慣性モーメント — 回転運動を支配するテンソル量

フィギュアスケーターがスピンの最中に腕を縮めると、回転が急に速くなります。腕を広げると回転はゆっくりになります。同じトルク(力のモーメント)が作用しているわけではないのに、なぜ回転速度が変わるのでしょうか。

答えは 慣性モーメント にあります。並進運動における「質量」に相当するのが、回転運動における「慣性モーメント」です。質量が大きい物体を加速しにくいように、慣性モーメントが大きい軸まわりに物体を回転させることは難しいのです。そして3次元の物体では、回転軸の方向によって慣性モーメントの値が異なるため、テンソル量 として記述する必要があります。

慣性モーメントとテンソルを理解すると、以下の応用に直結します。

  • 人工衛星の姿勢制御: 衛星の慣性テンソルが制御系の応答特性を決定する。設計段階で慣性テンソルを最適化することがミッション成功の鍵
  • ジャイロスコープ・リアクションホイール: 角運動量の保存則と慣性モーメントの関係がこれらのデバイスの動作原理
  • 機械設計: エンジンのフライホイール、車輪の設計で回転慣性の制御が不可欠
  • ロボットの動力学: 各リンクの慣性テンソルが関節トルク計算の基礎

本記事の内容

  • 慣性モーメントのスカラー定義からテンソルへの拡張
  • 慣性テンソルの対角化と主軸
  • 慣性積の物理的意味
  • 平行軸の定理(シュタイナーの定理)
  • 代表的な形状の慣性モーメント
  • 衛星の慣性テンソル設計
  • Pythonでの慣性テンソルの固有値・主軸可視化と慣性楕円体の3Dプロット

前提知識

この記事を読む前に、以下の知識があると理解がスムーズです。

  • 線形代数の基礎(行列の固有値・固有ベクトル、対称行列の対角化)
  • 積分の基本(体積積分)
  • 回転行列の基礎

慣性モーメントのスカラー定義

回転運動の基本法則

まず、1つの固定軸まわりの回転から始めましょう。ある軸まわりの並進運動の運動方程式 $F = ma$ に対応する回転の運動方程式は、

$$ \tau = I \alpha $$

ここで $\tau$ はトルク、$\alpha$ は角加速度、$I$ は慣性モーメントです。質量 $m$ が「動かしにくさ」を表すように、慣性モーメント $I$ は「回しにくさ」を表します。

スカラーの慣性モーメント

ある回転軸に対する慣性モーメントは、質量要素 $dm$ と回転軸からの距離 $r_{\perp}$ を使って次のように定義されます。

$$ I = \int r_{\perp}^2 \, dm $$

離散的な質点系の場合は、

$$ I = \sum_{i} m_i r_{\perp,i}^2 $$

$r_{\perp}^2$ という2乗の因子が入っていることに注目してください。回転軸から遠い質量要素ほど、慣性モーメントへの寄与が距離の2乗で大きくなります。だからこそ、スケーターが腕を広げると慣性モーメントが大きくなり、回転が遅くなるのです。

具体例: 2つの質点

質量 $m$ の2つの質点が距離 $L$ だけ離れているとき、中心を通る回転軸に対する慣性モーメントは、

$$ I = m \left(\frac{L}{2}\right)^2 + m \left(\frac{L}{2}\right)^2 = \frac{mL^2}{2} $$

同じ全質量 $2m$ でも、$L$ を大きくすれば慣性モーメントは増加します。

スカラーの慣性モーメントは、特定の1本の回転軸に対してしか定義できません。3次元空間では回転軸の方向は無数にあり、方向によって慣性モーメントの値は異なります。この方向依存性を統一的に記述するのが、次に説明する慣性テンソルです。

慣性テンソルの定義

角運動量からの動機づけ

慣性テンソルの必要性を理解するために、剛体の角運動量 $\bm{L}$ を考えましょう。角速度 $\bm{\omega}$ で回転する剛体の角運動量は、

$$ \bm{L} = \int \bm{r} \times (\bm{\omega} \times \bm{r}) \, dm $$

ベクトル三重積の公式 $\bm{a} \times (\bm{b} \times \bm{c}) = \bm{b}(\bm{a} \cdot \bm{c}) – \bm{c}(\bm{a} \cdot \bm{b})$ を適用すると、

$$ \bm{L} = \int \left[\bm{\omega}(\bm{r} \cdot \bm{r}) – \bm{r}(\bm{r} \cdot \bm{\omega})\right] dm $$

この式は $\bm{\omega}$ について線形なので、行列の形 $\bm{L} = \bm{J} \bm{\omega}$ で書けるはずです。成分表示で整理すると、

$$ L_i = \sum_{j} J_{ij} \omega_j $$

ここで $\bm{J}$ の各成分は以下のように定義されます。

慣性テンソルの成分

$$ J_{ij} = \int \left(|\bm{r}|^2 \delta_{ij} – r_i r_j\right) dm $$

具体的に書き下すと、

対角成分(慣性モーメント):

$$ J_{xx} = \int (y^2 + z^2) \, dm, \quad J_{yy} = \int (x^2 + z^2) \, dm, \quad J_{zz} = \int (x^2 + y^2) \, dm $$

対角成分はそれぞれ $x$, $y$, $z$ 軸まわりの慣性モーメントです。$J_{xx}$ の定義に $x^2$ が含まれず、$y^2 + z^2$($x$ 軸からの距離の2乗)になっていることに注意してください。

非対角成分(慣性積):

$$ J_{xy} = J_{yx} = -\int xy \, dm, \quad J_{xz} = J_{zx} = -\int xz \, dm, \quad J_{yz} = J_{zy} = -\int yz \, dm $$

慣性積にはマイナス符号が付いていることに注意してください。これは慣用的な定義で、角運動量の式 $\bm{L} = \bm{J}\bm{\omega}$ が自然な形になるようにするためです。

テンソルとしての記法

慣性テンソル $\bm{J}$ を行列として書くと、

$$ \bm{J} = \begin{pmatrix} J_{xx} & J_{xy} & J_{xz} \\ J_{yx} & J_{yy} & J_{yz} \\ J_{zx} & J_{zy} & J_{zz} \end{pmatrix} $$

$\bm{J}$ は 実対称行列 です($J_{ij} = J_{ji}$)。対称性は定義式から明らかで、$r_i r_j = r_j r_i$ であることから導かれます。

実対称行列であることは、線形代数の観点から非常に重要な意味を持ちます。固有値が全て実数であり、固有ベクトルが直交基底を構成するからです。この性質が慣性テンソルの対角化と主軸の概念につながります。

角運動量の表式

慣性テンソルを使えば、角運動量は簡潔に書けます。

$$ \bm{L} = \bm{J} \bm{\omega} $$

回転の運動エネルギーは、

$$ T = \frac{1}{2} \bm{\omega}^T \bm{J} \bm{\omega} = \frac{1}{2} \bm{\omega} \cdot \bm{L} $$

重要な事実: 一般に $\bm{L}$ と $\bm{\omega}$ は平行ではありません。角速度が $x$ 方向を向いていても、慣性積が存在すると角運動量は $x$ 方向からずれます。$\bm{L}$ と $\bm{\omega}$ が平行になる特別な方向が、次に説明する「主軸」です。

慣性テンソルの一般的な形を理解しました。では、非対角成分をゼロにできる特別な座標系は存在するのでしょうか。それが主軸座標系です。

慣性テンソルの対角化と主軸

主軸とは

主軸 とは、その軸まわりに回転したときに角運動量が回転軸と同じ方向を向く特別な方向のことです。数学的には、$\bm{L} = \bm{J}\bm{\omega}$ が $\bm{\omega}$ と平行になる方向、すなわち固有値問題

$$ \bm{J} \bm{e} = \lambda \bm{e} $$

の固有ベクトル $\bm{e}$ が主軸方向であり、対応する固有値 $\lambda$ がその軸まわりの 主慣性モーメント です。

実対称行列 $\bm{J}$ は必ず3つの実固有値 $I_1, I_2, I_3$ と、互いに直交する3つの固有ベクトル $\bm{e}_1, \bm{e}_2, \bm{e}_3$ を持ちます(スペクトル定理)。

対角化

固有ベクトルを列に並べた行列 $\bm{P} = [\bm{e}_1, \bm{e}_2, \bm{e}_3]$ を使うと、

$$ \bm{P}^T \bm{J} \bm{P} = \begin{pmatrix} I_1 & 0 & 0 \\ 0 & I_2 & 0 \\ 0 & 0 & I_3 \end{pmatrix} = \text{diag}(I_1, I_2, I_3) $$

$\bm{P}$ は直交行列($\bm{P}^T = \bm{P}^{-1}$)なので、これは座標系の回転に対応します。つまり、主軸を座標軸に選べば慣性テンソルは対角行列になるのです。

主軸座標系では、角運動量と回転エネルギーが特にシンプルな形になります。

$$ \bm{L} = (I_1 \omega_1, \; I_2 \omega_2, \; I_3 \omega_3)^T $$

$$ T = \frac{1}{2}(I_1 \omega_1^2 + I_2 \omega_2^2 + I_3 \omega_3^2) $$

各軸が独立に寄与する明快な形です。

対称性と主軸

物体の幾何学的対称性は、主軸の方向と主慣性モーメントの関係を決定します。

回転対称体(円柱、円錐、楕円体など): 回転対称軸が1つの主軸になり、対称軸に垂直な平面内の任意の方向が残り2つの主軸になります($I_1 = I_2 \neq I_3$ のような縮退が生じる)。

3つの対称面を持つ物体(直方体など): 対称面の交線が主軸方向になり、計算しなくても主軸がわかります。

非対称体: 一般に主軸は幾何学的に自明ではなく、固有値問題を解く必要があります。

慣性テンソルの非対角成分の意味をもう少し掘り下げておきましょう。

慣性積の物理的意味

慣性積 $J_{xy} = -\int xy \, dm$ がゼロでないとは、どういう状況でしょうか。

直感的に言えば、質量分布が座標平面に対して偏っていることを意味します。例えば、$J_{xy} \neq 0$ は質量が $x > 0, y > 0$ の象限と $x < 0, y < 0$ の象限に偏って分布していることを示します。

慣性積がゼロになる条件

慣性積 $J_{xy}$ がゼロになるための十分条件は、物体が $xz$ 平面($y = 0$)に関して対称であることです。対称性により、$y > 0$ の質量要素と $y < 0$ の対応する質量要素の寄与が打ち消し合い、$\int xy \, dm = 0$ になります。

衛星の設計では、質量分布をできるだけ対称にして慣性積を小さくすることが望まれます。慣性積が大きいと、1つの軸まわりの回転が他の軸にもトルクを生み、制御が複雑になるためです。

動バランスとの関係

工学では、回転する物体の慣性積をゼロにする操作を 動バランシング と呼びます。自動車のタイヤバランスや工作機械のスピンドルバランスがこれに該当します。慣性積がゼロでない回転体は、軸受けに周期的な力(動的アンバランス力)を生じ、振動や寿命低下の原因になります。

慣性テンソルの性質を理解したところで、実用上重要な定理を確認しましょう。重心を通らない軸まわりの慣性モーメントを計算する平行軸の定理です。

平行軸の定理(シュタイナーの定理)

定理の意味

重心(質量中心)を通る軸まわりの慣性モーメント $I_G$ がわかっているとき、それと平行な別の軸まわりの慣性モーメント $I$ は次の公式で計算できます。

$$ I = I_G + m d^2 $$

ここで $m$ は全質量、$d$ は2つの平行な軸の間の距離です。

この定理は「遠くの軸ほど慣性モーメントが大きい」という直感に合致します。重心から離れた軸のまわりに回すには、全質量が余分に $d$ だけ移動する分のエネルギーが必要だからです。

テンソル形式の平行軸の定理

テンソル全体に対する平行軸の定理は、重心を原点とする座標系での慣性テンソル $\bm{J}_G$ と、それから $\bm{d} = (d_x, d_y, d_z)^T$ だけ平行移動した座標系での慣性テンソル $\bm{J}$ の間の関係として書けます。

$$ J_{ij} = (J_G)_{ij} + m\left(|\bm{d}|^2 \delta_{ij} – d_i d_j\right) $$

行列形式で書くと、

$$ \bm{J} = \bm{J}_G + m\left(|\bm{d}|^2 \bm{I} – \bm{d}\bm{d}^T\right) $$

この公式は、複雑な形状の慣性テンソルを計算する際に非常に便利です。全体を単純な部品に分解し、各部品の重心まわりの慣性テンソルを計算してから、平行軸の定理で全体の重心まわりに変換して合算するのです。

導出

平行軸の定理を導出しましょう。新しい座標系の原点が重心から $\bm{d}$ だけずれているとき、質量要素の位置は $\bm{r}’ = \bm{r}_G + \bm{d}$ と書けます。ここで $\bm{r}_G$ は重心からの相対位置です。

$x$ 軸まわりの慣性モーメント $J’_{xx}$ を計算すると、

$$ J’_{xx} = \int (y’^2 + z’^2) \, dm = \int \left[(y_G + d_y)^2 + (z_G + d_z)^2\right] dm $$

展開すると3種類の項が現れます。

$$ = \int (y_G^2 + z_G^2) \, dm + 2d_y \int y_G \, dm + 2d_z \int z_G \, dm + (d_y^2 + d_z^2) \int dm $$

第1項は $J_{G,xx}$(重心まわりの慣性モーメント)、第2項・第3項は重心の定義 $\int \bm{r}_G \, dm = \bm{0}$ よりゼロ、第4項は $m(d_y^2 + d_z^2)$ です。

よって、

$$ J’_{xx} = J_{G,xx} + m(d_y^2 + d_z^2) $$

スカラーの場合に $d^2 = d_y^2 + d_z^2$($x$ 軸からの距離の2乗)なので、$I = I_G + md^2$ が得られます。非対角成分も同様に導出できます。

次に、代表的な形状の慣性モーメントを確認し、衛星設計への応用を見ていきましょう。

代表的な形状の慣性モーメント

球殻

半径 $R$、質量 $M$ の薄い球殻(中空の球)の、中心を通る任意の軸まわりの慣性モーメントは、

$$ I = \frac{2}{3} MR^2 $$

球の対称性から、どの軸でも同じ値になります。慣性テンソルは $\bm{J} = \frac{2}{3}MR^2 \bm{I}$ です。

中実球

半径 $R$、質量 $M$ の一様な中実球の慣性モーメントは、

$$ I = \frac{2}{5} MR^2 $$

球殻より小さいのは、質量が中心付近にも分布しており、回転軸からの平均距離が球殻より小さいためです。

導出を示しましょう。密度 $\rho$ の球を厚さ $dr$ の薄い球殻に分割します。半径 $r$ の球殻の質量は $dm = \rho \cdot 4\pi r^2 \, dr$、その慣性モーメントは $dI = \frac{2}{3}r^2 \, dm = \frac{2}{3}\rho \cdot 4\pi r^4 \, dr$ です。

$r = 0$ から $R$ まで積分すると、

$$ I = \frac{8\pi\rho}{3} \int_0^R r^4 \, dr = \frac{8\pi\rho}{3} \cdot \frac{R^5}{5} = \frac{8\pi\rho R^5}{15} $$

全質量 $M = \frac{4}{3}\pi R^3 \rho$ を代入して $\rho$ を消去すると、

$$ I = \frac{8\pi R^5}{15} \cdot \frac{3M}{4\pi R^3} = \frac{2}{5}MR^2 $$

円柱(回転対称軸まわり)

半径 $R$、質量 $M$ の一様な円柱の、中心軸まわりの慣性モーメントは、

$$ I_{\text{axis}} = \frac{1}{2} MR^2 $$

中心軸に垂直な軸(高さ $h$ の円柱の直径方向)まわりでは、

$$ I_{\perp} = \frac{1}{12} M(3R^2 + h^2) $$

$I_{\perp}$ が $R$ と $h$ の両方に依存するのは、直径方向まわりの回転では半径方向と高さ方向の両方の質量分布が効くためです。

直方体

辺の長さ $a$, $b$, $c$、質量 $M$ の一様な直方体の、重心を通り各辺に平行な軸まわりの慣性モーメントは、

$$ I_x = \frac{M}{12}(b^2 + c^2), \quad I_y = \frac{M}{12}(a^2 + c^2), \quad I_z = \frac{M}{12}(a^2 + b^2) $$

各辺に平行な軸が主軸(対称性から明らか)であり、慣性積はゼロです。

まとめ表

形状 回転軸 慣性モーメント
薄い球殻 中心通過の任意軸 $\frac{2}{3}MR^2$
中実球 中心通過の任意軸 $\frac{2}{5}MR^2$
円柱 中心軸 $\frac{1}{2}MR^2$
円柱 直径方向 $\frac{1}{12}M(3R^2+h^2)$
直方体 $x$ 軸(辺 $a$ 方向) $\frac{1}{12}M(b^2+c^2)$
細い棒 中央垂直 $\frac{1}{12}ML^2$
薄い円板 中心軸 $\frac{1}{2}MR^2$
薄い円板 直径 $\frac{1}{4}MR^2$

代表的な形状の慣性モーメントを確認しました。次に、これらの知識が衛星設計にどのように活かされるかを見ましょう。

衛星の慣性テンソル設計

衛星の慣性テンソルの重要性

人工衛星の慣性テンソルは、姿勢制御系の設計に直結する最も基本的なパラメータです。オイラーの回転方程式(衛星の姿勢のダイナミクスを記述する基本方程式)は、

$$ \bm{J}\dot{\bm{\omega}} + \bm{\omega} \times (\bm{J}\bm{\omega}) = \bm{T}_{\text{ext}} $$

ここで $\bm{J}$ は慣性テンソル、$\bm{\omega}$ は角速度、$\bm{T}_{\text{ext}}$ は外部トルクです。

$\bm{J}$ の値によって、同じトルクに対する角加速度の応答が変わります。具体的には、

  • 大きな慣性モーメント: 姿勢変更が遅いが、外乱に対して安定(角速度変化が小さい)
  • 小さな慣性モーメント: 姿勢変更が速いが、外乱に敏感

主軸設計の原則

衛星の設計では、以下の原則に従って慣性テンソルを設計します。

原則1: 機体座標系を主軸に揃える

衛星の機体座標系が慣性テンソルの主軸と一致するように設計します。これにより慣性積がゼロになり、3軸の姿勢制御が独立に(デカップルして)行えます。

原則2: 重力傾斜安定のための慣性モーメント比

地球指向衛星では、重力傾斜による安定化を利用するために、ナディア方向の慣性モーメントが最大になるように設計します。具体的には、

$$ I_z > I_x, \quad I_z > I_y $$

ここで $z$ 軸はナディア方向です。これにより、重力傾斜トルクが復元力として作用し、パッシブに姿勢が安定化します。

原則3: スピン安定の条件

スピン安定化衛星(スピン軸まわりに回転して安定化する衛星)では、

$$ I_{\text{spin}} > I_{\text{transverse}} $$

スピン軸の慣性モーメントが他の軸より大きい(最大主慣性モーメント軸まわりのスピン)と、内部のエネルギー散逸が安定化に寄与します。逆に最小主慣性モーメント軸まわりのスピンは、エネルギー散逸により不安定になります(フラットスピンへの遷移)。

複合体の慣性テンソル計算

実際の衛星は、複数の部品(本体、太陽電池パネル、アンテナ、燃料タンクなど)から構成されます。全体の慣性テンソルは、各部品の慣性テンソルを平行軸の定理で衛星重心まわりに変換し、足し合わせて求めます。

$$ \bm{J}_{\text{total}} = \sum_{k} \left[\bm{J}_{G,k} + m_k \left(|\bm{d}_k|^2 \bm{I} – \bm{d}_k \bm{d}_k^T\right)\right] $$

ここで $\bm{d}_k$ は部品 $k$ の重心から衛星重心までの位置ベクトルです。

燃料の消費に伴い質量分布が変化するため、慣性テンソルもミッション期間中に変化します。推進剤を多く搭載するミッションでは、初期と末期で慣性テンソルが大きく異なるため、制御ゲインの切り替えやオンライン推定が必要になります。

理論をPythonで実装し、慣性テンソルの性質を視覚的に確認しましょう。

Pythonでの実装と可視化

慣性テンソルの固有値と主軸の計算

import numpy as np

# 非対称な衛星の慣性テンソルの例
# (本体 + 太陽電池パネル + アンテナの寄与を合成した結果)
J = np.array([
    [150.0,  -20.0,  10.0],
    [-20.0,  200.0, -15.0],
    [ 10.0,  -15.0, 180.0]
])  # [kg m^2]

print("慣性テンソル J:")
print(J)
print(f"\n対称性の確認: ||J - J^T|| = {np.linalg.norm(J - J.T):.2e}")

# 固有値分解
eigenvalues, eigenvectors = np.linalg.eigh(J)

print(f"\n主慣性モーメント [kg m^2]:")
for i, val in enumerate(eigenvalues):
    print(f"  I_{i+1} = {val:.4f}")
    print(f"  主軸方向: {eigenvectors[:, i].round(4)}")

# 主軸座標系での慣性テンソル(対角化の確認)
J_principal = eigenvectors.T @ J @ eigenvectors
print(f"\n主軸座標系での慣性テンソル:")
print(np.round(J_principal, 8))

# 角運動量の計算例
omega = np.array([0.1, 0.0, 0.0])  # x軸まわりの回転
L = J @ omega
print(f"\n角速度 ω = {omega}")
print(f"角運動量 L = {L.round(4)}")
print(f"Lとωの角度: {np.degrees(np.arccos(np.dot(L, omega) / (np.linalg.norm(L) * np.linalg.norm(omega)))):.2f}°")

出力から以下のことが確認できます。

  1. 主慣性モーメント: 3つの固有値がすべて正であり、これは物理的に意味のある慣性テンソルの条件です。
  2. 対角化: 主軸座標系での慣性テンソルは対角行列になり、非対角成分はマシンイプシロン程度のゼロです。
  3. 角運動量の方向ずれ: $x$ 軸まわりに回転しても、慣性積が存在するため角運動量は $x$ 軸からずれています。このずれの角度は非対角成分の大きさに依存します。

慣性楕円体の3Dプロット

慣性テンソルの性質を直感的に理解するための強力なツールが 慣性楕円体 です。慣性楕円体は、各方向の慣性モーメントの大きさを3D表面で表したものです。

方向 $\hat{\bm{n}}$ まわりの慣性モーメントは $I(\hat{\bm{n}}) = \hat{\bm{n}}^T \bm{J} \hat{\bm{n}}$ で与えられます。慣性楕円体は、$\bm{x}^T \bm{J} \bm{x} = 1$ で定義される曲面です。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 慣性テンソル
J = np.array([
    [150.0,  -20.0,  10.0],
    [-20.0,  200.0, -15.0],
    [ 10.0,  -15.0, 180.0]
])

eigenvalues, eigenvectors = np.linalg.eigh(J)

# 球面上の点を生成
theta = np.linspace(0, np.pi, 60)
phi = np.linspace(0, 2*np.pi, 80)
theta_grid, phi_grid = np.meshgrid(theta, phi)

# 単位球面上の方向ベクトル
nx = np.sin(theta_grid) * np.cos(phi_grid)
ny = np.sin(theta_grid) * np.sin(phi_grid)
nz = np.cos(theta_grid)

# 各方向の慣性モーメント I(n) = n^T J n
I_n = np.zeros_like(nx)
for i in range(nx.shape[0]):
    for j in range(nx.shape[1]):
        n = np.array([nx[i,j], ny[i,j], nz[i,j]])
        I_n[i,j] = n @ J @ n

# 慣性楕円体: 半径 = 1/sqrt(I) で表示
r = 1.0 / np.sqrt(I_n)
x_ell = r * nx
y_ell = r * ny
z_ell = r * nz

fig = plt.figure(figsize=(16, 6))

# (a) 慣性楕円体
ax1 = fig.add_subplot(131, projection='3d')
surf = ax1.plot_surface(x_ell, y_ell, z_ell, facecolors=plt.cm.viridis(
    (I_n - I_n.min()) / (I_n.max() - I_n.min())), alpha=0.7)
# 主軸を表示
scale = 0.09
for i in range(3):
    ev = eigenvectors[:, i] * scale
    ax1.quiver(0,0,0,*ev, color=['r','g','b'][i], linewidth=2,
               arrow_length_ratio=0.15)
    ax1.quiver(0,0,0,*(-ev), color=['r','g','b'][i], linewidth=2,
               arrow_length_ratio=0)
    ax1.text(*(ev*1.3), f'$I_{i+1}$={eigenvalues[i]:.0f}',
             color=['r','g','b'][i], fontsize=9)

ax1.set_xlabel('X'); ax1.set_ylabel('Y'); ax1.set_zlabel('Z')
ax1.set_title('Inertia Ellipsoid (r = 1/√I)', fontsize=12)

# (b) 方向別の慣性モーメント(ヒートマップ)
ax2 = fig.add_subplot(132)
c = ax2.pcolormesh(np.degrees(phi), np.degrees(theta),
                   I_n.T, shading='auto', cmap='viridis')
plt.colorbar(c, ax=ax2, label=r'$I(\hat{n})$ [kg m²]')
ax2.set_xlabel('Azimuth φ [deg]', fontsize=12)
ax2.set_ylabel('Polar θ [deg]', fontsize=12)
ax2.set_title('Moment of Inertia Map', fontsize=12)

# (c) 主軸と機体軸のずれ
ax3 = fig.add_subplot(133, projection='3d')
# 機体座標系
for i, (c_ax, lbl) in enumerate(zip(['lightcoral','lightgreen','lightblue'],
                                     ['$x_b$','$y_b$','$z_b$'])):
    d = np.zeros(3); d[i] = 1.0
    ax3.quiver(0,0,0,*d, color=c_ax, linewidth=2, arrow_length_ratio=0.08)
    ax3.text(*(d*1.1), lbl, fontsize=11, color=c_ax)

# 主軸
for i, (c_ax, lbl) in enumerate(zip(['r','g','b'],
                                     ['$e_1$','$e_2$','$e_3$'])):
    ev = eigenvectors[:, i]
    ax3.quiver(0,0,0,*ev, color=c_ax, linewidth=3, arrow_length_ratio=0.08,
               linestyle='solid')
    ax3.text(*(ev*1.15), f'{lbl} ({eigenvalues[i]:.0f})',
             fontsize=10, color=c_ax)

ax3.set_xlabel('X'); ax3.set_ylabel('Y'); ax3.set_zlabel('Z')
ax3.set_title('Body Axes (light) vs Principal Axes (bold)', fontsize=11)
ax3.set_xlim([-1.3,1.3]); ax3.set_ylim([-1.3,1.3]); ax3.set_zlim([-1.3,1.3])

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

上の3つのグラフから、慣性テンソルの幾何学的性質が視覚的に理解できます。

  1. 左図(慣性楕円体): 楕円体の形状が慣性モーメントの方向依存性を反映しています。楕円体が膨らんでいる方向は $1/\sqrt{I}$ が大きい = 慣性モーメントが小さい方向で、回転させやすい方向です。赤・緑・青の矢印は3つの主軸を示し、楕円体の対称軸と一致しています。色は慣性モーメントの大きさを表しています(黄色が大きい方向)。
  2. 中図(方向別慣性モーメント): 全方向の慣性モーメントを球面座標でヒートマップにしたものです。最小値と最大値の比から、この衛星の異方性の程度がわかります。
  3. 右図(機体軸と主軸のずれ): 薄い色の軸が機体座標系、太い色の軸が主軸です。慣性積が存在するため、2つの座標系は完全には一致せず、わずかに傾いています。

平行軸の定理の数値検証

複合体の慣性テンソルを平行軸の定理で計算する例を示します。

import numpy as np

def parallel_axis(J_cm, mass, d):
    """平行軸の定理: 重心まわりJ_cmから距離dだけ離れた軸まわりのJを計算"""
    d = np.array(d)
    return J_cm + mass * (np.dot(d, d) * np.eye(3) - np.outer(d, d))

# 衛星本体(直方体: 1m x 1m x 1.5m, 300kg)
m_body = 300.0
a, b, c = 1.0, 1.0, 1.5  # [m]
J_body = m_body / 12 * np.diag([b**2+c**2, a**2+c**2, a**2+b**2])

# 太陽電池パネル1(薄い板: 2m x 1m, 15kg, 本体から+y方向に1.5m)
m_panel = 15.0
wp, hp = 2.0, 1.0  # 幅、高さ(厚さ無視)
J_panel_cm = m_panel / 12 * np.diag([hp**2, wp**2, wp**2+hp**2])
d_panel1 = np.array([0, 1.5, 0])  # 本体重心からの距離
J_panel1 = parallel_axis(J_panel_cm, m_panel, d_panel1)

# 太陽電池パネル2(対称配置: -y方向に1.5m)
d_panel2 = np.array([0, -1.5, 0])
J_panel2 = parallel_axis(J_panel_cm, m_panel, d_panel2)

# アンテナ(球殻: 半径0.3m, 5kg, +z方向に1.0m)
m_ant = 5.0
R_ant = 0.3
J_ant_cm = 2/3 * m_ant * R_ant**2 * np.eye(3)
d_ant = np.array([0, 0, 1.0])
J_ant = parallel_axis(J_ant_cm, m_ant, d_ant)

# 全体の慣性テンソル
J_total = J_body + J_panel1 + J_panel2 + J_ant

print("=== 各部品の慣性テンソル(重心まわり) ===")
print(f"本体:         diag = {np.diag(J_body).round(2)}")
print(f"パネル(CM):   diag = {np.diag(J_panel_cm).round(2)}")
print(f"アンテナ(CM): diag = {np.diag(J_ant_cm).round(4)}")

print(f"\n=== 全体の慣性テンソル ===")
print(np.round(J_total, 4))

# 主慣性モーメント
eigvals, eigvecs = np.linalg.eigh(J_total)
print(f"\n主慣性モーメント:")
for i in range(3):
    print(f"  I_{i+1} = {eigvals[i]:.2f} kg m², 主軸: {eigvecs[:,i].round(4)}")

# 質量合計
m_total = m_body + 2*m_panel + m_ant
print(f"\n全質量: {m_total} kg")

# 対称性の確認
print(f"\n対称配置のためJxy, Jyz ≈ 0の確認:")
print(f"  J_xy = {J_total[0,1]:.6f}")
print(f"  J_yz = {J_total[1,2]:.6f}")
print(f"  J_xz = {J_total[0,2]:.6f}")

出力から、衛星の複合慣性テンソルの特性が確認できます。

  1. 対称配置の効果: 太陽電池パネルを $\pm y$ 方向に対称配置しているため、$J_{xy}$ と $J_{yz}$ はほぼゼロです。これにより慣性テンソルはほぼ対角になり、姿勢制御の設計が簡単になります。
  2. 太陽電池パネルの寄与: パネルは質量が小さい(15kg×2)ですが、重心から1.5m離れているため、平行軸の定理の $md^2$ の項により $y$ 軸周り以外の慣性モーメントに大きく寄与します。
  3. 慣性モーメントの比: 衛星設計として重力傾斜安定の条件(ナディア軸の慣性モーメントが最大)を満たしているかどうかを確認できます。

代表的形状の慣性モーメントの比較可視化

import numpy as np
import matplotlib.pyplot as plt

# 各形状の慣性モーメント係数 I = k * M * R^2(または L^2)
shapes = {
    'Solid sphere': 2/5,
    'Thin spherical shell': 2/3,
    'Solid cylinder (axis)': 1/2,
    'Thin disk (axis)': 1/2,
    'Thin disk (diameter)': 1/4,
    'Thin ring': 1,
    'Thin rod (center)': 1/12,
    'Thin rod (end)': 1/3,
}

fig, ax = plt.subplots(figsize=(10, 6))
names = list(shapes.keys())
coeffs = list(shapes.values())
colors = plt.cm.viridis(np.linspace(0.2, 0.8, len(names)))

bars = ax.barh(range(len(names)), coeffs, color=colors, edgecolor='black', linewidth=0.5)
ax.set_yticks(range(len(names)))
ax.set_yticklabels(names, fontsize=11)
ax.set_xlabel(r'Coefficient $k$ in $I = k \cdot MR^2$', fontsize=13)
ax.set_title('Moment of Inertia Coefficients for Standard Shapes', fontsize=14)

# 数値をバーの上に表示
for bar, coeff in zip(bars, coeffs):
    ax.text(bar.get_width() + 0.01, bar.get_y() + bar.get_height()/2,
            f'{coeff:.4f}', va='center', fontsize=10)

ax.set_xlim([0, 1.15])
ax.grid(True, axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig('inertia_coefficients.png', dpi=150, bbox_inches='tight')
plt.show()

上のグラフから、形状による慣性モーメントの違いが一目でわかります。

  1. 薄いリング($k = 1$)が最大です。全質量が回転軸から最も遠い位置にあるためです。
  2. 中実球($k = 2/5$)は球殻($k = 2/3$)より小さいです。内部に質量が詰まっているため、平均的な回転軸からの距離が小さくなります。
  3. 棒の端を軸にした場合($k = 1/3$)は中央を軸にした場合($k = 1/12$)の4倍です。これは平行軸の定理から説明できます: $I_{\text{end}} = I_{\text{center}} + M(L/2)^2 = ML^2/12 + ML^2/4 = ML^2/3$。

慣性テンソルの回転変換

最後に、座標系を回転させたときに慣性テンソルがどう変換されるかを確認します。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 主軸座標系での対角慣性テンソル
I1, I2, I3 = 100.0, 200.0, 150.0
J_principal = np.diag([I1, I2, I3])

# 30度回転させた座標系での慣性テンソル
angle = np.radians(30)
# z軸まわり30度回転
Rz = np.array([
    [np.cos(angle), -np.sin(angle), 0],
    [np.sin(angle),  np.cos(angle), 0],
    [0, 0, 1]
])

# テンソルの回転変換: J' = R J R^T
J_rotated = Rz @ J_principal @ Rz.T

print("主軸座標系での慣性テンソル(対角):")
print(np.round(J_principal, 2))
print(f"\n30° 回転した座標系での慣性テンソル:")
print(np.round(J_rotated, 4))

# 固有値の確認(回転で不変)
eigvals_orig = np.linalg.eigvalsh(J_principal)
eigvals_rot = np.linalg.eigvalsh(J_rotated)
print(f"\n元の固有値:     {np.sort(eigvals_orig).round(4)}")
print(f"回転後の固有値: {np.sort(eigvals_rot).round(4)}")
print(f"固有値は座標系の回転に対して不変です")

# トレースの確認(回転で不変)
print(f"\ntr(J_principal) = {np.trace(J_principal):.4f}")
print(f"tr(J_rotated)   = {np.trace(J_rotated):.4f}")

出力から、テンソル変換の重要な性質が確認できます。

  1. 座標変換の公式: 回転行列 $R$ による座標変換は $J’ = R J R^T$ で行われます。回転した座標系では慣性積(非対角成分)が現れます。
  2. 固有値の不変性: 固有値(主慣性モーメント)は座標系の選び方によらず一定です。これは物理量として当然の性質です。
  3. トレースの不変性: $\text{tr}(\bm{J}) = I_1 + I_2 + I_3$ も座標系に依存しません。トレースの不変性は、$I_{xx} + I_{yy} + I_{zz} = 2\int |\bm{r}|^2 dm$ が座標系によらないことから理解できます。

まとめ

本記事では、回転運動の核となる物理量である慣性モーメントと慣性テンソルについて解説しました。

  • 慣性モーメント は回転に対する抵抗の大きさで、$I = \int r_\perp^2 \, dm$ で定義される。回転軸から遠い質量ほど大きく寄与する
  • 慣性テンソル $\bm{J}$ は3×3の実対称行列で、任意の方向の慣性モーメントを統一的に記述する。$\bm{L} = \bm{J}\bm{\omega}$ で角運動量と角速度を結ぶ
  • 主軸と対角化: 慣性テンソルの固有ベクトルが主軸、固有値が主慣性モーメント。主軸座標系では慣性テンソルが対角行列になる
  • 慣性積: 質量分布の非対称性を反映する非対角成分。ゼロでなければ角運動量と角速度の方向がずれる
  • 平行軸の定理: $\bm{J} = \bm{J}_G + m(|\bm{d}|^2\bm{I} – \bm{d}\bm{d}^T)$ で、重心以外の軸まわりの慣性テンソルを計算
  • 衛星設計: 慣性テンソルの主軸を機体座標系に揃え、重力傾斜安定の条件を満たすよう設計する

慣性テンソルは、オイラーの回転方程式 $\bm{J}\dot{\bm{\omega}} + \bm{\omega} \times (\bm{J}\bm{\omega}) = \bm{T}$ を通じて衛星の姿勢ダイナミクスを支配します。次のステップでは、この運動方程式を使った姿勢運動の解析(自由回転、トルクフリー運動)や、具体的な姿勢制御則の設計に進みます。

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