《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：A high-order stabilization-free Virtual Element Method for 3D linear elasticity
编辑推荐：
在计算力学领域，传统的虚拟单元法（VEM）依赖于数值稳定化项以确保收敛，这构成了其应用和性能的一个显著限制。为了解决这一问题，Timothée Bouchez 等人基于胡-鹫津（Hu-Washizu）变分框架，首次提出了一种适用于三维线性弹性问题的高阶无稳定化虚拟单元法（sf-VEM）。该方法采用多项式插值离散应变和应力场，并利用涉及偶然性投影算子的增强虚拟单元空间来离散位移，从而在保证稳定性的同时有效减少了自由度数量。研究通过在不同网格类型上进行的数值分析验证了其良好的收敛性能，为复杂几何结构的高效、高精度仿真提供了新工具。
在计算固体力学的广阔天地里，科学家和工程师们一直致力于发展能够高效、精确模拟材料行为的数值方法。其中，有限元法（FEM）长期占据主导地位，但其对网格形状（如需要凸多边形/多面体）的严格要求，有时会成为处理复杂几何结构（例如来自工业设计的异形零件或生物组织的多孔结构）的“绊脚石”。想象一下，当你试图用标准的、形状规则的“积木”去拼出一个极其不规则的物体时，会多么费力和低效。为了突破这一限制，虚拟单元法（Virtual Element Method, VEM）应运而生。它像一位“几何魔术师”，能够优雅地处理星形凸多面体网格，大大放宽了对网格几何形状的限制，从而可以使用更快速、更鲁棒的网格生成算法（如Voronoi镶嵌）。这为构建工业级仿真工具带来了巨大希望。
自2013年提出以来，VEM已在多个工程领域展现出潜力，从基础的拉普拉斯方程到复杂的线性弹性、超弹性、弹塑性乃至接触问题。然而，这位“魔术师”有一个与生俱来的“小瑕疵”：它需要依赖一个称为“稳定化项”的数值技巧来构建其离散的双线性形式。这个稳定化项本质上是数值性的，并非物理定律的直接体现，因此需要针对不同的问题精心选择和设计，以确保方法的收敛速度和精度。这无疑增加了VEM的复杂性和应用门槛，成为其发展的一个显著制约。
自2021年起，为了摆脱这个“包袱”，无稳定化虚拟单元法（stabilization-free VEM, sf-VEM）开始崭露头角。这类方法通过将感兴趣场的梯度高阶投影到一个多项式空间，并精心选择投影的阶次来确保稳定性，从而完全避免了额外稳定化项的需要。不过，相关研究大多集中于二维问题，例如二维泊松方程、平面弹性、有限应变等。而对于在实际工程中至关重要的三维问题，sf-VEM的应用则相对较新且有限，目前主要是一阶单元在线性弹性和超弹性中的应用。然而，在连续介质力学中，高阶单元往往是必需的，因为低阶单元在收敛到真实解时可能遇到困难，例如在薄壁结构分析中，一阶单元的“闭锁”行为就是一个重大挑战。因此，开发一种三维高阶无稳定化虚拟单元法，对于推动该方法的实际应用、提升复杂结构仿真能力具有迫切需求。
这正是由Timothée Bouchez、Anthony Gravouil、Nawfal Blal、Anthony Giacoma、Emmanuel Delor和Jean-Daniel Beley（来自法国INSA Lyon, CNRS, LaMCoS实验室）组成的研究团队在论文《A high-order stabilization-free Virtual Element Method for 3D linear elasticity》中所完成的工作。该论文发表在《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》上，首次在文献中提出了一种针对三维线性弹性问题的高阶无稳定化虚拟单元法。
为了开展这项研究，作者们主要运用了以下几个关键技术方法：首先，他们基于经典的胡-鹫津（Hu-Washizu）变分原理来表述线弹性问题，这为同时独立离散位移、应变和应力场提供了自然的框架。其次，在将计算域用非重叠的星形凸多面体单元进行划分后，他们将应变和应力的离散函数空间构建为多项式空间，并通过选择特定的多项式阶数来保证公式的稳定性。第三，他们构建了一个增强的虚拟单元空间用于离散位移，其中引入了一个“偶然性投影算子”（serendipity projection operator）。这个算子的选择旨在减少多项式积分步骤以及位于多面体各面上的自由度数量，从而优化计算效率。最后，他们通过数值实验，在多种网格类型上系统分析了所提方法的收敛行为。
研究结果
1. 方法构建与离散公式
研究基于胡-鹫津三场变分原理建立离散公式。该原理将线弹性问题的解表达为使一个包含位移(u)、应变(ε)和应力(σ)的泛函取驻值。通过对计算域Ω进行多面体网格划分，分别为这三个场构建了离散的近似空间：应变和应力空间被选为特定阶次的多项式空间，其阶次经过选择以确保整个公式的稳定性；位移空间则被构建为一个增强的虚拟单元空间，它依赖于一个定义的偶然性投影算子，该算子有助于减少位于单元面上的自由度总数，从而提升计算效率。最终，离散系统由平衡方程、本构关系和协调方程这三个生成子（generator）构成。
2. 网格与坐标系统
方法适用于满足特定弱条件的非重叠星形凸多面体网格。每个多面体单元E需要关于一个球是星形凸的，其每个多边形面F也需要关于一个圆盘是星形凸的。研究中引入了三套坐标系：全局坐标系X，基于全局正交基B0= {e1, e2, e3}；每个单元E的局部无量纲坐标系ξ = (X - XE)/hE；以及每个平面面F的二维局部无量纲坐标系ξ?，定义在其自身选择的正交基BF上。这些坐标系便于后续多项式的定义和计算。
3. 多项式空间与投影算子
研究定义了所需的多项式空间。对于单元的面F，定义了其上至k阶的三维向量值多项式空间[Pk(F)]3，并给出了基于局部坐标ξ?的显式基。对于单元E本身，定义了其上至k阶的三维向量值多项式空间[Pk(E)]3，其基由局部坐标ξ的单项式与全局基向量张量积构成。此外，还引入了用于应变和应力近似的对称张量值多项式空间。核心创新之一是为位移虚拟单元空间引入了一个“偶然性投影算子”ΠS。该算子将位移场投影到一个多项式子空间，其定义使得对于足够光滑的函数，投影与局部L2投影一致，从而减少了计算中对内部多项式模式进行积分的需求，并最终减少了位于单元面上的自由度。
4. 离散函数空间的具体构造
基于上述框架，研究具体构造了三个离散函数空间VEε, VEσ和 VEk,l，分别用于离散应变、应力和位移。
- •
应变空间VEε：被直接取为对称张量值多项式空间[Pmε(E)]3×3sym，其中mε是选定的应变多项式阶数。
- •
应力空间VEσ：被直接取为对称张量值多项式空间[Pmσ(E)]3×3sym，其中mσ是选定的应力多项式阶数。应变和应力多项式阶数的选择需满足稳定性条件。
- •
位移空间VEk,l：这是一个增强的虚拟单元空间。其函数在单元E内部没有显式表达式，但通过其在单元边界（面和边）上的值和矩来定义。具体来说，其自由度包括：在所有顶点处的位移值；在所有边e上，对至l-2阶的边多项式模式的矩；在所有面F上，对至k-2阶的面多项式模式（来自空间[Pk(F)]3）的矩。此外，位移函数还需要满足一个关键条件：其偶然性投影ΠSkuh在单元内部与一个由面矩确定的特定多项式相匹配。空间维数可以通过计数这些自由度并考虑约束关系得到。
5. 数值实现与收敛性分析
研究概述了所提公式的实现步骤，主要包括：构建每个单元的局部线性系统（对应于离散的平衡、本构和协调方程），然后组装成全局系统并施加边界条件进行求解。为了验证方法的有效性，研究在第7节进行了系统的数值收敛性分析。分析针对一个具有已知解析解的三维线弹性问题，在四种不同类型的网格上进行：规则立方体网格、扰动立方体网格、Voronoi网格以及由立方体和棱锥混合构成的非凸多面体网格。数值实验展示了当同时提升位移阶次k和应变/应力阶次(mε, mσ)时，位移误差和应变误差在能量范数和L2范数下均实现了最优收敛阶。结果证明了所提出的高阶无稳定化VEM对于多种复杂的单元几何形状都具有良好的鲁棒性和收敛性能。
研究结论与意义
本研究成功提出并数值验证了一种适用于三维线性弹性问题的高阶无稳定化虚拟单元法。该方法的核心贡献和意义在于：
- 1.
首次实现三维高阶sf-VEM：在作者已知的范围内，这是文献中首次针对三维线性弹性问题提出的高阶无稳定化虚拟单元法，填补了该领域的研究空白。
- 2.
创新的公式构建：通过采用胡-鹫津三场变分框架，为应变和应力独立地采用多项式插值，并为位移设计一个包含偶然性投影算子的增强虚拟单元空间，巧妙地绕过了对传统VEM稳定化项的依赖。应变和应力多项式阶数的选择直接用于保证数值稳定性。
- 3.
计算效率优化：引入的偶然性投影算子不仅具有清晰的数学意义，其设计还旨在减少多项式积分运算和位于单元面上的自由度数量，有助于控制总体计算规模。
- 4.
鲁棒性与通用性：数值实验表明，该方法在多种网格类型上，包括规则的、扰动的、Voronoi以及非凸多面体网格，都表现出了最优的收敛阶。这证明了该方法对复杂几何形状具有良好的适应性和鲁棒性，充分发挥了VEM在处理多面体网格方面的固有优势。
- 5.
重要的工程应用前景：该方法的成功开发，为分析薄壁结构、复合材料、生物组织等涉及复杂几何和可能出现的数值闭锁问题的工程领域，提供了一个新的、无需担忧稳定化项设计的高阶数值工具。它推动了无稳定化VEM从二维向更贴合工程实际的三维问题拓展，为发展更高效、更灵活的工业仿真软件奠定了重要的算法基础。
