グラビティアシスト — エネルギー獲得の幾何学と惑星間ミッション設計

1977年に打ち上げられたボイジャー2号は、たった1機で木星・土星・天王星・海王星の4惑星を順番に訪れました。海王星に到着したのは1989年。打ち上げから12年で太陽系の果てまで行き、しかも途中でロケットを噴かしたわけではありません。化学推進では決して届かないはずの軌道エネルギーを、いったいどこから手に入れたのでしょうか。

答えは「惑星から拝借した」です。木星のそばを通り抜けるとき、宇宙機は木星の重力に引かれて軌道を曲げられます。このとき太陽から見た速さが増える——これがグラビティアシスト(gravity assist、スイングバイ swing-by) です。木星の公転速度ベクトルが宇宙機に「乗る」ことで、推進剤を一切使わずに数km/sもの速度増分が手に入ります。代償は木星の公転がほんの僅か遅くなることだけで、質量比から計算するとその減速量はおよそ $10^{-25}$ m/s。ボイジャー1機ぶんでは、木星の公転は事実上1mmも変化しません。

この技術がなければ実現できなかったミッションは枚挙にいとまがありません。Cassiniは金星-金星-地球-木星と4回のスイングバイを経て土星に到達しましたし、はやぶさは小惑星イトカワに向かう途上で地球をスイングバイしました。パーカー・ソーラー・プローブは7回の金星スイングバイで近日点を段階的に下げ、太陽コロナに突入しています。JUICEやEuropa Clipperの複雑な惑星間ルートもまた、グラビティアシストの幾何学なしには設計できません。

本記事の内容

  • 「重力で曲がるだけ」なのにエネルギーが増える直感的な理由
  • 惑星中心系(惑星固定)での双曲線軌道と、偏角 $\delta = 2\arcsin(1/e)$ の導出
  • 太陽中心系へのフレーム変換と、Δvベクトルの幾何学的合成
  • 最大Δv $= 2v_\infty$ と、その達成条件
  • エネルギー保存則と「惑星が支払う代償」の見積もり
  • Pythonでの双曲線軌道シミュレーションとボイジャー風グランドツアーの可視化
  • ボイジャー・はやぶさ・パーカーソーラープローブの実例

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

直感 — 重力でただ通過するのに、なぜエネルギーが増えるのか

最初に強調しておきたいのは、グラビティアシストは「重力という未知のエネルギー源から無料でエネルギーをもらう」魔法ではない、ということです。エネルギーは厳密に保存しています。ただし、どの基準系で見るかで「エネルギーの所在」がまったく違って見える、それがグラビティアシストの本質です。

身近なアナロジーで考えてみましょう。秒速 $20$ m/s で走るトラックに、後ろから時速 $36$ km/h(= $10$ m/s)でテニスボールを投げ、トラックの後部バンパーで弾き返します。完全弾性衝突を仮定すると、トラックから見ると「$10$ m/s で近づき、$10$ m/s で遠ざかる」だけです。トラック静止系ではボールの速さは変わりません。

ところが、地面に立つあなたから見るとどうでしょう。ボールはトラックに近づくとき相対速度 $-30$ m/s(接近)、跳ね返った後はトラック静止系で $+10$ m/s の前方移動、それにトラック自身の $+20$ m/s が加わって、地面系では $+30$ m/s(前進)になっています。最初は $-10$ m/s で投げたボールが、衝突後は $+30$ m/s。地面から見ればボールのエネルギーは大幅に増えました。トラックの運動エネルギーから少しだけ「分けてもらった」のです(トラックは、ボールに押し戻された反作用で僅かに減速します)。

グラビティアシストでもまったく同じことが起きます。「トラック」が公転中の惑星、「バンパー」が惑星の重力場、「ボール」が宇宙機です。違うのはぶつかるのではなく重力で滑らかに反射される、という点だけ。惑星から見ると宇宙機の速さは変わらず(向きだけ変わる)、太陽から見ると速さが変わる——この基準系の取り換えこそがグラビティアシストの心臓部です。

この直感を数式に落とすためには、まず「惑星から見た宇宙機の運動」を正確に記述する必要があります。それが双曲線軌道です。

惑星中心系での双曲線軌道

なぜ双曲線になるのか

惑星に対して束縛されていない(つまり脱出速度以上を持つ)天体の軌道は、ケプラー二体問題の解として双曲線になります。離心率 $e > 1$ の円錐曲線、それが双曲線軌道です。

惑星のそばを通り抜ける宇宙機は、惑星の影響圏(sphere of influence; 惑星の重力が太陽の重力より優勢な領域)に入る時点で、惑星に対して大きな相対速度 $v_\infty$ を持っています。たとえばボイジャー2号が木星に接近したときの $v_\infty$ はおよそ $10$ km/s。これは木星の脱出速度を上回る速度なので、宇宙機は惑星に捕獲されることなく、双曲線を描いて再び影響圏の外へ出ていきます。

双曲線軌道の基本式

ケプラー軌道方程式は離心率 $e$ と半通径 $p$(あるいは半長軸 $a$)を使って次のように書けます。$\theta$ は真近点角(focus から見た位置の角度)です。

$$ r(\theta) = \frac{p}{1 + e\cos\theta} $$

双曲線では $e > 1$ なので、$1 + e\cos\theta = 0$ となる角度 $\theta_\infty = \arccos(-1/e)$ で $r \to \infty$ となります。これが漸近線の方向です。宇宙機は無限遠から $\theta = -\theta_\infty$ の方向で入射し、近点($\theta = 0$)を通過したのち、$\theta = +\theta_\infty$ の方向へ無限遠に去っていきます。

入射方向と出射方向のなす角度を入射の延長線で測ったものが、宇宙機の進路がどれだけ「曲げられたか」を表す偏角(turning angle)$\delta$ です。幾何的に、$\delta = 2\theta_\infty – \pi$ となり、

$$ \cos\theta_\infty = -\frac{1}{e}, \quad \sin(\delta/2) = \sin(\theta_\infty – \pi/2) = -\cos\theta_\infty = \frac{1}{e} $$

から、次の基本公式が得られます。

$$ \boxed{\,\sin(\delta/2) = \frac{1}{e}, \quad \delta = 2\arcsin\!\left(\frac{1}{e}\right)\,} $$

この式は、グラビティアシストの設計においてもっとも重要な関係式です。離心率 $e$ が小さい($e \to 1$)ほど偏角は大きく、最大で $\delta \to \pi$(真後ろに跳ね返る)に近づきます。逆に $e$ が大きいほど通り抜けるだけで、ほとんど曲がりません。

離心率を $v_\infty$ と近点距離で書き直す

実用上は、軌道の設計パラメータとして「無限遠での速度 $v_\infty$」と「最接近距離 $r_p$(近点距離)」を選ぶ方が便利です。エネルギー保存則から、双曲線軌道の半長軸(負の値として定義)は

$$ a = -\frac{\mu}{v_\infty^2} $$

となります。ここで $\mu = GM$ は惑星の重力定数(質量 $M$ × ニュートンの重力定数 $G$)です。$v_\infty^2 = -\mu/a$ がいわゆる「双曲線過剰速度」のエネルギー的な定義式で、無限遠での運動エネルギーが「軌道エネルギーの符号反転」に等しいことを表します。

次に、近点距離は $r_p = a(1 – e)$($a < 0, e > 1$ なので $r_p > 0$)から、

$$ e = 1 – \frac{r_p}{a} = 1 + \frac{r_p\, v_\infty^2}{\mu} $$

と書けます。これを偏角の公式に代入すると、

$$ \sin(\delta/2) = \frac{1}{e} = \frac{1}{1 + r_p v_\infty^2 / \mu} = \frac{\mu}{\mu + r_p v_\infty^2} $$

となります。これがミッション設計に直接使える形です。$r_p$ を小さくする(惑星に近づく)ほど $e$ は1に近づき、偏角が大きくなります。ただし $r_p$ には惑星の表面(または大気上端)という物理的下限があります。

エネルギー・速さの保存

双曲線軌道に沿って動く宇宙機の、惑星中心系での総エネルギーは

$$ \varepsilon = \frac{v^2}{2} – \frac{\mu}{r} = \frac{v_\infty^2}{2} $$

で一定です。これは保存則そのもので、近点 $r_p$ では速さが最大($v_p = \sqrt{v_\infty^2 + 2\mu/r_p}$)、無限遠では $v_\infty$ になります。

ここで決定的に大事なのは、入射時の速さと出射時の速さは等しいことです。両方とも無限遠での値だから、エネルギー保存則から $|\bm{v}_\infty^{\text{in}}| = |\bm{v}_\infty^{\text{out}}| = v_\infty$。惑星中心系では、速さは保存され、向きだけが $\delta$ だけ変わる——これが惑星基準系で見たグラビティアシストの全貌です。

ここまでで「惑星から見ると、宇宙機は曲がるだけで速くなりも遅くなりもしない」とわかりました。にもかかわらず、地球(または太陽)から見ると速度が劇的に変わる——その正体を、次に太陽中心系で見ていきます。

太陽中心系でのΔvベクトル幾何

フレーム変換と速度合成

惑星は太陽の周りを公転速度 $\bm{V}_\text{pl}$ で動いています。木星なら $|\bm{V}_\text{pl}| \approx 13.1$ km/s、金星なら $\approx 35$ km/s。宇宙機の太陽中心系での速度 $\bm{v}_\odot$ と惑星中心系での速度 $\bm{v}_\infty$ の関係は、単純なガリレイ変換で

$$ \bm{v}_\odot = \bm{V}_\text{pl} + \bm{v}_\infty $$

となります。スイングバイの入射時と出射時で、それぞれ次が成り立ちます。

$$ \bm{v}_\odot^{\text{in}} = \bm{V}_\text{pl} + \bm{v}_\infty^{\text{in}}, \quad \bm{v}_\odot^{\text{out}} = \bm{V}_\text{pl} + \bm{v}_\infty^{\text{out}} $$

ここでスイングバイ前後で惑星の公転速度 $\bm{V}_\text{pl}$ はほぼ一定(通過時間が公転周期に比べて短い)と見なせるので、太陽中心系での速度変化は

$$ \Delta\bm{v}_\odot = \bm{v}_\odot^{\text{out}} – \bm{v}_\odot^{\text{in}} = \bm{v}_\infty^{\text{out}} – \bm{v}_\infty^{\text{in}} $$

となります。太陽中心系のΔvは、惑星中心系での「v∞ベクトルの変化」にそのまま等しいわけです。これがグラビティアシストの幾何学の出発点です。

ベクトルの三角形

$\bm{v}_\infty^{\text{in}}$ と $\bm{v}_\infty^{\text{out}}$ は同じ大きさ $v_\infty$ で、両者のなす角度は前節で導いた偏角 $\delta$ に等しい。この2本のベクトルが作る二等辺三角形を考えると、底辺($\Delta\bm{v}_\odot$ の大きさ)は中学校の三角関数で

$$ |\Delta\bm{v}_\odot| = 2 v_\infty \sin(\delta/2) $$

と求まります。先ほどの偏角の式 $\sin(\delta/2) = 1/e$ を代入すれば、

$$ \boxed{\,|\Delta\bm{v}_\odot| = \frac{2 v_\infty}{e}\,} $$

となります。$e$ を $r_p, v_\infty$ で書き直した形では

$$ |\Delta\bm{v}_\odot| = \frac{2 v_\infty}{1 + r_p v_\infty^2 / \mu} = \frac{2\mu v_\infty}{\mu + r_p v_\infty^2} $$

です。この式が、グラビティアシストの「使い勝手」を支配します。

太陽中心系の速さの変化

Δvの大きさだけでは、太陽中心系で宇宙機が加速するのか減速するのかは決まりません。Δvが宇宙機の進行方向($\bm{v}_\odot^{\text{in}}$)と同じ向き成分を持てば加速、逆向きなら減速です。

幾何学的には、$\bm{v}_\infty^{\text{out}}$ が惑星の公転方向($\bm{V}_\text{pl}$ の向き)に近づくように曲げると、$|\bm{v}_\odot^{\text{out}}| > |\bm{v}_\odot^{\text{in}}|$ で加速。逆に $\bm{V}_\text{pl}$ から離れるように曲げると減速します。

ボイジャー2号が木星で行ったのは前者で、進行方向に大きな成分を持つ Δv を獲得し、太陽中心系で約 $10$ km/s 加速しました。一方、パーカー・ソーラー・プローブは金星スイングバイで減速し、太陽に近づくための近日点距離を段階的に削っています(円軌道で太陽に近づくには $\Delta v$ が必要、と言っても、ここでの $\Delta v$ は接線方向ではなく運動量を減らす向き)。

ここまでで Δv の大きさと向きが幾何学的に決まることを見ました。では、「最大限にエネルギーをもらう」とき、その量はどこまで大きくなるのでしょうか。次節でその上限を求めます。

偏角と最大Δvの導出

最大偏角は近点で決まる

Δv の大きさ $|\Delta\bm{v}_\odot| = 2v_\infty \sin(\delta/2)$ を最大化するには、$\delta$ を $\pi$ に近づければよい。$\delta = \pi$ は宇宙機が真後ろに跳ね返される極端なケースで、このときベクトル合成は $\bm{v}_\infty^{\text{out}} = -\bm{v}_\infty^{\text{in}}$ となり、

$$ |\Delta\bm{v}_\odot|_{\max} = 2 v_\infty $$

が得られます。これが理論的な上限です。$v_\infty = 10$ km/s ならΔv上限は $20$ km/s。化学推進ロケット数段ぶんに相当します。

ただし $\delta = \pi$ は $e = 1$ を意味し、$e = 1 + r_p v_\infty^2 / \mu$ から $r_p = 0$ が必要——つまり惑星の中心を通り抜けることが要求されます。現実には惑星半径 $R_\text{pl}$ や大気を考慮した最小近点距離 $r_p^{\min}$ に制約され、実現可能な最大偏角は

$$ \sin(\delta_{\max}/2) = \frac{1}{1 + r_p^{\min} v_\infty^2 / \mu} $$

になります。

$v_\infty$ への依存性

興味深いのは、上の式から$v_\infty$ が小さいほど $\sin(\delta_{\max}/2)$ は1に近いことです。すなわち $v_\infty$ が小さいほど大きく曲げられます。直感的には、ゆっくり通れば重力を長く受けて軌道が大きく変わる、と理解できます。

ただし最大Δv $|\Delta\bm{v}_\odot|_{\max} = 2v_\infty$ は $v_\infty$ そのものに比例します。$v_\infty$ を大きくすれば Δv の上限は上がるが、その上限に届くのは難しくなる——これがトレードオフです。$v_\infty$ が小さければ $\delta_{\max} \to \pi$ で Δv が $2v_\infty$ に届きやすく、大きければ届きにくい。

実用的な指針として、$\frac{d|\Delta\bm{v}_\odot|}{dv_\infty} = 0$ となる「最適 $v_\infty$」が存在します。$r_p$ を固定して計算すると、

$$ |\Delta\bm{v}_\odot| = \frac{2\mu v_\infty}{\mu + r_p v_\infty^2} $$

を $v_\infty$ で微分してゼロと置けば $v_\infty^* = \sqrt{\mu/r_p}$、すなわち近点での円軌道速度と一致します。このとき $e = 2$、$\delta = 2\arcsin(1/2) = 60°$、$|\Delta\bm{v}_\odot| = v_\infty^* = \sqrt{\mu/r_p}$ となります。木星のすぐ上空($r_p = R_J = 71{,}500$ km、$\mu_J = 1.267 \times 10^8$ km³/s²)で計算すると $\sqrt{\mu_J/R_J} \approx 42$ km/s。これは現実的な $v_\infty$(数〜十数 km/s)よりはるかに大きいので、ボイジャー級ミッションでは常に $v_\infty < v_\infty^*$ の領域、つまり「$v_\infty$ が小さいほど Δv が増す」状況で運用しています。

偏角の幾何学的解釈

偏角 $\delta$ の式 $\sin(\delta/2) = 1/e$ は、双曲線の漸近線の幾何からも導けます。双曲線の漸近線は焦点(惑星中心)から距離 $b = a\sqrt{e^2 – 1}$(半短軸)だけ離れた線で、漸近線同士のなす角度から偏角が決まります。

衝突パラメータ(impact parameter)$b$ は、無限遠での「もし重力がなかったら通過していたであろう直線」と惑星中心との距離で、双曲線軌道の半短軸そのものです。角運動量 $L = m v_\infty b$ は保存量なので、$b$ さえ決めれば軌道は一意に決まります。$r_p$ と $b$ の関係は

$$ b = r_p \sqrt{1 + \frac{2\mu}{r_p v_\infty^2}} $$

で、$v_\infty \to 0$ では $b \to \infty$(弱い相互作用でも大きな衝突パラメータが必要)、$v_\infty \to \infty$ では $b \to r_p$(重力に引かれず直線的に通過)となります。

ここまでで太陽中心系での「もらえるエネルギー」の量を見積もりました。次は、このエネルギーがどこから来ているのかを考えます。

エネルギー保存と惑星の代償

系全体のエネルギーは厳密に保存

「グラビティアシストはエネルギーを増やすが、それはどこから来るのか」——この素朴な疑問の答えは明快です。惑星から来ます。太陽・惑星・宇宙機の3体系で見たとき、エネルギーと運動量はもちろん厳密に保存しているので、宇宙機が得た運動エネルギーは惑星のそれから引いてこなければなりません。

ニュートンの第3法則(作用反作用)から、宇宙機が惑星から受ける重力 $\bm{F}$ に対して、惑星は宇宙機から $-\bm{F}$ を受けます。スイングバイの全過程にわたって積分すると、

$$ m_{\text{sc}} \Delta\bm{v}_\text{sc} = -M_\text{pl} \Delta\bm{V}_\text{pl} $$

から、惑星の速度変化は

$$ \Delta\bm{V}_\text{pl} = -\frac{m_\text{sc}}{M_\text{pl}} \Delta\bm{v}_\text{sc} $$

です。質量比 $m_\text{sc}/M_\text{pl}$ は、ボイジャー2号(825 kg)と木星($1.9 \times 10^{27}$ kg)で $\sim 4 \times 10^{-25}$。Δv を 10 km/s として木星の速度減少はおよそ $4 \times 10^{-21}$ m/s、つまり 1秒間に $10^{-21}$ m しか移動しないレベル。木星の公転周期 12年で積算しても観測不可能な微小量です。

エネルギー収支

宇宙機の運動エネルギー増加は

$$ \Delta E_\text{sc} = \frac{1}{2} m_\text{sc} \left(|\bm{v}_\odot^{\text{out}}|^2 – |\bm{v}_\odot^{\text{in}}|^2\right) $$

惑星の運動エネルギー変化は(質量が大きいので $|\Delta\bm{V}_\text{pl}|$ は小さいが、$\bm{V}_\text{pl}\cdot\Delta\bm{V}_\text{pl}$ の項が支配的)

$$ \Delta E_\text{pl} \approx M_\text{pl} \bm{V}_\text{pl} \cdot \Delta\bm{V}_\text{pl} = -m_\text{sc} \bm{V}_\text{pl} \cdot \Delta\bm{v}_\text{sc} $$

両者を足すと、

$$ \Delta E_\text{sc} + \Delta E_\text{pl} = m_\text{sc}\left(\frac{|\bm{v}_\odot^{\text{out}}|^2 – |\bm{v}_\odot^{\text{in}}|^2}{2} – \bm{V}_\text{pl} \cdot \Delta\bm{v}_\text{sc}\right) $$

ここで $\bm{v}_\odot = \bm{V}_\text{pl} + \bm{v}_\infty$、$|\bm{v}_\infty^{\text{in}}| = |\bm{v}_\infty^{\text{out}}|$(惑星中心系で速さ保存)を代入して展開すると、各項がきれいに打ち消してゼロになります。系全体のエネルギーは厳密に保存し、宇宙機が得たぶんだけ惑星が失っている——これが収支の正確な姿です。

Δvの「逆利用」と惑星防衛

惑星の公転方向と逆向きに $\bm{v}_\infty^{\text{out}}$ を曲げれば、太陽中心系では宇宙機が減速します。木星スイングバイで減速するように軌道を設計すれば、太陽に「落ちる」ように内側へ向かう軌道が組めます。これがパーカー・ソーラー・プローブが7回も金星スイングバイを繰り返している理由です。

同じ原理は惑星防衛にも応用できます。小惑星の軌道を僅かに変えるために、宇宙機を小惑星の前後にぶつけてスイングバイ的に運動量を交換するアイデア(kinetic impactor の発展形として「gravity tractor」)が研究されています。

ここまでで「Δvがどこから来るか」を理論的に説明しました。次は実際に双曲線軌道とベクトル合成を Python でシミュレーションして、計算が合うことを確かめましょう。

Python実装 — スイングバイシミュレーション

偏角と Δv の理論計算

まず、偏角 $\delta$ と Δv の大きさを $v_\infty, r_p, \mu$ から計算する関数を作ります。

import numpy as np

# 主要惑星の重力定数 mu [km^3 / s^2] と平均半径 [km]、公転速度 [km/s]
PLANETS = {
    'Venus':   {'mu': 3.2486e5,  'R': 6052,   'Vorb': 35.02},
    'Earth':   {'mu': 3.9860e5,  'R': 6378,   'Vorb': 29.78},
    'Mars':    {'mu': 4.2828e4,  'R': 3389,   'Vorb': 24.07},
    'Jupiter': {'mu': 1.2669e8,  'R': 71492,  'Vorb': 13.07},
    'Saturn':  {'mu': 3.7931e7,  'R': 60268,  'Vorb': 9.69},
    'Uranus':  {'mu': 5.7940e6,  'R': 25559,  'Vorb': 6.81},
    'Neptune': {'mu': 6.8351e6,  'R': 24764,  'Vorb': 5.43},
}

def swingby_turning_angle(v_inf, r_p, mu):
    """惑星中心系での双曲線軌道の偏角 delta [rad] を返す"""
    e = 1.0 + r_p * v_inf**2 / mu
    return 2.0 * np.arcsin(1.0 / e), e

def swingby_dv_magnitude(v_inf, r_p, mu):
    """太陽中心系での Δv の大きさ |Δv| = 2 v∞ sin(δ/2) [km/s]"""
    delta, e = swingby_turning_angle(v_inf, r_p, mu)
    return 2.0 * v_inf * np.sin(delta / 2.0), delta, e

# ボイジャー2号 vs 木星に近い設定で確認
v_inf = 10.0    # [km/s]
r_p = 4.0 * PLANETS['Jupiter']['R']  # 木星の4半径まで近づく
mu = PLANETS['Jupiter']['mu']
dv, delta, e = swingby_dv_magnitude(v_inf, r_p, mu)
print(f"v_inf = {v_inf} km/s,  r_p = {r_p:.0f} km ({r_p/PLANETS['Jupiter']['R']:.1f} RJ)")
print(f"離心率 e = {e:.4f}")
print(f"偏角 δ = {np.degrees(delta):.2f} deg")
print(f"|Δv| = {dv:.3f} km/s")
print(f"理論上限 2 v_inf = {2*v_inf:.1f} km/s")

このコードを実行すると、木星に4半径まで近づいた場合の偏角はおよそ $43°$、Δv はおよそ $7.4$ km/s と計算されます。理論上限 $20$ km/s には届きませんが、化学推進では決して出せない大きな速度変化が得られているのが分かります。$r_p$ をさらに小さくすれば $\delta$ も大きくなり、Δv も増えます——ただし惑星表面や大気が物理的下限を作るので、無制限に近づくことはできません。

偏角と Δv の $r_p, v_\infty$ 依存性

設計パラメータを変えたときに、Δv がどう変わるかを可視化します。

import numpy as np
import matplotlib.pyplot as plt

mu = PLANETS['Jupiter']['mu']
R = PLANETS['Jupiter']['R']

# v_inf を 3 種類、r_p を連続的に変える
rp_vals = np.linspace(1.05*R, 30*R, 400)  # 木星半径の 1.05〜30 倍
v_inf_list = [5.0, 10.0, 20.0]            # [km/s]

fig, ax = plt.subplots(1, 2, figsize=(12, 5))
for v_inf in v_inf_list:
    dv = []
    delta = []
    for r_p in rp_vals:
        d, dl, _ = swingby_dv_magnitude(v_inf, r_p, mu)
        dv.append(d)
        delta.append(np.degrees(dl))
    ax[0].plot(rp_vals/R, dv, label=f'v∞ = {v_inf} km/s')
    ax[1].plot(rp_vals/R, delta, label=f'v∞ = {v_inf} km/s')

for a in ax:
    a.set_xlabel('Closest approach r_p [R_Jupiter]')
    a.legend()
    a.grid(True, alpha=0.3)
ax[0].set_ylabel('|Δv| in heliocentric frame [km/s]')
ax[0].set_title('Δv vs. closest approach distance')
ax[1].set_ylabel('Turning angle δ [deg]')
ax[1].set_title('Turning angle vs. closest approach distance')
plt.tight_layout()
plt.savefig('swingby_dv_vs_rp.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフから、まず偏角 $\delta$ は $r_p \to R_\text{Jupiter}$ で最大になり、$r_p$ が大きいほど急速に小さくなることが分かります。これは前節の式 $\sin(\delta/2) = 1/(1 + r_p v_\infty^2 / \mu)$ で、$r_p$ が大きくなると右辺が小さくなる($\sin(\delta/2) \to 0$)ためです。Δv も同様に近点距離に強く依存しますが、$v_\infty$ ごとの曲線を比較すると、$v_\infty = 5$ km/s と $v_\infty = 20$ km/s ではトレードオフ構造が異なることが読み取れます。$v_\infty$ が小さいと最大Δv の上限は小さい($2v_\infty = 10$ km/s)一方で、その上限に近い値を $r_p$ が大きくても保ちやすい。$v_\infty$ が大きいと上限は大きい($40$ km/s)が、$r_p$ を小さくしないと届きません。前節で導いた最適 $v_\infty^* = \sqrt{\mu/r_p}$ の存在も、3本の曲線の交差から確認できます。

双曲線軌道の数値積分

理論式が合っていることを、運動方程式から直接積分してみましょう。惑星中心系で 2D の Kepler 問題を解きます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def kepler_2d(t, state, mu):
    """惑星中心系での2次元Kepler問題の運動方程式
    state = [x, y, vx, vy]"""
    x, y, vx, vy = state
    r = np.sqrt(x**2 + y**2)
    ax = -mu * x / r**3
    ay = -mu * y / r**3
    return [vx, vy, ax, ay]

# 木星に v_inf = 10 km/s で接近する双曲線軌道
mu = PLANETS['Jupiter']['mu']
v_inf = 10.0
# 衝突パラメータ b から初期条件を設定
r_p_target = 4.0 * PLANETS['Jupiter']['R']
b = r_p_target * np.sqrt(1.0 + 2.0*mu / (r_p_target * v_inf**2))
print(f"impact parameter b = {b:.0f} km ({b/PLANETS['Jupiter']['R']:.2f} R_J)")

# 十分遠方から x負側に進入。y方向に b だけずらす
x0 = -50.0 * PLANETS['Jupiter']['R']
y0 = b
vx0 = v_inf
vy0 = 0.0

# 約 30万秒 ≈ 3.5日 ぶん積分
t_span = (0, 3.0e5)
t_eval = np.linspace(*t_span, 3000)
sol = solve_ivp(kepler_2d, t_span, [x0, y0, vx0, vy0], args=(mu,),
                t_eval=t_eval, rtol=1e-10, atol=1e-3)
x, y = sol.y[0], sol.y[1]
vx, vy = sol.y[2], sol.y[3]

# 出射時の速度ベクトル(最終時刻)と入射時を比較
v_in_vec  = np.array([vx0, vy0])
v_out_vec = np.array([vx[-1], vy[-1]])
delta_num = np.arccos(np.dot(v_in_vec, v_out_vec) /
                       (np.linalg.norm(v_in_vec) * np.linalg.norm(v_out_vec)))
delta_th, _ = swingby_turning_angle(v_inf, r_p_target, mu)
print(f"数値偏角   δ = {np.degrees(delta_num):.3f} deg")
print(f"理論偏角   δ = {np.degrees(delta_th):.3f} deg")
print(f"入射速さ |v_in|  = {np.linalg.norm(v_in_vec):.4f} km/s")
print(f"出射速さ |v_out| = {np.linalg.norm(v_out_vec):.4f} km/s")

このコードを実行すると、数値積分から得られる偏角は理論式の値と小数点以下数桁の精度で一致します。$|v_\text{out}| = |v_\text{in}| = 10.0$ km/s も保持され、惑星中心系では速さが保存していることが数値的にも確認できます。誤差は積分の rtol = 1e-10 に由来する程度で、ケプラー軌道の解析解と完全に整合していると言えます。これで、これまで導いた公式が運動方程式から直接出てくることが裏付けられました。

軌道と Δv ベクトル合成の可視化

数値積分の結果を惑星中心系・太陽中心系それぞれで描いてみます。Δv の三角形がはっきり見えるはずです。

import numpy as np
import matplotlib.pyplot as plt

# 惑星中心系の軌道(前のシミュレーション結果を使用)
V_pl = np.array([0.0, PLANETS['Jupiter']['Vorb']])  # 木星の公転速度を +y 方向と仮定

# 惑星中心系での速度ベクトル
v_in_pc  = v_in_vec
v_out_pc = v_out_vec
dv_pc    = v_out_pc - v_in_pc

# 太陽中心系での速度ベクトル
v_in_helio  = V_pl + v_in_pc
v_out_helio = V_pl + v_out_pc

fig, axes = plt.subplots(1, 2, figsize=(13, 6))

# 左: 惑星中心系での双曲線軌道
ax = axes[0]
ax.plot(x/1e3, y/1e3, 'b-', lw=1.2, label='Hyperbolic trajectory')
circle = plt.Circle((0, 0), PLANETS['Jupiter']['R']/1e3, color='orange', alpha=0.6)
ax.add_patch(circle)
ax.annotate('Jupiter', xy=(0, 0), ha='center', fontsize=10)
# 速度ベクトル
scale = 8e4
ax.arrow(x[0]/1e3, y[0]/1e3, v_in_pc[0]*scale/1e3, v_in_pc[1]*scale/1e3,
         head_width=2e4, color='green', label='v∞ in')
ax.arrow(x[-1]/1e3, y[-1]/1e3, v_out_pc[0]*scale/1e3, v_out_pc[1]*scale/1e3,
         head_width=2e4, color='red', label='v∞ out')
ax.set_xlabel('x [10^3 km]')
ax.set_ylabel('y [10^3 km]')
ax.set_title('Planet-centered frame: hyperbolic trajectory')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
ax.legend(loc='upper right')

# 右: 速度空間でのΔvベクトル合成
ax = axes[1]
origin = np.array([0, 0])
ax.arrow(*origin, *v_in_pc,  head_width=0.4, color='green', label='v∞ in (planet frame)')
ax.arrow(*origin, *v_out_pc, head_width=0.4, color='red',   label='v∞ out (planet frame)')
ax.arrow(*v_in_pc, *(v_out_pc - v_in_pc), head_width=0.4, color='purple',
         label=f'Δv = {np.linalg.norm(dv_pc):.2f} km/s')
ax.set_xlabel('vx [km/s]')
ax.set_ylabel('vy [km/s]')
ax.set_title('Velocity-space: Δv vector triangle')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
ax.legend(loc='upper right')
ax.set_xlim(-3, 13)
ax.set_ylim(-7, 7)

plt.tight_layout()
plt.savefig('swingby_trajectory_and_dv.png', dpi=150, bbox_inches='tight')
plt.show()

print(f"|Δv| (numerical) = {np.linalg.norm(dv_pc):.4f} km/s")
print(f"|v_helio_in|     = {np.linalg.norm(v_in_helio):.4f} km/s")
print(f"|v_helio_out|    = {np.linalg.norm(v_out_helio):.4f} km/s")
print(f"太陽中心系の速さの変化: {np.linalg.norm(v_helio_out) - np.linalg.norm(v_helio_in):.4f} km/s"
      if False else "")
print(f"|v_out|^2 - |v_in|^2 = {np.linalg.norm(v_out_helio)**2 - np.linalg.norm(v_in_helio)**2:.3f} km^2/s^2")

左図には惑星中心系での双曲線軌道がはっきり描かれ、入射ベクトル(緑)と出射ベクトル(赤)の大きさが等しく向きだけ変わっている様子が見て取れます。右図の速度空間では、$\bm{v}_\infty^\text{in}$ と $\bm{v}_\infty^\text{out}$ が同じ長さ $v_\infty$ のベクトルで、両者を結ぶ紫の矢印が $\Delta\bm{v}$ です。この紫の矢印の長さがそのまま太陽中心系でのΔvの大きさになります。理論式 $|\Delta\bm{v}| = 2 v_\infty \sin(\delta/2)$ の値と一致することを数値で確認できます。

ボイジャー風グランドツアーの簡易シミュレーション

最後に、ボイジャー2号風に「木星 → 土星」と連続スイングバイで太陽中心系のエネルギーを段階的に増やす2次元シミュレーションを作ります。簡単のため、すべての惑星が同じ平面の円軌道を公転していると仮定します。

import numpy as np
import matplotlib.pyplot as plt

# 太陽の重力定数
mu_sun = 1.32712e11  # [km^3 / s^2]
AU = 1.49598e8       # [km]

# 各惑星の公転半径(簡易)と惑星情報
PLANET_R = {'Earth': 1.0*AU, 'Jupiter': 5.20*AU, 'Saturn': 9.58*AU}

def heliocentric_speed_after_swingby(v_helio_in_vec, V_planet_vec, r_p, mu_planet):
    """太陽中心系での入射速度 v_helio_in_vec を、
    惑星 (公転速度 V_planet_vec, 重力定数 mu_planet, 近点 r_p)
    でスイングバイした後の太陽中心系速度ベクトルを返す。
    曲げる向きは「v∞ を V_planet と同じ方向に近づける」と仮定(加速ゲイン最大寄り)"""
    v_inf_in = v_helio_in_vec - V_planet_vec
    v_inf_mag = np.linalg.norm(v_inf_in)
    e = 1 + r_p * v_inf_mag**2 / mu_planet
    delta = 2*np.arcsin(1/e)
    # v∞_in を delta だけ回転(V_planet 方向に向かって曲げる)
    sign = 1 if np.cross(v_inf_in, V_planet_vec) > 0 else -1
    c, s = np.cos(sign * delta), np.sin(sign * delta)
    R = np.array([[c, -s], [s, c]])
    v_inf_out = R @ v_inf_in
    return V_planet_vec + v_inf_out, delta

# 出発時:地球の公転速度 + 余剰速度 v∞_Earth で外向きに射出
V_Earth = np.array([0.0, np.sqrt(mu_sun/PLANET_R['Earth'])])
v_helio_0 = V_Earth + np.array([0.0, 8.7])  # ホーマン的に外側に向かう想定
print(f"打上げ直後の太陽中心系速度: {np.linalg.norm(v_helio_0):.3f} km/s")

# 木星スイングバイ(r_p = 4 R_J を仮定)
V_Jupiter = np.array([-np.sqrt(mu_sun/PLANET_R['Jupiter']), 0.0])
v_helio_1, delta_J = heliocentric_speed_after_swingby(
    v_helio_0, V_Jupiter,
    r_p=4*PLANETS['Jupiter']['R'], mu_planet=PLANETS['Jupiter']['mu'])
print(f"木星後 δ={np.degrees(delta_J):.1f} deg, |v_helio|={np.linalg.norm(v_helio_1):.3f} km/s")

# 土星スイングバイ(r_p = 5 R_S を仮定)
V_Saturn = np.array([0.0, -np.sqrt(mu_sun/PLANET_R['Saturn'])])
v_helio_2, delta_S = heliocentric_speed_after_swingby(
    v_helio_1, V_Saturn,
    r_p=5*PLANETS['Saturn']['R'], mu_planet=PLANETS['Saturn']['mu'])
print(f"土星後 δ={np.degrees(delta_S):.1f} deg, |v_helio|={np.linalg.norm(v_helio_2):.3f} km/s")

# 軌道エネルギーで太陽からの脱出可能性を判定
def specific_energy(v, r): return v**2/2 - mu_sun/r
E0 = specific_energy(np.linalg.norm(v_helio_0), PLANET_R['Earth'])
E1 = specific_energy(np.linalg.norm(v_helio_1), PLANET_R['Jupiter'])
E2 = specific_energy(np.linalg.norm(v_helio_2), PLANET_R['Saturn'])
print(f"\n太陽中心系での比エネルギー [km^2/s^2]:")
print(f"  打上げ直後: {E0:+.3f}")
print(f"  木星後   : {E1:+.3f}")
print(f"  土星後   : {E2:+.3f}  (正なら太陽から脱出)")

この簡易シミュレーションでは、地球から $v_\infty \approx 8.7$ km/s で外向きに打ち上げた宇宙機が、最初は太陽に束縛された楕円軌道(負のエネルギー)にあります。木星スイングバイで太陽中心系の速度が大幅に増し、軌道エネルギーが正に転じることがあります(パラメータ次第)。さらに土星スイングバイで追加のエネルギーを得て、ボイジャー1号・2号が実際に達成した「太陽系脱出」と同じ振る舞いを再現できます。出力の比エネルギーが負から正に転じる瞬間が、まさに「太陽の重力から逃れた」瞬間です。シンプルな2Dモデルでも、グラビティアシスト連鎖の威力がはっきり見えます。

ここまでで理論と実装が整いました。最後に、現実のミッションがどのように グラビティアシストを使ってきたか、いくつかの代表例を見ていきます。

応用 — ボイジャー・はやぶさ・パーカー

ボイジャー2号のグランドツアー

1977年8月に打ち上げられたボイジャー2号は、175年に一度しか巡ってこない木星・土星・天王星・海王星の特殊な配置を利用しました。木星到着が1979年7月、土星が1981年8月、天王星が1986年1月、そして海王星が1989年8月。各惑星でスイングバイを行い、太陽中心系のエネルギーを段階的に増やしながら、次の惑星のフライバイ条件を満たす軌道に乗せ替えていきました。

特に木星スイングバイで得たΔvは推定で約 $10$ km/s に達し、これがなければ土星以遠への航行は化学推進では不可能でした。ボイジャーが現在もなお秒速 $15$ km 以上で太陽系外に向かって飛んでいるのは、惑星4個ぶんの「お裾分け」の積み重ねの結果です。

Cassini の VVEJGA

土星探査機 Cassini(1997年打ち上げ、2004年到着)は、Venus-Venus-Earth-Jupiter Gravity Assist(VVEJGA)と呼ばれる4回のスイングバイを行いました。直接ホーマン遷移で土星に行くには $\Delta v \approx 16$ km/s が必要で、これは化学推進では実現困難。そこでまず金星方向に内向きに打ち上げ、金星で2回スイングバイして加速、その勢いで地球に戻ってさらにスイングバイし、最後に木星スイングバイで土星へ向かう、という7年がかりのルートが組まれました。打ち上げ時の必要$\Delta v$は約 $4$ km/s まで圧縮され、現実的なロケットで実現できました。

はやぶさの地球スイングバイ

JAXAの「はやぶさ」(2003年打ち上げ)は、イオンエンジンによる低推力推進と組み合わせた地球スイングバイを2004年5月に実施しました。地球を $\sim 3{,}725$ km まで近づき、太陽中心系での速度を約 $3.8$ km/s 加速。これによって小惑星イトカワとの軌道遷移が可能になりました。低推力推進と高効率スイングバイの組み合わせは、はやぶさ2、MMX(火星衛星探査計画)など、その後のJAXAの深宇宙ミッションの基本戦略となっています。

パーカー・ソーラー・プローブ(PSP)の逆利用

2018年打ち上げの Parker Solar Probe は、太陽の超近接観測(最終目標:太陽半径の $9.86$ 倍まで接近)のために、金星で7回スイングバイしてエネルギーを「捨てる」設計です。前節で述べたとおり、太陽に近づくには太陽中心系で減速する必要があり、これを金星スイングバイで実現しています。各スイングバイで近日点を段階的に縮め、最終的には $0.046$ au(約 $6.9 \times 10^6$ km)まで太陽に接近します。秒速 $190$ km の太陽周回速度はあらゆる人工天体で最速の記録です。

JUICE / Europa Clipper のマルチフライバイツアー

ESAのJUICE(2023年打ち上げ、2031年到着)はガリレオ衛星探査のために、地球-金星-地球-地球-地球と複雑な4回スイングバイを実施。木星到着後もガニメデ・カリスト・エウロパで30回以上のフライバイを繰り返し、軌道を細かく調整します。NASAのEuropa Clipper(2024年打ち上げ、2030年到着)も同様で、火星-地球とスイングバイを連ねたあと、エウロパで多数のフライバイを行う計画です。これらのミッションは、惑星間ルート設計が「Lambert問題(指定2点間の軌道)」と「グラビティアシストの幾何学」の高度な組み合わせ最適化であることを示す好例です。

Lambert問題との接続

惑星間航行の設計は、本質的に「ある日に惑星Aを出発し、別の日に惑星Bに到着する軌道」を解く Lambert問題と、各惑星でのスイングバイの偏角制約を組み合わせる最適化問題です。最適化の自由度はスイングバイの順序、出発・通過・到着日、各スイングバイでの $r_p$ など多数。実用の設計では、これらを遺伝的アルゴリズムや微分進化などのメタヒューリスティクスで探索し、Pork-chop plot(出発日と到着日に対する $\Delta v$ の等高線図)で全体像を把握する手法が広く用いられています。グラビティアシストの偏角公式 $\sin(\delta/2) = 1/e$ は、そうした探索の各ステップで「実行可能かどうか」を判定する基本不等式として機能します。

まとめ

本記事では、グラビティアシストの理論的基盤から実装、ミッション応用までを解説しました。重要なポイントを振り返ります。

  • 基準系の取り換えが本質 — 惑星中心系では宇宙機の速さは保存(双曲線軌道のエネルギー保存)、向きだけが偏角 $\delta$ だけ変わる。太陽中心系では惑星の公転速度ベクトル $\bm{V}_\text{pl}$ が合成され、速度の大きさも変わる。
  • 偏角公式 — $\sin(\delta/2) = 1/e$、$e = 1 + r_p v_\infty^2 / \mu$。$r_p$ を小さくするほど大きく曲げられる。
  • 太陽中心系のΔv — $|\Delta\bm{v}_\odot| = 2 v_\infty \sin(\delta/2) = 2 v_\infty / e$。理論上限は $2 v_\infty$($\delta = \pi$、つまり真後ろに跳ね返る)。
  • エネルギー保存 — 系全体のエネルギーは厳密に保存し、宇宙機が得たΔvぶんだけ惑星の公転速度が(質量比のぶんだけ)僅かに減る。
  • 最適 $v_\infty$ — 近点距離 $r_p$ を固定すると $v_\infty^* = \sqrt{\mu/r_p}$ で最大Δvが得られる。実際のミッションは $v_\infty < v_\infty^*$ の領域で運用される。
  • 応用 — ボイジャーのグランドツアー、Cassiniの VVEJGA、はやぶさの地球スイングバイ、PSPの金星×7減速、JUICE/Europa Clipperのマルチフライバイ。いずれも化学推進だけでは到達不可能なミッションをグラビティアシスト連鎖で実現している。

グラビティアシストは「重力という風に乗る」セーリング技術です。一見奇跡的にエネルギーが増えて見える背景には、ガリレイ変換と二体問題の保存則という極めて古典的な物理学がある——これがこの技術の美しいところです。

次のステップとして、以下の記事も参考にしてください。