Simulating Fluids, Fire, and Smoke in Real-Time

https://scrapbox.io/files/6951047ebc5ab2f4c03c35d0.png

炎を目的とした流体計算のチュートリアル

1.1

流体の空間を2Dとしてその空間の座標をx,流体の速度をu(x,t)とする。速度場は格子法で与える。

ここで流体上を流れる何らかの物理量をスカラー場Ψ(x,t)と表現し、流体上でこれがどのように輸送されるか(移流 advectionと呼ぶ)をシミュレーションする。

ナイーブに考えるとスカラー量は流体の速度場uによって移動されるのであるため、時刻Δt後のスカラー場には次の方程式が成り立つことが考えられる

https://scrapbox.io/files/659d0000557d4e0024406557.png

格子法では同じターゲット点が対象とされるし、そもそも移動先が格子点と合わない可能性があるのでこれをそのまま使うのは宜しくない（多分粒子法で解くんですかね?)。なので偏微分方程式(PED)を基に計算を考える。

第一法則として、ある点を注目してスカラー量が出ていく量とその点から消えた量というのは必ず等しいはずである。従って、ある固定の閉じた領域Wと境界Sを用いて以下の様に表される「質量保存の法則」が成り立つ。

https://scrapbox.io/files/659d01ca52aaee0025017562.png

この時、ガウスの発散定理からS面微分から立体積分の変換を行えば、次のような方程式が成立する。

https://scrapbox.io/files/659d02418279ee002aa9cdd5.png

これは極小のWでも成立するため、どの点においてもこの方程式が成立する(いつも思うけどこれ数学的に保証されるのかしら)。

https://scrapbox.io/files/659d0282701cfe0025c41068.png

よってスカラーの偏微分方程式として

https://scrapbox.io/files/659d02e9d4ea29002411b008.png

この方程式は非圧縮流れ移流方程式(incompressible flow advection equations)と呼ぶ。スカラー以外にも何らかのベクトル場vでも同様に導出される。

https://scrapbox.io/files/659d03cd54919c0025bf698c.png

一旦、uは時間に依存しないものとして考え、スカラー場の移流を考える。

https://scrapbox.io/files/659d0426bcd8140023b2b5ef.png

この方程式は次のような形で近似的に計算する方法がある。これはセミラグランジュ法と呼ぶ。（どうやら1次までの近似らしい、他にもCIP法とか色々あるらしい）

https://scrapbox.io/files/659d05d4525a200025386f3b.png

実装上は今いる場所の速度ベクトルを用いて前の位置x-uから値をその場所のスカラー量とすれば良い。

1.2 ナビエストークス方程式

uは現実的には時間に依存するはずである。このuが時間依存する際に従う方程式がナビエストークス方程式である(非圧縮流体が前提)。

https://scrapbox.io/files/659d06af7464910026b29085.png

μは流体の粘度を表す定数、Fは外力、pは多分圧力?(後で不要になるので特に書いてない)。ここでは流体は非粘性、外力はないものとして次のように考える。

https://scrapbox.io/files/659d0737099bf700240f2dc9.png

自己移流項 Self-advection

速度場の移流項である。これ自体はスカラー量のやつと同様に考えれば良い

https://scrapbox.io/files/659d07b791ddc700231e35f5.png

ここで移流後の速度場をu'と表す。(後で使う)

https://scrapbox.io/files/659d07c2e505390023ed25f4.png

圧力項 Pressure

非圧縮流体では圧力と速度の間には次の方程式が成り立つらしい。

https://scrapbox.io/files/659d09c0525a20002538d31f.png

ここから、次のようなポアソン方程式を導くことが出来る

https://scrapbox.io/files/659d0a6e2786be00232d31f1.png

なので、わざわざpを設定しなくてもu'から得ることが出来る。ただし、どうやらこの方程式を数値的に解くのは大変な事らしい。

右辺は以下の様に近似可能(1)

https://scrapbox.io/files/659d0b288279ee002aaac323.png

ラプラシアンも次のように書ける(2)

https://scrapbox.io/files/659d0b8da934c30024fde4cf.png

よって５つの未知数を持つ１次方程式としてポアソン方程式を考えることが出来る。ここで、グリッドNxN分だけこれら方程式を考えることによって、NxN連立方程式として圧力pを求めることが出来る。

Ap = b

なので連立方程式を高速に解きましょうねという話になる。GPU上で計算するのとしてはヤコビ法が簡単で良いとのこと、

ヤコビ法は連立方程式に対して、何度か反復して計算することで任意精度の近似を行う方法である。こうゆうx_1に対して方程式があった時、

https://scrapbox.io/files/659d0d1f557d4e002441c53c.png

k+1回目の反復時に計算する値は次のように求める。

https://scrapbox.io/files/659d0dbb557d4e002441d2f7.png

(1)と(2)を見比べれば、bというのは(1)式、A_11というのは-4, A_~は1として当てはめ、上下左右のAのみ考えれば良い。実装としては次のように書ける。

https://scrapbox.io/files/659d0e5e3807280022b0b183.png

これを何度かやるとpが得られるので、そこからgradientを求めれば圧力項pを得られる。

実装

全体的な実装としては次のよう

https://scrapbox.io/files/659d1026f8737200238bf026.png

Unityで実装するなら

advection計算 u'レイヤー

イテレーション k

p計算レイヤー(k枚)

Density更新レイヤー