高性能单元公式¶

三维实体单元的公式中所示的基于位移的标准公式，在应用于近似不可压缩材料或以弯曲为主的薄壁结构时，会表现出称为锁定（体积锁定、剪切锁定）的不合理过刚现象。为避免这一问题，FrontISTR 提供 B-bar 法和 F-bar 法，它们仅替换 B 矩阵或变形梯度的体积分量；还提供引入内部自由度的非协调单元、将压力作为独立未知场的 u-p 混合单元，以及面向板壳和梁结构的 MITC 壳单元与梁单元。

本章按单元类型整理这些高性能单元和结构单元的公式。

B-bar 法¶

将 8 节点线性六面体单元用于准不可压缩材料时，一个单元内的应变会与体积恒定约束产生矛盾，从而出现称为体积锁定的过刚现象。B-bar 法将 B 矩阵中对体积膨胀有贡献的分量替换为在单元中心处评估的值，从而缓解过约束 [Hughes1980]。

设在单元中心 \(\boldsymbol{r} = \boldsymbol{0}\) 处由形状函数的空间导数计算得到的 B 矩阵为 \(\bar{\boldsymbol{B}}\)，在积分点 \(\boldsymbol{r}\) 处计算得到的常规 B 矩阵为 \(\boldsymbol{B}(\boldsymbol{r})\)。在针对节点 \(\alpha\) 的自由度 \(i\) 的位移－应变关系中，对体积应变分量 \((\varepsilon_{11}, \varepsilon_{22}, \varepsilon_{33})\) 加上

\[ \Delta B_{i\alpha} = \tfrac{1}{3}\bigl(\bar{B}_{i\alpha}(\boldsymbol{0}) - B_{i\alpha}(\boldsymbol{r})\bigr) \]

而对剪切分量 \((\varepsilon_{12}, \varepsilon_{23}, \varepsilon_{31})\) 使用常规的 \(\boldsymbol{B}\)。然后使用所得的 B-bar 矩阵组装单元刚度和内力向量。

FrontISTR 仅针对 8 节点线性六面体（单元 ID 361，单元编号体系）提供该公式，并且可应用于小变形、Total Lagrange 和 Updated Lagrange。

F-bar 法¶

在有限变形下，体积变化会通过变形梯度 \(\boldsymbol{F}\) 以非线性方式起作用，因此 F-bar 法 [deSouzaNeto1996] 是在变形梯度层面实施与 B-bar 法等效的体积锁定对策。

设在单元中心 \(\boldsymbol{r} = \boldsymbol{0}\) 处评估的变形梯度体积比为 \(J_0 = \det \boldsymbol{F}(\boldsymbol{0})\)，积分点处的变形梯度体积比为 \(J = \det \boldsymbol{F}(\boldsymbol{r})\)，将积分点处的变形梯度替换为

\[ \bar{\boldsymbol{F}} = \left(\frac{J_0}{J}\right)^{1/3} \boldsymbol{F} \]

这样可得到 \(\det \bar{\boldsymbol{F}} = J_0\)，使整个单元内的体积比统一为单元中心值。应力评估和应变－位移矩阵的构建使用替换后的 \(\bar{\boldsymbol{F}}\)，在单元刚度中则包含由该替换产生的附加项来组装切线矩阵。

FrontISTR 仅针对 8 节点线性六面体实现 F-bar 法，并可应用于小变形以及 Total Lagrange / Updated Lagrange 的非线性几何分析。

非协调单元¶

8 节点线性六面体单元没有足够的应变自由度来表示弯曲模态，因此在弯曲主导的问题中会出现弯曲锁定。非协调单元 [Taylor1976] 为弥补这一缺陷，在单元内部引入附加位移模态。

除单元节点位移 \(\boldsymbol{u}^e\) 外，在每个单元中引入仅存在于单元内部的非协调模态自由度 \(\boldsymbol{\alpha} \in \mathbb{R}^{9}\)，即 3 个方向 × 3 个模态，并将位移场近似为

\[ \boldsymbol{u}(\boldsymbol{r}) = \sum_{\alpha=1}^{8} N_\alpha^e(\boldsymbol{r})\, \boldsymbol{u}^e_\alpha + \sum_{k=1}^{3} M_k(\boldsymbol{r})\, \boldsymbol{\alpha}_k \]

对于自然坐标 \(\boldsymbol{r} = (\xi, \eta, \zeta)\)，非协调形状函数取 \(M_1 = 1 - \xi^2\)、\(M_2 = 1 - \eta^2\)、\(M_3 = 1 - \zeta^2\)。这些函数虽然不能保证单元边界上的连续性，但可在单元内部增加用于再现弯曲模态应变的空间。

单元刚度按外部自由度－内部自由度进行分块，

\[ \begin{bmatrix} \boldsymbol{K}_{dd} & \boldsymbol{K}_{d\alpha} \\ \boldsymbol{K}_{\alpha d} & \boldsymbol{K}_{\alpha\alpha} \end{bmatrix} \begin{bmatrix} d\boldsymbol{u}^e \\ d\boldsymbol{\alpha} \end{bmatrix} = \begin{bmatrix} \boldsymbol{F}^e_{\text{ext}} \\ \boldsymbol{0} \end{bmatrix} \]

组装后，通过 \(d\boldsymbol{\alpha} = -\boldsymbol{K}_{\alpha\alpha}^{-1}\boldsymbol{K}_{\alpha d}\,d\boldsymbol{u}^e\) 消除内部自由度，执行静态缩聚，从而得到仅含外部自由度的单元刚度

\[ \boldsymbol{K}^e = \boldsymbol{K}_{dd} - \boldsymbol{K}_{d\alpha}\,\boldsymbol{K}_{\alpha\alpha}^{-1}\,\boldsymbol{K}_{\alpha d} \]

并将其传递给整体组装。

FrontISTR 仅针对 8 节点线性六面体（C3D8IC）实现非协调单元，并可应用于小变形、Total Lagrange 和 Updated Lagrange。

U-P 混合单元¶

与 B-bar 法、F-bar 法在基于位移的框架中修正体积分量不同，u-p 混合（U-P）单元是一种将压力 \(\lambda\) 作为独立于位移的未知场引入的混合公式 [Bathe1996]。对于准不可压缩材料（例如泊松比极接近 0.5 的橡胶类材料或塑性变形后的金属），如果仅用位移场来满足体积恒定约束，就会产生体积锁定；将压力作为独立变量可缓解这一约束。

将应力分解为偏差分量和压力，

\[ \boldsymbol{\sigma} = \boldsymbol{\sigma}_{\mathrm{dev}} + \lambda\,\boldsymbol{I}, \qquad \boldsymbol{\sigma}_{\mathrm{dev}} = \mathbf{D}_{\mathrm{dev}}\,\boldsymbol{\varepsilon} \]

其中 \(\mathbf{D}_{\mathrm{dev}}\) 是从弹性矩阵中除去与体积模量 \(K\) 成比例的体积部分后得到的偏差弹性矩阵。压力 \(\lambda\) 与体积应变 \(g = \mathrm{tr}\,\boldsymbol{\varepsilon}\) 通过压缩率 \(\alpha^{-1} = 1/K\) 的约束条件

\[ g - \alpha^{-1}\lambda = 0 \]

相联系。对位移 \(\boldsymbol{u}\) 和压力 \(\lambda\) 作为未知量进行离散化后，可得到单元联立方程

\[ \begin{bmatrix} \mathbf{K}_{uu} & \mathbf{K}_{up} \\ \mathbf{K}_{up}^{T} & \mathbf{K}_{pp} \end{bmatrix} \begin{bmatrix} d\boldsymbol{u} \\ d\lambda \end{bmatrix} = \begin{bmatrix} \boldsymbol{f}_{u} \\ \boldsymbol{f}_{p} \end{bmatrix} \]

其中，\(\mathbf{K}_{uu}\) 包含偏差弹性分量以及（在有限变形时）几何刚度，\(\mathbf{K}_{up}\) 是耦合体积应变与压力的矩阵，\(\mathbf{K}_{pp} = -\int \alpha^{-1}\,\boldsymbol{N}_p \boldsymbol{N}_p^{T}\,dV\) 是压力稳定化项（\(\boldsymbol{N}_p\) 为压力形状函数）。由于压力自由度封闭在单元内部，因此可按

\[ \mathbf{K}_{\mathrm{eff}} = \mathbf{K}_{uu} - \mathbf{K}_{up}\,\mathbf{K}_{pp}^{-1}\,\mathbf{K}_{up}^{T} \]

进行静态缩聚，并将仅包含外部（位移）自由度的有效刚度传递给整体组装。

FrontISTR 仅针对 8 节点线性六面体实现 U-P 单元，每个单元设置 1 个压力自由度（单元内为常数）。该单元可应用于小变形、Total Lagrange 和 Updated Lagrange；在 Updated Lagrange 中，先用客观应力率（Jaumann／Hughes-Winget 型）更新偏差应力，然后由静态缩聚得到的值施加压力 \(\lambda\,\boldsymbol{I}\)。

壳单元¶

薄壁板壳结构使用基于 Reissner-Mindlin 板壳理论的壳单元。低阶基于位移的板壳单元随着厚度减小，会通过与体积锁定类似的机制（剪切锁定）高估横向剪切应变，从而使针对弯曲模态的刚度发散式增大。MITC（Mixed Interpolation of Tensorial Components）法 [Dvorkin1984] [Bathe1986] 仅在单元内预先规定的 tying point 对剪切应变分量重新采样，再将这些采样值插值回单元内部，从而避免该问题。

MITC 壳单元的节点位于曲面中面上，每个节点具有 6 个自由度，包括 3 个平移分量和绕中面法线方向的 3 个转动分量。单元刚度通过中面自然坐标与厚度方向合计三维的高斯积分进行评估，板厚 \(h\) 在本构计算时作为单元属性给定。

FrontISTR 提供中面按 1 层处理的 MITC3（单元 ID 731）、MITC4（741）、MITC9（743），以及沿板厚方向配置两层节点的分层壳单元 MITC3-shell361（761，3\(\times\)2 节点、每节点 3 自由度）和 MITC4-shell361（781，4\(\times\)2 节点、每节点 3 自由度）。在分层壳单元中，节点自由度仅包括 3 个平移分量，通过两层节点配置表示相当于转动自由度的弯曲模态。

梁单元¶

梁、框架结构等线材结构使用梁单元离散化。FrontISTR 采用考虑剪切变形的 Timoshenko 梁公式，将弯曲和剪切都表示为位移自由度与转动自由度的函数。

梁单元的节点自由度共有 6 个，包括 3 个平移分量以及绕轴线和横轴方向的 3 个转动分量，单元刚度通过沿梁轴向的一维数值积分进行评估。截面面积 \(A\) 与弯曲、扭转方向的截面二次矩 \(I\) 作为梁的截面常数赋予单元属性，并与材料的纵向弹性模量 \(E\)、剪切模量 \(G\) 一起构成拉伸、弯曲、扭转和剪切的各刚度系数。

FrontISTR 提供 2 节点直线梁单元（单元 ID 611），以及用 3 节点表示的 4 节点四面体实体－梁混合单元（641，用于混合自由度）。