摘要:镁颗粒因其能量密度高、点火特性和燃烧效率好的优势, 作为燃料或添加剂应用于爆震燃烧动力系统具有广阔的应用前景. 本文建立了镁颗粒-空气混合物的一维稳态爆震波模型, 数值模拟爆震波稳态传播过程及其内部流场分布. 结果表明, 镁颗粒-空气混合物爆震波仅能以特征值速度稳定自维持传播, 特征值爆震速度的高低并不仅仅取决于反应放热多少, 两相间的相互作用也会影响热能向气相动能的转化效率. 当爆震波末端氧化镁处于熔化过程时, 满足一定的来流速度和镁颗粒密度条件, 爆震波仍能够稳定自维持传播. 气相吸收反应放热膨胀加速至声速的过程主要发生在镁颗粒纯蒸发反应阶段, 但在氧化镁熔化阶段由于熔化过程吸热量大, 使气相吸热膨胀过程近乎停止. 颗粒粒径变化主要影响爆震波尺寸, 而对特征值爆震速度以及波后声速面参数影响甚微. 在常温常压的初始条件下, 爆震波稳定自维持传播过程中波内不涉及氧化镁的汽化离解过程.
关键词: 异相爆震/
Zeldovich-von Neumann-D?ring模型/
镁颗粒燃料

English Abstract

全文HTML

--> --> -->

爆震燃烧具有能量释放速率快、热力循环效率高等优点, 在推进领域中一直都备受关注^[1]. 目前爆震发动机采用的燃料以气体和液体为主. 固体颗粒燃料(镁、铝、碳和硼等)因能量高、易存储、价格低廉, 不仅在常规的固体推进剂领域被广泛采用, 还可作为添加剂改善爆震波胞格质量^[2], 提高气相爆震波传播速度^[3], 提升脉冲爆震发动机性能^[4], 此外也有粉末燃料应用于连续旋转爆震发动机的探索研究^[5-9]. 虽然镁的能量密度低于铝和硼, 但镁较低的熔点和沸点使其点火特性和燃烧效率更优, 其燃烧过程以蒸发气相反应为主, 反应速度较铝、硼更快, 使得应用于爆震领域更有前景. 工业生产中, 与铝粉相比, 镁因其更高的活性更容易发生爆炸, 研究镁的爆震现象对工业生产安全也具有重要意义.
20世纪40年代, 由Zeldovich, von Neumann和D?ring提出的ZND模型^[10]是计算纯气相混合物一维稳态爆震波结构及稳定传播速度的经典模型. Fedorov等^[11]在此模型基础上, 建立了铝颗粒-氧气混合物一维稳态爆震模型, 分析了表征气固两相间的动量、速度和温度延迟的3个松弛参数${\tau _\xi }$, ${\tau _{u}}$和${\tau _{\rm T}}$对流场参数分布的影响, 当满足${\tau _u} < \max \left\{ {{\tau _\xi }, \sqrt {{\tau _{\rm T}}{\tau _\xi }}, 0.14{\tau _{\rm T}}} \right\}$时, 在前导激波后方会出现一个气相和固相颗粒同步的高密度区域; 当满足${\tau _{\rm T}} > {\tau _u}$时, 在颗粒受热点火过程中, 气相温度存在一个局部峰值. 洪滔^[12]也提出了一种适用于气固两相混合物爆震的ZND模型, 在求解方法中采用了一个自定义的参数——颗粒相马赫数, 由于文中没有给出具体算例, 因此上述模型是否合理有待进一步确定. Zhang^[13]在ZND模型的基础上, 针对铝粉-空气混合物建立了一维稳态爆震波模型, 结果表明, 两相混合物爆震与纯气相混合物爆震存在四点明显不同: 第一, 纯气相混合物在前导激波处压力即达到峰值, 两相混合物的压力峰值位于波后气相膨胀与拖曳力压缩作用的平衡处; 第二, 由于气体和颗粒惯性不同, 前导激波对于气相的压缩作用使得波后空燃比明显低于波前; 第三, 颗粒的熔化和蒸发等相变过程使两相混合物的爆震波结构呈现出不同特点(只是提出观点, 并无相关建模或算例结果); 第四, 在适当的速度松弛时间条件下, 在前导激波后的下游颗粒相存在一个高密度区域, 这与Federov等^[11]的结论一致.
由于粉末燃料或反应产物在熔化、蒸发等相变过程中会吸收大量热量, 一方面影响气相吸热膨胀做功, 另一方面会限制反应产物温度, 使其不至过高, 因此有必要针对相变过程对爆震波流场结构及传播过程的影响进行研究. 上述模型以及其他粉末燃料爆震相关的后续研究^[14-26], 大部分未考虑相变过程的影响, 另有少部分为了方便计算确保收敛, 将相变过程进行了简化处理, 如洪滔和秦承森^[17]直接把相变潜热折算在燃料热值中, 并通过设置逆向反应以确保气相温度不会超过反应产物的离解温度; Benkiewicz和Hayashi^[16]通过在高温条件下加入吸热反应的方式来限制反应产物温度等, 这与实际物理化学过程存在明显差距. 文献[13]与文献[27]考虑了反应物颗粒的相变过程影响, 但未考虑燃烧产物的相变影响. 此外, 关于气固两相爆震的研究, 采用的粉末燃料以铝和碳为主, 关于镁粉燃料的研究较少. 鉴于镁颗粒应用于爆震燃烧的优势, 本文建立了镁颗粒-空气混合物一维稳态模型, 分析研究来流速度、相变过程、颗粒相初始密度和半径对爆震波结构及稳定自维持传播过程的影响, 为镁燃料应用于爆震燃烧动力系统提供理论参考.

坐标系固定在前导激波上的一维稳态镁颗粒-空气两相混合物爆震波结构如图1所示. 经过前导激波, 反应物中的气相经绝热压缩过程温度升高, 在诱导区内通过两相间传热过程颗粒相温度升高达到着火温度. 在反应区内化学能快速释放, 气相膨胀加速远离前导激波, 膨胀产生的反向推力用于维持爆震波的传播. 当气相相对于前导激波加速至声速后, 下游流场对上游不再产生影响, 声速面后的气相膨胀过程对维持爆震波传播不再有贡献, 此时声速面又称为Chapman-Jouguet (CJ)面.
图 1 一维稳态爆震波结构示意图
Figure1. Schematic of the steady one-dimensionional detonation wave.

对于一维稳态镁颗粒-空气两相混合物, 其质量、动量和能量守恒方程如下:

$\frac{{\rm{d}}}{{{\rm{d}}x}}\left( {{\rho _1}{u_1} + {\rho _2}{u_2}} \right) = 0, $

$\frac{{\rm{d}}}{{{\rm{d}}x}}\left( {{\rho _1}u_1^2 + p + {\rho _2}u_2^2} \right) = 0, $

$\frac{{\rm{d}}}{{{\rm{d}}x}}\!\left[ {{\rho _1}{u_1}\left(\!{{e_1}\! +\! \frac{{u_1^2}}{2} \!+\! \frac{p}{{{\rho _1}}}}\!\right) \!+\! {\rho _2}{u_2}\left(\!{{e_2} \!+ \!\frac{{u_2^2}}{2} \!+\! {q_{\rm{c}} }}\!\right)}\!\right] \!=\! 0,$

其中下标1和2分别表示气相与颗粒相, ${q_{\rm{c}} }$是镁的热值. 为便于计算求解, 做出以下简化假设:
1)颗粒为球形, 初始粒径大小相同, 且均匀弥撒分布于气体内, 忽略颗粒间相互作用, 颗粒相压强为0;
2)颗粒内温度均匀分布, 颗粒达到熔点温度时若满足气相温度高于颗粒熔点的条件, 颗粒开始发生汽化和燃烧反应;
3)化学反应放能与压力无关, 反应放能仅被气相吸收, 镁只与空气中的氧气反应, 反应方程式为$2{\rm{Mg}} + {{\rm{O}}_2} \to 2{\rm{MgO}}$;
4)反应产物氧化镁算作气相组分^[28,29], 按准气相处理, 但仍考虑其相变过程, 气相热力参数计算方法参考文献[30]. 氧化镁的离解温度当作沸点处理, 物质沸点由Clausius-Clapeyron方程确定^[31], 仅气态工质对气相压强有贡献, 物质熔点不变;
5)相变潜热包含在内能之中, 在氧化镁熔化、离解过程中, 气相温度分别维持在氧化镁熔点和离解温度不变.
根据简化假设5)可知, 温度与内能之间的关系可表示为

${T_1} = \left\{ {\begin{aligned}& {{{{e_1}} / {{c_{\rm{v},1}},}}} \\& {{T_{1,\rm f}},} \\& {{{\!\left( {{e_1} \!-\! {e_{1,\rm f\left( \rm{l} \right)}}} \right)} / {{c_{\rm{v},1}} \!+\! {T_{1,\rm f}},}}} \\& {{T_{1,\rm b}},} \\& {{{\!\left( {{e_1} \!- \!{e_{1,\rm b\left( \rm{g} \right)}}} \right)} / {{c_{\rm{v},1}} \!+\! {T_{1,\rm b}},}}} \end{aligned}} \right.\;\;\begin{aligned} &{{e_1} < {e_{1,\rm f\left( s \right)}}}, \\& {{e_{1,\rm f\left( s \right)}} \leqslant {e_1} \leqslant {e_{1,\rm f\left( \rm{l} \right)}}}, \\ &{{e_{1,\rm f\left( \rm{l} \right)}} < {e_1} < {e_{1,\rm b\left( \rm{l} \right)}}}, \\& {{e_{1,\rm b\left( \rm{l} \right)}} \leqslant {e_1} \leqslant {e_{1,\rm b\left( \rm{g} \right)}}}, \\& {{e_1} > {e_{1,\rm b\left( \rm{g} \right)}}},\end{aligned} $

${T_2} = \left\{ \begin{aligned}& e_2/c_{\rm{v},2(s)}, \\& T_{2,\rm f}, \\& (e_2\!-\! e_{2,\rm f( \rm{l})})/c_{\rm{v},2(\rm{l})}\!+\!T_{2,\rm f}, \\& T_{2,\rm b},\end{aligned}\;\;\begin{aligned}& {e_2} < e_{2,\rm f(s)}, \\& e_{2,\rm f(s)} \!\leqslant\! {e_2} \!\leqslant\! e_{2,\rm f(\rm{l})}, \\& e_{2,\rm f(\rm{l})} \!<\! e_2 \!<\! e_{2,\rm b(\rm{l})}, \\& {e_2} \geqslant e_{2,\rm b(\rm{l})},\end{aligned} \right.$

其中下标f和b分别表示熔点和沸点, (s), (l)和(g)分别表示气态、液态和固态; ${T_{1, \rm f}}$和${T_{1, \rm b}}$分别表示氧化镁的熔点和沸点, ${e_{1, \rm f\left( s \right)}}$和${e_{1, \rm f\left( \rm{l} \right)}}$分别为温度达到氧化镁熔点时氧化镁全部为固态和液态所对应的气相比内能, ${e_{1, \rm b\left( \rm{l} \right)}}$和${e_{1, \rm b\left( \rm{g} \right)}}$分别为温度达到氧化镁沸点时氧化镁全部为液态和气态所对应的气相比内能. 颗粒相参数同理. 这里定压比热${c_{\rm{v}, 1}}$和${c_{\rm{v}, 2\left( s \right)}}$是一个等效平均值, 随温度升高而变化, 其与温度的值是在确定内能后通过多次迭代求解获得的.
将守恒方程(1)—(3)中气相与颗粒相分离, 得到关于气相、颗粒相以及气相中各组分的守恒方程如下:

$\frac{{\rm{d}}}{{{\rm{d}}x}}\left( {{\rho _1}{u_1}} \right) = {S_1}, $

$\frac{{\rm{d}}}{{{\rm{d}}x}}\left( {{\rho _1}u_1^2 + p} \right) = - {F_{\rm{d}}} + {S_1}{u_2}, $

$\begin{split} &\frac{{\rm{d}}}{{{\rm{d}}x}}\left[ {{\rho _1}{u_1}\left( {{e_1} + \frac{{u_1^2}}{2} + \frac{p}{{{\rho _1}}}} \right)} \right] \\ =\, & - {Q_{\rm{d}}} - {F_{\rm{d}}}{u_2} + {S_1}\left( {{e_2} + \frac{{u_2^2}}{2} + {q_{\rm{c}}}} \right),\end{split}$

$\frac{{\rm{d}}}{{{\rm{d}}x}}\left( {{\rho _2}{u_2}} \right) = - {S_1}, $

$\frac{{\rm{d}}}{{{\rm{d}}x}}\left( {{\rho _2}u_2^2} \right) = {F_{\rm{d}}} - {S_1}{u_2}, $

$\frac{{\rm{d}}}{{{\rm{d}}x}}\left[ {{\rho _2}{u_2}\left( {{e_2} + \frac{{u_2^2}}{2}} \right)} \right] = {Q_{\rm{d}}} + {F_{\rm{d}}}{u_2} - {S_1}\left( {{e_2} + \frac{{u_2^2}}{2}} \right), $

$\frac{{\rm{d}}}{{{\rm{d}}x}}\left( {{\rho _{{\rm{O} _2}}}{u_1}} \right) = {S_{{\rm{O} _2}}}, $

$\frac{{\rm{d}}}{{{\rm{d}}x}}\left( {{\rho _{\rm{MgO} }}{u_1}} \right) = {S_{\rm{MgO} }}, $

$\frac{{\rm{d}}}{{{\rm{d}}x}}\left( {{\rho _{{\rm{N} _2}}}{u_1}} \right) = 0, $

其中${S_1}$为单位体积内颗粒相对气相的质量转化速率. 根据文献[32], 镁颗粒燃烧初期异相反应影响明显, 采用经验公式^[13,33]

${S_1} = \frac{{3{\rho _2}}}{{{t_{\rm{b}}}}}\left( {1 + 0.276R{e^{0.5}}P{r^{0.33}}} \right), $

${t_{\rm{b}}} = \frac{{{{\left( {2{r_0}} \right)}^m}}}{{k{\xi _{{\rm{O} _2}}}\left( {1 + 0.217\sqrt {R{e_0}} } \right)}}, $

式中, $Pr$为Prandtl数; $Re = {{2{\rho _1}r\left| {{u_1} - {u_2}} \right|} / {{\mu _1}}}$为两相间的Renolds数; ${t_{\rm{b}} }$表示颗粒燃烧时间; $r$表示颗粒半径; ${\mu _1}$为气相黏性系数; 下标0表示初始态; ${\xi _{{\rm{O} _2}}}$表示气相中氧气质量分数; 参数$R{e_0}$, $m$和$k$的表达式为

$R{e_0} = {{2{\rho _1}\left| {{u_1} - {u_2}} \right|{r_0}} / {{\mu _1}}}, $

$m = \left\{ \begin{aligned}& 10{p^{ - 0.14}},\;\;\;\;\;\;\;{10^5}\;{\rm Pa} \leqslant p \leqslant 22 \times {10^5}\;{\rm Pa}, \\& 1.25,\;\;\;\;\;\;\;\;\;\;\;\;\;\; p > 22 \times {10^5}\;{\rm Pa},\end{aligned} \right.$

$k = \left\{ \begin{aligned} & 3.079 \times {10^{ - 10}}{p^{0.83}},\;\;\;\;\;{10^5}\;{\rm Pa} \!\leqslant\! p \!\leqslant\! 22\! \times \!{10^5}\;{\rm Pa} \\ & 5.658 \times {10^{ - 5}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;p > 22 \times {10^5}\;{\rm Pa},\end{aligned} \right.$

其中$m \leqslant 1.25 < 2$, 表明镁颗粒燃烧过程在一定程度上受表面反应动力学控制. 当镁颗粒温度达到沸点时, 镁颗粒燃烧过程与液滴纯蒸发燃烧类似, 由于气相反应特征时间远小于液滴纯蒸发燃烧时间, 液滴扩散燃烧过程反应速率可用纯蒸发速率代替. 根据文献[34]有

${S_1} = \frac{{12nr{\lambda _1}Nu\left( {{T_1} - {T_2}} \right)}}{{\Delta {h_{\rm{Mg},b}}}}, $

其中Nusselt数$Nu = 2 + 0.6 R{e^{0.5}}P{r^{0.33}}$, $\Delta {h_{\rm{Mg}, b}}$为单位质量镁的汽化潜热, ${\lambda _1}$为气相导热系数. 单位体积内, 气相对颗粒相的作用力为^[17]

${F_{\rm{d}}} = \frac{n}{2}{\text{π}}{r^2}{C_D}{\rho _1}\left| {{u_1} - {u_2}} \right|\left( {{u_1} - {u_2}} \right), $

式中,

${C_{\rm{D}}} = \left\{\begin{aligned}& \frac{{24}}{{Re\left( {1 + {{R{e^{{2 / 3}}}} / 6}} \right)}},\quad {Re < 1000},\\ & {0.424}, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{Re \geqslant 1000},\end{aligned} \right. $

$n$为单位体积内颗粒数目, 整个流动燃烧过程中颗粒总数守恒, 因此有

$n = {{{n_0}{u_{2,0}}} / {{u_2}}}, $

${n_0} = \frac{{3{\rho _{2,0}}}}{{4{\text{π}}r_0^3{\rho _{\rm{r} 2}}}}, $

其中${\rho _{\rm{r} 2}}$为颗粒的物质密度. 单位体积内, 气相对颗粒相的对流换热速率为

${Q_{\rm{d}}} = 2n{\text{π}}r{\lambda _1}Nu\left( {{T_1} - {T_2}} \right).$

根据文献[13], 辐射换热在两相混合物爆震过程中的影响可忽略不计. 为表征反应进行程度, 引入参数反应完成度$\eta $, 两相来流为贫氧工况时, 其定义为

$\eta = \frac{{{\varphi _{{\rm{O} _2},0}} - {\varphi _{{\rm{O} _2}}}}}{{{\varphi _{{\rm{O} _2},0}}}},$

富氧工况时, 其定义为

$\eta = \frac{{{\varphi _{2,0}} - {\varphi _2}}}{{{\varphi _{2,0}}}},$

其中$\varphi $为物质或相在两相中所占质量分数.

由于CJ面下游流场扰动对上游流场无影响, 默认计算求解域定为前导激波面至CJ平面处. 边界条件设定如下: 对于x < 0区域, 人为给定来流初始参数${\rho _{1, 0}}$, ${\rho _{2, 0}},u_0,T_0$和r₀. 在x = 0处, 对于气相, 根据Rankin-Hugoniot正激波关系式, 可得到前导激波面后的von Neumann状态; 对于颗粒相, 由于惯性作用, 认为通过前导激波面其状态参数保持不变. 若无特别说明, 本文的气相来流为常温常压下的空气, 即${\rho _{1, 0}} =$ 1.29 kg/m³和T₀ = 300 K. 常微分方程(6)—(14)在求解域内采用4阶Runge-Kutta法数值求解.
联立方程(6)—(14), 关于气相速度的常微分方程可表示为如下形式:

$\frac{{{\rm{d}}{u_1}}}{{{\rm{d}}x}} = \frac{\phi }{{1 - Ma_1^2}},$

其中$M{a_1}$为气相马赫数, $\phi $是一个与化学反应以及两相间质量、动量和能量转化有关的函数, 表征化学反应放热速率中实际转化为气相动能的部分, 且对于$\eta = 1$, 有$\phi = 0$. 当爆震波结构满足分子分母同时趋近于0的正则条件时, 所对应的初始来流速度为特征值爆震速度${u_{0, \rm{ZND} }}$. 此外, 在氧化镁熔化过程中, 结合假设5)以及理想气体状态方程, 可得到其气相速度的常微分方程为

$\frac{{{\rm{d}}{u_1}}}{{{\rm{d}}x}} = \frac{{{\phi _{\rm{MgO},f}}}}{{1 - \gamma Ma_1^2}}, $

其中$\gamma $为气相比热比. 同理可得氧化镁汽化离解过程中的气相速度常微分方程为

$\frac{{{\rm{d}}{u_1}}}{{{\rm{d}}x}} = \frac{{{\phi _{\rm{MgO},b}}}}{{1 - \dfrac{\gamma }{2}Ma_1^2}}.$

2

-->

3.1.网格无关性验证

选取${\rho _{2, 0}} = $ 0.445 kg/m³、来流速度为相应特征爆震速度u₀ = 1782.28 m/s以及颗粒半径r₀ = 5 μm的初始条件, 分别以1, 0.1和0.01 mm三种尺寸的网格计算得到前导激波后气相密度分布, 结果如图2所示. 通过图2内插图所示的局部放大图可知, 随着网格加密, 密度峰值略有提高, 峰值位置略有后移, 但在总体上三者计算结果无明显差别. 因此本文后续算例若无特别说明网格均为0.1 mm.
图 2 不同网格条件下气相密度分布(内插图为局部放大图)
Figure2. Spatial distribution of the gas-phase density with different grid sizes. Inset shows the view of partial enlargement.

由于镁粉爆震相关的文献较少, 此处借助铝粉爆震相关的实验和计算结果进行对比, 间接说明本文模型计算结果的合理性: 以文献[26]中验证算例的工况(燃料与空气当量比为1)为参考, 在燃料和空气当量比为1的条件下, 考虑到镁和铝自热热值及其理论空燃比, 理论上单位体积内镁反应释放的热量多于铝, 导致镁的爆震波速度应高于铝, CJ面流场参数也应高于铝, 表1列出了相关参数及结果, 表中数据说明本模型计算结果具有一定合理性.

燃料种类	质量热值/MJ·kg–1	理论空燃比	爆震波速度/m·s–1	CJ面密度/kg·m–3	CJ面速度/m·s–1	CJ面压力/MPa
镁	25.06	2.87	1782.28	3.27	942.7	2.87
铝	31.07	3.83	1650	2.43	673	2.04

表1镁和铝的结果对比
Table1.Comparison of results of magnesium and aluminum detonation.

2 -->

3.2.来流速度的影响

为讨论来流速度的影响, 在来流条件为${\rho _{2, 0}} = $ 0.445 kg/m³ (空燃比满足理论当量比)和r₀ = 5 μm的基础上, 以来流速度分别为0.8、1倍和1.2倍特征值爆震速度为例, 计算对应的爆震波结构, 截取末端两参数$1 - Ma_1^2$和$\phi $的空间分布, 如图3所示.
图 3 不同来流速度条件下爆震波末端$1 - Ma_1^2$与$\phi $的空间分布
Figure3. Spatial distribution of $1 - Ma_1^2$ and $\phi $ at the end of detonation wave with different inlet velocity.

由图3可知, 来流速度为特征值爆震速度时, 在$1 - Ma_1^2 \to 0$处同时有$\phi \to 0$, (26)式满足正则条件. 当来流速度高于特征值爆震速度(1.2倍)时, 在$\phi \to 0$处有$1 - Ma_1^2 \gg 0$, 说明在气相达到声速前反应已完成, 前导激波后的下游流场都是亚声速, 下游流场的扰动都能穿过反应区并减弱爆震强度, 这种情况即通常所说的过驱爆震, 在无外力维持的条件下爆震波无法稳定自维持传播. 当来流速度低于特征值爆震速度(0.8)时, 在$1 - Ma_1^2 \to 0$处有$\phi \gg 0$, 在气相加速至声速时出现奇点, 说明此时爆震波也无法稳定自维持传播. 综上所述, 特征值爆震速度是镁颗粒-空气混合物理论上唯一的稳定自维持传播速度.
2 -->

3.3.相变过程的影响

选取${\rho _{2, 0}} = $ 0.445 kg/m³、来流速度为相应特征爆震速度u₀ = 1782.28 m/s以及颗粒半径r₀ = 5 μm的初始条件, 计算得到爆震波流场参数如图4所示.
图 4 爆震波流场参数分布　(a) 密度和质量分数; (b) 速度和压力; (c) 温度
Figure4. Parameters distribution in detonation wave: (a) Density and mass fraction; (b) velocity and pressure; (c) temperature.

由图4可以看出, 气相通过前导激波面达到von Neumann状态, 与颗粒相之间存在较大的速度差, 之后通过拖曳力作用两相速度趋于相等. 在诱导区内通过两相间对流换热、辐射换热, 颗粒相温度上升直至镁熔点, 之后颗粒相进入反应区开始反应. 根据氧化镁质量分数${Y_{{\rm{Mg}}\rm{O} }}$曲线可以看出, 在颗粒相温度达到沸点(x = 0.2043 m处)前后反应速率有明显变化: 颗粒相温度达到沸点前反应速率较低, 在此阶段气相密度、速度、温度以及颗粒相密度、速度变化均无明显变化, 说明在镁熔化阶段时反应放热主要被镁的熔化过程吸收, 在镁完全熔化后反应放热主要用于提高颗粒相温度, 且反应放热速率分别与镁熔化时吸热速率及熔化后的相间传热速率大致相等; 颗粒相温度达到沸点后, 反应速率明显升高, 放热速率增加, 气相温度快速上升, 气相膨胀导致密度、压力降低, 速度增加. 在前导激波下游出现了气相和颗粒相的同步高密度区, 诱导区内气相温度存在局部峰值, 与文献[11]的结果一致.
在氧化镁熔化过程中(对应x = 0.2341 m至x = 0.2389 m之间的区域), 气相密度、速度和压力曲线发生明显变化: 密度升高, 速度和压力基本不变. 这是因为在假设5)的条件下, 气相温度保持不变, 此时反应放热大部分被非气态的氧化镁吸收, 气相中气态工质由吸热升温膨胀转变为等温膨胀, 膨胀速率降低, 同时颗粒燃烧持续消耗氧气, 膨胀速率进一步降低近乎中止, 可近似认为此时反应放热全部用于熔化吸热, 压力和速度基本不变, 气相密度升高主要是反应持续生成氧化镁所致.
在爆震波稳态传播条件下, 各相变区长度与相变过程时间尺度呈对应关系, 各区域长度具体数值如表2所列. 由于镁颗粒在燃烧过程中伴随着镁的蒸发, 燃烧反应区即对应镁颗粒的蒸发区(包括沸腾前的混合反应区和沸腾后的纯蒸发反应区). 由表2中数据可知, 氧化镁和镁的熔化过程持续时间不足整个爆震过程的5%, 爆震反应过程中的时间尺度由小到大依次为: 氧化镁熔化、镁熔化、点火延迟、燃烧反应. 由图4(c)可知, 在已知来流条件下CJ面处气相温度未达到氧化镁沸点, 不涉及氧化镁的汽化离解过程.

	总长度	诱导区	燃烧区 (Mg蒸发)	Mg 熔化	MgO 熔化	MgO 离解
空间长度/m	0.3324	0.0354	0.297	0.0116	0.0048	×
占比	1	0.106	0.894	0.035	0.014	×

表2诱导区及相变区长度
Table2.Length of induction region and phase-transition region.

为进一步探究相变过程对爆震波结构的影响程度, 在相同来流条件下(诱导区长度相同), 计算忽略部分相变后气相加速至0.99倍声速截面反应完成度${\eta _{0.99\rm{CJ} }}$和对应的燃烧区占比, 与全部考虑相变的工况进行对比, 结果如表3所列. 如前文所述, 镁的蒸发贯穿整个颗粒相燃烧过程, 镁蒸发过程所吸收的热量由气相通过相间换热过程提供, 根据模型假设, 蒸发后瞬时完成反应, 且反应放热全部被气相吸收, 镁通过蒸发燃烧由颗粒相转化为气相, 可以等效看作镁蒸发过程中吸收的热量仅是气相内部的能量转化, 对颗粒相未产生直接影响. 同时此工况不涉及氧化镁的汽化离解过程, 因此仅对比镁熔化和氧化镁熔化的影响.

$\Delta h_{\rm{Mg},f}$/MJ·kg–1	$\Delta h_{\rm{MgO},f}$/MJ·kg–1	$\eta_{0.99\rm{CJ}}$	燃烧区占比
0.354	1.946	0.995	0.894
0	1.946	0.952	0.885
0.354	0	0.573	0.872
0	0	0.549	0.865

表3相变过程对爆震波结构影响对比
Table3.Effect of phase transition on structure of detonation wave.

根据表3中数据可知, 忽略镁和氧化镁的熔化潜热后, 0.99倍声速截面处的反应完成度和燃烧区占比均降低, 这是因为在忽略熔化潜热后, 有更多反应放热通过气态工质膨胀做功转化为气相动能, 理论上气相加速至声速所需时间更短, 导致燃烧区缩短, 在0.99倍声速截面处颗粒相反应完成度降低, 此工况对应更高的特征值爆震速度. 此外, 氧化镁熔化的影响远大于镁熔化, 一方面是由于镁熔化过程吸收的热量会随着后续液态镁颗粒汽化后燃烧转化为气相内能, 参与气相膨胀做功, 而氧化镁熔化过程吸收的热量完全不参与气相膨胀做功; 另一方面, 根据氧化镁和镁的单位质量熔化潜热及其熔化过程中对应的密度, 可以得到单位体积内氧化镁熔化所需的能量大于镁熔化过程.
图5所示为气相中的气态工质自前导激波面至CJ面的压力-比体积曲线. 可知, 气态工质在镁和氧化镁的熔化过程中所做的膨胀功占总膨胀功的比例极小. 在诱导区内, 由于两相间的拖曳力作用, 气态工质经历先轻微膨胀后压缩的过程. 在燃烧区内, 气态工质先经历升压膨胀过程, 后转为降压膨胀过程, 转折处对应的颗粒相未达到沸点. 氧化镁完全熔化后的阶段与熔化前相比, 随比体积增加, 压降幅度更大, 这是由于液态氧化镁比热较固态更高, 熔化后反应放热转化为气态工质膨胀功的比例减少. 整体来看, 气态工质在纯蒸发反应阶段所做膨胀功最多, 对应图4(a)和图4(b)中气相膨胀加速过程主要发生在颗粒相达到沸点后. 此外, 在氧化镁熔化过程中(图5(b)点划线部分), 气态工质是伴随着氧气质量亏损的等温膨胀过程, 气态工质整体密度降低, 但氮气质量分数升高, 导致气态工质的平均气体常数升高, 气态工质密度和平均气体常数的综合作用, 使气相压力呈先降低后增加的趋势.
图 5 气态工质爆震燃烧过程压力-比体积曲线　(a) Mg熔化; (b) MgO熔化
Figure5. Pressure-specific volume curve of detonation combustion in gaseous phase:　(a) Mg fusion; (b) MgO fusion.

2 -->

3.4.颗粒相初始密度的影响

针对来流颗粒相初始半径r₀ = 5 μm, 计算得到爆震波特征值爆震速度和CJ面的流场参数随颗粒相初始密度的变化, 结果如图6所示. 在颗粒相初始密度小于0.148 kg/m³时, CJ面气相温度低于氧化镁熔点, 氧化镁为固态. 在颗粒相初始密度大于0.176 kg/m³时, CJ面气相温度高于氧化镁熔点, 氧化镁为液态. 在0.169—0.176 kg/m³范围内, CJ面处的氧化镁处于熔化过程中. 对于0.148—0.168 kg/m³范围内, 无法直接计算得到特征值爆震速度, 需要进一步讨论.
图 6 流场参数随颗粒相初始密度的变化　(a) 特征值速度和激波诱导下限速度; (b) 密度; (c) 速度和压力; (d) 温度
Figure6. Variation of parameters with different initial concentration of particle phase:　(a) Eigenvalue detonation velocity and minimum velocity of particle ignition induced by shock wave; (b) density; (c) velocity and pressure; (d) temperature.

在0.148—0.168 kg/m³范围内, 爆震波末端的氧化镁处于熔化过程中, 对应的来流速度有一个范围, 其下限为气相温度达到氧化镁熔点前气相速度出现奇点对应的最大来流速度, 上限为末端气相温度高于氧化镁熔点对应的最小来流速度. 通过计算可知在上述初始密度以及来流速度范围内, 氧化镁熔化阶段${\phi _{\rm{MgO}, f}} \leqslant 0$始终成立, 同时来流速度自下限增加至上限过程中, 对应的$1 - \gamma Ma_1^2$由负变正单调递增. 根据(27)式, $1 - \gamma Ma_1^2 < 0$是气相能够加速至声速的必要条件, 但在上述初始密度范围内, 气相加速至声速前反应已完成, 与3.2节中的过驱爆震波末端结构类似, 下游流场的扰动穿过反应区减弱前导激波强度, 存在一个使来流速度降低的趋势. 来流速度低于下限时, 爆震波内气相速度出现奇点, 存在一个使来流速度增加的趋势, 因此可认为在0.148—0.168 kg/m³范围内爆震波无法以一特定的速度稳定自维持传播, 取而代之的是在以下限为中心的速度范围内震荡传播. 造成震荡传播的主要原因还是如3.3节分析的那样, 氧化镁熔化过程吸收热量太多, 气相膨胀过程几乎停止, 无法正常加速至声速, 下游流场的扰动穿过反应区减弱前导激波强度, 当前导激波强度降低直至爆震波内温度均低于氧化镁熔点时, 这时无熔化过程吸热, 有充足的能量用于气相膨胀加速, 导致前导激波强度增加, 爆震波即在上述两综合作用下震荡传播. 故在图5中对应的CJ面参数无法确定, 仅在图6(a)中以虚线表示对应的均值速度.
图6(a)给出了激波诱导镁颗粒相点火的来流速度下限, 其随颗粒相初始浓度增加无明显变化, 特征值爆震速度曲线和激波诱导点火下限速度曲线在0.063和1.06 kg/m³处相交, 分别对应爆震波的稳定自维持传播的颗粒相初始密度上限和下限. 从图6整体来看, CJ面处的气相密度、压力和温度随颗粒相初始密度增加大致呈先增后减的变化趋势, 曲线峰值对应的颗粒相初始密度为理论当量比浓度0.445 kg/m³, 这是因为在理论当量比条件下燃料和氧化剂都完全反应, 燃烧反应释放能量最多且产物中没有多余颗粒相吸收反应放热影响气相膨胀做功. 特征值爆震速度和CJ面处气相速度曲线峰值对应的颗粒相初始密度为0.36 kg/m³, 明显小于理论当量比浓度, 这是因为在富氧工况下反应释放的总能量随颗粒相初始密度的增加而增加, 同时因两相间不平衡引起的做功、换热损失(如颗粒数目增加使得气相膨胀加速过程中克服拖曳力做功, 以及被氧化镁熔化吸收而无法参与气相膨胀做功的那部分能量等)也随之增加, 颗粒相初始密度增加导致的气相膨胀做功增益以及相间损失增益在0.36 kg/m³处达到平衡, 继续提高颗粒相初始密度, 相间损失增益大于气相膨胀做功增益, 导致特征值爆震速度和CJ面气相速度降低. 最高爆震速度对应的初始颗粒相浓度小于对应的理论当量比浓度, 这与文献[35]中针对铝粉-空气混合物的理论结算结果趋势一致.
在初始空燃比对应贫氧工况时, 随着颗粒相初始密度增加, CJ面处剩余颗粒相密度逐渐增加. 颗粒相因受气相拖曳力作用, CJ面处两相速度相差不大, 仅随颗粒相初始密度增加而略微增大. 由于颗粒相在CJ面达到沸点, 根据Clausius-Clapeyron方程, CJ面颗粒相温度随初始密度增加的变化趋势与气相压力一致.
2 -->

3.5.颗粒半径的影响

对于球形颗粒, 内部温度均匀分布的假设需满足Biot数条件^[36]

$Bi = \frac{{{h_1}r}}{{{\lambda _2}}} \leqslant 0.033,$

其中${h_1}$为气相对流换热系数. 根据相关特征数定义及经验公式, 可得颗粒半径与Biot数之间的关系式:

$r = \frac{1}{{0.36P{r^{0.67}}}}{\left( {Bi\frac{{{\lambda _2}}}{{{\lambda _1}}} - 2} \right)^2}\frac{{{\mu _1}}}{{2{\rho _1}\left| {{u_1} - {u_2}} \right|}}.$

在满足Biot数条件的前提下, 给定来流速度${u_0}$, 以von Neumann状态为参照, 可估算出适用于内部温度均匀假设的颗粒半径上限. 此处选取颗粒半径1—20 μm作为变量范围, 计算了颗粒相初始密度分别为0.8, 1.0倍和1.2倍的理论当量比密度条件下对应的特征值爆震速度、CJ面气相参数以及诱导区、燃烧区长度随颗粒半径的变化, 结果如图7所示.
图 7 流场参数随颗粒初始半径的变化　(a)速度; (b)密度和压力; (c)温度; (d)爆震波长度
Figure7. Variation of parameters with different initial particle radii:　(a) Velocity; (b) density and pressure; (c) temperature; (d) length of detonation wave.

根据图7(d)可知, 爆震波长度随颗粒半径增加单调递增, 这是由于随着粒径增加, 两相间界面积减小, 导致单位体积内的换热速率、反应速率降低, 相应的颗粒点火和燃烧时间延长. 根据图7(a)—(c), 特征值爆震速度以及CJ面气相参数随颗粒半径增加无明显变化, 说明在爆震波内, 随着粒径增加, 燃烧反应时间变长的同时, 气相吸热膨胀加速至声速的时间也随之变长, 且两者的变化规律近似几乎“同步”, 这才导致不同粒径条件下特征值爆震速度以及CJ面气相参数差别很小.
综合来流速度、颗粒相初始浓度以及颗粒半径的影响, 当颗粒相初始密度为满足当量空燃比(0.445 kg/m³)且爆震波能够稳定自维持传播的条件下, CJ面处气相温度达到理论最大值, 但仍然明显低于对应压力下的氧化镁汽化离解温度, 因此可认为在常温常压的初始条件下, 镁颗粒-空气混合物爆震燃烧过程不涉及氧化镁的汽化离解.

本文针对镁颗粒-空气混合物爆震建立了一维稳态模型, 通过数值模拟不同工况下的爆震波稳态自维持传播过程, 获得了爆震波内流场参数的分布以及相变过程、来流速度、颗粒半径以及颗粒初始密度等因素对爆震波结构的影响规律. 研究表明, 气相吸收反应放热通过膨胀做功逐渐加速至声速的过程主要发生在镁颗粒纯蒸发反应阶段, 镁和氧化镁熔化过程持续时间占整个点火燃烧过程比例较小. 在常温常压的初始条件下, 爆震波稳定自维持传播过程中波内理论最高温度低于产物氧化镁的离解温度, 因此不涉及氧化镁的汽化离解过程. 镁熔化过程吸收的热量伴随化学反应进入气相用于膨胀做功, 因此镁熔化对爆震波结构影响较小; 氧化镁熔化过程吸收的反应放热总量多, 致使气态工质膨胀做功过程近乎停止, 且所吸收的热量无法用于气相膨胀做功, 因此氧化镁熔化过程对爆震波结构影响较大.
镁颗粒-空气混合物爆震波仅能够以特征值速度稳定自维持传播, 低于此速度流场内出现奇点, 高于此速度, 波后无法加速至声速, 下游流场扰动可穿过反应燃烧区, 减弱爆震波强度. 当爆震波末端氧化镁处于熔化过程时, 满足一定的来流速度和镁颗粒密度条件, 爆震波仍能够稳定自维持传播, 否则爆震波理论上仅能以某一速度为均值震荡传播. 特征值爆震速度峰值对应的颗粒相初始密度小于密度、压力、温度三者峰值对应的理论当量比密度, 说明特征值爆震速度的高低并不是仅仅取决于反应放热多少, 两相间的相互作用也会影响热能向气相动能的转化效率. 在基于镁颗粒内部温度均匀分布的前提下, 颗粒粒径变化主要影响爆震波尺寸, 而对特征值爆震速度以及CJ面参数影响甚微.
本模型中涉及反应的部分是基于已有的经验模型, 计算结果由于缺乏相关实验数据难以直接验证, 本文通过与现有文献中铝的计算结果和研究成果进行对比分析, 间接证明了本模型结果的可靠性, 对采用粉末燃料的爆震动力装置设计具有一定的指导意义. 基于本文的工作, 下一步可开展镁颗粒-空气混合物非稳态爆震燃烧数值模拟的相关研究.