具有变系数的稳定型与自稳定型p-虚拟元素方法基准测试

《Computer Methods in Applied Mechanics and Engineering》：Benchmarking stabilized and self-stabilized p-virtual element methods with variable coefficients

　　帕奥拉·皮亚·福利尼奥|达尼埃莱·博菲|法比奥·克雷达利|里卡多·韦斯科维尼
米兰理工大学航空航天科学与技术系，意大利米兰20156

摘要
标准虚拟元素方法基于多项式投影，为评估离散空间中非多项式部分的贡献，需要引入稳定项。然而，该稳定项并非由底层变分形式唯一确定，通常以特定方式加入，可能会影响数值结果。为解决这一问题，人们提出了无稳定项和自稳定公式，但其理论分析仍不够成熟。本文对p版本虚拟元素方法的各类稳定及自稳定公式进行了深入的数值研究。结果表明，自稳定和无稳定公式虽能实现最优精度，但条件数较差。此外，本文在标准虚拟元素空间框架内引入了一种新的投影算子，可明确处理变系数问题。数值实验表明，对于较大的p值，这种新方法比现有方法更具鲁棒性。

引言
当前的研究致力于设计效率更高的新型航空航天结构。例如，变刚度面板就是一种此类结构，它由具有变向纤维的变刚度蒙皮与沿曲线路径分布的加劲肋组成。这类结构可通过具有非均匀弹性特性的偏微分方程来建模。针对这一应用，已有基于有限元法[1]和里兹法[2]的数值研究方法。但在有限元框架下，板材与加劲肋之间的节点共享需要精细的网格划分或特殊公式[3]。而里兹法则简化了板材与加劲肋之间的相容性约束，但在处理复杂几何形状时通常不够灵活，除非采用特定程序[4]。

在此背景下，虚拟元素方法[5]是一种很有前景的工具，因为它允许使用由一般多面体元素构成的网格，这些元素的边甚至可以是弯曲的[6]。因此，可轻松生成与加劲肋形状相匹配的网格，并结合高阶离散化方法使用。尽管虚拟元素方法已应用于计算力学中的诸多问题，如平面应力弹性问题[7][8][9][10][11]、屈曲问题[12][13][14]、基尔霍夫-洛夫问题[15]以及赖斯纳-明德林板模型[16][17]，但据作者所知，其在变刚度面板中的应用仍然有限。文献[18]中提出了对变系数进行分段常数近似的方法。在基于赫林格-赖斯纳变分原理的虚拟元素方法框架下，也有研究探讨了具有空间变化材料系数的弹性问题。在[19]的二维研究中，定义了一个局部能量投影算子，从而得到了在存在变弹性系数时不再是多项式的元素级近似空间。[20]也提出了类似方法，其公式为三维形式，且同样采用了基于投影的方法处理多面体网格。本文则旨在通过考虑变刚度面板带来的挑战，对p版本虚拟元素方法进行初步研究。

虚拟元素方法可视为在多面体网格上对有限元方法的推广。但由于试验函数和测试函数都是局部偏微分方程的解（在每个网格元素中定义），且在构建离散公式时并未明确求解这些方程，因此虚拟元素方法属于非一致化方法，这也正是“虚拟”一词的由来。由于仅能通过自由度对试验函数进行多项式投影，因此需引入稳定项来处理剩余的非多项式部分，确保离散问题的稳定性。不过，稳定项具有任意性，通常并不符合问题的物理规律。在变系数情况下，该稳定项通常通过常数近似进行缩放，而多项式投影则要么独立于系数构建[22]，要么基于该系数的低阶多项式近似来构建[23]。

稳定化过程需要引入用户定义的缩放参数，但这会降低方法的通用性，可能导致过度稳定化，进而降低解的质量。在各向异性问题[24]和特征值问题[25][26]中，若稳定参数选择不当，就可能出现虚假模态。关于稳定参数与虚拟元素方法性能之间的关系，已有相关研究[27][28]。尽管文献中存在多种稳定化形式[5][29][30][31][32]，但它们大多适用于低阶虚拟元素，而随着阶数的增加，其效果往往不再理想。

为解决稳定项带来的问题，有研究提出了无稳定项的虚拟元素[24][33][34]。这类元素是通过将虚拟函数投影到更大的多项式空间中得到的，因此所得离散问题具有自动稳定性。另一种策略是所谓自稳定元素，它是通过在近似空间中增加额外的内部自由度来实现的[35]。本文中，“自稳定虚拟元素方法”将指代这两种方法。虽然本文未涉及，但文献[36]提出了一种降阶方法，可用于低成本重构最低阶的虚拟元素方法基函数，从而避免使用稳定项。

本文共分为两个主要部分。首先，在p版本虚拟元素方法框架下，系统评估了几种稳定化及自稳定公式的性能。为此，研究了拉普拉斯问题、线性弹性问题和斯托克斯问题。目的是了解在结构化网格、扭曲网格甚至带弯曲边的网格上，这些方法的精度、稳定性、条件数以及对稳定参数的敏感性。据作者所知，目前尚无针对p版本虚拟元素方法的此类研究。

本文的第二个创新点在于提出了一种处理高阶近似中变系数问题的新方法。在现有的虚拟元素空间框架内，引入了一种新的多项式投影方式，以便考虑系数的特性。这种被称为VC-VEM的公式可轻松融入稳定化及自稳定框架中，从而提升精度。

本文的结构如下：第2节介绍几何框架；第3节回顾了利用稳定化及自稳定方法对拉普拉斯方程进行虚拟元素离散化的方法；第4节和第5节分别讨论线性弹性问题和斯托克斯问题；第6节通过大量数值测试对p版本虚拟元素方法进行性能评估；第7节介绍了用于二阶椭圆型问题和线性弹性问题的VC-VEM方法；第8节展示数值结果，证明所提元素的鲁棒性；第9节给出结论。

片段摘录

几何框架
研究域记为Ω，{Γi}i=1,?,N表示构成边界Γ??Ω的有限组光滑曲线，Ne为边界曲线总数。假设边界Γ是利普希茨连续的，且对于m≥0，每条曲线Γi都属于Cm+1类[6]，则可定义关于区间Ii=[ai,bi]的Cm+1类正则且可逆参数化γi: Ii→Γi。由于γi是可逆映射，对于任意函数g，都有g=g°γi?g∈Γi，以及g=g°γi?1?g∈Ii。在本研究中，弯曲边是通过贝塞尔曲线进行参数化的。

拉普拉斯问题
本节介绍了用于拉普拉斯问题的虚拟元素方法的稳定化及自稳定公式，其弱形式为：
寻找u∈H01(Ω)，使得：
a(δu,u)=(δu,f)Ω?δu∈H01(Ω)，
其中a(u, δu)表示内能，(δu, f)Ω表示外部作用力所做的功：
a(δu,u)=∫Ω?δu·?u dΩ，(δu,f)Ω=∫Ωδu·f dΩ。
文中使用紧凑符号(·,·)D表示L2(D)空间中的标量积。

随后定义离散空间Vh?H01(Ω)。

线性弹性问题
用于拉普拉斯问题的稳定化及自稳定公式可扩展到线性弹性问题。为简化问题，此处考虑二维情况。两个平面位移分量定义为：u(x,y)=[u(x,y)v(x,y)]，在微小位移假设下，应变-位移关系为：
??(u)=[εxxεyyγxy]=[u,xv,yu,y+v,x]，
其中(·),x和(·),y分别表示对坐标x和y的导数。

斯托克斯问题
斯托克斯问题是一个混合问题，未知量为速度u和压力ρ，其弱形式为：
寻找(u,ρ)∈[H01(Ω)]2×L02(Ω)，使得：
a(δu,u)+b(δu,ρ)=(δu,f)Ω?δu∈[H01(Ω)]2，b(δρ,u)=0?δρ∈L02(Ω)，
其中：
a(δu,u)=∫Ων?u·?δu dΩ，b(δu,ρ)=∫Ωdiv(δu)ρ dΩ，(δu,f)Ω=∫Ωδu·f dΩ，
ν为粘度，L02(Ω)是L2(Ω)中均值为零的函数所构成的子空间。在本研究中，采用文献[37][38]中描述的无散近似。因此，当p≥2且p≤m时，局部速度空间为……

数值结果：稳定化及自稳定虚拟元素方法的性能评估
本节从精度和条件数两个方面，比较了不同稳定化技术和自稳定公式的性能。众所周知，稳定项可能会影响高阶近似的计算结果，导致收敛性不佳[21]。由于本文的重点是p版本虚拟元素方法，因此本节旨在阐明稳定参数和稳定化技术对该方法精度和稳定性的影响。

含变系数的四元方程
许多实际问题中的系数会随空间位置发生变化。现代复合层压材料就是一个典型例子，其弹性特性取决于位置。对于大多数数值方法，通常会在多项式空间中对通用系数进行插值。如果采用足够高的多项式阶数或足够细的网格，这种近似方法是可行的。在虚拟元素方法中，关于如何处理变系数的问题已在[22]中得到探讨。

数值结果：含变系数的方程
本节针对具有变系数的一组问题进行了数值测试，旨在验证前一节中提出的方法的有效性，特别是分析系数变化对方法精度的影响。

结论
本文对p版本虚拟元素方法的稳定化及自稳定公式进行了深入比较。研究过程中充分考虑了在具有曲线形加劲肋的变刚度面板领域的潜在应用。在这种情况下，本构方程中的系数是变化的，为了准确描述几何形状，其离散化可能需要使用曲线形网格。在本文的第一部分，通过性能对比……

CRediT作者贡献声明
帕奥拉·皮亚·福利尼奥：写作——审稿与编辑，写作——初稿撰写，可视化，验证，软件，方法论，研究，概念设计。
达尼埃莱·博菲：写作——审稿与编辑，指导，资源获取，资金筹措，概念设计。
法比奥·克雷达利：写作——审稿与编辑，写作——初稿撰写，指导，研究，概念设计。
里卡多·韦斯科维尼：写作——审稿与编辑，指导，资源获取，概念设计。

利益冲突声明
作者声明不存在任何可能影响本文研究成果的已知财务利益或个人关系。