転がるコインの力学的シミュレーション
2009SE250 澤井優仁
指導教員：杉浦洋
重心 O の速度ベクトルを v とし，コインの半径を a とす

はじめに1る．以上の三式を用いると，速度ベクトルはコインの転がる様子には様々なものがある，この様子についてコンピュータのシミュレーションで詳しく調べる．

v = aθ̇k − aϕ̇i + aψ̇ cos θi

コインの運動を支配しているのは剛体力学である．本研究= −a(ϕ̇ + ψ̇ cos θ)i + aθ̇k(7)では剛体力学を詳しく学び，それによりコインの状態を表す微分方程式を導出する．その微分方程式を数値的に解くことにより，コンピュータ上にコインの運動を再現する．
国枝 [1] を参考に運動方程式を導いた．[1] では加速度計の力学を用いるが，我々は慣性系の力学のみで導出を行った．全ての式は Mathematica の数式処理機能を用いて検2.2運動方程式の導出

i，j ，k の時間微分を絶対座標で表すと
−ψ̇ sin ψdi =  ψ̇ cos ψ dt
0

算した．

= ψ̇ cos θj − ψ̇ sin θk,2

三次元空間でのコイン

2.1

パラメータの設定(8)θ̇ sin θ sin ψ − ψ̇ cos θ cos ψdj =  −θ̇ sin θ cos ψ − ψ̇ cos θ sin ψ dt
θ̇ cos θ

水平な床上を滑らない転がるコインの運動を解析する．
コインは O を中心とする半径 a の円盤で，その厚さは無= −ψ̇ cos θi + θ̇k,

視する．直交座標系 O − XY Z は，XY 軸は水平，Z 軸は(9)θ̇ cos θ sin ψ + ψ̇ sin θ cos ψdk =  −θ̇ cos θ cos ψ + ψ̇ sin θ sin ψ dt−θ̇ sin θ

鉛直の絶対座標である．また，O − xyz はコインに伴って動く直交座標系で，z 軸はコインの裏面から表面に抜ける垂直軸，x 軸は常に水平と取る．y 軸は座標系 O − xyz が右手系であることから自動的に決まる．Z 軸から z 軸への= ψ̇ sin θi − θ̇j.

角度をθ，X 軸から x 軸への角度をψ ，コインの z 軸周り(10)(7)，(8)，(9)，(10) より，速度を微分するとの回転角をϕとする．

x，y ，z 軸方向の単位ベクトルを i，j ，k として，それぞv̇ = (−aϕ̈ − aψ̈ cos θ + 2aθ̇ψ̇ sin θ)i+{−aθ˙2 − aψ̇ cos θ(ϕ̇ + ψ̇ cos θ)}j

れ絶対座標の数ベクトルで表すと．

()
cos ψ+{aθ̈ + aψ̇ sin θ(ϕ̇ + ψ̇ cos θ)}k(11)sin ψi=,0
ゆえに，重心の並進運動の方程式は以下のようになる．
()()− cos θ sin ψsin θ sin ψcos θ cos ψ− sin θ cos ψ . (1)  −ma(ϕ̈ + ψ̈ cos θ − 2θ̇ψ̇ sin θ) = Fx ,j=, k=sin θcos θ−ma{θ˙2 + ψ̇ cos θ(ϕ̇ + ψ̇ cos θ)} = Fy − mg sin θ, (12)ma{θ̈ + ψ̇ sin θ(ϕ̇ + ψ̇ cos θ)} = Fz − mg cos θ.
O − xyz 系の絶対座標系に対する角速度ベクトルΩはΩ = θ̇i − ψ̇ sin θj + ψ̇ cos θk.

(2)

すると

コインの絶対座標系に対する角速度ベクトルω は

ω = θ̇i + ψ̇ sin θj + (ϕ̇ + ψ̇ cos θ)k.

(3)

x，y ，z 軸周りのコインの慣性モーメントを A，A，C とす
ると，コインの角運動量 L は
L = Aθ̇i + Aψ̇ sin θj + C(ϕ̇ + ψ̇ cos θ)k.

(4)

系に入る外力は，鉛直下向きの重力と，コインと床との間に働く拘束力の二つである．コインの質量を m，拘束力を F とすると，二つの外力は−mge3 = −mg(sin θj + cos θk),
F = Fx i + Fy j + Fz k.

(4)，(8)，(9)，(10) より，同様にコインの角運動量を微分

(5)
(6)

L̇ = {Aθ̈ + Aψ˙2 sin θ cos θ + C ψ̇ sin θ(ϕ̇ + ψ̇ cos θ)}i+{Aψ̈ sin θ + 2Aθ̇ψ̇ cos θ − C θ̇(ϕ̇ + ψ̇ cos θ)}j+{C(ϕ̈ + ψ̈ cos θ − θ̇ψ̇ sin θ)}k(13)一方，外力のモーメントを N とするとL̇ = N = −aj × F = Fx ak − Fz ai(14)

重心回りの回転運動の方程式は，

 Aθ̈ + Aψ˙2 sin θ cos θ + C ψ̇ sin θ(ϕ̇ + ψ̇ cos θ) = −Fz a,(15)Aψ̈ sin θ + 2Aθ̇ψ̇ cos θ − C θ̇(ϕ̇ + ψ̇ cos θ) = 0,C(ϕ̈ + ψ̈ cos θ − θ̇ψ̇ sin θ) = Fx a.

(12)，(15) より，拘束力を消去して
(C + ma2 )(ϕ̈ + ψ̈ cos θ − θ̇ψ̇ sin θ) = ma2 θ̇ψ̇ sin θ,Aψ̈ sin θ + 2Aθ̇ψ̇ cos θ − C θ̇(ϕ̇ + ψ̇ cos θ) = 0,(16)(A+ ma2 )θ̈ − Aψ˙2 sin θ cos θ+(C+ma2 )ψ̇ sin θ(ϕ̇ + ψ̇ cos θ) + mga cos θ = 0.

3.2

ふらつきを伴う運動

以下の初期値を指定する．それ以外は 0 とする．

θ=π,2.001

θ̇ = 3.07,ψ̇ = −0.003,以上の計算式はMathematicaでプログラムを組んで検算

ϕ̇ = 1.805,
π2.001z = a sinした。

Mathematica による数値実験

3三次元でのコインの転がるシミュレーションを Mathematica 上に実現し，数値実験を行った．重心の絶対座標を
r = (x, y, z) とし，微分方程式 (7)，(16) を連立させて解く．
m=1，a=1，A=ma2 /4=1/4，C=ma2 /2 = 1/2 とする．

3.1

歳差運動

いくつかの実験の中から，2 つの典型例を紹介する．いずれの例でも，検算のためにエネルギーを計算したが，計算精度の範囲内で保存された．以下の初期値を指定する．
それ以外は 0 とする．

θ=π,
3ϕ̇ =

図 3 コインの運動モデル

3π,2ψ̇ = −0.868,

z = a sinπ343211234567-1-2-3

図 4 変数の値(図 4) はθ̇，ϕ̇，ψ̇ のグラフであり，それぞれ一定の周期で
図 1 コインの運動モデル変化をしているため，このコインの運動 (図 3) は周期的な運動を行っていることが分かる．

44終わりに
本研究では，コインの状態を表す微分方程式を導出した．

Mathematica でプログラムを組み，その微分方程式を数3値的に解くことにより，コンピュータ上にコインの運動を2再現した．初期値によって周期的に円運動を行うコインの1動きを再現したり，ふらつきがありながらも周期的な運動
1234567を行うコインの動きを再現することができた．

参考文献
図 2 変数の値(図 2) の変数は上からϕ̇，θ̇，ψ̇ である．この 3 変数の値は一定となっており，θ̇ の値は常に 0 であることから，コインはぶれがなく動いていることが分かる．コインの運動 (図 1) は周期的なものとなった．

[1] 国枝正春：『6 運動自由度系の力学』.
http://ds0.cc.yamaguchiu.ac.jp/ tsaito/Advanced/PDF/advibration01.pdf.