Riemann 問題：CFD の心臓#

有限体積法 (Finite Volume Method) では、セル境界の 数値フラックス (numerical flux) の置き方が解の正確さと安定性を分けます。鍵は Riemann 問題 です。左右が定数で、不連続は一つです。

圧力比のスライダーを 1 付近から 1000:1 へ動かすと、左右の波が膨張から衝撃へ切り替わります。

中央の接触面は圧力比を上げても残ります。左右の波が同時に衝撃になるのは、スター圧力 $p^*$ が左右の圧力よりどちらも大きくなったときです。

Riemann 問題とは？#

初期条件が一つの不連続面を境に、左右で異なる定数値を持つ初期値問題のことです。

$$\mathbf{U}(x, 0) = \begin{cases} \mathbf{U}_L & \text{if } x < 0 \\ \mathbf{U}_R & \text{if } x > 0 \end{cases}$$

1次元 Euler 方程式に対し、この問題の解析解は三つの波 (wave) で構成されます。

1. 左向きの波 (rarefaction または shock)
2. 接触不連続 (contact discontinuity)
3. 右向きの波 (rarefaction または shock)

Godunov のアイデア (1959)#

Godunov の核心となる洞察は、シンプルながらも強力なものでした。

“

セル境界におけるフラックスを求めるには、隣接する二つのセルの平均値を左右の状態量として、Riemann 問題を解けばよい。

有限体積法の更新公式は以下の通りです。

$$\mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t}{\Delta x}\left(\hat{\mathbf{F}}_{i+1/2} - \hat{\mathbf{F}}_{i-1/2}\right)$$

ここで $\hat{\mathbf{F}}_{i+1/2}$ こそが、Riemann 問題の解から決定される数値フラックスです。

Exact Riemann Solver#

厳密な Riemann 解法 (Exact Riemann Solver) では、Newton-Raphson 反復法を用いてスター領域 (star region) の圧力 $p^*$ を求めます。左右の波の種類 (shock/rarefaction) に応じて、Rankine-Hugoniot 関係式または等エントロピー関係式を適用します。

衝撃波 (Shock wave) の関係式 ($p^* > p_K$, $K = L$ or $R$):

$$f_K(p^*) = (p^* - p_K)\left[\frac{A_K}{p^* + B_K}\right]^{1/2}$$

ここで $A_K = \frac{2}{(\gamma+1)\rho_K}$、 $B_K = \frac{\gamma-1}{\gamma+1}p_K$ です。

膨張波 (Rarefaction wave) の関係式 ($p^* \leq p_K$):

$$f_K(p^*) = \frac{2c_K}{\gamma - 1}\left[\left(\frac{p^*}{p_K}\right)^{\frac{\gamma-1}{2\gamma}} - 1\right]$$

$p^*$ は以下の条件を満足します。

$$f_L(p^*) + f_R(p^*) + (u_R - u_L) = 0$$

分かれ目は $p^*$ と $p_K$ の大小だけです。$p^*$ が $p_K$ より大きいとその側は衝撃です。$f_K$ は圧力差に $[A_K/(p^*+B_K)]^{1/2}$ を掛けた値です。$A_K=2/((\gamma+1)\rho_K)$ なので、密度が大きいと同じ圧力ジャンプでも速度の変化は小さくなります。$B_K=(\gamma-1)p_K/(\gamma+1)$ は分母を $p^*$ より少し大きくします。$p^*$ が $p_K$ 以下なら膨張です。前の係数 $2c_K/(\gamma-1)$ は、音速が大きいほど、そして $\gamma$ が 1 に近いほど大きくなります。圧力比の指数は $(\gamma-1)/(2\gamma)$ です。$\gamma$ が下がると指数も下がり、圧力比は $f_K$ に効きにくくなります。

二つの $f$ の和に $u_R-u_L$ を足して 0 になる必要があります。$u_R$ が $u_L$ より小さいと流体は中央へ集まります。和 $f_L+f_R$ は正で、$p^*$ は左右の圧力より上へ押されます。そのとき左右の波は両方とも衝撃です。画面の 1000:1 はその端です。低い側の $p_K$ はスター圧力より小さくなり、その $f_K$ は衝撃の枝に残ります。高い側は圧力比次第で膨張のままもあり得ます。両側が衝撃になるには、中央へ寄る速度差が十分大きくなければなりません。セル一つをそのまま衝撃管にする計算は、Godunov がセルの中で解く衝撃管にあります。

Approximate Riemann Solver：なぜ必要なのか？#

厳密解法は反復計算が必要なため、計算コストが高くなります。実際の CFD コードでは、その多くが 近似 Riemann 解法 (Approximate Riemann Solver) を使用しています。

Roe Solver (1981)#

Roe は非線形な Riemann 問題を 線形化 (linearization) し、行列 $\hat{\mathbf{A}}$ に対する線形 Riemann 問題へと変換します。

$$\hat{\mathbf{F}}_{i+1/2} = \frac{1}{2}(\mathbf{F}_L + \mathbf{F}_R) - \frac{1}{2}\sum_{k=1}^{3}|\hat{\lambda}_k|\hat{\alpha}_k\hat{\mathbf{r}}_k$$

Roe 平均 (Roe average) を用いて $\hat{\mathbf{A}}$ を構成します。

$$\hat{u} = \frac{\sqrt{\rho_L}\,u_L + \sqrt{\rho_R}\,u_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}, \quad \hat{H} = \frac{\sqrt{\rho_L}\,H_L + \sqrt{\rho_R}\,H_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}$$

長所: 接触不連続を正確に捕捉。 短所: エントロピー条件に違反する可能性（エントロピー修正が必要）。

HLLC Solver (Toro, 1994)#

HLL 解法を改良し、 接触不連続 (Contact) を復元 したものが HLLC です。三つの波の速度 (wave speed) $S_L, S^*, S_R$ を推定します。

$$\hat{\mathbf{F}}_{i+1/2} = \begin{cases} \mathbf{F}_L & \text{if } S_L > 0 \\ \mathbf{F}^*_L & \text{if } S_L \leq 0 < S^* \\ \mathbf{F}^*_R & \text{if } S^* \leq 0 \leq S_R \\ \mathbf{F}_R & \text{if } S_R < 0 \end{cases}$$

スター領域のフラックスは以下の通りです。

$$\mathbf{F}^*_K = \mathbf{F}_K + S_K(\mathbf{U}^*_K - \mathbf{U}_K)$$

HLLC は実装がシンプルでありながら頑健 (robust) であるため、 圧縮性多相流コードで最も広く利用されています 。

AUSM+ (Liou, 1996)#

AUSM 系は、フラックスを対流項 (convective) と圧力項 (pressure) に分離します。

$$\hat{\mathbf{F}}_{i+1/2} = \dot{m}_{1/2}\boldsymbol{\Psi}_{1/2} + \hat{p}_{1/2}\mathbf{D}$$

ここで $\dot{m}_{1/2}$ はセル境界の質量流量、 $\boldsymbol{\Psi}$ は対流される物理量です。低マッハ数流でも安定しているため、全速度域スキーム (all-speed scheme) への拡張が容易です。

多相流への拡張

多相流では界面の両側で状態方程式 (EOS) が異なるため、Riemann 解法を マルチマテリアル (multi-material) 版へと一般化する必要があります。

代表的なアプローチには以下があります。

• Ghost Fluid Method: 界面付近で反対側の流体のゴーストセルを構成し、各相において独立したシングルマテリアルの Riemann 問題を解きます。
• HLLC for multi-material: スター領域において両側の EOS をそれぞれ適用し、 $p^*, u^*$ を決定します。

次の記事では、 界面捕捉法 (VOF, Level Set, Diffuse Interface) を比較します。

多相になって移流と圧力を切り離す分離は、Baer–Nunziato 流束の TV 分離へ続きます。