喬 野，聶萬勝，豐松江，蔡紅華，吳高楊

（中國人民解放軍裝備學院，北京101416）

0 引言

火箭發(fā)動機在燃燒室內消耗自身攜帶的推進劑通過噴管產生推力，這一過程伴隨大量發(fā)光發(fā)熱的燃氣從噴管噴出，并相繼形成第二、第三等若干個波陣面，即形成了火箭發(fā)動機的燃燒尾焰。尾焰具有高溫、高速、大流量的特點，在發(fā)射階段會對發(fā)射場產生很強的沖擊干擾和輻射干擾影響。

國外對火箭發(fā)動機尾焰的研究起步較早，但無論是數值模擬還是試驗測量，1990年前進展一直較為緩慢。隨著計算機技術和數值計算技術的發(fā)展，尾焰的研究取得了大量成果。文獻 [1-3]研究了尾焰流場對飛行器彈體受力的影響，文獻[4-6]研究了發(fā)動機尾焰流場形成與發(fā)展特點，文獻 [7-9]研究了發(fā)動機尾焰流場的沖擊效應，文獻 [10-11]研究了發(fā)動機尾焰流場的輻射特性。以上文獻為火箭發(fā)動機尾焰特性研究及發(fā)展提供了重要基礎。以液氫/液氧作為推進劑的液體火箭發(fā)動機反應能量遠遠大于液氧/煤油發(fā)動機和常規(guī)的偏二甲肼/四氧化二氮發(fā)動機。同時，由于氫的分子量極低，液氫/液氧發(fā)動機比推力比其他任何推進劑發(fā)動機的比推力都高，并具有無毒、無污染等優(yōu)點，因而也是世界各國爭相研究的對象，如美國的J-2發(fā)動機、日本的LE-7A發(fā)動機、蘇聯的RD-0120發(fā)動機[12]等都作為運載火箭芯一級動力系統發(fā)動機來使用。但目前，針對液氫/液氧發(fā)動機在地面發(fā)射階段的燃燒尾焰流場仿真計算還比較少，需要進一步研究。

為了分析航天發(fā)射中燃燒尾焰對發(fā)射場的影響，本文以液氫/液氧發(fā)動機為模型，考慮燃燒室內的燃燒反應，建立液氫/液氧發(fā)動機燃燒尾焰一體化仿真計算模型，得到發(fā)動機在地面發(fā)射階段燃燒尾焰流場的各項參數分布特點，為后期分析尾焰對航天發(fā)射的影響奠定基礎。

1 計算模型與數值方法

1.1 控制方程

控制方程采用二維N-S方程來描述，其質量、動量和能量方程基本形式為[13]：

式中：φ為通用變量，代表u，v，T等變量；ρ為流體密度；U為速度矢量；Гφ表示對應于φ的擴散系數；Sφ為相應的源項。

1.2 湍流模型

采用Realizable k-ε模型封閉控制方程，其對射流的仿真模擬可取得同試驗數據對比一致的結果[13]。

湍流動能k方程

湍流粘性系數方程

式中

湍流耗散率ε方程

1.3 燃燒模型

為簡化計算節(jié)省時間，本文考慮一步的全局化學反應：

發(fā)動機燃燒過程化學反應速率WCH采用湍流脈動機制REBU和Arrhenius機制RArr控制[14]：

1.4 計算方法

本文采用PISO算法進行壓力-速度耦合求解[15-16]。PISO算法基于對壓力-速度修正的高度近似，通過引入鄰近校正和傾斜校正克服了SIMPLE算法在完成壓力修正方程求解后由于節(jié)點速度以及相應通量不滿足動量平衡而不斷進行迭代的問題，提高了計算效率。

2 網格與邊界條件

圖1為尾焰計算網格。計算網格為以發(fā)動機及其尾焰遠場為一體的二維軸對稱計算區(qū)域。圖中OA為發(fā)動機入口，給定推進劑的質量流率，H2/O2氣氣噴注，混合比6.0，燃燒室室壓100 atm，噴管面積比5.0；ABC為發(fā)動機固壁面，采用壁面無滑移邊界條件；CDE為環(huán)境遠場，由于仿真主要研究地面發(fā)射階段尾焰流動特點，給定環(huán)境壓力1atm、環(huán)境溫度300 K，并且忽略來流速度；出口EF給定為壓力出口邊界條件，給定環(huán)境壓力以及環(huán)境溫度；OF為網格對稱軸，對稱軸上徑向速度為零，其余變量徑向梯度為零。

圖1 計算網格Fig.1 Mesh of plume calculation

3 計算結果分析

3.1 尾焰流場結構分析

圖2給出了尾焰流場的馬赫數分布。從圖中可以看出由于發(fā)動機處于欠膨脹工作狀態(tài)，在流場中可以捕捉到清晰的激波系結構。從局部放大的馬赫數分布中可以更加清晰的看到在發(fā)動機出口沿軸線外側首先形成了P-M膨脹波簇，用以實現尾焰內部高壓到環(huán)境低壓的順壓梯度；膨脹波由等壓流線反射為壓縮波，并匯聚成軸線內側的相交激波，相交激波同馬赫盤構成了高低壓場；相交激波穿過馬赫盤形成了反射激波，并以此構成典型的三波交點結構，反射激波同馬赫盤共同構成了局部高壓場；反射激波再經過尾焰射流邊界的反射，形成膨脹波，由此不斷循環(huán)，形成了典型的欠膨脹超音速射流近場激波系結構。這同文獻 [17]中給出的理論分析結果完全一致，如圖3所示，證明了計算方法的合理性和有效性。

圖2 馬赫數分布Fig.2 Mach-number distribution

圖3 典型欠膨脹燃氣射流流場結構[17]Fig.3 Typical structure of under-expansion gas flow field[17]

3.2 尾焰各項流場參數分析

圖4和圖5給出了不同時刻流場壓力云圖以及沿軸線上壓力分布。從圖4可以直觀看出尾焰壓力場的動態(tài)形成過程。在t=0.5～20 ms時，在發(fā)動機出口逐漸形成以發(fā)動機出口為球心的半球形波，此半球形波不斷干擾發(fā)動機出口外的靜止流體，所到之處會對該處的壓力場產生壓力突越（見圖5），由此可以判斷此半球形波為正激波，即尾焰流動過程中產生的初始沖擊波，且可近似認為此沖擊波在以恒定速度勻速傳播；隨著尾焰沖擊波對外界靜止流體的不斷作用，尾焰欠膨脹近場激波系結構逐步形成。隨著流動的進行，馬赫盤沿軸線不斷移動，其（正激波）強度不斷發(fā)展；t=20～150 ms時，尾焰欠膨脹近場激波系結構基本形成，且馬赫盤位置基本穩(wěn)定。

圖4 不同時刻流場壓力云圖Fig.4 Contour of pressure at different time

圖5 不同時刻軸線上壓力分布Fig.5 Pressure distribution on axial at different time

圖5中t=150 ms時，流場中最大壓力出現在燃燒室內，達到109.05 atm；燃氣流出發(fā)動機后，進入高低壓場，燃氣先膨脹后壓縮，壓力先降到低于環(huán)境壓力后，再在馬赫盤（正激波）的作用下急劇升高，進入高壓場；在高壓場中，壓力達到一個極大值（3.57 atm）后逐步降低；燃氣壓力通過以上方式不斷循環(huán)震蕩，隨著軸向距離的增大，激波系的強度逐漸降低，震蕩越來越小，最終衰減到1atm。

圖6(a)給出了流場的靜溫分布。圖中可以看出尾焰溫度場同壓力場具有相似的近場激波系結構，且流場高溫區(qū)集中在尾焰的高壓場中。隨著尾焰的流動，尾焰溫度逐漸降低。

圖6(b)為出口不同位置處流場溫度徑向分布。從圖中可以看出在X=2 m處為高低壓區(qū)，隨著徑向距離的增大，尾焰先在相交激波的作用下，溫度升高；隨后略微膨脹，溫度降低，但在邊界復燃作用的影響下，溫度有小幅上升；隨后在P-M膨脹波簇的作用下，溫度逐步降低到環(huán)境溫度。同理，X=4 m處在下一個波節(jié)的高低壓場，在相交激波和P-M膨脹波簇的作用下，溫度沿徑向先升高后降低，這里復燃的化學反應速率較低，對溫度產生的影響較小。在X=6～25 m時，隨著軸向距離的增加，尾焰中心區(qū)域的溫度“震蕩式”降低，而尾焰的厚度呈現先增加后減小的趨勢。

圖6 溫度分布Fig.6 Temperature distribution

圖7給出流場速度分布云圖。從圖中可以看到速度具有同壓力以及溫度相反的變化趨勢，這符合欠膨脹流動的基本規(guī)律。軸線上，燃氣在噴管中進行膨脹加速，在噴管出口外達到最大速度，為5 344 m/s；之后，燃氣在激波膨脹波的交替作用下，速度大小上下“震蕩”式降低，最終衰減為844 m/s。

圖8給出軸線上H2/O2/H2O體積分數分布，從圖中可以看出在發(fā)動機中H2和O2充分燃燒，O2余量很少。尾焰流場中，在X=1～12.8 m范圍內，尾焰燃氣的組分主要為H2O和H2；由于流動中的損耗以及和空氣中的氧發(fā)生復燃反應，H2的體積分數逐漸降低。在X=12.8～25 m范圍內，尾焰燃氣的主要成分為H2O和O2，且O2的體積分數逐漸增大，這是由于缺少H2而無法發(fā)生復燃反應，從而通過尾焰和空氣的摻混，增加了燃氣中O2含量。以上分析為之后研究發(fā)射階段尾焰輻射特性確定了研究對象，奠定了數據基礎。

圖7 速度分布Fig.7 Velocity distribution

圖8 軸線上H2/O2/H2O體積分數分布Fig.8 Volume fraction distribution of H2/O2/H2O on axis

4 結論

本文采用氣氣噴注、時均湍流模型對液氫液氧火箭發(fā)動機在發(fā)射階段的燃燒尾焰進行了一體化仿真計算，得到結論如下：

1）仿真結果同理論分析所得的欠膨脹燃氣射流近場激波系結構基本一致，證明了算法的有效性和正確性。

2）仿真計算實現了對尾焰壓力場形成過程的動態(tài)模擬，認為初始沖擊波是通過在發(fā)動機出口形成以發(fā)動機出口為球心的半球形波來對整個尾焰流場進行干擾的；這種沖擊波為正激波，會對壓力場產生壓力突越，且以恒定速度進行勻速傳播。發(fā)動機出口處，馬赫盤（正激波）隨時間的推移沿軸向發(fā)展運動，其激波強度逐漸增大，并最終趨于穩(wěn)定。

3）仿真得到了尾焰流場壓力、溫度、速度、馬赫數以及各個燃氣組分的分布情況，所得結果為進一步研究該型發(fā)動機發(fā)射階段尾焰的沖擊特性和輻射特性奠定了基礎。

[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]豐松江,聶萬勝.導彈尾焰對其飛行性能的影響研究[J].裝備學院學報,2006,17(5)：39-41.

[4]李軍,曹從詠,徐強.固體火箭燃氣射流近場形成與發(fā)展的數值模擬[J].推進技術,2003,24(5)：410-413.

[5]周松柏,郭正,高嵩,等.火箭發(fā)動機動態(tài)流場的數值模擬[J].推進技術,2007,28(2)：118-121.

[6]張光喜,周為民,張鋼錘,等.固體火箭發(fā)動機尾焰流場特性研究[J].固體火箭技術,2008,31(1)：18-23.

[7]周國儀,曹義華,胡繼忠.火箭尾噴流對帶孔平板沖擊流場的數值模擬[J].固體火箭技術,2001,24(2)：1-3.

[8]孫晉,曹從詠.火箭噴流對平板沖擊的數值模擬[J].南京理工大學學報,2002,26(4)：381-384.

[9]馬艷麗,姜毅,郝繼光,等.固體發(fā)動機燃氣射流對發(fā)射平臺沖擊效應研究[J].固體火箭技術,2010,33(4)：373-395.

[10]聶萬勝,楊軍輝,何浩波,等.液體火箭發(fā)動機尾噴焰紅外輻射特性[J].國防科技大學學報,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]丁兆波,孫繼國,路曉紅.國外典型大推力氫氧發(fā)動機推力室技術方案綜述[J].導彈與航天運載技術,2012(4)：28-30.

[13]陶文銓.數值傳熱學[M].西安：西安交通大學出版社,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]張福祥.火箭燃氣射流動力學[M].北京：國防工業(yè)出版社,1988.