[论文评述] 不让声速吞掉你的 CFL — Peluchon(2017) 声-传输 IMEX 分裂

一位跑大气再入仿真的博士生叹了口气。液-气界面处声速 1500 m/s，流速只有 5 m/s。CFL 完全由声速决定，每步只有微秒级。八小时的实际计算还没走完 1 秒物理时间。Peluchon、Gallice 与 Mieussens（2017）直击这一痛点 —— 在五方程多相模型里把声波单独隐式化，物质波依旧显式。

论文信息#

作者：S. Peluchon, G. Gallice, L. Mieussens
机构：CEA-CESTA / Univ. Bordeaux / INRIA
期刊：Journal of Computational Physics, 339, 2017
DOI：10.1016/j.jcp.2017.03.019
标题：A robust implicit-explicit acoustic-transport splitting scheme for two-phase flows

一句话总结 — 把声速从时间步分母里剪下来#

五方程 diffuse interface 模型用同一组守恒方程同时求解液体和气体。但液体声速（1500 m/s）是气体（340 m/s）的四倍多。显式 Godunov 受 $\Delta t < \Delta x / (|u|+c)$ 约束。在液体区 $|u| \ll c$ ，所以步长基本上是 $\Delta x / c$ 。1 mm 网格意味着每步只有 0.67 μs。

论文的核心是把那个 $c$ 从时间步分母里挑出来。只对声学步采用 backward Euler 隐式，传输步保持显式迎风。外层 CFL 就只看 $|u|/c$ ，即按 Mach 数比例提高 1/M 倍。Mach 0.01 时理论加速可达 100 倍。

这一思路源自 Coquel–Godlewski–Khalfallah（2016）的单相气动力学方法。Peluchon 把它移植到了 Allaire–Clerc–Kokh（2002）的五方程多相系。

把五个方程切成声学与传输#

1D 可压缩五方程系写作

\partial_t \rho + \partial_x (\rho u) = 0

\partial_t (\rho y) + \partial_x (\rho y u) = 0

\partial_t (\rho u) + \partial_x (\rho u^2 + p) = 0

\partial_t (\rho e) + \partial_x ((\rho e + p) u) = 0

\partial_t z + u \partial_x z = 0

其中 $z$ 是体积分数， $y = \rho_1 z / \rho$ 是第一相质量分数， $e = \epsilon + u^2/2$ 是比总能。只有体积分数方程是非守恒的。

论文把每个守恒变量的演化拆为 两部分。对 $\phi \in \{\rho, \rho y, \rho u, \rho e\}$ ，

\partial_t \phi + \phi \partial_x u + u \partial_x \phi = \partial_t \phi + \partial_x (\phi u)

左边前两项构成 声学部分（ $\phi \partial_x u$ 表示压缩/膨胀），第三项是 传输部分（ $u \partial_x \phi$ 表示对流）。含声速的压力梯度 $\partial_x p$ 归入动量与能量的声学部分。结果是：

声学系（传播速度 $0, \pm c$ ）：

\partial_t \rho + \rho \partial_x u = 0, \quad \partial_t (\rho u) + \rho u \partial_x u + \partial_x p = 0

传输系（仅传播速度 $u$ ）：

\partial_t \rho + u \partial_x \rho = 0, \quad \partial_t (\rho u) + u \partial_x (\rho u) = 0

声学系用 backward Euler 求解， $\pm c$ 就从时间步约束里脱离。传输系做显式迎风，CFL 条件仅需 $|u| \Delta t / \Delta x < 1$ 。

Riemann 解算器的正性条件#

还有一个关键细节。声学系是非守恒形式 $\partial_t V + \vartheta \partial_x G(V) = 0$ （ $\vartheta = 1/\rho$ ）。要做 Godunov 型方法就得构造 simple Riemann solver。论文采用左、右两个独立斜率 $\bar{a}_-,\,\bar{a}_+$ （扩展自 Gallice 2003 的解算器）。当斜率足够大时，中间状态的比体积和压力都能保持为正。

正性条件最终归为两个不等式（Appendix B）：

\bar{a}_-^2 \ge \rho_l c_l^2 - \rho_l (p_r - p_l) / \cdots

实操中取 $\bar{a}_\pm = K \max(\rho_l c_l, \rho_r c_r)$ ，把 $K$ 调到所需的阈值。下方的滑块可以直接调 $K$ 。当红色区域（出现负值的中间状态）消失时，就是正性保持的临界。

slope ā₋ / (ρc)_max1.60slope ā₊ / (ρc)_max1.60

Sod-like reference: ρ_L=1, p_L=1, ρ_R=0.125, p_R=0.1. Increase ā₋, ā₊ until the entire (Δu, Δp) square turns cyan — that's Peluchon's positivity-preserving slope choice.

把这个解算器隐式化只需把时间索引从 $n$ 换成 $n+1-$ 。非线性系统用 Picard 迭代求解即可。

自己实现 — 线性声-传输#

在搬运论文全套代码之前，可以先用 线性声-传输系 复现加速效果的核心。

\partial_t \rho + \rho_0 \partial_x u + u_0 \partial_x \rho = 0, \quad \partial_t u + \frac{c^2}{\rho_0} \partial_x \rho + u_0 \partial_x u = 0

特征值为 $u_0 - c,\,u_0,\,u_0 + c$ 。Unsplit 显式 Rusanov 要求 $\Delta t \le \Delta x / (|u_0| + c)$ 。Split 版本把声学（±c）隐式化，把传输（ $u_0$ ）保持显式。

import numpy as np
 
def step_unsplit_rusanov(rho, u, rho0, c, u0, dx, dt):
    """整体 Rusanov — CFL 受 |u0|+c 约束"""
    lam = abs(u0) + c
    F_rho = 0.5 * (rho0 * (u[:-1] + u[1:]) + u0 * (rho[:-1] + rho[1:])) \
            - 0.5 * lam * (rho[1:] - rho[:-1])
    F_u   = 0.5 * ((c*c/rho0) * (rho[:-1] + rho[1:]) + u0 * (u[:-1] + u[1:])) \
            - 0.5 * lam * (u[1:] - u[:-1])
    rho[1:-1] -= dt/dx * (F_rho[1:] - F_rho[:-1])
    u[1:-1]   -= dt/dx * (F_u[1:]   - F_u[:-1])
    return rho, u
 
def step_imex_split(rho, u, rho0, c, u0, dx, dt_outer):
    """声学(隐) + 传输(显) — 外层 CFL 只看 |u0|"""
    # 用子循环近似 backward Euler 的结果
    n_sub = max(1, int(np.ceil(c * dt_outer / dx)))
    dt_a = dt_outer / n_sub
    for _ in range(n_sub):
        F_rho = 0.5 * rho0 * (u[:-1] + u[1:]) - 0.5 * c * (rho[1:] - rho[:-1])
        F_u   = 0.5 * (c*c/rho0) * (rho[:-1] + rho[1:]) - 0.5 * c * (u[1:] - u[:-1])
        rho[1:-1] -= dt_a/dx * (F_rho[1:] - F_rho[:-1])
        u[1:-1]   -= dt_a/dx * (F_u[1:]   - F_u[:-1])
    # 传输 upwind (假设 u0 > 0)
    a = u0 * dt_outer / dx
    rho[1:-1] -= a * (rho[1:-1] - rho[:-2])
    u[1:-1]   -= a * (u[1:-1]   - u[:-2])
    return rho, u

让两个方案以同样的外层 $\Delta t = 0.5 \Delta x / |u_0|$ 运行，取 $M = u_0/c = 0.05$ 。Unsplit 在第一步就发散，Split 借助子循环保持稳定。从效率看，要达到相同物理时间，unsplit 需要多算 $1/M = 20$ 倍的步数。

直击发散现场#

下方仿真可以亲手调试。选 Explicit unsplit 并把 CFL 推到 1.2 以上，青色曲线立刻爆炸，红色警告亮起。切换到 IMEX，同样的外层 CFL 下依然稳定。把 Mach 数降到 0.05，IMEX 的有效时间步增益就扩大到 20 倍。

SchemeCFL (on |u|)1.40Mach number M0.100

IMEX time-step gain ≈ 1 + 1/M = 11.0×STATUS: stable

观察要点：

explicit：CFL 接近 1.0 已是走钢丝，1.05 就开始长出 spurious oscillation。
split：外层 CFL 调到 4 仍能让声学信息通过子循环完整传播，形状保持稳定。
M 越小，两者的时间步差距按比例放大。

批判 — 论文淡化的代价#

加速并非没有代价。

隐式声学求解的成本：backward Euler 是非线性的，需要 Newton 或 Picard。1D tridiagonal Thomas 算法很轻量；2D 结构网格要么走 ADI，要么面对 sparse solve。论文也推荐 ADI 风格的 IM3 变体。
分裂的一阶误差：算子分裂只有一阶精度。本文把二阶推广留到下一篇。2D 液-气相互作用算例里界面明显涂抹。
声波的隐式衰减：隐式声学步具有耗散性，声波会被阻尼。如果出口边界不正确，连真实的声波也会被非物理地抹掉，建议与 NSCBC（上周文章）搭配。
正性斜率的代价：为保证正性而增大 $\bar{a}_\pm$ 会让 Riemann 解算器本身更耗散。最终必须配合 low-Mach 修正（Appendix D）才能恢复精度。

OpenFOAM 的 compressibleInterFoam 通过 PIMPLE 迭代获得类似的声学隐式效果，但没有做特征级分裂。Peluchon 方案在大 CFL 下更稳健，正是因为这一点。

可复现性评分#

论文正文：算法 1、2 提供完整伪代码 → ★★★★☆
推导缺口：Appendix B 的正性不等式略紧，建议自行重新推导 → ★★★☆☆
代码公开：无，需自行实现。1D 约 200 行，2D ADI 约 600 行 → ★★★☆☆
验证算例：shock tube、droplet convection、液-气相互作用、low-Mach 声脉冲，全部贴近实际 → ★★★★★

如果你正打算在周六下午动手把五方程多相模型隐式化，这篇论文是最诚实的起点。一把剪刀就能让声速远离时间步分母。

接下来该读的论文#

Coquel–Godlewski–Khalfallah（2016）—— 单相声-传输分裂的鼻祖
Boscheri–Dumbser（2021）—— 将 IMEX 推广到 all-Mach Navier–Stokes
Bharate（2025）—— 同样的分裂思想用到六方程 Baer–Nunziato 模型