三维非定常水下超空泡射弹大涡模拟研究

张 木，谭俊杰，易文俊，蔡晓伟

（南京理工大学1.能源与动力工程学院；2.瞬态物理重点实验室，南京210094）

空化是流动中的常见现象，往往能导致空蚀、噪声和机械效率下降等负面现象的产生，因此，在早期工程设计中，人们对空化现象的分析与研究主要是为了避免空化现象的产生[1].但随着俄罗斯超空泡鱼雷“暴风”的问世，人们发现利用超空泡技术可使水下航行体被包裹在空泡内，仅头部与水体接触，从而减少航行体在水中的表面摩擦阻力与压差阻力，提高水下航行体的运动速度和航行距离[2，3].因此，国内外许多学者在空泡流动的数值模拟方面开展了大量的研究工作，如KINNAS和FINE[4]发展了一种基于速度势的边界元方法，用于二维水翼的局部和超空化流动分析；文献[5]采用回射流模型对三维水翼的局部空泡进行了非线性分析；文献[6]基于均相N－S方程，采用高雷诺数k－ξ湍流模型，对空化数为0.15条件下绕三维尖头回转体的空泡流动进行了数值模拟，其结果与实验数据吻合较好；文献[7]基于simple算法，采用带低雷诺数修正k－ξ湍流模型，模拟了在空化数为0.01条件下的绕三维轴对称航行体的空泡流动，得到航行体运动阻力系数与空化数之间的变化关系；文献[8]利用CFD商业软件对空化器线形与超空泡形态的关系及其阻力特性进行了数值仿真，分析了影响超空泡航行体减阻效果的因素.以上这些数值仿真研究主要针对空化数大于0.01的条件下，利用雷诺平均湍流模型模拟水下航行体的空泡流动.

针对此种情况，本文基于N－S方程及PISO算法，引入考虑气液间相变的空化模型，采用大涡模拟湍流模型和VOF方法，在空化数为0.78×10－3～4.88×10－3的条件下对绕三维水下自然超空泡射弹的空泡流动进行了数值模拟.

1 控制方程

1.1 大涡模拟

为求解非定常超空泡流动，近似认为气液两相流动是等温均相流动.另外根据文献[9]的研究，在等温均相流动的前提下，空化流动中气相与液相的可压缩性影响很小，而气液混合区域，其可压缩性影响需要考虑.因此，对于气液混合区域，建立速度、压力和密度的耦合关系式，考虑可压缩性的影响.对于气液两相采用VOF方法捕捉气液两相间的界面，以及采用大涡模拟湍流模型对流场的湍流运动进行模拟.关于气液两相之间相变问题，本文采用KUNZ[6]空化模型进行模拟.基于以上分析，控制方程可写为

式中，φl为液相的体积分数，ρl为液相密度，δij为Kronecker函数，μeff为有效动力粘性系数为相变质量传递源项，ρm为气液两相的混合密度为滤波速度，为滤波压力，ksgs为亚网格湍动能，g为重力.变形张量1，2，3.为了求解μeff＝μm＋μsgs，需要求解μm和ksgs.由于ρm＝φlρl＋（1－φl）ρv，ρv为气相密度，因此μm＝φlμl＋（1－φl）μv.关于ksgs求解，采用k－Δ模型进行封闭，其方程可以写为

1.2 空化模型

式中，为滤波压力，pv为空化压力，t∞为特征时间，u∞为来流速度，Cp和Cd为经验系数，本文中Cp和Cd分别为1 000和800，t∞＝L∞/u∞，L∞为特征长度.

2 离散方法

本文采用有限体积法对控制方程进行离散，其中时间项采用二阶精度的三点向后隐式差分格式，对流项采用二阶精度的线性迎风格式，扩散项采用二阶中心差分格式.另外基于PISO算法，建立压力、密度与速度的耦合算法，考虑气液混合相可压缩性的影响.在预测步骤中，将上一时刻的压力分布作为当前时刻动量方程中的压力分布p＊，然后求解如下方程：

式（8）中方程右端项为混合密度在p0点处的一阶泰勒展开式，而混合密度可通过式（1）、式（2）以及混合密度表达式导出.通过式（8）得出压力修正量p′，进而计算出速度修正量u′.之后利用关系式u＊＊＝u＊＋u′和p＊＊＝p＊＋p′，求解出速度与压力的修正值u＊＊和p＊＊.在第二校正步中，将速度与压力的修正值代入第二次压力修正方程得：

通过式（9）得出第二次压力修正量p″，进而计算出第二次速度修正量u″，然后根据关系式u＊＊＊＝u＊＊＋u″和p＊＊＊＝p＊＊＋p″，求解出改进的速度值和压力值.

3 计算结果与分析

采用图1所示的圆盘空化器多锥体射弹模型，进行了自然超空泡流动数值计算.

图1 圆盘空化器多锥体超空泡射弹模型

在图1中，空化器直径为3.5mm，射弹长径比为1∶10.计算域采用图2（a）所示的分块结构化网格对超空泡射弹进行三维网格建模，网格总数为1 602 000.为了保证1＜y＋＜3，第一层网格点距离壁面为6.5μm.在非定常计算中，为获得较好计算结果，通常以较小时间步长计算，文中计算时间步长为4.3ns，而且当流场变量的残差降到5×10－4以下时，认为计算过程收敛.由于网格数较大及时间步长较小，采用MPI并行算法提高计算效率.

超空泡流场计算边界条件如图2（b）所示.计算域顶部、底部与左端均为速度入口，右端为压力出口，射弹表面采用固壁边界条件.关于物性参数的选取，水的密度为1 000kg/m3，动力粘性系数为1×10－6m2/s；水蒸气的密度为23.08g/m3，动力粘性系数为4.27×10－4m2/s.定义无量纲量空化数σ＝（p∞－pv）/（0.5ρu2∞），其中，p∞为来流压力，pv为空化压力，ρ为水的密度，u∞为来流速度.在自然超空泡条件下，本文通过变化来流速度来调整空化数.

图2 计算网格与边界条件

3.1 超空泡形成过程

在空化数σ＝0.002 17的条件下，计算了绕超空泡射弹的空泡流动，并给出了不同时刻下超空泡形态的计算结果，如图3所示.

由图3可以看出，在图3（a）～图3（d）阶段，射弹头部空泡与尾部空泡随时间的增加而不断变大、变长；在图3（d）～图3（g）阶段，射弹头部形成的空泡与尾部空泡相连接，并继续向远处发展，形成完全包裹弹身的超空泡；在图3（h）阶段，在射弹表面出现稳定的空化区域，尾部被一个拖拽的巨大气泡包围着，形成一种稳态的开式尾流超空泡流型.

3.2 空化数对超空泡形态的影响

定义超空泡最大长度与空化器直径的比值c＝Lc/Dn，超空泡最大直径与空化器直径的比值c＝Dc/Dn.

为研究空化数对超空泡形态特性的影响，在空化数分别为4.88×10－3，3.12×10－3，2.17×10－3，1.59×10－3，1.22×10－3，0.96×10－3，0.78×10－3的条件下，对无量纲化空泡特征长度c和无量纲化空泡特征直径c进行了计算，计算结果如图4所示.

图4 超空泡特征长度与直径随空化数变化关系

由图4可以看出，超空泡特征长度与特征直径随空化数的增大而减小.为了验证计算结果的正确性，将其分别与LONGVINOVICH经验公式计算结果及SAVCHENKO经验公式计算结果进行对比，结果如图4所示.由对比可知，关于特征直径，本文数值计算结果与2种经验公式吻合较好；关于特征长度，本文计算结果与2种经验公式变化趋势大体一致，但本文计算结果与LONGVINOVICH[10]经验公式吻合较好，而与SAVCHENKO[11]经验公式的吻合较差，这是由于SAVCHENKO经验公式主要是适用于空化数在0.012～0.057的范围内，而LONGVINOVICH经验公式是适用于空化数在0～0.25的范围内.

3.3 攻角对超空泡形态的影响

在空化数为0.002 17的条件下，计算了攻角α分别为1°，3°，5°时水下超空泡射弹的空泡流动，计算结果如图5所示.

图5 不同攻角下超空泡射弹的空泡形态

从图5中可以看出，在有攻角的条件下，射弹表面上、下空泡厚度出现不对称性，且随着攻角的增加，其不对称性增强.为了分析攻角对空泡形态的影响规律，给出了不同攻角下的射弹上、下表面空泡厚度分布曲线，如图6所示，图中hu，hd分别为上、下表面空泡厚度.

从图6（a）可以看出，在攻角为0°和1°的情况下，射弹上表面的空泡厚度从射弹前缘起先逐渐增加，随后变小，再逐渐变大；在攻角为3°和5°的情况下，射弹上表面的空泡厚度从射弹前缘起逐渐变大；另外随着攻角的增加，射弹上表面的空泡厚度不断变大.从图6（b）可以看出，在攻角为0°和1°的情况下，射弹下表面的空泡厚度从射弹前缘起先逐渐增加，随后逐渐变小，再逐渐变大；在攻角为3°和5°的情况下，射弹下表面的空泡厚度由前缘起略有增大后逐渐变小，甚至在射弹的下表面，弹体有部分处于粘湿状态，之后再略有增加.综合上述分析可知，攻角对于超空泡射弹的空泡形态影响很大，特别是当攻角较大时，射弹有部分处于粘湿状态，这将影响超空泡对于射弹的减阻效果及射弹的水下弹道.

图6 不同攻角下射弹上、下表面空泡厚度分布曲线

4 结论

本文基于PISO算法及N－S方程，利用LES湍流模型和VOF方法，实现三维非定常水下超空泡射弹形态特性的数值模拟，得出如下结论：

①通过分析超空泡形成过程以及对比计算结果与经验公式，验证了本文计算结果的正确性；

②研究了不同空化数对超空泡形态的影响，计算结果表明，当空化器直径一定时，超空泡特征长度和特征直径随空化数的增加而减小；

③射弹航行攻角对水下超空泡形态特性影响较大，其显著改变超空泡形态的轴对称性.

[1]FAROUK M，NAYFEH H.A compressible multi－phase flow solver for the computation of the supercavitation over highspeed torpedo[J].AIAA，2002，14（17）：1－14.

[2]LOGVINOVICH G V.Some problems of supercavitating flows[C].Proceedings of NATO－AGARD.Ukraine：NAS－IHM，1997：1－5.

[3]SAVCHENKO Y N，SEMENENKO V N.The gas absorption into supercavity from liquid－gas bubble mixture[C].Proceedings of 3rd International Symposium on Cavitation.Grenoble，France：Universite Joseph Fourier，1998：49－53.

[4]KINNAS S A，FINE N E.A numerical nonlinear analysis of the flow around two－and three－dimensional partially cavitating hydrofoils[J].J Fluid Mech，1993，254：151－181.

[5]DANG J，KUIPER G.Re－entrant jet modeling of partial cavity flow on three－dimensional hydrofoils[J].Journal of Fluids Engineering，1999，121：781－787.

[6]KUNZ R F，BOGER D A，STINEBRING D R，et al.A preconditioned Navier－Stokes method for two－phase flows with application to cavitation prediction[J].Computers ＆ Fluids，2000，29：849－875.

[7]陈瑛，鲁传敬，吴磊.三维小空化数空泡流数值方法[J].计算物理，2008，25（2）：163－171.CHEN Ying，LU Chuan－jing，WU Lei.Numerical method for three－dimensional cavitation flows at small cavitation numbers[J].Chinese Journal of Computational Physics，2008，25（2）：163－171.（in Chinese）

[8]贾力平，张嘉钟，于开平.空化器线性与超空泡减阻效果关系研究[J].船舶工程，2006，28（2）：20－23.JIA Li－ping，ZHANG Jia－zhong，YU Kai－ping.Research on the relation of cavitator profile and supercavity drag reduction[J].Ship Engineering，2006，28（2）：20－23.（in Chinese）

[9]QIAO Qin.A virtual single－phase natural cavitation model and its application to CAV2003 hydrofoil[C].5th International Symposium on Cavitation.Japan：Nova Science，2003：1－8.

[10]LOGVINOVICH G V.Hydrodynamics of flow with free boundaries[M].Kiev：Naukova Dumka，1969.

[11]SAVCHENKO Y N，VLASENKO Y D，SEMENENKO V N.Experimental study of high－speed cavitation flow[J].Hydromechanics，1998，（72）：103－111.