关于有限元应力结果精度的几点讨论

作为一名汽车结构 CAE 工程师，使用有限元法进行应力分析是我多年的日常工作。但是关于有限元应力结果的一些技术点一直未能吃透，查到的相关文献也不多，文献内容也存在很多不尽不实之处。最近一段时间，抽空重新读了一遍王勖成老师的《有限单元法》，也查阅了一些其它资料，简单整理出了一些关于应力结果的笔记，再加上一点自己的理解，与大家分享交流。

1. 有限元位移解的下限性

利用最小位能原理求得的位移近似解所对应的弹性变形能是精确解变形能的下限，即位移近似解在总体上偏小，模型偏于刚硬。利用最小余能原理得到的应力近似解所对应的弹性余能是精确解余能的上限，即应力近似解在总体上偏大，结构的计算模型偏于柔软。

我们常用的结构单元大都是位移元，以位移为未知量，基于最小位能原理建立。位移元得到的位移解具有下限性质，在给定的载荷之下，计算模型的变形比实际要小。当单元网格分割得越来越细时，位移数值解将由下方收敛于精确解，即得到真实解的下界。

需要注意的是，位移数值解并不是在每个点上都小于精确解，数值解只是总体上小于精确解，更准确的说法应该是外载荷在数值解上做的功小于在精确解上做的功，即

$$\tilde{u}^TP \le u^TP$$

上式中，ũ和 u 分别为位移数值解和精确解，P 为外载荷。

位移解的下界性质从物理上非常好理解。连续体具有无限多个自由度，可以有无数种复杂的变形模式。划分单元后，位移场用有限个结点的位移和对应的形函数来表示，即用一些简单变形模式来逼近实际变形，如图 1 所示。这就意味着连续体的变形受到了约束和限制，即刚度较实际增加了。由刚度方程可知，在外力相同的情况下，所求得的位移近似解将在总体上偏小。

图 1 用特定的形函数逼近准确解

2. 应力结果的精度

使用位移有限元法进行结构分析时，未知的场函数是结构位移。利用最小位能原理建立的求解方程是系统的平衡方程，求解方程得到的是各结点的位移，但实际工程问题往往更关注结构应力分布。

位移有限元法求解应力的基本步骤如下：

• 引入位移边界条件；
• 求解方程组得到各结点的位移 a；

$$Ka=P$$

公式 1

• 根据单元各结点位移，通过导数运算求得应变和应力。

$$\varepsilon=Lu=Ba^e$$

公式 2

$$\sigma=D\varepsilon=DBa^e$$

公式 3

例如对于平面问题，

$$\varepsilon=\lbrack\varepsilon_x \varepsilon_y \varepsilon_z\rbrack$$

$$ L=\begin{bmatrix}{\partial}/{\partial{x} } & 0 & {\partial}/{\partial{y}}\\0 & {\partial}/{\partial{x} } & {\partial}/{\partial{y} }\end{bmatrix}^T$$

$$N=\begin{bmatrix} N_1& 0& N_ 2& 0& N_ 3&0&……\\ 0 &N_1&0&N_2&0&N_3 &…… \end{bmatrix}$$

$$\large B=\large L \large N$$

公式 4

应变矩阵 B 是对插值函数求导得到的矩阵，每求导一次，插值多项式的次数就降低一次。所以通过导数运算得到的应力解的精度较位移解降低了一阶。例如线性单元的应力是近似均布的，二次单元的应力是近似线性分布的，其精度都比位移解低（虽然应力分布函数还包含了一些高次项，但这些高次项都是非完全项，并不能提高应力精度）。

有限元应力解的近似性表现在：

• 单元内部一般不满足平衡方程；
• 单元与单元交界面上应力一般不连续；
• 力的边界上一般也不满足力边界条件。

只有单元尺寸无限趋近于 0 时，即自由度数趋近无穷时，才能精确满足平衡方程、力边界条件和单元交界面上的应力连续性。在单元数量有限时，这些条件只能近似满足，除非实际应力变化的阶次等于或低于所用单元的应力分布函数阶次。

3. 应力解的震荡性质

位移有限元法从力学上解释是求位移变分所引起的应变能为极小值的问题，从数学上解释是求解应力近似解与精确解差值的加权最小二乘问题。

与位移结果不同，位移元的应力结果并没有下限性质。应力结果是精确应力在加权最小二乘意义上的近似解，应力近似解必然在精确解上下震荡；并且在某些点上，近似解恰好等于精确解，即单元内存在最佳应力点，如图 2 所示。

图 2 有限元应力解的震荡性质

我们在有限元强度分析中，经常发现在应力集中部位的应力解低于精确解，所以有人误以为位移有限元法的应力解也有下限性质。但真正原因是单元数目有限，且应力近似解的阶次低，无法正确描述应力的剧烈变化。对于高应力梯度区域，有限元解通常给出的是一个比实际平滑的结果。如果结构上某部位的应力突然上升，该部位的应力解将低于精确解，但如果某个部位应力突然下降，该处的应力解反而会高于精确解。

位移元的应力解没有下限性质从物理上也容易解释，有限元相当于在结构内部增加了约束，即增加了整体刚度。在外载荷相同的情况下，刚度增加能导致位移数值解偏小，但解得的应力应该在总体上保持不变，否则就无法满足平衡方程。

4. 减缩积分和完全积分单元

完全积分是指当单元具有规则形状时，所用的高斯积分点的数目足以对单元刚度矩阵中的多项式进行精确积分。

减缩积分是指对单元刚度阵进行积分时，所用的高斯积分点数低于精确积分的要求，一般是按照下式来确定高斯积分点数，

$$n=p-m+1 $$

公式 5

式中，n 为高斯积分点数，p 是插值函数中完全多项式的方次，m 是微分算子 L 中的导数阶次（对于用于弹塑性力学分析的实体单元，m=1）。

按照式（5）确定积分点数的减缩积分单元比完全积分单元在每个方向上减少了一个积分点。

对于大部分商用有限元软件，二维实体单元中的线性与二次四边形单元有减缩积分形式，积分点数分别为 1 * 1 和 2 *2。三维实体单元中的线性和二次六面体单元有减缩积分形式，积分点数分别为 1 *1* 1 和 2 *2*2。三角形单元、四面体单元和楔形单元一般没有减缩积分形式。

图 3 完全积分和减缩积分单元

因为降低了刚度矩阵的积分精度，减缩积分单元看似会影响结果的准确性。但实际计算表明，采用减缩积分往往可以取得比精确积分更好的精度，主要原因如下：

• 精确积分常常是由插值函数中非完全多项式的最高方次所要求，而决定有限元精度的，通常是完全多项式的方次。这些非完全的最高方次项往往不能提高精度，反而有负面影响。采用式（5）的减缩积分方案，积分精度恰好保证完全多项式方次要求，实质是用一种新的插值函数代替原插值函数，从而改善单元精度。
• 在最小位能原理基础上建立的位移有限元，计算模型具有较实际结构偏大的整体刚度。选取减缩积分方案将使有限元模型的刚度有所降低，因此可能有助于提高精度。
• 减缩积分方案对于泛函中包含罚函数的情况也常常是必须的，以保证罚函数矩阵的奇异性。例如基于相对自由度的梁单元和板壳单元，采用完全积分会出现剪切锁死问题，导致出现完全歪曲的结果，改用减缩积分方案就可以有效解决。
• 减缩积分因为积分点数比全积分少很多，所以大幅减少了计算量，提高了计算效率。

实际工程应用中，完全积分单元容易出现剪切锁死和体积锁死等问题，即使划分很细的网格，精度依然很差，所以一般不推荐使用。