【游戏流体力学基础及Unity代码（一）】热传导方程

dice4dj

在游戏中模拟流体并不是什么新鲜事，但是我几乎就没看到什么好的入门文章。有些文章用尖峰波或者FFT模拟，但那毕竟是统计学方法，和流体力学还是不搭边。其余的文章倒是用了纳韦斯托克方程，但那也仅仅是把纳韦斯托克方程写了一遍，好一点的还给仍你一些居复杂的代码，对我等入门者来说实在是看不懂。
最近我在github上找到三个不错的python jupyter notebook库
1.首先是大名鼎鼎的12步纳韦斯托克方程，oschina上有中文翻译，油管上python版本和matlab版本的对应编程视频
barbagroup/CFDPython
2.然后是格子波兹曼方法，也是除了纳韦斯托克方法外的一个经典方方法
pylbm/pylbm
3.然后这个空气动力学python库，虽然说是空气动力学，但和流体力学相同的地方很多
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
{
    // Start is called before the first frame update

    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度，最后得到这样的图。
现在我们已经会搓小火球了，接下来我们正式开始学习如何搓水球，下一篇文章《用平流方程模拟染料流动》：