リーマン問題

By M.Sato

June 12, 2026

1. 概要

非粘性で、熱伝導が無い流れ場を対象として一次元のオイラー方程式を使い、リーマン問題の理論解を記述する。

衝撃波管問題を対象として、特性曲線に沿って保存されるリーマン不変量を利用することで、衝撃波、接触面、膨張波の位置と流速、圧力、密度の分布を求めることができることを示す。

また、いくつかのテストケースについて計算を行った結果を示す。

2. リーマン問題

非粘性で、熱伝導の無い圧縮性流体の流れ場を考える。

流体の密度を$\rho$、速度を$u$、圧力を$p$とする。この流れ場に対して、質量保存則と運動量保存則は次式で表される。

$$ \begin{align} & \frac{D \rho}{Dt} + \rho \nabla\cdot u = 0 \notag \\ & \frac{D u}{Dt} + \frac{1}{\rho} \nabla p =0 \notag \end{align} $$

ここで流れが等エントロピーであると仮定すると圧力変化$dp$は密度変化$d\rho$に対して$dp=a^2d\rho$と表される。ここで、$a$は音速である。従って質量保存式は以下のように表される。

$$ \begin{align} \frac{1}{\rho a} \frac{D p}{Dt} + a \nabla\cdot u = 0 \notag \end{align} $$

一次元流れ場の場合、$\rho,u,p$は、$x$と$t$の関数であることから、上の式は以下のように表される。

$$ \begin{align} & \frac{1}{\rho a} \left( \frac{\partial p}{\partial t} + u \frac{\partial p}{\partial x} \right) + a \frac{\partial u}{\partial x} = 0 \notag \\ & \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + \frac{1}{\rho} \frac{\partial p}{\partial x} =0 \notag \end{align} $$

これらの二式の和、および差をとると次式が得られる。

$$ \begin{align} \left[ \frac{\partial u}{\partial t} + (u+a) \frac{\partial u}{\partial x} \right] +\frac{1}{\rho a}\left[ \frac{\partial p}{\partial t} + (u+a) \frac{\partial p}{\partial x} \right] &= 0 \notag \\ \left[ \frac{\partial u}{\partial t} + (u-a) \frac{\partial u}{\partial x} \right] -\frac{1}{\rho a}\left[ \frac{\partial p}{\partial t} + (u-a) \frac{\partial p}{\partial x} \right] &= 0 \notag \end{align} $$

$\frac{\partial}{\partial t}+(u+a)\frac{\partial}{\partial x}$は特性方向$u+a$に対する物質微分を表しており、微小時間$dt$内での変化量は次式で表される。

$$ \begin{align} du &= \left[ \frac{\partial u}{\partial t} + (u+a) \frac{\partial u}{\partial x} \right]dt \notag \\ dp &= \left[ \frac{\partial p}{\partial t} + (u+a) \frac{\partial p}{\partial x} \right]dt \notag \end{align} $$

従って特性方向$u+a$に対する微分方程式が次式で表される。

$$ \begin{align} du + \frac{1}{\rho a} dp = 0 \notag \end{align} $$

等エントロピを仮定すると、$p/\rho^\gamma=C$（$C$は定数）と表されることから上式は積分することができ、次式が得られる。

$$ \begin{align} J_+ &= u + \frac{2a}{\gamma-1} = \mathrm{const.} \quad \text{along} \quad C_+ : u+a \notag \end{align} $$

これは、特性方向$C_+=u+a$に対して、リーマン不変量$J_+$が保存されることを示している。

同様に$\frac{\partial}{\partial t}+(u-a)\frac{\partial}{\partial x}$は特性方向$C_-=u-a$に対する物質微分を表しており、特性方向$C_-$に対するリーマン不変量$J_-$を導くことができる。

$$ \begin{align} J_- &= u - \frac{2a}{\gamma-1} = \mathrm{const.} \quad \text{along} \quad C_- : u-a \notag \end{align} $$

リーマン不変量$J_+,J_-$が決まれば、流速$u$と音速$a$は次式で求めることができる。

$$ \begin{align} a &= \frac{\gamma-1}{4}(J_+ - J_-) \notag \\ u &= \frac{1}{2}(J_+ + J_-) \notag \end{align} $$

初期状態（密度$\rho$、流速$u$、圧力$p$）が与えられた場合、特性方向に沿ってリーマン不変量を計算していくことで、流れの時間空間分布を計算することが可能となる。

💡 リーマン不変量は特性曲線に沿って保存されるため、初期条件さえ与えれば代数的な関係式だけで流れ場の厳密解を構成できる。この厳密解は、数値解法（リーマンソルバー）の検証用ベンチマークとしても広く使われている。

3. 衝撃波管の計算

衝撃波管の計算は、リーマン問題の解析を利用して行うことができる。

衝撃波管モデルでは、初期時刻において隔壁で分割された管の両側に異なる状態の気体が封入されている。
隔壁が取り除かれると、衝撃波(shock wave)、接触面(contact surface)、膨張波(expansion wave)などを伴う流れ場が形成される。ここでは$x=0$に設置された隔壁の左側に高圧気体、右側に低圧気体が入っている場合を考える。

図1 — 衝撃波管モデル

隔壁が取り除かれると左右の気体の接触面は速度$u_p$で移動することになる。また低圧気体側には衝撃波(shock wave)が移動速度$W$で右側に移動する。

一方左側の領域には膨張波(expansion wave)が伝播することになる。衝撃波は厚みが極めて小さいが膨張波は徐々に幅が広がっていく。

3.1. 波動線図

時刻ゼロで隔壁を除去した場合に、衝撃波、膨張波が時間とともに移動する様子を以下の波動線図に示す。

流体の領域を①、②、③、④で区別し、各領域における圧力$p$、温度$T$、流速$u$の典型的な挙動も合わせて示す。

図2 — 衝撃波管問題の波動線図

3.2. 衝撃波を通過する流れ

衝撃波は波速度$W$で右側に移動する。衝撃波が通過する前後では、物理量$(p,\rho,T)$が不連続的に変化することになる。衝撃波が通過した背後では流速は$u_p$であるとする。

図3 — 衝撃波前後の流れ（衝撃波と共に移動する座標系）

衝撃波と共に移動する座標系に基づき、質量、運動量、エネルギーの保存の式を書くと以下のようになる。

$$ \begin{align} \rho_1 W &= \rho_2 (W-u_p) \notag \\ p_1 + \rho_1 W^2 &= p_2 + \rho_2 (W-u_p)^2 \notag \\ e_1 + \frac{p_1}{\rho_1} + \frac{W^2}{2} &= e_2 + \frac{p_2}{\rho_2} + \frac{(W-u_p)^2}{2} \notag \end{align} $$

これらの式に基づき、ランキンウゴニオ式が導かれる。

$$ \begin{align} e_2 - e_1 = \frac{p_1+p_2}{2}\left(\frac{1}{\rho_1}-\frac{1}{\rho_2}\right) \notag \end{align} $$

内部エネルギーが$e=\frac{p}{\rho(\gamma-1)}$であり、理想気体の状態方程式より$\rho=\frac{p}{RT}$と表せることを利用すると次式が得られる。

$$ \begin{align} \frac{T_2}{T_1} &= \frac{p_2}{p_1}\left[ \frac{\frac{\gamma+1}{\gamma-1}+\left(\frac{p_2}{p_1}\right)}{1+\frac{\gamma+1}{\gamma-1}\left(\frac{p_2}{p_1}\right)} \right] \notag \\ \frac{\rho_2}{\rho_1} &= \frac{1+\frac{\gamma+1}{\gamma-1}\left(\frac{p_2}{p_1}\right)}{\frac{\gamma+1}{\gamma-1}+\left(\frac{p_2}{p_1}\right)} \notag \end{align} $$

また、衝撃波移動速度$W$、ピストン速度$u_p$は次式で表される。

$$ \begin{align} W &= a_1 \sqrt{\frac{\gamma+1}{2\gamma}\left(\frac{p_2}{p_1}-1\right)+1} \notag \\ u_p &= \frac{a_1}{\gamma}\left(\frac{p_2}{p_1}-1\right)\left[\frac{\frac{2\gamma}{\gamma+1}}{\frac{p_2}{p_1}+\frac{\gamma-1}{\gamma+1}}\right]^{\frac{1}{2}} \notag \end{align} $$

ここで$a_1=\sqrt{\gamma R T_1}$は領域①における音速である。

3.3. 膨張波の理論式

波動線図において膨張波は原点（$x=t=0$）から放射（ファン）状に広がっている。膨張波領域（領域③と④の間）における流れの特性を記述する。

図4 — 膨張波（ファン）領域と特性曲線

特性方向$C_+,C_-$については、それぞれリーマン不変量$J_+,J_-$が保存される。図中の青線は特性曲線$C_+$であり、$J_{+,A}=J_{+,C}$である。
また同様に$J_{+,B}=J_{+,D}$が成り立つ。 $A$点と$B$点は領域④に属しており、初期状態が同じであることから $J_{+,A}=J_{+,B}$が成り立つ。従って $J_{+,C}=J_{+,D}$が成り立つことになる。
一方図中の赤線は特性曲線$C_-$であり、この線上では$J_{-,C}=J_{-,D}$が成り立つ。
これらの関係を使うと結局、$u_C=u_D$および$a_C=a_D$が成り立つ事がわかる。

次に特性曲線$C_+$上でのリーマン保存量を考えると次式が成り立つ。

$$ \begin{align} u + \frac{2a }{\gamma-1} =u_4 + \frac{2a_4 }{\gamma-1} =\frac{2a_4 }{\gamma-1} \notag \end{align} $$

変形すると次式が得られる。これは膨張波（ファン）領域における流速と音速を関連付ける式である。

$$ \begin{align} \frac{a}{a_4} = 1 - \frac{\gamma-1}{2}\frac{u}{a_4} \notag \end{align} $$

一方、膨張波は特性方向$C_-$に沿って伸びている。従って特性線は次式で表される。

$$ \begin{align} \frac{dx}{dt}=u-a \notag \end{align} $$

膨張波の始点が$x=0,t=0$を通ることから、特性曲線は次式で表される。

$$ \begin{align} x = (u-a)t \notag \end{align} $$

上式より、流速$u$を$x,t$の関数として表すことができる。

$$ \begin{align} u = \frac{2}{\gamma+1}\left[a_4 + \left(\frac{x}{t}\right)\right] \notag \end{align} $$

膨張波（ファン）の領域では等エントロピー流れであることから、次式が得られる。

$$ \begin{align} \frac{T}{T_4} &= \left[ 1-\frac{\gamma-1}{2}\left(\frac{u}{a_4}\right)\right]^2 \notag \\ \frac{p}{p_4} &= \left[ 1-\frac{\gamma-1}{2}\left(\frac{u}{a_4}\right)\right]^{\frac{2\gamma}{\gamma-1}} \notag \\ \frac{\rho}{\rho_4} &= \left[ 1-\frac{\gamma-1}{2}\left(\frac{u}{a_4}\right)\right]^{\frac{2}{\gamma-1}} \notag \end{align} $$

3.4. 解析手順

上記で説明した関係式を使うことで衝撃波の流れを記述することができる。
接触面において$u_2=u_3$、および$p_2=p_3$であることを利用すると最終的に以下の非線形方程式が得られる。
なお、$p_1,p_4,a_1,a_4$は既知である。ここで未知数は$p_2/p_1$であり、例えばニュートン法などを使って$p_2/p_1$を求める。
この値が求まれば後は芋づる式に他の物理量を計算することができる。

$$ \begin{align} \frac{p_4}{p_1}=\frac{p_2}{p_1}\left[ 1 - \frac{(\gamma-1)\left(\frac{a_1}{a_4}\right)\left(\frac{p_2}{p_1}-1\right)}{\sqrt{2\gamma\left[2\gamma+(\gamma+1)\left(\frac{p_2}{p_1}-1\right)\right]}}\right]^{-2\gamma/(\gamma-1)} \notag \end{align} $$

4. 解析例

以下に示す5ケースについて、リーマン問題の解析を行った。なお比熱比$\gamma$はすべてのケースについて1.4とした。

これらの計算ケースについては、以下からダウンロードできるpythonコードを利用して計算した。これは元のpythonコードを参考にして若干修正作成したものになっている。

このコードは汎用的に作成されており、隔壁から左右に向かう衝撃波、膨張波の任意の組合せを計算できるようになっている。いくつかの解析条件はスライダーバーを使い対話的に変更することができるようになっている。

⬇️ ダウンロード

リーマン問題の計算コード

test	$\rho_L$	$u_L$	$p_L$	$\rho_R$	$u_R$	$p_R$	$t_{end}$
1	1.0	0.0	1.0	0.125	0.0	0.1	0.25
2	1.0	-2.0	0.4	1.0	2.0	0.4	0.15
3	1.0	0.0	1000.0	1.0	0.0	0.01	0.012
4	1.0	0.0	0.01	1.0	0.0	100.0	0.035
5	5.99924	19.5975	460.894	5.99242	-6.19633	46.0950	0.035

ここで、$\rho_L$、$u_L$、$p_L$は隔壁の左側の密度、速度、圧力を表し、 $\rho_R$、$u_R$、$p_R$は隔壁の右側の密度、速度、圧力を表している。また$t_{end}$は計算終了時間を表している。

なお今回の解析で使用したプログラムの理論式および例題の条件は、この書籍に詳細に記述されている。

4.1. test-1

以下にテスト1の計算結果を示す。これはSodの衝撃波管問題と呼ばれるものに対応している。衝撃波、接触面、膨張波の3つの波が発生していることがわかる。
なお図の横軸は空間座標$x$、縦軸は密度$\rho$、速度$u$、圧力$p$、内部エネルギー$e$を表している。
内部エネルギー$e$は、理想気体の状態方程式を利用して $e=\frac{p}{\rho(\gamma-1)}=\frac{R T}{\gamma-1}$で表すことができるため温度を表していると解釈することができる。

図5 — test-1（Sodの衝撃波管問題）の計算結果

4.2. test-2

このケースは、隔壁で分離された領域においてそれぞれが外側に向かって流動する条件になっている。
隔壁から外側に向かって膨張波が発生し、管の中央付近が低圧になるように変動する。

図6 — test-2の計算結果

4.3. test-3

このケースは、強い衝撃波を伴う流れの例である。右側に向かって衝撃波が移動する流れになっている。
これは後述のtest-5の左半分領域の初期条件としても利用されている。

図7 — test-3の計算結果

4.4. test-4

このケースは、左側に向かって衝撃波が移動する流れになっている。これは後述のtest-5の右半分領域の初期条件としても利用されている。

図8 — test-4の計算結果

4.5. test-5

このケースは、test-3とtest-4の初期条件を組み合わせたものである。

図9 — test-5の計算結果

4.6. インタラクティブデモ

上記のpythonコードと同じ計算をJavaScriptに移植し、ブラウザ上で実行できるようにした。
スライダーで初期条件（左右の密度・速度・圧力）と時刻を変更すると、厳密解がリアルタイムに更新される。
プリセットボタンでtest-1〜test-5の初期条件に切り替えることができる。
破線は初期条件、実線は厳密解を表している。なお比熱比は$\gamma=1.4$に固定している。

5. まとめ

非粘性で、熱拡散が無い流れ場を対象として一次元のオイラー方程式を使い、リーマン問題の解析を行った。

特性曲線に沿って保存されるリーマン不変量を利用することで、衝撃波、接触面、膨張波の位置と流速、圧力、密度の分布を求めることができることを示した。

また、いくつかの初期条件に対して、これらの波の位置と流速、圧力、密度の分布を計算した。

弊社の解析事例

弊社の流体解析事例については、下記のリンクからご覧ください。

👉 解析事例のご紹介 OpenFOAM

👉 流体解析実績及び事例