[论文综述] 低马赫下不杀死声音的方法 — AP-IMEX 波动格式

用可压缩程序求解低马赫（low Mach）流动，直觉上声速越快、解应越清晰。实际恰恰相反。马赫数越小，标准迎风（upwind）格式越是把一个完好的波从网格上抹掉。本文跟随 Arun、Das Gupta 与 Samantaray（2019）的工作，以线性波动方程组为模型诊断这种崩溃，并通过 IMEX-RK 时间积分构造一个与马赫数无关、稳定的渐近保持（AP）格式。读完之后，你能在代码里直接验证"为何必须隐式求解声学算子才能在低马赫下保持精度"。

当 ε 趋于零，网格吃掉波#

经过尺度化的等熵 Euler 方程里嵌着一个小参数：马赫数 $\varepsilon$ 。压力项以 $\nabla p / \varepsilon^2$ 进入动量方程。

\partial_t \rho + \nabla\cdot(\rho \mathbf{u}) = 0, \qquad \partial_t (\rho\mathbf{u}) + \nabla\cdot(\rho\mathbf{u}\otimes\mathbf{u}) + \frac{\nabla p}{\varepsilon^2} = 0

其中 $\rho$ 为密度， $\mathbf{u}$ 为速度， $p=\rho^\gamma$ 为压力， $\varepsilon$ 为参考速度与参考声速之比（马赫数）。

当 $\varepsilon \to 0$ ，解收敛到不可压极限。在数学上这是奇异摄动（singular perturbation）问题。为简化分析，作者使用线性波动方程组。

\partial_t \varrho + (\mathbf{u}\cdot\nabla)\varrho + \frac{a}{\varepsilon}\nabla\cdot\mathbf{u} = 0, \qquad \partial_t \mathbf{u} + (\mathbf{u}\cdot\nabla)\mathbf{u} + \frac{a}{\varepsilon}\nabla\varrho = 0

其中 $\varrho$ 为线性化的密度扰动， $a$ 为线性化状态的声速。有效声速为 $c = a/\varepsilon$ 。马赫数越小，这个声学速度就越爆炸。

well-prepared 空间 — 必须存活的解#

在极限 $\varepsilon\to 0$ 下，解并非随意漂移。它收敛到"已准备好（well-prepared）"的空间：密度为常数、速度无散的场的集合。

\varrho^{(0)} = \text{const}, \qquad \nabla\cdot\mathbf{u}^{(0)} = 0

这个空间是声学算子 $L(U)=(a\nabla\cdot\mathbf{u},\,a\nabla\varrho)$ 的核（kernel）。好的格式应让其离散解收敛到同一空间。Dellacherie（2010）证明，破坏这一不变性的格式会在低马赫下丧失精度。保持核即渐近精确（AA）；稳定约束与 $\varepsilon$ 无关即渐近保持（AP）。

迎风粘性随声速成比例增大#

问题的根源是迎风格式的数值粘性。把波动组对角化为特征变量，它分裂为以速度 $\pm c$ 传播的两个标量对流。对每个模态施加迎风（Rusanov）格式，得到每步的放大因子。

|g|^2 = \big(1 - \nu(1-\cos\theta)\big)^2 + \nu^2\sin^2\theta, \qquad \nu = \frac{c\,\Delta t}{\Delta x},\ \ \theta = k\Delta x

其中 $\nu$ 为声学 CFL 数， $\theta$ 为每格相位， $k$ 为波数。

关键在累积。到达固定物理时间 $T$ 需要 $M = T/\Delta t = T c /(\nu\Delta x)$ 个声学步。在小 $\theta$ 下每步衰减为 $\tfrac{1}{2}\nu(1-\nu)\theta^2$ ，累加 $M$ 次后对数衰减如下。

M\cdot\tfrac{1}{2}\nu(1-\nu)\theta^2 \;=\; \tfrac{1}{2}\,T\,c\,(1-\nu)\,k^2\Delta x \;\propto\; c \;\propto\; \frac{1}{\varepsilon}

马赫数每降一个数量级，有效粘性就升一个数量级。在下面的图中自己拖动看看。

Mach ε = 0.050

D_explicit ≈ 0.078

D_AP ≈ 7.8e-3

ratio ≈ 10.0×

把 $\varepsilon$ 向左移，迎风曲线（红）以斜率 $-1$ 飙升，而 AP 曲线（青）保持水平。在 $\varepsilon=10^{-3}$ 处两者之比拉开到数百倍。这就是低马赫精度崩溃的定量面貌。

算子分裂 — 对流显式，声学隐式#

处方从按时间尺度分裂方程开始。改写为发展形式，便出现两个算子。

\partial_t U + H(U) + \frac{1}{\varepsilon}L(U) = 0

其中 $H$ 为时间尺度 $O(1)$ 的对流算子， $L/\varepsilon$ 为时间尺度 $O(\varepsilon)$ 的声学算子。所有刚性都来自 $L/\varepsilon$ 。那么只把这部分隐式求解即可。对流不贵，便保持显式。这种显-隐混合就是 IMEX（Implicit-Explicit）。

隐式处理声学后，即使中心差分也稳定。迎风粘性消失。时间步由慢的对流尺度决定，而非快的声学。 $\varepsilon$ 再小，也无需缩小 $\Delta t$ 。

IMEX-RK 接住刚性#

作者用 IMEX 龙格-库塔（IMEX-RK）做时间积分。显式表 $(\tilde A,\tilde b)$ 处理 $H$ ，隐式表 $(A,b)$ 处理 $L/\varepsilon$ 。每一级中，声学部分凝缩为关于密度的椭圆方程（形似压力泊松）。

\varrho - \beta^2 \Delta t^2\, a^2 \nabla^2 \varrho = (\text{显式项})

其中 $\beta$ 为 IMEX 表的对角系数， $\nabla^2$ 为拉普拉斯算子。该椭圆问题的可逆性由变分鞍点（saddle point）理论保证。在离散层面，循环矩阵（circulant matrix）理论给出同样结果。由此得到三点：唯一解的存在、与 $\varepsilon$ 无关的一致稳定，以及 well-prepared 空间的不变性。

Python — 复现单一模态的振幅衰减#

下面的代码计算单一傅里叶模态的迎风放大因子，并按马赫数输出固定物理时间 $T$ 之后残留的振幅。隐式格式保存模态，故其振幅保持为 1。

import numpy as np
 
a  = 1.0        # 线性化声速
k  = 2 * np.pi  # [0,1] 上的单一傅里叶模态
T  = 1.0        # 固定物理终止时刻
nu = 0.5        # 声学 CFL 数
 
def upwind_amp_factor(eps, N):
    """波动组 Rusanov 单一特征模态的放大 |g|。"""
    c  = a / eps              # 声学速度 c = a/ε
    dx = 1.0 / N
    th = k * dx               # 每格相位 θ = kΔx
    re = 1.0 - nu * (1.0 - np.cos(th))
    im = nu * np.sin(th)
    return np.hypot(re, im)   # |g| <= 1
 
def retained_amplitude(eps, N):
    """物理时间 T 后的模态振幅（迎风 vs 隐式）。"""
    c  = a / eps
    dx = 1.0 / N
    dt = nu * dx / c          # 声学 CFL 步长
    M  = T / dt               # 显式步数
    explicit = upwind_amp_factor(eps, N) ** M
    implicit = 1.0            # 中心-隐式保存模态
    return M, explicit, implicit
 
if __name__ == "__main__":
    N = 64
    print(f"{'eps':>8} {'steps M':>10} {'explicit':>12} {'AP':>6}")
    for eps in [0.4, 0.2, 0.1, 0.05, 0.02, 0.01]:
        M, ex, ap = retained_amplitude(eps, N)
        print(f"{eps:8.3f} {M:10.0f} {ex:12.3e} {ap:6.2f}")

输出如下。

     eps    steps M     explicit     AP
   0.400        320    6.804e-01   1.00
   0.200        640    4.629e-01   1.00
   0.100       1280    2.143e-01   1.00
   0.050       2560    4.591e-02   1.00
   0.020       6400    4.430e-04   1.00
   0.010      12800    1.980e-07   1.00

每当 $\varepsilon$ 缩小十倍，步数翻倍，迎风振幅指数式地趋向零。在 $\varepsilon=0.01$ 处模态实质上消失。在同一网格上隐式格式保住了 100%。

在下面的模拟里直接操作。把 Mach $\varepsilon$ 滑块下移，红色曲线（迎风）在物理时间还未流逝前就跌到底。

Mach ε: 0.050acoustic CFL ν: 0.50grid cells N: 64

把 CFL $\nu$ 推向 1， $\nu(1-\nu)$ 变小，衰减会一时缓和。但再降 $\varepsilon$ ，这点收益很快被 $1/\varepsilon$ 的累积淹没。加大网格 $N$ 只有部分帮助，因为 $\theta^2\propto N^{-2}$ 与 $\Delta x \propto N^{-1}$ 只能部分抵消，根本问题仍在。

渐近保持（AP）与渐近精确（AA）#

这两个性质容易混淆。AP 问的是：离散格式的 $\varepsilon\to 0$ 极限是否为极限方程的正确离散，且稳定约束与 $\varepsilon$ 无关。AA 更进一步：在有限网格上，对小而非零的 $\varepsilon$ ，解是否停留在 well-prepared 空间附近。迎风格式两者皆败。IMEX-隐式格式两者皆满足。决定性因素是时间离散的选择，而非空间差分的精度。

记住的三行#

低马赫下迎风粘性随声速 $c=a/\varepsilon$ 以 $1/\varepsilon$ 增大，在固定物理时间上累积，把波抹掉。
是否保持 well-prepared 空间（常密度、无散速度）是精度的分水岭，而左右它的是时间离散。
只把声学算子 $L/\varepsilon$ 用 IMEX 隐式处理， $\Delta t$ 便从声学尺度解放到对流尺度，使格式与马赫数无关地稳定且精确（AP＋AA）。