Arbitrary order stationarity preserving stabilized finite elements for multidimensional nonlinear hyperbolic problems. Application to the Euler equations with gravity

这是一份关于论文《Arbitrary order stationarity preserving stabilized finite elements for multidimensional nonlinear hyperbolic problems. Application to the Euler equations with gravity》（多维非线性双曲问题的任意阶保稳态稳定有限元方法：应用于带重力的欧拉方程）的详细技术总结。

1. 研究背景与问题 (Problem)

• 核心挑战：在计算流体力学（CFD）中，数值模拟经常需要处理接近稳态（平衡态）的流动，例如大气中的静力平衡、天体物理中的吸积盘或低马赫数流动。传统的数值格式（如标准的有限元或有限体积法）在处理这些状态时，往往会产生非物理的数值误差（spurious numerical errors）。这些误差会破坏真实的物理平衡，导致虚假的波动，甚至掩盖真实的物理现象（如不稳定性或微弱的扰动）。
• 现有局限：
  • 现有的“保平衡”（Well-balanced）方法大多基于一维问题或特定的平衡态（如静止流体），难以推广到真正多维（genuinely multidimensional）的复杂平衡态（如旋转涡旋、非均匀重力场下的平衡）。
  • 对于高阶连续有限元方法（如 SUPG），其数值耗散结构通常与静力平衡条件（$\nabla \cdot F = S$）不兼容。特别是基于拉普拉斯算子的耗散项会破坏无散度条件（$\nabla \cdot v = 0$），导致无法在离散层面精确保持稳态。
  • 缺乏一种能够同时满足任意高阶精度、多维非线性平衡态保持以及低马赫数兼容性的统一框架。

2. 方法论 (Methodology)

本文提出了一种基于**全局通量积分（Global Flux Quadrature, GFQ）**的任意高阶保稳态稳定有限元方法（SUPG-GFQ）。

2.1 核心思想：通量势预处理

作者将 Barsukow 等人提出的“全局通量”概念推广到非线性平衡律和连续有限元框架中。

• 通量重构：不再直接离散散度算子 $\nabla \cdot F$，而是引入通量势（Flux Potentials）$F_1, F_2$ 和源项势 $S$。
  • 定义 $F_1 = \int F_1 d\eta$, $F_2 = \int F_2 d\xi$, $S = -\iint S d\xi d\eta$。
  • 利用恒等式 $\nabla \cdot F - S = \partial_{x_1}\partial_{x_2}(F_1 + F_2 + S)$。
• 离散实现：
  • 在笛卡尔网格上，利用 LobattoIIIA 配点法（与 Gauss-Lobatto 节点对应）对通量和源项进行逐行/逐列积分。
  • 这种预处理将连续的多项式向量场映射到一个具有 Raviart-Thomas (RT) 型结构的局部空间。
  • 这使得离散后的散度算子天然地满足保散度（solenoidal）性质，从而能够精确保持稳态。

2.2 数值格式构建

• 空间离散：采用连续有限元（Continuous FEM），使用张量积形式的 $Q_K$ 多项式空间。
• 稳定化：结合流线迎风 Petrov-Galerkin (SUPG) 方法。关键在于，稳定化项中的残差计算也基于上述预处理后的通量势，确保了稳定化项在稳态下自动消失（即 $\nabla \cdot F - S = 0 \implies \text{Stabilization} = 0$）。
• 时间离散：采用 缺陷修正（Defect Correction, DeC） 方法。这是一种显式迭代策略，用于求解由 SUPG 引入的非对角质量矩阵，从而获得任意高阶的时间精度，同时避免直接求逆矩阵的昂贵计算。
• 源项修正（Well-Balanced Correction）：针对带重力的欧拉方程，提出了一种特殊的源项离散化方法（基于参考状态 $\bar{\rho}, \bar{p}$），使得等温静力平衡态（Isothermal Hydrostatic Equilibria）可以被机器精度（Machine Precision）地保持。

2.3 理论分析

• 保稳态定义：论文形式化定义了“保稳态”（Stationarity Preserving, SP）离散格式，并提出了实现该性质的必要条件（对偶残差条件）。
• 超收敛性（Super-convergence）：证明了在稳态解附近，该方法的节点解具有比名义阶数更高的精度（对于 $K \ge 2$，达到 $K+2$ 阶或更高），这是由 LobattoIIIA 积分的内在性质决定的。
• 低马赫数兼容性：由于该方法在平衡态下退化为中心格式（无数值耗散），且数值耗散项具有正确的 $\nabla(\nabla \cdot v)$ 结构，因此天然具备低马赫数兼容性。

3. 主要贡献 (Key Contributions)

1. 理论框架的扩展：首次将保稳态的全局通量积分方法从线性声学/浅水方程扩展到非线性多维平衡律（如欧拉方程），并适配了任意高阶连续有限元。
2. 严格的保稳态定义与条件：给出了保稳态格式的严格数学定义，并提出了“对偶残差条件”作为判断格式是否保稳态的实用判据。
3. Raviart-Thomas 结构的隐式引入：揭示了通过通量势预处理，标准的 $Q_K$ 连续有限元空间在离散散度算子作用下，等价于具有 Raviart-Thomas 结构的离散空间，从而解决了连续元难以保持无散度条件的难题。
4. 机器精度的静力平衡保持：针对带重力欧拉方程，提出了一种简单的源项重构技术，实现了等温静力平衡态的机器精度保持，无需复杂的重构或分裂。
5. 低马赫数性能提升：证明了该方法在低马赫数极限下表现优异，能够避免标准 SUPG 格式在低马赫数下的数值耗散过大问题。

4. 数值结果 (Results)

论文通过广泛的基准测试验证了方法的有效性：

• 移动等熵涡旋（Moving Isentropic Vortex）：
  • 在非稳态情况下，SUPG-GFQ 的误差幅度比标准 SUPG 小 1.3 到 2 倍，且保持了预期的 $K+1$ 阶收敛精度。
• 稳态等熵涡旋（Steady Isentropic Vortex）：
  • 标准 SUPG 无法保持涡旋结构，产生虚假结构；而 SUPG-GFQ 即使在粗网格上也能完美保持涡旋。
  • 超收敛现象：对于 $K=2, 3$，SUPG-GFQ 展现出 $K+2$ 甚至更高的收敛阶，验证了理论预测。
• 低马赫数行为（Low Mach Behavior）：
  • 在 $Ma = 10^{-2}$ 到 $10^{-6}$ 的范围内，SUPG-GFQ 的误差远小于标准 SUPG 和传统的 HLLC 有限体积法。在 $Ma=10^{-6}$ 时，粗网格即可达到机器精度。
  • 长时间模拟（$t=2000s$）显示，SUPG-GFQ 能完美保持初始结构，而标准方法会因数值耗散导致涡旋消散。
• 流体不稳定性（Kelvin-Helmholtz & Rayleigh-Taylor）：
  • 在 Kelvin-Helmholtz 不稳定性测试中，SUPG-GFQ 在粗网格上就能捕捉到清晰的涡旋结构，且网格印记（mesh imprinting）极少。
  • 在 Rayleigh-Taylor 不稳定性中，该方法保持了径向对称性，而标准方法表现出明显的网格依赖性。
• 热气泡上升（Rising Thermal Bubble）：
  • 模拟大气中的热气泡上升，SUPG-GFQ 生成的蘑菇云形状更规则，无虚假振荡，显著优于标准 SUPG。
• 激波 - 涡旋相互作用：
  • 在 $M=1.1$ 的激波测试中，方法保持了稳定性，且高阶格式进一步减少了网格印记。

5. 意义与影响 (Significance)

• 通用性与高精度：该方法提供了一种通用的框架，能够以任意高阶精度处理多维非线性双曲守恒律，特别适用于需要极高精度的稳态或近稳态模拟。
• 物理保真度：通过消除非物理的数值耗散对平衡态的破坏，该方法能够更真实地模拟大气物理、天体物理中的微弱扰动和不稳定性演化。
• 低马赫数模拟的突破：为低马赫数流动（如燃烧、气象）提供了一种无需特殊低马赫预处理（Preconditioning）即可保持精度的离散方案。
• 计算效率：虽然引入了预处理步骤，但通过 DeC 方法实现了显式时间推进，避免了隐式求解的昂贵代价，使得高阶方法在实际应用中更具可行性。

总结：这篇论文通过引入全局通量积分和通量势预处理，成功地将保稳态能力赋予了任意高阶连续有限元方法。它不仅解决了长期存在的“稳态保持”与“高阶精度”之间的矛盾，还显著提升了低马赫数流动和复杂不稳定性问题的模拟能力，是计算流体力学数值格式发展的重要进展。