基于化学动力的强迫对流下硼粒子的燃烧特性

(1.火箭军研究院，北京 100091；2.国防科技大学 高超声速冲压发动机技术重点实验室，长沙 410073)

方传波1,2，夏智勋2，张旭荣1，余 勇1，袁天保1

0 引言

良好的着火燃烧性能是实现镁、铝、硼等固体粒子在发动机内高效释放能量的基础。当前，纳米粒子及其衍生品以其良好的着火和燃烧性能也开始引起广泛注意[1-2]。对固体火箭冲压发动机而言，受进气道冲压空气的影响，补燃室内的多相湍流掺混过程非常剧烈，气-固/液多相之间的速度、温度差异对粒子相的着火燃烧过程影响很大，必须对该过程展开深入研究。

针对硼粒子的着火燃烧性能，张鹏等[3]分析了温度和气相组分对气相硼燃烧的化学反应动力学参数影响。敖文等[4]利用热天平和激光点火装置，开展了气流速度对晶体硼粒子热氧化及点火燃烧特性的试验研究，但其研究对象为硼粉堆，在硼粉堆的着火燃烧过程中，气相组分的扩散和粉堆的热量累积均受硼粉粒径及粉堆的孔隙率影响，这与冲压发动机补燃室内的单个硼粒子燃烧有一定差别。胡建新等[5]基于薄火焰面假设研究了强迫对流下的硼粒子燃烧速率模型，但该模型属于扩散燃烧，未考虑有限化学反应动力的影响，且该模型不适用于粒子尾部出现分离流的情况。文献[6]基于单步总包反应假设，通过求解硼粒子周围化学反应流控制方程组，研究了强迫对流因素对硼粒子燃烧速率的作用机制，并拟合出了硼粒子的燃烧质量流率通量公式，为工程问题的应用提供了方便。但在强迫对流下，硼粒子的整个燃烧过程可能受到强迫对流的流体动力学和有限化学反应动力的耦合影响，基于单步总包反应的研究尚不能完全反映硼粒子在强迫对流下的燃烧特性。目前，国内外针对强迫对流下硼粒子燃烧特性的研究很少，正如文献[7]所言，仍需对此科学问题开展进一步深入研究。

1 物理数学模型表征方法

1.1 硼粒子燃烧的多步反应机理

Yetter等[8-9]经过多年研究，逐渐建立了相对完善的硼燃烧化学动力学机理，但其涉及的基元反应数目巨大，若计算涵盖所有气相组分，计算成本巨大。因此，本文以氧化性气体仅为氧气组分，并视氮气为惰性气体为例展开分析。气相反应表征为

(1)

则第k个气相组分的反应速率为

(2)

(3)

式中kfi、kri分别为正向、逆向反应速率常数。

硼粒子燃烧的表面反应机理和动力学参数如表1所示。其中，A为反应的指前因子，cm/s，n为级数，E为反应活化能，kcal/mol。

表面反应表征为

Bs/l+Z1g+Z2g+…+Zng→生成物

(4)

其中，Z表示气相组分，后缀s/l、g分别表示固相/液相和气相。气相组分的比热容等热物性参数采用JANAF数据库中的数值[10]。本文暂不考虑硼的融化过程，仅以液态硼的反应展开分析，表面组分的生成焓采用文献[11]中的数据。

若将表面反应也表征为

(5)

则第k个组分的反应速率为

(6)

相应的气相反应动力学反应机理如表2所示，其中，指前因子A单位为cm3/(mol·s)，反应活化能E单位为kcal/mol。

表1 表面反应动力学参数[7]

表2 气相反应动力学参数[7]

1.2 气相燃烧产物的凝结动力学模型

在实际的固体火箭冲压发动机工作过程中，补燃室内部环境很可能促使气相B2O3的凝结[12]，其凝结热值占据硼粒子完全燃烧所释放热值的相当大一部分，需要展开分析。气相B2O3凝结过程可表示为

B2O3g=B2O3l

经典的凝结形核理论经过Chiu等[13]的发展逐渐完善，得到单位体积内的凝相组分反应速率为

(7)

其中,IKin为粒子相通量，cm-3·s-1；r*为最大吉普斯自由能对应的凝相粒子半径，其具体计算方法见文献[13]，硼的氧化物在500～2100 ℃范围内，表面张力系数为[14]

σ= (72.11-33.38×10-3T1+70.57×

(8)

式中 下标l表示液态B2O3。

固体火箭冲压发动机补燃室内的气相温度可达1500～2500 K，若取2000 K作为典型代表温度，相应的气相B2O3饱和蒸气压约为0.048 atm，而固体火箭冲压发动机补燃室内环境压力大都高于此值。因此，气相B2O3的凝结过程不能忽视。

1.3 数学模型及仿真方法

1.3.1 基本假设

(1)粒子为球形，非球形可通过形状因子修正[15]；

(2)气相组分满足理想气体状态方程；

(3)硼粒子内部温度均匀，外部为层流流场；

(4)燃烧过程中粒子表面没有凝相氧化物覆盖；

(5)硼粒子燃烧过程远离熄火边界。

其中，形状因子FP=rPAP/VP,rP、AP、VP分别为粒子半径、外表面积和体积，测量粒子物理尺寸的方法较多，如筛分法、显微镜法、沉降法等，具体可见文献[15]。

1.3.2 控制方程组

控制方程组[16]：

式中ρ、φ、Γφ和Sφ等参数的具体取值见文献[6]和文献[16]，在此不再赘述。

气相满足状态方程：

(10)

式中p和T分别表示环境压力和温度；Ru为通用气体常数；Mi和Yi分别为组分i的摩尔质量和质量分数。

1.3.3 边界条件

(1)来流边界(r=r∞,0≤θ≤π/2)

uθ=U∞sinθ,ur=-U∞cosθ

T=T∞,p=p∞

YO2=YO2,∞,YB2O2,∞=0

YN2=1.0-YO2-YB2O2

(11)

式中U∞为来流相对速度。

(2)出口边界(r=r∞,π/2≤θ≤π)

(12)

出口法向速度ur利用连续方程求解得到。

(3)轴对称边界(θ=0,θ=π)

uθ=0

(13)

(4)粒子表面边界

粒子表面切向气流速度：

uθ,w=0

(14)

法向速度满足组分质量守恒：

(15)

式中 下标w表示粒子相和气相的界面，则

(16)

硼粒子的总燃烧速率为

(17)

1.3.4 输运和热物性参数计算方法

输运参数和热物性参数的计算采用文献[17]中的方法，参数选取见表3。

本文将凝相组分B2O3l处理为具有一定扩散能力的“准气相组分”[7,18]，考虑其扩散能力相对较弱，参考气相B2O3的输运参数，估算凝相B2O3l的输运参数σ=3.0,ε/kB=240，其他输运参数仍采用气相B2O3的对应值。气相组分的比热容等热物性参数采用热力学数据库中数值[10]。

表3 气相组分的输运参数取值

1.3.5 数值仿真方法

数值仿真方法借鉴文献[6]和文献[16]的处理方法，采用有限体积离散方法和SIMPLE算法开发了一套C++程序，采用交错网格、ADI隐式交错迭代和块修正等技术，依次求解动量方程、能量方程和组分方程。大量数值计算结果表明，数值仿真域外半径取粒子半径20倍时，仿真边界已对粒子燃烧流场影响很小,θ方向和r方向均匀划分为200个网格，此时能同时满足计算精度和仿真效率要求。

2 结果与讨论

考虑到非稳态计算耗时非常大，本文仅开展稳态仿真，令时间步长dt=1.0×1030s，且只考虑氧化性气体为氧气组分，共涉及11个气相反应和前文所述的8个表面异相反应、1个凝结动力反应。目前固体火箭冲压发动机所用硼粒子直径多为1～20 μm，燃烧室压力多为3.0～5.0 atm。因此，本文主要针对该范围内硼粒子燃烧特性展开研究。

2.1 强迫对流下硼粒子燃烧流场组分分布

以相对速度U∞=10.0m/s、环境温度T∞=2000 K、粒子半径rp=10 μm、粒子温度Tp=2500 K、环境压力p=1 atm、环境氧气质量分数YO2,∞=0.80为例，图2(a)～(f)依次给出了该工况下仿真得到的粒子周围气相组分分布云图。

由图2可见，在强迫对流作用下，硼粒子周围的多组分燃烧火焰呈轴对称结构，由于各气相组分的扩散能力和气相反应速率不同，使得各组分在硼粒子周围流场中的分布差异较大。硼粒子与O2组分反应的中间气相产物B2O2、BO等主要在硼粒子的表面反应中生成，且以B2O2为主；而气相产物BO2、B2O3g以及凝相B2O3l组分则主要在气相火焰中生成，且以BO2为主。其中，B2O2、O2的分布也与文献[6]中基于总包反应研究所得的分布也相互印证。由于硼粒子的粒径相对很小，两相之间燃烧传质过程的总质量很小，使得空间气相反应相对表面反应较弱，加上强迫对流作用和各气相组分扩散能力的差异，微量的气相B2O3g和凝相B2O3l主要集中在粒子尾部，形成了微弱的“气相尾焰”结构。

图3(a)～(d)给出了该工况下硼粒子外表面上不同位置处各气相组分的质量分数分布。为便于比较，采用双纵坐标轴形式，且将数据表示为а×10n的形式，纵坐标为a的值。由图3可见，在强迫对流作用下，硼粒子外表面上O2组分的质量分数随着迎风角度θ的增加而减小，而B2O2、BO2等燃烧产物的质量分数则随着迎风角度θ的增加而增加。可见，强迫对流作用降低了环境中O2分子和O原子反应基向粒子外表面上的输运阻力和气相生成物向周围环境的排放阻力，使得反应生成物在来流携带作用下向粒子背风面聚集。从燃烧反应动力学过程分析可知，气相火焰中的B2O3g、BO2等生成物通过对流和扩散作用传递至粒子表面，再与硼粒子发生表面异相反应，生成BO等。结合图2可知，在该工况下，硼粒子燃烧的主要产物为BO2和B2O2，这2种气相产物将在燃烧室内进一步反应生成B2O3l。

2.2 硼粒子燃烧特性影响因素

在分析强迫对流下硼粒子周围燃烧流场组分分布的基础上，本节将重点研究相对速度、环境压力和氧气质量分数等因素对强迫对流下硼粒子燃烧特性的作用规律。

2.2.1 相对速度对硼粒子燃烧特性的影响

以环境温度T∞=2000 K、压力p=1 atm、氧气质量分数YO2,∞=0.80、粒子温度TP=2500 K为例，图4给出了几组不同相对速度下半径10 μm的硼粒子表面上各个基元反应的摩尔消耗速率变化，考虑到不同反应的速率量级差异，将几组反应速率放在两个图中显示，并以箭头指向标示。由图4可见，表面反应4的反应速率最大，其次为表面反应7、2、5、6和1、3的反应。表面反应4的反应速率较其他反应的速率高几个量级，说明硼与O2组分反应生成B2O2的过程占据了硼消耗的主导，这一点与Zhou等[12]的研究结论亦相同。由图4还可看出，表面反应2、3、4的反应速率随着相对速度的增加而增大，而其他几个表面反应的反应速率则随着来流速度的增加而减小。

分析可知，从气液两相之间反应的整个过程来讲，当来流相对速度增加时，来流中O2组分向粒子表面和反应生成物向周围环境的输运过程随之得以强化，气相生成物B2O3g、BO2等组分在来流携带作用下向粒子尾部聚集，且这些组分向周围环境的输运过程也随着相对速度的增加得以强化。因此，总的效果使得O2组分参与的硼粒子反应速率均随来流相对速度的增加而增加，而BO2、B2O3等气相组分参与的基元反应速率则随之降低。由于表面反应4的基元反应的速率相比其他几个反应大几个量级，因此，硼粒子总的燃烧速率随着来流相对速度的增加而增加。

2.2.2 粒子半径对硼粒子燃烧特性的影响

以相对速度U∞=100.0 m/s、环境温度T∞=2000 K、压力p=1 atm、氧气质量分数YO2,∞=0.23、粒子温度TP=2500 K为例，图5给出了几组不同粒子半径对应的硼粒子表面上基元反应摩尔消耗速率的变化。由图5可见，在不同粒径下，表面反应4的基元反应速率均为最大，说明在本节研究的粒子半径范围内，硼与O2反应生成B2O2的过程仍占据主导。由图5还可看出，当其他环境条件均相同时，硼粒子的粒径越大，粒子表面上各个反应的速率均越大。

分析认为，如果仅从两相之间的物理输运过程上讲，粒子半径越大，围绕硼粒子的各气相组分的输运阻力越大；但硼粒子的半径越大，粒子相与气相间的强迫对流作用越显著，这反而会促进来流气体中的O2等组分相粒子表面的输运过程。因此，两相之间的输运过程是两者竞相作用的结果。仿真结果表明，硼粒子的燃烧速率随着粒子半径的增加而增加，这说明粒径增大对强迫对流作用下两相输运过程的促进作用起主导作用。

2.2.3 环境压力对硼粒子燃烧特性的影响

以相对速度U∞=100.0 m/s、环境温度T∞=2000 K、氧气质量分数YO2,∞=0.80、粒子温度TP=2500 K为例，图6给出了几组不同环境压力下半径10 μm的硼粒子表面上基元反应摩尔消耗速率的变化。由图6可见，在不同环境压力下，表面反应4的基元反应速率均为最大，说明在本节研究的环境压力范围内，硼粒子与O2反应生成气相B2O2的过程仍主导着硼粒子燃烧反应。由图6分析还可看出，硼粒子表面上各个反应的速率均随环境压力的增加而增大，且几乎呈线性增长趋势。

分析可知，环境压力的增加一方面使得环境中氧气浓度增大，另一方面将使得O2组分向硼粒子表面以及B2O3、BO2等燃烧生成物向周围环境的扩散输运阻力加大。因此，硼粒子外表面上气相产物B2O3、BO2等组分的质量分数随环境压力的增加而加大，同时使得这些组分参与的各基元反应的速率随之增加而增大，而O2组分参与的各基元反应速率也随着环境中氧气浓度的增加而加大。因此，硼粒子总的燃烧速率随环境压力的增加而加大。

2.2.4 环境中氧气质量分数对硼粒子燃烧特性的影响

以相对速度U∞=100.0 m/s、环境温度T∞=2000 K、环境压力p=1 atm、粒子温度TP=2500 K为例，图7给出了几组不同环境氧气质量分数下半径10 μm的硼粒子表面基元反应摩尔消耗速率的变化。由图7可见，硼粒子表面上各基元反应速率均随环境中氧气质量分数的增加几乎呈线性增加，且表面反应4的基元反应速率仍为最大，说明在本节研究的氧气质量分数范围内，硼粒子与O2组分反应生成气相B2O2的过程也占据了硼粒子反应消耗的主导。

分析可知，当环境中氧气质量分数增加时，通过强迫对流和扩散作用到达硼粒子外表面上的O2组分随之增加，使得O2组分参与的各基元反应速率均随之加大，而硼粒子外表面上的气相燃烧产物总量也随着O2组分反应速率的增加而加大，进而使得气相B2O3、BO2等组分参与的基元反应速率亦随之增加。因此，硼粒子总的燃烧速率随环境中氧气质量分数的增加而增加。

2.2.5 环境温度对硼粒子燃烧特性的影响

考虑到固体火箭冲压发动机燃烧室内温度分布的不均匀性，本节以来流相对速度U∞=100.0 m/s、环境压力p=1 atm、氧气质量分数YO2,∞=0.80、粒子温度TP=2500 K为例，图8给出了几组不同环境温度下半径10 μm的硼粒子表面反应摩尔消耗速率变化。

由图8可见，在本节研究的温度范围内，表面反应4的基元反应速率随环境温度的增加而增加，表面反应5和6的基元反应速率随着环境温度的增加而降低，而其他几个基元反应速率随环境温度的变化均不大。由于表面反应5和6的基元反应速率均相对较小，对硼粒子的整个燃烧过程影响相对很小，据此可认为,整体上环境温度对强迫对流下硼粒子的燃烧速率影响不大。

分析认为，环境温度主要影响各组分的输运过程和化学反应动力过程，硼粒子的燃烧消耗也是这些因素竞相作用的结果。由前文研究可知，环境温度对硼粒子的燃烧速率影响不大，尽管如此，从能量释放角度来讲，环境温度的升高并不利于放热反应的进行。因此，从这个角度上讲，一旦硼粒子被点燃，为促进硼能量的高效释放，往往需要富氧环境。

3 结论

(1)由于各气相组分反应速率和组分扩散能力不同，使得各组分在流场中的分布差异较大，气相中间产物B2O2、BO等主要在硼粒子表面反应中生成，而气相产物B2O3、BO2等以及凝相B2O3s组分则主要在气相火焰中生成。在强迫对流环境下，微量的气相B2O3和凝相B2O3s集中在粒子尾部，形成了“气相尾焰”。

(2)在本文研究的环境范围内，随着来流相对速度、粒子半径、环境压力、环境中氧气质量分数的增加，硼粒子总的燃烧反应消耗速率加大。

(3)当硼粒子远离熄火边界时，环境温度对硼粒子表面上各基元反应的反应速率影响较小。

[1] Young G,Sullivan K,Michael R Zachariah.Combustion characteristics of boron nanoparticles[J].Combustion and Flame,2009,156(2):322-333.

[2] Yanan Gan,Yi Syuen Lim,Li Qiao.Combustion of nanofluid fuels with the addition of boron and iron particles at dilute and dense concentrations[J].Combustion and Flame,2012,159(4): 1732-1740.

[3] 张鹏,洪延姬,丁小雨,等.气相硼燃烧的化学动力学分析[J].推进技术,2015,36(10):1582-1587.

[4] 敖文,杨卫娟,汪洋 等.气流速度对晶体硼粒子热氧化及点火燃烧特性的影响[J].固体火箭技术,2013,36(4):511-515.

[5] Hu Jian-xin,Xia Zhi-xun,Wang De-quan.Study on the burning rate of boron particles under forced convection conditions in secondary chamber of ducted rocket[J].Journal of Solid Rocket Technology,2007,30(1): 21-25.

[6] 方传波,夏智勋,胡建新,等.强迫对流下硼粒子燃烧特性影响因素研究[J].航空学报,2015,36(2):492-500.

[7] Zhou W.Numerical study of muti-phase combustion: ignition and combustion of an isolated boron particle in fluorinated environments[D].Princeton University,Princeton,1998.

[8] Yetter R A.Kinetics of high-temperature B/O/H/C chemistry [J].Combustion and Flame,1991,83(1-2):43-62.

[9] Pasternack L.Gas-phase modeling of homogeneous boron /oxygen /hydrogen /carbon combustion[J].Combustion and Flame,1992,90(3-4):259-368.

[10] Malcolm W C Jr.NIST-JANNAF thermochemical tables,4th ed,Parts I and II[J].Journal of Physical and Chemical Reference Date,1998,9(1):177-192.

[11] Brown R C.Kinetic model of liquid B2O3gasification in a hydrocarbon combustion environment:Heterogeneous surface reactions[J].International Journal of Chemical Kinetics,1991,23(11):957-970.

[12] Raymond B Edelman,Constantino Economos,John Boccio.An analytical and experimental study of some problems in two-phase flows involving mixing and combustion with application to B-O-H-N system[C]//AIAA Reaction Turbulent Flows Conference San Diego,United States,1970.

[13] Girshick S L,Chiu C P.Kinetic nucleation theory:A new expression for the rate of homogeneous nucleation from an ideal supersaturated vapor[J].Journal of Chemical Physics,1990,93(2): 1273-1277.

[14] Shpilrain E E,Yakimovich K A,Tsitsarkin A F.Surface tension of liquid boric oxide at up to 2100℃[J].High Temperature,1974,22(1):77-82.

[15] 葛庆仁.气固反应动力学[M].北京:原子能出版社,1991:57-62;228-243.

[16] Daniel N Pope.Numerical simulation of convective fuel droplet vaporization and combustion in a low pressure zero-gravuty environment[D].University of Nebraska,Nebraska,Lincoln,2001.

[17] Reid R C,Prausnitz J M,Poling B E.The properties of gases and liquids[M].McGraw Hill,New York,1987: 256-364.

[18] Davidchuk E L,Dimitrov V I,Rafalovich M L,et al.Kinetics of boron combustion in dry air[J].Combustion,Explosion,and Shock Waves,1991,27(1): 45-52.

[19] Atsushi Makino,Norie Umehara.Combustion rates of graphite rods in the forward stagnation field of the high-temperature,humid airflow[J].Proceedings of the Combustion Institute,2003,132(4):743-753.

[20] 陶文铨.数值传热学(第二版)[M].北京:高等教育出版社,2001:90-132.