保存形での剛体の回転直立現象

2006 年 2 月 17 日
福井大学 工学部 物理工学科
１４年度入学 ０２３８００９６
加藤友直

目次
第 1 章 回転直立現象

1.1
1.2
1.3

1.4

回転直立現象とは . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Moffatt-下村理論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 剛体回転のオイラー方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.2 座標系と接触点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.3 ジャイロスコピック近似 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.4 ジェレット定数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Moffatt-下村理論の問題点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
研究の歴史

第 2 章 剛体の運動の数値シミュレーション

2.1
2.2

剛体の運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.3
2.4
2.5

2.2.1 剛体の運動の力学的記述 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
剛体のエネルギー . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
初期条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
摩擦力の効果を調べる . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

剛体の形状. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 3 章 摩擦力の変化による回転直立の変化

3.1
3.2

3.3
3.4
3.5

動摩擦係数を変えた場合の影響 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
摩擦力の方向を人為的にずらした場合の影響 . . . . . . . . . . . . . . . . . . . . . . . .

3.2.1 |η| < 90◦ の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 |η| = 90◦ の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
|η| > 90◦ の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
η を変化させたときの直立時間の変化 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
時間反転解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 4 章 結論

2
2
3
4
4
6
9
10
11
12
12
14
14
16
17
19
20
21
21
24
24
26
27
29
30
32

参考文献

33

謝辞

34

付録

35

1

第 1 章 回転直立現象
1.1

回転直立現象とはゆで卵や碁石などの剛体を勢いよく回転させると直立する。この現象は１００年以上前から知られていたが、長い間理論的には解明されていなかった。２００２年に Moffatt と下村によって卵のような軸対称物体の回転に伴う立ち上がりの現象の近似解が示された。しかし、これは近似解であるため様々な問題点がある。そこで本卒業研究では広い状況下での回転直立現象の普遍性を調べるために、摩擦力が唯一絶対必要な条件かどうかを調べることにした。

図 1.1: 回転直立現象

回転直立現象をゆで卵を例にとって説明する。図 1.1 のように対称軸であるζを横に倒し勢いよくゆで卵を回転させると重心の位置が上昇し直立するという現象である。このことはゆで卵のようなプロレー
ト形 (対称軸方向に長く伸びた回転体) についてだけ起こるものではなく、碁石のようなオブレート形(対称軸方向が短い回転体) についても同様のことが言え、勢いよく回転させると直立する。ただし碁石の場合は初め対称軸を鉛直にした状態で回転させる。
なお回転直立という言葉は、文献 [1] のタイトルから我々がつくった造語であるので WEB サーチしてもヒットしないことを断っておく。

2

1.2

研究の歴史

本研究以前にも様々な研究が行われているが、そのうち最も重要なものとして以下の２つが挙げられる。
まず、１９５２年に Braams, Hugenholtz によって逆立ちごまの厳密解が示された [2, 3]。

図 1.2: 逆立ちごまの回転

逆立ちごまとは図 1.2 のような球の部分上の部分をとって棒を付けたものである。これは幾何学的には偏心球 (重心のずれた球) である。
次に２００２年、５０年ぶりに Moffatt, 下村によって任意の剛体の近似解が示された [4]。これについては次節で詳しく説明する。
そして、昨年、一昨年と卒業研究で数値シュミレーション法の開発、改良、計算結果の考察が村井と大森によって行われ [5, 6]、本年度木澤が様々な形状で普遍的に直立することを示した。

3

1.3

Moffatt-下村理論

この節では文献 [4, 8] をもとに Moffatt-下村理論について説明する。

1.3.1

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

次ページからの説明で使用するオイラー方程式を以下で導出する。
~~を剛体上のベクトル B(t)の微小時間 dt における変化 dB~ の変化(1) 剛体系 O-ξηζ からみたときの B~ の変化(2) 実験室系 O-xyz からみたときの Bの 2 つの場合に分けて考えてみる。
~(1) の場合の ξ, η, ζ 成分は dBξ , dBη , dBζ である。(2) の場合は (~ω × B)dtとなる。これらのベクトル和~ を dt で割るとである dB
~
dB
dt
~
dB
dt
~
dB
dt

!
=

dBξ
~ ξ ,
+ (~
ω × B)
dt

=

dBη
~ η ,
+ (~
ω × B)
dt

=

dBζ
~ ζ ,
+ (~
ω × B)
dt

ξ

!

η

!

ζ

(1.1)

~は
となる。ξ, η, ζ 軸を慣性主軸に一致させたときの O に関する全角運動量 L
Lξ = Iξ0 ωξ ,Lη = Iη0 ωη ,Lζ = Iζ0 ωζ

(1.2)

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

dLη
dωξ
dωη
dLξ
= Iξ0
,
= Iη0
,
dt
dt
dt
dt
~ として L~ を用い、(1.3) 式を使うととなる。(1.1) 式で B
~
dL
dt

!

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

dωξ
+ (ωη Lζ − ωζ Lη )
dt


Lζ
Lη
0 dωξ
= Iξ
+
−
ωη ωζ
dt
ωζ
ωη
dωξ
− (Iη0 − Iζ0 )ωη ωζ
= Iξ0
dt

(1.3)

= Iξ0

ξ

(1.4)

となる。同様に η, ζ 成分は

~
dL
dt
~
dL
dt

!
= Iη0

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

(1.5)

= Iζ0

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

(1.6)

η

!

ζ

4

となり角運動量に対する運動方程式

~
dL
dt

~
= N
X
=
(~ri × F~i )

(1.7)

i

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

(1.8)

となる。これを剛体の回転に対するオイラー方程式という。
本論文では以降からは ξ, η, ζ のかわりに X1 , X2 , X3 を用い、それぞれの座標軸を１軸, ２軸, ３軸と呼ぶことにする。

5

1.3.2座標系と接触点

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

(1.9)

~r
・ 垂直抗力 F(1.10)~f
・ Y 方向への摩擦力 F(1.11)

6である。
物体は X, Y, Z 系でみて、３軸のまわりに角速度 ω30 で回転しているとする。また、２軸 (= Y 軸) のま
わりには角速度 θ̇ で緩慢回転している。実験室系で見た剛体の角速度ベクトルは実験室系でみた X, Y, Z
系の角速度ベクトルと X, Y, Z 系で見た剛体の角速度ベクトルの和である。したがって物体の角速度ベ
クトル ω~ の各成分はω1= −Ω sin θ ,ω2= θ̇ ,

ω3= n,(1.12)となる。ただし n は ω の 3 軸成分で n = Ω cos θ + ω30 である。また水平面内回転系での各成分は

(~
ω )X

= ω1 (~e1 )X + ω2 (~e2 )X + ω3 (~e3 )X(1.13)= −ω1 cos θ + ω3 sin θ= −Ω sin θ cos θ + n sin θ
=
(~
ω )Y

(n − Ω cos θ) sin θ ,= ω1 (~e1 )Y + ω2 (~e2 )Y + ω3 (~e3 )Y(1.14)= ω2 ,= θ̇ ,(~ω )Z= ω1 (~e1 )Z + ω2 (~e2 )Z + ω3 (~e3 )Z(1.15)= ω3 cos θ + ω1 sin θ= n cos θ + Ω sin2 θ .
ただし e~1 , e~2 , e~3 はそれぞれの単位ベクトルである。

1,2,3 軸の主慣性モーメントをそれぞれ I1 , I2 , I3 とする。１, ２軸の主観性モーメントは等しいとする~ の各成分は (1.12) 式よりと (I1 = I2 )、角運動量 L

L1= −I1 Ω sin θ ,L2= I2 θ̇ (= I1 θ̇) ,L3= I3 n ,(1.16)となる。また

LX= L1 cos θ + L3 sin θ ,

LY

= L2

LZ

= −L1 sin θ + L3 cos θ ,(1.17)~はである。(1.16) 式を (1.17) 式に代入して水平面内回転系にうつると L~ = ((I3 n − I1 Ω cos θ) sin θ , I1 θ̇ , I1 Ω sin2 θ + I3 n cos θ)L7(1.18)

となる。
~ の時間に対する変化は，オイラーの方程式
この L

~
dL
~ =X
~ p × F~r + F~f+ (ω × L)
dt(1.19)~f は点 P にはたらく摩擦力、F~r は垂直抗力である。
で与えられる。ここで F~
位置ベクトル XP = (XP , 0, ZP ) の成分は

XP =

dh
,
dθ

ZP = −h(θ)(1.20)と表される [4]。
~r + F~f = (0, Ff , Fr ) であるから
また，F

XP × (F~r + F~f ) = (−ZP Ff , −Fr XP , Ff XP )(1.21)となり (1.19) 式を用いるとオイラーの方程式の各成分は

d[(I3 n − I1 Ω cos θ) sin θ] − I1 Ωθ̇ = −ZP Ff
dt
I1 θ̈ + Ω(I3 n − I1 Ω cos θ) sin θ = −Fr XPd
I1 Ω̇ + [(I3 n − I1 Ω cos θ) cos θ] = Ff XP
dtと表される。

8(1.22)(1.23)(1.24)

1.3.3

ジャイロスコピック近似

軸対称物体の立ち上がり現象を考えるとき，回転が遅いと立ち上がらないので回転が速い場合について考える。オイラー方程式の Y 成分は

I1 θ̈ − I1 Ω2 cos θ sin θ + I3 nΩ sin θ = −Fr XP(1.25)である。Ω2 が大きいとすると，Ω が含まれている左辺の第 2 項と第 3 項は右辺の −Fr XP よりも十分に大きいとすることができる。次に時間スケールに着目すると軸対称物体が立ち上がる時間は回転で決まる短いものではなく摩擦で決まる長いものであるので |θ̈|  Ω とすることができる。よって (1.25) 式
は近似的に(I3 n − I1 Ω cos θ)Ω sin θ = 0(1.26)とでき，立ち上がるまでの過程すなわち，Ω sin θ 6= 0 では

I3 n = I1 Ω cos θ(1.27)となる。これをジャイロスコピック近似とよぶ。
また，オイラー方程式の X 成分，(1.22) 式はI1 Ωθ̇ = Ff ZP(1.28)

I1 Ω̇ = Ff XP(1.29)となる。そして Z 成分，(1.24) 式はとなり (1.28) 式，(1.29) 式および (1.20) 式から

Ω̇
XP
ḣ
=
θ̇ = −Ω
ZP
h(1.30)

が導かれる。この (1.30) 式を積分すると

Z

Ω̇
dt
Ω
log Ω

Z

ḣ
dt ,
h
− log h + C .

= −(1.31)=(1.32)

よって log hΩ = C となり Ωh= 一定 となる。
この逆数に I1 をかけたものが次節で述べるジェレット定数である。

9

1.3.4

ジェレット定数

まずジェレット定数を説明するにあたり先ほど説明した卵や逆立ちごまが、どのような条件で重心が上昇するか説明する。重心が上昇するには・剛体と接触面との摩擦が不可欠である。
・接触点がすべって力学的エネルギーが減少しなければならない。
ということが分かっている。このようにエネルギーが保存しない系であるにもかかわらず、ひとつの運~ と接触点の位置ベクトル XP との内積として動定数が存在する。その運動定数は、角運動量ベクトル L~ ·X~p
J = −L

(1.33)

で定義される量であり、ジェレット定数と呼ばれる。
任意の軸対称物体に対するジェレット定数の時間微分を計算するとJ˙~˙ · X~p − L
~ ·X
~˙ p
= −L
(I3 n − I1 Ω cos θ)XP2=

d sin θdt Xp

(1.34)となる。このとき (1.18) 式 (1.20) 式を使用した。また、ここでジャイロスコピック近似が成立すると(1.34) 式は J˙ = 0 となり J は運動定数となる。

~はまた (1.18) 式で与えられる角運動量 L~ = (0 , I1 θ̇ , I1 Ω)
L(1.35)となる。したがって (1.35) 式から

J= −(−h(θ) · I1 Ω)= I1 Ωh(1.36)と表せることができる。このとき I1 , J は定数であるので摩擦力により Ω が減少すると h が増加する。
そのため、重心の位置が上昇し剛体が立ち上がる。
また (1.20) 式，(1.36) 式，(1.28) 式より任意の軸対称物体の θ に対するジャイロスコピック近似を用いた方程式はJ θ̇ = −Ff h2 (θ)(1.37)という簡単な 1 階の微分方程式となり、このとき h2 (θ) = a2 − (a2 − b2 ) sin2 θ (a.b：短軸, 長軸の半径)とすると一般解はaarctan
b



b
tan θ
a


=

である。

10

Ff a2
t+C
J

(1.38)

1.4

Moffatt-下村理論の問題点

Moffatt-下村理論の問題点として１つ目に、仮定するジャイロスコピックバランスは時間平均してはじめて成立することがあげられる [6]。ジャイロスコピックバランスは個々の時間だと大きなズレが生じているのである。
２つ目に自発的にジャンプするということがあげられる [4, 5, 8, 9]。図 1.4 のように剛体は回転させ立ち上がってくると、ジャンプしてしまう。そのためジャンプをして着地した後の解析的研究が行いづらい。

図 1.4: 剛体のジャンプ

そこで我々の研究室では Moffatt-下村理論が仮定したような滑らかな直立過程に対象を限ることなく、より広い状況下での回転直立現象の普遍性を調べることを目標とした。
特に本研究では摩擦力によるエネルギーの散逸が必要かどうかを調べることにする。

11

第 2 章 剛体の運動の数値シミュレーション
2.1

剛体の形状

剛体の形状は回転楕円体とする。慣性主軸系での座標を (X1 , X2 , X3 ) とすると表面方程式は



X1
R1

2


+

X2
R2

2


+

X3
R3

2
=1

(2.1)

である。本論文では形状を変えても体積は変わらないようにする。即ち３軸の比 R1 : R2 : R3 が与えられたとき体積 V = 4π
3 R1 R2 R3 が半径 R00 の球の体積に等しくなるように R1 , R2 , R3 を決めることにする。与えられた比を R10 : R20 : R30 とすると

R0 =

R00
0
(R1 R20 R30 )1/3

(2.2)

R1

= R0 R10

(2.3)

R2

= R0 R20

(2.4)

R3

R0 R30

(2.5)

として

=

と決めることで実現される。実際このとき体積は4π
4π 3 0 0 0
4π
R1 R2 R3 =R R R R =R00
3
3 0 1 2 33となっている。質量は密度を ρ として
V =

M =Vρ

(2.6)

(2.7)

である。主慣性能率は

I1

=

I2

=

I3

=

1
M (R22 + R32 ) ,51
M (R32 + R12 ) ,51
M (R12 + R22 ) ,5で与えられる。
剛体系で (n1 , n2 , n3 )(単位ベクトル) を外向き法線とするような、表面の点の座標 (t1 , t2 , t3 ) は

12

(2.8)
(2.9)
(2.10)

t1= t01 S ,(2.11)

t2= t02 S ,(2.12)

t3=

t03 S ,

(2.13)

t01= R12 n1(2.14)

t02=

R22 n2(2.15)=

R32 n3(2.16)となる。このとき t01 , t02 , t03 は

t03であり、また1
S=p0
t1 n1 + t02 n2 + t03 n3
である。

13

(2.17)

2.2

剛体の運動方程式

2.2.1

剛体の運動の力学的記述

図 2.1: 剛体の運動

図 2.1 で L 系 (Laboratory frame) は実験室系であり、その成分は x, y, z で表すことにする。B 系~ は重心の座標である。
(body − fixed frame) は剛体系であり、その成分は 1, 2, 3 とラベルする。また R~ei , ~ej は(~ei · ~ej = δij =01(i 6= j)(i = j)(2.18)

を満たし、かつ~e1 × ~e2

= ~e3 ,

(2.19)~e2 × ~e3

= ~e1 ,

(2.20)~e3 × ~e1

= ~e2 ,

(2.21)

14

が成り立つ。

ω~ と ω1 , ω2 , ω3 の間の関係は

ω1= ~e1 · ω
~ ,(2.22)

ω2= ~e2 · ω
~ ,(2.23)

ω3= ~e3 · ω~

(2.24)でありω~ = ω1~e1 + ω2~e2 + ω3~e3(2.25)となる。よって下記の等式が成立する

ω
~ × ~e1 = ω3~e2 − ω2~e3 ,(2.26)

ω~ × ~e1 = ω1~e3 − ω3~e1 ,(2.27)

ω~ × ~e1 = ω2~e1 − ω1~e2 .

(2.28)

15

2.2.2

運動方程式~e˙1 , ~e˙2 , ~e˙3 は

~e˙1

= ω
~ × ~e1 ,(2.29)~e˙1= ω~ × ~e1 ,(2.30)~e˙1= ω~ × ~e1 ,(2.31)(~e˙1 )x

= ωy (~e1 )z − ωz e1y ,(2.32)(~e˙1 )x

= ωy (~e1 )z − ωz e1y ,(2.33)(~e˙1 )x

= ωy (~e1 )z − ωz e1y ,(2.34)であり、成分表示するととなる。また

ω̇1

=

ω̇2

=

ω̇3

=

~˙ CM
R
~˙ CM
V

I2 − I3
ω2 ω3 + N1 ,I1
I3 − I1
ω3 ω1 + N2 ,I2
I1 − I2
ω1 ω2 + N3 ,
I3

~CM ,
= V
1 ~
=
F ,
M

(2.35)(2.36)(2.37)(2.38)(2.39)~ は外力 (ここでは床との接点に床が働く力) の重心のまわりのモーメント、F~となる。このとき M は質量、N~ , F~ は R~ CM , V~CM , ~e1 , ~e2 , ~e3 , ω
は外力のベクトル和 (重力と床からの力)、RCM は重心の座標である。なお、N~の関数として与えられる。

16

2.3

剛体のエネルギー

剛体のエネルギーには１つ目は回転の運動エネルギーがある。これはErot =

I1 ω12 + I2 ω22 + I3 ω32
2(2.40)で表され、剛体が角速度 ω~ で回転しているときの運動エネルギーである。ここで I1 , I2 , I3 は主慣性能率とする。
２つ目は並進の運動エネルギーである。これはEtra =~CM )2x + (V~CM )2y + (V~CM )2z ]M [(V2(2.41)~CM は重心の速で表され、剛体が並進運動しているときの運動エネルギーである。ここで M は質量 、V度とする
３つ目は重力のポテンシャルエネルギーがある。これは、床上に寝かして置かれた場合を 0 として~ CM )z − (R~ CM )z.min ]Egra = M g[(R(2.42)~ CMで表され、回転して重心が上昇することによりこのポテンシャルエネルギーも増加する。ここで R
は重心の位置ベクトルとする。
この３つの定義は一般的であるが抗力のポテンシャルエネルギーについては何らかの仮定が必要である。
本研究では抗力 Fr を以下のように仮定した。

Fr

= −

z
z
∂Vz
= M ge− D (1+α D )∂z(2.43)

ここで D 床の軟らかさを決めるパラメーターは α は減衰を早めるパラメーター z 方向の減衰を早くzzさせるためガウシアン e− D (1+α D ) を使い図 2.2 のように z が大きくなると速やかに Fr が大きくなるようにした。

図 2.2: 床の抗力 Fr の z 依存性17なお、この Fr を与えるポテンシャル Vz は下式で与えられる。





1 
√
π e 4α1z√ Erf
Epot = Vz = −M gD
α
+
−1
2
D 2α
α
√

18

(2.44)

2.4

初期条件

図 2.3: 初期条件

図 2.3 のように対称軸が鉛直軸となす角度 θ を 88.5 °とし、角速度ベクトルは鉛直上向きにとり、その大きさ ω を 50 回転/秒としてこの状態から回転させる。すなわち、このとき単位ベクトル ~e1 , ~e2 , ~e3 、
角速度 ω の初期条件を以下のようにする。

(e~a )x = 1.0(e~a )y = 0.0(e~a )z = 0.0(2.45)(e~b )x = 0.0(e~b )y = cos θ(e~b )z = − sin θ(2.46)(e~c )x = 0.0(e~c )y = sin θ(e~c )z = cos θ(2.47)ω1= ω(~ω1 )z ,

(2.48)ω2= ω(~ω2 )z ,

(2.49)ω3= ω(~ω3 )z .

(2.50)~CM と重心の位置ベクトル R~ CM の初期条件を以下のようにする。
また重心の速度 V~ CM )x = 0(R~ CM )y = 0(R~ CM )z = na ta + nb tb + nc tc(R(2.51)~CM )x = 0(V~CM )y = 0(V~CM )z = 0(V(2.52)

ここで n1 , n2 , n3 は鉛直下向きの単位ベクトル ~n の 1, 2, 3 成分であり、t1 , t2 , t3 は (2.11) 式〜(2.17) 式で与えられる。~n を法線とする剛体の表面上の点の重心に対する変位ベクトルの 1, 2, 3 成分である。

19

2.5

摩擦力の効果を調べる

本研究では摩擦力の効果を調べることにした。
摩擦力の大きさ Ff を摩擦係数 µ と抗力の大きさ Fr で Ff = µFr と表し、動摩擦係数を µ = 0.5 とする。摩擦力の働く方向を人為的にずらしてみることにした。

図 2.4: 剛体と床との接点での力図 2.4 は剛体を上から見た図で、剛体は角速度 Ω で左回りに回転している。
~ を床に対する接点の速度、図中の「•」を剛体と床との接点とする。
V~ の逆方向 (図中の F~f の方向) に働く。本研究では、その向きを角度 η だけずらし、
通常は摩擦力は V
そのとき回転する剛体にどのような現象が起きるかを調べた。
図 2.4 に示したように η は −180 °≤ η ≤ 180 °の範囲で動かし、η の符号は、η > 0 では軸を寝かせる方向に摩擦力がかかり、η < 0 では軸を起こす方向に摩擦力がかかるものとする。

20

第 3 章 摩擦力の変化による回転直立の変化
3.1

動摩擦係数を変えた場合の影響

摩擦力が回転直立にとって重要な役割を果たしていることを確かめるため、動摩擦係数を µ = 0 として計算した結果を図 3.1 に示した。縦軸は重心の高さ (床の高さをゼロとした z 座標) であり、横軸は経過時間である。対称軸を水平にして床においた状態で重心の高さは 1.8cm である。回転開始後 10 秒経
過してもこの高さからの変化は認められない。したがって摩擦力がないと直立はおきないと思われる。
1.8177mu

1.8176R[cm]1.81751.8174
1.8173
1.8172
1.81710

2

4
6time[sec]

810

図 3.1: 剛体の重心の高さの時間変化： 摩擦力のない場合

2.5
2.4

R[cm]

2.3
2.2
2.1
2
case1
case2
case3

1.9
1.8

0

10

20
30
time[sec]

4050

図 3.2: 剛体の重心の高さの時間変化： 摩擦係数の異なる場合の比較

21

図 3.2 には摩擦係数が µ = 0.01 (case 1)、µ = 0.5 (case 2)、µ = 0.7 (case 3) の３通りの場合を示した。縦軸と横軸は図 3.1 と同じである。卵の殻と同じオーダーの µ = 0.5 の場合と比較して、摩擦係数の非常に小さい場合 (µ = 0.01) は直立までに要する時間が一桁以上長くなる。この結果は摩擦が回転直立の原因であると考えれば理解できる。一方、摩擦係数をある程度以上大きくした場合 (µ = 0.7) も、や
はり所要時間が長くなる。この場合は直立の起き方が他の２つの場合と異なり、長い時間をかけて緩慢に直立して行く。

22

10
9
8

t[sec]

7
6
5
4
3
2
1
0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

mu
図 3.3: 回転直立に要する時間と動摩擦係数の関係

図 3.3 は回転直立に要する時間を動摩擦係数の関数としてプロットしたものである。直立の判定条件
としては「対称軸と鉛直方向のなす角度の余弦の絶対値が 0.9 以上であること」を課した。このグラフを見ると、µ < 0.53 では直立所要時間は摩擦係数 µ の増加とともに減少するが、µ = 0.53 を境にして挙動が一変し、µ > 0.53 では、µ とともに急速に増加するようになる。
この挙動の変化の原因として考えられることは、摩擦力が強くなり過ぎて回転の角速度が落ちてしまうことであるが、その当否は本論文では調べず将来の研究に譲ることにする。
これらの計算結果から、摩擦係数がある有限の値 (計算で用いた剛体の形状や回転の初期条件に対しては µ = 0.53) をとる場合に最も速やかに直立することがわかった。また、摩擦係数がゼロだと長時間経過しても直立しないこともわかった。しかしながら、これらの結果だけからでは、摩擦力が回転直立の不可欠の要因であると断言することはできない。
確かに Moffatt-下村の理論では摩擦力が回転直立を引き起こしている。しかしその結論は、摩擦がない場合に実現する定常的な歳差型の運動が緩慢に変化するという仮定の下で導き出されたものである。
実際に石ころなどの剛体を手で回してみると、そのような仮定にあてはまらない状況下でもほとんどの
場合に回転直立が起きることがわかる。例えば大抵の石ころは軸対称ではないので非常に激しい跳躍を伴う振動を起こしてから直立するが、跳躍運動は Moffatt-下村理論の適用範囲外である。
あるいは次のような状況も想像できる。摩擦のない滑らかな面の上で回転する場合でも、面に凹凸があるとしたらどうなるだろうか。長い対称軸を水平にして回転している状態は、対称軸を鉛直にして回転している状態より擾乱を受けやすいはずである。したがって、寝て回転している時間より立って回転している時間のほうが長くなる可能性があるのではないだろうか。もしそういう実例を見付けることができれば、それは保存系において回転直立に相当する現象であると言えるだろう。
このような回転直立現象をより広い立場からとらえるための研究の第一歩として、以下では、摩擦力の働き方を人為的に通常とは違えてみて、摩擦力の何が回転直立現象にとって本質的な役割を果たしているのか、通常の摩擦力以外の力でも直立がおきることはないのか、を調べることにする。

23

3.2

摩擦力の方向を人為的にずらした場合の影響

動摩擦力の働く方向は接地点の速度の方向とは正反対の方向であるが、この節では、この方向を人為的に違う方向にずらしてみたとき回転直立現象はどのような影響を受けるかを調べることにする。このずれの角度をηで表す。ηのとる値の範囲は −180◦ ≤ η ≤ 180◦ とし、η > 0 は軸を倒すような力の成分を持つ場合を、η < 0 は軸を起こすような力の成分を持つ場合を表すことにする。動摩擦係数は µ = 0.5に固定したまま角度 η を 0 からずらした場合の影響を調べた結果を、|η| < 90◦ の場合、|η| = 90◦ の場
合、|η| > 90◦ の場合に分けて以下に述べる。

3.2.1|η| < 90◦ の場合

この場合は、|η| の増大とともに摩擦力の「接地点の速度方向の成分」の大きさは cos η に比例して小さくなるが、向きは速度の反対方向のままである。|η| が増大すると、本来の摩擦力の大きさが減少すると同時に、「接地点の速度に垂直な方向の成分」の力が sin η に比例して発生する。この力は η > 0 (< 0)の場合は対称軸を倒す (起こす) ような方向を向いている。

2.5
2.4

R[cm]

2.3
2.2
2.1
2
case1
case2
case3

1.9
1.8

0

2

4
6
time[sec]

810

図 3.4: 剛体の重心の高さの時間変化:η = 0◦ , 30◦ , 60◦ の場合

図 3.4 に剛体の重心の z 成分 (重心の床からの高さ) の時間変化を、η が 0◦ (case 1)、30◦ (case 2)、60◦ (case 3) の３つの場合について示した。対称軸を水平 (鉛直) にして回転している状態では重心の高さは 1.82 cm (2.42cm) である。３つのケースともに、回転直立が起きているが、直立に要する時間は|η| が小さいほど短いことがわかる。

24

2.6
2.5
2.4
R[cm]

2.3
2.2
2.1
2
case1
case2
case3

1.9
1.8

0

2

4
6
time[sec]

810

図 3.5: 剛体の重心の高さの時間変化：η = 0◦ , ±45◦ の場合

3.1 節で述べた動摩擦係数を変化させた場合の計算結果から、|η| の増大とともに回転直立が遅くなる
ことの原因のひとつとして摩擦力の「接地点の速度方向の成分」の減少があることは間違いない。しかし、垂直方向の成分の増加も同時になんらかの影響を与えている可能性がある。これを調べたのが |η|の正・負で結果がどう変わるかを示した図 3.5 である。case 1 は η = 0◦ 、case 2 は η = 45◦ 、case 3 は

η = −45◦ の場合を示し、縦軸と横軸は図 3.4 と同じである。この図より、η = 0◦ の場合と比較して、
η = −45◦ の場合は直立に要する時間は 0.8 秒程度にしか増加していないが、η = 45◦ の場合は 2.0 秒に大幅に増加していることがわかる。η が負の場合には速度に垂直な方向の力は対称軸を起こす向きに働くことを考えると、この力が回転直立を促進していると考えるのが自然である。
なお、η = −45◦ の場合には直立直前に起きるジャンプの振幅が大きくなっているが、これも、速度に垂直な方向の力に帰することができる。即ち、この力が直立を助けるため、直立するペースが速くな
り過ぎてジャンプが早い段階から起き始めるからだと解釈できる。

25

3.2.2|η| = 90◦ の場合

前節の結果から、速度に垂直な方向の力だけでも回転直立が起きるのではないだろうかという疑問が湧く。この場合、力は仕事を行わないので力学的エネルギーは保存されることになる。もし保存系で回転直立が起きる例が見付かれば、回転直立現象の理解にとって画期的な発見となるであろう。

1.81765
1.81761.81755

R[cm]1.81751.81745
1.8174
1.81735
1.8173
1.81725
case1
case2

1.8172
1.817150

2

4
6time[sec]

810

図 3.6: 剛体の重心の高さの時間変化： η = ±90◦ の場合残念ながら、|η| = 90◦ の場合には、系の力学的発展を長時間にわたって追跡しても、回転直立は全く起きなかった。図 3.6 に示したのは、case 1 が η = 90◦ の場合, case 2 が η = −90◦ の場合である。縦軸と横軸は図 3.4 と同じである。どちらの場合も重心はほとんど上昇していない。また、+90◦ と −90◦との差もわずかしかない。
コマや地球などが行う定常的な歳差運動の最も単純な理論では、回転軸を倒したり起こしたりする方向の力は、回転軸が鉛直方向となす角度を変える効果は持たず、その代わりに回転軸を鉛直軸の回りに歳差運動させる効果を持つ。この観点からは、速度に垂直な力では回転直立が起きないことは常識的な結果であるかもしれない。しかしながら、通常の摩擦力と同時に働く場合には、速度に垂直な力が回転
直立を促進しうるという結果はかなり意外性のある発見であると言えるのではないだろうか。
この回転直立促進のメカニズムの理解法としては、次のような説明が考えられる。まず、歳差運動の一般的性質として [7]、η > 0 では、対称軸を倒す向きの力のモーメントが働くので、歳差運動の角速度Ω は増加することになる。実際には摩擦による減衰があるので、Ω は減少するが、その減少率が小さくなるはずである。したがって、もしジェレット定数 ∝ Ωh が保存されるなら h (重心の高さ) の増加は遅れるはずである。一方、η < 0 では、対称軸を起こす向きの力のモーメントが働くので、Ω の減少は促進されることになる。したがって、ジェレット定数 ∝ Ωh が保存されるなら h の増加は速められるはず
である。

26

3.3

|η| > 90◦ の場合

さらに、|η| が 90◦ を越える場合についても、理論的興味から調べてみた。図 3.7 は η = 150◦ の場合である。縦軸と横軸は図 3.4 と同じである。重心の高さは、対称軸を水平にして寝かせた状態にあたる

1.8cm から対称軸を鉛直にして立てた状態の 2.4cm に近い値にまで上昇していることが見て取れる。これは回転直立とほぼ同じ挙動である。摩擦力の向きが速度と鋭角をなす場合といえば、エネルギーが散逸でなく増加する場合であるから、これはエネルギーの散逸を伴う卵の回転直立とは全く違った力学的プ
ロセスである。それにもかかわらず、やはり回転直立は起きるとすれば、まさに画期的な発見であろう。

2.35
2.3
2.25
2.2
R[cm]

2.15
2.1
2.05
2
1.95
1.9
1.85
1.8

0

2

4
6
time[sec]8

図 3.7: 剛体の重心の高さの時間変化： η = 150◦ の場合

2710しかし残念ながら、この運動を詳しく調べてみると、回転直立とは別種の運動であることがわかった。
図 3.8 は重心の x 座標と y 座標の時間変化を示したものである。この図から剛体は x と y の間の方向へと加速度運動していることがわかる。この水平方向への加速度は、g を重力加速度、µ を動摩擦係数として |gµ sin η| ≒ 0.4g として見積もった結果と粗くは一致している。結論として、この場合の重心の高さの変化は、回転直立の結果ではなく、接地点の速度と鋭角をなす方向への「摩擦力」が水平方向の並進運動および剛体の対称軸に垂直な軸の回りの回転を引き起こした結果、図 3.9 に描いたような転がり運動が起きたためだとして理解できる。
14000

x
y

12000
10000
R[cm]

8000
6000
4000
2000
0
-2000

0

2

4

6

8

10

time[sec]

図 3.8: 剛体の重心の水平面上での移動の様子: η = 150◦ の場合

図 3.9: 剛体の転がり

28

3.4

η を変化させたときの直立時間の変化
100t[sec]

10

1

0.1

-80

-60

-40

-20

0
20
eta[dec]

40

60

80

図 3.10: 回転直立に要する時間と摩擦力の働く方向の関係

図 3.10 を見ると |η| が大きいと直立するのが遅く、|η| が小さくなるにつれ直立するのが早くなりη = −17 °で一番早く直立することが分かる。また η ≥ 90 °では直立しない。
また η > 0 と η < 0 を比べると η < 0 の方が t が小さいことが分かる。

29

3.5

時間反転解

やはり理論的興味から η = 180◦ という極端な場合も計算してみた。この場合の摩擦力は速度と同じ方向を向いているが、これは、通常方向の摩擦力が働く場合について、力学的状態を時間を遡って求めて行く場合に相当しているので、時間反転解と呼ぶことにする。
時間反転すると、回転直立の逆の現象が起きるはずである。即ち、図 3.11 に示したように、直立して回転している剛体が時間が経つと寝て回転する状態になるはずである。そこで、直立して回転している状態を初期状態にとって計算を行った。結果は図 3.12 に示したとおりである。縦軸と横軸は図 3.4 と同じである。この図を見ると経過時間が 0.1 秒程度までは重心が下がるのだが、その後、激しくジャンプをする期間があり、その後、重心は約 2.3cm まで上昇して止まることがわかる。この高さは 3.3 節で見た転がり運動と同じ高さであり、この場合もやはり、並進運動と「対称軸に垂直な軸の回りの回転運動」
が同時に起きる転がり運動が起きているのだと予想される。実際、図 3.13 を見ると水平方向への並進運動が起きていることが確認できる。図 3.13 の縦軸と横軸は図 3.4 と同じである。結論としては、回転直立を時間的に遡ることは 0.1 秒程度しかできず、その後は、図 3.14 に示したようなエネルギーの高い乱れた運動状態へ移ってしまったということである。

図 3.11: 時間反転解

2.8
2.7
2.6
R[cm]

2.5
2.4
2.3
2.2
2.1
2

0

2

4
6
time[sec]8

図 3.12: 時間に対する重心座標の z 成分の変化

30

10

0

R[cm]

-5000
-10000-15000-20000-25000

x
y
0

2

4

6

8

10

time[sec]
図 3.13: 時間に対する重心座標の x, y, z 成分の変化

図 3.14: 剛体の不安定さ

散逸のある力学系の時間反転に関しての一般的を言えば、エネルギーは下に有界なので、エネルギーが散逸する場合は、状態は時間発展によりエネルギー最小の解へと収束していく。しかし、エネルギーは上には有界でなく、いくらでもエネルギーの大きな状態が存在し、無限の体積の位相空間を占めているので、時間を逆行してエネルギーを増加させる場合には、限りなくエネルギーの大きな状態に移っていくばかりであって、特定の解に収束していく確率はゼロなのである。

31

第 4 章 結論

図 4.1: 摩擦の方向

本研究では摩擦力の働き方を変えることで剛体の回転直立現象にどのような影響があるかを調べた。
まず１つ目に摩擦力がないと直立は起きない。２つ目に、図 4.1 のように摩擦力の方向と床との接点の速度ベクトルの反対方向の角度を η とすると、|η| < 90 °のとき直立し、|η| = 90 °のときは直立しない。また |η| > 90 °のときは並進運動が起き、転がりながら進み、不安定である。
従って剛体の進む方向と逆の方向に摩擦力が働かないと剛体は直立しないが、それと垂直な方向の力
は直立を助けたり妨げたりするということが分かった。
タイトルに掲げた「保存系での回転直立」は η = 90 °の場合にあたるが、この場合は直立する例を見つけることができなかった。

32

参考文献
[1] 戸田盛和:「回転する卵はなぜ直立する」科学 (岩波書店),
72,932-939(2002).
[2] C.M.Braams,Physica, 18,503-514(1952).
[3] N.M.Hugenholtz,Physica, 18,515-527(1952).
[4] H.K.Moffatt and Y.Shimomura: Spinning eggs - a paradox resolved,
Nature 416,385-386(2002).
[5] 村井興平:「卵を回すとなぜ立つか」,
福井大学工学部 物理工学科 卒業研究 (2004).
[6] 大森英胤:「卵を回すとなぜ立つか２」,
福井大学工学部 物理工学科 卒業研究 (2005).
[7] 小出昭一郎:「力学」岩波書店,
pp.90-165(1987).
[8] 下村裕:「立ち上がる回転ゆで卵の解」パリティ,
18,52-56(2003).
[9] 日本物理学会年会講演,
27pYC10,(野田,2005).
[10] 戸田盛和:「力学」岩波書店,
pp.144-188(1982).

33

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

34

付録#include <stdio.h>#include <math.h>#include <limits.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] */double beta; // angle [degree] of the friction force
double cos_beta, sin_beta;
int kaiten5(){
Bvec radius_ratio = { 3.00, 3.00, 4.0}; /* relative lengths of principal axes */double R00=2.0 ;/* radius for shperical shape [cm] */double density=1;/* density [g/cm^3] */35

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;/* maximum time [sec] to calculate */double max_itime;/* the number of time steps */int
mprint;/* period to print the mechanical state */int
mprintTmp;
Bvec moi;/* moment of inertia */Bvec omg0;/* initial angular velocity vector, B-frame */Bvec omg;/* angular velocity vector, B-frame */Lvec omgL;
double omgsize;
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 */36

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;
FILE *log_res5;
FILE *log_res6;
Lvec angmom; /* angular momentum in L-frame */double angmomsize;
double F1,F2,F3;/* log_res1 = NULL;*//* log_res2 = NULL;*//* log_res3 = NULL; *//* log_res4 = NULL; *//* log_chk1 = NULL;*/

log_res1=fopen("rbm1.g","w");log_res2=fopen("rbm2.g","w");log_res3=fopen("rbm3.g","w");log_res4=fopen("rbm4.g","w");log_chk1=fopen("rbm5.g","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; /* volume of the rigid body [cm^3] */mass=density*volume;
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, max_t, beta[deg]");scanf("%lf %lf %lf",&dt, &max_t, &beta);
fprintf(stderr,"check dt=%e
max_t=%e beta=%e\n",dt,max_t,beta);
cos_beta=cos(beta/180*pi);
sin_beta=sin(beta/180*pi);if(max_t/dt > INT_MAX){
fprintf(stderr,"max_t/dt=%e > %e",max_t/dt,(double)INT_MAX);37

exit(1);}max_itime=ceil(max_t/dt);
mprint=0.005/dt;
mprintTmp=max_itime/ 5000.0; if(mprint>mprintTmp) mprint=mprintTmp;
mprintTmp=max_itime/50000.0; if(mprint<mprintTmp) mprint=mprintTmp;
if(mprint<1) mprint=1;
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);
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"
1
2
3
4
5
6
7
8
9
10
,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;

38

/*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(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
);

// 1
// 2
// 3
// 4
// 5
// 6
// 7

/*
}
if(itime >= max_itime) break;forces(&rc, &vc, &omg, &ea, &eb, &ec, &frc, &trq);

//
//
//

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"
1
2
3
4
5
6
,t,
angmom.x, angmom.y, angmom.z, angmomsize, omgsize
78910,eng_rot, eng_tra, eng_gra, eng_tot
11
12
13
14
15
16
17
18
19
39

,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);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 ;k1ax=dt*(omgL.y*ea.z-omgL.z*ea.y);k1ay=dt*(omgL.z*ea.x-omgL.x*ea.z);k1az=dt*(omgL.x*ea.y-omgL.y*ea.x);
k1bx=dt*(omgL.y*eb.z-omgL.z*eb.y);k1by=dt*(omgL.z*eb.x-omgL.x*eb.z);
k1bz=dt*(omgL.x*eb.y-omgL.y*eb.x);k1cx=dt*(omgL.y*ec.z-omgL.z*ec.y);k1cy=dt*(omgL.z*ec.x-omgL.x*ec.z);k1cz=dt*(omgL.x*ec.y-omgL.y*ec.x);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;
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);omgL.x = omgm.a*eam.x + omgm.b*ebm.x + omgm.c*ecm.x ;40

omgL.y = omgm.a*eam.y + omgm.b*ebm.y + omgm.c*ecm.y ;omgL.z = omgm.a*eam.z + omgm.b*ebm.z + omgm.c*ecm.z ;k2ax=dt*(omgL.y*eam.z-omgL.z*eam.y);k2ay=dt*(omgL.z*eam.x-omgL.x*eam.z);k2az=dt*(omgL.x*eam.y-omgL.y*eam.x);k2bx=dt*(omgL.y*ebm.z-omgL.z*ebm.y);k2by=dt*(omgL.z*ebm.x-omgL.x*ebm.z);k2bz=dt*(omgL.x*ebm.y-omgL.y*ebm.x);k2cx=dt*(omgL.y*ecm.z-omgL.z*ecm.y);k2cy=dt*(omgL.z*ecm.x-omgL.x*ecm.z);k2cz=dt*(omgL.x*ecm.y-omgL.y*ecm.x);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;
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);omgL.x = omgm.a*eam.x + omgm.b*ebm.x + omgm.c*ecm.x ;omgL.y = omgm.a*eam.y + omgm.b*ebm.y + omgm.c*ecm.y ;omgL.z = omgm.a*eam.z + omgm.b*ebm.z + omgm.c*ecm.z ;k3ax=dt*(omgL.y*eam.z-omgL.z*eam.y);k3ay=dt*(omgL.z*eam.x-omgL.x*eam.z);k3az=dt*(omgL.x*eam.y-omgL.y*eam.x);

41

k3bx=dt*(omgL.y*ebm.z-omgL.z*ebm.y);k3by=dt*(omgL.z*ebm.x-omgL.x*ebm.z);k3bz=dt*(omgL.x*ebm.y-omgL.y*ebm.x);k3cx=dt*(omgL.y*ecm.z-omgL.z*ecm.y);k3cy=dt*(omgL.z*ecm.x-omgL.x*ecm.z);k3cz=dt*(omgL.x*ecm.y-omgL.y*ecm.x);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;
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);k4z=dt*(f3*omgm.a*omgm.b+fn3*trq.c);omgL.x = omgm.a*eam.x + omgm.b*ebm.x + omgm.c*ecm.x ;omgL.y = omgm.a*eam.y + omgm.b*ebm.y + omgm.c*ecm.y ;omgL.z = omgm.a*eam.z + omgm.b*ebm.z + omgm.c*ecm.z ;k4ax=dt*(omgL.y*eam.z-omgL.z*eam.y);k4ay=dt*(omgL.z*eam.x-omgL.x*eam.z);k4az=dt*(omgL.x*eam.y-omgL.y*eam.x);k4bx=dt*(omgL.y*ebm.z-omgL.z*ebm.y);k4by=dt*(omgL.z*ebm.x-omgL.x*ebm.z);k4bz=dt*(omgL.x*ebm.y-omgL.y*ebm.x);k4cx=dt*(omgL.y*ecm.z-omgL.z*ecm.y);k4cy=dt*(omgL.z*ecm.x-omgL.x*ecm.z);k4cz=dt*(omgL.x*ecm.y-omgL.y*ecm.x);42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);
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;
eb.y=eb.y+kby;
eb.z=eb.z+kbz;
43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;}
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) {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 ;44

mez.a = - ea->z ; mez.b = - eb->z ;normal_point(&mez, &tp);

mez.c = - ec->z ;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 ;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)));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*cos_beta - vtp.y*sin_beta);
frctp.y = fct*(vtp.x*sin_beta + vtp.y*cos_beta);
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 ;return 0;}
int normal_point(Bvec *n, Bvec *tp){45

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;
return 0;}
void init2(Lvec *ea, Lvec *eb, Lvec *ec, Bvec *omg0, Lvec *rc, Lvec *vc){
double b;
double theta=(88.5/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;
eb->x=0.0;
ec->x=0.0;

ea->y=0.0;
eb->y=cos(theta);
ec->y=sin(theta);

ea->z=0.0;
eb->z=-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);}46