一维爆轰波在扰动来流中传播的数值研究

郗雪辰，杨鹏飞，滕宏辉

(1.北京理工大学宇航学院，北京 100081;2.山西警察学院治安系，山西太原 030401;3.中国科学院力学研究所高温气体动力学国家重点实验室，北京 100190;4.中国科学院大学工程科学学院，北京 100049)

引 言

与传统的爆燃燃烧相比，爆轰燃烧依靠强激波压缩进行点火并实现自持传播，具有能量转化速率快、波传播速度高的特点。基于爆轰燃烧的航空航天推进系统可望实现较高的热循环效率，且燃烧室尺寸小、结构质量小[1-2]，在高超声速推进领域具有潜在的工程应用价值。气相爆轰波涉及流动、燃烧及其耦合作用等复杂的物理-化学过程，爆轰波面结构复杂多变，波面失稳后形成的非定常波系对燃烧组织和发动机稳定工作构成挑战[3]。波面失稳源于化学反应和前导激波不断变化的耦合关系，宏观上表现为多道横波与入射波共同构成的复杂波面。弱耦合的入射激波与强耦合的Mach杆，通过横波联结，三者的交点被称为三波点。研究者利用烟迹技术记录了三波点的运动轨迹[4]，通常称为胞格结构，并用来表征爆轰波的动力学过程。过去几十年，研究者围绕爆轰的胞格结构开展了大量的实验研究，获取了氢气、甲烷、乙炔、乙烯等气体燃料在不同状态下的胞格特性[5]。随着计算机技术的发展，利用数值模拟围绕胞格结构开展了大量研究，深化了对气相爆轰的认识[6-9]。

虽然真实爆轰波均为三维，但是一维爆轰波在很多情况下也是一种很好的宏观近似，反映了激波和放热在传播方向上的耦合关系，可作为理解多维爆轰波不稳定的基础，因此得到了广泛的研究。Fickett等[10]针对一维爆轰波开展了开创性的研究，Erpenbeck[11]对一维和二维爆轰波进行了稳定性分析，利用渐进理论得到了爆轰波的稳定性边界。He等[12]研究了活化能对一维爆轰波的影响，发现随着活化能的增加，爆轰波从稳定传播发展为间歇爆轰波，最后由于振荡幅度过大，爆轰波不再能稳定传播。以上研究假设波前气体状态固定不变。事实上，在爆轰发动机研究的过程中发现，高速来流中爆轰波的波系结构对来流参数的变化非常敏感[13]，且来流状态的变化还可能引起波前组分不均匀，而这种不均匀性也被认为是造成发动机预测性能与实测性能之间差异的重要原因之一[14-16]。因此研究爆轰波在非定常来流中的传播就有了非常重要的实际意义。

非定常来流问题之前多围绕来流状态突变和瞬时离散扰动等进行[17-20]。近来，Kasimov等[21]采用Burgers模型研究了一维爆轰波对连续扰动的响应过程。结果显示，连续小扰动不会影响一维稳定爆轰波的动力学过程，但会激发脉冲爆轰波的内在不稳定性，使其呈现出更为复杂的动力学特征。Burgers模型获得的结果与单步反应模型的结果定性一致。Kim等[22]采用单步反应模型对一维脉冲爆轰波在连续密度扰动中的传播进行了数值模拟，研究了扰动波长对爆轰波振荡行为的影响。结果发现对于自身不稳定的爆轰波，波长在一定范围内的扰动，反而会抑制一维爆轰的不稳定性。这一结果不仅验证了非线性理论在爆轰动力学研究领域应用的可行性，也为研究利用微扰动技术控制动态爆轰参数等问题提供了新的思路，具有非常重要的启发意义。

此前的研究采用简化的单步反应模型，该模型在讨论爆轰波气体动力学和化学动力学之间的相互作用时，存在一定的局限性。两步反应模型可更好地模拟真实爆轰中的诱导-放热反应过程，且具有比单步反应更多的自由参数，更利于开展数值研究。本文采用两步诱导-放热反应模型，模拟一维爆轰波在连续正弦密度扰动中的传播，主要关注扰动波长对稳定爆轰波动力学特性的影响。模拟获得了一些新的现象，以此为基础对扰动爆轰波的传播机理进行了分析，深化了对扰动来流中一维爆轰波传播机理的认识。

1 数值方法

1.1 控制方程

控制方程采用的是忽略黏性的一维Euler方程

假设气体为具有固定比热比的理想气体，总能量和状态方程为

式中，ρ,u,p,e,q,T分别代表流体的密度、速度、压力、总能量、化学反应放热量和气体温度。化学反应放热量q的表达式为

q=ηQ

其中η表示放热反应进程变量。上述所有变量均使用波前未反应气体状态参数进行无量纲化

其中，下角标0表示波前的气体参数。上述方程与化学动力学模型耦合，用来描述爆轰波的结构。采用Ng等[23]在Short等[24]基础上修正的两步链分支反应模型，用两个化学速率控制方程模拟化学反应过程。第1步是热中性诱导区域点火过程，反应速率是对温度敏感的Arrhenius形式

其中，ξ表示诱导区反应过程的变量，其变化范围是1 → 0;诱导区定义为爆轰波前导激波与最大释热速率位置之间的距离，对于两步诱导-放热总包反应模型而言，化学反应进程变量ξ=0对应着最大释热速率的位置;kI表示诱导区的化学反应速率常数，令kI=-Uvn，Uvn是以一维CJ爆轰波为坐标系前导激波后气流的流动速度;EI表示诱导区活化能;Ts表示气体经过一维CJ正爆轰波前导激波之后的温度;H(1-ξ)是阶跃函数

第2步是缓慢放热之后的快速热量释放过程，其反应速率方程如下

其中，η是放热反应进程变量，其变化范围是0 → 1;kR表示放热区反应速率常数;ER表示放热区活化能。上述公式中的EI和ER都通过RT0无量纲化，Ts通过T0无量纲化。

1.2 计算模型及参数

本文采用基于Mach数分裂的AUSMPW+格式[25]求解通量，空间差分采用的是具有3阶精度的MUSCL格式，时间推进采用3阶Runge-Kutta方法。计算域网格尺寸为无量纲的 0.05，诱导区宽度内的平均网格密度为20。将稳态Zel′dovich-von Neumann-Döring(ZND)爆轰参数作为初始点火源设置在计算域最左侧，计算开始后爆轰波自左向右传播。模拟的总体长度不同算例有差别，最少为100倍扰动波长距离，以保证压力峰值统计结果的可靠性。为提高计算效率，采用了移动网格技术，计算域随着爆轰波面不断前移。计算域的长度取决于扰动波长和诱导区长度，保证不小于扰动波长且至少达到200倍的诱导区长度。根据计算经验，此长度能够有效避免动网格后边界截断对爆轰传播的影响。为了避免初始数值振荡的影响，爆轰波首先在定常流动中稳定传播一段时间之后，再施加连续的密度扰动，其扰动形式为

其中A和λ分别表示扰动的振幅和波长，振幅被设置为波前静止气体压力的10%，波长是利用诱导区无量纲化后的长度。计算采用的无量纲参数为:放热量Q=50，比热比γ=1.2，诱导区活化能EI=6.0Ts，放热区活化能ER=1.0Ts。放热反应速率常数kR作为可调参数，用来控制爆轰波的内在不稳定性。

2 结果与讨论

2.1 无扰动来流中的一维爆轰波

可燃气体经过爆轰波前导激波面的绝热压缩，压力和温度迅速升高，并引发快速的化学反应。反应过程包含诱导阶段和紧随其后的放热阶段。诱导阶段的压力和温度基本保持不变，放热阶段释放能量的同时驱使反应产物膨胀，表现为温度上升且压力下降。采用ZND模型可以更直观地表达上述过程。图1显示了利用两步诱导-放热反应模型获得的稳态ZND结构。爆轰波设置在流场最左侧，并自左向右传播，跨过激波面的压力和温度瞬间上升至von Neumann状态，随后进入诱导反应阶段，诱导结束后开始放热反应，压力逐渐下降，温度继续上升至CJ状态。两步反应模型中的放热反应速率由参数kR控制，从图1中可以看出，放热阶段的温度和压力变化会随kR的增大加剧，然而并不影响其最终状态。对于一维爆轰波而言，kR只控制放热进程中温度和压力等状态参数变化的快慢，对一维爆轰波的静态ZND结构影响较小。

图1 稳态ZND爆轰结构Fig.1 Steady ZND detonation structure

虽然一维爆轰的ZND结构受kR的影响很小，但其动态传播特性会随kR的调整发生显著变化。传播特性可以通过记录前导激波后的压力变化来描述，具体表现为压力在传播方向上的振荡。研究者根据不同的振荡模态，将一维爆轰波分为稳定爆轰波和不稳定爆轰波，后者也被称为脉冲爆轰波。稳定爆轰波的激波后压力在传播过程中不发生变化，而在不稳定爆轰波中，压力则表现为一定范围内的振荡，且随着kR的变化呈现出不同的振荡模态。图2显示了kR对前导激波后压力变化过程的影响。波后压力在kR= 1.4时始终保持von Neumann状态不变，此时为稳定爆轰波。将参数调整到kR=1.8后，转变为脉冲爆轰波，其特征是波后压力具有单一振幅的周期性振荡。继续增大kR，则演化成更加复杂的振荡模态，单个周期内有两个不同的压力极大值，即所谓的双周期振荡。进一步增加kR，压力振荡的振幅明显增加且周期性消失，具有高度不规则的特征，这意味着爆轰波在传播过程中的振荡行为变得不可预测。

(a) kR=1.4

利用相空间方法，可以描述压力振荡状态的演变过程。图3展示了根据图2计算结果绘制的脉冲爆轰波的平面相轨迹，它可以反映传播过程中前导激波后压力的变化趋势。kR=1.8时，压力的变化具有固定的振幅和循环周期。kR=1.9时，压力变化分裂成两个不同振幅的循环周期，且两个周期的长度不一致。kR=2.1时压力的变化轨迹变得复杂，无显著的周期性。

(a) kR=1.8，1.9

上述计算结果表明，在给定其他两步反应参数的情况下，随着kR逐渐增大，前导激波后压力的变化会从定值过渡为周期性振荡，最终演变成高度不规则振荡。为了更详细地反映演化的过程，进一步计算了不同kR参数的一维爆轰波，并对各参数下的压力峰值进行了统计，绘制出局部压力峰值随kR变化的分岔图谱，如图4所示。其中，局部峰值压力用Pmax表示。该图谱可以明确展示各振荡模态相应kR的取值区间。图中的点表示取不同kR时统计的压力峰值。一组kR参数内，出现两个点则表示双周期振荡，多个点即多周期或不规则振荡。

图4 压力峰值随kR的变化Fig.4 Evolutions of local maximum pressure with kR

在放热反应速率较低的区间内(kR≤ 1.45)，前导激波后压力在传播过程中不发生变化，始终等于von Neumann状态的压力(Pvn=42.1)，属于稳定模态。放热速率加快(1.45

2.2 一维稳定爆轰的扰动振荡特性

根据图4的统计结果，认为kR<1.45都属于稳定爆轰波。为了研究稳定爆轰在扰动影响下的振荡特性，首先模拟了kR=1.0时一维稳定爆轰波在连续变化密度扰动中的传播，主要关注扰动波长对爆轰传播动力学过程的影响，这种影响通过传播过程中压力振荡模态的变化来表达。为了反映振荡模态与扰动波长之间的关系，采用了与图4相同的方法，对压力峰值随扰动波长的变化进行了统计，结果如图5所示。扰动波长λ=0表示未施加扰动，此时稳定爆轰的波面压力不发生变化。根据振荡模态的演化情况，将扰动对爆轰动力学的影响过程分为以下几个阶段:无影响阶段(λ<90)，表现为正弦扰动引起的单周期振荡。调整波长只会略微改变振幅，并不会引起振荡模态的变化。此过程与图4中kR在1.45～1.86区间内的变化趋势相似;相互作用阶段(90≤λ<1 000)，原始的稳定爆轰波开始响应扰动波，因此表现出相对丰富的现象。扰动与爆轰波的相互作用过程，具有以下两个特征:(1)随着λ增大，振荡的周期数先经过倍周期分岔，后逐渐回归单周期振荡模态。此外，在λ=800附近，压力振荡的波峰处出现短暂的小幅高频振荡，因此在统计图谱中出现了两个峰值。小幅振荡并不影响爆轰波的振荡模态。(2)伴随扰动波长的增加，峰值压力的最大值经历了线性上升和下降;扰动支配阶段(λ≥1 000)，此阶段振荡模态被锁定且振幅不再随λ变化，扰动对爆轰传播动力学的影响失效。计算获得了与单步反应不同的结果，发现对于稳定的爆轰波，在特定波长范围内，扰动也可以导致波后压力振荡。

图5 扰动波长对局部压力峰值的影响Fig.5 Influence of perturbation wavelength on local maximum pressure

为进一步观察稳定爆轰波的参数区间内，扰动对爆轰波动力学过程的影响趋势是否一致，统计了kR=1.4时压力峰值随扰动波长的变化，如图6所示。扰动波长λ<100时，激波后压力只存在一个峰值，此时为单周期振荡;λ增大，发现压力峰值的数量显著增多，说明波后压力转变为不规则振荡。随着λ进一步增大，压力振幅越来越大，峰值压力的最大值线性升高。当λ>550之后，振幅开始下降，最终演化为扰动支配的单周期振荡。调整扰动波长，压力振荡的整体演化趋势与kR=1.0时基本一致，但很显然kR=1.4时扰动对爆轰动力学过程的影响更显著。结果表明，扰动对一维稳定爆轰波动力学过程的影响，与放热反应速率有关。反应速率越快，爆轰波对扰动的响应程度越大。

图6 扰动波长对局部压力峰值的影响Fig.6 Influence of perturbation wavelength on local maximum pressure

利用功率谱密度(power spectral density，PSD)方法，可以更直观地反映整体演化趋势的一致性。图7展示了kR分别在1.0和1.4条件下，扰动波长λ=40，400，800对应的功率谱密度。其中，λ分别为40和800时，有一个明显的主频信号，说明能量都集中在这个频率上，因此导致振荡具有很强的规则性，即只有一个峰值的单周期振荡。需要注意的是，在λ=800的情况下，除主频信号外，还存在一些倍数于基础频率的谐波，但这并不表示振荡有多个周期。kR=1.0，λ=400时，主频信号数量增加，此时为多周期振荡。相同扰动波长条件下将kR调整到1.4，没有明显的主频功率，这意味着能量在频谱范围内的分布相对更广。PSD结果也说明扰动对稳定爆轰波动力学过程的影响只在适当波长范围内生效。

图7 前导激波压力的功率谱密度Fig.7 Power spectral density of the leading shock pressure

与kR=1.0不同的是，在kR=1.4时峰值压力的最大值并未线性下降，而是在λ>600之后出现振荡幅度的突然下降，这说明当扰动波长大于该值后，波前扰动对稳定爆轰波的影响会被大幅削弱。图8展示了扰动波长分别为600和650时，波后压力的振荡历程。可以看到，λ=600时波面压力在振荡的过程中会间歇性下降至Pvn的一半后重新上升，说明爆轰波在传播的过程中，激波和化学反应面会在扰动的干预下经历短暂的解耦，且重新耦合后会随机出现压力的突跃。扰动波长增大至λ=650后，这种随机的压力突跃消失，因此在图6中表现为最大压力的突然下降。对图8的结果进行功率谱密度分析，如图9所示。可以看到λ=600时能量在频谱范围的分布相对更广。λ=650时频率在0.01 附近存在一个非常明显的主频信号，说明此时振荡的规则性更强。

(a) λ=600

(a) λ=600

2.3 扰动幅值对扰动振荡的影响

前文在给定扰动幅值A=10%的条件下，研究了扰动波长的影响，发现对于不同的放热反应速率，稳定爆轰波随扰动波长的演化过程呈现出不同的特征。本节初步研究了扰动幅值的影响，计算结果如图10所示。

(a) kR=1.0，A=20%

kR=1.0，A=20%的扰动图谱显示，扰动对爆轰波的影响与图5相比更加显著，局部压力峰值随λ的演化过程与图6相似。另外，注意到扰动波长从200增大到260时，不稳定性并没有增加，反而呈现出更稳定的状态。这与单步反应中获得的重新稳定区相似，但这种现象在两步反应模型中是否具有普遍性，尚待进一步验证。kR=1.0，A=5%时绘制的扰动图谱，相较于图6的结果，扰动对爆轰波动力学过程的影响被显著削弱，表现为响应区间缩小和振荡模态的周期性变化。上述结果说明，扰动幅值越大，对爆轰波在扰动来流中传播的影响越大。

2.4 一维不稳定爆轰的扰动振荡特性

现有的研究结果发现，稳定爆轰波和不稳定爆轰波对扰动的响应过程有明显的区别[21]。前文已经讨论了一维稳定爆轰的扰动振荡特性，观察到区别于单步反应模型的现象。接下来讨论扰动对一维不稳定爆轰波的影响。图11展示了给定扰动幅值A=10%的情况下，对不稳定爆轰波施加不同波长扰动后统计的压力变化结果。由图2可知，kR=1.8无扰动情况下，前导激波后压力的变化呈单一振幅的周期性振荡。施加扰动后，在爆轰波的内在不稳定性和扰动波的共同作用下，前导激波后压力波动增加，局部峰值压力分布在更宽的范围。扰动波长增大，压力的整体变化趋势与稳定爆轰的结果相似。区别在于持续增大扰动波长λ，压力振荡行为也不会转变为简单的周期振荡。kR=2.1时，压力随λ的变化趋势与低kR的结果不同。随着λ的增大，逐渐演化为高频振荡的过程中，随机出现一些振幅很高的低频振荡。

(a) kR=1.8

3 结论

本文采用两步诱导-放热反应模型分析了连续密度/温度扰动对一维爆轰波传播特性的影响，重点关注了扰动波长、扰动幅值对不同化学反应活性的稳态爆轰振荡特征的影响规律，并与连续扰动作用下单周期振荡和多周期振荡模态的不稳定爆轰波进行了对比。

对于稳定爆轰波，扰动波长会在一定范围内触发爆轰波的内在不稳定性，使前导激波后压力呈现出多模态振荡的特征，且振荡模态随着反应物化学反应活性的增加而趋于复杂。超过一定的波长(本文中波长临界值接近800～1 000)，爆轰波恢复为单峰振荡的模态。总体上，爆轰波的振荡峰值压力随着扰动波长的增加呈现先增加后减小的趋势。同时，观察到在某些扰动波长下，爆轰波的周期性解耦-再起爆过程被削弱，导致振荡峰值压力出现突降。

通过对扰动幅值影响的分析，发现扰动幅值越大越容易触发爆轰波的内在不稳定性。进一步对比扰动作用下单周期和多周期振荡模态的结果可知，其振荡模态更加复杂多变，外界施加的周期性扰动难以让其恢复到初始的振荡模态。因此，连续扰动作用下的爆轰波传播特性主要是由两方面决定:放热区/激波面的耦合关系，以及扰动幅值和波长。前者决定了爆轰波的基础振荡模态和稳定特性，后者通过对波前气体干扰影响燃烧-流动的耦合关系，实现对爆轰波传播模态的调整。

致谢本文得到了国家自然科学基金(11822202)资助。