基于有限體積法的液體火箭發動機燃氣紅外特性仿真

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

基于有限體積法的液體火箭發動機燃氣紅外特性仿真

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

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

運用有限體積法對某軌姿控發動機典型工況燃氣流場紅外輻射特性進行仿真計算。首先通過Fluent軟件對發動機典型工況燃氣流場進行計算，再利用最新的光譜數據庫對主要輻射組分吸收系數進行計算，最后對燃氣紅外輻射特性進行計算；重點分析不同波長、不同天頂角、不同波段下燃氣紅外輻射強度分布規律。研究結果表明：所提計算方法可靠，紅外特性計算結果準確，有利于掌握液體火箭發動機紅外輻射特性計算和分析方法。

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

0 引 言

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

本文采用有限體積法與燃氣流場數值計算相結合的方法對軌姿控發動機燃氣流場紅外輻射特性進行了計算和分析。通過Fluent軟件對發動機典型工況下軸對稱二維網格燃氣流場進行計算。基于HΙTRAN2012和HΙTEMP2010光譜數據庫，通過自編逐線法程序對H2O、CO2、CO主要輻射組分吸收系數進行了計算。最后利用流場計算溫度、壓力、吸收組分分布以及各組分吸收系數等數據，采用燃氣輻射計算程序對不同工況下的燃氣流場紅外輻射特性進行計算，并分析了不同波長、不同天頂角、不同波段下燃氣紅外輻射分布規律。

1 流場計算模型

1.1 流場控制方程

質量守恒方程為

式中 ρ為混合氣體密度；iY為組分i的質量分數；v為速度矢量；iD為組分i的擴散系數。

動量守恒方程為

式中 P為流體微元體上的靜壓力；τ為微元體上粘性應力張量，表達式為

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

能量守恒方程為

湍流模型采用RNGkε?（Renormalized Groupkε?）模型，該模型對于瞬變流和流線彎曲的影響較為敏感，適用于對燃氣結構的捕捉。

1.2 輻射傳輸方程

由于液體火箭發動機燃燒產物無顆粒，所以燃氣輻射傳遞方程可以寫為

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

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

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

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

其中，

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

2 燃氣流場計算

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

表1 平衡組分摩爾百分數

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

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

燃氣流場數值計算結果如圖2至圖5所示（圖中X、Y軸分別表示發動機軸向方向和徑向方向的相對距離）。

圖2 發動機燃氣流場溫度分布

圖3 發動機燃氣流場CO摩爾百分數分布

圖4 發動機燃氣流場H2O摩爾百分數分布

圖5 發動機燃氣流場CO2摩爾百分數分布

由圖2可知，燃氣從噴管噴出后，由于出口壓力小于環境壓力，因此在噴管出口不遠處出現正激波，溫度上升，當燃氣壓力高于大氣壓力時，形成膨脹波，使得壓力、溫度下降，如此往復，形成多個膨脹壓縮波系，溫度在每個波系下變化規律相同，且隨著湍流耗散和激波損失影響，燃氣溫度逐漸降低。

圖3至圖5為主要輻射組分CO、H2O、CO2摩爾百分數分布，3種組分分布相似，沿噴管軸向逐漸減小；圖2至圖5計算結果與文獻[5]流場計算結果相近，表明燃燒流場計算模型正確，計算結果可作為燃氣紅外輻射數值計算流場數據。燃氣紅外圖像反應示意見圖6。

圖6 燃氣紅外圖像

由圖6可知，燃氣流場模擬結果與試驗結果相符。

3 燃氣紅外輻射數值計算

3.1 燃氣吸收系數計算

燃燒產物中紅外輻射主要組分為H2O、CO2和CO，燃氣的光譜吸收系數為主要組分光譜吸收系數之和，本文采用逐線法對各氣體光譜吸收系數進行計算。

高溫下燃氣光譜譜線強度為

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

其中，

式中0P為標準狀態壓力，0P=0.101 325 MPa；Dγ為多普勒半寬；Lγ為碰撞半寬。

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

式中 m為分子質量。碰撞半寬為

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

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

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

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

3.2 燃氣紅外輻射計算結果和分析

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

圖8 燃氣在2～10μm波段內沿θ=3π/8光譜輻射強度分布

圖9 燃氣在2～10μm波段內不同天頂角輻射強度分布

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

由圖8可知，在3種典型發動機工況下，燃氣光譜輻射強度分布規律一致，燃氣流場光譜輻射強度存在3個輻射強度較大的峰值區域，第1個峰值區域為2.5～3.0 μm波段，這是因為2.5～3.0 μm波段存在CO2輻射帶（2.64～2.84 μm）以及H2O輻射帶（2.24～3.37 μm）；第2個峰值區域為4.2～4.7 μm波段，該波段存在CO2輻射帶（4.01～4.8 μm），在4.5 μm附近存在CO輻射帶（4.3～4.5 μm）；第3個峰值區域為5.05～5.2 μm波段，在5.15 μm附近存在H2O輻射帶（5.1～5.25 μm）。通過對比發現，室壓越高，混合比越大，輻射強度越強。這是因為此時燃氣核心區平均溫度最高，主要輻射組分H2O和CO2含量也最高。本文計算結果與文獻[10]、文獻[11]計算結果一致，表明輻射數值計算方法有效、結果可靠。

由圖9可知，不同工況下，輻射強度在不同天頂角下2～10 μm波段內分布規律相同，室壓越高，混合比越大，輻射強度越大；在0～180°天頂角范圍內，越接近90°，燃氣邊界投影面積越大，輻射強度越高，相反越接近0°或者180°，燃氣邊界投影面積越小，輻射強度越小。

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

4 結 論

本文采用有限體積法與燃氣流場數值計算相結合的方法對軌姿控發動機地面典型工況燃氣紅外輻射特性進行計算和分析，為掌握發動機燃氣紅外輻射特性以及飛行器紅外目標識別奠定了基礎。

a）燃氣紅外輻射特性在2～10 μm波段內具有光譜選擇性，在2.5～3 μm、4.2～4.7 μm以及5.05～5.2 μm波段輻射較強，主要由H2O、CO2和CO等選擇吸收性強的燃氣組分造成的。

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

c）發動機典型工況下，各波段輻射強度百分比基本相同，波段越短，所占百分比越大，燃氣輻射強度集中在2～5 μm波段內，占總波段輻射強度的70%以上。

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

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

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

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

[5] 王偉臣, 李世鵬, 張嶠, 王寧飛 ,何丹薇. 火箭發動機羽流紅外特性計算方法研究[J]. 推進技術, 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] 劉尊洋, 邵立, 汪亞夫, 孫曉泉. 飛行參數對液體火箭發動機尾焰紅外輻射特性的影響[J]. 光學學報, 2013, 33(4): 1-7.

[11] 金偉, 呂相銀. 衛星姿控推進器尾焰紅外輻射的數值模擬[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-），男，博士，工程師，主要研究方向為液體火箭發動機推力室設計及發動機狀態監測技術