周 康，李清廉，成 鵬，常一冰

(國防科技大學,湖南 長沙 410073)

0 引言

相對傳統(tǒng)噴注器，針栓噴注器結構簡單、可靠性高且燃燒較為穩(wěn)定，在應用上具有獨特的優(yōu)勢。1961年，首個針栓噴注器成功應用于MIRA 500火箭發(fā)動機上，推力變化為111.2～ 2 224 N。此后，針栓噴注器又先后應用到阿波羅登月下降級發(fā)動機、姿軌控發(fā)動機TR 312以及Space X的梅林1 D發(fā)動機[1-2]上，足見其應用十分廣泛。

目前針栓噴注器霧化研究大部分集中于對噴霧錐角的研究。對于徑向環(huán)縫式針栓噴嘴，Son等[3]發(fā)現噴霧錐角隨徑向/軸向動量比的增大而增大，并總結出了與氣液動量比TMR和We(韋伯數)有關的氣液針栓噴嘴半錐角的指數關系式，數值仿真工作[4]也得到了同樣的規(guī)律。方昕昕等[5]發(fā)現氣液比的增加會減小噴霧錐角，徑向環(huán)縫夾角、無量綱跳躍距離、液氧環(huán)縫厚度的增加均會增加噴霧錐角。成鵬[6]也研究了液液針栓噴嘴的噴霧錐角隨工況的變化，并通過理論分析提出了一個與徑向/軸向動量比TMR有關的半噴霧錐角的預測公式。數值模擬的結果還表明軸向液膜與徑向液膜相互撞擊在徑向環(huán)縫出口位置形成局部高壓區(qū)，該局部高壓區(qū)促進了噴霧錐角增大。

但是國內外在針栓噴注器的燃燒性能與燃燒過程基礎研究方面還比較有限。針栓發(fā)動機還存在燃燒效率較低的問題。國外普渡大學針對667.5 N的針栓小推力發(fā)動機，研究了總動量比、孔徑比以及針栓長度與針栓直徑之比對燃燒性能的影響，但并沒有得出燃燒性能與這些參數之間明確的相互作用規(guī)律[7-8]。2016年日本東京大學的Sakaki等[9-10]以LOX/Ethanol為推進劑，進行了平面針栓(開窗觀測)和軸對稱針栓的對比試驗，發(fā)現兩種情況下推力室燃燒效率均對TMR呈現反相關的關系，故而Sakaki認為可以用便于使用光學觀測等手段的平面針栓來進行觀察針栓燃燒過程的基礎研究。仿真方面，2017年國防科技大學方昕昕的進一步研究表明：跳躍距離Ls與針栓直徑之比為1時，特征速度效率最高(約0.96)。方昕昕[5]和韓國航天大學的Min Son等[11-12]研究表明了徑向環(huán)縫尺寸影響燃燒效率。方昕昕[5]的研究還表明燃燒室特征長度越長燃燒效率越高，這一點與常規(guī)噴注器發(fā)動機的結論是一致的?？偟膩碚f，對燃燒效率的基礎研究還不夠深入，需進一步明確結構工況參數對燃燒效率的影響規(guī)律。

在針栓發(fā)動機的燃燒過程基礎性研究上，仿真方面國內學者以及韓國學者一般采用Lagrange隨機軌道模型來模擬推進劑液滴的運動軌跡[13-16]，并以此作為霧化、蒸發(fā)乃至最終燃燒計算的初始條件。西北工業(yè)大學的李進賢[15]及北京航空航天大學的張連博等[14]對NTO/MMH針栓發(fā)動機燃燒室流場結構和推力室性能的影響因素進行了2D仿真研究，結果表明針栓深入燃燒室的長度對內流場結構有一定影響。2017年國防科技大學方昕昕的研究表明：跳躍距離過小時，徑向噴出的液氧射流將在甲烷氣流的撞擊下加速破碎，并與燃燒室壁面發(fā)生碰撞并反彈，然后向燃燒室中心移動，造成回彈位置附近和燃燒室中心富氧而燃燒不充分，出現相對低溫區(qū)。徑向射流環(huán)縫大小控制著液氧流速大小，開口愈大，則液氧流速愈小，對維持燃燒穩(wěn)定性起著非常重要的針栓尖部下方的回流區(qū)愈小，使得其對推進劑霧化混合的促進作用減弱。Min Son等[11-12]的研究也給出了類似結論，研究還表明相同結構下，噴注質量流率減小，火焰張角隨之減小，針栓尖部下方回流區(qū)變小。總而言之，現有資料缺乏燃燒流場深入的研究，應從基本流場分布和火焰結構方面來解釋結構工況參數對燃燒過程的作用規(guī)律。

縱觀國內外關于針栓發(fā)動機的燃燒性能和燃燒過程基礎研究還是比較有限的。跳躍距離和動量比等對燃燒過程包括流場結構、溫度分布和火焰結構等的影響機理缺乏更深層次的研究，因此還需要繼續(xù)開展仿真研究來明確燃燒性能和燃燒過程與針栓噴注器關鍵設計參數之間的規(guī)律。

1 物理模型與計算方法

1.1 物理模型與研究問題

本次研究的對象如圖1所示的二維發(fā)動機結構，發(fā)動機總流量為152 g/s，其中甲烷流量為36.2 g/s。發(fā)動機喉部半徑為9.55 mm，圓筒段半徑為22.5 mm，圓筒段長度為48.6 mm。In1(甲烷噴注口)噴進燃燒室的是甲烷，In2(氧氣噴注口)噴進燃燒室的是氧氣，箭頭方向表示燃料和氧化劑的流動方向。工作過程中，兩股氣流在燃燒室內混合燃燒，發(fā)生化學反應，產生高溫高壓氣體，再經過噴管加速之后，以高速氣流排出，在這次仿真中暫時不考慮推力室壁面?zhèn)鳠釋攘鲌龊蜏囟葓龅挠绊憽?/p>

研究采用商業(yè)軟件FLUENT進行數值仿真，研究的問題在于氣氧氣甲烷針栓發(fā)動機的燃燒特性，考察推進劑發(fā)生化學反應后燃燒室的內流場結構和溫度場分布，研究參數主要有不同針栓噴孔直徑、位置以及環(huán)縫寬度，研究的孔的位置有兩個，分別為位置A和B，如圖1所示。對應仿真算例如表1所示，對比算例1和2，可分析跳躍距離對內流場結構和溫度場分布的影響，對比算例1，3，4和5可分析動量比對內流場結構和溫度場分布的影響。

圖1 仿真模型的結構Fig. 1 Configuration of simulation model

算例序號孔的位置孔徑/mm環(huán)縫寬度/mm跳躍距離/mm徑向與軸向動量比備注Case 1B位置0.51.311.40.66Case 2A位置0.51.370.66動量比不變Case 3B位置1.11.311.40.47變徑向動量Case 4B位置0.51.511.40.81Case 5B位置0.5211.41.15變軸向動量

1.2 控制方程

處理含有化學反應的傳熱和流動問題需要求解質量、動量、能量守恒和化學組分輸運的聯立方程組，為更清楚地表達研究所涉及的控制方程及其物理意義，這里采用計算機程序計算時所采用的通用方程[17]：

(1)

式中：ρ為密度；v為速度向量;φ為待求變量;г為變量φ的擴散系數;S為源項。由該式可看出通用微分方程一般包含4項，即非穩(wěn)態(tài)項、對流項、擴散項和源項，本問題研究燃燒室內定常流場，所以非穩(wěn)態(tài)項直接忽略。

動量方程包括x和y兩個方向的方程，源項中的兩項分別代表單位質量流體的壓力差和體積力。

對于化學反應，這里采用最簡單的一步總包反應：

CH4+2O2→CO2+2H2O

(2)

反應共有4種組分，由于質量分數之和為1，故只需求解3個組分輸運方程，能量方程是以溫度T為變量來表達。

另由于流動是可壓縮流動，故需要補充狀態(tài)方程，假定流體滿足理想氣體狀態(tài)方程：

p=ρRT

(3)

式中：T為溫度；R為氣體常數。

1.3 計算模型

這里主要說明湍流模型、化學反應模型以及湍流和化學反應相互作用的模型。

湍流模型采用Standardk-εmodel，其可適用于航空發(fā)動機和火箭發(fā)動機燃燒室內的流動計算[18]。

化學反應模型選擇組分輸運模型(Species Transport)。本問題研究的是氣氣燃燒，所以選擇了組分輸運模型，該模型應用廣泛，能夠滿足本問題的仿真要求[19]，且計算成本也較小。

對于湍流和化學反應的相互作用模型，選用渦耗散模型(Eddy-Dissipation Model)，只要湍流出現(k/ε>0)，燃燒即可進行，不需要點火源來啟動燃燒。這通常對于非預混燃燒是可行的。在渦耗散模型中，每個反應都有同樣的湍流速率，因而模型能用于單步或是雙步整體反應，并適用于高雷諾數燃燒過程[20]。該研究中采用一步總包反應，且經計算本問題的雷諾數達105量級，所以綜合來看選擇渦耗散模型是較為合適的。

1.4 邊界條件與物性參數

仿真模型邊界條件設置如表2所示。

表2 邊界條件

1.5 算法使用

綜合考慮各種選項的精度、計算開銷、數值穩(wěn)定性等，本文采用SIMPLE算法，空間離散上除湍動能和湍流耗散率采用一階迎風各式外，其余各變量均采用二階迎風各式，壓力采用二階格式，梯度采用最小二乘法計算。在本次研究中，氧氣與甲烷的混合比為3.2，在該混合比下，根據蔡震宇[21]研究的不同室壓下氣氧/氣甲烷理論燃燒溫度，最高反應溫度不超過3 470 K，超出的溫度無研究意義，因此考慮算法中加入限制器，并且研究關注的是溫度分布規(guī)律，而不是最高溫度，這是合理化的處理。

2 網格無關性及計算方法驗證

基準算例的網格無關性驗證采用3套網格，網格數量級分別為6萬、10萬和20萬。仿真主要關注燃燒室內的溫度場和內流場，因而網格無關性驗證的參數選擇溫度和壓力。壓力和溫度的計算都是基于Area-weighted Average算法求值，原理為將所有網格面上的物理量φi與對應網格面的面積Ai乘積之和除以總面積A，即可求的對應物理量的平均值，公式如式(4)所示：

(4)

在軸向上選取某點處(該次選取x=63 mm)徑向直線，導出數據進行網格無關性驗證，結果如圖2所示。對于不同數量網格計算的最大壓力誤差為1.04%，在溫度方面，粗網格與密網格誤差為3.01%，中度密網格與密網格基本符合，最大誤差僅為0.85%，因此綜合考慮精確度與計算成本使用10萬數量級的中度密網格。

圖2 網格無關性驗證結果Fig. 2 Result of grid independent verification

對于計算方法的驗證，采用高玉閃[22]的同軸剪切式氣氧/氣甲烷單噴嘴物理模型，計算模型和本研究的保持一致，邊界條件設置和試驗參數設置一致。由于研究關注的燃燒特性與溫度和壓力相關，而試驗結果給出的是壓力結果，所以此次模型驗證指標為室壓。通過仿真得出的室壓結果為2.75 MPa，試驗結果為2.56 MPa，誤差為7.4%。論文關注溫度分布規(guī)律的差別，在誤差允許范圍內(小于10%)，用此計算方法來研究氣氧/氣甲烷針栓發(fā)動機的燃燒特性可行的。

3 結果討論與分析

3.1 基準算例結果與分析

以甲烷噴孔位于B位置，噴孔直徑d=0.5 mm為基準算例，對其仿真結果進行分析，分析壓力和組分分布云圖。

此前，定義初始燃燒反應帶狀區(qū)域氧氣濃度梯度為0.001(即氧氣濃度幾乎處于不變，如公式(5)所示)時的界面線與軸線方向形成的銳角為火焰錐角。燃燒效率使用特征速度效率表征。

(5)

燃燒仿真得出的基準算例燃燒室壓力分布如圖3所示，其平均室壓為0.922 MPa。甲烷和氧氣的組分分布云圖如圖4所示，甲烷和氧氣的混合主要在于兩組分的分界面之間，形成了一個混合層，化學反應也主要發(fā)生在這個混和層內。

圖3 壓力分布云圖Fig. 3 Contour of pressure distribution

圖4 組分分布云圖Fig. 4 Contour of species distribution

圖5的流線分布圖展示了針栓噴注器存在兩組明顯的回流區(qū)，并且角回流區(qū)和針栓頭部的回流區(qū)會影響其燃燒特性[8]。

聯創(chuàng)汽車電子有限公司（DIAS）由上海汽車集團股份有限公司和中聯汽車電子有限公司聯合投資，主要從事汽車電子產品的研發(fā)、制造、銷售和技術服務，致力于成為國內領先、國際一流的汽車電子系統(tǒng)制造商；DIAS發(fā)展迅速，已成為國內主要汽車制造公司（如上汽、上海通用和奇瑞）的供應商。

圖5 流線分布云圖Fig. 5 Contour of streamline distribution

定量分析溫度場(如圖6所示)的變化，在軸向方向上，截取x=50 mm和x=100 mm兩條徑向直線上溫度變化(如圖7所示)，溫度隨Y向距離增加呈現出先增加后減少的趨勢。這是因為反應主要集中在甲烷與氧氣混合的帶狀區(qū)，因此反應釋熱的高溫區(qū)在氧化劑和燃料之間混合的帶狀區(qū)?；鶞仕憷幕鹧驽F角為54°，如圖6所示，燃燒效率的結果為87.8%，平均排氣速度為2 122 m/s。

圖6 溫度分布云圖Fig. 6 Contour of temperature distribution

圖7 不同位置溫度隨徑向距離分布Fig. 7 Distribution of temperature with different radial distances

3.2 針栓噴孔位置對溫度場分布的影響

對針栓噴孔處于不同位置算例進行燃燒仿真，得到圖8中結果。定義溫度處于3 400 K與3 500 K之間的帶狀區(qū)域為高溫區(qū)，對應的溫度變化帶狀區(qū)域的寬度為帶寬。對溫度仿真結果進行定量分析，在Case 1和Case 2的溫度云圖中，取x= 50 mm如和x=100 mm徑向直線上的溫度進行分析，其結果見圖9所示。A位置(Case 2)燃燒室?guī)罡邷貐^(qū)域要比B位置(Case 1)的寬，但火焰錐角沒有隨著噴孔位置改變而改變，通過計算Case 1燃燒室內平均溫度為2 179.5 K，Case 2為2 398.5 K。

圖8 不同孔位置的溫度云圖Fig. 8 Contour of temperature with different positions of pintle orifice

圖9 不同孔位置的溫度隨徑向距離分布Fig. 9 Distribution of temperature with different positions of pintle orifice

對于以上結果，分析原因在于不同的噴注位置(Case 1和Case 2)，實質上改變了跳躍距離，Case 1的跳躍距離為11.4 mm，Case 2的跳躍距離為7 mm。跳躍距離減小，導致Case 2中氧氣與甲烷混合界面帶擴大，燃燒反應釋熱帶相對Case 1寬，穩(wěn)態(tài)狀態(tài)下，靠近壁面的溫度和對稱軸的溫度比Case 1溫度高。同時跳躍距離減小，使得氧氣甲烷在燃燒室內滯留時間比Case 1長，等效增加了特征長度，混合燃燒更加充分，導致Case 2 的平均溫度要高于Case 1，室壓同樣如此，Case 2平均室壓為0.97 MPa，Case 1為0.922 MPa，這表明Case 2的燃燒效率高于Case 1，即跳躍距離減小，燃燒效率增大，為92.3%，平均排氣速度為2 235 m/s。但是跳躍距離沒有改變火焰錐角。對兩個算例燃燒結果總結如表3所示。

表3 Case 1與Case 2結果總結

3.3 針栓噴孔直徑對溫度場分布的影響

不同針栓噴孔直徑的溫度和混合比仿真結果如圖10和圖11所示，從溫度云圖可看出，火焰錐角隨噴孔直徑增加而減小。對溫度進行分析，得到如圖12所示變化曲線，從曲線中可知，Case 1燃燒室?guī)罡邷貐^(qū)與Case 3的帶狀高溫區(qū)寬度接近，接近壁面處Case 3溫度高于Case 1，而靠近對稱軸附近剛好相反。并通過計算得出Case 3平均溫度為2 236 K，此時Case 3平均室壓為0.925 MPa，與Case 1相近。

對此進行分析，噴孔直徑增大，徑向動量變小，使得動量比減小，Case 1的動量比為0.66，Case 3為0.47。動量比減小，火焰錐角減小，如圖 10所示，所以實質上火焰錐角隨動量比減小而變小。仿真混合比為3.2，動量比減小，徑向氣流穿透軸向氣流能力減弱，相對的軸向氣流增強，擴散到壁面處氧氣增加，從圖11可知，3和4 之間的局部混合比區(qū)域面積Case 3大于Case 1，即動量比減小導致靠近壁面處氧氣與甲烷混合反應更加充分，而總體性能燃燒效率與Case 1接近，結果為88.1%，平均排氣速度為2 152 m/s。兩者燃燒性能總結如表4所示。

圖10 不同孔直徑的溫度云圖Fig. 10 Contour of temperature with different pintle orifice diameters

此外，從圖11可看出在帶狀高溫區(qū)域下方的混合情況很差，這是因為在二維情況下，氧氣穿透甲烷氣膜進入到帶狀高溫區(qū)域下方的量很少，導致混合變差。改善的方法在于減小徑向氣膜噴射速度，減小徑向動量，增加軸向氣膜的穿透能力，使得進入到帶狀高溫區(qū)域下方氧氣的量增加，提高氧氣與甲烷的摻混，從而改善其燃燒性能，Case 1和Case 3的仿真結果驗證了該方法的可行性。

圖11 不同孔直徑的混合比Fig. 11 Contour of mixture ratio with different pintle orifice diameters

算例燃燒效率/%排氣速度/(m·s-1)火焰錐角/(°)Case 187.82 12254Case 388.12 15235

圖12 不同位置不同孔徑溫度隨徑向距離分布Fig. 12 Distribution of temperature with different positions of pintle orifice diameter

3.4 針栓環(huán)縫寬度對溫度場分布的影響

圖13 不同環(huán)縫寬度的溫度分布Fig. 13 Contour of temperature with different annular gaps

圖14 不同位置不同環(huán)縫寬度溫度隨徑向距離分布Fig. 14 Distribution of temperature with different annular gaps

圖15 不同環(huán)縫寬度流線圖Fig. 15 Streamline with different annular gaps

分析以上結果，環(huán)縫寬度變大，使得軸向動量減小，動量比增加，混合帶狀區(qū)域向燃燒室頭部傾斜，導致燃燒時Case 4和Case 5火焰錐角變大。動量比增加，徑向氣流撞擊壁面，使得擴散到靠近壁面處氧氣減少，氧氣與甲烷混合變差，反應釋熱減少，溫度降低，如圖14曲線分布，并且動量比增加，靠近針栓頭部的回流區(qū)軸向尺寸減小，如圖15所示，對靠近喉部氣體卷吸作用減弱，氣體與高溫區(qū)交換熱量減小，軸向溫度降低。環(huán)縫寬度增加，導致動量比增加，氧氣與甲烷混合燃燒變差，使得燃燒室內平均溫度降低，平均室壓降低，燃燒效率降低，Case 4和Case 5的燃燒效率分別為86.3%和86.1%，平均排氣速度分別為2 055 m/s和2 036 m/s。以上燃燒結果可總結如表5所示。

表5 Case 1，Case 4 和Case 5結果總結

3.5 火焰錐角對燃燒效率的影響

在這次仿真中發(fā)現噴孔直徑、位置以及環(huán)縫寬度均會對燃燒效率產生影響，而其中火焰錐角也會發(fā)生相應變化，因此燃燒效率與火焰錐角其中規(guī)律可以通過曲線圖來說明，如圖16所示。燃燒效率隨著火焰錐角增加而減小。

圖16 不同火焰錐角下的燃燒效率Fig. 16 Combustion efficiency with different flame angles

對此現象可以從混合角度解釋，本文理論混合比為3.2，因此局部混合比處于3到4之間可認為是理想混合范圍，因而可從3.3節(jié)中圖11可以簡化建立一個理想混合區(qū)模型，如圖17所示。理想混合主要有兩個，一個靠近燃燒室頭部，另一個靠近燃燒室壁。在前面分析中可知，火焰錐角的變化是由動量比變化引起的，動量比增加，火焰錐角增加，而此時徑向射流由于動量相對軸向增加，使得射流撞擊壁面，而使得靠近壁面氧氣減小，而甲烷濃度沒有變化，因而混合變差，靠近壁面的理想混合區(qū)域被壓縮變小，燃燒不充分，燃燒效率也會相應降低。因而提高燃燒效率可以通過減小火焰錐角實現，但火焰錐角不能無限減小，推測存在一個變化范圍，這需要通過進一步研究來確定。

圖17 理想混合區(qū)模型Fig. 17 Model of ideal mixing zone

4 結論

仿真主要模擬了二維氣氧/氣甲烷針栓發(fā)動機的內流場結構和溫度場分布，研究了噴孔位置、直徑和環(huán)縫寬度對流場結構和溫度分布的影響，主要結論有；

1)穩(wěn)態(tài)下的二維氣氧氣甲烷燃燒仿真，溫度場在燃燒室中心一定區(qū)域內呈帶狀分布。

2)噴孔位置位于A位置，跳躍距離減小，反應釋熱區(qū)域變寬，燃燒室內平均溫度增加，室壓增加至0.97 MPa，燃燒效率提高至92.3%，而火焰錐角不變。

3)環(huán)縫寬度增加，動量比增加，火焰錐角增加，燃燒室內平均溫度減小，室壓從0.922 MPa降低到0.905 MPa，燃燒效率從87.8%降低至86.1%。

4)噴孔直徑增加，動量比減小，火焰錐角減小，燃燒效率變化很小，從87.8%增加到88.1%。

5)在35°到68°的錐角范圍之內，燃燒效率隨著火焰錐角增加而減小。