推进剂质量流量对三维非预混氢氧旋转爆轰波传播模态的影响

康朝辉,翁春生,武郁文,徐 高,丁陈伟,雷 特,邱彦铭

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

燃烧分为爆燃和爆轰两种模式,其中爆燃是一种等压燃烧,而爆轰近似为等容燃烧。研究表明,爆轰能够在短时间内释放出巨大能量,与等压燃烧模式相比,理论热循环效率最高可提升近50%[1]。爆轰推进技术可应用在航空航天领域,其中,连续旋转爆轰发动机(continuous rotating detonation engine,CRDE)由于结构简单紧凑、仅需单次点火起爆和能适用于从低速到高超声速来流范围的优势,成为了国内外学者研究的热点[2]。CRDE不仅适用于火箭基和冲压基,还能与涡轮等发动机组合起来使用,将来在飞机、导弹及临近空间飞行器等领域都有较好的发展前景。

关于CRDE的研究最早可追溯到上世界五六十年代,俄罗斯的VOITSEKHOVSKII[3]首次在盘式燃烧室中利用乙烯/氧气的推进剂组合实现了连续旋转爆轰,并研究了旋转爆轰波的波系结构。此后,研究人员针对不同燃料和氧化剂组合的推进剂开展了大量的研究。固体燃料发生爆轰所需条件较为苛刻,但是DUNN等[4]通过实验发现使用纳米级炭黑颗粒可以在空气中成功实现旋转爆轰。在液体燃料方面,在航空发动机领域使用最为广泛的是煤油,然而气液两相旋转爆轰过程存在燃料雾化和两相掺混等诸多难点,实验研究中主要通过在燃料中掺入微量氢气[5]、将富氧空气作为氧化剂使用[6]以及提高推进剂的喷注温度[7]等方法实现旋转爆轰波的触发。

在气态燃料方面,LIU等[8]在小尺寸燃烧室内实现了甲烷/空气的推进剂组合的旋转爆轰波稳定传播,但工作范围较窄。PENG等[9]分别使用甲烷、乙烯、氢气与空气的推进剂组合开展了CRDE实验。实验结果表明,旋转爆轰波随着当量比的增加呈现出由不稳定传播向稳定自持传播变化的趋势,其中,仅在氢气/空气的推进剂组合实验中出现了双波模态,且氢气/空气推进剂组合的工作可调范围比乙烯和甲烷更宽。MA等[10]使用氢气/空气推进剂组合开展了CRDE实验。实验结果表明,在燃烧室内同时存在爆燃和爆轰现象,并在实验中首次分析了单-双-单波的跃迁现象。李宝星等[11]开展氢气/空气CRDE实验研究,发现随推进剂总质量流量的增加,旋转爆轰波会发生传播方向的改变以及单波向多波转变等过程。刘世杰等[12]对氢气/空气推进剂在大范围质量流量条件下开展了实验研究,发现模态转变所需工作条件存在极限边界,在边界范围之内,工作条件的调节并不会改变旋转爆轰波的传播模态。空气作为氧化剂时,由于空气中的氮气作为惰性气体降低了化学反应速率,并且在相同的燃料流量条件下,空气作为氧化剂比氧气的流速更大,所以旋转爆轰波的稳定性更高。

以氧气为氧化剂时,化学反应活性的急剧提升会导致复杂的旋转爆轰波的行为模式。LIN等[13]对甲烷/氧气的推进剂组合进行了CRDE实验研究,发现在较低的质量流量下旋转爆轰波无法维持稳定传播,而在高流量条件下发现了旋转爆轰波的多波传播模态。由于氢气与纯氧燃烧获得的比冲比烃类燃料高约45%[14],将氢气/氧气的推进剂组合应用于火箭基连续旋转爆轰发动机中,能进一步提升火箭发动机的性能。SOSA等[15]通过在氢气中引入示踪剂,成功捕捉到了五波同向的工作模态,证明了氢气/氧气的推进剂组合在火箭基连续旋转爆轰发动机中的应用潜力。

由于在纯氧环境下氢气的反应速率快,可能对实验设备造成毁坏,现阶段主要开展氢氧旋转爆轰数值仿真研究。如GAILLARD等[16]通过限制喷注段的化学反应,对非预混氢氧旋转爆轰进行准三维仿真,详细地探究了氢氧旋转爆轰波的流场结构。孙健等[17]在此种方法的基础上对氢气/氧气预混气进行了二维数值模拟,发现由于氢气/氧气的推进剂组合活性较高,爆轰波传播的不稳定性更高,且传播速度和频率远高于氢气/空气组合的CRDE。FAN等[18]通过对氢氧CRDE的实验和仿真研究发现,在氢气/氧气的推进剂组合的CRDE燃烧室内极易形成多波模态,且爆轰波波头数目的增加会引起爆轰波波速的损失。

针对CRDE的仿真模拟研究通常采用预混喷注模型,即忽略燃料和氧化剂在喷注腔内的掺混过程,燃料和氧化剂提前以均匀混合的状态直接进入燃烧室。目前已有SWIDERSKI[19]、WU等[20]分别针对不同燃料和空气的推进剂组合进行了研究分析,揭示了旋转爆轰波的传播特性。但在实验条件下,为了防止回火,CRDE大多采用非预混喷注结构,燃料和氧化剂分别喷注进入混合腔,在燃烧室内边掺混边燃烧。孙健等[21]对非预混喷注的氢气/空气推进剂组合CRDE进行数值模拟,发现空气喷注喉道变大时,需要提高推进剂总质量流量来维持旋转爆轰波的持续传播。马虎等[22]针对氢气/空气的推进剂组合,研究了喷注压力等条件对非预混CRDE爆轰性能的影响,并阐述了旋转爆轰波向上游喷注扩张段传播的行为模式。考虑到氢气/氧气的推进剂组合化学反应极其剧烈,STECHMANN发现采用氢气/氧气预混气喷注的设置方式时,燃烧室内更倾向于形成爆燃燃烧[23]。

从国内外研究现状可以看出,目前的研究结果大多局限于预混喷注或二维流场,无法掌握发动机内真实流动条件下氢氧旋转爆轰波传播特性及传播模态,气动参数等对氢氧旋转爆轰波传播机理的影响尚未明晰。因此,本文拟针对三维非预混小孔-环缝喷注构型连续旋转爆轰发动机,采用详细化学反应机理,开展氢氧旋转爆轰流场数值模拟,获得旋转爆轰波稳定传播时燃烧室内流场结构,并在此基础上探索推进剂总质量流量对氢氧旋转爆轰波传播模态的影响,本文研究结果可为阐明氢氧旋转爆轰波传播机理提供一定的理论基础。

1 计算方法与物理模型

1.1 数值方法

在旋转爆轰数值模拟中,相对于对流项,耗散项的影响较小[24],为了兼顾计算效率,在本文计算中忽略黏性耗散项、热传导和扩散等输运效应的影响,基于理想气体假设,使用密度基求解器来求解三维带化学反应的多组分气体Euler方程。

在三维笛卡尔坐标系下的微分守恒形式的多组分气体化学非平衡流Euler方程的表达式为[25]

式中:U为守恒变量;F,G,H为无黏对流通量;S为化学反应源项。各符号的表达式为

U=(ρρuρvρwEρwi)T
F=(ρuρu2+pρuvρuw(E+p)uρwiu)T
G=(ρvρuvρv2+pρvw(E+p)vρwiv)T
H=(ρwρuwρvwρw2+p(E+p)wρwiw)T
S=(0 0 0 0 0ωi)T

式中:ρ为混合气体密度;u,v,w分别为3个方向速度;p为混合气体压力,wi(i=1,2,…,N-1)为i组分的质量分数;N为总组分数;E为混合气体单位体积的总能。

式中:h为混合气体单位质量焓;各组分的质量反应生成率ωi由基元反应的动力学模型给出:

式中:Mi为第i种组分的摩尔质量;βij,αij分别为第j个基元反应中生成物及反应物的化学计量数;Rf,j,Rb,j分别为摩尔浓度表示的第j个基元反应的正反应和逆反应的速率。

1.2 计算模型

图1给出了简化后的连续旋转爆轰发动机三维模型,燃烧室的内径为66 mm,外径为76 mm,计算域总长度为50 mm,喷注结构采用小孔环缝构型,其中氧气喷注环缝的喉道宽度为1 mm,90个氢气喷注孔均匀分布在拉法尔喷管扩张段中部,喷注孔直径为1 mm。计算域采用结构化网格划分,燃烧室主流区保持统一的网格精度,轴向方向上网格尺度为0.4 mm,径向方向上网格尺度为0.25 mm,周向方向上网格尺度为0.25 mm,在喷注段附近进行局部网格加密,氢气喷注孔附近及氧气喷注环缝喉部附近最小网格尺度均为0.05 mm,计算域的网格数量约为300万。

图1 计算域及网格加密Fig.1 Computational domain

氢气喷孔和氧气环缝的进口采用质量流量入口条件,出口采用压力出口边界,分两种情况:①当出口为超声速时,所有守恒变量由内部流场外推得到;②当出口为亚声速时,边界点压力等于外部背压,而其他守恒变量由内部流场外推得到,外部背压设置为0。采用AUSM失通量分裂法分解物理通量,化学反应有限速率模型使用Arrhenius9组分21步的基元反应模型[26]。

对于非预混进气方式下的旋转爆轰点火起爆过程,首先仿真获得稳定的氢氧冷态流场,然后仅保留燃烧室头部部分区域的氢氧推进剂,其余区域填充不参与反应的N2,并在燃烧室头部设置高温高压的点火区域,诱导燃烧室内形成沿单向传播的初始爆轰波。

1.3 模型验证

1.3.1 求解方法验证

为验证该求解方法的可行性,首先采用该求解方法计算二维旋转爆轰波(rotating detonation wave,RDW)结构,计算结果如图2所示。图2(a)中仿真计算所得旋转爆轰波结构与图中BYKOVSKII等[27]实验观测到的结构基本吻合;图2(b)为仿真计算中燃烧室入口压力及温度监测曲线,由计算可得爆轰波波速为1 827.0 m/s,与在同等静温260 K静压0.08 MPa条件下通过CEA软件[5]计算出来的CJ理论值1 965.1 m/s的误差为7.03%。因此可认为该计算方法能实现旋转爆轰波的计算。

图2 仿真计算结果与实验结果对比Fig.2 Comparison between calculated results and experimental results

1.3.2 网格无关性验证

当氢气质量流量为10 g/s,氧气质量流量为80 g/s,总温为300 K时,对网格的无关性进行检验。对该小孔-环缝喷注结构燃烧室计算域分别采用最大网格尺度为0.2 mm、0.4 mm和0.6 mm的网格进行计算,在喷注段小孔和环缝附近进行网格局部加密,燃烧室主流区保持统一的网格精度,最终计算域的网格量分别为469万、300万和263万。表1为不同网格分辨率下计算所得平均爆轰波(detonation wave,DW)波速以及与CJ理论值相比的速度损失,当网格分辨率大于300万时,继续提高网格分辨率对计算结果影响较小,可认为300万网格分辨率已满足计算需要。因此,本文所有算例均采用该网格分辨率进行。

表1 不同网格分辨率下爆轰波的速度损失Table 1 Velocity deficit of different grid resolutions

2 结果分析

保持当量比为化学恰当比,开展了氢氧旋转爆轰波传播特性研究,分析推进总质量流量对爆轰波传播模态的影响,计算工况如表2所示。

表2 计算工况Table 2 Calculation conditions

图3给出了工况1下950 μs时刻,燃烧室中心切面的温度二维展开云图。其中,旋转爆轰波沿圆周方向传播,爆轰产物经过膨胀沿轴向排出。在爆轰波波后附近,氢气和氧气的喷注受到阻塞,新鲜燃料无法喷入;在爆轰波波后较远的位置,新鲜燃料进入燃烧室,逐渐形成三角形的新鲜燃料填充区。其中,由于氢气和氧气化学反应的活性非常高,在燃烧室入口附近出现了氢气的预燃烧现象,在该区域部分氢气和氧气刚接触就被燃烧消耗掉。另外,爆轰波和下游斜激波之间的温度间断为滑移线,此线为两次循环中爆轰产物的接触面。爆轰波、斜激波、滑移线的交汇结构保持了流场的稳定性。

图3 温度分布云图(工况1,t=950 μs)Fig.3 Temperature distribution in Case 1 at 950 μs

如图4所示,在拉法尔喷管的扩张段上,爆轰波向喉道延伸并向波后方向弯曲,在喷管扩张段壁面的约束下形成反射激波,该反射激波在壁面作用以及氢气喷注出流的影响下发生变形。此外,在旋转爆轰波的传播过程中,集气腔内会形成一道上游斜激波,与旋转爆轰波同向传播。

图4 壁面压力分布云图(工况1,t=950 μs)Fig.4 Pressure distribution on combustion chamber wall in Case 1 at 950 μs

2.1 RDW稳定单波传播模态

如图5所示,在工况1下,点火起爆950 μs后,燃烧室内形成了稳定单波传播11～12个周期的RDW,爆轰压力为1.1 MPa,爆轰波波速为2 596.1 m/s,相较同等条件下CJ理论值的波速损失为8.48%,这是在非预混喷注条件下燃料和氧化剂的不均匀掺混导致的。从燃烧室入口处监测点的压力曲线局部放大图可见,图4中的反射激波在爆轰波波峰之后形成了第二个尖峰。

图5 燃烧室入口压力监测及局部放大图(工况1,t=0～950 μs)Fig.5 Pressure at combustor inlet in Case 1 from 0 to 950 μs and enlarged view

如图6所示,点火起爆的初始阶段,产生了沿顺时针方向传播的爆轰波(DW1),而在反方向上由于缺乏推进剂的支持,仅有一道激波(shock wave,SW)1沿逆时针方向传播。DW1再次传播至初始点火区域附近时与SW1发生对撞,此时新鲜推进剂已填充足够高度,使DW1在接触上一轮的高温燃烧产物后仍能持续传播。而SW1由于缺乏新鲜推进剂支持,并未能增强为新的爆轰波。最终燃烧室中形成了稳定单波传播模态。

图6 压力及温度分布云图(工况1,t=10～300 μs)Fig.6 Pressure and temperature distribution in Case 1 from 10 to 300 μs

2.2 RDW稳定双波同向传播模态

工况2及工况4下,经过一段时间燃烧室内形成了稳定的双波同向模态。如图7所示,在工况2中,燃烧室形成了沿顺时针方向传播的爆轰波同向双波模态,平均爆轰压力为0.9 MPa,单个爆轰波平均传播速度为2 458.3 m/s,与CJ理论值相比波速损失13.6%。

图7 燃烧室入口压力监测(工况2,t=0～2 700 μs)Fig.7 Pressure at combustor inlet in Case 2 from 0 to 2 700 μs

如图8所示,1 200 μs时燃烧室内DWA和DWC沿顺时方向同向传播。在1 286 μs时,燃烧室内发生了再起爆,在DWC前形成了沿顺时针方向传播的DWB,而沿反方向传播的激波与DWC对撞后,由于新鲜燃料及氧化剂被消耗殆尽(如图9所示),爆轰波前无法形成三角形新鲜混气区域,导致爆轰波无法维持传播,进而衰减并最终湮灭,燃烧室内又再度形成爆轰波双波同向的传播模态。此后,燃烧室内均维持爆轰波双波同向传播模态,但由图明显可见两爆轰波的强度并不相等。

图9 H2和O2组分分布云图(工况2,t=1 296 μs)Fig.9 H2 and O2 distribution in Case 2 at 1 296 μs

实际上该RDW稳定双波同向传播模态中,两爆轰波强度处于此消彼长的状态。为分析这种强弱交替现象,通过爆轰波波速vDW来代表爆轰波强度的相对大小。如图10所示,取50 μs时间内爆轰波的平均速度代表爆轰波波速。在爆轰波实现稳定传播后,DWA的波速在一段时间内大于DWB,而在下一时间段内DWA的波速又小于DWB。其中,两爆轰波速度差值最大为515 m/s。产生这种现象的原因可能为:在同向双波模态下,当爆轰波A传播速度较快,爆轰波B与爆轰波A之间的距离变长,使得爆轰波B的波前新鲜混气区域变大,从而增大了爆轰波B的强度;而爆轰波A的波前新鲜混气区域变小,反过来抑制了爆轰波A的强度。

图10 DWA和DWB波速分布(工况2,t=2 000～2 600 μs)Fig.10 Speed of DWA and DWB in Case 2 from 2 000 to 2 600 μs

在工况4下,推进剂总质量流量提高至157.5 g/s,相比于工况2,此时氢气和氧气的喷注速度较大。如图11所示,燃烧室内形成了沿逆时针方向传播的RDW稳定同向双波模态,爆轰压力为0.8 MPa,单个爆轰波平均传播速度为2 303.0 m/s,相较同等条件下CJ理论值的波速损失为19.9%。

图11 燃烧室入口压力监测(工况4,t=0～2 500 μs)Fig.11 Pressure at combustor inlet in Case 4 from 0 to 2 500 μs

如图12所示,燃烧室经历复杂的非稳态过程后,在1 300 μs时可以观察到燃烧室内存在一道沿逆时针方向传播的爆轰波DWB和数道向不同方向传播的激波,其中激波SWA经过发展,增强为爆轰波DWA,形成了两道较为稳定的逆时针方向传播的爆轰波DWA和DWB,而其他激波逐渐消散。最终,在1 570 μs时刻的云图中燃烧室内的弱激波均已消散,仅有两道稳定自持的旋转爆轰波DWA和DWB,且两个爆轰波的强度接近。

图12 压力及温度分布云图(工况4,t=1 300～1 750 μs)Fig.12 Pressure and temperature distribution in Case 4 from 1 300 to 1 750 μs

图13给出了工况4下爆轰波实现稳定传播后的波速随时间变化图,DWA和DWB的波速变化趋势相同,爆轰波波速的差值最大仅为37 m/s。

图13 DWA和DWB波速分布(工况4,t=1 650～2 050 μs)Fig.13 Speed of DWA and DWB in Case 4 from 1 650 to 2 050 μs

结合图11的压力曲线可知,在爆轰波实现稳定传播后,稳定双波同向传播模态中的两个爆轰波爆轰压力基本相等,爆轰波波速变化趋势一致,波速相近,说明该工况下两爆轰波的强度基本相等。在工况4更大的质量流量条件下,燃烧室内爆轰波波速变化幅度更小,说明该同向双波传播模态更加稳定。

2.3 RDW稳定四波同向传播模态

在工况5及工况6下,在燃烧室内均形成了4个同向传播的旋转爆轰波。以工况6为例,如图14所示,此时燃烧室内形成了4个沿逆时针方向传播的同向爆轰波,爆轰压力为0.70 MPa。从壁面压力云图中可见,爆轰波波头位置靠近喉道且明显向氧气喷注面倾斜,燃烧室内部存在伴随激波跟随爆轰波传播。

图14 壁面压力分布云图(工况6,t=2 440 μs)Fig.14 Pressure distribution on combustion chamber wall in Case 6 at 2 440 μs

如图15所示,在压力曲线图中可见明显的压力尖峰,以及在尖峰之后较低的压力脉动,即该伴随激波扫过测点时引起的压力脉动。在四波同向模态下,单个旋转爆轰波的平均传播速度为1 722.9 m/s,相较同等条件下CJ理论值的波速损失为40.5%。

图15 燃烧室入口压力监测(工况6,t=1 600～2 600 μs)Fig.15 Pressure at combustor inlet in Case 6 from 1 600 to 2 600 μs

图16为2 380～2 440 μs时刻内爆轰流场结构。由图可见4个沿逆时针方向传播的爆轰波DW1、DW2、DW3、DW4,爆轰波波头之间相位差基本相等。

图16 压力及温度分布云图(工况6,t=2 380～2 440 μs)Fig.16 Pressure and temperature distribution in Case 6 from 2 380 to 2 440 μs

如表3所示,随着推进剂质量流量的增加,爆轰波波数增加时,爆轰波波速损失增加。工况2和工况4在同一RDW传播模态下,随推进剂质量流量增加,爆轰波波速损失也会增加。这是由于推进剂流量的增加,喷注速度加快,导致爆轰波向喷注面倾斜,从而引起爆轰波波速在圆周方向上分量减小,即爆轰波波速损失增加。

表3 仿真计算结果Table 3 Calculation results

2.4 RDW不稳定传播模态

在工况3下,推进剂总质量流量为135 g/s,从图17可以看出,经过较长时间燃烧室内仍不能形成稳定的爆轰波传播模态,爆轰波压力和波速等仍在不断变化。实际上,在该工况下燃烧室内爆轰波处于单波和三波对撞交替的不稳定传播模态。

如图18所示,在1 750 μs时,燃烧室内存在一个沿顺时针方向传播的爆轰波DW1和一对向相反方向传播的爆轰波DW2、DW3。在1 770 μs时DW1与DW2发生对撞,由于在透射激波SW1和SW2的波前区域的新鲜燃料及氧化剂被消耗殆尽,爆轰波前无法形成三角形新鲜混气区域,导致爆轰波无法维持传播,进而衰减并最终湮灭。由于氢气和氧气的化学反应活性较高,经过一段时间会发生局部再起爆,又形成一对向相反方向传播的爆轰波DW4、DW5。在1 830 μs时DW3会与DW4产生新一轮的对撞过程,由于同样的原因使得爆轰波再次处于单波传播模态。

图18 压力及温度分布云图(工况3,t=1 750～1 840 μs)Fig.18 Pressure and temperature distribution in Case 3 from 1 750 to 1 840 μs

为探究不同推进剂质量流量下形成的RDW传播模态的稳定性,对各工况下爆轰波峰值压力进行统计,并进一步计算平均压力及其标准差以衡量旋转爆轰波传播稳定性。

如图19所示,在推进剂质量流量为90 g/s的工况下RDW平均压力最高,标准差最小,而在推进剂质量流量为135 g/s工况下RDW的平均压力最低,标准差最大。RDW传播模态稳定性呈先下降后上升的趋势,这是由于在推进剂质量流量增加的过程中,RDW传播模态首先由稳定单波转向双波同向传播模态,然后由双波同向传播模态转向多波对撞的不稳定传播模态,又转变成双波同向传播模态,最后转化成四波同向传播模态。其中,在多波对撞传播模态下爆轰波的稳定性最低。

图19 各工况平均压力及标准差分布Fig.19 Average pressure and standard deviation for various cases

3 结束语

本文以氢气为燃料,氧气为氧化剂在小孔/环缝喷注结构下对三维非预混旋转爆轰燃烧室进行数值模拟,在不同工况下获得多种RDW传播模态,对旋转爆轰流场进行分析,得到如下结论:

①当推进剂总质量流量为90 g/s时,在燃烧室内形成了稳定单波传播模态;当推进剂总质量流量为112.5 g/s和157.5 g/s时,均在燃烧室内形成了稳定双波同向传播模态,其中在112.5 g/s的入口条件下形成的两个爆轰波是强弱交替的,而推进剂总质量流量为157.5 g/s时两个爆轰波的强度接近;当推进剂总质量流量为180 g/s和270 g/s时,燃烧室内存在爆轰波稳定四波同向传播模态;推进剂总质量流量为135 g/s时,燃烧室内爆轰波处于单波和三波对撞交替的不稳定传播模态。

②小孔-环缝喷注结构下,爆轰波在传播过程中,会在燃烧室入口扩张段弯折并形成反射激波。此外,在旋转爆轰波的传播过程中,在爆轰波压力作用下,集气腔内会形成一道激波,与旋转爆轰波同向传播。随着推进剂质量流量的增加,爆轰波波数增多,会引起爆轰波波速损失增加。爆轰波波数不改变时,推进剂质量流量的增加会导致爆轰波向喷注端面的倾斜程度增加,使爆轰波速度沿圆周方向速度分量减少,引起速度亏损增加。

③通过对各工况平均压力及标准差的分析,发现工况1中获得的RDW稳定单波传播模态下RDW平均压力最高,工况3中不稳定传播模态的RDW平均压力最低。结果表明,随着推进剂质量流量增加,RDW传播模态稳定性呈先下降后上升的趋势,其中,在多波对撞传播模态下爆轰波的稳定性最低。