概要

　土木の分野において、粒子法の研究が盛んに行われております。特に土石流や津波、液状化現象などのシミュレーションで活用されています。今回の記事では、安定化完全陰解法非圧縮SPH（Smoothed Particle Hydrodynamics）を紹介したいと思います。

　本記事は以下の書籍の内容を参考にしています。

支配方程式

ラグランジュ記述の非圧縮ナビエ・ストークス方程式は次式で表せる。

\begin{align*}
\frac{D\mathbf{v}}{Dt} &= -\frac{1}{\rho}\nabla p + \frac{\mu}{\rho}\nabla^2\mathbf{v} + \mathbf{g} \\
\frac{D\rho}{Dt} &= -\rho \nabla \cdot \mathbf{v} = 0
\end{align*}

ここで、

$\mathbf{v}$ : 速度ベクトル
$\rho$ : 密度
$p$ : 圧力
$\mu$ : 粘性係数
$\mathbf{g}$ : 重力加速度ベクトル

第1式が運動方程式、第2式が連続の式（非圧縮条件）を表す。今回の記事では、上記の方程式をSPHを使い離散化する。

SPHの簡単な説明

　SPH法に関する記事はたくさんあり、詳細は他の記事に譲り、ここでは本記事に必要な部分を簡潔に説明する。

　偏微分方程式をコンピューターで解く際は、微分や積分を厳密に表現できないので離散化を行う必要がある。離散化とは微分や積分を四則演算（+-×÷）で表すことである。

　SPHでは、デルタ関数を重み関数$W$で近似的に表現し、物理量$\phi$を次式で示す。　

\begin{align*}
 \phi(\mathbf{x}) = \int_{\Omega}d^3\xi\phi(\mathbf{\xi}) \delta(\mathbf{\xi}-\mathbf{x}) \approx \int_{\Omega}d^3\xi \phi(\mathbf{\xi}) W(|\mathbf{\xi}-\mathbf{x}|)
\end{align*}

重み関数$W$は以下を満たさなければならない。

ユニティ条件
偶関数特性
コンパクトサポート
非負性条件
デルタ関数へ収束
単調減少条件
微分連続性

　今回の記事では、以下の3次スプラインカーネル関数を用いる。

\begin{align*}
W(r,h) = \frac{1}{\pi h^3} \begin{cases}
1 - \frac{3}{2}(\frac{r}{h})^2 + \frac{3}{4}(\frac{r}{h})^3 & (0 \leq \frac{r}{h} \leq 1) \\
\frac{1}{4}(2-\frac{r}{h})^3 & (1 < \frac{r}{h} \leq 2) \\
0 & ( \frac{r}{h} > 2)
\end{cases}
\end{align*}

$r$ : 粒子間距離
$h$ : 影響半径

　重み関数$W$を使うことで、$\phi$に作用する微分が重み関数に作用する。

\begin{align*}
\nabla_x \phi(\mathbf{x}) = \int_{\Omega}d^3\xi \phi(\mathbf{\xi}) \nabla_x W(|\mathbf{\xi}-\mathbf{x}|) 
\end{align*}

カーネル関数の微分は解析的に求めることができる。

　次に、積分を評価するために粒子を導入する。各点に粒子を配置し、以下の式のように粒子$i$の近傍にある粒子の物理量を足し上げることで積分を表現する。

\begin{align*}
\phi(\mathbf{x}) &=\int_{\Omega}d^3\xi \phi(\mathbf{\xi}) W(|\mathbf{\xi}-\mathbf{x}|) \\
& \approx \sum_j^N V_{(j)}  \phi(\mathbf{x}_{(j)})W(|\mathbf{x}_{(j)}-\mathbf{x}|) = \sum_j^N \frac{m_{(j)}}{\rho_{(j)}}  \phi(\mathbf{x}_{(j)})W(|\mathbf{x}_{(j)}-\mathbf{x}|) \\
&\rightarrow \langle \phi(\mathbf{x}) \rangle
\end{align*}

$(i)$ は粒子の表す添字である。また、$m$は質量、$V$は体積である。離散化された量は、$ \langle \rangle$ で表現する。

　ナビエ・ストークス方程式を解いて速度が求まれば、粒子$(i)$の位置は以下のように更新される。

\begin{align*}
\mathbf{x}_{(i)} \leftarrow \mathbf{x}_{(i)} + \mathbf{v}_{(i)}\Delta t
\end{align*}

$\Delta t$ は時間を離散化した際のステップ幅である。

空間の離散化

　上記の基本的な離散化の方法では、精度や保存則の観点から不十分であることが知られている。今回の記事では以下の和モデルと呼ばれる離散化手法を用いる。

　粒子間距離を以下で表す。

 {\bf r}_{ij} =|{\bf x'}_{(j)}-{\bf x}_{(i)}|

各演算子の離散化

密度

\begin{align*}
\langle \rho_{(i)} \rangle = \sum_{j=1}^N m_{(j)} W(r_{ij} ) 
\end{align*}

圧力勾配

\begin{align*}
\left\langle \frac{\nabla P_{(i)}}{\rho_{(i)}} \right\rangle = \sum_{j=1}^N m_{(j)} \left( \frac{P_{(i)}}{\rho_{(i)}^2} +\frac{P_{(j)}}{\rho_{(j)}^2}  \right) \nabla W(r_{ij} ) 
\end{align*}

ラプラシアン

\begin{align*}
\left\langle \nabla^2 {\bf v}_{(i)} \right\rangle = 2\sum_{j=1}^N  \frac{m_{(j)}}{\rho_{(j)}}\frac{{\bf v}_{(j)} -{\bf v}_{(i)} }{r_{ij}^2} \left({\bf x}_{(j)} -{\bf x}_{(i)} \right) \cdot \nabla W(r_{ij} ) 
\end{align*}

発散

\begin{align*}
\left\langle \nabla \cdot {\bf v}_{(i)} \right\rangle = \sum_{j=1}^N  \frac{m_{(j)}}{\rho_{(j)}} \left({\bf v}_{(j)} -{\bf v}_{(i)} \right)  \cdot \nabla W(r_{ij} ) 
\end{align*}

　影響半径$h$は計算する演算子によって異なり、以下のように設定した。

ラプラシアンモデル：$h = 4.1 r_0$
その他　　　　　　：$h = 2.1 r_0$

$r_0$は初期粒子間距離である。

時間の離散化

　非圧縮ナビエ・ストークス方程式を時間離散化すると以下となる。

\begin{align*}
&{\bf v}_{(i)}^{n+1} ={\bf v}_{(i)}^n + \Delta t \left\{ -\left\langle \frac{\nabla P_{(i)}}{\rho_{(i)}} \right\rangle  + \nu \left\langle \nabla^2 {\bf v}_{(i)} \right\rangle + {\bf g} \right\} \\
&\rho_{(i)}^{n+1} -\rho_{(i)}^{n} = - \rho_{(i)}^{n} \left\langle \nabla \cdot {\bf v}_{(i)} \right\rangle \Delta t \approx 0
\end{align*}

$\nu = \mu/\rho $は動粘性係数である。

　上記の方程式からは、圧力$P$は陽に求まらない。従って、今回の記事では、射影法を採用した。この方法は以下の手順で計算を行う。

圧力勾配項を除いた方程式から仮想速度$v^*$を求める
非圧縮条件を満たすように圧力を求める
圧力勾配により速度を修正する

　他の方法として、変分マルチスケーリング法を用いた一体型解法などがある。

完全陰解法 SPH

　$n$をステップ数として、$v^n$は既知とする。仮想速度は以下の式から計算できる。

\begin{align*}
{\bf v}_{(i)}^{*} ={\bf v}_{(i)}^n + \Delta t \left\{ \nu \left\langle \nabla^2 {\bf v}^*_{(i)} \right\rangle + {\bf g} \right\}
\end{align*}

ただし、この方程式は仮想速度$v^*$に関して陰的であり、連立方程式を解く必要がある。

　便宜上、位置・速度ベクトルの成分は次式で表現する。

\begin{align*}
{\bf x}_{(i)} = (x,y,z) \ \ , \ \  {\bf v}_{(i)} = (u_{(i)},v_{(i)},w_{(i)})
\end{align*}

　各方向の仮想速度は互いに独立に計算できる。以下、各方向について連立方程式を示す。

x方向の仮想速度に関する連立方程式

\begin{align*}
A^x_{ij} u_{(j)}^* = u_{(j)}^n 
\end{align*}

\begin{align*}
A^x_{ij} = \begin{cases}
  -2\Delta t \nu \frac{m_{(j)}}{\rho_{(j)}}\frac{x_{(j)} - x_{(i)} }{r_{ij}^2}  \frac{\partial r_{ij}}{\partial x_{(i)}}\frac{\mathrm{d}  W(r_{ij} )}{\mathrm{d} r_{ij}}  & (i \neq j) \\
 1+2\Delta t \nu \sum_{j=1}^N \frac{m_{(j)}}{\rho_{(j)}}\frac{x_{(j)} - x_{(i)} }{r_{ij}^2}  \frac{\partial r_{ij}}{\partial x_{(i)}}\frac{\mathrm{d}  W(r_{ij} )}{\mathrm{d} r_{ij}} & (i=j) \\
\end{cases}
\end{align*}

y方向の仮想速度に関する連立方程式

\begin{align*}
A^y_{ij} v_{(j)}^* = v_{(j)}^n 
\end{align*}

\begin{align*}
A^y_{ij} = \begin{cases}
  -2\Delta t \nu \frac{m_{(j)}}{\rho_{(j)}}\frac{y_{(j)} - y_{(i)} }{r_{ij}^2}  \frac{\partial r_{ij}}{\partial y_{(i)}}\frac{\mathrm{d}  W(r_{ij} )}{\mathrm{d} r_{ij}}  & (i \neq j) \\
 1+2\Delta t \nu \sum_{j=1}^N \frac{m_{(j)}}{\rho_{(j)}}\frac{y_{(j)} - y_{(i)} }{r_{ij}^2}  \frac{\partial r_{ij}}{\partial y_{(i)}}\frac{\mathrm{d}  W(r_{ij} )}{\mathrm{d} r_{ij}} & (i=j) \\
\end{cases}
\end{align*}

z方向の仮想速度に関する連立方程式

\begin{align*}
A^z_{ij} w_{(j)}^* = w_{(j)}^n -g \Delta t
\end{align*}

\begin{align*}
A^z_{ij} = \begin{cases}
  -2\Delta t \nu \frac{m_{(j)}}{\rho_{(j)}}\frac{z_{(j)} - z_{(i)} }{r_{ij}^2}  \frac{\partial r_{ij}}{\partial z_{(i)}}\frac{\mathrm{d}  W(r_{ij} )}{\mathrm{d} r_{ij}}  & (i \neq j) \\
 1+2\Delta t \nu \sum_{j=1}^N \frac{m_{(j)}}{\rho_{(j)}}\frac{z_{(j)} - z_{(i)} }{r_{ij}^2}  \frac{\partial r_{ij}}{\partial z_{(i)}}\frac{\mathrm{d}  W(r_{ij} )}{\mathrm{d} r_{ij}} & (i=j) \\
\end{cases}
\end{align*}

　
　連立方程式を解く際の解法は、クリロフ部分空間を用いるならば、共役勾配法（CG法）、安定化双共役勾配法（BiCGSTAB法）、一般化最小残差法（GMRES）が挙げられる。対称行列で正定値行列ならばCG法の収束性が良いと思われる。上記の連立方程式は、対称行列ではないためCG法は使えない。今回の記事では、全粒子の質量が等しいことを利用して、以下のように変換した。

\begin{align*}
2\frac{m_{(j)}}{\rho_{(j)}}\rightarrow \frac{m_{(j)}}{\rho_{(i)}\rho_{(j)}}(\rho_{(i)}+\rho_{(j)})
\end{align*}

この変換により、係数行列が対称となりCG法が適用できる。

　粒子間距離が非常に小さくなる($r \approx 0$)と、行列要素が発散する。これを防ぐため

\begin{align*}
\frac{1}{{r_{ij}^2}}\rightarrow \frac{1}{{r_{ij}^2}+\eta^2}  \ \ , \ \ \eta = 0.01h
\end{align*}

とする。$\eta$ を省くと、以下のアニメーションのように数値不安定性が発生する。

非圧縮 SPH

　仮想速度を用いて、ナビエ・ストークス方程式を表すと

\begin{align*}
{\bf v}_{(i)}^{n+1} ={\bf v}_{(i)}^* - \Delta t \left\langle \frac{\nabla P_{(i)}}{\rho_{(i)}} \right\rangle 
\end{align*}

となる。両辺に発散を作用させると、非圧縮条件を適用すると、次式のポアソン方程式が得られる。

\begin{align*}
\left\langle \nabla \cdot {\bf v}_{(i)}^{n+1} \right\rangle & = \left\langle \nabla \cdot {\bf v}_{(i)}^* \right\rangle  - \Delta t \left\langle \frac{\nabla^2 P_{(i)}}{\rho^0} \right\rangle = 0 \\
\therefore \left\langle \nabla \cdot {\bf v}_{(i)}^* \right\rangle&= -\frac{\rho^{0} -\langle \rho_{(i)} \rangle}{\rho^{0} \Delta t} 
\end{align*}

また、密度は変化しないとして、$\rho^{0} \approx \rho_{(i)}^{n+1} $と仮定した。

　圧力$P$を求めるため、次式の連立方程式を解く。

\begin{align*}
A^P_{ij} P_{(j)} = \frac{\rho^0 - \langle \rho_{(i)} \rangle}{\Delta t^2}
\end{align*}

\begin{align*}
C_{ij} &= -2\frac{m_{(j)}}{\rho_{(j)}}\left\{ \frac{x_{(j)} - x_{(i)} }{r_{ij}^2}  \frac{\partial r_{ij}}{\partial x_{(i)}}+\frac{y_{(j)} - y_{(i)} }{r_{ij}^2}  \frac{\partial r_{ij}}{\partial y_{(i)}}+\frac{z_{(j)} - z_{(i)} }{r_{ij}^2}  \frac{\partial r_{ij}}{\partial z_{(i)}} \right\} \frac{\mathrm{d}  W(r_{ij} )}{\mathrm{d} r_{ij}} 
\end{align*}

\begin{align*}
A^P_{ij} &= \begin{cases}
  C_{ij} & (i\neq j)\\
 -\sum_{j=1}^N C_{ij} & (i=j) \\
\end{cases}
\end{align*}

である。CG法を使う際は、対称化を行う。

安定化非圧縮 SPH

　非圧縮SPHは、非圧縮条件が完全に満たされ、密度一定であることが前提である。しかし、数値計算では密度は初期密度からズレてしまい、非圧縮条件は満たされない。

\begin{align*}
\left\langle \nabla \cdot {\bf v}_{(i)}^{n+1} \right\rangle  = \left\langle \nabla \cdot {\bf v}_{(i)}^* \right\rangle  - \Delta t \left\langle \frac{\nabla^2 P_{(i)}}{\rho^0} \right\rangle \neq 0
\end{align*}

　従って、初期密度からズレを密度補正項$\Delta \rho_{(i)}$として導入した。

\begin{align*} 
\rho^{0} \approx \langle \rho_{(i)}^{n+1} \rangle = \langle \rho_{(i)}^n \rangle + \Delta \rho_{(i)} 
\end{align*}

\begin{align*} 
\Delta \rho_{(i)} = \langle \rho_{(i)}^{n+1} \rangle- \langle \rho_{(i)}^n \rangle \approx \alpha_{\epsilon} (\rho^{0} - \langle \rho_{(i)}^n \rangle)
\end{align*}

緩和パラメータ$\alpha_{\epsilon}$は十分小さな値に設定し、今回の記事では$\alpha_{\epsilon}=1.0\times10^{-4}$とした。

　密度の時間発展を考慮して、以下の方程式を解く。

\begin{align*} 
\frac{\langle \rho_{(i)}^{n+1} \rangle -\langle \rho_{(i)}^n \rangle}{ \Delta t} +\rho^{0}\left\langle \nabla \cdot {\bf v}_{(i)}^{n+1} \right\rangle = 0
\end{align*}

最終的に、圧力に関する連立方程式は次式となり、これを解けば圧力$P$が求まる。

\begin{align*}
A^P_{ij} P_{(j)} = \frac{\rho^0}{\Delta t} \left\langle \nabla \cdot {\bf v}_{(i)}^* \right\rangle + \alpha_{\epsilon} \frac{\rho^0 - \langle \rho_{(i)}^n \rangle}{\Delta t^2}
\end{align*}

自由表面の探索

　水面の圧力はゼロ（または大気圧）になるため、自由表面上の粒子を特定し、ディリクレ境界条件$P=0$を課す必要がある。

　参考にした書籍で一番良いとされているモーメント行列の方法を用いる。モーメント行列は粒子分布の偏りを表現し、以下で定義される。

\begin{align*}
{\bf M}_{(i)} = \sum_{j=1}^N  \frac{m_{(j)}}{\rho_{(j)}} \left({\bf x}_{(j)} -{\bf x}_{(i)} \right)  \otimes \nabla W(r_{ij} ) 
\end{align*}

モーメント行列の固有値は、規則正しい粒子配置なら固有値が1になる。自由表面では影響半径内の粒子が少ないため、モーメント行列の固有値は１より小さくなる。

　モーメント行列の最小固有値$\lambda_{(i)}^{\min}$から自由表面を探索する。固有値を求める際は、モーメント行列が対称行列であることを利用してヤコビ法を用いた。

\begin{align*}
\begin{cases}
&\lambda_{(i)}^{\min} < 0.2 \ \ \qquad\mbox{自由表面} \\
0.2 \leq &\lambda_{(i)}^{\min}\leq 0.75  \qquad\mbox{グレーゾーン} \\
 &\lambda_{(i)}^{\min} >0.75  \qquad\mbox{流体内部} \\
\end{cases}
\end{align*}

グレーゾーンの場合

　グレーゾーンの場合、さらに詳細な幾何学的判定を行う。以下の条件を満たす粒子jがいなければ、粒子iを自由表面と判定する。

\begin{align*}
|{\bf x}_{(j)} - {\bf x}_{(i)}| \leq \sqrt{2}H  \ \ &\mbox{かつ} \ \  \cos^{-1} ({\bf n}_{(i)} \cdot ({\bf x}_{(j)} - {\bf x}_{(i)}) ) \leq \frac{\pi}{4}　\\
|{\bf x}_{(j)} - {\bf x}_{(i)}| > \sqrt{2}H  \ \ &\mbox{かつ} \ \  |{\bf x}_{(j)} - {\bf x}_{(i')}| \leq H 
\end{align*}

\begin{align*}
H = 1.33 d
\end{align*}

　法線ベクトル${\bf n}_{(i)}$は、最小固有値の勾配から以下のように計算される。

\begin{align*}
{\bf n}_{(i)} = -\frac{\left\langle \nabla \lambda^{\min}_{(i)}\right\rangle}{\left| \left\langle \nabla \lambda^{\min}_{(i)}\right\rangle\right|} 
\end{align*}

\begin{align*}
\left\langle \nabla \lambda^{\min}_{(i)}\right\rangle  = {\bf M}_{(i)} \sum_{j=1}^N \frac{m_{(j)}}{\rho_{(j)}}\left(\lambda_{(j)}^{\min}-\lambda_{(i)}^{\min} \right)\nabla W(r_{ij})
\end{align*}

粒子間の強制斥力

　越塚らの「粒子法入門」のプログラムでは、粒子間の距離が近くなると斥力が強制的に働くよう衝突力を計算する。

\begin{align*}
\delta{\bf v}_{(i)} = (1+e)\sum_{j=1}^N\frac{\rho_{(j)}}{\rho_{(i)}+\rho_{(j)}} \frac{({\bf v}_{(j)}-{\bf v}_{(i)})\cdot ({\bf x}_{(j)}-{\bf x}_{(i)})}{|{\bf x}_{(j)}-{\bf x}_{(i)} |} \frac{{\bf x}_{(j)}-{\bf x}_{(i)}}{|{\bf x}_{(j)}-{\bf x}_{(i)} |} 
\end{align*}

$e$は反発係数である。強制斥力を用いて、速度を次式のように更新する。

\begin{align*}
{\bf v}_{(i)} \leftarrow {\bf v}_{(i)} +\delta{\bf v}_{(i)} 
\end{align*}

強制斥力は壁面との相互作用に適用することで、粒子の壁すり抜けを効果的に防止できる。

静力学平衡の理論解との比較

　SPH法の実装が正しいかを検証するため、静止流体の圧力分布を理論解と比較する。静力学平衡の式は、次式の圧力と重力の釣り合い式である。

\begin{align*}
\Delta P = -\rho g \Delta z 
\end{align*}

以下の図のように、SPHで流体を静止させたシミュレーションを行い、理論解と一致するか検討した。粒子数は$108,045$個とし、図のように壁粒子やダミー粒子を配置した。密度は $1000 \mbox{kg/}\mbox{m}^3$ 、動粘性係数は$1.0\times 10^{-6} \mbox{m}^2/s$とした。

以下の図は$30$ステップ時の圧力分布である。縦軸は圧力、横軸は表面からの深さである。オレンジ線が理論解、青プロットはSPHの数値解である。平衡状態となり理論解と一致していることが確認できる。

しかし、$500$ステップ時では解が振動する。

水柱崩壊解析（ダムブレーク）

　以下は、ダムブレークのシミュレーション結果である。粒子数は$164,688$個、２秒間の解析である。ただし、長方形の箱から飛び出た粒子は計算対象外としている。

水柱が崩壊し、壁面に衝突して飛散する様子が再現されており、自由表面の追跡や大変形にも対応できていることが確認できる。

まとめ

　本記事では、安定化完全陰解法非圧縮SPH法の理論と実装について解説した。

SPH法の基礎理論：カーネル近似と粒子近似
和モデルによる改良された空間離散化
完全陰解法による時間積分
射影法による圧力と速度の分離計算
安定化手法：密度補正項の導入
自由表面探索：モーメント行列法
理論解との比較による精度検証

今後は、土石流や雪崩などの実現象に向けたシミュレーションを行いたい。ソースコードは githubの準備ができ次第公開する。