厳密なリーマン解を捨てて得たもの — Roe 近似リーマンソルバとエントロピー修正#

一次元 Euler のリーマン問題には厳密解があります。中心圧力の非線形方程式を Newton で一度解けば終わりです。実務の圧縮性コードは、その解を面ごとに使いません。答えは手の中にあります。それでも捨てます。コストです。理想気体の外では厳密解自体が無いこともあります。

左右の密度と速度のスライダーを動かすと、面上の √ρ 重みと三つの波速 $\hat{u}-\hat{c}$, $\hat{u}$, $\hat{u}+\hat{c}$ が計算し直されます。

$\hat{u}-\hat{c}<0<\hat{u}+\hat{c}$ のとき、波の扇は面にまたがります。亜音速です。左右の速度を同じ側へ大きく押すと、三つの波の符号がそろい、超音速に移ります。

厳密解はなぜ現場から姿を消したのか

厳密リーマンソルバ（Godunov）は、面ごとに非線形方程式の根を探します。セルが数百万個あれば、その根探索ループは毎ステップ数百万回まわります。さらに厄介なことに、厳密解は理想気体の状態方程式に縛られています。実在気体や二相混合物に進むと、厳密解そのものが存在しなくなります。

近似リーマンソルバはこの問題を回避します。非線形なリーマン問題を局所的に線形化するのです。反復なしに、一組の代数式でフラックスが得られます。その代わり、厳密性をわずかに手放します。

ヤコビアンを一つの平均に畳み込む

Roe のアイデアはシンプルです。面を挟んだ二つの状態 $q_L$, $q_R$ にかかるフラックス差を、定数行列 $\hat{A}$ 一つで近似します。

ここで $q$ は保存量ベクトル、$f$ は物理フラックスです。この条件は決定的です。もし二つの状態が単一の波（衝撃波あるいは接触不連続）で結ばれているなら、$\hat{A}$ は有限振幅であってもその波を厳密に伝播させます。つまり近似ソルバでありながら、純粋な衝撃波・接触面の前では厳密なのです。

問題は $\hat{A}$ をどう選ぶかです。いい加減な平均を使えば、上の条件は破れてしまいます。

√ρ 平均はどこから来るのか

Roe の答えは、密度の平方根で重み付けした平均です。

ここで $u$ は速度、$H$ は全比エンタルピー、キャレットは面平均値を表します。Roe 音速はこれに続いて $\hat{c}=\sqrt{(\gamma-1)(\hat{H}-\hat{u}^2/2)}$ で得られます。

なぜよりによって $\sqrt{\rho}$ なのでしょうか。状態ベクトル $q$ とフラックス $f$ を、パラメータベクトル $W=\sqrt{\rho}\,(H,u,1)^T$ の成分で書き下すと、どちらも完全な二次式になります。二次式であれば $q$ と $f$ の $W$ に関する微分が線形になり、$\hat{A}$ を組み立てる経路積分がぴったり評価できます。そこから落ちてくるのが、まさにこの √ρ 重み付き平均なのです。

重みの分母は $\sqrt{\rho_L}+\sqrt{\rho_R}$ です。一方の密度が4倍なら平方根は2倍で、その側の速度とエンタルピーが平均を2倍の強さで引きます。算術平均では、有限ジャンプで $\hat{A}(q_R-q_L)=f(q_R)-f(q_L)$ が合いません。√ρ でなければ純粋な衝撃波の前でも波速がずれます。音速は $\hat{c}=\sqrt{(\gamma-1)(\hat{H}-\hat{u}^2/2)}$ です。$\gamma=1.4$ なら $\gamma-1=0.4$ です。全エンタルピーから運動エネルギー $\hat{u}^2/2$ を引いた残りが音速の二乗になります。三つの波速の符号が分かれる場所は、任意法線で解いたオイラー固有構造の三つの波と同じです。

波に分解して、また足し合わせる

フラックスの組み立ては三段階です。まず状態のジャンプ $\Delta q = q_R - q_L$ を、三つの固有ベクトル $\hat{K}_k$ の和に分解します。各成分の大きさが波の強度 $\alpha_k$ です。次に、それぞれの波を固有値 $\hat{\lambda}_k$ の符号に従って風上方向へ流します。

前の項は中心平均、後ろの項は固有値の大きさで重み付けした風上散逸です。この構造はそのままコードへ移せます。検証には Shu–Osher 問題を使います。マッハ 3 の衝撃波が正弦波状の密度場へ突入し、衝撃波の後方に高周波の密度構造を残す問題です。Roe の低い数値散逸が真価を発揮する、まさにうってつけの試験です。

import numpy as np
 
gamma = 1.4
 
def phys_flux(U):                       # 保存量 U=(ρ, ρu, E) → 物理フラックス
    rho = U[0]; u = U[1] / rho; E = U[2]
    p = (gamma - 1) * (E - 0.5 * rho * u * u)
    return np.array([rho * u, rho * u * u + p, u * (E + p)])
 
def roe_flux(UL, UR, delta):            # Roe 近似リーマンフラックス (+ Harten エントロピー修正)
    rhoL, rhoR = UL[0], UR[0]
    uL, uR = UL[1] / rhoL, UR[1] / rhoR
    pL = (gamma - 1) * (UL[2] - 0.5 * rhoL * uL * uL)
    pR = (gamma - 1) * (UR[2] - 0.5 * rhoR * uR * uR)
    HL = (UL[2] + pL) / rhoL; HR = (UR[2] + pR) / rhoR
    sL, sR = np.sqrt(rhoL), np.sqrt(rhoR)      # √ρ 重み
    u = (sL * uL + sR * uR) / (sL + sR)        # Roe 平均速度
    H = (sL * HL + sR * HR) / (sL + sR)        # Roe 平均エンタルピー
    c = np.sqrt((gamma - 1) * (H - 0.5 * u * u))   # Roe 平均音速
    rr = sL * sR
    drho, dp, du = rhoR - rhoL, pR - pL, uR - uL
    alpha = np.array([(dp - rr * c * du) / (2 * c * c),   # 波の強度 α_k
                      drho - dp / (c * c),
                      (dp + rr * c * du) / (2 * c * c)])
    lam = np.array([u - c, u, u + c])          # 固有値 (波速)
    K = np.array([[1, u - c, H - u * c],       # 右固有ベクトル
                  [1, u,     0.5 * u * u],
                  [1, u + c, H + u * c]])
    al = np.abs(lam)
    small = al < delta                         # エントロピー修正: |λ| の下限値
    al[small] = (lam[small] ** 2 + delta ** 2) / (2 * delta)
    diss = (al * alpha) @ K
    return 0.5 * (phys_flux(UL) + phys_flux(UR)) - 0.5 * diss
 
def run_shu_osher(N=400, tmax=1.8, cfl=0.4, delta=0.1):
    x = np.linspace(0, 10, N); dx = x[1] - x[0]
    rho = np.where(x < 1, 3.857143, 1 + 0.2 * np.sin(5 * x))   # 衝撃波 + 正弦波密度
    u   = np.where(x < 1, 2.629369, 0.0)
    p   = np.where(x < 1, 10.33333, 1.0)
    U = np.array([rho, rho * u, p / (gamma - 1) + 0.5 * rho * u * u])
    t = 0.0
    while t < tmax:
        r = U[0]; v = U[1] / r; pp = (gamma - 1) * (U[2] - 0.5 * r * v * v)
        dt = cfl * dx / np.max(np.abs(v) + np.sqrt(gamma * pp / r))
        dt = min(dt, tmax - t)
        F = np.zeros((3, N + 1))
        for i in range(1, N):
            F[:, i] = roe_flux(U[:, i - 1], U[:, i], delta)
        F[:, 0] = phys_flux(U[:, 0]); F[:, N] = phys_flux(U[:, N - 1])
        U[:, 1:N - 1] -= dt / dx * (F[:, 2:N] - F[:, 1:N - 1])
        t += dt
    return x, U[0]
 
x, rho = run_shu_osher()
print(f"t=1.8  min rho={rho.min():.3f}  max rho={rho.max():.3f}")   # -> min~0.81  max~4.08

わずか 40 行あまりで、完結した圧縮性ソルバができあがります。Newton 反復も、厳密リーマン解もありません。delta が、この後で扱うエントロピー修正パラメータです。

Shu–Osher の初期場は $x<1$ で $\rho=3.857143$、$u=2.629369$、$p=10.33333$、右側は $\rho=1+0.2\sin(5x)$、$u=0$、$p=1$ です。$N=400$、区間は $[0,10]$、CFL は 0.4、$t=1.8$、$\delta=0.1$ です。出力は $\rho$ の最小がおよそ 0.81、最大がおよそ 4.08 です。左側の密度 3.857 が最大 4.08 付近まで少し上がり、正弦波の谷は 0.8 付近まで下がります。Newton なしでその範囲が残ったのが、この実行の結果です。

Roe はエントロピーを破る#

Roe ソルバは、すべての波をジャンプとして扱います。膨張波（rarefaction）でさえ、小さな衝撃波の連なりとして近似します。たいていはそれで問題ありません。ところが、膨張波が音速点（sonic point）を含むと事情が変わります。その点で固有値 $\hat{\lambda}_k$ が符号を変え、$0$ を通過するのです。

$|\hat{\lambda}_k|\to 0$ になると、その波の風上散逸が消えてしまいます。数値スキームは、なめらかな扇形を開く代わりに、静止した膨張衝撃波（expansion shock）をそのまま固定してしまいます。これは Rankine–Hugoniot 条件は満たすものの、エントロピー条件を破った非物理的な解です。

処方箋は Harten のエントロピー修正です。$0$ の近くで固有値の大きさに下限を敷きます。

$|\hat{\lambda}_k|=0$ のとき修正式は $\delta/2$ になります。$(0+\delta^2)/(2\delta)=\delta/2$ です。$\delta$ が 0 なら床も 0 で、風上散逸が消えます。コードの delta=0.1 はその床に 0.05 を敷きます。$(|\hat{u}|+\hat{c})$ の 5〜10% は、接触面がつぶれすぎる前の幅です。

$\delta$ のスライダーを 0 から上げると、$u_L=-0.5$、$u_R=1.0$ の Burgers 膨張 $u_t+(u^2/2)_x=0$ で、$x=0$ の青い折れが動くかが見えます。音速点 $u=0$ は区間のちょうど中央です。

$\delta=0$ のときは、$x=0$ に青い折れ — 静止した膨張衝撃波 — がそのまま残ります。$\delta$ を上げると、数値解が琥珀色の厳密膨張扇形の上へすべり降りてきます。局所の $(|\hat{u}|+\hat{c})$ の 5〜10% なら、たいてい安全です。大きく取りすぎると、接触面がつぶれてしまいます。

もっと安く：HLL と HLLC#

Roe が重いと感じるなら、波の数を減らせばよいのです。HLL（Harten–Lax–van Leer）は、左右の音波の二つだけを残します。中央の接触波は捨てます。

ここで $S_L$, $S_R$ は左右の最外側の波速の推定値です。HLL は頑健で、密度の正値性をよく保ちます。その代わり、接触面を鋭くとらえることはできません。中央波がないからです。

HLLC の C は中央波（Contact）を表します。HLL が捨てた波を復活させるのです。波は三つ、定常状態の領域は四つになります。結局この三つは、波の数という一つのスペクトル上に並びます。

現場の既定値は、たいてい HLLC です。接触面を生かしつつ、密度・圧力の正値性を守りやすいからです。Roe は分解能に優れますが、エントロピー修正が必須で、格子に平行な強い衝撃波ではカーバンクル（carbuncle）現象に注意しなければなりません。

δ=0のままなら膨張衝撃波が面に残る#

√ρ 重みは、$q$ と $f$ を二次式にするパラメータ化 $W=\sqrt{\rho}\,(H,u,1)^T$ から来ます。その経路積分がなければ、$\hat{A}$ のジャンプ条件は有限振幅で破れます。音速点で $|\hat{\lambda}_k|\to 0$ をそのまま置くと、Rankine–Hugoniot は合ってもエントロピーは破れます。Harten の床を敷かない Roe は、膨張衝撃波を静止解として固定します。

波をいくつ残すかが次の決定です。HLL は音波を二つだけ残して接触面をつぶし、HLLC は接触波を戻して状態領域を四つにします。Roe は三次元で波が五つまで行きますが、エントロピー修正が必須で、格子に平行な強い衝撃波ではカーバンクルを見る必要があります。低Machでは人工散逸のスケールが別に分かれます。Hope–Collinsが低MachのRoeを三つのスケーリングに分けた記録が、その分かれ目です。