剛体衝突の動力学（速度と角速度の変化）

## これは何？

２つの剛体の衝突によって、速度や角速度がどのように変化するのかを計算するための理論と計算です。参考に、この記事に、ChatGPT o1 mini による計算結果があります（対話しながら導きました）。

箱庭ドローンシミュレータを作っている時に使っています。ドローン同士の衝突、ドローンと固定壁の衝突、ARにおけるドローンと手の接触、などに使います。以下関連記事です。

https://qiita.com/kenjihiranabe/items/0fd58475d5a0fee29072

https://qiita.com/kenjihiranabe/items/498418ceb8a80f8242d0

## 定式化

衝突による力積(impulse)を計算します。これによって、速度の変化と角速度の変化を計算することができます。力積は力と時間の積であり、この量 $J$ を時間を分割して、数回にわたって徐々に速度に加算することもできますし、一回の更新で、どん、と速度を更新することもできます。

$$
\int_0^T F(t)dt = J \quad \text{（$J$は力積 - $T$は接触時間）}
$$

接触の力学モデルを立てます。まず、基本となる変数を定義します。図は2剛体の一般的なモデルとよりシンプルに一方が圧倒的に大きな質量と慣性モーメントを持っている場合（衝突に動じない）です。ここでは、一般的に話を進め、シンプルなモデルでの式も最後に示します。

![collision.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/51055/ffc920d3-dd1e-722f-6167-5aebbe9fd183.png)


- $m_1, m_2$ : それぞれの剛体の質量
- $I_1, I_2$ : それぞれの慣性モーメント（衝突位置中心）
- $v_1, v_2$ : 衝突前の並進速度ベクトル
- $\omega_1, \omega_2$ : 衝突前の角速度ベクトル
- $v'_1, v'_2$ : 衝突後の並進速度ベクトル
- $\omega'_1, \omega'_2$ : 衝突後の角速度ベクトル
- $r_1, r_2$ : 重心から接触点までのベクトル
- $F$ : 接触位置での力ベクトル
- $n$ : 接触面の法線ベクトル
- $e$ : 反発係数

摩擦は無視して力を法線方向のみとし、接触位置での力（作用・反作用）から、

$$
F = f(t) n,   \quad \text{ $f(t)$ は時間変化するスカラー関数}
$$

を仮定して、以下の式を立てます。

(1). 接触位置の力による並進運動方程式：

$$
\begin{array}{rl}
m_1\dot{v_1} &= F \\\\
m_2\dot{v_2} &= -F
\end{array}
$$

(2). 接触位置の力による回転運動方程式：

$$
\begin{array}{ll}
I_1\dot{\omega_1} = r_1 \times F \\\\
I_2\dot{\omega_2} = r_2 \times (-F)
\end{array}
$$

(3). 接触点での相対速度：

$$
v_{rel} = (v_1 + \omega_1 \times r_1) - (v_2 + \omega_2 \times r_2)
$$

(4). 衝突による力積：

$$
\int_0^T Fdt = J \quad \text{（$J$は力積 - $T$は接触時間）}
$$

(5). 反発係数の定義より：

$$
v_{rel}' \cdot n = -e(v_{rel} \cdot n)
$$

ここで $v_{rel}$ は相対ベクトル、 $n$ は接触面の法線ベクトル、 $'$ （ダッシュ）は衝突後を表します。

（参考：この4式は、運動量保存則と角運動量保存則と同値です。上2式を足し合わせて $J$ を消せば運動量保存。下に2式を接触位置における慣性モーメントに変換して足し合わせれば角運動量保存）

これらの式を連立して解くことで、衝突後の運動を求めることができます。具体的な解き方としては、(1),(2) を時間積分します：

反発の式(5)から $J$ について立式すると、解ける形になります。それを使って、(1),(2)k衝突後の運動量変化と角運動量変化を計算します。

$$
\begin{array}{ll}
v_1' &= v_1 + J/m_1 \\\\
v_2' &= v_2 - J/m_2 \\\\
\omega_1' &= \omega_1 + I_1^{-1}(r_1 \times J) \\\\
\omega_2' &= \omega_2 - I_2^{-1}(r_2 \times J)
\end{array}
$$

これらの式を式(3)に代入し、

$$
\begin{array}{ll}
v_{rel}' &= (v_1' + \omega_1' \times r_1) - (v_2' + \omega_2' \times r_2) \\\\
 &= (v_1 + J/m_1 + \omega_1 \times r_1 + I_1^{-1}(r_1 \times J) \times r_1) - (v_2 - J/m_2 + \omega_2 \times r_2 - I_2^{-1}(r_2 \times J) \times r_2) \\\\
 &= v_{rel} + (1/m_1 + 1/m_2)J + I_1^{-1}(r_1 \times J) \times r_1 + I_2^{-1}(r_2 \times J) \times r_2
\end{array}
$$

これを、反発の式(5)に代入て解くことで、$J$ を求めることができます。
$j_n = J \cdot n$ すなわち $J = j_n n$ ( $J$ は法線方向のみ、摩擦なし)として、

$$
   j_n = -(1+e) (v_{rel} \cdot n)/ \left( \frac{1}{m_1} + \frac{1}{m_2} +
   (r_1 \times n) \cdot I_1^{-1} (r_1 \times n) + (r_2 \times n) \cdot I_2^{-1} (r_2 \times n) \right)
$$

これで、 $J$ が求まりました。これを元に、速度と角速度の変化量を求めることができます。

$m_2, I_2$ が非常に大きい場合、この式は以下に簡略化されます。

$$
   j_n = -(1+e) (v_{rel} \cdot n)/ \left( \frac{1}{m_1} +
   (r_1 \times n) \cdot I_1^{-1} (r_1 \times n) \right)
$$

## 移動体の速度や角速度への反映方法

この計算で求まるのは力積です。剛体1については、ここから、速度変化として $J/m$ を加算することで、衝突後の速度が求まります。角速度については、$I^{-1}(r \times J)$ を加算します。

上記は、一瞬で速度・加速度が変化するモデルです。衝突中の運動を追跡したい場合には、 $F(t), (0 \le t \le T)$ をデザインします。この関数は、剛体表面のバネ係数なども関係します。そして、$\int^T_0 F(t) dt = J$ となるように $F(t)$ と接触時間 $T$ を決め、その間、速度・加速度変化を与え続けます。

## まとめ

二つの剛体の衝突による速度・角速度変化について定式化し、解を求めました。

この物理学的考察は、claude と chat GPT o1 mini によってヒントを得て、手計算と整合させました。

- [Chat GPT o1 mini による考察](https://chatgpt.com/share/6789c6a0-5d4c-8002-bc1e-e54169e4fdf4)