發(fā)動機熱噴流紅外輻射計算與仿真

發(fā)動機熱噴流紅外輻射計算與仿真李建勛，童中翔*，王超哲，童奇，李賀，張志波【摘要】摘要利用求每個小視場視線方向輻射亮度的方法計算噴流紅外輻射的光譜分布。以輻射傳遞方程數(shù)值和形式為基礎(chǔ)，采用Malkmus統(tǒng)計窄譜帶模型和Curtis-Godson(CG)近似求視線方向的輻射強度。采用CFD分析軟件FLUENT模擬流場和組分摩爾分數(shù)分布。建立噴流紅外成像仿真模型，仿真生成了液體火箭發(fā)動機熱噴流紅外圖像。結(jié)果表明，該方法可以很好地分辨出流場的細微結(jié)構(gòu)。該模型也適用于航空發(fā)動機噴流紅外輻射計算與仿真?！酒诳Q】光譜學(xué)與光譜分析【年(卷)，期】2013(033)001【總頁數(shù)】7【關(guān)鍵詞】關(guān)鍵詞熱學(xué)；發(fā)動機；紅外輻射；熱噴流；紅外圖像引言熱噴流是飛機后向主要紅外輻射來源之一，研究其紅外輻射強度和光譜輻射強度，對飛機紅外抑制技術(shù)、空空導(dǎo)彈及彈道導(dǎo)彈預(yù)警系統(tǒng)的紅外探測系統(tǒng)的優(yōu)化設(shè)計有重要意義。鑒于發(fā)動機熱噴流的紅外輻射特性研究的重要軍事意義，國內(nèi)外卜許多研究單位都開展了相關(guān)研究。國外已經(jīng)建立了比較成熟的紅外輻射特性計算軟件［1］。國內(nèi)的學(xué)者也在這方面進行了大量研究：鄧洪偉等［2］進行了發(fā)動機噴管紅外輻射數(shù)值模擬和實驗研究，潘丞雄等［3,4］計算了發(fā)動機尾噴管的紅外輻射特性，宗靖國等［5,6］研究了發(fā)動機排氣系統(tǒng)紅外特征。本工作以發(fā)動機熱噴流為對象，研究熱噴流紅外輻射強度的光譜分布，建立熱噴流紅外輻射特性的通用計算模型，開發(fā)一種發(fā)動機熱噴流紅外輻射計算程序，實現(xiàn)了對紅外熱像儀的紅外圖像的仿真。1建模計算1.1流場計算噴流輻射取決于燃氣的溫度、壓力和CO2，H2O等成分的組分濃度參數(shù)。因此，流場和組分濃度的計算是必要的。對于亞音速圓噴管熱噴流，可以采用CFD和理想湍流射流模型兩種方法。FLUENT商用軟件［7］,是目前功能最全面、適用性最強、使用最廣泛的CFD軟件之一。其求解器基于有限體積法（finitevolumemethod,FVM）,包含豐富的物理模型，可以精確模擬無粘流、層流、湍流、傳熱和傳質(zhì)以及多相流等復(fù)雜的流動現(xiàn)象，為用戶提供了二次開發(fā)接口，可開展從低速、跨音速到超聲速和高超聲速流動模擬。理想湍流射流模型［8］認為亞音速圓噴管熱噴流的截面包括主體段和初始段部分［9］。1.2視線方向氣體輻射亮度計算1.2.1吸收性介質(zhì)輻射傳遞方程設(shè)定向光譜輻射亮度為Lv的輻射能在吸收/發(fā)射性介質(zhì)內(nèi)傳遞，按照布格爾定律，該能量通過厚度為ds的氣體后，在原射線方向上輻射能的衰減量dLv正比于投射量Lv及厚度ds,即如果沿途介質(zhì)處于局部熱力學(xué)平衡之中，根據(jù)基爾霍夫定律，某方向厚度為ds的介質(zhì)沿著該方向發(fā)射的光譜輻射亮度為考慮到吸收和發(fā)射，由式（1）和（2）可知，該方向光譜輻射亮度的變化量是式（3）即為吸收/發(fā)射性介質(zhì)輻射傳遞微分方程。該式等價于X（S）=為光學(xué)厚度。式（4）兩邊乘以exp（X）從0積分到X從而有式（8）表明：定向光譜輻射亮度Lv（0）經(jīng)過光學(xué)厚度為X（S）的吸收和發(fā)射性介質(zhì)后，最終輻射亮度是兩項輻射亮度的疊加，第一項是Lv（0）經(jīng)過光學(xué)厚度X（L）,被衰減后的輻射亮度；第二項是X*處氣體介質(zhì)發(fā)射的輻射亮度Lvb（X*）經(jīng)X-X*光學(xué)厚度的氣體衰減后的光譜輻射亮度。式（8）是一個積分方程，為便于計算，還應(yīng)推導(dǎo)它的數(shù)值和形式。1.2.2輻射傳遞方程數(shù)值和形式對式（1）從Lv（0）到Lv（S）進行積分因為X（S）,經(jīng)計算，化簡得氣體透過率的定義是設(shè)氣體介質(zhì)厚度為L,以觀測端為起點，輻射傳遞的起始位置為m=N,將傳遞路徑分為N小段，如圖1所示。N的值足夠大，因此在每段內(nèi)氣體壓力，溫度和組分濃度可認為是相等的。則（8）式可以寫成式中X（m）表示從1到m這段氣體的光學(xué)厚度。將式（11）代入式（12）,得到吸收/發(fā)射性介質(zhì)輻射傳遞的數(shù)值和形式在該式中，輻射傳遞方向從1到NoTq（m）表示從點m到點N的透過率。數(shù)值和形式利用透過率來求解氣體輻射傳遞方程，便于計算機編程。對光譜輻射亮度在波長域或頻率域進行積分可以得到波段積分輻射強度。1.3譜帶模型利用輻射傳遞方程式（13）計算視線方向光譜輻射亮度，要用到光譜吸收系數(shù)，光譜吸收系數(shù)是各種分子所有譜線在該波長的吸收貢獻之和，與溫度和壓力（濃度）等有關(guān)，需要采用逐線法來計算，計算量非常大。在譜帶間隔內(nèi)，吸收系數(shù)的變化相對普朗克函數(shù)變化很快，將普朗克函數(shù)看作常數(shù)，利用譜帶平均透過率求譜帶平均輻射亮度，計算量將大幅減小，也能得到比較準確計算結(jié)果。譜帶模型按波數(shù)間隔的大小分為窄譜帶模型和寬譜帶模型。窄譜帶模型［10］的波數(shù)間隔一般為5~50cm-1,而寬譜帶模型則可包括整個振動-旋轉(zhuǎn)譜帶。根據(jù)譜線分布規(guī)律的假設(shè)，譜帶模型分為正規(guī)譜帶模型和統(tǒng)計譜帶模型。正規(guī)譜帶模型認為氣體的吸收帶是由一系列等譜線強度、等譜線半寬、等譜線間隔排列的洛倫茲譜線組成。正規(guī)譜帶模型適用于譜線分布規(guī)則的氣體分子譜帶輻射特性計算。統(tǒng)計譜帶模型認為：氣體分子在波數(shù)間隔內(nèi)分布著N條譜線，這些譜線的位置、強度的分布均是隨機的，即認為在△n內(nèi)任意位置上譜線出現(xiàn)的概率是相同的。按照譜線強度分布規(guī)律假設(shè)的類型不同，統(tǒng)計譜帶模型可分為等線強度分布，指數(shù)線強度分布（Goody），指數(shù)尾倒數(shù)線強度分布（Malkmus）的統(tǒng)計模型。Malkmus模型在大氣科學(xué)研究中應(yīng)用廣泛，后來也被氣體高溫輻射計算采用。對多原子分子氣體Malkmus模型被公認為最佳譜帶模型，當譜帶寬度為25cm-1時和逐線計算的吸收率間的最大誤差為10%，如果用最小二乘法從分子譜線庫推導(dǎo)模型參數(shù)，在帶寬為10cm-1、氣體壓力大于0.1個大氣壓時，Malkmus窄帶模型的誤差小于1%。在很多工程問題中，涉及非均勻氣體輻射傳遞。因此，由均勻氣體推導(dǎo)的譜帶模型處理非均勻氣體，要采用近似方法，即按照非均勻路徑的參數(shù)分布計算出〃等效均勻路徑”的參數(shù)，再采用譜帶模型來計算。以輻射傳遞方程數(shù)值和形式為基礎(chǔ)，本文采用Malkmus統(tǒng)計窄譜帶模型［11］和非均勻氣體Curtis-Godson近似［12］來求視線方向的輻射亮度。2噴流紅外成像仿真2.1噴流紅外圖像仿真方法紅外圖像反映探測陣列中每個像素的輻射照度分布，輻射照度與瞬時視場內(nèi)的輻射亮度分布是一樣的，因此只要計算出所有像素對應(yīng)瞬時視場的輻射亮度，并將其歸一化，按照偽彩色編碼關(guān)系，就能模擬紅外熱像儀拍攝噴流形成的熱像。紅外圖像仿真流程如圖2所示。2.2瞬時視場視線輻射傳遞路徑瞬時視場氣體輻射傳遞路徑指的是瞬時視場中心的視線穿過噴流包絡(luò)體的路徑，由于瞬時視場很小，該傳遞路徑的輻射亮度看作是瞬時視場內(nèi)的平均輻射亮度。將CFD的計算域作為輻射源（完全包裹熱噴流的圓柱形包絡(luò)體），瞬時視場輻射傳遞路徑起點與終點坐標由焦平面陣列探測元視線與圓柱包絡(luò)體的表面交點確定。設(shè)紅外熱像儀鏡頭中心和噴流水平中心在同一面上，鏡頭外法線平行于噴流正側(cè)向。紅外熱像儀視場（FOV）為ax&,其空間分辨率是mxn（即探測陣列有m行n列的探測元）。探測器位于焦平面上，設(shè)探測器陣列中心的坐標是（0,0,f）,噴流與熱像儀的相對位置如圖3所示。第i行，第j列的像素中心的坐標是系數(shù)a',b',c‘的表達式是設(shè)從像素中心經(jīng)透鏡中心引出視線的方向向量為（a,b,c）該視線方程也可寫成如下形式其中t是形參。而圓柱包絡(luò)面的方程是聯(lián)立方程得若該方程有兩個實數(shù)根t1，2,將兩個根分別代入式（17）,得到氣體輻射傳遞路徑起點和終點坐標（xA，yA，zA）和（xB，yB，zB）。如果無根，就說明該探測面視線與包絡(luò)面無交點，接收的是背景輻射。2.3輻射量偽彩色編碼設(shè)紅外熱像儀的動態(tài)響應(yīng)范圍是E0min到E0max，探測器焦平面接收到的輻照度變化范圍使探測元輸出在響應(yīng)范圍之內(nèi)，那么以下關(guān)系式應(yīng)成立E0min<Emin<Emax<E0maxo計算探測器陣列各元素的灰色值其中，color_number是灰色級別數(shù)，表示從最強到最弱分成的等份數(shù)。由于計算機能夠顯示的是8bit的視頻格式，所以范圍一般從0到255，可以取color_number為255。下面采用三原色（RGB）模式來研究灰度圖與偽彩色圖生成的方法?？梢姽獍嗽谝欢ǚ秶鷥?nèi)的不同波長的電磁波，但并沒有必要對每一種波長的可見光都用一種顏色來獨立表示。因為幾乎所有的顏色都可以用紅、綠、藍（RGB）這三種顏色的不同強度來組合而得。用I這個無量綱量來表示紅外熱像儀探測陣列各探測元接收的能量強弱的不同，再將I按一定函數(shù)關(guān)系映射成相應(yīng)的彩色，不同的灰度I對應(yīng)不同的彩色，因此圖像的偽彩色生成首先要設(shè)定紅、綠、藍三個變色函數(shù)，對應(yīng)每一個I都有相應(yīng)的RGB輸出，將這種變換關(guān)系稱為顏色編碼［13L3算例3.1液體火箭發(fā)動機噴流3.1.1流場參數(shù)采用文獻［14］中的算例檢驗本算法準確性，文獻中使用aerodyneradiationcode（ARC）計算了一臺推力為1.125x106N的液體火箭發(fā)動機在2-5pm的紅外輻射特性，發(fā)動機氧化劑為N2O4，還原劑為單甲基肼［14］,考慮噴流中H2O和CO2等組分的吸收和發(fā)射。將圓柱面包住噴焰，圓柱的半徑是11m,長度是112m。整個噴流的靜壓保持為1個大氣壓不變，溫度和組分摩爾分數(shù)隨軸向和徑向位置變化，如表1所示。3.1.2計算方法液體火箭發(fā)動機噴流不含固體顆粒，因此認為沒有散射，并且認為它處于化學(xué)和熱力學(xué)平衡狀態(tài)。根據(jù)表中的流場參數(shù)計算正側(cè)面（方位角為90°）噴流紅外卜光譜輻射強度，所用的方法是：首先，建立噴流的圓柱包絡(luò)面，劃分其表面網(wǎng)格，計算出各面元面積Ai和法線方向。其中，cosei是面元發(fā)現(xiàn)與觀測方向夾角的余弦。將光譜輻射強度對波長積分就能計算出某波段的輻射強度（W-Sr-1）。3.1.3計算結(jié)果圖5是文獻［14］計算的液體火箭噴焰沿正側(cè)向的光譜輻射強度，計算過程可參見文獻［14］,圖6是本方法計算結(jié)果。比較兩圖可以看出，本計算光譜輻射強度及隨波長的變化規(guī)按照觀測方向，判斷在觀測方向可見的面元。對于每個可見面元，以面元中心為起點，按照觀測方向弓I出射線與噴流相交得到氣體輻射傳遞路徑，計算出該路徑的輻射亮度日，該輻射亮度的計算方法與1.2節(jié)相同。通過下式計算出該觀測方向的光譜輻射強度［W-(Sr-pm)-1］律與文獻計算的結(jié)果符合較好。但是在2.5-3pm波段和4.6pm處兩種計算結(jié)果存在較大差異，分析原因可能是由于譜帶模型參數(shù)不同所導(dǎo)致的。3.2熱噴流偽紅外圖像仿真3.2.1流場計算設(shè)置計算域如圖7所示。噴管內(nèi)部從渦輪后開始建模，包括內(nèi)外涵、波瓣混合器和中心錐等部件。其尺寸來自某型發(fā)動機的測量結(jié)果。以某型發(fā)動機最大狀態(tài)作為噴管計算工況，其進口參數(shù)如表2所示。內(nèi)外涵進口燃氣為亞音速，設(shè)定為壓力入口，給定進口總壓、靜壓、總溫。出口截面條件由流動狀態(tài)決定，若為亞音速流，則給定反壓，其余變量通過外插給出；若為超音速流，則外插給出所有流動參數(shù)。壁面采用絕熱、無滑移物面條件。對稱邊界采用無流量通過條件，法向速度為0，其他參數(shù)外推。對于外流的入口，采用黎曼邊界條件。空氣和燃氣的組分參數(shù)如表3所示。最大狀態(tài)時噴管喉道半徑R8=29.3cm，出口半徑R9=30.6cm。其中噴流部分長度和半徑為：1217cmx400cm，網(wǎng)格節(jié)點的數(shù)目是151x111(軸向，徑向)。流場參數(shù)預(yù)處理：噴流輻射模擬需要的流場參數(shù)包括靜壓、靜溫、水蒸汽和CO2摩爾分數(shù)。熱噴流網(wǎng)格如圖8所示，劃分網(wǎng)格時做了局部加密處理以適應(yīng)實際流場參數(shù)的變化。這種非均勻網(wǎng)格的流場在使用時很不方便，先對該流場進行預(yù)處理，插值計算出新的均勻網(wǎng)格的流場作為輻射計算的輸入條件。這里新的流場節(jié)點數(shù)目是151x120。3.2.2噴流流場圖9所示為噴流靜壓等高線分布，內(nèi)外涵氣體摻混后經(jīng)噴管加速膨脹后噴管的出口壓力為104608Pa，由于大氣壓力為95930Pa，因此該噴管處于不完全膨脹狀態(tài)，從圖中可以看到膨脹波和壓縮波組。圖10所示為噴管出口局部靜壓等高線，可見由于噴管出口壓力大于外界大氣壓，在出口邊緣處產(chǎn)生兩道膨脹波，氣流經(jīng)過后氣壓下降并與外界氣流壓力相等。這兩道膨脹波相交后在交點處繼續(xù)產(chǎn)生膨脹波與自由邊界相交，氣流經(jīng)過后壓力繼續(xù)下降，以致壓力稍低于外界大氣，這時外界氣體將壓縮氣流，迫使氣流向內(nèi)轉(zhuǎn)折。由于壓力相差較少，超音速氣流向內(nèi)轉(zhuǎn)折時，在轉(zhuǎn)折點處將產(chǎn)生弱壓縮波。這兩道壓縮波相交后，會在交點處再產(chǎn)生兩道壓縮波。此后氣流壓力又稍大于外界大氣壓力，氣流將繼續(xù)膨脹，再度交替產(chǎn)生膨脹波和壓縮波組，但強度逐漸衰減，直至消失。由圖可知，在該工況下噴流在出現(xiàn)第二道膨脹波時，強度已經(jīng)大為衰減。噴流的溫度分布以及H2O,CO2的摩爾分數(shù)分布如圖11—圖13所示。3.2.3紅外仿真圖像紅外熱像儀的總視場(FOV):20°x15°;空間分辨率：320x240;波段：2000~5000cm-1。噴流橫截面的尺寸是800cmx1217cm。與噴流中心在同一水平面，考慮了兩種熱像儀位置：使噴流計算域軸向恰好充滿視場，距噴流中軸Z0=(1217/2)/tan(10/180xpi)=3451cm，熱像儀距噴管出口Y0=608.5cm。(2)距噴流中軸Z0=1300cm，距噴管出口Y0=Z0xtan(10/180xpi)=229.23cm。如圖15和圖16所示，熱像儀離得越近，噴流高溫區(qū)所占的熱像儀像素點越多，越能清楚的描述高溫區(qū)的輻亮度變化。從圖中可以看到，噴管出口(左側(cè))及下游不遠處分別出現(xiàn)了輻亮度明顯大于周邊的球形區(qū)域，而連接兩球形區(qū)域的地方，高輻射亮度區(qū)域向內(nèi)收縮形成頸部。由于該部分噴流的溫度和紅外活性氣體的摩爾分數(shù)幾乎相等，產(chǎn)生該現(xiàn)象的原因是該處存在交替出現(xiàn)的膨脹波和壓縮波，因此導(dǎo)致壓力高低交替變化，在高壓部位與附近低壓區(qū)域紅外活性氣體的摩爾分數(shù)幾乎相同，但是分子密度更大，因此同樣體積分子數(shù)越多，紅外輻射越大。4結(jié)論以發(fā)動機熱噴流為對象，研究了發(fā)動機熱噴流紅外輻射特性和熱噴流偽紅外圖像的生成方法。分析研究得出以下結(jié)論(1)由均勻氣體推導(dǎo)的譜帶模型處理非均勻氣體，要采用近似方法，即按照非均勻路徑的參數(shù)分布計算出〃等效均勻路徑”的參數(shù)，再采用譜帶模型來計算。提出的視線方向光譜輻射強度計算結(jié)果與文獻計算結(jié)果基本一致。但在2.5-3pm波段和4.6pm處本計算結(jié)果與文獻計算結(jié)果存在較大差異，可能是譜帶模型參數(shù)不同導(dǎo)致的。發(fā)動機噴管出口正后向氣流壓力與外界大氣壓力不相等，引起膨脹波和壓縮波組交替出現(xiàn)，但強度逐漸衰減，直至消失。熱噴流偽紅外圖像中，噴管出口及下游不遠處分別出現(xiàn)輻亮度明顯大于周邊的球形區(qū)域，而連接兩球形區(qū)域的地方，高輻射亮度區(qū)域向內(nèi)收縮形成頸部。紅外仿真圖中可以清晰地分辨出流場的細微結(jié)構(gòu)。References[1]Murtagh,DonalP.35thInternationalConferenceonInfrared,Millimeter,andTerahertzWaves,IRMMW-THz2010,Rome,Italy.DocumentNo.20104813432516,2010,9:445.[2]DENGHong-wei,SHAOWan-ren,ZHOUSheng-tian,etal（鄧洪偉，邵萬仁，周勝田，等）.Aeroengine（航空發(fā)動機），2010,36（1）:44.[3]PANCheng-xiong,ZHANGJing-zhou,SHAOWan-ren,etal（潘丞雄，張靖周，邵萬仁，等）.JournalofAerospacePower（航空動力學(xué)報），2010,25（7）:1518.[4]LIUChang-chun,JIHong-hu,LINa,etal（劉常春，吉洪湖，李娜，等）.JournalofEngineeringThermophysics（工程熱物理學(xué)報），2010,31（9）:1567.[5]ZONGJing-guo,ZHANGJian-qi,LIUDe-lian（宗靖國，張建奇，劉德連）.ActaPhotonicaSinica（光子學(xué)報）,2011,40（2）:289.[6]WEIJun-guang,YANGQing-zhen,LIYue-feng（未軍光，楊青真，李岳鋒）.ComputerSimulation（計算機仿真）,2011,28（4）:66.[7]ZHOUZhang-gen,MADe-yi（周章根，馬德毅）.MachineBuildingandAutomation（機械制造與研究），2010,39（1）:61.[8]WUHui-xia