輻射平衡壁溫下有限催化模型數值模擬

莫 凡, 王鎖柱, 高振勛

(1. 北京航空航天大學航空科學與工程學院, 北京 100191;2. 北京航天長征飛行器研究所, 北京 100076)

引 言

高超聲速飛行器在大氣中進行高Mach數飛行時, 來流經過頭部脫體激波的劇烈壓縮后往往能達到很高的溫度[1], 導致高溫空氣中發(fā)生復雜的化學反應. 一般情況下, 在2 500 K以上空氣中的氧氣開始離解, 4 000 K以上氮氣開始離解, 而到9 000 K以上氧原子和氮原子將進一步電離.

當離解之后的高溫空氣到達飛行器表面時, 表面材料會對空氣中的原子和離子產生顯著的催化復合作用, 這種催化復合作用將改變飛行器壁面附近各組分的質量分數梯度分布, 引起流場中輸運特性的改變, 進而顯著改變壁面熱流密度的分布. 已有研究表明, 在高Mach數情況下關鍵部位防熱材料表面的催化復合效應最高將帶來近50%以上的氣動熱載荷[2]. 因此在CFD模擬中, 催化邊界條件的確定往往非常重要, 完全非催化(non-catalytic wall, NCW)和完全催化(full-catalytic wall, FCW)[3]是高超聲速飛行器氣動熱環(huán)境數值模擬中最常用的兩種催化邊界條件. 完全非催化假設壁面催化復合速率為0, 即壁面組分質量分數的梯度為0. 而完全催化一般假設壁面來流原子的離子在壁面完全復合, 此時壁面組分質量分數的梯度也將最大.

然而, 以上兩種組分邊界條件均不能反映實際發(fā)生的壁面催化過程, 實際飛行中均為有限催化情況, 因此必須使用有限催化[3](finite-rate catalysis, FRC)邊界條件才能反映真實的壁面催化效應, 該邊界條件的關鍵在于催化復合系數的取值. 粟斯堯等[4]選取了催化復合系數為常數的有限催化邊界條件, 針對類聯盟號飛船返回艙外形研究了氣動熱環(huán)境隨壁面催化效率的變化規(guī)律, 并對壁面有限催化影響氣動加熱的物理機制進行了探討. 丁明松等[5]開展了在有限催化條件下, 表面材料催化特性差異性對熱環(huán)境的影響, 得出壁面催化復合系數的差異性會帶來局部熱流的跳變, 使局部區(qū)域熱流明顯高于全表面催化復合系數為0.1的情況, 且壁面復合系數差異越大, 熱流跳變的程度越大. 然而對于實際的高超聲速飛行器表面催化效應, 其表面的復合系數取常數是不恰當的. Inger[6]的實驗表明催化復合系數往往與壁面材料和壁面溫度有關, 并給出了某碳基材料催化復合系數隨表面溫度變化的擬合公式. Scott[7]開展了典型材料的催化復合系數研究, 并給出了2 000 K以下大部分碳基和硅基熱防護材料催化復合系數小于0.1的結論. Scott[8]和Zoby等[9]分別給出了實驗和飛行數據修正后氧原子和氮原子催化復合系數隨溫度擬合的公式. 國內關于采用溫度擬合催化復合系數研究有限催化的工作尚少, 董維中等[10]采用源于實驗數據擬合得到的與壁溫相關的O原子的催化復合系數, 開展了輻射平衡表面溫度分布與氣動熱的耦合計算, 但由于未考慮氮原子的催化復合效應, 因此模型并不完整.

本文在現有考慮離子組分的有限催化邊界條件的基礎上, 引入了由溫度擬合得到的氮氧原子催化復合系數, 研究了在輻射平衡壁溫下有限催化模型對典型再入高超聲速飛行器氣動熱環(huán)境的影響規(guī)律. 文中首先給出了流動控制方程及有限催化邊界條件的數學形式, 在此基礎上建立了高超聲速飛行器化學非平衡流場有限催化氣動熱環(huán)境數值計算方法及程序, 并通過與風洞實驗氣動熱結果對比驗證了數值模擬結果的可靠性. 在此基礎上考慮到有限催化中催化復合系數與壁面溫度相關, 采用了Scott和Zoby給出的氧原子和氮原子催化復合系數隨溫度擬合的公式, 并研究鈍錐外形在輻射平衡壁溫下有限催化對壁面熱環(huán)境的影響規(guī)律.

1 計算方法

1.1 化學非平衡流動的瞬時N-S方程

在Descartes坐標系下, 含化學反應流動的瞬時N-S方程可寫成如下守恒形式

(1)

式中,ρ,ui,p,T分別為氣體的密度、 速度、 壓強與靜溫;E為單位(體積)質量總能量

(2)

(3)

其中，μ為混合氣體的分子黏性系數. 各組分分子黏性系數、 熱傳導系數、 質量擴散系數等輸運系數和比熱、 比焓均采用Gupta的擬合公式得到[11], 而空氣混合物整體的輸運系數則根據Wilke的公式計算得到[12].

1.2 數值求解方法

在對方程(1)無量綱化處理后, 采用有限差分方法進行離散, 其中對流項采用Roe格式離散, 并采用Harten-Yee型熵修正, 且熵修正系數取0.25, 黏性項采用2階中心差分離散. 時間推進上為提高CFL步長以及求解的穩(wěn)定性, 采用LU-SGS隱式方法, 并采用當地時間步長加速收斂. 對于三維復雜多塊結構網格, 采用MPI分區(qū)并行求解技術以提高計算效率.

1.3 有限催化邊界條件

在高超聲速化學非平衡流動的CFD模擬中, 壁面催化效應通過壁面邊界條件影響對應的控制方程的解. 下面以7組分電離空氣為例, 給出了包含離子在內的組分質量分數邊界條件.

對于包含離子組分的高溫氣體, 在壁面上發(fā)生的催化復合反應主要有兩類, 即原子的復合和離子的復合. 現僅考慮以下3個催化復合反應

N+N→N2; O+O→O2; NO++e-→NO

(4)

有限催化壁邊界條件一般通過求解壁面附近質量守恒方程的方法構建. 不考慮質量引射時, 氣體組分由催化反應生成(或消耗)的質量與擴散出(或入)壁面微元的質量相等.

由催化反應引起的單位時間, 單位面積壁面上各組分質量消耗(生成)率為

(5)

其中，kws為反應速率常數. 假設壁面原子服從Maxwell-Boltzmann分布, 則kws與壁面催化復合系數γs具有如下關系[13]

(6)

式中,R0為普適氣體常數;Tw為壁面溫度;Ms為組分s的摩爾質量. 壁面催化復合系數γs表示壁面發(fā)生催化復合的原子數與入射到壁面總原子數之比. 其取值范圍在0到1之間, 取0時表示完全非催化(NCW), 取1時表示完全催化(FCW), 取0到1之間則為有限催化(FRC).

另一方面, 由Fick定律得壁面處由擴散作用產生的質量通量為

(7)

其中，n為壁面法向單位矢量, 方向由壁面指向流場內部. 根據質量守恒定律, 各組分單位時間單位面積由催化反應產生(消耗)的質量與擴散出(到)壁面的質量應相等, 即各組分在壁面處的凈質量通量為零. 故下述關系式成立

(8)

式(7)中存在導數項, 在CFD計算時需對其進行恰當的離散, 本文中消耗組分(s=O,N,NO+)采用以下離散方式

(9)

整理得

(10)

(11)

整理得

(12)

1.4 輻射平衡壁溫邊界條件

輻射平衡模型依據的主要原理是在壁面處能量守恒. 飛行器壁面一方面要接收流場的熱流(壁面熱流密度)Qflux, 另一方面要向內部結構傳熱Qinside, 另外還要通過電磁輻射向外發(fā)射能量(灰體輻射能量)Qradiation.當流場處于平衡狀態(tài)時, 上述3種能量的傳輸也應滿足以下的平衡方程

Qflux=Qinside+Qradiation

(13)

現對結構的內部傳熱忽略不計, 平衡方程則可以進一步化簡為

Qflux=Qradiation

(14)

即流場向壁面的傳熱量應與壁面出射的輻射能量達到平衡. 當流場氣體考慮組分時, 壁面熱流密度可以表示為對流項熱流密度Qconv和組分擴散項熱流密度Qdiff之和,可由式(15)進行計算

(15)

其中,n代表壁面的法向方向,ρ為壁面流場密度,ns為氣體組分個數,Ds代表組分s的擴散系數,hs是組分s的絕對焓,Ys代表組分s的質量分數.

若將壁面視為發(fā)射率等于ε的灰體, 則壁面出射的輻射能根據Stefan-Boltzmann定律, 可以寫作

(16)

其中,σ為Stefan-Boltzmann常數, 取值為σ=5.67×10-8W/(m2·K4);Tw為壁面溫度. 聯立式(13), (15), 并將Qflux離散化, 則最終可得平衡方程

(17)

其中,T1和Tw分別為第1層網格和壁面處的溫度;Ys1和Ysw分別為第1層網格和壁面處組分s的質量分數; Δy為第1層網格與壁面間的距離. 求解式(17)采用了Newton迭代法, 迭代關系如式(18)所示

(18)

2 等溫壁有限催化計算

卡爾斯本大學巴法羅研究中心(CUBRC)利用LENS系列激波風洞對標模返回艙進行了一系列的研究[14]. 其中幾何尺寸如圖1所示, 實驗工況如表1所示. 根據實驗來流靜溫和Mach數可計算得到來流總溫超過8 000 K, 因此采用了Gupta 5組分6反應化學反應動力學模型, 壁面上僅考慮O原子和N原子的催化復合反應. 計算網格如圖2所示. 為了更好地捕捉頭部脫體激波, 在返回艙頭部前體附近對網格加密至0.2 mm.

圖1 返回艙模型幾何尺寸(單位： in[mm])Fig. 1 Geometric dimensions of re-entry capsule model(unit: in[mm])

表1 返回艙實驗工況

LENS激波風洞實驗給出了標模返回艙對稱面壁面上測量點的壓力與熱流測量值. 針對Test-1的CFD計算中組分邊界條件采用完全催化(FCW)、 完全非催化(NCW)以及γ分別取0.015和0.1的有限催化(FRC), 在圖3中與實驗熱流進行了對比. 由于催化效應幾乎沒有改變壁面的壓強分布, 因此壓強曲線在圖2中僅給出一條并與實驗結果對比.

針對 Test-2的CFD 計算中組分邊界條件采用完全催化(FCW)、 完全非催化(NCW)以及γ取0.005 的有限催化(FRC), 在圖4中與實驗壁面熱流及壓強進行了對比. 圖5給出了壁面的壓強分布云圖. 圖6～8分別給出了完全非催化, 有限催化和完全催化邊界條件下的熱流密度分布云圖.

從Test-1與Test-2的CFD計算結果可以得到壁面熱流密度分布QFCW>QFRC>QNCW.這與以往經驗和其他文獻結論相符.

圖2 返回艙模型計算網格Fig. 2 Three-dimensional structural grids for capsule geometry

圖3 Test-1返回艙壁面熱流與壓強分布Fig. 3 Heat flux and pressure distribution on the capsule wall in Test-1

圖4 Test-2返回艙壁面熱流與壓強分布Fig. 4 Heat flux and pressure distribution on the capsule wall in Test-2

圖5 Test-2返回艙壁面壓強分布云圖Fig. 5 Contour of pressure on the capsule wall in Test-2

圖6 Test-2返回艙壁面完全非催化熱流密度分布云圖Fig. 6 Contour of heat flux at the NCW in Test-2

圖7 Test-2返回艙壁面有限催化熱流密度分布云圖Fig. 7 Contour of heat flux at the FRC in Test-2

圖8 Test-2返回艙壁面完全催化熱流密度分布云圖Fig. 8 Contour of heat flux at the FCW in Test-2

而在使用有限催化模型時, 當催化復合系數γ取合適值時可以得到與實驗結果更加吻合的結果. 這從側面反應了高超聲速化學非平衡流動的壁面催化反應是有限催化, 與以往人們的觀點相符. 此外, 當催化復合系數γ取0時, 組分質量分數邊界條件自動退化為完全非催化壁條件.

由圖9分析Test-1壁面熱流組成可得, 壁面熱流由熱傳導引起的熱流和各組分質量擴散引起的熱流組成, 且在本算例中熱傳導引起的熱流大于各組分質量擴散引起的熱流. 在質量擴散引起的熱流中, 氧原子質量擴散引起的熱流密度占主要成分, 其他組分(如氮原子)僅僅占很小一部分, 這與各組分質量分數在壁面的梯度有關. 以氧原子和氮原子為例, 如圖10所示, 由于激波后靜溫不到7 000 K, 氮氣分子離解程度遠小于氧氣分子, 因此在壁面處氮原子質量分數梯度遠小于氧原子質量分數梯度, 進而由氮原子質量擴散引起的熱流密度小于氧原子質量擴散引起的熱流密度.

圖9 Test-1返回艙壁面熱流組成Fig. 9 Composition of the heat flux in Test-1

圖10 Test-1有限催化駐點線氧原子和氮原子質量分數分布圖Fig. 10 Mass fraction distribution of oxygen and nitrogrn along the stagnation line at the FRC in Test-1

圖11～12給出了Test-1有限催化氧分子沿駐點線質量分數分布, 可得壁面氧原子質量分數下降, 氧氣上升, 表明在催化邊界的作用下, 氧原子在壁面復合為氧氣. 完全催化的復合程度最大且壁面氧原子質量分數降到接近于0, 小于來流的氧原子質量分數.

完全非催化條件下, 氧原子質量分數在壁面也有一定程度的降低. 這是由于壁面溫度僅有300 K, 激波后到壁面之間溫度逐漸下降, 在壁面附近溫度急劇下降, 而溫度影響了化學反應速率進而使當地流動趨向于新的化學平衡, 因此由激波后高溫離解產生的原子在壁面附近自動復合.

由以上可以更清晰地得出, 壁面復合效應有兩個組成部分, 第一是由壁面材料催化直接導致的復合, 第二是由壁面附近低溫導致的化學反應速率變化而間接引起的復合.

圖11 Test-1有限催化氧原子沿駐點線質量分數Fig. 11 Mass fraction distribution of oxygen along the stagnation line at different catalytic conditions in Test-1

圖12 Test-1有限催化氧分子沿駐點線質量分數Fig. 12 Mass fraction distribution of oxygen molecule along the stagnation line at different catalytic conditions in Test-1

3 輻射平衡模型驗證

北卡羅萊納州立大學的 Keenan 等對高超聲速球繞流問題進行了數值模擬研究[15], 且在文中考慮了壁面輻射的影響. 本節(jié)擬參考該數值模擬的結果, 對本文的輻射平衡模型進行驗證. 計算的具體參數如下: 球的半徑R=1 m, 壁面的發(fā)射率ε=0.9, 海拔高度為H=65 km, 來流速度取 8 km/s.

圖13給出了高超聲速球頭繞流計算所采用的結構網格, 并在激波附近對網格進行了加密, 以便提高對激波的捕捉能力. 圖14給出了高超聲速球頭繞流壓強的分布云圖, 可以發(fā)現來流在經過頭部脫體激波后壓強迅速增加, 其中在駐點附近由于激波強度最大, 因此壓強增加也最為顯著.

圖13 高超聲速球繞流網格Fig. 13 Structural grids for the hypersonic spherical body

圖14 高超聲速球繞流壓強分布云圖Fig. 14 Contour of pressure for the hypersonic spherical body

圖15, 16分別給出了球頭沿流向壁面上的溫度和熱流密度分布. 可以看出, 壁面溫度與文獻結果十分吻合, 熱流密度與文獻結果在中間部分有一定偏差, 但最大偏差不超過10%. 由此驗證了輻射平衡模型在計算壁面溫度分布時的正確性, 所以可以將此模型應用于之后的飛行器外形溫度分布計算.

圖15 溫度分布, v=8 km/sFig. 15 Distribution of temperature along the wall， v=8 km/s

圖16 熱流密度分布, v=8 km/sFig. 16 Distribution of heat flux along the wall， v=8 km/s

4 輻射平衡壁溫有限催化計算

本節(jié)借助三維高超聲速鈍頭雙錐繞流實驗模型[16], 將該模型等比例放大, 放大之后的模型外形參數如下: 鈍錐體頭部的前緣半徑為 0.153 4 m, 前錐體的半錐角為12.84°, 后錐體的半錐角為7°, 兩錐的連接處與頭部的距離為2.782 m, 全錐總長4.889 6 m. 計算采用的結構網格如圖17所示， 其中對雙錐頭部網格進行了加密處理以便更精確地捕捉頭部脫體激波. 計算工況采用STS-2飛行實驗來流參數, 如表2所示.

圖17 雙錐結構網格Fig. 17 Structural grids for bi-cone geometry

表2 雙錐計算工況

由于Mach數等于26.3, 靜溫202.72 K, 總溫達到了15 000 K以上, 因此需要考慮電離反應, 采用Gupta 7組分9反應化學反應動力學模型, 壁面催化效應考慮式(4)中的3個反應. CFD計算分別采用完全催化、 有限催化和完全非催化邊界條件. 其中， 有限催化模型的氧原子和氮原子催化復合系數采用Scott實驗數據擬合結果和Zoby飛行數據修正結果. 考慮到催化復合系數隨溫度變化的連續(xù)性, 因此本節(jié)最終采用式(19)的連續(xù)形式, 且γNO+取常數.

(19)

圖18給出了雙錐輻射平衡壁溫下有限催化流場溫度云圖. 圖19～21分別給出了完全催化、 有限催化和完全非催化條件下駐點線各組分質量分數分布, 其中NO+和e-由于電離程度較低, 組分質量分數含量極少, 最終這兩種組分在圖中幾乎與x軸重合. 對于氧氣、 氮氣質量分數, 壁面處完全催化大于有限催化大于完全非催化. NO質量分數, 有限催化大于完全催化大于完全非催化,這與壁面溫度及各組分質量分數有關.

圖18 雙錐輻射平衡壁溫下有限催化溫度云圖Fig. 18 Contour of temperature at the radiation balanced wall under FRC

圖19 完全催化駐點線組分質量分數分布Fig. 19 Mass fraction distribution of components along the stagnation line under FCW

圖20 有限催化駐點線組分質量分數分布Fig. 20 Mass fraction distribution of components along the stagnation line under FRC

圖21 完全非催化駐點線組分質量分數分布Fig. 21 Mass fraction distribution of components along the stagnation line under NCW

圖22給出了在輻射平衡邊界條件下, 完全催化, 有限催化和完全非催化駐點線溫度分布曲線. 表明在高空高Mach情況下, 來流空氣稀薄, 激波厚度較厚, 同時壁面催化復合效應也將改變激波位置. 其中完全非催化激波位置離壁面更遠, 且駐點溫度更低.

圖23, 24分別給出了完全催化、 有限催化和完全非催化壁面的溫度分布和熱流密度分布. 圖23表明有限催化下壁面溫度的峰值約為2 265 K, 而在大面積區(qū)域壁面溫度在1 150 K左右, 同時通過式(19)的催化復合系數計算公式求得在駐點附近氧原子的催化復合系數約為0.17, 氮原子的催化復合系數約為0.026, 而在大面積區(qū)域, 氧原子的催化復合系數約為0.005 3, 氮原子的催化復合系數約為 0.01. 另一方面對比了圖24中3種催化邊界的壁面熱流密度, 結果表明完全催化的壁面熱流峰值比有限催化高約21%, 而完全非催化的壁面熱流峰值比有限催化低約29%, 對于大面積的熱流分布, 完全催化的壁面熱流峰值比有限催化高約10%, 而完全非催化的壁面熱流峰值比有限催化低約8.7%.

圖22 完全催化、 有限催化和非催化駐點線溫度對比Fig. 22 Comparison of temperature along the stagnation line between FCW, FRC and NCW

圖23 完全催化、 有限催化和非催化壁面溫度對比Fig. 23 Comparison of wall temperature between FCW, FRC and NCW

圖24 完全催化、 有限催化和非催化壁面熱流對比Fig. 24 Comparison of heat flux between FCW, FRC and NCW

圖25給出了有限催化壁面總熱流、 熱傳導熱流與各組分熱流分布. 對比圖24， 25中有限催化和完全非催化的熱傳導熱流可以發(fā)現, 有限催化條件下壁面的熱傳導熱流顯著增加, 這是由于在有限催化和輻射平衡耦合作用下, 使得激波位置相較于完全非催化發(fā)生了改變, 且更靠近壁面(如圖22所示), 從而使得壁面溫度梯度上升, 同時在這種耦合作用下壁面溫度也將上升, 進而導致壁面熱傳導系數上升, 最終使得熱傳導熱流增強. 圖26給出了有限催化壁面熱傳導熱流與各組分熱流占總熱流比值. 結合圖25分析可知, 在此來流條件下, 熱流密度中熱傳導引起的熱流占主要成分, 質量擴散引起的熱流中氧原子擴散引起的熱流占主要成分, 且質量擴散引起的熱流在頭部駐點占比最大, 并隨著壁面離駐點的距離增大而下降.

圖25 有限催化壁面總熱流、 熱傳導熱流與各組分熱流對比Fig. 25 Comparison of total heat flux, thermal heat flux and heat flux of each component under FRC

圖26 有限催化壁面熱傳導熱流與各組分熱流占總熱流比值Fig. 26 Ratio of thermal heat flux and heat flux of each component to total heat flux under FRC

5 結論

本文發(fā)展了高超聲速飛行器輻射平衡壁溫下有限催化的數值方法, 基于返回艙外形的風洞實驗數據進行了數值模擬結果對比, 并進一步針對典型高超聲速飛行器鈍雙錐研究了輻射平衡壁溫下有限催化對氣動熱環(huán)境的影響規(guī)律, 主要結論如下:

(1)針對返回艙外形的數值實驗結果表明, 完全催化與完全非催化邊界條件下壁面熱流密度均與風洞實驗結果偏差較大, 不能表征真實壁面的催化復合效應, 而采用合適的有限催化模型, 則壁面熱流密度與風洞試驗結果符合良好, 表明高超聲速飛行器壁面的催化復合效應為有限催化.

(2)針對典型高超聲速飛行器鈍雙錐的研究表明, 在輻射平衡溫度邊界條件下駐點附近氧原子的催化復合系數約為0.17, 氮原子的催化復合系數約為0.026, 大面積區(qū)則分別降為0.005 3和0.01. 另一方面, 與完全催化和完全非催化相比, 有限催化對壁面溫度和熱流分布具有顯著影響. 駐點壁面溫度方面, 完全催化的壁面溫度峰值比有限催化高約5%, 而完全非催化的壁面溫度峰值比有限催化低約8%. 駐點熱流密度方面, 完全催化的壁面熱流峰值比有限催化高約21%, 而完全非催化的壁面熱流峰值比有限催化低約29%.

(3)高超聲速化學非平衡流場中飛行器壁面的組分復合效應有兩個組成部分, 第一是由壁面材料催化直接導致的復合, 第二是由壁面附近低溫導致的化學反應速率變化而間接引起的復合. 輻射平衡壁溫條件下, 與完全非催化相比, 有限催化效應將間接改變壁面溫度梯度及壁面熱傳導系數, 進而影響熱傳導熱流.

本文建立的高超聲速飛行器化學非平衡流場輻射平衡與有限催化邊界耦合計算方法具有較大的拓展空間, 未來可結合分子動力學計算得到的壁面催化復合系數, 實現跨尺度的CFD數值模擬.