剪切来流下超空泡射弹空化与水动力特性的数值研究

2023-04-06周标军赵子杰

弹道学报 2023年1期

袁馨,周标军,戴琪,赵子杰,张辉

(南京理工大学瞬态物理国家重点实验室,江苏南京 210094)

射弹在水下高速航行时,弹身附近的水压降低并达到饱和蒸汽压,周围液体发生相变而产生空泡[1-2]。当射弹速度继续增加时,空泡延伸并完全包裹弹身,产生超空泡,此时水蒸气代替水与射弹接触,会大大降低射弹受到的阻力。采用超空泡减阻的水下射弹具有高效的突防能力和较远的水下射程,因而具有广阔的应用前景[3]。

目前,很多学者对水下超空泡射弹的空化形态及水动力特性进行了数值模拟研究。张艳洲[4]对典型水下高速射弹模型进行数值模拟,研究射弹的超空泡生成和发展过程,发现射弹的头部最先发生空化,随后肩部和尾部产生空泡并逐渐包裹弹体,最后形成一个完整且对称的空泡;王瑞等[5]建立了带尾翼射弹的水下高速航行计算模型,结果表明射弹的阻力系数随着速度的增加而增大;LI等[6]构建了无附体射弹的跨声速超空化流场数值模型,研究了马赫数为0.42～1.68速度范围内射弹周围的流场压力分布和阻力特性,结果表明随着马赫数的增加,射弹头部压力急剧增加,导致射弹阻力系数增加且弹头附近空泡径向尺寸显著减小;熊天红等[7]通过数值模拟,研究了不同射弹模型在高速运动下的超空泡形态以及阻力系数,发现随着空化数的减小,空泡直径与长度均增大,阻力系数随之减小。

然而,实际海洋环境中受海岸面、目标物体表面等受限壁面的影响,射弹绕流速度随着与壁面距离的变化呈现出剪切分布的特性[8-12]。近年来,流动剪切性的作用受到了学者们的广泛关注。BRANDAO等[8]通过大涡模拟研究了后台阶剪切层流动中的空化特性,并取得了与实验一致的结果;张辉等[9]对剪切来流作用下圆柱绕流的水动力特性进行了实验和数值模拟研究,发现在剪切来流的绕流流场中,圆柱所受的升力偏向流速较小的一侧;WIMSHURST等[10]模拟了由斜坡引起的剪切来流作用下水轮机的水力性能,结果表明剪切流动会导致非对称性载荷;杨福昌等[11]对不同剪切率下的弹性管绕流进行了数值分析,总结了弹性管升力系数随剪切率的变化规律。

以往水下超空泡射弹的研究大都局限于均匀来流的情形,而在实际的航行过程中,射弹还会受到剪切来流的作用,但是相关研究到目前为止还是空白。因此,分析剪切来流对超空泡射弹的空化与水动力的影响规律具有十分重要的意义。本文针对典型超空泡射弹模型,采用Mixture多相流模型、Schnerr and Sauer空化模型和Realizablek-ε湍流模型,对剪切来流条件下的高速射弹绕流场进行数值模拟研究,并与均匀来流进行对比分析,研究剪切流作用下的空泡形态以及绕射弹流场动力特性,分析射弹阻力系数、升力系数随剪切率的变化规律,揭示来流剪切率对超空化的作用机理,以期为超空泡射弹构型设计提供理论依据。

1 数值方法

1.1 控制方程

本文采用不可压缩的Mixture多相流模型,该模型计算效率高,将水与水蒸气两相混合流作为密度可变的单一流体来处理,且水下运动温度变化较小,可忽略能量方程,故对混合物的连续性方程和动量方程进行模拟求解：

(1)

(2)

由于超空泡射弹在水下高速航行,因此湍流模型采用对高雷诺数复杂流场模拟效果较好的Realizablek-ε湍流模型。非平衡壁面函数适用于处理与压力梯度相关的复杂流动问题,所以采用该壁面函数处理近壁面区域的流动。气液两相之间的相变空化采用Schnerr and Sauer空化模型来计算[13],该模型数值稳定性强,计算效率高。对于控制方程的离散,采用有限体积法求解,二阶迎风差分格式计算空间导数,时间离散采用二阶隐式格式,压力与速度之间的耦合求解采用Coupled算法。气液两相体积各占50%(体积分数等于0.5)的等值面视为气液界面。

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

如图1所示,本文所选用的弹型为一个由截锥体与圆柱体组成的回转体,弹长L=170 mm,弹径D=12.65 mm,头部空化器为圆盘形,空化器直径Dn=4 mm。

图1 射弹模型Fig.1 Projectile model

为了避免“边界效应”[14-15]对射弹超空泡的影响,计算流域选取直径为20L、长度30L的圆柱形流域,如图2所示。以空化器中心为坐标系原点,射弹轴向为笛卡尔坐标系下x方向的正向,径向为y方向的正向,右手定则确定z方向的正向。流域入口距离弹头10L,为速度入口边界条件,加载剪切来流入口速度函数u(y),x轴上方为高速侧,下方为低速侧;流域出口距离弹尾19L,边界条件为压力出口;圆柱面远场边界设为压力出口。饱和蒸汽压强设为3 540 Pa。SAVCHENKO等[16]通过实验总结出Fr>30时可以忽略重力影响,We>2 000时可以忽略液体表面张力的影响,本文研究的Fr=232,We=1.4×106,因此,二者影响均可忽略。黄闯等[17-18]研究发现亚音速范围内的液体可压缩性对弹丸受力及空泡形态影响较小,故文中计算采用不可压缩流体。

图2 计算域构建及边界设置Fig.2 Establishment of computational domain and setup of boundary conditions

考虑到计算域大小限制,设置剪切来流入口速度为

(3)

式中：u0=(u1+u2)/2为剪切流平均(中心)速度;u1、u2分别为剪切流最大、最小速度;y∈[0,10L]为高速侧,y∈[-10L,0]为低速侧,剪切区域取为y∈[-0.1 m,0.1 m];Δδ为剪切区域的一半,即Δδ=0.1 m。弹体位置处的来流剪切率计算为

(4)

本文探究来流剪切率对超空泡射弹的影响,设置剪切流平均速度u0=600 m/s,模拟了剪切率K分别为0,2 000 s-1,4 000 s-1,6 000 s-1,6 500 s-1,7 000 s-1和7 500 s-1的7个算例。

采用ICEM软件的O-Block技术划分结构化网格,总体质量在0.6以上,过渡均匀。为精确捕捉空泡形态和弹丸受力,根据y+的计算结果优化近壁面处的边界层网格,并通过经验公式计算得到空泡范围,在该范围内加密网格。网格总数为167万,弹体附近的网格如图3所示。

图3 网格划分设置Fig.3 Setup of mesh generation

1.3 水动力系数定义

射弹阻力系数定义为CD=FD/(0.5ρv2S),FD为总阻力,包括压差阻力FDP和黏性阻力FDV[19],方向沿流向;升力系数定义为CL=FL/(0.5ρv2S),其中FL为总升力,可分解为压差升力FLP和黏性升力FLV[20],方向沿法向;ρ为流体密度,v为来流平均速度;S为射弹最大横截面积。射弹各位置处的阻力与升力也均可分解为垂直于弹体表面的压力与平行于弹体表面的黏性力。

为了探究剪切来流作用下弹体表面的阻力分布,本文对压差阻力FDP和黏性阻力FDV的曲面积分：

(5)

(6)

式中：px为单位面积的流向压力,τx为单位面积的流向剪切应力。

进行分解,计算得到不同流向位置处弹体表面的流向压力和黏性力,而流向合力为流向压力与黏性力之和：

(7)

(8)

fD,x=fDP,x+fDV,x

(9)

式中：Cx为该流向位置处弹体截面的周长,循环计算x=0～0.17 m中所有流向位置处的数据,最终可得到流向力沿流向的分布。

同样,为了探究剪切来流作用下弹体表面的升力分布,本文对压差升力FLP和黏性升力FLV的曲面积分：

(10)

(11)

式中：py为单位面积的法向压力,τy为单位面积的法向剪切应力。

进行分解,计算得到不同流向位置处弹体表面的法向压力和黏性力,而法向合力为法向压力与黏性力之和：

(12)

(13)

fL,x=fLP,x+fLV,x

(14)

式中：Cx为该流向位置处弹体截面的周长,循环计算x=0～0.17 m中所有流向位置处的数据,最终可得到法向力沿流向的分布。

1.4 数值验证

1.4.1 数值方法验证

众多国内外学者提出并修正了空泡轮廓线的经验公式,其中LOGVINOVICH[21]提出的经验公式应用广泛。依据LOGVINOVICH经验公式,超空泡最大截面半径Rc及最大空泡长度Lc分别为

式中：Rn为空化器半径,Cx0为空化数σ=0时的阻力系数,圆盘空化器的Cx0给定为0.83;系数ks与空化数大小有关,空化数越小,ks越大,可以由下式求得：

依据超空泡独立膨胀理论,得出超空泡的轮廓公式：

式中：x为超空泡轮廓距离空化器截面的水平位移,R(x)为空泡半径。

为验证文中所采用的数值方法,将均匀来流300 m/s条件下的仿真结果与LOGVINOVICH经验公式对比。该工况下液体的来流压力p0=111 307 Pa,饱和蒸汽压p1=3 540 Pa,因此空化数σ=(p0-p1)/(0.5ρv2)=0.002 39,满足LOGVINOVICH经验公式的空化数适用范围(0,0.25)。图4为仿真数据和经验公式对比图。由图4可知,液相体积分数为50%的气液交界面轮廓仿真结果与经验公式吻合较好,误差在3%左右。

图4 仿真数据和经验公式对比Fig.4 Comparison between simulation results and empirical formulas

进一步,本文应用HRUBES[22]实验中的弹型及工况,在此数值方法下进行仿真,结果对比如图5所示,可以看出仿真空泡与HRUBES实验中的空泡轮廓基本一致,空泡半径的相对偏差低于2%。综上说明,本文所采用的数值方法对于超空泡射弹具有较高的模拟精度。

图5 仿真与实验结果对比Fig.5 Comparison of simulation and experimental results

1.4.2 网格无关性验证

在不改变网格总体分布规律的基础上,通过调整线节点数,分别建立数量为100万、167万、340万的3种不同精度的网格,在速度为300 m/s的工况下进行模拟,所得空泡形态的网格无关性检验如图6所示。由图6可知,167万网格和340万网格下计算的空泡轮廓一致,而100万网格下计算的空泡轮廓较小。在保证足够计算精度的前提下,实现效率最大化,本文选取数量为167万的网格进行计算。

图6 不同网格密度下的空泡形态Fig.6 Shapes of supercavities under different mesh density

2 结果与讨论

2.1 超空泡形态分析

表1为不同来流条件下空泡发展过程。如表1所示,选取剪切率K分别为0,4 000 s-1和7 500 s-1的3个代表性算例,研究不同来流条件下的空泡形态,白色部分为射弹模型,深蓝色部分为空泡。按空泡形态随时间演化的过程分为6部分：弹头初生阶段;弹肩、弹尾初生阶段;空泡包裹弹肩阶段;空泡包裹弹体阶段;空泡继续发展阶段和空泡稳定阶段。均匀来流时,起初射弹头部产生空泡且空泡尺寸最大,肩部和尾部产生的空泡直径接近于弹丸的直径,头部产生的空泡继续向后发展形成一个完整的空泡,包裹住整个弹体。剪切来流下,起初射弹头部产生空泡,而随后肩部和尾部低速侧先于高速侧发生空化,并且低速侧空泡演化更快,体积更大。头部、肩部和尾部生成的空泡向后方同时发展,肩部产生的空泡迅速延伸至弹尾,随后头部产生的空泡继续向后发展,包裹住弹体,也形成一个完整的空泡。在空泡继续发展过程中,高速侧的空泡较快发展,略微缩小与低速侧空泡的差距。稳定后,包裹弹体的超空泡表现出明显的不对称性,空泡向低速侧偏斜,高速侧的空泡半径较小。此外,当剪切率K较大时,稳定后的空泡半径明显减小,弹肩高速侧存在一定程度的沾湿现象,弹体与水直接接触。

表1 不同来流条件空泡发展过程Table 1 Development of cavity under different inflow conditions

图7展现了不同剪切率下的空泡轮廓,空泡轮廓在不同剪切率下存在明显差异。均匀来流中空泡轮廓上下完全对称,而剪切来流下空泡高低速侧呈现出不对称性,随着剪切率增大,高速侧空泡半径减小,低速侧空泡半径无明显规律。当剪切率进一步增大时,弹肩高速侧发生沾湿。

图7 不同剪切率条件下的空泡轮廓(K=0～7 500 s-1)Fig.7 Cavity shapes for different shear rates (K=0-7 500 s-1)

空泡的不对称性可以通过空泡中心线位置δc来定量表示,如图8所示。剪切率较小弹体未沾湿时,随着剪切率增大,空泡的不对称性不断增大,空泡中心线逐渐向低速侧偏斜;而当剪切率较大弹体有沾湿时,空泡的中心在沾湿位置之后逐渐回升,空泡的不对称性明显减弱;随着剪切率进一步增加,沾湿增强,上空泡半径减小,下空泡半径增大,导致空泡的不对称性增强。

图8 不同剪切率条件下的空泡中心线位置Fig.8 Locations of cavity centerlines for different shear rates

2.2 射弹阻力特性分析

基于上述仿真算例,分析来流条件对射弹阻力系数的影响。图9给出了空泡稳定后不同来流剪切率的总阻力系数值、压差阻力系数值和黏性阻力系数值。随着剪切率的增加,阻力系数逐渐增大,而当剪切率K>6 000 s-1,射弹发生沾湿时,阻力系数更加显著地增大。剪切率较小时,弹体被空泡中的水蒸气完全包裹,黏性阻力很小,可以忽略不计,总阻力受压差阻力主导;而在剪切率较大时,射弹周围发生沾湿,弹体局部与水直接接触,导致黏性和压差阻力均显著增大,使得总阻力大幅增加。

图9 空泡稳定后剪切来流作用下的阻力系数Fig.9 Drag coefficients after cavitation stabilization in shear flow

为了进一步探究不同来流条件下弹体表面的阻力分布,采用1.3节的方法计算得到均匀来流中流向合力、压力和黏性力沿流向分布的曲线,如图10所示。因为弹头迎风面直接与水接触,弹头处受到的压力为一个极大值,而弹身区域被空泡包裹,空泡内部各点的相对压力大小均为饱和蒸汽压,故弹肩前的圆台部分压力的流向分量恒定为一固定值,弹肩后形状对称的圆柱部分压力为0,弹尾负压区的压力显著小于弹头。流向黏性力在弹头处为0,随后一直保持比较小的数值,经过弹肩后缓慢降低,在弹尾处重新为0。因此,流向合力主要分布在弹头位置处,弹肩后合力很小,弹尾处合力与弹头处方向相反,大小明显小于弹头,使得弹体受到的阻力主要由弹头和弹尾的压力差构成。

图10 均匀来流中弹体表面流向力的流向演化Fig.10 Streamwise evolution of the area-weighted average streamwise force in uniform flow

选取均匀来流、剪切率较小和剪切率较大的3个代表性算例,研究弹体表面流向力沿着流向方向的分布曲线(局部放大图对应主图中蓝色框线区域),如图11所示。剪切率较小时,弹体被空泡完全包裹,剪切来流下流向合力、压力和黏性力的分布趋势与均匀来流基本相同。弹头处压力随着剪切率的增加而增大,而弹尾上的压力与剪切率无关,如图12(a)所示,因此压差阻力显著增大。弹体圆台面压力流向分量为一固定负值,受剪切率影响较小,而黏性力随着剪切率的增加而略微增大;圆柱面流向压力为0,黏性力也随着剪切率的增加而略微增大。这是因为单位面积的剪切应力τ=μdu/dy,空泡内黏度μ不变,随着剪切率的增加,弹体高速侧流向速度梯度显著增大,低速侧速度梯度略微减小,如图12(b)所示,使得流向黏性力增大。剪切率较大时,弹体的压差阻力进一步增大,并且弹肩处发生沾湿,弹肩附近流向压力逐渐增大至0,弹体周围环境也由蒸汽变成了水,因而黏度增大,黏性力也随之显著增大,使得流向合力增大,既而导致阻力增大。

图11 不同剪切率的弹体表面流向力的流向演化Fig.11 Streamwise evolution of the area-weighted average streamwise force for different shear rates

图12 不同剪切率的展向切面流场云图Fig.12 Contours on the spanwise slice for different shear rates

2.3 射弹升力特性分析

基于上述仿真算例,分析来流条件对射弹升力系数的影响。图13给出了空泡稳定后不同来流剪切率下的升力系数值。均匀来流中,弹体被空泡完全包裹,且空泡对称,因而升力系数为0。然而,剪切来流作用下,升力系数变成负值,表明射弹受到朝向低速侧的法向力作用。剪切率较小时,随着剪切率增大,升力系数幅值略微增加,而当K>6 000 s-1时,弹体肩部沾湿,升力系数的幅值显著增加。

图13 空泡稳定后剪切来流作用下的升力系数Fig.13 Lift coefficients after cavitation stabilization in shear flow

进一步,图14比较了剪切率为0,4 000 s-1和7 500 s-1的3个代表性工况下法向合力、压力和黏性力的流向分布。均匀来流中,因为弹体上下侧的流速一致,空泡完全对称,所以法向的压力、黏性力和合力都恒为0。剪切来流下,当剪切率较小时,弹体仍然被空泡完全包裹,弹头、弹尾处的压力沿着流向方向,弹身上的压力高低速两侧相抵,因而压力的法向分量也都为0。由于弹体高速侧绕流更快,更多液体被排挤到低速侧,弹头处的法向黏性力为负值,而后先波动再沿着流向方向缓慢增加,并且剪切率越大,法向黏性力越大。

图14 不同剪切率的弹体表面法向力的流向演化Fig.14 Streamwise evolution of the area-weighted average normal force for different shear rates

图15展示了展向切面中弹体周围的流线分布,图16给出了流向切面中弹尾前端的速度矢量。均匀来流中,弹体周围的流线和速度矢量在高低速侧完全对称,并且在空化相变作用下,弹体周围产生的蒸汽向外均匀扩散,弹尾区域产生一对旋转方向相反的对称规则涡旋,上下侧卷吸的流体相等。而在剪切来流下,高速侧绕流作用更强,弹体周围的流线向低速侧偏斜,并且弹尾低速侧的涡旋作用更强,卷吸更多的流体从高速侧向低速侧运动[23-24],导致空泡不对称,并向低速侧偏斜,也使得弹体受到朝向低速侧的法向黏性力作用,因而总的升力也朝向低速侧方向。图17展示了剪切率为7 500 s-1的弹体表面法向合力、压力和黏性力沿着流向方向的分布。由于弹肩处发生沾湿,高速侧与水直接接触,黏度显著增大,黏性力因而明显增加,并显著大于低速侧,使得法向黏性力朝向高速侧方向。但是沾湿处的压力由饱和蒸汽压变为水压,高于低速侧的蒸汽压,压力的法向分量剧烈地指向低速侧方向,并且显著大于法向黏性力,导致朝向低速侧的升力进一步增大。