今回は、粒子法のプログラムでどのように力を計算しているかを説明します。
力の計算は次のコードになります。
void calculate_force( Particles* ps )
{
double pterm, vterm, r, c;
vec dr(3), force(3), fcurr(3);
FOR_EACH_PARTICLES( ps, p )
{
force(0) = 0.0; force(1) = 0.0; force(2) = 0.0;
FOR_EACH_PARTICLES( ps, pj )
{
if ( p == pj ) continue;
dr = ( p->r - pj->r ) * SPH_SIMSCALE;
r = norm( dr );
if ( H > r )
{
c = H - r;
pterm = -0.5 * c * SpikyKern * ( p->p + pj->p ) / r;
vterm = LapKern * SPH_VISC;
fcurr = pterm * dr + vterm * ( pj->v - p->v );
fcurr *= c * p->rho * pj->rho;
force += fcurr;
}
}
p->f = force;
}
}それぞれの粒子について、ナビエ=ストークス方程式の圧力項、粘性項にあたる力を計算しています。影響半径内にある自分以外の粒子に対して、距離に応じて重みを足し合わせています。SpikyKernとLapKernがカーネルの係数です。
力の計算も、密度と圧力の計算と同様に、近傍粒子探索を工夫したり対称性を利用することで、高速化できます。
次回は、どのように粒子の位置を更新するかという説明になります。
粒子法のプログラム
粒子法のプログラム第１回（概要）粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
その他の解説エントリ
SPHによる巻き波のシミュレーション１SPHによる巻き波のシミュレーション２
SPHによる巻き波のシミュレーション３
このあとやりたいこと
固液連成シミュレーションに関するエントリ
粒子法による固液連成シミュレーション流体シミュレーションに関するエントリ
【粒子法】粒子を流体としてレンダリング３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation
剛体シミュレーションに関するエントリ
粒子ベース剛体シミュレーション（プレビュー）粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies
動画
シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。--
粒子法(SPH)のプログラムを解説したシリーズです。ソースコードも公開しています。