11.自由な剛体

剛体の運動を求めることは、重心位置および回転行列の時間発展、それぞれ xg(t),R(t) を求めることに他ならない(回転行列 R の代わりにオイラー角 θ などでもよい)。これらは、運動方程式(12)、あるいはそれを数値計算に適した形に変形した式(21)により計算できる。ただし、予め、剛体の重心位置と慣性モーメントを計算しておく必要がある(式(17))
この章のシミュレーション(クリックで計算開始)

自由な速度 Ω に対する運動方程式(2)が欲しい

剛体の運動 xg(t),R(t) を計算したい。前章では、剛体の自由な速度として、「剛体上の基準点 x0 の速度 ˙x0 」と「角速度 ω 」という自然なものが存在することを見た。これらをまとめて Ω で表すことにする: Ω=[˙x0ω](1) 剛体の運動を計算するには、この自由な速度 Ω に対する運動方程式 ˙Ω=[](2) を求めればよい。
剛体の運動方程式(2)を得るには、第9章の議論に従えばよい。また、式(2)を用いて実際の計算を行う方法は自明ではないので、多少議論が必要である。そこでこの章では、以下の2つの節に分けて議論を行う: 11.1:11.2:
なお、式(2)の代わりに、自由な座標 x0,θ に対する運動方程式(2階微分方程式)を考えてもよい。しかし、 ω=˙θ となるような回転の自由な座標 θ は存在しないので(前章の【10.3-注1】)、特定の θ に限定した議論をすると一般性を欠くことになる。一方、式(2)を求めておけば、特定の θ を選んだ場合の運動方程式もすぐに導けるのでこれで十分である。実際、 ω との関係式は ω=A˙θ と言う形に書けるので(例えば前章の【10.3-注2】)、この関係式を式(2)に代入すれば、 ¨θ に対する方程式になる。

11.1 剛体の運動方程式の導出

この節では、剛体の運動方程式(12)と、それを数値計算に適した形に変形した式(21)を導く。

質点要素の速度 ˙X と自由な速度 Ω の関係式(5)が分かればよい

式(2)を求める方法は、自由な座標による運動方程式(第9章の【9.1-注1】)を求めた時と同様である。まず、全ての質点要素 xi に対する運動方程式は ⎢ ⎢m1m2⎥ ⎥M⎢ ⎢¨x1¨x2⎥ ⎥¨X=⎢ ⎢f1f2⎥ ⎥F+⎢ ⎢fc1fc2⎥ ⎥Fc(3) となる。
第8章の8.2節でみたように、拘束力 Fc は、ダランベールの原理により、拘束面と垂直である。即ち、(拘束が時間に依存しない場合) Fc は、拘束条件を満たす全ての速度 ˙X と「直交」する: ˙XTFc=0(4) ˙X は、自由な速度 Ω を用いて ˙X=AΩ(5) の形で書けるのでA は何らかの行列)、式(4)は、 ATFc=0 となる。従って、式(3)に AT を左乗すると、 Fc が消える: ATM¨X=ATF(6) 後は、この式に、式(5)を代入して ¨X を消去すれば、 ˙Ω に対する方程式の形にできる。
実際に代入を行うと、以下のようになる: ddt(μΩ)=ATF+˙ATMAΩμATMA(7) この式は、 ˙Ω=[] の形で書くこともできるが、 μΩ の微分の形で書いているのは、右辺第2項が消えて簡単になることが後で分かるからである(式(11))

質点要素の速度 ˙X を、自由な速度 Ω で表す:式(9)

この方針に従って、質点要素の速度 ˙X を、自由な速度 Ω (式(1))で表す。まず、 i 番目の質点要素の速度 ˙xi は、相対座標 xi0=xix0 を用いて ˙xi=˙x0+˙xi0=˙x0+ω×xi0=[1xi0×]Ω(8) となる(ベクトルのクロス x× の定義は第10章の【10.2-注2】)。これを ˙X に代入すれば、 ˙XΩ との関係式が得られる: ˙X=⎢ ⎢˙x1˙x2⎥ ⎥=⎢ ⎢1x10×1x20×⎥ ⎥ΩAΩ(9)

基準点を重心(10)に取った時の運動方程式:式(11)

後は、運動方程式(7)に式(9)の A を代入して、各項を計算していくだけである。実際の計算を行うに当たって、基準点 x0 を重心 xg に取ることにしよう:mmi は剛体の全質量) x0=xg=1mimixi(10) 重心を基準に取る理由は、剛体内の拘束力は作用・反作用の法則(第5章の【5.1-注1】)を満たすので、重心の速度 ˙xg は拘束力の影響を受けない(同【5.2-注1】)からである。つまり、拘束力は ˙xg を変化させず角速度 ω のみに影響を与えるようになるため、運動方程式が簡単になるわけである。
式(10)のもとで、運動方程式(7)の各項を変形していけばよく、その結果は以下の【11.1-注1】のようになる。それらを代入すると、以下の結果が得られる:xig=xixg ddt(μΩ)=[ifiixig×fi]μ[m00imi(xig2xTigxTig)](11) この式を見ると、外力 fi0 の場合、右辺が消えるので、 μΩ=const. となることが分かる。式が簡単になると式(7)で述べたのは、このことである。また、 μ がブロック対角行列になっているのは、基準点を x0=xg としたおかげである。

【11.1-注1】運動方程式(7)の各項

運動方程式(7)の各項を、式(9)の A 、および、 x0=xg のもとで計算すると、以下のようになる:xig=xixg μ=ATMA=[m00imi(xig×)2]ATF=[ifiixig×fi]˙ATMAΩ=0

導出

第1式: μATMA=[11x1g×x2g×]⎢ ⎢m1m2⎥ ⎥⎢ ⎢1x1g×1x2g×⎥ ⎥=[imiimixig×imixig×imi(xig×)2]‖ ‖ ‖ ‖ ‖ ‖imixig=imi(xixg)=m(xgxg)=0=[m00imi(xig2xTigxTig)] 赤字部分がうまく消えるのは、重心を基準にとったからである。また、緑字部分では、簡約公式 a×b×=baTaTb (第10章の【10.2-注2】)を用いてクロスを消去した。
第2式: ATF=[11x1g×x2g×]⎢ ⎢f1f2⎥ ⎥=[ifiixig×fi]
第3式: ˙ATMAΩ=[00˙x1g×˙x2g×]⎢ ⎢m1m2⎥ ⎥⎢ ⎢1x1g×1x2g×⎥ ⎥Ω=[00imi˙xig×imi˙xig×xig×]Ω‖ ‖ ‖ ‖ ‖ ‖imi˙xig=0=[000imi˙xig×xig×ω]‖ ‖ ‖ ‖˙xig×xig×ω=˙xig×(ω×xig)=˙xig×˙xig=0=0

剛体の運動方程式(11)の分離:式(12)

得られた結果をまとめておこう。式(11)を、重心速度 ˙xg に対するものと、角速度 ω に対するものに分けて書くと、以下のようになる: m¨xg=fddt(Igω)=τg(12) ただし、赤字部分は以下のように定義される: m=imif=ifiIg=imi(xig2xTigxTig)τg=ixig×fi(13)(14) 慣性モーメントとトルクに添え字 g がついているのは、重心を基準にしていることを表している。
式(12)の第1式を見ると、重心 xg の運動方程式は、質点の運動方程式と同じ形になっている。即ち、 xg の時間変化を知るだけであれば、剛体に働く外力の和さえ分かればよく、物体の形状を考慮する必要はない。キャッチボールや振り子を考える際に、物体の形状を考慮しなかったが、実際それでよかったわけである。
次に、式(12)の第2式についてである。簡単のため、慣性モーメント Ig がスカラー行列(=単位行列に実数をかけた行列)になる場合を考えると、その時間微分は、 ˙Ig=0 となる(以下の【11.1-注2】)。従って、式(12)の第2式は Ig˙ω=τg となる。これは、第1式 m¨xg=f と同じ形をしている。質量 m が大きくなるほど速度を変化させづらくなるのと同様に、 Ig は、大きいほど回転運動を変化させづらくなるような量(=回転の慣性を表す量)と見なせる。一方、トルク τg は、物体を回転させようとする「力」のようなものということになる。式(14)を見ると、トルクが最大になるのは、重心方向と外力が直交する時であることが分かる。例えば、ボールに力を加えて回転させる時、最も効率よく回転させることができるのは、球面に沿った方向に力を加える場合であることが直感的にわかる。実際この時、ちょうどトルクの大きさも最大になっている。逆に、ボールの重心に向かうような力がかかっている場合、トルクが 0 になり、ボールの回転には影響しない。
一方、 Ig がスカラー行列でない場合、式(12)の第2式を ˙ω=[] の形に変形すると、以下のようになる:(以下の【11.1-注2】を代入する) ˙ω=I1g(τgω×Igω)(15) これからも分かるように、 τg=0 であっても、右辺第2項が残るので、一般には ˙ω0 である。つまり、外力が働いていない場合であっても、回転軸が変化することになる。

【11.1-注2】慣性モーメント Ig の時間微分

慣性モーメント Ig の時間微分は、以下のようになる: ˙Ig=ω×IgIgω×(16)

導出

定義式(13)の微分を素直に計算すると以下のようになる:(見やすくするため xigxi と略記している)
˙Ig=ddtimi(|xi|2xixTi)‖ ‖ ‖ ‖ ‖ ‖ ‖ ‖ ‖ ‖ddt|xi|2=2xTi˙xi=2xTiω×xi=0ddtxixTi=˙xixTixi˙xTi=ω×xixTi+xixTiω×=imi(ω×xixTi+xixTiω×) 一方、式(16)の右辺も変形すれば同じ結果になる: ω×IgIgω×=ω×imi(|xi|2xixTi)=ω×imi(|xi|2xixTi)ω×=imi(ω×xixTi+xixTiω×)

11.2 剛体運動の計算

運動方程式(12)を用いれば、重心速度 ˙xg と角速度 ω の時間変化が計算できることになる。しかし、初期値をどのように設定するかなど、はっきりさせるべき点がある。この節では、それら実際の計算な議論を行う。特に、1階の正規形に変形すると式(21)のようになる。

剛体のモデリング

まず、剛体の形状を定義する必要がある。剛体の形状は変化しないので、適当な位置・向きに配置し、その時の各質点要素 mi,˜xi を定義すれば十分である。これを剛体のモデル位置と呼ぶことにする。このモデル位置での慣性モーメント ˜Ig と重心 ˜xg を計算しておく: ˜Ig=mi(˜xig2˜xTig˜xTig)˜xg=1mimi˜xi⎪ ⎪ ⎪⎪ ⎪ ⎪(17) モデル位置は、計算しやすいようにとればよい。例えば、 ˜Ig が対角行列になるようにとるといった具合である(以下の【11.2-注1】の式(19)のように対角行列にすることは、常に可能である)
剛体の位置・向きが変化した時、任意の質点要素の座標 xi は、重心 xg と回転行列 R を与えることによって以下のように決まる:(前章の10.1節) xi=xg+R(˜xi˜xg)(18) よって xg,R の時間変化を計算すれば、剛体の全ての質点要素の運動を計算できる、即ち、剛体の運動が計算できることになる。なお、任意の R での慣性モーメントについても、以下の【11.2-注1】の式(20)により計算できる。

【11.2-注1】慣性モーメントの性質

[1] どのような剛体であっても、適当な R をとることによって、即ち適当に回転させることによって、 I を以下のように対角化することができる: I=I1000I2000I3(19) 対角成分 I1,I2,I3 を主慣性モーメントという。これは、線形代数のスペクトル定理「任意の対称行列は回転行列で対角化できる」の帰結である。
[2] モデル位置での慣性モーメントが ˜I である剛体が、回転行列 R だけ回転したとする。回転後の慣性モーメント I は以下のようになる: I=R˜IRTI1=R˜I1RT}(20)
[3] 全ての質点要素が同一直線上に乗っている場合、主慣性モーメントのなかに0となるものが現れる。すると、 I に逆行列が存在しなくなり、運動方程式(12)は成立しなくなる。これは自然な結果である。というのも、全ての質点要素が同一直線上に乗っている場合、その直線の周りの回転角度が意味をなさなくなるためである。逆に、質点要素が、平面的あるいは立体的に分布している場合には、 I は常に逆行列を持つ。

解くべき運動方程式は、式(21)

上述の通り、剛体の運動を計算することは、重心位置 xg(t) および回転行列 R(t) の時間変化を計算することに他ならない。そのためには、運動方程式(12)を、1階の微分方程式(第1章の【1.3-注1】)にすると見通しがよい:(上2つが重心に関するもの、下2つが回転に関するもの) ddt⎢ ⎢ ⎢ ⎢xg˙xgRω⎥ ⎥ ⎥ ⎥=⎢ ⎢ ⎢ ⎢˙xgm1fω×RI1g(τgω×Igω)⎥ ⎥ ⎥ ⎥(21) (最後の成分は式(15)である。)時刻 t=0 における [] 部分の値を与えたうえで、微小な時間ステップ δt ずつ時刻を進めていけば、 xg(t) および R(t) を求めることができる。
ただし、質点要素に作用する外力が変化する場合は、 fτg は、数値計算のステップごとに与えられている必要がある。また、式(21)の第4成分では、各時刻での慣性モーメント Ig が必要になる。これは、式(20)により、 R とモデル座標での慣性モーメント ˜Ig から計算できる。
あるいは、 R の代わりに、第10章の10.3節で述べたオイラー角などの自由な座標 θ を用いることもできる。その場合、同章の【10.3-注2】や【10.3-注3】に示した ω˙θ の関係式を ˙θ=Aω と書くことにして、以下のようになる: ddt⎢ ⎢ ⎢ ⎢xg˙xgθω⎥ ⎥ ⎥ ⎥=⎢ ⎢ ⎢ ⎢˙xgm1fAωI1g(τgω×Igω)⎥ ⎥ ⎥ ⎥(22) θ の初期値は任意の値をとることができる。なお、 Ig の計算には R が必要だが、 θ から計算できる。

【例題】重力下ではトルクが0

例として、外力として重力のみが作用している場合を考える。力の総和 f およびトルク τg は、例外的に簡単になる:
f=imig=mgmimiτg=ixig×mig=0imixig=m(xgxg)=0 よって、運動方程式(12)の第1式より、重心 xg の運動方程式は m¨xg=mg となり、第1章の質点のキャッチボールの場合と同じになる。また、回転部分については、同第2式よりトルクが発生しないので、重力は回転には影響しないことも分かる。
そこでここでは、外力が働いていない場合に、剛体がどのように回転運動するのかについて、数値計算を行う。例えば右図は、立方体の各頂点に等しい質量のおもりを配置した場合である。以下の数値計算では、剛体の回転安定性について、いくつかの例を示す。