一種圓軌道航天器外熱流通用計算方法

易樺 黃興 江海

(北京空間飛行器總體設計部 空間熱控技術北京市重點實驗室，北京 100094)

航天器在空間軌道運行所受到的空間外熱流對航天器熱設計有重要影響，在航天器熱設計初期一般都必須進行外熱流分析，以決定航天器熱設計頂層方案和航天器熱分析的極端工況。空間外熱流包括太陽直射熱流，行星反照熱流和行星紅外熱流等。由于太陽為平行光，太陽直射熱流只需要考慮行星的遮擋，相對容易分析；但是行星紅外和行星反照需要考慮航天器與行星之間的位置關系，行星反照還需要額外考慮行星表面受照情況，計算過程非常復雜。早期由于計算機軟硬件的限制，國內外主要采用對不同積分域進行積分的方法，對不同輸入參數需進行查圖表來獲得行星紅外和行星反照[1]，但積分域的確定過程十分復雜，而圖表只能給出部分典型參數對應的結果，局限性較大。隨著計算機技術的高速發展，國內陸續提出了通過蒙特卡洛(Monte Carlo)隨機模擬法[2-4]等求行星紅外和行星反照的方法，這些方法求行星反照時還需要將行星中心到太陽矢量與航天器表面法線在局部球坐標系的經度差作為輸入，卻未給出求經度差的通用方法，無法滿足當前日益復雜的工程應用需求。文獻[5-7]應用積分方法，對太陽同步軌道、傾斜軌道在對地姿態的外熱流進行了分析，文獻[8]用平均模擬法對月球紅外輻射進行了精細計算，但這些方法都是針對特定的某種軌道和姿態進行分析，通用性不足，且部分方法仍采用基于相角的簡化方法[9]求行星反照，準確度不夠。

本文基于軌道坐標系到天球坐標系的轉換矩陣，提出了一種圓軌道航天器外熱流通用計算方法，采用解析法計算太陽直射熱流，將行星可視球冠劃分為等面積微元，采用數值積分法求解行星紅外和行星反照，無需求解行星中心到太陽矢量與航天器表面法線在局部球坐標系的經度差，能有效簡化計算過程；并將該方法應用于典型案例上，與熱分析軟件計算結果進行對比。

1 圓軌道航天器外熱流計算方法

首先定義天球坐標系為全局坐標系，軌道坐標系為局部坐標系。定義Ω為升交點赤經，i為軌道傾角，f為真近點角，ω為近地點幅角。

天球坐標系可以按以下步驟旋轉到軌道坐標系：先繞+Z軸旋轉Ω，再繞+X軸旋轉i，再繞+Z軸旋轉f+ω+π/2，最后繞+X軸旋轉-π/2。因此，軌道坐標系到天球坐標系的轉換矩陣為

(1)

1.1 太陽直射熱流計算

光照期(f≤fin或f≥fout)航天器表面太陽直照熱流qSUN為

(2)

式中：S為太陽輻照強度；I為黃赤交角；Φ為太陽黃經；fin和fout為從近地點算起的進陰影和出陰影位置的真近點角，可按文獻[9]求解。

陰影區(fin

1.2 行星紅外和行星反照熱流計算

圖1 行星紅外和行星反照熱流的計算Fig.1 Calculation of planet infrared and albedo heat flux

航天器表面行星紅外熱流qIR和行星反照熱流qAL為[1,4]

(3)

(4)

式中：ρ為行星反照率；α1為行星面積微元法線(行星中心到面積微元的矢量)和面積微元到航天器矢量的夾角；α2為航天器表面法線和航天器到行星面積微元矢量的夾角；η為行星面積微元法線和行星中心到太陽矢量的夾角；L為行星面積微元到航天器的距離。

考慮到隨機方法[4]在微元劃分數量較小的時候相對不準確，本文使用等面積平均數值積分方法，以每個行星可視局部球冠表面微元的面積相等為原則，按θ和φ分別劃分為m1和m2份，θ和φ取在面積微元的中心值。則

(5)

(6)

(7)

因此，航天器表面行星紅外熱流qIR為[1,4]

(8)

航天器表面行星反照熱流qAL為[1,4]

(9)

(10)

(11)

與文獻[1,4]需要將行星中心到太陽矢量與航天器表面法線在局部球坐標系的經度差作為輸入不同，本文提出cosα2和cosη的計算方法如下：

(12)

(13)

式中：γ為航天器到面積微元矢量與航天器到行星中心矢量的夾角。

(14)

需要說明的是，本方法同樣適用于對可視行星球冠進行隨機分割的情況，分割份數越大，計算結果越接近；另外，本方法只適用于凸表面，不適用于凹表面。

2 本方法典型應用案例

本文充分考慮到圓軌道的特性，為進一步簡化計算，定義一種特殊圓軌道如下：i=β(陽光與軌道面的夾角)，Ω=π/2，ω=-π/2，Φ=0。對于標準開普勒圓軌道，只需要求出β角，并將該特殊圓軌道計算出的外熱流結果進行-Λ(會日點到近地點角距)相位調整即可，對于周期平均值，由于和相位無關，標準開普勒圓軌道和特殊圓軌道計算結果相同。

考慮到大部分航天器姿態為對地定向、對日定向、偏航姿態或上述姿態的組合，本文針對上述姿態的特殊圓軌道，以六面體各個面為典型表面進行外熱流分析。

對于不同姿態，六面體各個表面法線方向不同，因此太陽直射熱流和cosα2有所不同，具體如下。

2.1 對地定向姿態

對地定向姿態定義與軌道坐標系相同，即+Z軸對地，+X軸為飛行方向，+Y為軌道面負法線方向。

各個面光照期qSUN和cosα2如表1所示。

表1 對地定向姿態六面體各個面的光照期太陽直射熱流和cosα2Table 1 Solar incident heat flux during illumination and cosα2 value of cube’s six sides in +Z pointing nadir attitude

2.2 對日定向姿態

對日定向姿態定義為：-Z軸對日，+X軸為行星中心到太陽矢量與軌道面法線方向的叉乘，+Y由右手定則確定。對于-Z面，光照期qSUN=S；對于其他面，光照期qSUN=0。各個面的cosα2如表2所示。

表2 對日定向姿態六面體各個面的cosα2Table 2 The cosα2 value of cube’s six sides in -Z pointing sun attitude

2.3 偏航姿態

偏航姿態定義為：+Z軸對地，+Y軸為+Z軸和行星中心到陽光矢量的叉乘，+X根據右手定則確定。按照上述定義求得偏航姿態局部坐標系到對地定向局部坐標系的轉換矩陣為

(15)

因此，各個面光照期qSUN和cosα2如表3所示。

表3 偏航姿態六面體光照期各個面的太陽直射熱流和cosα2Table 3 Solar incident heat flux during illumination and cosα2 value of cube’s six sides in yaw attitude

3 計算結果與驗證

根據上述計算方法，編制了計算程序。本文以β=60°，H=600 km的上述特殊圓軌道為例，取S=1367 W/m2，ρ=0.3，m1=m2=100，將一個周期劃分為24個點，對偏航姿態外熱流進行分析。

-X和±Y面太陽直射熱流始終為0，其他表面的太陽直射熱流見圖2。-Z面行星反照熱流始終為0，其他表面的行星反照熱流計算結果如圖3所示。行星紅外不隨時間變化，其計算結果如表4所示。上述圖2、圖3中，末尾的1和2分別代表采用本文方法編制的程序計算結果和采用蒙特卡洛方法的Thermal Desktop熱分析軟件計算結果。

圖2 太陽直射熱流分析結果對比Fig.2 Comparison of incident solar heat flux

圖3 行星反照熱流分析結果對比Fig.3 Comparison of the albedo heat flux

表4 行星紅外熱流分析結果對比Table 4 Comparison of the infrared heat flux W/m2

本方法計算結果與軟件計算結果的對比表明：

(1)太陽直射熱流絕對誤差不超過0.1 W/m2，相對誤差不超過0.02%；

(2)行星紅外和行星反照熱流絕對誤差不超過2 W/m2，相對誤差不超過3%(絕對值大于2 W/m2時)。

以上對比結果充分說明了本方法的準確性。本方法編制的程序運行時間為秒級，與熱分析軟件相當；此外，與熱分析軟件相比，本程序無需建模，只需要輸入少量參數和待求表面法線單位矢量，操作簡單，準確度和運行效率均能較好地滿足外熱流分析工作的需求。

4 結束語

本文提出了一種圓軌道航天器外熱流通用計算方法，采用解析法計算太陽直射熱流，采用數值積分法求解行星紅外和反照，無需求解行星中心到太陽矢量與航天器表面法線在局部球坐標系的經度差，能有效簡化計算過程。基于此方法，選取一種特殊圓軌道，對常用姿態(+Z對地，-Z對日，基于+Z對地的+X對日偏航)六面體衛星的瞬態外熱流進行分析，并將偏航姿態外熱流分析結果與熱分析軟件計算結果進行了對比。對比結果表明：本方法計算的瞬態外熱流結果與軟件計算結果相差不超過2 W/m2，相對誤差不超過3%(絕對值大于2 W/m2時)，充分說明了本方法的準確性。本方法能快速計算圓軌道航天器常用姿態外熱流，為航天器熱控設計和熱分析提供重要依據，可以應用于各領域航天器，通用性和擴展性好，準確度和運行效率高，能較好地滿足實際工程應用需求。