【游戏流体力学基础及Unity代码（二）】用平流方程模拟染料流动

maxueben · 发表于 2020-11-25 13:28

从现在开始，我们就以征服纳维·斯托克斯方程（Navier-Stokes equations）这个大Boss为最终目标之一。
在此之前，我们要从最基本的水系魔法开始学起。【误】
一维平流/对流

首先认识三个名词，平流，对流以及扩散，英文：Advection,Convection and Diffusion。
如果我们将一滴染料滴在清水中，那么染料会在水中扩散，也就是随机朝四面八方移动。如果清水原本也有速度，那么染料也会随着这个速度移动，这就是平流。对流与平流基本上是同义词，唯一区别是，平流是指在水中的其它物质的移动，对流是指流体本身的移动。
现在有一组一维网格，第0秒时第2到第5个网格的染料粒子的密度符合sin(i -1) + 1，其中i是所在网格的编号，其它网格染料粒子密度是零。那么其它时刻所有网格的粒子密度会怎么变化？我知道“染料粒子”这个说法既不魔法也不科学，但是为了方便先这么叫着吧。
我们现在假设这个染料粒子只能平流，速度向右，也就是某个网格下一时刻的染料粒子密度，只能由这一时刻此网格左侧网格决定，并且速度是每秒一网格，那么所有网格的水染料粒子密度将如何变化？我们可以直接这样写：
if((i-t-1) >= 1 && (i-t-1) <= 4)
height = sin(i-t-1) + 1;
else
height = 0;
但是每秒每个网格都算一次sin非常耗时，于是我们规定只有计算初始条件时可以使用sin，迭代时只能加减乘除。写出一个简洁的数学公式，也就是一维平流/对流方程

$\frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = 0 \\$
等式左边就是此网格下一时刻染料粒子密度相较于这一时刻的增长，而右边a就是速度，偏微分代表此时刻此网格与前一个网格的密度差。真实世界中，dt和dx都是无穷小的，所以上面这个公式是绝对精确的。
为了获取一维平流方程的精确离散形式，我们要用到泰勒展开。你可以在各种高等数学书上找到泰勒展开的定义。我们的一维平流方程经过泰勒展开后就是这样：

$\frac{U_{i,t+1} - U_{i,t}}{\Delta t} = a\frac{U_{i-1,t} - U_{i,t}}{\Delta x} + O(\Delta x)\\$
我知道这步跨度有点大，但其原理我们之后的篇章再讨论。上式虽然是准确无误的，但是为了算出上式的结果，我们需要算无穷多次。所以我们抛弃掉以后的二阶以及以后的项，得出下面这个式子。

$U_{i,t+1} \approx U_{i,t} - a\frac{\Delta t}{\Delta x}(U_{i,t} - U_{-1i,t}) \tag{1}$
其中a是速度，现在我们将其设置为1。dt和dx也都是1。c++代码如下，可以方便地看到各网格各时刻的值的变化。
int dx = 1, dt = 1, a = 1;//空间跨度，时间跨度，速度
const int length = 10;//10个网格
float PreValue[length], NowValue[length];

//初始条件
for (int i = 0; i < length; i++)
{
      if (i >= 1 && i <= 4)
         NowValue = sin(i - 1)  + 1;
      else
         NowValue = 0;
}

for (int t = 0; t < 100; t++)
{
      printf(&#34;第%d次迭代\n&#34;, t);
      for (int i = 0; i < length; i++)
      {
         printf(&#34;第%d个网格的值为%f\n&#34;,i,NowValue);
         PreValue = NowValue;
      }
      for (int i = 1; i < length; i++)
      {
         NowValue = PreValue - a * dt / dx * (PreValue - PreValue[i - 1]);
      }
}

如果想要可视化，可以使用python或者matlab，例如下面两份参考代码。不过这两份代码里 a * dt / dx 都小于1，所以它们最终画出的图有些奇怪，我们将在之后的篇章讨论这个。
不过就算不用可视化，我们也可以脑补出我们的染料粒子朴实无华而又枯燥的往右移动了过去，形状也没发生变化。不过这样实在太枯燥了，我们来干票大的。
Unity上模拟二维平流

二维平流方程只比一维平流多了一项：

$\frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} + b \frac{\partial u}{\partial y} = 0 \\$
其中a是x轴方向的速度，b是y轴方向的速度。离散化形式：

$\frac{U_{i,j,t+1} - U_{i,j,t}}{\Delta t} = a\frac{U_{i-1,j,t} - U_{i,j,t}}{\Delta x} + b\frac{U_{i,j-1,t} - U_{i,j,t}}{\Delta y} + O(\Delta x,\Delta y)\\$
你可以试着在c++或python上写一遍，因为之后我们在unity不会直接使用这个公式。
为了画出一张好看的二维平流图，我们首先需要画出一个黑白棋盘格纹理，代表染料的密度。新建一个Shader取名为InitDensity，修改像素着色器代码如下：
int _TexelNumber = 16;
int intuvx = floor(i.uv.x * _TexelNumber);
int intuvy = floor(i.uv.y * _TexelNumber);

if ((intuvx + intuvy) % 2)return float4<span class="p">(1.0f, 1.0f, 1.0f, 1.0f);
else  return float4(0.0f, 0.0f, 0.0f, 0.0f);
然后需要一个橙色三角形图片，尖端代表流动的方向
然后把它画到场景中上。因为我们要画出这个三角形，所以着色器不再是全屏效果，而是放在一块平板上。在摄像机上挂这么一段代码，我叫它Manager.cs。
public class Manager : MonoBehaviour
{
public GameObject prefab;
private float PlaneSize = 10.0f;
private const int QuadNumber = 16; //箭头在一个方向上的数量
private GameObject[,] Arrow = new GameObject[QuadNumber, QuadNumber];

public Material Advection;//平流
public Material InitDensity;//用于密度初始化

public RenderTexture RtPreFrame;//前一帧图像
public RenderTexture RtNowFrame;//本帧图像
void Start()
{
      Application.targetFrameRate = 10;
      float QuadSize = PlaneSize / QuadNumber;

      for (int j = 0; j < QuadNumber; j++)
      {
         for (int i = 0; i < QuadNumber; i++)
         {
            float posy = j * QuadSize - PlaneSize / 2.0f + QuadSize / 2.0f;
            float posx = i * QuadSize - PlaneSize / 2.0f + QuadSize / 2.0f;
            Arrow[i, j] = Instantiate(prefab, new Vector3(posx, 0.1f, posy), Quaternion.identity);
            Arrow[i, j].transform.localScale = new Vector3(2.0f, 2.0f, 2.0f);
            Arrow[i, j].transform.Rotate(new Vector3(90.0f, 0.0f, 0.0f));
            Arrow[i, j].transform.Rotate(new Vector3(0.0f, 0.0f, -90.0f));
         }
      }
      Graphics.Blit(null, RtPreFrame, InitDensity);//绘制初始的黑白棋盘格纹理
}
private void OnRenderImage(RenderTexture source, RenderTexture destination)
{
      Graphics.Blit(RtPreFrame, RtNowFrame, Advection);//移动之前的纹理
      Graphics.Blit(RtNowFrame, RtPreFrame);
      Graphics.Blit(source, destination);
}
}
为了防止采样误差，我们要将这些RenderTexture的尺寸设置为2的次方，比如2048 * 2048
然后新建一个用于模拟平流的着色器，我叫它Advection.shader。顶点着色器同样不用动，像素着色器代码如下：
float uvx, uvy;
uvx = i.uv.x - 1.0f / 256.0f;
uvy = i.uv.y;
if (uvx < 0)
{
uvx += 1.0f;
}
float4 col = tex2D(_MainTex, float2(uvx, uvy));
return col;

最终可以看到这样的画面，灰边是场景。
此时所有文件可在下面这次提交记录中找到。

一直往正右流太单调了，我们试着改变它的方向。首先是在cs文件中，设置默认的旋转角度为0，也就是向右，然后如果它处在右半部分并且不是在边界上，那么就逆时针旋转45度
ArrowRotation[j * QuadNumber + i] = 0.0f;
if (j != 0 && j < QuadNumber - 1 && i >= 6 && i <= 10)
{
ArrowRotation[j * QuadNumber + i] = 45.0f;
}
Arrow[i, j].transform.Rotate(new Vector3(0.0f, 0.0f, ArrowRotation[j * QuadNumber + i]));
新建一个Shader叫做UpdateVelocity.shader，由于根据箭头方向绘制速度纹理。这时候我们还需要在Manager.cs上新添加一个RenderTexture，大小依然为2048x2048，格式为R8G8B8A8_SNORM，S就是signed有符号的意思，因为我们要在这个速度图中存负一到正一的值，而一般的NORM只能存零到一的值。
int _TexelNumber = 16;
int intuvx = floor(i.uv.x * _TexelNumber);
int intuvy = floor(i.uv.y * _TexelNumber);
float rad = _ArrowRotation[intuvy * _TexelNumber + intuvx] * 3.1415927f / 180.0f;
float cosRes = cos(rad);
float sinRes = sin(rad);
return float4(cosRes, sinRes, 0.0f, 1.0f);
然后在Advection.shader文件中，获取速度图中对应位置的速度，用于改变密度。
float2 uv = i.uv - 1.0f / 256.0f * tex2D(_VelocityTex, i.uv).xy;
if (uv.x < 0)
{
uv.x += 1.0f;
}
float4 col = tex2D(_MainTex, uv);
return col;
上面这些代码已经用到了改过的平流方程。这个方程来源于Stam在1999年发表的论文“Stable Fluid”
并且2004年的《GPU Gems 3》第38章也提到了这个方法
其公式如下：
其中q代表某个物质量，在我们这里就是密度。如下图，右上角白色的点是现在我们需要求解的点，我们可以根据这一点此刻的速度算出它上一时刻在哪里，然后通过插值算出上一时刻的值，经过计算就可以得到这一点这个时刻的值。我目前看到的在unity上模拟流体平流的代码都是用这个方法。
另一个方法是，我们已经算出左下角黑点此时刻的值，为了算出下一刻的图像，我们可以计算这个黑点下一刻的位置，也就是灰点。这个灰点正好在第一行第二个网格里。然后我们“只要”算出所有下一时刻可能到达第一行第二个网格的点的值，然后平均就行了。
这种方法不仅计算量大，而且算出的数据要存起来也很麻烦。于是我们直接采用第一种方法，最终结果如下：
完整的改动可以查看git上这个提交记录，浅绿色部分是新添加的行，浅红色部分是删除的行，白色的表示无变动，比较起来很方便：
结果不怎么理想，首先由于左边界如果需要采样的点小于零，那么就会转而采样右边界的点，所以左边出现的不再是标准的棋盘格了。最方便的解决方法是再添加两个RenderTexture用于模拟棋盘格纹理稳定向右移动，然后如果对流着色器需要采样的点小于0，那么就采样这两个RenderTexture之中的一个。下面是Advection.shader的部分改动：
float2 uv = i.uv - 1.0f / 256.0f * tex2D(_VelocityTex, i.uv).xy;
float4 col = tex2D(_MainTex, uv);
if (uv.x < 0)
{
uv.x += 1.0f;
col = tex2D(_DensityLeft, uv);
}
然后新建一个UpdateVelocity.shader用于棋盘格的移动，其中像素着色器的代码如下：
float uvx = i.uv.x - 1.0f / 256.0f;
if (uvx < 0)  uvx += 1.0f;
fixed4 col = tex2D(_MainTex, float2(uvx, i.uv.y));
结果如下：
完整的改动可以查看git上这个提交记录：
另一个问题是，下图两个红圈部分，右边的明显要比左边更加模糊。
我们的RendetTexture的长度是2048，而着色器中采样的步长却是1/256.0f，而三角函数算出的值乘上1/256.0f并不总是1/2048.0f的倍数，而这会导致unity去计算插值，于是白色加灰色就变成了灰色，如此一来就模糊了。
这是因为unity的RenderTexture的采样格式默认是双线性过滤，也就是如果采样点不在像素正中间的话，就会取采样点周围四个像素加权平均起来。而根据Stam的这篇Stable Fluid论文，我们正需要这样的过滤方法。
现在我们可以试一下正弦波：
上面的结果有内味了，但是与真实的流体相比，缺点也很明显。
第一个缺点是我们计算时只有平流项，而没考虑流体的扩散和粘性。如果旁边网格的流体只要不是直接流入我这个网格，就和我一丝关系也没有。比如下图红框圈出的部分。这显然可以把牛顿爷爷气得从棺材里爬出来。
第二个是缺点是，我们的RenderTexture是2048x2048，但是用于表示流体方向的箭头以及数组却只有16x16个。
第三个缺点，也就是我们假设dt和dx，以及速度都是1，然而真实世界中不可能有这种理想的情况。不过我们暂时不管这个，毕竟可以简化很多东西。
现在我们把流体速度的方向和大小写死在代码里，能不能在运行时用鼠标点击改变这些箭头的方向，进而控制流体的方向呢？这是这篇文章留下的挑战。
然后安利一下一个9k星的webgl库，使用webgl的流体模拟，最初几次提交的代码量很少，适合阅读。不过这个库一开始就用了N-S方程。
下篇文章：（3）用波动方程模拟三维落雨池塘，连续性方程
上篇文章：（1）热传导方程

		自动登录	找回密码
密码			立即注册

[笔记] 【游戏流体力学基础及Unity代码（二）】用平流方程模拟染料流动

本帖子中包含更多资源