梁景奇,徐保成,王 瑞,祁晓斌,李瑞杰

(西北机电工程研究所, 陕西咸阳 712099)

0 引言

超空泡航行器大部分表面被空泡包覆,航行阻力显著减小,航行速度大幅提高。超空泡射弹就是一种利用了超空泡减阻技术的动能武器,具备水下高亚声速甚至超声速航行的能力,可以对鱼雷、水雷等水下目标实现快速毁伤[1]。

超空泡射弹的重要特征之一是射弹的流体动力受到超空泡与弹体之间相对位置耦合关系的影响。当航速在300～1 200 m/s时,受到初始扰动后,包裹于超空泡内的射弹开始周期性地与空泡壁面发生碰撞,即尾拍运动。Kulkarni[2]等人基于Milwitzky[3]的理论采用加入附加质量的方法研究超空泡运动体的“尾拍”力;Euteneuer[4]建立了超空泡运动体的动力学模型,将后体滑行力作为外界高频扰动看待,忽略了空泡的时间延迟效应;王茂励[5]基于弹道仿真技术研究了流体动力参数的摄动对弹道稳定性的影响规律;魏英杰[6]等提出了超空泡射弹流体动力特性的动态描述方法,建立了射弹尾拍运动的无约束动力学方程;张广[7]等运用CFX建立了超空泡航行体三维6DOF数值模型,发现无控高速超空泡航行体具有运动稳定性;时素果[8]等开展了自由飞试验,获得了不同预置舵角下超空泡航行体水平运动的弹道特性。水下高速超空泡射弹的弹道稳定性与其尾拍运动具有重要关系,弹道失稳会造成射击精度和射程不足,影响射弹的作战效能。因此,对射弹空泡流场与弹体耦合的尾拍特性的进一步研究极具必要性。

依托CFD软件FLUENT18.0及其二次开发,结合动网格移动计算域技术,忽略重力因素,在考虑液体压缩性的基础上建立了高速射弹刚体运动与空泡流场耦合数值模型,从理论及实验方面验证了数值模型的准确性;求解超空泡射弹水下纵平面运动,探究水下运动稳定规律,分析了初始速度对射弹尾拍特性的影响。

1 数值模型

1.1 控制方程

采用VOF多相流模型模拟相界面运动。VOF是一种在固定Euler网格下的界面捕捉法,常用于由两种及以上不相混液体组成的流体中,适用于多相间有清晰界面的流动。

1.1.1 连续性方程

(1)

(2)

(3)

式中:ρ是流体混合密度;u是混合速度;n1是相数,由于考虑了不凝气体,在这里,n1=3;αk、ρk和μk分别为第k相的体积分数、密度和速度。

1.1.2 动量方程

匀质平衡流模型认为各相间不存在速度差,在流体微元中基于平均密度和平均动力粘度的动量守恒方程可描述为:

(4)

式中：ρm、μm分别为流体微元的各相平均密度和动力粘度,且有:

(5)

(6)

1.2 空化模型

采用Schnerr and Sauer空化模型模拟超空泡射弹的水下稳定空化绕流,Schnerr and Sauer空化模型将汽相体积分数和单位体积流体含有的空泡数量联系起来,表达式为:

(7)

(8)

式中:ρl为液体密度;ρv为水蒸气密度;ρ为混合相密度；pv为饱和为蒸气压;p为当前压力;αnuc为气核体积分数;RB为气核空泡直径;n2为单位体积内的空泡数量。

1.3 湍流模型

Realizabledk-ε湍流模型主要针对充分发展的湍流,稳定性好,近壁面区域的流动状态使用壁面函数预测,对边界层网格要求较为宽松;结合尺度化壁面函数,其在不增加计算量和保持模型稳定性的前提下增加了模型的适用范围,模拟效果好。

湍流强度k方程:

(9)

湍流耗散率ε的方程:

(10)

式中:μt为湍动粘度;μ是流体的时均速度;σk、σε分别为k、ε方程的湍流能量普朗特数;C1、C2为经验常数;E为时均应变率;ν为运动粘度;xi、xj为各方向距离。

1.4 可压缩液体模型

Tait方程是通过采用非线性回归的方法,对能够反应p-v-T三者关系的试验数据进行拟合,而得到的液体状态方程,广泛应用于描述可压缩液体的物性。没有温度修正的简化Tait液体状态方程可描述为:

(11)

K=K0+n3ΔP

(12)

Δp=p-p0

(13)

(14)

式中:p0为参考压力;ρ0为参考压力下的液体密度;K0为参考压力下的液体体积弹性模型;n3为密度指数;p为当前压力;ρ为当前压力下的液体密度;K为当前压力下的液体体积弹性模量;c为水中声速。

1.5 动网格耦合运动

耦合运动即射弹运动与空泡流场的计算相互耦合。已知前一时刻刚体重心位置和偏转角,通过对物体表面压力和剪切应力积分得到流体动力和力矩,再根据刚体运动方程计算物体运动的平移速度和角速度,然后重新计算重心位置和偏转角。

(15)

(16)

式中:G为变换矩阵,刚体的位置与方向根据每时间步线速度和角速度的变化而变化,c.g.表示质心。

采用移动计算域技术模拟弹体运动,计算过程中仅涉及计算域的移动,不存在网格的变形与重构,计算效率高,结果一致性好。

2 计算模型及方法验证

2.1 模型、网格划分及边界条件

文中研究的射弹模型如图1所示,射弹采用圆盘空化器,前部为两段锥段,中间为圆柱段,柱段尾部安装有6片尾翼。前端空化器直径3.2 mm,圆柱段最大直径15 mm,质量为0.23 kg。

图1 射弹模型图

采用圆柱形计算域,直径取50倍尾截面空泡直径,计算域轴向长度为11倍弹长,前端边界距离空化器4倍弹长,后端边界距离弹尾6倍弹长,该计算域径向尺度可以忽略空泡阻塞效应[9]。

针对所建立的三维计算域,采用ICEM软件的O-Block技术划分全结构化网格,如图2所示在弹体周围3 mm范围内的流域划分外O-block用于设置边界层网格,近壁面添加边界层网格,并根据y+值对网格进行优化。划分网格时特别注意在空泡两相交界面位置进行网格加密,最终划分的网格总数约80万,网格质量均在0.6以上。

图2 剖面网格划分图

计算域四周边界均设置为压力入口条件(pressure-inlet),设置静压值且静压值随深度变化,射弹表面的边界条件设置为壁面(wall),并且壁面与临界网格相对静止。计算域和边界条件设置如图3所示。

图3 计算域及边界条件图

2.2 数值模型验证

在超空泡射弹运动中,运动特性取决于流体动力特性,流体动力取决于弹体与空泡的相对位置关系,因此空泡形态计算准确性是运动模拟准确性的前提。将文献[10]经验公式和文献[11]试验结果与数值模型的空泡外形作对比,验证数值模型准确性。

2.2.1 经验公式验证

Savechenko[10]等在乌克兰国家科学院水动力实验室进行了空泡截面独立扩张原理有效性验证试验,给出了空泡外形描述经验公式:

(17)

(18)

(19)

式中:Rn为空化器半径;Rc,max、Lc,max分别为超空泡最大截面半径和空泡全长;x为空泡截面到空化器的距离;σ为空化数,σ=0.012～1.057。

对比同工况下经验公式和数值模型计算的空泡轮廓,验证数值模型的准确性。取水深1 m,航速100 m/s,空化器半径1.6 mm,前段空泡轮廓对比如图4所示。

图4 空泡轮廓对比图

图4对比结果显示,同工况下经验公式得到的空泡轮廓在径向尺度略大于CFD计算的空泡轮廓,计算误差小于12%,可以证明数值模型的准确性。

2.2.2 试验结果验证

Hrubes[11]采用试验方法,在水深4 m,速度970 m/s条件下水平发射平头锥形射弹,利用高速摄像机捕捉射弹入水后的空泡形态。试验结果如图5所示,射弹完全被空泡包裹。运用建立的数值模型,采用与文献[11]相同的弹体模型及初始条件,模拟射弹水下直航运动,空泡轮廓对比如图6所示。

图5 970 m/s初速试验与数值模拟空泡形态图

图6 数值模拟与试验空泡轮廓对比图

图6给出了Hrubes试验数据点和数值模拟的空泡外形曲线。从中可以看出,数值模拟结果与试验数据点基本吻合,均能有效包裹射弹,最大误差不超过10%,表明文中数值模型合理可行。

3 计算结果及分析

采用变参数法探究不同初速下超空泡射弹尾拍的流体动力特性及运动特性,水深恒定1 m,忽略重力、海流和波浪影响,考虑液体压缩性,模拟定深绕质心旋转的纵平面尾拍运动[12]。先定常计算初始空泡流场,再基于三维动网格模型并结合二次开发技术模拟自由减速运动,探究运动稳定性及流体动力、弹道参数的变化规律,初始条件如表1所示。

表1 初始条件表

3.1 流体动力特性

探究1 m水深、30 rad/s角速度扰动条件下不同初速射弹的流体动力变化特性,以计算时长0.015 s为计算终止条件,时间步长根据计算残差在1e-08到1e-05之间逐渐增加,计算结果如图7、图8所示。

图7 尾拍升力对比图

图8 尾拍阻力对比图

图7、图8分别为3种速度工况尾拍升、阻力对比曲线。初始扰动使弹体轴线转动,并与空泡壁面发生周期性碰撞,弹体的尾拍升、阻力均呈现周期性振荡特性。根据曲线特征,可将受扰动后超空泡射弹的振荡过程分为发散阶段和收敛阶段。在第一阶段,由于弹体穿刺空泡深度不断增加,尾拍升、阻力曲线均呈发散特性,振荡幅值迅速增大并达到最大,振荡周期不断减小;在第二阶段,由于速度衰减,升、阻力曲线呈收敛特性,振荡幅值逐渐减小,振荡周期缓慢增大,两个阶段升、阻力曲线的幅值和周期均同步变化,同一时刻升力幅值约为阻力幅值的3.5倍。

尾拍力与速度的平方正相关,初速越大,振荡周期越小,尾拍升、阻力越快达到最大幅值,最大幅值越大,均值越大;随着速度的衰减,3种工况下阻力和升力均值的差距逐渐减小;0.01 s后时间步长增至1e-05,振荡曲线呈阶梯形,但依然能客观反映尾拍力振荡特性。

3.2 空泡形态

以1 200 m/s工况为例,分别取振荡发散阶段和振荡收敛阶段的一个尾拍周期T,分析尾拍过程中空泡形态变化及弹体与空泡的位置关系,结果如图9、图10所示。

图9 振荡发散阶段尾拍动态图(第一周期)

图9为振荡发散阶段(第一个尾拍周期)内的空泡动态变化过程。弹体向下偏转,尾翼局部穿刺空泡下壁面后产生二次空泡,升、阻力达到最大,弹体回弹并触碰上壁面,最终回到起始位置。在该阶段,尾拍角幅值随穿刺深度的增加呈增大趋势,尾拍角为0°时,弹体除空化器外不沾湿,空泡轴线与弹体轴线基本重合。

图10为振荡收敛阶段(第6个尾拍周期)的空泡动态变化过程,在该阶段,弹体穿刺空泡深度更大,空泡形态发生显著改变。到达最大尾拍角时,除尾翼外,弹体第二锥段和靠近尾翼的圆柱段也发生沾湿,且两次沾湿区域具有不对称性,而尾拍角为0°时弹体依然存在沾湿情况,这是由于空泡变化与弹体振荡周期存在相位差,空泡轴线变化滞后于射弹轴线变化,即空泡时间延迟效应。

3.3 运动特性

以地面系为参考坐标系,分析3种初速射弹受扰动后的尾拍运动特性,俯仰角、俯仰角速度、轴向速度、质心位移对比曲线分别如图11～图14所示。

图11、图12为俯仰角和俯仰角速度随时间变化曲线。不同速度下的俯仰角、俯仰角速度曲线均呈周期性振荡特性,同样可分为振荡发散和振荡收敛阶段,变化规律与流体动力曲线相似。由俯仰角曲线可知,初速越大,振荡周期越小,振荡发散阶段俯仰角幅值越小;在振荡收敛阶段,3种工况俯仰角幅值均接近4°,基本保持恒定,这表明,振荡收敛阶段弹体俯仰角幅值受初速影响极小。

图11 俯仰角随时间变化曲线

图12 俯仰角速度随时间变化曲线

由俯仰角速度曲线可知,初速越大,俯仰角速度曲线振荡周期越小,越快达到最大幅值,最大幅值越大,1 200 m/s工况最大俯仰角约为270 rad/s,增幅高达9倍;达到最大幅值后,俯仰角速度振荡幅值均随时间逐步衰减。通过与图7、图8对比可知,尾拍力曲线俯仰角曲线存在相位差,亦是由空泡延迟效应所致。

图13为X向速度随时间变化曲线,3条曲线均呈单调递减趋势,初速度越大,阻力越大,因此速度衰减越快,0.014 s内1 200 m/s、1 000 m/s、800 m/s工况的速度衰减幅度分别为260 m/s、190 m/s和130 m/s;随着时间的增加3种工况航速差距逐渐缩小,由此可知通过增加初速增加射程效果有限。

图13 X向速度随时间变化曲线

图14为0.015 s内3种工况射弹质心位移变化曲线,质心随着弹体的振荡均产生了侧向位移,且位移量不断增加。以1 200 m/s为例,轴向位移15 m时,侧向位移为0.02 m,与轴向位移量相比,侧向位移量可忽略不计,弹体基本保持了直线航行,这表明超空泡射弹的水下弹道具有动态稳定性。

图14 质心位移变化曲线

4 结论

文中基于建立的流场运动耦合数值模型,数值仿真了超空泡射弹的尾拍运动,分析了初速对尾拍力的特性及运动特性的影响,结论如下:

1)受扰动后超空泡射弹的水下运动具有动态稳定特性,尾拍力、俯仰角、俯仰角速度等参数均呈现准周期振荡特性,可将该过程分为振荡发散阶段和振荡收敛阶段；

2)尾拍运动使空泡轴线产生偏移,空泡轴线变化滞后于射弹轴线变化,同一尾拍周期的两次沾湿区域具有不对称性；

3)初速度越大,尾拍周期越小,弹体尾拍升、阻力及俯仰角速度振荡幅值越大,俯仰角振荡幅值受初速影响极小。