概要

　土木の分野において、粒子法の研究が盛んにおこなわれています。例えば、離散要素法（DEM）、Smoothed Particle Hydrodynamics（SPH）、Material Point Method（MPM）などがあります。

　DEMは、昔から地盤力学の分野で使われている手法で、地盤解析・土砂崩れ・落石解析・流体と組み合わせれば地盤の越流現象など広い範囲で適応されています。

今回は、球体 DEM について簡単に紹介したいと思います。

以下の書籍を参考にしています。

DEM で用いる方程式

　粒子 $i$ における球体 DEM の並進・回転を記述する方程式は

\begin{align}
m_i \frac{\mathrm{d} \vec{v}_i}{\mathrm{d} t} &= \sum_j \vec{F}^C_{ij} + m_i \vec{g} \\
\frac{\mathrm{d} \vec{\omega}_i}{\mathrm{d} t} &= I_i^{-1}  \sum_j \vec{T}_{ij} \\
\end{align}

である。添え字は、粒子の番号を表す。$m_i$ は粒子 $i$ の質量であり、$v_i$ は速度、$\omega$ は角速度、$g$ は重力、$F^C_{ij}$は粒子 $i$ と接触した粒子 $j$ との接触力、$T_{ij}$ はモーメントである。
　
　今回の記事では、3次元の球体 DEM について説明するの。慣性モーメントは、球体なので

\begin{align}
I_i = \frac{2}{5}m_i R_i^2
\end{align}

と書ける。$R_i$ は粒子 $i$ の半径である。粒子は、変形をしないことを仮定し、剛体として扱う。

フォークトモデル（接触力のモデル）

　
　上記の方程式を（数値的に）解くことで、粒子の運動を記述できるが、そのためには接触力とモーメントを求める必要がある。接触力に関するモデルは多数あり、まずは、フォークトモデルを紹介する。

　フォークトモデルは、弾性力と粘性減衰を組み合したモデルであり、接触力は法線方向と接線方向の成分に分離して表される。法線方向の接触力を $F^{C_n}$ 、接線方向の接触力を $F^{C_t}$ とすれば、粒子 $i$ と接触した粒子 $j$ との接触力$F^C_{ij}$ は

\begin{align}
\vec{F^C}_{ij} = \vec{F^{C_n}}_{ij} + \vec{F^{C_t}}_{ij}
\end{align}

と書ける。法線方向の接触力と接線方向の接触力は、弾性力と粘性減衰を組み合して

\begin{align}
\vec{F^{C_n}}_{ij} &= -K_n \vec{\delta_n}_{ij} - \eta_n \vec{v_n}_{ij}    \\
\vec{F^{C_t}}_{ij} &= -K_t \vec{\delta_t}_{ij} - \eta_t \vec{v_t}_{ij}     \\
\end{align}

と書くことができる。１項目は弾性力、２項目は粘性減衰である。$K,\eta$ は、ばね定数、粘性減衰係数である。${\delta},{v} $は、相対変位、相対速度である。そして、添え字 $n,t$ は、法線方向成分、接線方向成分に対応している。

　最初に、相対速度の法線・接線方向成分を求める。法線方向成分の相対速度は、粒子 $i$ の重心の位置 $x_i$ 、さらに並進速度に回転速度を加えた速度を使い、

\begin{align}
\vec{v_n}_{ij} &= \left\{\left\{(\vec{v}_i+R_i \vec{\omega}_i \times \vec{n}_{ij} ) -(\vec{v}_j+R_j \vec{\omega}_j \times \vec{n}_{ji} )  \right\}  \cdot \vec{n}_{ij} \right\}\vec{n}_{ij}  \\
&= \left\{(\vec{v}_i-\vec{v}_j)\cdot \vec{n}_{ij} \right\}\vec{n}_{ij}　\\
\vec{n}_{ij} &= \frac{\vec{x}_i-\vec{x}_j}{|\vec{x}_i-\vec{x}_j|}
\end{align}

と書ける。$n$ は法線ベクトルである。ベクトルの三重積の性質

\begin{align}
A\cdot B \times C = C\cdot A \times B = B\cdot C \times A 
\end{align}

と同じベクトルの外積はゼロであることを用いた。回転速度の成分は消えているので、モーメントは、接線方向成分の力から与えられる。

　接線方向成分の相対速度は、並進速度に回転速度を加えた速度から、法線方向成分を引き

\begin{align}
\vec{v_t}_{ij} &= \left\{(\vec{v}_i+R_i \vec{\omega}_i \times \vec{n}_{ij} ) -(\vec{v}_j+R_j \vec{\omega}_j \times \vec{n}_{ji} )  \right\} - \left\{(\vec{v}_i-\vec{v}_j)\cdot \vec{n}_{ij} \right\}\vec{n}_{ij} \\
& = (\vec{v}_i- \vec{v}_j) - \left\{(\vec{v}_i-\vec{v}_j)\cdot \vec{n}_{ij} \right\}\vec{n}_{ij} + (R_i \vec{\omega}_i + R_j \vec{\omega}_j)\times \vec{n}_{ij}
\end{align}

と書くことができる。もちろん、接線方向成分の相対速度と法線ベクトルの内積を計算すると、

\begin{align}
\vec{v_t}_{ij}\cdot \vec{n}_{ij} &= (\vec{v}_i- \vec{v}_j)\cdot \vec{n}_{ij} - \left\{(\vec{v}_i-\vec{v}_j)\cdot \vec{n}_{ij} \right\} + (R_i \vec{\omega}_i + R_j \vec{\omega}_j)\times \vec{n}_{ij}\cdot \vec{n}_{ij}=0
\end{align}

となる。

　次に、相対変位の法線・接線方向成分を求める。法線方向成分の相対変位は、幾何学的な関係から求まる。

\begin{align}
\vec{\delta_n}_{ij} &= \left\{L_{ij}-(R_i + R_j )  \right\} \vec{n}_{ij}　\\
L_{ij} &= |\vec{x}_i-\vec{x}_j|
\end{align}

$L_{ij}$ は、粒子 $i$ と粒子 $j$ の重心間の距離である。粒子 $i$ と粒子 $j$ の半径を足した値より、重心間の距離の方が小さい場合、つまり、

\begin{align}
L_{ij}-(R_i + R_j ) < 0
\end{align}

ならば、接触していると言える。接触判定は上記の不等式を満たすかどうかで決める。粒子 $i$ と粒子 $j$が、接触判定を満たさないなら、接触力はゼロなので

\begin{align}
\vec{F^C}_{ij} =\vec{ 0}
\end{align}

である。

　接線方向成分の相対変位は、幾何学的な関係から導出するのは難しく、

\begin{align}
\vec{\delta_t}_{ij} &= \int_{t_0}^{t_f} \mathrm{d}t \vec{v_t}_{ij} \\
\end{align}

から計算する。$t_0$ は粒子 $i$ と粒子 $j$ が接触した時間で、$t_f$ は粒子 $i$ と粒子 $j$ が離れた時間である。この式を離散化すると、$n$ ステップ時の接線方向成分の相対変位は、

\begin{align}
\vec{\delta_t^n}_{ij} &= \vec{\delta_t^{n-1}}_{ij} +  \vec{v^n_t}_{ij} \Delta t\\
\end{align}

となる。しかし、$n$ ステップ時に与えたい接線方向成分の力の方向は、１ステップ前の $\delta_t^{n-1} $ の方向と同一ではない。よって、$n$ ステップ時の接線ベクトル

\begin{align}
\vec{t}_{ij} = 
\begin{cases}
\vec{\delta_t^{n-1}}_{ij} / |\vec{\delta_t^{n-1}}_{ij}| & (v_{ij} = 0)  \\
\vec{v}_{ij}/|\vec{v_{ij}|}    & (v_{ij} \neq 0)
\end{cases}
\end{align}

を使い、$n$ ステップ時の接線方向成分の相対変位を

\begin{align}
\vec{\delta_t^n}_{ij} &= |\vec{\delta_t^{n-1}}_{ij}| \vec{t}_{ij} +  \vec{v^n_t}_{ij} \Delta t\\
\end{align}

と計算する。

　以上の結果を使い、ばね定数・粘性減衰係数を与えれば、粒子 $i$ と接触した粒子 $j$ との接触力 $F^C_{ij}$ が求まる。

粘性減衰係数と反発係数との関係

　DEM では、粘性減衰係数を反発係数 $e$ を使い

\begin{align}
\eta = -2 \ln e \sqrt{\frac{m K}{\pi^2+(\ln e)^2} }
\end{align}

が用いられる。反発係数は、$0\leq e \leq 1$ の値をとるので、粘性減衰係数を正の値とするとマイナスの符号が付く。反発係数がゼロのときは、

\begin{align}
\eta = 2 \sqrt{{m K}}
\end{align}

である。以下、粘性減衰係数の導出を行う。
　
　粒子が壁に衝突した際、運動方程式は

\begin{align}
m\frac{\mathrm{d}^2 x}{\mathrm{d} t^2} + \eta\frac{\mathrm{d} x}{\mathrm{d} t}+Kx =0 
\end{align}

で表せる。この式の特性方程式の解は、

\begin{align}
\lambda_{1,2} &= \frac{-\eta \pm \sqrt{\eta^2-4mK}}{2m} =-\frac{\eta}{2m} \pm \sqrt{\frac{\eta^2}{4m^2} -\frac{K}{m}}\\
&\equiv -\frac{\eta}{2m} \pm i \Gamma  \\
\Gamma&=\sqrt{\frac{K}{m}-\frac{\eta^2}{4m^2}}
\end{align}

であるので、運動方程式の解は、任意の定数 $A,B$ を使い

\begin{align}
x(t) = A e^{\lambda_1 t} + B e^{\lambda_2 t}
\end{align}

と書ける。初期条件として、

\begin{align}
x(0) = 0  \ \ , \ \ \frac{\mathrm{d} x}{\mathrm{d} t} \Big|_{t=0} = v_a
\end{align}

を課す。すると

\begin{align}
&A+B = 0  \ \ , \ \ A\lambda_1+B\lambda_2 = v_a \\
\therefore & A=-B=\frac{v_a}{\lambda_1-\lambda_2 } = \frac{v_a}{2i \Gamma}
\end{align}

となり、最終的に解は、

\begin{align}
x(t) &= \frac{v_a}{ \Gamma}\frac{e^{i\Gamma t}-e^{-i\Gamma t} }{2i} \exp\left\{ 
  -\frac{\eta}{2m} t\right\} \\
&=\frac{v_a}{ \Gamma}\sin(\Gamma t) \exp\left\{ 
  -\frac{\eta}{2m} t\right\}
\end{align}

と求まる。$v_a$ は壁との衝突前の速度であり、壁と衝突した後の速度を $v_b$ とすると、反発係数は、

\begin{align}
e = -\frac{v_b}{v_a}
\end{align}

と表せる。壁と衝突した粒子は、壁と接触し反発をもらい壁から離れるのに半周期費やす。つまり、初期条件から $t=0$ のとき $x(0)=0$ であり、接触して$x(t)=0$ に戻ってくる最短の時刻は、$t=\pi/\Gamma$ となる。よって、壁と衝突した後の速度 $v_b$ は

\begin{align}
v(t) &={v_a}\cos(\Gamma t) \exp\left\{ 
  -\frac{\eta}{2m} t\right\} -\frac{v_a\eta}{2m \Gamma}\sin(\Gamma t) \exp\left\{ -\frac{\eta}{2m} t\right\}\\
\therefore & v_b =v(t=\pi/\Gamma) = -v_a \exp\left\{-\frac{\eta\pi}{2m\Gamma} \right\}
\end{align}

そして、反発係数は、

\begin{align}
e = \exp\left\{-\frac{\eta\pi}{2m\Gamma} \right\}
\end{align}

と表せる。これを $\eta$ について解くと

\begin{align}
&\ln e = -\frac{\eta\pi}{2m\Gamma} \\
\therefore & \left(\frac{K}{m}-\frac{\eta^2}{4m^2} \right) (\ln e)^2 = -\frac{\eta^2\pi^2}{4m^2} \\
\therefore & \ \ \eta^2 = \frac{4mK}{\pi^2 +(\ln e)^2 }(\ln e)^2 \\
\therefore & \ \  \eta = -2 \ln e \sqrt{\frac{m K}{\pi^2+(\ln e)^2} } \ \ (\eta \geq  0  )
\end{align}

が求まる。

ばね係数と物性値について

　有効ヤング率 $E^{\ast}$ と有効粒子半径 $R^{\ast}$ 、そして有効質量 $m^{\ast}$を

\begin{align}
\frac{1}{E^{\ast}_{ij}} &= \frac{1-\nu^2_i}{E_i} + \frac{1-\nu^2_j}{E_j}  \\
\frac{1}{R^{\ast}_{ij}} &= \frac{1}{R_i} + \frac{1}{R_j} \\
\frac{1}{m^{\ast}_{ij}} &= \frac{1}{m_i} + \frac{1}{m_j} 
\end{align}

とする。$E_i,\mu_i$ は、粒子 $i$ におけるヤング率・ポアソン比を表す。ヘルツの理論によると法線方向の相対変位 $\delta_n$ は

\begin{align}
\delta_{ij} = \left\{ \frac{9F_n^2}{16R^{\ast}_{ij}(E_{ij}^{\ast})^2} \right\}^{1/3}
\end{align}

と与えられる。$F_n$ は法線方向にかかる力である。この式を変形すると

\begin{align}
F_n &= -\frac{4}{3}E_{ij}^{\ast}\sqrt{R^{\ast}_{ij}}\delta_{ij}^{3/2}\\
&= -\frac{4}{3}\frac{E_iE_j}{(1-\nu^2_i)E_j+(1-\nu^2_j)E_i}\sqrt{\frac{R_iR_j}{R_i+R_j} } \delta_{ij}^{3/2}
\end{align}

となる。単純に、フックの法則 $F_n=-K_{ij}^n\delta_{ij} $ と比較すれば、ばね係数は

\begin{align}
K_{ij}^n=\frac{4}{3}\frac{E_iE_j}{(1-\nu^2_i)E_j+(1-\nu^2_j)E_i}\sqrt{\frac{R_iR_j}{R_i+R_j} } \delta_{ij}^{1/2}
\end{align}

となる。また、ばね係数のなかにある$\delta_{ij}$ を参照食い込み量 $\delta^c$ に置き換えた式

\begin{align}
K_{ij}^n=\frac{4}{3}\frac{E_iE_j}{(1-\nu^2_i)E_j+(1-\nu^2_j)E_i}\sqrt{\frac{R_iR_j}{R_i+R_j} } \sqrt{\delta^c}
\end{align}

も使われる。参照食い込み量 $\delta^c$ は、平均粒子半径の $5$パーセントの値で設定される。このとき、粘性減衰係数は、係数の中にある質量を有効質量に置き換え

\begin{align}
\eta^n_{ij} = -2 \ln e \sqrt{\frac{m^{\ast}_{ij} K_{ij}^n}{\pi^2+(\ln e)^2} } = -\frac{2 \ln e }{\sqrt{\pi^2+(\ln e)^2}} \sqrt{\frac{m_im_j}{m_i+m_j}K_{ij}^n} 
\end{align}

とする。接線方向のばね係数は、ミンドリンの接触理論より、

\begin{align}
K_{ij}^t=\frac{8E_iE_j}{(1+\nu_i)(2-\nu_i)E_j+(1+\nu_j)(2-\nu_j)E_i}\sqrt{\frac{R_iR_j}{R_i+R_j} } \sqrt{\delta^c}
\end{align}

粘性減衰係数は、

\begin{align}
\eta^t_{ij} = -2 \ln e \sqrt{\frac{m^{\ast}_{ij} K_{ij}^t}{\pi^2+(\ln e)^2} } = -\frac{2 \ln e }{\sqrt{\pi^2+(\ln e)^2}} \sqrt{\frac{m_im_j}{m_i+m_j}K_{ij}^t} 
\end{align}

とする。これらを使うと、法線方向の接触力と接線方向の接触力は、

\begin{align}
\vec{F^{C_n}}_{ij} &= -K_{ij}^n \vec{\delta_n}_{ij} -\eta^n_{ij} \vec{v_n}_{ij}    \\
\vec{F^{C_t}}_{ij} &= -K_{ij}^t\vec{\delta_t}_{ij} - \eta^t_{ij} \vec{v_t}_{ij}     \\
\end{align}

と書くことができる。

クーロン摩擦

　
　摩擦の効果をモデルに組み込む場合を考える。クーロン摩擦モデルにおいて

\begin{align}
\left|\vec{F^{C_t}}_{ij} \right| > \mu \left|\vec{F^{C_n}}_{ij} \right|  \\
\end{align}

を満たすとき、粒子表面においてすべりが生じる。$\mu$ は摩擦係数である。これは、最大静止摩擦力 $\mu |F^{C_n}| $ までは力を加えるが、それ以上の力は伝わらずすべりが生じる。

この場合、接線方向成分の接触力は

\begin{align}
\vec{F^{C_t}}_{ij} = -\mu \left|\vec{F^{C_n}}_{ij} \right| t_{ij} 
\end{align}

と表せる。すべりが発生する場合、接線方向成分の変位は変化しないので

\begin{align}
\vec{\delta_t^n}_{ij} &= \vec{\delta_t^{n-1}}_{ij} 
\end{align}

と更新する。

速度と位置の更新

　最初に、並進方向について考える。並進方向の速度を表す運動方程式は、

\begin{align}
m_i \frac{\mathrm{d} \vec{v}_i}{\mathrm{d} t} &= \sum_j \vec{F}^C_{ij} + m_i \vec{g} 
\end{align}

なので離散化を行い、時間ステップが $n+1$ ステップにおける速度の更新式は

\begin{align}
\vec{v}_i^{n+1} = \vec{v}_i^n + \frac{\Delta t}{m_i} \sum_j \vec{F}^C_{ij} + \vec{g}\Delta t
\end{align}

と書け、位置の更新式は

\begin{align}
\vec{x}_i^{n+1} = \vec{x}_i^n + \vec{v}_i^{n+1}\Delta t
\end{align}

と書ける。

　次に、回転方向について考える。粒子 $i$ と接触した粒子 $j$ における力のモーメントは、接線方向成分の力のみ寄与するので

\begin{align}
\vec{T}_{ij} = -R_i \vec{n}_{ij} \times \vec{F^{C_t}}_{ij}
\end{align}

と書くことができる。回転方向の速度を表す運動方程式は、

\begin{align}
\frac{\mathrm{d} \vec{\omega}_i}{\mathrm{d} t} &= I_i^{-1}  \sum_j \vec{T}_{ij} \\
\end{align}

なので離散化を行い、時間ステップが $n+1$ ステップにおける角速度の更新式は

\begin{align}
\vec{\omega}_i^{n+1} = \vec{\omega}_i^n + \Delta t I_i^{-1} \sum_j \vec{T}_{ij} 
\end{align}

と書ける。時間増分は、

\begin{align}
\Delta t < 2\sqrt{\frac{m}{k}}
\end{align}

の条件を満たすように決める場合が多い。詳細は、前回の記事で紹介した。

回転摩擦がある場合

　以上の説明は、粒子が球体の場合における DEM について説明した。ただし、実際の砂粒や岩石は、歪な形をしており球体 DEM で表現することは難しい。従って、球体 DEM において回転摩擦を加えることで、歪な形状の物体を表現する。

　回転摩擦力は、法線方向成分の力 $| F^{C_n}| $ に比例し、モーメントと同様、力の作用点を $R_i$ とし、力の方向を角速度の方向として、

\begin{align}
\vec{M}_{ij} =-\mu_r R_i \left| \vec{F^{C_n}}_{ij} \right|  \frac{\vec{\omega}_i-\vec{\omega}_j}{ \left| \vec{\omega}_i-\vec{\omega}_j \right|}
\end{align}

とする。モーメントと反対方向の力（回転の抗力）になるように負の符号をつける。$\mu_r$ は転がり摩擦係数である。角速度は

\begin{align}
\vec{\omega}_i^{n+1} = \vec{\omega}_i^n + \Delta t I_i^{-1} \sum_j \left( \vec{T}_{ij}+\vec{M}_{ij} \right) 
\end{align}

と更新する。転がり摩擦係数は、実験に合うように調整する必要がある。

壁との接触

　上記は、粒子同士の接触について取り扱った。壁を粒子で表現する方法もあるが、この記事では、平面の方程式を使い壁を表現する。

　壁の座標を $(x_w,y_w,z_w)$ とする。平面の方程式は、

\begin{align}
ax_w+by_w+cz_w+d =0
\end{align}

を満たす。また、この平面における法線ベクトルは、

\begin{align}
\vec{n} = (a,b,c)^T
\end{align}

と書くことができる。粒子位置を $\vec{x}=(x,y,z)$ として、壁が衝突する座標 $\vec{x_w}=(x_w,y_w,z_w)$ を求めたい。このときパラメータ $t$ を用いて

\begin{align}
(x_w,y_w,z_w) =(x,y,z)+ t(a,b,c)
\end{align}

と表され、平面の方程式に代入すると

\begin{align}
& a(x+ta)+b(y+tb)+c(z+tc)+d =0 \\
\therefore & \ \ t = -\frac{ax+by+cz+d }{a^2+b^2+c^2} 
\end{align}

パラメータ $t$ が求まる。壁との接触判定は、粒子 $i$ と壁との距離より粒子半径が大きい場合に接触とする。つまり、

\begin{align}
&\sqrt{(x_w-x_i)^2+(y_w-y_i)^2+(z_w-z_i)^2 } - R_i < 0 \\
\therefore & \ \ |t|\sqrt{a^2+b^2+c^2 }- R_i < 0 \\
\therefore & \ \ \frac{|ax_i+by_i+cz_i+d |}{\sqrt{a^2+b^2+c^2 }}- R_i < 0
\end{align}

である。

壁と接触する粒子 $i$ の法線方向成分の相対変位は、

\begin{align}
\vec{\delta_n}_{iw} &= \left\{\frac{|ax_i+by_i+cz_i+d |}{\sqrt{a^2+b^2+c^2 }}- R_i  \right\} \vec{n}_{iw}　\\
\vec{n}_{iw} &= \frac{\vec{x}_i-\vec{x}_w}{|\vec{x}_i-\vec{x}_w|}
\end{align}

である。

粒子同士の衝突と同様に、壁に衝突時の接線方向成分の相対変位を

\begin{align}
\vec{\delta_t}_{iw} &= \vec{v^n_t}_{iw} \Delta t\\
\end{align}

と計算する。法線方向・接線方向の相対速度は

\begin{align}
\vec{v_n}_{iw} &= \left\{(\vec{v}_i-\vec{v}_w)\cdot \vec{n}_{iw} \right\}\vec{n}_{iw}  \\
\vec{v_t}_{iw} &= (\vec{v}_i- \vec{v}_w) - \left\{(\vec{v}_i-\vec{v}_w)\cdot \vec{n}_{iw} \right\}\vec{n}_{iw} + R_i \vec{\omega}_i \times \vec{n}_{iw}
\end{align}

である。$\vec{v}_w$ は、壁が動く速さである。また、壁は回転しないとした。

　ばね係数と粘性減衰係数については、壁の半径や質量は無限大とする。法線方向のばね係数は

\begin{align}
K_{iw}^n=\frac{4}{3}\frac{E_iE_w}{(1-\nu^2_i)E_w+(1-\nu^2_w)E_i}\sqrt{R_i } \sqrt{\delta^c}
\end{align}

となり、粘性減衰係数は、

\begin{align}
\eta^n_{i_W} =  -\frac{2 \ln e }{\sqrt{\pi^2+(\ln e)^2}} \sqrt{m_iK_{iw}^n} 
\end{align}

となる。接線方向も同様に、ばね係数は、

\begin{align}
K_{iw}^t=\frac{8E_iE_w}{(1+\nu_i)(2-\nu_i)E_w+(1+\nu_w)(2-\nu_w)E_i}\sqrt{R_i } \sqrt{\delta^c}
\end{align}

粘性減衰係数は、

\begin{align}
\eta^t_{iw} = -\frac{2 \ln e }{\sqrt{\pi^2+(\ln e)^2}} \sqrt{m_iK_{iw}^t} 
\end{align}

となる。$E_w,\mu_w$ は、壁におけるヤング率・ポアソン比を表す。あとは、粒子同士の接触と同じ計算をする。

検証

　上記で説明した理論の検証を行う。以下で用いる解析解は、前回の記事で紹介した。

　法線方向成分の計算の確認を行う。

　質量$0.01 \mbox{kg}$ の物体を高さ$0.1 \mbox{m}$ から落下させる。数値計算と解析解との比較を行った。（解析解（厳密解）と呼んでいるが、実際は、自由落下を数値計算で求めて、壁に衝突した時のみ、速度を変えている。）

　反発係数を、$e=0.5,0.8,1.0$ として計算を行った。時間ステップは、$\Delta t = 10^{-6}$ とした。

　シミュレーション結果は、概ね、数値計算の結果と解析解の結果は一致している。ただし、シミュレーション後半は、解析解と数値解とのずれが少しばかり生じている。参考にした書籍では、$\Delta t = 10^{-7}$ であるので、さらに時間増分を小さくとる必要がある。

　以下の図は、時間と高さのグラフである。
　
　以下の図は、時間と鉛直方向の速度のグラフである。

　次に接線方向成分の計算の確認を行う。

　角度が$45$度の斜面に、質量$0.01 \mbox{kg}$ の球体を転がす（または滑らす）。摩擦を $\mu = 0.0,0.1,0.2,0.3,0.5$ として計算を行った。解析解では、

\begin{align}
\mu_c =  0.2857
\end{align}

より大きい摩擦では、球体は滑らず転がるため、摩擦を増やしても速度は減少しない。シミュレーション結果は、初期の粒子を配置するときに、斜面に対して少し高さを持ってしまったため、シミュレーション初期では、自由落下してしまっている。その後は、概ね、数値計算の結果と解析解の結果は一致している

　以下の図は、時間と $x$ 方向の速度のグラフである。

　以下の図は、時間と $y$ 方向の速度のグラフである。

最後に

　ヘルツの接触理論やミンドリンの接触理論について、もう少し勉強する必要がある。

　今回の記事では、球体 DEM について紹介を行った。歪な形状で計算を行うために、回転摩擦を導入するが、実際のシミュレーションでは、剛体の形状を与えて計算を行いたい。

　機会があれば、剛体の形状を与え、さらに、物性値は反発係数と摩擦係数のみで計算ができる力積型 DEM の紹介をしたい。