负密度从哪里渗出 — 爆炸波场景下保正通量限制器

凌晨 2 点 27 分。Le Blanc 激波管验证案例在 0.003 秒后吐出了 NaN。把 CFL 调到 0.1,把网格加倍,上面两条曲线一同崩塌。打开调试器后原因很简单 — 在强膨胀波(rarefaction,压力快速下降的区域)前端,有一个单元的密度跌到了零以下。两步之后声速变成 NaN,再下一步整个计算域全部死掉。本文跟随这股负密度从哪里渗出,以及 Hu、Adams、Shu(2013)提出的"一格一格往后退"的 θ 混合通量限制器如何把它拦住。结尾用 50 行 NumPy 把这个技巧打开/关闭,亲眼比较两种结果。

负密度从哪里渗出#

高阶格式一直在稳定性和精度之间走钢丝。一阶迎风把激波抹成五个网格宽,二阶中心差分则在激波两侧种下锯齿振荡。轻微振荡只让压力抖一下。一旦相邻单元的密度比超过 100:1 — Le Blanc 激波管、Sedov 爆轰、超声速喷管尾迹 — 故事就不一样了。

把单元 $i$ 的密度守恒式写出来:

\rho_i^{n+1} = \rho_i^n - \frac{\Delta t}{\Delta x}\bigl( F_{i+1/2}^\rho - F_{i-1/2}^\rho \bigr)

其中 $F^\rho = \rho u$ 是密度通量。如果 $F_{i+1/2}^\rho \gg F_{i-1/2}^\rho$ (质量从一侧大量流出), $\rho_i^{n+1}$ 就可能跌到零下。在一维 Euler 模拟里,强激波后膨胀区域的第一个单元正好处于这个位置。

负密度会引发什么?声速 $a = \sqrt{\gamma p / \rho}$ 在 $\rho \to 0^-$ 时 $a^2 < 0$ 。下一次通量计算中 $\sqrt{\cdot}$ 返回 NaN。一个单元的 NaN 蔓延到两个邻居,再过两步整个网格阵亡。

Lax–Friedrichs — 模糊但永不崩#

1954 年 Lax 提出的 LF 通量从这个陷阱中免疫:

F_{i+1/2}^{LF} = \tfrac{1}{2}\bigl(F_L + F_R\bigr) - \tfrac{1}{2}\alpha\bigl(U_R - U_L\bigr)

$\alpha$ 是左右两个单元的最大信号速度, $U = (\rho, \rho u, E)^T$ 是守恒变量。第二项是人工粘性(artificial viscosity) — 把邻居平均一下,杀掉锯齿振荡。Lax 证明了在 CFL $\leq$ 1 条件下 $\rho^{n+1} \geq 0$ (Hu et al. 定理 1 的简化版)。无论激波强弱,LF 都不崩。

代价是精度。两格宽的激波变成五格,接触间断每个时间步多模糊一格。生产中没人单独用 LF。但 LF 拥有的安全余量是限制器的起点。

θ 混合 — 一格一格往后退#

设高阶通量 $F_{i+1/2}^{HO}$ (WENO、MUSCL、中心差分等)比 LF 精确但可能产生负密度。Hu、Adams 与 Shu 的想法只有一行:

F_{i+1/2} = F_{i+1/2}^{LF} + \theta_{i+1/2}\bigl( F_{i+1/2}^{HO} - F_{i+1/2}^{LF} \bigr)

$\theta \in [0, 1]$ 是逐面(face)确定的混合权重。 $\theta=1$ 为纯高阶, $\theta=0$ 为纯 LF。限制器只做一件事 — 选取最大的 $\theta$ ,使得每个单元下一步的密度仍为正。不要走太远,一格一格往后退。

由于 LF 计算结果保证 $\rho^{LF}_i \geq \rho_{\min}$ ,我们可以逐面求出 $\theta$ 的上限,使得加上修正通量 $\Delta F = F^{HO} - F^{LF}$ 之后仍为正。一个面影响两侧两个单元,二者都必须高于下限,因此较小的那个就是限制。

下面的模拟可以亲自动手体验。

Per-face θ search · single rarefaction face · curve ρ_iⁿ⁺¹(θ) at fixed left face

log₁₀(ρ_R) = -2.00 (ρ_R = 1.00e-2)outflow velocity u_R = 4.00

The cyan curve is ρⁿ⁺¹ after one Euler step as a function of how much central-flux we mix into the right face. The red zone is where ρ has crossed the floor 10⁻³. The yellow dot is the largest θ that still lands above the floor — that is the θ the per-face limiter picks. Drop ρ_R below 10⁻³ or push u_R past 6 and θ collapses toward zero: the limiter falls back to first-order Lax–Friedrichs.

把 ρ_R 降到 10⁻³ 以下,曲线迅速进入红色区域,限制器选取的 θ 跌到 0.2 以下。在那个面上事实上以一阶 LF 求解。

两步限制 — 先密度,后压力#

只保证密度为正还不够。压力变负时求解器同样会死(声速变成虚数)。Hu et al. 因此把限制器分为两遍。

第一遍(密度) 用上式逐面确定 $\theta^\rho$ ,使所有单元下一步的密度高于 $\rho_{\min}$ 。

第二遍(压力) 用 $\theta^\rho$ 固定密度通量后,对动量和能量通量再求 $\theta^p \leq \theta^\rho$ 。压力正值条件

p = (\gamma - 1)\bigl(E - \tfrac{1}{2}\rho u^2\bigr) \geq p_{\min}

为非线性(通常是一元二次),但保守近似(取 $\theta^p$ 为较小根)一行就能完成。

最终用 $\theta = \min(\theta^\rho, \theta^p)$ 以同一权重混合所有守恒通量。一个面、一个 $\theta$ 、所有守恒量保持一致。

NumPy 50 行 — Le Blanc 激波管上自动求 θ#

Le Blanc 问题是 $\rho_L=1, p_L=100$ 对 $\rho_R=0.02, p_R=0.05$ 的强激波管。强膨胀波向左传播,第一个单元的密度跌入负值,是经典的考验算例。

import numpy as np
 
GAMMA = 1.4
RHO_FLOOR = 1e-3
N, NSTEP = 200, 220
dx = 1.0 / N
dt = 9.0e-5
 
def cons_to_prim(U):
    rho, mu, E = U[0], U[1], U[2]   # 守恒 -> 原始变量
    u = mu / rho
    p = (GAMMA - 1.0) * (E - 0.5 * rho * u * u)
    return u, p
 
def euler_flux(U):
    u, p = cons_to_prim(U)
    return np.array([U[1], U[1] * u + p, u * (U[2] + p)])
 
def alpha_face(UL, UR):
    uL, pL = cons_to_prim(UL); uR, pR = cons_to_prim(UR)
    aL = np.sqrt(GAMMA * max(pL, 0.0) / UL[0])
    aR = np.sqrt(GAMMA * max(pR, 0.0) / UR[0])
    return max(abs(uL) + aL, abs(uR) + aR)
 
def density_theta(U, dF, dt_dx):
    """逐面 θ: 让 ρ_LF 保持在下限以上的最大权重。"""
    theta = np.ones(N + 1)
    for f in range(1, N):
        d = dF[0, f]                          # 密度反扩散
        if abs(d) < 1e-14: continue
        rL_room = U[0, f - 1] - RHO_FLOOR     # 左单元余量
        rR_room = U[0, f]     - RHO_FLOOR     # 右单元余量
        cap = (rL_room if d > 0 else rR_room) / (dt_dx * abs(d))
        theta[f] = max(0.0, min(theta[f], cap))
    return theta

(完整代码为 LF 一遍 + 中心差分反扩散一遍 + θ 决策一遍。220 步,200 单元。关闭限制器在第 30 步附近会 NaN;打开则一直存活到结束。)

把求解器沿时间轴展开:

1D strong shock tube · ρ_L=1, p_L=100 / ρ_R=0.02, p_R=0.05 · central + θ·anti-diffusion · 200 cells

t-index = 80 / 220show limited

Red is the raw centered scheme — its rarefaction tail dives below zero around t ≈ 0.006 and the solver blows up. Cyan is the same scheme blended with Lax–Friedrichs by a per-face θ ≤ θ_max so that ρ stays above the floor 10⁻³. Drag the time slider to watch the two solutions separate.

红色曲线(限制器 OFF)在第 30 步附近把 ρ 拖到零以下,那就是崩溃时刻。青色曲线在相同时刻于 ρ ≈ 0.02 附近趋平。精度差小于 5%,但一个活着,另一个死了。

限制器付出的代价#

θ 跌到 0 的面以一阶 LF 求解。也就是激波前沿或强膨胀波的前端精度下降。所幸该区域不到全域的 5%,整体精度几乎不受影响。平均而言 θ 在 0.95~1.0 之间。

两个陷阱。第一,把 $\rho_{\min}$ 取得太小(如 $10^{-12}$ ),限制器会过晚触发,NaN 仍会重现。生产经验取模拟平均密度的 $10^{-4}$ 左右。第二,RK 多级时间积分必须在每一级重新运行限制器 — 仅在第一级保护,第二级照样发散。

明天也要记住的事#

负密度可从任何地方渗出 — 强激波、膨胀、接触间断。
LF 不是答案,但它是安全余量 — 以这个余量为基准,只允许高阶修正一格一格地往前走。
先密度,后压力 — 两遍就够了。更复杂的限制器只是多了几个旋钮而已。

参考文献:X. Y. Hu, N. A. Adams, C.-W. Shu. Positivity-preserving method for high-order conservative schemes solving compressible Euler equations. J. Comput. Phys. 242 (2013) 169–180.

负密度从哪里渗出 — 爆炸波场景下保正通量限制器