平板近壁面双孔通气气液两相流场的数值模拟研究

2021-09-08吕亚飞张孟杰刘涛涛

数字海洋与水下攻防 2021年4期

吕亚飞，张孟杰，刘涛涛，*，陈杰，赵欣

（1.北京理工大学机械与车辆学院，北京 100081；2.中国北方车辆研究所，北京 100072）

0 引言

水下垂直发射航行体具有隐蔽性、机动性和突击性等优点，已成为各军事大国战略发展的重要方向。航行体出水载荷、弹道及姿态是水下发射技术的关键参数，由于穿越自由液面过程中存在复杂因素的干扰，如何确保航行体出水参数满足要求，已经成为亟待解决的水下发射技术难点之一[1]。近年来，国内外均发展了通气空泡技术来改善航行体在出水过程中的力学环境。作为一种特殊的通气空泡技术，多孔排气形成的覆盖航行体局部表面的等压气膜，可以有效改善航行体出水载荷、弹道及姿态[2]，该技术相关的科学问题近年来被提出并受到广泛关注。

多孔排气流场调控技术是利用发射前预置的高压气腔作为气源，在环境降压及贴体水流剪切下依靠向航行体边界层注入不可凝气体而形成，伴随着航行体的出水运动排气气泡持续空间扩张，逐渐由条带状气束发展成覆盖航行体较大部分外表面的周向融合的层状气膜[3]。出于研究问题的复杂性及研究对象的敏感性，国内外有关该技术的文献较少，主要集中在简单机理试验的初步探索，聚焦于气泡形态宏观发展演化过程的描述以及气膜形成后的整体流体动力特性的分析[4-7]，从一定程度上提高了人们对多孔排气流场结构演化过程的认知。

然而，从流动特征来看，多孔排气过程本质上为近壁面边界层同源多束射流流动问题，表现为典型的横射流流动过程。对于单相介质横射流流动，相关的研究表明[8-11]，流场发展过程中会呈现出反向涡对、马蹄涡、尾迹涡和环状剪切层涡等复杂的多尺度旋涡结构，这些旋涡结构在湍流形成、卷吸和流体掺混过程中起到了主导作用。相比于传统的单相横射流流动，水下发射多孔排气形成的流场会更加复杂，其主要原因在于气液两相介质物理属性的巨大差异使得旋涡结构在向下游发展的过程中呈现出更强的各向异性和非定常特性。目前针对水下多孔排气多相湍流流动的 CFD数值模拟，主要采用标准或修正 RANS模型进行求解[12-17]，加深了人们对大尺度空泡“回射流”推进、非定常空泡的流动分离、旋涡脱泻等宏观尺度流动的理解，但受限于商业软件在数学模型和求解器的灵活性，仍然无法做到精确捕捉气体射流与水体主流相互干扰形成的涡系结构呈现出的旋涡对并、合拢等复杂非定常现象，导致预测的流场结构与实验结果存在较大误差。

为了进一步探究水下航行体多孔排气流场演化机制，本文以平板绕流为研究对象，基于发展的近壁面边界层通气气液两相流动高精度算法，分析了近壁面双孔通气气泡的非定常发展演化过程，探讨了流场的旋涡结构分布与气泡发展之间的关系。

1 数值计算方法

1.1 基本控制方程

本文的研究对象为气液两相流动，考虑到均相流模型具有较好的收敛性和计算经济性，研究中将气液两相混合物当作一种密度可变的均匀介质，从而共同求解一套控制方程。同时，为了简化计算，将涉及到的气相和液相都视为不可压缩流体，即把流体的密度视为常量。对于不可压缩流体，笛卡尔坐标系下的连续性方程和动量方程可表述为

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

1.2 VOF方法

本文采用VOF方法对气液界面进行追踪捕捉。整个计算在笛卡尔坐标系下进行，并通过液相体积分数α进行界面追踪和重构，即满足以下定义：

液相体积分数α满足连续性方程

混相流体密度ρ和ν动力粘性系数为

但是，任何事物都存在矛盾体，即“福兮祸所伏”，要想得之先于与之，好东西的获得通常会要求一定的付出。经济的急速发展，伴随而至的是环境被破坏、河水被污染、生态失衡、气候变暖等等，人类终于为了现代化生活而付出了惨痛的代价。中国人民的幸福生活一半是痛苦的，生态破坏、环境恶化给人们带来的痛心的教训。

式中，ρwater和ρair分别为液相和气相密度，νwater和νair分别为液相和气相动力粘性系数。

1.3 求解方法及流程

在建立了控制方程的基础上，应用 VOF方法对两相界面进行追踪捕捉，并获得相界面法向量与曲率，建立的数值计算方法的求解计算流程如图 1所示。

图1 数值模拟计算顺序流程图Fig.1 Flow chart of numerical simulation

第3步根据不可压缩连续性方程可以得到公式（10）左端的速度微分项等于0，即

可以看出，本文采用的两步投影法将原来不可压缩的N-S方程进行了分步求解，通过计算过程流体速度 u*并且利用不可压缩连续性方程构造了一般形式的泊松方程，解除了对传统RANS湍流模型的依赖，实现了对湍流转捩过程的精确捕捉。

1.4 计算域及边界条件的设置

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

2 结果与讨论

2.1 数值结果验证

图3 给出了不同通气速度下数值模拟计算得到的气泡形态与实验结果的对比。

图3 不同通气速度下实验观测的气泡形态与数值预测结果的对比Fig.3 Comparison of experimental and numerical cavity shape under different ventilation speeds

从图中可以看出，本文建立的数值模拟方法计算得到的气泡发展过程与实验结果具有较好的一致性，即气体从通气孔进入流场后紧贴平板壁面形成连续条状气泡，气泡前端为连续透明气相区，气液界面较为光滑；随着流动向下游发展，气泡逐渐向两侧膨胀，并在 K-H不稳定性的影响下气液界面出现明显的脉动，气、液两相发生强烈掺混，导致气泡呈现出不透明的混相区。同时，随着通气速度的逐渐增大，气泡界面处气、液两相速度差增大，气液界面失稳位置逐渐向通气孔靠近，连续透明区逐渐缩小。由于实验中无法完全确保气体持续均匀通入，实验观测到的连续透明区的气泡界面与数值计算结果存在一定的误差。为了进一步定量说明数值计算方法的准确性，图4给出了实验测量得到的不同通气速度下气泡扩张角与数值计算的对比，其中实验结果为多次测量的平均值，可以看出，数值计算得到的气泡扩张趋势与实验结果也吻合较好。

图4 不同通气速度下实验测量的气泡扩张角与数值结果的对比Fig.4 Comparison of experimental and numerical cavity expansion angle sunder different ventilation speeds

2.2 通气气泡的发展演化过程

图5给出了来流速度 U∞=5 m/s、通气速度Ug=10 m/s下平板近壁面双孔通气气泡形态随时间的演变过程，其中气相体积分数取0.9。从图中可以看出，双孔通气气泡的发展呈现出明显的3个阶段：第1个阶段为相对稳定阶段，在此发展阶段，气体在水流的冲击下迅速沿流向弯折，由与水流正交变为与水流流动方向平行并逐渐向下游运动，形成贴附于平板表面且具有一定厚度和展向宽度的膜状气泡，此时气液界面清晰且光滑，气泡内部主要由气相介质构成，气泡形态较为稳定，两股气泡之间呈现出独立发展趋势，没有出现明显的混掺融合现象。第2阶段为气液界面失稳初步融合阶段，此阶段，在K-H不稳定性的影响下，气液界面开始出现失稳，且越靠近气泡尾部气液界面变形越大，同时气泡不断卷吸环境水体进入其内部，造成两股气泡间开始接触并发生混掺融合。由于气液两相间的掺混过程，气泡末端开始出现湍流结构，在湍流的作用下气泡末端呈现出大量离散的多种尺度气泡杂乱分布。第3阶段为充分接触融合阶段，在此阶段，气液界面失稳加剧，气泡末端湍流区不断扩大，湍流流场紊乱的结构则进一步加剧气液之间的掺混，引起两股气泡束充分接触并融合为单一气膜向下游发展，此时气泡界面模糊，呈现出小尺度的离散空泡。另外，值得注意的是，随着气泡间的掺混融合过程的不断发展，气液界面失稳位置和气泡间初步接触融合位置都发生了改变，气液界面失稳位置由初期的x/l=10提前至x/l=17，而气泡间初步接触融合位置则由初期的x/l=15提前至x/l=10，这说明气泡间的掺混融合过程对上游气泡的发展产生了明显的影响。

图5 平板近壁面双孔通气气泡形态随时间的演变过程Fig.5 Evolution of cavity shape for double-hole ventilated cavity near wall of the flat plate

图6给出了平板壁面（y/l=0）处的压力随时间演变过程，其中黑色轮廓线为气体体积分数为0.5时所表征的通气气泡形态。从图中可以看到，当气体从孔口通入流场后，来流发生绕流作用，在孔口附近的两股通气气泡束前端及其周围形成高压区。随着通气气泡的发展，两股通气气泡束的相对稳定区内的压力始终保持较低均匀分布，且比较稳定；而x/l=6位置后的水气掺混区域的压力分布随着通气气泡束界面的波动呈现波动性分布。随着气泡的进一步发展，两股通气气泡束发生掺混融合作用，两股气泡束间出现相对较强的高压、低压间隔分布，直至充分融合后，气泡内部的压力表现为较低压、稳定分布。

图6 平板表面压力随时间的演变过程Fig.6 Evolution of pressure of flat plate surface with time

图7进一步给出了y/l=0平板壁面处的周向速度随时间演变过程，其中黑色轮廓线为气体体积分数0.9所表征的通气气泡形态。速度的正负代表着其周向运动方向，从图中可以看到，气体通入流场后，再水流来流的扰流作用下，孔口位置的通气气泡内部产生均朝向孔中心位置的周向速度；同时，由于通气气泡束对周围环境水体的卷吸作用，在孔口附近x/l=0～5.5范围内所生成的稳定通气气泡外部产生均朝向通气空泡外侧的周向速度。随着气泡的发展，两股通气气泡束之间不断发生着掺混融合作用，两股通气气泡束内侧的方向相反的周向速度不断增强；而随着通气气泡束向周围环境水体的扩散作用，两者之间不断发生着动量交换，从而导致两股通气气泡束外侧的周向速度逐渐减弱直至趋于稳定。此外，当两股通气气泡束充分融合形成单一空泡后，通气气泡尾部的断裂脱落位置存在较强的周向速度分布。

图7 平板表面周向速度随时间的演变过程Fig.7 Evolution of circumferential velocity of flat plate surface with time

2.3 旋涡结构特性

旋涡结构的生成和演化对平板近壁面通气气泡流场的发展演变过程中起着关键性的作用。为了进一步分析流场的结构演化过程，图8给出了t5=42.5 ms时刻下流向速度u染色的Q=0.1等值面所表征的单孔通气气泡流场中的旋涡结构分布。从图中可以看出在近壁面通气气液流场中，分布着多种不同尺度的旋涡结构，主要包括马蹄涡、反向旋转涡对（CVP）、剪切层涡、发夹涡、Λ涡等典型旋涡结构。随着气泡的形成，孔口位置处的水流绕流作用相对较强，从而在通气孔前缘处生成明显的马蹄涡结构。马蹄涡结构分布在相对稳定区域，其与通气气泡在前缘分离开并表现为独立发展状态，分布在通气气泡两侧的涡腿比较长，达到x/l=4.5的位置处，并呈现为对称分布。同时在水流的冲击作用下，通入气体所形成的气泡迅速弯曲，其断面形状呈肾形，而在肾形截面内出现空间尺度、强度较大的反向旋转涡对，即CVP涡。剪切层涡卷则发生在空泡界面和平板近壁面区域，主要由流体间的速度差异所导致，其中气泡界面处剪切层涡主要分布在相对稳定区和转捩区相接区域，且随着气液界面的波动呈规则排列，而壁面处剪切层涡随着气泡的发展沿流向不断增大，分布相对比较规则。另外，值得注意的是，随着气泡的失稳破碎，其在流场下游开始断裂脱落，流动比较紊乱，从而在湍流区域生成大量的发夹涡，如Λ涡、Ω涡等。

图8 Q等值面表征的近壁面单孔通气气泡流场中的旋涡结构分布Fig.8 Vortex structure distribution presented by Q for single-hole ventilated cavity near wall of the flat plate

图9则进一步给出了t5=42.5 ms时刻下流向速度u染色的Q=0.1等值面所表征的双孔通气气泡流场中的旋涡结构的整体分布。从图中可以看出，随着气泡的发展，气泡两侧生成的壁面剪切层涡不断卷吸着周围环境水体向展向扩张，并在两孔中间一定位置发生接触生成新的剪切层涡，新的剪切层涡卷的高度要明显低于周围两束通气气泡界面处的剪切层涡卷。这说明气泡束之间的初步掺混融合过程是由近壁面处的剪切层涡卷主导的，而非气泡界面处的剪切层涡卷。同时，随着两股气泡间掺混融合的进一步发展，融合后产生的剪切层涡卷发生一定程度的抬升，并逐渐发展成为气泡界面处的剪切层涡卷。

图9 Q等值面表征的近壁面双孔通气气泡流场中的旋涡结构分布Fig.9 Vortex structure distribution presented by Q for double-hole ventilated cavity near wall of the flat plate

为了进一步探究反向旋转涡对在气泡束的掺混融合过程中的作用，图10给出了x/l=1.5、9、15.8、27.4等4个不同横截面内的瞬时流向涡量ωx和横截面内的流线分布。从图中可以看到，在孔中心x/l=1.5横截面内，2个通气孔周围的旋涡呈现完全规则的对称分布，反向旋转涡对也分别自平板壁面孔中心处产生并向内侧卷吸。在x/l=9横截面内，2股气泡束内部的旋涡分布保持相似性，根据流线分布可以观察到2束气泡束对应的反向旋转涡对在内侧开始发生相互作用，而空泡轮廓显示此时的两股气泡束的界面还没有发生接触。这说明反向旋转涡对也是诱导气泡间发生掺混融合的原因之一。在x/l=15.8横截面内，气泡内侧的壁面剪切层发生相互作用并出现接触交混现象，分布在内侧中间位置处的 CVP旋涡也相应发生接触并与形成的剪切层涡交汇融合在一起。随着2股气泡束的进一步融合，在流场下游的x/l=27.4横截面内，2股气泡束经过充分融合后表现为单一空泡形态，此时流场中的旋涡分布在空泡内部，而不再分布在空泡界面处，进一步说明所形成的单一空泡是相对稳定的。此外，单一空泡内部的反向旋转涡对变得杂乱无序，旋涡分布表现为紊乱、无序、小尺度涡聚集。