概要

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

この記事では、MPSについて、簡単に説明したいと思います。

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

重み関数

　
　粒子法は、基本的に、計算したい粒子 $i$ とその周りの粒子 $j$ の物理量を使い、粒子 $i$ に関する物理量の勾配を計算する。つまり、偏微分方程式に出てくる微分を計算する。

　MPS の場合、以下の重み関数を使い加重平均を行い、物理量の勾配を求める。

\begin{align}
\omega(|r_{ij}|) = 
\begin{cases}
\frac{r_e}{|r_{ij}|}-1 & (0 \leq |r_{ij}| \leq r_e)  \\
0             & (r_e \leq |r_{ij}|)
\end{cases}
\end{align}

また $r_{ij}$ は、粒子 $j$ における半径 $r_j$ と粒子 $i$ における半径 $r_i$ と
の差、 $r_{ij} = r_j - r_i $ である。半径 $r_e$ は、粒子間距離の2～4倍程度の値に設定する。

　粒子 $i$ における数密度は、重み関数を使い

\begin{align}
n_i = \sum_{j \neq i} \omega(|r_{ij}|)
\end{align}

とする。以下、勾配・ラプラシアンの具体的な計算方法を説明していく。

勾配の計算

　物理量 $\phi$ の勾配 ∇$\phi$ の計算方法について説明する。
　2点間、粒子 $i$ と粒子 $j$ における物理量 $\phi$ の勾配は、粒子 $i$ と粒子 $j$における物理量 $\phi_i,\phi_j$ および半径 $r_i,r_j$ を用いて

\begin{align}
(\nabla \phi)_{ij} = \frac{\phi_j -\phi_i}{|r_j - r_i|} \frac{r_j -r_i}{|r_j - r_i|}
\end{align}

と書ける。右辺の1つ目の分数は、傾きを意味しており、2つ目の分数は方向を意味している。

　粒子 $i$における勾配は、重み関数を使い加重平均を行い

\begin{align}
\langle \nabla \phi\rangle_{i} &= \frac{1}{n_i}\sum_{j \neq i} \omega(|r_{ij}|)(\nabla \phi)_{ij}  \\
&= \frac{1}{n_i} \sum_{j \neq i} \omega(|r_{ij}|) \frac{\phi_{ij}}{|r_{ij}|^2} r_{ij}
\end{align}

となる。$\phi_{ij} = \phi_j -\phi_i $ である。
この式は、1次元方向しか考えていないので、どの方向も等方的かつ均質であると仮定し（等方的なのでどの方向も上記の式、均質なので（質量は保たれるので）$n_i\approx n_0 $）、粒子 $i$における勾配を、$n_i \to n_0$として

\begin{align}
\langle \nabla \phi\rangle_{i} 
&= \frac{D}{n_0} \sum_{j \neq i} \omega(|r_{ij}|) \frac{\phi_{ij}}{|r_{ij}|^2} r_{ij}
\end{align}

とする。$D$ は次元数で、$n_0$は初期状態の数密度である。初期状態の数密度$n_0$ は、シミュレーションの初期に計算される。

ラプラシアンの計算

　次にラプラシアンの計算を行う。極座標表示においてラプラシアンは、

\begin{align}
 \nabla^2 \phi 
&= \frac{1}{r^2}\frac{\partial}{\partial r} \left(r^2\frac{\partial \phi}{\partial r}  \right) = \frac{1}{r^2}\left(2r\frac{\partial \phi}{\partial r} + r^2 \frac{\partial^2 \phi}{\partial r^2}  \right) 
\end{align}

である。$r$ の2階微分を離散化したいため、$\phi_j$ を $\phi_i$ でテイラー展開し、3次以上は無視するとして

\begin{align}
&\phi_j = \phi_j + \frac{\partial \phi_{ij}}{\partial r_{ij}}r_{ij} +\frac{1}{2}\frac{\partial^2 \phi_{ij}}{\partial r_{ij}^2}r_{ij}^2 + \cdots \\
&\frac{\partial^2 \phi_{ij}}{\partial r_{ij}^2}r_{ij}^2 = 2\phi_{ij} - 2\frac{\partial \phi_{ij}}{\partial r_{ij}}r_{ij}
\end{align}

と書けるので、2点間、粒子 $i$ と粒子 $j$における物理量$\phi$ のラプラシアンは、

\begin{align}
(\nabla^2 \phi)_{ij} &= \frac{1}{r_{ij}^2}\left(2r_{ij}\frac{\partial \phi_{ij}}{\partial r_{ij}} + 2\phi_{ij} - 2\frac{\partial \phi_{ij}}{\partial r_{ij}}r_{ij}  \right) \\
&=\frac{ 2\phi_{ij}}{r_{ij}^2}
\end{align}

と書けるので、等方性を仮定して、粒子 $i$におけるラプラシアンは、

\begin{align}
\langle \nabla^2 \phi\rangle_{i} 
&= \frac{2D}{n_i} \sum_{j \neq i} \omega(|r_{ij}|) \frac{\phi_{ij}}{|r_{ij}|^2} 
\end{align}

となる。また、上式の右辺を書き換えて、

\begin{align}
\frac{1}{n_i} \sum_{j \neq i} \omega(|r_{ij}|) \frac{\phi_{ij}}{|r_{ij}|^2} &=  
\frac{\sum_{j \neq i} \omega(|r_{ij}|) \frac{\phi_{ij}}{|r_{ij}|^2}}{\sum_{j \neq i} \omega(|r_{ij}|)}　\\
&\rightarrow \frac{\sum_{j \neq i} \omega(|r_{ij}|) \phi_{ij}}{\sum_{j \neq i} \omega(|r_{ij}|)|r_{ij}|^2}　\\
\end{align}

\begin{align}
\lambda  \equiv  \frac{\sum_{j \neq i} \omega(|r_{ij}|)|r_{ij}|^2}{\sum_{j \neq i} \omega(|r_{ij}|)}
\\ 
\end{align}

置けば、粒子 $i$におけるラプラシアンは、

\begin{align}
\langle \nabla^2 \phi\rangle_{i} 
&= \frac{2D}{\lambda n_0} \sum_{j \neq i} \omega(|r_{ij}|) \phi_{ij} 
\end{align}

となる。

　$\lambda$ は、$n_0$と同様にシミュレーションの初期に計算される。

圧力項の計算（ナビエストークス方程式と非圧縮条件）

　ナビエストークス方程式をMPSの方法で離散化し、流速・圧力を計算する方法を説明する。流体の運動を記述するナビエストークス方程式は、質量保存に関する式と共に

\begin{align}
&\frac{Du}{Dt} = -\frac{1}{\rho}\nabla P + \nu \nabla^2 u + g \\
&\frac{D\rho}{Dt} + \rho\nabla \cdot u = 0
\end{align}

と書ける。$u$ は流速、$P$ は圧力、$\rho$ は密度、$g$ は重力加速度である。非圧縮流体の場合、密度は時間に依らないため、

\begin{align}
\frac{D\rho}{Dt} = 0 \ \  \Rightarrow  \ \ \nabla \cdot u = 0
\end{align}

となり $\nabla \cdot u = 0$ は非圧縮条件である。非圧縮条件からポアソン方程式を導き出し圧力を計算する。以下では、projection法と呼ばれる半陰解法型アルゴリズムを説明する。

圧力項の計算（ポアソン方程式の導出）

　ナビエストークス方程式を離散化すると、$k+1$ ステップの流速は、$k$ ステップの流速と圧力を使い、

\begin{align}
u_i^{k+1} &= u_i^{k}-\frac{1}{\rho_0} \langle \nabla P \rangle_i^{k+1} \Delta t + \nu \langle \nabla^2 u  \rangle_i^k \Delta t + g \Delta t \\
\end{align}

と書ける。右辺の１項目は圧力項、２項目は粘性項であり、最後の項は外力項（重力項）である。時間微分は、

\begin{align}
\frac{Du_i}{Dt} \approx \frac{u_i^{k+1}-u_i^{k}}{\Delta t} 
\end{align}

と離散化した。上記の式を（添え字 $i$ を省いて）

\begin{align}
u^{k+1} & = u^{k} + \delta u^p + \delta u^c \\
\delta u^p & = \nu \langle \nabla^2 u  \rangle^k \Delta t + g \Delta t  \\
\delta u^c & = -\frac{1}{\rho_0} \langle \nabla P \rangle^{k+1} \Delta t
\end{align}

と分ける。$\delta u^p$ は求まるが、圧力 $P$ は未知量なので$\delta u^c$ は現時点では求まらない。$\delta u^c$ を除いた仮速度・仮位置を

\begin{align}
u^* &= u^k + \delta u^p \\
r^* &= r^k + u^* \Delta t \\
\end{align}

とする。
　
　圧力の計算は、質量保存に関する条件より計算される。仮速度により更新された仮の位置$r^{\ast}$ を使い、数密度を計算する。それを $n^{\ast}$ とすると、質量は保存されること、つまり数密度は一定であるので、$k+1$ ステップの数密度 $n^{k+1}$ は、

\begin{align}
n^{k+1} = n^{\ast} + n^c = n_0
\end{align}

となる。$n^c$は、$\delta u^c$ によって計算された数密度である。また密度と数密度は比例関係にあるので、

\begin{align}
\frac{\rho_0 -\rho^*}{\rho_0} = \frac{n_0-n^*}{n_0} 
\end{align}

が成立する。

　$\delta u^c$ における質量保存則は、

\begin{align}
\frac{1}{\rho_0}\left(\frac{D\rho}{Dt} \right)^c+ \nabla\cdot \delta u^c = 0
\end{align}

または、

\begin{align}
\frac{1}{n_0}\left(\frac{Dn}{Dt} \right)^c+ \nabla\cdot \delta u^c = 0
\end{align}

とかけ、この式を離散化すると、$\delta u^c= u^{k+1} - u^k$ であり、$k+1$ ステップの流速 $u^{k+1}$ は、非圧縮条件 $\nabla\cdot u^{k+1}=0$ が成立するので、

\begin{align}
&\frac{1}{\rho_0}\frac{\rho_0 -\rho^*}{\Delta t}+ \nabla\cdot (u^{k+1} - u^*)= 0 \\
\end{align}

または、

\begin{align}
&\frac{1}{n_0}\frac{n_0 -n^*}{\Delta t}+ \nabla\cdot (u^{k+1} - u^*)= 0 \\
\end{align}

より

\begin{align}
 \nabla\cdot u^* = -\frac{1}{\Delta t}\frac{\rho^*-\rho_0 }{\rho_0}
= -\frac{1}{\Delta t} \frac{n^*-n_0}{n_0} 
\end{align}

となる。さらに、

\begin{align}
\nabla\cdot\delta u^c & =\nabla\cdot (u^{k+1} - u^*)=-\nabla\cdot u^*= -\nabla\cdot \left(\frac{1}{\rho} \langle \nabla P \rangle^{k+1}\right) \Delta t
\end{align}

であるので、MPSにおける（密度が一様の場合の）ポアソン方程式は、

\begin{align}
\langle \nabla^2 P \rangle^{k+1} = -\frac{\rho_0}{\Delta t^2} \frac{n^*-n_0}{n_0} 
\end{align}

と求まる。

圧力項の計算（連立方程式）

　ポアソン方程式の左辺は、ラプラシアンの計算なので

\begin{align}
\langle \nabla^2 P \rangle^{k+1}_i
& = \frac{2D}{\lambda n_0} \sum_{j \neq i} \omega(|r_{ij}^*|) \left(P^{k+1}_j - P^{k+1}_i \right)\\
& = -\frac{2D}{\lambda n_0}P^{k+1}_i \sum_{j \neq i} \omega(|r_{ij}^*|) 
+ \frac{2D}{\lambda n_0} \sum_{j \neq i} \omega(|r_{ij}^*|) P^{k+1}_j
\end{align}

と離散化される。ポアソン方程式は、連立方程式 $Ax=b$ に帰着でき

\begin{align}
&Ax=b　\\
&x_i = P_i^{k+1} , \ \  b_i = \frac{\rho_0}{\Delta t^2} \frac{n^*_i-n_0}{n_0} 
\end{align}

\begin{align}
A_{ij} = 
\begin{cases}
\frac{2D}{\lambda n_0}\sum_{j \neq i} \omega(|r_{ij}^*|)  & (i=j)  \\
-\frac{2D}{\lambda n_0} \omega(|r_{ij}^*|)             & (i\neq j)
\end{cases}
\end{align}

となる。以上の連立方程式を解けば、圧力 $P_i^{k+1}$ が得られる。

　連立方程式を解く際、dirichlet 境界条件を考慮する必要がある。表面における粒子は、（流体の外は粒子がいないので）圧力がゼロ（または大気圧）になる。したがって、粒子が、自由表面粒子かどうか判断する必要がある。

　表面の粒子は、流体内部の粒子と比べ、数密度が少なくなるので（流体の外に粒子はいないから）、仮の位置$r^{\ast}$ から得られた数密度 $n^{\ast}$ が

\begin{align}
n^{\ast}_i < \beta n_0 \ \ , \ \ \beta = 0.95 \sim 0.97
\end{align}

の条件を満たす場合、自由表面粒子と判定し、$P_i^{k+1}=0$ とする。

圧力項の計算（圧力勾配）

　連立方程式を解くことで、圧力 $P_i^{k+1}$ が求まったので、その値を使い勾配を計算する。ただし、計算を安定させるため、${P}_i \to \hat{P}_i $ と置き換え

\begin{align}
\langle \nabla P \rangle_{i} ^{k+1}
&= \frac{D}{n_0} \sum_{j \neq i} \omega(|r_{ij}|) \frac{P_j^{k+1}-\hat{P}_i^{k+1}}{|r_{ij}|^2} r_{ij}
\end{align}

とする。 $\hat{P}_i $ は、

\begin{align}
\hat{P}_i &= \min_{j\in \Omega} (P_i,P_j) \\
\Omega &= \{j:\omega(|r_{ij}|) \neq 0\}
\end{align}

であり、半径 $r_e$ 内に存在する粒子の中で最小となる圧力である。
　${P}_i \to \hat{P}_i $ と置き換えることで、圧力差は非負

\begin{align}
P_j^{k+1}-\hat{P}_i^{k+1} \geq 0
\end{align}

であることが保証される。これは斥力（引力の逆、粒子同士を遠ざける力）が働くことを意味し、粒子が近づきすぎることによる計算の不安定性を回避することができる。

まとめ

　離散化されたナビエストークス方程式は

\begin{align}
u_i^{k+1} &= u_i^{k}-\frac{1}{\rho_0} \langle \nabla P \rangle_i^{k+1} \Delta t + \nu \langle \nabla^2 u  \rangle_i^k \Delta t + g \Delta t \\
\end{align}

であり、粘性項は

\begin{align}
\langle \nabla^2 u \rangle_{i} ^k 
&= \frac{2D}{\lambda n_0} \sum_{j \neq i} \omega(|r_{ij}|) (u_j^k  - u_i^k ) 
\end{align}

と離散化され、$P_i^{k+1}$ の値をポアソン方程式から計算し

\begin{align}
&Ax=b　\\
&x_i = P_i^{k+1} , \ \  b_i = \frac{\rho_0}{\Delta t^2} \frac{n^*_i-n_0}{n_0} 
\end{align}

\begin{align}
A_{ij} = 
\begin{cases}
\frac{2D}{\lambda n_0}\sum_{j \neq i} \omega(|r_{ij}^*|)  & (i=j)  \\
-\frac{2D}{\lambda n_0} \omega(|r_{ij}^*|)             & (i\neq j)
\end{cases}
\end{align}

圧力項は、ポアソン方程式から得られた $P_i^{k+1}$ の値を使い、

\begin{align}
\langle \nabla P \rangle_{i} ^{k+1}
&= \frac{D}{n_0} \sum_{j \neq i} \omega(|r_{ij}|) \frac{P_j^{k+1}-\hat{P}_i^{k+1}}{|r_{ij}|^2} r_{ij}
\end{align}

と離散化され、$k+1$ ステップの流速が求まる。

　以上が、MPSの説明である。MPSの離散化では、質量保存則とは別に運動量保存則や角運動量保存則を満たす必要がある。

　次回の記事では、それらの保存則について、または、高精度化について紹介する。