[论文评论] 为什么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系。

把 $M$ 滑到 $10^{-2}$ 和 $10^{-3}$。三个矩阵的 (1,1) 在同一个 Mach 上拉开多少,格子大小直接显示。

$M=10^{-2}$ 时 convective 的 (1,1) 是 $10^4$,mixed 是 $10^2$。一百倍。acoustic 的 (2,2) 按 $M^{-1}$,先把对流速度抹掉。

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方程为:

其中 $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$ 对应对流极限,$\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就是这两个陷阱之间的窄路。

格子大小在上面已经分成 $10^4$ 和 $10^2$。补音速点空洞是另一件事,在 给 Roe 线性化加上熵修正的地方。今天的抹平来自 (1,1) 的 $M$ 指数,不是那个修正。

用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滑块感受这一差距。

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

三行输出照字面读。$M=10^{-1}$ 时 convective 0.832,acoustic 0.978,mixed 0.974。间距不超过 0.15。$M=10^{-2}$ 时只有 convective 掉到 0.018,acoustic 仍是 0.974,mixed 仍是 0.971。把 0.018 和 0.974 并排放,一个周期里只有对流尺度丢掉了大部分振幅。$M=10^{-3}$ 时 convective 按三位小数印成 0.000,acoustic 0.973,mixed 0.971。mixed 与 acoustic 的差是 0.004、0.003、0.002,Mach 越小反而越窄。

格子 $N=128$ 时 $\Delta x=1/128$。声速 $a=1/M$,一个周期 $T=M$,时间步是 $0.25\,\Delta x/a$。$M$ 从分子分母里消掉,步数不随 Mach 变。$1/(0.25\Delta x)=512$。三个 Mach 都走同样的 512 步。差别不在步数,在 $a_{11}$。convective 用 $1/M^2$,在 $M=10^{-2}$ 是 $10^4$。mixed 和 acoustic 的 $a_{11}$ 是 $1/M=10^2$。矩阵 (1,1) 的一百倍原样进了阻尼系数。

初始压力是 $p=M\sin(2\pi x)$,扰动振幅本身就是 $O(M)$。表上是 $\max(p)/M$,也就是除以 $M$ 之后的数。0.018 是这个归一化振幅在一个周期后落到 0.018,0.000 是三位小数里消失了。并没有证明它在数学上精确为零。$\gamma=1.4$ 只进左端的压力-速度耦合,这个线性测试里密度固定为 $\rho_0=1$。

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网格模态振荡,本文准确地指出了根源。

对流尺度撞上声波的时候

$M=10^{-2}$ 时 convective 残幅是 0.018,mixed 是 0.971。声波若是答案的一部分,对流矩阵就不能用。若是定常、接近不可压、声波可以抹掉,则 $A^c_{11}\sim M^{-2}$ 反而对。表上分开的就是这两个数。

激波附近若出现网格模态,那是 mixed 的 (2,2) 停在 $O(M^0)$ 留下的痕迹。这时改用 acoustic 矩阵,或把 inf–sup 稳定化和 (1,1) 对齐。把声音放到单独的时间积分里,是 渐近保持 IMEX 留下声波的做法。只要压力场,mixed 就够。要保住 0.971 的残幅,convective 就丢掉。