尹振躍，朱定強(qiáng)，任泓帆

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

0 引言

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

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

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

1 計(jì)算方法

1.1 反向蒙特卡洛方法

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

如圖1所示，面元i或者體元j沿著θk角方向在微元立體角dΩk發(fā)射出來(lái)的能量被面元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ā)射出來(lái)的能量沿著θ角方向在立體角dΩ內(nèi)進(jìn)入面元0的能量的比例因子。

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

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

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

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

(3)

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

(4)

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

BMCM可得到輻射能量

(5)

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

1.2 遮擋算法

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

(6)

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

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

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

2 計(jì)算結(jié)果分析

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

采用《傳熱學(xué)》[10]標(biāo)準(zhǔn)算例9-8對(duì)所編寫(xiě)的反向蒙特卡洛程序進(jìn)行驗(yàn)證計(jì)算。該尾焰流場(chǎng)以x軸正方向?yàn)槲惭鎳娏鞣较?。參?shù)如表1所示。

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

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

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

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

在圓柱流場(chǎng)模型x=0 m，計(jì)算得到尾焰流場(chǎng)底面輻射強(qiáng)度如圖4所示。

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

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

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

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

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

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

圖5 尾焰流場(chǎng)溫度云圖Fig.5 Temperature cloud map of flame flow field

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

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

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

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

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

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

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

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

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

遮擋類別進(jìn)入流場(chǎng)射線數(shù)輻射強(qiáng)度無(wú)遮擋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

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

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

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

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

隨著圓柱圓心y坐標(biāo)增大，其距離探測(cè)面越來(lái)越近，輻射熱流密度數(shù)值逐漸減小，可見(jiàn)遮擋效應(yīng)增強(qiáng)。遮擋的y坐標(biāo)4～6.5 m段，即遮擋距離探測(cè)面3～0.5 m段，輻射強(qiáng)度出現(xiàn)驟降，該現(xiàn)象原因是：① 輻射強(qiáng)度與溫度最為相關(guān)。流場(chǎng)非均勻，其軸線附近是溫度最高區(qū)域，隨著遮擋與探測(cè)面距離越來(lái)越近，探測(cè)面發(fā)射的射線逐漸不能進(jìn)入高溫區(qū)，輻射強(qiáng)度出現(xiàn)驟降；② 距離對(duì)輻射驟降的影響不是線性的。以圖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中三角點(diǎn)曲線，可見(jiàn)隨著y坐標(biāo)增大，遮擋與探測(cè)面越來(lái)越近，即n值減少，遮擋半角增幅逐漸加大，也就是遮擋效應(yīng)驟增，與輻射強(qiáng)度下降趨勢(shì)相反。直徑截面規(guī)律與上述相似。上述兩種因素共同作用，導(dǎo)致了輻射強(qiáng)度—y坐標(biāo)曲線出現(xiàn)驟降。

2.3 火箭發(fā)動(dòng)機(jī)推力室上方紅外輻射計(jì)算

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

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

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

3個(gè)探測(cè)面與發(fā)動(dòng)機(jī)的相對(duì)位置如圖12所示。探測(cè)面的法向量均是指向發(fā)動(dòng)機(jī)的。

圖12 3個(gè)探測(cè)面與發(fā)動(dòng)機(jī)的相對(duì)位置Fig.12 Relative position of three detection surfaces and engine

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

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

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

表6 探測(cè)點(diǎn)輻射強(qiáng)度數(shù)值

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

3 結(jié)論

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

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

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

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

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