【粒子法 MPS】アルゴリズムの改善：HS法 ポアソン方程式生成項の高精度化

上記のnoteで解説した標準MPS法に加えて、下記の書籍で詳しく解説してある標準MPSの高精度化に着手していくよ。目次は次のとおりだよ。

書籍：粒子法:連続体・混相流・粒状体のための計算科学

P.83 HS法 ポアソン方程式生成項の高精度化
P.85 HL法 ラプラシアンモデルの高精度化
P.87 ECS法 ポアソン方程式生成項の誤差補償
P.92 GC法　勾配モデルの高精度化
P109 DS法　人工斥力導入による安定化
P118 SPP法　表面粒子の境界条件適正化
P124 WPP法　壁面における境界条件適正化

粒子法:連続体・混相流・粒状体のための計算科学

高精度化第一弾はHS法（ポアソン方程式生成項の高精度化）だよ。

HS法のアルゴリズム

HS法は各粒子に加わる圧力$${P}$$が満たすポアソン方程式

$$
\left\langle \nabla^2 P\, \right\rangle^{(k+1)}_i = \frac{\rho_0}{ \Delta t} \,\left\langle\nabla\cdot\boldsymbol{v}\right\rangle^{(*)}_i
$$

の右辺の高精度化を行う手法だよ。標準MPS法では仮時刻における連続の式

$$
\frac{ D\rho^{(*)} }{D t} +\rho^{(0)} \left\langle\nabla\cdot\boldsymbol{v}\right\rangle^{(*)}_i =0
$$

から

$$
\left\langle\nabla\cdot\boldsymbol{v}\right\rangle^{(*)}_i = -\frac{1}{\rho^{(0)} }\,\frac{ D\rho_i^{(*)} }{D t} \simeq -\frac{1}{\Delta t}\, \frac{\rho_i^{(*)}-\rho^{(0)}}{\rho^{(0)}} \simeq -\frac{1}{\Delta t}\,\frac{n_i^{(*)}-n^{(0)}}{ n^{(0)}}
$$

と密度の時間微分を１次で近似して、さらに粒子数密度で置き換える形で計算したね。HS法ではこの密度の時間微分を厳密に扱います。標準MPSと同様に仮時刻の密度を数密度で置き換えたあとに重み関数で表した

$$
\frac{\rho_i^{(*)}}{\rho^{0}} \simeq \frac{n_i^{(*)}}{n^{(0)}} = \frac{1}{n^{(0)}} \sum\limits_{j\ne i} \omega(|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|)
$$

を考えて、この時間微分を考えるよ。つまり、密度の時間微分を重み関数$${\omega(r)}$$の時間微分で表すわけだね。微分のチェーンルールを考慮すると

$$
\frac{ D \omega(r) }{D t} =\frac{\partial \omega(r)}{\partial r} \, \frac{\partial r}{\partial \boldsymbol{r}} \cdot \frac{d \boldsymbol{r} }{dt} =-\frac{r_e}{r^2}\,\frac{ \boldsymbol{r}}{ r } \cdot \boldsymbol{v} \ \ ,\ \ \frac{\partial r}{\partial \boldsymbol{r}}=\nabla r \ \ ,\ \ r=\sqrt{x^2+y^2+z^2}
$$

となるので、

$$
\frac{1}{\rho^{(0)} }\,\frac{ D\rho_i^{(*)} }{D t} \simeq\frac{1}{n^{(0)} }\,\frac{ Dn_i^{(*)} }{D t} =- \frac{1}{n^{(0)}} \sum\limits_{j\ne i} \frac{r_e}{|\boldsymbol{r}_j^{(*)}-\boldsymbol{r}^{(*)}_i|^3}\, (\boldsymbol{r}_j^{(*)}-\boldsymbol{r}^{(*)}_i) \cdot (\boldsymbol{v}_j^{(*)}-\boldsymbol{v}^{(*)}_i)
$$

となるので、速度の勾配は

$$
\left\langle\nabla\cdot\boldsymbol{v}\right\rangle^{(*)}_i =\frac{1}{n^{(0)}} \sum\limits_{j\ne i} \frac{r_e}{|\boldsymbol{r}_j^{(*)}-\boldsymbol{r}^{(*)}_i|^3}\, (\boldsymbol{r}_j^{(*)}-\boldsymbol{r}^{(*)}_i) \cdot (\boldsymbol{v}_j^{(*)}-\boldsymbol{v}^{(*)}_i)
$$

となるね。結果として圧力$${P}$$が満たすポアソン方程式は

$$
\left\langle \nabla^2 P\, \right\rangle^{(k+1)}_i = \rho_0\, \frac{1}{\Delta t}\,\frac{1}{n^{(0)}} \sum\limits_{j\ne i} \frac{r_e}{|\boldsymbol{r}_j^{(*)}-\boldsymbol{r}^{(*)}_i|^3}\, (\boldsymbol{r}_j^{(*)}-\boldsymbol{r}^{(*)}_i) \cdot (\boldsymbol{v}_j^{(*)}-\boldsymbol{v}^{(*)}_i)
$$

となるね。標準MPS法と合わせるために、両辺を$${-\rho^{(0)}}$$で割って圧力に関するラプラシアンモデルを適用すると$${i}$$番目の粒子の圧力に関する方程式が得られるよ。

$$
-\frac{1}{\rho^{(0)}}\,\frac{2{\rm dim}}{\lambda^{(0)}n^{(0)}} \sum\limits_{j=1,\, j\ne i}^N (P^{(k+1)}_j-P^{(k+1)}_i)\, \omega(|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|) \\
=- \frac{1}{\Delta t}\,\frac{1}{n^{(0)}} \sum\limits_{j\ne i} \frac{r_e}{|\boldsymbol{r}_j^{(*)}-\boldsymbol{r}^{(*)}_i|^3}\, (\boldsymbol{r}_j^{(*)}-\boldsymbol{r}^{(*)}_i) \cdot (\boldsymbol{v}_j^{(*)}-\boldsymbol{v}^{(*)}_i)
$$

右辺の式は複雑に見えるけれども、仮時刻の位置と速度はすでに計算済みなので、１回だけ和を計算するだけだね。ただし、注意は右辺の和を計算する際には、２粒子間の距離（$${|\boldsymbol{r}_j^{(*)}-\boldsymbol{r}^{(*)}_i|}$$）が影響半径$${r_e}$$よりも大きい場合はゼロになる点には注意が必要すべきだね。

HS法のまとめだけれども、ポアソン方程式の右辺を

$$
\frac{1}{(\Delta t)^2}\,\frac{n_i^{(*)}-n^{(0)}}{n^{(0)}} \to- \frac{1}{\Delta t}\,\frac{1}{n^{(0)}} \sum\limits_{j\ne i} \frac{r_e}{|\boldsymbol{r}_j^{(*)}-\boldsymbol{r}^{(*)}_i|^3}\, (\boldsymbol{r}_j^{(*)}-\boldsymbol{r}^{(*)}_i) \cdot (\boldsymbol{v}_j^{(*)}-\boldsymbol{v}^{(*)}_i)
$$

と置き換えるだけなので、実装は簡単だね。さっそく、水柱崩壊で効果を確かめてみよう！

計算結果比較

次の図は標準MPS法のみとHS法を加えたアルゴリズムの水柱崩壊後0.4秒後のスナップショットだよ。圧力によって色を変えているよ（赤：6000[Pa]、水色：0[Pa]）。

効果はてきめんだね。標準MPS法のままだと圧力がスパイク状に高い粒子が存在するけれども、HS法を加えると圧力スパイクは抑制されるだけでなく、圧力分布も自然な感じになっているね。

標準MPSのみ

標準MPS＋HS法

Pythonプログラム

以下は標準MPS法にHS法を加えたPythonプログラムソースを販売するよ。もし良かったら試してみてください。