带罚项的Morley型虚拟元方法求解四阶椭圆型奇异摄动问题

《Communications in Nonlinear Science and Numerical Simulation》：The Morley-type virtual element method with penalty for a fourth-order elliptic singular perturbation problem

　　

摘要
我们发展了一种带罚项的非协调虚拟单元方法（VEM），用于求解四阶椭圆奇异摄动问题的数值解。该方法采用Morley型虚拟单元离散化方案，能够自然地适配多边形和多面体网格，从而便于在几何复杂域上进行计算。由于标准的Morley型单元在低阶微分算子的逼近中已知存在失效问题，我们通过引入内部罚函数策略对离散格式进行了修正。我们的主要贡献包含两个互补部分：一种精心设计的罚函数机制，可增强方法对参数的鲁棒性；以及一种计算高效的格式，保留了虚拟单元框架的优良特性。我们建立了所提方案的最优收敛率和ε一致误差界。我们通过一系列数值试验验证了理论预测，并展示了带罚项的Morley型VEM的实际性能。

引言
我们考虑有界多边形域Ω??2上的四阶椭圆奇异摄动问题：

{??2Δ2???Δ??=??, 在Ω内,
??=????????=0, 在?Ω上,

其中n是边界?Ω的单位外法向量，f∈L2(Ω)为右端源项，ε>0为小的正常摄动参数。

问题(1.1)描述了二维薄屈曲板的挠曲[1]；它同时也是完全非线性Monge-Ampère方程的线性化形式[2]。当(1.1)中ε→0时，我们形式上恢复Poisson极限。在这种奇异情况下，精确解在靠近边界的一条薄带内表现出急剧的空间过渡，即文献[3]中报告的边界层效应。这种截然不同的尺度的共存，从根本上挑战了标准数值方案在离散层面以及后续收敛分析中的可靠性。

在有限元方法（FEM）框架内，H2协调单元，包括Argyris单元[4]及其三维推广[5]，以及Hsieh-Clough-Tocher单元[6]，能够同时逼近(1.1)中的两个算子，详见[7]。然而，这些单元依赖高阶多项式、复杂的空间结构和大量自由度（DoF），使得它们在处理以边界层现象为主导的问题时相当笨重。鉴于这些困难，非协调有限元方案为数值处理四阶奇异摄动椭圆方程提供了一种有吸引力的替代方案。其主要优势源于网格界面间松弛的光滑性约束，从而在保持整体逼近精度的同时大幅减少自由度数。对于问题(1.1)，文献[3]、[8]、[9]、[10]、[11]、[12]、[13]、[14]、[15]、[16]中构建了各种H2非协调但H1协调的单元，而完全H2非协调的单元也在[13]、[17]、[18]、[19]、[20]中得到研究。Brenner等人[21]发展了一种间断Galerkin离散化方案，其中包含内部罚函数（IPDG）并保持了C?连续性；在仅施加最小正则性假设的条件下，他们建立了参数鲁棒的误差界。该技术随后被Franz等人[22]推广至边界层适配网格。许多其他针对四阶椭圆问题的间断Galerkin方案[23]、[24]、[25]、[26]、[27]、[28]、[29]、[30]、[31]、[32]、[33]、[34]、[35]也可以自然地适配离散化问题(1.1)。

在早期为四阶边界值问题开发的非协调单元中，Morley-Wang-Xu（MWX）单元[17]、[18]以其在每个多边形上最少自由度的复杂性而闻名。然而，后续研究揭示了一个关键缺陷：该单元对二阶椭圆问题无法收敛，如文献[3]、[36]所记录。这一否定性结果立即排除了MWX单元对我们目标方程(1.1)的直接使用，尤其是在ε→0的极限情况下，此时方程变为奇异摄动形式。为弥补这一不足，出现了三类设计离散方案的策略，其本质差异在于低阶项和右端项的离散方式。

第一类引入一个辅助插值算子来修正与Laplace算子相关的双线性形式[37]、[38]，使得双调和项和低阶项在同一非协调框架内得到统一处理。尽管这些方法保证了收敛性，但它们在Poisson极限ε=0时退化为低阶方案。

第二类采用IPDG技术来离散化Laplace算子[39]，其中双调和项保持非协调离散化，而低阶项通过IP技术进行稳定化。这些方法在不引入额外插值的情况下实现了对ε的鲁棒收敛，但跳跃项使得刚度矩阵变为稠密矩阵。

第三类采用右端修正方法[40]：保持标准MWX离散双线性形式不变，将右端项(f, vh)替换为(f, Phvh)，其中Ph是到Lagrange单元空间的H1正交投影。该方法达到了最优收敛性，并进一步通过Nitsche技术改善了边界层附近的收敛行为，但其依赖于辅助协调空间的构建。

上述三类策略均在有限元框架内开发，其中网格受限于特殊剖分。这一情况促使我们探索一种新的离散框架，该框架既保留MWX单元的低自由度优势，又能适配任意多边形网格。

为克服这一困境，我们转向虚拟单元方法（VEM）。与MWX单元在奇异摄动设置下存在严重缺陷不同，VEM保留了标准FEM的变分框架，并消除了传统的几何约束：由任意多边形和多面体组成的网格自然可被接受[41]、[42]。这种固有的几何适应性，反过来为处理定义在复杂物理域上的问题建立了一个灵活的数值环境。其核心思想是通过投影算子将自由度映射到多项式空间，并通过稳定化项确保稳定性，从而避免了显式构建形函数。近年来，VEM已被成功应用于四阶椭圆方程。特别是对于奇异摄动问题，已有一些关于ε→0时退化行为及边界层现象处理的重要研究[43]、[44]、[45]。

值得注意的是，非协调VEM与IP（IPVEM）技术的结合，可以在保持计算精度的同时放宽连续性要求。

IPVEM文献中也出现了针对奇异摄动四阶问题的贡献。Zhang和Zhao[44]发展了一种带内部罚函数的非协调虚拟单元离散化，证明了最低阶单元在奇异摄动参数下的一致收敛性。此外，根据第1种策略，Feng等人在[43]中通过引入协调插值算子修改了IPVEM，获得了最低阶单元的一致误差估计。

受文献[39]中IP-MWX有限元的启发，在本文中我们构建了一种带罚项的Morley型虚拟单元方法，用于在任意多边形网格上求解(1.1)。离散双线性形式通过VEM投影算子ΠΔ??和Π???计算。通过采用经典定义的Morley型自由度，并对函数值在单元间的跳跃施加罚项，该方法保留了Morley型VEM的最小自由度优势，并在ε→0的退化情况下实现了鲁棒性，在边界层主导的情况下达到了h^(1/2)的一致收敛率。

我们的构造从文献[42]中引入的H2协调虚拟单元空间V????(K)出发。通过两个缩减操作获得修正的Morley型子空间Vh(K)：通过关系式r?g(v)=r?g(ΠΔ??v)消除顶点处的梯度相关自由度，并将内部矩基限制为阶数至多k?4的多项式。剩余的有效自由度为顶点求值、边上的矩（针对势函数及其法向通量），以及低阶内部矩。这种有效自由度的特定选择对投影算子的可计算性有两个直接后果。首先，由于条件r?g(v)=r?g(ΠΔ??v)被子空间继承，H2投影算子ΠΔ??可以从剩余自由度中唯一恢复。其次，法向导数矩χn的存在允许通过Green-Gauss恒等式直接计算H1型投影Π???，完全避免了任何辅助求积规则。这确保了所有刚度矩阵的组装完全依赖于我们定义的自由度，保证了算法在实际实现中的完全可计算性。特别是对于最低阶单元k=2，自由度仅为2n，仅由顶点值和边上的矩组成，完全独立于顶点梯度或边上插值点。这种构造确保了空间的维数是最小的，同时所有自由度都具有明确的物理解释并且在计算上是可行的。

由于Morley型空间是完全非协调的（函数在单元界面间不连续），受文献[39]中思想的启发，我们引入了一种IP格式，其中包括一个函数跳跃罚项J1和两个对称化项J2、J3，以弱方式强制C?连续性，从而为低阶Laplace算子提供必要的稳定性。

值得注意的是，我们对单元间断裂的处理与早期IPVEM工作[44]、[46]中采用的方法不同。在那些工作中，边上的法向导数跳跃受到罚项作用，而在本文的格式中，只有函数值跳跃[[v]]进入罚项。后者的处理方式直接受到协调有限元框架中Morley单元的经典IP对应版本的启发[39]；它同时也带来了一个实际优势——罚参数更容易调节而不会降低收敛行为。

我们证明了，对于足够大的罚参数σ，离散方案同时满足连续性和强制性，从而保证了适定性。在误差分析中，我们推导了最优收敛率O(h^(k?1))。对于边界层主导的情况，我们为最低阶情形k=2建立了ε一致的半阶收敛率O(h^(1/2))。最后，我们在非结构多边形网格上进行了数值实验以验证理论结果。我们通过验证相关双线性形式的强制性和连续性建立了离散格式的适定性，其中罚参数σ取足够大的值。在误差估计方面，我们的方法在能量范数下达到了最优收敛率O(h^(k?1))。更重要的是，当问题进入边界层主导的域时，最低阶版本（k=2）保持了不依赖于ε的鲁棒半阶精度O(h^(1/2))。这些分析结论得到了非规则多边形剖分上计算所得数值证据的支持，表明理论收敛阶与观测收敛阶之间存在良好的一致性。

本文其余部分安排如下。预备知识和网格假设收集在第2节。在第3节中，我们详细构造Morley型虚拟单元空间：通过约束消除顶点梯度和高阶内部矩，并建立自由度的单值性和可计算性。第4节引入IP格式——配备一个仅作用于函数跳跃的罚项——并给出适定性证明。第5节关注收敛分析，其中我们推导最优误差界以及边界层情形下的ε一致估计。第6节以多边形网格上的数值实验作为总结，旨在支持理论误差预测。

章节摘要
预备知识
在本节中，S?Ω表示任意开集。符号(·,·)S和‖·‖S表示通常的L2内积和诱导的范数，理解上这些适用于定义在S上的标量值函数和向量值函数的分量。对于整数m>0，我们使用标准的Sobolev空间H^m(S)和H?^m(S)，其范数和半范数分别记为‖·‖_{m,S}和|·|_{m,S}。当域S与整个Ω重合时，下标省略以便简化。

虚拟单元
在本节中，我们详细构造Morley型虚拟单元空间。我们的论述涵盖三个关键方面：所采用的自由度的单值性、插值误差界的推导以及相应基函数的有界性。

VEM方案
在本节中，我们构造内部罚函数方案并建立所得离散问题的适定性。特别关注仅作用于单元间边上函数值的跳跃罚项的设计。

误差分析
在本节中，我们建立所提方法的最优收敛性质，并推导在边界层主导域中成立的ε一致误差界。

数值结果
现在我们给出数值实验，以验证前面各节建立的理论收敛结果。计算在一系列非结构多边形剖分上进行；该系列中一个典型的网格如图1所示。这些多边形网格是使用文献[52]中记录的PolyMesher算法生成的。

在所有测试中，我们考虑所提方法的最低阶版本（即 k=2）应用于单位正方形域 Ω=(0,1)2 上的问题 (1.1)。

CRediT 作者贡献声明

赵济坤：写作 – 审阅与编辑，监督，软件，资源，概念化。党锋：验证，方法论，调查研究，形式化分析。刘琦：监督。杨永勤：监督。

竞争利益声明

作者声明他们没有已知的竞争性经济利益或个人关系，这些利益或关系可能被认为影响了本文所报道的工作。

赵济坤 | 党锋 | 刘琦 | 杨永勤

郑州大学数学与统计学院，中国郑州，450001