テニスラケットの定理をシミュレーションで確認する2022/12/11POSTB!はじめに「テニスラケットの定理」って知ってますか?3軸の慣性モーメントが異なる物体の回転運動を考えた際に、第二慣性主軸(慣性モーメントが大中小となるように並べたときの二番目の慣性主軸)周りの運動は不安定になるという定理のことです。effect」で検索すると無重力環境でT字形の物体(T形レンチか何か?)を回している動画が見つかると思います。見たことがない方はまず調べてみて下さい。動画の中で起きている不思議な運動は空力的な効果(乱流とか)とは関係なく発生するので、本記事では剛体の回転運動シミュレーションを行い、テニスラケットの定理(とそれにより発生する不思議な運動)を確認します。シミュレーションで回す物体は何でも良いのですが、折角なのでNASAの動画のものと同じ形状にします。適当に寸法を決めて以下のようにモデリングしました。シミュレーションにおいては、この物体を空気抵抗や重力といった外力が全く働かない状況で回転させます。限らず、物体の運動はダイナミクス(Dynamics)とキネマティクス(Kinematics)に分けて考えることができます。一般に、ダイナミクスは外力と速度の関係を表し、キネマティクスは速度と位置の関係を表します。両者の関係については文献[1]に以下のような図が掲載されています。<DynamicsとKinematicsの関係(文献[1]中の図6を参考に作成)>図中の上段が並進運動の場合、下段が回転運動の場合のダイナミクスとキネマティクスの関係を表しています。今回は回転運動のみのシミュレーションなので直接関係があるのは下段だけですね。<慣性座標系と物体座標系>また、式中の変数がどちらの座標系上の値であるかを上付きの添字で表します。キネマティクス説明の都合上キネマティクスの式を先に示します。回転運動の場合には姿勢表現の方法が複数あり、今回行うような物体がくるくる回るシミュレーションで特異点があると非常に厄介なので姿勢表現には四元数を用います。本記事中においては、四元数を表す変数は上にチルダを付けて表記します。\mathbb{R}^3$)は物体座標系上における角速度です。最終的にはこの微分方程式を解くことで物体の回転運動を計算します。すると剛体のダイナミクスは以下の式で表されます[4]。※(注:12月15日)以下の式に関して、トルクの座標系が怪しいので近いうちに修正します。\tag{4}$]の関係が成り立ち、この関係と式(3)より以下の式が得られます。$)に関する微分方程式ですから、これを数値的に解けばキネマティクスに必要な角速度が得られ、文献[4]においてもこの式を用いる方法が記載されています。この式を用いて実際にシミュレーションを行うと(特にオイラー法では)徐々に角運動量が変化してしまいます。$)は一定となるはずなのでこの結果は望ましくありません。そこで、式(2)と式(4)より角運動量を陽に扱います。$)を更新し、次に式(4)を用いて角速度を求めます。式(4)における角運動量は物体座標系上の値ですから、四元数を用いて物体座標系上の値に変換します。\tag{6}$]と書くことができるので[5]、式(4)より物体座標系上における角速度が以下のように求まります。$]これにより、角運動量保存則を成立させた状態で角速度を求めることができるようになりました。数値解法ここまでに示した微分方程式を数値解法で解いていきます。$)は常に一定となります。$)を計算しておくだけで済みます。仮に外力が加わる場合でも台形公式などを使って数値積分すれば計算できます。$)も簡単に計算できます。$)が含まれているので数値解法で解くために式(8)を式(1)に代入します。$)に関する微分方程式が得られたので、上式をルンゲ・クッタ法[5]で解くことで物体の回転運動を計算することができます。$]より、各計算ステップにおける更新式が以下のように得られます。$)(恒等四元数)とします。行列演算や四元数演算がやりやすいのでシミュレーション部分はRustで実装します。body_gyr.into())))}四元数の演算にはquaternion-coreクレートを使用しています。使いやすいと思うのでぜひ使って下さい(自分が作ったクレートなので使いやすいと感じるのは当然かもしれませんが...)。シンタックスハイライトがRustに対応していなくて見にくいのでGithubにも上げています。https://github.com/HamaguRe/tennis_racket_theoremシミュレーション結果のグラフを表示するPythonスクリプトとアニメーション表示用のProcessingコードもリンク先で見れます。$)としてシミュレーションを行いました。$)にしてしまうと全く遷移しなくなるので小さい値を設定しています。シミュレーション結果はCSVファイルに保存され、グラフにすると以下のようになります。$)軸周りの角速度を見ると何か規則的な動きをしていることが確認できると思います。全て物体座標系上の値)>0秒から10秒の範囲をアニメーション表示すると以下のようになります。NASAの動画と同じような挙動になっていますね!おわりに誰か宇宙に行く機会があったら本当にシミュレーション通りの挙動を示すか確認してみて下さい。Articleシンボリック計算を使わずに展開して多項式の係数を求める2022/10/15POSTAUTHORHamaguRe実験テーマは季節とともに移ろいゆく。