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

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

0 引 言

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

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

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

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

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

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

1 無黏外流場數(shù)值計算方法

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

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

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

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

2.1 物面熱流密度計算

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

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

(1)

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

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

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

(2)

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

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

(3)

初始條件為：

Tw|t=0=T0

(4)

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

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

表層：

(5)

最內(nèi)層：

(6)

中間層：

(7)

初始條件為：

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

(8)

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

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

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

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

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

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

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

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

(9)

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

4 高溫化學非平衡效應估算

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

(10)

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

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

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

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

(11)

(12)

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

(13)

(14)

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

(15)

(16)

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

(17)

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

(18)

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

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

5 算例與分析

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

圖2 高超聲速氣動加熱方法示意圖

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

圖3 求解流程示意圖

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

5.1 典型外形氣動加熱驗證

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

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

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

圖4 計算模型

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

(a)表面熱流密度分布

(b)子午面熱流對比

表2 駐點熱流值對比

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

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

(a)外表面溫度分布

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

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

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

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

圖7 計算模型與TPS方案

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

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

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

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

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

(d)第1000 s外表面熱流

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

(f)駐點熱流

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

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

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

表4 彈道參數(shù)

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

(a)t=0 s

(b)t=100 s

(c)t=350 s

(d)t=600 s

(e)t=1000 s

6 結(jié) 論

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

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

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

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