喬 野,聶萬(wàn)勝,豐松江,蔡紅華,吳高楊
(中國(guó)人民解放軍裝備學(xué)院,北京101416)
火箭發(fā)動(dòng)機(jī)在燃燒室內(nèi)消耗自身攜帶的推進(jìn)劑通過(guò)噴管產(chǎn)生推力,這一過(guò)程伴隨大量發(fā)光發(fā)熱的燃?xì)鈴膰姽車(chē)姵觯⑾嗬^形成第二、第三等若干個(gè)波陣面,即形成了火箭發(fā)動(dòng)機(jī)的燃燒尾焰。尾焰具有高溫、高速、大流量的特點(diǎn),在發(fā)射階段會(huì)對(duì)發(fā)射場(chǎng)產(chǎn)生很強(qiáng)的沖擊干擾和輻射干擾影響。
國(guó)外對(duì)火箭發(fā)動(dòng)機(jī)尾焰的研究起步較早,但無(wú)論是數(shù)值模擬還是試驗(yàn)測(cè)量,1990年前進(jìn)展一直較為緩慢。隨著計(jì)算機(jī)技術(shù)和數(shù)值計(jì)算技術(shù)的發(fā)展,尾焰的研究取得了大量成果。文獻(xiàn) [1-3]研究了尾焰流場(chǎng)對(duì)飛行器彈體受力的影響,文獻(xiàn)[4-6]研究了發(fā)動(dòng)機(jī)尾焰流場(chǎng)形成與發(fā)展特點(diǎn),文獻(xiàn) [7-9]研究了發(fā)動(dòng)機(jī)尾焰流場(chǎng)的沖擊效應(yīng),文獻(xiàn) [10-11]研究了發(fā)動(dòng)機(jī)尾焰流場(chǎng)的輻射特性。以上文獻(xiàn)為火箭發(fā)動(dòng)機(jī)尾焰特性研究及發(fā)展提供了重要基礎(chǔ)。以液氫/液氧作為推進(jìn)劑的液體火箭發(fā)動(dòng)機(jī)反應(yīng)能量遠(yuǎn)遠(yuǎn)大于液氧/煤油發(fā)動(dòng)機(jī)和常規(guī)的偏二甲肼/四氧化二氮發(fā)動(dòng)機(jī)。同時(shí),由于氫的分子量極低,液氫/液氧發(fā)動(dòng)機(jī)比推力比其他任何推進(jìn)劑發(fā)動(dòng)機(jī)的比推力都高,并具有無(wú)毒、無(wú)污染等優(yōu)點(diǎn),因而也是世界各國(guó)爭(zhēng)相研究的對(duì)象,如美國(guó)的J-2發(fā)動(dòng)機(jī)、日本的LE-7A發(fā)動(dòng)機(jī)、蘇聯(lián)的RD-0120發(fā)動(dòng)機(jī)[12]等都作為運(yùn)載火箭芯一級(jí)動(dòng)力系統(tǒng)發(fā)動(dòng)機(jī)來(lái)使用。但目前,針對(duì)液氫/液氧發(fā)動(dòng)機(jī)在地面發(fā)射階段的燃燒尾焰流場(chǎng)仿真計(jì)算還比較少,需要進(jìn)一步研究。
為了分析航天發(fā)射中燃燒尾焰對(duì)發(fā)射場(chǎng)的影響,本文以液氫/液氧發(fā)動(dòng)機(jī)為模型,考慮燃燒室內(nèi)的燃燒反應(yīng),建立液氫/液氧發(fā)動(dòng)機(jī)燃燒尾焰一體化仿真計(jì)算模型,得到發(fā)動(dòng)機(jī)在地面發(fā)射階段燃燒尾焰流場(chǎng)的各項(xiàng)參數(shù)分布特點(diǎn),為后期分析尾焰對(duì)航天發(fā)射的影響奠定基礎(chǔ)。
控制方程采用二維N-S方程來(lái)描述,其質(zhì)量、動(dòng)量和能量方程基本形式為[13]:
式中:φ為通用變量,代表u,v,T等變量;ρ為流體密度;U為速度矢量;Гφ表示對(duì)應(yīng)于φ的擴(kuò)散系數(shù);Sφ為相應(yīng)的源項(xiàng)。
采用Realizable k-ε模型封閉控制方程,其對(duì)射流的仿真模擬可取得同試驗(yàn)數(shù)據(jù)對(duì)比一致的結(jié)果[13]。
湍流動(dòng)能k方程
湍流粘性系數(shù)方程
式中
湍流耗散率ε方程
為簡(jiǎn)化計(jì)算節(jié)省時(shí)間,本文考慮一步的全局化學(xué)反應(yīng):
發(fā)動(dòng)機(jī)燃燒過(guò)程化學(xué)反應(yīng)速率WCH采用湍流脈動(dòng)機(jī)制REBU和Arrhenius機(jī)制RArr控制[14]:
本文采用PISO算法進(jìn)行壓力-速度耦合求解[15-16]。PISO算法基于對(duì)壓力-速度修正的高度近似,通過(guò)引入鄰近校正和傾斜校正克服了SIMPLE算法在完成壓力修正方程求解后由于節(jié)點(diǎn)速度以及相應(yīng)通量不滿(mǎn)足動(dòng)量平衡而不斷進(jìn)行迭代的問(wèn)題,提高了計(jì)算效率。
圖1為尾焰計(jì)算網(wǎng)格。計(jì)算網(wǎng)格為以發(fā)動(dòng)機(jī)及其尾焰遠(yuǎn)場(chǎng)為一體的二維軸對(duì)稱(chēng)計(jì)算區(qū)域。圖中OA為發(fā)動(dòng)機(jī)入口,給定推進(jìn)劑的質(zhì)量流率,H2/O2氣氣噴注,混合比6.0,燃燒室室壓100 atm,噴管面積比5.0;ABC為發(fā)動(dòng)機(jī)固壁面,采用壁面無(wú)滑移邊界條件;CDE為環(huán)境遠(yuǎn)場(chǎng),由于仿真主要研究地面發(fā)射階段尾焰流動(dòng)特點(diǎn),給定環(huán)境壓力1atm、環(huán)境溫度300 K,并且忽略來(lái)流速度;出口EF給定為壓力出口邊界條件,給定環(huán)境壓力以及環(huán)境溫度;OF為網(wǎng)格對(duì)稱(chēng)軸,對(duì)稱(chēng)軸上徑向速度為零,其余變量徑向梯度為零。
圖1 計(jì)算網(wǎng)格Fig.1 Mesh of plume calculation
圖2給出了尾焰流場(chǎng)的馬赫數(shù)分布。從圖中可以看出由于發(fā)動(dòng)機(jī)處于欠膨脹工作狀態(tài),在流場(chǎng)中可以捕捉到清晰的激波系結(jié)構(gòu)。從局部放大的馬赫數(shù)分布中可以更加清晰的看到在發(fā)動(dòng)機(jī)出口沿軸線(xiàn)外側(cè)首先形成了P-M膨脹波簇,用以實(shí)現(xiàn)尾焰內(nèi)部高壓到環(huán)境低壓的順壓梯度;膨脹波由等壓流線(xiàn)反射為壓縮波,并匯聚成軸線(xiàn)內(nèi)側(cè)的相交激波,相交激波同馬赫盤(pán)構(gòu)成了高低壓場(chǎng);相交激波穿過(guò)馬赫盤(pán)形成了反射激波,并以此構(gòu)成典型的三波交點(diǎn)結(jié)構(gòu),反射激波同馬赫盤(pán)共同構(gòu)成了局部高壓場(chǎng);反射激波再經(jīng)過(guò)尾焰射流邊界的反射,形成膨脹波,由此不斷循環(huán),形成了典型的欠膨脹超音速射流近場(chǎng)激波系結(jié)構(gòu)。這同文獻(xiàn) [17]中給出的理論分析結(jié)果完全一致,如圖3所示,證明了計(jì)算方法的合理性和有效性。
圖2 馬赫數(shù)分布Fig.2 Mach-number distribution
圖3 典型欠膨脹燃?xì)馍淞髁鲌?chǎng)結(jié)構(gòu)[17]Fig.3 Typical structure of under-expansion gas flow field[17]
圖4和圖5給出了不同時(shí)刻流場(chǎng)壓力云圖以及沿軸線(xiàn)上壓力分布。從圖4可以直觀(guān)看出尾焰壓力場(chǎng)的動(dòng)態(tài)形成過(guò)程。在t=0.5~20 ms時(shí),在發(fā)動(dòng)機(jī)出口逐漸形成以發(fā)動(dòng)機(jī)出口為球心的半球形波,此半球形波不斷干擾發(fā)動(dòng)機(jī)出口外的靜止流體,所到之處會(huì)對(duì)該處的壓力場(chǎng)產(chǎn)生壓力突越(見(jiàn)圖5),由此可以判斷此半球形波為正激波,即尾焰流動(dòng)過(guò)程中產(chǎn)生的初始沖擊波,且可近似認(rèn)為此沖擊波在以恒定速度勻速傳播;隨著尾焰沖擊波對(duì)外界靜止流體的不斷作用,尾焰欠膨脹近場(chǎng)激波系結(jié)構(gòu)逐步形成。隨著流動(dòng)的進(jìn)行,馬赫盤(pán)沿軸線(xiàn)不斷移動(dòng),其(正激波)強(qiáng)度不斷發(fā)展;t=20~150 ms時(shí),尾焰欠膨脹近場(chǎng)激波系結(jié)構(gòu)基本形成,且馬赫盤(pán)位置基本穩(wěn)定。
圖4 不同時(shí)刻流場(chǎng)壓力云圖Fig.4 Contour of pressure at different time
圖5 不同時(shí)刻軸線(xiàn)上壓力分布Fig.5 Pressure distribution on axial at different time
圖5中t=150 ms時(shí),流場(chǎng)中最大壓力出現(xiàn)在燃燒室內(nèi),達(dá)到109.05 atm;燃?xì)饬鞒霭l(fā)動(dòng)機(jī)后,進(jìn)入高低壓場(chǎng),燃?xì)庀扰蛎浐髩嚎s,壓力先降到低于環(huán)境壓力后,再在馬赫盤(pán)(正激波)的作用下急劇升高,進(jìn)入高壓場(chǎng);在高壓場(chǎng)中,壓力達(dá)到一個(gè)極大值(3.57 atm)后逐步降低;燃?xì)鈮毫νㄟ^(guò)以上方式不斷循環(huán)震蕩,隨著軸向距離的增大,激波系的強(qiáng)度逐漸降低,震蕩越來(lái)越小,最終衰減到1atm。
圖6(a)給出了流場(chǎng)的靜溫分布。圖中可以看出尾焰溫度場(chǎng)同壓力場(chǎng)具有相似的近場(chǎng)激波系結(jié)構(gòu),且流場(chǎng)高溫區(qū)集中在尾焰的高壓場(chǎng)中。隨著尾焰的流動(dòng),尾焰溫度逐漸降低。
圖6(b)為出口不同位置處流場(chǎng)溫度徑向分布。從圖中可以看出在X=2 m處為高低壓區(qū),隨著徑向距離的增大,尾焰先在相交激波的作用下,溫度升高;隨后略微膨脹,溫度降低,但在邊界復(fù)燃作用的影響下,溫度有小幅上升;隨后在P-M膨脹波簇的作用下,溫度逐步降低到環(huán)境溫度。同理,X=4 m處在下一個(gè)波節(jié)的高低壓場(chǎng),在相交激波和P-M膨脹波簇的作用下,溫度沿徑向先升高后降低,這里復(fù)燃的化學(xué)反應(yīng)速率較低,對(duì)溫度產(chǎn)生的影響較小。在X=6~25 m時(shí),隨著軸向距離的增加,尾焰中心區(qū)域的溫度“震蕩式”降低,而尾焰的厚度呈現(xiàn)先增加后減小的趨勢(shì)。
圖6 溫度分布Fig.6 Temperature distribution
圖7給出流場(chǎng)速度分布云圖。從圖中可以看到速度具有同壓力以及溫度相反的變化趨勢(shì),這符合欠膨脹流動(dòng)的基本規(guī)律。軸線(xiàn)上,燃?xì)庠趪姽苤羞M(jìn)行膨脹加速,在噴管出口外達(dá)到最大速度,為5 344 m/s;之后,燃?xì)庠诩げㄅ蛎洸ǖ慕惶孀饔孟?,速度大小上下“震蕩”式降低,最終衰減為844 m/s。
圖8給出軸線(xiàn)上H2/O2/H2O體積分?jǐn)?shù)分布,從圖中可以看出在發(fā)動(dòng)機(jī)中H2和O2充分燃燒,O2余量很少。尾焰流場(chǎng)中,在X=1~12.8 m范圍內(nèi),尾焰燃?xì)獾慕M分主要為H2O和H2;由于流動(dòng)中的損耗以及和空氣中的氧發(fā)生復(fù)燃反應(yīng),H2的體積分?jǐn)?shù)逐漸降低。在X=12.8~25 m范圍內(nèi),尾焰燃?xì)獾闹饕煞譃镠2O和O2,且O2的體積分?jǐn)?shù)逐漸增大,這是由于缺少H2而無(wú)法發(fā)生復(fù)燃反應(yīng),從而通過(guò)尾焰和空氣的摻混,增加了燃?xì)庵蠴2含量。以上分析為之后研究發(fā)射階段尾焰輻射特性確定了研究對(duì)象,奠定了數(shù)據(jù)基礎(chǔ)。
圖7 速度分布Fig.7 Velocity distribution
圖8 軸線(xiàn)上H2/O2/H2O體積分?jǐn)?shù)分布Fig.8 Volume fraction distribution of H2/O2/H2O on axis
本文采用氣氣噴注、時(shí)均湍流模型對(duì)液氫液氧火箭發(fā)動(dòng)機(jī)在發(fā)射階段的燃燒尾焰進(jìn)行了一體化仿真計(jì)算,得到結(jié)論如下:
1)仿真結(jié)果同理論分析所得的欠膨脹燃?xì)馍淞鹘鼒?chǎng)激波系結(jié)構(gòu)基本一致,證明了算法的有效性和正確性。
2)仿真計(jì)算實(shí)現(xiàn)了對(duì)尾焰壓力場(chǎng)形成過(guò)程的動(dòng)態(tài)模擬,認(rèn)為初始沖擊波是通過(guò)在發(fā)動(dòng)機(jī)出口形成以發(fā)動(dòng)機(jī)出口為球心的半球形波來(lái)對(duì)整個(gè)尾焰流場(chǎng)進(jìn)行干擾的;這種沖擊波為正激波,會(huì)對(duì)壓力場(chǎng)產(chǎn)生壓力突越,且以恒定速度進(jìn)行勻速傳播。發(fā)動(dòng)機(jī)出口處,馬赫盤(pán)(正激波)隨時(shí)間的推移沿軸向發(fā)展運(yùn)動(dòng),其激波強(qiáng)度逐漸增大,并最終趨于穩(wěn)定。
3)仿真得到了尾焰流場(chǎng)壓力、溫度、速度、馬赫數(shù)以及各個(gè)燃?xì)饨M分的分布情況,所得結(jié)果為進(jìn)一步研究該型發(fā)動(dòng)機(jī)發(fā)射階段尾焰的沖擊特性和輻射特性奠定了基礎(chǔ)。
[1]LEE Y K,RAGHUNATHAN S,BENARD E,et al.Plume interference effect on missile bodies and control,AIAA 2003-1241[R].Reston,USA:AIAA,2003.
[2]AVITAL G,POMPAN J,MACALES J,et al.Experimental and CFD study of rocket plume effects on missile longitudinalaerodynamicstability,AIAA2004-5196[R].Reston,USA:AIAA,2004.
[3]豐松江,聶萬(wàn)勝.導(dǎo)彈尾焰對(duì)其飛行性能的影響研究[J].裝備學(xué)院學(xué)報(bào),2006,17(5):39-41.
[4]李軍,曹從詠,徐強(qiáng).固體火箭燃?xì)馍淞鹘鼒?chǎng)形成與發(fā)展的數(shù)值模擬[J].推進(jìn)技術(shù),2003,24(5):410-413.
[5]周松柏,郭正,高嵩,等.火箭發(fā)動(dòng)機(jī)動(dòng)態(tài)流場(chǎng)的數(shù)值模擬[J].推進(jìn)技術(shù),2007,28(2):118-121.
[6]張光喜,周為民,張鋼錘,等.固體火箭發(fā)動(dòng)機(jī)尾焰流場(chǎng)特性研究[J].固體火箭技術(shù),2008,31(1):18-23.
[7]周?chē)?guó)儀,曹義華,胡繼忠.火箭尾噴流對(duì)帶孔平板沖擊流場(chǎng)的數(shù)值模擬[J].固體火箭技術(shù),2001,24(2):1-3.
[8]孫晉,曹從詠.火箭噴流對(duì)平板沖擊的數(shù)值模擬[J].南京理工大學(xué)學(xué)報(bào),2002,26(4):381-384.
[9]馬艷麗,姜毅,郝繼光,等.固體發(fā)動(dòng)機(jī)燃?xì)馍淞鲗?duì)發(fā)射平臺(tái)沖擊效應(yīng)研究[J].固體火箭技術(shù),2010,33(4):373-395.
[10]聶萬(wàn)勝,楊軍輝,何浩波,等.液體火箭發(fā)動(dòng)機(jī)尾噴焰紅外輻射特性[J].國(guó)防科技大學(xué)學(xué)報(bào),2005,27(5):91-94.
[11]NIE Wansheng,FENG Songjiang.Numerical simulation of liquid rocket exhaust plume radiation,AIAA 2007-4413[R].Reston,USA:AIAA,2007.
[12]丁兆波,孫繼國(guó),路曉紅.國(guó)外典型大推力氫氧發(fā)動(dòng)機(jī)推力室技術(shù)方案綜述[J].導(dǎo)彈與航天運(yùn)載技術(shù),2012(4):28-30.
[13]陶文銓.數(shù)值傳熱學(xué)[M].西安:西安交通大學(xué)出版社,1988.
[14]FENG Songjiang,CHENG Yufeng,NIE Wansheng.Combustion instability analysis in a hydrogen-oxygen rocketengine,AIAA2009-4864[R].Reston,USA:AIAA,2009.
[15]ISSA R I.Solution of implicitly discretized fluid flow equations by operator splitting[R].Journal of Computational Physics,1986(62):40-65.
[16]FERZIEGER J L,PERIC M.Computational methods for fluid dynamics[M].Heidelberg:Springer-Verlag,1996.
[17]張福祥.火箭燃?xì)馍淞鲃?dòng)力學(xué)[M].北京:國(guó)防工業(yè)出版社,1988.