ビザビバ方程式(vis-viva equation)を導出して軌道計算に活用する

国際宇宙ステーション(ISS)は、地上約420 kmの高さを秒速およそ7.7 kmで周回しています。これは音速の約22倍、東京から大阪まで約1分で到達する速さです。では、なぜISSはこの速度で飛んでいるのでしょうか? 速すぎれば地球から離れてしまい、遅すぎれば落下してしまいます。「ちょうどよい速度」は軌道の高さと形によって一意に決まるのです。

この「軌道上の任意の位置での速度」を一発で求めてくれるのが、ビザビバ方程式(vis-viva equation)です。ラテン語で「生きた力(vis viva)」を意味するこの名前は、17世紀にライプニッツが提唱した運動エネルギーの概念に由来しています。ライプニッツは $mv^2$ こそが運動の本質的な量だと主張しましたが、この議論は後にエネルギー保存則として結実しました。ビザビバ方程式はまさにその保存則から生まれた公式です。

ビザビバ方程式を理解すると、次のような幅広い問題に答えることができます。

  • 人工衛星の軌道設計: 衛星を特定の高度に投入するために必要な速度を計算できます。静止軌道(高度約36,000 km)の衛星がどのくらいの速さで周回しているかもすぐにわかります
  • 惑星間遷移軌道の設計: 地球から月へ、あるいは火星へ向かうときに必要な速度増分($\Delta v$)を見積もれます。ホーマン遷移軌道の計算はビザビバ方程式なしでは成り立ちません
  • 脱出速度の計算: 地球の重力圏を抜け出すために最低限必要な速度(第二宇宙速度)が、この方程式の特殊ケースとして自然に得られます
  • スイングバイの解析: 惑星の重力を利用して加速する技術の理解にも、双曲線軌道上でのビザビバ方程式が活躍します

本記事の内容

  • 軌道上で「速度が決まる仕組み」の直感的理解
  • エネルギー保存則と角運動量保存を用いたビザビバ方程式の完全な導出
  • 円軌道・楕円軌道・放物線軌道・双曲線軌道への適用と物理的意味
  • ISS・月遷移・火星遷移・静止軌道など実践的な具体例
  • Pythonによる軌道速度の可視化と数値計算

前提知識

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

軌道における速度の直感

ジェットコースターに乗ったことがあるでしょうか。高い場所ではゆっくりと進み、低い場所では猛スピードで駆け抜けます。これはエネルギー保存則の現れです。高い場所では位置エネルギーが大きく運動エネルギーが小さい。低い場所ではその逆です。両者の合計は常に一定に保たれます。

人工衛星の軌道でもまったく同じことが起きています。地球に近い場所(近地点)では衛星は速く動き、地球から遠い場所(遠地点)では遅く動きます。ジェットコースターにおける「高さ」が、軌道における「地球からの距離」に対応するのです。

ただし、ジェットコースターと衛星には大きな違いが一つあります。ジェットコースターの位置エネルギーは $mgh$ と高さに比例しますが、衛星の万有引力ポテンシャルは $-GMm/r$ と距離の逆数に比例します。この違いが、軌道力学を地上の力学よりも一段奥深くしている理由です。

とはいえ、エネルギー保存則という根本原理は共通しています。「運動エネルギー + ポテンシャルエネルギー = 一定」という関係を正しく書き下せば、軌道上のあらゆる位置での速度を求められるはずです。

では、この直感を数式で表現していきましょう。まず、衛星がもつ力学的エネルギーを正確に定義することから始めます。

軌道エネルギーの定義

中心天体(質量 $M$、たとえば地球)のまわりを運動する衛星(質量 $m$、ただし $m \ll M$)を考えます。二体問題の枠組みでは、中心天体は動かないと近似できるため、衛星の運動だけに注目すればよいことになります。

衛星は2種類のエネルギーをもっています。まず、速度 $v$ で運動することによる運動エネルギーです。

$$ K = \frac{1}{2}mv^2 $$

これは常に正(またはゼロ)で、衛星が速く動くほど大きくなります。

次に、万有引力によるポテンシャルエネルギーです。

$$ U = -\frac{GMm}{r} $$

ここで $G$ は万有引力定数($6.674 \times 10^{-11}$ m$^3$/(kg$\cdot$s$^2$))、$r$ は中心天体の中心から衛星までの距離です。ポテンシャルエネルギーが負の値をとるのは、「無限遠を基準($U = 0$)として、重力に引かれて落ちてきた分だけエネルギーが下がっている」ことを意味しています。地球に近づくほど $r$ が小さくなるので $U$ は大きな負の値になり、衛星はより深い重力の井戸にはまっている状態です。

万有引力は保存力であるため、力学的エネルギーの和は時間によらず一定に保たれます。

$$ \begin{equation} E = \frac{1}{2}mv^2 – \frac{GMm}{r} = \text{一定} \end{equation} $$

この式がビザビバ方程式の出発点です。衛星が軌道上のどの位置にいても、速度 $v$ と距離 $r$ の組み合わせは変化しますが、その「差し引きの合計」は変わりません。近地点で $r$ が小さく $U$ が深くなれば、そのぶん $v$ が大きくなって $K$ が増え、合計を一定に保つのです。

ここで、軌道力学で頻繁に登場する記号を導入しておきます。標準重力パラメータ(gravitational parameter)を次のように定義します。

$$ \mu = GM $$

地球の場合 $\mu_{\oplus} = 3.986 \times 10^{14}$ m$^3$/s$^2$ です。$G$ と $M$ は個別に精度よく測定することが難しいのですが、$\mu$ の値は天体観測から極めて高精度に求められます。そのため、軌道力学では $G$ と $M$ を分けずに $\mu$ をひとまとめに使うのが標準的です。

また、質量 $m$ で両辺を割った比力学的エネルギー(specific mechanical energy)を $\varepsilon$ と書くと、

$$ \varepsilon = \frac{E}{m} = \frac{1}{2}v^2 – \frac{\mu}{r} = \text{一定} $$

となります。比エネルギーは衛星の質量に依存しないため、「軌道そのもの」の性質を表す量です。同じ軌道に乗っている衛星なら、質量が1 kgでも1,000 kgでも比エネルギーは同じです。

ここまでで「軌道上のエネルギーは保存される」ことがわかりました。しかし、このままでは「一定」の値がいくらなのかがわかりません。次に、楕円軌道の幾何学的性質を使って、この定数を軌道パラメータ(半長軸 $a$)だけで表現します。

楕円軌道の全エネルギー — 半長軸で決まる

ケプラーの第一法則によると、中心力場での軌道は楕円(または円・放物線・双曲線)になります。楕円軌道の形は半長軸 $a$(楕円の長い方の半径)と離心率 $e$(楕円のつぶれ具合、$0 \leq e < 1$)で決まります。

ここで驚くべき事実があります。楕円軌道の全エネルギーは半長軸 $a$ だけで決まり、離心率 $e$ には依存しないのです。これは直感的には、楕円が細長かろうが丸かろうが、半長軸が同じならエネルギーは同じ、ということを意味します。

結論を先に述べると、比エネルギーは次のように表されます。

$$ \begin{equation} \varepsilon = -\frac{\mu}{2a} \end{equation} $$

$a > 0$(楕円軌道)のとき $\varepsilon < 0$ です。これは衛星が中心天体の重力に束縛されていることを意味します。$\varepsilon$ がちょうどゼロなら放物線軌道(束縛されていないが余裕もない状態)、正なら双曲線軌道(重力圏を脱出できる状態)です。

では、この重要な結果を証明しましょう。

証明: 近点と遠点での保存則を連立する

楕円軌道には2つの特別な点があります。中心天体に最も近い点(近点、$r = r_p$)と、最も遠い点(遠点、$r = r_a$)です。この2点では衛星の速度ベクトルが動径方向に垂直になるため、計算が大幅に簡単になります。

証明のゴールは「$\varepsilon = -\mu/(2a)$」を示すことです。使う道具はエネルギー保存則と角運動量保存則の2つです。

ステップ1: エネルギー保存則を近点と遠点で等置する

近点での速度を $v_p$、遠点での速度を $v_a$ として、比エネルギーが一定であることから、

$$ \frac{1}{2}v_p^2 – \frac{\mu}{r_p} = \frac{1}{2}v_a^2 – \frac{\mu}{r_a} $$

ステップ2: 角運動量保存則を立てる

中心力場では角運動量が保存されます。近点と遠点では速度が動径に垂直なので、比角運動量 $h$ は特に簡単な形をとります。

$$ r_p v_p = r_a v_a = h $$

これは「地球に近いほど速く、遠いほど遅い」というケプラーの第二法則(面積速度一定)の数学的表現にほかなりません。

ステップ3: 角運動量を使って速度を消去する

$v_p = h/r_p$、$v_a = h/r_a$ をエネルギー保存の式に代入します。

$$ \frac{1}{2}\frac{h^2}{r_p^2} – \frac{\mu}{r_p} = \frac{1}{2}\frac{h^2}{r_a^2} – \frac{\mu}{r_a} $$

ポテンシャルの項を右辺に、運動エネルギーの項を左辺に集めます。

$$ \frac{h^2}{2}\left(\frac{1}{r_p^2} – \frac{1}{r_a^2}\right) = \mu\left(\frac{1}{r_p} – \frac{1}{r_a}\right) $$

ステップ4: 左辺を因数分解する

左辺の括弧内は「二乗の差」なので、$(a^2 – b^2) = (a+b)(a-b)$ の形に因数分解できます。

$$ \frac{h^2}{2} \cdot \frac{(r_a + r_p)(r_a – r_p)}{r_p^2 r_a^2} = \mu \cdot \frac{r_a – r_p}{r_p r_a} $$

楕円軌道では $r_a \neq r_p$(近点と遠点は異なる位置)なので、両辺を共通因子 $(r_a – r_p)$ で割ることができます。

$$ \frac{h^2}{2} \cdot \frac{r_a + r_p}{r_p r_a} = \mu $$

ステップ5: 楕円の幾何学的関係を使って $h^2$ を求める

楕円の定義から $r_p + r_a = 2a$ が成り立ちます(軌道要素の定義を参照)。これを代入して $h^2$ について解きます。

$$ h^2 = \frac{2\mu r_p r_a}{r_p + r_a} = \frac{\mu r_p r_a}{a} $$

これで角運動量が軌道パラメータだけで表されました。

ステップ6: 比エネルギーを近点で計算する

得られた $h^2$ を使って、近点での比エネルギーを計算します。

$$ \varepsilon = \frac{1}{2}v_p^2 – \frac{\mu}{r_p} = \frac{h^2}{2r_p^2} – \frac{\mu}{r_p} $$

$h^2 = \mu r_p r_a / a$ を代入すると、

$$ \varepsilon = \frac{\mu r_a}{2 a r_p} – \frac{\mu}{r_p} $$

$\mu/r_p$ でくくります。

$$ \varepsilon = \frac{\mu}{r_p}\left(\frac{r_a}{2a} – 1\right) $$

ここで $2a = r_p + r_a$ なので、$r_a/(2a) = r_a/(r_p + r_a)$ です。これを代入すると、

$$ \varepsilon = \frac{\mu}{r_p} \cdot \frac{r_a – (r_p + r_a)}{r_p + r_a} = \frac{\mu}{r_p} \cdot \frac{-r_p}{2a} = -\frac{\mu}{2a} $$

以上により、

$$ \begin{equation} \boxed{\varepsilon = -\frac{\mu}{2a}} \end{equation} $$

が示されました。全エネルギーは $E = m\varepsilon = -GMm/(2a)$ です。

この結果の物理的意味を確認しておきましょう。半長軸 $a$ が大きいほど(軌道が大きいほど)$|\varepsilon|$ は小さくなります。つまり、大きな軌道ほどエネルギーのマイナス幅が小さく、「あと少しで脱出できる」状態に近づくのです。逆に $a$ が小さいほど衛星は深い重力の井戸に沈んでおり、脱出が困難です。

また、離心率 $e$ がこの式に登場しないことは重要です。$a = 10{,}000$ km の円軌道と、$a = 10{,}000$ km で $e = 0.5$ の楕円軌道は、全エネルギーがまったく同じです。違うのはエネルギーの「配分」(位置エネルギーと運動エネルギーの比率)だけであり、総量は半長軸のみで決まります。

これで「一定」の正体がわかりました。あとはエネルギー保存の式に代入して速度を求めれば、ビザビバ方程式の完成です。

ビザビバ方程式の導出

いよいよ本題です。ここまでに得られた2つの式を組み合わせます。

出発点1: 比エネルギーの定義

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

出発点2: 楕円軌道の比エネルギー

$$ \varepsilon = -\frac{\mu}{2a} $$

これらを等置すると、

$$ \frac{1}{2}v^2 – \frac{\mu}{r} = -\frac{\mu}{2a} $$

$v^2$ について解きます。まず $\mu/r$ を右辺に移項します。

$$ \frac{1}{2}v^2 = \frac{\mu}{r} – \frac{\mu}{2a} $$

右辺の $\mu$ をくくり出します。

$$ \frac{1}{2}v^2 = \mu\left(\frac{1}{r} – \frac{1}{2a}\right) = \mu\left(\frac{2a – r}{2ar}\right) $$

両辺に2を掛けると、

$$ v^2 = \mu \cdot \frac{2(2a – r)}{2ar} = \mu\left(\frac{2}{r} – \frac{1}{a}\right) $$

こうして得られたのがビザビバ方程式です。

$$ \begin{equation} \boxed{v^2 = \mu\left(\frac{2}{r} – \frac{1}{a}\right)} \end{equation} $$

この方程式は極めてシンプルでありながら、驚くほど強力です。軌道上の任意の位置 $r$ での速度 $v$ を、標準重力パラメータ $\mu$ と半長軸 $a$ だけから即座に計算できます。離心率も真近点角も軌道傾斜角も不要です。必要なのは「今どこにいるか」($r$)と「軌道の大きさ」($a$)だけなのです。

式の構造を眺めると、$2/r$ の項は「現在位置の重力ポテンシャルから得られる速度の寄与」、$-1/a$ の項は「軌道全体のエネルギーによる制限」と解釈できます。近い位置($r$ が小さい)では $2/r$ が大きくなるため速度が上がり、遠い位置では速度が下がります。これはまさに、冒頭のジェットコースターのアナロジーと合致しています。

それでは、この一般公式がさまざまな軌道形状でどのような特殊ケースを生むのか、一つずつ見ていきましょう。

特殊な軌道への適用

ビザビバ方程式 $v^2 = \mu(2/r – 1/a)$ は、半長軸 $a$ の値に応じて異なる軌道形状を表します。ここでは円軌道、放物線軌道、双曲線軌道の3つの重要な特殊ケースを取り上げ、それぞれの物理的意味を詳しく解説します。

円軌道 — 高さが決まれば速度が決まる($a = r$)

円軌道は離心率 $e = 0$ の楕円軌道であり、衛星は常に同じ距離 $r$ を保って中心天体を周回します。このとき半長軸 $a$ は軌道半径 $r$ に等しいので、ビザビバ方程式は、

$$ v_{\text{circ}}^2 = \mu\left(\frac{2}{r} – \frac{1}{r}\right) = \frac{\mu}{r} $$

したがって、

$$ \begin{equation} v_{\text{circ}} = \sqrt{\frac{\mu}{r}} \end{equation} $$

この式は「円軌道速度」と呼ばれ、第一宇宙速度の基本式でもあります。高度だけで速度が一意に決まるという、非常にすっきりした結果です。

地球表面($r = R_{\oplus} \approx 6{,}371$ km)に適用すると、

$$ v_1 = \sqrt{\frac{3.986 \times 10^{14}}{6.371 \times 10^6}} \approx 7.91 \text{ km/s} $$

これが地表ギリギリ(大気の抵抗を無視した場合)で周回するのに必要な速度です。実際の衛星は大気の影響を避けるためにもっと高い軌道を使いますが、高度が上がるにつれて $r$ が大きくなるため必要速度は減少します。

円軌道速度が $r$ の平方根に反比例する($v \propto 1/\sqrt{r}$)ことは重要な性質です。高度35,786 kmの静止軌道では $r \approx 42{,}157$ km となり、速度は約3.07 km/sまで下がります。地表近くの7.9 km/sと比べると半分以下です。「高い軌道ほどゆっくり回る」のです。

放物線軌道 — 脱出のギリギリ($a \to \infty$)

軌道のエネルギーがちょうどゼロ($\varepsilon = 0$)のとき、半長軸 $a \to \infty$ となります。このとき軌道は放物線になり、衛星は中心天体からの脱出にギリギリ成功する状態です。

ビザビバ方程式で $1/a \to 0$ とすると、

$$ v_{\text{esc}}^2 = \mu \cdot \frac{2}{r} $$

$$ \begin{equation} v_{\text{esc}} = \sqrt{\frac{2\mu}{r}} \end{equation} $$

これが脱出速度(escape velocity)、すなわち第二宇宙速度です。地球表面では、

$$ v_2 = \sqrt{\frac{2 \times 3.986 \times 10^{14}}{6.371 \times 10^6}} \approx 11.19 \text{ km/s} $$

円軌道速度との関係を見ると、

$$ v_{\text{esc}} = \sqrt{\frac{2\mu}{r}} = \sqrt{2} \cdot \sqrt{\frac{\mu}{r}} = \sqrt{2} \cdot v_{\text{circ}} $$

脱出速度は常に円軌道速度のちょうど $\sqrt{2} \approx 1.414$ 倍です。この関係は高度によらず成り立つ普遍的な性質であり、覚えておくと見積もりに便利です。

物理的な意味を考えましょう。放物線軌道では、衛星は無限遠($r \to \infty$)に到達したときにちょうど速度がゼロになります。エネルギーが余らず不足もしない、ぴったりの脱出です。現実のミッションでは「ぴったり」では効率が悪いため、少し余裕をもたせた双曲線軌道を使うのが一般的です。

双曲線軌道 — 余裕をもった脱出($a < 0$)

エネルギーが正($\varepsilon > 0$)のとき、$a = -\mu/(2\varepsilon)$ は負の値をとります。軌道は双曲線となり、衛星は中心天体のそばを通過した後、再び戻ってくることなく飛び去ります。

ビザビバ方程式に $a < 0$ を代入すると、

$$ v^2 = \mu\left(\frac{2}{r} + \frac{1}{|a|}\right) $$

$1/|a|$ の項が加わるため、同じ距離 $r$ で比較すると放物線軌道よりも速くなります。

特に、$r \to \infty$ としたときの速度(超過速度 $v_\infty$)は、

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

$$ \begin{equation} v_\infty = \sqrt{\frac{\mu}{|a|}} \end{equation} $$

これは「無限遠に到達しても残っている速度」です。放物線軌道では $v_\infty = 0$ だったのに対し、双曲線軌道では有限の速度が残ります。この超過速度 $v_\infty$ は惑星間航行のミッション設計で極めて重要なパラメータです。たとえば地球の重力圏を脱出した後の太陽に対する速度を決定します。

双曲線軌道の具体例としては、太陽系外に向かうボイジャー探査機や、他の惑星からの重力アシスト(スイングバイ)後の軌道があります。スイングバイでは、探査機が惑星に対して双曲線軌道で接近・通過し、$v_\infty$ の方向を変えることで太陽系内での軌道を効率よく変更します。

軌道タイプの判別まとめ

ここまでの結果を表にまとめます。軌道のタイプは比エネルギー $\varepsilon$ の符号で判別できます。

軌道タイプ 離心率 $e$ 半長軸 $a$ 比エネルギー $\varepsilon$ $v$ と $v_{\text{circ}}$ の関係
円軌道 $e = 0$ $a = r > 0$ $\varepsilon < 0$ $v = v_{\text{circ}}$
楕円軌道 $0 < e < 1$ $a > 0$ $\varepsilon < 0$ $v_{\text{circ}} < v < v_{\text{esc}}$
放物線軌道 $e = 1$ $a \to \infty$ $\varepsilon = 0$ $v = v_{\text{esc}}$
双曲線軌道 $e > 1$ $a < 0$ $\varepsilon > 0$ $v > v_{\text{esc}}$

この表は万有引力とケプラーの法則の内容とも深く関連しています。軌道のエネルギーが負なら束縛軌道(円・楕円)、ゼロ以上なら非束縛軌道(放物線・双曲線)です。

理論的な分類がわかったところで、次は具体的な数値を入れてビザビバ方程式を「使って」みましょう。実際のミッションでどのような速度が得られるかを確認します。

具体例: 実際の軌道での速度計算

例1: ISS(国際宇宙ステーション)

ISSは地球の低軌道(LEO: Low Earth Orbit)を周回しています。

  • 平均軌道高度: $h \approx 420$ km
  • 軌道半径: $r = R_{\oplus} + h = 6{,}371 + 420 = 6{,}791$ km $= 6.791 \times 10^6$ m
  • 離心率は非常に小さいため、ほぼ円軌道として扱えます

円軌道速度を計算します。

$$ v_{\text{ISS}} = \sqrt{\frac{\mu}{r}} = \sqrt{\frac{3.986 \times 10^{14}}{6.791 \times 10^6}} = \sqrt{5.870 \times 10^7} \approx 7{,}660 \text{ m/s} \approx 7.66 \text{ km/s} $$

ISSは約7.66 km/sで飛行しています。これは約90分で地球を1周する速さです。地上から見ると夜空を横切る明るい点として数分間見ることができますが、その間にISSは数百kmも移動しているのです。

この高度での脱出速度も計算しておきましょう。

$$ v_{\text{esc}} = \sqrt{2} \times 7.66 \approx 10.83 \text{ km/s} $$

ISSの速度は脱出速度よりも十分に小さいので、ISSは地球に束縛された円軌道を安定して維持しています。

例2: 静止軌道衛星(GEO)

静止軌道は、衛星の周期が地球の自転周期(約23時間56分)と一致する特別な円軌道です。気象衛星や通信衛星がこの軌道に配置されています。

  • 軌道高度: $h \approx 35{,}786$ km
  • 軌道半径: $r = 6{,}371 + 35{,}786 = 42{,}157$ km $= 4.216 \times 10^7$ m

$$ v_{\text{GEO}} = \sqrt{\frac{3.986 \times 10^{14}}{4.216 \times 10^7}} \approx 3{,}075 \text{ m/s} \approx 3.07 \text{ km/s} $$

ISSと比較すると約40%の速度しかありません。高い軌道では重力が弱いため、釣り合いに必要な速度も小さくなります。

例3: 月遷移軌道(TLI)

地球から月へ向かうミッション(アポロ計画など)では、地球周回の低軌道から月に到達する楕円軌道(遷移軌道)に乗り移ります。

  • 近地点: $r_p \approx 6{,}571$ km(地表から約200 kmの低軌道)
  • 遠地点: $r_a \approx 384{,}400$ km(月軌道の距離)
  • 半長軸: $a = (r_p + r_a)/2 = (6{,}571 + 384{,}400)/2 \approx 195{,}486$ km $= 1.955 \times 10^8$ m

近地点(出発点)での速度をビザビバ方程式で計算します。

$$ v_p = \sqrt{\mu\left(\frac{2}{r_p} – \frac{1}{a}\right)} = \sqrt{3.986 \times 10^{14}\left(\frac{2}{6.571 \times 10^6} – \frac{1}{1.955 \times 10^8}\right)} $$

$$ v_p = \sqrt{3.986 \times 10^{14} \times (3.044 \times 10^{-7} – 5.115 \times 10^{-9})} \approx \sqrt{3.986 \times 10^{14} \times 2.993 \times 10^{-7}} $$

$$ v_p \approx \sqrt{1.193 \times 10^{8}} \approx 10{,}920 \text{ m/s} \approx 10.92 \text{ km/s} $$

低軌道での円軌道速度は約7.73 km/s なので、月遷移に必要な速度増分は、

$$ \Delta v_{\text{TLI}} = v_p – v_{\text{circ}} \approx 10.92 – 7.73 \approx 3.19 \text{ km/s} $$

この $\Delta v$ がTLI(Trans-Lunar Injection)噴射で必要なエネルギーの指標です。アポロ計画ではサターンVロケットの第3段(S-IVB)がこの加速を担いました。

例4: 火星遷移軌道(ホーマン遷移)

地球から火星へ最もエネルギー効率のよい経路はホーマン遷移軌道です。太陽を中心とした楕円軌道で、近日点が地球軌道、遠日点が火星軌道に接する軌道を考えます。

  • 地球の太陽周回半径: $r_{\oplus} \approx 1.496 \times 10^{11}$ m(1 AU)
  • 火星の太陽周回半径: $r_{\text{Mars}} \approx 2.279 \times 10^{11}$ m(1.524 AU)
  • 太陽の標準重力パラメータ: $\mu_{\odot} = 1.327 \times 10^{20}$ m$^3$/s$^2$
  • 遷移軌道の半長軸: $a = (r_{\oplus} + r_{\text{Mars}})/2 = 1.888 \times 10^{11}$ m

出発点(地球軌道)での速度をビザビバ方程式で求めます。

$$ v_{\text{dep}} = \sqrt{\mu_{\odot}\left(\frac{2}{r_{\oplus}} – \frac{1}{a}\right)} = \sqrt{1.327 \times 10^{20}\left(\frac{2}{1.496 \times 10^{11}} – \frac{1}{1.888 \times 10^{11}}\right)} $$

$$ v_{\text{dep}} \approx \sqrt{1.327 \times 10^{20} \times 8.067 \times 10^{-12}} \approx \sqrt{1.071 \times 10^{9}} \approx 32{,}720 \text{ m/s} \approx 32.7 \text{ km/s} $$

地球の太陽周回速度は約29.8 km/sなので、地球に対する超過速度は $32.7 – 29.8 \approx 2.9$ km/s です。ただし、これは太陽中心系での値であり、実際のミッションでは地球の重力圏脱出のための $\Delta v$ も加わります。

これらの例から、ビザビバ方程式がいかに幅広い軌道計算に使えるかがわかります。次に、Pythonを使ってこれらの計算を体系的に可視化してみましょう。

Pythonでの実装: 軌道速度の可視化

ここでは3つの図を作成します。(1) 高度に対する円軌道速度と脱出速度の変化、(2) 楕円軌道上の速度変化、(3) 軌道形状の比較です。これにより、ビザビバ方程式が示す物理的な振る舞いを視覚的に確認します。

import numpy as np
import matplotlib.pyplot as plt

# 定数
mu = 3.986e14      # 地球の標準重力パラメータ [m³/s²]
R_E = 6.371e6      # 地球半径 [m]

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# --- 図1: 軌道速度 vs 高度 ---
h = np.linspace(200e3, 40000e3, 1000)   # 高度 200 km ~ 40,000 km
r = R_E + h

# 円軌道速度(vis-viva で a = r)
v_circ = np.sqrt(mu / r) / 1000  # km/s に変換
# 脱出速度(vis-viva で a → ∞)
v_esc = np.sqrt(2 * mu / r) / 1000

axes[0].plot(h / 1000, v_circ, 'b-', linewidth=2, label='Circular orbit $v_{circ}$')
axes[0].plot(h / 1000, v_esc, 'r--', linewidth=2, label='Escape $v_{esc}$')

# 代表的な衛星の位置をプロット
satellites = {
    'ISS (420 km)': 420,
    'GEO (35786 km)': 35786,
}
for name, h_sat in satellites.items():
    r_sat = R_E + h_sat * 1e3
    v_sat = np.sqrt(mu / r_sat) / 1000
    axes[0].plot(h_sat, v_sat, 'o', markersize=8, zorder=5)
    axes[0].annotate(name, (h_sat, v_sat), textcoords="offset points",
                     xytext=(10, 10), fontsize=9)

axes[0].set_xlabel('Altitude [km]', fontsize=12)
axes[0].set_ylabel('Velocity [km/s]', fontsize=12)
axes[0].set_title('Orbital velocity vs altitude', fontsize=14)
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim([0, 40000])

# --- 図2: 楕円軌道上の速度変化 ---
a = 15000e3    # 半長軸 15,000 km [m]
e = 0.5        # 離心率
theta = np.linspace(0, 2 * np.pi, 500)
# 軌道方程式 r = a(1-e²)/(1+e·cosθ)
r_orbit = a * (1 - e**2) / (1 + e * np.cos(theta))

# ビザビバ方程式で各位置の速度を計算
v_orbit = np.sqrt(mu * (2 / r_orbit - 1 / a)) / 1000

axes[1].plot(np.degrees(theta), v_orbit, 'b-', linewidth=2)
axes[1].axhline(y=np.sqrt(mu / a) / 1000, color='r', linestyle='--',
                alpha=0.7, label='Circular velocity at $r=a$')
axes[1].axvline(x=0, color='g', linestyle=':', alpha=0.5, label='Periapsis ($\\theta=0$)')
axes[1].axvline(x=180, color='orange', linestyle=':', alpha=0.5, label='Apoapsis ($\\theta=180°$)')
axes[1].set_xlabel('True anomaly $\\theta$ [deg]', fontsize=12)
axes[1].set_ylabel('Velocity [km/s]', fontsize=12)
axes[1].set_title(f'Velocity on elliptical orbit ($e={e}$, $a={a/1e3:.0f}$ km)', fontsize=14)
axes[1].legend(fontsize=9)
axes[1].grid(True, alpha=0.3)

# --- 図3: 軌道形状の比較 ---
theta_plot = np.linspace(-np.pi * 0.95, np.pi * 0.95, 1000)

orbits = [
    (15000e3, 0.0, 'Circle ($e=0$)', 'blue'),
    (15000e3, 0.3, 'Ellipse ($e=0.3$)', 'green'),
    (15000e3, 0.6, 'Ellipse ($e=0.6$)', 'orange'),
    (15000e3, 0.9, 'Ellipse ($e=0.9$)', 'red'),
]

for a_orb, e_orb, label, color in orbits:
    r_plot = a_orb * (1 - e_orb**2) / (1 + e_orb * np.cos(theta_plot))
    r_plot = np.where(r_plot > 0, r_plot, np.nan)
    x = r_plot * np.cos(theta_plot) / 1e6
    y = r_plot * np.sin(theta_plot) / 1e6
    axes[2].plot(x, y, linewidth=2, label=label, color=color)

# 中心天体(地球)を描画
circle = plt.Circle((0, 0), R_E / 1e6, color='deepskyblue', alpha=0.4, label='Earth')
axes[2].add_patch(circle)
axes[2].set_xlabel('x [$\\times 10^3$ km]', fontsize=12)
axes[2].set_ylabel('y [$\\times 10^3$ km]', fontsize=12)
axes[2].set_title('Orbit shapes (same $a$, different $e$)', fontsize=14)
axes[2].legend(fontsize=9, loc='upper right')
axes[2].set_aspect('equal')
axes[2].grid(True, alpha=0.3)
axes[2].set_xlim([-35, 35])
axes[2].set_ylim([-35, 35])

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

このコードで生成される3つのグラフから、ビザビバ方程式の物理的な意味を読み取ることができます。

図1(左: 速度 vs 高度) からは、2つの重要な特徴がわかります。第一に、高度が上がるほど円軌道速度・脱出速度ともに単調に減少します。これは「遠い軌道ほどゆっくり回る」という直感と一致しています。第二に、脱出速度(赤い破線)は常に円軌道速度(青い実線)の $\sqrt{2}$ 倍であり、2つの曲線は常に一定の比率を保っています。ISSは低軌道のため高い速度が必要ですが、GEOでは半分以下の速度で十分であることが視覚的にわかります。

図2(中央: 楕円軌道上の速度変化) では、近点($\theta = 0°$)で速度が最大、遠点($\theta = 180°$)で速度が最小になることが明確に見えます。離心率 $e = 0.5$ の場合、近点速度は遠点速度の約3倍にもなります。赤い破線は半長軸 $a$ に相当する距離での円軌道速度を示しており、楕円軌道の速度がこの値を上下に振動していることがわかります。

図3(右: 軌道形状の比較) では、半長軸 $a$ が同じでも離心率 $e$ が異なると軌道の形が大きく変わることが示されています。$e = 0$ の完全な円から $e = 0.9$ の極端に細長い楕円まで、すべて同じ半長軸をもち、したがってビザビバ方程式より同じ全エネルギーをもっています。形が変わってもエネルギーは変わらない — これが「全エネルギーは $a$ のみで決まる」ことの視覚的な確認です。

次に、さまざまな高度での具体的な速度値を数値で確認してみましょう。

Pythonでの実装: 各種軌道の速度一覧

以下のコードでは、代表的な軌道高度での円軌道速度・脱出速度・軌道周期を一覧表示し、さらに月遷移軌道と火星遷移軌道のパラメータも計算します。

import numpy as np

# 定数
mu_earth = 3.986e14     # 地球の標準重力パラメータ [m³/s²]
mu_sun = 1.327e20       # 太陽の標準重力パラメータ [m³/s²]
R_E = 6.371e6           # 地球半径 [m]

print("=" * 70)
print("  地球周回軌道の速度一覧(ビザビバ方程式による計算)")
print("=" * 70)
print(f"{'軌道':>12s} {'高度[km]':>10s} {'v_circ[km/s]':>14s} "
      f"{'v_esc[km/s]':>13s} {'周期':>12s}")
print("-" * 70)

orbits_earth = [
    ("LEO (低)", 200),
    ("ISS", 420),
    ("LEO (高)", 2000),
    ("MEO (GPS)", 20200),
    ("GEO (静止)", 35786),
]

for name, h_km in orbits_earth:
    r = R_E + h_km * 1e3
    v_c = np.sqrt(mu_earth / r) / 1e3        # 円軌道速度 [km/s]
    v_e = np.sqrt(2 * mu_earth / r) / 1e3     # 脱出速度 [km/s]
    T = 2 * np.pi * np.sqrt(r**3 / mu_earth)  # 周期 [s]
    T_h = T / 3600                             # 周期 [時間]
    if T_h < 24:
        period_str = f"{T_h:.2f} h"
    else:
        period_str = f"{T_h / 24:.2f} days"
    print(f"{name:>12s} {h_km:>10d} {v_c:>14.3f} {v_e:>13.3f} {period_str:>12s}")

print()
print("=" * 70)
print("  遷移軌道の計算(ビザビバ方程式)")
print("=" * 70)

# 月遷移軌道(地球中心)
r_p_moon = R_E + 200e3           # 近地点: 地表 + 200 km
r_a_moon = 384400e3              # 遠地点: 月軌道
a_moon = (r_p_moon + r_a_moon) / 2

v_p_moon = np.sqrt(mu_earth * (2 / r_p_moon - 1 / a_moon)) / 1e3
v_circ_leo = np.sqrt(mu_earth / r_p_moon) / 1e3
dv_tli = v_p_moon - v_circ_leo

print(f"\n--- 月遷移軌道 (TLI) ---")
print(f"  近地点高度:     200 km")
print(f"  遠地点距離:     384,400 km (月軌道)")
print(f"  半長軸:         {a_moon/1e3:.0f} km")
print(f"  近地点速度:     {v_p_moon:.3f} km/s")
print(f"  LEO円軌道速度:  {v_circ_leo:.3f} km/s")
print(f"  必要 Delta-v:   {dv_tli:.3f} km/s")

# 火星ホーマン遷移軌道(太陽中心)
r_earth = 1.496e11      # 地球の太陽周回半径 [m]
r_mars = 2.279e11       # 火星の太陽周回半径 [m]
a_mars = (r_earth + r_mars) / 2

v_dep = np.sqrt(mu_sun * (2 / r_earth - 1 / a_mars)) / 1e3
v_arr = np.sqrt(mu_sun * (2 / r_mars - 1 / a_mars)) / 1e3
v_earth_orbit = np.sqrt(mu_sun / r_earth) / 1e3
v_mars_orbit = np.sqrt(mu_sun / r_mars) / 1e3

T_transfer = np.pi * np.sqrt(a_mars**3 / mu_sun)  # 遷移時間(半周期)

print(f"\n--- 火星ホーマン遷移軌道 ---")
print(f"  地球軌道半径:     {r_earth/1e9:.3f} x 10^8 km  (1 AU)")
print(f"  火星軌道半径:     {r_mars/1e9:.3f} x 10^8 km  (1.524 AU)")
print(f"  遷移軌道半長軸:   {a_mars/1e9:.3f} x 10^8 km")
print(f"  出発速度(太陽系): {v_dep:.3f} km/s")
print(f"  地球公転速度:     {v_earth_orbit:.3f} km/s")
print(f"  出発時の超過速度: {v_dep - v_earth_orbit:.3f} km/s")
print(f"  到着速度(太陽系): {v_arr:.3f} km/s")
print(f"  火星公転速度:     {v_mars_orbit:.3f} km/s")
print(f"  到着時の超過速度: {v_mars_orbit - v_arr:.3f} km/s")
print(f"  遷移時間:         {T_transfer / (3600*24):.1f} days ({T_transfer / (3600*24*30):.1f} months)")

このコードの出力を確認しましょう。地球周回軌道の一覧表では、高度200 kmのLEOで約7.8 km/s、ISS軌道で約7.7 km/s、GEO(静止軌道)で約3.1 km/sと算出されます。高度が上がるにつれて速度が低下していく様子が数値で裏付けられます。周期もケプラーの第三法則と整合しており、GEOでは約24時間(正確には約23時間56分)になっていることが確認できます。

月遷移軌道の計算では、近地点での速度が約10.9 km/sと求まり、LEO円軌道速度(約7.8 km/s)との差から $\Delta v \approx 3.1$ km/s が必要であるとわかります。これはアポロ計画で実際に使われた値と良く一致しています。

火星遷移軌道の計算では、遷移時間が約259日(約8.6か月)と算出されます。出発時に地球公転速度を約2.9 km/s 上回る速度が必要であり、到着時には火星公転速度を約2.6 km/s 下回る速度になります。到着時の速度不足分は火星への減速(あるいは大気ブレーキ)で補う必要があります。

次に、ビザビバ方程式と密接に関連する2つの宇宙速度(第一・第二宇宙速度)について、改めて整理しておきましょう。

第一宇宙速度と第二宇宙速度

ビザビバ方程式の特殊ケースとして自然に現れる2つの速度は、宇宙開発の文脈で特別な名前をもっています。

第一宇宙速度 — 「落ち続ける」ための速度

地表ギリギリで地球を周回する円軌道速度が第一宇宙速度 $v_1$ です。

$$ v_1 = \sqrt{\frac{\mu}{R_{\oplus}}} \approx 7.91 \text{ km/s} $$

この速度は「水平に投げた物体が、落下しながらも地面のカーブに沿って周り続ける速度」と理解できます。ニュートンの思考実験で有名な「大砲の弾」のイメージです。山の上から水平に大砲を撃ったとき、初速が小さければ弾は放物線を描いて落下します。しかし、初速を十分に大きくすると、弾は落下しているにもかかわらず地球の丸みのために地面に到達せず、地球を一周してしまいます。その臨界速度が第一宇宙速度です。

「落ちているのに落ちない」というのは不思議に思えるかもしれませんが、ISSの宇宙飛行士が無重力状態にあるのはまさにこの現象です。ISSは自由落下しており、その中にいる人や物も同じ加速度で落下しているので、相対的には無重力に見えるのです。

第二宇宙速度 — 「二度と戻らない」ための速度

地表から打ち出した物体が地球の重力圏を脱出するための最小速度が第二宇宙速度(脱出速度)$v_2$ です。

$$ v_2 = \sqrt{\frac{2\mu}{R_{\oplus}}} = \sqrt{2} \cdot v_1 \approx 11.19 \text{ km/s} $$

第二宇宙速度は第一宇宙速度のちょうど $\sqrt{2}$ 倍です。この関係はビザビバ方程式から直ちに導かれます。円軌道速度 $v_1 = \sqrt{\mu/r}$ に対して脱出速度 $v_2 = \sqrt{2\mu/r}$ なので、比をとると $v_2/v_1 = \sqrt{2}$ です。

物理的な意味として、脱出速度で打ち出された物体は放物線軌道に乗り、無限遠で速度がちょうどゼロになります。つまり「ギリギリ戻ってこない」状態です。これより速い場合は双曲線軌道に乗り、無限遠でも有限の速度が残ります。

なお、第一・第二宇宙速度はあくまで「地球に対する速度」です。太陽系を脱出する(第三宇宙速度)にはさらに追加のエネルギーが必要で、地表基準で約16.7 km/sになります。これは地球の公転運動を利用するかどうかで変わるため、打ち出す方向にも依存する少し複雑な問題です。

これらの宇宙速度の概念が実際にどのように使われるかを、もう一つPythonコードで可視化してみましょう。

Pythonでの実装: エネルギーダイアグラムと宇宙速度

ビザビバ方程式の本質は「エネルギー保存則」です。ここでは、衛星の比エネルギーを距離の関数としてプロットし、各軌道タイプのエネルギーレベルを視覚的に比較します。

import numpy as np
import matplotlib.pyplot as plt

mu = 3.986e14
R_E = 6.371e6

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

# --- 図1: エネルギーダイアグラム ---
r = np.linspace(R_E, 60000e3, 1000)

# ポテンシャルエネルギー(比エネルギー)
U = -mu / r / 1e6   # MJ/kg に変換

# 各軌道の全エネルギーライン
a_leo = R_E + 400e3         # LEO (h=400km)
a_geo = R_E + 35786e3       # GEO
a_transfer = (a_leo + a_geo) / 2  # GTO (遷移軌道)

E_leo = -mu / (2 * a_leo) / 1e6
E_geo = -mu / (2 * a_geo) / 1e6
E_gto = -mu / (2 * a_transfer) / 1e6
E_escape = 0  # 放物線軌道

axes[0].plot(r / 1e6, U, 'k-', linewidth=2, label='$U(r) = -\\mu/r$')
axes[0].axhline(y=E_leo, color='blue', linestyle='--', alpha=0.7,
                label=f'LEO ($a$={a_leo/1e3:.0f} km)')
axes[0].axhline(y=E_gto, color='green', linestyle='--', alpha=0.7,
                label=f'GTO ($a$={a_transfer/1e3:.0f} km)')
axes[0].axhline(y=E_geo, color='orange', linestyle='--', alpha=0.7,
                label=f'GEO ($a$={a_geo/1e3:.0f} km)')
axes[0].axhline(y=E_escape, color='red', linestyle='--', alpha=0.7,
                label='Escape ($\\varepsilon=0$)')

# 運動エネルギー = E - U を塗りつぶしで表示(LEOの例)
K_leo = E_leo - U  # U は負なので E_leo - U > 0
K_leo_clipped = np.where(K_leo > 0, K_leo, 0)
axes[0].fill_between(r / 1e6, U, E_leo,
                     where=(E_leo > U.flatten()),
                     alpha=0.15, color='blue', label='Kinetic energy (LEO)')

axes[0].set_xlabel('Distance $r$ [$\\times 10^3$ km]', fontsize=12)
axes[0].set_ylabel('Specific energy [MJ/kg]', fontsize=12)
axes[0].set_title('Energy diagram for Earth orbits', fontsize=14)
axes[0].legend(fontsize=9, loc='lower right')
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim([0, 60])
axes[0].set_ylim([-65, 10])

# --- 図2: v_circ, v_esc, sqrt(2)の関係 ---
h_km = np.linspace(0, 50000, 500)
r_plot = R_E + h_km * 1e3

v_c = np.sqrt(mu / r_plot) / 1e3
v_e = np.sqrt(2 * mu / r_plot) / 1e3
ratio = v_e / v_c  # 常に sqrt(2)

ax2_left = axes[1]
ax2_right = ax2_left.twinx()

ax2_left.plot(h_km, v_c, 'b-', linewidth=2, label='$v_{circ}$')
ax2_left.plot(h_km, v_e, 'r-', linewidth=2, label='$v_{esc}$')
ax2_right.plot(h_km, ratio, 'g--', linewidth=2, label='$v_{esc}/v_{circ}$')

ax2_left.set_xlabel('Altitude [km]', fontsize=12)
ax2_left.set_ylabel('Velocity [km/s]', fontsize=12, color='black')
ax2_right.set_ylabel('Ratio $v_{esc}/v_{circ}$', fontsize=12, color='green')
ax2_right.set_ylim([1.0, 2.0])
ax2_right.axhline(y=np.sqrt(2), color='green', linestyle=':', alpha=0.5)
ax2_right.text(40000, np.sqrt(2) + 0.02, f'$\\sqrt{{2}} \\approx {np.sqrt(2):.3f}$',
               fontsize=11, color='green')

lines1, labels1 = ax2_left.get_legend_handles_labels()
lines2, labels2 = ax2_right.get_legend_handles_labels()
ax2_left.legend(lines1 + lines2, labels1 + labels2, fontsize=10, loc='center right')

ax2_left.set_title('Circular and escape velocities', fontsize=14)
ax2_left.grid(True, alpha=0.3)

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

この2つの図から得られる知見を確認しましょう。

図1(左: エネルギーダイアグラム) は、ビザビバ方程式の「エネルギー保存」という本質を視覚化しています。黒い曲線は重力ポテンシャル $U(r) = -\mu/r$ で、距離が離れるにつれて0に漸近します。水平の破線は各軌道の全エネルギー $\varepsilon$ です。あるエネルギーレベルの軌道が存在できるのは、$\varepsilon > U(r)$ の範囲のみです(運動エネルギーが正でなければならないため)。LEO(青い破線)は深いポテンシャルの井戸に沈んでおり、脱出には大量のエネルギーが必要であることがわかります。一方、GEO(オレンジ色)はゼロに近く、脱出までの差が小さいことが見てとれます。青い塗りつぶし部分は、LEO軌道での運動エネルギー($K = \varepsilon – U$)を示しており、地球に近いほど運動エネルギーが大きい(速い)ことが視覚的に理解できます。

図2(右: 速度と比率) では、円軌道速度(青)と脱出速度(赤)がともに高度とともに減少する一方で、その比(緑の破線)は高度によらず常に $\sqrt{2} \approx 1.414$ であることが明瞭に示されています。これはビザビバ方程式から解析的に証明した結果の数値的な裏付けです。

Pythonでの実装: 双曲線軌道と超過速度

最後に、双曲線軌道を含む全4タイプの軌道形状を同一の座標系に描画し、それぞれの軌道上でのビザビバ方程式による速度プロファイルを可視化します。

import numpy as np
import matplotlib.pyplot as plt

mu = 3.986e14
R_E = 6.371e6

fig, axes = plt.subplots(1, 2, figsize=(16, 7))

# --- 図1: 4タイプの軌道形状 ---
p = 15000e3  # 半直弦 [m](全軌道で共通)

orbit_params = [
    (0.0, 'Circle ($e=0$)', 'blue', (-0.99, 0.99)),
    (0.5, 'Ellipse ($e=0.5$)', 'green', (-0.99, 0.99)),
    (1.0, 'Parabola ($e=1$)', 'orange', (-0.85, 0.85)),
    (1.5, 'Hyperbola ($e=1.5$)', 'red', (-0.7, 0.7)),
]

for e, label, color, (t_min_frac, t_max_frac) in orbit_params:
    theta = np.linspace(t_min_frac * np.pi, t_max_frac * np.pi, 1000)
    denom = 1 + e * np.cos(theta)
    r_orb = np.where(denom > 0.05, p / denom, np.nan)  # 発散を避ける
    r_orb = np.where(r_orb > 0, r_orb, np.nan)

    x = r_orb * np.cos(theta) / 1e6
    y = r_orb * np.sin(theta) / 1e6
    axes[0].plot(x, y, linewidth=2, label=label, color=color)

# 焦点(原点)にマーカー
axes[0].plot(0, 0, 'k*', markersize=15, zorder=5)
axes[0].text(1, -3, 'Focus\n(Earth)', ha='center', fontsize=10)

axes[0].set_xlabel('x [$\\times 10^3$ km]', fontsize=12)
axes[0].set_ylabel('y [$\\times 10^3$ km]', fontsize=12)
axes[0].set_title('Four types of Keplerian orbits\n(same semi-latus rectum $p$)', fontsize=13)
axes[0].legend(fontsize=10)
axes[0].set_aspect('equal')
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim([-50, 80])
axes[0].set_ylim([-60, 60])

# --- 図2: 速度 vs 距離(各軌道タイプ) ---
r_range = np.linspace(R_E + 200e3, 100000e3, 500)

# 円軌道
v_circ = np.sqrt(mu / r_range) / 1e3

# 楕円軌道 (a = 25000 km)
a_ell = 25000e3
r_ell = r_range[r_range <= 2 * a_ell]  # 楕円の範囲内のみ
v_ell = np.sqrt(mu * (2 / r_ell - 1 / a_ell)) / 1e3

# 放物線軌道 (a → ∞)
v_par = np.sqrt(2 * mu / r_range) / 1e3

# 双曲線軌道 (a = -30000 km)
a_hyp = -30000e3
v_hyp = np.sqrt(mu * (2 / r_range - 1 / a_hyp)) / 1e3

axes[1].plot(r_range / 1e6, v_circ, 'b-', linewidth=2, label='Circle ($a=r$)')
axes[1].plot(r_ell / 1e6, v_ell, 'g-', linewidth=2, label='Ellipse ($a=25000$ km)')
axes[1].plot(r_range / 1e6, v_par, 'orange', linewidth=2, label='Parabola ($a \\to \\infty$)')
axes[1].plot(r_range / 1e6, v_hyp, 'r-', linewidth=2, label='Hyperbola ($|a|=30000$ km)')

# v_∞ の漸近線
v_inf = np.sqrt(mu / abs(a_hyp)) / 1e3
axes[1].axhline(y=v_inf, color='red', linestyle=':', alpha=0.5)
axes[1].text(85, v_inf + 0.1, f'$v_\\infty = {v_inf:.2f}$ km/s',
             fontsize=10, color='red')

axes[1].set_xlabel('Distance $r$ [$\\times 10^3$ km]', fontsize=12)
axes[1].set_ylabel('Velocity [km/s]', fontsize=12)
axes[1].set_title('Velocity profiles for different orbit types', fontsize=14)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim([0, 100])
axes[1].set_ylim([0, 12])

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

この図からは、4種類の軌道におけるビザビバ方程式の振る舞いの違いが明確に読み取れます。

図1(左: 軌道形状) では、同じ半直弦 $p$ をもつ4タイプの軌道を重ねて描いています。円軌道(青)は閉じた曲線、楕円軌道(緑)はやや細長い閉曲線、放物線軌道(オレンジ)は開いた曲線で無限遠に向かい、双曲線軌道(赤)はさらに急な角度で開いています。離心率 $e$ が1を超えると軌道は閉じなくなり、衛星は二度と中心天体に戻ってこないことが形状から直感的に理解できます。

図2(右: 速度プロファイル) が特に重要です。全ての軌道で、距離 $r$ が大きくなるにつれて速度は減少しますが、その減り方が軌道タイプによって異なります。円軌道(青)は各距離での最小速度を示し、これが「その距離で周回するのに必要な速度」です。楕円軌道(緑)は近点で円軌道より速く、遠点で遅くなります。放物線軌道(オレンジ)は遠方で速度がゼロに漸近します。そして双曲線軌道(赤)は $r \to \infty$ でも有限の速度 $v_\infty$ が残り、これが赤い点線で示されています。$v_\infty \approx 3.65$ km/s は惑星間航行で重要な超過速度に対応します。

ビザビバ方程式の応用上の注意点

ビザビバ方程式は非常に汎用的ですが、使う際にはいくつかの前提条件を意識しておく必要があります。

二体問題の仮定

ビザビバ方程式は二体問題の枠組みで導出されています。つまり、中心天体と衛星の2つの質点だけを考え、他の天体からの重力は無視しています。実際の太陽系では月や太陽の重力が摂動として作用し、長期間では軌道が変化します。しかし、短期間の軌道計算や初期見積もりには十分な精度を提供します。

非球対称の効果

地球は完全な球ではなく、赤道方向にわずかに膨らんだ回転楕円体です。この非球対称性はJ2摂動として知られ、軌道要素の長期的な変動を引き起こします。ビザビバ方程式は点質量の仮定に基づくため、J2摂動の効果は含まれていません。精密な軌道予報にはこれらの摂動を考慮した数値計算が必要です。

推力がある場合

ビザビバ方程式はエネルギー保存則に基づくため、衛星に外力(推力)が加わっている間は適用できません。ロケットエンジンを噴射している間は力学的エネルギーが変化するためです。ただし、噴射前後の「自由飛行」区間ではビザビバ方程式が成り立つので、$\Delta v$ の計算(噴射前の軌道と噴射後の軌道をそれぞれビザビバ方程式で評価し、その差をとる)には非常に有用です。ホーマン遷移軌道の $\Delta v$ 計算はまさにこの方法で行います。

相対論的効果

GPS衛星のような高精度が求められるシステムでは、一般相対論的な効果(重力による時間の遅れ)を考慮する必要があります。しかし、軌道速度の計算自体に対する相対論的補正は非常に小さく($v/c \sim 10^{-5}$ 程度)、通常の軌道力学では無視できます。

これらの制約を理解した上で使えば、ビザビバ方程式は軌道力学における最も強力な「万能ツール」と言えるでしょう。

まとめ

本記事では、ビザビバ方程式をエネルギー保存則から丁寧に導出し、さまざまな軌道への応用を解説しました。

  • ビザビバ方程式 $v^2 = \mu(2/r – 1/a)$ は、軌道上の任意の位置 $r$ での速度 $v$ を、標準重力パラメータ $\mu$ と半長軸 $a$ だけから直ちに計算できる強力な公式です
  • エネルギー保存則と角運動量保存則から導出され、楕円軌道の全エネルギーは半長軸 $a$ のみで決まる($\varepsilon = -\mu/(2a)$)という深い結果に基づいています
  • 円軌道($a = r$)では $v = \sqrt{\mu/r}$(第一宇宙速度の一般化)
  • 脱出速度($a \to \infty$)では $v = \sqrt{2\mu/r}$(第二宇宙速度、円軌道の $\sqrt{2}$ 倍)
  • 双曲線軌道($a < 0$)では無限遠でも速度 $v_\infty = \sqrt{\mu/|a|}$ が残り、惑星間航行のミッション設計に不可欠です
  • ISS、静止軌道、月遷移軌道、火星ホーマン遷移軌道など、具体的なミッションへの適用を通じて、方程式の実用性を確認しました

ビザビバ方程式は軌道力学の「スイスアーミーナイフ」とも言える存在です。ケプラーの法則が軌道の「形」と「時間」を記述するのに対し、ビザビバ方程式は軌道の「速度」を記述します。この2つを組み合わせることで、軌道力学の問題の多くに対処することができます。

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