2次元流体シミュレーション: 自由表面と剛体との相互作用

自由表面、剛体との相互作用のある、2次元の格子法流体シミュレーションをBevyエンジンを用いてGPUで開発しました。

開発したシミュレーションを利用したデモがこちらにありますので、ぜひご覧ください。ソースコードも公開しています: https://github.com/narasan49/bevy_eulerian_fluid

動機

ゲームで用いられる流体表現には様々ありますが、ユーザーが触れて遊ぶことのできるものはあまりないなと思い、作ってみたくなりました。流体に触れた結果さらにほかの要素に影響を与えてほしいので剛体との相互作用を、流体として身近なものはやはり水なので水面(自由表面)を含むシミュレーションを開発しました。

実装の概要

流体シミュレーションは次のNavier-Stokes方程式と非圧縮の条件を解くことで実現できます。

$$ \begin{align*} \frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u}\cdot\nabla)\boldsymbol{u} &= -\frac{1}{\rho} \nabla p + \boldsymbol{F} \\ \nabla \cdot \boldsymbol{u} &= 0 \end{align*} $$

GPU GemsやFluid Simulation for Computer Graphicsが詳しいですが、Navier-Stokes方程式の各項を分離して順に解くことで、解を近似することができます。

$$ \begin{align*} &\frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u}\cdot\nabla)\boldsymbol{u} = 0 & (移流) \\ &\frac{\partial \boldsymbol{u}}{\partial t} = \boldsymbol{F} & (外力) \\ &\frac{\partial \boldsymbol{u}}{\partial t} = -\frac{1}{\rho} \nabla p & (圧力勾配) \end{align*} $$

移流項は速度場によって物理量が流される効果を表しますが、この場合は速度自身が速度場によって流されることを表しています。移流項を解くにはセミラグランジアンと呼ばれる手法を用い、注目しているセルを速度場によって流される粒子とみなし、前の時刻でその粒子がどこにいたかを推定します。前の時刻の位置における速度が、注目しているセルの次の時刻における速度となります。

外力項では、ユーザー入力などによって与えられる局所的な力や重力を扱っています。外力項は移流後の速度$u'$を用いて前進オイラー差分で更新できます。 $$ \boldsymbol{u''} = \boldsymbol{u'} + \Delta t \boldsymbol{F} $$

最後に圧力勾配力の式を解きます。非圧縮の条件($\nabla \cdot \boldsymbol{u} = 0$)と同時に解くことで、流体の非圧縮性を維持できます。後程詳しく説明します。

剛体との相互作用については、既存の剛体シミュレーション利用し、剛体の情報をGPUへの転送とReadbackにより実現しています。これらに加え、剛体シミュレーションとの相互作用や流体の流入流出などを加えると、シミュレーションのGPUパイプライン全体の構成は以下のようになります。

剛体のレベルセット更新
速度場の移流
外力の適用
圧力勾配力と非圧縮性による速度の更新
速度場の外挿
レベルセットの移流
流体の流入流出の計算
レベルセットの再初期化
流体から剛体へ与える力の評価

以下では、工夫した点などをかいつまんで紹介します。

レベルセット

自由表面における流体と空の領域の区別、剛体との相互作用における流体と剛体の区別には、レベルセットと呼ばれるSDFを用います。

自由表面おいてレベルセット$\phi^f$は$\phi^f < 0$を流体の内側、$\phi^f\ge0$を流体の外側と表現します。
剛体においても同様で$\phi^s$は$\phi^s < 0$を剛体の内側、$\phi^s \ge 0$を剛体の外側と表現します。

剛体のレベルセットはシーン上の剛体の形状、TransformなどをGPUに転送して構築しています。

離散化

今回は格子法という、計算領域を等間隔の格子に分割してシミュレーションを行う方法をとっており、格子の構造には、Staggered Gridと呼ばれる構造を採用しています。

Staggered Gridは速度がグリッド幅の半分だけずれた位置にあり、それ以外の量はグリッドの中央に位置している格子構造です。x方向の速度はx方向に半分ずれた点$u_{i-1/2, j}$に、y方向の速度はy方向に半グリッドずれた点$u_{i, j-1/2}$に位置します。

グリッド中央の発散場(中間的な速度場は必ずしも $\nabla \cdot \boldsymbol{u} = 0$ ではありません)を計算する際、すべての物理量がグリッド中央にある格子(Collocatedといいます)の中心差分は$(u_{i+1, j} - u_{i-1, j}) / 2\Delta x$となる一方、Staggered Gridの場合、$(u_{i+1/2, j} - u_{i-1/2, j}) / \Delta x$となります。 $u_{i+1/2, j}$と$u_{i-1/2, j}$はデータ上は隣り合う要素のため、Collocatedの離散化と同じ精度で解像度が2倍になるメリットがあります。

有限体積法

プロジェクション

ここでは圧力勾配力による速度の更新についてみていきます。 $$ \begin{align*} \frac{\partial \boldsymbol{u}}{\partial t}= -\frac{1}{\rho} \nabla p \\ \nabla \cdot \boldsymbol{u} = 0 \end{align*} $$

圧力勾配は以下のように離散化されます。 $$ u_{i-1/2, j}^{n+1} = u_{i-1/2, j}'' - \frac{\Delta t}{\rho\Delta x}(p_{i, j} - p_{i - 1, j}) $$

これを離散化した非圧縮の式に代入することで圧力$p$に関する式として解くことができます。非圧縮の式を有限差分法によって差分化すると、

$$ \frac{u_{i+1/2, j}^{n + 1} - u_{i-1/2, j}^{n + 1}}{\Delta x} + \frac{v_{i, j+1/2}^{n + 1} - v_{i, j-1/2}^{n + 1}}{\Delta x} = 0 $$

となり、圧力勾配力の式を代入すると、

$$ u_{i+1/2, j}'' - u_{i-1/2, j}'' + u_{i, j+1/2}'' - u_{i, j-1/2}'' + \frac{\Delta t}{\rho\Delta x}(4p_{i, j} - p_{i - 1, j} - p_{i + 1, j} - p_{i, j - 1} - p_{i, j + 1}) = 0 $$

となります。この式により圧力場を更新し、圧力勾配力の式に再度代入することで、最終的な速度場を得られます。圧力勾配力と非圧縮の条件を同時に解くことで最終的な速度場が非圧縮性を維持することができます。

自由表面、剛体とのカップリングがある場合、セルの状態によって修正が必要になります。 $(i, j)$の隣接セル$(i-i, j)$が空気の場合は$p_{i-1,j}=0$とします。隣接セル$(i-i, j)$が剛体だった場合、その境界の速度$u_{i-1/2, j}^{n+1}$は剛体の速度$u_{i-1/2, j}^{solid}$となります。その分、上の式においては$p_{i, j}$の係数は1小さくなり、$p_{i-1,j}$の項も0になります。

この圧力勾配と非圧縮の式から速度を更新する過程はプロジェクションと呼ばれます。

有限体積法

有限差分法では、各セルに単一の状態(剛体、空気、流体)が割り当てられていましたが、境界付近では通常、一つのセルで複数の状態が混在しています。このような状態を正しく扱わないと、剛体の周辺の流れにアーティファクトが現れたり、水面が正しく追跡できず、流体の体積の損失を引き起こしたりします。

このような問題を改善する方法に有限体積法があります。有限体積法は非圧縮条件をセルの境界に沿って積分することにより、より精度の高い評価を可能にしています。

$$ \int_{\partial C}{\nabla \cdot \boldsymbol{u}} = 0 $$

$\partial C$ はセルの境界に沿った経路での積分です。 $\nabla \cdot \boldsymbol{u}= 0$なので、積分しても0になります。2次元の場合、セル$(i, j)$について以下のように離散化されます。

$$ \begin{align*} F_{i+1/2, j} u_{i+1/2, j}^{n + 1} &- F_{i-1/2, j} u_{i-1/2, j}^{n + 1} \\ + F_{i, j+1/2} v_{i, j+1/2}^{n + 1} &- F_{i, j-1/2} v_{i, j-1/2}^{n + 1} \\ + (1 - F_{i+1/2, j}) u_{i+1/2, j}^{solid} &- (1 - F_{i-1/2, j}) u_{i-1/2, j}^{solid} \\ + (1 - F_{i, j+1/2}) v_{i, j+1/2}^{solid} &- (1 - F_{i, j-1/2} )v_{i, j-1/2}^{solid} = 0\\ \end{align*} $$

ここで、$F_{i-1/2, j}$はセルの左側の境界が占める非剛体の割合を表しています。$F_{i-1/2, j}=1$なら$(i-1/2, j-1/2)$から$(i-1/2, j+1/2)$は完全に流体または空気で、$F_{i-1/2, j}=0$なら、完全に剛体です。$F_{i-1/2, j}$の値は剛体のレベルセットから計算可能です。これによりサブグリッドスケールで剛体と流体の境界を扱えるようになりました。

自由表面の圧力についても、単に$p_{i-1,j}=0$とするのではなく、$\phi^f=0$の位置が$p=0$となるように$p_{i-1,j}$の値を外挿しています。

$$ \begin{equation*} p_{i-1, j} = \begin{cases} p_{i-1, j} & (\phi_{i-1,j}^f < 0) \\ \frac{\phi_{i-1,j}^f}{\phi_{i,j}^f}p_{i, j} & (\phi_{i-1,j}^f \ge 0) \end{cases} \end{equation*} $$

プロジェクションの解法

プロジェクションをGPUで解くには反復解法が一般的で、Jacobi法、Gauss-Seidel法といった解法が一般的で、解が十分収束するためには解像度にもよりますが、100回単位の反復が必要になり、シミュレーションのボトルネックとなります。煙のシミュレーションなどではある程度のところで打ち切ることでパフォーマンスを優先する選択を取ることができます。

しかし、水面のある水のシミュレーションではプロジェクションの収束が十分でないと、質量が保存せず全体の水位が上下に振動したり、質量が失われたりと好ましくない状況になってしまいます。

この問題を解決できる方法に、Multigridという方法があります。プロジェクションは変数を平滑化するような操作で、一回のイテレーションで隣り合うセルとの平滑化を行うので高周波の誤差はすぐに収束する一方で、低周波の誤差は収束するのに多くのイテレーションを要します。そこでグリッドの解像度を段階的に粗くしながらGauss-Seidelなどのイテレーションを回すことで低周波の誤差も効率よく収束させることができます。

Multigridにはいくつかの型がありますが、今回はV-Cycleという実装を採用しており、pre smoothと呼ばれる工程で数回のGauss-Seidelを実施します。その後残差を計算し、Ristrictionという工程で解像度を落とします(私の実装では、x, yともに1/2にしています)。これを最も粗い解像度(512x256のシミュレーションでは16x8まで落としています)まで繰り返します。最も粗い解像度では、pre smoothよりも多く、収束するまでイテレーションを回します。その後Prolongationと呼ばれる工程で、粗い解像度で求めた結果を一段細かいグリッドで加算します。その後、Post smoothの工程で、数回のGauss-Seidelのイテレーションを回して均します。

これがMultigridの一連の流れです。低解像度にすることで一回のイテレーションに要する時間が短くなり、さらに収束までの全体のイテレーション数も削減することができ、高速化につながります。 Multigridについてはこちらの動画シリーズがとても丁寧に解説してくれています。

以下にそれぞれの手法における水面の挙動についての例を掲載しています。Jacobi法は大幅に水面が上下しています。反復100回のGauss-Seidel法は水面の振動は小さいです。Multigridはほとんど振動していません。実行環境はSnapdragon X-Eliteの内蔵GPUです。

反復100回のJacobi法

反復100回のGauss-Seidel法

Multigrid法

Multigridは精度が高いだけでなく、Gauss-Seidelよりも速く収束します。512x256解像度における各手法の所要時間をPIXキャプチャで取得すると、以下のようになり、圧倒的にMultigridが高速です。

手法	Projectionの所要時間 (ms)	反復回数
Jacobi	6.0	100
Gauss-Seidel	9.0	100
Multigrid	0.83	pre/post smooth 2回, 最粗グリッドでの反復: 3回, 最粗解像度: 16x8

レベルセットの移流

レベルセットによって、流体と空領域の境界、剛体と流体の境界を表現できますが、レベルセットを移流することによって、水面の変位を表現できます。しかし、移流操作に含まれる補間などにより、レベルセットが距離場ではなくなってしまうため、精度よくシミュレーションを継続させるためには、再初期化という操作を定期的に行う必要があります。

再初期化とは、現在の境界付近のレベルセットを正しいものとして、SDFを再構築する操作で、SDFの再構築の手法としてJump Flooding AlgorithmとFast Iterative Methodを試しました。

Jump Flooding Algorithmは与えられたシード点に関するボロノイ図を構築する手法で並列性能が高くGPUでの実装に向いています。

Fast Iterative Methodは境界付近から離れた方向に向かって、反復的な操作でSDFを構築する手法で、Jump Flooding Algorithmよりも精度が良いです。

水面の移流では、境界付近のみで精度のよいレベルセットがあれば十分なので、Fast Iterative Methodを採用しました。

それぞれの手法を以下の動画で比較してみます。Jump Flooding Algorithmの場合、水位が時間とともに下がっていっていますが、Fast Iterative Methodの場合、水位はほぼ一定でした。さらに、Fast Iterative Methodは境界付近のSDFを再構築すればよいので、パフォーマンス面でも優位なことが多いです。

Jump Flooding Algorithm

Fast Iterative Method

まとめ

以上、剛体との相互作用、自由表面のある2次元の流体シミュレーションを開発した際の工夫を簡単に紹介しました。 MultigridやFast Iterative Methodにより、内蔵GPUのlaptopでもリアルタイムに高品質な水面のあるシミュレーションを実現できました。これからも剛体の情報のCPU-GPU間のやり取りの最適化や3Dシミュレーションなどにも取り組みたいと思っています。リポジトリ(https://github.com/narasan49/bevy_eulerian_fluid)も公開しているので、実装の詳細などもご覧いただけます。