什么是涡（vortex），如何识别涡？捕捉涡的方法

什么是涡

涡在自然界中无处不在，大到星云变换，小到分子运动，当然我们最常见的还是龙卷风，台风，倒入茶水时产生的各种各样的涡旋。尤其是在湍流系统中存在着大量大小不一的涡旋，我们通常用尺度来描述这些涡旋的大小。

因为湍流是一个需要能量持续输入的强耗散系统（没有能量输入，湍流是不会持续出现的），由于湍流的这种特性，所以可以很轻易地就想到，湍流中各种各样的涡在其能量耗散中起到了非常重要的作用。这是因为涡的出现意味着在空间上存在速度梯度，在能量耗散率中，速度梯度是其构成部分，因为速度梯度的存在使得流体之间在剪切变形，流体的粘性使得流体的机械能通过类似摩擦的作用变成内能。

这种涡结构在湍流的生成和维持过程中起到了非常重要的作用。通常来讲，涡是流体的旋转。但这种定义是唯象的定义也并不清晰，对于我们如何提取涡结构并无帮助。因此，我们需要一些数学手段，在过去，通常将速度的旋度或者涡量作为涡的数学定义，也是一种提取涡的数学手段，这也被一些科研工作者定义为第一代识别涡的方法。

尽管在各种资料中都有对涡的定义，但我认为对于涡的定义至今仍然是不清晰的，仍然是在一个发展完善的过程中。但涡是流体的旋转仍将是涡定义中关键的一部分。

如何提取涡

第一代识别涡的方法

正如上文所提到的，涡量最先被用于提取涡的结构特征，也就是第一代识别涡的方法，但是这种识别方法有比较大的缺陷，因为涡量并不等于涡，涡量与流体旋转没有很强的关联性。一个非常经典的例子就是在低速平面Poiseuille流动(管道流)中靠近壁面的部分，涡量较大但是完全没有肉眼可见的流体旋转。

涡量 \omega=\frac{\partial{u}}{\partial{x}}-\frac{\partial{v}}{\partial{y}} 在壁面处不为0，却没有明显的旋转。

所以，在壁面处涡量与涡的关联性很弱。显然，用涡量捕捉涡的方法是不合适的。

第二代识别涡的方法

鉴于用涡量提取涡这种方法，有很明显的缺陷，所以产生了新的识别方法，如 \Omega，Q，\Delta，\lambda_{2} 等方法，它们是基于，速度梯度 \frac{\partial{u_{i}}}{\partial{x_{j}}} ，的特征值演化得来的。

\frac{\partial{u_{i}}}{\partial{x_{j}}}=A+B

A=\begin{bmatrix} \frac{\partial u}{\partial x} & \frac{1}{2}(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}) & \frac{1}{2}(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z})\\\\ \frac{1}{2}(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}) & \frac{\partial v}{\partial y} & \frac{1}{2}(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y})\\\\ \frac{1}{2}(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}) & \frac{1}{2}(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}) & \frac{\partial w}{\partial z}\\\\ \end{bmatrix}

B=\begin{bmatrix} 0 & \frac{1}{2}(\frac{\partial u}{\partial y}-\frac{\partial v}{\partial x}) & -\frac{1}{2}(\frac{\partial w}{\partial x}-\frac{\partial u}{\partial z})\\\\ -\frac{1}{2}(\frac{\partial u}{\partial y}-\frac{\partial v}{\partial x}) & 0 & \frac{1}{2}(\frac{\partial v}{\partial z}-\frac{\partial w}{\partial y})\\\\ \frac{1}{2}(\frac{\partial w}{\partial x}-\frac{\partial u}{\partial z}) & -\frac{1}{2}(\frac{\partial v}{\partial z}-\frac{\partial w}{\partial y}) & 0\\\\ \end{bmatrix}

速度梯度张量可以拆分为一个对称张量A和一个反对称张量B，其中对称张量A表示流体的变形效果，反对称张量B表示流体的旋转。

例如 \Omega，Q 两种识别方法，都是要求在涡结构中B要存在，也就是指涡量要存在，并且B的旋转效应要大于A的变形效应。

\Delta，\lambda_{2} 则是由另外的基于特征值的理论推导得到，下面将简要介绍这四种方法的具体内容和推导。

Q方法：

|\frac{\partial{u_{i}}}{\partial{x_{j}}}-\lambda I|=\begin{bmatrix} \frac{\partial u}{\partial x}-\lambda & \frac{\partial u}{\partial y} & \frac{\partial u}{\partial z}\\\\ \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y}-\lambda & \frac{\partial v}{\partial z}\\\\ \frac{\partial w}{\partial x} & \frac{\partial w}{\partial y} & \frac{\partial w}{\partial z}-\lambda\\\\ \end{bmatrix}

也就是 \lambda^{3}+P\lambda^{2}+Q\lambda+R=0

其中

P=-(\lambda_{1}+\lambda_{2}+\lambda_{3})

Q=\lambda_{1}\lambda_{2}+\lambda_{2}\lambda_{3}+\lambda_{1}\lambda_{3}

R=-\lambda_{1}\lambda_{2}\lambda_{3}

而且Q可以用A和B重写为

Q=\frac{1}{2}(||B||_{F}^{2}-||A||_{F}^{2})

其中Frobenius范数是该矩阵各个元素的平方加和再开方。

所以直观上来讲，要求Q大于0就是要求在涡结构中B要存在，也就是指涡量要存在，并且B的旋转效应要大于A的变形效应。具体这个阈值如何选择，根据实际情况选择。

\Omega 方法：

\Omega 方法一定程度上是Q的扩展，在特殊情况下，可以退化为Q方法。

从上面Q的公式可以看出，他是B的旋转效应要大于A的变形效应，仅仅是一个差值，具体涡结构的提取依赖于阈值的选择，这样会存在一个问题，假设B与A的F范数较小，且B大于A，但由于阈值的限制，有些细小的涡结构会被忽略掉。

所以就需要重新定义一下，以避免这个问题。

所以 \Omega 不是以差值来寻找涡结构，而是以B的旋转效应在整个速度梯度张量中所占的比重来判断涡结构。

公式如下：

\Omega=\frac{||B||_{F}^{2}}{||B||_{F}^{2}+||A||_{F}^{2}}

为了避免分母较小时，带来较大的误差，可以在分母上加上一个小的正数。

\Omega=\frac{||B||_{F}^{2}}{||B||_{F}^{2}+||A||_{F}^{2}+\varepsilon}

自然而然 \Omega 的取值是0-1，基于Q，我们可以选择 \Omega>0.5 作为阈值，具体阈值如何选择，读者可以参考相应的文献。

\Delta 方法：

该方法来源于A E Perry, M S Chong, 1987&1990 所提出的临界点理论，在他的临界点理论中，用特征值来辨别和分离流动模式。

对于一个矩阵而言，可以看作是对一个向量进行旋转和放缩。

\frac{\partial{u_{i}}}{\partial{x_{j}}} 实际上是对一个带有方向的流体质点产生作用（我个人觉得这里不如理解为一小段带有方向的流线）

实特征值，代表对该向量的放大，不涉及旋转。

复特征值，带表对该向量具有旋转作用。

所以对于三维流动，要求至少有两个复特征值（复特征值共轭），可以认为这里流线发生了弯曲。

因此我们要找到一个量来判断 \frac{\partial{u_{i}}}{\partial{x_{j}}} 是否有两个复特征值。

对于QPR满足如下公式：

27R^{2}+(4P^{3}-18PQ)R+(4Q^{3}-P^{2}Q^{2})=0

整理得到，若要求有两个复根，则 \Delta>0 （具体可以参考上面两篇文章）

\Delta=(Q'/3)^{3}+(R'/2)^{2}

Q'=Q-P^{2}/3

R'=R+2P^{3}/27-PQ/3

对于不可压缩流动，P=0

所以Q>0时，也满足 \Delta>0 ，也就是说Q方法也是它的子集。而且该方法有个优势，由于它是严格的数学推导，对于涡流区和非涡流区的判断是非常严谨和准确的。

（2023.02.09更新）

\lambda_{ci} 方法：

\lambda_{ci} 是 \Delta 方法的扩展，按照 \Delta 的正负可以判断出涡区和非涡区，但其实并不能很好的判断出涡的强度。所以需要对复数特征值进一步考察。

按照 \Delta ，我们有两个复数特征值，记为 \lambda=\lambda_{cr}\pm\lambda_{ci}

按照线性代数的理论，对于复数特征值，虚数部分代表旋转的角度大小。因此用其虚数部分来表示涡的强度。

第三代识别涡的方法

第三代涡识别方法以正在发展和完善的Liutex方法为主。

Liutex 涡定义：

待续……