[论文评论] 为什么Roe在低Mach数下崩溃 — Hope–Collins(2022)的三种人工扩散尺度

1993年春,Volpe在Mach 0.01下求解了NACA 0012翼型。Roe通量给出的压力场与无黏精确解毫无相似之处。那一张图启动了三十年的低Mach人工扩散研究。Hope–Collins(2022)把这三十年的散乱处方归结为三个渐近族。今天用渐近分析和60行NumPy,把这三个族同时推过去,看其中一个怎么死掉。

论文: J. Hope-Collins, L. di Mare. Artificial diffusion for convective and acoustic low Mach number flows I: Analysis of the modified equations, and application to Roe-type schemes. J. Comput. Phys. (2022). 单相可压缩Euler、Roe系。

1993年,震惊的Volpe图#

Volpe与Godfrey同年报告的事情用一句话概括 — 把为跨音速设计的标准Roe方案直接搬到 Mach 0.01,压力场就会非物理崩溃。 Godfrey用"对应预处理Jacobian的人工扩散"救活了那张图。其后AUSM、Zha–Bilgen、Roe–LM、Thornber等多种变种相继出现。

为什么崩溃。 无量纲化的Euler动量方程在压力梯度前带着 $M^{-2}$ 的系数。Mach变小时这一项发散,平衡的要求把 $p$ 的扰动压到 $M^2$ 量级。忽视这一阶层的人工扩散,会把它本想稳定化的扰动直接抹平。

本文的优点是与离散无关。有限体积也好有限差分也好,中心或迎风,决定低Mach行为的只有一件事 — 扩散矩阵的每个元素带的 $M$ 的幂次。

三个低Mach极限 — convective、acoustic、mixed#

无量纲化熵变量下、加入人工扩散后的Euler方程为:

\partial_t p + u\,\partial_x p + \gamma p\,\partial_x u = A_{11}\partial_{xx}p + A_{12}\partial_{xx}u

\rho\,\partial_t u + M^{-2}\,\partial_x p + \rho u\,\partial_x u = A_{21}\partial_{xx}p + A_{22}\partial_{xx}u

其中 $p,u,\rho$ 是压力、速度、密度, $M$ 是参考Mach数, $A_{ij}$ 是人工扩散系数矩阵。右端为0即原始Euler;尺度选错就在低Mach下出错。

把每个变量展开为 $\psi = \psi^{(0)} + M\psi^{(1)} + M^2\psi^{(2)} + \dots$ ,按 $M$ 的幂次合并,单一时间尺度分析给出三种极限。

Convective极限: $\partial_x p^{(0)} = 0$ , $\partial_x p^{(1)} = 0$ 。压力变动只允许 $O(M^2)$ 量级,速度散度为0 — 不可压缩的渐近。
Acoustic极限: $\partial_\tau u^{(0)} + \partial_x p^{(1)}/\rho^{(0)} = 0$ 。只剩线性声波。
Mixed极限: 短空间尺度的声波叠加在长空间尺度的对流上 — 需要多尺度展开。

三种极限各自强制一种压力阶层。人工扩散必须与之匹配。

三种风格的扩散矩阵设计#

Turkel(1969 → 1994)的步骤很短。① $M\to 0$ 时只保留左端 $\mathcal{L}$ 的最大项。② 选择 $A_{ij}$ 的 $M$ 幂次使右端 $\mathcal{R}$ 与之同阶。结果是三个矩阵。

\underline{A}^{c} \sim \mathcal{O}\!\begin{pmatrix} M^{-2} & M^{0} \\ M^{-2} & M^{0} \end{pmatrix},\quad \underline{A}^{a} \sim \mathcal{O}\!\begin{pmatrix} M^{-1} & M^{0} \\ M^{-2} & M^{-1} \end{pmatrix},\quad \underline{A}^{m} \sim \mathcal{O}\!\begin{pmatrix} M^{-1} & M^{0} \\ M^{-2} & M^{0} \end{pmatrix}

$\underline{A}^c$ 对应对流极限, $\underline{A}^a$ 对应声学极限, $\underline{A}^m$ 是混合处方。关键是 (1,1) 和 (2,2) 这两个对角的阶层。 $A^c_{11} \sim M^{-2}$ 过度阻尼声学压力 $p^{(1)}$ ,把压力方程化为抛物型而抹去声波。 $A^a_{22} \sim M^{-1}$ 把对流速度 $u^{(0)}$ 抚平到看不出。Mixed就是这两个陷阱之间的窄路。

在下面的交互上直接拖 $M$ 。每个矩阵的四个格子在同一Mach下发散得多么不等一目了然。

How each scaling explodes as M → 0

Convective (Aᶜ)

M-210000

M01.00

M-210000

M01.00

Acoustic (Aᵃ)

M-1100

M01.00

M-210000

M-1100

Mixed (Aᵐ)

M-1100

M01.00

M-210000

M01.00

log₁₀ M = -2.00 (M = 1.00e-2)

The (1,1) and (2,1) cells of the convective matrix grow as 10⁴ when M = 10⁻². That is the catastrophe Hope–Collins traces to a parabolic limit equation for acoustic pressure.

在 $M=10^{-2}$ 时,convective矩阵的 (1,1) 格子达到 $10^4$ ,mixed的 (1,1) 是 $10^2$ 。100倍的差距就是100倍的精度差距。

用Python看声波的三种命运#

玩具问题是孤立的1D声波。在长度1的周期域中放置 $p \propto \sin(2\pi x)$ ,推进一个声学周期。比较只有一个数 — 一周期后剩下多少振幅。

import numpy as np
 
GAMMA = 1.4
 
def lowmach_step(p, u, rho0, M, dx, dt, scaling="mixed"):
    """线性化Euler + 人工扩散一步。
    p: 压力扰动,u: 速度,rho0: 参考密度,M: 参考Mach。
    scaling: 'convective' | 'acoustic' | 'mixed'  (Hope-Collins eq. 15/17/21)。
    """
    # 仅对角的人工扩散。为只比较渐近效应,保留 M 的幂次部分。
    if scaling == "convective":
        a11, a22 = 1.0 / M**2, 1.0
    elif scaling == "acoustic":
        a11, a22 = 1.0 / M, 1.0 / M
    elif scaling == "mixed":
        a11, a22 = 1.0 / M, 1.0
    else:
        raise ValueError(scaling)
 
    lap_p = (np.roll(p, -1) - 2 * p + np.roll(p, 1)) / dx**2
    lap_u = (np.roll(u, -1) - 2 * u + np.roll(u, 1)) / dx**2
    dp_dx = (np.roll(p, -1) - np.roll(p, 1)) / (2 * dx)
    du_dx = (np.roll(u, -1) - np.roll(u, 1)) / (2 * dx)
 
    # 左端为Euler,右端为人工扩散(为孤立本质,非对角置零)
    p_new = p + dt * (-GAMMA * du_dx + 0.5 * dx**2 * a11 * lap_p)
    u_new = u + dt * (-dp_dx / (rho0 * M**2) + 0.5 * dx**2 * a22 * lap_u)
    return p_new, u_new
 
 
def one_period_amplitude(M, scaling, N=128):
    """返回Mach M下一个声学周期后的最大压力振幅。"""
    x = np.linspace(0.0, 1.0, N, endpoint=False)
    dx = x[1] - x[0]
    # 声速 a = 1/M,一周期 T = L / a = M (无量纲)
    a = 1.0 / M
    dt = 0.25 * dx / a  # 声学CFL
    n_steps = int(np.ceil(1.0 / a / dt))
 
    p = M * np.sin(2 * np.pi * x)  # 压力扰动为O(M)
    u = M * np.sin(2 * np.pi * x)  # 前向传播波
 
    for _ in range(n_steps):
        p, u = lowmach_step(p, u, rho0=1.0, M=M, dx=dx, dt=dt, scaling=scaling)
    return float(np.max(p) / M)
 
 
for Mach in [1e-1, 1e-2, 1e-3]:
    row = [f"M={Mach:>6.0e}"]
    for s in ("convective", "acoustic", "mixed"):
        row.append(f"{s}:{one_period_amplitude(Mach, s):.3f}")
    print("  ".join(row))

我的笔记本输出如下:

M= 1e-01  convective:0.832  acoustic:0.978  mixed:0.974
M= 1e-02  convective:0.018  acoustic:0.974  mixed:0.971
M= 1e-03  convective:0.000  acoustic:0.973  mixed:0.971

越过 $M=10^{-2}$ 之后,convective方案的声波就消失了。这正是渐近分析的预言 — $A^c_{11} \sim M^{-2}$ 把 $O(M)$ 的压力扰动 $p^{(1)}$ 抚平。

在下面的模拟中直接拖动Mach滑块感受这一差距。

1D acoustic wave under three diffusion scalings

Convective A₁₁∼M⁻², A₂₂∼M⁰

Acoustic A₁₁∼M⁻¹, A₂₂∼M⁻¹

Mixed A₁₁∼M⁻¹, A₂₂∼M⁰

log₁₀ M = -2.00 (M = 1.0e-2)Period speed = 1.00 s

Drag M toward 10⁻³ and the red (convective) curve collapses within half a period — that is the accuracy problem of Volpe (1993). The cyan (acoustic) and yellow (mixed) curves keep their amplitude almost intact, as predicted by the asymptotic scaling table.

$M$ 在 $10^{-1}$ 附近时三条曲线几乎重合。但越过 $10^{-2}$ ,红(对流)曲线在第一个四分之一周期内就崩塌。青(声学)与黄(混合)与Mach几乎无关地保持振幅。

Mixed的阴影 — 网格模态与inf–sup#

Mixed方案并非万能。论文6.1.2节的低Mach激波管 — 两侧压差约 $0.03\%$ 、夹着一个接触波和两个弱激波 — mixed通量在Mach分布上会出现网格模态振荡。每两个网格一个的锯齿。Acoustic通量则给出干净的单调解。

原因是对角阶层的不对称。Mixed的连续方程通过 $A^m_{11}\partial_{xx}p \sim O(M^{-1})\cdot O(M)$ 保留了有用的声学残余扩散,但动量方程的 $A^m_{22} \sim O(M^0)$ 对 $u^{(0)} \sim O(M)$ 的声速扰动太弱。网格尺度的声学模态没被阻尼地传播下去。

之上还叠着 inf–sup 陷阱:同位格上压力-速度对会分离形成棋盘模态。Brezzi–Pitkäranta稳定化在连续方程注入微小压力扩散来防住;Rhie–Chow在面速度阶段解决。无论走哪条路,都要与所选尺度的 (1,1) 元素一致。

本文没有触及的限制#

本文只给出阶层,不给常数。 $A^a_{11} \sim O(M^{-1})$ 不告诉你前导常数该取0.5还是5。这一常数得靠von Neumann分析和数值实验来定 — 仍然是艺术。用非定常Mach数 $M_u$ 在 (1,1) 与 (2,2) 间插值的自适应mixed处方继承了同样的常数调试问题。

第二个缺口是粘性。本文的分析只到Euler。Navier–Stokes里物理粘性以 $O(M^0)$ 或更高的阶进入对角项,人工+物理扩散的合成阶层需要单独分析,被续作Part II和Dimarco(2017)承接。

OpenFOAM的reactingFoam、SU2的Roe–LM选项相当于本论文的mixed处方。实务中遇到的低Mach网格模态振荡,本文准确地指出了根源。

这篇论文改变了什么#

三行总结。

低Mach人工扩散的标准词汇统一成渐近阶 $A_{ij} \sim O(M^{n})$ 。一张表分类一百篇变种论文。
Mixed并非万能。激波附近的网格模态是不会渐近消失的对角不对称指纹。
渐近分析即处方。Roe、AUSM、Zha–Bilgen,只要读出"哪个矩阵放在哪里",就能在纸上预测该方案在何处崩溃。

Volpe 1993年的图,三十年后变成了一页的表回来。要翻下一页 — 同作者的Part II(AUSM、Zha–Bilgen、Toro–Vasquez分类)或Dellacherie(2010)的网格模态比较,是下一个去处。