はじめに
初めまして。Craft Eggの3D開発室というところで描画・TA業務を行っている江口大喜と申します。
今回のブログはCraft Eggアドカレ16日目になります。
このブログについて
このブログもこの機会に開設しましたが、会社とは関係ない個人ブログになります。
これから自分がInputしたことをOutputしていく場にできればと思います。
Unityの流体シミュレーションについて
Unityの流体シミュレーションについて、概念や理論、そして実装まで説明しようとするとかなりの分量になってしまいます。そこで本記事では、概念および格子法をしっかり説明していきたいと思います。
また、今回示す実装に関しては以下URLのリポジトリにありますUnityGraphicsProgrammingを参考にしていきます。 github.com
また、理論内容は以下のReal-Time Fluid Dynamics for Gamesを参考にします。
http://graphics.cs.cmu.edu/nsp/course/15-464/Fall09/papers/StamFluidforGames.pdf
どちらも以下の1999年発表の論文、Stable Fluidsを元にしています。
https://d2f99xq7vri1nk.cloudfront.net/legacy_app_files/pdf/ns.pdf
今回の記事の目標
今回の記事で、上記に挙げたReal-Time Fluid Dynamics for Games、UnityGraphicsProgrammingを理解できるようになれればいいなと思います。
また、自分も勉強中ですので間違いなどあったらご教示いただきたいです。
ナビエストークス方程式
最終アウトプット
まず、今回実装していく流体シミュレーション/格子法とはどのようなものになるのでしょうか。最終アウトプットの動画を見てください。
ナビエ・ストークス方程式
まず、流体シミュレーションを始めるにあたってナビエ・ストークス方程式を少し理解する必要があります。
ナビエ・ストークス方程式は以下になります。
$$ ρ(\frac{∂\overrightarrow{u}}{∂t} +(\overrightarrow{u}・\nabla)\overrightarrow{u})=-\nabla{P}+μ\nabla{^2}\overrightarrow{v}+ρ\overrightarrow{f} $$
また、実装するに当たっては両辺をρを割り、 $$ \frac{∂\overrightarrow{u}}{∂t} +(\overrightarrow{u}・\nabla)\overrightarrow{u}=-\frac{\nabla{P}}{ρ}+\frac{μ\nabla{^2}\overrightarrow{v}}{ρ}+\overrightarrow{f} $$
として、左辺の第二項を右辺に移項した $$ \frac{∂\overrightarrow{u}}{∂t} =-(\overrightarrow{u}・\nabla)\overrightarrow{u} -\frac{\nabla{P}}{ρ}+\frac{μ\nabla{^2}\overrightarrow{v}}{ρ}+\overrightarrow{f} $$
という形で使いました。以下の導出では一つ目のナビエ・ストークスを使いますが、実装の観点では最後に挙げた式の左辺にある$$\frac{∂\vec{u}}{∂t}$$を算出することが目標とします。ある粒子が動いた様が分かれば、流体を記述することができそうです。
そして単語として、最後に挙げた式の第1項を移流項、第2項を圧力項、第3項を粘性項、第4項を外力項とします。
移流項について
まず移流項についてです。
以下のように流体を微小区間に分けたボックスを考えます。それがΔtの間に移動したとします。この時、移動後の位置は古典力学と違い$\vec{r}+\vec{r}Δt$となります。そのため、位置が$\vec{r}+\vec{v}Δt$で時間がt+Δt。従って速度$\vec{v}$は$v(\vec{r}+\vec{r}Δt, t+Δt)$になります。
これはオイラー表現というものが絡んでいるそうです。以下の記事2つを参考にしましたが、オイラー表現とは空間座標を固定しておいた上で、当該時刻にその位置にいる物質粒子を見る方法です。
http://gpvjma.ccs.hpcc.jp/~tanaka/ugomeku/lingo/euler_lagrangian.pdf
https://www.gfd-dennou.org/library/riron/renzoku/kijutu/pub/kijutu.pdf
そのため、位置は$\vec{r}$ではなく$\vec{r}+\vec{v}Δt$になります。
ナビエ・ストークス方程式は流体力学の支配方程式で、力学でいうところの$ma=F$になります。また、ナビエ・ストークス方程式を見ればこれが単位体積あたりのものだということがわかります。
また、方程式の下図で表されたところは$ma = F$の加速度に対応しそうです。
では加速度を算出してみます。
$$ \lim_{Δt \to 0}{\frac{\vec{v(\vec{r}+\vec{v}Δt, t+Δt)-\vec{v(\vec{r}, t)}}}{Δt}} $$
これを求めるわけですが、xyzと3成分になっていますので1方向に限定して記述します。ここで$\vec{r}=(x, y, z)$、$\vec{v}=(vx, vy, vz)$とします。
$$ \lim_{Δt \to 0}{\frac{vx(x+vxΔt + y+vyΔt + z+vzΔt, t+Δt)-v(x, t)}{Δt}} $$
これを計算するために、1次テイラー展開します。
aの周りでのオイラー展開は以下のようなものです。
$$ f(x) = f(a) + f^{\prime}(x-a)+\frac{1}{2!} f^{\prime\prime}(x-a){^2}+\frac{1}{3!} f^{\prime\prime\prime}(x-a){^3}+............ $$
Δtは微小量なのでvxΔt、vyΔt、vzΔtも微少量として展開します。
$$ \lim_{Δt \to 0}{\frac{vx(x+vxΔt + y+vyΔt + z+vzΔt, t+Δt)-v(x, t)}{Δt}} $$
もう一度展開前の式を見てみると最後にΔtを0に飛ばしているので、分子にΔtが残れば極限の段階でその項は0になります。分母にはΔtがいるのでテイラー展開後にΔtの二乗まで展開すると先の話で極限をとって0になるので一次までのテイラー展開で大丈夫です。
$$ \lim_{Δt \to 0}{\frac{vx(x+vxΔt + y+vyΔt + z+vzΔt, t+Δt)-v(x, t)}{Δt}} =\frac{∂vx}{∂t}+(\vec{v}・\nabla)vx $$
x成分に着目した時yやzに関係ない変数で表すことができたのでy、zにおいても同様の結果が得られます。 よってこれでナビエストークス方程式の左辺を出すことができ、同時に実装の際に使用する移流項を導出できました。
圧力項
移流項の時と同じように微小な六面体で考えると、六面体の右面を押す圧力はP(x+dx/2, y, z)になリます。
左面も同じように六面体の左面を押す圧力はP(x-dx/2, y, z)になリます。
では圧力項を出してみます。 $$ {P(x-\frac{dx}{2}, y, z) - P(x+\frac{dx}{2}, y, z)}dydz = (-\frac{dx}{2}\frac{∂P}{∂x}-frac{dx}{2}\frac{∂P}{∂x})dydz = -\frac{∂P}{∂x}dxdydz $$
dxdydzは体積Vになるので単位体積あたりの圧力は
$$ -\frac{∂P}{∂x} $$ になります。x成分と同じことがy、z成分でも計算されるため圧力項は
$$ -\nabla{\vec{P}} $$ になります。
外力項
外力に関してはその名の通り、流体に及ばされる力をいい、主にコリオリ力や重力を意味します。ここは流体に及ぶ力のベクトルの和をとっていけばいいと思います。
今回はUIやマウス操作などで流体にインタラクションしますので、その分の力などはここに寄与します。
粘性項
クエット流れ
以下で書いているように大きな板二枚を用意します。この上側の板を引っ張ると図で表しているように速度勾配は上の板の近付くほど大きくなります。
板の中にある流体を取り出した時、その微小区間での剪断応力を考えてると以下のような力が働いています。
ニュートン力学において実験的に得られている剪断応力は
$$ τ=μ\frac{∂vx}{∂t} $$
と表されます。
対称性
上記のように微小区間での六面体を考えます。この六面体において力のモーメントの釣り合いから
$$ τxy = τyx $$
として与えられます。
剛体回転
剛体回転を考えます。なぜ、ここで急に剛体を考えるかというと、剛体は大きさを持ち変形しない物体です。流体要素の位置関係が変化しないものなので$τ$で表される粘性抵抗が0になるはずだからです。
粘性抵抗が0になることから粘性抵抗の一般系を導き出したいという意図になります。
上記のような剛体の回転を考えるとき剛体の速度場は
$$ \vec{v}=\vec{Ω}×\vec{r}=(-Ωy, Ωx, 0) $$ となります。これに関しては以下の記事の3.5 固定軸の周りの回転の章を参考にしました。
https://www.sci.u-hyogo.ac.jp/material/electro_phys/member/kohara/text3.pdf
では粘性抵抗を考えてみると $$ τxy=μ(\frac{∂vx}{∂y} + a) (a:未知数) $$
となり、この式を恒等的に0にするアプローチとして
$$ τxy=μ(\frac{∂vx}{∂y} + \frac{∂vy}{∂x}) $$
とすればゼロになることがわかります。これが粘性応力の一般系です。
では導出のためにまずx成分の合力を考えます。
$$ \frac{∂τxx}{∂x} dxdydz + \frac{τyz}{∂y}dydxdz+\frac{τzx}{∂z}dzdydx = (\frac{∂τxx}{∂x}+\frac{∂τyz}{τy}+\frac{∂τzx}{∂z})dxdydz $$
となります。dxdydzは微小六面体の体積なのでさらに
$$ μ\frac{∂}{∂x}(\frac{∂vx}{∂x}+\frac{∂vx}{∂x})+μ\frac{∂}{∂y}(\frac{∂vy}{∂x}+\frac{∂vx}{∂y})+μ\frac{∂}{∂z}(\frac{∂vz}{∂x}+\frac{∂vx}{∂z}) = μ(\frac{∂{^2}vx}{∂x{^2}}+\frac{∂{^2}vx}{∂y{^2}}+\frac{∂{^2}vx}{∂z{^2}})+μ\frac{∂}{∂x}(\frac{∂vx}{∂x}+\frac{∂vy}{∂y}+\frac{∂vz}{∂z}) $$
この右辺の第2項は非圧縮性から導かれる$\nabla・\vec{v}=0$により0になります。
従って粘性項は
$$ μ\frac{∂{^2}{\vec{v}}}{∂x{^2}} $$
となります。
格子法と粒子法
粒子シミュレーションにおいて、大きく3種類に分別できるそうです。
・格子法
・粒子法
・格子法+粒子法
特徴として格子法は「格子状の場をつくり、時間で微分した際にその各格子がどういった速度になっているか」をシミュレーションし、粒子法は「粒子に着目し、粒子自体の移流をシミュレーション」する方法になっています。
格子法
参考Gitにも記載がありますが、流体のシミュレーションにおいて圧力、粘性、拡散等の計算は得意ですが移流の計算は苦手だそうです。
また、格子法の内容は以下の記事を参考にしました。
格子法では格子の配置の仕方は以下の3つ種類があるそうです。
・Arakawa-B型スタッガード格子
以下の配置で格子を定義する方法です。
速度に対する物体境界条件の設定が容易というメリットがあるようです。
・コロケート格子
以下の配置で格子を定義する方法です。
全て同じ格子点上で定義する方法はコロケート格子配置と呼ばれ、計算精度は高いが速度と圧力のカップリングが悪く、計算が不安定になり易い
・Arakawa-C型スタッガード格子
以下の配置で格子を定義する方法です。
速度と圧力を以下のように交互に定義する配置方法を(Arakawa-C型)スタッガード格子と呼び、標準的なスタッガード格子配置である。カップリングが良く、安定な計算できる
今回の格子
引用:Real-Time Fluid Dynamics for Games
今回の実装ではReal-Time Fluid Dynamics for Gamesは格子法、そして上記のコロケート格子を採用しています。格子の中心に圧力と速度を持たせる方法になります。
粒子法
こちらは次の記事で詳細に書こうと思います。
実装
以上でナビエストークス方程式に関して少しわかったところで、ではどうやって実装を行おうかということになります。
格子法における方程式
格子法におけるナビエ・ストークスの方程式は以下の3式を扱いました。
$$ \frac{∂\vec{u}}{∂t}=-(\vec{u}・\nabla)\vec{u}+ν\nabla{^2}\vec{u}+\vec{f} $$
$$ \frac{∂ρ}{∂t}=-(\vec{u}・\nabla)ρ+κ\nabla{^2}ρ+S $$
$$ \nabla・\vec{u} = 0 $$
一つ目は速度場について、二つ目は密度場、三つ目は質量保存則になります。
一つ目の速度場と二つ目の密度場はかなり似ている式であることは見てわかると思います。スカラー表現になっている点と動粘性係数が拡散係数κになってる点です。 こちらは参考のUnityGraphicsProgramming.Vol1を引用します。
密度場は必ずしも必要ではありませんが、速度場を求めた際の各ベクトルに対し、密度場で拡散された画面上のピクセルを載せることで、溶けながら流れるような、より流体らしい表現が可能になります。
Real-Time Fluid Dynamics for Gamesを参考にすると、煙の粒子などを全てモデル化することは難しく、コストがかかってしまいます。なのでそれを解決するため、全てモデル化する代わりに煙の密度場に変換するというアイデアを提案しています。
以下のCGにおける流体映像の効率的な生成に関する研究、においてもこの「スカラー場の解析方法」として密度場が格子法において取り上げられています。
CGにおける流体映像の効率的な生成に関する研究 [論文内容及び審査の要旨] : HUSCAP
質量保存則
参考Gitでは「連続の式(質量保存則)」として記述がある以下の式ですが、
$$ \nabla ・ \vec{v} = 0 $$
これは今回見る流体の非圧縮性を示しています。上記に挙げたナビエストークスの方程式は仮定として密度ρと粘性係数μが一定であるという非圧縮であり等質な流体に着目したものでした。 その非圧縮性から速度場の発散が0という式が導出されます。非圧縮性とは砕けた言い方で言うと「密度が場所に変わらず一定」と言うことになりますがそれを表した式になります。
上記の図で表したように微小区間で囲われたセルを考えます。このセルへの流入、流出を考えると
$$ \frac{i(x+Δx, y)Δy+i(x, y+Δy)Δx-i(x, y)Δy-i(x,y)Δx}{ΔxΔy} = \frac{i(x+Δx, y)-i(x, y)}{Δx}+\frac{i(x, y+Δy)-i(x, y)}{Δy} $$
ここから極限をとって
$$ {\lim_{Δx \to 0}\frac{i(x+Δx, y)-i(x, y)}{Δx}}+\lim\frac{i(x, y+Δy)-j(x, y)}{Δy} = \frac{∂i}{∂x} + \frac{∂j}{∂y} $$
となり、微小区間における流体の流出流入は発散で表されることがわかります。よって非圧縮性の流体では微小区間セル内の流入量と流出量はイコールになります。 そのため、流入量と流出量の総和が0 = 速度場の発散が0 と言う式になるというわけだそうです。
実装全体のフロー
実装は以下のようなフローになります。速度場も密度場も外力が関係するのでUIやマウス等のインテラクションによる力の発生を二つに伝えなければいけません。
処理順としては大きく、速度場から密度場に渡します。これは、密度場の計算での「密度の移流」で速度場の速度を用いるためです。そのため密度場よりも速度場のほうが早く処理させる必要がありそうです。
もう一つ取り上げたいのが密度場の計算に質量保存がないことです。これは密度場は速度場と違い非圧縮の必要がないためです。非圧縮性から導かれた質量保存ですから、非圧縮である必要がない密度場には処理を入れる必要はありません。
密度場の考え方
格子分けと変数定義
先に述べたように処理順は速度場から密度場ですが、先に密度場から説明します。これは、密度場で考えた処理のフローやアイデアが速度場にも応用が効くためです。
引用:http://graphics.cs.cmu.edu/nsp/course/15-464/Fall09/papers/StamFluidforGames.pdf
流体を先に述べたように格子状に分割します。圧力と速度をセルの中心で表すのがコロケート格子なのですが、この考えは実装にも役に立ちました。各セルの中心で流体をサンプリングすればいいのです。
RWTexture2D<float4> solver; //xy = velocity, z = density RWTexture2D<float>density; //density field RWTexture2D<float2> velocity; //velocity field RWTexture2D<float3> prev; //xy = prev velocity field, z = prev density field. when project step x = p, y = div Texture2D source;
このような考え方から上記のようにTextureやRWTexture2Dのような変数定義を行いid管理で格納していきます。
参考Gitでいうところで以下の表記になっています。
density[id] = (id.x == 0) ? density[id + uint2(1,0)] : density[id]; // idはuint2
引用:http://graphics.cs.cmu.edu/nsp/course/15-464/Fall09/papers/StamFluidforGames.pdf
これは上記の図で言うとdensity[uint2(1,2)]で以下で表すような意味合いです。
密度場のフロー
$$ \frac{∂ρ}{∂t}=-(\vec{u}・\nabla)ρ+κ\nabla{^2}ρ+S $$
密度場は先に示すように、上記の式で表されます。密度場は右辺の3つの要素で決めることができるということになります。
最初の項は密度場が速度場に影響されると言うこと、2項目は密度が割合的に拡散されることを示し、3項目はUI等のインプットから密度場が上昇することを示します。
引用:http://graphics.cs.cmu.edu/nsp/course/15-464/Fall09/papers/StamFluidforGames.pdf
実装では上記で示されるように右辺3項を逆の順番で処理させます。3項目のInput、2項目の拡散、1項目の速度場影響という感じで、さらにこれらを繰り返します。
外力項
外力は先に述べているように単純にInput(マウスの移動方向)ベクトルを重ね合わせるだけで処理します。
拡散項
引用:http://graphics.cs.cmu.edu/nsp/course/15-464/Fall09/papers/StamFluidforGames.pdf
拡散を考えるときはまず、上記のように最小の拡散影響範囲で考えます。自分のセルからみて、隣接するセル4つに自分のセルにある密度が拡散し(流出)、隣接するセル4つから自分のセルに密度が拡散(流入)します。
$$ denisity[i, j] = t * density[i-1,j] + t * denisity[i+1,j] + t * density[i,j-1] + t * density[i,j+1] - 4 * t * denisity[i,j] $$
という形でかけます。
指しているところが分かりづらいですが、上記のような形です。tは拡散率のようなもので仮で置いています。
しかし、この方法はうまくいきません。数式通りに実装するとシミュレーションが発散してしまうのです。
この理由としては実装において粘性係数・微分時間・格子数を乗算させた拡散率が高くなってしまうためです。
if (id.x < w && id.y < h) { float a = dt * visc * w * h; [unroll] for (int k = 0; k < GS_ITERATE; k++) { density[id] = ( prev[id].xy + a * (velocity[int2(id.x - 1, id.y)] + velocity[int2(id.x + 1, id.y)] + velocity[int2(id.x, id.y - 1)] + velocity[int2(id.x, id.y + 1)])) / (1 + 4 * a); SetBoundaryVelocity(id, w, h); } }
for文で書いている部分はシミュレーション発散を回避するやり方なのですが、その上の
float a = dt * visc * w * h
が大きくなりすぎるということです。 viscは動粘性係数、dtは微小時間、w*hは格子数に対応します。
dtにはTime.deltaTimeを渡しますがそれよりもwhviscの寄与が多い影響になります。
そのため、拡散を安定させる方法として定常反復法を使用します。
定常反復法
定常反復法について参考にしたのが以下二つになります。
https://ocw.tsukuba.ac.jp/wp/wp-content/uploads/2018/10/compmath2-2018-slide-04.pdf
http://nkl.cc.u-tokyo.ac.jp/13n/SolverIterative.pdf
定常反復法には大きく3つ以下のようになります。
・ヤコビ法
・ガウスザイダル法
・SOR法
参考のGit、Real-Time Fluid Dynamics for Gamesではガウスザイダル法を採用しています。
以下の記事を参考にしましたが、(Pythonでの計測ですが)SOR法の処理が一番早く、その次にガウスザイダル法、ヤコビ法と続くみたいです。 qiita.com
どこまでシミュレーションの速さに寄与するかは調べたいと思います。
また、
density[id] = ( prev[id].xy + a * (velocity[int2(id.x - 1, id.y)] + velocity[int2(id.x + 1, id.y)] + velocity[int2(id.x, id.y - 1)] + velocity[int2(id.x, id.y + 1)])) / (1 + 4 * a);
で、最後に(1+4*a)で1を足したもので割っているのは0で割らないためになります。
密度場への移流
引用:http://graphics.cs.cmu.edu/nsp/course/15-464/Fall09/papers/StamFluidforGames.pdf
移流も実装に工夫が必要です。上記引用図のaで示されるような速度場あったとき、この格子は速度をセルの中心で表すのでbのようにセルの中心に粒子があるように考えることができます。
ただ、今度は動いた後粒子がグリッド値に戻す必要があります。その解決策として、逆に1ステップ前にセルの中心にいる、現在グリッド上にいる粒子を見つけるというものがあります。(図c)
この粒子が運ぶ密度量はまず二つのグリッドを考えます。一つ目は前のタイムステップの密度値をもち、もう一つは最新の密度値をもちます。最新の密度値を持っているグリッドセルごとにセルの中心位置を速度の逆方向にバックトレースします。
次に、1ステップ前の密度値をもつグリッドセルの密度値と線形補完させてこの値を現在のグリッドセルの密度値にします。
つまり、最も近い 4 つの隣接点から開始位置の密度を線形補間することによって粒子が運ぶ密度量が得られます。
参考Gitのコードは以下のようになっています。
if (id.x < w && id.y < h) { int ddx0, ddx1, ddy0, ddy1; float x, y, s0, t0, s1, t1, dfdt; dfdt = dt * (w + h) * 0.5; //バックトレースポイント割り出し. x = (float)id.x - dfdt * velocity[id].x; y = (float)id.y - dfdt * velocity[id].y; //ポイントがシミュレーション範囲内に収まるようにクランプ. clamp(x, 0.5, w + 0.5); clamp(y, 0.5, h + 0.5); //バックトレースポイントの近傍セル割り出し. ddx0 = floor(x); ddx1 = ddx0 + 1; ddy0 = floor(y); ddy1 = ddy0 + 1; //近傍セルとの線形補間用の差分を取っておく. s1 = x - ddx0; s0 = 1.0 - s1; t1 = y - ddy0; t0 = 1.0 - t1; //バックトレースし、1step前の値を近傍との線形補間をとって、現在の速度場に代入。 density[id] = s0 * (t0 * prev[int2(ddx0, ddy0)].z + t1 * prev[int2(ddx0, ddy1)].z) + s1 * (t0 * prev[int2(ddx1, ddy0)].z + t1 * prev[int2(ddx1, ddy1)].z); SetBoundaryDensity(id, w, h); }
速度場の実装
$$ \frac{∂\vec{u}}{∂t}=-(\vec{u}・\nabla)\vec{u}+ν\nabla{^2}\vec{u}+\vec{f} $$
速度場は上で示したようにこの式で表されます。密度場と同じように考えると三項目のInputからの速度場上昇、2項目の粘性拡散、1項目の自己移流によって速度場がきまるということです。
実装は密度場の考え方を流用できます。ただ、密度場と違い、質量保存が成り立ちます(非圧縮性)。
これを実装するためのフローとしては、速度場から発散の算出、求めた発散からPoisson方程式をガウス・ザイデル法で解く、$∇・u = 0$を行うというフローで実装します。
境界条件
流体は壁から流れ出してはいけないため、垂直方向の壁での流体の水平成分の速度/密度は0、水平方向の壁での流体の水平方向の速度/密度は0にならないといけません。
または、0ではなくて逆側から出るようにしたりする式を実装することもできます。
参考Gitでは以下のように速度、密度ともに境界処理が入っています。
void SetBoundaryDensity(uint2 id, uint w, uint h) { density[id] = (id.x == 0) ? density[id + uint2(1,0)] : density[id]; density[id] = (id.x == w-1) ? density[uint2(w-2, id.y)] : density[id]; density[id] = (id.y == 0) ? density[id + uint2(0,1)] : density[id]; density[id] = (id.y == h-1) ? density[uint2(id.x, h-2)] : density[id]; density[id] = (id.x == 0 && id.y == 0) ? 0.5 * (density[uint2(1,0)] + density[uint2(0,1)]) : density[id]; density[id] = (id.x == 0 && id.y == h-1) ? 0.5 * (density[uint2(1,h-1)] + density[uint2(0,h-2)]) : density[id]; density[id] = (id.x == w-1 && id.y == 0) ? 0.5 * (density[uint2(w-2,0)] + density[uint2(w-1,1)]) : density[id]; density[id] = (id.x == w-1 && id.y == h-1) ? 0.5 * (density[uint2(w-2,h-1)] + density[uint2(w-1,h-2)]) : density[id]; }
最後に
ここまでで流体シミュレーションの実装に入るまでの「ナビエストークスの方程式」と「格子法」、そして実装の大まかな流れについてまとめれたので次はComputeShaderを使いながらさらに深い実装に入っていきたいと思います。
せっかくアドカレで始められたので、これで終わらず実装のところまでしっかりと書き続けたいと思います。