【游戏流体力学基础及Unity代码（六）】用NavierStokes方程模拟粘性染料流动
粘度推荐看这篇文章之前先看看这个生动的的视频https://www.bilibili.com/video/BV1pW411X7xE?p=2。
很久很久以前，牛顿爵士闲来无聊，去做了一个实验。他将流体放在两个平行的平板之间，固定下层平板而以速度U移动上层平板。如下图，
U是上层平板的速度，u是流体的速度，L是平板间的距离。牛顿发现，越靠近上层平板的流体，它的速度越接近U，简单来说，就是流体的速度随到上层平板的距离发生变化：
https://www.zhihu.com/equation?tex=u%28y%29+%3D+%5Cfrac%7BU%7D%7BL%7Dy%5C%5C 上式y就是到下层平板的距离。显然是有什么力导致了流体的移动，这就是粘性力(viscous force)，和固体间的动摩擦力很相似。如何计算这种力呢？
我们之前提过压力(pressure)，它能让可压缩流体产生线变形，如下图：
压力是怎么算的呢？用力的总大小除以面积就算到了压力，压力的单位是Pa，也就是(N/m^2)牛顿每平方米。压力是沿着流体表面的法向量作用的，所以也叫做正应力和法向应力(normal stress)。
然而牛顿的这个实验力，力是沿着表面的切向量产生作用了。我们也用这个力的大小除以面积，算出来的东西就叫切向应力或剪切力(shear stress)，它和平板的速度和平板之间的距离有关：
https://www.zhihu.com/equation?tex=shear+%5C+stress+%3D+%5Cmu+%5Cfrac%7BU%7D%7BL%7D+%5C%5C 上式那个miu是个比例系数，叫做动力粘性系数(dynamic viscosity)。而(U/L)就是剪切变形速率(rate of shear deformation)。同样写成微分形式：
https://www.zhihu.com/equation?tex=%5Ctau+%3D+%5Cmu+%5Cfrac%7B%5CDelta+u%7D%7B%5CDelta+y%7D+%5C%5C 上式说明如果剪切力不变，如果流体越粘，那么移动的距离越小。比如清水粘性很小，甘油粘度很大。这就是牛顿粘性定律(Newton`s Law of Viscosity)。牛顿做实验的那种流体，粘度可以看做是常数，因此可以叫牛顿流体(Newtonian fluid)。
然而流体和流体之间是不能一概而论的，比如由玉米淀粉和水调制而成的玉米糊(oobleck)，如果我们将手慢慢伸进去，我们会感觉它没那么粘，然而如果用力快速给它一锤子，我们会觉得锤子在锤一个很硬的东西。粘度会随着的剪切力的大小变化而变化的流体，就是非牛顿流体(Non-Newtonian fluid)。
比如下图，如果塑料盒子里装的是清水，那么子弹可以直接打穿盒底。然而遇上非牛顿流体，子弹可能只好洗洗睡了。
上面的动态图出自https://www.youtube.com/watch?v=bLiNHqwgWaQ。关于非牛顿流体的一篇很好的科普请看https://www.zhihu.com/question/267108394/answer/559562147。
粘性力
对于xy平面上的二维流体，它会受到X轴方向的力，导致X轴方向上的速度沿着Y轴变化。同理也会受到Y轴方向的力，也就说：
https://www.zhihu.com/equation?tex=%5Ctau_%7Bxy%7D+%3D+%5Ctau_%7Byx%7D+%3D+%5Cmu+%28%5Cfrac%7B%5CDelta+v%7D%7B%5CDelta+x%7D+%2B+%5Cfrac%7B%5CDelta+u%7D%7B%5CDelta+y%7D%29%5C%5C 而如果力作用在X轴方向，导致X轴方向上的速度沿着X轴变化，那么这就是压力。这时候再看看第四章说到的速度分解，由线变形和角变形导致的第二项，你是否想到了什么呢？
https://www.zhihu.com/equation?tex=%5Cbegin%7Bpmatrix%7D+%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+x%7D%26+%5Cfrac%7B1%7D%7B2%7D%28%5Cfrac%7B+%5Cpartial+u%7D%7B%5Cpartial+y%7D+%2B+%5Cfrac%7B+%5Cpartial+v%7D%7B%5Cpartial+x%7D%29+%5C%5C+++%5Cfrac%7B1%7D%7B2%7D%28%5Cfrac%7B+%5Cpartial+u%7D%7B%5Cpartial+y%7D+%2B+%5Cfrac%7B+%5Cpartial+v%7D%7B%5Cpartial+x%7D%29%26+%5Cfrac%7B%5Cpartial+v%7D%7B%5Cpartial+y%7D%5C%5C++%5Cend%7Bpmatrix%7D%5C%5C 动力粘性系数mu经常和密度rho作为未知数一起出现，并且mu的大小不能直接反映不同流体的易流动性程度，因此引入运动粘性系数(kinematic viscosity)nu。按照开头所给视频的说法，动力粘性系数可以看作是相同变形率时的粘性大小，而运动粘性系数可以看做是相同加速度时粘性的大小：
https://www.zhihu.com/equation?tex=%5Cnu+%3D+%5Cfrac%7B%5Cmu%7D%7B%5Crho%7D+%5C%5C 下面是一些常见的物质在一个标准大气压下的运动粘性系数参考：
粘性力大部分时候都是非常小的，为了简化计算，有时候必须是被忽略的。但是是什么时候呢？我们至今讨论的力中，可以分为惯性力(Inertia Force)和粘性力(Viscous Force)，前者是让流体保持前进速度的力，而后者是让流体速度稳定下来的力。给两者取个比值，就是雷诺数(Reynold Number)。
https://www.zhihu.com/equation?tex=%5Cfrac%7B%E6%83%AF%E6%80%A7%E5%8A%9B%7D%7B%E7%B2%98%E6%80%A7%E5%8A%9B%7D+%3D+%5Cfrac%7B%5Crho+uu%7D%7B%5Cmu+u%2F+L%7D+%3D+%5Cfrac%7B%5Crho+u+L%7D%7B%5Cmu%7D+%3D+%5Cfrac%7Bu+L%7D%7B%5Cnu%7D+%3D+Re+%5C%5C 当雷诺数很小的，也就是粘性力相比惯性力很大的时候，流体微团如果想突然改变速度，都会被附近的流体微团把速度改回来一些，此时流体速度很稳定，这就是平滑的层流(laminar flow)。而当雷诺数很大，粘性力远远小于惯性力，甚至可以忽略掉。此时流体微团随机游走，就变成了复杂的紊流(turbulent flow)。
Naiver-Stokes方程
先回顾第四章关于动量推导的部分。在Euler方程中，我们考虑的表面力只有压力，现在我们还要加上粘性力，并且由散度定理将它从面积分变成体积分。这里的上三角符号是Delta的意思，下三角符号是散度符号。
你可以把粘度想象成动量的扩散。动量多的地方会把自己的动量送给周围动量少的地方。然后把这项加在Euler方程后面，就得了著名的不可压缩的Navier-Stokes公式。注意动力粘性系数被密度除后就变成了运动粘性系数。
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+t%7D%2Bu%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+x%7D%2Bv%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+y%7D+%3D+-%5Cfrac%7B1%7D%7B%5Crho%7D%5Cfrac%7B%5Cpartial+p%7D%7B%5Cpartial+x%7D%2B%5Cnu+%5Cleft%28%5Cfrac%7B%5Cpartial%5E2+u%7D%7B%5Cpartial+x%5E2%7D%2B%5Cfrac%7B%5Cpartial%5E2+u%7D%7B%5Cpartial+y%5E2%7D+%5Cright%29%5Ctag%7B1%7D
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+v%7D%7B%5Cpartial+t%7D%2Bu%5Cfrac%7B%5Cpartial+v%7D%7B%5Cpartial+x%7D%2Bv%5Cfrac%7B%5Cpartial+v%7D%7B%5Cpartial+y%7D+%3D+-%5Cfrac%7B1%7D%7B%5Crho%7D%5Cfrac%7B%5Cpartial+p%7D%7B%5Cpartial+y%7D%2B%5Cnu%5Cleft%28%5Cfrac%7B%5Cpartial%5E2+v%7D%7B%5Cpartial+x%5E2%7D%2B%5Cfrac%7B%5Cpartial%5E2+v%7D%7B%5Cpartial+y%5E2%7D%5Cright%29%5Ctag%7B2%7D
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+u%7D%7B%5Cpartial+x%7D%2B%5Cfrac%7B%5Cpartial+v%7D%7B%5Cpartial+y%7D+%3D+0+%5C%5C 这里我们还是忽略了重力和其它力的影响。如果你对上面这个式子不熟悉，那么请先看看本系列第四章的动量方程推导。大体计算步骤和Euler方程差不多，先算平流和鼠标拖拽产生的力
https://www.zhihu.com/equation?tex=%5Cfrac%7BV%5E%7B%2A%7D+-+V%5En%7D%7B%5CDelta+t%7D+%3D+F%5En%5C%5C+%5C%5C 然后我们现在再计算粘性力。
https://www.zhihu.com/equation?tex=%5Cfrac%7BV%5E%7B%2A%2A%7D+-+V%5E%2A%7D%7B%5CDelta+t%7D+%3D+%5Cnu+%28%5CDelta+V%5E%2A%29+%5C%5C 粘性会使流体微团的速度受到周围流体微团的影响，对于二维平面，离散化就是：
https://www.zhihu.com/equation?tex=%7B%7Bu_%7Bi%2Cj%7D%5E%7Bn%2B1%7D+-+u_%7Bi%2Cj%7D%5En%7D+%5Cover+%7B%5CDelta+t%7D%7D+%3D+%5Cnu+%7B%7Bu_%7Bi-1%2Cj%7D%5En+-+2u_%7Bi%2Cj%7D%5En+%2B+u_%7Bi%2B1%2Cj%7D%5En%7D+%5Cover+%5CDelta+x%5E2%7D+%2B+%5Cnu+%7B%7Bu_%7Bi%2Cj-1%7D%5En+-+2u_%7Bi%2Cj%7D%5En+%2B+u_%7Bi%2Cj%2B1%7D%5En%7D+%5Cover+%5CDelta+y%5E2%7D+%5C%5C
https://www.zhihu.com/equation?tex=%7B%7Bv_%7Bi%2Cj%7D%5E%7Bn%2B1%7D+-+v_%7Bi%2Cj%7D%5En%7D+%5Cover+%7B%5CDelta+t%7D%7D+%3D+%5Cnu+%7B%7Bv_%7Bi-1%2Cj%7D%5En+-+2v_%7Bi%2Cj%7D%5En+%2B+v_%7Bi%2B1%2Cj%7D%5En%7D+%5Cover+%5CDelta+x%5E2%7D+%2B+%5Cnu+%7B%7Bv_%7Bi%2Cj-1%7D%5En+-+2v_%7Bi%2Cj%7D%5En+%2B+v_%7Bi%2Cj%2B1%7D%5En%7D+%5Cover+%5CDelta+y%5E2%7D%5C%5C 上面步骤可以不用迭代直接计算。但按照GPU GEM3第三十八章介绍的方法，我们还是把它当成一个泊松方程来多次迭代来计算出速度。最后一步就是求解压力，算出无散度的最终速度：
代码实现
继续在第五章的基础上改动。cs文件里，添加粘性迭代相关的代码
//第三步：计算受到粘性影响的速度，得到中间速度
for (int i = 0; i < 10; i++)
{
ViscosityMat.SetTexture(&#34;_VelocityTex&#34;, VelocityRT2);
Graphics.Blit(VelocityRT2, VelocityRT, ViscosityMat);
Graphics.Blit(VelocityRT, VelocityRT2);
}然后是相应的Viscosity.Shader。
float2 Top = tex2D(_VelocityTex, i.uv + float2(0.0f, _VelocityTex_TexelSize.y)).xy;
float2 Bottom = tex2D(_VelocityTex, i.uv + float2(0.0f, -_VelocityTex_TexelSize.y)).xy;
float2 Right = tex2D(_VelocityTex, i.uv + float2(_VelocityTex_TexelSize.x, 0.0f)).xy;
float2 Left = tex2D(_VelocityTex, i.uv + float2(-_VelocityTex_TexelSize.x, 0.0f)).xy;
float2 Center = tex2D(_VelocityTex, i.uv).xy;
float nu = 0.1;//运动粘性系数
float VelocityX = Center.x + nu * ((Right.x - 2 * Center.x + Left.x) + (Top.x - 2 * Center.x + Bottom.x));
float VelocityY = Center.y + nu * ((Right.y - 2 * Center.y + Left.y) + (Top.y - 2 * Center.y + Bottom.y));
return float4(VelocityX, VelocityY, 0.0f, 0.0f);然后我们就能继续毁灭这个世界了【误】
完整代码在下面，部分着色器我直接使用之前的了
再现牛顿当年的实验
虽说是再现，但是稍微有一点不同。同样假设上平板以速度U向右走，下平板不动。但是现在我们还要加上左平板和右平板，左右平板也不会移动。这就是模拟流体中经典的Lid-Driven Cavity Flow问题。
现在我们试试CPU模拟。github上这种代码一搜一大把，都是python和matlab的，但初次接触可能仍然有些困惑。比如当你看12步NS方程的第14章的时候ttps://nbviewer.jupyter.org/github/barbagroup/CFDPython/blob/master/lessons/14_Step_11.ipynb。我把它转写成unity上的代码如下，然后慢慢分析。
这个问题的解决的方法和我们之前说的差不多，唯一不同的就是边界条件。由于我们模拟的是粘性流体，因此边界平板的速度会带动边界处的流体，边界流体的速度等于这里平板上的速度，这就是无滑移(No Slip)边界条件，写成公式如下：
https://www.zhihu.com/equation?tex=V_%7B%E8%BE%B9%E7%95%8C%7D+%3D+a+%5C%5C a就是平板的速度。这又被称作是狄利克雷边界条件(Dirichlet boundary condition)。Y方向的速度在边界上永远都是零，因为左右边界并未移动，而X方向上，下边界速度应该是零，而上边界速度应该是上平板移动的速度，我们设置为1，对应代码在106~119行。
而压力的边界条件是什么呢？看看求解的最后一步，也就说下面这个式子。
在边界处，流体的速度保持不变，因此上式左边是零，那么右边也应该是零。也就说边界处压力的梯度是零。写成公式如下：
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+p%7D%7B%5Cpartial+x%7D+%3D+%5Cfrac%7B%5Cpartial+p%7D%7B%5Cpartial+y%7D+%3D+0+%5C%5C 上面又被叫做黎曼边界条件(Neumann boundary conditions)。对应代码在第84到第87行。
上面那个问题的最后求解结果是流体在旋转，在真实世界中它看起来是这样子的，图源自维基https://en.wikipedia.org/wiki/Wingtip_vortices
至于为什么会这样，又是一个很长的故事了，我们之后的篇章再讨论。
接下来我再展示一下如何把公式写成代码，公式的推导请看12步NS方程的第14章。这里我把之前GPU模拟中缺失的很多项都补上来了，比如dx不再是1了。首先是计算压力的公式，相应代码在79~80行。这里的压力同样是猜的，所以需要迭代很多次。
https://www.zhihu.com/equation?tex=%5Cbegin%7Bsplit%7D+p_%7Bi%2Cj%7D%5E%7Bn%7D+%3D+%26+%5Cfrac%7B%5Cleft%28p_%7Bi%2B1%2Cj%7D%5E%7Bn%7D%2Bp_%7Bi-1%2Cj%7D%5E%7Bn%7D%5Cright%29+%5CDelta+y%5E2+%2B+%5Cleft%28p_%7Bi%2Cj%2B1%7D%5E%7Bn%7D%2Bp_%7Bi%2Cj-1%7D%5E%7Bn%7D%5Cright%29+%5CDelta+x%5E2%7D%7B2%5Cleft%28%5CDelta+x%5E2%2B%5CDelta+y%5E2%5Cright%29%7D+++-%5Cfrac%7B%5Crho%5CDelta+x%5E2%5CDelta+y%5E2%7D%7B2%5Cleft%28%5CDelta+x%5E2%2B%5CDelta+y%5E2%5Cright%29%7D+%5Ctimes+b_%7Bi%2Cj%7D%5E%7Bn%7D+%5Cend%7Bsplit%7D
p = (((pn + pn) *
(dx * dx) + (pn + pn) * (dy * dy)) / (2 * (dx * dx + dy * dy)) -
rho * (b * ((dy * dx) * (dy * dx <span class="o">/ (2 * (dx * dx + dy * dy))));其中b和我们之前说过的散度很相似，是为了减少不必要的计算的，代码在第63~66行
https://www.zhihu.com/equation?tex=b_%7Bi%2Cj%7D%5E%7Bn%7D+%3D+%5Cleft%5B%5Cfrac%7B1%7D%7B%5CDelta+t%7D%5Cleft%28%5Cfrac%7Bu_%7Bi%2B1%2Cj%7D-u_%7Bi-1%2Cj%7D%7D%7B2%5CDelta+x%7D%2B%5Cfrac%7Bv_%7Bi%2Cj%2B1%7D-v_%7Bi%2Cj-1%7D%7D%7B2%5CDelta+y%7D%5Cright%29-%5Cfrac%7Bu_%7Bi%2B1%2Cj%7D-u_%7Bi-1%2Cj%7D%7D%7B2%5CDelta+x%7D%5Cfrac%7Bu_%7Bi%2B1%2Cj%7D-u_%7Bi-1%2Cj%7D%7D%7B2%5CDelta+x%7D+%5C%5C-2%5Cfrac%7Bu_%7Bi%2Cj%2B1%7D-u_%7Bi%2Cj-1%7D%7D%7B2%5CDelta+y%7D%5Cfrac%7Bv_%7Bi%2B1%2Cj%7D-v_%7Bi-1%2Cj%7D%7D%7B2%5CDelta+x%7D-%5Cfrac%7Bv_%7Bi%2Cj%2B1%7D-v_%7Bi%2Cj-1%7D%7D%7B2%5CDelta+y%7D%5Cfrac%7Bv_%7Bi%2Cj%2B1%7D-v_%7Bi%2Cj-1%7D%7D%7B2%5CDelta+y%7D%5Cright%5D
b = (u - u) / (2 * dx) + (v - v) / (2 * dy);
b -= ((((u - u) / (2 * dx)) * ((u - u) /
(2 * dx)))+ 2 * (((u - u) / (2 * dy)) * ((v - v) /
(2 * dy))) + (((v - v) / (2 * dy)) * ((v - v) /
(2 * dy))));我们把1式写成离散形式，代码在95~98行
https://www.zhihu.com/equation?tex=%5Cbegin%7Bsplit%7D+u_%7Bi%2Cj%7D%5E%7Bn%2B1%7D+%3D+u_%7Bi%2Cj%7D%5E%7Bn%7D+%26+-+u_%7Bi%2Cj%7D%5E%7Bn%7D+%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+x%7D+%5Cleft%28u_%7Bi%2Cj%7D%5E%7Bn%7D-u_%7Bi-1%2Cj%7D%5E%7Bn%7D%5Cright%29+-+v_%7Bi%2Cj%7D%5E%7Bn%7D+%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+y%7D+%5Cleft%28u_%7Bi%2Cj%7D%5E%7Bn%7D-u_%7Bi%2Cj-1%7D%5E%7Bn%7D%5Cright%29+%5C%5C+%26+-+%5Cfrac%7B%5CDelta+t%7D%7B%5Crho+2%5CDelta+x%7D+%5Cleft%28p_%7Bi%2B1%2Cj%7D%5E%7Bn%7D-p_%7Bi-1%2Cj%7D%5E%7Bn%7D%5Cright%29+%5C%5C+%26+%2B+%5Cnu+%5Cleft%28%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+x%5E2%7D+%5Cleft%28u_%7Bi%2B1%2Cj%7D%5E%7Bn%7D-2u_%7Bi%2Cj%7D%5E%7Bn%7D%2Bu_%7Bi-1%2Cj%7D%5E%7Bn%7D%5Cright%29+%2B+%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+y%5E2%7D+%5Cleft%28u_%7Bi%2Cj%2B1%7D%5E%7Bn%7D-2u_%7Bi%2Cj%7D%5E%7Bn%7D%2Bu_%7Bi%2Cj-1%7D%5E%7Bn%7D%5Cright%29%5Cright%29+%5Cend%7Bsplit%7D
float advection = (un - (un * (dt / dx) * (un - un)) -
vn * (dt / dy) * (un - un));//平流项
float diffuse = (nu * dt / (dx * dx)) * (un - 2 * un + un) +
(nu * dt / (dy * dx)) * (un - 2 * un + un);//粘性项
float density = (dt / (2 * rho * dx)) * (p - p);//密度项
u = advection + diffuse - density;2式也写成离散形式，代码在第99到103行
https://www.zhihu.com/equation?tex=%5Cbegin%7Bsplit%7D+v_%7Bi%2Cj%7D%5E%7Bn%2B1%7D+%3D+v_%7Bi%2Cj%7D%5E%7Bn%7D+%26+-+u_%7Bi%2Cj%7D%5E%7Bn%7D+%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+x%7D+%5Cleft%28v_%7Bi%2Cj%7D%5E%7Bn%7D-v_%7Bi-1%2Cj%7D%5E%7Bn%7D%5Cright%29+-+v_%7Bi%2Cj%7D%5E%7Bn%7D+%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+y%7D+%5Cleft%28v_%7Bi%2Cj%7D%5E%7Bn%7D-v_%7Bi%2Cj-1%7D%5E%7Bn%7D%29%5Cright%29+%5C%5C+%26+-+%5Cfrac%7B%5CDelta+t%7D%7B%5Crho+2%5CDelta+y%7D+%5Cleft%28p_%7Bi%2Cj%2B1%7D%5E%7Bn%7D-p_%7Bi%2Cj-1%7D%5E%7Bn%7D%5Cright%29+%5C%5C+%26+%2B+%5Cnu+%5Cleft%28%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+x%5E2%7D+%5Cleft%28v_%7Bi%2B1%2Cj%7D%5E%7Bn%7D-2v_%7Bi%2Cj%7D%5E%7Bn%7D%2Bv_%7Bi-1%2Cj%7D%5E%7Bn%7D%5Cright%29+%2B+%5Cfrac%7B%5CDelta+t%7D%7B%5CDelta+y%5E2%7D+%5Cleft%28v_%7Bi%2Cj%2B1%7D%5E%7Bn%7D-2v_%7Bi%2Cj%7D%5E%7Bn%7D%2Bv_%7Bi%2Cj-1%7D%5E%7Bn%7D%5Cright%29%5Cright%29+%5Cend%7Bsplit%7D
advection = (vn - (un * dt / dx * (vn - vn)) -
vn * dt / dy * (vn - vn));
diffuse = (nu * dt / (dx * dx)) * (vn - 2 * vn + vn) +
(nu * dt / (dx * dx)) * (vn - 2 * vn + vn);
density = (dt / (2 * rho * dy)) * (p - p);
v = advection + diffuse - density;最终速度U如下，效果非常丑，因为横竖只有10个网格，而且配色没有python和matlab的库或自带函数配得好看。这里红色代表正数，蓝色代表负数。
更漂亮应该看看下图。
至此，Navier-Stokes方程的入门之旅就结束了。不过先别急着完结撒花，我们能在unity上玩的东西还有很多。
下一篇：光影帽子：【游戏流体力学基础及Unity代码（七）】车流量问题，非线性水波以及burgers方程
上一篇：光影帽子：【游戏流体力学基础及Unity代码（五）】用欧拉方程模拟无粘性染料之代码实现
页:
[1]