[有限元]从结点开始谈稀疏矩阵

结点在有限元方法中可以说是最简单的一部分了: 编号+空间坐标. 而单元, 仅从几何上说, 就是一系列结点按某种次序所组成的有序数组. 然而, 仅仅是简单的单元的结点数组信息(也就是先不谈形函数积分点等), 我们就可以确定一个非常重要的关系: 整体刚度矩阵的非零元素位置.

我们知道有限元中整体刚度矩阵的稀疏性是很高的, 也就是非零元素的个数远小于零元素的个数. 因此, 在Newton-Raphson方法中求解自由度的位移的修正量的时候, 我们得用稀疏矩阵去求解线性方程组, 否则, 无论是从存储, 还是浮点运算量, 都是不可忍受的.

在常见的有限元问题中, 稀疏矩阵有三类存储方式: 1) 对称稀疏矩阵, 用于材料模量是对称的情况; 2) 结构对称稀疏矩阵, 也就是非零元素的位置是对称的, K(i,j)和K(j,i)同时为零元素或者非零元素, 常见于材料模量是非对称的情况; 3) 非对称稀疏矩阵, 非零元素的位置不对称, 常见于完全耦合多物理场问题. 其实第三种情况是第二种情况的一种特殊形式: 把非对称稀疏矩阵中所有非零元素及其关于对角线对称位置的元素通通标记为非零元素, 即可用结构对称稀疏矩阵表达.

而对于有限元问题, 非零元素的位置可以直接从单元的结点数组信息获得. 我们称结点i和结点j是相邻的, 如果存在一个单元使得这个单元结点数组包含结点i和j. 那么, 结点i和j的所有自由度配对在整体刚度矩阵中是非零的. 举个例子, 二维问题中, 结点3的自由度在整体刚度矩阵中是5, 6, 结点7与结点3相邻, 其自由度编号是13, 14, 那么在整体刚度矩阵中, K(5,5), K(5,6), K(5,13), K(5,14), K(6,5), K(6,6), K(6,13), K(6,14), K(13,5), K(13,6), K(13,13), K(13,14), K(14,5), K(14,6), K(14,13), K(14,14)位置为非零元素.

以上是一些简单的关于有限元内稀疏矩阵的介绍. 注意这个稀疏矩阵只是在稳态以及隐式求解方法中需要, 而在显示求解方法中是不需要的.

下面说一些在具体问题中关于稀疏矩阵的一些思考.

第一个要说的是约束场变量边界条件的处理. 约束场变量的边界条件是说在某个分析步内, 将某个自由度的场变量(例如机械场里的位移, 温度场里面的问题)约束为零, 或者让其按照某个时间函数变化. 也就是说, 对于任意的分析步时间, 其场变量的值是可以立即获得, 而非在Newton-Raphson方法中求得. 在这种情况下, 不少书本以及实践都是将这个自由度在整体刚度矩阵中的所在行列清零, 然后将对角值设置成1, 变成 1 * x = r 的形式. 这样这个自由度的值和其余自由度没有任何的耦合关系.

这样处理诚然简单直接, 然而, 间接的增加了未知数的个数(当然计算量肯定不会随着未知数个数成比例的增加, 这些未知数的求解代价不会很高). 而且在某些稀疏矩阵存储方法下, 将某个自由度对应的行列清零也不是一个简单的事情.

其实我们直接将这些约束自由度完全从整体刚度矩阵中剔除. 我们可以建立一个这样的映射关系:

自由度编号: 1 2 3 4 5 6 7 8 ...

未知数编号: 0 1 2 3 0 4 5 0 ...

这个例子中, 自由度1, 5, 8为约束自由度, 而2, 3, 4, 6, 7分别为第1, 2, 3, 4, 5个未知数. 在从单元刚度矩阵组成整体刚度矩阵的时候, 我们可以先将单元刚度矩阵内自由度编号映射到未知数编号, 非零未知数配对才加入整体刚度矩阵.

第二个要说的是自由度对自由度约束的问题. 我们考虑这样一个简单的自由度i(成为从自由度)对自由度j(称为主自由度)的约束: u_i = u_j, 这里u是某个场变量. 举个例子, 在ABAQUS中, 我们给一个模型施加一个线载荷, 就可以做这样一个步骤: 1) 计算载荷总量, 也就是线载荷*长度; 2) 在所施加的线上选取一个结点, 在上面施加一个集中力载荷, 其大小等于载荷总量; 3) 将线上所有的其他结点和这个结点设置自由度对自由度的约束. 在这种情况下, 我们相当于在整体刚度矩阵当中增加了一些简单的约束方程. 那么, 这种情况下, 我们应该怎么处理呢.

首先要提的是, 这是一个最普遍意义下的约束方程下的一个最简单的情况. 这种约束可以是线性, 可以是非线性; 可以是一个自由度对一个自由度, 也可以是多个自由度对一个自由度. 最通常的情况下, 我们可以使用类似罚函数或者Lagangian乘子法等. 然而, 对于这里这个最简单的情况, 我们可以做一些简单的处理, 使得我们并不需要增加额外的负担. 回到上面提到的自由度-未知数编号的映射上, 假设我们现在令自由度2必须等于自由度7, 那么编号映射就得修改成

自由度编号: 1 2 3 4 5 6 7 8 ...

未知数编号: 0 4 1 2 0 3 4 0 ...

自由度2就不再是一个未知数. 然而, 这并没有结束, 因为这样的映射表会影响稀疏矩阵的非零元素的个数. 简单来说, 就是与从自由度相邻的结点, 也将与主自由度相邻. 按照这个准则, 去修改稀疏矩阵非零元素位置.

第三个要说的是扩充方程的问题. 在有些方法中, 例如arclength方法, 需要在整体刚度矩阵中增加一些与自由度有关的未知数, 也就是要扩充方程. 这类问题是有一定风险的. 因为哪怕扩充后的整体刚度矩阵仍然是稀疏的, 但经过类似LU分解之后, 可能结果就不再是稀疏的了.