なぜ回っているコマは倒れないの？という疑問は多くの方が持っていると思います。
私が「コマが倒れない理由」でネット検索しても、「ジャイロ効果があるから」「歳差運動するから」といった解説しかされないことが多く、「そのジャイロ効果ってそもそもなんなんだ」という疑問に答えてくれませんでした。
そのため、「ジャイロ効果」を詳しく調べる準備として、コマの動きをシミュレーションを試みました。コマの運きを記述する運動方程式をたてて、pythonに解かせることで、確かに倒れないコマのシミュレーションが意外にも簡単にできたので紹介します。
オイラー方程式とオイラー角を使う必要がある
コマの動きを知るには剛体の「オイラー方程式」が必要
コマのような大きさのある物体の運動を知るためには、いわゆる「剛体の力学」を使う必要があります。
高校までの力学は「質点の力学」だよね
このページでは導出を割愛しますが、形のある物体の動きをは、$\frac{\text{d}\boldsymbol{L}}{\text{d}t}=\boldsymbol{N}$で表される「オイラー方程式」を解けば知ることができます。
このオイラー方程式は以下のようになります。
$\begin{eqnarray}
\left\{
\begin{array}{l}
A\frac{\text{d}\omega_1}{\text{d}t}-(B-C)\omega_2\omega_3=N_1\\
B\frac{\text{d}\omega_2}{\text{d}t}-(C-A)\omega_3\omega_1=N_2\\
C\frac{\text{d}\omega_3}{\text{d}t}-(A-B)\omega_1\omega_2=N_3
\end{array}
\right.
\end{eqnarray}$
ここで、$\omega_1,\omega_2,\omega_3$は角速度$\boldsymbol{\omega}$のそれぞれ慣性主軸$\xi,\eta,\zeta$成分のことです。
オイラー角を使う
「オイラー方程式」という「ひな形」に対して、$A$や$\omega$、$N$などの値を具体的に入れなければ、コマの動きはシミュレーションできません。
コマは3次元で複雑に動くのに、角運動量$\omega$はどうやって記述したらいいんだ？
３次元での物体の姿勢を表現するのによく用いられるのが「オイラー角」です。
「オイラー角」を使うと、「コマが垂直からどれだけ倒れているか」「コマがこれだけ自転、公転しているか」を定量的に書き下すことができます。
オイラー角はよく使われるので、角速度などの表記は文献が多いです
文献が多いと、細かい部分は覚えていなくてもすぐ調べられるから便利だね
オイラー角には回転の順番によって12個もの種類があります。文献値を使う場合は、「どの順番で回転させて得られたオイラー角か」に十分注意しましょう。
微分方程式を解けばコマの動きが分かる
「剛体のオイラー方程式」に「オイラー角」で表した$\omega_1,\omega_2,\omega_3$を代入しましょう。あわせて、重力によるモーメント（要するにコマを床に倒そうとする力）も加味すると、次のようになります。
$\begin{eqnarray}
\left\{
\begin{array}{l}
A\ddot{\theta}-A\dot{\phi}^2\sin\theta\cos\theta+Cn\dot{\phi}\sin\theta=Mgh\sin\theta\\
A\ddot{\phi}\sin\theta+2A\dot{\phi}\dot{\theta}\cos\theta-Cn\dot{\theta}=0\\
\dot{\phi}\cos\theta+\dot{\psi}=n
\end{array}
\right.
\end{eqnarray}$
この式では、コマの対称性（$A=B$）を利用しており、$h$は床からコマの重心までの長さ、$n$は（要するに）コマの回転数($[\text{rad}/\text{s}]$)です。
回るコマのシミュレーションに必要な運動方程式を導出から詳しく説明してくれているのが「原田 鮮著 力学Ⅰ 質点・剛体の力学」です。
本書では剛体の運動の他、コリオリ力や強制振動など、大学の力学で取り扱うテーマを重点的に取り扱っています。
シミュレーションしてみた
さて、実際に上の微分方程式をpythonでプログラムしてシミュレーションしたものが↓になります。
- コマの大きさは$r_1=5\text{mm}$、$r_2=50\text{cm}$、$r_3=70\text{cm}$、$r_4=10\text{cm}$、$h_1=5\text{cm}$、$h_2=50\text{cm}$、$h_3=30\text{cm}$としている。
- 初期のコマの傾き$\theta$は$60$度とする。
- 初期のコマの回転数は$2.5$回転/秒とする。
プログラミングコードはこちら
import math
import numpy as np
import numpy.linalg as LA
np.set_printoptions(precision=6,floatmode='fixed')
PI = 3.1415926
theta = 1/6*PI
phi = 0
psi = 0
eta = 0
gzai = 0
d_psi = 0
rho = 1
r1 = 0.005
r2 = 0.5
r3 = 0.7
r4 = 0.1
h1 = 0.05
h2 = 0.5
h3 = 0.3
Dt = 0.005
g = 9.80665
n = 2.5*2*PI
alpha = (r3-r2)/h2
lz1 = PI*rho/2*(r1**4)*h1
lz2 = PI*rho/(10*alpha)*((r2+alpha*h2)**5-r2**5)
lz3 = PI*rho/2*(r4**4)*h3
lx1 = PI*rho*(r1**2)*h1*((h1**2)/3+(r1**2)/4)
lx2 = PI*rho/(20*alpha)*((r2+alpha*h2)**5-r2**5)+PI*rho*((alpha**2)/5*(h2**5)+(r2*alpha+h1*(alpha**2))/2*(h2**4)+((r2**2)+4*h1*r2*alpha+(h1**2)*(alpha**2))/3*(h2**3)+(h1*(r2**2)+(h1**2)*r2*alpha)*(h2**2)+(h1**2)*(r2**2)*h2)
lx3 = PI*rho*(r4**2)*(((h1+h2+h3)**3)/3+(r4**2)*(h1+h2+h3)/4-((h1+h2)**3)/3-(r4**2)*(h1+h2)/4)
A = lx1 + lx2 + lx3
C = lz1 + lz2 + lz3
alpha = (r3-r2)/h2
total_V = PI*((r1**2)*h1+1/3*(r2**2+r2*r3+r3**2)*h2+(r4**2)*h3)
for ite3 in range(1000):
x = ite3*h2/1000
half_V = PI*((r1**2)*h1+(r2**2)*x+r2*alpha*(x**2)+1/3*(alpha**2)*(x**3))
if(half_V>0.5*total_V):
hg = h1 + x
break
Mgh = rho * total_V * g * hg
e = open('spin.txt','a')
f = open('spin_koma.txt','a')
g = open('spin_theta.txt','a')
ite = 0
for ite in range(20000):
k_theta1 = eta * Dt
k_phi1 = gzai * Dt
k_eta1 = ((A*(gzai**2)*math.cos(theta)-C*n*gzai+Mgh)*math.sin(theta))/A*Dt
k_gzai1 = -(2*A*eta*gzai*math.cos(theta)-C*n*eta)/(A*math.sin(theta))*Dt
k_psi1 = (n-gzai*math.cos(theta))*Dt
k_theta2 = (eta+k_eta1/2) * Dt
k_phi2 = (gzai+k_gzai1/2) * Dt
k_eta2 = ((A*((gzai+k_gzai1/2)**2)*math.cos(theta+k_theta1/2)-C*n*(gzai+k_gzai1/2)+Mgh)*math.sin(theta+k_theta1/2))/A*Dt
k_gzai2 = -(2*A*(eta+k_eta1/2)*(gzai+k_gzai1/2)*math.cos(theta+k_theta1/2)-C*n*(eta+k_eta1/2))/(A*(math.sin(theta+k_theta1/2)))*Dt
k_psi2 = (n-(gzai+k_gzai1/2)*math.cos(theta+k_theta1/2))*Dt
k_theta3 = (eta+k_eta2/2) * Dt
k_phi3 = (gzai+k_gzai2/2) * Dt
k_eta3 = ((A*((gzai+k_gzai2/2)**2)*math.cos(theta+k_theta2/2)-C*n*(gzai+k_gzai2/2)+Mgh)*math.sin(theta+k_theta2/2))/A*Dt
k_gzai3 = -(2*A*(eta+k_eta2/2)*(gzai+k_gzai2/2)*math.cos(theta+k_theta2/2)-C*n*(eta+k_eta2/2))/(A*(math.sin(theta+k_theta2/2)))*Dt
k_psi3 = (n-(gzai+k_gzai2/2)*math.cos(theta+k_theta2/2))*Dt
k_theta4 = (eta+k_eta3) * Dt
k_phi4 = (gzai+k_gzai3) * Dt
k_eta4 = ((A*((gzai+k_gzai3)**2)*math.cos(theta+k_theta3)-C*n*(gzai+k_gzai3)+Mgh)*math.sin(theta+k_theta3))/A*Dt
k_gzai4 = -(2*A*(eta+k_eta3)*(gzai+k_gzai3)*math.cos(theta+k_theta3)-C*n*(eta+k_eta3))/(A*(math.sin(theta+k_theta3)))*Dt
k_psi4 = (n-(gzai+k_gzai3)*math.cos(theta+k_theta3))*Dt
theta = theta + (k_theta1+2*k_theta2+2*k_theta3+k_theta4)/6
phi = phi + (k_phi1 +2*k_phi2 +2*k_phi3 +k_phi4 )/6
eta = eta + (k_eta1 +2*k_eta2 +2*k_eta3 +k_eta4 )/6
gzai = gzai + (k_gzai1 +2*k_gzai2 +2*k_gzai3 +k_gzai4 )/6
psi = psi + (k_psi1 +2*k_psi2 +2*k_psi3 +k_psi4 )/6
d_psi = (k_psi1 +2*k_psi2 +2*k_psi3 +k_psi4 )/6
u0x = math.sin(theta) * math.cos(phi)
u0y = math.sin(theta) * math.sin(phi)
u0z = math.cos(theta)
u1x = math.sin(theta) * math.cos(phi)
u1y = math.sin(theta) * math.sin(phi)
u1z = math.cos(theta)
u2x = math.cos(theta) * math.cos(phi)
u2y = math.cos(theta) * math.sin(phi)
u2z = -math.sin(theta)
u3x = -math.sin(phi)
u3y = math.cos(phi)
u3z = 0
ab = 0.6
cd = 0.3
ef = 0.2
gh = 0.2
if (ite % 10) == 0:
co2 = math.cos(psi)
si2 = math.sin(psi)
print(ite,ite*Dt,0,0,0,1.5*u1x,1.5*u1y,1.5*u1z,gh*u1x,gh*u1y,gh*u1z,gh*u1x+cd*(co2*u2x+si2*u3x), gh*u1y+cd*(co2*u2y+si2*u3y),gh*u1z+cd*(co2*u2z+si2*u3z), file=e)
for ite2 in range(8):
co = math.cos(2*PI*ite2/8)
si = math.sin(2*PI*ite2/8)
print(ab*u1x+ef*(co*u2x+si*u3x),ab*u1y+ef*(co*u2y+si*u3y),ab*u1z+ef*(co*u2z+si*u3z), file=f)
print(ite,ite*Dt,theta*180/PI,phi*180/PI,psi*180/PI,d_psi/(2*PI)/Dt,file=g)
e.close()
f.close()
g.close()
プログラムの正しさを確認
私が作成したプログラムの正しさを確認するため、「コマ シミュレーション」でインターネット検索し、コマの歳差運動の数値計算結果が掲載されているページを探しました。
すると、「高校物理で探検するコマのワンダーランド」というホームページ（URL：https://yamauo1945.sakura.ne.jp/komamove.html）が掲載の前提条件と計算結果をわかりやすく掲載していたので、計算条件を同じにして同じ結果が得られるか確認しました。
- コマの慣性モーメントは$A=B=0.135\text{kgm}^2/\text{s}$、$C=0.09\text{kgm}^2/\text{s}$とする。
- 重力加速度は$9.80665\text{m}/\text{s}^2$とする。
- 床から重心までの距離は$0.3\text{m}$とする。
- 初期のコマの傾き$\theta$は$60$度とする。
- 初期のコマの回転数は$2.73$回転/秒とする。
計算の結果、本ホームページのプログラムによる結果と「高校物理で探検するコマのワンダーランド」による結果はぴたりと一致することが確認できました。
まとめ
このページでは、コマの運きを記述する運動方程式をたてて、pythonに解かせることで、コマが倒れないシミュレーション結果を紹介しました。
つまり「剛体の運動方程式を解けば回るコマが倒れないシミュレーションはできる」というところまで明らかにすることができました。しかし、残念ながら、コマが倒れない「定性的な理解」までは落とし込めていないように感じます。
さらに踏み込んで「回るコマが倒れない理由を定性的に理解する」ことについては別のページで紹介したいと思います。。
