鄭鴻儒,馬 巖 ,范林東,任 翔
(1.北京跟蹤與通信技術(shù)研究所, 北京 100094;2.長(zhǎng)光衛(wèi)星技術(shù)有限公司, 長(zhǎng)春 130000;3.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所, 北京 100094)
隨著空間技術(shù)的不斷發(fā)展和進(jìn)步,空間資源對(duì)國(guó)防的戰(zhàn)略意義日益凸顯,如何在未來(lái)空天一體化戰(zhàn)爭(zhēng)中占領(lǐng)制高點(diǎn)、掌握主動(dòng)權(quán)已成為各國(guó)空間技術(shù)發(fā)展的重點(diǎn)。天基紅外預(yù)警衛(wèi)星通過(guò)在衛(wèi)星上搭載紅外探測(cè)設(shè)備,可監(jiān)測(cè)空間目標(biāo)散發(fā)的紅外輻射,進(jìn)行早期的預(yù)警及跟蹤,為實(shí)施導(dǎo)彈攔截提供可靠的彈道信息和時(shí)間差[1-2]。
高超聲速飛行器在大氣層內(nèi)飛行時(shí),往往需要?jiǎng)恿ο到y(tǒng)以調(diào)整姿態(tài)或者速度保持。此時(shí)高超聲速來(lái)流與動(dòng)力系統(tǒng)的噴流相互作用,產(chǎn)生復(fù)雜的干擾流場(chǎng)。從目標(biāo)探測(cè)的角度來(lái)說(shuō),離軌發(fā)動(dòng)機(jī)工作時(shí)產(chǎn)生的復(fù)雜流場(chǎng)會(huì)對(duì)紅外探測(cè)提供額外的特征信號(hào),因此,開展紅外特性研究具有重要意義。
由于涉及領(lǐng)域的特殊性,在軌試驗(yàn)研究較少,試驗(yàn)研究主要以地面模型試驗(yàn)為主,且成本較高。如余曉婭等人[3]使用光譜儀對(duì)半圓柱模型以再入速度8 km/s的高超聲速流動(dòng)試驗(yàn)進(jìn)行了測(cè)量,得到了繞流激波的三維效應(yīng),說(shuō)明了激波層內(nèi)輻射存在動(dòng)態(tài)非平衡效應(yīng)。在數(shù)值計(jì)算方面,楊霄、牛青林等人[4-5]采用求解Navier-Stokes(NS)方程方法獲得了類HTV-2高超聲速滑翔飛行器的流場(chǎng)參數(shù),采用視在光線法(LOS)計(jì)算得到其紅外輻射特性,并分析了天基和地基平臺(tái)的可觀測(cè)性。江濤等人[6]針對(duì)印度烈火-Ⅱ?qū)棧捎谜瓗л椛淠P烷_展了助推段和再入段的輻射特性計(jì)算分析,考慮了CO2、H2O和CO組分的紅外輻射。王少平等人[7]對(duì)助推滑翔高超聲速導(dǎo)彈CAV-L的紅外輻射強(qiáng)度進(jìn)行了仿真,并通過(guò)計(jì)算背景輻射強(qiáng)度給出了其被SBIRS預(yù)警衛(wèi)星探測(cè)到的探測(cè)距離。在采用反向蒙特卡洛法(BMC)求解輻射傳輸方程的應(yīng)用方面,包醒東等人[8]采用直接模擬蒙特卡洛方法(DSMC)計(jì)算了某0.9 N模型發(fā)動(dòng)機(jī)羽流場(chǎng),并采用反向蒙特卡洛方法對(duì)輻射特性進(jìn)行了分析,結(jié)果表明,由于高空羽流溫度整體較低,輻射峰值向中波移動(dòng)。到目前為止,國(guó)內(nèi)對(duì)再入體流場(chǎng)輻射仿真分析方面已有研究或關(guān)注再入時(shí)的大氣來(lái)流,或關(guān)注發(fā)動(dòng)機(jī)燃燒產(chǎn)生的尾焰紅外可探測(cè)性,對(duì)來(lái)流與噴流干擾流場(chǎng)的分析較少,尤其是對(duì)不同攻角、不同來(lái)流速度下的干擾流場(chǎng)紅外輻射特性研究鮮有報(bào)道。
本文首先使用數(shù)值求解NS方程獲得某型高超聲速飛行器離軌發(fā)動(dòng)機(jī)噴流與不同攻角、不同速度條件下來(lái)流形成的干擾流場(chǎng),獲得流場(chǎng)溫度、壓力及組分分布,然后采用基于三維非結(jié)構(gòu)化網(wǎng)格開發(fā)的低軌衛(wèi)星紅外探測(cè)可見性分析軟件,建立高超聲速來(lái)流與噴流相互干擾流場(chǎng)的紅外輻射特性數(shù)值模擬模型,采用反向蒙特卡洛方法計(jì)算流場(chǎng)對(duì)給定軌道高度的紅外輻射特性,并結(jié)合大氣衰減效應(yīng)計(jì)算分析了不同位置處低軌衛(wèi)星紅外探測(cè)的可見性。
本文使用的物理模型及發(fā)動(dòng)機(jī)參數(shù)與文獻(xiàn)[9]一致。飛行器本體為鈍錐形,尾部安裝有10 000 N的MMH/N2O4組合發(fā)動(dòng)機(jī),混合比為1.65±0.05,燃燒室壓強(qiáng)為4 MPa,燃燒室溫度為3 160 K。依據(jù)常見的飛行彈道狀態(tài),流場(chǎng)數(shù)值計(jì)算中分別考慮飛行體在高度為94 km的高空無(wú)來(lái)流和以攻角為90°來(lái)流、135°來(lái)流飛行時(shí)的來(lái)流與發(fā)動(dòng)機(jī)噴流的相互干擾作用,其中,來(lái)流大氣速度除考慮了標(biāo)準(zhǔn)飛行速度7 000 m/s外,還考慮了5 000 m/s和6 000 m/s的不同來(lái)流速度以供對(duì)比,大氣組分為79%的N2和21%的O2,溫度和壓強(qiáng)設(shè)置參考1976年美國(guó)標(biāo)準(zhǔn)大氣。發(fā)動(dòng)機(jī)的燃燒產(chǎn)物及其比例根據(jù)對(duì)噴管的熱力學(xué)計(jì)算得到,流場(chǎng)計(jì)算時(shí)考慮燃?xì)獾闹饕煞譃镃O2、CO、H2O、O2、N2、H2,而來(lái)流條件通過(guò)設(shè)置為遠(yuǎn)場(chǎng)邊界條件實(shí)現(xiàn)。在計(jì)算紅外特性時(shí)主要考慮CO2和H2O兩種成分。
發(fā)動(dòng)機(jī)噴流與來(lái)流干擾流場(chǎng)的數(shù)值仿真通過(guò)數(shù)值求解Navier-Stokes方程方法實(shí)現(xiàn),控制方程如公式(1)所示:

其中,U為守恒量,F(xiàn)i,inv為無(wú)粘通量、Fi,vis為粘性通量,W為源項(xiàng)。
采用Fluent軟件求解基于密度的二階CFD穩(wěn)態(tài)求解器,耦合求解連續(xù)性方程、組分方程、動(dòng)量方程、能量方程。湍流模型采用sst k-ω兩方程湍流模型。
本文主要考慮的輻射組分CO2和H2O的氣體輻射參數(shù)基于HITRAN 2008[10]及HITEMP 2010[11]數(shù)據(jù)庫(kù),采用Voigt線型函數(shù)的逐線積分法(LBL)進(jìn)行求解。對(duì)于混合氣體的光譜吸收系數(shù)通過(guò)各組分吸收系數(shù)加權(quán)求和[12],如公式(2)所示:

其中C是組分的摩爾分?jǐn)?shù),κ是標(biāo)準(zhǔn)溫度下組分的光譜吸收系數(shù)(cm-1),T是組分的實(shí)際溫度(K),p是組分的實(shí)際壓強(qiáng)(Pa)。
圖1(彩圖見期刊電子版)中給出了本文計(jì)算的CO2和H2O在1 000 K溫度下 2.7 μm譜帶的吸收系數(shù)與文獻(xiàn)[13]和文獻(xiàn)[14]中給出的結(jié)果對(duì)比。其中藍(lán)色三角曲線為文獻(xiàn)參考數(shù)值,紅色方形曲線為本文計(jì)算數(shù)值??梢钥闯鑫章视?jì)算結(jié)果一致性較好,僅在局部峰值存在一定差異,這與選用的數(shù)據(jù)庫(kù)參數(shù)更新等因素有關(guān),偏差在可接受范圍,驗(yàn)證了氣體輻射參數(shù)計(jì)算的正確性。

圖1 吸收系數(shù)對(duì)比Fig.1 Comparison of absorption coefficients
本文求解輻射傳輸方程采用的是反向蒙特卡洛法(BMCM)。反向蒙特卡洛法是蒙特卡洛法的逆過(guò)程,基于輻射傳遞的互易性原理,從探測(cè)點(diǎn)開始逆向追蹤光線,在傳輸過(guò)程中判斷光線是否被吸收,直到光線被吸收或者超出計(jì)算邊界,然后將光線的吸收點(diǎn)作為歸宿點(diǎn),計(jì)算歸宿點(diǎn)對(duì)探測(cè)點(diǎn)的輻射貢獻(xiàn)量。由于不需要對(duì)發(fā)射點(diǎn)進(jìn)行空間離散,BMCM在處理各向異性輻射、復(fù)雜形狀、小立體角問(wèn)題時(shí)非常有效,提高了計(jì)算效率,保證了計(jì)算精度,適用于本文中假設(shè)的低軌遙感衛(wèi)星紅外探測(cè)過(guò)程。
光線從觀測(cè)點(diǎn)發(fā)射,在給定的視場(chǎng)角內(nèi)隨機(jī)分布,統(tǒng)計(jì)所有光線的熱流密度為:

其中I為方向輻射強(qiáng)度,F(xiàn)OV為紅外輻射源對(duì)觀測(cè)點(diǎn)所張的立體角。
將發(fā)射的N條光線中進(jìn)入流場(chǎng)計(jì)算得到的Na條光線的熱流密度除以總光線數(shù)得到輻射強(qiáng)度Icos?在視場(chǎng)內(nèi)的平均值為:

最終到達(dá)觀測(cè)點(diǎn)的輻射熱流密度為:

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

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

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

圖4 不同來(lái)流速度下尾焰流場(chǎng)溫度云圖Fig.4 Temperature nephogram of exhaust plume at different inlet velocities

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

圖6 不同攻角下尾焰流場(chǎng)H2O質(zhì)量分?jǐn)?shù)云圖Fig.6 H2O mass fraction nephogram of exhaust plume at different attack angles
圖4給出了不同風(fēng)速條件下的流場(chǎng)溫度分布,可以看出,在135°攻角下,隨著來(lái)流速度的逐漸升高,噴流核心區(qū)域被壓縮在更小的區(qū)域內(nèi),且噴流與來(lái)流間的弓形激波的溫度不斷升高,最高溫度分別為7 940、10 570、13 460 K。
由圖5~圖6可以看出,在噴流核心區(qū)的H2O和CO2濃度與發(fā)動(dòng)機(jī)出口處的燃?xì)鉂舛缺WC一致,而在高溫剪切層內(nèi),由于物質(zhì)摻混和擴(kuò)散作用,氣體濃度存在梯度。相較于無(wú)風(fēng)工況,有來(lái)流條件下噴流核心區(qū)域的減小造成H2O和CO2的高濃度區(qū)域減小,不過(guò)在弓形激波和高溫剪切層位置處存在濃度變化。135°攻角下,5 000 m/s和6 000 m/s來(lái)流的H2O和CO2組分分布情況與7 000 m/s來(lái)流分布相近,此處不再贅述。
圖7~圖9(彩圖見期刊電子版)分別給出了2.65~3 μm波段、3.7~4.95 μm波段和8~13.5 μm波段在無(wú)來(lái)流、攻角90°、攻角135°流場(chǎng)條件下的輻射分布對(duì)比。可以看出,各波段輻射熱流密度受觀測(cè)角度的影響很小,而不同工況下影響較大,且3.70~4.95 μm波段輻射強(qiáng)于其他波段。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波段輻射熱流密度對(duì)比Fig.7 Comparison of radiative heat flux in 2.65~3.00 μm wave band

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

圖9 8.00~13.5 μm波段輻射熱流密度對(duì)比Fig.9 Comparison of radiative heat flux in 8.00~13.5 μm wave band
在無(wú)風(fēng)條件下,大氣與飛行器相對(duì)靜止,此時(shí)產(chǎn)生的紅外特征信號(hào)強(qiáng)度較低,3.7~4.95 μm波段在10-9W/m2量級(jí),2.65~3 μm和8~13.5 μm波段在10-10W/m2量級(jí)。
當(dāng)大氣來(lái)流以90°和135°攻角與尾焰流場(chǎng)作用后產(chǎn)生的紅外信號(hào)強(qiáng)度升高了2個(gè)數(shù)量級(jí)以上,且攻角越大,干擾流場(chǎng)影響范圍越大,紅外信號(hào)強(qiáng)度也隨之增大。當(dāng)存在攻角時(shí),2.65~3 μm和3.7~4.95 μm波段紅外輻射強(qiáng)度在10-7W/m2量級(jí),8~13.5 μm波段在10-8W/m2量級(jí),可見,攻角對(duì)2.65~3 μm波段輻射強(qiáng)度提升效果最為明顯。
圖10給出了135°攻角下,5 000 m/s和6 000 m/s來(lái)流速度下的輻射熱流密度與7 000 m/s來(lái)流的對(duì)比。圖中w1、w2、w3分別代表2.65~3 μm、3.7~4.95 μm和8~13.5 μm波段。從圖中可以看出,5 000 m/s和6 000 m/s來(lái)流工況的輻射熱流密度在各個(gè)波段都與7 000 m/s來(lái)流工況較為相近,同時(shí)也存在速度越大、輻射熱流密度越大的趨勢(shì)。

圖10 不同來(lái)流速度下各波段輻射熱流密度對(duì)比Fig.10 Comparison of radiative heat flux at different wave bands and different inlet velocities
本文設(shè)定的94 km飛行高度,由于處于稠密大氣邊緣,采用528 km軌道衛(wèi)星天頂方向觀測(cè)時(shí)(圖11中位置1),大氣衰減效應(yīng)可以忽略不計(jì)。

圖11 低軌衛(wèi)星觀測(cè)方向示意圖Fig.11 Detection diagram of a low orbit satellite
當(dāng)衛(wèi)星位于位置2時(shí),由于穿過(guò)大氣需要考慮沿程衰減,因此采用大氣輻射傳輸軟件Modtran計(jì)算各波數(shù)下大氣透過(guò)率,疊加至輻射光譜中進(jìn)行積分計(jì)算,如公式(6)所示[4]。Modtran中使用的大氣條件為1976年美國(guó)標(biāo)準(zhǔn)大氣,氣溶膠模型為海軍海洋模型。

其中,R為衛(wèi)星與飛行器之間的距離,Iλ為λ1至λ2波段內(nèi)的輻射強(qiáng)度,τ(λ,R)和Lpath,λ分別為大氣透過(guò)率和沿程輻射亮度,由Modtran軟件計(jì)算得到,At為流場(chǎng)在觀測(cè)方向上的面積。此處僅考慮大氣透過(guò)率,得到在位置1和位置2處3種來(lái)流工況下的紅外輻射熱流密度如表1所示。
根據(jù)文獻(xiàn)和實(shí)際在軌工作傳感器過(guò)往數(shù)據(jù),在計(jì)算時(shí)將2.65~3 μm波段和3.7~4.95 μm波段紅外傳感器靈敏度取為10-12W/cm2[5],8~13.5 μm波段靈敏度取為4×10-12W/cm2??梢钥闯?,在位置1處,不同攻角、不同速度的幾種工況對(duì)于3個(gè)波段均可探測(cè),而無(wú)風(fēng)條件下則均無(wú)法探測(cè)。在位置2處,由于探測(cè)距離明顯增大,探測(cè)路徑通過(guò)稠密大氣距離較長(zhǎng),衰減效應(yīng)明顯,對(duì)于無(wú)風(fēng)、不同攻角及不同速度來(lái)流條件,3個(gè)波段均不可探測(cè)。對(duì)于可探測(cè)的位置1處,可以看出2.65~3 μm波段和3.7~4.95 μm波段強(qiáng)度高于8~13.5 μm波段,這與已知的紅外預(yù)警普遍采用的2.7 μm和4.3 μm 敏感器件的情況相符。

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