概要
土木の分野において、粒子法の研究が盛んにおこなわれています。例えば、離散要素法（DEM）、Smoothed Particle Hydrodynamics（SPH） 、Material Point Method（MPM）などがあります。
（球形）DEMは、昔から地盤力学の分野で使われている手法で、地盤解析・土砂崩れ・落石解析・流体と組み合わせれば地盤の越流現象など広い範囲で適応されています。
前回の記事で、球形DEM の紹介をしました。
実際の砂粒や岩石は、歪な形をしており球体 DEM で表現することは難しいことが知られています。球体 DEM においては、回転摩擦を加えることで、歪な形状の物体を表現できます。しかし、回転摩擦にも限界があり、任意の形状を表現したシミュレーションを行う方が良いと考えられます。
今回の記事では、任意の形状（粒子を繋げて表現できる範囲）で解析でき、さらに、物性値は反発係数と摩擦係数のみで計算ができる力積型個別要素法（iDEM）の紹介をします。
以下の書籍を参考にしています。
接触条件
接触・非接触の条件について述べる。
2つの物体（剛体）A,B を考える。物体は粒子で構成されているとして、物体 A に属する粒子 $i$ と物体 B に属する粒子 $j$ が衝突することを考える。接地・貫入条件は、粒子の位置を $x$ とすると、
\begin{align}
C_p(x_{iA},x_{jB},t) = (\vec{x}_{iA}(t)-\vec{x}_{jB}(t))\cdot \vec{n}_{jB} \leq 0
\end{align}
と書くことができる。 $n_{jB} $ は、物体 B における粒子 $j$ の法線ベクトルである。貫入条件は、$C_p(x_{iA},x_{jB},t) < 0$ である。
重心位置ベクトル $X$、相対位置ベクトル $r$ （重心から粒子までのベクトル $r=x-X$）とすれば、
\begin{align}
C_p(x_{iA},x_{jB},t) = \{(\vec{X}_A+\vec{r}_{iA})-(\vec{X}_B+\vec{r}_{jB})\}\cdot \vec{n}_{jB} \leq 0
\end{align}
と書くことができる。同様に速度に関しても考えることができ、
\begin{align}
C_v(x_{iA},x_{jB},t) & =\frac{\mathrm{d} C_p}{\mathrm{d} t} = \{\dot{\vec{x}}_{iA}(t)-\dot{\vec{x}}_{jB}(t)\}\cdot \vec{n}_{jB} \\
&=\{ (\dot{\vec{X}}_A+\vec{\omega}_A \times \vec{r}_{iA})-(\dot{X}_B+\vec{\omega}_B \times \vec{r}_{jB})\}\cdot \vec{n}_{jB} \leq 0
\end{align}
となる。$\omega$ は角運動量である。法線ベクトルは、接触の際は変化しないと仮定した。貫入条件は、$C_v(x_{iA},x_{jB},t) < 0$ である。
位置に関する接地・貫入条件 $C_p(x_{iA},x_{jB},t) \leq 0$ を満たしているが、 速度に関する接地・貫入条件を満たしていない、つまり、$C_v(x_{iA},x_{jB},t) > 0$ のとき、以前の計算ステップで衝突し接触力が与えられ、お互いの物体が離れる過程だと考えることができる。速度に関する接地条件は、以上の条件のとき、接触力を与えないために必要な条件である。
まとめると
\begin{align}
\left\{\begin{array}{ll}
C_p(x_{iA},x_{jB},t) >0 & \mbox{非接触}\\
C_p(x_{iA},x_{jB},t) \leq 0 & \mbox{接地・貫入状態}
\Longrightarrow
\left\{\begin{array}{ll}
C_v(x_{iA},x_{jB},t) &>0 & \mbox{離反過程・非接触}\\
C_v(x_{iA},x_{jB},t) &\leq 0 & \mbox{接触}
\end{array}\right.
\end{array}\right.
\end{align}
となる。
力積
ニュートンの運動方程式と角運動量の時間微分は、
\begin{align}
m\frac{\mathrm{d} \vec{V}}{\mathrm{d} t} &= \vec{F} \\
\frac{\mathrm{d} (I\vec{\omega})}{\mathrm{d} t} &= \vec{N} = \vec{r} \times \vec{F}
\end{align}
と書ける。$I$ は慣性モーメントであり、$F,N$ は力とモーメントである。$V$ は、並進方向の速度である。
微小時間 $\Delta t$ に働く力が力積なので、力積を$J$ とすると、
\begin{align}
\vec{J} = \vec{F} \Delta t
\end{align}
と書ける。
力積を使い速度の増分を表すと、ニュートンの運動方程式から、
\begin{align}
& m\frac{\Delta \vec{V}}{\Delta t} = \vec{F} \\
\therefore \ \ & \Delta \vec{V} = \frac{\vec{F}}{m} \Delta t = \frac{\vec{J} }{m}\\
\end{align}
並進方向の速度の増分が得られる。
角運動量の時間微分から、
\begin{align}
\frac{\mathrm{d} (I\vec{\omega})}{\mathrm{d} t} &= I\frac{\mathrm{d} \vec{\omega}}{\mathrm{d} t}+\frac{\mathrm{d} I}{\mathrm{d} t}\vec{\omega} = \vec{r} \times \vec{F} \\
\frac{\Delta \vec{\omega}}{\Delta t} &= I^{-1} (\vec{r} \times \vec{F}) - I^{-1}\frac{\mathrm{d} I}{\mathrm{d} t}\vec{\omega}
\\
\therefore \ \ \Delta \vec{\omega} &= I^{-1}(\vec{r} \times \vec{F} \Delta t) = I^{-1}(\vec{r} \times \vec{J})
\end{align}
回転方向の速度の増分が得られる。衝突時に力積を計算する際、慣性モーメントは変化しない（瞬間的に起こる）として$\mathrm{d} I/ \mathrm{d} t =0$ とした。力積が求まれば、位置や速度を更新しシミュレーションを行うことができる。
2体衝突問題の場合、接地条件
\begin{align}
C_v(x_{iA},x_{jB},t)=0
\end{align}
と反発係数 $e$ の関係
\begin{align}
v(t+\Delta t) = - e v(t)
\end{align}
から、力積 $J$ が具体的な形で求めることができる。多体衝突問題の場合は、具体的な形で力積が求まらないので、反復法を用いて力積を計算する。
エネルギー保存型力積法
以下では、エネルギー保存型力積法は以下の4点を前提とする。
- 反発係数は、Stronge の仮定
\begin{align}
e^2 = \frac{W_{rel}^n}{W_{comp}^n}
\end{align}
から定義される。$W_{comp}^n$ は圧縮エネルギーの法線方向成分、$W_{rel}^n$ は解放エネルギーの法線方向成分である。法線方向の成分である理由は、衝突時にエネルギーが散逸するのは相対速度の法線方向成分だからである（とのこと）。
2. 速度に関する接地・貫入条件を満たす粒子を探す。すべての粒子において、速度に関する非接触条件 $C_v(x_{iA},x_{jB},t) > 0$ を満たすまで、力積を与え反復的に速度の更新を行う。ただし、位置の更新は行わない。
3. 力積を与え速度の更新すると同時に、圧縮エネルギーを蓄積する。蓄積されたエネルギーを解放することで反発させる。
4. 粒子の速度を更新する際の優先度順位は、相対速度が小さい粒子から順に力積を与える。
エネルギー保存型力積法の具体的な計算方法
エネルギー保存型力積法の具体的な計算方法について述べる。エネルギー保存型力積法は、大きく2つの過程からなり、圧縮過程とエネルギー解放過程からなる。
時間ステップを $n$ 、圧縮過程の反復法における計算ステップを $k$ 、エネルギー解放過程の反復法における計算ステップを $l$ とする。
物体 A の粒子 $i$ と物体 B の粒子 $j$ が衝突したとする。衝突した粒子の中で、最も相対速度が小さい粒子の組み合わせが、物体 A の粒子 $i$ と物体 B の粒子 $j$ だとする。
時間ステップ $n$ 、 圧縮過程の $k$ ステップにおいて、そのとき生じた力積を$J_{ij}^{n(k)}$とする。物体 A,B における並進速度・角速度の増分は、
\begin{align}
\Delta \vec{V}_A^{n(k)} &= \frac{\vec{J}_{ij}^{n(k)} }{M_A} \ \ , \ \
\Delta \vec{\omega}_A^{n(k)} = I^{-1}_A(\vec{r}_{iA} \times \vec{J}_{ij}^{n(k)} ) \\
\Delta \vec{V}_B^{n(k)} &= -\frac{\vec{J}_{ij}^{n(k)} }{M_B} \ \ , \ \
\Delta \vec{\omega}_B^{n(k)} = -I^{-1}_B(\vec{r}_{jB} \times \vec{J}_{ij}^{n(k)} ) \\
\end{align}
と書ける。作用・反作用の関係を使っている。
この2点間の力積によって生じる相対速度の変化は、
\begin{align}
\Delta \vec{v}_{ij}^{n(k)} &\equiv \Delta \vec{v}_{iA}^{n(k)} - \Delta \vec{v}_{jB}^{n(k)} \\
& = (\Delta \vec{V}_A^{n(k)}+\Delta \vec{\omega}_A^{n(k)}\times \vec{r}_{iA} ) - (\Delta \vec{V}_B^{n(k)}+\Delta \vec{\omega}_B^{n(k)}\times \vec{r}_{jB} ) \\
& = \left(\frac{1}{M_A}+ \frac{1}{M_B}\right) \vec{J}_{ij}^{n(k)} + I^{-1}_A\left( \vec{r}_{iA} \times \vec{J}_{ij}^{n(k)}\right)\times \vec{r}_{iA}+ I^{-1}_B\left( \vec{r}_{jB} \times \vec{J}_{ij}^{n(k)}\right)\times \vec{r}_{jB}
\end{align}
と書くことができる。求めたいのは、力積 $J_{ij}^{n(k)}$ である。
外積に関する公式
\begin{align}
A \times (B \times C) +B \times (C \times A) + C \times (A \times B) = 0
\end{align}
を使えば、同じベクトルの外積はゼロなので、
\begin{align}
I^{-1} \vec{r} \times \left(\vec{J}\times \vec{r}\right) =- \vec{r} \times \left( I^{-1}\ \vec{r} \times \vec{J}\right)
\end{align}
となる。さらに、
\begin{align}
\vec{r} &= (r_x \ \ r_y \ \ r_z)^T \ \ , \ \ \vec{J} = (J_x \ \ J_y \ \ J_z)^T \\
\vec{r} \times \vec{J} &\equiv [r]_{\times} \vec{J} =
\begin{pmatrix}
0 & -r_z & r_y \\
r_z & 0 & -r_x \\
-r_y & r_x & 0 \\
\end{pmatrix}
\begin{pmatrix}
J_x\\
J_y \\
J_z \\
\end{pmatrix}
\end{align}
を使うと、
\begin{align}
\Delta \vec{v}_{ij}^{n(k)}
& = \left\{ \left(\frac{1}{M_A}+ \frac{1}{M_B}\right) {\bf{1}} - [r_{iA}]_{\times}I^{-1}_A[r_{iA}]_{\times}- [r_{jB}]_{\times}I^{-1}_B[r_{jB}]_{\times}
\right\}\vec{J}_{ij}^{n(k)}\\
& = K_{AB} \vec{J}_{ij}^{n(k)} \\
\therefore \ \ \vec{J}_{ij}^{n(k)} &= K_{AB}^{-1}\Delta \vec{v}_{ij}^{n(k)}
\end{align}
となり、力積 $J_{ij}^{n(k)}$ が計算できる。
後々必要となるので、力積を法線方向と接線方向に分解しておく。相対速度と力積は、
\begin{align}
\vec{v}_{ij}^{n(k)} = {v^n_{ij}}^{n(k)} \vec{n}_{jB} + {v^t_{ij}}^{n(k)} \vec{t}_{v} \\
\vec{J}_{ij}^{n(k)} = {J^n_{ij}}^{n(k)} \vec{n}_{jB} + {J^t_{ij}}^{n(k)} \vec{t}_{v} \\
\end{align}
と分解できる。新たに表れた添え字 $n,t$ は、法線・接線成分を意味している。
法線方向の速度成分は、相対速度に法線ベクトル $n_{jB}$ を掛ければ、
\begin{align}
{v^n_{ij}}^{n(k)} = \vec{v}_{ij}^{n(k)} \cdot \vec{n}_{jB}
\end{align}
と書けるので、残りの成分が接線方向に対応する。
\begin{align}
{v^t_{ij}}^{n(k)} \vec{t}_{v} = \vec{v}_{ij}^{n(k)} - {v^n_{ij}}^{n(k)} \vec{n}_{jB}
\end{align}
接線ベクトルは、規格化されているので二乗したら $1$ になる。よって、接線方向の速度成分は
\begin{align}
({v^t_{ij}}^{n(k)} \vec{t}_{v} )^2 &=( \vec{v}_{ij}^{n(k)} - {v^n_{ij}}^{n(k)} \vec{n}_{jB} )^2 \\
\therefore \ \ {v^t_{ij}}^{n(k)} &=\| \vec{v}_{ij}^{n(k)} - {v^n_{ij}}^{n(k)} \vec{n}_{jB} \| \\
\end{align}
と書ける。よって、速度成分における接線ベクトル $t_v$ は、
\begin{align}
\vec{t}_{v} =\frac{ \vec{v}_{ij}^{n(k)} - {v^n_{ij}}^{n(k)} \vec{n}_{jB} }{{v^t_{ij}}^{n(k)}} = \frac{ \vec{v}_{ij}^{n(k)} - {v^n_{ij}}^{n(k)} \vec{n}_{jB} }{\| \vec{v}_{ij}^{n(k)} - {v^n_{ij}}^{n(k)} \vec{n}_{jB} \|}
\end{align}
と書ける。力積成分における接線ベクトルは、速度成分における接線ベクトル $t_v$ と同一とする。
同様に、力積の法線・接線成分は、
\begin{align}
{J^n_{ij}}^{n(k)} &= \vec{J}_{ij}^{n(k)} \cdot \vec{n}_{jB} \\
{J^t_{ij}}^{n(k)} &=\left\| \vec{J}_{ij}^{n(k)} - {J^n_{ij}}^{n(k)} \vec{n}_{jB} \right\|
\end{align}
と求まる。
圧縮過程
上記の計算において力積を求める際に、相対速度の変化 $\Delta v_{ij}^{n(k)}$ を求める必要がある。接触時の相対速度を、ある一定の割合 $\alpha_{comp}$ で力積を与え解消する。割合は、接触時の法線方向の相対速度 ${v^n_{ij}}^{n(k)}$ とすれば、
\begin{align}
\alpha_{comp} \equiv \frac{{v^n_{ij}}^{n(k)}}{\Delta {v^n_{ij}}^{n(k)}} = \alpha_{0} N_{comp}^{n(k)}
\end{align}
で与えられる。$N_{comp}^{n(k)}$ は、$k$ 回目の圧縮過程の反復において、相対速度が最小となる粒子の組み合わせが複数発生した時の数である。接触時の相対速度は既知なので、相対速度の変化の法線方向成分 $\Delta {v_{ij}^n}^{n(k)}$は、
\begin{align}
\Delta {v_{ij}^n}^{n(k)} = -\frac{{v_{ij}^n}^{n(k)}}{\alpha_{0} N_{comp}^{n(k)}}
\end{align}
と計算する。負の符号は、非接触条件を満たすには反対方向に力積を加える必要があるからである。（相対速度は、接触条件から必ず負になる。反対方向に力積を加えるため、以下説明する、法線・接線方向の速度は正の符号になる。）
力積は、法線方向成分のみとして計算を行うと、
\begin{align}
\Delta \vec{v}_{ij}^{n(k)} &= K_{AB} \vec{J}_{ij}^{n(k)} \\
&=K_{AB} {J^n_{ij}}^{n(k)} \vec{n}_{jB}
\end{align}
両辺に法線ベクトルの内積をとれば
\begin{align}
\Delta {v_{ij}^n}^{n(k)} &= \vec{n}_{jB} \cdot \Delta \vec{v}_{ij}^{n(k)} \\
&={J^n_{ij}}^{n(k)}\vec{n}_{jB} K_{AB} \vec{n}_{jB} \\
\therefore \ \ {J^n_{ij}}^{n(k)} &= \frac{\Delta {v_{ij}^n}^{n(k)}}{\vec{n}_{jB} K_{AB} \vec{n}_{jB}}
\end{align}
となり、法線方向成分の力積が求まる。
接線方向成分の力積は、静止摩擦状態の場合と摩擦すべり状態のときの2通りを考える。
静止摩擦状態のときは、接線方向に運動しないことを意味する。よって、接線方向成分の相対速度を打ち消すように、反対方向に力積を与えればよいので、
\begin{align}
{\vec{v}_{ij}}^{n(k)} &= K_{AB} {J^t_{ij}}^{n(k)}\vec{t}_{v} \\
{v_{ij}^t}^{n(k)} &= {\vec{v}_{ij}}^{n(k)} \cdot \vec{t}_{v} \\
&={J^t_{ij}}^{n(k)}\vec{t}_{v} K_{AB} \vec{t}_{v} \\
\therefore \ \ {J^t_{ij}}^{n(k)} &= \frac{{v_{ij}^t}^{n(k)}}{\vec{t}_{v} K_{AB} \vec{t}_{v}}
\end{align}
となる。計算の途中で接線ベクトルの内積を計算した。
摩擦すべり状態のときは、クーロン摩擦モデルによると、
\begin{align}
\mu {J^n_{ij}}^{n(k)} < {J^t_{ij}}^{n(k)}
\end{align}
の不等式を満たす状態である。$\mu$ は摩擦係数である。このとき、接線方向成分の力積は、
\begin{align}
{J^t_{ij}}^{n(k)} = \mu {J^n_{ij}}^{n(k)}
\end{align}
と書くことができる。摩擦すべり状態のときは、再度、法線方向成分の力積を計算し
\begin{align}
\Delta \vec{v}_{ij}^{n(k)} &= K_{AB} \vec{J}_{ij}^{n(k)} \\
&=K_{AB}\left\{ {J^n_{ij}}^{n(k)} \vec{n}_{jB} + {J^t_{ij}}^{n(k)}\vec{t}_{v} \right\} \\
&={J^n_{ij}}^{n(k)} K_{AB}\left(\vec{n}_{jB} + \mu \vec{t}_{v} \right) \\
\Delta {v_{ij}^n}^{n(k)} &= \vec{n}_{jB} \cdot \Delta \vec{v}_{ij}^{n(k)} \\
&={J^n_{ij}}^{n(k)}\vec{n}_{jB} K_{AB}\left(\vec{n}_{jB} + \mu \vec{t}_{v} \right) \\
\therefore \ \ {J^n_{ij}}^{n(k)} &= \frac{\Delta {v_{ij}^n}^{n(k)}}{\vec{n}_{jB} K_{AB}\left(\vec{n}_{jB} + \mu \vec{t}_{v} \right)}
\end{align}
となる。
まとめると、
\begin{align}
{J^n_{ij}}^{n(k)} =& \frac{\Delta {v_{ij}^n}^{n(k)}}{\vec{n}_{jB} K_{AB} \vec{n}_{jB}} \ \ , \ \ {J^t_{ij}}^{n(k)} = \frac{{v_{ij}^t}^{n(k)}}{\vec{t}_{v} K_{AB} \vec{t}_{v}} \\
\vec{J}_{ij}^{n(k)} =&{J^n_{ij}}^{n(k)} \vec{n}_{jB} + {J^t_{ij}}^{n(k)} \vec{t}_{v} \\
=&
\left\{\begin{array}{ll}
\frac{\Delta {v_{ij}^n}^{n(k)}}{\vec{n}_{jB} K_{AB} \vec{n}_{jB}} \vec{n}_{jB}+\frac{{v_{ij}^t}^{n(k)}}{\vec{t}_{v} K_{AB} \vec{t}_{v}} \vec{t}_{v} & \mu {J^n_{ij}}^{n(k)} \geq {J^t_{ij}}^{n(k)}\mbox{ （静止摩擦状態）} \\
\frac{\Delta {v_{ij}^n}^{n(k)}}{\vec{n}_{jB} K_{AB}\left(\vec{n}_{jB} + \mu \vec{t}_{v} \right)} \left(\vec{n}_{jB} + \mu \vec{t}_{v} \right) & \mu {J^n_{ij}}^{n(k)} < {J^t_{ij}}^{n(k)} \mbox{ （摩擦すべり状態）}
\end{array}\right.
\end{align}
となる。
以上の力積を使い、
\begin{align}
\vec{V}_A^{n(k+1)} &= \vec{V}_A^{n(k)} + \Delta \vec{V}_A^{n(k)} \ \ , \ \
\Delta \vec{V}_A^{n(k)} = \frac{\vec{J}_{ij}^{n(k)} }{M_A} \\
\vec{\omega}_A^{n(k+1)} &=\vec{\omega}_A^{n(k)} + \Delta \vec{\omega}_A^{n(k)} \ \ , \ \
\Delta \vec{\omega}_A^{n(k)} = I^{-1}_A(\vec{r}_{iA} \times \vec{J}_{ij}^{n(k)} ) \\
\vec{V}_B^{n(k+1)} &= \vec{V}_B^{n(k)} + \Delta \vec{V}_B^{n(k)} \ \ , \ \
\Delta \vec{V}_B^{n(k)} = -\frac{\vec{J}_{ij}^{n(k)} }{M_B} \\
\vec{\omega}_B^{n(k+1)} &=\vec{\omega}_B^{n(k)} + \Delta \vec{\omega}_B^{n(k)} \ \ , \ \
\Delta \vec{\omega}_B^{n(k)} = -I^{-1}_B(\vec{r}_{jB} \times \vec{J}_{ij}^{n(k)} )
\end{align}
物体の並進・回転速度を更新する。
以上が、圧縮過程における力積の計算方法と粒子の速度の更新方法である。後に説明するエネルギー解放過程において、必要となる圧縮エネルギーを蓄積する必要がある。
圧縮エネルギーの法線方向成分は、物体を法線方向の力 $F^n$ で、法線方向の速度 $v^n$ で動かしたとき計算される量（仕事率）を時間で積分すれば得られ
\begin{align}
{W^n_{ij}}_{comp} = - \int_{t^{n(0)}}^{t^{n(k_{end})}} \mathrm{d}t F^n_{ij} v^n_{ij} = - \int_0^{{J^n_{ij}}^{n(k_{end})}} \mathrm{d}J_{ij} v^n_{ij}
\end{align}
と書ける。$k_{end}$ は、圧縮過程が終了した計算ステップである。負の符号は、非接触条件を満たすために力積が反対方向に作用するためである。離散化すれば、
\begin{align}
{W^n_{ij}}_{comp} &= - \int_0^{{J^n_{ij}}^{n(k_{end})}} \mathrm{d}J_{ij} v^n_{ij}\\
&\approx \sum_{k=1}^{k_{end}} \frac{1}{2} \left({v_{ij}^n}^{n(k)}+{v_{ij}^n}^{n(k+1)} \right) {J^n_{ij}}^{n(k)} \\
&= \sum_{k=1}^{k_{end}} \frac{1}{2} \left({v_{ij}^n}^{n(k)}+({v_{ij}^n}^{n(k)} +\Delta {v_{ij}^n}^{n(k)} )\right) {J^n_{ij}}^{n(k)}\\
&= \sum_{k=1}^{k_{end}} \frac{1}{2} \left(2{v_{ij}^n}^{n(k)}+\Delta {v_{ij}^n}^{n(k)} \right) {J^n_{ij}}^{n(k)}
\end{align}
と書ける。（相対速度 ${v_{ij}^n}^{n(k)}$ は、反対方向に力積を与えたいため、負の符号をつけていた。マイナスはそれと打ち消した。）
各ステップ $k$ における圧縮エネルギー
\begin{align}
{\Delta{W^n_{ij}}}^{n(k)}_{comp} &\equiv \frac{1}{2} \left(2{v_{ij}^n}^{n(k)}+\Delta {v_{ij}^n}^{n(k)} \right) {J^n_{ij}}^{n(k)}
\end{align}
を計算し
\begin{align}
{W^n_{ij}}_{comp} &=\sum_{k=1}^{k_{end}} {\Delta{W^n_{ij}}}^{n(k)}_{comp}
\end{align}
蓄積する。
圧縮過程における計算手順は、以下の通りである。
- 接触と判断された粒子の中で、最も相対速度が小さい組み合わせの粒子に力積を与え、粒子の速度を更新する。その際に、圧縮エネルギー ${\Delta{W^n_{ij}}}^{n(k)}_{comp}$ も蓄積する。
- 最も相対速度が小さい組み合わせの粒子の速度を更新した後、接触と判断された粒子に対して、速度に関する非接触条件 $C_v(x_{iA},x_{jB},t) > 0$ が満たされているか計算を行う。
- 接触と判断された粒子に対して、速度に関する非接触条件が満たされていないなら、最も相対速度が小さい組み合わせの粒子に力積を与え、粒子の速度を更新し、圧縮エネルギーを蓄積。これを速度に関する非接触条件が満たされるまで、1.、 2. を反復する。
エネルギーの散逸
エネルギーは、Stronge の仮定より
\begin{align}
W_{rel}^n = e^2 W_{comp}^n
\end{align}
が成立する。散逸したエネルギーは、
\begin{align}
W_{comp}^n - W_{rel}^n = W_{comp}^n-e^2 W_{comp}^n = (1-e^2)W_{comp}^n
\end{align}
である。
解放されるエネルギーは、
\begin{align}
{W^n_{ij}}_{rel} = \int_{t^{n(0)}}^{t^{n(l_{end})}} \mathrm{d}t F^n_{ij} v^n_{ij} = \int_{{{J^n_{ij}}^{n(0)}}}^{{J^n_{ij}}^{n(l_{end})}} \mathrm{d}J_{ij} v^n_{ij}
\end{align}
と書ける。$l_{end}$ は、エネルギー解放過程が終了した計算ステップである。離散化すれば、形式的に
\begin{align}
{W^n_{ij}}_{rel} &= \int_{{{J^n_{ij}}^{n(0)}}}^{{J^n_{ij}}^{n(l_{end})}} \mathrm{d}J_{ij} v^n_{ij}\\
&\approx \sum_{l=1}^{l_{end}} \frac{1}{2} \left({v_{ij}^n}^{n(l)}+{v_{ij}^n}^{n(l+1)} \right) {J^n_{ij}}^{n(l)} \\
&= \sum_{l=1}^{l_{end}} \frac{1}{2} \left({v_{ij}^n}^{n(l)}+({v_{ij}^n}^{n(l)} +\Delta {v_{ij}^n} )\right) {J^n_{ij}}^{n(l)}\\
&= \sum_{l=1}^{l_{end}} \frac{1}{2} \left(2{v_{ij}^n}^{n(l)}+\Delta {v_{ij}^n} \right) {J^n_{ij}}^{n(l)}\\
&=\sum_{l=1}^{l_{end}} {W^n_{ij}}^{n(l)}_{rel}
\end{align}
と書ける。
エネルギー解放過程
物体 A の粒子 $i$ と物体 B の粒子 $j$ が衝突しており、衝突した粒子の中で、最も解放エネルギーが大きい粒子の組み合わせが、物体 A の粒子 $i$ と物体 B の粒子 $j$ だとする。時間ステップ $n$ 、 エネルギー解放過程の $l$ ステップにおいて、そのとき生じた力積を$J_{ij}^{n(l)}$とする。
圧縮過程と同様に、相対速度の変化の法線方向成分 $\Delta {v_{ij}^n}^{n(l)}$ を求めたい。蓄積されたエネルギーは、
\begin{align}
{W^n_{ij}}^{n(l)}_{rel} &= \frac{1}{2} \left(2{v_{ij}^n}^{n(l)}+\Delta {v_{ij}^n} \right) {J^n_{ij}}\\
&= \frac{1}{2} \left(2{v_{ij}^n}^{n(l)}+\Delta {v_{ij}^n} \right)\frac{\Delta {v_{ij}^n}}{\vec{n}_{jB} K_{AB} \vec{n}_{jB}} \\
\end{align}
と書き直せる。蓄積されたエネルギーと相対速度は既知なので、これを解けば、
\begin{align}
&{W^n_{ij}}^{n(l)}_{rel} = \frac{1}{2} \left(2{v_{ij}^n}^{n(l)}+\Delta {v_{ij}^n}\right)\frac{\Delta {v_{ij}^n}}{\vec{n}_{jB} K_{AB} \vec{n}_{jB}} \\
\therefore \ \ & \left(\Delta {v_{ij}^n}\right)^2+2{v_{ij}^n}^{n(l)}\Delta {v_{ij}^n}-2{W^n_{ij}}^{n(l)}_{rel}\vec{n}_{jB} K_{AB} \vec{n}_{jB} =0 \\
\therefore \ \ & \Delta {v_{ij}^n} = - {v_{ij}^n}^{n(l)}+\sqrt{\left( {v_{ij}^n}^{n(l)}\right)^2 +2{W^n_{ij}}^{n(l)}_{rel}\vec{n}_{jB} K_{AB} \vec{n}_{jB} }
\end{align}
$\Delta {v_{ij}^n} $ が求まる。
圧縮過程と同様に、相対速度をある一定の割合 $\alpha_{rel}$ で力積を与え解消する。割合は、
\begin{align}
\alpha_{rel} \equiv \frac{\Delta{v^n_{ij}}}{\Delta {v^n_{ij}}^{n(l)}} = \alpha_{0}
\end{align}
で与えられる。よって、相対速度の変化の法線方向成分 は、
\begin{align}
\Delta {v^n_{ij}}^{n(l)}= \frac{1}{\alpha_{rel} } \Delta{v^n_{ij}}
\end{align}
と求まる。
力積の計算方法や速度の更新方法は、圧縮過程と同様に行う。
\begin{align}
{J^n_{ij}}^{n(l)} =& \frac{\Delta {v_{ij}^n}^{n(l)}}{\vec{n}_{jB} K_{AB} \vec{n}_{jB}} \ \ , \ \ {J^t_{ij}}^{n(l)} = \frac{{v_{ij}^t}^{n(l)}}{\vec{t}_{v} K_{AB} \vec{t}_{v}} \\
\vec{J}_{ij}^{n(l)} =&{J^n_{ij}}^{n(l)} \vec{n}_{jB} + {J^t_{ij}}^{n(l)} \vec{t}_{v} \\
=&
\left\{\begin{array}{ll}
\frac{\Delta {v_{ij}^n}^{n(l)}}{\vec{n}_{jB} K_{AB} \vec{n}_{jB}} \vec{n}_{jB}+\frac{{v_{ij}^t}^{n(l)}}{\vec{t}_{v} K_{AB} \vec{t}_{v}} \vec{t}_{v} & \mu {J^n_{ij}}^{n(l)} \geq {J^t_{ij}}^{n(l)}\mbox{ （静止摩擦状態）} \\
\frac{\Delta {v_{ij}^n}^{n(l)}}{\vec{n}_{jB} K_{AB}\left(\vec{n}_{jB} + \mu \vec{t}_{v} \right)} \left(\vec{n}_{jB} + \mu \vec{t}_{v} \right) & \mu {J^n_{ij}}^{n(l)} < {J^t_{ij}}^{n(l)} \mbox{ （摩擦すべり状態）}
\end{array}\right.
\end{align}
\begin{align}
\vec{V}_A^{n(l+1)} &= \vec{V}_A^{n(l)} + \Delta \vec{V}_A^{n(l)} \ \ , \ \
\Delta \vec{V}_A^{n(l)} = \frac{\vec{J}_{ij}^{n(l)} }{M_A} \\
\vec{\omega}_A^{n(l+1)} &=\vec{\omega}_A^{n(l)} + \Delta \vec{\omega}_A^{n(l)} \ \ , \ \
\Delta \vec{\omega}_A^{n(l)} = I^{-1}_A(\vec{r}_{iA} \times \vec{J}_{ij}^{n(l)} ) \\
\vec{V}_B^{n(l+1)} &= \vec{V}_B^{n(l)} + \Delta \vec{V}_B^{n(l)} \ \ , \ \
\Delta \vec{V}_B^{n(l)} = -\frac{\vec{J}_{ij}^{n(l)} }{M_B} \\
\vec{\omega}_B^{n(l+1)} &=\vec{\omega}_B^{n(l)} + \Delta \vec{\omega}_B^{n(l)} \ \ , \ \
\Delta \vec{\omega}_B^{n(l)} = -I^{-1}_B(\vec{r}_{jB} \times \vec{J}_{ij}^{n(l)} )
\end{align}
最後に、蓄積されたエネルギーを解放する。法線方向成分の力積と相対速度、相対速度の変化を使い、蓄積されたエネルギーを、
\begin{align}
{W^n_{ij}}^{n(l+1)}_{rel} &={W^n_{ij}}^{n(l)}_{rel} -\frac{1}{2} \left(2{v_{ij}^n}^{n(l)}+\Delta {v_{ij}^n}^{n(l)} \right) {J^n_{ij}}^{n(l)}\\
\end{align}
解放（更新）する。
エネルギー解放過程における計算手順は、以下の通りである。
- 接触と判断された粒子の中で、最も蓄積されたエネルギーが大きい組み合わせの粒子に力積を与え、粒子の速度を更新する。
- 蓄積されたエネルギーを解放する。
- 接触と判断された粒子に対して、蓄積されたエネルギーが残っているなら、最も蓄積されたエネルギーが大きい組み合わせの粒子に力積を与え、粒子の速度を更新し、蓄積されたエネルギーを解放する。これを蓄積されたエネルギーがゼロになるまで、1.、 2. を反復する。
位置の更新
圧縮過程とエネルギー解放過程のセットを実行し、すべての粒子に対して、速度に関する非接触条件 $C_v(x_{iA},x_{jB},t) > 0$ が満たされていないなら、圧縮過程とエネルギー解放過程のセットを繰り返し実行する。最終的に物体 A の速度は、
\begin{align}
\vec{V_A}^{n+1} &= \vec{V_A}^{n(l_{end})} \\
\vec{\omega_A}^{n+1} &=\vec{\omega_A}^{n(l_{end})}
\end{align}
と更新する。
すべての粒子に対して、速度に関する非接触条件 $C_v(x_{iA},x_{jB},t) > 0$ が満たされたなら、位置の更新を行う。
重心の位置の更新
物体 A の重心の位置は、エネルギー保存型力積法で求まった速度 $V_A^{n+1}$ を使い
\begin{align}
\vec{X_A}^{n+1} & = \vec{X_A}^{n} + \vec{V_A}^{n+1}\Delta t
\end{align}
と更新する。
重力がある場合、速度を
\begin{align}
\vec{V_A}^{n+1} & \leftarrow \vec{V_A}^{n+1} + \vec{g}\Delta t \ \ , \ \ \vec{g}=(0\ \ 0\ \ -g)
\end{align}
と更新する必要がある。場合によっては、エネルギー保存型力積法を行う前に重力による速度の更新を行う。
慣性テンソルと粒子の位置の更新
物体 A の回転表現を行うため、クォータニオンを用いる。クォータニオンの記事は多数あるので、詳細は省く。
クォータニオンの初期値は、
\begin{align}
\vec{q}^0=(1\ \ 0\ \ 0\ \ 0)
\end{align}
とする。
角速度 $\omega_A^{n+1}$ を用いると、クォータニオンの増分
\begin{align}
\Delta \vec{q}&=\left(\cos\left(\frac{\Delta \theta }{2} \right)\ \
\sin\left(\frac{\Delta \theta }{2} \right)n^{\Delta}_x\ \
\sin\left(\frac{\Delta \theta }{2} \right)n^{\Delta}_y\ \
\sin\left(\frac{\Delta \theta }{2} \right)n^{\Delta}_z \right) \\
n^{\Delta} &= \frac{\vec{\omega_A}^{n+1}}{\|\vec{\omega_A}^{n+1} \|} \ \ , \ \
\Delta \theta = \|\vec{\omega_A}^{n+1} \| \Delta t
\end{align}
が計算でき、クォータニオンの更新は、
\begin{align}
\vec{q}^{n+1}&= \Delta \vec{q} \wedge \vec{q}^{n} \\
\end{align}
と書ける。各成分は、
\begin{align}
q_0^{n+1} &= \Delta q_0 q_0^{n}-\Delta q_1 q_1^{n}-\Delta q_2 q_2^{n}-\Delta q_3 q_3^{n} \\
q_1^{n+1} &= \Delta q_0 q_1^{n}+\Delta q_1 q_0^{n}+\Delta q_2 q_3^{n}-\Delta q_3 q_2^{n} \\
q_2^{n+1} &= \Delta q_0 q_2^{n}+\Delta q_2 q_0^{n}-\Delta q_1 q_3^{n}+\Delta q_3 q_1^{n} \\
q_3^{n+1} &= \Delta q_0 q_3^{n}+\Delta q_3 q_0^{n}+\Delta q_1 q_2^{n}-\Delta q_2 q_1^{n} \\
\end{align}
である。クォータニオンのノルムは、$|\vec{q}|=1$ となる。そうでないときは、
\begin{align}
\vec{q} \leftarrow \frac{\vec{q}}{\|\vec{q}\|}
\end{align}
とする。クォータニオンを使えば、回転行列は
\begin{align}
R= &
\begin{pmatrix}
1-2(q_2^2+q_3^2) & 2(q_1q_2-q_0q_3) & 2(q_1q_3+q_0q_2) \\
2(q_1q_2+q_0q_3) & 1-2(q_1^2+q_3^2) & 2(q_2q_3-q_0q_1) \\
2(q_1q_3-q_0q_2) & 2(q_2q_3+q_0q_1) & 1-2(q_1^2+q_2^2) \\
\end{pmatrix}
\end{align}
と書ける。
回転行列を使えば、重心からの相対距離
\begin{align}
\vec{r_{iA}}^{0}= \vec{x}_{iA}^{0}-\vec{X_A}^{0}
\end{align}
として、物体 A における粒子 $i$ の位置は、
\begin{align}
\vec{x}_{iA}^{n+1} &= \vec{X_A}^{n+1} + \vec{r_{iA}}^{n+1}\\
&=\vec{X_A}^{n+1} + R^{n+1}\vec{r_{iA}}^{0}\\
\end{align}
と更新される。慣性テンソルは、
\begin{align}
I_A^{n+1} = R^{n+1}I_A^{0}{R^{n+1}}^T
\end{align}
と計算でき、法線ベクトルも
\begin{align}
\vec{n}_{iA}^{n+1} &= R^{n+1}\vec{n_{iA}}^{0}\\
\end{align}
と計算できる。
補足
重心の位置と慣性モーメントの計算法
今回の記事で必要となる、重心の位置と慣性モーメントの計算法について説明する。
物体を構成する粒子の密度を $\rho$ とし、粒子の位置を $x$ とすると、物体の重心 $X$ は
\begin{align}
\vec{X} & = \frac{\int_V \mathrm{d}^3x \rho \vec{x} }{\int_V \mathrm{d}^3x \rho }
\approx \frac{\sum_{i\in V} m_i \vec{x}_i }{M}
\end{align}
粒子の質量を $m$ 、物体の質量を $M$ とした。$i$ は粒子番号である。
慣性モーメントは、重心からの相対距離を $r=x-X$ として、
\begin{align}
I_{xx} &= \int_V \mathrm{d}^3x \rho (r^2_y +r^2_z) \approx \sum_{i\in V} m_i (r^2_{iy} +r^2_{iz}) \\
I_{yy} &= \int_V \mathrm{d}^3x \rho (r^2_z +r^2_x) \approx \sum_{i\in V} m_i (r^2_{iz} +r^2_{ix}) \\
I_{zz} &= \int_V \mathrm{d}^3x \rho (r^2_x +r^2_y) \approx \sum_{i\in V} m_i (r^2_{ix} +r^2_{iy}) \\
I_{yz} &=I_{zy} = -\int_V \mathrm{d}^3x \rho r_y r_z \approx -\sum_{i\in V} m_i r_{iy} r_{iz} \\
I_{xz} &=I_{zx} = -\int_V \mathrm{d}^3x \rho r_x r_z \approx -\sum_{i\in V} m_i r_{ix} r_{iz} \\
I_{yx} &=I_{xy} = -\int_V \mathrm{d}^3x \rho r_y r_x \approx -\sum_{i\in V} m_i r_{iy} r_{ix} \\
\end{align}
と計算される。
法線ベクトルを粒子群から求める方法
力積型個別要素法では、法線ベクトルを計算する必要がある。単純な形状や簡単な検証（衝突する粒子がある面に限られていたりする場合）ならば、プログラムで手動で設定すればよい。しかし、複雑な形状においては、数値的に求める必要がある。
粒子 $i$ における曲面は、接平面に近似できるので、粒子 $i$ の法線ベクトル $n$ は接平面の法線ベクトルとすればよい。
粒子 $i$ に近い曲面上の粒子 $j$ に対して、粒子 $i$ における接平面の法線ベクトルを $n$ とすれば、粒子 $j$ と接平面との距離は
\begin{align}
(\vec{x}_j - \vec{x}_i)\cdot \vec{n}=(x_j-x_i)n_x+(y_j-y_i)n_y+(z_j-z_i)n_z
\end{align}
と書ける。$x_i$ は粒子 $i$ の位置である。粒子 $i$ における曲面を接平面に近似したいため、粒子 $i$ に近い曲面上の粒子 $j$ は、接平面と十分近い距離のはずである。よって、求めたい法線ベクトルは、粒子 $j$ と接平面との距離が近くなるベクトルである。
粒子 $i$ に近い曲面上の粒子を$N$ 個集め、それらと接平面との距離が最小となるように目的関数を構築し最小化問題を解く。
法線ベクトルのノルムは、$1$ となる制約があるので、ラグランジュ未定乗数$\lambda$ を使い、以下の目的関数
\begin{align}
L(n,\lambda)=\sum_{j}^N \|(\vec{x}_j - \vec{x}_i)\cdot \vec{n} \|+\lambda\left(1-\|\vec{n} \| \right)
\end{align}
を最小化する $n$ を求める。
$\alpha,\beta,\gamma,\ldots$ を $x,y,z$ の指標とする。目的関数を微分して、極値を求めると
\begin{align}
\frac{\partial L}{\partial n_{\alpha}} &= \sum_{j}^N \frac{\partial }{\partial n_{\alpha}} \left\{(x_{j\beta} - x_{i\beta})n_{\beta}(x_{j\gamma} - x_{i\gamma})n_{\gamma} \right\}+\lambda\frac{\partial }{\partial n_{\alpha}} \left(1-n_{\beta}n_{\beta} \right)\\
&=2\sum_{j}^N (x_{j\beta} - x_{i\beta})n_{\beta}(x_{j\alpha} - x_{i\alpha})-2\lambda n_{\alpha} =0\\
\therefore \ \ &\sum_{j}^N (x_{j\alpha}-x_{i\alpha})(x_{j\beta}-x_{i\beta})n_{\beta}=\lambda n_{\alpha} \\
\end{align}
となり、行列を
\begin{align}
M= &
\begin{pmatrix}
x_{1}-x_{i} & x_{2}-x_{i} & \cdots & x_{N}-x_{i} \\
y_{1}-y_{i} & y_{2}-y_{i} & \cdots & y_{N}-y_{i} \\
z_{1}-z_{i} & z_{2}-z_{i} & \cdots & z_{N}-z_{i} \\
\end{pmatrix}
\end{align}
とすれば、
\begin{align}
M^TM \vec{n} = \lambda \vec{n}
\end{align}
とまとまる。よって、求めたい法線ベクトルは、$M^TM $ の固有ベクトルである。しかし、3つの固有ベクトルが求まってしまう。
目的関数に戻って計算すると、固有ベクトルは、規格化されているので$ |n|=1 $、
\begin{align}
L(n,\lambda)&=\sum_{j}^N \|(\vec{x}_j - \vec{x}_i)\cdot \vec{n} \|\\
&=\|M \vec{n} \| = (M \vec{n})^TM \vec{n} \\
&= \vec{n}^T M^T M \vec{n} \\
&= \vec{n}^T \lambda \vec{n} \\
&= \lambda
\end{align}
となる。よって、$M^TM $ の最小となる固有値に対応する固有ベクトルが、最終的に求めたい法線ベクトルとなる。
$M^TM $ は（自明に）対称行列なので、ヤコビ法を使い、固有値と固有ベクトルが求まる。
検証
上記で説明した理論の検証を行う。以下で用いる解析解は、前回の記事で紹介した。
法線方向成分の確認
法線方向成分の計算の確認を行う。
$2.0m \times 1.0m \times 2.0m$ の物体を高さ$0.6m$ から落下させる。数値計算と解析解との比較を行った。（解析解（厳密解）と呼んでいるが、実際は、自由落下を数値計算で求めて、壁に衝突した時のみ、速度を変えている。）
反発係数を $e=0.5$ として計算を行った。時間ステップは、$\Delta t = 5.0\times 10^{-4}$とした。反復法に必要なパラメータは、$\alpha_{0}=15$ とした。
以下の図は、時間と鉛直方向の速度のグラフである。
シミュレーションの後半になると、誤差が大きくなっているが、概ね、解析解と一致している。
接線方向成分の確認
次に接線方向成分の計算の確認を行う。
角度が $30$ 度の斜面に $2.0m \times 1.0m \times 2.0m$ の物体を滑らす。摩擦を $\mu =0.0,0.1,0.3,0.5$ として計算を行った。
他のパラメータは、反発係数は $e=0.0$ 、時間ステップは$\Delta t = 5.0\times 10^{-4}$、反復法に必要なパラメータは、$\alpha_{0}=15$ とした。
以下の図は、時間と $x$ 方向の速度のグラフである。
以下の図は、時間と $y$ 方向の速度のグラフである。
概ね、解析解と一致している。
物体の衝突
解析解
初速度 $v_A^0,v_B^0$ で滑っている2つの物体 A,B があり、その2つの物体が衝突したとする。衝突後、物体の速度は、$v_A,v_B$ となった。 質量は $m_A,m_B$ とする。
運動量保存則を適応すると
\begin{align}
m_Av_A^0+m_Bv_B^0=m_Av_A+m_Bv_B
\end{align}
となり、反発係数の定義から
\begin{align}
e = - \frac{v_B-v_A}{v_B^0-v_A^0 }
\end{align}
となる。これらの方程式から、衝突後の物体 A の速度 $v_A$ は、
\begin{align}
V_A = \frac{(m_A-m_B e)v_A^0+ m_B(1+e)v_B^0}{m_A+m_B}
\end{align}
衝突後の物体 B の速度 $v_B$ は、
\begin{align}
V_B = \frac{m_A(1+e)v_A^0+(m_B-m_A e)v_B^0 }{m_A+m_B}
\end{align}
となる。
$n-1$ 個の物体が横一列に並んで接触し、静止している。そこに、初速度 $v^0$ で横一列に並んでいる物体に衝突させる。もし、反発係数が $e=1$ の場合、衝突した物体の速度はゼロになる。最終的に、横一列に並んでいる最後の物体は、速度 $v^0$ で動く。
2個の物体の衝突
2つの物体を用意し、片方の物体は静止している。もう片方の物体は、初速度 $2.0m/s$ で運動しており、静止した物体に衝突させる。質量は同じにした。
摩擦は $\mu =0.0$ 、反発係数は $e=0.5$ 、時間ステップは$\Delta t = 5.0\times 10^{-4}$、反復法に必要なパラメータは、$\alpha_{0}=150$ とした。
以下の図は、2つの物体の時間と速度のグラフである。
概ね、解析解と一致している。
10個の物体の衝突
$9$ 個の物体が横一列に並んで接触し、静止している。そこに、初速度 $1.0m/s$ で横一列に並んでいる物体に衝突させる。質量は同じにした。
摩擦は $\mu =0.0$ 、反発係数は $e=1.0$ 、時間ステップは$\Delta t = 5.0\times 10^{-4}$、反復法に必要なパラメータは、$\alpha_{0}=150$ とした。床との接触計算に時間がかかるので、重力はゼロとした。
以下の動画は、シミュレーションの結果である。
以下の図は、衝突した物体と横一列に並んでいる最後の物体の、時間と速度のグラフである。
衝突してから、中間にある物体が力を伝え、最後の物体に力が伝わるまでに、時間がかかっている。衝突を繰り返しているので、誤差も大きくなっているが、概ね、解析解と一致している。
最後に
力積型個別要素法（iDEM）は、球形 DEM と違い、複雑な形状を再現させてシミュレーションを行うことができる。さらに、必要な物性値は、反発係数と摩擦係数である。
今後は、もう少し土木よりな現象をシミュレーションしてみたい。