一种全局无散度的虚拟元素方法，结合IMEX-SAV格式，用于Kelvin-Voigt粘弹性流体流动模型

《Communications in Nonlinear Science and Numerical Simulation》：A globally divergence-free virtual element method with IMEX-SAV scheme for the Kelvin-Voigt viscoelastic fluid flow model

　　白艳红 | 冯敏夫 | 高婷怡 | 魏华艺
西安华科大学理学院，成都，610039，中国

摘要
我们提出了一种稳定的、无散度的虚拟元素方法，用于一般多边形网格上的Kelvin-Voigt粘弹性流体流动模型。速度和压力分别通过类Stokes虚拟元素空间和任意阶的不连续分段多项式空间来近似。对于时间离散化，采用了IMEX-SAV方案来有效处理非线性项。通过结合适当设计的稳定策略，所提出的方法能够有效抑制在对流主导的流动状态下通常会出现的虚假振荡。在合理的假设下，我们建立了一个L∞→L2的逆不等式，并推导出了任意空间维度下类Stokes虚拟元素空间的插值误差估计。此外，获得了与模型参数无关的统一误差界限，强调了该方法在不同参数范围内的稳健性。最后，通过数值实验展示了所提出的VEM方法的准确性、稳健性和参数不敏感性。

引言
Kelvin-Voigt粘弹性流体流动模型描述了弱浓度水-聚合物溶液的运动。该模型最初由Pavlovskii在[1]中引入，随后由Oskolkov在[2]中命名。对于受周期性或无滑移Dirichlet边界条件约束的三维Kelvin-Voigt模型，已经在[3]中建立了正则化结果。我们现在回顾Kelvin-Voigt模型：
\[u_t - \kappa \Delta u_t - \nu \Delta u + (u \cdot \nabla)u + \nabla p = f\]
在\([0,T]\times\Omega\)上，\(\nabla \cdot u = 0\)，\(u = g\)在\([0,T]\times\partial\Omega\)上，\(u(0,x) = u_0(x)\)在\(\Omega\)内。
这里，\(\Omega\subset R^d\)（\(d=2,3\)）是一个有界的多边形域，\(\partial\Omega\)是Lipschitz边界，\([0,T]\)表示时间间隔。未知数\(u\)和\(p\)分别表示速度场和压力场；\(f\)是给定的体力，\(u_0\)表示规定的初始速度。参数\(\nu=1/Re>0\)是运动粘度，\(Re\)是雷诺数，而\(\kappa \geq 0\)表示表征变形弹性松弛的延迟时间。为了简便起见，我们仅考虑均匀的Dirichlet边界条件，即\(g = 0\)。

从数学角度来看，Kelvin-Voigt模型通过引入一个由\(\kappa\)加权的延迟项，扩展了不可压缩Navier-Stokes方程，从而考虑了材料响应中的记忆效应。在\(\kappa=0\)的极限情况下，系统简化为经典的Navier-Stokes方程。由于对流项的强非线性和由于不可压缩性约束\(\nabla \cdot u = 0\)而产生的速度\(u\)和压力\(p\)之间的耦合，解析解通常不可得。因此，开发稳健和精确的Kelvin-Voigt数值方法引起了广泛关注。在三角形和四边形网格上，代表性的方法包括连续和不连续Galerkin方法[4]以及双网格方法[5]。

近年来，由于计算域的几何复杂性以及对网格灵活性和适应性的需求增加，基于多边形网格离散化的数值方法受到了广泛关注。代表性的方法包括模仿有限差分（MFD）方法[6][7]、不连续Galerkin方法[8][9]、混合不连续Galerkin方法[10][11][12]、弱Galerkin方法[13][14][16]、混合高阶方法[17][18]和虚拟元素方法（VEM）[19][20][21][22][23][24]。在这些技术中，VEM作为MFD方法的推广[25]特别吸引人，因为它能够在Galerkin框架内处理一般的多边形和多面体网格。该方法避免了显式构造基函数，而是依赖于精心设计的自由度和可计算的投影算子，从而在数学分析和工程模拟中获得了广泛的应用。

对于与Kelvin-Voigt模型密切相关的不可压缩流动问题，包括Stokes方程和Navier-Stokes方程，已经提出了各种VEM公式并进行了分析；例如，参见[19][20][21][22][23][24][26][27]。为了确保质量守恒，已经提出了多种无散度VEM公式，包括H1-符合空间的[28][29][30]、非符合变体[31][32]以及弱VEM公式[33][34]。特别是，黄和王[30]将[28][29]中引入的无散度虚拟元素空间扩展到了任意空间维度。与标准公式相比，无散度VEM能够精确执行不可压缩性约束，产生压力稳健的速度近似，保持重要的物理结构，并且通常不需要额外的稳定化。为了减轻对流主导的不利影响，已经将多种稳定技术结合到VEM框架中，例如流线前向/Petrov-Galerkin（SUPG）方法[35][36]、基于涡度的分稳技术[37]、最小二乘稳定[38]和基于投影的方法[39][40]。特别是，[40]中提出的基于投影的稳定技术，结合类Stokes虚拟元素空间[28]，已经应用于二维非稳态Navier-Stokes方程。然而，由此产生的Crank-Nicolson时间离散化导致了一个需要在每个时间步长进行牛顿迭代的非线性方案，从而增加了计算成本，并且没有完整的理论分析。

隐式-显式（IMEX）方案通过隐式处理刚性项和显式处理非刚性项，结合了隐式和显式时间离散化的优势，从而提高了稳定性并降低了计算成本。因此，IMEX方案允许相对较大的时间步长，为时变问题提供了一种高效且稳定的方法（例如参见[41][42][43]）。标量辅助变量（SAV）方法是一种相对较新的数值策略，用于高效和稳定地处理梯度流和耗散系统[44]。其核心思想是通过引入一个标量辅助变量来重新制定原始能量函数，从而得到一个满足修改后的离散能量律的等效系统。SAV方法提供了两个主要优势：无条件能量稳定性和降低的计算复杂性。在非线性流体问题的背景下，已经报告了几个成功的应用和分析；例如，参见[16][45][46][47][48][49][50]。在[16]中，提出了一种结合IMEX-SAV方案的全局无散度弱Galerkin方法，用于二维Kelvin-Voigt模型。然而，这种方法的L∞(L2)范数误差估计仅为O(hk+τ)，并且数值实验仅限于三角形网格。

据我们所知，将虚拟元素方法应用于对流主导的Kelvin-Voigt模型尚未被研究。在这项工作中，我们为Kelvin-Voigt问题（1.1）开发了一种稳定的、无散度的虚拟元素方法。主要贡献总结如下：
(1) 任意阶离散化。我们考虑了二维和三维空间情况。速度由k阶的类Stokes符合虚拟元素空间[29][30]近似，而压力由k-1阶的不连续分段多项式空间[30]近似。结合了基于投影的稳定项以有效处理对流主导的情况。对于时间离散化，我们采用了IMEX方案与SAV方法相结合，从而在每个时间步长产生具有恒定系数的线性系统，从而提高了计算效率。
(2) 参数稳健的稳定化。所提出的稳定化在合理的稳定参数范围内对选择是稳健的，通过仅调整VEM公式中内在稳定项的系数，有效减轻了对流主导效应，而不会增加额外的计算成本。数值实验表明，即使对于依赖于\(\nu\)的解，该方法在对流主导的情况下以及在一般多边形网格上仍然稳健，稳定参数\(rK\)只需简单选择为O(1)。
(3) 插值和逆不等式。为了推导速度的均匀L∞(L2)误差估计，一个关键贡献是建立了任意空间维度下类Stokes虚拟元素空间的L∞→L2逆不等式[30]，以及插值算子IkK的稳定性和误差估计，包括L∞范数的界限。
(4) 统一误差估计。在这些理论结果的基础上，通过递归分析，我们推导出了速度的均匀L∞(L2)误差估计，其常数与粘度\(\nu\)和延迟时间\(\k\)的负幂无关。此外，在形式为\(\tau \leq Ch^2\)的二次CFL条件下，我们建立了速度的无条件能量估计，与[48]中报告的结果一致。

本文的其余部分组织如下。第2节介绍了模型问题、其变分公式和几个初步结果。第3节介绍了具有IMEX-SAV时间离散化的无散度VEM，并建立了完全离散方案的稳定性。第4节进行了完全离散方案的误差分析。最后，第5节报告了验证理论发现的数值实验。

在整篇文章中，符号\(a \sim b\)（分别是\(a \geq b\)）表示存在一个通用常数\(C > 0\)，使得\(a \leq Cb\)（分别是\(a \geq Cb\)）。常数\(C\)可能取决于精确解\(u, p\)和力项\(f\)的规则性，但与\(1/h, 1/\tau, 1/\nu, 1/\k\)无关。我们写\(a \sim b\)表示\(a \leq b \leq a\)。C的值可能因行而异。此外，\(C(\cdot)\)表示取决于括号中指示的量的常数。

**变分问题**
设\(M\)是所有\(d \times d\)张量的空间。对于任意有界域\(E\subset R^d\)，我们用\(W_{m,p}(E)\)表示标准Sobolev空间，其中函数的弱导数直到\(m\)阶属于\(L^p(E)\)，\(p \in [1, \infty)\)且\(m\)是非负整数。此外，我们用\(W_{m,p}(E;X)\)表示在有限维向量或张量空间\(X\)中取值的函数的相应Sobolev空间，其中\(X=R^d\)或\(M\)。当\(p=2\)时，我们写\(W_{m,p}(E)=H_m(E)\)和\(W_{m,p}(E;X)=H_m(E;X)\)。相关的范数和...

**稳定化半离散方案**
为了简化空间离散化的分析复杂性，我们将讨论限制在Kelvin-Voigt模型的二维或三维虚拟元素空间上。遵循标准的VEM构造[29]，我们引入了局部双线性和三线性形式。对于\(u,w,v \in H_1(K;R^d)\)和\(p, q \in L^2(K)\)，它们定义如下：

**半离散方案的误差估计**
定理4.1
设\((J, u, p) \in L^2([0, T]) \times L^2([0, T]; V) \times L^2([0, T]; Q)\)和\((J_h, uh, q_h) \in L^2([0, T]) \times L^2([0, T]; V_h) \times L^2([0, T]; Q_h)\)分别是变分问题(2.3)和半离散方案(3.34)的解。以下误差估计成立：
\[ \max_{0 \leq t \leq T} \{ \sum_{K \in Th} \| uh - u \|_0 \|_K^2 + hK^2 shK((I - \Pi_kK)(uh - u), (I - \Pi_kK)(uh - u)) + \kappa \| \nabla(uh - u) \|_0 \|_K^2 + |J_h - J|^2 + \int_{0}^T (\nu \| \nabla(uh - u) \|_0^2 + sh_1(uh - u, uh - u)) ds \leq C(u, ut, f, T)h^2 \], \int_{0}^T \sum_{K \in Th} \| ph - p \|_0 \|_K ds \leq C(u, ut, f, T)h^k. \]

**备注4.2**
为了节省空间，证明被省略了；可以通过...

**数值实验**
我们进行了一系列数值实验来验证二维设置（\(d=2\)）下完全离散方案(3.37)的理论结果。特别强调了验证预测的收敛率、评估参数稳健稳定性的有效性，以及证明误差常数不依赖于参数\(\nu\)和\(\k\)。

**结论**
在这项工作中，我们为一般多边形网格上的非稳态Kelvin-Voigt问题开发了一种稳定且全局无散度的虚拟元素方法。通过结合类Stokes虚拟元素离散化和适当的稳定策略，所提出的方法能够有效处理对流主导的情况。对于时间离散化，采用了IMEX-SAV方案，从而获得了对模型参数稳健的速度误差估计。

**声明**
在准备这项工作期间，作者使用ChatGPT来纠正语法错误并改进手稿的英语语言。使用此工具/服务后，作者根据需要审阅和编辑了内容，并对已发表文章的内容承担全部责任。

**作者贡献声明**
白艳红：写作——审阅与编辑、撰写——原始草稿、验证、方法论。
冯敏夫：项目管理、方法论。
高婷怡：验证。
魏华艺：验证。

**利益冲突声明**
作者声明他们没有已知的竞争性财务利益或个人关系可能会影响本文报告的工作。