引 言
气相介质通常存在非均匀性, 区别于经典的气相爆轰理论, 非均匀气相介质的爆轰机理研究面临许多挑战性的基础理论问题. 研究非均匀气相介质的爆轰过程具有极为重要的实际应用价值, 可以为工业安全防护和爆轰推进设计提供重要的理论和技术支持.
在可燃气体泄露后的流场中, 非均匀性以各个热力学状态参数的空间分布梯度的形式体现, ****们研究了气体组分浓度梯度[1-3]或者温度梯度[4-5]对爆轰波波阵面和传播极限的影响. 当前对于非均匀气相介质爆轰波的研究主要集中于热力学状态参数的空间分布存在梯度这一模型, 进而研究梯度分布对爆轰波波阵面和传播极限的影响. 在可燃气体泄露后的流场和爆轰发动机的工作过程中, 存在着多方向的组分浓度梯度或者温度梯度. 国内外很多****对垂直于传播方向存在组分浓度梯度的气体混合物的爆轰传播过程, 特别是波阵面的演化过程, 进行了实验[6-8]和数值模拟[9-11]研究. Ishii和Kojima[6]研究了波前介质非均匀对爆轰波阵面形状和爆轰不稳定性的影响, 结果表明, 波阵面会形成一偏转角, 并且该角度随局部梯度的增加而增加. Ettner等[12]研究了直管中不同浓度梯度对爆轰波阵面、爆轰稳定性和压力分布等的影响, 结果表明, 壁面的约束使得波阵面出现弯曲或者马赫反射, 但不会出现规则反射. 蔺伟等[13]研究了浓度梯度对瓦斯爆炸影响的数值模拟, 以及爆轰波形成的胞格大小和三波点随梯度的变化规律, 发现由于浓度梯度的存在, 波阵面上会出现马赫反射现象, 同时研究了马赫反射对整个传播过程的影响. Han等[14]通过高精度数值模拟和一步反应模型研究了纵向存在组分梯度情况下爆轰波的波阵面结构、流场细节和传播极限问题.
同时很多****也对平行于爆轰波传播方向存在组分浓度梯度的情况进行了研究. Kuznetsov等[15-16]便对组分活性梯度降低的爆轰传播进行了类似的研究, 并研究了驱动段长度的影响. Thomas等[17]则发现小的组分浓度梯度有助于爆轰波的透射, 而大的梯度则可能使激波阵面和反应区分离, 导致爆轰波熄爆. 王春等[18]对爆轰波穿越惰性气团时的透射激波进行了参数分析. 他们发现当爆轰波穿越惰性气团再进入可燃混合气体时, 透射激波的参数主要受到3个基本过程的影响.
上述研究中, 热力学参数的非均匀性是来自整体空间分布上的梯度. 这种梯度分布模型虽然更接近实际的情况, 但是并不利于建立非均匀介质爆轰波的稳态传播机理. 另外一种思路则是局部上添加扰动, 爆轰波此时会以一种近似稳态的模式传播, 波阵面近似平面(虽然存在小尺度的三波结构). 例如, Mi等[19]通过在均匀介质中添加随机分布的离散单元实现非均匀性, 离散单元的存在使得爆轰波波阵面更加扭曲、离散, 同时也促进了爆轰波的传播, 使其不容易熄爆失效. Mi等[20-21]和Lam等[22]则利用一个一维热方程研究了具有无限刘易斯数和空间离散放热源的系统中火焰传播的动力学问题, 该方程的源项为反应过程变量, 该变量的扩散率为零, 模型揭示了火焰传播模式的谱系. Li等[23]通过大振幅、二维密度正弦波动引入非均匀性, 对薄炸药层中爆轰波传播极限附近的爆速进行了计算. 结果表明, 非均匀介质中的爆轰波具有与均匀介质中类似的胞格结构. 爆轰能够在一个非常薄的炸药层传播而不熄爆, 并且比相应的均匀介质呈现出更大地传播速度. Higgins[24]通过在化学反应速率模型的基础上乘以一个与当地坐标有关的周期函数来实现人工热点, 作为均匀反应速率的替代方法, 分析了在含能材料中离散效应对爆轰波传播的影响. 但是这些研究存在一个问题, 即他们使用的反应模型为压力敏感型方程, 这使得胞格结构不易显现出来, 这有利于研究波阵面曲率的影响, 但是忽略了胞格不稳定性和横波结构的影响, 这种假设更适用于凝聚相爆轰.
另外, 在爆轰发动机的工作过程中, ****针对于非均匀性对于爆轰发动机的影响也开展了相应的数值仿真[25]与实验[26]研究. 例如, 分析了不同的非预混过程对组分浓度分布的影响规律[27], 以及湍流效应在这一过程中的作用[28]. 此类研究通常不考虑点火以及其后的爆轰过程, 也有很多研究涉及发动机非预混爆轰的整个过程, 如点火过程[29]、波头演化[30], 进而研究非均匀性对发动机推力的影响[31].
研究非均匀气相介质的爆轰过程具有极为重要的实际应用价值, 可以为工业安全防护和爆轰推进系统设计提供重要的理论支持. 在非均匀性存在的前提下, 气相介质的爆炸过程, 包括直接起爆, 稳态传播和失效机制, 都极其复杂, 很多物理机制尚不明确. 当前的很多研究中, 由于实验手段的限制, 很难对周期性非均匀介质中气相爆轰波传播机制进行实验研究, 而数值模拟可以容易地实现这一目的, 便于对物理机制的理解. 基于此,本文使用反应欧拉方程和两步反应模型对气相爆轰波在非均匀介质中的传播过程进行数值模拟, 主要探究非均匀介质横向均匀温度扰动对气相爆轰波的传播机理进行系统的参数研究, 特别是非均匀性将以何种形式影响爆轰波的直接起爆、传播和失效以及波阵面多波精细结构, 这种影响又是否与爆轰波自身的不稳定性和动力学参数存在某种定量的关系, 将是本文研究的重点. 相对于复杂的基元反应模型, 两步法可以容易地实现爆轰参数的解耦, 方便系统地进行单一参数研究, 而又不失一般性.
1.
控制方程和数值方法
1.1
控制方程
在本文中控制方程使用耦合了两步化学反应模型的欧拉方程. 忽略了热传导、分子扩散以及黏性效应的影响, 控制方程表示为
$$ frac{{partial {{boldsymbol{U}}}}}{{partial t}} + frac{{partial {{boldsymbol{F}}}}}{{partial x}} + frac{{partial {{boldsymbol{G}}}}}{{partial y}} = {{boldsymbol{S}}} $$  |
其中为
$$ begin{gathered}{{boldsymbol{U}}} = left[ { ho, ho u, ho v,E, ho {y_{{ m{I}}}}, ho {y_{{ m{R}}}}} ight] hfill {{boldsymbol{F}}} = left[ { ho u, ho {u^2} + p, ho uv,u(E + p), ho u{y_{{ m{I}}}}, ho u{y_{{ m{R}}}}} ight] hfill end{gathered} $$  |
$$ begin{gathered}{{boldsymbol{G}}} = left[ { ho v, ho vu, ho {v^2} + p,v(E + p), ho v{y_{{ m{I}}}}, ho v{y_{{ m{R}}}}} ight] hfill {{boldsymbol{S}}} = left[ {0,0,0,0,{omega _{{ m{I}}}},{omega _{{ m{R}}}}} ight] hfill end{gathered} $$  |
对于理想气体, 总能和状态方程可表示为
$$ begin{gathered}E = frac{p}{{gamma - 1}} + frac{ ho }{2}left( {{u^2} + {v^2}} ight) - ho {y_{{ m{R}}}}Q hfill p = ho RT hfill end{gathered} $$  |
其中p为压力, ρ为密度, E为单位体积的能量, u和v分别为x和y方向上的速度, Q为反应热, yI与yR分别为两步反应模型中表征诱导和放热过程的进程参数.
1.2
两步化学反应模型
化学反应采用两步反应模型[32], 即先进行诱导过程, 诱导过程结束后再进行放热反应过程. 第一步诱导过程速率可表示为
$$ {omega _{{ m{I}}}} = - Hleft( {1 - {y_{{ m{I}}}}} ight){tilde k_{{ m{I}}}} ho exp left[ {{E_{{ m{I}}}}left( {frac{1}{{R{T_{{ m{s}}}}}} - frac{1}{{RT}}} ight)} ight] $$  |
其中
m{I}}}} $
m{s}}}} $
m{I}}}}
ight) $
$$ Hleft( {1 - {y_{{ m{I}}}}} ight) = left{ {begin{array}{*{20}{l}}{1,}&{;;;{y_{{ m{I}}}} < 1} {0,}&{;;;{y_{{ m{I}}}} geqslant 1} end{array}} ight. $$  |
诱导过程结束后, 第二阶段为放热反应阶段, 反应速率表示为
$$ {omega _{{ m{R}}}} = left[ {1 - Hleft( {1 - {y_{text{I}}}} ight)} ight] ho {tilde k_{{ m{R}}}}left( {1 - {y_{{ m{R}}}}} ight)exp left( {frac{{{E_{{ m{R}}}}}}{{RT}}} ight) $$  |
其中ER为反应阶段的活化能, 指前因子
m{I}}}} $
m{R}}}}$
m{I}}}} $
m{I}}}}$
m{R}}}}$
$$ {M_{{text{CJ}}}} = sqrt { {1 + dfrac{{{gamma ^2} - 1}}{gamma }Q} + sqrt {{{left( {1 + dfrac{{{gamma ^2} - 1}}{gamma }Q} ight)}^2} - 1} } $$  |
为了便于与文献比较, 这里定义新的活化能参数, 即
$$ begin{split} &{varepsilon _{{ m{I}}}} = frac{{{E_{{ m{I}}}}}}{{R{T_{{ m{S}}}}}} &{varepsilon _{{ m{R}}}} = frac{{{E_{{ m{R}}}}}}{{R{T_{{ m{S}}}}}} &{T_{{ m{S}}}} = frac{{(2gamma M_{{{ m{CJ}}}}^{2} - gamma + 1)[2 + (gamma - 1)M_{{{ m{CJ}}}}^{2}]}}{{{{(gamma + 1)}^2}M_{{{ m{CJ}}}}^{2}}} end{split}$$  |
一般来说, 对于常见的碳氢燃料, 活化能
m{I}}}}$
m{I}}}} gg {varepsilon _{{
m{R}}}}$
表
1
爆轰参数
Table
1.
Detonation parameters
table_type1 ">
| Parameter | Value | Unit |
| R | 218.79 | J/kg· K |
| p0 | 50 | kPa |
| T0 | 295 | K |
| c0 | 304.86 | m/s |
| ρ0 | 0.775 | kg/m3 |
| Q/(RT0) | 19.7 | ? |
| γ | 1.44 | ? |
| MCJ | 5.6 | ? |
| ${varepsilon _{ m{I}}}$ | 4.8 | ? |
| ${varepsilon _{ m{R}}}$ | 1.0 | ? |
| kI = ${tilde k_{{ m{I}}}}{x_{{text{ref}}}}/{c_0}$ | 1.3875 | ? |
| kR = ${tilde k_{text{R}}}{x_{{text{ref}}}}/{c_0}$ | 2.0~5.0 | ? |
下载: 导出CSV
|显示表格
1.3
计算方法
为了分辨爆轰波阵面精细的多波结构, 必须使用足够高的网格分辨率. 为了满足这一要求, 在本研究中使用了基于块结构化自适应网格方法(SAMR)的AMROC程序[33-34]. AMROC是一个基于有限体积法开发的开源框架, 广泛应用于燃烧、爆轰和空气动力学的模拟研究, 而且用户可以在该框架下嵌入自定义的算法和反应模型. 本研究在AMROC中耦合了两步化学反应模型, 使其可以进行相关的爆轰研究. 在整个计算区域, 使用粗网格, 在流场剧烈变化的区域使用细网格, 使之具有更高的网格分辨率, 以减少计算时间和内存需求. 针对反应欧拉方法, 使用两阶Strang分裂方法解耦对流项和化学反应源项, 然后采用基于Roe-HLL求解器的显式二阶MUSCL-Hancock方法求解对流项. 反应源项采用半隐式的4阶通用Runge-Kutta刚性常微分求解器(GRK4A)[35]进行数值积分.
1.4
初始设定
在本研究中, 爆轰波的比热比
m{I}}}} $
m{R}}}}$
m{I}}}}{varDelta _{{
m{I}}}}/{varDelta _{{
m{R}}}}$
本研究针对非均匀介质中爆轰波发展分别进行了一维和二维数值模拟. 一维计算域设定如图1(a)所示, 左侧为一ZND爆轰波, 右侧为温度T1的热区或温度为T2的冷区. 当待起爆区为热区时, 表达式为T1 = T0+A, 为冷区时, 表达式为T2 = T0?A, T0为基准温度295 K, A为温度幅值. 二维周期性扰动初始设定则如图1(b)所示, 横向长度为
m{I}}}}$
m{I}}}}$
m{I}}}}$
m{I}}}}$
m{I}}}}$
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-1.jpg'" class="figure_img
figure_type1 bbb " id="Figure1" />
图
1
计算域设置
Figure
1.
Setup of the simulation domain
下载: 全尺寸图片
幻灯片
1.5
收敛性测试
如图2所示, 首先对于一维爆轰波的计算进行网格收敛性测试. 相同初始条件下, 对于5种不同的网格分辨率进行了对比, 分别为4pts/ΔI, 8pts/ΔI, 16pts/ΔI, 32pts/ΔI和64pts/ΔI. 可以从图2中清晰地观察到, 无论是爆轰波传播进冷区或热区, 在网格分辨率为32pts/ΔI时均收敛, 因此本文一维算例网格分辨率均为32pts/ΔI. 之后对于二维ZND爆轰波向胞格爆轰波的转变过程进行了测试, 结果如图3所示, 网格分辨率分别为16pts/ΔI, 32pts/ΔI和64pts/ΔI. 对比可得, 随着网格分辨率的增加, 可以捕捉到更为清晰的爆轰波胞格结构, 与分辨率为64pts/ΔI结果相比较, 分辨率32pts/ΔI下的计算结果并未有太大不同, 考虑到计算成本, 本文针对于二维爆轰波的数值模拟, 网格分辨率采用32pts/ΔI.
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-2.jpg'" class="figure_img
figure_type1 bbb " id="Figure2" />
图
2
一维爆轰波网格分辨率测试
Figure
2.
Numerical resolution test for one-dimensional detonations
下载: 全尺寸图片
幻灯片
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-3.jpg'" class="figure_img
figure_type1 bbb " id="Figure3" />
图
3
二维胞格爆轰波网格分辨率测试(无扰动, kR = 2.0)
Figure
3.
Numerical resolution test for transition of ZND to cellular detonations (without disturbance, kR = 2.0)
下载: 全尺寸图片
幻灯片
2.
结果与讨论
为了研究爆轰波在非均匀介质中的传播问题, 做了以下3方面的工作. 首先, 通过对平面CJ爆轰波进入“热”或“冷”区转变过程进行简单的一维数值模拟研究, 揭示了爆轰波在不同类型扰动下的传播机理, 特别是在 “冷”区内的失效极限和重新起爆过程进行了分析. 其次, 将研究拓展到二维, 揭示了爆轰波阵面结构在存在间断温度扰动时的结构和演化过程, 在这一部分中将温度扰动限定为一个周期. 最后, 将扰动设定为多个周期, 并且改变扰动的振幅, 讨论了爆轰波在温度扰动下的传播过程和波阵面结构. 通过从简单到复杂的研究思路, 希望可以更清晰地展示爆轰波波阵面结构的变化, 更容易解释传播机理.
2.1
一维爆轰波转变研究
首先对于平面CJ爆轰波进入“热”或“冷”区进行了一维研究, 讨论了爆轰波的转变过程和随后的传播过程. 计算设定如图1(a)所示, 左侧为一初始温度为T0的CJ爆轰波, 右侧基准温度为T0, 通过添加或减小一个温度幅值变化A表征热区和冷区.
具有不同稳定性系数(kR = 2.0, kR = 3.0, kR = 4.0)的爆轰波的一维计算结果如图4(a) ~ 图4(c)所示, 最终的稳态爆轰波波速如图4(d)所示. 稳态传播的CJ爆轰波遇到新的温度界面时, 由于两侧的热力学性质不同, 会在界面处产生右行的透射波和左行的反射波. 研究发现, 平面CJ爆轰波在遇到初始温度较高的界面时, 会产生一个压力较低但速度较高的透射波向前传播, 如图4(d)所示, 与此同时, 会产生一个反向传播的反射波. 而当平面CJ爆轰波遇到初始温度较低的界面时, 会产生与热界面完全相反的现象, 即透射波的压力增加, 速度降低.
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-4.jpg'" class="figure_img
figure_type1 bbb " id="Figure4" />
图
4
压力历史和速度曲线图
Figure
4.
Pressure history curve
下载: 全尺寸图片
幻灯片
爆轰波完成转变过程后, 根据传播界面温度的不同, 透射爆轰波会以不同的稳态速度和压力传播. 另外, 计算结果表明, 热界面温度大大降低了爆轰波稳态传播时的峰值压力, 略微加快了爆轰波传播速度. 与此相反的是, 冷界面增大了爆轰波稳态传播的峰值压力, 减慢了爆轰波的传播速度. 然而, 随着冷界面温度的继续下降, 爆轰波会经历一个熄爆再起爆的过程, 即为临界状态, 如图4(c)中冷区温度为280 K时曲线所示. 当温度再次降低时, 爆轰波会完全熄爆, 如图4(a) ~ 图4(c)中温度为190 K曲线所示. 在一维计算中, 冷界面的临界温度在240 K左右, 当低于此温度时, 爆轰波会熄爆.
对于不稳定爆轰波和稳定爆轰波, 前者的化学反应率对温度更为敏感, 在遇到冷区扰动时, 反应率下降更快, 更容易衰减熄爆. 图4(d)中色块表示透射过程中出现熄爆和重新起爆的范围. 可以看到, 随着冷区温度的下降, 不稳定爆轰波(kR = 4.0)熄爆更早出现. 但是对3种爆轰波而言, 完全熄爆的冷区温度基本都在240 K左右, 没有明显的区别. 这与实验中的多维爆轰波的结论有区别, 原因主要是一维爆轰波和多维爆轰波传播机制的不同.
2.2
二维激波/爆轰波单周期转变过程研究
在这一部分中, 首先使用激波然后再使用爆轰波作为入射波传入非均匀介质中, 通过这样的对比研究, 可以更容易的对爆轰波的传播机制进行研究. 在二维情况下, 非均匀区域由中心的热区域和上下边界处两个冷区组成的阶梯分布构成, 温度扰动设定为一个半周期. 与本文的一维研究的类似, 激波透射后发展到稳定态经历了一个转变过程. 不同之处在于, 由于界面两侧的初始温度不同, 导致流动速度不同, 从而在界面处产生剪切层, 因此在界面附近处可以观察到Kevin-Helmholtz不稳定性. 由于本研究中使用了无扩散效应的反应欧拉方程, 因此流动中显示的漩涡并不是真正的物理黏性现象, 而是由于数值黏性造成的结果. 尽管如此, 如Han等[14]所分析的, 目前所使用的方法在某种程度上依然能够体现Navier-Stokes方程的某些特性.
如图5(a) ~ 图5(f)所示, 随着激波在横向非均匀介质中的传播, 激波前沿演化成一个复杂的结构. 在温度较高的中心区域存在一轻微前凸的激波波阵面, 而在上下温度较低的冷区, 可以发现由入射波、马赫杆、反射激波和剪切层组成的三波结构. 入射激波与中心弯曲的激波面相连, 并被剪切层隔开. 因此, 存在两个剪切层, 即由温度界面演化成的剪切层和三波结构造成的马赫反射剪切层, 而且这两个剪切层在斜入射波后形成一个薄通道. 激波面后的流场由三部分组成, 第一部分为中心处前凸激波面后的高温低压膨胀区, 第二部分为斜入射激波后两剪切层之间的狭窄通道, 最后一部分为边界附近两正入射激波后的高压低温区.
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-5.jpg'" class="figure_img
figure_type1 bbb " id="Figure5" />
图
5
二维激波波阵面演化图, 温度扰动振幅A = 105 K
Figure
5.
Evolution of shock front structures during the two-dimensional with A = 105 K
下载: 全尺寸图片
幻灯片
图5(e)的局部放大图中给出了波阵面的精细结构以及流线分布, 设置流线时将坐标系建在了波阵面上. 斜入射激波的存在使得流线向边界处弯曲, 导致中心区流场的膨胀以及边界处流场的收缩. 相邻的两对剪切层在激波面后一定距离处合二为一, 在靠近上边界的另一对剪切层沿与中心线方向相反的方向卷曲. 随着涡流尺寸的增大, 两剪切层在距离激波面后较远距离处重合, 并且在非均匀介质中形成了独特的激波复合结构体.
图5(g)给出了最大压力历史云图(对应于爆轰波模拟中的数值胞格图), 可以看出两侧临边界处压力大, 中心压力低, 但是中心低压区(初始高温区)面积扩张, 两侧高压区(初始低温区)面积收缩.
图6比较了不同温度扰动幅值下的激波波阵面的结构特征及其压力历史云图, 其中扰动幅度A分别为15 K, 35 K和105 K. 可以发现在温度扰动幅度较小时并没有典型的三波点结构生成, 而是在中心区域产生了一个微凸的激波以及上下边界两个微凹的激波. 这是因为较低的温度幅度使得中心与上下边界流动状态相差不大, 这也延缓了剪切层中涡的演化. 从波阵面上的最大压力历史中可以发现不同扰动幅度下, 马赫杆的高度和微凹激波的高度发生变化, 温度扰动幅值越大, 马赫杆越短, 幅值越小, 马赫杆越长. 同时发现较大幅度的温度扰动下, 波后湍流产生的压缩波使得中心区域出现较大的压力震荡.
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-6.jpg'" class="figure_img
figure_type1 bbb " id="Figure6" />
图
6
不同温度扰动幅值下激波波阵面纹影图和最大压力历史
Figure
6.
Schlieren photography of the shock front at the same instant for cases with different temperature perturbations and pressure history
下载: 全尺寸图片
幻灯片
下面使用ZND爆轰波作为入射波传入存在单周期温度扰动的非均匀介质中. 与激波相比, 因为透射和传播过程涉及到化学反应, 使得问题更为复杂. 图7给出了扰动幅值A为105 K时的爆轰波前沿的演化过程, 其中kR = 2.0. 从图中可以观测到, 在上下边界冷区内, 爆轰波波阵面出现一定程度的解耦现象, 即诱导区长度不断变长. 此外, 在解耦区出现不断试图再次起爆的热点, 但均以失败告终. 即使在中心区域处前导激波与反应面也稍有解耦. 在未有化学反应的冷区也发现了三波结构. 由于冷区存在周期性的熄爆与再起爆过程以及流场的湍流特性, 产生了未反应的区域. 而此区域的燃烧产生了不同强度的压缩波. 可以发现这些压缩波的相互作用可能在冷区产生局部的热点, 即使此时冷区初始温度为190 K. 联想到一维数值模拟中冷区温度为190 K时, 透射爆轰波最终熄爆, 可以说明一维和二维的情况存在不同, 这也表明了第二个维度上的梯度对于此问题的影响.
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-7.jpg'" class="figure_img
figure_type1 bbb " id="Figure7" />
图
7
二维爆轰波波阵面演化图, 温度扰动振幅A = 105 K, kR = 2.0
Figure
7.
Evolution of detonation front structures during the two-dimensional with A = 105 K and kR = 2.0
下载: 全尺寸图片
幻灯片
图8给出了两个时刻沿中心线和边界处压力与温度的分布曲线. 与激波的例子相比, 爆轰波波阵面由于化学反应的影响, 从未到达与激波波阵面相同的稳态, 而是一个不稳定的演化过程. 沿中心线和边界处压力与温度历史曲线表明, 冷区波阵面前沿的压力高于中心波阵面前沿的压力. 图8(a)中压力曲线的第二个峰值清晰地显示了重新起爆(热点)的位置(t =0.25 ms). 观察温度曲线图, 可以发现在x = 0.25 ~ 0.35 m之间温度出现高频率的震荡, 表明在涡流中存在大量未反应完全的气团. 这是未完全反应可燃气体沿着上下壁面通道被涡卷入波后的中部区域造成的.
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-8.jpg'" class="figure_img
figure_type1 bbb " id="Figure8" />
图
8
不同时刻下沿中心线和边界线压力与温度分布和流线图, 温度扰动振幅A = 105 K, kR = 2.0
Figure
8.
Pressure, temperature distribution and streamlines along the centering and boundary lines at two instants with A = 105 K, kR = 2.0
下载: 全尺寸图片
幻灯片
而后研究了不同稳定性下的爆轰波演变过程, 其中kR分别为2.0, 3.0和4.0, 温度幅值分别为15 K, 35 K和105 K, 其纹影图与数值胞格图如图9所示. 从图中可以清晰地观测到均匀介质中典型的二维胞胳爆轰波结构, 唯一的区别在于由于横向温度分布不均产生的剪切层. 由于计算域的宽度仅为诱导区长度的30倍左右, 所以横向只容许存在一个胞胳. 另外, 爆轰波的不稳定性仅影响了胞胳爆轰波的相位, 并未显著改变胞胳的大小. 如图9(a)所示, 在温度扰动幅值较小时(15 K), 可以清楚的观察到规则的胞格结构, 温度扰动并没有对胞格形状产生明显的影响. 如图9(b)所示, 在温度扰动幅值增大到35 K时, 可以看出胞格的形状发生了一定程度的改变, 主要体现在长宽比的增加, 这说明横波的速度相对于前导波阵面的速度减小. 而且在图9(b)的数值烟膜中, 爆轰波存在一定程度的解耦?失效, 然后重新起爆过程, 这体现了爆轰不稳定性的影响. 如图9(c)所示, 在温度扰动幅值增大到105 K时, 可以看出数值烟膜中不存在典型的胞格结构, 这说明横波发展受到了温度扰动的抑制, 波阵面上不断变化的爆轰波三波结构基本消失, 取而代之的是与激波例子中类似地稳态的波阵面结构, 只有在中心处能够看到非常弱的横波产生的痕迹. 同时也可以看到, 由于边界两侧存在熄爆和重新起爆不断重复的过程, 马赫杆的长度也发生了周期性的变化.
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-9.jpg'" class="figure_img
figure_type1 bbb " id="Figure9" />
图
9
不同稳定性下爆轰波波阵面纹影、温度扰动幅度以及胞格结构
Figure
9.
Schlieren photography of the detonation front at the same instant for cases with different instability
下载: 全尺寸图片
幻灯片
2.3
非均匀介质中波阵面结构
当在二维管道中布置一维ZND爆轰波作为初始条件时, 均匀介质中平面爆轰波的形态在向前传播的过程中能够维持一段时间或者距离. 但是数值扰动会诱发波阵面局部上小尺度的扭曲或者褶皱, 并在爆轰不稳定性的作用下, 波阵面的扭曲放大, 进而在横向上产生压力、温度和密度的梯度, 导致弱横波的产生. 弱横波放大受不稳定性因子的控制, 最终会达到平衡态, 即横波的强度以及横波的间距相对保持不变, 而且该平衡态取决于可燃介质的热力学状态和化学动力学特性, 这也说明横波间距是一个特征尺寸, 此时便形成了一个完全发展的自维持胞格爆轰波. 图10给出了
m{R}}}} = 2 $
m{I}}}} $
m{I}}}} $
m{R}}}} $
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-10.jpg'" class="figure_img
figure_type1 bbb " id="Figure10" />
图
10
均匀介质中无扰动时ZND爆轰波到胞格爆轰波的演化
Figure
10.
Transition from ZND to cellular detonation without artificial perturbations
下载: 全尺寸图片
幻灯片
图11 ~ 图13给出了不同扰动波长s情况下ZND爆轰波演化发展过程的数值胞格. 图11对应的压力历史曲线(取自中心线)如图14所示. 这一部分的研究中使用了小振幅的温度扰动, 即A = 15 K. 从图11中可以看出, 引入温度扰动后, 只经过了几个诱导区的距离波阵面上便出现了扭曲变形, 并形成非常弱且非常规则的横波, 这可以从图中规则的小胞格中看出, 也可以从图14中小振幅震荡的规则压力曲线中反应出来. 图11中规则小胞格的尺寸几乎等于扰动的波长10ΔI. 但是, 经过了一段传播距离达到位置x = 400ΔI时, 横向不稳定性开始放大增长, 最终在x = 700ΔI处达到平衡状态, 形成了充分发展的稳态胞格爆轰波. 压力历史曲线的震荡也开始不规则, 并且震荡增加, 与无扰动时压力的震荡特性相一致(如图10(b)所示). 这说明, 人工温度扰动可以在一定的范围内抑制胞格不稳定性的发展, 但是并不能够终止这一进程. 这主要是因为人工扰动的波长(s = 10ΔI)远小于特征横波间距(40ΔI ~ 50ΔI), 并且人工扰动的振幅也较小.
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-11.jpg'" class="figure_img
figure_type1 bbb " id="Figure11" />
图
11
存在人工温度扰动时的数值胞格, A = 15 K, s = 10ΔI
Figure
11.
Numerical smoked foils for cases with artificial perturbations of A = 15 K, s = 10ΔI
下载: 全尺寸图片
幻灯片
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-12.jpg'" class="figure_img
figure_type1 bbb " id="Figure12" />
图
12
存在人工温度扰动时的数值胞格, A = 15 K, s = 20ΔI
Figure
12.
Numerical smoked foils for cases with artificial perturbations of A = 15 K, s = 20ΔI
下载: 全尺寸图片
幻灯片
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-13.jpg'" class="figure_img
figure_type1 bbb " id="Figure13" />
图
13
存在人工温度扰动时的数值胞格, A = 15 K, s = 40ΔI
Figure
13.
Numerical smoked foils for cases with artificial perturbations of A = 15 K, s = 40ΔI
下载: 全尺寸图片
幻灯片
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-14.jpg'" class="figure_img
figure_type1 bbb " id="Figure14" />
图
14
存在人工温度扰动时爆轰波的压力历史曲线
Figure
14.
Pressure history for cases with artificial perturbations
下载: 全尺寸图片
幻灯片
如图12所示, 当人工温度扰动的波长增大至s = 20
从上述的分析中也可以发现, 随着不稳定性参数的增大, 上述ZND爆轰波的发展演化过程可被大大提前. 因此总结来说, ZND爆轰波在人工温度扰动下的响应(发展演化过程)主要受制于两个竞争性的因素: 一是爆轰波内在的不稳定性(主要受热力学和化学动力学特性控制), 二是人工扰动的波长, 前者是内因, 后者是外因. 但是, 由于上述算例中, 人工温度扰动均设置了相同的、相对较小的振幅(15 K), 即便存在温度扰动, 最终也能够得到一个充分发展的胞格爆轰波, 并且这些受扰动胞格爆轰波的主要特征与在均匀介质中的未受扰动胞格爆轰波基本一致.
由于人工扰动的存在, 压力震荡的发散程度减小, 但是频率增大. 人工温度扰动可以看作是横向周期性排列的热点和冷点区域. 热点区域能够通过增大反应速率而增强爆轰波的不稳定性, 其中冷点通过减小反应速度而限制爆轰波的发展. 但是冷点的负面效应要强于热点的正面效应, 这是由于阿伦尼乌斯反应率与温度不是线性关系. 因此人工温度扰动的存在, 一方面能够限制波阵面的发展, 另一方面也能够限制爆轰波横波的发展, 使得波阵面更加的离散, 压力震荡的波动幅度减小.
上述的研究中, 重点关注较小温度扰动幅度下ZND爆轰波的发展演化过程. 下面研究在较大温度扰动时, 波阵面的结构的变化. 图15中给出了4种不同温度扰动时ZND爆轰转变为胞格爆轰波时的数值胞格结构, 管道中线处的压力历史曲线见图16. 在x = 1.2 m处波阵面结构及波后的压力分布如图17所示. 从数值胞格中可以看出, 随着温度扰动幅度的增加, 胞格的形态发生了很大的变化, 当A = 15 K时, 温度扰动的影响非常有限, 胞格的形状与无扰动时基本一致. 图17(a)中能够看到扰动使得波阵面结构局部上发生细微的影响, 波阵面存在小的间断, 但是并没有对波阵面整体的结构产生大的影响, 横波依然很清晰, 压力曲线的震荡也在正常的范围内. 但在A = 35 K和A = 55 K时, 数值胞格的形态发生了较大的变化, 三波的轨迹局部上发生了较大的扭曲, 这说明横波在穿越温度扰动时的速度发生了变化, 产生了震荡, 使得胞格形状发生了一定的变化. 从图17(b)和图17(c)中可以看出, 温度间断使得波阵面更加的离散化, 除了存在很长的自然横波之外, 还存在局部上温度间断造成的短的三波结构, 使得波阵面表现出更多的湍流特性. 在温度振幅A = 105 K的例子中, 非常强的温度扰动几乎完全抑制了波阵面的演化, 只能观察到强横波消失的湍流结构的波阵面, 这也可以从压力历史曲线图17(d)中看出, 从波后横向的压力分布中看出横波强度得到了大大地削弱. 上述讨论说明适度的温度扰动能够增加爆轰波自身的不稳定性, 促进爆轰波的传播, 但是过大的温度扰动则会通过抑制横波的发展而限制爆轰波的传播. 除此之外, 由于波阵面更多的湍流特性, 产生较多小的未反应气团.
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-16.jpg'" class="figure_img
figure_type1 bbb " id="Figure16" />
图
16
存在不同振幅人工温度扰动时(kR = 4.0, s = 20ΔI)的压力历史曲线图
Figure
16.
Pressure history for cases with artificial perturbations of kR = 4, s = 20ΔI
下载: 全尺寸图片
幻灯片
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-17.jpg'" class="figure_img
figure_type1 bbb " id="Figure17" />
图
17
人工温度扰动s = 20ΔI, kR = 4.0时,波阵面的纹影结构
Figure
17.
Schlieren photography of detonation front under artificial perturbations with different amplitude with s = 20ΔI and
下载: 全尺寸图片
幻灯片
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/10//lxxb2021-069-15.jpg'" class="figure_img
figure_type1 bbb " id="Figure15" />
图
15
存在不同振幅人工温度扰动时(kR = 4.0, s = 20ΔI)的数值胞格
Figure
15.
Numerical smoked foils for cases with artificial perturbations of kR = 4, s = 20ΔI
下载: 全尺寸图片
幻灯片
3.
结 论
使用反应欧拉方程和两步反应模型对气相爆轰波在非均匀介质中的发展演化过程进行了系统的数值模拟研究, 重点分析了不同波长、幅度的温度扰动对波阵面波系结构的影响, 特别是对自维持传播起到关键作用的横波强度和间距的影响. 为更深入地分析问题, 本研究遵循从简到繁的思路, 先后进行了一维和二维数值研究, 同时计算了激波算例, 进行了对比研究.
一维研究表明, 平面CJ爆轰波在遇到初始温度较高的界面时, 会产生一压力较低但速度较高的透射波向前传播, 同时会产生一个反向传播的反射波. 而当平面CJ爆轰波遇到初始温度较低的界面时, 会产生与热界面完全相反的现象, 即透射波压力增加, 速度降低. 当温度继续降低时, 根据自身不稳定性的区别, 爆轰波透射后会熄爆?重新起爆或完全熄爆. 在本研究中, 二维数值模拟使用平面ZND爆轰波而不是充分发展的胞格爆轰波作为初始条件, 这样做的目的是避免在一开始就引入横波不稳定性, 使得问题过于复杂. 当ZND爆轰波作为入射波传入横向存在单周期温度扰动的非均匀介质中时, 由于在冷区和热区爆轰波传播的速度不同, 温度界面附近会产生剪切层, 导致不同区波阵面的演化存在差异, 冷区甚至存在熄爆和重新起爆的现象, 使得整体上波阵面的结构与在均匀介质中有所不同. 在温度扰动幅值较小时(15 K), 可以清楚地观察到规则的胞格结构, 且温度扰动并没有对胞格形状产生明显的影响. 但当温度扰动幅值增大到35 K时, 胞格的形状发生了一定程度的改变, 主要体现在胞格长宽比的增加. 而当温度扰动幅值增大到105 K时, 数值烟膜中不存在典型的胞格结构, 这说明横波发展受到了温度扰动较强的影响, 波阵面上不断变化的爆轰波三波结构, 只有在中心处能够看到非常弱的横波产生的痕迹.
对于不存在人为温度扰动的均匀介质中, 平面ZND爆轰波可以在传播一定距离后逐步演变为充分发展的胞格爆轰波, 这是由于数值截断误差会诱发波阵面的局部扭曲, 并可在内在不稳定性的控制下发展出横波和波阵面复杂的三波结构. 在ZND爆轰波如此的发展模式下, 若一开始便考虑温度扰动, 在早期阶段, 温度扰动可通过抑制胞格不稳定性的放大来控制横波间隔, 并且维持一段时间的稳态发展. 但随后, 胞格不稳定性再次开始放大, 并最终转变为完全发展的胞格爆轰波. 这表明人为温度扰动的存在通过抑制横波的发展, 延迟了ZND爆轰波向胞格爆轰波的演化, 并且内在不稳定性的增加可以减缓此延迟现象. 另外, 人工的温度扰动可以在一定的范围内抑制胞格不稳定性的发展, 但是并不能够终止这一过程. ZND爆轰波在人工温度扰动下的响应主要受制于两个竞争性的因素: 一是爆轰波内在的不稳定性, 二是人工扰动的波长, 前者是内因, 后者是外因. 但是, 由于上述算例中, 人工温度扰动均设置了相对较小的振幅(15 K), 即便存在温度扰动, 最终也能够得到一个充分发展的胞格爆轰波, 并且这些受扰动胞格爆轰波的主要特征与在均匀介质中的未受扰动胞格爆轰波基本一致. 同时还发现, 如果人工扰动的波长接近均质介质中的胞格爆轰波的横波间隔, 则这两个因素发生同步, 即立即形成具有极其规则胞格图案的胞格爆轰波, 仅在经过很长的距离后, 相位变化才使胞格模式变得稍微不规则.
非均匀介质中温度的不连续性使得波阵面波更加离散, 表现出更多的湍流特性, 除了自然的横波外, 还存在局部弱的三波结构. 在温度扰动幅度较小时, 会增加胞格不稳定性, 从而改变爆轰波阵面的传播机制. 相反, 较大的人工扰动会抑制横波的发展, 此时仅能观察到强横波消失的湍流结构的波阵面, 从而降低爆轰波的不稳定性. 综上, 人工温度扰动的特性和爆轰波内在不稳定性的竞争支配着爆轰前沿的演化过程和传播机制.
