用歐拉方程模擬無粘性染料之代碼實現

2021-03-02 Unity3D遊戲開發精華教程乾貨
一維壓力求解

回顧上一篇第二步要求解的方程，也就是

V星是中間速度，Vn+1是下一時刻速度，rho是密度，p是壓力。V星可能有散度，而我們必須讓下一時刻速度無散度從而符合不可壓縮的性質。我們知道密度和中間速度，而壓力和下一時刻速度都是未知的。將它寫成離散化形式是

然後繼續回顧我們討論的連續方程，也就是不可壓縮流體必須要滿足散度為0的條件，不僅要在時刻n滿足，也要在時刻n+1也滿足，寫成一維形式如下：

然後令dx = 1,dt = 1,就可求得本時刻中間速度與壓力的關係

繼續變化，就能根據中間速度求得壓力

而密度需要乘的括號裡的數又是中間速度的散度。現在來看一個可以手算的簡單例子，假如求得的中間的速度u = 2x的話，那麼就希望最終的速度為 u = 0，而需要被減去的速度為u = -2x，繼而求得壓力的梯度和壓力。

我們看看一維情況下的速度分解就知道了，等式右邊那個大大的零就是我們希望求得的無散度的場的速度，也就是下一時刻的速度的微分要是0。而等式右邊第二項就是我們希望減去的那項，第二項是無旋度的，因此它可以表示為一個標量場的梯度，而壓力正好可以來充當這個標量場。

下面是一個可以手算的例子。

我們這裡避開了邊界情況，因為邊界情況介紹起來又是上千字qwq，下一篇再說如何處理邊界速度不為0的情況。

但這僅僅是一維的簡單例子，並且還事先知道了速度函數這個bug級的東西。到了二維三維，並且速度函數很複雜的時候，就沒這麼好算了。也就是，一般情況下，我們只能直接算出中間速度和中間速度的散度，而計算希望的最終的速度和壓力場就很難麻煩。

但這個解決我們在第一章就說過了，我們雖然一開始不知道確切值，但我們可以猜嘛，多猜幾次就和正確結果差不多了。所以1式右邊是壓力就是上次猜的結果，左邊的壓力就是這次的結果。

上面那個例子寫成C++代碼如下，貼到main函數裡即可運行。下面這段代碼實際上是錯的，並不能正確求得壓力，因為沒考慮邊界情況。但是看看下面的代碼至少能讓你對大體求解方法是怎樣的有個直觀的印象。

const int length = 8;double v[length], div[length], pressure[length], pressure2[length];for (int i = 0; i < length; i++){    pressure[i] = pressure2[i] = 0;    if (i >= 2 && i <= 5)        v[i] = 2 * i; else v[i] = 0;}for (int i = 0; i < length; i++){    printf("第%d網格速度wei%f\n", i, v[i]);}for (int j = 1; j < length - 1; j++){    div[j] = v[j + 1] - v[j];    printf("第%d個網格的散度是%f\n", j, div[j]);}for (int i = 0; i < 100; i++){    if (i % 10 == 0)        printf("\n第%d次迭代\n", i);    for (int j = 0; j < length; j++)    {        if (i % 10 == 0)            printf("第%d網格壓力wei%f\n", j, pressure[j]);        pressure2[j] = pressure[j];    }    for (int j = 1; j < length - 1; j++)    {        pressure[j] = 0.5 * (pressure2[j + 1] + pressure2[j - 1] - div[j]);    }}for (int i = 1; i < length - 1; i++){    double diff = pressure[i] - pressure[i - 1];    v[i] -= diff;
}for (int i = 0; i < length; i++){    printf("網格%d最終速度為%lf\n", i, v[i]);}

最後求得除了第一個網格和最後一個網格速度為0外，其它網格速度為4，雖然整個散度不為0，但現在只有邊界上有誤差，因此可以忽略掉這個誤差。以後我們使用unity模擬時將繼續忽略邊界誤差。
二維壓力求解
二維形式的散度寫成微分為
 不過上面這個我寫使用的中心差分(Central  Difference )，也就是這個網格的下一個格子減去上一個格子再除以格子長度的兩倍。而之前使用的是前向差分(Upwind Difference)，也就是這個網格減去上一個格子再除以格子長度的一倍。兩種方法各有優缺點。除此之外還有很多不同的差分法，比如Leafrog，MacCormack方法等。你可以很容易在各種計算流體力學的書上找到這些方法。
然後令dx = dy = dt = 1，求得中間速度與壓力的關係
 
改變一下各項位置，得出迭代需要的方程，同樣，等號左邊的p是本次迭代的結果，右邊的p是上次迭代的結果
 
上式密度rho所乘的括號內的式子，正好是中間速度的散度，因此將其寫為divergence，它在迭代過程中並不變化，所以在迭代之前就應該計算好防止重複計算
 
2式又被稱作壓力泊松方程(Pressure Possion Equation)，或者說是雅可比迭代。泊松方程和第一章介紹的拉普拉斯方程的區別是，前者有源，後者無源。這裡的源就是密度乘以散度。
壓力泊松方程也可從速度分解的角度考慮，下面是一個速度分解：
 無旋度場可以視為標量場的梯度，這裡的標量場就是壓力。將它們同乘以散度
 那麼無散度的散度是零，因此這項刪掉。散度與梯度的聯繫請看這篇https://zhuanlan.zhihu.com/p/136836187。上面這個式子我們就建立起了壓力和速度的關係，兩個倒三角形就是拉普拉斯運算符。寫成微分形式就是
 
繼續寫成離散形式就是，然後加上密度項，變換一下即可得到1式
 
unity上的代碼實現
首先看看是附著到攝像機上的代碼，與之前c++代碼中的步驟大體相同。
Graphics.Blit(VelocityRT, VelocityRT2);AdvectionMat.SetTexture("_VelocityTex", VelocityRT2);Graphics.Blit(VelocityRT2, VelocityRT, AdvectionMat);
SplatMat.SetTexture("_VelocityTex", VelocityRT);SplatMat.SetFloat("PointerX", (float)MouseX / Screen.width);SplatMat.SetFloat("PointerY", (float)MouseY / Screen.height);SplatMat.SetFloat("PointerDX", (float)MouseDX / Screen.width);SplatMat.SetFloat("PointerDY", (float)MouseDY / Screen.height);SplatMat.SetInt("MouseDown", MouseDown);Graphics.Blit(null, VelocityRT2, SplatMat);
DivergenceMat.SetTexture("_VelocityTex", VelocityRT2);Graphics.Blit(VelocityRT2, DivergenceRT, DivergenceMat);
PressureMat.SetTexture("_DivergenceTex", DivergenceRT);for (int i = 0; i < 20; i++){    Graphics.Blit(PressureRT, PressureRT2);    PressureMat.SetTexture("_PressureTex", PressureRT2);    Graphics.Blit(DivergenceRT, PressureRT, PressureMat);}
SubtractMat.SetTexture("_VelocityTex", VelocityRT2);SubtractMat.SetTexture("_PressureTex", PressureRT);Graphics.Blit(VelocityRT2, VelocityRT, SubtractMat);
Graphics.Blit(DyeRT, DyeRT2);DisplayMat.SetTexture("_DyeTex", DyeRT2);DisplayMat.SetTexture("_VelocityTex", VelocityRT);Graphics.Blit(DyeRT2, DyeRT, DisplayMat);
 Graphics.Blit(DyeRT, destination);

然後分別是各步的著色器代碼。第一步平流的，和本系列第二節的基本上相同，最後乘一個消散比例，讓速度不至於一直存在。
float Speed = 2.0f;float uv = i.uv - Speed * _VelocityTex_TexelSize.x * tex2D(_VelocityTex, i.uv).xy;float disspation = 0.999f;float4 col = disspation * tex2D(_VelocityTex, uv);

第二步，也就是滑鼠拖動造成的力，這裡的函數我使用的是這個webgl庫https://github.com/PavelDoGreat/WebGL-Fluid-Simulation所使用的，效果還很不錯~
float2 pointeruv = float2(PointerX, PointerY);float2 p = i.uv - pointeruv;float radius = 0.001;//圓的半徑float3 color = float3(PointerDX, PointerDY, 0.0f) * 50.0f;//根據滑鼠的方向產生不同方向的速度float3 splat = pow(2.1, -dot(p, p) / radius) * color;float3 base = tex2D(_VelocityTex, i.uv).xyz;if (MouseDown == 1)base += splat;float4 col = float4(base, 1.0f);return col;

第三步，計算散度的，方法和計算高斯模糊差不多，使用的公式為3式
float Top = tex2D(_VelocityTex, i.uv + float2(0.0f, _VelocityTex_TexelSize.y)).y;float Bottom = tex2D(_VelocityTex, i.uv + float2(0.0f, -_VelocityTex_TexelSize.y)).y;float Right = tex2D(_VelocityTex, i.uv + float2(_VelocityTex_TexelSize.x, 0.0f)).x;float Left = tex2D(_VelocityTex, i.uv + float2(-_VelocityTex_TexelSize.x, 0.0f)).x;float divergence = 0.5f * (Right - Left + Top - Bottom);return float4(divergence, 0.0f, 0.0f, 0.0f);

第四步：迭代計算壓力的，方法也和高斯模糊差不多，使用的公式為2式
float Top = tex2D(_PressureTex, i.uv + float2(0.0f, _PressureTex_TexelSize.y)).y;float Bottom = tex2D(_PressureTex, i.uv + float2(0.0f, -_PressureTex_TexelSize.y)).y;float Right = tex2D(_PressureTex, i.uv + float2(_PressureTex_TexelSize.x, 0.0f)).x;float Left = tex2D(_PressureTex, i.uv + float2(-_PressureTex_TexelSize.x, 0.0f)).x;float div = tex2D(_DivergenceTex, i.uv).x;float alpha = 0.5f;float pressure = (Top + Bottom + Right + Left - alpha * div) * 0.25f;return float4(pressure, 0.0f, 0.0f, 1.0f);

第五步，減去壓力梯度的，得到最終無散度，符合不可壓縮性質的下時刻速度，也很簡單
float Top = tex2D(_PressureTex, i.uv + float2(0.0f, _PressureTex_TexelSize.y)).y;float Bottom = tex2D(_PressureTex, i.uv + float2(0.0f, -_PressureTex_TexelSize.y)).y;float Right = tex2D(_PressureTex, i.uv + float2(_PressureTex_TexelSize.x, 0.0f)).x;float Left = tex2D(_PressureTex, i.uv + float2(-_PressureTex_TexelSize.x, 0.0f)).x;float2 velocity = tex2D(_VelocityTex, i.uv).xy;float factor = 0.5f;velocity.xy -= factor * float2(Right - Left, Top - Bottom);return float4(velocity, 0.0f, 1.0f);

第六步：用最終速度平流顏色
float2 coord = i.uv - _VelocityTex_TexelSize.x * tex2D(_VelocityTex, i.uv).xy;float4 col = tex2D(_DyeTex, coord);return col;

最後特別需要注意的是RenderTexture的格式，除了Dye以外，其它的都是R16G16B16A16_SNORM，如果不用16位那麼精度不夠，會出現各種奇怪的結果。SNORM是為了讓能RenderTexture存儲負數。
最後的坑爹的導入圖片問題，要允許讀寫，並且更改到合適的格式
完整的代碼在這個倉庫的第五次提交：
https://github.com/clatterrr/FluidSimChineseTutorialsWithUnity/tree/b66b5b15ca4ff46741c53502a0da85c33f7474ae
然後我們就可以破壞世界，扭曲現實了！【誤】
速度圖如下：
可視化
我們可以分析一下剛才提到的那個webgl庫，修改最終的輸出圖像代碼如下
gl.viewport(0, 0, gl.drawingBufferWidth, gl.drawingBufferHeight);displayProgram.bind();gl.uniform1i(displayProgram.uniforms.uTexture, density.first[2]);blit(null);


中間速度可視化如下：
中間速度的散度可視化如下，滑鼠從左往右，此時在那個大紅亮圓的右側
壓力可視化如下，滑鼠從左往右，此時在那個大黑圓的右側
最後
如果你搞定了本系列第四，五篇，特別是速度分解和壓力泊松方程，那麼大部分不可壓縮流體模擬代碼的主幹部分你就已經懂了，畢竟NS方程只是在歐拉方程上加個粘性項。
最後，介紹流體模擬的資料真是很稀少，因此我將本系列第四五章參考的主要資料放上來
https://developer.nvidia.com/sites/all/modules/custom/gpugems/books/GPUGems/gpugems_ch38.html首先是大名鼎鼎的GPU GEMS3，當然GPU GEM3把所有能省略的數學公式推導都省略了，直接看至少我沒看懂...不過GPU GEM3也是參考1999年的一篇論文《Stable Fluid》的http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.15.9203&rank=1，不過不建議直接看這篇論文
《Fluid Simulation for Computer Graphics》, Second Edition by Bridson, Robert ，也是很有名的一本書，前半本書基本上都在介紹模擬NS流體。這本書是2015年寫的，但這本書的前身似乎是一次notehttps://www.cs.ubc.ca/~rbridson/fluidsimulation/fluids_notes.pdf。這本書總體也很不錯，值得一讀https://b-ok.global/book/2871841/46dc7e
《Computational Fluid Dynamics Incompressible Turbulent Flows》by Takeo Kajishima Kunihiko Taira的第三章，這本書總體上也很不錯
《Essential Computational Fluid Dynamics》 by Oleg Zikanov 第10章。
《Guide To CFD》一篇很短的文章https://www.montana.edu/mowkes/research/source-codes/GuideToCFD.pdf
https://nbviewer.jupyter.org/github/barbagroup/CFDPython/blob/master/lessons/14_Step_11.ipynb就是那個NavierStokes的第11章
http://www.thevisualroom.com/02_barba_projects/barba_cfd_projects.html用python解決CFD問題的例子合集
http://folk.ntnu.no/leifh/teaching/tkt4140/._main000.html在線的有關數值分析的資料
http://jamie-wong.com/2016/08/05/webgl-fluid-simulation/，這是一個個人博客，使用Webgl模擬的，並做了簡單介紹。也省略了很多東西。這網頁經常打不開，因此我把它緩存成一個pdf了，當然圖片和文字位置有點錯亂。https://wwe.lanzous.com/iy4fthw8cva
https://wwe.lanzous.com/iy4fthw8cva，也是一個個人博客，比上面那個稍微詳細一點。同樣打開速度巨慢，因此緩存版本的PDF下載https://wwe.lanzous.com/idl0Vhw8cwb
然後是幾個github庫
https://github.com/PavelDoGreat/WebGL-Fluid-Simulation
https://github.com/kodai100/Unity_EulerianFluidSimulation，在CPU上模擬的
https://github.com/candycat1992/2DFluidSim
https://github.com/Scrawk/GPU-GEMS-2D-Fluid-Simulation
https://github.com/keijiro/StableFluids這個使用了computeShader
上一篇：光影帽子：【遊戲流體力學基礎及Unity代碼（四）】用歐拉方程模擬無粘性染料之公式推導
下一篇：光影帽子：【遊戲流體力學基礎及Unity代碼（六）】用NavierStokes方程模擬粘性染料流動

聲明：發布此文是出於傳遞更多知識以供交流學習之目的。若有來源標註錯誤或侵犯了您的合法權益，請作者持權屬證明與我們聯繫，我們將及時更正、刪除，謝謝。
作者：光影帽子
原文：https://zhuanlan.zhihu.com/p/270531017
More：【微信公眾號】 u3dnotes

相關焦點

  • 【遊戲流體力學基礎及Unity代碼(四)】用歐拉方程模擬無粘性染料之公式推導
    所以下面的內容是完全適合NS方程的。各位請準備好！模擬流體的時候會遇到許多數學公式。為了深刻理解這些數學公式，我們先從最簡單的複習起，比如說二十以內的加減法。 歐拉方程與NS方程唯一的區別就是，歐拉方程沒有計算流體的粘性，而NS方程計算了流體的粘性。
  • 【遊戲流體力學基礎及Unity代碼(二)】用平流方程模擬染料流動
    另一個方法是，我們已經算出左下角黑點此時刻的值，為了算出下一刻的圖像，我們可以計算這個黑點下一刻的位置，也就是灰點。這個灰點正好在第一行第二個網格裡。第一個缺點是我們計算時只有平流項，而沒考慮流體的擴散和粘性。如果旁邊網格的流體只要不是直接流入我這個網格，就和我一絲關係也沒有。比如下圖紅框圈出的部分。這顯然可以把牛頓爺爺氣得從棺材裡爬出來。
  • 科學網—帶自由邊界歐拉方程的幾何分析與先驗估計
    帶自由邊界歐拉方程的幾何分析與先驗估計
  • 流體中失效的歐拉方程
    歐拉方程是一種理想化的對流體運動的數學描述，它們在一定的假設範圍內，模擬流體的運動。更確切地說，歐拉方程描述了流體中無窮小的粒子的瞬時運動。這個描述包括一個粒子的速度和它的渦量（即旋轉的速度和方向）。顯然，真正的流體內部是有摩擦的，歐拉方程的描述將流體的運動規定在了一個特定的理想化世界中。若要模擬更真實的流體運動，則需要使用納維-斯託克斯方程（NS方程）。
  • 數學上最複雜的公式——納維斯託克斯方程,到底困難在哪兒?
    1775年，歐拉大神決定換個口味，去研究了一個與力學相關的數學領域——流體力學。他從最基本的無粘性流體的特性開始，仔細研究了無粘性流體的運動與動量變化的關係，於是寫成一本《流體運動的一般原理》。書裡留下了一個無粘性流體力學領域最重要的基礎方程。這本書也是流體力學的開山之作。
  • 北京電影學院發了滿是數學公式的計算機頂會論文,並開源了其代碼
    而諸如洪水、煙霧、爆炸等特效計算的背後，實際上是用電腦程式在求解已有百年歷史的「納維-斯託克斯方程」這個方程，對於做流體動力學的讀者一定不陌生，數十年來科學家們為了計算機翼升力，已將其研究了百千萬遍。然而基於影視製作的特別需求，影視科技工作者們對這個方程的求解提出了新的需求。我們需要能夠處理更大的時間步長以及不損失精度細節！！
  • 基於部分浸潤效應的歐拉壁膜流動形態演變模擬仿真
    ■ 質量守恆和動量守恆：歐拉壁膜模型的核心基礎■ 能量守恆：涉及壁膜的傳熱需要啟用■ 由於歐拉壁膜模型適用的液膜很薄，故採用潤滑近似(平行流動)，因此這些控制方程是在平行於壁面的局部當地坐標中求解的。■ 在無外部流體流動的重力場下，如壁膜生長流動問題中，壁膜速度完全受重力驅動。
  • 尋找「最好」(2)歐拉-拉格朗日方程
    歐拉-拉格朗日方程(Euler -Lagrange equation) 為變分法中的一條重要方程。
  • 《歐拉方程及微分方程建模》思路與方法
    一、歐拉方程及其求解方法具有結構的變係數線性微分方程稱之為歐拉方程.令x=eu，則u=lnx，於是有將原歐拉方程中xky(k)全部用上式代入，則可以將原方程轉化為以y為函數，u為自變量的常係數線性微分方程
  • 變分法——歐拉-拉格朗日方程
    然後根據"變分法基本引理"(參見變分法基本引理)就可以導出歐拉-拉格朗日方程啦:需要注意的是,歐拉方程是泛函極值的必要條件,但不是充分條件,>在處理實際泛函極值問題時,一般不去考慮充分條件,而是從實際問題的性質出發,間接地判斷泛函極值的存在性,直接使用歐拉方程求出極值曲線【往期精選】●1+1/2+...+1/n  ,可能是整數嗎?
  • 39.積分、泛函 + 歐拉-拉格朗日方程、實數、標量、變分法、極值、弧微分、範數(數學篇)
    2.2.1 歐拉-拉格朗日方程這個方程是泛函中非常重要的方程，也是非常經典的能量泛函極小化的方法，不論在物理還是計算機領域，應用非常廣泛。所謂能量泛函，是指微分的範數平方再積分。這個方程是泛函中非常重要的方程，也是非常經典的能量泛函極小化的方法，不論在物理還是計算機領域，應用非常廣泛。
  • 數學第一家族和「伯努利方程」
    況且伯努利家族一向以擅長微積分聞名於世，歐拉是一個合格的助理，他自作主張把伯努利方程升級為伯努利方程2.0：再對比一下伯努利方程，有沒有感覺到一下子從中學到了碩士的水平。此時的歐拉剛剛20歲，不愧是著名數學家兼教育家約翰.伯努利的弟子，他把動量和質量守恆植入到流體力學之中，開創了無黏流體的統一方程，雖然後世流體力學研究者萬千，優秀的數學家也不可勝數，但兩大守恆無人敢背叛……不過，丹尼爾自己是怎麼想的，我們就不得而知了，不過作為歐拉一生的摯友
  • 世界級千禧難題「納維-斯託克斯方程」:數學史上最複雜的公式
    這個方程並不是一個人提出來的，1775年，著名數學家歐拉，對，沒有錯就是數學界四大天王歐拉，他如今又來摻和流體力學了，他在《流體運動的一般原理》一書中根據無粘性流體運動時流體所受的力和動量變化從而推導出了一組方程。
  • 【遊戲流體力學基礎及Unity代碼(一)】熱傳導方程
    有些文章用尖峰波或者FFT模擬，但那畢竟是統計學方法，和流體力學還是不搭邊。其餘的文章倒是用了納韋斯託克方程，但那也僅僅是把納韋斯託克方程寫了一遍，好一點的還給仍你一些居複雜的代碼，對我等入門者來說實在是看不懂。
  • 困擾人類200年,數學史最難最複雜的公式之一:納維-斯託克斯方程
    這個方程並不是一個人提出來的，1775年，著名數學家歐拉，對，沒有錯就是數學界四大天王歐拉，他如今又來摻和流體力學了，他在《流體運動的一般原理》一書中根據無粘性流體運動時流體所受的力和動量變化從而推導出了一組方程。
  • 歐拉是如何用根式解的形式表示方程sinx=0的
    這是一個含有立方的多項式，如果它在-1,0,1三個點的位置等於0那麼這個立方多項式就可以寫成如下形式：每個0點位置都有一個乘積因式通過在垂直方向非零常量進行縮放，可以輕鬆解決問題如果您知道一個有n個0點的n次方程，該又該如何確定呢？
  • 湍流的描述方程LES-NS(大渦模擬) 和RANS (雷諾平均) 的區別
    1　　基礎　　用幾個詞來概括湍流的本質：三維，非定常（隨時間變化），多尺度。這就導致了直接模擬湍流計算代價非常大。為了在有限的計算機資源下模擬湍流，各種前輩大牛提出了幾種方法，包括了LES 和RANS。　　LES，中文名大渦模擬，基本思想是對NS 方程進行某種過濾，然後只計算大尺度的湍流，而將小於過濾尺度的湍流用模型加以刻畫。數學上，小於過濾尺度的湍流表現為額外的應力項，稱為亞網格應力。現有的湍流理論已經有結論，幾乎所有的湍流在足夠小的尺度上都具有一定的相似性。
  • 此「歐拉」非彼「歐拉」 你可真正了解歐拉?
    人生最後7年，歐拉的雙目完全失明，他還是以驚人的速度產出了生平一半的著作。歐拉著作的驚人多產並不是偶然的，他可以在任何不良的環境中工作，他常常抱著孩子在膝上完成論文，也不顧孩子在旁邊喧譁。他那頑強的毅力和孜孜不倦的治學精神，歐拉的記憶力和心算能力是罕見的，他能夠複述年青時代筆記的內容，心算並不限於簡單的運算，高等數學一樣可以用心算去完成。
  • 建築工程概述之水力學
    二十世紀初的重要突破是普朗特的邊界層理論，它把無粘性理論和粘性理論在邊界層概念的基礎上聯繫起來。　　二十世紀蓬勃發展的經濟建設提出了越來越複雜的水力學問題：高濃度泥沙河流的治理；高水頭水力發電的開發；輸油幹管的敷設；採油平臺的建造；河流湖泊海港汙染的防治等。
  • 看得懂的數學之美:從青年歐拉對巴塞爾問題的解法說起
    機器之心報導參與：思、肖清什麼是數學之美？就是思考的時候忘記時間的流逝，解答或證明後無與倫比的快樂。歐拉，歷史上最重要的數學家之一，也是最高產的數學家，平均每年能寫八百多頁論文。我們經常能見到以他名字命名的公式與定理，可能最廣為人知的便是「世界上最美的公式」歐拉公式。