级间热分离流场干扰特性及动力学数值研究
0 引言
对于一个两级的导弹,级间热分离过程是在两级连接处切割解锁之前,上面级发动机就开启点火程序,一级弹体在二级发动机喷焰的作用下与二级弹体分离。级间热分离过程涉及到的流动现象十分复杂,级间间隙附近的横向喷流与外流场自由来流相互之间会形成强烈的干扰现象[1],这种干扰现象成为分离过程中上面级弹体姿态不稳定变化的重要原因[2]。
为了研究级间热分离过程的复杂流动现象,国内外开展了相关风洞试验和数值计算。Hohn等[3]利用红外成像技术分析了一级反推火箭羽流对二级气动热载荷的影响。Applebaum等[4]采用3种模型进行建模,讨论了不同模型对远场冲击和空腔流动等情况的适用性。林敬周等[5]在φ1米高超声速风洞中对多喷流干扰条件下级间分离流场及其热环境变化规律进行了风洞试验研究,清晰观测到级间包含了分离激波、弓形激波和边界层干扰等现象的流场结构。Luchinsky等[6]提出了一种基于物理的级间分离再接触故障建模方法,重建了发动机启动过程中受损喷嘴在羽流载荷下的屈服和熔化等流程。Stahl等[7]利用纹影图像技术和油流技术,在三声速风洞中对高速导弹控制中侧向射流对超声速横流的干扰效应,总结了影响测射流前方超压分离区范围大小的主要因素,但未能通过壁面压力测量获得法向力,更无法对法向力振幅进行定量分析。Stahl等[8]对冷热喷流与超声速来流的干扰作用进行了试验研究,结果表明,相较于冷喷流,热喷流前方的高压分离区作用范围更大,而低压区则小得多。王志坚等[9]采用冷喷流模拟技术对助推器分离过程进行了风洞试验研究,分析了喷流干扰对二级姿态控制的影响情况。Lungu等[10]在超声速风洞中利用火花纹影技术可视化流场,发现增加连续阶段出口压力能够增加剪切层的厚度,而随着级间距离的增加,各种流动特征均逐渐减弱。王常悦等[11]分别采用风洞试验和数值模拟2种方式,对分离过程中横向喷流和外部来流的干扰情况进行研究,统计了来流攻角对二级弹体气动特性的影响。梁杰等[12-13]利用DSMC数值算法,对高空条件下稀薄大气与三维平板模型侧向射流的干扰流场,探究了来流稀薄程度、射流压比以及飞行姿态等参数对流动分离区域范围和分离涡强度、力特性等干扰流场特性的作用规律。朱中根等[14]利用数值模拟研究了稠密大气下带斜喷动力导弹在喷流干扰效应中的涡流结构和波系结构,总结了攻角和来流马赫数等参数对干扰效应的影响规律,结果表明,喷流位于迎风侧和背风侧时喷流干扰对喷流推力和干扰力矩起的作用各异。李易等[15]采用结构化Chimera技术对分离体的运动进行了建模并进行了飞行力学分析。文献[16-17]分别对化学非平衡流模型、单组份以及双组分气体模型对固体火箭级间分离进行仿真,结果表明,使用双组分模型可以避免非平衡流模型的巨大计算开销,又能确保相对较高的保真度,其中文献[16]讨论了一级前封头导流板形状对干扰流场的影响,发现锥形和球形导流板可以减小二级附近流动分离区域的长度。傅德彬等[18]采用非定常流动和刚体六自由度耦合模型对级间热分离过程中的分离载荷和喷流状态进行了数值模拟,总结了远场冲击干扰的作用规律。
本文中利用Fluent软件,对非定常流场和刚体运动进行耦合仿真计算,根据计算结果探究级间分离过程中横向喷流与超声速自由来流耦合作用的变化规律,分析弹体所受载荷及运动情况,计算分析结果可以用于预测弹体运动轨迹,为实现级间安全分离提供一定的参考价值。
1 数学模型与计算方法
1.1 几何模型
由于级间分离过程只有二级发动机点火,一级弹体被适当简化,略去了一级弹体的发动机喷管,导弹模型由上下两级弹体和二级发动机喷管组成,如图1所示。弹体直径D=1 200 mm,一级、二级弹体长度分别为6.89D和6.63D,α、β分别为攻角和侧滑角。
图1 计算模型示意图
Fig.1 Schematic diagram of calculation model
在发出控制指令使得两级解锁以前,上下两级弹体之间的级间间隙为0 mm,处于密封状态,但在数值模拟过程中如果让两级初始的级间间隙为0,基于现有的数值模拟方法,无法正常划分计算网格,无法使得在发出级间解锁指令后保证两级正常分离,所以为了仿真过程的顺利进行,需要在级间段建模时留有一定尺寸的微小间隙,在计算模型中将初始间隙设置为10 mm。
1.2 控制方程和耦合计算模型
级间热分离过程存在级间间隙横向喷流与外部来流的耦合等典型非定常流动问题,流动状态复杂,需要采用非定常可压缩的Navier-Stokes方程作为流动控制方程,方程表示为
(1)
式(1)中:Q为守恒变量;F、G、H为无粘通量;Fv、Gv、Hv为粘性通量。
在级间分离初期的较短时间内,忽略燃料消耗造成弹体质量的变化,两级弹体均可视为质量和转动惯量恒定的刚体,刚体在六自由度运动中,外力由发动机推力、重力和空气动力组成,空气动力分解为升力、阻力和侧向力。将质心运动的动力学方程展开为标量形式,表示为
(2)
根据动量矩定理: M=dh/dt。由于模型在纵平面上是对称的,因此Ixy=Ixz=0。将质心旋转动力学方程展开为标量形式,表示为
(3)
欧拉角约束方程为:
(4)
式(3)、式(4)中:(u,v,w)为飞行速度在弹体坐标系中的投影;(p,q,r)为角速度在弹体坐标系中的投影;(L,D,C)为速度坐标系中定义的气动分量;(α,β)为攻角和侧滑角;mg为重力;(L,M,N)为转矩在弹体坐标系的投影;(Ix,Iy,Iz)为转动惯量;Iyz为惯性积;(φ,θ,ψ)表示弹体的姿态坐标。
对非定常可压缩Navier-Stokes方程和刚体运动方程进行耦合求解,求解流程如图2所示。在流固耦合计算中利用结合(用户自定义函数)程序,将两级弹体设定为刚体,刚体运动造成了流场边界位置和速度的变化,而对流动的求解获得了分布作用在刚体上的载荷,对分布的载荷进行表面积分获得刚体所受合力,进而得出刚体加速度,并对时间进行积分获得刚体速度和位移。
图2 耦合求解流程图
Fig.2 Flow chart for coupled solution
1.3 计算域网格划分
计算域网格划分情况如图3所示,在喷管内部采用了结构网格,其他区域则采用非结构网格。级间段附近的流动情况十分复杂,需要网格进行局部加密,加密区域的网格最小尺度保持在2 mm。计算域网格数量在400万左右。级间分离过程两级弹体之间存在相对运动,因此必须采用动网格技术来模拟级间分离运动。利用弹簧式光顺方法以及局部网格重构法实现动态网格的变形,每次迭代后,根据边界位置的变化,网格出现变形,变形规律与弹簧伸缩类似,每次变形后均进行网格质量检查,当质量不合格时进行局部重构,直到生成新质量合格的网格。图3(b)和图3(c)分别为分离进行到10 ms和36 ms时的级间段网格。
图3 计算网格示意图
Fig.3 Schematic diagram of calculation mesh
1.4 离散格式
数值计算中,先要对控制方程中的微分项进行离散处理。时间离散格式采用隐式算法,空间离散格式采用一阶迎风格式和Roe-FDS格式。
迎风型格式本身具有格式黏性,能将涨落和流量等特性较好地表现出来,在提高间断分辨率和计算效率上效果显著。迎风格式精度不高,但是对于对流主导的流动,迎风格式是最稳定的数值格式,可以使用这一格式来产生一个初始稳定的解,然后再把数值格式改成其他精度更高的格式来使得流场更加精确,通量差分分裂(Roe-FDS)格式则对激波和间断均具有较高的辨识率。所以本文中采用一阶迎风格式和Roe-FDS格式进行空间离散。
1.5 湍流模型
级间分离过程中的流动属于湍流流动,数值仿真过程中一般采用湍流模型对湍流运动进行模化。本研究中,应用S-A单方程湍流模型求解运动湍流粘度的模型输运方程。S-A模型能够有效承受不利压力梯度的边界层,广泛应用于涉及壁面有界流动的航空航天问题。在Ansys Fluent中对S-A模型进行了扩展,采用了不敏感壁面处理,使得模型的应用不依赖于近壁分辨率。
1.6 计算状态
飞行高度为H= 20～35 km,来流马赫数Ma∞=4,攻角α=2°,侧滑角β=2°,二级发动机入口总压7.0 MPa,总温3 500 K。不同飞行高度工况如表1所示。
表1 不同高度工况列表
Table 1 List of working conditions at different altitudes
在0 ms时,二级导弹发动机启动点火程序,当二级喷管喉部憋压至2 MPa时,将喉部堵片打开,2 ms后发出指令,级间段切割解锁,两级开始分离。
计算采用有限体积法,湍流模型选用S-A模型,空间离散采用一阶迎风格式和Roe-FDS格式,时间离散采用隐式算法,时间步长设置为1×10-5 s,残差收敛条件设置为1×10-3,计算时间持续100 ms。
2 计算结果分析
2.1 网格无关性分析
通常情况下,计算域网格划分越精细,数值模拟的可靠性越高,计算结果越趋近于真实情况,但也会造成计算的规模和占用的存储空间增大,不利于高效率计算,因此,需要选取合适精度的网格,兼顾计算效率和结果的准确性。
为验证选取网格的合理性,划分4套不同网格,最小网格尺寸和网格总数见表2。分别利用每套网格计算飞行高度为H=30 km,来流马赫数Ma∞=4,攻角α=2°,侧滑角β=2°的条件下级间段内部监测点处压强变化情况,监测点位置如图4所示,监测点位于xOy截面上二级弹体后方内壁面,距分离面0.2D。
表2 不同网格最小尺寸及网格总数
Table 2 Minimum size and total number of different meshes
图4 监测点位置
Fig.4 Location of monitoring point
利用不同网格计算得到监测点处压强变化曲线如图5所示,在100 ms的计算时间内,级间段内部监测点处压强呈现先上升后下降的趋势,在20 ms时达到峰值。网格最小尺寸为5 mm时的监测点处压强变化曲线误差较大,峰值处误差超过30%;随着网格最小尺寸减小,不同网格的计算结果逐渐接近,网格Ⅲ与网格Ⅳ压强峰值仅相差1.7%,网格Ⅲ的网格数目则比网格Ⅳ少43%,因此为了节约计算开销同时保证计算精度,将2 mm(网格Ⅲ)作为本研究网格最小尺度。
图5 监测点处压强曲线对比
Fig.5 Comparison of pressure curves at monitoring points
2.2 级间段内部流场结构
本文中二级发动机尾喷管采用的是收敛-扩张喷管,对应着完全膨胀、欠膨胀与过膨胀3种工作状态。喷管会处于何种工作状态主要由其可用压力比和设计压力比之间的大小关系决定,当喷管出口背压等于设计压力时,可用压力比与设计压力比相等,喷管出口气体完全膨胀;当出口背压小于设计压力时,可用压力比大于设计压力比,喷管处于欠膨胀工作状态;当出口背压大于设计压力时,可用压力比小于设计压力比,喷管处于过膨胀工作状态。
收敛-扩张喷管过膨胀状态如图6所示,在本文中的计算状态下,喷管扩张段内的燃气由于过度膨胀,在到达出口前气体压力便下降到一个很低的值,为了使燃气压力能够与外界环境压力相平衡,最终会在扩张段内形成一道正激波,同时激波位置会随出口背压的改变而移动,背压升高,激波向靠近喉部截面方向移动;背压降低,激波向喷管出口截面方向移动。
图6 收敛-扩张喷管过膨胀状态
Fig.6 Overexpansion state of convergent-divergent nozzle
图6(a)为过膨胀状态的收缩扩张喷管示意图和喷管内部轴线上马赫数和压强分布曲线,图6(b)为仿真获得的对应云图和曲线,仿真结果与理论曲线吻合较好。在过膨胀状态下,喷管出口背压大于设计压力,可用压力比小于设计压力比。从图6可以看出,气流穿过正激波后,马赫数瞬间降至1以下,气流速度为亚声速,静压强瞬间上升。
图7、图8给出了20～90 ms时间段级间段内马赫数云图和压强云图,图9为20～90 ms内喷管中心轴线上压强分布图,由于分离距离的增大,喷管中心轴线横坐标也不断延长。
图7 级间段内部流场马赫数云图
Fig.7 Mach number cloud map of internal flow field in the interstage
图8 级间段内部流场压强云图
Fig.8 Pressure cloud map of internal flow field in the interstage
图9 喷管中心轴线压强分布
Fig.9 Pressure distribution at the central axis position of the nozzle
从图7—图9可以看出,随着分离过程的进行,级间距离增大,喷管出口背压整体上呈现不断下降的趋势,进而导致喷管内正激波强度不断降低,正激波位置向出口截面方向移动,在分离进行到70 ms时,正激波已达到出口截面位置。
2.3 级间段间隙喷流干扰流场结构
文献[2]给出了级间分离过程中级间段内外典型流场结构示意图,如图10所示,级间间隙的横向喷流与超声速自由来流产生非常复杂的流动干扰效应。外界超声速自由来流在横向喷流的阻碍作用下,在外流场无粘区域形成一道高强度的弓形激波,较高的逆压梯度会造成来流上游的湍流边界层与上游外壁面发生分离。超声速来流遇到分离的边界层后发生转向,在流动分离点上方产生分离激波。间隙处的高温高压横向喷流在离开间隙后迅速膨胀加速,喷流压力迅速降低。经过在主流中的膨胀过程后,高温高压气流转变成为低压高速射流,最终会形成能够用于平衡外界环境动量和压力的羽流激波。
图10 级间段内外典型流场结构
Fig.10 Typical flow field structures inside and outside the interstage
图11是在来流马赫数Ma∞=4、飞行高度H= 30 km条件下,计算获得的在36 ms时级间间隙横向喷流与外界超声速来流形成的干扰流场结构。从图11中可以看出,计算获得的流场结构与图10吻合较好,具有弓形激波、分离激波以及羽流激波、回流区、剪切层等主要流动特征,在流动分离点的下游方向、排气羽流的后方、喷管出口下面级发动机前封头附近的低速区域为流动分离区,均出现了旋涡结构。
图11 级间段内外流场云图和流线图
Fig.11 Cloud map and streamline diagram of the flow field inside and outside the interstage
在级间段外壁面附近设置一条监测线,用于获取外壁面压强分布,监测线位置如图12所示。在36、44、60 ms等3个时刻监测线上压强分布曲线如图13所示,3个时刻级间段附近的马赫数云图如图14所示。
图12 监测线位置
Fig.12 Location of monitoring line
图13 监测线上压强分布曲线
Fig.13 Pressure distribution curves at monitoring lines
图14 不同时刻级间段附近马赫数云图
Fig.14 Mach number cloud map near interstage at different times
超声速来流在遇到分离的边界层后会发生转向,形成分离激波,从而导致流动分离区域的压力要比分离激波上游区域的压力高,因此上面级表面流动分离点对应级间段外壁面压力突然上升的位置,在36、44、60 ms等3个时刻流动分离点的位置分别为x/D=5.85、x/D=5.75、x/D=6.15。在流动分离点的右侧,各条监测线上的压力均呈现为水平直线,这是因为对应位置为流动分离区,流动分离区内压力处处相等,近似等于分离点处压力大小。
进一步分析图13、图14可以发现,当分离时间从36 ms增加到44 ms时,随着级间间隙的增大,上面级表面流动分离点向上游移动,从x/D=5.85的位置逐渐向上游移动至x/D=5.75的位置,但流动分离区域的压力大小基本保持一致,均在13 000 Pa左右,说明此过程级间间隙的变化主要使流动分离点位置发生改变,而对分离激波的强度影响较小。当分离时间从44 ms增加到60 ms时,上面级表面流动分离点向下游移动至x/D=6.15的位置,并且分离流动区域压力明显降低至不足11 000 Pa,这是因为当级间间隙大到一定程度时,上面级发动机尾喷管排出的燃气在下面级发动机前封头处的反溅效果变差,导致大部分高温高压气体离开间隙后沿下游方向膨胀加速,只有少部分气体沿侧向膨胀,导致横向喷流干扰强度减弱。因此当级间间隙达到一定值时,继续增大级间间隙会促进二级弹体表面流动分离点向下游移动,并减小分离激波的强度,造成二级弹体表面分离区压强降低。
在36 ms,受横向喷流的影响,外流场监测线上的压强在上面级表面流动分离区后方的排气羽流的作用范围内(6.6<x/D<6.8)急剧上升,监测线在排气羽流后方的位置处于一级弹体表面的流动分离区,监测线上压强再次下降,一级弹体表面流动分离区内压力接近于0。随着级间间隙增大,分离激波、排气羽流作用范围均向下游移动,在60 ms时,羽流激波与监测线相交于x/D=6.9和x/D=7.1处,在该范围内监测线上压强峰值仅为约45 000 Pa。
2.4 飞行高度对横喷干扰流场作用分析
分离进行到40 ms时,两级弹体间的距离已达28 mm,此时不同飞行高度条件下干扰流场在该时刻的流线分布如图15所示。
图15 干扰流场流线分布
Fig.15 Streamline distribution of interference flow field
在飞行高度为20 km时,级间间隙上游未出现明显的流动分离现象;随着飞行高度的增加,来流动压降低,间隙处横向喷流膨胀效应增强,当飞行高度增大至30 km时,能观察到二级弹体后侧(级间间隙上游)壁面附近出现流动分离和再附,形成一对驻涡,前者为流动分离区域内顺时针的主上游涡,壁面附近流动方向与自由来流相反,后者为再附区内逆时针的次上游涡,壁面附近流动方向与自由来流相同;飞行高度增大至35 km时,上游分离区向上游延伸,面积不断增大。
一级弹体前侧(级间间隙下游)的壁面附近也出现了流动分离现象,流动流动分离区域内形成了顺时针的主下游涡,随着飞行高度的增加,一级弹体前侧附近的流动分离区域不断增大,并逐渐向下游扩张。
2.5 分离载荷与动力学分析
在飞行高度H=30 km、来流马赫数Ma∞=4、攻角α=2°、侧滑角β=2°的条件下,计算获得一级弹体沿轴向的受力曲线和位移曲线如图16所示。
图16 一级弹体轴向载荷和位移曲线
Fig.16 Axial load and displacement curve of the first stage
由图16可知,在100 ms的计算时间内,一级弹体首先受到二级发动机燃气的作用,在初始阶段,级间间隙较小,间隙喷流流量小于喷管出口流量,级间段内处于憋压状态,一级弹体轴向载荷不断上升,在20 ms时刻轴向受力上升至峰值903 kN左右,一级弹体在载荷作用下沿x方向加速运动。随着分离距离的增大,级间间隙流量超过喷管出口流量,级间处于泄压状态,一级弹体所受轴向载荷持续下降。当分离进行到36 ms时,一级弹体轴向载荷出现了小幅度波动,此时分离距离已达0.22 m,级间段间隙内外压强差减小,造成在36～40 ms的时间段内泄压速度低于憋压速度,一级弹体轴向载荷在小幅度上升之后继续下降。在分离运动进行到100 ms时,两级轴向距离达到1.5 m左右。
2.6 一级弹体侧向运动分析
图17、图18为一级弹体在x、y、z等3个方向所受载荷曲线和位移曲线。一级弹体侧向(y、z方向)载荷来源主要为一级弹体外表面所受压力和摩擦力的合力,而轴向(x方向)载荷主要来源于所受二级发动机燃气作用力,因此一级弹体侧向载荷远小于轴向载荷,侧向位移也远小于轴向位移。
图17 一级弹体载荷曲线
Fig.17 Load curve of the first stage
图18 一级弹体位移曲线
Fig.18 Displacement curve of the first stage
在不同飞行高度条件下,一级弹体y方向和z方向位移曲线分别如图19、图20所示。由于重力的存在,一级弹体y方向位移始终为负,随着飞行高度从20 km增加至35 km,来流越来越稀薄,升力逐渐降低,抵消重力效果减弱,一级弹体在y方向的位移绝对值从0.035 m不断增大至0.047 m。在飞行高度为20 km时,一级弹体在z方向的位移为正,在100 ms达到0.011 m,在分离后期,随着飞行高度增加至30 km以上,弹体逐渐向z轴负方向移动。
图19 一级弹体y方向位移曲线
Fig.19 The y-direction displacement curve of the first stage
图20 一级弹体z方向位移曲线
Fig.20 The z-direction displacement curve of the first stage
3 结论
本文中以多级导弹的级间热分离过程为研究对象,利用CFD软件Fluent结合动网格技术对非定常可压缩N-S方程和刚体运动方程进行了耦合求解,对级间热分离过程的流场结构与一级弹体的动力学情况进行了数值模拟,通过数值仿真还原了级间段内部喷管过膨胀状态的喷流流场以及间隙横向喷流与外界自由来流的干扰流场,针对流场中激波结构、排气羽流、流动分离等三维非定常流场结构的形成原因和变化规律进行深入分析,并对一级弹体的动力学情况进行了研究。研究结果表明:
1) 在本文中的计算条件下,喷管处于过膨胀状态,随着分离距离的增大,出口背压不断降低,喷管内的正激波不断向喷管出口截面方向移动。
2) 级间分离过程中级间间隙处的横向喷流与超声速自由来流形成复杂的流动干扰现象,在二级弹体的尾部和一级弹体前侧附近均形成流动分离区域,随着分离过程的进行,级间间隙气流逐渐向下游膨胀,而侧向的横喷干扰强度减弱,分离激波和羽流激波强度降低并向下游移动。
3) 在20～35 km的飞行高度范围内,随着飞行高度增大,随着飞行高度的增加,来流动压降低,级间间隙处横向喷流膨胀效应增强,一级弹体、二级弹体表面流动分离区域面积均不断扩大。
4) 分离初期,级间段处于憋压状态,一级弹体轴向所受载荷迅速上升;随着分离过程的进行,级间间隙达到一定距离,间隙喷流流量超过二级喷管出口流量,级间段内压强下降,一级弹体轴向受力达到峰值并迅速下降,之后经历小幅度波动,继续缓慢下降。
5) 一级弹体侧向受力和位移远小于轴向受力和位移,随着飞行高度的增加,y、z方向位移均向负方向增大。