粒子法を使った解析
By M.Sato
1. 概要
連続体を解析するための手法として、格子法（Finite Volume Methodなど）が広く使われているが、 粒子法は格子法とは異なるアプローチで連続体を表現し、 特に自由表面の激しい変動や大規模な変形を伴う現象のシミュレーションに適している。
ここでは主に代表的な手法であるSPH法とMPS法の理論的な背景、検証結果について解説する。
粒子法を使った解析結果のポスト処理については以下の記事が詳しいので参考にしてほしい。
👉 SPH法の可視化
2. 粒子法の特徴
粒子法は、流体や固体などの連続体を、格子（メッシュ）で分けるのではなく、 「粒子の集まり」として表現して計算するシミュレーション手法である。 以下にその特徴を簡潔にまとめる。
-
格子法との違い： 従来のシミュレーション（格子法）では、空間を固定されたメッシュで分割する必要がある。 一方、粒子法では計算対象そのものを動く点（粒子）として扱う。
- 格子法： 空間を固定し、その中を流体が通り過ぎる様子を計算する（オイラー座標系に基づく計算）。
- 粒子法： 粒子一つ一つが物理量（質量、速度、圧力など）を持って空間を移動する（ラグランジュ座標系に基づく計算）。
-
主なメリット
- 自由表面の扱いが得意：水しぶきが飛ぶ、波が崩れる、液体が混ざるといった、形が激しく変化する現象を容易に扱うことができる。
- メッシュ作成が不要：複雑な形状の容器や部品の間を流れる計算でも、面倒なメッシュを切る作業が不要である。
-
主なデメリット
- 計算コストが大きい：粒子数が多いと計算時間が長くなる。
- 境界条件の取り扱いが難しい：特に固体境界との相互作用を正確に扱うのが困難である。
-
代表的な手法：現在、主に以下の2つの手法が広く使われている。
- SPH法 (Smoothed Particle Hydrodynamics)： 天体物理学から発展した手法。主に圧縮性流体や高速衝突、樹脂の成形などに使われる。
- MPS法 (Moving Particle Semi-implicit)： 日本で開発された手法。水などの縮まない液体（非圧縮性流体）の計算に特化しており、 土木や造船、自動車のオイル挙動解析などで多用される。
-
オープンソースソフトウェア：代表的なオープンソースソフトウェアには以下のようなものがある。
| 手法 | ソフトウェア |
|---|---|
| SPH法 | DualSPHysics, SPlisHSPlasH, TrixiParticles.jl, PySPH |
| MPS法 | openMPS |
3. SPH法
3.1. カーネル関数
SPH法では、カーネル関数と呼ばれる重み関数を用いて、粒子の物理量を平滑化して連続体として近似する。 粒子は質量を持ち、位置や速度などの物理量を持って空間を移動する。 粒子間の相互作用は、カーネル関数を通じて計算される。
解析領域を$\Omega$とする。領域$\Omega$上で定義される任意の関数$\phi$は、領域内の 任意点$\mathbf{x}$において、次のように表される。
$$ \phi(\mathbf{x}) = \int_{\Omega} \phi(\mathbf{\xi}) \delta(\mathbf{\xi}-\mathbf{x}) d\mathbf{\xi} $$
ここで、$\delta$はディラックのデルタ関数である。 デルタ関数は、$\mathbf{\xi}=\mathbf{x}$のときに無限大の値をとり、それ以外のときには$0$の値をとる関数である。 SPH法ではデルタ関数を、カーネル関数$w(\mathbf{\xi}-\mathbf{x},h)$で近似する。
$$ \phi(\mathbf{x}) \approx \int_{\Omega} \phi(\mathbf{\xi}) w(|\mathbf{\xi}-\mathbf{x}|) d\mathbf{\xi} $$
ここで、$h$は平滑化長と呼ばれるパラメータである。 カーネル関数は、$\mathbf{\xi}$が$\mathbf{x}$から離れるにつれて値が減少するような関数である。 カーネル関数として、スプライン関数などが使われる。 スプライン関数の場合、半径$r$に対して、区分的に関数式が定義 されることになるため、連続的な関数として表現可能なWendland関数が使われることもある。
以下に二次元領域におけるカーネル関数の例を示す。
図1 — SPH法で使われるカーネル関数（2次元）
詳細は後述するが、SPH法の計算では、カーネル関数の空間微分（一階微分、二階微分）が必要となる。 カーネル関数の一階微分、二階微分を下図に示す。
図2 — SPHカーネル関数の一階微分
図3 — SPHカーネル関数の二階微分
3.2. 空間の離散化
SPH法では離散的な粒子が持っている物理量を、カーネル関数を用いて平滑化することで、 連続体として近似する。例えば、粒子$i$が持っている物理量$\phi(\mathbf{r}_i)$は、次のように表される。
$$ \begin{align} \phi(\mathbf{r}_i) = \sum_j \phi(\mathbf{r}_j) w(|\mathbf{r}_j-\mathbf{r}_i|) V_j = \sum_j \phi(\mathbf{r}_j) w(|\mathbf{r}_j-\mathbf{r}_i|) \frac{m_j}{\rho_j} \notag \end{align} $$
ここで、$j$は粒子のインデックス、$V_j$は粒子$j$の体積、$m_j$は粒子$j$の質量、$\rho_j$は粒子$j$の密度である。
図4 — カーネル関数による粒子の離散化
次に物理量の勾配を求めてみる。$\nabla_\xi \phi(\mathbf{\xi})$をカーネル関数を用いて近似することを考える。
$$ \begin{align} \nabla \phi(\mathbf{r}_i) &= \int_{\Omega} \nabla_\xi \phi(\mathbf{\xi}) w(|\mathbf{\xi}-\mathbf{r}_i|) d\mathbf{\xi} \notag \\ &= \int_{\Omega} \nabla_\xi \left[ \phi(\mathbf{\xi}) w(|\mathbf{\xi}-\mathbf{r}_i|) \right] d\mathbf{\xi} - \int_{\Omega} \phi(\mathbf{\xi}) \nabla_\xi w(|\mathbf{\xi}-\mathbf{r}_i|) d\mathbf{\xi} \notag \\ &= \int_{\Gamma} \phi(\mathbf{\xi}) w(|\mathbf{\xi}-\mathbf{r}_i|) \hat{n} dS - \int_{\Omega} \phi(\mathbf{\xi}) \nabla_\xi w(|\mathbf{\xi}-\mathbf{r}_i|) d\mathbf{\xi} \notag \\ &= - \int_{\Omega} \phi(\mathbf{\xi}) \nabla_\xi w(|\mathbf{\xi}-\mathbf{r}_i|) d\mathbf{\xi} \notag \\ &\approx -\sum_j \phi(\mathbf{r}_j) \nabla_j w(|\mathbf{r}_j-\mathbf{r}_i|) V_j \notag \\ &= \sum_j \phi(\mathbf{r}_j) \nabla_i w(|\mathbf{r}_i-\mathbf{r}_j|) \frac{m_j}{\rho_j} \notag \end{align} $$
$w$は既知の関数であり、物理量$\phi$の空間勾配は、$\phi$と$\nabla w$の線形和で表されることになる。
なお、ここでは、$\nabla \phi$を計算する一例を示したが、この手法は実際にはあまり使われない。 勾配を計算する手法として他にも$\phi$関数をテーラー展開して近似する方法もある。 以下に、SPH法で最も一般的に使われる関数の内挿近似、勾配モデル、ラプラシアンモデルを示す。
$$ \begin{align} \langle \phi_i \rangle^{SPH} &= \sum_j \phi_j w_{ij} \frac{m_j}{\rho_j} \notag \\ \langle \nabla \phi_i \rangle^{SPH} &= \sum_j (\phi_j-\phi_i) \nabla w_{ij} \frac{m_j}{\rho_j} \notag \\ \langle \nabla^2 \phi_i \rangle^{SPH} &= 2 \sum_j (\phi_j-\phi_i) \frac{(\mathbf{r}_j-\mathbf{r}_i) \cdot \nabla w_{ij}}{|\mathbf{r}_j-\mathbf{r}_i|^2} \frac{m_j}{\rho_j} \notag \end{align} $$
ここで $w_{ij}$は$w(|\mathbf{r}_j-\mathbf{r}_i|)$、$\nabla w_{ij}$は$\nabla_i w(|\mathbf{r}_i-\mathbf{r}_j|)$である。
3.3. 計算手法
SPH法は圧縮性流体を陽解法で解く手法として開発されたが、その後、非圧縮性流体の計算に適用する ために多数の解析モデルが提案されている。ここではいくつかの代表的な手法を紹介する。
非圧縮性のナビエ・ストークス方程式は次式で表される。
$$ \begin{align} \frac{D \mathbf{u}}{D t} &= -\frac{\nabla p}{\rho} + \nu \nabla^2 \mathbf{u} + \mathbf{g} \notag \end{align} $$
ここで、$\rho$は密度、$\mathbf{u}$は流速、$p$は圧力、$\nu$は動粘度、$\mathbf{g}$は重力加速度である。 左辺は物質微分であり、流体の粒子が移動する際の時間変化を表す。 粒子法では物質の移動（流れ）をラグランジュ記述で表すため、格子法で必要となる移流項を表す必要はない。
💡 粒子法ではラグランジュ記述を採用するため、ナビエ・ストークス方程式の移流項を扱う必要がなく、移流項に起因する数値拡散の問題を回避できる。
ナビエ・ストークス方程式には未知変数として、流速$\mathbf{u}$、圧力$p$が含まれているため、 さらに連続方程式を解く必要がある。
$$ \begin{align} \frac{D \rho}{D t} &= -\rho \nabla \cdot \mathbf{u} \notag \end{align} $$
SPH法で計算を行う場合、時間積分法と圧力の計算方法からいくつかの解析手法がある。 まずこれらの解析手法で使われる粒子離散化モデルを記述しておく。
$$ \begin{align} \langle \rho_i \rangle^{SPH} &= \sum_j m_j w_{ij} \notag \\ \langle \frac{\nabla p_i}{\rho_i} \rangle^{SPH} &= \sum_j m_j \left( \frac{p_i}{\rho_i^2}+\frac{p_j}{\rho_j^2} \right) \nabla w_{ij} \notag \\ \langle \nabla^2 \mathbf{u}_i \rangle^{SPH} &= 2 \sum_j m_j \left( \frac{\mathbf{u}_j-\mathbf{u}_i}{\rho_j} \right) \frac{(\mathbf{r}_j-\mathbf{r}_i) \cdot \nabla w_{ij}}{|\mathbf{r}_j-\mathbf{r}_i|^2} \notag \end{align} $$
3.3.1. 微圧縮性SPH法（WCSPH）
WCSPH (Weakly Compressible SPH) 法は流体の微小な圧縮性を許容し、陽解法を使って密度を計算する。
$$ \begin{align} \rho_i^{n+1} &= \rho_i^n + \frac{D\rho}{Dt}\Delta t = \rho_i^n - \rho_i^n \langle \nabla \cdot \mathbf{u}_i^n \rangle^{SPH} \Delta t \notag \end{align} $$
方程式の個数（4個）に対して未知数（u,v,w,p,$\rho$）の数（5個）が多いため、 さらに状態方程式を用いて密度から圧力を計算する必要がある。 テイト方程式を使うと状態方程式は次式で表される。
$$ \begin{align} p^n = B \left[ \left(\frac{\rho^n}{\rho_0}\right)^\gamma - 1 \right] \notag \end{align} $$
$\gamma$は一般的に7が使われる。$B$は圧力スケーリングパラメータである。 $B$は、流体の音速を$c_s$とすると、次式で表される。
$$ \begin{align} B = \frac{\rho_0 c_s^2}{\gamma} \notag \end{align} $$
上記の方程式を陽解法で計算する方法が微圧縮性SPH法（WCSPH）である。 この計算方法はクーラン条件により時間刻み幅が制限される。
3.3.2. 非圧縮性SPH法（ISPH）
ISPH (Incompressible SPH) 法では、速度については陽的に時間積分し、 圧力を連続方程式を満たすように陰的に求める手法である。 これは格子法における非圧縮流体の解析で使用されるプロジェクション（射影）法をSPH法に適用したものと考えることができる。
ナビエ・ストークス方程式の左辺項を以下のように近似する。
$$ \begin{align} \left( \frac{D \mathbf{u}}{D t} \right)_i^n \approx \frac{\mathbf{u}_i^{n+1}-\mathbf{u}_i^n}{\Delta t} = \frac{\mathbf{u}_i^{n+1}-\mathbf{u}_i^{*}}{\Delta t} +\frac{\mathbf{u}_i^\ast-\mathbf{u}_i^n}{\Delta t} \notag \end{align} $$
ここで $\mathbf{u}_i^*$は、圧力項を除いたナビエ・ストークス方程式の右辺を陽的に計算した仮速度であり次式で計算される。
$$ \begin{align} \mathbf{u}_i^* = \mathbf{u}_i^n + \Delta t \left( \nu \langle \nabla^2 \mathbf{u}_i^n \rangle^{SPH} + \mathbf{g} \right) \notag \end{align} $$
この仮速度$\mathbf{u}_i^*$は、非圧縮性条件を満たさないため、 連続方程式を使い非圧縮性条件を満たすための圧力を求める。
$$ \begin{align} \mathbf{u}_i^{n+1} = \mathbf{u}_i^* - \Delta t \left( \frac{\nabla p_i^{n+1}}{\rho_0} \right) \notag \end{align} $$
両辺の発散を取り、$\nabla\cdot\mathbf{u}_i^{n+1}=0$と仮定すると、次式が得られる。
$$ \begin{align} \nabla^2 p_i^{n+1} = \rho_0 \frac{ \nabla\cdot\mathbf{u}_i^* }{\Delta t} \notag \end{align} $$
粒子を使った空間離散化を行うと、次式が得られる。
$$ \begin{align} 2\sum_j m_j \frac{p_j^{n+1}-p_i^{n+1}}{\rho_j} \frac{(\mathbf{r}_j-\mathbf{r}_i) \cdot \nabla w_{ij}}{|\mathbf{r}_j-\mathbf{r}_i|^2} = \frac{ \sum_j m_j (\mathbf{u}_j^\ast-\mathbf{u}_i^\ast) \cdot \nabla w_{ij} }{\Delta t} \notag \end{align} $$
この連立方程式を解くことで圧力を求めることができる。 次に次式で速度を更新する。
$$ \begin{align} \mathbf{u}_i^{n+1} = \mathbf{u}_i^* - \Delta t \langle \frac{1}{\rho_0} \nabla p_i^{n+1} \rangle^{SPH} \notag \end{align} $$
次に次式で粒子位置を更新する。
$$ \begin{align} \mathbf{r}_i^{n+1} = \mathbf{r}_i^n + \Delta t \mathbf{u}_i^{n+1} \notag \end{align} $$
ここで示した半陰解法は、さまざまな改良モデル（例えばIISPH (Implicit Incompressible SPH) 法など）が提案されている。 例えば、圧力のポアソン方程式の右辺項を次式で表すモデルもある。
$$ \begin{align} \nabla^2 p_i^{n+1} = \frac{\rho_0 - \langle \rho_i^* \rangle}{\Delta t^2} \notag \end{align} $$
3.3.3. エントロピ減衰人工圧縮性法 (EDAC)
EDAC法は、流体の連続方程式を解く代わりに、次の圧力の時間発展方程式を解くことで非圧縮性条件を満たす方法である。
$$ \begin{align} \frac{D p}{D t} = -\rho c_s^2 \nabla \cdot \mathbf{u} + \nu_p \nabla^2 p \notag \end{align} $$
この手法は圧力場に発生する高周波ノイズを減衰させることができ、また圧力は陽的に求めるため、 ISPH法などの半陰解法と比較して、低コストで計算を行うことができる。
4. MPS法
4.1. カーネル関数
MPS法はSPH法とは異なり、非圧縮性流体を対象として開発された手法である。 MPS法では、SPH法と同様にカーネル関数を使い物理量の空間分布を表現する。 ただしSPH法におけるカーネル関数とは異なり粒子間の相互作用を重み付き平均することで空間分布を求める。
MPS法で使われるカーネル関数は、次式で表される。
$$ \begin{align} w(r) &= \frac{r_e}{r} - 1 & 0 < r < r_e \notag \\ &= 0 & r \ge r_e \notag \end{align} $$
ここで、$r$は粒子間の距離、$r_e$は相互作用半径である。
図5 — MPS法で使われるカーネル関数
4.2. 計算手法
MPS法で一般的に使用される関数の内挿近似、勾配モデル、ラプラシアンモデルを示す。
$$ \begin{align} \langle \phi_i \rangle^{MPS} &= \sum_{j\ne i} \phi_j w(|\mathbf{r}_j-\mathbf{r}_i|) \notag \\ \langle \nabla \phi_i \rangle^{MPS} &= \frac{D}{n_0} \sum_{j\ne i} (\phi_j-\phi_i) \frac{\mathbf{r}_j-\mathbf{r}_i}{|\mathbf{r}_j-\mathbf{r}_i|^2} w(|\mathbf{r}_j-\mathbf{r}_i|) \notag \\ \langle \nabla^2 \phi_i \rangle^{MPS} &= \frac{2D}{\lambda n_0} \sum_{j\ne i} (\phi_j-\phi_i) w(|\mathbf{r}_j-\mathbf{r}_i|) \notag \end{align} $$
ここで$n_0$は粒子の数密度、$D$は次元（1,2,3のいずれか）、$\lambda$はスケーリングパラメータである。 SPH法とは異なり、MPS法ではカーネル関数の微分は使われない。
$$ \begin{align} n_0 &= \sum_{j\ne i} w(|\mathbf{r}_j-\mathbf{r}_i|) \notag \\ \lambda &= \frac{\sum_{j\ne i} |\mathbf{r}_j-\mathbf{r}_i|^2 w(|\mathbf{r}_j-\mathbf{r}_i|)}{\sum_{j\ne i} w(|\mathbf{r}_j-\mathbf{r}_i|)} \notag \end{align} $$
MPS法では、半陰解法を使い、圧力のポアソン方程式を解くことで圧力を求める。 計算アルゴリズムはSPH法のISPH法と同様であるのでここでは詳細を記述しない。
5. 計算例
ここでは、SPH法を対象として複数の計算モデルを比較した結果を示す。
なおここでは、PySPHというPythonで書かれたSPH法のオープンソースコードを使って計算を行った。
PySPHには多数の例題が同梱されている。例題を実行する際に
コマンドライン引数に--schemeオプションを指定することで計算スキームを指定できる。
5.1. dam break (2D)
ダムブレイクは、流体が重力の影響を受けて自由落下する現象である。 ダムブレイクの解析は、流体の運動や衝突、波の伝播などを理解するために重要な問題である。
図6 — ダムブレイク問題の初期条件
PySPHに同梱されている例題（dam_break_2d.py）を使い計算を行った。
以下に、WCSPH法、IISPH法、EDAC法を使った場合の流速分布のアニメーションを示す。
IISPH法の計算では水滴が細かく飛び散る様子が見られることが分かった。
動画1 — WCSPH法によるダムブレイク解析
動画2 — IISPH法によるダムブレイク解析
動画3 — EDAC法によるダムブレイク解析
5.2. テイラーグリーン渦
テイラーグリーン渦は、非圧縮性流体の解析でよく使われるベンチマーク問題である。 テイラーグリーン渦の理論解は次式で表される。
$$ \begin{align} u &= -U \cos(2\pi x)\sin(2\pi y) e^{b t} \notag \\ v &= U \sin(2\pi x)\cos(2\pi y) e^{b t} \notag \\ p &= -\frac{U^2}{4}(\cos(4\pi x)+\cos(4\pi y)) e^{2b t} \notag \end{align} $$
計算条件は以下の通りである。
| パラメータ | 値 |
|---|---|
| 代表流速 $U$ | $1$ [m/s] |
| 代表長さ $L$ | $1$ [m] |
| 動粘度 $\nu$ | $0.01$ [m²/s] |
| 減衰係数 $b$ | $-8\pi^2/Re$（$Re=UL/\nu$） |
最大流速の時間変化は次式で表される。
$$ \begin{align} u_{max} = U e^{b t} \notag \end{align} $$
PySPHに同梱されている例題（taylor_green.py）を使い計算を行った。
以下にEDAC法を使った場合の流速分布のアニメーションを示す。
動画4 — EDAC法によるテイラーグリーン渦の流速分布
以下に、最大流速の経時変化を示す。WCSPH法、IISPH法、TVF法、EDAC法、および理論解の結果を比較している。 この計算結果によればEDAC法が理論解に最も近い結果を示していることがわかる。
図7 — 最大流速の時間履歴（各SPHスキーム vs 理論解）
💡 テイラーグリーン渦の減衰挙動を見ると、EDAC法が理論解に最も近い結果を示している。EDAC法は陽的に圧力を解くため計算コストも低く、非圧縮性流れの解析において有力な選択肢となる。
テイラーグリーン渦は時間と共に流速が減衰していく。 初期時刻（t=0）と最終時刻（t=2）における流速の分布を以下に示す。 なお色分布はそれぞれの時刻での最大最小値で正規化したものになっている。
最終時刻でのコンター図を見ると、EDAC法は理論解に類似した流速分布、圧力分布になっていることがわかる （すなわち初期値と比較すると強度は異なるが相似な空間分布形状）。
図8 — 初期時刻 (t=0) における流速・圧力分布
図9 — 最終時刻 (t=2) におけるWCSPH法・EDAC法の流速・圧力分布
図10 — 最終時刻 (t=2) におけるIISPH法・TVF法の流速・圧力分布
6. まとめ
粒子法を使った連続体の解析は、格子法とは異なるアプローチで連続体を表現し、 特に自由表面の激しい変動や大規模な変形を伴う現象のシミュレーションに適している。
粒子法の代表的な手法であるSPH法とMPS法の理論的な背景、検証結果について解説した。 SPH法は圧縮性流体を陽解法で解く手法として開発されたが、 その後、非圧縮性流体の計算に適用するために多数の解析モデルが提案されている。 MPS法は非圧縮性流体を対象として開発された手法である。
SPH法を対象として複数の計算モデル（WCSPH, IISPH, TVF, EDAC）を比較した結果、 ダムブレイク・テイラーグリーン渦のいずれにおいてもスキームによる挙動の違いが顕著に現れ、 EDAC法が理論解に最も近い結果を示すことが確認できた。
弊社の解析事例
弊社の流体解析事例については、下記のリンクからご覧ください。