一、前言

游戲中的場(chǎng)景交互一直是人們追求的，比如人物與水的交互、物體與煙霧交互、可交互雪地等。它們都有一個(gè)共性，就是這些物體都是具有流體效果的（雪地可以看做是流體粘度無(wú)限大）。因此我想嘗試一下流體的模擬是如何做出來(lái)的。這個(gè)難度是遠(yuǎn)高于預(yù)期的，其中充斥著無(wú)數(shù)的數(shù)學(xué)方程。與此同時(shí)，將這些數(shù)學(xué)公式離散化轉(zhuǎn)換為分步求解的代碼也是一大挑戰(zhàn)。

我看了網(wǎng)上很多的文獻(xiàn)以及許多大神的做法，但大神們普遍言簡(jiǎn)意賅，而有些在求解過(guò)程中出現(xiàn)了不必要的求解方式。我則在這些基礎(chǔ)上，根據(jù)自己的理解整理了一份模擬思路。

然后我簡(jiǎn)要介紹一下NS方程，這個(gè)方程是流體模擬的關(guān)鍵。NS方程全程是納維-斯托克斯偏微分方程組（Naiver-Stokes），使用這個(gè)方程基本可以描述各種流體在各種情況下的運(yùn)動(dòng)方式。求解三維空間中的N-S方程組光滑解的存在性問(wèn)題是千禧年十大問(wèn)題之一，它的難度可想而知。下面展示的是不可壓縮流體的NS方程

?

在前言的最后，我們要明確以下幾個(gè)問(wèn)題，首先NS方程的求解是舍棄了很多項(xiàng)，并且簡(jiǎn)化很多條件后的數(shù)值解；我們始終是在2D維度求解NS方程；我所解釋的更為具體化，盡管我已經(jīng)盡可能的求證，但也可能存在錯(cuò)誤的解釋；下文涉及到的公式我都盡可能以普通高數(shù)水平和較高高中數(shù)學(xué)水平進(jìn)行解釋。一定要注意區(qū)分方程中矢量場(chǎng)和標(biāo)量場(chǎng)。

?

二、NS方程推導(dǎo)

矢量微積分算子

不可壓縮流體的NS方程中，很多量都是矢量，因此對(duì)于矢量的計(jì)算，我們需要引入一些算子，用于簡(jiǎn)化表示矢量的計(jì)算。下述的都是2維的算子公式。

梯度算子 ▽ ：，物理含義是矢量場(chǎng)在各個(gè)方向上的偏導(dǎo)的矩陣。更通俗來(lái)說(shuō)，就是一個(gè)向量場(chǎng)上的某個(gè)點(diǎn)，在此位置的某個(gè)方向上的變化率，以速度場(chǎng)求梯度▽u，就可以理解為某點(diǎn)在x方向的加速度和在y方向的加速度所組成的向量。

散度算子 ▽· ：，物理含義是矢量場(chǎng)在此點(diǎn)偏導(dǎo)的求和，得到的結(jié)果是一個(gè)標(biāo)量。也就是說(shuō)，散度是當(dāng)前點(diǎn)的物質(zhì)增加還是減少。

旋度算子 ▽× ：，物理含義是矢量場(chǎng)在此點(diǎn)偏導(dǎo)的差，得到的也是一個(gè)標(biāo)量。通俗來(lái)說(shuō)就是此點(diǎn)的環(huán)流情況，也就是是否會(huì)產(chǎn)生渦旋現(xiàn)象。

拉普拉斯算子 ?2 = ▽ · ▽：，拉普拉斯算子不太好進(jìn)行直觀的解釋，大致來(lái)說(shuō)它提供了標(biāo)量場(chǎng)的曲率和平滑度的信息。

?

重要的一點(diǎn)，速度場(chǎng) = 速度的散度場(chǎng) + 速度的旋度場(chǎng)

?

?

基礎(chǔ)物理公式

首先先列出幾個(gè)最基礎(chǔ)的物理公式：

動(dòng)量方程：p = m*v

動(dòng)量方程另一形式：F = p / t

牛二定理：F = ma

密度公式：ρ = m/V

這些基礎(chǔ)公式有了，再非?？焖俚慕榻B一下歐拉方法和拉格朗日方法，兩個(gè)方法都是假設(shè)空間中有一個(gè)微元i，歐拉法相對(duì)參考系是整個(gè)世界，拉格朗日法相對(duì)參考系是以微元自身。我們將使用歐拉方法進(jìn)行NS方程的推導(dǎo)，下面開(kāi)導(dǎo)：

動(dòng)量守恒

現(xiàn)在假設(shè)一個(gè)方程，u = u(t, X)，它用于描述一個(gè)粒子在時(shí)間 t 和位置X時(shí)的速度方程。方程u對(duì)時(shí)間求導(dǎo)，物理意義為一個(gè)二維平面上的點(diǎn)元在t時(shí)的加速度。這個(gè)式子被稱為拉格朗日變量。注意這里是點(diǎn)乘而不是乘積。

?

x是二維的，?X/ ?t 就是該點(diǎn)的速度 u，然后我們將?u / ?X = (?u / ?x , ?u / ?y) 記作▽u，得到：

?????????（1）

?

再換一個(gè)角度，考慮一個(gè)質(zhì)點(diǎn)受到的力，F(xiàn)總 = ?接觸力+黏性力（剪應(yīng)力），這里是根據(jù)牛頓流體進(jìn)行推導(dǎo)的。（牛頓流體（英語(yǔ)：Newtonian fluid）指應(yīng)力與應(yīng)變率成正比的流體。此比例系數(shù)為流體的黏度。）

?

就可以寫(xiě)成 F總 = F(壓強(qiáng)) + τ + F(額外力)，我們分別求出F(壓強(qiáng)) 和 τ

F(壓強(qiáng)) = ?VdP / dX = (?p / ?x, ?p / ?y)。我們將?P/ ?x記作▽p，得到：

F(壓強(qiáng)) = V▽p ????????????????（2）

然后請(qǐng)看牛頓流體特性的最后一句話，因?yàn)檫@里的速度場(chǎng)是矢量場(chǎng)，因此需要先對(duì)速度場(chǎng)求梯度，得到方向速率變化，得到的值點(diǎn)乘梯度才能得到速度場(chǎng)在各方向的梯度，得到：

τ = μ*▽·(▽u)= μ*?2u ??????（3）

下面，我們綜合（1） = （2）+ （3）+F(額外力)

m(?u / ?t + u·▽u) = V▽p + ?μ*?2u + f ， 兩邊同時(shí)除m 并移項(xiàng)

=> ??u / ?t = -u·▽u - 1/ρ* ▽p + v*?2u + F???(常數(shù)項(xiàng)v = μ/m， F = f/m)

?

質(zhì)量守恒：

m0 - m1 = ρ*△V = u*△t*ρ =ρ* ▽ · u = 0

=>?▽ · u = 0 ?

?

到此為止，令人驚嘆的，我們推導(dǎo)出來(lái)了不可壓縮流體的NS方程！

?

三、數(shù)學(xué)方法求解方程

現(xiàn)在我們假設(shè)一個(gè)情況，我們有一個(gè)當(dāng)前時(shí)刻的速度場(chǎng)，那么在接下來(lái)的一段時(shí)間內(nèi)，這個(gè)場(chǎng)將如何變化呢？也就是說(shuō)，我們已知u，如何去求解上面的NS方程，得到加速度，進(jìn)而得到下一刻的速度場(chǎng)。以u(píng) = at為例，這是連續(xù)的方程，由于是在電腦中模擬，我們還需要離散化。

很明顯這是一個(gè)歐拉視角的求解。

亥姆霍茲-霍奇分解

用公式表示為

這個(gè)在講解算子的時(shí)候說(shuō)到了，任何向量場(chǎng)都可以分解為兩個(gè)其他向量場(chǎng)的和，

速度場(chǎng) = 速度的散度場(chǎng) + 速度的旋度場(chǎng)。

回看一下動(dòng)量守恒項(xiàng)，我們簡(jiǎn)化一下NS方程，假設(shè)這個(gè)流體是無(wú)粘度的理想液體，公式的黏性項(xiàng)v*?2u 就可以舍去了：

?

也就是說(shuō)力只剩下了F = -1/ ρ ▽p

又有：F = ma = m * (u(n+1)-un) / △t,帶入得：

??????（1）

又因?yàn)椋ㄖ罢f(shuō)到的質(zhì)量守恒公式）：

（2）

?

現(xiàn)在令△x = △y = △t = 1，將（1）帶入（2）得：

?

?

這個(gè)式子被稱為壓力泊松方程(Pressure Possion Equation)

?Δt/ρ???p=?Δ??u

?

四、轉(zhuǎn)換為代碼

上面主要是說(shuō)了如何將▽p轉(zhuǎn)換為關(guān)于u的離散化方程，速度場(chǎng)也是要迭代的，這個(gè)迭代就是平流項(xiàng)。那么大致來(lái)說(shuō)呢，先初始化一下速度場(chǎng)，求個(gè)平流項(xiàng)（平流項(xiàng)是速度場(chǎng)點(diǎn)乘速度場(chǎng)的梯度），然后加一個(gè)壓力場(chǎng)（壓力場(chǎng)是通過(guò)速度場(chǎng)迭代求解出來(lái)的），最后再保證一下速度場(chǎng)的散度為0。就這幾步就ok了！

1、初始化一個(gè)速度場(chǎng)和密度場(chǎng)（也可以叫染料場(chǎng)）。

兩個(gè)場(chǎng)的的區(qū)別只有顏色是不同的。下圖是我自己提前畫(huà)了一個(gè)速度場(chǎng)存儲(chǔ)在了RenderTexture中，（這一篇我希望更偏學(xué)術(shù)向，我不再介紹什么是RT，什么是shader等等了。不懂的可以簡(jiǎn)單理解為RT就是圖片，shader就是如何計(jì)算像素（體素）的。）為什么顏色是這樣的呢，畫(huà)圖解釋：顏色是rgb，我們只取rg通道作為速度向量的xy，黑色是由于負(fù)數(shù)，由于0也會(huì)顯示黑色，造成了顏色變成了這樣。

?

模擬后可以發(fā)現(xiàn)確實(shí)如我們所想的，鼠標(biāo)初始化速度場(chǎng)也就是這么做的。密度場(chǎng)相當(dāng)于對(duì)有速度場(chǎng)的區(qū)域上色，比如空氣流動(dòng)是看不到的，但我們噴一些煙上去，也就看到了整個(gè)速度場(chǎng)的形狀。換成鼠標(biāo)上色的方式。

2、計(jì)算平流項(xiàng)（對(duì)流項(xiàng)）。

平流項(xiàng)也就是計(jì)算：-（u·▽）u = -（▽u）·u，速度場(chǎng)剛剛初始化了，對(duì)它求梯度再相乘就好了，這里能轉(zhuǎn)化是因?yàn)樘荻人阕訚M足分配律、分解率，但不滿足交換律。那么對(duì)速度場(chǎng)求梯度，我們只需要代入之前的離散化方程就可以了。

看代碼，我們?cè)趕hader中對(duì)每一個(gè)像素點(diǎn)進(jìn)行平流，簡(jiǎn)化代碼只剩下平流項(xiàng)，畫(huà)圓圈可以發(fā)現(xiàn)顏料都好像是沿著切線方向移動(dòng)了。看shader代碼，先來(lái)一個(gè)簡(jiǎn)化版本的，看公式這不就是求平均值嘛！。采樣速度場(chǎng)，得到的是當(dāng)前點(diǎn)的速度向量，x = vt, 采用前向歐拉法，i.uv=> i.uv-△x(i - (i-1))。運(yùn)行一下，沒(méi)問(wèn)題！But！注意步長(zhǎng)，我們選的很小，如果稍微大一點(diǎn)呢？

?

爆炸原因是：零空間問(wèn)題，有興趣可以自己搜索一下，簡(jiǎn)單解釋就是這種方法是有偏的，并且精度是O(△x)。短步長(zhǎng)會(huì)讓離散解較為接近，增加步長(zhǎng)后由于沒(méi)有考慮中間時(shí)刻的變化，導(dǎo)致數(shù)值解出現(xiàn)高頻振蕩。反映在屏幕上，遠(yuǎn)處還好一點(diǎn)，鼠標(biāo)位置爆炸最嚴(yán)重。這就是在鼠標(biāo)地方速度有一個(gè)突變。

看來(lái)需要一種好一些的方法來(lái)計(jì)算步長(zhǎng)。這里我們引入一種網(wǎng)格結(jié)構(gòu)：MAC網(wǎng)格，速度分量記錄在網(wǎng)格邊界上，然后我們采用一種半拉格朗日法（拉格朗日視角解決歐拉視角的問(wèn)題），保證了函數(shù)解最終會(huì)收斂。GPUGems圖中展示了一個(gè)向量如何反向求解，得到的是質(zhì)點(diǎn)在△t后的新的速度。為什么要雙線性差值，向量很大概率不指向uv中心，而雙線差實(shí)際上就是對(duì)周圍取均值，精度也提高到了O(△x^2)。并且好就好在，tex2D函數(shù)，或者UE中的Sample2D采樣Texture都是先二線差再返回顏色，因此從代碼角度來(lái)說(shuō)簡(jiǎn)化了我們的編程。

?

coord對(duì)應(yīng)實(shí)際應(yīng)該偏移的距離（目標(biāo)位置的uv坐標(biāo)），將這個(gè)點(diǎn)對(duì)應(yīng)到MAC網(wǎng)格上，求MAC網(wǎng)格附近線性插值后的顏色。

?

3、求壓力項(xiàng)▽p

亥姆霍茲-霍奇分解，得到的公式，這個(gè)公式很明顯是一個(gè)迭代的形式，這個(gè)就被稱為雅可比（Jacobi）迭代。

這一項(xiàng)不就是▽·u嘛！因此我們先求出速度場(chǎng)的散度，再代入散度值進(jìn)行迭代。乘0.9999是浮點(diǎn)數(shù)精度問(wèn)題。

?

還有一種迭代法叫高斯-賽德?tīng)柕?，也很不錯(cuò)，而且更為易懂，我就不詳解了。

4、組合起來(lái)

上面的步驟我們得到了迭代后的壓力場(chǎng)，還有速度場(chǎng)，我們?cè)賾?yīng)用一下NS方程的第二項(xiàng)就是那個(gè)不可壓縮項(xiàng)（u無(wú)散度 = u-p），最終得到的是無(wú)散度速度場(chǎng)，將這個(gè)場(chǎng)應(yīng)用到密度場(chǎng)上，就可以了。

舉一反三，如果把速度擴(kuò)散項(xiàng)設(shè)置為比較小的值，那么是不是可以做一個(gè)雪地效果，或者是書(shū)法筆跡？如果在最后不對(duì)密度場(chǎng)平流，而對(duì)壓力場(chǎng)平流，好像水波效果有了哦?。梢钥匆幌耮ames103 P10課程，思路是一樣的。）再想想風(fēng)場(chǎng)呢……這些只是一個(gè)大概思路，有興趣的可以繼續(xù)探索哦！

5、擴(kuò)展計(jì)算 ???粘度項(xiàng)（渦度）

把他放在最后是因?yàn)槠鋵?shí)目前的模擬效果已經(jīng)可以接受了，而從原理上來(lái)講，我們之前計(jì)算壓力項(xiàng)時(shí)，實(shí)際是舍棄了這一項(xiàng)，現(xiàn)在我們將這一項(xiàng)加回來(lái)。

還記得這個(gè)公式吧，剪應(yīng)力的表達(dá)式，仔細(xì)看看，你會(huì)發(fā)現(xiàn)為什么這么巧，這不就是對(duì)速度場(chǎng)求旋度嘛

我們令 ：

w = ▽×u

則 f = ε (N×w)， 其中N = ▽|w| / （| ▽|w| |），w就表示渦量，對(duì)渦量歸一化再乘w的大小，實(shí)質(zhì)就是讓渦量方向統(tǒng)一。求得的f = ma = mvt，m = 1，則△v = f*△t。讓速度場(chǎng)疊加上去就可以了。

?

五、總結(jié)

用物理推導(dǎo)公式 -》從公式中可以看出來(lái)我們要求的是平流項(xiàng)、壓力項(xiàng)、擴(kuò)散項(xiàng)、質(zhì)量守恒項(xiàng) -》離散化后發(fā)現(xiàn)很多項(xiàng)和提到的算子息息相關(guān) -》 離散化求解平流項(xiàng)用半拉格朗日法求 -》 壓力項(xiàng)先求散度再用雅各比迭代求解 -》 粘度項(xiàng)先求旋度再歸一化求解 -》 得到的結(jié)果還需要滿足質(zhì)量守恒 -》最后將速度場(chǎng)應(yīng)用到密度場(chǎng)就行了。

?

六、參考文獻(xiàn)

https://yangwc.com/2019/05/01/fluidSimulation/

https://zhuanlan.zhihu.com/p/374280832

https://zhuanlan.zhihu.com/p/270531017 https://developer.nvidia.com/gpugems/gpugems/part-vi-beyond-triangles/chapter-38-fast-fluid-dynamics-simulation-gpu https://www.bilibili.com/video/BV1i54y1F73W/?spm_id_from=333.788&vd_source=e883ef8d644cde54ca3470fb7bdc548d

?

源碼地址： https://github.com/clatterrr/FluidSimulationTutorialsUnity

?

b站 邱丑丑啊