この記事ではMPS法のラプラシアンモデルをまず導出し、その次にSPH法のラプラシアンモデルと比較します。

MPS法のラプラシアンモデル

\nabla^2_{\rm M} F_i = A \sum_j \left( F_j - F_i \right) w_{ij} \tag{1}

ラプラシアンモデルを式(1)の形で記述できるとします (後ほど出てくる真のラプラシアン $\nabla^2$ と区別するために、$\nabla^2_{\rm M}$ としています)。ここで、$F$ は任意の物理量、 $w_{ij}$ は重み関数です。$A$ はこれから求めます。
$F_j$ を $\boldsymbol{r}_i$ 周りでテイラー展開すると

F_j = F_i + \nabla F_i \cdot \boldsymbol{r}_{ji} + \frac{1}{2} \left( \boldsymbol{r}_{ji} \cdot \nabla \right)^2 F_i + \mathcal{O}(r_{ji}^3)

となります。ここで、 $\boldsymbol{r}_{ji} = \boldsymbol{r}_j - \boldsymbol{r}_i$ としています。これを式(1) に代入します。

\nabla^2_{\rm M} F_i = A \sum_j \nabla F_i \cdot \boldsymbol{r}_{ji} w_{ij} + \frac{A}{2} \sum_j \left[ \left( \boldsymbol{r}_{ji} \cdot \nabla \right)^2 F_i \right] w_{ij} + \mathcal{O}(A r_{ji}^3) \tag{2}

まず、式(2)の右辺第1項は

 A \sum_j \nabla F_i \cdot \boldsymbol{r}_{ji} w_{ij} = A \left[ \sum_j \boldsymbol{r}_{ji} w_{ij} \right] \cdot \nabla F_i

で、格子状の粒子配置を仮定すると、対称性から

\sum_j \boldsymbol{r}_{ji} w_{ij} = 0

となるため無視することができます。

次に右辺第2項を展開します。3次元だと書くのが面倒なので2次元とします。

\frac{A}{2} \sum_j \left[ \left( \boldsymbol{r}_{ji} \cdot \nabla \right)^2 F_i \right] w_{ij}
= \frac{A}{2} \left[ \sum_j x_{ji}^2 w_{ij} \right] \partial_x^2 F_i + \frac{A}{2} \left[ \sum_j y_{ji}^2 w_{ij} \right] \partial_y^2 F_i + A \left[ \sum_j x_{ji} y_{ji} w_{ij} \right] \partial_x \partial_y F_i \tag{3}

この式(3)の右辺第3項については格子状の粒子配置の仮定から

\sum_j x_{ji} y_{ji} w_{ij} = 0

となるため無視できます。式(3)の右辺第1項と第2項について、またまた粒子配置の仮定を使うと

\sum_j x_{ji}^2 w_{ij} = \sum_j y_{ji}^2 w_{ij} = \frac{1}{d} \sum_j \boldsymbol{r}_{ji}^2 w_{ij}

となります。ここで、$d$ は次元です (2次元なので2)。これらを使って式(3)を整理すると、

\frac{A}{2} \sum_j \left[ \left( \boldsymbol{r}_{ji} \cdot \nabla \right)^2 F_i \right] w_{ij}
= \frac{A}{2d} \left[ \sum_j \boldsymbol{r}_{ji}^2 w_{ij} \right] \nabla^2 F_i

が得られます。3次元の場合でも同様です。

ここまでの議論から、式(2)は

\nabla^2_{\rm M} F_i = \frac{A}{2d} \left[ \sum_j \boldsymbol{r}_{ji}^2 w_{ij} \right] \nabla^2 F_i + \mathcal{O}(A r_{ji}^3)

と変形でき、係数 $A$ は

A= \frac{2d}{\sum_j \boldsymbol{r}_{ji}^2 w_{ij}} \tag{4}

とすればよいということがわかります。この $A$ は実際にMPSのラプラシアンモデルで使われてる係数と同じです (ただし、実際のシミュレーションでは、初期粒子配置から計算した定数値を $A$ として使います)。また、$A$ は $\mathcal{O}(r_{ij}^{-2})$ なので、ラプラシアンモデルの誤差は $\mathcal{O}(r_{ij})$ 、空間1次精度であることもわかります。

ただし、これは粒子配置の等方性が成立する場合の話で、粒子配置が乱雑になっている場合は、式(2)の右辺第1項が $\mathcal{O}(r_{ij}^{-1})$、第2項が $\mathcal{O}(1)$ の誤差を発生させることに注意してください。この誤差に関してはTamai et al. (2013)などに詳しく書かれています。

SPH法のラプラシアンモデル

SPH法ではラプラシアンはカーネル関数 $W_{ij}$ を使って

\nabla^2_{\rm S} F_i = 2 \sum_j \frac{m_j}{\rho_j} (F_j - F_i) \frac{\boldsymbol{r}_{ji}}{r_{ji}^2} \cdot \nabla W_{ij} \tag{5}

のように書かれます (Price 2012)。ここで重み関数を

w_{ij} = \frac{m_j}{\rho_j} \frac{\boldsymbol{r}_{ji}}{r_{ji}^2} \cdot \nabla W_{ij} \tag{6}

として、式(1)のMPS法のラプラシアンモデルを考えます。係数 $A$ は式(4)から

A= \frac{2d}{\sum_j \boldsymbol{r}_{ji}^2 w_{ij}} = \frac{2d}{\sum_j \frac{m_j}{\rho_j} \boldsymbol{r}_{ji} \cdot \nabla W_{ij}}

となります。SPH法では物理量の発散を

\nabla \cdot \boldsymbol{G}_i = \sum_j \frac{m_j}{\rho_j} (\boldsymbol{G}_{j} - \boldsymbol{G}_{i}) \cdot \nabla W_{ij}

のように計算します (他の発散モデルもあります)。したがって、$A$ の分母は

\sum_j \frac{m_j}{\rho_j} \boldsymbol{r}_{ji} \cdot \nabla W_{ij} = \nabla \cdot \boldsymbol{r}_i = d  \tag{7}

とすることができ、$A = 2$ になります。
重み関数を式(6)としてMPS法のラプラシアンモデルに代入すると、式(5)のSPH法のラプラシアンモデルが導出されました。

係数 $A$ が粒子間の距離に依存しない定数になるのは、圧縮性流体のシミュレーションのために作られたSPH法にとって非常に都合が良いことです。
MPS法では式(4)の係数 $A$ を初期粒子配置から計算し、そのまま使い続けるのですが、そうすると圧縮・膨張によって粒子間距離が大きく変化したときに係数 $A$ に大きな誤差が生じてしまいます。したがって、MPS法のラプラシアンモデルは粒子間の距離がそれほど変化しない非圧縮性流体のシミュレーションにのみ使用できることがわかります。
SPH法のラプラシアンモデルではこの問題は発生しません。となると、SPH法のラプラシアンモデルのほうが粒子配置に起因する誤差に強いのかと思うかもしれませんが、式(7)の発散モデルの誤差などもあり、一概にそうとは言えません。

そんなこんなで、ラプラシアンモデルに現れる係数 $A$ に着目してみると、それぞれのスキームが生まれた背景がわかるという話でした。

参考文献

Tamai, T., Shibata, K., & Koshizuka, S. (2013). Development of the Higher-order MPS Method Using the Taylor Expansion. 20130003. https://doi.org/10.11421/jsces.2013.20130003
Price, D. J. (2012). Smoothed particle hydrodynamics and magnetohydrodynamics. Journal of Computational Physics, 231(3), 759–794. https://doi.org/10.1016/j.jcp.2010.12.011