1,5 平面问题有限元分析整体刚度矩阵,曹国华,整体刚度矩阵 整体刚度矩阵的特点 边界条件处理 计算结果整理,2,前文对单元体进行了分折，得到了单元刚度方程 ，但要解决问题，还必须进一步建立整个计算模型的整体刚度方程。完成这一步的关键，在于怎样将单元的刚度矩阵和节点荷载列阵，分别“组装”成整体刚度矩阵和整体节点荷载列阵。 通过研究任意节点的平衡来建立整体刚度矩阵，该方法不但比较直观、易懂，而且对怎样编写计算机程序是很有帮助。,整体刚度矩阵,3,整体节点载荷列阵：由各节点所受等效节点力按节点号码以从小到大的顺序排列组成的列阵。等效节点力是由集中力、表面力和体积力共同移置构成的，其中集中力包括直接作用在弹性体上的外力和边界约束力，如支座反力。,为了研究整体刚度矩阵的组装过程，先引入两个概念。,整体节点位移列阵：由各节点位移按节点号码以从小到大的顺序排列组成的列阵。,整体刚度矩阵,4,式中： ，,不失一般性，仅考虑计算模型中有4个单元，如图所示。四个单元的整体节点位移列阵为,整体刚度矩阵,5,对每个单元都可以写出相应的单元刚度方程 ，即单元节点平衡方程。例如，对号单元，有,为了便于组装整体刚度矩阵，将上式以整体节点位移 表示，即,号单元的扩大刚度矩阵或称为单元贡献矩阵。,整体刚度矩阵,6,同理，对于单元，有,整体刚度矩阵,7,对于3单元，有,整体刚度矩阵,8,对于单元，有,整体刚度矩阵,9,对于任意一个节点，可能承受两种力的作用，一种是其它单元给予该节点的反作用力；另一种是作用在节点上的等效节点力。对整体而言，前者属于内力，后者属于外力，每个节点在两种力的作用下处于平衡。 将各单元刚度方程左边相加，即将各节点所受力相加，由于对于整体而言，单元给予节点的反作用力属于内力，在相加过程中相互抵消，所以各节点所受力相加的结果只有外力，即等效节点力，从而得到整体节点荷载列阵，如下,整体刚度矩阵,10,将各单元刚度方程右边相加，从而得到整体刚度矩阵，如下,整体刚度矩阵,11,通过以上分析得，整体节点载荷与整体节点位移之间的关系式，即结构整体有限元方程，如下,整体刚度矩阵,12,整体刚度矩阵组装的基本步骤： 1）将单元刚度矩阵中的每个子块放到在整体刚度矩阵中的对应位置上，得到单元的扩大刚度矩阵。注意对于单元刚度矩阵是按照局部编码排列的，即对应单元刚度矩阵中的i、j、m；对于整体刚度矩阵是按照整体编码排列的，即按节点号码以从小到大的顺序排列。在组装过程中，必须知道单元节点的局部编码与该节点在整体结构中的整体编码之间的关系，才能得到单元刚度矩阵中的每个子块在整体刚度矩阵中的位置。将单元刚度矩阵中的每个子块按总体编码顺序重新排列后，可以得到单元的扩大矩阵。例如在图中，单元的局部编码为i、j、m，对应整体编码为1、3、4，然后将单元刚度矩阵中的每个子块按总体编码顺序重新排列后，可以得到单元的扩大矩阵。注意有些书籍中将局部编码表示为1、2、3或1，2，3等； 2）将全部单元的扩大矩阵相加得到整体刚度矩阵。,整体刚度矩阵,13,通过以上组装过程可以得到组装整体刚度矩阵的一般规则： 1）结构中的等效节点力是相关单元结点力的叠加，整体刚度矩阵的子矩阵是相关单元的单元刚度矩阵子矩阵的集成；,2）当整体刚度矩阵中的子矩阵 中r=s时，该节点（节点r或s）被哪几个单元所共有，则 就是这几个单元的刚度矩阵中的子矩阵 的相加。如 应该是单元中对应子矩阵的集成，即,整体刚度矩阵,14,3）当 中 时，若rs边是组合体的内边，则 就是共用该边的两相邻单元刚度矩阵中的子矩阵 的相加。如13边为单元和的共用边，则,4）当 中r和s不同属于任何单元时，则 =0。如节点r1和 s5不同属于任何单元，此时 =0。,整体刚度矩阵,15,16,上述组装基本步骤和规则具有普遍意义 对于不同类型、不同形式的单元，只是相应节点的子矩阵的阶数（节点自由度节点自由度）可能不同，至于组装整体刚度矩阵的规律仍是相同的。正是因为有了这种组装规律，使得有限元法能够很方便地应用电子计算机进行计算。,17,例： 如图所示有限元模型，弹性模量为 ，厚度为 ，为简化计算取 ，求整体刚度矩阵。,18,解：该模型中共有6个节点，4个单元，各单元的信息如表所示。,各单元信息,19,同上例类似的分析，得,根据单元刚度矩阵的性质，可知 ，若3单元5,3,2，则 整体刚度矩阵中的各子块是对所有单元相应的子块求和得到的（实际只是对相关单元求和），其中各子块矩阵均为2行2列，整体刚度矩阵用子块矩阵可以表示为,20,上式中任意一子块矩阵均为2行2列，在计算过程中，无需将每个单元刚度矩阵进行扩大，只需判断整体刚度矩阵子块的下标，然后利用组装整体刚度矩阵的一般规则进行计算，如 ，由图形可知，节点2由单元、和所共有，则,21,，由图形可知，25边为单元和的共用边，则,，由图形可知，节点1、5不同属于任何单元，则,采用同样的方法进行计算，得到整体刚度矩阵为,22,23,整体刚度矩阵简单算例2,已知平面矩形结构，边长为1，E=1，t=1，NU=0.25,约束及外载形式如图所示。假定划分为两个单元。,24,单元1节点总码：1，2，4；单元2节点总码：3，4，2；逆时针排列。单元1刚度矩阵为：,K22,25,单元2刚度矩阵与单元1相同，注意单刚中的k33块在整体中为K22。,K22,单刚来说为K33,26,单元2节点：2，3，4顺序，单元刚度矩阵为,27,单元1刚度矩阵扩充到整体刚度矩阵为,28,单元2刚度矩阵扩充到整体刚度矩阵为,K22,29,30,作业： 如图所示有限元模型，弹性模量为 ，厚度 ，为简化计算取 ，E=1，t=1 ，求整体刚度矩阵。,31,4. 是奇异矩阵，在排除刚体位移后，它是正定阵,3. 是稀疏矩阵，非零元素呈带状分布,用有限元方法分析复杂工程问题时，节点的数目比较多，整体刚度矩阵的阶数通常也是很高的。那么在进行计算时，如果存储整体刚度矩阵的全部元素，将会浪费较大的资源、降低计算效率。如果根据整体刚度矩阵的特点进行编写程序，可以大大节省资源、并提高计算效率。因此有必要了解和掌握整体刚度矩阵的特点，整体刚度矩阵具有以下几个显著的特点：,1. 是对称矩阵,2. 中主对角元素总是正的,整体刚度矩阵的特点,32,1. 是对称矩阵,由单元刚度矩阵的对称性和整体刚度矩阵的组装过程，可知整体刚度矩阵必为对称矩阵。利用对称性，在计算机编写程序时，只存储整体刚度矩阵上三角或下三角部分即可。,2. 中主对角元素总是正的,例如，刚度矩阵 中的元素 表示节点2在x方向产生单位位移，而其它位移均为零时，在节点2的x方向上必须施加的力； 表示节点2在y方向产生单位位移，而其它位移均为零时，在节点2的y方向上必须施加的力。很显然在此情况下力的方向应该与位移方向一致，故应为正号。,33,3. 是稀疏矩阵，非零元素呈带状分布,如果遵守一定的节点编号规则，就可使矩阵的非零元素都集中在主对角线附近呈带状。整体刚度矩阵中的子矩阵 只有当下标s等于r或者s与r同属于一个单元时才不为零，这就说明，在第r双行中非零子矩阵的块数，应该等于节点r周围直接相邻的节点数目加1。可见， 中元素一般都不是填满的，是稀疏矩阵，且非零元素呈带状分布。,以下图所示的单元网格为例，其整体刚度矩阵中的非零子块（每个子块为2行2列）的分布情况如下图所示，图中阴影部分表示该子块不为零，其它子块部位均为零。,34,35,36,显然，带状刚度矩阵的带宽取决于单元网格中相邻节点号码的最大差值D。把半个斜带形区域中各行所具有的非零元素的最大个数叫做刚度矩阵的半带宽（包括主对角元），用B表示，如下 B=2(D+1) 通常的有限元程序，一般都利用刚度矩阵的对称和稀疏带状的特点，在计算求解中，只存储上半带的元素，即所谓的半带存储。因此，在划分完有限元网格进行节点编号时，要采用合理的编码方式，使同一单元中相邻两节点的号码差尽可能小，以便节省存储空间、提高计算效率。,37,对于同样的有限元单元网格，按照图（a）的结点编码，最大的半带宽为 按照图（b）的结点编码，最大的半带宽为B 按照图（c）的结点编码，最大的半带宽为B,(a) (b) (c),2(106+1)10,B = 2(104+1)=14,2(102+1)18,38,4. 是奇异矩阵，在排除刚体位移后，它是正定阵,无约束的弹性体（或结构物）的整体刚度矩阵是奇异的，不存在逆矩阵，即关于位移的解不唯一。这是因为弹性体在外力的作用下处于平衡，外力的分量应该满足三个静力平衡方程。这反映在整体刚度矩阵中就意味着存在三个线性相关的行或列，所以是个奇异阵，不存在逆矩阵。例如：设弹性体在外力的作用下处于平衡，这时相应的解为 ，然后在给予弹性体以刚体位移而相应的节点位移 ，这时， 仍是问题的解，因为刚体位移不会破坏平衡。,注：当排除刚体位移后，整体刚度矩阵是正定矩阵。,39,前文已经提到在排除刚体位移后，整体刚度矩阵是正定的，方程才可求得唯一解。排除刚体位移可以通过引入边界约束条件来实现，这里介绍两种比较简单的引入已知位移的方法。,1代入法,2乘大数法,边界条件,40,1代入法,该方法保持方程组仍为2n2n系统，仅对整体刚度矩阵 和整体载荷列阵进行修正。下面以一个只有四个方程的简单例子加以说明，方程如下,假定系统中节点位移 、 ，则当引入这些节点的已知位移之后，方程就变成,41,若 ，则,然后，用这组维数不变的方程来求解所有的节点位移。显然，其解答仍为原方程的解答。 在手算时，可直接将零位移约束所对应的整体刚度矩阵中的行和列直接划去，使得整体刚阵的维数变小，更便于手算。,42,2乘大数法,将 中与指定的节点位移相对应的主对角元素乘上一个大数，同时将 中的对应元素换成指定的节点位移值、该大数与节点位移相对应的主对角元素三者的乘积。若把此方法用于上面的例子，则方程就变成,该方程组的第一个方程为,解得 ，这种方法就是使 中相应行的修正项远大于非修正项。,43,1代入法,2乘大数法,在以上的两种方法中，代入法接近人工解法，虽然该方法比较直观，但该方法对刚度矩阵改变较多，程序效率不高。乘大数法对刚度矩阵改变较少，工作量较小，但相乘的“大数”若取得过大，求解时会发生“溢出”、若取得太小则会引起较大的误差。,44,对于三节点三角形单元，单元内各点的应力值相等，算出的应力一般作为单元形心处的应力。由于单元应力为常数，整个结构的应力场呈阶梯状，在单元之间不连续。而工程上往往更加关心边界和节点上的受力情况，因此，必须对所得到的应力再次进行处理，得到更加合理的应力场，并得到所需点上的应力值。 这里介绍两种简单的方法，一种方法称为节点平均法，即把环绕某一节点的各单元的应力加以平均作为该节点的应力值。例如图中节点3的应力为,计算结果整理,45,为了使通过这样平均得来的应力比较接近实际情况，要求环绕节点的各单元尺寸不应相差太大。这种做法，对内节点比较好，对边界点则可能很差。 因此，边界节点处的应力不宜直接由单元应力平均来获得，而要根据内节点的应力构造插值函数推算出来。例如图中边界点1的应力，可以先用节点平均法求得节点2、3、4处的应力，在构造相应的插值函数推算边界点1的应力，如常用的抛物线插值公式如下,计算结果整理,46,另一种方法称为单元平均法，即把两相邻单元的应力加以平均，用以表示公共边界中点的应力。为了由这样平均所得到的应力具有较好的精度，两相邻单元的面积不应相差太大。如图中单元和边界的中点处的应力为,在不同的有限元软件中均具有各自的后处理方法，但无论后期怎样处理，应力场来源于应力的计算，应力的精度主要依赖于单元的尺寸、单元