减缩积分单元、沙漏控制与自定义单元：与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.

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

基本思想

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

其中

定义均匀    矩阵为

因此可得

沙漏位移为

用沙漏基向量表示

其中    为沙漏基向量，    沙漏模态位移

沙漏形状向量为

单元内力

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

刚度矩阵推导

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

均匀应变刚度

沙漏刚度

当    时，

当    时，

均匀B矩阵

在上述公式中可以发现，要想得到刚度阵，只要得到沙漏形状向量    和均匀    矩阵即可，而沙漏形状向量    也是由均匀    矩阵计算得到的，因此关键的一步是如何得到均匀    矩阵。论文中（22）和（23）式给出了均匀    矩阵的计算方法，我们可以通过符号计算软件直接求解出均匀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}];

这里给出计算得到    的表达式，其中    的结果和论文给出的结果（式（79））一致。

验证算例

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

 

模型尺寸及边界

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

结果对比

C3D8I单元结果

 

C3D8R单元结果

单元设置如下：

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

结果对比如下：

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

 

 

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

 

 

总结

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

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