第四章 结构刚度矩阵的存贮方式和组集程序

1. 第四章 结构刚度矩阵的存贮方式和组集程序

2. 前面几章，我们介绍了有限元法的基本原理，较深入地讨论了单元分析和系统分析理论。有限元法的重要价值在于方便编写程序，借助于计算机完成全部计算工作。 在整个编程和计算中，影响最大的是结构刚度矩阵的存贮和组集。本章就对这一问题进行专门讨论。 结构刚度矩阵的组集方法、程序和它的存贮方式密切相关，还依赖于约束条件的处理方式。本章介绍三种存贮方式：方阵存贮、等带宽存贮和一维数组存贮。

3. 约束条件采用第三章思想处理：从无约束的结构刚度矩阵中删去与受约束位移号对应的行和列，再将矩阵压缩排列成n×n阶方阵。实际操作时，根本就不去存贮这些行和列。约束条件采用第三章思想处理：从无约束的结构刚度矩阵中删去与受约束位移号对应的行和列，再将矩阵压缩排列成n×n阶方阵。实际操作时，根本就不去存贮这些行和列。 本章讨论的结构刚度矩阵存贮方式和组集的基本内容，适合于各类结构，包括杆系结构、弹性连续结构等。适合于结构静力分析，也是结构动力分析以及稳定性分析的重要基础。 下面，我们就从方阵存贮方式开始。

4. 5.1 结构刚度矩阵的方阵存贮 （3-8） 为叙述简便起见，本章中所用“单元刚度矩阵”一词都指结构坐标单元刚度矩阵。只有特别需要时才加“结构坐标”定语。 由第三章式（3-8）知 结构刚度矩阵[K]等于所有膨胀后的单元刚度矩阵[k]相累加。下面详细考察如何在程序中实现这个操作。 1、单元刚度矩阵元素的两类地址 单刚地址和总刚地址。

5. （1）单刚地址 单刚地址是指单元刚度矩阵元素在膨胀前的单元刚度矩阵中的地址，即行、列号。式（4-1）矩阵中每个元素的两个下标就是该元素的单刚地址。第一个下标指行号，第二个下标指列号。例如k15的说明它是矩阵中第1行第5列位置的元素。 （4-1）

6. （2）总刚地址 总刚地址是指单元刚度矩阵元素在膨胀后的单元刚度矩阵中的地址（行、列号）。这个地址与单元的节点号有关，而节点号是离散结构时编好的，还与节点自由度有关。 Y i j X i3 j3 i1 j2 l i2 j1 图4-1 设式（4-1）是图4-1所示梁单元的单元刚度矩阵。该单元两个节点的节点号为i、j，节点自由度是3。 单元刚度矩阵元素的总刚地址应为：i1、i2、i3、j1、j2、j3，如式（4-2）。

7. i1 i2 j3 i3 j1 j2 i1 i2 i3 （4-2） j1 j2 j3 （4-3） 有

8. m2 Y m1 j2 i2 j1 i1 m j i X （4-4） 如果式（4-1）是图4-2所示平面问题三角形单元的单元刚度矩阵。该单元三个节点的节点号为i、j、m，节点自由度是2。 图4-2 单元刚度矩阵元素的总刚地址应为：i1、i2、j1、j2、m1、m2，如式（4-4）、（4-5）。 有

9. i1 i2 m2 j1 j2 m1 i1 i2 j1 j2 m1 m2 （4-5）

10. （3）两类地址的连接数组 （4-6） （4-7） [LM]=[ 3i-2 3i-1 3i 3j-2 3j-1 3j ] [LM]=[ 2i-1 2i 2j-1 2j 2m-1 2m ] 建立一个连接数组供计算机识别两类地址，以便正确地从单元刚度矩阵中取出元素，并放入结构刚度矩阵之中。 给连接数组取个名，例如：LM（NDF2），它是一维数组，数组长度NDF2是单元自由度的最大值。 先不考虑约束情况： 对于平面梁单元组成的结构 对于平面问题三角形单元组成的结构

11. （4）建立连接数组的条件 从式（4-6）、（4-7）知，建立连接数组应先知道单元的节点号i、j、m 。这些数据应作为结构计算的初始数据在程序开始后即读入的，因此是已知的。通常，这些数据读入后，按照单元号顺序由小到大依次存放在若干个一维数组中，例如，数组IO、JO、MO分别存放i、j、m节点号；数组IO、JO、MO的长度为单元总数。当然也可用一个多维数组存放它们。

12. 2、组集无约束结构刚度矩阵 ^ NE [K]=  [k]e e=1 （1）方法 现在我们指出一个明显而重要的事实。所谓单元刚度矩阵的膨胀只不过为了说明问题的概念。事实上，执行时并不去膨胀单元刚度矩阵。所谓把所有单元膨胀后的单元刚度矩阵累加起来，即算式 的实现操作方法非常简单： 只要把每个单元的单元刚度矩阵[k]的每个元素按照它的总刚地址送入总刚度矩阵[K]，并简单地累加起来就组集起了总刚度矩阵[K]。

13. （2）步骤 组集无约束结构总刚度矩阵的具体步骤如下： ① 矩阵[K]充零 即把[K]的N1×N1个元素全部置为零元素。 N1=NJ（节点总数）×NDF（一个节点自由度数） ② 从第1号单元开始，依次对所有单元执行下述操作：(下面符号“e”代表被执行的单元号) (a) 计算e号单元单元坐标系下的单元刚度矩阵[k]e。 (b) 计算e号单元的坐标变换矩阵[R]。 (c) 形成e号单元结构坐标系下的单元刚度矩阵[k]。计算公式为: [k]=[R]T[k]e[R]

14. (d) 建立e号单元的连接数组[LM] (e) 把e号单元结构坐标单元刚度矩阵[k]的元素按连接数组给出的总刚地址送入结构刚度矩阵[K]并进行累加。赋值语句为： K(i,j)=K(i,j)+k(p,q) （4-8） K(i,j) ——结构刚度矩阵中第i行j列元素 k(p,q) ——单元刚度矩阵中第p行q列元素 p, q ——单刚地址 i, j ——总刚地址 为了节省存贮、节省运算时间，执行第(e)步时应利用单元刚度矩阵和结构刚度矩阵的对称性。用式（4-8）只须将单元刚度矩阵的上三角部分（含对角线元素）装入结构刚度矩阵中，形成结构刚度矩阵的上三角部分。

15. 执行命令 K( j, i )=K( i, j ) ( i＜j ) 即得结构刚度矩阵[K]的下三角部分元素。 单元刚度矩阵中的上三角元素 组成 结构刚度矩阵中的上三角元素

17. 执行上述第②步时，对每个单元皆重复同样规律的计算，因而应编写独立的子程序。 [K] ← 0 M←（1，1，NE） ELEM（M） II ← IO[M]， JJ ← JO[M] I1←（1，1, NDF） （3）程序框图

18. LM(I1) ← NDF*(II-1)+I1 LM(I1+NDF) ← NDF*(JJ-1)+I1 R ← LM(I) I←(1，1，NDF2) J←(I+1, 1, NDF2) K(R,R)←K(R,R)+k(I,I) C ← LM [ J ] R < C K [R,C ] ← K [ R,C ] + k [I,J ] K [ C,R ] ← K [C,R ] + k [ J,I ]

19. 3、组集约束结构总刚度矩阵 所谓约束是指刚性约束。弹性约束可按弹性单元处理，而在刚性支承处引入约束。 组集起无约束结构的总刚度矩阵[K]之后，形成受约束结构的结构刚度矩阵Kff的原则是很简单的。例如，如果第L号位移分量受有刚性约束，即L=0。那么，只要把无约束结构刚度矩阵[K]的第L行和第L列删去即可。对所有受刚性约束的位移分量都作这样处理后，就得到约束结构的结构刚度矩阵[Kff]。

20. 下面我们就按照这个原则来组集[Kff]。但不是在形成[K]之后再从其中删去一些行和列，而是根本不去存贮应该删去的那些行和列。 （1）约束信息数组 假定所讨论的结构受刚性支承约束的节点共有NRJ个。我们用一维数组 KRJ(NRJ) 来存贮这些受刚性支承约束的节点的号码。 每个受约束的节点，它的各个自由度可能全受到约束，也可能只有一个自由度被约束。为了表明这种情况，还需一个二维数组: KRL(NDF,NRJ)

21. 数组[KRL]的第i列，说明数组[KRJ]中的第i个节点的每个自由度是否受约束。我们规定：元素为1受约束，0表示不受约束。数组[KRL]的第i列，说明数组[KRJ]中的第i个节点的每个自由度是否受约束。我们规定：元素为1受约束，0表示不受约束。 4 5 6 1 2 3 2 Y 1 3 图4-3 ○ X ○ ○ ○ ○ 例 图4-3所示框架结构NRJ=3 1号节点1、 2受约束， 3不 受约束； 2号节点3个位移均受约束；3号节点只2受约束。 数组KRJ、KRL给出了结构的全部约束信息。这些信息应作为初始数据交由计算机读入的。

22. （2）节点位移信息数组 2 Y 1 3 ○ X ○ ○ ○ ○ 节点位移信息数组可命名为ID(NDF,NJ)。NJ是节点总数。数组ID是一个NDF行, NJ列的二维数组。 ID数组先用来表明每个节点（无论它是否受约束）的各个自由是否受到约束。对图4-3情形，NDF=3，NJ=6，ID数组的内容是： 这一步可借助于KRL数组实现。 4 5 6 （4-9） 1 2 3

23. 再对ID数组进行如下改造：从第1列开始，自上而下检查，遇1改0，遇0则用上一个非0元素加1；然后，对第2列、第3列、…依次按同样办法检查处理，完成后，得最终形式的ID数组。仍以图4-3为例，对式（5-9）改造后，得: 4 5 6 1 2 3 2 Y 1 3 图4-3 ○ X ○ ○ ○ ○ 0 0 2 4 7 10 11 0 5 8 0 0 3 12 0 1 6 9 （4-10）

24. 改造后的ID数组的性质是: a) 如果 ID（i, j ) = 0 则表明j号节点第i个自由度受有约束。 b) 如果 ID( i, j ) ≠ 0 则j号节点第i个自由度不受约束。并且， j号节点第i个位移分量在非约束节点位移列向量 f的序号就是: ID( i, j )

25. ID数组形成框图 [ID] ← 0 J ←KRJ(M) M←（1，1，NRJ） I←（1，1, NDF） ID[ I,J ] ← KRL[ I,M ] KK←0 受约束节点总数 受约束节点号 受约束节点信息

26. J ←（1，1, NJ） I ←（1，1, NDF） K K ← K K + 1 ID(I,J) ← K K ID(I,J) ← 0 I D [ I，J ]=1

27. （3） 受约束结构总刚度矩阵的组集 组集受约束结构总刚度矩阵[Kff]和组集无约束结构总刚度矩阵[K]的区别在于前者引入了约束，而后者没有。因此，组集受约束结构总刚度矩阵[Kff]的步骤与组集无约束结构总刚度矩阵[K]的步骤完全相同。只须在组集无约束结构总刚度矩阵[K]的方法中，做以下三点修改，以反应约束的引入。 ① 结构刚度矩阵的体积从n1×n1改为n×n n1——不考虑约束的节点总自由度 n1=NDF×NJ

28. n——删去约束后，结构的节点自由度总数 n=n1-nr nr——受约束节点的自由度数 ② 根据ID数组形成连接数组 ③ 在单元刚度[k]的总刚地址中只要为0，则根本不把[k]中的相应行和列送入结构刚度矩阵。 执行上述过程，得到的就是受约束结构刚度矩阵[Kff]。

29. [K] ← 0 ELEM（M） M←（1，1，NE） II ← IO[M]， JJ ← JO[M] I1←（1，1, NDF） LM(I1) ← ID(I1,II) LM(I1+NDF) ← ID(I1,JJ) I←(1，1，NDF2) （4）程序框图

30. IR ← LM(I) K(IR,IR)←K(IR,IR)+k(I,I) J←(I+1, 1, NDF2) JC ← LM [ J ] K[IR,JC] ← K[IR,JC] + k [I,J ] IR < JC IR = 0 JC = 0 K [JC,IR] ← K[JC,IR] + k [ J,I ]

31. 4.2 结构刚度矩阵的等带宽存贮 1、结构刚度矩阵的稀疏性 在第3章里，我们曾经指出，结构刚度矩阵是具有大量零元素的稀疏矩阵。 现在对这一性质作进一步的考察。

32. 在对结构系统进行平衡分析时知，如果i节点发生单位位移，只有与节点i有直接联系的那些节点才会受到影响，产生节点力。反之， i节点是否产生节点力也只会受与它直接相连的节点的影响。 设结构中共有100个节点，其中第i（i = 4）号节点及其与i号节点直接联系的节点如图4-4所示。 j=7 5 i=4 图4-4 平面框架中的节点 2 则结构刚度矩阵中的3i-2, 3i-1, 3i行的形式为：

33. 10 21 1 2 12 ⋮ 10 11 12 ⋮ ⋮ 100 4 图4-5 上图中每个小方块是一个三阶子块阵，所示出的与节点i相应的100个三阶子块中，只有4个非零块，其余96个皆为零元素组成的三阶子块。通常，对于任何一个节点，与它直接联系的节点个数和问题中的节点总数比较起来总是小得多。因而，呈现出结构刚度矩阵的高度稀疏性。

34. UBW 0 UBW 对 称 0 图4-6 由于刚度矩阵为对称矩阵，只要对节点适当编号，就能使刚度矩阵的非零元素聚集在主对角线两侧，呈现为图4-6所示带状矩阵。

35. 一般，刚度矩阵各行的非零元素个数和排列是不同的。对第i行而言，从主对角线元素Kii算起，直到它最右边的一个非零元素为止，其间所有元素（包括可能有的零元素）的个数，称为该行的右带宽。以图4-5为例，第10行的右带宽为12，即包括 K10，10，K10，11，…，K10，21。 这12个元素。 把相互连接的两个节点i和j的号码差∣j – i∣称为相邻点号差，当i<j时从主对角线子块[Kii]到最右边的非零块[Kij] ，共有∣j – i∣+ 1个三阶子块；当i>j时，从[Kjj]到[Kji]也共有∣j – i∣+ 1个三阶子块。因而：

36. 10 21 1 2 12 ⋮ 10 11 12 ⋮ ⋮ 100 图中所示三行中最大的右带宽为： （∣j – i∣+ 1 ）* 3 4 图4-7

37. 刚度矩阵各行右带宽中之最大者，称为最大右带宽。因刚度矩阵为对称矩阵，最大右带宽和最大左带宽是相同的，故简称最大右带宽为最大半带宽，记为UBW（Upper Band Width ）。 当结构的节点编号确定后，刚度矩阵的最大半带宽随之确定。找出所有相邻点号差中的最大者，记为rmax，则刚度矩阵最大半带宽为 UBW =（1 + rmax ）* NDF （4-9） 其中，NDF为节点自由度。对平面刚架，NDF = 3。对弹性平面问题，NDF=2。

38. 2、结构刚度矩阵的等带宽存贮 鉴于工程结构都是几何不变结构，从现在起，我们只讨论受约束结构的刚度矩阵。为书写简便起见，用[K]代表受约束结构的刚度矩阵，即以前的[Kff]。 因为刚度矩阵是对称的，只须存贮它的一半即可。约定：只存贮[K]的上三角部分。带区之外全是零元素，自然不必存贮，因此只需存贮包括主对角元素在内的上半带宽区域内的元素。这样，大大节约了需占用的计算机内存。称这种存贮方式为等带宽存贮。 设N为独立结点位移分量总数，则原始形式的刚度矩阵[K]是一个N行N列的方阵。所谓等带宽存贮，就是把[K]的上半带宽区域内全部元素用一个N行UBW列的矩形数组予以存贮：

39. UBW UBW 行 号 行 号 1 → 1 → IR→ IR → N→ N→ JC-(IR-1) 1 1 JC 列 号 方阵形式 等带宽形式 图4-8

40. 对比上图知： （1）等带宽存贮实际上是矩形数组存贮，矩形数组的行数和方阵行数相同，列数为半带宽数； （2）方阵存贮中第IR行第JC列（JC≥IR）元素Kir,jc在矩形数组中的行号与方阵存贮相同； （3）由于方阵存贮中的主对角线上元素在矩形数组中都在第一列，第IR行主对角线元素左移了（IR-1）列，且主对角线以右的所有元素都左移了（IR-1）列。所以，方阵存贮中第IR行第JC列（JC≥IR）元素Kir,jc在矩形数组中的列号是: 列号= JC-(IR-1)

41. 方阵存贮和矩形数组存贮地址关系 综上所述，可得到等带宽存贮下组集受约束结构总刚度矩阵的计算框图。与按方阵存贮时的计算框图相比，区别仅在于J循环应修改为如下形式：

42. JC ← LM [ J ] JC ← JC-IR+1 J←(I+1, 1, NDF2) K[IR,JC] ← K[IR,JC] + k [I,J ] IRR ← JC JCC ← IR-JC+1 K [IRR,JCC] ← K[IRR,JCC] + k [ J,I ] IR < JC JC = 0

43. 3、关于节点编号问题 与方阵存贮所需的存贮量相比，等带宽存贮大大减少了所需占用的计算机内存。例如，对于具有100个节点的平面刚架而言，总刚度矩阵共有300×300即9万个元素。按等带宽存贮时，假定其最大半带宽UBW = 18，则所需存贮量为300×18 = 5400个元素。 按等带宽方式存贮刚度矩阵也存在优化问题，这就是最大半带宽UBW最小。欲使最大半带宽UBW最小，必须注意节点编号方法，使直接联系的相邻节点的最大点号差最小。

44. 5.3 结构刚度矩阵的变带宽存贮 等带宽存贮虽然已经节省了不少内存，但认真研究半带宽内的元素，还有相当数量的零元素。在平衡方程求解过程中，有些零元素只增加运算工作量而对计算结果不产生影响。如果这些零元素不存、不算，更能节省内存和运算时间。这就是本章要介绍的变带宽存贮，也称一维数组存贮。 1、存贮方式 假定结构刚度矩阵[K]在方阵存贮中有如下形式：

45. 顶 线 UBW=4 对 称 图4-9 方阵形式的刚度矩阵[K]

46. （1）半带宽内零元素的两种情况 ① 零元素所在列的上方全是零元素。 ② 零元素所在列的上方还有非零元素。

47. （2）顶线以上零元素无须存贮 因为在求解平衡方程组 Kδ= P 的时候，它们根本不参加运算。例如采用自上而下地消去法求解或以后采用的其它解法都有是如此。 （3）顶线以下零元素须要存贮 例如，把第一个方程的若干倍加到第二个方程上去的时候，图4-9 中第四列顶线以下的那个零元素将成为非零元素，在之后的消元过程中要起作用，因而仍须存贮。

48. MAXA 1 2 4 6 10 12 16 18 22 图4-10一维数组[A]存贮刚度矩阵[K]

49. 变带宽存贮可采用按列存贮方式：从左到右，逐列存放，对每一列，先存主对角线元素，然后由下而上顺序存放，直到顶线下第一个元素为止。为避免混淆，我们把存贮[K]的一维数组称为[A]。如图4-10中，A(7)里存放的就是K34等。变带宽存贮可采用按列存贮方式：从左到右，逐列存放，对每一列，先存主对角线元素，然后由下而上顺序存放，直到顶线下第一个元素为止。为避免混淆，我们把存贮[K]的一维数组称为[A]。如图4-10中，A(7)里存放的就是K34等。

50. 为实现变带宽存贮，需要解决的标本问题是：刚度矩阵中i行j列的元素Kij在一维数组A中的地址是什么。为实现变带宽存贮，需要解决的标本问题是：刚度矩阵中i行j列的元素Kij在一维数组A中的地址是什么。 2、列高及其计算 把刚度矩阵[K]的第j列上方第1个非零元素的行号记为mj。例如，在图4-9中，m6 = 3， m8 = 5。