有限元基础编程 | “选择性减缩积分”理论推导&数值实现

本次给大家分享的是：有限元数值编程中面对几乎不可压问题时，如何选择性减缩积分？

低阶单元完全积分在分析几乎不可压缩问题（泊松比接近0.5）时，易出现体积自锁，精度是非常差的！为适应此类情况，前辈们拓展了积分方案（减缩积分）以及各种数值技术来适应不可压问题，其中选择性减缩积分和B-BAR技术比较具有代表性。

本节主要分享选择性减缩积分技术，B-BAR的具体数值实现留作下一篇推文介绍，本次主要内容有：

• 何为体积自锁？
• 选择性减缩积分理论公式推导
• 数值案例
• 精度对比

何为体积自锁？

如下图所示的平面应变问题（${\epsilon }_z=0$），体积应变应为${\epsilon }_x+{\epsilon }_y$，对于不可压缩或者几乎不可压缩问题（如橡胶材料）来说，其泊松比无限接近于0.5，体积模量接近于无穷大，单元的体积在受力过程中几乎保持不变。

上图中的1、2、5号节点的自由度被固定，若要使单元1的体积保持不变，节点6不能沿着x方向移动，同时为了保持单元2的体积不变，节点6同样不能沿着y方向移动，因此就造成了节点6被固定，同样的道理，延伸至其他节点也是被固定，模型在受力过程中，网格各个节点均不能发生移动，产生了"闭锁"现象，有限元世界里常称之为"体积自锁"（volumetric
locking）或"网格自锁"（mesh locking）。

以上是从物理现象角度出发，解释何为体积自锁，其实背后有其数学物理方程来解释，本文档不打算记录该物理过程，感兴趣的读者可翻阅张雄老师的《有限元法基础》。

选择性减缩积分

公式推导

从这个名字来看，顾名思义就是要选择性的进行减缩积分计算，不要被名字所吓到，静下心一点一点学习这个概念，学习前辈们的智慧。

单元刚度矩阵$K$可以分解为偏刚度矩阵$K_{dev}$（deviatoric stiffness matrix）和体积刚度矩阵$K_{dil}$（volumetric/dilatational
stiffness matrix），类似于弹性力学中将应力分解为偏应力和静水压力。

$$K=K_{dev}+K_{dil}$$

对于几乎不可压缩问题，体积模量趋于无穷大，体积刚度矩阵$K_{dil}$远远大于偏刚度矩阵$K_{dev}$，若要避免自锁（体积自锁或剪切自锁）现象，必须使得体积刚度矩阵阵$K_{dil}$发生奇异时（有非零解），相关理论解释可翻阅《有限元法基础》P275。

在数值积分时，对于偏刚度矩阵$K_{dev}$可采用完全积分，对于体积刚度矩阵$K_{dil}$可采用减缩积分，保证其矩阵奇异。

对于体积自锁现象，用户也可使用全局减缩积分，但可能会出现虚假的沙漏模态，此时需要添加沙漏刚度避免沙漏现象，选择性减缩积分不必需要添加额外的沙漏刚度就可以保证解的收敛性。

接下来就详细展开，这两个刚度矩阵具体该如何求解。应力将其分解为体积应力（${\sigma }_{dil}$）和偏应力（${\sigma }_{dev}$）：
$$\{\begin{array}{l}\sigma ={\sigma }_{dil}+{\sigma }_{dev}\\ {\sigma }_{dil}=D_{dil}{\epsilon }_{dil}\\ {\sigma }_{dev}=D_{dev}{\epsilon }_{dev}\end{array}$$
其中，$D_{dil}$表示体积应变部分的弹性矩阵，$D_{}$