平板液相横流中气体射流的多尺度旋涡结构分析

黄彪，黄润航，吕亚飞，刘涛涛，郝亮

（1.北京理工大学 机械与车辆学院，北京 100081；2.北京理工大学 重庆创新中心，重庆 401120）

横向射流(Jet in cross flow, JICF)广泛存在于自然界和工业应用中，如涡轮叶片的气膜冷却[1-2]、航空发动机的燃油雾化[3]、火箭发动机推力矢量控制[4-5]以及叶片边界层流动控制[6]等.射流在横流的作用下逐渐发生偏转和弯曲，并与横流之间产生剧烈的能量和动量交换，形成分离、再附、绕流、卷吸等复杂流动现象，进而诱导产生一系列复杂的多尺度旋涡流动结构.湍流结构的演变是射流动量扩散的主要机制，对实际工程应用有着非常重要的影响，一直是横射流研究中学者们关注的核心问题[1-3].

KAMOTANI 等[7]是最早对射流孔下游流场结构做出描述的学者：射流孔下游的横截面上可以观察到一对反旋转涡对(counter rotating vortex, CVP)的存在，同时在下游流场占据主导地位，其影响范围直至射流孔的远下游区域.MARGASON[8]于20 世纪90年代梳理总结了横射流近 50 余年的研究结果，形成了一个理论体系，并将横射流流场的特点概括为反旋转涡对、马蹄涡和尾迹涡等3 个方面.FRIC 等[9]则进一步运用流场显示技术对射流孔下游的流场结构进行了实验观测，指出射流孔附近的涡系还存在自由射流发展过程中的剪切层涡特征，完善形成了四涡系流动结构.至此，以反旋转涡对和剪切层涡构成的主要结构和以马蹄涡和尾迹涡构成的从属结构等四涡系流动结构成为了横射流研究的共识.随着实验技术和数值模拟方法的不断提升，人们开始探究横射流发展过程中各种复杂涡系结构的成因.KOZLOV 等[10]指出剪切层涡的诱因在于射流与横向来流之间速度差引发的K-H 不稳定性，进而诱导交界面的准周期性波动和破碎.同时大量的学者通过研究发现，横射流流场存在强烈的三维和非定常特性，射流初期的流场旋涡结构在发展演化过程中彼此间会发生复杂的相互作用[11-13].

经过70 多年的研究，人们对横射流的湍流结构演化已经有了较为深入的认识，但这些研究多聚焦于气-气单相横射流.美国密西根大学CECCIO 教授研究团队[14-15]基于气层减阻对液相横流中的气体射流开展了卓有成效的研究工作，相关实验通过空气注入的方式在平板模型上构建气相结构，结果表明在不同的气体注入量的条件下会形成气泡减阻、空气层减阻以及介于两者之间的过渡状态减阻等3 种不同的状态.出于研究问题的复杂性，研究者对液相横流中的气体射流多聚焦于气体射流形成的气泡形态发展演化的宏观描述[16-18]，对各种旋涡结构在流场不同发展阶段扮演的作用以及如何影响气泡的演化规律尚未形成系统的认知.近年来，在工程应用的推动下，部分学者对液相横流中气体射流流场中的涡系结构与气泡演化之间的联系进行了初步的探究.王国玉等[19-22]利用不同模型对平板通气进行了数值模拟，捕捉到了射流孔附近的马蹄涡和射流气体内部的反旋转涡对，发现漩涡结构在气层形成不同阶段具有明显的特征.ZHENG 等[23]对回转体通气空化开展了数值计算，结果表明通气空泡气液界面失稳与旋涡结构演化密切相关.马贵辉等[24-25]基于数值结果分析了水下航行体表面气体射流的发展演化机制，捕捉到了气体射流内部的反旋转涡对，并指出气体射流初期阶段的纯气相区主要由反旋转涡对主导.

相较于气-气单相横射流，液相横流中的气体射流流场结构为典型的多相流场，气液两相介质物理属性的显著差异，势必会造成单相横射流中的涡系结构演化过程发生改变，进而造成射流气体和主流水体之间动量交换和传递机制发生改变.为了进一步研究液相横流中气体射流流场结构的发展演化过程，文中基于两步投影法和VOF 界面捕捉方法对水下气体横射流多相流动进行了数值模拟，详细分析了两相横射流流场中旋涡结构的演化特征.

1 数值计算方法

1.1 基本控制方程

文中的研究涉及气液两相，采用均相流模型求解控制方程，采用VOF 方法追踪气液界面.不考虑气液两相介质的可压缩性.

①连续性方程：

②动量方程：

式中：ρ为混合介质密度；ν为混合动力黏度系数；u为速度矢量；p为压力；σ为表面张力系数；k为气液界面的局部曲率；m为自由液面的法向量；δs为气液交界面的Dirac 函数；g为重力加速度.式(2)右侧第3 项[26]为考虑了气液两相界面表面张力加入的表面张力项.

③曲率k：

④法向量m：

式中，α为液相体积分数.

1.2 两步投影方法

文中选用了两步投影方法[27]更新整个流场的压力和速度.首先，将控制方程式(2)代入，当前时间步的速度为un，忽略流场中压力的影响，利用对流项、扩散项、重力项和表面张力项，得到下一个时间步的预估速度u*：

第2 步，即投影步，增加压力的梯度项，来更新下一时间步的速度un+1：

为了求解新时刻的压力场，我们将散度算子运用于式(6)，并使用连续性方程消掉un+1，得到压力的泊松方程:

通过式(7)得到新时刻的压力后，可以利用式(6)获得新时刻的速度场.N-S 方程需要将偏微分方程离散为代数方程，因此需要选用恰当的离散格式.在文中中的对流项使用了二阶迎风格式，扩散项使用了中心差分格式作离散化处理.

1.3 流域边界条件

图1 给出了文中的计算域及边界条件设置示意图，计算域长度为80L、宽度为15L、高度为10L，L=0.01 m.计算域左侧为速度入口，右侧为开放性出口，四周均为固壁边界条件.圆形通气孔为气相速度入口，位于计算域下侧，孔口中心距离速度入口8L，距离前后两侧均为7.5L，直径d=0.453L.其中U∞=1 m/s，Uj=3 m/s，定义射流速度比r=Uj/U∞.空气介质和水介质取25℃时的性质参数.将全局库朗数设置为0.1.

图1 计算域及边界条件设置示意图Fig.1 Schematic diagram of computational domain and boundary conditions

为了达到计算精度和计算资源消耗的最佳平衡，文中采用Richardson 外推法[28]对网格无关性进行分析，采用的网格信息如表1 所示.图2 分别给出了在气体射流充分发展之后，3 种网格数下预测和基于Richardson 外推法获得的孔口中心上方0.01mm 处的垂向速度，可以看出，网格3 的外推相对误差在0.96%以内，满足计算精度的要求.

表1 计算网格划分方案Tab.1 Computational grid division scheme

图2 孔口中心上方0.01 mm 处的垂向速度分布(U∞=1 m/s, Uj=3 m/s)Fig.2 Vertical velocity distribution 0.01 mm above the center of the hole(U∞=1 m/s, Uj=3 m/s)

1.4 数值计算方法验证

为了验证数值计算方法的准确性，图3 分别给出了U∞=1 m/s,Uj=2 m/s 和4 m/s 时，数值模拟获得的特定时刻下气体射流形态与实验结果[29-30]的对比.左侧是实验结果，右侧是数值模拟结果，可以看出气液界面脉动现象较为一致.图4 给出了数值计算获得的气体射流径向扩张角与实验结果的对比，二者变化规律一致，平均误差在3%以内.说明文中建立的数值模型捕捉到的射流气泡形态特征与实验具有较好的相似性.

图3 数值模拟预测的气体射流形态与实验结果的对比Fig.3 Comparison of experimental and numerical gas jet bubbles morphology

图4 数值模拟预测的气体射流径向扩张角与实验结果的对比Fig.4 Comparison of experimental and numerical bubbles expansion angle

2 结果与讨论

2.1 气体射流的发展演化过程

图5 给出了U∞=1 m/s,Uj=3 m/s 时数值预测的气体射流形态随时间的演变过程，图6 给出了射流气泡演化过程中流动区域的分布示意图.从图5 可以看出，射流的发展过程可以分为3 个特征阶段，即射流气泡初步形成阶段、射流气液界面失稳阶段和射流气泡充分发展阶段.随着气体开始通入流场，如图5(a)所示，在横向水流的推动下，射流气体向来流方向弯曲.在沿来流方向运动的同时，气体向横向两侧发展，这时气泡形态较为稳定.当射流气泡沿流向发展一段时间后，如图5(b)所示，由于气液流体密度不同且具有速度差，气液界面在K-H 不稳定性的诱导下开始发生形变并逐渐发生失稳，导致气液两相之间开始发生掺混，进而造成原有的膜状气泡发生局部断裂.断裂后的气泡团在与主流分离后会进一步拉伸变形破碎成为小尺度泡团向下游运动，气泡尾部开始出现湍流结构，可以认为此阶段是层流向湍流进行转变的转捩阶段.此时射流气泡呈现出两个明显的区域，气泡前端为形态稳定的气相区，尾部为界面失稳引起的湍流转捩区.随着射流气泡进一步发展后，如图5(c)所示，扰动不断发展扩大，气液界面变形收缩截断，截断后的气泡团不断破碎分散成许多小气泡，流场结构则从层流转变为湍流.此时气泡流动呈现出3 个典型的流动区域，分别是相对稳定区、转捩区以及湍流区，如图6 所示.

图5 瞬时气泡发展演化过程(U∞=1m/s, Uj=3m/s)Fig.5 Development and evolution of instantaneous bubble

图6 气泡演化过程不同区域示意图Fig.6 Schematic diagram of different regions during bubble evolution

为了进一步说明射流气泡的演化特性，图7 给出了射流气泡不同发展区域气液界面处监测点流向速度随时间的变化.P1点位于相对稳定区，气泡内部始终为纯气相区，流向速度基本不随时间变化，保持为0.7Uj；P2点位于转捩区，在K-H 不稳定性的诱导下气液界面开始出现脉动，导致流向速度趋于接近射流中心速度，并呈低频高幅值脉动，其平均值为0.8Uj；P3点位于湍流区，气液界面失稳导致气液两相间掺混过程加剧，气液界面破碎不断分裂为细碎的小气泡，流向速度呈高频低幅值脉动.

图7 不同监测点处流向速度随时间的演变Fig.7 Evolution of flow velocity with time at different monitoring points

图8 给出了不同截面内垂向速度分布.可以发现，垂向速度沿流向方向基本呈现震荡分布.在通气孔口附近，即在z/d=1 以及x/d=0 处，由于射流气体在刚出孔时具有较高的速度以及动量，并与来流之间发生动量交换，垂向速度脉动幅度较大.在靠近孔口后方位置，射流气体在来流的作用下逐渐与来流方向平行并向下游发展，此时气液界面相对稳定，垂向速度震荡频率较小.在转捩区，气液界面在K-H 不稳定性的诱导下开始出现脉动，导致垂向速度呈高频脉动.随着下游距离x/d的增加，流动进入湍流区，气液两相之间存在较强的混掺作用以及非稳定发展过程，致使速度亏损的值逐渐增大，垂向速度呈高频低幅脉动.

图8 不同截面内垂向速度分布Fig.8 Vertical velocity distribution at different sections

图9 给出了不同位置处流向速度分布.图中蓝色虚线为U/Uj=0.33 速度等值线，红色实线为提取的该位置z/d=0～20 范围内流向速度，云图为y/d=0 截面内以液相体积分数F=0.95 表征的等值面.当速度位于蓝色虚线左边时表示射流气体减速横流，位于虚线右边时表示射流气体加速液体横流.从速度分布中可以看出，在孔口附近处存在一个速度大于主流速度的区域，称之为高速区.高速区速度分布相对稳定，随着截面位置向孔口后方移动，高速区逐渐出现明显的波动，表现为在靠近壁面处呈现低速流动，在接近气液界面处呈现明显的射流气体加速横流的过程.当x/d=13.2 时，由于水气交界面的脉动，在气液界面附近，速度波动更为明显，产生了小尺度旋涡结构，气液界面失稳.随着空泡进一步发展，扰动不断发展扩大，气液界面变形收缩截断，截断后的气泡团不断破碎分散成许多小气泡，表现为基本与来流速度保持一致并呈小范围脉动.

图9 不同位置处流向速度分布Fig.9 Flow velocity distribution at different sections

2.2 多尺度旋涡结构演化过程

液相横流中的气体射流是既包括不同微小尺度的相互作用行为，又具有大尺度气泡演变过程的复杂多尺度流动，这样的流动必然诱导多尺度的旋涡结构.为了进一步阐明流场中旋涡结构发展演化与气体射流形态发展之间的关系，图10～图12 分别给出了气体射流3 个特征阶段的流场中的典型涡系结构，其中旋涡结构以Q涡识别方法进行表征.

图10 射流气泡初步形成阶段Q 分布云图Fig.10 Vortex structure distribution presented by Q at the initial formation stage of jet bubbles

从图10 可以看出，在射流气泡初步形成阶段，由于水流的绕流作用，在孔口前缘处形成了明显的马蹄涡结构，此马蹄涡结构在射流气泡前缘分开，呈螺旋状随气泡向下游发展.在气液两相速度差引起的剪切力作用下，一方面在气泡内部诱导产生反旋转涡对，此涡对起点位于通气孔后缘附近；另一方面在气泡表面诱导产生来源于平板表面的剪切层涡，即近壁面剪切层涡，此涡卷结构尺度沿流向方向不断增大.通过不同流向位置截面的Q分布云图可以看出，反旋转涡对在通气孔附近呈现较为规整的对称形状，随着向下游发展，在平板壁面射流气泡两侧的马蹄涡和沿流向尺度不断增大的近壁面剪切层涡的共同作用下，反旋转涡对逐渐演变为不对称形态，同时伴随着明显的抬升作用，这也说明了气泡在发展过程中会呈现一定的三维效应.同时，反旋转涡对在演变过程中影响范围会不断扩大，卷吸气体沿展向和垂向运动，形成具有一定展向宽度和厚度的贴壁膜状气泡，这说明在气泡的相对稳定区，反旋转涡对主导着流场结构的发展.值得注意的是，反旋转涡对上方，即射流气泡顶端位置处，气液界面的强剪切作用诱导产生了一对旋转方向相反的涡管，此涡管起点位于通气孔中心，并逐渐延伸至气泡尾部，有别于气气横射流流场结构.

随着射流气泡的进一步发展，如图11 所示，由KH 不稳定性引起的气液界面失稳现象开始凸显，主要体现在2 个方面：①近壁面剪切层涡尺度首先不断发展扩大，与气泡内部的反旋转涡对进一步相互作用，造成反旋转涡对不断变形、扭曲，直至消失；其次，大尺度近壁面剪切层涡发生破碎并逐渐演变成大量小尺度涡卷结构，并在水流的夹带作用下向下游运动，在移动过程中会不断发生变形、融合与破碎.近壁面剪切层涡在形成和发展过程中与环境水流体间存在较强的卷吸作用，引起气泡射流束两侧的气液界面发生形变；②射流气泡顶端的双涡管结构遭到破坏，诱发产生新的剪切层涡，即射流剪切涡.这种射流剪切涡在气液交界面处持续作用加强了气液之间的动量和质量交换，诱导气液界面发生准周期性波动和破碎.总而言之，气液界面失稳将引起速度剪切流逐渐向湍流发生转变，气液之间开始发生掺混，气泡形态也随之开始发生改变.

图11 射流气液界面失稳阶段Q 分布云图Fig.11 Vortex structure distribution presented by Q at the instability stage of gas-liquid interface

图12 和图13 分别给出了射流气泡充分发展阶段Q分布云图及不同截面处的涡量分布.结合图12和图13 可以看出，充分发展的射流气泡呈现出3 个典型的流动区域：相对稳定区、转捩区和湍流区.双涡管、反旋转涡对等大尺度旋涡结构主要存在于相对稳定区，大尺度旋涡结构断裂、脱落以及破碎主要发生在转捩区；而湍流区则主要是杂乱无章的小尺度发卡涡结构，这些小尺度旋涡结构之间会发生随机的融合以及再破碎，引起气液之间发生强烈的掺混，进而造成气泡尾部发生一定尺度的断裂脱落.

图12 射流气泡充分发展阶段Q 分布云图Fig.12 Vortex structure distribution presented by Q at the full development stage of jet bubbles

图13 射流气泡充分发展阶段不同截面处涡量分布Fig.13 Vorticity distribution at different sections in the full development stage of jet bubbles

3 结 论

文中基于数值计算对平板液相横流中气体射流进行了研究，分析了气体射流形态的非稳态演化过程，并结合Q涡识别方法，探讨了水下气体横射流流场中旋涡结构的多尺度演化特征，得到的主要结论如下：

①文中基于两步投影法和VOF 界面捕捉方法，结合高分辨率网格策略，建立了平板液相横流中气体射流的高精度数值模拟方法，该方法可以捕捉到流场中气体射流的多尺度流动特征，数值模拟结果与实验结果吻合较好.

②平板液相横流中气体射流的演化过程可以分为射流气泡初步形成、射流气液界面失稳和射流气泡充分发展3 个典型阶段.在Kelvin-Helmholtz 不稳定性的诱导下，气液界面发生形变并逐渐失稳产生破碎，导致射流形态由大尺度连续气泡逐渐转变为小尺度离散气泡，气液界面处的压力和速度则由低频高幅值脉动转变为高频小幅度脉动.

③平板液相横流中气体射流的涡系结构主要包括射流内部和气液界面的反旋转涡对、近壁面和射流界面的剪切层涡以及通气孔前缘的马蹄涡.气泡的相对稳定区主要由射流内部的反旋转涡对所主导；随着气体射流的发展，大尺度近壁面剪切层涡不断与反旋转涡对相互作用，造成反旋转涡对不断变形、扭曲，直至消失，并诱发产生射流界面处的剪切涡，进而导致气液界面发生准周期性波动和破碎.随着大尺度旋涡结构逐渐演变成大量小尺度涡卷结构，气体射流液也由速度剪切流逐渐转变为湍流.