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

（北京理工大學(xué) 宇航學(xué)院，北京100081）

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

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

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

1 翼面氣動加熱熱流

運(yùn)用牛頓公式，氣動加熱熱流為

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

1.1 飛行器翼面駐點(diǎn)熱流

采用Kemp－Riddell公式計(jì)算翼面前緣駐點(diǎn)熱流：

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

1.2 飛行器翼面非駐點(diǎn)熱流

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

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

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

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

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

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

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

2 翼面輻射熱流

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

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

3 翼面熱傳導(dǎo)

熱傳導(dǎo)微分方程為

熱流邊界條件為

式中：k為導(dǎo)熱系數(shù)，ρ為物體的密度，c為比熱容，W為內(nèi)部熱源強(qiáng)度。

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

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

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

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

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

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

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

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

圖2 翼面?zhèn)鳠狁詈嫌?jì)算流程圖

5 算例分析

5.1 氣動加熱工程算法驗(yàn)證

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

圖3 CFD二維翼面網(wǎng)格圖

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

圖4 熱流數(shù)值計(jì)算與工程計(jì)算對比

5.2 氣動加熱與熱傳導(dǎo)耦合分析

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

圖5 高超聲速飛行器三維翼面結(jié)構(gòu)圖

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

5.2.1 氣動加熱、輻射換熱、熱傳導(dǎo)耦合分析與非耦合分析對比

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

圖6 非耦合各部分輸入輸出關(guān)系

圖7 耦合100時(shí)間步長翼面溫度分布

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

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

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

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

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

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

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

6 結(jié)束語

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

［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］.南京航空航天大學(xué)學(xué)報(bào)，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］黃唐，毛國良，姜貴慶，等.二維流場、熱、結(jié)構(gòu)一體化數(shù)值模擬［J］.空氣動力學(xué)學(xué)報(bào)，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］楊榮，王強(qiáng).高超聲速旋轉(zhuǎn)體氣動加熱、輻射換熱與結(jié)構(gòu)熱傳導(dǎo)的耦合數(shù)值分析［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］張志成.高超聲速氣動熱和熱防護(hù)［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ì)算方法［J］.彈箭與制導(dǎo)學(xué)報(bào)，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.