【游戏流体力学基础及Unity代码（十二）】卡门涡街，边界层，涡方法

六翼天使494 · 发表于 2021-1-19 11:51

模拟卡门涡街

除了圆，还有正方形以及机翼的版本，完整视频请看https://www.bilibili.com/video/BV1u5411H7hr/。你可以试试各种其它奇奇怪怪的形状。
关于卡门涡街的知识可看看李永乐老师的视频https://www.bilibili.com/video/BV1BT4y1u7DA?from=search&seid=6057852290300366231。用着色器模拟卡门涡街非常简单，简单到令人发指的地步。相比与第五章的Euler方程，重要改动只有两处。第一处就是平流密度时的像素着色器部分，与第二章介绍的差不多，最左侧的染料图像保持不变。
float4 frag(v2f i) :SV_Target{
float2 vel = tex2D(VelocityTex, i.uv).xy;
float4 col = tex2D(DensityTex, i.uv - 0.01f*vel);
if (i.uv.x < 0.1f)col = tex2D(InitDyeTex, i.uv);
return col;
}
下面是平流速度时的像素着色器部分，BlockTex就是障碍物，如果在障碍物上速度就设置为零。最左侧速度一直为1。
float4 frag(v2f i) :SV_Target{
float4 col = tex2D(VelocityTex, i.uv - 0.01f*tex2D(VelocityTex, i.uv));
if (tex2D(BlockTex, i.uv).x > 0.99f)col.xy = float2(0.0f, 0.0f);
if (i.uv.x < 0.01f)col.xy = float2(1.0f, 0.0f);
return col;
}
上面这两个平流时，计算出的上一步位置几乎不可能位于像素正中间，所以着色器内部会通过插值计算值，这样就导致了人工粘性的产生。所以无需特别添加计算粘性的着色器，也能保证上面模拟的流体有粘性。unity代码仓库地址如下
上图的锯齿是由用到着色器的jpg图片本身质量问题导致的。如果平流速度时用如下代码而不是采样图片的画就不会有锯齿产生clatterrr/FluidSimulationTutorialsUnity上图的锯齿是由用到着色器的jpg图片本身质量问题导致的。如果平流速度时用如下代码而不是采样图片的画就不会有锯齿产生
float2 dir = i.uv - float2(0.3f, 0.5f);
float dis = dir.x*dir.x + dir.y*dir.y;
if(dis < 0.001f)col.xy = float2(0.0f, 0.0f);
效果如下：
边界层

边界层分离(Boundary Layer Separation)看起来就是这样的。
关于边界层的入门理论直接看看王洪伟老师的视频就挺不错。https://www.bilibili.com/video/BV1pW411X7xE?p=15。更深入推荐一本书《Boundary-Layer Theory》 by Hermann Schlichting, Klaus Gersten。很多流体力学书比如《Fundamentals of Fluid Mechanics, 6th Edition 》by Bruce R. Munson, Donald F. Young, Theodore H. Okiishi, Wade W. Huebsch 的第8章或是B站知乎油管也有与边界层有关的东西，很容易找到。弄懂边界层为什么会发生分离就是本章读者的任务。这些到处都很容易找的知识是不会出现在本系列中的。
要抑制边界层分离的办法之一是增加流体粘性，这在我们的代码中很容易办到，如下，但是现实中无法这么做。
增大粘性会让卡门涡街消失，圆柱绕流变成如下状态。
抑制边界层分离的一个方法是涡流发生器，在本篇稍后介绍。设计高尔夫球也需要抑制边界层分离，所以把球表面做的坑坑洼洼。具体可看https://www.zhihu.com/question/23890724/answer/1179979539。足球中的香蕉球或电梯球也与边界层有关，具体可看https://zhuanlan.zhihu.com/p/37894502。
信天翁(Albatross)这种鸟类也是边界层理论大师。信天翁有时候会有好几天连续在海面上飞翔，而中途无法休息。而信天翁借助海风的力量，可以沿着海风的垂直方向飞行很久而几乎不消耗一点自身的能量。根据边界层理论，海风在靠近海面的地方，由于之前提到的边界层效果，速度非常慢接近于零，但是当离海面越来越远，海风的速度也会加快。如下图，信天翁一开始在2号位置，距离海面有一定距离，海风的速度很大。然后信天翁沿着234的轨迹飞行，一边顺风飞行，一边向下飞行，重力让信天翁下降得更快。到了4号位置，离海面很近，海风速度接近于零，信天翁的速度已经很快了，当然动能也很大。此时4号位置相当于0号位置，于是信天翁再沿着012的路径飞行，一边逆风一边克服重力的影响，到了2号位置，信天翁动能消耗殆尽，速度接近于零。此时再沿着234路径飞行，并调整自己的方向，又可以让海风为自己补充足够的动能。如此循环往复，这就是Dynamic Soaring。信天翁借助这种方法可以一天飞行1000公里甚至更远，整个过程中甚至连翅膀都不需要扇动几下。
德国流体学家Ludwig Prandtl在1904年发现了边界层。关于这位大佬有一本传记可读《Ludwig Prandtl A Life for Fluid Mechanics and Aeronautical Research》 by Michael Eckert。而且他的学生和冯卡门发现了卡门涡街。就像《上帝会掷骰子吗》那样八卦量子力学科学家的故事一样，把这些研究流体力学的大佬的故事八卦一下也很意思，可以写一本小说了，不过我在国内外都只看到零散的故事，如果有人能写个系统的小说我肯定买买买...
一点空气动力学

飞机附近的涡流主要如下
1处的漩涡也叫BoundVortex，这个漩涡与飞机的升力密切相关，这在本系列第11章已说过。由于机翼下方为高压区，机翼上方为低压区，这样的压力差别可以产生升力。
但是在三维情况下，高压气气流还是会流向低压区，这就是机翼尖涡流(wingtip trailing Vortex)。如上图的2号涡流。如下图左所示。下图是飞机正面图。
这个漩涡会把气流砸在机翼上方，从而减小升力。为了避免这个情况，可以使用翼梢小翼(winglets)，减小这个漩涡以减小升力损失。如上图右侧。但这个漩涡也有好的一面。如上图左的机翼延伸处，这个漩涡实际上产生了向上的升力，也就是红色箭头。嗯，这个升力自己虽然用不了，但是也不能白白浪费了。这就是为什么大雁要排成V字形飞行。前面大雁飞行时产生的额外升力可以留给后面的大雁用。
3处的StartingVortex，其环量正好是机翼上的BoundVortex的负数。StartingVortex其实就是那些卡门涡街。同样将粘性调大，也就是将粘性着色器循环次数增大，可以减小卡门涡街的大小和频率。StartingVortex如下图。
由于边界层分离现象，机翼表面会产生漩涡扰乱气流，造成飞机无法稳定飞行。其解决办法就是让这些漩涡远离机翼表面，这需要用到涡流发生器(Vortex Generator)。虽然用fluent或openfoam这类软件分析起来更严谨一些，但现在咱也可以使用着色器，而且涡流发生器以及各种结构画起来挺方便的。涡流发生器的作用大概就是下图这样：
上图上面是有涡流发生器的，下面则没有。上面气流经过涡流发生器后，漩涡形成在远离机翼表面的地方，机翼附近的气流较为稳定。下图漩涡形成在靠近机翼表面的地方，导致机翼表面气流不稳定。也许观看速度图或压力图更直观一些。下面是X方向速度。
下面是压力。有涡流发生器的机翼上表面附近压力稳定，而没有涡流发生器则压力忽大忽小。注意下图的锯齿是由用到着色器的jpg图片本身质量问题导致的。
老是模拟挺没意思，咱也可以折纸飞机做实验。折纸飞机也需要考虑升力，翼尖扰流等问题。可参考下面这个搬运视频。https://www.bilibili.com/video/BV1XK411u7bS?from=search&seid=2343044064008772604
飞机的升力当然越多越好，阻力越小越好。但是地上跑的赛车并不需要升力，向下的力越大越好，也就是Downwash Force，这样就能尽量抵消掉升力。咋办才能减小赛车的升力呢？
第一个方法是把能产生升力的形状反过来安装在后面，这就是赛车后面的扰流板(spoiler)
第二个方法是，由于气流撞到流线型的车身前方后后速度变小，从而压力增大。为了防止高压气流进入车辆下面把车辆抬起来，就在车前方安装前铲(Splitter)。
要减小车辆的阻力，最好使用流线型形的车身。比如下图，上面是公交车的形状，公交车前后面积都很大，导致前面高压区大，车身后低压区或者叫尾流区(wake region)大。这样大量的高压气流会将公交车推向车后的低压区，这是不容忽视的阻力。但如果设计合适形状的车身，可以减小尾流区的面积，进而减小阻力。
当然赛车上的空气动力学不止这些，感兴趣的可看KYLE.ENGINEERS的油管频道 https://www.youtube.com/user/Kyleengineers或EngineersExplained的油管频道https://www.youtube.com/channel/UClqhvGmHcvWL9w3R48t9QXQ。
另外与空气动力学比较相关的就是风力发电机(wind turbine)了。可以阅读《Introduction to Wind Turbine Aerodynamics 》by A. P. Schaffarczyk。各种动物，尤其是鱼类与鸟类上也有很多空气动力学现象。可以阅读《Experimental Hydrodynamics of Fast-Floating Aquatic Animals》 by Babenko Viktor Vitaliiovych。
涡方法

我们可以用NS方程模拟流体，但是用这个方程很难模拟特殊漩涡现象，这时候就要用到流函数涡方法(StreamFunction-Vorticity)或简称涡方法(Vortex Method)。首先再写一遍涡量和速度的关系，以及流函数与速度的关系。u是X方向上的速度，v是Y方向上的速度。

$\omega = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}\qquad u = \frac{\partial \psi}{\partial y} \qquad v = -\frac{\partial \psi}{\partial x}\\$
将两个式子代入，就得到了流函数和涡量的关系。

$\frac{\partial^2 \psi}{\partial x^2} + \frac{\partial^2 \psi}{\partial y^2} + \omega = 0\tag{1}$
上面也是泊松方程，与压力泊松方程类似，不过把压力换成了流函数，散度换成了涡量。其解法也类似，写成python代码如下，代码一共17行。压力泊松方程会在第13章讲SIMPLE算法时详细介绍。
import numpy as np
nmax = 8
psi = np.zeros((nmax,nmax))
v = np.zeros((nmax,nmax))
u = np.zeros((nmax,nmax))
omega = np.zeros((nmax,nmax))
dx = dy = 1
omega[3,3] = 1#涡量为正则逆时针旋转
for k in range(0,100):
for i in range(1,nmax-1):
      for j in range(1,nmax-1):
         term = (psi[i+1,j] + psi[i-1,j])*dy*dy + (psi[i,j+1] + psi[i,j-1])*dx*dx
         psi[i,j] = (term + omega[i,j])/(2*(dx*dx + dy*dy))
for i in range(1,nmax-1):
for j in range(1,nmax-1):
      u[i,j] = (psi[i,j+1] - psi[i,j-1])/(2*dy)
      v[i,j] = (psi[i-1,j] - psi[i+1,j])/(2*dx)打开变量查看器，下图左侧是u，右侧是v。由上到下是X轴，由左到右是Y轴，蓝色是正数，红色是负数。可以看到这些速度正好组成一个逆时针漩涡。这里压图太严重了，自己运行一遍看最好。
很好，涡方法至此我们已经实现了一半了，无论是从内容上还是代码行数上。这时候我们再写一遍涡量输运方程，就是普通NS方程的旋度。

$\frac{\partial \omega}{\partial t} + u\frac{\partial \omega}{\partial x} + v\frac{\partial \omega}{\partial y} = \frac{1}{Re}(\frac{\partial^2 \omega}{\partial x^2} + \frac{\partial^2 \omega}{\partial y^2})\\$
用流函数代替速度：

$\frac{\partial \omega}{\partial t} = - \frac{\partial \psi}{\partial y}\frac{\partial \omega}{\partial x} + \frac{\partial \psi}{\partial x}\frac{\partial \omega}{\partial y} + \frac{1}{Re}(\frac{\partial^2 \omega}{\partial x^2} + \frac{\partial^2 \omega}{\partial y^2})\tag{2}$
这样涡方法的全部两个方程就推导完了，很简单吧。注意流函数和涡量有关系，说明上面是个非线性方程，而非线性方程并不好解，需要用到各种格式。但现在简便起见就先用中心差分了。使用中心差分离散化后如下

$\frac{\omega_{i,j}^{t+1} - \omega_{i,j}}{\Delta t} = - (\frac{\psi_{i,j+1} - \psi_{i,j-1}}{2\Delta y})(\frac{\omega_{i+1,j} - \omega_{i-1,j}}{2\Delta x}) + (\frac{\psi_{i+1,j} - \psi_{i-1,j}}{2\Delta x})(\frac{\omega_{i,j+1} - \omega_{i,j-1}}{2\Delta y})\\ + \frac{1}{Re}(\frac{\omega_{i+1,j} -2\omega_{i,j} + \omega_{i-1,j}}{\Delta x^2} + \frac{\omega_{i,j+1} -2\omega_{i,j} + \omega_{i,j-1}}{\Delta y^2})$
剩下就是处理边界条件了。这里先处理顶盖驱动流(Lid Driven Cavity)的情况。如果左右边界为墙，那么X轴的速度u为零，可得到下面的式子。上下边界同理。

$左右边界为墙 \omega_{wall} = -\frac{\partial^2 \psi}{\partial x^2}\\ 上下边界为墙 \omega_{wall} = -\frac{\partial^2 \psi}{\partial y^2}\\$
流函数边界情况，首先是下边界

$\psi_{i,j=1} = \psi_{i,j=0} + \frac{\psi_{i,j=0}}{\partial y}\Delta y + \frac{\psi^2_{i,j=0}}{\partial y^2}\frac{\Delta y^2}{2} + O(\Delta y^3)\\ \omega_{i,j=0} = -\frac{\psi^2_{i,j=0}}{\partial y^2} \qquad u_{i,j=0} = \frac{\partial \psi_{i,j=0}}{\partial y} = 0\\ \omega_{i,j=0} = -(\psi_{i,j=1} - \psi_{i,j=0})\frac{2}{\Delta y^2}\\$
上边界，顶盖驱动流上边界的u不为零。

$\psi_{i,j=n-2} = \psi_{i,j=n-1} - \frac{\psi_{i,j=n-1}}{\partial y}\Delta y + \frac{\psi^2_{i,j=n-1}}{\partial y^2}\frac{\Delta y^2}{2} + O(\Delta y^3)\\ \omega_{i,j=n-1} = -\frac{\psi^2_{i,j=n-1}}{\partial y^2} \qquad u_{i,j=n-1} = \frac{\partial \psi_{i,j=n-1}}{\partial y} = U_{wall}\\ \omega_{i,j=n-1} = -(\psi_{i,j=n-2} - \psi_{i,j=n-1})\frac{2}{\Delta y^2} - U_{wall}\frac{2}{\Delta y}$
左右边界读者自己推导吧。这样用涡方法写顶盖驱动流的python代码如下，总共34行。
import numpy as np
nmax = 8
psi = np.zeros((nmax,nmax))
v = np.zeros((nmax,nmax))
u = np.zeros((nmax,nmax))
omega = np.zeros((nmax,nmax))
w = np.zeros((nmax,nmax))
dx = dy = dt = 1
uwall = 1
Re = 10
for k1 in range(0,200):
for i in range(1,nmax-1):#用中心差分计算涡量运输方程
      for j in range(1,nmax-1):
         w[i,j] = -(psi[i,j+1] - psi[i,j-1])/(2*dy)*(omega[i+1,j] - omega[i-1,j])/(2*dx)
         w[i,j] += (psi[i+1,j] - psi[i-1,j])/(2*dx)*(omega[i,j+1] - omega[i,j-1])/(2*dy)
         w[i,j] += 1/Re*(omega[i+1,j] -2*omega[i,j] + omega[i-1,j])/(dx*dx)
         w[i,j] += 1/Re*(omega[i,j+1] -2*omega[i,j] + omega[i,j-1])/(dy*dy)
for i in range(1,nmax-1):
      for j in range(1,nmax-1):
         omega[i,j] += dt*w[i,j]
omega[0,:] = -2*psi[1,:]/(dx*dx)#左面
omega[nmax-1,:] = -2*psi[nmax-2,:]/(dx*dx)#右面
omega[:,nmax-1] = -2*(psi[:,nmax-2])/(dy*dy) - uwall*2/dy#上面
omega[:,0] = -2*psi[:,1]/(dy*dy)#下面
for k2 in range(0,100):#泊松方程求解流函数
      for i in range(1,nmax-1):
         for j in range(1,nmax-1):
            term = (psi[i+1,j] + psi[i-1,j])*dy*dy + (psi[i,j+1] + psi[i,j-1])*dx*dx
            psi[i,j] = (term + omega[i,j])/(2*(dx*dx + dy*dy))
for i in range(1,nmax-1):
      for j in range(1,nmax-1):
         u[i,j] = (psi[i,j+1] - psi[i,j-1])/(2*dy)
         v[i,j] = (psi[i-1,j] - psi[i+1,j])/(2*dx)
u[:,nmax-1] = uwall最后是顶盖驱动流的X轴速度u与Y轴速度v的数据，嗯，很完美。
还可以用涡方法模拟2dChannel Flow和后台阶流动，比较折腾的就是边界条件。代码也可在相同仓库找到。
首先是2DChannelFlow，下边界的流函数为零，上边界的流函数为固定值，左边界是流体入口，速度不变，所以流函数为线性，也就是从零增长到上边界的那个固定值。右边界是流体出口则是梯度为零。
#初始条件
for j in range(0,nmax):
psi[0,j] = uwall*j
psi[1:nmax-1,nmax-1] = psi[0,nmax-1]
#每次泊松迭代的边界条件
psi[nmax-1,:] = psi[nmax-2,:]涡流的边界条件是左边界入口的涡量为零，右边界出口涡量梯度为零。后台阶流动差不多。
最后注意，涡量运输方程是非线性方程，而我们所用的中心差分方法不怎么精确，所以如果把时间步调大，或者迭代次数调多，数值最后会爆炸。除了中心差分，还可以选择LaxWendroff和LaxFriedriches或更高阶更精确的格式。但是即使选择了这些格式，仍然不能很好地解非线性问题，至少在2021年是如此。
上面涡方法使用的网格为直角坐标下的均匀网格，但也可以使用极坐标系下的网格，请看参考代码。当然涡方法也可以模拟卡门涡街，但我暂时还没试验成功。
除了这种涡方法，还有一种VortexInCell方法也会用到涡量运输方程。这个方法与ParticleInCell方法相似。而GAMES201中介绍过这种ParticleInCell混合网格粒子方法。
三篇漩涡到此结束。下一篇是泊松压力方程与SIMPLE算法。敬请期待！
上一篇：光影帽子：【游戏流体力学基础及Unity代码（十一）】理想流体机翼绕流和升力原理
参考代码

https://github.com/shohei/simulation 第四章是用涡方法模拟各种算例，包括极坐标系网格下的卡门涡街。
https://github.com/amandaghassaei/VortexShedding webgl模拟卡门涡街
A Finite Difference Code for the Navier-Stokes Equations in Vorticity/Stream Function Formulation  有关涡方法非常棒的讲义！我看了几十篇有关涡方法的介绍文献，都不如这个讲义清楚
Streamfunction-Vorticity Formulation  by A. Salih  Department of Aerospace Engineering
https://curiosityfluids.com/2016/03/14/streamfunction-vorticity-solution-lid-driven-cavity-flow/

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

[笔记] 【游戏流体力学基础及Unity代码（十二）】卡门涡街，边界层，涡方法

本帖子中包含更多资源