减缩积分单元、沙漏控制与自定义单元：与Abaqus C3D8R单元的精度对比之旅

为了提高计算效率，有限元中经常用到减缩积分单元，但是减缩积分单元会有沙漏现象。
因此减缩积分单元需要引入沙漏控制来保证计算精度，常见商业软件中（如ABAQUS）一阶单元的减缩积分通常采用一种均匀应变的方式，可参考：

Flanagan D P, Belytschko T. A uniform strain hexahedron and quadrilateral with orthogonal hourglass control[J]. International journal for numerical methods in engineering, 1981, 17(5): 679-706.

这种应变不是在一阶高斯点处直接计算得到，而是在体积上平均得到的。本文主要复现上述论文中的方法。

基本思想

均匀应变实际上是构造了一个线性的位移场，同时定义了和该线性位移场正交的沙漏形状向量，并用沙漏形状向量定义沙漏位移场，从而保证位移场的完全性。（个人理解，如有误欢迎指正）
大概的推导过程如下（详细推导请见论文）：

uiI=uiILIN+uiIHGu_{iI}=u_{iI}^{LIN}+u_{iI}^{HG}

uiLIN=uˉi+uˉi,j(xj−xˉj)u_{i}^{LIN}=\bar{u}_i+\bar{u}_{i,j}(x_j-\bar{x}_j)

其中

uˉi,j=1V∫Vui,jdV\bar{u}_{i,j}=\frac{1}{V}\int_{V}u_{i,j}dV

定义均匀Bˉ\bar{B}矩阵为

BˉiI=∫VϕI,idV\bar{B}_{iI}=\int_{V}\phi_{I,i}dV

因此可得

uˉi,j=1VuiIBˉjI\bar{u}_{i,j}=\frac{1}{V}u_{iI}\bar{B}_{jI}

xiIBˉjI=Vδijx_{iI}\bar{B}_{jI}=V\delta_{ij}

沙漏位移为

uiIHG=uiI−uiILINu_{iI}^{HG}=u_{iI}-u_{iI}^{LIN}

用沙漏基向量表示

uiIHG=18qiαΓαIu_{iI}^{HG}=\frac{1}{\sqrt{8}}q_{i \alpha}\Gamma_{\alpha I}

其中ΓαI\Gamma_{\alpha I}为沙漏基向量，qiαq_{i \alpha}沙漏模态位移

qiα=18uiIγαIq_{i \alpha}=\frac{1}{\sqrt{8}}u_{i I}\gamma_{\alpha I}

沙漏形状向量为

γαI=ΓαI−1VBˉiIXiJΓαJ\gamma_{\alpha I}=\Gamma_{\alpha I}-\frac{1}{V}\bar{B}_{iI}X_{iJ}\Gamma_{\alpha J}

单元内力

单元内力由俩部分构成，线性位移场（均匀应变）内力和沙漏力

fiIIN=fiIHG+fiIMEANf_{iI}^{IN}=f_{iI}^{HG}+f_{iI}^{MEAN}

fiIHG=18QiαγαI=18κλ+2μ3BˉiIBˉiIVqiαγαI=κλ+2μ24BˉiIBˉiIVuiJγαJγαIf_{iI}^{HG}=\frac{1}{\sqrt{8}}Q_{i \alpha}\gamma_{\alpha I}\\ =\frac{1}{\sqrt{8}} \kappa \frac{\lambda+2 \mu}{3} \frac{\bar{B}_{iI}\bar{B}_{iI}}{V}q_{i \alpha}\gamma_{\alpha I}\\ =\kappa \frac{\lambda+2 \mu}{24}\frac{\bar{B}_{iI}\bar{B}_{iI}}{V}u_{i J}\gamma_{\alpha J}\gamma_{\alpha I}

fiIMEAN=TˉijBˉjI=(λDˉkkδij+2μDˉij)BˉjI=λDˉkkBˉiI+2μDˉijBˉjIf_{iI}^{MEAN}=\bar{T}_{ij}\bar{B}_{jI}\\ =(\lambda\bar{D}_{kk}\delta_{ij}+2\mu\bar{D}_{ij})\bar{B}_{jI}\\ =\lambda\bar{D}_{kk}\bar{B}_{iI}+2\mu\bar{D}_{ij}\bar{B}_{jI}

Dˉij=12(uˉi,j+uˉj,i)\bar{D}_{ij}=\frac{1}{2}(\bar{u}_{i,j}+\bar{u}_{j,i})

刚度矩阵推导

和单元内力一样，刚度矩阵也由2部分构成,可由内力对位移求偏导得到

均匀应变刚度

ε=[uˉ1,1uˉ2,2uˉ3,3uˉ2,3+uˉ3,2uˉ1,3+uˉ3,1uˉ1,2+uˉ2,1]=1V[u1IBˉ1Iu2IBˉ2Iu3IBˉ3Iu2IBˉ3I+u3IBˉ2Iu1IBˉ3I+u3IBˉ1Iu1IBˉ2I+u2IBˉ1I]=1V[Bˉ1I000Bˉ2I000Bˉ3I0Bˉ3IBˉ2IBˉ3I0Bˉ1IBˉ2IBˉ1I0][u1Iu2Iu3I]=BU\varepsilon = \begin{bmatrix} \bar{u}_{1,1}\\ \bar{u}_{2,2}\\ \bar{u}_{3,3}\\ \bar{u}_{2,3}+\bar{u}_{3,2}\\ \bar{u}_{1,3}+\bar{u}_{3,1}\\ \bar{u}_{1,2}+\bar{u}_{2,1}\\ \end{bmatrix}=\frac{1}{V} \begin{bmatrix} u_{1I}\bar{B}_{1I}\\ u_{2I}\bar{B}_{2I}\\ u_{3I}\bar{B}_{3I}\\ u_{2I}\bar{B}_{3I}+u_{3I}\bar{B}_{2I}\\ u_{1I}\bar{B}_{3I}+u_{3I}\bar{B}_{1I}\\ u_{1I}\bar{B}_{2I}+u_{2I}\bar{B}_{1I}\\ \end{bmatrix}\\ =\frac{1}{V} \begin{bmatrix} \bar{B}_{1I}&0&0\\ 0&\bar{B}_{2I}&0\\ 0&0&\bar{B}_{3I}\\ 0&\bar{B}_{3I}&\bar{B}_{2I}\\ \bar{B}_{3I}&0&\bar{B}_{1I}\\ \bar{B}_{2I}&\bar{B}_{1I}&0\\ \end{bmatrix} \begin{bmatrix} u_{1I}\\ u_{2I}\\ u_{3I}\\ \end{bmatrix}=BU

K=BTDBVK=B^TDBV

沙漏刚度

fiIHG=κλ+2μ24BˉiIBˉiIVuiJγαJγαIf_{iI}^{HG}=\kappa \frac{\lambda+2 \mu}{24}\frac{\bar{B}_{iI}\bar{B}_{iI}}{V}u_{i J}\gamma_{\alpha J}\gamma_{\alpha I}

(KIJHG)ij=∂fiIHG∂ujJ(K^{HG}_{IJ})_{ij}=\frac{\partial{f_{iI}^{HG}}}{\partial{u_{jJ}}}

当i≠ji \neq j时，

(KIJHG)ij=0(K^{HG}_{IJ})_{ij}=0

当i=ji = j时，

(KIJHG)ij=κλ+2μ24BˉiIBˉiIVγαJγαI(K^{HG}_{IJ})_{ij}=\kappa \frac{\lambda+2 \mu}{24}\frac{\bar{B}_{iI}\bar{B}_{iI}}{V}\gamma_{\alpha J}\gamma_{\alpha I}

均匀B矩阵的计算

在上述公式中可以发现，要想得到刚度阵，只要得到沙漏形状向量γ\gamma和均匀Bˉ\bar{B}矩阵即可，而沙漏形状向量γ\gamma也是由均匀Bˉ\bar{B}矩阵计算得到的，因此关键的一步是如何得到均匀Bˉ\bar{B}矩阵。
论文中（22）和（23）式给出了均匀Bˉ\bar{B}矩阵的计算方法，我们可以通过符号计算软件直接求解出均匀B矩阵的表达式。计算均匀B矩阵的 mathmatica 的代码如下：

mlambda={{-1,1,1,-1,-1,1,1,-1},
          {-1,-1,1,1,-1,-1,1,1},
          {-1,-1,-1,-1,1,1,1,1}};         
MatrixForm[%];
mgamma={{1,1,-1,-1,-1,-1,1,1},
         {1,-1,-1,1,-1,1,1,-1},
         {1,-1,1,-1,1,-1,1,-1},
         {-1,1,-1,1,1,-1,1,-1}};
MatrixForm[%];
e=Table[100,{i,3},{j,3},{k,3}];
For[i=1, i <= 3, ++i,
  For[j=1, j <= 3, ++j,
    For[k=1, k <= 3, ++k,
        Which[
            ((i==1 && j ==2 && k ==3)||(i==2 && j ==3 && k ==1)||(i==3 && j ==1 && k ==2)),e[[i,j,k]]=1,
            ((i==3 && j ==2 && k ==1)||(i==2 && j ==1 && k ==3)||(i==1 && j ==3 && k ==2)),e[[i,j,k]]=-1,
            True,e[[i,j,k]]=0
            ]     
    ]
  ]
]
MatrixForm[e];
(*(23)*)
Crst=Table[
    Sum[
        Rationalize[1/192]*e[[i,j,k]]*(3*mlambda[[i,r]]*mlambda[[j,s]]*mlambda[[k,t]]+
                                        mlambda[[i,r]]*mgamma[[k,s]]*mgamma[[j,t]]+
                                        mgamma[[k,r]]*mlambda[[j,s]]*mgamma[[i,t]]+
                                        mgamma[[j,r]]*mgamma[[i,s]]*mlambda[[k,t]]),
    {i,3},{j,3},{k,3}],
    {r,8},{s,8},{t,8}];
MatrixForm[Crst];

xddd={Subscript[X,1],Subscript[X,2],Subscript[X,3],Subscript[X,4],Subscript[X,5],Subscript[X,6],Subscript[X,7],Subscript[X,8]};
yddd={Subscript[Y,1],Subscript[Y,2],Subscript[Y,3],Subscript[Y,4],Subscript[Y,5],Subscript[Y,6],Subscript[Y,7],Subscript[Y,8]};
zddd={Subscript[Z,1],Subscript[Z,2],Subscript[Z,3],Subscript[Z,4],Subscript[Z,5],Subscript[Z,6],Subscript[Z,7],Subscript[Z,8]};
(*(22)*)
bone=Table[Collect[Sum[12*yddd[[j]]*zddd[[k]]*Crst[[i,j,k]],{j,8},{k,8}],{Subscript[Y,1],Subscript[Y,2],Subscript[Y,3],Subscript[Y,4],Subscript[Y,5],Subscript[Y,6],Subscript[Y,7],Subscript[Y,8]}],{i,8}];
btwo=Table[Collect[Sum[12*zddd[[j]]*xddd[[k]]*Crst[[i,j,k]],{j,8},{k,8}],{Subscript[Z,1],Subscript[Z,2],Subscript[Z,3],Subscript[Z,4],Subscript[Z,5],Subscript[Z,6],Subscript[Z,7],Subscript[Z,8]}],{i,8}];
bthree=Table[Collect[Sum[12*xddd[[j]]*yddd[[k]]*Crst[[i,j,k]],{j,8},{k,8}],{Subscript[X,1],Subscript[X,2],Subscript[X,3],Subscript[X,4],Subscript[X,5],Subscript[X,6],Subscript[X,7],Subscript[X,8]}],{i,8}];

这里给出计算得到Bˉ1\bar{B}_{1}的表达式，其中Bˉ11\bar{B}_{11}的结果和论文给出的结果（式（79））一致。

Bˉ1I=112[Y3(Z2−Z4)+Y8(Z4−Z5)+Y6(Z5−Z2)+Y2(−Z3−Z4+Z5+Z6)+Y4(Z2+Z3−Z5−Z8)+Y5(−Z2+Z4−Z6+Z8)Y4(Z3−Z1)+Y5(Z1−Z6)+Y1(Z3+Z4−Z5−Z6)+Y7(Z6−Z3)+Y6(Z1−Z3+Z5−Z7)+Y3(−Z1−Z4+Z6+Z7)Y1(Z4−Z2)+Y6(Z2−Z7)+Y2(Z1+Z4−Z6−Z7)+Y8(Z7−Z4)+Y7(Z2−Z4+Z6−Z8)+Y4(−Z1−Z2+Z7+Z8)Y2(Z1−Z3)+Y8(−Z1+Z3−Z5+Z7)+Y7(Z3−Z8)+Y3(Z1+Z2−Z7−Z8)+Y5(Z8−Z1)+Y1(−Z2−Z3+Z5+Z8)Y2(Z6−Z1)+Y8(Z1+Z4−Z6−Z7)+Y4(Z1−Z8)+Y1(Z2−Z4+Z6−Z8)+Y7(Z8−Z6)+Y6(−Z1−Z2+Z7+Z8)Y1(Z2−Z5)+Y8(Z5−Z7)+Y3(Z7−Z2)+Y2(−Z1+Z3−Z5+Z7)+Y5(Z1+Z2−Z7−Z8)+Y7(−Z2−Z3+Z5+Z8)Y2(Z3−Z6)+Y8(−Z3−Z4+Z5+Z6)+Y6(Z2+Z3−Z5−Z8)+Y5(Z6−Z8)+Y4(Z8−Z3)+Y3(−Z2+Z4−Z6+Z8)Y1(Z5−Z4)+Y7(Z3+Z4−Z5−Z6)+Y3(Z4−Z7)+Y4(Z1−Z3+Z5−Z7)+Y6(Z7−Z5)+Y5(−Z1−Z4+Z6+Z7)]\bar{B}_{1I}=\frac{1}{12} \begin{bmatrix} Y_3\left(Z_2-Z_4\right)+Y_8\left(Z_4-Z_5\right)+Y_6\left(Z_5-Z_2\right)+Y_2\left(-Z_3-Z_4+Z_5+Z_6\right)+Y_4\left(Z_2+Z_3-Z_5-Z_8\right)+Y_5\left(-Z_2+Z_4-Z_6+Z_8\right)\\ Y_4\left(Z_3-Z_1\right)+Y_5\left(Z_1-Z_6\right)+Y_1\left(Z_3+Z_4-Z_5-Z_6\right)+Y_7\left(Z_6-Z_3\right)+Y_6\left(Z_1-Z_3+Z_5-Z_7\right)+Y_3\left(-Z_1-Z_4+Z_6+Z_7\right)\\ Y_1\left(Z_4-Z_2\right)+Y_6\left(Z_2-Z_7\right)+Y_2\left(Z_1+Z_4-Z_6-Z_7\right)+Y_8\left(Z_7-Z_4\right)+Y_7\left(Z_2-Z_4+Z_6-Z_8\right)+Y_4\left(-Z_1-Z_2+Z_7+Z_8\right)\\ Y_2\left(Z_1-Z_3\right)+Y_8\left(-Z_1+Z_3-Z_5+Z_7\right)+Y_7\left(Z_3-Z_8\right)+Y_3\left(Z_1+Z_2-Z_7-Z_8\right)+Y_5\left(Z_8-Z_1\right)+Y_1\left(-Z_2-Z_3+Z_5+Z_8\right)\\ Y_2\left(Z_6-Z_1\right)+Y_8\left(Z_1+Z_4-Z_6-Z_7\right)+Y_4\left(Z_1-Z_8\right)+Y_1\left(Z_2-Z_4+Z_6-Z_8\right)+Y_7\left(Z_8-Z_6\right)+Y_6\left(-Z_1-Z_2+Z_7+Z_8\right)\\ Y_1\left(Z_2-Z_5\right)+Y_8\left(Z_5-Z_7\right)+Y_3\left(Z_7-Z_2\right)+Y_2\left(-Z_1+Z_3-Z_5+Z_7\right)+Y_5\left(Z_1+Z_2-Z_7-Z_8\right)+Y_7\left(-Z_2-Z_3+Z_5+Z_8\right)\\ Y_2\left(Z_3-Z_6\right)+Y_8\left(-Z_3-Z_4+Z_5+Z_6\right)+Y_6\left(Z_2+Z_3-Z_5-Z_8\right)+Y_5\left(Z_6-Z_8\right)+Y_4\left(Z_8-Z_3\right)+Y_3\left(-Z_2+Z_4-Z_6+Z_8\right)\\ Y_1\left(Z_5-Z_4\right)+Y_7\left(Z_3+Z_4-Z_5-Z_6\right)+Y_3\left(Z_4-Z_7\right)+Y_4\left(Z_1-Z_3+Z_5-Z_7\right)+Y_6\left(Z_7-Z_5\right)+Y_5\left(-Z_1-Z_4+Z_6+Z_7\right) \end{bmatrix}

验证算例

C3D8在弯曲载荷下精度较差，主要是因为在弯曲载荷下会产生parasitic shear stresses，因此通常在弯曲方向画分多层网格，或者使用C3D8I单元来提高精度。同时弯曲变形也是沙漏模态之一，因此特选如下模型来说明沙漏控制的有效性。这里我们暂且把C3D8I的结果作为对比的参考解。

沙漏模态

模型尺寸及边界

模型尺寸及边界

模型厚度为1，一端固定，一端对称约束，p为侧面压力载荷。

结果对比

C3D8I单元结果

C3D8I单元位移值

C3D8R单元结果
单元设置如下：

单元设置
这里displacement hourglass factors可以认为是上述刚度矩阵里的κ\kappa。（其实κ\kappa和displacement hourglass factors相差了个系数，同时测试时发现ABAQUS的默认系数比较小)

结果对比如下：

C3D8R displacement hourglass factors为1.0时，ABAQUS结果和自编程序结果对比：

Abaqus结果

自编程序结果

C3D8R displacement hourglass factors为50.0时，ABAQUS结果和自编程序结果对比：

ABAQUS结果

自编程序结果

总结

ABAQUS C3D8IABAQUS C3D8R factor=1.0自编 C3D8R factor=1.0ABAQUS C3D8R factor=50.0自编 C3D8R factor=50.0
最大位移0.092368.2088.2080.16450.1645
  1. 可见在displacement hourglass factors在较小时，沙漏现象严重，位移严重失真，增加displacement hourglass factors可以大幅提高精度，说明沙漏控制的有效性。
  2. 同时也说明按照论文方法自编单元和ABAQUS C3D8R带displacement hourglass控制选项单元基本完全一致。

以上内容是个人对于减缩积分和沙漏控制的一点浅显理解，如有疑问，欢迎在下方留言区参与讨论，共勉。

有限元
#有限元编程
减缩积分单元、沙漏控制与自定义单元：与Abaqus C3D8R单元的精度对比之旅
http://yimumu.top/2023/08/18/C3D8R沙漏控制/
作者
易公子
发布于
2023年8月18日
更新于
2023年9月1日
许可协议