有限催化對返回艙氣動熱環(huán)境影響

粟斯堯, 石義雷, 柳 森, 彭治雨, 黎作武

(1. 中國空氣動力研究與發(fā)展中心 超高速空氣動力研究所, 四川 綿陽 621000； 2. 國家計算流體力學(xué)實驗室, 北京 100191)

0 引 言

飛船返回艙等高超聲速飛行器在大氣層中飛行時，頭部前方空氣經(jīng)激波強(qiáng)烈壓縮而急劇升溫，并將發(fā)生離解、電離等化學(xué)反應(yīng)，產(chǎn)生相應(yīng)原子和離子。而飛行器表面材料往往對流場中的原子、離子具有催化復(fù)合作用，使其在飛行器表面附近發(fā)生復(fù)合反應(yīng)形成相應(yīng)分子，由此改變飛行器表面附近的組分分布。

另一方面，催化復(fù)合反應(yīng)是一種放熱過程，它也會改變飛行器表面附近流場的溫度分布。顯然這兩方面因素都將直接影響飛行器的氣動熱環(huán)境。這一現(xiàn)象已被美國NASA航天飛機(jī)防熱瓦催化特性飛行試驗所證實[1]。因此研究壁面催化對飛行器熱環(huán)境影響規(guī)律，發(fā)展壁面催化條件下氣動熱環(huán)境精確預(yù)測方法對高超聲速飛行器熱防護(hù)系統(tǒng)設(shè)計和可靠性驗證具有重大意義。

對于高超聲速飛行器，壁面催化反應(yīng)往往發(fā)生在高溫非平衡流動中。受地面設(shè)備能力及經(jīng)費的限制，目前難以依靠地面或飛行試驗系統(tǒng)研究壁面催化對飛行器氣動熱環(huán)境影響。地面試驗主要用于研究材料的催化特性，為催化計算建模和校核提供依據(jù)[2]。隨著CFD技術(shù)的發(fā)展和計算機(jī)硬件能力的提升，CFD數(shù)值模擬已成為解決存在壁面催化等高溫?zé)峄瘜W(xué)非平衡效應(yīng)條件下高超聲速飛行器氣動熱環(huán)境預(yù)測的重要手段[3]。在CFD方法中，通常以催化邊界條件的形式模擬壁面催化效應(yīng)，即需給出壁面上的組分分布。完全非催化和完全催化是兩種最容易實現(xiàn)的催化邊界條件，其中完全非催化假定壁面附近組分分布不存在梯度，壁面組分分布等于流場內(nèi)點組分分布；而完全催化認(rèn)為流場中的原子、離子在壁面完全復(fù)合，在冷壁條件下，可認(rèn)為壁面組分與來流組分分布相同(超級催化)。由于完全非催化條件下壁面組分分布梯度為零，不會形成擴(kuò)散熱流，故完全非催化壁條件下的熱流計算結(jié)果將低于完全催化壁[4]。但實際上，壁面一般處于有限催化狀態(tài)，上述兩種極限情況只能用于定性分析。因此要實現(xiàn)高超聲速飛行器氣動熱環(huán)境精確預(yù)測，必須采用有限催化邊界條件。

根據(jù)催化反應(yīng)速率常數(shù)的計算方法，有限催化邊界條件可分為兩類[5]：一是指定催化效率(催化復(fù)合系數(shù))，再根據(jù)催化效率計算催化反應(yīng)速率常數(shù)[6]；二是采用壁面有限速率化學(xué)反應(yīng)動力學(xué)方法，對氣固復(fù)相催化反應(yīng)過程進(jìn)行分析和建模，直接計算催化反應(yīng)速率常數(shù)[7]，得到催化反應(yīng)速率常數(shù)后，再根據(jù)表面質(zhì)量守恒原理確定壁面各組分分布，即實現(xiàn)了有限催化邊界條件。由于第一種實現(xiàn)方法相對簡單，且有大量防熱材料催化系數(shù)實驗數(shù)據(jù)的支持[8-10]，因此得到了廣泛的應(yīng)用。

返回艙等飛行器再入過程中將因電離反應(yīng)形成等離子體鞘套，嚴(yán)重時還會導(dǎo)致通訊黑障，因此壁面處也會發(fā)生離子組分的催化反應(yīng)。而目前已發(fā)表的研究工作主要是針對氧原子、氮原子以及碳原子等中性原子的有限催化反應(yīng)[11-13]，存在離子組分時壁面有限催化對氣動加熱的影響規(guī)律還不清晰。為此本文首先采用指定催化效率的方法，發(fā)展了包含離子組分的有限催化邊界條件，并結(jié)合多組分化學(xué)非平衡N-S方程數(shù)值求解，建立了有限催化條件下的高超聲速飛行器氣動熱環(huán)境計算方法。然后針對類聯(lián)盟號飛船返回艙外形，采用7組分電離空氣化學(xué)模型和有限催化邊界條件對其典型再入工況氣動熱環(huán)境開展了計算分析，研究了存在離子組分時其氣動熱環(huán)境隨壁面催化效率的變化規(guī)律，并對壁面有限催化影響氣動加熱的物理機(jī)制進(jìn)行了探討。

1 計算方法

1.1 流動控制方程

本文旨在考察分子離解、電離以及催化復(fù)合等現(xiàn)象對氣動熱環(huán)境的影響，因此只考慮了化學(xué)非平衡效應(yīng)，暫時不模擬熱力學(xué)非平衡和輻射等現(xiàn)象。流動控制方程為三維守恒形式的多組分化學(xué)非平衡氣體N-S方程組[14]：

式中，Q為守恒變量，Sr為化學(xué)反應(yīng)源項，F(xiàn)、G、H為無黏通量，F(xiàn)v、Gv、Hv為黏性通量。采用采用有限體積法，在網(wǎng)格控制體單元內(nèi)對控制方程進(jìn)行積分，結(jié)合Gauss定理，得：

式中，Σ和n分別為控制體單元表面面積和外法向單位向量，V為控制體單元體積。f=Fi+Gj+Hk和fv=Fvi+Gvj+Hvk分別是無黏和黏性通量矢量。

黏性應(yīng)力張量、熱流矢量、組分質(zhì)量擴(kuò)散通量分別滿足牛頓應(yīng)力關(guān)系、傅里葉熱傳導(dǎo)定律和費克擴(kuò)散定律這些線性輸運本構(gòu)關(guān)系。對飛行器的氣動加熱由下式計算：

式中，Qc、Qd分別為熱傳導(dǎo)和組分?jǐn)U散對總氣動加熱的貢獻(xiàn)，可稱為傳導(dǎo)熱流和擴(kuò)散熱流。

本文采用了7組分Dunn-Kang[15]電離空氣模型,組分比焓hs等熱力學(xué)函數(shù)由多項式擬合給出[16]，黏性系數(shù)、熱傳導(dǎo)系數(shù)及擴(kuò)散系數(shù)等輸運系數(shù)的具體計算方法可參考文獻(xiàn)[17]。

1.2 有限催化邊界條件

在CFD數(shù)值模擬中，壁面催化效應(yīng)通常以催化邊界條件的形式出現(xiàn)，即需給出壁面上的組分分布。下面以7組分電離空氣為例，給出包含離子組分的有限催化邊界條件。

高溫空氣包含離子組分時，壁面催化反應(yīng)可包括原子復(fù)合和離子復(fù)合兩種類型。對7組分電離空氣，考慮如下3個催化復(fù)合反應(yīng)：

由催化反應(yīng)引起的單位時間，單位面積壁面上組分質(zhì)量消耗為：

式中，ks為催化反應(yīng)速率常數(shù)，ρw為壁面處密度，fws為壁面處組分質(zhì)量分?jǐn)?shù)，負(fù)號表示質(zhì)量損失。

給定催化效率后，壁面無滑移時ks可由下式計算[18]：

式中γs為催化效率，也被稱為催化復(fù)合系數(shù)。它等于在壁面發(fā)生催化復(fù)合的原子數(shù)與入射到壁面總原子數(shù)之比。顯然其取值在0到1之間，取0時表示完全非催化，取1時表示完全催化。

中性分子為催化反應(yīng)的生成物，根據(jù)質(zhì)量守恒原理，催化反應(yīng)引起的壁面單位面積組分質(zhì)量增加為：

式中，MNO、MNO+為對應(yīng)組分的摩爾質(zhì)量。另一方面，壁面處由擴(kuò)散作用產(chǎn)生的質(zhì)量通量為：

其中n為壁面法向單位矢量，方向由壁面指向流場內(nèi)部。壁面無質(zhì)量引射時，各組分單位時間單位面積由催化反應(yīng)產(chǎn)生(消耗)的質(zhì)量與擴(kuò)散出(到)飛行器表面的質(zhì)量應(yīng)相等，即各組分在壁面處的凈質(zhì)量通量為零。故下述關(guān)系式成立：

由于邊界條件中涉及組分質(zhì)量分?jǐn)?shù)的壁面法向?qū)?shù)，實際計算時還需要恰當(dāng)離散。本文綜合計算效率和計算穩(wěn)定性兩方面考慮，發(fā)展了一種近似隱式處理方法：壁面組分質(zhì)量分?jǐn)?shù)采用隱式離散，其它系數(shù)項采用顯示離散。例如，s=N,O,NO+時，式(9)可離散為：

式中，上標(biāo)n表示用第n時間步的量計算，dn為壁面第一層網(wǎng)格中心到壁面的距離，f1s為該處組分質(zhì)量分?jǐn)?shù)。流場迭代計算時，已知n時間步流場變量后，根據(jù)下式更新壁面組分質(zhì)量分?jǐn)?shù)：

其他組分處理方法類似，不再贅述。在給定催化效率后，由前述關(guān)系式可以確定所有重粒子組分在壁面的分布，最后再根據(jù)電中性假設(shè)得出電子在壁面的質(zhì)量分?jǐn)?shù)。

1.3 求解方法

對控制方程進(jìn)行無量綱化處理后，采用有限體積法進(jìn)行離散。為了克服剛性問題，同時兼顧計算穩(wěn)定性和計算效率，對無黏項和化學(xué)反應(yīng)源項采用隱式格式離散，而黏性項的離散采用顯式二階中心格式離散。最后采用LU-SGS方法[19]進(jìn)行時間推進(jìn)求解，并使用局部時間步長加速收斂。為解決復(fù)雜外形網(wǎng)格生成困難問題，采用了多塊結(jié)構(gòu)網(wǎng)格MPI分區(qū)并行求解技術(shù)。

氣動熱計算好壞關(guān)鍵在于能否準(zhǔn)確模擬邊界層，特別是壁面附近的溫度梯度、組分質(zhì)量分?jǐn)?shù)梯度。高超聲速流場中含有強(qiáng)激波等間斷，如果計算格式耗散太小可能引起計算發(fā)散或出現(xiàn)紅玉現(xiàn)象等非物理解，而計算格式耗散太大又會影響邊界層的刻畫，最終都不能得到好的熱流計算效果。因此要想準(zhǔn)確計算熱流，需要在激波處引入適當(dāng)耗散以穩(wěn)定激波，同時在邊界層內(nèi)要盡量減小耗散。為解決這一問題，本文在流場變量梯度大的方向(壁面法向)采用耗散小的Godunov格式[20]，而在流場變量梯度小的方向采用耗散大的Steger-Warming矢通量分裂格式[21]，以提高計算魯棒性，同時不降低熱流計算精度[22]。

對于邊界條件的處理，壁面可采用無滑移完全非催化/催化，以及有限催化邊界條件，壁面溫度可按等溫壁設(shè)定或由輻射平衡關(guān)系計算。由于非規(guī)則外形流場存在激波/激波,激波/邊界層干擾現(xiàn)象，通常使用的壁面法向零壓力梯度條件不一定適用，故本文采用由內(nèi)點外插的方法計算壁面壓力。超聲速來流邊界則直接采用來流條件，對于超聲速出口邊界采用一階外插的方法。

2 方法驗證

標(biāo)模ELECTRE幾何外形是半錐角為4.6°的球錐體，頭部半徑為0.175 m，彈身總長為2.0 m，且擁有飛行試驗測熱數(shù)據(jù)[23]。其中發(fā)射后第293 s時刻(高度53.3 km，迎角近似為0°，馬赫數(shù)13)公開報道的熱流數(shù)據(jù)最完整，包含沿標(biāo)模母線所有測點數(shù)據(jù)；且該工況飛行速度在4000 m/s以上，高溫氣體效應(yīng)已不能忽略；加之飛行高度在50 km以上，可以避免轉(zhuǎn)捩及湍流模擬，更易于考察高溫氣體效應(yīng)。基于上述原因，已有大量文獻(xiàn)把ELECTRE第293 s時刻工況作為高溫非平衡流場及氣動熱計算方法的考核算例。本文為驗證所建方法和計算程序的可靠性，也對該算例進(jìn)行了考核計算。計算時壁面溫度取為343 K,且采用了完全非催化和完全催化兩種壁面催化條件。

網(wǎng)格分布，特別是壁面法向網(wǎng)格間距對氣動熱計算影響較大[24]。為考核計算方法對氣動熱計算的網(wǎng)格無關(guān)性，并由此確定合適的壁面法向網(wǎng)格間距，共劃分了Rec=1、5、10、20、50五種不同網(wǎng)格雷諾數(shù)的計算網(wǎng)格。網(wǎng)格雷諾數(shù)的定義如下：

式中，ρ∞、V∞、μ∞分別為自由來流密度、速度和黏性系數(shù)，Δn為壁面法向第一層網(wǎng)格間距。為同時驗證程序的三維問題計算能力，五套網(wǎng)格均為三維網(wǎng)格，且拓?fù)浣Y(jié)構(gòu)一致，壁面法向網(wǎng)格數(shù)均為81，如圖1所示。

圖2給出了不同網(wǎng)格和壁面催化條件下ELECTRE彈身母線熱流密度計算結(jié)果與飛行試驗數(shù)據(jù)比較情況。圖例中“nc”表示壁面為完全非催化，“fc”表示壁面為完全催化。本文的計算方法具有較好的網(wǎng)格無關(guān)性，網(wǎng)格雷諾數(shù)小于20后熱流計算結(jié)果基本相同。熱流分布規(guī)律計算與試驗基本一致，除尾部外飛行試驗數(shù)據(jù)處于完全非催化壁和完全催化壁預(yù)測結(jié)果之間。

圖2 計算與飛行試驗結(jié)果比較Fig.2 Comparison among the CFD results and Flight data

3 返回艙熱環(huán)境計算分析

利用上節(jié)所述計算方法，對類聯(lián)盟號飛船返回艙氣動熱環(huán)境開展了數(shù)值計算。由于返回艙通常采用碳基或硅基防熱材料，其催化效率通常小于0.1，因此本文只考慮了0.0、0.02、0.05、0.1四種壁面催化效率，且假設(shè)各組分的催化效率相同。計算工況高度為70 km，馬赫數(shù)為24.9，迎角為30°，壁面溫度為350 K。由于該工況飛行高度較高，均按層流流態(tài)計算。

圖3是返回艙流場溫度分布云圖。返回艙激波壓縮強(qiáng)烈，過激波后流場溫度急劇升高，并在返回艙大底、迎風(fēng)面與弓形激波之間形成了高溫區(qū)。與量熱完全氣體總靜關(guān)系式估計結(jié)果比較，波后溫度明顯降低，可以預(yù)見波后高溫激波層內(nèi)發(fā)生離解等吸熱化學(xué)反應(yīng)。這一點可從圖4給出的流場組分質(zhì)量分?jǐn)?shù)云圖得到證實：過激波后O2分子和N2分子發(fā)生離解并生成相應(yīng)O原子與N原子，因此原子質(zhì)量分?jǐn)?shù)明顯增加而對應(yīng)分子質(zhì)量分?jǐn)?shù)顯著下降。 此外流場中還生成了NO分子，并進(jìn)一步電離產(chǎn)生了NO+離子，因此有必要研究包含離子組分時有限催化對返回艙氣動加熱的影響。

圖5給出了不同催化效率條件下返回艙表面熱流密度分布情況。當(dāng)催化效率等于零，即完全非催化壁時熱流密度最低；隨著催化效率增大，返回艙迎風(fēng)面熱流密度明顯增加?？梢姳诿嬗邢薮呋瘜鈩蛹訜嵊兄匾绊憽２捎玫痛呋实姆罒岵牧峡捎行Ь徍头祷嘏摎鈩訜岘h(huán)境。

駐點是返回艙熱環(huán)境最為嚴(yán)酷的部位，圖6給出了返回艙駐點峰值熱流隨壁面催化效率的變化情況。

圖3 返回艙流場溫度云圖Fig.3 Temperature contour of capsule flowfield

圖4 組分質(zhì)量分?jǐn)?shù)云圖Fig.4 Mass fraction contours

圖5 返回艙表面熱流密度云圖Fig.5 Heat-flux contour of capsule surface

如圖6所示，在本文考慮的范圍內(nèi)，返回艙駐點總熱流和擴(kuò)散熱流均隨著催化效率單調(diào)遞增，但熱流增加率并不是線性，而是隨催化效率逐漸減小。擴(kuò)散熱流對催化效率更加敏感，且量值上可以超過傳導(dǎo)熱流。此外通過比較總熱流與擴(kuò)散熱流的變化規(guī)律可以看出，在駐點區(qū)域，催化效應(yīng)主要通過擴(kuò)散機(jī)制影響氣動加熱。

圖6 返回艙駐點熱流密度比較Fig.6 Heat-flux contour of capsule surface

圖7 返回艙大底中心線熱流密度分布Fig.7 Heat-flux distribution along the center line of capsule ’s base

除駐點外，返回艙大底也是氣動加熱嚴(yán)重的部位，為分析有限催化對其熱環(huán)境的影響，圖7給出了不同催化效率條件下沿大底中心線的熱流密度分布情況?？梢钥闯鲇邢薮呋瘲l件下熱流密度明顯升高，但中心線不同位置處熱流密度隨催化效率的變化規(guī)律不盡相同：駐點附近區(qū)域熱環(huán)境受催化作用影響顯著，熱流密度隨催化效率的絕對增幅也最大；而在遠(yuǎn)離駐點的流動膨脹區(qū)，催化效率大于0.02后熱流密度變化不大。

為進(jìn)一步厘清壁面有限催化對氣動熱環(huán)境的影響機(jī)制，圖8和圖9分別給出了沿返回艙大底中心線的傳導(dǎo)熱流和擴(kuò)散熱流密度分布。壁面催化效率大于0.02后，傳導(dǎo)熱流基本不受壁面催化效率影響，而且量值上反而低于完全非催化壁。而擴(kuò)散熱流受壁面催化系數(shù)影響明顯，這再次說明有限催化條件下，擴(kuò)散熱流是影響氣動加熱的主要機(jī)制。

從物理上分析，催化效率反映了原子、離子在壁面發(fā)生催化復(fù)合反應(yīng)的強(qiáng)弱，對壁面及壁面附近組分質(zhì)量分?jǐn)?shù)分布有直接影響。由擴(kuò)散熱流表達(dá)式(3)可以看出，擴(kuò)散熱流與壁面法向組分質(zhì)量分?jǐn)?shù)梯度直接相關(guān)，因此對催化效率更加敏感。另一方面催化復(fù)合反應(yīng)屬于放熱反應(yīng)，對壁面附近溫度梯度也會產(chǎn)生影響。但從本文計算結(jié)果看，這種影響相對較小，與壁面溫度梯度直接相關(guān)的傳導(dǎo)熱流對催化效率并不敏感。這兩方面原因共同造成了傳導(dǎo)與擴(kuò)散熱流的相對大小會隨催化效率發(fā)生變化。同時也說明壁面催化效應(yīng)主要是通過原子、離子在壁面發(fā)生催化復(fù)合反應(yīng)，改變壁面及壁面附近組分分布，進(jìn)而導(dǎo)致擴(kuò)散熱流發(fā)生變化，對總氣動加熱產(chǎn)生影響。

圖8 返回艙大底中心線傳導(dǎo)熱流密度分布Fig.8 Conductive heat-flux distribution along the center line of capsule’s base

圖9 返回艙大底中心線擴(kuò)散熱流密度分布Fig.9 Diffusive heat-flux distribution along the center line of capsule’s base

需要指出的是，在駐點附近區(qū)域和遠(yuǎn)離駐點區(qū)域，催化效率對擴(kuò)散熱流的影響規(guī)律明顯不同。駐點附近擴(kuò)散熱流隨催化效率增大而增大，而在遠(yuǎn)離駐點的流動膨脹區(qū)擴(kuò)散熱流則可能隨催化系數(shù)增大反而減小。這說明壁面催化對氣動加熱的影響不僅與表征材料催化能力的催化效率有關(guān)。

考慮到實際計算中，壁面法向組分梯度是以離散形式計算，即：

通過上述分析，我們可以得到如下初步推論：壁面有限催化對氣動熱的影響不僅與壁面材料催化效率有關(guān)，也與流場離解電離程度、壁面密度、溫度等當(dāng)?shù)亓鲃訁?shù)相關(guān)。這也解釋了本文返回艙算例中，在駐點附近區(qū)域和遠(yuǎn)離駐點區(qū)域，催化效率對熱流的影響規(guī)律明顯不同的原因。

4 結(jié) 論

初步建立了包含離子組分有限催化條件下的高超聲速飛行器氣動熱環(huán)境計算方法和軟件。針對類聯(lián)盟號飛船返回艙外形，研究了壁面有限催化對其氣動熱環(huán)境的影響規(guī)律，主要結(jié)論如下：

1) 表征材料催化能力的催化效率對氣動加熱影響顯著，隨著催化效率增大，返回艙迎風(fēng)面熱流密度明顯增加。采用低催化效率壁面材料可有效緩和返回艙氣動熱環(huán)境。

2) 催化作用主要通過擴(kuò)散機(jī)制影響氣動加熱。擴(kuò)散熱流對壁面催化效率更加敏感，且量值上可以超過傳導(dǎo)熱流，但熱流并不隨催化效率增加而線性增大。

3) 催化作用對氣動加熱的影響不僅與壁面材料催化效率本身有關(guān)，也與流場離解電離程度、壁面密度、溫度等當(dāng)?shù)亓鲃訁?shù)相關(guān)。

本文的工作假定了所有催化反應(yīng)的催化效率都為相同常數(shù)，且只考慮了0.1以下的催化效率，因此所得結(jié)論還需要進(jìn)一步研究確認(rèn)。本文采用了等溫壁條件，沒有研究壁溫對催化的影響。而實際上催化效率與壁面溫度是密切相關(guān)的，在以后的研究中應(yīng)考慮隨壁溫變化的催化效率模型。另外對燒蝕熱防護(hù)系統(tǒng)，壁面氧化、氮化反應(yīng)也對氣動熱環(huán)境有重要影響，因此除壁面催化反應(yīng)外，還需要發(fā)展同時考慮其它壁面反應(yīng)的邊界條件。

熱力學(xué)非平衡現(xiàn)象如振動-化學(xué)反應(yīng)耦合，會改變相關(guān)反應(yīng)的化學(xué)反應(yīng)速率，并由此對流場組分分布造成影響。而壁面催化主要是通過組分?jǐn)U散機(jī)制影響氣動加熱，對流場組分分布相對敏感，因此有必要評估熱力學(xué)非平衡效應(yīng)對壁面催化氣動加熱的影響。本文作者針對此問題已開展了初步研究，結(jié)果表明：對返回艙大底這類大鈍頭體外形，飛行速度在第一宇宙速度以下時，熱力學(xué)非平衡效應(yīng)對其氣動熱環(huán)境影響較小，不會改變本文的研究結(jié)論，相關(guān)情況將另文詳細(xì)探討。

致謝:感謝中國空氣動力研究與發(fā)展中心超高速所李志輝研究員對本文工作的幫助。