黃 杰，姚衛(wèi)星，陳 炎,3，孔 斌,3

(1. 南京航空航天大學(xué)飛行器先進(jìn)設(shè)計(jì)技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，南京 210016；2. 南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京 210016；3. 中國(guó)航空工業(yè)集團(tuán)公司成都飛機(jī)設(shè)計(jì)研究所，成都 610091)

0 引 言

空天飛行器再入階段受到巨大的氣動(dòng)加熱作用，為保證內(nèi)部機(jī)體結(jié)構(gòu)在可承受的溫度范圍內(nèi)，需要在機(jī)體表面附加熱防護(hù)系統(tǒng)(Thermal protection system,TPS)[1-3]，熱防護(hù)系統(tǒng)的設(shè)計(jì)為高超聲速飛行器關(guān)鍵技術(shù)之一，而陶瓷防熱瓦是應(yīng)用最廣泛的隔熱結(jié)構(gòu)之一，它通過(guò)應(yīng)變隔離墊粘接在機(jī)體表面，其熱控性能直接影響到飛行器的安全，必須進(jìn)行熱防護(hù)系統(tǒng)的熱分析。

傳統(tǒng)的熱防護(hù)系統(tǒng)熱控分析方法將氣動(dòng)熱與結(jié)構(gòu)傳熱人為地分割開(kāi)來(lái)。文獻(xiàn)[4-7]在已知熱流密度的情況下進(jìn)行了飛行器頭部防熱瓦、蓋板式熱防護(hù)系統(tǒng)以及多層隔熱結(jié)構(gòu)的熱分析，獲得了熱防護(hù)系統(tǒng)和機(jī)體結(jié)構(gòu)的溫度-時(shí)間歷程。其求解過(guò)程為先進(jìn)行氣動(dòng)加熱的分析，再將熱流密度作為邊界條件直接分析傳熱方程求解熱防護(hù)系統(tǒng)的溫度場(chǎng)，而忽略氣動(dòng)熱與結(jié)構(gòu)傳熱之間的耦合效應(yīng)，這必然會(huì)影響熱防護(hù)系統(tǒng)熱控性能的分析精度。

由于空天飛行器再入過(guò)程中產(chǎn)生巨大的氣動(dòng)加熱會(huì)造成結(jié)構(gòu)溫度急劇升高，而結(jié)構(gòu)溫度升高后，激波層和邊界層內(nèi)氣體與熱防護(hù)系統(tǒng)外表面之間的溫度梯度將會(huì)減小，造成壁面熱流密度的降低，即氣動(dòng)加熱與結(jié)構(gòu)溫度場(chǎng)存在強(qiáng)烈的耦合效應(yīng)，應(yīng)給予足夠的重視。

耦合分析的關(guān)鍵在于氣動(dòng)熱的準(zhǔn)確計(jì)算、耦合迭代的策略以及耦合面之間的數(shù)據(jù)插值。以往的研究常常采用Eckert[8]提出的參考焓法等工程算法進(jìn)行熱流密度的分析，其分析精度有限，現(xiàn)代氣動(dòng)熱分析采用計(jì)算流體力學(xué)(Computational fluid dynamics,CFD)方法求解Navier-Stokes方程[9-10]，可借助CFD++或FLUENT等軟件或自編程序準(zhǔn)確分析駐點(diǎn)及其他區(qū)域的熱流密度分布。而耦合分析常采用迭代求解法，其核心技術(shù)是要保證瞬態(tài)耦合的時(shí)間精度和迭代的協(xié)調(diào)性。此外由于結(jié)構(gòu)傳熱網(wǎng)格單元的尺度遠(yuǎn)大于流體單元，耦合面節(jié)點(diǎn)非一一對(duì)應(yīng)，故耦合面之間的數(shù)據(jù)交換需要采用插值法完成，其關(guān)鍵技術(shù)是要保證插值精度和耦合面熱流通量的守恒性。

本文在充分考慮氣動(dòng)熱與結(jié)構(gòu)傳熱耦合的基礎(chǔ)上，提出了熱防護(hù)系統(tǒng)熱控性能分析的分區(qū)協(xié)調(diào)耦合推進(jìn)方法，對(duì)比了圓管算例的計(jì)算結(jié)果與試驗(yàn)結(jié)果，最后進(jìn)行了防熱瓦與主動(dòng)冷卻系統(tǒng)同時(shí)存在時(shí)的熱防護(hù)系統(tǒng)熱控性能及其影響因素分析。

1 流體和結(jié)構(gòu)傳熱的數(shù)值算法

1.1 流體分析數(shù)值算法

不考慮體積力和內(nèi)熱源情況下，直角坐標(biāo)系下的流體動(dòng)力學(xué)N-S方程的積分形式為：

(1)

式中：W為守恒向量，F(xiàn)c為對(duì)流通量，F(xiàn)v為黏性通量，dS為控制體V的邊界面，n為邊界面dS的外法線單位向量。將式(1)按有限體積法(Finite volume method,FVM)進(jìn)行空間離散可得：

(2)

式中：Wi和Vi分別為控制體i的守恒向量和體積，NF為控制體邊界面的數(shù)目，ΔSN為第N個(gè)邊界面的面積。由于對(duì)流通量Fc具有高度非線性特點(diǎn)，并且集中體現(xiàn)了流場(chǎng)的對(duì)流特征，采用具有TVD(Total Variation Diminishing)性質(zhì)的NND格式[11]對(duì)其進(jìn)行空間離散，其中半離散化的上風(fēng)型NND格式為：

(3)

(4)

(5)

(6)

式中：sgn為符號(hào)函數(shù)，min表示取較小值。

控制方程中的黏性項(xiàng)采用中心差分格式離散，湍流模型采用Menter’s SST兩方程模型[12]，而針對(duì)非定常問(wèn)題的時(shí)間離散，在n+1時(shí)刻采用時(shí)間二階精度的隱式三點(diǎn)后差離散可得到二階精度的離散方程：

(7)

(8)

式中：Δτ和Δt分別為虛擬時(shí)間步長(zhǎng)和物理時(shí)間步長(zhǎng)，稱(chēng)為雙時(shí)間步長(zhǎng)法[13]；p和n分別為虛擬時(shí)間迭代步和物理時(shí)間迭代步。虛擬時(shí)間步上的內(nèi)迭代可采用LU-SGS格式[14]求解，當(dāng)p→∞時(shí)虛擬時(shí)間項(xiàng)趨近于零，式(8)的定常解即為二階精度的非定常解。

1.2 結(jié)構(gòu)傳熱數(shù)值算法

在無(wú)體積熱源的假設(shè)下結(jié)構(gòu)瞬態(tài)熱傳導(dǎo)的控制方程為：

(9)

式中：ρ0為結(jié)構(gòu)材料密度，c為材料比熱容，kx，ky和kz分別為材料三個(gè)方向的導(dǎo)熱系數(shù)。其中比熱容和導(dǎo)熱系數(shù)一般為溫度的函數(shù)。

針對(duì)本文熱防護(hù)系統(tǒng)的熱分析問(wèn)題，其外表面邊界條件為壁面熱流密度Qaero和壁面熱輻射量Qrad，其表達(dá)式分別為：

(10)

(11)

式中：?T/?n為壁面法向溫度梯度，ε為壁面熱輻射率，玻爾茲曼常數(shù)σ=5.67×10-8W/(m2·K4)，Twall為壁面溫度，Tat為大氣環(huán)境溫度。

對(duì)式(9)進(jìn)行有限元離散可得到總體合成矩陣求解方程：

(12)

(13)

求解式(13)即可得到結(jié)構(gòu)各個(gè)時(shí)刻節(jié)點(diǎn)的瞬時(shí)溫度。

2 耦合分析模型和策略

空天飛行器再入過(guò)程中熱防護(hù)系統(tǒng)外壁面巨大的氣動(dòng)加熱效應(yīng)會(huì)造成結(jié)構(gòu)溫度急劇升高，而結(jié)構(gòu)溫度升高后，激波層和邊界層內(nèi)氣體與壁面之間的溫度梯度將會(huì)減小，造成壁面熱流密度的降低，即氣動(dòng)加熱與結(jié)構(gòu)溫度場(chǎng)存在強(qiáng)烈的相互耦合效應(yīng)，如圖1所示。

本文采用分區(qū)協(xié)調(diào)耦合推進(jìn)方法進(jìn)行氣動(dòng)熱與結(jié)構(gòu)傳熱的分析，其耦合推進(jìn)策略如圖2所示，即將流場(chǎng)與結(jié)構(gòu)分為兩個(gè)不同的區(qū)域，分開(kāi)建模與求解，其中流場(chǎng)采用FVM求解，而結(jié)構(gòu)熱傳導(dǎo)采用有限元法(Finite element method,FEM)求解，其基本假設(shè)和特點(diǎn)為：

1) 在任意n到n+1時(shí)間步內(nèi)FVM(或FEM)求解過(guò)程中壁面溫度(或壁面熱流密度)不變，即各子學(xué)科求解時(shí)邊界條件凍結(jié)。

2) 各子學(xué)科在相同的時(shí)間節(jié)點(diǎn)進(jìn)行數(shù)據(jù)交換，以保證耦合分析的協(xié)調(diào)性以及耦合時(shí)間精度。

根據(jù)圖2所示的耦合推進(jìn)策略，本文提出了適用于瞬態(tài)和穩(wěn)態(tài)的熱防護(hù)系統(tǒng)熱分析流程(見(jiàn)圖3)，主要步驟為：

1) 首先建立相應(yīng)的FVM和FEM數(shù)值分析模型，定義來(lái)流條件和初始溫度T0，進(jìn)行定常氣動(dòng)熱分析，結(jié)果作為耦合分析的熱流密度初始條件。

2) 進(jìn)行第n到n+1時(shí)間步的流場(chǎng)求解，將計(jì)算獲得的熱流密度Qn+1通過(guò)插值方法傳遞給結(jié)構(gòu)FEM模型外表面。

3) 進(jìn)行第n到n+1時(shí)間步的結(jié)構(gòu)FEM傳熱分析，輸出分析獲得的結(jié)構(gòu)溫度場(chǎng)Tn+1。

4) 判斷時(shí)間是否達(dá)到分析的總時(shí)間ttotal(針對(duì)瞬態(tài)分析)，或壁面熱流密度和結(jié)構(gòu)溫度場(chǎng)是否收斂(針對(duì)穩(wěn)態(tài)分析)，若達(dá)到ttotal或收斂，結(jié)束耦合分析，否則進(jìn)行熱流密度和壁面溫度的數(shù)據(jù)傳遞。

5) 返回2)，進(jìn)行下一個(gè)迭代步的求解，直到分析時(shí)間達(dá)到ttotal或壁面熱流密度和結(jié)構(gòu)溫度場(chǎng)收斂，結(jié)束分析。

3 耦合面數(shù)據(jù)傳遞

以上氣動(dòng)熱和結(jié)構(gòu)傳熱的耦合分析中涉及到熱流密度和壁面溫度的數(shù)據(jù)傳遞，若這兩個(gè)學(xué)科之間的數(shù)值模型網(wǎng)格節(jié)點(diǎn)一一對(duì)應(yīng)，那么通過(guò)節(jié)點(diǎn)之間的對(duì)應(yīng)關(guān)系就可以方便地進(jìn)行熱流密度和壁面溫度的數(shù)據(jù)傳遞。但實(shí)際分析中結(jié)構(gòu)網(wǎng)格尺度通常遠(yuǎn)大于流體網(wǎng)格，即FVM和FEM耦合面節(jié)點(diǎn)不一致，如圖4所示，故需要在耦合面進(jìn)行插值以實(shí)現(xiàn)數(shù)據(jù)傳遞。

本文的分析涉及到熱流密度和結(jié)構(gòu)壁面溫度的數(shù)據(jù)插值，其中的關(guān)鍵技術(shù)是要保證耦合面上熱流密度通量的守恒性。

(14)

式中：Q為FVM網(wǎng)格面的熱流密度，q為與FVM網(wǎng)格節(jié)點(diǎn)對(duì)應(yīng)的結(jié)構(gòu)傳熱FEM網(wǎng)格面的熱流密度。當(dāng)Scoup為耦合面時(shí)，式(14)體現(xiàn)了熱流密度通量的總體守恒性。本文采用如下基于控制面的雙向映射插值方法進(jìn)行熱流密度的數(shù)據(jù)傳遞，算法主要步驟為：

1) 將FVM和FEM耦合面上節(jié)點(diǎn)從物理空間(x,y,z)通過(guò)坐標(biāo)轉(zhuǎn)換映射到參數(shù)空間(u,v)中，即u=u(x,y,z)和v=v(x,y,z)，它將三維曲面上的空間節(jié)點(diǎn)轉(zhuǎn)換到二維的控制面上。

2) 在物理空間中采用Kdtree算法搜索任意FEM網(wǎng)格節(jié)點(diǎn)ζi(x,y,z)附近的FVM網(wǎng)格節(jié)點(diǎn)ηi(x,y,z)，并將其轉(zhuǎn)換到二維的控制面得到ζi(u,v)和ηi(u,v)。

3) 將FVM網(wǎng)格節(jié)點(diǎn)坐標(biāo)ηi(u,v)和相應(yīng)的熱流Qi(u,v)代入到如下三次插值函數(shù)Q(u,v)，采用最小二乘法求解插值函數(shù)的系數(shù)bj(j=1,2,…,10)。

Q(u,v)=b1u3+b2v3+b3u2v+b4v2u+

b5u2+b6v2+b7uv+b8u+b9v+b10

(15)

4) 將所有FEM網(wǎng)格節(jié)點(diǎn)ζi(u,v)代入到已知系數(shù)bj(j=1,2,…,10)的插值函數(shù)Q(u,v)中，即可求得FEM網(wǎng)格節(jié)點(diǎn)插值熱流密度qi(u,v)，并通過(guò)式(14)進(jìn)行熱流密度通量的守恒性檢驗(yàn)。

4 算例1：耦合方法的驗(yàn)證

本文選取NASA高超聲速圓管風(fēng)洞試驗(yàn)?zāi)Ｐ蚚15]進(jìn)行算例分析。圓管內(nèi)徑R1=25.4 mm，外徑R2=38.1 mm，圓管材料為不銹鋼，其密度ρ=8030 Kg/m3，導(dǎo)熱系數(shù)k=16.72 W/(m·K)，比熱容c=502.48 J/(kg·K)，此外圓管初始壁面溫度T0=294.4 K。高超聲速來(lái)流馬赫數(shù)Ma=6.47，來(lái)流溫度T=241.5 K，來(lái)流壓強(qiáng)P=648.1 Pa，攻角α=0°。建立了二維分析模型，如圖5所示，其中流場(chǎng)網(wǎng)格量為60×100，壁面第一層網(wǎng)格高度Δh=1×10-5m，并且在激波處進(jìn)行了局部加密。圓管結(jié)構(gòu)有限元網(wǎng)格量為40×25，采用四節(jié)點(diǎn)平面單元模擬。此外分析類(lèi)型為瞬態(tài)分析，時(shí)間步長(zhǎng)Δt=1×10-4s，總時(shí)間ttotal=2 s，并且將定常流場(chǎng)作為耦合分析的初始條件。

分析獲得了圓管駐點(diǎn)溫度隨時(shí)間的變化曲線，如圖6所示，從圖中可觀察到曲線斜率在初始階段很大，隨時(shí)間的推移逐漸減小，說(shuō)明結(jié)構(gòu)開(kāi)始階段升溫很快，隨后溫度上升速度逐漸減緩，最后趨于一個(gè)穩(wěn)定的溫度值(穩(wěn)態(tài)解)。

圖7為2 s時(shí)刻圓管外壁面的相對(duì)溫度分析情況,從圖中可知壁面耦合分析結(jié)果分布曲線與試驗(yàn)結(jié)果分布曲線吻合良好，其中駐點(diǎn)處溫度耦合分析結(jié)果Tstag為442 K，而試驗(yàn)結(jié)果為465 K，相對(duì)誤差為4.95%。從而驗(yàn)證了分區(qū)協(xié)調(diào)耦合推進(jìn)方法的正確性與分析精度，其可用于分析熱防護(hù)系統(tǒng)氣動(dòng)熱與結(jié)構(gòu)傳熱的耦合問(wèn)題。

5 算例2：空天飛行器頭錐熱防護(hù)系統(tǒng)耦合分析

現(xiàn)代空天飛行器同時(shí)采用防熱瓦及主動(dòng)冷卻系統(tǒng)以保證蒙皮內(nèi)表面溫度穩(wěn)定在設(shè)定值，如圖8所示。本文進(jìn)行了空天飛行器頭錐熱防護(hù)系統(tǒng)的穩(wěn)態(tài)耦合分析。其中防熱瓦外表面涂層熱輻射率ε=0.85，LI-900防熱瓦厚度h1=60 mm，應(yīng)變隔離墊(SIP)厚度h2=4.5 mm，鋁合金蒙皮厚度h3=1.5 mm，蒙皮導(dǎo)熱系數(shù)k=180 W/(m·K)。防熱瓦和SIP的導(dǎo)熱系數(shù)均隨溫度呈非線性變化趨勢(shì)[16]。飛行馬赫數(shù)Ma=7，高度H=60 km，攻角α=25°，初始結(jié)構(gòu)溫度T0=300 K，大氣環(huán)境溫度Tat=247.1 K。建立了三維分析模型，流場(chǎng)和結(jié)構(gòu)網(wǎng)格如圖9所示，結(jié)構(gòu)熱分析采用八節(jié)點(diǎn)六面體單元模擬，蒙皮內(nèi)表面溫度約束在325 K，為此主動(dòng)冷卻系統(tǒng)需要持續(xù)的提供功率以保證蒙皮內(nèi)表面維持在此溫度。

利用本文的耦合方法進(jìn)行了空天飛行器頭錐熱防護(hù)系統(tǒng)的分析。圖10為非耦合和耦合方法分析獲得的對(duì)稱(chēng)壁面處的熱流密度分布曲線，可知非耦合分析方法獲得的駐點(diǎn)熱流密度比耦合方法高了53.8%。這是由于非耦合方法未考慮壁面溫度升高造成激波層和邊界層與壁面溫度梯度減小。

圖11和和12分別為非耦合和耦合方法分析得到的防熱瓦及SIP外表面溫度云圖。其中非耦合和耦合方法獲得的防熱瓦最高溫度分別為1108.5 K和994.1 K，相差了114.4 K；而非耦合和耦合方法獲得的SIP最高溫度分別為472.9 K和440.3 K，相差了32.6 K。故非耦合方法獲得的防熱瓦和SIP溫度偏高，這是由于非耦合方法考慮壁面溫度升高對(duì)氣動(dòng)熱的反饋?zhàn)饔迷斐傻摹?/p>

研究了飛行馬赫數(shù)Ma對(duì)防熱瓦和SIP最高溫度、主動(dòng)冷卻系統(tǒng)最大功率的影響，如圖13所示，從圖中可觀察到馬赫數(shù)從5增加到9，防熱瓦最高溫度上升了69.2%，SIP最高溫度上升了36.1%，而最大冷卻功率上升了266.1%。圖14為以上響應(yīng)隨熱防護(hù)系統(tǒng)外表面涂層熱輻射率ε的變化情況，熱輻射率從0.25增加到1，防熱瓦最高溫度降低了24.5%，SIP最高溫度降低了18.2%，而最大冷卻功率降低了50.1%，故應(yīng)盡量采用高輻射率的涂層以降低熱防護(hù)系統(tǒng)的厚度和主動(dòng)冷卻系統(tǒng)的功率，達(dá)到減重的目的。

在進(jìn)行熱防護(hù)系統(tǒng)的設(shè)計(jì)時(shí)，防熱瓦的導(dǎo)熱系數(shù)和厚度的選取直接影響到熱控性能，是熱防護(hù)系統(tǒng)最重要的設(shè)計(jì)參數(shù)。圖15和圖16分別為防熱瓦和SIP最高溫度、主動(dòng)冷卻系統(tǒng)最大功率隨防熱瓦導(dǎo)熱系數(shù)比k/k0以及厚度h1的變化曲線，k0代表防熱瓦原始導(dǎo)熱系數(shù)。從圖中可觀察到防熱瓦導(dǎo)熱系數(shù)比從0.128增加到8時(shí)，防熱瓦最高溫度降低了3.2%，而防熱瓦厚度從20 mm增加到80 mm，防熱瓦最高溫度上升了1.5%，即防熱瓦導(dǎo)熱系數(shù)和厚度對(duì)其最高溫度影響很小。隨著防熱瓦導(dǎo)熱系數(shù)的增加和厚度的降低，外部氣動(dòng)熱能量能夠更容易地傳遞到SIP及蒙皮，造成SIP溫度升高以及主動(dòng)冷卻系統(tǒng)功率增加，并且影響顯著。故在防熱瓦重量約束條件下應(yīng)盡量采取較厚和較低導(dǎo)熱系數(shù)的防熱瓦，以提高熱防護(hù)系統(tǒng)的隔熱性能和降低主動(dòng)冷卻系統(tǒng)的重量。

6 結(jié) 論

1) 本文針對(duì)空天飛行器熱防護(hù)系統(tǒng)熱控性能研究，提出了分區(qū)協(xié)調(diào)耦合推進(jìn)分析方法，其中氣動(dòng)熱和結(jié)構(gòu)傳熱分別采用FVM和FEM求解，且在耦合面進(jìn)行數(shù)據(jù)插值以實(shí)現(xiàn)數(shù)據(jù)傳遞。

2) 進(jìn)行了圓管算例分析，圓管壁面熱流密度和溫度分布與試驗(yàn)結(jié)果吻合良好，從而驗(yàn)證了本文的耦合推進(jìn)方法的正確性與分析精度。

3) 研究了空天飛行器頭錐熱防護(hù)系統(tǒng)的熱控性能，非耦合方法獲得的防熱瓦和SIP溫度相對(duì)耦合方法偏高，這是由于非耦合方法未考慮壁面溫度升高對(duì)氣動(dòng)熱的反饋?zhàn)饔?，而耦合方法充分考慮了此影響。

[1] Olynick D. Trajectory-based thermal protection system sizing for an X-33 winged vehicle concept [J]. Journal of Spacecraft and Rockets, 1998, 35(3): 249-257.

[2] Shideler J L, Webb G L, Pittman C M. Verification tests of durable thermal protection system concepts [J]. Journal of Spacecraft and Rockets, 1985, 22(6): 598-604.

[3] Oscar A M, Anurag S, Bhavani V S, et al. Thermal force and moment determination of an integrated thermal protection system [J]. AIAA Journal, 2010, 48(1): 119-128.

[4] 劉雙, 張博明. 發(fā)汗式主動(dòng)冷卻金屬熱防護(hù)系統(tǒng)主動(dòng)冷卻效率研究[J]. 宇航學(xué)報(bào), 2011, 32(2):433-438. [Liu Shuang, Zhang Bo-ming. Investigation on transpiration active cooling metallic thermal protection systems [J]. Journal of Astronautics, 2011, 32(2): 433-438.]

[5] 孟松鶴, 楊強(qiáng), 霍施宇, 等. 一體化熱防護(hù)技術(shù)現(xiàn)狀和發(fā)展趨勢(shì)[J]. 宇航學(xué)報(bào), 2013, 34(10):1295-1302. [Meng Song-he, Yang Qiang, Huo Shi-yu, et al. State-of-arts and trend of integrated thermal protection system [J]. Journal of Astronautics, 2013, 34(10): 1295-1302.]

[6] 魏東, 杜雁霞, 石友安, 等. 非均勻氣動(dòng)加熱下隔熱層結(jié)構(gòu)的優(yōu)化設(shè)計(jì)[J]. 宇航學(xué)報(bào), 2015, 36(10):1108-1113. [Wei Dong, Du Yan-xia, Shi You-an, et al. Optimization design of heat insulation layer with non-uniform heat flow loads [J]. Journal of Astronautics, 2015, 36(10): 1108-1113.]

[7] 朱言旦, 劉偉, 曾磊, 等. 飛行器結(jié)構(gòu)部件導(dǎo)熱/輻射耦合傳熱特性預(yù)測(cè)方法[J]. 宇航學(xué)報(bào), 2016, 37(11):1371-1377. [Zhu Yan-dan, Liu Wei, Zeng Lei, et al. Method of predicting conduction-radiation coupled heat transfer characteri-stics for vehicle structural component [J]. Journal of Astron-autics, 2016, 37(11): 1371-1377.]

[8] Eckert E R G. Engineering relations for friction and heat transfer to surfaces in high velocity flow [J]. Journal of Aerospace Sciences, 1955, 22(8): 585-587.

[9] Ju Y. Lower-upper scheme for chemically reacting flow with finite rate chemistry [J]. AIAA Journal, 1985, 33(8): 1418-1425.

[10] 楊肖峰, 唐偉, 桂業(yè)偉, 等. 火星環(huán)境高超聲速催化加熱特性[J]. 宇航學(xué)報(bào), 2017, 38(2):205-211. [Yang Xiao-feng, Tang Wei, Gui Ye-wei, et al. Hypersonic catalytic aeroheating characteristics for Mars entry process [J]. Journal of Astronautics, 2017, 38(2): 205-211.]

[11] 沈清, 顧鋼民, 高樹(shù)椿, 等. NND格式在航天飛機(jī)頭部段N-S方程求解中的應(yīng)用[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 1989, 7(2):145-155. [Shen Qing, Gu Gang-min, Gao Shu-chun, et al. Applications of the NND-scheme to the solving of Navier-Stokes equations on the fore-head of space-shuttle-like body [J]. Atca Aerodynamica Sinica, 1989, 7(2): 145-155.]

[12] Menter F R. Two-equation eddy-viscosity turbulence models for engineering applications [J]. AIAA Journal, 1994, 32(8): 1598-1605.

[13] 王保國(guó), 劉淑艷, 張雅, 等. 雙時(shí)間步長(zhǎng)加權(quán)ENO-強(qiáng)緊致高分辨率格式及在葉輪機(jī)械非定常流動(dòng)中的應(yīng)用[J]. 航空動(dòng)力學(xué)報(bào), 2005, 20(4):534-539. [Wang Bao-guo, Liu Shu-yan, Zhang Ya, et al. Pseudo-time method of weighted ENO-strongly compact high resolution schemes and its application to turbomachinery unsteady flow [J]. Journal of Aerospace Power, 2005, 20(4): 534-539.]

[14] Yoon S, Jameson A. Low-upper Gauss-Sediel method for the Euler and Navier-Stokes equations [J]. AIAA Journal, 1988, 26(9): 1025-1026.

[15] Wieting A R. Experimental study of shock wave interference heating on a cylindrical leading edge [R]. NASA-TM-100484, 1987.

[16] Ng W H, Friedmann P P, Waas A. Thermomechanical analysis of a thermal protection system with defects and heat shorts [R]. AIAA-2006-2212, 2006.