连续旋转爆轰发动机气液两相爆轰波传播特性二维数值研究

连续旋转爆轰发动机气液两相爆轰波传播特性二维数值研究李宝星;翁春生
【摘要】为了研究液体燃料连续旋转爆轰发动机( Continuous Rotating Detonation Engine,CRDE)中爆轰波形成与传播过程,采用二维CE/SE方法,对汽油、富氧空气两相连续旋转爆轰发动机爆轰过程进行数值模拟,分析了连续旋转爆轰发
动机气液两相爆轰流场和爆轰波结构及入口和出口处的流场变化规律,揭示了CRDE自持传播机理. 计算结果表明,燃料以时段阶梯填充方式来起爆旋转爆轰,可快速有效地形成单方向稳定传播的爆轰波;在周向方向上出口处的流场间断面要延后
于入口处的间断面,出口流场间断面主要是由斜激波和接触间断面造成的,而入口流
场间断面是由爆轰波引起的. 通过对气液两相CRDE的二维数值模拟,可更好地了解液体燃料CRDE的工作过程,为液体燃料CRDE研究提供指导.%In order to discuss the formation and propagation of detonation wave of the liquid fuel continuous rotating detonation engine( CRDE) ,the two-dimensional CE/SE method is used to simulate the detonation process of gasoline and oxygen-enriched air two-phase CRDE.The gas-liquid two-phase detonation flow field of CRDE,the structure of detonation wave,and the variation of flow field at inlet and exit were analyzed,then the self-sustaining mechanism of CRDE was revealed.The results show that rotating detona-tion is initiated by the method of multistep filling fuel,forming a stable detonation wave along one direction propagation rapidly and effectively.In circumferential direction,the discontinuity of exit flow field appears behind the discontinuity of inlet. Oblique shock wave and contact discontinuity
are the main reasons of the formation of discontinuity in the exit flow field,
while the discontinuity of inlet flow field is caused by detonation wave.The two-dimensional numerical simulation of gas-liquid two-phase CRDE gives us better understanding of the liquid fuel CRDE and provides guidance for the research of liquid fuel CRDE.
【期刊名称】《固体火箭技术》
【年(卷),期】2015(038)005
【总页数】7页(P646-652)
【关键词】CE/SE方法;连续旋转爆轰发动机;气液两相;流场;自持机理
【作者】李宝星;翁春生
【作者单位】南京理工大学瞬态物理国家重点实验室,南京 210094;南京理工大学瞬态物理国家重点实验室,南京 210094
【正文语种】中文
【中图分类】V435
连续旋转爆轰发动机是利用爆轰波在环形燃烧室内的一端连续旋转传播，爆轰产物从另一端开口处排出，爆轰产物在出口处高速排出时产生推力新概念发动机。

与脉冲爆轰发动机相比，连续旋转爆轰发动机具有以下优点：可连续性的填充高速燃料，相对脉冲间歇式，可填充更高的燃料流量，获得稳定较大的推力；只需要一次点火起爆，爆轰频率可高达几千Hz，而PDE每次都需要点火起爆，频率会受到限制。

CRDE是极具吸引力的新型推进系统。

关于连续爆轰的实验研究，俄罗斯的Bykovskii等[1-4]用不同的燃料(如氢气、煤油、丙酮等)，在不同大小以及不同形状的燃烧室内进行实验，形成旋转爆轰。

国
内的刘世杰、刘卫东等[5-6]利用H2和空气在环形燃烧室内做了不同模态下的连
续旋转爆震试验研究。

对连续旋转爆轰数值模拟方面的研究，法国Davidekol[7]
利用高精度的WENO和半隐的四阶龙格库塔的方法，对氢-氧爆轰进行了二维数
值模拟；Takayuki Yamada等[8]采用二阶精度的Harten-Yee non-MUSCL修正通量型的迎风TVD方法，对氢-氧爆轰进行了二维数值模拟。

国内的范宝春、姜孝海等[9-10]采用二阶半隐的龙格库塔法和五阶WENO格式，对连续爆轰进行了相关的二维及三维数值模拟。

北京大学的邵业涛、王健平等[11]利用五阶保单调MPWENO的格式，对连续旋转爆轰发动机进行了二维数值模拟。

国内外学者大多采用气体燃料对CRDE进行数值模拟，对于液体燃料的研究较少。

守恒元和求解元(the Method of Conservation Element and Solution Element，简称CE/SE方法)方法[12]被广泛用于脉冲爆轰数值模拟中。

目前，尚未发现
CE/SE方法在连续旋转爆轰数值模拟中的应用。

本文尝试应用CE/SE方法来模拟气液两相CRDE的二维爆轰现象，研究以液态汽
油为燃料，富氧空气为氧化剂的气液两相爆轰波特性；研究控制时间的阶梯填充方式对起爆过程的影响；分析连续爆轰波在环形燃烧室内的传播与循环过程，得到气液两相CRDE爆轰流场变化规律，揭示了气液两相连续旋转爆轰波的自持机理。

实际上连续旋转爆轰过程为三维爆轰过程，为了问题的简化，在此只考虑沿环形燃烧室的二维爆轰过程。

图1是由一个没有厚度的圆环，将圆环沿着一条母线ab剪开，得到矩形计算域，左边界ab和右边a′b′界是通过周期边界相连的。

下端为填充燃料端面，上端为排气端面。

在CRDE内气液两相的连续爆轰过程是非常复杂的。

为了简化计算，提出以下假设：
(1)旋转爆轰过程为二维的；
(2)气液两相爆轰过程为无粘过程；
(3)液滴为球形，并且温度均匀分布；
(4)液滴间互不影响；
(5)爆轰波经过液滴时，液滴仍保持球状，在气动力作用下发生剥离；
(6)液滴通过剥离蒸发成为气体，立即与富氧空气瞬间均匀混合。

根据上面基本假设，得到气液两相连续爆轰发动机的控制方程[12-15]：
式(1)中，各物理量含义解释如下：下标1表示气相，2表示液滴相；x，y分别为CRDE周向与轴向坐标；变量φ、ρ、u、v、E和p分别为体积分数、密度、周向速度、轴向速度、总能和气相压力，其中，φ1+φ2=1；Yk(k=1，2，3，4)依次
为气相中氧气、二氧化碳、水蒸气及汽油蒸气的质量分数，而氮气不参加任何反应，其质量分数可归一化处理；N为单位体积内液滴的个数；m21为液相转化为气相
的质量贡献率；分别为各组分的生成速率与消耗速率；Fx21，Fy21分别为气相与液滴群间的轴向作用力和径向作用力；Qcomb为化学反应释放的热量；Qconv
为气相与液滴群间的对流传热[13]。

燃料液滴剥离和蒸发对气相质量的贡献率m21计算为[16]
其中，r为燃料液滴半径。

燃料液滴半径变化率由气动剥离与蒸发两部分组成[16]：其中，通过剥离引起半径变化率为
式中μ、η分别为气体粘度和液滴粘度。

通过蒸发引起的半径变化率为[16]
式中λ为气体热传导系数；Nu为努赛尔(Nusselt)数；T为温度；L为燃料液滴的蒸发潜热。

燃料液滴剥离和蒸发过程中，液滴周围都是由燃料蒸汽包围，假设当液滴温度低于液滴沸点时，液滴只通过气动产生剥离，此时不考虑蒸发；只有液滴温度达到液滴沸点时，两者皆考虑。

气相与燃料液滴群间的对流传热Qconv[13]：
汽油的主要成分为辛烷，在这里采用辛烷的化学反应方程式替代汽油的燃烧反应过程。

辛烷一步总包反应方程式为
化学反应燃料消耗速率计算为
式中 A为化学反应指前因子；m，n为反应级数；Ea为活化能；Ru为普适气体常数。

其余化学反应组分的消耗或生成速率氧气二氧化碳)和水蒸气)，通过化学反应方程式(9)的当量系数进行计算。

CE/SE方法[12，17]是求解含强间断问题的一种新的计算方法，它最初是由NASA科学家Chang提出。

CE/SE方法的基本思想为将空间与时间作为同等变量，在整个空间与时间的计算区域内定义守恒元和求解元。

CE/SE方法与其他CFD计
算方法不同，将时间和空间统一进行处理，设立求解元和守恒元，保证了计算格式在整个计算域内满足物理上守恒；计算格式上简单、精度高、捕捉爆轰波等强间断能力强，不用黎曼分解，在计算时空间上的通量是不用方向分裂。

采用CE/SE方
法计算连续旋转爆轰波等强间断拥有独特的优势。

通过对守恒元每个边界面的时间-空间密度矢量的积分通量进行计算，得到的
CE/SE方法计算格式[12]：
2.1 源项的处理
在连续旋转爆轰过程中，爆轰波所引起的化学反应特征时间要远小于对流的特征时间。

因此，其源项可看作为刚性的，本文应用四阶龙格库塔法处理式(1)中的化学
反应源项[12-14]；先不考虑源项，用CE/SE方法对进行求解，然后将求解出的作为初值，求解常微分方程组dU/dt=R，最后得出四阶龙格库塔法中时间步长可取：式中ΔtR-K为四阶龙格库塔的时间步长；ΔtCE为CE/SE方法计算的时间步长，
一般取N=5～20。

2.2 初始条件及边界条件
在数值模拟过程中，周向长度L=300 mm、轴向长度H=100 mm的矩形区域，
计算网格数取300×100。

初始条件：图1左下角的红色区域1为点火区域，大小为20×30个网格数，点火条件为2.0 MPa和20×288.15 K以及气相的周向初速度为1 200 m/s；蓝色区域2为预混合好新鲜燃料；其他区域处为富氧空气。

边界条件：下端为入口边界，与进气管道相连的，总压为p0=0.6 MPa，总温为
T0=288.15 K。

假设气相、液相是以相同速度进入燃烧室，设燃烧室内边界临近
处的计算压力为p，进气边界分3种情况：
(1)当p≥p0，此时预混燃料不能进入燃烧室内，将入口边界按照固壁边界处理；
(2)当pcr<p<p0，此时为亚声速喷射。

入口边界各参数为
(3)当p<pcr，预混的燃料以声速喷射进入燃烧室，此时的边界参数为
式中 R为气体常数。

计算域的上边aa′为出口边界，使用无反射自由边界条件，分为2种：当出口边界为超声速时，出口边界状态根据二阶外推得到；当出口为亚声速时，出口压力等于环境压力(0.1 MPa)。

左右边界即为周期边界(左边物理参数与右边物理参数相同)。

3.1 阶梯填充方式下气液两相连续旋转爆轰波形成过程分析
试验中发现，CRDE燃烧室内会出现双波对撞现象，碰撞之后，可能会导致连续旋转爆轰波湮灭，不能形成稳定连续旋转爆轰。

为了点火之后，仅仅产生沿一个方向传播的爆轰波。

所以，在点火之后的第一个周期内，采取按时间进行阶梯填充新鲜的预混燃料(如图1所示)，来起爆并形成单方向传播的连续旋转爆轰波。

图2为点火之后不同时刻压力云图。

初始阶段在0≤x≤0.02 m、0≤y≤0.03 m该区域，填充标准状态下(压力0.1 MPa、温度288.15 K)气液两相预混燃料，液滴半径为50 μm，其他区域填充富氧空气
(氧气的质量分数为40%)。

t=0时刻，在0≤x≤0.0 2m、0≤y≤0.03 m区域，用
高温高压来模拟点火条件。

从图2看出，点火区域的右边有填充的新鲜燃料，点火之后，很快就形成向x+方向传播的爆轰波，在t=13.2 μs时刻，最高压力约为2.2 MPa，爆轰波在传播的
同时，前面不断填充新鲜燃料来维持其继续传播。

在t=153 μs时刻，爆轰波压力峰值已经到达了3.6 MPa，其相对应的温度峰值在2 100 K左右，爆轰波传播速度已经达到1 400 m/s，表明成功达到爆轰状态。

爆轰波在传播过程中，波后的压力会慢慢衰减，随着距离爆轰波波阵面越远，压力会降得越低，当压力低于填充总压时，即开始填充新鲜燃料。

在t=203 μs时刻，爆轰波传播到右端周期边界；同时，在区域已经填形成一定高度的新鲜燃料层，来维持爆轰波继续传播。

当t=307 μs时刻，形成了较为0≤x≤0.01 m稳定的爆轰波，波速保持在1 400
m/s，爆轰压力为4.5 MPa。

计算结果表明，点火之后，采取时段性阶梯填充燃料的方式，有利于在燃烧室内快速有效形成单方向的稳定爆轰波。

图3为周向上x=0.08 m的入口处压力和温度随时间变化曲线，从点火开始一共
循环了8个周期。

由压力和温度变化曲线可看出，爆轰波波阵面处高压间断面与
高温间断面是相互耦合在一起的，突显出了爆轰波的基本特征。

前2个周期内还
没有达到稳定状态，在第1个周期内波动较大，压力峰值只有3.5 MPa，相应的
温度峰值高达2 150 K，是受到高温高压点火条件的影响，从第3个周期开始压力峰值保持在4.6 MPa，温度峰值保持在2 100 K左右。

由图3可得，第3个压力
峰值和第3个温度峰值都出现t1=489 μs时刻，而第8个压力峰值和第8个温度峰值都是在t2=1 561 μs时刻，这段时间内一共循环了5个周期，通过计算公式
得到的爆轰频率f=4 664 Hz，通过计算公式f×L得到爆轰波波速为1 399.2 m/s，表明已经达到稳定爆轰状态。

3.2 气液两相连续旋转爆轰流场与爆轰参数分析
图4为t=1 218 μs时刻的稳定连续旋转爆轰的温度场和压力场云图，与Bykovkii[3]利用煤油-空气-氧气在环形燃烧室的爆轰试验中所揭示的连续旋转爆
轰流场定性的一致。

图4(a)中，1是爆轰波后的爆轰燃烧产物，2是连续旋转爆轰波，3是新鲜的预混燃料，4是斜激波，5是接触间断面，6是爆轰燃烧产物与新鲜预混燃料的接触面。

稳定的爆轰温度为2 100 K，爆轰压力为4.6 MPa。

图5为t=1218 μs时刻的出口与入口处的压力、温度、轴向速度以及周向速度沿
周向分布曲线。

从图5(a)可看出，进气总压为0.6 MPa，入口处的压力峰值为4.6 MPa，出现的
位置是在x=0.02 m处。

在0.1m≤x≤0.2 m这段范围内是大于进气总压的，没有
预混新鲜燃料填充进入燃烧室内部，随着离爆轰波的距离越远，压力会慢慢衰减，而在0.05 m≤x≤0.1 m入口压力值在总压值上下波动，该区域内部分压力已经减
小到比进气总压小，表明已经开始填充新鲜燃料。

从图5(a)中得知，其填充比(即
可填充燃料的长度与总入口长度之比)为0.58。

出口处的压力波动较大，该波动主要是由燃烧室内爆轰波头部的斜激波传到出口处引起的；沿着爆轰波的传播方向，出口处压力出现波动的位置要滞后于入口处的压力峰值出现的位置。

出口的压力值一直大于环境压力，由出口的边界条件判断，可得到出口处的流速均为超音速，这样有利于发动机出口处产生较大推力性能。

图5(b)中的入口处温度峰值为2 100 K，与压力峰值出在同一位置，体现出爆轰波高温高压断面相互耦合的基本特征，在0.05 m≤x≤0.2 m这段范围内温度都高于
进气总温，在0.05 m≤x≤0.1 m这段范围内温度处出现了较大波动，表明开始有
少量的混合燃料进入燃烧室内，进入的部分燃料的边界层会被之前的高温燃烧产物引燃，导致高低温波动较大；而其他区域温度较低，该区域是形成好的新鲜预混燃
料区，但在重新填充的新鲜燃料与高温的爆轰产物接触面处引起的燃烧，会影响新鲜燃料填充的同时，还会影响发动机出口的爆轰参数以及发动机的爆轰性能。

出口处的平均温度均在1 200 K以上，在斜激波处温度出现较大波动，主要也是由斜
激波膨胀引起的。

从图5(c)、(d)可看出，填充新鲜燃料时，气相与液相是以相同的速度进入燃烧室，速度最大值为200 m/s，该值出现在0.2 m处，正好也是入口压力和温度峰值位置；而在入口处压力大于进气总压的范围内进入速度都为0 m/s。

出口处的轴向气相速度均值都在600 m/s以上，而峰值可高达950 m/s，周向速度峰值只有380 m/s，两者的波动都较大，主要是受到斜激波的影响；而液相的轴向速度为零，说明出口处液相很早就反应完全。

而入口处的气相周向速度最大处出现在爆轰波处，随着距离爆轰波的越远，压力慢慢下降，当压力下降到与进气总压相等时，入口处的气相周向速度达到反向最大值，随着压力继续下降，反向的气相周向速度会渐渐趋向于零。

在周向速度从反向最大速度变为0 m/s这段范围内有出现一些波动，
这些波动主要是由已经进入燃烧室的燃料造成的，阻碍了稀疏波的膨胀。

在爆轰波波后区域的液相存在一定的周向速度，因为在爆轰波扫过燃料混合区域，液滴会同时发生剥离和蒸发，液滴在没有完全被剥离和蒸发之前，会在气相的作用力下，使液滴运动随着爆轰波周向运动，但液滴的周向速度要比气相的周向速度小很多；而波前形成稳定的新鲜燃料区域内气相和液相的周向速度都趋向于零，它们的运动规律基本相同。

(1)CE/SE方法能够很好地对连续旋转爆轰过程进行计算，有效地捕捉到爆轰波等
强断面。

(2)点火后，采用控制时间进行阶梯式的燃料填充方法，有利于连续旋转爆轰发动
机快速起爆，产生沿一个方向传播的爆轰波，进而达到稳定爆轰状态。

(3)分析达到稳定状态时连续旋转爆轰波结构，得出入口处和出口处流场变化规律，
入口处流场间断是由爆轰波引起的，随着爆轰波的传播会周期性变化；周向上出口处流场间断出现的位置要延后于爆轰波的位置，出口处流场参数的波动主要是受爆轰波波头产生的斜激波和接触断面的影响。

【相关文献】
[1] Bykovkii F A,Vedemikov E F.Continuous detonation of a subsonic flow of a propellant[J].Combustion,Explosion,and Shock Wave,2003,39(3):323-334.
[2] Bykovkii F A,Zhdan S A,Vedemikov E F.Continuous Spin Detonation[J].Journal of Propulsion and Power,2006,22(6):1204-1216.
[3] Bykovkii F A,Zhdan S A,Vedemikov E F.Continuous spin detonation of fuel-Air mixtures[J].Combustion,Explosion and Shock,2006,42(4):463-471.
[4] Bykovkii F A,Zhdan S A,Vedemikov E F.Continuous spin detonation of hydrogen-oxygen mixtures 2 combustion with an expanding annular
channel[J].Combustion,Explosion and Shock,2008,44(3):330-342.
[5] 刘世杰,刘卫东,林志勇,等.连续旋转爆震波传播过程研究(Ⅰ):同向传播模式[J].推进技
术,2014,35(1):138-144.
[6] 刘世杰,刘卫东,林志勇,等.连续旋转爆震波传播过程研究(Ⅱ):双波对撞传播模式[J].推进技
术,2014,35(2):269-275.
[7] Davidekol D M.Numerical simulation of H2/O2 continuous spin detonation with a detailed chemical mechanism[C]//21st ICDERS(CD-ROM),2007.
[8] Takayuki Yamada,Koichi A Hayashi.Numerical analysis of threshold of limit detonation in rotating detonation engine[C]//48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition.4-7 January 2010,Orlando,Florida.
[9] 姜孝海,范宝春,董刚.旋转爆轰流场的数值模拟[J].推进技术,2007,28(4):403-407.
[10] 张旭东,范宝春,归明月,等.旋转爆轰的三维结构和侧向稀疏波的影响[J].爆炸与冲
击,2010,30(4):37-340.
[11] 邵业涛,王健平.连续爆轰发动机的二维数值模拟研究[C]//第十三届全国激波与激波管会议,长沙:2008.
[12] 翁春生,王浩.计算内弹道学[M].北京:国防工业出版社,2006.
[13] 马丹花,翁春生.爆震管内扰流片对爆震波影响的数值分析[J].推进技术,2011,32(3):425-430.
[14] 马丹花,翁春生.脉冲爆轰发动机内三维两相爆轰的数值计算[J].推进技术,2010,31(4):503-507.
[15] Wang G,Zhang D,Liu K,et al.An improve CE/SE scheme for numerical simulation of gaseous and two-phase detonations[J].Computers and Fluids,2010,39(1):168-177.
[16] 洪滔,秦承森.气体-燃料液滴两相系统爆轰的数值模[J].爆炸与冲击,1999,19(4):335-342.
[17] Chang S C.A new approach for constructing highly sTablehigh order CSSE schemes[R].AIAA 2010-543.。