概要
土木の分野において、粒子法の研究が盛んにおこなわれています。前回の記事で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は、それら両方を満たし保存則が成立する。