11. 2段ロケット上昇

推進剤を使い終えたタンクやエンジンは、もう役に立たない重りです。途中で切り離せば、残りの推進剤で軽い機体を加速できます。 これが多段ロケットの考え方です。例題10と同じ推進剤 3000 kg で、1段目の構造を切り離すと、到達高度はどこまで伸びるでしょうか。

項目	内容
状態	地球中心からの距離 $R$、鉛直速度 $v$、質量 $m$
制御	推力 $T$
目的	最終的な距離 $R(t_f)$（到達高度）を最大にする
区間	2区間（1段目・2段目、各25ノード）、分離時刻と終了時刻は自由
新しく学ぶこと	段分離で質量が不連続に変わること、接続条件を自分で書くこと

状況

全体の質量 5000 kg のロケットを地表から真上へ打ち上げます。物理のモデル（大気、重力、エンジンの性能、最大推力）は 例題10とまったく同じです。違うのは、途中で段分離をすることです。

• 1段目（全質量 3000 kg、うち推進剤 1800 kg）で上昇し、ある時刻 $t_1$ に1段目の構造 1200 kg を切り離します。
• その後は2段目（全質量 2000 kg、うち推進剤 1200 kg）で上昇を続け、推進剤を使い切って最終質量 800 kg になります。
• 上昇が止まった点が到達高度です。

各段の空の質量は、その段の全質量に比 0.4 を掛けた値です（Mc = 0.4）。使える推進剤は合計 $1800 + 1200 = 3000$ kg で、例題10と同じです。 違いは、最後まで運ぶ「推進剤以外の質量」が例題10の 2000 kg から 800 kg に減ることです。 各段が使える推進剤はその段のタンクに入っている分だけで、1段目の推進剤を2段目で使うことはできません。分離の時刻 $t_1$ は最適化で決めます。

分離の瞬間に質量は 1200 kg 不連続に減ります。1本の多項式ではこの段差を表せないので、時間を1段目と2段目の2つの区間に分けます。

問題設定

運動方程式

両方の区間で、例題10と同じ式を使います。

$$\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}$$ $$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}$$

$R$ は地球中心からの距離、$v$ は鉛直上向きの速度、$m$ は質量、$T$ は推力、$D$ は空気抵抗、$\rho$ は高度 $h = R - R_e$ での空気の密度、 $g$ は重力加速度です。各項の意味は例題10を参照してください。

区間と接続条件

区間1（1段目）を $0 \le t \le t_1$、区間2（2段目）を $t_1 \le t \le t_f$ とします。分離時刻 $t_1$ で次の条件を課します。

$$\begin{aligned} R^{(1)}(t_1) &= R^{(2)}(t_1) \\[2pt] v^{(1)}(t_1) &= v^{(2)}(t_1) \\[2pt] m^{(2)}(t_1) &= 2000\ \mathrm{kg} \\[2pt] m^{(1)}(t_1) &= m^{(2)}(t_1) + 1200\ \mathrm{kg} \end{aligned}$$

上付きの $(1), (2)$ は区間の番号です。位置と速度は分離の前後でつながりますが、質量は 1200 kg だけ跳びます。 分離直後の質量は2段目の全質量 2000 kg に固定し、分離直前の質量は 3200 kg（1段目の空の質量 1200 kg + 2段目）になります。 これで1段目が燃やす推進剤は $5000 - 3200 = 1800$ kg、2段目は $2000 - 800 = 1200$ kg と、どちらも各段の推進剤の量に決まります。

例題08・09では「すべての状態が境目でつながる」という条件をライブラリが自動で課していました。 この例題では質量が跳ぶので自動の接続を外し、上の4つの式を自分で書いています （区間をつなぐ）。

境界条件

$$\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) = 800\ \mathrm{kg} \end{aligned}$$

地表で静止した状態から出発し、推進剤を使い切って（質量 800 kg）上昇が止まった瞬間を終端とします。 800 kg は2段目の全質量 2000 kg に比 0.4 を掛けた値です。

制約

両方の区間で、次の範囲を守ります。

$$\begin{aligned} &0 \le T \le T_{\max}, \qquad R \ge R_e, \qquad v \ge 0 \\[2pt] &3200 \le m \le 5000\ \mathrm{kg} \quad (\text{区間1}) \\[2pt] &800 \le m \le 2000\ \mathrm{kg} \quad (\text{区間2}) \end{aligned}$$

質量の範囲は、各段が自分の推進剤の量を超えて燃やさないことを表します（接続条件で端の値は決まっているので、途中のノードのための範囲です）。 区間の終わりの時刻には、区間が潰れないように小さな下限をつけます。

$$t_1 \ge 10\ \mathrm{s}, \qquad t_f \ge 10\ \mathrm{s}$$

最大推力 $T_{\max} = 2\, m(0)\, g_0 = 98066.5$ N は1段目・2段目で共通です。

目的関数

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

例題10と同じく、到達点の距離を最大にします。

パラメータ

記号	意味	値
$m(0)$	初期質量（1段目 3000 kg + 2段目 2000 kg）	5000 kg
—	各段の空の質量 ÷ 全質量 Mc	0.4
$\Delta m$	分離で切り離す質量（1段目の空の質量）	1200 kg
$m(t_f)$	2段目が推進剤を使い切ったときの質量（2段目の空の質量）	800 kg
$C_d$	抵抗係数	0.2
$A$	断面積	10 m²
$I_{sp}$	比推力	300 s
$T_{\max}$	最大推力	初期重量の 2 倍（98066.5 N）
$\rho_0,\ H$	地表の空気の密度、スケールハイト	1.225 kg/m³、8500 m
$GM,\ R_e,\ g_0$	地球の重力定数・半径・標準重力加速度	例題10と同じ

単位のスケーリングも例題10と同じです（距離 $R_e$、速度 $\sqrt{GM/R_e}$、質量 5000 kg など）。

初期推定

区間の時刻は $[0, 300]$ s と $[300, 600]$ s から始めます。

• 距離は全体の時間にわたって $R_e$ から $R_e + 50$ km へなめらかに変わる3次式、速度は 0。
• 質量は区間ごとの直線で、区間1は 5000 kg から 3200 kg へ、区間2は 2000 kg から 800 kg へ減らします。
• 推力は全体の時間にわたって最大推力から 0 へ下がる3次式。

結果

高度 R − Re・速度・質量・推力の時間履歴。青が1段目（区間1）、橙が2段目（区間2）で、点は各区間25個のノードでの値、線は区間ごとの多項式です。縦の破線は分離時刻、赤い点線は最大推力 98 kN です。

上昇の様子。左の高度は1段目を青、2段目を橙で描いています。右は速度と推力の時間変化で、縦の破線が分離時刻、縦の黒線が現在の時刻です。

• 到達高度は 184.1 km で、同じ推進剤量の単段（例題10）の 65.5 km の約2.8倍です。
• 1段目は約 64 s で推進剤を使い切り、少し惰性で上がってから時刻 69.2 s に分離します。分離時刻は下限の 10 s から離れており、制約で決まった値ではありません。
• 質量は分離の直前に 3,200 kg、直後に 2,000 kg で、ちょうど 1200 kg 減っています。
• 2段目は分離直後に最大推力で点火し、推進剤は約 152 s で尽き、 307 s に最高点に達します。

解説

重りを捨てると伸びる理由

ロケットが得られる速度の増分は、ツィオルコフスキーの式

$$\Delta v = g_0 I_{sp} \ln\frac{m_{\text{start}}}{m_{\text{end}}}$$

で決まります（重力と空気抵抗を無視した場合）。$m_{\text{start}}/m_{\text{end}}$ は燃焼前後の質量の比です。 単段では 5000 kg が 2000 kg になるだけなので、比は 2.5 です。2段式では、分離で重りを捨てることで、 2段目の燃焼の間は機体がずっと軽く、同じ推進剤からより大きな速度が得られます。 速度の最大値は例題10の約 780 m/s に対して約 1420 m/s になり、到達高度は大きく伸びます。

1段目の燃焼と分離

1段目は最初の約 22 s だけ最大推力で、その後は推力を絞りながら約 64 s で推進剤 1800 kg を使い切ります。 例題10と同じく、空気の濃い低空で速く飛びすぎると抵抗で損をするので、推力を抑えて速度を調整しています。

推進剤が尽きたあと、すぐには分離せず、数秒だけ惰性で上がってから分離しています。この間、機体は1段目の空の構造を抱えたままですが、 推進剤を使っていないので損はありません。分離の時刻は、実質的には「2段目に点火する時刻」の選択になっています。

段ごとの推進剤の量を決めないと（試してみようの1つ目）、ソルバーは推進剤をなるべく2段目に回して早く分離しようとします。 分離の前に燃やした推進剤は、捨てる予定の 1200 kg まで一緒に加速することになるからです。 しかしそれは「2段目のタンクに入りきらない推進剤を2段目で燃やす」ことになり、実際の機体では実現できません。

2段目の推力

2段目では、例題10と同じく「最大推力 → 中間の推力 → ゼロ」の形が現れます。分離時刻ですでに一定の速度があり、 機体も軽いので、最大推力の時間は短く、すぐに中間の推力（特異弧）に移ります。 ここでも推力の波打ちや燃焼終了付近のなめらかな落ち方は、1本の多項式で不連続を表そうとしたためです。

試してみよう

• equality から2段目の初期質量の条件 result.equal(m2[0], obj.M0_2nd) を外し、区間ごとの質量の範囲を set_states_bounds(2, 0, obj.Mdry_1st + obj.Mdry_2nd, obj.M0)、 set_states_bounds(2, 1, obj.Mdry_2nd, obj.M0 - obj.Mdry_1st) に広げて、段ごとの推進剤の量を自由にする。 分離は下限の 10 s に張り付き、到達高度は約 261 km まで伸びますが、2段目の初期質量は約 3530 kg と、2段目の全質量 2000 kg を超えてしまいます。
• Rocket の Mc を 0.4 から 0.3 にして構造を軽くし、到達高度がどれだけ伸びるかを見る（約 476 km になります）。

実行する

uv run python examples/11_rocket_ascent_two_stage.py

examples/11_rocket_ascent_two_stage.py を GitHub で見る