楊 樺，姜 毅，陶倩楠，張曼曼，黃陽陽

(北京理工大學(xué) 宇航學(xué)院，北京 100081)

0 引言

固體火箭發(fā)動機(jī)作為一種高性能動力裝置具有結(jié)構(gòu)簡單、可靠性強(qiáng)、使用方便等諸多優(yōu)點，被廣泛用于火箭、導(dǎo)彈的發(fā)射和推進(jìn)。由于固體火箭發(fā)動機(jī)在工作過程中燃燒室內(nèi)的推進(jìn)劑燃燒不充分，未燃燒完全的燃?xì)馀懦龊笈c空氣接觸產(chǎn)生二次燃燒，使尾焰流場更加復(fù)雜，進(jìn)而影響安全范圍的確定以及隔熱、排導(dǎo)裝置的合理設(shè)計。

為模擬真實的燃?xì)饬鲌觯芯空邆儼l(fā)展了關(guān)于多種不同推進(jìn)劑組分的化學(xué)反應(yīng)理論，以研究尾焰流場復(fù)燃現(xiàn)象。黃振宇等[1]研究了基于非平衡化學(xué)反應(yīng)流動的火箭燃?xì)饬鲌?，發(fā)現(xiàn)化學(xué)反應(yīng)集中在混合層內(nèi)。姜毅等[2]根據(jù)基于一維等熵流動獲得的噴管出口參數(shù)和燃?xì)饨M分含量，列寫組分輸運方程，考慮化學(xué)反應(yīng)動力過程，計算得出了復(fù)燃尾噴焰流場的結(jié)構(gòu)。馬艷麗等[3]采用有限速率化學(xué)反應(yīng)模型模擬了燃?xì)馍淞鲗Πl(fā)射平臺的熱沖擊和動力沖擊效應(yīng)，發(fā)現(xiàn)平臺離發(fā)動機(jī)噴口距離不同時壓力分布變化不大，發(fā)射平臺越遠(yuǎn)，平臺上溫度變化越小。賀衛(wèi)東等[4]基于9組分10計算反應(yīng)的H2/CO氧化反應(yīng)體系分析了不同海拔高度下發(fā)射時燃?xì)馍淞鲗Φ孛婕鞍l(fā)射裝置的影響，研究發(fā)現(xiàn)，發(fā)射海拔增加，射流速度降低，欠膨脹度增加，溫度升高。

本文借助計算流體力學(xué)技術(shù)對燃?xì)鈱?dǎo)流器流場進(jìn)行數(shù)值模擬，將模擬發(fā)射過程中燃?xì)饣瘜W(xué)反應(yīng)的方法應(yīng)用于實際的流場分析，參考文獻(xiàn)[5]把基元反應(yīng)總數(shù)量為20的C—H—O—N化學(xué)反應(yīng)體系引入流場計算過程中，對發(fā)射過程中安全范圍的確定及導(dǎo)流器的合理設(shè)計有著重要的參考意義。

1 計算模型

1.1 湍流模型

Realizablek-ε模型[6]是由Shih等于1995年提出的一種帶旋流修正的模型，在處理旋轉(zhuǎn)流動、強(qiáng)逆壓梯度的邊界層流動以及流動分離等現(xiàn)象時具有非常好的表現(xiàn)。Realizablek-ε模型中湍動能k的輸運方程為

Gb-ρε-YM+Sk

(1)

模型中湍動能ε的輸運方程為

(2)

其中

C1=max[0.43，η/η+5]

η=Sk/ε

式中μt為湍流粘性系數(shù)；Gk為由于平均速度梯度產(chǎn)生的湍動能；Gb為由于浮力產(chǎn)生的湍動能；YM為可壓縮湍流中的波動擴(kuò)張對總耗散率的貢獻(xiàn)；Sk和Sε為自定義源項；C2和C1ε為常數(shù)，C2=1.9，C1ε=1.44；σk、σε分別為k和ε的普朗特數(shù)，σk=1.0，σε=1.2。

1.2 化學(xué)反應(yīng)模型

(1)組分輸運方程

使用質(zhì)量守恒方程來描述尾焰中的復(fù)雜組分：

(3)

(4)

式中Di,m為第i種組分的質(zhì)量擴(kuò)散系數(shù)；DT,i為該組分的溫度擴(kuò)散系數(shù)；Sct為湍流施密特數(shù)；T為溫度。

(2)基元反應(yīng)

化學(xué)反應(yīng)過程中，一個或多個組分僅通過一個過渡態(tài)，一步反應(yīng)直接生成產(chǎn)物的過程稱為基元反應(yīng)[7]?；磻?yīng)的反應(yīng)速率與產(chǎn)物濃度、環(huán)境的溫度直接相關(guān)。一般用阿雷尼烏斯公式來描述基元反應(yīng)的化學(xué)反應(yīng)速率，對于雙基元反應(yīng)，其溫度為T時，正向化學(xué)反應(yīng)速率為

(5)

式中 [A]和[B]為反應(yīng)物A、B的摩爾濃度；A稱為前指因子，其量綱與反應(yīng)級數(shù)有關(guān)；b為溫度指數(shù)，無量綱；Ea為反應(yīng)的活化能；Ru為通用氣體常數(shù)。

以吉布斯自由能梯度等于零作為判斷某時刻反應(yīng)達(dá)到平衡態(tài)的標(biāo)志，即平衡狀態(tài)下吉布斯自由能達(dá)到最小。根據(jù)該判據(jù)在已知化學(xué)反應(yīng)表達(dá)式的情況下，計算出化學(xué)動力學(xué)中化學(xué)平衡常數(shù)，由正向化學(xué)反應(yīng)速率可以推知逆向化學(xué)反應(yīng)速率，反之亦可。

(3)有限速率反應(yīng)模型

式(3)中化學(xué)反應(yīng)的凈產(chǎn)生率Ri可由有限速率化學(xué)反應(yīng)模型描述：

(6)

對于第r個化學(xué)反應(yīng)，其當(dāng)量表達(dá)式為

(7)

由此可得出在第r個基元反應(yīng)中第i個組分單位體積物質(zhì)的量的生成率[8]。

(8)

其中，Cj,r為該物質(zhì)的摩爾濃度；Г為第三體的作用對于反應(yīng)生成(消耗)速率的影響，其表達(dá)式為

(9)

式中Ck為第三體中第k種組分的摩爾濃度；γk,r為在此基元反應(yīng)中該組分的化學(xué)計量數(shù)。

(10)

(11)

所有反應(yīng)的單位體積反應(yīng)焓生成率之和取負(fù)，即為流場中的由于化學(xué)反應(yīng)產(chǎn)生的能量SR：

(12)

將該能量項作為化學(xué)反應(yīng)源項帶入能量守恒方程，實現(xiàn)流場與動態(tài)化學(xué)反應(yīng)的耦合計算。

2 湍流模型驗證

本文的主要研究對象為固體火箭發(fā)動機(jī)尾部燃?xì)馍淞鳑_擊壁面形成的高溫高壓流場，選用的湍流模型為Realizablek-ε模型。本節(jié)通過經(jīng)典的燃?xì)饬髌桨鍥_擊問題對選用的湍流模型進(jìn)行驗證。

根據(jù)文獻(xiàn) [9]中數(shù)據(jù)建立沖擊模型，如圖 1所示。燃燒室壓強(qiáng)為2.1 MPa，總溫為3000 K，環(huán)境壓強(qiáng)為101 325 Pa，環(huán)境溫度為300 K，射流軸線與平板夾角為55°，平板壁面選為無滑移絕熱壁面，噴管出口中心到平板的距離為噴管出口直徑的2倍，選擇1/2模型進(jìn)行計算。

圖1 沖擊模型示意圖

圖 2為文獻(xiàn)[9]中實驗結(jié)果與本文流場仿真結(jié)果對比。射流在沖擊平板處產(chǎn)生燃?xì)饷芏确逯狄约榜R赫數(shù)谷值。對比仿真結(jié)果與實驗結(jié)果，仿真結(jié)果較好地反映了流場物理量的真實分布情況。

圖2 仿真計算結(jié)果與實驗對比

本節(jié)還采用了標(biāo)準(zhǔn)k-ε模型、Realizablek-ε模型對該沖擊流場進(jìn)行進(jìn)一步的計算比較。對稱面處平板上與噴管軸線距離不同的各點的壓強(qiáng)分布如圖 3所示，橫坐標(biāo)是對稱面處平板上各點到噴管中軸的水平距離與噴管出口半徑的比值，縱坐標(biāo)是相應(yīng)位置的壓強(qiáng)與總壓之比。

圖3 平板與對稱面的交線上的壓強(qiáng)分布

基于兩種湍流模型都計算得到了平板上的壓強(qiáng)分布情況，在x軸正向處，出現(xiàn)壓強(qiáng)兩個峰值，壓強(qiáng)最大值點離噴管出口更近。Realizablek-ε模型為湍流粘性增加了計算公式，而不作常值處理；該模型中的耗散率輸運方程更為精確，更加符合湍流的物理特性。對比結(jié)果說明該模型更好地計算出了壓強(qiáng)的兩個峰值點，數(shù)值也與實驗數(shù)據(jù)更接近，因此采用Realizablek-ε模型對燃?xì)鉀_擊流場進(jìn)行數(shù)值計算是合理的。

3 數(shù)值計算及結(jié)果分析

研究使用的固體推進(jìn)劑為RDX-CMDB[10]，該推進(jìn)劑主要由硝化纖維素(NC)、硝化甘油(NG)和黑索金(RDX)組成。

引燃推進(jìn)劑后，固相組元在預(yù)熱區(qū)內(nèi)即率先由于受熱而氣化，氣化形成的氣體進(jìn)入表面分解區(qū)。表面分解區(qū)部分溫度較高，氣體成分在高溫作用下分解，組分中的N—O化學(xué)鍵和N—N化學(xué)鍵首先發(fā)生斷裂，重新組合生成二氧化氮和醛基生成物。嘶嘶區(qū)處由于反應(yīng)放熱的作用，溫度繼續(xù)升高。表面分解區(qū)的反應(yīng)產(chǎn)物進(jìn)入嘶嘶區(qū)之后繼續(xù)反應(yīng)。分解反應(yīng)中的醛基產(chǎn)物被二氧化氮氧化，生成一氧化氮，一氧化碳和水。這些燃燒產(chǎn)物進(jìn)入發(fā)光區(qū)之后進(jìn)一步發(fā)生氧化與還原反應(yīng)，生成H2、CO2以及N2等生成物，在反應(yīng)過程中繼續(xù)放出大量熱量。反應(yīng)區(qū)溫度持續(xù)升高，高溫作用提高了這些氧化還原反應(yīng)的正向及逆向反應(yīng)速率，使得反應(yīng)的中間產(chǎn)物濃度也較高。由于壓強(qiáng)升高導(dǎo)致表面分解區(qū)、嘶嘶區(qū)和暗區(qū)的厚度大大降低，甚至薄到難以分辨，因此可以忽略大分子的醛基產(chǎn)物，認(rèn)為其存在的范圍和時間都小到可不計，推進(jìn)劑的燃燒產(chǎn)物在燃面附近很快就可以達(dá)到化學(xué)平衡狀態(tài)，因此可以根據(jù)最小吉布斯自由能理論計算燃燒過程的產(chǎn)物組成與絕熱溫度。詳細(xì)的反應(yīng)體系組成和反應(yīng)正向反應(yīng)速率系數(shù)見表 1。表中正向反應(yīng)速率包括前指因子A(其量綱[m,kmol,s]與反應(yīng)級數(shù)相關(guān))、溫度指數(shù)b以及反應(yīng)活化能Ea。

表1 反應(yīng)體系及正向反應(yīng)速率表

因為基元反應(yīng)R2和R9的反應(yīng)速率系數(shù)在對數(shù)坐標(biāo)系中非線性特質(zhì)較強(qiáng)，為方便進(jìn)行描述，采用多條曲線的疊加來表示，由每一條曲線得出相應(yīng)的反應(yīng)速率系數(shù)，進(jìn)行加和計算得出該溫度下的反應(yīng)速率系數(shù)。在一定的溫度范圍內(nèi)可以通過擬合得出能夠符合工程應(yīng)用精度的一組表達(dá)式?；磻?yīng)R5、R6、R7、R8、R11、R12和R20是三分子基元反應(yīng)，其中第三體的組分構(gòu)成如式(13)所示。

{[M1]=[H2]+6.5[H2O]+0.4[O2]+0.4[N2]+0.75[CO]+1.5[CO2]

[M2]=2[H2]+6[H2O]+[N2]+1.5[CO]+2[CO2]

[M3]=[H2]+2.5[M2O]+0.4[O2]2+0.4[N2]+0.75[CO]+1.5[CO2]

(13)

本節(jié)基于上表所描述的基元反應(yīng)總數(shù)量為20的C—H—O—N化學(xué)反應(yīng)體系對尾焰的化學(xué)反應(yīng)過程進(jìn)行分析，數(shù)據(jù)來源于文獻(xiàn) [11-13]。氣體組分包括CO2、CO、H2O、H2、NO2、NO、N2和O2，反應(yīng)中間產(chǎn)物氣態(tài)自由基包括OH、H、O、N，一共12種基本反應(yīng)組分。

3.1 三維發(fā)射試驗平臺流場數(shù)值計算說明

本節(jié)的研究對象是沖擊導(dǎo)流器的燃?xì)饬鲌?。?dǎo)流器模型如圖 4所示。導(dǎo)流器放置于發(fā)射裝置后部，在導(dǎo)彈發(fā)動機(jī)尾噴管下方。具體發(fā)射架整體模型如圖 5所示。壁面邊界設(shè)置為絕熱，流場入口設(shè)置為壓力入口，參數(shù)按照某型發(fā)動機(jī)進(jìn)行設(shè)置，發(fā)動機(jī)內(nèi)部氣體設(shè)置為混合燃?xì)?。具體如表 2所示。流場出口設(shè)置為壓力出口，參數(shù)按照標(biāo)準(zhǔn)大氣環(huán)境設(shè)置，即壓力為101 325 Pa，溫度為300 K。

圖4 導(dǎo)流器模型

圖5 發(fā)射架整體模型

表2 流場入口邊界條件

3.2 三維發(fā)射試驗平臺有導(dǎo)流器的流場結(jié)果分析

圖 6給出了化學(xué)反應(yīng)射流流場中的化學(xué)反應(yīng)熱分布?？煽闯觯鲌鲋写嬖谖鼰岱磻?yīng)和放熱反應(yīng)，放熱反應(yīng)主要發(fā)生在噴管喉部、導(dǎo)流器上的燃?xì)舛逊e區(qū)、地面上的燃?xì)舛逊e區(qū)以及燃?xì)馀c空氣的混合區(qū)域。

圖6 化學(xué)反應(yīng)熱

圖 7對比了計算結(jié)果中化學(xué)凍結(jié)流流場與化學(xué)反應(yīng)射流流場中幾種主要氣體組分在模型對稱面上的質(zhì)量分?jǐn)?shù)分布?？梢?，相對于凍結(jié)流流場，化學(xué)反應(yīng)射流流場中在射流與空氣接觸的混合區(qū)域以及導(dǎo)流器導(dǎo)流面、地面等燃?xì)舛逊e區(qū)域CO、H2的質(zhì)量分?jǐn)?shù)降低，CO2、H2O的質(zhì)量分?jǐn)?shù)增加。這是由于在燃?xì)馀c空氣的混合區(qū)域內(nèi)CO、H2與空氣中的O2發(fā)生化學(xué)反應(yīng)，產(chǎn)生CO2、H2O，在導(dǎo)流面、地面等燃?xì)舛逊e區(qū)域內(nèi)燃?xì)饨M分之間相互作用發(fā)生化學(xué)反應(yīng)產(chǎn)生CO2、H2O。

(a) H2 (b) CO

通過以上分析表明，在化學(xué)反應(yīng)射流流場中，在燃?xì)馀c空氣的混合區(qū)域燃?xì)饨M分與空氣成分發(fā)生化學(xué)反應(yīng)產(chǎn)生復(fù)燃現(xiàn)象(與文獻(xiàn)[14-15]所得結(jié)論一致)，在導(dǎo)流面的燃?xì)舛逊e區(qū)域、地面的燃?xì)舛逊e區(qū)域燃?xì)饨M分之間發(fā)生化學(xué)反應(yīng)產(chǎn)生復(fù)燃現(xiàn)象。

圖 8圖 9所示分別對比了計算結(jié)果中化學(xué)凍結(jié)流流場與化學(xué)反應(yīng)射流流場在模型對稱面上的溫度分布以及在導(dǎo)流面、地面上的溫度分布。圖 10給出了兩種工況計算結(jié)果在模型對稱面處燃?xì)鉀_擊壁面上的溫度分布曲線，橫坐標(biāo)是模型對稱面處燃?xì)鉀_擊壁面上的各點到噴管中軸的水平距離，射流在對稱面處地面上的流動方向為正方向，縱坐標(biāo)為溫度。

圖8 對稱面處燃?xì)饬鲌鰷囟仍茍D

圖9 地面及導(dǎo)流器導(dǎo)流面流場溫度云圖

從圖 8～圖 10可看出，考慮化學(xué)反應(yīng)后，在燃?xì)舛逊e區(qū)域、燃?xì)馀c空氣混合區(qū)域的射流溫度以及導(dǎo)流面、地面的溫度都顯著升高，高溫區(qū)域范圍明顯增大，這是由于在所述區(qū)域中發(fā)生了放熱的復(fù)燃現(xiàn)象；在燃?xì)鉀_刷的壁面上最高溫度位置無變化，出現(xiàn)在導(dǎo)流面的燃?xì)舛逊e處，最高溫度升高約200 K；地面上相同位置處的溫度差較大，并且在一定范圍內(nèi)隨著與噴管中軸的距離增大而增大，最大溫度差約為1250 K。

圖10 對稱面處燃?xì)鉀_刷壁面上的溫度分布

除了對比溫度與燃?xì)饨M分分布的差異，如圖 11、圖 12所示，對比了凍結(jié)流與化學(xué)反應(yīng)射流中導(dǎo)流面上的壓強(qiáng)分布以及對稱面處導(dǎo)流面、地面上的壓強(qiáng)分布。通過對比發(fā)現(xiàn)，在考慮化學(xué)反應(yīng)的計算結(jié)果中，導(dǎo)流面上的高壓區(qū)域分布與凍結(jié)流場相同，最大壓強(qiáng)升高約0.3 MPa，最大壓強(qiáng)的位置一致，地面上相同位置處的壓強(qiáng)明顯升高，壓強(qiáng)變化約為0.1 MPa。

圖11 地面及導(dǎo)流器導(dǎo)流面流場溫度云圖

圖12 地面及導(dǎo)流器導(dǎo)流面流場溫度云圖

4 結(jié)論

(1)復(fù)燃現(xiàn)象主要發(fā)生在射流與空氣的混合邊界層以及導(dǎo)流器和地面等高溫燃?xì)獾亩逊e處。

(2)化學(xué)反應(yīng)射流流場中燃?xì)鉀_刷的壁面上的最高溫度位置與凍結(jié)射流流場相同，為導(dǎo)流器的燃?xì)舛逊e處，最高溫度升高約200 K。

(3)加入化學(xué)反應(yīng)后，燃?xì)鉀_刷地面的溫度顯著升高，最多升高約1250 K，地面上的高溫區(qū)域明顯增大。

(4)加入化學(xué)反應(yīng)后，導(dǎo)流器表面的壓強(qiáng)升高，最大壓強(qiáng)出現(xiàn)的位置不變，最大壓強(qiáng)升高約0.3 MPa。