MPS法による計算アルゴリズム

粒子法とはナビエ・ストークス方程式のような偏微分方程式を数値的に解くための手法の一つです。粒子法は解析対象とする流体を細かな粒子の集合体とみなして、粒子同士が相互作用しながら時間発展していくという描像だね。粒子法の利点は有限要素法のように格子をあらかじめに用意する必要がないので、自由表面流れや場合によっては水の飛沫のような極端な境界が存在する場合でも、特段の処置が必要ない点が最大の利点として知られているね。また、ナビエ・ストークス方程式は空間に固定された各地点における物理量の関係を記述するオイラー記述よりも、粒子の運動とみなすラグランジュ記述のほうが簡単なことも大きな利点だね。本稿では、非圧縮流体の解析法としてスタートした粒子法の一つであるMPS法の基本を解説するよ。

参考書籍

粒子法入門 流体シミュレーションの基礎から並列計算と可視化まで C/C++ソースコード付

流体→粒子

MPS法では流体を$${N}$$個の粒子で表現するよ。$${i}$$番目の粒子の位置と速度を $${\boldsymbol{ r}_i, \boldsymbol{ v}_i }$$ と表し、粒子同士の相互作用の大きさを重み関数

$$
\omega (|\boldsymbol r_j-\boldsymbol r_i|):= \left\{
\begin{matrix}
\left(\frac{r_e}{r}\right) -1 &  ( r < r_e) \\  0 \ \ \ \ & (r_e \leq r)
\end{matrix}
\right.
$$

で表すよ。$${r_e}$$ は影響半径と呼ばれ、粒子の大きさを表すよ。この重み関数を用いて、粒子数密度と流体密度が定義されるよ。

粒子数密度

$$
n_i :=\sum\limits_{j\ne i} \omega(|\boldsymbol{r}_j - \boldsymbol{r}_i|)
$$

流体密度

$$
\frac{\rho_i^{(k)}-\rho^{(0)}}{\rho^{(0)}} \simeq \frac{n_i^{(k)}-n^{(0)}}{ n^{(0)}}
$$

$${ n^{(0) }$$ は粒子を影響半径の敷き詰めたときの粒子数密度だよ。$${ \rho^{(0)} }$$ は基準となる流体密度で、水の場合は $${\rho^{(0)}=1000{\rm[kg/m^3]}}$$ と与えられるね。また、重み関数を用いて、任意の物理量に対するグラディエントとラプラシアンを次のように定義されるよ。

グラディエントモデル

$$
\langle \nabla \phi \rangle_i :=\frac{{\rm dim}}{n^{(0)}} \sum\limits_{j=1}^N \frac{ \phi_j-\phi_i}{|\boldsymbol{r}_j - \boldsymbol{r}_i|^2}\,(\boldsymbol{r}_j - \boldsymbol{r}_i)\, \omega(|\boldsymbol{r}_j - \boldsymbol{r}_i|)
$$

ラプラシアンモデル

$$
\langle \nabla^2 \phi \rangle_i :=\frac{2{\rm dim}}{\lambda^{(0)}n^{(0)}} \sum\limits_{j=1}^N (\phi_j-\phi_i)\, \omega(|\boldsymbol{r}_j - \boldsymbol{r}_i|)
$$

$${\lambda^{(0)}}$$ は表面ではない粒子の初期配置に対して

$$
\lambda^{(0)}:=\frac{\sum\limits_{j\ne i}|\boldsymbol{r}_j^{(0)} - \boldsymbol{r}_i^{(0)}|^2 \omega(|\boldsymbol{r}_j^{(0)} - \boldsymbol{r}i^{(0)}|) }{\sum\limits{j\ne i} \omega(|\boldsymbol{r}_j^{(0)} - \boldsymbol{r}_i^{(0)}|)}
$$

で与えられる量だよ。以上のグラディエントとラプラシアンを対象とするナビエ・ストークス方程式に適用することで、流体の運動を粒子の運動に置き換えることを考えるのがMPS法だよ。

ナビエ・ストークス方程式と連続の式

非圧縮流体のナビエ・ストークス方程式は次のように与えられるね。

$$
\frac{D \boldsymbol{v}}{Dt} = -\frac{1}{\rho_0} \nabla P + \nu\nabla^2\boldsymbol{v} + \boldsymbol{g}
$$

$${\boldsymbol{v}}$$ は速度ベクトル、$${P}$$ は圧力（$${\rm [Pa/m^2]}$$）、$${\rho_0}$$ は密度（$${\rm [kg/m^3]}$$）、$${\nu}$$ は動粘性係数（粘性係数$${ \mu }$$を$${\rho}$$で割った値）、$${\boldsymbol{g}}$$は重力加速度（$${\rm [m/s^2]}$$）だよ。$${D}$$はラグランジュ記述の微分演算子だね。もう一つ重要な方程式は質量保存則を表す連続の式

$$
\frac{D \rho}{Dt} +\rho\nabla\cdot\boldsymbol{v} =0
$$

だね。非圧縮流体の場合、第１項目がゼロなので、任意の時刻で速度は

$$
\nabla\cdot\boldsymbol{v} =0
$$

を満たすことになるね。この式とナビエ・ストークス方程式で流体の運動を計算することができるよ。粒子法ではこのナビエ・ストークス方程式を質点のような粒子の加速度を与える式と考えるので、ラジュランジュ記述の微分演算子は通常の微分演算子と同様に扱うことができて、$${i}$$ 番目の粒子の加速度は

$$
\frac{d \boldsymbol{v}_i}{dt} = -\frac{1}{\rho_0}\left\langle \nabla P \right\rangle_i + \nu\left\langle \nabla^2\boldsymbol{v} \right\rangle _i + \boldsymbol{g}
$$

で与えられるね。式中の $${{\langle \sim \rangle}_i}$$ は $${i}$$ 番目の粒子の周りの環境によって与えられる量で、先に解説したグラディエントモデルとラプラシアンモデルで定義されるね。MPSでは速度の時間発展を１次近似

$$
\frac{d \boldsymbol{v}_i}{dt} \simeq \frac{ \boldsymbol{v}^{(k+1)}_i - \boldsymbol{v}^{(k)}_i }{\Delta t}
$$

で表して（上付き添字$$(k)$$は時間ステップ数だよ）、右辺の速度の時刻は$$(k)$$として置きつつ、圧力は$$(k+1)$$という１ステップ未来の時刻の値を用いるよ。式で表すと次のとおりだね。

$$
\boldsymbol{v}^{(k+1)}_i =\boldsymbol{v}^{(k)}_i+\left[ -\frac{1}{\rho_0} \left\langle \nabla P \right\rangle^{(k+1)}_i + \nu \left\langle \nabla^2\boldsymbol{v} \right\rangle^{(k)}_i + \boldsymbol{g}\right] \Delta t
$$

このように、未来の時刻の速度を計算するのに未来の時刻の量（$${\left\langle \nabla P \right\rangle^{(k+1)}_i }$$）を用いるアルゴリズムは陰的数値解法と呼ばれるね。　MPS法ではこの未知の圧力を計算するために、次の手順をとるよ。

圧力に関するポアソン方程式の導出と計算方法

１．仮速度の導入と時刻$${k+1}$$の速度の表式

先の式で圧力が存在しないと仮定したときの速度（仮速度）を

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

と定義しておいて、この式に速度に対するラプラシアンモデルを適用すると

$$
\boldsymbol{v}^{(*)}_i = \boldsymbol{v}^{(k)}_i + \left[ \frac{2 \nu{\rm dim}}{ \lambda^{(0)} n^{(0)} } \sum\limits_{j=1,, j\ne i}^N (\boldsymbol{v}_j^{(k)} - \boldsymbol{v}_i^{(k)})\, \omega(|\boldsymbol{r}_j^{(k)} - \boldsymbol{r}_i^{(k)}|) + \boldsymbol{g} \right]\Delta t
$$

となるね。右辺の値は時刻 $${k}$$ の粒子の位置と速度ですべて既知なので、$${\boldsymbol{v}^{(*)}_i }$$は計算することができるね。ちなみに、時刻$${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
$$

２．仮位置と時刻$${k+1}$$の位置を定義

この速度を用いて圧力が存在しないと仮定したときの位置（仮位置）を前方差分

$$
\boldsymbol{r}^{(*)}_i = \boldsymbol{r}^{(k)}_i +\boldsymbol{v}^{(*)}_i\Delta t
$$

で定義するよ。$${\boldsymbol{v}^{(*)}_i }$$ が与えられれば、$${\boldsymbol{r}^{(*)}_i }$$ も直ちに得られるね。一方、時刻$${k+1}$$の位置は後方差分

$$
\boldsymbol{r}^{(k+1)}_i = \boldsymbol{r}^{(k)}_i +\boldsymbol{v}^{(k+1)}_i\Delta t
$$

で定義するよ。

４．圧力が満たす条件

時刻$${k+1}$$の速度の式の両辺に発散をとると

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

となるね。左辺の量は質量保存則からゼロになることは先に示したね。つまり、

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

となるね。一方、仮時刻における速度と密度は元の連続の式を満たすので

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

となるわけなので、

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

と粒子数密度で近似することができるね。つまり、圧力が満たす方程式は次のとおりだよ。

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

この方程式はポアソン方程式の形になっているね。

５．圧力の計算方法

圧力の条件式の右辺の量は仮位置からすべて計算できるね。一方、左辺はラプラシアンモデルに置き換えると

$$
\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}^{(*)}_j - \boldsymbol{r}^{(*)}_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)^2}\,\frac{n_i^{(*)}-n^{(0)}}{n^{(0)}}
$$

という連立方程式が得られるね。これを行列で表すと次のようになるよ。

$$
\left(\begin{matrix} a_{11} & a_{12} & a_{13} & \cdots \\ a_{21} & a_{22} & a_{23} & \cdots \\ a_{31} & a_{32} & a_{33} & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{matrix} \right)\left(\begin{matrix} P_1\\ P_2\\ P_3\\ \vdots \end{matrix} \right)=\left(\begin{matrix} b_1 \\ b_2\\ b_3 \\ \vdots \end{matrix} \right)
$$

$$
a_{ij} = \left\{\begin{matrix}- \frac{1}{\rho_0}\,\frac{2{\rm dim}}{\lambda^{(0)}n^{(0)}}\ \omega(|\boldsymbol{r}^{(*)}_j- \boldsymbol{r}^{(*)}_i|) & (i,\ne, j) \\  \frac{1}{\rho_0}\,\frac{2{\rm dim}}{\lambda^{(0)}n^{(0)}}\sum\limits_{j'=1,, j'\ne i}^N \omega(|\boldsymbol{r}^{(*)}_{j'}- \boldsymbol{r}^{(*)}_i|)& (i= j)\end{matrix}\right.
$$

$$b_i = \frac{1}{(\Delta t)^2}\,\frac{n_i^{(*)}-n^{(0)}}{n^{(0)}}$$

行列で簡単に表すと

$$
\boldsymbol{A}\boldsymbol{P}=\boldsymbol{b}
$$

となるわけだけれども、$${\boldsymbol{A}}$$ の逆行列を計算することができれば、圧力が得られるね。

$$\boldsymbol{P} = \boldsymbol{A}^{-1} \boldsymbol{b}$$

位置と速度の更新

先に計算した圧力を用いると時刻$${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
$$

$$
\boldsymbol{r}^{(k+1)}_i = \boldsymbol{r}^{(k)}_i +\boldsymbol{v}^{(k+1)}_i\Delta t
$$

以上だよ。

MPS法の Pythonプログラムソース

僕が作成したPythonプログラムソースを販売するよ。今後、機能を追加するたびに更新していくね。販売するプログラムソースの動作チェック動画は次のとおりだよ。

公開日：2024年9月6日（ver.1.0）