SPH法について

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

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

SPH法による離散化式

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

sph_eq_01.gif

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

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

sph_eq_05.gif

と表され,物理量sph_eq_phi.gifのラプラシアンは,

sph_eq_06.gif

と表される.

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

sph_eq_07.gif

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

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

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

sph_eq_08.gif
sph_eq_09.gif

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

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

sph_eq_12.gif
sph_eq_13.gif

sph_eq_Vf_i.gifの各項について以下で述べる.

粘性拡散項

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

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

sph_eq_14.gif

圧力項

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

sph_eq_15.gif

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

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

sph_eq_16a.gif

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

sph_eq_16.gif

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

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

sph_eq_19.gif

ここで,sph_eq_20.gif

Bは圧力定数で,以下の式となる.

sph_eq_21.gif

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

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

sph_eq_24.gif

と計算される. ここで,sph_eq_Vf_a.gifは外力,sph_eq_alpha.gifは粘性定数で[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.


*1 カーネル関数にしたい式sph_eq_fr.gif積分した値sph_eq_03.gifを求めることで,sph_eq_04.gifとして設定できる

添付ファイル: filesph_eq_rho.gif 1692件 [詳細] filesph_eq_Vu.gif 1567件 [詳細] filesph_eq_Vf_a.gif 1521件 [詳細] filesph_eq_p_i.gif 1423件 [詳細] filesph_eq_nabla_phi.gif 1526件 [詳細] filesph_eq_Vu_i.gif 1456件 [詳細] filesph_eq_phi.gif 1494件 [詳細] filesph_eq_Vf_ext.gif 1532件 [詳細] filesph_eq_c_s.gif 1471件 [詳細] filesph_eq_alpha.gif 1479件 [詳細] filesph_eq_fr.gif 1185件 [詳細] filesph_eq_Vf_i.gif 1418件 [詳細] filesph_eq_eta.gif 1293件 [詳細] filesph_eq_rho_0.gif 1362件 [詳細] filesph_eq_13.gif 1408件 [詳細] filesph_eq_22.gif 1354件 [詳細] filesph_eq_16.gif 1427件 [詳細] filesph_eq_17.gif 1482件 [詳細] filesph_eq_09.gif 1864件 [詳細] filesph_eq_24.gif 1798件 [詳細] filesph_eq_01.gif 1838件 [詳細] filesph_eq_11.gif 1491件 [詳細] filesph_eq_16a.gif 1411件 [詳細] filesph_eq_12.gif 1540件 [詳細] filesph_eq_23.gif 1531件 [詳細] filesph_eq_05.gif 1708件 [詳細] filesph_eq_07.gif 1600件 [詳細] filesph_eq_19.gif 1634件 [詳細] filesph_eq_04.gif 1227件 [詳細] filesph_eq_08.gif 1580件 [詳細] filesph_eq_18.gif 1530件 [詳細] filesph_eq_14.gif 1779件 [詳細] filesph_eq_15.gif 1546件 [詳細] filesph_eq_03.gif 1236件 [詳細] filesph_eq_02.gif 1580件 [詳細] filesph_eq_06.gif 1612件 [詳細] filesph_eq_10.gif 1577件 [詳細] filesph_eq_20.gif 1556件 [詳細] filesph_eq_21.gif 1611件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2012-12-12 (水) 17:56:19 (2839d)