概要
時間発展する方程式を陽解法で解くとき、時間増分を上手く決める必要があります。
個別要素法（DEM）において、時間増分 $\Delta t$ は
\begin{align}
\Delta t < 2\sqrt{\frac{m}{k}}
\end{align}
を満たすように設定されます。今回の記事では、この条件の導出を行います。
運動方程式と差分化
DEMの場合、フォークトモデルと呼ばれる、弾性力 $(\propto x)$ と粘性減衰 $(\propto dx/dt)$ を組み合されたモデルが使われる。よって、運動方程式は
\begin{align}
m\frac{\mathrm{d}^2 x}{\mathrm{d} t^2} + c\frac{\mathrm{d} x}{\mathrm{d} t}+kx =0
\end{align}
の形で書ける。$m,c,k$ は質量、粘性減衰係数、ばね定数である。
この方程式を、単純に、離散化すると
\begin{align}
m\frac{x(n+2)-2x(n+1)+x(n)}{\Delta t^2} + c\frac{x(n+1)-x(n)}{\Delta t}+kx(n+1) =0
\end{align}
となる。$n$ は時間ステップ数である。 この式を整理すれば、
\begin{align}
\begin{pmatrix}
x(n+2) \\
x(n+1) \\
\end{pmatrix}= &
\begin{pmatrix}
2-\frac{c}{m}\Delta t-k\frac{\Delta t^2}{m} & \frac{c}{m}\Delta t-1 \\
1 & 0 \\
\end{pmatrix}
\begin{pmatrix}
x(n+1) \\
x(n) \\
\end{pmatrix} \\
= &
\begin{pmatrix}
2-\frac{c}{m}\Delta t-k\frac{\Delta t^2}{m} & \frac{c}{m}\Delta t-1 \\
1 & 0 \\
\end{pmatrix}^2
\begin{pmatrix}
x(n) \\
x(n-1) \\
\end{pmatrix} \\
& \vdots \\
=& \begin{pmatrix}
2-\frac{c}{m}\Delta t-k\frac{\Delta t^2}{m} & \frac{c}{m}\Delta t-1 \\
1 & 0 \\
\end{pmatrix}^n
\begin{pmatrix}
x(1) \\
x(0) \\
\end{pmatrix} \\
\end{align}
と書ける。
初期条件として、$x(1),x(0)$ の値を与えれば、時間ステップ $n$ における解 $x(n)$ が求まる。
解の安定条件
以下の行列
\begin{align}
\begin{pmatrix}
2-\frac{c}{m}\Delta t-k\frac{\Delta t^2}{m} & \frac{c}{m}\Delta t-1 \\
1 & 0 \\
\end{pmatrix}
\end{align}
の対角化を考える。固有値が縮退しなければ（２つの固有値が異なれば）対角化可能である。
上記の行列の固有方程式は、
\begin{align}
\lambda \left(\lambda-2+\frac{c}{m}\Delta t+k\frac{\Delta t^2}{m}\right)-\frac{c}{m}\Delta t+1 =0
\end{align}
であり、これを解くと
\begin{align}
\lambda_{\pm} = \frac{1}{2}\left\{2-\frac{c}{m}\Delta t-k\frac{\Delta t^2}{m} \pm \sqrt{\left(-2+\frac{c}{m}\Delta t+k\frac{\Delta t^2}{m}\right)^2+4\left(\frac{c}{m}\Delta t-1\right) } \right\}
\end{align}
固有値が求まる。平方根の中がゼロ以上だとすれば、縮退はしていないので、この行列は対角化できる。直交行列を$R$ とすれば、
\begin{align}
\begin{pmatrix}
2-\frac{c}{m}\Delta t-k\frac{\Delta t^2}{m} & \frac{c}{m}\Delta t-1 \\
1 & 0 \\
\end{pmatrix} = R^{-1}
\begin{pmatrix}
\lambda_{+} & 0 \\
0 & \lambda_{-} \\
\end{pmatrix} R
\end{align}
と書くことができ、離散化された方程式は、
\begin{align}
\begin{pmatrix}
x(n+2) \\
x(n+1) \\
\end{pmatrix}= &
R^{-1}
\begin{pmatrix}
\lambda_{+} & 0 \\
0 & \lambda_{-} \\
\end{pmatrix} R
\begin{pmatrix}
x(n+1) \\
x(n) \\
\end{pmatrix} \\
= &
R^{-1}
\begin{pmatrix}
\lambda_{+}^2 & 0 \\
0 & \lambda_{-}^2 \\
\end{pmatrix} R
\begin{pmatrix}
x(n) \\
x(n-1) \\
\end{pmatrix} \\
& \vdots \\
=&R^{-1}
\begin{pmatrix}
\lambda_{+}^n & 0 \\
0 & \lambda_{-}^n \\
\end{pmatrix} R
\begin{pmatrix}
x(1) \\
x(0) \\
\end{pmatrix} \\
\end{align}
と書くことができる。従って、以下の条件
\begin{align}
-1 < \lambda_{-} \ \ , \ \ \lambda_{+} < 1
\end{align}
を満たすなら、解は安定すると言える。
最初に $\lambda_{+} < 1 $ を満たす$\Delta t$ を求める。計算をしていくと
\begin{align}
&\frac{1}{2}\left\{2-\frac{c}{m}\Delta t-k\frac{\Delta t^2}{m} + \sqrt{\left(-2+\frac{c}{m}\Delta t+k\frac{\Delta t^2}{m}\right)^2+4\left(\frac{c}{m}\Delta t-1\right) } \right\} < 1 \\
& \therefore \sqrt{\left(-2+\frac{c}{m}\Delta t+k\frac{\Delta t^2}{m}\right)^2+4\left(\frac{c}{m}\Delta t-1\right) } < \frac{c}{m}\Delta t+k\frac{\Delta t^2}{m}
\\
& \therefore
\left(-2+\frac{c}{m}\Delta t+k\frac{\Delta t^2}{m}\right)^2+4\left(\frac{c}{m}\Delta t-1\right) < \left(\frac{c}{m}\Delta t+k\frac{\Delta t^2}{m} \right)^2
\\
& \therefore 4-4\left(\frac{c}{m}\Delta t+k\frac{\Delta t^2}{m} \right)+4\left(\frac{c}{m}\Delta t-1\right) < 0 \\
& \therefore \ \ k\frac{\Delta t^2}{m} > 0
\end{align}
条件が求まる。質量とばね定数は、正の値に設定するので、常に $\lambda_{+} < 1 $ が満たされる。
次に $-1 < \lambda_{-} $ を満たす$\Delta t$ を求める。計算をしていくと
\begin{align}
&-1 < \frac{1}{2}\left\{2-\frac{c}{m}\Delta t-k\frac{\Delta t^2}{m} - \sqrt{\left(-2+\frac{c}{m}\Delta t+k\frac{\Delta t^2}{m}\right)^2+4\left(\frac{c}{m}\Delta t-1\right) } \right\} \\
& \therefore -\left(4-\frac{c}{m}\Delta t-k\frac{\Delta t^2}{m} \right) < -\sqrt{\left(-2+\frac{c}{m}\Delta t+k\frac{\Delta t^2}{m}\right)^2+4\left(\frac{c}{m}\Delta t-1\right) } \\
& \therefore \left(-2+\frac{c}{m}\Delta t+k\frac{\Delta t^2}{m}\right)^2+4\left(\frac{c}{m}\Delta t-1\right) < \left(4-\frac{c}{m}\Delta t-k\frac{\Delta t^2}{m} \right)^2 \\
& \therefore 4+4\left(\frac{c}{m}\Delta t-1\right)<16-4\left(\frac{c}{m}\Delta t+k\frac{\Delta t^2}{m} \right) \\
& \therefore \ \ \Delta t^2 +2\frac{c}{k} \Delta t-4\frac{m}{k}<0 \\
& \therefore \ \ \Delta t < -\frac{c}{k} + \sqrt{\frac{c^2}{k^2}+4\frac{m}{k}} \ \ (\Delta t > 0)
\end{align}
条件が求まる。
以上より、
\begin{align}
\Delta t < -\frac{c}{k} + \sqrt{\frac{c^2}{k^2}+4\frac{m}{k}}
\end{align}
を満たすように、時間増分 $\Delta t $ を決めればよい。
最後に
減衰定数 $\beta$ は、
\begin{align}
\beta = \frac{c}{2\sqrt{mk}}
\end{align}
と定義され、減衰定数を使い時間増分の条件式に代入すれば、
\begin{align}
\Delta t &< -2\beta\sqrt{\frac{m}{k}} + \sqrt{4\beta^2\frac{m}{k}+4\frac{m}{k}}\\
& = 2\sqrt{\frac{m}{k}} \left(\sqrt{\beta^2+1} -\beta \right)
\end{align}
が得られる。$\beta=0$ とすれば
\begin{align}
\Delta t &< 2\sqrt{\frac{m}{k}}
\end{align}
概要に記載した式が得られる。
（実際のDEMシミュレーションにおいて、この条件式を満たすように時間増分を設定しても、粒子が小さすぎたりすると壁を通り抜けたりするので、さらに小さい値に設定する場合がある。）