人工衛星が宇宙空間を飛んでいるとき、その機体はゆっくりと回転していることがあります。地上の私たちがフリスビーを投げると、回転しながら安定して飛ぶこともあれば、本を空中に投げ上げるとぐるぐると不安定に回転することもあります。同じ「回転」なのに、なぜこれほど振る舞いが異なるのでしょうか？
この問いに答えるのがオイラーの運動方程式です。ニュートンの運動方程式 $\bm{F} = m\bm{a}$ が質点の並進運動を記述するのと同様に、オイラーの運動方程式は剛体の回転運動を支配する基本法則です。
オイラーの運動方程式を理解すると、以下のような幅広い問題に取り組めるようになります。
- 人工衛星の姿勢制御: 衛星がどのように回転するかを予測し、望みの姿勢に制御する設計の出発点となります
- ジャイロスコープの歳差運動: 地球ゴマやジャイロセンサの首振り運動を数学的に記述できます
- 天体力学: 小惑星や彗星核のタンブリング運動を解析できます
- ロボティクス: 多関節ロボットの各リンクの回転ダイナミクスを定式化する基盤になります
本記事の内容
- 角運動量の時間変化とトルクの関係
- 慣性座標系から機体座標系への変換
- オイラーの運動方程式の導出
- 軸対称体における簡約化
- トルクフリー運動（自由回転）と歳差運動
- 回転安定性の解析（長軸回転 vs 短軸回転）
- Pythonによるトルクフリー運動のシミュレーションと3D可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
角運動量とトルク — 回転運動の基本法則
並進運動では「力 = 質量 × 加速度」が基本法則ですが、回転運動ではこれに対応する法則として角運動量の変化がトルクに等しいという関係が成り立ちます。まずはこの基本関係を確認しましょう。
角運動量の定義
剛体の角運動量 $\bm{H}$ は、慣性テンソル（慣性モーメント行列）$\bm{I}$ と角速度ベクトル $\bm{\omega}$ を用いて次のように定義されます。
$$ \bm{H} = \bm{I} \bm{\omega} $$
ここで $\bm{I}$ は $3 \times 3$ の対称行列であり、物体の質量分布と回転軸の関係を表します。対角成分は各軸まわりの慣性モーメント、非対角成分は慣性乗積です。
直感的には、角運動量は「回転の勢い」を表す量です。重いフライホイールが高速で回転しているほど、止めるのに大きな力が必要になります。これが角運動量の大きさに対応しています。
角運動量の変化則
慣性座標系（空間に固定された座標系）において、角運動量の時間変化はトルク（力のモーメント）に等しいという法則が成り立ちます。
$$ \left(\frac{d\bm{H}}{dt}\right)_{\text{inertial}} = \bm{N} $$
ここで $\bm{N}$ は外部トルク（重力傾度トルク、磁気トルク、スラスタなど）の合計です。
この式はニュートンの第2法則 $\bm{F} = \frac{d\bm{p}}{dt}$ の回転版です。並進運動で力が運動量を変化させるのと同様に、回転運動ではトルクが角運動量を変化させます。外部トルクが $\bm{N} = \bm{0}$ のとき、角運動量は保存されます。
この式は慣性座標系で書かれていますが、実際に衛星の回転を扱うときには機体座標系（衛星に固定された座標系）で方程式を記述する方が圧倒的に便利です。なぜなら、慣性テンソル $\bm{I}$ は機体座標系では定数行列になるのに対し、慣性座標系では物体の回転とともに時々刻々変化してしまうからです。
次のセクションでは、この慣性座標系の式を機体座標系に変換することで、オイラーの運動方程式を導出します。
機体座標系への変換 — オイラーの運動方程式の導出
回転座標系における時間微分
オイラーの運動方程式を導出する鍵は、慣性座標系と回転座標系（機体座標系）で時間微分が異なるという事実にあります。
任意のベクトル $\bm{A}$ について、慣性座標系での時間微分と機体座標系での時間微分の間には次の関係が成り立ちます。
$$ \left(\frac{d\bm{A}}{dt}\right)_{\text{inertial}} = \left(\frac{d\bm{A}}{dt}\right)_{\text{body}} + \bm{\omega} \times \bm{A} $$
ここで $\bm{\omega}$ は機体座標系の角速度ベクトルです。
この式は「輸送定理」と呼ばれ、直感的には次のように理解できます。回転する座標系に乗っている観測者が見るベクトルの変化（第1項）に加えて、座標系自体が回転していることによる見かけの変化（第2項 $\bm{\omega} \times \bm{A}$）を足し合わせると、慣性座標系から見た真の変化率になるのです。
たとえばメリーゴーラウンドに乗っている人が、手に持ったボールを「動いていない」と感じても、地面に立つ人から見ればボールは回転しています。この違いを記述するのが $\bm{\omega} \times \bm{A}$ の項です。
角運動量への適用
輸送定理を角運動量 $\bm{H}$ に適用しましょう。
$$ \left(\frac{d\bm{H}}{dt}\right)_{\text{inertial}} = \left(\frac{d\bm{H}}{dt}\right)_{\text{body}} + \bm{\omega} \times \bm{H} $$
左辺は先ほどの角運動量変化則から外部トルク $\bm{N}$ に等しいので、
$$ \bm{N} = \left(\frac{d\bm{H}}{dt}\right)_{\text{body}} + \bm{\omega} \times \bm{H} $$
機体座標系では慣性テンソル $\bm{I}$ は時間に依存しない定数行列です。これは機体に固定された座標系なので、質量分布と座標系の相対関係が変わらないためです。したがって、
$$ \left(\frac{d\bm{H}}{dt}\right)_{\text{body}} = \left(\frac{d(\bm{I}\bm{\omega})}{dt}\right)_{\text{body}} = \bm{I}\dot{\bm{\omega}} $$
ここで $\dot{\bm{\omega}}$ は機体座標系で見た角速度の時間変化率（角加速度）です。
オイラーの運動方程式
上の結果を代入すると、オイラーの運動方程式が得られます。
$$ \begin{equation} \bm{I}\dot{\bm{\omega}} + \bm{\omega} \times (\bm{I}\bm{\omega}) = \bm{N} \end{equation} $$
これが剛体の回転ダイナミクスを支配する基本方程式です。
第1項 $\bm{I}\dot{\bm{\omega}}$ は「角速度の変化による角運動量の変化」、第2項 $\bm{\omega} \times (\bm{I}\bm{\omega})$ は「回転座標系であることによるジャイロスコピック項（見かけの項）」に対応します。
並進運動のニュートンの第2法則 $m\dot{\bm{v}} = \bm{F}$ と比べると、回転版には $\bm{\omega} \times (\bm{I}\bm{\omega})$ という非線形項が加わっています。この項が回転運動を並進運動よりもはるかに豊かで複雑なものにしているのです。特に、慣性モーメントが異方的（3軸で異なる）な場合に、歳差運動やニューテーション運動といった興味深い現象が生じます。
次に、機体座標系として主軸系を選ぶことで、この方程式をさらに具体的に書き下してみましょう。
主軸座標系での成分表示
主軸座標系の選択
機体座標系として慣性テンソルの主軸系を採用すると、慣性テンソルは対角行列になります。
$$ \bm{I} = \begin{pmatrix} I_1 & 0 & 0 \\ 0 & I_2 & 0 \\ 0 & 0 & I_3 \end{pmatrix} $$
ここで $I_1, I_2, I_3$ はそれぞれの主軸まわりの主慣性モーメントです。
主軸系は慣性テンソルを固有値分解したときの固有ベクトル方向に対応します。多くの人工衛星は設計段階で質量分布が対称になるように作られるため、機体の幾何学的な対称軸がそのまま主軸に一致します。
角速度ベクトルを主軸系の成分で $\bm{\omega} = (\omega_1, \omega_2, \omega_3)^T$ と表すと、オイラーの運動方程式は次の3本の連立常微分方程式になります。
$$ \begin{align} I_1 \dot{\omega}_1 + (I_3 – I_2) \omega_2 \omega_3 &= N_1 \\ I_2 \dot{\omega}_2 + (I_1 – I_3) \omega_3 \omega_1 &= N_2 \\ I_3 \dot{\omega}_3 + (I_2 – I_1) \omega_1 \omega_2 &= N_3 \end{align} $$
各項の物理的意味
1つ目の方程式を例に取ると、
- $I_1 \dot{\omega}_1$: 第1軸まわりの角加速度による角運動量変化
- $(I_3 – I_2)\omega_2 \omega_3$: ジャイロスコピック結合（第2軸と第3軸の回転が第1軸に影響を与える項）
- $N_1$: 第1軸まわりの外部トルク
ジャイロスコピック結合項 $(I_3 – I_2)\omega_2 \omega_3$ は、慣性モーメントが等方的（$I_1 = I_2 = I_3$）でない限り現れます。この非線形結合こそが、回転運動の複雑さの源です。
もし慣性モーメントが全て等しければ、ジャイロスコピック項は消えて $I\dot{\omega}_i = N_i$ という単純な式になります。これは各軸が独立に振る舞い、並進運動の $ma = F$ と同じ構造です。しかし現実の物体はほぼ確実に異方的な慣性モーメントを持つため、3軸は常に結合しています。
特に重要なのは、外部トルクがゼロの場合（$\bm{N} = \bm{0}$）でも角速度が変化しうるという点です。並進運動では力がなければ速度は一定ですが、回転運動ではトルクがなくても角速度の方向と大きさが変化します。これがトルクフリー運動であり、次に詳しく解析します。
軸対称体の場合 — 方程式の簡約化
多くの人工衛星や回転機械は、1つの軸に対して回転対称な形状をしています。たとえばスピン安定衛星は円筒形やドラム形であり、回転軸まわりに対称です。このような軸対称体では、オイラーの運動方程式が大幅に簡略化されます。
軸対称体の慣性モーメント
軸対称体の対称軸を第3軸に取ると、残りの2軸の慣性モーメントが等しくなります。
$$ I_1 = I_2 = I_T, \quad I_3 = I_S $$
ここで $I_T$ は横軸（transverse）の慣性モーメント、$I_S$ は対称軸（spin/symmetry）の慣性モーメントです。
このとき、オイラーの運動方程式は次のように簡約化されます。
$$ \begin{align} I_T \dot{\omega}_1 + (I_S – I_T) \omega_2 \omega_3 &= N_1 \\ I_T \dot{\omega}_2 – (I_S – I_T) \omega_1 \omega_3 &= N_2 \\ I_S \dot{\omega}_3 &= N_3 \end{align} $$
注目すべきは第3式です。対称軸まわりの回転には $(I_2 – I_1) = 0$ から結合項が消え、対称軸まわりのスピン角速度 $\omega_3$ は外部トルク $N_3$ のみで決まります。つまり、対称軸のスピンは他の2軸の運動と独立しています。
この性質は、スピン安定衛星にとって極めて重要です。対称軸まわりに一定速度でスピンさせれば、そのスピン速度はスピン軸まわりのトルクがない限り変化しません。
さらに、第1式と第2式に現れる $(I_S – I_T)\omega_3$ をひとまとめにすると、横軸の運動を複素平面で記述できるようになります。 $\Omega = \omega_1 + i\omega_2$ と定義すると、トルクフリーの場合に
$$ \dot{\Omega} + i\lambda \Omega = 0, \quad \lambda = \frac{I_S – I_T}{I_T}\omega_3 $$
という1次の線形常微分方程式になり、解は
$$ \Omega(t) = \Omega_0 e^{-i\lambda t} $$
です。これは横軸の角速度成分が一定の振幅で回転することを意味し、歳差運動（プリセッション）を記述しています。
それでは、トルクフリー運動の一般的な解析に進みましょう。
トルクフリー運動 — 外部トルクがない剛体の自由回転
基本設定
宇宙空間で外部トルクを受けない剛体（$\bm{N} = \bm{0}$）の回転は、剛体力学の古典的問題であり、オイラーの運動方程式から多くの知見が得られます。
トルクフリーの場合、オイラーの運動方程式は次のようになります。
$$ \begin{align} I_1 \dot{\omega}_1 &= (I_2 – I_3) \omega_2 \omega_3 \\ I_2 \dot{\omega}_2 &= (I_3 – I_1) \omega_3 \omega_1 \\ I_3 \dot{\omega}_3 &= (I_1 – I_2) \omega_1 \omega_2 \end{align} $$
保存量
トルクフリー運動では、2つの重要な保存量が存在します。
1. 角運動量の大きさ
外部トルクがゼロなので、慣性座標系での角運動量ベクトル $\bm{H}$ は時間変化しません（大きさも方向も保存）。
$$ |\bm{H}|^2 = I_1^2 \omega_1^2 + I_2^2 \omega_2^2 + I_3^2 \omega_3^2 = \text{const} $$
2. 回転の運動エネルギー
外部トルクがなく、剛体なのでエネルギー散逸もないため、回転の運動エネルギーも保存されます。
$$ 2T = I_1 \omega_1^2 + I_2 \omega_2^2 + I_3 \omega_3^2 = \text{const} $$
これら2つの保存量は、$(\omega_1, \omega_2, \omega_3)$ の3次元空間において楕円体の表面として表されます。角運動量の保存は「角運動量楕円体」、エネルギーの保存は「エネルギー楕円体（慣性楕円体）」を定義し、角速度ベクトルの先端はこれら2つの楕円体の交線上を動きます。この幾何学的な描像をポアンソの構成と呼びます。
軸対称体のトルクフリー運動
軸対称体（$I_1 = I_2 = I_T$, $I_3 = I_S$）の場合、前セクションで導出したように解析解が得られます。
第3式から $\dot{\omega}_3 = 0$、つまりスピン角速度は一定です。
$$ \omega_3 = \omega_{30} = \text{const} $$
横軸の角速度は、
$$ \omega_1(t) = A \cos(\lambda t + \phi_0), \quad \omega_2(t) = A \sin(\lambda t + \phi_0) $$
$$ \lambda = \frac{I_S – I_T}{I_T}\omega_{30} $$
ここで $A$ は横軸角速度の振幅、$\phi_0$ は初期位相です。
$\lambda$ は機体ニューテーション周波数と呼ばれます。角速度ベクトル $\bm{\omega}$ の先端は機体座標系では第3軸（スピン軸）を中心とする円を描きます。
2つの円錐運動
トルクフリー運動では、角速度ベクトルが描く軌跡を「円錐運動」として幾何学的に理解できます。
- ボディコーン（機体円錐）: 角速度ベクトルが機体座標系のスピン軸のまわりを回転する円錐。機体から見た運動で、周波数は $\lambda$ です。
- スペースコーン（空間円錐）: 角速度ベクトルが慣性空間での角運動量ベクトル $\bm{H}$ のまわりを回転する円錐。宇宙から見た運動です。
機体のスピン軸（第3軸）は、角運動量ベクトル $\bm{H}$（慣性空間で固定）のまわりを歳差運動します。日常で目にするジャイロスコープの首振り運動はまさにこの現象です。
- $I_S > I_T$（扁平な物体、例: 円盤）の場合: ボディコーンがスペースコーンの外側を転がる（外接転がり）。歳差は角速度ベクトルと同方向。
- $I_S < I_T$（細長い物体、例: ロケット）の場合: ボディコーンがスペースコーンの内側を転がる（内接転がり）。歳差は角速度ベクトルと逆方向。
トルクフリー運動では、角速度ベクトルの大きさは一定（エネルギー保存）ですが、その方向は機体座標系でも慣性座標系でも変化します。唯一固定なのは、慣性座標系における角運動量ベクトル $\bm{H}$ の方向です。
それでは次に、一般の非対称剛体（$I_1 \neq I_2 \neq I_3$）の場合に、どの軸まわりの回転が安定で、どの軸が不安定かを解析しましょう。
回転安定性の解析 — テニスラケットの定理
問題設定
3軸の慣性モーメントが全て異なる剛体（$I_1 < I_2 < I_3$ と仮定）が、主軸のいずれかのまわりにほぼ純粋に回転しているとき、その回転は安定でしょうか？
日常で実験できます。本やスマートフォンを空中に投げ上げてみてください。
- 最も短い軸（最大慣性モーメント $I_3$）まわり: 安定に回転する
- 最も長い軸（最小慣性モーメント $I_1$）まわり: 安定に回転する
- 中間の軸（中間慣性モーメント $I_2$）まわり: 不安定にタンブリングする
この現象はテニスラケットの定理（またはジャニベコフ効果、中間軸の定理）と呼ばれます。
線形安定性解析
第1軸まわりの回転を例に解析します。定常回転を $\bm{\omega}_0 = (\Omega, 0, 0)^T$ とし、微小摂動 $\bm{\omega} = (\Omega + \delta\omega_1, \delta\omega_2, \delta\omega_3)^T$ を加えます。
オイラーの運動方程式に代入し、摂動の2次以上の項を無視すると、
$I_1$ 軸の式からは $\delta\dot{\omega}_1 \approx 0$ となり、摂動は成長しません。
$I_2$ 軸と $I_3$ 軸の式から連立方程式が得られます。
$$ I_2 \delta\dot{\omega}_2 = (I_3 – I_1) \Omega \, \delta\omega_3 $$
$$ I_3 \delta\dot{\omega}_3 = (I_1 – I_2) \Omega \, \delta\omega_2 $$
$\delta\omega_2$ を時間微分して $\delta\dot{\omega}_3$ を代入すると、
$$ \delta\ddot{\omega}_2 = \frac{(I_3 – I_1)(I_1 – I_2)}{I_2 I_3}\Omega^2 \, \delta\omega_2 $$
$\sigma^2 = -\frac{(I_3 – I_1)(I_1 – I_2)}{I_2 I_3}\Omega^2$ と置くと、
$$ \delta\ddot{\omega}_2 + \sigma^2 \delta\omega_2 = 0 $$
安定性条件
$\sigma^2 > 0$ なら $\delta\omega_2$ は正弦振動（安定）、$\sigma^2 < 0$ なら指数的に発散（不安定）します。
$I_1 < I_2 < I_3$ のとき、各軸について $\sigma^2$ の符号を調べます。
第1軸まわり（最小慣性モーメント $I_1$）:
$$ \sigma^2 = -\frac{(I_3 – I_1)(I_1 – I_2)}{I_2 I_3}\Omega^2 $$
$I_3 – I_1 > 0$, $I_1 – I_2 < 0$ より $(I_3-I_1)(I_1-I_2) < 0$。したがって $\sigma^2 > 0$ で安定です。
第3軸まわり（最大慣性モーメント $I_3$）:
同様の解析で $\sigma^2 > 0$ となり安定です。$(I_1 – I_3)(I_3 – I_2)$ を計算すると、$I_1 – I_3 < 0$, $I_3 - I_2 > 0$ より積は負、$-$ を掛けて正。
第2軸まわり（中間慣性モーメント $I_2$）:
$(I_3 – I_2)(I_2 – I_1)$ は両方とも正なので積は正、$-$ を掛けて $\sigma^2 < 0$ で不安定です。
まとめ: テニスラケットの定理
| 回転軸 | 慣性モーメント | 安定性 |
|---|---|---|
| 第1軸 | 最小 $I_1$ | 安定 |
| 第2軸 | 中間 $I_2$ | 不安定 |
| 第3軸 | 最大 $I_3$ | 安定 |
中間軸まわりの回転のみが不安定という結果が得られました。
この定理は人工衛星の設計に直結します。スピン安定衛星はスピン軸の慣性モーメントが最大（短軸回転、扁平体）または最小（長軸回転）になるように設計しなければなりません。中間軸まわりにスピンさせると不安定になり、制御を失う危険があります。
ただし、実際の衛星にはエネルギー散逸（内部の液体、柔軟部材の振動など）が存在します。エネルギー散逸がある場合、角運動量は保存されますがエネルギーは減少するため、最終的に最大慣性モーメント軸まわりの回転に遷移します。したがって、長期的に安定なスピン安定化には $I_S > I_T$（短軸回転、メジャーアクシス・スピン）が必要です。初期の Explorer I 衛星は長軸回転（$I_S < I_T$）で打ち上げられましたが、柔軟アンテナのエネルギー散逸により短軸回転に遷移してしまった歴史的事例があります。
理論的な解析を一通り終えたところで、次はPythonでトルクフリー運動をシミュレーションし、これらの現象を視覚的に確認しましょう。
Pythonによるトルクフリー運動のシミュレーション
オイラー方程式の数値解法
まず、非対称剛体のトルクフリー運動を数値的に解きます。3つの主軸の慣性モーメントを異なる値に設定し、各軸まわりの初期回転に対する安定性を確認します。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# 主慣性モーメント (I1 < I2 < I3)
I1, I2, I3 = 1.0, 2.0, 3.0
def euler_equations(t, omega):
"""トルクフリーのオイラー方程式"""
w1, w2, w3 = omega
dw1 = (I2 - I3) / I1 * w2 * w3
dw2 = (I3 - I1) / I2 * w3 * w1
dw3 = (I1 - I2) / I3 * w1 * w2
return [dw1, dw2, dw3]
t_span = (0, 30)
t_eval = np.linspace(0, 30, 3000)
# ケース1: 第1軸まわり（最小I）+ 微小摂動 → 安定
omega0_axis1 = [5.0, 0.1, 0.1]
sol1 = solve_ivp(euler_equations, t_span, omega0_axis1, t_eval=t_eval, rtol=1e-10, atol=1e-12)
# ケース2: 第2軸まわり（中間I）+ 微小摂動 → 不安定
omega0_axis2 = [0.1, 5.0, 0.1]
sol2 = solve_ivp(euler_equations, t_span, omega0_axis2, t_eval=t_eval, rtol=1e-10, atol=1e-12)
# ケース3: 第3軸まわり（最大I）+ 微小摂動 → 安定
omega0_axis3 = [0.1, 0.1, 5.0]
sol3 = solve_ivp(euler_equations, t_span, omega0_axis3, t_eval=t_eval, rtol=1e-10, atol=1e-12)
fig, axes = plt.subplots(3, 1, figsize=(10, 10))
titles = [
f'Rotation about axis 1 (I={I1}, minimum) - Stable',
f'Rotation about axis 2 (I={I2}, intermediate) - Unstable',
f'Rotation about axis 3 (I={I3}, maximum) - Stable'
]
solutions = [sol1, sol2, sol3]
colors = ['#e74c3c', '#2ecc71', '#3498db']
for ax, sol, title in zip(axes, solutions, titles):
ax.plot(sol.t, sol.y[0], label=r'$\omega_1$', color='#e74c3c', linewidth=1.2)
ax.plot(sol.t, sol.y[1], label=r'$\omega_2$', color='#2ecc71', linewidth=1.2)
ax.plot(sol.t, sol.y[2], label=r'$\omega_3$', color='#3498db', linewidth=1.2)
ax.set_title(title, fontsize=12)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Angular velocity [rad/s]')
ax.legend(loc='upper right')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('euler_torque_free.png', dpi=150, bbox_inches='tight')
plt.show()
このシミュレーション結果から、テニスラケットの定理が鮮明に確認できます。
- 第1軸まわり（最小慣性モーメント）: $\omega_1 \approx 5$ を維持しながら、$\omega_2, \omega_3$ は微小振動にとどまります。初期摂動が増幅されず、安定な回転です。
- 第2軸まわり（中間慣性モーメント）: 初期に $\omega_2 \approx 5$ だったのが、時間とともに大きく変動し、角速度ベクトルが3軸全体に広がるタンブリング運動を起こします。中間軸まわりの回転が不安定であることが明瞭です。
- 第3軸まわり（最大慣性モーメント）: $\omega_3 \approx 5$ を維持し、他の2成分は微小振動にとどまります。第1軸と同様に安定です。
保存量の確認
トルクフリー運動では角運動量とエネルギーが保存されるはずです。数値誤差の確認も兼ねて、保存量をプロットします。
fig, axes = plt.subplots(2, 1, figsize=(10, 6))
labels = ['Axis 1 (stable)', 'Axis 2 (unstable)', 'Axis 3 (stable)']
linestyles = ['-', '--', ':']
for sol, label, ls in zip(solutions, labels, linestyles):
w1, w2, w3 = sol.y
# 角運動量の大きさ
H = np.sqrt((I1*w1)**2 + (I2*w2)**2 + (I3*w3)**2)
H_rel = (H - H[0]) / H[0]
axes[0].plot(sol.t, H_rel, label=label, linestyle=ls, linewidth=1.5)
# 運動エネルギー
T = 0.5 * (I1*w1**2 + I2*w2**2 + I3*w3**2)
T_rel = (T - T[0]) / T[0]
axes[1].plot(sol.t, T_rel, label=label, linestyle=ls, linewidth=1.5)
axes[0].set_title('Conservation of angular momentum (relative error)')
axes[0].set_ylabel(r'$(|H| - |H_0|)/|H_0|$')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[1].set_title('Conservation of kinetic energy (relative error)')
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel(r'$(T - T_0)/T_0$')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('euler_conservation.png', dpi=150, bbox_inches='tight')
plt.show()
保存量の相対誤差はいずれのケースでも $10^{-10}$ 以下のオーダーに収まっていることが確認できます。これは solve_ivp の厳しい許容誤差設定（rtol=1e-10）の効果です。不安定な中間軸の場合でも保存量は正確に維持されており、タンブリング運動は数値誤差ではなく、物理的に正しい不安定性を反映しています。
角速度ベクトルの3D軌跡（ポアンソの構成）
角速度ベクトルの先端がどのような軌跡を描くかを3Dで可視化します。角運動量楕円体とエネルギー楕円体の交線として表現されるこの軌跡は、ポアンソの構成の核心です。
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(12, 5))
for idx, (sol, title) in enumerate(zip(solutions,
['Axis 1 (stable)', 'Axis 2 (unstable)', 'Axis 3 (stable)'])):
ax = fig.add_subplot(1, 3, idx+1, projection='3d')
w1, w2, w3 = sol.y
ax.plot(w1, w2, w3, linewidth=0.8, alpha=0.8)
ax.scatter(w1[0], w2[0], w3[0], color='red', s=40, label='Start')
ax.set_xlabel(r'$\omega_1$')
ax.set_ylabel(r'$\omega_2$')
ax.set_zlabel(r'$\omega_3$')
ax.set_title(title, fontsize=10)
ax.legend(fontsize=8)
plt.tight_layout()
plt.savefig('euler_3d_trajectory.png', dpi=150, bbox_inches='tight')
plt.show()
3Dプロットから、各ケースの角速度ベクトルの軌跡がはっきりと異なることがわかります。
- 第1軸まわり（安定）: 角速度ベクトルは第1軸の近傍で小さな閉曲線を描きます。摂動が増幅されず、元の回転状態の近くにとどまっています。
- 第2軸まわり（不安定）: 角速度ベクトルは $\omega_1$-$\omega_3$ 平面を大きく横切る軌跡を描き、第2軸から遠く離れた領域まで到達します。これがタンブリング運動の3D的な姿です。
- 第3軸まわり（安定）: 第1軸と同様に、スピン軸近傍の小さな閉曲線にとどまります。
軸対称体の歳差運動
最後に、軸対称体（$I_1 = I_2$）のトルクフリー運動で歳差運動（プリセッション）を可視化します。角速度ベクトルがスピン軸のまわりを回転する様子と、解析解との一致を確認します。
# 軸対称体 (IT = IT, IS)
IT = 2.0 # 横軸慣性モーメント
IS = 3.0 # スピン軸慣性モーメント (扁平体: IS > IT)
def euler_axisymmetric(t, omega):
w1, w2, w3 = omega
dw1 = (IT - IS) / IT * w2 * w3
dw2 = (IS - IT) / IT * w1 * w3
dw3 = 0.0 # 対称軸のスピンは一定
return [dw1, dw2, dw3]
# 初期条件: スピン + 横軸摂動
omega_s0 = 10.0 # スピン角速度
omega_t0 = 2.0 # 横軸角速度の振幅
omega0 = [omega_t0, 0.0, omega_s0]
t_span = (0, 20)
t_eval = np.linspace(0, 20, 2000)
sol = solve_ivp(euler_axisymmetric, t_span, omega0, t_eval=t_eval, rtol=1e-12, atol=1e-14)
# 解析解
lam = (IS - IT) / IT * omega_s0 # ニューテーション周波数
omega1_ana = omega_t0 * np.cos(lam * t_eval)
omega2_ana = omega_t0 * np.sin(lam * t_eval)
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
# 時間波形
axes[0, 0].plot(sol.t, sol.y[0], 'b-', label=r'$\omega_1$ (numerical)', linewidth=1.5)
axes[0, 0].plot(t_eval, omega1_ana, 'r--', label=r'$\omega_1$ (analytical)', linewidth=1.0)
axes[0, 0].set_xlabel('Time [s]')
axes[0, 0].set_ylabel(r'$\omega_1$ [rad/s]')
axes[0, 0].set_title('Transverse angular velocity (axis 1)')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
axes[0, 1].plot(sol.t, sol.y[1], 'b-', label=r'$\omega_2$ (numerical)', linewidth=1.5)
axes[0, 1].plot(t_eval, omega2_ana, 'r--', label=r'$\omega_2$ (analytical)', linewidth=1.0)
axes[0, 1].set_xlabel('Time [s]')
axes[0, 1].set_ylabel(r'$\omega_2$ [rad/s]')
axes[0, 1].set_title('Transverse angular velocity (axis 2)')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# omega3 の安定性
axes[1, 0].plot(sol.t, sol.y[2], 'g-', linewidth=1.5)
axes[1, 0].set_xlabel('Time [s]')
axes[1, 0].set_ylabel(r'$\omega_3$ [rad/s]')
axes[1, 0].set_title(f'Spin angular velocity (constant = {omega_s0})')
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].set_ylim([omega_s0 - 0.5, omega_s0 + 0.5])
# 横軸の位相面 (歳差の円)
axes[1, 1].plot(sol.y[0], sol.y[1], 'b-', linewidth=1.0, label='Numerical')
axes[1, 1].plot(omega1_ana, omega2_ana, 'r--', linewidth=0.8, label='Analytical')
axes[1, 1].set_xlabel(r'$\omega_1$ [rad/s]')
axes[1, 1].set_ylabel(r'$\omega_2$ [rad/s]')
axes[1, 1].set_title('Precession in transverse plane')
axes[1, 1].set_aspect('equal')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.suptitle(f'Axisymmetric body: $I_T$={IT}, $I_S$={IS}, Nutation freq $\\lambda$={lam:.2f} rad/s',
fontsize=13, y=1.02)
plt.tight_layout()
plt.savefig('euler_precession.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"ニューテーション周波数 λ = {lam:.4f} rad/s")
print(f"ニューテーション周期 T = {2*np.pi/abs(lam):.4f} s")
print(f"数値解と解析解の最大誤差 (ω1): {np.max(np.abs(sol.y[0] - omega1_ana)):.2e}")
print(f"数値解と解析解の最大誤差 (ω2): {np.max(np.abs(sol.y[1] - omega2_ana)):.2e}")
このシミュレーション結果から、軸対称体の歳差運動の特徴が明確に読み取れます。
- 横軸角速度 $\omega_1, \omega_2$ は正弦波: 数値解（青実線）と解析解（赤破線）が完全に一致しています。最大誤差は $10^{-11}$ のオーダーであり、解析解が正しいことが数値的にも裏付けられています。
- スピン角速度 $\omega_3$ は厳密に一定: 軸対称体の性質として、スピン角速度はジャイロスコピック結合の影響を受けません。
- 位相面上の円軌道: $\omega_1$-$\omega_2$ 平面での軌跡は完全な円を描きます。これは横軸の角速度が一定振幅で回転していることを示し、ニューテーション周波数 $\lambda = (I_S – I_T)/I_T \cdot \omega_3$ で円運動します。
- ニューテーション周波数: $\lambda = 5.0$ rad/s、周期 $T \approx 1.26$ sです。スピン角速度とモーメント比 $(I_S – I_T)/I_T = 0.5$ で決まるこの値は、衛星設計における姿勢安定性の評価に直結します。
エネルギー散逸と衛星設計への示唆
ここまでの議論は剛体（エネルギー散逸なし）を前提としていましたが、実際の衛星にはエネルギー散逸の機構が存在します。液体燃料の残留、柔軟な太陽電池パドル、ケーブルハーネスの振動などが典型的な散逸源です。
角運動量は外部トルクがない限り厳密に保存されますが、運動エネルギーは内部散逸により減少していきます。
$$ |\bm{H}|^2 = I_1^2\omega_1^2 + I_2^2\omega_2^2 + I_3^2\omega_3^2 = \text{const} $$
$$ T = \frac{1}{2}(I_1\omega_1^2 + I_2\omega_2^2 + I_3\omega_3^2) \to \text{min} $$
角運動量が一定のまま運動エネルギーを最小化する回転状態を求めると、ラグランジュの未定乗数法から最大慣性モーメント軸まわりの純粋回転が解になります。
$$ T_{\min} = \frac{|\bm{H}|^2}{2I_{\max}} $$
したがって、エネルギー散逸がある場合:
- 短軸回転（メジャーアクシス・スピン: $I_S = I_{\max}$）: エネルギーが最小の状態であり、散逸があっても安定。衛星は最終的にこの状態に落ち着きます。
- 長軸回転（マイナーアクシス・スピン: $I_S = I_{\min}$）: エネルギーが最大の状態であり、散逸によってエネルギーが減少すると中間軸・最大軸方向への遷移が起こります。
この知見は、1958年に打ち上げられた Explorer I 衛星の事例で劇的に実証されました。細長い形状（$I_S < I_T$）の衛星を長軸スピンで打ち上げたところ、柔軟アンテナの振動によるエネルギー散逸のため、数時間で横軸回転（タンブリング）に遷移してしまったのです。この教訓から、スピン安定衛星は $I_S > I_T$（短軸回転）で設計するのが鉄則となりました。
以上の解析を踏まえると、衛星の姿勢制御設計において重要なのは、オイラーの運動方程式を基盤として角運動量の管理戦略を構築することです。次の記事では、この角運動量の管理問題を衛星全体のシステムレベルで扱います。
まとめ
本記事では、剛体の回転ダイナミクスを記述するオイラーの運動方程式について、基礎理論から導出、解析、シミュレーションまでを解説しました。
- オイラーの運動方程式: $\bm{I}\dot{\bm{\omega}} + \bm{\omega} \times (\bm{I}\bm{\omega}) = \bm{N}$ — 機体座標系（回転座標系）で記述することで慣性テンソルが定数行列になり、扱いやすくなります
- ジャイロスコピック項: $\bm{\omega} \times (\bm{I}\bm{\omega})$ は慣性モーメントの異方性に起因する非線形結合項であり、回転運動を並進運動よりも複雑にする本質的な要素です
- 軸対称体の簡約化: $I_1 = I_2$ のとき方程式は大幅に簡約化され、スピン軸角速度は保存され、横軸ではニューテーション周波数 $\lambda$ による歳差運動が生じます
- テニスラケットの定理: 最大・最小慣性モーメント軸まわりの回転は安定、中間軸まわりの回転は不安定です
- エネルギー散逸の効果: 実際の衛星ではエネルギー散逸により最大慣性モーメント軸まわりの回転に遷移するため、スピン安定衛星は短軸回転で設計する必要があります
次のステップとして、以下の記事も参考にしてください。