吳 睿，聶萬勝，蔡紅華，喬 野，馮 偉

（中國人民解放軍裝備學(xué)院航天裝備系，北京，101416）

偏二甲肼/四氧化二氮火箭發(fā)動機尾焰流場特性三維仿真研究

吳 睿，聶萬勝，蔡紅華，喬 野，馮 偉

（中國人民解放軍裝備學(xué)院航天裝備系，北京，101416）

以偏二甲肼/四氧化二氮（UDMH/NTO）火箭發(fā)動機為研究對象，采用κ-ε湍流模型，運用PISO算法分別對發(fā)動機內(nèi)流-場和尾焰流場進行三維仿真。采用相同方法計算液氫/液氧（LH2/LOX）火箭發(fā)動機尾焰，仿真結(jié)果和試驗結(jié)果吻合得較好，證明了計算模型的正確性與有效性；同時對比分析了UDMH/NTO發(fā)動機與LH2/LOX發(fā)動機尾焰流場特性。結(jié)果表明，兩者具有相似的溫度和馬赫數(shù)變化趨勢，但是UDMH/NTO發(fā)動機尾焰核心區(qū)溫度相對較低，而LH2/LOX發(fā)動機尾焰將更快衰減至亞聲速射流。

液體火箭發(fā)動機；三維數(shù)值仿真；尾焰流場特征

0 引 言

火箭發(fā)動機尾焰是由燃燒室中高溫、高壓的化學(xué)反應(yīng)產(chǎn)物經(jīng)過噴管加速后達到超聲速形成的。這些高溫、高速產(chǎn)物噴入靜止大氣后，壓縮周圍的空氣形成特征鮮明的波系結(jié)構(gòu)[1,2]。研究火箭發(fā)動機尾焰波系結(jié)構(gòu)對火箭和反導(dǎo)系統(tǒng)的研發(fā)具有重要意義。

現(xiàn)有的大型運載火箭仍有許多采用雙組元自燃推進劑，如歐空局的阿里安4系列運載火箭、美國的大力神系列運載火箭，以及中國的長征二號、長征三號系列運載火箭。同時液氫/液氧作為目前比沖最高的化學(xué)推進劑也廣泛應(yīng)用于大型運載火箭中。

尾焰研究方法主要有數(shù)值模擬和試驗兩種，國外在20世紀對尾焰流場進行了大量數(shù)值模擬和試驗。Vitkin等[3]詳細論述了火箭發(fā)動機用的物理數(shù)學(xué)模型，考慮了尾流中復(fù)燃、化學(xué)反應(yīng)、羽煙等的作用；George等[4]耦合CFD-DSMC方法對小推力發(fā)動機尾流流場進行了仿真。由于試驗成本高、模擬的條件有限[5]，因此中國針對尾焰的研究主要集中于數(shù)值模擬[6～9]；文獻[9]中則指出，尾焰三維模型計算結(jié)果較二維軸對稱模型計算結(jié)果有較大區(qū)別，其原因是二維軸對稱模型無法準確模擬燃氣流動，而采取三維模型才能得到較為準確的結(jié)果。

本文針對偏二甲肼/四氧化二氮（UDMH/NTO）火箭發(fā)動機，考慮液體火箭發(fā)動機噴嘴及隔板對發(fā)動機內(nèi)流場以及尾焰流場的影響，采用一種分塊計算三維尾焰流場的方法。即對于燃燒較為劇烈的發(fā)動機燃燒室及噴管收縮段，采用較小時間步長和精細網(wǎng)格以得到一個比較準確的內(nèi)流場，并將喉部參數(shù)作為尾焰流場質(zhì)量入口的邊界條件。而對于噴管擴張段以及外流場這個燃燒相對緩慢的區(qū)域，采用較大的時間步長和更加稀疏的網(wǎng)格進行計算，具有較好的針對性，合理分配了計算資源，同時可以得到較為準確的三維尾焰流場，具有一定的實際意義。

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

1.1 控制方程

控制方程采用三維N-S方程來描述，其質(zhì)量、動量和能量方程的基本形式為

式中 φ為通用變量；ρ為流體密度；U為速度矢量；Γφ為對應(yīng)于φ的擴散系數(shù)；Sφ為相應(yīng)的源項。

1.2 湍流模型

對于發(fā)動機內(nèi)流場，采用標準k-ε雙方程模型處理；對于近壁區(qū)采用標準壁面函數(shù)處理[10]；對于尾焰流場，根據(jù)文獻[11]，采用標準的k-ε雙方程模型處理。因為其相對于標準k-ε雙方程模型，可以更好地模擬超聲速噴流與周圍大氣的摻混與燃燒。對于近壁區(qū)仍采用標準壁面函數(shù)處理。

1.3 化學(xué)反應(yīng)模型

本文采用有限速率/渦耗散模型對尾焰流場特性進行計算，其中化學(xué)反應(yīng)速率由Arrhenius公式給出[10]：

式中 kf,r為正向化學(xué)反應(yīng)速率；Ar為指前因子；βr為溫度系數(shù)；Er為活化能；R為氣體常數(shù)；T為化學(xué)反應(yīng)溫度。

本文采用的化學(xué)反應(yīng)方程式為

1.4 求解方法

采用PISO算法求解。PISO算法適用于各種流動速度的定常/非定常、可壓/不可壓及兩相燃燒流動問題。同時對于瞬態(tài)問題，PISO算法總體效率較高，具有較明顯的優(yōu)勢[12,13]。

1.5 物性參數(shù)

燃燒過程中，各組分的物性參數(shù)隨著溫度變化而變化，為了提高數(shù)值計算精度，必須考慮物性參數(shù)隨溫度的變化。

各氣相組分的定壓比熱CP,i可采用溫度的分段擬合多項式計算，即：

式中 R為理想氣體常數(shù)；Mi為各氣相組分摩爾質(zhì)量；A0～A5為定壓比熱組分多項式各溫度項系數(shù)。

各氣體粘性系數(shù)和導(dǎo)熱系數(shù)采用分子動理論計算得到。

2 熱力計算模型

2.1 基本假設(shè)

a）推進劑燃燒過程是絕熱的，燃燒產(chǎn)物與外界沒有熱交換；

b）推進劑燃燒產(chǎn)物處于化學(xué)平衡狀態(tài)；

c）燃燒產(chǎn)物及其混合氣體都認為是完全氣體，符合完全氣體狀態(tài)方程；

d）燃燒產(chǎn)物在噴管內(nèi)是絕熱、等熵的膨脹過程。

2.2 最小吉布斯自由能法

最小吉布斯自由能法是在給定溫度和壓強條件下計算燃燒產(chǎn)物平衡組分最常用的方法之一。等溫、等壓條件下系統(tǒng)達到平衡狀態(tài)的判據(jù)[14]為

式中 k為推進劑中含有的不同元素的編號；Nk為1 kg推進劑中含有第k個元素的摩爾原子數(shù)；Akj為1 molj組分中含k個元素的摩爾原子數(shù)；nj為各組分摩爾數(shù)。

利用拉格朗日數(shù)乘法建立目標函數(shù)為

函數(shù)F的極值條件為

對方程組中對數(shù)項線性化處理后迭代求解，得到平衡組分摩爾數(shù)。

3 發(fā)動機內(nèi)流場網(wǎng)格及邊界條件

本文只考慮氣相燃燒，UDMH/NTO發(fā)動機的基本參數(shù)如表1所示，內(nèi)流場計算的邊界條件如表2所示，內(nèi)流場網(wǎng)格如圖1所示。

表1 發(fā)動機參數(shù)

由表1可知，熱力計算得到UDMH/NTO發(fā)動機燃燒室絕熱燃燒溫度為3339.871 K。

表2 發(fā)動機內(nèi)流場邊界條件

圖1 UDMH/NTO發(fā)動機網(wǎng)格

4 發(fā)動機內(nèi)流場計算結(jié)果及分析

計算得到的發(fā)動機內(nèi)流場如圖2所示。

采用單步化學(xué)反應(yīng)模型計算得到的UDMH/NTO發(fā)動機燃燒室X=0.15 m截面平均壓力值為671 kPa；截面平均溫度為3650 K，熱力計算偏差為9.29%。由于單步化學(xué)反應(yīng)模型的溫度高于實際溫度，因此發(fā)動機內(nèi)流場的計算結(jié)果可信。

圖2 UDMH/NTO發(fā)動機內(nèi)流場云圖

5 尾焰流場網(wǎng)格及邊界條件

為了節(jié)約計算資源，尾焰流場從喉部截面開始，選取1/6區(qū)域進行計算，其網(wǎng)格及邊界條件如圖3、表3所示。

圖3 尾焰流場網(wǎng)格及邊界條件

表3 尾焰場邊界條件

6 計算結(jié)果及分析

計算得到的UDMH/NTO發(fā)動機地面試車時的尾焰結(jié)構(gòu)，如圖4、圖5所示。

圖4 UDMH/NTO發(fā)動機尾焰溫度場云圖及等值線

圖5 UDMH/NTO發(fā)動機尾焰馬赫數(shù)云圖及等值線

從圖5可以看出，超聲速射流與大氣介質(zhì)相互作用，形成湍流混合區(qū)，當(dāng)其衰減至聲速，可以視作無激波平行射流。

6.1 模型有效性驗證

針對文獻[16]中給出的LH2/LOX一級發(fā)動機參數(shù)，采用上述模型進行三維仿真，得到的地面試車時尾焰結(jié)構(gòu)如圖6所示。

圖6 LH2/LOX發(fā)動機尾焰溫度云圖及其等值線

從圖6可以看出，LH2/LOX發(fā)動機尾焰射流噴出后壓縮較為劇烈，而UDMH/NTO發(fā)動機尾焰射流噴出后壓縮則相對溫和，這是由于2臺發(fā)動機噴口靜壓比Pe/Pa不同造成的，但二者都屬于過膨脹超聲速射流。從圖4～6中可以清楚地辨識尾焰射流邊界混合區(qū)、核心區(qū)、馬赫盤和倒錐形沖波等結(jié)構(gòu)。

本文得到的溫度場云圖與試驗結(jié)果比對效果如圖7所示。

圖7 試驗與仿真效果對比

由圖7可以看出，試驗與仿真結(jié)果吻合較好，因此上述計算模型有效。

6.2 尾焰形態(tài)結(jié)構(gòu)特征分析

運載火箭一級發(fā)動機地面工作時發(fā)動機噴口靜壓低于外界環(huán)境壓力(Pe＜Pa)，發(fā)動機工作在過膨脹狀態(tài)。理論上，欠膨脹超聲速射流的形態(tài)如圖8所示，過膨脹超聲速射流的形態(tài)如圖9所示[17]。

圖8 欠膨脹超聲速射流理論形態(tài)

圖9 過膨脹超聲速射流理論形態(tài)

發(fā)動機工作在過膨脹狀態(tài)時，燃氣在噴口內(nèi)過度膨脹，出口壓力小于外界壓力，噴出后遇到高于它的反壓作用而產(chǎn)生錐形波，氣流方向向內(nèi)折轉(zhuǎn)一個角度，使壓力提高與周圍介質(zhì)壓力平衡。由于周邊氣流都向內(nèi)折轉(zhuǎn)，在射流軸線上發(fā)生相撞，于是錐形波的頂點形成一個倒錐形沖波，使氣流流動平行于軸線。經(jīng)過這個波形后，氣流壓力又提高到大于介質(zhì)壓力，氣流流動方向平行于軸線。這時氣流已經(jīng)變成與欠膨脹射流相類似，所以隨后的流動就形成由膨脹-壓縮波段所組成的周期性震蕩結(jié)構(gòu)。

尾焰流場中心線上的溫度和馬赫數(shù)曲線如圖10、圖11所示。

圖10 中心線溫度曲線

圖11 中心線馬赫數(shù)曲線

氣流在噴管擴張段內(nèi)溫度逐漸降低，馬赫數(shù)逐漸升高，噴入大氣，溫度震蕩式上升后開始逐漸下降，而馬赫數(shù)的變化趨勢與溫度變化趨勢相反。從圖10中可知，UDMH/NTO發(fā)動機尾焰核心區(qū)溫度較低。

火箭燃氣射流流場由超聲速流和亞聲速流組成。超聲速燃氣射流從噴管噴出后，經(jīng)過膨脹和壓縮過程及復(fù)雜的激波系，最終使氣流流動變?yōu)閬喡曀倭?。燃氣軸線上的聲速點則是超聲速流和亞聲速流的轉(zhuǎn)折點。Piesit[18]建議利用下列經(jīng)驗公式計算燃氣射流中心線上聲速點的位置：

LH2/LOX發(fā)動機噴管出口燃氣中的主要成分是H2O，其平均比熱比為1.379；UDMH/NTO發(fā)動機噴管出口的主要成分是CO2、H2O和N2，其平均比熱比為1.222。采用上述經(jīng)驗公式計算發(fā)現(xiàn)，LH2/LOX發(fā)動機聲速尖位置距離噴口更近，即LH2/LOX火箭尾焰將更快衰減至亞聲速射流，這與圖11結(jié)果一致。

7 結(jié) 論

本文對UDMH/NTO和LH2/LOX火箭地面試車的尾焰流場進行了三維仿真計算，并進行比對，得到如下結(jié)論：

a）仿真結(jié)果分別與熱力計算結(jié)果、過膨脹燃氣射流理論基本一致，證明了模型的有效性和正確性；

b）LH2/LOX火箭一級發(fā)動機地面試車尾焰和UDMH/NTO火箭一級發(fā)動機地面試車尾焰都屬于過膨脹超聲速射流，都存在一個高溫的核心區(qū)，具有相似的溫度和馬赫數(shù)變化趨勢；

c）UDMH/NTO火箭尾焰核心區(qū)溫度較低，尾焰形態(tài)細長；

d） LH2/LOX火箭尾焰將更快衰減至亞聲速射流。

[1] 劉尊洋, 邵立, 汪亞夫,等. 復(fù)燃對固體火箭尾焰紅外輻射特性的影響[J]. 光學(xué)學(xué)報, 2013(6):32-39.

[2] 聶萬勝, 楊軍輝, 何浩波, 等. 液體火箭發(fā)動機尾噴焰紅外輻射特性[J]. 國防科技大學(xué)學(xué)報, 2005, 27(5):91-94.

[3] Vitkin E I, Karelin V G, Kirillov A A, et al. A physico-mathematical model of rocket exhaust plumes[J]. International Journal of Heat & Mass Transfer, 1997, 40(5):1227-1241.

[4] Boyd G.Simulation of nozzle plume flows using a combined CFDDSMC approach[R]. AIAA-99-3454, 1999.

[5] 徐春光. 復(fù)雜噴流流場數(shù)值模擬及應(yīng)用研究[D]. 長沙: 國防科學(xué)技術(shù)大學(xué), 2002.

[6] 郭景晶, 朱小蘭. 考慮化學(xué)反應(yīng)的真空羽流的數(shù)值模擬[J]. 宇航學(xué)報, 2001, 22(4):13-19.

[7] 帥永, 董士奎, 談和平. 數(shù)值模擬噴焰2.7微米紅外輻射特性[J]. 航空學(xué)報, 2005, 26(4):402-405.

[8] 蔡紅華，聶萬勝，鄭剛.液氧/煤油發(fā)動機噴管內(nèi)型面對尾焰特性影響[J].導(dǎo)彈與航天運載技術(shù), 2015(5):35-38.

[9] 李茂, 陳世哲, 陳春富. 火箭發(fā)動機地面水平試車尾流溫度場仿真分析[J]. 火箭推進, 2012, 38(6):29-34.

[10] 王振國.液體火箭發(fā)動機燃燒過程建模與數(shù)值仿真[M]. 北京: 國防工業(yè)出版社, 2012.

[11] 李茂，王占林.氫氧火箭發(fā)動機射流仿真與試驗臺熱防護[J]. 載人航天, 2014, 20(5): 421-426.

[12] Issa R I. Solution of the implicitly discretised fluid flow equations by operator-splitting[J].Journal of Computational Physics, 1986, 62(1):40-65.

[13] Ferzieger J L, Peric M. Computational methods forfluid dynamics[M]. Heidelberg:Springer-Verlag, 1996.

[14] Gordon S, Mcbride B J. Comprter program for calculation and fitting thermodynamic functions[R]. NASA RP-1271, 1992.

[15] Krueger A J, Minzner R A. A midlatitude ozone model for the 1976 U.S. standard atmosphere[J]. Journal of Geophysical Research Atmospheres, 1976,81(24): 4477-4481.

[16] 姜毅, 魏晰林, 陳苗. 發(fā)射動力學(xué)[M]. 北京: 北京理工大學(xué)出版社, 2015.

[17] 張福祥. 火箭燃氣射流動力學(xué)[M]. 北京: 國防工業(yè)出版社, 1988.

Three-dimensional Simulation Study of UDMH/NTO Rocket Engine Plume Flow Field Characteristics

Wu Rui, Nie Wan-sheng, Cai Hong-hua, Qiao Ye, FengWei
(Equipment Academy of PLA, Department of Space Equipment, Beijing, 101416)

This passage is modeled with UDMH/NTO rocket engine, and use the κ-ε turbulence model, PISO algorithm to simulate the rocket flow filed and exhaust plume flow field. With the compare between the simulation andexperiment of LH2/LOX rocket exhaust plume, the model is proved to be effective and correct. What’s more, the characteristic analysis of plume flow field between UDMH/NTO engine and LH2/LOX engine shows that tendency of temperature and mach number are the same, but the UDMH/NTO engine exhaust plume has a core of lower temperature and LH2/LOX engine exhaust plume attenuation to the subsonic flow faster.

Liquid rocket engine; Three-dimensional simulation; Plume folw field characteristics

1004-7182(2016)05-0074-06

10.7654/j.issn.1004-7182.20160516

V434

A

2016-05-01；

2016-08-14

國家自然科學(xué)基金項目（51206185，91441123）

吳 睿（1994-），男，碩士研究生，主要研究方向為發(fā)動機尾焰與紅外輻射特性