12. 極座標・単段軌道投入

これまでのロケットは真上に上がるだけでした。人工衛星になるには、高く上がるだけでなく、 地球を回り続けられるだけの横向きの速度が必要です。 上へ向けるか横へ向けるか、推力の向きをどう変えていけば、最も多くの質量を軌道に残せるでしょうか。

項目	内容
状態	地球中心からの距離 $R$、角度 $\theta$、鉛直速度 $V_r$、水平速度 $V_t$、質量 $m$
制御	推力の鉛直成分 $T_r$ と水平成分 $T_t$
目的	軌道投入時の質量 $m(t_f)$ を最大にする
区間	1区間、30ノード、終了時刻は自由
新しく学ぶこと	2次元の極座標で書く運動方程式、推力を2成分で表す方法、多数の経路制約、最終質量の最大化

状況

地上で静止している質量 100 t のロケットを、高度 400 km の円軌道に入れます。円軌道に乗るには、終端で

• 高度がちょうど 400 km であること
• 上下方向の速度が 0 であること
• 水平方向の速度が、その高度で地球を回り続けられる速さ（円軌道速度、約 7.67 km/s）であること

の3つが必要です。

機体の 99 % は推進剤で、空になったときの質量は 1 t です。推進剤をどれだけ残して軌道に着けるかが性能の指標になり、 これを最大にします。推力の大きさだけでなく向きも自由に選べるのが、例題10 との違いです。

途中では、空気抵抗を受けること、推力に上限があること、機体にかかる加速度が 8 G を超えないことを守らせます。 地球の自転は考えません。

問題設定

座標系

ロケットは赤道面のような1つの平面の中を飛ぶとし、位置を地球中心からの極座標で表します。

• $R$：地球中心からの距離。高度は $h = R - R_e$（$R_e$ は地球半径）です。
• $\theta$：射点から地球中心のまわりに回った角度。地表に沿った距離（ダウンレンジ）は $R_e\theta$ です。
• $V_r$：$R$ が増える向き、つまり真上向きの速度（鉛直速度）。
• $V_t$：それに直角で、進行方向を向く速度（水平速度）。
• $T_r,\ T_t$：推力の鉛直成分と水平成分。推力の大きさは $|T| = \sqrt{T_r^2 + T_t^2}$、 水平から測った推力の向き（仰角）は $\arctan(T_r / T_t)$ です。

推力を「大きさと角度」ではなく「2つの成分」で持つのは、角度を使うと $\sin$・$\cos$ が入って最適化が難しくなり、 推力 0 のとき角度が決まらないという問題も生じるからです。

運動方程式

$$\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}} \end{aligned}$$

極座標で書くと、速度の式に2つの項が加わります。

• $V_t^2/R$ は遠心力にあたる項です。水平速度が大きいほど、機体は地球から離れる向きに押されます。 $V_t^2/R = g$ となる速度が円軌道速度で、このとき重力と遠心力がつり合って高度が変わらなくなります。
• $-V_r V_t/R$ は、上昇すると水平速度が自然に減ることを表します。推力も抵抗もなければ角運動量 $R V_t$ が保存されるので、 $R$ が増えれば $V_t$ は減ります。

空気抵抗は速度ベクトルと逆向きにかかり、その成分は次のとおりです。

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

空気密度と重力は高度の関数です。

$$\rho(h) = 1.225\, e^{-h/8500}, \qquad g(R) = \frac{\mu}{R^2}$$

密度は高度 8.5 km ごとに $1/e$ 倍になる指数関数で近似しています。重力は地球の重力定数 $\mu$ を距離の2乗で割った形で、終端の円軌道速度と同じ $\mu$ を使います。 質量は推力の大きさに比例して減り、比推力 $I_{sp}$ が大きいほど推進剤を節約できます。

境界条件

$$\begin{aligned} &R(0) = R_e,\quad \theta(0) = 0 \\ &V_r(0) = 0,\quad V_t(0) = 0 \\ &m(0) = 100\ \mathrm{t} \\[4pt] &R(t_f) = R_e + 400\ \mathrm{km} \\ &V_r(t_f) = 0 \\ &V_t(t_f) = \sqrt{\mu / R(t_f)} \end{aligned}$$

地表で静止した状態から出発します。終端では高度 400 km、鉛直速度 0、水平速度は円軌道速度です。 $\mu$ は地球の重力定数 $3.986 \times 10^{14}\ \mathrm{m^3/s^2}$ で、円軌道速度は約 7672.6 m/s になります。 終端の角度 $\theta(t_f)$（どこで軌道に入るか）と終了時刻 $t_f$ は自由です。

制約

$$\begin{aligned} &R \ge R_e \\ &0 \le T_r \le T_{\max},\quad 0 \le T_t \le T_{\max} \\ &\sqrt{T_r^2 + T_t^2} \le T_{\max} \\ &1\ \mathrm{t} \le m \le 100\ \mathrm{t} \\ &\left|\frac{\boldsymbol{T} - \boldsymbol{D}}{m}\right| \le 8\, g_0 \end{aligned}$$

• 地面より下には行けません。
• 推力の成分はどちらも 0 以上です。つまり推力を下向きや後ろ向きには出しません。
• 推力の大きさは最大推力 $T_{\max}$（初期重量の 1.5 倍、約 1471 kN）以下です。
• 質量は空の機体（1 t）を下回れません。推進剤の残量に対する制約です。
• 推力と抵抗による加速度の大きさを 8 G 以下にします。乗っているものが感じる加速度（加速度計が測る値）なので、重力は含めません。

これらは各ノードで課されます。

目的関数

$$J = -\frac{m(t_f)}{m(0)}$$

最終質量を最大にしたいので、マイナスをつけて最小化します。初期質量で割って 1 程度の大きさにそろえています。

パラメータ

記号	意味	値
$R_e$	地球半径	6371 km
$g_0$	標準重力加速度（比推力・最大推力・加速度の上限の換算に使う）	9.80665 m/s²
$\mu$	地球の重力定数（重力と目標速度の計算に使う）	$3.986 \times 10^{14}$ m³/s²
$m(0)$	初期質量	100 t
—	推進剤の質量（空の機体は 1 t）	99 t
$C_D A$	抵抗係数 × 断面積	0.6 × 4 m²
$I_{sp}$	比推力	300 s
$T_{\max}$	最大推力（$1.5\, m(0)\, g_0$）	約 1471 kN
—	加速度の上限	8 G
—	目標高度	400 km

$R$ は数千 km、質量は 100 t、推力は MN の大きさで、そのままでは桁がそろいません。そこで距離は $R_e$、速度は $\sqrt{\mu/R_e}$、 質量は $m(0)$ を単位にして、最適化の変数を 1 程度にそろえています（単位をそろえる）。

初期推定

最初の区間の長さは 200 秒とし、次のように与えます。

• 距離 $R$：地表から目標の高さへなめらかに増える3次式。
• 角度 $\theta$：0 から 25° へなめらかに増える3次式。
• 鉛直速度 $V_r$ は 0、水平速度 $V_t$ は 0 から目標速度へ直線的に増やします。
• 質量：100 t から 1 t へ減る3次式。
• 推力：鉛直・水平成分とも $T_{\max}/2$ から 0 へ下がる3次式。

結果

終端の高度は 400.0 km（目標 400.0 km）、 終端の速さは 7,672.6 m/s（目標 7,672.6 m/s）で、 境界条件を満たしています。

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

時間履歴。高度、水平速度（実線）と鉛直速度（破線）、質量、推力の大きさ、推力の仰角（水平から測った角度、エンジン停止中は非表示）、推力と抵抗による加速度。赤い点線は最大推力と加速度の上限です。

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

• 打上げから 263 秒ごろまで燃焼し、いったんエンジンを止めます。
• その後は惰性で高度 400 km まで上がり、時刻 886 秒ごろから 901 秒まで再び燃焼して水平速度を円軌道速度まで上げます。
• 軌道に残る質量は 4,370 kg で、初期質量の約 4 % です。 そのうち 1 t は空の機体で、残りは使い残した推進剤です。
• 軌道に入る位置は、射点から地表に沿って 5,105 km 先です。
• 加速度は機体が軽くなる燃焼の終わりごろに大きくなり、最大は円軌道化の燃焼の終わりで 8.00 G（ノードでの値）です。 上限の 8 G には達しておらず、この解では加速度の制約は効いていません。

解説

立ち上がって、倒れる

推力の仰角は最初ほぼ真上（約 80°）を向き、時間とともに水平へ倒れていきます。 空気の濃い低空では、速く飛ぶと抵抗で損をするので、まず上へ出て空気の薄いところへ抜けます。 一方、軌道に必要なのは水平速度なので、上へ向け続けると重力に逆らうことに推進剤を使いすぎます（重力損失）。 両者の兼ね合いで、早めに機体を倒して水平方向に加速するのが得になります。実際のロケットの上昇（重力ターン）とよく似た形です。

燃焼・慣性飛行・燃焼

最適な飛び方は「加速 → 慣性飛行 → 円軌道化の燃焼」の3段階になっています。最初の燃焼で、高度 400 km を最高点とする 楕円の軌道に乗る速度まで低い高度で加速し、あとは推進剤を使わずに最高点まで上がります。そこで水平に加速して円軌道にします。

低い円軌道から高い円軌道へ移るホーマン遷移と同じ考え方です。高く上がるための仕事は慣性飛行が「ただで」してくれるので、 燃焼は低い高度での加速と最高点での加速に集中させるのが効率的です。慣性飛行中、上昇するにつれて水平速度が少しずつ減っていくのは、 運動方程式の $-V_r V_t/R$ の項（角運動量の保存）の効果です。

図を読むときの注意

• 制約はノードでしか課していないので、ノードの間を結ぶ多項式は推力の上限を守るとは限りません。加速度の線にも同じことが言えます。
• 燃焼と停止の切り替わりを1本の多項式で表しているので、切り替わり付近の推力の細かい形はノードの数に左右されます （区間に分ける）。結果を確かめる手順も参照してください。

試してみよう

• ノード数（Problem([0.0, 200], [30], ...) の 30）や初期推定の区間長 200 を変え、ソルバーが収束を報告するか、最終質量がどう変わるかを調べる。
• Rocket の MaxG を 8 から 4 に下げ、加速度の制約が効いて推力を絞るようになる様子と、最終質量の変化を見る。
• Htarget を 400 km から 300 km や 600 km に変え、慣性飛行の長さと最終質量を比べる。

実行する

uv run python examples/12_rocket_ascent_polar_ssto.py

examples/12_rocket_ascent_polar_ssto.py を GitHub で見る