C3D8R单元沙漏刚度如何添加？理论&有限元数值实现
在学习本节之前建议读者先移步节自研有限元程序的减缩积分单元如何添加沙漏刚度（理论解释+数值实现）进行平面单元的沙漏刚度学习，三维线性等参单元（C3D8）的沙漏刚度理论相对来讲比较难，不过本节还是从初学者的角度出发，逐渐揭开C3D8R单元神秘的面纱，力争让看到此文档的读者都可以编制出带有沙漏控制的C3D8R单元程序。
主要的参考文献是Ted《A uniform strain hexahedron and quadrilateral with orthogonal hourglass control》，以下简称“文献”和之前的公 众号推文：减缩积分单元、沙漏控制与自定义单元：与Abaqus C3D8R单元的精度对比之旅。
由于木木的理论水平有限，仅分享我在学习时的个人感悟，对程序编制有帮助的部分记录于此，若读者想要全面了解沙漏刚度的理论框架，可参照“文献”。
为了提高计算效率，有限元中经常用到减缩积分单元，但是减缩积分单元会有沙漏现象。因此减缩积分单元需要引入沙漏控制来保证计算精度，常见商业软件中（如ABAQUS）一阶单元的减缩积分通常采用一种均匀应变的方式。这种应变不是在一阶高斯点处直接计算得到，而是在体积上平均得到的。本节主要复现该数值方法。
基本思想
均匀应变实际上是构造了一个线性的位移场，同时定义了和该线性位移场正交的沙漏形状向量，并用沙漏形状向量定义沙漏位移场，从而保证位移场的完全性。
C3D8单元的位移模式加入沙漏基向量后可表示为：
C3D8单元位移模式各项系数见下表：
| Node | ξ | η | ζ | ΣI | Λ₁I | Λ₂I | Λ₃I | Γ₁I | Γ₂I | Γ₃I | Γ₄I |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | -1/2 | -1/2 | -1/2 | 1 | -1 | -1 | -1 | 1 | 1 | 1 | -1 |
| 2 | 1/2 | -1/2 | -1/2 | 1 | 1 | -1 | -1 | 1 | -1 | -1 | -1 |
| 3 | 1/2 | 1/2 | -1/2 | 1 | 1 | 1 | -1 | -1 | -1 | 1 | -1 |
| 4 | -1/2 | 1/2 | -1/2 | 1 | -1 | 1 | -1 | -1 | 1 | -1 | -1 |
| 5 | -1/2 | -1/2 | 1/2 | 1 | -1 | -1 | 1 | 1 | 1 | 1 | 1 |
| 6 | 1/2 | -1/2 | 1/2 | 1 | 1 | -1 | 1 | 1 | -1 | -1 | 1 |
| 7 | 1/2 | 1/2 | 1/2 | 1 | 1 | 1 | 1 | -1 | -1 | 1 | 1 |
| 8 | -1/2 | 1/2 | 1/2 | 1 | -1 | 1 | 1 | -1 | 1 | -1 | 1 |
六面体单元由上式中的八个基向量表示，其中， 表示刚体平移， 表示体积基向量（拉压变形、剪切变形、旋转）， 导致沙漏变形，称为沙漏基向量。接下来直接沿用“文献”中的插图，较为形象地解释各个基向量所代表的含义。
节点位移 被分解为两部分：[线性部分]{style="color: blue"} ，描述刚体运动和均匀应变场和[沙漏部分]{style="color: blue"} ，包含与线性模式正交的高阶零能量模式（hourglass模式），即：
线性部分表示为：
其中， 表示均匀应变， 表示节点相对于单元质心位置的偏移量。
均匀应变可通过对整个单元体积 进行积分得到，体现了均匀应变模式的特点：
引入均匀应变矩阵 矩阵
因此可得
其中 是Kroneckerdelta，用于表示节点间的正交性，确保了 矩阵与单元体积 的一致性
沙漏位移 是总位移 与线性部分 的差值：
沙漏位移可以用沙漏基向量表示：
为沙漏基向量，定义了单元内hourglass模式的形状， 为沙漏模态位移：
沙漏模态向量在表达上相对于CPS4单元较为复杂，C3D8单元的沙漏模态向量 可表示为：
其中， 为均匀应变矩阵， 为节点坐标， 为单元体积。
当六面体单元是平行六面体单元时， ，原文：Theyare identical if and only if the hexahedron is a parallelepiped.
刚度矩阵推导
单元刚度矩阵由两部分组成：均匀应变刚度矩阵和沙漏刚度矩阵。
沙漏引起的单元内力表示为：
沙漏刚度可以进一步表示为：
当 时，
当 时，
在上述公式中可以发现，要想得到单元刚度矩阵，只要得到沙漏形状向量 和均匀应变矩阵 即可，而沙漏形状向量 也是由均匀应变矩阵 计算得到的，因此关键的一步是如何得到均匀应变矩阵。
论文中（22）和（23）式给出了均匀应变矩阵 的计算方法，我们可以通过符号计算软件直接求解出均匀应变 矩阵的表达式。计算均匀应变 矩阵的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}
];
通过上面的mma代码可进一步推导出均匀应变 矩阵，看似很复杂，实则也不是很简单。
用户在matlab中进行编程时，可以使用符号语言进行公式推导，但是如果把符号语言写进程序里会降低该单元的刚度组装效率，建议直接将均匀应变 矩阵显式写出，由于将 完全展开需要占据较大篇幅，这里仅给出 的表达式：
将之编制进matlab中，由于编程能力有限，编程风格仅供参考。更多计算细节详看MFEA单元库中的C3D8R单元计算程序。
C = elemNodeCoordinate; % 单元节点坐标
bbar = zeros(3, 8); % 初始化 bbar 矩阵大小为 3x8
bbar(1,1) = (1/12) * (
C(3,2) * (C(2,3) - C(4,3)) +
C(8,2) * (C(4,3) - C(5,3)) +
C(6,2) * (-C(2,3) + C(5,3)) + ...
C(2,2) * (-C(3,3) - C(4,3) + C(5,3) + C(6,3)) +
C(4,2) * (C(2,3) + C(3,3) - C(5,3) - C(8,3)) + ...
C(5,2) * (-C(2,3) + C(4,3) - C(6,3) + C(8,3))); % 第1列
bbar(1,2) = (1/12) * (
C(4,2) * (-C(1,3) + C(3,3)) +
C(5,2) * (C(1,3) - C(6,3)) +
C(1,2) * (C(3,3) + C(4,3) - C(5,3) - C(6,3)) + ...
C(7,2) * (-C(3,3) + C(6,3)) +
C(6,2) * (C(1,3) - C(3,3) + C(5,3) - C(7,3)) + ...
C(3,2) * (-C(1,3) - C(4,3) + C(6,3) + C(7,3))); % 第2列
bbar(1,3) = (1/12) * (
C(1,2) * (-C(2,3) + C(4,3)) +
C(6,2) * (C(2,3) - C(7,3)) +
C(2,2) * (C(1,3) + C(4,3) - C(6,3) - C(7,3)) + ...
C(8,2) * (-C(4,3) + C(7,3)) +
C(7,2) * (C(2,3) - C(4,3) + C(6,3) - C(8,3)) + ...
C(4,2) * (-C(1,3) - C(2,3) + C(7,3) + C(8,3))); % 第3列
bbar(1,4) = (1/12) * (
C(2,2) * (C(1,3) - C(3,3)) +
C(8,2) * (-C(1,3) + C(3,3) - C(5,3) + C(7,3)) +
C(7,2) * (C(3,3) - C(8,3)) + ...
C(3,2) * (C(1,3) + C(2,3) - C(7,3) - C(8,3)) +
C(5,2) * (-C(1,3) + C(8,3)) + ...
C(1,2) * (-C(2,3) - C(3,3) + C(5,3) + C(8,3))); % 第4列
bbar(1,5) = (1/12) * (
C(2,2) * (-C(1,3) + C(6,3)) +
C(8,2) * (C(1,3) + C(4,3) - C(6,3) - C(7,3)) +
C(4,2) * (C(1,3) - C(8,3)) + ...
C(1,2) * (C(2,3) - C(4,3) + C(6,3) - C(8,3)) +
C(7,2) * (-C(6,3) + C(8,3)) + ...
C(6,2) * (-C(1,3) - C(2,3) + C(7,3) + C(8,3))); % 第5列
bbar(1,6) = (1/12) * (
C(1,2) * (C(2,3) - C(5,3)) +
C(8,2) * (C(5,3) - C(7,3)) +
C(3,2) * (-C(2,3) + C(7,3)) + ...
C(2,2) * (-C(1,3) + C(3,3) - C(5,3) + C(7,3)) +
C(5,2) * (C(1,3) + C(2,3) - C(7,3) - C(8,3)) + ...
C(7,2) * (-C(2,3) - C(3,3) + C(5,3) + C(8,3))); % 第6列
bbar(1,7) = (1/12) * (
C(2,2) * (C(3,3) - C(6,3)) +
C(8,2) * (-C(3,3) - C(4,3) + C(5,3) + C(6,3)) +
C(6,2) * (C(2,3) + C(3,3) - C(5,3) - C(8,3)) + ...
C(5,2) * (C(6,3) - C(8,3)) +
C(4,2) * (-C(3,3) + C(8,3)) + ...
C(3,2) * (-C(2,3) + C(4,3) - C(6,3) + C(8,3))); % 第7列
bbar(1,8) = (1/12) * (
C(1,2) * (-C(4,3) + C(5,3)) +
C(7,2) * (C(3,3) + C(4,3) - C(5,3) - C(6,3)) +
C(3,2) * (C(4,3) - C(7,3)) + ...
C(4,2) * (C(1,3) - C(3,3) + C(5,3) - C(7,3)) +
C(6,2) * (-C(5,3) + C(7,3)) + ...
C(5,2) * (-C(1,3) - C(4,3) + C(6,3) + C(7,3))); % 第8列
数值算例
C3D8单元在弯曲载荷下精度较差，主要是因为在弯曲载荷下会产生parasitic shear stresses，因此通常在弯曲方向画分多层网格，或者使用C3D8I单元来提高精度。同时弯曲变形也是沙漏模态之一，因此特选如下模型来说明沙漏控制的有效性。这里我们暂且把C3D8I的结果作为对比的参考解。
有限元模型见下图，厚度设置为1，左端固定，上表面施加向下的表面压强荷载，大小为10。试用Abaqus内置的C3D8R（沙漏默认设置）、Abaqus内置的C3D8R（沙漏增强设置）、Abaqus内置的C3D8I单元以及MFEA的C3D8R单元分别进行有限元分析，将位移场进行对比分析，一方面用于研究单元性质另一方面验证自研程序的合理性。
位移云图以及各单元类型最大位移对比如下所示，
| C3D8I | C3D8R-默认沙漏控制 | C3D8R-增强沙漏控制 | MFEA-C3D8R | |
|---|---|---|---|---|
| 最大位移 | 3.752e-1 | 3.574e1 | 3.753e-1 | 3.756e-1 |
从上面结果可以看出，Abaqus内置的C3D8R默认沙漏控制似乎并不是很有效，比C3D8I单元的位移场大了两个数量级，经过增强设置之后的位移精度与C3D8I单元几乎一样。
注：自研程序的C3D8R单元位移精度与C3D8I精度相当，程序中不是一次计算完成的，需要根据实际情况来调节比例系数，来逐渐达到C3D8I的精度。该单元的程序有待改进，目前的水平就只能设计于此了，望读者朋友见谅！
