✨

Position Based Dynamics で作るソフトボディ&流体シミュレーション

に公開
Rust
TypeScript
WebAssembly
シミュレーション
並列処理
tech

Position Based Dynamics という手法を用いて、ブラウザ上で動作する物理シミュレーションを実装しました。ぶよぶよしたソフトボディが水に浮かびます。

マウスでソフトボディをドラッグしたり、ウィンドウのリサイズでシミュレーション領域の大きさを変えられます。いろいろといじって遊んでみてください。

Demo :
https://floaty-fluid.netlify.app/

Repository :
https://github.com/matsuoka-601/Floaty

シミュレーションのソルバーは大半を TypeScript で実装しており、ボトルネックとなる流体計算の部分だけ Rust で書いています。

本記事では、このシミュレーションに関する以下の話題について説明します。

  • シミュレーションの概要
  • Position Based Dynamics について
  • ソフトボディと流体のカップリング
  • シミュレーションの高速化

シミュレーションの概要

このシミュレーションで実装したのは、Position Based Dynamics (PBD) という手法です。かねてから PBD を使ったシミュレーションを作ってみたいと思っていたのですが、以下の動画を見たことがきっかけで実装に踏み切ることにしました。

https://www.youtube.com/watch?v=fuvfpveaacY

この動画の 13 秒から始まるシーンでは、布・剛体・ロープが水に浮かぶ様子が、リアルタイムで見事にシミュレーションされています。このような素晴らしいシミュレーションを自分でも再現してみたいと思ったことが、今回の作品を作ったきっかけです。

今回の作品の最大のポイントは、ソフトボディと流体の相互作用をシミュレーションすることにより、浮力を再現したことです。流体やソフトボディ単体のシミュレーションは既にたくさん存在しますが、その 2 つを組み合わせたシミュレーションで、かつブラウザ上で動作するものは今回の作品が初めてだと思います（自分の知る限りでは）。

PBD のソルバーの大半は TypeScript で書きましたが、ボトルネックとなる流体計算の部分には Rust を使っており、rayon という並列処理のライブラリを用いたマルチスレッド並列化をしています。並列化の恩恵はかなり大きく、手元のラップトップでは、シングルスレッドの場合と比べて 3 倍程度の高速化を達成しています。

Position Based Dynamics について

ここでは、Position Based Dynamics が何者なのかについて非常にざっくりと説明します。本当にざっくりとした説明なので、もっと詳しく知りたい方は 参考文献 の部分を参照してください。

Position Based Dynamics (PBD) とは、以下のような特徴を備えた物理シミュレーションの手法です。

  • タイムステップを多少大きくとっても爆発したりしないため、リアルタイム性が求められる場面に強い
  • 流体・剛体・気体・ソフトボディ・布など、様々な物質を統一的にシミュレーションできる
  • 実装が比較的容易

「多少の物理的精度は失われてもいいので、視覚的にそれっぽく見えるシミュレーションを高速に行いたい」ようなユースケース（例: ゲーム・VFX）に向いている手法として知られています。

PBD では、各粒子に対して位置に関する制約をかけ、その制約を満たすように各粒子を少しずつ動かしていくことを繰り返してシミュレーションを進めます。もう少し具体的には、各制約について、いわば「その制約を満たさない度」を表す制約関数 C(\bm{p}_1,\dots,\bm{p}_n) を定義し、その制約関数が 0 に近づくように各粒子を動かします（\bm{p}_i は i 番目の粒子の位置ベクトル）。ここで、粒子 i の変位 \Delta\bm{p}_i は以下のように計算されます（導出は 参考文献 の元論文をご覧ください）。

\begin{equation} \Delta\bm{p}_i = \frac{C(\bm{p}_1,\dots,\bm{p}_n)} {\sum_j|\nabla_{\bm{p}_j}C(\bm{p}_1,\dots,\bm{p}_n)|^2} (-\nabla_{\bm{p}_i}C(\bm{p}_1,\dots,\bm{p}_n)) \end{equation}

ややこしく見えるかもしれませんが、「各粒子を制約関数が小さくなる方向（=勾配の逆方向）へ動かす」ことと、「制約関数の大きさが大きい（=制約を満たさない度が大きい）ほど粒子を大きく動かす」ことが分かれば、この記事を理解するうえでは OK です。

このように、各粒子の位置を直接修正することが "Position Based"（位置ベース）と銘打っている理由であり、粒子にかかる力から変位を求める "Force Based" な手法に比べて高い安定性を持つとされています。

PBD における制約には様々なものがありますが、ここでは例として最も基本的な距離制約についてとりあげてみます。距離制約とは、「\bm{p}_1 と \bm{p}_2 が距離 l でなければならない」という制約です。自然長 l のバネを想像してもらえればわかりやすいかと思います。距離制約の制約関数は以下のように定義されます。

C(\bm{p}_1, \bm{p}_2) = |\bm{p}_1-\bm{p}_2|-l

この制約関数から、\bm{p}_1 と \bm{p}_2 の変位を式 (1) にしたがってそれぞれ求めると、

\begin{aligned} \Delta\bm{p}_1 &= -\frac{1}{2}(|\bm{p}_1-\bm{p}_2|-l)\bm{d}\\ \Delta\bm{p}_2 &= +\frac{1}{2}(|\bm{p}_1-\bm{p}_2|-l)\bm{d}\\ \end{aligned}

となります（\bm{d} = \frac{\bm{p}_1-\bm{p}_2}{|\bm{p}_1-\bm{p}_2|} は \bm{p}_2 から \bm{p}_1 へ向かう単位ベクトル）。要は、\bm{p}_1 と \bm{p}_2 が l より離れていれば 2 点が近づく方向へ、逆に l よりも近ければ 2 点が遠ざかる方向へと \bm{p}_1 と \bm{p}_2 を修正するわけです。また、「制約関数が 0 から離れているほど変位の大きさが大きくなる」こともわかります。

このように、各制約を満たすように粒子を少しずつ動かしていくことでシミュレーションを進めるのが、PBD の大まかな流れになります。

ソフトボディ

今回のシミュレーションでは、粒子をリング状につなげることでソフトボディをシミュレーションしています。リングで隣接する粒子どうしは、すでに説明した距離制約によってつながれています。

ただ、距離制約だけだとリング状のソフトボディはぺしゃんこにつぶれてしまうので、それを防ぐために面積制約を各ソフトボディの粒子に対して課しています。面積制約は、粒子で囲まれた領域が一定の面積を保つようにするための制約であり、この制約によって適度に「ぷよぷよ」したソフトボディがシミュレーションできます。

また、異なるソフトボディの粒子同士がぶつかったときには、その 2 粒子に対して距離制約をかけます。具体的には、2 粒子の距離が一定距離（r とする）以下になったときに、その 2 粒子に対して距離 r の距離制約を追加します。これにより、近づきすぎた粒子同士が反発して、ソフトボディどうしがボヨンと跳ね返ります。

流体

流体のシミュレーションには、Position Based Fluids (PBF) という手法を用いています。PBFは、その名前が示唆するように、流体シミュレーションを行うための位置ベースな手法です。同じく位置ベースの手法である PBD と相性がよく、PBD のフレームワークに組み込んで流体シミュレーションをすることが可能です。

PBF では、「粒子周りの密度が一定になるようにする」という制約を各粒子に課すことで、非圧縮性を持つ流体（例: 水）のシミュレーションを行います。具体的には、i 番目の粒子に対して以下の制約関数を考えます。

C_i(\bm{p}_1,\dots,\bm{p}_n) = \frac{\rho_i}{\rho_0}-1

\rho_i は i 番目の粒子周りの密度であり、\rho_0 が通常状態での流体の密度です。各粒子は、できるだけ C=0 となるように、すなわち各粒子まわりの密度が \rho_0 に近づくようにふるまうことになります。

\rho_i を推定する方法として PBF で用いられているのは、粒子 i の近傍粒子の質量を、カーネル関数という関数で重みづけして足し合わせるというものです（W がカーネル関数）。

\rho_i = \sum_{j} m_j\times W(\bm{p}_i-\bm{p}_j)

これは、SPH 法というシミュレーション手法で標準的に用いられる密度推定手法であり、以前私が書いた ブラウザ上でヌルヌル動作する流体シミュレーションを Rust + wasm-bindgen-rayon で実装するでも説明しているので、より深く知りたい方はこちらもあわせてご覧ください。

さて、これで制約関数を具体的に計算できるので、あとは式 (1) に従って各粒子の変位を求めてやれば、PBD の場合と同様に各粒子の運動をシミュレーションできます(*)。

（* とはいいつつ、このまま実装すると流体の挙動が安定しないので、シミュレーションを安定させるための追加の工夫をいくつか実装する必要があります。これらの工夫については、PBF の元論文の 4 章 (Tensile Instability) や 5 章（Vorticity Confinement and Viscosity）をご覧ください。）

ソフトボディと流体のカップリング

ここでは、今回のシミュレーションの肝である、ソフトボディと流体の相互作用の計算と浮力の再現について説明します。とはいいつつ、浮力の再現にあたって新たな理論が必要になるわけではなく、以下の 2 つの変更を加えれば PBD のフレームワーク内で自然に浮力が再現されます。

  • PBF の計算に、流体粒子だけではなくソフトボディ粒子も含める
    • \rho_i の計算でソフトボディの粒子も含めた密度を計算します。
  • 各粒子の質量をすべて 1 としている式 (1) の代わりに、各粒子の質量を考慮した式 (2) を用いる（w_i は粒子 i の逆質量 1/m_i）
\begin{equation} \Delta\bm{p}_i = w_i\times\frac{C(\bm{p}_1,\dots,\bm{p}_n)} {\sum_j w_j|\nabla_{\bm{p}_j}C(\bm{p}_1,\dots,\bm{p}_n)|^2} (-\nabla_{\bm{p}_i}C(\bm{p}_1,\dots,\bm{p}_n)) \end{equation}

このような変更を加えると、ソフトボディ粒子の質量を流体粒子の質量よりも小さくすることで、浮力が再現されます。

なお、浮力の実装にあたって Unified Particle Physics for Real-Time Applications という論文を参考にしました。ここに書いたのはこの論文をかいつまんだ内容なので、より詳しく知りたい場合は論文の 7 章を参照してください。

シミュレーションの高速化

このシミュレーションのパフォーマンスのボトルネックは、PBF の流体計算です。幸いにも、PBF は並列化しやすいアルゴリズムなので、マルチスレッドによる並列化を行うことにしました。

今回のシミュレーションのソルバーの大半は TypeScript で実装しましたが、PBF の部分は Rust で実装し、並列化に rayon というクレートを用いています。今回のシミュレーションでは並列化が効果てきめんで、手元のラップトップでは 3 倍近い高速化が達成できています。


1スレッドの場合のパフォーマンス


7スレッドの場合のパフォーマンス

今回のデモは比較的低いスペックの端末でもスムーズに動くことを確認していますが、並列化なしではこのレベルのパフォーマンスを達成することは難しかっただろうと考えています。

まとめ

以下が本記事のまとめです。

  • Position Based Dynamics (PBD) と Position Based Fluids (PBF) に基づくソフトボディ&流体シミュレーションを実装
  • 以下の 2 つの変更を PBF のソルバーに加えることで浮力を再現
    • PBF の計算にソフトボディの粒子も含める
    • 重み付きの更新式を用いる
  • rayon クレート により PBF ソルバーを並列化

今回は、ソフトボディと流体という目新しい組み合わせのシミュレーションを実装したわけですが、比較的スムーズに浮力の再現まで漕ぎつけられたのは PBD の持つシンプルさや柔軟性ゆえだと感じています。学習するための資料も豊富にあるので、インタラクティブな物理シミュレーションをはじめて作ってみたい方に PBD はおすすめです。

今後は、シミュレーションをさらに高速化しようと目論んでいます。そのために有効だと考えているのが、WebAssembly SIMD の活用です。実は、WebAssembly SIMD を用いて高速化した流体シミュレーションはすでにいくつか存在しており、saharan さんの Drops や Grant Kot さんの Liquid Layers がその実例です。これらのシミュレーションでは、いずれも WebAssembly SIMD による 2~3 倍の高速化が報告されているため、今回実装したシミュレーションも SIMD による高速化の余地が大いにあると考えています。最終的には、ローエンドのスマホでも 60FPS で動作させることを目標に、開発を続けていきたいと考えています。

参考文献

Discussion