旋成体高速入水可压缩性影响研究

李国良， 尤天庆， 孔德才， 李静， 周伟江

(1.中国航天空气动力技术研究院 一所， 北京 100074；2.北京宇航系统工程研究所， 北京 100076)

0 引言

航行体入水运动是一个具有强非定常特性的复杂多相流动过程，这一过程具有强瞬时及高载荷特性。航行体在入水过程中会形成一个由空气和空化产生的蒸汽混合空腔，形成过程涵盖多相之间的掺混、液体可压缩、自由面与空泡间相互作用等复杂物理问题。现阶段航行体高速入水问题研究更具有工程意义，它涉及武器装备如空投鱼雷、超空泡射弹、多次出入水导弹等的关键技术问题。

航行体入水问题的实验研究国外开展得较早[1-5]，国内时素果等[6]采用高速摄影方法对细长体以速度60 m/s入水过程开展实验研究，观察入水过程空泡形态演化规律，采用压力传感器获得了物体表面压力变化规律。蒋运华等[7]采用实验方法，在实验水箱中开展约束和不约束运动体姿态的入水实验，分析了有无扰动下空泡尺寸与弗劳德数和欧拉数之间的关系。何春涛等[8]通过垂直和倾斜两种入水方式开展实验研究，分析了入水过程中的一系列流动现象，以及空泡的生成、发展和闭合，并且开展了串联和并联入水时空泡的生成和相互之间的影响。

由于实验研究成本较高及周期偏长的原因，入水问题研究的早期更多采用数值仿真手段。夏维学等[9]采用流体体积函数(VOF)多相流模型、k-ω湍流模型(k为湍流动能，ω为比耗散率)及动网格技术，数值模拟了旋转角速度为199 rad/s、入水速度为5.45 m/s条件下入水空泡的形成、发展、闭合和溃灭的演化过程。孙钊等[10]采用VOF多相流模型和连续表面力(CSF)模型，对亲水性及疏水性球体垂直入水过程开展了数值模拟研究，分析了亲水性及疏水性球体低速入水对空化初生以及流体动力的影响。马庆鹏等[11]基于有限体积法和VOF方法，对小型锥头圆柱体高速垂直入水开展数值模拟，研究了入水空泡演化及深闭合过程空泡流场的流动特性和压力分布等问题。Marrone等[12]采用光滑粒子流体动力学(SPH)模型对航空器和直升机水上迫降涉及的机理问题开展研究，把飞行器简化为平板，对其在低速小攻角条件下水击载荷进行数值模拟，并与实验数据与势流理论进行了对比。Sean等[13]采用浸没边界法建立了多相流数值求解方法，克服了水、空气、蒸汽多相之间相界面大密度比的难题，并且模拟了二维方柱和圆柱低速入水问题。Neaves等[14]忽略湍流对空泡发展的影响采用预处理方法和低耗散迎风格式，对二维旋成体入水问题开展了数值模拟，验证了该方法的适用性。陈晨等[15]采用数值模拟方法研究了小型运动体高速倾斜入水时的空泡流动特性，并深入研究了不同入水角度对空泡形态、力学特性的影响。

目前航行体入水问题研究通常不考虑液体的可压缩性。本文基于压力耦合方程组的半隐式(SIMPLE)算法，在考虑液体可压缩性条件下，对锥柱旋成体以50 m/s、100 m/s、200 m/s、400 m/s、800 m/s速度入水的过程进行数值仿真，研究湍流模型对入水过程的影响、液体可压缩性对入水冲击载荷及空泡形态的影响，以及入水速度对入水后运动特性的影响。

1 计算方法

1.1 空化多相流方程

混合相连续方程：

(1)

混合相动量方程：

(2)

蒸汽相输运方程：

(3)

式中：ρm为混合相密度，ρm=αlρl+αvρv+αgρg，αl、αv、αg分别为液相、蒸汽相、气相的体积分数，ρl、ρv、ρg分别为液相、蒸汽相、气相的密度；t为时间；u为速度，x为坐标，下标i和j分别为1，2，3，表示3个坐标方向；p为压力；μm为层流混合相黏性系数，μm=αlμl+αvμv+αgμg，μl、μv、μg分别为液相、蒸汽相、气相的黏性系数；Re、Rc为空化源项；μt为混合相湍流黏性系数；g为重力加速度。

1.2 空化模型

采用Zwart-Gerber-Belamri空化模型：

(4)

式中：Fv为蒸发系数，Fv=50；rn为成核区的体积分数，rn=5×10-4；Rb为气泡半径，Rb=1 μm；pv为饱和蒸汽压；Fc为凝结系数，Fc=0.01.

1.3 湍流模型

湍流模型对于空泡形态的演化及运动过程有较大影响，本文采用剪应力传递(SST)k-ω、标准(Standard)k-ε(ε为湍流动能耗散率)、重整化群(RNG)k-ε、可实现(Realizable)k-ε4种湍流模型开展对比分析研究。

1.4 可压缩状态方程

对于高速入水，入水载荷及空泡界面主要受液相影响，因此本文中气相和蒸汽相仍按照不可压缩来处理。液相要考虑其可压缩性，在实现液相可压缩性时应用Tait状态方程。不含温度修正的Tait状态方程为

(5)

式中：pr和ρr为参考压力和参考密度；K和n分别为温度和压力的弱函数，这里设为常数，K=3×108，n=7.

利用Tait方程定义流体密度变化的同时，还要定义流体中的声速c方程：

(6)

1.5 数值求解方法

流场采用基于压力的SIMPLE算法进行求解。动量方程中的对流项采用2阶迎风格式，压力项离散采用PRESTO差分格式。多相流模型采用VOF+Level Set方法，时间项采用显示格式。本文主要模拟高速入水问题，在入水过程中会产生空泡，需要多相流与动网格耦合求解。动网格计算采用动态层法，为保证射弹在运动过程中空泡产生的区域网格状态始终不变，采用全域网格随射弹运动的动网格方法，网格更新分别在流场入口及出口边界处，在流场入口与出口新生网格高度始终与原网格高度相同。计算模型如图1所示，模型采用锥柱旋成体，具体参数如表1所示，其中L为模型长度，D为模型直径，θ为前部圆锥体锥角，m为模型质量。

图1 计算模型Fig.1 Computational model

表1 模型参数
Tab.1 Model parameters

L/mmD/mmθ/(°)m/kg50101270.01

图2为计算域示意图，由于计算模型为旋成体外形，采用二维轴对称网格。图2中，上部为压力进口边界条件，设置为标准大气压力，底部为压力出口边界条件，按照公式p=p∞+ρlgh设置，本文中p∞为1个标准大气压，h表示水深，空气域与液相域的内侧为轴对称边界。

图2 计算域示意图Fig.2 Computational domain

2 结果分析

2.1 湍流模型比较分析

湍流采用SSTk-ω、Standardk-ε、RNGk-ε、Realizablek-ε4种模型，比较湍流模型对入水速度衰减及入水深度的影响。根据文献[16]提出的入水速度和空泡半径随时间变化的理论公式，得到瞬时加速度和空泡半径分别为

(7)

(8)

式中：vp为瞬时速度；A0为模型横截面；Cd为阻力系数，Cd=Cd0+σ,σ为空化数，Cd0为σ=0时的阻力系数，本文中Cd0=0.498；r(x,t)为空泡半径；tb为前一运动时刻；r0为tb时刻空泡半径；A(x)、B(x)的表达式见参考文献[16]。

图3为旋成体入水时间为0.006 s时的速度衰减曲线，可以看出4种湍流模型的数值总体趋势上与理论解一致，但是在量值上有差异。在0.006 s时刻，理论解瞬时速度为39.08 m/s，Standardk-ε湍流模型模拟瞬时速度为33.57 m/s，SSTk-ω湍流模型模拟瞬时速度为38.45 m/s，RNGk-ε湍流模型模拟瞬时速度为37.24 m/s，Realizablek-ε湍流模型模拟瞬时速度为37.92 m/s. 从数值比较来看，使用SSTk-ω湍流模型得到的速度衰减，无论从整体趋势上还是从瞬时点上都与理论解误差最小。

图3 旋成体入水速度随时间衰减曲线(入水速度100 m/s)Fig.3 Velocity during water-entry vs. time(water-entry velocity of 100 m/s)

图4 无量纲入水深度随时间变化(入水速度100 m/s)Fig.4 Nondimensional water-entry depth vs. time (water-entry velocity of 100 m/s)

图5 旋成体入水1 ms时数值模拟空泡形态及与理论解比较(Standard k-ε模型，入水速度100 m/s)Fig.5 Numerically simulated cavitation configuration and its comparison with theoretical solution at 1ms of water-entry(Standard k-ε model，water-entry velocity of 100 m/s)

图5～图8分别为采用4种湍流模型计算入水时间为1 ms时刻空泡界面及与理论解的比较图。由图5～图8可见，采用Standardk-ε、Realizablek-ε、RNGk-ε3种湍流模型数值模拟的空泡界面在水气界面处为向内收缩，采用SSTk-ω湍流模型数值模拟的空泡界面在自由面处为向外扩张，造成这一现象的原因是因为SSTk-ω湍流模型能够更好地模拟空泡界面剪应力的传递作用，因此采用SSTk-ω湍流模型得到的空泡界面与理论值符合最好。为了验证本文方法的准确性，进一步在高速入水理论解及实验结果比对方面开展工作。

图6 旋成体入水1 ms时数值模拟空泡形态及与理论解比较(Realizable k-ε模型，入水速度100 m/s)Fig.6 Numerically simulated cavitation configuration and its comparison with theoretical solution at 1 ms of water-entry(Realizable k-ε model，water-entry velocity of 100 m/s)

图7 旋成体入水1 ms时数值模拟空泡形态及与理论解比较(RNG k-ε模型，入水速度100 m/s)Fig.7 Numerically simulated cavitation configuration and its comparison with theoretical solution at 1 ms of water-entry (RNG k-ε model，water-entry velocity at 100 m/s)

图8 旋成体入水1 ms时数值模拟空泡形态及与理论解比较(SST k-ω模型，入水速度100 m/s)Fig.8 Numerically simulated cavitation configuration and its comparison with theoretical solution at 1ms of water-entry (SST k-ω model，water-entry velocity of 100 m/s)

采用SSTk-ω湍流模型对入水速度200 m/s、400 m/s 2种情况开展数值模拟，图9为入水衰减速度数值模拟与理论值的对比，图10为无量纲入水深度数值模拟与理论值的对比。从图9和图10中可以看出，数值模拟与理论值具有很好的一致性，整体误差在2%以内。

图9 旋成体入水速度随时间衰减曲线Fig.9 Velocity during water-entry vs. time

参考文献[15，17]的实验模型、实验状态、网格总数及设置，采用本文建立的数值模拟方法与该高速入水实验进行对比研究，图11所示为入水速度衰减数值模拟与实验值的比较。由图11可见，初始入水速度为141.15 m/s条件下，数值模拟与实验结果一致性较好，整体误差控制在1%以内。

图12所示为数值模拟空泡形态发展与文献[15，17]实验结果对比，可见数值模拟基本能够得到实验中空泡发展的过程，在空泡形态上数值模拟与实验结果基本一致。因此，从数值模拟与理论结果的比较以及数值模拟与参考实验的比较来看，本文建立的数值方法是正确的。

2.2 流体可压缩性的影响分析

高速入水会产生较大的冲击载荷，液体可压缩性对冲击载荷有较大影响，在入水速度分别为50 m/s、 100 m/s、200 m/s、400 m/s、800 m/s 5种条件下，分析考虑液体可压缩和不可压缩两种处理方式对入水力学特性的影响。表2为数值模拟结果。

根据表2的结果可知：在入水速度为50 m/s、100 m/s时，两种处理方式最大载荷出现的时间相同，冲击载荷基本相当；在入水速度为200 m/s时，不可压缩流体入水冲击载荷最大值出现在入水1.6×10-5s时，峰值为3 765 N，可压缩流体入水冲击载荷最大值出现在入水1.7×10-5s时，峰值为3 474 N，二者相差7.7%，最大密度提高2%；在入水速度为400 m/s时，不可压缩流体入水冲击载荷最大值出现在入水0.95×10-5s时，峰值为13 323 N，可压缩流体入水冲击载荷最大值出现在入水1.1×10-5s时，峰值为12 211 N，二者相差9.1%，最大密度提高6.5%；在入水速度为800 m/s时，不可压缩流体入水冲击载荷最大值出现在入水0.55×10-5s时，峰值为55 478 N，可压缩流体入水冲击载荷最大值出现在入水0.875×10-5s时，峰值为37 718 N，二者相差47.1%，最大密度提高32.0%.

图13为入水速度为100 m/s时两种处理方式冲击载荷的变化历程，图中Fm为冲击载荷；图14为冲击载荷最大时密度等值线图；图15为入水速度为800 m/s时两种处理方式冲击载荷的变化历程；图16为冲击载荷最大时密度等值线图。由图13～图16可知，入水速度在100 m/s时液体可压缩性对入水冲击载荷基本无影响，入水速度≥200 m/s时，液体可压缩性对冲击载荷有影响，冲击载荷会变小，而且最大冲击载荷出现的时间会滞后，随着入水速度不断增加，滞后时间会加大。

图17、图18分别为入水速度800 m/s两种液体处理方式入水0.5 ms时刻的空泡形态图，其中图17为按照不可压缩方式数值模拟的空泡形态图，图18为按照可压缩方式数值模拟的空泡形态图。为了比较液体可压缩性对空泡形态的影响，对空泡轮廓线提取出来以作比较，如图19所示，可见考虑液体可压缩特性时空泡会有向内收缩现象。从图20可以看出，考虑液体可压缩性时，在空泡界面附近水的密度降低，在压差作用下周围液体向空泡内部挤压，造成空泡收缩现象。

图10 旋成体无量纲入水深度随时间变化曲线Fig.10 Nondimensional depth during water-entry vs. time

图11 旋成体入水速度数值模拟与实验结果[15，17]对比(入水速度141.15 m/s)Fig.11 Numerically simulated results vs. experimental results[15,17] of velocity(water-entry velocity of 141.15 m/s)

图12 空泡形态数值模拟与实验结果[15，17]对比(入水速度141.15 m/s)Fig.12 Numerically simulated results vs. experimental results in Refs.[15] and [17] of cavitation configuration(water-entry velocity of 141.15 m/s)

图13 冲击载荷随时间变化(入水速度100 m/s)Fig.13 Impact load vs. time(water-entry velocity of 100 m/s)

图14 最大载荷时刻密度云图(入水速度100 m/s)Fig.14 Density contour under maximun impact load(water-entry velocity of 100 m/s)

图15 冲击载荷随时间变化(入水速度800 m/s)Fig.15 Impact load vs. time(water-entry velocity of 800 m/s)

图16 最大载荷时刻密度云图(入水速度800 m/s)Fig.16 Density contour under maximun impact load(water-entry velocity of 800 m/s)

表2 旋成体不同入水速度冲击载荷数值模拟
Tab.2 Numerically simulated results of water-entry impact load at different water-entry velocities

入水速度/(m·s-1)最大载荷发生时间/10-5s最大载荷/N最大密度/(kg·m-3)不可压缩可压缩不可压缩可压缩不可压缩可压缩504.44.431531410001001100 2.42.4853851100010042001.61.737653474100010204000.951.101332312211100010658000.5500.875554783771810001320

图17 旋成体入水0.5 ms时刻空泡形态图(入水速度800 m/s，不可压缩液体)Fig.17 Cavitation configuration of revolutionary body at 0.5 ms (water-entry velocity of 800 m/s, incompressible water)

图18 旋成体入水0.5 ms时刻空泡形态图(入水速度800 m/s，可压缩液体)Fig.18 Cavitation configuration of revolutionary body at 0.5 ms(water-entry velocity of 800 m/s, compressible water)

图19 空泡轮廓图的比较(入水速度800 m/s)Fig.19 Comparison of cavitation outlines (water-entry velocity of 800 m/s)

图20 旋成体入水0.5 ms时刻液体密度云图(入水速度800 m/s，可压缩流体)Fig.20 Density contour of water at 0.5 ms during water entry of revolutionary body (water-entry velocity of 800 m/s, compressible water)

2.3 入水速度对运动过程的影响

图21 旋成体无量纲入水速度随时间变化Fig.21 Nondimensional velocity during water-entry vs. time

图22 旋成体入水加速度随时间变化Fig.22 Acceleration during water-entry vs. time

3 结论

本文基于SIMPLE算法，在考虑液体可压缩性的条件下，对锥柱旋成体在50 m/s、100 m/s、200 m/s、400 m/s、800 m/s入水速度下开展数值模拟研究。得出以下主要结论：

1)建立了高速入水问题的数值求解方法，采用SSTk-ω湍流模型进行数值模拟的结果在入水运动过程及空泡形态方面与理论解一致性最好，同时与参考实验的结果开展了比较研究，验证了本文方法的准确性。

2)在入水速度≤100 m/s时，液体可压缩性对冲击载荷基本没有影响，在入水速度≥200 m/s时，需要考虑液体可压缩性；随着入水速度增加，液体可压缩性对入水冲击载荷影响越大，其会弱化入水冲击载荷及延缓最大载荷出现的时间。

3)入水速度越大，入水过程速度衰减越快，加速度值在入水初期较大；在计算模型周围被超空泡包裹的航行阶段，随着入水深度增加，加速度逐渐减小，而且加速度的变化逐渐平缓。