王偉臣，魏志軍，張 嶠，王寧飛，王 宏，孫 美

(1.北京理工大學 宇航學院，北京 100081;2.西安近代化學研究所，西安 710065))

0 引言

固體發(fā)動機羽流紅外輻射特性是飛行器的重要特征信號，研究羽流紅外輻射機理對于飛行器底部部件設計，探測系統(tǒng)設計和隱身性能研究，推進劑配方改進等都具有重要意義。固體推進劑在燃燒過程中產生貧氧富燃的燃氣，噴入大氣后，逐漸從周圍空氣中卷吸一部分氧氣進行二次燃燒(后燃)，這會使得羽流近場區(qū)域的溫度有較大升高，極大地加強了羽流的輻射強度。對于含鋁粉的固體推進劑，羽流中含有較多的Al2O3顆粒，顆粒相具有較強的發(fā)射能力，是羽流紅外特性的重要因素。將流場計算和輻射傳輸過程耦合求解以顯著提高輻射計算的精度，在以往研究中，流場計算和輻射傳輸計算大多是獨立進行。董士奎［1］、樊士偉和張小英等［2］使用半經驗法計算了流場物性參數，然后對羽流場幾個特定波段內的輻射特性進行了計算。詹光［3］和豐松江等［4］使用 CFD方法計算得到了較為準確的羽流場，然后求解了1～5 μm波段內羽流紅外輻射信號，但沒有將輻射傳輸與流場計算進行耦合求解。豐松江［5］等在計算得到的氣-固兩相流場的基礎上，考慮Al2O3顆粒輻射作用，使用有限增量法對固體推進劑羽流紅外特性進行了計算。在發(fā)表的文獻中，大多是對單個觀測點處接收到來自羽流的輻射能量進行研究，對輻射強度在羽流近場內各點處的分布情況研究較少，沒有得到較細致的輻射強度在羽流場內的分布圖。

本文建立了固體推進劑羽流紅外特性計算模型，通過在流場能量方程中引入輻射源項，實現(xiàn)了流場計算與輻射傳輸的耦合求解。使用FLUENT軟件進行計算，基于歐拉-拉格朗日方法對考慮后燃化學反應的氣-固兩相羽流流場進行了計算，使用有限速率化學反應模型模擬羽流后燃現(xiàn)象，使用離散顆粒模型DPM對羽流中顆粒運動進行了跟蹤計算。采用離散坐標法求解羽流輻射傳輸方程，得到了羽流紅外輻射強度在1 000～4 500 cm－1內的分布情況及輻射強度在近場內的分布云圖，通過與試驗數據進行對比，對耦合求解方法的可靠性和計算精度進行了驗證;就不同鋁粉含量對推進劑羽流紅外特性的影響程度進行了研究。

1 物理模型

1.1 羽流流場能量耦合模型

通過在能量方程中引入顆粒源項和化學反應源項，可在流場計算中考慮顆粒運動和化學反應的影響;通過在能量方程中引入輻射能量源項，可實現(xiàn)流場與輻射換熱的耦合計算，獲得較高輻射傳輸的求解精度。由此得到流場能量守恒方程為

方程右邊后3項分別為由于導熱、組分擴散和粘性耗散所引起的能量傳遞。Sp為顆粒能量源項;Sh為化學源項，代表化學反應引起的能量增量;Sr是輻射源項，代表輻射發(fā)射與傳輸引起的能量增量。

1.2 顆粒運動模型

Al2O3顆粒在羽流中所占的體積分數較小，與氣相之間存在熱量和動量交換，不考慮顆粒的蒸發(fā)和質量變化，顆粒與流體相互作用時忽略重力、熱泳力、Brownian力和Saffman升力，顆粒受到湍流作用時采用隨機游走模型。顆粒運動方程為

式中 up是顆粒速度為氣體拖曳力，ρp為顆粒密度，dp為顆粒直徑，Re為相對雷諾數，CD為流體拖曳系數［6］。

1.3 有限速率化學反應模型

羽流中的化學反應主要發(fā)生在近場區(qū)域，該區(qū)域基本為超音速流，使用Arrhenius定律描述后燃現(xiàn)象中詳細的化學反應機理［7］。用Arrhenius定律表示的正向化學反應速率常數如式(3)所示:

式中 Ar是指前因子;n是溫度指數;Er是活化能;R是通用氣體常量;Ar與kf的量綱相同，為cm3/(mol/s)。

1.4 輻射傳輸模型

對于含顆粒相的羽流流場，輻射傳輸方程為

式中→r為位置向量;→s為方向向量;a為氣體吸收系數;σp為顆粒散射系數;n為氣體折射率;σ為黑體輻射常數代表單元體內所有顆粒發(fā)射的輻射能量，εp，n為單個顆粒的發(fā)射率，Tp，n為單個顆粒的溫度;I(r→，s→)為位置 r→處的微元體在 s→ 方向上的輻射強度;s→＇為散射方向向量;Φ(s→，s→＇)為散射相函數;Ω＇為空間立體角。

2 數值計算方法

2.1 流場計算方法

對于連續(xù)氣相流場，使用歐拉方法對其進行計算;對于離散顆粒的運動軌跡，使用拉格朗日方法進行計算。進行流場計算時，采用有限體積法對流場控制方程進行離散，使用基于密度算法的隱式求解器進行求解，湍流模型選用RNG k-ε模型［8］，壁面附近采用標準壁面函數。

2.2 輻射計算方法

由于輻射傳輸具有方向性，使用天頂角θ和方位角φ對每個微元體內的空間傳輸方向進行離散，在每個傳播方向上采用二階迎風格式的離散坐標法［9］對輻射傳輸方程進行離散和求解。

利用高分辨率氣體光譜數據庫HITRAN 2004［10］，使用數值平均方法［11］計算每種組分氣體在各個波帶內的吸收系數;使用 Lorentz-Mie散射理論［12］計算Al2O3顆粒的散射系數σp和發(fā)射率εp。

3 羽流紅外特性測試

為研究固體推進劑中鋁粉含量對羽流紅外特性的影響，本文選取了2種配方的改性雙基推進劑，配方二在配方一的基礎上添加了5%的鋁粉，使用紅外光譜儀(FTIR)對羽流的紅外特性進行測試。圖1為羽流紅外測試儀器放置示意圖。FTIR距發(fā)動機軸線距離為4.3 m，測試對象為羽流軸線上1個直徑為60 mm的光斑，該光斑距發(fā)動機噴口0.5 m。詳細測試方法及原理見文獻［13］。本次試驗對該光斑在1 000～4 500 cm－1范圍內的紅外輻射強度進行測試，得到紅外強度的光譜分布。根據CO2、H2O和CO等主要氣體的發(fā)射峰的分布特性，將1 000～4 500 cm－1范圍劃分為 7個波數區(qū)間(1 000～1 600 cm－1，1 600 ～1 950 cm－1，1 950 ～2 150 cm－1，2 150 ～2 400 cm－1，2 400 ～ 3 100 cm－1，3 100 ～ 3 800 cm－1，3 800～4 500 cm－1)，將紅外強度譜圖在每個區(qū)間內進行光譜積分，得到紅外強度在每個區(qū)間內的分布情況。

圖1 羽流紅外輻射測試儀器放置示意圖Fig.1 Schematic diagram of infrared radiation test device of exhaust plume

4 計算模型及結果分析

4.1 計算模型

本文使用二維軸對稱計算區(qū)域對羽流流場和輻射傳輸過程進行計算，噴管出口直徑為23 mm，噴喉直徑為10 mm，收縮半角為45°，擴張半角為15°，羽流計算區(qū)域長10 m，高1 m，網格節(jié)點劃分為300×130。在噴管軸線和壁面附近采用較密的網格，遠離軸線的區(qū)域由于流場參數變化不大而采用稀疏網格。整個計算域內網格總數為4.4×104，通過對多次計算結果進行對比表明，此網格劃分策略具備較好的無關性。計算域輪廓及噴管出口處網格劃分情況如圖2所示。

圖2 計算域輪廓及噴管出口處網格劃分Fig.2 Computation region and the grid around the nozzle exit

噴管入口邊界采用壓力入口，配方一和配方二的燃燒室的實測平均壓強分別為7 MPa和7.4 MPa，燃氣溫度分別為2 646 K和2 884 K。羽流外邊界條件為壓力出口，環(huán)境壓力為101.325 kPa，溫度為300 K。忽略外部熱源的輻射加熱作用。

計算中考慮9種主要氣體組分:H2O、CO、CO2、H2、N2、O2、OH、H 和 O，使用最小自由能法［14］對 2 種推進劑進行熱力學計算，得到2種配方推進劑的燃氣組分在噴管入口處的含量如表1所示。各組分的熱導率kj和比定壓熱容cp，j由文獻［15］中提供的數據進行插值得到。

羽流中的后燃現(xiàn)象主要由H2和CO的氧化反應組成，本文使用H2/CO氧化反應體系［16］進行化學計算，所用反應機理數據如表2所示。

配方二推進劑羽流中含有鋁粉，其燃燒產物為Al2O3顆粒。在噴管入口處加入 Al2O3顆粒，使用Reed等人推薦的Hermsen公式［17］計算顆粒平均粒徑為2.02 μm，顆粒粒徑分布采用Rosin-Rammler分布規(guī)律。

為研究羽流流場在7個波數區(qū)間內的紅外特性，分別對每個波數區(qū)間內的羽流的流場及紅外特性進行計算。通過對比發(fā)現(xiàn)，二維算例和三維算例的計算結果差別不大，為節(jié)約計算時間，本文均采用二維算例進行計算。

表1 各燃氣組分在燃燒室中的含量Table 1 Mass fractions of species in combustion chamber of two kind of propellants

表2 羽流中化學反應模型Table 2 Chemistry reaction model in exhaust plume

4.2 計算結果分析

圖3為配方一和配方二推進劑羽流的近場溫度云圖，圖4為2種配方推進劑羽流軸線上的溫度分布對比。由圖3可見，羽流核心區(qū)的中后部溫度較高，這是由于富燃燃氣與空氣中的氧氣在此混合，發(fā)生了后燃反應。由圖4可知，配方二羽流軸線附近的溫度高于配方一，這可能是由于配方二推進劑的能量特性較高，羽流中后燃程度較為劇烈;配方二羽流中的Al2O3顆粒具有一定的溫度滯后，這會進一步提高流場溫度。

羽流中輻射能力較強的氣體組分為H2O和CO2，這2種氣體的質量分數亦較大，限于篇幅，只給出配方二推進劑羽流中這2種氣體組分的分布情況。圖5所示為羽流近場H2O和CO2質量分數云圖。

圖3 配方一和配方二推進劑羽流近場溫度云圖Fig.3 Temperature contour of plume near field from case 1 and case 2

圖4 配方一和配方二推進劑羽流軸線上溫度分布對比Fig.4 Comparison of temperature distribution in axis line of plume from case 1 and case 2

圖5 羽流近場H2O和CO2組分質量分數分布云圖Fig.5 Mass fractions contour of H2O and CO2 in plume near field

由圖5可見，此2種組分氣體的分布隨流場波系結構而變化，在軸線附近分布較多，在羽流混合層內分布較少。對于配方二推進劑，Al2O3顆粒是羽流中的重要組分。圖6為直徑2 μm的Al2O3顆粒的運動軌跡圖。圖6表明，多數粒子離開噴管出口后在羽流核心區(qū)內運動，少數粒子在運動一定距離后偏離羽流核心區(qū)域。

圖6 直徑2 μm的顆粒運動軌跡Fig.6 Motion trajectory of 2 μm particles

圖7為配方一推進劑羽流軸線上x=0.5 m處紅外強度的計算結果與試驗數據的對比情況。為與試驗數據進行對比，將計算所得的各點輻射強度在0.5 m處光斑區(qū)域內進行了積分平均計算，得到該光斑處的輻射強度。由圖7可見，觀測點上的紅外強度在2 150～2 400 cm－1區(qū)間內最大，主要由 CO2和 CO的發(fā)射引起的;在1 000～1 600 cm－1區(qū)間內，紅外強度亦較為明顯，主要由 H2O的長波發(fā)射導致的;在3 800～4 500 cm－1區(qū)間內，H2O和CO有少量發(fā)射譜線，試驗結果表明這2種氣體在該區(qū)間的發(fā)射能力較弱;在2 400～3 100 cm－1區(qū)間內，輻射強度極其微弱，是因為該區(qū)間內無明顯的輻射組分。本文計算結果與試驗數據的平均絕對誤差為21%，吻合較好。

圖7 配方一推進劑羽流軸線上觀測點處各區(qū)間內輻射強度的分布情況Fig.7 Spectral radiation intensity of plume at monitor point on axis line from case 1

圖8為配方二推進劑羽流軸線上x=0.5 m處紅外強度的計算結果與試驗數據的對比情況。由圖8可見，配方二推進劑羽流的輻射強度在不同波數區(qū)間內的變化趨勢與配方一基本一致，各個區(qū)間內的輻射強度比配方一均有大幅提高。配方二推進劑中增加的5%鋁粉會提高燃燒室內燃氣溫度，從而提高羽流流場的溫度;配方二推進劑羽流中含有的Al2O3顆粒具有較強的發(fā)射能力，會提高羽流的整體發(fā)射能力。這2個因素可能是推進劑中引入鋁粉后，羽流紅外強度增加的主要原因。3 100～3 800 cm－1區(qū)間內的輻射強度相對上升，1 600～1 950 cm－1和1 950 ～2 150 cm－1區(qū)間內的輻射強度相對下降。本文計算結果與試驗數據的平均絕對誤差為20%，吻合較好。

圖8 配方二推進劑羽流軸線上觀測點處各區(qū)間內輻射強度的分布情況Fig.8 Spectral radiation intensity of plume at monitor point on axis line from case 2

圖9為配方一和配方二推進劑羽流近場在2 150～2 400 cm－1區(qū)間內輻射強度分布云圖。

圖9 配方一和配方二推進劑羽流近場在2 150～2 400 cm－1區(qū)間內輻射強度分布圖云圖Fig.9 Radiation intensity contour of plume in 2 150～2 400 cm band near field－1 from case 1 and case 2

由圖9可見，輻射強度變化與流場溫度的變化基本一致，輻射強度最高值出現(xiàn)在羽流核心區(qū)軸線上的中后部。這是由于該區(qū)域內化學反應較為充分，導致流場溫度較高，氣體輻射能力較強。在羽流軸線方向上，紅外強度較高，傳播較遠;在垂直于軸線的方向上，紅外強度衰減較快。相對于以往的計算方法，將流場和輻射傳輸進行耦合求解能夠捕捉到輻射特性隨流場物性變化的情況，表明此種方法能提高輻射計算結果的精度。

4.3 鋁粉含量對推進劑羽流紅外特性的影響

由圖7和圖8可知，2 150～2 400 cm－1區(qū)間內的輻射能量占羽流紅外輻射的主要部分，為研究推進劑中不同鋁粉含量對羽流紅外特性的影響，對鋁粉含量分別為10%、15%、20%的推進劑羽流在2 150～2 400 cm－1區(qū)間內的紅外特性進行了計算。

圖10和圖11分別為不同鋁粉含量的推進劑羽流軸線上溫度和2 150～2 400 cm－1區(qū)間內輻射強度分布的對比情況。

圖10 不同鋁粉含量的推進劑羽流軸線上溫度分布對比Fig.10 Comparison of temperature distribution of propellants with different aluminum in axis line

圖11 不同鋁粉含量的推進劑羽流軸線上2 150～2 400 cm－1區(qū)間內輻射強度分布對比Fig.11 Comparison of radiation intensity distribution of propellants with different aluminum in 2 150～2 400 cm－1in axis line

由圖10可見，推進劑中鋁粉含量增加會提高其能量特性，燃燒室內燃氣溫度上升，羽流流場的溫度也隨之上升。由圖11可見，羽流輻射強度隨推進劑中鋁粉含量增加而顯著上升，輻射強度峰值出現(xiàn)在0.7 m位置附近。0.5 m后軸線上輻射強度呈現(xiàn)鋸齒形波動，這可能是離散的Al2O3顆粒進入流場后導致相鄰流體微元的輻射呈現(xiàn)一定的不連續(xù)性。結果表明，Al2O3顆粒會顯著提高推進劑羽流的紅外輻射能力。

5 結論

(1)建立了固體推進劑羽流紅外傳輸的計算模型，基于歐拉-拉格朗日方法對考慮后燃化學反應的氣-固兩相羽流流場進行了計算，通過在流場能量方程中引入輻射源項，實現(xiàn)了流場計算與輻射傳輸的耦合求解。

(2)對推進劑羽流紅外強度在1 000～4 500 cm－1范圍內的分布情況進行了試驗和數值研究，得到了紅外輻射強度在典型區(qū)間內的變化曲線和在羽流近場內的分布云圖，計算結果與試驗數據吻合較好，證明了計算模型的可靠性。

(3)在羽流近場內能捕捉到輻射特性隨流場物性變化的情況，表明耦合求解能提高輻射計算結果的精度。

(4)對不同鋁粉含量對推進劑羽流紅外特性的影響程度進行了研究，Al2O3顆粒會顯著提高推進劑羽流的紅外輻射能力，羽流輻射強度隨推進劑中鋁粉含量增加而大幅度上升。

［1］董士奎，于建國，李東輝，等.貼體坐標系下離散坐標法計算尾噴焰輻射特性［J］上海理工大學學報，2003，25(2):159-162.

［2］樊士偉，張小英，朱定強，等.用FVM法計算固體火箭羽流的紅外特性［J］.宇航學報，2005，26(6):793-797.

［3］詹光，李椿萱.發(fā)動機燃氣噴流紅外輻射場的數值模擬［J］.北京航空航天大學學報，2005，31(8):829-833.

［4］Feng S，Nie W，Xie Q，et al.Numerical simulation of flow field and radiation of an aluminized solid propellant rocket multiphase exhaust plume［R］.AIAA 2007-4415，2007.

［5］豐松江，聶萬勝，宋豐華，等.固體發(fā)動機尾焰紅外輻射特性預估研究［J］.固體火箭技術，2009，32(2):183-187.

［6］Morsi S A，Alexander A J.An investigation of particle trajectories in two-phase flow systems［J］.J.Fluid Mech.，1972，55(2):193-208.

［7］苗瑞生.發(fā)射氣體動力學［M］.北京:國防工業(yè)出版社，2006.

［8］Yakhot V，Orszag S A.Renormalization group analysis of turbulence:I.Basic theory［J］.Journal of Scientific Computing，1986，1(1):1-51.

［9］談和平，夏新林，劉林華，等.紅外輻射特性與傳輸的數值計算［M］.哈爾濱:哈爾濱工業(yè)大學出版社，2006.

［10］Rothman L S，Jacquemart D，Barbe A，et al.The HITRAN 2004 molecular spectroscopic database［J］.Journal of Quantitative Spectroscopy ＆ Radiative Transfer，2005，96(2):139-204.

［11］Young S J.Nonisothermal band model theory［J］.Journal of Quantitative Spectroscopy and Radiative Transfer，1976，18:1-28.

［12］張小英，朱定強，蔡國飆.固體火箭羽流中Al2O3粒子的輻射特性［J］.固體火箭技術，2006，29(4):247-250.

［13］張小玲，李春迎，王宏，等.FTIR光譜遙測固體火箭發(fā)動機推進劑燃燒輻射能［J］.廣西師范大學學報(自然科學版)，21(2):205-207.

［14］張平，孫維申，眭英.固體火箭發(fā)動機原理［M］.北京:北京理工大學出版社，1992:75-129.

［15］鄭亞，陳軍，鞠玉濤，等.固體火箭發(fā)動機傳熱學［M］.北京:北京航空航天大學出版社，2006:256-260.

［16］Jensen D E，Jones G A.Reaction rate coefficients for flame calculations［J］.Combustion and Flame，1978，32:1-34.

［17］Reed R A，Calia V S.Review of aluminum oxide rocket exhaust particles［R］.AIAA 93-2819，1993.