李佳偉, 王江峰, 程克明, 伍貽兆

(南京航空航天大學(xué) 航空學(xué)院, 江蘇 南京 210016)

0 引 言

高超聲速飛行器在大氣層中持續(xù)飛行時(shí)，飛行器表面將經(jīng)受嚴(yán)峻的氣動(dòng)加熱[1]。高超聲速氣動(dòng)熱會(huì)導(dǎo)致飛行器結(jié)構(gòu)溫度急劇升高，嚴(yán)重時(shí)局部高溫還可能導(dǎo)致飛行器結(jié)構(gòu)材料的物理屬性和剛度發(fā)生改變，對(duì)飛行器的飛行安全造成極大隱患。因此快速預(yù)測(cè)氣動(dòng)加熱與結(jié)構(gòu)傳熱的物理過程，對(duì)高超聲速飛行器設(shè)計(jì)的初步選型和熱防護(hù)系統(tǒng)輕量化設(shè)計(jì)，具有重要作用[2]。

高超聲速氣動(dòng)加熱問題[3]的主要研究方法包括：數(shù)值模擬方法[4]、工程計(jì)算方法[5-6]、地面風(fēng)洞試驗(yàn)[7]及飛行試驗(yàn)等，其中后兩種方法因代價(jià)高昂等因素不適于工程設(shè)計(jì)的初期選型階段。

數(shù)值模擬方法主要是對(duì)N-S(Navier-Stokes)方程及相關(guān)簡化形式的控制方程進(jìn)行求解，具有計(jì)算精度高、可處理復(fù)雜流動(dòng)和全機(jī)外形等優(yōu)點(diǎn)，但在氣動(dòng)熱與結(jié)構(gòu)傳熱耦合問題的求解方面對(duì)計(jì)算資源需求巨大且非常耗時(shí)。在這方面的研究，相關(guān)學(xué)者做了大量工作，主要工作集中在計(jì)算效率與精度等方面。董維中[8]等開展了高超聲速飛行器氣動(dòng)加熱與飛行器表面溫度耦合的數(shù)值模擬研究，結(jié)果表明在氣動(dòng)熱環(huán)境的預(yù)測(cè)中，應(yīng)采用表面溫度分布與氣動(dòng)熱耦合計(jì)算的方法，以減小表面溫度分布對(duì)氣動(dòng)熱計(jì)算結(jié)果的影響。楊愷[9]等采用基于點(diǎn)隱式的二階精度迎風(fēng)TVD格式的N-S方程有限體積法，開展了高溫?zé)峄瘜W(xué)非平衡條件下球錐和RAM-CII模型的氣動(dòng)熱數(shù)值模擬計(jì)算，分析了不同非平衡模型在熱化學(xué)非平衡條件下流場(chǎng)的影響。計(jì)算結(jié)果表明建立的數(shù)值模擬方法具有較高的精度。夏剛[10]等依據(jù)流場(chǎng)特征時(shí)間步長遠(yuǎn)小于結(jié)構(gòu)傳熱時(shí)間步長的特點(diǎn)，提出了一種松耦合的流-固-熱耦合計(jì)算方法，并采用二維高超聲速圓管繞流與結(jié)構(gòu)傳熱的非定常算例對(duì)方法進(jìn)行了驗(yàn)證。耿湘人[11]等基于Levelset方法和有限差分方法，將流場(chǎng)與固體結(jié)構(gòu)統(tǒng)一到同一組控制方程中進(jìn)行求解，計(jì)算結(jié)果與試驗(yàn)值吻合良好，探索了流-固-熱一體化計(jì)算的可行性。張勝濤[12]等針對(duì)高超聲速流場(chǎng)-熱-結(jié)構(gòu)松耦合分析策略框架，采用自適應(yīng)耦合計(jì)算時(shí)間步長、混合插值策略等方法，建立了多場(chǎng)耦合分析平臺(tái)。由于流-固-熱物理問題復(fù)雜，該問題的求解在數(shù)值計(jì)算技術(shù)及計(jì)算設(shè)備硬件條件等方面要求甚高，采用全流場(chǎng)黏性數(shù)值模擬方法難以快速得到工程設(shè)計(jì)初期選型所需的參考數(shù)據(jù)。

氣動(dòng)加熱工程計(jì)算方法[13]對(duì)簡單外形的氣動(dòng)熱求解具有高效、準(zhǔn)確的特點(diǎn)，因此在實(shí)際工程應(yīng)用中率先得到發(fā)展。但在復(fù)雜氣動(dòng)外形的氣動(dòng)熱分析方面適應(yīng)性較差，并且需要基于大量試驗(yàn)數(shù)據(jù)對(duì)計(jì)算方法與結(jié)果進(jìn)行人工修正。在這方面比較著名的是NASA蘭利研究中心研發(fā)的一套氣動(dòng)加熱預(yù)測(cè)程序(MINIVER)[14]，該程序中駐點(diǎn)區(qū)域使用經(jīng)典的Fay-Riddle公式，采用參考焓方法來計(jì)算高速流動(dòng)壓縮性效應(yīng)，另外可對(duì)過渡流區(qū)氣動(dòng)熱進(jìn)行計(jì)算。Hamilton[15]針對(duì)三維鈍頭體發(fā)展了一種適用于空氣平衡氣體的高超聲速流動(dòng)熱流密度計(jì)算方法，該方法可計(jì)算不同高速流動(dòng)狀態(tài)(層流、轉(zhuǎn)捩、湍流)下的熱流密度。Zoby[16]等人開發(fā)了LATCH方法，該方法基于參考焓和修正雷諾比擬計(jì)算熱流密度，可用于有化學(xué)反應(yīng)參與的鈍頭體氣動(dòng)熱計(jì)算。樂嘉陵[17]針對(duì)高超聲速流動(dòng)的黏性效應(yīng)，利用邊界層近似方法和工程方法對(duì)高超聲速飛行器的傳熱系數(shù)和表面摩擦阻力進(jìn)行了相應(yīng)計(jì)算。李建林[18]等人對(duì)氣動(dòng)熱工程計(jì)算方法進(jìn)行拓展應(yīng)用，對(duì)乘波體構(gòu)型高速飛行器進(jìn)行氣動(dòng)熱快速估算，得到較好結(jié)果?？傮w而言，工程計(jì)算方法的缺陷在于需要大量人工干涉以及龐大試驗(yàn)數(shù)據(jù)，在處理復(fù)雜外形飛行器方面的通用性上仍需完善。

綜上，本文依據(jù)高超聲速邊界層[19]理論，利用數(shù)值計(jì)算方法與工程算法的各自優(yōu)勢(shì)，將氣動(dòng)熱的計(jì)算簡化為繞飛行器的無黏外流(邊界層以外)數(shù)值解和邊界層內(nèi)熱流求解兩個(gè)部分，同時(shí)耦合防熱系統(tǒng)結(jié)構(gòu)傳熱計(jì)算模型，考慮了高溫化學(xué)非平衡效應(yīng)，發(fā)展了一種可用于快速計(jì)算三維復(fù)雜外形飛行器氣動(dòng)加熱與結(jié)構(gòu)傳熱特性的方法，實(shí)現(xiàn)了復(fù)雜飛行條件下飛行器全機(jī)表面熱流密度與防熱結(jié)構(gòu)溫度場(chǎng)時(shí)變特性的快速計(jì)算。

該方法的優(yōu)點(diǎn)在于：綜合了全流場(chǎng)數(shù)值模擬與工程算法各自的優(yōu)勢(shì)，可以快速、高效地對(duì)三維復(fù)雜外形高速飛行器的多種飛行狀態(tài)下的氣動(dòng)加熱與結(jié)構(gòu)耦合傳熱特性進(jìn)行計(jì)算與分析，給出熱流密度與防熱結(jié)構(gòu)層內(nèi)溫度等參數(shù)分布，彌補(bǔ)了全流場(chǎng)黏性數(shù)值模擬方法計(jì)算代價(jià)高、效率低、周期長等缺陷，同時(shí)拓展了氣動(dòng)熱工程算法的應(yīng)用范圍。

1 無黏外流場(chǎng)數(shù)值計(jì)算方法

氣動(dòng)加熱與結(jié)構(gòu)傳熱耦合數(shù)值計(jì)算技術(shù)考慮了瞬態(tài)氣動(dòng)加熱、巡航狀態(tài)長時(shí)傳熱和彈道狀態(tài)長時(shí)傳熱三種模型。根據(jù)不同飛行狀態(tài)完成無黏流場(chǎng)的求解，取無黏流場(chǎng)解結(jié)果中的物面參數(shù)(如表1所示)作為邊界層外緣參數(shù)提供給工程方法中的邊界層簡化算法作為輸入條件，之后可得到用于防熱系統(tǒng)中結(jié)構(gòu)傳熱計(jì)算所需的飛行器表面熱流及傳熱系數(shù)等參數(shù)。

表1 無黏外流解物面輸出參數(shù)

在無黏外流場(chǎng)求解方面，本文采用基于塊結(jié)構(gòu)網(wǎng)格技術(shù)及分布式并行計(jì)算技術(shù)的三維流場(chǎng)數(shù)值計(jì)算方法，為氣動(dòng)熱工程估算方法提供如表1所列的物面流動(dòng)參數(shù)。

2 物面熱流與結(jié)構(gòu)傳熱工程估算方法

2.1 物面熱流密度計(jì)算

物面熱流的計(jì)算采用工程中簡化求解邊界層方程的方法。將該工程方法所需的邊界層外緣參數(shù)設(shè)定為無黏流場(chǎng)解的物面流動(dòng)參數(shù)(見表1)，計(jì)算輸出結(jié)果為物面熱流密度。

本文研究針對(duì)的是全機(jī)外形，因此根據(jù)熱流密度工程計(jì)算方法將飛行器表面熱流的計(jì)算劃分為駐點(diǎn)區(qū)和非駐點(diǎn)區(qū)兩個(gè)區(qū)域。駐點(diǎn)區(qū)熱流計(jì)算采用目前廣泛使用的Fay-Riddell公式[20]：

(1)

式中ρw、μw、hw分別表示物面上氣體的密度、黏性系數(shù)和焓，ρs和μs分別表示駐點(diǎn)處的氣體密度及氣體黏性系數(shù)，hD為平均空氣電離焓。計(jì)算中取普朗特?cái)?shù)Pr=0.71，路易斯數(shù)Le=1.0。非駐點(diǎn)區(qū)的熱流密度計(jì)算公式參見文獻(xiàn)[21]。

2.2 結(jié)構(gòu)傳熱耦合計(jì)算方法

采用工程中常用的平板傳熱模型進(jìn)行結(jié)構(gòu)傳熱計(jì)算。傳熱方程為一維熱傳導(dǎo)方程，熱防護(hù)結(jié)構(gòu)內(nèi)表面采用絕熱壁邊界條件，采用差分法進(jìn)行求解。在傳熱計(jì)算模型的選取上，根據(jù)式(2)定義的防熱結(jié)構(gòu)材料的畢奧數(shù)Bi的取值來選取[22]。

(2)

式中α為傳熱系數(shù)，δ為材料結(jié)構(gòu)厚度，λδ為當(dāng)前材料熱傳導(dǎo)系數(shù)。

當(dāng)Bi≤0.1時(shí)，采用式(3)所示的熱薄壁傳熱模型：

(3)

初始條件為：

Tw|t=0=T0

(4)

式中ρδ、cδ和ε分別表示防熱材料的密度、比熱容和表面輻射系數(shù)。采用差分方法對(duì)式(3)進(jìn)行時(shí)間t的推進(jìn)求解，即可得到熱薄壁條件下熱防護(hù)結(jié)構(gòu)層溫度隨時(shí)間的變化。

當(dāng)Bi>0.1時(shí)，采用熱厚壁傳熱模型。在具體計(jì)算時(shí)，將熱厚壁由內(nèi)向外分為j層進(jìn)行逐層推進(jìn)求解。計(jì)算方程為：

表層：

(5)

最內(nèi)層：

(6)

中間層：

(7)

初始條件為：

Tj|t=0=T1|t=0=Tn|t=0=T0

(8)

式中下標(biāo)m為材料類型，n表示結(jié)構(gòu)材料的分層數(shù)，λm為材料m的熱傳導(dǎo)系數(shù)，該計(jì)算方法可用于計(jì)算多種材料組成的熱防護(hù)結(jié)構(gòu)的傳熱情況。將式(5)～式(7)進(jìn)行差分離散并按時(shí)間步長推進(jìn)求解，可以計(jì)算得出各層材料溫度隨時(shí)間的變化。氣動(dòng)加熱與結(jié)構(gòu)傳熱過程的耦合計(jì)算[23-27]比較復(fù)雜，在求解氣動(dòng)熱參數(shù)時(shí)，需要以物面溫度Tw作為物面邊界輸入條件，然后通過簡化邊界層工程方法計(jì)算得到表面熱流密度qn；而在計(jì)算結(jié)構(gòu)傳熱時(shí)，又以表面熱流密度qn作為熱邊界輸入條件，計(jì)算得到壁面溫度。因此，本文采用圖1所示的耦合迭代計(jì)算思路[23]。

圖1 氣動(dòng)加熱與結(jié)構(gòu)傳熱耦合求解示意圖

3 彈道狀態(tài)動(dòng)態(tài)插值方法

文中所述的彈道狀態(tài)是指飛行器在真實(shí)飛行中隨飛行時(shí)間改變飛行高度、迎角、馬赫數(shù)三個(gè)參數(shù)的狀態(tài)(暫不考慮側(cè)滑角)。

根據(jù)第2小節(jié)中所述的物面熱流與結(jié)構(gòu)傳熱工程計(jì)算方法的特點(diǎn)，將其擴(kuò)展到隨時(shí)間變化非定常彈道(飛行包線)狀態(tài)下的氣動(dòng)加熱與結(jié)構(gòu)傳熱耦合計(jì)算，此時(shí)需要足夠多的無黏外流的數(shù)值解作為輸入?yún)?shù)。為了在滿足工程設(shè)計(jì)誤差要求前提下提高耦合算法的計(jì)算效率，提出了一種無黏外流解動(dòng)態(tài)插值方法：即通過已經(jīng)預(yù)先完成的有限個(gè)無黏外流解的流場(chǎng)結(jié)果，插值得到當(dāng)前彈道時(shí)間點(diǎn)上飛行條件下的流場(chǎng)解，以供氣動(dòng)加熱與結(jié)構(gòu)傳熱耦合計(jì)算使用。

由于高超聲速飛行器飛行高度跨度大，通常涵蓋整個(gè)大氣層高度，若針對(duì)飛行高度H進(jìn)行插值，必須要先構(gòu)建飛行器在各個(gè)飛行高度條件下的流場(chǎng)數(shù)據(jù)庫，其計(jì)算量非常龐大，難以達(dá)到快速計(jì)算的目的。針對(duì)該問題，參考國內(nèi)外有關(guān)大氣參數(shù)隨大氣高度的變化特征及擬合方法，采用如下方法對(duì)飛行高度參數(shù)的插值進(jìn)行簡化。

首先，在飛行高度H對(duì)流場(chǎng)參數(shù)的影響規(guī)律方面，研究發(fā)現(xiàn)飛行器物面上的當(dāng)?shù)亓鲌?chǎng)參數(shù)與來流參數(shù)的比值在不同飛行高度下幾乎保持不變[28]。根據(jù)這一規(guī)律，對(duì)同一飛行器，若已知某高度H0下的邊界層外緣參數(shù)(以P0表示)以及該高度下的流動(dòng)參數(shù)Pinf,0，且已知任意高度下的大氣參數(shù)Pinf,x，那么就可以通過Px=Pinf,xP0/Pinf,0計(jì)算獲得該飛行器在任意高度上的邊界層外緣參數(shù)[24,28]。

由此一來，無黏外流解的計(jì)算可僅考慮馬赫數(shù)和迎角兩個(gè)參數(shù)。具體方法為：在對(duì)整個(gè)彈道狀態(tài)進(jìn)行計(jì)算時(shí)，無黏外流的求解只需計(jì)算設(shè)定參考飛行高度下的兩個(gè)馬赫數(shù)和兩個(gè)迎角，共四個(gè)飛行狀態(tài)，其他飛行狀態(tài)可根據(jù)插值方法快速得到。這四個(gè)計(jì)算狀態(tài)分別為設(shè)定飛行高度下的兩個(gè)馬赫數(shù)和兩個(gè)迎角的組合，即：(Ma1,α1),(Ma1,α2),(Ma2,α1),(Ma2,α2)。與這四個(gè)狀態(tài)所對(duì)應(yīng)的飛行器表面上任一點(diǎn)處的流動(dòng)參數(shù)的值分別記作P11、P12、P21和P22，則某個(gè)彈道時(shí)間點(diǎn)上飛行條件(Max,αx)狀態(tài)下的飛行器表面同一位置處的流動(dòng)參數(shù)Pxx可通過如下插值算法得到：

(9)

構(gòu)造該插值算法的目的在于減少彈道狀態(tài)無黏外流解的計(jì)算次數(shù)以提高數(shù)值仿真效率，該算法具有簡單可行、易實(shí)現(xiàn)等特點(diǎn)。由式(9)表示的插值算法可看出，在流場(chǎng)參數(shù)的插值精度方面主要依賴于已經(jīng)求得的無黏外流解的精度及插值變量間距(即馬赫數(shù)與迎角的增量值)。因此在實(shí)際計(jì)算中根據(jù)彈道參數(shù)特征合理建立插值狀態(tài)數(shù)據(jù)庫及采用高精度無黏外流求解器將有利于提高流場(chǎng)參數(shù)的插值精度，能夠最大程度地滿足工程設(shè)計(jì)中的誤差要求。該方法在插值精度與誤差分析方面的詳細(xì)分析參見文獻(xiàn)[23]。

4 高溫化學(xué)非平衡效應(yīng)估算

計(jì)算模型中，將無黏外流解得到的物面流動(dòng)參數(shù)作為邊界層外緣參數(shù)，即該無黏流的物面流動(dòng)參數(shù)不是實(shí)際黏性流場(chǎng)中飛行器的物面參數(shù)。由前述可知，實(shí)際物面的熱流密度是由Fay-Riddell熱流密度計(jì)算公式(1)得到的[29]。需要說明的是，該公式是在平衡邊界層的前提條件下推導(dǎo)得出的，因而未考慮高溫條件下的空氣化學(xué)非平衡效應(yīng)。在實(shí)際高超聲速飛行條件下，空氣的高溫非平衡效應(yīng)對(duì)氣動(dòng)熱影響顯著，因此對(duì)于實(shí)際高溫化學(xué)非平衡流動(dòng)，需要對(duì)Fay-Riddell駐點(diǎn)熱流密度計(jì)算公式進(jìn)行修正。式(10)給出了高溫空氣非平衡邊界層駐點(diǎn)傳熱與平衡邊界層駐點(diǎn)傳熱表面熱通量的關(guān)系式[30]：

(10)

式中，下標(biāo)O、N分別表示氧原子和氮原子，φ為壁面催化因子，h0為生成焓，CO,s和CN,s分別為氧原子和氮原子的質(zhì)量濃度，Le為路易斯數(shù)。

通過式(10)可看出，求解高溫化學(xué)非平衡流狀態(tài)下的駐點(diǎn)熱流，需要確定駐點(diǎn)位置氧原子和氮原子的質(zhì)量濃度。本文主要發(fā)展的是一種面向工程應(yīng)用的快速氣動(dòng)加熱計(jì)算技術(shù)，因此在組元參數(shù)特性的獲取方面，使用目前國內(nèi)外廣泛認(rèn)可的組元參數(shù)簡化計(jì)算模型，而不使用耗費(fèi)資源巨大的精確多組元N-S方程數(shù)值解。依據(jù)Dunn和Kang[31]發(fā)展的各組元濃度計(jì)算模型，在流場(chǎng)溫度9000 K以下時(shí)，可以只考慮O2、O、O+、N2、N、N+、NO、NO+、e-等九種組元及以下六個(gè)化學(xué)反應(yīng)式：

其中Ki為摩爾密度平衡常數(shù)，對(duì)于不同的化學(xué)反應(yīng)式，取值可由文獻(xiàn)[31]附表查得。

由高溫空氣化學(xué)反應(yīng)特性可知：當(dāng)高溫空氣在9000 K以下時(shí)化學(xué)反應(yīng)以離解反應(yīng)為主，電離反應(yīng)可忽略，在混合氣體組成中，NO、NO+、e-、N+、O+的濃度比O2、O、N2、N的濃度小得多，即：

(11)

(12)

因此，在化學(xué)反應(yīng)效應(yīng)對(duì)氣動(dòng)加熱的影響特性分析方面，暫不考慮NO、NO+、N+、O+、e-等5種組元。同時(shí)，氧原子和氮原子的壁面催化因子有如下關(guān)系成立：

(13)

(14)

其中[O2]0、[N2]0分別表示為空氣中氧分子和氮分子的摩爾密度。聯(lián)立式(13)、式(14)與反應(yīng)方程①②求解，可以分別得到氧原子和氮原子的摩爾密度為：

(15)

(16)

由于φO、φN遠(yuǎn)小于1，因此在一級(jí)近似時(shí)可以取φO=φN=0。然后聯(lián)立反應(yīng)方程①～③，可得到[O2]、[N2]和[NO]的摩爾密度計(jì)算式為：

(17)

最后聯(lián)立反應(yīng)方程④～⑥和電荷守恒方程[19]可以得到電子的摩爾密度為：

(18)

在氣動(dòng)加熱計(jì)算中，各組元質(zhì)量分?jǐn)?shù)需根據(jù)外流流動(dòng)特征采用時(shí)間歷程相關(guān)的迭代方法計(jì)算獲得。由此，考慮高溫化學(xué)非平衡效應(yīng)下的氣動(dòng)加熱計(jì)算方法可修正為：首先由Fay-Riddell駐點(diǎn)熱流密度計(jì)算式(1)，得到平衡條件下熱流密度qeq，然后根據(jù)化學(xué)反應(yīng)計(jì)算模型得到各組元摩爾密度和質(zhì)量分?jǐn)?shù)，最后由駐點(diǎn)化學(xué)非平衡邊界層與平衡邊界層熱流密度關(guān)系式(10)，求得高溫非平衡條件下的熱流密度qne。

研究中將如上所述的高溫空氣化學(xué)非平衡多組元效應(yīng)嵌入到所發(fā)展的全機(jī)外形飛行器彈道飛行狀態(tài)下的氣動(dòng)加熱與結(jié)構(gòu)傳熱數(shù)值計(jì)算方法中，使得計(jì)算模型更接近于高超聲速飛行器真實(shí)飛行條件下的熱環(huán)境。

5 算例與分析

根據(jù)以上理論分析與數(shù)值建模過程，發(fā)展了一種復(fù)雜飛行器氣動(dòng)加熱與結(jié)構(gòu)傳熱耦合計(jì)算技術(shù)。該技術(shù)主要由無黏外流解、氣動(dòng)熱工程計(jì)算和結(jié)構(gòu)傳熱三部分構(gòu)成，如圖2所示。

圖2 高超聲速氣動(dòng)加熱方法示意圖

在實(shí)際計(jì)算中，無黏外流解模塊需根據(jù)計(jì)算任務(wù)要求完成所需無黏流動(dòng)狀態(tài)的計(jì)算，并形成后續(xù)計(jì)算所需的外流解數(shù)據(jù)庫。邊界層加熱和結(jié)構(gòu)傳熱的計(jì)算可根據(jù)飛行器實(shí)際飛行時(shí)間進(jìn)行耦合求解。圖3給出了計(jì)算技術(shù)總體流程示意圖。

圖3 求解流程示意圖

下面采用三類典型的高超聲速氣動(dòng)加熱與結(jié)構(gòu)傳熱算例來對(duì)所發(fā)展的耦合計(jì)算方法的效率、功能及精度進(jìn)行具體分析。

5.1 典型外形氣動(dòng)加熱驗(yàn)證

5.1.1 鈍頭體瞬態(tài)氣動(dòng)加熱

本算例計(jì)算模型(圖4)選取NASA報(bào)告[32]中的三維鈍錐模型，模型長為0.352 m，頭部半徑0.027 94 m，錐角30°。無黏外流求解所用的計(jì)算網(wǎng)格為塊結(jié)構(gòu)網(wǎng)格，網(wǎng)格總單元數(shù)約32萬，物面網(wǎng)格單元約1.1萬。來流條件設(shè)定為：馬赫數(shù)Ma∞=10.6，迎角α=20°，壓強(qiáng)P∞=132.1 Pa，溫度T∞=47.3 K；設(shè)定鈍頭體壁溫Tw=294.44 K。

該算例在普通微機(jī)(3.3 GHz/4 GB，下同)上約在8 s CPU機(jī)時(shí)(不包含無黏外流解計(jì)算時(shí)間，下同)內(nèi)即可完成瞬態(tài)氣動(dòng)加熱的計(jì)算并輸出飛行器表面熱流分布。

圖4 計(jì)算模型

圖5為該算例物面熱流密度計(jì)算結(jié)果。圖5(a)為模型表面總熱流密度分布云圖?？梢娫谟杏乔闆r下，計(jì)算所得的沿子午面熱流密度分布具有很好的對(duì)稱性，而且熱流密度等值線圖分布均勻，沒有出現(xiàn)明顯的鋸齒、跳躍等現(xiàn)象。圖5(b)為采用駐點(diǎn)熱流值進(jìn)行歸一化處理后的子午面(截面方位角0°與180°)上熱流密度分布與試驗(yàn)值[32]的比較，可以看到計(jì)算結(jié)果與試驗(yàn)值相符。由于飛行器大面積氣動(dòng)熱數(shù)據(jù)的試驗(yàn)值難以獲得，因此在熱流具體數(shù)值的對(duì)比方面，取駐點(diǎn)作為對(duì)比驗(yàn)證點(diǎn)。表2給出了當(dāng)前流動(dòng)狀態(tài)下駐點(diǎn)熱流密度的計(jì)算結(jié)果與參考試驗(yàn)值對(duì)比，相對(duì)誤差小于10%，基本滿足氣動(dòng)熱工程設(shè)計(jì)的誤差要求。

(a)表面熱流密度分布

(b)子午面熱流對(duì)比

表2 駐點(diǎn)熱流值對(duì)比

5.1.2 RAM-CII巡航狀態(tài)長時(shí)加熱

計(jì)算外形采用典型高速飛行器RAM-CII球錐體試驗(yàn)?zāi)Ｐ?，幾何參?shù)為：頭部曲率半徑0.152 m，半頂角9°，模型長度1.3 m。無黏外流場(chǎng)解的計(jì)算網(wǎng)格為塊結(jié)構(gòu)網(wǎng)格，總單元數(shù)約40萬，物面單元約1.5萬。設(shè)定的巡航狀態(tài)為：飛行高度H=60 km，迎角α=0°，Ma∞=6.0，初始壁溫Tw=247 K，飛行時(shí)間t=1000 s。長時(shí)加熱中考慮結(jié)構(gòu)傳熱，飛行器熱防護(hù)結(jié)構(gòu)材料為金屬Ti，厚度2 mm，表面發(fā)射率為0.8。結(jié)構(gòu)傳熱計(jì)算中，時(shí)間步長取為0.05 s(即總時(shí)間推進(jìn)步數(shù)為2萬步)，采用熱薄壁傳熱模型。該算例在普通微機(jī)上計(jì)算耗時(shí)約1020 s(17 min)CPU機(jī)時(shí)。

(a)外表面溫度分布

(c)子午線溫度分布對(duì)比

通過圖6(c)中的對(duì)比可以看出，在子午線溫度分布上，計(jì)算結(jié)果與輻射平衡溫度吻合較好，當(dāng)時(shí)間足夠長時(shí)，表面溫度計(jì)算結(jié)果愈趨近于輻射平衡溫度，與理論分析一致。

5.2 復(fù)雜飛行器巡航狀態(tài)長時(shí)傳熱

計(jì)算模型為如圖7所示的類X-37B外形，圖中同時(shí)給出了計(jì)算所用的表面網(wǎng)格。計(jì)算狀態(tài)為巡航狀態(tài)長時(shí)傳熱：Ma∞=5，迎角α=3°，飛行高度H=40 km，飛行時(shí)間t=1000 s。用于無黏外流求解的塊結(jié)構(gòu)網(wǎng)格總單元數(shù)約512萬，物面網(wǎng)格單元約25.1萬。在熱防護(hù)結(jié)構(gòu)方面，該算例進(jìn)行了更接近于工程設(shè)計(jì)的復(fù)雜方案，即飛行器頭部區(qū)域設(shè)定為熱防護(hù)區(qū)1，以 TPS1表示，全機(jī)其他部位設(shè)為熱防護(hù)區(qū)2，以TPS2表示(見圖7)，不同的熱防護(hù)區(qū)域使用表3給出的不同的熱防護(hù)方案，其中熱防護(hù)結(jié)構(gòu)材料的最外層(飛行器表面)和最內(nèi)層溫度初值設(shè)為飛行高度上的大氣環(huán)境溫度245 K。根據(jù)該算例的飛行條件，不考慮高溫非平衡效應(yīng)影響特性。該算例耗時(shí)約1680 s(28 min)CPU機(jī)時(shí)。

圖7 計(jì)算模型與TPS方案

表3 熱防護(hù)系統(tǒng)參數(shù)設(shè)置

圖8給出了計(jì)算結(jié)果，其中圖8(a)和圖8(b)分別為第1000 s時(shí)飛行器防熱層外表面(Tw,surf)和內(nèi)表面(Tw,in)溫度分布，飛行器在計(jì)算初始時(shí)刻(第0 s)和計(jì)算結(jié)束時(shí)(第1000 s)的防護(hù)層外表面熱流密度分布如圖8(c)、圖8(d)所示。由計(jì)算結(jié)果可見，巡航1000 s時(shí)，類X-37B熱防護(hù)層內(nèi)表面溫度在不同的熱防護(hù)區(qū)域出現(xiàn)明顯差異(圖8(b))，而熱防護(hù)層整個(gè)外表面溫度分布均勻。這是由于TPS1區(qū)使用的是熱薄壁與隔熱層的組合防熱結(jié)構(gòu)，并且隔熱層采用熱沉性好的二氧化硅(SiO2)材料，故在該熱防護(hù)區(qū)域內(nèi)表面溫度出現(xiàn)顯著降低，大部分區(qū)域溫度在410 K左右，說明防熱層隔熱效果良好。而TPS2區(qū)僅采用了熱薄壁結(jié)構(gòu)，防熱材料為厚度0.002 m的金屬Ti，其熱傳導(dǎo)性很好，故防護(hù)層內(nèi)外表面溫度很快達(dá)到平衡，與外表面溫度差別很小，且防護(hù)區(qū)域熱流密度很小，該區(qū)域溫度大部分在600 K以上。圖8(e)、8(f)分別為駐點(diǎn)內(nèi)外層溫度與熱流密度隨時(shí)間變化曲線，可以發(fā)現(xiàn)TPS1區(qū)域內(nèi)外表面溫度變化差別較大，計(jì)算結(jié)果顯示駐點(diǎn)區(qū)域防熱層外表面最高溫度約為979 K，內(nèi)表面最高溫度約為521 K，造成此溫度值差異的原因仍然是由于在這個(gè)區(qū)域使用了SiO2隔熱層。算例計(jì)算與分析結(jié)果表明，發(fā)展的復(fù)雜飛行器氣動(dòng)加熱與結(jié)構(gòu)傳熱耦合計(jì)算方法可用于復(fù)雜外形及多種熱防護(hù)方案問題的求解。

(a)第1000 s 時(shí)外表面溫度

(b)第1000 s 時(shí)內(nèi)表面溫度

(c)初始時(shí)刻外表面熱流

(d)第1000 s外表面熱流

(e)駐點(diǎn)內(nèi)/外層溫度

(f)駐點(diǎn)熱流

5.3 復(fù)雜飛行器彈道狀態(tài)長時(shí)加熱

文中的彈道狀態(tài)長時(shí)加熱指的是飛行器按照給定的彈道狀態(tài)(變高度、迎角、馬赫數(shù))飛行時(shí)的氣動(dòng)加熱與結(jié)構(gòu)傳熱特性，彈道計(jì)算狀態(tài)更接近飛行器實(shí)際飛行包線狀態(tài)。計(jì)算模型采用算例5.2的模型及計(jì)算網(wǎng)格。由于考慮的是彈道狀態(tài)，因此使用了前文中的彈道狀態(tài)動(dòng)態(tài)插值方法。具體計(jì)算方法流程如圖9所示。

由于缺乏相關(guān)彈道參數(shù)數(shù)據(jù)，本文根據(jù)相關(guān)文獻(xiàn)設(shè)定的彈道參數(shù)如表4所示。彈道參數(shù)考慮了馬赫數(shù)、飛行高速及迎角的變化，彈道飛行時(shí)間1000 s，在此期間飛行高度由60 km 降至30 km，飛行馬赫數(shù)由Ma=6降至Ma=4，迎角變化范圍為±5°。該算例全機(jī)表面采用同一種組合熱防護(hù)方案模型，具體參數(shù)與算例5.2中 TPS1區(qū)相同。該算例耗時(shí)約1860 s(31 min)CPU機(jī)時(shí)。

表4 彈道參數(shù)

圖10給出了該算例計(jì)算結(jié)果。從圖10中可以清晰地看到飛行器表面沿彈道飛行時(shí)逐漸加熱到平衡狀態(tài)的過程，并且隨迎角從-5°至5°變化過程中,飛行器頭部與機(jī)翼前緣處的高溫區(qū)往下表面移動(dòng)，第1000 s時(shí)，最高溫出現(xiàn)在飛行器頭部下表面，約為1050 K，這與高超聲速飛行器防熱結(jié)構(gòu)材料特性理論分析相符。由于沒有在公開發(fā)表文獻(xiàn)中找到與此類似的算例進(jìn)行對(duì)比，因此這里僅給出本文算例計(jì)算結(jié)果的分析。該算例表明所發(fā)展的計(jì)算方法可對(duì)復(fù)雜高超聲速飛行器在彈道飛行狀態(tài)下的氣動(dòng)加熱與結(jié)構(gòu)傳熱耦合進(jìn)行快速高效計(jì)算與分析，可為高速飛行器的熱防護(hù)設(shè)計(jì)與初期選型提供技術(shù)參考。

(a)t=0 s

(b)t=100 s

(c)t=350 s

(d)t=600 s

(e)t=1000 s

6 結(jié) 論

針對(duì)復(fù)雜外形飛行器在復(fù)雜飛行條件下的氣動(dòng)加熱與結(jié)構(gòu)傳熱耦合問題，發(fā)展了一種基于耦合邊界層外無黏流場(chǎng)解與氣動(dòng)熱工程方法的高超聲速全機(jī)外形氣動(dòng)加熱與結(jié)構(gòu)傳熱快速計(jì)算方法。對(duì)典型外形及典型狀態(tài)進(jìn)行了數(shù)值計(jì)算與分析，得到了與理論分析及文獻(xiàn)相符的計(jì)算結(jié)果。結(jié)論如下：

1)以邊界層理論為基礎(chǔ)，綜合數(shù)值計(jì)算方法與工程算法各自的優(yōu)勢(shì)，將高超聲速飛行器氣動(dòng)熱的計(jì)算簡化為飛行器無黏外流歐拉數(shù)值解和邊界層內(nèi)熱流求解，并考慮了熱化學(xué)非平衡效應(yīng)。同時(shí)在計(jì)算方法中嵌入不同防熱模型(熱薄壁、熱厚壁)的結(jié)構(gòu)傳熱計(jì)算方法，可用于快速數(shù)值計(jì)算分析熱防護(hù)系統(tǒng)中防熱材料的溫度場(chǎng)時(shí)變特性。

2)所發(fā)展的計(jì)算方法避開了全流場(chǎng)黏性數(shù)值模擬所需的巨大計(jì)算量，同時(shí)拓展了氣動(dòng)熱工程估算方法的應(yīng)用范圍，計(jì)算精度可滿足工程設(shè)計(jì)初期方案論證要求，另外嵌入彈道狀態(tài)動(dòng)態(tài)插值方法，可方便地應(yīng)用于高超聲速復(fù)雜全機(jī)外形在復(fù)雜彈道飛行條件下的氣動(dòng)加熱與結(jié)構(gòu)傳熱多場(chǎng)耦合問題的快速計(jì)算分析。

3)所發(fā)展的方法具有較好的可移植性，在后續(xù)研究中將考慮不同防熱材料及材料燒蝕效應(yīng)對(duì)高超聲速飛行器氣動(dòng)加熱特性的影響。