【粒子法 MPS法】アルゴリズムの改善：SPP法 表面粒子の境界条件適正化
上記のnoteで解説したDS法に加えて、下記の書籍で詳しく解説してある標準MPSの高精度化に着手しているよ。目次は次のとおりだよ。
書籍：粒子法:連続体・混相流・粒状体のための計算科学
P.83 HS法 ポアソン方程式生成項の高精度化
P.85 HL法 ラプラシアンモデルの高精度化
P.87 ECS法 ポアソン方程式生成項の誤差補償
P.92 GC法 勾配モデルの高精度化
P109 DS法 人工斥力導入による安定化
P118 SPP法 表面粒子の境界条件適正化
P124 WPP法 壁面における境界条件適正化
書籍の順番に従って、表面粒子の境界条件適正化を目指すSPP法の導入を行うよ。ここまでのMPS法では圧力に関するポアソン方程式をディリクレ境界条件のもとで解く場合、表面粒子の圧力を０として計算していたね。この結果、隣り合う表面粒子は互いに力を及ぼせないため、物理的に不自然な挙動が生じてしまうことが知られているね。そこで、表面粒子間に相互作用を与えるために表面粒子の近傍に仮想的な粒子（SPP粒子）を配置して、このSPP粒子にディリクレ境界条件を課すことで、表面粒子にも圧力計算（ポアソン方程式）に加わってもらうよ。これはSPP法って呼ばれるよ。
上図は着目した$${i}$$番目の粒子が表面粒子の場合に配置するSPP粒子の模式図だよ。SPP粒子は着目する粒子が表面粒子の場合だけ１個だけ配置され、着目する粒子にのみ作用するものとするよ。そして、$${i}$$番目の粒子のSPP粒子の粒子数密度を$${ n^{(*)}_{i{\rm SPP}} }$$と表した場合、$${i}$$番目の粒子の粒子数密度$${ n^{(*)}_i }$$を用いて
$$
n^{(*)}_{i{\rm SPP}} = n^{(0)} - n^{(*)}_i
$$
で与えられると定義するね。時刻が仮時刻$$(*)$$なのは、圧力のポアソン方程式を計算する直前の数密度状態で境界条件を判定するためだよ。もともと粒子数密度の定義は
$$
n^{(*)}_i = \sum\limits_{j\ne i} \omega(|\boldsymbol{r}^{(*)}_{j} - \boldsymbol{r}^{(*)}_i|)
$$
なので、SPP粒子の粒子数密度$${ n^{(*)}_{i{\rm SPP}} }$$は、和の項が１項のみとなるので
$$
n^{(*)}_{i{\rm SPP}} = \omega(|\boldsymbol{r}^{(*)}_{i{\rm SPP}} - \boldsymbol{r}^{(*)}_i|) = \frac{r_e}{|\boldsymbol{r}^{(*)}_{i{\rm SPP}} - \boldsymbol{r}^{(*)}_i|} -1
$$
で表されるね。$${\boldsymbol{r}^{(*)}_{i{\rm SPP}}}$$は$${i}$$番目のSPP粒子の位置ベクトルだよ。上記の２つの式をつなげると、表面粒子$${i}$$とSPP粒子の距離は次のようになるね。
$$
|\boldsymbol{r}^{(*)}_{i{\rm SPP}} - \boldsymbol{r}^{(*)}_i|= \frac{r_e}{ n^{(0)} -n_i^{(*)} + 1 }
$$
次にSPP粒子の位置$${\boldsymbol{r}^{(*)}_{i{\rm SPP}}}$$を決定するために、$${i}$$番目の粒子の近傍粒子の重心を次の式で計算するよ。
$$
\boldsymbol{r}^{(*)}_{ig} = \frac{1}{n^{(0)}} \sum\limits_{j=1, j\ne i}^N \boldsymbol{r}^{(*)}_{ij}\, \omega (|\boldsymbol{r}^{(*)}_{ij}|)
$$
先に挙げた図のように$${\boldsymbol{r}^{(*)}_{i{\rm SPP}}}$$を重心から$${i}$$番目の粒子位置へ引いた直線上に配置するように考えると、
$$
\boldsymbol{r}^{(*)}_{i{\rm SPP}} = \boldsymbol{r}^{(*)}_{i} - |\boldsymbol{r}^{(*)}_{i{\rm SPP}} - \boldsymbol{r}^{(*)}_i|\, \frac{\boldsymbol{r}^{(*)}_{ig}}{|\boldsymbol{r}^{(*)}_{ig}|} = \boldsymbol{r}^{(*)}_{i} - \frac{r_e}{ n^{(0)} -n_i^{(*)} + 1 }\, \hat{\boldsymbol{t}}_{ig}
$$
と決定することができるね。なお、SPP粒子は流体粒子にせん断応力を与えないので速度$${\boldsymbol{v}_{{i{\rm SPP}}}^{(*)}}$$は流体粒子と同じで相対速度が０となるよ。
$$
\boldsymbol{v}_{{i{\rm SPP}}}^{(*)} = \boldsymbol{v}_{i}^{(*)} \ \to \ \boldsymbol{v}_{i{i{\rm SPP}}}^{(*)}=0
$$
実装まとめ
この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
$$
の右辺をラプラシアンモデルで置き換えた
$$
\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} |)\\
=\frac{2{\rm dim}}{\lambda^{(0)}n^{(0)}}\left[ \sum\limits_{j=1,\, j\ne i}^N P^{(k+1)}_j\, \omega(|\boldsymbol{r}^{(*)}_{ij}|) -P^{(k+1)}_i\sum\limits_{j=1,\, j\ne i}^N \omega(|\boldsymbol{r}^{(*)}_{ij}|)\right]
$$
を次のように書き換えるよ。
$$
\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{\rm SPP}}\, \omega_{i{\rm SPP}} -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]\\
=\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 n^{(0)}\right] + \frac{2{\rm dim}}{\lambda^{(0)}n^{(0)}}P^{(k+1)}_{i{\rm SPP}}\, \omega_{i{\rm SPP}}
$$
最後の項はポアソン方程式の左辺に移項すべきなのだけれども、ディリクレ条件として
$$
P^{(k+1)}_{i{\rm SPP}} = 0
$$
を与えるならば、消えてなくなる項だね。なお、上記の導出の過程で次の関係式を用いているよ。
$$
\sum\limits_{j=1,\, j\ne i, i{\rm SPP}}^N \omega(|\boldsymbol{r}^{(*)}_{ij}|) + \omega_{i{\rm SPP}} = n^{(0)}
$$
以上の結果、$${i}$$番目の粒子が表面粒子である場合、圧力を計算するための行列要素を次のように置き換えるよ。
$$
a_{ii} = \frac{1}{\rho_0}\,\frac{2{\rm dim}}{\lambda^{(0)}n^{(0)}}\sum\limits_{k=1,\, k\ne i}^N \omega(|\boldsymbol{r}^{(*)}_{ik}|)\ \to\ a_{ii} = \frac{1}{\rho_0}\,\frac{2{\rm dim}}{\lambda^{(0)}}
$$
２．時刻$${k+1}$$の速度を計算する際の圧力に関するグラディエントモデルの差し替え。
$$
\boldsymbol{v}^{(k+1)}_i = \boldsymbol{v}^{(*)}_i -\frac{1}{\rho_0} \left\langle \nabla P\, \right\rangle^{(k+1)}_i \Delta t
$$
圧力のグラディエントを次のように置き換えるよ。
$$
\langle \nabla P\, \rangle_i^{(k+1)} =\frac{{\rm dim}}{n^{(0)}} \sum\limits_{j=1,\, j\ne i, i{\rm SPP}}^N \frac{ P^{(k+1)}_j-P^{(k+1)}_i}{|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|^2}\,(\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i)\, \omega(|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|)
$$
→
$$
\langle \nabla P\, \rangle_i^{(k+1)} =\frac{{\rm dim}}{n^{(0)}}\left[ \sum\limits_{j=1,\, j\ne i, i{\rm SPP}}^N \frac{ P^{(k+1)}_j-P^{(k+1)}_i}{|\boldsymbol{r}^{(*)}_{ij}|^2}\,\boldsymbol{r}^{(*)}_{ij}\, \omega(|\boldsymbol{r}^{(*)}_{ij}|) + \frac{ P^{(k+1)}_{i{\rm SPP}}-P^{(k+1)}_i}{|\boldsymbol{r}^{(*)}_{{i{\rm SPP}}j}|^2}\,\boldsymbol{r}^{(*)}_{{ii{\rm SPP}}}\, \omega(|\boldsymbol{r}^{(*)}_{{ii{\rm SPP}}}|) \right]\\
=\frac{{\rm dim}}{n^{(0)}}\left[ \sum\limits_{j=1,\, j\ne i, i{\rm SPP}}^N \frac{ P^{(k+1)}_j-P^{(k+1)}_i}{|\boldsymbol{r}^{(*)}_{ij}|^2}\,\boldsymbol{r}^{(*)}_{ij}\, \omega(|\boldsymbol{r}^{(*)}_{ij}|) \right] + \frac{{\rm dim}}{n^{(0)}}\left[\frac{ P^{(k+1)}_{i{\rm SPP}}-P^{(k+1)}_i}{|\boldsymbol{r}^{(*)}_{{ii{\rm SPP}}}|^2}\,\boldsymbol{r}^{(*)}_{{ii{\rm SPP}}}\, \omega(|\boldsymbol{r}^{(*)}_{{ii{\rm SPP}}}|) \right]
$$
以上だよ。今回はアルゴリズムの導出のみだよ。次回、Pythonによる実装と検証を行っていくね。これからも応援よろしくお願いしまーす。