シミュレータのための剛体力学理論

Mitsuru Nishibe (mitsuman)

[email protected]

1999/03/10

1.はじめに

近年の

(中略)

だろう．

そこで現実の物体に近い性質を持つ剛体を扱ったシミュレーションの研究を行う．

具体的にはシミュレータは物体を多面体で定義し，全ての物体間の作用は衝突により処理することにより汎用性の高いシステムを目指した．

ここでは，シミュレーションに必要な剛体力学の式の導出及び解析と，姿勢制御に用いられるquaternionの解析を行う．参考文献と併せて読むと良いと思う．

2.理論

2.1.剛体力学による運動方程式

剛体を扱う場合，質点の移動にニュートンの運動方程式，回転にオイラーの運動方程式を用いるのが一般的である．移動に関しては扱い方は1通りしかないが[1]，回転は幾つかの表現方法がある．まずはニュートンの運動方程式から導かれる形から考える．

使う記号の意味は以下に準ずる．

2.1.1.並進運動

まず質点の移動について考えると，ニュートンの運動方程式より

(2.1)

であり，速度の変化量が得られる．

(2.2)

位置の変化量は速度の定義より明らかなので，これで質点の移動に関する運動方程式は十分である．

2.1.2.角速度による姿勢変化

姿勢の変化量は位のように単純には得られない．

まず，姿勢を表す行列はその慣性系の基底ベクトルを束ねたものであるのでベクトルを使うと，次のようになる．

(2.3)

が重心周りの角速度を表しているので，重心からの位置にある点の速度はであるから，これを各基底ベクトルにあてはめると，

(2.4)

表現が冗長なので行列を定義する．

(2.5)

この行列は外積に相当する事がわかる．

これを用いるとは次のように簡潔に表記できる．

(2.6)
2.1.3. 回転モーメント

角速度の変化量を求める前に角運動量について解析しよう．剛体を質点の集まりと考え，それぞれの中心からの位置を，質量をとすると，角運動量の定義は

(2.7)

である．これを微分しよう．

ここで，

(2.8)

として原点周りの外力の回転モーメント(トルク)をと定義する．

これより「角運動量の時間変化はトルクに等しい」ことが言える．

(2.9)

2.1.4.. 慣性テンソル

さてここで剛体が角速度を持っているとして角運動量を展開してみよう．

このように角速度以外の部分を２階のテンソルとしてまとめることができる．これを慣性テンソルという．また対角の部分を慣性モーメント，それ以外を慣性相乗モーメントと呼ぶ．

(2.10)

この慣性テンソルは姿勢の関数であるため姿勢が変更される度にこの膨大な計算が必要である．姿勢を用いた式に書き改めて検討しよう．

まず，とを利用して，

(2.11)

は姿勢の関数なので基準の姿勢の位置をとすれば

姿勢に依存しない定数が得られたので基準姿勢の慣性テンソルとして定義する．

(2.12)

すると任意の姿勢の慣性テンソルは

(2.13)

と表記することができる．を事前に求めておけば容易に慣性テンソルが求まることが分かる．なお，慣性テンソルは対称行列であるから任意の直交行列で対角化が可能である．即ち任意の正規直交行列つまり姿勢で対角化される．このことを主軸変換といい，新しい座標系を慣性主軸と呼ぶ．最終的に

(2.14)

を得る．

2.1.5.トルクから角速度の計算

トルクから角速度の変化を求めよう．角運動量の式の左辺を展開すると，

(2.15)

(2.13)より，

よって

(2.16)

角運動量の式(2.9)の右辺と合わせると

一見複雑になってしまったが，系全体を回転させれば簡単になる．

(2.17)

をワールド座標系と捉えず物体の座標系と考えて置きなおせば

(2.18)

展開するとオイラーの運動方程式になる．

(2.19)

さて，角速度の変化を求める式を得たが，もっとシンプルに求める方法がある．それは単にLを求めた後に

(2.20)

とするだけである．さらに系の変換をする必要も無い点でもオイラーの運動方程式より扱いが簡単であるので，角速度を求めるのにオイラーの運動方程式を使う必要は無い[2]．

2.1.6.剛体の状態ベクトル

並進運動と回転運動が揃ったので，分かりやすく剛体の状態ベクトルとその微分を定義しよう．

(2.21)

2.2.quaternion

2.2.1.quaternionによる姿勢管理

ここまでで導いた式により数値積分を繰り返すと，計算機の有効桁数の問題で誤差が蓄積されていく．それによりが正規直交行列でなくなり，物体が変形したことになる．これを防ぐため正規直交化を定期的に行う必要がある．そこで１つの実数とi,j,k３つの虚数単位からなるquaternionと呼ばれる4元数を導入する．このquaternionを正規化すると姿勢が表せる．Quaternionは行列のようにアフィン変換を全て表現可能でない分，姿勢管理に向いている[3]．

虚数部をまとめて１つのベクトルとしても表記できる．

(2.22)

絶対値，規格(Norm)，正規化，共役をそれぞれ，

(2.23)

(2.24)

(2.25)

(2.26)

と定める．を満たす逆quaternionは，

(2.27)

となる．正規化されている場合はであることが分かる．

虚数単位には以下の性質がある．

(2.28)

このことから幾つかの演算が導かれる．とすると，

(2.29)

(2.30)

ところで，

(2.31)

を定義し(は正規化されているとする)，ベクトルに対して下の演算を行う．

(2.32)

するとを周りにθ回転させたベクトルが得られる[4]．これは姿勢行列Rと同じ役目をする．つまり姿勢をquaternionで表すことが可能である．

との関係はを時間微分すると次のようになる．

(2.33)

以上の結果よりquaternionを使い剛体状態ベクトルを再定義しよう．

(2.34)

2.2.2.回転quaternionの行列表記

姿勢をquaternionで表すと様々な利点があるが，座標変換時には行列にしてスケーリング行列等とまとめると効率がよい．また姿勢行列が無いとIが計算できない．このことよりquaternionを使った回転を行列表記する必要がある．

(2.28)の性質を満たすquaternionの行列表記を示す．

(分かり易くするためとする)

(2.35)

変換前のベクトルと変換後のベクトルを定義する．

これを用いて(2.32)を計算しよう．

より，

は正規化されているのでであるから，

故に，

(2.36)

を得る．

2.3.衝突

多面体同士の接触を検出し，その時点の速度，角速度，質量，慣性テンソル，反発係数から撃力を与える．接触点の座標を，接触点での相対速度を，接触面の法線をとすると接触店における撃力は

(2.37)

である．

3.実装

3.1. 環境

Windows95/98上でWin32アプリケーションとして作成した．3D表示にはDirect3Dを使用した．C++で記述しVisualC++6.0でコンパイルした．

開発環境のハードウェアはK6-2/300MHz メモリ96MBである．

3.2. 剛体

理論で解析した剛体の状態ベクトルの微分方程式をオイラー法で実装した．

3.3.衝突検出の高速化

n個の物体の衝突チェック回数は，つまりである．多面体同士の衝突判定は軽い処理ではないのでチェック回数を減らす必要がある．まず考えうるのが物体同士の外接球による判定を行い，完全に外れている物体同士では衝突判定を行わないとする方法を試した．単純な方法ではあるが，非常に効果があり状態にもよるが判定回数を抑える働きがあった．それに加え凹面が無い物体に限定すれば物体の裏にあたる部分では衝突は起きないので，そのことを利用するとさらに半分近くの無駄な判定を回避できた．

違ったアプローチとしては毎回判定するのではなく一度いつ頃までは衝突しないかを予想し，それまでは衝突判定を行わないといったアルゴリズムもあるが，地面に接触して常に衝突が起きている場合を考えると，全くの無駄，むしろその計算の分だけ重くなりかねないので試してはいない．その代案として空間Aの中に物体がn個あるとして，仮に空間BとCに分割しそれぞれn/2個ずつの物体に分割し判定回数を半分に減らす事を思いつき試したが，物体数が少ないときは分割するコストだけ遅くなってしまった．しかし物体数が多くなる場合劇的とは言えないが速度の向上が確認できた．

4.結果

摩擦力を考慮していない以外はかなり自然にシミュレートできた．具体的には斜めに積み上げた板の崩壊，ドミノ倒し，コマ等の運動がシミュレートできた．

4.1.曲面への物体の落下

曲面は三角関数より作り100ポリゴンに分割したものを使用．

4.2.斜めに積み上げた板の崩壊

5.今後の課題

まず挙げられるのが摩擦の考慮である．これは一見簡単のように思えるが，面と面の接触の場合，正確な摩擦力の発生位置は接触面積を計算する必要がある．接触面積の計算は非常に重いため，考慮中である．

他には剛体同士の拘束力，例えばバネやダンパーの実装が挙げられる．また，マウスにより物体を直接動かすことが出来ればよりシミュレータの価値が上がるだろう．

画面表示では影を付ければ物体の位置が把握しやすくなるだろう．実際影がないと物体の位置がよくわからない．球などは姿勢を視覚的に見るために模様をつけることも必要だろう．

全般的にいえることは高速化である．多面体同士の衝突判定の処理が全体の処理の殆どを占めていることからその部分の高速化を重点的に研究していきたい．

参考文献

[1] David Baraff. An Introduction to Physically Based Modeling:Rigid body simulation I/II,1997

[2] 神谷 周,撃力ベースによる３次元剛体運動シミュレーション,東京大学精密機械工学科 卒業論文,1998

[3] 宇治社中,クウォータニオン,1998

---

[1] ここでは極座標，円柱座標といった直交座標以外の座標系は扱わないことにする

[2] 計算量はオイラー法が慣性テンソルを対角化していることを前提にしている分若干少ない．しかしが後々必要になることを考慮するべきである．

[3] 9個のパラメータで姿勢を表現する行列よりも4個のパラメータで表現するquaternionの方が冗長さが無く誤差が少なくなると経験的に知られている．またquaternionが球面線形補間等を容易にする事も姿勢を扱う点ではポイントが高い．

[4] 参考文献[3]を参照せよ