基于拉格朗日方法的通气云状空泡三维非定常脱落特性研究

2020-06-15段忠平刘涛涛张孟杰王国玉

宇航总体技术 2020年3期

段忠平，刘涛涛，张孟杰，黄彪，王国玉

(北京理工大学机械与车辆学院, 北京 100081)

0 引言

对于水下高速航行体而言，为了有效利用通气超空泡的减阻作用，必须在水下航行体表面形成特定尺度、形态可控、且相对稳定的超空泡。然而在大多数情况下，物体是从全湿流过渡到局部空泡，最后才发展为超空泡，因此通气空泡的发展是一个渐近稳定过程：由低频振荡云状空泡、雾状空泡，逐渐转变成半透明的空泡，并在通气量越过某一临界值后，生成包围航行体的稳定超空泡[1]。

通气云状空泡作为通气超空泡发展过程中的中间形态，由于空泡尺度较小，空泡发展过程容易受到流动参数、通气参数、重力效应以及相应的流体属性等多种因素影响而出现空泡团的断裂脱落现象，导致流场发生强烈变化，这种空泡的脉动和断裂等非定常流动特性一直是通气空泡研究的热点问题。Rashidi等[2]运用高速摄影技术观测了通气空泡的流动特性，指出反向射流对空泡界面、空泡的断裂以及空泡内部的流动具有较大影响。Kinzel等[3]采用数值计算对绕回转体通气空泡流动特性进行了研究，阐明了反向射流是引起大尺度空泡脱落的主要原因。王国玉等[4-5]结合实验和数值计算，基于对空泡尾部闭合区域处的压力分析发现，空泡内外压差形成的逆压梯度引起空泡闭合处出现流动分离，形成回射流，回射流与空泡内部气体相互作用产生复杂的旋涡结构，此旋涡结构与主流相互作用造成空泡断裂脱落。王一伟等[6-7]基于反向射流的运动对通气云状空泡脱落机制进行了详细分析，发现通入气体造成反向射流呈蛇形向上游运动，在运动过程中于壁面处产生高压，随着空泡的继续发展和反向射流的持续运动，两个高压区扩大并逐渐连接，促使空泡在高压区的连线形成断裂。这些研究基于对回射流的分析，促进了对通气空泡非定常特性的认识。然而对于轴对称回转体而言，物体往往存在很强的三维效应，因此形成的通气空泡流动也具有三维特性，进而造成反向射流的运动特征也发生相应的改变。Lindau等[8]基于数值计算对绕三维平头回转体的空化流动特性进行了研究，分析了空泡产生的起始点位置，并对空泡尾部的反向射流运动过程进行了描述，说明了其运动的三维特性。Hu等[9-11]基于混合密度分域的PANS湍流模型对平头回转体的云状空化流动进行了数值计算研究，发现云状空穴形态沿回转体周向的分布始终呈现不对称，导致反向射流的推进速度在圆周方向上存在差异，进而造成周向上反向射流到达回转体肩部的时间不一致，引起空泡发生不规则断裂。

由于通气空泡流动的复杂性，限制了人们对其非定常流动机理的有效认识。因此，为了进一步详细分析通气空泡的三维非定常脱落特性，本文基于CFX软件对绕轴对称回转体通气云状空泡发展演化过程进行了数值模拟，并基于FTLE和LCS方法对非定常通气空泡流动进行描述，采用粒子追踪方法对特征流动结构进行分析，揭示通气云状空泡的三维非定常脱落特性。

1 数值计算方法

1.1 基本控制方程

研究基于Navier-Stokes方程的计算框架，采用均质平衡流模型，则Favre平均的N-S方程连续性方程和动量方程可表述为

(1)

(2)

ρm=ρlαl+ρgαg

(3)

式中，下标i和j分别代表坐标方向;ρl，ρg分别为液相和气相密度;αl，αg分别为液相和气相体积分数;ρm，u和p分别为混合介质的密度、速度和压强;μ和μt分别为混合介质的层流和湍流黏性系数。

1.2 湍流模型

由Johansen等[12]提出的滤波器湍流模型(Filter based model，FBM)通过在流域中加入滤波器函数，湍流尺度小于滤波尺度时计算采用k-ε湍流模型进行求解；当滤波尺度小于RANS模型中的湍流特征尺度时采用LES进行计算，能较好地捕捉非定常流动引起的大尺度旋涡结构。

(4)

(5)

湍流黏性由式(6)确定

(6)

式中，k，ε分别为湍动能和湍流耗散率，Pt为湍动能生成项，λ为滤波尺度。相关模型常数分别为：Cε1=1.44，Cε2=1.92，σε=1.3，σk=1.0。

1.3 网格及边界条件设置

为了与实验结果对比，模型的几何尺寸及空间相对位置与实验[13]保持一致。采用的模型为40°钝锥头型回转体，其长度L=126 mm，直径D=20 mm，计算域入口距回转体头部约为12D，出口距回转体尾部为17D，上下壁面距离为4.75D。计算中，基于CFX平台，采用速度入口和压力出口作为边界条件，通气孔设置为质量流量入口，如图1所示。回转体表面采用绝热、无滑移的固壁条件。根据实验条件，入口速度设定为U∞=6.93 m/s，Fr=15.66(Fr=U∞/(gD)1/2)，出口压力设定为p∞=101.325 kPa，空化数σ=4.52(σ=(p∞-pv)/0.5ρv2),通气孔设置为质量流量入口，CQ=0.06(CQ=Qin/(U∞D2))。采用包裹回转体的C型网格拓扑结构来控制回转体周围的网格分布，根据实验结果中通气空泡的区域来确定C型网格拓扑区域的范围。为了可以准确捕捉通气空泡区域内的流场信息以及减少整个计算域的网格数量，包裹回转体C型网格拓扑区域的网格尺度在0.1～0.2 mm，其他区域的网格尺度逐渐增大，网格总数为1.3×106。

图1 计算域边界条件与网格划分Fig.1 The boundary conditions and mesh used in the CFD computations

2 拉格朗日拟序结构

拉格朗日拟序结构(Lagrangian Coherent Structures, LCS)指的是有限时间李雅普诺夫指数(Finite-Time Lyapunov Exponent, FTLE)代表的标量场中出现的“脊”结构，它代表该区域中无限相邻的质点，经过流场的拉扯、卷曲和交错等，将发散至较大的距离。该方法以流体质点在流场中对流特性的时间积分为主要的考虑因素，具有拉格朗日性质的物质线显示旋涡边界，对流场涡的边界和位置有整体的客观描述，可以描述流体质点的真实运动轨迹[14]。

假设流体质点在系统空间内的初始时刻t=t0所处位置为(x0,y0,z0)，经过T时间的运动之后所处位置为(xT,yT,zT)，则有限时间李雅普诺夫指数可以定义为

(7)

(8)

(9)

本次研究中，首先基于CFX商业软件获得通气云状空泡非定常流动流场数据，然后基于Fortran程序对计算所得到的速度场进行处理，获得FTLE结果，然后通过Tecplot软件对计算结果进行显示。

3 结果与分析

3.1 通气云状空泡的发展演化过程

图2分别给出了数值计算得到的通气云状空泡

形态随时间的演化过程，并和实验结果的对比情况。从图中可以看出，数值计算结果和实验结果吻合较好，可以较为准确地描述通气云状空泡发展演化过程中反向射流的推进、空泡断裂、断裂空泡团旋转翻滚以及大尺度空泡团脱落等非定常特性。在该工况下，通气云状空泡的发展周期约为T=17 ms。在t0时刻，反向射流运动到空泡中间位置，此时空泡前端为连续透明的气相区，尾部为水气混合的云雾状，并伴随着上个周期脱落的大尺度云状空泡团。从t0到t0+0.32T时刻，反向射流持续向上游运动，空泡尾部的水气混相区域逐渐扩大，而透明状的气相区域则逐渐缩小，同时空泡尾部出现小尺度空泡团脱落。在空泡前缘附近，反向射流与通入的气流相遇，两股方向相反的流动相互作用，在t0+0.74T时刻，空泡表面出现凹陷和涨起，发生不规则的断裂，涨起空泡在向下游运动的过程中不断翻滚，直至t0+T时刻在空泡尾部再次形成大尺度U型空泡团的脱落现象。值得注意的是，通气云状空泡的整个发展过程中，空泡尾部始终伴随着随机性较强的小尺度云状空泡团的脱落现象，同时由数值结果中的三维效果图可以看出，绕回转体的通气云状空泡呈现很强的不对称性，且周向不同位置出现凹陷裂痕的时刻存在差异。

图2 空泡形态随时间的演化过程
Fig.2 The time-evolution process of ventilated cavity

为了进一步说明通气云状空泡的发展演化过程，图3给出了不同的轴向位置监测点的表面压力系数随时间的变化曲线，其中监测点的位置如图4所示。从图3中可以看出，监测点x/D=0.5，x/D=0.75和x/D=1的压力保持相对稳定，与其他监测点的变化趋势存在显著差异，这是由于这3个监测点位于通气孔附近，受到尾部空泡脱落影响较小，始终位于空泡内部。随着监测点向空泡尾部移动，空泡末端周期性大尺度空泡团的断裂、脱落影响逐渐加剧，对应的压力也出现了低频、高幅度的变化。基于上述的分析可知，对于绕回转体的通气空泡而言，反向射流存在强烈的三维特性，反向射流与主流相互作用后并未造成空泡完全断裂，而是表现为空泡局部断裂，当裂痕发展至空泡尾部后出现脱落，因此图3中只有空泡尾部监测点x/D=3.5，x/D=4压力出现高幅度的变化。

图3 不同轴向位置压力随时间的变化曲线Fig.3 Predicted pressure distribution at different axial positions

图4 回转体表面监测线位置示意图Fig.4 Sketch map of line-monitors at different positions of the axisymmetric body

为了研究通气空泡发展过程的三维特性，图5给出了不同时刻下不同轴向监测线上的表面压力分布，其中4条监测线(A，B,C,D)的位置如图4所示。从图5可以看出，不同监测线上空泡内部低压区的压力分布基本一致。压力的最大峰值P1主要位于流动分离的再附着区空泡末端，且其轴向位置会随着时间的推移发生改变，这说明空泡的脱落周期内空泡末端会发生小尺度空泡团的脱落。同时不同监测线上最大峰值P1的轴向位置在相同时刻也存在着较大差异，这表明空泡末端小尺度空泡团脱落周向位置具有一定的随机性。由于反向射流的作用，在空泡内部低压区会出现压力次高峰P2，其位置会随着时间的推移逐渐向上游移动。对比同一时刻不同监测线上的压力分布可以看出，次高峰P2的轴向位置存在显著差异，这说明反向射流在沿回转体壁面向上游运动的过程中周向方向上存在着差异。

(a) t0

(b) t0+0.25T

(d) t0+0.75T图5 不同时刻下监测线上回转体表面压力系数分布Fig.5 Instantaneous pressure distribution of different monitoring lines

图6统计了不同时刻下P2的位置，其中纵坐标表示P2所处位置与回转体头部之间的无量纲距离。从图中可以看出，曲线的斜率表示各监测线上反向射流的运动速度。在t0时刻，C监测线上P2的位置较A，B和D的靠下游。随着时间的推移，由曲线的斜率可以看出，不同监测线上的反向射流以不同的速率向上游运动，并且同一监测线上不同时刻反向射流的推进速度也不尽相同，其中D监测线上反向射流的平均推进速度最快。结合图2可以看出，当反向射流与通气气流相遇后被推向外侧，反向射流开始远离回转体壁面朝空泡壁面运动，从而引起空泡发生涨起。由于沿回转体表面圆周方向上反向射流的推进速度存在差异，空泡外侧高压造成空泡中部边界出现凹陷，形成局部断裂；随着时间的推移，断裂位置逐渐向两侧扩散，从而U型空泡团脱落。

图6 次高压P2位置随时间的变化Fig.6 Instantaneous positions of P2

3.2 三维非定常空泡断裂特性分析

为了进一步详细地分析通气云状阶段空泡的非定常断裂特性，本节采用拉格朗日方法进行了研究。FTLE可表示流场中各质点的动态趋势，而非限于当前的流场，高FTLE值区域会与低FTLE值区域卷曲交错，进而形成“脊”状结构，即拉格朗日拟序结构，可精确捕捉流场中的旋涡边界。由于FTLE主要用于描述两个相邻粒子在有限时间T范围内的平均变动率，其值与选定的参考时间尺度或周期性运动的周期有着直接联系，在进行相应计算时，采用时间向后积分的方法，积分时间为一个周期T=17 ms。图7给出了通气云状空泡在一个脱落周期内典型时刻下纵截面内的FTLE云图以及示踪粒子的分布情况，图8为相应的含气量云图。初始时刻，如图7(a)所示，在回转体的上下表面空泡内部的5个特征位置：通气孔处、通气孔前、通气孔后、空泡中部、空泡尾部，呈对称状分别布置了5组示踪粒子，图中以不同颜色加以区分并依次编号为A，B，C，D和E。在一段时间内，每个示踪粒子的位置会发生相应的改变，根据粒子的分布情况可以分析得到流场结构的变化，进而可以分析空泡的脱落特性。

观察图7中不同时刻下的FTLE云图可以发现，空泡覆盖区域的LCS结构整体呈椭球状分布，与图8中的空泡外形基本一致；椭球状内部为不规则的复杂LCS结构，与空泡内部形成的若干小尺度旋涡结构一一对应，这说明LCS可以描绘出通气云状空泡内复杂的流动结构。同时由粒子的始末位置可以看出，示踪粒子基本沿着LCS结构运动，回转体上下表面初始位置呈对称分布的粒子经过相同时间的运动后呈现出非对称性，说明绕回转体的通气云状空泡内部流动结构是非对称的。

(a) t0

(b) t0+20.6%T

(d) t0+55.9%T

(e) t0+79.4%T

(f) t0+94.1%T图7 不同时刻纵截面内的FTLE云图及示踪粒子的始末位置Fig.7 The contours of FTLE and the distributions of Lagrangian tracers on the longitudinal section

A，B，C这3组粒子主要沿水气界面向下游运动，并且随着时间的发展发生一定的混合，同时由原来的团状分布变化为带状分布，如图7(b)和图8(a)所示，这说明相同轴向位置不同径向位置处粒子的运动速度存在显著差异。另外，回转体上下表面粒子的运动方式存在显著差异，回转体上表面的粒子主要聚集在通气孔附近且位置基本不随时间改变，下表面的粒子以断面形式逐渐向下游移动，粒子经过后的区域表现为透明的气相区。

D和E两组粒子的运动轨迹较为复杂。在t0+20.6%T时刻，D组粒子在反向射流与主流的剪切拉伸作用下，靠近回转体表面的粒子反向朝空泡前缘移动，远离回转体表面的粒子随主流向下游移动，回转体上表面中反向移动的粒子呈现出整体的反向移动，而下表面反向移动的粒子呈现出先聚集后移动的趋势。到t0+32.4%T时刻时，如图7(c)所示，回转体下表面反向移动的部分D组粒子则进一步聚集并开始偏离回转体表面，说明此处的反向射流可能遭遇前方气流被推向外侧。随着时间的推移，下表面的D组粒子在偏离回转体表面向空泡界面运动的过程中受到空泡界面处C组粒子的影响逐渐转向空泡内部，进而向下游移动，说明反向射流被推向外侧遭遇空泡壁面后再次向内转向，反向射流与空泡界面的多次接触极易引起空泡的断裂脱落。

在空泡的整个脱落周期内，回转体上表面部分E组粒子沿其壁面反向朝回转体肩部移动，部分粒子贴壁向下游发展。回转体下表面反向移动的E组粒子与上表面运动方式基本一致，也呈现出贴壁移动，但由于D组粒子的影响，E组粒子在t0+55.9%T时刻开始出现偏离回转体表面，表明反向射流在此处也被推向外侧。而回转体下表面向下游移动的粒子与上表面的运动方式存在很大差异。在t0+20.6%T时刻，如图7(a)所示，回转体下表面已经出现部分E组粒子偏离回转体表面，说明空泡尾部的反向射流有将空泡抬起的趋势，从而空泡发生卷起并形成小尺度脱落。随着时间的推移，E组粒子与A，B，C，D这4组粒子发生混合并逐渐向下游移动，直至t0+94.1%T时刻，A，B和部分C，D，E组粒子组成的混合粒子团完全脱离空泡，空泡完成大尺度的脱落过程。

(a) t0+20.6%T

(b) t0+55.9%T

(d) t0+94.1%T图8 不同时刻纵截面内的气体含量云图和示踪粒子的始末位置细节Fig.8 The contours of air fraction and the distributions of Lagrangian tracers

为了进一步研究其三维特性，图9给出了绕回转体不同横截面上的FTLE云图和部分示踪粒子的运动轨迹。其中，x表示横截面到回转体头部的距离，分别用回转体直径进行无量纲化处理。在初始时刻t0时，分别在空泡内部沿y=0和z=0的坐标轴线两侧布置A，B，C，D这4组粒子，由于LCS的分布通常与旋涡边界相对应，可以用来界定具有不同流动结构的区域，因此每组粒子依据所处LCS结构位置的不同又细分为3小组，即A1，B1，C1，D1组粒子位于圆环LCS结构内部，A2，B2，C2，D2组粒子位于圆环LCS结构边界上，A3，B3，C3，D3组粒子位于圆环LCS结构外部，如图9所示。

图9 x/D=1截面处初始时刻FTLE云图、示踪粒子位置以及含气量云图Fig.9 The contours of FTLE and air fraction, the tracking particle initial positions at x/D=1 section

对比不同截面的LCS结构分布可以看出，与回转体头部的距离越近，LCS的分布越简单清晰。如图10(a)所示，在靠近回转体头部x/D=1的截面上，LCS整体呈现较为规则的圆环状，从粒子的运动轨迹可知，此截面上的流动主要分为3个区域：1) 圆环LCS结构内部靠近回转体壁面的区域，如A1，B1，C1，D1组粒子所示，此区域具有明显的周向流动，方向或呈顺时针或呈逆时针或呈螺旋状；2) 圆环LCS结构外部远离回转体壁面的区域，粒子主要位于空泡壁面处，如A3，B3，C3，D3组粒子所示，粒子并没有明显的周向运动，此区域的流动方向基本与主流方向(x方向)一致；3) 圆环LCS结构边界区域，如A2，B2，C2，D2组粒子所示，此区域的流动兼具上述两种运动，粒子首先呈现一定的周向运动，在周向运动过程中逐渐靠近空泡壁面，到达空泡壁面后沿主流向下游运动。

随着截面位置向下游移动，如图10(b)所示，在空泡中部x/D=2截面上，LCS的分布区域逐渐扩大同时变得不规则，从粒子的运动轨迹可以看出，此截面上的流动与x/D=1的截面上流动基本一致，不同之处在于x/D=2截面上第一个流动区域中粒子的运动变得更加复杂且无规律性，同时周向运动的粒子越来越多，说明周向流动的区域逐渐扩大。当流动发展至空泡尾部，如图10(c)所示，LCS的分布更加杂乱不规则，已经不再呈现为圆环状，几乎所有粒子都表现出周向的运动，说明空泡尾部的三维运动特性最为明显且复杂，也进一步说明在一个空泡脱落准周期内，空泡尾部容易出现脉动诱发多次的小尺度空泡发生脱落。