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

Riemann ファンと厳密解 Lax 衝撃波管から見る Godunov 法の原理と限界。

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

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$ から左右どちら側か、波が衝撃か膨張かで状態を抽出します。

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

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 つのファンが触れた瞬間に自己相似性は壊れ、保存も壊れます。

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|λ|.

スライダーを 1.0 を超えて動かすと隣のファンが衝突し、赤い警告が出ます。1.0 未満であれば、境界フラックスは時間に対して定数だという仮定がそのまま通ります。

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)として残っています。

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

下のシミュレーションは同じ 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)で表して単調性を保ちつつ精度を上げます。骨組みはそのまま生き残っています。

今日の持ち帰り#

各セル境界で厳密な Riemann ファンを 1 回解けば、衝撃波を特別扱いせずに自然と捕獲できる。
時間ステップはファンが触れない条件 $\Delta t \le \Delta x / \max|\lambda|$ に縛られる。
1 次精度の天井は Godunov 自身の定理で確定した。以後のすべての高次衝撃捕獲スキームは、この骨格の上での再構成の戦いです。