【粒子法 MPS法】アルゴリズムの改善：DS法 人工斥力導入による安定化
上記のnoteで解説したGC法に加えて、下記の書籍で詳しく解説してある標準MPSの高精度化に着手しているよ。目次は次のとおりだよ。
書籍：粒子法:連続体・混相流・粒状体のための計算科学
P.83 HS法 ポアソン方程式生成項の高精度化
P.85 HL法 ラプラシアンモデルの高精度化
P.87 ECS法 ポアソン方程式生成項の誤差補償
P.92 GC法 勾配モデルの高精度化
P109 DS法 人工斥力導入による安定化
P118 SPP法 表面粒子の境界条件適正化
P124 WPP法 壁面における境界条件適正化
書籍の順番に従って、人工斥力導入による安定化であるDS法の導入を行うよ。MPS法では粒子同士の重なりが生じると瞬間的な圧力の異常上昇を引き起こすのだけれども、これは計算不安定性の原因となるね。標準MPS法ではこのような粒子の重なりを排除するために、①仮時刻における位置を計算した時点で粒子同士に重なりが生じた場合に例外的衝突撃力を与える、②圧力ポアソン方程式を解いた後の速度更新
$$
\boldsymbol{v}^{(k+1)}_i = \boldsymbol{v}^{(*)}_i -\frac{1}{\rho_0} \left\langle \nabla P\, \right\rangle^{(k+1)}_i \Delta t
$$
のグラディエントモデルを$${i}$$番目の粒子に加わる圧力を生の$${ P^{(k+1)}_i }$$から$${i}$$番目の粒子の周囲の最も小さな圧力$${\hat{P}^{(k+1)}_i }$$に置き換えた
$$
\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}^{(*)}_j - \boldsymbol{r}^{(*)}_i|^2}\,(\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i)\, \omega(|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|)
$$
を用いて計算することで圧力勾配による粒子に加わる力をすべて斥力とする、という以上２つの粒子重なりの排除のアルゴリズムが導入されているね（これまで詳しく解説していなかったけれども…）。どちらともナビエ・ストークスとは無関係に導入された非物理的な効果だけれども、②はとくに圧力勾配が斥力しか存在できないことで粒子密度が低い領域への粒子が流れ込む負圧域の計算ができないことは問題だね。
この２つの施策を改善するのがDS法だよ。まず圧力勾配を負圧も扱えるように $${\hat{P}^{(k+1)}_i \to P^{(k+1)}_i }$$ とした
$$
\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-P^{(k+1)}_i}{|\boldsymbol{r}^{(*)}_{ij}|^2}\,\boldsymbol{r}^{(*)}_{ij}\, \omega(|\boldsymbol{r}^{(*)}_{ij}|)
$$
を用いるよ。この圧力勾配から時刻$${k+1}$$の候補とした仮速度$${ \boldsymbol{v}^{(k+1)^*}_i }$$と仮位置$${ \boldsymbol{r}^{(k+1)^*}_i }$$を次のように計算するよ。
$$
\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{ F}_{i}^{\, ({\rm DS)}}}$$を加えることで重なりを取り除くよ。
上図は人工斥力$${\boldsymbol{ F}_{i}^{\, ({\rm DS)}}}$$を計算する際に必要なベクトル量を表しているよ。$${\boldsymbol r_i^{(*)}}$$と$${\boldsymbol r_j^{(*)}}$$はこれまでと同様、圧力を無視して計算した仮位置だよ。$${\boldsymbol{r}^{(k+1)^*}_i }$$と$${\boldsymbol{r}^{(k+1)^*}_j }$$は負圧を考慮した圧力勾配で計算した次の時刻$${k+1}$$の仮位置と仮速度だよ。この図のように$${i}$$番目と$${j}$$番目の粒子が重なった場合（青色の丸）、引き離してちょうど接するところまでに相当する人工斥力$${\boldsymbol{ F}_{i}^{\, ({\rm DS)}}}$$をそれぞれの粒子に加えてるね。DS法では$${\boldsymbol{ F}_{i}^{\, ({\rm DS)}}}$$の方向を仮時刻$$(*)$$時の２粒子間の方向ベクトル（$${\hat{\boldsymbol {t}}_{ij} =\boldsymbol r^{(*)}_{ij}/|\boldsymbol r^{(*)}_{ij}| }$$）とするよ。そして、$${i}$$番目の粒子に加わる人工斥力はすべての粒子に関しての和
$$
\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}|)
$$
で与えるよ。$${\rho_i}$$は$${i}$$番目の粒子の密度だけれども今は密度一定で考えているから$${\rho_i \simeq \rho^{0} }$$と考えるよ。そして、この２粒子間の人工斥力は次の式で定義されるね。
$$
\boldsymbol{ F}_{ij}^{\, ({\rm DS)}}=\left\{ \begin{matrix} 0 & (|\boldsymbol r^{(k+1)^*}_{ij}| \geq d_{ij}) \cr \Pi _{ij} \,\hat{\boldsymbol {t}}_{ij} & (|\boldsymbol r^{(k+1)^*}_{ij}| < d_{ij}) \end{matrix} \right. \ , \ \ \hat{\boldsymbol {t}}_{ij}= \frac{\boldsymbol r^{(*)}_{ij}}{|\boldsymbol r^{(*)}_{ij}|}
$$
$${ d_{ij}}$$は人工斥力を発生させる距離、$${\Pi _{ij}}$$は重なり状態から接する状態までに必要な人工斥力の大きさを表すよ。具体的な導出は割愛（書籍を参照）するけれども、次のように与えられるよ。
$$
d_{ij} = \alpha_{\rm DS}\, \frac{d_i+d_j}{2} =\left( 1-\frac{v_{\rm max} \Delta t}{d_i} \right) \frac{d_i+d_j}{2}
$$
$$
\Pi _{ij} = \frac{ \rho_j }{\Delta t^2(\rho_i + \rho_j)} \left(\sqrt{ d^2_{ij}-| \boldsymbol r^{(k+1)^*}_{ij} \cdot \hat{\boldsymbol{n}}_{ij} |^2 } -|\boldsymbol r^{(k+1)^*}_{ij} \cdot \hat{\boldsymbol{t}}_{ij} | \right)
$$
$${\hat{\boldsymbol{n}}_{ij}}$$は２次元の場合は単純に$${\hat{\boldsymbol{t}}_{ij}}$$と垂直な単位ベクトルだよ。また、$${d_i, d_j}$$は粒子半径、$${ v_{\rm max} }$$は計算対象の流速の最大値で因子$${ v_{\rm max} \Delta t / d_i}$$はクーランド数で、$${\alpha_{\rm DS}}$$が例えば0.9の場合、10%の重なりで斥力が発生させるという意味を持つよ。ここで計算した$${\boldsymbol{ F}_{i}^{ ({\rm DS)}} }$$を用いて時刻$${k+1}$$の速度ベクトルと位置ベクトルは次のように決定されるね。
$$
\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
$$
今回の導出過程では圧力勾配のグラディエントモデルを標準MPS法で表したけれども、前回導出したCG法におけるグラディエントモデルにも適用できるよ。先と同様に$${\hat{P}^{(k+1)}_i \to P^{(k+1)}_i }$$とした
$$
\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-P^{(k+1)}_i}{|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|^2} \,\boldsymbol{C}^{(*)}_{i} \cdot (\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i)\, \omega(|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|)
$$
$$
\boldsymbol{C}^{(*)}_{i}= \left[ \frac{1}{n^{(0)}}\sum\limits_{j=1,\, j\ne i}^N\frac {\omega(|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|)}{|\boldsymbol{r}^{(*)}_j - \boldsymbol{r}^{(*)}_i|^2}\left( \boldsymbol r^{(*)}_{ij}\otimes\boldsymbol r^{(*)}_{ij} \right) \right]^{-1}
$$
を用いるだけだね。
今回はアルゴリズムの導出のみだよ。次回、Pythonによる実装と検証を行っていくね。これからも応援よろしくお願いしまーす。