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

《第四章结构刚度矩阵的存贮方式和组集程序.ppt》由会员分享，可在线阅读，更多相关《第四章结构刚度矩阵的存贮方式和组集程序.ppt（64页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、第四章 结构刚度矩阵的存贮方式和组集程序,前面几章，我们介绍了有限元法的基本原理，较深入地讨论了单元分析和系统分析理论。有限元法的重要价值在于方便编写程序，借助于计算机完成全部计算工作。,在整个编程和计算中，影响最大的是结构刚度矩阵的存贮和组集。本章就对这一问题进行专门讨论。,结构刚度矩阵的组集方法、程序和它的存贮方式密切相关，还依赖于约束条件的处理方式。本章介绍三种存贮方式：方阵存贮、等带宽存贮和一维数组存贮。,本章讨论的结构刚度矩阵存贮方式和组集的基本内容，适合于各类结构，包括杆系结构、弹性连续结构等。适合于结构静力分析，也是结构动力分析以及稳定性分析的重要基础。,约束条件采用第三章思想处

2、理：从无约束的结构刚度矩阵中删去与受约束位移号对应的行和列，再将矩阵压缩排列成nn阶方阵。实际操作时，根本就不去存贮这些行和列。,下面，我们就从方阵存贮方式开始。,为叙述简便起见，本章中所用“单元刚度矩阵”一词都指结构坐标单元刚度矩阵。只有特别需要时才加“结构坐标”定语。,5.1 结构刚度矩阵的方阵存贮,（3-8）,由第三章式（3-8）知,结构刚度矩阵K等于所有膨胀后的单元刚度矩阵k相累加。下面详细考察如何在程序中实现这个操作。,1、单元刚度矩阵元素的两类地址,单刚地址和总刚地址。,（1）单刚地址 单刚地址是指单元刚度矩阵元素在膨胀前的单元刚度矩阵中的地址，即行、列号。式（4-1）矩阵中每个元

3、素的两个下标就是该元素的单刚地址。第一个下标指行号，第二个下标指列号。例如k15的说明它是矩阵中第1行第5列位置的元素。,（4-1）,（2）总刚地址 总刚地址是指单元刚度矩阵元素在膨胀后的单元刚度矩阵中的地址（行、列号）。这个地址与单元的节点号有关，而节点号是离散结构时编好的，还与节点自由度有关。,设式（4-1）是图4-1所示梁单元的单元刚度矩阵。该单元两个节点的节点号为i、j，节点自由度是3。,单元刚度矩阵元素的总刚地址应为：i1、i2、i3、j1、j2、j3，如式（4-2）。,有,（4-3）,图4-2,如果式（4-1）是图4-2所示平面问题三角形单元的单元刚度矩阵。该单元三个节点的节点号为

4、i、j、m，节点自由度是2。,单元刚度矩阵元素的总刚地址应为：i1、i2、j1、j2、m1、m2，如式（4-4）、（4-5）。,有,（4-4）,（4-5）,（3）两类地址的连接数组,建立一个连接数组供计算机识别两类地址，以便正确地从单元刚度矩阵中取出元素，并放入结构刚度矩阵之中。,给连接数组取个名，例如：LM（NDF2），它是一维数组，数组长度NDF2是单元自由度的最大值。,先不考虑约束情况：,对于平面梁单元组成的结构,对于平面问题三角形单元组成的结构,从式（4-6）、（4-7）知，建立连接数组应先知道单元的节点号i、j、m 。这些数据应作为结构计算的初始数据在程序开始后即读入的，因此是已知的

5、。通常，这些数据读入后，按照单元号顺序由小到大依次存放在若干个一维数组中，例如，数组IO、JO、MO分别存放i、j、m节点号；数组IO、JO、MO的长度为单元总数。当然也可用一个多维数组存放它们。,（4）建立连接数组的条件,现在我们指出一个明显而重要的事实。所谓单元刚度矩阵的膨胀只不过为了说明问题的概念。事实上，执行时并不去膨胀单元刚度矩阵。所谓把所有单元膨胀后的单元刚度矩阵累加起来，即算式 的实现操作方法非常简单：,只要把每个单元的单元刚度矩阵k的每个元素按照它的总刚地址送入总刚度矩阵K，并简单地累加起来就组集起了总刚度矩阵K。,2、组集无约束结构刚度矩阵,（1）方法,（2）步骤,组集无约束

6、结构总刚度矩阵的具体步骤如下： 矩阵K充零 即把K的N1N1个元素全部置为零元素。N1=NJ（节点总数）NDF（一个节点自由度数） 从第1号单元开始，依次对所有单元执行下述操作：(下面符号“e”代表被执行的单元号) (a) 计算e号单元单元坐标系下的单元刚度矩阵ke。 (b) 计算e号单元的坐标变换矩阵R。 (c) 形成e号单元结构坐标系下的单元刚度矩阵k。计算公式为: k=RTkeR,(d) 建立e号单元的连接数组LM (e) 把e号单元结构坐标单元刚度矩阵k的元素按连接数组给出的总刚地址送入结构刚度矩阵K并进行累加。赋值语句为： K(i,j)=K(i,j)+k(p,q) （4-8） K(i

7、,j) 结构刚度矩阵中第i行j列元素 k(p,q) 单元刚度矩阵中第p行q列元素 p, q 单刚地址 i, j 总刚地址,为了节省存贮、节省运算时间，执行第(e)步时应利用单元刚度矩阵和结构刚度矩阵的对称性。用式（4-8）只须将单元刚度矩阵的上三角部分（含对角线元素）装入结构刚度矩阵中，形成结构刚度矩阵的上三角部分。,执行命令 K( j, i )=K( i, j ) ( ij )即得结构刚度矩阵K的下三角部分元素。,结构刚度矩阵中的上三角元素,执行上述第步时，对每个单元皆重复同样规律的计算，因而应编写独立的子程序。,（3）程序框图,3、组集约束结构总刚度矩阵,组集起无约束结构的总刚度矩阵K之后

8、，形成受约束结构的结构刚度矩阵Kff的原则是很简单的。例如，如果第L号位移分量受有刚性约束，即L=0。那么，只要把无约束结构刚度矩阵K的第L行和第L列删去即可。对所有受刚性约束的位移分量都作这样处理后，就得到约束结构的结构刚度矩阵Kff。,所谓约束是指刚性约束。弹性约束可按弹性单元处理，而在刚性支承处引入约束。,假定所讨论的结构受刚性支承约束的节点共有NRJ个。我们用一维数组 KRJ(NRJ)来存贮这些受刚性支承约束的节点的号码。 每个受约束的节点，它的各个自由度可能全受到约束，也可能只有一个自由度被约束。为了表明这种情况，还需一个二维数组:,（1）约束信息数组,KRL(NDF,NRJ),下面

9、我们就按照这个原则来组集Kff。但不是在形成K之后再从其中删去一些行和列，而是根本不去存贮应该删去的那些行和列。,例 图4-3所示框架结构NRJ=3,数组KRL的第i列，说明数组KRJ中的第i个节点的每个自由度是否受约束。我们规定：元素为1受约束，0表示不受约束。,1号节点1、 2受约束， 3不受约束； 2号节点3个位移均受约束；3号节点只2受约束。,数组KRJ、KRL给出了结构的全部约束信息。这些信息应作为初始数据交由计算机读入的。,（2）节点位移信息数组,节点位移信息数组可命名为ID(NDF,NJ)。NJ是节点总数。数组ID是一个NDF行, NJ列的二维数组。,ID数组先用来表明每个节点（

10、无论它是否受约束）的各个自由是否受到约束。对图4-3情形，NDF=3，NJ=6，ID数组的内容是：,（4-9）,1,2,3,4,6,5,这一步可借助于KRL数组实现。,再对ID数组进行如下改造：从第1列开始，自上而下检查，遇1改0，遇0则用上一个非0元素加1；然后，对第2列、第3列、依次按同样办法检查处理，完成后，得最终形式的ID数组。仍以图4-3为例，对式（5-9）改造后，得:,（4-10）,0,0,1,0,0,0,2,0,3,4,5,6,7,8,9,10,11,12,a) 如果 ID（i, j ) = 0则表明j号节点第i个自由度受有约束。,b) 如果 ID( i, j ) 0则j号节点第

11、i个自由度不受约束。并且， j号节点第i个位移分量在非约束节点位移列向量 f的序号就是: ID( i, j ),改造后的ID数组的性质是:,ID数组形成框图,受约束节点总数,受约束节点号,受约束节点信息,组集受约束结构总刚度矩阵Kff和组集无约束结构总刚度矩阵K的区别在于前者引入了约束，而后者没有。因此，组集受约束结构总刚度矩阵Kff的步骤与组集无约束结构总刚度矩阵K的步骤完全相同。只须在组集无约束结构总刚度矩阵K的方法中，做以下三点修改，以反应约束的引入。,（3） 受约束结构总刚度矩阵的组集, 结构刚度矩阵的体积从n1n1改为nn n1不考虑约束的节点总自由度 n1=NDFNJ, 根据ID数

12、组形成连接数组 在单元刚度k的总刚地址中只要为0，则根本不把k中的相应行和列送入结构刚度矩阵。,执行上述过程，得到的就是受约束结构刚度矩阵Kff。,n删去约束后，结构的节点自由度总数 n=n1-nr nr受约束节点的自由度数,（4）程序框图,4.2 结构刚度矩阵的等带宽存贮,1、结构刚度矩阵的稀疏性,在第3章里，我们曾经指出，结构刚度矩阵是具有大量零元素的稀疏矩阵。,现在对这一性质作进一步的考察。,在对结构系统进行平衡分析时知，如果i节点发生单位位移，只有与节点i有直接联系的那些节点才会受到影响，产生节点力。反之， i节点是否产生节点力也只会受与它直接相连的节点的影响。 设结构中共有100个节

13、点，其中第i（i = 4）号节点及其与i号节点直接联系的节点如图4-4所示。,则结构刚度矩阵中的3i-2, 3i-1, 3i 行的形式为：,上图中每个小方块是一个三阶子块阵，所示出的与节点i相应的100个三阶子块中，只有4个非零块，其余96个皆为零元素组成的三阶子块。通常，对于任何一个节点，与它直接联系的节点个数和问题中的节点总数比较起来总是小得多。因而，呈现出结构刚度矩阵的高度稀疏性。,图4-5,4,由于刚度矩阵为对称矩阵，只要对节点适当编号，就能使刚度矩阵的非零元素聚集在主对角线两侧，呈现为图4-6所示带状矩阵。,一般，刚度矩阵各行的非零元素个数和排列是不同的。对第i行而言，从主对角线元素

14、Kii算起，直到它最右边的一个非零元素为止，其间所有元素（包括可能有的零元素）的个数，称为该行的右带宽。以图4-5为例，第10行的右带宽为12，即包括 K10，10，K10，11，K10，21。这12个元素。,把相互连接的两个节点i和j的号码差j i称为相邻点号差，当ij时，从Kjj到Kji也共有j i+ 1个三阶子块。因而：,图中所示三行中最大的右带宽为： （j i+ 1 ）* 3,图4-7,4,刚度矩阵各行右带宽中之最大者，称为最大右带宽。因刚度矩阵为对称矩阵，最大右带宽和最大左带宽是相同的，故简称最大右带宽为最大半带宽，记为UBW（Upper Band Width ）。,当结构的节点编号

15、确定后，刚度矩阵的最大半带宽随之确定。找出所有相邻点号差中的最大者，记为rmax，则刚度矩阵最大半带宽为 UBW =（1 + rmax ）* NDF （4-9）,其中，NDF为节点自由度。对平面刚架，NDF = 3。对弹性平面问题，NDF=2。,2、结构刚度矩阵的等带宽存贮,鉴于工程结构都是几何不变结构，从现在起，我们只讨论受约束结构的刚度矩阵。为书写简便起见，用K代表受约束结构的刚度矩阵，即以前的Kff。,因为刚度矩阵是对称的，只须存贮它的一半即可。约定：只存贮K的上三角部分。带区之外全是零元素，自然不必存贮，因此只需存贮包括主对角元素在内的上半带宽区域内的元素。这样，大大节约了需占用的计算

16、机内存。称这种存贮方式为等带宽存贮。,设N为独立结点位移分量总数，则原始形式的刚度矩阵K是一个N行N列的方阵。所谓等带宽存贮，就是把K的上半带宽区域内全部元素用一个N行UBW列的矩形数组予以存贮：,方阵形式,等带宽形式,图4-8,对比上图知： （1）等带宽存贮实际上是矩形数组存贮，矩形数组的行数和方阵行数相同，列数为半带宽数； （2）方阵存贮中第IR行第JC列（JCIR）元素Kir,jc在矩形数组中的行号与方阵存贮相同； （3）由于方阵存贮中的主对角线上元素在矩形数组中都在第一列，第IR行主对角线元素左移了（IR-1）列，且主对角线以右的所有元素都左移了（IR-1）列。所以，方阵存贮中第IR行

17、第JC列（JCIR）元素Kir,jc在矩形数组中的列号是: 列号= JC-(IR-1),方阵存贮和矩形数组存贮地址关系,综上所述，可得到等带宽存贮下组集受约束结构总刚度矩阵的计算框图。与按方阵存贮时的计算框图相比，区别仅在于J循环应修改为如下形式：,3、关于节点编号问题,与方阵存贮所需的存贮量相比，等带宽存贮大大减少了所需占用的计算机内存。例如，对于具有100个节点的平面刚架而言，总刚度矩阵共有300300即9万个元素。按等带宽存贮时，假定其最大半带宽UBW = 18，则所需存贮量为30018 = 5400个元素。,按等带宽方式存贮刚度矩阵也存在优化问题，这就是最大半带宽UBW最小。欲使最大半

18、带宽UBW最小，必须注意节点编号方法，使直接联系的相邻节点的最大点号差最小。,5.3 结构刚度矩阵的变带宽存贮,等带宽存贮虽然已经节省了不少内存，但认真研究半带宽内的元素，还有相当数量的零元素。在平衡方程求解过程中，有些零元素只增加运算工作量而对计算结果不产生影响。如果这些零元素不存、不算，更能节省内存和运算时间。这就是本章要介绍的变带宽存贮，也称一维数组存贮。,1、存贮方式,假定结构刚度矩阵K在方阵存贮中有如下形式：,对 称,图4-9 方阵形式的刚度矩阵K,（1）半带宽内零元素的两种情况, 零元素所在列的上方全是零元素。 零元素所在列的上方还有非零元素。,例如，把第一个方程的若干倍加到第二个

19、方程上去的时候，图4-9 中第四列顶线以下的那个零元素将成为非零元素，在之后的消元过程中要起作用，因而仍须存贮。,（3）顶线以下零元素须要存贮,（2）顶线以上零元素无须存贮,因为在求解平衡方程组 K= P 的时候，它们根本不参加运算。例如采用自上而下地消去法求解或以后采用的其它解法都有是如此。,图4-10 一维数组A存贮刚度矩阵K,变带宽存贮可采用按列存贮方式：从左到右，逐列存放，对每一列，先存主对角线元素，然后由下而上顺序存放，直到顶线下第一个元素为止。为避免混淆，我们把存贮K的一维数组称为A。如图4-10中，A(7)里存放的就是K34等。,2、列高及其计算,把刚度矩阵K的第j列上方第1个非

20、零元素的行号记为mj。例如，在图4-9中，m6 = 3， m8 = 5。,为实现变带宽存贮，需要解决的标本问题是：刚度矩阵中i行j列的元素Kij在一维数组A中的地址是什么。,从第j列的主对角线元素起到该列上方第一个非零元素为止，所含元素的个数称为第j列的列高，记为hj 。有 hj = j - mj + 1显然，hj其实也就是第j行的左带宽，因而必有 UBW= max(hj) j=1,2, ,N,利用4.1所述中的节点位移信息数组（ID），可容易地确定刚度矩阵K任何一列的列高。,例如，对图4-11所示框架结构，利用ID数组可得各单元的连接数组LM：（假定小号为i）,图5-11,1号单元：LM=0

21、, 0, 1, 0, 0, 2 2号单元：LM=0, 0, 2, 3, 4, 5 3号单元：LM=3, 4, 5, 6, 7, 8 4单元：LM=0, 0, 1, 6, 7, 8,确定任何一列例如第7列列高的办法是：从1号单元起，对所有单元逐个进行检查。其中，与7号位移分量同在一个连接数组中的最小非零号码就是m7。显然有 m7 = 1,于是第7列的列高为 h7 = 7 1 + 1 = 7用同样方法即可确定其余各列的列高。,计算列高的程序框图,3、主对角线元素地址,为了计算Ki,j在一维数组A排列中的地址，除数组A之外，再设置一个一维数组MAXA来存放K的主对角线元素在一维数组A中的地址。数组M

22、AXA的长度是K的行或列数加1（N+1）。对于图4-10情形，数组MAXA的内容给予该图最右边。元素MAXA(N+1)用来存放假想的第N+1个主对角元的地址。,知道了K的各列列高之后，K的任何一个主对角元在一维数组A中的地址是容易确定的：K矩阵第j列主对角线元素在一维数组A中的地址等于前(j-1)列的列高之和加1，即,MAXA(j) = h1 + hj-2 - hj-1 + 1 =（h1 + hj-2 + 1）+ hj-1 = MAXA（j - 1）+ hj-1 因为永远有 MAXA(1)= 1， MAXA(2)= 2,故计算主元地址的公式可写为 MAXA(j+1)= MAXA(j)+ hj

23、（5-12）式中， j = 2，3，N； hj刚度矩阵K第j列的列高。,因为假想的第N+1个主元地址MAXA(N+1)等于K的所有N列列高之和再加1，所以，一维数组A的总长度(S) ，即刚度矩阵K按变带宽存贮的总存贮量 S = MAXA(N+1)- MAXA(1) （5-13）,主对角线元素地址计算框图,4、K矩阵中Ki,j在一维数组A中的地址,记Ki,j在一维数组A中的地址为AIJ。则由图4-12知，,A I J = MAXA J + J I （4-14） 其中，I = mj，mj+1，J。,图4-12,j列,MAXA(J),AIJ,5.4 组集结构刚度矩阵小结,1、压缩存贮 结构刚度矩阵是

24、具有高度稀疏性的对称矩阵。元素的数量很大，采用压缩存贮形式是绝对必要的。我们介绍了等带宽存贮和变带宽存贮这两种最常用的存贮方式。它们都适用于线性方程组的直接解法，也适用于迭代解法。,2、只存贮刚度矩阵的上三角部分 无论是等带宽存贮还是变带宽存贮，都利用结构刚度矩阵的对称性，都只存贮了刚度矩阵的上三角部分，即主对角线元素以及上半带区（或顶线以下）内的元素。,3、不去存贮与刚性支承约束相应的行和列 处理约束时，我们采用的处理方式是根本不去存贮与刚性支承约束相应的行和列。为此，须根据给定的支承约束条件形成ID数组。根据ID数组可得到单元连接数组LM，它给出了单元两端结点位移分量在独立结点位移列矢量中

25、的序号。其中，与受刚性支承约束的位移分量相应的序号为零。,4、节点编号应得当 无论采用等带宽存贮或变带宽存贮，都要求结点号码编得恰当，使得刚度矩阵的非零元素都聚集在主对角线附近狭窄的带形区内。,对于大型结构，应该特别注意选取恰当的结点编号方式。这是因为总存贮量大体上与带宽成正比的原故（对等带宽存贮，则准确地与最大半带宽成正比）。,4、程序中应通过试算验证 编好的程序，必须经过试算，以检查刚度矩阵是否被正确地形成。通常，可采用简单但具有典型性的算例进行计试算，输出形成后的刚度矩阵。也可以在组集刚度矩阵的程序段内设置一些输出语句，输出少量信息以便检查。导致刚度矩阵未被正确形成的因素是多方面的，可能是子程序编写有错，也可能计算程序完全正确而输入数据本身有错，例如支承约束条件填写有错等等。,总之，计算机铁面无情，不允许在任何环节上发生任何差错。另一方面，它也是很公正的，对于所有认真而细心的使用者，都迅速准确地提供所要求的计算结果。,本 章 完,