Godunov 法 — 各セルの中に小さな衝撃波管を一つずつ

星圧が左右の圧力のどちら側かで波を分け、隣の扇が重なるステップは捨てます。

100 セルを切っても、1 次精度の風上スキームは衝撃波を 5 セル幅に塗りつぶしてしまいます。1959 年、Sergei Godunov は逆を選びました。各セル境界に小さな衝撃波管(Riemann 問題)を一つずつ置き、その厳密な自己相似解を 1 時間ステップ流したあと、セル平均を取り直す。本記事では Riemann 問題の 5 領域ファン構造を厳密解で追い、Lax 衝撃波管を 12 行の Python で解き、Godunov 自身が証明した 1 次精度の壁まで歩きます。

スライダーを 1.0 の前後で動かします。隣り合う扇が重なると赤い警告が点き、1.0 未満では界面フラックスが時間について定数のままです。

time step Δt / Δt_CFL = 0.70max wave speed = 1.00

Each cell interface spawns a self-similar Riemann fan. Push the time step past 1 and the fans collide — Godunov requires Δt ≤ Δx / max|λ|.

警告が点いたステップは、自己相似がすでに崩れたステップです。 $\Delta t$ が $\Delta x/\max|\lambda|$ を越えた絵と同じです。

1959 年 — 各セル境界に衝撃波管を一つずつ#

これまでの有限体積法はセル内部を一定と仮定します。すると境界 $x_{i+1/2}$ には左状態 $\mathbf{q}_i$ と右状態 $\mathbf{q}_{i+1}$ の間に跳びが必ず立ちます。その跳びを左に $-\infty$ 、右に $+\infty$ まで伸ばせば、それがそのまま Riemann 問題になります。

Godunov のアイデアは率直です。すべてのセル境界でこの小さな Riemann 問題を解き、その厳密な自己相似解を 1 ステップ進めてからセル平均を取り直す。圧縮性 Euler の固有構造(特性速度 $u$ および $u \pm c_s$ )を分けずに一度に扱う点が肝心です。

Riemann ファン — 5 領域の自己相似解#

$\gamma$ -気体は次に従います。

\partial_t \mathbf{q} + \partial_x \mathbf{f}(\mathbf{q}) = 0, \quad \mathbf{q} = (\rho, \rho u, \rho E)^\top

$\rho$ は密度、 $u$ は速度、 $E$ は総エネルギー、 $\mathbf{f}$ は保存形フラックスです。

初期条件 $(\rho_L, u_L, p_L)$ と $(\rho_R, u_R, p_R)$ が与えられると、解は $\xi = x/t$ のみの関数になります。隔膜から開くファンは正確に 5 領域に分かれます。

左の未擾乱領域 — 左側の波の先頭(head)がまだ到達していない場所。
左の波そのもの — 膨張波なら連続変化、衝撃波なら跳び。
星付き左(star-L)領域 — 速度と圧力は星の値、密度は左由来。
星付き右(star-R)領域 — 同じ星圧・星速度、別の星密度。
右の未擾乱領域。

3 と 4 の境界が接触不連続(contact discontinuity)です。左右で圧力と速度は一致しますが、密度とエントロピーだけが跳びます。

単一変数 $p^{\star}$ の関数方程式#

未知数は 4 つに見えますが、星圧 $p^{\star}$ ひとつを解けば残りは追従します。Toro の古典的な縮約は次の関数方程式に閉じます。

f_L(p^{\star}; \mathbf{q}_L) + f_R(p^{\star}; \mathbf{q}_R) + (u_R - u_L) = 0

各 $f_K$ は衝撃波と膨張波で分岐します。

f_K(p) = \begin{cases} (p - p_K)\sqrt{\dfrac{A_K}{p + B_K}} & p > p_K \;\;\text{(shock)} \\[6pt] \dfrac{2 c_K}{\gamma - 1}\left[\left(\dfrac{p}{p_K}\right)^{\!\frac{\gamma-1}{2\gamma}} - 1\right] & p \le p_K \;\;\text{(rarefaction)} \end{cases}

$A_K = 2/[(\gamma + 1)\rho_K]$ 、 $B_K = \frac{\gamma - 1}{\gamma + 1} p_K$ 、 $c_K = \sqrt{\gamma p_K / \rho_K}$ です。星速度は

u^{\star} = \tfrac{1}{2}(u_L + u_R) + \tfrac{1}{2}\bigl[f_R(p^{\star}) - f_L(p^{\star})\bigr].

Newton 反復で $p^{\star}$ を求めたあと、 $\xi = x/t$ から左右どちら側か、波が衝撃か膨張かで状態を抽出します。その抽出が数値フラックスになる道筋は Riemann 問題から Godunov へ続くフラックスに別にあります。

時間ステップの上限 — ファンが触れる瞬間

Godunov 法の本質的な制約は、隣り合うファンが 1 ステップ内で衝突しないことです。

\Delta t \le \min_i \frac{\Delta x}{\max_k |\lambda_k^{(i+1/2)}|}

固有値 $\lambda_k$ は $u$ および $u \pm c_s$ です。2 つのファンが触れた瞬間に自己相似性は壊れ、保存も壊れます。

上の赤い警告が、この不等式が崩れた場所です。

Python — Lax 衝撃波管の厳密解#

Toro 本の標準手続きを短く写します。Lax 衝撃波管 $(\rho_L, u_L, p_L) = (0.445, 0.698, 3.528)$ 、 $(\rho_R, u_R, p_R) = (0.5, 0, 0.571)$ 。

import numpy as np
 
GAMMA = 1.4
 
def f_wave(p, rho_K, p_K, c_K, gamma=GAMMA):
    # 単一波の圧力-速度関数(衝撃 vs 膨張で分岐)
    if p > p_K:
        A = 2.0 / ((gamma + 1.0) * rho_K)
        B = (gamma - 1.0) / (gamma + 1.0) * p_K
        return (p - p_K) * np.sqrt(A / (p + B))
    return (2.0 * c_K / (gamma - 1.0)) * ((p / p_K) ** ((gamma - 1.0) / (2.0 * gamma)) - 1.0)
 
def f_prime(p, rho_K, p_K, c_K, gamma=GAMMA):
    if p > p_K:
        A = 2.0 / ((gamma + 1.0) * rho_K)
        B = (gamma - 1.0) / (gamma + 1.0) * p_K
        return np.sqrt(A / (p + B)) * (1.0 - (p - p_K) / (2.0 * (p + B)))
    return (1.0 / (rho_K * c_K)) * (p / p_K) ** (-(gamma + 1.0) / (2.0 * gamma))
 
def solve_star(left, right, gamma=GAMMA, tol=1e-9, max_iter=80):
    rho_L, u_L, p_L = left
    rho_R, u_R, p_R = right
    c_L = np.sqrt(gamma * p_L / rho_L)
    c_R = np.sqrt(gamma * p_R / rho_R)
    # 初期推定: 圧力平均を速度跳びで補正
    p = max(1e-6, 0.5 * (p_L + p_R) - 0.125 * (u_R - u_L) * (rho_L + rho_R) * (c_L + c_R))
    for _ in range(max_iter):
        f  = f_wave(p, rho_L, p_L, c_L) + f_wave(p, rho_R, p_R, c_R) + (u_R - u_L)
        df = f_prime(p, rho_L, p_L, c_L) + f_prime(p, rho_R, p_R, c_R)
        dp = f / df
        p_new = max(1e-6, p - dp)
        if 2.0 * abs(p_new - p) / (p_new + p) < tol:
            p = p_new
            break
        p = p_new
    u = 0.5 * (u_L + u_R) + 0.5 * (f_wave(p, rho_R, p_R, c_R) - f_wave(p, rho_L, p_L, c_L))
    return p, u, c_L, c_R
 
# Lax 衝撃波管
left  = (0.445, 0.698, 3.528)
right = (0.5,   0.0,   0.571)
p_star, u_star, c_L, c_R = solve_star(left, right)
print(f"p* = {p_star:.5f}, u* = {u_star:.5f}")
# -> p* = 3.40948, u* = 0.62556

この 12 行の Newton ループが Godunov 1 ステップの心臓です。厳密解はやがて 1981 年 Roe の線形化に置き換わり、後に HLL や HLLC が標準になりました。それでも厳密解は新スキームの検証における基準解(reference solution)として残っています。

印刷された左右の状態だけで星の値を分けます。左圧力 3.528、右 0.571、星圧 3.40948。左からは 0.11852 だけ下がり、右からは 2.83848 上がります。星速度 0.62556 は左の 0.698 より 0.07244 低く、右の 0 より上です。左右速度の算術平均は 0.349。星速度の式の第2項が 0.349 を 0.62556 まで押し上げます。

分岐だけで足ります。星圧が左圧力より小さいので左の $f_L$ は膨張の枝、右圧力より大きいので右の $f_R$ は衝撃の枝です。許容誤差は $10^{-9}$ 、反復上限は 80、圧力の床は $10^{-6}$ 。初期推定は $0.5(p_L+p_R)$ に速度跳びを足します。 $u_R-u_L=-0.698$ なので補正は正です。二圧力の平均は 2.0495。推定はその平均より上から出て 3.40948 に着きます。この引き算はコードに書いた式と出力コメントだけを使っています。

100 セルで 1 次風上が衝撃を 5 セルに伸ばすと、その幅は格子長さの 5% です。スライダー 1.0 は $\Delta t=\Delta x/\max|\lambda|$ の境界です。1 を越えると隣の扇が重なり、界面フラックスが 1 ステップのあいだ定数だという仮定が崩れます。 $\gamma=1.4$ は $A_K$ 、 $B_K$ 、 $c_K$ の三か所に一緒に入ります。

ファンを自分で開いてみる

下のシミュレーションは同じ Newton 反復をブラウザで再実行し、時間スライダーで自己相似ファンを開きます。プリセットを 123 problem(両側へ逃げる真空近くの流れ)にすると中央に深い圧力の谷ができ、Strong-blast ではほぼ単一の衝撃波の形に潰れます。

presettime t = 0.140

Three waves emerge from the origin: a left-going rarefaction or shock, a contact discontinuity at u*, and a right-going wave. The vertical pink dashed line marks the initial diaphragm at x = 0.

t を小さく取り、徐々に大きくしてみてください。左右の波は一定速度で離れていきます。Lax では左側が膨張波(rarefaction)として滑らかに開き、右側が衝撃波(shock)として跳びます。中央の小さな密度段が接触不連続です。

1 次の壁 — Godunov の定理#

Godunov 自身は同じ 1959 年の論文で 線形かつ単調なスキームは 1 次精度より高くなれない ことを証明しています。自分のアイデアが 1 次に縛られるという冷ややかな帰結でした。

この壁を迂回する道筋が 1970–80 年代の高解像度衝撃波捕獲スキームの歴史です。van Leer の MUSCL、Harten の TVD、Roe の線形化、Colella–Woodward の PPM、Shu–Osher の WENO — どれも Godunov の骨格を残し、セル内部を区分一定ではなく非線形再構成(reconstruction)で表して単調性を保ちつつ精度を上げます。骨組みはそのまま生き残っています。

厳密な扇を界面に置くとき

$p^{\star}=3.40948$ は左の 3.528 より小さく、右の 0.571 より大きいので、左は膨張、右は衝撃です。 $u^{\star}=0.62556$ は左速度 0.698 から 0.07244 下がった値です。この一組が界面ごとに要るなら厳密な扇を回します。100 セルで衝撃が 5 セルに広がってもよいなら、1 次風上で足ります。

厳密解は高いです。1981 年以降に Roe、HLL、HLLC がその席を取った理由です。検証用に残る基準は、この $p^{\star}$ と $u^{\star}$ です。衝撃をわざと厚くする反対側は 1950 年の人工粘性が衝撃をぼかしたやり方です。Godunov はぼかさずに捕獲し、その代わり $\Delta t$ が扇に縛られます。