基于化學動力的強迫對流下硼粒子的燃燒特性

(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.