那佳琪， 楊 錄， 李文強， 張 明， 崔夢君

(中北大學(xué) 信息與通信工程學(xué)院 電子測試國家重點實驗室，山西 太原 030051)

0 引 言

固體火箭發(fā)動機測試是整個發(fā)動機實驗研究和工程發(fā)展階段的重要環(huán)節(jié)，也是推進(jìn)發(fā)動機技術(shù)的重要支撐。其中，尾焰溫度場的測試對火箭發(fā)動機來說至關(guān)重要。西安電子科技大學(xué)的趙文娟對固體火箭發(fā)動機尾焰流場進(jìn)行了研究分析，使用經(jīng)驗公式計算發(fā)動機尾焰流場，結(jié)果表明,含氧量越高，最小點火能量越低[1]。北京理工大學(xué)的王偉臣等人通過建立尾焰的計算模型研究了固體火箭發(fā)動機尾焰規(guī)律，使用歐拉和拉格朗日相結(jié)合的方法對尾焰流場進(jìn)行數(shù)值計算，結(jié)果表明不同點火劑量的增加對尾焰有著顯著影響[2]。綜上所述，目前,國內(nèi)外對于航空發(fā)動機燃燒室方面的技術(shù)研究,基本都是關(guān)注結(jié)構(gòu)和流動參數(shù)對尾焰流場的宏觀影響，而有關(guān)燃燒過程的動態(tài)模擬等方面比較少?；鸺l(fā)動機尾焰流場結(jié)構(gòu)十分復(fù)雜，其噴射出的氣流與環(huán)境中的氣體劇烈摻混，固體火箭發(fā)動機點火過程時間十分短暫，所以,化學(xué)試驗的方法無法對進(jìn)行精確分析，通過仿真的方法能夠?qū)Πl(fā)動機點火過程進(jìn)行數(shù)值研究。

本文針對固體火箭發(fā)動機尾焰溫度場,設(shè)計了一種基COMSOL Multiphysics軟件的尾焰流場數(shù)值仿真模型。開展發(fā)動機燃燒過程數(shù)值計算，研究成果對工程應(yīng)用具有重要的指導(dǎo)意義。

1 尾焰流場的數(shù)值模型

1.1 湍流方程

固體火箭發(fā)動機尾焰中的雷諾數(shù)Re大于2 300，因此,尾焰的流體處于湍流流體狀態(tài)[3,4]。

在求解流場中運用湍流模型時，選取了k-ω湍流模型

(1)

(2)

式中Gk為湍動能，Gω為ω方程，Yk和Yω為k和ω的發(fā)散項，τk和τω為k和ω的有效擴散項，Dω為正交發(fā)散項，Sk和Sω為用戶定義項[5]。

1.2 湍流反應(yīng)速率

使用湍流模型時，反應(yīng)產(chǎn)生的物質(zhì)產(chǎn)率模擬Rij定義為平均值閉合反應(yīng)速率和渦耗散模型速率的最小值[6]

Rij=VijMi×min[rMVC,jrED,j]

(3)

式中rMVC,j為平均值閉合反應(yīng)速率，rED,j為渦耗散模型速率。平均值閉合反應(yīng)速率是用平均質(zhì)量分?jǐn)?shù)表示的動力學(xué)反應(yīng)速率，對應(yīng)于反應(yīng)速率比湍流混合慢的反應(yīng)的特征反應(yīng)速率，或是湍流度可忽略不計的區(qū)域的反應(yīng)速率。該反應(yīng)速率可通過達(dá)姆科勒數(shù)Da來量化，達(dá)姆科勒數(shù)為湍流時間尺度與化學(xué)反應(yīng)時間尺度之比。平均值閉合適用于低達(dá)姆科勒數(shù)

(4)

渦耗散模型定義的反應(yīng)速率為

(5)

式中τT為湍流混合時間尺度，τc為化學(xué)反應(yīng)時間尺度，ρ為混合物密度，ω為質(zhì)量分?jǐn)?shù)，V為化學(xué)計量系數(shù)，M為摩爾質(zhì)量[7]。

當(dāng)Re和達(dá)姆科勒數(shù)都足夠高，導(dǎo)致反應(yīng)速率受到湍流混合時間尺度的限制。在分子層面，由于湍流的存在，全局反應(yīng)最多以新反應(yīng)物混合的速率進(jìn)行。當(dāng)反應(yīng)速率受到反應(yīng)物不足的限制時，反應(yīng)物具有最低局部濃度。模型參數(shù)指定反應(yīng)需要產(chǎn)物物質(zhì)，對活化能進(jìn)行建模。對于氣態(tài)非預(yù)混燃燒，取α=4，β=4。在當(dāng)前模型中，反應(yīng)的分子反應(yīng)速率無限快時，在模型中通過為反應(yīng)指定相當(dāng)高的反應(yīng)速率常數(shù)來實現(xiàn)。

1.3 點火藥質(zhì)量的計算

以點火壓力為依據(jù)，具體用量的計算公式為

(6)

式中ρ為火藥試樣密度，fb為點火藥的火藥力，pb為點火藥壓力，V為密閉爆發(fā)器燃燒室體積，mp為火藥試樣裝藥量，αb為點火藥氣體的余容。

2 尾焰溫度場有限元模型

2.1 數(shù)值計算方法

利用軟件求解二維Navier-Stokes方程組，使用基于密度算法的求解器，采用有限體積法對方程組進(jìn)行離散并通過二階迎風(fēng)格式進(jìn)行重構(gòu)，同時,對連續(xù)方程、動量方程、能量方程和組分運輸方程進(jìn)行求解。壁面附近采用標(biāo)準(zhǔn)壁面函數(shù)。

2.2 溫度場幾何模型與網(wǎng)格剖分

為方便計算，可將溫度場幾何模型采用二維軸對稱模型。以入口作為進(jìn)氣口，燃?xì)鈴娜肟谧杂缮淞鞫觯⒑雎宰陨淼暮穸?。具體尺寸和網(wǎng)格模型如圖1(a)所示。

圖1 溫度場幾何模型與網(wǎng)格剖分

為了在仿真過程中不影響網(wǎng)格密度，最好通過局部網(wǎng)絡(luò)加密以及結(jié)構(gòu)化和非結(jié)構(gòu)化網(wǎng)絡(luò)的組合來優(yōu)化網(wǎng)格分離。該模型對核心區(qū)域和火焰周圍區(qū)域的網(wǎng)格進(jìn)行了不同的剖分[8]。為進(jìn)行網(wǎng)格獨立性驗證，分別劃分為8，15，25萬 三組網(wǎng)格，并選取溫度、CO2質(zhì)量分?jǐn)?shù)作為主要求解參數(shù)進(jìn)行對比。結(jié)果表明，當(dāng)網(wǎng)格數(shù)量增加到一定程度時，溫度和CO2質(zhì)量分?jǐn)?shù)變化很小。為減小計算量，選用15萬的網(wǎng)格進(jìn)行火焰特征仿真計算。經(jīng)過不斷實驗，選取147 778萬的網(wǎng)格進(jìn)行仿真計算。尾焰的網(wǎng)格模型如圖1(b)所示。

模型網(wǎng)格劃分采用自由4面體網(wǎng)格，網(wǎng)格最大單元為0.039 9 m，最小網(wǎng)格尺寸為1.76×10-8m，時間步長的選擇受到穩(wěn)定性條件的限制，為了保證計算結(jié)果為收斂，選擇時間步長為(1.36×10-6)s。網(wǎng)格數(shù)量約147 778，為最佳解。由于仿真模擬是隨著時間的變化而變化，故研究選擇瞬態(tài)[9]。

2.3 溫度場邊界條件的設(shè)置

給定入口和軸向速度350 m/s；出口取外插值邊界條件。外插值的邊界條件，即所有變數(shù)值均使用外插計算[10]。壓力值為82 atm(l atm=101.325 Pa)，溫度為1 000 K。為計算發(fā)動機尾焰溫度場，首先，需要確定燃燒室內(nèi)燃燒產(chǎn)物及組分含量，利用基于吉布斯最小自由能方法的熱力計算可獲得相應(yīng)結(jié)果，本文使用文獻(xiàn)[10]所介紹的吉布斯最小自由能方法對入口射流組分的含量分布情況進(jìn)行熱力計算，為流場計算提供數(shù)據(jù)。入口射流含有6種化學(xué)成分，如表1所示。

表1 燃?xì)馍淞鞒煞?/p>

3 溫度場尾焰數(shù)值計算驗證與分析

3.1 數(shù)值模型計算尾焰溫度驗證

如圖2所示，實驗值和計算值溫度均在0.6 m左右的區(qū)域位置開始上升，并在1.0 m左右到達(dá)溫度最高值3 000 K。將計算結(jié)果與參考文獻(xiàn)[11]實驗數(shù)據(jù)進(jìn)行對比。實驗值與計算值整體趨勢基本一致，說明數(shù)值模型計算尾焰溫度驗證是可靠的。

圖2 溫度場尾焰溫度驗證

3.2 溫度場尾焰燃燒過程

點火過程與整個發(fā)動機工作過程相對比較短暫，為了得到更精確的分析結(jié)果將尾焰燃燒過程假設(shè)為4部分：誘導(dǎo)期、傳播期、充滿期和穩(wěn)定期，各期溫度云圖如圖3所示。

圖3 各期溫度云圖

在發(fā)出發(fā)動機點火指令后，在火箭發(fā)動機后迅速產(chǎn)生一個高溫區(qū)域，燃燒產(chǎn)物進(jìn)入高溫區(qū)域并充填自由容積，擠壓高溫區(qū)域原有的空氣，造成局部壓力開始升高從而引起振蕩，并產(chǎn)生燃?xì)?，如圖3(a)所示。由圖可知，在t=3.8 ms時，燃?xì)忾_始向外擴散，此時氣體內(nèi)能轉(zhuǎn)換成動能，以很高的速度向后噴射，從而獲得更大的速度，相應(yīng)溫度也就會較低，這個時候的溫度和周圍的環(huán)境一致。高溫區(qū)域的壓強還不穩(wěn)定，正在緩慢上升，擴大流場燃燒區(qū)域，溫度場存在振蕩。

剛產(chǎn)生燃?xì)鈺r，因溫度不高，沒有達(dá)到其臨界值。隨著燃燒反應(yīng)的進(jìn)行，產(chǎn)生的燃?xì)庠絹碓蕉?，高溫區(qū)域也逐漸增大，如圖3(b)所示。由圖可知，在t=16.1 ms時產(chǎn)生的高溫燃?xì)獬熟F狀向兩側(cè)迅速發(fā)展，火焰繼續(xù)迅速發(fā)展，因此，傳播期在整個點火過程中時間占比非常大。

當(dāng)用火點燃的藥劑開始釋放大量高溫氣體時，火焰開始蔓延到藥劑的末端，此時開始燃燒，然后進(jìn)入充滿期，如圖3(c)所示。由圖可知，在t=156 ms時，此時已經(jīng)全部點火成功，但是尾焰核心區(qū)域溫度仍然沒有達(dá)到平衡狀態(tài)，燃?xì)獬掷m(xù)不斷的向后半段推進(jìn)，同時前半段核心區(qū)域溫度迅速上升。

隨著尾焰釋放的燃?xì)庠蕉啵l(fā)動機高溫區(qū)域充滿期的時間越短，發(fā)動機能夠以更快的時間達(dá)到工作狀態(tài)，隨著充滿期深入的進(jìn)行，發(fā)動機內(nèi)部流場將更趨于穩(wěn)定，如圖3(d)所示。由圖可知，發(fā)動機已經(jīng)進(jìn)入穩(wěn)定燃燒狀態(tài)，此時溫度分布處于穩(wěn)定階段，發(fā)動機現(xiàn)已完全燃燒，釋放了大量燃燒氣體，火箭將繼續(xù)運行，氣體到達(dá)尾部，釋放的燃?xì)廨^少。

本文模型的傳播時間非常短，從計算結(jié)果上看，大約196.7 ms后，固體火箭發(fā)動機進(jìn)入穩(wěn)定燃燒階段，高溫區(qū)域內(nèi)的溫度趨于穩(wěn)定，最高溫度能夠達(dá)到3 000 K。

3.3 不同質(zhì)量點火藥對溫度場影響

對點火藥質(zhì)量為2，5，8 g三種狀態(tài)對發(fā)動機尾焰的溫度影響特性的分析，研究了點火藥性質(zhì)對溫度場特征的影響作用，如圖4所示?？梢钥闯觯粫r刻2 g點火劑量產(chǎn)生的高溫區(qū)域最小，而8 g最大。

圖4 同一時刻不同點火劑量溫度云圖

隨著點火藥量的增加，產(chǎn)生的高溫燃燒氣體越多，尾焰核心區(qū)及中后部高溫區(qū)的長度和寬度均有所增加，如圖5所示。可以看出，隨著點火藥劑量的增加，溫度場最高溫度不斷的在增加。因此，在保持藥柱結(jié)構(gòu)完整性的同時，火藥劑量增加，從而釋放更多燃燒氣體，同時產(chǎn)生高溫區(qū)域的面積較大，發(fā)動機燃燒充滿期變短，發(fā)動機運行速度越快。

圖5 同一時刻不同點火劑量溫度云圖

由此可以推出，燃燒室內(nèi)溫度隨著點燃劑量的增加而逐漸升高，如果點燃劑量過大并超過其極限，導(dǎo)致藥柱破碎，影響其燃燒效率，使發(fā)動機發(fā)生故障，甚至可能導(dǎo)致發(fā)動機爆炸；而點火劑量太小，會出現(xiàn)缺點、點火過度延遲和間歇燃燒等現(xiàn)象。

3.4 不同壓強對溫度場的影響

對工作環(huán)境壓強等級各設(shè)為6，8，10，12，14 MPa，如圖6所示。通過比較可發(fā)現(xiàn)，由于工作壓強的上升，在尾焰軸上溫度峰值的位置逐漸向后方移動，這是由于燃燒室工作壓強的上升使得傳播到周圍空氣的范圍逐漸擴大，因此導(dǎo)致噴嘴出口處附件的密度降低更明顯，該處的溫度也降低一定程度。

圖6 同一時刻不同壓強溫度云圖

4 結(jié) 論

本文建立了固體火箭發(fā)動機尾焰模型，對燃燒過程和不同條件下溫度場特性進(jìn)行研究，得到溫度場的計算結(jié)果。仿真結(jié)果表明，在保證藥柱結(jié)構(gòu)完整性的前提下，點火劑量越大，工作壓強越大，產(chǎn)生的高溫區(qū)域面積越大，從而縮短點火時間，加快固體火箭的反應(yīng)速度。