李志松，馬昌健，毛云杰，郭 濤，汪 行，牟旭娜

（1. 上海衛(wèi)星工程研究所，上海 201109； 2. 上海航天技術研究院，上海 201109）

0 引言

微納衛(wèi)星具有體積小、重量輕、結構簡單等優(yōu)點，其研發(fā)周期短，經(jīng)濟成本低。隨著外形尺寸的縮小，衛(wèi)星的體積比表面積縮減得更快；如果器件的密度變化不大，則微納衛(wèi)星自身的熱容和熱慣性將隨著整星質量的下降而減小。因此，微納衛(wèi)星比普通大中型衛(wèi)星更容易受到內、外熱流的影響而出現(xiàn)溫度波動[1]。

傳統(tǒng)的衛(wèi)星熱設計方法從平臺整體到局部單機的建模和仿真分析，細節(jié)繁多而復雜，并不適用于微納衛(wèi)星快速設計的思路。因此，有必要針對微納衛(wèi)星開發(fā)新的簡化熱分析方法。Tsai[2]總結了一系列的衛(wèi)星熱分析數(shù)學模型，將現(xiàn)有的熱分析模型劃分為通用型和簡化型：通用型的分析模型過于復雜，難以應用于概念設計、迭代優(yōu)化和地面測試模擬等階段；而簡化型的分析模型僅考慮輻射傳熱，精度無法滿足工程設計的需要。Totani 等[3]綜合了輻射熱傳遞和單個零部件內的傳導熱傳遞兩方面因素，提出了單節(jié)點和雙節(jié)點的微納衛(wèi)星簡化熱分析數(shù)學模型，且計算結果對比表明，雙節(jié)點模型能達到與多節(jié)點模型類似的計算精度。Reiss 等[4]指出，對于目前高度模塊化的微納衛(wèi)星，熱分析模型對于不同的微納衛(wèi)星應該具備大部分重用的能力，而無須對每個型號都建立一個新的模型。此外，Elhady[5]在其微型衛(wèi)星設計研究中指出，接觸熱傳導也是熱控系統(tǒng)設計中重要的影響因素。

在現(xiàn)有的熱分析商業(yè)軟件中，基于蒙特卡羅光線追蹤方法的角系數(shù)計算消耗了很大的建模和運算工作量，而微納衛(wèi)星上規(guī)則分布的電路板疊層有利于簡化星內熱輻射中角系數(shù)的計算，可大幅提高建模和仿真的效率。

本文將綜合內、外熱流，內、外輻射及接觸導熱等因素，以集總參數(shù)法對微納衛(wèi)星上的各主要部件進行建模，通過有限復雜度的熱網(wǎng)絡節(jié)點的構建來完成整星的熱分析，以期快速實現(xiàn)熱控設計，并確保模型具有較好的通用性。

1 理論方法

1.1 軌道參數(shù)分析

在衛(wèi)星外熱流的計算中，為確定衛(wèi)星的光照條件，首先要確定衛(wèi)星的運行軌道情況。以軌道6 根數(shù)（半長軸a、偏心率e、傾角i、升交點赤經(jīng)Ω、近地點幅角ω、真近點角θ）作為輸入條件[6]，需要求得衛(wèi)星軌道的周期T0，運行周期內的真近點角θ及其對應時刻t，每個軌道周期內衛(wèi)星進、出地影時的真近點角θ1、θ2及其對應時刻t1、t2，以及衛(wèi)星-地心連線與太陽光線之間的夾角ρ。再結合衛(wèi)星的形狀和姿態(tài)信息，就可以確定星上各表面的外熱流情況。在短期的分析里，不考慮軌道升交點赤經(jīng)和近地點幅角的攝動變化，也不考慮軌道傾角和偏心率的攝動變化。根據(jù)開普勒第二定律，設衛(wèi)星過近地點的時刻為t0、軌道周期為T0，可求得在半個周期內任一時刻t與真近點角θ之間的對應關系為

計算衛(wèi)星表面的太陽外熱流時，需要確定太陽光照與衛(wèi)星運行軌道之間的關系，2 個關鍵的參數(shù)是太陽光線與軌道面的最小夾角β（可以看作衛(wèi)星軌道平面與黃道平面之間的夾角）和軌道近地點與日-地連線之間的夾角Λ。Λ作為軌道近地點到會日點的地心角距，Λ-ω相當于會日點從升交點起量得的地心角距（幅角），Λ+θ相當于衛(wèi)星在軌道面上與會日點之間的地心角距，如圖1 所示。

圖 1 Λ 角示意Fig. 1 The Λ angle

對Λ值的求解，可以利用Cunningham 的坐標轉換方法[7]。相對于以赤道平面為XY平面、春分點方向為+X軸、地球北極方向為+Z軸的不轉動的地球坐標系，設太陽在給定時刻的等效位置坐標為[XS,YS,ZS]，根據(jù)太陽所在的赤經(jīng)ΩS和赤緯δS，可以算得：

可以確定在該地球坐標系中，日-地連線的單位矢量為[SX,SY,SZ]，其中：

圖 2 坐標變換示意Fig. 2 The coordinate transformation

利用Ω、ω和i所確定的空間對應關系，對XYZ坐標進行轉換，可以得到一個以地心為原點、衛(wèi)星軌道平面為xy平面、地心指向衛(wèi)星近地點方向為+x軸的新的地球坐標系（如圖2 所示）。

在新坐標系中，日-地連線的單位矢量為[Sx,Sy,Sz]。新舊坐標間的轉換關系為

根據(jù)Sx和Sy的值，按具體情況選取Λ的角度為0、π/2、arctan(±Sy/Sx)或π-arctan(±Sy/Sx)。獲得Λ值以后，可以算出在真近點角θ，衛(wèi)星-地球連線與太陽光線之間的夾角ρ（相角，取0≤ρ＜π）為

對衛(wèi)星進、出地影的計算，當Sx和Sy均為0 時，β=π/2，衛(wèi)星在整個軌道周期內均接受太陽外熱流；其他情況下，β范圍設為0～π/2，計算式為

假定地球的陰影區(qū)為一圓柱體空間，衛(wèi)星軌道平面與該圓柱體空間相交的區(qū)域為半個橢圓面，則該橢圓的半短軸為地球半徑RE、半長軸為RE/sinβ，不考慮衛(wèi)星軌道非常高和位于半陰影區(qū)時的情況。衛(wèi)星進、出地影的具體時刻主要通過Sx、Sy和衛(wèi)星所在位置的真近點角θ來判斷[8]，當陰影函數(shù)

且其中的(Sxcosθ+Sysinθ)＞0 時，即認為衛(wèi)星處于地影之中。

1.2 外熱流計算

衛(wèi)星外熱流計算主要包括太陽輻射、地球紅外輻射和地球反照，暫不考慮衛(wèi)星自身對外熱流的二次輻射。

在地球附近，太陽輻射強度可以表示為

式中：aE為地球與太陽之間的平均距離，即1 個天文單位（AU）；eE為地球繞太陽運行的軌道偏心率；f為太陽在天球坐標系中的真近點角；CS為太陽常數(shù)。

地球紅外輻射所產生的外熱流為

式中FE為地球與衛(wèi)星之間的角系數(shù)。

僅考慮衛(wèi)星運行在光照區(qū)內和軌道高度不太高的情況，地球反照的熱流密度為

式中：定義(ab)為地球反照系數(shù)，由地球反照角系數(shù)和地球反照率共同決定，計算中取其值為0.35；η為衛(wèi)星-地心連線與太陽光線之間的夾角，與ρ一致，系數(shù)0.9 是考慮了當衛(wèi)星-地心連線與太陽光線之間的夾角超過π/2，但仍然可能受到地球反照的情形。

式(9)、(10)中，F(xiàn)E的近似擬合計算式[9]為

其中：λ為衛(wèi)星的某個表面相對于衛(wèi)星-地心連線（天頂方向）的夾角，例如對地定向的衛(wèi)星，其對地的那一面的λ值恒為π；r為地球半徑RE與衛(wèi)星到地心距離（RE+h）之比，即

S的擬合經(jīng)驗表達式為

根據(jù)以上的理論和近似分析結果，可在MatLab SimuLink 中建立模型，計算外熱流隨衛(wèi)星真近點角的變化。

微小型衛(wèi)星多采用六面體外形，其姿態(tài)控制采用三軸穩(wěn)定并對地定向時，假定+X面為指向飛行方向，+Z面指向地心，+Y面的方向通過右手定則確定。以此衛(wèi)星質心為坐標系原點，衛(wèi)星指向地心的單位矢量為[0, 0, 1]。衛(wèi)星姿態(tài)發(fā)生變化時，假定偏航角為φ，俯仰角為γ，滾動角為φ，其姿態(tài)矩陣為

到達該平面的太陽輻射熱流為Ei=CScosζ。

根據(jù)以上論述，在加入衛(wèi)星的姿態(tài)角信息以后，可以獲得外形為六面體的衛(wèi)星各面上所受的外熱流的值Ei。

1.3 衛(wèi)星溫度場計算

衛(wèi)星溫度場的確定，除了軌道外熱流的計算，還需考慮衛(wèi)星內外表面的吸收與反射，星內熱源以及衛(wèi)星內部的輻射、傳導等耦合傳熱關系。針對小型衛(wèi)星的結構和布局特點，其整體溫度場可采用整星的熱平衡方程進行分析，即

在前面的論述中，外熱流Ei是軌道真近點角θ的函數(shù)。而式(1)中軌道真近點角θ卻是時間t的隱函數(shù)，從任意t求出對應的θ需要采用數(shù)值方法。因此，Ei難以表示為時間t的函數(shù)?，F(xiàn)有的一個解決方法[3]是，先將軌道按真近點角θ從0 到2π平均細分成很多份（如8000～10 000 份）；然后將各θ值與具體的時間點t以及對應的Ei通過列表進行記錄，從而得到離散的Ei(t)；最后對t進行數(shù)值積分，以求出溫度T的變化。這種方法雖然可行，但是操作煩瑣，難以對任意參數(shù)的軌道進行自動求解，且計算精度受到細分份數(shù)的限制。本文提出一個簡單的連續(xù)積分的方法，在式(1)的兩端對t求導，可以求得

從而能夠通過對θ求積分，求出整星溫度T隨真近點角θ變化的情況。通過θ與t之間的連續(xù)對應關系，可獲得溫度T隨時間t變化的連續(xù)曲線。

2 算例

目前大量的微納衛(wèi)星采用了如圖3 所示的印制電路板-機箱式結構布局，內熱源的分布較大中型衛(wèi)星更為簡單和規(guī)律化，使得采用集總參數(shù)法進行簡化分析成為可能。設：微納衛(wèi)星對地定向，偏航角、俯仰角和滾動角均為0，幾何外形為正六面體，尺寸為0.3 m×0.3 m×0.3 m，單個外表面的面積A為0.09 m2，整星質量m為15 kg，整體比熱容c為720 J/(kg·K)，星內熱源Qi為50 W；各表面上鋪設太陽電池的面積比為0 或0.6，太陽電池的太陽吸收比αSs為0.92、發(fā)射率εs為0.85；衛(wèi)星各部件的初始溫度均為20 ℃（293 K）。

圖 3 典型微納衛(wèi)星的內部結構Fig. 3 Typical internal structure of a nanosatellite

設微納衛(wèi)星外殼為鋁合金蜂窩板，厚0.015 m，單個面板面積較小，面內近似等溫。星內熱源主要在印制電路板上，發(fā)熱元器件通過電路連接和導熱硅橡膠等與電路板整體形成熱耦合，故可近似認為電路板本身就是均溫體熱源，從而假定每片電路板就是熱網(wǎng)絡中的一個節(jié)點，并可進一步簡化為沒有厚度的薄平板。各平板相互平行，面積相等，間距遠小于板長和板寬，板與板和板與衛(wèi)星外殼之間的輻射及傳導換熱可通過分析的方法求解。衛(wèi)星的六面體外殼可以進一步分解為6 個單獨的外殼節(jié)點。按電路板堆疊分布的特點，假設前方的外殼板僅與排列在最前方的電路板和四側的外殼板換熱，后方的外殼板僅與排列在最后方的電路板和四側的外殼板換熱。并假設，6 片外殼板之間的相互接觸熱傳導很弱；正前方、正后方的外殼板不與側面的4 片外殼板輻射換熱，同時由于電子元器件的遮擋，4 片側面外殼板之間也不考慮互相輻射換熱。星上各部件之間的輻射換熱，其角系數(shù)可以通過理論公式[10]算出。假設星內有5 片大小相同的電路板，均與X軸垂直，正面均為單一朝向，則整星的傳熱關系可以等效于圖4 所示的部件節(jié)點劃分和圖5 所示的熱網(wǎng)絡節(jié)點及熱阻關系。將整星連同宇宙空間共劃分為16 個節(jié)點和40 個熱阻關系，在SimuLink中進一步構建數(shù)學模型求解。

圖 4 微納衛(wèi)星上等溫化假設的簡化部件Fig. 4 The simplified isothermal components of a nanosatellite

圖 5 多節(jié)點熱網(wǎng)絡節(jié)點及熱阻關系Fig. 5 The multi-nodal network of thermal resistances

圖5 的熱網(wǎng)絡中存在一系列的傳導熱阻需要進行計算評估。衛(wèi)星通過聚酰亞胺（PI）薄膜隔熱，各接觸面用硅橡膠粘貼太陽電池片，故其總傳熱系數(shù)U可以估算為

式中：t為厚度，聚酰亞胺薄膜厚度tPI設為1.0 mm，硅橡膠涂覆厚度tadh設為0.08 mm，太陽電池片厚度tcell設為1.5 mm；k為熱導率，kPI約為0.12 W/(m·K)，kadh約為0.39 W/(m·K)，kcell約為14 W/(m·K)。

據(jù)此估算的傳熱系數(shù)約為1653 W/(m2·K)，對于0.054 m2大小的面積，其熱阻僅約0.011 2 K/W。因此可以認為太陽電池片與外殼是接近等溫的，在分析中可作為同一節(jié)點處理，總節(jié)點數(shù)可由16 個減少到12 個。電路板相鄰間距設為0.05 m，相互之間的導熱連接主要依靠電連接器、支撐柱和在衛(wèi)星殼體上固定點，相鄰兩片電路板之間的傳導熱阻估算為2.32 K/W。電路板與±Y面和±Z面的外殼接觸固定，根據(jù)國外的經(jīng)驗數(shù)據(jù)，每片電路板與外殼面板之間的熱阻設為1.4 K/W。非相鄰的電路板之間，以及電路板與沒有接觸固定的外殼面板之間不考慮直接熱傳導。

各節(jié)點的熱物理性質[11]和建模參數(shù)如表1 所示。表中帶上劃線的為面內平均值。

基于以上設定，對于星內各節(jié)點的溫度變化，可以計算如下：

對于電路板節(jié)點，

3 算例結果與討論

在此，以常見的800 km 高度太陽同步軌道和地球同步軌道為例。兩軌道偏心率均為0.0，太陽同步軌道的傾角為95.597°、半長軸為7178 km，地球同步軌道的傾角為0°、半長軸為42 164.17 km，初始真近點角均取為0°。在春分日，兩個軌道的外熱流分別如圖6(a)、(b)所示。為了對建模計算結果的正確性進行校驗，圖6 還給出采用Thermal Desktop 5.0商業(yè)軟件的仿真數(shù)據(jù)進行對比，假設表面太陽吸收比為1。可以看到，建模計算結果與軟件仿真數(shù)據(jù)在具體的熱流強度和進入地影時刻等方面均符合良好；在最大熱流密度方面兩者稍有偏差，應該是計算時太陽常數(shù)的設置存在差異所致。

圖 6 春分日六面體衛(wèi)星各面所受軌道外熱流隨真近點角變化的情況Fig. 6 The external heat flux vs. the true anomaly for each surface of an orbiting hexahedron satellite

分別對高溫工況和低溫工況進行熱控設計數(shù)值模擬評估，模擬結果如圖7 和圖8 所示。軌道參數(shù)依照前述的太陽同步軌道設置。

由圖7 可見，高溫工況（星上熱源全部啟動）下，由于初始溫度比較接近穩(wěn)態(tài)值，衛(wèi)星在經(jīng)歷了約4 個軌道周期后就基本達到穩(wěn)定的溫度循環(huán)。圖7(a)顯示，各外殼面板之間溫度差異較大。這是由于，±Y和±Z面與星上載荷（電路板）有熱耦合，且面外側做了隔熱多層和噴漆等熱控設計，加之地球紅外輻射和反照的影響，使得這些面的溫度波動較小（13～35 ℃）；而±X面受外熱流影響較大，且與衛(wèi)星內部熱耦合程度不高，故溫度變化劇烈（-10～75 ℃）。圖7(b)顯示，星內各電路板間溫差較小，在衛(wèi)星外殼遮擋和內部熱耦合的作用下，平均溫度約為22 ℃，在1 個軌道周期內的變化范圍僅約15 ℃，有利于電路板上電子設備的可靠工作。

圖 7 高溫工況下衛(wèi)星各部件溫度隨時間的變化Fig. 7 Temperature vs. time for each satellite component under high temperature condition

由圖8 可見，低溫工況（星上熱源全部關閉）下，從初始到穩(wěn)態(tài)的溫度差距較大，衛(wèi)星實現(xiàn)穩(wěn)定溫度循環(huán)所需的時間較長，約為8 個軌道周期。各面板的溫度波動現(xiàn)象與高溫工況類似，但溫度下降，星內平均溫度維持在-25 ℃左右，軌道周期內變化范圍也約15 ℃。通過與Thermal Desktop 5.0 的仿真結果進行對比，可以發(fā)現(xiàn)本研究之建模計算的瞬態(tài)溫度曲線在時間和變化趨勢上與仿真結果吻合良好，溫度值的最大偏差僅約2～3 ℃，表明建模計算的精度已可滿足工程應用的需求。建模計算與仿真結果間的差異主要表現(xiàn)在電路板的溫度預測上，其原因是星內的熱耦合建模較為簡化，一些輻射換熱關系被忽略，不如Thermal Desktop 5.0 中的詳細完整。若模型內繼續(xù)添加熱網(wǎng)絡內各節(jié)點間的換熱關系細節(jié)，則計算精度有望進一步提高。

圖 8 低溫工況下衛(wèi)星各部件溫度隨時間的變化Fig. 8 Temperature vs. time for each satellite component under low temperature condition

4 結束語

本文在分析微納衛(wèi)星熱控系統(tǒng)需求特性和結構布局特點的基礎上，對一般性橢圓軌道衛(wèi)星的外熱流進行建模，并利用集總參數(shù)法建立微納衛(wèi)星的簡化多節(jié)點分析熱阻網(wǎng)絡模型。在建模過程中，突破了衛(wèi)星進入地影、熱平衡方程對真近點角變化的連續(xù)積分和星內熱耦合分析等技術問題。該熱阻網(wǎng)絡模型與真實的微納衛(wèi)星具有較高的熱物理相似性，其對衛(wèi)星的高溫和低溫工況計算示例結果與專業(yè)航天器熱分析軟件的仿真數(shù)據(jù)對比顯示，該模型能有效預測軌道外熱流和衛(wèi)星的瞬態(tài)溫度場，并具有較好的精度，而計算耗時僅為幾秒鐘。該分析模型可以根據(jù)用戶的具體需求選擇不同的軌道參數(shù)，調整衛(wèi)星的外形尺寸、飛行姿態(tài)、電路板數(shù)量、熱負荷工況以及各部件的熱物理性質。在現(xiàn)有模型的基礎上，還可以進一步細化結構，如加入連接各外殼面板和電路板的金屬框架，通過固體導熱提升系統(tǒng)整體的等溫性等，而數(shù)學模型本身的復雜程度不會大幅增加。因此，該模型將是微納衛(wèi)星研制過程中實現(xiàn)熱控設計快速迭代優(yōu)化和方案比較評估的有力工具。