黄宝珠,李代金,黄 闯,古鉴霄,罗 凯
(西北工业大学 航海学院,陕西 西安,710072)
超空泡技术是可以使水下航行器获得减阻90%的革命性技术[1]。利用该技术,水下高速射弹航行时弹体表面出现低压区,当压力低于当地水的饱和蒸汽压时,出现局部空化,随着局部空泡进一步发展成超空泡,射弹表面几乎被完全包裹,此时射弹水下航行的摩擦阻力大大减小,实现增速增程[2]。
水下超空泡射弹在300~1 000 m/s 航速下的运动稳定模式主要是尾拍稳定形式,即在运动过程中伴随着尾拍现象的发生[3],如图1 所示。射弹水下航行时,由于初始扰动的存在,射弹后体发生上下摆动而与超空泡壁面发生反复触碰,形成动态稳定[4]。
图1 水下超空泡射弹尾拍运动示意图Fig.1 Schematic diagram of underwater supercavitating projectile tail-slapping
水下超空泡射弹的尾拍运动对其弹道稳定性和作战效能具有重要影响,引起了国内外学者的广泛关注。May[5]早在1975 年就阐述了尾拍现象形成的原因,分析了尾拍对稳定航行器弹道产生的影响;Rand等[6]研究了航行器后体与空泡壁面的尾拍特性,指出考虑航行器与空泡的相对位置时超空泡航行器可能稳定在尾拍运动模式下;Savchenko[7]给出了水下超空泡航行器以尾拍运动保持稳定的速度范围;梁景奇等[8]在建立了射弹多自由度数值模型的基础上研究了弹体不同初速下的尾拍特性;陈伟善等[9]在2 种长径比下比较了平头弹、凹口弹和锥头弹的尾拍运动,并分析了其水动力特性;赵成功等[10]研究了不同质心位置的超空泡射弹尾拍运动及超空泡形态,对比分析了不同质心位置的超空泡射弹尾拍力的变化特征,张浩[11]仿真分析了在不同锥头角、平头面积弹丸、长径比、质心位置及空化槽形状等条件下射弹入水的弹道稳定性和减阻特性;Zhao 等[12]建立了一种基于数学推论的向前运动和尾拍动力学模型,对比分析不同角速度下尾拍次数和反射角的变化特征;蔡涛等[13]研究分析了空化槽的形状与尺寸对弹丸水下运动的弹道稳定性和减阻特性的影响;魏平等[14]研究了超空化射弹的稳定性模型及当弹丸尾部接触腔壁时的尾拍频率和尾拍力。综合来看,国内外学者针对弹体材料属性对尾拍运动的影响研究较少。
文中采用计算流体力学(computational fluid dynamics,CFD)软件平台,结合动网格移动计算域技术,在考虑液体压缩性的基础上开展数值仿真,在相同出口动能的条件下,研究铝合金、结构钢及钨合金3 种不同材料密度对射弹尾拍特性的影响,获得材料属性对超空泡形态及力学特性的变化规律。
水下超空泡射弹的外流场属于汽-液两相流问题,涉及相间质量传递。针对该问题,文中采用流体体积函数(volume of fluid,VOF)多相流模型模拟相界面运动。其基本控制方程包括连续性方程、动量方程和能量方程,考虑到超空泡射弹运动过程中温度变化不明显[15],可忽略能量方程。
1)连续方程
水和水蒸气的两相流动连续方程为
式中:αi、ρi和vi分别表示流体微元中第i相的体积分数、密度和速度,各相体积分数之和为1;表示从第j相到第i相的质量传递率;表示从第i相到第j相的质量传递率;n表示总相数,对于水下超空泡射弹n=2。
2)动量方程
VOF 多相流模型在整个流域内仅求解一个动量方程,并将所得速度场赋予各相,动量方程为
式中:p和v分别表示控制体的静压和速度;ρm和µm分别表示混合相的密度和动力粘度,取各相的密度和动力粘度按体积分数加权平均,即
选用工程适用性较强的基于雷诺平均纳维-斯托克斯(Reynolds equation Navier-Stokes,RANS)方程的Realizablek-ε湍流模型,该模型考虑涡流粘度的影响[16],适用于时均应变率特别大的流动问题的求解[17],具有较高的仿真精度和数值稳定性。湍流方程为
式中:C1=为经验常数;µt=ρCuk2/ε为湍流粘性系数;S=为时均应变率;Gk为速度梯度引起的湍流动能;Gb为浮力引起的湍动能;YM为可压缩湍流脉动膨胀对耗散率的影响,不可压缩流动时该项为 0;C2=1.9,C1ε=1.44为模型常数;k、ε的湍流普朗特数 σk、σε分别为1.0和1.2。
由于射弹运动的超空化流场在近壁面区域湍流未得到充分发展,k-ε湍流模型无法精确模拟,因此研究引入壁面函数对近壁面区域进行处理,以捕捉湍流边界层到过渡层区域的特性,并配合Realizablek-ε模型使用,以提高对于模拟复杂流动的适用性和计算精度[15]。壁面函数在y*>11 的网格区域满足平均速度壁面法则
式中:κ=0.42 为冯卡门常数;E=9.81 为实验系数;Uq为Q点平均流体速度;kq为Q点的湍流动能;yq为Q点到壁面距离;µ为动力粘度系数;y*为标准壁面函数在该区域的无量纲距离,y*在30~60 的区域内平均速度符合对数分布率。
Scalable 壁面函数是在标准壁面函数的基础上改进的壁面处理方法,其在y*≥11时与标准壁面函数完全一致,但为了避免在y*<11时出现数值恶化,Scalable 壁面函数引入限制量来强制使用对数分布率[18],其网格分布律可表示为
水下射弹超空化流场空化数极小且空泡较为稳定,为尽量提高计算效率和数值稳定性,文中研究采用Schnerr-Sauer 空化模型,该模型将汽相体积分数和单位体积流体含有的空泡数量联系起来,计算效率高,数值稳定性高[19]。其相间质量传递表达式为
式中,n为单位体积内的空泡数。
Tait 方程是通过采用非线性回归的方法,对能够反应p-v-T三者关系的试验数据进行拟合,而得到的液体状态方程,广泛用于可压缩液体的物性描述[15]。由于在建模时没有考虑射弹的绕流引起的温度变化,简化后的Tait 液体状态方程可描述为
式中:a和b是描述压力与密度关系的系数,基于体积模量与压力线性相关假设,其值可根据参考状态的压力、密度和模量求得;p0为参考压力;ρ0为参考压力下的液体密度;K为当前压力下的液体体积弹性模量;K0为参考压力下的液体体积弹性模量;n为密度指数;p为当前压力。
声速是介质中微弱压力扰动的传播速度,根据定义可由以下公式计算
式中:c为水中声速;ρ为当前压力下的液体密度。
针对以上控制方程组,采用有限体积法进行离散,基于分离式解法和压力耦合方程组的半隐式解法(Simple)求解流场动量方程、质量守恒方程等,其稳定性较好,可以大大提高计算速度,同时为了确保数值稳定性和结果准确性,压力求解采用Body Force Weighted 离散格式,通过假设压力和体积力之间差异的标准梯度是常数来计算面压力。对于动量,湍流方程采用2 阶迎风格式,稳定性好且收敛较快,体积分数采用改进的高分辨率界面捕捉(modified HRIC),用于生成清晰的空泡界面[20]。同时采用动网格移动计算域技术,将整个区域设置为动区域,通过用户自定义函数[21](user defined function,UDF)定义射弹的3 自由度运动,包括x、y方向的平动以及z方向的转动。
文中研究的超空泡射弹外形如图2 所示,由空化器、锥段和柱段组成。主要外形参数见表1,其中射弹总长为L;空化器直径为Dn;锥段长度La;柱段长度为Lc;圆柱段直径为Dc;质心距空化器距离为LG。
图2 超空泡射弹模型Fig.2 Supercavitating projectile model
表1 超空泡射弹主要外形参数Table 1 The main shape parameters of the supercavitating projectile
计算模型采用结构化网格,为提高计算精度,对射弹周围及柱段进行网格加密,添加流动边界层;同时为节省计算成本,提高计算效率,对远离射弹的远场网格进行稀疏处理,最终获得的网格数量为125 万,如图3 所示。
图3 弹身附近网格Fig.3 Local mesh near the body of the projectile
在流体域选取时采用有限流域来模拟无限远空间,并保证流域边界不受影响。选取圆柱形计算域,计算域布局及边界条件设置如图4 所示。
图4 计算域和边界条件设置Fig.4 Computational domain and boundary condition setting
为避免“空泡阻塞”效应,计算域直径选取50 倍射弹柱段直径,入口边界距离空化器4 倍弹长,出口边界距离航行器尾部6 倍弹长。研究中认为远场流体是静止的,射弹按照既定的规律运动,计算域四周的边界条件主要设定静压,根据射弹航行深度改变静压值;航行器表面的边界条件设置为壁面,并且壁面与临界网格相对静止,入口边界条件为压力入口,出口边界条件为压力出口。
在超空泡射弹运动过程中,其运动特性取决于流体动力特性,流体动力取决于弹体与空泡的相对位置关系,因此空泡形态计算准确性是运动模拟准确性的前提。
Hurbes[22]对锥台模型在水深4 m、速度970 m/s的条件下进行了试验,得出了高速超空化流场的空泡形态。Hurbes 试验所用模型的几何特征如图5 所示。基于Hurbes 试验所采用的模型,在同样工况下开展数值计算,将所得空泡外形和Hurbes试验结果进行对比,如图6 所示。图中分别以空化器直径Dn、射弹长度L0为参考值对空泡外形数据进行无量纲化。
图5 Hurbes 试验模型几何特征Fig.5 Geometric characteristics of the Hurbes test model
图6 970 m/s 空泡外形的数值仿真与试验数据的对比Fig.6 Comparison of numerical simulation and experimental data of 970 m/s bubble profile
图6 显示数值计算与Hurbes 试验所得超空泡外形吻合,数值仿真所得空泡轮廓略大于试验结果,如表2 所示相对偏差不超过5%,表明文中数值模型合理可行。图中,Lb为射弹全长,x为轴向位置,r为空化器径向位置。横坐标是以Lb为参考值对x进行无量纲化后获得的数据,纵坐标是以Dn为参考值对r进行无量纲化后获得的数据。
表2 数值仿真与试验数据点偏差值表Table 2 Table of the deviations of numerical simulation and experimental data points
针对图3 所示射弹网格模型,在保证边线节点分布率不变的前提下,仅改变节点数量,建立3 种不同精度的网格模型,实现网格尺寸按照倍的规律变化,网格数量分别约为95 万、125 万和160 万。基于相同的航行工况、边界设置和计算方法,分别使用这3 种不同的网格开展数值仿真,分析网格单元总量对计算结果的影响。
图7 为不同精度网格模型射弹阻力系数随时间t的变化对比。由图7 可以看出,3 种网格模型所得阻力系数变化规律结果几乎一致,网格数量在95 万时,在0.1 s 内的尾拍运动后期阻力系数有微小波动,计算精度较小;网格数量在160 万时,与网格数量125 万的数值结果一致性较好。综合考虑计算精度与经济性,文中选取适中的网格数量,认为125 万网格模型满足网格无关性要求,便于开展后续研究。
图7 不同精度网格模型射弹阻力系数对比图Fig.7 Comparison of drag coefficient of different precision grid model
1)材料属性
分别选用铝合金、结构钢及钨合金作为超空泡射弹的材料,射弹参数见表3。
表3 各材料密度射弹参数表Table 3 Table of the parameters of the projectiles under different material densities
2)尾拍过程中的空泡形态
文中研究的超空泡射弹主要利用自身动能摧毁目标,在出口动能一定的情况下,折算铝合金、结构钢和钨合金3 种射弹的初速分别为1 200、700 和450 m/s。数值仿真初始条件设定为:水深1 m,角速度5 rad/s,攻角0°。以钨合金射弹为例,提取在1 个弹道周期内的流场状态,图8 为该弹道周期内射弹与超空泡之间相对位置的演变过程,t0为初始时刻。
图8 钨合金射弹尾拍动态过程Fig.8 Motion process of tungsten alloy projectile tail-slapping
钨合金射弹在受到初始扰动5 rad/s 后,弹体首先绕质心逆时针转动,弹体柱锥段在约t=t0+12 ms时刻碰撞空泡下壁面且穿刺最深,受到尾拍作用力后射弹的扰动角速度逐渐衰减至0 rad/s;之后,弹尾在尾拍力作用下,弹回空泡内,在射弹转回至初始位置后,在惯性作用下继续绕质心转动直至撞击上侧空泡壁面,如此反复,射弹在航行过程中处于动态平衡。对于铝合金射弹,由于自身质量较小,在做尾拍运动时速度衰减更快,空泡形态由开始的开口泡逐渐变为闭口泡,图9 为尾拍运动后期铝合金射弹的闭口泡形态图。
图9 铝合金射弹尾拍动态过程Fig.9 Motion process of aluminum alloy projectile tailslapping
采用变参数法研究不同材料密度时超空泡射弹尾拍的流体动力特性及运动特性,设定水深1 m,忽略重力、海流和波浪的影响,模拟定深绕质心旋转的纵平面尾拍运动。以地面系为参考坐标系,获得铝合金、结构钢及钨合金 3 种不同材料射弹受扰动后的攻角、俯仰角速度、尾拍周期、垂直方向速度和运动位移变化曲线。
图10 为3 种材料密度的射弹攻角变化曲线,可以看出,材料密度越小,射弹在初始段的尾拍振荡频率越高,随着航速的衰减,尾拍周期缓慢增大;相同弹型下,3 种射弹的攻角变化范围均稳定在4°以内,在-2°~2°之间周期性变化,攻角变化幅值随时间略有减小,铝合金射弹由于速度衰减较快,攻角随时间变化其幅值减小较为明显。
图10 不同材料密度射弹攻角随时间变化曲线Fig.10 Attack angle under different material densities versus time
图11 为俯仰角速度变化曲线,材料密度改变对射弹尾拍运动俯仰角速度的影响主要体现在初期振荡周期和幅值变化上。材料密度越小的射弹形成尾拍的时间越早,弹体在空泡内的运动时间越短,俯仰角速度变化率就越大,尾拍频率越高,约0.2 s 后,振荡趋于稳定;材料密度越大的射弹初始俯仰角速度幅值越小,3 种射弹俯仰角速度幅值均收敛在±5 rad/s 以内。
图11 不同材料密度下射弹俯仰角速度随时间变化曲线Fig.11 Angular velocities of pitch angle under different material densities versus time
提取射弹运动速度与尾拍运动周期,形成对应曲线如图12 所示。由图可见,在相同速度时,材料密度越大的射弹,尾拍运动周期越大;同一材质的射弹,随着速度的衰减,尾拍周期呈现增大的趋势,其中钨合金射弹的尾拍周期最大。
图12 不同材料密度下射弹速度与尾拍运动周期对应关系曲线Fig.12 The curves of correspondence between projectile velocity and tail-slapping period under different material densities
如图13 为3 种射弹动能随时间的变化曲线。由图可见,曲线均呈下凹形,密度越大的射弹在相同航行时间内动能衰减越小,即材料密度最大的钨合金射弹在相同时间内动能衰减最小,材料密度最小的铝合金射弹在0.1 s 内动能迅速衰减,随后由于尾拍周期的增大,阻力较初期尾拍减小,速度衰减减慢,动能衰减趋于平缓。
图13 射弹动能随时间变化曲线Fig.13 The variation curves of kinetic energy of projectile versus time
图14 和图15 分别为射弹垂直方向速度和质心位移变化曲线。由图可见,3 种射弹垂直方向速度均呈现周期性振荡衰减趋势,质心纵平面的运动轨迹呈现一定波长的稳定振荡特性。铝合金射弹尾拍初期,尾拍周期较短(见图10),沾湿面积较大,较大的升力使其产生y轴偏移,并积累了较高的垂向速度,随后射弹速度衰减,尾拍周期增大,升力减小,垂向速度趋于平稳;结构钢射弹垂向速度在0 附近振荡;钨合金射弹的垂向速度总体为正,当弹体完全被空泡包裹时,几乎不受侧向力,因此垂向速度几乎不变,当弹体下表面沾湿时,升力方向向上,垂向速度增加,反之则垂向速度减小。初始段,铝合金射弹的质心位移最小,结构钢次之,钨合金最大。在相同出口动能和初始扰动下,不同材料密度射弹的尾拍运动质心运动轨迹均近似呈直线,当质心的纵向偏移量相对于射弹的航程可以忽略不计时,射弹具有良好的弹道特性。所以结合各材料射弹的水平及垂向速度变化规律可知,以结构钢为材料的射弹有良好的弹道特性,此时射弹航程可达200 m,纵向偏移量为0.15 m。
图14 垂直方向速度随时间变化曲线Fig.14 The variation curves of vertical velocity versus time
图15 不同材料密度下质心位移变化曲线Fig.15 The variation curves of centroid trajectories under different material densities
根据仿真计算结果,提取水下射弹的全弹道流体动力特性,将阻力、升力及俯仰力矩进行无量纲化处理,得
式中:CF为升力和阻力系数;CM为俯仰力矩系数;F为射弹阻力;ρ为水的密度;vl为射弹速度;S为圆盘空化器的面积;M为射弹俯仰力矩;L为射弹全长。
3 种射弹阻力系数(见表4)、升力系数以及俯仰力矩系数均呈现出显著的周期性变化规律,如图16~18 所示。在运动初期,振荡频率随着射弹材料密度的增大而增大,振荡幅值随射弹材料密度的增大而减小,在运动后期,不同材料密度的射弹流体动力曲线的振荡频率值和振荡幅值均趋于稳定。在形成更稳定的尾拍运动后,超空泡射弹的流体动力特性取决于头部空化器和周期性尾拍碰撞空泡壁面所形成的稳定沾湿区域。不同材料密度的射弹俯仰力矩系数均随时间的变化保持稳定的周期性变化,其中铝合金材料射弹随着时间的变化其俯仰力矩系数略有增大,由于随着速度的衰减,射弹在穿刺空泡之后的滑水力数值变化不大,但是滑水力作用点逐渐后移,使得滑水力的力矩增大,因此,以实时速度和射弹全长为参考值计算得到的俯仰力矩系数的幅值逐渐增大。根据图18 曲线特征,当射弹在尾拍过程中尾部产生沾湿时,俯仰力矩系数发生一定波动,材料密度越大的射弹波动幅值越大。
图16 不同材料密度下阻力系数变化曲线Fig.16 The variation curves of drag coefficient under different material densities
图17 不同材料密度下升力系数变化曲线Fig.17 The variation curves of lift coefficient under different material densities
图18 不同材料密度下俯仰力矩系数变化曲线Fig.18 The variation curves of pitching moment coefficient under different material densities
文中针对超空泡射弹在不同材料密度影响下形成尾拍运动的过程进行数值仿真,分析了不同材料密度下射弹的弹道特性和流体动力特性的变化规律,得出如下结论。
1)不同材料密度的超空泡射弹的尾拍运动攻角、俯仰角速度、垂直方向速度呈周期性变化规律。其攻角均在-2°~2°之间周期性规律变化,随着速度的衰减,攻角的幅值逐渐减小。
2)材料密度越大的射弹,尾拍周期越长。钨合金材料射弹的尾拍周期最长,尾拍初始阶段一次尾拍周期为0.5 s,随着速度的衰减,尾拍周期逐渐增大,且增大的幅度最小,仅增大0.15 s。
3)材料密度越大的射弹在相同时间内速度与动能衰减越慢,余速越大,俯仰角速度振荡幅值越小,最终3 种材料密度的射弹俯仰角速度均收敛在±5 rad/s 以内。对于垂向速度,铝合金材料的射弹初始尾拍阶段累积较高的垂向速度,随着尾拍周期的增大,垂向速度逐渐趋于平稳,结合其质心位移变化规律,不同材料密度射弹的尾拍运动质心运动轨迹均近似呈直线,结构钢材料射弹纵向偏移量最小,具有良好的弹道特性。
4)3 种射弹的阻力系数、升力系数以及俯仰力矩系数均呈现出显著的周期性变化规律。铝合金材料射弹由于材料密度较小,初始阻力值最大,随着速度的衰减,阻力减小的幅度最大,随后逐渐趋于稳定,升力随着时间的变化先增大后减小,俯仰力矩由于滑水力作用点的后移逐渐增大。
未来工作将选取更多合适的材料密度,对超空泡射弹的弹道特性开展研究,考虑材料密度对超空泡射弹阻尼力特性的影响,为优选水下超空泡射弹材质提供参考。