誘餌球鏡面反射內壁的輻射傳遞系數計算*

厲夫兵，冷俊敏，李紅蓮

(北京信息科技大學 信息與通信工程學院，北京 100101)

誘餌球鏡面反射內壁的輻射傳遞系數計算*

厲夫兵，冷俊敏，李紅蓮

(北京信息科技大學 信息與通信工程學院，北京 100101)

針對內表面具備鏡面反射特性的空間誘餌球，基于蒙特卡羅方法建立漫發射模型，推導光線的最大反射次數，利用光線在球腔的傳播特性求解光線與球壁的多次反射交點，計算球面微元吸收的光束能量，獲取不同反射率情況下單一球面微元對全部球面微元的輻射傳遞系數。計算結果表明，單一微元對球面各微元的輻射傳遞系數中存在2個峰值，分別為微元對其自身的輻射傳遞系數，以及微元對其球對稱微元的輻射傳遞系數。隨著反射率的增大，微元對其自身的輻射傳遞系數逐漸接近于微元對其球對稱微元的輻射傳遞系數；當反射率接近于1時，兩者的值近似相等。

輻射傳遞系數；輻射換熱；鏡面反射；蒙特卡羅；角系數；球

0 引言

在彈道目標飛行中段釋放氣球誘餌的方法可以降低目標的探測概率[1-2]，增加天基紅外監視系統對目標的識別與跟蹤難度[3-4]。公開文獻表明，空間氣球多采用鋁膜-聚酯薄膜-鋁膜的超薄層結構[5]。在日光照射下，氣球向陽面與背陰面的溫差會加劇氣球內部的輻射換熱，而輻射換熱又反過來影響著目標的溫度分布。因此，提高輻射換熱的計算精度對提高目標溫度的計算準確度具有重要意義。

近年來，國內外許多學者對影響彈頭、誘餌等空間目標的溫度與紅外輻射特性的主要因素進行了研究[6-8]，獲得了目標在日照區及地球陰影區的整體溫度或溫度分布隨時間的變化[9-10]。傳統方法計算目標內部輻射換熱時，通常將目標內壁視為灰體面元[11-12]，利用輻射角系數建立傳熱方程，通過多次迭代求解灰體面元間的輻射換熱，整個求解過程較費時間。張偉清等先利用輻射角系數求解灰體面元間的輻射傳遞系數[13]，再直接利用輻射傳遞系數求解灰體面元間的輻射換熱，避免了直接利用角系數求解輻射換熱時的迭代計算過程，加速了計算進程。值得注意的是，當氣球內壁鋁膜足夠光滑時，氣球內表面將具備良好的鏡面反射特性而非漫反射特性，此時利用基于面元灰體假設的計算方法將會給輻射換熱計算引入新的誤差，因此需要計算鏡反射面元間的輻射傳遞系數。

由于鏡面反射的存在，面元的輻射熱量不可避免地在球腔內多次反射，因此優先采用蒙特卡羅方法計算面元間的輻射傳遞系數。經典的蒙特卡羅方法利用面元的發射率隨機吸射光線，通過統計面元吸收光線數目的方法計算輻射傳遞系數[14-15]，這種隨機吸收光線的方法在面元發射率ε較大時會導致反射光線數量的急劇減少。例如，當系統發射率ε=0.9、面元發射光線數量N=10 000條時，經過系統的1，2，3次反射后，剩余光線數量分別約為1 000，100，10。多次反射發生后，由于剩余光線的樣本數量太少，難以對某面元的反射光線在半空間的分布情況進行分析。通過大量增加發射光線數量的方法雖然能夠緩解上述問題，但同時也引入了大量的碰撞與求交計算，限制了蒙特卡羅方法的計算速度。

本文通過建立鏡反射內壁面元的蒙特卡羅漫發射模型，根據光線的最大反射次數計算每次反射事件發生時光線損失的能量份額，利用光線在球腔內的傳播特性直接求解光線與內壁的多次反射交點，通過統計各次反射事件中面元吸收的光線能量份額計算面元間的輻射傳遞系數，在此基礎上分析鏡反射內壁的輻射傳遞系數分布特性。

1 蒙特卡羅漫發射模型

1.1 目標建模與微元剖分

氣球目標被建模為半徑為R0的球面，并在天頂角和方位角方向進行等間隔剖分，獲得一系列球面微元，并利用微元在天頂角、方位角方向的位序確定微元的位置。由球的性質可知，內壁微元表面任意位置點的法矢量為該點指向球心的單位矢量。

1.2 光線的隨機發射位置、方向與能量

在蒙特卡羅方法中，面元向半空間的輻射能量被視為從面元隨機位置、以隨機方向發射的光線的集合，每束光線都具備一定的能量[14-15]。假定面元向半空間發射的光線總數為N，則每束光線攜帶的初始能量為

(1)

式中：T為面元的瞬態溫度；ρ為面元的反射率；σ為Steven-Boltzmann常數。

在球坐標系中，光線在球面微元內表面的隨機發射位置可表示為

(2)

式中：θa，θb分別是球面微元天頂角的上、下界；φa，φb分別是微元方位角的上、下界；Rθ，Rφ是在[0，1]區間滿足均勻分布的隨機變量。

光線在隨機發射位置點上的隨機發射方向的求解過程可參考文獻[14]。通過分析面元向半空間范圍內投射光束的二維概率密度，可推導光束的二維分布函數并得到光束的發射方向與隨機數的關系為

(3)

式中：θf，φf分別為發射光束的天頂角、方位角；Rθ，Rφ為[0，1]區間內均勻分布的隨機數。

2 輻射傳遞系數的計算

2.1 光線在球腔中的傳播

圖1所示為參考坐標系OXYZ中球面微元A上光束的某次隨機發射位置p(θp，φp)。為計算方便，建立與參考坐標系OXYZ重合的本體坐標系Oxyz，并將本體坐標系Oxyz先沿z軸逆時針旋轉角度φp，再沿y軸順時針旋轉角度π-θp，使該坐標系的z軸負半軸恰好通過光線的發射位置p。因此，從本體坐標系Oxyz向參考坐標系OXYZ的坐標轉換矩陣為

(4)

式中：Ly(β)，Lz(γ)分別為坐標系沿y軸、z軸逆時針旋轉β，γ角度對應的旋轉矩陣，即

圖1 面元隨機發射位置p在坐標系OXYZ與Oxyz中的位置Fig.1 Random position ‘p’ on a spherical facet A in OXYZ coordinates and Oxyz coordinates

如圖2所示，在坐標系Oxyz中，光線從點 p[0，0，-R0] 以角度 (θf，φf) 發射。光線與球壁第1次相交于點d (1次反射點)，并遵循鏡面反射定律繼續向前傳播。由于球內表面上任意點的法矢量都指向球心，故過點d的反射光線必定在由球心、點p和點d所確定的大圓切平面內。容易推得，發射光線的多次反射光線也全部落在該圓內。從圖2中可以看出，發射光線與球壁的1次反射點d在大圓中的圓心角α1=π-2θf，光線與球壁的k次反射點對應的圓心角為

(5)

式中：%表示求余操作。

容易求出，k次反射點在坐標系Oxyz中的天頂角θk和方位角φk為

(6)

將坐標系Oxyz中的k次反射點坐標 (R0，θk，φk) 先轉換為直角坐標，再利用式(4)的坐標轉換矩陣L求反射點在坐標系OXYZ中的直角坐標，并將該直角坐標轉換為球坐標(R0，θr，φr)，最后利用反射點在坐標系OXYZ中的天頂角θr和方位角φr求出反射點所在的球面微元。

圖2 光線在球腔中的多次反射 (坐標系Oxyz)Fig.2 Multi-reflections of a light ray inside the sphere in Oxyz coordinates

2.2 輻射傳遞系數的計算

輻射傳遞系數是離開一個表面到達另一個表面的能量份額，其中包含了從任意面元反射后到達目標表面的能量。因此，球面微元的輻射傳遞系數需要考慮光線的多次反射貢獻。

本文中，單束光線攜帶的初始能量e0是被面元依吸收率逐漸吸收而衰減的。當光線發生第1次反射時，反射光束的能量變為ρe0，被交點所屬微元吸收的能量為(1-ρ)e0。當光線發生第2次反射時，反射光束的能量變為ρ2e0，被交點所屬微元吸收的能量為(1-ρ)ρe0。以此類推，第k次反射時，反射光束的能量變為ρke0，被交點所屬微元吸收的能量為

(7)

給定閾值η，當第k次反射發生時，若反射光束的能量ρke0與初始能量e0之比低于閾值η時視為光束能量被完全吸收，即ρk≤η，則最小反射次數nr的取值可表示為

(8)

式中：「?表示上取整。

綜上所述，對于球面微元的任意一束初始發射光線，利用式(5)～(7)可一次求得其nr個反射點的所屬微元及微元吸收能量。對微元的全部N束發射光線進行遍歷，統計每束光線的全部nr個反射點的所屬微元，并計算微元吸收能量，可求得該微元對其它微元的輻射傳遞系數為

(9)

式中：Ej為微元i的輻射能量中最終被微元j吸收的部分，包括微元i對微元j的直接入射能量和被其它面元反射后到達微元j的能量。

3 數值計算結果與分析

本例的參數設置如下。氣球半徑R0=1 m，內壁鏡面鋁膜的反射率為ρ，球面沿天頂角剖分18份 (dθ=10°)，沿方位角方向剖分36份 (dφ=10°)。對特定的球面微元，隨機發射光線的數目N=100萬條。光束被多次吸收、反射后，當反射光束的能量與初始能量之比低于閾值η=0.01時，輻射能量被視為全部吸收，光線的多次反射計算過程結束。

3.1 輻射角系數

當球內壁為黑體時 (ε=1，ρ=0)，輻射傳遞系數等于輻射角系數。球內壁任意兩球面微元m，n之間的角系數理論值為[15]

(10)

式中：Sn為微元n的面積；θa，θb分別為微元n的天頂角上、下界；φa，φb分別為微元的方位角上、下界。

由式(10)可知，球面微元間的角系數僅與微元n的面積有關，與發射微元m的位置和面積無關。

圖3a)為球面發射微元A(天頂角、方位角中心為 (135°，115°)) 對球內壁全部微元的角系數計算值二維展開圖，垂直方向為球內壁微元的天頂角，水平方向為方位角。圖3b)為3個特定的方位角下 (對應于圖3a)方框) 的輻射角系數計算值與理論值對比，圖3c)為計算值與理論值的相對誤差。

圖3 球面微元的角系數理論值與計算值Fig.3 Theoretical and calculated shape factors for sphere facets

從圖3a)可以看出，該微元對球面各微元的輻射角系數在方位角方向保持相對穩定，僅在天頂角方向變化。從圖3b)～c)可以看出，輻射角系數計算值與理論值很接近。在本例的參數設置下，相對誤差大部分處于±5%的范圍內，最大相對誤差不超過10%。由理論式(10)可知，當球面沿方位角方向均勻剖分時，φb，φa為常數，角系數與方位角無關，由此證明了本文方法計算輻射角系數的正確性。

3.2 輻射傳遞系數

圖4a)～c)分別為鏡面反射率ρ=0.40，0.70，0.99時 (反射次數nr=6，13，459)，發射微元A(135°，115°) 對全體微元的輻射傳遞系數矩陣。從圖中可以看出，與圖3a)的黑體面元不同，鏡面反射時，輻射傳遞系數沿方位角方向也發生變化，并呈現出幾個特性。①隨著反射率的增大，輻射傳遞系數存在兩個較為明顯的峰值區域，一個恰為微元A(135°，115°)自身，另一個是其球心對稱微元B(45°，295°) (注：球面兩點a，b關于球心對稱時，其天頂角滿足θa+θb=180°，方位角滿足 |φa-φb|=180°)。②當反射率較低時，微元A的輻射傳遞系數小于微元B，隨著反射率的增加，其輻射傳遞系數逐漸接近微元B。當反射率接近于1時，兩者的值近似相等。

圖4 發射微元A對球面微元的輻射傳遞系數Fig.4 Radiative heat transfer coefficient between a manua-lly selected facet A and the whole spherical facet

圖5a)～b)為反射率ρ=0.99時 (nr=459)，2個不同位置的發射微元A(5°，55°)，(85°，85°) 對全體球面微元的輻射傳遞系數矩陣。從圖5中可以看出，與圖4c)類似，輻射傳遞系數矩陣中的兩個峰值所在的微元在空間上關于球心對稱。由此可以推知，當反射率較大時，任意位置的發射微元對全體球面微元的輻射傳遞系數均存在2個峰值，且分別位于微元自身及其球心對稱微元。

接下來首先分析圖5中輻射傳遞系數存在2個關于球心對稱的峰值的原因。在圖2所示的坐標系Oxyz中，假設微元A，B分別為發射微元及其球心對稱微元。當反射率ρ增大時，單束光線的反射次數將增加；當ρ接近于1時，反射次數nr將非常大，各次反射點將散布于大圓圓周上。雖然單一光線的多次反射點在大圓圓周上的散布位置是確定的；但對于發射天頂角θf∈[0，π/2]的大量隨機發射光線而言，其多次反射點落在大圓圓周各處的概率將近似相等。由于大圓切割球心對稱微元A，B形成的弧長相等，因此，在方位角φ=φf上，反射點落在微元A上弧線的概率等于其落在微元B上弧線的概率。特別地，由圖2可知，由于任意方位角φ∈[0，2π)的大圓都經過微元A，B，因此反射點落在微元A，B上的概率遠高于其他微元。即反射率接近于1時，微元A，B的輻射傳遞系數近似相等且遠高于其他球面微元。

圖5 ρ=0.99時不同位置的發射微元A對應的輻射傳遞系數矩陣Fig.5 Radiative heat transfer coefficient matrix corresponding to two different manually selected facets A when reflectivity ρ is 0.99

對圖4a)～c)中微元A的輻射傳遞系數隨反射率ρ的增大而逐漸逼近微元B的原因進行分析時，需要對微元A，B吸收的光線能量成份進行研究。通過對微元A的每條發射光線進行追蹤，記錄微元上的各個反射點都源于光線的第幾次反射行為，然后統計微元上各次反射行為發生的次數。圖6所示為在圖4參數設置下，反射率ρ=0.40，0.70，0.83時，微元A，B上光線的第k次反射行為 (k=1，2，…，nr) 發生的次數統計。

圖6 不同反射率時面元A，B上各次反射發生的次數Fig.6 Statistics of multi-reflection of rays on facet A and its spherically symmetric facet B

當η=0.01，ρ=0.4時，由式(8)可知光線發生6次反射即被完全吸收。如圖6a)所示，當k=1時微元A，B上光線的第1次反射行為發生的次數近似相等，即微元A發出的光線中，與微元A，B直接發生碰撞的光線數量近似相等。由于微元A，B是關于球心對稱的微元，面積相等，因此光線與兩面元發生碰撞的概率相等，這與節3.1中輻射角系數的計算結論是一致的。

k=2表示光線的第2次反射行為，即發射微元A發出的光線經過球面微元的1次鏡面反射后再與微元A或B相交。由圖2可知，光線由微元A發出后，只有那些能夠與微元B(或其附近微元) 發生直接碰撞的光線才有可能經1次反射后與微元A相交，此部分光線數量較少。但是，光線由微元A發出后，所有能夠與“赤道”微元 (與Oxy平面相交的球面微元) 發生直接碰撞的光線經過1次鏡面反射后，都有可能與微元B相交。因此，2次反射點落在微元B上的光線數量要遠大于落在微元A上的光線數量，即微元B吸收的光線的2次碰撞能量遠大于微元A，直接導致此時微元B的輻射傳遞系數較微元A大。

由圖6a)～c)可以看出，當反射率ρ增大時，反射次數增多，且光線的第k次反射行為 (k≥3) 發生在微元A，B上的數目趨于相等。此時，雖然微元B吸收的光線的2次碰撞能量仍較微元A高，但相對于總能量而言，其所占的能量份額隨反射率的增大而逐漸減小。因此，如圖4所示，隨著反射率的增大，微元A的輻射傳遞系數逐漸接近于微元B。

在本例中，光線被多次反射后，其攜帶的能量會逐漸降低，若反射光射的能量與初始能量之比低于閾值η=0.01，則光線能量被視為全部吸收，即光束被吸收的能量份額為0.99。當該閾值減小時，計算精度將會提高。圖7所示為η=0.000 1、其他設置與圖5相同時，得到的輻射傳遞系數矩陣。

圖7 η=0.000 1，ρ=0.99時不同位置的微元A對應的輻射傳遞系數矩陣 Fig.7 Radiative heat transfer coefficient matrix corresponding to two different manually selected facets when η=0.000 1，ρ is 0.99

從圖7與圖5可以看出計算結果非常相近。這是因為當閾值η=0.000 1時，光束被吸收的能量份額為0.999 9；而當η=0.01時，光束被吸收的能量份額為0.99。由此可見，兩者的能量變化不足1%，因此圖7得到的結果與圖5類似，但數據精度更高。

4 結束語

本文利用蒙特卡羅方法計算了球內壁的輻射角系數及不同鏡面反射率情況下的輻射傳遞系數。輻射角系數的計算值與理論值一致。輻射傳遞系數的結果表明，首先，微元對全部球面微元的輻射傳遞系數矩陣中存在2個峰值，即微元對其自身、以及微元對其球心對稱微元的輻射傳遞系數。其次，隨著反射率的增大，微元對其自身的輻射傳遞系數逐漸接近于微元對其球心對稱微元的輻射傳遞系數。當反射率接近于1時，兩者的輻射傳遞系數值近似相等。由此可知，對于誘餌球鏡反射內壁微元而言，當某微元的輻射出射能量一定時，反射率越大，輻射能量越傾向于向該微元及其球對稱微元集中。本文的結論有助于分析空間誘餌球向陽面高溫微元與背陰面低溫微元之間的輻射換熱與溫度變化過程。

[1] 酈蘇丹，任萱.大氣層外誘餌釋放研究[J].宇航學報，2001，22(2)：100-104. LI Su-dan，REN Xuan.Research of Releasing Decoy Outside Atmosphere[J].Journal of Astronautics，2001，22(2)：100-104.

[2] 林兩魁，謝愷，徐暉，等.中段彈道目標群的紅外成像仿真研究[J].紅外與毫米波學報，2009，28(3)：218-223. LIN Liang-kui，XIE Kai，XU Hui，et al.Research on Infrared Imaging Simulation of Midcourse Ballistic Object Target Complex[J].Journal of Infrared and Millimeter Waves，2009，28(3)：218-223.

[3] 王碩，張奕群，孫冰巖.紅外點目標跟蹤方法綜述[J].現代防御技術，2016，44(2)：124-134. WANG Shuo，ZHANG Yi-qun，SUN Bing-yan.Review of Point Target Tracking Methods Based on Infrared Sensor[J].Modern Defence Technology，2016，44(2)：124-134.

[4] 浦甲倫，崔乃剛，郭繼峰.天基紅外預警衛星系統及其探測能力分析[J].現代防御技術，2008，36(4)：68-72. PU Jia-lun，CUI Nai-gang，GUO Ji-feng.Space-Based Infrared System and the Analysis of Its Detecting Capability[J].Modern Defence Technology，2008，36(4)：68-72.

[5] SESSLER A M，CORNWALL J M，DIETZ B，et al.Countermeasures：a Technical Evaluation of the Operational Effectiveness of the Planned US National Missile Defense System[R].Cambridge，MA：Union of Concerned Scientists and MIT Security Studies Program，April 2000.

[6] 姚連興，侯秋萍，羅繼強.彈道導彈中段目標表面溫度與紅外突防研究[J].航天電子對抗，2005，21(2)：5-16. YAO Lian-xing，HOU Qiu-ping，LUO Ji-qiang.Research on Surface Temperature and Infrared Signature of the Ballistic Warhead in Midcourse[J].Aerospace Electronic Warefare，2005，21(2)：5-16.

[7] 汪洪源，陳赟.天基空間目標紅外動態輻射特性建模與仿真[J].紅外與激光工程，2016，45(5)：0504002. WANG Hong-yuan，CHEN Yun.Modeling and Simulation of Infrared Dynamic Characteristics of Space-Based Targets[J].Infrared and Laser Engineering，2016，45(5)：0504002.

[8] 敬韓博，王英瑞.臨近空間高超聲速目標天基紅外探測技術研究[J].現代防御技術，2016，44(6)：80-84. JING Han-bo，WANG Ying-rui.Space-Based Infrared Detection for Near Space Hyptrsonic Targets[J].Modern Defence Technology，2016，44(6)：80-84.

[9] ZHU D Q，SHEN W T，CAI G B，et al.Numerical Simulation and Experimental Study of Factors Influencing the Optical Characteristics of a Spatial Target[J].Applied Thermal Engineering，2013(50)：749-762.

[10] 張駿，楊華，凌永順，等.彈道導彈中段彈頭表面溫度場分布理論分析[J].紅外與激光工程，2005，34(5)：583-586. ZHANG Jun，YANG Hua，LING Yong-shun，et al.Theoretical Analysis of Temperature Field on the Surface of Ballistic Missile Warhead in Midcourse[J].Infrared and Laser Engineering，2005，34(5)：583-586.

[11] LI F B，XU X J.Modeling Time-Evolving Infrared Characteristics for Space Object with Mircromotions[J].IEEE Transactions on Aerospace and Electronic Systems，2012，48(4)：3567-3577.

[12] 厲夫兵，許小劍.計算超薄多層介質誘餌溫度的快速算法[J].紅外與激光工程，2011，40(6)：986-991.

LI Fu-bing，XU Xiao-jian.Fast Algorithm for Temperature Calculation of Multi-Layered Ultra-Thin Decoy[J].Infrared and Laser Engineering，2011，40(6)：986-991.

[13] 張偉清，宣益民，韓玉閣.單元表面間輻射傳遞系數的新型計算方法[J].宇航學報，2005，26(1)：77-81. ZHANG Wei-qing，XUAN Yi-min，HAN Yu-ge.A New Method of Radiative Transfer Coefficient between Unit Surfaces[J].Journal of Astronautics，2005，26(1)：77-81.

[14] HOWELL J R，PERLMUTTER M.Monte Carlo Solution of Thermal Transfer Through Radiant Media between Gray Walls[J].Journal of Heat Transfer，1964，86(1)：116-122.

[15] 卞伯繪.輻射換熱的分析與計算[M].北京：清華大學出版社，1988：81-83，144-145. BIAN Be-hui.Analysis and Calculation of Radiative Heat Transfer[M].Beijing:Tsinghua Universtiy Press,1988:81-83,144-145.

Radiative Heat Transfer Coefficient Calculation of Inflatable Sphere with Specular Inner Surface

LI Fu-bing，LENG Jun-min，LI Hong-lian

(Beijing Information Science & Technology University，School of Information & Communication Engineering，Beijing 100101，China)

Monte Carlo technique is used to calculate the radiative heat transfer coefficient for an inflatable sphere with specular reflection inner surface. Light rays are emitted diffusely from a spherical facet and being traced. The maximum number of times a ray can be reflected is derived, and the multi-intersections of a ray with the inner side of spherical surface are obtained. Radiative heat transfer coefficients between a manually selected facet and the whole spherical facet are calculated at different reflectivities. The results indicate that there are two peaks in the calculated coefficients, the smaller one is the coefficient between the selected facet and itself (i.e., the self-to-self coefficient), and the larger one is between the selected facet and its spherically symmetric facet (i.e., the self-to-spherical-symmetry coefficient). The self-to-self radiative heat transfer coefficient gradually approaches and finally equals the self-to-spherical-symmetry coefficient as the reflectivity increases to 1.

radiation transfer coefficient；radiation heat transfer；specular reflection；Monte Carlo；shape factor；sphere

2016-12-01；

2017-02-11

北京市自然科學基金資助項目(3164043)

厲夫兵(1982-)，男，山東日照人。講師，博士，主要從事目標紅外輻射特性建模研究。

10.3969/j.issn.1009-086x.2017.03.030

TK124；TP391.9

A

1009-086X(2017)-03-0193-07

通信地址：100101 北京市朝陽區北四環中路35號北京信息科技大學信息與通信工程學院

E-mail：lifubing@bistu.edu.cn