孫 創(chuàng)，夏新林，戴貴龍

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

0 引言

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

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

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

1 軌道計算方法

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

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

2 環(huán)境映射面

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

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

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

圖2 環(huán)境映射面結(jié)構(gòu)Fig.2 The structure of mapping plane with environment

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

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

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

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

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

2. 2 環(huán)境映射面熱流計算實例

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

2.2.1 太陽輻射熱流

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

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

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

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

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

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

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

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

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

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

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

3 輻射傳遞因子計算

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

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

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

圖4 雙射反射分布函數(shù)Fig.4 Bidirectional reflectance distribution function

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

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

4 算例及結(jié)果分析

根據(jù)上述介紹的環(huán)境熱流計算方法，仍以圖5所示的飛行器外結(jié)構(gòu)為例，分析其不同表面的熱流計算結(jié)果，結(jié)果如圖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附近需經(jīng)過陰影區(qū)，此時不存在環(huán)境熱流。由圖5的太陽照射方向可知，位置1是最容易受到太陽照射的，所以在飛行周期內(nèi)其外熱流均保持較高的值。對于位置3，由于在部分情況下存在遮擋，所以其表面熱流波動較大。而對于位置2，在多數(shù)情況下，其表面熱流僅能通過飛行器其它表面的反射才能得到，所以整體熱流值較低。另外，由于飛行器的部分結(jié)構(gòu)存在轉(zhuǎn)動情況，造成飛行周期內(nèi)各表面的熱流值呈非規(guī)律性變化。

5 結(jié)論

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

［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］.中國空間科學技術(shù)，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］ 安敏杰，程惠爾，李鵬.空間對接機構(gòu)太陽外熱流的計算與分析［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］ 夏新林.空間光學系統(tǒng)的雜散輻射計算與熱分析研究［D］.哈爾濱工業(yè)大學博士學位論文，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］ 夏新林，任德鵬，郭亮.求解介質(zhì)內(nèi)熱輻射傳遞的雙向統(tǒng)計蒙特卡羅法［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］ 談和平，夏新林，劉林華，等.紅外輻射特性與傳輸?shù)臄?shù)值計算—計算熱輻射學［M］.哈爾濱:哈爾濱工業(yè)大學出版社，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.］