第九章有限元线性方程组的解法课件.ppt

1、 前面，我们主要讨论了静力平衡问题的有限元格式。在确定了离散所需要的单元形式后，需要进行单元特性矩阵的计算，最后由单元特性矩阵集合而成的有限元求解方程组RK 是一组联立的线性代数方程组。这组方程在静力平衡问题中就是以结点位移为基本未知量的系统结点平衡方程。有限元求解的效率很大程度上取决于这组线性代数方程组的解法。有限元分析可以通过细分单元的网格来提高解的精度。因此，当有限元分析采用越来越多的单元离散模型来近似实际的结构时，线性联立方程组的阶数越来越高。曾经有相当一部分的研究工作是围绕如何有效地求解这组庞大的线性方程组的。 （9-1） 在线性静力分析中，解代数方程组的时间在整个解题时间中占很大的

2、比重。而在动力分析和非线性分析中这部分比重也是相当大的。若采用不适当的求解技术，不仅计算费用大量增加，更严重的是有可能导致求解过程的不稳定和求解的失败。 在有限单元法中，需求解的线性代数方程组（6-1）的系数矩阵K，即刚度矩阵，具有大型、对称、稀疏、带状分布以及正定、主元占优势的特点。在求解方程组时必须利用上述特点，以提高方程求解的效率，否则将无谓的提高计算费用。 线性联立方程组的解法可以分作两大类：直接解法和迭代法。 直接解法以高斯消元法为基础，求解效率高。在方程组的阶数不是特别高时（例如不超过10000阶），通常采用直接解法。当方程组的阶数过高时，由于计算机有效位数的限制，直接求解法中的舍

3、入误差，消元中有效位数的损失等将会影响方程求解的精度，此时可用迭代解法。 本章主要讨论几种常用的比较有效的直接解法，LU分解法和波前法。 有限元法中，线性代数方程组的系数矩阵是对称的，因此可只存储一个上三角（或下三角）矩阵。但是由于矩阵的稀疏性，仍然会发生零元素占绝大多数的情况。考虑到非零元素的分布呈带状特点，在计算机中系数矩阵的存储一般采用二维等带宽存储或一维变带宽存储，后者更为常用。 一维变带宽存储就是把变化的带宽内的元素按一定的顺序存储在一维数组中。按照解法可分为按行一维变带宽存储及按列一维变带宽存储。这里主要进行按列一维变带宽存储。887877676658565546454436343

4、3232214121100000000000000000KKKKKKKKKKKKKKKKKKK对称)18()19()16()20()17()12()21()13()10()14()11()6()15()7()4()8()5()2()9() 3() 1 (AAAAAAAAAAAAAAAAAAAAA 按列一维变带宽存储是按列依次存储元素，每列应从主对角线元素直至最高的非零元素，即该列中符号最小的非零元素为正，即图中突线所包括的元素。由图可以看出，这种存储对夹在非零元素内的零元素，如K24，K58等则必须存储。上图表示的是这些元素按列在一维数组中的排列。 把系数矩阵中的元素紧凑存储在一维数组中，必须

5、有辅助的数组帮助记录原系数矩阵的情况，例如对角元素的位置、每列元素的个数等。辅助数组M（n + 1)，用以记录主对角元素在一维数组中的位置。对于下图的一维数组，它的 （8 + 1) 数组是： M： 1 2 4 6 10 12 16 18 22 前n个数记录的是主对角元素的位置，最后一个数是一维数组刚度加1。 利用辅助数组M，除了知道各主元在一维数组中的位置以外，还可以用以计算每列元素的列高Ni，即每列元素的个数，以及每列元素的起始行号mi。Ni = M ( i + 1 ) M ( i ) （9-2） mi = i - Ni + 1 例如求第四列元素个数N4N4 = M ( 5 ) M ( 4

6、) = 10 6 = 4 求第六列元素的个数及非零元素的起始行号N6 = M ( 7 ) M ( 6 ) = 16 12 = 4 m6 = 6 4 + 1 = 3 有了辅助数组M后，可以找到一维数组中相应的元素进行各数组的求解。 一维变带宽存储是最节省内存的一种存储方法。 TLDLK1 nnninjnniiijiilllllllllllllllL21213332312221110 数学上已证明，对于对称正定的总刚度矩阵K，可以唯一地分解为其中（9-4）（9-3） nniillllD002211 TL 1D RLDLT1而 是L的转置， 是D的逆矩阵。将式（9-3）代入总刚方程后有（9-5） T

7、LD1 RL 若令则上式就变为这样，求解总刚方程组的通解，现在就转化为求解两个三角形方程组（9-6）和（9-7），即先由式（9-7）求出中间变量 ，然后以 为已知的右端项，由式（9-6）求得结点位移 。求解上面两个三角形方程组是很好解的，具体过程如下：（9-6） （9-7） 将式（9-3）的右端按矩阵乘法规则乘开后有nnnjnninijiinjnjKKKKKKKKKKKKKKKK2121222221111211一分解总刚度矩阵，求出L阵的lij11112112111111121121112111212111212211212121112111npnnppnpnpjpnjppjpnpnnnnpi

8、nppnpipjpijppjpipiiinnjjnjlllllllllllllllllllllllllllllllllllllllllll11iiKlijjpijppjpipKllll11jinjilllKjiljpppjpipijij, 2 , 1,011令对应元素相等，即得分解公式； 由此求的L阵的元素lij为（9-9） （9-8） 讨论： 1 从式（99）看出，在按行列由Kij计算lij时，计算完lij后，Kij就失去存在的作用，同时所用到lip、ljp和lpp排列顺序都在Kij之前，因此可将分解后得到的元素lij存贮在Kij单元中，即原来存贮K的内存单元，现在可用来存贮L矩阵，以减少对

9、内存贮量的要求。 2 由于这里只存贮下三角形带内元素，所以在利用式（99）由Kij计算lij时，求和号内各元素的列号应从第i行和第j列上第一个非零元素所在列号（i1和j1）中最大的列号开始。 3 从式（98）看出，在分解K时，每行的第一个非零元素其值保持不变，因此在分解总刚时，每行可从第二个非零元素的列号开始，这样lij的最后递推公式为iiijnilllKjiljjipppjpipijij, 2, 1, 3 , 20111),max(11（9-10） ), 2 , 1(11nillRiiikkikii 令式（97）两端对应元素相等得（9-11）二顺追赶求中间列向量iiii1i ), 3 , 2

10、(11nillRiiiikkikii 讨论： 1 从式（911）看出， 与Ri相对应，Ri只在求中间变量 时有用，求出 和Ri后就失去作用，因此可把求得的中间变量 存贮在Ri所用的内存单元中。 2 由于这里只存贮下三角形带内元素，因此求和号内的lik的列号k应从第i行的第一个非零元素的列号i1开始。 3 从上式还看出，当i = 1时， = R1，因此求中间量 时，可从第2行开始。考虑上面的原因，求中间列向量 的最后递推公式应为：（9-12） ) 1 , 2, 1(1nnilxliinijjjiii 将上面求得的中间列向量 作为常数项，代入式（96）即可求得未知结点位移列向量 。将（96）式左端

11、两矩阵相乘并令等式两端对应元素相等，得到如下的递推公式（9-13）iiiiii 讨论： 1因为 与 相对应，而且一旦求出 后， 就失去作用，因此把求得的 存贮在 的内存单元中，即存贮在结点荷载的内存单元中。 2. lij必须是带内元素，因此它的列号i必不小于该行的第一个非零元素的列号j1。三逆赶追求解未知结点位移列向量 采用高斯消去法和三角分解（LU分解）法时，方程一般都按结点自然编号排列。在有些情况下按自然顺序的带宽D很大，而中间夹有很多零元素，如复连通域的问题。造成工作三角形很大，这时可采用波前法。 波前（阵）法（Frontal Technique）是七十年代初提出并获得广泛应用的有限元法

12、所特有的一种程序方法，它是高斯消去法的一种特殊形式，也包括变量的消元和回代两个过程。与高斯消去法不同的是，总刚的组集和变量的消元是交错进行的，在内存中从不形成完整的总刚；组集时，变量是按单元和单元自由度的顺序进入总刚方程的；变量的消元是按变量成熟的先后次序进行的。下面通过一个例子说明它的基本思想。123456 图中连续体共划分4个单元，有6个结点。假定每一个结点有一个自由度，因此自由度编码与结点编码相同。单元扫描次序按单元编码顺序进行。 计算过程如下： 1. 按单元顺序扫描计算单元刚度矩阵及等效结点荷载列阵并送入内存进行组装。如首先扫描等元，进入内存的元素见图（a） 在K和R中集成尚未完毕的自

13、由度称为活动变量；集成完毕的自由度称为不活动变量。不活动变量可以作为主元行对其它非主元行的元素进行消元。 内存中存储刚度矩阵K元素的三角形称为波前三角形。在波前三角形中活动变量和主元按一定顺序构成波前。波前中变量的个数叫做波前数，记为W。 图(a)中，波前为2，4，5，波前数W = 3。 2检查哪些自由度已集成完毕，以集成完毕的自由度i作为主元对其它行列的元素进行消元修正。 图（b）中，自由度4已等成完毕，是不活动变量，现在作为主元，用 表示。主元行元素 ，不再变化，对其它行列元素进行消元修正。自由度 K P 2 扫描单元 4 5 波前 波前三角形 波前数W = 3 （a）紧凑2 2 5 3

14、集成完毕，为主元 紧凑波前区，将自由度5前移 主元行元素 被消元修正的其它元素（b） （c）扫描单元 5 442 2 扫描单元 3 自由度5，6同时集成完毕，分两步进行，每次一个主元（d）53 6 6 a1扫描单元 a2 a3全部扫描完毕依次消元，回代求解 （e）231iiiibibaiaRKKK,aba 3对其它行列元素进行消元修正后，主元已完成消元作用，将主元行有关元素Kij，Ri送入外存共有W+1个数。此后的紧凑波前区见图（c） 送入外存的元素是代数方程组中一个方程的系数和自由项 方程由哪些自由度 组成是由主元行 行未送出内存时的波前决定，因此在把主元行的元素送入外存前需记录有关信息：(

15、1)主元号A,(2)主元在波前中的位置I，(3)波前数W以后需要用这组信息来恢复波前。 送出主元行 行前,这组信息为A = 4；I = 2；W = 3。i4 4重复步骤13，将全部单元扫描完毕。全部扫描过程内存中的情况见图(a)(e)。 在最后一个单元扫描完成后，内存中全部自由度都已集成完毕，如图（e），此时，可直接消元后回代求解，得到最后内存中n个未知量 的解。i 在消元过程中得到一组A I W信息如下： A I W (送入外存元素) 4 2 3 （K4i i = 14） 5 2 4 （K5i i = 15） 6 3 3 （K6i i = 14）12364663362261KKKK 5回代求

16、解 按消元的顺序，由后向前逐个恢复波前，调入送到外存的元素，依次回代求解。 恢复波前需要利用A I W信息。在现有内存基础上根据A、I可将自由度A插在位置I上，然后根据W取出前n个自由度即为恢复的上一个波前。恢复一个波前就顺序有外存调入一组元素（后调出的元素先调入），一次求出方程组的解。 本例中最后内存中的波前为2，3，1，消元后解出未知量 ， 。利用最后一组AIW信息6，3，3恢复上一个波前：增加自由度6在第3的位置上，波前数为3，即2，3，1 2，3，6，1 2，3，6调入K6i ( i = 14 )，相当于得到方程236132654 ， 已经在上一个波前解得，由上列方程可解得 。然后再推

17、出前一个波前，用相同的方法求解一个新进入的自由度，由后向前直至全部求解完毕。过程如下所示 A I W 波前 调入 求解 （最后内存中） 2，3，1 ， ， 6，3，3 2，3，6 K6i（4个） 5，2，4 2，5，3，6 K5i（5个） 4，2，3 2，4，5 K4i（4个） 由上述解题过程可见，保留在内存中的波前区，包括波前三角形与自由项列阵，它的大小与结点的编码无关而与单元扫描的次序有关。为了确定不活动变量，需先记下每个结点的相关单元数，这只需事先把单元扫描一下即可解决。一般波前法扫描可以有两种次序： （1）按单元编码顺序扫描。上述例题就是如此。 （2）按自由度进入内存的顺序，保证先进入内存的自由度首先集成完毕作为主元退出波前。即按先进先出的次序进行扫描。 前者扫描方便，后者恢复波前较简单，波前区较小，但需多耗扫描时间。