はじめに

　一般的な流体シミュレーションでは空間をグリッドに分割して流体を表現します．これを「格子法」と呼びます．
　一方，流体を粒子で表現しようとするのが「粒子法」です．粒子法にはいろんなメリットがあります（こちらのサイトで詳しく説明されています）．何よりも，個人的には粒子法の方が物理的に自然なモデリングかなぁと思っています．
　この記事は「MPS法解説シリーズ」の一部です．このシリーズでは粒子法の一種である「MPS法」を解説しています．初めましての方は，ぜひシリーズの全体像からご覧ください．

目次

・やりたいこと
・計算アルゴリズム
・圧力のポアソン方程式を導出する
・まとめ

やりたいこと

　本記事ではMPS法の計算アルゴリズムを紹介し，その中で用いる「圧力のポアソン方程式」を導出します．
　論文などを読んでいると，圧力のポアソン方程式の導出には様々な方法があることが分かります．途中式や考え方の省略が多くて，理解するのに苦労しました...
　この記事では（私にとって）一番納得のいく導出を目指します．

#計算アルゴリズム
　MPS法の計算アルゴリズムには様々なものがありますが，一番基礎的なアルゴリズムは半陰解法を用いています．すなわち，圧力以外を陽的に計算したのち，非圧縮条件から圧力を陰的に計算します．この時，圧力の計算に用いるのが「圧力のポアソン方程式」です．
　第 $k$ ステップから第 $k+1$ ステップへの計算アルゴリズムは以下のようになっています．

・ナビエ・ストークス方程式の圧力勾配項以外の項を計算し，仮の速度ベクトル $\boldsymbol{u}_i^*$ を求める．

\boldsymbol{u}_i^*
=\boldsymbol{u}_i^k
+ \nu \langle \boldsymbol{\nabla}^2\boldsymbol{u} \rangle_i^k \Delta t
+ \rho \boldsymbol{f} \Delta t

・仮の速度ベクトル $\boldsymbol{u}_i^*$ を用いて仮の位置 $\boldsymbol{r}_i^*$ に粒子を移動させる．

\boldsymbol{r}_i^*
=\boldsymbol{r}_i^k
+\boldsymbol{u}_i^*\Delta t

・仮の粒子数密度 $n_i^*$ を計算する．

n_i^*
=\sum_{j\ne i}w_{ij}

・圧力のポアソン方程式（後述）を用いて，次の時刻における圧力 $P_i^{k+1}$ を計算する．

\langle \boldsymbol{\nabla}^2 P \rangle_i^{k+1}
=-\frac{\rho^0}{(\Delta t)^2} \frac{n_i^*-n^0}{n^0}

・圧力勾配項による速度増分を足し合わせて，真の速度ベクトル $\boldsymbol{u}_i^{k+1}$ を計算する．

\boldsymbol{u}_i^{k+1}
=\boldsymbol{u}_i^*
-\frac{1}{\rho^0} \langle \boldsymbol{\nabla} P \rangle _i^{k+1} \Delta t

・真の速度ベクトル $\boldsymbol{u}_i^{k+1}$ を用いて，真の位置ベクトル $\boldsymbol{r}_i^{k+1}$ に粒子を移動させる．

\boldsymbol{r}_i^{k+1}
=\boldsymbol{r}_i^k
+\boldsymbol{u}_i^{k+1}\Delta t

圧力のポアソン方程式を導出する

　圧力のポアソン方程式とは流体の非圧縮条件を考慮して得られる粒子 $i$ が次の時刻に持つ圧力 $P_i^{k+1}$ に関する方程式です，なお，ポアソン方程式とは一般に2階の楕円型偏微分方程式のことを言います．
　粒子 $i$ が仮の位置 $\boldsymbol{r}_i^*$ に移動した時点では，非圧縮条件は満たされません．そこで，次の時刻の圧力 $P_i^{k+1}$ によって速度の修正が生じ，非圧縮条件を満たす真の速度 $\boldsymbol{u}_i^{k+1}$ が得られると考えます．

math \boldsymbol{u}_i^{k+1} =\boldsymbol{u}_i^* -\frac{1}{\rho^0} \langle \boldsymbol{\nabla} P \rangle _i^{k+1} \Delta t
真の速度 $\boldsymbol{u}_i^{k+1}$ は非圧縮流体の連続の式 $\boldsymbol{\nabla}\cdot\boldsymbol{u}=0$ を満たすことを考慮して，上式の発散を取ってみましょう．

math \begin{align} &0=\boldsymbol{\nabla}\cdot\boldsymbol{u}_i^* -\frac{1}{\rho^0} \langle \boldsymbol{\nabla}^2 P \rangle _i^{k+1} \Delta t \\ &\therefore \langle \boldsymbol{\nabla}^2 P \rangle _i^{k+1} =\frac{\rho^0}{\Delta t}\boldsymbol{\nabla}\cdot\boldsymbol{u}_i^* \end{align}
　この式でも次の時刻の圧力 $P_i^{k+1}$ を現時点で分かっている値を用いて表すことができていますので，この式を圧力のポアソン方程式と呼ぶこともできます．しかし，この式のように速度の発散を用いると計算誤差に起因する密度の誤差が検出されず，体積が保存されにくい場合があります（この点については越塚先生がこちらで詳しく説明されていますので，ぜひご覧ください）．したがって，上式の右辺を変形し，（粒子数）密度の条件を用いて表せないか考えていきます．
　流体密度 $\rho_i^k(=\rho^0)$ から仮の流体密度 $\rho_i^*$ への密度変化は，仮の速度 $\boldsymbol{u}_i^*$ の影響で生じています．このことを，一般的な連続の式

math \frac{D\rho}{Dt} +\rho \boldsymbol{\nabla}\cdot\boldsymbol{u} =0
を用いて表現すると，

math \begin{align} \frac{\rho_i^*-\rho^0}{\Delta t} +\rho^0 \boldsymbol{\nabla} \cdot \boldsymbol{u}_i^* =0 \\ \therefore \boldsymbol{\nabla} \cdot \boldsymbol{u}_i^* =-\frac{1}{\Delta t} \frac{\rho_i^*-\rho^0}{\rho^0} \end{align}
となります．流体密度と粒子数密度の関係より，流体密度の代わりに粒子数密度を用いると，

math \boldsymbol{\nabla} \cdot \boldsymbol{u}_i^* =-\frac{1}{\Delta t} \frac{n_i^*-n^0}{n^0}
となり，速度の発散を粒子数密度で表現することができました．
　これを先ほどのポアソン方程式に代入すれば，

math \langle \boldsymbol{\nabla}^2 P \rangle_i^{k+1} =-\frac{\rho^0}{(\Delta t)^2} \frac{n_i^*-n^0}{n^0}
となり，MPS法で用いる圧力のポアソン方程式が導出されました．

まとめ

　この記事では，MPS法の計算アルゴリズムの紹介と圧力のポアソン方程式の導出を行いました．
　MPS法の計算アルゴリズムには今回紹介した半陰解法だけでなく陽解法も存在します．また，計算精度を上げるために圧力のポアソン方程式の形式も様々なものが提案されています．
　今回はとりあえずシンプルな形を導出しましたが，今後はより深くまとめられたらなぁと考えています．

参考文献

・粒子法（越塚先生）
・粒子法入門
・粒子法（後藤先生）