1、 前面，我们主要讨论了静力平衡问题有限元格式。在确定了离散所需要单元形式后，需要进行单元特征矩阵计算，最终由单元特征矩阵集合而成有限元求解方程组 是一组联立线性代数方程组。这组方程在静力平衡问题中就是以结点位移为基本未知量系统结点平衡方程。有限元求解效率很大程度上取决于这组线性代数方程组解法。有限元分析能够经过细分单元网格来提升解精度。所以，当有限元分析采取越来越多单元离散模型来近似实际结构时，线性联立方程组阶数越来越高。曾经有相当一部分研究工作是围绕怎样有效地求解这组庞大线性方程组。第九章第九章 有限元线性方程组的解法有限元线性方程组的解法（9-1）第1页 在线性静力分析中，解代数方程组时间
2、在整个解题时间中占很大比重。而在动力分析和非线性分析中这部分比重也是相当大。若采取不适当求解技术，不但计算费用大量增加，更严重是有可能造成求解过程不稳定和求解失败。在有限单元法中，需求解线性代数方程组（6-1）系数矩阵K，即刚度矩阵，含有大型、对称、稀疏、带状分布以及正定、主元占优势特点。在求解方程组时必须利用上述特点，以提升方程求解效率，不然将无谓提升计算费用。线性联立方程组解法能够分作两大类：直接解法和迭代法。直接解法以高斯消元法为基础，求解效率高。在方程组阶数不是尤其高时（比如不超出10000阶），通常采取直接解法。当方程组阶数过高时，因为计算机有效位数限制，直接求解法中舍入误差，消元中
3、有效位数损失等将会影响方程求解精度，此时可用迭代解法。本章主要讨论几个惯用比较有效直接解法，LU分解法和波前法。第2页 有限元法中，线性代数方程组系数矩阵是对称，所以可只存放一个上三角（或下三角）矩阵。不过因为矩阵稀疏性，依然会发生零元素占绝大多数情况。考虑到非零元素分布呈带状特点，在计算机中系数矩阵存放普通采取二维等带宽存放或一维变带宽存放，后者更为惯用。一维变带宽存放就是把改变带宽内元素按一定次序存放在一维数组中。按照解法可分为按行一维变带宽存放及按列一维变带宽存放。这里主要进行按列一维变带宽存放。9-1 9-1 系数矩阵的一维变带宽存储系数矩阵的一维变带宽存储第3页对称第4页第5页 按列
4、一维变带宽存放是按列依次存放元素，每列应从主对角线元素直至最高非零元素，即该列中符号最小非零元素为正，即图中突线所包含元素。由图能够看出，这种存放对夹在非零元素内零元素，如K24，K58等则必须存放。上图表示是这些元素按列在一维数组中排列。把系数矩阵中元素紧凑存放在一维数组中，必须有辅助数组帮助统计原系数矩阵情况，比如对角元素位置、每列元素个数等。辅助数组M（n+1)，用以统计主对角元素在一维数组中位置。对于下列图一维数组，它（8+1)数组是：M：1 2 4 6 10 12 16 18 22 前n个数统计是主对角元素位置，最终一个数是一维数组刚度加1。第6页 利用辅助数组M，除了知道各主元在一
5、维数组中位置以外，还能够用以计算每列元素列高Ni，即每列元素个数，以及每列元素起始行号mi。Ni=M(i+1)M(i)（9-2）mi=i-Ni+1 比如求第四列元素个数N4N4=M(5)M(4)=10 6=4 求第六列元素个数及非零元素起始行号N6=M(7)M(6)=16 12=4 m6=6 4+1=3 有了辅助数组M后，能够找到一维数组中对应元素进行各数组求解。一维变带宽存放是最节约内存一个存放方法。第7页 数学上已证实，对于对称正定总刚度矩阵K，能够唯一地分解为其中（9-4）（9-3）8-28-2修改的平方根法（即修改的平方根法（即LULU分解法）分解法）第8页而 是L转置，是D逆矩阵。将
6、式（9-3）代入总刚方程后有（9-5）第9页若令则上式就变为这么，求解总刚方程组通解，现在就转化为求解两个三角形方程组（9-6）和（9-7），即先由式（9-7）求出中间变量 ，然后以 为已知右端项，由式（9-6）求得结点位移 。求解上面两个三角形方程组是很好解，详细过程以下：（9-6）（9-7）第10页 将式（9-3）右端按矩阵乘法规则乘开后有一分解总刚度矩阵，求出L阵lij第11页第12页令对应元素相等，即得分解公式；由此求L阵元素lij为（9-9）（9-8）第13页 讨论：1 从式（99）看出，在按行列由Kij计算lij时，计算完lij后，Kij就失去存在作用，同时所用到lip、ljp和l
7、pp排列次序都在Kij之前，所以可将分解后得到元素lij存贮在Kij单元中，即原来存贮K内存单元，现在可用来存贮L矩阵，以降低对内存贮量要求。2 因为这里只存贮下三角形带内元素，所以在利用式（99）由Kij计算lij时，求和号内各元素列号应从第i行和第j列上第一个非零元素所在列号（i1和j1）中最大列号开始。3 从式（98）看出，在分解K时，每行第一个非零元素其值保持不变，所以在分解总刚时，每行可从第二个非零元素列号开始，这么lij最终递推公式为第14页（9-10）令式（97）两端对应元素相等得（9-11）二顺追赶求中间列向量第15页 讨论：1 从式（911）看出，与Ri相对应，Ri只在求中间
8、变量 时有用，求出 和Ri后就失去作用，所以可把求得中间变量 存贮在Ri所用内存单元中。2 因为这里只存贮下三角形带内元素，所以求和号内lik列号k应从第i行第一个非零元素列号i1开始。3 从上式还看出，当i=1时，=R1，所以求中间量 时，可从第2行开始。考虑上面原因，求中间列向量 最终递推公式应为：（9-12）第16页 将上面求得中间列向量 作为常数项，代入式（96）即可求得未知结点位移列向量 。将（96）式左端两矩阵相乘并令等式两端对应元素相等，得到以下递推公式（9-13）讨论：1因为 与 相对应，而且一旦求出 后，就失去作用，所以把求得 存贮在 内存单元中，即存贮在结点荷载内存单元中。
9、2.lij必须是带内元素，所以它列号i必大于该行第一个非零元素列号j1。三逆赶追求解未知结点位移列向量第17页 采取高斯消去法和三角分解（LU分解）法时，方程普通都按结点自然编号排列。在有些情况下按自然次序带宽D很大，而中间夹有很多零元素，如复连通域问题。造成工作三角形很大，这时可采取波前法。波前（阵）法（Frontal Technique）是七十年代初提出并取得广泛应用有限元法所特有一个程序方法，它是高斯消去法一个特殊形式，也包含变量消元和回代两个过程。与高斯消去法不一样是，总刚组集和变量消元是交织进行，在内存中从不形成完整总刚；组集时，变量是按单元和单元自由度次序进入总刚方程；变量消元是按
10、变量成熟先后次序进行。下面经过一个例子说明它基本思想。83 83 求解总刚分担的波前法求解总刚分担的波前法第18页123456 图中连续体共划分4个单元，有6个结点。假定每一个结点有一个自由度，所以自由度编码与结点编码相同。单元扫描次序按单元编码次序进行。第19页 计算过程以下：1.按单元次序扫描计算单元刚度矩阵及等效结点荷载列阵并送入内存进行组装。如首先扫描等元，进入内存元素见图（a）在K和R中集成还未完成自由度称为活动变量；集成完成自由度称为不活动变量。不活动变量能够作为主元行对其它非主元行元素进行消元。内存中存放刚度矩阵K元素三角形称为波前三角形。在波前三角形中活动变量和主元按一定次序组
11、成波前。波前中变量个数叫做波前数，记为W。图(a)中，波前为2，4，5，波前数W=3。第20页 2检验哪些自由度已集成完成，以集成完成自由度i作为主元对其它行列元素进行消元修正。图（b）中，自由度4已等成完成，是不活动变量，现在作为主元，用 表示。主元行元素 ，不再改变，对其它行列元素进行消元修正。自由度 K P 2 扫描单元 4 5 波前 波前三角形 波前数W=3 （a）第21页紧凑2 2 5 3 集成完成，为主元 紧凑波前区，将自由度5前移 主元行元素 被消元修正其它元素（b）（c）扫描单元 5 44第22页2 2 扫描单元 3 自由度5，6同时集成完成，分两步进行，每次一个主元（d）53
12、 6 6第23页 a1扫描单元 a2 a3全部扫描完成依次消元，回代求解 （e）231第24页 3对其它行列元素进行消元修正后，主元已完成消元作用，将主元行相关元素Kij，Ri送入外存共有W+1个数。今后紧凑波前区见图（c）送入外存元素是代数方程组中一个方程系数和自由项 方程由哪些自由度 组成是由主元行 行未送出内存时波前决定，所以在把主元行元素送入外存前需统计相关信息：(1)主元号A,(2)主元在波前中位置I，(3)波前数W以后需要用这组信息来恢复波前。送出主元行 行前,这组信息为A=4；I=2；W=3。i4第25页 4重复步骤13，将全部单元扫描完成。全部扫描过程内存中情况见图(a)(e)
13、。在最终一个单元扫描完成后，内存中全部自由度都已集成完成，如图（e），此时，可直接消元后回代求解，得到最终内存中n个未知量 解。在消元过程中得到一组A I W信息以下：A I W (送入外存元素)4 2 3 （K4i i=14）5 2 4 （K5i i=15）6 3 3 （K6i i=14）第26页 5回代求解 按消元次序，由后向前逐一恢复波前，调入送到外存元素，依次回代求解。恢复波前需要利用A I W信息。在现有内存基础上依据A、I可将自由度A插在位置I上，然后依据W取出前n个自由度即为恢复上一个波前。恢复一个波前就次序有外存调入一组元素（后调出元素先调入），一次求出方程组解。本例中最终内存
14、中波前为2，3，1，消元后解出未知量 ，。利用最终一组AIW信息6，3，3恢复上一个波前：增加自由度6在第3位置上，波前数为3，即2，3，1 2，3，6，1 2，3，6调入K6i(i=14)，相当于得到方程第27页 ，已经在上一个波前解得，由上列方程可解得 。然后再推出前一个波前，用相同方法求解一个新进入自由度，由后向前直至全部求解完成。过程以下所表示 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个）第28页 由上述解题过程可见，保留在内存中波前区，包含波前三角形与自由项列阵，它大小与结点编码无关而与单元扫描次序相关。为了确定不活动变量，需先记下每个结点相关单元数，这只需事先把单元扫描一下即可处理。普通波前法扫描能够有两种次序：（1）按单元编码次序扫描。上述例题就是如此。（2）按自由度进入内存次序，确保先进入内存自由度首先集成完成作为主元退出波前。即按先进先出次序进行扫描。前者扫描方便，后者恢复波前较简单，波前区较小，但需多耗扫描时间。第29页
