Smoothed Particle Hydrodynamics

Smoothed Particle Hydrodynamics (SPH、SPH法)は流体力学や材料力学にて用いられる微分方程式の数値解析手法（CFD）の一つ。

対象となる１個の連続体を有限個の粒子の集団に置き換えて計算する粒子法の、その代表的なもの。連続体を計算粒子集団に置き換え、連続体の支配方程式を元に導かれる粒子間相互作用を粒子集団に与えて移動させる。

Smoothed-Particle（SP、滑らか粒子）とは、中心で最大、周辺に向けてゼロへ連続変化する滑らかな密度分布をもつ粒子である。このSPはオーバーラップするように配置され、個々のSPの密度分布の重ね合わせで連続体の密度場や他の物理量の分布を表すことが根本のコンセプトである。

いわゆるフルラグランジュと言われるタイプの手法であり、基本的に移流項がなく、保存性がよい、メッシュフリー近似を実現できる、といったメリットがあり、物体や界面の大変形/大移動をともなう対象では有効である。

天体物理学、構造力学（主に破壊や亀裂進展など）、流体力学など多くの分野で利用されている。しかし、メッシュ手法に比べれば適用数は桁違いに少ない。計算妥当性や境界条件の扱い方など、十分に確立されていない部分もある。

SPHは1977年のGingoldとMonaghanによる天体物理のシミュレーションに関する論文^[1]が起源とされており、手法の名前もこの論文が元となっている。しかし、この論文の謝辞でも述べているようにSPHの離散化はLucyのアイディアであり、Lucyも同年に少し遅れて同様の離散化を用いた論文を公開している^[2]。

SPHでは粒子１個の密度分布として補間関数を定義する。補間関数は伝統的にカーネルと呼ばれる。このカーネルという語は離散データの補間関数の呼び名として一般的に用いられる。

SPHでは、このカーネルの重ね合わせで物理量を置き換える。すなわち物理量分布はカーネルの足し合わせで表され、空間微係数はカーネルの微係数の足し合わせで表される。この空間離散化コンセプトがSPHのベースである。

　SPHカーネル ${\displaystyle w_{h}}$ は例えば次のような条件を満たす。（ここでxとyは同次元であることに注意、修正求む）

${\displaystyle {\begin{alignedat}{3}&w_{h}(r){\begin{cases}>0,&0\leq r<h,\\=0,&r\geq h,\end{cases}}\qquad \qquad &{\rm {(1)}}\\&\int _{\mathbb {R} ^{d}}w_{h}(|x|)dx=1.\qquad &{\rm {(2)}}\end{alignedat}}}$

hはカーネルのサイズを決める変数でありこれがゼロに近い極限を考えるとき、その重み関数 ${\displaystyle w_{h}}$ を用いて、任意の関数 ${\displaystyle \phi }$ を

${\displaystyle \phi (x)\approx \int _{\Omega }\phi (y)w_{h}(|x-y|)dy\qquad \qquad {\rm {(3)}}}$

と無限個のカーネル足し合わせに置き換えられる。（ここで、${\displaystyle \Omega }$ は対象領域であるがあまり意味はない）。

(3)に則り、スカラ関数${\displaystyle \phi }$ に対するhを現実的な値としたときの有限個のSPHカーネル重ね合わせによる近似を

${\displaystyle \langle \phi (x)\rangle =\sum _{i=1}^{N}{\frac {m_{i}}{\rho _{i}}}\phi (x_{i})w_{h}(|x-x_{i}|)\qquad \qquad {\rm {(5)}}}$

と定義する。

また、${\displaystyle \Omega }$ 内に有限の${\displaystyle N}$ 個の粒子${\displaystyle x_{i}}$ を配置したとき、関数${\displaystyle \phi }$ の${\displaystyle \Omega }$ 内での積算は領域内粒子の個々の値の総和で表される。すなわち、

${\displaystyle \int _{\Omega }\phi (y)dy\approx \sum _{i=1}^{N}{\frac {m_{i}}{\rho _{i}}}\phi (x_{i})\qquad \qquad {\rm {(4)}}}$

と近似される。ここに、${\displaystyle m_{i}}$ と${\displaystyle \rho _{i}}$ はそれぞれ、粒子${\displaystyle x_{i}}$ が代表する質量と密度である。

勾配作用素${\displaystyle \nabla }$ （１階の微分作用素）は、カーネルの勾配の積算で近似される。

近似関数 ${\displaystyle \langle \phi (x)\rangle }$ を微分して

${\displaystyle \langle \nabla \phi (x)\rangle =\sum _{i=1}^{N}{\frac {m_{i}}{\rho _{i}}}\phi (x_{i})\nabla w_{h}(|x-x_{i}|)\qquad \qquad {\rm {(6)}}}$

と定義される。

また、${\displaystyle \nabla \phi }$ が

${\displaystyle \nabla \phi (x)\approx \int _{\Omega }\{\phi (y)-\phi (x)\}\nabla w_{h}(|x-y|)dy\qquad \qquad {\rm {(7)}}}$

と積分式で近似できることから、(4)の近似を用い、すなわち差分の重ね合わせによる置き換えもある。

${\displaystyle \langle \nabla \phi (x)\rangle =\sum _{i=1}^{N}{\frac {m_{i}}{\rho _{i}}}\{\phi (x_{i})-\phi (x)\}\nabla w_{h}(|x-x_{i}|)\qquad \qquad {\rm {(8)}}}$

　拡散を表す微分作用素であるラプラシアン${\displaystyle \Delta }$ は、

${\displaystyle \Delta \phi (x)\approx 2\int _{\Omega }{\frac {\phi (x)-\phi (y)}{|x-y|^{2}}}(x-y)\cdot \nabla w_{h}(|x-y|)dy\qquad \qquad {\rm {(9)}}}$

と近似できることから、

${\displaystyle \langle \Delta \phi (x)\rangle =2\sum _{i=1}^{N}{\frac {m_{i}}{\rho _{i}}}{\frac {\phi (x)-\phi (x_{i})}{|x-x_{i}|^{2}}}(x-x_{i})\cdot \nabla w_{h}(|x-x_{i}|)\qquad \qquad {\rm {(10)}}}$

と定義する。他にカーネル重み付け平均に基づいて定義する方法もある。

　弱圧縮性SPH(Weakly compressible SPH)の場合、エネルギー保存の考察により、次のような近似作用素が用いられる。連続の式に現れる ${\displaystyle \rho \nabla \cdot u}$ ( ${\displaystyle u}$ は流速)は

${\displaystyle \langle \rho \nabla \cdot u(x_{i})\rangle =\sum _{j=1}^{N}m_{j}(u(x_{j})-u(x_{i}))\cdot \nabla w_{h}(|x_{i}-x_{j}|)\qquad \qquad {\rm {(11)}}}$

と定義され、運動方程式に現れる圧力勾配項${\displaystyle \rho ^{-1}\nabla p}$ ( ${\displaystyle p}$ は圧力)は

${\displaystyle \langle \rho ^{-1}\nabla p(x_{i})\rangle =\sum _{j=1}^{N}m_{i}\left({\frac {p(x_{j})}{\rho (x_{j})^{2}}}+{\frac {p(x_{i})}{\rho (x_{i})^{2}}}\right)\nabla w_{h}(|x_{i}-x_{j}|)\qquad \qquad {\rm {(12)}}}$

と定義される。 また、拡散率が場によって異なる場合の拡散${\displaystyle \nabla \cdot \mu \nabla \phi }$ (${\displaystyle \mu }$ は場の拡散率を表す関数)に対する近似は

${\displaystyle \langle \nabla \cdot \mu \nabla u(x)\rangle =\sum _{i=1}^{N}{\frac {m_{i}}{\rho _{i}}}{\frac {(\mu (x)+\mu (x_{i}))(u(x)-u(x_{i}))}{|x-x_{i}|^{2}}}(x-x_{i})\cdot \nabla w_{h}(|x-x_{i}|)\qquad \qquad {\rm {(13)}}}$

または、

${\displaystyle \langle \nabla \cdot \mu \nabla u(x)\rangle =2\sum _{i=1}^{N}{\frac {m_{i}}{\rho _{i}}}{\frac {\mu (x_{i})(u(x)-u(x_{i}))}{|x-x_{i}|^{2}}}(x-x_{i})\cdot \nabla w_{h}(|x-x_{i}|)\qquad \qquad {\rm {(14)}}}$

と定義される^[3]。

SPH法は元来は単純な陽解法であり密度変動をともなう圧縮性流れ場のための手法であった。その後CFDソルバとして対象を広げ、非圧縮流れへの種々の拡張が試みられている。

特に、行列式ソルバを利用してポアソン式を反復解法で解き圧力場と同時に非圧縮速度場を得るタイプのものがISPHと呼ばれることがある。他に粒子移動を反復して速度場を非圧縮場に修正するPCISPH法がある。

ISPHと呼ばれるスキームの多くは、射影法を適用したものである。圧力勾配項と他の項とを分けて、速度発散ゼロの回転成分を取り出すという、差分法などでは古典的なアイディアを粒子法に適用したもの。粒子法の場合の最初の適用として1996年のMoving Particle Semi-implicit (MPS、MPS法)^[4]があり、その後、SPH空間離散化と組み合わせられた^[5]。

粒子法の場合は原則的に粒子移動後の粒子分布密度を元に移動前の速度を修正することになる。必然的に、いわゆる日本国内でフラクショナルステップ法と言われる、圧力のみ陰的に1ステップ先を見る方式となる。

ただし粒子法の速度発散と密度変動は定義が一通りではなく、粒子モデルに由来する変動が顕著に含まれるため、厳密に密度変動ゼロないしは速度発散ゼロを追求せず、ソースタームに0.3~0.7といった程度の緩和係数が乗じられることがある。また、非圧縮条件を近似的にでも満たすのは粒子移動前の速度場か移動後の粒子分布密度のどちらか一方だけとなるのが標準的である。

• 数値解析
• 粒子法
• メッシュフリー法
• 差分法
• 有限要素法
• 数値流体力学

• R.A., Gingold; J.J. Monaghan (1977), “Smoothed particle hydrodynamics: theory and application to non-spherical stars”, Mon. Not. R. Astron. Soc. 181: 375–389 

• L.B., Lucy (1977), “A numerical approach to the testing of the fission hypothesis”, Astron. J. 82: 1013–1024 
• J.P., Morris; P. J. Fox, and Y. Zhu (1997), “Modeling low Reynolds number incompressible flows using SPH”, J. Comput. Phys, 136: 214-226 
• S., Koshizuka; Y. Oka (1996), “Moving-particle semi-implicit method for fragmentation of incompressible fluid”, Nuclear Sci. Eng. 123: 421-434 
• S.J., Cummins; M. Rudman (1999), “An SPH projection method”, J. Comput. Phys. 152: 584-607