近自由面超空泡航行体多相流特性分析

2023-09-02谢泽阳赵小宇

空气动力学学报 2023年7期

向敏，刘波，谢泽阳，赵小宇

(1.国防科技大学空天科学学院，长沙 410073；2.西安现代控制技术研究所，西安 710000)

0 引言

超空泡技术为通过生成包裹弹体表面的空泡实现高效减阻的技术，已在鱼雷、射弹、水下导弹等武器装备中得到成功应用。超空泡技术应用于水面航行器，可充分利用外界空气作为超空泡气源，在大幅提高水面艇航速的同时保证航行时间和航程，因此具有突出的优势。然而近自由面超空泡由于受到自由界面扰动和空气卷吸掺混效应的影响，其演化动力学机制与传统对称超空泡存在显著差异：一方面，空化将产生兴波效应；另一方面，自由界面将干扰空泡的产生、发展及闭合泄气，因此近自由面空化现象的研究与航行体水动力特性预示充满挑战。目前针对空化现象和自由界面效应相互耦合的研究仍较为有限。

实验研究方面，国内外针对超空泡动力学机制的研究目前主要集中于无自由界面干扰情况下的空泡演化机理研究[1-3]，高速稳定的自由界面空泡流实验研究难度较大。施红辉等[4]采用高速相机，通过对试验数据的分析，研究了水深和自由面对自然超空泡轮廓、弹道特性和速度衰减的影响。随着水深减小，自由面与超空泡相互作用增强，由于大气的补气作用，通常卷气后空泡尺寸增加，稳定性增强。Chung[5]对静止流体中运动的近自由面通气超空泡形成进行了一系列实验研究，证明柔性自由表面对空泡尺寸的影响非常显著。中科院力学所[6-7]利用自行研制的SHPB高速弹射装置，模拟了圆柱射弹在水下不同深度处的航行状态，研究了射弹表面空化与自由界面的相互干扰特性。Dawson[8]研究了自由界面附近楔形水翼的超空化流动，获得了水翼的法向力、俯仰力矩、压心、空泡长度与水翼攻角、空化数、弗劳德数（Fr）和浸没深度的关系。在自由界面附近，空气易被卷入空泡进而形成超空泡，而当水翼远离水面时，则需要额外通气形成超空泡。

数值计算方面，初期基于线性和非线性的边界元方法是研究主流。挪威科技大学的Faltinse[9]利用非线性分析方法，提出了改进的闭合空泡尾流理论模型，推导了复杂流动势能函数，并用数值方法求解方程组，讨论了自由界面Fr对稳态水翼空化现象的影响。自由界面和重力场都导致空泡长度缩短，水翼升力系数与浸没深度成反比，与Fr成正比。相较于基于传统势流理论的方法，CFD 方法考虑了黏性，对复杂几何外形和复杂物理问题（如空化）更适用。中国科学院王一伟等[10]采用VOF 界面捕捉、Zwart 空化模型和LES 湍流模型，模拟了锥形、平头两种细长圆柱在自由界面下的非定常空泡演化过程，研究表明当空泡与自由界面距离过近以至于空气与空泡内气体掺混后，将导致空泡长度增加，并增加其稳定性。韩国釜山国立大学的Mei-Shan Jin[11]基于Mixture 多相流模型、Merkle 空化模型和k-ε湍流模型建立了多相流仿真程序，计算了自由界面附近的NACA0012 空化水翼通气空化流场，并证明其与已有试验数据的一致性，讨论了自由界面高度变化随通气速率的变化规律。Brizzolara[12]针对高速SWATH构型的水面航行体设计了新的超空化跨介质水翼，利用数值仿真方法分析了三维水翼附近的流场结构。研究表明自由界面锐化因子的取值，直接影响多组分气体掺混以及回射流现象的模拟，对水翼水动力预测的准确性有一定影响。张亚涛[13]研究了水深对空泡形态的影响，自由液面上出现的先导波浪高度随浸深增大而减小。程闻[14]采用数值水池技术，分析了波高和水翼浸深对水翼非稳态空泡的形态演化和脱落频率的影响，结果表明液面对脱落涡系有吸引作用。

本文针对轴对称和弹-翼结合近水面超空泡航行体（如图1 所示），开展空泡演化动力学机制和流体动力特性分析。图中空化水翼用于提供平衡重力的升力，并作为水上和水下部件的连接机构。水下轴对称部件采用通气超空泡技术实现高效减阻。首先基于OpenFOAM 平台构建多相空化流求解器，针对通气空泡和跨介质水翼算例开展求解器验证，进一步开展超空泡近水面航行器流体动力仿真分析，揭示自由界面影响下空泡形态瞬时演化规律，分析自由界面效应对空泡界面非对称性、空泡尺度、空泡尾部闭合特性与泄气机制的影响，分析跨介质水翼对气液界面波动和弹体表面压力分布影响规律，为近水面超空泡航行器外形优化设计提供支撑。

图1 近水面超空泡航行器示意图Fig.1 Schematic of the near free-surface supercavitating vehicle

1 多相流数值计算方法及验证

1.1 数学模型

1.1.1 控制方程

近自由面空化流涉及界面生成与空化相变、多组分气体掺混、多相湍流流动等复杂物理现象。本文基于均质多相流模型构建空化流仿真方法，对不同相求解同一速度场、压力场和温度场，不考虑由于相间速度差、温度差等导致的热力不平衡状态。采用VOF方法追踪空泡界面和自由界面。本文基于OpenFOAM开源平台构建多相不可压求解器。首先针对混合相求解质量和动量守恒方程。

混合介质的连续方程为：

式中：U表示混合物速度，下标m 表示混合相，局部混合物的密度 ρm和黏度系数 µm由各相体积分数αi加权：

式中：N为相的个数，i=l 或g，αl、αg分别表示液体和气体的体积分数。混合介质的动量方程：

式中：τ为剪应力，f为表面张力和重力引起的动量源项，表达式如下：

式中：Fσ为表面张力，g为重力加速度。连续表面力（CSF）方法将不连续的压降视为连续函数，满足界面条件：

式中：C为表面张力系数，κ为界面处的曲率，α为体积分数。

包含相变的多相系统各相守恒方程表示为：

式中：ρi和αi为第i相的密度和体积分数，m˙i定义为各相质量源项，且有：

压力-速度耦合由PIMPLE 算法处理，其基本思想是采用SIMPLE 算法将每个时间步长内按照稳态流动求解，采用PISO 算法时间推进，并根据流动特性，在求解过程中应用亚松弛方法实现高效计算。

1.1.2 VOF 方法

流体体积（VOF）法属于欧拉方法中的一类，用于处理自由表面流的界面跟踪问题，如溃坝、液体晃动等，是运用最为广泛基于均质平衡流假设的界面追踪方法。OpenFOAM 中多个多相求解器采用VOF 方法，如interFoam 类求解器（interFoam、compressibleInterFoam、interDyMFoam、multiPhaseInterFoam 等）。

相界面通过计算控制单元内的体积分数场而后处理得到。使用VOF 方法中为了保证相界面更加锐利，加入人工压缩项，对相界面附近的体积分数进行压缩。对不考虑相间质量传递的体积分数守恒方程，添加人工压缩项，得到：

式中第三项为人工压缩项，保证在α=0 或α=1 的非界面处为0。Uc为压缩速度，与界面法向同向：

式中，c为压缩因子，取值越大，压缩效应越明显。

1.1.3 湍流模型

本文采用k-ωSST 湍流模型来模化湍流方程，k-ωSST 模型结合了k-ε模型与k-ω模型的特点，且在边界层流动和空化流动中具有更好的性能。

以下给出湍动能k和比耗散率ω的输运方程：

其中：β*、β、σk、σω、γ 均为模型常数；为湍动能生成项；µt为湍流黏度；S为应变率；a1为模型常数；F1、F2为混合函数；σω2为模型常数；d为与壁面距离。

1.2 求解器应用及验证

1.2.1 通气空化算例

首先针对水洞中通气空化现象开展数值模拟，验证求解器精度。模型为头部带圆盘空化器的圆柱体，空化器直径15 mm，弹体直径10 mm，弹体长度100 mm，流速5 m/s，通过调整头部通气量和尾部喷管射流流量得到不同的空泡流型。图2 为计算域和边界条件，图3 为局部网格划分，图4 为空泡界面数值模拟与实验的对比。Qc和Qj分别为头部和尾部气体的体积流量，无量纲化后为：射流作用下通气空泡闭合于约2 倍模型长度处，尾部呈周期性震荡泄气模式。求解器能较好地捕捉空泡界面（如图5），空泡直径的相对位置误差小于9%，空泡长度误差小于3%，说明数值计算方法精度足够。

图2 计算域和边界条件Fig.2 Computational domain and boundary conditions

图3 局部网格划分Fig.3 Local mesh distributions

图4 数值模拟空泡界面与实验结果对比（Fr=13）Fig.4 Comparison between the predicted cavity surface with the experimental results (Fr=13)

图5 实验得到的通气超空泡轮廓与数值模拟结果对比Fig.5 Comparison between the predicted cavity profile with the experimental data

1.2.2 跨介质水翼算例

半浸式水翼在切割自由表面航行时，涉及自由面波动和空气卷吸等多种现象，因此应用半浸式水翼数值模拟开展求解器验证。水翼端面形状采用钝后缘、半椭圆形的Bi-ogival 对称翼型（如图6）。水翼弦长为279.4 mm，攻角α为10°，图7 为网格划分，图8为实验得到的垂直表面穿透水翼的尾部卷气与数值模拟结果的对比。弗劳德数Fr与展弦比AR特征长度均选取浸湿的展长。如图可知，水翼的压力侧出现大片的喷流幕，而卷气主要集中钝体尾部，形成倾斜的空腔，在分离的尾迹处产生大范围界面变形。数值模拟较为准确地复现了卷气位置、卷气空泡尺度与空气-水自由面曲率。

图6 Bi-ogival 翼型尺寸（单位：mm）Fig.6 Size for the Bi-ogival airfoil (Unit: mm)

图7 Bi-ogival 翼型网格划分Fig.7 Mesh distribution for the Bi-ogival airfoil

图8 垂直表面穿透水翼的尾部卷气（α=10°,Fr=2,AR=1.0）Fig.8 Tail air entrainment for the vertical surface piercing wing(α=10°,Fr=2, AR=1.0)

1.3 近自由面超空泡航行体构型及网格划分

近自由面超空泡航行体由升力体艇身、跨介质水翼、超空泡水下部件、空气引流系统、测控系统和动力系统等组成，其流体动力特性主要受跨介质水翼和轴对称水下部件影响。本文首先针对单独轴对称体开展流体动力特性分析，进一步针对水翼-轴对称体复合构型（简称为翼-弹组合体），研究水翼与轴对称体耦合效应。

轴对称水下部件设计为流线型，头部放置圆盘空化器，轴对称体直径为70 mm，空化器直径23 mm。计算网格由ICEM 生成，在航行体壁面附近及自由面附近加密，全部为六面体结构网格，总网格量约为223 万，如图9 所示，其计算域及边界条件设置如图10 所示。

图9 轴对称航行体附近结构网格划分Fig.9 The structural mesh distribution near the axisymmetricvehicle surface

图10 航行体计算域及边界条件Fig.10 Computational domain and boundary conditions for the axisymmetric vehicle

跨介质水翼整体外形为梯形，根梢比为2，由低速翼型与平板翼复合设计而成。水翼截面前后缘均采用低速NACA66 翼型。水翼尾部的半个NACA66翼攻角为5°，水翼下挂轴对称航行体，下反角为45°。复合体结构和网格分布如图11 所示，其计算域及边界条件设置如图12 所示。

图12 复合体计算边界条件Fig.12 Computational boundary conditionsfor the combined body

2 自由界面与通气空泡耦合作用机理研究

2.1 通气空化流与波动自由面耦合研究

针对不同弗劳德数（Fr）、空化数（σc）和通气流量系数（Cq）开展CFD 计算。引入无量纲参数来描述超空泡性能，特别地，引入浸深比（由浸没深度H除以空化器直径D得到）描述自由表面的影响：

其中：p∞为环境压力，pc为空泡内压力，ρ∞和U∞分别为来流密度和速度，Q为通气质量流量。

图13 对比了不同弹体长径比对气液界面形态的影响。由图可知，近自由面航行体导致自由面出现显著的兴波效应，最大波幅位置接近航行体中段。自由面影响下超空泡形态具有明显不对称性，上表面受到自由面“吸引”导致曲率增大，下表面相对比较平缓。原因是自由面和空泡表面间形成局部狭窄通道，导致液体流速增加，压力减小，从而空泡向自由面方向扩张。图13（a）中工况自由面最大波幅与空泡径向尺度相当，而空泡上表面最大直径，即上表面与弹体中轴线的最大距离，相对于下表面增加了约30%。相同通气和航行工况下，航行体长径比增加时，空泡闭合位置向后移动，空泡上表面曲率减小，且空化兴波波长增加，波幅减小。因此近自由面空泡形态与弹体外形密切相关。由图中流场结构分布可知，小长径比航行体在艉部产生流动分离和小范围的回流区，将导致压差阻力增加。图13 两种工况压差阻力系数分别为0.075 4 和0.069 6，大长径比弹体相对于小长径比弹体压差阻力系数减小了7.69%。对于近自由面航行体，在形成超空泡状态下提高航行体长细比对空泡稳定性和减阻性能有利，但是过长的航行体表面更容易发生沾湿。因此航行体外形必须与考虑自由面影响下的通气空泡流型进行匹配设计。

图13 内置航行体长度对通气空泡内部流动的影响（Cq=9.259,Fr=31.37,=3.48）Fig.13 Effect of the inside body length on the internal flow structure of the ventilated cavity (Cq=9.259,Fr=31.37,=3.48)

图14 对比了相同Fr和航行深度、不同通气量下的空泡界面演化图。随着通气量增加，空泡尺度增大，自由界面波动范围和幅度均显著增加。因此在生成超空泡的前提下，应尽量减小空泡最大直径，减小和自由界面之间的扰动。图15 比较分析了Fr对气液界面形态的影响。Fr=10.53 时，由于浮力影响显著，空泡呈现出强非对称性，航行体尾段1/3 处完全沾湿，尾流与自由面产生强烈的相互作用，自由面最大抬升位置在主空泡后方。Fr=62.75 时，形成稳定双涡管泄气超空泡，且由于空泡尺度增加，对自由面的扰动增强，自由面最大波幅位置向弹体尾端移动，且由于自由面作用，空泡中轴线向下偏转，与无穷水域中空泡形态存在显著差异。原因是高Fr状态下，随着自由界面波峰向后移动，空泡尾部压力显著增加，高压作用下空泡尾部向下偏转。

图14 不同通气量下气液界面形态对比Fig.14 Comparison of the gas-liquid interface under different ventilation rate

图15 不同弗劳德数下气液界面形态对比Fig.15 Comparison of the gas-liquid interface under different Fr number

表1 给出了在相同Fr和Cq下，不同浸深比下的摩擦阻力系数和压差阻力系数，可以看到两者的差别较小。小浸深比下摩擦阻力和压差阻力均较小，主要是由于接近自由面导致空化数减小，泡内压力增大，使得压差阻力减小，而回流区尺度增加也导致摩擦阻力有所减小。然而过于靠近自由面势必导致界面波动增强，因此存在最佳工作深度。

表1 不同浸深比下的阻力系数Table 1 Drag coefficient under different immersion depth

表2 对比了相同浸深比和通气系数下、Fr对阻力系数的影响。小Fr下出现大面积沾湿，压差阻力系数和摩擦阻力系数均大幅增加，Fr=62.75 时，空泡压差阻力系数较小。因此近自由面航速不宜过小，以避免浮力导致壁面沾湿。

表2 不同弗劳德数下的阻力系数Table 2 Drag coefficient under different Fr number

2.2 翼-弹组合体耦合效应分析

针对翼-弹组合体开展流体动力特性分析。如图16，在Cq=9.259、Fr=31.37、=3.48 时，连接跨介质水翼后自由表面出现较大起伏，总影响范围达到3 倍航行体长度区域，最大自由面抬升分别出现在主空泡兴波区与压力面对应的尾流，峰值达到1.74D。空化兴波使得水翼沾湿面积增加，水翼表面出现小范围的片状射流；最大自由面凹陷出现在吸力面尾流，达到-2.17D。整个尾流区是水翼尾流、通气空泡尾流、各相湍流相互作用的复杂区域，非定常特性显著。

图16 翼-弹组合体通气空泡流场自由界面起伏（Cq=9.259,Fr=31.37,=3.48）Fig.16 Free-surface fluctuation induced by the ventilated cavitating flow around the hydrofoil-projectile combined body(Cq=9.259,Fr=31.37,=3.48)

在航行体对称轴处建立局部柱坐标系，定义吸力侧和压力侧如图17 所示。图中云图表示横截面处压力分布，蓝色实线表示模型截面线，黑色实线表示空气体积分数等值线。跨介质水翼导致的压力分布变化引起空泡形态非对称性加剧，主体空泡流型并没有被破坏。对整个组合体，空化器前缘驻点是压力最高点。整个水翼压力分布几乎对称，水翼前后缘沾湿部分压力明显高于周围压力。

图17 水翼-航行体复合构件横截面处的通气空化流场压力分布及界面特征（Cq=9.259,Fr=31.37,=3.48）Fig.17 Pressure distributions and the interface characteristics on the cross sections induced by the ventilated cavitating flow around the hydrofoil-vehicle combined body(Cq=9.259,Fr=31.37,=3.48)

在截面1 至截面3 处未受水翼影响，前段主体空泡形态与2.1 节中类似，空泡非对称主要由重力效应引起，并且由于自由面抬升较少，自由面对空泡影响不显著。

在截面4 至截面6 处，连接水翼显著改变了主空泡流场压力分布，吸力面处低压导致翼-弹连接处的吸力侧通气空泡直径增大、压力侧空泡直径减小。

受水翼穿透影响，在截面7 至截面9 处，压力侧部分主空泡受压，向内凹陷。水翼切割了主空泡尾部，使得空泡没有完整包裹航行体，水翼下游出现壁面的间歇沾湿，称为水翼穿刺空泡，空泡特征与航行体控制面穿刺空泡的特征类似。主空泡在压力侧受压明显，在截面9 处出现压力最高点，由于后缘水翼存在攻角，吸力侧和压力侧两边的压差达到最大值，提供了主要升力。此时发生部分卷气，卷气方向沿翼展向下。这种卷气机制属于吸力面诱导卷气，不同的是，由于没有稳定的分离空间，因此只形成一个凹陷的空腔，空腔深度没有达到水翼的整个浸入展长，未出现不稳定的回射流，因此不满足维持完全卷气流动的稳定条件。

在截面10 至截面12 处，自由面最大起伏出现在水翼尾流处，水翼尾流与船体尾流类似，使得下游区域液面下沉，如果航行体尾部安装推进装置，会导致推进效率下降。如果航行体长度过长使控制舵面位于尾流卷气位置，可能造成上、下控制面沾湿面积不同，会引起控制性能急剧变化，设计时应当考虑卷气和自由面波动的综合影响，合理调整浸没深度。

相比于单独水下航行体，连接水翼的复合构件主空泡上方波幅最高降低35%，如图18 所示，吸力面处的卷气导致主波位置压力和波幅下降，压力面尾迹区产生次波。由于主波对空泡影响显著，因此连接水翼有减少主波扰动和增加主空泡稳定性的效果。

图18 连接水翼对对通气空泡主波的影响（Cq=9.259,Fr=31.37,=3.48）Fig.18 Effect of connecting with hydrofoil on the primary wave induced by the ventilated cavity(Cq=9.259,Fr=31.37,=3.48)

由图19 的流线分布可知，流线由水翼后缘的压力侧向吸力侧收缩汇集，空泡界面附近的流动存在周向流动。由于航行体连接水翼，使空泡后半段周围压力分布不均，高压流动循着空泡界面由下往上，滚卷流动到较低压的吸力侧，与流体向后直线流动的速度矢量相加，最终形成螺旋式的运动，类似于水翼或机翼产生的翼尖涡。对称航行体通气空泡的双涡管由一对大小相等，方向相反的环量产生。相比于控制面穿刺空泡，水翼穿刺空泡对空泡周围压力分布改变更明显，速度分布同样发生变化，产生额外环量，方向与双涡管其中一个涡方向相反，与另一个方向相同，导致两涡强度和大小非对称，影响双涡管生成，空泡泄气机制发生改变。对于单个水翼-航行体复合部件，航行体通气空泡尾涡与水翼作用产生大小不同、方向相反的横滚力矩；对于整个近水面航行器，镜像对称的水翼和航行体产生的力矩相互平衡。

图19 翼-弹组合体通气空泡流线图（Cq=9.259,Fr=31.37,=3.48）Fig.19 Streamline distributions for the ventilated cavitating flow around the hydrofoil-projectile combined body(Cq=9.259,Fr=31.37,=3.48)

由于跨介质水翼与通气空泡相互作用，呈现出后缘的压力侧的间歇沾湿现象，如图20 所示。在相同浸深比、相同通气系数和不同Fr下，航行体艏部至水翼平板翼周围的通气空泡无明显差异，而穿刺空泡的尾部沾湿存在显著区别。各个工况下吸力侧周围完全浸湿，而为了提供主要升力，水翼尾部的半个NACA66翼型存在攻角，压力侧呈现出相似的沾湿面。随着Fr增加，沾湿面积增大，因此高速航行时的减阻性能可能受限。

图20 不同Fr 下水翼穿刺空泡的尾部沾湿对比（Cq=9.259,=3.48）Fig.20 Comparison of the tail wet region for the surfacepiercing hydrofoil under different Fr number(Cq=9.259,=3.48)