李志杰 果琳麗

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

0 引言

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

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

1 耦合處理方法原理分析

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

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

圖1 中，在直角坐標系下，流場區(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為橫向坐標；y為縱向坐標；nx為沿x方向的單位矢量；ny為沿y方向的單位矢量。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 計算模型的選取

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

圖2 計算模型示意Fig.2 Numerical simulation model

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)束語

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

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ù)值模擬與熱舒適性評價[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]鐘奇,劉強.航天器密封艙流動和傳熱的數(shù)值研究[J].宇航學(xué)報,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)對流換熱數(shù)值研究進展[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]張政.傳熱與流體流動的數(shù)值計算[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.