執筆者：廣瀬
映像業界を経てモノリスソフトへ入社。 以来、テクニカルアーティストとして主にエフェクト関連の業務を担当。 好きな食べものはソフトクリーム。
TECH BLOG
宇宙でスパナをぶん回そう
はじめに
こんにちは。モノリスソフト テクニカルアーティストの廣瀬です。
今回は、宇宙空間でよく見るスパナの回転をHoudiniでシミュレーションしてみたいと思います。
Houdiniは最後に実装の解説で少し出てくるだけですので、Houdiniに詳しくなくても問題ありません。
この運動を再現するには、歳差運動を考慮したシミュレーションをする必要があります。
歳差運動は、オイラーの回転運動方程式を解くことでシミュレーションすることができます。
ですので、まずはオイラーの回転運動方程式について解説し、次にオイラーの回転運動方程式を解くために必要な慣性テンソルなどの算出方法を解説します。
最後に、シミュレーション方法と、Houdiniでの実装の解説します。
（本記事では行列は列優先で記述します）
サンプルファイル（monolithtech_inertia.zip 45.8 KB）はこちらです。
オイラーの回転運動方程式
歳差運動のシミュレーションをするためには、オイラーの回転運動方程式を解く必要があります。
剛体のローカル座標系におけるオイラーの回転運動方程式は以下です。
ここで N はトルク、 ω は剛体の角速度ベクトル、L は角運動量です。
\(N = \dfrac{dL}{dt} + \omega \times L \)
角運動量は以下の式で表すことができます。
ここで、I は慣性モーメントで3×3の対称行列であり、2階のテンソルです。そのことから慣性テンソルなどと呼ばれたりもします。
\(L = I\omega \)
よって、角加速度ベクトルは以下の式で求めることができます。
\(\dfrac{d\omega}{dt} =I^{-1} (N - \omega \times I\omega) \)
ですので、歳差運動のシミュレーションをするためには、まず、シミュレーション対象のメッシュの慣性テンソルを求める必要があります。
また、剛体は重心を中心に回転運動をするため、メッシュの重心を求める必要もあります。
それらの計算をする際にメッシュの体積も求める必要があるため、まずは体積の求め方から解説していきます。
体積
まず、メッシュの体積の求め方です。
メッシュの体積を求めるためには、初めに全てのフェースを三角化します。
三角化したフェースのうちの１つのフェースに注目して考えます。
フェースと原点とを結んだ四面体を構築します。メッシュ全体の体積はこの四面体の合計と考えることができます。
\( volume_{mesh} = \sum_{face\in{mesh}}{} volume_{face} \)
NOTE:
下の画像のようなフェースを選択した場合、メッシュの表面から飛び出した四面体が構築されてしまいますが、
このような場合、飛び出た分だけ別のフェースの四面体が負の体積として計算されるため、メッシュ全体の体積として辻褄が合うようになっています。
フェースの各頂点の位置から3×3行列を構築して行列式（スカラー三重積）を計算します。
行列式の値は各ベクトルが作る立体の体積です。ですので、行列式を6で割った値が四面体の体積となります。
\( \begin{eqnarray} A_{face} &=& \begin{bmatrix} r_{1xface} & r_{2xface} & r_{3xface} \\ r_{1yface} & r_{2yface} & r_{3yface} \\ r_{1zface} & r_{2zface} & r_{3zface} \end{bmatrix} \\ volume_{face} &=& \dfrac{\det A_{face}}{6} \end{eqnarray} \)
よって、メッシュ全体の体積は以下で求まります。
\( volume_{mesh} = \sum_{{face}\in{mesh}} \dfrac{\det A_{face}}{6} \)
重心
続いて、重心の求め方です。
メッシュ全体の重心は体積と同様に、四面体の重心を合計することで計算することができます。
四面体の重心は以下の式で求まります。
\( \vec{center\:of\:mass}_{face} = \dfrac{\vec{r}_{1face} + \vec{r}_{2face} + \vec{r}_{3face} + O}{4} \)
よって、メッシュ全体の重心は四面体の体積（質量）を加味した以下の式で求まります。
\( \vec{center\:of\:mass}_{mesh} = \dfrac{1}{{volume}_{mesh}} \sum_{{face}\in{mesh}} \dfrac{\det A_{face}}{6} \dfrac{\vec{r}_{1face} + \vec{r}_{2face} + \vec{r}_{3face} + O}{4} \)
慣性テンソル
慣性テンソルの求め方です。
慣性テンソルも体積・重心と同様に、四面体の慣性テンソルを合計することで計算することができます。
\( I_{mesh} = \sum_{{face}\in{mesh}} I_{face} \)
慣性テンソルの式に関しては、天下り的に記載しておきます。詳しく知りたい方は以下の文献をご覧ください。
参考
慣性テンソルは3×3対称行列です。
\( I_{face} = \begin{bmatrix} a_{face} & -c^{\prime}_{face} & -b^{\prime}_{face} \\ -c^{\prime}_{face} & b_{face} & -a^{\prime}_{face} \\ -b^{\prime}_{face} & -a^{\prime}_{face} & c_{face} \end{bmatrix} \)
行列の各成分は以下の式で求まります。ρ は剛体の密度です。
\( \begin{eqnarray} \vec{v}_1 &=& \vec{r}_{1face} - \vec{center\:of\:mass}_{mesh} \\ \vec{v}_2 &=& \vec{r}_{2face} - \vec{center\:of\:mass}_{mesh} \\ \vec{v}_3 &=& \vec{r}_{3face} - \vec{center\:of\:mass}_{mesh} \end{eqnarray} \)
\( \begin{eqnarray} a_{face} &=& \rho \det A_{face} \dfrac{1}{60}( \\ && \quad v_{1y}^2 + v_{1y}v_{2y} + v_{1y}v_{3y} + v_{2y}^2 + v_{2y}v_{3y} + v_{3y}^2 \\ && + v_{1z}^2 + v_{1z}v_{2z} + v_{1z}v_{3z} + v_{2z}^2 + v_{2z}v_{3z} + v_{3z}^2) \\ \\ b_{face} &=& \rho \det A_{face} \dfrac{1}{60}( \\ && \quad v_{1z}^2 + v_{1z}v_{2z} + v_{1z}v_{3z} + v_{2z}^2 + v_{2z}v_{3z} + v_{3z}^2 \\ && + v_{1x}^2 + v_{1x}v_{2x} + v_{1x}v_{3x} + v_{2x}^2 + v_{2x}v_{3x} + v_{3x}^2) \\ \\ c_{face} &=& \rho \det A_{face} \dfrac{1}{60}( \\ && \quad v_{1x}^2 + v_{1x}v_{2x} + v_{1x}v_{3x} + v_{2x}^2 + v_{2x}v_{3x} + v_{3x}^2 \\ && + v_{1y}^2 + v_{1y}v_{2y} + v_{1y}v_{3y} + v_{2y}^2 + v_{2y}v_{3y} + v_{3y}^2) \\ \\ a^{\prime}_{face} &=& \rho \det A_{face} \dfrac{1}{120}( \\ && \quad 2v_{1y}v_{1z} + v_{1y}v_{2z} + v_{1y}v_{3z} \\ && + v_{2y}v_{1z} + 2v_{2y}v_{2z} + v_{2y}v_{3z} \\ && + v_{3y}v_{1z} + v_{3y}v_{2z} + 2v_{3y}v_{3z}) \\ \\ b^{\prime}_{face} &=& \rho \det A_{face} \dfrac{1}{120}( \\ && \quad 2v_{1z}v_{1x} + v_{1z}v_{2x} + v_{1z}v_{3x} \\ && + v_{2z}v_{1x} + 2v_{2z}v_{2x} + v_{2z}v_{3x} \\ && + v_{3z}v_{1x} + v_{3z}v_{2x} + 2v_{3z}v_{3x}) \\ \\ c^{\prime}_{face} &=& \rho \det A_{face} \dfrac{1}{120}( \\ && \quad 2v_{1x}v_{1y} + v_{1x}v_{2y} + v_{1x}v_{3y} \\ && + v_{2x}v_{1y} + 2v_{2x}v_{2y} + v_{2x}v_{3y} \\ && + v_{3x}v_{1y} + v_{3x}v_{2y} + 2v_{3x}v_{3y}) \\ \end{eqnarray} \)
シミュレーション
オイラーの回転運動方程式と、各変数の値についての解説が済みましたので、シミュレーション方法について解説していきます。
シミュレーションは簡単で、
-
初期値のトルク (N) を与える。
-
剛体のローカル座標系でオイラーの回転運動方程式を解いて角加速度を求める。
-
角加速度を元に、角速度 (ω) と剛体の回転 (姿勢) を求める。
-
3. で計算した角速度 (ω) と剛体の回転 (姿勢) を新たな値として設定 (積分) する。
-
トルク (N) をリセットする。
-
2. に戻る。
と、フレーム毎に 2 ~ 6 の工程を繰り返すだけです。
今回のシミュレーションでは 2. と 3. の計算をする際、オイラー法では精度が不足する可能性がある為、4次のルンゲ・クッタ法を使用しました。
オイラー法とルンゲ・クッタ法に関しては、今回詳しい説明は省かせていただきます。
Houdiniでの実装
最後に、Houdiniでの実装について解説していきます。
- まずは、スパナをモデリングします。
- 全てのフェースを三角化します。
- 各フェース毎の四面体の行列式・体積・重心を計算します。
行列式の値は、今後慣性テンソルを計算するときにも使用するため、一時的に値を保存しておきます。
- 体積・重心それぞれの総和を取り、メッシュ全体の体積・重心を計算します。
- 四面体の慣性テンソルを計算します。
HoudiniのRBDソルバで用いられている慣性テンソルアトリビュートの値は質量を除外した値のため、それにならってここでは密度の乗算は行っていません。
慣性テンソルへの質量の反映はシミュレーション時に、慣性テンソルアトリビュートの値と、質量アトリビュートの値を乗算することで行います。
- メッシュ全体の慣性テンソルを計算します。
上記で解説したHoudiniと同じ慣性テンソルの値にするため、体積で除算しています。
- 体積と密度と乗算して質量を計算します。
- メッシュ全体を1つの剛体として扱うためにパックを行います。
また、剛体の中心が重心となるようにピポットを移動します。
- 初期トルク (N) を作成します。
- Solver SOP で毎フレーム角速度と剛体の回転 (姿勢) を蓄積 (積分) しながら、4次のルンゲ・クッタ法でオイラーの回転運動方程式を解いていきます。
function vector dw(const vector torque, w; const matrix3 I) { // Euler's equations // N = (dL/dt)space = (dL/dt)body + cross(w, L) return invert(I) * (torque - cross(w, I * w)); } vector t = v@torque; vector w = @w; vector4 q = @orient; matrix3 I = 3@itensor * @mass; w = qrotate(qinvert(q), w); vector kw1 = w; vector4 kq1 = q; vector kdw1 = dw(qrotate(qinvert(kq1), t), kw1, I); vector kw2 = w + kdw1 * @TimeInc * 0.5; vector4 kq2 = qmultiply(q, quaternion(kw1 * @TimeInc * 0.5)); vector kdw2 = dw(qrotate(qinvert(kq2), t), kw2, I); vector kw3 = w + kdw2 * @TimeInc * 0.5; vector4 kq3 = qmultiply(q, quaternion(kw2 * @TimeInc * 0.5)); vector kdw3 = dw(qrotate(qinvert(kq3), t), kw3, I); vector kw4 = w + kdw3 * @TimeInc; vector4 kq4 = qmultiply(q, quaternion(kw3 * @TimeInc)); vector kdw4 = dw(qrotate(qinvert(kq4), t), kw4, I); vector kdw = (kdw1 + kdw2 * 2 + kdw3 * 2 + kdw4) * @TimeInc / 6; vector kw = (kw1 + kw2 * 2 + kw3 * 2 + kw4) * @TimeInc / 6; v@torque = set(0, 0, 0); @orient = qmultiply(q, quaternion(kw)); @w = qrotate(@orient, w + kdw);
- おわり。
まとめ
今回は歳差運動のシミュレーションをしてみました。
入力するメッシュを違うものにしてシミュレーションしてみても面白いと思います。また、今回実装はHoudiniで行いましたが、Mayaなど別のソフトでも実装できます。興味がある方は色々なソフトで実装してみて、宇宙にいる感覚を味わってみてください。