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

snail · 发表于 2020-11-25 09:26

粘度

推荐看这篇文章之前先看看这个生动的的视频https://www.bilibili.com/video/BV1pW411X7xE?p=2。
很久很久以前，牛顿爵士闲来无聊，去做了一个实验。他将流体放在两个平行的平板之间，固定下层平板而以速度U移动上层平板。如下图，
U是上层平板的速度，u是流体的速度，L是平板间的距离。牛顿发现，越靠近上层平板的流体，它的速度越接近U，简单来说，就是流体的速度随到上层平板的距离发生变化：

$u(y) = \frac{U}{L}y\\$ 上式y就是到下层平板的距离。显然是有什么力导致了流体的移动，这就是粘性力(viscous force)，和固体间的动摩擦力很相似。如何计算这种力呢？
我们之前提过压力(pressure)，它能让可压缩流体产生线变形，如下图：
压力是怎么算的呢？用力的总大小除以面积就算到了压力，压力的单位是Pa，也就是(N/m^2)牛顿每平方米。压力是沿着流体表面的法向量作用的，所以也叫做正应力和法向应力(normal stress)。
然而牛顿的这个实验力，力是沿着表面的切向量产生作用了。我们也用这个力的大小除以面积，算出来的东西就叫切向应力或剪切力(shear stress)，它和平板的速度和平板之间的距离有关：

$shear \ stress = \mu \frac{U}{L} \\$ 上式那个miu是个比例系数，叫做动力粘性系数(dynamic viscosity)。而(U/L)就是剪切变形速率(rate of shear deformation)。同样写成微分形式：

$\tau = \mu \frac{\Delta u}{\Delta y} \\$ 上式说明如果剪切力不变，如果流体越粘，那么移动的距离越小。比如清水粘性很小，甘油粘度很大。这就是牛顿粘性定律(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轴方向的力，也就说：

$\tau_{xy} = \tau_{yx} = \mu (\frac{\Delta v}{\Delta x} + \frac{\Delta u}{\Delta y})\\$ 而如果力作用在X轴方向，导致X轴方向上的速度沿着X轴变化，那么这就是压力。这时候再看看第四章说到的速度分解，由线变形和角变形导致的第二项，你是否想到了什么呢？

$\begin{pmatrix} \frac{\partial u}{\partial x}& \frac{1}{2}(\frac{ \partial u}{\partial y} + \frac{ \partial v}{\partial x}) \\ \frac{1}{2}(\frac{ \partial u}{\partial y} + \frac{ \partial v}{\partial x})& \frac{\partial v}{\partial y}\\ \end{pmatrix}\\$ 动力粘性系数mu经常和密度rho作为未知数一起出现，并且mu的大小不能直接反映不同流体的易流动性程度，因此引入运动粘性系数(kinematic viscosity)nu。按照开头所给视频的说法，动力粘性系数可以看作是相同变形率时的粘性大小，而运动粘性系数可以看做是相同加速度时粘性的大小：

$\nu = \frac{\mu}{\rho} \\$ 下面是一些常见的物质在一个标准大气压下的运动粘性系数参考：
粘性力大部分时候都是非常小的，为了简化计算，有时候必须是被忽略的。但是是什么时候呢？我们至今讨论的力中，可以分为惯性力(Inertia Force)和粘性力(Viscous Force)，前者是让流体保持前进速度的力，而后者是让流体速度稳定下来的力。给两者取个比值，就是雷诺数(Reynold Number)。

$\frac{惯性力}{粘性力} = \frac{\rho uu}{\mu u/ L} = \frac{\rho u L}{\mu} = \frac{u L}{\nu} = Re \\$ 当雷诺数很小的，也就是粘性力相比惯性力很大的时候，流体微团如果想突然改变速度，都会被附近的流体微团把速度改回来一些，此时流体速度很稳定，这就是平滑的层流(laminar flow)。而当雷诺数很大，粘性力远远小于惯性力，甚至可以忽略掉。此时流体微团随机游走，就变成了复杂的紊流(turbulent flow)。
Naiver-Stokes方程

先回顾第四章关于动量推导的部分。在Euler方程中，我们考虑的表面力只有压力，现在我们还要加上粘性力，并且由散度定理将它从面积分变成体积分。这里的上三角符号是Delta的意思，下三角符号是散度符号。
你可以把粘度想象成动量的扩散。动量多的地方会把自己的动量送给周围动量少的地方。然后把这项加在Euler方程后面，就得了著名的不可压缩的Navier-Stokes公式。注意动力粘性系数被密度除后就变成了运动粘性系数。

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right)\tag{1}$

$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial y}+\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right)\tag{2}$

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} = 0 \\$ 这里我们还是忽略了重力和其它力的影响。如果你对上面这个式子不熟悉，那么请先看看本系列第四章的动量方程推导。大体计算步骤和Euler方程差不多，先算平流和鼠标拖拽产生的力

$\frac{V^{*} - V^n}{\Delta t} = F^n\\ \\$ 然后我们现在再计算粘性力。

$\frac{V^{**} - V^*}{\Delta t} = \nu (\Delta V^*) \\$ 粘性会使流体微团的速度受到周围流体微团的影响，对于二维平面，离散化就是：

${{u_{i,j}^{n+1} - u_{i,j}^n} \over {\Delta t}} = \nu {{u_{i-1,j}^n - 2u_{i,j}^n + u_{i+1,j}^n} \over \Delta x^2} + \nu {{u_{i,j-1}^n - 2u_{i,j}^n + u_{i,j+1}^n} \over \Delta y^2} \\$

${{v_{i,j}^{n+1} - v_{i,j}^n} \over {\Delta t}} = \nu {{v_{i-1,j}^n - 2v_{i,j}^n + v_{i+1,j}^n} \over \Delta x^2} + \nu {{v_{i,j-1}^n - 2v_{i,j}^n + v_{i,j+1}^n} \over \Delta y^2}\\$ 上面步骤可以不用迭代直接计算。但按照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)边界条件，写成公式如下：

$V_{边界} = a \\$ a就是平板的速度。这又被称作是狄利克雷边界条件(Dirichlet boundary condition)。Y方向的速度在边界上永远都是零，因为左右边界并未移动，而X方向上，下边界速度应该是零，而上边界速度应该是上平板移动的速度，我们设置为1，对应代码在106~119行。
而压力的边界条件是什么呢？看看求解的最后一步，也就说下面这个式子。

在边界处，流体的速度保持不变，因此上式左边是零，那么右边也应该是零。也就说边界处压力的梯度是零。写成公式如下：

$\frac{\partial p}{\partial x} = \frac{\partial p}{\partial y} = 0 \\$ 上面又被叫做黎曼边界条件(Neumann boundary conditions)。对应代码在第84到第87行。
上面那个问题的最后求解结果是流体在旋转，在真实世界中它看起来是这样子的，图源自维基https://en.wikipedia.org/wiki/Wingtip_vortices
至于为什么会这样，又是一个很长的故事了，我们之后的篇章再讨论。
接下来我再展示一下如何把公式写成代码，公式的推导请看12步NS方程的第14章。这里我把之前GPU模拟中缺失的很多项都补上来了，比如dx不再是1了。首先是计算压力的公式，相应代码在79~80行。这里的压力同样是猜的，所以需要迭代很多次。

$\begin{split} p_{i,j}^{n} = & \frac{\left(p_{i+1,j}^{n}+p_{i-1,j}^{n}\right) \Delta y^2 + \left(p_{i,j+1}^{n}+p_{i,j-1}^{n}\right) \Delta x^2}{2\left(\Delta x^2+\Delta y^2\right)} -\frac{\rho\Delta x^2\Delta y^2}{2\left(\Delta x^2+\Delta y^2\right)} \times b_{i,j}^{n} \end{split}$

p[i, j] = (((pn[i, j - 1] + pn[i, j + 1]) *
(dx * dx) + (pn[i - 1, j] + pn[i + 1, j]) * (dy * dy)) / (2 * (dx * dx + dy * dy)) -
rho * (b[i, j] * ((dy * dx) * (dy * dx <span class="o">/ (2 * (dx * dx + dy * dy))));

复制代码

其中b和我们之前说过的散度很相似，是为了减少不必要的计算的，代码在第63~66行

$b_{i,j}^{n} = \left[\frac{1}{\Delta t}\left(\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}+\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right)-\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} \\-2\frac{u_{i,j+1}-u_{i,j-1}}{2\Delta y}\frac{v_{i+1,j}-v_{i-1,j}}{2\Delta x}-\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right]$

b[i, j] = (u[i + 1, j] - u[i - 1, j]) / (2 * dx) + (v[i, j + 1] - v[i, j - 1]) / (2 * dy);
b[i, j] -= ((((u[i + 1, j] - u[i - 1, j]) / (2 * dx)) * ((u[i + 1, j] - u[i - 1, j]) /
(2 * dx)))+ 2 * (((u[i, j + 1] - u[i, j - 1]) / (2 * dy)) * ((v[i + 1, j] - v[i - 1, j]) /
(2 * dy))) + (((v[i, j + 1] - v[i, j - 1]) / (2 * dy)) * ((v[i, j + 1] - v[i, j - 1]) /
(2 * dy))));

复制代码

我们把1式写成离散形式，代码在95~98行

$\begin{split} u_{i,j}^{n+1} = u_{i,j}^{n} & - u_{i,j}^{n} \frac{\Delta t}{\Delta x} \left(u_{i,j}^{n}-u_{i-1,j}^{n}\right) - v_{i,j}^{n} \frac{\Delta t}{\Delta y} \left(u_{i,j}^{n}-u_{i,j-1}^{n}\right) \\ & - \frac{\Delta t}{\rho 2\Delta x} \left(p_{i+1,j}^{n}-p_{i-1,j}^{n}\right) \\ & + \nu \left(\frac{\Delta t}{\Delta x^2} \left(u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}\right) + \frac{\Delta t}{\Delta y^2} \left(u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}\right)\right) \end{split}$

float advection = (un[i, j] - (un[i, j] * (dt / dx) * (un[i, j] - un[i - 1, j])) -
vn[i, j] * (dt / dy) * (un[i, j] - un[i, j - 1]));//平流项
float diffuse = (nu * dt / (dx * dx)) * (un[i + 1, j] - 2 * un[i, j] + un[i - 1, j]) +
(nu * dt / (dy * dx)) * (un[i, j - 1] - 2 * un[i, j] + un[i, j + 1]);//粘性项
float density = (dt / (2 * rho * dx)) * (p[i + 1, j] - p[i - 1, j]);//密度项
u[i, j] = advection + diffuse - density;

复制代码

2式也写成离散形式，代码在第99到103行

$\begin{split} v_{i,j}^{n+1} = v_{i,j}^{n} & - u_{i,j}^{n} \frac{\Delta t}{\Delta x} \left(v_{i,j}^{n}-v_{i-1,j}^{n}\right) - v_{i,j}^{n} \frac{\Delta t}{\Delta y} \left(v_{i,j}^{n}-v_{i,j-1}^{n})\right) \\ & - \frac{\Delta t}{\rho 2\Delta y} \left(p_{i,j+1}^{n}-p_{i,j-1}^{n}\right) \\ & + \nu \left(\frac{\Delta t}{\Delta x^2} \left(v_{i+1,j}^{n}-2v_{i,j}^{n}+v_{i-1,j}^{n}\right) + \frac{\Delta t}{\Delta y^2} \left(v_{i,j+1}^{n}-2v_{i,j}^{n}+v_{i,j-1}^{n}\right)\right) \end{split}$

advection = (vn[i, j] - (un[i, j] * dt / dx * (vn[i, j] - vn[i - 1, j])) -
vn[i, j] * dt / dy * (vn[i, j] - vn[i, j - 1]));
diffuse = (nu * dt / (dx * dx)) * (vn[i + 1, j] - 2 * vn[i, j] + vn[i - 1, j]) +
(nu * dt / (dx * dx)) * (vn[i, j - 1] - 2 * vn[i, j] + vn[i, j + 1]);
density = (dt / (2 * rho * dy)) * (p[i, j + 1] - p[i - 1, j - 1]);
v[i, j] = advection + diffuse - density;

复制代码

最终速度U如下，效果非常丑，因为横竖只有10个网格，而且配色没有python和matlab的库或自带函数配得好看。这里红色代表正数，蓝色代表负数。
更漂亮应该看看下图。
至此，Navier-Stokes方程的入门之旅就结束了。不过先别急着完结撒花，我们能在unity上玩的东西还有很多。
下一篇：光影帽子：【游戏流体力学基础及Unity代码（七）】车流量问题，非线性水波以及burgers方程
上一篇：光影帽子：【游戏流体力学基础及Unity代码（五）】用欧拉方程模拟无粘性染料之代码实现

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

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

本帖子中包含更多资源