ジェット戦闘機が音速を超えた瞬間に機体の周囲にコーン状の白い雲が現れる光景を見たことがあるでしょうか。あるいは、スペースシャトルが大気圏に再突入する際、機体前方に数千度にも達する高温領域が形成される現象を耳にしたことがあるかもしれません。これらの劇的な現象の正体が衝撃波（shock wave）です。
衝撃波は、流れの速度が音速を超えたときに発生する、圧力・温度・密度が不連続的に急変する極めて薄い領域です。亜音速の日常世界ではほとんど意識しませんが、超音速の世界では衝撃波を理解しなければ何も設計できません。衝撃波は流体力学のなかでも特に美しい理論体系を持ち、わずか数本の保存則から驚くほど豊かな物理が導かれます。
衝撃波と膨張波を理解すると、以下のような幅広い分野で応用が見えてきます。
- 航空宇宙工学: 超音速航空機の翼型設計、ラバルノズルの流れ解析、ロケットエンジンの排気流解析など、超音速飛行体の空力設計には衝撃波理論が不可欠です
- 再突入体の熱防護: 宇宙機が大気圏に再突入する際、機体前方に形成される弓状衝撃波（bow shock）が空力加熱の主要因であり、衝撃波前後の状態量変化を正確に予測することが熱防護設計の出発点です
- 超音速風洞と実験流体力学: 風洞の試験部で所望のマッハ数を実現するためには、ノズル形状の設計にプラントル・マイヤー膨張理論を用います
- 爆発・衝撃工学: 爆発によって生じるブラスト波（blast wave）は球面衝撃波であり、建築物の耐爆設計にランキン・ユゴニオ関係式が応用されます
- 天体物理学: 超新星残骸の膨張、太陽風と地球磁気圏の相互作用（バウショック）、星間ガスの衝撃波加熱など、宇宙規模のスケールでも衝撃波は中心的な役割を果たします
本記事の内容
- 衝撃波の物理的原因 ── 情報伝播速度と流速の関係、マッハ線とマッハ角
- 垂直衝撃波のランキン・ユゴニオ関係式の導出（保存則から一歩ずつ）
- 垂直衝撃波前後の状態量変化とエントロピー増大
- 斜め衝撃波の理論と $\theta$-$\beta$-$M$ 関係
- プラントル・マイヤー膨張扇の理論と膨張角の計算
- Python による 4 つの数値計算と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 圧縮性流体とマッハ数 ── マッハ数の定義、音速、等エントロピー関係
- 連続の式 ── 質量保存の微分形・積分形
- 運動量保存則 ── 運動量の式（オイラー方程式）
- 熱力学第1法則 ── 内部エネルギーとエンタルピー
- エントロピー ── エントロピーの熱力学的定義と増大則
- 理想気体の状態方程式 ── $p = \rho R T$ の関係
- 比熱と比熱比 ── $c_v$, $c_p$, $\gamma$ の定義と関係
衝撃波はなぜ生じるのか
音速という情報伝播の壁
衝撃波の本質を理解するには、まず「気体中で情報がどのように伝わるか」を考える必要があります。
池に石を投げ込むと波紋が同心円状に広がっていきます。気体中でも同じことが起こっています。ある一点で圧力がわずかに上がると、その微小な圧力擾乱（音波）が周囲に球面状に広がっていきます。この擾乱が伝わる速度が音速 $a$ です。理想気体では、音速は次のように表されます。
$$ \begin{equation} a = \sqrt{\gamma R T} \end{equation} $$
ここで $\gamma$ は比熱比、$R$ は気体定数、$T$ は絶対温度です。空気（$\gamma \approx 1.4$, $R \approx 287 \, \mathrm{J/(kg \cdot K)}$）の場合、標準状態（$T = 288 \, \mathrm{K}$）で $a \approx 340 \, \mathrm{m/s}$ になります。
音速は気体中の「情報伝播の上限速度」です。流れの中に障害物（たとえば翼）があると、障害物は自分の存在を圧力擾乱として上流に知らせます。亜音速流では流速 $u$ が音速 $a$ より遅い（マッハ数 $M = u/a < 1$）ので、擾乱は上流に伝わることができ、流れは障害物の存在を「事前に知って」滑らかに迂回します。
超音速流で何が起きるか
ところが、流速が音速を超える（$M > 1$）と状況は一変します。流れの速度が情報の伝播速度を上回るので、上流の流体は障害物の存在を「知らない」まま高速で突っ込んでくることになります。流体には障害物を滑らかに避ける余裕がなく、急激な状態変化を強いられます。この急激な状態変化が衝撃波として現れるのです。
大雑把に言えば、衝撃波は超音速流が障害物に「出会い頭で衝突」した結果生じる現象であり、流体の圧力・密度・温度が不連続的に跳ね上がる極めて薄い領域です。衝撃波の厚さは分子の平均自由行程のオーダー（標準大気で約 $10^{-7}$ m）であり、巨視的にはゼロ厚さの不連続面として扱えます。
マッハ線とマッハ角
超音速流の中を移動する微小擾乱源を考えてみましょう。擾乱源は時刻ごとに球面波を放射しますが、擾乱源自身が音速より速く移動しているため、球面波は擾乱源に追いつけません。過去の各時刻に放射された球面波の包絡面がコーン状の面になり、これをマッハ円錐（Mach cone）、その母線をマッハ線（Mach line）と呼びます。
マッハ円錐の半頂角 $\mu$（マッハ角）は、単純な幾何学から求まります。時間 $\Delta t$ の間に、音波は距離 $a \Delta t$ だけ広がり、擾乱源は距離 $u \Delta t$ だけ進みます。直角三角形の関係から、次の式が得られます。
$$ \begin{equation} \sin \mu = \frac{a \Delta t}{u \Delta t} = \frac{a}{u} = \frac{1}{M} \end{equation} $$
したがって、
$$ \begin{equation} \mu = \arcsin\!\left(\frac{1}{M}\right) \end{equation} $$
$M = 1$ でマッハ角は $90°$（擾乱が全方向に広がる）、$M \to \infty$ でマッハ角は $0°$（円錐が極端に鋭くなる）になります。マッハ線は微小擾乱の包絡線であり、流れの状態量に不連続はありません。一方、有限の強さを持つ擾乱が集積すると衝撃波になり、状態量が不連続的に変化します。
マッハ線が「微小擾乱の限界」であるのに対し、衝撃波は「有限擾乱の集積結果」という関係を押さえておきましょう。次に、最もシンプルな衝撃波である垂直衝撃波の理論に入ります。
垂直衝撃波とランキン・ユゴニオ関係式
垂直衝撃波の問題設定
垂直衝撃波（normal shock）は、流れに対して垂直に立つ衝撃波で、衝撃波理論の最も基本的な形態です。衝撃波の上流（添字 1）と下流（添字 2）の状態量の関係を求めることが目標です。
衝撃波に固定した座標系で考えます。上流から速度 $u_1$ で超音速の気体が衝撃波に流入し、下流では速度 $u_2$ で流出します。衝撃波は定常とし、衝撃波の厚さは無視します。
衝撃波を横切る過程では粘性散逸が内部で起こるため、過程は不可逆です。しかし、衝撃波の上流面と下流面に注目すれば、粘性力や熱伝導の影響は衝撃波内部に閉じ込められており、衝撃波の外側では無視できます。そこで、衝撃波を挟んだ制御体積に対して、質量保存・運動量保存・エネルギー保存の 3 つの保存則を適用することで、上流と下流の状態量の関係が得られます。
3 つの保存則（ランキン・ユゴニオ関係式の出発点）
衝撃波を挟む一次元定常流に対して、3 つの保存則は次のように書けます。
質量保存（連続の式）:
$$ \begin{equation} \rho_1 u_1 = \rho_2 u_2 \end{equation} $$
運動量保存（運動量の式）:
$$ \begin{equation} p_1 + \rho_1 u_1^2 = p_2 + \rho_2 u_2^2 \end{equation} $$
エネルギー保存（熱力学第1法則から導かれるエネルギー方程式）:
$$ \begin{equation} h_1 + \frac{u_1^2}{2} = h_2 + \frac{u_2^2}{2} \end{equation} $$
ここで $\rho$ は密度、$u$ は速度、$p$ は圧力、$h$ は比エンタルピーです。理想気体では $h = c_p T$ であり、$c_p$ は定圧比熱です。
エネルギー方程式の左辺と右辺はそれぞれ全エンタルピー（total enthalpy）$h_0 = h + u^2/2$ に等しく、衝撃波前後で全エンタルピーが保存されることを意味しています。これは衝撃波が断熱過程（外部との熱の授受がない）であることの反映です。なお、断熱ではあっても等エントロピーではない点が重要です。
ランキン・ユゴニオ関係式の導出
上記の 3 つの保存則と理想気体の状態方程式 $p = \rho R T$ から、衝撃波前後の状態量比を上流マッハ数 $M_1$ のみの関数として表すのがランキン・ユゴニオ（Rankine-Hugoniot）関係式です。導出を一歩ずつ進めましょう。
まず、理想気体の音速 $a^2 = \gamma p / \rho$ と 比熱比 $\gamma = c_p / c_v$ を使って、比エンタルピーを次のように書き換えます。
$$ \begin{equation} h = c_p T = \frac{\gamma}{\gamma – 1} \frac{p}{\rho} \end{equation} $$
これは状態方程式 $p = \rho R T$ と $c_p = \gamma R / (\gamma – 1)$ を使えばすぐに確認できます。
エネルギー保存式に代入すると、次のようになります。
$$ \begin{equation} \frac{\gamma}{\gamma – 1} \frac{p_1}{\rho_1} + \frac{u_1^2}{2} = \frac{\gamma}{\gamma – 1} \frac{p_2}{\rho_2} + \frac{u_2^2}{2} \end{equation} $$
ここで、質量保存式から $u_2 = \rho_1 u_1 / \rho_2$ を代入できますが、最も見通しよく導出するために、マッハ数を使った表現に書き換えましょう。
$M_1 = u_1 / a_1$ および $a_1^2 = \gamma p_1 / \rho_1$ より $u_1^2 = M_1^2 \gamma p_1 / \rho_1$ です。同様に $u_2^2 = M_2^2 \gamma p_2 / \rho_2$ です。
運動量保存式を変形していきます。$\rho u^2 = \rho u^2 \cdot (\gamma p / \rho) / (\gamma p / \rho) = \gamma p M^2$ と書けるので、
$$ \begin{equation} p_1 + \gamma p_1 M_1^2 = p_2 + \gamma p_2 M_2^2 \end{equation} $$
両辺をまとめると、圧力比が得られます。
$$ \begin{equation} p_1 (1 + \gamma M_1^2) = p_2 (1 + \gamma M_2^2) \end{equation} $$
$$ \begin{equation} \frac{p_2}{p_1} = \frac{1 + \gamma M_1^2}{1 + \gamma M_2^2} \end{equation} $$
次に、エネルギー保存式をマッハ数で表現します。$h + u^2/2 = \text{const}$ を音速で書き換えると、$a^2/(\gamma – 1) + u^2/2 = \text{const}$ となります。これを上流と下流で等置すると、
$$ \begin{equation} \frac{a_1^2}{\gamma – 1} + \frac{u_1^2}{2} = \frac{a_2^2}{\gamma – 1} + \frac{u_2^2}{2} \end{equation} $$
$u = aM$ を使って $a$ で割ると、臨界音速 $a^*$（$M = 1$ での音速）を介した関係が得られますが、ここではより直接的に $M_2$ を $M_1$ で表す関係式を導きましょう。
質量保存、運動量保存、エネルギー保存を連立して整理すると（詳細な代数的操作は以下の通り）、衝撃波下流のマッハ数 $M_2$ は上流マッハ数 $M_1$ だけで決まります。
質量保存より $\rho_2/\rho_1 = u_1/u_2$ です。運動量保存を $p_2 – p_1 = \rho_1 u_1^2 – \rho_2 u_2^2 = \rho_1 u_1 (u_1 – u_2)$ と書き直します（最後の等号で質量保存 $\rho_1 u_1 = \rho_2 u_2$ を使いました）。
一方、エネルギー保存を $h_2 – h_1 = (u_1^2 – u_2^2)/2 = (u_1 + u_2)(u_1 – u_2)/2$ と書きます。ここで $h = a^2/(\gamma – 1)$ を使うと、
$$ \frac{a_2^2 – a_1^2}{\gamma – 1} = \frac{(u_1 + u_2)(u_1 – u_2)}{2} $$
また、運動量保存から $p_2 – p_1 = \rho_1 u_1 (u_1 – u_2)$ で、$p = \rho a^2/\gamma$ を使って
$$ \frac{\rho_2 a_2^2 – \rho_1 a_1^2}{\gamma} = \rho_1 u_1(u_1 – u_2) $$
これらの関係を組み合わせ、$u_1 – u_2 \neq 0$（非自明解）の条件下で整理すると、プラントル関係式（Prandtl relation）が得られます。
$$ \begin{equation} u_1 u_2 = a^{*2} \end{equation} $$
ここで $a^*$ は臨界音速です。この関係式は、衝撃波上流と下流の速度の積が臨界音速の二乗に等しいことを示しています。$u_1 > a^*$ ならば $u_2 < a^*$ でなければならず、超音速の流れは衝撃波を通過すると必ず亜音速になるという重要な結論が得られます。
衝撃波前後の状態量比
プラントル関係式と 3 つの保存則を組み合わせると、全ての状態量比を上流マッハ数 $M_1$ のみの関数として表せます。これがランキン・ユゴニオ関係式の最終形です。
下流マッハ数:
$$ \begin{equation} M_2^2 = \frac{M_1^2 + \frac{2}{\gamma – 1}}{\frac{2\gamma}{\gamma – 1} M_1^2 – 1} \end{equation} $$
圧力比:
$$ \begin{equation} \frac{p_2}{p_1} = 1 + \frac{2\gamma}{\gamma + 1}(M_1^2 – 1) \end{equation} $$
密度比（=速度比の逆数）:
$$ \begin{equation} \frac{\rho_2}{\rho_1} = \frac{u_1}{u_2} = \frac{(\gamma + 1) M_1^2}{(\gamma – 1) M_1^2 + 2} \end{equation} $$
温度比:
$$ \begin{equation} \frac{T_2}{T_1} = \frac{p_2}{p_1} \cdot \frac{\rho_1}{\rho_2} = \frac{\left[2\gamma M_1^2 – (\gamma – 1)\right]\left[(\gamma – 1) M_1^2 + 2\right]}{(\gamma + 1)^2 M_1^2} \end{equation} $$
これらの式から、いくつかの重要な物理的特徴が読み取れます。
- $M_1 = 1$ で全ての比が 1 になる: 衝撃波の強さがゼロ、すなわち擾乱のない場合に対応します
- $M_1 > 1$ で $p_2 > p_1$, $\rho_2 > \rho_1$, $T_2 > T_1$: 衝撃波を通過すると圧力・密度・温度はすべて増加します。これは衝撃波が「圧縮波」であることを意味しています
- $M_1 > 1$ で $M_2 < 1$: 超音速流は衝撃波を通過すると亜音速になります（ただし強い衝撃波の極限を除く）
- 密度比には上限がある: $M_1 \to \infty$ のとき $\rho_2/\rho_1 \to (\gamma + 1)/(\gamma – 1)$。空気（$\gamma = 1.4$）では密度比の上限は 6 です。一方、圧力比や温度比には上限がなく、$M_1 \to \infty$ で無限大に発散します
密度比に上限がある理由は直感的に理解できます。衝撃波で気体が圧縮されて密度が上がると温度も上がり、内部エネルギーが増えます。内部エネルギーの増加は圧縮に対する「抵抗」として働くため、どんなに強い衝撃波でも密度を無限に上げることはできないのです。
ここまでで、垂直衝撃波前後の状態量変化を理論的に導きました。これらの関係式が正しいことを Python で確認し、さらにグラフで視覚的に理解を深めましょう。
Python: 垂直衝撃波の状態量比
上流マッハ数 $M_1$ を横軸に、各状態量比を縦軸にプロットして、ランキン・ユゴニオ関係式の振る舞いを確認します。
import numpy as np
import matplotlib.pyplot as plt
# 比熱比
gamma = 1.4
# 上流マッハ数
M1 = np.linspace(1.0, 5.0, 500)
# ランキン・ユゴニオ関係式
# 下流マッハ数
M2_sq = (M1**2 + 2 / (gamma - 1)) / (2 * gamma / (gamma - 1) * M1**2 - 1)
M2 = np.sqrt(M2_sq)
# 圧力比
p_ratio = 1 + 2 * gamma / (gamma + 1) * (M1**2 - 1)
# 密度比
rho_ratio = (gamma + 1) * M1**2 / ((gamma - 1) * M1**2 + 2)
# 温度比
T_ratio = p_ratio / rho_ratio
# 可視化
fig, axes = plt.subplots(2, 2, figsize=(12, 9))
axes[0, 0].plot(M1, p_ratio, "b-", linewidth=2)
axes[0, 0].set_xlabel("Upstream Mach number $M_1$")
axes[0, 0].set_ylabel("$p_2 / p_1$")
axes[0, 0].set_title("Pressure ratio")
axes[0, 0].grid(True, alpha=0.3)
axes[0, 1].plot(M1, rho_ratio, "r-", linewidth=2)
axes[0, 1].axhline(
y=(gamma + 1) / (gamma - 1), color="r", linestyle="--", alpha=0.5,
label=r"$(\gamma+1)/(\gamma-1) = {:.1f}$".format((gamma + 1) / (gamma - 1)),
)
axes[0, 1].set_xlabel("Upstream Mach number $M_1$")
axes[0, 1].set_ylabel(r"$\rho_2 / \rho_1$")
axes[0, 1].set_title("Density ratio")
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
axes[1, 0].plot(M1, T_ratio, "g-", linewidth=2)
axes[1, 0].set_xlabel("Upstream Mach number $M_1$")
axes[1, 0].set_ylabel("$T_2 / T_1$")
axes[1, 0].set_title("Temperature ratio")
axes[1, 0].grid(True, alpha=0.3)
axes[1, 1].plot(M1, M2, "m-", linewidth=2)
axes[1, 1].axhline(y=1.0, color="k", linestyle="--", alpha=0.5, label="$M = 1$")
axes[1, 1].set_xlabel("Upstream Mach number $M_1$")
axes[1, 1].set_ylabel("$M_2$")
axes[1, 1].set_title("Downstream Mach number")
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
fig.suptitle("Normal Shock Relations ($\\gamma = 1.4$)", fontsize=14)
plt.tight_layout()
plt.show()
上の 4 つのグラフから、ランキン・ユゴニオ関係式の重要な特徴を視覚的に確認できます。
- 圧力比と温度比は $M_1$ の増加とともに単調に増加する: 衝撃波が強くなるほど、圧力と温度のジャンプは大きくなります。$M_1 = 5$ で圧力比は約 29、温度比は約 5.8 に達し、高マッハ数飛行での空力加熱の深刻さがわかります
- 密度比は上限 $(\gamma + 1)/(\gamma – 1) = 6$ に漸近する: 赤の破線が上限値を示しています。$M_1 = 5$ でも密度比は約 5 であり、上限にかなり近づいていることがわかります
- 下流マッハ数は常に 1 未満: $M_1$ がどんなに大きくても $M_2$ は 1 を下回り、衝撃波通過後は必ず亜音速になることが確認できます。$M_1 \to \infty$ の極限では $M_2 \to \sqrt{(\gamma – 1)/(2\gamma)} \approx 0.378$ に収束します
これらの結果は、先ほど理論的に議論した内容と完全に一致しています。次に、衝撃波がエントロピーにどのような影響を与えるかを定量的に調べましょう。
衝撃波前後のエントロピー変化
エントロピー変化の式
エントロピーは不可逆過程で増大するという性質を持ちます。衝撃波は断熱だが不可逆な過程なので、エントロピーは衝撃波を通過すると増大するはずです。これを定量的に確認しましょう。
理想気体のエントロピー変化は次のように表されます。
$$ \begin{equation} \Delta s = s_2 – s_1 = c_v \ln\frac{T_2}{T_1} – R \ln\frac{\rho_2}{\rho_1} \end{equation} $$
あるいは、圧力と温度を使って次のようにも書けます。
$$ \begin{equation} \Delta s = c_p \ln\frac{T_2}{T_1} – R \ln\frac{p_2}{p_1} \end{equation} $$
$c_p / R = \gamma/(\gamma – 1)$ の関係を使うと、無次元化されたエントロピー変化を次のように書くことができます。
$$ \begin{equation} \frac{\Delta s}{c_v} = \ln\frac{T_2}{T_1} – (\gamma – 1)\ln\frac{\rho_2}{\rho_1} \end{equation} $$
この式に先ほど導出した温度比と密度比のランキン・ユゴニオ関係式を代入すれば、エントロピー変化を $M_1$ のみの関数として表すことができます。
弱い衝撃波の近似
$M_1 = 1 + \epsilon$（$\epsilon \ll 1$）の弱い衝撃波では、エントロピー変化をテイラー展開できます。結果は次のようになります。
$$ \begin{equation} \frac{\Delta s}{c_v} \propto (M_1^2 – 1)^3 + O\!\left((M_1^2 – 1)^4\right) \end{equation} $$
エントロピー変化が $(M_1^2 – 1)$ の三乗に比例するということは、弱い衝撃波ではエントロピー変化が非常に小さいことを意味しています。これは重要な結論で、弱い衝撃波は近似的に等エントロピーとみなせるということです。マッハ線（無限に弱い衝撃波の極限）でエントロピーが変化しないことと整合しています。
一方、$M_1^2 – 1 < 0$、すなわち $M_1 < 1$ の場合、エントロピー変化は負になります。これはエントロピー増大則に反するため、亜音速流では衝撃波は物理的に存在できないことがわかります。衝撃波は超音速流でのみ発生するという事実が、熱力学の第二法則から自動的に導かれるのです。
Python: 衝撃波前後のエントロピー変化
エントロピー変化の $M_1$ 依存性を可視化して、上で述べた性質を確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
gamma = 1.4
M1 = np.linspace(1.0, 5.0, 500)
# 圧力比
p_ratio = 1 + 2 * gamma / (gamma + 1) * (M1**2 - 1)
# 密度比
rho_ratio = (gamma + 1) * M1**2 / ((gamma - 1) * M1**2 + 2)
# 温度比
T_ratio = p_ratio / rho_ratio
# エントロピー変化 (Δs / cv)
ds_cv = np.log(T_ratio) - (gamma - 1) * np.log(rho_ratio)
# エントロピー変化 (Δs / R)
ds_R = gamma / (gamma - 1) * np.log(T_ratio) - np.log(p_ratio)
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
# 左: エントロピー変化
axes[0].plot(M1, ds_cv, "b-", linewidth=2, label=r"$\Delta s / c_v$")
axes[0].plot(M1, ds_R, "r--", linewidth=2, label=r"$\Delta s / R$")
axes[0].set_xlabel("Upstream Mach number $M_1$")
axes[0].set_ylabel("Entropy change")
axes[0].set_title("Entropy change across normal shock")
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 右: 弱い衝撃波領域の拡大
M1_weak = np.linspace(1.0, 1.5, 300)
p_w = 1 + 2 * gamma / (gamma + 1) * (M1_weak**2 - 1)
rho_w = (gamma + 1) * M1_weak**2 / ((gamma - 1) * M1_weak**2 + 2)
T_w = p_w / rho_w
ds_w = np.log(T_w) - (gamma - 1) * np.log(rho_w)
# 3乗近似との比較
coeff = 2 * gamma / (gamma + 1)**2 * (gamma**2 - 1 + 3) / 3
ds_cubic = (2 * gamma * (gamma - 1)) / (3 * (gamma + 1)**2) * (M1_weak**2 - 1)**3
axes[1].plot(M1_weak, ds_w, "b-", linewidth=2, label=r"Exact $\Delta s / c_v$")
axes[1].plot(M1_weak, ds_cubic, "r--", linewidth=2, label=r"Cubic approx $\propto (M_1^2-1)^3$")
axes[1].set_xlabel("Upstream Mach number $M_1$")
axes[1].set_ylabel(r"$\Delta s / c_v$")
axes[1].set_title("Weak shock: entropy change (zoomed)")
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のグラフから、エントロピー変化は $M_1 = 1$ でゼロから出発し、$M_1$ の増大とともに単調に増加することがわかります。$\Delta s / c_v$ と $\Delta s / R$ の両方の表現でプロットしていますが、どちらも定性的な振る舞いは同じです。$M_1 = 5$ では $\Delta s / c_v \approx 1.6$ 程度に達し、無視できないエントロピー増大であることが見てとれます。
右のグラフでは、弱い衝撃波の領域を拡大して三乗近似と比較しています。$M_1 \lesssim 1.2$ の範囲では三乗近似が厳密解とよく一致しており、弱い衝撃波でのエントロピー変化が $(M_1^2 – 1)^3$ に比例するという理論的予測が確認できます。
重要な帰結として、エントロピー変化が常に正（$M_1 > 1$ のとき）であるという事実は、衝撃波が超音速から亜音速への「一方通行」であることの熱力学的証拠です。逆方向の亜音速から超音速への衝撃波はエントロピー減少を伴うため、孤立系では実現しません。
衝撃波に伴うエントロピー増大は全圧（よどみ点圧力）の低下を引き起こし、実用上は推進システムの効率低下として現れます。たとえば超音速インテークでは、衝撃波によるエントロピー増大を最小限に抑えるために、複数の弱い斜め衝撃波を使って段階的に減速する設計が採用されます。このような工学的工夫を理解するために、次に斜め衝撃波の理論に進みましょう。
斜め衝撃波
斜め衝撃波の幾何学
超音速流が楔（ウェッジ）や圧縮コーナーで方向を変えられるとき、流れの偏向角に応じた斜め衝撃波（oblique shock）が発生します。斜め衝撃波は垂直衝撃波の一般化であり、流れに対して斜めに傾いた平面上の不連続面です。
斜め衝撃波の幾何学的パラメータは次の 3 つです。
- 衝撃波角 $\beta$: 上流の流れ方向と衝撃波面のなす角
- 偏向角（転向角） $\theta$: 衝撃波通過前後の流れ方向の変化量
- 上流マッハ数 $M_1$
斜め衝撃波の解析のポイントは、速度を衝撃波面に垂直な成分と平行な成分に分解することです。衝撃波面に平行な速度成分は衝撃波を通過しても変化しません（衝撃波面に沿った方向には圧力勾配がないため）。したがって、衝撃波面に垂直な成分だけが垂直衝撃波と同じランキン・ユゴニオ関係式に従います。
上流速度の衝撃波面に垂直な成分は $u_{1n} = u_1 \sin\beta$ であり、この成分に対するマッハ数は次のようになります。
$$ \begin{equation} M_{1n} = M_1 \sin\beta \end{equation} $$
衝撃波が存在するための条件は $M_{1n} > 1$、つまり $M_1 \sin\beta > 1$ です。マッハ角 $\mu = \arcsin(1/M_1)$ を使えば、$\beta > \mu$ が衝撃波の存在条件です。
$\theta$-$\beta$-$M$ 関係
垂直成分に対して垂直衝撃波の関係式を適用し、幾何学的な関係（偏向角 $\theta$ は上流と下流の速度ベクトルの方向差）を使うと、$\theta$-$\beta$-$M$ 関係式が得られます。
下流の垂直成分に対するマッハ数 $M_{2n}$ はランキン・ユゴニオ関係式で計算でき、下流速度の衝撃波面に対する角度から偏向角が求まります。具体的には、上流では $\tan\beta = u_{1n}/u_{1t}$、下流では $\tan(\beta – \theta) = u_{2n}/u_{2t}$ です。平行成分は不変（$u_{1t} = u_{2t}$）なので、次の関係が得られます。
$$ \frac{\tan(\beta – \theta)}{\tan\beta} = \frac{u_{2n}}{u_{1n}} = \frac{\rho_1}{\rho_2} $$
ここで最後の等号は質量保存（垂直成分に対して）$\rho_1 u_{1n} = \rho_2 u_{2n}$ から来ています。密度比のランキン・ユゴニオ関係式を $M_{1n} = M_1 \sin\beta$ で書けば、整理して次の $\theta$-$\beta$-$M$ 関係式 が得られます。
$$ \begin{equation} \tan\theta = 2\cot\beta \cdot \frac{M_1^2 \sin^2\beta – 1}{M_1^2(\gamma + \cos 2\beta) + 2} \end{equation} $$
この式は、与えられた $M_1$ と $\beta$ に対して偏向角 $\theta$ を一意に与えます。逆に、$M_1$ と $\theta$ を指定すると、一般に 2 つの $\beta$ の解が存在します。
- 弱い衝撃波解（小さい $\beta$）: 超音速インテークや翼型周りの流れで通常観察される解
- 強い衝撃波解（大きい $\beta$）: 通常の外部流れでは稀だが、ダクト内流れなどで生じうる
最大偏向角と離脱衝撃波
各マッハ数 $M_1$ に対して、偏向角 $\theta$ には最大値 $\theta_{\max}$ が存在します。$\theta > \theta_{\max}$ の場合、斜め衝撃波は物体に付着した状態を維持できず、物体前方に離れた位置に弓状衝撃波（bow shock / detached shock）が形成されます。
$\theta = \theta_{\max}$ のときの衝撃波角 $\beta$ は、弱い解と強い解の境界に対応します。$\theta_{\max}$ を超えると、斜め衝撃波の解が存在しなくなり、流れ場の構造が本質的に変わるのです。鈍頭物体（たとえばスペースシャトルの先端やカプセル型宇宙機）の前方には常に離脱衝撃波が形成されます。
$\theta$-$\beta$-$M$ 関係式は非常に有用で、超音速流れの設計において中心的な役割を果たします。これを Python でプロットして、特徴を視覚的に確認しましょう。
Python: $\theta$-$\beta$-$M$ 関係図
$\theta$-$\beta$-$M$ 関係をいくつかのマッハ数について描画します。この図は超音速空気力学の「地図」とも呼ばれる重要な図です。
import numpy as np
import matplotlib.pyplot as plt
gamma = 1.4
def theta_beta_M(beta_deg, M1):
"""θ-β-M関係式: 偏向角θ(度)を返す"""
beta = np.radians(beta_deg)
sin_b = np.sin(beta)
cos_b = np.cos(beta)
tan_b = np.tan(beta)
numer = M1**2 * sin_b**2 - 1
denom = M1**2 * (gamma + np.cos(2 * beta)) + 2
# 物理的に意味のある範囲のみ
tan_theta = 2 * cos_b / sin_b * numer / denom
theta = np.degrees(np.arctan(tan_theta))
return theta
fig, ax = plt.subplots(figsize=(10, 7))
Mach_numbers = [1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 10.0]
colors = plt.cm.viridis(np.linspace(0, 0.9, len(Mach_numbers)))
for M1, color in zip(Mach_numbers, colors):
mu = np.degrees(np.arcsin(1 / M1)) # マッハ角
beta_arr = np.linspace(mu + 0.1, 89.9, 2000)
theta_arr = theta_beta_M(beta_arr, M1)
# θが正の部分のみプロット
mask = theta_arr > 0
ax.plot(
theta_arr[mask], beta_arr[mask],
color=color, linewidth=1.8, label=f"$M_1 = {M1}$",
)
# 最大偏向角の位置にマーカー
if np.any(mask):
idx_max = np.argmax(theta_arr[mask])
ax.plot(
theta_arr[mask][idx_max], beta_arr[mask][idx_max],
"o", color=color, markersize=5,
)
# M2 = 1 の境界（強い解と弱い解の境界の参考）
ax.set_xlabel("Deflection angle $\\theta$ [deg]", fontsize=12)
ax.set_ylabel("Shock wave angle $\\beta$ [deg]", fontsize=12)
ax.set_title("$\\theta$-$\\beta$-$M$ diagram ($\\gamma = 1.4$)", fontsize=14)
ax.set_xlim(0, 50)
ax.set_ylim(0, 90)
ax.legend(fontsize=10, loc="upper left")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
この $\theta$-$\beta$-$M$ 図から、斜め衝撃波の重要な特徴が読み取れます。
- 各マッハ数の曲線は逆 U 字型を描く: 偏向角 $\theta$ はある $\beta$ で最大値をとり、それ以上の $\beta$ では $\theta$ が減少して $\beta = 90°$（垂直衝撃波）で $\theta = 0$ に戻ります。丸いマーカーが最大偏向角 $\theta_{\max}$ の位置を示しています
- 同じ $\theta$ に対して 2 つの $\beta$ が存在する: 曲線の上側が強い衝撃波解、下側が弱い衝撃波解です。自然界では通常、弱い衝撃波解が選ばれます
- マッハ数が大きいほど $\theta_{\max}$ が大きい: 高マッハ数の流れほど大きな偏向角に耐えられます。$M_1 = 1.5$ では $\theta_{\max} \approx 12°$ 程度ですが、$M_1 = 10$ では $\theta_{\max} \approx 45°$ に近づきます
- $\beta \to \mu$（マッハ角）の極限で $\theta \to 0$: 弱い衝撃波解の下限はマッハ線であり、マッハ線は無限に弱い衝撃波と解釈できます
この図は超音速インテークの設計で極めて重要です。たとえば、インテークのランプ角（偏向角）を決めれば、飛行マッハ数に対応する衝撃波角が図から直ちに読み取れます。
斜め衝撃波は流れを圧縮して方向を変える波でした。では逆に、超音速流が膨張する（方向が広がる）場合には何が起きるのでしょうか。次に、衝撃波の対概念である膨張波について解説します。
プラントル・マイヤー膨張扇
膨張波の物理
超音速流が凸コーナー（外向きコーナー）で方向を変えるとき、流れは膨張して加速されます。この膨張過程で生じるのがプラントル・マイヤー膨張扇（Prandtl-Meyer expansion fan）です。
膨張波は衝撃波とは本質的に異なる性質を持っています。
| 性質 | 衝撃波（圧縮） | 膨張波 |
|---|---|---|
| 不連続面 | あり（有限の状態量ジャンプ） | なし（連続的に変化） |
| エントロピー | 増大（不可逆） | 一定（等エントロピー） |
| マッハ数 | 減少（超音速→亜音速） | 増大（超音速→さらに超音速） |
| 圧力 | 増大 | 減少 |
| 温度 | 増大 | 減少 |
衝撃波が「圧縮波の集積」で不連続面を形成するのに対し、膨張波は多数のマッハ線（微小擾乱波）が扇状に広がった形をしています。各マッハ線を横切るたびに、流れはわずかに加速・膨張・偏向します。マッハ線は微小擾乱なので等エントロピーであり、膨張扇全体も等エントロピー過程です。
なぜ膨張波が「扇」状になるかを直感的に理解しましょう。膨張が進むにつれてマッハ数が増大し、マッハ角 $\mu = \arcsin(1/M)$ が減少します。そのため、膨張扇の先頭のマッハ線（膨張が始まる前、マッハ数が最も小さい）は角度が大きく、末尾のマッハ線（膨張が終わった後、マッハ数が最も大きい）は角度が小さくなります。結果として、マッハ線は扇の骨のように広がるのです。
プラントル・マイヤー関数
膨張扇を通過する前後の流れの偏向角とマッハ数の関係は、プラントル・マイヤー関数 $\nu(M)$ で記述されます。この関数は、等エントロピーの仮定のもとで運動量とエネルギーの保存から導かれます。
微小偏向角 $d\theta$ とマッハ数の変化 $dM$ の関係は、等エントロピー流れの微小擾乱解析から次のように得られます。
$$ d\theta = \sqrt{M^2 – 1} \cdot \frac{dM/M}{1 + \frac{\gamma – 1}{2}M^2} $$
この式の右辺は $M$ のみの関数なので、$M = 1$（$\theta = 0$ と定義）から任意の $M$ まで積分できます。
$$ \begin{equation} \nu(M) = \int_0^\theta d\theta = \int_1^M \frac{\sqrt{M’^2 – 1}}{1 + \frac{\gamma – 1}{2}M’^2} \frac{dM’}{M’} \end{equation} $$
この積分は解析的に実行でき、置換 $M’^2 – 1 = t^2$ を使って計算すると、次の閉じた形になります。
$$ \begin{equation} \nu(M) = \sqrt{\frac{\gamma + 1}{\gamma – 1}} \arctan\!\sqrt{\frac{\gamma – 1}{\gamma + 1}(M^2 – 1)} – \arctan\!\sqrt{M^2 – 1} \end{equation} $$
プラントル・マイヤー関数 $\nu(M)$ の物理的意味は、「マッハ数 $M$ の超音速流を $M = 1$ の音速流まで等エントロピー的に圧縮するのに必要な偏向角」です。あるいは逆に、「音速流を等エントロピー的に膨張させてマッハ数 $M$ にするのに必要な偏向角」とも解釈できます。
膨張扇の計算方法
マッハ数 $M_1$ の超音速流が偏向角 $\theta$ の凸コーナーで膨張するとき、下流のマッハ数 $M_2$ は次の関係式で求まります。
$$ \begin{equation} \nu(M_2) = \nu(M_1) + \theta \end{equation} $$
右辺の $\nu(M_1) + \theta$ を計算し、プラントル・マイヤー関数の逆関数を使って $M_2$ を求めるという手順です。プラントル・マイヤー関数は単調増加関数なので、逆関数は一意に存在します。
プラントル・マイヤー関数の最大値は $M \to \infty$ のときで、次のようになります。
$$ \begin{equation} \nu_{\max} = \frac{\pi}{2}\left(\sqrt{\frac{\gamma + 1}{\gamma – 1}} – 1\right) \end{equation} $$
空気（$\gamma = 1.4$）では $\nu_{\max} \approx 130.45°$ です。これは、流体を完全真空（$p = 0$, $T = 0$, $M \to \infty$）まで膨張させるのに必要な最大偏向角に対応します。実際にはこのような極限的な膨張は実現しませんが、理論上の上限として重要な値です。
膨張扇は等エントロピー過程なので、圧縮性流体の等エントロピー関係をそのまま使って膨張扇前後の圧力比・温度比・密度比を計算できます。
$$ \frac{p_2}{p_1} = \left(\frac{1 + \frac{\gamma-1}{2}M_1^2}{1 + \frac{\gamma-1}{2}M_2^2}\right)^{\frac{\gamma}{\gamma-1}}, \quad \frac{T_2}{T_1} = \frac{1 + \frac{\gamma-1}{2}M_1^2}{1 + \frac{\gamma-1}{2}M_2^2} $$
衝撃波と違い、膨張扇では全圧が保存されます。これが超音速風洞のノズル設計でプラントル・マイヤー膨張理論が用いられる理由です。ノズル壁面の形状を膨張扇の理論に基づいて設計すれば、損失なく所望のマッハ数を実現できるのです。
それでは、プラントル・マイヤー関数を Python で計算し、その性質を確認しましょう。
Python: プラントル・マイヤー関数
プラントル・マイヤー関数 $\nu(M)$ をプロットし、さらに膨張扇の計算例を示します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq
gamma = 1.4
def prandtl_meyer(M):
"""プラントル・マイヤー関数 ν(M) [rad]"""
g = gamma
return (
np.sqrt((g + 1) / (g - 1))
* np.arctan(np.sqrt((g - 1) / (g + 1) * (M**2 - 1)))
- np.arctan(np.sqrt(M**2 - 1))
)
def inverse_prandtl_meyer(nu_target):
"""プラントル・マイヤー関数の逆関数: ν → M"""
func = lambda M: prandtl_meyer(M) - nu_target
return brentq(func, 1.0, 50.0)
# --- プラントル・マイヤー関数のプロット ---
M_arr = np.linspace(1.0, 20.0, 1000)
nu_arr = np.degrees(prandtl_meyer(M_arr))
# 理論的最大値
nu_max = np.degrees(
np.pi / 2 * (np.sqrt((gamma + 1) / (gamma - 1)) - 1)
)
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
axes[0].plot(M_arr, nu_arr, "b-", linewidth=2)
axes[0].axhline(y=nu_max, color="r", linestyle="--", alpha=0.6,
label=f"$\\nu_{{\\max}} = {nu_max:.1f}°$")
axes[0].set_xlabel("Mach number $M$", fontsize=12)
axes[0].set_ylabel("$\\nu(M)$ [deg]", fontsize=12)
axes[0].set_title("Prandtl-Meyer function", fontsize=13)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
# --- 膨張扇の計算例 ---
M1_example = 2.0
theta_expansion = np.linspace(0, 40, 200) # 膨張角 [deg]
M2_result = []
for th in theta_expansion:
nu1 = prandtl_meyer(M1_example)
nu2_target = nu1 + np.radians(th)
if nu2_target < np.pi / 2 * (np.sqrt((gamma + 1) / (gamma - 1)) - 1):
M2_result.append(inverse_prandtl_meyer(nu2_target))
else:
M2_result.append(np.nan)
M2_result = np.array(M2_result)
# 膨張扇前後の圧力比（等エントロピー）
p_ratio_exp = (
(1 + (gamma - 1) / 2 * M1_example**2)
/ (1 + (gamma - 1) / 2 * M2_result**2)
) ** (gamma / (gamma - 1))
ax2 = axes[1]
ax2.plot(theta_expansion, M2_result, "b-", linewidth=2, label="$M_2$")
ax2_twin = ax2.twinx()
ax2_twin.plot(theta_expansion, p_ratio_exp, "r--", linewidth=2, label="$p_2/p_1$")
ax2.set_xlabel("Expansion angle $\\theta$ [deg]", fontsize=12)
ax2.set_ylabel("Downstream Mach number $M_2$", fontsize=12, color="b")
ax2_twin.set_ylabel("Pressure ratio $p_2/p_1$", fontsize=12, color="r")
ax2.set_title(f"Expansion fan: $M_1 = {M1_example}$", fontsize=13)
ax2.tick_params(axis="y", labelcolor="b")
ax2_twin.tick_params(axis="y", labelcolor="r")
lines1, labels1 = ax2.get_legend_handles_labels()
lines2, labels2 = ax2_twin.get_legend_handles_labels()
ax2.legend(lines1 + lines2, labels1 + labels2, fontsize=11)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
左のグラフから、プラントル・マイヤー関数 $\nu(M)$ は $M = 1$ でゼロから出発し、$M$ の増大とともに単調に増加して、理論的最大値 $\nu_{\max} \approx 130.5°$（赤い破線）に漸近することがわかります。$M = 2$ で $\nu \approx 26.4°$、$M = 5$ で $\nu \approx 76.9°$ 程度です。関数の増加率は低マッハ数領域で大きく、高マッハ数領域で小さくなります。これは、高マッハ数ほどマッハ数を 1 だけ増やすのに必要な偏向角が小さくなることを意味しています。
右のグラフは、$M_1 = 2.0$ の超音速流が膨張扇を通過する場合の、膨張角 $\theta$ に対する下流マッハ数 $M_2$ と圧力比 $p_2/p_1$ の変化です。膨張角が大きくなるほどマッハ数は増大し、圧力比は急激に低下します。たとえば $\theta = 20°$ で $M_2 \approx 2.83$, $p_2/p_1 \approx 0.26$ となり、圧力が約 4 分の 1 に下がります。
これらの結果から、プラントル・マイヤー膨張が流れを効率よく加速する手段であることが実感できます。ラバルノズルの超音速部分はまさにこの原理を利用して流れを加速しているのです。
衝撃波と膨張波の対比
ここまでの内容を整理して、衝撃波（圧縮波）と膨張波の違いを体系的にまとめましょう。
基本的な対称性
超音速流は壁面の凹コーナー（流れが内側に曲がる）で衝撃波を、壁面の凸コーナー（流れが外側に曲がる）で膨張扇を発生させます。これは超音速空気力学における最も基本的な原理の一つであり、複雑な流れ場を理解する出発点です。
超音速翼型の周りの流れを考えると、この原理がよくわかります。翼の前縁では、上面と下面にそれぞれ斜め衝撃波が生じ、流れが圧縮されて翼面に沿います。翼の最大厚さの位置では、壁面の傾斜が変わり膨張扇が発生します。翼の後縁では再び衝撃波が生じます。このように、衝撃波と膨張扇が交互に現れて超音速翼型の圧力分布を決めるのです。
エントロピー生成の観点
工学的に最も重要な違いはエントロピー生成の有無です。衝撃波ではエントロピーが増大するため、全圧（よどみ点圧力）が低下します。全圧の低下は推進システムの効率損失に直結します。
たとえば超音速インテークの設計では、できるだけ弱い衝撃波を使って段階的に流れを減速させることで、全圧回復率（下流全圧と上流全圧の比）を最大化します。1 本の強い垂直衝撃波で一気に減速するよりも、複数の弱い斜め衝撃波で段階的に減速する方が全圧損失は小さくなります。これは弱い衝撃波のエントロピー変化が $(M_{1n}^2 – 1)^3$ に比例するという三乗則の帰結です。
一方、膨張扇は等エントロピーなので全圧の損失がなく、ラバルノズルやロケットノズルの設計に理想的な膨張過程を提供します。
数学的構造の違い
衝撃波は双曲型偏微分方程式（ナビエ・ストークス方程式のオイラー方程式近似）の弱解として数学的に記述されます。保存則を弱い意味で満たす不連続解は一般に一意ではなく、エントロピー条件（衝撃波通過でエントロピーが増大する）という追加条件で物理的に正しい解を選択します。
膨張扇は同じ双曲型方程式のリーマン不変量に沿った滑らかな解であり、特性曲線法で自然に構成できます。プラントル・マイヤー膨張は中心化された単純波（centered simple wave）の代表例です。
数学的な観点からは、衝撃波は「特性曲線が交差して不連続面が必要になる」状況に対応し、膨張扇は「特性曲線が広がって滑らかな解が構成できる」状況に対応しています。圧縮波では特性曲線が収束して衝突するため、有限時間で解が破綻し（多価になり）、衝撃波という不連続解の導入が必要になるのです。
ここまでの理論をまとめる意味で、最後に垂直衝撃波の状態量を一括計算する Python プログラムを作成しましょう。
Python: 垂直衝撃波の状態量一括計算
任意の上流マッハ数 $M_1$ を入力すると、ランキン・ユゴニオ関係式から全ての状態量比を計算し、加えてよどみ点圧力比（全圧比）とエントロピー変化も出力する関数を作成します。
import numpy as np
gamma = 1.4
def normal_shock(M1):
"""
垂直衝撃波の状態量比を計算する
入力: 上流マッハ数 M1 (M1 >= 1)
"""
g = gamma
M1sq = M1**2
# 下流マッハ数
M2 = np.sqrt((M1sq + 2 / (g - 1)) / (2 * g / (g - 1) * M1sq - 1))
# 圧力比 p2/p1
p_ratio = 1 + 2 * g / (g + 1) * (M1sq - 1)
# 密度比 rho2/rho1
rho_ratio = (g + 1) * M1sq / ((g - 1) * M1sq + 2)
# 温度比 T2/T1
T_ratio = p_ratio / rho_ratio
# 全圧比 p02/p01
p0_ratio = (
rho_ratio ** (g / (g - 1))
* (1 / p_ratio) ** (1 / (g - 1))
* ((2 * g * M1sq - (g - 1)) / (g + 1)) ** (-1 / (g - 1))
)
# 等エントロピー関係を使った全圧比の別表現
p0_ratio2 = (
((g + 1) * M1sq / ((g - 1) * M1sq + 2)) ** (g / (g - 1))
/ (1 + 2 * g / (g + 1) * (M1sq - 1)) ** (1 / (g - 1))
)
# エントロピー変化 Δs/R
ds_R = -np.log(p0_ratio2)
return {
"M1": M1,
"M2": M2,
"p2/p1": p_ratio,
"rho2/rho1": rho_ratio,
"T2/T1": T_ratio,
"p02/p01": p0_ratio2,
"ds/R": ds_R,
}
# いくつかの代表的なマッハ数で計算
print(f"{'M1':>5s} {'M2':>7s} {'p2/p1':>8s} {'rho2/rho1':>10s} "
f"{'T2/T1':>8s} {'p02/p01':>8s} {'ds/R':>8s}")
print("-" * 60)
for M1 in [1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 10.0]:
result = normal_shock(M1)
print(
f"{result['M1']:5.1f} {result['M2']:7.4f} {result['p2/p1']:8.4f} "
f"{result['rho2/rho1']:10.4f} {result['T2/T1']:8.4f} "
f"{result['p02/p01']:8.4f} {result['ds/R']:8.4f}"
)
この計算結果の表から、いくつかの重要な数値的特徴が確認できます。
- $M_1 = 1.0$ ですべての比が 1、エントロピー変化がゼロ: 衝撃波の強さがゼロの極限に対応し、理論と整合しています
- $M_1 = 2.0$ で $p_2/p_1 = 4.5$, $T_2/T_1 \approx 1.69$: マッハ 2 の衝撃波でも圧力が 4.5 倍に跳ね上がり、温度も約 70% 上昇します。航空機の構造設計において衝撃波荷重を考慮する必要性がわかります
- 全圧比 $p_{02}/p_{01}$ はマッハ数の増大とともに急速に低下: $M_1 = 2.0$ で約 0.72、$M_1 = 5.0$ で約 0.06 です。マッハ 5 の衝撃波では全圧の 94% が失われることになり、スクラムジェットエンジンの設計がいかに困難であるかがわかります
- $M_1 = 10.0$ で密度比が 5.71: 理論上の上限 $(\gamma+1)/(\gamma-1) = 6$ に非常に近い値です
このコードは超音速空気力学の基本ツールとして活用でき、たとえばインテーク設計や風洞実験の事前予測に使えます。
応用: 超音速インテークと衝撃波系
ここまでの理論の実用的な応用例として、超音速インテーク（吸気口）の設計思想を簡単に紹介します。
超音速航空機のジェットエンジン（ラムジェットやターボジェット）は、入口の空気を亜音速まで減速させてから燃焼器に導く必要があります。この減速を行うのがインテークの役割です。
最も単純なインテークは、入口面に垂直衝撃波を立てて一気にマッハ数を 1 以下に下げるピトー式インテークです。しかし、前述のように強い垂直衝撃波は大きな全圧損失を伴います。
そこで、多段圧縮インテークが採用されます。ランプ（傾斜面）を複数段設け、各段で弱い斜め衝撃波を発生させて段階的に減速した後、最後に弱い垂直衝撃波で亜音速にするという設計です。SR-71 偵察機のインテークはこの方式の代表例であり、マッハ 3.2 で飛行しながら高い全圧回復率を実現していました。
この設計最適化には、$\theta$-$\beta$-$M$ 関係図とランキン・ユゴニオ関係式が不可欠であり、本記事で学んだ理論がそのまま実用に結びつきます。
Python: 衝撃波とマッハ角の可視化
最後に、超音速流中の物体（楔型）周りの衝撃波とマッハ線を概念的に可視化するコードを示します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq
gamma = 1.4
def theta_from_beta_M(beta_rad, M1):
"""β, M1 から θ を計算"""
sin_b = np.sin(beta_rad)
cos_b = np.cos(beta_rad)
numer = M1**2 * sin_b**2 - 1
denom = M1**2 * (gamma + np.cos(2 * beta_rad)) + 2
tan_theta = 2 * cos_b / sin_b * numer / denom
return np.arctan(tan_theta)
def find_beta_weak(M1, theta_rad):
"""弱い衝撃波の β を求める"""
mu = np.arcsin(1 / M1)
func = lambda b: theta_from_beta_M(b, M1) - theta_rad
return brentq(func, mu + 1e-6, np.pi / 2 - 1e-6)
# パラメータ
M1 = 3.0
theta_deg = 15.0
theta_rad = np.radians(theta_deg)
mu_deg = np.degrees(np.arcsin(1 / M1))
# 衝撃波角を求める
beta_rad = find_beta_weak(M1, theta_rad)
beta_deg = np.degrees(beta_rad)
fig, ax = plt.subplots(figsize=(12, 6))
# 楔型の描画
wedge_length = 4.0
wedge_x = [0, wedge_length * np.cos(theta_rad), wedge_length]
wedge_y_upper = [0, -wedge_length * np.sin(theta_rad), -wedge_length * np.sin(theta_rad)]
wedge_y_lower = [0, wedge_length * np.sin(theta_rad), wedge_length * np.sin(theta_rad)]
ax.fill_between(
[0, wedge_length * np.cos(theta_rad)],
[0, -wedge_length * np.sin(theta_rad)],
[0, wedge_length * np.sin(theta_rad)],
color="gray", alpha=0.5, label="Wedge",
)
# 斜め衝撃波（上下対称）
shock_len = 3.5
shock_x_upper = [-shock_len * np.cos(beta_rad), 0]
shock_y_upper = [shock_len * np.sin(beta_rad), 0]
shock_x_lower = [-shock_len * np.cos(beta_rad), 0]
shock_y_lower = [-shock_len * np.sin(beta_rad), 0]
ax.plot(shock_x_upper, shock_y_upper, "r-", linewidth=2.5,
label=f"Oblique shock ($\\beta = {beta_deg:.1f}°$)")
ax.plot(shock_x_lower, shock_y_lower, "r-", linewidth=2.5)
# マッハ線（上流のマッハ角で参考表示）
mach_len = 3.0
mu_rad = np.arcsin(1 / M1)
for y_offset in [1.5, 2.5]:
x_start = -5.0
ax.plot(
[x_start, x_start + mach_len * np.cos(mu_rad)],
[y_offset, y_offset - mach_len * np.sin(mu_rad)],
"b--", linewidth=0.8, alpha=0.5,
)
ax.plot(
[x_start, x_start + mach_len * np.cos(mu_rad)],
[-y_offset, -y_offset + mach_len * np.sin(mu_rad)],
"b--", linewidth=0.8, alpha=0.5,
)
# 流れの矢印
for y_pos in [-2.0, -1.0, 0.0, 1.0, 2.0]:
ax.annotate(
"", xy=(-3.0, y_pos), xytext=(-5.0, y_pos),
arrowprops=dict(arrowstyle="->", color="navy", lw=1.5),
)
# テキスト
ax.text(-5.5, 2.8, f"$M_1 = {M1}$", fontsize=14, color="navy")
ax.text(-5.5, 2.2, f"Mach angle $\\mu = {mu_deg:.1f}°$", fontsize=11, color="blue")
ax.text(-4.0, -2.8, f"Wedge half-angle $\\theta = {theta_deg:.0f}°$",
fontsize=11, color="gray")
ax.set_xlim(-6, 5)
ax.set_ylim(-3.5, 3.5)
ax.set_aspect("equal")
ax.set_xlabel("$x$", fontsize=12)
ax.set_ylabel("$y$", fontsize=12)
ax.set_title(
f"Oblique shock around a wedge: $M_1={M1}$, $\\theta={theta_deg:.0f}°$, "
f"$\\beta={beta_deg:.1f}°$",
fontsize=13,
)
ax.legend(loc="upper right", fontsize=11)
ax.grid(True, alpha=0.2)
plt.tight_layout()
plt.show()
この図から、超音速流が楔型物体に衝突したときの流れ場の構造が視覚的に理解できます。
- 衝撃波角 $\beta$ は偏向角 $\theta$ より大きい: $M_1 = 3.0$, $\theta = 15°$ に対して $\beta \approx 32.2°$ です。衝撃波は楔面よりも急な角度で立っています
- マッハ角 $\mu$ は衝撃波角 $\beta$ より小さい: $\mu \approx 19.5°$ であり、マッハ線は衝撃波よりも浅い角度です。これは、マッハ線が無限に弱い衝撃波の極限であることと整合しています
- 上流の流れは衝撃波に到達するまで一様: 超音速流では擾乱が上流に伝播できないため、衝撃波の上流側の流れは完全に乱されていません
このような衝撃波の幾何学的構造は、超音速機の空力設計の基礎になっています。機首の先端角度やインテークのランプ角度を変えると、衝撃波の角度と強さが変わり、抗力や圧力分布が大きく変化します。
まとめ
本記事では、衝撃波と膨張波の物理的原理から数学的記述、Python による数値計算まで、体系的に解説しました。
-
衝撃波の物理的原因: 超音速流では流速が音速（情報伝播速度）を上回るため、流体は障害物の存在を事前に「知る」ことができず、急激な状態変化（衝撃波）を強いられます。マッハ線とマッハ角 $\mu = \arcsin(1/M)$ が微小擾乱の包絡線であり、有限擾乱の集積が衝撃波を形成します
-
ランキン・ユゴニオ関係式: 質量・運動量・エネルギーの保存則から、垂直衝撃波前後の全ての状態量比が上流マッハ数 $M_1$ のみの関数として表されます。衝撃波通過で圧力・密度・温度は増大し、マッハ数は減少します。密度比には上限 $(\gamma+1)/(\gamma-1)$ があります
-
エントロピー増大: 衝撃波は不可逆過程であり、エントロピーは常に増大します。エントロピー変化は弱い衝撃波で $(M_1^2 – 1)^3$ に比例し、亜音速での衝撃波はエントロピー減少を伴うため物理的に禁止されます
-
斜め衝撃波: $\theta$-$\beta$-$M$ 関係式により、流れの偏向角と衝撃波角、上流マッハ数が結びつけられます。各マッハ数に最大偏向角が存在し、それを超えると離脱衝撃波が形成されます
-
プラントル・マイヤー膨張扇: 超音速流の膨張は等エントロピー的に進行し、プラントル・マイヤー関数 $\nu(M)$ で記述されます。衝撃波と対照的に、全圧損失がなく効率的な加速が可能です
-
工学的応用: 超音速インテークの多段圧縮設計、ラバルノズルの超音速部設計、再突入体の空力加熱予測など、衝撃波と膨張波の理論は超音速・極超音速の航空宇宙工学に不可欠です
衝撃波の理論は、圧縮性流体力学のなかでも特にエレガントで実用的な体系です。本記事で扱った一次元・二次元の理論を基盤として、三次元の衝撃波干渉、衝撃波/境界層干渉、非定常衝撃波、極超音速域でのリアルガス効果など、より高度なトピックへと発展していきます。
次のステップとして、以下の記事も参考にしてください。
- 圧縮性流体とマッハ数 ── 衝撃波理論の前提となる等エントロピー関係の詳細
- ラバルノズル ── 膨張波理論を応用した超音速ノズルの設計
- ベルヌーイの定理 ── 非圧縮性流体における類似の保存則（衝撃波理論との違いを理解する参考に）
- ナビエ・ストークス方程式 ── 衝撃波内部の構造を記述するために必要な粘性流体の基礎方程式