【游戏流体力学基础及Unity代码（十）】漩涡和模拟二维烟雾

123456833

接下来几篇都和漩涡有关。
涡量，旋度和环量

物体会旋转是由于两端受力不均匀导致的。比如固定住木条的一端，而给另一端一定的速度，那么这个木条肯定会旋转起来。假如我们把一个正方形流体微团放在水流之中，它什么时候会旋转呢？如下图左边那个橙色正方形，流体微团左侧上部分的x轴速度为2，左侧下部分受到的x轴速度为1。流体微团下侧左部分y轴方向速度为2，下侧右部分y轴方向速度为3。这个正方形是不会旋转的。如果你没想明白，可以把同时把y轴方向速度减去一，就变成了下图右侧的情况。
虽然上面这个流体微团在xy轴速度都不一致，运动后甚至会线变形，但却并没有旋转，因为这个流体微团在xy方向的速度梯度都是一样的，为1。我们可以用一个公式来计算流体微团是否会旋转，得到的结果就是涡量(Vorticity)，用omega表示。二维涡量的计算公式如下。

这样涡量恰好是角速度的两倍。注意角速度的符号也是omega，因此有些书籍使用xi来表示涡量。但在本系列文章中，仍然使用omega表示涡量。来一道思考题，如果涡量为正，那么流体微团会顺时针旋转还是逆时针旋转？
仅仅旋转流体微团怎么够呢？于是地球自己玩了票大的。地球在由西向东旋转的时候，每一点的角速度都是一样的，但是赤道所在的纬度线比较长，而靠近南北极的地方的纬度线短，所以赤道上的线速度要比南北极快，所以自然产生了大漩涡，导致南北半球气旋的旋转方向不一样，从太空中看北半球逆时针，南半球顺时针。
如果已知了一个三维速度场，要计算涡量，就要计算速度场的旋度(curl)，如下：


上面是只有一个流体微团的情况。如果有很多小流体微团组成了一个大的区域，将它们的涡量相加，然后取负号，就得到了环量(Circulation)。用Gamma表示

这样当流体微团都在顺时针旋转时，它们的总环量就是正的。然而很多时候，计算每个小流体微团的涡量太麻烦，能不能只在大流体微团的边界计算就算出环量？答案是肯定的。这个方法与第四章的散度定理的本质一样，我们算每个小流体微团的的速度梯度再相加，与直接计算大流体微团的速度积分再相加，其结果是一样的。如下图左右的环量是一样的，假设X轴方向速度为y + 1，而Y轴方向速度为2x + 1。
那么计算公式如下

嗯，又是一个小学二年级的知识【误】。这就是旋度定理，又称斯托克斯定理(Stokes' theorem)，用公式表示如下。

上式V是区域C边界处的速度，l是边界。等式左边就是边界的速度梯度。omega是涡量，n是平面法向量，S是整个区域。旋度定理在高等数学中又叫做格林公式。格林公式如下，你只要把P看成是X轴方向的速度，把Q看成是Y轴方向的速度，就可以在格林公式和旋度定理间反复横跳。

最终环量的计算方法如下

有时候流体微团看起来像是在旋转但其实涡量为零，并没有旋转。比如你沿着400米一圈的塑胶跑道跑步，那么你在沿着跑道前进的时候，还要将自己的目光对准跑道的切线方向。一边调整自己身体的旋转角度，某些时候面朝北，某些时候面朝南。如下图中间那样，这是正常的跑法，相当于下图中间的有旋流(rotational flow)。但如果你不肯调整自己的旋转角度，无论跑到哪里，都只朝一个方法比如北方，那么你自身并未旋转，相当于下图右的无旋流(irrotational flow)。
Kelvin环量定理

假如舞台上有一些演员正在做旋转动作，演员A每分钟旋转10次，演员B每分钟旋转8次，并且他们的位置不发生变换。如果用记号笔在舞台的地面上画一条线将演员A圈起来，那么线内所有演员的旋转次数就是10次，如下图左面那样，橙色点是演员，蓝色线就是笔画的线，那么线内所有演员的环量就是360度乘以10。如果画的这条线没有将任何演员圈起来，那么旋转次数为零，总旋转次数也为零。这个和复变函数中的路径积分很相似。
圈内总旋转次数可以通过将各个演员分别的旋转次数相加起来进行计算，也可以通过路径积分来算。这是之前说到的旋度定理。但如果路径复杂导致路径积分太过复杂，我们可以将复杂的路径分解成围绕演员的简单路径，也就是围绕奇点的简单路径，所得到的环量是相同的，这就是复变函数中的留数定理。和旋度定理有异曲同工之妙？
演员可能移动，你也可以任意改变线的形状，但只要保证不把新的演员圈进来，或者把让原来的演员移到线圈外面去，那么圈里所有演员的总旋转次数也不会发生变化。但如果把演员B也圈起来，那么圈里所有演员的总旋转次数就会变化了，不再是10而是变成了18。
把上面的演员看成流体微团，把单个演员旋转次数看成涡量，把总旋转次数换成环量。就得到了流体力学中的一个定理，即在理想流体或者无粘性流体中，环量的物质导数沿着围线C不会发生改变。这是因为在无粘性流体中，没有剪切力的存在，流体微团的角速度不可能变化。这就是即Kelvin's circulation theorem。


涡量输运方程

当我们把一堆小木块丢入河流中，木块会由于对流和扩散而改变自己的位置，但总木块数量不会改变。涡量也是如此，对流和扩散作用仅仅是改变涡量的分布，不会改变涡量的总大小。如下的动量方程，最后那个q的散度表示外部体积力。对这个方程取旋度

就得到了涡量运输方程


Helmholtz漩涡定理

根据Kelvin环量定理，如果没有粘性作用，无旋的流体仍然会一直是无旋的。
漩涡由于对流运动的时候，漩涡走过的路径就是涡线(Vortex Line)。如果有粘性扩散效果，速度梯度会逐渐减小，大漩涡就会被分解成小漩涡。但如果在理想无粘性流体中，大漩涡仍能保持住自己的形状。这样顺着漩涡走过的路径，它的涡量也不会变。
由多条通电导线组合成一起，就成了电缆。同样让多条涡线组合在一起，就组成了涡管(Vortex tube)。由于牛顿第一定律，如果没有粘性所带来的剪切力，那么涡量就不会变化，漩涡就不会产生或消失，涡管在流体中就必须一直延续下去，直到形成闭合的曲线或碰到边界。如果涡管的横截面积无限小，就成了Vortex filament。
如果楼下烧烤摊卖的烧烤串也能不花钱一直延续下去就好了。
涡线的形状类似上图那样。涡线的准确定义是，涡线的切向方向与涡量的方向一致。如果一个速度场仅在xy方向有速度，z速度方法为零，那么涡量则在xy方向上为零，z方法上不为零，因此这时候的涡线是沿着z轴的。
Bio-Savart Law

水箱里原来有一些静止的水，可以伸手去搅动它让水旋转起来。搅动导致了涡量的产生，继而让水有了速度。如下图。
更准确地说，在一条vortex filament上，由单位环量Gamma产生在某一处的速度计算公式如下。

如果这个涡管无限长直，那么速度可以如下表示。

这就是Bio-Savart Law，我们更经常在电磁学中见到它。如果你想知道无限长直涡管产生的速度是怎么算的，推荐看看格里菲斯的《电动力学导论》。电磁学中的Bio-Savart Law表述如下，在一根通电长直导线上，由单位导线dl产生在某一点的磁通密度B计算公式公式如下。

上式B是磁场强度，mu0是真空磁导率，I是电流强度，r是由导线到空间某一点的距离。如果这个导线是无限长的话，那么磁通强度可以用下式表示

unity模拟二维烟雾
在烟雾中，通常有大大小小的漩涡和湍流。然而精度不高的数值算法会造成人工粘性，让漩涡在模拟中难以看见，就像在本系列第五篇那样。总涡量并没有减少，只是由于粘性扩散将大涡量分解成很小的涡量了。于是我们可以找到有这些小涡量的地方，将原来因为数值误差造成的能量耗散补充回来，也就是施加力的影响，将这些小漩涡重新变成清晰可见的大漩涡。这就是Vorticity Confinement。
第一步计算涡量，这样在应该有漩涡的地方，涡量就不为零。

第二步是归一化涡量位置向量，让它从低涡量的地方指向高涡量的地方

第三步，再让这个涡量方向乘上涡量本身，然后再乘上一个控制参数控制涡量放大倍数，就得到了为了形成漩涡需要额外添加的力。最后将这个力乘上时间，就可以加到速度中去了。

在c#代码中，为了让漩涡更加清晰，所以把粘性扩散着色器去掉了。然后处理好鼠标拖动的影响后，添加如下两个调用着色器的代码
//第四步：计算Curl
CurlMat.SetTexture("VelocityTex", VelocityRT2);
Graphics.Blit(VelocityRT2, CurlRT, CurlMat);

//第五步：计算旋度，更新速度，得到有散度的速度场
VorticityMat.SetTexture("VelocityTex", VelocityRT2);
VorticityMat.SetTexture("CurlTex", CurlRT);
VorticityMat.SetFloat("curl", Vorticity);
VorticityMat.SetFloat("dt", dt);
Graphics.Blit(VelocityRT2, VelocityRT, VorticityMat);
Graphics.Blit(VelocityRT, VelocityRT2);
对于第一个Curl.Shader，旋度或者说是涡量的计算方法如下
float4 frag(v2f i) :SV_Target{
  float L = tex2D(VelocityTex,i.uv - float2(VelocityTex_TexelSize.x,0.0f)).y;
  float R = tex2D(VelocityTex,i.uv + float2(VelocityTex_TexelSize.x, 0.0f)).y;
  float T = tex2D(VelocityTex, i.uv + float2(0.0f, VelocityTex_TexelSize.y)).x;
  float B = tex2D(VelocityTex, i.uv - float2(0.0f, VelocityTex_TexelSize.y)).x;
  float vorticity = R - L - T + B;
  return float4(vorticity,0.,0.,1.);
}
对于第二个Vorticity.Shader，代码如下
float4 frag(v2f i) :SV_Target{
  float L = tex2D(CurlTex,i.uv - float2(CurlTex_TexelSize.x,0.0f)).y;
  float R = tex2D(CurlTex,i.uv + float2(CurlTex_TexelSize.x, 0.0f)).y;
  float T = tex2D(CurlTex, i.uv + float2(0.0f, CurlTex_TexelSize.y)).x;
  float B = tex2D(CurlTex, i.uv - float2(0.0f, CurlTex_TexelSize.y)).x;

  float C = tex2D(CurlTex,i.uv).x;

  float2 force = float2(abs(T)-abs(B),abs(R)-abs(L));
  force *= 1./length(force + 0.00001) * curl * C;//为了防止除零
  float2 vel = tex2D(VelocityTex, i.uv).xy;
  return float4(vel + force * dt,0.0f,1.0f);
}
除了漩涡效果外，密度和温度都会影响流体的速度。较重的烟雾颗粒会逐渐下沉，而热空气则会上升。可以使用如下公式，其中z = (0,0,1)，也就是竖直向上的单位向量。公式：


代码地址：
上一篇：光影帽子：【游戏流体力学基础及Unity代码（九）】用浅水波方程模拟雨落池塘和DamBreak
下一篇：光影帽子：【游戏流体力学基础及Unity代码（十一）】理想流体机翼绕流和升力原理
参考

https://github.com/TNWX-Z/EnhanceSmokeSimulationPro
《Fundamentals of Aerodynamics》 by John D. Anderson
《Elements of Vorticity Aerodynamics 》by James C. Wu
Modification of the Euler equations for "vorticity confinement": Application to the computation of interacting vortex rings
Visual Simulation of Smoke