绕过带凹槽的圆盘空化器的超空泡流的研究

，，

(浙江理工大学机械与自动控制学院，杭州 310018)

0 引 言

航行体在水下高速运动时，头部空化器能够诱导产生超空泡，包裹整个航行体，大幅度降低航行体的阻力。但水下航行体超空泡受到扰动流场时，容易溃灭，极其不稳定。Savchenko等[1]提出一种带有浅凹槽的、多边形的空化器，该空化器能够产生稳定力矩，且具有抗扰动的作用。在此基础上，施红辉等[2]研究了水下航行体表面开槽数量对超空泡流场影响，发现随着航行体表面的凹槽数的增加，超空泡的最大无量纲直径和长度会减小，但无法确定航行体表面开槽对无量纲尺寸的影响，也未对开槽航行体超空泡水动力学特性进行分析等。所以针对超空泡航行体的主要沾湿部位——头部空化器，进行了深入的研究，能够达到超空泡的最大减阻效果。超空泡减阻增速的优势，已经被美俄等国用于高速水中兵器的设计中，但我国水下超空泡流动的研究仍亟待提高[3-6]，只有夯实超空泡流动技术，才能进入空(气)中到水下连续机动的超空泡兵器开发阶段。

超空泡流动是通过空化器的脱体绕流发生的，因此空化器的几何形状对超空泡流动有着十分重要的影响[7]。先前研究主要针对圆盘、圆锥、矩形板、卵巢型等空化器形状的空化器，其中圆盘型空化器是应用最为广泛[8]。针对以上空化器形状，国内外许多专家学者也进行了许多实验、理论和数值模拟方面的研究。如罗格维诺维奇等[9]在20世纪60年代通过理论和实验研究给出了圆盘空化器阻力系数、空化器直径和空化数与自然超空泡形态的函数关系。Savchenko[10]在实验基础上研究了锥形空化器，通过外推法给出了空化数趋近于零时的阻力系数，并与近似公式进行了比较。栗夫园等[11]在水洞试验中研究了锥形空化器，分析了锥形空化器的锥角对航行体流体动力学特性的影响，对圆盘空化器超空泡流体力学特性有待进一步分析。徐清沐等[12]基于雷诺平均N-S方程和均质平衡流理论对带凹槽和凸起的圆盘空化器进行了数值模拟，并给出了空化数对空泡凹凸形态和阻力系数的影响。金大桥等[13]分析了水下航行体结构参数对自然超空泡减阻特性的影响，获得了超空泡航行体的阻力系数变化趋势。牟斌等[14]利用空化流动隐式求解方法对三维尖锥空化器形成的自然空泡的流场特征进行研究。

为了进一步研究超空泡航行体最大限度的减阻机理，本文选取不同凹槽数的圆盘空化器航行体为研究对象，应用商业软件Fluent对其所形成的三维超空泡进行数值模拟，验证三维超空泡的自恢复特性，还分析了空化器凹槽数对航行体流体动力特性的影响。

1 数值模拟方法

1.1 控制方程

对超空泡射弹的数值模拟是基于连续性方程和动量方程，以及湍流模型采用标准k-ε模型，多相流模型采用VOF模型，空化模型采用Schnerr-Sauer空化模型。

空化数表征超空泡的最为重要的无量纲参数之一，其定义为：

(1)

其中：P∞表示境的压强；P表示空泡内压强；ρ表示流体密度；V∞表示自由来流速度或者是航行体的航行速度。连续方程为：

(2)

其中：t表示时间；xj(j=1，2，3)表示笛卡尔坐标；uj表示笛卡尔坐标系下的速度分量；ρ表示混合介质的密度。动量方程为：

(3)

其中：ui和fi表示在i方向的速度和质量力；Pm表示混合压力；μ表示混合相的动力粘度；ρ表示混合相的密度，被定义为：

ρ=αvρv+αgρg+(1-αv-αg)ρl

(4)

其中：α表示某一相的体积分数；v，g，l分别表示蒸汽相、空气相和液态水相。

标准k-ε模型湍流模型是两方程湍流模型中最具代表性的模型，也被认为是工程应用中最为普遍的模式[15]。采用湍流动能k和湍流动能耗散率ε来表示涡粘性系数vt，如式(5)所示：

(5)

k的输运方程为：

(6)

模化后的ε输运方程为：

(7)

其中各个经验系数取值分别为：Cε=0.07～0.09;Cε1=1.41～1.45;Cε2=1.91～1.92．

采用Schnerr-Sauer空化模型[16]直接从汽/液净质量传输率的表达式为切入口，对表达式中的体积分数项进行了计算，最后得到相变率的表达式如下：

(8)

(9)

(10)

其中：Me、Mc和RB分别表示冷凝的质量、蒸发的质量和气泡的半径。

1.2 计算域、网格划分及边界条件

数值模拟采用复杂表面形状带有圆盘形空化器的航行体，以Yang等[17]的航行体为原型，又根据试验中高压气枪的口径为参考设计，总长为48 mm，最大直径为6 mm，空化器是直径5 mm、厚度1 mm的圆盘；在空化器边缘处开槽数分别为0、2、4、6，凹槽为1 mm×1 mm×1 mm的正方体，如图1所示；航行体模型的几何尺寸参数如图2所示。

图1 航行体模型示意

图2 航行体的几何尺寸参数示意(长度单位：mm)

计算域是一个400 mm×200 mm×200 mm的长方形水箱，如图3所示。航行体头部空化器到前方来流的速度入口的距离为80 mm，位于计算域中间处。航行体的速度为100 m/s，压力值为一个标准大气压，四周壁面设为无滑移壁面边界条件，近壁面处采用标准壁面函数，见图3。通过Fluent中的Patch功能，可以设置计算域的初始相。

图3 计算区域边界条件设置示意

网格的划分如图4所示，网格全部采用六面体结构网格划分，为了能更清晰地显示出超空泡壁面处的汽液交界面，在航行体模型的四周进行了局部网格加密，最终得到的0、2、4、6道槽航行体计算域网格的节点数分别为727812、706307、774576、896792个。在数值计算中采用VOF的多相流方法捕获汽液交界面，并且采用标准的k-ε湍流模型，与连续性方程、动量方程等形成封闭的方程组，进行数值模拟的求解。此外使用压力基求解器，流动为非定常流动，流动中设置空气、水、水蒸气三相，并考虑重力的影响。

图4 计算域网格划分示意

通过改变航行体头部圆盘空化器上的凹槽数量，逐个进行数值模拟的计算，分析不同凹槽数空化器的航行体对超空泡流场的影响，并进行了水动力学分析。不同凹槽数空化器航行体网格划分情况示意图如图5所示。

图5 不同凹槽数空化器航行体网格划分情况示意

1.3 网格无关性验证

选取0道槽航行体的计算域进行网格划分，得到节点数分别为50万、70万和100万左右三种密度的网格，分别进行数值模拟计算。将计算结果进行后处理分析分析，绘制不同网格密度下航行体阻力系数Cd随时间变化的曲线图，如图6所示。

图6 不同网格密度下航行体的阻力系数随时间的变化曲线

对图6分析可知，中等程度的网格质量(70万左右)与较高程度的网格质量(100万左右)的计算精度相差不大，最大偏差为4.9%，最小偏差为1.2%，故选取节点数为70万左右的网格能够满足计算的要求。

2 计算结果分析

2.1 计算结果

头部空化器开0、2、4、6道凹槽的航行体在100 mm水深处诱导产生的超空泡数值模拟的三维和截面超空泡的水相图如图7—图10所示。每幅图由两部分构成，分别是超空泡三维形态图和超空泡纵向截面图，每一幅表示的时刻T依次为0.2、0.4、0.8、1.2、1.6、2.0、2.4 ms。

图7 0道槽空化器航行体数值模拟 超空泡的三维和截面水相图

图8 2道槽空化器航行体数值模拟 超空泡的三维和截面水相图

图9 4道槽空化器航行体数值模拟 超空泡的三维和截面水相图

图10 6道槽空化器航行体数值模拟 超空泡的三维和截面水相图

从图7—图10可知，当T=0.2 ms时，在航行体的空化器处、肩部以及尾部同时形成局部空泡；当T=0.4 ms时，空化器、肩部以及尾部空泡逐渐融为一体；当T=0.8 ms时，形成一个完全包裹航行体的超空泡。综合超空泡的截面水相图来看：在速度相同的情况下，不同凹槽数空化器航行体的超空泡的初生、发展过程比较接近。但是从空泡的三维形态图，可以看出，不同的凹槽数空化器航行体形成的超空泡会根据相应的凹槽数形成凸凹形状。然而，随着空泡长度的增加，凸凹形状就会很快被修复成光滑椭球形，超空泡所具有的这种特性被称作超空泡的自修复特性。

图11给出不同凹槽航行体的距离圆盘空化器的距离为1 mm超空泡的横截面图；图12给出了4、6道槽航行体的距离圆盘空化器的距离为1、3、5、7、9、11 mm处超空泡的横截面图。为了节约篇幅，在最大空泡直径Dc处(约在80 mm处)的结果没有给出。

图11 在x=1 mm处超空泡横截面水相图

图12 不同凹槽数航行体超空泡不同 位置处横截面水相图

将图11和图12对比分析后可知：a) 圆盘空化器有诱导产生超空泡的作用，且在超空泡的初生阶段，其截面形状会根据凹槽数的不同，形成不同凸凹形状的超空泡，但最终都变成圆形。b) 图12给出了4、6道槽航行体超空泡横截面从四角形、六角星形变化为圆形的过程。可知，在较短距离范围内，超空泡的横截面发生了较为显著的变化，体现出超空泡有很强的自修复能力，超空泡所具有的这种自动修复自身形态，我们称之为超空泡的自修复特性。这与超空泡所具有的自动缝合空洞的能力相似，均是由于超空泡处于高度非平衡状态，具有很大的熵差才驱动了这种自修复特性[18]。

图13给出了超空泡横截面自修复成圆形时所处的无量纲位置x/Dn，其中横坐标代表不同的凹槽数。图中对比分析了文献[1，19]中航行体超空泡的自修复特性结果，文献中航行体速度为40 m/s，在航行体表面开槽，开槽数同上。由图13可知，尽管计算工况的速度不同，但都出现了超空泡的自修复现象。

图13 不同凹槽数和速度的超空泡自修复曲线注：三角代表空化器上开槽，速度为100 m/s的工况；圆代表航行体表面开槽，速度为40 m/s的工况。

2.2 与半经验公式的比较

Logvinovich通过动量定理结合大量的实验数据，对超空泡的外形提出如下半经验公式[6]：

(11)

(12)

其中：Cx0为当空化数σ=0时的阻力系数，对于圆盘空化器而言，Cx0一般取值为0.827；k是一个经验常数，取值为0.96～1.00。根据超空泡独立膨胀原理得出的超空泡的外形轮廓的公式：

(13)

随机选取t=1.2 ms到t=2.1 ms这两个时刻的数值模拟结果，并与式(13)进行对比分析，得到超空泡外形轮廓的对比曲线图，如图14所示。

图14 0道凹槽航行体超空泡外形轮廓与 Logvinovich公式对比

分析图14可知：上述两个时刻的超空泡前沿外形轮廓与Logvinovich公式得出的轮廓曲线具有较好的一致性，这表明了数值模拟的准确性；而后沿是由于计算域限制，数值模拟并未达到超空泡的最大直径和长度，但从选取的两个时刻可以看出，数值模拟结果与经验公式结果是一个逐渐趋向的过程，在足够大的计算域中的某一时刻，数值模拟结果将会与半经验公式完全重合。

3 航行体阻力系数分析

航行体超空泡在水下航行时，所受到的阻力F主要包括压差阻力Fp和粘性阻力Ff两部分，压差阻力主要取决于航行体头部与航行体尾部的流场压力，而粘性阻力与航行体表面积及其沾湿介质的密度有关[20]。综上所述，本文中考虑了压差阻力系数、粘性阻力系数以及总阻力系数。如下式：

(14)

(15)

(16)

其中：S为特征面积，这里取圆柱段的横截面积。式(14)—(16)中的压差阻力、粘性阻力和总阻力可以在通过对数值模拟的后处理计算求得。

图15—图17是超空泡航行体的压差阻力系数、粘性阻力系数以及总阻力系数随时间的变化曲线。

图15 航行体在超空泡形成过程中的压差阻力系数变化

图16 航行体在超空泡形成过程中的粘性阻力系数变化

图17 航行体在超空泡形成过程中的总阻力系数变化

根据图15—图17，并结合超空泡的形成、发展的水相图(图7—图10)可知：a) 在T=0～0.6 ms时间内，航行体的阻力系数急剧下降；这是由于超空泡逐渐形成，航行体的头部与尾部压差以及沾湿面积都急剧下降，致使航行体所受到的阻力急剧降低的缘故。b) 在T=0.6～2.4 ms时间内，航行体的阻力系数趋于一条平滑的直线；这是因为在T=0.6 ms时间以后，超空泡完全包裹航行体，航行体的头部与尾部压差以及沾湿面积都保持不变，故其阻力系数均是直线。c)对单幅图来分析可知，对于不同的凹槽数航行体而言，随着凹槽数的增加，航行体的压差阻力、粘性阻力和总阻力系数均呈减小趋势。

4 结 论

通过对不同凹槽数的圆盘空化器航行超空泡的三维流场的数值研究，将计算得到的数据进行后处理分析、与Logvinovich半经验公式进行对比以及对其流体动力学进行分析，得出如下结论：

a) 随着航行体空化器的凹槽数的增加，超空泡无量纲直径和长度呈现一定的减小趋势。

b) 尽管不同凹槽数的空化器造成了非轴对称的三维超空泡形状，但在沿流向不长的距离处，超空泡能够很快自动修复形成光滑椭球形超空泡，这即是超空泡的自修复特性。

c) 随着空化器凹槽数的增加，其受到的压差阻力和粘性阻力均呈减小趋势。

参考文献：

[1] Savchenko Y N， Semenenko V N， Putilin S I. Unsteady supercavitated motion of bodies[J]. International Journal of Fluid Mechanics Reaearch，2000，27(1):109-137.

[2] 施红辉，孙亚亚，周杨洁，等.水下航行体表面凹槽数对超空泡流场影响的数值模拟[J].弹道学报，2017，29(1):44-50．

[3] 颜开，褚学森，许晟，等.超空泡流体动力学研究进展[J].船舶力学，2006，10(4):148-155．

[4] 施红辉，周素云，张晓萍，等.水下超空泡流体机械的机理和技术研究综述[C]//第十五届全国激波与激波管学术会议论文集.杭州，2012:554-564．

[5] 周后村，向敏，张为华.超空泡航行体控制面与主空泡相互作用分析[J].兵工学报，2017，38(5):949-958．

[6] 胡晓.水下航行体空泡流数值模拟研究[D].哈尔滨:哈尔滨工程大学，2012:39-50．

[7] 施红辉，胡俊辉，周浩磊.完全超空泡出水的实验研究及理论分析[J].空气动力学学报，2014，32(4):544-550．

[8] 陈瑛.自然空泡流数值模拟方法研究[D].上海:上海交通大学，2009:12-30．

[9] 罗格维诺维奇Γ Β.自由边界流动的水动力学[M].施红辉，译.上海:上海交通大学出版社，2012:97-126．

[10] Savchenko Y N. Perspectives of hydrodynamic drag reductoin methods [J]. Internation Journal of Fluid Mechanics Research，2001，28(4):526-540．

[11] 栗夫园，张宇文，党建军，等.锥形空化器流体动力特性研究[J].兵工学报，2014，35(7):1040-1044．

[12] 徐清沐，薛雷平.带凹槽和带凸起的圆盘空化器超空泡流数值研究[J].水动力学研究与进展:A辑，2013，28(5):518-525．

[13] 金大桥，王聪，魏英杰，等.水下射弹自然超空泡减阻特性的数值模拟[J].工程力学，2010，27(6):202-208．

[14] 牟斌，江雄，王建涛.空化流动隐式求解方法研究[J].空气动力学学报，2017，35(1):27-32．

[15] 张明远，景思睿，李国君.高等工程流体力学[M].西安:西安交通大学出版社，2006:310-320．

[16] 魏海鹏，符松.不同多相流模型在航行体出水流场数值模拟中的应用[J].振动与冲击，2015，34(4):48-52．

[17] Yang D， Xiong Y L， Guo X F. Drag reduction of a rapid vehicle in supercavitating flow[J].International Journal of Naval Architecture & Ocean Engineering，2016，9(1):35-44．

[18] 施红辉，罗喜胜.可压缩性和高速多相流动[M].北京：中国科学技术大学出版社，2014:186-199．

[19] 叶少东.水下航行体空化器形状对超空泡流影响的数值模拟[D].杭州:浙江理工大学，2017:20-35．

[20] 易文俊，王中原，熊天红，等.水下高速射弹超空泡减阻特性研究[J].弹道学报，2008，20(4):1-4．