高 山, 施 瑶, 潘 光
(1. 西北工业大学 航海学院,西安 710072; 2. 无人水下运载技术工信部重点实验室,西安 710072)
航行体水下发射过程中近似可分为三个阶段,即出筒阶段、水中航行阶段以及出水阶段。当航行体以较高的速度出筒时,其肩部低压区导致附着空泡生成。在水中航行阶段中,随着航行体速度减小,空化数增大,附着空泡在回射流的影响下发生脱落及断裂现象。在出水阶段,由于自由液面附近的介质密度差别极大,空泡在接触到外界大气后迅速产生相变而发生空泡溃灭,形成指向航行体壁面的射流,造成强烈的冲击载荷,严重影响到航行体出水姿态和结构安全性[1]。
在水下垂直发射过程中,由于航行体周围环境压力不断变化、表面的空化数随着航行体所处水深的变化而变化,导致相关的水洞试验[2]和数值模拟方法[3-4]难以适用水下发射特殊环境相关研究。关于水下发射过程航行体肩部附着空泡演化主要涉及大尺度空泡群生成、发展以及溃灭等复杂现象。在理论方面,目前主要集中在单空泡相关研究,但研究成果难以实际应用。王一伟等[5]从目前已有的数值模拟方法与试验测试手段,系统地总结了高速航行体水下发射水动力相关研究进展,特别是水下航行阶段空泡稳定性、出水空泡溃灭等关键问题。另外,王一伟等[6]对航行体出水附着空泡溃灭过程开展了研究,建立了空泡溃灭压力的物理模型,给出了相关参数的影响规律。陈玮琪等[7]基于势流理论、细长体理论和奇点分布法,对附着空泡、自由液面以及筒口气团的相互影响开展了理论研究,并采用试验结果验证了建立模型的合理性。施红辉等[8]采用试验测试手段研究了完全超空泡形态下航行体出水过程,建立了超空泡崩溃次数的相关模型。权晓波等[9]针对不同发射参数下空泡形态开展了试验测试,当在大攻角下,其周围空泡分布不对称加强,迎流侧和背流侧的压差范围增大。Zhang等[10]采用高速摄影技术研究了碎冰影响下航行体出水空化特性演变规律,获取了碎冰条件下航行体附着空泡以喷发状态溃灭机制。在数值模拟方面,魏英杰等[11]采用动网格技术结合混合介质RANS方程,研究了航行体垂直发射过程空化特性研究,给出了航行体阻力系数与空化数之间的关系。刘元清等[12]对海流对肩部空泡演变过程的影响进行了数值模拟研究,发现在不同海流状态下泡内的水气分布存在差异,因而溃灭过程中高压产生次序也不同。李卓越等[13]基于雷诺平均数值模拟方法和Schnerr-Sauer空化模型,建立了航行体水下垂直发射三维数值模型,分析了多种初始因素影响下壁面载荷特性。孔德才等[14]采用Singhal空化模型模拟了空泡界面对航行体头部附近流场压力的影响,并与试验结果进行了对比验证。刘涛涛等[15]基于均相流模型对垂直发射水下航行体周围的通气空化流动进行了三维数值模拟,并讨论了通气时序对通气空化流场的影响。施瑶等[16]研究了双发回转体齐射过程空化流场结构演变规律,发现齐射过程中由于流动干扰区域的存在,双体肩部空泡形态演变过程由不对称演变为对称状。
综上所述,虽然目前国内外针对水下发射非定常空化特性研究取得了一些成果。然而,关于出水过程湍流相干结构[17]与溃灭载荷内在演变机理的研究相对较少,且模拟方法主要集中在混合介质RANS。另外,由于水下发射试验测试难度高,成本高以及难以在航行体表面设置多个压力传感器,导致很难充分测量其压力脉动规律。因此本文采用改进型分离涡模型(improved delayed detached eddy simulation,IDDES)研究出水空泡形态演变与涡旋流场之间的内在联系,同时设置多个探针点来监测空化区域的脉动压力并分析其演变规律。
水下发射过程航行体带空泡出水过程中的力学环境变化剧烈,包含的物理现象和机制更为复杂。本文采用多相流模型(volume of fluid,VOF)对混合流场进行描述,基本控制方程形式如下:
连续性方程
(1)
动量方程
(2)
能量方程
(3)
水蒸气相方程:
(4)
考虑水蒸气相为可压缩理想气体,补充状态方程
p=ρgRgT
(5)
式中:t为时间;ui为速度分量;xi为坐标方向;p为压力;T为温度;gj为重力加速度分量;Ek为第k相流体总能;ρk为第k相流体密度;keff为混合物有效传热系数;Rg为气体常数。
本文采用的延迟IDDES模型基于SSTk-ω模型进行构造,引入湍流长度尺度lIDDES,对模型中湍动能耗散项进行了修正
(6)
其中,lIDDES的基本形式如下
式中:Sij为应变率张量;km和τ分别为湍流动能和湍流剪切力;CDES为模型系数,通常取为0.65;Δmesh为网格尺度;dw为计算点到壁面的距离;hmax为网格最大边长;hwn为垂直壁面方向的网格尺度。
本文采用的Schnerr-Sauer空化模型是第一个不需要经验常数的质量输运空化模型,同时考虑了气泡加速增长影响、黏性效应以及表面张力效应等,其质量源表达式为
(7)
(8)
式中:单位体积液体内的气核数n=1012/m3;pv为与空化相变相关的饱和蒸气压。另外,本文采用商用软件STAR-CCM开展计算相关计算。
航行体模型与边界条件示意图,如图1所示。其中直径D=0.04 m,长径比L/D=6,采用半球头型。航行体质心距离头部顶点的距离为0.5L,初始发射位置距离水面的深度为L;航行体以初始速度(V=5~6 m/s)和初始攻角(AOA0=0°~8°)向上作无动力自由运动,其中初始空化数范围是σ=0.198 0~0.023 1。水面上部区域是一个压力可调的空气域,以满足模拟与真实条件下空化数相似(Po=437 1 Pa),其中25℃下水的饱和蒸汽压为Pv=317 0 Pa。左侧、右侧、上部均为压力出口边界条件,其截面位置设置为对称边界。
图1 航行体模型与边界条件Fig.1 Vehicle models and boundary conditions
网格细节划分,如图2所示。背景加密区域和重叠域中基本网格单元大小为2.5%D。另外,为了保证满足IDDES湍流模型计算要求,靠近航行体第一层网格高度Y+取1,边界层共20层。重叠域网格数4.7×106,背景网格数1.63×107,总网格数2.1×107。
本文所采用的数值模型是在前期研究[18]基础上进一步考虑了自然空化的影响。在不考虑空化作用下该模型已得到了验证。因此,本文通过模拟局部自然空化流动来验证所采用的空化模型,从而验证所采用数值模型的有效性。其中,数值模拟的航行体直径为10 mm,长径比为10,与试验模型[19]保持一致;速度恒定为50 m/s,无穷远处水流场压力为378 955 Pa。另外,根数值模拟环境为了保持和试验环境一致,须满足雷诺数106和空化数0.3。
由图3可知,模拟数据和试验结果吻合度较好,其中空泡闭合区域的突增均发生在X/D=2附近;如图4所示为获取的航行体物面压力系数分布与Rouse等文献中试验数据对比图,其中X表示空泡末端距离航行体头部顶点轴向距离。图4表明压力系数演变基本一致,表明航行体在空化流场中所受的流体动力和试验结果也是近似相同的。因此,本文采用的数值模拟方法是有效的。
图3 空泡演变过程气相体积分数αv分布云图Fig.3 Cloud diagram of vapor phase αv distribution
图4 压力系数分布Fig.4 Distribution of pressure coefficients
当考虑水下发射自然空化现象时,应满足以下方程
f(t,V,μw,p∞,ρw,g,D,L)=0
(9)
式中:基本物理参数包括D为航行体长度;ρw为水密度;V为航行体发射速度。因此,可以获得式(9)的无量纲形式,如式(10)所示
(10)
式(10)进一步可写为如式(11)形式:
(11)
如图5和图6所示为攻角AOA0=0°,空化数σ=0.198下航行体出水空泡脱落及溃灭与压力云图演变。从图5中可以发现,航行体头部还未接触水面之前,回射流已到达航行体肩部空泡中间区域。如图6所示,空泡末端回射流区域出现较高的逆压梯度。当航行体继续运动过程中,其头部逐渐顶起自由液面,头部高压区消失,肩部空泡与液态水存在一个明显的交界面。然而,自由液面之上的高压区域会导致此交界面附近的水蒸气出现凝结现象,引起局部水蒸气发生液化。同时,回射流顶点附近的空泡形态自上而下发生了溃灭现象。随着航行体继续运动,回射流顶点下面空泡继续发生溃灭现象,直至大尺度空泡变成孤立的小空泡。在外界高压作用下,孤立空泡发生了强烈的坍塌溃灭现象。同时,形成高速壁面射流,从而冲击航行体壁面。
图5 AOA0=0°,σ=0.198下航行体出水空泡脱落及溃灭演变Fig.5 Evolution of the cavity shedding and collapse at AOA0=0° and σ=0.198 during vehicle exiting-water
图7所示为AOA0=0°,σ=0.231下航行体出水空泡脱落及溃灭演变细节图。与图5相比可知,随着空化数增大,空泡界面形态变得薄而短,肩部空化区域体积明显减小。分析发现,在航行体肩部空泡溃灭开始溃灭之前,泡内回射流顶点与航行体之间近似保持静止状态,这主要是因为随着回射流相对于航行体向上移动,空泡区域本身也会向下发展。在回射流和界面高压的共同作用下,空泡发生坍塌溃灭现象。然而,需要注意的是航行体尾流区域的空泡出现反复的收缩和膨胀,直到航行体完全离开自由液面。图8和图9分别是AOA0=4°和AOA0=8°下出水等值面空泡与速度矢量演变。从图中可以发现,航行体出水时刻的攻角状态对空泡形态演化具有显著的影响。在高速回射流的冲击下,背流侧空泡发生大尺度空泡界面断裂和脱落现象。随着出水攻角增大,迎流侧空泡和背流侧空泡演变行为不对称加剧,对空化区的涡旋结构演变机理和溃灭载荷将产生剧烈的影响。
图7 AOA0=0°,σ=0.231下航行体出水空泡脱落及溃灭演变Fig.7 Evolution of the cavity shedding and collapse at AOA0=0° and σ=0.231 during vehicle exiting-water
图8 AOA0=4°,σ=0.198下出水等值面空泡与速度流线演变Fig.8 Evolution of iso-surface cavity and velocity streamlines at AOA0=4° and σ=0.198 during the exiting-water process
图9 AOA0=8°,σ=0.198下出水等值面空泡与速度流线演变Fig.9 Evolution of iso-surface cavity and velocity streamlines at AOA0=8° and σ=0.198 during the exiting-water process
航行体以一定攻角出水时,由于空化区域剧烈演变,其周围流场的湍流相干结构受到显著的影响。图10给出了AOA0=4°,σ=0.198下出水速度矢量场演变。水下航行阶段,在回射流区域涡旋结构的作用下,航行体背流侧发生了小尺度空泡脱落现象。脱落空泡在远离航行体肩部附着空泡时,发生了溃灭现象,值得注意的是脱落空泡的溃灭现象与航行体肩部空泡出水溃灭完全不同。在出水阶段中,航行体两侧水面附近出现了反向旋转涡对(counter-rotating vortex pair,CVP),在反向旋转涡对的作用下出水部分空泡快速发生了大规模溃灭,直至水下部分空泡缩小至孤立空泡后,产生巨大的溃灭载荷。如图11所示,随着出水攻角增大,背流侧空泡末端附近的涡旋结构强度增大,从而导致脱落的空泡尺度越大。同样地,在大攻角状态下,水面附近的反向旋转涡对导致出水空泡发生了大规模溃灭现象。
图10 AOA0=4°,σ=0.198下出水速度矢量场演变Fig.10 Evolution of the velocity vector field at A AOA0=4° and σ=0.198 during the exiting-water process
图11 AOA0=8°,σ=0.198下出水速度矢量场演变Fig.11 Evolution of the velocity vector field at AOA0=8° and σ=0.198 during the exiting-water process
图12所示为AOA0=0°,σ=0.198下空化区涡旋结构细节图。由于λci方法在水下发射过程中涡识别过程表现较好。因此,本文采用λci准则对空化区域附近的涡旋结构进行了识别,其数学定义为
图12 空化流场涡结构细节Fig.12 Cavitation flow field vortex structure detail
(12)
(13)
(14)
当Δ>0,其特征值为λ1=λr,λ2,3=λcr±iλci;其中,
(15)
P=-(λ1+λ2+λ3)=-tr(∇V)
(16)
(17)
R=-det(∇V)
(18)
式中:P,Q,R为速度梯度张量∇V的三个伽利略不变量;tr为矩阵的迹;det为矩阵的行列式。值得注意的是:在本文中航行体头部触及自由液面之前,空化涡的流体介质是水和水蒸汽的混合物;随后在自由表面附近,水汽混合介质会掺混空气,此时空化涡的流体介质是水蒸汽、水和空气的混合物。本文重点关注了自由液面之下涡结构的演化机制。
由于肩部空化区域的存在,以发卡涡为代表的壁面涡旋结构发展明显受到了抑制。在壁面涡结构演变过程中,主涡环从空泡闭合区域脱落,涡管沿着流向发生拉伸形成涡腿,涡环与涡腿共同构成典型的发卡涡结构。随着涡结构的进一步发展,壁面发展并衍生出“发卡涡”相互锁定的“发卡涡包”,其形态具体表现为不规则状。另外,发现多个发卡涡沿轴向间隔排列,组成发卡涡包存在于航行体尾流中,而尾涡结构的演变将直接决定连续发射领域中后一发航行体的出水姿态稳定性[20]。
如图13和图14分别为AOA0=4°和AOA0=8°下出水空化涡结构演变过程。从图13、图14中可以发现,随着攻角增大,壁面发卡涡表现生成时间早、尺度大以及发展规则性更差。对比A区和C区发现,随着攻角增大,涡环的数量和尺度明显增加,涡腿长度变得细而长。对比B区和D区发现,随着攻角增大,尾涡中发卡涡的尺度也明显增加。
图13 AOA0=4°,σ=0.198下出水空化涡结构演变Fig.13 Cavitation flow field vortex structure detail at AOA0=4° and σ=0.198
图14 AOA0=8°,σ=0.198下出水空化涡结构演变Fig.14 Cavitation flow field vortex structure detail at AOA0=8° and σ=0.198
在上述空化流场结构和涡结构分析中,发现出水空泡在自由液面附近会发生强烈的坍塌溃灭行为。如图15所示为出水空泡溃灭过程示意图,首先航行体肩部附着空泡在自由液面上方反向旋转涡对作用下,空泡顶端溃灭现象开始发生。当溃灭位置到达回射流顶端位置时,水下空泡快速收缩并在回射流涡旋结构作用下发生大尺度脱落。当缩小至一个尺度较小的孤立空泡时,该空泡整体溃灭将直接冲击航行体表面,将形成大范围的较高压力脉动峰值。因此,出水空泡溃灭压力特性研究将具有重要的意义。如图16所示,在航行体肩部附近布置多个探针点监测溃灭脉动压力演变。
图15 出水空泡溃灭过程示意图Fig.15 Diagram of cavity collapse process
图16 航行体表面监测点分布Fig.16 Distribution of monitoring points on the surface of the vehicle
如图17和图18所示分别为AOA0=0°下,σ=0.198和σ=0.231不同监测点压力演变曲线。从图17(a)和图18(a)中可以发现,压力曲线前期出现了短暂而快速地下降过程,主要是航行体肩部空泡开始生成。在航行体还未出水之前,大部分监测点压力曲线基本平稳。然而,在出水阶段空泡溃灭导致监测点压力均出现剧烈的脉动变化。图17(b)和图18(b)所示分别为σ=0.198和σ=0.231下监测点压力曲线局部放大图。在出水空泡溃灭过程中,出现了面积范围大、脉动幅值高的压力峰值。特别需要注意的是,当空泡缩小至孤立空泡随后发生溃灭行为,表现出脉宽小和峰值高特点。然而,在此瞬态时刻下产生了巨大的冲击压力峰值,约为185个大气压和155个大气压。当此溃灭压力冲击航行体壁面,将会给其结构安全性带来极大的破坏。
图17 AOA0=0°, σ=0.198下不同监测点压力演变Fig.17 Pressure evolution with different monitoring points at AOA0=0° and σ=0.198
图18 AOA0=0°, σ=0.231下不同探针压力演变Fig.18 Pressure evolution with different monitoring points at AOA0=0° and σ=0.231
如图19所示为σ=0.198下不同攻角下背流面监测点压力演变曲线。从图19(a)可知,4°攻角下出水过程中未出现大范围的脉动压力峰值。相比0攻角状态下,出水空泡溃灭压力峰值范围和最大冲击压力峰值均减小;从图19(b)可知,8°攻角下的最大冲击压力峰值范围相比4°攻角下有所增大。图15所示为σ=0.198下不同攻角下迎流侧监测点压力演变曲线。从图20(a)和图20(b)可知,迎流侧与背流侧压力演变规律有所不同。在空泡生成阶段,不同攻角下迎流侧压力均出现了较大的脉动峰值。综上所述,航行器出水姿态与空泡溃灭载荷幅值紧密相关。
图19 σ=0.198下不同攻角下背流侧监测点压力演变Fig.19 Pressure evolution at the back flow monitoring point with different angles of attack at σ=0.198
图20 σ=0.198下不同攻角下迎流侧监测点压力演变Fig.20 Pressure evolution at the face flow monitoring point with different angles of attack at σ=0.198
对航行体出水过程空化结构演变与溃灭载荷特性进行了数值模拟研究,得到的主要结论如下:
(1)在自由液面较大的介质密度差和回射流共同作用下,附着空泡迅速自上而下发生坍塌溃灭现象;在溃灭末期空泡收缩为较小的孤立空泡时,在外界高压作用下孤立空泡发生了溃灭现象,形成的高速射流冲击结构表面。
(2)航行体出水过程中,水面附近的反向旋转涡对导致附着空泡发生了溃灭现象;由于附着空泡存在,以发卡涡为代表的壁面涡旋结构发展明显受到了抑制。随着涡结构的进一步发展,航行体壁面发展出“发卡涡”相互锁定的“发卡涡包”,其形态具体表现为不规则状。
(3)0攻角状态下,当附着空泡缩小至孤立空泡随后发生溃灭行为时,冲击载荷表现出脉宽小和峰值高特点。同时产生了巨大的冲击压力峰值,将会给其结构安全性带来极大的破坏。当航行体出水时带有一定攻角下(0°~8°),空泡溃灭压力的范围和峰值均减小。但是,在大攻角状态下,初始时刻迎流侧空泡脱落以及溃灭行为会导致较大的压力峰值出现。