SPH法の理論

SPH法について†

SPH(Smoothed Particle Hydrodynamics)はもともと 銀河系の衝突や天体の形成などの宇宙物理学におけるシミュレーションのために提案された手法である[Lucy1977]が， 近年流体や熱など他の多くの現象に応用されている． ただし，元々宇宙物理学の問題に多くみられる圧縮性流体を取り扱うために提案された手法であったため， 水や速度の遅い空気(流速<<音速)のような非圧縮性流体には向いていない (非圧縮性を強制する手法はいくつか考案されている[Becker2007,Solenthaler2009])． 計算時間がかかってもいいならば最初から非圧縮流体を対象として開発されたMPSの方がよい場合もある．

CG分野に初めてSPHを導入したのは[Desbrun1996]で， 3次元の粘性流体に適用したのは[Muller2003]である． 以下では，主に[Muller2003]に基づいた方法を述べる．

SPH法による離散化式†

SPH法による物理量の離散化式を以下に示す．

ここで，Nは近傍パーティクルの集合，mはパーティクル質量，はパーティクル密度， Wはカーネル関数である． カーネル関数は有効半径hを持ち，有効半径外では0となる(コンパクトサポート)． また，その積分が1となる()ように設定される *1．

上式のようにSPHでは物理量を周囲のパーティクルの重み付き和で近似する． そして，物理量の勾配はカーネル関数の導関数を用いて表す． 物理量の勾配は，

と表され，物理量のラプラシアンは，

と表される．

上式を用いることで流体の密度は以下で計算される．

シミュレーションに用いられるカーネル関数については， SPH法の重み関数 を参照．

SPH法によるNS方程式の解法†

支配方程式である非圧縮のナビエ・ストークス方程式を以下に示す．

ここで，は流体速度，は動粘性係数，は流体の密度，pは圧力，は外力で重力，表面張力などが含まれる． 上の式が質量保存式，下が運動量保存式である．

支配方程式をパーティクルで離散化し，SPH法を用いて解く． パーティクル自体が液体を表しているため，パーティクル質量が変化しないかぎり質量保存性が常に保持され，質量保存式を解く必要がない(あくまで質量保存であり，体積保存ではないことに注意)． また，粒子法ではグリッド法のような移流項の計算を行う必要がなく，粒子を移流させるだけでよい． そのため，移流項は離散化する必要がない． よって，各パーティクルiの速度を以下で更新すればよい．

の各項について以下で述べる．

粘性拡散項†

粘性項は流体の速度により構成される． 速度のラプラシアンをそのまま離散化すると(離散化式のをで置き換える)， あるパーティクルiからjに与える力とjからiへの力が異なる(asymmetric)． このままでは速度の速いパーティクルは遅いパーティクルに影響を与えるが，その逆の影響力が小さくパーティクル速度が発散してしまう．

粘性力は流体の速度差に依存すると考えて，symmetricに作用するようにしたのが以下の式である．

圧力項†

圧力項は圧力値の勾配で構成される． 勾配値用の離散化式を用いることで圧力項を離散化することができるが，粘性拡散項と同様にasymmetricの問題が出るので，[Muller2003]では，

と圧力の平均値を用いている．

の計算には圧力と密度の関係を表す状態方程式が用いられる．理想気体の状態方程式は，

ここで，kはガス定数である． 多くのSPHシミュレーションでは，[Desbrun1996]の提案した式，

が用いられている．ここで，はシミュレーションする対象の流体の密度である． 例えば，水ならばであるが，実装では初期配置での最大密度を用いてもよい(もちろん設定した密度になるように適切にhなどを設定した方がよい)．この式は単純でわかりやすいが，かなり流体が圧縮してしまうという欠点がある． 非圧縮性を強制するもっともポピュラーな方法はポアソン方程式を解くことであるが，そのためには巨大な線形システムを解く必要があり，非常に計算時間がかかる(MPS法ではこれをやっている)．

[Becker2007]ではTait方程式を用いることである程度の非圧縮性を確保している．Tait方程式は以下．

ここで，．

Bは圧力定数で，以下の式となる．

ここで，は音速である． として，の関係から算出する． 密度の変化率であるには0.01がよく用いられる(これはそのまま流体の非圧縮性を1%以下に強制することを示す)．

この方法は高い非圧縮性(1%以下)を保てるものの，安定した計算のためには，音速の大きさに依存してタイムステップ幅を非常に小さくする必要がある．タイムステップ幅はCFL(Courant-Friedrichs-Lewy)条件から，

と計算される． ここで，は外力，は粘性定数で[Becker2007]では0.08から0.5の値が用いられている． これで計算するとタイムステップ幅が0.0005[s]ぐらいになってしまうが， 各ステップの計算が速いのでMPSなどと比べても問題ない速度にはなる(リアルタイムには向かないが...)．

参考文献†

[Becker2007] M. Becker and M. Teschner, Weakly Compressible SPH for Free Furface Flows, In Proc. SCA2007, pp.209-217, 2007.

[Desbrun1996] M. Desbrun and M.-P. Cani, Smoothed Particles: A New Paradigm for Animating Highly Deformable Bodies, Eurographics Workshop on Computer Animation and Simulation (EGCAS), pp.61-76. 1996.

[Lucy1977] L. B. Lucy, A Numerical Approach to the Testing of the Fission Hypothesis, The Astronomical Journal, Vol.82, No.12, pp.1013-1024, 1977.

[Muller2003] M. Muller, D. Charypar and M. Gross, Particle-based Fluid Simulation for Interactive Applications, In Proc. SCA2003, pp.154-159, 2003.

[Solenthaler2009] B. Solenthaler and R. Pajarola, Predictive-Corrective Incompressible SPH, In Proc. SIGGRAPH2009, pp.1-6, 2009.