1 引言

镁作为一种金属燃料, 具有熔点低、能量高的特性, 已在固体火箭推进剂和烟火药剂中得到广泛应用^[1]。使用镁基燃料的火箭发动机, 其燃烧性能与镁颗粒着火燃烧特性密切相关。镁颗粒均匀散布在推进剂中, 其燃烧过程十分复杂, 属于非均相固、液、气三相反应, 包括金属颗粒的结团、着火、燃烧及燃烧产物的弥散等系列复杂过程^[2-3]。国内外许多学者都曾对单颗粒镁金属的燃烧特性进行了大量的实验和理论研究, 确定了镁颗粒的燃烧是气相燃烧过程^[3]。Derevyaga^[4]和Kashireninov^[5]对镁颗粒燃烧气相反应区进行了实验测试与分析, 获得了相关动力学参数。Brzustowski等^[6]基于Godsave的碳氢燃料液滴燃烧理论, 建立了简化的镁金属颗粒蒸汽相燃烧模型, 也即所谓快速反应模型, 他假设反应过程受扩散控制, 气相组分按当量比在“薄面区”内快速反应完成, 且反应速率无限快。然而在实际遇到的镁金属微粒燃烧过程中, 由于扩散速率大于气相氧化反应速率, 故气相中不存在“薄反应区”。Edward等^[7-8]选取了直径180~250 mm的镁颗粒, 通过中止燃烧实验发现:未燃的镁颗粒表面含有大量的氧原子, 从而证明反应区不是“薄面区”, 化学反应速率不是无限快。因此, 快速反应模型已经不能精确描述镁金属微粒的实际反应特性。基于以上背景, 本研究建立了镁金属微粒的一维球对称准稳态非快速反应燃烧模型(下文简称为非快速反应燃烧模型), 针对单颗粒镁粒子在大气环境中的燃烧过程进行数值模拟, 分析其燃烧场特性, 并与现有实验结果(文献值)进行对比。研究结果对深入认识镁金属微粒的燃烧机理有一定的参考价值。

2 物理模型

镁金属颗粒(半径r₀)在大气环境中点火后, 颗粒呈现熔融液态, 表面温度可达到沸点(1361 K)^[8]。由于点火时间很短, 对总体燃烧时间影响很小, 故忽略点火过程。镁颗粒在高温下蒸发出的镁蒸汽, 与空气中氧气相互扩散, 在气相中一定区域发生燃烧反应, 此区域即为“反应区(reaction zone)”。对于扩散速率与化学反应速率的快慢关系, 分别有图 1和图 2所示两种物理模型。其中图 1为本研究所建立的非快速反应模型, 反应区为一有限厚度薄壳, 由于反应速率的有限性, 在颗粒表面与反应区之间, 仍然存在扩散进来的氧气分子。图 2为传统的快速反应模型^[6], 反应区仅为薄火焰面, 镁蒸汽与氧气在这个燃烧面全部瞬间反应完全。

针对本研究所建立的非快速反应模型, 提出如下基本假设:

(1) 假设镁颗粒的燃烧过程是准稳态和球对称的, 不考虑颗粒表面的内移效应。且颗粒与大气环境无相对速度, 只有斯蒂芬流引起的球对称径向一维流动;

(2) 在镁颗粒燃烧过程中, 由于热辐射的影响率很小^[6], 故忽略热辐射;

(3) 镁蒸汽与氧气相互扩散, 在一定区域发生气相化学反应, 反应速率遵循Arrhenius定律;

(4) 各气相组分导热系数、比定压热容均为常数。

3 数学模型及计算方法

根据上述非快速反应物理模型(Model 1), 建立镁颗粒球对称一维准定常带化学反应的燃烧场基本方程组^[9]。

(1) 连续方程

$ 4\pi r_0^2{\rho _0}{\upsilon _0} = 4\pi {r^2}\rho \upsilon = 常数 = G $

(1)

式中, G为总质量流量, g·s^-1; r₀为镁颗粒半径, μm; ρ₀为镁颗粒表面处气流密度, g·m^-3; υ₀为镁颗粒表面处气流速度, m·s^-1; ρ为任意r处气流密度, g·m^-3; υ为任意r处气流速度, m·s^-1。

(2) 扩散方程

$ \rho \upsilon = \frac{{{\rm{d}}{f_{{\rm{ox}}}}}}{{{\rm{d}}r}} = \frac{1}{{{r^2}}}\frac{{\rm{d}}}{{{\rm{d}}r}}\left( {{r^2}D\rho \frac{{{\rm{d}}{f_{{\rm{ox}}}}}}{{{\rm{d}}r}}} \right) - {\omega _{{\rm{ox}}}} $

(2)

$ \rho \upsilon \frac{{{\rm{d}}{f_{\rm{F}}}}}{{{\rm{d}}r}} = \frac{1}{{{r^2}}}\frac{{\rm{d}}}{{{\rm{d}}r}}\left( {{r^2}D\rho \frac{{{\rm{d}}{f_{\rm{F}}}}}{{{\rm{d}}r}}} \right) - {\omega _{\rm{F}}} $

(3)

式中, f_ox和f_F分别代表氧气和镁蒸汽的组分质量分数, %; ω_ox表示氧气反应速率, ω_F代表镁蒸汽反应速率; 参数下标ox和F分别表示氧气组分和镁蒸汽组分; D为扩散系数, cm²·s^-1。

(3) 能量守恒方程

$ \rho \upsilon {c_p}\frac{{{\rm{d}}T}}{{{\rm{d}}r}} = \frac{1}{{{r^2}}}\frac{{\rm{d}}}{{{\rm{d}}r}}\left( {{r^2}\lambda \frac{{{\rm{d}}T}}{{{\rm{d}}r}}} \right) + {\omega _{{\rm{ox}}}}{Q_{{\rm{ox}}}} $

(4)

反应速率项ω_ox和ω_F由Arrhenius定律表示:

$ {\omega _{{\rm{ox}}}} = - k\exp \left( { - \frac{E}{{RT}}} \right){f_{{\rm{ox}}}}{f_{\rm{F}}}\frac{{{\rho ^2}}}{{{W_{{\rm{ox}}}}{W_{\rm{F}}}}} $

(5)

$ {\omega _{\rm{F}}} = \frac{{{\omega _{{\rm{ox}}}}}}{\beta } $

(6)

式中, β为化学当量比, k为反应速率常数, m³·(mol·s)^-1; E为活化能, J·mol^-1; W_ox和W_F分别为相应组分分子量; Q_ox为氧气组分的当量反应热, J·g^-1; c_p为比热容, J·(g·K)^-1; λ为导热系数, W·(m·K)^-1。

将(1)式代入(2)、(3)、(4)式, 可得如下微分方程组:

$ \frac{{{\rm{d}}\left( {G{f_{{\rm{ox}}}}} \right)}}{{{\rm{d}}r}} = \frac{{\rm{d}}}{{{\rm{d}}r}}\left( {4{\rm{ \mathsf{ π} }}{r^2}D\rho \frac{{{\rm{d}}{f_{{\rm{ox}}}}}}{{{\rm{d}}r}}} \right) - 4{\rm{ \mathsf{ π} }}{r^2}{\omega _{{\rm{ox}}}} $

(7)

$ \frac{{{\rm{d}}\left( {G{f_{\rm{F}}}} \right)}}{{{\rm{d}}r}} = \frac{{\rm{d}}}{{{\rm{d}}r}}\left( {4{\rm{ \mathsf{ π} }}{r^2}D\rho \frac{{{\rm{d}}{f_{\rm{F}}}}}{{{\rm{d}}r}}} \right) - 4{\rm{ \mathsf{ π} }}{r^2}{\omega _{\rm{F}}} $

(8)

$ \frac{{{\rm{d}}\left( {G{c_p}T} \right)}}{{{\rm{d}}r}} = \frac{{\rm{d}}}{{{\rm{d}}r}}\left( {4{\rm{ \mathsf{ π} }}{r^2}\lambda \frac{{{\rm{d}}T}}{{{\rm{d}}r}}} \right) + 4{\rm{ \mathsf{ π} }}{r^2}{\omega _{{\rm{ox}}}}{Q_{{\rm{ox}}}} $

(9)

对应边界条件为:

$ \left\{ \begin{array}{l} r = {r_0}\;处,T = {T_{\rm{b}}},{f_{{\rm{ox}}}} = {f_{{\rm{ox,0}}}},{f_{\rm{F}}} = {f_{{\rm{F,0}}}}\;\;\;\;\;\;\;\left( {10} \right)\\ r = {r_\infty }\;处,T = {T_\infty },{f_{{\rm{ox}}}} = {f_{{\rm{ox,}}\infty }},{f_{{\rm{F,}}\infty }} = 0\;\;\;\;\;\;\;\left( {11} \right) \end{array} \right. $

式中, T_b和T_∞分别代表镁沸点温度和无穷远处大气环境温度, K; 参数下标b表示镁金属沸点, ∞表示无穷远处大气环境; f_{ox, 0}和f_{F, 0}分别为表面处氧气组分质量分数和镁蒸汽组分质量分数, 同样, f_{ox, ∞}和f_{F, ∞}为氧气组分和镁蒸汽组分在无穷远处大气环境中的质量分数。

将(7)、(8)、(9)式由r₀到r上积分一次, 可得如下关系式:

$ G\left( {{f_{{\rm{ox}}}} - 1} \right) = 4{\rm{ \mathsf{ π} }}{r^2}D\rho \left( {\frac{{{\rm{d}}{f_{{\rm{ox}}}}}}{{{\rm{d}}r}}} \right) - \int_{{r_0}}^r {4{\rm{ \mathsf{ π} }}{r^2}{\omega _{{\rm{ox}}}}{\rm{d}}r} $

(12)

$ G\left( {{f_{\rm{F}}} - 1} \right) = 4{\rm{ \mathsf{ π} }}{r^2}D\rho \left( {\frac{{{\rm{d}}{f_{\rm{F}}}}}{{{\rm{d}}r}}} \right) - \int_{{r_0}}^r {4{\rm{ \mathsf{ π} }}{r^2}{\omega _{\rm{F}}}{\rm{d}}r} $

(13)

$ G\left[ {{c_p}\left( {T - {T_0} + {q_{\rm{e}}}} \right)} \right] = 4{\rm{ \mathsf{ π} }}{r^2}\lambda \frac{{{\rm{d}}T}}{{{\rm{d}}r}} + \int_{{r_0}}^r {4{\rm{ \mathsf{ π} }}{r^2}{\omega _{{\rm{ox}}}}{Q_{{\rm{ox}}}}{\rm{d}}r} $

(14)

联立(12)和(14)式消去反应源项, 可得:

$ G\left[ {{c_p}\left( {T - {T_0}} \right) + {q_{\rm{e}}} + {f_{{\rm{ox}}}}{Q_{{\rm{ox}}}}} \right] = 4{\rm{ \mathsf{ π} }}{r^2}\frac{\lambda }{{{c_p}}}\frac{{\rm{d}}}{{{\rm{d}}r}}\left( {{c_p}T + {f_{{\rm{ox}}}}{Q_{{\rm{ox}}}}} \right) $

(15)

对(15)式在边界上积分, 即可得到蒸发速率(总质量流量)G:

$ G = 4{\rm{ \mathsf{ π} }}{r_0}\frac{\lambda }{{{c_p}}}\ln \left[ {1 + \frac{{{c_p}\left( {{T_\infty } - {T_0}} \right) + \left( {{f_{{\rm{ox,}}\infty }} - {f_{{\rm{ox,}}0}}} \right){Q_{{\rm{ox}}}}}}{{{q_{\rm{e}}} + {f_{{\rm{ox,}}\infty }}{Q_{{\rm{ox}}}}}}} \right] $

(16)

式中, q_e为镁金属的蒸发潜热, J·mol^-1。根据质量守恒定律, 镁颗粒质量m_d的减少速率等于其蒸发速率G, 于是可得到镁颗粒半径r(或直径d)随时间的变化规律:

$ \frac{{{\rm{d}}\left( {{m_{\rm{d}}}} \right)}}{{{\rm{d}}t}} = \frac{{{\rm{d}}\left( {{\rho _{\rm{l}}}{\rm{ \mathsf{ π} }}{d^3}/6} \right)}}{{{\rm{d}}t}} = - G $

(17)

由(17)式可变为:

$ \frac{{{\rm{d}}\left( {{d^2}} \right)}}{{{\rm{d}}t}} = - \frac{{2G}}{{{\rho _{\rm{l}}}{\rm{ \mathsf{ π} }}r}} $

(18)

上式通过积分可计算出镁颗粒完全燃烧反应时间τ_f, 即:

$ {\tau _f} = d_0^2{\rho _{\rm{l}}}{\rm{ \mathsf{ π} }}{r_0}/\left( {2G} \right) = d_0^2/{K_f} $

(19)

式中, K_f为燃烧速率常数, K_f=2G/(ρ_lπr₀)；d₀为颗粒直径, μm；ρ_l为液态镁金属密度, kg·m^-3。

联立(12)和(13)式, 并利用反应速率之间的当量关系, 可消去反应源项, 由r₀到r积分一次, 即可得f_F、f_ox之间关系:

$ \frac{G}{{4{\rm{ \mathsf{ π} }}D\rho }}\left( {\frac{1}{r} - \frac{1}{{{r_\infty }}}} \right) = \ln \left( {\frac{{\beta {f_{{\rm{F,}}\infty }} - {f_{{\rm{ox,}}\infty }} - \beta }}{{\beta {f_{\rm{F}}} - {f_{{\rm{ox}}}} - \beta }}} \right) $

(20)

同样方法联立(13)和(14)式, 可得T、f_ox之间关系:

$ \ln \left[ {\frac{{{c_p}\left( {T - {T_0}} \right) + {f_{{\rm{ox}}}}{Q_{{\rm{ox}}}} + {q_{\rm{e}}}}}{{{q_{\rm{e}}} + {f_{{\rm{ox}},0}}{Q_{{\rm{ox}}}}}}} \right] = \frac{{G{c_p}}}{{4{\rm{ \mathsf{ π} }}\lambda }}\left( {\frac{1}{{{r_0}}} - \frac{1}{r}} \right) $

(21)

将(9)式进行变换, 并令Z=r₀/r, 联立(5)式, 可推导出T、r关系的微分方程:

$ \frac{{{{\rm{d}}^2}T}}{{{\rm{d}}{Z^2}}} + \frac{{G{c_p}}}{{4{\rm{ \mathsf{ π} }}\lambda {r_0}}}\frac{{{\rm{d}}T}}{{{\rm{d}}Z}} + \frac{{r_0^2{Q_{{\rm{ox}}}}k{\rho ^2}{f_{\rm{F}}}{f_{{\rm{ox}}}}{{\rm{e}}^{ - \frac{E}{{RT}}}}}}{{\lambda {Z^4}}} = 0 $

(22)

联立(20)、(21)和(22)式, 结合已知边界条件(10)和(11)式, 通过数值差分可求解出镁颗粒燃烧场温度分布和镁蒸汽、氧气组分质量分数分布规律。

4 数值结果模拟及分析 4.1 计算参数

由于微小单颗粒的粒径范围没有准确界定, 而相关文献实验所取镁颗粒粒径多在120 μm左右, 故本研究也选择直径120 μm的镁颗粒, 对其在大气环境中的燃烧过程进行模拟。大气环境条件为: T_∞=300 K, f_{ox, ∞}=23.2%。非快速反应模型(Model 1)计算中所用的热力学、动力学相关参数均来自文献[4]、^[5]、^[10]和^[11], 见表 1。为验证和对比模型合理性, 又采用传统的快速反应模型^[6](Model 2)进行同样条件下的计算, 所用参数取自文献[12]。

4.2 数值计算结果分析

(1) 燃烧场温度分布

两种燃烧模型算出的镁颗粒燃烧场温度分布如图 3所示。由图 3可见, 两种模型燃烧场的温度分布规律基本一致。流场中的温度随着距离镁颗粒中心O点径向距离的增大, 先升高后逐渐降低。Model 1温度曲线顶部较平滑, 结合图 1物理模型假设, 说明反应区有一定厚度, 其反应区最高温度为3568 K, 对应位置为距离镁颗粒中心O点2.7r₀处。Model 2温度曲线顶部为一尖峰, 结合图 2物理模型假设, 说明反应区为一几何面, 其反应区最高温度为3577 K, 对应位置为距离镁颗粒中心O点3.7r₀处。根据镁金属单颗粒燃烧实验^{[8, 13]}, 测出的燃烧区最高温度约为3373 K, 对应位置约为距离镁颗粒中心O点3r₀处, 图 3结果表明非快速反应模型反应区最高温度所在位置比快速反应模型更接近实测值。两种反应模型火焰最高温度计算值略高于实测值, 主要原因是燃烧实验中存在热损失, 而在计算中未考虑, 所以火焰温度实测值必然低于理论值。

(2) 反应物组分分布

两种燃烧模型算出的镁颗粒燃烧场反应物组分分布如图 4所示。对比图 4a与图 4b可见, 两种模型计算结果差别较大, Model 2中镁蒸汽和氧气在某点相遇后, 瞬时反应完全, 为快速反应过程, 反应组分间不存在共存现象。而Model 1中氧气组分扩散到了镁颗粒表面, 原因在于:气体的扩散速率大于反应速率, 而非快速反应模型正是考虑到镁蒸汽与氧气反应速率的有限性, 故不存在完全反应区, 有部分氧气、镁蒸汽组分会扩散到反应区外部, 这个结果和Edward等人^[7-8]所观察到的镁金属微粒燃烧中止实验后氧原子存在于未燃镁金属中的现象相吻合。

(3) 颗粒燃烧反应时间

针对120 μm的镁金属颗粒, 基于上述公式(19), 可算出镁颗粒燃烧反应时间τ_f。两种燃烧模型的计算结果与实测结果^[14-16]的对比见表 2。由表 2可见, Model 1的计算值与实验值吻合较好, 平均误差约为: 8.75%, 而Model 2的计算值与实验值相差很大。原因在于: Model 2假设化学反应速率无限快, 间接促使镁蒸发速率加快, 导致了镁颗粒燃烧反应时间大幅缩短。因此在计算镁颗粒燃烧反应时间时, 非快速反应模型更为精确合理。

5 结论

建立了镁金属微粒燃烧过程的非快速燃烧反应模型, 并进行了数值模拟。根据模拟结果可得出如下结论:

(1) 非快速燃烧反应模型计算出的镁颗粒燃烧场特征参数分布规律比传统的快速燃烧反应模型更接近实际工况。

(2) 针对120 μm的镁颗粒, 非快速燃烧反应模型算出的反应区最高温度为3568 K, 所在位置距离镁颗粒中心O点2.7r₀处, 颗粒完全燃烧时间为21.9 ms。快速燃烧反应模型算出的反应区最高温度为3577 K, 所在位置距离镁颗粒中心O点3.7r₀处, 颗粒完全燃烧时间为8.64 ms。与实测结果相比, 非快速燃烧反应模型计算结果精度更高。

(3) 非快速燃烧反应模型得到的氧气组分空间分布特性, 与镁颗粒燃烧中止实验后观察到的氧原子存在于未燃镁金属中的现象定性吻合。