01

前言

众所周知，在使用有限元法进行静力分析时，大多数情况下求解这样一个矩阵方程，我们需要求解一个大的矩阵方程：

在RFEM的【静态分析设置】中，我们可以选择求解矩阵方程的方法。程序提供两种求解方法：直接法和迭代法。

图1 选择方程的求解方法

那么什么是直接法，什么是迭代法呢？这是求解线性方程组的两类方法。在与数值计算相关的书籍中，我们可以找到如下说明：

直接法可以通过有限步的算术运算得到方程组的精确解，但由于舍入误差的存在刚度矩阵计算公式，该方法只能得到线性方程组的近似解。

迭代法是利用一定的极限过程逐渐逼近线性方程组精确解的方法，适用于求解大型矩阵方程组。

在RFEM6程序中，默认使用“直接法”求解矩阵方程。如果二维有限元节点数超过100,000个或三维有限元节点数超过30,000个，建议使用迭代法。

图 2 查看网格统计信息

02

矩阵的条件数

如前所述，无论是用直接法还是迭代法求解矩阵方程，都是数值计算方法，得到的解是数值解，不是精确解。那么如何评价数值方法求解矩阵方程的误差呢？

为了说明这个问题，我们使用教科书[1]中的计算例子。

我们遵循机械问题方程组的格式：

考虑以下矩阵方程组：

它的精确解是x1=2，x2=0。当等式右边稍有变化时，如方程组变为：

此时，其精确解为x1=1，x2=1。可见，方程系统的一个微小扰动都会带来解的巨大变化。我们称此时的方程组为病态方程组，刚度矩阵K为病态矩阵。对于力学解，我们希望我们的刚度矩阵K不是病态矩阵，以免出现数值解不收敛而误差过大的情况。

为了评估K对这个小扰动的放大刚度矩阵计算公式，我们定义矩阵“条件数”cond(A)，矩阵条件数是矩阵范数||K||的范数。和逆矩阵 ||K^-1|| 产品：

矩阵的条件数越大，表示病态越严重。

03

评估错误

在大多数情况下，我们很难直接计算矩阵的“条件数”来评估误差，因为矩阵的条件数需要求矩阵逆的范数。

我们可以使用两种方法来判断方程的病态程度。从上面的例子可以看出，系数矩阵的两行距离很近，这个矩阵对应的行列式很小，方程近似奇异。在这种情况下，方程相对病态。

即第一种方法是求解刚度矩阵的行列式大小。当刚度矩阵 K 的行列式较小时，方程可能是病态的。完成静态分析后，我们可以在计算概览表中查看刚度矩阵行列式的大小。

图 3 计算概览

对于刚度矩阵的条件数，一般采用估计的方法而不是直接计算矩阵条件数，因为求解矩阵条件数往往比求解原矩阵方程更难。

如果结构网格和材料是均匀的，可以采用Fried[2]给出的估计方法，矩阵条件数的估计公式如下：

式中hmin为有限元单元的最小边长，2m为所求微分方程的阶数。该估算公式没有考虑网格划分不均匀、材料不均匀、元素个数、网格形状对矩阵条件数的影响。进一步的估计公式在[3][4]中给出：

式中，Emax为模型中材料的最大弹性模量，Emin为模型中材料的最小弹性模量。

hmax为有限元单元的最大边长，hmin为有限元单元的最小边长。

对于三角形单元，θ0 是最小内角。四面体单元的θ0是四面体单元的最小θ角。Ne为有限元单元总数，2m为微分方程的阶数。

图4 四面体单元θ角示意图

参考：

[1] 王绪成．有限元法基本原理及数值方法

[2] I.Fried，来自

[3] 邵秀敏．关于有限元矩阵的条件数[J]．应用数学杂志，1983（04）：413-419。

[4] 周树子．有限元矩阵条件数的估计[C]．湖南：湖南大学，1984：7。

过往回顾

下载：

手动的：

教程：丨丨丨丨丨丨丨丨丨丨丨丨丨|||||

在线课程：丨丨丨丨丨丨丨丨||