ジェットエンジンはなぜ航空機を前に押し出すことができるのでしょうか。消防士がホースを握りしめるとき、なぜあれほど大きな反力を感じるのでしょうか。これらの問いに共通する鍵は、流体が運ぶ運動量です。流体が加速したり方向を変えたりすると、ニュートンの運動法則に従って力が発生します。この力を系統的に計算する枠組みが流体の運動量の法則です。
運動量の法則は流体力学のなかで最も実用的な道具の一つであり、以下のような幅広い工学分野で日常的に使われています。
- 航空宇宙工学: ジェットエンジンやロケットモータの推力を算出する際、噴流が運ぶ運動量の変化から力を直接計算します。ロケット方程式の土台にもなっています
- 配管・プラント設計: 管路の曲がり部、バルブ、ノズルに作用する力を求めて配管サポートを設計するとき、運動量の法則は欠かせません
- 土木・水工学: 水門に作用する力、河川の堰における衝撃力、ダムの放水路設計など、開水路の問題にも運動量の法則が適用されます
- 風力・翼型設計: 翼に作用する揚力と抗力を検査体積で評価する手法は、風洞実験データの解析において基本的な方法です
本記事では、ニュートンの第2法則から出発して流体の運動量方程式を導出し、検査体積アプローチの考え方を丁寧に解説します。その上で、ジェットの推力、管路の曲がり部に作用する力、衝突噴流の力学、翼周りの検査体積解析という4つの代表的な応用をPythonで実装します。
本記事の内容
- ニュートンの第2法則から流体の運動量方程式への拡張
- 検査体積アプローチとレイノルズ輸送定理の導出
- 定常流における運動量保存式の簡略化
- ジェットエンジンの推力計算
- 管路の曲がり部に作用する力の導出
- 衝突噴流が平板に及ぼす力
- 翼周りの検査体積解析（揚力・抗力の評価）
- Python による4種類の数値計算と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 連続の式 — 質量保存の微分形・積分形と検査体積の基本
- ベルヌーイの定理 — エネルギー保存の観点から流れを理解する
- ニュートンの運動法則 — 質点の運動量と力の関係
- 運動量保存則 — 質点系の運動量保存の基礎
- 流体の性質 — 密度・粘性などの基本的な流体物性
ニュートンの第2法則から流体の運動量方程式へ
質点系の運動量の法則
流体の運動量方程式を理解するために、まず質点系のニュートンの運動法則を振り返りましょう。質量 $m$ の質点に力 $\bm{F}$ が作用するとき、運動量 $\bm{p} = m\bm{v}$ の時間変化率は力に等しいという法則です。
$$ \bm{F} = \frac{d\bm{p}}{dt} = \frac{d(m\bm{v})}{dt} $$
質量が一定であれば $\bm{F} = m\bm{a}$ というおなじみの形になりますが、流体では検査体積を出入りする質量が変化しうるため、$d(m\bm{v})/dt$ の形のほうが本質的です。
この法則を多数の流体粒子（流体要素）の集まりに拡張することを考えます。流体塊（物質体積）全体に作用する外力の合計は、その流体塊全体の運動量の時間変化率に等しいはずです。
$$ \bm{F}_{\text{ext}} = \frac{d}{dt}\int_{\text{物質体積}} \rho \bm{v} \, dV $$
ここで $\rho$ は流体の密度、$\bm{v}$ は速度ベクトル、積分は同一の流体粒子を追跡する物質体積（ラグランジュ的体積）にわたって行います。この式が流体に対するニュートンの第2法則そのものです。
物質体積と検査体積の違い
しかし、流体力学では「同一の流体粒子を追い続ける」という物質体積の考え方は不便なことが多いのです。なぜなら、流体粒子は絶えず動き回り、形を変え、分裂し、混合するからです。
そこで、空間に固定された箱（検査体積、Control Volume: CV）を考え、その箱を出入りする流体の運動量を追跡するほうがはるかに実用的です。たとえば、エンジンの周りに仮想的な箱を設定し、箱に入ってくる空気と箱から出ていく排気ガスの運動量の差を見れば、エンジンが発生する推力がわかります。
この「物質体積で記述された法則を検査体積で記述し直す」ための数学的な道具が、レイノルズ輸送定理です。次節でこの定理を導出します。
レイノルズ輸送定理
直感的な理解
レイノルズ輸送定理は、「流体粒子を追跡する観点（ラグランジュ的）」と「空間に固定された箱で観測する観点（オイラー的）」を結びつける橋渡しの定理です。
身近なアナロジーで考えてみましょう。部屋（検査体積）の中にいる人の総数の変化を考えます。部屋の中の人数が変わる理由は2つあります。(1) 部屋の中にいる人が増えたり減ったりする（部屋内部の変化）、(2) ドアから人が出入りする（境界面を通じたフラックス）。レイノルズ輸送定理はまさにこの2つの寄与を数学的に分離するものです。
定理の導出
任意の流体の物理量（運動量、エネルギーなど）を単位質量あたり $\phi$ で表しましょう。物質体積 $V_m(t)$（同一の流体粒子を含む時間変化する領域）にわたる $\rho\phi$ の総量の時間変化率は次のように書けます。
$$ \frac{d}{dt}\int_{V_m(t)} \rho \phi \, dV $$
ある瞬間 $t = t_0$ において物質体積と検査体積が一致するとします。微小時間 $\Delta t$ 後、物質体積は流れに沿って移動しますが、検査体積 $V_{CV}$ は空間に固定されたままです。この差を考慮すると、次の関係が導かれます。
$$ \frac{d}{dt}\int_{V_m} \rho \phi \, dV = \frac{\partial}{\partial t}\int_{V_{CV}} \rho \phi \, dV + \oint_{CS} \rho \phi \, \bm{v} \cdot \bm{n} \, dA $$
右辺の第1項は、検査体積内部の $\rho\phi$ の時間変化率（蓄積項）です。第2項は、検査面（Control Surface: CS）を通じて流出入する $\rho\phi$ のフラックスです。$\bm{n}$ は検査面の外向き単位法線ベクトルであり、$\bm{v} \cdot \bm{n} > 0$ なら流出、$\bm{v} \cdot \bm{n} < 0$ なら流入を意味します。
この式がレイノルズ輸送定理です。左辺はラグランジュ的な時間微分、右辺はオイラー的な記述に変換されています。
導出のステップ
もう少し丁寧に導出を追いましょう。時刻 $t$ で物質体積 $V_m(t)$ と検査体積 $V_{CV}$ が一致しているとします。時刻 $t + \Delta t$ で、物質体積は流れに乗って動き、一部は検査体積から流出し（体積 $\Delta V_{\text{out}}$）、検査体積の外にあった流体が検査体積内に入ってきます（体積 $\Delta V_{\text{in}}$）。
物質体積内の総量の変化は次のように分解できます。
$$ \frac{d}{dt}\int_{V_m} \rho\phi \, dV = \lim_{\Delta t \to 0} \frac{1}{\Delta t}\left[\int_{V_m(t+\Delta t)} \rho\phi \, dV – \int_{V_m(t)} \rho\phi \, dV \right] $$
$V_m(t) = V_{CV}$ であり、$V_m(t+\Delta t) = V_{CV} – \Delta V_{\text{in}} + \Delta V_{\text{out}}$ と書けるので、上式を整理すると次のようになります。
$$ \frac{d}{dt}\int_{V_m} \rho\phi \, dV = \lim_{\Delta t \to 0}\frac{1}{\Delta t}\left[\int_{V_{CV}} \rho\phi \Big|_{t+\Delta t} dV – \int_{V_{CV}} \rho\phi\Big|_t \, dV\right] + \lim_{\Delta t \to 0}\frac{1}{\Delta t}\left[\int_{\Delta V_{\text{out}}} \rho\phi \, dV – \int_{\Delta V_{\text{in}}} \rho\phi \, dV\right] $$
第1項の極限は検査体積内の偏時間微分 $\frac{\partial}{\partial t}\int_{V_{CV}} \rho\phi \, dV$ に他なりません。第2項について、微小時間 $\Delta t$ の間に検査面の微小面積 $dA$ を通って流出する体積は $(\bm{v}\cdot\bm{n})\,dA\,\Delta t$（外向き正）ですから、第2項は面積分 $\oint_{CS} \rho\phi\,\bm{v}\cdot\bm{n}\,dA$ となります。
こうしてレイノルズ輸送定理が得られました。この定理は運動量だけでなく、質量（$\phi = 1$ とすれば連続の式が得られる）やエネルギー（$\phi$ をエネルギー密度とすればベルヌーイの定理につながる）など、あらゆる保存則の変換に使える汎用的な道具です。
次に、この定理を運動量に適用して、流体の運動量方程式を導きましょう。
検査体積に対する運動量方程式
一般形の導出
レイノルズ輸送定理で $\phi = \bm{v}$（速度ベクトル、つまり単位質量あたりの運動量）と置きます。すると、物質体積に対するニュートンの第2法則は次のように検査体積の言葉で書き直せます。
$$ \bm{F}_{\text{ext}} = \frac{\partial}{\partial t}\int_{V_{CV}} \rho \bm{v} \, dV + \oint_{CS} \rho \bm{v} (\bm{v} \cdot \bm{n}) \, dA $$
この式の各項の意味を確認しましょう。
- 左辺 $\bm{F}_{\text{ext}}$: 検査体積に作用する全ての外力の合計です。これには体積力（重力）と表面力（圧力、せん断応力）が含まれます
- 右辺第1項: 検査体積内の運動量の時間変化率です。非定常な流れ（時間とともに速度場が変化する流れ）でのみゼロでない値を持ちます
- 右辺第2項: 検査面を通じて正味流出する運動量のフラックスです。流体が検査体積から出ていくとき運動量を持ち出し、入ってくるとき運動量を持ち込みます
この式は流体力学における運動量の積分形であり、ナビエ・ストークス方程式の微分形と等価な情報を含んでいます。微分形は流れ場の各点での力のバランスを記述しますが、積分形は有限の領域全体での力のバランスを記述します。工学的な力の計算では、内部の流れの詳細を知らなくても、検査面での入口・出口条件さえわかれば力が求まるという点で、積分形のほうが強力です。
外力の分類
左辺の外力 $\bm{F}_{\text{ext}}$ をもう少し具体的に書き下しましょう。検査体積に作用する力は次の3種類に分類できます。
体積力（ボディフォース）: 検査体積内の流体全体に作用する力です。最も一般的なのは重力 $\rho\bm{g}$ です。
$$ \bm{F}_{\text{body}} = \int_{V_{CV}} \rho \bm{g} \, dV $$
圧力: 検査面に垂直に作用する表面力です。外向き法線 $\bm{n}$ に対して内側に押す力なので、符号は負になります。
$$ \bm{F}_{\text{pressure}} = -\oint_{CS} p \, \bm{n} \, dA $$
固体壁からの反力: 検査体積内に固体壁（管壁、翼面など）がある場合、壁が流体に及ぼす力 $\bm{R}$ があります。これが求めたい力の反作用です。
これらをまとめると、運動量方程式は次の形になります。
$$ \bm{R} + \bm{F}_{\text{body}} + \bm{F}_{\text{pressure}} = \frac{\partial}{\partial t}\int_{V_{CV}} \rho \bm{v} \, dV + \oint_{CS} \rho \bm{v}(\bm{v} \cdot \bm{n}) \, dA $$
ここで $\bm{R}$ は固体が流体に及ぼす力です。逆に流体が固体に及ぼす力は $-\bm{R}$ です。この符号の扱いは初学者がつまずきやすい点なので注意してください。
ここまでで一般的な運動量方程式を導出しました。次に、最も多く使われる定常流の場合に簡略化しましょう。
定常流の運動量保存式
定常流への簡略化
工学的な多くの問題では、流れが時間的に変化しない定常流を扱います。定常流では速度場や密度場が時間に依存しないため、検査体積内の運動量の時間変化率はゼロになります。
$$ \frac{\partial}{\partial t}\int_{V_{CV}} \rho \bm{v} \, dV = 0 \quad (\text{定常流}) $$
したがって、定常流の運動量方程式は次のように簡略化されます。
$$ \bm{F}_{\text{ext}} = \oint_{CS} \rho \bm{v}(\bm{v} \cdot \bm{n}) \, dA $$
つまり、定常流では検査面を通じて正味流出する運動量フラックスが、検査体積に作用する全外力に等しいのです。
離散的な入口・出口がある場合
さらに、多くの工学的応用では、流体は有限個の入口と出口を通じて検査体積を出入りします。そして、各入口・出口では速度が断面内でほぼ一様であると仮定できる場合が多いです。このとき、面積分は各入口・出口での離散的な和に置き換えられます。
連続の式より、定常流では質量流量 $\dot{m} = \rho A V$（$A$ は断面積、$V$ は断面平均速度）が保存されます。入口の質量流量を $\dot{m}_{\text{in}}$、出口の質量流量を $\dot{m}_{\text{out}}$ とし、入口が1つ、出口が1つの最も単純な場合を考えると、連続の式から $\dot{m}_{\text{in}} = \dot{m}_{\text{out}} = \dot{m}$ です。
このとき運動量方程式は、各成分（$x$, $y$, $z$）ごとに次のように書けます。
$$ \sum F_x = \dot{m}(v_{x,\text{out}} – v_{x,\text{in}}) $$
$$ \sum F_y = \dot{m}(v_{y,\text{out}} – v_{y,\text{in}}) $$
$$ \sum F_z = \dot{m}(v_{z,\text{out}} – v_{z,\text{in}}) $$
この式は非常に直感的です。力の合計は「出口の運動量フラックス $-$ 入口の運動量フラックス」、すなわち運動量の変化率に等しいのです。これは質点に対する運動量保存則を流体系に拡張したものに他なりません。
ここまでで理論的な枠組みが整いました。いよいよ具体的な応用問題に入りましょう。まず、最も代表的かつ重要な応用であるジェットの推力計算から始めます。
ジェットエンジンの推力計算
問題設定
ジェットエンジンは、前方から空気を取り込み、燃料と混合して燃焼させ、高速の排気ガスを後方に噴出します。この排気ガスの運動量変化がエンジンの推力を生み出します。
エンジンを囲む検査体積を設定しましょう。入口（エンジン前方の空気取り入れ口）では、空気が速度 $v_{\text{in}}$、圧力 $p_{\text{in}}$、密度 $\rho_{\text{in}}$、断面積 $A_{\text{in}}$ で流入します。出口（排気ノズル）では、燃焼ガスが速度 $v_{\text{out}}$、圧力 $p_{\text{out}}$、密度 $\rho_{\text{out}}$、断面積 $A_{\text{out}}$ で流出します。エンジンには燃料が質量流量 $\dot{m}_f$ で供給されます。
推力の導出
検査体積に対して $x$ 方向（推力の方向）の運動量方程式を適用します。ここでは、外力として圧力差とエンジンが受ける推力 $T$ を考えます。
入口の質量流量を $\dot{m}_a = \rho_{\text{in}} A_{\text{in}} v_{\text{in}}$（空気）、出口の質量流量を $\dot{m}_a + \dot{m}_f$（空気＋燃料）とすると、運動量方程式は次のようになります。
$$ T + p_{\text{in}}A_{\text{in}} – p_{\text{out}}A_{\text{out}} = (\dot{m}_a + \dot{m}_f)v_{\text{out}} – \dot{m}_a v_{\text{in}} $$
左辺第1項がエンジンの推力、第2項と第3項が入口・出口の圧力による力です。右辺は出口と入口の運動量フラックスの差です。
$T$ について解くと、ジェットエンジンの推力公式が得られます。
$$ T = (\dot{m}_a + \dot{m}_f)v_{\text{out}} – \dot{m}_a v_{\text{in}} + (p_{\text{out}} – p_{\text{in}})A_{\text{out}} $$
この式の第1項と第2項は運動量の変化に基づく推力（運動量推力）、第3項は圧力差による推力（圧力推力）です。多くの場合、排気ノズルが十分に膨張して出口圧力が大気圧と等しい（$p_{\text{out}} = p_{\text{in}} = p_{\text{atm}}$、完全膨張ノズル）場合、圧力推力の項はゼロになり、推力は運動量項のみで決まります。
$$ T = (\dot{m}_a + \dot{m}_f)v_{\text{out}} – \dot{m}_a v_{\text{in}} \quad (\text{完全膨張ノズル}) $$
さらに、燃料の質量流量が空気に比べて小さい（$\dot{m}_f \ll \dot{m}_a$）場合は、次のように簡略化されます。
$$ T \approx \dot{m}_a (v_{\text{out}} – v_{\text{in}}) $$
この結果は物理的に明快です。推力は「排気の速度と吸気の速度の差」に質量流量を掛けたもの、すなわち運動量の変化率そのものなのです。
ロケットとの違い
ロケットエンジンでは、吸気がありません（$v_{\text{in}} = 0$、$\dot{m}_a = 0$）。したがって推力は次のように表されます。
$$ T_{\text{rocket}} = \dot{m}_f v_{\text{out}} + (p_{\text{out}} – p_{\text{atm}})A_{\text{out}} $$
これはロケット方程式の出発点となる式です。宇宙空間では $p_{\text{atm}} = 0$ なので、ノズル出口圧力も推力に寄与します。
では、Python で様々な条件でのジェット推力を計算し、パラメータ依存性を可視化してみましょう。
Python: ジェット推力の計算
排気速度と空気質量流量がジェット推力にどのように影響するかを計算し、可視化します。
import numpy as np
import matplotlib.pyplot as plt
# --- ジェットエンジン推力の計算 ---
# パラメータ設定
v_in = 250.0 # 飛行速度 (m/s) ≈ 900 km/h
m_dot_a = 80.0 # 空気質量流量 (kg/s)
m_dot_f = 2.0 # 燃料質量流量 (kg/s)
# 排気速度を変化させる
v_out = np.linspace(300, 900, 200) # 排気速度 (m/s)
# 推力の計算（完全膨張ノズル: 圧力推力 = 0）
T_full = (m_dot_a + m_dot_f) * v_out - m_dot_a * v_in
T_approx = m_dot_a * (v_out - v_in) # 燃料流量無視の近似
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# (a) 推力 vs 排気速度
axes[0].plot(v_out, T_full / 1000, 'b-', linewidth=2, label='Full formula')
axes[0].plot(v_out, T_approx / 1000, 'r--', linewidth=2, label='Approx ($\\dot{m}_f \\approx 0$)')
axes[0].set_xlabel('Exhaust velocity $v_{out}$ (m/s)', fontsize=12)
axes[0].set_ylabel('Thrust $T$ (kN)', fontsize=12)
axes[0].set_title('Jet Thrust vs Exhaust Velocity', fontsize=13)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
# (b) 飛行速度の影響
v_out_fixed = 600.0 # 固定排気速度
v_flight = np.linspace(0, 350, 200) # 飛行速度
T_vs_flight = (m_dot_a + m_dot_f) * v_out_fixed - m_dot_a * v_flight
axes[1].plot(v_flight, T_vs_flight / 1000, 'g-', linewidth=2)
axes[1].set_xlabel('Flight velocity $v_{in}$ (m/s)', fontsize=12)
axes[1].set_ylabel('Thrust $T$ (kN)', fontsize=12)
axes[1].set_title('Thrust vs Flight Velocity', fontsize=13)
axes[1].grid(True, alpha=0.3)
axes[1].axhline(y=0, color='k', linewidth=0.5)
plt.tight_layout()
plt.savefig('jet_thrust.png', dpi=150, bbox_inches='tight')
plt.show()
# 数値例の出力
v_out_example = 600.0
T_example = (m_dot_a + m_dot_f) * v_out_example - m_dot_a * v_in
print(f"=== ジェット推力の数値例 ===")
print(f"空気質量流量: {m_dot_a} kg/s")
print(f"燃料質量流量: {m_dot_f} kg/s")
print(f"飛行速度: {v_in} m/s")
print(f"排気速度: {v_out_example} m/s")
print(f"推力: {T_example:.1f} N = {T_example/1000:.2f} kN")
左のグラフから、推力は排気速度に対してほぼ線形に増加することがわかります。また、完全な式（青実線）と燃料流量を無視した近似（赤破線）の差はごく小さく、空燃比が大きい一般的なジェットエンジン（空気流量が燃料の40倍程度）では、この近似が十分に実用的であることが確認できます。
右のグラフは、飛行速度が増加すると推力が減少する様子を示しています。これは直感と一致します。飛行速度が排気速度に近づくと、入口と出口の速度差が小さくなり、運動量の変化が減るからです。飛行速度がある値を超えると推力がゼロになる点があり、これがそのエンジンの理論的な最高速度に相当します。
次に、管路の曲がり部に作用する力の計算に移りましょう。配管設計では、曲がり部のサポートに必要な力を見積もることが非常に重要です。
管路の曲がり部に作用する力
問題設定
管路が角度 $\theta$ だけ曲がっている部分を考えます。流体が管路の曲がり部を通過する際、流れの方向が変わるため、運動量が変化し、管壁に力が作用します。
入口では流れが $x$ 方向正、出口では流れの方向が $x$ 軸から角度 $\theta$ だけ傾いているとします。入口の断面積を $A_1$、圧力を $p_1$、速度を $v_1$、出口のそれぞれを $A_2$、$p_2$、$v_2$ とします。
検査体積として管路の曲がり部を囲む領域を取ります。管壁が流体に及ぼす力を $\bm{R} = (R_x, R_y)$ とします。
力の導出
$x$ 方向と $y$ 方向それぞれに運動量方程式を適用します。
連続の式より、定常・非圧縮性流体では $\rho A_1 v_1 = \rho A_2 v_2 = \dot{m}$ です。
入口では速度ベクトルが $(v_1, 0)$、出口では $(v_2 \cos\theta, v_2 \sin\theta)$ です。$x$ 方向の運動量方程式は次のようになります。
$$ R_x + p_1 A_1 – p_2 A_2 \cos\theta = \dot{m}(v_2 \cos\theta – v_1) $$
左辺で、$p_1 A_1$ は入口の圧力が検査体積に押し込む力（$x$ 正方向）、$p_2 A_2 \cos\theta$ は出口の圧力が検査体積を押す力の $x$ 成分（出口の法線が斜めなので $\cos\theta$ がかかる）です。
$R_x$ について解くと、次のようになります。
$$ R_x = \dot{m}(v_2 \cos\theta – v_1) – p_1 A_1 + p_2 A_2 \cos\theta $$
同様に $y$ 方向の運動量方程式は次のとおりです。
$$ R_y – p_2 A_2 \sin\theta = \dot{m}(v_2 \sin\theta – 0) $$
ここで入口の $y$ 方向速度はゼロ、出口の $y$ 方向速度は $v_2\sin\theta$ です。$R_y$ について解くと、次のようになります。
$$ R_y = \dot{m} v_2 \sin\theta + p_2 A_2 \sin\theta $$
管壁が流体に及ぼす力の合力は次のとおりです。
$$ |\bm{R}| = \sqrt{R_x^2 + R_y^2} $$
逆に、流体が管壁に及ぼす力は $-\bm{R}$ です。配管サポートはこの力に耐えられるように設計しなければなりません。
特殊な場合: 90度曲がり
実際の配管で最も多い90度曲がり（$\theta = 90°$）を考えてみましょう。断面積が一定（$A_1 = A_2 = A$）の場合、連続の式から $v_1 = v_2 = v$ です。
$\cos 90° = 0$、$\sin 90° = 1$ を代入すると、次のようになります。
$$ R_x = \dot{m}(0 – v) – p_1 A + 0 = -\dot{m}v – p_1 A $$
$$ R_y = \dot{m}v + p_2 A $$
圧力が入口と出口で等しい（$p_1 = p_2 = p_g$、ゲージ圧）場合は、さらに簡略化されます。
$$ R_x = -(\dot{m}v + p_g A), \quad R_y = \dot{m}v + p_g A $$
合力は次のとおりです。
$$ |\bm{R}| = \sqrt{2} \cdot (\dot{m}v + p_g A) $$
この力の方向は、$x$ 軸から135度の方向（曲がりの外側方向）です。つまり、流体は管路の曲がりの外側に向かって管壁を押すことがわかります。これは日常の経験とも一致します。ホースを曲げるとき、曲がりの外側が膨らもうとする力を感じるのはこのためです。
特殊な場合: 180度曲がり（Uターン）
$\theta = 180°$ の場合（U字管）は、$\cos 180° = -1$、$\sin 180° = 0$ なので、次のようになります。
$$ R_x = \dot{m}(-v – v) – p_1 A – p_2 A = -2\dot{m}v – (p_1 + p_2)A $$
$$ R_y = 0 $$
180度曲がりでは、力は $x$ 方向のみに作用し、その大きさは90度曲がりの場合より大きくなります。運動量が完全に反転するため、運動量変化が最大になるからです。
では、Python で曲がり角度を変化させたときの力の変化を計算し、可視化しましょう。
Python: 管路曲がり部の力の計算
管路の曲がり角度 $\theta$ を 0度から180度まで変化させたときの管壁に作用する力を計算します。
import numpy as np
import matplotlib.pyplot as plt
# --- 管路曲がり部に作用する力の計算 ---
# パラメータ設定
rho = 1000.0 # 水の密度 (kg/m^3)
D = 0.10 # 管直径 (m)
A = np.pi * D**2 / 4 # 断面積
v = 5.0 # 流速 (m/s)
p_g = 200e3 # ゲージ圧 (Pa) = 200 kPa
m_dot = rho * A * v # 質量流量
# 曲がり角度を変化させる
theta_deg = np.linspace(0, 180, 200)
theta_rad = np.radians(theta_deg)
# 断面積一定の場合（v1 = v2 = v, p1 = p2 = p_g）
Rx = m_dot * v * (np.cos(theta_rad) - 1) - p_g * A * (1 - np.cos(theta_rad))
Ry = m_dot * v * np.sin(theta_rad) + p_g * A * np.sin(theta_rad)
R_mag = np.sqrt(Rx**2 + Ry**2)
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# (a) 力の各成分と合力
axes[0].plot(theta_deg, Rx / 1000, 'b-', linewidth=2, label='$R_x$ (kN)')
axes[0].plot(theta_deg, Ry / 1000, 'r-', linewidth=2, label='$R_y$ (kN)')
axes[0].plot(theta_deg, R_mag / 1000, 'k--', linewidth=2, label='$|R|$ (kN)')
axes[0].set_xlabel('Bend angle $\\theta$ (degrees)', fontsize=12)
axes[0].set_ylabel('Force (kN)', fontsize=12)
axes[0].set_title('Forces on Pipe Bend vs Angle', fontsize=13)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[0].axhline(y=0, color='k', linewidth=0.5)
# (b) 力のベクトル図（代表的な角度）
angles_show = [30, 60, 90, 120, 150, 180]
ax2 = axes[1]
for theta_d in angles_show:
theta_r = np.radians(theta_d)
rx = m_dot * v * (np.cos(theta_r) - 1) - p_g * A * (1 - np.cos(theta_r))
ry = m_dot * v * np.sin(theta_r) + p_g * A * np.sin(theta_r)
# 流体が管壁に及ぼす力は -R
fx, fy = -rx, -ry
scale = 1000 # kN単位
ax2.arrow(theta_d, 0, 0, fy / scale, head_width=3, head_length=0.05,
fc='red', ec='red', alpha=0.7)
ax2.arrow(theta_d, 0, 0, fx / scale, head_width=3, head_length=0.05,
fc='blue', ec='blue', alpha=0.4)
ax2.set_xlabel('Bend angle $\\theta$ (degrees)', fontsize=12)
ax2.set_ylabel('Force on pipe wall (kN)', fontsize=12)
ax2.set_title('Force Components at Selected Angles', fontsize=13)
ax2.grid(True, alpha=0.3)
ax2.axhline(y=0, color='k', linewidth=0.5)
plt.tight_layout()
plt.savefig('pipe_bend_forces.png', dpi=150, bbox_inches='tight')
plt.show()
# 数値結果の出力
print(f"=== 管路曲がり部の力（数値例）===")
print(f"管直径: {D*100:.0f} cm, 流速: {v} m/s, ゲージ圧: {p_g/1000:.0f} kPa")
print(f"質量流量: {m_dot:.2f} kg/s")
for theta_d in [45, 90, 135, 180]:
theta_r = np.radians(theta_d)
rx = m_dot * v * (np.cos(theta_r) - 1) - p_g * A * (1 - np.cos(theta_r))
ry = m_dot * v * np.sin(theta_r) + p_g * A * np.sin(theta_r)
r_mag = np.sqrt(rx**2 + ry**2)
print(f" θ={theta_d:3d}°: Rx={rx:.1f} N, Ry={ry:.1f} N, |R|={r_mag:.1f} N ({r_mag/1000:.2f} kN)")
左のグラフから、曲がり角度 $\theta$ が大きくなるほど合力 $|\bm{R}|$ が増大することが明確に読み取れます。$\theta = 0°$ では力はゼロ（直管）、$\theta = 180°$（U字管）では力が最大となります。$R_x$ は常に負（入口方向に押し戻される力）、$R_y$ は $0° < \theta < 180°$ で正（曲がりの外側方向）です。
数値的にも、直径10cmの管に水が毎秒5mで流れ、ゲージ圧200kPaという一般的な条件でも、90度曲がりで数kNの力が作用することがわかります。これは配管サポート設計において無視できない大きさです。
次に、噴流が平板に衝突する問題を考えます。これは運動量の法則の典型的な応用であり、タービンのバケットへの噴流衝突や消防のノズル設計にも関連する重要な問題です。
衝突噴流の力学
問題設定
断面積 $A_j$、速度 $v_j$ の噴流（ジェット）が、傾斜角 $\alpha$ の平板に衝突する場合を考えます。噴流は平板に衝突した後、平板に沿って上下（または左右）に分かれて流れていきます。
この問題で重要な仮定は以下のとおりです。
- 流体は非圧縮性で、密度 $\rho$ は一定
- 摩擦は無視できる（ベルヌーイの定理が適用可能）
- 噴流と周囲の圧力はともに大気圧（ゲージ圧ゼロ）
- 重力の影響は無視する
垂直衝突（$\alpha = 90°$）
まず、最も単純な垂直衝突から考えましょう。噴流が速度 $v_j$ で平板に垂直に衝突します。
検査体積として噴流が平板に衝突する部分を囲む領域を設定します。入口では噴流が平板に向かう方向（$x$ 方向とします）に速度 $v_j$ で流入し、出口では平板に沿って（$y$ 方向）上下対称に流出します。
$x$ 方向の運動量方程式を適用します。出口での $x$ 方向の速度成分はゼロ（平板に沿って流れるため）、入口での $x$ 方向速度は $v_j$ です。圧力はすべてゲージ圧ゼロなので、運動量方程式は次のようになります。
$$ -F_x = 0 – \dot{m}v_j = -\rho A_j v_j^2 $$
ここで $F_x$ は平板が噴流に及ぼす力（$x$ 負方向、噴流を止める方向）、$\dot{m} = \rho A_j v_j$ は質量流量です。符号を整理すると、噴流が平板に及ぼす力は次のとおりです。
$$ F = \rho A_j v_j^2 $$
この力は速度の2乗に比例することに注目してください。速度を2倍にすると力は4倍になります。これは流体の性質に起因する非線形性であり、流体力学の多くの現象に見られる特徴です。
斜め衝突（一般の角度 $\alpha$）
噴流が平板に角度 $\alpha$ で斜めに衝突する場合はどうなるでしょうか。ここで $\alpha$ は噴流の軸と平板のなす角度とします。
噴流の速度を平板に垂直な成分 $v_j \sin\alpha$ と平板に平行な成分 $v_j \cos\alpha$ に分解します。平板が滑らかで摩擦がないとき、平板に平行な方向の運動量は保存されます（平板は平行方向の力を及ぼさないため）。したがって、平板に垂直な方向のみ運動量が変化します。
平板に垂直な方向の運動量方程式を適用すると、平板に及ぼす力の垂直成分は次のようになります。
$$ F_n = \rho A_j v_j^2 \sin\alpha $$
$\alpha = 90°$ のとき $F_n = \rho A_j v_j^2$ となり、垂直衝突の結果と一致することが確認できます。
移動平板への衝突
ペルトン水車のバケットのように、平板が速度 $v_p$ で噴流と同じ方向に移動している場合はどうなるでしょうか。この場合、平板から見た相対速度は $v_j – v_p$ になります。
移動平板に及ぼされる力は次のとおりです。
$$ F = \rho A_j (v_j – v_p)^2 $$
この力が平板に対してなす仕事率（パワー）は次のようになります。
$$ P = F \cdot v_p = \rho A_j (v_j – v_p)^2 v_p $$
パワーが最大になる平板速度を求めるため、$dP/dv_p = 0$ とすると $v_p = v_j/3$ が得られます。このときの最大パワーは次のとおりです。
$$ P_{\max} = \frac{4}{27}\rho A_j v_j^3 $$
一方、噴流が持つ運動エネルギーフラックスは $\frac{1}{2}\rho A_j v_j^3$ ですから、平板1枚の場合の最大効率は $8/27 \approx 29.6\%$ にとどまります。実際のペルトン水車では、バケットがU字形をしていて噴流をほぼ180度転向させるため、効率は大幅に向上します（理論的にはバケットが噴流を180度転向させ、かつ $v_p = v_j/2$ で最大効率100%）。
では Python でこれらの衝突噴流の力学を計算・可視化しましょう。
Python: 衝突噴流の力学
噴流の衝突角度と速度に対する力の変化、および移動平板の効率を可視化します。
import numpy as np
import matplotlib.pyplot as plt
# --- 衝突噴流の力学 ---
rho = 1000.0 # 水の密度 (kg/m^3)
D_j = 0.05 # 噴流直径 (m)
A_j = np.pi * D_j**2 / 4 # 噴流断面積
v_j = 20.0 # 噴流速度 (m/s)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (a) 垂直衝突: 力 vs 噴流速度
v_range = np.linspace(1, 40, 200)
F_normal = rho * A_j * v_range**2
axes[0, 0].plot(v_range, F_normal, 'b-', linewidth=2)
axes[0, 0].set_xlabel('Jet velocity $v_j$ (m/s)', fontsize=12)
axes[0, 0].set_ylabel('Force $F$ (N)', fontsize=12)
axes[0, 0].set_title('Normal Impact: Force vs Velocity', fontsize=13)
axes[0, 0].grid(True, alpha=0.3)
# (b) 斜め衝突: 力 vs 衝突角度
alpha_deg = np.linspace(0, 90, 200)
alpha_rad = np.radians(alpha_deg)
F_oblique = rho * A_j * v_j**2 * np.sin(alpha_rad)
axes[0, 1].plot(alpha_deg, F_oblique, 'r-', linewidth=2)
axes[0, 1].set_xlabel('Impact angle $\\alpha$ (degrees)', fontsize=12)
axes[0, 1].set_ylabel('Normal force $F_n$ (N)', fontsize=12)
axes[0, 1].set_title(f'Oblique Impact: Force vs Angle ($v_j={v_j}$ m/s)', fontsize=13)
axes[0, 1].grid(True, alpha=0.3)
# (c) 移動平板: 力とパワー vs 平板速度
v_p = np.linspace(0, v_j, 200)
F_moving = rho * A_j * (v_j - v_p)**2
P_moving = F_moving * v_p
P_max = 4/27 * rho * A_j * v_j**3
ax_c = axes[1, 0]
ax_c2 = ax_c.twinx()
line1, = ax_c.plot(v_p / v_j, F_moving, 'b-', linewidth=2, label='Force $F$')
line2, = ax_c2.plot(v_p / v_j, P_moving, 'r-', linewidth=2, label='Power $P$')
ax_c2.axhline(y=P_max, color='r', linestyle='--', alpha=0.5, label=f'$P_{{max}}$={P_max:.1f} W')
ax_c2.axvline(x=1/3, color='gray', linestyle=':', alpha=0.5)
ax_c.set_xlabel('Plate velocity ratio $v_p / v_j$', fontsize=12)
ax_c.set_ylabel('Force $F$ (N)', fontsize=12, color='blue')
ax_c2.set_ylabel('Power $P$ (W)', fontsize=12, color='red')
ax_c.set_title('Moving Plate: Force and Power', fontsize=13)
lines = [line1, line2]
labels = [l.get_label() for l in lines]
ax_c.legend(lines, labels, fontsize=11, loc='upper right')
ax_c.grid(True, alpha=0.3)
# (d) 効率 vs 平板速度
eta = P_moving / (0.5 * rho * A_j * v_j**3)
axes[1, 1].plot(v_p / v_j, eta * 100, 'g-', linewidth=2)
axes[1, 1].axvline(x=1/3, color='gray', linestyle=':', alpha=0.5, label='$v_p/v_j = 1/3$')
axes[1, 1].axhline(y=8/27*100, color='gray', linestyle='--', alpha=0.5, label=f'$\\eta_{{max}}$={8/27*100:.1f}%')
axes[1, 1].set_xlabel('Plate velocity ratio $v_p / v_j$', fontsize=12)
axes[1, 1].set_ylabel('Efficiency $\\eta$ (%)', fontsize=12)
axes[1, 1].set_title('Moving Plate: Efficiency', fontsize=13)
axes[1, 1].legend(fontsize=11)
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('jet_impingement.png', dpi=150, bbox_inches='tight')
plt.show()
# 数値例の出力
print(f"=== 衝突噴流の数値例 ===")
print(f"噴流直径: {D_j*100:.1f} cm, 速度: {v_j} m/s")
print(f"質量流量: {rho * A_j * v_j:.2f} kg/s")
print(f"垂直衝突の力: {rho * A_j * v_j**2:.1f} N")
print(f"最大パワー (v_p = v_j/3): {P_max:.1f} W")
print(f"最大効率: {8/27*100:.1f}%")
4つのグラフから重要な知見が読み取れます。左上のグラフでは、力が速度の2乗に比例して増加する様子が明確に見えます。噴流速度を2倍にすると力は4倍になるという非線形性は、消防ホースのノズル速度が大きいときに消防士が受ける反力の急増を説明します。
右上のグラフでは、衝突角度に対する垂直力が $\sin\alpha$ に従って変化し、垂直衝突（$\alpha = 90°$）で最大になることがわかります。
左下と右下のグラフは特に工学的に重要です。移動平板のパワーは $v_p/v_j = 1/3$ で最大になり、最大効率は約29.6%にとどまります。これは平板1枚で噴流のエネルギーを回収する場合の理論的な限界を示しています。実際の水力タービンではバケット形状の工夫によってこの限界を大幅に超えることができます。
最後に、流体の運動量の法則を翼の揚力・抗力の評価に応用する方法を見ていきましょう。
翼周りの検査体積解析
揚力と抗力の運動量的解釈
飛行機の翼はなぜ揚力を発生するのでしょうか。翼の形状が上面と下面で異なるため、流れが曲げられ、翼の周りで下向きの運動量が生み出されます。ニュートンの運動法則により、流体の運動量が変化するなら、翼は流体に力を及ぼしているはずです。その反作用が揚力です。
この考え方は、ベルヌーイの定理に基づく圧力差の説明と矛盾するものではありません。ベルヌーイの定理は翼面上の圧力分布から揚力を説明しますが、運動量の法則は検査体積の境界での速度場から同じ揚力を導きます。両者は同じ物理現象の異なる見方です。
検査体積の設定
翼を囲む大きな矩形の検査体積を設定します。流れが左から右（$x$ 方向）に一様流速度 $U_\infty$ で入ってくるとします。検査体積は十分に大きく取り、上面と下面では流れがほぼ乱されていないとします。
検査体積の4つの面を考えます。
- 入口面（左面）: 一様流 $U_\infty$ が $x$ 方向に流入
- 出口面（右面）: 翼の影響で速度分布が変化している
- 上面・下面: 流れが多少曲げられているため、$y$ 方向の速度成分を持つ
抗力の計算
$x$ 方向の運動量方程式から抗力 $D$ を求めます。翼の後方（出口面）では、翼の後流（ウェイク）によって流速が一様流より低下しています。この速度欠損を $\Delta v(y)$ とすると、出口面での速度は $u(y) = U_\infty – \Delta v(y)$ です。
入口と出口の $x$ 方向運動量フラックスの差から抗力が求まります。入口面の高さを $H$（検査体積の高さ）、出口面の速度分布を $u(y)$ として、単位スパン長さあたりの抗力は次のようになります。
$$ D = \int_{-H/2}^{H/2} \rho u(y) \left[U_\infty – u(y)\right] dy $$
この式は、出口面での速度欠損が大きいほど抗力が大きいことを表しています。実際の風洞実験では、翼の後方で速度分布を測定し、この式から抗力を算出するウェイクサーベイ法が広く使われています。速度分布が既知であれば翼面の圧力分布を知る必要がないという点で、この方法は非常に実用的です。
揚力の計算
$y$ 方向の運動量方程式から揚力 $L$ を求めます。翼は流れを下向きに曲げるため、検査体積の上面と下面を通じて $y$ 方向の運動量が流出します。
翼の循環（circulation）$\Gamma$ を用いると、翼から十分離れた位置での流れの曲がりは循環に関連づけられます。クッタ・ジュコフスキーの定理により、単位スパンあたりの揚力は次のとおりです。
$$ L = \rho U_\infty \Gamma $$
この式は驚くほどシンプルですが、非常に深い内容を持っています。揚力は一様流速度と循環の積に密度を掛けたものであり、翼の形状の詳細には直接依存しません（循環の値を通じて間接的に依存します）。
検査体積アプローチでは、揚力は翼が流体に与える下向き運動量（ダウンウォッシュ）の反作用として理解できます。翼の上方では流れが下向きに曲げられ、下方でも下向きの流れが誘起されます。この下向き運動量のフラックスの積分が揚力に等しくなるのです。
Python: 翼周りの検査体積解析
翼の後流の速度分布から抗力を計算し、揚力との関係を可視化します。後流の速度分布はガウス分布で近似します。
import numpy as np
import matplotlib.pyplot as plt
# --- 翼周りの検査体積解析 ---
U_inf = 50.0 # 一様流速度 (m/s)
rho = 1.225 # 空気密度 (kg/m^3)
chord = 1.0 # 翼弦長 (m)
# 後流の速度分布モデル（ガウス分布近似）
# ウェイク幅と速度欠損の最大値をパラメータとする
def wake_profile(y, delta_v_max, wake_width):
"""ガウス分布で後流の速度欠損をモデル化"""
return delta_v_max * np.exp(-y**2 / (2 * wake_width**2))
# パラメータ
delta_v_max = 10.0 # 最大速度欠損 (m/s)
wake_width = 0.05 # ウェイク幅 (m)
y = np.linspace(-0.5, 0.5, 1000)
delta_v = wake_profile(y, delta_v_max, wake_width)
u_wake = U_inf - delta_v
# 抗力の計算（ウェイクサーベイ法）
# D = ∫ ρ u(y) [U∞ - u(y)] dy
D_integrand = rho * u_wake * delta_v
D = np.trapz(D_integrand, y)
# 抗力係数
q_inf = 0.5 * rho * U_inf**2 # 動圧
Cd = D / (q_inf * chord)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (a) 後流速度分布
axes[0, 0].plot(u_wake, y, 'b-', linewidth=2)
axes[0, 0].axvline(x=U_inf, color='gray', linestyle='--', alpha=0.5, label=f'$U_\\infty$={U_inf} m/s')
axes[0, 0].fill_betweenx(y, u_wake, U_inf, alpha=0.2, color='blue')
axes[0, 0].set_xlabel('Velocity $u(y)$ (m/s)', fontsize=12)
axes[0, 0].set_ylabel('$y$ (m)', fontsize=12)
axes[0, 0].set_title('Wake Velocity Profile', fontsize=13)
axes[0, 0].legend(fontsize=11)
axes[0, 0].grid(True, alpha=0.3)
# (b) 抗力の被積分関数
axes[0, 1].plot(D_integrand, y, 'r-', linewidth=2)
axes[0, 1].fill_betweenx(y, 0, D_integrand, alpha=0.2, color='red')
axes[0, 1].set_xlabel('Drag integrand $\\rho u \\Delta v$ (N/m$^2$)', fontsize=12)
axes[0, 1].set_ylabel('$y$ (m)', fontsize=12)
axes[0, 1].set_title(f'Drag Integrand (D = {D:.2f} N/m)', fontsize=13)
axes[0, 1].grid(True, alpha=0.3)
# (c) 速度欠損の大きさ vs 抗力
dv_range = np.linspace(1, 25, 100)
D_vs_dv = []
for dv in dv_range:
delta_v_temp = wake_profile(y, dv, wake_width)
u_temp = U_inf - delta_v_temp
D_temp = np.trapz(rho * u_temp * delta_v_temp, y)
D_vs_dv.append(D_temp)
D_vs_dv = np.array(D_vs_dv)
axes[1, 0].plot(dv_range, D_vs_dv, 'g-', linewidth=2)
axes[1, 0].set_xlabel('Max velocity deficit $\\Delta v_{max}$ (m/s)', fontsize=12)
axes[1, 0].set_ylabel('Drag per unit span $D$ (N/m)', fontsize=12)
axes[1, 0].set_title('Drag vs Wake Deficit', fontsize=13)
axes[1, 0].grid(True, alpha=0.3)
# (d) 揚力と循環の関係（クッタ・ジュコフスキー）
Gamma_range = np.linspace(0, 15, 200)
L_range = rho * U_inf * Gamma_range
# 揚力係数
Cl_range = L_range / (q_inf * chord)
ax_d = axes[1, 1]
ax_d2 = ax_d.twinx()
line1, = ax_d.plot(Gamma_range, L_range, 'm-', linewidth=2, label='Lift $L$')
line2, = ax_d2.plot(Gamma_range, Cl_range, 'c--', linewidth=2, label='$C_L$')
ax_d.set_xlabel('Circulation $\\Gamma$ (m$^2$/s)', fontsize=12)
ax_d.set_ylabel('Lift per unit span $L$ (N/m)', fontsize=12, color='m')
ax_d2.set_ylabel('Lift coefficient $C_L$', fontsize=12, color='c')
ax_d.set_title('Kutta-Joukowski: Lift vs Circulation', fontsize=13)
lines = [line1, line2]
labels = [l.get_label() for l in lines]
ax_d.legend(lines, labels, fontsize=11)
ax_d.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('wing_analysis.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"=== 翼の検査体積解析 ===")
print(f"一様流速度: {U_inf} m/s")
print(f"抗力 (per unit span): {D:.2f} N/m")
print(f"抗力係数 Cd: {Cd:.5f}")
print(f"動圧: {q_inf:.1f} Pa")
左上のグラフは翼の後方における速度分布を示しています。一様流速度 $U_\infty = 50$ m/s から中心部で速度が低下しているのが後流（ウェイク）です。青く塗られた領域が運動量の欠損を表しており、この面積が抗力に対応します。
右上のグラフは抗力の被積分関数 $\rho u(y)\Delta v(y)$ を示しており、この分布を $y$ 方向に積分したものが単位スパンあたりの抗力です。被積分関数がウェイクの中心に集中していることから、抗力は後流の狭い領域で主に生じていることがわかります。
左下のグラフは最大速度欠損と抗力の関係です。速度欠損が小さい範囲ではほぼ線形に見えますが、欠損が大きくなると $\Delta v^2$ の効果で増加率が鈍化します（$u = U_\infty – \Delta v$ が小さくなるため）。
右下のグラフはクッタ・ジュコフスキーの定理に基づく揚力と循環の線形関係を示しています。循環が大きいほど揚力が大きく、この関係は密度と一様流速度にのみ依存するシンプルなものです。
運動量補正係数
速度分布の不均一性の扱い
ここまでの解析では、管路の各断面で速度が一様であると仮定してきました。しかし実際の管内流れでは、壁面近くの速度はゼロ（滑りなし条件）、中心部で最大という速度分布を持ちます。
速度分布が一様でない場合、断面を通過する運動量フラックスは $\dot{m}V$（$V$ は断面平均速度）とは一致しません。これを補正するために運動量補正係数 $\beta$ を導入します。
$$ \beta = \frac{1}{A}\int_A \left(\frac{v}{V}\right)^2 dA $$
ここで $v$ は断面内の局所速度、$V = Q/A$ は断面平均速度（$Q$ は体積流量）です。$\beta$ は常に1以上であり、速度分布が一様なら $\beta = 1$ です。
代表的な値
層流の場合、管内の速度分布は放物線状です。円管の場合は次のように計算できます。
$$ v(r) = 2V\left(1 – \frac{r^2}{R^2}\right) $$
ここで $R$ は管の半径です。この分布を運動量補正係数の定義に代入すると、次のようになります。
$$ \beta_{\text{laminar}} = \frac{1}{\pi R^2}\int_0^{2\pi}\int_0^R \left(\frac{2V(1-r^2/R^2)}{V}\right)^2 r\,dr\,d\theta $$
$u/V = 2(1-r^2/R^2)$ と $\xi = r/R$ の置換を行い積分を実行すると、次の結果を得ます。
$$ \beta_{\text{laminar}} = 2\int_0^1 4(1-\xi^2)^2 \xi \, d\xi = 8\int_0^1 (1-2\xi^2+\xi^4)\xi\,d\xi = 8\left[\frac{1}{2}-\frac{2}{4}+\frac{1}{6}\right] = \frac{4}{3} $$
つまり、層流では $\beta = 4/3$ であり、一様速度と仮定した場合に比べて運動量フラックスが33%大きくなります。
乱流の場合、速度分布はべき乗則 $v/V_{\max} = (1-r/R)^{1/n}$ で近似され、$n$ はレイノルズ数に依存しますが、典型的には $n \approx 7$ です。この場合 $\beta \approx 1.02$ となり、速度分布が層流よりはるかに一様に近いため、$\beta \approx 1$ の近似が良い精度で成り立ちます。
運動量方程式への組み込み
運動量補正係数を考慮すると、定常流の運動量方程式は次のように修正されます。
$$ \sum \bm{F} = \dot{m}(\beta_{\text{out}} \bm{V}_{\text{out}} – \beta_{\text{in}} \bm{V}_{\text{in}}) $$
ここで $\bm{V}$ は断面平均速度ベクトルです。高精度が必要な解析では $\beta$ の補正が重要ですが、乱流（$\beta \approx 1$）の場合は多くの実用問題で補正を省略できます。
運動量補正係数の理解は、理論と実験を結びつける際に重要です。次節ではPythonを使って層流と乱流の速度分布から運動量補正係数を直接計算し、上記の理論値と比較してみましょう。
Python: 運動量補正係数の計算
層流（放物線分布）と乱流（べき乗則分布）の速度分布から運動量補正係数を数値的に計算します。
import numpy as np
import matplotlib.pyplot as plt
# --- 運動量補正係数の計算 ---
R = 0.05 # 管半径 (m)
N = 500 # 積分の分割数
r = np.linspace(0, R, N)
xi = r / R # 無次元半径
# 層流: 放物線分布 v(r) = 2V(1 - r^2/R^2)
# 平均速度 V で規格化
v_laminar = 2 * (1 - xi**2) # v/V
# 乱流: べき乗則 v(r)/V_max = (1 - r/R)^(1/n)
# V/V_max = 2n^2 / ((n+1)(2n+1)) なので v/V を計算
n_values = [6, 7, 8, 9, 10]
beta_turb = []
for n in n_values:
V_ratio = 2 * n**2 / ((n + 1) * (2 * n + 1)) # V / V_max
v_over_Vmax = (1 - xi)**(1/n)
v_over_Vmax[0] = 1.0 # r=0 での値を修正
v_over_V = v_over_Vmax / V_ratio # v/V
# β = (1/A) ∫ (v/V)^2 dA = 2 ∫_0^1 (v/V)^2 ξ dξ
integrand = v_over_V**2 * xi
beta = 2 * np.trapz(integrand, xi)
beta_turb.append(beta)
# 層流の β を計算
integrand_lam = v_laminar**2 * xi
beta_lam = 2 * np.trapz(integrand_lam, xi)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# (a) 速度分布の比較
axes[0].plot(v_laminar, xi, 'b-', linewidth=2, label=f'Laminar ($\\beta$={beta_lam:.4f})')
for i, n in enumerate([7]):
V_ratio = 2 * n**2 / ((n + 1) * (2 * n + 1))
v_over_Vmax = (1 - xi)**(1/n)
v_over_V = v_over_Vmax / V_ratio
axes[0].plot(v_over_V, xi, 'r-', linewidth=2,
label=f'Turbulent n={n} ($\\beta$={beta_turb[n_values.index(n)]:.4f})')
# 一様分布
axes[0].axvline(x=1, color='gray', linestyle='--', alpha=0.5, label='Uniform ($\\beta$=1)')
axes[0].set_xlabel('$v / V$ (normalized velocity)', fontsize=12)
axes[0].set_ylabel('$r / R$ (normalized radius)', fontsize=12)
axes[0].set_title('Velocity Profiles in Pipe Flow', fontsize=13)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim([0, 2.2])
# (b) 乱流のβ vs べき乗指数 n
axes[1].plot(n_values, beta_turb, 'ro-', linewidth=2, markersize=8, label='Turbulent (power law)')
axes[1].axhline(y=4/3, color='b', linestyle='--', linewidth=1.5, label=f'Laminar $\\beta$=4/3={4/3:.4f}')
axes[1].axhline(y=1.0, color='gray', linestyle=':', linewidth=1.5, label='Uniform $\\beta$=1')
axes[1].set_xlabel('Power law exponent $n$', fontsize=12)
axes[1].set_ylabel('Momentum correction factor $\\beta$', fontsize=12)
axes[1].set_title('Momentum Correction Factor', fontsize=13)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('momentum_correction.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"=== 運動量補正係数 ===")
print(f"層流 (放物線分布): β = {beta_lam:.4f} (理論値: {4/3:.4f})")
print(f"乱流 (べき乗則):")
for n, b in zip(n_values, beta_turb):
print(f" n={n}: β = {b:.4f}")
左のグラフは管内の速度分布を比較しています。層流（青線）では放物線状の分布で中心速度が平均速度の2倍になるのに対し、乱流（赤線、$n=7$）では分布がはるかに平坦で、壁面近くを除けば速度がほぼ一様に近いことがわかります。灰色の破線は完全に一様な分布（$\beta = 1$）を示しています。
右のグラフは乱流のべき乗指数 $n$ に対する運動量補正係数 $\beta$ の変化を示しています。$n$ が増加する（レイノルズ数が大きくなる）ほど速度分布が平坦になり、$\beta$ は1に近づきます。層流の $\beta = 4/3 \approx 1.333$ と比べて、乱流の $\beta$ は1.01〜1.06程度であり、多くの工学的問題で $\beta \approx 1$ の仮定が妥当であることが数値的にも確認できます。
数値計算の結果、層流の $\beta$ は理論値 $4/3$ と一致しており、数値積分の精度が十分であることも確認できます。
運動量の法則の応用に関する注意点
検査体積の選び方
運動量の法則を正しく適用するためには、検査体積を適切に設定することが重要です。以下の指針に従いましょう。
1. 求めたい力を含む: 管壁の力を求めたいなら、管壁を検査体積の内部に含めます。翼の力を求めたいなら、翼を検査体積の内部に含めます。
2. 検査面での情報が得やすい場所に設定する: 検査面は速度・圧力がわかる（または合理的に仮定できる）位置に設定します。たとえば、管の入口・出口の断面や、十分遠方の乱されていない流れの位置が適しています。
3. 検査面に垂直な速度成分を確認する: 検査面の法線方向の速度成分 $\bm{v}\cdot\bm{n}$ の符号によって流入・流出が決まります。外向き法線を正とするので、流出なら正、流入なら負です。
符号の取り扱い
運動量の法則を適用する際に最も間違いやすいのが符号の処理です。以下の3点を常に確認してください。
1. 外向き法線の定義: 検査面の法線 $\bm{n}$ は必ず外向きに取ります。入口面では流入方向と逆向き（$\bm{v}\cdot\bm{n} < 0$）、出口面では流出方向と同じ向き（$\bm{v}\cdot\bm{n} > 0$）になります。
2. 力の方向: 管壁が流体に及ぼす力 $\bm{R}$ と、流体が管壁に及ぼす力 $-\bm{R}$ は大きさが等しく方向が逆です。求める対象がどちらかを明確にしましょう。
3. 圧力の扱い: 圧力は検査面に垂直で内向きに作用します。外向き法線を基準とすると $-p\bm{n}$ です。ゲージ圧を使うと大気圧の寄与が自動的にキャンセルされるため、多くの問題でゲージ圧を使うのが便利です。
エネルギーの法則との使い分け
運動量の法則とベルヌーイの定理（エネルギーの法則）は、どちらも流体の問題を解くための強力なツールですが、適用すべき場面が異なります。
運動量の法則が適している場面: – 力を求めたいとき（管路の壁面力、推力、翼の揚力・抗力） – 流れの内部で損失がある場合でも適用可能（損失は内力であり、運動量の法則には影響しない） – 混合、渦、衝撃波など、不可逆過程を含む流れ
ベルヌーイの定理が適している場面: – 速度や圧力を求めたいとき – 流れが定常・非粘性・非圧縮で、損失がない場合 – 流線に沿った関係を利用するとき
多くの工学問題では、連続の式とベルヌーイの定理で未知の速度・圧力を求めてから、運動量の法則で力を計算するという組み合わせが使われます。
微分形との関係
オイラー方程式とナビエ・ストークス方程式
本記事で扱った運動量の法則は積分形ですが、これを微小な検査体積に対して適用し、体積をゼロに近づけると微分形が得られます。
非粘性流体の場合、ベクトル解析の発散の定理を用いて面積分を体積積分に変換し、被積分関数がゼロであることから、オイラー方程式が得られます。
$$ \rho\left(\frac{\partial \bm{v}}{\partial t} + (\bm{v}\cdot\nabla)\bm{v}\right) = -\nabla p + \rho\bm{g} $$
さらに粘性応力を含めると、ナビエ・ストークス方程式になります。
$$ \rho\left(\frac{\partial \bm{v}}{\partial t} + (\bm{v}\cdot\nabla)\bm{v}\right) = -\nabla p + \mu\nabla^2\bm{v} + \rho\bm{g} $$
積分形と微分形は等価な情報を持っていますが、工学的な力の計算では積分形（本記事の内容）が、流れ場の詳細な計算では微分形（CFDで解く形式）が使われます。両者の関係を理解しておくことは、流体力学の全体像を把握する上で重要です。
Python: 総合演習 — 噴流による推進装置の設計
最後に、これまでの理論を統合した総合演習として、水ジェット推進装置の推力とエネルギー効率をPythonで計算します。水ジェット推進は船舶やジェットスキーで実用化されている技術であり、ポンプで取り込んだ水を後方に高速で噴射することで推力を得る仕組みです。
import numpy as np
import matplotlib.pyplot as plt
# --- 水ジェット推進装置の設計解析 ---
rho = 1025.0 # 海水密度 (kg/m^3)
# 設計パラメータ
D_inlet = 0.30 # 取水口直径 (m)
D_nozzle = 0.15 # ノズル直径 (m)
A_inlet = np.pi * D_inlet**2 / 4
A_nozzle = np.pi * D_nozzle**2 / 4
# 船速の範囲
v_ship = np.linspace(1, 20, 200) # 船速 (m/s)
# ポンプの体積流量 Q (m^3/s) を一定とする
Q = 0.5 # 体積流量 (m^3/s)
# 連続の式から入口・出口速度を計算
v_inlet = Q / A_inlet # 取水口では船速の影響で相対速度が変わる
v_nozzle = Q / A_nozzle # ノズル出口速度
# 推力の計算
# 取水口での速度は船速 v_ship（船の座標系で考える）
# ノズルでの速度は v_nozzle
# T = ρQ(v_nozzle - v_ship)
T = rho * Q * (v_nozzle - v_ship)
# 有効パワー（推力 × 船速）
P_useful = T * v_ship
# ポンプが流体に与えるパワー
# ジェットの運動エネルギーフラックス - 取り込み水の運動エネルギーフラックス
P_pump = 0.5 * rho * Q * (v_nozzle**2 - v_ship**2)
# 推進効率
eta_prop = np.where(P_pump > 0, P_useful / P_pump, 0)
# 理論最大効率（フルード効率）
# η_F = 2 / (1 + v_nozzle/v_ship)
eta_froude = 2 / (1 + v_nozzle / v_ship)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (a) 推力 vs 船速
axes[0, 0].plot(v_ship, T / 1000, 'b-', linewidth=2)
axes[0, 0].set_xlabel('Ship velocity $v_{ship}$ (m/s)', fontsize=12)
axes[0, 0].set_ylabel('Thrust $T$ (kN)', fontsize=12)
axes[0, 0].set_title('Water Jet Thrust vs Ship Velocity', fontsize=13)
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].axhline(y=0, color='k', linewidth=0.5)
# (b) パワー
axes[0, 1].plot(v_ship, P_useful / 1000, 'g-', linewidth=2, label='Useful power $Tv_{ship}$')
axes[0, 1].plot(v_ship, P_pump / 1000, 'r--', linewidth=2, label='Pump power')
axes[0, 1].set_xlabel('Ship velocity $v_{ship}$ (m/s)', fontsize=12)
axes[0, 1].set_ylabel('Power (kW)', fontsize=12)
axes[0, 1].set_title('Power vs Ship Velocity', fontsize=13)
axes[0, 1].legend(fontsize=11)
axes[0, 1].grid(True, alpha=0.3)
# (c) 推進効率
valid = P_pump > 0
axes[1, 0].plot(v_ship[valid], eta_prop[valid] * 100, 'b-', linewidth=2, label='Propulsive efficiency')
axes[1, 0].plot(v_ship[valid], eta_froude[valid] * 100, 'r--', linewidth=2, label='Froude efficiency')
axes[1, 0].set_xlabel('Ship velocity $v_{ship}$ (m/s)', fontsize=12)
axes[1, 0].set_ylabel('Efficiency (%)', fontsize=12)
axes[1, 0].set_title('Propulsive Efficiency', fontsize=13)
axes[1, 0].legend(fontsize=11)
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].set_ylim([0, 100])
# (d) ノズル径の影響
D_nozzles = [0.10, 0.12, 0.15, 0.18, 0.20]
colors = plt.cm.viridis(np.linspace(0.2, 0.8, len(D_nozzles)))
for i, Dn in enumerate(D_nozzles):
An = np.pi * Dn**2 / 4
vn = Q / An
T_temp = rho * Q * (vn - v_ship)
axes[1, 1].plot(v_ship, T_temp / 1000, color=colors[i], linewidth=2,
label=f'$D_n$={Dn*100:.0f} cm ($v_n$={vn:.1f} m/s)')
axes[1, 1].set_xlabel('Ship velocity $v_{ship}$ (m/s)', fontsize=12)
axes[1, 1].set_ylabel('Thrust $T$ (kN)', fontsize=12)
axes[1, 1].set_title('Effect of Nozzle Diameter on Thrust', fontsize=13)
axes[1, 1].legend(fontsize=10)
axes[1, 1].grid(True, alpha=0.3)
axes[1, 1].axhline(y=0, color='k', linewidth=0.5)
plt.tight_layout()
plt.savefig('waterjet_propulsion.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"=== 水ジェット推進装置の解析結果 ===")
print(f"ノズル出口速度: {v_nozzle:.1f} m/s")
print(f"体積流量: {Q} m^3/s = {Q*1000:.0f} L/s")
for vs in [5, 10, 15]:
t = rho * Q * (v_nozzle - vs)
p_u = t * vs
p_p = 0.5 * rho * Q * (v_nozzle**2 - vs**2)
eta = p_u / p_p * 100 if p_p > 0 else 0
print(f" 船速 {vs} m/s: 推力={t:.0f} N ({t/1000:.1f} kN), "
f"有効パワー={p_u/1000:.1f} kW, 推進効率={eta:.1f}%")
この総合演習から、水ジェット推進装置の特性が4つのグラフで明確に把握できます。
左上のグラフから、推力は船速の増加とともに線形に減少することがわかります。これはジェットエンジンと同じ傾向であり、入口と出口の速度差が小さくなるためです。推力がゼロになる船速がこの装置の理論最高速度に相当します。
右上のグラフでは、有効パワー（推力と船速の積）がある船速で最大となる一方、ポンプパワーは船速が小さいほど大きいことがわかります。低速域では運動エネルギーが噴流に消費され、推進に有効に使われていないことを意味しています。
左下の効率グラフは特に重要です。推進効率は船速が増加するほど向上し、噴流速度に近づくと効率は100%に近づきます。これはフルードの効率式 $\eta = 2/(1 + v_j/v_s)$ と一致しています。つまり、推力は大きくなるが効率は低い（低速） vs 効率は高いが推力が小さい（高速） というトレードオフがあります。
右下のグラフはノズル径を変えた場合の影響を示しており、ノズル径が小さいほど噴流速度が高くなり、低速域での推力は増加しますが、高速域では推力がゼロになる速度（最高速度）も高くなります。設計者はこのトレードオフを考慮してノズル径を選定する必要があります。
まとめ
本記事では、流体の運動量の法則について、ニュートンの第2法則から出発して系統的に解説しました。
- ニュートンの第2法則の拡張: 質点の $\bm{F} = d(m\bm{v})/dt$ を流体塊に適用し、物質体積における運動量方程式を定式化しました
- レイノルズ輸送定理: 物質体積での記述を検査体積での記述に変換する数学的な橋渡しの定理を導出しました。蓄積項と境界面のフラックスの2つの寄与に分離できることを示しました
- 定常流の運動量保存式: 定常流では蓄積項が消え、$\sum\bm{F} = \dot{m}(\bm{v}_{\text{out}} – \bm{v}_{\text{in}})$ というシンプルな形に帰着することを示しました
- ジェットの推力: 推力が排気と吸気の運動量差に等しいことを導出し、飛行速度との関係を可視化しました
- 管路曲がり部の力: 曲がり角度に依存する壁面力を成分ごとに導出し、90度・180度の特殊ケースを確認しました
- 衝突噴流の力学: 垂直・斜め・移動平板への衝突を解析し、力が速度の2乗に比例することを示しました
- 翼の検査体積解析: ウェイクサーベイ法による抗力算出とクッタ・ジュコフスキーの定理による揚力の運動量的解釈を解説しました
- 運動量補正係数: 速度分布が一様でない場合の補正方法を示し、層流（$\beta = 4/3$）と乱流（$\beta \approx 1$）を比較しました
運動量の法則は、連続の式（質量保存）、ベルヌーイの定理（エネルギー保存）と並んで、流体力学の三大基礎方程式の一つです。これら3つを組み合わせることで、工学的な流体問題の大部分を解くことができます。
次のステップとして、以下の記事も参考にしてください。
- ナビエ・ストークス方程式 — 運動量の法則の微分形を粘性流体に拡張した支配方程式
- 管内流れ — 運動量の法則と粘性を組み合わせた具体的な応用
- ロケット方程式 — 運動量の法則を推進工学に発展させた理論
- レイノルズ数 — 流れの性質を決定する無次元数と運動量補正への影響