ロケットが地球の重力を振り切って宇宙へ飛び立つとき、エンジンの噴射ガスは音速をはるかに超える速度で噴き出しています。スペースXのマーリンエンジンでは排気速度がマッハ3を超え、スペースシャトルのメインエンジン（SSME）ではマッハ5にも達していました。ここで素朴な疑問が浮かびます――ノズルの形をどのように設計すれば、ガスをこれほどの超音速まで加速できるのでしょうか。
日常の感覚では、ホースの先を指でつまんで断面積を狭くすると水が勢いよく飛び出します。流路を絞れば流速が上がる、というベルヌーイの定理の直感です。しかし、この「絞れば速くなる」という常識は、実は音速を超えると正反対になります。超音速の流れでは、むしろ流路を広げなければ加速できないのです。
この驚くべき物理を巧みに利用したのがラバルノズル（Laval nozzle、de Laval nozzle）です。先細（収縮）部で亜音速のガスを音速まで加速し、スロート（最小断面積部）で音速に到達させ、その後の末広（拡大）部で超音速まで一気に加速する――この独特な形状が、圧縮性流体力学の基本原理から自然に導かれます。
ラバルノズルの原理を理解すると、以下のような幅広い応用分野が見えてきます。
- ロケット推進: ロケットノズルの設計においてラバルノズルは不可欠であり、比推力の最大化に直結します
- ジェットエンジン: ターボジェットやターボファンエンジンの排気ノズルにラバルノズルが使われています
- 超音速風洞: 航空機や弾道体の空力試験に用いる風洞では、試験部のマッハ数をラバルノズルの面積比で制御します
- 蒸気タービン: スウェーデンの技術者ド・ラバルが最初にこのノズルを発明したのは、蒸気タービンの効率向上のためでした
- 産業応用: 超音速ジェット洗浄、ロケットの噴射による消火、超音速スプレーコーティングなど
本記事では、ノズル内の流れを支配する面積-速度関係式を連続の式と圧縮性流体の基本式から丁寧に導出し、スロートでの臨界条件（マッハ数 = 1）がなぜ成り立つのかを明らかにします。さらに、チョーク流れと最大質量流量、ノズル内の圧力分布パターン（完全膨張・過膨張・不足膨張）、そしてロケットノズルとの関連を解説します。最後に、Pythonによる4つの数値計算・可視化を通じて、理論の理解を定着させます。
本記事の内容
- 先細ノズルとラバルノズルの違い
- 面積-速度関係式（$dA/A = (M^2 – 1)\,dv/v$）の導出
- スロートでの臨界条件（$M = 1$）の物理的意味
- チョーク流れと最大質量流量
- ノズル内の圧力分布パターン（完全膨張・過膨張・不足膨張）
- ロケットノズルとの関連と推力係数
- Pythonによる4つの数値計算・可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 圧縮性流体とマッハ数 — マッハ数の定義、等エントロピー関係式、よどみ点状態量の概念
- 連続の式 — 質量保存から導かれる流体の基本拘束条件
- ベルヌーイの定理 — 流体のエネルギー保存則（非圧縮性の場合）
- エントロピー — 等エントロピー過程の意味と条件
- 理想気体の状態方程式 — $pV = nRT$ と理想気体の基本性質
先細ノズルとラバルノズル
先細ノズル（収束ノズル）
まずは最もシンプルなノズルである先細ノズル（convergent nozzle）から考えましょう。これは入口から出口に向かって断面積が徐々に小さくなるだけの形状です。タンクに高圧ガスを溜めて出口から噴き出す場面を想像してください。
非圧縮性の流体（水など）であれば、連続の式 $\rho A v = \text{const}$ から、断面積 $A$ が小さくなれば流速 $v$ は単調に増加します。ホースの先を指でつまむと水が速く出る、というあの現象です。
しかし、空気のような圧縮性流体では話が複雑になります。密度 $\rho$ も変化するため、断面積を小さくしても流速が単調に増えるとは限りません。先細ノズルでは、出口で到達できるマッハ数の上限は $M = 1$（音速）です。タンクの圧力をいくら上げても、出口のマッハ数は1を超えることができません。これが後述するチョーク（choking）現象です。
では、音速を超えてさらに加速するにはどうすればよいのでしょうか。ここで登場するのがラバルノズルです。
ラバルノズル（収束-発散ノズル）
ラバルノズル（convergent-divergent nozzle, C-D nozzle）は、先細部（収束部）の後にスロート（最小断面積部）を経て、末広部（発散部）が続く形状をしています。1890年代にスウェーデンの技術者カール・グスタフ・パトリック・ド・ラバル（Carl Gustaf Patrik de Laval）が蒸気タービン用に発明しました。
ラバルノズルの動作を直感的に理解するために、3つの領域に分けて考えましょう。
- 先細部（収束部）: 亜音速の流れが加速される領域です。断面積が減少するにつれて流速が増加し、圧力と密度は低下します。ここは日常の直感に合います。
- スロート: 断面積が最小になる位置で、流れがちょうど音速（$M = 1$）に達します。後で示すように、面積-速度関係式からスロートでは必ず $M = 1$ となります。
- 末広部（発散部）: 超音速の流れがさらに加速される領域です。断面積が増加するにつれて流速がさらに増加します。これは日常の直感とは正反対であり、圧縮性効果ならではの現象です。
なぜ「超音速では広がると加速する」のでしょうか。この直感に反する振る舞いの根拠を数学的に理解するために、次に面積-速度関係式を導出します。
面積-速度関係式の導出
出発点: 3つの基本方程式
面積-速度関係式を導くために、以下の3つの基本方程式を用います。いずれも定常・一次元の等エントロピー流れを仮定しています。
1. 連続の式（質量保存）
連続の式は、定常流では流路に沿って質量流量が一定であることを表します。
$$ \begin{equation} \rho A v = \text{const} \end{equation} $$
ここで $\rho$ は密度、$A$ は流路断面積、$v$ は流速です。この両辺の微分（対数微分）を取ると、次のようになります。
$$ \begin{equation} \frac{d\rho}{\rho} + \frac{dA}{A} + \frac{dv}{v} = 0 \end{equation} $$
この式は、密度・断面積・流速のそれぞれの変化率の和がゼロであること、すなわち一方が増えれば他方が減って辻褄が合うことを意味しています。
2. オイラーの運動方程式（運動量保存）
摩擦のない定常1次元流れに対するオイラーの運動方程式は次の通りです。
$$ \begin{equation} dp + \rho v \, dv = 0 \end{equation} $$
これはベルヌーイの定理の微分形に相当し、圧力の増加が流速の減少を引き起こすことを表しています。
3. 音速の定義
圧縮性流体の理論から、音速 $a$ は次のように定義されます。
$$ \begin{equation} a^2 = \frac{dp}{d\rho} \end{equation} $$
等エントロピー過程では $a^2 = \gamma p / \rho = \gamma R T$ です（$\gamma$ は比熱比、$R$ はガス定数、$T$ は温度）。この式を変形すると、密度変化と圧力変化の関係が得られます。
$$ \begin{equation} dp = a^2 \, d\rho \end{equation} $$
導出のステップ
ここからが面積-速度関係式の導出の核心です。ゴールは、断面積の変化 $dA$ と流速の変化 $dv$ の関係を、マッハ数 $M = v/a$ を使って表すことです。
まず、オイラーの運動方程式 (3) から $dp$ を消去します。式 (3) より次が得られます。
$$ dp = -\rho v \, dv $$
この右辺を音速の定義式 (5) と等しいとおくことで、$d\rho$ と $dv$ の関係が得られます。
$$ a^2 \, d\rho = -\rho v \, dv $$
両辺を $a^2 \rho$ で割ると、密度の変化率が流速の変化率で表されます。
$$ \begin{equation} \frac{d\rho}{\rho} = -\frac{v^2}{a^2} \frac{dv}{v} = -M^2 \frac{dv}{v} \end{equation} $$
ここでマッハ数の定義 $M = v/a$ を用いました。この式は重要な物理的意味を持っています。マッハ数が大きいほど、同じ流速の変化に対する密度の変化が大きくなるのです。つまり、超音速の流れでは圧縮性の効果が支配的になります。
次に、この結果を連続の式の微分形 (2) に代入します。
$$ -M^2 \frac{dv}{v} + \frac{dA}{A} + \frac{dv}{v} = 0 $$
$dv/v$ を括り出して整理すると、最終的に以下の面積-速度関係式が得られます。
$$ \begin{equation} \frac{dA}{A} = (M^2 – 1) \frac{dv}{v} \end{equation} $$
これがノズル内の流れの振る舞いを決定する最も重要な式です。
面積-速度関係式の物理的解釈
式 (7) は一見シンプルですが、非常に深い物理を含んでいます。$M^2 – 1$ の符号によって、流れの振る舞いが劇的に変わるからです。場合分けして考えましょう。
場合1: 亜音速（$M < 1$）のとき
$M^2 – 1 < 0$ なので、$dA$ と $dv$ は逆符号になります。
- 断面積が減少する（$dA < 0$）ならば、流速は増加する（$dv > 0$）
- 断面積が増加する（$dA > 0$）ならば、流速は減少する（$dv < 0$）
これは日常の直感に合致します。水道のホースの先をつまめば水が速く出る、というのがまさにこれです。亜音速では圧縮性効果が小さく、密度がほぼ一定のため、断面積の減少は流速の増加に直結します。
場合2: 超音速（$M > 1$）のとき
$M^2 – 1 > 0$ なので、$dA$ と $dv$ は同符号になります。
- 断面積が増加する（$dA > 0$）ならば、流速も増加する（$dv > 0$）
- 断面積が減少する（$dA < 0$）ならば、流速も減少する（$dv < 0$）
これは亜音速の場合と正反対です。超音速で加速するには、流路を広げなければならないのです。なぜでしょうか。超音速の流れでは圧縮性効果が非常に強く、流速が増すと密度が急激に低下します。連続の式 $\rho A v = \text{const}$ において、$v$ の増加と $\rho$ の急激な低下を同時に満たすためには、$A$ を増加させて補う必要があるのです。式 (6) で見たように、$M > 1$ では $|d\rho/\rho| > |dv/v|$ であり、密度の変化率が流速の変化率を上回っています。
場合3: 音速（$M = 1$）のとき
$M^2 – 1 = 0$ なので、$dv$ の値に関係なく $dA = 0$ でなければなりません。つまり音速に達する位置では断面積が極値（最大または最小）を取ります。物理的に意味のある解は断面積が最小の位置、すなわちスロートです。
この結論は非常に重要です。亜音速から超音速へ連続的に加速するためには、まず断面積を狭めて亜音速の流れを加速し（先細部）、スロートで音速に到達させ（$M = 1$, $dA = 0$）、その後断面積を広げて超音速の流れをさらに加速する（末広部）という手順が必要になるのです。これがまさにラバルノズルの形状です。
以上で、面積-速度関係式の導出とその物理的解釈が完了しました。次に、スロートでの臨界条件をもう少し詳しく見ていきましょう。等エントロピー関係式を使って、スロートでの状態量と面積比の関係を定量的に求めます。
スロートでの臨界条件
よどみ点状態量と等エントロピー関係
ノズル流れの解析では、よどみ点（stagnation point、滞留点）の状態量が重要な役割を果たします。よどみ点とは、流れを等エントロピー的に速度ゼロまで減速させたときの仮想的な状態です。圧縮性流体の理論から、よどみ点の温度 $T_0$、圧力 $p_0$、密度 $\rho_0$ と、流れの局所的な温度 $T$、圧力 $p$、密度 $\rho$ の間には、以下の等エントロピー関係式が成り立ちます。
$$ \begin{equation} \frac{T_0}{T} = 1 + \frac{\gamma – 1}{2} M^2 \end{equation} $$
$$ \begin{equation} \frac{p_0}{p} = \left(1 + \frac{\gamma – 1}{2} M^2\right)^{\frac{\gamma}{\gamma – 1}} \end{equation} $$
$$ \begin{equation} \frac{\rho_0}{\rho} = \left(1 + \frac{\gamma – 1}{2} M^2\right)^{\frac{1}{\gamma – 1}} \end{equation} $$
これらの式は、理想気体の状態方程式 $p = \rho R T$ と等エントロピー関係 $p / \rho^\gamma = \text{const}$ から導かれます。直感的には、流れのマッハ数が大きいほど運動エネルギーが大きく、その分だけ温度・圧力・密度が低い状態になっていることを表しています。
臨界状態量（スロートでの値）
スロートではマッハ数がちょうど1に達するので、$M = 1$ を上の等エントロピー関係式に代入すれば、スロートでの状態量が得られます。スロートでの値には上付きの $*$（アスタリスク）を付けて表記する慣習があります。
$M = 1$ を式 (8) に代入すると、臨界温度比が求まります。
$$ \begin{equation} \frac{T^*}{T_0} = \frac{2}{\gamma + 1} \end{equation} $$
同様に、$M = 1$ を式 (9) に代入すると、臨界圧力比が得られます。
$$ \begin{equation} \frac{p^*}{p_0} = \left(\frac{2}{\gamma + 1}\right)^{\frac{\gamma}{\gamma – 1}} \end{equation} $$
$M = 1$ を式 (10) に代入すると、臨界密度比が得られます。
$$ \begin{equation} \frac{\rho^*}{\rho_0} = \left(\frac{2}{\gamma + 1}\right)^{\frac{1}{\gamma – 1}} \end{equation} $$
空気（$\gamma = 1.4$）の場合の具体的な値を計算してみましょう。
$$ \frac{T^*}{T_0} = \frac{2}{2.4} = 0.8333 $$
$$ \frac{p^*}{p_0} = \left(\frac{2}{2.4}\right)^{3.5} = 0.5283 $$
$$ \frac{\rho^*}{\rho_0} = \left(\frac{2}{2.4}\right)^{2.5} = 0.6340 $$
スロートでの圧力はよどみ点圧力の約53%、温度は約83%にまで低下しています。これらの値は理想気体の比熱比 $\gamma$ だけで決まり、ノズルの具体的な形状には依存しません。
面積比-マッハ数関係
ノズル設計において最も実用的な関係式の一つが、面積比-マッハ数関係です。任意の断面でのマッハ数 $M$ と断面積 $A$、スロートの断面積 $A^*$ の比を求めましょう。
連続の式 $\rho A v = \rho^* A^* v^*$ から、面積比は次のように表されます。
$$ \frac{A}{A^*} = \frac{\rho^* v^*}{\rho v} $$
右辺の密度比と速度比を等エントロピー関係式で書き換えます。$v = Ma$、$v^* = a^*$、$a = \sqrt{\gamma R T}$ を用いて整理していきましょう。
まず、速度比については次のようになります。
$$ \frac{v^*}{v} = \frac{a^*}{Ma} = \frac{1}{M}\sqrt{\frac{T^*}{T}} = \frac{1}{M}\sqrt{\frac{T^*/T_0}{T/T_0}} $$
式 (8) と式 (11) を代入すると、次のようになります。
$$ \frac{v^*}{v} = \frac{1}{M}\sqrt{\frac{\frac{2}{\gamma+1}}{1 / \left(1 + \frac{\gamma-1}{2}M^2\right)}} = \frac{1}{M}\sqrt{\frac{2}{\gamma+1}\left(1 + \frac{\gamma-1}{2}M^2\right)} $$
密度比についても同様に式 (10) と式 (13) から求まります。
$$ \frac{\rho^*}{\rho} = \frac{\rho^*/\rho_0}{\rho/\rho_0} = \frac{\left(\frac{2}{\gamma+1}\right)^{\frac{1}{\gamma-1}}}{\left(1 + \frac{\gamma-1}{2}M^2\right)^{-\frac{1}{\gamma-1}}} = \left[\frac{2}{\gamma+1}\left(1 + \frac{\gamma-1}{2}M^2\right)\right]^{\frac{1}{\gamma-1}} $$
これらを掛け合わせて $A/A^*$ を求めると、最終的に以下の面積比-マッハ数関係式が得られます。
$$ \begin{equation} \frac{A}{A^*} = \frac{1}{M}\left[\frac{2}{\gamma+1}\left(1 + \frac{\gamma-1}{2}M^2\right)\right]^{\frac{\gamma+1}{2(\gamma-1)}} \end{equation} $$
この式は超音速ノズルの設計における最も基本的な式であり、所望のマッハ数を出口で達成するために必要な面積比（出口面積/スロート面積）を求めるのに使われます。
式 (14) の重要な性質として、同じ面積比 $A/A^*$ に対して2つのマッハ数の解（1つは亜音速、もう1つは超音速）が存在することが挙げられます。$M = 1$ のとき $A/A^* = 1$（スロート）で最小値を取り、$M$ が1から離れるにつれて $A/A^*$ は単調に増加します。ノズルの先細部では亜音速の解が、末広部では超音速の解が実現されます。
ここまでで、スロートでの臨界条件と面積比-マッハ数関係式が得られました。Pythonで面積比-マッハ数関係を可視化し、理論の理解を深めましょう。
Pythonで面積比-マッハ数関係を可視化する
面積比-マッハ数関係式 (14) の振る舞いを確認するために、$M$ に対する $A/A^*$ をプロットします。
import numpy as np
import matplotlib.pyplot as plt
# 比熱比（空気）
gamma = 1.4
# マッハ数の範囲
M = np.linspace(0.01, 4.0, 1000)
# 面積比-マッハ数関係式
exponent = (gamma + 1) / (2 * (gamma - 1))
A_ratio = (1 / M) * ((2 / (gamma + 1)) * (1 + (gamma - 1) / 2 * M**2))**exponent
# 可視化
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(M, A_ratio, 'b-', linewidth=2)
ax.axvline(x=1.0, color='r', linestyle='--', alpha=0.7, label='M = 1 (throat)')
ax.axhline(y=1.0, color='gray', linestyle=':', alpha=0.5)
ax.set_xlabel('Mach number $M$', fontsize=13)
ax.set_ylabel('Area ratio $A / A^*$', fontsize=13)
ax.set_title('Area ratio vs Mach number (isentropic flow, $\\gamma = 1.4$)', fontsize=14)
ax.set_xlim(0, 4)
ax.set_ylim(0, 12)
ax.legend(fontsize=12)
ax.grid(True, alpha=0.3)
# 亜音速と超音速の領域を注記
ax.annotate('Subsonic\n(convergent section)', xy=(0.4, 2.5), fontsize=11, ha='center',
bbox=dict(boxstyle='round,pad=0.3', facecolor='lightyellow', alpha=0.8))
ax.annotate('Supersonic\n(divergent section)', xy=(2.8, 4.5), fontsize=11, ha='center',
bbox=dict(boxstyle='round,pad=0.3', facecolor='lightyellow', alpha=0.8))
plt.tight_layout()
plt.show()
上のグラフから、面積比-マッハ数関係の重要な特徴が読み取れます。
- $M = 1$ で面積比が最小値 $A/A^* = 1$ を取ります。これがスロートに対応し、音速はスロートでのみ達成されることが視覚的にも確認できます。
- 同じ面積比に対して亜音速と超音速の2つの解が存在することがわかります。例えば $A/A^* = 2$ に対して $M \approx 0.31$ と $M \approx 2.20$ の2つのマッハ数が対応します。先細部では亜音速の解が、末広部では超音速の解が実現されます。
- 超音速側では面積比が急速に増大します。マッハ4のノズルでは $A/A^* \approx 10.7$ もの面積比が必要です。これがロケットエンジンのベルノズルが大きく開いた形状をしている理由です。
次に、この面積比関係を使って、ノズル全体にわたるマッハ数・圧力・温度の分布を計算してみましょう。
チョーク流れと最大質量流量
チョーク現象とは
チョーク（choking）とは、ノズルのスロートでマッハ数が1に達し、それ以上背圧を下げても質量流量が増加しなくなる現象です。日本語では「閉塞」とも呼ばれます。
イメージとしては、高速道路のボトルネック（車線が減少する区間）を想像してください。車の流れがボトルネックで渋滞すると、上流側の車がいくら増えてもボトルネックを通過する車の台数は一定の上限に達します。これと同様に、ノズルのスロートを通過する質量流量にも上限が存在するのです。
先細ノズルを使ったタンクからの噴出を考えましょう。タンク内の圧力（よどみ点圧力 $p_0$）が一定で、出口の背圧 $p_b$ を徐々に下げていく場合を考えます。
- $p_b$ がわずかに $p_0$ より低い場合: ノズル全体で亜音速の流れが生じ、出口でのマッハ数は1未満です。質量流量は $p_b$ を下げるほど増加します。
- $p_b$ が臨界圧力 $p^*$ に達した場合: 出口でちょうど $M = 1$ になります。これがチョーク状態です。
- $p_b$ をさらに下げた場合: ノズル内の流れは変化せず、質量流量は一定のままです。出口の下流で膨張波が生じて圧力を $p_b$ まで下げますが、ノズル内部の流れはこの情報を受け取れません。なぜなら、スロートで流速が音速に達しているため、下流からの圧力波がスロートより上流に伝わることができないからです。
最後の点は非常に重要です。音速とは圧力の微小擾乱が伝搬する速度のことです。流速が音速に達すると、下流からの圧力波は流れに逆らって上流に伝わることができなくなります。そのため、下流の条件（背圧）の変化がノズル内部の流れに影響を与えることができないのです。
最大質量流量の計算
チョーク状態での質量流量（最大質量流量）を定量的に求めましょう。スロートでの質量流量は次の式で表されます。
$$ \dot{m} = \rho^* A^* v^* $$
スロートでは $M = 1$ なので $v^* = a^* = \sqrt{\gamma R T^*}$ です。理想気体の状態方程式 $\rho^* = p^* / (RT^*)$ を使い、式 (11) と (12) の臨界状態量を代入します。
まず $v^* = \sqrt{\gamma R T^*}$ に式 (11) を用いて $T^*$ をよどみ点温度で表します。
$$ v^* = \sqrt{\gamma R \cdot \frac{2 T_0}{\gamma + 1}} $$
密度は $\rho^* = p^*/(RT^*)$ であり、式 (11) と (12) を用いると次のようになります。
$$ \rho^* = \frac{p_0}{RT_0} \cdot \left(\frac{2}{\gamma+1}\right)^{\frac{\gamma}{\gamma-1}} \cdot \frac{T_0}{T^*} = \frac{p_0}{RT_0} \cdot \left(\frac{2}{\gamma+1}\right)^{\frac{\gamma}{\gamma-1}} \cdot \frac{\gamma+1}{2} $$
整理すると、$\rho^*$ は次のようになります。
$$ \rho^* = \frac{p_0}{RT_0} \left(\frac{2}{\gamma+1}\right)^{\frac{1}{\gamma-1}} $$
これらを掛け合わせて質量流量を計算すると、最終的に以下のチョーク質量流量の式が得られます。
$$ \begin{equation} \dot{m}_{\max} = A^* p_0 \sqrt{\frac{\gamma}{R T_0}} \left(\frac{2}{\gamma + 1}\right)^{\frac{\gamma + 1}{2(\gamma – 1)}} \end{equation} $$
この式から、チョーク状態での質量流量を決定する因子が明確にわかります。
- スロート面積 $A^*$: スロートが大きいほど質量流量は大きくなります。これは直感に合います。
- よどみ点圧力 $p_0$: 圧力に比例します。タンク圧力を2倍にすれば質量流量も2倍です。
- よどみ点温度 $T_0$: $\sqrt{T_0}$ に反比例します。温度が高いほど密度が下がるため、質量流量は減少します。
- ガスの種類（$\gamma$ と $R$）: 分子量が小さいガス（$R$ が大きい）ほど音速が大きいですが、密度が小さいため質量流量は減少します。
ロケットエンジンの推力は質量流量と排気速度の積に比例するので（ロケット方程式参照）、スロート面積とチャンバー圧力の設計がエンジン性能を直接左右することがわかります。
チョーク条件と最大質量流量が理解できたところで、次にノズル全体の圧力分布がどのように変化するかを見ていきましょう。特に、背圧の値によって流れのパターンが大きく変わる点が重要です。
ノズル内の圧力分布パターン
背圧による流れの分類
ラバルノズルの動作を完全に理解するには、下流の背圧（back pressure）$p_b$ とよどみ点圧力 $p_0$ の比がどのように流れを支配するかを知る必要があります。背圧を徐々に下げていくと、ノズル内の流れパターンは以下のように変化します。
パターン1: 全体が亜音速（$p_b$ が十分高い場合）
背圧がよどみ点圧力に近い場合、ノズル全体で亜音速の流れが維持されます。先細部で加速され、スロートで最大速度に達しますが、$M < 1$ のままです。末広部では流れが減速され、圧力は回復します。この状態ではノズルは事実上ベンチュリ管と同じ動作をしています。
パターン2: スロートで $M = 1$、末広部で衝撃波が発生
背圧をさらに下げると、スロートで $M = 1$ に達します。末広部では超音速の流れが一時的に発生しますが、背圧が十分に低くないため、末広部のどこかで垂直衝撃波が発生し、超音速の流れが一気に亜音速に戻ります。衝撃波の後方では圧力が急上昇します。背圧の値によって衝撃波の位置が変化し、背圧が低いほど衝撃波はノズル出口に近づきます。
パターン3: 過膨張（over-expanded）
衝撃波がノズル出口の外に出てしまう場合です。ノズル内部では完全に等エントロピーな超音速流れが実現されていますが、出口圧力 $p_e$ が背圧 $p_b$ よりも低い状態です。ノズル出口の直後で斜め衝撃波が発生し、ジェットが「圧縮」されます。過膨張ノズルでは、衝撃波のダイヤモンド模様（マッハディスクやバレルショック）が噴流中に観察されることがあります。
パターン4: 完全膨張（design condition）
出口圧力 $p_e$ がちょうど背圧 $p_b$ に等しい場合です。これが設計条件（design condition）であり、ノズルが最も効率よく動作する状態です。ノズル内部では入口から出口まで完全に等エントロピーな流れが実現され、出口で衝撃波も膨張波も発生しません。ノズルの面積比は、この設計条件でのマッハ数に合わせて決定されます。
パターン5: 不足膨張（under-expanded）
出口圧力 $p_e$ が背圧 $p_b$ よりも高い場合です。ノズル内の流れは完全膨張と同じですが、出口から噴出したガスがさらに膨張波によって広がります。ロケットが高高度で噴射する場合にこの状態になることが多く、ジェットが大きく広がるのが特徴です。
各パターンの判定条件
これらのパターンを定量的に判定するには、ノズル出口面積比 $A_e/A^*$ から等エントロピー関係を使って出口のマッハ数 $M_e$ と圧力 $p_e$ を計算し、背圧 $p_b$ と比較します。
等エントロピー流れでノズル全体が超音速の場合、出口のマッハ数 $M_e$ は面積比-マッハ数関係式 (14) の超音速解で与えられ、出口圧力は次の式で求まります。
$$ \begin{equation} p_e = \frac{p_0}{\left(1 + \frac{\gamma-1}{2}M_e^2\right)^{\frac{\gamma}{\gamma-1}}} \end{equation} $$
- $p_b = p_e$: 完全膨張（設計条件）
- $p_b > p_e$: 過膨張（出口の外で衝撃波）
- $p_b < p_e$: 不足膨張（出口の外で膨張波）
実際のロケットエンジンでは、打ち上げ時（地上での高い大気圧）から高高度（真空に近い大気圧）まで背圧が大きく変化するため、1つのノズルで常に完全膨張を維持することは不可能です。このため、設計では最も重要な高度（例えば推力が最大に必要な瞬間）で完全膨張になるように面積比を選び、それ以外の高度では過膨張や不足膨張を許容します。
Pythonでノズル内のマッハ数・圧力・温度分布を計算する
ラバルノズルの形状を定義し、ノズル内のマッハ数・圧力・温度の分布を数値計算で求めてみましょう。面積比からマッハ数を求めるには、式 (14) の逆関数が必要ですが、解析的には求まらないためニュートン法を使います。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq
# 比熱比（空気）
gamma = 1.4
def area_mach_relation(M, gamma=1.4):
"""面積比-マッハ数関係式: A/A* を計算"""
exp = (gamma + 1) / (2 * (gamma - 1))
return (1 / M) * ((2 / (gamma + 1)) * (1 + (gamma - 1) / 2 * M**2))**exp
def mach_from_area_ratio(A_ratio, gamma=1.4, supersonic=False):
"""面積比からマッハ数を求める（Brent法）"""
def func(M):
return area_mach_relation(M, gamma) - A_ratio
if supersonic:
return brentq(func, 1.0 + 1e-10, 50.0)
else:
return brentq(func, 1e-6, 1.0 - 1e-10)
# ノズル形状の定義（放物線プロファイル）
N = 500
x = np.linspace(-1, 1, N) # 正規化座標（-1: 入口, 0: スロート, 1: 出口）
# スロートで最小、入口と出口で面積が大きい放物線形状
A_ratio_profile = 1.0 + 1.5 * x**2 # A/A*
# マッハ数分布の計算
M_sub = np.zeros(N) # 全体亜音速
M_sup = np.zeros(N) # スロートから超音速
for i in range(N):
ar = A_ratio_profile[i]
# 亜音速解
M_sub[i] = mach_from_area_ratio(ar, gamma, supersonic=False)
# 超音速解（スロートより後方）
if x[i] >= 0:
M_sup[i] = mach_from_area_ratio(ar, gamma, supersonic=True)
else:
M_sup[i] = mach_from_area_ratio(ar, gamma, supersonic=False)
# 等エントロピー関係式による圧力比・温度比
def p_ratio(M, gamma=1.4):
return (1 + (gamma - 1) / 2 * M**2)**(- gamma / (gamma - 1))
def T_ratio(M, gamma=1.4):
return (1 + (gamma - 1) / 2 * M**2)**(-1)
# 可視化（3段パネル + ノズル形状）
fig, axes = plt.subplots(4, 1, figsize=(10, 12), sharex=True)
# (a) ノズル形状
axes[0].fill_between(x, -np.sqrt(A_ratio_profile), np.sqrt(A_ratio_profile),
color='lightblue', alpha=0.5, label='Nozzle wall')
axes[0].plot(x, np.sqrt(A_ratio_profile), 'b-', linewidth=1.5)
axes[0].plot(x, -np.sqrt(A_ratio_profile), 'b-', linewidth=1.5)
axes[0].set_ylabel('$\\sqrt{A/A^*}$', fontsize=12)
axes[0].set_title('Laval nozzle: isentropic flow properties', fontsize=14)
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)
# (b) マッハ数分布
axes[1].plot(x, M_sub, 'b-', linewidth=2, label='Fully subsonic')
axes[1].plot(x, M_sup, 'r-', linewidth=2, label='Supersonic (design)')
axes[1].axhline(y=1.0, color='gray', linestyle=':', alpha=0.5)
axes[1].set_ylabel('Mach number $M$', fontsize=12)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)
# (c) 圧力分布
axes[2].plot(x, p_ratio(M_sub), 'b-', linewidth=2, label='Fully subsonic')
axes[2].plot(x, p_ratio(M_sup), 'r-', linewidth=2, label='Supersonic (design)')
axes[2].set_ylabel('$p / p_0$', fontsize=12)
axes[2].legend(fontsize=10)
axes[2].grid(True, alpha=0.3)
# (d) 温度分布
axes[3].plot(x, T_ratio(M_sub), 'b-', linewidth=2, label='Fully subsonic')
axes[3].plot(x, T_ratio(M_sup), 'r-', linewidth=2, label='Supersonic (design)')
axes[3].set_ylabel('$T / T_0$', fontsize=12)
axes[3].set_xlabel('Normalized position $x/L$', fontsize=12)
axes[3].legend(fontsize=10)
axes[3].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
上の4段パネルのグラフから、以下のことが読み取れます。
-
ノズル形状（最上段）: 入口とスロートの間で断面積が減少し、スロート（$x = 0$）で最小になり、出口に向かって再び増加する、典型的なラバルノズルの形状です。
-
マッハ数分布（2段目）: 青線（全体亜音速）では、スロートで最大マッハ数に達した後、末広部で減速して出口マッハ数は入口と対称的な値に戻ります。赤線（超音速設計条件）では、スロートで $M = 1$ に達した後、末広部でマッハ数が単調に増加し続けます。出口マッハ数は約2.0に達しています。
-
圧力分布（3段目）: 超音速設計条件の場合、圧力はノズル全体を通じて単調に減少し、出口では $p/p_0 \approx 0.09$（よどみ点圧力の約9%）まで低下しています。一方、全体亜音速の場合は末広部で圧力が回復します。
-
温度分布（4段目）: 超音速設計条件では、出口温度は $T/T_0 \approx 0.56$ まで低下しています。この温度低下はガスの運動エネルギーへの変換に伴うもので、熱力学第1法則（エンタルピー保存）と整合しています。
これらの分布はノズルの形状（面積比プロファイル）が決まれば一意に定まることが、等エントロピー流れの重要な特徴です。次に、背圧を変化させたときにノズル内の流れパターンがどのように変わるかをPythonで可視化します。
Pythonで背圧変化による流れパターンを可視化する
背圧 $p_b / p_0$ を変化させたときのノズル内圧力分布を計算し、先に述べた5つのパターンを一枚のグラフで確認しましょう。末広部で衝撃波が発生する場合は、垂直衝撃波の関係式を使って衝撃波前後の状態量の跳びを計算します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq
gamma = 1.4
def area_mach(M, gamma=1.4):
"""面積比 A/A*"""
exp = (gamma + 1) / (2 * (gamma - 1))
return (1 / M) * ((2 / (gamma + 1)) * (1 + (gamma - 1) / 2 * M**2))**exp
def mach_from_area(A_ratio, gamma=1.4, supersonic=False):
"""面積比からマッハ数を求める"""
def func(M):
return area_mach(M, gamma) - A_ratio
if supersonic:
return brentq(func, 1.0 + 1e-10, 50.0)
else:
return brentq(func, 1e-6, 1.0 - 1e-10)
def p_over_p0(M, gamma=1.4):
"""等エントロピー圧力比 p/p0"""
return (1 + (gamma - 1) / 2 * M**2)**(- gamma / (gamma - 1))
def normal_shock_p02_p01(M1, gamma=1.4):
"""垂直衝撃波後のよどみ点圧力比 p02/p01"""
term1 = ((gamma + 1) * M1**2 / ((gamma - 1) * M1**2 + 2))**(gamma / (gamma - 1))
term2 = ((gamma + 1) / (2 * gamma * M1**2 - (gamma - 1)))**(1 / (gamma - 1))
return term1 * term2
def normal_shock_M2(M1, gamma=1.4):
"""垂直衝撃波後のマッハ数"""
return np.sqrt(((gamma - 1) * M1**2 + 2) / (2 * gamma * M1**2 - (gamma - 1)))
# ノズル形状（面積比プロファイル）
N = 500
x = np.linspace(-1, 1, N)
A_ratio = 1.0 + 1.5 * x**2 # A/A*
# 出口面積比
Ae_At = A_ratio[-1] # = 2.5
# 各パターンの圧力分布を計算
fig, ax = plt.subplots(figsize=(10, 6))
# パターン1: 全体亜音速（背圧が高い）
M_all_sub = np.array([mach_from_area(ar, supersonic=False) for ar in A_ratio])
p_pattern1 = p_over_p0(M_all_sub)
ax.plot(x, p_pattern1, 'b-', linewidth=2, label='Pattern 1: Fully subsonic')
# パターン4: 完全膨張（設計条件）
M_design = np.zeros(N)
for i in range(N):
if x[i] < 0:
M_design[i] = mach_from_area(A_ratio[i], supersonic=False)
else:
M_design[i] = mach_from_area(A_ratio[i], supersonic=True)
p_design = p_over_p0(M_design)
ax.plot(x, p_design, 'r-', linewidth=2, label='Pattern 4: Fully expanded (design)')
# パターン2: 末広部で衝撃波（衝撃波位置を x_shock = 0.4 に設定）
x_shock = 0.4
idx_shock = np.argmin(np.abs(x - x_shock))
M_shock_pattern = np.zeros(N)
p_shock_pattern = np.zeros(N)
for i in range(N):
if x[i] < 0:
M_shock_pattern[i] = mach_from_area(A_ratio[i], supersonic=False)
elif i <= idx_shock:
M_shock_pattern[i] = mach_from_area(A_ratio[i], supersonic=True)
else:
# 衝撃波後: p02/p01 を使って亜音速解を計算
M1_shock = mach_from_area(A_ratio[idx_shock], supersonic=True)
p02_p01 = normal_shock_p02_p01(M1_shock)
M2_shock = normal_shock_M2(M1_shock)
# 衝撃波直後の面積比 (A/A*_new)
A_star_new_ratio = A_ratio[idx_shock] / area_mach(M2_shock)
A_ratio_new = A_ratio[i] / A_star_new_ratio
M_shock_pattern[i] = mach_from_area(A_ratio_new, supersonic=False)
# 圧力の計算（衝撃波前後でよどみ点圧力が変わる）
for i in range(N):
if i <= idx_shock:
p_shock_pattern[i] = p_over_p0(M_shock_pattern[i])
else:
M1_shock = mach_from_area(A_ratio[idx_shock], supersonic=True)
p02_p01 = normal_shock_p02_p01(M1_shock)
p_shock_pattern[i] = p02_p01 * p_over_p0(M_shock_pattern[i])
ax.plot(x, p_shock_pattern, 'g-', linewidth=2,
label=f'Pattern 2: Normal shock at x={x_shock:.1f}')
ax.plot(x[idx_shock], p_shock_pattern[idx_shock], 'g^', markersize=10)
# パターン5: 不足膨張（背圧が設計値より低い）
p_under = np.copy(p_design)
pb_under = p_design[-1] * 0.5 # 設計出口圧力の半分
ax.plot(x, p_under, 'r--', linewidth=1.5, alpha=0.5)
ax.annotate('Pattern 5:\nUnder-expanded\n$p_b < p_e$', xy=(1.0, pb_under),
xytext=(0.6, 0.02), fontsize=10, arrowprops=dict(arrowstyle='->', color='purple'),
color='purple', bbox=dict(boxstyle='round', facecolor='lavender', alpha=0.8))
# パターン3: 過膨張
pb_over = p_design[-1] * 2.0
ax.annotate('Pattern 3:\nOver-expanded\n$p_b > p_e$', xy=(1.0, pb_over),
xytext=(0.6, 0.35), fontsize=10, arrowprops=dict(arrowstyle='->', color='orange'),
color='orange', bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))
ax.axvline(x=0, color='gray', linestyle=':', alpha=0.5, label='Throat')
ax.set_xlabel('Normalized position $x/L$', fontsize=13)
ax.set_ylabel('$p / p_0$', fontsize=13)
ax.set_title('Pressure distributions in a Laval nozzle for various back pressures', fontsize=14)
ax.legend(fontsize=10, loc='upper right')
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 1.05)
plt.tight_layout()
plt.show()
上のグラフから、背圧の変化に対するノズル内の流れパターンの違いが明確に見えます。
-
全体亜音速（青線）: 圧力はスロートで最小値を取った後、末広部で回復します。このパターンではノズルは効率的な加速装置として機能していません。末広部がディフューザ（減速・圧力回復装置）として働いています。
-
末広部で衝撃波（緑線）: スロートから衝撃波の位置（緑の三角形）まで超音速の等エントロピー膨張が起こり、圧力は急激に低下します。衝撃波の位置で圧力が不連続的に跳び上がり、それ以降は亜音速の流れとなって緩やかに圧力が回復します。衝撃波はエントロピーを増大させるため、この跳びの後ではよどみ点圧力が低下しています。
-
完全膨張（赤線）: 圧力はノズル全体を通じて単調に減少し、出口で設計背圧にちょうど一致します。これが最も効率的な動作条件です。
-
過膨張と不足膨張: ノズル内部の流れは完全膨張と同じですが、出口の下流で衝撃波（過膨張）または膨張波（不足膨張）が発生して圧力を背圧に調整します。
ここまでで、ラバルノズルの流れの基本メカニズムを理解しました。次に、ロケットエンジンの文脈でラバルノズルがどのように利用されるかを見ていきましょう。推力の計算と推力係数の概念が実用上重要です。
ロケットノズルとの関連
ロケットエンジンにおけるラバルノズルの役割
ロケットエンジンは、燃焼室で高温・高圧のガスを生成し、それをラバルノズルで超音速まで加速して噴射することで推力を得ます。ロケット方程式から、ロケットの推力 $F$ は次のように表されます。
$$ \begin{equation} F = \dot{m} v_e + (p_e – p_a) A_e \end{equation} $$
ここで $\dot{m}$ は質量流量、$v_e$ は出口での排気速度、$p_e$ はノズル出口圧力、$p_a$ は周囲大気圧、$A_e$ は出口断面積です。
第1項 $\dot{m} v_e$ は運動量推力（momentum thrust）と呼ばれ、排気ガスの運動量変化に起因する推力です。第2項 $(p_e – p_a) A_e$ は圧力推力（pressure thrust）と呼ばれ、出口圧力と大気圧の差によって生じる推力です。
完全膨張条件（$p_e = p_a$）では圧力推力はゼロになり、推力は運動量推力のみとなります。過膨張（$p_e < p_a$）では圧力推力が負となり、不足膨張（$p_e > p_a$）では圧力推力が正となります。ただし、不足膨張では排気速度 $v_e$ が完全膨張のときより小さくなるため、総推力が最大になるのは必ずしも不足膨張のときではありません。
推力係数
ロケットノズルの性能を表す無次元パラメータとして、推力係数 $C_F$ が広く使われています。推力係数は、推力を「チャンバー圧力 $\times$ スロート面積」で無次元化したものです。
$$ \begin{equation} C_F = \frac{F}{p_0 A^*} \end{equation} $$
式 (17) と等エントロピー関係を用いて $C_F$ を展開すると、次の表現が得られます。
$$ \begin{equation} C_F = \sqrt{\frac{2\gamma^2}{\gamma-1}\left(\frac{2}{\gamma+1}\right)^{\frac{\gamma+1}{\gamma-1}} \left[1 – \left(\frac{p_e}{p_0}\right)^{\frac{\gamma-1}{\gamma}}\right]} + \frac{p_e – p_a}{p_0}\frac{A_e}{A^*} \end{equation} $$
第1項は運動量推力からの寄与で、出口圧力比 $p_e/p_0$ が小さい（高膨張比）ほど大きくなります。第2項は圧力推力からの寄与です。
推力係数はノズルの幾何形状（面積比）とガスの性質（$\gamma$）だけで決まり、チャンバー圧力やスロート面積に依存しません。このため、異なるエンジンのノズル性能を公平に比較するのに適した指標です。
有効排気速度と比推力
ロケットの性能指標として比推力 $I_{sp}$ がよく使われます。比推力は「単位重量流量あたりの推力」であり、次のように定義されます。
$$ \begin{equation} I_{sp} = \frac{F}{\dot{m} g_0} = \frac{v_{\text{eff}}}{g_0} \end{equation} $$
ここで $v_{\text{eff}} = F / \dot{m}$ は有効排気速度（effective exhaust velocity）、$g_0 = 9.80665 \, \text{m/s}^2$ は標準重力加速度です。有効排気速度は運動量推力と圧力推力の両方を含んでおり、完全膨張の場合は実際の排気速度 $v_e$ に等しくなります。
チョーク質量流量の式 (15) を推力係数の定義 (18) に代入すると、有効排気速度は次のように表されます。
$$ v_{\text{eff}} = C_F \cdot \frac{p_0 A^*}{\dot{m}} = C_F \cdot \frac{1}{\sqrt{\frac{\gamma}{RT_0}}\left(\frac{2}{\gamma+1}\right)^{\frac{\gamma+1}{2(\gamma-1)}}} $$
この式から、比推力を最大化するには高い燃焼温度 $T_0$（大きな $v_e$）と適切なノズル膨張比（大きな $C_F$）が重要であることがわかります。水素-酸素エンジンが高い比推力を達成できるのは、水素の分子量が小さく（$R$ が大きく）、燃焼温度が高いためです。
実際のロケットノズル設計
実際のロケットノズルの設計では、以下のような追加的な考慮が必要になります。
-
ノズル輪郭の最適化: 理想的なラバルノズルは面積比の変化が滑らかであれば機能しますが、実際にはノズルの長さを短くしつつ効率を維持する必要があります。ベルノズル（釣鐘型）やラオの方法による最適輪郭が使われます。
-
境界層の影響: 実際のノズルでは壁面に境界層が発達し、有効断面積が幾何学的断面積より小さくなります。これにより流量係数が1未満になります。
-
化学反応の凍結と平衡: ロケット燃焼ガスは膨張中に化学平衡が「凍結」（frozen）することがあり、理想的な等エントロピー膨張からのずれが生じます。
-
二相流: 固体ロケットの排気には固体粒子（アルミナなど）が含まれ、二相流の効果を考慮する必要があります。
-
ノズル冷却: 高温の燃焼ガスによるノズル壁面の加熱に対処するため、再生冷却（推進剤をノズル壁面に流して冷却する方式）やアブレーション冷却が用いられます。
ここまでで、ラバルノズルとロケット推進の関連を理解しました。最後に、推力係数をPythonで計算し、ノズル膨張比と周囲圧力がロケット性能に与える影響を定量的に確認しましょう。
Pythonで推力係数を計算する
推力係数 $C_F$ をノズル膨張比（出口/スロート面積比）の関数として計算し、大気圧の影響も含めて可視化します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq
gamma = 1.4
def area_mach(M, gamma=1.4):
"""面積比 A/A*"""
exp = (gamma + 1) / (2 * (gamma - 1))
return (1 / M) * ((2 / (gamma + 1)) * (1 + (gamma - 1) / 2 * M**2))**exp
def mach_from_area(A_ratio, gamma=1.4):
"""面積比から超音速マッハ数を求める"""
def func(M):
return area_mach(M, gamma) - A_ratio
return brentq(func, 1.0 + 1e-10, 50.0)
def p_over_p0(M, gamma=1.4):
"""等エントロピー圧力比 p/p0"""
return (1 + (gamma - 1) / 2 * M**2)**(- gamma / (gamma - 1))
def thrust_coefficient(pe_p0, Ae_At, pa_p0, gamma=1.4):
"""推力係数 C_F"""
momentum_term = np.sqrt(
(2 * gamma**2 / (gamma - 1))
* (2 / (gamma + 1))**((gamma + 1) / (gamma - 1))
* (1 - pe_p0**((gamma - 1) / gamma))
)
pressure_term = (pe_p0 - pa_p0) * Ae_At
return momentum_term + pressure_term
# 膨張比の範囲
epsilon = np.linspace(1.01, 80, 1000) # Ae/A*
# 各膨張比に対する出口マッハ数と圧力比
Me = np.array([mach_from_area(e) for e in epsilon])
pe_p0 = p_over_p0(Me)
# 異なる大気圧条件での推力係数
pa_p0_values = [0.0, 0.01, 0.05, 0.1]
labels = ['Vacuum ($p_a = 0$)', '$p_a/p_0 = 0.01$',
'$p_a/p_0 = 0.05$', '$p_a/p_0 = 0.1$']
colors = ['red', 'orange', 'green', 'blue']
fig, ax = plt.subplots(figsize=(10, 6))
for pa_p0, label, color in zip(pa_p0_values, labels, colors):
CF = thrust_coefficient(pe_p0, epsilon, pa_p0)
ax.plot(epsilon, CF, color=color, linewidth=2, label=label)
# 最適膨張比（完全膨張条件: pe = pa）の位置にマーカー
if pa_p0 > 0:
idx_opt = np.argmin(np.abs(pe_p0 - pa_p0))
ax.plot(epsilon[idx_opt], CF[idx_opt], 'o', color=color, markersize=8)
ax.set_xlabel('Nozzle expansion ratio $A_e / A^*$', fontsize=13)
ax.set_ylabel('Thrust coefficient $C_F$', fontsize=13)
ax.set_title('Thrust coefficient vs expansion ratio ($\\gamma = 1.4$)', fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_xlim(1, 80)
ax.set_ylim(0.5, 2.2)
plt.tight_layout()
plt.show()
上のグラフから、推力係数に関する重要な知見が読み取れます。
-
真空中（赤線）では膨張比を大きくするほど推力係数が単調に増加します。これは、膨張比が大きいほど排気ガスのエネルギーをより多く運動エネルギーに変換できるためです。理論上は膨張比を無限大にすれば $C_F$ は理論上限に漸近しますが、実際にはノズルの重量とサイズの制約があります。
-
大気圧が存在する場合（他の線）には、推力係数に最適な膨張比が存在します。丸印が完全膨張条件の位置を示しており、これより膨張比が大きいと過膨張となり圧力推力が負になるため推力係数が低下します。地上用エンジンでは比較的小さな膨張比が最適であり、上段エンジン（真空用）では大きな膨張比が有利です。
-
大気圧が低い環境ほど推力係数が高いことが確認できます。これがロケットの上段エンジンに大きなノズル（高膨張比）が使われる理由です。例えば、スペースXのファルコン9のマーリンエンジンは地上版と真空版で膨張比が大きく異なり、真空版はノズルベルが大きく広がった形状をしています。
数値計算のまとめ: 面積比からマッハ数を求める逆関数
ここまでの計算で繰り返し使ってきた「面積比からマッハ数を求める」操作は、ノズル設計の基本ツールです。面積比-マッハ数関係式 (14) は陽にマッハ数について解くことができないため、数値的に逆関数を求める必要があります。
先ほどのコードではSciPyの brentq（ブレント法）を使いましたが、ここではもう少し詳しく、ニュートン法による実装も紹介しておきましょう。ニュートン法では、面積比-マッハ数関数の導関数を解析的に計算して収束を高速化できます。
$$ \frac{d}{dM}\left(\frac{A}{A^*}\right) = \frac{A}{A^*}\left(-\frac{1}{M} + \frac{(\gamma+1)M/2}{1 + \frac{\gamma-1}{2}M^2}\right) $$
以下のPythonコードで、ブレント法とニュートン法の比較を行います。
import numpy as np
import time
gamma = 1.4
def area_mach(M, gamma=1.4):
"""面積比 A/A*"""
exp = (gamma + 1) / (2 * (gamma - 1))
return (1 / M) * ((2 / (gamma + 1)) * (1 + (gamma - 1) / 2 * M**2))**exp
def dAdM(M, gamma=1.4):
"""面積比の導関数 d(A/A*)/dM"""
A = area_mach(M, gamma)
return A * (-1/M + (gamma + 1) * M / (2 * (1 + (gamma - 1)/2 * M**2)))
def mach_newton(A_ratio, gamma=1.4, supersonic=False, tol=1e-12, max_iter=50):
"""ニュートン法で面積比からマッハ数を求める"""
# 初期値の設定
if supersonic:
M = 1.0 + np.sqrt(A_ratio - 1.0) # 超音速側の初期推定
else:
M = 0.5 # 亜音速側の初期推定
for _ in range(max_iter):
f = area_mach(M, gamma) - A_ratio
fp = dAdM(M, gamma)
dM = -f / fp
M = M + dM
# マッハ数の範囲を制限
if supersonic:
M = max(M, 1.0 + 1e-10)
else:
M = min(max(M, 1e-6), 1.0 - 1e-10)
if abs(dM) < tol:
break
return M
# テスト: 面積比 2.0, 5.0, 10.0 に対する亜音速・超音速解
test_ratios = [2.0, 5.0, 10.0]
print("面積比-マッハ数逆関数の計算結果")
print("-" * 60)
print(f"{'A/A*':>8s} {'M (subsonic)':>14s} {'M (supersonic)':>14s} {'A/A* check':>12s}")
print("-" * 60)
for ar in test_ratios:
M_sub = mach_newton(ar, supersonic=False)
M_sup = mach_newton(ar, supersonic=True)
# 検算
ar_check_sub = area_mach(M_sub)
ar_check_sup = area_mach(M_sup)
print(f"{ar:8.1f} {M_sub:14.6f} {M_sup:14.6f} {ar_check_sub:12.6f}")
# 速度比較
N_test = 10000
A_test = np.linspace(1.01, 20.0, N_test)
start = time.time()
for ar in A_test:
mach_newton(ar, supersonic=True)
t_newton = time.time() - start
from scipy.optimize import brentq
start = time.time()
for ar in A_test:
brentq(lambda M: area_mach(M) - ar, 1.0 + 1e-10, 50.0)
t_brentq = time.time() - start
print(f"\n速度比較 (N = {N_test}):")
print(f" Newton法: {t_newton:.4f} 秒")
print(f" Brent法: {t_brentq:.4f} 秒")
print(f" 速度比: Newton法は Brent法の {t_brentq/t_newton:.1f} 倍速い")
上の出力から確認できる重要なポイントは次の通りです。
-
同じ面積比に対して亜音速と超音速の2つの解が存在することが数値的に確認できます。例えば $A/A^* = 2.0$ に対して $M \approx 0.306$ と $M \approx 2.197$ の2解があります。
-
逆関数の検算: 求めたマッハ数を面積比-マッハ数関係式に代入すると、元の面積比が正しく再現されていることが確認できます。これにより数値解の精度が保証されます。
-
ニュートン法はブレント法より高速です。導関数の解析式が利用可能な場合、ニュートン法の二次収束性が活きて計算効率が大幅に向上します。大規模な設計最適化やCFD（数値流体力学）の初期化で面積比-マッハ数変換を大量に実行する場合に有用です。
まとめ
本記事では、ラバルノズルの原理を圧縮性流体力学の基礎から導出し、ノズル内の流れの振る舞いを解説しました。
- 面積-速度関係式 $dA/A = (M^2 – 1)\,dv/v$ は、亜音速では「絞ると加速」、超音速では「広げると加速」という対照的な振る舞いを表し、音速（$M = 1$）がスロート（$dA = 0$）でのみ達成されることを示します
- 等エントロピー関係式とスロートでの臨界条件から、ノズルの面積比が所望の出口マッハ数を一意に決定することがわかります。面積比-マッハ数関係式は超音速ノズル設計の出発点です
- チョーク現象はスロートで $M = 1$ に達すると質量流量が上限値に固定される現象であり、最大質量流量はスロート面積・よどみ点圧力・よどみ点温度で決まります
- 背圧の変化によって、ノズル内の流れは全体亜音速・末広部での衝撃波・過膨張・完全膨張・不足膨張の5つのパターンに分類されます。ロケットエンジンでは飛行高度によって背圧が変化するため、設計条件の選択が重要です
- 推力係数 $C_F$ はノズルの幾何形状と作動ガスの性質だけで決まる無次元性能指標であり、真空中では膨張比を大きくするほど増加しますが、大気圧下では最適な膨張比が存在します
ラバルノズルの理解は、ロケットノズルの設計や比推力の最適化に直結するだけでなく、衝撃波の理論や超音速空気力学を学ぶ際の基礎にもなります。
次のステップとして、以下の記事も参考にしてください。
- 衝撃波と膨張波 — ノズル内外で発生する衝撃波のメカニズムと跳び条件
- ロケットノズルの設計 — ベルノズル輪郭の最適化と実際の設計手法
- 比推力 — ロケットエンジンの性能評価指標
- ロケット方程式 — 推力と軌道変更の関係