基于旋转的第三介质接触正则化方法

《Computer Methods in Applied Mechanics and Engineering》：A rotation-based approach to third medium contact regularization

　　

第三介质接触方法利用虚构的“第三”介质隐式模拟接触相互作用。当接触发生时，虚构介质发生严重变形，这需要数值正则化以确保数值稳定性。最终，该正则化应促进良好的单元质量，但不应干扰接触力学建模。一种方法使用位移Hessian（Hessian矩阵）惩罚拉伸和旋转变形模式，这需要高阶单元。另一种方法引入额外自由度以惩罚旋转梯度的近似，这大大增加了系统规模。研究人员提出了一种基于虚构介质中旋转梯度近似的新正则化方法，该方法不惩罚拉伸变形模式，并且可用于一阶单元。通过多个数值示例（包括基准测试、第三介质中寄生力的研究以及一般加载条件的新应用）展示了所提出方法的有效性。

**研究背景与问题**
接触力学是固体力学中的核心问题，但接触现象固有的不连续性和未知的接触力与变形分布使其建模极具挑战性。传统接触方法（如节点-节点、节点-面或面-面方法）需预先定义接触面并依赖搜索算法追踪变形，自接触问题则需特殊处理。第三介质接触（TMC）方法通过引入虚构的“第三”介质隐式传递接触力，避免了搜索算法且简化了自接触建模，但面临数值不稳定问题：当第三介质被高度压缩时，有限元网格严重扭曲，导致算法终止。现有正则化方法如HuHu方法（惩罚位移Hessian矩阵）和HuHu-LuLu方法（修正拉普拉斯项）需高阶单元，而基于辅助场的旋转梯度惩罚方法（如Wriggers等人提出的）虽可避免拉伸惩罚，但引入额外自由度，显著增加系统规模。为此，研究人员提出一种基于局部旋转变化的新正则化方法，旨在仅惩罚旋转梯度而不惩罚拉伸变形模式，且兼容一阶单元，以降低计算成本并提高数值稳定性。

**研究内容与结论**
研究人员在论文《Computer Methods in Applied Mechanics and Engineering》中提出了一种基于旋转度量的正则化方法，通过惩罚单元内旋转与参考旋转的偏差来模拟旋转梯度，避免了拉伸变形模式的惩罚。该方法无需辅助场或高阶单元，可直接用于一阶单元。通过多个数值示例（包括平面应变C-Beam基准测试、C-Box对称加载、一般加载条件及三维球体接触）验证了其有效性。结果表明，与HuHu-LuLu方法相比，所提方法在拉伸变形模式下寄生效应更小，且与经典接触方法（如Abaqus中的罚函数法）一致性良好。该方法在一般加载条件下也能准确捕捉接触行为，并具有较低的计算开销（MFLOPs值仅为HuHu-LuLu方法的约1/13）。

**主要技术方法**
研究人员使用以下关键技术方法：
（1）基于局部旋转度量的正则化函数：通过度量两个旋转张量之间的Frobenius范数距离（?(R₁, R₂) = ‖R₁ - R₂‖_F），惩罚单元内各点旋转与单元质心参考旋转（取自上一收敛步）的偏差，从而近似旋转梯度。
（2）可压缩neo-Hookean材料模型：第三介质采用缩放后的体积应变能（?_vol = ?K(log_eJ)²）和等容应变能（?_iso = ?G(?₁ - 3)），并通过小参数γ_K、γ_G降低非接触时的刚度。
（3）一阶有限元离散：使用一阶单元（8节点六面体或4节点四面体），结合对称数值积分规则（如Gauss-Lobatto或Hammer积分），并采用自适应步长控制（初始步长Δλ₀，收敛失败时减半至Δλ_min）。
（4）并行求解框架：基于Gmsh生成网格，PETSc并行求解，MUMPS进行LU分解。

**研究结果**
**5.1 C-Beam in Plane Strain**
通过平面应变条件下的C-Beam基准测试（左端夹紧，顶部右端施加向下位移），对比了所提正则化方法与HuHu-LuLu方法。结果表明：所提方法（k_R = 8×10^-4 N/m²）可达到完全接触（位移1.0 mm），而HuHu-LuLu方法（k_F = 32×10^-7 N）最大位移约0.68 m。力-位移曲线显示，所提方法在接触前对结构刚度影响极小，而HuHu-LuLu方法随惩罚参数增大显著增加寄生力，导致提前接触。向上拉伸时（位移900 mm），所提方法引起的下臂上抬效应更小，且单元形状允许较大变化。参数研究显示，所提方法在宽参数范围内均能实现高负载水平，而HuHu-LuLu方法的最佳参数需权衡寄生效应与最大位移。计算效率方面，所提方法（3.60 MFLOPs）远低于HuHu-LuLu方法（46.8 MFLOPs），且对初始步长不敏感。

**5.2 C-Box in Plane Strain**
采用对称的C-Box示例（2 m×0.5 m矩形含1.8 m×0.3 m空腔，底部两端固定，中心施加向下1 m位移），使用所提正则化方法（γ_K=0, γ_G=2×10^-6, k_R=4×10^-4 N/m²）。结果显示，所提方法能准确捕捉对称加载下的接触行为，变形网格合理，验证了方法在对称工况下的有效性。

**5.3 C-Beam in general loading conditions**
将C-Beam扩展至一般加载条件（非平面，节点(0.35 m, 0.7 m, 1 m)施加u_y=-2 m, u_x=0.5 m），使用一阶4节点四面体单元（30,697个）。所提方法（k_R=4×10^-3 N/m²）成功模拟了复杂接触，尽管第三介质单元出现显著拉伸变形，但接触行为仍准确。与文献中六面体单元的结果相比，由于四面体单元弯曲刚度不同，略有差异，但整体接触行为一致。

**5.4 Spheres**
三维示例：两个球体（半径0.9 m和0.95 m，中心(0,0,0)和(0,2 m,0)）及半个椭圆圆柱（半长轴1.8 m，半短轴0.9 m）嵌入矩形域（2 m×4 m×3 m），顶部施加u_z=-0.9 m位移。使用379,527个一阶四面体单元，所提方法（γ_G=γ_K=10^-4, k_R=4×10^-3 N/m²）。结果显示，接触区域（靠近墙的球体）产生高静水压力，非接触区域（左侧球体）压力接近零，与Wriggers等人结果吻合。接触时左侧球体与圆柱间存在微小间隙（因第三介质刚度非零），但整体响应真实，表明所提方法在接触区域压缩响应与非接触区域张力响应中均表现良好。

**讨论与结论**
讨论部分强调，所提正则化方法仅基于局部旋转变化，避免了拉伸惩罚，因此第三介质网格可能严重扭曲，但数值示例表明接触行为仍精确捕捉，网格畸变仅为视觉假象。方法具有三大优势：（1）适用于非结构化网格；（2）直接兼容一阶单元；（3）避免低效的高阶积分方案和额外自由度。结论部分翻译如下：
研究人员提出的基于局部旋转变化的正则化方法，在平面和一般加载条件下的多个第三介质接触基准问题中均表现出有效性。此外，该方法与经典接触公式结果吻合良好。在平面加载条件下，数值示例表明，与HuHu-LuLu方法相比，所提方法寄生效应更小，尤其是拉伸变形模式，这符合其仅惩罚旋转的设计预期。最终数值示例显示，所提方法能达到与Wriggers等人相同的负载水平。所提方法的直接后果是有限元网格在第三介质中可能高度扭曲，但数值示例证明接触行为仍被准确捕捉，即网格畸变仅为视觉假象。该方法具有以下优点：1）适用于非结构化网格；2）直接兼容一阶单元；3）有效避免低效且昂贵的积分方案及额外自由度的增加。