矩形噴管外尾焰紅外輻射特性的數(shù)值計(jì)算

馮云松，李曉霞，路遠(yuǎn)，金偉

(解放軍電子工程學(xué)院 安徽省紅外與低溫等離子體重點(diǎn)實(shí)驗(yàn)室，安徽合肥230037)

0 引言

飛機(jī)發(fā)動(dòng)機(jī)的排氣系統(tǒng)是飛機(jī)的主要紅外輻射源，采用矩形噴管是一種有效抑制排氣系統(tǒng)紅外輻射的手段。由于矩形噴管在實(shí)現(xiàn)隱身性上的優(yōu)勢(shì)，轟炸機(jī)、戰(zhàn)斗機(jī)以及無(wú)人機(jī)都廣泛采用了矩形噴管[1]。因此，對(duì)矩形噴管外尾焰紅外輻射特性研究頗有必要。

目前，國(guó)外在飛機(jī)排氣系統(tǒng)紅外輻射特性研究方面的理論、技術(shù)和方法都進(jìn)行了研究，如Varney[2]，Jim[3]對(duì)飛機(jī)排氣系統(tǒng)的紅外輻射特性進(jìn)行了較詳細(xì)的數(shù)學(xué)建模和數(shù)值計(jì)算。國(guó)內(nèi)的學(xué)者從20 世紀(jì)90 年代開(kāi)始已采用多種方法對(duì)火箭、導(dǎo)彈和飛機(jī)等發(fā)動(dòng)機(jī)尾焰的紅外輻射特性進(jìn)行了計(jì)算。如董士奎等[4]用貼體坐標(biāo)系下的離散坐標(biāo)法(DOM)研究了不同工況下尾噴焰的紅外光譜輻射特性;樊士偉等[5]用有限體積法(FVM)計(jì)算了固體火箭發(fā)動(dòng)機(jī)羽流的紅外特性;帥永等[6]用反向蒙特卡羅法(BMCM)模擬了高溫含粒子自由流的紅外輻射特性;阮立明等[7]提出了用源項(xiàng)六流法(SSFM)求解導(dǎo)彈尾噴焰的紅外輻射特性。綜合比較針對(duì)尾焰紅外輻射的多種算法，其中有限體積法具有原理簡(jiǎn)單、計(jì)算精度高和可以求解任意方向的輻射特性等優(yōu)點(diǎn)，因此，本文針對(duì)矩形噴管外尾焰，采用有限體積法數(shù)值計(jì)算其紅外輻射特性。

1 噴管外尾焰流場(chǎng)的計(jì)算

1.1 噴管幾何模型

矩形噴管的幾何模型如圖1 所示，噴管入口取在XOY 面上，Z 軸正向?yàn)樯淞鞣较颉＞匦螄姽艿膶抴=0.9 m，高h(yuǎn) =0.45 m，長(zhǎng)l =0.5 m[1]，噴管出口的當(dāng)量直徑D*=(4wh/π)1/2=0.718 m.

圖1 噴管模型Fig.1 Nozzle model

1.2 數(shù)學(xué)模型

為了簡(jiǎn)化尾焰的流動(dòng)模擬和分析，假設(shè):

1)噴管燃?xì)饬鲃?dòng)為定常流動(dòng);

2)忽略燃?xì)庠谂蛎涍^(guò)程中組分的變化，并認(rèn)為燃?xì)獾亩▔罕葻崾浅Ａ?

3)假設(shè)噴管中燃?xì)馀蛎浟鲃?dòng)過(guò)程是理想等熵流動(dòng)過(guò)程，忽略燃?xì)鈱?duì)噴管壁的傳熱;

4)噴管內(nèi)流動(dòng)為純氣相流動(dòng)，燃?xì)鉃槔硐霘怏w，服從理想氣體狀態(tài)方程;

5)不考慮輻射作用，忽略重力的影響。

根據(jù)以上假設(shè)，包括連續(xù)方程、動(dòng)量方程、能量方程的通用形式[8]為

式中:ρ 為氣體密度;v 為氣體速度矢量;φ 為通用求解變量;Γφ為有效擴(kuò)散系數(shù);Sφ為廣義源項(xiàng)。

湍流采用RNGκ－ε 模型，湍流動(dòng)能κ 方程和湍流耗散率ε 方程[8]分別為

式中:ueff為流體在i 方向上的有效速度;Gκ為由層流速度梯度而產(chǎn)生的湍流動(dòng)能;Gb為由浮力而產(chǎn)生的湍流動(dòng)能;YM為在可壓縮湍流中，過(guò)渡的擴(kuò)散產(chǎn)生的波動(dòng);C1z、C2z、C3z為常量;ακ和αε分別為湍流動(dòng)能κ 和耗散率ε 對(duì)應(yīng)的湍流Prandtl 數(shù)。

1.3 尾焰流場(chǎng)數(shù)值計(jì)算

流場(chǎng)計(jì)算區(qū)域在Z 方向的長(zhǎng)度為14D*，在X和Y 方向的長(zhǎng)度均為5D*.圖2 給出了流場(chǎng)計(jì)算區(qū)域的整體示意圖，噴管進(jìn)口位于ABCD 平面上，出口位于EFGHIJKL 組合面上，區(qū)域MNOPQRST 為整個(gè)計(jì)算區(qū)域邊界。計(jì)算中采用六面體結(jié)構(gòu)化網(wǎng)格，流場(chǎng)計(jì)算區(qū)域網(wǎng)格數(shù)目均在13 萬(wàn)左右。

假定噴管入口氣體為完全燃燒的燃?xì)?，基于燃燒理論，確定主流進(jìn)口處N2、CO2和H2O 蒸氣質(zhì)量百分比分別為0.706，0.209 和0.085;外場(chǎng)邊界和進(jìn)口中引射的氣體均為空氣，成分為N2和O2，其質(zhì)量百分比為0.756 和0.244[9].

設(shè)置飛機(jī)的飛行高度為3 km，飛行馬赫數(shù)為1.5，發(fā)動(dòng)機(jī)工作狀態(tài)為非加力狀態(tài)[1]。入口邊界條件:噴管入口為壓力入口，給定的總溫度為800 K，總壓力為0.16 MPa，入口流動(dòng)角θ = 0°[9－10];計(jì)算IJKLMNOP 區(qū)域的邊界條件:設(shè)定該長(zhǎng)方體的幾個(gè)面為壓力出口，溫度為280 K，總壓力為71 kPa;壁面邊界條件:采用無(wú)滑移固壁邊界條件，壁面設(shè)定為流、固耦合面，溫度場(chǎng)計(jì)算時(shí)不考慮壁面之間的輻射傳熱。

圖2 流場(chǎng)計(jì)算區(qū)域Fig.2 Computational area of flow field

計(jì)算中，設(shè)外界氣流靜止，尾噴管壁面絕熱，忽略尾噴管厚度。流場(chǎng)計(jì)算使用Fluent 6.3 軟件進(jìn)行，采用Couple 隱式算法，壁面附近采用標(biāo)準(zhǔn)壁面函數(shù)進(jìn)行修正，收斂精度為10－4.當(dāng)Fluent 對(duì)尾焰流場(chǎng)的迭代計(jì)算收斂時(shí)，獲得尾焰的靜壓力p 和靜溫度T 云圖，如圖3 所示。

圖3 Y=0 m 處XOZ 切面的尾焰靜壓力和靜溫度分布Fig.3 Distribution of static pressure and static temperature of exhaust plume in XOZ face in Y=0 m

2 氣體的輻射傳遞方程

結(jié)合介質(zhì)的發(fā)射、吸收和散射的相互關(guān)系，文獻(xiàn)[11 －12]建立了求解光譜輻射亮度的輻射傳遞方程，其表達(dá)式如下:

式中:λ 為窄帶l 的中心波長(zhǎng);r 為位置向量;s 為方向向量為氣體平均吸收系數(shù);σsλ(s)為氣體散射系數(shù);Lλ(s，r)為位置r 處微元體在s 方向上的光譜輻射亮度;Lbλ(s)為微元體在s 方向上的黑體光譜輻射亮度;s'為散射方向向量;Φλ(s'，s)為散射相函數(shù);Ωi為空間立體角。

航空發(fā)動(dòng)機(jī)在非加力狀態(tài)下工作，且燃燒較充分時(shí)，其尾焰幾乎不含固體顆粒，因此，可不考慮散射對(duì)輻射的影響。由此，輻射傳遞方程簡(jiǎn)化為

根據(jù)(5)式，若要求解輻射亮度，必須首先獲得傳輸介質(zhì)的平均吸收系數(shù)為了求解平均吸收系數(shù)αλ(s)，先采用洛倫茲線型的統(tǒng)計(jì)窄帶模型(6)式計(jì)算某個(gè)窄帶l 內(nèi)的平均透過(guò)率[13]

式中:ΔC、k 和d 分別為碰撞增寬的半寬、平均吸收因子和譜線間距;u 為光學(xué)厚度。

由于尾焰中含有多種吸收氣體，可以先求出每種吸收氣體的吸收系數(shù)，然后對(duì)所有吸收系數(shù)求和，即可得尾焰的吸收系數(shù)。

3 尾焰紅外輻射的計(jì)算

3.1 確定核心計(jì)算區(qū)域

根據(jù)對(duì)矩形噴管外尾焰流場(chǎng)核心區(qū)域的分析，當(dāng)尾焰的密度和壓強(qiáng)在噴口附近取得最大值后，隨著噴口徑向和軸向距離的增加，尾焰對(duì)紅外輻射吸收和發(fā)射作用的影響越小。可以用指定的CO2百分比含量ηCO2作為設(shè)定尾焰輻射計(jì)算區(qū)域的閾值，并對(duì)區(qū)域進(jìn)行簡(jiǎn)單的填補(bǔ)式修正，得到長(zhǎng)方體形尾焰[10]。具體的方法為:在X 軸方向上找到尾焰流場(chǎng)的CO2百分比含量大于等于ηCO2時(shí)的最大坐標(biāo)Xmax，同理在Y 軸方向上和Z 軸方向上找到最大坐標(biāo)Ymax和Zmax.把Zmax，2Ymax，2Xmax分別作為輻射計(jì)算區(qū)域的長(zhǎng)寬高，從而得到一個(gè)長(zhǎng)方體形的核心輻射計(jì)算區(qū)域，如圖4 所示。對(duì)計(jì)算區(qū)域進(jìn)行網(wǎng)格化，NX=8，NY=12，NZ=20，如圖5 所示，則每個(gè)網(wǎng)格(微控制體)的編碼為(nX，nY，nZ).

圖4 尾焰核心區(qū)域示意圖Fig.4 Sketch of plume core area

圖5 計(jì)算區(qū)域網(wǎng)格化示意圖Fig.5 Gridding of computational area

3.2 光譜亮度計(jì)算

FVM 的基本思想是保證微控制體在每個(gè)立體角內(nèi)的輻射能量守恒[11]。為此，需要對(duì)計(jì)算區(qū)域和4π 空間分別進(jìn)行空間離散和角度離散。前文已經(jīng)對(duì)尾焰輻射區(qū)域進(jìn)行了空間離散，得到的單個(gè)微控制體P，其體積為VP.

對(duì)4π 空間進(jìn)行角度離散，離散為互不重疊的立體角Ωmn，如圖6 所示。Ωmn中心對(duì)應(yīng)的天頂角和方位角分別為θm和φn.

對(duì)(5)式在控制體P 和立體角Ωmn積分，得到[15]

圖6 立體角離散示意圖Fig.6 Discrete solid angle

式中:re為輻射亮度沿Ωmn中心方向的單位矢量，在直角坐標(biāo)系中的計(jì)算表達(dá)式為

經(jīng)過(guò)分析與計(jì)算，尾焰計(jì)算區(qū)域中的每個(gè)控制體都會(huì)得到類似于(8)式的方程，從而形成以控制體中心節(jié)點(diǎn)沿離散立體角中心方向上的光譜輻射亮度為未知數(shù)的方程組，聯(lián)合迭代求解該方程組，便可求得任一控制體在任何離散立體角中心方向上的光譜輻射亮度。

3.3 紅外輻射強(qiáng)度計(jì)算

在得到尾焰某一微控制體的方向光譜輻射亮度Lλ(θm，φn)后，還必須利用(10)式和(11)式計(jì)算尾焰的某一微控制體外表面在某一方向的光譜輻射強(qiáng)度ΔIλ(θm，φn)和波段輻射強(qiáng)度ΔIλ1－λ2(θm，φn).通過(guò)對(duì)尾焰所有控制體外表面同一方向的輻射強(qiáng)度求和，如(12)式，進(jìn)一步得到尾焰在該方向的波段總輻射強(qiáng)度Iλ1－λ2(θm，φn).

式中:φ 為探測(cè)方向與微控制體外表面法線的夾角;ΔA 為微控制體外表面的面積。

式中:λ1－λ2為某一波段;A 為尾焰核心計(jì)算區(qū)域所有外表面。

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

4.1 紅外輻射光譜特性

隨機(jī)取核心計(jì)算區(qū)域表面的一編碼為(8，5，12)的微控制體，其外側(cè)表面中心沿θ=5π/28，φ=π/24 rad 和θ=13π/28 rad，φ =11π/24 rad 方向上的光譜輻射亮度曲線如圖7 所示。在圖7 的兩條光譜輻射亮度曲線上，都存在幾處較為明顯的輻射峰，輻射峰值波長(zhǎng)分別在2.7、4.3、5.97 μm 和6.55 μm附近，CO2在2.7 μm、4.3μm 處附近產(chǎn)生了2 個(gè)強(qiáng)吸收帶，H2O在2.7、5.97 μm 和6.55 μm 處附近產(chǎn)生了3 個(gè)強(qiáng)吸收帶，其中，左起第一處輻射峰2.7 μm的出現(xiàn)是H2O 與CO2共同作用的結(jié)果。在圖7(a)中，輻射亮度的最大值為6 ×10－2W/(cm2·μm·sr);而在圖7(b)中，輻射亮度的最大值為9 ×10－2W/(cm2·μm·sr)，后者約是前者1.5 倍，產(chǎn)生這種現(xiàn)象的原因是:隨著輻射方向的變化，輻射傳輸經(jīng)過(guò)尾焰區(qū)域的溫度和光學(xué)程長(zhǎng)是變化的，導(dǎo)致光譜輻射亮度發(fā)生變化。

圖7 控制體(8，5，12)在不同方向的光譜輻射亮度Fig.7 Spectral radiant luminance of control volume(8，5，12)in different directions

4.2 紅外輻射強(qiáng)度分布

對(duì)整個(gè)核心計(jì)算區(qū)域在寬邊對(duì)稱面XOZ 和窄邊對(duì)稱面YOZ 面內(nèi)的探測(cè)范圍均為－13π/28 ～13π/28 rad、探測(cè)間隔為π/14 rad 的各個(gè)方向利用(12)式，計(jì)算得到尾焰在3 ～5 μm 波段紅外輻射強(qiáng)度分布，如圖8 所示。圖8(a)和圖8(b)為尾焰分別在XOZ 與YOZ 平面內(nèi)的紅外輻射強(qiáng)度分布。在XOZ 平面內(nèi)，尾焰在θ =π/28 rad 方向的紅外輻射強(qiáng)度IXOZ僅為267 W/sr，隨著θ 的增大，尾焰紅外輻射強(qiáng)度迅速增大，當(dāng)在θ =13π/28 rad 方向時(shí)，紅外輻射強(qiáng)度IXOZ增大到1 965 W/sr，紅外輻射強(qiáng)度增加了7 倍;在YOZ 平面內(nèi)，尾焰在θ=π/28 rad 方向的紅外輻射強(qiáng)度IYOZ僅為142 W/sr，隨著θ 的增大，尾焰紅外輻射強(qiáng)度迅速增大，當(dāng)在θ =13π/28 rad 方向時(shí)，紅外輻射強(qiáng)度IYOZ增大到69 W/sr，紅外輻射增加了約5 倍;YOZ 平面內(nèi)最大紅外輻射強(qiáng)度IYOZ，max僅是XOZ 平面內(nèi)的最大紅外輻射強(qiáng)度IXOZ，max的35%，這主要因?yàn)榫匦螄姽艿膶捀弑華R 為2∶1，導(dǎo)致尾焰為扁平狀，所以造成了寬邊對(duì)稱面內(nèi)的IXOZ，max大于窄邊對(duì)稱面內(nèi)的IYOZ，max.

圖8 尾焰3 ～5 μm 波段總輻射強(qiáng)度分布Fig.8 Total intensity distribution in 3 ～5 μm of exhaust plume

5 結(jié)論

本文在只考慮燃?xì)庵蠬2O、CO2的光譜吸收與發(fā)射影響的情況下，針對(duì)矩形噴管外尾焰流場(chǎng)進(jìn)行了數(shù)值模擬計(jì)算，并采用FVM 和氣體輻射的窄譜帶模型，計(jì)算了矩形噴管尾焰的紅外輻射的光譜特性與在3 ～5 μm 波段的總強(qiáng)度分布。最終的模擬計(jì)算結(jié)果與文獻(xiàn)[9 －10]所得實(shí)驗(yàn)數(shù)據(jù)相吻合，證實(shí)了模型建立的合理性與數(shù)值計(jì)算結(jié)果的有效性。并由結(jié)果分析，初步得出以下結(jié)論:

1)在大氣窗口3 ～5 μm 內(nèi)尾焰輻射出現(xiàn)了2.7 μm和4.3 μm 兩個(gè)輻射峰，且4.3 μm 處的光譜輻射亮度較2.7 μm 處的要大，因此一般采用中心工作波長(zhǎng)為4.3 μm 的紅外探測(cè)器探測(cè)飛機(jī)的尾焰輻射。

2)矩形噴管的尾焰為扁平狀，且寬邊對(duì)稱面內(nèi)的紅外輻射強(qiáng)度要大于窄邊對(duì)稱面內(nèi)的紅外輻射強(qiáng)度，因此紅外探測(cè)器在尾焰寬邊對(duì)稱面內(nèi)將更容易探測(cè)到飛機(jī)的中波輻射。

3)本文的計(jì)算方法與結(jié)果可為實(shí)現(xiàn)飛機(jī)紅外隱身和紅外探測(cè)隱身飛機(jī)提供理論依據(jù)。

在實(shí)際的飛機(jī)尾焰中，必然有固體粒子存在(如未完全燃燒的碳粒子)，也就存在固體粒子對(duì)紅外輻射的散射，但是本文在求解輻射傳遞方程時(shí)，忽略了散射項(xiàng)，這必將造成一定的計(jì)算誤差，因此，本文數(shù)值計(jì)算結(jié)果精確度有進(jìn)一步提高的空間。

References)

[1]魏鋼.F－22“猛禽”戰(zhàn)斗機(jī)[M].北京:航空工業(yè)出版社，2008:5 －8，80 －81.WEI Gang.F－22“Accipiter”fighter plane[M].Beijing:Aviation Industry Press，2008:5 －8，80 －81.(in Chinese)

[2]Varney G E.Infrared signature measurement techniques and simulation methods for aircraft survivability，AIAA－1979－1186[R].Las Vegas:AIAA，SAE and ASME Joint Propulsion Conference，1979.

[3]Jim C.F/A－22 IR signature flight test model validation[J].Aircraft Survivablity，2003，29(4):9212 －9216.

[4]董士奎，于建國(guó)，李東輝，等.貼體坐標(biāo)系下離散坐標(biāo)法計(jì)算尾噴焰輻射特[J].上海理工大學(xué)學(xué)報(bào)，2003，25(2):159－162.DONG Shi－kui，YU Jian－guo，LI Dong－h(huán)ui，et al.Numerical modeling of infrared radiation properties of exhaust plume by the discrete ordinates method in body－fitted coordinates[J].Journal of University of Shanghai for Science and Technology，2003，25(2):159 －162.(in Chinese)

[5]樊士偉，張小英，朱定強(qiáng)，等.用FVM 法計(jì)算固體火箭羽流的紅外特性[J].宇航學(xué)報(bào)，2005，26(6):793 －797.FAN Shi－wei，ZHANG Xiao－ying，ZHU Ding－qiang，et al.Calculation of the infrared characteristics of the solid rocket plume with FVM method[J].Journal of Astronautics，2005，26(6):793 －797.(in Chinese)

[6]帥永，董士奎，劉林華.高溫含粒子自由流紅外輻射特性的反向蒙特卡羅法模擬[J].紅外與毫米波學(xué)報(bào)，2005，24(2):100－104.SHUAI Yong，DONG Shi－kui，LIU Lin－h(huán)ua.Simulation of infrared radiation characteristics of high temperature free－stream flow including particles by using backward Monte－Carlo method[J].Journal Infrared Millimeter and Waves，2005，24(2):100 －104.(in Chinese)

[7]阮立明，齊宏，王圣剛，等.導(dǎo)彈尾噴焰目標(biāo)紅外特性的數(shù)值仿真[J].紅外與激光工程，2008，37(6):959 －962.RUAN Li－ming，QI Hong，WANG Sheng－gang，et al.Numerical simulation of the infrared characteristic of missile exhaust plume[J].Infrared and Laser Engineering，2008，37(6):959 －962.(in Chinese)

[8]李進(jìn)良，李承曦.精通Fluent 6.3 流場(chǎng)分析[M].北京:化學(xué)工業(yè)出版社，2009:55 －58.LI Jin－liang，LI Cheng－xi.Mastery of flow field analysis in Fluent 6.3[M].Beijing:Chemical Industry Press，2009:55 －58.(in Chinese)

[9]范仁鈺，張靖周，單勇.遮擋板結(jié)構(gòu)對(duì)二元噴管紅外輻射特性的影響[J].航空動(dòng)力學(xué)報(bào)，2011，26(2):333 －348.FAN Ren－yu，ZHANG Jing－zhou，SHAN Yong.Effects of sheltering baffles on the infrared radiation characteristics of two－dimensional nozzles[J].Journal of Aerospace Power，2011，26(2):333－348.(in Chinese)

[10]侯曉春，季鶴鳴，劉慶國(guó)，等.高性能航空燃?xì)廨啓C(jī)燃燒技術(shù)[M].北京:國(guó)防工業(yè)出版社，2002:375 －378.HOU Xiao－chun，JI He－ming，LIU Qing－guo，et al.Combustion technology for high performance aviation gas turbine[M].Beijing:National Defense Industry Press，2002:375 － 378.(in Chinese)

[11]Siegel R，Howell J R.Thermal radiation heat transfer[M].Washington DC:Hemisphere and McGraw－Hill，1981.

[12]王偉臣，李世鵬，張嶠，等.工作條件對(duì)固體發(fā)動(dòng)機(jī)羽流溫度場(chǎng)的影響[J].兵工學(xué)報(bào)，2011，32(12):1493 －1497.WANG Wei－chen，LI Shi－peng，ZHANG Qiao，et al.Influence of operating conditions on temperature distributions of exhaust plume of solid rocket motor[J].Acta Armamentarii，2011，32(12):1493 －1497.(in Chinese)

[13]Ludwing C B，Malkus W，Reardon J E，et al.Hand－book of infrared radiation from combustion[R].US:NASA－SP－3080，1973:222 －230.

[14]Liu F，Gülder ? L，Smallwood G J.Three－dimensional non－grey gas radiative transfer analyses using the statistical narrow－band model[J].International Journal of Heat Mass Transfer，1998，37(9):759 －768.

[15]談和平，夏新林，劉林華，等.紅外輻射特性與傳輸?shù)臄?shù)值計(jì)算[M].哈爾濱:哈爾濱工業(yè)大學(xué)出版社，2006:194 －198.TAN He－ping，XIA Xin－lin，LIU Lin－h(huán)ua，et al.Numerical calculation of transmission and characteristic of infrared radiation[M].Harbin:Harbin Institute of Technology Press，2006:194－198.(in Chinese)