Page 52 - 《爆炸与冲击》2026年第10期

第 46 卷 陆 乐，等： 活塞驱动的病态爆轰波纵向稳定性研究 第 10 期

Å ã
Y A Y B (1−δλ)RT
pv = + RT = (8)
W A W B W A
式中：Y 和 A Y 分别为组分 A 和 B 的质量分数。注意式（8）中的第 2 个等式使用了组分质量分数和反应
B
度 λ 之间的关系：Y =1−λ，Y =λ。因此，根据式（7）～（8），比热比 Γ(λ) 的表达式可以写为：
B
A
c pA +λ(c pB −c pA )
c p
Γ(λ) = = (9)
c V c pA +λ(c pB −c pA )−(1−δλ)R/W A
式中：c 为等容比热。由式（9）可以看出，化学反应摩尔数减少的二元体系中，比热比不仅是 λ 的函数，
V
也是 δ 的函数。在本文中，固定 δ 不变（δ=0.5），因此在化学反应进行的过程中，比热比 Γ 仅随着反应的
进行而发生变化。
对式 (1) 中的对流通量项，采用 Lax-Friedrichs 通量分裂方法与五阶 WENO 格式 [20] 耦合的离散方法
进行求解。对化学源项 S，采用变阶变步长的 VODE 程序包求解。关于时间推进项，采用了显式三阶
Runge-Kutta 时间离散格式，设置库朗数（Courant Friedrichs Lewy, CFL）为 0.4，以确保对流-反应耦合系统
的时间精度与计算鲁棒性。
采 用 实 验 室 坐 标 系 建 立 计 算 模 型 ， 设 置 一 维 计 算 域 长 度 为 1 000Δ 1/2 ， 这 里 Δ 1/ 2 为 相 应 的 爆 轰 波
ZND 结构中的半反应区长度，即反应完成一半时的反应区长度。初始时刻，在计算域内的左侧置入给定
条件下预先计算得到的 ZND 结构的流场，设置其长度为 10Δ ，以保证 ZND 结构下游已基本达到反应结
1/2
束的状态。此时，ZND 结构的下游状态便可作为计算域左端边界条件，而计算域的右端则设定为零梯度
边界条件。数值模拟时，爆轰波从零时刻的稳态 ZND 结构出发，从左向右运动，在数值误差的扰动下爆
轰波很快变为非线性失稳状态，通过记录不同时刻下爆轰波阵面压力峰值 p ma x 来考察爆轰波的一维非
线性纵向不稳定特性。
1.2 病态爆轰的 ZND 模型和计算方法
在采用上述数值模拟方法求解反应摩尔数减小的病态爆轰的非线性稳定性特性之前，首先需要对
其 ZND 结构进行求解，以便作为数值模拟的初始条件。求解 ZND 结构的控制方程组是式 (1) 的稳态
版，即不考虑时间推进项（∂U/∂t）的影响，因此控制方程变为常微分方程组的形式，对稳态形式的式 (1) 的
前 3 个式子进行重排，并采用激波坐标系，则可以得到如下的形式：

dρ ρΣ
 = −

 dt η



du uΣ

= (10)
 dt η


2
 dp ρu Σ


 = −
dt η
2
式中：u 为激波坐标系下的流体质点运动速度，η=1−Ma ，Ma 为流体质点的运动马赫数，可表达为 Ma=u/a，
这里 a 为流体质量的当地声速，a=[Γ(λ)p/ρ] 。Σ 代表反应体系的热性（thermicity），定义 [19] 为：
1/2
dλ
Σ = σ (11)
dt
式中：σ 为热性系数。对变摩尔数的化学反应，σ 可以表示为 [21] ：
α f q δ
σ = − (12)
c pA +λ(c pB −c pA ) 1−δλ
式中：α 为膨胀系数。对于量热完全气体，α =ρ(∂v/∂T) =1/T [19] 。假设反应前后系统的等压比热保持不
f
p,λ
f
变，即，c =c ，于是式 (12) 进一步简化为：
pA
pB
q δ
σ = − (13)
c pA T 1−δλ
102201-4