卵を回すとなぜ立つか—– 剛体の力学 —–福井大学 工学部 物理工学科
2004 年 2 月

村井興平

目次
第 1 章 剛体の運動方程式
3
1.1 剛体とそのつりあい . . . . . . . . . . . . . . . . . . . . . 3
1.1.1 剛体とその自由度 . . . . . . . . . . . . . . . . . . . 3
1.1.2 剛体のつりあい . . . . . . . . . . . . . . . . . . . . 4
1.2 固定軸のある剛体の運動 . . . . . . . . . . . . . . . . . . . 6
1.3 慣性モーメント . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 固定点のある剛体の運動 . . . . . . . . . . . . . . . . . . . 9
1.4.1 慣性テンソル . . . . . . . . . . . . . . . . . . . . . 9
1.4.2 オイラーの方程式 . . . . . . . . . . . . . . . . . . . 11
第 2 章 先行する研究の概略

13

第 3 章 偏心球の回転直立の理論

14

第 4 章 軸対称物体の回転直立の理論

16

第 5 章 運動方程式の数値積分による回転直立の再現

22

第 6 章 数値計算の結果
25
6.1 軸対称物体のオイラー法による解 . . . . . . . . . . . . . . 25
6.2 ルンゲクッタ法により求めた回転楕円体の直立運動 . . . . 29
第 7 章 結論

35

1

はじめに

最近 H.K.Moﬀatt と下村裕による速報論文 “Spinning eggs a paradox
resolved“が Nature 誌に掲載された [1]。これはゆで卵のように軸対称の物体を机上で回転させると立ち上がる現象を統一的に扱った論文である。
これに先立って、重心のずれた球体の回転については長い研究の歴史がある。具体的には、重心の低い姿勢で回転させると自然と傾きはじめ、やがて重心の最も高い姿勢で回転するという不思議な運動をする「 逆さゴマ」を力学的に説明しようとする一連の試みがあり、1950 年代に最終的解決をみたのであった。Moﬀatt 、下村は逆さゴ マを軸対称ゴ マに拡張することで研究に新局面をひらいたわけである。
彼らは近似的にこの理論を説明したわけであるが、私は近似のない数値
計算によって剛体の直立運動を再現することを本研究の目的とする。

2

第1章

剛体の運動方程式

この章の執筆にあたっては文献 [6] を参考にした。

1.1

剛体とそのつりあい

1.1.1

剛体とその自由度

我々の身辺には、形がほぼ一定の物体が多い。これを理想化し 、形が全く変わらない物体を考えて、これを剛体と呼ぶ。もっと正確には、質点間の距離 |r i − r j | がすべていつも一定に保たれている質点系である。
n 個の質点からできている系の運動を指定するには、x1 , y1 , z1 , x2 , y2 , z2
, .., xn , yn , zn という 3n 個の変数の時間変化を知る必要がある。このとき
系は 3n 個の自由度を持つと言われる。物体を剛体と見なすと、6 個の変数で剛体の運動は完全に記述することができる。
剛体がその形を保持しているのは、各部分の間に働く内力によるのであるが 、内力は全く表面にでてこないで話がすむ。というのは、自由度が 6の剛体の運動は、重心の運動（ 自由度 3 ）と、重心の回りの相対運動として許される回転（ 自由度 3 ）で記述されるが、質点系の一般論で調べたよ
うに、これらは外力だけで決まるからである。一様な重力は、重心の回りのモーメントの和が 0 で 、剛体を重心の回りで回転させようとする働きを持たないから、重心運動だけをきめるものとして、重心に働く大きさがM G の一つの力として扱ってよい。重心が剛体に固定した点であることは、言うまでもない。

3

1.1.2

剛体のつりあい

剛体がつりあいを保つためには、重心が静止を続けなければならない。
そこで質量 m1 , m2 , m3 , ... をもった多数の質点からできている質点系を考える。それら質点の位置を r 1 , r 2 , r 3 ,... とし 、i 番目の質点が、j 番目におよぼす力を F ij , これら内力のほかに質点系の外部から働く外力もあってよいから、i 番目の質点が受けている外力を F i としめすことにする。そうすると、各質点の運動方程式は

m1 r¨1 = F 1 + F 21 + F 31 + ....
m2 r¨2 = F 2 + F 22 + F 32 + ....
..................

(1.1)となる。これらを全部加え合わせると、内力の部分は作用反作用の法則
F ij + F ji = 0 によってことごとく打ち消し合ってしまうから、右辺の和には外力だけが残り

m1 r¨1 + m2 r¨2 + .... = F 1 + F 2 + ...

(1.2)

が得られる。質点系の全質量をM = m1 + m2 + ...

(1.3)とし 、重心または質量中心を

M R = m1 r 1 + m2 r 2 + ...

(1.4)

によって定義すると、上の式は

M

∑
d2 R
=
F
+
F
+
...
=
Fi
1
2
dt2
i

(1.5)
∑

となる。これは、質量が M の 1 個の質点（その位置を R とする）が i F iという力を受けて運動する場合の運動方程式の形をしている。つまり、質点系の重心はそこに系の全部の質量と外力のすべてが集中した 1 個の質点のように運動する。
（ 1.1 ）式を用いると∑dP=
Fi = 0dtiが得られる。

4(1.6)

次に、たとえ重心が静止していても、そのまわりで回転が起っては静止とは言えないから 、重心の回りの角運動量 L0 が 0 のままで変化しない、ということを保証する条件として、
∑
d ∑ 0
d ∑ 0
ri =
r i × mi r˙i 0 =r 0i × F i = 0
dt i
dt i
i

(1.7)

が要求される。しかし 、重心が動かないのであるから LG も 0 のままであり、LG + L0 = 0 である。したがって

dL ∑=ri × F i = 0
dt
i

(1.8)

でも同じことになる。結局 (1.2) と (1.3) の 2 式を成分に分けて書いた合計 6 個の式が 、つりあいの条件ということになる。

5

1.2

固定軸のある剛体の運動

剛体が固定した軸のまわりで回転する場合には、自由度は 1 になる。今
固定軸を z 軸にとり、剛体の適当なところにつけた目印から z 軸に下した垂線が zx 面とつくる角 θ が t とともにど う変化するかを知ればよいからである。
剛体を細分して、その各小片を質点とみなすことにし 、i 番目の質点の質量を mi , 座標を xi , yi , zi とすると、√

ri =x2i + yi2 = 一定, zi = 一定

(1.9)

であり、

xi = ri cos θi ,yi = ri sin θi(1.10)

であるが 、θi − θ は t によらないので

θ˙i = θ̇,θ¨i = θ̈,

(1.11)となる。したがってẋi = −ri sin θi × θ˙i = −yi θ̇,

同様に y˙i = xi θ̇

(1.12)がわかる。以上を使って θ(t) を求めるには 、質点系の角運動量に関する式（ 1.3 ）式を利用すればよい。それは∑d ∑mi (r i × r˙i ) =rj × F j
dt ij(1.13)

であるが 、左辺のベクトル積の z 成分を考えると(r i × r˙i )z = xi y˙i − yi ẋi = x2i θ̇ + yi2 θ̇ = r 2i θ̇

(1.14)であって ri2 は t によらないから∑d ∑mi (r i × r˙i )z = ( mi r 2i )θ̈
dt i
i(1.15)

が得られる。

Iz =∑mi (x2i + yi2 )iを、この剛体の z 軸のまわりの慣性モーメントという。

6(1.16)

これから、今の場合には、剛体の全角運動量の z 成分の時間変化が Iz θ̈と表されることがわかる。θ̈ は角加速度である。以上によって（ 1.3 ）式の
z 成分は

Iz

d2 θ ∑=(xj Fiy − yi Fjx )dt2
j(1.17)となる。右辺は、剛体に働く力のモーメントのうちで実際に回転を左右す
る成分であるが 、これを特にトルクと呼ぶことがある。

7

1.3

慣性モーメント

密度が一様な球の慣性モーメントは、質量 を M とし 、半径を R とす
ると2
I = M R2
5

(1.18)

で与えられる。
密度が一様な楕円体の慣性モーメントは質量を M とし 、
x 軸方向の半径 Rx 、y 軸方向の半径 Ry 、z 軸方向の半径 Rz とすると

1
Ix = M (Ry2 + Rz2 )5
1
Iy = M (Rz2 + Rx2 )5
1
Iz = M (Rx2 + Ry2 )5で与えられる。

8(1.19)(1.20)(1.21)

1.4

固定点のある剛体の運動

1.4.1

慣性テンソル

固定点を O とし 、これを原点として空間に固定した座標系 O − xyz, 剛
体に固定した座標系 O − ξηδ をとる。
各瞬間における剛体の運動は、O を通る回転軸とその周りの角速度 ω で指定される。これを、その回転軸と方向が一致し 、大きさが ω で , そのときの剛体のように回転させた右ねじの進む向きを持つベクトル ω であらわすことにする。ω は、空間に固定されていないだけでなく、一般には剛体に対しても動くベクトルである。
角速度が ω のとき剛体は時間 dt のあいだに ω のまわりに角 ωdt だけ回転するから、剛体上の点 r i は、ω と r i の両方に垂直な方向に (|r i | sin γi )ωdt
だけ動く。γi は r i と ω の間の角である。これから向きまで含めて、この瞬間における速度がṙ i = ω × r i(1.22)と表されることがわかる。
この式に mi を掛けて i についての和をとれば 、剛体の運動量はP =∑mi r˙i = ω ×∑

mi r i = ω × M R(1.23)ii∑と表される。M = mi , R は重心の位置ベクトルである。つぎに、O 点のまわりの角運動量を考えるとL=∑

mi (r i × r˙i ) =∑imi r i × (ω × r i )(1.24)iであるが(r × (ω × r))x = y(ωx y − ωy x) − z(ωz x − ωx z)= (y 2 + z 2 )ωx − xyωy − zxωz

(1.25)

を用いれば

Lx =

∑

mi (yi2 + x2i )ωx −i∑

m i xi y i ω y −i∑

mi zi xi ωz(1.26)iと書けることがわかる。剛体に固定した座標軸方向の成分で記せば

Lξ =

∑
i

mi (ηi2 + ζi2 )ωξ −∑i

mi ξi ηi ωη −

∑
i

となる。

9

mi ζi ξi ωζ

(1.27)

慣性モーメントと慣性乗積の定義を用いると、Lξ = Iξ ωξ − Iξη ωη − Iζξ ωη

(1.28)となる。ほかの成分も同様にして

Lη = −Iξη ωη + Iζ ωζ − Iζξ ωζ

(1.29)

Lζ = −Iζξ ωξ − Iηζ ωη + Iζ ωζ(1.30)のように書けるから、まとめて行列で表せば








Lξ
Iξ
−Iξη −Iζξ

 

Iη
−Iηζ 
 Lη  =  −IξηLζ−Iζξ −Iηζ
Iζ





ωξ


 ωη 
ωζ

(1.31)となる。この右辺にある 3 行 3 列の行列で表されている量のことを慣性テンソルという。すでに述べたように座標軸の方向を適当に選べば 、慣性乗積がすべて 0 になるようにすることができるので、そのようにすれば








L0
I0 0 0 ξ   ξ
 L0  =  0 I 0 0 η η  
L0ζ
0 0 Iζ0





ω0
 ξ 
 ω0  η ωζ0(1.32)

このときの Iξ0 , Iη0 , Iζ0 を主慣性モーメント、座標軸の方向を慣性主軸という。

10

1.4.2オイラーの方程式

いま、ベクトル B(t) の微少時間 dt における変位 dB を ,(1) 動いているO − ξηζ 系でみた B の変化と、(2)B をこの間だけ O − ξηζ 系に固定させてそれといっしょに回転させた場合の B の変化とに分けて考えてみる。
r i の場合には前者 (1) はなくて、時間変化は後者 (2) のみであり、それが(ω × r i ) であった。一般のベクトル B でも同様に (2) が (ω × B)dt で表されることはすぐわかるであろう。(1) の ξ, η, ζ 成分は (dB ξ , dB η , dB ζ )である。したがって、これらのベクトル和である dB を dt で割ったもの
として

dBξ
dB
)ξ =
+ (ω × B)ξ
dt
dt
dB
dBη
(
)η =
+ (ω × B)η
dt
dt
dBζ
dB
)ζ =
+ (ω × B)ζ
(
dt
dt
(

(1.33)(1.34)(1.35)

が得られる。B = r i のときには右辺の第一項が消えるのである。
質点系一般の場合に成り立つ全角運動量の時間変化は外力のモーメン
トの総和に等しいという関係は剛体に対しても成り立つ。固定点 O に関する剛体の全角運動量 L は (1.31) あるいは (1.32) で与えられるが、ξ, η, ζ
軸を慣性主軸に一致させてあるとすれば (1.32) が使える。L と ω の肩符
0 を省略して書けばLξ = Iξ0 ωξ ,L
Ã η = Iη0 ωη ,Lξ = Iζ0 ωζ

であって、主慣性モーメント Iξ0 , Iη0 , Iζ0 は定数であるから

dLξ
dωξ
dηξ
dLζ
dωζ
dLη
= Iξ0
,
= Iη0
,
= Iζ0
dt
dt
dt
dt
dt
dt
と書かれる。そこで (1.4) (1.6) を L に対して適用してやるとdωξ
dL)ξ = Iξ0+ (ωη Lζ − ωζ Lη )
dt
dt
dωξ
= Iξ0
− (Iη0 − Iζ0 )ωη ωζ
dtなどが得られる。したがって
(

Iξ0

∑
dωξ
− (Iη0 − Iζ0 )ωη ωζ = ( r i × F i )ξ
dt
i

11

Iη0

∑
dωη
− (Iζ0 − Iξ0 )ωζ ωξ = ( r i × F i )η
dt
i

Iζ0

∑
dωζ
− (Iξ0 − Iη0 )ωζ ωη = ( r i × F i )ζdtiとなる。この式を、剛体の回転に対するオイラーの方程式と呼ぶ。

12

第2章

先行する研究の概略

1 、J.H Jellet(1872) は、重心のずれた球体の運動では、H を重心の高さのベクトル、X を重心から見た床との接点の位置ベクトルとすると
J = −H · X p(2.1)

は運動定数であることを最初に近似的に発見したので、J はジェレット定数と呼ばれている。

2 、E.J.Routh(1905)
すべる場合と転がる場合の両方について、ジェレット定数が運動定数であることを証明した。
3 、C.M.Brams, N.M.Hugonholts(1952)
逆立ゴ マの問題を最終的に解決した。
4 、M.K. Moﬀatt and Y.Shimomura(2002)
回転楕円体の場合を下村が１年かけて解いたが、その後任意の軸対称物体に拡張できることがわかり、卵の回転直立をはじめて説明する理論となっ
た。

13

第3章

偏心球の回転直立の理論

この章の執筆にあたっては文献 [7] を参考にした。
偏心球の例として逆立ゴマを例にあげる。滑らかな机の上で逆立ゴマをまわすとど うなるだろうか？
逆立ゴマの重心は曲率中心の近くにある。このことから、重力のモーメ
ントは無視してもよい。そのうえ、摩擦力のモーメントは水平方向である。この水平方向のモーメントも同様に角周波数 ω で回転を行っている
L ≈ 0 となり、から、時間平均をとれば 0 になる。この結果、平均して ddt
コマの角運動量はほとんど保存されている。逆立ゴマをはじめ上向きの回
転速度で回せば 、L の方向は近似的に机に対して常に鉛直な方向を向くことになる。この運動をコマとともに回転する座標系からみると、モーメン
ト N の時間変化が 0 である運動方程式は

N=

dL
+ω×L
dt

(3.1)

となる。いま重心の回りの３つの主慣性モーメントが互いに等しいものとすると

L ≈ Iω

(3.2)となり、(3.1) の中の ω × L の項は 0 となる。運動の間 L は終始鉛直方向
を、N は水平方向を向いている。L と N とは互いに垂直であるから、物体に固定した座標系からみるとき、(3.1) より、モーメントのために角運動量が一様な歳差運動を始めることになる。いま対称軸 k の方向の (3.1)の成分をとれば

k·N =k·dLd(k · L)
=
dt
dt

(3.3)

となる。N を近似的に水平方向とした図を用いれば 、(3.3) の２つのスカ
ラー積はそれぞれ

k · N = −N sin θk · L = L cos θ

14

(3.4)

となる。ここで θ は L と k の間の角である。(3.3) 、(3.4) から、R をコマの半径とすると

θ̇ =

N
µM gR
≈
L
Iω

(3.5)

が得られる。したがって、θ はたしかに時間が経つにつれて増加し 、逆立
ゴ マは引っくり返るわけである。

15

第4章

軸対称物体の回転直立の理論

この章の執筆に際しては文献 [3,8] を参考にした。
軸対称の物体を回転させると立ち上がる。H.K.Moﬀatt と下村は共同研究の結果この現象を表現する数学解を得た。ここではその理論の紹介をする。ここで用いる座標系を図に示した。

図：水平な机の上で重心 G のまわりに回転する物体
参考文献 [3] の図 3 を転載した
回転する物体の対称軸と重心 G 、そして机との接点 P を考え、その中で鉛直座標軸 GZ と水平座標軸 GX をとり、これらと直交し紙面を向こうへ抜ける座標軸を GY とする。机は十分滑らかであるとすると、物体が重心 G を通る Z 軸のまわりに回転するとき、接点 P は机上をすべって移動する。そのため座標系 GXY Z は運動座標系である。
対称軸を GZ 0 ,Z 軸とのなす角を θ とし 、Gz 0 に垂直な軸を Gx0 ,Gy 0 とする。また重心 G の高さを h(θ) とし 、接点 P の座標を GXYZ 系で XP =16

(XP , 0, ZP ) とする。
XP =

dh
,
dθ

ZP = −h(θ)(4.1)(4.1) で第２式は明らかであるが 、第 1 式は下村が研究の最終段階で発見した。この理論は任意の軸対象物体に適用できるという優れたものである。ここではそれを紹介する。

図は参考文献 [8] の < 図 2 > を転載した

17

図は角度 θ が微小角 ∆θ だけ変化した状況を、物体を固定して描いた図である。X ∗ , Y ∗ は変化後の座標軸で、h∗ は変化後の重心の高さである。
下の図は上の図の点線で囲んだ楕円部分を拡大したもので 、ここで線分
LN に注目する。
LN ≈ −XP ∆θ
(4.2)
また

LN

= OL − ONh∗= h−cos(∆θ)≈ h − h∗ = −∆h

ゆえに

XP =

dh
dh
XP =
dθ
dθ

(4.3)

(4.4)

となる。
また物体の全質量を M , 重力加速度を g とすると接点 P で上向きに働く力は R = M g であり、ここで摩擦力を F とすると、すべりの摩擦係数 µは十分小さく、F は Y 方向を向き F = (0, F, 0)

F = µM g(4.5)としてよい。この式は物体の上向の加速度を 0 として、摩擦によって回転が止まらないということを言いたいのである。
物体の角速度 ω は、鉛直軸 GZ のまわりの角速度 Ω と物体の軸 Gz 0 の
まわりの角速度 ψ̇ を合成したものと見ることができる。
これを物体の主軸方向に射影し 、これらの主軸方向の角速度を求めると

ωx0 = −Ω sin θ(4.6)

ωy0 = θ̇ωz0 = Ω cos θ + ψ̇ = n

(4.7)
(4.8)

となる。主軸方向の慣性モーメントを

Ix0 = Iy0 = A,

Iz0 = C

(4.9)

とすると角運動量は

L0x = Aωx0 = −AΩ sin θ(4.10)

L0y = −AΩ0y = Aθ̇
L0z = Cn

(4.11)

18(4.12)

であたえられる。また GXY Z 系の角運動 Ω の成分は

Ω = (0, 0, Ω)(4.13)である。
角速度 Ω で動く回転座標系 GXY Z における運動方程式をオイラーの方程式で表すとd0 L=Ω×L=Ndt(4.14)である。
ここで N は物体にはたらく力のモーメントである。この座標系に対する角運動量の成分は

Lx = L0x cos θ + L0z sin θ= (Cn − AΩ cos θ) sin θ(4.15)

Ly = L0y = Aθ̇Lz =(4.16)−L0x sin2 θ + L0z cos θ= AΩ sin θ + Cn cos θ= AΩ + (Cn − AΩ cos θ) cos θ(4.17)となり、ベクトル積 Ω × L の成分は(Ω × L)X= ΩY Lz − ΩZ LY = −AΩθ̇(Ω × L)Y= ΩZ LX − ΩX LZ(4.18)= −Ω(Cn − AΩ cos θ) sin θ(Ω × L)Z= ΩX LY − ΩY LX = 0(4.19)

R + F = (0, F, R) = (0, µM g, −M g)(4.20)となる。また、力はであり、接触点の座標は(

X P = (XP , YP , ZP ) =

)

dh(θ)
, 0, −h(θ)
dθ(4.21)

で与えられる。
これらを (4.14) に代入すると運動方程式はd[(Cn − AΩ cos θ) sin θ] − AΩθ̇ = −F ZP
dtAθ̈ + Ω(Cn − AΩ cos θ) sin θ = −RXPdAΩ̇ + [(Cn − AΩ cos θ) cos θ] = F XP
dt
19

(4.22)(4.23)(4.24)

これらの式は sin θ,cos θ について非線型となるので、このままでは解けない。そこで近似を使って解くことにする。
この方程式の y 成分に注目する。ここで Ω2 が十分大きいとすると、Ωを含む左辺第 2 項と第 3 項は右辺の −RXP を凌駕する。またゆで卵が起き上がっていく時間スケールは、回転で決まる短いものではなく、摩擦で決まる非常に長いものである。その時間スケールのみに注目すると、方程式 (4.23) は近似的に (Cn − AΩ cos θ) cos θ=0 となり、sin θ 6= 0 では

Cn = AΩ cos θ(4.25)となる。この近似を gyroscopic balance という。これを用いると運動方程式の X 成分と Y 成分はそれぞれAΩθ̇ = F ZP

(4.26)

AΩ̇ = F XP

(4.27)という簡単な式になる。(4.27) を (4.26) で割ると
dh
Ω̇
XP
dθ
ḣ
=
θ̇ = − dθ
=−
Ω
ZP
h dt
h

(4.28)という１階の変微分方程式になり、両辺を t で積分するとΩh = ecAΩh = Aec = J(4.29)(4.30)

が導き出される。この J を Jellet 定数という。即ち、回転の角速度 Ω と重心の高さ h の積はジェレット定数と呼ばれる運動の定数になる。Ω も h
も一定のまま定常的な歳差運動をするが 、摩擦があると Ω が減少し 、そのため h が増加するのである。
さらに (4.26) の第 1 式 AΩθ̇ = F ZP と (4.1) の第２式 ZP = −h(θ) を援用すれば (4.30) からJ θ̇ = −F h2 (θ)(4.31)を得る。これは物体の形 h(θ) が与えられれば直ちに積分できる式である。
すなわち、θ が t の関数として求まるということである。

20

一つの例として回転楕円体

z 2 x2 + y 2
+
=1
a2
b2

(4.32)の形をした物体を考えると、これを与える関数 h(θ) は

h2 (θ) = a2 − (a2 − b2 ) sin2 θ(4.33)であることが確かめられる。これが一様な密度を持っているとすると

A=

M 2(a + b2 ),52
C = M b2
5

(4.34)µM ga2b2(1 − (1 − 2 ) sin2 θ)Ja(4.35)

このとき (4.31) はθ̇ = −
これを積分するとbµM ga2aarctan( tan θ) = −(t − t0 )baJ(4.36)(to は積分定数) を得る。J の符合を考えると
atan θ = − tan µq(t − t0 )b(4.37)となる。ただしq = M gaba−b|a − b||J|(4.38)

とおく。細長い楕円体 (a > b すなわち g > 0) の場合 θ が π2 まで変化するのに要する時間は

∆t =

π
2µq

(4.39)となる。楕円体でも θ = 0 から θ = π2 になるまでには ∆t = π2 |µq| の時間を要する。
なお、θ = π4 のときの Ω を Ω∗ とすると Jellet 定数 J はJ = AΩ∗ (a2 + b2 1
)2
2

と書ける。

21

(4.40)

第5章

運動方程式の数値積分による
回転直立の再現

x 、y 、z 成分で実験室系の座標、ξ, η, ζ 成分で剛体系の座標を表示する。
以下の記号を与える• M : 質量• Ii : i 軸の主慣性モーメント• F : 外力
• N : 外力の重心のまわりのモーメント• R : 重心の座標
ここで Ṙ を V と書く。eξ , eζ , eη を主軸方向の単位ベクトルとする。
ei · ej = δij (i, j = ξ, η, ζ) を満たさなければならない。また eξ × eη =eζ ,eη × eζ = eξ ,eζ × eξ = eη が成り立つ。

ω の剛体系での成分を ωξ , ωη , ωζ とすると
ωξ = eξ · ω(5.1)ωη = eη · ω(5.2)= eζ · ω

(5.3)

ωζ

ω = ωξ eξ + ωη eη + ωζ eζ(5.4)

が成り立つ。
また以下の等式が成り立つ。

ω × eξ = ωζ eη − ωη eζ

(5.5)

ω × eη = ωξ eζ − ωζ eξ

(5.6)

ω × eζ = ωη eξ − ωξ eη

(5.7)

22

したがって剛体の運動は以下の１階の連立常微分方程式で決定される。

e˙ξ = ω × eξ

(5.8)

e˙η = ω × eη(5.9)

= ω × eζ(5.10)e˙ζω̇ = f (ω, eξ , eη , eζ ) + N (R, V , eξ , eη , eζ , ω)
(5.11)
Ṙ = V
1V̇ =F (R, V , eξ , eη , eζ , ω)M

(5.12)(5.13)(5.14)N は重心のまわりの外力のモーメント、F は外力である。以下に (5.11)の具体的な計算手順を示す。f の定義も同時に示す。
ωξ = ω · eξ

(5.15)

ωη = ω · eη

(5.16)= ω · eζ(5.17)

Nξ = N · eξ

(5.18)

Nη = N · eη

(5.19)= N · eζ(5.20)ω˙ξ = fξ ωη ωζ + Nξ

(5.21)ω˙η = fη ωζ ωξ + Nη

(5.22)ω˙ζ

(5.23)

ωζNζ

= fζ ωξ ωη + Nζω̇ = ω˙ξ eξ + ω˙η eη + ω˙ζ eζ+ ωξ e˙ξ + ωη e˙η + ωζ e˙ζ

(5.24)

但しfξ =Iη − IζIζ − IξIξ − Iη, fη =, fζ =
Iξ
Iη
Iζ

である。

23

(5.25)eξ , eη , eζ が従うべき条件は 6 つあり|eξ | = |eη | = |eζ | = 1, ex · ey = ey · ez = ez · ex = 0(5.26)

これらの条件は、初期状態でみたされていれば 、自動的にみたされる理由は下記の等式が成り立つからである。

d(eξ · eξ ) =dt=d(eξ · eη ) =dt=

2e˙ξ · eξ = 2(ω × eξ ) · eξ

(5.27)

2ω · (eξ × eξ ) = 0(5.28)e˙ξ · eη + eξ · e˙η = (ω × eξ ) · eη + eξ · (ω × eη )
ω · (eξ × eη + eη × eξ ) = 0(5.29)[力学変数の個数についての考察]
並進 3 ＋回転 3 の 6 変数が時間についての 2 階の微分方程式に従うということは 12 変数が時間についての１階の微分方程式に従うということである。
上の定式化では６ベクトルの合計１８成分即ち１８個の力学変数が時間の１階の微分方程式にしたがっている。18 − 12 ＝ 6 個の拘束条件が自動的に満たされているから自由度の数はあっている。

24

第6章

6.1

数値計算の結果

軸対称物体のオイラー法による解

まず物体に外力が全く働かない場合の軸対称物体の解を求めた。これは基本的なことなので割と簡単にプログラミングできる。付録としてプログ
ラムリストを載せる。オイラー法を簡単に式にして表すと

dω(t)
dt

= f (t)(6.1)

ω(t + dt) = ω(t) + ∆t ·dω(t) + θ((∆t)2 )
dt

(6.2)

で ω の値が導き出される。オイラー法を用いた数値計算の結果は以下のようになった。

まず初期条件を与える。
dt = 0.001 秒とした。これは 1000 個の値をプロットしたことを意味する。
他の条件を次ページのように与える。

25

記号

記号の意味

数値

変数名

R0ζ
Pξ
Pη
Pζ
ρ
Pξ ’
Pη ’
Pζ ’
Iξ
Iη
Iζ
ωξ0
ωη0
ωζ0
ωξ
ωη
ωζ
eξx
eξy
eξz
eηx
eηy
eηz
eζx
eζy
eζz

ζ 座標の長さ
ξ 座標の形の比η 座標の形の比ζ 座標の形の比密度
ξ 座標のベクトルの比
η 座標のベクトルの比
ζ 座標のベクトルの比ξ 軸の慣性モーメントη 軸の慣性モーメントζ 軸の慣性モーメント最初の ξ 軸の角速度
最初の η 軸の角速度
最初の ζ 軸の角速度
ξ 軸の角速度η 軸の角速度ζ 軸の角速度ξ 軸の x 方向の角速度単位ベクトル
ξ 軸の y 方向の角速度単位ベクトル
ξ 軸の z 方向の角速度単位ベクトル
η 軸の x 方向の角速度単位ベクトル
η 軸の y 方向の角速度単位ベクトル
η 軸の z 方向の角速度単位ベクトル
ζ 軸の x 方向の角速度単位ベクトル
ζ 軸の y 方向の角速度単位ベクトル
ζ 軸の z 方向の角速度単位ベクトル

2
3.1
2.9
5
1
3
2
1

R0ζ
radius − ratio − xradius − ratio − yradius − ratio − zdensityangular-velocity-ratio-x
angular-velocity-ratio-yangular-velocity-ratio-z
moi1
moi2
moi3
omg01
omg02
omg03
omg1
omg2
omg3
e1x
e1y
e1z
e2x
e2y
e2z
e3x
e3y
e3z

26

0
1
0
0
0
-1
1
0
0

前ページで掲載しなかった値をあげる。

R0ξ =
R0η =

Pξ
× R0ζ
Pζ
Pη
× R0ζ
Pζ

(6.3)
(6.4)

楕円体の体積 V はV=

4 Pξ Pη R0ζ 3
π
3
Pζ 2

(6.5)

となり慣性モーメントは

Iξ =
Iη =
Iζ

=

1
ρV (R02η + R02ζ )5
1
ρV (R02ζ + R02ξ )5
1
ρV (R02ξ + R02η )
5

(6.6)
(6.7)
(6.8)

角速度の初期値は
√

ωξ =

Pξ0

ωη =

Pη0

ωζ

Pζ0

√
√

=

|ωo |2Pξ 2 Pη 2 Pζ 2

(6.9)

|ωo |2Pξ 2 Pη 2 Pζ 2

(6.10)|ωo |2
Pξ Pη 2 Pζ 2

(6.11)

2となる。
これらの初期値を用いてオイラー法による数値計算を行う。

27

omega-1
omega-2
omega-380

angular velocity (1/sec)

60
40
20
0
-20
-40
-60
-80
0

0.2

0.4

0.6

0.8

time (sec)
図 6.1: 横軸 時間 t[sec] 、縦軸 角速度ベクトルの比 [1/s] を表したグラフ

図 6.1 は、外力の働かない剛体の楕円体の回転を計算した結果である。ここで、主軸方向の半径の比は、ξ 軸：η 軸：ζ 軸＝ 3.1 : 2.9 : 5.0 とした。
角速度の初期値は毎秒１０回転である。図の横軸は時刻（ 秒）、縦軸は角
速度ベクトルの３つの主軸方向の成分の値である（１／秒）。実線は ξ 軸成分、長い破線は η 軸成分、短い破線は ζ 軸成分である。図で示した計算例と同じように、初期条件をどのように変化させても角速度の主軸成分は周期的な運動を行うことが確かめられた。つまり外力の働かない軸対称物
体は、定常的な周期運動を続けることになり物体は直立せずに重心の高さも変化しない。Moﬀatt 、下村の理論からもこうなるはずであると推論できたが 、それを確かめることができた。

28

1

6.2

ルンゲクッタ法により求めた回転楕円体の直立運動

次に 4 次のルンゲクッタ法を用いて平らな摩擦のある床の上で回転する
剛体の運動の解を求めた。
先程オイラー法で用いた記号はそのまま使う。新たな条件を以下に記す。物体にかかる力は、摩擦力、重力、垂直効力があげられる。垂直抗力にはポテンシャルエネルギーを仮定すると、抗力 N は

N =0(z À 0)(6.12)

N = −kz(z < 0)(6.13)となる。ここで z は剛体表面の最下点の座標である。剛体表面の最下点を床との接点であると定義する。重心の速度を V 、重心から見た接点の位置ベクトルを r とする。このとき接点の速度、摩擦力はそれぞれ

v =VG+ω×rv
F = −µ N|v|(6.14)

Rξ = Rη = 1.8cm(6.16)

Rζ = 2.4cm(6.17)(6.15)である。楕円体の形状をで与え |ω| = 50(回転/秒) 、摩擦係数を 0.5 とする。
初期値をこれらで与え 4 次のルンゲクッタ法での数値積分の結果を以下のグラフで表示する。

29

400

ω1
ω2
ω3

300

angular velocity[1/sec]

200

100

0

-100

-200

-300

-400
0

1

2

3

4

time[sec]
図 6.2: 横軸 時間 t[sec] 、縦軸 角速度ベクトルの各成分 [1/sec] を表したグラフ

図 6.2 は横軸が時刻 t(秒) 、縦軸が 3 つの角速度ベクトルの主軸成分の比 (1/sec) の値である。ω3(式では ωζ と書いた) が下がると同時に楕円体が立ち上がりはじめ、楕円体が直立すると ω3 は一定値となる。ω1(ωξ ) 、
ω2(ωη ) は、波形となり楕円体が直立していくにつれ徐々に小さくなり、楕
円体が直立すると 0 に近い値になる。

30

5

2.5

2.4

2.3

z[cm]

2.2

2.1

2

1.9

1.8
0

1

2

3

4

time[sec]
図 6.3: 横軸 時間 t[sec] 、縦軸 z[cm] を表したグラフ

図 6.3 は横軸に時刻 t(秒), 縦軸は重心の高さ z(cm) の値である。図は ζ軸のみのグラフである。楕円体の重心は 2 秒から 3 秒の間だけ上がること
がわかる。

31

5

6.2e+06

6e+06

5.8e+06

energy[erg]5.6e+06

5.4e+06

5.2e+06

5e+06

4.8e+06

4.6e+06
0

1

2

3

4

time[sec]
図 6.4: 横軸 時間 [sec] 、縦軸 energy[erg] を表したグラフ

図 6.4 は横軸は時刻 t(秒) 、縦軸はエネルギー E(erg) の値である。楕円体が立ち上がるとき、摩擦によりエネルギーが減少する。このエネルギーはほとんど回転の運動エネルギーである。並進運動のエネルギーや、重力の位置エネルギーの変化量は 3 桁以上小さい値である。

325

20000

Lx
Ly
Lz

18000
16000
14000

L[g.cm2/s]

12000
10000
8000
6000
4000
2000
0
-2000
0

1

2

3

4

time[sec]
図 6.5: 横軸 時間 t[sec] 、縦軸 角運動量 L[g·cm2 /s] を表したグラフ

図 6.5 は横軸は時刻 t(秒) 、縦軸は角運動量 L の値である。鉛直上向きに立ち上がるとき、z 軸の角運動量が摩擦で減少している。x 軸、y 軸方向の成分は、常にほぼ 0 である。即ち角運動量は、回転直立の前後を通じて、
常にほとんど 鉛直方向を向いていることがわかる。

33

5

328

326

ω[rad/sec]

324

322

320

318

316

314
0

1

2

3

4

time[sec]
図 6.6: 横軸 時間 t[sec] 、縦軸 角速度の大きさ ω[rad/s] を表したグラフ

図 6.6 は横軸は時刻 t(秒) 、縦軸は ξ 軸の角速度の大きさ ω(rad/s) の値である。ζ 軸の角速度は楕円体が立ち上がると増える。これは ζ 軸のまわりの慣性能率が小さいためである。

34

5

第7章

結論

1, オイラー法による剛体の自由回転の数値計算の結果
初期条件をどのように変えても楕円体は剛体型でみて周期的な歳差運動を繰り返す。物体に働く力、とりわけ摩擦力を考慮しない場合、回転楕円体が直立することはない。このプログラムは自力で完成させたので付録に掲載した。

2, ルンゲクッタ法による数値計算の結果
4 次のルンゲクッタの場合、初期条件の値をうまく選んでやると楕円体が直立することがわかった。楕円体の対称軸を床と水平に回転を与えた場
合は直立しない。この場合は物体に摩擦力が働かないからである。また対称軸以外の角速度を少しでも大きい数字にしてしまうと、物体は床から離れる (跳びはねる) という結果になった。
今回の研究で、私はルンゲクッタ法によるプログラムをつくることができず、代わりに指導教員の作成したものを使用した。

35

参考文献表[1] H.K. Moﬀatt and Y. Simomura: Spiningeggs-a paradox resolved,
Nature 416,385-386(2002).
[2] H.K. Moﬀatt: Euler’s disk and its finite-time singularity,
Nature 404 ,833-834(2000).
[3] 戸田盛和:「回転する卵はなぜ直立する」科学 (岩波書店),
72 ,932-939 (2002).
[4] 今井 功:「魔法のコマ」日本物理学会誌,
8, 288-293(1953).
[5] 伏見康治:「逆立ゴ マ」数学セミナー (日本評論社),
3, 28-32(1967); 4, 38-42 (1967). 6, 32-34 (1967).
[6] 小井出昭一郎:「力学」 岩波書店,
pp.90-165 (1987).
[7] V.D. バーガー, M.G. オルソン : 力学「新しい視点に立って」,
pp.221-223 (1997).
[8] 下村裕: 「立ち上がる回転ゆで卵の解」パリティ,
18,52-56 (2003).

36

謝辞

本論文を作成するにあたり、田嶋 直樹 先生には 、終始御厚いご 指導をして頂いたことに誠に感謝し 、お礼申し上げます。また 鈴木 敏男 先生、林 明久 先生にも本研究及び日常的なことにおいても、実に丁寧なご指導、お世話をして頂きました。
本研究に対してご意見を頂いた、多くの物理工学科の先生方のみなさまにもお礼申し上げ、謝辞の言葉とさせて頂きます。

37

付録 : プログラムリスト
以下は、外力の働かない場合の剛体の回転運動をオイラー法で数値的に解く C 言語のプログラムのソースリストである。このプログラムは、卒業研究の一環として自作したものである。変数名などの説明はコード 中のコメントとして与えられている。
#include<stdio.h>#include<math.h>main(){
kaiten4();}
int kaiten4(){double radius_ratio_x=3.1;
double radius_ratio_y=2.9;
double radius_ratio_z=5;
double R00=2.0 ;double density=1;/* relative length of prinxipal axis *//* radius for shperical shape [cm] *//* density [g/cm^3] */double angnlar_velocity_ratio_x=3; /* relative size of initial ang. vel. */double angnlar_velocity_ratio_y=2;
double angnlar_velocity_ratio_z=1;
double moi1, moi2, moi3 ; /* moment of inertia */double omg01, omg02, omg03; /* initial angular velocity vector, B-frame */double omg1, omg2, omg3 ; /* angular velocity vector, B-frame */double R0x,R0y,R0z,R0 ;double omg0;
double a1,a2,a3 ;double dt,t,pi, b,V;
double e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z;
double e1dx,e1dy,e1dz,e2dx,e2dy,e2dz,e3dx,e3dy,e3dz;
double L1x,L1y,L1z,L2x,L2y,L2z,L3x,L3y,L3z,Lx,Ly,Lz;
pi=4*atan(1.0); fprintf(stderr,"check : pi = %20.16f\n",pi);dt=0.0001;/* time step size [sec] */omg0=20*pi; /* size of initial ang. vel. [radian/sec] */R0=R00/pow(radius_ratio_x*radius_ratio_y*radius_ratio_z,1.0/3.0);
R0x=R0*radius_ratio_x; R0y=R0*radius_ratio_y; R0z=R0*radius_ratio_z;
fprintf(stderr,"R0=(%f %f %f)\n",R0x,R0y,R0z);
V=4*pi*R0x*R0y*R0z/3; /* volume of the rigid body [cm^3] */moi1=density*V*(R0y*R0y+R0z*R0z)/5; /* moment inertia [g cm^2] */moi2=density*V*(R0z*R0z+R0x*R0x)/5;
moi3=density*V*(R0x*R0x+R0y*R0y)/5;
fprintf(stderr,"MOI=(%f %f %f)\n",moi1,moi2,moi3);
b=1/sqrt(angnlar_velocity_ratio_x*angnlar_velocity_ratio_x+angnlar_velocity_ratio_y*angnlar_velocity_ratio_y38

+angnlar_velocity_ratio_z*angnlar_velocity_ratio_z);
omg01=omg0*angnlar_velocity_ratio_x*b;
omg02=omg0*angnlar_velocity_ratio_y*b;omg03=omg0*angnlar_velocity_ratio_z*b;omg1=omg01;
omg2=omg02;
omg3=omg03;
e1x=1, e1y=0, e1z=0;
e2x=0, e2y=-1, e2z=0;
e3x=0, e3y=0, e3z=1;
for(t=0;t<=1;t=t+dt) {L1x=moi1*omg1*e1x;L2x=moi2*omg2*e2x;L3x=moi3*omg3*e3x;
L1y=moi1*omg1*e1y;L2y=moi2*omg2*e2y;L3y=moi3*omg3*e3y;
L1z=moi1*omg1*e1z;
L2z=moi2*omg2*e2z;
L3z=moi3*omg3*e3z;
Lx=L1x+L2x+L3x;
Ly=L1y+L2y+L3y;
Lz=L1z+L2z+L3z;/*printf("%f %f %f %f\n",t,Lx,Ly,Lz);*/a1=(moi2-moi3)*omg2*omg3/moi1;
a2=(moi3-moi1)*omg3*omg1/moi2;
a3=(moi1-moi2)*omg1*omg2/moi3;
e1dx=dt*(omg3*e2x-omg2*e3x);e1dy=dt*(omg3*e2y-omg2*e3y);
e1dz=dt*(omg3*e2z-omg2*e3z);
e2dx=dt*(omg1*e3x-omg3*e1x);e2dy=dt*(omg1*e3y-omg3*e1y);e2dz=dt*(omg1*e3z-omg3*e1z);e3dx=dt*(omg2*e1x-omg1*e2x);e3dy=dt*(omg2*e1y-omg1*e2y);e3dz=dt*(omg2*e1z-omg1*e2z);
printf("%f %f %f %f %f %f %f %f %f %f %f %f\n",t,omg1,omg2,omg3,
e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z);
omg1=omg1+dt*a1;
omg2=omg2+dt*a2;

39

omg3=omg3+dt*a3;e1x=e1x+e1dx;
e1y=e1y+e1dy;
e1z=e1z+e1dz;
e2x=e2x+e2dx;
e2y=e2y+e2dy;
e2z=e2z+e2dz;
e3x=e3x+e3dx;
e3y=e3y+e3dy;
e3z=e3z+e3dz;}}

40