---- #contents ---- *MPS法 [#g6ff33bf] 圧縮性流体を取り扱うために提案されたSPH法に対して, 粒子法で非圧縮性流体を取り扱えるように提案されたのが Moving Particle Semi-implicit (MPS)法¬e{S. Koshizuka, H. Tamako and Y. Oka, A particle method for incompressible viscous flow with fluid fragmentation, Comput. Fluid Dynamics J., 29(4), 1996.};である. 粒子法で非圧縮性流体を取り扱えるように提案されたのが Moving Particle Semi-implicit (MPS)法¬e{Koshizuka1996:S. Koshizuka, H. Tamako and Y. Oka, A particle method for incompressible viscous flow with fluid fragmentation, Comput. Fluid Dynamics J., 29(4), 1996.};である. MPS法では粒子数密度を導入し,この密度が一定となるように陰解法で圧力のポアソン方程式を解く. これにより,グローバルに非圧縮性を強制する. 一部分に陰解法を使っているので,Semi-implicitである. **非圧縮性流れの支配方程式 [#ode99163] 非圧縮流れの支配方程式であるナビエ・ストークス方程式は以下である. #ref(mps.eq1.gif,nolink,70%) ここで,&ref(mps.eq2.gif,nolink,70%);は流体の速度,&ref(mps.eq3.gif,nolink,70%);は液体の密度,&ref(mps.eq4.gif,nolink,70%);は粘度,&ref(mps.eq5.gif,nolink,70%);は外力を表す. **粒子間相互作用モデル [#vbee22b2] MPS法では勾配,発散,ラプラシアン等の微分演算子を,粒子間相互作用モデルに基づき離散化する. まず重み関数&ref(mps.eq6.gif,nolink,70%);を以下のように定義する. #ref(mps.eq7.gif,nolink,70%) ここで,&ref(mps.eq8.gif,nolink,70%);は粒子間距離,&ref(mps.eq9.gif,nolink,70%);は有効半径を表す. この式は粒子間距離が&ref(mps.eq9.gif,nolink,70%);より小さい粒子同士でのみ相互作用を及ぼすことを表している. 重み関数&ref(mps.eq6.gif,nolink,70%);を用いて,非圧縮性を強制するために次式で粒子数密度を定義する. #ref(mps.eq10.gif,nolink,70%) ここで,&ref(mps.eq11.gif,nolink,70%);は粒子&ref(mps.eq12.gif,nolink,70%);とその近傍粒子&ref(mps.eq13.gif,nolink,70%);の位置ベクトルを表している. 粒子数密度の初期値&ref(mps.eq14.gif,nolink,70%);を特に初期粒子数密度と呼ぶ. 非圧縮性流れにおいて液体の密度は一定である. MPS法では各粒子の粒子数密度&ref(mps.eq15.gif,nolink,70%);が&ref(mps.eq14.gif,nolink,70%);となるようにすることで非圧縮条件を強制する. 粒子間相互作用モデルによる勾配,発散,ラプラシアンの離散化は以下. -勾配 勾配は粒子&ref(mps.eq12.gif,nolink,70%);の位置に対して勾配ベクトルを与える演算子である. #ref(mps.eq16.gif,nolink,70%) ここで,&ref(mps.eq17.gif,nolink,70%);は次元数&ref(mps.eq14.gif,nolink,70%);は初期粒子数密度である. -発散 発散はベクトルに作用しスカラーを得る演算子である. #ref(mps.eq18.gif,nolink,70%) -ラプラシアン ラプラシアンは粒子&ref(mps.eq12.gif,nolink,70%);の変数量(ここでは速度,圧力等の物理量)を, 近傍の粒子&ref(mps.eq13.gif,nolink,70%);に重み関数の分布で分配する演算子である. #ref(mps.eq19.gif,nolink,70%) ここで,&ref(mps.eq20.gif,nolink,70%);は以下. #ref(mps.eq21.gif,nolink,70%) **MPS法による粒子位置更新 [#nb1c0cc7] ナビエ・ストークス方程式(質量保存式と運動量保存式)を離散化する. #ref(mps.eq22.gif,nolink,70%) 上付は用いる物理量のステップ数を示す. 圧力勾配項は陰的に計算するため,&ref(mps.eq23.gif,nolink,70%);の値を用いる. 粘性項,重力項については陽的に計算する. まず,粘性項と重力項を計算し,粒子の仮の速度&ref(mps.eq24.gif,nolink,70%);,仮の位置&ref(mps.eq25.gif,nolink,70%);を得る. #ref(mps.eq26.gif,nolink,70%) #ref(mps.eq27.gif,nolink,70%) 粘性項&ref(mps.eq28.gif,nolink,70%);はラプラシアンであるので以下のように離散化する. #ref(mps.eq29.gif,nolink,70%) 粒子位置&ref(mps.eq30.gif,nolink,70%);での粒子数密度を&ref(mps.eq31.gif,nolink,70%);とし,以下の圧力のポアソン方程式を解くことで&ref(mps.eq32.gif,nolink,70%);を得る. #ref(mps.eq33.gif,nolink,70%) 粘性項と同様に右辺を離散化すると, &ref(mps.eq34.gif,nolink,70%);に関する線形システムを得る. #ref(mps.eq35.gif,nolink,70%) Aは対称疎行列であり,不完全コレスキー分解付共役勾配法(Incomplete Cholesky Decomposition Conjugate Gradient method)などで解くことができる.求めた&ref(mps.eq32.gif,nolink,70%);から,&ref(mps.eq36.gif,nolink,70%);,&ref(mps.eq37.gif,nolink,70%);を求める. #ref(mps.eq38.gif,nolink,70%) #ref(mps.eq39.gif,nolink,70%) **計算アルゴリズム [#f2a5a1de] MPS法では以下の手順でシミュレーションを進める. +流れの初期条件を設定 +仮の速度と位置を陽的に計算 +粒子数密度を計算 +圧力のポアソン方程式を解く +圧力勾配項で速度と位置を修正 +2に戻る