今回は、粒子の境界条件と粒子位置の更新について説明します。
境界条件と粒子位置の更新のコードは以下のようになります。
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;
}
}それぞれの粒子について、粒子にかかる力から加速度を計算したあと、各壁面について境界条件を適用し、重力加速度を加えます。そうやって求まった加速度から、粒子の速度を求め、粒子の位置を更新します。境界条件には、ペナルティ法を用いています。ペナルティ法とは、粒子が壁面へめりこんだ場合に、そのめり込みの程度によって壁面との衝突による反発力を与える方法です。
このようにして、境界条件を処理し、粒子の位置を更新しています。
次回は計算結果の出力です。最終回です。
粒子法のプログラム
粒子法のプログラム第１回（概要）粒子法のプログラム第２回（プログラムの大枠）
粒子法のプログラム第３回（データ構造）
粒子法のプログラム第４回（密度と圧力の計算）
粒子法のプログラム第５回（力の計算）
粒子法のプログラム第６回（境界条件と粒子位置の更新）
粒子法のプログラム最終回（粒子の出力）
その他の解説エントリ
SPHによる巻き波のシミュレーション１SPHによる巻き波のシミュレーション２
SPHによる巻き波のシミュレーション３
このあとやりたいこと
固液連成シミュレーションに関するエントリ
粒子法による固液連成シミュレーション流体シミュレーションに関するエントリ
【粒子法】粒子を流体としてレンダリング３次元の粒子法シミュレーション
粒子法のシーンを２倍のサイズにしてみたが…
粒子法のシーンを２倍のサイズにしてみた
Haskell、OCamlでSPH法
カメラ位置を変えて流体をレンダリング
Bunny-shaped fluid simulation
剛体シミュレーションに関するエントリ
粒子ベース剛体シミュレーション（プレビュー）粒子ベース多体衝突シミュレーション
引き続き、粒子ベース剛体シミュレーション
Falling Rigid Bunnies
動画
シミュレーションの結果をレンダリングして作った動画です。流体シミュレーションや剛体シミュレーションの動画を見ることができます。--
