【遊戲流體力學基礎及Unity代碼(四)】用歐拉方程模擬無粘性染料之公式推導

先放一張動態圖吊一下胃口~下面就是最終的效果

不可壓縮的歐拉方程只比NS方程少一個粘性項。所以下面的內容是完全適合NS方程的。各位請準備好！

模擬流體的時候會遇到許多數學公式。為了深刻理解這些數學公式，我們先從最簡單的複習起，比如說二十以內的加減法。

如下圖，ABCD四個區域原本有一些人，此時刻這些人的離開速度為3xi + yj，也就是說位於（Xa = 1，Ya = 1）的A區域的人將有3 * Xa個人從x軸正方向離開A區域進入B區域，將有1 * Ya個人往y軸正方向離開A區域進入C區域。那麼經過一個時刻後，一共有多少人離開了ABCD四個小區域組成的大區域？請寫出這道題的兩種解法，其中第一種方法分別計算四個小區域的離開人數，第二種方法則把四個區域當成一個整體。

如果你正確寫出了兩種解法，那麼趕緊給自己戴一朵小紅花吧！如果沒有也沒關係，我們來看看這兩種解法的不同之處。如下圖，第一種解法，由十以內的加減法顯而易見，ABCD四個小區域人數的淨離開量都是4，相加得到總淨離開人數是16。第二種解法，我們只計算從大區域邊界上離開的人數，畢竟對於小區域之間的人互相移動並不影響大區域的人數離開量。最後結果也就是2+2+6+6=16。

於是我們得出一個一般的結論，對二維空間來說，如果將一個大區域分成很多小區域，那麼

從大區域的邊界上的總離開量 = 各個小區域的淨離開量相加

上面這個公式還有一個聽起來高大上的名字，散度定理(divergence theorem)，或者高斯定理(Gauss's theorem)。所以畢導說很多東西小學二年級我們就學過真的不是瞎吹的。

上面這個公式寫成積分形式就是

其中F是我們的向量場，在我們這道題中就是速度。n就是邊界的法向量，l是邊界的長度，倒三角就是向量場F的梯度。C就是curve曲線，也就是大區域的邊界，S是surface表面，也就是整個大區域的面積。至於那個積分符號，在計算機離散的世界裡，把它看成加號就行了。

然後三維散度定理，說人話就是，將三維坐標系上的一個三維大空間分成許多三維小空間，那麼

從大空間的邊界上的總離開量 = 各個小空間的淨離開量相加

這個公式的積分形式在各大教科書上長這樣，右邊那個V是指Volume體積的意思：

散度定理在本系列中主要用於推導動量方程。

現在我們來研究固體和流體有什麼區別。

比如說一塊磚頭，第一種運動是平移(Translation)，也就是從一個地方到另一個地方，但不改變朝向。第二種運動是旋轉(Rotation)。這兩種情況下磚頭的形狀都不發生改變，是不可壓縮的(incompressible)，因此也叫剛體(Rigid Body)。

但對於一隻貓來說，情況就不同了。它不僅可以平移旋轉，還可以任意改變自己的形狀，這也就是為什麼有「一杯貓」，「一灘貓」，「一盒貓」的說法。

因此研究人員Marc-Antoine Fardin用流變函數證實了廣為流傳的說法「貓是液體」，或者更準確的「貓是可壓縮流體」，也因此獲得了2017年搞笑諾貝爾物理學獎。

因此為了研究貓，或者說是可壓縮流體，我們先來看看它能做什麼類型的運動。如下圖了，除了平移和旋轉外。還包括線變形(linear deformation)和角變形(Angular Deformation)。

如果一個力垂直於可壓縮物體的表面的話，那麼這個力就會造成線變形。比如用力將枕頭壓扁，如果這個力與表面的切線方向一致，那麼這個力就會造成角變形。

比如對滑鼠滑輪來說，如果用力向下按，而不讓滑輪滑動，那麼這個力算垂直於滑鼠滑輪的切線方向。如果前後滑動滑輪，那麼就算是與表面切線方向一致。

對於由流體粒子組成流體微團來說，各點的平移速度不應該和坐標有關係，也就是一個常數，也就是

而各個點的線變形率應該是一個一次方程，也就是

對於角變形來說，角變形率用如下的方法求：

如上圖，當由橙色正方形變成藍色平行四邊形的時候，取一個非常短的時間，那麼點B變成了點B星，它在y軸上的速度為v，在x軸上的速度很小所以忽略不計。點D點在x軸上的速度為u，在y軸上的速度很小所以忽略不計。

所以角度B星AB上的的正切值在單位時間上的變化量應該為y軸上的速度除以AB之間的距離

又因為AB之間的距離很小，夾角B星AB也很小，所以得到下面這個式子。點D同理。

所以我們要求的第三種速率，也就是角變形速率如下：

然後對於旋轉來說，如下圖：

B在y軸上的速度為v，在x軸上的速度很小所以忽略不計。點D點在x軸上的速度為u，在y軸上的速度很小所以忽略不計。由於線速度等於角速度乘以半徑，那麼B點和D點的角速度：

因為B點和D點的角速度實際上是相等的，所以當我們考慮流體微團的x軸速度u和y軸速度v與流體微團整體角速度的關係時，會得到下面這個式子，也就是我們要求的第四種速度，旋轉角速度：

等式最左邊括號裡的式子，就是旋度(curl)的一部分。符號如下。由上式看出，旋度剛好等於角速度的兩倍。如果流體的速度的旋度等於零，那麼流體微團的角速度也是零，那麼這個流體就是無旋的(irrotational)。比如說如果地球僅僅繞太陽公轉，而自身不旋轉的話，那麼地球的運動就是無旋的。但如果地球不僅公轉，還自轉的話，那麼地球的運動就是有旋的。旋度在後面我們討論渦流的時候還會遇到。

研究可壓縮的「一灘貓」非常麻煩，所以我們先從不可壓縮的流體研究起。之前我們說過，我們的不可壓縮流體只能有平移和旋轉，而不能有線變形和角變形，所以我們把這兩種速度分開，也就是

速度  = （平移和旋轉） + (線變形和角變形)

寫成二維微分形式就是下面這樣。

這個也叫亥姆霍茲速度分解(Helmholtz decomposition)。在張量分析中，這叫把一個張量分解成一個反對稱張量和一個對稱張量。這個等式的性質可以看看《A Mathematical Introduction to Fluid Mechanics》. 3rd ed. Springer.Chorin, A.J., and J.E. Marsden. 1993.的1.3章。

等式右邊第一項就是我們想要的，它的正對角線上兩項都是0，代表平移，反對角線上就是它的旋轉角速度。它們互為相反數，這是旋轉導致的。它也符合無散的性質。

等式右邊第二項是我們不想要的，它的正對角線是線變形，而反對角線上的兩項是一樣的，它是角變形導致的。第二項符合無旋的性質。我們之後會用某個方法把第二項去掉。

我們再回顧一下簡單的加減法。例如下圖左邊16個小區域原本都有有一些人，他們的離開速度是V = {2y,4x}。那麼單位時間裡人員的流動可以分解可以如下分解：

可以看到中間那個無散度場，如果16個小區域組成的大區域的邊界上，進入的人數等於離開的人數，因此是無散度的，但它可能是有旋度的，比如有一個繞圈的人正在沿著EIMNOKGFE走。按照亥姆霍茲速度分解，也就是下面這樣 關於四種速度的分解還有一些不錯的資料，比如：

http://users.metu.edu.tr/csert/me305/ME%20305%20Part%206%20Differential%20Formulation%20of%20Fluid%20Flow.pdf

http://book.ucdrs.superlib.net/views/specific/2929/bookDetail.jsp?dxNumber=000006086061&d=F62E39880A1BF4A158F629BA10D59031&fenlei=18210203 的第3.4章，文獻傳遞50頁至70頁即可。

http://book.ucdrs.superlib.net/views/specific/2929/bookDetail.jsp?dxNumber=000012787629&d=4D3FB1D67093435A9FFB7DD6A7C3926D&fenlei=13020807 《粘性流體力學》的第3.2章，大約是在89頁到99頁。而且我覺得這本書整體非常不錯。

牛頓爺爺的第二定律告訴我們，外力是導致一個物體速度改變的原因，也就是

上面F就是外力，m是物體的質量，a是加速度。但如果這個力持續作用一個時間，那麼這個物體就會從這個時間開始時的速度V0增加到時間結束後的速度V1。這裡的mv，就是動量(momentum)吧。動量的概念在高中物理上就有。而流體中動量守恆的概念為

單位時間內控制體動量的增加量通過控制面流入控制體的動量總和通過流出控制面流出控制體內的動量總和外界作用的力

其中控制體(control volume)就是空間中一個固定的地方，比如一個立方體。而控制面(control suface)就是這個立方體的六個面。流體微團可以自由從控制面流進流出。控制體和控制面實際上並不存在，只是為了方便研究而已。

我們先看看把((流入動量) - (流出動量))簡化為(-(淨流出動量))。我們假設控制面的法向量指向控制面外面。那麼流出動量的速度方向就和控制面的方向相同。

如果控制體很小，那麼速度的變化就可以忽略，因此動量的變化主要是質量的變化。而質量的變化量等於密度乘以離開控制體的流體體積。密度在這裡是常數，而離開控制體的體積等速度乘以控制面的面積。也就是

上式中，m是質量，rho是密度，V是體積，S是表面的面積。那麼最終動量的變化是

現在讓表面S細分為無窮小的區域，然後將其結果相加，也就是積分起來，就得到了淨流出動量，而(-(淨流出動量))還要在下面這個式子前加個負號。

然後是流體受到的力。比如說一個水箱破了一個洞，那麼這個洞附近的水就會受到外面空氣壓力的作用而流向外面。但壓力只能讓在流體表面的水受到影響，所以壓力也被稱為表面力(Surface Force)。我們現在需要考慮的表面力只有壓力。

由於壓力的方向由外面指向控制體裡面，與控制面的方向相反，所以這個項的前面應該加個負號。然後讓壓力在表面上積分，得到下面這個式子：

除了表面力，流體整體也會受到力的作用。比如重力，電磁力。這種力作用在所有流體微團上，不管它是不是在控制體的表面，又叫做體積力(body force)。不過，我們這裡模擬的流體會忽略掉重力的影響，而加上滑鼠拖動的影響。這時候體積力項就是：

為了將體積與速度的符號區分開，上面這個式子我用Volume來代表體積。

最後，1式等式左邊動量可以用下面這個式子表示：

單位時間內動量的變化就是

集齊2~5這四個式子，就召喚出了流體力學三大基本守恆方程之一——動量方程。

它的意思就和1式差不多，不過得稍微修改一下：

單位時間內控制體動量的增加量 = 通過流出控制面淨流出控制體內的動量總和 + 外界作用的表面力 + 外界作用的體積力

不過上面這個式子有個缺點，就是它的積分符號多得嚇人。因此我們用散度定理把面積分換成體積分，然後換一下流出動量的位置。於是式子中多了兩個梯度符號，也就是倒三角形：

然後默念咒語，讓體積分這些妖魔鬼怪快點一起離開，那麼就簡潔多了

當然也可以同除以密度，就變成了下面這種更常見的形式：

以上的動量定理的推導在各大教科書上都有，很容易查到。但是我目前覺得推導得最好的是《Fundamentals of Aerodynamics》by John D. Anderson, Jr的第2.5節，這本書本身也非常不錯。不過好像只有譯註版的。

世界萬物都遵循牛頓定律，流體也不例外。不過牛頓第一定律在流體力學中就是連續性方程，對不可壓縮流體而言，上篇我們已經推導過了，也就是速度的散度為0。牛頓第二定律則變成了動量方程。聯合起來，就得了歐拉方程:

歐拉方程與NS方程唯一的區別就是，歐拉方程沒有計算流體的粘性，而NS方程計算了流體的粘性。也就是

歐拉方程 = 重力 + 壓力

NS方程 = 重力 + 壓力 + 粘性力

但這裡我們把重力省略，而加上了來自滑鼠拖動產生的力。這裡我們不能直接把5式代入6式，因為6式裡的那項實際上是個非線性方程

直接求解歐拉方程太過困難，所以我們將其分成兩步，第一步，由上一時刻的速度，以及已知的力，在這裡就是滑鼠的拖動力，求解出中間速度(intermediate velocity)，也就是V*。

第二步，由中間速度，和壓力，求解出下一時刻的速度。

這種方法也被稱為預測校正步(predictor - corrector)。第一步並不困難，就是本系列第二篇介紹的平流方程。但這時候求解出來的中間速度，它的散度並不是零，因此第二步我們需要一點技巧將其的散度變為零，讓下一時刻的流體繼續符合不可壓縮的性質。

這篇字實在有點太多了，剩下的公式推導和代碼就放在下一篇了。

下一篇：光影帽子：【遊戲流體力學基礎及Unity代碼（五）】用歐拉方程模擬無粘性染料之代碼實現

上一篇光影帽子：【遊戲流體力學基礎及Unity代碼（三）】用波動方程模擬三維落雨池塘，連續性方程

聲明：發布此文是出於傳遞更多知識以供交流學習之目的。若有來源標註錯誤或侵犯了您的合法權益，請作者持權屬證明與我們聯繫，我們將及時更正、刪除，謝謝。

作者：光影帽子

原文：https://zhuanlan.zhihu.com/p/270530827

More：【微信公眾號】 u3dnotes