負の密度が漏れる場所 — 爆風波のためのPositivity保存フラックスリミッター

更新後の密度が床より上に残る範囲で、最大のθだけを残します。

午前2時27分です。Le Blanc衝撃管が0.003秒でNaNを吐きました。CFLを0.1まで下げても、格子を倍にしても、二本の曲線が一緒に崩れます。原因は単純でした。強い膨張波(rarefaction、圧力が急に下がる領域)の先頭で、一つのセルの密度が負に落ちていました。二ステップ後に音速が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}$ が負になりえます。1D 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}$ なので、補正 $\Delta F = F^{HO} - F^{LF}$ を足しても正のままになる $\theta$ の上限を面ごとに解きます。一つの面は両側のセルに触れます。両方とも床の上にいる必要があるので、小さい方が限界です。

先ほどρ_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]                          # 密度 anti-diffusion
        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が一パス、中心差分のanti-diffusionが一パス、θの決定が一パスです。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.