不同头型旋成体入水初期流场特性数值分析

吕红庆,许 磊,王振清

(1.烟台哈尔滨工程大学研究院， 山东 烟台 264000； 2.哈尔滨工程大学 航天与建筑工程学院， 哈尔滨 150001)

1 引言

入水过程是跨介质武器从空中飞行转入水下航行的一个重要环节，旋成体以一定初速度撞击自由液面，将自身一部分动能传递给液体，液体发生复杂流动并与航行体表面分离。随着入水深度的增加，旋成体两侧形成了稳定的空泡区域，空泡区域主要由液体相变产生的水蒸气和高速诱导进入的空气组成。此外，在入水过程中，旋成体头部产生具有强瞬时性、非定常性及高载荷性的冲击载荷，同时伴随着相变、湍动等非定常流体动力学现象以及弹体弹道难可控性等问题[1]。因此研究旋成体不同头型在入水过程的空泡形态演化与瞬间冲击载荷特性对跨介质兵器的发展具有指导意义。

早在十九世纪末期，学者们就已针对入水问题开展了试验和理论研究，Worthington[2]利用刚兴起的闪光摄影技术对小球入水进行了试验研究，通过改变液体介质，观察到入水过程中液体飞溅与空泡演化现象，同时分析了液体属性、表面张力等因素对液体飞溅现象的影响，为之后的入水试验奠定了基础。Watanabe[3]通过对尖锥体进行入水试验研究，测量出尖锥体在不同入水时刻所受到冲击载荷大小。在早期获得的试验数据基础上，最早进行入水问题理论研究的是T.von Karman[4]，他假定流体无黏、无旋和不可压缩，物体在穿越气-水界面时液面不发生变化，同时引入速度势理论，设定自由液面速度势为零，提出附加质量法，基于动量守恒定律推导出入水冲击载荷计算公式，该计算公式对于二维小底升角(α)楔形体入水冲击载荷预测较为准确。Wagner[5]提出渐进匹配方法，针对物体入水过程中液面抬升现象，将流体分为内外流域，内流域中进行局部射流分析，外流域提出近似平板理论，使物体入水过程的数学描述更加符合实际问题。Wagner的理论研究思想为后续入水冲击理论奠定了坚实的基础。

随着入水问题数学模型的完善，Logvinovich[6]在他的专著中提出了空泡独立膨胀原理，根据相似理论给出了计算空泡尺寸的半经验公式，成为空泡理论的奠基人。Lee等[7]在空泡理论基础之上，推导了普适性的数学模型，能够解释部分简单的入水现象。与此同时，随着高性能计算机和先进计算方法的出现，对入水问题的研究迎来了黄金阶段。Park等[8]提出了计算高速入水冲击载荷以及跳弹现象的数值方法，假定入水冲击在极短的时间内发生，入水区域的流动为无粘势流，采用源面板法进行求解，计算获得的数值数据与试验数据具有良好的一致性。Neaves等[9]采用有限体积法，引入自然空化模型和Tait状态方程模拟了射弹高速垂直入水空泡演化规律，同时考虑了流体介质的压缩性和空化潜热效应，获得了高速射弹流场结构，数值结果与实验结果一致。Abraham等[10]对小球在不同雷诺数、入水速度和液体表面张力情况下展开了数值模拟，认为动量传递的原因主要是拖曳力而非摩擦阻力。

蒋运华等[11]研究了弹体在不同通气量和不同倾斜角度入水情况下近水表面空泡特性，包括表面封闭、颈缩、深闭合以及壁面波动等现象，同时针对超音速弹体串并联入水过程中空泡几何形状和阻力特性展开了研究[12]。王聪等[13]对不同角度锥头弹体高速垂直入水进行了数值模拟，得到不同头型条件下高速入水运动参数及空泡形态发展规律、流场的压力分布及速度分布规律，分析了头型对入水空泡流场的影响。段文洋等[14]采用试验与数值模拟相结合的办法，研究了弹体入水初期的空泡流动与演化过程，分析了头部形状对该过程的影响。潘光等[15]提出一种非对称机头的AUV，采用数值模拟的方法分析了在不同的入水速度、机头的形状对AUV入水过程的弹道特性的影响。

目前大多数的入水问题研究主要集中于旋成体低速入水过程，对于旋成体由低速到高速入水的研究还较少。因此本文对跨介质旋成体低速到高速入水过程开展数值模拟研究，基于雷诺平均N-S方程、VOF多相流模型、realizablek-ε湍流模型，引入Schnerr and Sauer空化模型，建立能较好描述入水过程的数值计算模型，采用重叠网格技术进行仿真计算。对不同头型旋成体在不同工况下的入水过程展开研究，分析旋成体在入水过程中流场参数、空泡演化、弹体动力学特性与瞬间冲击载荷特性。所得研究成果可为跨介质兵器设计提供一定参考。

2 数值方法与计算模型量

2.1 物理模型

物理模型采用旋成体构型，以空投鱼雷、跨介质射弹等几何特征的跨介质武器作为目标模型。弹体长度设定为航行体直径的3倍且不考虑尾翼的影响，计算所采用的具体模型如图1所示。

图1 旋成体物理模型示意图Fig.1 Physical model of the vehicle

旋成体的直径为324 mm，弹身长度为972 mm，球头半径为162 mm，航行体材料为普通碳素钢，密度为7.85 g/cm3。具体几何尺寸如表1所示。由于针对旋成体入水过程流场参数变化与空泡演化特性开展研究，所以将模型设定为刚体结构不考虑旋成体变形与破坏。

表1 旋成体参数

2.2 控制方程

本文采用均质平衡多相流理论，选取Volume of Fluid(VOF)模型[16]，模拟两相或多相具有明显分界面的流体运动问题。将气相与液相作为混合的单一介质，通过各组分的所占混合物的体积分数比值来描述各组分的流动状态，不同相共用同一套基本控制方程组，根据质量守恒方程，连续性方程可表示为：

(1)

(2)

(3)

(4)

式中：P为压力；S为源项；τij为剪切应力，表达式如下：

(5)

湍流模型选取为realizablek-ε模型[17]，该模型引入了旋转和曲率相关的数值，采用涡旋粘性各向同性假设，在模拟流动分离、旋转流动、圆孔射流和射流撞击等问题中具有良好适应性[18]。特别是对于射流曲率较大的情况，realizablek-ε模型具有很好的表现。该模型湍动能k和湍动能的耗散率ε的输运方程表达为：

(6)

式(6)中，我们通常取C1ε=1.44、C2=1.9、σk=1、σε=1.2。

realizablek-ε模型是对standardk-ε模型的修正，也被称之为可实现的k-ε模型。在计算过程中，为了保证计算结果的可实现性(Realizability)，Cμ不再视为常数，来保证正应力数学约束关系的实现。

预测常数Cμ的公式如下：

(7)

本文为确保仿真更加符合物理实际，所以采用Schnerr and Sauer模型[19]对弹体入水过程中的空化现象进行求解。水蒸气相的质量输运方程如下式所示：

(8)

式中:RB为气穴半径，RB=1×10-6m。不可凝结气体体积分数αnuc=5×10-4、预测常数Fvap、Fcond的取值分别为50与0.001。

本文在计算过程中采用有限体积法对构建的控制方程进行离散化处理，利用PISO算法，对压力和速度耦合问题求解，该算法对SIMPLE算法进行了改进，是一种基于压力的隐式算子分裂方法，计算收敛性健壮且计算效率高。控制方程中对流项采用二阶离散格式，扩散项采用二阶中心差分格式，混合相体积率的离散采用CICSAM格式，时域采用隐式离散方法进行求解。

2.3 计算域与网格划分

网格划分采用重叠网格技术[20]，该技术能够确保在弹体入水过程中网格不发生变形，重叠区域网格尺寸与背景网格处尺寸基本保持一致，这样物体在运动过程中运动网格与背景网格插值区域的网格质量保持最优。本研究的计算域由弹体运动域和背景区域两部分构成，如图2所示。

图2 计算域示意图Fig.2 Schematic diagram of calculation domain

图2(a)所示，运动域为围绕弹体的钝头体圆柱，圆柱直径为弹体直径的2倍，计算域长度为1 200 mm,约为4倍弹体直径。图2(b)所示，为计算背景区域，该区域采用立方体结构，长为6 000 mm，宽为3 000 mm，高为6 000 mm，上半部分为空气域，高度为1 500 mm，下半部分为水域，水深设定为4 500 mm，弹体的初始位置设定在空气域中，距离水面高度为300 mm。由于本次计算的模型的为轴对称模型，所以采用一半的模型进行计算以节约计算的资源和时间。

采用结构网格划分运动区域与背景域，如图3所示。计算域y-z平面的网格的剖分视图，z轴的负方向为重力加速度的方向，坐标的原点为旋成体模型质心位置。在旋成体周围进行局部的加密处理，以更好捕捉入水过程中空泡的演化与瞬间的冲击载荷。临近壁面处，法向速度会存在较大的梯度，对于该区域的计算采用壁面函数法。本次计算过程中雷诺数较高且存在分离与射流流动的现象，我们使用非平衡壁面函数法，采用结构化网格进行划分，运动域网格数量为60万，背景域网格数量为240万。背景计算域的四周边界条件与底面设定为壁面(wall)，它的上表面设定为压力出口(pressure-outlet)，在旋成体周围的区域边界设定为重叠网格可以识别的over-set边界条件，下落过程中，采用六自由度运动方式，背景域网格与运动域网格进行插值计算传递数据。

图3 y-z平面剖分网格示意图Fig.3 y-z plane mesh generation

3 计算结果与分析

3.1 计算方法验证

有必要验证本文中所采用的数值计算方法的可靠性，以确保研究结果的可信性。南京理工大学瞬态物理国家重点实验室[21]关于不同头型回转体入水的试验数据为参考，将数值模拟得到的空泡演化和回转体入水速度变化与试验结果进行对比。

验证计算的试验装置模型如图4所示，试验中回转体为长度40 mm、直径8 mm的实心圆柱体，材料为普通碳素钢，密度为7.85 g/cm3。液体介质为室温下的水，密度为1 g/cm3，水温为25 ℃，通过高速摄像机记录回转体在入水过程中对自由液面的扰动。

图4 试验装置模型示意图Fig.4 Schematic diagram of the test equipment

试验在高度为500 mm、底部尺寸为250 mm×250 mm且设有防护层的水槽中进行。在水槽的侧面贴有坐标纸，记录回转体运动的位置变化。回转体的初始位置在水槽的上方，被电磁铁吸附固定，释放后自由落体垂直入水，接触水面的入水速度为1.67 m/s。通过图5分析，可以清楚地看到，入水仿真所获得的回转体被空泡包裹、空泡的形成与演化、空泡的颈缩与闭合以及产生液体的射流等现象与试验结果基本一致，验证了本文数值计算方法的适用性。

图5 数值仿真与试验现象过程示意图Fig.5 Comparison between numerical simulation and experimental phenomena

3.2 无关性验证

在数值求解的过程中，离散点分布在网格的节点上，对于物理量过渡平缓的位置，网格节点布置可以较为稀疏，对于梯度变化较大区域，需要增加网格节点数量。选择适当的网格疏密程度不仅可以保证仿真结果的精确度，还能够节约数值计算资源与时间。使用2.1～2.3节中给定的计算模型与边界条件，以半球头旋成体100 m/s的初速进行垂直入水计算，选取出三套不同密度的网格模型，网格数量见表2。通过验证计算，最终确定出适合本次计算的网格尺度。

表2 网格划分数量

计算过程中设置监视点，得到旋成体入水过程中头部最大压力和速度随时间变化的曲线，以一个标准大气压P0=101 325 Pa作为冲击压力峰值的基准，对头部压力峰值进行无量纲化的处理。如图6与图7所示。可以看到，在计算过程中三套网格压力与速度随时间的变化趋势基本一致，在入水初期，三套网格对于压力峰值与后期压力变化的计算曲线基本一致。观察旋成体入水速度变化，可以看到稀疏网格的速度曲线虽然趋势与另外2种网格密度保持一致，但在不同时刻时与其他2种情况的航行体速度存在较大差异。通过对计算的数据差异以及所需的计算资源进行比较，最终确定在后续模拟中选用运动域60万网格、背景域240万网格作为本次计算的网格数量。

图6 旋成体无量纲压力峰值随时间变化曲线Fig.6 Dimensionless pressure peaks of axisymmetric body varying with time

图7 旋成体速度随时间变化曲线Fig.7 Variation of speed of axisymmetric body with time

3.3 旋成体头部外形影响

影响旋成体的入水流场特性的因素有很多，本小节主要研究相同长径比、相同质量密度情况下，在100 m/s垂直入水速度下，头型对入水过程的影响。图8给出了不同头型在相同入水深度时空泡形态，横坐标为空泡半径(R)与旋成体半径(R0)的比值，纵坐标为空泡长度(H)与旋成体直径(D)的比值，我们观察可知3种头型旋成体空泡均在旋成体的肩部产生，空泡的整体发展轮廓基本一致。

空泡沿着径向逐渐向外扩张，空泡的半径随之增大，自由液面处的空泡半径最大。不同头型旋成体入水初期，在相同深度的情况下，平头体产生的空泡半径最大，半球头体次之，尖锥头体的空泡尺寸最小。对于不同头型入水空泡大小差异这一现象，我们可以根据文献[7]，假设流体介质为无黏、无旋、不可压缩的理想流体，根据能量守恒定律，弹体在入水过程中，自身动能的损失都转化为液体压力势能和动能，由此可以得到入水空泡半径预测公式，经过一系列推导，并略去高阶小量，可知空泡半径的大小主要与弹体头型在空化数为零时的阻力系数及入水时间有关，阻力系数越大那么在自由液面处形成的空泡尺寸就越大，在未发生空泡闭合时，旋成体入水时间越长，自由液面处空泡半径越大。

图8 不同入水深度下空泡形态曲线Fig.8 Comparison of cavitation shape curves at different water entry depths

入水初期，旋成体接触水面后头部产生巨大压力差，尖锥头斜面可以更好地将压差转化为水的动能产生液体的飞溅，所以尖锥头体对自由液面扰动最大，溅起的水花区域最大，平头体对水面扰动最小，所溅起的水花区域也最小，半球头体介于两者之间。

旋成体头部外形影响着入水初期瞬间的压力载荷峰值的大小，如图9(a)所示为3种头型在入水初速100 m/s时的压力峰值曲线。由图可知压力峰值出现在入水的初期，即当旋成体头部接触到水面时，冲击载荷达到峰值，但压力峰值持续时间较短，随着入水深度增加，旋成体被空泡包裹，头部压力逐渐趋于稳定。在相同的入水速度下，可以看出不同的旋成体头型的局部载荷具有不同的特点，平头体的瞬间冲击压力峰值要明显大于另外2种头型的旋成体，半球头体次之，尖锥头体的冲击压力峰值最小。平头体的峰值宽度最窄，半球头次之，尖锥头峰值宽度最大，由此可以得出旋成体头部压力峰值与峰值宽度成反比例的关系。入水一段时间后，3种头型旋成体头部的峰值逐渐趋于稳定，入水冲击载荷的波动情况表现基本一致。

如图9(b)所示为不同头型旋成体撞击水面后速度变化曲线，当未接触到水面时，3种头型受重力加速度作用，三条速度曲线的变化趋势基本一致。当旋成体接触到水面时，可以很明显的看到平头体在撞水的瞬间速度曲线出现明显的拐点，受到的阻力和加速度最大，在接触水面的初期，尖锥头体的速度要比半球头体的速度衰减的慢，但随着入水深度的增加，尖锥头是速度衰减率要大于半球头体。在整个计算时间域内半球头旋成体速度曲线过渡较为平滑，速度逐渐下降，近似于线性速度下降趋势。随着时间推移，旋成体逐渐被空泡所包裹，3种头型旋成体的速度下降速率趋于相似。

图9 不同头型入水参数曲线Fig.9 Comparison of water entry parameters of different head types

3.4 入水速度影响

本小节探究不同头型旋成体在不同入水速度下空泡演化及冲击载荷特性，入水初速度由低速向高速过渡，具体计算工况参数见表3所示。

表3 入水速度工况参数

如图10所示为不同入水初速度下同一时刻的无量纲入水空泡演化过程。由图可知，相同入水速度下，在入水时间相同时，不同头型的入水深度存在明显差别。随着入水初速度的提高，自由液面溅起的水花速度方向逐渐偏离弹体轴线向外扩展。入水初速度越大，液体射流速度越快，自由液面受到的扰动越大，溅起的水花面积越大。因为航行体入水初速度越大，自身的动能越大，撞击水面后液体获得的动能与压力势能也越大，所以液体飞溅范围扩大。同一入水时刻下半球头旋成体的入水深度要大于尖锥头与平头旋成体的入水深度。

图10 不同初速度下相同时刻空泡演化过程曲线Fig.10 Comparison of cavitation evolution at different initial velocities at the same time

对比旋成体在不同初速度下入水过程的衰减速率，由于初始速度相差较大，不易对比观察，故采用归一化方法。通过对比各个时刻旋成体速度与入水初速度的比值，绘制出不同初速度下旋成体相对速度随时间的衰减变化率，如图11所示半球头旋成体在不同入水初速度下相对速度变化率曲线。可以看出，在旋成体接触水面后，不同初始速度下旋成体相对速度变化率出现明显的拐点，入水初速度越大，相对速度的衰减率越大，呈现正比例的关系。航行体接触水面的瞬间，所受到的反向加速度急剧增加，这个过程的持续时间很短，且入水初速度越快，这个过程持续的时间就越短，旋成体更快被空泡所包裹。随着入水深度的增加，在浮力、惯性力与重力的共同作用下，旋成体相对衰减速率逐渐趋于稳定，速度按一定比例的下降。半球头航行体在初速度50 m/s和100 m/s时的相对速度衰减曲线间隔要大于 100 m/s和150 m/s时曲线间隔，并且变化明显。

图11 半球头旋成体相对速度变化率曲线Fig.11 Relative velocity change rate of a blunt head body

对于同一头型的旋成体而言，入水初速度决定了瞬间冲击压力峰值的大小。从图12可以看出，半球头旋成体瞬间压力峰值随着入水初速度的增加而变大；对比初速度为50 m/s、100 m/s、150 m/s时不同头型旋成体压力峰值，可知相同头型在150 m/s速度下的压力峰值约为50 m/s速度下的8～10倍，100 m/s速度下的压力峰值约为50 m/s速度下的3～5倍，可知旋成体在入水所受到的压力峰值随着初始入水速度的增加呈现近似平方关系的正相关增长。将压力大小为峰值80%以内载荷分布的宽度定义为冲击压力峰值的幅值宽度，可以看到，随着入水速度的增大，压力峰值的幅值宽度也会变大。

图12 半球头旋成体无量纲压力变化曲线Fig.12 Curve of dimensionless pressure of a blunt head body

图13为不同入水初速度下半球头旋成体的阻力系数变化曲线，由图可得当旋成体接触自由液面时阻力系数达到峰值且3种速度下峰值基本相近。随着旋成体入水深度的增加，气相与液相的流场基本趋于稳定，同时不同初速度时的阻力系数趋于稳定。可以看出入水初速度对旋成体的阻力系数影响不大。

为了更好的分析不同头型旋成体在不同初速度下垂直入水时的冲击压力大小，提取不同头型旋成体不同初速度下的冲击压力峰值及其在旋成体进入空泡后的稳定值，提取所得具体数值如表4所示。

图13 不同入水速度阻力系数变化曲线Fig.13 Variation curve of resistance coefficient at different water entry velocities

表4 不同入水速度下冲击压力峰值与稳定值Table 4 Peak value and stable value of pressure at different water entry velocities

对于不同头型旋成体而言，初速度为50 m/s时，旋成体所受到的冲击压力峰值及其入水之后的稳定值最小。不同头型在不同初速度下头部冲击压力峰值的无量纲变化情况如图14所示。由图可知，随着入水初速度的增大，平头体的冲击压力峰值与稳定值始终大于半球头体和尖锥头体。随着旋成体入水速度增大，尖锥头体与半球头体头部冲击压力峰值的变化幅度相较平头体变化不明显。

图14 不同头型旋成体无量纲压力峰值随速度的变化曲线Fig.14 Variation of dimensionless pressure peak value with velocities for different head-shaped bodies

4 结论

采用数值模拟方法，研究了3种不同头型旋成体初速度从低速到高速入水过程中空泡形态与冲击载荷特性，得到以下结论：

1) 平头旋成体入水过程形成的空泡半径最大，半球头次之，尖锥头空泡半径最小；旋成体入水接触自由液面瞬间，头部产生巨大压力差，尖锥头斜面可以更好将压差转化为水的动能将水排开，所以尖锥头旋成体在入水过程中溅起的水花面积最大，半球头旋成体次之，平头旋成体溅起的水花面积最小。

2) 同一入水时刻下，半球头旋成体的入水深度最深，尖锥头次之，平头体入水深度最浅；不同头型旋成体从低速到高速的入水过程中，将旋成体不同初速度下入水速度变化进行归一化处理，3种头型旋成体的速度衰减特性表现基本一致，旋成体初速度越大，入水过程中速度衰减率就越大。

3) 不同头型旋成体所受到的冲击压力峰值随着入水初速度增加呈现近似平方关系的正相关增长，平头旋成体冲击压力峰值变化幅度受入水初速的影响较另外2种头型旋成体变化更加剧烈。随着入水初速度的升高，冲击压力峰值的幅值宽度也会变宽，且入水初速度由低速到高速的变化对旋成体的阻力系数影响不大。