俞南嘉，鮑啟林，張 洋，戴 健

(北京航空航天大學(xué) 宇航學(xué)院，北京 100191)

0 引言

目前，液體火箭發(fā)動(dòng)機(jī)是航天運(yùn)輸及空間飛行器的主要?jiǎng)恿ρb置，液氧/煤油火箭發(fā)動(dòng)機(jī)因其制造成本低、發(fā)動(dòng)機(jī)性能高、環(huán)保無污染、推力可控以及安全系數(shù)高等特點(diǎn)，在航天運(yùn)輸領(lǐng)域得到廣泛應(yīng)用[1]；同時(shí)，液氧/煤油火箭發(fā)動(dòng)機(jī)也是未來可重復(fù)使用發(fā)動(dòng)機(jī)發(fā)展的重要方向，可重復(fù)使用火箭發(fā)動(dòng)機(jī)意味著火箭的子級(jí)可以回收并且重復(fù)使用，其實(shí)現(xiàn)的重要途徑就是研發(fā)推力可調(diào)的火箭發(fā)動(dòng)機(jī)[2]。

針?biāo)ㄊ絿娮⑵髯鳛樽兺屏Πl(fā)動(dòng)機(jī)應(yīng)用的典型，被廣泛應(yīng)用到火箭發(fā)動(dòng)機(jī)中。最著名的應(yīng)用就是美國TRW公司登月艙下降段發(fā)動(dòng)機(jī)，該發(fā)動(dòng)機(jī)先后將12名宇航員送上月球[3]；1995年美國又研制了世界上最大的針?biāo)òl(fā)動(dòng)機(jī)TR-106，推力達(dá)到了290噸級(jí)[4]；美國SpaceX公司于2015年12月完成了獵鷹9號(hào)火箭的首次陸上垂直降落回收，其中，獵鷹9號(hào)火箭的梅林1D發(fā)動(dòng)機(jī)即采用針?biāo)ㄊ絿娮⑵鱗5]，推力可實(shí)現(xiàn)50 %～100 %的調(diào)節(jié)；2016年，Min Son等人提出了針?biāo)ㄊ絿娮⑵髯钚碌睦碚撛O(shè)計(jì)方法[6]；我國于上世紀(jì)70年代開始針?biāo)ㄊ絿娮⑵鞯难芯浚?992年研制成功機(jī)械定位的變推力發(fā)動(dòng)機(jī)，該發(fā)動(dòng)機(jī)采用自燃推進(jìn)劑N2O4/UDMH[6]；2013年嫦娥三號(hào)探測(cè)器采用雙調(diào)針?biāo)ㄊ絿娮⑵鲗?shí)現(xiàn)了月球軟著陸，該發(fā)動(dòng)機(jī)依然采用常溫自燃推進(jìn)劑[7]。但是，目前我國對(duì)于低溫非自燃推進(jìn)劑的針?biāo)ㄊ桨l(fā)動(dòng)機(jī)的研究很少，開展針?biāo)ㄊ揭貉?煤油發(fā)動(dòng)機(jī)的研究很有意義。

1 計(jì)算模型

1.1 物理模型

針?biāo)ㄊ阶兺屏σ后w火箭發(fā)動(dòng)機(jī)的原理是通過調(diào)節(jié)噴注器中的環(huán)形套筒來調(diào)節(jié)兩個(gè)同軸環(huán)縫的噴注通道的流通面積，實(shí)現(xiàn)對(duì)燃料和氧化劑噴注壓降、噴注液膜厚度等參數(shù)的調(diào)節(jié)，其中，環(huán)形調(diào)節(jié)套筒稱為針?biāo)ㄊ絿娮⑵鞯尼橀y，兩個(gè)環(huán)形噴注通道稱為氧化劑噴嘴和燃料噴嘴[8]。針?biāo)ㄊ絿娮⑵髟韴D如圖 1所示。

圖 1 針?biāo)ㄊ絿娮⑵髟韴D

1.2 數(shù)學(xué)模型

液體火箭發(fā)動(dòng)機(jī)的噴注燃燒一般是氣液兩相的噴注燃燒，推進(jìn)劑通過氧化劑噴嘴和燃料噴嘴以液體射流或者液膜的形態(tài)噴入燃燒室，氧化劑和燃料之間強(qiáng)烈的相互作用，使得液體推進(jìn)劑發(fā)生碰撞、霧化、混合等過程，并在對(duì)流、輻射等傳熱過程中蒸發(fā)和燃燒。因此，數(shù)學(xué)模型涉及到控制方程、湍流模型、燃燒模型等[9]。

1.2.1 控制方程

液體火箭發(fā)動(dòng)機(jī)噴霧燃燒的控制方程包含氣相控制方程和液相控制方程，氣相控制方程一般用歐拉坐標(biāo)系中的N-S方程來描述，液相控制方程則是由拉格朗日坐標(biāo)系中N-S方程對(duì)單個(gè)液滴進(jìn)行描述，氣液兩相之間的耦合采用氣液兩相的耦合源項(xiàng)來考慮[8]。

對(duì)于有N個(gè)組元的混合物來說，控制方程由連續(xù)方程、動(dòng)量方程、能量方程、組元方程以及狀態(tài)方程組成封閉方程組，其中，連續(xù)方程為

(1)

動(dòng)量方程為

(2)

能量方程為

(3)

組元方程為

(4)

狀態(tài)方程為

p=ρRT

式中：ρ為混合物密度；Sm為質(zhì)量守恒方程中的源項(xiàng)，代表單位體積中所有液滴蒸發(fā)導(dǎo)致的氣相質(zhì)量增長率；Sdi,Sdm為動(dòng)量方程中相互作用的源項(xiàng)，分別為單位體積的液滴反作用力使氣相的動(dòng)量的增長率和單位體積液滴蒸發(fā)帶給氣相的動(dòng)量；Se為能量守恒方程中的相互作用源項(xiàng)；Sg為組分方程中的源項(xiàng)；h為混合物靜焓，τij為應(yīng)力張量，μ為動(dòng)力粘性系數(shù)，wi為組分i的化學(xué)反應(yīng)速度，R為通用氣體常數(shù)。

1.2.2 湍流模型

湍流模型采用k-ε雙方程模型來計(jì)算液體火箭發(fā)動(dòng)機(jī)的燃燒，該模型廣泛應(yīng)用于各種工程流動(dòng)和熱交換的數(shù)值模擬中，并得到了較好工程驗(yàn)證，k-ε模型方程為：

(5)

(6)

式中：C1ε和C2ε為經(jīng)驗(yàn)常數(shù)，C1ε=1.44，C2ε=1.92；σk和σε分別是與湍流動(dòng)能k和耗散率ε對(duì)應(yīng)的普讓特?cái)?shù)；模型常數(shù)取值：Cμ=0.09，σk=1.0，σε=1.3；μt為湍流粘性；Gk是由于速度梯度引起的湍流動(dòng)能k的產(chǎn)生項(xiàng)。

1.2.3 組分輸運(yùn)模型

燃燒的組分輸運(yùn)模型采用有限速率-渦耗散模型，反應(yīng)考慮了組分之間的相互擴(kuò)散過程和推進(jìn)劑進(jìn)口的擴(kuò)散作用；該模型的化學(xué)反應(yīng)速率由大渦混合時(shí)間尺度k/ε控制，并對(duì)Arrhenius和渦耗散反應(yīng)速率都進(jìn)行計(jì)算，并取兩者之間化學(xué)反應(yīng)速率較小的一個(gè)[9]；液氧和煤油的化學(xué)反應(yīng)機(jī)理采用單步總包反應(yīng)，其化學(xué)反應(yīng)關(guān)系式為：

C12H23+17.75O2=12CO2+11.5H2O

由阿倫尼烏斯計(jì)算的燃燒反應(yīng)速率為：

(7)

由渦耗散控制的化學(xué)反應(yīng)速率為：

(8)

(9)

2 設(shè)計(jì)工況

所設(shè)計(jì)發(fā)動(dòng)機(jī)的總質(zhì)量流量為381 g/s，理論比沖為2 806.59 m/s，特征速度為1 776.6 m/s，噴注方案為燃料中心式，即燃料在內(nèi)側(cè)噴注，噴注角度為90o；氧化劑在外側(cè)噴注，噴注角度為52.36o，發(fā)動(dòng)機(jī)的參數(shù)如表 1所示。

2.1 網(wǎng)格劃分

由于發(fā)動(dòng)機(jī)計(jì)算域結(jié)構(gòu)的對(duì)稱性，本文采用1/8模型進(jìn)行網(wǎng)格的劃分，使用ICEM軟件劃分計(jì)算域網(wǎng)格。由于燃燒室內(nèi)存在劇烈的化學(xué)反應(yīng)，并伴隨著復(fù)雜的流動(dòng)過程，因此為更準(zhǔn)確地捕捉流場信息，對(duì)噴嘴、燃燒室頭部和燃燒室壁面進(jìn)行網(wǎng)格加密，計(jì)算域網(wǎng)格數(shù)量為118萬，圖 2為燃燒室頭部的計(jì)算域網(wǎng)格。

圖2 燃燒室頭部計(jì)算域網(wǎng)格

2.2 邊界條件設(shè)置

1)入口邊界條件：對(duì)于液體火箭發(fā)動(dòng)機(jī)的數(shù)值仿真，一般采用離散相模型作為入口邊界。本文設(shè)置液滴每隔20步化學(xué)反應(yīng)噴注一次，同時(shí)考慮顆粒對(duì)連續(xù)相的影響，顆粒軌道的迭代步數(shù)設(shè)置為1 000，離散物模型選擇壓力梯度模型、隨機(jī)碰撞模型以及裂解模型；液氧和煤油顆粒噴射源的噴注類型選擇SURFACE，顆粒形式為DROPLET，同時(shí)分別設(shè)置液滴的入口溫度、入口速度、質(zhì)量流量、液滴尺寸等，其中，煤油和液氧的液滴尺寸均設(shè)置為0.4 mm。

2)出口邊界條件：發(fā)動(dòng)機(jī)噴管出口一般為超音速的氣流，各參數(shù)計(jì)算采用外插公式得到。出口邊界條件設(shè)置湍流粘度和水利直徑，出口壓力為環(huán)境壓力。

3)壁面邊界條件：壁面邊界采用絕熱無滑移的壁面條件，即u=w=v=0；并且壓力梯度導(dǎo)數(shù)、組分質(zhì)量分?jǐn)?shù)梯度導(dǎo)數(shù)為零。

4)對(duì)稱面邊界：流場是對(duì)稱分布的，因此對(duì)稱面采用對(duì)稱面邊界條件，對(duì)稱面處參數(shù)的法向?qū)?shù)為零[10-11]。

3 計(jì)算結(jié)果

仿真計(jì)算采用Fluent軟件完成，計(jì)算獲得了發(fā)動(dòng)機(jī)的燃燒流場分布。圖 3為發(fā)動(dòng)機(jī)燃燒室內(nèi)頭部附近流線分布，根據(jù)圖可知，在燃燒室頭部以及針?biāo)^部附近存在兩塊回流區(qū)，即外側(cè)回流區(qū)和中心回流區(qū)，這是針?biāo)ㄊ絿娮⑵鲊婌F混合模式的典型特征，這種形式是實(shí)現(xiàn)高效燃燒和保證燃燒穩(wěn)定性的基礎(chǔ)。外側(cè)回流區(qū)主要是推進(jìn)劑向頭部的回流；中心回流區(qū)主要是推進(jìn)劑向針?biāo)^部回流，加劇推進(jìn)劑液滴的二次破碎和摻混燃燒效果，提高燃燒效率。

圖 3 燃燒室內(nèi)頭部附近流線分布

圖4為動(dòng)量比1工況下的溫度分布以及Y=100 mm位置截面的溫度分布，根據(jù)軸向的溫度分布可知，燃燒室壁面上游區(qū)域存在由氧化劑形成的液膜低溫區(qū)，這是因?yàn)橐貉鹾兔河蛷膰娮⑼ǖ懒鞒龊蟀l(fā)生碰撞，氧化劑獲得徑向速度分量，向上運(yùn)動(dòng)碰撞到燃燒室壁面后，分成兩部分，一部分向低壓區(qū)(燃燒室頭部)流動(dòng)，形成燃燒室頭部回流區(qū)；一部分沿燃燒室壁面向下游流動(dòng)，形成燃燒室壁面的液膜，該液膜對(duì)發(fā)動(dòng)機(jī)壁面起到冷卻的作用，對(duì)發(fā)動(dòng)機(jī)冷卻是有利的。流動(dòng)過程中煤油與氧的不斷燃燒消耗，使得液膜厚度變得越來越小，到燃燒室中下游附近消失；中心回流區(qū)內(nèi)的主要組分是未參與燃燒的煤油，在主流燃?xì)獾膸?dòng)下，未燃燒的煤油不斷向下游流動(dòng)，使得中心軸線附近出現(xiàn)低溫區(qū)，這個(gè)低溫區(qū)有利于降低針?biāo)^部的溫度，起到保護(hù)針?biāo)ǖ淖饔谩?/p>

圖 4 溫度分布

圖5為燃燒室中心軸線的溫度曲線，根據(jù)曲線可知，在針?biāo)^部附近溫度大約為1 445 K，在噴管的喉部附近溫度達(dá)到最高，約為2 500 K；圖 6為Y=100 mm，Z=0 mm位置處X方向的湍流強(qiáng)度分布，據(jù)圖可知，在徑向X=8 mm左右位置處湍流強(qiáng)度最大，該位置對(duì)應(yīng)溫度流場中的燃燒高溫區(qū)，該區(qū)域存在劇烈的化學(xué)反應(yīng)，湍流脈動(dòng)速度大，導(dǎo)致推進(jìn)劑混合劇烈，可以提高燃燒效率。

圖 5 軸向溫度曲線

圖 6 湍流強(qiáng)度分布

3.1 動(dòng)量比的影響

動(dòng)量比是指推進(jìn)劑的徑向射流動(dòng)量值與軸向射流動(dòng)量值的比，即

動(dòng)量比作為針?biāo)ㄊ絿娮⑵鞯闹匾阅軈?shù)，代表了噴注器的霧化混合程度，推進(jìn)劑的霧化混合主要是通過徑向射流液膜與軸向射流液膜碰撞實(shí)現(xiàn)，針?biāo)ㄊ絿娮⑵骺梢酝ㄟ^改變動(dòng)量比的大小，提高或者降低射流的穿透能力，增強(qiáng)或者減弱推進(jìn)劑的霧化混合程度，從而改變?nèi)紵矢叩屯瑫r(shí)。因此，本文分別設(shè)計(jì)動(dòng)量比為0.45，0.8，1.0，1.2和1.5的工況進(jìn)行數(shù)值計(jì)算，得出如下結(jié)果。

圖7為不同動(dòng)量比的溫度分布，根據(jù)圖可知，隨著動(dòng)量比增加，燃燒室壁面液膜厚度越來越小，壁面溫度越來越高；動(dòng)量比為1.5時(shí)，燃燒高溫區(qū)已經(jīng)貼近燃燒室壁面，說明隨著徑向射流動(dòng)量的增加，燃料和氧化劑的燃燒面向壁面移動(dòng)，使得燃燒高溫區(qū)靠近壁面。根據(jù)圖 8和圖 9的壓力分布可知，隨著動(dòng)量比增加燃燒室壓力呈現(xiàn)先增加后減小的趨勢(shì)，在動(dòng)量比為1時(shí)，燃燒室壓力最高，這是因?yàn)閯?dòng)量比小于1時(shí)，隨著動(dòng)量比增加，推進(jìn)劑徑向相對(duì)射流動(dòng)量增加，一定程度上提高了煤油和氧的混合程度，有利于燃燒，使的燃燒室壓力有增加的趨勢(shì)，但是，動(dòng)量比大于1時(shí)，隨著動(dòng)量比的增加，燃燒室壓力有減小的趨勢(shì)，說明動(dòng)量比大于1時(shí)，徑向射流動(dòng)量比較大，軸向射流動(dòng)量相對(duì)減小，使得軸向推進(jìn)劑組元穿透徑向射流的難度增加，導(dǎo)致碰撞、混合的效果減弱，發(fā)動(dòng)機(jī)燃燒效果略有減小，發(fā)動(dòng)機(jī)燃燒室壓力有所降低。圖10為不同動(dòng)量比的煤油質(zhì)量分?jǐn)?shù)分布云圖，根據(jù)圖可知，隨著動(dòng)量比增加，徑向射流動(dòng)量增強(qiáng)，使得煤油組分在燃燒室內(nèi)占據(jù)的區(qū)域增加，氧和煤油摻混的區(qū)域向燃燒室壁面靠近，導(dǎo)致燃燒的高溫區(qū)向燃燒室壁面靠近，提高了發(fā)動(dòng)機(jī)壁面的熱載荷，不利于發(fā)動(dòng)機(jī)燃燒室的熱防護(hù)。

圖 7 不同動(dòng)量比溫度分布

圖 8 不同動(dòng)量比壓力分布

圖 9 不同動(dòng)量比軸向壓力分布曲線

圖 10 不同動(dòng)量比煤油質(zhì)量分?jǐn)?shù)分布

3.2 燃燒效率分析

為了評(píng)估燃燒室內(nèi)能量轉(zhuǎn)換過程中的損失大小，引入燃燒效率的概念，火箭發(fā)動(dòng)機(jī)的燃燒效率是指實(shí)際特征速度與理論特征速度的比值[11]，即：

根據(jù)數(shù)值仿真結(jié)果計(jì)算出不同動(dòng)量比工況下的燃燒效率如表 2所示。根據(jù)表可知，動(dòng)量比為1時(shí)，發(fā)動(dòng)機(jī)的燃燒效率最高為96.65%。為了分析發(fā)動(dòng)機(jī)燃燒效率和動(dòng)量比之間的變化關(guān)系，探究發(fā)動(dòng)機(jī)具有最佳燃燒性能的工況點(diǎn)，根據(jù)表 2數(shù)據(jù)采用多項(xiàng)式擬合的方法擬合出燃燒效率與動(dòng)量比的關(guān)系曲線如圖 11所示，根據(jù)圖可知，隨著動(dòng)量比增加，燃燒效率呈現(xiàn)先增加后減小的趨勢(shì)，并且動(dòng)量比為1時(shí)燃燒效率最高，說明該狀態(tài)下推進(jìn)劑混合效果最好，具有最佳的燃燒性能。因此，在設(shè)計(jì)針?biāo)ㄊ絿娮⑵鲿r(shí)，為了提高發(fā)動(dòng)機(jī)性能，盡量選擇動(dòng)量比為1的工況點(diǎn)進(jìn)行設(shè)計(jì)。

表 2 不同動(dòng)量比燃燒效率

圖 11 燃燒效率與動(dòng)量比關(guān)系曲線

4 結(jié)論

本文針對(duì)液氧/煤油推進(jìn)劑的針?biāo)ㄊ桨l(fā)動(dòng)機(jī)進(jìn)行數(shù)值仿真，得出如下結(jié)果：

1)針?biāo)ㄊ桨l(fā)動(dòng)機(jī)在燃燒室內(nèi)形成兩個(gè)回流區(qū)，即燃燒室頭部回流區(qū)和針?biāo)^部回流區(qū)，回流區(qū)的存在有利于發(fā)動(dòng)機(jī)噴注器的熱防護(hù)；

2)針?biāo)ㄊ桨l(fā)動(dòng)機(jī)的噴注形式能夠在燃燒室壁面形成液膜，有利于燃燒室冷卻；

3)隨著動(dòng)量比增加，燃燒該高溫區(qū)向壁面靠近，增加了壁面熱載荷，對(duì)發(fā)動(dòng)機(jī)的熱防護(hù)是不利的；

4)動(dòng)量比接近1時(shí)，霧化混合效果最好，燃燒效率最高。