卵を回すとなぜ立つか II

2005 年 2 月 10 日
福井大学 工学部 物理工学科

11 年度入学 16 番 大森英胤

目次
序章

2

第 1 章 剛体の力学

3
3
3
4
5
6
6
6
7
9
9
11
12

1.1

1.2
1.3

1.4

1.5

剛体の定義とつりあい . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.1.1 剛体の定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.2 剛体のつりあい . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
固定軸まわりの剛体の回転運動 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
慣性モーメントと慣性乗積 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 慣性モーメント . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.2 慣性乗積 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.3 慣性楕円体 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
剛体回転のオイラー方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4.1 慣性テンソル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4.2 剛体回転のオイラー方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ポアンソーの表現法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 2 章 軸対称物体の回転直立

2.1
2.2
2.3
2.4
2.5

–Moffatt 下村理論の紹介–
座標系と接触点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ジェレット定数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ジャイロスコピック近似 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
回転楕円体の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Moffatt 下村理論のまとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 3 章 数値シミュレーションによる回転直立現象の再現

3.1
3.2
3.3
3.4

プログラムの改良と検証 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
重心位置の上昇の様子 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ジェレット定数の時間変化について . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ジャイロスコピック近似 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 4 章 結論

14
14
19
20
22
24
25
25
27
28
29
30

関連図書

31

謝辞

32

付録

Program List

33

1

序章
ゆで卵や碁石などの剛体を机の上において，勢いよく回転させると立ち上がる。この現象の不思議なと
ころは，重心が重力に逆らって上昇する，ということである。この現象は 100 年以上前から知られていたが，これまで解明されていなかった。2002 年に Moffatt と下村によって卵のような軸対称物体の回転に伴う立ち上がりの現象を説明する理論が発表された。
この理論によると, 回転する剛体の運動方程式に対して C を対称軸，A をそれ以外の軸まわりの慣性モーメント，Ω, n をそれぞれ鉛直軸，対称軸方向の角速度，θ を鉛直軸と対称軸のなす角度としたとき

Cn = AΩ cos θ という近似を用いると，運動方程式を積分することができる。この近似をジャイロスコ
ピック近似とよび，回転が速い場合においては, 自然な仮定であるが，定量的な検討はまだなされたことがない。
軸対称物体の立ち上がり現象と関係の深いものとして，逆立ちごまがある。逆立ちごまのような偏心球の回転に伴う立ち上がりの現象は，1952 年 Braams および Hugenholtz によって厳密解が得られていて，そこにはジェレット定数とよばれる運動定数が存在することがわかっている。Moffatt と下村の理論では一般の軸対称物体に対しても, このジェレット定数が保存される。上述の記号および重心の高さ hを用いると J = AΩh と表される。このジェレット定数を用いると運動方程式は非常に簡単な形で書くことができ, 容易に積分できるようになる。
本卒業研究では剛体の並進 · 回転運動を表す微分方程式を 4 次のルンゲクッタ法を用いて解き，上記の近似を用いずに軸対称物体の回転直立現象を再現し，その際にジャイロスコピック近似がどの程度の精度で成立しているのか や ジェレット定数の時間変化について調べた。

2

第 1 章 剛体の力学
この章の執筆にあたり，文献 [4, 5] を参考にした。

1.1

剛体の定義とつりあい

1.1.1

剛体の定義

剛体とは形が全く変化しない物体のことをいう。言い替えると，物体内の質点間の距離 |r i − r j | がい
つでも全て一定に保たれている質点系であるということができる。
剛体の 3 点に A,B,C と目印を付けると，この 3 点の位置が決まると剛体の位置も向きも決まる。まずA の位置を決めるには，xA , yA , zA の 3 個の座標を決めなければならない。次に B の位置を決めるには，AB の長さは決まっているので，極座標と同様に線分 AB の方向 θ と向き φ を決めればよい。C は ABを軸として回転する角 ψ によって決まる。すなわち剛体の運動は 6 個の変数によって記述されることが分かる。このことを剛体の自由度は 6 であるという。
一般に拘束条件が 1 個増えると自由度は 1 減るということができる。

3

1.1.2

剛体のつりあい

剛体がつりあいを保つということは，次の 2 つが成り立っていることである。

(1) 剛体の重心が静止している。
(2) 剛体の重心まわりで回転が起こっていない。
剛体の重心の位置をR, 剛体を質点の集まりとして見たときの各質点の位置と質量をそれぞれ r 1 , r 2 , · · · , m1 , m2 , · · ·とし，全質量を M とすると重心は次のように定義される。
R=

m1 r 1 + m2 r 2 + · · ·
M

(1.1)

(1.1) 式の両辺を時間で 2 回微分して変形すると
M

d2 R
dt2

= m1 r¨1 + m2 r¨2 + · · · =Xmi r̈ i

(1.2)

i

また i 番目の質点が受けている外力を F i とする。剛体内の質点が及ぼしあう内力は打ち消し合うことを考慮すると (1.2) 式はM

d2 R X
=
Fi
dt2
i

(1.3)

となる。すなわち (1) の条件はX

Fi = 0(1.4)

iとかくことができる。

2) の条件については重心まわりの角運動量 L0 が常に 0 であると言えればよい。すなわち
X
d X 0
d X 0
L =
r i × mi r˙0i =r 0i × F i = 0
dt i
dt i
i

(1.5)

であればよい。また重心も動かないので LG = 0 であり，

X

ri × F i = 0iとかくことができる。

(1.4) 式と (1.6) 式をそれぞれ成分に分けた 6 個の式が，つり合いの条件式である。

4

(1.6)

1.2

固定軸まわりの剛体の回転運動

剛体が固定した軸のまわりで回転するときは，p3 でも述べたように自由度は 1 である。
固定軸を z 軸にとり角速度を ω とする。また剛体を質点の集まりとみなすとして i 番目の質点の質量を mi ，軸からの距離を ri とする。質点の x 座標からの偏角を θi とすると

dθi
dt

ω=

(1.7)

であり，質点 i の運動量は mi ri ω ，軸 z に関する角運動量は mi ri2 ω である。したがって全角運動量 Lz は

Lz =

X

mi ri2 ω(1.8)

iとなる。ここでIz=

Xmi ri2i=X

mi (x2i + yi2 )

(1.9)

i

を固定軸 z まわりの剛体の慣性モーメントという。これを用いると，軸のまわりの角運動量は

Lz = Iz ωとなる。

5(1.10)

1.3

慣性モーメントと慣性乗積

1.3.1

慣性モーメント

剛体に固定した直交座標軸 ξ, η, ζ をとり，剛体を質点の集まりとしてみたときの質点の質量を mi , 位
置を (ξ, η, ζ) とする。ξ, η, ζ 軸まわりの慣性モーメントはそれぞれ

Iξ =

X

mi (ηi2 + ζi2 )

i

Iη =

X

mi (ζi2 + ξi2 )(1.11)iIζ =X

mi (ξi2 + ηi2 )iとなる。剛体を密度が ρ(ξ, η, ζ) の連続体として扱える場合には，慣性モーメントは次のように求めることができる。

ZZZ
ρ(ξ, η, ζ)(η 2 + ζ 2 )dξdηdζ

Iξ =
ZZZ

ρ(ξ, η, ζ)(ζ 2 + ξ 2 )dξdηdζIη =(1.12)ZZZIζ =ρ(ξ, η, ζ)(ξ 2 + η 2 )dξdηdζ

また M を剛体の質量として，慣性モーメントを

I = M k2q
とおいたときの k ，すなわち

1.3.2(1.13)

I
M のことを回転半径と言う。

慣性乗積

剛体に固定した座標系 ξ, η, ζ に関して次の量を慣性乗積という。

Iξη =

X

mi ξi ηi

i

Iηζ =

X

mi ηi ζi

i

Iζξ =

X

mi ζi ξii

この量は対称性がよい場合 (例えば (ξ, η, ζ) ⇒ (−ξ, η, ζ) が成り立つ場合) は全て 0 になる。

6

(1.14)

1.3.3

慣性楕円体

剛体に固定した座標軸 O − ξηζ に関する慣性モーメントと慣性乗積がわかっているとする。このとき

O を通る任意の軸 OX と ξ, η, ζ 軸とのなす角の cos をそれぞれ λ, µ, ν とする。
ζ
Pi
µ
O

X

θi Qi
υ

η

λ
ξ

剛体内の一点 Pi (ξi , ηi , ζi ) から OX におろした垂線を Pi Qi , OPi と OXi のなす角を θi , OPi = ri とす
ると

cos θi =1(λξi + µηi + νζi )ri(1.15)とかけるからPi Q2i= ri2 sin2 θi= ri2 − ri2 cos2 θi= (ξi2 + ηi2 + ζi2 ) − (λξi + µηi + νζi )2

ここで λ2 + µ2 + ν 2 = 1 であるから

Pi Q2i= (ξi2 + ηi2 + ζi2 )(λ2 + µ2 + ν 2 ) − (λξi + µηi + νζi )2= (ηi2 + ζi2 )λ2 + (ζi2 + ξi2 )µ2 + (ξi2 + ηi2 )ν 2
−2µνηi ζi − 2νλζi ξi − 2λµξi ηi

これに mi をかけて加えX
2
mi Pi Qi
i

=

X
i

−2

(ηi2 + ζi2 )λ2 +XX(ζi2 + ξi2 )µ2 +iµνηi ζi − 2iXνλζ iξi − 2

i

X

X

(ξi2 + ηi2 )ν 2i

λµξi ηiiとすると，これは OX まわりの慣性モーメント IOX となるが，ここで (1.11) 式と (1.14) 式を用いて書きかえるとIOX = Iξ λ2 + Iη µ2 + Iζ ν 2 − 2Iηζ µν − 2Iζξ νλ − 2Iξη ξη

(1.16)となる。

(1.16 ) 式の λ, µ, ν にそれぞれ ξ, η, ζ を代入して，IOX = 1 とするとIξ ξ 2 + Iη η 2 + Iζ ζ 2 − 2Iηζ ηζ − 2Iζξ ζξ − 2Iξη ξη = 1(1.17)となる。この式を満足する点 (ξ, η, ζ) は一つの楕円体面をつくる。これを慣性楕円体という。この楕円体と OX との交点を R(ξX , ηX , ζX ) とすると，OR = r とすると

ξX = rλ,ηX = rµ,7

ζX = rνであり，これらを (1.17) 式の ξ, η, ζ に代入すると

r2 (Iξ ξ 2 + Iη η 2 + Iζ ζ 2 − 2Iηζ ηζ − 2Iζξ ζξ − 2Iξη ξη) = 1(1.18)

ここで (1.18) 式の括弧の中は IOX に他ならないので

r2 IOX
⇒

IOX

= 1
1
=
r2

(1.19)となる。このことは原点 O をきめて，それに対応した慣性楕円体を描き O を通る任意の軸と慣性楕円体との交点を R とするとその軸のまわりの慣性モーメントはIOX = 1/ OR(1.20)

で与えられることを意味している。
慣性楕円体は，原点のとり方には依存するが，ξ, η, ζ 軸のとり方には依存せず，座標軸の方向を適当に選ぶことで (1.18) 式は Aξ 2 + Bη 2 + Cζ 2 の様な形にすることができる。このように選んだ ξ, η, ζ 軸を，慣性主軸に対する慣性モーメントを主慣性モーメントと言う。

8

1.4

剛体回転のオイラー方程式

1.4.1

慣性テンソル
ζ

ω

z

rii

η
y

O

x
ξ

図 1.1: ω で回転する剛体

図 1.1 のように O を固定点とし，実験室系 O − xyz と重心系 O − ξηζ をとる。剛体の点 O に角速度ω を与えると剛体上の任意の点 r i は時間 dt に ω, r i の両ベクトルに対して垂直な方向に (|r i | sin θi )ωdt
動くことになる。すなわちこの瞬間における r i の速度は
r˙i = ω × r i(1.21)と表すことができる。この (1.2) 式に mi を掛けて i について和をとると，剛体の運動量 P は

P

=

X

mi r˙i = ω ×iX

mi r ii= ω × MR

(1.22)とかける。ここで M は剛体の全質量で，R は剛体の重心である。
点 O のまわりの角運動量 L は r と P との外積でL=X

mi (r i × r˙i ) =X

imi r i × (ω × r i )(1.23)iとなる。ここで L の x 成分は

Lx

=

X

mi [yi (ωx yi − ωy xi ) − zi (ωz xi − ωx zi )]

i

=

X

mi [(yi2 + zi2 )ωx − xi yi ωy − zi xi ωz ]iとかける。重心系で表すと

Lξ =

X

mi [(ηi2 + ζi2 )ωξ − ξi ηi ωη − ζi ξi ωζ ](1.24)i

慣性モーメント Iξ =P

2
2i mi (ηi + ζi ) と慣性乗積 Iξη =Pi mi ξi ηi ,

Iζξ =

P

i mi ζi ξi を用いて (2.4) 式を書き換えるとLξ = Iξ ωξ − Iξη ωη − Iζξ ωζ

9(1.25)

となる。η ，ζ の成分も同様に

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

(1.26)

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

(1.27)となる。(1.25) 式，(1.26) 式，(1.27) 式をまとめて行列で書くと



 
Lξ
Iξ

 
 Lη  =  −Iξη
Lζ
−Iζξ

−Iξη
Iη
−Iηζ



−Iζξ
ωξ


−Iηζ   ωη Iζωζ

(1.28)となる。この (1.28) 式の右辺の 3 行 3 列の行列で表される量を慣性テンソルという。また慣性乗積は対称性がよくなるように座標系を選んでやると，すべて 0 にすることができることから，そのとき (1.28)
式は



  0
L0ξ
Iξ
 0  
 Lη  =  0
L0ζ
0

0
Iη0
0

 0 0ωξ
0   ωη0 Iζ0ωζ0(1.29)とすることができる。この Iξ0 , Iη0 , Iζ0 を主慣性モーメント，この時の座標軸の方向を慣性主軸という。

10

1.4.2

剛体回転のオイラー方程式

剛体上のベクトル B(t) の微小時間 dt における変位 dB を

(1) 重心系 O − ξηζ からみたときの B の変化(2) 実験室系 O − xyz からみたときの B の変化の 2 つの場合に分けて考えてみる。
(1) の場合の ξ, η, ζ 成分は dBξ , dBη , dBζ である。(2) の場合は (ω × B)dt となる。これらのベクトル和である dB を dt で割ると


dBξ
dB
=
+ (ω × B)ξ
dt ξ
dt


dB
dBη
=
+ (ω × B)η
(1.30)
dt η
dt


dBζ
dB
=
+ (ω × B)ζ
dt ζdtとなる。ξ, η, ζ 軸を慣性主軸に一致させたときの O に関する全角運動量 L は (1.29) 式

Lξ = Iξ0 ωξ ,Lη = Iη0 ωη ,Lζ = Iζ0 ωζ

で与えられ，時間微分すると

dLξ
dωξ
= Iξ0
,
dt
dt

dLη
dωη
= Iη0
,
dt
dt

dLζ
dωζ
= Iζ0
dt
dt

となる。(1.31) 式と L を (1.30) 式に用いると，ξ 成分については


dL
dωξ
= Iξ0
+ (ωη Lζ − ωζ Lη )
dt ξ
dt


Lζ
Lη
0 dωξ
= Iξ
+
−
ωη ωζ
dt
ωζ
ωη
dωξ
− (Iη0 − Iζ0 )ωη ωζ= Iξ0dtとなる。同様に η, ζ 成分は




dL
dt
dL
dt


= Iη0

dωη
− (Iζ0 − Iξ0 )ωζ ωξ
dt

= Iζ0

dωζ
− (Iξ0 − Iη0 )ωξ ωη
dt

η



ζ

(1.31)(1.32)となり角運動量に対する運動方程式

dL
dt

= N
X
=
(r i × F i )iに代入すると，重心系の各成分はXdωξ− (Iη0 − Iζ0 )ωη ωζ = (r i × F i )ξ
dt
i
X
dωη
− (Iζ0 − Iξ0 )ωζ ωξ = (r i × F i )η
Iη0
dt
i
X
dωζ
Iζ0
− (Iξ0 − Iη0 )ωξ ωη = (r i × F i )ζ
dt
i
Iξ0

これを剛体の回転に対するオイラー方程式という。

11(1.33)

1.5

ポアンソーの表現法

外力が 0 の場合のオイラー方程式は

dωξ
dt
0 dωη
Iη
dt
0 dωζ
Iζ
dt
Iξ0

= (Iη0 − Iζ0 )ωη ωζ
= (Iζ0 − Iξ0 )ωζ ωξ

(1.34)= (Iξ0 − Iη0 )ωξ ωη

となるが，(1.34) 式にそれぞれ ωξ , ωη , ωζ を掛けて加えると

Iξ ωξ

dωξ
dωη
dωζ
+ Iη ωη+ Iζ ωζ
=0
dt
dt
dt

(1.35)となる。この (1.35) 式を積分すると

1(Iξ ωξ2 + Iη ωη2 + Iζ ωζ2 ) = const
2(1.36)となる。この (1.36) 式の左辺は剛体の力学的エネルギーになっていて，力学的エネルギーが保存されていることが確認できる。
また，(1.34) 式にそれぞれ Iξ ωξ , Iη ωη , Iζ ωζ を掛けて加えるとIξ2 ωξdωξdωηdωζ+ Iη2 ωη
+ Iζ2 ωζ
=0
dt
dt
dt

(1.37)となり，この (1.37) 式を積分するとIξ2 ωξ2 + Iη2 ωη2 + Iζ2 ωζ2 = const(1.38)となる。Iξ ω = Lξ であり，(1.38) 式から角運動量が保存されることが確認できる。
慣性楕円体で点 O から引いた ω に平行な直線を考える。この直線と慣性楕円体との交点を W として，

OW = l とすると，点 W の座標 (ξW , ηW , ζW ) は
ξW = lωξ,ωηW = lωη,ωζW = lωζω(1.39)となる。これらは慣性楕円体上にあるので Iξ ξ 2 + Iη η 2 + Iζ ζ 2 = 1 を満たすので

l2(Iξ ω 2 + Iζ ω 2 + Iη ω 2 ) = 1
ω2(1.40)(1.36) 式の右辺を K とすると (1.40) 式は
ω=

√

2Kl

(1.41)となる。
点 W で慣性楕円体に接する平面上の点を r 0 = (ξ 0 , η 0 , ζ 0 ) とするとこれらはIξ ξW ξ 0 + Iη ηW η 0 + Iζ ζW ζ 0 = 1(1.42)

を満たすことが解析幾何学で知られている。これに (1.39) 式を代入すると

l(Iξ ωξ ξ 0 + Iη ωη η 0 + Iζ ωζ ζ 0 ) = 1
ωω √⇒ Lξ ξ 0 + Lη η 0 + Lζ ζ 0 = = 2Kl

12

(1.43)(1.44)

となる。左辺は L · r 0 であるが，これが一定であるということは，r 0 の乗る平面が L に垂直であることを示している。原点 O からこの接平面への垂線の長さを h とすると，L · r 0 = |L|h であり，これが√2K に等しいから√
2K(1.45)h=Lである。K, L は一定であることから h も一定である。このことからこの接平面は空間に固定した一定不変の面であるとわかる。これを不変面と言う。点 O を固定した慣性楕円体は，不変面に接しながら，滑べること無く，転がっていく。その回転の角速度は，各瞬間における OW の長さに比例する。
このとき，接点 W が不変面上に描く曲線をハーポールホード，慣性楕円体上にのこす軌跡をポール
ホードという。このように，慣性楕円体の運動として剛体の運動を表す方法をポアンソーの表現法と
いう。

L

ハーポール
ホード

ω

ζ

ポール
ホード

O

ξ
図 1.3: ポールホード
図 1.2: ハーポールホードとポールホード

13η

第 2 章 軸対称物体の回転直立–Moffatt 下村理論の紹介–
この章では文献 [1, 2, 3] をもとに Moffatt と下村による軸対称物体の回転直立の理論の詳細な解説を
行う。

2.1

座標系と接触点

Z

ζ

g
θ
X

O
h(θ)
R
P

ξ

図 2.1: 回転する軸対称物体
点 O を重心とする任意の軸対称物体が, 水平面上を滑りながら点 O のまわりに回転運動をしているとする。考えている瞬間における物体と水平面との接点を P とする。点 P において摩擦力が働かなければ, 運動方程式の解として定常的な歳差運動が存在することは, よく知られている [4,5]。そこで, さしあたり定常的な歳差運動が実現していて, したがって, 点 O の速度はゼロと見なして, 以下の議論を展開する。その上で, 非常に弱い摩擦力が働いた場合に, その効果により運動にどのような変化が起きるか, 具
体的には点 O が上昇するかどうか, を考察することにする。
まず座標軸の取り方を説明する。X, Y, Z 系 は 点 O を原点として水平面内で回転する座標系である。
Z 軸は鉛直上向きであり, X 軸は水平面内に, 点 P が XZ 平面にのるようにとる。Y 軸は X 軸と Z 軸とに垂直で, X, Y, Z 軸が右手系をなす向きにとる。X, Y 軸は点 P の運動に追随して水平面内で回転するが, 考えている瞬間におけるその角速度を Ω とする。ξ, η, ζ 系は物体に固定された座標系 (剛体系) であり, 点 O を原点とする。ζ 軸は物体の対称軸方向を, ξ, η 軸はそれに垂直な２方向を向くようにとる。
ξ, η, ζ 軸も右手系をなす。形状が軸対称であるため, ξ, η 軸の取り方に任意性はない。ξ, η, ζ 軸は, 物体の−→角速度ベクトル ω で回転している。Xp = OP は O を原点とした接点 P の位置ベクトルとする。この

14

とき物体が受ける力は・ 重力 M g = (0, 0, −Mg)
・ 垂直抗力 R
・ Y 方向への摩擦力 Fである。摩擦力が Y 方向を向いていることは, 物体の対称軸 (ζ) 軸の傾き (図中の θ) の変化が, 回転角速度に比べて無視できるほど小さいという仮定の帰結である。
物体の角速度ベクトル ω の各成分は

ωξ= −Ω sin θ

ωη

= θ̇

ωζ

= n

(2.1)

となる。
また水平面内回転系での各成分は

ωX= −ωξ cos θ + ωζ sin θ= −Ω sin θ cos θ + n sin θ= (n − Ω cos θ) sin θωY= ωη = θ̇

(2.2)

ωZ

= ωζ cos θ + ωξ sin θ= n cos θ + Ω sin2 θ

ただし n は ω の ζ 軸成分で

n = ψ̇ + Ω cos θである。ただし ψ̇ は図 (2.2) に示すような関係にある。

Z

Ω θ

ζ

ω

.

Ωcosθ

ψ

X

O

h(θ)

R

ξ

P
図 2.2: 2 つの系での ω と n

15

主慣性モーメントをIξ = Iη = A, Iζ = C とすると，角運動量 L の各成分はLξ

= −AΩ sin θ

Lη

= Aθ̇

Lζ

= Cn

(2.3)

となる。また

LX= Lξ cos θ + Lζ sin θ

LY

= Lη

LZ

= −Lζ sin θ + Lζ cos θ(2.4)(2.3) 式を (2.4) 式に代入して水平面内回転系にすると L は
L = ((Cn − AΩ cos θ) sin θ, Aθ̇, AΩ sin2 θ + Cn cos θ)

(2.5)

となる。
この L の時間に対する変化は，オイラー方程式∂L+ (ω × L) = XP × R + F
∂t

(2.6)

で与えられる。ここで F は点 P にはたらく摩擦力である。
図 (2.3) は角度 θ が微小角 ∆θ 変化したときを描いてある。X’,Z’ は変化後の水平面内回転系の座標軸
で h’ は変化後の重心の高さである。図 (2.3) の点 P 付近を拡大したものが図 (2.4) である。ここで線分

LN に着目するとLN ' −XP ∆θ

(2.7)

であり，また

LN= OL − ONh0= h−cos(∆θ)である。∆θ ' 0 のとき，cos(∆θ) = 1 となることから

LN ' h − h0 = −∆h

(2.8)

ゆえに

XP =

dh
dθ

(2.9)

となる。これを用いると位置ベクトル X P = (XP , 0, ZP ) の成分は

XP =

dh
,
dθ

ZP = −h(θ)(2.10)となる。
また，R + F = (0, F, R) であるから

XP × (R + F ) = (−ZP F, −RXP , F XP )16

(2.11)

Z’ Z

ζ

∆θ
θ
X’
X

O
h(θ)

h’(θ)

ξ
P

-Xp

図 2.3: 微小な角度変化 1

N
∆θ
-Xp
図 2.4: 微小な角度変化 2

17Lとなり (2.6) 式を用いるとオイラー方程式の各成分はd[(Cn − AΩ cos θ) sin θ] − AΩθ̇ = −ZP F
dt
Aθ̈ + Ω(Cn − AΩ cos θ) sin θ = −RXPdAΩ̇ + [(Cn − AΩ cos θ) cos θ] = F XP
dtと表される。

18(2.12)(2.13)(2.14)

2.2

ジェレット定数

図 2.5: 逆立ちごま

軸対称物体の回転直立の現象に関係のあるものとして，逆立ちごまがある。逆立ちごまとは球の上部を切り取って，棒を付けたもので，偏心球となっている。これを回転させると重心が上昇し，上下逆になり心棒でまわり続ける (図 2.5)。これについては 1952 年にブラームスとヒュヘンホルツによって理論的に解明されている。それによると逆立ちするためには・こまと接触面との摩擦が不可欠である。
・接触点がすべって力学的エネルギーが減少しなければならない。
ということが分かっている。このようにエネルギーが保存しない系であるにもかかわらず，ひとつの運動定数が存在する。その運動定数は，角運動量ベクトル L と接触点の位置ベクトル XP との内積としてJ = −L · X P(2.15)

で定義される量である。
任意の軸対称物体に対するジェレット定数の時間微分を計算するとJ˙= −L̇ · X P − L · X˙ P
= (Cn − AΩ cos θ)XP2

d sin θdt XP

(2.16)arとなる。図 2.6 のように軸対称物体の形が半径 r の球であるとき重心が球の中心から a だけずれている

θ

h(θ)
: 重心

図 2.6: 球の場合の h(θ)

とするとh(θ) = r − d cos θ(2.17)dh= a sin θdθ(2.18)

となるので

XP =

となることから，J˙ = 0 となり J は運動定数となる。
しかし，一般の軸対称物体の場合は J は運動定数とはならない。

19

2.3

ジャイロスコピック近似

軸対称物体の立ち上がり現象を考えるとき，回転が遅いと立ち上がらないので回転が速い場合について考える。オイラー方程式の Y 成分はAθ̈ − AΩ2 cos θ sin θ + CnΩ sin θ = −RXP

(2.19)である。Ω2 が大きいとすると，Ω が含まれている左辺の第 2 項と第 3 項は右辺の −RXP よりも十分に大きいとすることができる。次に時間スケールに着目すると軸対称物体が立ち上がる時間は回転で決まる短いものではなく摩擦で決まる長いものであるので |θ̈|  Ω とすることができる。よって (2.19) 式は近似的に

Cn − AΩ cos θ = Ω sin θ(2.20)とでき，立ち上がるまでの過程すなわち，sin θ 6= 0 では

Cn = AΩ cos θ(2.21)となる。これをジャイロスコピック近似とよぶ。
このジャイロスコピック近似を取り入れると (2.16) 式は J˙ = 0 となり J は運動定数となる。また (2.5)
式で与えられる角運動量 L は

L = (0, Aθ̇ AΩ)(2.22)となる。したがって (2.15) 式から

J= −(−h(θ) · AΩ)= AΩh(2.23)となる。
また，オイラー方程式の X 成分，(2.12) 式はAΩθ̇ = F ZP

(2.24)

AΩ̇ = F XP

(2.25)XPḣΩ̇=θ̇ = −Ω
ZP
h(2.26)となる。そして Z 成分，(2.14) 式はとなり (2.24)，(2.25) 式からが導かれる。この (2.26) 式を積分すると Ωh= 一定 となる。また (2.8) 式，(2.23) 式，(2.24) 式より任意の軸対称物体の θ に対するジャイロスコピック近似を用いた方程式はJ θ̇ = −F h2 (θ)(2.27)という簡単な 1 階の微分方程式となる。
ここで軸対称物体と床との接触点 P での摩擦係数を µ，すべり速度を Vp としてクーロン摩擦力を仮定すると，摩擦力 F は

F = −µM g

20

Vp
|Vp |

(2.28)

で与えられる。

VP は剛体の角運動量 ω と接触点 P の位置ベクトル XP との外積で表すことができる。すなわち
VP= ω × XP=(Ω sin2 θ + n cos θ)dh+ (n − Ω cos θ) h(θ) sin θdθ(2.29)である。β(θ) = (sin2 θ + (A/C) cos2 θ)−1/2 として (2.21) 式，(2.23) 式をもちいて (2.29) 式の VP を書き換えると次のようになる。

VP =J −3 −1 d(βh)
β h
A
dθ

(2.30)

F と VP の関係は分かっているものとし，h(θ) も幾何学的な考察から分かっているとするとき (2.27)
式を積分することで θ は t の関数として求まる。

21

2.4

回転楕円体の場合

ここでは密度が一様な回転楕円体の場合に (2.25) 式を積分してみる。回転楕円体の対称軸方向の半径を a，それと垂直な方向の半径を b，質量を M ，密度を ρ とする。この回転楕円体は

a2 (x2 + y 2 ) + b2 z 2 = a2 b2と表すことができる。このときの対称軸 (ζ 軸) まわりの慣性モーメントはZIζ =ρ(ξ 2 + η 2 )dξdηdζ
V
Z
ZZ= 2ρ dη dζ ξ 2 dξ

Z b
= 2ρ

r2 dr

0

Z b
= 2ρ

Z π

Z 2π
sin θdθ

0

r4 dr

Z π

0

dφ(r sin θ cos φ)20

sin3 θdθ

Z 2π0cos2 φdφ0= 2M b2 /5(2.31)

同様にして= Iη = M (a + b)2 /5

Iξ(2.32)となる。楕円体上の点 (x0 , 0, z0 ) において

X
b /x
2

0

(x , 0, z )

h(θ

)

0

-a

0

θ

o

a

a /z
2

0

Z

-b
図 2.7: 重心の高さ h(θ)

 z 20 x 20=1(2.33)zx
dz + 2 dx = 02ab(2.34)a+bであるが，両辺を微分するととなるから点 (x0 , 0, z0 ) における接線の式はx0z0(z − z0 ) + 2 (x − x0 ) = 0
a2b z 2  x 2
z0
x0
0
0
⇒ 2z + 2 x =
+
=1
a
b
a
b

(2.35)

であり, この接線と原点 O との距離が h であるから図 2.7 からわかるように

h(θ) =b2a2
cos θ =
sin θ
z0
x0
22

(2.36)

である。(2.33) 式と (2.36) 式から

h(θ) =

p
a2 cos2 θ + b2 sin2 θ(2.37)となる。このときの P 点すべり速度は
VP = −

J(a2 − b2 ) sin 2θ
4Ah2


(2.38)と簡単になり，クーロン摩擦力を仮定する式を積分すると

tan θ = (a/b) tan µq(t − t0 )となる。ただし q = M gab(a − b)/|a − b||J| であり，t0 は積分定数である。

23(2.39)

2.5

Moffatt 下村理論のまとめ

回転ゆで卵が立ち上がることのもっと直観的な説明はないのだろうか? ジェレット定数が本質的であったが，それが回転の十分大きい場合に運動定数になることを数式なしでいうことはまだ誰によってもな
されていない。しかし, もしそれを認めれば，後の説明は極めて簡単である。すべり摩擦によって角速度Ω は時間とともに減少していくであろう。ここで (2.23) 式を用いると，h が増加しなければならない。
つまり，重心が上昇していくことになる。卵も碁石も重心が上昇するように状態を変化させ，最終的にもっとも高い位置に重心がある状態に落ち着くことがこうして理解できる。

24

第 3 章 数値シミュレーションによる回転直立現
象の再現
3.1

プログラムの改良と検証

4 次のルンゲクッタ法を用いて，剛体の直立現象を再現してみた。プログラムは昨年度の村井興平
氏の卒業研究 [6] に使用したものをもとに以下の改良を加えた。
元のプログラムでは水平な床の上に対して軸対称剛体の対称軸 (ζ 軸) が平行になっていて，角速度 ωを鉛直上向きから僅かにずらした方向に与えていたため，回転が始まってすぐに剛体が跳ねてしまっていた。これを防ぐために ζ 軸を水平に平行な状態からずらして，ω を鉛直上向き方向に与えるように改良した。
プログラムを用いて最初に行った計算目的は, プログラム自体の検証である。f (x, ~y ) = d~y /dx を 4 次のルンゲクッタ法で近似するとき

x → x + h で |y~1 | に生じる誤差は O(h5 ) であり
xa → xb で |y~1 | に生じる誤差は, ステップ数が (xb − xa )/h であるから
xb − xa× O(h5 ) → O(h4 )h~ に生じる誤差が h4 に比例していとなるはずである。したがって区間 [xa , xb ] を固定して, その区間で |y|ればよいといえる。また，|∆~y | = ah4 を対数にとるとlog |∆~y | = log a + 4 log hとなることから対数をとったとき傾きが 4 となればよい。
以上のことから今回は h = ∆t = 1.0 × 10−10 のときのエネルギー値 [erg] を真値であると仮定して，ステップ幅 ∆t を横軸に，エネルギー値の誤差を縦軸にとって対数プロットしたときの傾きが 4 次のル
ンゲクッタ法に従っているかどうかを確かめた。

25

100001000d(Energy)[erg]

100
10
1
0.1
0.01
0.001
0.0001
1e-05
1e-06

1e-05

0.0001
∆t

0.001

0.01

図 3.1: プログラムの検証

図は t = 0.2[sec] の時のものである。+は各 ∆t に対する誤差であり，直線の傾きは 4 である。この図から ∆t に対する誤差の傾きは，ほぼ 4 になっているということがわかる。

∆t ≤ 1 × 10−5 ではこの性質がうまく成り立っていないのだが，その原因として可能性が高いのは, プ
ログラム中で ξ, η, ζ 軸間の直交性が僅かずつ失われことだと考えられる。どのようにして直交性を回復させるのがよいかは今後の課題である。

26

3.2

重心位置の上昇の様子

2.5
2.4

z[cm]

2.3
2.2
2.1
2
1.9
1.8

0

0.5

1

1.5

2
2.5
time[sec]

3

3.5

4

図 3.2: 重心の高さの時間変化

次に重心座標の z 成分を調べた。図 3.2 の横軸は時刻 [sec], 縦軸は重心の高さ [cm] の値である。
この時の諸条件は以下のとおりである。

· 剛体の形状楕円体· 半径rξ = 6.0cm, rη = 6.0cm, rζ = 8.0cm· 慣性モーメントIξ = 61.471469, Iη = 61.471469, Iζ = 44.259458· 楕円体の質量m = 33.510322· 回転角速度100π [rad/s]· ステップ幅∆ t=0.0001· ζ 軸の水平方向とのなす角度 0.35π [rad]楕円体の重心は 1.1 秒 ∼1.8 秒の間だけ上昇することが見てとれる。

27

3.3

ジェレット定数の時間変化について

Jellett constant [g· cm2 · rad/sec]

35160
35150
35140
35130
35120
35110
35100
35090
35080
0

0.5

1

1.5

2

2.5

time [sec]図 3.3: ジェレット定数の時間変化

次にジェレット定数の時間変化について調べた。図 3.3 で横軸は時刻 [sec] で, 縦軸はジェレット定数[g · cm2 · rad/sec] である。ここで，軸対称物体が立ち上がる時間は 1.1 秒 ∼1.8 秒であるが，その間にジェレット定数は 0.2% 程度変化しているということがわかった。

28

3.4

ジャイロスコピック近似

Cn, AΩcosθ [g· cm2· rad/sec]CnAΩcosθ
10000

1000

100

0

0.5

1

1.5

2

2.5

time [sec]図 3.4: ジャイロスコピック近似の成立について

ジャイロスコピック近似が成立しているかどうかを調べた。図 3.4 においてここでも軸対称物体が立
ち上がるまで時間は 1.1 秒 ∼1.8 秒であるが，その間に左辺 Cn に対して右辺 AΩ cos θ は最大で約 4%ずれていた。しかし短時間の変動をならして考えるとジャイロスコピック近似はよく成り立っている。
ジャイロスコピック近似の破れが 4%程度あるのに対し, ジェレット定数保存則の破れが 0.2%程度しかないのはこのためであると考えられる。

29

第 4 章 結論
数値計算によるシミュレーションを行い，Moffatt 下村理論の仮定したジャイロスコピック近似が成立しているか，また, ジェレット定数が定数となっているのか，について調べた。その結果，ジャイロス
コピック近似に関しては，立ち上がる過程で左辺の値に対する右辺の値は最大で約 4%程度ずれていた。
これに対してジェレット定数は立ち上がる過程で，約 0.2%程度の変化しかなかった。
このことは，ジャイロスコピック近似は，個々の時刻においては, 大きなずれが生じうるが，短時間の変動をならしてみたときには高い精度で成り立っており，このことによりジェレット定数の変化が約0.2%しかないと解釈することができる。

30

関連図書
[1] H.K.Moffatt and Y.Shimomura: Spinning eggs - a paradox resolved,
Nature 416,385-386(2002).
[2] 下村裕:「立ち上がる回転ゆで卵の解」パリティ,
18,52-56(2003).
[3] 戸田盛和:「回転する卵はなぜ直立する」科学 (岩波書店),
72,932-939(2002).
[4] 小井出昭一郎:「力学」岩波書店,
pp.90-165(1987).
[5] 戸田盛和:「力学」岩波書店,
pp.144-188(1982).
[6] 村井興平:「卵を回すとなぜ立つか」,
福井大学工学部 物理工学科 卒業研究 (2004).

31

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

32

Program List

付録

私はこのプログラム中の初期条件を与える部分を改良した。
A の部分は元々, 関数 kaiten5 内にあったものを分離した。しかしこの条件では回転角速度 ! を鉛直上向きから少しずらした
方向に与えていたために, 物体が床にめり込むようになるため回転が始まってすぐに跳ねてしまうので B のように改良した。
B の部分の説明
回転角速度 ! を鉛直上向きに与えるようにした。このままだと物体は立ち上がらずにその場で回転し続けるだけなので, 物体を傾けた (右図)。このときの実験室系 (x, y, z) から見た重心系 (a, b, c) の単位ベクトルはつぎのようになる。

z

ea

=

(1, 0, 0)eb=(0, cos θ, − sin θ)

ec=(0, sin θ, cos θ)cθ~ と, 重心からの鉛直下向
また, 重心 O の z 座標 rc.z は重心 O から接点へ引いたベクトル tp
きの単位ベクトル ~n との内積で与えられる。
/*
kaiten6h.c*/#include <stdio.h>#include <math.h>/*method = 0 : Euler method, 1 : 4th-order Runge-Kutta method */#define METHOD 1
typedef struct { double x; double y; double z; } Lvec; /* vector in L-frame */typedef struct { double a; double b; double c; } Bvec; /* vector in B-frame */typedef struct { double x; double y; double z; double a; double b; double c; }LBvec; /* vector whose components both in L- and B-frames are necessary */int forces(Lvec *rc, Lvec *vc, Bvec *omg, Lvec *ea, Lvec *eb, Lvec *ec,
Lvec *frc, LBvec *trq);
int normal_point(Bvec *n, Bvec *tp);/*void init1(Lvec *ea, Lvec *eb, Lvec *ec, Bvec *omg0, Lvec *rc, Lvec *vc);*/void init2(Lvec *ea, Lvec *eb, Lvec *ec, Bvec *omg0, Lvec *rc, Lvec *vc);
const double pi = 3.141592653589793;
main(){
kaiten5();}
Bvec R;/* radius in the principal axes */double mass;/* total mass of the rigid body */double volume; /* volume of the rigid body */const double g_gravity = 980.0; /* gravitational acceleration [cm/s^2] */int kaiten5(){
Bvec radius_ratio = { 3.00, 3.00, 4.0};

33

tp

n

y
a

/* relative lengths of principal axes *//* radius for shperical shape [cm] *//* density [g/cm^3] */double R00=2.0 ;double density=1;

Lvec ea /*= { 0.0, 0.0,-1.0}*/ ; /* unit vector for a-axis of B-frame */Lvec eb /*= {-1.0, 0.0, 0.0}*/ ; /* unit vector for b-axis of B-frame */Lvec ec /*= { 0.0, 1.0, 0.0}*/ ; /* unit vector for c-axis of B-frame */double dt;/* time step size [sec] *//*double max_t=30.0;*/
double max_itime=300000;/* the number of time steps */int

mprint=100;

Bvec moi;
Bvec omg0;
Bvec omg;
double omgsize;/* period to print the mechanical state *//* moment of inertia *//* initial angular velocity vector, B-frame *//* angular velocity vector, B-frame */double omg_gosa;
Bvec omgm;/* for Runge Kutta */double R0 ;double t,fct;
double minimum_height ;double eng_tot, eng_rot, eng_tra, eng_gra;
Lvec rc ; /* center of mass coordinate, =(0,0,0) in B-frame */Lvec vc ; /* velocity of center of mass, =(0,0,0) in B-frame */Bvec tp ; /* tangential point */Lvec dtp ; /* displacement from center of mass to tangential point */Lvec vtp ; /* velocity of tangential point */Lvec evtp ; /* unit vector parallel to velocity of tangential point */Lvec frc ; /* total force */LBvec trq ; /* Toruque as for center of mass */double d_eng_tot;
double d_omgsize;
double f1,f2,f3,fn1,fn2,fn3;
Lvec eam, ebm, ecm; /* for Runge Kutta */Lvec rcm,vcm; /* for Runge Kutta */double e1dx,e1dy,e1dz, e2dx,e2dy,e2dz, e3dx,e3dy,e3dz;
double kx,ky,kz, kax,kay,kaz, kbx,kby,kbz, kcx,kcy,kcz;
double krx,kry,krz, kvx,kvy,kvz;
double k1x ,k1y ,k1z , k2x ,k2y ,k2z , k3x ,k3y ,k3z , k4x ,k4y ,k4z ;double k1ax,k1ay,k1az, k2ax,k2ay,k2az, k3ax,k3ay,k3az, k4ax,k4ay,k4az;
double k1bx,k1by,k1bz, k2bx,k2by,k2bz, k3bx,k3by,k3bz, k4bx,k4by,k4bz;
double k1cx,k1cy,k1cz, k2cx,k2cy,k2cz, k3cx,k3cy,k3cz, k4cx,k4cy,k4cz;
double k1rx,k1ry,k1rz, k2rx,k2ry,k2rz, k3rx,k3ry,k3rz, k4rx,k4ry,k4rz;
double k1vx,k1vy,k1vz, k2vx,k2vy,k2vz, k3vx,k3vy,k3vz, k4vx,k4vy,k4vz;
int itime;/* counter of time step, taking 0..max_itime */double Jellett;/* Jellett constant */Bvec mez;
FILE *log_res1; /* output for detailed graphs */FILE *log_res2; /* output for detailed graphs */FILE *log_chk1; /* output for check */FILE *log_res3;
FILE *log_res4;
Lvec angmom; /* angular momentum in L-frame */double angmomsize;
double F1,F2,F3;

34

/* log_res1 = NULL;*/
log_res1=fopen("kaiten.g1","w");/* log_res2 = NULL;*/
log_res2=fopen("kaiten.g2","w");log_chk1 = NULL;/*log_chk1=fopen("kaiten.chk","w");*//* log_res3 = NULL; */ log_res3=fopen("kaiten.g3","w");/* log_res4 = NULL; */log_res4=fopen("kaiten.g4","w");fprintf(stderr,"check : pi = %20.16f\n",pi);
R0=R00/pow(radius_ratio.a*radius_ratio.b*radius_ratio.c,1.0/3.0);
R.a=R0*radius_ratio.a; R.b=R0*radius_ratio.b; R.c=R0*radius_ratio.c;
if(R.a < R.b) minimum_height=R.a; else minimum_height = R.b;if(R.c < minimum_height) minimum_height = R.c;fprintf(stderr,"R=(%f %f %f) min=%f\n",R.a,R.b,R.c,minimum_height);
volume=4*pi*R.a*R.b*R.c/3;
mass=density*volume;/* volume of the rigid body [cm^3] */fprintf(stderr,"volume=%f density=%f mass=%f\n",volume,density,mass);moi.a=mass*(R.b*R.b+R.c*R.c)/5;/* moment inertia [g cm^2] */moi.b=mass*(R.c*R.c+R.a*R.a)/5;
moi.c=mass*(R.a*R.a+R.b*R.b)/5;
fprintf(stderr,"MOI=(%f %f %f)\n",moi.a,moi.b,moi.c);f1=(moi.b-moi.c)/moi.a; f2=(moi.c-moi.a)/moi.b; f3=(moi.a-moi.b)/moi.c;
fn1=1/moi.a; fn2=1/moi.b; fn3=1/moi.c;
fprintf(stderr,"input dt=");scanf("%lf",&dt);
fprintf(stderr,"check dt=%e\n",dt);/*max_itime=max_t/dt;
mprint=max_itime/3000.0; if(dt<=0.00001&&mprint%10!=0) mprint++;
fprintf(stderr,"max_itime=%f mprint=%d\n",max_itime,mprint);*//*init1(&ea, &eb, &ec, &omg0, &rc ,&vc);*/init2(&ea, &eb, &ec, &omg0, &rc ,&vc);
fprintf(stderr,"rc.z=%f\n",rc.z);

fprintf(stderr,"%s\n ec - (ea x eb) =(%e %e %e)\n","check of right-handedness of vectors {ea,eb,ec}:",ec.x - (ea.y*eb.z-ea.z*eb.y),ec.y - (ea.z*eb.x-ea.x*eb.z),ec.z - (ea.x*eb.y-ea.y*eb.x));
omg.a=omg0.a;

omg.b=omg0.b;

omg.c=omg0.c;

for(itime=0;;itime++){ t=itime*dt;

angmom.x=moi.a*omg.a*ea.x + moi.b*omg.b*eb.x + moi.c*omg.c*ec.x ;angmom.y=moi.a*omg.a*ea.y + moi.b*omg.b*eb.y + moi.c*omg.c*ec.y ;angmom.z=moi.a*omg.a*ea.z + moi.b*omg.b*eb.z + moi.c*omg.c*ec.z ;angmomsize=sqrt(angmom.x*angmom.x + angmom.y*angmom.y + angmom.z*angmom.z);
omgsize = sqrt(omg.a*omg.a + omg.b*omg.b + omg.c*omg.c);eng_rot = (moi.a*omg.a*omg.a + moi.b*omg.b*omg.b + moi.c*omg.c*omg.c)*0.5 ;eng_tra = (0.5*mass)*(vc.x*vc.x + vc.y*vc.y + vc.z*vc.z) ;eng_gra = (mass*g_gravity)*(rc.z-minimum_height);

35eng_tot = eng_rot + eng_tra + eng_gra;if(itime==0) fprintf(stderr,"eng_tot=%f\n",eng_tot);if(log_res1 != NULL && itime % mprint == 0)
fprintf(log_res1,"%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n",t,omg.a,omg.b,omg.c,ea.x,ea.y,ea.z,eb.x,eb.y,eb.z,ec.x,ec.y,ec.z,rc.x,rc.y,rc.z,vc.x,vc.y,vc.z);d_eng_tot=eng_tot-2841751.291223551612347364;/*if(t==0.2) {fprintf(stderr,"eng=%30.18f,omega=%f\n",eng_tot,omgsize);fprintf(log_res3,"%30.18f %20.18f\n",dt,fabs(d_eng_tot));
fprintf(log_res4,"%30.18f %30.18f\n",dt,fabs(d_omgsize));}*//*jellett constant*/
mez.a=-ea.z; mez.b=-eb.z; mez.c=-ec.z;/*fprintf(stderr,"mez(%f %f %f)\n",mez.a,mez.b,mez.c);*/normal_point(&mez,&tp);
Jellett=-(moi.a*omg.a*tp.a + moi.b*omg.b*tp.b + moi.c*omg.c*tp.c);if(itime%mprint==0) fprintf(log_res3,"%f %f\n",t,Jellett);/*Gyroscopic balance*/
F1= moi.c * omg.c;
F2 = moi.a * ec.z*(ec.x*(omg.b*ea.y-omg.a*eb.y)-ec.y*(omg.b*ea.x-omg.a*eb.x))/(ec.x*ec.x+ec.y*ec.y );F3=F1-F2;
if (itime%mprint==0) fprintf(log_res4,"%f %f %f\n",t,F1,F2);/*if (itime%mprint==0) fprintf(log_res4,"%f %f\n",t,fabs(F3));*//*if(itime%mprint==0) fprintf(log_res4,"%f %f\n",t,fabs((2*F3)/(F1+F2)));*//*Omega Size*/if(t==0) fprintf(stderr,"init OMG=%f\n",omgsize);/*if (t>2.4228) {fprintf(stderr,"omg=%f\n",omgsize); break;}*/if(log_chk1 != NULL && itime % mprint == 0){/*
fprintf(log_chk1,"%f %e %e %e %e %e %e\n",t,ea.x*ea.x + ea.y*ea.y + ea.z*ea.z -1.0,
eb.x*eb.x + eb.y*eb.y + eb.z*eb.z -1.0,
ec.x*ec.x + ec.y*ec.y + ec.z*ec.z -1.0,
ea.x*eb.x + ea.y*eb.y + ea.z*eb.z,
eb.x*ec.x + eb.y*ec.y + eb.z*ec.z,
ec.x*ea.x + ec.y*ea.y + ec.z*ea.z);*/
fprintf(log_chk1,"%f %e\n",t,fabs(ea.x*eb.x + ea.y*eb.y + ea.z*eb.z)+fabs(eb.x*ec.x + eb.y*ec.y + eb.z*ec.z)+fabs(ec.x*ea.x + ec.y*ea.y + ec.z*ea.z));}
if(itime >= max_itime) break;forces(&rc, &vc, &omg, &ea, &eb, &ec, &frc, &trq);36

if(log_res2 != NULL && itime % mprint == 0)
fprintf(log_res2,"%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n",t,angmom.x,angmom.y,angmom.z,angmomsize,omgsize,eng_rot,eng_tra,eng_gra,eng_tot,frc.x,frc.y,frc.z, trq.x,trq.y,trq.z, trq.a,trq.b,trq.c);k1x=dt*(f1*omg.b*omg.c+fn1*trq.a);k1y=dt*(f2*omg.c*omg.a+fn2*trq.b);k1z=dt*(f3*omg.a*omg.b+fn3*trq.c);k1ax=dt*(omg.c*eb.x-omg.b*ec.x);k1ay=dt*(omg.c*eb.y-omg.b*ec.y);k1az=dt*(omg.c*eb.z-omg.b*ec.z);k1bx=dt*(omg.a*ec.x-omg.c*ea.x);
k1by=dt*(omg.a*ec.y-omg.c*ea.y);k1bz=dt*(omg.a*ec.z-omg.c*ea.z);k1cx=dt*(omg.b*ea.x-omg.a*eb.x);k1cy=dt*(omg.b*ea.y-omg.a*eb.y);k1cz=dt*(omg.b*ea.z-omg.a*eb.z);
k1rx=dt*vc.x;
k1ry=dt*vc.y;
k1rz=dt*vc.z;
k1vx=(dt/mass)*frc.x;
k1vy=(dt/mass)*frc.y;
k1vz=(dt/mass)*frc.z;
omgm.a=omg.a+k1x/2; omgm.b=omg.b+k1y/2; omgm.c=omg.c+k1z/2;
eam.x=ea.x+k1ax/2; eam.y=ea.y+k1ay/2; eam.z=ea.z+k1az/2;
ebm.x=eb.x+k1bx/2; ebm.y=eb.y+k1by/2; ebm.z=eb.z+k1bz/2;
ecm.x=ec.x+k1cx/2; ecm.y=ec.y+k1cy/2; ecm.z=ec.z+k1cz/2;
rcm.x=rc.x+k1rx/2; rcm.y=rc.y+k1ry/2; rcm.z=rc.z+k1rz/2;
vcm.x=vc.x+k1vx/2; vcm.y=vc.y+k1vy/2; vcm.z=vc.z+k1vz/2;/*

printf("%f %f %f %f\n",vc.z,frc.z,k1vz,vcm.z);*/forces(&rcm, &vcm, &omgm, &eam, &ebm, &ecm, &frc, &trq);
k2x=dt*(f1*omgm.b*omgm.c+fn1*trq.a);k2y=dt*(f2*omgm.c*omgm.a+fn2*trq.b);k2z=dt*(f3*omgm.a*omgm.b+fn3*trq.c);k2ax=dt*(omgm.c*ebm.x-omgm.b*ecm.x);k2ay=dt*(omgm.c*ebm.y-omgm.b*ecm.y);k2az=dt*(omgm.c*ebm.z-omgm.b*ecm.z);k2bx=dt*(omgm.a*ecm.x-omgm.c*eam.x);k2by=dt*(omgm.a*ecm.y-omgm.c*eam.y);k2bz=dt*(omgm.a*ecm.z-omgm.c*eam.z);k2cx=dt*(omgm.b*eam.x-omgm.a*ebm.x);k2cy=dt*(omgm.b*eam.y-omgm.a*ebm.y);k2cz=dt*(omgm.b*eam.z-omgm.a*ebm.z);37

k2rx=dt*vcm.x;k2ry=dt*vcm.y;k2rz=dt*vcm.z;k2vx=(dt/mass)*frc.x;
k2vy=(dt/mass)*frc.y;
k2vz=(dt/mass)*frc.z;
omgm.a=omg.a+k2x/2; omgm.b=omg.b+k2y/2; omgm.c=omg.c+k2z/2;
eam.x=ea.x+k2ax/2; eam.y=ea.y+k2ay/2; eam.z=ea.z+k2az/2;
ebm.x=eb.x+k2bx/2; ebm.y=eb.y+k2by/2; ebm.z=eb.z+k2bz/2;
ecm.x=ec.x+k2cx/2; ecm.y=ec.y+k2cy/2; ecm.z=ec.z+k2cz/2;
rcm.x=rc.x+k2rx/2; rcm.y=rc.y+k2ry/2; rcm.z=rc.z+k2rz/2;
vcm.x=vc.x+k2vx/2; vcm.y=vc.y+k2vy/2; vcm.z=vc.z+k2vz/2;/*

printf("%f %f %f %f %f %f\n",rcm.z,k2rz,vc.z,frc.z,k2vz,vcm.z);*/forces(&rcm, &vcm, &omgm, &eam, &ebm, &ecm, &frc, &trq);
k3x=dt*(f1*omgm.b*omgm.c+fn1*trq.a);k3y=dt*(f2*omgm.c*omgm.a+fn2*trq.b);k3z=dt*(f3*omgm.a*omgm.b+fn3*trq.c);k3ax=dt*(omgm.c*ebm.x-omgm.b*ecm.x);k3ay=dt*(omgm.c*ebm.y-omgm.b*ecm.y);k3az=dt*(omgm.c*ebm.z-omgm.b*ecm.z);k3bx=dt*(omgm.a*ecm.x-omgm.c*eam.x);k3by=dt*(omgm.a*ecm.y-omgm.c*eam.y);k3bz=dt*(omgm.a*ecm.z-omgm.c*eam.z);k3cx=dt*(omgm.b*eam.x-omgm.a*ebm.x);k3cy=dt*(omgm.b*eam.y-omgm.a*ebm.y);k3cz=dt*(omgm.b*eam.z-omgm.a*ebm.z);k3rx=dt*vcm.x;k3ry=dt*vcm.y;k3rz=dt*vcm.z;
k3vx=(dt/mass)*frc.x;
k3vy=(dt/mass)*frc.y;
k3vz=(dt/mass)*frc.z;
omgm.a=omg.a+k3x; omgm.b=omg.b+k3y; omgm.c=omg.c+k3z;
eam.x=ea.x+k3ax; eam.y=ea.y+k3ay; eam.z=ea.z+k3az;
ebm.x=eb.x+k3bx; ebm.y=eb.y+k3by; ebm.z=eb.z+k3bz;
ecm.x=ec.x+k3cx; ecm.y=ec.y+k3cy; ecm.z=ec.z+k3cz;
rcm.x=rc.x+k3rx; rcm.y=rc.y+k3ry; rcm.z=rc.z+k3rz;
vcm.x=vc.x+k3vx; vcm.y=vc.y+k3vy; vcm.z=vc.z+k3vz;/*

printf("%f %f %f %f %f %f\n",rcm.z,k3rz,vc.z,frc.z,k3vz,vcm.z);*/forces(&rcm, &vcm, &omgm, &eam, &ebm, &ecm, &frc, &trq);
k4x=dt*(f1*omgm.b*omgm.c+fn1*trq.a);k4y=dt*(f2*omgm.c*omgm.a+fn2*trq.b);38

k4z=dt*(f3*omgm.a*omgm.b+fn3*trq.c);k4ax=dt*(omgm.c*ebm.x-omgm.b*ecm.x);k4ay=dt*(omgm.c*ebm.y-omgm.b*ecm.y);k4az=dt*(omgm.c*ebm.z-omgm.b*ecm.z);k4bx=dt*(omgm.a*ecm.x-omgm.c*eam.x);k4by=dt*(omgm.a*ecm.y-omgm.c*eam.y);k4bz=dt*(omgm.a*ecm.z-omgm.c*eam.z);k4cx=dt*(omgm.b*eam.x-omgm.a*ebm.x);k4cy=dt*(omgm.b*eam.y-omgm.a*ebm.y);k4cz=dt*(omgm.b*eam.z-omgm.a*ebm.z);k4rx=dt*vcm.x;k4ry=dt*vcm.y;k4rz=dt*vcm.z;k4vx=(dt/mass)*frc.x;
k4vy=(dt/mass)*frc.y;
k4vz=(dt/mass)*frc.z;
kx=(k1x+2*(k2x+k3x)+k4x)*(1.0/6.0);
ky=(k1y+2*(k2y+k3y)+k4y)*(1.0/6.0);
kz=(k1z+2*(k2z+k3z)+k4z)*(1.0/6.0);
kax=(k1ax+2*(k2ax+k3ax)+k4ax)*(1.0/6.0);
kay=(k1ay+2*(k2ay+k3ay)+k4ay)*(1.0/6.0);
kaz=(k1az+2*(k2az+k3az)+k4az)*(1.0/6.0);
kbx=(k1bx+2*(k2bx+k3bx)+k4bx)*(1.0/6.0);
kby=(k1by+2*(k2by+k3by)+k4by)*(1.0/6.0);kbz=(k1bz+2*(k2bz+k3bz)+k4bz)*(1.0/6.0);
kcx=(k1cx+2*(k2cx+k3cx)+k4cx)*(1.0/6.0);
kcy=(k1cy+2*(k2cy+k3cy)+k4cy)*(1.0/6.0);
kcz=(k1cz+2*(k2cz+k3cz)+k4cz)*(1.0/6.0);
krx=(k1rx+2*(k2rx+k3rx)+k4rx)*(1.0/6.0);
kry=(k1ry+2*(k2ry+k3ry)+k4ry)*(1.0/6.0);
krz=(k1rz+2*(k2rz+k3rz)+k4rz)*(1.0/6.0);
kvx=(k1vx+2*(k2vx+k3vx)+k4vx)*(1.0/6.0);
kvy=(k1vy+2*(k2vy+k3vy)+k4vy)*(1.0/6.0);
kvz=(k1vz+2*(k2vz+k3vz)+k4vz)*(1.0/6.0);/*

printf("%f %f %f %f\n",krz,frc.z,k4vz,kvz);*/omg.a=omg.a+kx;
omg.b=omg.b+ky;
omg.c=omg.c+kz;
ea.x=ea.x+kax;
ea.y=ea.y+kay;
ea.z=ea.z+kaz;
eb.x=eb.x+kbx;

39

eb.y=eb.y+kby;
eb.z=eb.z+kbz;
ec.x=ec.x+kcx;
ec.y=ec.y+kcy;
ec.z=ec.z+kcz;
rc.x=rc.x+krx;
rc.y=rc.y+kry;
rc.z=rc.z+krz;
vc.x=vc.x+kvx;
vc.y=vc.y+kvy;
vc.z=vc.z+kvz;
fct=1/sqrt(ea.x*ea.x+ea.y*ea.y+ea.z*ea.z);ea.x=fct*ea.x;ea.y=fct*ea.y;ea.z=fct*ea.z;
fct=1/sqrt(eb.x*eb.x+eb.y*eb.y+eb.z*eb.z);eb.x=fct*eb.x;eb.y=fct*eb.y;eb.z=fct*eb.z;
fct=1/sqrt(ec.x*ec.x+ec.y*ec.y+ec.z*ec.z);ec.x=fct*ec.x;ec.y=fct*ec.y;ec.z=fct*ec.z;}
if (log_res1 != NULL) fclose(log_res1);
if (log_res2 != NULL) fclose(log_res2);
if (log_chk1 != NULL) fclose(log_chk1);}
int forces(Lvec *rc, Lvec *vc, Bvec *omg, Lvec *ea, Lvec *eb, Lvec *ec,
Lvec *frc, LBvec *trq) {/*uses global variables : mass, g_gravity,*/
const double mu_friction = 0.5 ;const double tspeed0 = 1.0e-3 ;const double d_fr = 0.05; /* [cm] */Lvec omgL ; /* angular velocity vector in L-frame */Bvec tp ; /* tangent point in B-frame */Bvec mez ; /* unit vector parallel to the gravity */Lvec dtp ; /* displacement from center of rotor to tangent point */Lvec vtp ; /* velocity of tangent point */Lvec frctp ; /* force operating at the tangent point */double htp ; /* height of tangent point */double fr ; /* normal reaction force */double ff ; /* tangential friction force */double tspeed ; /* tangential speed */double fct,t1 ;double mufrict ; /* friction coefficient mu */omgL.x = omg->a*ea->x + omg->b*eb->x + omg->c*ec->x ;omgL.y = omg->a*ea->y + omg->b*eb->y + omg->c*ec->y ;omgL.z = omg->a*ea->z + omg->b*eb->z + omg->c*ec->z ;/*fprintf(stderr,"omgL.z=%f\n",omgL.z);*/mez.a = - ea->z ;mez.b = - eb->z ;mez.c = - ec->z ;normal_point(&mez, &tp);

40

dtp.x = tp.a*ea->x + tp.b*eb->x + tp.c*ec->x ;dtp.y = tp.a*ea->y + tp.b*eb->y + tp.c*ec->y ;dtp.z = tp.a*ea->z + tp.b*eb->z + tp.c*ec->z ;
vtp.x
vtp.y
vtp.z

= vc->x + omgL.y * dtp.z - omgL.z * dtp.y ;= vc->y + omgL.z * dtp.x - omgL.x * dtp.z ;= vc->z + omgL.x * dtp.y - omgL.y * dtp.x ;/*printf("%f %f %f %f %f %f\n",vc->x,vc->y,vc->z,vtp.x,vtp.y,vtp.z);*/
htp = rc->z + dtp.z ;t1=htp*(1.0/d_fr); if(t1 < -5.0) t1=-5.0;
fr = (mass*g_gravity)*exp(-t1*(1.0+0.1*t1))* (1- 0.1*tanh(vc->z * (1.0/15.0)));/*t1=htp*(1.0/d_fr); if(t1 < -4.6) t1=-4.6;
fr = (mass*g_gravity)*exp(-t1);*//*
printf("%f %f %f %f %f\n",rc->z, dtp.z, htp, t1, fr); if(1) exit(1);*/tspeed = sqrt(vtp.x * vtp.x + vtp.y * vtp.y);mufrict = mu_friction * tanh(tspeed * (1.0/tspeed0));
ff = fr * mufrict ;if(tspeed > 1.0e-32) fct = - ff / tspeed ; else fct = 0.0;
frctp.x = fct* vtp.x ;frctp.y = fct* vtp.y ;frctp.z = fr ;frc->x = frctp.x ;frc->y = frctp.y ;frc->z = frctp.z - mass * g_gravity ;trq->x = dtp.y * frctp.z - dtp.z * frctp.y ;trq->y = dtp.z * frctp.x - dtp.x * frctp.z ;trq->z = dtp.x * frctp.y - dtp.y * frctp.x ;
trq->a = trq->x * ea->x + trq->y * ea->y + trq->z * ea->z ;trq->b = trq->x * eb->x + trq->y * eb->y + trq->z * eb->z ;trq->c = trq->x * ec->x + trq->y * ec->y + trq->z * ec->z ;/*printf("%f %f %f %f %f %f F\n",rc->z,dtp.z,htp,vc->z,fr,frc->z);*//*printf("%f %f %f %f %f %f %f\n",mufrict,fr,ff,frctp.x,frctp.y,frctp.z,frc->z);if(1) exit(1);*/return 0;}
int normal_point(Bvec *n, Bvec *tp){/*41

For ellipsoid.
uses global variables : Bvec R*/
double fct;tp->a=R.a*R.a*n->a; tp->b=R.b*R.b*n->b; tp->c=R.c*R.c*n->c;
fct = 1.0/sqrt(tp->a*n->a + tp->b*n->b + tp->c*n->c) ;tp->a = tp->a*fct;tp->b = tp->b*fct;
tp->c = tp->c*fct;/*
printf("%f %f %f\n",R.a,R.b,R.c);printf("%f %f %f\n",n->a,n->b,n->c);printf("%f %f %f\n",tp->a,tp->b,tp->c);if(1) exit(1);*/return 0;}/*
If call the init1 , rc and vc are necessary toestablish the initial value.
*//*
void init1(Lvec *ea, Lvec *eb, Lvec *ec, Bvec
*omg0, Lvec *rc, Lvec *vc){
double b;
Bvec angular_velocity_ratio = {-1.0, 0.0, -0.01};
double angular_velocity = 100*pi;
ea->x=0.0; ea->y=0.0;
eb->x=-1.0; eb->y=0.0;
ec->x=0.0; ec->y=1.0;

ea->z=-1.0;
eb->z=0.0;
ec->z=0.0;

b=1/sqrt(angular_velocity_ratio.a*angular_velocity_ratio.a+angular_velocity_ratio.b*angular_velocity_ratio.b+angular_velocity_ratio.c*angular_velocity_ratio.c);
omg0->a=angular_velocity*angular_velocity_ratio.a*b;
omg0->b=angular_velocity*angular_velocity_ratio.b*b;
omg0->c=angular_velocity*angular_velocity_ratio.c*b;
rc->x=0;
vc->x=0;

rc->y=0;
vc->y=0;

rc->z=R.a;vc->z=0;

}
*/

42

A

void init2(Lvec *ea, Lvec *eb, Lvec *ec, Bvec *omg0,
Lvec *rc, Lvec *vc){
double b;
double theta=(85.0/90.0)*0.5*pi;
Bvec angular_velocity_ratio = {0.0, -sin(theta), cos(theta)};/*relative size of body-frame components
of initial angular velocity */Bvec n ={0.0, sin(theta), -cos(theta)};
Bvec tp;
double angular_velocity = 100*pi;/*size fo initial ang. vel.[radian/sec]*/
ea->x=1.0;

ea->y=0.0;

ea->z=0.0;

eb->x=0.0;

eb->y=cos(theta);

eb->z=-sin(theta);

ec->x=0.0;

ec->y=sin(theta);

ec->z=cos(theta);

b=1/sqrt(angular_velocity_ratio.a*angular_velocity_ratio.a+angular_velocity_ratio.b*angular_velocity_ratio.b+angular_velocity_ratio.c*angular_velocity_ratio.c);
omg0->a=angular_velocity*angular_velocity_ratio.a*b;
omg0->b=angular_velocity*angular_velocity_ratio.b*b;
omg0->c=angular_velocity*angular_velocity_ratio.c*b;
normal_point(&n, &tp);
rc->x=0; rc->y=0; rc->z=n.a*tp.a + n.b*tp.b + n.c*tp.c;
vc->x=0; vc->y=0; vc->z=0;
fprintf(stderr,"%f,%f,%f,%f,%f,%f\n",tp.a,tp.b,tp.c,omg0->a,omg0->b,omg0->c);}43

B