基于King模型的含硼富燃燃氣燃燒模擬①

馮喜平，劉 洋，任全彬，李進賢

(1.西北工業大學燃燒、熱結構與內流場重點實驗室，西安 710072 2.中國航天科技集團公司第四研究院，西安 710025)

0 引言

硼以高體積熱值優勢已成為理想的推進劑材料，在固沖發動機中表現出巨大的應用潛力[1]。但硼的高熔點(2 450 K)與高沸點(3 931 K)，使其點火困難。

國外自20世紀70～90年代就對硼粒子的點火機理進行了大量研究[2－6]，并建立半經驗公式，其中硼粒子點火模型以 King模型[3]和 L-W 模型[4]為代表。King點火模型基于O擴散至B-B2O3表面與B反應生成B2O3(l)的假設，考慮了氧化層的產生、消除對B點火和燃燒的影響。L-W模型認為，在高溫情況下，硼在氧化層界面上溶解，與B2O3發生反應生成聚合物(BO)n;該聚合物在氧化層中由內向外擴散，并在外界面上通過蒸發和化學反應消耗。Yeh和 Kuo[5－6]也研究了硼粒子點火過程，發現B顆粒在1 213 K就開始熔化，遠低于其熔化溫度2 450 K;在溫度升高的過程中，B不斷溶解于B2O3，并與B2O3反應生成了一種新的聚合物(BO)n。

國內霍東興等[7]采用King點火模型和L-W燃燒模型研究了某發動機補燃室內硼粒子粒徑對點火位置和點火效率的影響，得出2μm以下的硼粒子燃燒效率較高的結論。吳婉娥[8]建立了單個硼粒子的點火燃燒數學模型，通過數值方法研究環境溫度、氧分壓、水蒸氣分壓、粒徑對硼顆粒點火的影響，得出高溫環境、高水蒸氣分壓、小粒徑有利于硼粒子點火，但過高氧分壓會延長硼粒子點火時間的結論。李海波[9]對比了不同湍流模型、氣相燃燒模型對模擬結果的影響，并給出沖壓發動機補燃室中適用的一些模型[9]。

對于固沖發動機而言，含硼富燃燃氣二次燃燒是其關鍵技術。對于二次燃燒過程的計算，燃燒模型是關鍵。盡管國內外研究者進行了大量相關研究，但主要集中于硼粒子的點火燃燒模型等細節性和局部性研究環節，許多研究采用的數值模擬方法過分依賴于商業軟件自帶模型，對燃燒過程過分簡化，這些燃燒模擬普遍存在缺乏完整數值模型的問題(包括顆粒相燃燒模型、氣相燃燒模型等)。

針對上述問題，本文以燃燒實驗裝置為物理模型，選用King模型作為硼粒子點火燃燒模型，選用有限速率/渦耗散作為氣相燃燒模型，構建含硼富燃燃氣流動與燃燒過程的數學模型，結合典型實驗工況進行數值模擬，通過模擬結果與實驗結果對比，驗證模擬方法。

1 物理模型與計算方法

1.1 物理模型

燃燒室結構如圖1所示。燃燒室由燃氣發生器、燃燒室和尾噴管3部分組成。

圖1 燃燒室結構Fig.1 Structure of the chamber

為方便實現燃燒流場的測量，燃燒室采用橫截面為100 mm×100 mm、長635 mm的長方體結構;數值計算時，為減少計算量，將燃燒室簡化為等流通面積的圓柱體，直徑113 mm。其他參數見圖1。整個燃燒室采用軸對稱結構。

實驗裝置工作原理如下:含硼富燃燃氣由燃氣入口進入燃燒室，空氣由位于燃燒室前部的環形空氣入口噴入燃燒室，二者進入燃燒室后沿各自流線運動;在燃氣和空氣組分梯度的作用下相互擴散、燃燒，在燃燒室形成擴散火焰;燃燒后產物由噴管噴出。燃氣流量通過燃氣發生器噴管音速喉道控制，空氣流量通過空氣管道中的音速喉道控制。

1.2 數學模型

1.2.1 基本假設和控制方程

基本假設:

(1)燃燒室內燃氣流動為準定常流動，與外界無熱交換;

(2)忽略燃氣各組分之間的輻射作用，忽略重力的影響;

(3)燃燒室燃氣為理想氣體，服從理想氣體狀態方程;

(4)不可燃氣體假設為N2。

控制方程:

式中 各項表達式見文獻[10]。

1.2.2 湍流模型

文獻[9]通過對比不同湍流模型對模擬結果的影響，得出RNG k-ε湍流模型更適用于本物理模型內二次燃燒的結論。因此，湍流模型選取RNG k-ε模型。具體方程及說明見參考文獻[11]。

1.2.3 燃燒模型

(1)氣相燃燒模型

氣相燃燒模型采用通用有限速率/渦耗散模型，分別計算Arrhenius和渦耗散反應速率，凈反應速率取兩個速率中較小的。Arrhenius反應速率作為一種動力學開關，阻止反應在火焰穩定之前發生。一旦火焰被點燃，渦耗散速率通常會小于Arrhenius反應速率，這樣可延遲計算點火的開始，提供點火時間，更加符合實際[9]。

(2)碳顆粒燃燒模型

碳顆粒燃燒模型采用動力學/擴散控制反應速率模型。該模型假定表面反應速率同時受到擴散過程和反應動力學的影響。該模型的擴散速率常數為

動力學反應速率常數為

二者加權得反應速率為

式中 pox為顆粒周圍的氣相氧化劑分壓;R為考慮了焦碳內表面的反應及其擴散的化學(動力學)反應速率常數。

(3)硼顆粒相燃燒模型

硼顆粒相點火燃燒模型選用King模型[3]。

硼粒子表面覆蓋有致密的氧化層，阻止潔凈硼顆粒與氧氣直接接觸，見圖2。

硼粒子點火經歷3個過程:

(1)表面B2O3蒸發;

(2)B2O3與H2O反應生成HBO2氣體;

(3)O2擴散至B-B2O3邊界，與潔凈B發生反應，生成新的B2O3。

圖2 King模型中硼粒子點火機理Fig.2 IgnitionofboroninKingmodel

King指出，只有B2O3的消耗速率大于生成速率時，才能認定B點火成功，進入燃燒階段。

其數學描述[3]如下:

式中 dB、δ、f、TP分別為硼粒子直徑、氧化層厚度、熔化分數、硼粒子溫度;為點火過程中的能量變化，包括顆粒反應釋放的熱量及其與周圍環境之間的對流換熱和輻射換熱;T、TRAD分別為環境溫度和周圍輻射溫度;RB、RE、RH分別為 B消耗速率、B2O3蒸發速率、B2O3與水的反應速率;pO2和pH2O分別為環境氣體中氧氣和水蒸氣的分壓;ΔHM、p、XH2O、XO2分別為 B 的熔化熱、環境壓力、H2O的摩爾分數、O2的摩爾分數。

燃燒階段，潔凈硼粒子反應速率遵從有限速率模型:

式中 k為Arrhenius定律計算出的反應速率;dB為當前潔凈硼粒子粒徑。

為了完整、準確描述硼顆粒相點火燃燒過程，采用Fluent擴展應用UDF編寫King模型程序。程序編制中，采用如下假設:

(1)假定硼粒子表面存在1%厚度的氧化層;

(2)認定H2O(g)存在能加快點火過程，但不考慮HBO2的生成反應;

(3)認為顆粒相中只存在一種成分B，且未點火前B顆粒質量、密度、粒徑不發生變化，點火完成前B2O3的生成和消耗速率進行單獨計算，不加入迭代方程中;

(4)點火過程完成后，硼粒子遵從Fluent自帶顆粒表面反應模型。

1.3 網格劃分

圖3為計算網格。

圖3 計算網格Fig.3 Meshofthechamber

采用四邊形與三角形網格相結合，部分加密的方法對模型進行網格劃分。其中，1區域為燃氣通道，采用三角形網格并加密;2區域為燃燒室，采用四邊形網格，中部區域加密;3區域為尾噴管，采用四邊形網格，整體加密。總網格數量53854。

1.4 邊界條件

邊界條件見表1。

表1 邊界條件Table1 Boundaryconditions

燃氣流量為 0.03kg/s，空燃比為 10。

顆粒相采用離散相、面入射、隨機軌道模型，初始粒徑為10－5m，速度為100m/s，溫度1800K。

燃氣組分通過燃氣發生器熱力計算得到，結果見表2。燃氣中，氣相組分32.2%，顆粒相67.8%。

表2 富燃燃氣主要組分Table 2 M ain component of fuel-rich gas

2 計算結果與討論

2.1 燃燒室流場分布

圖4為模擬得到的燃燒室流場分布。

圖4(a)為燃燒室壓力等值線分布。由圖可見，整個燃燒室壓力均勻，平均壓力為0.48 MPa;燃燒室出口處，由于拉瓦爾噴管作用，壓力逐漸降低。

圖4(b)給出了燃燒室內溫度分布。最低溫度出現在空氣入口處，為573 K;高溫區域(2 800 K)分別出現在燃燒室前端燃氣與空氣接觸面(210mm處)、燃燒室1/3(280～370 mm處)和燃燒室中后部靠近軸線區域(400 mm～尾部)。燃燒室前端燃氣與空氣接觸面處，燃氣與氧氣發生劇烈的化學反應。因此，鄰近微小空間內溫度急速上升，形成第一個2 800 K高溫區域;前部區域空氣、燃氣速度梯度與溫度梯度較大，在接觸面附近形成嚴重湍流，熱量很快向周圍傳遞，2 800 K高溫區消失。

圖4 燃燒室流場分布Fig.4 Distribution of flow field in the chamber

燃燒室280 mm處，氣相燃料基本燃燒完畢，少量硼粒子點火成功，開始與氧氣反應，形成第2個高溫區域;由于僅有少量硼粒子在此區域點火成功，且很快完全燃燒，其他硼粒子尚未開始反應。因此，第2個2 800 K高溫區持續90 mm后就消失。燃燒室400 mm處，氧氣向中心軸線區域擴散，大量硼粒子點火成功，開始劇烈的燃燒反應，與前兩個高溫區域不同，第3個2 800 K高溫區由空-燃接觸面擴展至中心軸線處，表明氧氣已擴散至中心區域，與大量處于軸線附近的硼粒子發生反應，釋放巨大熱量;高溫區域一直延展到燃燒室尾部，說明硼粒子燃燒反應持續進行。

從整個溫度場結構來看，隨著流動進行，燃氣向壁面方向擴散，與外部氧氣發生反應，火焰寬度逐漸增大;尾部區域，溫度場分層現象仍然嚴重，中心溫度2 800 K，壁面溫度僅為1 400 K，二者差值達1 400 K。溫度場呈錐形擴散火焰形狀。

圖4(c)為燃燒室速度分布圖。速度分布呈現射流特征。在頭部區域靠近軸線處，中心射流速度高達140 m/s，而外部流場速度僅10 m/s左右。二者接觸面速度梯度大，分層明顯。由于氣相的粘性作用、燃氣與空氣之間分子擴散產生的動量交換，中心氣流速度逐漸降低，外部氣流速度逐漸升高。燃燒室200mm左右位置，射流核心區消失，射流主體區域出現，此時軸線中心速度依然高達80 m/s。而在燃燒室中后部510 mm處，速度梯度已經很小，流場速度趨于均勻化。如果燃燒室足夠長，流動到一定階段后，整個流場速度將完全均勻，變為層流運動。

圖4(d)給出了O2的分布。在燃燒室前部區域，氧化劑由空氣入口處進入，O2質量分數最大;燃氣中不含O2，故O2質量分數最小。隨著流動進行，燃氣向外擴散，二者在火焰面處進行燃燒，消耗氧氣，高含氧區域縮小，氧組分等值線呈現向外彎曲形狀;隨著流動向下游進行，氧氣進一步消耗。在燃燒室軸線附近區域，盡管存在氧氣向中心區域的擴散，但由于燃燒的進行，氧氣在火焰面位置大量消耗，因此燃燒室軸線附近區域氧含量非常小。

圖4(e)給出了B2O3的分布。在燃燒室頭部不存在B2O3。這是由于硼粒子剛進入燃燒室，尚未點火成功。210 mm處，B2O3開始生成，標志部分B粒子點火成功，并開始與O2發生反應。B2O3生成位置并非在中心軸線處，而是處于氧氣與燃氣的接觸面附近區域。這是由于燃氣中可燃氣體成分在頭部區域大量消耗O2，處于中心區域的硼粒子幾乎很難接觸到氧氣。只有少量硼粒子偏離中心，射入燃氣與氧氣反應的區域中，在高溫、高氧氣分壓、高水蒸氣分壓條件下，很快點火成功，生成B2O3。400 mm處，B2O3濃度明顯提高，即大量硼粒子開始與O2反應。這一點在圖4(b)溫度分布中也有所反映:400mm處，燃燒室溫度第三次達到峰值2 800 K。在這一階段，中心區域也生成大量B2O3，這是由于燃氣已經消耗完畢，O2擴散至中心，與軸線附近的B粒子發生反應。B2O3最高濃度出現在燃燒室尾部，說明B顆粒與O2的劇烈反應一直持續到燃燒室尾部;O2質量分數分布也支持這一結論，中后部區域中心軸線附近O2濃度仍為0，而靠近壁面處O2濃度基本達到一致，說明中心軸線區域持續發生燃燒反應。

2.2 硼粒子點火距離

硼粒子進入燃燒室到硼粒子劇烈燃燒的距離為硼粒子的點火距離，模擬結果見圖5～圖7。

圖5 B顆粒群質量變化曲線Fig.5 M ass-length of boron particles

圖6 B2O3質量分數分布圖Fig.6 M ass fraction of B2 O 3

圖7 硼粒子軌跡圖Fig.7 Boron particle tracks

圖5為燃燒室內不同入射位置的硼粒子質量變化曲線。圖6為燃燒室頭部B2O3的質量分數分布圖(出現B2O3即認定點火成功)。圖7為不同入射位置硼粒子軌跡。由于入射位置不同，硼顆粒在燃燒室中的軌跡也不同。靠近軸線方向的硼粒子周圍氧氣與水蒸氣濃度較低，溫度上升慢，點火距離較長。靠近邊界層的硼粒子能更早地接觸氧氣與水蒸氣，且溫度較高，因此點火距離短。以圖5中硼顆粒1、2為例，其軌跡如圖7中所示。顆粒1最先點火成功，但燃速很慢，300 mm處突然燃速加快，并且很快完全燃燒，對應圖7中顆粒1此前一直處于中心軸線處，溫度、氧分壓、水蒸氣分壓低下;在300 mm處開始偏離中心軸心，射入高溫反應區(溫度圖見圖4(b))，因此很快燃燒完畢，軌跡線消失。顆粒2也是由于偏離中心軸線，得以快速點火燃燒。

由圖5和圖6得出，硼粒子在210mm處質量發生變化，即成功點火，除去燃氣入口及延長段147 mm距離，得出硼粒子點火距離為63 mm。大部分硼粒子在400 mm處才點火成功，開始劇烈燃燒，這種現象也可說明圖4(b)中400 mm處2 800 K高溫區產生的原因。由圖5和圖7得出，大量硼粒子沒有完全燃燒，隨氣體一起流出尾噴管。

圖8是實驗者捕捉的穩定燃燒階段燃燒室中硼與氧氣反應時明亮火焰圖像，即硼粒子點火位置。經過處理計算，可得出穩定階段點火距離約為50 mm。數值模擬中硼粒子點火距離為63 mm，二者數值接近。

圖8 硼粒子燃燒火焰圖Fig.8 Flame of boron combustion

2.3 模擬結果與實驗結果對比

為了驗證King點火燃燒模型的準確性，將數值模擬結果與實驗參數對比，結果見表3。

表3 模擬結果與實驗結果對比Table 3 Com parison of simulation and experim ent result

由表3可得出，使用King模型進行數值模擬得出的燃燒室壓力為 0.48MPa，實驗結果為 0.43MPa，誤差11.6%;點火距離為63 mm，實驗結果為 50 mm，誤差26%。可見，使用King模型的模擬結果較好地預測了實驗工況。

3 結論

(1)進行含硼富燃燃氣的二次燃燒時，粒子相采用King模型點火與燃燒，氣相燃燒采用有限速率/渦耗散模型，湍流采用RNG k-ε模型，模擬結果與實驗結果對比說明，該模擬方法較好預測了燃燒過程。

(2)采用Fluent軟件UDF功能，依據粒子相點火燃燒模型自編程序，克服了目前數值模擬中燃燒過程過于簡化問題，較準確地表述了硼粒子點火燃燒機理和燃燒過程，模擬結果更加可靠。

(3)數值模擬表明，粒子相在距燃氣出口一定距離點火并開始燃燒，燃燒火焰頭部呈錐形火焰結構。

[1] 湯中權，銳鑫.含硼固體推進劑[J].化學推進劑與高分子材料，1998，3:40-42.

[2] 王寧飛，關大林，范紅杰.硼顆粒點火和燃燒研究進展[J].含能材料，2001，9:86-89.

[3] King M K.Ignition and combustion of boron particles and clouds[J].Journal of Spacecraft and Rockets，1982，19(4):294-306.

[4] Li SC，Williams F A.Ignition and combustion of boron particles[M].Combustion of Boron-Based Solid Propellants and Solid Fuels.Boca Raton，CRC Press，1993:248-271.

[5] Yeh C L，Kuo K K.Ignition and combustion ofboron particles[J].Progress in Energy and Combustion Science，1996，22(6):511-541.

[6] Kuo K K.Boron ignition model[M].NPU，1998.

[7] 霍東興，陳林泉，劉霓生，等.硼粒子直徑對點火位置及燃燒效率的影響研究[J].固體火箭技術，2004，27(4):272-275.

[8] 吳婉娥，裴明敬，郭耳鈴，等.硼粒子在固沖環境中點火過程影響因素的數值模擬[J].火炸藥學報，2008，6:79-82.

[9] 李海波.含硼燃氣與空氣擴散燃燒模型研究[D].西安:西北工業大學，2013.

[10] 方丁酉.兩相流動力學[M].長沙:國防科技大學出版社，1988.

[11] Yakhot，Orszag S A.Renormalization group analysis of turbulent[J].J.Sci.Comput.，1986，1:3-9.