概要
数値計算において、プログラムを構築する際、計算結果が正しいかどうか検証する必要があります。
その際、解析解（または厳密解）があるモデルと比較することで、理論やプログラムの妥当性を確認できます。今回の記事では、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}
となる。球体の場合、摩擦係数を大きく取っても、静止することなく転がることを意味する。
最後に
今回の記事は、基礎物理の復習に役立った。