崔胜楠,白桥栋,翁春生,孟豪龙,吴明亮,张世健,韩家祥
(南京理工大学 瞬态物理国家重点实验室,江苏 南京 210094)
连续旋转爆轰发动机(continuous rotating detonation engine,CRDE)是一种基于爆轰燃烧方式的新型发动机,通过在环形燃烧室中形成稳定的旋转爆轰波(rotating detonation wave,RDW)而产生推力。CRDE具有热效率高、结构简单、工况范围宽等优点,国内外针对气态燃料连续旋转爆轰发动机的点火起爆、旋转爆轰波的自持机理和传播模态等问题进行了大量研究。
煤油是CRDE工程应用的理想燃料。近年来,基于液态煤油为燃料的连续旋转爆轰发动机成为旋转爆轰发动机研究的一个热点,但以煤油/空气为推进剂的旋转爆轰发动机存在点火起爆困难等难题。作为煤油裂解产物最主要的成分,乙烯的化学性质较为活泼,是由氢气向煤油过渡研究的首要选择,但乙烯与氢气相比活性较低,爆轰波的稳定性相对较差。目前国内外对以乙烯为燃料的旋转爆轰发动机进行了较多实验研究。
HAN等采用乙烯/氧气为推进剂,ANAND等采用乙烯/空气为推进剂,研究了当量比、质量流量、喷注结构、背压等对CRDE的影响,发现只有在特定质量流量和当量比下连续旋转爆轰发动机才可稳定运行;发动机推力随质量流量的增加呈线性增加,爆轰波传播速度略有增加。KASAHARA等采用乙烯/氧气为推进剂开展了CRDE的火箭橇实验和地面推力测试实验。BURR等采用乙烯/氧气为推进剂研究了爆轰波在通过反应物射流阵列的基本特性,探究了乙烯/氧气燃烧室的流场特性。王宇辉等以甲烷/乙烯/空气为推进剂进行了连续旋转爆轰实验,研究了当量比对爆轰波传播速度和比冲等的影响。彭皓阳等采用乙烯为燃料,通过实验研究了燃烧室宽度、氧化剂活性、燃料活性、喷注方案和凹腔等对连续旋转爆轰波的影响。ZHONG等采用乙烯/乙炔/氢气为燃料开展了CRDE实验研究,结果表明,随着混合气中乙炔和氢气的增加,稳定工作范围扩大,爆轰波速度增大。
关于乙烯旋转爆轰传播特性的数值模拟研究目前还偏少。SATO等采用乙烯/空气为推进剂,使用氢稀释燃料,采用简化的8组分2步化学模型和21组分38步化学模型进行了详细的数值模拟,研究表明,加入氢气后经过弱爆燃到强爆燃之间的过渡,爆轰波稳定性增强,最终旋转爆轰波稳定传播。FAN等以乙烯/空气为推进剂,采用二维模型进行数值模拟,研究了喷嘴对乙烯旋转爆轰的影响。
从以上的研究可以看出,目前对乙烯燃料CRDE特性开展了大量实验研究,大多集中在乙烯/氧气为推进剂,采用乙烯/空气为推进剂进行实验研究也较少。数值模拟采用多组分反应机理的研究多集中于氢气/空气,乙烯/空气旋转爆轰采用多组分反应机理研究较少。
乙烯/空气连续旋转爆轰发动机工况范围、传播特性及其不稳定传播特性研究对于液体煤油CRDE的机理研究具有重大意义。本文采用详细反应机理对以乙烯/空气作为推进剂的旋转爆轰过程进行数值模拟研究,研究乙烯/空气旋转爆轰的传播特性和当量比等对其影响,研究结果可为乙烯和以煤油为燃料的CRDE实验研究提供理论指导。
在爆轰波的数值模拟中,不影响流场特征的情况下一般可忽略黏性、热扩散和热传导的影响,使用带化学反应源项的二维欧拉控制方程,其表达式为
(1)
(2)
(3)
(4)
(5)
由理想状态方程得:
=
(6)
采用乙烯/空气21组分(O、HO、N、CO、H、O、H、OH、CO、HO、HO、CH、CH、CHCHO、CHCO、CH、CHO、CH、CHO、CHO、CH)36步基元反应生成化学源项。反应速率常数采用Arrhenius公式计算:
(7)
式中:,和分别为指前因子、温度指数和活化能。
基于上述控制方程,本文使用开源CFD软件OpenFOAM中rhoReactingCentralFoam求解器对CRDE流场进行数值模拟。该求解器将可压缩流求解器rhoCentralFoam和反应流求解器rhoReactingFoam结合起来,对流项采用Kurganov的二阶Godunov型中心格式和迎风中心格式求解,化学源项采用欧拉隐式方法进行求解。孟豪龙等基于rhoReactingCentralFoam计算得到稳定爆轰波内流场,其准确性得到验证。
CRDE通常采用环形燃烧室,一端为燃料和氧化剂的进口,一端为爆轰产物的出口。当燃烧室周向长度远大于燃烧室宽度时可忽略宽度影响,将燃烧室沿母线展开得到二维旋转爆轰燃烧室模型,如图1所示。燃烧室模型的长和宽分别为300 mm和100 mm,下端为燃料和氧化剂入口边界,上端为爆轰产物出口边界,左右为周期性边界。
初始时刻燃烧室中填充条件如图1所示,在燃烧室左下角0≤≤20 mm,0≤≤40 mm处设置高温高压点火区域,为防止产生两道爆轰波,温度和压力以梯度增加,温度300 K≤≤3 000 K,压力0.1 MPa≤≤2 MPa;20 mm≤≤200 mm,0≤≤40 mm为预填充区,填充当量比为1的预混燃料,温度300 K,压力0.1 MPa;其余空间为空气。
图1 二维计算模型
计算域下端为入口边界,采用一维等熵入流边界条件,不考虑气流从集气腔进入燃烧室的流动损失。填充总压为,温度为,边界处计算压力为,入口边界分为3种情况:
①当<时,燃料不能进入燃烧室,按固壁边界处理。
②当<<时(为声速填充条件下的临界压强),此时按亚声速条件填充。
=
(8)
(9)
(10)
③当<时,预混的燃料按声速条件填充。
=
(11)
(12)
(13)
计算区域上端为出口边界,采用无反射边界条件,分为两种情况:当出口边界气相速度为超声速,出口边界状态二阶外推得到;当出口速度为亚声速时,出口压力等于环境压力。
考虑到计算精度、计算成本和反应机理的适用性,本文在爆轰管中对网格无关性和反应机理进行了验证。爆轰管直径20 mm,长300 mm,充满当量比为1的乙烯/空气混合气体,初始压力和温度分别为一个标准大气压和300 K。爆轰管左端设置高温高压点火区域,物理模型如图2所示。
图2 爆轰管计算模型
图3 不同网格下中心轴线处压力分布
为验证数值模型的准确性,在喷注总压0.6 MPa,喷注总温300 K,当量比1的条件下采用详细反应机理对旋转爆轰过程进行数值模拟。将乙烯/空气CRDE流场温度云图数值结果与BYKOVSKII实验中观测到的乙炔/氧气CRDE流场结构进行了对比,如图4所示。图中,为爆轰波,为斜激波,为滑移线,为爆轰产物膨胀区域,为新鲜燃料预混区域与燃烧产物接触断面,为阻塞段,为新鲜预混气。通过对比可知,本文计算结果与BYKOVSKII的实验结果所揭示内部流场结构定性一致。数值模拟结果计算得到爆轰波周向传播速度为1 772.0 m/s。利用CEA程序计算得到理论CJ速度为1 850.5 m/s,数值模拟结果与理论CJ速度相比爆轰波速度亏损为4.24%。
图4 模型验证对比图
本文数值模拟中喷注压力和总温分别为0.6 MPa和300 K,改变预混气各组分质量分数实现当量比从0.4~1.4逐渐增大,研究当量比对旋转爆轰波传播特性的影响,表1为计算得到的不同当量比下爆轰特性参数。由表1可看出,随着当量比的增加,爆轰波速度、压力增益以及燃料比冲均随之增大,且压力增益均保持在30%以上,速度亏损Δ逐渐减小。
表1 不同当量比下各工况爆轰特性参数
预混气体在喷注压力0.6 MPa、总温300 K和当量比1的条件下喷注进燃烧室,计算得到了稳定传播的旋转爆轰波。通过监测入口边界上的参数变化了解爆轰波的传播特性。图5给出了=100 mm,=0监测点处压力和温度随时间变化的曲线图,图6为爆轰波传播过程中压力和温度云图。由图5可知,=52 μs时由高温高压点火区域诱导形成的爆轰波第一次经过监测点;由图6(a)可知,=200 μs时由于燃料未能及时供应,一部分爆轰波与上一轮燃烧产物相交退化为斜激波,爆轰波和弱压力波经过碰撞均得到增强,继续向相反方向传播,阻塞燃料进入形成阻塞段;爆轰波扫过后逐渐形成稳定三角形新鲜燃料填充区域,燃烧产物从出口流出,爆轰波随后通过=100 mm,=0处监测点,出现图5中的第二个尖峰。由图5可知,在=222 μs后两条曲线呈现规律的周期性变化,当爆轰波扫过该监测点时温度和压力都迅速上升并存在明显的峰值点。爆轰波在燃烧室内连续稳定地传播,直至计算结束旋转爆轰波在燃烧室内共传播了11个周期,根据相邻尖峰之间的时间差计算得到稳定传播的爆轰波,平均传播速度为1 772.0 m/s。
图5 工况5下x=100 mm,y=0处压力和温度曲线图
图6(b)为爆轰波稳定传播时压力和温度云图。从图中可以观察到稳定的流场结构,由于波前预混气具有一定轴向喷注速度,爆轰波面和来流混合速度垂直,因此靠近入口壁面处的段略有倾斜,其波后峰值压力变化由预混气进入流场后进一步膨胀引起。由于斜激波后压力较低,爆轰波受到来自斜激波侧的一系列稀疏波的影响而衰减。段离斜激波较近,受影响程度较大,爆轰波强度明显减弱,爆轰波阵面落后于的延伸线,波面发生弯曲,呈现弧形,破坏了爆轰波原有的倾斜角度。
图6 工况5爆轰波传播过程中温度和压力云图
图7为工况5中=480 μs时刻化学组分云图,图7(a)为新鲜燃料CH分布图,图7(b)为中间产物OH浓度图。可看出燃烧产物与新鲜预混区存在明显接触断面;化学反应中间产物OH主要集中在爆轰波后方和斜激波下游,反映了化学反应阵面的形状。图8为=480 μs,=10 mm处沿轴方向部分组分曲线图。CH和O在经过新鲜燃料预混区域与燃烧产物接触断面时,燃烧室内反应物质量分数急剧下降,中间产物OH、CO和最终燃烧产物CO、HO含量升高;进入爆轰产物膨胀区域时,CH和O发生反应几乎转换为中间产物OH、CO和最终燃烧产物CO、HO,但O下降速度相对于CH较慢;随着反应进行,中间产物OH、CO和最终燃烧产物CO、HO含量趋于稳定。在爆轰波阵面处由于诱导区的存在,中间产物OH、CO含量小幅度升高后下降,燃烧产物CO、HO含量降低,CH和O随后进入燃烧室,含量升高。
图7 t=480 μs化学组分云图
图8 t=480 μs,y=10 mm处部分组分曲线图
图9和图10为=1 490 μs时温度、压力云图和曲线图。由图5黑色线框中=1 490 μs时温度和压力曲线图发现温度曲线突然升高,结合图9和图10 温度压力云图和曲线图可以看出,爆轰波扫过后波后反应产物温度和压力都较高,新鲜预混气三角区变形形成阻塞段,新鲜预混气未能及时供入,导致图5中线框温度曲线有所波动,对应压力曲线图也出现一小尖峰,当新鲜燃料供入后温度随之降低。
图9 t=1 490 μs压力和温度云图
图10 入口边界处t=1 490 μs压力和温度曲线图
图11为不同当量比下爆轰波传播速度及速度亏损曲线图。由图可知,随着当量比的增加,爆轰波传播速度和CJ速度呈增大趋势,增大趋势逐渐减小,当量比每增加0.1,平均增长率为3.98%;随着当量比的增加,速度亏损逐渐减小,减小趋势放缓,当量比每增加0.1,平均增长率为-23.15%。当量比为1.1时爆轰波传播速度为1 798.7 m/s,当量比下降到0.7时爆轰波传播速度为1 590.9 m/s,衰减了11.55%。在爆轰波成功起爆情况下,随当量比的增加,燃料活性增加,提高了爆轰燃烧释放的能量,爆轰波后产物的膨胀做工能力增加,导致爆轰波速度增加。
图11 爆轰波速度和速度亏损曲线图
当量比分别为0.7和0.8时压力和温度曲线图如图12和图13所示。由图12和13可知,当量比略低,爆轰波传播更加稳定,出现如图5中温度曲线剧烈波动的情况较少,这是因为当量比略低,爆轰燃烧化学反应剧烈程度不高,爆轰波后燃烧产物压力较低,膨胀做工能力减小,导致工质向下游传播能力下降,新鲜燃料可成功供入燃烧室,未形成图5温度曲线的波动。
图12 Er=0.7时压力和温度曲线图
图13 Er=0.8时压力和温度曲线图
由于CRDE工作过程中流场参数如压力和温度等变化剧烈,给出了进出口压力沿周向分布曲线图。图14和图15是=480 μs时刻不同当量比下,出口平面和进口平面压力沿周向分布情况。由图14可知,出口平面静压均大于外界反压()0.1 MPa。由图15和图6(b)可知,出口平面压力峰值是由斜激波引起,进口平面压力峰值是由爆轰波引起,沿爆轰波传播方向,出口平面压力峰值滞后于入口压力峰值;有部分区域进口静压高于给定喷注总压0.6 MPa,此段入口处于阻塞状态,没有新鲜气体喷入燃烧室。在稳定时刻计算了不同当量比下阻塞比(被阻塞长度与总入口长度之比)。图16为当量比与阻塞比曲线图。由图16可知,随着当量比的提高阻塞比逐渐增大,这是因为当量比的提高增加了化学反应剧烈程度,从图14和图15可以看出,随当量比提高爆轰波峰值压力和燃烧室内压力均有提高。当量比每增加0.1,出口压力峰值平均增长率为9.23%,进口压力峰值平均增长率为4.48%,阻塞比平均增长率为7.33%。
图14 出口压力曲线图
图15 进口压力曲线图
图16 阻塞比变化曲线图
为了比较当量比对燃烧室性能的影响,引入燃烧室压力增益和燃料比冲两个指标,分析不同当量比下燃烧室压力增益和燃料比冲变化,和表达式为
=(-)
(14)
(15)
式中:为出口截面平均总压,为入口总压,为出口密度,为出口轴向速度,-为出口与环境压差,为进口密度,进口轴向速度。
图17为进气总压为0.6 MPa时不同当量比条件下压力增益和燃料比冲变化曲线图。由图17可知,随着当量比的增大,压力增益和比冲呈增大趋势,压力增益保持在30%以上。分析认为,当量比的增大提高了化学反应剧烈程度,导致爆轰波的强度增强,爆轰波后的反射激波强度增加,燃烧室压力增益和燃料比冲随之增大。当量比每增加0.1,压力增益平均增长率为15.38%,燃料比冲平均增长率为7.27%。
图17 pt-g和If曲线图
本文对乙烯/空气的旋转爆轰过程进行了数值模拟,研究了当量比对旋转爆轰的影响和复杂反应机理在旋转爆轰反应计算中的适用性,分析了旋转爆轰的形成和传播过程,得出以下结论:
①采用乙烯/空气21组分36步基元反应可对旋转爆轰过程进行数值模拟,本文计算条件下旋转爆轰成功自持传播的当量比范围为0.7~1.1。
②在来流总压和总温一定的条件下,改变预混气的当量比,随着当量比的增加,化学反应剧烈程度增加,提高了燃烧室内温度,爆轰波传播速度、阻塞比、爆轰波峰值、燃烧室内压力和压力增益都随之增加。燃烧室内压力增益均保持在30%以上,比冲也随之增加。
③爆轰波稳定传播时,中间产物OH集中在爆轰波后方和斜激波下游,其反映了化学反应阵面的形状。传播过程中通过监测点检测温度,发现温度曲线出现波动,这是由于爆轰波扫过后压力过高,新鲜燃料无法供应引起的。
④出口平面压力峰值是由斜激波引起,进口平面压力峰值是由爆轰波引起,沿爆轰波传播方向,出口平面压力峰值滞后于入口压力峰值;有部分区域入口静压高于给定喷注总压0.6 MPa,形成阻塞段。