HOME > natural science Laboratory > コンピュータ・シミュレーション講座 > 仮想物理実験室

【物理シミュレーションに挑戦!】古典力学
斜面を転がる剛体球運動の計算アルゴリズム16:力学的エネルギー保存則を満たした剛体球同士の衝突(回転あり)

文責:遠藤 理平 (2016年10月20日) カテゴリ:仮想物理実験室(247)計算物理学(126)

本項は古典力学の様々な系の物理現象を解析的に扱うのではなく、数値計算による物理シミュレーションを実行するために必要な計算アルゴリズムを示すことを目的とします。 様々な初期条件に対する物理シミュレーションを実現するために、最も汎用的な直交座標系を用います。
斜面を転がる剛体球運動の計算アルゴリズム1:摩擦力無しの場合
斜面を転がる剛体球運動の計算アルゴリズム2:摩擦力無しの場合の例(曲面上の運動)
斜面を転がる剛体球運動の計算アルゴリズム3:静止摩擦力のみの場合
斜面を転がる剛体球運動の計算アルゴリズム4:静止摩擦力のみの場合の例(曲面上の運動)
斜面を転がる剛体球運動の計算アルゴリズム5:滑りながら転がる剛体球
斜面を転がる剛体球運動の計算アルゴリズム6:滑りながら転がる剛体球(曲面上の運動)
斜面を転がる剛体球運動の計算アルゴリズム7:空気抵抗力と転がり摩擦抵抗力
斜面を転がる剛体球運動の計算アルゴリズム8:斜面との衝突する剛体球(回転なし)
斜面を転がる剛体球運動の計算アルゴリズム9:斜面との衝突する剛体球(回転あり)
斜面を転がる剛体球運動の計算アルゴリズム10:剛体球同士の衝突(回転なし)
斜面を転がる剛体球運動の計算アルゴリズム11:複数剛体球の同時衝突(回転なし)
斜面を転がる剛体球運動の計算アルゴリズム12:接触中の剛体球に衝突した場合(回転なし)
斜面を転がる剛体球運動の計算アルゴリズム13:剛体球同士の衝突(回転あり)
斜面を転がる剛体球運動の計算アルゴリズム14:固定された剛体球に衝突した場合(回転あり)
斜面を転がる剛体球運動の計算アルゴリズム15:力学的エネルギー保存則を満たした斜面との衝突(回転あり)
斜面を転がる剛体球運動の計算アルゴリズム16:力学的エネルギー保存則を満たした剛体球同士の衝突(回転あり)

斜面を転がる剛体球運動の計算アルゴリズム13:剛体球同士の衝突(回転あり)」で示した回転による推進力は、衝突力から見積もられる摩擦力として導出しました。この場合、回転運動を含めた力学的エネルギーを満たす条件を全く考慮していませんでした。本項では、回転する剛体球の衝突前後の力学的エネルギーが保存する条件から推進力を計算します。

運動量保存則から、衝突前後の重心運動量と回転運動量の変化量は力積(力×時間)で与えられます。そのため、2つの剛体球の衝突後の重心速度と角加速度は、衝突力を \mathbf{f}_{\rm col} 、回転による推進力を \mathbf{f}_{\rm d} と表した場合、

\mathbf{v'}_1 = \mathbf{v}_1+\frac{\Delta t}{m_1}\, (\mathbf{f}_{\rm col12} + \mathbf{f}_{\rm d12})
\mathbf{v'}_2 = \mathbf{v}_2+\frac{\Delta t}{m_2}\, (\mathbf{f}_{\rm col21} + \mathbf{f}_{\rm d21})
\mathbf{\omega'}_1=\mathbf{\omega}_1+\frac{a_1\Delta t}{I_1} (\mathbf{f}_{\rm d12} \times\hat{\mathbf{n}}_{12})
\mathbf{\omega'}_2=\mathbf{\omega}_2+\frac{a_2\Delta t}{I_2} (\mathbf{f}_{\rm d21} \times\hat{\mathbf{n}}_{21})

となります。ただし、 \mathbf{f}_{\rm col12} は1番目の剛体球に加わる衝突力を表し、「斜面を転がる剛体球運動の計算アルゴリズム10:剛体球同士の衝突(回転なし)」で示したとおり

\mathbf{f}_{\rm col12}  =\frac{1}{\Delta t}\left[ \frac{(1+e)m_1m_2}{m_1+m_2}\,(\mathbf{v}_1- \mathbf{v}_2)\cdot\hat{ \mathbf{n} }_{12} \right]\hat{ \mathbf{n}}_{12}

です。 また、\mathbf{f}_{\rm col12}\mathbf{f}_{\rm col21}は1番目と2番目の剛体球に加わる角運動量と並進運動量の交換を媒介する力(以後、運動量交換力と呼ぶことにします)です。 なお、\mathbf{f}_{\rm col12}\mathbf{f}_{\rm col21}は作用・反作用の法則から \mathbf{f}_{d21} = -\mathbf{f}_{d12} です。この\mathbf{f}_{\rm col12}をエネルギー保存則から導き出します。

 まず、2つの剛体球が衝突時に角運動量と並進運動量の交換力が生じる条件について考えます。交換力は2つの剛体球の回転を含めた衝突面における接点の相対速度が一致している場合には交換力は生じません。そこで2つの剛体球の接点の相対速度を並進運動と回転運動によるから考えます。接点の相対速度は必ず衝突面と平行成分しか持たないことを考慮すると並進運動に対する接点の相対速度は

並進運動に対する接点の相対速度 =\mathbf v_{1 \parallel } - \mathbf v_{2 \parallel } = (\mathbf v_{1} - (\mathbf v_{1} \cdot\hat{\mathbf{n}}_{12})\hat{\mathbf{n}}_{12} ) -  (\mathbf v_{2} -( \mathbf v_{2} \cdot\hat{\mathbf{n}}_{21})\hat{\mathbf{n}}_{21} )
= (\mathbf v_{1} -\mathbf v_{2} ) - \left\{(\mathbf v_{1} -\mathbf v_{2} )\cdot\hat{\mathbf{n}}_{12} \right\} \hat{\mathbf{n}}_{12}

となります。ただし、 \hat{\mathbf{n}}_{21}=-\hat{\mathbf{n}}_{12}} を考慮しています。同様に回転運動に対する接点の相対速度は

回転運動に対する接点の相対速度 =a_1(\mathbf{\omega }_1\times \hat {\mathbf n}_{12})-a_2(\mathbf{\omega }_2\times \hat {\mathbf n}_{21})
=(a_1\mathbf{\omega }_1+a_2\mathbf{\omega }_2) \times \hat {\mathbf n}_{12}

となります。よって、並進運動と回転運動を考慮した接点の相対速度ベクトルは上記の2つの相対速度ベクトルは

\mathbf D\equiv \left[ (\mathbf v_{1} -\mathbf v_{2} ) - \left\{(\mathbf v_{1} -\mathbf v_{2} )\cdot\hat{\mathbf{n}}_{12} \right\} \hat{\mathbf{n}}_{12} \right]-\left[ (a_1\mathbf{\omega }_1+a_2\mathbf{\omega }_2) \times \hat {\mathbf n}_{12} \right]

と定義することができます。 \mathbf {D} が有限の大きさを保つ場合、衝突時に並進運動の衝突力に加え、角運動量と並進運動量の交換が発生することになり、運動量交換力の方向は \mathbf {D} と同じ方向(衝突面の節ベクトル)になります。つまり、運動量交換力\hat{\mathbf{n}}_{21}=-\hat{\mathbf{n}}_{12}}は衝突面の接方向ベクトル

\hat{\mathbf{t}} \equiv \mathbf {D}/|\mathbf {D}|

を用いて

\mathbf{f}_{d12}=f_d\, \hat{\mathbf{t}}

と表されます。あとはこの運動量交換力の大きさf_dが未知の量です。f_dを力学的エネルギー保存則から導きましょう。

衝突後の力学的エネルギーは

E'=\frac{1}{2}m_1\mathbf{v'}_1^2+\frac{1}{2}m_2\mathbf{v'}_2^2+\frac{1}{2}I_1\mathbf{\omega'}_1^2+\frac{1}{2}I_2\mathbf{\omega'}_2^2

と表されます。衝突前後の力学的エネルギーの差を弾性衝突(反発係数e=1)を仮定して注意深く計算すると次のとおりになります。

E'-E =\frac{1}{2} \left[ \frac{1}{m_1}+\frac{1}{m_2} + \frac{a_1^2}{I_1} +\frac{a_2^2}{I_2}\right] \Delta t^2f_d^2
+ \left[(\mathbf{v}_1-\mathbf{v}_2)\cdot \hat{\mathbf{t}} +(a_1 \mathbf{\omega}_1 + a_2 \mathbf{\omega}_2)\cdot( \hat{\mathbf{t}} \times \hat{\mathbf{n}}_{12} )\right] \Delta t f_d

左辺を0として、推進力f_dに関する2次方程式が得られました。一つ目の解はf_d=0は自明な解、もう一つは

f_d \Delta t= \frac{2}{ \frac{1}{m_1}+\frac{1}{m_2} + \frac{a_1^2}{I_1} +\frac{a_2^2}{I_2} } \left[ (a_1 \mathbf{\omega}_1 + a_2 \mathbf{\omega}_2)\cdot( \hat{\mathbf{n}}_{12}\times\hat{\mathbf{t}} )  -(\mathbf{v}_1-\mathbf{v}_2)\cdot \hat{\mathbf{t}}\right]

です。この推進力を与えることで、衝突前後の剛体球の回転を計算することができます。

回転する剛体球同士の衝突シミュレーション(力学的エネルギーが保存)


物理シミュレーションについては「HTML5による物理シミュレーション」を参照ください。 数式の表示は「Tex表記によるHTML文書への式の埋め込み」をご覧ください。



▲このページのトップNPO法人 natural science トップ

関連記事

仮想物理実験室







計算物理学







▲このページのトップNPO法人 natural science トップ