概要

　数値計算において、プログラムを構築する際、計算結果が正しいかどうか検証する必要があります。
　その際、解析解（または厳密解）があるモデルと比較することで、理論やプログラムの妥当性を確認できます。今回の記事では、DEM（個別要素法）において計算結果の確認でよく使われる、簡単な解析解の紹介をします。
　DEM では、法線方向と接線方向で力を分解するので、それぞれの方向で検討できる解を紹介します。

自由落下からの反発

　ここで説明する解は、個別要素法において、法線方向成分の計算結果の確認に使用できる。

　高さ $h_1$ から初速ゼロで物体を落下させ、壁に衝突させることを考える。壁に衝突し反発した物体は、$h_2$ の高さまで到達したとする。$h_2$ の高さまで到達した時、速度はゼロとする。反発係数 $e$ と初期の物体の高さ $h_1$ が与えられたとし、$h_2$ を求めたいとする。

　エネルギー保存則が成立するので、

\begin{align}
\frac{1}{2} mv_1^2 = mgh_1
\end{align}

が成立する。$v_1$ は、壁に衝突寸前の速度である。$m$ は物体の質量で、重力を $g$ とするが負の符号はつけない。速度 $v_1$ は、

\begin{align}
v_1 = -\sqrt{2gh_1}
\end{align}

と求まる。落下しているので、負の符号をつけた。壁に衝突し反発した物体の速度を$v_2$とすると、反発係数 $e$ の定義から

\begin{align}
e &= - \frac{v_2}{v_1}　\\
\therefore v_2 & = \sqrt{2gh_1} e
\end{align}

となる。エネルギー保存則の式を作り、$v_2$ を代入すると

\begin{align}
&\frac{1}{2} mv_2^2 = mgh_2 \\
\therefore \ \ & h_2  = e^2h_1
\end{align}

高さ $h_2$ が求まる。
　
　高さ $h_1$ から初速ゼロで物体を落下させ、物体が壁に $n$ 回衝突を繰り返し、反発して高さ $h_{n+1}$ まで到達したとする。$h_{n+1}$ の高さまで到達した時、速度はゼロとする。高さ $h_{n+1}$ は、

\begin{align}
h_{n+1}  = e^2 h_{n} = \cdots =  e^{2n} h_{1} 
\end{align}

と求まる。

　高さ $h_{n+1}$ と数値計算の結果が一致しているかどうかにより、理論やプログラムの確認ができる。

摩擦がある斜面を滑る物体

　ここで説明する解は、個別要素法において、接線方向成分の計算結果の確認に使用できる。

　摩擦がある斜面を滑る物体の運動について考える。初速をゼロとする。斜面の角度を $\theta$ とすれば、斜面に対して物体にかかる力は、

\begin{align}
F &=m\frac{\mathrm{d} v}{\mathrm{d} t } = mg\sin \theta - \mu N \\
N &= mg\cos \theta 
\end{align}

と書ける。$m$ は物体の質量で、$\mu$ は摩擦定数、$g$ は重力である。$N$ は垂直抗力である。物体の速度 $v$ は、初速度がゼロなので

\begin{align}
\frac{\mathrm{d} v}{\mathrm{d} t } &= g\sin \theta - \mu  g\cos \theta \\
 \therefore v(t) &=  gt(\sin \theta - \mu \cos \theta)
\end{align}

位置 $x$ は、

\begin{align}
\frac{\mathrm{d} x}{\mathrm{d} t } &= v(t) \\
 \therefore x(t) &= x_0 + \frac{1}{2}gt^2(\sin \theta - \mu \cos \theta)
\end{align}

と求まる。$x_0$ は初期の物体の位置である。上記から得られた速度と位置が、数値計算の結果と一致しているかどうかにより、理論やプログラムの確認ができる。

静止摩擦係数 $\mu_s$ において、物体が斜面を滑る力と摩擦力が釣り合うので

\begin{align}
F &=0 \\
\therefore  \mu_s &= \tan \theta  
\end{align}

となる。摩擦係数を $\mu > \tan \theta $ を満たすように設定してしまうと、物体は斜面を滑らない。

斜面を転がる物体が球体の場合

　斜面を転がる物体が球体の場合を考える。並進方向と回転方向の運動方程式は

\begin{align}
m\frac{\mathrm{d} v}{\mathrm{d} t } &= mg\sin \theta - f_s \\
I \frac{\mathrm{d} \omega}{\mathrm{d} t} &= rf_s 
\end{align}

と書ける。$f_s$ は摩擦力、$I$ は慣性モーメント、$\omega$ は角速度、$r$ は球体の半径である。慣性モーメントは球体なので

\begin{align}
I = \frac{2}{5}mr^2
\end{align}

と書ける。速度と角速度の関係は

\begin{align}
v = r \omega 
\end{align}

なので、微分をして、回転方向の運動方程式に代入すれば

\begin{align}
\frac{\mathrm{d} v}{\mathrm{d} t} &= r \frac{\mathrm{d}\omega }{\mathrm{d} t}　\\
\therefore f_s &= \frac{I}{r}\frac{\mathrm{d}\omega }{\mathrm{d} t} =\frac{I}{r^2}\frac{\mathrm{d}v }{\mathrm{d} t} 
\end{align}

摩擦力が求まる。この式を並進方向の運動方程式に代入すれば、

\begin{align}
m\frac{\mathrm{d} v}{\mathrm{d} t } &= mg\sin \theta -\frac{I}{r^2}\frac{\mathrm{d}v }{\mathrm{d} t} \\
\therefore \frac{\mathrm{d}v }{\mathrm{d} t} &= \frac{mgr^2\sin\theta}{mr^2 + I } \\
\therefore v(t) &= \frac{mgr^2\sin\theta}{mr^2 + I }t \\
\therefore x(t) &=x_0+ \frac{1}{2}\frac{mgr^2\sin\theta}{mr^2 + I }t^2 \\
\end{align}

速度と位置が求まる。$x_0$ は初期の物体の位置であり、物体の初速はゼロとした。

斜面を転がる球体が回転せず滑る条件

　斜面を転がる球体が回転せず滑るには、摩擦力が

\begin{align}
f_s =  \mu N \ \ , \ \ N = mg\cos \theta 
\end{align}

とかけ、回転モーメントと釣り合うときである。そのときの摩擦係数 $\mu$ は、

\begin{align}
 f_s &= \frac{I}{r^2}\frac{\mathrm{d}v }{\mathrm{d} t} =\frac{I}{r^2}\frac{mgr^2\sin\theta}{mr^2 + I } \\
& = \mu mg\cos \theta \\
\therefore \mu &= \frac{I}{mr^2 + I } \tan\theta
\end{align}

である。つまり、摩擦係数が

\begin{align}
\mu \leq \frac{ I }{mr^2 + I }\tan\theta
\end{align}

を満たすとき、球体は回転せず斜面をすべる。つまり、物体の速度は、

\begin{align}
v(t) &=  gt(\sin \theta - \mu \cos \theta)
\end{align}

である。逆に、摩擦係数が

\begin{align}
\mu > \frac{ I }{mr^2 + I }\tan\theta
\end{align}

を満たす場合、速度は摩擦係数に依らず

\begin{align}
v(t) &= \frac{mgr^2\sin\theta}{mr^2 + I }t \\
\end{align}

となる。球体の場合、摩擦係数を大きく取っても、静止することなく転がることを意味する。

最後に

　今回の記事は、基礎物理の復習に役立った。