位置ベース流体(Position Based Fluids)

粒子法を用いた流体シミュレーションにおいて重要な要素となるのが 非圧縮性の強制である. 非圧縮性とは流体の密度が流れによって空間的にも時間的にも変化しないことを表し, 水のような液体はもちろん, 流速が音速より十分小さいならば空気も非圧縮性流体と考えることができる (ちなみに音速に近くなると空気は圧縮する, 音速に達した飛行機が発する衝撃波いわゆるソニックブームは空気が圧縮することで引き起こされる). 非圧縮性のための様々な方法(MPS,WCSPH,PCISPH,IISPHなど)が提案されているが, それらに比べて非常に高速かつ安定した計算が可能なのが位置ベース流体法(Position Based Fluids, 以下PBF)(1) である.

PBFは位置ベース法(Position-based method)(2) の一種である. 位置ベース法では従来の力学ベースの方法と異なり,位置(物体の空間座標値)を満たすべき制約に従って直接変更する方法である 力学ベースでは支配方程式に従い加速度->速度->位置と変更していくため, それぞれで積分に伴う誤差の蓄積が発生するが位置に直接制約をかける位置ベース法はそれがない. ただし,頂点や粒子の位置関係に基づく制約をかけるためどちらかというと幾何学的な方法で, あまり物理的ではないのではないかという議論もある.

位置ベース法の基本

PBFの前に基本となる位置ベース法の考え方を説明する. まず,移動する計算点pbf.eq1.gifに対して,その座標値そのものに制約条件pbf.eq2.gifを掛けることを考える. pbf.eq3.gifの条件が常に満たされているとすると, 粒子が微少量pbf.eq4.gif移動した後もこの条件は満たされているということになるので, pbf.eq5.gifとなる. これをテイラー展開してやれば,

pbf.eq6.gif

となる.2次以降の項を無視すればpbf.eq4.gifに関する方程式が得られる.

pbf.eq7.gif

pbf.eq8.gifは制約として定義しているのでそこから計算できる. いま,pbf.eq1.gifを粒子の中心座標とすると求めるべき未知数の数は 粒子数をpbf.eq9.gifとしたら3次元空間の場合pbf.eq10.gifとなる. 一方でpbf.eq11.gifが各粒子ごとに計算される(つまりpbf.eq12.gif)ならば, 上記の方程式の数はpbf.eq9.gifなので未知数pbf.eq10.gifに対して少なすぎて解が定まらない (疑似逆行列などを使ってできないことはない). そこでpbf.eq4.gifに対して別の制約をかける.

計算点pbf.eq1.gifが運動方程式に従って移動していると考えると, ここでのpbf.eq4.gifはその運動とは別の制約pbf.eq11.gifを満たすための移動となる. そうすると元の運動(平行移動と回転)を阻害しないような移動方向, 言い換えると平行移動と回転運動量を保存するような移動方向を定める必要がある. 制約pbf.eq11.gifがこれらの運動の状態に全く作用されない(依存しない)と考えるならば, その傾きpbf.eq13.gifはこれらの運動方向に対して垂直な方向となる. つまり,pbf.eq4.gifの方向をpbf.eq13.gifとなるような制約を掛けてやればよいことになる. これはちょうど制約付きの最適化問題になるのでラグランジュ乗数pbf.eq14.gifを導入して,

pbf.eq15.gif

とする. これを上のpbf.eq4.gifに関する方程式に代入して整理すると,

pbf.eq16.gif

となる.求めたpbf.eq14.gifpbf.eq4.gifの式に代入すれば制約を満たすための移動ベクトルが分かるので, それで計算点位置を修正してやればよい. 複数の計算点が絡む場合は1回の修正では制約を満たすことは難しいので, 上記の式による制約を条件を満たすまで反復処理する.

位置ベース法の粒子法への適用

さて,話を粒子法の場合に戻そう. まず,非圧縮性条件が満たすべき制約(Constraint)を考える. 流体が圧縮しないということはその密度pbf.eq17.gifが初期密度pbf.eq18.gifから変化しないということであり, つまり,条件としてはpbf.eq19.gifである.これを変形して位置pbf.eq1.gifに関する制約pbf.eq11.gifとして定式化すると,
  pbf.eq20.gif    ・・・(1)
となる.

ここでは粒子法の一種であるSPH法を使ってある粒子pbf.eq21.gifの密度pbf.eq22.gifを求める.

pbf.eq23.gif

ここでpbf.eq24.gifは粒子pbf.eq21.gifの近傍粒子の集合,pbf.eq25.gifは粒子質量,pbf.eq26.gifは有効半径である. 近傍粒子数をpbf.eq9.gifとすると粒子pbf.eq21.gifの密度に関する制約は近傍粒子座標pbf.eq27.gifの 関数となるので,pbf.eq14.gifの計算式は以下のようになる.
  pbf.eq28.gif    ・・・(2)

pbf.eq29.gifは密度pbf.eq22.gifの式を式(1)に代入すれば計算できる. pbf.eq30.gifも同様で,SPH法における勾配値は単に重み関数pbf.eq31.gifの勾配をとればよいので,

pbf.eq32.gif

となる. この式で注意しないといけないのはpbf.eq31.gifpbf.eq33.gifpbf.eq34.gif(とpbf.eq26.gif)に関する関数なので, pbf.eq35.gifで微分したときに0にならないのはpbf.eq36.gifpbf.eq37.gifの場合だけである. そうすると上記の式はさらに書き下すことができて,
  pbf.eq38.gif    ・・・(3)

となる.pbf.eq36.gifの時はpbf.eq39.gifに関わらずpbf.eq40.gifpbf.eq35.gifを含むのでそのまま計算すればよい. pbf.eq41.gifの場合はpbf.eq42.gifとなったときだけpbf.eq40.gifを含むのでそれだけ計算する (マイナスの符号がついているのはpbf.eq43.gifなので). また,粒子質量pbf.eq25.gifがすべての粒子で一定ならばpbf.eq44.gifの外に出せて, 最終的にpbf.eq4.gifの計算ではすべて打ち消されるので計算上は省略してもよい.

式(2)と式(3)により各粒子のラグランジュ乗数pbf.eq45.gifが計算されるので, これらから各粒子の位置修正量pbf.eq46.gifを求める. このとき近傍粒子への影響と近傍粒子からの影響が非対称にならないようする(SPH法の圧力項等と同じ).
  pbf.eq47.gif    ・・・(4)

先ほども説明したように質量一定ならばこちらのpbf.eq48.gifも省略できる. pbf.eq45.gifpbf.eq49.gifの平均をとるならばpbf.eq50.gifをつけるところだが, SOR法の加速緩和係数と同じ役割を果たすようで,ない方が収束は速くなる.

張力の安定化

SPH法では近傍粒子数が少ない場合に負の圧力が発生して粒子同士が過度に集まってしまう (particle clusteringやparticle stackingと呼ばれている)現象が発生する. PBFではこれを解決するためにMonaghan (3) が提案した以下のような項を式(4)に追加している.

pbf.eq51.gif

式中のパラメータの値としてはpbf.eq52.gif,pbf.eq53.gif,pbf.eq54.gifあたりがよいようである. 式(4)のpbf.eq55.gifの部分をpbf.eq56.gifとするのだが, Monaghanの元の論文ではpbf.eq57.gifは加速度なので位置修正に直接使うのはおかしくなるので, 下記の私の実装ではpbf.eq58.gifを掛けてある.

境界粒子

前述のparticle stackingの問題は固体境界付近でも発生する. 粒子と固体境界の単純な衝突処理だけだと固体内に粒子が存在しないため, 境界付近に粒子が集まってしまう. それを解決する一つの方法が固体境界内に位置固定の粒子(境界粒子)を散布することである. 通常,境界粒子は有効半径に合わせて複数層配置するが, Akinciらの方法(4) を使えば,仮想的な体積を事前に計算しておくことで1層だけで済ますことができる. 境界に沿って粒子を配置した後,境界粒子だけを使って以下の仮想体積を計算して格納しておく.

pbf.eq59.gif

そして,PBFでも計算時に境界粒子の質量をpbf.eq60.gifとしてやればよい (詳しい式などはこちらの論文を参照).

実装例

Visual Studio 2012 + CUDA7.5の環境で作成したコードは以下.

  • 簡単な説明&注意事項
    • GUIツールキットとしてFLTKを用いている.
    • rx_fltk_glcanvas.hの55行目
      #define RX_USE_GPU
      をコメントアウトするとCPUで計算する.
    • OpenVDBはポリゴンモデルを読み込んでその内部に境界粒子を配置する場合に必要.使用する場合はプロプロセッサ定義にRX_USE_OPENVDBを追加する.
    • モデルの読み込み自体は3Dモデルファイルの入出力のrx_model.libを使っている.
    • Visual Studioから実行する場合は,作業ディレクトリを"../../bin"にしておく.
    • "s"キー or Start/Stopボタンでシミュレーションスタート/ストップ.
    • F1〜F9キー or Sceneメニューからシーン切替可能.binフォルダに"sph_scene_*.cfg"という名前のファイルがあり, これがシーンを記述するファイル.デフォルトでは4つのシーンが入っている.
    • CUDAはv7.5で動作を確認している.
    • ポリゴンモデルを固体境界とした場合に境界粒子配置のためにOpenVDBを用いている(正確には境界粒子を配置するための陰関数場生成のため).デフォルトではOpenVDB使用をOFFにしているので,ポリゴン境界に粒子を配置したい場合はプリプロセッサでRX_USE_OPENVDBを定義するか,rx_sph_solid_poly.cppの先頭でdefineすること.

(1)  M. Macklin and M. Muller, "Position based fluids" ACM Trans. Graph., 32, pp.104:1-104:12, 2013.
(2)  M. Muller, B. Heidelberger, M. Hennix and J. Ratcliff, "Position based dynamics", J. Vis. Comun. Image Represent., 18, pp.109-118, 2007.
(3)  J. Monaghan, "SPH without a Tensile Instability", Journal of Computational Physics, 159, pp.290-311, 2000.
(4)  N. Akinci, M. Ihmsen, G. Akinci, B. Solenthaler and M. Teschner, "Versatile rigid-fluid coupling for incompressible SPH", ACM Trans. Graph., ACM, 31, pp.62:1-62:8, 2012.

添付ファイル: filerx_pbf.zip 7307件 [詳細] filepbf.eq46.gif 643件 [詳細] filepbf.eq44.gif 704件 [詳細] filepbf.eq40.gif 575件 [詳細] filepbf.eq1.gif 639件 [詳細] filepbf.eq49.gif 663件 [詳細] filepbf.eq6.gif 640件 [詳細] filepbf.eq39.gif 594件 [詳細] filepbf.eq58.gif 564件 [詳細] filepbf.eq16.gif 584件 [詳細] filepbf.eq60.gif 586件 [詳細] filepbf.eq51.gif 580件 [詳細] filepbf.eq4.gif 540件 [詳細] filepbf.eq13.gif 568件 [詳細] filepbf.eq22.gif 588件 [詳細] filepbf.eq3.gif 603件 [詳細] filepbf.eq17.gif 583件 [詳細] filepbf.eq24.gif 539件 [詳細] filepbf.eq59.gif 554件 [詳細] filepbf.eq23.gif 562件 [詳細] filepbf.eq14.gif 560件 [詳細] filepbf.eq30.gif 539件 [詳細] filepbf.eq47.gif 565件 [詳細] filepbf.eq10.gif 588件 [詳細] filepbf.eq2.gif 538件 [詳細] filepbf.eq9.gif 537件 [詳細] filepbf.eq42.gif 613件 [詳細] filepbf.eq57.gif 583件 [詳細] filepbf.eq55.gif 541件 [詳細] filepbf.eq53.gif 562件 [詳細] filepbf.eq38.gif 521件 [詳細] filepbf.eq52.gif 597件 [詳細] filepbf.eq35.gif 524件 [詳細] filepbf.eq48.gif 558件 [詳細] filepbf.eq20.gif 513件 [詳細] filepbf.eq45.gif 573件 [詳細] filepbf.eq50.gif 505件 [詳細] filepbf.eq36.gif 521件 [詳細] filepbf.eq31.gif 607件 [詳細] filepbf.eq33.gif 541件 [詳細] filepbf.eq25.gif 583件 [詳細] filepbf.eq27.gif 549件 [詳細] filepbf.eq21.gif 548件 [詳細] filepbf.eq15.gif 537件 [詳細] filepbf.eq5.gif 570件 [詳細] filepbf.eq26.gif 539件 [詳細] filepbf.eq29.gif 586件 [詳細] filepbf.eq37.gif 577件 [詳細] filepbf.eq19.gif 564件 [詳細] filepbf.eq32.gif 557件 [詳細] filepbf.eq54.gif 559件 [詳細] filepbf.eq8.gif 540件 [詳細] filepbf.eq11.gif 544件 [詳細] filepbf.eq7.gif 541件 [詳細] filepbf.eq41.gif 544件 [詳細] filepbf.eq56.gif 568件 [詳細] filepbf.eq34.gif 566件 [詳細] filepbf.eq18.gif 574件 [詳細] filepbf.eq12.gif 643件 [詳細] filepbf.eq43.gif 570件 [詳細] filepbf.eq28.gif 667件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2015-11-12 (木) 16:06:10 (2214d)