ロケットの打ち上げ映像を見ると、ロケットは最初まっすぐ上に飛び立ちますが、しばらくすると徐々に傾き始め、やがて地平線と平行に近い方向に飛んでいきます。なぜ最初から斜めに飛ばないのでしょうか？ そして、いつ、どのくらい傾けるのでしょうか？
この「最初は垂直、途中から徐々に寝ていく」飛行プロファイルが重力ターン（gravity turn）です。重力ターンは、ロケットの推力軸を常に速度ベクトルの方向に合わせることで、重力の作用で自然に飛行経路が曲がっていく打ち上げ方式です。この方法の巧みな点は、ロケットに横向きの力（構造に負担がかかる）をほとんどかけずに、効率的に軌道投入ができることです。
重力ターンを理解すると、以下の応用が見えてきます。
- ロケットの打ち上げ設計: 推力、質量比、打ち上げ角度の最適化
- 軌道投入の効率計算: 重力損失と空力損失を定量化し、必要なΔVを見積もる
- ミッション計画: 異なる軌道（LEO、SSO、GTO）への投入に必要な条件を理解する
本記事の内容
- ロケットはなぜ斜めに飛ぶか — 軌道速度の必要性
- 垂直上昇から軌道投入までの飛行フェーズ
- 重力ターンの運動方程式の導出
- ピッチオーバーキック
- 大気中の空力損失と重力損失
- Pythonで2次元重力ターンのシミュレーション
- 高度・速度プロファイルの可視化と最適な投入角度
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ロケットはなぜ斜めに飛ぶか
軌道に「のる」とは何か
人工衛星を軌道に投入するとは、地球のまわりを周回するのに十分な水平速度を与えることです。高度200 kmの円軌道に必要な軌道速度は約7.78 km/sで、これはほぼ水平方向の速度です。
もしロケットがまっすぐ上に飛び続けたらどうなるでしょうか？ 推進が終わった瞬間、ロケットは上向きの速度を持っていますが、水平方向の速度はゼロです。したがって、ロケットは放物線を描いて地球に戻ってきます — 花火と同じです。
軌道に乗るためには、垂直方向の上昇で高度を稼ぎつつ、水平方向に約7.8 km/sの速度を獲得する必要があります。問題は、この「垂直→水平」の遷移をどうやって効率的に行うかです。
理想的なΔVと実際のΔV
もし大気と重力が存在しなければ（宇宙空間でのマヌーバ）、必要なΔVは軌道速度そのもの（約7.78 km/s）です。しかし現実の打ち上げでは、大気と重力の影響で追加のΔVが必要です。
$$ \Delta V_{\text{total}} = \Delta V_{\text{orbital}} + \Delta V_{\text{gravity loss}} + \Delta V_{\text{drag loss}} – \Delta V_{\text{Earth rotation}} $$
各項を概算すると
| 項目 | 概算値 [km/s] |
|---|---|
| 軌道速度（高度200 km） | 7.78 |
| 重力損失 | 1.0〜1.5 |
| 空力損失 | 0.1〜0.5 |
| 地球自転の利得（赤道、東向き） | -0.46 |
| 合計 | 約8.5〜9.5 |
実際のロケットの総ΔVは約9〜9.5 km/sで、軌道速度の約1.2倍です。重力損失と空力損失の合計が1.1〜2.0 km/sに達するため、これらを最小化する飛行プロファイルの設計が重要です。
重力損失と空力損失がどのように生じるかを理解するために、打ち上げの各フェーズを順に見ていきましょう。
打ち上げの飛行フェーズ
フェーズ1: 垂直上昇（Vertical Rise）
ロケットは発射台からまっすぐ上に飛び立ちます。この段階の目的は3つです。
- 発射台からの安全な離脱: 発射台や周辺施設との衝突を避ける
- 稠密大気の通過: 最も空気が濃い最初の数kmを高速で通過する
- 初期速度の獲得: ピッチオーバーのための基本速度を確保する
垂直上昇は通常10〜20秒程度で、高度数百mから1 km程度まで続きます。
フェーズ2: ピッチオーバーキック（Pitch-Over Kick）
垂直上昇の後、ロケットの姿勢をわずかに（1〜3°程度）傾けます。この小さな傾きがピッチオーバーキックです。打ち上げ軌道全体を決定する極めて重要なパラメータです。
ピッチオーバー角度が大きすぎると大気中での空力荷重が増大し、小さすぎると十分な水平速度が得られません。最適な角度はロケットの推力重量比、大気モデル、目標軌道によって決まります。
フェーズ3: 重力ターン（Gravity Turn）
ピッチオーバーの後、ロケットの推力軸を速度ベクトルの方向に維持します。この状態では、推力は常に速度の方向に作用するため、ロケットに横方向の力がほとんどかかりません。
重力の鉛直下向きの成分が速度ベクトルを「引き下げる」ことで、飛行経路角（水平面となす角度）が徐々に減少していきます。これが「重力ターン」の名前の由来です — 重力が自然にロケットの飛行方向を水平に向けてくれるのです。
この方式の大きな利点は、ロケットの機体に横向きの空力荷重（迎角による荷重）がほとんどかからないことです。大気中での高速飛行時に横方向の力がかかると、機体構造に大きな負担がかかり、最悪の場合は構造破壊につながります。
フェーズ4: 軌道投入（Orbit Insertion）
大気圏を抜けた後、飛行経路角がほぼゼロ（水平飛行）になった時点で、エンジンの残りの燃焼で軌道速度を達成します。多段式ロケットの場合、上段エンジンがこの役割を担います。
各フェーズの物理がイメージできたところで、重力ターンの運動方程式を数学的に導出しましょう。
重力ターンの運動方程式
座標系と変数の定義
地球の表面を基準とした2次元の極座標系を使います。地球中心からロケットまでの距離を $r$、水平方向から測った飛行経路角を $\gamma$（flight path angle）とします。$\gamma = 90°$ が垂直上昇、$\gamma = 0°$ が水平飛行です。
ロケットの速度を $v$、質量を $m(t)$（燃料消費で時間変化）、推力を $T$、大気抵抗を $D$ とします。
速度方向の方程式
推力は速度方向に作用し、重力と抗力は速度方向に反対の成分を持ちます。速度方向のニュートンの運動方程式は
$$ m\dot{v} = T – D – mg\sin\gamma $$
ここで $g = \mu_E/r^2$ は高度に依存する重力加速度です。
右辺の各項を解釈すると
- $T$: 推力によるロケットの加速
- $-D$: 大気抵抗による減速
- $-mg\sin\gamma$: 重力の速度方向成分（飛行経路角が正のとき、重力は減速方向に作用する）
$m$ で割って
$$ \dot{v} = \frac{T}{m} – \frac{D}{m} – g\sin\gamma $$
経路角方向の方程式
速度に垂直な方向（経路角を変化させる方向）の方程式は
$$ v\dot{\gamma} = -g\cos\gamma + \frac{v^2\cos\gamma}{r} $$
右辺第1項 $-g\cos\gamma$ は重力の法線方向成分で、経路角を減少させる方向に作用します。これが重力ターンの本質です — 重力が飛行方向を「水平方向に引き下げる」効果です。
右辺第2項 $v^2\cos\gamma/r$ は遠心力の効果で、経路角を増加させる方向に作用します。軌道速度に達すると、この項が重力項と釣り合い、$\dot{\gamma} = 0$（一定の高度で水平飛行）が実現します。
位置の方程式
ロケットの位置変化は
$$ \dot{r} = v\sin\gamma $$
$$ \dot{\theta} = \frac{v\cos\gamma}{r} $$
ここで $\theta$ は地球中心から見た角度（ダウンレンジ角）です。
質量の方程式
ロケットの質量変化は推進剤の消費で決まります。
$$ \dot{m} = -\frac{T}{I_{\text{sp}} g_0} $$
ここで $I_{\text{sp}}$ は比推力（秒）、$g_0 = 9.80665$ m/s² は標準重力加速度です。
大気抵抗
$$ D = \frac{1}{2}\rho(h) v^2 C_D A $$
$\rho(h)$ は高度 $h = r – R_E$ での大気密度（指数大気モデル）です。
これらの連立微分方程式を数値的に解くことで、ロケットの打ち上げ軌道を計算できます。ただし、その前に重力損失と空力損失の物理をもう少し詳しく理解しておきましょう。
重力損失と空力損失
重力損失
重力損失（gravity loss）は、ロケットが重力に逆らって上昇するために「消費」するΔVです。速度方向の運動方程式の $-g\sin\gamma$ 項を時間積分したものです。
$$ \Delta V_{\text{gravity}} = \int_0^{t_{\text{burn}}} g\sin\gamma \, dt $$
重力損失を最小化するためには
- 燃焼時間を短くする: 高い推力で素早く加速すれば、重力に逆らう時間が短い
- 経路角を早く下げる: $\sin\gamma$ が小さいほど重力損失は少ない
- 推力重量比（TWR）を大きくする: TWRが大きいほど加速が速く、重力損失が少ない
極端な例として、推力が重力とちょうど釣り合う（TWR = 1）ロケットが垂直に上昇する場合を考えます。推力のすべてが重力に対抗するために使われ、加速は一切できません。ΔV全体が重力損失です。逆にTWRが非常に大きければ（瞬間的に加速できれば）、重力損失はほぼゼロになります。
空力損失
空力損失（drag loss）は、大気抵抗に打ち勝つために消費するΔVです。
$$ \Delta V_{\text{drag}} = \int_0^{t_{\text{burn}}} \frac{D}{m} \, dt = \int_0^{t_{\text{burn}}} \frac{\rho v^2 C_D A}{2m} \, dt $$
空力損失を最小化するためには
- 低空で速度を上げすぎない: 密度が高い低空で高速にすると抗力が大きい
- 機体を細くする: 断面積 $A$ を小さくして抗力係数を下げる
- 早く大気圏を抜ける: 高いTWRで素早く高度を稼ぐ
マックスQ（最大動圧）
打ち上げ中に動圧 $q = \frac{1}{2}\rho v^2$ が最大になる瞬間をマックスQ（Max Q）と呼びます。打ち上げ初期は速度が低く $v^2$ が小さいですが、高度が上がると大気密度 $\rho$ が急激に減少します。この2つの効果が拮抗する高度（通常10〜15 km付近）でマックスQが発生します。
マックスQは機体構造への空力荷重が最大になる瞬間で、ロケットの設計においてクリティカルな制約条件です。SpaceXのFalcon 9はマックスQ付近でエンジンの推力を一時的に絞るスロットリングを行い、空力荷重を低減しています。
重力損失と空力損失のトレードオフが理解できたところで、Pythonで重力ターンのシミュレーションを実装しましょう。
Pythonで2次元重力ターンのシミュレーション
ロケットのパラメータ設定
Falcon 9に近い2段式ロケットのパラメータを設定します。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# --- 定数 ---
mu_E = 3.986004418e14 # 地球の重力パラメータ [m^3/s^2]
R_E = 6378.137e3 # 地球半径 [m]
g0 = 9.80665 # 標準重力加速度 [m/s^2]
rho0 = 1.225 # 海面の大気密度 [kg/m^3]
H_scale = 8500.0 # スケールハイト [m]
# --- ロケットパラメータ（Falcon 9風） ---
# 第1段
m_total_1 = 540000.0 # 第1段の総質量 [kg]（ペイロード+第2段含む）
m_prop_1 = 395000.0 # 第1段の推進剤質量 [kg]
T_1 = 7607000.0 # 第1段推力 [N]（海面、9基Merlin 1D）
Isp_1 = 282.0 # 第1段比推力 [s]（海面）
# 第2段
m_total_2 = 115000.0 # 第2段の総質量 [kg]（ペイロード含む）
m_prop_2 = 92000.0 # 第2段の推進剤質量 [kg]
T_2 = 934000.0 # 第2段推力 [N]（真空）
Isp_2 = 348.0 # 第2段比推力 [s]（真空）
# 空力パラメータ
C_D = 0.3 # 抗力係数
A_ref = 10.75 # 基準断面積 [m^2]（直径3.7mの円）
# 目標軌道
h_target = 200e3 # 目標高度 [m]
v_orbital = np.sqrt(mu_E / (R_E + h_target)) # 軌道速度
print(f"目標軌道速度: {v_orbital:.0f} m/s ({v_orbital/1e3:.2f} km/s)")
print(f"第1段TWR: {T_1 / (m_total_1 * g0):.2f}")
print(f"第1段燃焼時間: {m_prop_1 / (T_1/(Isp_1*g0)):.0f} s")
運動方程式の実装
def atmospheric_density(h):
"""指数大気モデル"""
if h < 0:
return rho0
return rho0 * np.exp(-h / H_scale)
def gravity_turn_eom(t, state, params):
"""
重力ターンの運動方程式
state = [v, gamma, r, theta, m]
v: 速度 [m/s]
gamma: 飛行経路角 [rad]（水平=0, 垂直=pi/2）
r: 地球中心からの距離 [m]
theta: ダウンレンジ角 [rad]
m: ロケット質量 [kg]
"""
v, gamma, r, theta, m = state
T, Isp, phase = params['T'], params['Isp'], params['phase']
# 重力加速度
g = mu_E / r**2
# 高度
h = r - R_E
# 大気抵抗
rho = atmospheric_density(h)
D = 0.5 * rho * v**2 * C_D * A_ref
# マックスQでのスロットリング（簡易モデル）
q_dyn = 0.5 * rho * v**2
throttle = 1.0
if phase == 1 and q_dyn > 30000:
throttle = 0.7 # 動圧が高いときスロットリング
T_actual = T * throttle
# 質量流量
m_dot = -T_actual / (Isp * g0)
# 速度が非常に小さい場合の処理
if v < 1.0:
dv = T_actual / m - D / m - g
dgamma = 0.0
else:
dv = T_actual / m - D / m - g * np.sin(gamma)
dgamma = -(g / v - v / r) * np.cos(gamma)
dr = v * np.sin(gamma)
dtheta = v * np.cos(gamma) / r
return [dv, dgamma, dr, dtheta, m_dot]
打ち上げシミュレーションの実行
def simulate_launch(kick_angle_deg):
"""
重力ターンの打ち上げシミュレーション
Parameters:
kick_angle_deg: ピッチオーバーキック角度 [度]
Returns:
results: 辞書（時間、速度、高度など）
"""
# 初期条件
v0 = 0.1 # 初期速度 [m/s]
gamma0 = np.pi / 2 # 垂直上昇
r0 = R_E
theta0 = 0.0
m0 = m_total_1
state = [v0, gamma0, r0, theta0, m0]
# シミュレーション結果の格納
t_all = [0.0]
v_all = [v0]
gamma_all = [gamma0]
r_all = [r0]
theta_all = [theta0]
m_all = [m0]
phase_all = [0]
# --- フェーズ0: 垂直上昇（10秒間） ---
params = {'T': T_1, 'Isp': Isp_1, 'phase': 0}
sol = solve_ivp(
lambda t, s: gravity_turn_eom(t, s, params),
(0, 10), state,
method='RK45', max_step=0.5,
rtol=1e-8, atol=1e-10
)
for k in range(1, len(sol.t)):
t_all.append(sol.t[k])
v_all.append(sol.y[0, k])
gamma_all.append(sol.y[1, k])
r_all.append(sol.y[2, k])
theta_all.append(sol.y[3, k])
m_all.append(sol.y[4, k])
phase_all.append(0)
state = [sol.y[i, -1] for i in range(5)]
# --- ピッチオーバーキック ---
kick_rad = np.radians(kick_angle_deg)
state[1] = np.pi / 2 - kick_rad # 経路角をキック
# --- フェーズ1: 第1段重力ターン ---
t_burn1 = m_prop_1 / (T_1 / (Isp_1 * g0))
t_remaining_1 = t_burn1 - 10 # 垂直上昇分を引く
params = {'T': T_1, 'Isp': Isp_1, 'phase': 1}
t_start = t_all[-1]
sol = solve_ivp(
lambda t, s: gravity_turn_eom(t, s, params),
(0, t_remaining_1), state,
method='RK45', max_step=0.5,
rtol=1e-8, atol=1e-10,
events=lambda t, s: s[4] - (m_total_1 - m_prop_1)
)
for k in range(1, len(sol.t)):
t_all.append(t_start + sol.t[k])
v_all.append(sol.y[0, k])
gamma_all.append(sol.y[1, k])
r_all.append(sol.y[2, k])
theta_all.append(sol.y[3, k])
m_all.append(sol.y[4, k])
phase_all.append(1)
state = [sol.y[i, -1] for i in range(5)]
# --- 段間分離（2秒のコースト） ---
state[4] = m_total_2 # 第2段の質量に切り替え
params_coast = {'T': 0, 'Isp': Isp_2, 'phase': 2}
t_start = t_all[-1]
sol = solve_ivp(
lambda t, s: gravity_turn_eom(t, s, params_coast),
(0, 2), state,
method='RK45', max_step=0.5,
rtol=1e-8, atol=1e-10
)
for k in range(1, len(sol.t)):
t_all.append(t_start + sol.t[k])
v_all.append(sol.y[0, k])
gamma_all.append(sol.y[1, k])
r_all.append(sol.y[2, k])
theta_all.append(sol.y[3, k])
m_all.append(sol.y[4, k])
phase_all.append(2)
state = [sol.y[i, -1] for i in range(5)]
# --- フェーズ3: 第2段燃焼 ---
t_burn2 = m_prop_2 / (T_2 / (Isp_2 * g0))
params = {'T': T_2, 'Isp': Isp_2, 'phase': 3}
t_start = t_all[-1]
# 目標高度・速度に達するか、燃料が尽きるまで
def event_orbit(t, s):
h = s[2] - R_E
v_orb = np.sqrt(mu_E / s[2])
return s[0] - v_orb # 軌道速度に達したら停止
event_orbit.terminal = True
event_orbit.direction = 1
sol = solve_ivp(
lambda t, s: gravity_turn_eom(t, s, params),
(0, t_burn2), state,
method='RK45', max_step=0.5,
rtol=1e-8, atol=1e-10,
events=event_orbit
)
for k in range(1, len(sol.t)):
t_all.append(t_start + sol.t[k])
v_all.append(sol.y[0, k])
gamma_all.append(sol.y[1, k])
r_all.append(sol.y[2, k])
theta_all.append(sol.y[3, k])
m_all.append(sol.y[4, k])
phase_all.append(3)
return {
't': np.array(t_all),
'v': np.array(v_all),
'gamma': np.array(gamma_all),
'r': np.array(r_all),
'theta': np.array(theta_all),
'm': np.array(m_all),
'phase': np.array(phase_all),
}
結果の可視化
# --- シミュレーション実行 ---
results = simulate_launch(kick_angle_deg=2.0)
t = results['t']
v = results['v']
gamma = results['gamma']
h = (results['r'] - R_E) / 1e3 # [km]
m = results['m']
phase = results['phase']
downrange = results['theta'] * R_E / 1e3 # [km]
# 動圧
rho_arr = np.array([atmospheric_density(hi * 1e3) for hi in h])
q_dyn = 0.5 * rho_arr * v**2
# --- 4パネルプロット ---
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 高度 vs 時間
ax = axes[0, 0]
for p, color, label in [(0, 'gray', 'Vertical'), (1, 'blue', 'Stage 1'),
(2, 'orange', 'Coast'), (3, 'red', 'Stage 2')]:
mask = phase == p
if np.any(mask):
ax.plot(t[mask], h[mask], color=color, linewidth=2, label=label)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Altitude [km]')
ax.set_title('高度プロファイル')
ax.legend()
ax.grid(True, alpha=0.3)
# 速度 vs 時間
ax = axes[0, 1]
for p, color, label in [(0, 'gray', 'Vertical'), (1, 'blue', 'Stage 1'),
(2, 'orange', 'Coast'), (3, 'red', 'Stage 2')]:
mask = phase == p
if np.any(mask):
ax.plot(t[mask], v[mask] / 1e3, color=color, linewidth=2, label=label)
ax.axhline(y=v_orbital / 1e3, color='green', linestyle='--',
alpha=0.7, label=f'Orbital velocity ({v_orbital/1e3:.1f} km/s)')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Velocity [km/s]')
ax.set_title('速度プロファイル')
ax.legend()
ax.grid(True, alpha=0.3)
# 飛行経路角 vs 時間
ax = axes[1, 0]
for p, color, label in [(0, 'gray', 'Vertical'), (1, 'blue', 'Stage 1'),
(2, 'orange', 'Coast'), (3, 'red', 'Stage 2')]:
mask = phase == p
if np.any(mask):
ax.plot(t[mask], np.degrees(gamma[mask]), color=color,
linewidth=2, label=label)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Flight Path Angle [deg]')
ax.set_title('飛行経路角')
ax.legend()
ax.grid(True, alpha=0.3)
# 動圧 vs 時間
ax = axes[1, 1]
ax.plot(t, q_dyn / 1e3, 'purple', linewidth=2)
max_q_idx = np.argmax(q_dyn)
ax.plot(t[max_q_idx], q_dyn[max_q_idx] / 1e3, 'r*', markersize=15,
label=f'Max Q = {q_dyn[max_q_idx]/1e3:.1f} kPa\n(t = {t[max_q_idx]:.0f} s)')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Dynamic Pressure [kPa]')
ax.set_title('動圧（Max Q）')
ax.legend()
ax.grid(True, alpha=0.3)
plt.suptitle('重力ターン打ち上げシミュレーション（ピッチキック = 2°）',
fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('gravity_turn_profile.png', dpi=150, bbox_inches='tight')
plt.show()
# --- 統計の出力 ---
print(f"\n=== シミュレーション結果 ===")
print(f"最終高度: {h[-1]:.1f} km")
print(f"最終速度: {v[-1]/1e3:.2f} km/s")
print(f"最終経路角: {np.degrees(gamma[-1]):.2f}°")
print(f"ダウンレンジ距離: {downrange[-1]:.0f} km")
print(f"総飛行時間: {t[-1]:.0f} s")
print(f"Max Q: {q_dyn.max()/1e3:.1f} kPa (t = {t[max_q_idx]:.0f} s)")
print(f"最終質量: {m[-1]:.0f} kg")
この4パネルプロットから、重力ターンの打ち上げの全体像が読み取れます。
-
高度プロファイル（左上）: 第1段燃焼中に急速に高度を稼ぎ、第2段でさらに高度を上げながら軌道投入します。コースト期間（段間分離）では慣性飛行で高度がわずかに変化します。
-
速度プロファイル（右上）: 第1段で約3 km/s、第2段でさらに約5 km/sを獲得し、合計で軌道速度に達します。第1段の速度獲得が遅いのは、大気抵抗と重力損失が大きいためです。
-
飛行経路角（左下）: 90°（垂直）から始まり、重力ターンにより徐々に0°（水平）に近づきます。第1段燃焼終了時点で30°程度、軌道投入時にはほぼ0°になります。
-
動圧（右下）: Max Qは打ち上げ後60〜80秒付近で発生し、約30 kPa程度に達します。実際のFalcon 9では約35 kPaで、スロットリングで抑制しています。
ピッチオーバーキック角度の最適化
ピッチキック角度を変えた場合の影響を比較します。
# --- ピッチキック角度の比較 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
kick_angles = [0.5, 1.0, 2.0, 3.0, 5.0]
colors = plt.cm.viridis(np.linspace(0.2, 0.9, len(kick_angles)))
for angle, color in zip(kick_angles, colors):
res = simulate_launch(angle)
axes[0].plot(res['t'], (res['r'] - R_E) / 1e3, color=color,
linewidth=1.5, label=f'Kick = {angle}°')
axes[1].plot(res['t'], res['v'] / 1e3, color=color,
linewidth=1.5, label=f'Kick = {angle}°')
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Altitude [km]')
axes[0].set_title('高度プロファイル（ピッチキック角度比較）')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Velocity [km/s]')
axes[1].set_title('速度プロファイル（ピッチキック角度比較）')
axes[1].axhline(y=v_orbital / 1e3, color='green', linestyle='--', alpha=0.5)
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('kick_angle_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、ピッチキック角度の影響が明確にわかります。
-
キック角度が小さすぎる（0.5°）場合: ロケットがほぼ垂直に上昇し続けるため、高度は稼げますが水平速度の獲得が遅くなります。重力損失が大きくなり、結果として軌道投入に必要な総ΔVが増加します。
-
キック角度が大きすぎる（5°）場合: 早い段階で経路が水平に近づくため、大気中での水平速度が大きくなります。しかし、低高度で水平速度を上げると空力損失が増大し、また十分な高度に達する前に燃料を使い切るリスクがあります。
-
最適なキック角度（1〜3°）: 重力損失と空力損失のバランスが取れた角度です。実際のロケットでも1〜2°程度のキック角度が使われることが多いです。
損失の内訳
# --- 損失の計算 ---
res = simulate_launch(2.0)
dt_arr = np.diff(res['t'])
v_final = res['v'][-1]
# 重力損失
g_arr = mu_E / res['r']**2
gravity_loss = np.sum(g_arr[:-1] * np.sin(res['gamma'][:-1]) * dt_arr)
# 空力損失
rho_sim = np.array([atmospheric_density((r - R_E)) for r in res['r']])
D_arr = 0.5 * rho_sim * res['v']**2 * C_D * A_ref
drag_loss = np.sum(D_arr[:-1] / res['m'][:-1] * dt_arr)
# 理想ΔV（ツィオルコフスキー）
dv_ideal_1 = Isp_1 * g0 * np.log(m_total_1 / (m_total_1 - m_prop_1))
dv_ideal_2 = Isp_2 * g0 * np.log(m_total_2 / (m_total_2 - m_prop_2))
dv_ideal = dv_ideal_1 + dv_ideal_2
print(f"\n=== ΔVバジェット ===")
print(f"理想ΔV（ツィオルコフスキー）: {dv_ideal/1e3:.2f} km/s")
print(f" 第1段: {dv_ideal_1/1e3:.2f} km/s")
print(f" 第2段: {dv_ideal_2/1e3:.2f} km/s")
print(f"最終速度: {v_final/1e3:.2f} km/s")
print(f"重力損失: {gravity_loss/1e3:.2f} km/s")
print(f"空力損失: {drag_loss/1e3:.2f} km/s")
print(f"合計損失: {(gravity_loss+drag_loss)/1e3:.2f} km/s")
print(f"損失率: {(gravity_loss+drag_loss)/dv_ideal*100:.1f}%")
この結果から、打ち上げのΔVバジェットの内訳が定量的にわかります。理想的なΔV（ツィオルコフスキーの式から得られる値）のうち、約15〜20%が重力損失と空力損失で失われています。重力損失が空力損失を大きく上回ることも確認できます。
まとめ
本記事では、ロケットの重力ターン軌道の理論とシミュレーションを解説しました。
- 軌道投入の本質: 水平方向に約7.8 km/sの速度を獲得すること。大気と重力による損失を加えた実際の必要ΔVは約9〜9.5 km/s
- 重力ターンの仕組み: 推力を速度方向に維持し、重力の法線成分で飛行経路を自然に水平方向に曲げる。ロケットに横方向の荷重がかからない効率的な方式
- ピッチオーバーキック: 垂直上昇後のわずかな傾き（1〜3°）が飛行プロファイル全体を決定する重要パラメータ
- 重力損失: 重力に逆らう上昇で失うΔV。TWRが大きいほど小さい（約1.0〜1.5 km/s）
- 空力損失: 大気抵抗で失うΔV。低空での高速飛行が多いほど大きい（約0.1〜0.5 km/s）
- Max Q: 動圧最大の瞬間（打ち上げ後60〜80秒、高度10〜15 km）。構造設計のクリティカルポイント
次のステップとして、以下の記事も参考にしてください。