【粒子法 MPS法】アルゴリズム高精度化まとめ（HS法・HL法・ECS法・GC法・DS法・SPP法）

下記の書籍に基づいてMPS法のアルゴリズム高精度化を進めてきて、ようやく一通り完成したのでまとめるよ。この書籍にはプログラミング時の実装方法について具体的な言及がないので、僕の理解不足もあるのだけれども解釈に不明確な点があったね。とくにHL法の導出過程とGC法との整合性、そもそもGC法をどこにどのように適用するのかがわからない、SPP法の適用条件が不明確だったので、今回ためした実装方法をメモがてら残しておくね。

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

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

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

１．仮時刻の速度計算のための速度のラプラシアンモデル

MPS法序盤の仮時刻における速度計算に必要な速度ラプラシアンのモデル化だよ。

$$
\boldsymbol{v}^{(*)}_i = \boldsymbol{v}^{(k)}_i + \left[\nu \left\langle \nabla^2\boldsymbol{v}\, \right\rangle^{(k)}_i + \boldsymbol{g} \right] \Delta t
$$

標準MPS法

$$
\left\langle \nabla^2\boldsymbol{v}\, \right\rangle^{(k)}_i = \frac{2 {\rm dim}}{ \lambda^{(0)} n^{(0)} } \sum\limits_{j=1,\, j\ne i}^N \boldsymbol{v}^{(k)}_{ij} \, \omega(|\boldsymbol{r}^{(k)}_{ij} |)
$$

HS法

$$
\left\langle \nabla^2\boldsymbol{v}\, \right\rangle^{(k)}_i = \frac{\nu(5-{\rm dim})}{n^{(0)}}\sum\limits_{j\ne i}\frac{ r_e }{|\boldsymbol{r}^{(k)}_{ij} |^3}\, \boldsymbol{v}_{ij}^{(k)}
$$

GC法

$$
\left\langle \nabla^2\boldsymbol{v}\, \right\rangle^{(k)}_i = {\rm tr}\left( \boldsymbol{B}^{(k)}_i\right)\left[ \frac{\nu}{n^{(0)}}\sum\limits_{j\ne i}\frac{\boldsymbol{v}^{(k)}_{ij} }{|\boldsymbol{r}^{(k)}_{ij} |^2}\,\omega(|\boldsymbol{r}^{(k)}_{ij} |) \right]
$$

$$
\boldsymbol{B}^{(k)}_i = \left[\frac{1}{n^{(0)}}\sum\limits_{j\ne i} \left(\frac {\boldsymbol r^{(k)}_{ij}\otimes\boldsymbol r^{(k)}_{ij}}{|\boldsymbol r^{(k)}_{ij}|^2}\right) \omega(|\boldsymbol r^{(k)}_{ij}|) \right]^{-1}
$$

２．圧力計算のためのポアソン方程式右辺の速度の発散モデル

次はMPS法の要である圧力計算のためのポアソン方程式の右辺、速度の発散のモデル化だよ。書籍ではHS法とGC法は独立に実装可能との記述があるのだけれども、ここのHS法をそのままとして、他のGC法を適用すると計算不安定となってしまったよ。

$$
\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法

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

HS法

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

GC法

$$
\langle\nabla\cdot \boldsymbol{v}\rangle^{(*)}_i=\frac{1}{n^{(0)}} \sum\limits_{j\ne i}\,\frac{ \omega(|\boldsymbol{r}^{(*)}_{ij} |) }{|\boldsymbol{r}^{(*)}_{ij} |^2} \, {\rm tr}\!\!\left[ \left( \boldsymbol{v}^{(*)}_{ij} \otimes \boldsymbol{r}^{(*)}_{ij} \right) \cdot \boldsymbol{B}^{(*)}_i\right]
$$

ECS法

ECS法は独立ではなく、標準MPS法・HS法・GC法に加える形で実装されるね。

$$
\langle\nabla\cdot \boldsymbol{v}\rangle^{(*)}_i +\!\!= -\frac{\gamma}{\Delta t}\,\frac{n_i^{(k)}-n^{(0)}}{ n^{(0)}} - \frac{\beta}{n^{(0)}} \sum\limits_{j\ne i} \frac{r_e}{|\boldsymbol{r}^{(k)}_{ij} |^3}\, \boldsymbol{r}^{(k)}_{ij}\cdot \boldsymbol{v}^{(k)}_{ij}
$$

３．圧力計算のためのポアソン方程式左辺の圧力のラプラシアンモデル

次はMPS法の要である圧力計算のためのポアソン方程式の左辺、圧力ラプラシアンのモデル化だよ。書籍ではHL法とSPP法は独立である記述があるのだけれども、HL法には$${\omega}$$が存在しないので、そもそもSPP法の導出と整合がつくのかわからなかったね。

$$
\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法

$$
\langle \nabla^2 P\, \rangle_i^{(k+1)} =\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}^{(*)}_{ij}|)
$$

HL法

$$
\langle \nabla^2 P\,\rangle_i^{(k+1)} =\frac{5-{\rm dim}}{n^{(0)}}\sum\limits_{j=1,\, j\ne i}^N (P^{(k+1)}_j-P^{(k+1)}_i)\,\frac{ r_e }{|\boldsymbol{r}^{(*)}_{j}-\boldsymbol{r}^{(*)}_{i} |^3}
$$

GC法

$$
\langle \nabla^2 P\,\rangle_i^{(k+1)} = \frac{{\rm tr}\left( \boldsymbol{B}^{(*)}_i\right)}{ n^{(0)}}\left[ \sum\limits_{j=1,\, j\ne i}^N(P^{(k+1)}_j-P^{(k+1)}_i)\,\frac{ \omega(|\boldsymbol{r}^{(*)}_{ij}|) }{|\boldsymbol{r}^{(*)}_{ij} |^2} \right]
$$

標準MPS法＋SPP法

$$
\langle \nabla^2 P\, \rangle_i^{(k+1)} =\frac{2{\rm dim}}{\lambda^{(0)}n^{(0)}}\left[ \sum\limits_{j=1,\, j\ne i, i{\rm SPP}}^N P^{(k+1)}_j\, \omega(|\boldsymbol{r}^{(*)}_{ij}|)-P^{(k+1)}_i \left\{\sum\limits_{j=1,\, j\ne i, i{\rm SPP}}^N \omega(|\boldsymbol{r}^{(*)}_{ij}|) + \omega_{i{\rm SPP}}\right\}\right]
$$

GC法＋SPP法

$$
\langle \nabla^2 P\, \rangle_i^{(k+1)} =\frac{{\rm tr}\left( \boldsymbol{B}^{(*)}_i\right)}{ n^{(0)}} \left[ \sum\limits_{j=1,\, j\ne i, i{\rm SPP}}^N P^{(k+1)}_j\,\frac{ \omega(|\boldsymbol{r}^{(*)}_{ij}|) }{|\boldsymbol{r}^{(*)}_{ij} |^2} -P^{(k+1)}_i \left\{\sum\limits_{j=1,\, j\ne i, i{\rm SPP}}^N \frac{\omega(|\boldsymbol{r}^{(*)}_{ij}|)}{|\boldsymbol{r}^{(*)}_{ij}|^2} + \frac{\omega_{i{\rm SPP}}}{R^2_{ i{\rm SPP}}}\right\}\right]
$$

$$
R_{i{\rm SPP}}= \frac{r_e}{ n^{(0)} -n_i^{(*)} + 1 } \ , \ \omega_{i{\rm SPP}} = \frac{r_e}{ R_{i{\rm SPP}}} -1
$$

SPP法は完全に孤立した粒子の場合は$${\omega_{i{\rm SPP}}=0}$$となるよ。

４. 速度更新時の圧力のグラディエントモデル

次は速度更新時の圧力グラディエントのモデル化だよ。$${\hat{P}^{(k+1)}_i}$$は$${i}$$番目の粒子自身を含む周囲の圧力の最小値だね。なお、次のDS法を実装する場合はこの速度が最終的な次時刻の速度とはならないよ。

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

標準MPS法

$$
\langle \nabla P\, \rangle_i^{(k+1)} =\frac{{\rm dim}}{n^{(0)}} \sum\limits_{j=1,\, j\ne i }^N \frac{ P^{(k+1)}_j-\hat{P}^{(k+1)}_i}{|\boldsymbol{r}^{(*)}_{ij} |^2}\,\boldsymbol{r}^{(*)}_{ij}\, \omega(|\boldsymbol{r}^{(*)}_{ij}|)
$$

GC法

$$
\langle \nabla P\, \rangle_i^{(k+1)} =\frac{{\rm dim}}{n^{(0)}} \sum\limits_{j=1,\, j\ne i}^N \frac{ P^{(k+1)}_j-\hat{P}^{(k+1)}_i}{|\boldsymbol{r}^{(*)}_{ij}|^2} \,\boldsymbol{C}^{(*)}_{i} \cdot \boldsymbol{r}^{(*)}_{ij}\, \omega(|\boldsymbol{r}^{(*)}_{ij}|)
$$

$$
\boldsymbol{C}^{(*)}_{i}= \frac{1}{{\rm dim}}\left[ \frac{1}{n^{(0)}}\sum\limits_{j=1,\, j\ne i}^N\frac {\omega(| \boldsymbol{r}^{(*)}_{ij}|)}{| \boldsymbol{r}^{(*)}_{ij}|^2}\left( \boldsymbol r^{(*)}_{ij}\otimes\boldsymbol r^{(*)}_{ij} \right) \right]^{-1}
$$

SPP法

ディリクレ境界条件として$${ P^{(k+1)}_{i{\rm SPP}}}=0$$を考慮しているよ。

$$
\langle \nabla P\, \rangle_i^{(k+1)} +\!\!\!= \frac{{\rm dim}}{n^{(0)}}\left[\frac{-\hat{P}^{(k+1)}_i}{R_{{i{\rm SPP}}}^2}\,\boldsymbol{r}^{(*)}_{{ii{\rm SPP}}}\,\omega_{i{\rm SPP}} \right]
$$

$$
\boldsymbol{r}^{(*)}_{ii{\rm SPP}} = -\frac{ R_{i{\rm SPP}} }{|\boldsymbol{r}^{(*)}_{ig}|}\,\boldsymbol{r}^{(*)}_{ig} \ , \ \boldsymbol{r}^{(*)}_{ig} = \frac{1}{n^{(0)}} \sum\limits_{j=1, j\ne i}^N \boldsymbol{r}^{(*)}_{ij}\, \omega (|\boldsymbol{r}^{(*)}_{ij}|)
$$

ちなみに$${|\boldsymbol{r}^{(*)}_{ig}|>0}$$であれば近傍に粒子が存在することになるので、$${|\boldsymbol{r}^{(*)}_{ig}|=0}$$の場合に$${\omega_{i{\rm SPP}}=0}$$を与えれば良いね。

５. 人工斥力の導入（DS法）

最後はMPS法の宿命である粒子同士の異常接近を阻止するために導入する人工斥力だよ。DS法を導入する場合は、$${\hat{P}^{(k+1)}_i} \to P^{(k+1)}_i}$$と置き換え、先の次時刻の速度と位置を仮のものと考えて、特定の条件を満たしたときに人工斥力$${\boldsymbol{ F}_{i}^{\, ({\rm DS)}}}$$を粒子に加えるよ。

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

$$
\boldsymbol{v}^{(k+1)}_i = \boldsymbol{v}^{(k+1)^*}_i -\frac{1}{\rho_0}\, \boldsymbol{ F}_{i}^{\, ({\rm DS)}} \Delta t\\
\boldsymbol{r}^{(k+1)}_i = \boldsymbol{r}^{(k+1)^*}_i -\frac{1}{\rho_0}\, \boldsymbol{ F}_{i}^{\, ({\rm DS)}} \Delta t^2
$$

$$
\boldsymbol{ F}_{i}^{\, ({\rm DS)}} = \frac{\rho_i }{n^{(0)}} \sum\limits_{j=1,\, j\ne i}^N \boldsymbol{F}_{ij}^{\,(\rm DS)} \omega(|\boldsymbol{r}^{(k+1)^*}_{ij}|)
$$

人工斥力の発生条件と$${\boldsymbol{F}_{ij}^{\,(\rm DS)} }$$の具体的な表式は以下のページを参考にしてください。

実行結果

Python プログラムソース

標準MPS法・HS法・ECS法・DS法・GC法・SPP法の導入を切り替えることの機能をつけたPythonプログラムソースを販売するよ。もしよかったら購入してみてください。