14. 大樹からの2段軌道投入

例題13 の2段式ロケットの問題を、実際の小型ロケットに近い条件に置き換えます。 大気は標準大気の表、抵抗係数はマッハ数ごとの表で与え、地球の自転による初速と、外の気圧で変わるエンジンの性能も入れます。 問題の形は例題13と同じでも、モデルを詳しくすると、何が変わり何が変わらないでしょうか。

項目	内容
状態	地球中心からの距離 $R$、角度 $\theta$、鉛直速度 $V_r$、水平速度 $V_t$、質量 $m$
制御	推力の鉛直成分 $T_r$ と水平成分 $T_t$
目的	軌道投入時の2段目の質量 $m(t_f)$ を最大にする
区間	2区間（20ノードと30ノード）、分離時刻も終了時刻も自由
新しく学ぶこと	表で与えた大気と抵抗係数、地球の自転による初速、外の気圧による比推力と推力の変化

状況

北海道大樹町（北緯 42.5°）から、直径 1.8 m の2段式ロケットで 100 kg の荷物（ペイロード）を、高度 561 km、傾斜角 96.7° の円軌道へ運びます。 傾斜角が 90° を超える軌道は、地球の自転と逆向きに少し西へ傾いて南北に回る軌道で、地球観測衛星によく使われる太陽同期軌道がこの種類です。

例題13からの主な違いは次の4つです。

• 大気：空気密度・気圧・音速を、高度 0〜86 km の米国標準大気の表から線形補間で求めます。86 km より上では密度と気圧を 0 とします。
• 抵抗係数：マッハ数（速さ ÷ 音速）ごとの表から求めます。音速付近で大きくなる、実際の機体に近い形です。
• 地球の自転：射点は地球とともに東へ動いているので、静止していても宇宙から見ると速度を持っています。
• エンジンの性能：ノズルの出口に外の気圧がかかるので、空気の薄い上空ほど推力と比推力が上がります。

問題設定

座標系

例題12 と同じです。$R$ は地球中心からの距離（高度 $h = R - R_e$）、$\theta$ は射点から地球中心のまわりに回った角度、 $V_r$ は真上向きの速度、$V_t$ は飛ぶ向きの水平速度、$T_r$ と $T_t$ は推力の鉛直成分と水平成分です。速度は地球の自転とともに回らない座標で測ります。

区間と接続条件

区間1（1段目）を $0 \le t \le t_s$、区間2（2段目）を $t_s \le t \le t_f$ とし、$t_s$ と $t_f$ を最適化で決めます。 分離では位置と速度が連続で、質量は2段目の初めの質量 $m_2$ から再出発します。

$$\begin{aligned} R(t_s^+) &= R(t_s^-) \\ \theta(t_s^+) &= \theta(t_s^-) \\ V_r(t_s^+) &= V_r(t_s^-) \\ V_t(t_s^+) &= V_t(t_s^-) \\ m(t_s^+) &= m_2 \end{aligned}$$

運動方程式

形は例題13と同じです。段 $i$（$i = 1, 2$）の値を使います。

$$\begin{aligned} \dot{R} &= V_r \\[2pt] \dot{\theta} &= \frac{V_t}{R} \\[2pt] \dot{V}_r &= \frac{T_r - D_r}{m} - g(R) + \frac{V_t^2}{R} \\[2pt] \dot{V}_t &= \frac{T_t - D_t}{m} - \frac{V_r V_t}{R} \\[2pt] \dot{m} &= -\frac{\sqrt{T_r^2 + T_t^2}}{g_0\, I_{sp,i}(h)} \end{aligned}$$

$V_t^2/R$ は遠心力、$-V_r V_t/R$ は上昇すると水平速度が減る効果（角運動量の保存）です。重力は $g(R) = \mu/R^2$ で、終端の円軌道速度と同じ重力定数 $\mu$ を使います。

大気と抵抗

空気密度 $\rho(h)$、気圧 $p(h)$、音速 $a(h)$ は標準大気の表を線形補間した関数です。マッハ数と抵抗は次のとおりです。

$$\begin{aligned} V &= \sqrt{V_r^2 + V_t^2} \\[2pt] M &= V / a(h) \\[2pt] D_r &= \tfrac12 \rho(h)\, V_r V\, C_D(M) A \\[2pt] D_t &= \tfrac12 \rho(h)\, V_t V\, C_D(M) A \end{aligned}$$

抵抗係数 $C_D(M)$ は表を線形補間したもので、亜音速で約 0.28、マッハ 1.2〜1.4 付近で最大の約 0.68、 マッハ 4.5 付近で約 0.34 になります。表の範囲外は直線で延長します。断面積は直径 1.8 m から $A \approx 2.54$ m² です。

比推力と推力の変化

ロケットエンジンの推力は、ノズル出口の圧力と外の気圧の差の分だけ変わります。この例題では、表で与えた比推力 （1段目は地上での値 261 s、2段目は真空での値 322 s）からのずれを次のように表します。

$$\begin{aligned} \Delta T_1(h) &= \{p(h_0) - p(h)\}\, A_{e,1} \\[2pt] \Delta T_2(h) &= -p(h)\, A_{e,2} \end{aligned}$$

$A_{e,i}$ はノズル出口の面積、$h_0 = 10$ m は射点の高度です。1段目では、上昇して外の気圧が下がるほど推力が増えます。 2段目は真空での値が基準なので、外の気圧 $p(h)$ がノズル出口を押し返す分だけ推力が減ります。コードでは Rocket.delta_thrust がこの $\Delta T_i(h)$ を返します。 比推力はこの推力の増分を基準の流量で割った分だけ変わります。

$$\begin{aligned} \dot{m}_{\mathrm{ref},i} &= \frac{T_{\max,i}}{g_0\, I_{sp,i}} \\[2pt] I_{sp,i}(h) &= I_{sp,i} + \frac{\Delta T_i(h)}{g_0\, \dot{m}_{\mathrm{ref},i}} \end{aligned}$$

推力の上限も同じだけ変えます（制約）。最大推力 $T_{\max,i} + \Delta T_i$ で吹かしたときの流量は、 $\{T_{\max,i} + \Delta T_i\}/\{g_0 I_{sp,i}(h)\} = \dot{m}_{\mathrm{ref},i}$ となり高度によらず一定です。 「エンジンが送り出す推進剤の流量は同じで、外の気圧が低いほど推力が大きい」という実際のエンジンの性質を表しています。 1段目の場合、真空では推力が約 41.6 kN 増え、比推力は約 294 s になります。 2段目を使う高度では気圧がほぼ 0 なので、2段目の補正は結果にほとんど影響しません。

地球の自転と初速

大樹町の地表は、地球の自転で東向きに

$$v_{\mathrm{site}} = \omega R_e \cos\phi = 464.6\ \mathrm{m/s} \times \cos 42.5^\circ \approx 342.5\ \mathrm{m/s}$$

で動いています。$\omega = 7.292 \times 10^{-5}$ rad/s は自転の角速度、$\phi$ は射点の緯度、464.6 m/s は赤道での自転速度です。 傾斜角 $i$ の軌道面が射点を通る向き（北から測った方位角 $A_z$）は、球面三角法から $\sin A_z = \cos i / \cos\phi$ で決まります。 自転の速度のうち、この向きの成分が水平速度の初期値です。

$$V_0 = v_{\mathrm{site}} \sin A_z = \omega R_e \cos i \approx -54.2\ \mathrm{m/s}$$

緯度が打ち消し合い、結果は傾斜角だけで決まります。傾斜角が 90° を超える軌道へは西寄り（$A_z \approx -9^\circ$、北北西）に打ち上げるので、 自転の速度は助けにならず、わずかに逆向きの初速から出発することになります。軌道面と直角な成分（約 338 m/s）は、2次元のモデルでは無視しています。 抵抗もこの速度から計算しており、大気が地球と一緒に回っていることは考えていません。

境界条件

$$\begin{aligned} &R(0) = R_e + 10\ \mathrm{m},\quad \theta(0) = 0 \\ &V_r(0) = 0,\quad V_t(0) = V_0 \\ &m(0) = m_1 + m_2 \\[4pt] &R(t_f) = R_e + 561\ \mathrm{km} \\ &V_r(t_f) = 0 \\ &V_t(t_f) = \sqrt{\mu / R(t_f)} \end{aligned}$$

$m_1 = 13000$ kg は1段目、$m_2 = 1566.7$ kg は2段目（ペイロード込み）の質量で、打上げ時の質量は 14566.7 kg です。 終端の水平速度は高度 561 km の円軌道速度（約 7583.0 m/s）です。

制約

$$\begin{aligned} &R \ge R_e \\ &m_{d,1} + m_2 \le m \le m_1 + m_2 \quad (\text{区間1}) \\ &m_{d,2} + m_{\mathrm{pl}} \le m \le m_2 \quad (\text{区間2}) \\ &|\boldsymbol{T}| = \sqrt{T_r^2 + T_t^2} \le T_{\max,i} + \Delta T_i(h) \\ &q = \tfrac12 \rho V^2 \le 500\ \mathrm{kPa} \\ &\left|\frac{\boldsymbol{T} - \boldsymbol{D}}{m}\right| \le 20\, g_0 \end{aligned}$$

• 区間1の質量の下限 $m_{d,1} + m_2 = 2866.7$ kg は、1段目が空になり2段目が載っている状態です。
• 区間2の質量の下限は、2段目の空の質量 $m_{d,2} = 220$ kg とペイロード $m_{\mathrm{pl}} = 100$ kg の合計 320 kg です。荷物を軌道まで運べない解は認めません。
• 推力は大きさだけを制限し、向きは自由です。上限は高度によって変わります（比推力と推力の変化）。
• $q$ は動圧で、機体が空気から受ける力の目安です。上限は 500 kPa と大きく、実際にはほとんど効きません。
• 推力と抵抗による加速度（重力を含めない値）は 20 G 以下です。

目的関数

$$J = -\frac{m(t_f)}{13000\ \mathrm{kg}}$$

軌道に入ったときの2段目の質量を最大にします。1段目の質量で割って大きさをそろえています。

パラメータ

記号	意味	1段目	2段目
$m_{d,i}$	空の質量	1300 kg	220 kg
—	構造の比率（空の質量 ÷ 段の質量）	10 %	15 %
—	段の質量（推進剤込み）	13000 kg	1466.7 kg
—	ペイロード	—	100 kg
$T_{\max,i}$	基準の最大推力	33.3 tf（約 326.6 kN）	4.2 tf（約 41.2 kN）
$I_{sp,i}$	基準の比推力	261 s（地上）	322 s（真空）
—	エンジン	9基、ノズル面積比 20	1基、ノズル面積比 140
$A_{e,i}$	ノズル出口の面積（スロート径 53.9 mm から）	約 0.411 m²（9基分）	約 0.319 m²

地球半径 $R_e = 6371$ km、$g_0 = 9.80665$ m/s²、重力と目標速度の計算に使う重力定数 $\mu = 3.986 \times 10^{14}$ m³/s²、 射点の緯度 42.506°、地球の自転角速度 $7.2921159 \times 10^{-5}$ rad/s、傾斜角 96.7°、目標高度 561 km です。 単位のそろえ方は例題12・13と同じで、質量の単位は 13000 kg です。

初期推定

区間1を 0〜200 秒、区間2を 200〜800 秒として与えます。

• 距離 $R$ は地表から目標の高さへ、角度 $\theta$ は 0 から 25° へなめらかに増える3次式。
• 鉛直速度は 0、水平速度は $V_0$ から目標速度へ直線的に増やします。
• 質量は、区間1で打上げ質量から下限 2866.7 kg へ、区間2で $m_2$ から下限 320 kg へ減る3次式。
• 推力は各段で、鉛直成分を $0.9\,T_{\max,i}$、水平成分を $0.1\,T_{\max,i}$ から 0 へ下げる3次式。

結果

飛行経路。横軸は射点からの地表距離、縦軸は高度で、高度方向を拡大しています。青は1段目、橙は2段目の区間です。点はノードでの値、線はそれらを通る多項式、橙色の矢印はエンジンが動いているノードでの推力の向き（間引いて表示）です。

時間履歴。高度、水平速度（実線）と鉛直速度（破線）、質量、推力の大きさ、推力の仰角（エンジン停止中は非表示）、動圧。青は1段目、橙は2段目、縦の破線は段分離の時刻です。赤い点線は各段の基準の最大推力で、1段目の上限は上昇とともにこの線より大きくなります。

上昇の様子。右側は速度・推力の大きさ・推力の仰角で、縦の実線が現在の時刻です。エンジンが止まっている間は早送りしています。

• 段分離は時刻 139.0 秒です。分離直前の質量は 2,866.7 kg で下限に達し、 1段目の推進剤を使い切っています。分離時の水平速度は 3,356 m/s です。
• 2段目は分離直後から時刻 293 秒ごろまで燃焼し、長い慣性飛行の後、 2,302 秒ごろから 2,311 秒まで短く燃焼して円軌道に入ります。
• 軌道に残る質量は 342.5 kg です。区間2の質量の下限 320 kg（空の質量 220 kg + ペイロード 100 kg）をわずかに上回り、 この下限の制約は効いていません。100 kg の荷物を運べることがわかります。余りは使い残した推進剤です。
• 終端の高度は 561.0 km、速さは 7,583.0 m/s で目標どおりです。
• 動圧の最大はノードで 45.9 kPa、加速度の最大は 6.99 G で、 どちらも上限（500 kPa、20 G）よりずっと小さく、この解では効いていません。

解説

詳しいモデルでも形は同じ

飛び方は例題13とよく似ています。ほぼ真上に上がってから水平へ倒れ、2段目で軌道に近い速さまで加速し、 地球のまわりを大きく回る慣性飛行で目標高度まで上がってから、最高点で円軌道にします（ホーマン遷移の形）。 大気や抵抗係数を表にしても、最適な飛び方の骨格は変わりません。簡単なモデルで性質をつかみ、詳しいモデルで数値を詰めるという使い分けができます。

詳しいモデルで変わるもの

• 抵抗：抵抗係数は音速付近で約 2.4 倍に増えます。動圧のグラフは打上げ後しばらくで最大になり、その後は空気が薄くなって急に下がります。
• 自転：西寄りの打上げでは自転が約 54 m/s 足を引っ張ります。東向きに打ち上げる軌道なら逆に数百 m/s 得をします。
• エンジン：1段目の推力と比推力は上昇とともに増えます。推力のグラフで1段目のノードが基準の最大推力（赤い点線）を上回っているのは、 上空ほど上限が大きくなるためです。2段目は空気の薄い高度で使うので、この補正はほとんど効きません。

1段目の推力のジグザグ

ソルバーは収束を報告していますが、1段目の推力はノードごとに大きく上下しています。最終質量は1段目の推進剤の総量でほぼ決まり、 1段目の中で推力をいつ強めいつ弱めるかには鈍感なため、推力の形がはっきり決まらないと考えられます。 1段目の推力の細かい形は、意味のある特徴として読まないでください。ノード数を増やしたり、推力の変化の速さに制約を加えたりすると改善することがあります。

試してみよう

• Rocket の MaxQ を 500000 Pa から 30000 Pa（30 kPa）に下げると、動圧の制約が効くようになります。動圧のグラフの頭が抑えられる様子と、最終質量の減り方を見る。
• inclination を 96.7 から 45 に変えると、初速 $V_0$ が東向き（正）になります。最終質量がどれだけ増えるかを確かめる。
• Mpayload を 100 から 150 kg に増やし、最終質量から空の質量とペイロードを引いた「余り」がまだ正かどうかを調べる。

実行する

uv run python examples/14_polar_tsto_taiki.py

examples/14_polar_tsto_taiki.py を GitHub で見る