基于有限體積法的近空間高超聲速飛行器輻射特性仿真研究

1熱輻射源紅外特性的研究近距離飛機(jī)是一種新型的戰(zhàn)略遏制和戰(zhàn)術(shù)應(yīng)用武器平臺(tái)。這是世界上所有主要軍事能力的集合和研究的熱點(diǎn)。隨著紅外探測(cè)、跟蹤和制導(dǎo)技術(shù)的發(fā)展,空中軍事目標(biāo)面臨日益嚴(yán)峻的紅外威脅,具備紅外隱身能力已成為現(xiàn)代先進(jìn)飛行器提高自身生存和突防能力的必然要求。開展近空間高超聲速飛行器紅外輻射特性研究對(duì)于其紅外隱身設(shè)計(jì)具有十分重要的意義。獲取目標(biāo)的紅外特性數(shù)據(jù)可通過實(shí)驗(yàn)測(cè)量和數(shù)值仿真來實(shí)現(xiàn)。由于實(shí)驗(yàn)測(cè)量耗資巨大,實(shí)驗(yàn)結(jié)果容易受諸多環(huán)境因素的影響,建立紅外特性理論模型,通過數(shù)值仿真來獲取目標(biāo)的紅外輻射特征被廣泛采用。國(guó)外這方面的研究開展得較早,并已開發(fā)出相關(guān)計(jì)算軟件,如北約部分國(guó)家聯(lián)合開發(fā)的空中目標(biāo)的紅外輻射模型NIRATAMNATO;近年來,國(guó)內(nèi)開始重視飛行器隱身及相關(guān)方面的研究工作,羅明東等研究了非加力渦輪發(fā)動(dòng)機(jī)排氣系統(tǒng)紅外輻射強(qiáng)度的數(shù)值計(jì)算方法;張小英等研究了軸對(duì)稱矢量噴管在光譜2100～5260cm-1范圍內(nèi)紅外特性的數(shù)值計(jì)算方法;郝金波等在給定濃度場(chǎng)、溫度場(chǎng)、粒子和氣體組分輻射物性參數(shù)的條件下,將貼體坐標(biāo)系下有限體積法用于主動(dòng)段尾噴焰紅外光譜特性計(jì)算。但研究重點(diǎn)主要集中在各類常規(guī)飛行器的發(fā)動(dòng)機(jī)熱部件及尾噴焰。對(duì)于近空間高超聲速飛行器來說,氣動(dòng)加熱引起的熱輻射是其主要的紅外輻射源。本文從紅外隱身設(shè)計(jì)的角度,針對(duì)飛行馬赫數(shù)Ma=4～7、飛行高度H=15km～30km、具有長(zhǎng)時(shí)間熱載荷和非燒蝕構(gòu)型等特點(diǎn)的近空間高超聲速飛行器氣動(dòng)熱紅外特性進(jìn)行了理論建模和數(shù)值仿真研究。希望以此為基礎(chǔ)能夠進(jìn)一步建立一套完整的數(shù)值預(yù)測(cè)體系,可以為近空間高超聲速飛行器紅外隱身設(shè)計(jì)提供參考依據(jù)。2輻射場(chǎng)-流場(chǎng)當(dāng)飛行器在近空間大氣層內(nèi)以高超聲速飛行時(shí),強(qiáng)烈的氣動(dòng)熱會(huì)引起周圍介質(zhì)(空氣)各點(diǎn)的溫度、壓力等熱力學(xué)狀態(tài)在空間上的劇烈變化,從而使介質(zhì)表現(xiàn)為強(qiáng)烈的參與性和不均勻性。但由于介質(zhì)熱力學(xué)狀態(tài)只在飛行器近域空間變化比較劇烈,在無窮遠(yuǎn)處則趨向于外界環(huán)境。為此,在遠(yuǎn)場(chǎng)處引入一“假想面”將飛行器周圍空間劃分為兩部分,如圖1所示?！凹傧朊妗蓖鈬橘|(zhì)認(rèn)為是常物性的或透明的,只要知道“假想界”處輻射參量,就可得到任意遠(yuǎn)處的輻射參量值,這樣問題便轉(zhuǎn)化為求解封閉系統(tǒng)內(nèi)參與性介質(zhì)的輻射換熱問題。通常,輻射場(chǎng)和流場(chǎng)是耦合在一起的,但這里所研究的對(duì)象,速度不大于3km/s,流場(chǎng)最高溫度一般不超過3000K。根據(jù)目前文獻(xiàn),在此情況下可以不考慮它們之間的干擾效應(yīng),認(rèn)為輻射場(chǎng)與流場(chǎng)是可以解耦的。采用解耦計(jì)算的具體步驟為:首先基于氣體動(dòng)力學(xué)方程計(jì)算流場(chǎng)參量;然后以已知流場(chǎng)參量為基礎(chǔ),計(jì)算介質(zhì)輻射物性;最后通過求解介質(zhì)內(nèi)輻射傳輸方程獲得相應(yīng)輻射參量。3流場(chǎng)數(shù)值法3.1navier-st對(duì)流通量方程的數(shù)學(xué)模型假設(shè)空氣是由多種組分氣體均勻混合而成的完全氣體模型。為了能夠獲得較為精確的流場(chǎng)參量,考慮可壓縮粘性流動(dòng)的Navier-Stokes方程組。在直角坐標(biāo)系下,其守恒形式為:?U?t+?(FU+FV)?x+?(GU+GV)?y+?(ΗU+ΗV)?z=0(1)?U?t+?(FU+FV)?x+?(GU+GV)?y+?(HU+HV)?z=0(1)式中,U為守恒變量;FU、GU、HU為對(duì)流通量;FV、GV、HV為粘性通量。3.2類tvd模型本文采用Jameson中心差分格式求解上述控制方程。人工粘性項(xiàng)采用二階與四階差分的組合形式;針對(duì)高超聲速流動(dòng)問題,引入類TVD變量修正;為了加速求解過程,采用當(dāng)?shù)貢r(shí)間步長(zhǎng)和隱式殘值光順等加速收斂措施。湍流模型采用Baldwin-Lomax代數(shù)模型。4輻射場(chǎng)數(shù)值法4.1熱流密度qkr的基本原理氣體分子遵循Rayleigh散射,一般工程中,其散射能量很小,故忽略不計(jì)?；诓ǘ桅う藘?nèi)的發(fā)射、吸收性非灰介質(zhì)內(nèi)輻射傳輸方程表達(dá)式為:dΙk(→r,?s)ds=-κk(→r)Ιk(→r,?s)+κk(→r)Ιbk(→r)(2)dIk(r?,s?)ds=?κk(r?)Ik(r?,s?)+κk(r?)Ibk(r?)(2)式中,下標(biāo)k表示波段Δλ區(qū)域(如3～5μm);?ss?為單位方向矢量;→rr?為空間位置矢量;Ik(→rr?,?s)s?)為位置→rr?處、?ss?方向上Δλ內(nèi)的輻射強(qiáng)度;Ibk(→rr?)為位置→rr?處黑體輻射強(qiáng)度;κk(→r)為位置→r處介質(zhì)Δλ內(nèi)吸收系數(shù)平均當(dāng)量。在已知位置→r處各方向輻射強(qiáng)度值后,位置→r處的輻射熱流密度→qK(→r)的表達(dá)式為:→qk(→r)=∫Ω=4πΙk(→r,?s)?sdΩ(3)其中,Ω為空間立體角。認(rèn)為飛行器壁面為不透明、漫發(fā)射、漫反射面,其邊界條件為:Ιk,w(→rw,?s)=εk,wΙbk,w(→rw)+1-εk,wπ∫nw?si＜0Ιk?w(→rw??si)|?nw??si|dΩi(4)式中,εk,w為壁面Δλ波段內(nèi)發(fā)射率平均當(dāng)量;?nw為壁面單位法向矢量。4.2輔助網(wǎng)格單元在求解RTE的方法中,有限體積法是一種非常強(qiáng)有力的方法,它可確保輻射能量整體守恒,對(duì)不規(guī)則邊界適應(yīng)性強(qiáng),可以使用與流場(chǎng)相同的計(jì)算網(wǎng)格,方便地與流場(chǎng)求解的程序相結(jié)合。因此本文采用FVM求解輻射傳輸方程。采用FVM求解輻射傳輸方程時(shí),需要進(jìn)行空間離散和角度離散?？臻g離散與流場(chǎng)計(jì)算相同;角度離散就是將整個(gè)4π空間劃分為互不重疊的若干微元立體角,本文采用文獻(xiàn)提出的角度離散方法,如圖2所示。在控制體Vp和微元立體角Ωm(如圖3所示)內(nèi)對(duì)輻射傳輸方程積分得到其在→sm方向的離散形式為:∑j=e,w,s,n,t,b[Ιmk,c,jSjDmj]≈[-κk,ΡΙmk,Ρ+κk,ΡΙbk,Ρ]ΩmVΡ(5)式中,Ωm為第m個(gè)微元立體角;Imk,c,j為控制體j面上Ωm方向的k區(qū)域內(nèi)輻射強(qiáng)度;Imk,Ρ為控制體中心P處Ωm方向的k區(qū)域內(nèi)輻射強(qiáng)度;Ibk,P為中心P處k區(qū)域內(nèi)黑體輻射強(qiáng)度;Sj為控制體j面的面積;Dmj為控制體j面的方向權(quán)值,其表達(dá)式為:Dmj=∫Ωm(?sm??nj)dΩm(6)式中,→sm表示微元立體角Ωm中心的單位方向矢量;→nj表示表面Sj的單位外法矢量。為了將控制體表面上的輻射強(qiáng)度和體內(nèi)節(jié)點(diǎn)處的輻射強(qiáng)度聯(lián)系起來,本文選用階梯格式(thestepscheme),即:Imk,c,j=Imk,p(7)則方程(5)可寫成如下形式:amΡΙmk,Ρ+∑J=E,W,S,Ν,Τ,BamJΙmk,J=bmk,Ρ(8)式中各項(xiàng)系數(shù)分別為:amΡ=∑j=e,w,s,n,t,bmax[SjDmj,0]+κk,ΡVΡΩmamJ=min[SjDmj,0]bmk,p=κk,pBk,ΤpσΤ4pπVpΩm式中,下標(biāo)J=E,W,S,N,T,B分別表示與控制體VP相鄰的六個(gè)控制體中心節(jié)點(diǎn);Bk,TP表示溫度TP下Δλ波段內(nèi)的黑體輻射能量占其總輻射能量的份額。壁面邊界條件(4)的離散形式為:Ιmk,w=εk,wBk,ΤwσΤ4wπ+1-εk,wπ∑Dlw＜0Ιlk,w|Dlw|(9)對(duì)于“假想面”邊界的處理,不考慮環(huán)境輻射源影響,則就沒有外部輻射能量進(jìn)入封閉系統(tǒng)內(nèi)部,假設(shè)系統(tǒng)通過“假想面”輻射到外界去的輻射能量,將會(huì)被為外界大氣完全吸收掉,不再反射回系統(tǒng)內(nèi)部,也即系統(tǒng)內(nèi)輻射能“只出不進(jìn)”。為了實(shí)現(xiàn)這一點(diǎn),引入一層輔助網(wǎng)格單元作為最外層邊界,如圖4所示。設(shè)DmΡ為系統(tǒng)內(nèi)部P處控制體“單元交界面”上的方向權(quán)值,Imk,Ρ為P處的輻射強(qiáng)度,Imk,Μ為M處的輻射強(qiáng)度,則當(dāng)DmΡ>0時(shí),Imk,Μ=Imk,Ρ;當(dāng)DmΡ≤0時(shí),Imk,Μ保持原值不變。對(duì)于得到的離散方程組,本文采用CGSTAB算法求解。4.3氣體輻射吸收系數(shù)計(jì)算在求解輻射傳輸方程之前需要先知道氣體的輻射吸收系數(shù)。把每個(gè)控制體內(nèi)的氣體看作是等物性的,在已知其熱力學(xué)參數(shù)的基礎(chǔ)上,利用逐線計(jì)算法計(jì)算氣體輻射吸收系數(shù),逐個(gè)計(jì)算所有控制體內(nèi)氣體的吸收系數(shù)即得封閉系統(tǒng)內(nèi)氣體吸收系數(shù)的空間分布。在計(jì)算中考慮光譜線的壓力增寬和多普勒增寬混合效應(yīng),所需譜線參數(shù)由HITRAN數(shù)據(jù)庫(kù)獲得。對(duì)于波段Δλ區(qū)域的吸收系數(shù)平均當(dāng)量,本文采用普朗克平均法計(jì)算,即:κk=∫λ2λ1κλΙbλdλ∫λ2λ1Ιbλdλ(10)此外,假設(shè)組成空氣的各組分氣體之間互不影響,空氣介質(zhì)的吸收系數(shù)等于各組分氣體吸收系數(shù)之和。5計(jì)算結(jié)果和分析5.1計(jì)算程序的初步驗(yàn)證1壁溫ma模型為鈍錐體:球頭半徑為27.94mm,半錐角為15°,底部圓半徑為152.4mm。計(jì)算條件:Ma∞=10.6,P∞=131.97pa,T∞=47.34K,給定壁溫為Tw=294.44K。圖5為零攻角時(shí)鈍錐體表面熱流密度分布與文獻(xiàn)提供的試驗(yàn)數(shù)據(jù)的比較,計(jì)算值與試驗(yàn)值符合得很好(圖中q0為駐點(diǎn)熱流密度值)。2輻射熱流密度q為了在一定程度上驗(yàn)證RTE計(jì)算方法和程序的正確性和可靠性,設(shè)計(jì)了一個(gè)簡(jiǎn)單球形黑體模型:球體半徑為1m,表面為黑體,溫度為500K;周圍為透明介質(zhì)。設(shè)R為到球體中心的距離,圖6表示輻射熱流密度q隨R的變化情況。從圖中可以看出,計(jì)算結(jié)果比理論值稍低些,最大相對(duì)誤差約為30%,這主要是由空間離散和角度離散引起的,但計(jì)算結(jié)果與理論值的變化趨勢(shì)是一致的,完全能反映出球形黑體周圍空間熱輻射的分布規(guī)律,這表明方法和程序是正確可靠的。5.2空氣流場(chǎng)模擬結(jié)果乘波構(gòu)型被認(rèn)為是高超聲速飛行器的最佳構(gòu)型。以乘波構(gòu)型為實(shí)例,利用上述建立的方法和計(jì)算程序數(shù)值仿真了其某一飛行狀態(tài)下的3～5μm波段內(nèi)紅外輻射特性的空間分布?；诟叱曀偌?xì)長(zhǎng)圓錐繞流的近似解,根據(jù)文獻(xiàn)中的設(shè)計(jì)方法構(gòu)建了簡(jiǎn)單的錐導(dǎo)乘波構(gòu)型?；緟?shù)為:總長(zhǎng)度為L(zhǎng)=5m,設(shè)計(jì)馬赫數(shù)為Ma=6,壓縮角為θc=β-ωK=4°,外形如圖7所示。利用Gridgen軟件進(jìn)行網(wǎng)格劃分?？臻g網(wǎng)格的數(shù)量為121×51×61,在物面和前緣附近進(jìn)行了加密,計(jì)算網(wǎng)格如圖8所示。流場(chǎng)計(jì)算條件:攻角為0°,來流馬赫數(shù)Ma=6,飛行高度H=20km。數(shù)值邊界條件:給定壁面溫度Tw=300K;在入口處取自由來流條件;在出口處采用一階外推;遠(yuǎn)場(chǎng)條件為無反射邊界條件。圖9為乘波構(gòu)型流場(chǎng)壓強(qiáng)和溫度的分布圖,圖中錐形激波清晰明顯;壓縮角θc值大約為5°左右,與設(shè)計(jì)值非常接近。以已知的流場(chǎng)數(shù)據(jù)為基礎(chǔ),計(jì)算介質(zhì)的吸收系數(shù)。在這里只考慮了H2O和CO2兩種組分氣體,設(shè)H2O的體積比為2%,CO2的體積比為0.033%。圖10為乘波構(gòu)型流場(chǎng)空氣介質(zhì)的收系數(shù)分布。在已知流場(chǎng)數(shù)據(jù)和介質(zhì)吸收系數(shù)數(shù)據(jù)的基礎(chǔ)上,計(jì)算乘波構(gòu)型周圍的輻射強(qiáng)度場(chǎng)。在計(jì)算中,假設(shè)飛行器壁面在3～5μm波段內(nèi)平均當(dāng)量發(fā)射率為εk,w=0.8;空間角度離散數(shù)量為N(N+2)=8(8+2)=80。下面給出本次計(jì)算的部分計(jì)算結(jié)果。乘波構(gòu)型在3～5μm波段內(nèi)的輻射熱流密度q(W/m2)的空間分布如圖11所示;圖12和圖13分別為截取的輻射熱流密度空間分布的側(cè)視圖、后視圖。由圖11、12、13可以看出:①乘波構(gòu)型紅外輻射場(chǎng)的空間分布基本上為一“錐形”,這與錐導(dǎo)乘波構(gòu)型流場(chǎng)特征相一致;②在乘波構(gòu)型上、下表面附近的空間區(qū)域,紅外輻射分布強(qiáng)烈,這是由于粘性摩擦使物面附近空氣介質(zhì)具有較高的溫度所致;同時(shí)由于下表面還存在強(qiáng)激波,空氣溫度相對(duì)更高,最高可達(dá)約700K,因而下方的紅外輻射更為強(qiáng)烈,傳輸?shù)靡哺h(yuǎn);③乘波構(gòu)型前向迎頭方向上的紅外輻射比較微弱,很快就與環(huán)境基本上沒有差別。這是符合物理實(shí)際的,因?yàn)樵诟叱曀贍顟B(tài)下,飛行器對(duì)波前空間的流場(chǎng)是沒有影響的。綜合以上分析可知,仿真結(jié)果完全能夠反映出乘波構(gòu)型紅外輻射特性的分布規(guī)律。6輻射模型的建立本文以紅