概要

　土木の分野において、粒子法の研究が盛んにおこなわれています。前回の記事でMPSについて説明しました。

　今回の記事では、保存則を満たすように圧力項と修正する方法（cMPS）、ポアソン方程式の高精度生成項（HS）、高精度ラプラシアン（HL）について紹介したいと思います。

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

保存則

　質点系の全運動量 $G$ と全角運動量 $H$ は、

\begin{align}
 G = \sum_i m_iu_i \ \ , \ \ H = \sum_i r_i \times m_i u_i
\end{align}

とかける。$m_i$ は粒子 $i$ の質量、$u_i$ は粒子 $i$ の速度、$r_i$ は粒子 $i$ の半径である。

　ここでは、外力はゼロとして、粒子 $i$ の内力を $T_i$

\begin{align}
T_i = \sum_j T_{ij}
\end{align}

とすれば、ニュートンの第２法則より

\begin{align}
m_ia_i=m_i\frac{du_i}{dt} = -T_i
\end{align}

と書ける。全運動量または全角運動量が保存するとは、それらが時間変化しないこと、つまり、微分してゼロになる必要がある。

\begin{align}
\frac{dG}{dt} &=  \sum_i m_i \frac{du_i}{dt} = -\sum_i T_i = 0 \\
\frac{dH}{dt} &=  \sum_i r_i \times m_i \frac{du_i}{dt} = -\sum_i r_i \times T_i = 0 \\
\end{align}

つまり、

\begin{align}
&\sum_i T_i =\sum_i \sum_j T_{ij} = 0 \\
&\Rightarrow T_{ij}+T_{ji} = 0 
\end{align}

\begin{align}
&\sum_i r_i \times T_i =\sum_i \sum_j r_i \times T_{ij} = 0 \\
&\Rightarrow r_i \times T_{ij}+ r_j \times T_{ji} = -(r_j-r_i)\times T_{ij}=
-r_{ij}\times T_{ij}= 0 
\end{align}

を満たす必要がある。

　これを満たすには、$T_{ij}$ が $i$ と $j$ に関して反対称であること

\begin{align}
T_{ij} = -T_{ji} 
\end{align}

$T_{ij}$ が $r_{ij}$ に比例すること（同じベクトル同士の外積はゼロなので）

\begin{align}
T_{ij}  \propto r_{ij}
\end{align}

である。

MPS の圧力項について

　MPSの圧力項は、内力 $T_{ij}$ で書くと

\begin{align}
T_{ij}= - \frac{D}{n_0} \omega(|r_{ij}|) \frac{P_j^{k+1}-\hat{P}_i^{k+1}}{|r_{ij}|^2} r_{ij}
\end{align}

と書ける。$T_{ij}$ が $r_{ij}$ に比例しているが、

\begin{align}
P_j^{k+1}-\hat{P}_i^{k+1} \neq -( P_i^{k+1}-\hat{P}_j^{k+1})
\end{align}

であるので反対称でない。

　今回の記事では、粘性項が保存則を満たすようには修正はしないが（修正の仕方を勉強していない）、粘性項も内力 $T_{ij}$ で書くと

\begin{align}
T_{ij}= - \frac{2D}{\lambda n_0}  \omega(|r_{ij}|) (u_j^k  - u_i^k ) 
\end{align}

と書ける。、$T_{ij}$ が $i$ と $j$ に関して反対称であるが、$T_{ij}$ が $r_{ij}$ に比例していないため保存則を満たさない。

　この記事では、保存則を満たすように圧力項を修正する。この方法をcMPSと呼ばれる。

cMPS

　cMPSについて説明する。cMPSは、粒子 $i,j$ の間に仮想粒子 $k$ を考えることにより、MPSで満たされていなかった反対称性を満たすように修正を行う。

　粒子 $i,j$ の間にある仮想粒子 $k$ は、

\begin{align}
P_k = \frac{P_i+P_j}{2}  \ \ , \ \ r_{ik}=r_k-r_i = \frac{r_{ij}}{2}
\end{align}

と計算できるとする。重み関数は、

\begin{align}
\omega(|r_{ik}|) = \frac{r_{eik}}{r_{ik}}-1=  \frac{r_{eij}/2}{r_{ij}/2}-1
=\omega(|r_{ij}|)
\end{align}

となるので、数密度は

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

となり添え字を $k$ を $j$ に変えても、変化しない。

　MPSの圧力勾配を、添え字 $k$ で足し上げて、$P_k , r_{ik}$ を代入し、

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

圧力の最小値

\begin{align}
\hat{P}_i \rightarrow \frac{\hat{P}_i +\hat{P}_j}{2}
\end{align}

とすれば、

\begin{align}
\langle \nabla P \rangle_{i}
&= \frac{D}{n_0} \sum_{j \neq i} \omega(|r_{ij}|) \frac{\left(P_i+P_j \right)-\left(\hat{P}_i +\hat{P}_j\right) }{|r_{ij}|}\frac{r_{ij}}{|r_{ij}|} \\
\end{align}

圧力勾配が求まる。 $i$ と $j$ に関して反対称であり、$T_{ij}$ が $r_{ij}$ に比例しているので保存則を満たしている。

ポアソン方程式の復習

　上記でcMPSについて説明したが、圧力勾配の精度は、圧力 $P$ を求める精度、つまりポアソン方程式の解の精度に依ってくる。

　ポアソン方程式を再度導出すると、$\delta u^c$ における質量保存則は、

\begin{align}
\frac{1}{n_0}\left(\frac{Dn}{Dt} \right)^c+ \nabla\cdot \delta u^c = 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}

が成立するため、ポアソン方程式は

\begin{align}
\nabla\cdot \left(\frac{1}{\rho} \langle \nabla P \rangle^{k+1}\right)  = \frac{1}{n_0 \Delta t} \left(\frac{Dn}{Dt} \right)^c
\end{align}

と書ける。ポアソン方程式の右辺の高精度化の方法を高精度生成項と呼ばれ、左辺の高精度化の方法を高精度ラプラシアンと呼ばれる。

　この記事では、高精度ラプラシアンの方法を粘性項の計算でも使う。

ポアソン方程式の高精度生成項（HS）

　ここでは、ポアソン方程式の高精度生成項、連立方程式の右辺の項の高精度化について説明する。

　数密度の微分は、重み関数の微分

\begin{align}
\frac{d}{d|r|}\omega(|r|)  = - \frac{r_e}{|r|^2}
\end{align}

と $r$ に関する微分

\begin{align}
\frac{d}{dr}|r|  = - \frac{r}{|r|}　\ \ , \ \ \frac{d}{dt}r_{ij}  = u_j-u_i \equiv u_{ij}
\end{align}

を使い、

\begin{align}
\frac{Dn_i}{Dt} &= \sum_{j \neq i} \frac{D}{Dt} \omega(|r_{ij}|) \\
& = \sum_{j \neq i} \frac{dr_{ij} }{dt} \frac{d|r_{ij}|}{dr_{ij}} \frac{d\omega(|r_{ij}|) }{d|r_{ij}|} \\
&=- \frac{r_e}{|r_{ij}|^3} r_{ij}\cdot u_{ij}
\end{align}

と求まる。（本来は式の途中でδ関数の和を取る。）

　質量保存から

\begin{align}
\frac{Dn}{Dt} &=\frac{D}{Dt}(n^p+n^c)= 0
\end{align}

であるので、

\begin{align}
\left(\frac{Dn}{Dt}\right)^c &= - \left(\frac{Dn}{Dt}\right)^p\\
& = - \sum_{j \neq i} \left(\frac{D\omega(|r^*_{ij}|) }{Dt}\right)^p \\
& = \sum_{j \neq i} \frac{r_e}{|r_{ij}^*|^3} r_{ij}^*\cdot u_{ij}^*
\end{align}

ともとまる。最終的に、ポアソン方程式の右辺は

\begin{align}
\nabla\cdot \left(\frac{1}{\rho} \langle \nabla P \rangle^{k+1}\right)  = \frac{1}{n_0 \Delta t} \sum_{j \neq i} \frac{r_e}{|r_{ij}^*|^3} r_{ij}^*\cdot u_{ij}^*
\end{align}

と書ける。

高精度ラプラシアン（HL）

　ここでは、ラプラシアンの高精度化、ポアソン方程式の左辺の項の高精度化について説明する。勾配モデルを

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

と書くことにする。（このモデルはSPHの勾配モデルであり、単純には本来のMPSの勾配モデルと異なる。）

　ラプラシアンは、

\begin{align}
\langle \nabla^2 \phi \rangle_i &= \nabla\cdot \langle \nabla \phi \rangle_i  \\
&= \frac{1}{n_0}\sum_{j \neq i} \left(\nabla\phi_{ij} \cdot \nabla\omega(|r_{ij}|) +\phi_{ij} \nabla^2 \omega(|r_{ij}|) \right)
\end{align}

と書ける。連鎖律を用いれば

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

となるので、ラプラシアンは、

\begin{align}
\langle \nabla^2 \phi \rangle_i 
&= \frac{1}{n_0}\sum_{j \neq i} \left\{\frac{\partial \phi_{ij}}{\partial |r_{ij}|} \frac{\partial \omega(|r_{ij}|)}{\partial |r_{ij}|} + \phi_{ij}\left(\frac{\partial^2 \omega(|r_{ij}|)}{\partial |r_{ij}|^2} + \frac{D-1}{|r_{ij}|}\frac{\partial \omega(|r_{ij}|)}{\partial |r_{ij}|} \right) \right\}
\end{align}

と書ける。

密度が非一様な高精度ラプラシアン（HL）

　今までの説明は、密度が一様であることを前提であったが、固体と液体が混合状態で存在するとき、領域内の密度の変化も考慮する必要がある。
　そのとき、ラプラシアンは、

\begin{align}
\nabla\cdot \left(\frac{\nabla \phi}{\rho} \right)  
&= \frac{1}{\rho}\nabla\cdot \left( \nabla \phi \right) +  \nabla\left(\frac{1}{\rho} \right) \cdot  \nabla \phi
\end{align}

と書ける。１項目は、上記で説明した高精度ラプラシアン（HL）の結果を使う。２項目の計算は、連鎖律を用いて

\begin{align}
\nabla\left(\frac{1}{\rho} \right) &= \frac{\partial }{\partial |r|}\left(\frac{1}{\rho} \right) \nabla |r|  = \frac{\partial }{\partial |r|}\left(\frac{1}{\rho} \right) \frac{r}{|r|} \\
\nabla\omega(|r|) &= \frac{\partial \omega(|r|)}{\partial |r|} \frac{r}{|r|} 
\end{align}

粒子$i,j$ において、密度に作用する微分の離散化を

\begin{align}
\left\{ \frac{\partial }{\partial |r|}\left(\frac{1}{\rho} \right) \right\}_{ij} &= 
\frac{1}{r_{ij}} \left\{\frac{1}{\rho_j} - \frac{1}{\rho_i} \right\} \\
\end{align}

とすれば

\begin{align}
\left\langle \left\{\nabla\left(\frac{1}{\rho} \right) \right\} \cdot \nabla \phi  \right\rangle_i &= \frac{1}{n_0}\sum_{j \neq i} \phi_{ij}\left\{\nabla\left(\frac{1}{\rho} \right) \right\}_i \cdot \nabla\omega(|r_{ij}|) \\
&=  \frac{1}{n_0}\sum_{j \neq i} \phi_{ij}\frac{1}{r_{ij}} \left\{\frac{1}{\rho_j} - \frac{1}{\rho_i} \right\} \frac{\partial \omega(|r_{ij}|)}{\partial |r_{ij}|}\\
\end{align}

となるので、密度が非一様なラプラシアンは、

\begin{align}
&\left\langle \nabla\cdot \left(\frac{\nabla \phi}{\rho} \right) \right\rangle_i = \\
& \frac{1}{n_0}\sum_{j \neq i} \left\{ \frac{1}{\rho_i} \frac{\partial \phi_{ij}}{\partial |r_{ij}|} \frac{\partial \omega(|r_{ij}|)}{\partial |r_{ij}|} +\frac{\phi_{ij}}{\rho_i}\frac{\partial^2 \omega(|r_{ij}|)}{\partial |r_{ij}|^2} + \left(\frac{1}{\rho_j} + \frac{D-2}{\rho_i} \right) \frac{\phi_{ij}}{|r_{ij}|} \frac{\partial \omega(|r_{ij}|)}{\partial |r_{ij}|}\right\}
\end{align}

と書ける。
　次に、具体的にMPSの重み関数を代入すると

\begin{align}
&\frac{d}{d|r|}\omega(|r|) = - \frac{r_e}{|r|^2} \\
&\frac{d^2}{d|r|^2}\omega(|r|) = 2\frac{r_e}{|r|^3} \\
&\frac{d\phi_{ij}}{d|r_{ij}|}= -2\frac{\phi_{ij}}{|r_{ij}|} 
\end{align}

なので、

\begin{align}
&\left\langle \nabla\cdot \left(\frac{\nabla \phi}{\rho} \right) \right\rangle_i = \\
& \frac{1}{n_0}\sum_{j \neq i} \left\{ \frac{1}{\rho_i} \left(\frac{-2\phi_{ij}}{|r_{ij}|} \right) \left(\frac{-r_e}{|r_{ij}|} \right)
+\frac{\phi_{ij}}{\rho_i} \left(\frac{2r_e}{|r_{ij}|^3} \right)
+ \left(\frac{1}{\rho_j} + \frac{D-2}{\rho_i} \right) \frac{\phi_{ij}}{|r_{ij}|}  \left(\frac{-r_e}{|r_{ij}|} \right) \right\} \\
& = \frac{1}{n_0}\sum_{j \neq i}  \left(\frac{6-D}{\rho_i} -\frac{1}{\rho_j} \right) \frac{r_e}{|r_{ij}|^3}\phi_{ij}
\end{align}

と求まる。高精度ラプラシアンは、$i$ と $j$ に関して反対称であるが、$r_{ij}$ に比例していないため保存則を満たさない。

　一様な密度の場合、$\rho_i=\rho_j=\rho_0$ であるので、

\begin{align}
\frac{1}{\rho_0} \left\langle  \nabla^2 \phi \right\rangle_i = \frac{5-D}{\rho_0n_0}\sum_{j \neq i}  \frac{r_e}{|r_{ij}|^3}\phi_{ij} \\
\therefore  \left\langle  \nabla^2 \phi \right\rangle_i = \frac{5-D}{n_0}\sum_{j \neq i}  \frac{r_e}{|r_{ij}|^3}\phi_{ij}
\end{align}

と書ける。粘性項の計算ではこの式を使う。

ポアソン方程式のまとめ

　高精度生成項と高精度ラプラシアンを用いた場合、ポアソン方程式は

\begin{align}
\nabla\cdot \left(\frac{1}{\rho} \langle \nabla P \rangle^{k+1}\right) & = \frac{1}{n_0 \Delta t} \left(\frac{Dn}{Dt} \right)^c　\\
\end{align}

\begin{align}
\therefore \ \  \frac{1}{n_0}\sum_{j \neq i}  \left(\frac{6-D}{\rho_i} -\frac{1}{\rho_j} \right) \frac{r_e}{|r_{ij}^*|^3}(P_j^{k+1} - P_i^{k+1}) &= \frac{1}{n_0 \Delta t} \sum_{j \neq i} \frac{r_e}{|r_{ij}^*|^3} r_{ij}^*\cdot u_{ij}^*
\end{align}

となる。

　連立方程式$Ax=b$ の形にすると

\begin{align}
&Ax=b　\\
&x_i = P_i^{k+1} , \ \  b_i = \frac{1}{n_0 \Delta t} \sum_{j \neq i} \frac{r_e}{|r_{ij}^*|^3} r_{ij}^*\cdot u_{ij}^*
\end{align}

\begin{align}
A_{ij} = 
\begin{cases}
 -\frac{1}{n_0}\sum_{j \neq i}  \left(\frac{6-D}{\rho_i} -\frac{1}{\rho_j} \right) \frac{r_e}{|r_{ij}^*|^3}& (i=j)  \\
 \frac{1}{n_0}  \left(\frac{6-D}{\rho_i} -\frac{1}{\rho_j} \right) \frac{r_e}{|r_{ij}^*|^3} & (i\neq j)
\end{cases}
\end{align}

となる。密度が一様な場合は、対称行列になるが、非一様な場合は、非対称行列となり、連立方程式をBi-CGSTAB 法などの解法で解く必要がある。

MPS と cMPS-HS-HL の比較

　単純に、MPS と cMPS-HS-HL の比較を行う。

　以下は、MPS の damBreak の結果であり

以下は、cMPS-HS-HL の damBreak の結果である。

まとめ

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

\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}

である。

　以下、MPSとcMPS-HS-HLの差異をまとめる。

　粘性項は、MPSでは

\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}

であったが、高精度ラプラシアン（HL）を使えば

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

と離散化される。

　MPSでは、$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}

から計算したが、高精度生成項と高精度ラプラシアンを使うと（密度を一様として）

\begin{align}
&Ax=b　\\
&x_i = P_i^{k+1} , \ \  b_i = \frac{\rho_0}{n_0 \Delta t} \sum_{j \neq i} \frac{r_e}{|r_{ij}^*|^3} r_{ij}^*\cdot u_{ij}^*
\end{align}

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

から計算する。

　MPSにおいて、圧力項は、ポアソン方程式から得られた $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}

と離散化する。cMPSでは、

\begin{align}
\langle \nabla P \rangle_{i}
&= \frac{D}{n_0} \sum_{j \neq i} \omega(|r_{ij}|) \frac{\left(P_i+P_j \right)-\left(\hat{P}_i +\hat{P}_j\right) }{|r_{ij}|}\frac{r_{ij}}{|r_{ij}|} \\
\end{align}

と離散化する。MPSは、$i$ と $j$ に関して反対称であるが $r_{ij}$ に比例していないので保存則が成立しない。しかしCMPSは、それら両方を満たし保存則が成立する。