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