0 引言
为了能够以最低成本将有效载荷送入轨道,必须提高喷管性能。只有当喷管出口处压力等于环境压力时,喷管效率最大。传统喷管由于面积比固定,具有恒定的出口压力,因此只能在设计高度下提供最佳性能。喷管的性能损失来源有化学非平衡损失、摩擦损失、喷管出口流动不均匀造成的损失、燃料混合和燃烧不充分造成的损失以及喷管流动的非适应性损失。其中由于气流对外界环境压力变化不适应而产生的非适应性损失占比最大,约为全部损失的15%[1]。膨胀偏流喷管(expansion deflection nozzle,ED)是一种具有高度补偿能力的新型喷管,因其长度短、质量轻、没有可移动部件[1]等优点而受到广泛关注。
1960年,Rao[2]首先提出了ED喷管这一概念,采用无粘流数值计算方法分析了喉部流动,并按照最大推力型面设计ED喷管,进行了试验验证。然而随着多级火箭技术的广泛应用,ED喷管的研究陷入停滞。进入21世纪,为了能够减轻发动机质量、降低发射成本,以及开发可重复使用的单级入轨飞行器(SSTO),研究者们对ED喷管的流场特性和高度补偿性能进行了一系列数值和试验研究。Taylor提出了ED喷管喉部的一般设计方法[3],并总结了不同的设计原则对ED喷管性能的影响[4-5]。Wagner等[6-7]等设计了一个平面ED喷管,并进行冷流试验,研究了ED喷管出口附近的流场结构和流动模式的转变。Schomberg等[8-9]研究了几何结构参数对推力性能的影响,在试验测试条件下,所设计的平面ED喷管性能均高于钟形喷管。John Paul P等[10-11]采用数值模拟的方法研究了Wagner设计的平面ED喷管中心体基底形状对喷管流动和性能的影响,结果表明,具有最高基底曲率的弯曲基底整体性能最好。Wang等[12]对面积比100的ED喷管内的流动进行了数值模拟,以了解在工作模态转变过程中ED喷管推力效率(η)下降背后的机理。张琦[13]、李明肤[14]、张渴欣[15]、王革[16]等利用数值模拟研究了扩张比、中心体型面等结构参数对ED喷管高度补偿性能的影响。
为了能更好的理解ED喷管内的流动行为,探究ED喷管高度补偿的产生机理,在本文中通过数值仿真手段得到平面ED喷管在不同工作压比下的流场特性和推力性能,与同扩张比传统收敛-发散喷管(conventional convergent-divergent nozzle,CD)进行比较,并研究压力比、推力以及工作模式转换之间的关系。
1 计算模型和数值方法
1.1 计算模型
ED喷管喉部几何型面参考Taylor[3]提出的设计方法进行设计,给定喉部最小距离Gt=3 mm,喉部各圆弧半径均为Gt的倍数。根据文献[4,8]的参数优化结果,选择较大的喉部偏转角θt,选择合适的中心体壁面半径使得喉部面积缓慢增大,并且中心体最后一段圆弧半径
较大,与文献[8]中性能最佳的结构15一样为10Gt。扩张段型面为3次抛物线,设计了一个面积比为30的平面ED喷管。各结构参数符号表示如图1所示,几何模型采用的详细几何参数值列在表1中。
图1 ED喷管几何模型
Fig.1 Sketch of ED nozzle geometry
表1 ED喷管计算模型几何参数
Table 1 ED nozzle design parameters
几何参数符号数值喷管壁面半径R-w/Gt2R+w/Gt2中心体壁面半径R-p/Gt5R+p/Gt10中心体半径yd/Gt1.5喷管出口半径Re/mm90扩张段长度L/mm130喉部偏转角θt/(°)60中心体弯曲角θi/(°)30喉部出口角θe/(°)10喷管出口扩张半角αe/(°)12
设定ED喷管的入口为压力入口;外流场长度20Re,宽度6Re,3个边界均设置为压力出口;其余壁面的边界条件均为绝热且无滑移,计算域如图2所示。压力出口条件为0.1 MPa,温度为300 K;通过改变入口总压来改变喷管的工作压力比,以表示不同的工作高度,入口总压在1～15 MPa之间变化,温度仍为300 K。定义喷管压力比(NPR)为进口总压P0与环境压力Pa之比,则进行仿真计算的压力比范围为10～150。利用ICEM对模型进行结构网格划分,对喷管喉部以及喷管壁面处进行加密处理,最终结果如图3所示。使用Fluent 2021 R2对模型进行仿真计算。
图2 ED喷管计算域
Fig.2 Computational domain of ED nozzle
图3 网格划分结果
Fig.3 Meshing of computational model
ED喷管的推力是喷管工作时内、外表面所受气体压力的合力,计算公式为
F=F内+F外
(1)
式(1)中,F内包括喷管入口处燃气贡献的推力、喷管内壁面所受燃气压力在轴线方向上的积分和中心体壁面所受燃气压力在轴线方向上的积分,F外为外界背压贡献的推力。Pw表示喷管壁面压力,Ae表示喷管出口面积。
为了消除喷管厚度对推力大小的影响,通过比较推力系数的大小来比较ED喷管与同扩张比CD喷管的推力性能,推力系数计算公式为
(2)
式(2)中,At表示喷管喉部面积。
1.2 数值方法
计算时假定流场为二维稳态模型,不考虑化学反应,不考虑辐射传热作用,假定燃气为理想气体且流场中不含固体颗粒。其中粘度采用三系数的Sutherland定律[17]定义。利用雷诺平均Navier-Stokes方程[18]求解流场,控制方程为
连续方程:
(3)
动量守恒方程:
(4)
能量守恒方程:
(5)
式(3)—式(5)中:
是速度矢量,p是静压,ρ是密度,μ是动力粘度,I是单位张量,kT是热导率,E是总能,T是温度,μt是湍流粘度,Prt是湍流普朗特数。能量守恒方程的最后一项表示粘性热,且μeff=μ+μt。
采用Menter[19]提出的两方程剪应力输运(SST)k-ω湍流模型进行仿真计算。SST k-ω模型是基于湍流动能k和比耗散率ω的模型方程,具体表示如下
(6)
(7)
式(6)、式(7)中,Gk表示由平均速度梯度产生的湍流动能,Gω表示由ω方程产生的湍流动能,Γk和Γω分别表示k和ω的有效扩散系数,Yk和Yω分别表示k和ω由于湍流作用引起的耗散,Dω表示交叉扩散项。
选择基于压力的耦合求解器求解控制方程,采用基于最小二乘单元的空间离散化,其中假设解呈线性变化。对流项采用二阶上风插值方案进行求解。计算分析是在稳态条件下进行的,初始化方式为混合初始化。在计算过程中,对喷管出口质量流量进行监控,当残差下降至10-3以下,出口质量流量无明显变化时,判定计算收敛。
2 结果与讨论
利用数值分析方法模拟了平面ED喷管在不同工作压力比下的流场特征,获得了流动模式、和推力系数变化规律,并与CD喷管进行比较,具体结果如下。
2.1 ED喷管流动行为分析
ED喷管在不同工作压力比下具有不同的工作状态,对应的流场特征也不相同。图4(a)、(b)、(c)展示了ED喷管在NPR分别为30、60、100时流场的马赫数云图和喷管内壁面压力分布曲线。NPR=30时,气流经喉部膨胀,在中心体圆弧末端产生一道入射激波,到达喷管扩张段壁面后发生反射,气流的流动方向发生偏转。此时喷管内流动被剪切层一分为二,剪切层上方是超声速气流,其流动范围始终被限制在内壁面和剪切层之间;剪切层下方,即中心体后方是开放的粘性再循环区,ED喷管处于“开尾流”工作模式。从内壁面压力分布曲线中可以看出,在扩张段中部,激波反射导致超声速流动发生压缩转向,引起壁面压力小幅度升高,在反射点之后,压力又逐渐下降至略低于环境压力,在出口处产生一道斜激波,使壁面压力与环境压力相适应。
图4 ED和CD喷管不同压力比下流场和壁面压力分布
Fig.4 Flow field and pressure distributions of the ED nozzle and CD nozzle in different NPRs
NPR=60时,入射激波在扩张段壁面的反射点移动到喷管外部,超声速气流的流动面积增加,粘性再循环区的宽度变窄,喷管内壁面压力在喉部前后发生骤降,在扩张段部分压力缓慢下降,ED喷管仍处于“开尾流”工作模式。NPR=100时,超声速气流在中心体末端发生转向,剪切层附着在喷管中心轴线所在壁面上,产生再附着激波。此时,粘性再循环区面积很小,被限制在中心体后方和超声速气流之间,与外界环境隔离,ED喷管此时为“闭尾流”工作模式。喷管内壁面压力分布曲线与NPR=60相似。此后,随着NPR继续增大,流场将始终处于闭合尾流状态,内壁面压力分布曲线也与NPR=60相似。
图4(d)、(e)、(f)分别为同扩张比CD喷管在NPR=30、60和100时流场的马赫数云图和喷管内壁面压力分布曲线。CD喷管的几何尺寸如图5所示。其中喉部半径Rt=Gt=3 mm,收敛段入口半径Rc=3Rt,收敛半角β=30°,喉部收敛半径Rd=3Rt,喉部扩张半径Rm=Rt,初始膨胀半角αm=29°,出口膨胀半角αe=17°,扩张段长度为ED喷管的2倍。
图5 CD喷管几何尺寸
Fig.5 Conventional convergent-divergent nozzle geometry
从图4(d)、(e)、(f)中可以看出,CD喷管在工作压力比较低时会发生壁面流动分离,随着NPR不断增大,分离点位置越靠近喷管出口,分离点压力也越小。NPR=30时,喷管内壁面压力持续下降至低于环境压力,直到发生壁面分离,分离点后压力上升至与环境压力相当,产生了一个压力抬升的台阶。随着NPR逐渐增大,压力抬升值逐渐降低,台阶高度减小。
比较ED喷管和CD喷管的流场特性和壁面压力分布曲线,可以发现:由于ED喷管中的超声速气流被限制在喷管内壁面和剪切层之间,因此没有像CD喷管一样出现壁面流动分离;增加中心体使得ED喷管内流动状态发生显著变化,流动结构更加复杂。
2.2 ED喷管流动模式转变过程分析
从图4中ED喷管流场的流线分布可以看出,剪切层下方的粘性再循环区会将外界大气“吸入”喷管内部,进而导致中心体后方的压力(简称“基底压力”,用Pb表示)略低于环境压力。随着NPR逐渐增大,剪切层逐渐向ED喷管中心线移动,再循环区气流被加速,剪切层的吸入效应越强,基底压力逐渐下降,当降低至某一数值时,超声速气流在中心体附近的膨胀状态就会发生变化,剪切层与轴线相交,流动模式发生转变。因此,基底压力的变化与ED喷管流动模式的转变密切相关。
图6显示了基底压力随NPR的变化规律。从图中可以看出,在NPR很小时,基底压力略低于环境压力,随着NPR逐渐增大,基底压力也逐渐下降;当NPR增大到90之后,降低速率变缓,且在NPR=96～97之间,Pb发生骤降,NPR=97时的Pb比NPR=96时的Pb降低超过3/4。结合图7所示2个压力比下的马赫数云图,可以说明Pb的突降对应着流动模式从开尾流转变为闭尾流。当流动模式转变为闭尾流后,基底压力又会随NPR增加而缓慢增大。这是因为此时中心体后方的再循环区被超声速气流包围,与外界环境隔离,此时基底压力仅受超声速气流影响,再循环区随着超声速气流的膨胀而不断被压缩,基底压力升高。
图6 基底压力随NPR的变化曲线
Fig.6 Base pressure at different NPRs
图7 NPR=96和97时ED喷管马赫数云图
Fig.7 Mach number contours of ED nozzle at NPR=96 and NPR=97
2.3 ED喷管推力性能分析
分别计算ED、CD喷管在压力比10～150范围内的推力系数,绘制得到的曲线如图8所示。从图中可以看出,在进行仿真的22个不同工况下,CD喷管的推力系数随NPR增加而增大;ED喷管的推力系数均高于CD喷管,但在NPR=96～97之间推力系数有所下降。ED喷管的推力系数相对于CD喷管的提升值随NPR的变化曲线如图9所示。在所有计算条件下,ED喷管的推力系数均有所提升,且提升值随压力比的增加整体呈现先降后升的趋势。NPR小于60之前,推力系数提升值逐渐下降,NPR=10时性能提升最大,推力系数提升了约0.29;NPR=60时性能提升最小,推力系数仅提升0.051 6。NPR大于60之后,推力系数提升值开始增大,但在NPR=96～97之间,推力系数提升值出现明显下降,随后上升幅度变缓。
图8 不同压力比下ED、CD喷管的推力系数比较
Fig.8 Comparison of the thrust coefficient of ED and CD nozzle at different NPRs
图9 不同压力比下ED喷管的推力系数提升值
Fig.9 The thrust coefficient increase of ED nozzle at different NPRs
图10展示了NPR=50、60和70时ED喷管的马赫数云图。
图10 NPR=50、60和70时ED喷管马赫数云图
Fig.10 Mach number contours of ED nozzle at NPR=50、60 and 70
从图10中可以发现,NPR=50时,入射激波会撞击喷管内壁面,并发生反射,使得超声速流动压缩转向;NPR=60时,入射激波恰好到达喷管内壁面的末端;NPR=70时,入射激波倾斜角度更小,没有在喷管内壁面发生撞击,而是在喷管出口后方与出口斜激波相交。结合图9可知,入射激波的撞击点恰好在喷管出口末端会导致一定的性能损失,使得推力系数增加值最低。
比较图7中NPR=96和97的马赫数云图可以发现,流动模式发生转变后,在中心体后方产生了一道朝向轴线的强激波,且流场内的最大马赫数显著增大。结合图8和图9分析可知,NPR=97时由于该激波的产生,导致喷管性能损失较大,推力系数有所下降,推力系数提升值明显下降。由此可见,ED喷管推力系数下降现象与流动模式发生转变有关。
总结ED喷管性能提升的原因有以下2点:
1) 由3.1节的流场分析可知,ED喷管内的超声速气流被限制在剪切层和内壁面之间,不发生流动分离,大大减小了ED喷管的过膨胀损失,因此性能优于CD喷管。
2) 压力比较低时,从中心体发出的入射激波撞击喷管内壁面,驱动气流发生压缩转向,引起壁面压强升高,使得ED喷管低空性能提高。
3 结论
利用数值仿真方法计算了平面ED喷管在压力比10～150间的流场马赫数、壁面压力和推力系数,分别与CD喷管比较,得到以下结论:
1) ED喷管有2种不同的工作模式,分别为“开尾流”和“闭尾流”模式。在“开尾流”模式下,流场分为超声速流动区、剪切层和开放粘性再循环区3部分;在“闭尾流”模式下,剪切层与中心轴线相交,产生再附着激波,粘性再循环区被超声速气流封闭。
2) ED喷管基底压力Pb的变化与流动模式的转变密切相关。随着NPR不断增大,剪切层的吸入效应越强,Pb不断减小,达到某一个压力比时Pb发生突降,与此同时流动模式发生转变。
3) 本文所设计的ED喷管在所有计算条件下的推力系数均高于CD喷管,且在较低压力比时优势更明显。ED喷管性能提升的原因为:没有因壁面流动分离而产生的过膨胀损失;从中心体发出的入射激波撞击喷管内壁面,驱动气流发生压缩转向,引起壁面压强升高。
[1]HAGEMANN G,IMMICH H,NGUYEN T V,et al.Advanced rocket nozzles[J].Journal of Propulsion and Power,1998,14(2):62-64.
[2]RAO GVR.Analysis of a new concept rocket nozzle[J].Journal of Liquid Rockets and Propellants,1960,2:669-682.
[3]TAYLOR N V,HEMPSELL C M.Throat flow modelling of expansion deflection nozzles[J].Journal of the British Interplanetary Society,2004,57:242-250.
[4]TAYLOR N V.An integrated approach to expansion deflection nozzle analysis[D].Bristol:University of Bristol,2002.
[5]TAYLOR N V.Experimental and CFD analysis of the expansion deflection nozzle[C]//7th European Symposium on Aerothermodynamics.2011,692:101.
[6]WAGNER B,STARK R,SCHLECHTRIEM S.Experimental study of a planar expansion-deflection nozzle[J].Progress in Propulsion Physics,2011,2:641-654.
[7]WAGNER B,SCHLECHTRIEM S.Numerical and experimental study of the flow in a planar expansion-deflection nozzle[C]//47th AIAA/ASME/SAE/ASEE Joint Propulsion Conference &Exhibit.2011:5942.
[8]SCHOMBERG K A,DOIG G,OLSEN J,et al.Geometric analysis of the linear expansion-deflection nozzle at highly overexpanded flow conditions[C]//50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference.2014:4001.
[9]SCHOMBERG K,OLSEN J,NEELY A,et al.Experimental analysis of a linear expansion-deflection nozzle at highly overexpanded conditions[C]//19th Australasian Fluid Mechanics Conference.2014:74-77.
[10]PAUL P.J,PAUL J,NAIR P P,et al.Numerical analysis of effects of pintle shape on flow in planar expansion-deflection nozzles[C]//ISABE2019.2019.
[11]PAUL P.J,NAIR P P,SURYAN A,et al.Numerical simulation on optimization of pintle base shape in planar expansion-deflection nozzles[J].Journal of Spacecraft and Rockets,2020,57(3):539-548.
[12]WANG G,CHEN L,ZHOU B,et al.Numerical investigation on thrust efficiency dropping phenomenon of annular expansion-deflection nozzles[J].Physics of Fluids,2021,33(12):126107.
[13]张琦,李明肤,王革,等.膨胀偏流喷管高度补偿机制数值研究[J].推进技术,2019,40(1):44-52.ZHANG Qi,LI Mingfu,WANG Ge,et al.Numerical study on altitude compensation mechanism of expansion-deflection nozzle[J].Journal of Propulsion Technology,2019,40(1):44-52.
[14]李明肤.膨胀偏流喷管的性能研究[D].哈尔滨:哈尔滨工程大学,2017.LI Mingfu.Performance research of expansion-deflection nozzle[D].Harbin:Harbin Engineering University,2017.
[15]张渴欣,周博成,陈磊,等.环喉式膨胀偏流喷管结构参数研究[C]//第六届空天动力联合会议暨中国航天第三专业信息网第四十二届技术交流会暨2021航空发动机技术发展高层论坛论文集(第2册).成都:中国科协航空发动机产学联合体,2022:327-334.ZHANG Kexin,ZHOU Bocheng,CHEN Lei,et al.Research on structural parameters of the floating nozzle[C]//Collection of the 6th Joint Conference on Space Power and the 42nd Technical Exchange Meeting and 2021 Aero-Engine Technology Development High-level Forum (Volume 2).Chengdu:Aero-engine Production and Student Consortium of China Association for Science and Technology,2022:327-334.
[16]王革,张渴欣,周博成,等.基于混合正交法的环喉式膨胀偏流喷管性能影响因素研究[J].固体火箭技术,2023,46(3):410-418.WANG Ge,ZHANG Kexin,ZHOU Bocheng,et al.Investigation of influence of design factors on the annular throat expansion-deflection nozzle performance based on the mixed orthogonal method[J].Journal of Solid Rocket Technology,2023,46(3):410-418.
[17]SUTHERLAND W.The viscosity of gases and molecular force[J].The London,Edinburgh,and Dublin Philosophical Magazine and Journal of Science,1893,36(223):507-531.
[18]阎超.计算流体力学方法及应用[M].北京:北京航空航天大学出版社,2006.YAN Chao.Methods and applications of computational fluid mechanics[M].Beijing:Beihang University Press,2006.
[19]MENTER F R.Two-equation eddy-viscosity turbulence models for engineering applications[J].AIAA Journal,1994,32(8):1598-1605.