一种流场旋涡涡核核点提取和连接方法

技术领域

本申请涉及流体力学领域，具体为一种流场旋涡涡核核点提取和连接方法。更具体地，本申请提供一种基于欧拉框架下的流体旋涡涡核识别和涡核线连接方法。

背景技术

在飞行过程中，飞机周围空气受到强烈扰动，会在机翼后方形成复杂的尾流场。飞机尾流由一系列复杂的涡系组成，包括机翼、舵面、增升部件产生的翼尖涡、翼面后缘空域产生的边界层湍流，以及发动机的喷流、机身扰流等。在进一步的演化中，涡系卷起融合为一对反向旋转的涡对。大型飞机的尾流持续时间通常较长，在不同的飞机型号和尺寸、飞行速度、高度、气象条件下，尾流持续时间通常在几分钟到十几分钟不等。飞机尾流的影响问题与工程应用紧密相关。

对于起降间隔问题，进入前机尾流区域的后机而言，尾流产生的横向滚转力矩可能超过所受影响航空器的滚转操作限制，使进入后机发生俯仰、翻滚、机身剧烈抖震、飞机引擎喘振、翼尖过载突变等严重状况，甚至会使后机发生失速等危险。在编队飞行中，后机可以在前飞机尾流的上洗气流区域飞行。在定升力飞行情况下，这意味着诱导阻力减小，进而使得燃油消耗减少。空降期间，飞机构型、降落伞系统和空降策略旨在最大限度保证伞降的安全性，空降者本身和飞机周围的不稳定流动(舱门打开)之间总是存在一些相互作用，特别是在空投的最初几秒内，完全由飞行器的尾流场所决定。因此，对于近场尾流场的精确预测对于空投系统的轨迹预测非常重要。进一步，对流场中旋涡的提取和分析变得尤为重要。

传统的涡核提取方法存在如下两方面的缺陷：

(1)由于旋涡动力学特征较为模糊，包括衰减速度法、Lambda2最小值法、压力最小值法等在内的传统涡核提取方法，旋涡识别准则就没有完全分离剪切污染，最终提取出的涡核线不能较好表征旋涡的涡核；

(2)传统的涡核线提取方法需要给定一个种子点，或者采用种子点生成方法生成一系列种子点，需要研究者人为给定一个可能存在涡核点的范围；但是对于未知全貌的流场，就需要不断尝试以获得完整的涡核线。

为此，迫切需要一种不需要预先知晓流场中旋涡的大致分布，就可以提取出流场中所有涡核点，并连接成涡核线的方法。

发明内容

本申请的发明目的在于：针对上述存在的问题，提供一种流场旋涡涡核核点提取和连接方法。本发明提出一种涡核点的定义，通过在每个网格面上遍历获得满足定义的涡核点，再通过连接算法将其连接起来。本申请不需要预先知晓流场中旋涡的大致分布，就能够提取出流场中所有涡核点，并连接成涡核线。同时，采用本申请的方法，涡核点定义足够准确，能够解决传统涡核提取方法存在的两方面问题。

为了实现上述目的，本申请采用如下技术方案：

一种流场旋涡涡核核点提取和连接方法，包括如下步骤：

S1、对于已有流场数据，根据速度场计算出后续用于识别旋涡与涡核的场量；

S2、将涡核点定义转化为涡识别场量；

S2.1计算旋转强度的梯度

S2.2根据涡核点的定义，在涡核点处，满足以下公式(14)：

因此，所有涡核点可以表征为以下涡核点集S，如下公式(15)所示：

公式(15)中，p代表点的空间坐标，

S2.3根据平行算子方法，涡核点集S可以转化为以下定义，如公式(16)所示：

S3、在网格面上迭代的雅可比矩阵J

S3.1、将公式(16)中的定义转化函数c，采用如下公式(17)表示：

S3.2、假设空间中的网格点分布为结构化网格点，其节点沿着空间三个方向X,Y,Z的索引分别为I,J,K，那么在某一切面K＝k上，存在由四个网格点(i,j,k)、(i+1,j,k)、(i,j+1,k)、(i+1,j+1,k)构成的一个网格面，(i,j,k)、(i+1,j,k)、(i,j+1,k)、(i+1,j+1,k)、(i,j,k+1)、(i+1,j,k+1)、(i,j+1,k+1)、(i+1,j+1,k+1)八个网格点构成一个网格单元；

记沿着网格面索引方向上的偏导数为ξ

公式(18)中，

S4、在网格面上采用牛顿迭代法，求解满足公式(16)的解；

S4.1、对于已有的函数f(x)，采用牛顿迭代法有：

x

公式(19)中的x

J

公式(20)中，Δξ为两个空间坐标ξ

Δξ＝ξ

S4.2、采用最小二乘法计算Δξ；

以网格面中间点ξ

ξ

其中，Δξ

S4.3、判断解的有效性；

对求解出来的距离做出判断，确定涡核点是否在网格面上，即：

ξ

公式(24)中，S(ξ

如果流场的网格为均匀网格，且索引方向I,J,K分别与空间三个方向X,Y,Z平行，公式(24)的表达式转化为：

公式(25)中，p

如果ξ

S5、确定涡核点的连接准则；

假设涡核点集S中的任意两个点P

插值获得P

计算r

计算

若D小于相邻两个网格单元中心之间的距离且Δ小于45°，则满足连接条件；

S6、获得涡核点集S；

根据步骤S4，对流场中的所有网格面遍历，获得整个流场的涡核点集S；

S7、将S6获得的涡核点集连接成线段；

记涡核点集合为S，涡核线集合为L；

S7.1遍历流场中的每一个网格单元，具体操作如下：

(a)当网格单元中的多个网格面上只存在1个涡核点，则继续遍历；

(b)当网格单元中的多个网格面上存在2个涡核点，则将两个涡核点进行连接，将连接的涡核线归入到涡核线集合L中，并从涡核点集S中删除这两个涡核点；

(c)当网格单元中的多个网格面上存在2个以上的涡核点，则执行以下操作：

d)根据排列组合方式对可能的连接方式进行排列，如网格单元上存在3个涡核点，则一共有

e)计算每种排列组合方式的连接方向矢量

f)计算涡核点的平均旋转轴方向

g)计算每种排列方向矢量V

h)选择最大的乘积之模max(||V

S7.2、将不同网格单元上的同一个网格面上的涡核点进行连接；

S7.2.1、将结构化面索引按以下公式转化为单调增函数的一维数值索引集合

公式(25)中，exp代表指数函数，cos代表余弦函数；i

S7.2.2将根据一维数值索引集合n递增顺序排列对应的涡核点集合S，记为S

S7.2.3遍历重新排序的涡核点集合S

S7.3对涡核点集S以及涡核线集合L进行临近搜索连接；

S7.3.0令h的初始值＝1；

S7.3.1令i的初始值＝1，对于涡核线集合L执行以下操作：

选择涡核线集合中的第h条涡核线的第一个点P

对涡核点集S，执行以下操作：

N1.对于涡核点集的第i个点，记为P

N2.如果两个条件均满足或其中一个都满足，则将对应的涡核点进行连接，将连接后的新涡核线更新到涡核线集合L中，并从涡核点集S中删除这些涡核点，对新的循环重新执行步骤S7.3的操作；

N3.如果两个条件均不满足，则i＝i+1；

S7.3.2在步骤S7.3.1中，当i＞涡核点集S中涡核点总数时，h＝h+1，再返回步骤S7.3.1；

S7.3.3在步骤S7.3.1中，当h＞涡核线集合L中的涡核线总数时，迭代中止。

已有流场数据包括x,y,z,u,v,w；其中，x,y,z分别对应空间相互正交方向上的三个坐标，u,v,w分别对应x,y,z三个方向上的速度；

所述步骤S1包括如下步骤：

S1.1计算速度梯度张量与涡量，计算公式如下公式(1)、公式(2)所示：

式中，U为速度向量，包含速度场的三个分量u,v,w，

式中，

S1.2根据获得的速度梯度张量,求解其特征方程；

S1.2.1列出特征方程表达式，如下公式(3)所示：

式(3)中，λ是待求解的特征值,I代表单位对角矩阵，将公式(3)写成公式(4)的形式：

λ

公式(4)的方程中，设解为λ

式(5)-式(7)中，tr表示矩阵的迹，det表示矩阵的行列式；

S1.2.2根据获得的P,Q,R，计算特征方程的判别式，如下公式(8)所示：

公式(8)中，Δ为特征方程的判别式，

S1.3当特征方程的判别式Δ>0，根据临界点理论，可知该区域有旋，求解公式(4)的特征根，此时得到一个实特征根λ

公式(11)中，λ

S1.4根据以上步骤获得的涡量ω、局部旋转轴r以及复数特征根的实数部分λ

式(12)中，·表示点积。

附图说明

图1为网格面示意图。

图2为网格单元示意图。

图3为本申请的方法应用到梯形翼旋涡提取的效果图。

图4本发明的方法与传统方法的对比图。

具体实施方式

本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

本说明书中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中的一个例子而已。

实施例1

本实施例提供一种流场旋涡涡核核点提取和连接方法，其包括如下步骤。

S1、对于已有流场数据，根据速度场计算出后续用于识别旋涡与涡核的场量。

已有流场数据包括x,y,z,u,v,w；其中，x,y,z分别对应空间相互正交方向上的三个坐标，u,v,w分别对应x,y,z三个方向上的速度。

所述步骤S1包括如下步骤：

S1.1计算速度梯度张量与涡量，计算公式如下公式(1)、公式(2)所示：

式中，U为速度向量，包含速度场的三个分量u,v,w，

式中，

S1.2根据获得的速度梯度张量,求解其特征方程，具体如下。

S1.2.1列出特征方程表达式，如下公式(3)所示：

式(3)中，λ是待求解的特征值,I代表单位对角矩阵。

将公式(3)写成公式(4)的形式：

λ

公式(4)的方程中，设解为λ

式(5)-式(7)中，tr表示矩阵的迹，det表示矩阵的行列式。

S1.2.2根据获得的P,Q,R，计算特征方程的判别式，如下公式(8)所示：

公式(8)中，Δ为特征方程的判别式，

S1.3当特征方程的判别式Δ>0，根据临界点理论，可知该区域有旋，求解公式(4)的特征根，此时得到一个实特征根λ

公式(11)中，λ

S1.4根据以上步骤获得的涡量ω、局部旋转轴r以及复数特征根的实数部分λ

式(12)中，·表示点积。

S2、将涡核点定义转化为涡识别场量，具体操作如下。

S2.1计算旋转强度的梯度

S2.2根据涡核点的定义，在涡核点处，满足以下公式(14)：

因此，所有涡核点可以表征为以下涡核点集S，如下公式(15)所示：

公式(15)中，p代表点的空间坐标，

S2.3根据平行算子方法，涡核点集S可以转化为以下定义，如公式(16)所示：

S3、在网格面上迭代的雅可比矩阵J

S3.1、将公式(16)中的定义转化函数c，采用如下公式(17)表示：

S3.2、假设空间中的网格点分布为结构化网格点，其节点沿着空间三个方向X,Y,Z的索引分别为I,J,K，那么在某一切面K＝k上，存在由四个网格点(i,j,k)、(i+1,j,k)、(i,j+1,k)、(i+1,j+1,k)构成的一个网格面，(i,j,k)、(i+1,j,k)、(i,j+1,k)、(i+1,j+1,k)、(i,j,k+1)、(i+1,j,k+1)、(i,j+1,k+1)、(i+1,j+1,k+1)八个网格点构成一个网格单元，如图1所示。

记沿着网格面索引方向上的偏导数为ξ

公式(18)中，

S4、在网格面上采用牛顿迭代法，求解满足公式(16)的解。

S4.1、对于已有的函数f(x)，采用牛顿迭代法有：

x

公式(19)中的x

J

公式(20)中，Δξ为两个空间坐标ξ

Δξ＝ξ

S4.2、采用最小二乘法计算Δξ。

由于流场中的节点之间直接采用线性插值，因此采用牛顿迭代并不能提高精度，因此在求解满足公式(16)解的涡核点时，直接采用最小二乘法获得结果。此时，以网格面中间点ξ

ξ

其中，Δξ

S4.3、判断解的有效性。

对求解出来的距离做出判断，确定涡核点是否在网格面上，即：

ξ

公式(24)中，S(ξ

公式(25)中，p

S5、确定涡核点的连接准则。

假设涡核点集S中的任意两个点P

插值获得P

计算r

计算

若D小于相邻两个网格单元中心之间的距离且Δ小于45°，则满足连接条件。

S6、获得涡核点集S。

根据步骤S4，对流场中的所有网格面遍历，获得整个流场的涡核点集S。

记录网格单元所在网格单元索引以及网格单元面的索引。如图2以结构化网格六面体为例，以I′,J′,K′记录其在空间上三个方向的索引，如图2所示，假设网格单元为结构化索引为(i,j,k)，则网格单元的六个面索引分别为：S1、S2、S3、S4、S5、S6。为降低数据存储量，只记录网格单元结构化索引以及一连串0/1代表是否有涡核点，同时记录对应网格面涡核点所在的位置。典型例子，如表1所示。

表1

表1中，(i,j,k)代表索引为(i,j,k)的网格单元上，0,1,0,0,0,1分别代表S1(i-1/2,j,k)、S2(i+1/2,j,k)、S3(i,j-1/2,k)、S4(i,j+1/2,k)、S5(i,j,k-1/2)、S6(i,j,k+1/2)。[x

S7、将S6获得的涡核点集连接成线段。

记涡核点集合为S，涡核线集合为L。

S7.1遍历流场中的每一个网格单元，具体操作如下。

(a)当网格单元中的多个网格面上只存在1个涡核点，则继续遍历。

(b)当网格单元中的多个网格面上存在2个涡核点，则将两个涡核点进行连接，将连接的涡核线归入到涡核线集合L中，并从涡核点集S中删除这两个涡核点。

(c)当网格单元中的多个网格面上存在2个以上的涡核点，则执行以下操作：

d)根据排列组合方式对可能的连接方式进行排列，如网格单元上存在3个涡核点，则一共有

e)计算每种排列组合方式的连接方向矢量

f)计算涡核点的平均旋转轴方向

g)计算每种排列方向矢量V

h)选择最大的乘积之模max(||V

S7.2、将不同网格单元上的同一个网格面上的涡核点进行连接。

S7.2.1、将结构化面索引按以下公式转化为单调增函数的一维数值索引集合

公式(25)中，exp代表指数函数，cos代表余弦函数；i

S7.2.2将根据一维数值索引集合n递增顺序排列对应的涡核点集合S，记为S

S7.2.3遍历重新排序的涡核点集合S

记第i个和第i+1个涡核点为P

S7.3对涡核点集S以及涡核线集合L进行临近搜索连接。

S7.3.0令h的初始值＝1。

S7.3.1令i的初始值＝1，对于涡核线集合L执行以下操作：

选择涡核线集合中的第h条涡核线的第一个点P

对涡核点集S，执行以下操作：

N1.对于涡核点集的第i个点，记为P

N2.如果两个条件均满足或其中一个都满足，则将对应的涡核点进行连接，将连接后的新涡核线更新到涡核线集合L中，并从涡核点集S中删除这些涡核点，对新的循环重新执行步骤S7.3的操作；

N3.如果两个条件均不满足，则i＝i+1。

S7.3.2在步骤S7.3.1中，当i＞涡核点集S中涡核点总数时，h＝h+1，再返回步骤S7.3.1。

S7.3.3在步骤S7.3.1中，当h＞涡核线集合L中的涡核线总数时，迭代中止。

流体力学中，流场通常包含数百万个以上的数据，为便于理解，选取十个节点作为典型进行说明，节点数据如下表2所示。

表2节点及三个方向速度的数据

根据步骤s1，获得识别旋涡和涡核的场量，如下表3所示。

表3根据三个方向速度获得涡识别量数据

根据步骤S2，获得识别旋涡和涡核的场量，如下表4所示。

表4

根据步骤S3进行计算，结果如下表5所示。

表5

根据步骤S6，基于整个流场(包含数百万个以上的数据)，得到涡核点所在的网格单元、网格面以及涡核位置，如下表6所示。

表6

根据步骤S7，获得涡格线集，一共787条涡核线，如表7所示。

表7

/>

将本申请的方法与现有技术对比，相关结果如下。

图3给出了该方法应用到梯形翼旋涡提取的效果。其中，图3左侧图是遍历整个流场获得的涡核点集合，图3右侧图是采用发明提出的线连接算法连接而成的涡核线。该方法将梯形翼的翼尖涡、襟翼涡、机身涡提取出来。

图4给出了本申请提出的涡核点定义与著名后处理软件Tecplot之间的差异。从图4中可以看出：如图4左侧图所示的旋涡中，采用Tecplot内置算法提出的涡核点缺少连接细节(如图4中间图所示)，而采用本发明提出算法提取的涡核点则将旋涡连接细节也提取出来了(如图4右侧图所示)。

本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。