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

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

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

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

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

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

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

液體火箭發(fā)動(dòng)機(jī)；三維數(shù)值仿真；尾焰流場(chǎng)特征

0 引 言

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

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

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

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

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

1.1 控制方程

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

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

1.2 湍流模型

對(duì)于發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)，采用標(biāo)準(zhǔn)k-ε雙方程模型處理；對(duì)于近壁區(qū)采用標(biāo)準(zhǔn)壁面函數(shù)處理[10]；對(duì)于尾焰流場(chǎng)，根據(jù)文獻(xiàn)[11]，采用標(biāo)準(zhǔn)的k-ε雙方程模型處理。因?yàn)槠湎鄬?duì)于標(biāo)準(zhǔn)k-ε雙方程模型，可以更好地模擬超聲速噴流與周圍大氣的摻混與燃燒。對(duì)于近壁區(qū)仍采用標(biāo)準(zhǔn)壁面函數(shù)處理。

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

本文采用有限速率/渦耗散模型對(duì)尾焰流場(chǎng)特性進(jìn)行計(jì)算，其中化學(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算法適用于各種流動(dòng)速度的定常/非定常、可壓/不可壓及兩相燃燒流動(dòng)問題。同時(shí)對(duì)于瞬態(tài)問題，PISO算法總體效率較高，具有較明顯的優(yōu)勢(shì)[12,13]。

1.5 物性參數(shù)

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

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

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

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

2 熱力計(jì)算模型

2.1 基本假設(shè)

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

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

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

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

2.2 最小吉布斯自由能法

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

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

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

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

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

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

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

表1 發(fā)動(dòng)機(jī)參數(shù)

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

表2 發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)邊界條件

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

4 發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)計(jì)算結(jié)果及分析

計(jì)算得到的發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)如圖2所示。

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

圖2 UDMH/NTO發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)云圖

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

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

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

表3 尾焰場(chǎng)邊界條件

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

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

圖4 UDMH/NTO發(fā)動(dòng)機(jī)尾焰溫度場(chǎng)云圖及等值線

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

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

6.1 模型有效性驗(yàn)證

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

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

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

本文得到的溫度場(chǎng)云圖與試驗(yàn)結(jié)果比對(duì)效果如圖7所示。

圖7 試驗(yàn)與仿真效果對(duì)比

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

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

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

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

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

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

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

圖10 中心線溫度曲線

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

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

火箭燃?xì)馍淞髁鲌?chǎng)由超聲速流和亞聲速流組成。超聲速燃?xì)馍淞鲝膰姽車姵龊螅?jīng)過膨脹和壓縮過程及復(fù)雜的激波系，最終使氣流流動(dòng)變?yōu)閬喡曀倭鳌Ｈ細(xì)廨S線上的聲速點(diǎn)則是超聲速流和亞聲速流的轉(zhuǎn)折點(diǎn)。Piesit[18]建議利用下列經(jīng)驗(yàn)公式計(jì)算燃?xì)馍淞髦行木€上聲速點(diǎn)的位置：

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

7 結(jié) 論

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

a）仿真結(jié)果分別與熱力計(jì)算結(jié)果、過膨脹燃?xì)馍淞骼碚摶疽恢拢C明了模型的有效性和正確性；

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

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

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

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

[2] 聶萬勝, 楊軍輝, 何浩波, 等. 液體火箭發(fā)動(dòng)機(jī)尾噴焰紅外輻射特性[J]. 國(guó)防科技大學(xué)學(xué)報(bào), 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ù)雜噴流流場(chǎng)數(shù)值模擬及應(yīng)用研究[D]. 長(zhǎng)沙: 國(guó)防科學(xué)技術(shù)大學(xué), 2002.

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

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

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

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

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

[11] 李茂，王占林.氫氧火箭發(fā)動(dòng)機(jī)射流仿真與試驗(yàn)臺(tái)熱防護(hù)[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ā)射動(dòng)力學(xué)[M]. 北京: 北京理工大學(xué)出版社, 2015.

[17] 張福祥. 火箭燃?xì)馍淞鲃?dòng)力學(xué)[M]. 北京: 國(guó)防工業(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

國(guó)家自然科學(xué)基金項(xiàng)目（51206185，91441123）

吳 睿（1994-），男，碩士研究生，主要研究方向?yàn)榘l(fā)動(dòng)機(jī)尾焰與紅外輻射特性