耿利寅，韓東陽，張傳強，孟恒輝，李國強

（北京空間飛行器總體設(shè)計部，北京 100094）

0 引言

熱分析計算是航天器熱設(shè)計驗證和在軌溫度預(yù)示的重要手段。目前，一般航天器熱設(shè)計中熱分析計算的步驟是：建立航天器在軌或試驗狀態(tài)的熱物理模型，并將其離散化為一系列描述節(jié)點間熱交換過程的傳熱學(xué)方程（組），即熱網(wǎng)絡(luò)方程；再通過對熱網(wǎng)絡(luò)方程的求解，獲得所關(guān)注節(jié)點和時段的瞬態(tài)溫度結(jié)果。熱網(wǎng)絡(luò)方程因其中的輻射換熱項包含溫度的4 次冪而隸屬于典型的非線性方程。對于非線性方程，通常采用迭代法進行求解，而目前尚無完善的理論方法能直接判斷有限次迭代是否可以獲得收斂解。

研究結(jié)果表明，對于一個包含輻射換熱項的確定的熱物理模型，其瞬態(tài)計算的收斂性與時間步長密切相關(guān)：若設(shè)置的時間步長過大，在迭代計算過程中可能引起部分節(jié)點溫度振蕩失穩(wěn)、不能收斂，造成計算失??；反之則使計算效率低下，耗時過長。因此，對瞬態(tài)熱計算收斂性與時間步長及其優(yōu)化選定等開展分析研究，對優(yōu)化建模、提高計算效率，以及指導(dǎo)模型構(gòu)建中的合理簡化、節(jié)點劃分和換熱關(guān)系設(shè)定等具有重要意義。

本文對Gauss-Seidel 迭代法求解熱網(wǎng)絡(luò)方程的收斂性問題進行分析，結(jié)合航天器熱控設(shè)計的工程實際，提出一種近似線性化方法，即，將含有輻射項的非線性熱網(wǎng)絡(luò)方程轉(zhuǎn)化為線性熱擴散方程。在此基礎(chǔ)上進行分析以確定計算收斂的條件和最優(yōu)時間步長，并輔以算例驗證。

1 熱網(wǎng)絡(luò)方程（組）及其節(jié)點分類

運行在宇宙空間中的在軌航天器通常無對流換熱過程，僅需考慮輻射和傳導(dǎo)兩種換熱方式。對于一個由

n

個節(jié)點組成的航天器熱物理模型，描述其能量平衡的熱網(wǎng)絡(luò)方程為

式中：下標

i

,

j

代表節(jié)點編號，

i

,

j=

0, 1, …,

n

?1；

C

為節(jié)點

i

的熱容，J/K；

T

為溫度，K；

t

為時間，s；

Q

為外熱流功率，W；

q

為內(nèi)熱源功率，W；

D

為節(jié)點

i

、

j

間的導(dǎo)熱系數(shù)，W/K；

B

為節(jié)點

i

、

j

間的輻射換熱系數(shù)（GebHart 系數(shù)）；

A

為節(jié)點

i

的輻射換熱面積，m；

ε

為全波長發(fā)射率；

σ

為Stefan-Boltzman 常量。

在實際計算中，通常將一個熱分析模型的節(jié)點分成3 類——擴散節(jié)點、算術(shù)節(jié)點和邊界節(jié)點。其中，邊界節(jié)點溫度是已知的，不需求解。擴散節(jié)點是具有一定熱容，并與其他全部或部分節(jié)點有換熱關(guān)系的節(jié)點。在瞬態(tài)計算中，給定起始溫度，則擴散節(jié)點的溫度可由式(1)進行迭代計算獲得。算術(shù)節(jié)點通常用以模擬覆蓋在實體外的多層隔熱組件、薄膜等熱容非常小的表面，其熱容在模型中設(shè)為0，其能量方程與時間無關(guān)，實際退化為顯式格式，求解不存在收斂性問題。因此，本文主要針對擴散節(jié)點開展討論。

2 求解收斂性與時間步長分析

2.1 方程的線性化近似

考查節(jié)點

i

在迭代過程中的情況，因研究對象為節(jié)點

i

，則首先可假設(shè)節(jié)點

i

以外的其他節(jié)點都可收斂，并已經(jīng)收斂，其溫度為

T

（

j≠i

）。

對于擴散節(jié)點，由于式(1)右側(cè)第4 項的輻射換熱項為非線性項，在計算傳熱學(xué)相關(guān)領(lǐng)域公開發(fā)表的文獻中尚未發(fā)現(xiàn)其迭代求解收斂性的相關(guān)研究結(jié)論；但對于線性擴散方程，收斂性研究的結(jié)論是明確的。為此，本文考慮首先進行方程的線性化。在迭代法求解方程組的過程中，每個節(jié)點的起始溫度都是給定的，為了便于分析，將式(1)中的輻射換熱項簡化為導(dǎo)熱形式，即

式(4)與典型的熱傳導(dǎo)擴散方程形式完全相同，說明通過上述近似轉(zhuǎn)換，已經(jīng)將式(1)的非線性方程轉(zhuǎn)換為線性方程。需要說明的是，上述線性化近似只用于計算收斂性的分析，實際方程求解時，仍須按照式(1)進行迭代計算。

2.2 節(jié)點收斂性分析

在每個特定時間步長的迭代計算中，式(4)中的

Q

、

q

（即內(nèi)、外熱源）為恒定值，由此引起的節(jié)點溫度變化在每次迭代計算中是恒定的。造成節(jié)點溫度不能收斂的因素來源于式(4)中的第3 項，即節(jié)點與其他節(jié)點之間的等效熱導(dǎo)項。下面分析迭代過程中的等效熱導(dǎo)，為求簡便，將研究的節(jié)點定義為節(jié)點0，其他節(jié)點為1, 2, …,

n

?1。

1）第1 次迭代

在此迭代步驟中，節(jié)點0 向其他節(jié)點釋放的熱量為

3）收斂性分析

若節(jié)點0 為擴散節(jié)點，在迭代求解的過程中，其溫度結(jié)果收斂的條件為：每一次迭代，其溫度相對上一步的變化量小于上一步迭代的值。即，若令

K

為相鄰兩步迭代的溫度變化量之比，則有

由式(12)顯見

K

與迭代次數(shù)無關(guān)，說明每一次溫度的迭代結(jié)果與上一次的差值是一個與迭代次數(shù)無關(guān)的等比數(shù)列。因此，|

K

|＜1 是節(jié)點0 溫度結(jié)果收斂的充要條件。同時，對節(jié)點

i

，記

D

+

D

+ …+

D

=

D

，有節(jié)點

i

溫度結(jié)果收斂的充要條件為

式(13)可敘述為：在迭代法求解過程中，某個節(jié)點收斂的充要條件是該節(jié)點與其他所有（與之有換熱關(guān)系的）節(jié)點的等效熱導(dǎo)之和與計算時間步長的積小于該節(jié)點的熱容。上述“等效熱導(dǎo)之和與時間步長的積”的單位為[(W/K)?s]即J/K，與節(jié)點熱容的單位相同。熱容的物理意義為“節(jié)點溫度每變化1 ℃所吸收或釋放的熱量”，即單位溫升熱流量。類似地，“等效熱導(dǎo)與時間步長的積”的物理意義可以理解為節(jié)點與其他節(jié)點溫差為1 ℃時的熱流量，可以將其命名為“單位溫差等效熱流量”。那么，式(13)可簡明地敘述為：某節(jié)點收斂的條件是其與其他節(jié)點單位溫差下的等效熱流量小于該節(jié)點的熱容。

2.3 模型的收斂性與時間步長約束

以上分析是選定模型中的1 個節(jié)點進行的，但分析結(jié)果對全部擴散節(jié)點有效。則，當(dāng)一個熱物理模型的全部擴散節(jié)點均滿足式(13)時，該模型計算結(jié)果收斂。

另外，上述結(jié)論也可用于模型時間步長的選定。對某個節(jié)點

i

，其時間步長應(yīng)滿足

可見，只要采用式(14)對所有擴散節(jié)點計算其允許時間步長，選擇其中的最小值，即為整個模型收斂允許的最大時間步長。

2.4 時間步長優(yōu)化選定

由式(12)可知每一次迭代溫度結(jié)果與上一次的差值是一個與迭代次數(shù)無關(guān)的等比數(shù)列，因此節(jié)點

i

第

m

次迭代的結(jié)果為

進行迭代計算時，需指定迭代的截斷誤差

δ

。截斷誤差通常需就每個實際算例綜合考量而確定：設(shè)置較小的截斷誤差將有助于提高計算精度，但計算量大、收斂過程緩慢、計算耗時長；設(shè)置較大的截斷誤差可以減少計算量、縮短計算時間，但計算精度不高。指定截斷誤差后，當(dāng)

時，認為該節(jié)點溫度滿足收斂條件，停止迭代。

參考式(9)和式(6)，對節(jié)點

i

，

根據(jù)式(9)可知，時間步長越長，

q

的絕對值越接近于1，則式(15)、(16)、(19)清晰表明迭代次數(shù)將增多。因此，時間步長的增大，在加快模型時間推進的同時，會增多單個時間步上的迭代次數(shù)，對計算效率的影響存在2 種趨勢相反的影響，需權(quán)衡優(yōu)化之。式(19)兩邊同時取自然對數(shù)，得到節(jié)點

i

滿足收斂條件的最少迭代次數(shù)為

其中時間步長Δ

t

越長，模型進展越快；同時，迭代的次數(shù)

m

越多，計算量越大，消耗的計算時間越長。因此，Δ

t

/

m

是一個表征計算效率的物理量，其物理意義是完成一定計算量所獲得的模型進展時間，或可更具體地表達為：完成1 次迭代計算獲得的模型進展時間，其單位為s/次。記

當(dāng)研究的物理模型確定時，則

D

、

D

、

C

等均為定值；

T

、

T

為計算的起始溫度，已知并確定；

δ

為計算設(shè)定的截斷誤差，也是確定值。因此，式(21)右側(cè)實際上是一個以Δ

t

為自變量的形如

y

=[

x?

ln(

bx

)]/

a

的函數(shù)。基于以上分析，以Δ

t

為自變量，

t

為因變量，對式(21)求導(dǎo)，可確定

t

取最大值時對應(yīng)的最優(yōu)Δ

t

為

式中e 為自然常數(shù)。

通過式(22)可以確定瞬態(tài)計算中獲得最高計算效率的最優(yōu)時間步長。對確定的熱物理模型，時間步長Δ

t

的選取決定了計算的效率；對所有滿足收斂條件的時間步長Δ

t

來說，存在一個最優(yōu)值，使得時間步長取該值時模型計算的效率最高，即

t

最大，則此時間步長為該模型該步計算的最優(yōu)時間步長Δ

t

。結(jié)合式(14)可以看出，若定義滿足收斂條件的最大允許時間步長為Δ

t

，即Δ

t

=

C

/

D

，則有

即，迭代法求解航天器熱網(wǎng)絡(luò)方程時，最優(yōu)時間步長是計算收斂允許的最大時間步長的1/e。

3 算例驗證

3.1 收斂條件驗證

為驗證上述方法的正確性，選取一個典型的航天器熱物理模型進行分析。該模型用于模擬一個衛(wèi)星的小艙，由艙板、設(shè)備、多層隔熱組件等組成，艙板外側(cè)開設(shè)有散熱面。模型被劃分為40 個節(jié)點，其中邊界節(jié)點2 個、算術(shù)節(jié)點13 個、擴散節(jié)點25 個；這些節(jié)點之間有傳導(dǎo)和輻射換熱關(guān)系，共有傳導(dǎo)項49 個、輻射項151 個。模型置于空間軌道環(huán)境下，其熱源包括外熱流和設(shè)備熱耗。

對所有擴散節(jié)點的傳導(dǎo)和輻射換熱關(guān)系進行統(tǒng)計，發(fā)現(xiàn)其中編號為3015 的節(jié)點，按式(14)計算所得的收斂允許時間步長最小。節(jié)點3015 代表艙板內(nèi)表面的一部分，其表面積

A

=0.04 m，太陽吸收比

α

=0.20，紅外發(fā)射率

ε

=0.87，熱容為66.47 J/K。該節(jié)點與相鄰艙板節(jié)點、同一艙板的外側(cè)節(jié)點以及其上安裝的設(shè)備均有傳導(dǎo)換熱關(guān)系；與模型小艙的其他艙板內(nèi)表面和設(shè)備有輻射換熱關(guān)系。按照上述與其他節(jié)點的傳導(dǎo)和輻射換熱關(guān)系，計算節(jié)點3015 的等效總熱導(dǎo)為

D

=4.578 7 W/K，按照式(14)計算，可得確保計算收斂的最大時間步長Δ

t

=14.52 s?，F(xiàn)設(shè)置截斷誤差

δ

=1.0×10℃，并采用不同時間步長Δ

t

對節(jié)點3015 進行迭代計算，其收斂情況見圖1，圖中曲線表示迭代過程中節(jié)點3015 的溫度結(jié)果變化過程。

圖1 節(jié)點3015 不同時間步長下的溫度迭代計算過程Fig. 1 Iterative calculation processes for node 3015 with different time steps

從圖1 可明顯看出，在時間步長Δ

t

＜Δ

t

時，節(jié)點3015 的溫度迭代過程很快收斂；隨著Δ

t

的增大，結(jié)果收斂需要的迭代次數(shù)逐漸增加；當(dāng)Δ

t

=14.3 s時，結(jié)果收斂需經(jīng)過410 次迭代；當(dāng)Δ

t=

14.5 s（接近Δ

t

）時，經(jīng)過500 次迭代，溫度結(jié)果仍呈振蕩狀態(tài)，不能收斂；當(dāng)Δ

t

＞Δ

t

時，迭代結(jié)果發(fā)散，且Δ

t

越大，發(fā)散速度越快。這樣的數(shù)據(jù)和結(jié)果與式(13)理論分析的結(jié)論高度一致，驗證了理論分析的正確性。

3.2 最優(yōu)時間步長的驗證

仍以3.1 節(jié)所述模型進行驗證。選取不同的時間步長進行迭代計算，結(jié)果如表1 所示。表1 列出了選取不同時間步長Δ

t

時，通過式(20)和式(21)計算得到的理論迭代次數(shù)

m

和理論計算效率

t

，以及實際計算的迭代次數(shù)

m

和計算效率

t

。圖2 以曲線形式較為直觀地對比了

t

和

t

隨時間步長變化的情況。

圖2 實際計算效率與理論計算效率的對比Fig. 2 Comparison between actual calculation efficiency and theoretical calculation efficiency

表1 采用不同時間步長的迭代計算效率對比Table 1 Results of calculation efficiency with different time steps

從表1 和圖2 的結(jié)果可見，計算效率最高時的最優(yōu)時間步長Δ

t

在5.0～5.5 s 之間。這與通過式(22)和式(23)計算所得的Δ

t

=14.52/e=5.34 s的結(jié)果吻合，驗證了2.4 節(jié)關(guān)于最優(yōu)時間步長選定方法的正確性。

4 結(jié)束語

本文采用理論分析的方法對瞬態(tài)熱物理模型迭代求解過程中的收斂性問題進行研究，獲得了節(jié)點和模型收斂的條件；并進一步分析時間步長與計算效率的關(guān)系，提出了以計算效率為優(yōu)化目標的最優(yōu)時間步長確定方法。主要結(jié)論如下：

1）對熱物理模型的某擴散節(jié)點

i

，若其與其他節(jié)點單位溫差下的等效熱流量小于該節(jié)點的熱容，則該節(jié)點溫度在迭代計算中可收斂，否則不收斂；對于整個模型，若其所有擴散節(jié)點滿足上述條件，則模型收斂，否則不收斂。2）以上結(jié)論的另一種形式為，某擴散節(jié)點

i

迭代收斂的允許最大時間步長為該節(jié)點的熱容與等效熱流量之比。對模型的所有擴散節(jié)點逐一計算其確保收斂的允許最大時間步長，選擇其中的最小值，即為保證整個模型收斂允許的最大時間步長。

3）在所有滿足收斂條件的時間步長中，存在一個最優(yōu)值，采用該值時計算的效率最高，即1 次迭代計算獲得的模型時間進展最大。該最優(yōu)時間步長為計算收斂允許的最大時間步長的1/e。

需要說明的是，在結(jié)果收斂的情況下，即模型所有節(jié)點的最終溫度結(jié)果恒定或接近恒定時，時間步長會影響從初始值走向最終穩(wěn)定結(jié)果的瞬態(tài)過程，但不影響最終溫度結(jié)果。實際在軌航天器的溫度往往做周期性波動，在這類航天器的瞬態(tài)熱物理模型中，時間步長不僅影響瞬態(tài)過程，也對最終溫度波動范圍有影響。關(guān)于時間步長對計算精度的影響，受目前研究所限，本文暫不涉及，后續(xù)將對該問題開展深入探討。