【遊戲流體力學基礎及Unity代碼(一)】熱傳導方程

在遊戲中模擬流體並不是什麼新鮮事，但是我幾乎就沒看到什麼好的入門文章。有些文章用尖峰波或者FFT模擬，但那畢竟是統計學方法，和流體力學還是不搭邊。其餘的文章倒是用了納韋斯託克方程，但那也僅僅是把納韋斯託克方程寫了一遍，好一點的還給仍你一些居複雜的代碼，對我等入門者來說實在是看不懂。

最近我在github上找到三個不錯的python jupyter notebook庫

1.首先是大名鼎鼎的12步納韋斯託克方程，oschina上有中文翻譯，油管上python版本和matlab版本的對應編程視頻

https://github.com/barbagroup/CFDPython

2.然後是格子波茲曼方法，也是除了納韋斯託克方法外的一個經典方方法

https://github.com/pylbm/pylbm

3.然後這個空氣動力學python庫，雖然說是空氣動力學，但和流體力學相同的地方很多

https://github.com/barbagroup/AeroPython

於是這個「遊戲流體力學基礎」系列將在這些庫的幫助下，從最基本的方程開始推導，儘量保證貓都都能看懂，來完成遊戲中的流體模擬吧！預計會有十幾篇。

好了，我們開始學習一些基本的火系魔法知識吧！【誤】

首先我們看一個這樣的問題，對於一個一維物體，如果它的左端點溫度是700度，右端點溫度是300度，並且這個物體上的溫度是連續均勻變化的，那麼這個物體的其它位置的溫度是多少？

當然，由於計算機只能求解離散的物體，我們需要將問題轉化成如下形式，將一個物體分成八格，每格都有一個自己的溫度值。如果最左邊一格的溫度是700度，最右邊一格的溫度是300度，並且每格之間的溫度變化也是均勻連續的，那麼中間六格的溫度各是多少？

你會說，這還不簡單嘛，不就是解一元一次方程嘛。但是這僅僅是一維的情況，並且我們確定的溫度值只有兩個而已。如果到了三維，並且確定的溫度值有很多個的時候，要求解的溫度值也很多的時候，要通過解析方法求解就很麻煩了。我們得換個法子。

既然我們不知道中間六格的溫度，我們就瞎猜嘛，說不定猜對了呢？

我們猜中間六格的溫度都是0度。嗯，結果顯然看起來不太對，但我們可以改一下。除了我第一格和第八格的溫度已經給定不能再改變外，其它六格的溫度我們都可以繼續改，讓結果看起來更加合理一些。

牛頓爺爺說過，當兩個互相接近的物體溫度不相等的時候，它們的溫度肯定會被互相影響。也就說，第二格的溫度會因第一格和第三格的溫度的影響而改變，第三格的溫度會因為第二格和第四格的溫度的影響而改變....但具體會如何被影響呢？

某個方格的值會被附近方格的值影響，是不是很熟悉？對啦，高斯模糊就是用的這種方法！

簡單起見，我們就設第二格的溫度等於（第一格的溫度 加上 第三格的溫度）除以二。於是經過一次計算後，每個方格的溫度就變成了下面的樣子。

顯然結果看起來還是不太對。但是，熱量會持續傳導，時間會見證一切。

可以看到在第100秒，方格之間的溫度差是57度左右，拿出除法一算，400 / 7 大約也是57度，很不錯的。實際上第50秒的結果已經很接近第100秒的結果了，所以如果對精度要求不那麼高，那麼就算100次就結束吧。

既然我們的瞎猜法很不錯，那麼我們就要總結經驗，把瞎猜法發揚光大。不過我們之前說的：「第二格的溫度等於（第一格的溫度 加上 第三格的溫度）除以二」 這句話太長了，我們將其轉換為數學語言：

其中T就是溫度，下標n代表第n個方格，t代表第t秒。然後我們就可以得到第n個方格第t+1秒的溫度和第t秒溫度之間的關係：

不過上面這個公式還不夠簡潔，讓我們將其寫為人看不懂的形式吧：

dt即時間的跨度，我們之前將其設置為1秒，也就是1。dx就是空間的跨度，也是不同網格之間的距離，我們之前將其設置為1格，也就是1。這裡，時間t只有一階導數是因為時刻t的溫度只受到t-1的溫度的影響。而x為二階導數是因為，第n個格子會受到n-1和n+1兩個格子的影響。

然後我們將這個公式寫得更加簡潔一些，然後取個名字，就叫一維熱傳導方程吧：

這裡的a就是一個迭代參數，我們之前將其設置為1/2。實際上可以設置得更小，也就是熱傳導速度更慢，防止出現一些莫名奇妙的結果。當然這樣求解起來也就更慢了。如果大於1/2，那麼就會出現不穩定的狀態，比如將其設置為0.7，僅僅第三次迭代後第二個網格的溫度就高於730度了，產生不穩定狀態的原因也很好理解，第三格實際上不能把0.7的熱量同時分配給第二格和第四格，第三格沒那麼高的溫度。

好了，現在你可以拿這個公式去欺負那些沒見過世面的麻瓜了！記住把式子改成下面的樣子讓他們更加看不懂，其中三角形是拉普拉斯算子，u是待求解的量，我們這裡就是溫度。

不過這個的一維熱傳導方程是否還能寫成別的形式呢？

比如我們要測量b處的溫度變化速度，而a處的溫度會影響b處，它們之間的關係應該是怎樣的呢？首先它們之間會有一個介質，介質面積越大，熱傳導速度越快，熱量變化速度也越快，這就是面積A。介質厚度越大，熱量變化越慢，這就是厚度dx。a處和b處的溫度差越大，溫度變化也越快，這就是溫度差Ta - Tb。除了這個三個因素外，介質的材質參數k也會影響溫度變化速率，如下圖所示。

根據傅立葉法則，我們得到熱傳導方程為如下形式

其中等式左邊是熱量的變化速度。所以對於之前的迭代參數a來說，它的實際物理意義是材質參數乘上介質面積。

先廢話兩句。其實這種模擬流體的代碼通常寫在Python和Matlab和Foratan上，上github上一搜一大把，當然看不看得懂就是另一回事了。但是由於看這篇文章的讀者對unity更加熟悉一些，而且很方便移植到opengl/directx上，所以就用unity寫了。順便吐槽一句unreal寫著色器太麻煩了。

看懂下面的內容需要一點unity著色器的知識，如果你沒有，現在立馬去把那本《UnityShader入門精要》吃掉。

我們在著色器上實現的方法現在很簡單，就是在攝像機上掛一個全屏處理腳本，直接把著色器生成的材質放到整塊屏幕上。

我們需要三個RenderTexture，第一個_rtNull是全黑的紋理，也就是初始狀態溫度全為0，此時還沒有設置邊界條件。

第二個_rt0代表上一幀的結果，經過計算後得到本幀的結果_rt1，再將後者輸出的屏幕上。最後再將rt1複製給rt0，用於下一幀的計算。

下面是要掛到攝像機上的代碼，之後我們會將_TexelNumber設置為8，也就是方格的數量。同時將每秒的幀數設置為2，方便我們動態觀察。

public class FullScreen : MonoBehaviour{    
    public Material _mat;    public RenderTexture _rtNull;    public RenderTexture _rt0;    public RenderTexture _rt1;    [Range(2,100)]    public int _TexelNumber;
    private void Start()    {        Graphics.Blit(_rtNull, _rt0);        Application.targetFrameRate = 2;        _TexelNumber = 8;    }    private void Update()    {        _mat.SetInt("_TexelNumber", _TexelNumber);        _mat.SetTexture("_MainTex", _rt0);    }
    private void OnRenderImage(RenderTexture source, RenderTexture destination)    {        Graphics.Blit(_rt0, _rt1,_mat);        Graphics.Blit(_rt1, destination);        Graphics.Blit(_rt1, _rt0);    }}

然後在資源管理器中創建這三個RenderTexture，然後拖到腳本上即可。注意RenderTexture大小要改成和輸出屏幕一樣大小。
我們現在不用修改頂點著色器，所以下面是像素著色器的代碼，原理和高斯模糊很類似，採集上一幀圖像中左邊網格的值和右邊網格的值再平均，就得到新的溫度值。
float  _Temprature  = 0.0f;float  _TexelSize  = 1.0f / _TexelNumber;if ( i.uv.x < _TexelSize )  _Temprature = 0.7f;else if ( i.uv.x >= 1 - _TexelSize )  _Temprature = 0.3f;else{  float2  leftuv  = float2( i.uv.x - _TexelSize, i.uv.y );  float  leftT  = tex2D( _MainTex, leftuv ).g;  float2  rightuv = float2( i.uv.x + _TexelSize, i.uv.y );  float  rightT  = tex2D( _MainTex, rightuv ).g;  _Temprature = (leftT + rightT) / 2.0f;}return(float4( 1.0f, _Temprature, 0.0f, 1.0f ) );
最後Return語句的原理如下。我們假設我們模擬的世界中，溫度的顏色如下變化
於是，我們就可以將RGB中的R設置為1，B設置為0，僅將G作為變量。當G從0到1變化時，我們就假設溫度從0度到1000度變化，於是顏色就也如上圖變化。點擊，播放，現在場景看起來如下，大約過了十秒左右顏色變化就不太明顯了。如果把_TexelNumber改得更大一些，那麼迭代需要的時間也更久。
是不是很簡單？現在你已經初步掌握了火系魔法的要點，快用爆裂魔法去和魔王軍首領對決吧！相信你一定可以的！【誤】


之前說過我們的瞎猜法，在一維的情況下，就是左邊格子和右邊格子的平均值。那麼在二維情況下怎麼辦？你肯定已經想到了，就是上下左右四個格子的平均值！三維就是上下左右前後六個格子的平均值，就這麼簡單。
但不僅僅能取上下左右四個格子，對角線上四個其實也可以取，就連要求解的格子本身也可以取上一幀的值，不過要乘上一些不同的權重。這就是格子波茲曼方法，lattice boltzmann method。這個我們之後的篇章再討論。
寫成離散化形式就是如下：
二維代碼和三維方程我就不貼了。就把它當作一個挑戰留下來吧。值得注意的是我們這裡的dt，dx和dy都是1，而a也僅僅是簡單的權重，所以算是大大簡化了這個方程。三維方程代碼編寫方式更加複雜一些，我們之後的篇章再討論。
我將左下角溫度設置為700度，右上角300度，右下角0度，最後得到這樣的圖。
現在我們已經會搓小火球了，接下來我們正式開始學習如何搓水球，下一篇文章《用平流方程模擬染料流動》：
宣傳：我創建了一個模擬流體交流Q群，歡迎各位喜歡自己編寫代碼模擬流體的小夥伴們入群互相交流學習！群號：1001290801。

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

More：【微信公眾號】 u3dnotes