陳 鑫，劉 莉，李昱霖，楊 武

（北京理工大學 宇航學院，北京100081）

隨著飛行馬赫數的提高，氣動加熱效應對高超聲速飛行器的影響愈來愈不可忽略。準確地預測高超聲速飛行器的熱環(huán)境是飛行器設計及飛行試驗成敗的關鍵難點之一。根據國外在NASA研究開發(fā)計劃中的經驗，利用分層求解［1］的思路可以滿足工程的需求，即進行氣動加熱計算時，只考慮與飛行器表面輻射換熱達到熱平衡，暫不考慮表面向結構內部的傳導。

國外關注氣動熱、熱傳導及結構的耦合問題的工作開展較早，1989年文獻［2］完成了流場、傳熱與結構的一體化耦合計算；文獻［3］從數值模擬角度對氣動加熱與瞬態(tài)熱傳導耦合性的影響進行了分析，給出了部分模擬計算結果，指出了耦合性在分析持續(xù)氣動加熱問題中的必要性；文獻［4］給出了一種把流場、熱、結構耦合起來進行一體化數值模擬的方法，以簡單的二維圓柱繞流為例，驗證了該方法的可行性。文獻［5］采用流場－熱－結構耦合的計算方法，耦合氣動加熱，輻射換熱及結構熱傳導，數值模擬了旋轉體的外流場和壁面結構溫度場的非定常過程。

上述參考文獻指出了考慮氣動加熱、輻射換熱與熱傳導耦合的必要性，但均采用了數值模擬氣動加熱和結構熱傳導過程，計算效率較低，計算模型多為二維簡單模型。利用氣動加熱工程算法能極大地減少求解流場的計算時間，比傳統(tǒng)的數值計算效率更高。本文利用以經驗公式為基礎的工程算法計算氣動加熱熱流，考慮輻射換熱熱流，利用有限單元法計算瞬態(tài)熱傳導過程，并通過實時更新邊界條件，實現(xiàn)了三維翼面的氣動加熱、輻射換熱與結構熱傳導的耦合求解。

1 翼面氣動加熱熱流

運用牛頓公式，氣動加熱熱流為

式中：αh為焓值熱交換系數，hr為恢復焓值，hw為飛行器壁面焓值。將飛行器升力面分為前緣與非前緣兩部分分別進行氣動加熱熱流計算。

1.1 飛行器翼面駐點熱流

采用Kemp－Riddell公式計算翼面前緣駐點熱流：

式中：vc＝7 900m／s，ρsl＝1.225kg／m3為大氣海平面密度，RN為駐點曲率半徑，ρ∞、v∞分別為無窮遠來流密度及速度，h0為駐點總焓，h300為300K時空氣的焓值。

1.2 飛行器翼面非駐點熱流

采用片條理論，翼面弦向氣動加熱熱流［6］為

式中：ρ＊、μ＊分別為參考焓下氣體密度及粘性系數，ρe、μe分別為邊界層外緣氣體密度及粘性系數，Rex為參考焓下的氣體雷諾數，帶上標“＊”的參數為參考焓下的氣流參數。

參考焓法的基本思想是利用不可壓縮邊界層理論求得的公式來估算可壓縮邊界層中的摩擦和傳熱。Eckert參考焓公式［7］最常用，在層流及湍流中均與試驗吻合得較好。即：

式中：h＊為參考焓值，he為邊界層外緣焓值。

1.3 邊界層外緣及翼面壁面氣流參數的確定

利用布澤曼理論可以考慮不同飛行參數及翼型參數對翼面表面溫度分布的影響［8］。利用布澤曼理論得到迎風面邊界層外緣氣流參數：

式中：帶下標e的參數為邊界層外緣參數，帶下標c的參數為翼面前緣參數，Cp為翼面表面壓力系數，Cp0為翼面前緣壓力系數，φ為迎風面某點處翼型切線與來流方向的夾角，λ為比熱比，pe、p∞分別為翼面前緣和來流壓強，Ma∞為來流馬赫數，ρe、ρ∞為翼面前緣和來流密度，Te、T∞分別為翼面前緣和來流溫度。

2 翼面輻射熱流

飛行器壁面向外輻射的熱流量為

式中：波爾茲曼常數σ＝57.6×10－9W／（m2·K4），ε為蒙皮表面輻射系數，Tw為壁面溫度。

3 翼面熱傳導

熱傳導微分方程為

熱流邊界條件為

式中：k為導熱系數，ρ為物體的密度，c為比熱容，W為內部熱源強度。

4 翼面氣動加熱、輻射換熱與瞬態(tài)熱傳導耦合分析

高超聲速飛行器飛行過程中，翼面壁面氣動加熱熱流qw、輻射熱流qr和壁面溫度Tw相互耦合，通過邊界條件相互作用，如圖1所示。氣動加熱熱流qw求解需要翼面壁面溫度Tw，Tw求解需要以熱流qw為邊界條件；輻射換熱熱流qr求解需要翼面溫度Tw，Tw求解需要以熱流qr為邊界條件求解。數學上表現(xiàn)為各部分的求解對象是其它求解對象的定解條件，即三部分求解對象之間是互為定解條件的。

圖1 翼面?zhèn)鳠岣鞑糠州斎胼敵鲫P系

傳統(tǒng)處理氣動加熱、輻射換熱、熱傳導需要聯(lián)立三部分方程耦合求解，這樣處理求解復雜，計算繁瑣，實際操作困難。本文嘗試采用一種準定常求解方法，其計算思路如下：

①給定高超聲速飛行器典型彈道下隨時間變化的飛行條件、初始時刻的翼面溫度；

②由初始時刻的飛行條件、翼面溫度求解翼面氣動加熱與輻射換熱熱流，并以此為熱流邊界條件求解結構熱傳導方程，得到單位時間步長的新的翼面溫度分布；

③若新的溫度分布與原有翼面溫度分布滿足收斂條件，則更新到下一步時間步長，回到①繼續(xù)進行下一時間步迭代，直到整個彈道計算完畢，否則更新翼面溫度返回②直到計算收斂。

給定典型彈道，考慮傳熱耦合性，氣動加熱、輻射換熱與熱傳導耦合算法計算流程圖如圖2所示。

圖2 翼面?zhèn)鳠狁詈嫌嬎懔鞒虉D

5 算例分析

5.1 氣動加熱工程算法驗證

對某一高超聲速二維翼面進行氣動加熱數值及工程算法計算，飛行高度H＝30km，馬赫數Ma＝5.0，來流溫度T＝226.51K，攻角為0，翼面弦長為0.5m。CFD－FASTRAN是目前廣泛應用于高超聲速流的專用CFD商業(yè)軟件［9］。本文高超聲速氣動熱計算采用CFD－FASTRAN求解器計算熱流分布。數值計算邊界條件為壓力遠場、壓力出口和壁面邊界。翼面邊界表面考慮輻射，輻射系數為0.85，時間離散采用隱式Jacobi點迭代方法進行推進，空間離散采用一階Roe’s FDS格式。圖3為利用數值計算的CFD－FASTRAN計算結果網格圖。

圖3 CFD二維翼面網格圖

圖4 為熱流數值計算與工程計算的對比圖，由圖4可知，氣動加熱工程算法與數值計算熱流在駐點及弦向上均較為接近，在駐點熱流求解上，工程算法與數值求解結果吻合較好，約為60.7kW／m2。在x＝0.25m處，數值求解考慮翼型時出現(xiàn)熱流較小的跳變，在弦向0.25m以后熱流均略微小于工程算法熱流。

圖4 熱流數值計算與工程計算對比

5.2 氣動加熱與熱傳導耦合分析

某一高超聲速三維翼面如圖5所示，翼面根弦長為0.5m，根梢比為0.45，展長為0.287 5m翼面表面蒙皮采用鈦合金，厚度2.5mm。

圖5 高超聲速飛行器三維翼面結構圖

飛行器助推飛行時間t＝0～50s，馬赫數Ma＝1.5～5.5，飛行助推段H＝15～30km。

5.2.1 氣動加熱、輻射換熱、熱傳導耦合分析與非耦合分析對比

不考慮氣動加熱與熱傳導耦合的輸入輸出關系，如圖6所示。以0.001s為時間步長，在典型彈道上進行100步長時間后，考慮氣動加熱、輻射換熱與熱傳導耦合的方法得到翼面溫度，如圖7所示，圖中x為翼面弦長坐標，y為翼面展長坐標。非耦合翼面溫度分布如圖8所示，耦合與非耦合前緣最大相對溫度差約為40%，且耦合方法求得的翼面分布更均勻，更能反映翼面真實的溫度分布。非耦合求解時氣動加熱熱流計算不能利用實時變化的翼面溫度，隨著飛行時間增加，氣動加熱熱流將明顯增大，不能反映真實的翼面熱流。隨著飛行時間的增加，非耦合求解的氣動加熱熱流在翼面的積累效應將越來越明顯，非耦合翼面分布溫度越來越高于耦合翼面溫度。故預測翼面熱環(huán)境時必須考慮氣動加熱與瞬態(tài)熱傳導的耦合效應。

圖6 非耦合各部分輸入輸出關系

圖7 耦合100時間步長翼面溫度分布

圖8 非耦合100時間步長翼面溫度分布

5.2.2 氣動加熱、輻射換熱、熱傳導耦合分析與分層求解方法對比

分層求解思想指求解翼面溫度時僅考慮氣動加熱熱流與輻射換熱熱流瞬時平衡，即qw＝qr，從而得到表面輻射熱平衡溫度Tw。由美國NASA研究經驗［1］可以看出，分層求解方法可以滿足一定的工程需求。

針對上述高超聲速翼面，以0.001s為時間步長進行耦合算法迭代，為了與分層求解對比研究，耦合求解初始溫度為零時刻分層求解翼面的穩(wěn)態(tài)溫度，整個典型彈道求解需時約為55h。

耦合與分層求解根部弦向溫度對比如圖9所示，耦合與分層求解根部前緣溫度－時間曲線如圖10所示。如圖9、圖10所示，考慮耦合時整個翼面弦向溫度分布與分層求解溫度分布相當，前緣溫度分布和弦向溫度分布相差不大，但考慮耦合求解的溫度分布低于分層求解的溫度分布，對比發(fā)現(xiàn)耦合求解翼面弦向溫度比分層求解溫度降低約5%。同時，隨著飛行時間的增加，耦合求解的溫度越來越低于分層求解溫度。通過與分層求解結果對比，可以看出，本文的耦合求解方法能夠反映翼面的溫度分布，是一種可行的計算氣動加熱、輻射換熱、熱傳導耦合分析的方法，同時比分層求解所得的溫度分布更低?？紤]氣動加熱、輻射換熱與瞬態(tài)熱傳導耦合后，僅考慮氣動加熱熱流與輻射換熱流平衡的分層求解方法求得的翼面溫度將高于耦合求解翼面的溫度。

圖9 耦合與分層求解根部弦向溫度對比

圖10 耦合與分層求解根部前緣溫度－時間圖

6 結束語

本文給出了一種求解氣動加熱、輻射換熱與瞬態(tài)熱傳導耦合求解的方法，首先保證了氣動加熱的有效性，其次實現(xiàn)了高超聲速飛行器在典型彈道飛行條件下三維機翼的氣動加熱、輻射換熱、瞬態(tài)熱傳導的耦合求解，與非耦合的求解方法對比，指出了考慮氣動加熱、輻射換熱與瞬態(tài)熱傳導耦合的必要性。與工程上的分層求解方法相比，該耦合求解方法能獲得相近的翼面溫度分布，且耦合求解的翼面溫度更低。

［1］SPAIN C，SOISTMANN D，PARK E.An overview of selected NASP aeroelastic studies at the NASA Langley Research Center，AIAA－90－5218［R］.1990.

［2］DECHAUMPHAI P，THORNTON E A，WIETING A R.Flow－thermal－structural study of aerodynamically heated leading edges［J］.J Spacecraft，1989，26（4）：201－209.

［3］程克明，呂英偉.飛行器持續(xù)氣動加熱的耦合性分析［J］.南京航空航天大學學報，2000，32（3）：150－155.CHENG Ke－ming，LV Ying－wei.An analysis of coupling feature in continuing gasdynamic heating over flight vehicle［J］.Journal of Nanjing University of Aeronautics and Astronautics，2000，32（2）：150－155.（in Chinese）

［4］黃唐，毛國良，姜貴慶，等.二維流場、熱、結構一體化數值模擬［J］.空氣動力學學報，2000，18（1）：115－119.HUANG Tang，MAO Guo－liang，JIANG Gui－qing，et al.Two dimensional coupled flow－thermal－structural numerical simulation［J］.Acta Aerodynamica Sinica，2000，18（1）：115－119.（in Chinese）

［5］楊榮，王強.高超聲速旋轉體氣動加熱、輻射換熱與結構熱傳導的耦合數值分析［J］.上海航天，2009（4）：25－29.YANG Rong，WANG Qiang.Coupled numerical study on aroheating，radiative heat transfer and structure heat conduction for hypersonic bodies of revolution［J］.Aerospace Shanghai，2009（4）：25－29.（in Chinese）

［6］張志成.高超聲速氣動熱和熱防護［M］.北京：國防工業(yè)出版社，2003.ZHANG Zhi－cheng.Hypersonic heating and thermal protection［M ］.Beijing：National Defense Industry Press，2003.（in Chinese）

［7］ECKERT R G.Engineering relations for friction and heat transfer to surfaces in high velocity flow［J］.AIAA Journal，1955，22（8）：585－587.

［8］陳鑫，劉莉，李昱霖，等.高超聲速飛行器翼面氣動加熱的工程計算方法［J］.彈箭與制導學報，2013，33（3）：133－137.CHEN Xin，LIU Li，LI Yu－lin，et al.Engineering calculation of aerodynamic heating for airfoils of hypersonic vehicles［J］.Journal of Projectiles，Rockets，Missiles and Guidance，2013，33（3）：133－137.（in Chinese）

［9］LEI Jian，ZHAO Xiang，ZHANG S J.Hypersonic non－equilibrium computations for ionizing air，AIAA 2009－1591［R］.2009.