旋转爆轰发动机内煤粉-氢气两相爆轰流场数值研究

翁春生 倪晓冬 续晗

摘 要：旋转爆轰发动机因其较高的热循环效率在航空航天领域备受关注。 为了研究气固混合燃料在空气氛围中的旋转爆轰特性， 在圆柱坐标系中建立了气固混合相燃料爆轰理论模型， 并基于三维守恒元与求解元方法， 对圆盘形燃烧室内的爆轰过程进行三维数值模拟。 计算获得了气固混合相旋转爆轰波稳定传播时的流场结构， 分析了爆轰流场的热力学参数分布特征、 化学反应区分布以及旋转爆轰波后的波系分布特点。 研究结果表明， 微米级煤粉在氢气辅助作用下可快速反应并支持旋转爆轰波的稳定传播， 但因盘形燃烧室的扁平结构特征和非预混喷注方式， 导致爆轰波波阵面呈现为不规则的曲面。 煤粉和氢气射流穿透空气层的能力差异导致气固混合燃料的化学反应区发生分离， 最终以双反应区的爆轰组织形式支持旋转爆轰波稳定传播。 旋转爆轰波不规则的曲面结构使其在扁平的燃烧室中发生多次反射， 进而在爆轰波后有规律地出现多道反射激波。

关键词：  气固混合相爆轰； 爆轰波； 粉末燃料； 煤粉； 非预混； 数值模拟

中图分类号：TJ760； V231

文献标识码： A

文章編号：  1673-5048（2024）02-0060-11

DOI： 10.12132/ISSN.1673-5048.2024.0036

0 引  言

旋转爆轰发动机（Rotating detonation engine， RDE）是一种利用爆轰波在环形［1-2］、 空桶形［3-5］或圆盘形［6-10］等具有圆柱面约束的燃烧室中连续旋转传播所产生的高温高压气体产物由发动机出口高速排出， 进而产生推力的新概念发动机。 与定压燃烧相比， 爆轰过程接近于定容燃烧， 因此RDE相较于传统发动机具有更高的热循环效率。 此外， RDE凭借其一次点火便可连续稳定运行的工作模式， 以及推力稳定、 推重比大、 宽范围的流入速度等优势［11］受到世界各国的广泛关注。 在军机、 导弹、 超高声速飞行器等领域具有广阔的应用前景。

目前， 旋转爆轰推进技术的开发应用主要以气体和液体燃料为主， 而随着军事应用场景的不断扩大， 以含能材料等固体推进剂为燃料的旋转爆轰发动机的研制被提上日程。 相较于气体和液体燃料， 固体燃料具有储存方便、 能量密度高、 抗冲击性能好等优点。 将固体燃料制成粉末状应用于RDE， 使其继承固体燃料优点的同时， 又兼备了气体、 液体燃料燃烧效率高、 流量可调等优点。 因此， 近几年作为航天及军事强国的俄罗斯和美国， 先后开展了固体粉末旋转爆轰发动机的实验研究。 俄罗斯的Bykovskii等人［12-16］率先在不同规格大小的圆盘形燃烧室中对各类型煤粉的旋转爆轰过程进行了相关实验研究， 在不严格控制燃料流量的情况下验证了各类型煤粉末旋转爆轰的可行性。 美国的Dunn等人［17-19］随后对煤粉在环形燃烧室中的旋转爆轰过程进行了实验研究， 并通过理论及数值研究进一步分析了煤粉爆轰的可行性。 国内方面， 与Dunn等人同期， Xu等人［6-8］在控制燃料流量稳定的情况下于直径150 mm的盘形燃烧室中对煤粉、 铝粉等燃料的爆轰特性进行了实验研究， 并对爆轰波波头高度、 强度、 传播速度以及发动机推力和燃料比冲进行了系统的分析。

粉末旋转爆轰发动机作为近几年刚开始研究的课题， 研究初期， 由于固体燃料的旋转爆轰特性尚不明确， 许多关键技术需要探索与验证。 因此， 选用一种取材方便、 成本较低、 安全可靠的实验材料非常重要。 煤炭， 作为生活中常见的固体燃料， 与旋转爆轰中经常用到的气态甲烷、 液态煤油等同为碳氢化石燃料，  在世界范围内分布广、 储量大、 开采成本低， 且燃烧产物以无害的CO2和H2O为主， 因此成为现阶段粉末旋转爆轰发动机的理想燃料。 但煤的化学反应活性较低， 研究其旋转爆轰特性需借助于氢气等活性气体的辅助。 因此， 上述研究工作中关于煤粉的爆轰特性研究均是在氢气等活性气体辅助的情况下进行的。 然而， 因为研究时间较短， 且实验测试手段的局限， 煤粉-氢气-空气旋转爆轰燃烧室中大量的物理参数无法观测和记录， 因此对于燃烧室中的流场参数分布及变化规律缺乏直观的认识和理解， 使得粉末旋转爆轰的研究留下巨大的空白。 为此， 需借助于数值模拟的方式进一步展开研究。 Zhu等人［20-21］对环形燃烧室中煤粉、 铝粉的旋转爆轰过程进行了二维模拟研究， 围绕爆轰波前颗粒三角区与空气三角区不完全重叠的现象进行了一系列研究。 然而， 现实中的旋转爆轰燃烧室均为三维非预混喷注结构， 煤粉-氢气-空气旋转爆轰波的组织结构也呈现出三维的结构特点。 为了更好地了解煤粉-氢气-空气旋转爆轰波的传播特性， 且对前期实验研究工作［6， 8］进行补充， 本文建立了三维煤粉-氢气-空气混合相爆轰理论模型， 并基于时空守恒元与求解元（The Space-Time Conservation Element and Solution Element， CE/SE）方法［22-27］， 模拟了盘形燃烧室中非预混喷注条件下煤粉-氢气-空气旋转爆轰波的传播过程， 进一步分析了煤粉氢气混合燃料爆轰波的三维结构特征、 化学反应区分布特征以及热力学参数、 波系结构在不同方向上的分布特点。

1 混合相爆轰理论模型与计算方法

1.1 混合相爆轰理论模型

实践中， 圆盘形RDE与离心式压气机和径流涡轮表现出良好的匹配特性， 同时也具有便于可视化观测和布置监测点的优势， 因此对处于初始研究阶段的粉末旋转爆轰， 其多数实验研究工作是在盘形燃烧室中进行的。 为了衔接并进一步完善前期研究工作， 本文的数值研究工作也围绕盘形燃烧室开展。 为此， 在考虑了圆盘形燃烧室主要结构特征的基础上， 建立了扁平的煤粉-氢气-空气非预混圆盘形计算域， 如图1所示。 计算域上下端面分别为燃烧室上壁面和下壁面； 外圆柱面为喷注面板， 喷注面板上阵列分布有燃料及空气喷注入口， 位于上方的入口阵列为燃料入口， 位于下方的入口阵列为空气入口， 其余部分为固壁（为便于描述， 喷注面板的固壁部分称为喷注面）； 内圆柱面为出口。 该计算域出口半径为R1、 喷注面半径为R2、 厚度为L。

1.1.1 基本假设

由于煤粉-氢气-空气旋转爆轰过程既包含气相燃料的爆轰， 也包含固相粉末燃料的爆轰， 整个过程极其复杂， 为了方便研究与计算， 需对该过程进行简化。

在开展粉末爆轰特性研究的工作中， 通常以微米或亚微米级的超细粉为研究对象。  由于颗粒尺寸极小且在粉末发动机流量供给条件下颗粒数目极大， 因此， 本文数值模型主要用于描述大量超细固体颗粒在气体中的宏观运动以及气固两相之间的相互作用过程。 旋转爆轰发动机工作时， 固相燃料颗粒群与气相燃料充分混合后经小孔喷入燃烧室， 整个过程颗粒群悬浮于气相环境， 且颗粒数量足够大， 因此在宏观上可以使用连续体假设来将固相颗粒群近似为“准流体”， 而单个颗粒的各物理特征则均匀一致。 数值计算时， 气相与固相颗粒群可被处理为共享空间区域的双流体， 且用流体动力学方程来描述这两相的运动， 而质量、 动量与能量的源项则被用来描述两相流体之间的热-机械-化学相互作用。 由于固相的真实密度远大于气相， RDE中固相燃料的体积分数始終不超过总体的1%， 因此， 固相体积的变化对流场的影响十分有限。 此外， 爆轰过程是瞬态变化过程， 爆轰波的传播速度通常达到几千米/秒， 因此燃烧室流场与壁面的传热作用对爆轰波的传播影响较小， 分析爆轰流场特征时可忽略其与壁面的热交换过程。

目前， 大量的旋转爆轰发动机数值研究工作基于Euler方程开展。 Bruce等人［28］指出， 粘性对非定常的激波运动过程几乎没有影响。 林玲等人［26］和Oran等人［29］则通过对比Euler方程和Navier-Stokes方程计算得到的爆轰流场， 发现粘性对爆轰波的传播及结构特征也没有明显影响。 因此， Euler方程完全满足对RDE中旋转爆轰波传播过程及流场结构特征研究的需求。

基于上述分析和论证， 本文在使用Euler方程的基础上总结出如下几条假设： （1）爆轰过程中燃烧室壁面绝热； （2）粉末燃料视为固体颗粒群， 颗粒及其反应后产生的气体在每个网格内与气体氧化剂瞬间混合均匀； （3）颗粒在整个爆轰过程中均保持球状， 且固体颗粒间互不影响， 颗粒内部温度均匀分布； （4）忽略气固混合物中颗粒体积对流场的影响。

1.1.2 气固混合相爆轰控制方程

基于圆盘形燃烧室的结构特征与前文假设， 在柱坐标系中建立如下控制方程组：

Ut+Fr+Grθ+Hz=R－Fr （1）

式中： U为守恒项； F， G， H为对流项； R为源项。

U=ρgYO2ρgYH2ρgYCO2ρgYH2OρgYN2ρgugρgvgρgwgρgEgρsρsusρsvsρswsρsEs ；

F=ρgugYO2ρgugYH2ρgugYCO2ρgugYH2OρgugYN2ρgu2g+pρgugvgρgugwg（ρgEg+p）ugρsusρsu2sρsusvsρsuswsusρsEs；

G=ρgvgYO2ρgvgYH2ρgvgYCO2ρgvgYH2OρgvgYN2ρgugvgρgv2g+pρgvgwg（ρgEg+p）vgρsvsρsusvsρsv2sρsvswsvsρsEs； H=ρgwgYO2ρgwgYH2ρgwgYCO2ρgwgYH2OρgwgYN2ρgugwgρgvgwgρgw2g+p（ρgEg+p）wgρswsρsuswsρsvswsρsw2swsρsEs；

R=－8m·s/3－16ω·－2ω·11m·s/318ω·0m·sus－Frm·svs－Fθm·sws－Fz－Qconv－（Frus+Fθvs+Fzws）+m·s（Es+p/ρs）－m·s－m·sus+Fr－m·svs+Fθ－m·sws+FzQconv+（Frus+Fθvs+Fzws）－m·s（Es+p/ρs）。 其中： 下标g， s分别表示气相和固相； 变量ρ， u， v， w， E， p分别表示密度、 径向速度、 周向速度、 轴向速度、 总能和气相压力； YO2， YH2， YCO2， YH2O， YN2分别表示气相中氧气、 氢气、 二氧化碳、 水蒸汽和氮气的质量分数， 且YO2+YH2+YCO2+YH2O+YN2=1； m·s为固相燃料质量消耗速率， ω·为氢气反应摩尔浓度消耗速率； Fr， Fθ， Fz分别为固体颗粒群在径向、 周向和轴向上受到的气相作用力； Qconv表示气相与固体颗粒群之间的对流换热。

能量方程中， 气相和固相总能分别满足如下关系：

Eg=h－pρ+12（u2g+v2g+w2g） （2）

Es=cvTs+12（u2s+v2s+w2s） （3）

式中： h为气相的比焓； cv为固相的比热容； Ts为固相温度。

本文所用固相燃料为煤粉， 鉴于前期所用煤粉固定碳含量为95%以上， 故在计算时主要考虑碳颗粒群的反应。 目前的研究中， 碳颗粒的燃烧主要受到表面动力学或物质扩散的限制。 扩散限制的反应速率依赖于氧化剂的浓度， 而动力学限制的反应速率则取决于颗粒表面处的反应速率。 因此， 鉴于Salvadori等人［19］的分析与验证， 碳颗粒群的总反应速率被定义为

m·s=4Ncπr2cpO21/ks+1/kd（4）

式中： Nc为颗粒数密度； rc为颗粒半径； pO2为氧气分压； ks和kd分别代表煤粉颗粒的动力学反应速率和扩散反应速率（单位为kg·m-2·s-1·Pa-1）。

ks=0.86exp－1.495×108RgTs（5）

kd=9.72DrcRg（Tg+Ts） （6）

式中： Tg为气相温度； Rg为气体常数； D为扩散系数。

D=4.24×10－7T 3/2g（1/MO2+1/MCO2）1/2p（V1/3O2+V1/3CO2）2（7）

式中： MO2， MCO2分别为氧气和二氧化碳的摩尔质量； VO2， VCO2分别为氧气和二氧化碳的摩爾体积。

本文所用气相燃料为氢气， 其主要作用为流化并辅助微米级煤粉爆轰， 因此， 这里采用一步反应［20］来描述其化学反应过程， 相应的反应速率由Arrhenius公式计算：

ω·=A［H2］m［O2］nexp（－EaRgTg）（8）

式中： A为化学反应指前因子， 取值为9.87×108； ［H2］， ［O2］分别为氢气、 氧气的摩尔浓度； m， n为反应级数， 这里均取值为1； Ea为活化能， 取值为5×107J/kmol。

气固两相流场中， 气体与固体颗粒间存在拖曳力［24］， 结合粉末颗粒的浓度， 可以在宏观上得到气相在各方向上对固相的作用力为

Fr=12Ncπr2cCDρg|Vg－Vs|（ug－us）

Fθ=12Ncπr2cCDρg|Vg－Vs|vg－vs）

Fz=12Ncπr2cCDρg|Vg－Vs|（wg－ws）  （9）

式中： Vg， Vs分别为气相和固相的速度矢量； CD为阻力系数。

CD与气固两相之间的相对雷诺数Re存在如下关系［24］：

CD=24Re1+16Re2/3（10）

Re=2ρgrc|Vg－Vs|μg（11）

式中： μg为气体粘度。

两相流动过程中， 由于温差和速度差等因素的影响， 气相与固相之间通过对流换热传递的热量为［24］

Qconv=2πrcNcλgNu（Tg－Ts） （12）

式中： λg为气体的导热系数； Nu为努塞尔数， 且

Nu=2+0.459Re0.55Pr0.33（13）

此处， 普朗特数Pr=0.72。

1.2 计算方法

爆轰过程的一个主要特征是有激波等强间断的传播。 为此， 爆轰过程的数值模拟对强间断的捕捉能力具有较高的要求。 故本文采用适用于强间断问题求解的时空守恒元与求解元方法（CE/SE方法）来对上述控制方程组进行求解。

CE/SE方法于1995年被Chang正式提出以来已发展30年左右， 它是基于空间通量与时间通量的守恒性原理推导而来， 目前被大量应用于爆轰等强间断问题的数值计算［1， 23-25］。 该方法作为一种全新的守恒型方程计算格式， 与传统数值方法（如有限差分法和有限体积法）不同， 它将时间视为第四个维度， 与空间中的三个维度进行统一划分处理。 在求解计算时， 该方法将整个时间-空间区域划分为若干个求解元。 在每个求解元内， 假设流场的变量是连续的， 并可用泰勒级数展开， 穿过相邻求解元的边界， 流场的变量可以是不连续的。 对于每个网格点对应的守恒元， 空间-时间密度矢量的积分通量是守恒的。

CE/SE方法计算精度高， 且无需采用算子分裂即可保证局部和整体的通量守恒， 无需Riemann分解即可实现激波捕捉， 在保证守恒性的同时， 能够准确地获取流场的激波和其他接触间断现象。

1.2.1 计算格式

在采用CE/SE方法对控制方程组进行求解时， 首先将整个计算域划分为交替的六面体网格， 网格点（i， j， k）位于六面体的中心位置， 且i， j， k分别是柱坐标系中的r， θ， z方向的网格点数。 通常， 在一个完整的时间步Δt内， 间隔分布的网格点分两步计算完成， 即前半个时间步求解一半网格点， 后半个时间步计算另一半网格点。 三维CE/SE方法守恒元与求解元的定义涉及到四维结构的构建， 其计算格式的推导过程详见文献［27］， 这里直接给出其在柱坐标系下的推导结果：

（U）li， j， k－Δt2（R－F/r）li， j， k=16 U+Δr4Ur－Δt4（Gθ/r+Hz）+  3ΔtΔrF+Δr6Fr+Δt4Ftl－1/2i－1/2，  j，  k+

U－Δr4Ur－Δt4（Gθ/r+Hz）－3ΔtΔrF－Δr6Fr+Δt4Ftl－1/2i+1/2，  j，  k+U+Δθ4Uθ－Δt4（Fr+Hz）+

3ΔtΔθG+Δθ6Gθ/r+Δt4Gtl－1/2i， j－1/2， k+U－Δθ4Uθ－Δt4（Fr+Hz）－ 3ΔtΔθG－Δθ6Gθ/r+Δt4Gtl－1/2i， j+1/2， k+U+Δz4Uz－Δt4（Fr+Gθ/r）+3ΔtΔzH+Δz6Hz+Δt4Htl－1/2i，  j，  k－1/2+

U－Δz4Uz－Δt4（Fr+Gθ/r）－

3ΔtΔzH－Δz6Hz+Δt4Htl－1/2i，  j，  k+1/2 （14）

式中： l为时间点。 Fr， Gθ， Hz， Ft， Gt， Ht与Ur， Uθ， Uz之间存在如下关系：

Fr=Fr=FU Ur=AUr

Gθ=Gθ=GU Uθ=BUθ

Hz=Hz=HU Uz=CUz （15）

Ft=A（R－F/r－AUr－BUθ/r－CUz）

Gt=B（R－F/r－AUr－BUθ/r－CUz）

Ht=C（R－F/r－AUr－BUθ/r－CUz）  （16）

其中， A=FU， B=GU和C=HU均为雅克比系数矩阵。 因此， 在使用该计算格式时只需再提供Ur， Uθ和Uz的算法即可。 因旋转爆轰波等强间断的存在， 本文在计算Ur， Uθ和Uz时宜采用加权函数法求解。 为此， 定义

（U±r）li，  j， k=±［（U′）li±1/2， j， k－（U）li， j， k］Δr/2（U±θ）li， j， k=±［（U′）li， j±1/2， k－（U）li， j， k］Δθ/2（U±z）li， j， k=±［（U′）li， j， k±1/2－（U）li， j， k］Δz/2  （17）

式中： （U′）li±1/2， j， k=U+Δt2Utl－1/2i±1/2， j， k；

（U′）li， j±1/2， k=U+Δt2Utl－1/2i， j±1/2， k；

（U′）li， j， k±1/2=U+Δt2Utl－1/2i， j， k±1/2。

进一步定义U（κ）r， U（κ）θ和U（κ）z（κ = 1， 2， 3， …， 8）分别对应于（U±r）li， j， k， （U±θ）li， j， k和（U±z）li， j， k的8种组合， 即

U（1）r=（U+r）li， j， k， U（1）θ=（U+θ）li， j， k， U（1）z=（U+z）li， j， k；

U（2）r=（U+r）li， j， k， U（2）θ=（U+θ）li， j， k， U（2）z=（U－z）li， j， k；

U（3）r=（U+r）li， j， k， U（3）θ=（U－θ）li， j， k， U（3）z=（U+z）li， j， k；

U（4）r=（U+r）li， j， k， U（4）θ=（U－θ）li， j， k， U（4）z=（U－z）li， j， k；

U（5）r=（U－r）li， j， k， U（5）θ=（U+θ）li， j， k， U（5）z=（U+z）li， j， k；

U（6）r=（U－r）li， j， k， U（6）θ=（U+θ）li， j， k， U（6）z=（U－z）li， j， k；

U（7）r=（U－r）li， j， k， U（7）θ=（U－θ）li， j， k， U（7）z=（U+z）li， j， k；

U（8）r=（U－r）li， j， k， U（8）θ=（U－θ）li， j， k， U（8）z=（U－z）li， j， k。

基于此， Ur， Uθ和Uz可通过如下加权函数求解：

（Ur）li， j， k=0， 若ψ（ζ）=0 （ζ=1， 2， 3，…， 8）∑8κ=1［（W（κ））αU（κ）r］/∑8κ=1（W（κ））α， 其他

（Uθ）li， j， k=0， 若ψ（ζ）=0 （ζ=1， 2， 3， …， 8）

∑8κ=1［（W（κ））αU（κ）θ］/∑8κ=1（W（κ））α， 其他

（Uz）li， j， k=0， 若ψ（ζ）=0 （ζ=1， 2， 3， …， 8）

∑8κ=1［（W（κ））αU（κ）z］/∑8κ=1（W（κ））α， 其他  （18）

式中： α是一個可调常数， 通常取为1或者2；  W（κ）=∏8ζ=1， ζ≠κψ（ζ） ， ψ（ζ）=（U（ζ）r）2+（U（ζ）θ）2+（U（ζ）z）2。

1.2.2 时间步长及源项的处理

数值计算过程中， CE/SE方法时间步长Δt的选取通常应满足稳定性条件：

Δt=minεΔr|u|+c， εrΔθ|v|+c， εΔz|w|+c （19）

式中： c为声速； ε为小于1的数， 弱波的情况下可以取大一些， 而对于强波则可以取小一些。

由于爆轰过程中， 化学反应的特征时间远小于物理特征时间， 因此其源项是刚性的。 这里采用四阶龙格库塔法对刚性源项进行求解， 且求解的时间步长按如下关系选取：

Δtr－k=Δt/Nr－k （20）

式中： 迭代次数Nr－k通常取5 ～ 20。

本文刚性问题的处理步骤为： 在不考虑源项的情况下， 采用CE/SE方法先求解（U）li， j， k， 然后将所得结果作为初值， 再对常微分方程组dU/dt=R－F/r进行求解。

1.3 初始条件与边界条件

初始条件： 初始时刻燃烧室内填充有常温常压下总当量比为1.0的煤粉-氢气-空气可燃混合物， 其中氢气/空气当量比为0.7， 煤粉/空气当量比为0.3。 图2所示为点火区域分布图， 为快速获得稳定传播的旋转爆轰波， 在图2中角度为15°的红色区域设置高温高压高速的周向气流作为点火条件， 并在指定时间内将红色区域左侧的计算域缝合面设置为固壁。

边界条件： 如图1所示， 上壁面、 下壁面和喷注面为固壁边界； 喷注面上层入口阵列为煤粉-氢气混合燃料的压力入口边界， 下层入口阵列为空气的压力入口边界； 作为出口的内圆柱面为无反射自由边界［3］； 圆周缝合处在初始爆轰波传播270°后由初始的固壁边界转变为周期边界。 对于压力入口边界， 参数的计算与设置可参考文献［10］。 同时， 为便于燃料与空气的混合， 并增加燃料在燃烧室中的停留时间， 参照前期实验研究工作［6-8］， 燃料喷孔处煤粉-氢气混合燃料的喷注速度向空气喷孔一侧偏转， 而空气喷孔处气流的喷注速度则向圆周切线方向发生偏转。 对于出口无反射自由边界， 参数的具体计算与设置可参考文献［30］， 即出流为亚声速时， 采用外展出流边界条件， 而出流为超声速时， 采用自由出流边界条件。

2 计算结果与分析

网格大小的选择既要考虑计算精度， 同时也需要兼顾计算成本。 在进行数值模拟之前， 首先进行网格无关性验证。 本文分别在Mesh-1（1 040×124×44）， Mesh-2 （720×90×32）和Mesh-3（560×56×20） 三种网格规模下进行爆轰流场的数值计算， 得到各网格尺寸条件下同一时刻燃烧室外圆压力沿周向的分布情况如图3所示。 可知以上三种网格均能有效捕捉到爆轰波的强间断面， 且各网格规模下爆轰波压力差异均不超过2%。 因此， 本文选取720×90×32网格开展数值研究。

本文对以氢气-煤粉气固混合物为燃料、 空气为氧化剂的盘形旋转爆轰发动机进行三维流场数值计算。 前期的实验研究在文献［8］中已经证实， 在半径R2=75 mm， 厚度L=16 mm的圆盘形燃烧室中， 由260 g/s的空气、 5.3 g/s的氢气和6.7 g/s的微米级多孔球状无烟煤粉末组成的可燃混合物， 经高能点火可实现单个旋转爆轰波的稳定传播。 为进一步探究实验现象的本质， 同时验证数值计算的准确性， 本文数值计算均按上述条件进行设置， 根据实验中所用粉末材料的粒径检测结果， 计算中的球状颗粒粒径设为5 μm。

图4所示为旋转爆轰波在盘形燃烧室稳定传播阶段， 距离喷注面35 mm（即距燃烧室中轴线R=40 mm）处下壁面监测点的实验与计算压力时程曲线。 可以看出， 两曲线的波形及幅值大小基本一致， 证明文中所建计算模型可准确模拟非预混喷注条件下煤粉-氢气-空气在盘形燃烧室中的旋转爆轰过程。

2.1 气固混合相旋转爆轰结构特征

从高温高压点火开始， 经过数十圈的传播， 盘形燃烧室中最终形成一个沿逆时针传播的旋转爆轰波。 图5所示为下壁面上紧贴喷注面（距燃烧室中轴线R=75 mm）处监测点压力及爆轰波传播频率时程曲线， 可以看到该测点处所测压力峰值周期性振荡且平均值为5.54 MPa。 显然， 该值远大于图4监测点（R=40 mm）处压力峰值。 这主要是因为该点处于爆轰波波头可覆盖范围， 且位于两壁面交汇处， 壁面对激波的反射作用使其压力成倍增加， 而图4监测点距离喷注面35 mm， 距爆轰波波头相对较远， 所测压力峰值为斜激波压力， 压力峰值相对较低。 然而， 各测点处压力峰值大小的差异并不影响爆轰波传播频率的获取。 通过求取压力时程曲线上爆轰波或斜激波两相邻峰值时间差的倒数， 可近似得到爆轰波的瞬时传播频率。 图5显示， 该爆轰波的平均传播频率为4.211 kHz， 若以测点所在圆周（R=75 mm）计算， 则该旋转爆轰波的传播速度为1 984 m/s。 图6所示为8.25 ms时刻煤粉-氢气混合相燃料在空气中稳定爆轰时盘形燃烧室的压力及温度分布图， 可以看到爆轰波在喷注面的约束作用下沿着圆盘外缘稳定传播。

图6（a） 显示， 非预混喷注条件下， 旋转爆轰波波阵面因壁面约束和非预混喷注等原因呈现为不规则的曲面， 从喷注面方向看爆轰波阵面主要由上下两个弧面交汇而成， 交汇点位于盘形燃烧室下半层， 即靠近空气喷孔一侧。 喷注面上两弧面交汇处压力最大， 喷注面下半层整体压力明显高于上半层。 这一方面是因为激波交汇导致压力升高， 另一方面是交汇点位于喷注流量较大的空气一侧， 该区域密度较大， 相应地， 爆轰波经过时强度较大。 此外， 爆轰波沿圆周传播过程中， 爆轰波与上壁面之间存在一层很薄的已燃气体层， 由于该薄层无燃料及氧化剂的喷注， 因此爆轰波在该薄层的燃烧产物中产生一透射激波随爆轰波波阵面一起向前传播。 而沿径向方向， 不规则的爆轰波波阵面在爆轰波与燃烧室出口之间又透射產生一系列强度较弱的激波。

图6（b）显示， 爆轰波前的区域在低温燃料及空气的高速喷注作用下， 燃烧室外缘靠近喷注面板的中间区域形成低温燃料填充区域， 而在靠近上下壁面的喷注盲区， 因上一循环留下的高温产物驻留而保持相对较高的温度。 喷注面附近， 爆轰波扫过后， 在稀疏膨胀波的作用下， 燃烧室压力迅速下降， 喷孔在压差的作用下向燃烧室中喷注低温燃料及空气， 使得波后喷孔附近混合气体的温度相较于贴近上下壁面的喷注盲区明显下降。 从喷注面方向看， 旋转爆轰波附近的温度场呈现出上下两个“V”型高温区分布的特点。

2.2 气固混合相爆轰化学反应分布特征

为了进一步探究煤粉-氢气-空气旋转爆轰波的组织形式， 图7给出了煤粉-氢气-空气旋转爆轰过程中各燃料组分的反应放热速率透视图。 从图中可以清晰地看到， 混合燃料反应放热速率的空间分布特征与旋转爆轰波的结构高度一致， 表明此处混合燃料（煤粉-氢气）与空气反应放热支持了旋转爆轰波的持续传播。

与实验喷注情形一致， 混合均匀的气固混合相燃料经燃料喷孔阵列向空气的喷注区域喷注。 但因为气相燃料与固相粉末燃料在物理性质上的巨大差异（主要表现为固体颗粒相较于气体具有更大的密度与惯性）， 混合燃料经空气的“过滤”作用后在燃烧室中会呈现出分层分布的特点。 进而导致各燃料的反应放热区域也呈现出分层分布的特点。 因此， 从图7中各燃料的放热速率空间分布来看， 煤粉-氢气混合相燃料旋转爆轰波在燃料与空气分开喷注的情形下依靠双反应区的爆轰组织形式自持传播。 其中， 煤粉与空气中的氧气主要在燃烧室贴近下壁面的区域反应放热以支持旋转爆轰波的传播， 而氢气与氧气发生反应的区域则主要位于煤粉反应区上方。 图7中煤粉化学反应区分布与图6中爆轰波处压力、 温度空间分布的高度耦合， 证明了煤粉参与爆轰并支持混合相旋转爆轰的连续传播。

2.3 流场三维分布特征

煤粉-氢气混合燃料与空气在轴向上分开喷注的方式， 使得燃料与氧化剂的掺混在轴向上存在较大差异， 进而导致流场剖面在轴向上发生变化， 而轴向上燃料喷孔向空气喷孔一侧倾斜喷注的方式也使得最佳掺混区域向下壁面偏移。 因此， 根据上一节中燃烧室压力温度在爆轰波处的轴向分布特征， 在轴向上选取若干切面及特征点展开分析与讨论。

图8所示为盘形燃烧室轴向不同切面上的压力与温度分布图。 因喷注掺混区域整体向下壁面偏移， 径向喷注速度的不均匀分布导致旋转爆轰波在上、 下壁面处的激波前沿在径向偏转趋势上存在明显差异。 中心切面处因受上下波阵面交汇的影响， 此处波阵面在径向上发生明显转折， 且转折处因激波汇聚而压力最大。 与中心切面不同， 上下壁面处的压力最大区域紧贴喷注面， 这主要是圆柱喷注面的收敛压缩作用导致的结果。 因燃料空气掺混区域向下壁面偏移， 所以爆轰波经过时靠近下壁面处化学反应更为剧烈， 可以看到， 爆轰波经过时下壁面处压力梯度较大， 上壁面处压力梯度较小。

图8（b）中各切面的温度分布显示， 盘形燃烧室中的高温燃气呈旋涡状流出燃烧室。 这一方面是因为喷注空气的切向分速度引起的旋转流动； 另一方面是旋转爆轰波传播过程中， 沿周向高速运动的斜激波压缩气体， 增加了气体的切向速度分量。 从轴向中心切面的温度分布可以看到， 低温燃料与空气的喷注掺混过程在径向上形成一环形的新鲜燃料预混区， 该区域的稳定存在保证了盘形燃烧室中旋转爆轰波的稳定传播。 因空气喷孔所在切面距离下壁面较近， 在高速的低温空气对流和新鲜的煤粉颗粒群吸热作用下， 下壁面上靠近喷注面的区域与整个切面其他区域相比， 在爆轰波到来之前也存在一温度较低的环形带。 爆轰波扫过时， 经过一段时间预热的煤粉与附近氧气快速反应释放热量， 支持旋转爆轰波的传播。 相比之下， 燃烧室上半区因属于新鲜来流喷注的盲区， 且距离喷孔相对较远， 在上壁面处并未形成明显的低温环形带。

图9所示为径向坐标R=75 mm时， 下壁面、 中心切面和上壁面处监测点的压力与温度时程曲线。 在掺混区向燃烧室下半区偏移以及圆柱形喷注面对流场的收敛压缩作用下，  下壁面处监测点平均压力峰值最高，  达到5.54 MPa。 其次是距离掺混区较近的中心切面处测点， 平均压力峰值为2.42 MPa； 相比之下， 离掺混区最远的上壁面处， 监测点平均压力峰值最低， 为2.18 MPa。 因处于喷注死角， 上壁面贴近喷注面处监测点的温度值基本维持在1  500 K以上，  且在爆轰波传播至对应相位角时， 平均温度峰值可达到3 063 K。 与之相比， 距离喷注区域更近的下壁面处， 贴近喷注面的位置虽位于喷注盲区， 但强烈的对流作用使得监测点温度在每一个旋转周期内都有一段时间维持在1 000 K以下。 喷注区域的中心切面处， 因新鲜低温氧化剂与燃料的喷注， 监测点温度在一半以上的时间内维持在500 K以下， 仅在爆轰波扫过时温度短暂超过1 500 K， 平均温度峰值可达到2 464 K。

盘形燃烧室扁平的结构特征， 使其在轴向上对爆轰波的传播产生了明显的约束。 为了进一步探究上下壁面对爆轰波传播的影响以及爆轰波后压力场的分布特点， 图10给出了距圆盘形燃烧室中心轴线不同距离的圆周切面压力分布展开图。 从各切面的压力分布间断来看， 因非预混喷注条件下爆轰波波阵面为不规则的曲面结构， 使得其与上下壁面之间形成不同的倾斜角度， 进一步导致了爆轰波波后区域激波在上下壁面处的交替反射。 圆周切面的半径越小， 距离出口越近， 距喷注面越远， 压力场受到爆轰波以及喷注面收敛压缩作用的影响越小。 因此， 紧贴喷注面的切面处因反射激波而产生的压力变化最为明显。 而随着圆周切面半径的减小， 同时考虑到侧向膨胀作用的影响， 因爆轰波传播引起的斜激波及波后一系列反射激波逐渐减弱， 切面上压力变化逐渐减缓。

从图8中压力和温度沿径向的分布来看， 新鲜燃料及氧化剂的喷注混合区域高度有限， 导致爆轰波在径向上的高度有所限制。 随着与喷注面距离的增加， 燃烧室中的强压缩扰动由爆轰波转变为斜激波， 其中， 斜激波以“弓形激波系”的形式呈现。 图8中各轴向圆盘切面的压力分布显示， 因爆轰波阵面在不同切面上形状有所差异， 爆轰波与燃烧室出口之间产生了一系列强度较弱的弓形激波。 这些弓形激波均汇聚于不规则的爆轰波阵面， 而发散于燃烧室出口。 因燃烧室为圆盘形， 激波沿圆周切向传播相同的距离， 距燃烧室中轴线越近的位置相位差越大， 因此由爆轰引起的“弓形激波系”前沿在相位上始终领先于爆轰波前沿。 这些规律在图例取值范围较小的图10中得到进一步体现。 从图10中可以看到， 随着与喷注面距离的增加， 爆轰波波阵面受非预混喷注的影响变小， 波阵面逐渐趋于平整。 并且从R=65 mm切面处开始， 各弓形激波之间的相位差随径向坐标R的减小快速增加。

图11所示为下壁面上“弓形激波系”的主要分布区域内， 同一相位角度不同径向位置的压力时程曲线。 从图中各监测点捕捉的第一道弓形激波压力扰动的上升沿开始时间点， 可以看出“弓形激波系”率先覆盖距燃烧室中轴线较近的位置。 而对以上各时间点作差， 从R=40 mm处开始， 到R=65 mm， R每增加5 mm， 测点捕捉到压力上升沿的时间间隔Δt依次为4.2 μs， 4.2 μs， 4.1 μs， 4.1 μs和4.2 μs。 因此， 在角速度基本不变的情况下， “弓形激波系”前沿的相位角随径向坐标R的增大线性减小。 结合2.1节中旋转爆轰波的传播频率， 可以计算得到径向坐标每间隔5 mm， “弓形激波系”前沿相位差Δθ=2πf·Δt= 0.108 ～ 0.111 rad。 因弓形激波强度较弱， 前文提到的因上下壁面反射而产生的反射激波在该区域内也相对较弱， 图11中各反射激波峰值压力的大小对比再次证明， 随着径向坐标R的减小， 反射激波逐渐减弱， 在R小于55 mm时， 弓形激波后无明显的反射激波存在。

3 结  论

为给固体燃料旋转爆轰发动机的研制提供技术支持与验证， 在国内外关于煤粉-氢气-空气旋转爆轰实验研究持续开展的背景下， 本文基于双流体模型， 在圆柱坐标系下建立了三维气固混合相燃料爆轰理论模型， 采用CE/SE方法计算求解了盘形燃烧室中煤粉-氢气-空气混合相旋转爆轰流场； 证实了微米级煤粉可支持旋转爆轰波的传播， 同時分析总结了粉末旋转爆轰的三维分布特征。 计算结果表明：

（1） 煤粉-氢气-空气旋转爆轰波波阵面因壁面约束和非预混喷注等原因呈现为不规则的曲面， 从喷注面方向看， 波阵面主要由上下两个弧面交汇而成， 喷注面上两弧面交汇处压力最大。 因喷注方式采用非预混小孔喷注， 喷注面与上下壁面之间形成喷注盲区， 旋转爆轰波附近的温度场呈现出上下两个“V”型高温区分布的特点。

（2） 因煤粉颗粒与氢气穿透空气区域的能力存在较大差异， 煤粉氢气混合流在穿越空气喷注区域时分层分布， 煤粉-氢气混合相燃料旋转爆轰波在非预混喷注的情形下依靠双反应区的爆轰组织形式自持传播。 其中， 煤粉与空气中的氧气主要在燃烧室贴近下壁面的区域反应放热以支持旋转爆轰波的传播， 而氢气与氧气发生反应的区域则主要位于煤粉反应区上方。

（3） 因非预混喷注条件下爆轰波波阵面为不规则的曲面结构， 使得其与上下壁面之间形成不同的倾斜角度， 进一步导致了爆轰波波后区域出现多重反射激波。

（4） 盘形燃烧室中爆轰波与燃烧室出口之间产生了一系列强度较弱的弓形激波， 这些弓形激波均汇聚于不规则的爆轰波阵面， 发散于燃烧室出口， 形成“弓形激波系”； 而该“弓形激波系”前沿的相位角坐标随燃烧室中径向坐标的增加线性减小。

参考文献：

［1］ 李宝星， 翁春生. 气体与液体两相连续旋转爆轰发动机爆轰波传播特性三维数值模拟研究［J］. 兵工学报， 2017， 38（7）： 1358-1367.

Li Baoxing， Weng Chunsheng. Three-Dimensional Numerical Simulation on the Propagation Characteristics of Detonation Wave in Gas-Liquid Two-Phase Continuous Rotating Detonation Engine［J］. Acta Armamentarii， 2017， 38（7）： 1358-1367.（in Chinese）

［2］ Meng H L， Xiao Q， Feng W K， et al. Air-Breathing Rotating Detonation Fueled by Liquid Kerosene in Cavity-Based Annular Combustor［J］. Aerospace Science and Technology， 2022， 122： 107407.

［3］ Tang X M， Wang J P， Shao Y T. Three-Dimensional Numerical Investigations of the Rotating Detonation Engine with a Hollow Combustor［J］. Combustion and Flame， 2015， 162（4）： 997-1008.

［4］ Lin W， Zhou J， Liu S J， et al. An Experimental Study on CH4/O2 Continuously Rotating Detonation Wave in a Hollow Combustion Chamber［J］. Experimental Thermal and Fluid Science， 2015， 62： 122-130.

［5］ Liu S J， Huang S Y， Peng H Y， et al. Characteristics of Methane-Air Continuous Rotating Detonation Wave in Hollow Chambers with Different Diameters［J］. Acta Astronautica， 2021， 183： 1-10.

［6］ Xu H， Ni X D， Su X J， et al. Experimental Investigation on the Application of the Coal Powder as Fuel in a Rotating Detonation Combustor［J］. Applied Thermal Engineering， 2022， 213： 118642.

［7］ 续晗， 罗永晨， 倪晓冬， 等. 铝粉燃料连续旋转爆轰发动机工作特性［J］. 兵工学报， 2022， 43（5）： 1046-1053.

Xu Han， Luo Yongchen， Ni Xiaodong， et al. Operating Characte-ristics of Aluminum Powder Rotating Detonation Engine［J］. Acta Armamentarii， 2022， 43（5）： 1046-1053.（in Chinese）

［8］ Ni X D， Xu H， Su X J， et al. Effects of Different Physical Properties of Anthracite Powder Fuel on Detonation Characteristics of a Rotating Detonation Engine［J］. Physics of Fluids， 2023， 35（5）： 053325.

［9］ 夏鎮娟， 张义宁， 马虎， 等. 点火位置对圆盘结构下旋转爆震波起爆过程的影响［J］. 推进技术， 2021， 42（4）： 805-814.

Xia Zhenjuan， Zhang Yining， Ma Hu， et al. Effects of Ignition Position on Initiation Process of Rotating Detonation Wave in Plane-Radial Structure［J］. Journal of Propulsion Technology， 2021， 42（4）： 805-814.（in Chinese）

［10］ 夏镇娟， 武晓松， 马虎， 等. 圆盘结构下旋转爆震波的二维数值研究［J］. 推进技术， 2017， 38（6）： 1409-1418.

Xia Zhenjuan， Wu Xiaosong， Ma Hu， et al. Two-Dimensional Numerical Simulation of Rotating Detonation Wave in Plane-Radial Structure［J］. Journal of Propulsion Technology， 2017， 38（6）： 1409-1418.（in Chinese）

［11］ Ma J Z， Luan M Y， Xia Z J， et al. Recent Progress， Development Trends， and Consideration of Continuous Detonation Engines［J］. AIAA Journal， 2020， 58（12）： 4976-5035.

［12］ Bykovskii F A， Zhdan S A， Vedernikov E F， et al. Continuous and Pulsed Detonation of a Coal-Air Mixture［J］. Doklady Phy-sics， 2010， 55（3）： 142-144.

［13］ Bykovskii F A， Zhdan S A， Vedernikov E F， et al. Detonation Combustion of Coal［J］. Combustion， Explosion， and Shock Waves， 2012， 48（2）： 203-208.

［14］ Bykovskii F A， Zhdan S A， Vedernikov E F， et al. Continuous Spin Detonation of a Coal-Air Mixture in a Flow-Type Plane-Radial Combustor［J］. Combustion， Explosion， and Shock Waves， 2013， 49（6）： 705-711.

［15］ Bykovskii F A， Zhdan S A， Vedernikov E F， et al. Detonation Burning of Anthracite and Lignite Particles in a Flow-Type Radial Combustor［J］. Combustion， Explosion， and Shock Waves， 2016， 52（6）： 703-712.

［16］ Bykovskii F A， Vedernikov E F， Zholobov Y A. Detonation Combustion of Lignite with Titanium Dioxide and Water Additives in Air［J］. Combustion， Explosion， and Shock Waves， 2017， 53（4）： 453-460.

［17］ Dunn I B， Malik V， Flores W， et al. Experimental and Theoretical Analysis of Carbon Driven Detonation Waves in a Heteroge-neously Premixed Rotating Detonation Engine［J］. Fuel， 2021， 302： 121128.

［18］ Dunn I B， Flores W， Morales A， et al. Carbon-Based Multi-Phase Rotating Detonation Engine［J］. Journal of Energy Resources Technology， 2022， 144（4）： 042101.

［19］ Salvadori M， Dunn I B， Sosa J， et al. Numerical Investigation of Shock-Induced Combustion of Coal-H2-Air Mixtures in a Unwrapped Non-Premixed Detonation Channel［C］∥AIAA Scitech 2020 Forum， 2020.

［20］ Zhu W C， Wang Y H， Wang J P. Flow Field of a Rotating Detonation Engine Fueled by Carbon［J］. Physics of Fluids， 2022， 34（7）： 073311.

［21］ 李世全， 楊帆， 王宇辉， 等. 当量比对铝粉/空气旋转爆轰发动机流场影响的数值模拟［J/OL］. 航空动力学报， doi： 10.13224/j.cnki.jasp. 20210560.

Li Shiquan， Yang Fan， Wang Yuhui， et al. Numerical Simulations on Effect of Equivalence Ratio on Flow Field in Aluminum/Air Rotating Detonation Engines［J/OL］. Journal of Aerospace Power， doi： org/10.13224/j.cnki.jasp.20210560. （in Chinese）

［22］ Chang S C. The Method of Space-Time Conservation Element and Solution Element-A New Approach for Solving the Naver-Stokes and Euler Equations［J］. Journal of computational physics， 1995， 119（2）： 295-324.

［23］ Wang G， Zhang D L， Liu K X， et al. An Improved CE/SE Scheme for Numerical Simulation of Gaseous and Two-Phase Detonations［J］. Computers and Fluids， 2010， 39（1）： 168-177.

［24］ Zhang Z J， Wen C Y， Liu Y F， et al. Application of CE/SE Method to Gas-Particle Two-Phase Detonations under an Eulerian-Lagrangian Framework［J］. Journal of Computational Physics， 2019， 394： 18-40.

［25］ 馬丹花， 翁春生. 二维守恒元和求解元方法在两相爆轰流场计算中的应用［J］. 燃烧科学与技术， 2010， 16（1）： 85-91.

Ma Danhua， Weng Chunsheng. Application of Two-Dimensional CE/SE Method Calculation of Two-Phase Detonation Flow Field［J］. Journal of Combustion Science and Technology， 2010， 16（1）： 85-91.（in Chinese）

［26］ 林玲， 翁春生. 等离子体射流点火对燃烧转爆轰影响的二维数值计算［J］. 兵工学报， 2014， 35（9）： 1428-1435.

Lin Ling， Weng Chunsheng. Two-Dimensional Numerical Calculation for the Influence of Plasma Jet Ignition on Deflagration-to-Deto-nation Transition［J］. Acta Armamentarii， 2014， 35（9）： 1428-1435.（in Chinese）

［27］ 翁春生， 王浩. 计算内弹道学［M］. 北京： 国防工业出版社， 2006： 290-306.

Weng Chunsheng， Wang Hao. Computational Interior Ballistics［M］. Beijing： National Defense Industry Press， 2006： 290-306.（in Chinese）

［28］ Bruce P J K， Babinsky H. Unsteady Shock Wave Dynamics［J］. Journal of Fluid Mechanics， 2008， 603： 463-473.

［29］ Oran E S， Weber J W， Stefaniw E I， et al. A Numerical Study of a Two-Dimensional H2-O2-Ar Detonation Using a Detailed Chemical Reaction Model［J］. Combustion and Flame， 1998， 113（1/2）： 147–163.

［30］ 王健平， 姚松柏. 连续爆轰发动机原理与技术［M］. 北京： 科学出版社， 2018： 68-72.

Wang Jianping， Yao Songbai. Principle And Technology Of Continuous Detonation Engine［M］. Beijing： Science Press， 2018： 68-72.（in Chinese）

Numerical Study on the Two-Phase Detonation Flow Field of

Coal Powder-Hydrogen in a Rotating Detonation Engine

Weng Chunsheng， Ni Xiaodong*， Xu Han

（National Key Laboratory of Transient Physics， Nanjing University of Science and Technology， Nanjing 210094， China）

Abstract： Rotating detonation engines （RDEs） have attracted significant attention in the aerospace field due to their high thermal cycle efficiency. This study explores the characteristics of rotating detonation in gas-solid mixed fuel within an air atmosphere. A theoretical detonation model for gas-solid mixed-phase fuel is developed in a cylindrical coordinate system. Using the three-dimensional conservation element and solution element （CE/SE） method， numerical simulations of the detonation process in a disc combustor are conducted in three dimensions. The flow field structure during the stable propagation of gas-solid mixed-phase rotating detonation waves is calculated. Subsequently， the analysis is conducted on the distribution characteristics of thermodynamic parameters in the detonation flow field， the distribution of chemical reaction zone， and the characteristics of wave system following the rotating detonation wave. The research results demonstrate that micron-sized coal powder can undergo a rapid reaction and support the stable propagation of rotating detonation waves with the assistance of hydrogen. However， due to the flat structural characteristics of the disc combustor and the non-premixed injection mode， the wave front of detonation wave appears as an irregular surface. The difference in the ability of coal powder and hydrogen jet to penetrate the air layer leads to the separation of the chemical reaction zone of the gas-solid mixed fuel， ultimately supporting the stable propagation of rotating detonation wave in the form of a dual reaction zone. The irregular surface structure of the rotating detonation wave leads to multiple reflections in the flat combustor， consequently causing the regular appearance of multiple reflected shock waves after the detonation wave.

Key words： gas-solid mixed-phase detonation； detonation wave； powder fuel； coal powder； non-premixed； numerical simulation