自研有限元程序的减缩积分单元如何添加沙漏刚度（理论解释+数值实现）

今天分享的是减缩积分单元如何在原有单元刚度矩阵的基础上添加沙漏刚度，以防止虚假沙漏模态的出现。

主要包括以下内容：

• 历史背景
• 添加沙漏刚度（数值实现）
• 数值案例
• 精度讨论

历史背景

线性单元在有限元计算过程中每个单元的刚度矩阵形成之前需要遍历4个高斯积分点，前辈们为了尽可能缩短有限元计算时间，将单元每个方向上减少一个积分点（减缩积分）来减少计算成本。对于原先$2\times 2$分布的高斯积分点，减缩积分后每个单元只保留1个积分点，单元刚度计算时只需计算该积分点的几何矩阵$B$即可，大大减少了计算时间。

但是在减少时间成本的同时，相应的弊端也显现了出来。如下图所示，为1/4圆环有限元模型承受内壁压力荷载作用下的位移云图，采用用平面等参四节点减缩积分计算，缩放因子为0.001。

从上图可以看出单元在受到荷载作用后，出现严重的"软化"现象，有限元中常常称之为沙漏现象，显然是不符合我们对结果的预期的。

这时候应该怎么办呢？从物理现象角度来看，既然单元太软，那就让单元变硬点！具体如何实施呢？请读者耐心往下阅读学习，且看木木如何在原有单元刚度基础上添加刚度使得减轻单元的沙漏现象。

沙漏刚度

本小节的内容主要想带着读者了解沙漏刚度的公式是如何数值实现，重要是将结论写进代码中！等后面自己对概念相对熟悉一些后，可回过头来翻阅公式的演变过程，本节的主要参考资料为张雄老师的*《有限元法基础》P234*。

Q4（平面等参四节点）单元的位移模式可写为：

$$N_I^{Q4}(\xi ,\eta )=\frac{1}{4}\left({\Sigma }_I+{\Lambda }_{1I}\xi +{\Lambda }_{2I}\eta +{\Gamma }_I\xi \eta \right)$$

上式中的四个基向量分别代表刚体平移（${\Sigma }_I$）、拉压变形（${\Lambda }_{1I}$）、剪切变形（${\Lambda }_{2I}$）、沙漏变形（${\Gamma }_I$），如下图所示。若单元发生沙漏变形，则单元形状呈现梯形，单元中心应变为0。

相关的具体公式的推导因本人水平有限不能以大白话讲之，遂弃之。可参照推文：https://mp.weixin.qq.com/s/4WF2adiWb-ZAINNmIEsVcg或张雄老师的*《有限元法基础》P234*。

沙漏模态向量$\gamma$可表示为：
$$\gamma =\frac{1}{A}\left[\begin{array}{c}{x}_2{y}_{34}+x_3{y}_{42}+x_4{y}_{23}\\ x_3{y}_{14}+x_4{y}_{31}+x_1{y}_{43}\\ x_4{y}_{12}+x_1{y}_{24}+x_2{y}_{41}\\ x_1{y}_{32}+x_2{y}_{13}+x_3{y}_{21}\end{array}\right]$$

其中，$x$、$y$指的是单元节点坐标，$A$指的是单元面积。引入单元沙漏刚度矩阵：

$$K_{ab}^{\mathrm{HG}}=\kappa V_e{\gamma }_I{\gamma }_J$$
式中
$$\{\begin{array}{l}a=\left(I-1\right)n_{dof}+i\\ b=\left(J-1\right)n_{dof}+i,i=1:n_{dof}\end{array}$$
，$V_e$为单元体积，系数$\kappa$一般可取为：

$$\kappa =0.01G{N}_{I,i}{N}_{I,i}$$
其中$G$为剪切模量，$N_{I,i}$表示物理坐标系下形函数对坐标轴的偏导数，用户可将调为形参用于控制沙漏刚度大小。

以上的理论就是程序中所用到的公式，比如在求沙漏刚度时，程序可以编写循环嵌套结构，我用到了四重循环，其实是一种比较笨的遍历方法，不过可以初步实现该公式，用户可以依据自己能力进行代码优化。

for a  = 1 : elenode
    for i  = 1 : ndof
      for b  = 1 : elenode
        for j  = 1 : ndof
          row = ndof*(a-1)+i;
          col = ndof*(b-1)+j;
          hs = hs_vectors(a)*hs_vectors(b)*hourglass_k * Ae;
          k(col,row) = k(col,row) + hs;
        end
      end
    end
end

系数$\kappa$在程序可使用二重循环，可体现为：

for a  = 1 : elenode
    for i  = 1 : ndof
      DNNorm = DNNorm + XYderivatives(a,i)*XYderivatives(a,i);
    end
end
hourglass_k = 0.01*G*DNNorm;

数值案例

现考虑如下图所示的$\frac{1}{4}$圆环，下边界约束$y$方向自由度，左边界约束$x$方向自由度。内壁承受大小为1的压强，弹性模量取13，泊松比取0.3，单元厚度为1，将之作为平面应变问题看待，径向位移的解析解可表示为：

$$u_r=\frac{(1+\nu )a^2{b}^2}{E\left(b^2-a^2\right)}\left[\frac{p_a-p_b}{r}+(1-2\nu )\frac{p_a{a}^2-p_b{b}^2}{a^2{b}^2}r\right]$$

程序分别按照完全积分和减缩积分进行有限元计算，位移结果罗列于下图，

左图为不采用沙漏控制计算得到的位移变形云图（变形放大系数为0.01），右图在原有单元刚度矩阵的基础上添加沙漏刚度计算得到的位移变形云图（变形放大系数为1）。这两个图可以看出添加沙漏刚度控制后，极大程度上降低了虚假的沙漏模态，保证了有限元解的收敛性。

精度讨论

内壁上径向的位移解析解为0.1173，精度对比罗列于下表
。从下表可以看出，由减缩积分得到的径向位移比由完全积分得到的结果更接近解析解，精度更高！

加入沙漏刚度后的有限元程序也可以通过拼片试验，精度良好！

---

本次分享的理论解释和matlab代码均上传至《有限元基础编程百科全书》最新版PDF中，获取方式：自研有限元程序的减缩积分单元如何添加沙漏刚度（理论解释+数值实现）