CN111414720B - 一种基于神经网络的流场旋涡检测方法 - Google Patents

一种基于神经网络的流场旋涡检测方法 Download PDF

Info

Publication number
CN111414720B
CN111414720B CN202010097688.5A CN202010097688A CN111414720B CN 111414720 B CN111414720 B CN 111414720B CN 202010097688 A CN202010097688 A CN 202010097688A CN 111414720 B CN111414720 B CN 111414720B
Authority
CN
China
Prior art keywords
vorticity
grid
field
grid points
flow field
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010097688.5A
Other languages
English (en)
Other versions
CN111414720A (zh
Inventor
王岳青
邓亮
王昉
赵丹
杨志供
刘杨
喻杰
杨文祥
陈呈
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Original Assignee
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Computational Aerodynamics Institute of China Aerodynamics Research and Development Center filed Critical Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority to CN202010097688.5A priority Critical patent/CN111414720B/zh
Publication of CN111414720A publication Critical patent/CN111414720A/zh
Application granted granted Critical
Publication of CN111414720B publication Critical patent/CN111414720B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06NCOMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00Computing arrangements based on biological models
    • G06N3/02Neural networks
    • G06N3/04Architecture, e.g. interconnection topology
    • G06N3/045Combinations of networks
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06NCOMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00Computing arrangements based on biological models
    • G06N3/02Neural networks
    • G06N3/08Learning methods
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00Geometric image transformations in the plane of the image
    • G06T3/60Rotation of whole images or parts thereof
    • G06T3/604Rotation of whole images or parts thereof using coordinate rotation digital computer [CORDIC] devices

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • Biomedical Technology (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Health & Medical Sciences (AREA)
  • Image Analysis (AREA)
  • Complex Calculations (AREA)

Abstract

本发明一种基于神经网络的流场旋涡检测方法,先对流场内所有物理网格点打标签,标记其是否属于旋涡区域,得到标签数据,对涡量场进行网格转换,将物理网格中的网格点一一映射到计算网格的网格点,抛弃物理坐标信息,得到计算网格下的涡量场。根据计算网格下的的涡量场,将每个网格上的涡量值减去涡量场的均值后除以标准差,得到归一化后的涡量场;对得到的归一化后的涡量场和标签数据同时随机采样,采样后对区域内的网格点对应的标签进行判断,构建神经网络,使用随机采样得到的带标签局部涡量场数据训练该神经网络,固定神经网络参数,得到固定参数的神经网络;利用第五步得到的固定参数的神经网络,对待检测涡区的流场进行涡区检测。

Description

一种基于神经网络的流场旋涡检测方法
技术领域
本发明涉及一种基于神经网络的流场旋涡检测方法,属于流场物体检测技术领域。
背景技术
旋涡是流体在流动过程产生的一种螺旋状形态,常见的旋涡包括飞机机翼后的涡、***等。旋涡在工程中发挥重要的作用,在流体力学界被认为是流体运动中的肌腱,飞机飞行时,旋涡往往产生在机翼后方,旋涡的发生会减少飞机的升力,造成飞机失速,产生消极作用;在流体混合的过程中,旋涡能够加速流体的混合,发挥着积极的作用;旋涡还与气动噪声的产生和消除有着密切的联系。计算流体力学(CFD)通过计算机和数值方法来求解流体力学的控制方程,对流体力学问题进行模拟和分析,CFD是迭代计算的,因此一次计算会产生数以万计的流场。流场是指运动流体所占有的空间区域,包括速度场、压力场、密度场等多个不同的物理变量场。流场往往通过网格点划分,本发明中所说的流场内的点就是指流场内的网格点,CFD计算就是通过对初始场进行计算,不断迭代,产生多个中间流场,模拟出流场的变化。CFD一次计算会产生海量的流场数据,流场内产生旋涡的位置很重要,因此需要对所有流场数据进行判断其是否包含旋涡。
目前尚无现有技术针对流场进行涡结构检测,现有技术只有针对图像领域的物体检测,由于流场与图像完全不同,因此无法直接适用于流场旋涡检测。
关于现有技术中已有的涡检测方法包括全局检测方法、局部检测方法、机器学习方法和深度学习方法。全局方法误检和漏检较少,结果较为可信,但是耗时较长,计算流体力学计算过程往往会产生海量的流场,如果用全局方法一一检测,时间太长,无法接受;局部方法简单快速,易于实现,但存在大量的误检和漏检,结果不可靠;机器学习算法高效精确,但通用性不高,泛化性较差,可扩展性太弱;基于深度学习的方法很少,已有的方法结合全局和局部的优势,针对流场内每个点逐一进行判断其是否属于涡区域,但这些方法得到的结果精度比全局方法低,速度比局部方法慢,只是对两种方法性能的折中。已有的深度学习方法具有以下缺点:(1)针对流场内的点进行逐一判断,这会导致数据的重复计算;(2)检测时间与流场内网格点数目成正比;(3)深度学习模型中使用了全连接层,导致网络参数太多,深度神经网络计算一次的时间较长;(4)预处理时,需要将不规则的物理网格转换为规则的计算网格,使用的输入是计算网格下的速度场,丢失了物理网格的特性。
发明内容
本发明解决的技术问题为:克服现有技术不足,提供一种基于神经网络的流场旋涡检测方法,解决流场内涡结构的自动快速检测问题,实现流场内旋涡的高效精确检测。
本发明解决的技术方案为:一种基于分割网络的流场旋涡检测方法,步骤如下:
第一步:对流场内所有物理网格点打标签,标记其是否属于旋涡区域,若属于则标记为1,否则标记为0,得到标签数据,即标记后的物理网格点;
第二步:根据第一步的流场确定出涡量场,对涡量场进行网格转换,将物理网格中的网格点一一映射到计算网格的网格点,抛弃物理坐标信息,得到计算网格下的涡量场。流场的物理网格和计算网格如图3(a)、(b)所示。
第三步,根据第二步得到的计算网格下的的涡量场,计算该涡量场的均值和标准差;将每个网格上的涡量值减去涡量场的均值后除以标准差,再进行归一化操作,得到归一化后的涡量场;
第四步,对第三步得到的归一化后的涡量场和第一步得到的标签数据同时随机采样,如图4(a)所示,对归一化后的涡量场按照设定的区域大小进行采样,采样后对区域内的网格点对应的标签进行判断,如果区域内标签全部为0或者全部为1,则丢弃该采样数据,剩余的数据为带标签局部涡量场数据;
第五步,构建卷积神经网络CNN,使用第四步的随机采样得到的带标签局部涡量场数据训练该CNN,固定CNN网络参数,得到固定参数的CNN;
第六步,利用第五步得到的固定参数的CNN,对待检测涡区的流场进行涡区检测,步骤如下:
优选的,第一步:对流场内所有物理网格点打标签,标记其是否属于旋涡区域,若属于则标记为1,否则标记为0,得到标签数据,即标记后的物理网格点,具体如下:
步骤1.1,对流场的物理网格中的每个网格点(i,j,k)计算其涡量场,假设用x,y,z表示笛卡尔坐标系中的三个坐标轴方向,则i表示网格点在x方向上处于第i列,j表示网格点在y方向上处于第j行,k表示网格点在z方向处于第k层,i∈{1,2,…,Nx},j∈{1,2,…,Ny},k∈{1,2,…,Nz},Nx,Ny,Nz分别表示在x,y,z方向上网格点的数量。在网格点(i,j,k)处的速度值为vi,j,k=(vxvyvz),vx,vy,vz分别表示在该网格点处的速度在三个方向上的值,网格点处的速度值是流场中已知的变量,无需计算。使用公式(1),根据流场内网格点(i,j,k)的速度vi,j,k,计算网格点(i,j,k)处的涡量值ωi,j,k
Figure BDA0002385753390000031
式中,
Figure BDA0002385753390000032
表示速度分量vz在y方向的偏导数,使用公式(1-1)计算,
Figure BDA0002385753390000033
表示速度分量vy在z方向的偏导数,使用公式(1-2)计算,
Figure BDA0002385753390000034
表示速度分量vx在z方向的偏导数,使用公式(1-3)计算,
Figure BDA0002385753390000035
表示速度分量vz在x方向的偏导数,使用公式(1-4)计算,
Figure BDA0002385753390000041
表示速度分量vy在x方向的偏导数,使用公式(1-5)计算,
Figure BDA0002385753390000042
表示速度分量vx在y方向的偏导数,使用公式(1-6)计算。下面分别说明其计算公式。
Figure BDA0002385753390000043
Figure BDA0002385753390000044
Figure BDA0002385753390000045
Figure BDA0002385753390000046
Figure BDA0002385753390000047
Figure BDA0002385753390000048
公式(1-1)~(1-6中)ξ,η,ζ是曲面坐标系下的三个方向。
Figure BDA0002385753390000049
表示vx,vy,vz对ξ方向求偏导数,计算公式为(1-7)~(1-9);
Figure BDA00023857533900000410
Figure BDA00023857533900000411
表示vx,vy,vz对η方向求偏导数,计算公式为(1-10)~(1-12);
Figure BDA00023857533900000412
Figure BDA00023857533900000413
表示vx,vy,vz对ζ方向求偏导数,计算公式为(1-13)~(1-15);ξxyz表示曲面坐标系下ξ方向对笛卡尔坐标系三个方向x,y,z的偏导数,计算公式为公式(1-16)~(1-18)。ηxyz表示曲面坐标系下η方向对笛卡尔坐标系三个方向x,y,z的偏导数,计算公式为公式(1-19)~(1-21)。ζxyz表示曲面坐标系下ζ方向对笛卡尔坐标系三个方向x,y,z的偏导数,计算公式为公式(1-22)~(1-24)。
Figure BDA0002385753390000051
Figure BDA0002385753390000052
Figure BDA0002385753390000053
Figure BDA0002385753390000054
Figure BDA0002385753390000055
Figure BDA0002385753390000056
Figure BDA0002385753390000057
Figure BDA0002385753390000058
Figure BDA0002385753390000059
Figure BDA00023857533900000510
Figure BDA00023857533900000511
Figure BDA00023857533900000512
Figure BDA00023857533900000513
Figure BDA00023857533900000514
Figure BDA00023857533900000515
Figure BDA00023857533900000516
Figure BDA0002385753390000061
Figure BDA0002385753390000062
公式(1-7)~(1-15)中,(i-1,j,k),(i+1,j,k)表示当前网格点(i,j,k)在x方向上的前一个和后一个网格点,(i,j-1,k),(i,j+1,k)表示当前网格点(i,j,k)在y方向上的前一个和后一个网格点,(i,j,k-1),(i,j,k+1)表示当前网格点(i,j,k)在z方向上的前一个和后一个网格点.
Figure BDA0002385753390000063
表示网格点(i-1,j,k)处在x,y,z三个方向的速度值,
Figure BDA0002385753390000064
表示网格点(i+1,j,k)处在x,y,z三个方向的速度值。
Figure BDA0002385753390000065
表示网格点(i,j-1,k)处在x,y,z三个方向的速度值,
Figure BDA0002385753390000066
表示网格点(i,j+1,k)处在x,y,z三个方向的速度值。
Figure BDA0002385753390000067
表示网格点(i,j,k-1)处在x,y,z三个方向的速度值,
Figure BDA0002385753390000068
表示网格点(i,j,k+1)处在x,y,z三个方向的速度值。上述18个速度值均为流场内已知量,无需计算。Δξ表示网格点(i-1,j,k)和网格点(i+1,j,k)在曲面坐标系下ξ方向的距离,Δη表示网格点(i,j-1,k)和网格点(i,j+1,k)在曲面坐标系下η方向的距离,Δζ表示网格点(i,j,k-1)和网格点(i,j,k+1)在曲面坐标系下ζ方向的距离,Δξ,Δη,Δζ都默认设置为1。
公式(1-16)~(1-24)中,J-1可以采用公式(1-25)表示。
Figure BDA0002385753390000069
公式(1-16)~(1-25)中,xξ,yξ,zξ,xη,yη,zη,xζ,yζ,zζ使用公式(1-26)~(1-34)计算。
Figure BDA00023857533900000610
Figure BDA00023857533900000611
Figure BDA0002385753390000071
Figure BDA0002385753390000072
Figure BDA0002385753390000073
Figure BDA0002385753390000074
Figure BDA0002385753390000075
Figure BDA0002385753390000076
Figure BDA0002385753390000077
公式(1-26)~(1-34)中,xi-1,j,k,yi-1,j,k,zi-1,j,k表示网格点(i-1,j,k)在笛卡尔坐标系的物理坐标值。xi+1,j,k,yi+1,j,k,zi+1,j,k表示网格点(i+1,j,k)在笛卡尔坐标系的物理坐标值。xi,j-1,k,yi,j-1,k,zi,j-1,k表示网格点(i,j-1,k)在笛卡尔坐标系的物理坐标值。xi,j+1,k,yi,j+1,k,zi,j+1,k表示网格点(i,j+1,k)在笛卡尔坐标系的物理坐标值。xi,j,k-1,yi,j,k-1,zi,j,k-1表示网格点(i,j,k-1)在笛卡尔坐标系的物理坐标值。xi,j,k+1,yi,j,k+1,zi,j,k+1表示网格点(i,j,k+1)在笛卡尔坐标系的物理坐标值。上述物理坐标值均可以直接得到,无需计算。Δξ表示网格点(i-1,j,k)和网格点(i+1,j,k)在曲面坐标系下ξ方向的距离,Δη表示网格点(i,j-1,k)和网格点(i,j+1,k)在曲面坐标系下η方向的距离,Δζ表示网格点(i,j,k-1)和网格点(i,j,k+1)在曲面坐标系下ζ方向的距离,Δξ,Δη,Δζ都默认设置为1。
至此,计算得到网格点(i,j,k)的涡量值ωi,j,k,逐个网格点计算后得到所有网格点的涡量值,进而形成涡量场。
步骤1.2,根据步骤1.1得到的涡量场计算平均涡量场。具体计算方式如公式(2)所示。
Figure BDA0002385753390000081
公式(2)中,ωμ表示待计算的平均涡量值,ωi,j,k表示步骤1.1中得到的网格点(i,j,k)的涡量值。Nx,Ny,Nz分别表示流场在x,y,z方向上网格点的维度,是已知的量。
使用公式(3)计算瞬时涡量偏差IVD值。
IVDi,j,k=|ωi,j,kμ| 公式(3)
公式中ωμ表示待计算的平均涡量值,ωi,j,k表示步骤1.1中得到的网格点(i,j,k)的涡量值,IVDi,j,k表示网格点(i,j,k)处的瞬时涡量偏差值。逐个网格点计算后得到所有网格点的IVD值。
步骤1.3,找到流场中IVD值的所有局部最大值,并得到这些局部最大值的位置,获取包含IVD局部最大值的所有等值线,得到这些等值线的凸度偏差值。假设第i条等值线所包围的区域面积为Si,该等值线包围的区域中必然存在凸包区域,这是已知的定理,假设该凸包区域面积为Scvi,则可计算出该等值线的凸度偏差值为
Figure BDA0002385753390000082
计算出所有等值线的凸度偏差值后,得到最大的凸度偏差值对应的等值线,则认为该等值线包围的区域属于涡区。将凸度偏差值最大的等值线包围的所有网格点标记为1,其余位置标记为0,至此获得流场内所有网格点处的标签信息即标签数据。
优选的,第二步:对步骤1.1中得到的涡量场进行网格转换,将物理网格中的网格点一一映射到规则的计算网格的网格点,抛弃物理坐标信息,得到计算网格下的涡量场。流场的物理网格和计算网格如图3(a)、(b)所示。
优选的,第三步,根据第二步得到的计算网格下的的涡量场,计算该涡量场的均值和标准差;将每个网格上的涡量值减去涡量场的均值后除以标准差,再进行归一化操作,得到归一化后的涡量场,具体如下:
根据第二步得到的计算网格下的的涡量场,使用公式(6)与公式(7)计算涡量场的均值ωμ和标准差ωδ,将涡量场减去均值后除以标准差,如公式(8)所示,之后使用公式(9)进行归一化操作;
Figure BDA0002385753390000091
Figure BDA0002385753390000092
其中ωi,j,k表示步骤1.1中得到的网格点(i,j,k)的涡量值,Nx,Ny,Nz分别表示流场在x,y,z方向上网格点的维度,是已知的量。
根据公式(6)和公式(7)得到的均值和标准差,使用公式(8)计算其去中心化后网格点(i,j,k)处的正则化涡量值
Figure BDA0002385753390000093
其中i={1,2,…,Nx},j={1,2,…,Ny},k={1,2,…,Nz}。
Figure BDA0002385753390000094
之后使用公式(9)计算归一化后的涡量值
Figure BDA0002385753390000095
遍历所有网格点,计算归一化后的涡量值,形成归一化后的涡量场。
优选的,第四步,对第三步得到的归一化后的涡量场和第一步得到的标签数据同时随机采样,对归一化后的涡量场按照设定的区域大小进行采样,采样后对所有采样块内的网格点对应的标签进行判断,如果采样块内标签全部为0或者全部为1,则丢弃该采样数据,剩余的采样块为带标签局部涡量场数据。随机采样如图4(a)所示。
优选的,第五步,构建卷积神经网络CNN,使用第四步的随机采样得到的带标签局部涡量场数据训练CNN,得到固定参数的CNN,具体如下:
假设CNN共有n层,则CNN的函数表达为:
O=fn(fn-1(fn-2(…f1(ω,W1)…),Wn-1),Wn),其中f1…fn表示第一层到第n层的运算,如果第i层为卷积层,则fi表示卷积运算;W1,…,Wn表示第一层到第n层的待训练的参数矩阵,O表示CNN的输出,输出维度与采样块大小相同,ω表示输入,即采样涡量块。在确定了层数及对应层的类型后,通过带标签的采样块可以不断调整各层的参数矩阵。长时间训练后,可以固定网络参数,得到固定参数的CNN。
优选的,第六步,利用第五步得到的固定参数的CNN,对待检测涡区的流场进行涡区检测,步骤如下:6.1、根据待检测涡区的流场,确定出涡量场,对该涡量场进行网格转换,将物理网格中的网格点一一映射到计算网格的网格点,抛弃物理坐标信息,得到计算网格下的涡量场。
6.2根据步骤6.1得到的计算网格下的的涡量场,计算该涡量场的均值和标准差;将每个网格上的涡量值减去涡量场的均值后除以标准差,再进行归一化操作,得到归一化后的涡量场;
6.3对归一化后的涡量场按照设定的区域大小和位置进行顺序采样,如图4(b)所示,得到采样后的局部涡量场数据;
6.4将步骤6.3采样的局部涡量场数据,送入第五步的固定参数的CNN中,得到CNN输出,CNN输出为数据块;根据神经网络输出中每个网格点的值,判断每个网格点否属于涡,若网格点的值为1,则判定属于涡,打上标签,否则判定为不属于涡,打上标签,形成标签块。
6.5将所有标签块按照位置顺序合并,得到整个流场内所有网格点的标签。
其中,6.1、根据待检测涡区的流场的速度值,确定出待检测流场的涡量场,优选方案如下:
对待检测流场的物理网格中的每个网格点(i,j,k)计算其涡量场,假设用x,y,z表示笛卡尔坐标系中的三个坐标轴方向,则i表示网格点在x方向上处于第i列,j表示网格点在y方向上处于第j行,k表示网格点在z方向处于第k层。在网格点处的速度值为v'i,j,k=(v'xv'yv'z),v'x,v'y,v'z分别表示在该网格点处的速度在三个方向上的值,网格点处的速度值是流场中已知的变量,无需计算。使用公式(10),根据流场内网格点(i,j,k)的速度v'i,j,k,计算网格点(i,j,k)处的涡量值ω'i,j,k
Figure BDA0002385753390000111
式中,
Figure BDA0002385753390000112
表示速度分量v'z在y方向的偏导数,使用公式(10-1)计算,
Figure BDA0002385753390000113
表示速度分量v'y在z方向的偏导数,使用公式(10-2)计算,
Figure BDA0002385753390000114
表示速度分量v'x在z方向的偏导数,使用公式(10-3)计算,
Figure BDA0002385753390000115
表示速度分量v'z在x方向的偏导数,使用公式(10-4)计算,
Figure BDA0002385753390000116
表示速度分量v'y在x方向的偏导数,使用公式(10-5)计算,
Figure BDA0002385753390000117
表示速度分量v'x在y方向的偏导数,使用公式(10-6)计算。下面分别说明其计算公式。
Figure BDA0002385753390000118
Figure BDA0002385753390000119
Figure BDA00023857533900001110
Figure BDA00023857533900001111
Figure BDA00023857533900001112
Figure BDA00023857533900001113
公式(10-1)~(10-6中)ξ,η,ζ是曲面坐标系下的三个方向。
Figure BDA00023857533900001114
Figure BDA0002385753390000121
表示v'x,v'y,v'z对ξ方向求偏导数,计算公式为(10-7)~(10-9);
Figure BDA0002385753390000122
Figure BDA0002385753390000123
表示v'x,v'y,v'z对η方向求偏导数,计算公式为(10-10)~(10-12);
Figure BDA0002385753390000124
表示v'x,v'y,v'z对ζ方向求偏导数,计算公式为(10-13)~(10-15);ξxyz表示曲面坐标系下ξ方向对笛卡尔坐标系三个方向x,y,z的偏导数,计算公式为公式(10-16)~(10-18)。ηxyz表示曲面坐标系下η方向对笛卡尔坐标系三个方向x,y,z的偏导数,计算公式为公式(10-19)~(10-21)。ζxyz表示曲面坐标系下ζ方向对笛卡尔坐标系三个方向x,y,z的偏导数,计算公式为公式(10-22)~(10-24)。
Figure BDA0002385753390000125
Figure BDA0002385753390000126
Figure BDA0002385753390000127
Figure BDA0002385753390000128
Figure BDA0002385753390000129
Figure BDA00023857533900001210
Figure BDA00023857533900001211
Figure BDA00023857533900001212
Figure BDA00023857533900001213
Figure BDA0002385753390000131
Figure BDA0002385753390000132
Figure BDA0002385753390000133
Figure BDA0002385753390000134
Figure BDA0002385753390000135
Figure BDA0002385753390000136
Figure BDA0002385753390000137
Figure BDA0002385753390000138
Figure BDA0002385753390000139
公式(10-7)~(10-15)中,(i-1,j,k),(i+1,j,k)表示当前网格点(i,j,k)在x方向上的前一个和后一个网格点,(i,j-1,k),(i,j+1,k)表示当前网格点(i,j,k)在y方向上的前一个和后一个网格点,(i,j,k-1),(i,j,k+1)表示当前网格点(i,j,k)在z方向上的前一个和后一个网格点.
Figure BDA00023857533900001310
表示网格点(i-1,j,k)处在x,y,z三个方向的速度值,
Figure BDA00023857533900001311
表示网格点(i+1,j,k)处在x,y,z三个方向的速度值。
Figure BDA00023857533900001312
表示网格点(i,j-1,k)处在x,y,z三个方向的速度值,
Figure BDA00023857533900001313
表示网格点(i,j+1,k)处在x,y,z三个方向的速度值。
Figure BDA00023857533900001314
表示网格点(i,j,k-1)处在x,y,z三个方向的速度值,
Figure BDA00023857533900001315
表示网格点(i,j,k+1)处在x,y,z三个方向的速度值。上述18个速度值均为流场内已知量,无需计算。Δξ表示网格点(i-1,j,k)和网格点(i+1,j,k)在曲面坐标系下ξ方向的距离,Δη表示网格点(i,j-1,k)和网格点(i,j+1,k)在曲面坐标系下η方向的距离,Δζ表示网格点(i,j,k-1)和网格点(i,j,k+1)在曲面坐标系下ζ方向的距离,Δξ,Δη,Δζ都默认设置为1。
公式(10-16)~(10-24)中,J-1可以采用公式(10-25)表示。
Figure BDA0002385753390000141
公式(10-16)~(10-25)中,xξ,yξ,zξ,xη,yη,zη,xζ,yζ,zζ使用公式(10-26)~(10-34)计算。
Figure BDA0002385753390000142
Figure BDA0002385753390000143
Figure BDA0002385753390000144
Figure BDA0002385753390000145
Figure BDA0002385753390000146
Figure BDA0002385753390000147
Figure BDA0002385753390000148
Figure BDA0002385753390000149
Figure BDA00023857533900001410
公式(10-26)~(10-34)中,xi-1,j,k,yi-1,j,k,zi-1,j,k表示网格点(i-1,j,k)在笛卡尔坐标系的物理坐标值。xi+1,j,k,yi+1,j,k,zi+1,j,k表示网格点(i+1,j,k)在笛卡尔坐标系的物理坐标值。xi,j-1,k,yi,j-1,k,zi,j-1,k表示网格点(i,j-1,k)在笛卡尔坐标系的物理坐标值。xi,j+1,k,yi,j+1,k,zi,j+1,k表示网格点(i,j+1,k)在笛卡尔坐标系的物理坐标值。xi,j,k-1,yi,j,k-1,zi,j,k-1表示网格点(i,j,k-1)在笛卡尔坐标系的物理坐标值。xi,j,k+1,yi,j,k+1,zi,j,k+1表示网格点(i,j,k+1)在笛卡尔坐标系的物理坐标值。上述物理坐标值均可以直接得到,无需计算。Δξ表示网格点(i-1,j,k)和网格点(i+1,j,k)在曲面坐标系下ξ方向的距离,Δη表示网格点(i,j-1,k)和网格点(i,j+1,k)在曲面坐标系下η方向的距离,Δζ表示网格点(i,j,k-1)和网格点(i,j,k+1)在曲面坐标系下ζ方向的距离,Δξ,Δη,Δζ都默认设置为1。
至此,计算得到待检测流场在网格点(i,j,k)处的涡量值ω'i,j,k,逐个网格点计算后得到所有网格点的涡量值,进而形成待检测流场的涡量场。对该涡量场进行网格转换,将物理网格中的网格点一一映射到计算网格的网格点,抛弃物理坐标信息,得到计算网格下的涡量场。
6.2根据步骤6.1得到的计算网格下的的涡量场,计算该涡量场的均值和标准差;将每个网格上的涡量值减去涡量场的均值后除以标准差,再进行归一化操作,得到归一化后的涡量场,优选方案如下:
根据步骤6.1得到的计算网格下的的涡量场,使用公式(11)与公式(12)计算涡量场的均值ω'μ和标准差ω'δ,将涡量场减去均值后除以标准差,如公式(13)所示,之后使用公式(14)进行归一化操作;
Figure BDA0002385753390000151
Figure BDA0002385753390000152
其中ω'i,j,k表示步骤1.1中得到的网格点(i,j,k)的涡量值,N'x,N'y,N'z分别表示流场在x,y,z方向上网格点的维度,是已知的量。
根据公式(11)和公式(12)得到的均值和标准差,使用公式(13)计算其去中心化后网格点(i,j,k)处的正则化涡量值
Figure BDA0002385753390000161
其中i={1,2,…,Nx’},j={1,2,…,Ny’},k={1,2,…,Nz’}。
Figure BDA0002385753390000162
之后使用公式(14)计算归一化后的涡量值
Figure BDA0002385753390000163
遍历所有网格点,计算归一化后的涡量值,形成归一化后的涡量场。
6.3对归一化后的涡量场按照设定的区域大小m*m和位置进行顺序采样,得到采样后的局部涡量场数据(不带标签),顺序采样方式如图4(b)所示。
6.4将步骤6.3采样的局部涡量场数据,送入第五步的固定参数的神经网络中,得到神经网络输出,大小与输入大小一致,均为m*m;根据神经网络输出的值,判断采样块中每个网格点是否属于涡,若网格点的输出值为1,则判定属于涡,打上标签,否则判定为不属于涡,打上标签,形成局部标签块。
6.5将所有标签块按照位置顺序合并,得到整个流场内所有网格点的标签,优选方案如下:
按照采样局部涡量场数据时的顺序,将局部涡量场对应的局部标签块放置到整个流场中,对于重叠的部分使用覆盖策略判断,即总是使用后面标签块中的标签覆盖前面标签块中的标签。最终得到待检测流场内所有网格点的标签。
本发明与现有技术相比的优点在于:
(1)本发明设计了流场中涡量场的计算方式。相比于速度场,涡量场融合了速度信息和物理网格信息,可以保持物理网格特性;
(2)本发明逐块判断减少了数据重复计算,减少了单个流场的检测时间;
(3)本发明提出了涡量场的归一化方式,通过这种归一化方式,能够保证将涡量场去中心化,得到无偏差归一化后的涡量场,有利于神经网络的训练;
(4)本发明首次提出用采样后的局部涡量场做神经网络输入,这使神经网络的输入与流场本身大小无关,因此对于任意大小的流场,均可以采用本发明的方法进行涡区检测,具有较强的泛化性。
附图说明
图1为本发明整体架构图;
图2为本发明训练和测试流程图,其中(a)为训练过程流程示意图,(b)为测试过程流程示意图;
图3为本发明不规则的物理网格图和转换后的规则计算网格图,其中(a)为不规则的物理网格图,(b)为转换后的规则计算网格图;
图4为本发明随机采样和顺序采样示意图,其中(a)为随机采样示意图;(b)为顺序采样示意图;
图5为本发明的方法流程图。
具体实施方式
下面结合附图和具体实施例对本发明做进一步详细描述。
本发明一种基于神经网络的流场旋涡检测方法,对流场内所有物理网格点打标签,标记其是否属于旋涡区域,得到标签数据;根据第一步的流场确定出涡量场,对涡量场进行网格转换,将物理网格中的网格点一一映射到计算网格的网格点,抛弃物理坐标信息,得到计算网格下的涡量场。根据第二步得到的计算网格下的的涡量场,计算该涡量场的均值和标准差;将每个网格上的涡量值减去涡量场的均值后除以标准差,再进行归一化操作,得到归一化后的涡量场;对第三步得到的归一化后的涡量场和第一步得到的标签数据同时随机采样,对归一化后的涡量场按照设定的区域大小进行采样,采样后对区域内的网格点对应的标签进行判断,得到带标签局部涡量场数据;构建神经网络,使用第四步的随机采样得到的带标签局部涡量场数据训练该神经网络,固定神经网络参数,得到固定参数的神经网络;利用第五步得到的固定参数的神经网络,对待检测涡区的流场进行涡区检测。
本发明所涉及的计算流体力学(CFD)中的旋涡结构非常重要,及时的发现和跟踪能够有效减少一些负面影响,例如飞机飞行时,产生了旋涡会减少飞机的升力。如果能在CFD计算中分析得到在何种情况下产生旋涡,就可以有效避免飞机失速等消极作用的产生。在CFD计算中,为防止丢掉关键信息,总是间隔一定时间步就保存一次结果,因此一次CFD计算可能需要保存几万到几十万个流场。同时,每个流场数据比较大,单个流场数据量最大可达10GB,在如此大量的数据中找哪些时间步存在旋涡是很困难的。已有的分析工具采用传统算法从海量数据中提取旋涡,算法不够智能,时间开销大,结果不够精确。要想快速准确的掌握海量数据规律,必须结合机器学习技术,由机器代替人去挖掘信息,获取知识。本发明结合机器学习技术,采用卷积神经网络CNN,能够高效精确地从大量流场中检测出旋涡,从而大大解放了劳动力,节省了时间。
如图1所示,为本发明采用卷积申请网络进行流场旋涡检测的原理,预处理是在训练过程中步骤1至步骤4的方案,在测试过程中步骤6.1到步骤6.3的的方案,后处理是指在测试过程中步骤6.5的方案,可视化是指将流场的标签数据形成可视的表达方式。
如图5所示,本发明的一种基于分割网络的流场旋涡检测方法,进一步的优方案如下
第一步:对流场内所有物理网格点打标签,标记其是否属于旋涡区域,若属于则标记为1,否则标记为0,得到标签数据,即标记后的物理网格点;
步骤1.1,对流场的物理网格中的每个网格点(i,j,k)计算其涡量场,假设用x,y,z表示笛卡尔坐标系中的三个坐标轴方向,则i表示网格点在x方向上处于第i列,j表示网格点在y方向上处于第j行,k表示网格点在z方向处于第k层,i∈{1,2,…,Nx},j∈{1,2,…,Ny},k∈{1,2,…,Nz},Nx,Ny,Nz分别表示在x,y,z方向上网格点的数量。在网格点(i,j,k)处的速度值为vi,j,k=(vxvyvz),vx,vy,vz分别表示在该网格点处的速度在三个方向上的值,网格点处的速度值是流场中已知的变量,无需计算。使用公式(15),根据流场内网格点(i,j,k)的速度vi,j,k,计算网格点(i,j,k)处的涡量值ωi,j,k
Figure BDA0002385753390000191
式中,
Figure BDA0002385753390000192
表示速度分量vz在y方向的偏导数,使用公式(15-1)计算,
Figure BDA0002385753390000193
表示速度分量vy在z方向的偏导数,使用公式(15-2)计算,
Figure BDA0002385753390000194
表示速度分量vx在z方向的偏导数,使用公式(15-3)计算,
Figure BDA0002385753390000195
表示速度分量vz在x方向的偏导数,使用公式(15-4)计算,
Figure BDA0002385753390000196
表示速度分量vy在x方向的偏导数,使用公式(15-5)计算,
Figure BDA0002385753390000197
表示速度分量vx在y方向的偏导数,使用公式(15-6)计算。下面分别说明其计算公式。
Figure BDA0002385753390000198
Figure BDA0002385753390000199
Figure BDA00023857533900001910
Figure BDA00023857533900001911
Figure BDA00023857533900001912
Figure BDA00023857533900001913
公式(15-1)~(15-6)中ξ,η,ζ是曲面坐标系下的三个方向。
Figure BDA0002385753390000201
Figure BDA0002385753390000202
表示vx,vy,vz对ξ方向求偏导数,计算公式为(15-7)~(15-9);
Figure BDA0002385753390000203
Figure BDA0002385753390000204
表示vx,vy,vz对η方向求偏导数,计算公式为(15-10)~(15-12);
Figure BDA0002385753390000205
Figure BDA0002385753390000206
表示vx,vy,vz对ζ方向求偏导数,计算公式为(15-13)-(15-15);ξxyz表示曲面坐标系下ξ方向对笛卡尔坐标系三个方向x,y,z的偏导数,计算公式为公式(15-16)~(15-18)。ηxyz表示曲面坐标系下η方向对笛卡尔坐标系三个方向x,y,z的偏导数,计算公式为公式(15-19)-(15-21)。ζxyz表示曲面坐标系下ζ方向对笛卡尔坐标系三个方向x,y,z的偏导数,计算公式为公式(15-22)-(15-24)。
Figure BDA0002385753390000207
Figure BDA0002385753390000208
Figure BDA0002385753390000209
Figure BDA00023857533900002010
Figure BDA00023857533900002011
Figure BDA00023857533900002012
Figure BDA00023857533900002013
Figure BDA00023857533900002014
Figure BDA0002385753390000211
Figure BDA0002385753390000212
Figure BDA0002385753390000213
Figure BDA0002385753390000214
Figure BDA0002385753390000215
Figure BDA0002385753390000216
Figure BDA0002385753390000217
Figure BDA0002385753390000218
Figure BDA0002385753390000219
Figure BDA00023857533900002110
公式(15-7)~(15-15)中,(i-1,j,k),(i+1,j,k)表示当前网格点(i,j,k)在x方向上的前一个和后一个网格点,(i,j-1,k),(i,j+1,k)表示当前网格点(i,j,k)在y方向上的前一个和后一个网格点,(i,j,k-1),(i,j,k+1)表示当前网格点(i,j,k)在z方向上的前一个和后一个网格点.
Figure BDA00023857533900002111
表示网格点(i-1,j,k)处在x,y,z三个方向的速度值,
Figure BDA00023857533900002112
表示网格点(i+1,j,k)处在x,y,z三个方向的速度值。
Figure BDA00023857533900002113
表示网格点(i,j-1,k)处在x,y,z三个方向的速度值,
Figure BDA00023857533900002114
表示网格点(i,j+1,k)处在x,y,z三个方向的速度值。
Figure BDA00023857533900002115
表示网格点(i,j,k-1)处在x,y,z三个方向的速度值,
Figure BDA00023857533900002116
表示网格点(i,j,k+1)处在x,y,z三个方向的速度值。上述18个速度值均为流场内已知量,无需计算。Δξ表示网格点(i-1,j,k)和网格点(i+1,j,k)在曲面坐标系下ξ方向的距离,Δη表示网格点(i,j-1,k)和网格点(i,j+1,k)在曲面坐标系下η方向的距离,Δζ表示网格点(i,j,k-1)和网格点(i,j,k+1)在曲面坐标系下ζ方向的距离,Δξ,Δη,Δζ都默认设置为1。
公式(15-16)~(15-24)中,J-1可以采用公式(15-25)表示。
Figure BDA0002385753390000221
公式(15-16)~(15-25)中,xξ,yξ,zξ,xη,yη,zη,xζ,yζ,zζ使用公式(15-26)~(15-34)计算。
Figure BDA0002385753390000222
Figure BDA0002385753390000223
Figure BDA0002385753390000224
Figure BDA0002385753390000225
Figure BDA0002385753390000226
Figure BDA0002385753390000227
Figure BDA0002385753390000228
Figure BDA0002385753390000229
Figure BDA00023857533900002210
公式(15-26)~(15-34)中,xi-1,j,k,yi-1,j,k,zi-1,j,k表示网格点(i-1,j,k)在笛卡尔坐标系的物理坐标值。xi+1,j,k,yi+1,j,k,zi+1,j,k表示网格点(i+1,j,k)在笛卡尔坐标系的物理坐标值。xi,j-1,k,yi,j-1,k,zi,j-1,k表示网格点(i,j-1,k)在笛卡尔坐标系的物理坐标值。xi,j+1,k,yi,j+1,k,zi,j+1,k表示网格点(i,j+1,k)在笛卡尔坐标系的物理坐标值。xi,j,k-1,yi,j,k-1,zi,j,k-1表示网格点(i,j,k-1)在笛卡尔坐标系的物理坐标值。xi,j,k+1,yi,j,k+1,zi,j,k+1表示网格点(i,j,k+1)在笛卡尔坐标系的物理坐标值。上述物理坐标值均可以直接得到,无需计算。Δξ表示网格点(i-1,j,k)和网格点(i+1,j,k)在曲面坐标系下ξ方向的距离,Δη表示网格点(i,j-1,k)和网格点(i,j+1,k)在曲面坐标系下η方向的距离,Δζ表示网格点(i,j,k-1)和网格点(i,j,k+1)在曲面坐标系下ζ方向的距离,Δξ,Δη,Δζ都默认设置为1。
至此,计算得到网格点(i,j,k)的涡量值ωi,j,k,逐个网格点计算后得到所有网格点的涡量值,进而形成涡量场。
步骤1.2,根据步骤1.1得到的涡量场计算平均涡量场。具体计算方式如公式(16)所示。
Figure BDA0002385753390000231
公式(16)中,ωμ表示待计算的平均涡量值,ωi,j,k表示步骤1.1中得到的网格点(i,j,k)的涡量值。Nx,Ny,Nz分别表示流场在x,y,z方向上网格点的维度,是已知的量。
使用公式(17)计算瞬时涡量偏差IVD值。
IVDi,j,k=|ωi,j,kμ| 公式(17)
公式中ωμ表示待计算的平均涡量值,ωi,j,k表示步骤1.1中得到的网格点(i,j,k)的涡量值,IVDi,j,k表示网格点(i,j,k)处的瞬时涡量偏差值。逐个网格点计算后得到所有网格点的IVD值。
步骤1.3,找到流场中IVD值的所有局部最大值,并得到这些局部最大值的位置,获取包含IVD局部最大值的所有等值线,得到这些等值线的凸度偏差值。假设第i条等值线所包围的区域面积为Si,该等值线包围的区域中必然存在凸包区域,这是已知的定理,假设该凸包区域面积为Scvi,则可计算出该等值线的凸度偏差值为
Figure BDA0002385753390000241
计算出所有等值线的凸度偏差值后,得到最大的凸度偏差值对应的等值线,则认为该等值线包围的区域属于涡区。在具体实现过程中,可以直接调用Matlab仿真软件中的convhull函数计算凸度偏差值。将凸度偏差值最大的等值线包围的所有网格点标记为1,其余位置标记为0,至此获得流场内所有网格点处的标签信息即标签数据。
第二步:第二步:对步骤1.1中得到的涡量场进行网格转换,将物理网格中的网格点一一映射到规则的计算网格的网格点,抛弃物理坐标信息,得到计算网格下的涡量场。
第三步,根据第二步得到的计算网格下的的涡量场,计算该涡量场的均值和标准差;将每个网格上的涡量值减去涡量场的均值后除以标准差,再进行归一化操作,得到归一化后的涡量场,具体如下:
根据第二步得到的计算网格下的的涡量场,使用公式(18)与公式(19)计算涡量场的均值ωμ和标准差ωδ,将涡量场减去均值后除以标准差,如公式(20)所示,之后使用公式(21)进行归一化操作;
Figure BDA0002385753390000242
Figure BDA0002385753390000243
其中ωi,j,k表示步骤1.1中得到的网格点(i,j,k)的涡量值,Nx,Ny,Nz分别表示流场在x,y,z方向上网格点的维度,是已知的量。
根据公式(18)和公式(19)得到的均值和标准差,使用公式(20)计算其去中心化后网格点(i,j,k)处的正则化涡量值
Figure BDA0002385753390000244
其中i={1,2,…,Nx},j={1,2,…,Ny},k={1,2,…,Nz}。
Figure BDA0002385753390000251
之后使用公式(21)计算归一化后的涡量值
Figure BDA0002385753390000252
遍历所有网格点,计算归一化后的涡量值,形成归一化后的涡量场。
第四步,第四步,对第三步得到的归一化后的涡量场和第一步得到的标签数据同时随机采样,对归一化后的涡量场按照设定的区域大小进行采样,采样后对所有采样块内的网格点对应的标签进行判断,如果采样块内标签全部为0或者全部为1,则丢弃该采样数据,剩余的采样块为带标签局部涡量场数据。
第五步,构建卷积神经网络CNN,使用第四步的随机采样得到的带标签局部涡量场数据训练CNN,得到固定参数的CNN,如图2(a)所示,优选方案如下:
假设CNN共有n层,则CNN的函数表达为:
O=fn(fn-1(fn-2(…f1(ω,W1)…),Wn-1),Wn),其中f1…fn表示第一层到第n层的运算,如果第i层为卷积层,则fi表示卷积运算;W1,…,Wn表示第一层到第n层的待训练的参数矩阵,O表示CNN的输出,输出维度与采样块大小相同,ω表示输入,即采样涡量块。在确定了层数及对应层的类型后,通过带标签的采样块可以不断调整各层的参数矩阵。长时间训练后,可以固定网络参数,得到固定参数的CNN。
第六步,利用第五步得到的固定参数的CNN,对待检测涡区的流场进行涡区检测,如图2(b)所示,步骤如下:
6.1、根据待检测涡区的流场,确定出涡量场,对该涡量场进行网格转换,将物理网格中的网格点一一映射到计算网格的网格点,抛弃物理坐标信息,得到计算网格下的涡量场。
对待检测流场的物理网格中的每个网格点(i,j,k)计算其涡量场,假设用x,y,z表示笛卡尔坐标系中的三个坐标轴方向,则i表示网格点在x方向上处于第i列,j表示网格点在y方向上处于第j行,k表示网格点在z方向处于第k层。在网格点处的速度值为v'i,j,k=(v'xv'yv'z),v'x,v'y,v'z分别表示在该网格点处的速度在三个方向上的值,网格点处的速度值是流场中已知的变量,无需计算。使用公式(10),根据流场内网格点(i,j,k)的速度v'i,j,k,计算网格点(i,j,k)处的涡量值ω'i,j,k
Figure BDA0002385753390000261
式中,
Figure BDA0002385753390000262
表示速度分量v'z在y方向的偏导数,使用公式(22-1)计算,
Figure BDA0002385753390000263
表示速度分量v'y在z方向的偏导数,使用公式(22-2)计算,
Figure BDA0002385753390000264
表示速度分量v'x在z方向的偏导数,使用公式(22-3)计算,
Figure BDA0002385753390000265
表示速度分量v'z在x方向的偏导数,使用公式(22-4)计算,
Figure BDA0002385753390000266
表示速度分量v'y在x方向的偏导数,使用公式(22-5)计算,
Figure BDA0002385753390000267
表示速度分量v'x在y方向的偏导数,使用公式(22-6)计算。下面分别说明其计算公式。
Figure BDA0002385753390000268
Figure BDA0002385753390000269
Figure BDA00023857533900002610
Figure BDA00023857533900002611
Figure BDA00023857533900002612
Figure BDA00023857533900002613
公式(22-1)~(22-6中)ξ,η,ζ是曲面坐标系下的三个方向。
Figure BDA0002385753390000271
Figure BDA0002385753390000272
表示v'x,v'y,v'z对ξ方向求偏导数,计算公式为(22-7)~(22-9);
Figure BDA0002385753390000273
Figure BDA0002385753390000274
表示v'x,v'y,v'z对η方向求偏导数,计算公式为(22-10)~(22-12);
Figure BDA0002385753390000275
表示v'x,v'y,v'z对ζ方向求偏导数,计算公式为(22-13)~(22-15);ξxyz表示曲面坐标系下ξ方向对笛卡尔坐标系三个方向x,y,z的偏导数,计算公式为公式(22-16)-(22-18)。ηxyz表示曲面坐标系下η方向对笛卡尔坐标系三个方向x,y,z的偏导数,计算公式为公式(22-19)-(22-21)。ζxyz表示曲面坐标系下ζ方向对笛卡尔坐标系三个方向x,y,z的偏导数,计算公式为公式(22-22)-(22-24)。
Figure BDA0002385753390000276
Figure BDA0002385753390000277
Figure BDA0002385753390000278
Figure BDA0002385753390000279
Figure BDA00023857533900002710
Figure BDA00023857533900002711
Figure BDA00023857533900002712
Figure BDA00023857533900002713
Figure BDA0002385753390000281
Figure BDA0002385753390000282
Figure BDA0002385753390000283
Figure BDA0002385753390000284
Figure BDA0002385753390000285
Figure BDA0002385753390000286
Figure BDA0002385753390000287
Figure BDA0002385753390000288
Figure BDA0002385753390000289
Figure BDA00023857533900002810
公式(22-7)~(22-15)中,(i-1,j,k),(i+1,j,k)表示当前网格点(i,j,k)在x方向上的前一个和后一个网格点,(i,j-1,k),(i,j+1,k)表示当前网格点(i,j,k)在y方向上的前一个和后一个网格点,(i,j,k-1),(i,j,k+1)表示当前网格点(i,j,k)在z方向上的前一个和后一个网格点。
Figure BDA00023857533900002811
表示网格点(i-1,j,k)处在x,y,z三个方向的速度值,
Figure BDA00023857533900002812
表示网格点(i+1,j,k)处在x,y,z三个方向的速度值。
Figure BDA00023857533900002813
表示网格点(i,j-1,k)处在x,y,z三个方向的速度值,
Figure BDA00023857533900002816
表示网格点(i,j+1,k)处在x,y,z三个方向的速度值。
Figure BDA00023857533900002814
表示网格点(i,j,k-1)处在x,y,z三个方向的速度值,
Figure BDA00023857533900002815
表示网格点(i,j,k+1)处在x,y,z三个方向的速度值。上述18个速度值均为流场内已知量,无需计算。Δξ表示网格点(i-1,j,k)和网格点(i+1,j,k)在曲面坐标系下ξ方向的距离,Δη表示网格点(i,j-1,k)和网格点(i,j+1,k)在曲面坐标系下η方向的距离,Δζ表示网格点(i,j,k-1)和网格点(i,j,k+1)在曲面坐标系下ζ方向的距离,Δξ,Δη,Δζ都默认设置为1。
公式(22-16)~(22-24)中,J-1可以采用公式(22-25)表示。
Figure BDA0002385753390000291
公式(22-16)~(22-25)中,xξ,yξ,zξ,xη,yη,zη,xζ,yζ,zζ使用公式(22-26)-(22-34)计算。
Figure BDA0002385753390000292
Figure BDA0002385753390000293
Figure BDA0002385753390000294
Figure BDA0002385753390000295
Figure BDA0002385753390000296
Figure BDA0002385753390000297
Figure BDA0002385753390000298
Figure BDA0002385753390000299
Figure BDA00023857533900002910
公式(22-26)~(22-34)中,xi-1,j,k,yi-1,j,k,zi-1,j,k表示网格点(i-1,j,k)在笛卡尔坐标系的物理坐标值。xi+1,j,k,yi+1,j,k,zi+1,j,k表示网格点(i+1,j,k)在笛卡尔坐标系的物理坐标值。xi,j-1,k,yi,j-1,k,zi,j-1,k表示网格点(i,j-1,k)在笛卡尔坐标系的物理坐标值。xi,j+1,k,yi,j+1,k,zi,j+1,k表示网格点(i,j+1,k)在笛卡尔坐标系的物理坐标值。xi,j,k-1,yi,j,k-1,zi,j,k-1表示网格点(i,j,k-1)在笛卡尔坐标系的物理坐标值。xi,j,k+1,yi,j,k+1,zi,j,k+1表示网格点(i,j,k+1)在笛卡尔坐标系的物理坐标值。上述物理坐标值均可以直接得到,无需计算。Δξ表示网格点(i-1,j,k)和网格点(i+1,j,k)在曲面坐标系下ξ方向的距离,Δη表示网格点(i,j-1,k)和网格点(i,j+1,k)在曲面坐标系下η方向的距离,Δζ表示网格点(i,j,k-1)和网格点(i,j,k+1)在曲面坐标系下ζ方向的距离,Δξ,Δη,Δζ都默认设置为1。
至此,计算得到待检测流场在网格点(i,j,k)处的涡量值ω'i,j,k,逐个网格点计算后得到所有网格点的涡量值,进而形成待检测流场的涡量场。对该涡量场进行网格转换,将物理网格中的网格点一一映射到计算网格的网格点,抛弃物理坐标信息,得到计算网格下的涡量场。
6.2.根据步骤6.1得到的计算网格下的的涡量场,计算该涡量场的均值和标准差;将每个网格上的涡量值减去涡量场的均值后除以标准差,再进行归一化操作,得到归一化后的涡量场,优选方案如下:
根据步骤6.1得到的计算网格下的的涡量场,使用公式(23)与公式(24)计算涡量场的均值ω'μ和标准差ω'δ,将涡量场减去均值后除以标准差,如公式(25)所示,之后使用公式(26)进行归一化操作;
Figure BDA0002385753390000301
Figure BDA0002385753390000302
其中ω'i,j,k表示步骤1.1中得到的网格点(i,j,k)的涡量值,N'x,N'y,N'z分别表示流场在x,y,z方向上网格点的维度,是已知的量。
根据公式(23)和公式(24)得到的均值和标准差,使用公式(25)计算其去中心化后网格点(i,j,k)处的正则化涡量值
Figure BDA0002385753390000311
其中i={1,2,…,Nx’},j={1,2,…,Ny’},k={1,2,…,Nz’}。
Figure BDA0002385753390000312
之后使用公式(26)计算归一化后的涡量值
Figure BDA0002385753390000313
遍历所有网格点,计算归一化后的涡量值,形成归一化后的涡量场。
6.3对归一化后的涡量场按照设定的区域大小和位置进行顺序采样,位置选择与区域大小相关,由人为指定,得到采样后的局部涡量场数据;
6.4将步骤6.3采样的局部涡量场数据,送入第五步的固定参数的神经网络中,得到神经网络输出,大小与输入大小一致,均为m*m;根据神经网络输出的值,判断采样块中每个网格点是否属于涡,若网格点的输出值为1,则判定属于涡,打上标签,否则判定为不属于涡,打上标签,形成局部标签块。
6.5将所有标签块按照位置顺序合并,得到整个流场内所有网格点的标签,优选方案如下:
按照采样局部涡量场数据时的顺序,将局部涡量场对应的局部标签块放置到整个流场中,对于重叠的部分使用覆盖策略判断,即总是使用后面标签块中的标签覆盖前面标签块中的标签。最终得到待检测流场内所有网格点的标签。
本发明进一步优选的方案为:对步骤一中的公式(2)进行改进。流场中一般存在物体,例如机翼等,这些物体的存在对流场的影响巨大。因此在设计物理网格的时候,这些物体的边界处网格密度较大,即相邻物理网格点的距离较小,远离物体边界的空间网格密度较小,也就是相邻网格间距较大。这种改进使整个方法集中检测重要性更大的区域中的旋涡,具备更强的实际应用性。
考虑到网格点的这种不平等性质,可以在计算平均涡量场时加入贡献值。如公式(27)所示。
Figure BDA0002385753390000321
αi,j,k表示网格点Pi,j,k的贡献值,也称为权重。
Figure BDA0002385753390000322
d(Pi,j,k,B)表示网格点Pi,j,k和边界点集合B的距离,该距离定义如公式(29)。
Figure BDA0002385753390000323
边界点集合B在流场和物体给定的情况下可以直接得到。公式(29)表示,网格点Pi,j,k与物体边界集合的距离为该网格点与所有物体边界点距离的最小值。
本发明进一步优选的方案为:对步骤三中的公式(6)~公式(8)进行改进。在计算平均涡量场时加入网格点的贡献值。
Figure BDA0002385753390000324
Figure BDA0002385753390000325
Figure BDA0002385753390000326
αi,j,k表示网格点Pi,j,k的贡献值,也称为权重。
Figure BDA0002385753390000327
d(Pi,j,k,B)表示网格点Pi,j,k和边界点集合B的距离,该距离定义如公式(29)。
Figure BDA0002385753390000331
边界点集合B在流场和物体给定的情况下可以直接得到。公式(34)表示,网格点Pi,j,k与物体边界集合的距离为该网格点与所有物体边界点距离的最小值。
本发明的一种基于分割网络的流场旋涡检测方法,进行了仿真,充分体现了本发明的技术效果,以双三角翼流场为例。该双三角翼流场大小为301*201*101,流场的雷诺数为1000000,双三角翼的攻角为22.5°。对该案例进行CFD仿真,可以得到上万个中间流场,我们选择最后一个流场进行旋涡检测。使用精度、召回率和执行时间作为判断依据,精度和召回率如公式(35)和(36)所示:
Figure BDA0002385753390000332
Figure BDA0002385753390000333
p表示精度,其物理含义为我们的方法判断为涡区的网格点中判断正确的比例是多少,比例越大越好。r表示召回率,物理含义为,本发明找到的涡区网格点数量占所有属于涡区的网格点数量的比例,该比例越大越好。
公式(35)和(36)中,TP表示真实标签为1,且使用本发明判断为1的网格点数量;FP表示真实标签为0,本发明判断为1的网格点数量;FN表示真实标签为0,本发明判断为1的网格点数量。这些数量可以统计出来。
本发明在上述双三角翼流场中涡检测精度达到90%,与背景技术中的工作相比提升了6个百分点,召回率r与相关工作比较提升了5个百分点,达到86%。同时,本发明执行时间较短,仅为150秒,与相关工作比较,加速比达到10倍。
本发明设计了流场中涡量场的计算方式。相比于速度场,涡量场融合了速度信息和物理网格信息,可以保持物理网格特性,而且本发明逐块判断减少了数据重复计算,减少了单个流场的检测时间;进一步的本发明提出了涡量场的归一化方式,通过这种归一化方式,能够保证将涡量场去中心化,得到无偏差归一化后的涡量场,有利于神经网络的训练;
本发明首次提出用采样后的局部涡量场做神经网络输入,这使神经网络的输入与流场本身大小无关,因此对于任意大小的流场,均可以采用本发明的方法进行涡区检测,具有较强的泛化性。

Claims (10)

1.一种基于神经网络的流场旋涡检测方法,其特征在于步骤如下:
第一步:对流场内所有物理网格点打标签,标记其是否属于旋涡区域,若属于则标记为1,否则标记为0,得到标签数据,即标记后的物理网格点;
第二步:根据第一步的流场确定出涡量场,对涡量场进行网格转换,将物理网格中的网格点一一映射到计算网格的网格点,抛弃物理坐标信息,得到计算网格下的涡量场;
第三步,根据第二步得到的计算网格下的涡量场,计算该涡量场的均值和标准差;将每个网格上的涡量值减去涡量场的均值后除以标准差,再进行归一化操作,得到归一化后的涡量场;
计算去中心化后网格点(i,j,k)处的正则化涡量值
Figure FDA0003731010760000011
Figure FDA0003731010760000012
其中ωμ,ωδ分别为涡量场的均值和标准差,ωi,j,k表示网格点(i,j,k)的涡量值,i={1,2,…,Nx},j={1,2,…,Ny},k={1,2,…,Nz},Nx,Ny,Nz分别表示流场在x,y,z方向上网格点的维度,是已知的量;
计算归一化后的涡量值:
Figure FDA0003731010760000013
遍历所有网格点,计算归一化后的涡量值,得到归一化后的涡量场;
第四步,对第三步得到的归一化后的涡量场和第一步得到的标签数据同时随机采样,对归一化后的涡量场按照设定的区域大小进行采样,采样后对区域内的网格点对应的标签进行判断,得到带标签局部涡量场数据;
第五步,构建卷积神经网络CNN,使用第四步的随机采样得到的带标签局部涡量场数据训练该CNN,从而固定CNN网络参数,得到固定参数的CNN;
第六步,利用第五步得到的固定参数的CNN,对待检测涡区的流场进行涡区检测。
2.根据权利要求1所述的一种基于神经网络的流场旋涡检测方法,其特征在于:第一步:对流场内所有物理网格点打标签,标记其是否属于旋涡区域,若属于则标记为1,否则标记为0,得到标签数据,即标记后的物理网格点,具体如下:
步骤1.1,对流场的物理网格中的每个网格点(i,j,k)计算其涡量场;
步骤1.2,根据步骤1.1得到的涡量场,计算平均涡量场;根据平均涡量场,计算流场IVD值;
步骤1.3,找到流场中IVD值的所有局部最大值,并得到这些局部最大值的位置,获取包含IVD局部最大值的所有等值线,得到这些等值线的凸度偏差值,将凸度偏差值最大的等值线包围的所有网格点进行标记,至此获得流场内所有网格点处的标签信息即标签数据。
3.根据权利要求2所述的一种基于神经网络的流场旋涡检测方法,其特征在于:步骤1.3中,将凸度偏差值最大的等值线包围的所有网格点进行标记,至此获得流场内所有网格点处的标签信息即标签数据,具体为:
将凸度偏差值最大的等值线包围的所有网格点标记为1,其余位置标记为0,至此获得流场内所有网格点处的标签信息即标签数据。
4.根据权利要求1所述的一种基于神经网络的流场旋涡检测方法,其特征在于:第三步中,根据第二步得到的计算网格下的涡量场,计算该涡量场的均值ωμ,公式如下:
Figure FDA0003731010760000021
5.根据权利要求1所述的一种基于神经网络的流场旋涡检测方法,其特征在于:第三步中,根据第二步得到的计算网格下的涡量场,计算该涡量场的标准差,公式如下:
Figure FDA0003731010760000031
6.根据权利要求1所述的一种基于神经网络的流场旋涡检测方法,其特征在于:第四步中,采样后对区域内的网格点对应的标签进行判断,得到带标签局部涡量场数据,具体为:
采样后对所有采样块内的网格点对应的标签进行判断,如果采样块内标签全部为0或者全部为1,则丢弃该采样块,剩余的采样块为带标签局部涡量场数据。
7.根据权利要求1所述的一种基于神经网络的流场旋涡检测方法,其特征在于:第五步,构建卷积神经网络CNN,使用第四步的随机采样得到的带标签局部涡量场数据训练CNN,得到固定参数的CNN,具体如下:
首先确定CNN层数、CNN的函数表达式以及确定CNN的输入和输出,对于CNN的输出,输出维度与采样块大小相同,在确定了CNN层数及对应层的类型后,通过带标签的采样块不断调整各层的参数矩阵,训练后,固定网络参数,得到固定参数的CNN。
8.根据权利要求1所述的一种基于神经网络的流场旋涡检测方法,其特征在于:第六步,利用第五步得到的固定参数的CNN,对待检测涡区的流场进行涡区检测,步骤如下:
6.1、根据待检测涡区的流场,确定出涡量场,对该涡量场进行网格转换,将物理网格中的网格点一一映射到计算网格的网格点,抛弃物理坐标信息,得到计算网格下的涡量场;
6.2根据步骤6.1得到的计算网格下的涡量场,计算该涡量场的均值和标准差;将每个网格上的涡量值减去涡量场的均值后除以标准差,再进行归一化操作,得到归一化后的涡量场;
6.3对归一化后的涡量场按照设定的区域大小和位置进行顺序采样,得到采样后的局部涡量场数据;
6.4将步骤6.3采样的局部涡量场数据,送入第五步的固定参数的CNN中,得到CNN输出,CNN输出为数据块;根据神经网络输出中每个网格点的值,判断每个网格点否属于涡,若网格点的值为1,则判定属于涡,打上标签,否则判定为不属于涡,打上标签,形成标签块;
6.5将所有标签块按照位置顺序合并,得到整个流场内所有网格点的标签。
9.根据权利要求8所述的一种基于神经网络的流场旋涡检测方法,其特征在于:6.4中,判断每个网格点否属于涡,打上标签,形成标签块,具体如下:
判断每个网格点否属于涡,若网格点的值为1,则判定属于涡,打上标签,否则判定为不属于涡,打上标签,形成标签块。
10.根据权利要求8所述的一种基于神经网络的流场旋涡检测方法,其特征在于:6.5中,将所有标签块按照位置顺序合并,得到整个流场内所有网格点的标签,具体如下:
按照采样局部涡量场数据时的顺序,将局部涡量场对应的局部标签块放置到整个流场中,对于重叠的部分使用设定的覆盖策略判断,即总是使用后面标签块中的标签覆盖前面标签块中的标签;最终得到待检测流场内所有网格点的标签。
CN202010097688.5A 2020-02-17 2020-02-17 一种基于神经网络的流场旋涡检测方法 Active CN111414720B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010097688.5A CN111414720B (zh) 2020-02-17 2020-02-17 一种基于神经网络的流场旋涡检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010097688.5A CN111414720B (zh) 2020-02-17 2020-02-17 一种基于神经网络的流场旋涡检测方法

Publications (2)

Publication Number Publication Date
CN111414720A CN111414720A (zh) 2020-07-14
CN111414720B true CN111414720B (zh) 2023-01-20

Family

ID=71490820

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010097688.5A Active CN111414720B (zh) 2020-02-17 2020-02-17 一种基于神经网络的流场旋涡检测方法

Country Status (1)

Country Link
CN (1) CN111414720B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113743577B (zh) * 2021-06-25 2023-11-21 上海大学 用于中尺度涡识别的精细化网格数据分区构建方法及***
CN113609596B (zh) * 2021-09-29 2021-12-14 中国空气动力研究与发展中心计算空气动力研究所 一种基于神经网络的飞行器气动特性预测方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102921915A (zh) * 2012-10-23 2013-02-13 杭州谱诚泰迪实业有限公司 基于钢水表面旋涡图像识别的下渣检测方法及装置
CN109859311A (zh) * 2019-01-29 2019-06-07 河海大学 一种基于Liutex-Omega涡识别理论的空化流动数值模拟方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102921915A (zh) * 2012-10-23 2013-02-13 杭州谱诚泰迪实业有限公司 基于钢水表面旋涡图像识别的下渣检测方法及装置
CN109859311A (zh) * 2019-01-29 2019-06-07 河海大学 一种基于Liutex-Omega涡识别理论的空化流动数值模拟方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Liang Deng 等.A CNN-based vortex identification method.《Journal of Visualization》.2018,第32卷第65-78页. *
郝滢洁.基于卷积神经网络的海洋中尺度涡旋检测算法研究.《中国优秀硕士学位论文全文数据库 (信息科技辑)》.2018,(第03期),第I140-205页. *

Also Published As

Publication number Publication date
CN111414720A (zh) 2020-07-14

Similar Documents

Publication Publication Date Title
CN111414720B (zh) 一种基于神经网络的流场旋涡检测方法
CN107247259B (zh) 基于神经网络的k分布海杂波形状参数估计方法
CN112966853B (zh) 基于时空残差混合模型的城市路网短时交通流预测方法
CN106970375A (zh) 一种机载激光雷达点云中自动提取建筑物信息的方法
CN110456355B (zh) 一种基于长短时记忆和生成对抗网络的雷达回波外推方法
CN105258906B (zh) 一种风洞自由飞试验模型飞行轨迹预估方法
CN113777931B (zh) 结冰翼型气动模型构造方法、装置、设备及介质
CN110082738B (zh) 基于高斯混合和张量循环神经网络的雷达目标识别方法
CN108332756B (zh) 一种基于拓扑信息的水下航行器协同定位方法
CN108564600A (zh) 运动物体姿态跟踪方法及装置
CN111199105B (zh) 一种扑翼运动参数优化方法
CN114997296A (zh) 一种基于gru-vae模型的无监督航迹异常检测方法及***
CN113627075B (zh) 基于自适应粒子群优化极限学习的弹丸气动系数辨识方法
CN112329627B (zh) 一种高空抛掷物判别方法
CN109255361B (zh) 一种考虑不可行区域的潮汐流能发电场机组布局方法
Du et al. Super resolution generative adversarial networks for multi-fidelity pressure distribution prediction
CN107220589A (zh) 一种基于elm与hmm的序列飞机目标识别方法
CN111651930B (zh) 一种基于极限学习机的流场涡区域检测方法
Ghassemi et al. Adaptive model refinement with batch bayesian sampling for optimization of bio-inspired flow tailoring
CN115620082B (zh) 模型训练方法、头部姿态估计方法、电子设备及存储介质
CN116628886A (zh) 基于有限工程数据的盾构机穿越施工参数实时优化方法
CN116757321A (zh) 太阳直接辐射量预测方法、***、设备及存储介质
KR102480382B1 (ko) 인공지능 기반 풍하중 산정 시스템
CN112347915B (zh) 一种高空抛掷物判别***
CN114492176B (zh) 一种基于生成对抗网络的动力学模型参数辨识方法及***

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant