火箭發(fā)動機有遮擋情況的尾焰紅外輻射計算

尹振躍，朱定強，任泓帆

(北京航空航天大學(xué) 宇航學(xué)院，北京 100083)

0 引言

火箭發(fā)動機尾焰具有高速、高溫、大流量的特點[1]，會產(chǎn)生強烈的紅外輻射特性。對于結(jié)構(gòu)復(fù)雜的火箭發(fā)動機系統(tǒng)，其主發(fā)動機和多個姿控游機會產(chǎn)生多個尾流流場，這些發(fā)動機和渦輪排氣出口組件等結(jié)構(gòu)會對一些位置產(chǎn)生遮擋效應(yīng)，使得該位置的輻射量數(shù)值變小，影響了整體的熱流密度分布[2]。同時，由于實驗測量尾焰紅外輻射強度成本較高且不易實現(xiàn)[3]，故數(shù)值仿真尾焰紅外輻射強度成為了一種成本低并且可靠的研究方式。

運用反向蒙特卡洛方法(Backward Monte Carlo Method, BMCM)的火箭發(fā)動機尾焰流場輻射特性研究起始于上世紀(jì)90年代，取得了一系列成果。H.F.Nelson[4]1992年在計算火箭發(fā)動機尾焰對發(fā)動機底部加熱問題中首次提出了反向蒙特卡洛法計算模型，并通過計算結(jié)論證明BMCM可以很好地計算熱輻射傳遞問題。D.V.Walters和R.O.Backius[5]研究了在非均勻并具有吸收、發(fā)散、散射介質(zhì)中反向蒙特卡洛法的應(yīng)用，其中熱輻射相互作用為反向蒙特卡洛法提供了理論基礎(chǔ)。張濤[6]的博士論文中提到了光線發(fā)射過程中遮擋判斷的基本數(shù)學(xué)算法。以上研究成果為復(fù)雜有遮擋情況下的火箭發(fā)動機尾焰流場的紅外輻射特性研究提供了重要幫助。

目前，針對復(fù)雜狀況下尾焰的輻射特性研究并不完善，因此需要進(jìn)一步研究。本文把遮擋算法的數(shù)學(xué)模型與反向蒙特卡洛輻射計算模型結(jié)合在一起，編寫基于C++語言的計算機程序，研究了有遮擋情況下，火箭發(fā)動機尾焰流場的紅外輻射特性。

1 計算方法

1.1 反向蒙特卡洛方法

反向蒙特卡洛法是終點開始，逆向追蹤光束路徑，從而確定沿路徑發(fā)射的輻射能進(jìn)入探測器的數(shù)量，其主要依據(jù)為輻射傳遞因子的互易性[7]。BMCM在處理小面元或者小立體角接收到能量輻射量問題時非常有效[8]，利用發(fā)射、反射和吸收等傳遞過程的概率模型，統(tǒng)計到達(dá)被測目標(biāo)的各個體元和面元的熱射線數(shù)量[9]。

如圖1所示，面元i或者體元j沿著θk角方向在微元立體角dΩk發(fā)射出來的能量被面元0沿著θ角方向在微元立體角dΩ吸收的能量可以表示為

Qi→0=ΔAiεicosθkdΩkIb(Ti)Di,k0

(1)

Qj→0=ΔVjkjdΩkIb(Tj)Dj,k0

(2)

式中：kj為體元j的吸收系數(shù)；ΔAi為面元i的面積；ΔVj為體元j的體積；輻射傳遞因子Di,k0或者Dj,k0為從面元i或者體元j發(fā)射出來的能量沿著θ角方向在立體角dΩ內(nèi)進(jìn)入面元0的能量的比例因子。

圖1 閉合空間內(nèi)的反向光束追蹤Fig.1 Backward beam tracing in closed space

若T0=Ti或T0=Tj，則面元0與面元i或體元j之間無凈熱交換，即Q0→i=Qi→0或Q0→i=Qi→0。

則可計算得到以下關(guān)系，即

A0ε0cosθdΩD0,jk= ΔAiεicosθkdΩkDi,k0

(3)

ΔA0ε0cosθdΩD0,jk=ΔVjkjdΩkDj,k0

(4)

上述兩式即為輻射傳遞因子的相對性。由于傳遞因子只與幾何形狀及物性參數(shù)有關(guān)，而與工況無關(guān)，該公式在特殊環(huán)境條件下推導(dǎo)出來，卻在所有環(huán)境下都適用。

BMCM可得到輻射能量

(5)

為計算輻射能量，正向蒙特卡洛法需要跟蹤所有面元i和體元j發(fā)射的能量束，而反向蒙特卡洛方法只需要跟蹤面元0發(fā)射的能量束即可，效率遠(yuǎn)高于正向蒙特卡洛法。

1.2 遮擋算法

本文采用逼近效果最好的三角形網(wǎng)格統(tǒng)一表示物體表面。判斷線段O1mO2m與三角形Es的位置關(guān)系。如圖2所示，單元Es的3個頂點為i，j，m,則

(6)

式中：AijO12，AjmO12，AmiO12分別為三角形ijO12，jmO12,，miO12的面積;ζ為判斷閾值，一般取小于10-3的值即可滿足計算精度。若式(6)成立，則線段O1mO2m與遮擋單Es相交。

遮擋判斷時，判斷式(6)是否成立，若存在一個遮擋單元Es對計算單元E1和E2產(chǎn)生遮擋，則兩者之間輻射熱流為0,并且中止其余遮擋單元的遮擋判斷。

圖2 線段與三角形相交Fig.2 Segments intersect with triangles

2 計算結(jié)果分析

2.1 反向蒙特卡洛方法計算紅外輻射正確性

采用《傳熱學(xué)》[10]標(biāo)準(zhǔn)算例9-8對所編寫的反向蒙特卡洛程序進(jìn)行驗證計算。該尾焰流場以x軸正方向為尾焰噴流方向。參數(shù)如表1所示。

表1 等溫流場參數(shù)

計算得到尾焰流場圓柱側(cè)表面的輻射強度，如圖3所示。

圖3 柱面輻射強度結(jié)果Fig.3 Results of cylindrical radiation intensity

圓柱側(cè)面的總輻射熱量為1.5×105W。總輻射熱量與標(biāo)準(zhǔn)驗證算例的數(shù)據(jù)相對誤差1.06%。從柱面熱流密度分布規(guī)律上來看，柱面探測面中心區(qū)域出現(xiàn)熱流密度最大值，以柱面中心端面為中心向底面呈對稱趨勢減小，并且相同x截面上的輻射強度數(shù)值大小十分相近。

在圓柱流場模型x=0 m，計算得到尾焰流場底面輻射強度如圖4所示。

圖4 底面輻射強度結(jié)果Fig.4 Results of bottom radiation intensity

底面的總輻射熱量為1.8×104W，總輻射熱量與標(biāo)準(zhǔn)驗證算例相對誤差5%。從輻射分布規(guī)律上來看，底面中心區(qū)域出現(xiàn)熱流密度最大值，并以此為中心沿徑向呈逐漸降低趨勢，同心圓上的輻射強度數(shù)值大小十分相近，例題與程序計算結(jié)果如表2所示。

表2 例題與程序結(jié)果對比

綜上所述，本文的輻射計算結(jié)果與教材算例相比，輻射強度分布規(guī)律符合基本理論，總輻射熱量的誤差在可接受范圍內(nèi)，可以認(rèn)為本程序的計算結(jié)果是正確的。

2.2 遮擋對復(fù)燃尾焰流場紅外輻射的影響

基礎(chǔ)流場出自文獻(xiàn)[11]，其溫度云圖如圖5所示。

圖5 尾焰流場溫度云圖Fig.5 Temperature cloud map of flame flow field

流場、遮擋物和探測面三者間的相對位置如圖6所示。計算探測面的紅外輻射并記錄射線路徑。結(jié)果中僅顯示進(jìn)入流場區(qū)域的射線，被遮擋和溢出流場等未與流場相交的射線則未顯示。這樣的好處是可以清楚地展示遮擋物體的遮擋效應(yīng)。計算了平板、管、桿等典型的遮擋物結(jié)構(gòu)時的輻射。流場、遮擋和探測面的幾何參數(shù)如表3所示。計算時發(fā)射5 000條射線。

圖6 流場、遮擋面和探測面三者之間相對位置關(guān)系Fig.6 Relative positional relationship among flow field, shading surface and detection surface

類別幾何參數(shù)尺寸及說明流場長/m48 半徑/m3下底面圓心位置原點軸線方向x軸探測面法向量(0,0,-1)計算點坐標(biāo)(20,0,7)板遮擋法向量(0,0,1)中心點坐標(biāo)(20,0,6.1)面積/m20.2×0.4圓柱遮擋中心點坐標(biāo)(20,0,6.1)長度/m2半徑/m0.05軸線方向x軸

無遮擋物體時，射線發(fā)射情形如圖7所示。圖中的圓柱是流場區(qū)域，射線從探測面元外法向出發(fā)，進(jìn)入流場區(qū)域，最終被吸收。外法向與探測平面所形成的半球空間內(nèi)，射線的分布對稱于外法向，除半球空間部分低緯度的射線由于未進(jìn)入流場區(qū)域沒被記錄外，射線充滿了半球空間。符合物理規(guī)律。

圖7 無遮擋時射線發(fā)射情形Fig.7 Radiation emission situation without shielding

增加平板類遮擋，其幾何參數(shù)如表3所示。結(jié)果如表4所示，射線減少18.47%，輻射強度減少24.57%，遮擋效應(yīng)明顯。

僅增加圓柱類遮擋，模仿真實系統(tǒng)的管、桿等結(jié)構(gòu)。結(jié)果如表4所示，遮擋效應(yīng)明顯。射線被圓柱體遮擋情形如圖8所示，與無遮擋情況對比，圖中間部分形成一個沒有射線的區(qū)域，被遮擋之后，射線無法穿透遮擋物進(jìn)入流場區(qū)域。

對于上述的板和圓柱類遮擋，遮擋物的存在阻止了射線進(jìn)入流場中心高溫區(qū)，而輻射強度大小與溫度最為相關(guān)，所以輻射強度減少的幅度要大于射線減少的幅度。

圖8 射線被圓柱體遮擋Fig.8 Rays blocked by the cylinder

遮擋類別進(jìn)入流場射線數(shù)輻射強度無遮擋2 1501.180 6×103(W·m-2)板遮擋1 7538.905 7×102(W·m-2)下降百分?jǐn)?shù)/%18.4724.57圓柱遮擋1 9157.158 6×102(W·m-2)下降百分?jǐn)?shù)/%11.4639.36

針對圓柱類遮擋，僅改變圓柱半徑，所得輻射數(shù)值變化折線圖如圖9所示。隨著圓柱尺寸變大，輻射熱流密度數(shù)值逐漸減小，遮擋效應(yīng)增強。該變化規(guī)律符合理論。

圖9 圓柱遮擋物不同半徑下的輻射強度Fig.9 Radiation intensity at different radius of cylindrical shield

針對圓柱類遮擋，固定探測面和流場位置，僅改變圓柱的位置，使得圓柱與探測面之間的距離變化。所得輻射數(shù)值變化折線圖如圖10中的方塊點曲線。

圖10 圓柱遮擋物不同位置下的輻射強度Fig.10 Radiation intensity at different positions of cylindrical shield

隨著圓柱圓心y坐標(biāo)增大，其距離探測面越來越近，輻射熱流密度數(shù)值逐漸減小，可見遮擋效應(yīng)增強。遮擋的y坐標(biāo)4～6.5 m段，即遮擋距離探測面3～0.5 m段，輻射強度出現(xiàn)驟降，該現(xiàn)象原因是：① 輻射強度與溫度最為相關(guān)。流場非均勻，其軸線附近是溫度最高區(qū)域，隨著遮擋與探測面距離越來越近，探測面發(fā)射的射線逐漸不能進(jìn)入高溫區(qū)，輻射強度出現(xiàn)驟降；② 距離對輻射驟降的影響不是線性的。以圖11遮擋圓截面為例。

圖11 遮擋半角變化Fig.11 Change of shielding half angle

圖11中y坐標(biāo)每次增加Δx,可以得到第n次的遮擋半角

(7)

設(shè)Δθ=θn-θn+1，n取0～12，得到Δθ-y曲線，如圖10中三角點曲線，可見隨著y坐標(biāo)增大，遮擋與探測面越來越近，即n值減少，遮擋半角增幅逐漸加大，也就是遮擋效應(yīng)驟增，與輻射強度下降趨勢相反。直徑截面規(guī)律與上述相似。上述兩種因素共同作用，導(dǎo)致了輻射強度—y坐標(biāo)曲線出現(xiàn)驟降。

2.3 火箭發(fā)動機推力室上方紅外輻射計算

本文針對上述的復(fù)燃尾焰流場的推力室，計算了該液體火箭發(fā)動機推力室的頭部上方若干個位置的輻射。此情形下，推力室作為遮擋物。推力室?guī)缀螀?shù)如表5所示。

表5 推力室設(shè)計參數(shù)

采用結(jié)構(gòu)化網(wǎng)格劃分發(fā)動機結(jié)構(gòu)。按照其排布規(guī)律，計算出的發(fā)動機的子平面法向量是外法向量，滿足遮擋算法要求。計算射線5 000條。

3個探測面與發(fā)動機的相對位置如圖12所示。探測面的法向量均是指向發(fā)動機的。

圖12 3個探測面與發(fā)動機的相對位置Fig.12 Relative position of three detection surfaces and engine

對于探測面，存在發(fā)動機作為遮擋物時，射線發(fā)射情況如圖13所示，可以看出射線與遮擋存在交匯處，射線被吸收，而終止行進(jìn)，圖13(b)發(fā)動機后出現(xiàn)光線缺口；其他進(jìn)入流場的部分當(dāng)滿足公式時終止。

圖13 有遮擋時探測面1發(fā)射射線情形Fig.13 Emitting rays of detecting surface 1 with shielding

探測面1、2和3計算所得輻射結(jié)果如表6所示。

表6 探測點輻射強度數(shù)值

有遮擋時，3個探測面的輻射強度數(shù)值很接近，均在3 600 W/m2；探測面2與探測面3相比，探測面2更接近推力室軸線，被遮擋程度更甚，所以輻射強度下降更多，約為28%，探測面3輻射強度略大；無遮擋時，探測面1輻射強度數(shù)值最大，約為11 180 W/m2，探測面2的輻射強度稍大于探測面3，約為5 100 W/m2和4 490 W/m2。探測面2相比于探測面3距離流場軸線更近，所以輻射強度數(shù)值更大。

3 結(jié)論

本文編寫了一個可以計算流場任意處紅外輻射特性的計算程序，得到了遮擋對于流場紅外輻射特性的影響。獲得如下結(jié)論:

1)流場的紅外輻射會使附近結(jié)構(gòu)表面的熱流密度升高，對于結(jié)構(gòu)的熱環(huán)境分析很重要。對于本文中的復(fù)燃流場，表面熱流密度可達(dá)104W/m2量級。

2)遮擋對于紅外輻射的傳輸起到了阻礙作用。本文圓柱類和板類遮擋的計算結(jié)果，與無遮擋的情況對比，有遮擋的情況下紅外輻射熱流密度降幅可達(dá)50%以上。

3)相同條件下，遮擋物的尺寸越大，對輻射傳輸?shù)南鳒p越大；遮擋物位于探測面和流場中間時，改變遮擋物位置，使其與探測面距離逐漸變小時，輻射強度逐漸減少，并在一定范圍內(nèi)出現(xiàn)驟降。

4)對于本文采用的復(fù)燃流場主流推力室上方的紅外輻射計算，大小約為3 600 W/m2。與無遮擋的情況相比下降范圍為10%～70%。