概要

　時間発展する方程式を陽解法で解くとき、時間増分を上手く決める必要があります。

　個別要素法（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シミュレーションにおいて、この条件式を満たすように時間増分を設定しても、粒子が小さすぎたりすると壁を通り抜けたりするので、さらに小さい値に設定する場合がある。）