王 革,張渴欣,周博成,陳 磊,關(guān) 奔*,汪根來

(1.哈爾濱工程大學(xué) 航天與建筑工程學(xué)院,哈爾濱 150001;2.中國航天科工集團(tuán)有限公司第六研究院41所,呼和浩特 010076)

0 引言

火箭發(fā)動機(jī)是一種可以直接產(chǎn)生推力的噴氣推進(jìn)動力裝置,其自帶燃料和氧化劑,可以在大氣層以外工作,成為宇宙航行和大氣層外飛行的主要?jiǎng)恿ρb置。相關(guān)研究[1]表明,發(fā)動機(jī)的性能指標(biāo)實(shí)際上總是低于理論上可達(dá)到的值,這是由于在發(fā)動機(jī)工作過程中,會存在燃料混合、燃燒及燃?xì)馀蛎洸怀浞值葐栴},造成發(fā)動機(jī)性能損失。在發(fā)動機(jī)所有損失中,噴管的非適應(yīng)性損失占據(jù)較大比例,最大可達(dá)到15%。非適應(yīng)性損失[2]是指傳統(tǒng)噴管在非設(shè)計(jì)高度工作時(shí),由于噴管內(nèi)部氣流的過膨脹或者欠膨脹而帶來的損失,為了減小這種性能損失,多種先進(jìn)的火箭噴管概念相繼出現(xiàn)[3-5],包括雙鐘型噴管、塞式噴管、膨脹偏流噴管和延伸噴管等,通常將它們稱為高度補(bǔ)償噴管。在這些高度補(bǔ)償噴管中,膨脹偏流噴管具有長度短、質(zhì)量輕、型面變化連續(xù)等優(yōu)點(diǎn),可以削減火箭發(fā)動機(jī)30%的體積和15%的結(jié)構(gòu)質(zhì)量[6],因此受到了廣泛的關(guān)注。

膨脹偏流噴管(ED噴管)的概念由RAO[7]于1960年提出。典型的ED噴管結(jié)構(gòu)包括外壁輪廓和中心塞錐兩個(gè)部分,這種結(jié)構(gòu)使發(fā)動機(jī)燃?xì)饬鹘?jīng)噴管喉部時(shí)的方向朝向噴管外壁輪廓,燃?xì)庠谌F拐角處發(fā)生膨脹,并在外壁輪廓的限制作用下偏轉(zhuǎn)為噴管軸向噴出。ED噴管存在兩種工作模態(tài),即低空開放模態(tài)和高空閉合模態(tài)。此兩種工作模態(tài)的轉(zhuǎn)換為ED噴管提供了優(yōu)秀的高度補(bǔ)償能力。為了獲得ED噴管內(nèi)更多的流動細(xì)節(jié),學(xué)者們進(jìn)行了大量的實(shí)驗(yàn)和數(shù)值研究[8-14]。研究發(fā)現(xiàn)在兩種工作模態(tài)下,ED噴管內(nèi)部均存在激波和膨脹波,在“開放”工作模態(tài),超音速核心流經(jīng)過一系列激波-膨脹波偏轉(zhuǎn)以保證波后壓力能夠適應(yīng)大氣環(huán)境壓力;在閉合工作模態(tài),氣流在膨脹波束作用下,流動方向發(fā)生了較大偏轉(zhuǎn),超、亞聲速區(qū)域之間的剪切層與噴管軸線相交,形成再附激波,使氣流以平行于軸線方向排出,噴管中激波的位置及強(qiáng)度決定了噴管的流動特性,進(jìn)而決定了噴管的性能。WANG等[15]數(shù)值研究了大擴(kuò)張比ED噴管內(nèi)的流動規(guī)律,重點(diǎn)分析了推力效率下降的流動機(jī)理。SCHOMBERG[16]和TAYLOR[17]等實(shí)驗(yàn)對比了ED噴管和傳統(tǒng)噴管、雙鐘型噴管的性能,發(fā)現(xiàn)ED噴管的效率要高于傳統(tǒng)收斂-擴(kuò)張型噴管,與雙鐘型噴管對比后發(fā)現(xiàn),在高壓比條件下,ED噴管推力性能更具優(yōu)勢,驗(yàn)證了ED噴管的先進(jìn)性。同時(shí),學(xué)者們也開展了結(jié)構(gòu)參數(shù)對噴管性能的影響研究。TAYLOR[18]和SCHOMBERG[19]等研究了喉部幾何形狀對ED噴管性能的影響,發(fā)現(xiàn)ED噴管允許被大幅縮短長度,從而節(jié)省結(jié)構(gòu)重量。PARK等[20]研究了不同壓比下具有不同塞錐折轉(zhuǎn)角的ED噴管的流動特性,JOHN等[21]數(shù)值研究了塞錐基底曲率半徑對噴管性能的影響。張琦等[22]采用數(shù)值手段,分析了塞錐尾椎角、初始擴(kuò)張角和收斂角對環(huán)喉式膨脹偏流噴管性能的影響,并獲得了其最佳值。

本研究擬采用數(shù)值手段,對環(huán)喉式膨脹偏流噴管(ATEDN)的內(nèi)流場進(jìn)行仿真計(jì)算,利用混合正交法[23]設(shè)計(jì)數(shù)值模擬實(shí)驗(yàn),探究擴(kuò)張比、入口壓力和塞錐基底半徑三個(gè)影響因素共同作用下噴管性能的變化規(guī)律,為后續(xù)ATEDN型面設(shè)計(jì)提供參考。

1 數(shù)值方法

1.1 幾何模型及邊界條件

ATEDN的計(jì)算模型主要由噴管外壁和中心塞錐兩部分組成,噴管外壁擴(kuò)張段母線采用三次樣條曲線進(jìn)行設(shè)計(jì),與喉部過渡圓弧保持相切,如圖1(a)所示。圖中,RH和Re分別為塞錐基底半徑和出口半徑,Gt為等效喉部寬度,噴管設(shè)計(jì)直徑為30 mm,初始擴(kuò)張角∠C=53.51°,塞錐下游延伸段長度L=120 mm,其余參數(shù)取值參照文獻(xiàn)[24]。ATEDN的計(jì)算區(qū)域如圖1(b)所示,噴管入口采用壓力入口邊界,將遠(yuǎn)場設(shè)置為壓力出口邊界,噴管壁面為絕熱無滑移壁面,考慮本文所使用的ATEDN為二維軸對稱結(jié)構(gòu),因此采用軸對稱邊界。外部流場沿噴管軸向方向長度設(shè)置為20Re,徑向長度設(shè)置為6Re。

(a)Calculation model for ATEDN

(b)Boundary conditions圖1 ATEDN計(jì)算模型和邊界條件設(shè)置示意圖Fig.1 Sketch of the calculation model for ATEDN and boundary conditions

1.2 數(shù)值模型

針對噴管地面實(shí)驗(yàn)工況,不考慮噴管外部高速來流,忽略ATEDN內(nèi)的化學(xué)反應(yīng)及多相流動,不考慮其內(nèi)部輻射傳熱、型面燒蝕等影響,忽略重力的影響,控制方程為可壓縮流的N-S方程。

質(zhì)量守恒方程:

(1)

動量守恒方程:

(2)

能量守恒方程:

(λT)+S

(3)

理想氣體狀態(tài)方程:

p=ρRT

(4)

其中,

Φ=·(τ·v)-(·τ)·v

式中ρ為流體的密度;t為時(shí)間;v為速度矢量;p為壓強(qiáng);τ為粘性應(yīng)力張量;μ為動力黏度;I為單位張量;S為變形速率張量;λ為熱導(dǎo)率;T為溫度;Φ為耗散函數(shù);R為特定氣體常數(shù)。

采用剪切應(yīng)力輸運(yùn)k-ω(SSTk-ω)湍流模型[25]將雷諾平均N-S方程進(jìn)行封閉,SSTk-ω模型是基于湍流動能k和比耗散率ω的兩方程模型,能夠較好地預(yù)測由于激波引起的分離流動特性[26-27],其表達(dá)式如下:

(5)

(6)

式中Gk為由層流速度梯度產(chǎn)生的湍流動能;Gω為ω方程產(chǎn)生的湍流動能;Гk和Гω為k和ω的擴(kuò)散率;Yk和Yω為擴(kuò)散產(chǎn)生的湍流。

選取基于密度的隱式求解器求解控制方程,差分格式使用AUSM格式[28],空間離散使用二階迎風(fēng)格式。計(jì)算在穩(wěn)態(tài)條件下進(jìn)行,庫朗數(shù)取值為2,在計(jì)算過程中,對噴管出口質(zhì)量流量進(jìn)行監(jiān)控,當(dāng)殘差下降至10-3以下,監(jiān)控參數(shù)無明顯變化時(shí),判定計(jì)算收斂。

1.2.1 網(wǎng)格無關(guān)性驗(yàn)證

采用結(jié)構(gòu)網(wǎng)格對圖1(a)所示計(jì)算域進(jìn)行網(wǎng)格劃分,并對貼壁網(wǎng)格進(jìn)行加密處理,網(wǎng)格生成結(jié)果如圖2所示。在噴管入口壓強(qiáng)為4 MPa,出口狀態(tài)為海平面、壓力為101325 Pa的工作條件下,利用40 000、80 000、160 000、260 000網(wǎng)格對ATEDN流動進(jìn)行數(shù)值模擬,繪制不同網(wǎng)格數(shù)量情況下的噴管壁面壓力變化曲線,如圖3所示。可以看出,隨著網(wǎng)格數(shù)量增加,曲線之間差異逐漸減小,當(dāng)網(wǎng)格數(shù)量達(dá)到160 000后,網(wǎng)格數(shù)量的繼續(xù)增加對壁面壓力曲線變化幾乎不再產(chǎn)生影響,因此本文采用數(shù)量為160 000的網(wǎng)格進(jìn)行數(shù)值計(jì)算。

圖2 計(jì)算區(qū)域的網(wǎng)格劃分Fig.2 Mesh layout of the computational domain

圖3 四種不同數(shù)量網(wǎng)格下噴管壁面壓強(qiáng)分布Fig.3 Pressure distributions along the nozzle shroud with four different mesh layouts

參照WAGNER和SCHLECHTRIEM的實(shí)驗(yàn)結(jié)果[12]對當(dāng)前使用的數(shù)值方法進(jìn)行驗(yàn)證,實(shí)驗(yàn)與數(shù)值模擬圖像對比見圖4。

(a) Experimental images (b) Numcrical images圖4 實(shí)驗(yàn)結(jié)果[12]與當(dāng)前數(shù)值結(jié)果對比Fig.4 Comparison of experimental images[12] with the present numerical results

由圖4可以看到,該數(shù)值方法可以很好地模擬出WAGNER和SCHLECHTRIEM的實(shí)驗(yàn)結(jié)果。在該實(shí)驗(yàn)條件下,當(dāng)NPR=10時(shí),噴管流動處于“開放”工作模態(tài),而在NPR=30時(shí),噴管流動處于“閉合”工作模態(tài)。數(shù)值方法很好地模擬了這兩種工作模態(tài)下的流動特征,激波和膨脹波清晰,從而證明了數(shù)值方法的準(zhǔn)確性。

1.2.2 參數(shù)設(shè)置

噴管工質(zhì)采用熱流燃?xì)?其總溫T0=2500 K、比熱容比γ=1.161、氣體常數(shù)Rg=408.35 J/(kg·K)。將燃?xì)庖暈槎▔罕葻崛莸目蓧嚎s理想氣體,熱導(dǎo)率按動能理論給定,粘性系數(shù)按薩瑟蘭粘性定律[29]給定。為了分析ATEDN在寬高度下的性能,選取0～31 km之間的8個(gè)典型高度工況,不同高度工況對應(yīng)的大氣參數(shù)通過《標(biāo)準(zhǔn)大氣參數(shù)公式》[30]進(jìn)行計(jì)算,計(jì)算結(jié)果如表1所示。

表1 大氣參數(shù)Table 1 Atmospheric parameters

1.3 噴管推力計(jì)算

ATEDN推力為其內(nèi)、外表面全部作用力的合力,計(jì)算公式為

F=Fin+Fex

(7)

其中,Fin由ATEDN入口、噴管外壁面的燃?xì)庖粋?cè)(圖1(a)中綠線)和中心塞錐壁面(圖1(a)中藍(lán)線)提供,表示燃?xì)鈱TEDN內(nèi)表面的作用力的合力;Fex由噴管外壁面外側(cè)提供,表示外界大氣對ATEDN外表面的作用力。

1.4 混合型正交表設(shè)計(jì)

現(xiàn)有火箭發(fā)動機(jī)的燃燒室壓強(qiáng)一般都保持在4～10 MPa[31]的水平上,因此選取噴管入口壓力的4個(gè)水平分別為4、6、8、10 MPa;考慮到ATEDN能夠滿足大擴(kuò)張比要求的結(jié)構(gòu)特性[13,22],選取擴(kuò)張比的4個(gè)水平分別為40、60、80、100;而考慮到塞錐基底半徑RH過大會降低發(fā)動機(jī)的有效載荷,在可接受范圍內(nèi)選取塞錐基底半徑的3個(gè)水平分別為30、34、38 mm。使用上述參數(shù)設(shè)計(jì)用于數(shù)值仿真的混合正交表,如表2所示,在之后分析中,將編號1、2、3、…所代表的組合方式稱為Case 1、Case 2、Case 3、…。

表2 混合型正交表Table 2 Mixed orthogonal table

2 結(jié)果與討論

2.1 ATEDN典型工作模態(tài)分析

為了分析ATEDN在不同高度下的流場特征,此處以Case 5為例,將其典型工作模態(tài)下的流場和壁面壓力分布進(jìn)行對比,如圖5所示?？梢钥吹?在低空時(shí)(0～1.95 km),ATEDN尾流處于“開放”模式(圖5(a)),高速燃?xì)庠谌F末端發(fā)生流動分離形成剪切層。剪切層兩側(cè)分別為沿噴管外壁的超聲速核心流和中心塞錐后方的亞聲速回流區(qū)。由塞錐基底尖角處發(fā)出的入射激波撞擊到噴管外壁面發(fā)生發(fā)射。從噴管外壁面壓力分布曲線可以看到,噴管收斂段壁面壓強(qiáng)有小幅下降,而進(jìn)入喉部后則壓強(qiáng)驟降,之后保持穩(wěn)定。在噴管中后段,激波反射造成壁面壓力的小幅度上升,隨后逐漸降低至環(huán)境壓力。將ATEDN的此工作模態(tài)稱為開放-固壁反射模態(tài)(OW模態(tài))。隨著工作高度增加 (3.01～4.21 km,見圖5(b)),噴管尾流雖然仍處于“開放”模式,但噴管內(nèi)部超聲速核心區(qū)面積增加,入射激波在噴管外壁面的反射點(diǎn)被推出噴管。此時(shí)激波發(fā)生自由邊界反射,噴管外壁面壓力在經(jīng)過喉部之后保持恒定。將ATEDN的此工作模態(tài)稱為開放-自由邊界反射模態(tài)(OF模態(tài))。當(dāng)工作高度進(jìn)一步增加(達(dá)5.57 km后,見圖5(c)),ATEDN尾流“閉合”,使得塞錐基底處的小面積回流區(qū)無法受到外界環(huán)境壓力變化影響,噴管壁面壓力分布曲線與OF模態(tài)相同且不隨工作高度的繼續(xù)增加而改變。此處將ATEDN的這種閉合工作模態(tài)簡稱為C模態(tài)。

圖5 Case 5在三種典型工作模態(tài)下的流場(上)及壁面壓力(下)分布Fig.5 Flow fields (upper row) and pressure distributions(lower row) of Case 5 at three typical operating modes

對表2所示各組合方式下ATEDN的工作過程進(jìn)行數(shù)值仿真,發(fā)現(xiàn)在0～31 km高度范圍內(nèi),Case 5的三種典型工作模態(tài)并不總是出現(xiàn)在所有的噴管工作過程中。比如Case 2不存在OW模態(tài),Case 3僅存在C模態(tài),經(jīng)過必要的算例補(bǔ)充后發(fā)現(xiàn),Case 10和14存在非常短暫的OF模態(tài),而Case 16的C模態(tài)發(fā)生在 7.2 km。將不同ATEDN所經(jīng)歷的典型工作模態(tài)及相對應(yīng)的工作高度范圍進(jìn)行歸納,可將其工作過程分為以下三類:

(1)噴管全程處于閉合工作模態(tài)(Case 3、4、7、8、9、11、13);

(2)噴管先后經(jīng)歷OF和C模態(tài)(Case 2、12、15);

(3)噴管先后經(jīng)歷OW、OF和C三種工作模態(tài)(Case 1、5、6、10、14、16),歸納結(jié)果可見表3。

表3 ATEDN典型工作模態(tài)及對應(yīng)的高度范圍Table 3 Typical operating modes and corresponding altitude ranges of ATEDNs

選取Case 3對第一種工作過程(即全程C模態(tài))的流場進(jìn)行分析,圖6為其在0、3.01、5.57 km高度工況下對應(yīng)的馬赫數(shù)分布圖?？梢钥吹?噴管尾流在海平面即處于“閉合”模態(tài),這主要與噴管的入口壓力和擴(kuò)張比有關(guān)。通過多個(gè)噴管流場的對比結(jié)果可以看出,當(dāng)擴(kuò)張比和塞錐基底半徑相同時(shí),入口壓力越大,噴管達(dá)到閉合模態(tài)所對應(yīng)的工作高度越低;當(dāng)入口壓力和塞錐基底半徑相同時(shí),擴(kuò)張比越小,噴管達(dá)到閉合模態(tài)時(shí)對應(yīng)的工作高度越低。處于第一種工作過程的噴管或是入口壓力足夠大,或是擴(kuò)張比足夠小,使其在海平面即處于“閉合”工作模態(tài)。全程C模態(tài)使ATEDN失去了高度補(bǔ)償能力,所以在噴管設(shè)計(jì)時(shí)應(yīng)該盡量避免噴管在海平面即處于閉合模態(tài)的情況發(fā)生。

圖6 Case 3在H取0、3.01、5.57 km的馬赫數(shù)分布圖Fig.6 Mach number distributions of Case 3 at H=0,3.01,5.57 km

圖7為Case 2、12、15(均處于第二種工作過程,即由OF轉(zhuǎn)變?yōu)镃模態(tài))在0、0.99、1.95 km工作高度的對應(yīng)的馬赫數(shù)分布云圖?？梢钥闯?在海平面工作條件下,激波反射點(diǎn)位于噴管外部,激波發(fā)生自由邊界反射,工作高度增加后,噴管尾流處于“閉合”模式。

圖7 Case 2、Case 12、Case 15(c)在H取0、0.99、1.95 km的馬赫數(shù)分布圖Fig.7 Mach number distributions of Case 2, Case 12 and Case 15 at H=0,0.99,1.95 km

繪制0～31 km高度范圍內(nèi)Case 2、Case 12和Case 15對應(yīng)的噴管推力變化曲線,如圖8所示?？梢钥闯?隨著工作高度的增加,Case 2、Case 12和Case 15對應(yīng)的噴管推力逐漸增加。由式(7)可知,噴管推力由兩部分組成,分別是燃?xì)鈱姽軆?nèi)表面的作用力Fin和外界大氣對噴管外表面的作用力Fex。隨著噴管工作高度的增加,噴管推力的Fin值變化很小,但是Fex的值卻顯著降低,此時(shí)噴管推力主要由噴管出口壓差所產(chǎn)生,因此,其推力逐漸增加。

圖8 第二種工作過程的推力隨高度變化曲線Fig.8 Thrust histories with altitude variation of the second operating process

選取Case 5對第三種工作過程(即由OW轉(zhuǎn)變?yōu)镺F再轉(zhuǎn)變?yōu)镃模態(tài))流場進(jìn)行分析,圖9為典型高度工況下其對應(yīng)的馬赫數(shù)分布圖?？梢钥闯?在0～ 1.95 km 處,噴管處于開放-固壁反射工作模態(tài),入射激波在噴管壁面發(fā)生反射,當(dāng)高度增加至3.01 km時(shí),剪切層下移,噴管內(nèi)高速核心區(qū)域增加,激波反射點(diǎn)移動至噴管出口外,剪切層上方的入射激波發(fā)生自由邊界反射。當(dāng)高度增加至5.57 km時(shí)剪切層與噴管軸線相交,塞錐下游的亞聲速回流區(qū)不再與外界環(huán)境連通,噴管進(jìn)入“閉合”工作模態(tài)。

圖9 Case 5在寬高度工況下的馬赫數(shù)分布圖Fig.9 Mach number distributions of Case 5 within wide height condition

圖10為第三種工作過程的噴管推力隨高度變化曲線?？梢钥吹?Case 1、5、6對應(yīng)的噴管推力先增加后減小后又增加;Case 16經(jīng)歷了2次推力增加和2次推力減小,Case 14經(jīng)歷了3次推力增加和2次推力減小。在該工作過程中,噴管推力變化情況較為復(fù)雜,由WANG等[15]的研究可知,這種推力的變化極大地受到噴管內(nèi)激波位置和激波強(qiáng)度的影響。

圖10 第三種工作過程的推力隨高度變化曲線Fig.10 Thrust histories with altitude variation of the third operating process

在由開放-固壁反射到開放-自由邊界反射的模態(tài)轉(zhuǎn)換過程中,噴管均出現(xiàn)推力驟減現(xiàn)象,這一現(xiàn)象可利用壁面壓強(qiáng)變化進(jìn)行解釋。以Case 5為例,其噴管所對應(yīng)的Fin和Fex隨工作高度變化曲線如圖11所示。

圖11 Case 5的Fin和Fex隨高度變化曲線Fig.11 Fin and Fex histories with altitude variation of Case 5

由圖11可見,隨著工作高度增加,Fin和Fex均逐漸減小。由于在模態(tài)轉(zhuǎn)換過程中(1.95～3.01 km),ATEDN噴管內(nèi)部的高速核心流區(qū)域擴(kuò)大,激波反射點(diǎn)被逐漸推出噴管,噴管擴(kuò)張段壁面的局部高壓區(qū)消失,使得Fin驟減,且其減小值要顯著大于Fex的減小值。因此,在該模態(tài)轉(zhuǎn)換過程中,噴管的總推力出現(xiàn)短暫的驟減現(xiàn)象。

2.2 正交法分析影響因素

(8)

表4為正交模擬計(jì)算結(jié)果。正交設(shè)計(jì)法所采集的數(shù)據(jù)常用極差分析法進(jìn)行處理。極差表示某一因素在所有水平條件下性能參數(shù)最大值與最小值之間的差值,此處的極差大小反映各因素對噴管性能影響的主次關(guān)系。極差的計(jì)算公式為

表4 正交模擬結(jié)果Table 4 Orthogonal simulation results

(9)

圖12 三個(gè)因素不同水平時(shí)高度積分平均比沖變化曲線Fig.12 Altitude-integral average specific impulse of the three factors at different levels

3 結(jié)論

本文對環(huán)喉式膨脹偏流噴管的流動過程進(jìn)行了數(shù)值仿真計(jì)算,采用正交法和極差分析法研究了擴(kuò)張比、入口壓力和塞錐基底半徑對噴管性能的影響及其重要程度,得出如下結(jié)論:

(1)ATEDN在不同高度工況工作時(shí),存在開放-固壁反射(OW模態(tài))、開放-自由邊界反射(OF模態(tài))和閉合(C模態(tài))三種典型工作模態(tài)。研究發(fā)現(xiàn),噴管在0～31 km高度范圍內(nèi)存在三類不同的工作過程:噴管全程處于C模態(tài);噴管先后經(jīng)歷OF和C模態(tài);噴管先后經(jīng)歷OW、OF和C三種工作模態(tài)。

(2)ATEDN的不同工作過程表現(xiàn)出不同的推力變化特性。對于第二種工作過程,ATEDN推力逐漸增加;對于第三種工作過程,噴管推力變化情況較為復(fù)雜,且在由OW模態(tài)到OF模態(tài)的轉(zhuǎn)換過程中,噴管推力出現(xiàn)驟減。