1 引言
潜射航行体的发射是一个涉及多因素干扰，目前为应对各种复杂环境需求，航行体朝着变深度发射、简化结构、提高发射可靠性等趋势发展，其中发射可靠性一直是其重点研究方向之一[1]。以往的研究往往只关注水下过程[2-4]，忽略了出水过程的阻力动态演变。因此出水过程的阻力来源对于水下航行体的减阻设计是非常必要的。利用实验方法与重叠网格技术结合，揭示了航行体出水过程中水冢现象、浮力变化以及附连水质量变化对于水动力及水动力系数影响的内在规律，为工程设计提供了指导。
2 数学模型
2.1 VOF模型
VOF模型中动量方程为求解整个区域内统一的动量方程，作为结果的速度场由各相共享[5]，即：
ρg+F
(1)
相间界面跟踪通过求解一项或多项容积比率的连续方程来完成。对第q相方程为：
(2)
式(2)中：
为第q相对第p相的质量传输；
为p相对第q相的质量传输；Sαq为源项；αq为容积比率。
2.2 水动力系数
水动力系数是水下航行体结构设计的重要参数，一般而言，水动力系数大多通过水动力测量试验和理论方法计算相结合获得，由于航行体水动力系数很多，想要全部辨识是不切实际的，本文选取最常用的无量纲化水动力系数来建立数学模型[6]，有：
Cd=F/(0.5ρV2S)
(3)
式(3)中：ρ为水的密度；V为水深z处航行体最大纵向速度，即为质心速度在纵向方向的分量；S为航行体特征横截面积；L为航行体轴向长度。
航行体质心的流体力合力为：
(4)
式(4)中：
为表面压力积分项；
为黏性力积分项。
3 实验方案
实验前航行体与发射平台在室内进行吊装，吊装过程如图1所示，吊装完成后,将航行体和发射平台移至湖内水下，竖直固定在湖内水下。
图1 航行体发射平台室内吊装工作示意图
Fig.1 Schematic diagram of indoor hoisting work of vehicle launch
为获得航行体表面压力，实验时在航行体表面布置了多个压力监测点，测点1为航行体头部驻点，测点2分布于头部母线段，测点3分布在柱段，测点4位于尾段。航行体几何尺寸及测点具体位置如图2所示。
图2 航行体几何尺寸及测点位置
Fig.2 Geometrical dimensions of the vehicle and location of measuring points
试验当天无风，湖面基本无浪，航行体从水下15m处垂直起射直至射出水面。其水下弹射出筒过程如图3所示，射出水面过程如图4所示。
图3 航行体水下垂直弹射过程
Fig.3 The underwater vertical ejection process of the vehicle
图4 航行体出水过程
Fig.4 The process of exiting the water
实验中航行体从初始体到完全出水等过程时序如图5所示。航行体尾部离开发射筒时，高压发射气体在冲出筒口的过程中将会对流场产生了剧烈扰动，形成筒口效应[7]，这将会使得此时的数据误差较大、可信度低，数据分析时要避免这一段的筒口效应影响，因此本文数值模拟重点为模拟筒口效应消失到航行体完全出水这一段过程，并对这一过程的流场以及水动力系数进行分析讨论。
图5 出水过程时序图
Fig.5 Timing diagram of water outlet process
4 数值模拟
采用重叠网格技术，利用k-epsilon(2eqn)湍流模型、VOF模型计算当前时刻瞬态流场，通过对航行体表面压力积分，获得当前时刻航行体所受水动力，在此基础上计算航行体的水动力系数等; 利用迭代获得下一时刻的航行体边界和网格，进而求解下一时刻的瞬态流场，接着计算下一时刻航行体的加速度、速度和位移等得到航行体运动特性。上述过程不断循环，直至航行体完全出水为止。
为简化模型和计算，做如下假设：
1)流体常物性，不可压缩；
2)流体在航行体运动之前是静止的；
3)航行器是刚体，其外形关于xOz、yOz平面对称；
4)不考虑其偏航运动；
5)不考虑空泡的影响。
4.1 流体控制方程
流体力学的控制方程为3个基本的守恒方程，即连续性方程(质量守恒方程)、动量守恒方程(运动方程)和能量守恒方程，其中能量守恒方程是包含有热交换的流动问题必须满足的基本方程，因为本文所做数值计算中涉及到的流体(海水)是不可压缩流体，故将其忽略，数值计算时只需要考虑质量守恒方程和动量守恒方程即可。
考虑到本文中的流体处于湍流状态，所以还需遵守附加的湍流输运方程。Menter[8]和Rattanasiri[9]等人分别将k-ω和k-ε模型应用于近壁面和远流场区。本文应用的湍流模型为SST k-ω湍流输运方程。
连续性方程：
非定常N-S方程[10]：
(5)
湍流动能方程[11]：
(6)
湍流动能耗散率方程：
(7)
式(5)～(6)中：μ为流体粘度；σk为湍流动能普朗特数；σω为湍流耗散普朗特数；Gk和Gω分别为k和ω的耗散；Sω为湍流耗散率的源项；Dω为交叉扩散项；μt为湍流黏度。
4.2 几何模型及边界和初始条件
将整个计算域划分为2个区域：一个是数值波浪水池区域(背景区域)；另一个为航行体区域(子区域)。航行体区域设置相对比较简单，对于包裹体(航行体的重叠区域)的外部边界设置为重叠网格边界。由于没有波浪和海流，因此不需要考虑波浪和海流的进出口，边界设置相对简单，波浪水池四周设置为对称面，底部设置为滑移壁面,顶部设置为压力出口。边界条件设定如图6所示。对于模拟输入的初始条件，我们采用了实验测定的0 s时的速度值，即：Vz0=25 m/s。
图6 边界条件设定
Fig.6 Boundary condition setting
4.3 重叠网格模型
重叠网格的划分主要分为2个区域，即背景区域和航行体区域(子区域)。背景区域在自由表面附近划分很细的网格，用来捕捉航行体穿越自由液面时自由液面的细节变化，在竖直方向划分出与子区域尺寸一致的切割体网格，用来更好地与子网格装配。背景计算域网格划分完成后，对子区域进行网格划分，并对其近壁面网格细化。子网格与背景网格重叠，子网格通过插值向背景网格传递计算结果，加快计算效率和提高计算精度[12]。网格划分如图7所示。
图7 背景区域和航行体区域网格划分
Fig.7 Meshing of background area and vehicle area
5 结果与分析
将实验数据与仿真计算的结果整理，绘制速度曲线如图8所示，航行体表面压力曲线图如图9所示。从速度曲线可知，航行体在水下速度下降显著，水下0.88 s内速度实验值由25 m/s下降至9.5 m/s，计算值由25 m/s下降至9.7 m/s。从压力曲线分布可知，头部驻点压强随着航行体速度减小而降低；头锥段压力，因为绕流作用和环境压力的下降而呈现出较稳定的下降趋势，0.88 s后头部出水，驻点压强呈现明显的下降趋势，出水前达到最大压力为0.3 MPa。仿真计算结果与实验吻合较好，验证了用重叠网格方法处理航行体运动的准确性和可靠性。
图8 航行体运动速度曲线
Fig.8 Velocity velocity curve
图9 航行体表面压力曲线
Fig.9 Surface pressure curve of vehicle
取航行体在水下0.7 s截面压力云图如图10所示，由图10可以看出，航行体运动对整个流场压力的分布影响很小，其主要影响集中在航行体周围，可以看出明显的高压和低压区。其中最明显的高压区包括航行体顶部高压区，达到0.386 MPa，是阻力的主要来源，低压区主要表现在航行体肩部附近的低压区和尾端，分别达到0.184 MPa和0.132 MPa。
图10 0.7 s时航行体表面压力云图
Fig.10 Surface pressure cloud map of the vehicle at 0.7 s
航行体出水过程数值仿真如图11所示，其中t0指的是航行体头部出水时刻；t1是肩部出水时刻；t2是因为肩部到主段之间承接处出水时刻；t3指的是航行体尾部刚离开自由液面时刻；t4时刻拐点是航行体出水后完全摆脱附连水时刻。由图11可知，航行体出水产生水冢，还会穿越水面而形成水冢破碎[13]，之后随着航行体出水的距离增加，水冢还会增长一段时间，在达到最大水冢之后，水面开始沿着航行体表面回落，完全出水后仅有少量附连水附着于航行体表面，附连水引起的粘滞力是航行体完全出水后主要的流体阻力来源[14]。
图11 航行体出水过程数值仿真
Fig.11 Numerical simulation of the water exit process of the vehicle
同时，绘制水动力曲线如图12所示，水动力系数曲线如图13所示。其中纵轴代表力的大小和方向。由图12、图13可以看出，从航行体头部到完全出水过程中，航行体水动力及系数在t0时刻达到最大值，分别为0.65 kN和0.184。在整个过程中，水动力及力系数主要是沿Z轴方向变化明显，X方向波动很小，基本可以忽略。在航行体出水过程中，水动力和力系数波动变化拐点主要出现在整个航行体出水过程，水下段水动力和力系数变化相对稳定。但从筒口效应到出水过程中，水动力及系数在0 s最大，这主要是因为水的静压所导致的。
图12 水动力曲线
Fig.12 Hydrodynamic curve
图13 水动力系数曲线
Fig.13 Hydrodynamic coefficient curve
结合图11，分析航行体出水过程水力和力系数变化曲线上各个拐点出现的原因。t0时刻出现的拐点是因为航行体头部出水导致的，这是航行体出水过程的开始时刻，水动力以及水动力系数达到最大值(从头部接触水面到完全出水)；t1时刻出现的拐点是由于肩部出水导致受力发生变化；t2时刻拐点是因为肩部到尾部之间有段承接处几何形状改变导致的；t3时刻拐点是由于航行体尾部出水，这是航行体出水过程的结束时刻。此时航行体已经出水，所受浮力变为零，出水时刻速度较低，其空气阻力又特别小，可以忽略不计，但此时航行体水动力不为零，这是由于航行体出水时其表面所附带的水沾湿产生的粘滞力导致的，粘滞力是航行体刚出水后主要的流体阻力来源[15]。t4时刻拐点是因为航行体出水后完全摆脱水的影响，而空气阻力又特别小，可以忽略不计，因此所受到的流体力及力系数趋于零。
6 结论
1)航行体在水下时顶部压力是航行体阻力的主要来源力。水动力及系数呈现较为稳定的上升趋势，而速度下降显著，计算值由25 m/s下降9.7 m/s。
2)出水过程水动力及力系数主要受水冢现象、浮力变化以及附连水质量变化影响而呈现波动并逐渐趋向于零。
3)粘滞力是航行体完全出水后主要的流体阻力来源。
[1] 倪火才.潜载导弹水下发射技术的发展趋势分析[J].舰载武器,2001(01):8-16.
Ni H C.Analysis on the development trend of underwater launch technology of submarine-borne missiles[J].Shipborne Weapons,2001(01):8-16.
[2] Dai Y,Xue C,Su Q,et al.Numerical analysis on hydrodynamic characteristics of a deep-sea mining vehicle under three typical motions[J].Ocean Engineering,2021(235):109446.
[3] Hanssen F C W,Greco M.A potential flow method combining immersed boundaries and overlapping grids:Formulation,validation and verification[J].Ocean Engineering,2021(227):108841.
[4] Shah S M H,Mustaffa Z,Matínez-Gomariz E,et al.Hydrodynamic effect on non-stationary vehicles at varying froude numbers under subcritical flows on flat roadways[J].Journal of Flood Rrisk Management,2020,13(04):12657.
[5] Akhlaghi M,Mohammadi V,Nouri N M,Taherkhani M,et al.Multi-fluid VOF model assessment to simulate the horizontal air-water intermittent flow[J].Chemical Engineering Research and Design,2019(152):48-59.
[6] Panda J P,Mitra A,Warrior H V.A review on the hydrodynamic characteristics of autonomous underwater vehicles[J].Proceedings of the Institution of Mechanical Engineers Part M.Journal of Engineering for The Maritime Environment,2020,235(01):15-29.
[7] 崔可.航行体出筒以及水下运动的水动力研究[D].哈尔滨：哈尔滨工业大学，2012.
Cui K.Hydrodynamic research on the exit of the vessel and the underwater movement[D].Harbin:Harbin Institute of Technology,2012.
[8] Menter F R.Two-equation eddy-viscosity turbulence models for engineering applications[J].AIAA J,1994,32(08):1598-1605.
[9] Rattanasiri P,Wilson P A,Phillips AB.Numerical investigation of a fleet of towed AUVs[J].Ocean Engineering,2014(80):25-35.
[10] Mittal K,Dutta S,Fischer P.Multiratetimestepping for the incompressible Navier-Stokes equations in overlapping grids[J].Journal of Computational Physics 2021(437):110335.
[11] Wilcox D C.Turbulence modeling for CFD[M].LaCanada CA:DCW Industries,1998.
[12] 周超杰,洪亮,张周康,等.重叠网格在多浮体结构CFD中的应用[J].兵器装备工程学报,2018,39(11):199-204.
Zhou C J,Hong L,Zhang Z K,et al.Application of overlay grid in multi-floating structure CFD[J].Chinese Journal of Ordnance Equipment Engineering,2018,39(11):199-204.
[13] 孙士丽,孙义龙,胡竞中,等.圆柱体出水运动的自由面效应及水冢现象分析[J].计算物理,2013,30(02):187-193.
Sun S L,Sun Y L,Hu J Z,et al.Analysis of free surface effect and water mound phenomenon of cylinder water outlet movement[J].Computational Physics,2013,30(02):187-193.
[14] Pan D,Xu X,Liu B.Influence of flanks on resistance performance of high-speed amphibious vehicle[J].Journal of Marine Science and Engineering,2021,9(11):1260.
[15] 杨晓光,陈焕龙,刘华坪,等.导弹水下运动及出水过程的三维流场仿真[J].弹道学报,2010,22(03):107-112.
Yang X G,Chen H L,Liu H P,et al.3D Flow field simulation of Missile’s underwater movement and water exit process[J].Journal of Ballistics,2010,22(03):107-112.