李志杰 果琳麗

（中國空間技術(shù)研究院載人航天總體部，北京 100094）

0 引言

與非載人航天器相比，載人航天器均含有密封艙，而密封艙內(nèi)會(huì)有供人生存的大氣存在[1]。因此在對(duì)艙內(nèi)熱環(huán)境進(jìn)行預(yù)測分析計(jì)算時(shí)，不可避免的需要同時(shí)對(duì)流場與溫度場進(jìn)行耦合計(jì)算[2]。國內(nèi)外學(xué)者針對(duì)此問題的緊耦合計(jì)算方法已經(jīng)展開了很多研究工作[3]。文獻(xiàn)[4]中Garimella等人完成了艙內(nèi)流場與傳熱一體化耦合計(jì)算，得到了較為理想的結(jié)果。文獻(xiàn)[5]中Hofacker等人利用PHOENICS軟件模擬了歐空局赫爾墨斯（Hermes）航天飛機(jī)艙內(nèi)對(duì)流換熱溫度場、流場和狹縫集中送風(fēng)的通風(fēng)效果。文獻(xiàn)[6]對(duì)航天器密封艙內(nèi)流動(dòng)和傳熱問題進(jìn)行了數(shù)值計(jì)算研究，得出了艙內(nèi)對(duì)流換熱系數(shù)與平均風(fēng)速的經(jīng)驗(yàn)關(guān)系式。文獻(xiàn)[7]利用CFD數(shù)值計(jì)算方法研究了載人空間站艙內(nèi)對(duì)流換熱的問題。但上述工作中均采取的是CFD緊耦合數(shù)值計(jì)算方法[8-10]，即計(jì)算時(shí)采用流場的特征時(shí)間為計(jì)算步長，特點(diǎn)是流場與結(jié)構(gòu)傳熱實(shí)時(shí)耦合，反映了物理實(shí)際，但計(jì)算需要耗費(fèi)大量時(shí)間，在實(shí)際工程中的應(yīng)用價(jià)值大打折扣，因此需要尋求一種快速便捷的仿真計(jì)算方法解決以上問題。

本文首先討論了數(shù)值計(jì)算中的耦合傳熱問題，分析了流固界面耦合求解原理，指出了不同求解方法的優(yōu)缺點(diǎn)。針對(duì)普遍使用的緊耦合方法計(jì)算時(shí)間過長的問題，結(jié)合一般松耦合求解思想提出了一種針對(duì)載人航天器艙內(nèi)流場/溫度場松耦合計(jì)算方法，并對(duì)其進(jìn)行了詳細(xì)的數(shù)值仿真驗(yàn)證。結(jié)果證明在流場隨時(shí)間變化不大的情況下，本文提出的流場/溫度場松耦合計(jì)算方法在求解瞬態(tài)耦合換熱問題時(shí)不但具有很好的魯棒性，而且在與緊耦合方法的計(jì)算結(jié)果進(jìn)行對(duì)比時(shí)，精度上也只存在較小的誤差。

1 耦合處理方法原理分析

緊耦合計(jì)算的特點(diǎn)是流場與固體計(jì)算域在接觸面實(shí)時(shí)交互數(shù)據(jù)、耦合迭代，主要體現(xiàn)在流體與固體接觸面上的邊界條件處理，即各自的溫度邊界條件[11-12]。流/固接觸面的耦合原理，如圖1所示。

圖1 流/固接觸面的耦合原理Fig.1 Coupling principle of fluid/solid interface

圖1 中，在直角坐標(biāo)系下，流場區(qū)域和固體結(jié)構(gòu)區(qū)域分別用FΩ和SΩ表示；Γ為流、固接觸面；sq為接觸面Γ傳至結(jié)構(gòu)內(nèi)部的熱流；xq和yq分別為沿x和y方向由流場傳入至結(jié)構(gòu)內(nèi)部的熱流。

熱傳導(dǎo)方程在邊界Γ上的定解條件是[10]：

式中 k為固體熱傳導(dǎo)系數(shù)；T為溫度；x為橫向坐標(biāo)；y為縱向坐標(biāo)；nx為沿x方向的單位矢量；ny為沿y方向的單位矢量。

從能量守恒出發(fā)，在忽略輻射效應(yīng)的假設(shè)前提下，流體區(qū)域傳至固壁的熱流應(yīng)等于固壁傳至結(jié)構(gòu)體內(nèi)部的熱流，因此有[10]：

由上式建立起流場與結(jié)構(gòu)傳熱計(jì)算的耦合關(guān)系，先給定時(shí)間為0、流場FΩ的初場和Γ上的溫度分布，通過求解Navier-Stokes方程[13]得到tΔ時(shí)間后壁面Γ上的熱流分布，此處tΔ為一個(gè)計(jì)算時(shí)間步長。然后以該熱流分布為邊界條件，用有限單元法計(jì)算tΔ時(shí)間后的結(jié)構(gòu)溫度分布，并得到壁面Γ上新的溫度分布，至此完成耦合計(jì)算的一個(gè)周期。如此反復(fù)迭代，直到計(jì)算要求的截止時(shí)刻。

而在松耦合方法中，假定從t=0時(shí)刻開始，Γ上的熱流分布為恒定值，即上述緊耦合計(jì)算的0時(shí)刻熱流分布，結(jié)構(gòu)區(qū)域SΩ從此作為邊界條件進(jìn)行熱傳導(dǎo)分析。利用松耦合方法求解流場與結(jié)構(gòu)傳熱耦合的非定常問題，可按“準(zhǔn)定?！狈椒ㄌ幚?。首先，把整個(gè)時(shí)間歷程劃分成一定的時(shí)間序列0t，1t，2t，…，1nt-，nt時(shí)間間隔tΔ的大小可視精度要求進(jìn)行控制。在每一時(shí)間間隔tΔ內(nèi)，認(rèn)為固體表面接收到一個(gè)固定的外表面熱流，具體步驟如下[14-15]：

1）給定初始t0時(shí)刻外表面壁溫和結(jié)構(gòu)內(nèi)溫度分布初場，算出穩(wěn)態(tài)流場，得到外表面溫度分布第一個(gè)時(shí)間步Δt1內(nèi)的外表面壁溫；

2）將外表面溫度分布作為邊界條件，求解得到t1時(shí)刻結(jié)構(gòu)溫度分布；

3）以t1時(shí)刻固體邊界溫度為邊界條件，求得t1時(shí)刻的穩(wěn)態(tài)流場；

4）凍結(jié)t1時(shí)刻的流場，繼續(xù)在t1～t2時(shí)間段內(nèi)單獨(dú)求解結(jié)構(gòu)瞬態(tài)熱傳導(dǎo)，獲得t2時(shí)刻的結(jié)構(gòu)溫度分布；

5）轉(zhuǎn)向步驟3），繼續(xù)計(jì)算直到整個(gè)時(shí)間推算完畢轉(zhuǎn)向執(zhí)行步驟2），直至整個(gè)時(shí)間推算完畢。

緊耦合方法中，雖然流場與結(jié)構(gòu)傳熱實(shí)時(shí)耦合，反映了物理實(shí)際，但其計(jì)算效率較為低下，不滿足實(shí)際應(yīng)用中的高效性要求。而針對(duì)此類耦合問題，如果計(jì)算對(duì)象的流場隨時(shí)間變化較為緩慢，即流場流動(dòng)特征時(shí)間與結(jié)構(gòu)傳熱特征時(shí)間相差較大，松耦合方法在結(jié)構(gòu)體每一時(shí)間步長的時(shí)間間隔內(nèi)可以認(rèn)為流場是穩(wěn)態(tài)，從而大大減少計(jì)算量。

2 流場/溫度場松耦合求解思想

在前面所介紹的松耦合求解方法中，在計(jì)算固體導(dǎo)熱的每一個(gè)時(shí)間步長內(nèi)都將流場視為穩(wěn)態(tài)，從而大大節(jié)省了整體的計(jì)算時(shí)間提高了效率。由此受到啟發(fā)，提出一種新的流場/溫度場松耦合求解方法，將其應(yīng)用到載人航天器艙內(nèi)對(duì)流換熱的CFD數(shù)值計(jì)算當(dāng)中，具體求解方法與步驟如下：

1）在一個(gè)時(shí)間步長1tΔ內(nèi)同時(shí)計(jì)算流場+溫度場，也就是說在1tΔ的時(shí)間步長內(nèi)對(duì)能量方程和動(dòng)量方程同時(shí)進(jìn)行迭代求解；

2）固定流場，在下一個(gè)時(shí)間步長2tΔ內(nèi)不再進(jìn)行計(jì)算流場變化，即不再對(duì)流體動(dòng)量方程進(jìn)行計(jì)算，只對(duì)能量方程進(jìn)行迭代求解一個(gè)步長2tΔ，得出溫度場；

3）固定溫度場不變，在1tΔ時(shí)間內(nèi)繼續(xù)同時(shí)求解流場；

4）重復(fù)步驟1）到步驟3），直到求解結(jié)束。

需要說明的是，一般的松耦合計(jì)算方法中是將1tΔ與2tΔ作為進(jìn)行迭代計(jì)算的步長，在2tΔ的計(jì)算步長內(nèi)將流場視為瞬態(tài)固定，而在本文中，1tΔ與2tΔ將作為進(jìn)行計(jì)算的時(shí)間步長，而在每一個(gè)時(shí)間步長內(nèi)會(huì)進(jìn)行若干次相應(yīng)的迭代計(jì)算，進(jìn)一步提高整體計(jì)算效率。

在流場/溫度場松耦合計(jì)算方法中，1tΔ與2tΔ時(shí)間步長的選取方式至關(guān)重要。根據(jù)非穩(wěn)態(tài)導(dǎo)熱問題數(shù)值解法原理[10]，時(shí)間計(jì)算步長越小，所需計(jì)算機(jī)內(nèi)存及計(jì)算時(shí)間就大大增加；時(shí)間步長選取過大，計(jì)算精度則相應(yīng)降低，綜合考慮以上因素，本文選取9種1tΔ與2tΔ的組合，如表1所示。

表1 不同時(shí)間步長組合的選取Tab.1 Different combination of time steps

3 計(jì)算模型的選取

本文選取的計(jì)算對(duì)象為某升力體載人航天器設(shè)備艙的艙段模型，如圖2所示。為方便進(jìn)行數(shù)值計(jì)算，提高計(jì)算效率，對(duì)模型進(jìn)行相應(yīng)簡化，簡化后的模型及艙內(nèi)7個(gè)設(shè)備的編號(hào)如圖3所示，模型尺寸單位為m。

圖2 計(jì)算模型示意Fig.2 Numerical simulation model

圖3 簡化后的計(jì)算模型及艙內(nèi)設(shè)備編號(hào)Fig.3 Simplified model and serial numbers of devices in cabin

數(shù)值計(jì)算的邊界條件設(shè)定如下：艙內(nèi)設(shè)備材料為鋁，艙蒙皮內(nèi)壁面（弧面）給定溫邊界條件100℃，其它壁面設(shè)為絕熱條件，CFD數(shù)值計(jì)算數(shù)學(xué)模型采取k-ε湍流計(jì)算模型[10]。

4 計(jì)算結(jié)果與分析

首先考察不同時(shí)間步長組合下基于本文松耦合求解思想的計(jì)算結(jié)果，與緊耦合計(jì)算結(jié)果對(duì)比，選取一組計(jì)算效率和計(jì)算精度均較優(yōu)的時(shí)間步長組合，然后對(duì)此時(shí)間步長組合下不同時(shí)刻的計(jì)算結(jié)果進(jìn)行對(duì)比分析。

4.1 不同時(shí)間步長下的計(jì)算結(jié)果

表2所示為采用本文松耦合方法在9種時(shí)間步長組合下進(jìn)行瞬態(tài)計(jì)算1000s所消耗的計(jì)算時(shí)間，并與緊耦合方法的計(jì)算結(jié)果進(jìn)行了對(duì)比。表中t為計(jì)算消耗時(shí)間，單位為h。

表2 不同時(shí)間步長組合下計(jì)算1000s所消耗的時(shí)間Tab.2 Time expended when calculating 1 000s at different time steps combination

從表2中可以看出，相比緊耦合計(jì)算方法，流場/溫度場松耦合計(jì)算方法在進(jìn)行瞬態(tài)計(jì)算時(shí)節(jié)省了大量的計(jì)算時(shí)間。結(jié)合表1中時(shí)間步長的選取可以看出，1tΔ越大，整體的計(jì)算時(shí)間就會(huì)越長，相反2tΔ的增大會(huì)減少計(jì)算時(shí)間，提高計(jì)算效率。也就是說，實(shí)際計(jì)算時(shí)間基本是與1tΔ/2tΔ的大小成正比關(guān)系。

很明顯，不同時(shí)間步長的組合對(duì)數(shù)值計(jì)算所用時(shí)間的影響很大，但在實(shí)際計(jì)算中，只一味地選擇計(jì)算耗時(shí)少的情況也是不可取的，因?yàn)橐惨瑫r(shí)考察不同計(jì)算方式下計(jì)算結(jié)果的精度問題。

取艙內(nèi)設(shè)備3～4的溫度分布結(jié)果作為分析對(duì)象，在1000s時(shí)，表1中所示的9種不同時(shí)間步長組合下的計(jì)算結(jié)果如圖4～6所示。圖中T為溫度；橫坐標(biāo)為沿艙段橫向的距離，參照?qǐng)D3所示。

圖4 計(jì)算時(shí)間1000s時(shí)算例1～3的溫度計(jì)算結(jié)果Fig.4 Temperature results of case 1～3 in 1 000s

圖5 計(jì)算時(shí)間1000s時(shí)算例4～6的溫度計(jì)算結(jié)果Fig.5 Temperature results of case 4～6 in 1 000s

圖6 計(jì)算時(shí)間1000s時(shí)算例7～9的溫度計(jì)算結(jié)果Fig.6 Temperature results of case 7～9 in 1 000s

從圖4中不難發(fā)現(xiàn)，當(dāng)流場與溫度場同時(shí)求解的時(shí)間步長1tΔ固定時(shí)，計(jì)算結(jié)果與緊耦合結(jié)果之間的誤差會(huì)隨著求解溫度場的時(shí)間步長2tΔ的增大而增大，當(dāng)1tΔ較大時(shí)，計(jì)算結(jié)果與緊耦合結(jié)果較為符合，這意味著相應(yīng)的計(jì)算時(shí)間也會(huì)較大，對(duì)計(jì)算效率沒有太多提高。

從圖5和圖6中可以看出，減小1tΔ的大小，同時(shí)選取合適的2tΔ的值，計(jì)算結(jié)果與緊耦合結(jié)果也會(huì)符合得很好。另一方面，在保證精度的同時(shí)也應(yīng)將計(jì)算時(shí)間盡量減少，結(jié)合表2可以發(fā)現(xiàn)在算例5的情況下計(jì)算1000s耗時(shí)為3h，比算例4情況下的4.3h和算例7情況下的3.5h要小，而三者的誤差均在1%～1.5%（見表3）。結(jié)合上述分析，最終選取算例5即1tΔ=10s和2tΔ=10s作為時(shí)間步長應(yīng)用到本文的流場/溫度場松耦合計(jì)算方法中，作進(jìn)一步計(jì)算分析。

上述分析結(jié)果與1000s時(shí)艙內(nèi)設(shè)備的平均溫度計(jì)算結(jié)果見表3，從中進(jìn)一步可以看出，與緊耦合計(jì)算方法相比，雖然算例1、算例4和算例7的計(jì)算誤差較算例5要小，但計(jì)算效率較差，因此算例5的時(shí)間步長組合選擇為最優(yōu)。設(shè)備表面的在流場中的平均對(duì)流換熱系數(shù)[13]見表4，從表4中發(fā)現(xiàn)不同時(shí)間步長組合對(duì)對(duì)流換熱系數(shù)的計(jì)算結(jié)果影響不大，當(dāng)算例9中引起的溫度誤差為7.51%時(shí)換熱系數(shù)的誤差卻僅為2.87%。

表3 瞬態(tài)1000s時(shí)刻，設(shè)備平均溫度在不同時(shí)間步長情況下的計(jì)算結(jié)果Tab.3 Tab. 3 Average temperature at t=1 000s at different time steps

表4 設(shè)備表面平均對(duì)流換熱系數(shù)的計(jì)算結(jié)果Tab.4 Average heat convection coefficient on the surface of devices

4.2 同時(shí)間步長、不同時(shí)刻計(jì)算結(jié)果

根據(jù)上一節(jié)中得出的結(jié)論，選取算例5中的Δt1=10s，Δt2=10s作為流場/松耦合計(jì)算方法中的時(shí)間步長。在瞬態(tài)計(jì)算時(shí)長為200、500、1000和1500s時(shí)，對(duì)比分析溫度場的（設(shè)備1到設(shè)備2之間）計(jì)算結(jié)果，如圖7～10所示。

圖7 t=200s時(shí)刻溫度場計(jì)算結(jié)果Fig.7 Temperature results at t = 200s

圖8 t=500s時(shí)刻溫度場計(jì)算結(jié)果Fig.8 Temperature results at t = 500s

圖9 t=1000s時(shí)刻溫度場計(jì)算結(jié)果Fig.9 Temperature results at t = 1 000s

圖10 t=1500s時(shí)刻溫度場計(jì)算結(jié)果Fig.10 Temperature results at t = 1 500s

可以看出本文流場/溫度場松耦合計(jì)算方法的結(jié)果與緊耦合方法符合得較好，但隨著計(jì)算時(shí)間的增大，二者計(jì)算結(jié)果的誤差也相應(yīng)增大，4個(gè)圖中設(shè)備的計(jì)算溫度出現(xiàn)最大誤差的時(shí)刻在1500s，但也僅為2%左右。艙內(nèi)設(shè)備在不同時(shí)刻下的溫度計(jì)算結(jié)果見表5所示。從表中可看出，在選取算例5進(jìn)行計(jì)算同時(shí)間步長下不同時(shí)刻的設(shè)備溫度時(shí)，與緊耦合方法相比，計(jì)算結(jié)果均能夠滿足精算精度要求，進(jìn)一步說明了本文流場/溫度場松耦合計(jì)算方法的合理性。

艙內(nèi)設(shè)備表面平均對(duì)流換熱系數(shù)的計(jì)算值見表6所示。與表4中的結(jié)果類似，艙內(nèi)設(shè)備表面的瞬態(tài)對(duì)流換熱系數(shù)基本不受計(jì)算方法和時(shí)間的影響。究其原因，一方面從對(duì)流換熱系數(shù)的定義分析，物體與流體接觸面的對(duì)流換熱系數(shù)h的大小主要取決于表面的幾何形狀、流體的運(yùn)動(dòng)特性及流體的熱力學(xué)與輸運(yùn)性質(zhì)，而受其它因素的影響較小[10,13]。

表5 不同時(shí)刻下設(shè)備平均溫度的計(jì)算結(jié)果Tab.5 Average temperature at different time

表6 不同時(shí)刻下設(shè)備表面平均對(duì)流換熱系數(shù)的計(jì)算結(jié)果Tab.6 Average heat convection coefficient at different time

5 結(jié)束語

文中首先討論了數(shù)值計(jì)算中的耦合對(duì)流換熱問題，分析了流固界面耦合傳熱求解原理，指出了不同求解方法的優(yōu)缺點(diǎn)。針對(duì)普遍使用的緊耦合方法計(jì)算時(shí)間過長的問題，結(jié)合一般松耦合求解思想提出了一種針對(duì)艙內(nèi)流場/溫度場的松耦合計(jì)算方法，并對(duì)其進(jìn)行了詳細(xì)的數(shù)值仿真驗(yàn)證。結(jié)果證明此種方法在求解瞬態(tài)耦合換熱問題時(shí)不但具有很好的魯棒性，與緊耦合方法的計(jì)算結(jié)果進(jìn)行對(duì)比，精度上也只存在較小的誤差。尤其是在計(jì)算對(duì)流換熱系數(shù)的時(shí)候發(fā)現(xiàn)，其計(jì)算結(jié)果不受到計(jì)算時(shí)間步長和求解時(shí)間的影響，這樣進(jìn)一步驗(yàn)證了本文提出的流場/溫度場松耦合計(jì)算方法的合理性。

References)

[1]于瑩瀟,田佳林.美國新型載人火星探測技術(shù)方案[J].航天返回與遙感,2009,30(3):1-7,15.YU Yingxiao, TIAN Jialin. USA

[2]黃家榮,范含林.載人航天器生活艙內(nèi)熱濕環(huán)境的數(shù)值模擬[J].中國空間科學(xué)技術(shù),2004(6):7-13.HUANG Jiarong, FAN Hanlin. Numerical Simulation of Thermal/Humidity in Manned Spacecraft Cabin[J]. Chinese Space Science and Technology, 2004(6): 7-13. (in Chinese)

[3]鄭忠海.空間站艙內(nèi)氣流分布特性數(shù)值模擬與熱舒適性評(píng)價(jià)[D].哈爾濱:哈爾濱工業(yè)大學(xué),2005.ZHENG Zhonghai. Numerical Simulation of Air Distribution in Space Station and Thermal Evaluation[D]. Harbin: Harbin Institute of Technology, 2005. (in Chinese)

[4]Garimella S. Flow and Heat Transfer in Space Vehicle Tile Gaps[R]. AIAA Paper No. 92-4061, Houston: American Institute of Aeronautics and Astronautics, 1992.

[5]Hofacker W, Lucas J. Thermal Analysis of Spacecraft by Combining Different Analysis Tools[R]. SAE911584, Detroit:Society of Automotive Engineers, 2009.

[6]鐘奇,劉強(qiáng).航天器密封艙流動(dòng)和傳熱的數(shù)值研究[J].宇航學(xué)報(bào),2002,23(5):44-48.ZHONG Qi, LIU Qiang. A Numerical Investigation on Heat Transfer and Flow in a Pressured Cabin of Spacecraft[J]. Journal of Astronautics, 2002, 23(5): 44-48. (in Chinese)

[7]張吉禮,梁珍,鄭忠海,等.載人航天空間站艙內(nèi)通風(fēng)對(duì)流換熱數(shù)值研究進(jìn)展[J].暖通空調(diào),2006,36(1):28-34.ZHANG Jili, LIANG Zhen, ZHENG Zhonghai, etal. Numerical Simulation Progresses of Convection Heat Transfer in Manned Spacecraft Cabin[J]. Heating, Ventilating and Air Conditioning, 2006, 36(1): 28-34. (in Chinese)

[8]Patankar S V. Numerical Heat Transfer and Fluid Flow[M]. New York: McGraw-Hill, 1980: 125-129.

[9]張政.傳熱與流體流動(dòng)的數(shù)值計(jì)算[M].北京:科學(xué)出版社,1984:25-28.ZHANG Zheng. Numerical Simulation of Heat Transfer and Fluid Flow[M]. Beijing: Science Press, 1984: 25-28. (in Chinese)

[10]陶文銓.數(shù)值傳熱學(xué)[M].第2版.西安:西安交通大學(xué)出版社,2001:112-174.TAO Wenquan. Numerical Heat Transfer[M]. Second Edition. Xi’an: Xi’an Jiaotong University Press, 2001: 112-174. (inChinese)

[11]Zabrodskiy S S. Several Specific Terms and Concepts in Russian Heat Transfer Notation[J]. Heat Transfer-soviet Research,1972, 4(2): 1-4.

[12]Versteeg H K, Malasekera W. An Introduction to Computational Fluid Dynamics: The Fine Volume Method[M]. Essex:Longman Scientific Technical, 1995: 35-68.

[13]Bejan A. Convection Heat Transfer [M]. New York: Wiley, 2004: 334-384.

[14]Sun Y S, Emery A F. Multigrid Computation of Natural Convection in Enclosures with a Conductive Baffle[J]. Numerical Heat Transfer, Part A, 1994(25): 575-592.

[15]Anand N K, Chin C D, Mcmath J G. Heat Transfer in Rectangular Channels with a Series of Normally in-line Positioned Plates [J]. Numerical Heat Transfer, Part A, 1995(27): 19-34.