部分裂解煤油的旋转爆震发动机数值模拟

王 丹，周晨初，陈宏玉，严 宇，洪 流

(西安航天动力研究所液体火箭发动机技术重点实验室，陕西 西安 710100)

0 引言

旋转爆震发动机(rotating detonation engine，RDE)是一种利用爆震波在燃烧室内连续旋转传播，并产生稳定推力的发动机[1-5]。图 1为典型旋转爆震发动机的燃烧室结构和爆震波三维传播结构。

图1 旋转爆震发动机燃烧室结构Fig.1 Combustion chamber structure of RDE

由于旋转爆震发动机具有结构简单，在亚声速和超声速入流条件下均可工作，只需要一次起爆的优势，具备广阔的应用前景，美国、俄罗斯、法国、波兰和日本等国家已经开展了大量的研究工作[6-13]。在数值仿真方面，Zhdan等针对RDE结构特点对其流场进行了简化，采用两步化学反应模型进行了二维数值模拟[14]。Hishida采用加密网格对RDE开展了二维数值模拟，精确捕捉到了爆震波面上的三波点结构，数值胞格结果表明，受侧向膨胀影响，靠近入口处的胞格尺寸较小，而靠斜激波处的胞格尺寸较大[15]。Davidenk等采用氢氧基元反应模型对RDE流场结构进行了二维分析，并对比研究了燃烧室长度对发动机性能的影响[16]。Eude等针对RDE开展了三维数值模拟，研究结果表明受曲率影响，爆震波面后燃烧室内、外壁处压力不等，此压力梯度在流场内引起了横向扰动[17]。

国内方面，国防科技大学刘世杰等人采用二维数值模拟，验证了RDE的推力矢量调节能力，采用三维数值模拟，对比研究了燃烧室长度、宽度、尾喷管等几何参数对发动机性能的影响[18]。北京大学王健平等使用两步化学反应模型对RDE进行了二维数值模拟研究。数值计算获得了同轴圆管腔中间层曲面上连续爆震的多个循环过程，分析了燃料入射、提前燃烧、爆震波结构和波传播速度等几个关键问题[19]。南京理工大学马虎等[20]采用 9 组分 19 步的基元反应模型，对以氢气、氧气、氮气为反应混合物的RDE进行了二维模拟，研究结果表明，采用分段填充方法可以有效地形成沿同一方向传播的爆震波。

国内外采用氢气作为燃料的旋转爆震发动机燃烧室研究取得了众多进展，但是氢气存在贮存困难、价格昂贵、危险性高等缺点，采用常温且无毒的煤油作为燃料将极大地拓展旋转爆震发动机的应用领域，且能够实现低成本和长时间贮存。然而，采用液态煤油作为燃料存在起爆困难的问题，本文以煤油—空气作为推进剂的旋转爆震发动机为研究对象，提出一种采用预先加热模式使得液态煤油发生部分裂解，剩余被加热的煤油及裂解产物再与空气混合进行点火起爆的方案。通过裂解产生氢气、乙烯等高活性成分，降低煤油起爆的难度。采用数值模拟的手段开展部分裂解煤油旋转爆震发动机的仿真分析，分析流场内燃料组分分布对爆震波传播方向的影响，以及煤油裂解率对发动机性能的影响。

1 数值方法及计算模型

采用部分裂解方式实现点火起爆的煤油—空气旋转爆震发动机基本原理为:将发动机设计为两部分，第一部分为加热裂解部分，第二部分为爆震燃烧部分。其中第一部分通过加热的方式使得液态煤油发生汽化，其中小部分在高温作用下产生裂解，裂解过程将使大分子长碳链碳氢化合物变为小分子短链碳氢化合物，从而提高了燃料的化学活性。对煤油的加热方式有多种，其一可采用电加热模式，即通过将电加热丝盘绕于煤油管路之上，在传输过程中实现汽化和热裂解。其二可采用预先燃烧的方式，取全部煤油和少量空气进行富燃燃烧，以此过程作为煤油加热、裂解的途径。然而，采用预先燃烧方式所获得的裂解煤油组分与单纯电加热时组分不同，因此仅针对电加热模式的产物进行分析。

表1 化学反应组分Tab.1 Chemical reaction group

第二部分的爆震燃烧由第一部分所产生的部分裂解混合物与主燃烧室的空气掺混，通过起爆管点火起爆。由于有第一部分加热裂解过程的存在，进入第二部分的燃料活性明显提升，大幅降低了点火起爆的难度。与此同时，发生部分裂解的混合物与液态煤油相比，组分更加复杂，化学反应过程也更加复杂，需要通过数值模拟的方式着重对第二部分的爆震燃烧过程进行深入分析。将部分裂解煤油作为燃料边界条件，空气为氧化剂，模拟燃烧室的起爆过程。利用Fluent求解器，基于理想气体假设。采用隐式格式求解二维欧拉控制方程，忽略黏性、热传导和扩散等输运效应。用三阶MUSCL格式对对流项进行离散，时间项采用具有二阶精度的四步龙格—库塔法。计算采用压力基、非稳态计算，化学反应模型为有限速率模型，反应速率常数采用Arrhenius公式计算。

由于煤油成分中的长链化合物在高温下不稳定，很容易裂解为碳链更短的烷烃、烯烃等，详细的化学反应机理十分复杂，产物有上百种，NASA Langley提出了17组分，30步反应的简化机理[21]，在此机理的基础上，选取主要裂解产物：H2、CH4、C2H4、C3H6，其余为不参与化学反应的N2，忽略占比小于1%的其他产物。以上5种裂解产物加上被加热汽化的C12H23，共6种成分的混合物作为燃料，空气(21%O2+79%N2)作为氧化剂进行燃烧，化学反应采用多组分模型，各组分的化学反应均采用总包模型，化学反应模型来源于CHEMKIN软件[22]，组分列于表1，表中PEF(pre-exponential factor)为指前因子，AE(activation energy)为活化能。

计算针对的燃烧室为环形燃烧室，由于其厚度相对于其周长较小，计算中忽略其厚度的影响(见图1)，将其沿其母线拉直后得到如图 2所示的二维计算区域。计算域为尺寸为100 mm×260 mm 的矩形，对应图 1将圆柱沿x方向剪开，沿y方向拉直(圆柱圆形截面周长为260 mm，圆柱x方向长度L为100 mm)。

计算域的上下边界为周期性边界，模拟爆震波在周向不间断的旋转过程。计算域的左边界为入口边界(入口1和入口2)，为部分裂解煤油，划分为两个入口是为了分析初始爆震波的传播过程，实际上两入口的物理意义相同，都代表环腔头部环缝向燃烧室不间断供应的燃料(对应图 1左侧的环缝)，初始过程中，在入口区域初始化一个高温高压区域作为点火起爆源，该区域为0≤x≤10,40≤y≤50。计算域右边界为压力出口边界。

图2 计算模型Fig.2 Simulation model

2 网格划分及模型校验

采用ICEM软件绘制二维结构化网格，为了保证对爆震波的准确捕捉，估算爆震波波速在2 000 m/s附近，计算时间步长10-7s，每个时间步爆震波通过的距离为0.2 mm，因此将网格尺寸选择为0.1 mm，保证对爆震波的可靠捕捉，网格数量为260万。

为了验证计算模型的可靠性，进行典型爆震管起爆过程计算。爆震管初始结构为直径30 mm，长度600 mm，采用二维轴对称模型，取一半进行计算。初始计算时，在爆震管区域填充常温的、压力为0.1 MPa的煤油—空气当量混合比预混燃气，外界区域填充常温常压的空气。在爆震管入口侧设计一个矩形的点火起爆源，位置在x：0～5 mm，y:0～15 mm，点火能量为1 000 J，保证充足的点火能量使得发生爆震。起爆源在0时刻开始工作，持续时间10-5s。计算结果列于图3。仿真模型能够准确模拟点火起爆过程，模拟获得的爆震波波速与理论C—J速度计算值误差不超过15%，计算模型可靠。

表2 边界条件Tab.2 Boundary condition

图3 爆震管仿真Fig.3 Detonation tube simulation

3 起爆过程计算及分析

随着裂解率的升高，混合物中小分子烯烃的比例不断提升，化学活性增强，更加容易起爆。考虑到盘管电热丝加热能力以及煤油在管路中的停留时间，无法使得煤油发生充分的裂解，因此从工程设计角度考虑，着重针对裂解率20%以下的工况进行计算，研究爆震波的形成和发展规律，获得可靠的煤油爆震发动机起爆条件。裂解组分含量由CEA(chemical equilibrum with applications)软件计算获得。

表2列出了计算的边界条件，二维计算忽略燃料与氧化剂的掺混过程，进入燃烧室的工作介质是部分裂解煤油与空气完全掺混的混合气，入口1及入口2为速度入口，需要给定组分、温度及初始压力。根据煤油发生裂解的起始温度，选择850 K作为燃料混合物温度。速度边界通过自定义udf给定。计算方法为：当边界处压力p大于预混燃料喷射总压p0=0.35 MPa时，流动被阻塞，速度为0，临界喷射压力

当p0>p>pcr时，速度通过等熵膨胀计算，即

(2)对动漫行业典型人物进行访谈。典型人物访谈侧重对其个人成长与奋斗经历的了解与认识，并与自己的职业生涯规划相联系，撰写2000字左右的访谈报告。

当p

式中：R为预混燃料气体常数；T0为总温；γ为比热比。

上下边界采用周期性边界，模拟发动机的旋转过程。点火源初始化为温度为3 000 K，压力2 MPa。出口为空气，0.1 MPa压力边界，温度300 K。在旋转爆震试验过程中，为了保证爆震波单向传播，往往采用切向入射的方式。而在二维数值仿真计算中，难以模拟切向入射的过程，高温高压的点火区域将会向外产生球形爆震波，向上下两个方向传播，当两道爆震波相遇时，局部压力高于喷射压力，入口喷射速度为0，预混燃气无法进入燃烧室，此时缺乏维持爆震的预混燃料，可能导致火焰熄灭，因此要求在点火处产生的爆震波最好仅沿一个方向传播。因此在计算中通过分区初始化来模拟该过程，初始时刻在0≤x≤10，50≤y≤260区域填充850 K的预混燃料，其余部分则填充常温空气(模拟部分裂解煤油从头部环缝进入一段时间后充填于燃烧室头部区域的状态)，在0≤x≤10，0≤y≤40区域填充空气而不填充燃料，是为了人为控制爆震波向上传播(下方无燃料，仅能传播压力波而难以形成爆震波)，选择850 K是对裂解煤油的一个温度预估，为了对比裂解率的影响而将温度保持恒定，方便对比分析。图 4为计算初始化示意图。

图5 算例1起爆过程Fig.5 Initiation process of test 1

图4 初始条件Fig.4 Initial condition

首先针对算例1(煤油裂解率5%)进行分析。图 5为裂解率5%工况下压力及C12H23质量分数随时间的变化情况，其中C12H23的质量分数代表了主要燃料的占比。点火时间为0时刻，25 μs时，形成了正向的爆震波(y轴正方向为正向)，反方向仅有压力波。50 μs时，由于入口燃料的加入，反向也形成了一道爆震波，随后两个爆震波相遇对撞，对撞后的透射激波沿原方向继续传播。同时计算区域下方还产生了一个较弱的向上传播的压力波，150 μs时向上传播的压力波与向下传播的透射激波对撞，而此刻计算区域上半部分填入的燃料较多，二次对撞后透射的正向压力波得到了燃料的补充，迅速发展，形成了稳定的爆震波，而其余压力波由于未得到充分的燃料补充，逐渐耗散，最终在流场内只剩下了正向传播的一道爆震波，240 μs后该爆震波一直保持稳定自持传播，460 μs时经过调整过程和3个完整的旋转爆震周期，爆震波状态未发生明显改变。

为了监测爆震波的形成和发展规律，在计算域中取了7个点来分析该点温度、压力等参数的变化情况，各点的坐标如表3所示，图 2中示意画出了监控点P1～P7的位置。

对于爆震发动机的性能评估，最重要的参数就是发动机比冲和爆震频率，由于二维计算与实际发动机的情况略有差异，采用P1点压力的变化曲线来反映爆震频率，因为爆震波集中在入口侧传播，当调整阶段结束后可认为每次在P1点形成一个高压就意味着爆震波经过一次，即为一个周期。

表3 监控点坐标分布Tab.3 Coordinate distribution of monitoring points 单位:mm

图6和图7分别列出了入口附近沿y方向间隔50 mm均布的P1、P2、P3三点的压力及y方向速度的变化情况。从图中可以看出，不考虑200 μs前爆震波的形成和调整过程，仅分析之后较稳定的3个周期。可以看出，沿y轴方向，不同点的压力峰值有所差异，最低2 MPa，最高3 MPa以上，且P1压力小于P2，P2压力小于P3，具有一定的规律性。各点y方向的气流速度峰值基本保持一致，基本在1 000 m/s附近。可见，爆震波在环形燃烧室中传播过程中，流场内气流的传播速度基本保持不变，压力则呈现规律性变化，存在一个由弱变强再由强变弱的过程，可见爆震波在传播过程中受到入射预混气速度及压力的影响、波后燃烧产物分布的影响，压力并非时刻保持稳定，存在起伏现象。用P1、P2、P3三点出现压力峰的时间间隔来估算爆震波传播的平均速度。例如，P1和P2点y方向距离为50 mm，相继出现压力峰的时间间隔为25.5 μs，因此爆震波从P1点传播到P2点过程中的平均速度为距离与时间的比值，即1 940 m/s。经计算，不同周期不同监测位置爆震波平均传播速度基本相同。

图6 算例1监控点P1、P2、P3的压力变化Fig.6 Pressure history in monitor point P1,P2 and P3 of test 1

图8列出了沿x轴方向间隔20 mm均布的P4、P5、P6三点压力变化情况。从图2中可以看出，虽然P4、P5、P6与P1在同一条直线上，但是其压力峰值明显偏低，这是由于稳定后的旋转爆震波是由爆震波和斜激波组合而成的，具有侧向膨胀的性质，爆震波峰面宽度仅有20 mm左右。从稳定后的压力分布可以看出，P4压力大于P5，P5压力大于P6，可见随着向出口的发展，斜激波逐渐膨胀，压力降低。

图7 算例1监控点P1、P2 、P3的y方向速度变化Fig.7 Velocity history of y direction in monitor point P1,P2 and P3 of test 1

图8 算例1监控点P4、P5、P6的压力变化Fig.8 Pressure history in monitor point P4,P5 and P6 of test 1

4 裂解率对发动机性能的影响

以上分析均针对算例1(裂解率5%)，由于裂解率不同将影响发动机的燃料组分，为了分析组分对发动机性能的影响，增加算例2(裂解率10%)和算例3(裂解率20%)2种工况开展计算。3种工况仅燃料组分不同，其余计算条件相同。

算例2在50 μs时，两道面对面传播的压力波中间区域发生缓燃，图 9为此时刻的压力和组分云图，与图 5不同，裂解率升高后预混燃料活性更高，发生了期望外的缓燃，将消耗掉一定量的预混燃料，不利于稳定爆震波的快速形成。压力波对撞后形成了一些峰值较小的分散的压力波，在120 μs时甚至看不出存在明显的压力波，然而在170 μs时其中一个压力波得到了燃料的支持并不断增强，最终形成了正向传播的爆震波。

图9 算例2初始过程Fig.9 Initial process of test2

算例3在50 μs时，两道面对面传播的压力波中间区域同样发生缓燃，对撞后90 μs形成了一个比较明显向下传播的压力波，但是由于预混燃料并未填充进去，压力又衰减了。在120 μs到240 μs的时间段内流场里并没有明显的压力波，280 μs时产生了一个向下传播的压力波，随着燃料的不断填充，最终形成了稳定的反向传播的爆震波。图10为该工况稳定工作时的压力及组分云图。

图10 算例3稳定传播结果Fig.10 Stable propagation result of test 3

对比了不同组分的爆震情况可以看出，燃料组分对爆震波的形成和发展过程有非常重要的影响，3种工况除了燃料组分外其余条件完全相同，但是其传播方向有所差异，形成稳定爆震波的时间也不同。图 11列出了P1点压力变化情况。从图11中可以看出，算例1和算例2第一个稳定的爆震波形成于200 μs，而算例3的第一个稳定的爆震波则形成于400 μs，中间有很长一段时间并未形成爆震波(相当于稳定工作爆震周期的1.5倍)。这是由于裂解率提升后预混燃料活性提升，容易在燃烧室产生缓燃，消耗预混燃料，不利于稳定爆震波的形成，流场内的微小波系需要较长的调整时间，等待足够的预混燃料填入，方能形成稳定的爆震波，导致调整阶段时间延长。调整阶段过后3种工况的爆震波频率基本相同，压力峰值基本不变。

图11 不同裂解率工况P1的压力变化Fig.11 Pressure history in P1 under different crack rate

挑选一个周期来对比燃料组分对发动机比冲的影响，为了方便比较，将3种裂解率放在同一横坐标起点进行对比。图12列出了3种裂解率工况在稳定工作后某1个周期P1点的压力变化情况，从图中可以看出，3种裂解率情况下，1个周期的时间均为134 μs，旋转爆震波频率为7 500 Hz左右，且P1点的压力峰值均为2.3 MPa左右。

为了比较裂解率对性能的影响，通过作用在出口平面的流体性质来计算比冲。比冲由出口区域上产生的推力与流过的流量之比来确定，具体由下式计算

图12 一个周期内P1点压力变化情况(x轴平移为同一零时刻)Fig.12 Pressure history in P1 during one period (translate x axis to the same zero time)

计算中分子分母均在1个周期内对时间进行积分，可以获得3种工况在1个周期内的平均比冲(基于燃料、氧化剂及稀释气体混合物的比冲)。同时为了排除计算误差及爆震波传播稳定性的影响，计算了3个周期内的平均比冲，计算结果列于表4。

表4 比冲计算结果Tab.4 Impulse simulation result

可以看出，1个周期内，裂解率5%和10%时，比冲非常接近，然而裂解率20%时，比冲略有提升；3个周期内，不同裂解率下比冲基本相同。可见，裂解率对发动机稳定工况的平均比冲无明显影响。

5 结论

1)在爆震波形成的初始阶段，具有一定的不确定性，该不确定性与流场的初始状态有紧密的联系，最终稳定爆震波的传播方向取决于调整阶段燃烧场的波系变化及组分分布。

2)燃料组分对爆震波的形成和发展过程有非常重要的影响，裂解率提升导致预混燃料活性提升，容易在燃烧室产生缓燃从而消耗预混燃料，流场内的微小波系需要较长的调整时间等待足够的预混燃料填入，方能形成稳定的爆震波，导致形成稳定爆震波的时间延后。

3)调整阶段过后在不同裂解率预混燃料的工况下爆震波频率基本相同，压力及速度峰值也基本不变。

4)煤油裂解率20%以下时，改变裂解率对发动机稳定工作后的性能无明显影响。