Riemann 問題：CFD の心臓#

有限体積法 (Finite Volume Method) において、セル境界の 数値フラックス (numerical flux) をどのように計算するかは、解の正確さと安定性を左右する極めて重要な要素です。この数値フラックスを決定する鍵となるのが Riemann 問題 です。

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$$

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) を比較します。

左右の圧力比を 1000:1 まで上げると、両側の波が shock になる瞬間が見える。