水下航行体垂直出筒流体动力特性研究

刘海军，王 聪，邹振祝，王本利，王柏秋

(哈尔滨工业大学航天学院，150001哈尔滨)

水下航行体的发射方式可分为有动力发射和无动力发射.根据布置方式可分为垂直发射、水平发射和倾斜发射;根据发射筒内航行体周围的环境分为干式发射和湿式发射;按其发动机的点火时机分为筒内点火、水下点火和水面点火等［1］.无论采用哪种水下发射方式，航行体都要经历出筒、水中航行、出水和空中飞行4个阶段，尽管前3个阶段的运动距离和时间比较短，但会形成复杂多相流场，对航行体流体动力、弹道及载荷等产生直接影响，从而影响航行体运动的稳定性［2］.

潜射航行体在水下运动的时间很短，却存在复杂的流场情况，殷崇一等［3］运用有限体积法和SIMPLE算法计算了航行体水下发射内流场的数值研究.陈雄州等［4］根据固体发动机推力最大原则，并结合水下航行体外弹道运动方程，求出了不同使用环境下的发动机最佳出口压强.卢立秀等［5］运用MATLAB对飞机内弹道弹射机构进行了研究.王汉平等［6］采用多相流模型和动网格技术研究了均压气体对航行体射后效的影响.仲维国等［7］基于鱼雷的航行力学导出水下航行体的弹道方程.Acosta等［8］采用数值模拟的方法对重力场中的轴对称细长体沿轴线的定常空泡长度、空泡截面面积和阻力系数进行了研究，得到了重力场中弗洛德数对空泡的长度和面积影响较大，对轴对称细长体的阻力系数影响较小.Basharova等［9］采用基于细长体空泡理论的近似方程组计算了有重力影响无限流场中沿轴线方向上空泡的形态和基本尺寸.J.J.Yagla等［10］对水下垂直发射系统的动力环境进行了研究，提出了筒内发动机点火的方式可以为航行体的出水提供1个从出筒到出水整个过程的气幕，并进行了相关的实验.Chris J.Weiland等［11］研究了不同水深和横向流对气幕的影响.刘筠乔等［12］研究了航行体垂直发射出筒过程的通气空泡流.张红军等［13］对不带空化模型的潜射航行体垂直出筒过程进行了研究.刘乐华等［14］采用TVD有限体积法研究了深水发射的内流场，得到燃气流喷出之后快速膨胀，进而形成高温燃气射流，对反射板有极大的破坏作用.朱明骏等［15］对潜射航行体的充气均压系统进行了研究，并利用PID控制方法对充压系统进行了设计.

综上所述，国内外只是对潜射航行体的垂直发射出筒过程的内流场进行了初步研究，没有对航行体出筒过程中的流体动力进行研究，由于航行体在出筒过程中受到的流体动力对航行体的出筒及出筒后飞行的弹道都有影响，故对航行体在出筒过程中的流体动力特性进行研究是必要的.本文采用多相Mixture模型并结合动网格技术对运动流场中混合介质RANS方程进行求解，实现了重力场中运动固体边界与气(汽)水流场的耦合求解，针对采用热发射方式的航行体垂直发射出筒过程的多相流场进行了数值模拟，研究了重力影响下的水下航行体垂直出筒过程和流体动力特性.

1 数学模型

1.1 控制方程

描述水下航行体垂直出筒过程的控制方程主要包括混合物的连续方程、动量方程、能量方程、湍流模型和空化模型.

混合物连续性方程为

式中:ui为i方向的速度分量;ρm为混合物的密度，ρm=αlρl+αvρv+αngρng，ρi和αi(i=l，v，ng)分别为水、水蒸气及不可结气体的密度和体积分数;fi=αiρi/ρm(i=l，v，ng)分别为液体、水蒸汽及不可凝结气体的质量分数;不可凝结气体为空气;m-是水的汽化源相，m+是水蒸汽液化过程的源相.

混合物动量方程为

式中:p是混合介质的压力;μm是运动粘性系数，其中n是构成混合介质的相数，μk是每个相的湍流粘性系数;ρm为混合物的密度;g是重力加速度.

1.2 湍流模型

标准的k-ε湍流模型为

湍动能k和湍动耗散率ε的输运方程如下:

式中:k是混合介质的湍动能;ε是混合介质湍流耗散率;Gb是浮力产生湍动能项;Gk是速度梯度产生的湍动能项;μt为湍流粘性系数;ρm为混合物的密度;σk和σε分别为k和ε的普朗特数;可以调整的经验常数取值分别为:Cε1=1.44，Cε2=1.92，σk=1.0，σε=1.3，Cμ=0.09.

1.3 空化模型

基于均质平衡流理论的Singhal空化模型，表达式如下:

式中:Vch是液相和汽相之间转化的特征速度;τ是液体表面张力;pv是饱和蒸汽压;fl=1-fv-fng，其中不可凝结气体分数fng=1.5×10-5;Cc=0.01和Ce=0.02是两个经验常数.

1.4 能量方程

能量方程为

式中:keff为传导率，keff=(k+kt);kt为湍流的热传导率为物质j的扩散通量.上式中等号右边的最后3项分别代表由热传导、物质扩散和粘性耗散引起的能量交换.Sh包含化学反应热在内的一切热源.

1.5 航行体的运动学方程

航行体的运动学方程由牛顿第二定律得到，对航行体运动学方程使用下式进行离散:

式中:Vi是航行体在第i个时间步的速度;Vi-1是航行体在第i-1个时间步的速度;Δt是时间步长;Fi-1是航行体在第i-1个时间步合外力，m是航行体的质量.合力由下式算得:

式中:FG为航行体受到的重力;FB为航行体的底部受到的推力;FG和FB由初始条件给出;FP是航行体前后压强差产生阻力;FV是航行在粘性流体运动受到的粘性力阻力，FP和FV均是通过对流场的计算得到.

2 数值模拟

2.1 数值求解及动网格

使用SIMPLEC方法进行计算，采用PRESTO对控制方程的压力项进行离散，采用二阶迎风格式对控制方程中的能量方程和湍流耗散项及对流项进行离散.利用动网格技术解决运动边界和计算域变化的问题.编制控制航行体运动的程序，通过自定义函数(UDF)嵌入FLUENT，求解出航行的速度和位移.

采用动网格技术进行计算时，动网格区域中的网格具有一定对流运动，在实际计算过程中要去掉动网格对流运动的影响，动网格守恒方程统一形式如下:

式中:ρm为混合物的密度;u为流体的速度矢量;ug为网格的运动速度;Γ为扩散系数;Sφ为标量φ的源相;∂V为控制体积V的边界.

2.2 数值求解方法及正确性

本文采用间接证明的方法，对比分析数值模拟与实验结果，得到数值求解方法的正确性.验证算例采用模型为:半球头型航行体直径为50 mm、总长度为200 mm.计算条件是:未受扰动水流场压强p∞=378 955 Pa，采用非定常过程计算.以航行体直径D计算雷诺数为Re=106，空化数σ=0.31，图1给出了算例中的航行体在尾部离开筒口时刻的物面压力系数分布曲线与文献［16］中实验数据对比曲线图.H为距离航行体的头部的轴向距离，从图中看出仿真结果和实验数据吻合较好.

图1 实验结果与仿真结果的比较

在计算过程中采用网格重构法，结合多相流Mixture模型和湍流模型，对流场混合介质RANS方程进行求解，进而求解固体运动边界与汽水流场的流-固耦合问题，采用此方法对水下航行体垂直发射出筒过程流体动力特性展开数值模拟研究.

3 计算模型及边界条件

航行体的长细比L/D=10，L为航行体的长度，D为航行体的直径，头为半球型，发射深度为50 m.计算域由筒内的气体和筒外的水域组成.图2给出了流场计算区域、网格和边界条件的示意图，在航行体附近的网格进行了加密处理.图3给出了静网格和动网格示意图，整个计算域分为两块:中间为动态区域，网格随时间变化，周围为静态区域，网格固定不动.流场的长度是航行体长度的5倍，宽度是航行体直径的12倍.

航行体喷管入口是本模型的压力入口，喷管入口的压力采用图4给出随时间变化燃气总压曲线.

图2 计算区域与网格

图3 计算区域与模型

图4 燃气总压随时间变化曲线

4 计算结果及分析

4.1 流场分析

图5给出了水下航行体不同时刻的流体相体积分数分布云图.I至III是航行体头部没有穿透气团之前的气团的生成过程，此时气团是半球型.III至VI是航行体头部穿透气体时气团的演化过程，此时气团演变成椭球型，因航行体的速度逐渐增大和环境的静压力逐渐减小，空化数逐渐减小，故在航行体的肩部产生空化，并逐渐发展成包裹航行体肩部的空泡.航行体离开发射筒后，筒内的压力减小，但是仍高于周围环境的压力，筒口的气团在海水的作用下开始径向回缩.

图5 不同时刻流场体积分数ai云图

4.2 流场压力分布

图6给出了不同时刻流场的压力云图，所选的时间点与图5一一对应，从压力云图中可以看出，水域流场的压力分布主要受重力的影响.初始时刻，筒口附近水域压力比较大，发动机点火筒内的压力开始增加.随着出筒高度的增加周围环境压力开始减小.从IV时刻后，在航行体的肩部出现了绕流引起的低压区，随出筒高度的增加，航行体的速度增大，在肩部绕流引起的低压区继续增大，这与图5的压力云图中的肩部空泡的产生是对应的.航行体尾部离开发射筒后，筒口附近的海水压力沿径向往筒口逐渐增加，但筒口气团的压力仍高于环境压力，此时筒口周围的水暂时没有向筒内回灌.

图6 不同时刻流场压力云图

4.3 航行体表面压力分析

图7给出了不同时刻航行体表面的环向压力沿轴向的分布，横坐标代表的是航行体的特征长度，纵坐标代表的是表面压力与最大压力的比值，h是距离航行体头部的高度，L是航行体的长度，T=L/V0，V0是出筒的速度，T是特征时间.左侧是航行体的头部，右侧是航行体的尾部.从图中看出，随着出筒时间的增加，航行体表面的压力系数是逐渐降低的，在航行体的头部出现了由绕流引起的低压区，在航行体的尾部由于处在发射筒内高压燃气之中，出现了高压区，随出筒时间的增加航行体尾部压力逐渐减小并趋缓，这与图6中压力场的分布一致.

图7 不同时刻航行体表面压力比值沿轴向分布曲线

4.4 航行体出筒弹道分析

图8 航行体速度和位移变化曲线

图9是加速度变化曲线.加速度曲线在初始刻有上下波动，这是由于发动机的推力在逐渐增加至峰值时，加速度也迅速增大，然而在筒口，由于气团的膨胀和收缩导致加速度又有所波动，最后加速度抖动下降.

4.5 航行体出筒过程受力分析

图10给出了压差力系数、粘性力系数和总流体动力系数变化曲线.由图9可得出，在t＜0.375T时压差力系数、粘性力系数和总流体动力系数的都在对应的时间上有波动.压差系数和总流体动力系数波动的趋势一样，产生的原因是:1)发动机的推力随着时间的增加逐渐增大，并且到达峰值;2)气团在没有被航行体穿透时，气团在水静压的作用下收缩和膨胀，产生瞬间压力的变化.t＞0.375T时，压差系数曲线逐渐下降且趋缓.粘性力系数在t＜0.375T时的波动，是由于海水静压作用下气团的收缩和膨胀影响高压燃气的速度，此时气体和航行体相对速度发生改变，进而影响粘性力系数，粘性力系数为正值.航行体受到的粘性力由高速气体对航行体的粘性力和水对航行体的粘性力组成，两者的作用力方向相反，水的粘性力随着沾湿面积的增加而增大，因此总粘性力系数逐渐减小.这说明了压差系数是影响总流体动力系数的主要因素，粘性系数对总系数的影响是次要因素.

图9 航行体加速度变化曲线

图10 总系数变化曲线

5 结论

采用数值模拟的方法研究了有重力影响的水下航行体垂直发射出筒过程流体动力特性，得到如下结论:

1)航行体头部距离筒口0.01D处，气团的膨胀和收缩变化剧烈，加速度出现剧烈波动，距离筒口高度的增加，筒口气团的膨胀和收缩变化趋缓;

2)气团的形状从半球型变化到椭球型，最后变成细长型;航行体头部在海水静压作用下出现了高压和在肩部出现了绕流引起的低压区;

3)航行体头部穿透气团之前，气团的膨胀和收缩引起航行体的流体动力系数较大的波动，缩短航行体穿透气团的时间可以减弱航行体载荷的波动，进而优化航行体的水弹道.

［1］颜开，王宝寿.出水空泡流动的一些研究进展［C］//第二十一届全国水动力学研讨会.北京:海洋出版社，2008:9-16.

［2］黄寿康.流体动力弹道载荷环境［M］.北京:宇航出版社，1991:165-182.

［3］殷崇一，张宇文，刘乐华，等.导弹水下发射内流场的数值研究［J］.弹箭与制导学报，2003，23(3):56-58.

［4］陈雄州，舒旭光，廉斌，等.水下火箭发动机喷管出口压强研究［J］.舰船科学技术，2004，26(6):38-40.

［5］卢立秀，汤军社，刘永超，等.导弹弹射机构的建模及优化设计［J］.航空计算技术，2007，37(6):58-61.

［6］王汉平，吴友生，程栋，等.潜射模拟弹弹射后效分析［J］.船舶力学，2010，14(10):1122-1128.

［7］仲国维，张嘉钟.潜射航行体的水下弹道模拟［J］.弹道学报，2005，17(1):8-12.

［8］ACOSTA A J.The effect of a longitudinal gravity field on the supercavitating flow over a wedge［J］.Trans ASME，1961，28(2):188-192.

［9］BASHAROVA V N，BUIVOL V N，SEREBRYAKOV V V.Slender axisymmetric cavities in the flow around bodies in a longitudinal gravity force field［J］.International Applied Mechanics，1983，19(4):359-366.

［10］YAGLA J J.Launch dynamics environment of a water piercing missile launcher［C］//Proceedings of the 24th international symposium on ballistics.New Orleans，LA:NDIA，2008:1-17.

［11］WEILAND C J，VLACHOS P P.Concept analysis and laboratory observations on a waterpiercing missile launcher［J］.Ocean Engineering，2010，3(9):959-965.

［12］刘筠乔，鲁传敬，李杰，等.导弹垂直发射出筒过程通气空泡流研究［J］.水动力学研究与进展(A辑)，2007，22(5):549-554.

［13］张红军，陆宏志.潜射导弹出筒过程三维非定常数值模拟研究［J］.水动力学研究与进展(A辑)，2010，25(3):405-415.

［14］刘乐华，张宇文，袁绪龙.水下大深度垂直发射内流场的数值研究［J］.水动力学研究与进展(A辑)，2005，20(1):90-94.

［15］朱明骏，杨军.潜射导弹充气均压系统建模与控制系统设计［J］.战术导弹控制技术，2010，27(1):26-30.

［16］ROUSE H，MCNOWN J S.Cavitation and pressure distribution，head forms at zero angle of yaw［R］.Iowa City:State University of Iowa，1948:23-26.