跨介质航行器能够顺利跨越空水界面，同时具备导弹武器高速远程和鱼雷武器隐蔽性好的优势，是对抗大型水面舰船的杀手锏。为了追求空中弹道良好的气动性能，跨介质航行器一般不采用传统水中兵器的回转体外形，更倾向于采用适合高超声速飞行的外形。这种非回转体的航行器在入水及水下运动过程中表现出的流体动力特性和运动特性较回转体外形有很大不同。随着跨介质航行器相关研究的不断深入，具有非回转体外形的航行器在跨越自由液面时的运动及受力问题逐渐受到关注。

跨介质航行器的研究始于20世纪30年代，苏联最早提出了跨介质航行器的概念，但相关研究由于技术受限而搁浅^[1]。21世纪后，随着大型舰艇防御体系的日趋完备，导弹从空中突防的难度增大，而水下超空泡武器作战范围有限，于是跨介质航行器的突防优势逐渐凸显，引发各国进一步开展研究^[2]。以麻省理工学院的Lincoln实验室为代表的研究机构^[3-5]围绕仿生跨介质飞行器开展了大量研究，实现了空/水介质跨越。近年来高超声速飞行器研究表明，高升阻比气动布局形式主要包括翼身组合体、翼身融合体、升力体以及乘波体等非回转体外形^[6]，而跨介质航行器由于具有较长的空中弹道，故倾向于采用上述面对称外形以获得更好的气动性能。

跨介质航行器的介质跨越问题始终是制约其发展的瓶颈^[7]，其入水过程涉及复杂的相变过程，流场、流体动力和弹道表现出强烈的瞬态特性，国内外对入水问题已展开了大量研究并取得一系列成果。Worthington等^[8]用闪光照相机对小球垂直入水过程进行了试验研究，获得了入水空泡图片。Truscott等^[9-10]研究了球体和射弹入水时的流体动力特性。侯宇等^[11]通过试验研究了射弹高速入水过程中空泡形态的变化以及空泡对弹道的影响。袁绪龙等^[12-13]开展了模型实验，利用高速摄像机拍摄了回转体入水过程的空泡形态演变，分析了头型、入水角度和入水速度对水下弹道的影响。郭子涛等^[14]开展了水平回转体入水实验，研究了回转体入水过程的弹道稳定性和空泡拓展特性。黄鸿鑫等^[15]研究了回转体射弹的头部形状和质心位置对高速入水稳定性的影响。马庆鹏等^[16]研究了带有不同角度锥头圆柱体高速入水时的流场特性和运动特性。胡青青^[17]讨了4种头型、3种长径比射弹高速入水时的稳定性和空泡形态，并分析了其对入水过程的影响。唐楚淳等^[18]对非轴对称头型弹体的入水弹道特性开展了弹丸低速垂直入水实验研究。华扬等^[19]通过高速摄影法对不同构型的非对称头型航行器低速入水时的入水弹道和空泡形态开展了实验研究。

关于入水问题的研究主要集中在回转体方面，关于非回转体的研究则是主要以空投鱼雷、水下飞机等为应用背景开展的，入水速度一般较低。随着高速跨介质航行器研究逐渐深入，非回转体外形在高速斜入水过程中的载荷和运动特性问题亟需解决。然而，从公开发表的文献来看，围绕非回转体外形开展的高速入水研究鲜见。为解决非回转体高速入水问题，本文采用数值仿真方法，建立了非回转体航行器高速斜入水模型，研究了非回转体外形对航行器入水过程中空泡形态、流体动力和弹道特性的影响，研究结果可为跨介质航行器的设计提供依据。

1 非回转体航行器高速入水数值模型 1.1 研究对象

本文以某非回转体跨介质航行器为研究对象，为了适应空中飞行工况，采用具有椭圆截面的外形。几何结构如图 1所示，航行器由空化器段、收缩段和柱段组成。

航行器总长为L, 头部空化器直径为D_n, 空化器锥角为60°。通过等比缩小回转体航行器收缩段和圆柱段截面在航行器立轴方向上的尺寸, 可以得到一系列横截面为椭圆形的非回转体航行器。若截面椭圆的长轴长度为D, 短轴长度为d, 则短轴与长轴之比为d/D, 记为k。对同一个航行器, 其收缩段和柱段上各截面的k是相同的。参考某跨介质航行器的入水工况, 在入水角为20°、初速度为200 m/s的典型入水工况下, 对不同截面形状(k值分别取0.8, 0.7, 0.6)的非回转体航行器入水过程进行数值仿真, 研究其空泡演化规律、流体动力特性和弹道特性, 选择回转体模型(k=1.0)作为对照组。

1.2 数值模型

非回转体跨介质航行器入水过程是一个典型的存在清晰界面的非定常多相流问题, 其中还涉及自然空化问题。本文基于商用流场仿真软件平台STAR-CCM+, 采用VOF多相流模型、Realizable k-ε湍流模型、Schnerr-Sauer空化模型和重叠网格技术进行数值仿真。

1) 多相流模型

航行器高速入水过程涉及空气、水、水蒸气三相, 属于非定常多相流问题。基于均质多相流理论, 本文采用VOF模型建立数值模型。

VOF模型要求在每个控制体中, 所有相的体积分数之和为1, 即各相体积分数满足关系式

(1)

式中, α_l, α_a, α_v表示液相、气相及水蒸气相的体积分数。

混合介质的连续性方程为

(2)

式中, i=1, 2, 3, 下同。

动量守恒方程为

(3)

(2)~(3)式中: u_i为速度分量; ρ_m和μ_m表达式分别为

(4)

(5)

对于VOF多相流模型, 空气相体积分数的输运方程为

(6)

2) 湍流模型

航行器高速入水过程中流场变化剧烈、雷诺数大、水气大量掺混, 为了精准模拟流动参数分布和航行器的流体动力特性, 需要对湍流进行建模求解。Realizable k-ε湍流模型是一种雷诺时均的湍流模型, 具有较高的精度和数值稳定性, 在涉及多相耦合、非定常问题的流场模拟中表现稳健, 可用于解决本文所研究的问题。该模型输运方程为

(7)

(8)

式中：; G_k为速度梯度引起的湍流动能；G_b为浮力引起的湍动能; Y_M为可压缩湍流脉动膨胀对耗散率的影响, 对于不可压缩流动该项为0; 为湍流黏性系数；Realizable k-ε湍流模型中C_μ是平均应变率与旋度的函数, 在平衡边界层惯性底层C_μ=0.09;模型常数C₂=1.9；C_1ε=1.44；σ_k=1.0。

3) 空化模型

在高速入水过程伴随着自然空化现象, 空化对流场特性有着显著的影响。Schnerr-Sauer空化模型具有形式简单、计算效率高、数值稳定性强的特点, 故以该模型模拟入水过程中的空化现象, 其对相间质量传递的描述如下

(9)

(10)

式中: 为蒸发率; 为凝结率; α_v为水蒸气相体积分数; n为单位体积; R_B为气泡半径。

4) 重叠网格技术

入水问题涉及到物体的大幅度运动, 需要把运动与流场耦合求解。重叠网格的优势在于可以实现物体的大幅度运动, 在运动过程中网格不会发生变形, 且可以清晰显示自由液面, 所以采用重叠网格技术将航行器自由运动与非定常多相流动耦合起来。重叠网格原理如图 2所示, 计算域网格分为背景网格区域和重叠网格区域, 各网格区域在空间上重叠, 但相互独立存在。背景网格随地面坐标系保持不动, 重叠区域可以随雷体坐标系在背景网格内部进行移动, 利用重叠部分的网格以一定的插值形式进行流场信息交互, 最终实现全流场计算。

1.3 边界条件与网格划分

图 3展示了背景计算域尺寸及其边界条件。背景计算域关于x₀O₀y₀平面对称, 长32L, 宽4L, 高19L, 水深15L, 空气域4L, 重力沿y₀轴负向。初始时刻航行器轴线与水平面夹角为20°, 航行体初速度为200 m/s。计算域左侧边界为速度入口, 前后边界为对称平面, 上、下、右侧边界为压力出口, 压强设置为随水深变化的压强。

入水点、前驻点、空化器等流动参数变化剧烈的区域全部位于重叠区域内, 因此重叠区域网格较密, 合理选择重叠区域范围能够提高对流场的模拟精度。经对比分析, 最终确定重叠区域范围如图 4所示, 其长6.6L, 宽4D, 高10D, 其中航行器头部距右边界0.6L。航行器位于重叠区域右视图中央, 其椭圆截面的短轴与x₀O₀y₀平面(铅垂面)重合。

采用切割体网格划分策略, 并在近壁面处添加棱柱层网格进行细化, 以期获得更精准的湍流模拟结果。在航行体头部、尾部附近流场设置加密区, 并根据理论计算公式所得的空泡轮廓对两相交界区域进行周向加密, 以获得更好的流场界面。网格划分结果如图 5所示。

1.4 数值方法验证

1) 与试验对比

为了验证数值计算的可靠性, 采用所建立的数值方法, 对文献中的试验工况开展仿真计算。首先对文献[20]的平头圆柱垂直入水试验进行数值模拟计算, 试验采用的平头圆柱长25.4 mm, 直径为12.56 mm, 质量为25.1 g, 初速度为603 m/s, 仿真计算结果与试验结果对比如图 6所示。

对陈晨等人^[21]的小型运动体斜入水试验进行数值模拟计算, 试验所用模型空化器直径为2.2 mm, 弹身直径6 mm, 质量8.5 g, 初速度141.15 m/s, 入水角为43.6°, 仿真计算结果与试验结果对比如图 7所示。

可以看出, 数值计算结果与试验结果相比有较好的一致性, 说明数值计算方法是可信的。

2) 结果无关性验证

为保证仿真结果独立性, 在入水速度200 m/s、入水角20°的工况下, 使用相同的数值模型进行计算, 分别进行网格无关性验证和时间步长无关性验证。在最小时间步长为1×10^-6 s的条件下, 选取不同的网格数进行仿真, 如图 8所示。

根据图 8a)可得, 网格节点总数为238万的算例满足网格无关性要求; 在238万网格数的条件下, 选取不同的最小时间步长进行数值仿真, 根据图 8b)可得, 最小步长为1×10^-6 s的算例满足步长无关性要求。因此, 在研究非回转体入水问题时, 选用238万网格和1×10^-6 s最小时间步长。

2 仿真结果与分析

使用第1节建立的数值模型, 在入水角20°、入水速度200 m/s, 无攻角、初始角速度和初始扰动的工况下, 对不同截面形状的非回转体航行器高速入水过程进行数值模拟研究。其中, 同一航行器的收缩段、柱段各横截面椭圆形状相似, 即截面椭圆长短轴之比相同, 非回转体航行器横截面的椭圆短轴与长轴之比分别取0.8, 0.7, 0.6, 对照组回转体的该比例为1.0。

2.1 入水过程非定常空泡特性

根据仿真计算结果提取非定常超空泡形态, 如图 9所示, 后处理时将水蒸气体积分数为50%的等值面作为两相交界。入水空泡经历了入水撞击、空泡形成、开口空泡和空泡闭合阶段。在入水撞击阶段, 航行器撞击液面, 形成不对称的空泡, 随后空泡迅速发展并包裹大部分航行器。

当入水9 ms时, k为1.0的航行器收缩段下表面刺破空泡壁面, 使得航行器下表面沾湿, 航行器未能稳定入水，发生弹跳翻转现象。当入水11 ms时, k为0.8的航行器尾部刺破空泡壁面, 沾湿区域扩大到一定程度后几乎不再变化, 此后航行器一直保持着在空泡壁面上滑行的姿态。k为0.7和0.6的航行器, 航行器尾部刺破空泡壁面, 随后在沾湿区域升力作用下, 航行器与空泡壁面分离, 航行器被空泡完全包裹, 空化器产生的抬头力矩使航行器尾部端面下沿再次刺破空泡壁面; 自此在流体动力综合作用下, 射弹发生了尾拍现象。

由图 9局部放大图可以看出, k为1.0的航行器, 收缩段率先沾湿, 升力作用点过于靠前, 沾湿区域无法提供恢复力矩, 故而航行器攻角会不断增大, 直至失稳; 当k为0.8时, 航行器柱段最早沾湿, 航行器柱段始终未能脱离空泡壁面, 在此过程中, 虽然收缩段也有所沾湿, 但沾湿区升力作用点仍在质心之后; 当k为0.6时, 航行器尾端最早沾湿, 且此时柱段距离空泡交界面还有相当的距离, 沾湿区升力作用点一直处于质心之后。

综上所述, 航行器收缩段及柱段截面的椭圆形状的变化会在很大程度上影响航行器刺穿空泡壁面的部位, 即改变了航行器的泡体位置关系。随着k值减小, 航行体收缩段愈发难以沾湿。

2.2 非回转体航行器的流体动力特性

根据仿真计算结果, 取不同航行器模型在入水过程中的流体动力参数进行对比如图 10所示。

入水初期, 不同短长轴之比的航行器所受的入水载荷相似。在0.16 ms时, 航行器空化器与水面开始接触, 阻力系数绝对值骤然升高。同时, 空化器受到的法向力使得航行器受到正向俯仰力矩作用, 引起航行器抬头。在3.2 ms时, 航行器尾部入水, 升力系数略微增大。尾部入水后, 空泡完全包裹航行器, 阻力系数趋于稳定, 约为-0.035, 此时航行器所受俯仰力矩很小。

在9 ms时, k为1.0的航行器收缩段发生沾湿, 航行器所受的力和力矩均急剧增大, 过大的流体动力载荷将对航行器结构安全产生不利影响。

在10 ms之前, 3种椭圆截面航行器所受的流体动力较为接近。在11.5 ms时, k为0.8的航行器尾部触水, 造成升阻力系数绝对值增大, 俯仰力矩系数跌至谷值。k为0.7和0.6的航行器先后触水, 触水后升阻力系数绝对值增大, 俯仰力矩系数跌至谷值, 随后航行器在沾湿区域升力作用下与空泡两相交界面分离, 升阻力系数逐渐回归触水前的数值。其中, 较之于k为0.7的航行器, k为0.6的航行器尾拍周期更短, 所受的力与力矩的振幅更高。相较于k为1.0的回转体航行器, 椭圆截面航行器的阻力系数绝对值均值更小, 升力系数和俯仰力矩系数峰值更小。

2.3 非回转体航行器的弹道特性

图 11给出了各截面形状航行器在典型工况下入水过程中的速度衰减曲线和质心位移曲线。

由图可知, k为1.0的航行器的速度衰减曲线在0.01 s内近似为一条斜直线, 但在0.01 s以后速度快速衰减, 这是由于航行器在这段时间里持续与空泡壁面大面积接触, 动能严重损失。而k为0.8, 0.7和0.6时, 航行器速度衰减趋势大体相同, 近似为斜直线, 这是由于k值减小使得航行器沾湿面积在较小范围内变化, 阻力系数因此受限。其中, 在入水初期, k为1.0的非回转体航行器速度衰减最慢, k为0.6的非回转体航行器速度衰减最快, 这是由于k越小航行器质量就越小, 故而速度衰减越快。k为1.0的航行器在5L位移之内, 质心轨迹近似为斜直线, 此后则不断向上偏转, 发生了跳弹现象; k为0.8, 0.7和0.6时, 航行器在入水后0.03 s内质心轨迹基本一致, 近似为一条斜线。

图 12分别为不同截面形状下非回转体航行器的俯仰角、俯仰角速度和攻角随时间变化曲线。

由图可知k为1.0的航行器, 俯仰角先减小后增大且增速越来越快, 俯仰角速度在5 ms后也表现出相同的趋势; 攻角在入水初期近似为线性增大, 在收缩段沾湿后, 升力作用点过于靠前, 难以提供足够的恢复力矩, 攻角增长没有被遏止, 航行器最终失稳。

k为0.8的航行器在尾部沾湿后, 其俯仰角速度显著下降, 俯仰角增速变慢, 攻角持续减小。k为0.7和0.6的航行器, 由于射弹尾拍, 故而俯仰角、俯仰角速度和攻角呈规律性的振荡状态, 总体来看, 其俯仰角速度和攻角幅值呈下降态势, 俯仰角增速也得到了抑制。其中, k为0.7的航行器在同一个振荡周期内俯仰角和俯仰角速度均值都更低。对比来看, 截面为一定椭圆形的非回转体航行器, 尾部沾湿区域的升力作用点始终在质心之后, 可以持续提供恢复力矩, 使得航行器攻角减小; k减小到一定程度时, 会发生尾拍现象, 此时俯仰角增速放缓程度更高; 但k继续减小则会导致升力系数和俯仰力矩系数峰值增大, 从而使俯仰角和俯仰角速度振荡的幅值增大, 导致尾拍频率加快、尾拍幅度增大, 这将影响航行器的结构安全。同时俯仰角增速放缓的程度有所减弱, 使得航行器姿态变化加快。

结合以上分析可知, 相比于回转体航行器, 具有椭圆截面的航行器尾部沾湿区域更为靠后, 故而滑水力作用点更易位于质心之后, 从而能产生恢复力矩维持航行姿态。航行器椭圆截面的k值越小, 其尾部距离空泡两相交界面的距离越远, 需要更大的攻角才会刺穿空泡, 故而刺穿的时间更晚; 同时, 由于航行器刺穿空泡前的角速度随时间增加, 故而k值越小会导致航行器拍击水的角速度更大、升力峰值更高、沾湿区提供的恢复力矩更大, 从而使尾拍周期更短、攻角振幅更大。总之, 减小航行器截面的短轴与长轴之比k能改变泡体位置关系, 影响航行器入水后下表面的沾湿区域, 进而增大恢复力矩, 有效改善该航行器跳弹倾覆的问题, 获得稳定的入水航迹, 但k值过小则会导致尾拍幅度过大, 俯仰角和攻角变化加快且幅度增大, 速度衰减加剧。

3 结论

建立了非回转体跨介质航行器斜入水过程仿真计算模型，研究了高速斜入水工况下，航行器横截面短长轴之比对多相流场、入水载荷、非定常运动特性的影响规律。主要结论如下：

1) 在不改变空化器的条件下，航行器横截面短长轴之比对入水空泡的影响有限，但是显著改变了泡体相对位置关系，短长轴之比越小航行体收缩段越难以穿刺空泡，短长轴之比不超过0.8时，航行器收缩段不会沾湿。

2) 选择合适的短长轴之比，可改善航行器在入水过程中的流体动力特性。适当减小短长轴之比，可降低流体动力载荷。但是，当短长轴之比小于0.6时，流体动力载荷明显增大，会对减阻效果和结构安全带来不利影响。

3) 随着航行器横截面短长轴之比的减小，入水弹道稳定性逐渐增强，主要原因在于航行器越扁，滑水力的作用点越靠后，尾部沾湿区域升力能够为航行器提供恢复力矩。短长轴之比不超过0.8时，入水弹道是稳定的；短长轴之比不超过0.7时入水弹道具有明显的“尾拍”特性。