王大銳，梁 薇，馬麗娜，姚世強

基于有限體積法的液體火箭發(fā)動機燃?xì)饧t外特性仿真

王大銳，梁 薇，馬麗娜，姚世強

（北京航天動力研究所，北京，100076）

運用有限體積法對某軌姿控發(fā)動機典型工況燃?xì)饬鲌黾t外輻射特性進行仿真計算。首先通過Fluent軟件對發(fā)動機典型工況燃?xì)饬鲌鲞M行計算，再利用最新的光譜數(shù)據(jù)庫對主要輻射組分吸收系數(shù)進行計算，最后對燃?xì)饧t外輻射特性進行計算；重點分析不同波長、不同天頂角、不同波段下燃?xì)饧t外輻射強度分布規(guī)律。研究結(jié)果表明：所提計算方法可靠，紅外特性計算結(jié)果準(zhǔn)確，有利于掌握液體火箭發(fā)動機紅外輻射特性計算和分析方法。

軌姿控發(fā)動機；紅外輻射；仿真計算；逐線法

0 引 言

液體火箭發(fā)動機是衛(wèi)星、火箭等航天器的重要輻射源，因此掌握發(fā)動機紅外輻射特性對航天器紅外輻射研究具有重要意義，中國采用多種輻射計算方法對該問題進行研究。Surzhikov等[1]運用蒙特卡羅法對發(fā)動機尾焰輻射特性和底部加熱效應(yīng)進行研究；阮立明等[2]提出了源項六流法對導(dǎo)彈尾焰紅外輻射特性的研究方法；張小英等[3]利用離散坐標(biāo)法，結(jié)合新的光譜吸收系數(shù)計算方法，對不同飛行高度下羽流紅外輻射特性的影響進行了研究；劉尊洋等[4]利用有限體積法（Finite Volume Method，F(xiàn)VM）研究了復(fù)燃對固體火箭發(fā)動機尾焰紅外輻射特性的影響。

本文采用有限體積法與燃?xì)饬鲌鰯?shù)值計算相結(jié)合的方法對軌姿控發(fā)動機燃?xì)饬鲌黾t外輻射特性進行了計算和分析。通過Fluent軟件對發(fā)動機典型工況下軸對稱二維網(wǎng)格燃?xì)饬鲌鲞M行計算?；贖ΙTRAN2012和HΙTEMP2010光譜數(shù)據(jù)庫，通過自編逐線法程序?qū)2O、CO2、CO主要輻射組分吸收系數(shù)進行了計算。最后利用流場計算溫度、壓力、吸收組分分布以及各組分吸收系數(shù)等數(shù)據(jù)，采用燃?xì)廨椛溆嬎愠绦驅(qū)Σ煌r下的燃?xì)饬鲌黾t外輻射特性進行計算，并分析了不同波長、不同天頂角、不同波段下燃?xì)饧t外輻射分布規(guī)律。

1 流場計算模型

1.1 流場控制方程

質(zhì)量守恒方程為

式中 ρ為混合氣體密度；iY為組分i的質(zhì)量分?jǐn)?shù)；v為速度矢量；iD為組分i的擴散系數(shù)。

動量守恒方程為

式中 P為流體微元體上的靜壓力；τ為微元體上粘性應(yīng)力張量，表達(dá)式為

式中 μ為湍流粘性系數(shù)；iu，ju，ku分別為流體在3個垂直方向的速度分量；ijδ為流體變形速率張量。

能量守恒方程為

湍流模型采用RNGkε?（Renormalized Groupkε?）模型，該模型對于瞬變流和流線彎曲的影響較為敏感，適用于對燃?xì)饨Y(jié)構(gòu)的捕捉。

1.2 輻射傳輸方程

由于液體火箭發(fā)動機燃燒產(chǎn)物無顆粒，所以燃?xì)廨椛鋫鬟f方程可以寫為

式中 λ為某一波長；κλ為光譜吸收系數(shù)；Lλ(r, s )為r位置s方向的光譜輻射亮度；Ib,λ為光譜定向輻射強度。

輻射計算域是以二維對稱流場繞噴管軸向方向旋轉(zhuǎn)一周所劃定的區(qū)域，將計算域進行空間離散和立體角離散，最終劃分為互不重疊的控制體和立體角，控制體和立體角的離散示意如圖1所示。

圖1 空間離散和角度離散示意

利用高斯公式和階梯差分公式將式（5）離散為

其中，

式中 c為光速；h為普朗克常數(shù)，h=6.622 617 6× 10-34J·s；BK為波爾茲曼常數(shù)，BK=1.38×10-23J/K。

2 燃?xì)饬鲌鲇嬎?/h2>

以某軌姿控發(fā)動機為研究對象，發(fā)動機采用MMH/NTO推進劑組合，設(shè)計點壓力為2.0 MPa，混合比為1.65；另外，選定室壓為1.8 MPa、混合比為1.5及室壓為2.1 MPa、混合比為1.76兩種工況進行流場計算。以設(shè)計點工況為例，根據(jù)熱力學(xué)計算平衡組分摩爾百分?jǐn)?shù)計算結(jié)果如表1所示。

表1 平衡組分摩爾百分?jǐn)?shù)

流場入口、出口條件如表2所示。

表2 流場計算入口、出口條件

燃?xì)饬鲌鰯?shù)值計算結(jié)果如圖2至圖5所示（圖中X、Y軸分別表示發(fā)動機軸向方向和徑向方向的相對距離）。

圖2 發(fā)動機燃?xì)饬鲌鰷囟确植?/p>

圖3 發(fā)動機燃?xì)饬鲌鯟O摩爾百分?jǐn)?shù)分布

圖4 發(fā)動機燃?xì)饬鲌鯤2O摩爾百分?jǐn)?shù)分布

圖5 發(fā)動機燃?xì)饬鲌鯟O2摩爾百分?jǐn)?shù)分布

由圖2可知，燃?xì)鈴膰姽車姵龊?，由于出口壓力小于環(huán)境壓力，因此在噴管出口不遠(yuǎn)處出現(xiàn)正激波，溫度上升，當(dāng)燃?xì)鈮毫Ω哂诖髿鈮毫r，形成膨脹波，使得壓力、溫度下降，如此往復(fù)，形成多個膨脹壓縮波系，溫度在每個波系下變化規(guī)律相同，且隨著湍流耗散和激波損失影響，燃?xì)鉁囟戎饾u降低。

圖3至圖5為主要輻射組分CO、H2O、CO2摩爾百分?jǐn)?shù)分布，3種組分分布相似，沿噴管軸向逐漸減??；圖2至圖5計算結(jié)果與文獻(xiàn)[5]流場計算結(jié)果相近，表明燃燒流場計算模型正確，計算結(jié)果可作為燃?xì)饧t外輻射數(shù)值計算流場數(shù)據(jù)。燃?xì)饧t外圖像反應(yīng)示意見圖6。

圖6 燃?xì)饧t外圖像

由圖6可知，燃?xì)饬鲌瞿M結(jié)果與試驗結(jié)果相符。

3 燃?xì)饧t外輻射數(shù)值計算

3.1 燃?xì)馕障禂?shù)計算

燃燒產(chǎn)物中紅外輻射主要組分為H2O、CO2和CO，燃?xì)獾墓庾V吸收系數(shù)為主要組分光譜吸收系數(shù)之和，本文采用逐線法對各氣體光譜吸收系數(shù)進行計算。

高溫下燃?xì)夤庾V譜線強度為

譜線線形選為佛奧特譜線，歸一化后得：

其中，

式中0P為標(biāo)準(zhǔn)狀態(tài)壓力，0P=0.101 325 MPa；Dγ為多普勒半寬；Lγ為碰撞半寬。

具體計算方法采用Humlcek算法的改進算法[6,7](,)K x y進行計算。多普勒半寬為

式中 m為分子質(zhì)量。碰撞半寬為

式中zP為氣體總壓力；sP為組分氣體分壓；n為空氣增寬溫度依賴指數(shù)；airγ為空氣增寬半寬；selfγ為自增寬半寬；n，airγ和selfγ的參數(shù)由分子譜線數(shù)據(jù)庫提供。

通過自編程逐線法對標(biāo)準(zhǔn)狀況下CO2在4.3 μm譜帶吸收系數(shù)進行計算，并將結(jié)果與HΙTRAN-Web[8]計算結(jié)果進行對比，如圖7所示。

圖7 逐線法程序計算CO24.3μm譜帶光譜吸收系數(shù)

由圖7可知，自編程光譜吸收系數(shù)計算結(jié)果與HΙTRAN-Web計算結(jié)果符合較好，說明該計算程序正確，計算精度高。

3.2 燃?xì)饧t外輻射計算結(jié)果和分析

以3種典型的發(fā)動機燃燒工況流場計算結(jié)果為基礎(chǔ)，以溫度為1 000 K為閾值，將流場中燃?xì)夂诵膮^(qū)進行提取，并將該區(qū)域按照圓柱體補齊，立體角劃分為24個，其中天頂角4個，每個天頂角下圓周角分布分別有4、8、8、4個。利用輻射能量傳輸離散方程和共軛梯度穩(wěn)定法[9]得到該燃?xì)夂诵膮^(qū)輻射分布情況如圖8～10所示。

圖8 燃?xì)庠?～10μm波段內(nèi)沿θ=3π/8光譜輻射強度分布

圖9 燃?xì)庠?～10μm波段內(nèi)不同天頂角輻射強度分布

圖10 不同波段輻射強度百分比

由圖8可知，在3種典型發(fā)動機工況下，燃?xì)夤庾V輻射強度分布規(guī)律一致，燃?xì)饬鲌龉庾V輻射強度存在3個輻射強度較大的峰值區(qū)域，第1個峰值區(qū)域為2.5～3.0 μm波段，這是因為2.5～3.0 μm波段存在CO2輻射帶（2.64～2.84 μm）以及H2O輻射帶（2.24～3.37 μm）；第2個峰值區(qū)域為4.2～4.7 μm波段，該波段存在CO2輻射帶（4.01～4.8 μm），在4.5 μm附近存在CO輻射帶（4.3～4.5 μm）；第3個峰值區(qū)域為5.05～5.2 μm波段，在5.15 μm附近存在H2O輻射帶（5.1～5.25 μm）。通過對比發(fā)現(xiàn)，室壓越高，混合比越大，輻射強度越強。這是因為此時燃?xì)夂诵膮^(qū)平均溫度最高，主要輻射組分H2O和CO2含量也最高。本文計算結(jié)果與文獻(xiàn)[10]、文獻(xiàn)[11]計算結(jié)果一致，表明輻射數(shù)值計算方法有效、結(jié)果可靠。

由圖9可知，不同工況下，輻射強度在不同天頂角下2～10 μm波段內(nèi)分布規(guī)律相同，室壓越高，混合比越大，輻射強度越大；在0～180°天頂角范圍內(nèi)，越接近90°，燃?xì)膺吔缤队懊娣e越大，輻射強度越高，相反越接近0°或者180°，燃?xì)膺吔缤队懊娣e越小，輻射強度越小。

由圖10可知，不同工況下各波段輻射強度百分比排列規(guī)律一致，2～5 μm波段最大，占總輻射強度的70%以上；5～7 μm波段次之，占總輻射強度的19%以上；7～10 μm波段最小，約占總輻射強度的9%左右。

4 結(jié) 論

本文采用有限體積法與燃?xì)饬鲌鰯?shù)值計算相結(jié)合的方法對軌姿控發(fā)動機地面典型工況燃?xì)饧t外輻射特性進行計算和分析，為掌握發(fā)動機燃?xì)饧t外輻射特性以及飛行器紅外目標(biāo)識別奠定了基礎(chǔ)。

a）燃?xì)饧t外輻射特性在2～10 μm波段內(nèi)具有光譜選擇性，在2.5～3 μm、4.2～4.7 μm以及5.05～5.2 μm波段輻射較強，主要由H2O、CO2和CO等選擇吸收性強的燃?xì)饨M分造成的。

b）對于相同天頂角，室壓越大，混合比越大，輻射強度越大；對于不同天頂角，各工況輻射強度分布相似，天頂角越接近90°，輻射強度越大，天頂角越接近0°或者180°，輻射強度越小。

c）發(fā)動機典型工況下，各波段輻射強度百分比基本相同，波段越短，所占百分比越大，燃?xì)廨椛鋸姸燃性?～5 μm波段內(nèi)，占總波段輻射強度的70%以上。

[1] Surzhikov S T. Direct simulation monte aarlo algorithms for the rocket exhaust plumes emissivity prediction[R]. AΙAA-2002-0795, 2002.

[2] 阮立明, 齊宏, 王圣剛, 等. 導(dǎo)彈尾噴焰目標(biāo)紅外特性的數(shù)值仿真[J].紅外與激光工程, 2008, 37(6): 959-962.

[3] 張小英, 朱定強, 蔡國飆. 固體火箭羽流紅外特性的DOM 法模擬及高度影響研究[J]. 宇航學(xué)報, 2007, 28(3): 702-706.

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

[5] 王偉臣, 李世鵬, 張嶠, 王寧飛 ,何丹薇. 火箭發(fā)動機羽流紅外特性計算方法研究[J]. 推進技術(shù), 2010, 31(4): 423-427.

[6] Kuntz M. A new implementation of Humlicek algorithm for the calculation of the voigt profile function[J]. Journal of Quantit Spectrosc and Radiative Transfer, 1997, 57(6): 819-824.

[7] Ruyten W. A new implementation of Humlicek algorithm for the calculation of the voigt profile function[J]. Journal of Quantit Spectrosc and Radiative Transfer, 2004, 86(2): 231-233.

[8] HΙTRAN on the Web[EB/OL]. (2014-03-15)[2014-03-15]. http:// hitran. iao. ru/home.

[9] Ferziger H J, Peric M. Computational methods for fluid dynamics[M]. Berlin: Springer, 1996.

[10] 劉尊洋, 邵立, 汪亞夫, 孫曉泉. 飛行參數(shù)對液體火箭發(fā)動機尾焰紅外輻射特性的影響[J]. 光學(xué)學(xué)報, 2013, 33(4): 1-7.

[11] 金偉, 呂相銀. 衛(wèi)星姿控推進器尾焰紅外輻射的數(shù)值模擬[J]. 紅外與激光工程, 2011, 40(4): 595-599.

Simulation on Liquid Rocket Engine Gas Infrared Characteristic Base on Finite Volume Method

Wang Da-rui, Liang Wei, Ma Li-na, Yao Shi-qiang
(Beijing Aerospace Propulsion Ιnstitute, Beijing, 100076)

Using method of finite volume, we simulate infrared radiation properties of the attitude and orbit control engine gas flow field. First for the engine gas flow field is calculated by Fluent in typical working conditions. Used the latest spectrum database for the major components of radiation absorption coefficient to calculate. Finally, the gas infrared radiation characteristic is calculated. Analyzed infrared radiation intensity distribution in different wavelengths, different zenith angle, different wave band. The results show that the proposed method of calculation is reliable and calculation results of infrared characteristics is accurate, conducive to liquid rocket engineers further grasp the infrared radiation characteristics calculation and analysis methods.

Attitude and orbit control engine; Ιnfrared radiation; Simulation calculation; Line by line method

V439+.7

A

1004-7182(2017)04-0039-05

DOΙ：10.7654/j.issn.1004-7182.20170410

2016-11-30；

2017-03-04

王大銳（1986-），男，博士，工程師，主要研究方向為液體火箭發(fā)動機推力室設(shè)計及發(fā)動機狀態(tài)監(jiān)測技術(shù)