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

把低马赫（low Mach）流动放进可压缩程序，声速越快，画面似乎该越清楚。实际相反。马赫数越小，标准迎风（upwind）越是把一个完好的波从网格上抹掉。Arun、Das Gupta 与 Samantaray（2019）用线性波动方程组把这次崩溃拆开，再用 IMEX-RK 搭了一个与马赫数无关、仍然稳定的渐近保持（AP）格式。声学算子为什么必须隐式求解，下面的数字写得很直接。

把 $\varepsilon$ 滑块向左推，红色迎风曲线以斜率 $-1$ 爬升，青色 AP 曲线保持水平。

在 $\varepsilon=10^{-3}$ 处，两条曲线的比拉开到数百倍。马赫数下降时变大的不是波，是迎风数值粘性。

当 ε 趋于零，网格吃掉波

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

$\rho$ 是密度，$\mathbf{u}$ 是速度，$p=\rho^\gamma$ 是压力。$\varepsilon$ 是参考速度除以参考声速。

$\varepsilon \to 0$ 时，解走向不可压极限。这是奇异摄动（singular perturbation）。作者把分析收短，用的是线性波动组。

$\varrho$ 是线性化密度扰动，$a$ 是该状态的声速。有效声速是 $c=a/\varepsilon$。马赫数变小，这个声学速度就爆开。

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

极限 $\varepsilon\to 0$ 下，解并不随便漂。它收进已准备好（well-prepared）的空间：密度为常数，速度无散。

这个空间是声学算子 $L(U)=(a\nabla\cdot\mathbf{u},\,a\nabla\varrho)$ 的核（kernel）。离散解也得落到同一处。Dellacherie（2010）指出，打破这个不变性的格式会在低马赫下丢掉精度。保住核，就是渐近精确（AA）。稳定限制与 $\varepsilon$ 无关，就是渐近保持（AP）。

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

问题出在迎风的数值粘性。波动组对角化成特征变量后，变成速度为 $\pm c$ 的两个标量对流。每个模态套上 Rusanov 迎风，就有一步的放大因子。

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

一步的损失很小。堆起来才是问题。走到固定物理时间 $T$，需要 $M = T/\Delta t = T c /(\nu\Delta x)$ 个声学步。$\theta$ 小时，一步衰减是 $\tfrac{1}{2}\nu(1-\nu)\theta^2$，$M$ 次相乘得到这样的对数衰减。

马赫数降一个数量级，有效粘性也升一个数量级。上面那张图里红线斜率 $-1$，就是这个比例。

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

按时间尺度把方程拆开，算子有两个。

$H$ 是对流，时间尺度 $O(1)$。$L/\varepsilon$ 是声学，时间尺度 $O(\varepsilon)$。刚性全部在 $L/\varepsilon$ 上。只把这一部分隐式求解。对流保持显式。这种混合就是 IMEX（Implicit-Explicit）。同一分裂拿到马赫数 0 时，时间步还站不站得住，接在 IMEX TVD 的复现 里。

声学改为隐式之后，中心差分也稳定。迎风粘性退出。$\Delta t$ 由慢的对流决定，不由快的声学决定。在这个模型里，$\varepsilon$ 变小并不是缩小步长的理由。

IMEX-RK 接住刚性#

时间积分用 IMEX 龙格-库塔（IMEX-RK）。显式表 $(\tilde A,\tilde b)$ 负责 $H$，隐式表 $(A,b)$ 负责 $L/\varepsilon$。每一级里，声学凝成关于密度的椭圆方程，样子接近压力泊松。

$\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

把 Mach $\varepsilon$ 往下拨，红色振幅会在物理时间 $T=1$ 结束前掉到地板上。把 CFL $\nu$ 推向 1，$\nu(1-\nu)$ 变小，衰减会缓一口气。再把 $\varepsilon$ 降低，这点好处就被 $1/\varepsilon$ 的累积盖住。加大 $N$ 也只抵消一部分，因为 $\theta^2\propto N^{-2}$ 与 $\Delta x\propto N^{-1}$ 并不完全对消。

把代码里的数代进公式，步数有封闭形式。$a=1$、$T=1$、$\nu=0.5$、$N=64$ 时，$\Delta x=1/64$，$c=1/\varepsilon$。于是 $\Delta t=\nu\Delta x/c$，所以 $M=T/\Delta t=128/\varepsilon$。表里的 320、640、1280、2560、6400、12800，就是 $\varepsilon=0.4,0.2,0.1,0.05,0.02,0.01$ 代进这个式子的结果。$\varepsilon$ 减半，$M$ 翻倍。$\varepsilon$ 变成十分之一，$M$ 变成十倍。

对数衰减 $\tfrac12 T c(1-\nu)k^2\Delta x$ 代进同一组数，指数是 $0.154/\varepsilon$。系数来自 $k=2\pi$。$\varepsilon=0.4$ 时指数为 0.386，$\mathrm{e}^{-0.386}=0.680$，与输出的 $6.804\times 10^{-1}$ 在三位小数上一致。$\varepsilon=0.1$ 时 $\mathrm{e}^{-1.542}=0.214$，对上输出的 $2.143\times 10^{-1}$。$\varepsilon=0.01$ 时指数是 15.4。$\mathrm{e}^{-15.4}$ 在 $2.0\times 10^{-7}$ 附近，输出则是 $1.980\times 10^{-7}$。$\theta=k\Delta x=2\pi/64\approx 0.098$，小角度假设到表的末行并不完全准。表的底端比近似死得稍多，原因在这个 $\theta$。$\nu=0.5$ 时 $\nu(1-\nu)=0.25$，是这个乘积的最大值。所以这一列记的是单步衰减最强的那个 CFL。

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

两件事不是同一个问题。AP 问的是离散格式的 $\varepsilon\to 0$ 极限是不是极限方程的正确离散，稳定限制必须与 $\varepsilon$ 无关。AA 再问窄一点：网格有限、$\varepsilon$ 小但不为零时，解是否留在 well-prepared 空间附近。迎风两者都失败。IMEX 隐式两者都过。分叉点是时间离散，不是空间差分的阶。

即使 ν=0.5 是最狠的一步，AP 列仍是 1.00#

这张表的 AP 列六行都是 1.00。隐式步拿掉了迎风粘性，指数 $0.154/\varepsilon$ 不会再乘到振幅上。$\nu=0.5$ 已经是单步衰减因子最大的 CFL，堆完 $M=128/\varepsilon$ 步之后 AP 仍是 1.00。对这个线性波动组，只把 $L/\varepsilon$ 做成隐式，是对的选择。

这个选择不会原样搬到每一个低马赫系统。有的方程已经把声学做成隐式，时间步还是再次大幅缩短。相对速度里藏着第二个 $1/M^2$ 的位置，写在焓里剩下的刚性里。