鄭鴻儒，馬 巖 ，范林東，任 翔

（1.北京跟蹤與通信技術(shù)研究所, 北京 100094；2.長光衛(wèi)星技術(shù)有限公司, 長春 130000；3.北京應(yīng)用物理與計算數(shù)學(xué)研究所, 北京 100094）

1 引 言

隨著空間技術(shù)的不斷發(fā)展和進步，空間資源對國防的戰(zhàn)略意義日益凸顯，如何在未來空天一體化戰(zhàn)爭中占領(lǐng)制高點、掌握主動權(quán)已成為各國空間技術(shù)發(fā)展的重點。天基紅外預(yù)警衛(wèi)星通過在衛(wèi)星上搭載紅外探測設(shè)備，可監(jiān)測空間目標(biāo)散發(fā)的紅外輻射，進行早期的預(yù)警及跟蹤，為實施導(dǎo)彈攔截提供可靠的彈道信息和時間差[1-2]。

高超聲速飛行器在大氣層內(nèi)飛行時，往往需要動力系統(tǒng)以調(diào)整姿態(tài)或者速度保持。此時高超聲速來流與動力系統(tǒng)的噴流相互作用，產(chǎn)生復(fù)雜的干擾流場。從目標(biāo)探測的角度來說，離軌發(fā)動機工作時產(chǎn)生的復(fù)雜流場會對紅外探測提供額外的特征信號，因此，開展紅外特性研究具有重要意義。

由于涉及領(lǐng)域的特殊性，在軌試驗研究較少，試驗研究主要以地面模型試驗為主，且成本較高。如余曉婭等人[3]使用光譜儀對半圓柱模型以再入速度8 km/s的高超聲速流動試驗進行了測量，得到了繞流激波的三維效應(yīng)，說明了激波層內(nèi)輻射存在動態(tài)非平衡效應(yīng)。在數(shù)值計算方面，楊霄、牛青林等人[4-5]采用求解Navier-Stokes（NS）方程方法獲得了類HTV-2高超聲速滑翔飛行器的流場參數(shù)，采用視在光線法（LOS）計算得到其紅外輻射特性，并分析了天基和地基平臺的可觀測性。江濤等人[6]針對印度烈火-Ⅱ?qū)?，采用窄帶輻射模型開展了助推段和再入段的輻射特性計算分析，考慮了CO2、H2O和CO組分的紅外輻射。王少平等人[7]對助推滑翔高超聲速導(dǎo)彈CAV-L的紅外輻射強度進行了仿真，并通過計算背景輻射強度給出了其被SBIRS預(yù)警衛(wèi)星探測到的探測距離。在采用反向蒙特卡洛法（BMC）求解輻射傳輸方程的應(yīng)用方面，包醒東等人[8]采用直接模擬蒙特卡洛方法（DSMC）計算了某0.9 N模型發(fā)動機羽流場，并采用反向蒙特卡洛方法對輻射特性進行了分析，結(jié)果表明，由于高空羽流溫度整體較低，輻射峰值向中波移動。到目前為止，國內(nèi)對再入體流場輻射仿真分析方面已有研究或關(guān)注再入時的大氣來流，或關(guān)注發(fā)動機燃燒產(chǎn)生的尾焰紅外可探測性，對來流與噴流干擾流場的分析較少，尤其是對不同攻角、不同來流速度下的干擾流場紅外輻射特性研究鮮有報道。

本文首先使用數(shù)值求解NS方程獲得某型高超聲速飛行器離軌發(fā)動機噴流與不同攻角、不同速度條件下來流形成的干擾流場，獲得流場溫度、壓力及組分分布，然后采用基于三維非結(jié)構(gòu)化網(wǎng)格開發(fā)的低軌衛(wèi)星紅外探測可見性分析軟件，建立高超聲速來流與噴流相互干擾流場的紅外輻射特性數(shù)值模擬模型，采用反向蒙特卡洛方法計算流場對給定軌道高度的紅外輻射特性，并結(jié)合大氣衰減效應(yīng)計算分析了不同位置處低軌衛(wèi)星紅外探測的可見性。

2 紅外輻射計算方法

2.1 物理模型

本文使用的物理模型及發(fā)動機參數(shù)與文獻[9]一致。飛行器本體為鈍錐形，尾部安裝有10 000 N的MMH/N2O4組合發(fā)動機，混合比為1.65±0.05，燃燒室壓強為4 MPa，燃燒室溫度為3 160 K。依據(jù)常見的飛行彈道狀態(tài)，流場數(shù)值計算中分別考慮飛行體在高度為94 km的高空無來流和以攻角為90°來流、135°來流飛行時的來流與發(fā)動機噴流的相互干擾作用，其中，來流大氣速度除考慮了標(biāo)準(zhǔn)飛行速度7 000 m/s外，還考慮了5 000 m/s和6 000 m/s的不同來流速度以供對比，大氣組分為79%的N2和21%的O2，溫度和壓強設(shè)置參考1976年美國標(biāo)準(zhǔn)大氣。發(fā)動機的燃燒產(chǎn)物及其比例根據(jù)對噴管的熱力學(xué)計算得到，流場計算時考慮燃氣的主要成分為CO2、CO、H2O、O2、N2、H2，而來流條件通過設(shè)置為遠場邊界條件實現(xiàn)。在計算紅外特性時主要考慮CO2和H2O兩種成分。

2.2 數(shù)值求解NS方程

發(fā)動機噴流與來流干擾流場的數(shù)值仿真通過數(shù)值求解Navier-Stokes方程方法實現(xiàn)，控制方程如公式(1)所示：

其中，U為守恒量，F(xiàn)i,inv為無粘通量、Fi,vis為粘性通量，W為源項。

采用Fluent軟件求解基于密度的二階CFD穩(wěn)態(tài)求解器，耦合求解連續(xù)性方程、組分方程、動量方程、能量方程。湍流模型采用sst k-ω兩方程湍流模型。

2.3 氣體輻射參數(shù)計算

本文主要考慮的輻射組分CO2和H2O的氣體輻射參數(shù)基于HITRAN 2008[10]及HITEMP 2010[11]數(shù)據(jù)庫，采用Voigt線型函數(shù)的逐線積分法（LBL）進行求解。對于混合氣體的光譜吸收系數(shù)通過各組分吸收系數(shù)加權(quán)求和[12]，如公式(2)所示：

其中C是組分的摩爾分?jǐn)?shù)，κ是標(biāo)準(zhǔn)溫度下組分的光譜吸收系數(shù)（cm-1），T是組分的實際溫度（K），p是組分的實際壓強（Pa）。

圖1（彩圖見期刊電子版）中給出了本文計算的CO2和H2O在1 000 K溫度下 2.7 μm譜帶的吸收系數(shù)與文獻[13]和文獻[14]中給出的結(jié)果對比。其中藍色三角曲線為文獻參考數(shù)值，紅色方形曲線為本文計算數(shù)值?？梢钥闯鑫章视嬎憬Y(jié)果一致性較好，僅在局部峰值存在一定差異，這與選用的數(shù)據(jù)庫參數(shù)更新等因素有關(guān)，偏差在可接受范圍，驗證了氣體輻射參數(shù)計算的正確性。

圖1 吸收系數(shù)對比Fig.1 Comparison of absorption coefficients

2.4 輻射計算方法

本文求解輻射傳輸方程采用的是反向蒙特卡洛法（BMCM）。反向蒙特卡洛法是蒙特卡洛法的逆過程，基于輻射傳遞的互易性原理，從探測點開始逆向追蹤光線，在傳輸過程中判斷光線是否被吸收，直到光線被吸收或者超出計算邊界，然后將光線的吸收點作為歸宿點，計算歸宿點對探測點的輻射貢獻量。由于不需要對發(fā)射點進行空間離散，BMCM在處理各向異性輻射、復(fù)雜形狀、小立體角問題時非常有效，提高了計算效率，保證了計算精度，適用于本文中假設(shè)的低軌遙感衛(wèi)星紅外探測過程。

光線從觀測點發(fā)射，在給定的視場角內(nèi)隨機分布，統(tǒng)計所有光線的熱流密度為：

其中I為方向輻射強度，F(xiàn)OV為紅外輻射源對觀測點所張的立體角。

將發(fā)射的N條光線中進入流場計算得到的Na條光線的熱流密度除以總光線數(shù)得到輻射強度Icos?在視場內(nèi)的平均值為：

最終到達觀測點的輻射熱流密度為：

由于壁面輻射遠小于流場輻射，因此，本文僅考慮流場輻射的計算和可觀測性分析。

2.5 計算域及觀測視角

本文中使用的坐標(biāo)系定義發(fā)動機出口平面中心為原點，出口方向為+X方向，豎直向上方向為+Y方向，+Z遵循右手定則，來流攻角定義為與X軸的夾角。計算域為以50 m半徑在出口平面上形成的半球空間。觀測視角采用方位角 θ和俯仰角 φ兩個角度進行定義，如圖2所示，其中方位角在Y0Z平面上以X軸為旋轉(zhuǎn)軸順時針變化，俯仰角在XOY平面上以Z軸為旋轉(zhuǎn)軸順時針變化。文中設(shè)置方位角度間隔為30°，變化范圍為-90°～+90°，俯仰角設(shè)置為0°。觀測距離參考528 km軌道衛(wèi)星對94 km飛行高度流場進行觀測，在各工況下均設(shè)置為434 km。由于再入段反導(dǎo)攔截常使用3～5 μm和8～12 μm紅外焦平面陣列進行觀測[6]。參考吉林一號光譜星紅外載荷參數(shù)，本文選擇計算波段為2.65～3.00、3.70～4.95和8.00～13.50 μm，波長間隔為0.05 μm。

圖2 計算域及觀測視角示意圖Fig.2 Schematic diagram of computational domain and observation angle

3 結(jié)果與討論

3.1 流場計算結(jié)果分析

根據(jù)噴管擴張段型面及發(fā)動機熱力學(xué)參數(shù)，使用 Fluent 計算了高度為94 km時，無風(fēng)條件以及攻角分別為90°、135°來流條件下的流場分布[9]，來流速度設(shè)置為7 000 m/s。為了進一步分析同一攻角下不同來流速度對尾焰紅外輻射特性的影響，對135°攻角狀態(tài)下，5 000 m/s和6 000 m/s來流情況進行了仿真分析，并與7 000 m/s工況進行了對比。圖3～圖6（彩圖見期刊電子版）分別給出了幾種工況下的溫度、H2O和CO2質(zhì)量分布云圖。由圖3可以看出，94 km高度無風(fēng)條件下，發(fā)動機尾焰流場呈軸對稱結(jié)構(gòu)，充分膨脹并向后方流動，由于發(fā)動機燃氣的內(nèi)能大部分轉(zhuǎn)化為動能，此時噴流核心區(qū)域的靜溫會降到100 K以下。另一方面，在稀薄大氣約束下，發(fā)動機燃氣與大氣形成拋物線型的高溫剪切層，其溫度峰值為1 300 K左右。在攻角分別為90°和135°的來流條件下，迎風(fēng)側(cè)的高溫剪切層轉(zhuǎn)變?yōu)楣渭げǎ瑴囟燃眲∩仙?，峰值? 000 K以上，并且噴流核心區(qū)域被壓縮在小范圍空間內(nèi)。且從圖中可以看出，來流攻角越大，干擾流場影響范圍越大。

圖3 不同攻角下尾焰流場溫度云圖Fig.3 Temperature fraction of exhaust plume at different attack angles

圖4 不同來流速度下尾焰流場溫度云圖Fig.4 Temperature nephogram of exhaust plume at different inlet velocities

圖5 不同攻角下尾焰流場CO2質(zhì)量分?jǐn)?shù)云圖Fig.5 CO2 mass fraction nephogram of exhaust plume at different attack angles

圖6 不同攻角下尾焰流場H2O質(zhì)量分?jǐn)?shù)云圖Fig.6 H2O mass fraction nephogram of exhaust plume at different attack angles

圖4給出了不同風(fēng)速條件下的流場溫度分布，可以看出，在135°攻角下，隨著來流速度的逐漸升高，噴流核心區(qū)域被壓縮在更小的區(qū)域內(nèi)，且噴流與來流間的弓形激波的溫度不斷升高，最高溫度分別為7 940、10 570、13 460 K。

由圖5～圖6可以看出，在噴流核心區(qū)的H2O和CO2濃度與發(fā)動機出口處的燃氣濃度保證一致，而在高溫剪切層內(nèi)，由于物質(zhì)摻混和擴散作用，氣體濃度存在梯度。相較于無風(fēng)工況，有來流條件下噴流核心區(qū)域的減小造成H2O和CO2的高濃度區(qū)域減小，不過在弓形激波和高溫剪切層位置處存在濃度變化。135°攻角下，5 000 m/s和6 000 m/s來流的H2O和CO2組分分布情況與7 000 m/s來流分布相近，此處不再贅述。

3.2 流場輻射特性分析

圖7～圖9（彩圖見期刊電子版）分別給出了2.65～3 μm波段、3.7～4.95 μm波段和8～13.5 μm波段在無來流、攻角90°、攻角135°流場條件下的輻射分布對比?？梢钥闯?，各波段輻射熱流密度受觀測角度的影響很小，而不同工況下影響較大，且3.70～4.95 μm波段輻射強于其他波段。2.65～3 μm波段包括了H2O在2.66 μm和CO2在2.7 μm附近的輻射，3.7～4.95 μm波段包括了CO2在4.3 μm左右的輻射，8～13.5 μm波段包括了H2O在12 μm 左右的輻射。

圖7 2.65～3.00 μm波段輻射熱流密度對比Fig.7 Comparison of radiative heat flux in 2.65～3.00 μm wave band

圖8 3.70～4.95 μm波段輻射熱流密度對比Fig.8 Comparison of radiative heat flux in 3.70～4.95 μm wave band

圖9 8.00～13.5 μm波段輻射熱流密度對比Fig.9 Comparison of radiative heat flux in 8.00～13.5 μm wave band

在無風(fēng)條件下，大氣與飛行器相對靜止，此時產(chǎn)生的紅外特征信號強度較低，3.7～4.95 μm波段在10-9W/m2量級，2.65～3 μm和8～13.5 μm波段在10-10W/m2量級。

當(dāng)大氣來流以90°和135°攻角與尾焰流場作用后產(chǎn)生的紅外信號強度升高了2個數(shù)量級以上，且攻角越大，干擾流場影響范圍越大，紅外信號強度也隨之增大。當(dāng)存在攻角時，2.65～3 μm和3.7～4.95 μm波段紅外輻射強度在10-7W/m2量級，8～13.5 μm波段在10-8W/m2量級，可見，攻角對2.65～3 μm波段輻射強度提升效果最為明顯。

圖10給出了135°攻角下，5 000 m/s和6 000 m/s來流速度下的輻射熱流密度與7 000 m/s來流的對比。圖中w1、w2、w3分別代表2.65～3 μm、3.7～4.95 μm和8～13.5 μm波段。從圖中可以看出，5 000 m/s和6 000 m/s來流工況的輻射熱流密度在各個波段都與7 000 m/s來流工況較為相近，同時也存在速度越大、輻射熱流密度越大的趨勢。

圖10 不同來流速度下各波段輻射熱流密度對比Fig.10 Comparison of radiative heat flux at different wave bands and different inlet velocities

3.3 可探測性分析

本文設(shè)定的94 km飛行高度，由于處于稠密大氣邊緣，采用528 km軌道衛(wèi)星天頂方向觀測時（圖11中位置1），大氣衰減效應(yīng)可以忽略不計。

圖11 低軌衛(wèi)星觀測方向示意圖Fig.11 Detection diagram of a low orbit satellite

當(dāng)衛(wèi)星位于位置2時，由于穿過大氣需要考慮沿程衰減，因此采用大氣輻射傳輸軟件Modtran計算各波數(shù)下大氣透過率，疊加至輻射光譜中進行積分計算，如公式(6)所示[4]。Modtran中使用的大氣條件為1976年美國標(biāo)準(zhǔn)大氣，氣溶膠模型為海軍海洋模型。

其中，R為衛(wèi)星與飛行器之間的距離，Iλ為λ1至λ2波段內(nèi)的輻射強度，τ(λ,R)和Lpath,λ分別為大氣透過率和沿程輻射亮度，由Modtran軟件計算得到，At為流場在觀測方向上的面積。此處僅考慮大氣透過率，得到在位置1和位置2處3種來流工況下的紅外輻射熱流密度如表1所示。

根據(jù)文獻和實際在軌工作傳感器過往數(shù)據(jù)，在計算時將2.65～3 μm波段和3.7～4.95 μm波段紅外傳感器靈敏度取為10-12W/cm2[5]，8～13.5 μm波段靈敏度取為4×10-12W/cm2?？梢钥闯?，在位置1處，不同攻角、不同速度的幾種工況對于3個波段均可探測，而無風(fēng)條件下則均無法探測。在位置2處，由于探測距離明顯增大，探測路徑通過稠密大氣距離較長，衰減效應(yīng)明顯，對于無風(fēng)、不同攻角及不同速度來流條件，3個波段均不可探測。對于可探測的位置1處，可以看出2.65～3 μm波段和3.7～4.95 μm波段強度高于8～13.5 μm波段，這與已知的紅外預(yù)警普遍采用的2.7 μm和4.3 μm 敏感器件的情況相符。

表1 不同觀測位置的輻射熱流密度Tab.1 Radiative heat fluxes at different observation positions (W/cm2)

4 結(jié) 論

本文針對高超聲速飛行器離軌發(fā)動機工作時產(chǎn)生的尾焰與來流干擾形成的流場紅外輻射特性進行了仿真分析。仿真結(jié)果表明，相比無來流的情況，攻角為90°和135°時由于產(chǎn)生弓形激波，流場內(nèi)溫度急劇升高，產(chǎn)生的干擾流場紅外輻射強度顯著升高2～3個數(shù)量級以上，相同工況下長波波段較3.70～4.95 μm波段低1個數(shù)量級以上。對于同一135°攻角下的不同速度工況開展了紅外輻射特性分析。結(jié)果表明，隨著來流速度的增大，紅外輻射熱流密度也越來越高。同時，還結(jié)合大氣衰減分析了低軌衛(wèi)星在不同位置處對流場紅外信號的可觀測性。結(jié)果表明，針對2.65～3 μm波段和3.7～4.95 μm波段紅外傳感器靈敏度取為10-12W/cm2，8～13.5 μm波段靈敏度取為4×10-12W/cm2，當(dāng)衛(wèi)星對地且飛行器位于視場中心時除無風(fēng)條件不可探測，不同攻角、不同速度的來流條件下均可探測，而當(dāng)衛(wèi)星側(cè)擺通過臨邊路徑進行探測時由于穿過稠密大氣，不同來流條件不同波段下均不可探測。