2009年6月20日土曜日

粒子法のプログラム第6回(境界条件と粒子位置の更新)


今回は、粒子の境界条件と粒子位置の更新について説明します。

境界条件と粒子位置の更新のコードは以下のようになります。

void calculate_position( Particles* ps )
{
vec accel(3), g(3), norm(3);
double speed, diff, adj;

g(0) = 0.0; g(1) = -9.8; g(2) = 0.0;
FOR_EACH_PARTICLES( ps, p )
{
accel = p->f * SPH_PMASS;

speed = norm2( accel );
if ( speed > SPH_LIMIT*SPH_LIMIT ) {
accel *= SPH_LIMIT / sqrt(speed);
}

// Z-axis walls
diff = 2.0 * RADIUS - ( p->r(2) - MIN(2) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 0.0; norm(1) = 0.0; norm(2) = 1.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}
diff = 2.0 * RADIUS - ( MAX(2) - p->r(2) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 0.0; norm(1) = 0.0; norm(2) = -1.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}

// X-axis walls
diff = 2.0 * RADIUS - ( p->r(0) - MIN(0) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 1.0; norm(1) = 0.0; norm(2) = 0.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}
diff = 2.0 * RADIUS - ( MAX(0) - p->r(0) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = -1.0; norm(1) = 0.0; norm(2) = 0.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}

// Y-axis walls
diff = 2.0 * RADIUS - ( p->r(1) - MIN(1) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 0.0; norm(1) = 1.0; norm(2) = 0.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}
diff = 2.0 * RADIUS - ( MAX(1) - p->r(1) ) * SPH_SIMSCALE;
if ( diff > EPSILON )
{
norm(0) = 0.0; norm(1) = -1.0; norm(2) = 0.0;
adj = SPH_EXTSTIFF * diff - SPH_EXTDAMP * inner_prod( norm, p->v );
accel += adj * norm;
}

accel += g;
p->v += accel * DT;
p->r += p->v * DT / SPH_SIMSCALE;

// stick on x-y plane
p->r(2) = 0.0;

}
}
それぞれの粒子について、粒子にかかる力から加速度を計算したあと、各壁面について境界条件を適用し、重力加速度を加えます。そうやって求まった加速度から、粒子の速度を求め、粒子の位置を更新します。

境界条件には、ペナルティ法を用いています。ペナルティ法とは、粒子が壁面へめりこんだ場合に、そのめり込みの程度によって壁面との衝突による反発力を与える方法です。

このようにして、境界条件を処理し、粒子の位置を更新しています。

次回は計算結果の出力です。最終回です。

粒子法のプログラム
粒子法のプログラム第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

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

動画の一覧


--

1 件のコメント:

hanahanasama さんのコメント...

質問なのですが、9行目のaccelはもし加速度なら、運動方程式は 「ma = F」なので
「accel = p->f * SPH_PMASS」でなく
「accel = p->f / SPH_PMASS」のような気がするのですがどうでしょうか
また次の速度を求めるときも
「speed = norm2( accel )」でなく純粋に
「speed = accel * DT」にしてもいいでしょうか。よろしくお願いします