曲面を滑らずに転がる剛体球運動の計算アルゴリズム【HTML, WebGL】

前回は斜面を滑らずに転がる計算アルゴリズムの導出を行ったね。並進運動と回転運動に関する運動方程式

$$
M\,\frac{d^2\boldsymbol{R}}{dt^2} =M\boldsymbol{g}+\boldsymbol{D}+\boldsymbol{f}
$$

$$
I\,\frac{d^2\boldsymbol{\theta}}{dt^2} =-a\hat{\boldsymbol{n}}\times\boldsymbol{f}
$$

は斜面に限らず任意の曲面上を滑る場合でも成り立つので、今回は半球内を滑らずに転がる剛体球の計算アルゴリズムの導出を行うよ。

前回と同様に図のような球面を滑る剛体球の運動の拘束条件を考えるよ。

並進運動の拘束条件

曲面と球面との接点を$${\boldsymbol{R}_0}$$と表した場合、剛体球の重心位置は

$$
\boldsymbol{R} =\boldsymbol{R}_0+a\,\hat{\boldsymbol{n}}
$$

を満たすね。この表式自体は斜面の場合と同じだけれども、斜面の場合には一定だった法線ベクトルが剛体球の位置に依存するね。位置ベクトルの原点を球面（半径$${A}$$）の中心をとした場合、法線ベクトルは

$$
\hat{\boldsymbol{n}} = -\frac{ \boldsymbol{R} }{ |\boldsymbol{R}| } = -\frac{ \boldsymbol{R} }{ A-a }
$$

となるね。法線ベクトルの内積はどんなときでも１なので、球面上を運動する剛体球の運動の拘束条件は

$$
\boldsymbol{R}\cdot \boldsymbol{R} = (A-a)^2
$$

となるね。さらにこの表式を時間で微分した

$$
\frac{d\boldsymbol{R}}{dt}\cdot \boldsymbol{R} = 0 ,     \frac{d^2\boldsymbol{R}}{dt^2}\cdot \boldsymbol{R} + \left|\frac{d\boldsymbol{R}}{dt}\right|^2= 0
$$

も拘束条件となるね。これらの条件を運動方程式に課すことで、垂直抗力$${\boldsymbol{S}}$$が決定できるよ。

回転運動の拘束条件

速度と角速度と法線の関係は斜面のときと一緒だね。

$$
\frac{d\boldsymbol{R}}{dt} = a\boldsymbol{\omega}\times \hat{\boldsymbol{n}} \ , \ \boldsymbol{\omega}=\frac{ d\boldsymbol{\theta}}{dt}
$$

斜面の場合と違うのは法線ベクトルが位置に依存するので、

$$
\frac{d\boldsymbol{R}}{dt} =-\frac{a}{A-a}\,\boldsymbol{\omega}\times \boldsymbol{R}
$$

となるね。これは回転運動に関する拘束条件の１つ目だね。さらに両辺を時間で微分した

$$
\frac{d^2\boldsymbol{R}}{dt^2} = a\,\frac{d^2\boldsymbol{\theta}}{dt^2}\times \hat{\boldsymbol{n}}
$$

も拘束条件となるね。これらの拘束条件を用いることで摩擦力$${\boldsymbol{f}}$$を決定できるよ。

垂直抗力の表式

斜面の場合、垂直抗力$${\boldsymbol{S}}$$は位置によらず一意に決定できたのだけれども、曲面の場合は見かけの力である遠心力が働くので、速度にも依存することになるね。並進運動の運動方程式の両辺に\hat{\boldsymbol{R}}との内積をとった

$$
M\,\frac{d^2\boldsymbol{R}}{dt^2}\cdot \boldsymbol{R} =M\boldsymbol{g}\cdot \boldsymbol{R}+\boldsymbol{S} \cdot \boldsymbol{R}
$$

に先の関係式を代入して整理すると、垂直抗力は

$$
\boldsymbol{S} = \left[\frac{M}{A-a}\left|\frac{d\boldsymbol{R}}{dt}\right|^2-M\boldsymbol{g}\cdot \hat{\boldsymbol{n}}\right] \hat{\boldsymbol{n}}
$$

が得られるね。速度の２乗に比例した第１項が遠心力に対抗する垂直抗力、第２項は重力に対する垂直抗力を表すね。この導出には$${\boldsymbol{S} = S \hat{\boldsymbol{n}}}$$であることを考慮しているよ。

摩擦力の表式

摩擦力$${\boldsymbol{f}}$$を計算するために、まず回転運動に関する運動方程式の両辺に右から$${\boldsymbol{R}}$$との外積をとってベクトル三重積の変換を行うよ。

$$
I \,\frac{d^2\boldsymbol{\theta}}{dt^2} \times \boldsymbol{R} = a(A-a) \boldsymbol{f}
$$

この式に回転に関する拘束条件を課すと、

$$
\frac{d^2\boldsymbol{R}}{dt^2} =- \frac{a^2}{I} \,\boldsymbol{f} -\frac{a}{A-a}\, \boldsymbol{\omega} \times \frac{d\boldsymbol{R}}{dt}
$$

が得られるので、これと先に導出した$${\boldsymbol{S}}$$を並進運動の運動方程式に代入するすることで$${\boldsymbol{f}}$$に関する条件式が得られるよ。するとなんと斜面の際と全く同じ表式

$$
\boldsymbol{f}=\frac{M}{1+\frac{Ma^2}{I}} \left[(\boldsymbol{g}\cdot \hat{\boldsymbol{n}})\, \hat{\boldsymbol{n}}-\boldsymbol{g}\right]
$$

が得られるね。導出過程で$${|\boldsymbol{V}|^2 = a^2 |\boldsymbol{\omega}|^2}$$の関係を考慮しているよ。それにしても不思議なのは垂直抗力に速度依存性があるのに、摩擦力が速度に依存しないのはちょっと不思議だね。

曲面を転がる剛体球の計算アルゴリズム

以上の結果を踏まえて、曲面を滑らずに転がる剛体球の計算アルゴリズムは次のとおりだよ。　

$$
\frac{d^2\boldsymbol{R}}{dt^2} = \frac{1}{1+\frac{I}{Ma^2}} \left[\boldsymbol{g} -(\boldsymbol{g}\cdot \hat{\boldsymbol{n}})\, \hat{\boldsymbol{n}}\right] + \frac{1}{A-a} \,|\boldsymbol{V}|^2 \hat{\boldsymbol{n}}
$$

$$
\frac{d^2\boldsymbol{\theta}}{dt^2} = -\frac{1}{a}\,\frac{d^2\boldsymbol{R}}{dt^2}\times \hat{\boldsymbol{n}}
$$

並進運動の第１項目は斜面のときと同じで、第２項目が曲面によって生じる向心力に相当するね。回転運動は斜面と全く一緒となるね。剛体球の慣性モーメント

$$
I = \frac{2}{5}\, Ma^2
$$

を考慮すると、並進運動の第１項目は質量にも半径にも依存しないけれども、第２項目は半径に依存するね。つまり、斜面のときとは異なり、運動は半径に依存することがわかるね。

シミュレーション結果

次のアニメーションは斜面の傾きの向きからずらした方向へ初速度（角速度）を与えたときの計算結果だよ。表示している矢印は角速度ベクトル$${\boldsymbol{\omega}}$$ を表しているよ。また、左上の数値は計算誤差の大きさを把握するために計算した力学的エネルギーの変化だよ。意図通りの運動結果が得られていそうだね。

プログラムソース（HTML）

上記シミュレーションの計算と動画作成のプログラムソースを販売するよ。もし良かったら試してみてください。