制御エンジニアのための剛体の回転運動

制御の教科書には、振り子、モーター、自転車、自動車などさまざまな回転する機械が、解析手法や制御手法を説明するための例題として取り上げられ、数式を使ってモデリングされる。そこではしばしば、運動方程式の回転バージョン

I\frac{d\omega}{dt}={N}
\tag{0}

が登場する。$I$は物体の慣性モーメント、$\omega$は角速度（回転数）、$N$はトルクである。これは、電子工学や化学工学や生物学から制御の世界に足を踏み入れた学習者にとっては曲者である。われわれは、運動方程式との類似から、なんとなく、トルクによって回転が速くなることや、慣性モーメントは回転数の変化しにくさを表していることを推し量ることができる。しかし、そもそもこの式がどこから出てきたのかとか、慣性モーメントの具体的な計算とかいったことになるとおぼつかないし、回転軸が固定されていない3次元の運動となるともはやお手上げである。そこで本稿では、なじみ深い運動方程式からスタートして、3次元の回転運動の方程式と慣性モーメントの3次元バージョンである慣性モーメントテンソルを導出する。

質点の回転運動

剛体の話をするまえに、まず質点の話から始める。静止した座標系（慣性座標系）から見た質量$m$の質点$P$の運動方程式は

m\frac{d^2\boldsymbol{r}}{dt^2} = \boldsymbol{f} \tag{1}

である。$\boldsymbol{f}$は質点$P$に働く力、$\boldsymbol{r}$は質点$P$の位置ベクトルである。時間微分$d/dt$について、ここでは慣性座標系を採用しているために、特に注意する必要はなく、通常通り3つの成分をそれぞれ微分してやればよい（これは、後の節で物体とともに回転する座標系を採用する場合に重要となる）。式(1)の左から$\boldsymbol{r}$を外積でかけると

m\boldsymbol{r} \times \frac{d^2\boldsymbol{r}}{dt^2} 
=
\boldsymbol{r} \times \boldsymbol{f} \tag{2}

となる。ここで同じベクトルどうしの外積はゼロベクトルであることを利用すると

\boldsymbol{r} \times \frac{d^2\boldsymbol{r}}{dt^2} 
= 
\boldsymbol{r} \times \frac{d^2\boldsymbol{r}}{dt^2}
+
\frac{d\boldsymbol{r}}{dt} \times \frac{d\boldsymbol{r}}{dt}
=
\frac{d}{dt}\left(
    \boldsymbol{r} \times \frac{d\boldsymbol{r}}{dt}
\right)

と変形できる（ベクトルの外積についても積の微分公式が利用できる）。これを式(2)に代入すると

\frac{d}{dt}\left(
    m\boldsymbol{r} \times \frac{d\boldsymbol{r}}{dt}
\right)
=
\boldsymbol{r} \times \boldsymbol{f} 
\tag{3}

となる。左辺の時間微分の中身が角運動量ベクトル$\boldsymbol{L}$と定義される量である。回転の勢いを表す量といわれることもあるが、後の節で角運動量ベクトルの向きと回転（角速度ベクトルの向き）が異なるという事実が出てきたときに混乱を招くことになるので、ここでは形式的に受け入れることにする。右辺は力のモーメント、もしくはトルクと定義され$\boldsymbol{N}$と表記される。これは明らかに、座標原点に対して質点$P$を回転させようとする作用を表している。これらの記号を使って式(3)を書き直すと

\frac{d\boldsymbol{L}}{dt}
=
\boldsymbol{N} 
\tag{4}

となる。これが質点の回転運動の方程式である。

大きさのある物体の回転運動

つぎに、大きさのある物体の回転について考える。大きさのある物体は、多数の質点$P_1, \dots,P_n$を質量の無視できる軽い棒でつないだものとして近似できる。ただし、ここではまだこれらの棒は絶対的な硬さを持っている必要はない。言い換えると、この物体は必ずしも剛体でなくてもよいとする。この節の結論を先に書くと、この物体の回転運動の方程式は、質点の方程式(3)の総和を取ったもの

\frac{d}{dt}\left(
    \sum_{i=1}^n m_i\boldsymbol{r}_i \times \frac{d\boldsymbol{r}_i}{dt}
\right)
=
\sum_{i=1}^n \boldsymbol{r}_i \times \boldsymbol{f}_i 
\tag{5}

ということになる。$m_i$と$\boldsymbol{r}_i$はそれぞれ質点$P_i$の質量および位置ベクトルである。$\boldsymbol{f}_i$は質点$P_i$に物体の外から働く力である。左辺の微分の中身と右辺をそれぞれ物体の全角運動量$\boldsymbol{L}$および物体にかかる全トルク$\boldsymbol{N}$として改めて定義すると

\frac{d\boldsymbol{L}}{dt}
=
\boldsymbol{N}
\tag{5'}

と書ける。結果的に質点のときと同じ方程式が得られる。当たり前のように思えるが、注意すべきは、物体を構成する各質点は、外力のほかに、周囲の他の質点からも力（内力）も受けるという点である。以下では、なぜ最終的にこの内力が方程式(5)の右辺から消えたのかを説明する。

周囲の質点から受ける力を含めて、物体を構成する質点$P_i$について運動方程式を立てると

m_i\frac{d^2\boldsymbol{r}_i}{dt^2} = \boldsymbol{f}_i + \sum_{j=1}^n \boldsymbol{F}_{ij}
\tag{6}

である。$\boldsymbol{F}_{ij}$は質点$P_i$が他の質点$P_j$から受ける力であり$\boldsymbol{F} _{ii} = \boldsymbol{0}$である。また、作用反作用の法則より$\boldsymbol{F} _{ij} = - \boldsymbol{F} _{ji}$である。先ほどと同じように、運動方程式(6)の両辺に左から$\boldsymbol{r}_i$との外積をとり、すべての質点について総和を取ると

\frac{d}{dt}\left(
    \sum_{i=1}^n m_i\boldsymbol{r}_i \times \frac{d\boldsymbol{r}_i}{dt}
\right)
=
\sum_{i=1}^n \boldsymbol{r}_i \times \boldsymbol{f}_i 
+
\sum_{i=1}^n \sum_{j=1}^n \boldsymbol{r}_i \times \boldsymbol{F} _{ij}
\tag{7}

となる。問題の右辺第2項について

\sum_{i=1}^n \sum_{j=1}^n \boldsymbol{r}_i \times \boldsymbol{F} _{ij}
=
\sum_{i=1}^n \sum_{j=1}^n \boldsymbol{r}_j \times \boldsymbol{F} _{ji}
=
\sum_{i=1}^n \sum_{j=1}^n \boldsymbol{r}_j \times (-\boldsymbol{F} _{ij})

である。この式の左辺を左辺自身および右辺に足すと

2\sum_{i=1}^n \sum_{j=1}^n \boldsymbol{r}_i \times \boldsymbol{F} _{ij}
=
\sum_{i=1}^n \sum_{j=1}^n (\boldsymbol{r}_i \times \boldsymbol{F} _{ij}
-
\boldsymbol{r}_j \times \boldsymbol{F} _{ij})
=
\sum_{i=1}^n \sum_{j=1}^n (\boldsymbol{r}_i - \boldsymbol{r}_j) \times \boldsymbol{F} _{ij}

である。$\boldsymbol{r}_i - \boldsymbol{r}_j$は質点$P_j$からみた質点$P_i$の方向であるが、これは質点$P_j$から質点$P_i$に働く力$\boldsymbol{F} _{ij}$の方向と同じはずである。同じ方向を向いたベクトルの外積はゼロベクトルとなるので$(\boldsymbol{r}_i - \boldsymbol{r}_j) \times \boldsymbol{F} _{ij} = \boldsymbol{0}$である。よって、式(7)の右辺第2項が消え、大きさのある物体の回転運動の方程式(5)が得られる。まだかなり抽象的な感じがするが、これに剛体（形の変わらない物体）という仮定を加えると、直観的な意味をつかめるようになる。

軸が固定された剛体の回転運動

すでに述べたとおり、方程式(5)は一般の物体について成り立つ。つぎに、特殊なケースとして、$z$軸まわりに回転数$\omega$で回転している剛体を考える。ただし、今回も座標系は慣性座標系（静止している）であり、剛体とともに回転しているわけではないことに注意する。まず、剛体を構成する質点のうち、ある時刻に位置$\boldsymbol{r}_i := (x_i\ y_i\ z_i)^T$に存在する質点$P_i$について考える。回転軸（$z$軸）からの距離を$r_i$とおくと、$x$座標と$y$座標は

\boldsymbol{r}_i
=
\begin{bmatrix}
r_i\cos\theta_i \\ 
r_i\sin\theta_i  \\
z_i
\end{bmatrix}

と書き直せる。$\theta_i$は$zx$平面からの回転角度である。回転軸からの距離$r_i$は一定なので、速度$d\boldsymbol{r}/dt$は

\frac{d\boldsymbol{r}_i}{dt}
=
\begin{bmatrix}
-\frac{d\theta_i}{dt}r_i\sin\theta_i \\ 
\frac{d\theta_i}{dt}r_i\cos\theta_i  \\
0
\end{bmatrix}
=
\begin{bmatrix}
-\frac{d\theta_i}{dt}y_i \\ 
\frac{d\theta_i}{dt}x_i  \\
0
\end{bmatrix}

である。ところで、剛体であるということは、すべての質点が同じ回転数で回転するということである。つまり位置によらず$d\theta_i/dt = \omega$であり

\frac{d\boldsymbol{r}_i}{dt}
=
\begin{bmatrix}
-\omega y_i \\ 
\omega x_i  \\
0
\end{bmatrix}

である。よって質点$P_i$の角運動量ベクトル$\boldsymbol{L}_i$は

\boldsymbol{L}_i
:=
m_i \boldsymbol{r} \times \frac{d\boldsymbol{r}_i}{dt}
=
m_i
\begin{bmatrix}
x_i \\ 
y_i  \\
z_i
\end{bmatrix}
\times
\begin{bmatrix}
-\omega y_i \\ 
\omega x_i  \\
0
\end{bmatrix}
=
m_i
\begin{bmatrix}
-x_i z_i \\ 
-y_i z_i \\
r_i^2
\end{bmatrix}
\omega

である。これをすべての質点について足し合わせた剛体の全角運動量ベクトル$\boldsymbol{L}$は

\boldsymbol{L}
:=
\sum_{i=1}^n \boldsymbol{L}_i
=
\begin{bmatrix}
-\sum_{i=1}^n m_i x_i z_i \\ 
-\sum_{i=1}^n m_i y_i z_i \\
\sum_{i=1}^n m_i r_i^2
\end{bmatrix}
\omega

となる。これを回転運動の方程式(5')に代入すると

\begin{bmatrix}
-\sum_{i=1}^n m_i x_i z_i \\ 
-\sum_{i=1}^n m_i y_i z_i \\
\sum_{i=1}^n m_i r_i^2
\end{bmatrix}
\frac{d\omega}{dt}
+
\begin{bmatrix}
-\frac{d}{dt}\sum_{i=1}^n m_i x_i z_i \\ 
-\frac{d}{dt}\sum_{i=1}^n m_i y_i z_i \\
0
\end{bmatrix}
\omega
=
\boldsymbol{N}
\tag{8}

となる。第2項の$z$成分がゼロであるのは、剛体が回転しても、各点の中心軸からの距離$r_i$は変化しないためである。この方程式の意味を考えてみる。まず、この剛体が一定の回転数で回っているとき、$d\omega/dt = 0$より左辺第1項はゼロになる。よって、たとえ回転数が一定であっても、$x$軸と$y$軸まわりにはトルクが生じるということになる。実際に、いびつな物体をある軸のまわりに一定の回転数で回すことを想像してみると、軸をベアリングでしっかり固定しておかないと、軸がガタガタと動いてしまうだろう。つまり、この方程式の$x$成分と$y$成分は軸を固定して回転を保つために必要なトルクを表していると捉えることができる。残りの$z$成分を取り出すと

\left(
    \sum_{i=1}^n m_i r_i^2
\right)
\frac{d\omega}{dt}
=
N_z \\
\tag{8'}

である。$N_z$は$\boldsymbol{N}$の$z$成分である。スカラー定数$\sum_{i=1}^n m_i r_i^2$は固定軸まわりの慣性モーメントと定義され、$I$と表記される。これを使って式(8')を書き直すと、最初に挙げた固定軸まわりの回転運動の方程式

I_z \frac{d\omega}{dt}
=
N_z \\
\tag{9}

が得られる。最初に書いたとおり、$I_z$は回転の変化のしにくさを表しているが、その定義から、質量が同じ物体を比べると、軸から遠い部分に質量が集まっているほど回転数が変化しにくい、つまり回りにくく止まりにくい物体であることがわかる。

円盤の回転運動：慣性モーメントの計算例

前節で導いた固定軸まわりの回転運動の方程式(8)を使って、均質な素材でできた円盤の回転について考えてみる。なぜなら、われわれは経験からこの円盤が中心軸まわりになめらかに回ると知っているからである。円盤の半径を$l$、厚みを$2h$、質量を$M$とする。円盤の中心軸を$z$軸とし、厚みの真ん中に原点があるとする。まず方程式(8)の$x$成分に現れる総和を計算する。円盤の密度を$\rho=M/(2\pi l^2 h)$とおくと

\sum_{i=1}^n m_i x_i z_i 
\approx
\int xz\ dm
=
\rho\int \int \int xz\ dxdydz

である。円盤は$yz$平面に関してつねに対称であり、被積分関数$xz$は$x$に関して奇関数である。したがって$yz$平面の両側で被積分関数が打ち消し合うことになり、$\sum_{i=1}^n m_i x_i z_i = 0$である。同様に$\sum_{i=1}^n m_i y_i z_i = 0$である。 固定軸まわりの慣性モーメント$I_z$は、円盤を幅$dr$、厚さ$2h$の径の異なるたくさんのリングに分割することで求めることができる。これらのリングのうち、半径が$r$のリングの質量$dm$は、リングの周長が$2\pi r$であることから

dm = \rho 2\pi r 2h dr = \frac{2Mrdr}{l^2}

である。よって

I_z := \sum_{i=1}^n m_i r_i^2
\approx
\int r^2 dm
=
\frac{2M}{l^2}\int_0^l r^3dr = \frac{Ml^2}{2} 

と求まる。以上より回転運動の方程式は

\begin{bmatrix}
0 \\ 
0 \\
Ml^2/2
\end{bmatrix}
\frac{d\omega}{dt}
=
\boldsymbol{N}

となる。この式からはふたつのことがわかる。まず、理想的な円盤は経験的にも明らかなように、極めてなめらかに回転するため、軸を保持するための力（トルク）が不要であるということ。もうひとつは、円盤の中心軸まわりの慣性モーメントは、質量が同じであれば厚みには依存しないということである。

軸が固定されていない剛体の回転運動

剛体の運動は、並進運動と回転運動に分解することができる。並進運動とは、剛体を構成する各質点が、固定点と同じ方向に平行に移動する運動である。一方、回転運動とは、剛体が固定点に対して回転する運動である。ここでの固定点とは、剛体に対する相対的な位置が固定されている点という意味である。並進運動を表すうえで便利であるために、重心を固定点として定義することが多いが、実は固定点は必ずしも重心でなくてもよい。固定点まわりの回転というと複雑な感じがしてイメージがわきにくいが、瞬間ごとにみると、固定点を通るある軸のまわりの回転にすぎない。これまでの軸が固定された回転と違うのは、この軸の方向が時々刻々と変化していることである。

角速度ベクトル

まず、軸の向きが変化する回転をあらわすために角速度ベクトル$\boldsymbol{\omega}$というものを導入する。角速度ベクトルは、その向きが回転軸の方向（回転に対して右ねじの向き）として、その大きさが回転数として定義されるベクトルである。例えば角速度ベクトル

\boldsymbol{\omega}
=
\begin{bmatrix}
1 \\ 
1 \\
0
\end{bmatrix}

は$xy$平面上の直線$y = x$を軸とする回転数$\sqrt{2}$の回転を表す。このように、軸の向きと回転の速さをまとめて表すことができるという点で角速度ベクトルは便利であるが、角速度ベクトルを足し合わせるということは、それらが表す回転運動を足し合わせる、つまり同時に行うことと同じだろうか？

軸の方向が同じふたつの回転運動を同時に行うと、軸の方向は変わらず、回転数が足し合わされるというのは理解しやすい。問題は軸の方向が異なる場合である。物体を異なるふたつの軸のまわりに同時に回転させると、結果として現れる新しい回転の軸や回転数はどうなるのか、直観的にはやや理解しがたい。そこでこの項では、回転の合成と、それらを表す角速度ベクトルの合成が整合していることを確認する。方針としては、一般の回転について直接確認するのは難しいので、$x,y,z$軸方向の3つの回転の合成について確かめることにする。任意の方向の回転どうしを合成する場合は、それらをいったん$x,y,z$軸方向の3つの回転に分解したあと、同じ軸どうしで個別に合成すればよい。

いま、原点を中心とする半径$1$の球がある。球の$z=0$における円周を「赤道」とよぶことにする。まず、球を$x$軸まわりに回転数$\omega_x > 0$で軸に対して右ねじの方向に回転させるとする。このとき、瞬間的には赤道上のすべての点は$z$方向にのみ動く。赤道上のある点の位置を$x$軸からの角度$\psi$で表すと、この点が$z$方向に動く速度$v_x(\psi)$は

v_x(\psi) = \omega_x \sin \psi

である。今度は、球を$y$軸まわりに回転数$\omega_y > 0$で回した場合を考える。$x$軸まわりの回転と同じく、赤道上のすべての点は$z$方向にのみ動く。$x$軸から角度$\psi$にある点の速度$v_y(\psi)$は

v_y(\psi) = - \omega_y \cos \psi

である。これら$x$軸まわりと$y$軸まわりのふたつの回転運動が同時に生じると、赤道上では場所によってこれらの回転による$z$方向速度が強め合ったり打ち消し合ったりすることになる。速度がちょうど打ち消し合ってゼロとなる点では

\omega_x \sin \psi - \omega_y \cos \psi = 0

すなわち$\tan \psi = \omega_y/\omega_x$を満たす。$\tan\psi$は周期$\pi$の周期関数なので、赤道上に速度ゼロとなる点は原点対称な位置にふたつある（それぞれ角度を$\psi_0, \psi_0 + \pi$とおく）。これら2点を結んだ直線こそがふたつの回転が同時に生じてできた新しい回転の軸である（$xy$合成回転軸とよぶことにする）。つまり、新しい回転軸は$x$軸を$y$軸に向かって原点まわりに$\psi_0 = \tan^{-1}(\omega_y/\omega_x)$だけ回した直線である。つぎに回転数を求める。赤道上で$z$方向速度が最大となるのは、この新しい回転軸からもっとも遠い点、すなわち$\psi_0 + \pi/2$と$\psi_0 + 3\pi/2$である。よって、この最大速度は

v_{x+y,\rm{max}}
= 
\omega_x \sin \left(
    \psi_0 + \frac{\pi}{2}
\right)
-
\omega_y \cos \left(
    \psi_0 + \frac{\pi}{2}
\right) \\
=
\omega_x \cos \left(
    \psi_0
\right)
+
\omega_y \sin \left(
    \psi_0
\right) \\

である。両辺を2乗して$\omega_x \sin \psi_0 - \omega_y \cos \psi_0 = 0$を利用すると

v_{x+y,\rm{max}}^2
= 
\omega_x^2
+
\omega_y^2

と変形できる。球の半径は$1$なので、この最大速度$v_{x+y,\rm{max}}$は、新しい回転数$\omega_{x+y}$に等しい。したがって

\omega_{x+y}
= 
\sqrt{
\omega_x^2
+
\omega_y^2
}

である。さて、これら$x$軸まわりと$y$軸まわりの回転に加え、さらに同時に$z$軸まわりにも回転数$\omega_z>0$で球を回すことを考える。$xy$合成回転軸は$z$軸と直交するので、これらふたつの回転軸で張られる平面上に新しい赤道を定めることで、ここまでの議論をそのまま繰り返すことができる。その結果、3つの軸まわりの回転が同時に生じると、回転数は

\omega_{x+y+z}
= 
\sqrt{
\omega_x^2
+
\omega_y^2
+
\omega_z^2
}

となり、回転軸の向きは$xy$合成回転軸から$z$軸に向かって原点まわりにさらに$\phi_0 = \tan^{-1}(\omega_z/\sqrt{\omega_x^2 + \omega_y^2})$だけ回したものになることがわかる。少しわかりにくいが、要するにこれは原点と点$(\omega_x,\omega_y,\omega_z)$を通る直線である。これらはまさに角速度ベクトル

\omega_x \boldsymbol{i}
+
\omega_y \boldsymbol{j}
+
\omega_z \boldsymbol{k}
=:
\begin{bmatrix}
\omega_x \\
\omega_y \\
\omega_z
\end{bmatrix}

の絶対値および方向と同じである（ただし$\boldsymbol{i},\ \boldsymbol{j},\ \boldsymbol{k}$は単位ベクトル）。以上より、回転の合成と、それらを表す角速度ベクトルの合成が整合していることが確認できた。

慣性モーメントテンソル

角速度ベクトルの定義に納得がいったところで、これを使って角運動量ベクトルを表してみる。いま、慣性座標系の原点を固定点として剛体が回転している。ある瞬間の剛体の角速度ベクトルを$\boldsymbol{w}$とする。このとき、剛体のすべての質点は、$\boldsymbol{w}$で表される軸のまわりを、回転数$|\boldsymbol{w}|$で回転している。位置$\boldsymbol{r}_i$にある質点$P_i$と回転軸の距離$h_i$は、$\boldsymbol{r}_i$と回転軸のなす角を$\theta_i$とすると$h_i = |\boldsymbol{r}_i|\sin\theta_i$である。よって、質点$P_i$の速度の大きさ$|d\boldsymbol{r}_i/dt|$は

\left|
    \frac{d\boldsymbol{r}_i}{dt}
\right|
= 
|\boldsymbol{w}| h_i
=
|\boldsymbol{w}||\boldsymbol{r}_i|\sin\theta_i

である。また速度の向きは、$\boldsymbol{r}_i$と回転軸$\boldsymbol{w}$のどちらに対しても垂直である。これはベクトルの外積の定義そのものである。したがって質点$P_i$の速度は外積を使って

\frac{d\boldsymbol{r}_i}{dt}
=
\boldsymbol{w} \times \boldsymbol{r}_i
\tag{10}

と表すことができる。ところで、上記の考え方は位置ベクトルに限らずどんなベクトルの回転にも適用できる。よって

\frac{d\boldsymbol{a}}{dt}
=
\boldsymbol{w} \times \boldsymbol{a}
\tag{10'}

である。さて、式(10)を角運動量ベクトルの定義式に代入すると

\boldsymbol{L}_i
:=
m_i \boldsymbol{r}_i \times \frac{d\boldsymbol{r}_i}{dt}
=
m_i \boldsymbol{r}_i \times (\boldsymbol{w} \times \boldsymbol{r}_i)
\tag{11}

である（質点$P_i$の質量を$m_i$とした）。$\boldsymbol{r}_i \times (\boldsymbol{w} \times \boldsymbol{r}_i)$に対してベクトル3重積の公式を用いると

\boldsymbol{L}_i
=
m_i(
    (\boldsymbol{r}_i^T \boldsymbol{r}_i)\boldsymbol{w}
    -
    (\boldsymbol{r}_i^T\boldsymbol{w})\boldsymbol{r}_i
)
\tag{12}

と変形できる。ここで$\boldsymbol{r}_i := (x_i\ y_i\ z_i)^T,\ \boldsymbol{\omega} := (\omega_x\ \omega_y\ \omega_z)^T$とおいて、右辺第2項の$(\boldsymbol{r}_i^T\boldsymbol{w})\boldsymbol{r}_i$を成分で書き下すと

(\boldsymbol{r}_i^T\boldsymbol{w})\boldsymbol{r}_i
=
(
x_i \omega_x 
+
y_i \omega_y 
+
z_i \omega_z
)
\begin{bmatrix}
x_i \\
y_i \\
z_i
\end{bmatrix}
=
\begin{bmatrix}
x_i^2 \omega_x + x_i y_i \omega_y + z_i x_i \omega_z\\
x_i y_i \omega_x + y_i^2 \omega_y + y_i z_i \omega_z\\
z_i x_i \omega_x + y_i z_i \omega_y + z_i^2 \omega_z
\end{bmatrix}
=
\begin{bmatrix}
x_i^2 & x_i y_i & z_i x_i\\
x_i y_i & y_i^2 & y_i z_i \\
z_i x_i & y_i z_i & z_i^2
\end{bmatrix}
\boldsymbol{w}

これを、$(\boldsymbol{r}_i^T \boldsymbol{r}_i)\boldsymbol{w} = (\boldsymbol{r}_i^T \boldsymbol{r}_i) E \boldsymbol{w}$と書けることに注意して（ただし$E$は単位行列）、式(12)に代入すると

\boldsymbol{L}_i
=
m_i
\begin{bmatrix}
y_i^2 + z_i^2 & -x_i y_i & -z_i x_i\\
-x_i y_i & z_i^2 + x_i^2 & -y_i z_i \\
-z_i x_i & -y_i z_i & x_i^2 + y_i^2
\end{bmatrix}
\boldsymbol{w}
\tag{13}

となる。よって剛体の全角運動量ベクトル$\boldsymbol{L}$は$\boldsymbol{L}_i$の総和を取って

\boldsymbol{L}
=
\begin{bmatrix}
\sum_i m_i (y_i^2 + z_i^2) & -\sum_i m_i x_i y_i & -\sum_i m_i z_i x_i\\
-\sum_i m_i x_i y_i & \sum_i m_i (z_i^2 + x_i^2) & -\sum_i m_i y_i z_i \\
-\sum_i m_i z_i x_i & -\sum_i m_i y_i z_i & \sum_i m_i (x_i^2 + y_i^2)
\end{bmatrix}
\boldsymbol{w}
\tag{14}

となる。右辺に現れた対称行列こそが3次元バージョンの慣性モーメント、すなわち慣性モーメントテンソル$I$である。また、$I$の非対角成分を慣性乗積とよぶ。この式をみると、$I$が単位行列のスカラー倍である特殊な場合を除き、角運動量ベクトルの方向と角速度ベクトルの方向は一致しないことがわかる。これが、角運動量ベクトルを定義したときに、これを「回転の勢いを表す量」として捉えることを避けた理由である。

回転する座標系におけるベクトルの時間微分

あとは、$\boldsymbol{L}$の時間微分をとれば3次元バージョンの剛体の回転運動の方程式の出来上がりだが、ここでひとつの問題が生じる。ここまでずっと慣性座標系（静止している座標系）を採用してきたわけだが、慣性座標系では剛体の座標が時間とともに変わっていくので、結果として慣性モーメントテンソル$I$も変わってしまう。これは剛体の運動を考えるうえで非常に面倒くさい。ちなみに、前の節で求めた固定軸まわりの慣性モーメントは、剛体の各点と固定軸との距離にしか依存しないため、剛体の向きによらず一定であった。

この問題を回避するため、剛体とともに原点まわりに回転する座標系を使って、回転運動を考えなおすことにする。剛体とともに回転する座標系では、当然ながら剛体の各点の座標は変化しないので、慣性モーメントテンソルもやはり不変である。ただし、回転する座標系で物体の運動を考える際には時間微分に気を付ける必要がある。

静止した座標系では、ベクトルそのものの時間微分と、ベクトルの成分の時間微分は等しい。つまり、あるベクトル$\boldsymbol{a}$が単位ベクトルを使って$\boldsymbol{i},\boldsymbol{j},\boldsymbol{k}$

\boldsymbol{a}
=
x\boldsymbol{i}
+
y\boldsymbol{j}
+
z\boldsymbol{k}

と表せるとき

\frac{d\boldsymbol{a}}{dt}
=
\frac{dx}{dt}\boldsymbol{i}
+
\frac{dy}{dt}\boldsymbol{i}
+
\frac{dz}{dt}\boldsymbol{k}

である。ベクトルそのものの時間微分$d\boldsymbol{a}/dt$と区別するために、右辺を$\delta\boldsymbol{a}/\delta t$と表すことにする。これに対して、回転する座標系では、単位ベクトルが変化するため$d\boldsymbol{a}/dt\neq \delta\boldsymbol{a}/\delta t$である。すなわち、回転する座標系の座標を$\tilde{x},\ \tilde{y},\ \tilde{z}$、単位ベクトルを$\tilde{\boldsymbol{i}},\ \tilde{\boldsymbol{j}},\ \tilde{\boldsymbol{k}}$とすると

\frac{d\boldsymbol{a}}{dt}
=
\frac{\delta \boldsymbol{a}}{\delta t}
+
\tilde{x}\frac{d\tilde{\boldsymbol{i}}}{dt}
+
\tilde{y}\frac{d\tilde{\boldsymbol{j}}}{dt}
+
\tilde{z}\frac{d\tilde{\boldsymbol{k}}}{dt}

である。いま、座標系が原点まわりに角速度ベクトル$\boldsymbol{\omega}$で回転しているとすると、式(10')より、単位ベクトルを$\tilde{\boldsymbol{i}},\ \tilde{\boldsymbol{j}},\ \tilde{\boldsymbol{k}}$の時間微分は

\frac{d\tilde{\boldsymbol{i}}}{dt}
=
\boldsymbol{\omega} \times \tilde{\boldsymbol{i}},\

\frac{d\tilde{\boldsymbol{j}}}{dt}
=
\boldsymbol{\omega} \times \tilde{\boldsymbol{j}},\

\frac{d\tilde{\boldsymbol{k}}}{dt}
=
\boldsymbol{\omega} \times \tilde{\boldsymbol{k}},\

である。よってベクトル$\boldsymbol{a}$そのものの時間微分$d\boldsymbol{a}/dt$は

\frac{d\boldsymbol{a}}{dt}
=
\frac{\delta \boldsymbol{a}}{\delta t}
+
\boldsymbol{\omega} \times \boldsymbol{a}
\tag{15}

となる。この公式を使ってもう一度角運動量ベクトルと慣性モーメントテンソルを求めてみる。まず回転する座標系で質点の座標は変化しないので$\delta \boldsymbol{r}_i/\delta t = \boldsymbol{0}$である。よって質点の位置ベクトルの時間微分は

\frac{d\boldsymbol{r}_i}{dt}
=
\frac{\delta \boldsymbol{r}_i}{\delta t}
+
\boldsymbol{\omega} \times \boldsymbol{r}_i
=
\boldsymbol{\omega} \times \boldsymbol{r}_i

である。原点が移動しない限り、位置ベクトルの時間微分は座標の取り方によらないので、これは慣性座標系での結果と同じである。したがって、前の節と同じ式展開を繰り返すと、角運動量ベクトル$\boldsymbol{L}$と慣性モーメントテンソル$\tilde{I}$は

\boldsymbol{L}
=
\sum_{i} m_i \boldsymbol{r}_i \times (\boldsymbol{\omega} \times \boldsymbol{r}_i)
=
\tilde{I}
\boldsymbol{w}
\tag{16}

\tilde{I}
=
\begin{bmatrix}
\sum_i m_i (\tilde{y}_i^2 + \tilde{z}_i^2) & -\sum_i m_i \tilde{x}_i \tilde{y}_i & -\sum_i m_i \tilde{z}_i \tilde{x}_i\\
-\sum_i m_i \tilde{x}_i \tilde{y}_i & \sum_i m_i (\tilde{z}_i^2 + \tilde{x}_i^2) & -\sum_i m_i \tilde{y}_i \tilde{z}_i \\
-\sum_i m_i \tilde{z}_i \tilde{x}_i & -\sum_i m_i \tilde{y}_i \tilde{z}_i & \sum_i m_i (\tilde{x}_i^2 + \tilde{y}_i^2)
\end{bmatrix}
\tag{17}

となる。慣性モーメントテンソルの式の形は慣性座標系の場合と同じだが、使われている座標はすべて不変である。角運動量ベクトル$\boldsymbol{L}$を時間微分すると

\frac{d\boldsymbol{L}}{dt}
=
\frac{\delta\boldsymbol{L}}{\delta t}
+
\boldsymbol{\omega} \times \boldsymbol{L}
=
\tilde{I}
\frac{\delta\boldsymbol{\omega}}{\delta t}
+
\boldsymbol{\omega} \times \tilde{I}\boldsymbol{\omega}

である。ここで

\frac{d\boldsymbol{\omega}}{d t}
=
\frac{\delta\boldsymbol{\omega}}{\delta t}
+
\boldsymbol{\omega} \times \boldsymbol{\omega}
=
\frac{\delta\boldsymbol{\omega}}{\delta t}

より

\frac{d\boldsymbol{L}}{dt}
=
\tilde{I}
\frac{d\boldsymbol{\omega}}{d t}
+
\boldsymbol{\omega} \times \tilde{I}\boldsymbol{\omega}
\tag{18}

である。なお

\frac{d\boldsymbol{L}}{dt}
=
\frac{d}{dt}
(\tilde{I} \boldsymbol{\omega})
=
\tilde{I} \frac{d\boldsymbol{\omega}}{dt}

とするのは誤りである。なぜならベクトルに行列をかけるという操作は、ベクトルの成分の表し方、つまり座標の取り方に依存するからである。以上より、3次元の回転の運動方程式は

\tilde{I}
\frac{d\boldsymbol{\omega}}{d t}
+
\boldsymbol{\omega} \times \tilde{I}\boldsymbol{\omega}
=
\boldsymbol{N}
\tag{19}

となる。

慣性モーメントテンソルの対角化とオイラーの方程式

3次元の回転の運動方程式(19)には、慣性モーメントテンソル$\tilde{I}$が2回出てくる。もし$\tilde{I}$が対角行列であれば、つまり慣性乗積がすべてゼロであれば、計算がかなり簡単になるはずである。$\tilde{I}$は実対称行列であるので、正規直交行列$U$によって$U^T\tilde{I}U = \Lambda := \rm{diag}(\lambda_1,\lambda_2,\lambda_3 )$と対角化可能である。よって式(16)より$U^T\boldsymbol{L} = \Lambda U^T \boldsymbol{\omega}$である。そこで剛体とともに回転する新しい座標を

\tilde{\boldsymbol{i}}'
=
U_{11}\tilde{\boldsymbol{i}}
+
U_{21}\tilde{\boldsymbol{j}}
+
U_{31}\tilde{\boldsymbol{k}}

\tilde{\boldsymbol{j}}'
=
U_{12}\tilde{\boldsymbol{i}}
+
U_{22}\tilde{\boldsymbol{j}}
+
U_{32}\tilde{\boldsymbol{k}}

\tilde{\boldsymbol{k}}'
=
U_{13}\tilde{\boldsymbol{i}}
+
U_{23}\tilde{\boldsymbol{j}}
+
U_{33}\tilde{\boldsymbol{k}}

ととると、（成分表示された）角運動量ベクトルおよび角速度ベクトルはそれぞれ
$\boldsymbol{L}' = U^T \boldsymbol{L},\ \boldsymbol{\omega}' = U^T \boldsymbol{\omega}$と変換される。ただし$U$の各成分を$U_{ij}\ (1\leq i\leq 3,\ 1\leq j\leq 3)$とおいた。この変換により$\boldsymbol{L}' = \Lambda \boldsymbol{\omega}'$となり、慣性モーメントテンソルは対角行列となる。このような基底ベクトル$\tilde{\boldsymbol{i}}',\ \tilde{\boldsymbol{j}}',\ \tilde{\boldsymbol{k}}'$を剛体の慣性主軸とよぶ。

いま、最初から慣性主軸を基底ベクトルとして選んで回転運動の方程式を立てると

\begin{bmatrix}
\tilde{I}_x & 0 & 0 \\
0 & \tilde{I}_y & 0 \\
0 & 0 & \tilde{I}_z \\
\end{bmatrix}
\frac{d\boldsymbol{\omega}}{d t}
+
\boldsymbol{\omega} \times \begin{bmatrix}
\tilde{I}_x & 0 & 0 \\
0 & \tilde{I}_y & 0 \\
0 & 0 & \tilde{I}_z \\
\end{bmatrix}
\boldsymbol{\omega}
=
\boldsymbol{N}
\tag{20}

と書ける。$\tilde{I}_x,\ \tilde{I}_y,\ \tilde{I}_z$は慣性モーメントテンソル$\tilde{I}$の対角成分である。これは成分ごとに

\tilde{I}_x \frac{d\omega_x}{d t}
+
(\tilde{I}_z - \tilde{I}_y)\omega_y \omega_z
=
N_x

\tilde{I}_y \frac{d\omega_y}{d t}
+
(\tilde{I}_x - \tilde{I}_z)\omega_z \omega_x
=
N_y

\tilde{I}_z \frac{d\omega_z}{d t}
+
(\tilde{I}_y - \tilde{I}_x)\omega_x \omega_y
=
N_z

と簡潔に書き下すことができる。これらの式はオイラーの方程式とよばれる。

参考文献

• 吉田春夫. (1996). キーポイント力学（物理のキーポイント1 ）. 岩波書店.
• 物理のかぎプロジェクト. 慣性モーメント. (2010). 物理のかぎしっぽ. https://hooktail.sub.jp/mechanics/momentOfInertia/
• 物理のかぎプロジェクト. 速度座標系と慣性力. (2010). 物理のかぎしっぽ. http://hooktail.sub.jp/mechanics/acCoordinates/