飛行器復雜外結構的環境熱流計算方法

孫 創，夏新林，戴貴龍

(哈爾濱工業大學能源科學與工程學院，哈爾濱150001)

0 引言

對于空間飛行器的熱分析計算而言，真實可靠的環境熱流模擬是獲得精確溫度場計算結果的前提保證。而環境熱流的計算需要飛行器軌道動力學和熱力學結合，是一個跨學科研究。目前，國外的外熱流計算軟件中，NEVADA是基于Monte-Carlo方法和電磁輻射理論模擬輻射換熱過程的商業軟件，STK為軌道分析軟件，I－DEAS利用其軌道分析結果計算環境熱流。以軌道動力學為基礎，在熱力學領域內計算飛行器復雜外結構的環境熱流，可保證環境熱流的準確性，并使數據結果更好的與熱分析計算耦合，提高熱分析計算的精度。

環境熱流求解過程中，輻射傳遞計算采用的方法通常有角系數法［1］和蒙特卡羅法［2－3］，由于需考慮飛行器表面輻射特性、表面間的多次反射和相互遮擋等問題，所以計算主要以蒙特卡羅法為主［4］。同時，在采用蒙特卡羅法的基礎上，為了更合理的對熱流進行統計，趙立新采用了“假想地球等效面”思想對地球紅外熱流進行統計［5］，安敏杰等利用了“太陽窗口”概念對空間對接結構的太陽熱流進行分析［6］。

空間飛行器復雜外結構主要有以下特點:結構眾多、空間位置復雜、設備尺寸差異大、表面輻射特性多樣，這些因素增加了飛行器環境熱流準確計算的難度。本文采用環境映射面結構對空間飛行器進行包覆后，根據飛行器的軌道特性，利用矢量轉換和矢量點積等方法確定不同時刻映射面接收環境熱流的大小及方向，然后通過求解映射面與飛行器表面間的輻射傳遞因子，最終得到飛行器表面的環境熱流。該方法的優勢在于具有一定的通用性，并在求解輻射傳遞過程中考慮飛行器表面的鏡、漫射特性，結構間的遮擋及多次反射問題。計算中采用發射能束份額修正、正反蒙特卡羅法及雙向統計等改進措施［7－8］，解決飛行器結構中不同設備尺寸差異大等問題。

1 軌道計算方法

對于空間運行的飛行器，了解其軌道六要素即可確定軌道面位置及飛行器的空間位置。六要素包括:飛行器軌道面對地球赤道面的傾角i;由春分點到升交點的地心角距αΩ;飛行器軌道半長軸a;飛行器軌道偏心率e;升交點到軌道近地點角距ω;飛行器通過近地點時刻τp。

以地球為中心，圖1給出黃道面，地球赤道及軌道面的空間位置關系。圖中，坐標系XYZ為地球赤道坐標系，XeYeZe為太陽黃道面坐標系，X0Y0Z0為飛行器軌道坐標系。

2 環境映射面

環境映射面的設計原理:映射面對于外部投射的環境熱流全部吸收，并記錄其投射位置、方向和大小。映射面本身并不積累能量，而是在映射面內表面與飛行器結構外表面組成的封閉系統中，按同樣的位置、方向及大小進行輻射換熱模擬計算，以此得到飛行器外表面環境熱流。由映射面組成的封閉系統形狀依飛行器外部結構形狀而定，以圖2為例，映射面封閉系統為長方體結構。為了避免在計算映射面的熱流時考慮遮擋等影響因素，環境映射面多以平面或圓面為主。

圖1 黃道面、赤道及飛行軌道關系圖Fig.1 The relationship among the ecliptic plane，the equator plane and orbit plane

在確定了飛行器的軌道六要素后，根據飛行器的瞬時位置，可確定飛行器外部映射面上的熱流大小及方向。圖3中，隨著飛行器位置的變化，計算中需了解軌道坐標系 X0Y0Z0、飛行器初始坐標系X1Y1Z1及瞬時坐標系XrYrZr之間的轉換關系。下面主要針對太陽輻射熱流、地球紅外及地球反照熱流的計算進行介紹。

圖2 環境映射面結構Fig.2 The structure of mapping plane with environment

2. 1 太陽投射方向的推導

圖3 飛行器瞬時位置及受照示意圖Fig.3 Schematic for illuminated spacecraft

在圖3所示的三個坐標系中，太陽光矢量→■S是在軌道坐標系OX0Y0Z0描述的，而飛行器結構是在瞬時坐標系OXrYrZr中描述的，所以需通過矢量坐標轉換及乘積的方法，將二者在同一坐標系內表達，再進行相關計算。假設飛行器的瞬時坐標系為三軸穩定，經過τ時刻坐標系旋轉的角度為θ=πτ/2T。則經過τ時刻，則飛行器的瞬時坐標系(OXrYrZr)與飛行器初始坐標系(OX1Y1Z1)的轉換矩陣T1為:

由圖3知，飛行器初始時刻坐標系(OX1Y1Z1)到軌道坐標系(OX0Y0Z0)的轉換矩陣T2為:

至此，可計算出飛行器瞬時坐標系中，太陽光矢量與三坐標軸間的夾角余弦，即確定了太陽矢量在飛行器瞬時坐標系中的表達式。

2. 2 環境映射面熱流計算實例

環境熱流的統計計算中，以太陽投射方向較難確定，尤其針對轉移軌道飛行器。在前述的理論基礎上，舉例說明環境映射面的各熱流統計。

2.2.1 太陽輻射熱流

映射面結構及軌道特性如圖2與圖3所示，在計算各映射面的太陽輻射熱流前，需首先判斷太陽光能否直接照射在該映射面上，判斷依據為太陽光矢量與各面法線夾角應大于90°，小于等于180°。

下面以映射面1為例，列出其表面的太陽輻射熱流計算公式。該面的法向量為(0 0 －1)，太陽光在動態坐標系中的夾角余弦已由前式求得，為(cos?1，cos?2，cos?3)，則太陽光與面1 夾角的法向余弦為:

映射面1接收的太陽輻射熱流表達式為:

式中，S為太陽輻射常數，1353W/m2;A1為映射面1的面積，m2。

2.2.2 地球反照及地球紅外熱流

相比于太陽輻射熱流的計算，映射面上的地球反照太陽及地球紅外熱流計算要簡單一些。下面將以文獻［9］為參考，提出本文的處理方法。

假設飛行器飛行過程中三軸穩定，+X軸指向飛行器飛行方向，+Z軸始終指向地心。則計算中僅需考慮映射面3(與地球相對的虛擬面)所接受的地球紅外輻射及地球反照熱流即可。

式中，Eio為將地球作為250K黑體時輻射強度，值為220W/m2，Re為地球半徑，值為6370km，h為飛行器高度。

地球反照太陽輻射熱流記為:

式中，ρ為地球反照率，取為0.35，φ為日地連線與微元面與地心連線夾角。

無論是地球自身紅外輻射還是地球反照太陽，映射面3發射光束的天頂角均有限制，設其最大角度為θ，則:

3 輻射傳遞因子計算

映射面結構中各面的輻射熱流大小及方向統計完成后，采用蒙特卡羅法計算各映射面與飛行器設備表面組成的封閉系統中各面間的輻射傳遞因子，最終獲得飛行器不同位置表面處的軌道外熱流。

計算中，引入雙向反射分布函數(BRDF)來考慮飛行器結構表面的鏡、漫反射特性［10］。

式(9)定義為:在波長λ下，由(θi，Ψi)投射的光束在(θr，Ψr)方向上引起的光譜輻射強度與投射光譜能量之比。另外，Ωi投射立體角，T為表面溫度。

圖4 雙射反射分布函數Fig.4 Bidirectional reflectance distribution function

計算中考慮發射能束的方向性，以圖5所示的飛行器結構為例，在外表面環境熱流計算中，設備表面間的遮擋及多次反射現象明顯，且部分環境熱流通過開口處進入飛行器內部。在由映射面結構與設備表面形成的封閉系統內，太陽輻射模擬能束以一定方向發射，利用自動搜索功能完成能束在各表面間傳遞，避免因外結構復雜為計算帶來的困難。圖中，①線代表了蒙特卡羅法計算中，能束在設備表面間的多次反射情況;② 線代表了部分計算中，能束在飛行器內部傳遞情況;③線代表了部分區域由于結構遮擋，接收不到輻射熱流的情況;④ 線代表地球反照及紅外能束在模擬過程中，其發射天頂角需保證在θ范圍內。

復雜飛行器外結構中，不同設備間的尺寸差距大，而個別較小設備又是熱分析中比較重視的部分，如飛行器外部直徑為0.6mm的管路。為了提高其計算精度，在環境熱流計算中，可以在能束發射過程中，單獨提高該類設備表面的抽樣密度，并通過改進的蒙特卡羅法進行統計計算，以達到精確求解目的。

4 算例及結果分析

根據上述介紹的環境熱流計算方法，仍以圖5所示的飛行器外結構為例，分析其不同表面的熱流計算結果，結果如圖6。圖中曲線1、曲線2和曲線3所對應的位置在圖5中已有標識。

圖5 映射面與飛行器表面間輻射傳遞Fig.5 Radiation transfer between mapping plane and surfaces of spacecraft

圖6 部分位置外熱流Fig.6 Heat flux for some surfaces

圖6 中，飛行器在3000s附近需經過陰影區，此時不存在環境熱流。由圖5的太陽照射方向可知，位置1是最容易受到太陽照射的，所以在飛行周期內其外熱流均保持較高的值。對于位置3，由于在部分情況下存在遮擋，所以其表面熱流波動較大。而對于位置2，在多數情況下，其表面熱流僅能通過飛行器其它表面的反射才能得到，所以整體熱流值較低。另外，由于飛行器的部分結構存在轉動情況，造成飛行周期內各表面的熱流值呈非規律性變化。

5 結論

針對空間飛行器復雜外結構，各種因素給精確計算其表面的環境熱流帶來困難。本文引入環境映射面的概念，并與蒙特卡羅法結合，提出了一種較通用的環境熱流計算方法。環境映射面的引入為輻射熱流的統計帶來方便，并為蒙特卡羅法的實施提供支持。計算中考慮了復雜結構下飛行器表面輻射特性不同，相互遮擋及多次反射問題。同時，蒙特卡羅法中發射能束份額修正、正反蒙特卡羅法及雙向統計等改進措施可為復雜結構的環境熱流計算提供幫助。

［1］ O’Neill R F，Lorentz D R.A vector approach to numerical computation of view factors and application to space heating［J］.AIAA，1983，1－5.

［2］ Howell J R.The monte carlo method in radiative heat transfer［J］.Journal of Heat Transfer，1998，120:547 －560.

［3］ Howell JR，Perlmutter M.Monte carlo solution of thermal transfer trough radiant media between gray walls［J］.Journal of Heat Transfer，1964，86(1):547 －560.

［4］ 翁建華，潘增富，閔桂榮.鏡反射凹面及相互可視表面的軌道外熱流計算［J］.中國空間科學技術，1994，4:7－14.［Weng Jian-hua，Pan Zeng-fu.Calculating orbital heat fluxes of concave surfaces and mutually visible surfaces with specular reflection［J］.Chinese Space Science and Technology，1994，4:7－14.］

［5］ 趙立新.Monte-Carlo法在空間外熱流計算中的應用—空間光學遙感器光學窗口的地球紅外輻射分析計算［J］.光學精密工程，1995，3(6):71－79.［Zhao Li-xin.Earth infrared radiation analysis and calculation for the optical window of high resolution space optical remote senser［J］.Optics and Precision Engineering，1995，3(6):71 －79.］

［6］ 安敏杰，程惠爾，李鵬.空間對接機構太陽外熱流的計算與分析［J］.上海航天，2006，1:22－26.［An Min-jie，Cheng Hui-er，Li Peng.Computation and analysis of solar external heat flux in the docking mechanism［J］.Aerospace Shanghai，2006，1:22 －26.］

［7］ 夏新林.空間光學系統的雜散輻射計算與熱分析研究［D］.哈爾濱工業大學博士學位論文，1997.［Xia Xin-lin.A study of the stray radiation calculation and thermal analysis for spaceborne optical systems［D］.Harbin Institute of Technology，1997:1 －11.］

［8］ 夏新林，任德鵬，郭亮.求解介質內熱輻射傳遞的雙向統計蒙特卡羅法［J］.工程熱物理學報，2006，27(2):21－24.［Xia Xin-lin，Ren De-peng，Guo liang.Bidirectionally statistical monte carlo method for radiative heta transfer in medium［J］.Journal of Engineering Thermophysics，2006，27(2):21 －24.］

［9］ 閔桂榮，郭舜.航天器熱控制(第二版)［M］.北京:科學出版社，1998.［Min Gui-rong，Guo Shun.Spacecraft thermal control［M］.Second Edition.Beijing:Science Press，1998.］

［10］ 談和平，夏新林，劉林華，等.紅外輻射特性與傳輸的數值計算—計算熱輻射學［M］.哈爾濱:哈爾濱工業大學出版社，2006.［Tan He-ping，Xia Xin-lin，Liu Lin-hua，et al.Numerical calculation of infrared radiative transfer［M］.Harbin Press of Harbin Institute of Technology，2006.］