解决体积闭锁现象：用于接触问题的无体积闭锁混合虚拟单元方法

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：Volumetric locking-free mixed Virtual Element Methods for contact problems

本文针对二维无摩擦弹性接触问题，在近不可压缩材料中出现的“体积闭锁”（volumetric locking）数值现象，提出了一种基于混合变分公式（位移/压力，u/p）的虚拟单元法（Virtual Element Method, VEM）。研究团队开发了通用的误差分析框架，考察了含短边的网格，并分别构造了一阶和二阶VEM格式。数值结果表明，两种格式均能有效避免体积闭锁，对短边网格鲁棒，并能达到理论预期的收敛速率，为复杂网格下接触问题的精确、稳健求解提供了有效工具。

　　

在计算固体力学领域，弹性体的接触问题是一个经典且具有挑战性的课题。从简单的机械咬合到复杂的生物组织相互作用，理解接触过程中的应力分布和变形至关重要。然而，当使用有限元法等数值工具进行模拟时，工程师和科学家们常会遇到一个棘手的“拦路虎”——体积闭锁（volumetric locking）。这种现象尤其容易发生在模拟近乎不可压缩的材料时，比如橡胶、生物软组织或处于塑性流动状态的金属。简单来说，体积闭锁会导致数值解“过度僵硬”，计算出的位移场严重失真，即使网格不断加密，精度也难以提升，这极大地限制了数值模拟的可靠性和应用范围。另一方面，实际工程问题中的接触界面往往形状复杂，生成规则、匹配的网格（如全部由三角形或四边形组成）非常困难且耗时。虚拟单元法（Virtual Element Method, VEM）的出现为这一困境带来了曙光。作为有限元法（FEM）和模仿有限差分法（Mimetic Finite Differences）的“进化版”，VEM继承了FEM的变分框架，同时能极其灵活地处理多边形甚至多面体网格，轻松应对悬挂节点和非匹配网格，这使其在处理复杂接触界面时具有天然优势。但是，将VEM应用于接触问题，特别是在近不可压缩 regime 下，如何同时保证方法的准确性（避免体积闭锁）和稳健性（应对网格中的短边等不利情况），仍需深入的理论探索和算法设计。为此，来自米兰大学（University of Milano）的C. Lovadina和L. Molinari在期刊《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》上发表了他们的研究成果，致力于发展一种能够“免疫”体积闭锁的混合虚拟单元法，用于求解二维无摩擦接触问题。

为了攻克体积闭锁难题并利用VEM处理复杂网格的能力，研究人员在本研究中主要运用了几个关键技术方法。首先是混合变分公式的采用。研究放弃了传统的、只以位移为未知量的公式，转而采用了位移/压力（u/p）混合形式，将压力（与体积应变相关）显式地作为一个独立的未知量引入。这种方法是处理近不可压缩材料的经典策略之一。其次是虚拟单元法（VEM）的离散框架。研究在VEM的框架下构建离散的位移和压力空间，其中位移空间包含直到k阶的多项式。关键步骤是设计了离散的双线性形式a_h(·,·)来近似表示连续的双线性形式a(·,·)（与剪切变形能相关），该设计满足k阶一致性（consistency）和稳定性（stability）条件。最后是节点插入算法与接触匹配条件的处理。为了便于在接触边界上施加离散的、点对点的接触条件，研究采用了来自文献的节点插入算法来处理可能非匹配的接触边界网格，从而得到一个统一的、节点匹配的接触边界离散。这个过程会自然产生包含“短边”（small edges）的网格，研究在理论分析中明确考虑了这种情形的影响。

本研究的主要内容和结论通过以下部分展开：

1. 接触问题的混合变分形式

研究首先从二维线弹性力学的基本方程出发，考虑两个弹性体的无摩擦接触（Signorini条件）。为了规避体积闭锁，作者没有采用标准的位移法公式，而是引入了以位移u和压力p（定义为pⁱ= λⁱdiv uⁱ）为独立变量的混合变分形式。该形式最终归结为求解一个定义在凸集K（由不可侵入接触条件定义）上的变分不等式系统。研究证明了该混合问题解的存在唯一性及其关于Lamé参数λ的一致稳定性。

2. 虚拟单元离散化

研究详细描述了VEM对上述混合变分系统的离散化。对于位移场，在每一个多边形单元K上构建了VEM空间V_hⁱ，其包含直到k阶的向量多项式。对于压力场，则采用分片常数或高阶多项式空间Q_hⁱ。对于接触条件的处理，研究假设通过节点插入算法后，接触边界Γ_C拥有一个统一的、节点匹配的离散化，从而可以方便地定义离散的凸集K_h，其由在接触边界节点上满足不可侵入条件的离散函数构成。离散的核心在于构造局部双线性形式a_h^K(·,·)，它由两部分组成：一是对位移能量投影算子的精确计算，二是对一个稳定项的设计，以确保离散格式满足一致性（对多项式精确）和稳定性（其常数α_{?, hK}和α^?_hK可能依赖于网格尺寸）。

3. 一般收敛性分析与两种示例格式

研究建立了一个通用的误差分析框架。关键假设包括离散的inf-sup条件成立（常数β_h可能依赖于h），以及离散空间具有良好的逼近性质。在此框架下，论文给出了位移误差和压力误差的抽象估计。作为该框架的具体实现，论文提出了两种VEM格式：一种具有一阶（k=1）一致性和逼近性质；另一种则具有二阶（k=2）性质，其中位移空间采用了更复杂的二阶VEM空间。论文指出，对于接触问题，即使解足够光滑，二阶格式的理论收敛阶也存在一个5/2的上限，这是由于接触凸集离散化引入的额外一致性误差导致的，这是接触问题的典型特征。

4. 对二阶格式的详细分析

论文第五部分专门针对二阶格式进行了深入分析。重点估计了在网格可能包含“短边”的情况下，离散双线性形式a_h(·,·)的稳定性常数α_?,h、α^?_h以及离散inf-sup常数β_h。分析表明，这些常数可能以对数形式依赖于网格的最小边尺寸h_m。尽管存在这种依赖，但由于对数增长缓慢，在实际计算中其导致的退化并不明显。

5. 数值实验

最后的数值部分通过算例验证了理论预测。结果表明：（1）无论是二阶格式还是一阶格式，其误差均与Lamé参数λ无关，即完全消除了体积闭锁现象；（2）即使在人为引入的、包含极短边的网格上，方法依然表现稳健，误差并未出现显著恶化；（3）两种格式均达到了其理论预期的收敛速率，一阶格式关于位移和压力的L²误差收敛阶约为2，能量误差收敛阶约为1；二阶格式的位移和压力L²误差收敛阶约为3，能量误差收敛阶约为2，与理论分析相符。

结论与意义

本研究成功地将混合变分公式与虚拟单元法相结合，为二维无摩擦弹性接触问题发展了一套稳健、高效的数值离散方案。其最主要的贡献在于彻底解决了近不可压缩材料模拟中的体积闭锁难题，同时充分发挥了VEM处理任意多边形网格的灵活性，使得对具有复杂接触界面和非匹配网格的问题进行高精度分析成为可能。理论分析保持了高度的严谨性，明确追踪了各稳定性常数对网格尺寸（特别是短边）的依赖关系，为方法在实际复杂场景中的应用提供了可靠性保证。论文中发展的分析框架具有一般性，为后续研究更复杂的接触问题（如摩擦、大变形、三维问题）奠定了基础。这项工作不仅推动了虚拟单元法在计算接触力学领域的发展，也为工程中广泛存在的、涉及软材料或大变形接触的实际问题提供了更强大的数值模拟工具。