10. 単段ロケット上昇

例題07の Goddard 問題を、メートル・キログラム・秒の実際の単位で解きます。空気の密度は高度とともに指数関数的に薄くなり、重力は地球中心からの距離の2乗に反比例します。桁の大きく違う量を同時に扱うため、ソルバーが解きやすいように単位をそろえる工夫も学びます。

項目	内容
状態	地球中心からの距離 $R$ 、鉛直速度 $v$ 、質量 $m$
制御	推力 $T$
目的	最終的な距離 $R(t_f)$ （到達高度）を最大にする
区間	1区間、50ノード、終了時刻は自由
新しく学ぶこと	実際の単位、指数関数の大気モデル、逆2乗の重力、単位のスケーリング

状況

全体の質量 5000 kg、そのうち 3000 kg が推進剤のロケットを、地表から真上へ打ち上げます。推進剤を使い切ったあとも惰性で上がり、上昇が止まった点が到達高度です。この到達高度を最大にする推力の出し方を求めます。

考え方は例題07と同じです。全開で吹かすと低空で速度が出すぎて空気抵抗で損をし、ゆっくり吹かすと重力に逆らう時間が長くなって損をします。違いは、物理量を実際の単位で書くことです。位置は地表からの高度ではなく地球中心からの距離 $R$ で表します（地表では $R = R_e = 6371$ km）。機体は真上にしか動かない1次元の運動です。

問題設定

運動方程式

\begin{aligned} \dot{R} &= v \\[2pt] \dot{v} &= \frac{T - D}{m} - g(R) \\[2pt] \dot{m} &= -\frac{T}{g_0 I_{sp}} \end{aligned}

$R$ は地球中心からの距離 [m]、 $v$ は鉛直上向きの速度 [m/s]、 $m$ は機体の質量 [kg]、 $T$ は推力 [N] です。
質量の式の $g_0 I_{sp}$ は排気速度 [m/s] です。比推力 $I_{sp}$ は「1 kg 重の推進剤で何秒間 1 kg 重の推力を出せるか」を表すエンジンの燃費の指標で、大きいほど同じ推力で推進剤の消費が少なくなります。

右辺に現れる空気抵抗 $D$ 、空気の密度 $\rho$ 、重力加速度 $g$ は次のとおりです。

D = \tfrac{1}{2}\, \rho(R - R_e)\, v^2\, C_d\, A

\rho(h) = \rho_0 \exp\!\left(-\frac{h}{H}\right), \qquad g(R) = \frac{GM}{R^2}

空気抵抗は動圧 $\tfrac12 \rho v^2$ に抵抗係数 $C_d$ と断面積 $A$ を掛けたものです。
密度 $\rho$ は高度 $h = R - R_e$ の指数関数で、高度が $H = 8500$ m 上がるごとに $1/e$ （約 37 %）に減ります。 $H$ をスケールハイトと呼びます。
重力加速度は地球中心からの距離の2乗に反比例します。 $GM$ は地球の重力定数です。地表では約 9.82 m/s² になります。

境界条件

\begin{aligned} &R(0) = R_e,\quad v(0) = 0,\quad m(0) = 5000\ \mathrm{kg} \\[2pt] &v(t_f) = 0,\quad m(t_f) = 2000\ \mathrm{kg} \end{aligned}

地表で静止した状態から出発します。終端では推進剤 3000 kg を使い切り（質量 2000 kg）、上昇が止まった（速度 0）瞬間とします。終端の質量 2000 kg は、初期質量に比 $0.4$ を掛けた値です。

制約

\begin{aligned} &0 \le T \le T_{\max}, \qquad R \ge R_e, \qquad v \ge 0 \\[2pt] &2000 \le m \le 5000\ \mathrm{kg}, \qquad t_f \ge 10\ \mathrm{s} \end{aligned}

最大推力は初期重量の2倍です。

T_{\max} = 2\, m(0)\, g_0 = 98066.5\ \mathrm{N}

目的関数

J = -\frac{R(t_f)}{R_e}

到達点の距離を最大にするため、マイナスをつけて最小化します。 $R_e$ で割っているのは値を 1 程度にするためで、答えには影響しません。

パラメータ

記号	意味	値
$GM$	地球の重力定数	$3.986004418 \times 10^{14}$ m³/s²
$R_e$	地球の半径	6371 km
$g_0$	標準重力加速度	9.80665 m/s²
$m(0)$	初期質量	5000 kg
$m(t_f)$	推進剤を使い切ったときの質量	2000 kg（初期質量の 0.4 倍）
$C_d$	抵抗係数	0.2
$A$	断面積	10 m²
$I_{sp}$	比推力	300 s
$T_{\max}$	最大推力	初期重量の 2 倍（98066.5 N）
$\rho_0$	地表の空気の密度	1.225 kg/m³
$H$	スケールハイト	8500 m

最大推力で吹かすと、推進剤は毎秒 $T_{\max}/(g_0 I_{sp}) \approx 33.3$ kg 減ります。全開なら 3000 kg を 90 s で使い切る計算です。

単位のスケーリング

この問題では、距離は $6.4 \times 10^6$ m、質量は $10^3$ kg、時間は $10^2$ s、推力は $10^5$ N と、量の桁が大きく違います。そのまま最適化すると、ソルバーは変数ごとの動かし方の感覚がつかめず、収束しにくくなります。そこで各量を次の代表的な大きさで割り、1 程度の値にしてから最適化します（単位をそろえる）。

量	単位にとった値	大きさ
距離	$R_e$	6371 km
速度	$\sqrt{GM/R_e}$ （地表での円軌道速度）	約 7.91 km/s
質量	$m(0)$	5000 kg
時間	距離の単位 ÷ 速度の単位	約 805 s
推力	質量 × 距離 ÷ 時間²	約 49.1 kN

運動方程式や条件は物理単位のまま書き、割り算はライブラリが内部で行います。

初期推定

時間区間を $[0, 600]$ s とし、次のように与えます。

距離は $R_e$ から $R_e + 50$ km へなめらかに変わる3次式（両端で傾き 0）。
速度は 0。
質量は 5000 kg から 2000 kg へ変わる3次式（出発時の傾き −0.6 kg/s、終端で傾き 0）。
推力は最大推力から 0 へ下がる3次式（両端で傾き 0）。

結果

高度 R − Re・速度・質量・推力の時間履歴。点は50個のノードでの値、線はそれらを通る多項式、赤い点線は最大推力 98 kN です。

上昇の様子。左は高度、右は速度と推力の時間変化で、縦の黒線が現在の時刻です。

到達高度は 65.5 km です。
推進剤は約 114 s で尽き、その後は惰性で上がって 188 s に最高点に達します。
推力は「最大推力 → 中間の推力 → ゼロ」と変わります。例題07と同じ形です。

解説

例題07と同じ形の答え

最初の約 20 s は最大推力で吹かして速度を稼ぎます。速度が上がると空気抵抗が大きくなるので、推力を 70〜80 kN 程度に落とし、抵抗による損失と重力による損失のつり合う速度を保ちながら上昇します（特異弧）。高度が上がって空気が薄くなると、残りの推進剤を燃やし切って速度を最大にし、あとは惰性で上がります。

最大推力のまま吹かせば 90 s で燃え尽きるところを、推力を絞ることで燃焼が約 114 s まで延びています。空気の濃い低空では急がず、空気が薄くなってから速度を上げる、というのが最適な使い方です。

例題07は量を無次元にした問題でしたが、物理の構造は同じなので、実際の単位で書いても同じ形の答えになります。

推力の波打ち

中間の推力の区間で推力が波打ち、燃焼終了の前後でなめらかにゼロへ落ちています。本当の最適解では推力はゼロへ不連続に切り替わります。例題07と同じく、不連続な関数を1本の多項式で表そうとしたためです（区間に分ける）。例題08のように燃焼終了で区間を分ければ、この振動は小さくできます。

試してみよう

抵抗係数 self.Cd を 0（空気抵抗なし）にする。特異弧がなくなり、最大推力で吹かして燃やし切る答えに近づくか、到達高度がどれだけ上がるかを見る。
比推力 self.Isp を 300 から 450 に上げる。同じ推進剤で長く吹かせるので、到達高度が大きく伸びることを確かめる。
最大推力 self.max_thrust（初期重量に対する比）を 1.5 に下げると、最大推力で吹かす時間が長くなり、特異弧が短くなるかを見る。

実行する

uv run python examples/10_rocket_ascent_single_stage.py

examples/10_rocket_ascent_single_stage.py を GitHub で見る