低频爆轰不稳定性形成机理的数值模拟研究*

2021-11-19张允祯荣光耀王健平

爆炸与冲击 2021年9期

张允祯，程杪，荣光耀，王健平

（北京大学工学院力学与工程科学系，北京 100871）

爆轰是一种激波与波后化学反应区强烈耦合的物理现象。爆轰现象自被发现以来，已受到了广泛的关注，并已形成了相关的研究理论。目前而言，被普遍认可并被广泛应用的爆轰理论主要有两种，分别为C-J理论和ZND模型[1-2]。由于爆轰相较于燃烧具有熵增小、热效率高等特点[3]，因此研究人员也一直致力于将爆轰引入到化学推进系统中，以期实现发动机性能的大幅提升。而近年来，连续爆轰发动机（rotating detonation engine）更是成为国际上研究的热点。

图1为连续爆轰发动机的原理示意图。在发动机运行过程中，在头部进气壁面附近有一道或多道沿圆周方向连续旋转传播的爆轰波头（detonation front），推进剂从头部进气壁面连续不断地喷入燃烧室内，在爆轰波头前形成一个新鲜气体层（fresh gas layer）。在爆轰波头的斜后方附着有一道斜激波（oblique shock wave）和一道接触间断（contact surface），接触间断两侧分别为上周和本周的爆轰产物（burnt gas），高温高压的爆轰产物从尾部高速喷出，从而获得所需的推力。

图1 连续爆轰发动机原理（二维圆柱流场及其展开结构）Fig.1 Schemetic of RDE(2D cylindrical flow field and its unfolding structure)

20世纪50年代末至60年代，Voitsekhovskii[4]和Nicholls等[5]分别对连续爆轰发动机进行了验证性研究。此后由于研究条件的限制，连续爆轰发动机的相关研究陷入近乎停滞的状态，直到Bykovskii等[6]重新在实验上实现了较长时间的连续旋转爆轰。此后世界各地学者们[7-17]对相关领域开展了大量的理论和实验研究。

对于一种推进系统，其能够稳定可靠地运行是研究人员追求的最终目标，因此连续爆轰发动机的稳定性研究是一个不可避免的课题。在以往针对连续爆轰发动机的研究中，人们发现，在连续爆轰发动机运行过程中存在各式各样的不稳定性，并针对这些不稳定性开展了一定的研究。在这个过程中，一些研究人员对这些不稳定性现象进行了分类。Wang等[18]根据实验测得的不稳定信号的振荡频率将这些不稳定性分为低频、中频、高频不稳定性，并指出低频爆轰不稳定性是由燃烧室和进气室的声学耦合引起的。Anand等[19-21]按照压力信号的特征将所发现的不稳定现象分为4类，并针对这些不稳定性开展了持续的、系统的实验研究。针对低频爆轰不稳定性（low frequency instability,LFI），Anand等[21]通过实验开展了深入研究，他们认为产生这种低频爆轰不稳定性的一个原因是由于燃烧室头部的激波与爆轰波相互作用所致。此外，Zhang等[22]的研究表明，燃烧室头部的激波会引起新鲜气体层的不规则分布，进而导致了低频爆轰不稳定性，但是研究未能揭示不规则新鲜气体层与压力峰值低频振荡之间的直接关联。此外，实验和数值上均未给出激波的来源、发展及导致低频爆轰不稳定性的具体作用过程。本文将通过数值模拟的方式针对低频爆轰不稳定性现象开展进一步的研究。

1 数值方法和物理模型

1.1 数值方法

本文采用H2/Air 9组分19步基元化学反应模型针对同轴圆环腔燃烧室开展二维的数值模拟研究。爆轰波的传播速度非常快，一般都在千米每秒的量级，因此在对爆轰波的数值模拟中，一般忽略组分的扩散、黏性以及热传导[23]的影响，采用二维圆柱坐标系下的带源项的守恒型的Euler方程作为控制方程[22,24]：

其中

式中：U为守恒量，F和G为通量向量，S为化学反应源项；ρ为混合气体总密度，ρi为各组分气体的密度，需满足归一化条件；E为混合气体单位体积的总能，h为混合气体总比焓；uθ为周向速度分量，uz为轴向速度分量；为第i个组分的化学生成速率，由基元化学反应模型得到；hi为第i个组分的比焓，与混合气体的温度有关。由于温度、压力、能量强烈耦合，因此联立式(3)、(5)，采用Newton-Rapson迭代法求解温度。热力学参数采用温度的多项式进行拟合。

在空间上采用Steger-Warming矢通量分裂法对通量向量进行分裂，采用五阶WENO格式对通量进行重构，时间上采用2阶TVD Runge-Kutta方法推进求解。在对化学反应的处理上，采用解耦算法将流动和化学反应进行解耦，并采用刘君等[25]提出的半隐式格式求解化学反应源项以解决求解源项时产生的刚性问题，并同时提高计算效率。

1.2 物理模型

本文将对同轴圆环腔燃烧室内存在的低频爆轰不稳定性开展数值模拟研究。同轴圆环腔燃烧室的环腔内外径之差远小于半径，因此在计算中可以将其简化为一个没有厚度的圆柱面，并忽略径向的影响，于是可将此圆柱面沿母线展开，在二维平面上进行计算。图2为一轴向长度（L）为0.04 m、半径（R）为0.01 m的燃烧室的平面展开计算域。其中下边界为进气壁面，上边界为燃烧室尾部出口。按化学恰当比预混好的H2/Air混合气从进气壁面喷入燃烧室，入流总压ptot=5×105Pa，入流总温Ttot=360 K。

图2 初始时刻计算域Fig.2 Computational domain at the initial time

初始时刻，燃烧室内充满按照化学恰当比预混好的新鲜可燃气体，初始压强p0和初始温度T0分别为1×105Pa和300 K。在头部壁面附近覆盖一个典型的一维爆轰波，用以起爆周向传播的旋转爆轰波。

边界条件采用如下设置。

（1）下边界为入流边界，采用收敛喷管入流边界。将每个网格点视为一个收敛喷管，在每个网格点上对物理量按照收敛喷管的条件进行设置。进气壁面上的压强pw为收敛喷管的背压，入流总压为ptot，入流总温为Ttot。收敛喷管的计算中有两个关键参数，入流总压ptot和临界压强pcr：(a)若pw＞ptot，新鲜气体无法喷入燃烧室，此时采用固壁反射边界条件，速度矢量在壁面两侧沿法向相反，沿切向相等，密度、压力、温度等标量值在壁面两侧对称相等；(b)若pw＞ptot＞pcr，新鲜气体以亚声速喷入燃烧室；(c)若pw＜pcr，喷管壅塞，新鲜气体以声速喷入燃烧室。

（2）上边界为尾部出流边界，采用亚声速/超声速出流边界条件。(a)若出口处的马赫数Ma＜1，气体以亚音速出流，边界处的压力取为周围环境的压力，其余物理量由边界内部的网格点上对应的物理量插值得到；(b)若Ma＞1，气体以超声速喷出燃烧室，下游的流场不能影响上游的流场，因此虚拟网格的值可随意设定，一般设定与出口处的值相等。

（3）左右边界为圆柱壁面上同一条母线，设定周期性边界条件，即将一侧的边界内侧节点上的值赋值给另一侧的边界外侧的虚拟网格点上。

由于采用基元化学反应模型，计算量较大，因此采用MPI并行计算，以提高计算效率。图3为采用0.1 mm和0.2 mm网格计算所得流场的温度云√图和特征压力梯度对数云图。特征压力梯度对数定义为[26]：dp∗=0.8Exp[−100|∇p|/max|∇p|]，其中|∇p|=(∂p/∂x)2+(∂p/∂y)2。从温度云图可以看到，两种尺寸的网格均能够模拟出在燃烧室头部周向传播的旋转爆轰波，且二者得到的流场基本一致。但从特征压力梯度对数云图来看，0.2 mm的网格得到的流场的波系结构较粗糙，不如0.1 mm的网格得到的流场展示出来的细节丰富。由于本文的研究中需要针对流场中的细微结构进行分析，因此采用0.1 mm的网格进行计算。

图3 不同网格尺寸的计算结果（温度云图和特征压力梯度对数云图）Fig.3 Calculation results of different grid sizes(temperature nephogram and pressure gradient logarithmic nephogram)

2 结果与讨论

在上述的初始及边界条件下，最终在燃烧室内形成了一个沿周向传播的爆轰波，记录燃烧室进气壁面上某采样点的压强随时间的变化，得到p-t曲线，如图4(a)所示。图4(b)为Anand等[21]的H2/Air连续爆轰发动机实验结果，3种不同颜色为沿环向均布的3个位置处的压力时程曲线。数值模拟结果与实验结果高度类似，爆轰波峰值均呈现出周期性起伏振荡的现象。计算爆轰波在燃烧室内传播的频率fD（会受到燃烧室尺寸的影响）和压力振荡频率fosc，数值模拟的结果分别为fD≈30 kHz和fosc≈1.5 kHz，两者之比fosc/fD≈0.05；Anand的实验结果两者分别约为2 kHz和150 Hz，两者之比。可以看到，数值和实验结果均非常小且十分相近，在同一量级。因此，一般将这种现象称为低频爆轰不稳定性。

图4 采样点（传感器处）的压强-时间曲线Fig.4 The pressure-time track of sampling points(or the points where the sensors are placed)

2.1 爆轰波头上压强的周期性变化

为探究这种低频爆轰不稳定性产生的原因，针对燃烧室内爆轰波的传播过程进行重点研究。图5为340～402 µs进气壁面上最大压强值随时间的变化。可以发现，进气壁面上的最大压强也呈现出了周期性起伏振荡的现象。选取黑色虚线框内（374～385 µs）的一个周期，详细观察燃烧室内的流场。图6为这段时间内爆轰波头在燃烧室头部的传播过程，其中白色曲线为接触间断，即已燃气与未燃气的分界线（因此白色曲线与进气壁面之间即为新鲜气体层）。图7为这段时间内爆轰波头上的压强沿轴向的分布情况，横轴代表轴向的距离，纵轴代表相应距离上爆轰波头上的压强。在这种工况下，一般认为压强超过4.5 MPa为爆轰波。

图5 进气壁面上最大压强随时间的变化Fig.5 Track of the peak pressure on the inlet wall

图6 374～384 µs爆轰波头的传播情况Fig.6 Propagation of detonation front between 374−384 µs

从图6和图7可以看到，爆轰波头上的压强分布在时间上并不一致，会存在一些向进气壁面运动的高压区。在374 µs时，爆轰波头的高压区正与进气壁面接触（图6），因此此时进气壁面最大压强达到最大值（图5）。在374～378 µs之间，爆轰波越过新鲜气体层上的凹点（以下将之称为进气阻滞点(injet blocking point,IBP)，在2.3节中会解释原因），波前的新鲜气体层的宽度急剧变宽（图6）。此时，在爆轰波头与新鲜气/已燃气分界线相交的地方产生了一道局部爆炸。这个现象可以直观地从图7中377～379 µs的波头压强分布曲线上看出。随后，这道局部爆炸产生的环形激波使得爆轰波头上产生一个高压区，这个高压区向进气壁面运动，与进气壁面相遇（如图7所示）从而使得进气壁面上的最大压强值起伏变化。由于爆轰波周期性与进气阻滞点相遇，爆轰波头上压强的轴向分布周期性变化，从而进气壁面的最大压强便产生了周期性的振荡。

图7 374～385 µs爆轰波头上的压强轴向分布Fig.7 Axial distribution of the pressure on the detonation front between 374−385 µs

像这种爆轰波头上产生局部爆炸的现象，其他学者在研究中也有发现。Hishida等[27]在他们的数值研究中曾经发现过类似的现象并采用了较精细的网格进行了研究，认为爆轰波越过新鲜气体层上的凹点时，爆轰波在新鲜气/已燃气分界线处会产生未燃气包，从而会形成局部爆炸。最近，Athmanathan等[28]在开展的连续爆轰发动机可视化实验中观察到，爆轰波头上存在所谓的强爆轰区，且这些强爆轰区会向燃烧室的头部方向运动。他们认为这种现象是由于爆轰波头前新鲜气体层高度的变化导致的。

2.2 低频爆轰不稳定性产生的机理

图8为340～402 µs之间进气壁面上最大压强（爆轰波头与进气壁面相交处的压强）的变化，蓝色虚线为爆轰波头经过进气阻滞点的时刻。将爆轰波经过两个进气阻滞点之间的时间记为 ∆t，相邻的两个进气阻滞点间的距离为 ∆x，爆轰波速D视为恒定的。由以上研究可以看到，进气壁面上相邻两个进气阻滞点间的距离 ∆x基本一致且保持不变，因此爆轰波经过两个进气阻滞点间的时间 ∆t=∆x/D也基本恒定，这一点可以很直观地从图8中看出。可以看到，爆轰波每经过任意相邻的两个进气阻滞点的时间里，进气壁面上最大压强的历程基本一致，总体上呈现出较规则的正弦振荡。

图8 进气壁面上最大压强的变化（数值模拟结果）Fig.8 Track of the peak pressure on the inlet wall(numerical result)

如图9(a)所示，当爆轰波处于两进气阻滞点（红色方块处）之间的某个位置时，设爆轰波头距离上个进气阻滞点（设为第n个）的距离为l。此时，爆轰波距离经过第n个进气阻滞点时已经过了t′=l/D的时间，于是此时进气壁面上的最大压强即为如图9(b)所示的 (n∆t+t′)时刻黑色虚线所指示的压强。

图9 采样点压强振荡的形成机理Fig.9 The schematic diagram of mechanism of the pressure oscillation at the sampling point

若此时爆轰波头正好与进气壁面上的采样点（sampling point，SP，图9(a)中为黑色方块处）相遇，则此刻采样点的压强正好达到压强峰值，且该压强峰值的大小为图9(b)所示的n∆t+t′时刻黑色虚线所指示的压强。在进一步的研究中发现，这些进气阻滞点产生的位置会在进气壁面上沿圆周方向缓慢移动，如图10所示，其中w(H2)为H2的质量分数。因此对于进气壁面上某固定的采样点来说，爆轰波每次与之相遇时，爆轰波头相对于上一个进气阻滞点的位置l会缓慢变化。对于每个确定的l(0≤l≤∆x)，都能够在图9(b)上有一个对应的压强，此压强即为本周爆轰波旋转传播过程中采样点的压强峰值。随着l从 ∆x逐渐减小到零，采样点的压强峰值便完成了一个周期的振荡。又由于进气阻滞点所产生的位置的移动速度v非常小，因此采样点的压强-时间曲线上压强峰值的振荡频率f=v/∆x也就很低，从而产生了所谓的低频爆轰不稳定性。

图10 爆轰波每旋转两周进气阻滞点产生的位置Fig.10 The position of IBP in every two cycles of detonation wave rotation

2.3 不规则新鲜气体层的形成

从上述研究中可以看到，新鲜气体层的不规则分布对于低频爆轰不稳定性的产生起着关键作用。本节将阐释不规则新鲜气体层的形成原因。

图11为进气壁面上进气阻滞点产生时该点周围H2质量分数的分布云图。从图11中可以看出，391～393 µs之间，进气壁面上Y=0.03 m附近的进气被阻滞，该区域两侧的进气均先于该区域的进气，其中0.029 m处的进气被阻滞的时间最长。分别在Y=0.027,0.029,0.032 m处设置采样点，记录进气被阻滞处前后的压强变化，并将之与稳定的进气情况进行对比，如图12所示。进气壁面压强低于进气总压（红色虚线，500 kPa）被认为气体可以喷入燃烧室。

图11 进气阻滞点的产生Fig.11 Generation of injet blocking points

图12 进气被阻滞处及左右各一点在进气阻滞点产生时的压强变化和稳定进气情况Fig.12 Pressure track of the place where intake process is interrupted and one point on each side when IBP is being generated,and the stable intake situation

与稳定情况进行对比，可以发现，发生进气阻滞的情况下，开始进气的前后进气壁面上的压强在时间上分布并不平滑，有一系列的激波沿着与爆轰波相反的方向从0.032 m处向0.029 m处运动（图12中绿色圆圈所示）。这些激波传到0.032 m处时，该点的压强还远高于进气总压，因此只是在压强曲线位于进气总压上面的部分引起了一定的振荡，在一定程度上使得该区域的进气过程滞后了一些，但对进气的阻滞效果十分有限；当这些激波传播到0.029 m处时，该点的压强刚好下降到进气总压，新鲜气体即将喷入燃烧室，但是这些激波的到来引起了即将进气的地方压强产生强烈振荡，并使压强保持在进气总压以上很长一段时间，即将喷入的新鲜气体被阻止进入燃烧室，进气被强行打断，于是在391～393 µs形成了图11中0.029 m处的进气阻滞点；同时，当这些激波经过0.029 m并引起进气阻滞的时候，上游（0.032 m附近）的压强却已经下降到进气总压以下，新鲜气体开始喷入燃烧室；随后，当这些激波越过0.029 m（进气被长时间阻滞处）到达下游的0.027 m处时，如前所述，该点的压强也早已下降到了进气总压以下并开始进气，并且这些激波随着传播距离的增加其强度也逐渐衰减，因此在0.027 m处不能阻碍进气过程，只是在该点的压强曲线上位于进气总压以下的部分引起了一定的压强回升。因此，这些激波仅在与新鲜气体层顶点（即压强刚下降到进气总压以下即将开始进气的点）相遇的地方间断性地导致了进气阻滞点的产生。

接下来，将从流场内波系结构的发展出发，揭示这些激波的来源、发展与导致进气阻滞点产生的详细过程。首先，由Chen等[29]的研究可知，当反传激波与爆轰波相遇后，反传激波的强度将得到大幅增强。图13为362～393 µs燃烧室头部波系结构的发展过程。从图13(a)中可以看到，爆轰波后存在一些强度较强的激波（黑色箭头所示），这些激波不断地向后传播，其强度不断衰减，最终进入新鲜气体层（371～374 µs），在下周的爆轰波前形成了一个横波区。在381～384 µs之间（图13(b)），这些波前横波（黑色箭头所示）与爆轰波相遇后被大幅增强，在波后重新形成了一个更强的横波区（红色箭头所示）。之后，在385～393 µs之间（图13(c)），这些被增强后的横波（黑色箭头所示）向进气壁面运动，与进气壁面相撞，发生反射。反射点与新鲜气体层的顶点（白色箭头所示）相遇，使得该点处的压强长时间维持在进气总压以上，该点的进气过程被长时间打断，于是形成了新鲜气体层上的进气阻滞点。图13(c)的白色圆圈内可以看到一道较强的激波在进气壁面上的反射点与新鲜气体层顶点相遇的情景。