论坛 › 代码库 Code › 【游戏流体力学基础及Unity代码（五）】用欧拉方程模拟无粘性染料之代码实现

fly_027@163.com 发表于 2020-11-29 11:06

【游戏流体力学基础及Unity代码（五）】用欧拉方程模拟无粘性染料之代码实现

一维压力求解

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

https://www.zhihu.com/equation?tex=%5Cfrac%7BV%5E%7Bn%2B1%7D+-V%5E%7B%2A%7D%7D%7B%5CDelta+t%7D+%3D%5Cfrac%7B1%7D%7B%5Crho%7D%5Cnabla+p+%5C%5C V星是中间速度，Vn+1是下一时刻速度，rho是密度，p是压力。V星可能有散度，而我们必须让下一时刻速度无散度从而符合不可压缩的性质。我们知道密度和中间速度，而压力和下一时刻速度都是未知的。将它写成离散化形式是

https://www.zhihu.com/equation?tex=u_%7Bi%7D%5E%7Bn%2B1%7D+%3D+u_%7Bi%7D+-+%5CDelta+t+%5Cfrac%7B1%7D%7B%5Crho%7D%5Cfrac%7Bp_%7Bi%2B1%7D+-+p_%7Bi%7D%7D%7B%5CDelta+x%7D+%5C%5C 然后继续回顾我们讨论的连续方程，也就是不可压缩流体必须要满足散度为0的条件，不仅要在时刻n满足，也要在时刻n+1也满足，写成一维形式如下：

https://www.zhihu.com/equation?tex=%5Cnabla+%5Ccdot+V+%3D+%5Cfrac%7Bu_%7Bi%7D%5E%7Bn%2B1%7D+-+u_%7Bi-1%7D%5E%7Bn%2B1%7D%7D%7B%5CDelta+x%7D+%3D+0+%5C%5C 然后令dx = 1,dt = 1,就可求得本时刻中间速度与压力的关系

https://www.zhihu.com/equation?tex=%28u_%7Bi%7D+-+%5Cfrac%7B1%7D%7B%5Crho%7D%28p_%7Bi%2B1%7D+-+p_%7Bi%7D%29%29+-+%28u_%7Bi-1%7D+-+%5Cfrac%7B1%7D%7B%5Crho%7D%28p_%7Bi%7D+-+p_%7Bi-1%7D%29%29+%3D+0+%5C%5C 继续变化，就能根据中间速度求得压力

https://www.zhihu.com/equation?tex=p_%7Bi%7D+%3D+%5Cfrac%7B1%7D%7B2%7D%28p_%7Bi%2B1%7D+%2B+p_%7Bi-1%7D+-+%5Crho%28u_%7Bi%7D+-+u_%7Bi-1%7D%29%29+%5Ctag%7B1%7D
而密度需要乘的括号里的数又是中间速度的散度。现在来看一个可以手算的简单例子，假如求得的中间的速度u = 2x的话，那么就希望最终的速度为 u = 0，而需要被减去的速度为u = -2x，继而求得压力的梯度和压力。
我们看看一维情况下的速度分解就知道了，等式右边那个大大的零就是我们希望求得的无散度的场的速度，也就是下一时刻的速度的微分要是0。而等式右边第二项就是我们希望减去的那项，第二项是无旋度的，因此它可以表示为一个标量场的梯度，而压力正好可以来充当这个标量场。

https://www.zhihu.com/equation?tex=%5Cbegin%7Bpmatrix%7D+%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+x%7D+%5Cend%7Bpmatrix%7D+%3D++%5Cbegin%7Bpmatrix%7D+0+%5Cend%7Bpmatrix%7D+%2B+%5Cbegin%7Bpmatrix%7D+%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+x%7D+%5Cend%7Bpmatrix%7D++%5C%5C
下面是一个可以手算的例子。
我们这里避开了边界情况，因为边界情况介绍起来又是上千字qwq，下一篇再说如何处理边界速度不为0的情况。
但这仅仅是一维的简单例子，并且还事先知道了速度函数这个bug级的东西。到了二维三维，并且速度函数很复杂的时候，就没这么好算了。也就是，一般情况下，我们只能直接算出中间速度和中间速度的散度，而计算希望的最终的速度和压力场就很难麻烦。
但这个解决我们在第一章就说过了，我们虽然一开始不知道确切值，但我们可以猜嘛，多猜几次就和正确结果差不多了。所以1式右边是压力就是上次猜的结果，左边的压力就是这次的结果。
上面那个例子写成C++代码如下，贴到main函数里即可运行。下面这段代码实际上是错的，并不能正确求得压力，因为没考虑边界情况。但是看看下面的代码至少能让你对大体求解方法是怎样的有个直观的印象。
const int length = 8;
double v, div, pressure, pressure2;
for (int i = 0; i < length; i++)
{
    pressure = pressure2 = 0;
    if (i >= 2 && i <= 5)
      v = 2 * i; else v = 0;
}
for (int i = 0; i < length; i++)
{
    printf("第%d网格速度wei%f\n", i, v);
}
for (int j = 1; j < length - 1; j++)
{
    div = v - v;
    printf("第%d个网格的散度是%f\n", j, div);
}
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);
      pressure2 = pressure;
    }
    for (int j = 1; j < length - 1; j++)
    {
      pressure = 0.5 * (pressure2 + pressure2 - div);
    }
}
for (int i = 1; i < length - 1; i++)
{
    double diff = pressure - pressure;
    v -= diff;

}
for (int i = 0; i < length; i++)
{
    printf("网格%d最终速度为%lf\n", i, v);
}最后求得除了第一个网格和最后一个网格速度为0外，其它网格速度为4，虽然整个散度不为0，但现在只有边界上有误差，因此可以忽略掉这个误差。以后我们使用unity模拟时将继续忽略边界误差。
二维压力求解

二维形式的散度写成微分为

https://www.zhihu.com/equation?tex=%5Cnabla+%5Ccdot+V+%3D+%5Cfrac%7Bu_%7Bi%2B1%2Cj%7D%5E%7Bn%2B1%7D+-+u_%7Bi-1%2Cj%7D%5E%7Bn%2B1%7D%7D%7B2%5CDelta+x%7D+%2B+%5Cfrac%7Bv_%7Bi%2Cj%2B1%7D%5E%7Bn%2B1%7D+-+v_%7Bi%2Cj-1%7D%5E%7Bn%2B1%7D%7D%7B2%5CDelta+y%7D%3D+0+%5C%5C 不过上面这个我写使用的中心差分(CentralDifference )，也就是这个网格的下一个格子减去上一个格子再除以格子长度的两倍。而之前使用的是前向差分(Upwind Difference)，也就是这个网格减去上一个格子再除以格子长度的一倍。两种方法各有优缺点。除此之外还有很多不同的差分法，比如Leafrog，MacCormack方法等。你可以很容易在各种计算流体力学的书上找到这些方法。
然后令dx = dy = dt = 1，求得中间速度与压力的关系

https://www.zhihu.com/equation?tex=%28u_%7Bi%2B1%2Cj%7D+-+%5Cfrac%7B1%7D%7B%5Crho%7D%28p_%7Bi%2B1%2Cj%7D+-+p_%7Bi%2Cj%7D%29%29+-+%28u_%7Bi-1%2Cj%7D+-+%5Cfrac%7B1%7D%7B%5Crho%7D%28p_%7Bi%2Cj%7D+-+p_%7Bi-1%2Cj%7D%29%29+%5C%5C+%2B+%28v_%7Bi%2Cj%2B1%7D+-+%5Cfrac%7B1%7D%7B%5Crho%7D%28p_%7Bi%2Cj%2B1%7D+-+p_%7Bi%2Cj%7D%29%29+-+%28v_%7Bi%2Cj%2B1%7D+-+%5Cfrac%7B1%7D%7B%5Crho%7D%28p_%7Bi%2Cj%7D+-+p_%7Bi%2Cj-1%7D%29%29%3D+0
改变一下各项位置，得出迭代需要的方程，同样，等号左边的p是本次迭代的结果，右边的p是上次迭代的结果

https://www.zhihu.com/equation?tex=p_%7Bi%2Cj%7D+%3D+%5Cfrac%7B1%7D%7B4%7D%28p_%7Bi%2B1%2Cj%7D%2Bp_%7Bi-1%2Cj%7D%2Bp_%7Bi%2Cj%2B1%7D%2Bp_%7Bi%2Cj-1%7D-%5Crho%28u_%7Bi%2B1%2Cj%7D+-+u_%7Bi-1%2Cj%7D+%2B+v_%7Bi%2Cj%2B1%7D+-+v_%7Bi%2Cj-1%7D%29%29%5Ctag%7B2%7D
上式密度rho所乘的括号内的式子，正好是中间速度的散度，因此将其写为divergence，它在迭代过程中并不变化，所以在迭代之前就应该计算好防止重复计算

https://www.zhihu.com/equation?tex=divergence+%3D+%5Cfrac%7B1%7D%7B2%7D%28u_%7Bi%2B1%2Cj%7D+-+u_%7Bi-1%2Cj%7D+%2B+v_%7Bi%2Cj%2B1%7D+-+v_%7Bi%2Cj-1%7D%29+%5Ctag%7B3%7D
2式又被称作压力泊松方程(Pressure Possion Equation)，或者说是雅可比迭代。泊松方程和第一章介绍的拉普拉斯方程的区别是，前者有源，后者无源。这里的源就是密度乘以散度。
压力泊松方程也可从速度分解的角度考虑，下面是一个速度分解：

https://www.zhihu.com/equation?tex=V+%3D+V_%7B%E6%97%A0%E6%95%A3%E5%BA%A6%7D+%2B+V_%7B%E6%97%A0%E6%97%8B%E5%BA%A6%7D+%5C%5 无旋度场可以视为标量场的梯度，这里的标量场就是压力。将它们同乘以散度

https://www.zhihu.com/equation?tex=%5Cnabla+%5Ccdot+V+%3D+%5Cnabla+%5Ccdot+V_%7B%E6%97%A0%E6%95%A3%E5%BA%A6%7D+%2B+%5Cnabla+%5Ccdot+%5Cnabla+p+%3D+%5Cnabla+%5Ccdot+%5Cnabla+p+%5C%5C 那么无散度的散度是零，因此这项删掉。散度与梯度的联系请看这篇https://zhuanlan.zhihu.com/p/136836187。上面这个式子我们就建立起了压力和速度的关系，两个倒三角形就是拉普拉斯运算符。写成微分形式就是

https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+V%7D%7B%5Cpartial+x%7D+%2B+%5Cfrac%7B%5Cpartial+V%7D%7B%5Cpartial+y%7D+%3D+%5Cfrac%7B%5Cpartial%5E2+p%7D%7B%5Cpartial+x%5E2%7D+%2B++%5Cfrac%7B%5Cpartial%5E2+p%7D%7B%5Cpartial+y%5E2%7D+%5C%5C
继续写成离散形式就是，然后加上密度项，变换一下即可得到1式

https://www.zhihu.com/equation?tex=%5Cfrac%7Bu_%7Bi%2B1%2Cj%7D+-+u_%7Bi-1%2Cj%7D%7D%7B2%5CDelta+x%7D+%2B+%5Cfrac%7Bv_%7Bi%2Cj-1%7D+-+v_%7Bi%2Cj-1%7D%7D%7B2%5CDelta+y%7D+%3D+%5C%5C++%5Cfrac%7B%28p_%7Bi%2B1%2Cj%7D+-+p_%7Bi%2Cj%7D%29+%2B+%28p_%7Bi%2Cj%7D+-+p_%7Bi-1%2Cj%7D%29%7D%7B2%5CDelta+x%5E2%7D+%2B+%5Cfrac%7B%28p_%7Bi%2Cj%2B1%7D+-+p_%7Bi%2Cj%7D%29+%2B+%28p_%7Bi%2Cj%7D+-+p_%7Bi%2Cj%2B1%7D%29%7D%7B2%5CDelta+y%5E2%7D
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;
//x存x轴速度，y存y轴速度
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 + loat2(-_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存储负数。
最后的坑爹的导入图片问题，要允许读写，并且更改到合适的格式
完整的代码在这个仓库的第五次提交：
然后我们就可以破坏世界，扭曲现实了！【误】
速度图如下：
可视化

我们可以分析一下刚才提到的那个webgl库，修改最终的输出图像代码如下
// display result
gl.viewport(0, 0, gl.drawingBufferWidth, gl.drawingBufferHeight);
displayProgram.bind();
gl.uniform1i(displayProgram.uniforms.uTexture, density.first);//最终颜色
//gl.uniform1i(displayProgram.uniforms.uTexture, divergence)//散度
//gl.uniform1i(displayProgram.uniforms.uTexture, pressure.first);//压力
//gl.uniform1i(displayProgram.uniforms.uTexture, velocity.second);//中间速度
//gl.uniform1i(displayProgram.uniforms.uTexture, velocity.first);//最终速度
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方程模拟粘性染料流动

查看完整版本: 【游戏流体力学基础及Unity代码（五）】用欧拉方程模拟无粘性染料之代码实现