2009年5月30日土曜日

粒子法のプログラム第5回(力の計算)


今回は、粒子法のプログラムでどのように力を計算しているかを説明します。

力の計算は次のコードになります。


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がカーネルの係数です。

力の計算も、密度と圧力の計算と同様に、近傍粒子探索を工夫したり対称性を利用することで、高速化できます。

次回は、どのように粒子の位置を更新するかという説明になります。

粒子法のプログラム
粒子法のプログラム第1回(概要)
粒子法のプログラム第2回(プログラムの大枠)
粒子法のプログラム第3回(データ構造)
粒子法のプログラム第4回(密度と圧力の計算)
粒子法のプログラム第5回(力の計算)
粒子法のプログラム第6回(境界条件と粒子位置の更新)
粒子法のプログラム最終回(粒子の出力)

その他の解説エントリ
SPHによる巻き波のシミュレーション1
SPHによる巻き波のシミュレーション2
SPHによる巻き波のシミュレーション3
このあとやりたいこと

固液連成シミュレーションに関するエントリ
粒子法による固液連成シミュレーション

流体シミュレーションに関するエントリ
【粒子法】粒子を流体としてレンダリング
3次元の粒子法シミュレーション
粒子法のシーンを2倍のサイズにしてみたが…
粒子法のシーンを2倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation

剛体シミュレーションに関するエントリ
粒子ベース剛体シミュレーション(プレビュー)
粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies

動画
シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。

動画の一覧


--

4 件のコメント:

匿名 さんのコメント...

力は粒子間距離に応じて線形なのでしょうか

kamonama さんのコメント...

力は、粒子間距離を引数とするカーネル関数を使っています。

下記ページの
http://www.fluid.mech.tohoku.ac.jp/public_html/sph/sph.html
図にあるような釣り鐘型のものです。

匿名 さんのコメント...

圧力、粘性、それぞれの項を数値として表示させたいのですが、うまくいきません。

Functions for outputとmain functionをいじればよいと思うのですが。。。
宜しければ教えていただきたいです。

星 さんのコメント...

いつも参考にさせていただいております。

fcurr *= c * p->rho * pj->rho;

の部分ですが、SPHの文献や参考書では、密度で割り算をしているかと存じます。
上記のような実装になるのはなぜでしょうか。
ご教授いただけますと幸いです。