はじめに,発生時間での仮想粒子の位置を初期条件として式(3)を解く。 説明のため,b を速度場 (v),s を時間 (t) とする。 ある粒子l の位置を (xl (t), yl(t), zl (t)) とすると,速度場と初期条件(xl (0), yl (0), zl (0)) が決定すれば,解析的または数値的に式(3)を解くことができる。 数値計算によって求められた速度場は空間を埋め尽くした計算セルの頂点で与えられているものとする(与えられる位置は計算手法によって異なるが頂点とすることにより一般性は失われない)。 速度場は離散点上で与えられているために粒子を追跡する(式(3)を解く)ためには,粒子が属する格子セルを探すための方法と局所的な補間によって近似されたベクトル場を求める方法が必要となる。 これらの方法は計算手法によっても異なり,それぞれの計算手法に対し効率的な方法が提案されているので,詳細については触れない。 例えば,一般曲線座標系,構造格子を用いた計算を対象とした場合の方法は文献1を参照してほしい。
粒子密度はオイラー的に計算する。 粒子l の離散時間 tn における位置を (xl , yl , zl ) とする。 粒子l は仮想的なものであるが粒子密度を考えるときは仮の質量を与える。 この粒子の質量を Ml とすると,図3に示すようにその粒子が存在する格子セルでのセル頂点への寄与を考えることにより粒子密度を求めることができる。 粒子密度を計算する方法としては,NGP法,CIC法,SPM法等がある。 また,粒子密度は,流れ場の計算を行った計算格子,または,補助的な直交格子を用いて計算される(粒子密度が算出される格子を補助格子と呼ぶこととする)。 これらの具体的な手続きについては Appendix を参照されたい。
なお,この論文の範疇では発生点(初期条件)は前もって決められているものとする。
図3 粒子密度の補助格子での算出