连续旋转爆轰发动机外流场三维数值模拟①

邱玉杰，翁春生，武郁文，李宝星

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

0 引言

连续旋转爆轰发动机(Continuous Rotating Detonation Engine，CRDE)是一种利用爆轰波在燃烧室内沿周向高速旋转，产生的高温高压爆轰产物从出口排出形成推力的新型发动机。因具有结构简单、热循环效率高[1]、工作范围宽和推力矢量可调节等优点，CRDE近年来受到了广泛关注。

20世纪60年代，前苏联的Voitsekhovskii[2]首先开展了旋转爆轰的研究，通过圆盘实验装置实现了乙炔/氧气旋转爆轰的短暂传播。此后，Bykovskii等[3-4]采用CH4、H2、煤油和丙酮等多种气体、液体燃料，对不同燃烧室构型和尺寸大小的连续旋转爆轰发动机进行了系统的实验研究，并通过速度补偿技术观测到了旋转爆轰波的基本结构。Fotia等[5]研究了不同背压条件对爆轰波起爆过程的影响，并分析了不同来流条件和尾喷管构型对发动机推进性能的影响。Wolanski等[6]采用CH4/O2成功实现了旋转爆轰波的稳定自持传播，并在内径为140 mm、外径为150 mm的带有塞式喷管的发动机上，获得了250～300 N的推力。南京理工大学翁春生课题组[7]实现了液态燃料连续旋转爆轰波的成功起爆和稳定传播，分析了当量比对爆轰波传播速度的影响。林伟等[8]对单波和双波模态下的H2/Air连续旋转爆轰发动机进行了推力测试，详细分析了两种模态下爆轰波自持传播的典型波形特征和时域、频域特征。

由于爆轰波具有高温、高压和高速的特点，通过实验方法很难获得发动机流场内部详细的物理参数，需要通过数值计算揭示连续旋转爆轰复杂流场结构。2007年Zhdan等[9]首次实现了H2/O2连续旋转爆轰波平面内的二维数值模拟，获得的流场结构与实验结果定性一致。Yi等[10]对以H2/Air为燃料的连续旋转爆轰发动机进行了二维和三维数值模拟，研究了来流参数及燃烧室结构等对发动机爆轰流场及推进性能的影响，分析不同类型尾喷管对发动机出口参数及性能的影响。Eude等[11]利用网格自适应技术对三维连续旋转爆轰发动机燃烧室进行了数值模拟研究，将三维结果与二维进行对比，发现二者的流场结构和爆轰波传播特性基本相同。Schwer等[12]对燃烧室的外部流场进行了数值研究，分析了外流场对燃烧室内部流场的结构和性能影响。邵业涛等[13]通过三维数值模拟对实现连续爆轰的机理问题进行了分析，并探究了不同质量流量对发动机推进性能的影响。张旭东等[14]分析了H2/Air连续旋转爆轰波的三维结构和侧向膨胀波对爆轰波阵面的影响。李宝星等[15]采用CE/SE方法对气液两相连续旋转爆轰过程开展了数值研究，分析了气液两相爆轰波的流场结构以及燃烧室流场参数的变化规律。马虎等[16]对分开喷注下的连续旋转爆轰发动机进行了三维数值模拟，并对氧化剂和燃料的混合效果及三维连续旋转爆轰波的发展过程进行了分析。

总体而言，目前国内外学者对CRDE的研究主要集中在爆轰波传播特性和燃烧室内流场结构等方面，而对燃烧室外流场分布特性的探讨较少，出口流场的流动过程尚不清楚。为探究爆轰产物出燃烧室后外流场的分布规律及其对发动机推进性能的影响，本文对连续旋转爆轰发动机内流场及爆轰产物出燃烧室后的非定常排气过程进行了三维数值模拟。对爆轰产物出燃烧室后的流场发展过程以及充分发展羽流流场结构和主要参数的分布特性进行了分析，并对比分析了不同进气总压条件对流场结构及发动机性能的影响。

1 物理模型与计算方法

1.1 物理模型

连续旋转爆轰发动机通常采用柱状环形燃烧室，燃料和氧化剂从燃烧室进气端通过环缝或小孔进行喷注，经高能点火起爆后，会形成一个或多个爆轰波在燃烧室头部沿圆周方向传播，并燃烧波前燃料，高温高压爆轰产物从燃烧室尾部高速排出，从而产生推力。

连续旋转爆轰发动机内外流场计算域如图1所示，坐标系采用空间直角坐标系，分(a)、(b)、(c)计算域。其中，区域(a)为CRDE环形燃烧室，其内径为40 mm、外径为60 mm、长度为50 mm；区域(b)和(c)为外流场区域，分别为长200 mm×宽度200 mm×高30 mm、长200 mm×宽度200 mm×高210 mm的长方体。采用网格前处理软件GAMBIT进行网格划分，内外流场计算域内均采用六面体结构化网格，网格总数约为110万。

图1 计算域Fig.1 Computational domain

1.2 数值方法

采用FLUENT求解器，基于理想气体假设进行三维数值模拟。计算中，忽略粘性、扩散和热传导等输运效应，利用密度基显式求解器求解三维非稳态欧拉控制方程。采用三阶MUSCL格式对对流项进行离散，物理通量采用AUSM格式进行分解，时间项采用四步龙格-库塔法，选取恰当量比的H2/Air混合气体作为反应物，化学反应模型为一步反应的有限速率模型，反应速率常数由Arrhenius公式计算。

1.3 边界条件和初始条件

计算域的上边界为预混气入口边界，预混气进气总压为p0，总温T0=300 K，采用FLUENT的用户自定义函数(user-define function，UDF)进行设置，边界上每个网格单元的流动情况由该单元处的压力p决定，分为三种情况：

(1)当p≥p0，此时预混气不能进入燃烧室内，喷注速度为

v=0

(2)当pcr

(3) 当p≤pcr，预混气以声速喷注，即

燃烧室内外壁面为绝热固壁边界，燃烧室内柱端面处中心平面采用绝热固壁边界。外部流动区域上各出口平面为压力出口边界，分为两种情况：当出口为超声速时，边界点上所有流动参数由流场内部通过插值外推得到；当出口为亚声速时，边界点压力等于外界反压，而其他流动参数由内部流动外推得到，外界反压为0.1 MPa。由于化学反应主要在燃烧室内进行，将外流场区域(b)和(c)设置为无化学反应区域。

初始时刻，将爆轰波传播稳定的三维内流场的温度、压力、组分、密度、速度等参数插值到环形燃烧室区域(a)，对内流场进行初始化；外部流场区域(b)和(c)充满静止空气，温度为300 K，压力为0.1 MPa。

1.4 计算方法验证

1.4.1 网格无关性验证

本文分别采用0.5 mm×0.5 mm、1 mm×1 mm、2 mm×2 mm 三种网格，在三种网格条件下的计算结果如图2所示。可见，网格尺寸为1 mm×1 mm时能很好地捕捉到爆轰波间断面，进一步减小网格尺寸，对计算结果的影响不大，考虑到计算成本，在三维连续旋转爆轰发动机流场计算过程中采用尺度为1 mm的网格。

1.4.2 计算准确性

为验证该数学模型的可靠性，将二维CRDE爆轰流场参数的计算结果与理论值及文献[4]中的实验结果进行对比。图3(a)、(b)分别给出了计算获得二维流场温度分布云图与Bykovskii通过速度补偿法观测到的旋转爆轰波结构图，可看出，数值模拟的结果与实验定性一致。图3(c)给出了燃烧室入口附近某点压力和温度随时间变化曲线，将计算得到的爆轰参数和运用NASA CEA计算的同等条件下CJ理论值进行了对比，如表1所示。两者相对偏差都在5%以下，仿真结果与理论值吻合较好，表明了仿真结果的准确性。

图2 不同网格尺寸下入口处压力分布情况Fig.2 Distribution of pressure at inlet with different size of the grid

(a)计算所得流场静温分布 (b)实验扫描图片

(c)压力、温度曲线计算值

参数峰值压力pCJ/MPa峰值温度TCJ/K传播速度vCJ/(m/s)计算值1.241 02911.31950.0理论值1.303 32925.431961.9相对偏差/%4.80.50.6

2 结果与讨论

2.1 旋转爆轰三维内流场结构及参数分析

图4为预混气进气总压0.35 MPa条件下，爆轰波达到稳定传播状态时，燃烧室外壁面上温度和压力分布云图。其中，I为爆轰波，II为新鲜预混气体与上一循环爆轰产物间的接触面，III为与爆轰波相连的斜激波，IV为新产生的爆轰产物与上一循环爆轰产物间的接触面。接触面和爆轰波以及斜激波将流场分为4个区域：A区域为波前新鲜预混气体区，B区域为燃烧产物区，C区域为经斜激波压缩过的爆轰产物区，D区为爆轰波后的爆轰产物区。由图4可见，爆轰波阵面处是流场中温度和压力最高的区域，波头温度高达3000 K，波头处压力为3.1 MPa；爆轰波阵面前的新鲜预混气体区域，压力低于进气压力，燃料可由入口进入该区域，从而为旋转爆轰波的周向持续传播提供燃料。

图4 燃烧室外壁面上温度和压力分布云图Fig.4 Temperature and pressure distribution at outer wall of chamber

图5为爆轰波稳定传播时，燃烧室外壁面上监测点(x=30 mm,y=0 mm,z=5 mm)处的压力和温度随时间的变化曲线，其中共包含5个旋转周期。

图5 (x=30 mm,y=0 mm,z=5 mm)处压力和温度随时间变化曲线Fig.5 Temporal variation of pressure and temperature at a location(x=30 mm,y=0 mm,z=5 mm)

由图5可看出，压力和温度几乎同时到达峰值，反映了爆轰波阵面上化学反应界面和前导激波的耦合特性。压力峰值变化范围为2.99～3.17 MPa，温度峰值变化范围为2872.8～2908.3 K，说明此时爆轰波达到稳定传播状态，平均峰值压力(所有周期的峰值压力之和与总周期数的比值)为3.06 MPa，平均峰值温度为2884.0 K；相邻波峰之间的平均时间间隔为82.5 μs，由于燃烧室外径为60 mm，内径为40 mm，选取二者的平均值50 mm，获得爆轰波在燃烧室内的平均传播速度为1904 m/s，其工作频率高达12.1 kHz。爆轰波前混合气体的压力为0.22 MPa，温度为300 K，由CEA计算得到的C-J温度为2997.3 K，C-J压力为3.45 MPa，C-J速度为1980.5 m/s，数值计算所得的爆轰波温度、压力和速度与CEA计算的理论值的相对误差分别为4%、11%和4%。

2.2 外流场分析

2.2.1 外流场发展过程分析

图6为爆轰产物从燃烧室出口排出后，外流场x=0截面上不同时刻温度分布云图，描述了外流场的变化过程。初始时刻，高温、高压、高速的爆轰产物以射流的形式从环形燃烧室排出，由于失去了环形壁面的限制，爆轰产物迅速向四周膨胀，在燃烧室出口两侧附近形成球形高温射流区域，如图6(a)所示；在爆轰产物前方存在呈球形向静止空气中传播的激波阵面，不断压缩波前静止空气，图中爆轰产物与激波之间的高亮区域为受激波扫过的区域，该区域温度明显升高，上升到400 K左右，而激波前为未受扰动的空气，此处温度仍为300 K。至t=0.25 ms时刻，左右两侧的爆轰产物主要沿轴线方向，向下游流动，此时左右两股燃气射流并未交汇在一起，二者中间仍存在低温区，这主要是由于爆轰产物具有较高的轴向速度，而径向速度相对较低。t=0.8 ms时刻，左右两侧射流向中心膨胀，汇聚后共同沿中心轴线向下游流动，而在靠近中心壁面区域温度仍较低，约为500 K左右，这是由于中心壁面的存在，使此处流动受到限制，高温的爆轰产物未传播到此处。至t=1.8 ms时刻，爆轰产物已流出下部边界，外流场已达到充分发展状态，由图6(d)可看出，爆轰产物从喷管内排出后温度迅速下降，在出口附近区域流场温度约为1800 K，这是因为高度欠膨胀、超声速的爆轰产物出燃烧室后由于失去了环形壁面的限制迅速膨胀，在出口附近区域温度迅速降低；羽流中心区域爆轰产物的温度明显高于外围爆轰产物的温度，温度高达2700 K，这是由于出口爆轰产物速度比周围气体高，在引射作用下，在中心轴线附近形成低压区域，将周围高温爆轰产物卷吸至此处，导致此处流场温度升高；在流场下游区域，随着距喷管出口截面距离的增加，流场温度呈下降趋势，在外场下边界附近，流场温度下降到1000 K左右。

(a)t=0.16 ms (b)t=0.25 ms

(c)t=0.8 ms (d)t=1.8 ms

2.2.2 充分发展外流场结构及参数分布特性分析

图7为1.8 ms时刻，流场达到稳定状态时，环形燃烧室出口下游z=51 mm截面上的压力和温度分布云图。由图7可看出，与环形燃烧室对应的圆环区域处压力和温度较高，其最高压力出现在斜激波波面处，高达0.49 MPa；流场中存在呈螺旋状向外扩散的高压区域，这是由于燃烧室下游周期性传播的斜激波，传播至出口进入外流场后，没有内外壁面的限制，沿径向方向快速向周围大气中传播导致的。此外，可看出，羽流中心区域的压力较低，由于燃烧室结构的特殊性，中心圆形壁面附近区域无射流入射，爆轰产物在出环形燃烧室后沿径向方向，向周围大气中膨胀的同时，也会向射流中心膨胀，导致中心圆形壁面附近区域压力迅速下降。由图7(b)中可看出，羽流中心存在呈圆形的高温区域，最高温度达到2700 K，对应压力云图可知，由于该区域的压力低于周围区域压力，导致高温燃气不断流到此处，从而形成高温区域。

图8给出了t=1.8 ms时刻，外流场x=0平面上马赫数、压力分布云图。由图8(a)可看出，高温高压的爆轰产物以欠膨胀的方式从燃烧室内排出，使出口附近流场速度迅速增加，最大速度增加到Ma=3.2；而在距燃烧室出口截面约26 mm处马赫数发生突变，甚至出现了Ma<1的区域，对比压力云图可看出，在燃烧室出口附近，燃气射流的压力在很短的距离内下降到环境压力以下，为了协调出口附近膨胀流场和下游流场的压差，使其满足等压条件，在距燃烧室出口截面约26 mm处出现了一道斜激波，跨过激波后压力和温度迅速上升，而速度迅速下降。由图8(b)可看出，羽流外围存在呈圆弧状高压区域，对应于图7(a)可知，这是由周期性循环的斜激波出燃烧室后，同时沿径向和轴向传播，在周围大气中形成螺旋状的高压区域。在羽流流场下游区域，流场速度有逐渐减小的趋势，在外流场下边界附近羽流的马赫数下降到1.4左右。

(a)静压分布云图 (b) 温度分布云图

2.2.3 沿轴向方向流场压力分布

为了进一步研究爆轰产物出燃烧室后在出口附近的流动情况，在y=0平面(x=25 mm，y=0直线)上沿z轴方向设置6个监测点，各点相对位置如图9所示，分别为燃烧室内点1(25,0,30)、出口平面上点2(25,0,50)、出口近场点3(25,0,55)、点4(25,0,60)、点5(25,0,80)、点6(25,0,100)(单位为mm)。

(a)Ma分布云图 (b)压力分布云图

图10给出了1.364～1.884 ms时间内各监测点的压力时程曲线。

图9 y=0平面上监测点位置Fig.9 Monitoring point on the section y=0

(a)点1 (b)点2

(c)点3 (d)点4

(e)点5 (f)点6

由图10可看出，各监测点处的压力波形呈周期性变化，压力峰值总体稳定，说明此时流场达到较为稳定的状态。各监测点处压力峰值之间的平均时间间隔相同，约为0.082 ms，与爆轰波在燃烧室内的传播周期相同。由图10(a)～(d)可看出，点1～4各点压力变化规律基本一致，斜激波经过各点时，压力迅速上升到最大值，各点处的压力峰值分别为0.60、0.46、0.33、0.24 MPa，表明燃烧室内各点压力明显高于外流场中各点的压力，且随着距离燃烧室出口平面距离的增加，出口附近各监测点平均压力峰值逐渐减小；斜激波过后，爆轰产物开始迅速膨胀，各点压力开始迅速下降达到达极小值，随后在接触间断的作用下压力又细微上升，到达到极大值后开始缓慢下降，直到下一次斜激波到达后又开始上升，呈现周期性变化。5、6两点处的压力仍呈周期性变化，但变化规律并不十分明显，随着距出口平面距离的增加，各点压力峰值呈缓慢下降趋势，由0.21 MPa下降到0.16 MPa。

2.2.4 沿径向方向外流场压力分布

为研究出口附近流场压力沿径向方向分布规律，在y=0，z=60 mm直线上设置6个监测点，各点的编号和位置如图9所示，分别为点7(0,0,60)、点8(10,0,60)、点9(35,0,60)、点10(45,0,60)、点11(60,0,60)、点12(90,0,60)(单位mm)。图11给出了1.424～1.924 ms时间内各点压力变化曲线。

(a)点7 (b)点7

(c)点8 (d)点10

(e)点11 (f)点12

由图11可见，各点对应的压力曲线均呈周期性变化，点7的压力值在0.032～0.077 MPa之间变化，点8处的压力变化范围是0.010～0.082 MPa，二者的压力始终为负压(压力小于环境压力)，出现负压主要因素是出燃烧室的爆轰产物快速沿径向，向射流中心膨胀，导致压力快速降低到负压；中心轴线上点7的压力峰值间的时间间隔为28.5 μs，约为爆轰波传播周期的1/3，这主要是因为中心轴线上各点处的压力，同时受燃烧室出口处各方向上斜激波的影响，导致此处压力振荡频率较高。

9、10两点处压力曲线变化基本一致，斜激波扫过后压力迅速衰减到负压，随后在接触间断的作用下压力又略有上升后迅速下降，直到下一次循环的斜激波经过出现新的压力峰值；随着距中心轴线距离的增加，压力峰值略有下降，由0.21 MPa变化到0.14 MPa，表示爆轰产物沿径向方向，向流场外围膨胀，导致压力下降。11、12两点处的压力曲线呈现出单调的下降和上升，这是由于在流场外围区域的空气，受到斜激波的扰动，各点压力呈周期性单调变化；随着距中心轴线距离的增加，各点压力值的波动逐渐减小，压力峰值逐渐降低，由0.124 MPa下降到0.118 MPa，压力极小值逐渐上升，由0.079 MPa上升到0.087 MPa，二者同时接近于环境压力。

2.3 不同进气总压对流场分布及发动机性能的影响

在其他参数保持不变的情况下，研究不同进气总压力0.15、0.35、0.55 MPa条件对连旋转爆轰发动机流场特性及推进性能的影响。

图12给出了不同进气总压条件下，爆轰流场达到稳定状态时，某一时刻内流场以及外流场出口近场y=0平面上的压力分布云图及流线分布图，此时斜激波传播至图中燃烧室左侧出口处。总体上来看，不同进气条件下各平面上压力分布规律类似，爆轰产物出燃烧室后迅速膨胀，出口下游流场静压迅速下降，均在出口中心平面附近形成负压区域；在流场下游靠近中心轴线附近存在激波使流场静压回升；当进气总压升高时，出口处爆轰产物压力升高，出口下游流场的膨胀更剧烈，激波距燃烧室出口距离增加。由图12中流线的分布可看出，燃烧室内爆轰产物几乎都是沿轴向方向运动；而在燃烧室出口处，内壁面附近的爆轰产物沿轴向方向运动的同时，也会沿径向方向，向压力较低的射流中心区域流动，在轴线附近交汇在一起，共同沿轴向方向，向流场下游流动。由于出口左侧斜激波处的爆轰产物压力高于另一侧的爆轰产物压力，在中心壁面附近，左侧爆轰产物沿径向方向膨胀更为强烈，使得此处爆轰产物径向方向速度较大。在进气压力为0.15 MPa时，在中心平面附近区域形成了回流区域，随着进气总压的升高，斜激波处爆轰产物沿径向膨胀作用增强，径向方向流动速度增加，回流区域逐渐消失。由图12(c)可看出，当进气总压升高到0.55 MPa时，由于爆轰产物压力较高，使得交汇处爆轰产物压力高于环境压力，交汇后的爆轰产物向下游膨胀，压力迅速下降到环境压力以下，并在z=112 mm处形成一道激波。出口处爆轰产物沿径向方向流动将会导致轴向速度损失，给发动机推力带来了不利影响。

(a)p0=0.15 MPa (b)p0=0.35 MPa (c)p0=0.55 MPa

为进一步分析爆轰产物沿中心平面膨胀对发动机推进性能的影响。图13给出了不同进气总压条件下，燃烧室出口中心平面上等压线分布图，此时出口面上斜激波传播至y轴正方向处，并沿逆时针方向传播。由图13可见，各平面上压力分布规律基本一致，平面右侧靠近斜激波处壁面上的压力最高，沿y轴负方向压力呈递减趋势，在平面左侧边缘区域附近压力最低，且中心平面上大部分区域压力均在环境压力以下。这是由于在燃烧室出口处，爆轰产物压力在周向上分布不均，斜激波处压力最大。爆轰产物出燃烧室后在沿轴向方向流动的同时，也会由斜激波一侧沿中心平面向压力低的一侧膨胀，导致在中心平面上的压力迅速下降到环境压力以下。当进气总压升高时，燃烧室出口处斜激波压力增加，导致爆轰产物在中心平面附近膨胀作用增强，各平面上的压力最小值呈下降趋势。

(a)p0=0.15 MPa (b)p0=0.35 MPa (c)p0=0.55 MPa

发动机的推力F主要由两部分组成：一是气体流动产生的推力F1；二是中心平面所受的推力F2。各推力计算公式：

(1)

(2)

F=F1+F2

(3)

式中ρ和vz分别为出口处爆轰产物密度和轴向速度；p为爆轰产物压力；pb为出口环境压力；pc为中心平面上的压力；exit为出口平面；center为中心平面。

计算得到发动机推力与进气总压力的关系曲线如图14所示。由气体流动产生的推力F1随着进气总压的增大而增大，增大趋势较为明显；但中心平面上产生的推力F2始终为负值，给发动机的推力带来负的增益。

图14 发动机推力与进气总压关系曲线Fig.14 Relationship diagram between engine thrust and inlet stagnation pressure

总体上来看，在其他条件不变的情况下，进气压力越高，发动机推力越大，与进气总压力几乎成线性关系。当进气总压为0.55 MPa时，发动机的推力达到了1160 N。在后续研究中，可以通过在燃烧室尾部安装中心锥来改善出口附近流场分布，消除出口附近羽流中心的负压区域，实现发动机的推进性能增益。

3 结论

(1)本文所用的数值方法能够真实模拟连续旋转爆轰波三维流场，计算所得结果与实验结果基本一致。爆轰波阵面处是内流场中温度和压力最高的区域，进气总压为0.35 MPa时，波头温度高达3000 K，波头处压力为3.1 MPa。

(2)爆轰产物排出燃烧室后静压迅速降低到环境压力以下，在流场下游产生激波使静压上升，随着进气总压增大，激波距燃烧室出口距离增加。羽流外围的空气受出口处斜激波的扰动，压力呈周期性变化，随着距出口距离增加，压力值的波动减小，压力峰值逐渐接近于环境压力。

(3)引射作用使得在出口附近羽流中心形成低压高温区域；中心平面上产生了负的推力，降低了发动机的推进性能。

(4)在给定的模型下，当出口环境压力不变时，发动机推力随进气总压的升高而近似呈线性增加。当进气总压为0.55 MPa时，发动机推力达到了1160 N。