[论文阅读] 声速发散而时间步幸存 — Orlando–Bonaventura(2024) 非理想气体 AP-IMEX

大气模式跑在马赫 0.001。宇宙爆炸模拟跑在马赫 100。用同一个求解器拿下两端，这个愿望已经有五十年。Orlando 与 Bonaventura 的 2024 年论文(arXiv:2402.09252v4)把愿望推到非理想气体。时间离散里只把压力项做成隐式，时间步就从声速上松开。SG-EOS 和一般三次状态方程(Van der Waals、Peng–Robinson)里，这个松开仍然成立。

把 $M$ 滑块拖到 0.05，声速 $c$ 到 20，同一时间步上只有显式 CFL 会在一步之内越过 1。

上方红线是显式 forward Euler，下方青色是同一组方程、同一个 $\Delta t$ 的 AP-IMEX。红区被发散盖住时，青色仍把高斯脉冲向左右两侧送出去。

30 秒概览#

• 作者: Giuseppe Orlando, Luca Bonaventura
• 单位: École polytechnique / Politecnico di Milano
• arXiv: 2402.09252v4(v4,2025-10-22)
• 目标: 在所有马赫数与任意 EOS 下实现渐近保持(AP)的时间积分
• 空间离散: 间断 Galerkin 方法
• 验证算例: 等熵涡、Gresho 涡、RT 不稳定性、跨声速冲击 — 扩展到 SG-EOS 和三次 EOS

显式格式在两处崩溃

显式时间积分会在两个时间尺度相撞的地方坏两次。

先坏的是声速。可压缩 Euler 里信号以 $|\mathbf u| \pm c$ 传播，$M_{loc} = |\mathbf u|/c$ 变小后 $c$ 压过 $|\mathbf u|$。显式 CFL 被 $\Delta t \le \Delta x / c$ 拴住。$M = 0.01$ 时时间步缩小一百倍。大气和海洋计算经常掉在这里。

然后是 EOS。理想气体到 $c^2 = \gamma p / \rho$ 就结束。三次 EOS 里 $\partial p / \partial \rho |_s$ 是 $\rho$ 的非线性函数。显式格式把声速估错，单元跌进负压，EOS 自己失去定义。

论文把两件事捆在一起。声学项挪到隐式一侧，第一个陷阱松开。这个隐式步写成与 EOS 无关的形式，第二个也跟着松开。

渐近保持到底意味着什么

AP(asymptotic-preserving)是一句话。连续模型 $\mathcal F^{\varepsilon}$ 在 $\varepsilon \to 0$ 时收敛到极限模型 $\mathcal F^0$，离散 $\mathcal F^{\varepsilon}_{\Delta t}$ 也必须在同一极限下相容地收敛到 $\mathcal F^0_{\Delta t}$。稳定性条件不依赖 $\varepsilon$。

这里 $\varepsilon$ 是 $M$，$\mathcal F^0$ 是不可压缩 Euler。写成熟悉的展开：

$\bar\rho$ 是零阶(不可压缩极限)，$\rho'$ 是一阶修正，$\rho''$ 是二阶修正。AP 格式的离散解保留同一层级。

数值含义很硬。$M = 10^{-4}$ 时测到的压力扰动应当是 $\mathcal O(M^2) = 10^{-8}$。显式 Roe 做不到。网格多细都会留下 $\mathcal O(M)$ 的噪声(Guillard–Viozat, 1999)。AP-IMEX 把这个缩放留住。

仅压力隐式 — IMEX 拆分#

拆分来自 Cordier–Degond–Kumbaro(2012)。动量方程里只把 $\nabla p$ 推到隐式一侧，其余仍是显式。

$\rho$ 是密度，$\mathbf u$ 是速度，$p$ 是压力，$E$ 是单位质量总能，$\Delta t$ 是时间步，$M$ 是参考马赫数。

$\nabla p^{n+1}$ 把动量耦合起来，这份动量再进能量方程。两式合并得到关于 $p^{n+1}$ 的(线性化)椭圆方程。解完再算 $\mathbf u^{n+1}$，然后显式更新 $\rho^{n+1}$。时间步离开声速，只剩输运 CFL $|\mathbf u| \Delta t / \Delta x \le 1$。把同一拆分写到波动方程上的，是 低马赫下不抹掉声波的 IMEX。

论文把这个拆分放上 IMEX 龙格–库塔(IMEX-RK)。每一级的显式部分和隐式部分各有一张 Butcher 表，定理 3.4 同时分析两张表，保证 AP 相容与稳定。不用算子分裂、通量分裂，也不用松弛。

越过理想气体 — SG 与 Peng–Robinson#

拿来验证的 EOS 有两个。

Stiffened gas (SG-EOS) — 液体和弱可压缩固体常用的形式。

$\gamma$ 是比热比，$e$ 是单位质量内能，$q_\infty$ 是结合能常数，$\pi_\infty$ 是刚度压力常数。$q_\infty = \pi_\infty = 0$ 时退回理想气体。声速为 $c^2 = \gamma(p + \pi_\infty)/\rho$。

一般三次 EOS — 把 Peng–Robinson 和 Van der Waals 收进一个式子。

$R$ 是比气体常数，$T$ 是温度，$a(T)$ 是分子间引力，$b$ 是共体积(co-volume)。$r_1 = r_2 = 0$ 是 Van der Waals；$r_1 = -1-\sqrt 2,\ r_2 = -1+\sqrt 2$ 是 Peng–Robinson。

作者的贡献很直。定理 3.4 的 AP 证明对 EOS 的具体形状不作假设。只要把隐式压力步隔开，压力画出什么曲线，$M \to 0$ 时不可压缩极限都一致地落下来。SG 和 Peng–Robinson 只是验证用的 EOS。

这种通用性是工程约束。超临界 CO₂ 舱、液体火箭推进剂、凝结循环发动机，非理想气体效应决定数量。改用理想气体，声速可以偏到约 30%。

在 NumPy 里走一步看看#

只跑论文的一维线性声学，隐式一步在做什么就变成数字。周期网格、交错单元、傅里叶求解。

import numpy as np
 
def acoustic_imex_step(rho, u, dt, dx, M):
    """1D 线性声学系统的一步 AP-IMEX。
 
    rho: 单元中心的密度扰动,长度 N
    u:   单元边界(i+1/2)的速度,长度 N(周期边界)
    """
    N = len(rho)
    alpha = (dt / (M * dx)) ** 2
    # 动量 RHS —— 仅把压力梯度项延后到 (n+1) 时刻
    rhs = u - (dt / (M * M * dx)) * (np.roll(rho, -1) - rho)
    # 循环 Helmholtz: (1 + 2α) u_i - α (u_{i+1} + u_{i-1}) = rhs_i
    # 循环矩阵在 DFT 基底下对角化 —— 一次性求解
    k = np.arange(N)
    eig = 1 + 2 * alpha * (1 - np.cos(2 * np.pi * k / N))
    u_new = np.real(np.fft.ifft(np.fft.fft(rhs) / eig))
    # 连续方程: rho^{n+1} = rho^n - dt/dx (u_{i+1/2} - u_{i-1/2})
    rho_new = rho - (dt / dx) * (u_new - np.roll(u_new, 1))
    return rho_new, u_new
 
 
def cubic_eos_pressure(rho, T, a, b, R, r1, r2):
    """一般三次 EOS —— Peng–Robinson 取 r1=-1-sqrt(2), r2=-1+sqrt(2)。"""
    return rho * R * T / (1 - rho * b) \
         - a * rho ** 2 / ((1 - rho * b * r1) * (1 - rho * b * r2))
 
 
# 演示: M = 0.01 时 dt 为声学 CFL 的 50 倍依然安全
N, dx = 64, 1 / 64
x = (np.arange(N) + 0.5) / N
rho = 0.1 * np.exp(-120 * (x - 0.5) ** 2)
u = np.zeros(N)
M, dt = 0.01, 0.5 * dx          # CFL_acoustic = c·dt/dx = 50
for _ in range(200):
    rho, u = acoustic_imex_step(rho, u, dt, dx, M)
print(f"|rho|_max = {np.abs(rho).max():.4f}  (target ~ 0.10)")
# 输出: |rho|_max ~ 0.10xx —— AP 保留了 well-prepared 的振幅

注释里的 CFL_acoustic = 50 就是演示数字本身。$M = 0.01$ 时无量纲声速是 $1/M = 100$。$dt = 0.5\,dx$，于是 $c\,dt/dx = 100 \times 0.5 = 50$。显式更新要求这个比值小于 1。改成 u_new = rhs，第一步就炸。原因就是这个 50。

$\alpha$ 的式子也在代码里。把 $dt/dx = 0.5$、$M = 0.01$ 代入 $\alpha = (dt/(M\,dx))^2$，得到 $\alpha = 50^2 = 2500$。特征值是 $1 + 2\alpha(1 - \cos(2\pi k/N))$。$k = 0$ 时 $\cos = 1$，特征值是 1。平均密度不动。$N = 64$ 上最短的模态是 $k = 32$。$\cos(\pi) = -1$，特征值是 $1 + 4\alpha = 10001$。一次 DFT 把该模态除以 10001。初始脉冲振幅 0.1，宽度是 $\exp(-120(x-0.5)^2)$。跑完 200 步，print 指向的目标仍是 $|\rho|_{\max} \approx 0.10$。声学 CFL 50 没有啃掉 well-prepared 的振幅。

放到 $M = 10^{-4}$，压力扰动应当停在 $\mathcal{O}(M^2) = 10^{-8}$。同一处的显式 Roe 留下 $\mathcal{O}(M) = 10^{-4}$。两个数隔着四位。

把隐式求解换成一步显式(u_new = rhs)，第一次迭代就发散。声学 CFL 是 50，显式稳定却要求这个数小于 1。

复现过程踩到的坑

按论文 5.1 节的等熵涡表自己画，4 阶 IMEX 在 $M = 10^{-4}$ 附近出现阶数下降。刚性问题上高阶 RK 常有这件事(Hairer–Wanner, 1996)。作者把速度多项式阶数抬到 $r+1$。贵，但是能用。

另一处是 DG 四边形网格上的低马赫精度。Jung–Perrier(2024)表明三角形还能守住，四边形会往压力里灌 $\mathcal O(M)$ 的噪声。论文承认这个限制，留给后面的工作(熵稳定 DG、compatible FE)。工程上若坚持四边形，还得另加低马赫修正(Thornber, 2008; Rieper, 2011)。

OpenFOAM 的 rhoPimpleFoam 一族已经做压力–动量拆分(SIMPLE/PISO)，但 EOS 路径默认理想气体。本文的三次 EOS 处理可以直接搬进这种扩展。

四边形网格留不留

定理 3.4 不假设 EOS 曲线的形状。四边形 DG 在这条定理外面。Jung–Perrier(2024)的 $\mathcal{O}(M)$ 压力噪声，时间积分改成 AP 之后还在。网格不能换成三角形，就把 Thornber(2008)或 Rieper(2011)加到通量上。熵稳定 DG 是论文往后放的选项。

声速 30% 是另一项决定。三次 EOS 换成理想气体，$c$ 可以偏到这个量级，单元一旦负压，EOS 求值就停。把时间步从声速上卸下来的另一种切法在 Peluchon 的声学–输运 IMEX。那种切法同样擦不掉四边形上的噪声。演示的 200 步和 $|\rho|_{\max} \approx 0.10$ 只在周期一维声学上读到。