平面分析-整体刚度矩阵

1、2022-7-4平面问题有限元分析-总刚15 平面问题有限元分析平面问题有限元分析整体刚度矩阵整体刚度矩阵曹国华曹国华5.1 整体刚度矩阵整体刚度矩阵5.2 整体刚度矩阵的特点整体刚度矩阵的特点5.3 边界条件边界条件5.4 计算结果整理计算结果整理5.5 收敛准则收敛准则 前文对单元体进行了分折，得到了单元刚度方程 ，但要解决问题，还必须进一步建立整个计算模型的整体刚度方程。完成这一步的关键，在于怎样将单元的刚度矩阵和节点荷载列阵，分别“组装”成整体刚度矩阵和整体节点荷载列阵。 这里通过研究任意节点的平衡来建立整体刚度矩阵，该方法不但比较直观、易懂，而且对怎样编写计算机程序是很有帮助的。ee

2、eFK 2022-7-42平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵 整体节点载荷列阵：由各节点所受等效节点力按节点号码以从小到大的顺序排列组成的列阵。等效节点力是由集中力、表面力和体积力共同移置构成的，其中集中力包括直接作用在弹性体上的外力和边界约束力，如支座反力。为了研究整体刚度矩阵的组装过程，先引入两个概念。 整体节点位移列阵：由各节点位移按节点号码以从小到大的顺序排列组成的列阵。2022-7-43平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵式中： ，TTTTTT12345 Tiiiuv (1,2,5)i 不失一般性，仅考虑计算模型中有4个单元，如图所示。四个单

3、元的整体节点位移列阵为2022-7-44平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵 对每个单元都可以写出相应的单元刚度方程 ，即单元节点平衡方程。例如，对号单元，有 eeeFK (1)(1)(1)(1)11112131(1)(1)(1)(1)22122232(1)(1)(1)(1)33132333FKKKFKKKFKKK 式中：(1)iF号单元中第i（i=1,2,3）节点所受力。为了便于组装整体刚度矩阵，将上式以整体节点位移 表示，即 (1)(1)(1)(1)11111213(1)(1)(1)(1)22212223(1)(1)(1)(1)(1)333132334500000000

4、0000000000FKKKFKKKKFKKK (1)K号单元的扩大刚度矩阵或称为单元贡献矩阵。2022-7-45平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵 同理，对于单元，有(2)(2)(2)(2)111113142(2)(2)(2)(2)(2)33313334(2)(2)(2)(2)444143445000000000000000000FKKKKFKKKFKKK (2)iF(2)K式中：号单元中第i（i=1,3,4）节点所受力； 号单元的扩大刚度矩阵。2022-7-46平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵对于3单元，有1(3)(3)(3)(3)222232

5、52(3)(3)(3)(3)(3)332333534(3)(3)(3)(3)55253555000000000000000000 FKKKFKKKKFKKK(3)iF(3)K式中：3号单元中第i（i=3,4,5）节点所受力； 3号单元的扩大刚度矩阵。2022-7-47平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵对于单元，有12(4)(4)(4)(4)(4)33334353(4)(4)(4)(4)44344454(4)(4)(4)(4)55354555000000000000000000FKKKKFKKKFKKK (4)iF(4)K式中：号单元中第i（i=3,4,5）节点所受力； 号

6、单元的扩大刚度矩阵。2022-7-48平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵 对于任意一个节点，可能承受两种力的作用，一种是其它单元给予该节点的反作用力；另一种是作用在节点上的等效节点力。对整体而言，前者属于内力，后者属于外力，每个节点在两种力的作用下处于平衡。 将各单元刚度方程左边相加，即将各节点所受力相加，由于对于整体而言，单元给予节点的反作用力属于内力，在相加过程中相互抵消，所以各节点所受力相加的结果只有外力，即等效节点力，从而得到整体节点荷载列阵，如下(1)(2)111(3)(1)222(3)(4)(1)(2)33333(4)(2)444(3)(4)555000000

7、00FFFFFFF =FFFFFFFFFFF2022-7-49平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵 将各单元刚度方程右边相加，从而得到整体刚度矩阵，如下(1)(2)(3)(4)11121314152122232425313233343541424344455152535455(1)(2)(1)(1)(2)(2)111112131314(1)(1)(3)(1)(3)(3)212222232325(1)3100KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK(2)(1)(3)(1)(2)(3)(4)(2)(4)(3)(4)31323233

8、33333334343535(2)(2)(4)(2)(4)(4)414343444445(3)(3)(4)(4)(3)(4)52535354555500KKKKKKKKKKKKKKKKKKKKKK2022-7-410平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵 通过以上分析得，整体节点载荷与整体节点位移之间的关系式，即结构整体有限元方程，如下FK 式中：K整体刚度矩阵。4( )1iiKK2022-7-411平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵 整体刚度矩阵组装的基本步骤：1）将单元刚度矩阵中的每个子块放到在整体刚度矩阵中的对应位置上，得到单元的扩大刚度矩阵。注意

9、对于单元刚度矩阵是按照局部编码排列的，即对应单元刚度矩阵中的i、j、m；对于整体刚度矩阵是按照整体编码排列的，即按节点号码以从小到大的顺序排列。在组装过程中，必须知道单元节点的局部编码与该节点在整体结构中的整体编码之间的关系，才能得到单元刚度矩阵中的每个子块在整体刚度矩阵中的位置。将单元刚度矩阵中的每个子块按总体编码顺序重新排列后，可以得到单元的扩大矩阵。例如在图中，单元的局部编码为i、j、m，对应整体编码为1、3、4，然后将单元刚度矩阵中的每个子块按总体编码顺序重新排列后，可以得到单元的扩大矩阵。注意有些书籍中将局部编码表示为1、2、3或1，2，3等；2）将全部单元的扩大矩阵相加得到整体刚度

10、矩阵。2022-7-412平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵通过以上组装过程可以得到组装整体刚度矩阵的一般规则： 1）结构中的等效节点力是相关单元结点力的叠加，整体刚度矩阵的子矩阵是相关单元的单元刚度矩阵子矩阵的集成； 2）当整体刚度矩阵中的子矩阵 中r=s时，该节点（节点r或s）被哪几个单元所共有，则 就是这几个单元的刚度矩阵中的子矩阵 的相加。如 应该是单元中对应子矩阵的集成，即 rsKrsKersK33K(1)(2)(3)(4)3333333333KKKKK2022-7-413平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵 3）当 中 时，若rs边是组合体

11、的内边，则 就是共用该边的两相邻单元刚度矩阵中的子矩阵 的相加。如13边为单元和的共用边，则 r srsKrsKersK(1)(2)131313K= KK 4）当 中r和s不同属于任何单元时，则 =0。如节点r1和 s5不同属于任何单元，此时 =0。rsKrsK15K2022-7-414平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵 上述组装基本步骤和规则具有普遍意义 对于不同类型、不同形式的单元，只是相应节点的子矩阵的阶数（节点自由度节点自由度）可能不同，至于组装整体刚度矩阵的规律仍是相同的。正是因为有了这种组装规律，使得有限元法能够很方便地应用电子计算机进行计算。2022-7-4

12、15平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵 例：例： 如图所示有限元模型，弹性模量为 ，厚度为 ，为简化计算取 ，求整体刚度矩阵。t0E2022-7-416平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵单元编号单元编号整体编码整体编码1、2、32、4、55、3、23、5、6局部编码局部编码i、j、mi、j、mi、j、mi、j、m以整体编码表以整体编码表示的单元刚度示的单元刚度矩阵子块矩阵子块解：解：该模型中共有6个节点，4个单元，各单元的信息如表所示。 各单元信息(1)(1)(1)111213(1)(1)(1)212223(1)(1)(1)313233KKKKKKK

13、KK(2)(2)(2)222425(2)(2)(2)424445(2)(2)(2)525455KKKKKKKKK(3)(3)(3)555352(3)(3)(3)353332(3)(3)(3)252322KKKKKKKKK(4)(4)(4)333536(4)(4)(4)535556(4)(4)(4)636566KKKKKKKKK2022-7-417平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵同上例类似的分析，得(2)(2)(2)222425(2)(2)(2)(2)424445(2)(2)(2)5254551011010202001031211213014002020101101EtK

14、KKKKKKKKK 根据单元刚度矩阵的性质，可知 ，若3单元5,3,2，则 整体刚度矩阵中的各子块是对所有单元相应的子块求和得到的（实际只是对相关单元求和），其中各子块矩阵均为2行2列，整体刚度矩阵用子块矩阵可以表示为(2)(1)(4)KK= K(2)(3)KK2022-7-418平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵111213141516212223242526313233343536414243444546515253545556616263646566KKKKKKKKKKKKKKKKKKK =KKKKKKKKKKKKKKKKKK上式中任意一子块矩阵均为2行2列，在计算

15、过程中，无需将每个单元刚度矩阵进行扩大，只需判断整体刚度矩阵子块的下标，然后利用组装整体刚度矩阵的一般规则进行计算，如 ，由图形可知，节点2由单元、和所共有，则 22K(1)(2)(3)22222222KKK+ K2022-7-419平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵25K(2)(3)252525KKK，由图形可知，25边为单元和的共用边，则15K150K，由图形可知，节点1、5不同属于任何单元，则采用同样的方法进行计算，得到整体刚度矩阵为2022-7-420平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵10110100000002020000000010614

16、11101001216120210000041610021011012160014000010003121004001200130100000121206121001014111601000000002020000010001101EtK2022-7-421平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵MATHCAD例子例子(1)(2)(1)(1)(2)(2)111112131314(1)(1)(3)(1)(3)(3)212222232325(1)(2)(1)(3)(1)(2)(3)(4)(2)(4)(3)(4)313132323333333334343535(2)(2)(4)(2)

17、(4)(4)414343444445520000KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK(3)(3)(4)(4)(3)(4)5353545555KKKKK公共邮箱：公共邮箱：密码：密码： cncumtcncumt 作业：作业： 采用采用MATHCAD求解求解如图所示有限元模型，弹性模量为 ，厚度为 ，为简化计算取 ，E=1，t=1 ，求整体刚度矩阵。t0E2022-7-424平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵每个班由负责人压缩成一每个班由负责人压缩成一个压缩包后发到个压缩包后发到每个人提交的文件名为：每个人提交的文件名为：班级班级-姓名姓名-序号（序

18、号（1）比如：比如：机自机自10（1）班）班-某某某某-12（1）4. 4. 是奇异矩阵，在排除刚体位移后，它是正定阵是奇异矩阵，在排除刚体位移后，它是正定阵 K3. 3. 是稀疏矩阵，非零元素呈带状分布是稀疏矩阵，非零元素呈带状分布K 用有限元方法分析复杂工程问题时，节点的数目比较多，整体刚度矩阵的阶数通常也是很高的。那么在进行计算时，如果存储整体刚度矩阵的全部元素，将会浪费较大的资源、降低计算效率。如果根据整体刚度矩阵的特点进行编写程序，可以大大节省资源、并提高计算效率。因此有必要了解和掌握整体刚度矩阵的特点，整体刚度矩阵具有以下几个显著的特点：1. 1. 是对称矩阵是对称矩阵K2. 2.

19、 中主对角元素总是正的中主对角元素总是正的K2022-7-425平面问题有限元分析-总刚5.2 整体刚度矩阵的特点整体刚度矩阵的特点1. 是对称矩阵K 由单元刚度矩阵的对称性和整体刚度矩阵的组装过程，可知整体刚度矩阵必为对称矩阵。利用对称性，在计算机编写程序时，只存储整体刚度矩阵上三角或下三角部分即可。2. 中主对角元素总是正的K 例如，刚度矩阵 中的元素 表示节点2在x方向产生单位位移，而其它位移均为零时，在节点2的x方向上必须施加的力； 表示节点2在y方向产生单位位移，而其它位移均为零时，在节点2的y方向上必须施加的力。很显然在此情况下力的方向应该与位移方向一致，故应为正号。K33K44K

20、2022-7-426平面问题有限元分析-总刚5.2 整体刚度矩阵的特点整体刚度矩阵的特点3. 是稀疏矩阵，非零元素呈带状分布K 如果遵守一定的节点编号规则，就可使矩阵的非零元素都集中在主对角线附近呈带状。整体刚度矩阵中的子矩阵 只有当下标s等于r或者s与r同属于一个单元时才不为零，这就说明，在第r双行中非零子矩阵的块数，应该等于节点r周围直接相邻的节点数目加1。可见， 中元素一般都不是填满的，是稀疏矩阵，且非零元素呈带状分布。rsKK 以下图所示的单元网格为例，其整体刚度矩阵中的非零子块（每个子块为2行2列）的分布情况如下图所示，图中阴影部分表示该子块不为零，其它子块部位均为零。2022-7-

21、427平面问题有限元分析-总刚5.2 整体刚度矩阵的特点整体刚度矩阵的特点2022-7-428平面问题有限元分析-总刚5.2 整体刚度矩阵的特点整体刚度矩阵的特点1011010000000202000000001061411101001216120210000041610021011012160014000010003121004001200130100000121206121001014111601000000002020000010001101EtK2022-7-429平面问题有限元分析-总刚5.1 整体刚度矩阵整体刚度矩阵 显然，带状刚度矩阵的带宽取决于单元网格中相邻节点号码的最大差值D

22、。把半个斜带形区域中各行所具有的非零元素的最大个数叫做刚度矩阵的半带宽（包括主对角元），用B表示，如下 B=2(D+1) 通常的有限元程序，一般都利用刚度矩阵的对称和稀疏带状的特点，在计算求解中，只存储上半带的元素，即所谓的半带存储。因此，在划分完有限元网格进行节点编号时，要采用合理的编码方式，使同一单元中相邻两节点的号码差尽可能小，以便节省存储空间、提高计算效率。2022-7-430平面问题有限元分析-总刚5.2 整体刚度矩阵的特点整体刚度矩阵的特点 对于同样的有限元单元网格，按照图（a）的结点编码，最大的半带宽为 按照图（b）的结点编码，最大的半带宽为B 按照图（c）的结点编码，最大的半带

23、宽为B (a) (b) (c)2022-7-431平面问题有限元分析-总刚2(106)105.2 整体刚度矩阵的特点整体刚度矩阵的特点B = 2(104)=142(102)184. 是奇异矩阵，在排除刚体位移后，它是正定阵 K 无约束的弹性体（或结构物）的整体刚度矩阵是奇异的，不存在逆矩阵，即关于位移的解不唯一。这是因为弹性体在外力的作用下处于平衡，外力的分量应该满足三个静力平衡方程。这反映在整体刚度矩阵中就意味着存在三个线性相关的行或列，所以是个奇异阵，不存在逆矩阵。例如：设弹性体在外力的作用下处于平衡，这时相应的解为 ，然后在给予弹性体以刚体位移而相应的节点位移 ，这时， 仍是问题的解，因

24、为刚体位移不会破坏平衡。1 2 12 注：当排除刚体位移后，整体刚度矩阵是正定矩阵。2022-7-432平面问题有限元分析-总刚5.2 整体刚度矩阵的特点整体刚度矩阵的特点 前文已经提到在排除刚体位移后，整体刚度矩阵是正定的，方程才可求得唯一解。排除刚体位移可以通过引入边界约束条件来实现，这里介绍两种比较简单的引入已知位移的方法。1代入法2乘大数法2022-7-433平面问题有限元分析-总刚5.3 边界条件边界条件1代入法 该方法保持方程组仍为2n2n系统，仅对整体刚度矩阵 和整体载荷列阵进行修正。下面以一个只有四个方程的简单例子加以说明，方程如下11121314112122232412313

25、23334234142434424KKKKuFKKKKvFKKKKuFKKKKvF 假定系统中节点位移 、 ，则当引入这些节点的已知位移之后，方程就变成11uc22uc2022-7-434平面问题有限元分析-总刚5.3 边界条件边界条件1112131411212223241231323334234142434424KKKKuFKKKKvFKKKKuFKKKKvF 11221 123 22224122441 143 242442100000001000cuFK cK cKKvcuFK cK cKKv若 ，则120cc1222412242442410000000010000uKKvFuKKvF 然

26、后，用这组维数不变的方程来求解所有的节点位移。显然，其解答仍为原方程的解答。 在手算时，可直接将零位移约束所对应的整体刚度矩阵中的行和列直接划去，使得整体刚阵的维数变小，更便于手算。2022-7-435平面问题有限元分析-总刚5.3 边界条件边界条件2乘大数法 将 中与指定的节点位移相对应的主对角元素乘上一个大数，同时将 中的对应元素换成指定的节点位移值、该大数与节点位移相对应的主对角元素三者的乘积。若把此方法用于上面的例子，则方程就变成KF151511112131411112122232421515231323334233241424344410101010uKKKKc KvKKKKFuKK

27、KKc KvKKKKF该方程组的第一个方程为151511112 11321421111010KuK vK uK vc K解得 ，这种方法就是使 中相应行的修正项远大于非修正项。11uc2022-7-436平面问题有限元分析-总刚5.3 边界条件边界条件1代入法代入法2乘大数法乘大数法 在以上的两种方法中，代入法接近人工解法，虽然该方法比较直观，但该方法对刚度矩阵改变较多，程序效率不高。乘大数法对刚度矩阵改变较少，工作量较小，但相乘的“大数大数”若取得过大，求解时会发生“溢出溢出”、若取得太小则会引起较大的误差误差。2022-7-437平面问题有限元分析-总刚5.3 边界条件边界条件 对于三节点

28、三角形单元，单元内各点的应力值相等，算出的应力一般作为单元形心处的应力。由于单元应力为常数，整个结构的应力场呈阶梯状，在单元之间不连续。而工程上往往更加关心边界和节点上的受力情况，因此，必须对所得到的应力再次进行处理，得到更加合理的应力场，并得到所需点上的应力值。 这里介绍两种简单的方法，一种方法称为节点平均法，即把环绕某一节点的各单元的应力加以平均作为该节点的应力值。例如图中节点3的应力为(1)(2)(3)(4)(5)(6)316 2022-7-438平面问题有限元分析-总刚5.4 计算结果整理计算结果整理 为了使通过这样平均得来的应力比较接近实际情况，要求环绕节点的各单元尺寸不应相差太大。

29、这种做法，对内节点比较好，对边界点则可能很差。 因此，边界节点处的应力不宜直接由单元应力平均来获得，而要根据内节点的应力构造插值函数推算出来。例如图中边界点1的应力，可以先用节点平均法求得节点2、3、4处的应力，在构造相应的插值函数推算边界点1的应力，如常用的抛物线插值公式如下3424231234232432344243xxxxxxxxxxxxxxxxxxxxxxxx 5.4 计算结果整理计算结果整理2022-7-439平面问题有限元分析-总刚 另一种方法称为单元平均法，即把两相邻单元的应力加以平均，用以表示公共边界中点的应力。为了由这样平均所得到的应力具有较好的精度，两相邻单元的面积不应相差

30、太大。如图中单元和边界的中点处的应力为(2,3)(1)(2)12 在不同的有限元软件中均具有各自的后处理方法，但无论后期怎样处理，应力场来源于应力的计算，应力的精度主要依赖于单元的尺寸、单元的类型（位移模式）。2022-7-440平面问题有限元分析-总刚5.4 计算结果整理计算结果整理 对于有限元这种数值计算方法，一般总是希望随着网格的逐步细分所得到的解能够收敛于问题的精确解。根据前面的分析，可知在有限元分析中，一旦确定了单元的形状之后，位移模式的选择将是非常关键的。由于载荷的移置、应力矩阵和刚度矩阵的建立等等，都依赖于单元的位移模式，所以，如果所选择的位移模式与真实的位移分布有很大的差别，那

31、么就很难获得良好的数值解。为了保证解答的收敛性，要求位移模式必须满足以下三个条件，即 1）位移模式必须包含单元的刚体位移。 2）位移模式必须包含单元的常应变。 3）位移模式在单元内要连续、且在相邻单元之间的位移必须协调。 5.5 收敛准则收敛准则2022-7-441平面问题有限元分析-总刚 1）位移模式必须包含单元的刚体位移。也就是说，当节点位移是由某个刚体位移所引起时，弹性体内将不会产生应变。所以，位移模式不但要具有描述单元本身形变的能力，而且还要具有描述由于其它单元形变而通过节点位移引起单元刚体位移的能力。 2）位移模式必须包含单元的常应变。每个单元的应变一般都是包含着两个部分：一部分是与该单元中各点的坐