负密度从哪里渗出

只留下让下一步每个单元密度仍高于下限的最大θ，再大就是渗漏。

凌晨2点27分。Le Blanc激波管在0.003秒后吐出NaN。CFL降到0.1，网格加倍，两条曲线还是一起塌了。原因很小。强膨胀波头部，有一个单元的密度掉到了零以下。两步之后声速变成NaN，再下一步整个区域死去。Hu、Adams、Shu(2013)的θ混合，只在那个位置把高阶通量退回一格。

把右侧密度ρ_R降到10⁻³以下，允许曲线进入红色区域，那个面上的θ掉到0.2以下。

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.

红色区域是下一步密度会穿破下限的一侧。θ若停在这条边界上，该面实际上就是一阶Lax–Friedrichs。

高阶格式在稳定性和精度之间晃。一阶迎风把激波抹成大约五个格子。二阶中心差分在激波两侧种下锯齿。振荡轻的时候，只是压力抖一下。相邻单元密度比超过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$ 是守恒变量。被减去的那一项是人工黏性。它把邻居平均掉，锯齿就没了。这和1950年故意把黏性涂在激波上的做法是同一个位置。Lax在CFL $\leq$ 1下证明了 $\rho^{n+1} \geq 0$ ，是Hu et al. 定理1的简化。激波强或弱，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]$ 是每个面的权重。 $\theta=1$ 是纯高阶， $\theta=0$ 是纯LF。限制器只做一件事：在每个单元下一步密度仍为正的前提下，取最大的 $\theta$ 。不往整个区域退。只退一格。

LF的结果已经满足 $\rho^{LF}_i \geq \rho_{\min}$ ，所以每个面都求一个最大的 $\theta$ ，使得加上修正 $\Delta F = F^{HO} - F^{LF}$ 之后仍然为正。一个面碰到两侧单元。两边都得在下限之上，于是较小的那个胜出。

前面把ρ_R拉低时，θ掉到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$ 。强膨胀向左走，若没有保护，第一个单元的密度会变成负数。

左右相除，密度比是 $1/0.02 = 50$ ，压力比是 $100/0.05 = 2000$ 。前面把100:1当作危险线的例子，可这一对的密度跳跃还在那条线下面。猛的是压力。 $\gamma = 1.4$ 时声速 $a = \sqrt{\gamma p / \rho}$ ，左边是 $\sqrt{140} \approx 11.83$ ，右边是 $\sqrt{3.5} \approx 1.871$ 。两者之比正好是 $\sqrt{40} \approx 6.32$ 。

把网格数字代进更新式，就知道不是时间步长的错。 $N = 200$ ，所以 $\Delta x = 1/200 = 0.005$ ， $\Delta t = 9.0 \times 10^{-5}$ 。于是 $\Delta t/\Delta x = 0.018$ 。只用左态声速估的库朗数是 $a_L \Delta t/\Delta x \approx 0.213$ 。这和凌晨把CFL降到0.1仍然死掉的情况对得上。由 $\rho^{n+1} = \rho^n - (\Delta t/\Delta x)\Delta F$ ，通量差一旦超过密度余量，步长再小也会穿破下限。

右态初始密度 $0.02$ 到下限 $10^{-3}$ 的余量是 $0.019$ 。要让 $\theta$ 停在1，密度修正 $|d|$ 不能超过 $0.019/0.018 \approx 1.06$ 。 $|d| = 2$ 时上限是 $0.019/(0.018 \times 2) \approx 0.53$ 。余量缩到 $5 \times 10^{-4}$ 时，同样的 $|d| = 1$ 把 $\theta$ 上限压到 $0.0005/0.018 \approx 0.028$ 。那个面几乎就是LF。

按权重把混合写成 $F = (1-\theta)F^{LF} + \theta F^{HO}$ 。 $\theta = 0.53$ 时高阶占53%。 $\theta = 0.028$ 时高阶只占2.8%，其余是LF。代码里的上限就是更新式的同一条不等式。修正搬动密度的量是 $(\Delta t/\Delta x)\,\theta\,|d|$ ，不能超过右侧余量 $0.019$ ，所以 $\theta \le 0.019 / (0.018\,|d|)$ 。 $|d| \approx 1.06$ 时等号落在 $\theta = 1$ 。

压力公式是同一种读法。 $p = (\gamma - 1)(E - \tfrac{1}{2}\rho u^2)$ ， $\gamma - 1 = 0.4$ 。固定 $E$ 和 $\rho$ ，1单位动能把压力削掉0.4。右态 $p = 0.05$ 在动能 $0.05/0.4 = 0.125$ 时碰到0。左态 $p = 100$ 按同一算法可以撑到 $250$ 。只守住密度下限，渗漏会转到密度仍为正、声速却变成虚数的那一侧。

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；打开，就一直算到结束。

关掉限制器开关，红色曲线会在第30步附近把密度拖到0以下。

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.

那一次穿越就是崩溃。青色在同一时刻于 $\rho \approx 0.02$ 附近变平。精度差不到5%。一个活着，另一个死了。

限制器付出的代价

$\theta$ 落到0的面按一阶LF求解。激波前沿和强膨胀的头部在那里变钝。这块区域不到全域的5%，整体精度几乎不动。θ的平均落在0.95到1.0之间。

把 $\rho_{\min}$ 放到 $10^{-12}$ 这种过低的值，限制器会醒得太晚，同样的NaN会回来。实际计算里，下限大约取平均密度的 $10^{-4}$ 。多级RK必须在每一级重新跑限制器。只保住第一级，下一级照样发散。

把下限放到1e-12，限制器会在NaN之后才醒#

$\rho_{\min} = 10^{-12}$ 不是下限，是装饰。对着右态密度0.02，余量几乎还全在，上面那个θ上限仍停在1附近。密度已经掉到 $10^{-8}$ 以下，上限才开始收紧。在那之前， $\gamma p / \rho$ 已经变号。限制器是在单元死掉之后才退回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.

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