はじめに
一般的な流体シミュレーションでは空間をグリッドに分割して流体を表現します．これを「格子法」と呼びます．
一方，流体を粒子で表現しようとするのが「粒子法」です．粒子法にはいろんなメリットがあります（こちらのサイトで詳しく説明されています）．何よりも，個人的には粒子法の方が物理的に自然なモデリングかなぁと思っています．
この記事は「MPS法解説シリーズ」の一部です．このシリーズでは粒子法の一種である「MPS法」を解説しています．初めましての方は，ぜひシリーズの全体像からご覧ください．
目次
・やりたいこと
・微分演算モデルについて
・勾配モデル
・圧力勾配項の計算に用いる勾配モデル
・ベクトルの発散
・まとめ
・参考文献
やりたいこと
のちほど述べるように，流体計算には特定の微分演算モデルが必要です．この記事ではMPS法の勾配モデルを導出します．
スカラー $\phi$ の勾配モデル
\langle\boldsymbol{\nabla}\phi\rangle_i=\frac{d}{n^0}\sum_{j\ne i}\frac{\phi_{ij}}{|\boldsymbol{r}_{ij}|^2}\boldsymbol{r}_{ij}w_{ij}
この記事は「数学的な厳密さは置いといて，物理的な感覚を知りたいんだ！」という工学者気質な方向けです．特に，空間の次元数 $d$ をかける意味や粒子数密度の基準値 $n^0$ で割る意味を直観的に理解して頂けるように頑張ります．
数学的には，これらのモデルは「Taylar展開＋粒子配置の均等性」から厳密に導出できます．数学者気質な方はそちらの記事をご覧ください．（まだ作成中です...ごめんなさい）
なお，ナビエ・ストークス方程式の圧力勾配項の計算には次の勾配モデルを用います．詳しくはのちほど述べます．
\langle\boldsymbol{\nabla}P\rangle_i=\frac{d}{n^0}\sum_{j\ne i}\frac{P_j-\hat{P_i}}{|\boldsymbol{r}_{ij}|^2}\boldsymbol{r}_{ij}w_{ij}
微分演算モデルについて
まず，流体シミュレーションに必要な微分演算モデルについて説明しておきます．支配方程式を確認しておきましょう．
ナビエ・ストークス方程式（粒子法では実質微分が用いられることにご注意ください）
\frac{D\boldsymbol{u}}{Dt}
=-\frac{\boldsymbol{\nabla}P}{\rho}+\frac{\mu}{\rho}\boldsymbol{\nabla}^2\boldsymbol{u}+\boldsymbol{f}
連続の式
\boldsymbol{\nabla}\cdot\boldsymbol{u}=0
この2式にはスカラーとベクトルに関する1階偏微分と2階偏微分が登場するので，これらを離散化したモデルが必要です．
ベクトルの場合は各成分（スカラー）ごとに偏微分すればよいので，スカラーの偏微分モデルがあれば大丈夫です．例えば，ベクトル $\boldsymbol{u}$ の発散は次式のようにスカラーの偏微分のみで表せます．
\boldsymbol{\nabla \cdot u} = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z}
したがって，必要なのはスカラーの1階偏微分（勾配）と2階偏微分（ラプラシアン）のモデルということになります．この2つのうち，この記事では勾配モデルを導出します．残るラプラシアンについてはこちらで解説します．
勾配モデル
粒子 $i$ と粒子 $j$ の間のスカラー $\phi$ の局所的な勾配は，
(\boldsymbol{\nabla}\phi)_{ij}=\frac{\phi_{ij}}{|\boldsymbol{r}_{ij}|}\frac{\boldsymbol{r}_{ij}}{|\boldsymbol{r}_{ij}|}
と書けます．この式の重み平均を取ります．
\langle\boldsymbol{\nabla}\phi\rangle_i=\frac{1}{n_0}\sum_{j\ne i}\frac{\phi_{ij}}{|\boldsymbol{r}_{ij}|^2}\boldsymbol{r}_{ij}
最後に，次元数 $d$ 倍すればMPS法の勾配モデルが導出されます．
\langle\boldsymbol{\nabla}\phi\rangle_i=\frac{d}{n^0}\sum_{j\ne i}\frac{\phi_{ij}}{|\boldsymbol{r}_{ij}|^2}\boldsymbol{r}_{ij}w_{ij}
「重み平均を取る」とは？
「式の重み平均を取る」とは，「その式に重み関数をかける⇒すべての粒子に関する和を取る⇒粒子数密度の基準値 $n^0$ で割る」ことをいいます．ほかにも，「kernel積分する」といった言い方があります．
重み関数をかけて和を取ることで，近傍粒子の影響を大きく，遠方粒子の影響を小さく反映させることができます．
粒子数密度の基準値 $n^0$ で割るのは，重み関数をかけたことで元の値が変化してしまうのを防ぐ（規格化する）ためです．本来は局所的な粒子数密度 $n_i$ を用いるべきですが，非圧縮性の観点から粒子数密度は一定とみなせるので，簡単のため基準値 $n^0$ を用いています．
次元数倍するわけ
次元数 $d$ 倍する前の勾配モデルで足し合わせているのは，相対位置ベクトル $\boldsymbol{r_{ij}}$ の方向に1次元化された方向成分です．しかし，勾配ベクトルの計算には空間の次元数分の方向を考慮する必要があります．例えば，坂道の勾配を計算する際に坂を横切る方向のみを考慮しても正しい勾配は得られません．
したがって，例えば2次元空間における勾配モデルは
\langle
\boldsymbol{\nabla}\phi
\rangle_i
=\frac{1}{n_0}\sum_{j\ne i}
\Biggl\{
\frac{\phi_{ij}}{|\boldsymbol{r}_{ij}|^2}
\boldsymbol{r}_{ij}w_{ij}
+\frac{\phi_{ij_{⊥}}}{|\boldsymbol{r}_{ij_{⊥}}|^2}
\boldsymbol{r}_{ij_{⊥}}w_{ij_{⊥}}
\Biggl\}
とするのが妥当です．ここで，粒子 $j_{⊥}$ は粒子 $i$ から見て粒子 $j$ と等距離でかつ $ij$ 間と $ij_{⊥}$ 間の相対位置ベクトルが直交するような粒子です．
|\boldsymbol{r}_{ij_{⊥}}|=|\boldsymbol{r}_{ij}|\\
\boldsymbol{r}_{ij_{⊥}}\cdot\boldsymbol{r}_{ij}=0
3次元空間の場合，もう1つの垂直方向成分を足し合わせればOKです．
粒子 $i$ と粒子 $j$ のすべての組に対して，上の条件を満たす粒子 $j_{⊥}$ を見つけるのは困難です．というか，不可能です．そこで，計算の簡略化のために粒子 $i$ 近傍の粒子配置が均等であると仮定します．すると，ある粒子 $j$ は空間の次元数に応じて2回または3回，他の粒子 $j'$ に対する粒子 $j'_{⊥}$ となります．これが，次元数倍する理由です．
圧力勾配項の計算に用いる勾配モデル
上で定義したMPS法の勾配モデルを用いると，圧力 $P$ の勾配は
\langle\boldsymbol{\nabla}P\rangle_i=\frac{d}{n^0}\sum_{j\ne i}\frac{P_{ij}}{|\boldsymbol{r}_{ij}|^2}\boldsymbol{r}_{ij}w_{ij}
となるはずですが，ナビエストークス方程式中で圧力を求める際に用いる勾配計算には
\langle\boldsymbol{\nabla}P\rangle_i=\frac{d}{n^0}\sum_{j\ne i}\frac{P_j-\hat{P_i}}{|\boldsymbol{r}_{ij}|^2}\boldsymbol{r}_{ij}w_{ij}
を用います．すなわち， $P_i$ の代わりに $\hat{P_i}$ を用いています． $\hat{P_i}$ は粒子 $i$ の近傍で最も圧力が小さい粒子の圧力です．なお，ここでいう近傍は重み関数に用いる影響半径 $r_e$ 内のことを指します．
このように変形すると圧力勾配によって生じる力は常に斥力になります．これにより粒子が近づきすぎることを防ぎ，計算の安定化を図っています．
計算のためとはいえ，勾配モデルをこのように変形してしまうと物理的直観から離れてしまいますし，計算誤差が気になりますね．この点については別の記事で詳しくお話するつもりです．
ベクトルの発散
スカラーの勾配モデルが求まったので，そこからベクトルの発散を導出しておきます（いつか使うかもしれませんし，使わないかもしれません）．すでに述べたように，ベクトルの発散は次のように表されます．
\boldsymbol{\nabla \cdot u} = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z}
ここで，スカラーである $u$ に上で求めた勾配モデルを適用してみましょう．
\begin{align}
\langle\boldsymbol{\nabla}u\rangle_i=\frac{d}{n^0}\sum_{j\ne i}\frac{u_{ij}}{|\boldsymbol{r}_{ij}|^2}\boldsymbol{r}_{ij}w_{ij}
\end{align}
分かりやすいように成分ごとに見てみます．
math \begin{align} \begin{pmatrix} \langle \frac{\partial u}{\partial x} \rangle_i \\ \langle \frac{\partial u}{\partial y} \rangle_i \\ \langle \frac{\partial u}{\partial z} \rangle_i \end{pmatrix} &=\frac{d}{n^0}\sum_{j\ne i}\frac{u_{ij}}{|\boldsymbol{r}_{ij}|^2} \begin{pmatrix} x_{ij} \\ y_{ij} \\ z_{ij} \end{pmatrix} w_{ij} \\ \\ \therefore \bigg\langle\frac{\partial u}{\partial x}\bigg\rangle_i &=\frac{d}{n^0}\sum_{j\ne i}\frac{u_{ij}x_{ij}}{|\boldsymbol{r}_{ij}|^2}w_{ij} \end{align}
$v$， $w$ についても同様に考えて，発散の式に代入します．
`math
\begin{align}
\langle\boldsymbol{\nabla \cdot u}\rangle_i
&=\frac{d}{n^0}\sum_{j\ne i}\frac{u_{ij}x_{ij}}{|\boldsymbol{r}{ij}|^2}w{ij}
- \frac{d}{n^0}\sum_{j\ne i}\frac{v_{ij}y_{ij}}{|\boldsymbol{r}{ij}|^2}w{ij}
- \frac{d}{n^0}\sum_{j\ne i}\frac{w_{ij}z_{ij}}{|\boldsymbol{r}{ij}|^2}w{ij}\
&=\frac{d}{n^0}\sum_{j\ne i}\frac{u_{ij}x_{ij}+v_{ij}y_{ij}+w_{ij}z_{ij}}{|\boldsymbol{r}{ij}|^2}w{ij} \
&=\frac{d}{n^0}\sum_{j\ne i}\frac{\boldsymbol{u}{ij}\cdot\boldsymbol{r}{ij}}{|\boldsymbol{r}{ij}|^2}w{ij}
\end{align}
`
以上でベクトルの発散モデルが求まりました．
まとめ
MPS法の勾配モデルと発散モデルの導出を行いました．
私は特に次元数の話で直観的に理解できずに苦しみましたが，坂の例でだいぶしっくりきました．
皆さんのお役に立てれば幸いです．質問やご意見お待ちしてます．