C3D8 中的选择性缩减积分：从体积自锁到 B-bar 方法

简介

八节点六面体单元 C3D8 是有限元中最常见的三维实体单元之一。但对于近似不可压缩材料，或者问题逐渐进入不可压缩状态时，标准全积分 C3D8 很容易出现一个经典问题：体积自锁。在《有限元中的 Locking：为什么单元会变得异常刚硬》已有介绍。

解决这一问题有多种方法，例如混合 (u-p) 单元、选择性缩减积分（Selective Reduced Integration, SRI）以及 B-bar（\(\bar B\)）方法等。这一节主要讨论 C3D8 中常见的一类处理：选择性缩减积分。

Why SRI?

如果要满足材料不可压，标准全积分 C3D8 单元需要在每个积分点上都满足局部不可压缩条件。但是一个低阶 C3D8 位移场的自由度有限，这会导致单元完全无法动弹，于是表现得异常硬，这就是 volumetric locking。

这时候，我们自然而然就会想到：如何不把单元约束得太死呢？

现在我们不动 deviatoric 部分，只处理 volumetric 部分。原来体积刚度为：

\[K_{\mathrm{vol}} \approx \sum_{g=1}^{8} B_v(\xi_g)^T K B_v(\xi_g)J_gw_g \]

其中 \(K\) 是体积模量（Bulk Modulus）。现在对体积项使用 \(1\times1\times1\) 积分：

\[\boxed{ K_{\mathrm{vol}} \approx B_v(0)^T K B_v(0)V } \]

也就是只在单元中心计算体积项。

于是当 \(K\rightarrow\infty\)时，原来要求的是：\(\varepsilon_v(\xi_1) = \varepsilon_v(\xi_2) = \cdots = \varepsilon_v(\xi_8) = 0\)，现在主要约束变成：\(\varepsilon_v(0)=0\)。这样，对离散位移场施加的体积约束显著减少了，原本被错误抑制的等体积 deformation modes 就可以重新出现，这就是 SRI 能解决 locking 的本质。

为什么不能把整个 C3D8 都改成 1 点积分？

这马上引出另一个问题。既然减缩积分这么好，那干脆整个 \(K = K_{\mathrm{dev}}+K_{\mathrm{vol}}\)，全部使用 \(1\times1\times1\) 积分不就行了吗？答案是不行，因为这样会产生新的问题：hourglass modes。某些非零节点位移模式可能在单元中心恰好满足 \(B(0)u=0\)，于是积分点认\(\varepsilon(0)=0\)，从而\(U_e=0\)，但节点实际上已经发生了非刚体变形。这就是零能 hourglass mode。

这里就不展开说了，先埋个坑，我会在后面的文章里介绍如何避免 hourglass modes。

几何线性情况下的推导

从应变分解到 \(\bar B\)

如果只考虑几何线性问题，那么整个思路其实相当清楚。标准有限元使用 \(\varepsilon=\boldsymbol B\boldsymbol u\)，现在重新定义 \(\boldsymbol B\)：

• 偏应变部分仍然使用原来的局部 \(\boldsymbol B\)；
• 体积应变部分不再完全使用积分点处的局部值，而使用单元平均体积应变。

因此构造新的应变—位移矩阵

\[\boldsymbol B = \boldsymbol B_\mathrm{dev} + \bar{\boldsymbol B}_\mathrm{vol} \]

其中 \(\bar{\boldsymbol B}_{\mathrm{vol}}\) 来自单元体积平均，使用平均形函数梯度 \(\frac{1}{V}\int_V N_{I,i},dV\) 来表示。

这就是整个方法最重要的一步，也叫做B-Bar 方法。随后仍然按照普通有限元流程计算\(\int_V \bar{\boldsymbol B}^{T} \boldsymbol D \bar{\boldsymbol B} \, dV\)。

因此从程序结构上看，它与普通实体单元并没有本质区别：

\[\boxed{ \text{shape function} \rightarrow \nabla N \rightarrow B \rightarrow \varepsilon \rightarrow \sigma \rightarrow K } \]

真正发生变化的是：

\[\boxed{ B \rightarrow \bar B } \]

或者更准确地说：

\[\boxed{ B_{\mathrm{vol}} \rightarrow \bar B_{\mathrm{vol}}} \]

平均梯度的计算

可以利用 Flanagan 与 Belytschko [1] 给出的 uniform strain formulation，直接根据节点坐标计算相应的几何量。例如其中一个分量可以写成

\[\frac{1}{12}\Big[ y_2((z_6-z_3)-(z_4-z_5)) +y_3(z_2-z_4)\ +y_4((z_3-z_8)-(z_5-z_2))\ +y_5((z_8-z_6)-(z_2-z_4))\ +y_6(z_5-z_2) +y_8(z_4-z_5) \Big]. \]

其他分量可以通过节点以及坐标轴的循环置换得到。不过从算法理解的角度，并不需要首先记住这些展开式。真正重要的是它代表的积分量：

\[B_{iI} \sim \int_V N_{I,i},dV. \]

显式坐标公式更多是这一积分的高效实现方式。

如何得到单元体积？

得到上述几何量之后，可以利用恒等关系\(x_{iI} B_{jI} = V \delta_{ij}\)，例如取 (i=j=1)，得到

\[\boxed{ V=x_{1I}B_{1I} } \]

这里采用 Einstein 求和约定，即节点指标 (I) 自动求和。因此\(\sum_{I=1}^{8}x_I B_{1I}\)。

类似地，理论上也可以从另外两个方向得到相同的单元体积。这一关系非常有用，因为它把节点坐标、平均形函数梯度和单元体积联系了起来。

几何线性问题很简单，但是几何非线性下的推导会复杂得多，感兴趣的读者可以试着推导看看😏。

总结

最后总结一下，C3D8 的体积自锁，本质上来自低阶位移离散对近不可压缩条件施加了过强的局部约束。选择性缩减积分以及\(\bar B\)一类方法的共同思想，是降低这种不合理的体积约束。

References

[1] Flanagan, D. P., and T. Belytschko. ‘A Uniform Strain Hexahedron and Quadrilateral with Orthogonal Hourglass Control’. International Journal for Numerical Methods in Engineering 17, no. 5 (May 1981): 679–706. https://doi.org/10.1002/nme.1620170504.