火箭發(fā)動機(jī)斜切噴管流場與推力特性的數(shù)值模擬研究①

金賀龍，王 浩，江 坤

(南京理工大學(xué) 能源與動力工程學(xué)院，南京 210094)

0 引言

噴管作為發(fā)動機(jī)動力系統(tǒng)的能量轉(zhuǎn)化裝置，其結(jié)構(gòu)設(shè)計的優(yōu)劣對火箭發(fā)動機(jī)性能具有重要的影響。發(fā)動機(jī)噴管按其結(jié)構(gòu)形式的不同，可分為直噴管和斜切噴管兩種形式。直噴管結(jié)構(gòu)作為發(fā)動機(jī)噴管的基本形式，許多學(xué)者已開展大量的研究工作，并取得一系列的研究成果；然而，針對斜切噴管的研究工作相對較少。斜切噴管在導(dǎo)彈的級間分離、推力終止及推力矢量控制技術(shù)中發(fā)揮著重要的作用。由于斜切噴管通常不具有常規(guī)的幾何外形，因此斜切噴管的設(shè)計要比常規(guī)直噴管復(fù)雜。通過對噴管型面的設(shè)計，斜切形式與安裝方式的選取，可實現(xiàn)不同的噴管性能，滿足不同的工程應(yīng)用需求。近年來，一些學(xué)者逐步開展了對斜切噴管的研究工作，并取得一定的進(jìn)展。鮑福廷等[1]針對斜切噴管推進(jìn)系統(tǒng)，提出一種優(yōu)化算法，以參數(shù)分析為基礎(chǔ)確定目標(biāo)函數(shù)的最優(yōu)值。劉君等[2]應(yīng)用計算流體力學(xué)軟件PHOENICS，通過求解二維湍流Navier-Stokes(N-S)方程，對有/無斜切的雙噴管固體火箭發(fā)動機(jī)內(nèi)流場進(jìn)行了數(shù)值模擬研究。樂貴高等[3]采用3階MUSCL TVD格式求解三維可壓縮平均Navier-Stokes，對兩種斜切噴管在噴口壓比為2.0、2.8、3.4和4.0條件下的發(fā)動機(jī)噴管燃?xì)馍淞鞯牧鲌鎏匦赃M(jìn)行數(shù)值模擬研究，得到噴管射流場流譜和流場特征參數(shù)的分布規(guī)律。徐瑋等[4]應(yīng)用有限體積法求解三維守恒型N-S方程組，對某非軸對稱斜切噴管的內(nèi)流場特性開展了數(shù)值模擬研究，重點分析噴管斜切角度變化對發(fā)動機(jī)推力的影響情況。伊進(jìn)寶等[5]也應(yīng)用數(shù)值方法對魚雷燃?xì)鉁u輪機(jī)的斜切噴管內(nèi)流場特性開展研究工作。上述學(xué)者研究主要集中在斜切噴管的流場特性方面，在斜切噴管推力等性能參數(shù)方面，一些學(xué)者也開展了一定的研究工作。林霞等[6]采用壓力積分法推導(dǎo)斜切噴管推力、推力偏轉(zhuǎn)角和力矩的計算公式，同時論證噴管斜切角度和擴(kuò)張半角對上述特征參數(shù)的影響。邢鵬濤等[7]采用數(shù)值分析與試驗論證相結(jié)合的方式，對某多斜切噴管結(jié)構(gòu)發(fā)動機(jī)的流場與推力特性進(jìn)行研究。劉沛等[8]、王立武等[9]也對斜切噴管的推力特性開展了相關(guān)的研究工作。由于斜切噴管的幾何形狀不再對稱，這將使得噴管內(nèi)燃?xì)鈪?shù)分布規(guī)律更為復(fù)雜。并且斜切噴管發(fā)動機(jī)存在推力偏轉(zhuǎn)現(xiàn)象，以上這些現(xiàn)象與常規(guī)噴管有很大不同。因此，有必要對發(fā)動機(jī)斜切噴管的流場與推力特性進(jìn)行深入研究。

本文針對斜切噴管火箭發(fā)動機(jī)的特點，通過數(shù)值模擬手段，采用三維非定?？蓧嚎sN-S方程組與Realizablek-ε湍流模型相結(jié)合的方法，并運(yùn)用混合網(wǎng)格技術(shù)，對斜切噴管火箭發(fā)動機(jī)噴管流場與推力特性進(jìn)行研究，同時分析了噴管斜切角度變化對發(fā)動機(jī)性能的影響。

1 計算模型與數(shù)值方法

1.1 幾何模型

本文所研究的固體火箭發(fā)動機(jī)噴管出口采用斜切結(jié)構(gòu)設(shè)計(如圖1所示)，發(fā)動機(jī)燃燒室軸線與噴管軸線之間成45°夾角。同時，由于空間結(jié)構(gòu)限制，發(fā)動機(jī)噴管入口省略收縮段結(jié)構(gòu)。此發(fā)動機(jī)燃燒室直徑22 mm，噴管喉部直徑6 mm，噴管擴(kuò)張半角15°。定義噴管斜切角度為發(fā)動機(jī)噴管軸線與噴管出口截面之間的夾角(噴管斜切角用β表示)，文中四種工況噴管斜切角度分別為45°、60°、75°和90°。不同角度斜切噴管軸線長度同為26.4 mm，噴管入口安裝角度相同。斜切噴管出口面積和喉部面積之比分別為14.21、10.78、 9.40和9。推進(jìn)劑選用雙基管狀推進(jìn)劑，推進(jìn)劑性能參數(shù)如表1所示。

(a) Structural diagram of 45° scarfed nozzle motor

表1 推進(jìn)劑性能參數(shù)

1.2 基本假設(shè)

為簡化計算，同時又能夠真實模擬斜切噴管內(nèi)的燃?xì)饬鲃?，本文引入如下假設(shè)條件：

(1)忽略化學(xué)反應(yīng)的具體過程，裝藥燃燒所生成的燃?xì)鉃橐环N混合氣體，并且該混合氣體滿足理想氣體狀態(tài)方程。

(2)燃?xì)饬鲃影牧鳡顟B(tài)。

(3)不考慮燃?xì)夂凸腆w壁面之間的輻射換熱現(xiàn)象及對熱結(jié)構(gòu)型面燒蝕帶來的影響。

(4)忽略發(fā)動機(jī)擋藥板結(jié)構(gòu)對流場計算的影響。

1.3 計算方法

1.3.1 控制方程

流場計算采用非定?？蓧嚎sN-S方程，以連續(xù)、動量和能量方程為基礎(chǔ)，考慮氣體比熱系數(shù)和導(dǎo)熱系數(shù)隨溫度的變化。在直角坐標(biāo)系中，可壓縮流體守恒型N-S方程可以表示成如下通用形式：

(1)

式中U=[ρ,ρu,ρv,ρw,e]T；ρ、u、v、w、e分別為燃?xì)獾拿芏取向速度、y向速度、z向速度和內(nèi)能；E、F和H為對流通量；Ev、Fv和Hv為黏性通量；S為源項。

1.3.2 湍流模型[10-12]

湍流模型選用Realizablek-ε模型：

Gk+Gb-ρε-YM+Sk

(2)

(3)

式中Gk為由層流速度梯度而產(chǎn)生的湍流動能；Gb為由浮力而產(chǎn)生的湍流動能；YM為由于在可壓縮湍流中，過渡擴(kuò)散而產(chǎn)生的波動；其他參數(shù)為相關(guān)特征常量與自定義源項。

1.3.3 火箭發(fā)動機(jī)噴管流量系數(shù)和推力[13-15]

定義發(fā)動機(jī)噴管流量系數(shù)φ為噴管實際質(zhì)量流量與噴管理論質(zhì)量流量之比。

噴管理論質(zhì)量流量計算公式為

(4)

式中Γ為氣動參數(shù)，是與燃?xì)獗葻岜萲有關(guān)的單值函數(shù)；pc為發(fā)動機(jī)燃燒室平均壓強(qiáng)值，Pa；Tc為發(fā)動機(jī)燃燒室燃?xì)鉁囟龋琄；At為發(fā)動機(jī)噴管喉部面積，m2。

火箭發(fā)動機(jī)推力為發(fā)動機(jī)內(nèi)、外表面全部作用力的合力。因此，發(fā)動機(jī)推力計算公式為

F=?AinpindA+?AexpexdA

(5)

式中pin和pex分別為作用于發(fā)動機(jī)內(nèi)、外表面上的壓強(qiáng)；dA為發(fā)動機(jī)內(nèi)外表面的微元面積。

由于發(fā)動機(jī)燃燒室內(nèi)的壓強(qiáng)變化較小，假設(shè)計算模型入口處燃?xì)鈮簭?qiáng)等于發(fā)動機(jī)前封頭處壓強(qiáng)；同時，作用在燃燒室圓筒段的壓強(qiáng)互相抵消。因此，進(jìn)行推力計算時，只對本文所建立的計算模型壁面壓強(qiáng)積分即可。

1.4 計算網(wǎng)格與邊界條件

取圖1(a)中斜切噴管火箭發(fā)動機(jī)后部作為發(fā)動機(jī)流場與推力特性研究的計算模型區(qū)域(如圖1(b)所示)，并對此區(qū)域進(jìn)行網(wǎng)格劃分。由于計算區(qū)域具有明顯的三維特性，并且噴管喉部入口區(qū)域幾何形狀較為復(fù)雜，因此對計算區(qū)域分段劃分網(wǎng)格是一種很好的方法。本文將計算區(qū)域分為3段：發(fā)動機(jī)燃燒室圓柱段、尾部收縮段和斜切噴管段。在燃燒室圓柱段和斜切噴管段建立結(jié)構(gòu)網(wǎng)格，對于幾何形狀較為復(fù)雜的尾部收縮段則建立非結(jié)構(gòu)網(wǎng)格。不同計算網(wǎng)格區(qū)域通過交界面進(jìn)行參數(shù)傳遞，網(wǎng)格劃分結(jié)果如圖2所示。

(a) Full grid

為獲取火箭發(fā)動機(jī)工作過程基本性能參數(shù)，采用4階龍哥庫塔法，自編計算程序求解火箭發(fā)動機(jī)內(nèi)彈道常微分方程組[16]，得到發(fā)動機(jī)內(nèi)彈道基本參數(shù)，同時通過平衡壓強(qiáng)工程計算公式進(jìn)行校正，發(fā)動機(jī)燃燒室平衡壓強(qiáng)約為7.15 MPa，程序計算結(jié)果如圖3所示。通過此壓強(qiáng)值、燃速方程和推進(jìn)劑燃燒面積，得到燃?xì)赓|(zhì)量流量為0.147 3 kg/s，并以此流量值作為數(shù)值計算的入口流量條件。計算區(qū)域邊界：在發(fā)動機(jī)燃燒室尾部燃?xì)馊肟谶吔缟?，給定燃?xì)馊肟谫|(zhì)量流量0.147 3 kg/s，燃?xì)饪倻?100 K。發(fā)動機(jī)噴管出口為超音速出口，此處出口參數(shù)按一階外推確定。對于發(fā)動機(jī)燃燒室與噴管壁面則采用絕熱、無滑移固壁邊界條件進(jìn)行處理。

圖3 發(fā)動機(jī)燃燒室壓強(qiáng)曲線

1.5 網(wǎng)格無關(guān)性

為確保數(shù)值模擬的準(zhǔn)確性，排除各個工況因網(wǎng)格數(shù)量不同而造成的結(jié)果失真，需對已建立的計算網(wǎng)格模型的無關(guān)性進(jìn)行驗證。驗證參量為發(fā)動機(jī)噴管軸線速度分布曲線和噴管質(zhì)量流量數(shù)值，取45°工況的計算網(wǎng)格模型，將網(wǎng)格單元數(shù)為2 340 000、2 510 000、 2 770 000三種方案進(jìn)行比較，計算結(jié)果如圖4和表2所示。由圖4可知，三種方案發(fā)動機(jī)噴管軸線速度分布曲線基本重合；同時，根據(jù)噴管質(zhì)量流量計算結(jié)果，工況1與工況2、3噴管質(zhì)量流場差異小于0.01%。表明當(dāng)計算網(wǎng)格數(shù)達(dá)到2 340 000以后，計算網(wǎng)格數(shù)量的增加對計算結(jié)果沒有影響。因此，為節(jié)省計算資源，提升計算效率，選取計算網(wǎng)格數(shù)2 340 000方案進(jìn)行數(shù)值計算。

表2 網(wǎng)格無關(guān)性

(a)Full curves (b)Local curves

2 計算結(jié)果與分析

2.1 斜切噴管火箭發(fā)動機(jī)內(nèi)部流場特性分析

圖5給出了45°斜切噴管火箭發(fā)動機(jī)燃燒室尾部與噴管流場的靜壓和馬赫數(shù)分布云圖。由圖5可知，在火箭發(fā)動機(jī)燃燒室尾部前端，燃燒室內(nèi)氣體壓強(qiáng)變化很小，燃?xì)饬鲃虞^為平緩，發(fā)動機(jī)燃燒室壓強(qiáng)約為7 MPa；當(dāng)燃?xì)饬鞯桨l(fā)動機(jī)噴管附近時，燃?xì)獠粩嗯蛎浖铀?，氣體速度增加，壓強(qiáng)減小。

(a)Pressure distribution contour (b)Mach number distribution contour

由于本文所研究的發(fā)動機(jī)沒有收縮段噴管結(jié)構(gòu)，這將使得氣流在噴管喉部入口附近產(chǎn)生較大的收縮圓弧形結(jié)構(gòu)，最大收縮弧形區(qū)域到達(dá)180°；并且此處壓強(qiáng)梯度和速度梯度較大，這種結(jié)構(gòu)設(shè)計不可避免導(dǎo)致噴管喉部附近附面層的厚度較大，即所謂的“縮頸現(xiàn)象”，造成較大的噴管流量損失。在發(fā)動機(jī)燃燒室壓強(qiáng)為7 MPa時，通過數(shù)值計算得到發(fā)動機(jī)噴管實際質(zhì)量流量為0.142 7 kg/s，而由理論性能參數(shù)公式(4)得到發(fā)動機(jī)噴管理論質(zhì)量流量為0.152 1 kg/s，噴管實際質(zhì)量流量為理論流量的0.938，即噴管流量系數(shù)為0.938，小于常規(guī)直噴管發(fā)動的噴管流量系數(shù)0.98。通過觀察圖5斜切噴管內(nèi)燃?xì)怦R赫數(shù)場和壓強(qiáng)場的分布特性發(fā)現(xiàn)，在斜切噴管內(nèi)流場存在兩個區(qū)域；在噴管結(jié)構(gòu)對稱區(qū)域，燃?xì)鈿鈩訁?shù)的流動沿噴管軸線左右對稱，而在非對稱區(qū)域，則存在氣體的單邊流動特性。非對稱區(qū)域和單邊流動特性將產(chǎn)生推力偏轉(zhuǎn)現(xiàn)象。

圖6給出了45°斜切噴管軸線速度和壓強(qiáng)分布曲線。由圖6可知，沿噴管軸線方向，燃?xì)馑俣认瓤焖僭黾樱_(dá)到第一個極值點后，速度下降；但很快速度到達(dá)一個極小值點，然后速度繼續(xù)快速上升達(dá)到第二個極值點；而后速度略有下降，再繼續(xù)增大直至噴管出口。而噴管軸線壓強(qiáng)分布與速度分布趨勢完全相反，并且兩曲線極值位置十分接近。燃?xì)馀蛎涀龉κ侨細(xì)馑俣仍黾拥闹鲃恿?，燃?xì)馀蛎涀龉?，速度增加，同時燃?xì)庾陨韷簭?qiáng)下降；而燃?xì)獗粔嚎s，壓強(qiáng)上升，氣體速度則會下降。產(chǎn)生上述流場參數(shù)變化的原因，可以根據(jù)膨脹-壓縮波理論進(jìn)行解釋。

圖6 45°斜切噴管發(fā)動機(jī)噴管軸線速度和壓強(qiáng)分布

本文所研究的發(fā)動機(jī)噴管結(jié)構(gòu)，噴管入口與噴管喉部采用直接連接形式進(jìn)行過度。這種過度形式結(jié)構(gòu)簡單，加工方便，但會導(dǎo)致“縮頸現(xiàn)象”，使噴管喉部附近邊界層厚度增大。噴管喉部邊界層呈現(xiàn)弧形形狀，在噴管喉部入口處氣體流動通道逐漸變窄，在喉部與擴(kuò)張段連接處附近氣體流動通道最小，之后氣體流動通道逐漸變寬。當(dāng)燃?xì)饬魅雵姽芎?，由于弧形邊界層的存在，燃?xì)饬鲃訉⑦B續(xù)內(nèi)折，由膨脹波與壓縮波理論可知，氣體連續(xù)內(nèi)折將會產(chǎn)生一系列的壓縮波，同時壓縮波匯聚將形成一道弱激波；不同方向的弱激波在噴管軸線附近交匯形成一個區(qū)域，燃?xì)饬鹘?jīng)此區(qū)域壓強(qiáng)上升，速度下降，形成局部高壓區(qū)。此區(qū)域?qū)?yīng)于噴管軸線速度下降的第一個極小值點。同時這些弱激波經(jīng)過噴管擴(kuò)張段壁面的反射作用，在下游噴管軸線上交匯，再次形成一個較強(qiáng)的激波交匯區(qū)域，導(dǎo)致燃?xì)馑俣鹊诙蜗陆担鴫簭?qiáng)增加。

2.2 不同角度斜切噴管內(nèi)流場特性與推力特性分析

為研究不同角度斜切噴管的流場與推力特性，取噴管安裝角度相同，噴管出口斜切角度分別為45°、60°、75°和90°四種工況進(jìn)行對比，不同角度噴管流場速度分布云圖如圖7所示。

(a)45° (b)60°

由圖7可知，當(dāng)發(fā)動機(jī)噴管安裝角度相同而噴管出口斜切角度不同時，噴管出口斜切角度不同并不會對發(fā)動機(jī)燃燒室與噴管內(nèi)速度場的整體分布情況產(chǎn)生明顯的影響。四種工況噴管喉部與噴管擴(kuò)張段對稱結(jié)構(gòu)部分的速度場分布狀況基本相同，而在噴管擴(kuò)張段非對稱部分，噴管速度分布存在一定的單邊現(xiàn)象。在非對稱區(qū)域，噴管軸線附近速度分布較為近似；而在噴管兩邊邊界處，四種工況速度分布存在差異，并且噴管斜切角度越小，速度分布差異越大。圖8為不同角度斜切噴管軸線壓強(qiáng)和速度分布曲線?？梢?，不同角度斜切噴管軸線壓強(qiáng)分布和速度分布基本重合。四種工況噴管內(nèi)燃?xì)饬鲃右?guī)律都是氣體先膨脹做功，燃?xì)饬魉僭黾?，壓?qiáng)數(shù)值逐漸降低；當(dāng)氣體流經(jīng)噴管擴(kuò)張段入口處第一個激波交匯區(qū)后，燃?xì)饬魉傧陆?，壓?qiáng)上升；之后氣體繼續(xù)加速膨脹直到燃?xì)饬鹘?jīng)第二個激波交匯區(qū)域，氣流速度再次下降，隨后燃?xì)饫^續(xù)膨脹，直到流出噴管。同時，不同角度斜切噴管軸線壓強(qiáng)和速度分布基本一致，表明噴管斜切角度變化并不會對噴管軸線氣動參數(shù)分布產(chǎn)生明顯影響。

(a)Pressure distribution curves (b)Velocity distribution curves

為研究不同角度斜切噴管的推力特性，同樣選用噴管斜切角度45°、60°、75°和90°這四種工況進(jìn)行對比分析。通過數(shù)值仿真與噴管特征參量的計算得到噴管特征參數(shù)，計算結(jié)果如表3所示。定義噴管軸向推力Fx為平行于發(fā)動機(jī)燃燒室軸線方向作用力，噴管側(cè)向推力Fy為垂直于燃燒室軸線方向作用力。定義推力偏轉(zhuǎn)角θ為發(fā)動機(jī)推力作用線與噴管軸線之間的夾角。對于不同角度的斜切噴管，隨著噴管斜切角度的增加，噴管中燃?xì)獾馁|(zhì)量流量數(shù)值相差不大，表明當(dāng)斜切噴管入口安裝角度相同時，僅改變噴管出口的斜切角度并不會對噴管中燃?xì)赓|(zhì)量流量大小產(chǎn)生明顯的影響。說明對于本文所研究的斜切噴管結(jié)構(gòu)，噴管入口安裝角度與連接方式對噴管中燃?xì)赓|(zhì)量流量起決定作用。此外，對于不同角度的斜切噴管，隨著噴管斜切角度的增加，發(fā)動機(jī)噴管軸向推力逐漸增加，側(cè)向推力逐漸減小。

表3 不同角度斜切噴管性能參數(shù)

圖9為不同角度斜切噴管推力特征參數(shù)分布圖。由圖9(a)可知，隨著噴管斜切角度的增加，發(fā)動機(jī)噴管軸向推力Fx近似于線性增大，而側(cè)向推力Fy則線性減?。划?dāng)噴管斜切角度為90°時(噴管擴(kuò)張段為直噴管結(jié)構(gòu))，法向推力和側(cè)向推力近似相等。但對于不同角度的斜切噴管，噴管總推力數(shù)值之間相差不大。

(a)Axial thrust and lateral thrust (b)Thrust deflection angle

通過比較不同角度斜切噴管發(fā)現(xiàn)，噴管采用斜切結(jié)構(gòu)設(shè)計時，45°、60°和75°三種工況同時存在推力偏轉(zhuǎn)現(xiàn)象。隨著噴管斜切角度的增加，發(fā)動機(jī)的推力偏轉(zhuǎn)角度不斷減小，當(dāng)噴管斜切角度為90°時，推力偏轉(zhuǎn)現(xiàn)象基本消失。對于斜切噴管的推力偏轉(zhuǎn)現(xiàn)象可做如下解釋：根據(jù)火箭發(fā)動機(jī)推力計算原理，發(fā)動機(jī)推力為發(fā)動機(jī)內(nèi)外表面全部作用力的合力。當(dāng)發(fā)動機(jī)噴管采用斜切結(jié)構(gòu)設(shè)計時，斜切噴管存在幾何對稱結(jié)構(gòu)與幾何非對稱結(jié)構(gòu)兩部分區(qū)域。在噴管幾何對稱部分，垂直于噴管軸線方向的作用力是可以相互抵消的；在幾何非對稱部分，將產(chǎn)生垂直于噴管軸線方向的作用合力，正是這部分作用力，導(dǎo)致推力偏轉(zhuǎn)現(xiàn)象的產(chǎn)生。并且隨著噴管斜切角度的增加，斜切噴管的不對稱程度逐漸降低。因此，推力偏轉(zhuǎn)角度不斷減小。

3 結(jié)論

(1)對于本文所研究的斜置斜切噴管發(fā)動機(jī)，由于噴管入口采用傾斜安裝形式，同時省略收縮段結(jié)構(gòu)，此種結(jié)構(gòu)設(shè)計將會帶來噴管質(zhì)量流量損失，噴管實際質(zhì)量流量為理論流量的0.938，明顯小于常規(guī)直噴管發(fā)動機(jī)的噴管流量系數(shù)。

(2)對于不同角度的斜切噴管，噴管喉部與噴管擴(kuò)張段對稱結(jié)構(gòu)部分的速度場分布狀況基本相同，而在噴管擴(kuò)張段非對稱部分，噴管速度分布存在一定的單邊現(xiàn)象。同時，燃?xì)庠趪姽軆?nèi)運(yùn)動時，會兩次穿過激波交匯的高壓強(qiáng)區(qū)域，通過此區(qū)域時，燃?xì)鈮簭?qiáng)增加，速度下降。

(3)隨著噴管斜切角度的增加，發(fā)動機(jī)噴管中燃?xì)獾馁|(zhì)量流量數(shù)值相差不大；噴管軸向推力Fx線性增大，側(cè)向推力Fy線性減小，推力偏轉(zhuǎn)角則從2.323°減小到0.063°，但噴管總推力數(shù)值之間相差不大。發(fā)動機(jī)噴管的非對稱結(jié)構(gòu)區(qū)域是產(chǎn)生推力偏轉(zhuǎn)現(xiàn)象的原因，并且斜切噴管不對稱結(jié)構(gòu)區(qū)域越大，則推力偏轉(zhuǎn)角度越大。