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

尹振躍，朱定強，任泓帆

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

0 引言

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

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

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

1 計算方法

1.1 反向蒙特卡洛方法

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

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

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

(1)

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

(2)

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

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

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

則可計算得到以下關系，即

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

(3)

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

(4)

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

BMCM可得到輻射能量

(5)

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

1.2 遮擋算法

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

(6)

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

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

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

2 計算結果分析

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

采用《傳熱學》[10]標準算例9-8對所編寫的反向蒙特卡洛程序進行驗證計算。該尾焰流場以x軸正方向為尾焰噴流方向。參數如表1所示。

表1 等溫流場參數

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

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

圓柱側面的總輻射熱量為1.5×105W。總輻射熱量與標準驗證算例的數據相對誤差1.06%。從柱面熱流密度分布規律上來看，柱面探測面中心區域出現熱流密度最大值，以柱面中心端面為中心向底面呈對稱趨勢減小，并且相同x截面上的輻射強度數值大小十分相近。

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

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

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

表2 例題與程序結果對比

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

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

基礎流場出自文獻[11]，其溫度云圖如圖5所示。

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

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

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

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

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

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

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

僅增加圓柱類遮擋，模仿真實系統的管、桿等結構。結果如表4所示，遮擋效應明顯。射線被圓柱體遮擋情形如圖8所示，與無遮擋情況對比，圖中間部分形成一個沒有射線的區域，被遮擋之后，射線無法穿透遮擋物進入流場區域。

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

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

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

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

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

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

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

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

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

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

(7)

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

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

本文針對上述的復燃尾焰流場的推力室，計算了該液體火箭發動機推力室的頭部上方若干個位置的輻射。此情形下，推力室作為遮擋物。推力室幾何參數如表5所示。

表5 推力室設計參數

采用結構化網格劃分發動機結構。按照其排布規律，計算出的發動機的子平面法向量是外法向量，滿足遮擋算法要求。計算射線5 000條。

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

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

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

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

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

表6 探測點輻射強度數值

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

3 結論

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

1)流場的紅外輻射會使附近結構表面的熱流密度升高，對于結構的熱環境分析很重要。對于本文中的復燃流場，表面熱流密度可達104W/m2量級。

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

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

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