剛体の衝突について考えてみます(2回目)。
本を片手に書いているので、おかしな所がありましたらぜひ教えてください。
一部説明を省略したりするので、前回の記事を呼んでおくことをお勧めします。
前回では、剛体の回転について全く考えていませんでした。
球体などは摩擦がなければ回転しないので前回の式がそのまま使えますが、
箱など(というか、ほとんどの形状で)は力を与える場所によっては回転が生じます。
剛体の回転運動に関する公式がこちらです。
τ=Iα
よく分からない記号が並んでいますが、F=maの回転版とも言えるもので、意味は大体同じです。
τはトルクといい、ひねる力を表します。
I は慣性モーメントといい、回転に対する抵抗の強さを表します。
αは角加速度といい、回転の加速度を表します。
ひねる力は、回転に対する抵抗と角加速度の積に等しい、ということです。
剛体に力を与えた時も、ひねる力は発生します。
ですが、与えた力=ひねる力ではありません。
ひねる力は、力を与えた場所によって変化します。
重心に向かって力を与えても、ひねる力は発生しません。
力を与えた位置が重心から離れているほど、力の向きが重心への方向と違うほど、
ひねる力は強くなります。
与えた力と、ひねる力の関係は、次の式で表せます。
τ=r×F
×はベクトルの外積で、r は重心から見て、力を与えた位置です。
外積は2ベクトルのなす平行四辺形の面積に比例する
さて、これらを元に回転も考慮した力積を求めてみましょう。
まずは前回と同じく反発係数の式です。
e = -n・(v'2 - v'1) / n・(v2 - v1)
おや、前回と少し式が違いますね。
衝突の反発においては、法線以外の方向への速度は関係ないので、
法線との内積を取って、その方向の力だけを取り出しています。
衝突で発生した力積をJとすると、発生する力Fは
F = nJ
です。
Fによる速度変化を考えると、衝突後の各剛体の速度は次のように表せます。
vL'1 = vL1 + F/m1
vA'1 = vA1 + r1×F/I1
vL'2 = vL2 - F/m2
vA'2 = vA2 - r2×F/I2
ここで vL は並進速度、vA は角速度を表します。
また、詳細は省略しますが、ある相対位置 r における各剛体の衝突前・衝突後の速度 v は、
v1 = vL1 + vA1×r1
v2 = vL2 + vA2×r2
v'1 = vL'1 + vA'1×r1
v'2 = vL'2 + vA'2×r2
と表せます。
後はひたすら式をこねくり回していけばJが求まりそうですね。
一度に計算すると非常にややこしくなるので、分割して計算します。
まずは v2 - v1、v'2 - v'1 を整理します。
v2 - v1 = (vL2 + vA2×r2) - (vL1 + vA1×r1)
= (vL2 - vL1) + (vA2×r2 - vA1×r1)
v'2 - v'1 = (vL'2 + vA'2×r2) - (vL'1 + vA'1×r1)
= (vL'2 - vL'1) + (vA'2×r2 - vA'1×r1)
= ([vL2 - F/m2] - [vL1 + F/m1]) + ([vA2 - r2×F/I2]×r2 - [vA1 + r1×F/I1]×r1)
= ([vL2 - vL1] - [F/m2 + F/m1]) + ([vA2×r2 - vA1×r1] - [{r2×F/I2}×r2 + {r1×F/I1}×r1])
= (vL2 - vL1) + (vA2×r2 - vA1×r1) - ([F/m2 + F/m1] + [{r2×F/I2}×r2 + {r1×F/I1}×r1])
= (v2 - v1) - ([nJ/m2 + nJ/m1] + [{r2×nJ/I2}×r2 + {r1×nJ/I1}×r1])
= (v2 - v1) - (nJ * [1/m1 + 1/m2] + J[{r2×n/I2}×r2 + {r1×n/I1}×r1])
ではこの式を反発係数の式に代入します。
e = -n・(v'2 - v'1) / n・(v2 - v1)
e = -n・( (v2 - v1) - (nJ * [1/m1 + 1/m2] + J[{r2×n/I2}×r2 + {r1×n/I1}×r1]) ) / n・(v2 - v1)
e * n・(v2 - v1) = -n・( (v2 - v1) - (nJ * [1/m1 + 1/m2] + J[{r2×n/I2}×r2 + {r1×n/I1}×r1]) )
e * n・(v2 - v1) = -n・(v2 - v1) + n・(nJ * [1/m1 + 1/m2] + J[{r2×n/I2}×r2 + {r1×n/I1}×r1])
e * n・(v2 - v1) + n・(v2 - v1) = n・(nJ * [1/m1 + 1/m2] + J[{r2×n/I2}×r2 + {r1×n/I1}×r1])
n・(v2 - v1) * (e + 1) = n・(nJ * [1/m1 + 1/m2] + J[{r2×n/I2}×r2 + {r1×n/I1}×r1])
n・(v2 - v1) * (e + 1) = n・(nJ * [1/m1 + 1/m2]) + n・J([r2×n/I2]×r2 + [r1×n/I1]×r1)
n・(v2 - v1) * (e + 1) = J * (1/m1 + 1/m2) + nJ・([r2×n/I2]×r2 + [r1×n/I1]×r1)
n・(v2 - v1) * (e + 1) = J * (1/m1 + 1/m2 + n・[{r2×n/I2}×r2 + {r1×n/I1}×r1])
J = n・(v2 - v1) * (e + 1) / (1/m1 + 1/m2 + n・[{r2×n/I2}×r2 + {r1×n/I1}×r1])
Jが求まりました。
n・(v2 - v1)を法線方向の相対速度 vR とすれば、
J = vR * (e + 1) / (1/m1 + 1/m2 + n・[{r2×n/I2}×r2 + {r1×n/I1}×r1])
と書くこともできます。一般的にはこちらの式が使われています。
・・・これで角度を考慮した衝突の計算ができるようになりました。
いざコードに書き起こすとなると、なかなか大変ですが、
やる気のある方はぜひ挑戦してみてください!
P.S.しばらく出かけていますので、その間コメントの返信などできなくなります。
≪ 氷追加 | HOME | 木材追加&剛体の衝突 ≫
≪ 氷追加 | HOME | 木材追加&剛体の衝突 ≫
Author:saharan