聶春生,楊光,聶亮,周禹,趙良

(中國運(yùn)載火箭技術(shù)研究院空間物理重點(diǎn)實(shí)驗(yàn)室,北京 100076)

0 引言

臨近空間飛行器高馬赫數(shù)飛行時,高溫空氣發(fā)生電離,形成一個非均勻等離子體包裹流場[1-2],同時飛行器外表面的防熱材料會發(fā)生燒蝕,燒蝕產(chǎn)物通過引射和擴(kuò)散效應(yīng)進(jìn)入流場,影響等離子體流場特性,因此準(zhǔn)確獲取等離子流場參數(shù),必須考慮表面防熱材料燒蝕產(chǎn)物對等離子體流場的影響[3-4]。

國外在高溫等離子體流場方面的研究較早,針對球錐再入體開展了一系列研究,最有代表的是電磁衰減測量項目(RAM)[5],整個項目進(jìn)行了8次飛行試驗(yàn),從理論和試驗(yàn)兩方面研究了等離子鞘套對再入體通信系統(tǒng)的影響,其中第1次和第2次飛行試驗(yàn)端頭采用碳酚醛燒蝕材料,開展了表面材料燒蝕產(chǎn)物對等離子體流場的影響研究,但相關(guān)結(jié)果公開發(fā)表的很少。數(shù)值研究方面,Bhutta等[6]最先將數(shù)值方法引入燒蝕產(chǎn)物對高溫非平衡流場的影響研究;Lewis等[7]利用拋物化Narier-Stokes(N-S)方程數(shù)值研究了泰氟隆材料燒蝕產(chǎn)物對高溫流場的影響;Bianchi等[8-9]對電弧風(fēng)洞中的球錐模型石墨燒蝕流場進(jìn)行了計算分析;Keenan等[10]針對半球和鈍錐體外形,數(shù)值研究了碳碳材料燒蝕與流場的耦合的影響;Thomas等[11]建立了高超聲速平衡燒蝕再入三維模擬程序。

國內(nèi)還沒有開展專門針對高溫等離子體流場研究的飛行試驗(yàn),但在地面風(fēng)洞環(huán)境中開展了一些研究工作。袁野等[12]在高頻等離子體風(fēng)洞中研究了碳/碳和碳/炭化鋯兩種復(fù)合防熱材料燒蝕情況下對下游流場的影響;聶春生等[13]在高頻等離子體風(fēng)洞中開展了純碳-碳材料燒蝕對球錐體等離子體流場的影響研究。數(shù)值研究方面,高鐵鎖等[14]進(jìn)行了包含和不包含碳-酚醛燒蝕產(chǎn)物的再入小鈍錐流場計算,分析了燒蝕產(chǎn)物對流場電子數(shù)密度、溫度等流動參數(shù)的影響;魏叔如等[15]求解邊界層方程,開展了碳-碳材料燒蝕對電離邊界層影響的研究;董維中等[16]建立了碳-碳材料燒蝕流場的數(shù)值計算方法,研究了碳-碳材料燒蝕對再入體頭部區(qū)域溫度和熱流的影響,沒有對等離子體參數(shù)進(jìn)行分析。

總之,國內(nèi)外針對高溫等離子體流場開展的試驗(yàn)研究相對較少,主要通過數(shù)值模擬方法研究表面材料燒蝕產(chǎn)物對等離子體流場的影響,但數(shù)值方法存在多個難點(diǎn):1)不同材料的燒蝕反應(yīng)過程不同,燒蝕產(chǎn)物對等離子體流場的影響也有較大差別,需要針對不同材料建立相應(yīng)的燒蝕反應(yīng)模型;2)由于表面催化、氧化、燒蝕與熱解等復(fù)雜的氣體表面相互作用,氣體表面邊界條件會非常復(fù)雜,表面邊界方程在流場中的求解算法十分復(fù)雜,難以建立高效、魯棒的耦合計算方法,使得針對復(fù)雜外形應(yīng)用難度較大;3)由于進(jìn)入流場的燒蝕組分與高溫空氣組分差別較大,流場中的化學(xué)反應(yīng)過程會變得更加復(fù)雜,剛性問題和求解難度大大增加[17]。因此,目前文獻(xiàn)中的研究對象主要針對簡單的球錐再入體,材料主要針對傳統(tǒng)的泰氟隆材料和碳-碳材料。與球錐再入體不同,臨近空間飛行器一般會采用復(fù)雜的升力體外形,并長時間高馬赫數(shù)飛行,其熱環(huán)境特點(diǎn)為中高焓、低熱流,防隔熱設(shè)計往往會采用低密度的熱解炭化燒蝕材料[4],其燒蝕機(jī)理主要靠材料熱解吸熱和熱解氣體注入邊界層的熱阻塞效應(yīng)。調(diào)研發(fā)現(xiàn),目前針對升力體外形的等離子體流場,以及熱解炭化材料的燒蝕產(chǎn)物注入對等離子體流場的影響研究鮮有報道。

本文針對采用低密度熱解炭化防熱材料的臨近空間升力體飛行器,建立表面材料燒蝕與高溫等離子體流場數(shù)值計算方法,研究不同馬赫數(shù)下燒蝕產(chǎn)物注入流場對等離子體的影響規(guī)律。

1 計算方法

1.1 控制方程與求解方法

求解三維熱化學(xué)非平衡N-S方程[18],守恒形式的方程如下:

式中:Q為守恒變量;F、G、H與Fv、Gv、Hv分別對應(yīng)不同方向的對流項與黏性項;W為非平衡源項。

采用有限體積法求解控制方程,無黏通量離散求解采用Roe平均方法與Yee的對稱TVD限制器,并對特征值采用了Harten的熵修正[19];黏性通量離散求解采用中心差分格式;時間推進(jìn)采用全隱式LU-SGS方法,非平衡源項采用隱式處理[20]。

1.2 化學(xué)模型

當(dāng)飛行器表面溫度高于表面材料的燒蝕溫度門檻時就會發(fā)生燒蝕,此時流場中要考慮從壁面注入的熱解燒蝕產(chǎn)物的氣體組分;本文等離子體流場計算考慮了32種組分,包括N、O、N2、O2、NO、NO+、C2H、C2H2、CN、H、H2、HCN、CH、Si、SiO、SiO2、SiC、Si2、SiN,化學(xué)反應(yīng)模型見表1,其中組分的熱力學(xué)參數(shù)見文獻(xiàn)[21],輸運(yùn)參數(shù)見文獻(xiàn)[22-23],化學(xué)反應(yīng)速率常數(shù)見文獻(xiàn)[24-25]。

表1 高溫空氣/材料燒蝕組分的化學(xué)反應(yīng)模型Tab.1 Chemical reaction model of material ablative components

1.3 燒蝕壁面邊界條件

本文研究定常的燒蝕流場,對應(yīng)于壁面溫度最高,固體材料內(nèi)部溫度分布達(dá)到穩(wěn)態(tài)的情況,避免了耦合熱傳導(dǎo)方程的求解,并且不考慮燒蝕引起的壁面外形變化對流場的影響。

假設(shè)固體碳處于與表面氣體達(dá)到平衡的狀態(tài),存在如下平衡約束條件[26]:

式中:ρw為壁面處氣體密度;cw,C為壁面碳原子的質(zhì)量濃度;MC為碳的摩爾質(zhì)量;KC,C為異相反應(yīng)C(固相)→C(氣相)的化學(xué)平衡常數(shù);Tw為壁面溫度。(2)式的求解需要知道壁面碳原子的質(zhì)量濃度cw,C,此碳原子質(zhì)量分?jǐn)?shù)假設(shè)在壁面溫度、壓力和給定的元素比例條件下處于平衡。

化學(xué)平衡常數(shù)KC,C采用擬合關(guān)系式計算:

式中:T為氣體溫度。

對于某一元素k壁面處的質(zhì)量濃度可以從如下表面元素質(zhì)量平衡方程中獲得元素質(zhì)量分?jǐn)?shù):

壁面壓力pw采用如下壁面法向動量獲得:

式中:vw為壁面法向速度,vw來自如下質(zhì)量連續(xù)方程,

式中:qC為壁面熱流;α為輻射吸收系數(shù);qr為壁面向氣體的輻射熱流;ε為材料輻射系數(shù);σ為斯蒂芬-玻爾茲曼常數(shù);hw為壁面焓值。

式中:ρv、ρC分別為原始材料密度和炭化層密度。

2 方法驗(yàn)證

采用公開發(fā)表的電磁衰減測量項目C方案(RAM-C)飛行試驗(yàn)等離子體測試結(jié)果,對本文計算方法進(jìn)行驗(yàn)證。RAM-C模型為鈍球錐體,球頭半徑Rn=0.154 2 m,半錐角θ=9°,模型長度L=1.295 m。

飛行條件:速度U∞=7 650.0 m/s,高度H=71 km,壁面溫度Tw=1 500 K。本文主要考核流場中電子數(shù)密度分布,空氣化學(xué)反應(yīng)模型采用Park模型,考慮了N2、O2、NO、NO+、N、O、e-共計7種組分,熱力學(xué)模型采用兩溫度模型,壁面考慮完全非催化和完全催化兩種條件。

圖1給出了沿軸向位置電子數(shù)密度N峰值分布與飛行試驗(yàn)結(jié)果對比,橫軸為軸向位置與頭部半徑的無量綱距離,縱軸為電子數(shù)密度。由圖1可以看出,本文電子數(shù)密度峰值分布計算結(jié)果與飛行數(shù)據(jù)符合較好,誤差在半個數(shù)量級以內(nèi)。圖2給出了典型截面處法向電子數(shù)密度分布與飛行試驗(yàn)結(jié)果對比。由圖2可以看出:完全非催化和完全催化壁面條件會影響壁面附近的電子數(shù)密度,完全非催化壁結(jié)果高于完全催化壁結(jié)果,遠(yuǎn)離壁面的空間等離子體流場受影響較小?？傮w上,完全非催化壁結(jié)果更接近飛行試驗(yàn)測量值,表明本文等離子體流場計算方法在電子數(shù)密度分布計算中有較好的準(zhǔn)確性。

圖1 電子數(shù)密度峰值分布比較Fig.1 Comparison of peak distributions of electron number density

圖2 法向電子數(shù)密度分布比較Fig.2 Comparison of normal electron number density distributions

3 燒蝕產(chǎn)物的影響分析

3.1 計算模型及防熱材料屬性

針對某升力體氣動布局的臨近空間飛行器前體部分,開展燒蝕產(chǎn)物對等離子體流場的影響分析。圖3給出了計算外形和網(wǎng)格示意圖,采用六面體結(jié)構(gòu)化網(wǎng)格,在近壁面附近進(jìn)行了加密處理,網(wǎng)格總量150萬。計算來流工況為高度60 km、馬赫數(shù)分別為20、25、28,飛行攻角0°。

圖3 計算模型和網(wǎng)格Fig.3 Computational model and grid

為研究不同熱解炭化材料燒蝕產(chǎn)物對等離子體流場的影響,選取表2所示成分存在差異的同一類輕質(zhì)燒蝕防熱材料[27-28],其中材料1為文獻(xiàn)[27-28]中給出的阿波羅飛船采用的輕質(zhì)燒蝕材料;材料2和材料3分別在材料1基礎(chǔ)上進(jìn)行了工藝改進(jìn),分別去除材料中的氫元素和硅元素,分別記為輕質(zhì)燒蝕材料2和輕質(zhì)燒蝕材料3。3種材料密度均取544.63 kg/m3,熱解燒蝕后的炭化層密度取256.29 kg/m3。

表2 材料參數(shù)Tab.2 Material parameters

3.2 計算結(jié)果分析

3.2.1 純空氣等離子體流場規(guī)律分析

圖4給出了高度60 km、馬赫數(shù)27狀態(tài)的流場溫度和電子數(shù)密度分布。由圖4可以看出,高馬赫數(shù)來流在飛行器周圍形成了脫體激波,激波對等離子體流場的影響主要取決于激波后的流場溫度;對于頭部區(qū)域,激波強(qiáng)度大、壓縮性強(qiáng),波后溫度超過了10 000 K,空氣發(fā)生大規(guī)模離解和電離反應(yīng),流場中的O2分子幾乎完全離解,O、N原子的電離反應(yīng)是整個流場中等離子體產(chǎn)生的主要源頭,該區(qū)域的電子數(shù)密度最高,峰值超過了1×1014cm-3;飛行器身部區(qū)域的激波強(qiáng)度較弱,激波后的流場溫度約為5 000 K,該區(qū)域的空氣離解和電離程度較弱,飛行器身部流場的等離子體主要由頭部流場輸運(yùn)擴(kuò)散而來;受激波層影響,飛行器身部區(qū)域的等離子體層主要貼近飛行器壁面,其中前緣部位的等離子體層厚度相對迎風(fēng)面較薄,隨著流動向下游發(fā)展,流場溫度逐漸降低,飛行器身部電子數(shù)密度比頭部低2～3個量級。

圖4 流場溫度和電子數(shù)密度分布Fig.4 Temperature and electron number density distribution of flow field

圖5給出了頭部駐點(diǎn)線流場溫度分布。由圖5可以看出,駐點(diǎn)線上的振動溫度與平動溫度差別較大,流場的熱力學(xué)非平衡效應(yīng)顯著,由于激波層內(nèi)高溫區(qū)的氣體電離反應(yīng)是等離子體層形成的主要原因,需要采用兩溫度模型進(jìn)行等離子體流場計算。

圖5 沿駐點(diǎn)線溫度分布Fig.5 Temperature distribution along stagnation line

圖6給出了駐點(diǎn)線流場氣體電離組元的數(shù)密度分布。由圖6可以看出:馬赫數(shù)為20時,NO+離子是流場中電子數(shù)密度的主要貢獻(xiàn)源;馬赫數(shù)增大為25、27時,除了NO+離子,N+對電子數(shù)密度也有重要貢獻(xiàn),其中NO+主要由O與N原子的締合電離反應(yīng)機(jī)制產(chǎn)生,N+主要由N與其他組分的碰撞電離反應(yīng)機(jī)制產(chǎn)生。綜合上述分析可知,對流場中電子數(shù)密度貢獻(xiàn)最重要的離子是NO+和N+,其次依次為這些離子對電子數(shù)密度的貢獻(xiàn)很小。

圖6 沿駐點(diǎn)線電離組元的數(shù)密度分布Fig.6 Number density distribution of ionized components along stagnation line

3.2.2 燒蝕產(chǎn)物對等離子流場影響分析

圖7給出了材料3在高度60 km、馬赫數(shù)為20狀態(tài)的表面溫度和燒蝕質(zhì)量流率分布。由圖7可以看出:飛行器頭部和前緣溫度超過了2 000 K,該區(qū)域發(fā)生明顯燒蝕,其中頭部燒蝕最嚴(yán)重,燒蝕質(zhì)量流率超過了0.2 kg/(m2·s);飛行器身部迎風(fēng)面的溫度約為1 000 K,燒蝕質(zhì)量流率相對較小,沿流向燒蝕效應(yīng)逐漸減弱。

圖7 表面溫度和燒蝕質(zhì)量流率(高度60 km、馬赫數(shù)20,材料3)Fig.7 Surface temperature and ablation mass flow rate(height of 60 km,Mach number of 20,and Material 3)

圖8給出了流場中主要燒蝕組分的數(shù)密度。由圖8可見:注入流場的主要熱解燒蝕產(chǎn)物是CO和H2,其他燒蝕氣體組分的質(zhì)量分?jǐn)?shù)相對較小,并主要集中在頭部區(qū)域;下游飛行器身部的燒蝕效應(yīng)逐漸減弱,但隨著流動向下游發(fā)展,頭部燒蝕產(chǎn)物會隨流動的擴(kuò)散效應(yīng)輸運(yùn)到下游,燒蝕效應(yīng)對流場等離子分布的影響范圍擴(kuò)大。

圖8 對稱面燒蝕產(chǎn)物分布(高度60 km、馬赫數(shù)20,材料3)Fig.8 Distribution of ablated products on symmetry plane(height of 60 km,Mach number of 20,and Material 3)

圖9給出了頭部駐點(diǎn)線的流場參數(shù)對比。由圖9可以看出:燒蝕產(chǎn)物注入流場后,流場基本結(jié)構(gòu)沒有發(fā)生明顯變化,溫度和電子數(shù)密度峰值未受影響;由于壁面燒蝕的引射效應(yīng),使得頭部激波脫體距離增大,駐點(diǎn)線上的溫度、電子數(shù)密度和組分質(zhì)量分?jǐn)?shù)的分布也發(fā)生了一定變化,流場中高溫區(qū)域變大,O2和N2分子更多地分解成O和N原子,O原子參加氧化反應(yīng)被消耗而質(zhì)量分?jǐn)?shù)明顯減小,消耗的O原子主要生成了CO。

圖9 頭部軸線流場參數(shù)分布比較(高度60 km、馬赫數(shù)20,材料3)Fig.9 Distribution of axial flow field parameters on aircraft head(height of 60 km,Mach number of 20,and Material 3)

圖10和圖11分別給出了材料3在不同馬赫數(shù)下駐點(diǎn)線的溫度和電子數(shù)密度曲線。由圖10和圖11可以看出,對于材料3,隨著飛行馬赫數(shù)變大,表面燒蝕速率增強(qiáng),燒蝕產(chǎn)物進(jìn)入流場使得激波脫體距離增大,但駐點(diǎn)線上的電子密度峰值變化很小,這主要是因?yàn)榱鲌鲋须娮訑?shù)密度貢獻(xiàn)最重要的離子是NO+。圖12給出了材料3在不同馬赫數(shù)駐點(diǎn)線NO+離子數(shù)密度分布。從圖12中可以看出,不同馬赫數(shù)下燒蝕效應(yīng)對頭部駐點(diǎn)線上的NO+離子數(shù)密度影響相對較小導(dǎo)致的。

圖10 不同馬赫數(shù)下駐點(diǎn)線溫度分布(材料3)Fig.10 Distribution of temperature along stagnation line at different Mach numbers(Material 3)

圖11 不同馬赫數(shù)下駐點(diǎn)線電子數(shù)密度分布(材料3)Fig.11 Distribution of electron number density along stagnation line at different Mach numbers(Material 3)

圖12 不同馬赫數(shù)駐點(diǎn)線NO+離子數(shù)密度分布(材料3)Fig.12 Distribution of NO+ion number density along stagnation line at different Mach numbers(Material 3)

圖13給出了材料3的對稱面電子數(shù)密度峰值對比曲線。由圖13可以看出:隨著馬赫數(shù)增大,表面燒蝕速率增大,由于進(jìn)入流場的燒蝕產(chǎn)物增多,使得流場中的電子數(shù)密度有一定的增大;隨著流動從頭部向下游發(fā)展,飛行器身部的熱流會逐漸降低,表面材料的燒蝕速率逐漸減小,但燒蝕產(chǎn)物對飛行器身部流場電子數(shù)密度的影響逐漸增強(qiáng);下游0.5 m的位置,材料3的燒蝕產(chǎn)物對電子數(shù)密度的影響增量約為20%～30%。

圖13 對稱面電子數(shù)密度峰值對比(材料3)Fig.13 Comparison of electron number density peaks on symmetry plane(Material 3)

圖14給出了表2中不同成分的低密度防熱材料在高度60 km、馬赫數(shù)25狀態(tài)下飛行器表面中心線的燒蝕質(zhì)量流率對比曲線。由圖14可以看出:由于飛行器頭部的氣動加熱最為嚴(yán)酷,使得燒蝕質(zhì)量流率最大,材料1和材料3接近,約為0.55 kg/(m2·s),材料2為0.33 kg/(m2·s);與頭部相比,飛行器身部的燒蝕質(zhì)量流率會迅速降低約2～3個量級,表2中3種材料的燒蝕質(zhì)量流率均低于0.013 kg/(m2·s),其中材料2要略高于材料1和材料3。

圖14 不同材料飛行器表面對稱線的燒蝕質(zhì)量流率Fig.14 Ablation mass flow rate of surface center lines of aircraft with different materials

圖15、圖16分別為高度60 km、馬赫數(shù)25狀態(tài)下駐點(diǎn)線溫度和電子數(shù)密度對比曲線。由圖15和圖16可以看出:該狀態(tài)不同材料的燒蝕效應(yīng)存在差異,從而會影響頭激波的脫體距離,但駐點(diǎn)線的峰值溫度和峰值電子數(shù)密度幾乎不受影響;駐點(diǎn)線脫體距離從大到小依次為材料1、材料2和材料3;材料1和材料2的燒蝕質(zhì)量流率接近,但由于材料的成分不同,注入流場的燒蝕產(chǎn)物不同,使得兩種材料燒蝕后對流場激波脫體距離的影響程度存在差異。

圖15 不同材料條件下駐點(diǎn)線溫度分布Fig.15 Distribution of temperatures of different materials along stagnation line

圖16 不同材料條件下駐點(diǎn)線電子數(shù)密度分布Fig.16 Distribution of electron number density along stagnation line under different material conditions

圖17給出了高度60 km、馬赫數(shù)25狀態(tài)下對稱面電子數(shù)密度峰值對比曲線。由圖17可以看出:對稱面材料的燒蝕產(chǎn)物對流場峰值電子數(shù)密度影響沿流動向下游發(fā)展逐漸增強(qiáng),并且不同材料燒蝕產(chǎn)物對流場中的電子數(shù)密度峰值影響量不同,其中材料2的影響相對最大;在0.5 m位置,電子數(shù)密度的影響增量約為300%～320%。結(jié)合圖15表面燒蝕質(zhì)量流率分析可知,飛行器身部流場的電子數(shù)密度一方面受頭部區(qū)域燒蝕產(chǎn)物流動擴(kuò)散的影響,另一方面還要受當(dāng)?shù)氐奈⒘繜g產(chǎn)物的影響,二者均會使得流場中的電子數(shù)密度增大,其影響機(jī)理為燒蝕產(chǎn)物中CH、CO會發(fā)生反應(yīng)產(chǎn)生O原子,O原子與N原子發(fā)生締合電離反應(yīng),產(chǎn)生更多NO+。圖18所示為對稱面x=0.4 m位置流場中的NO+離子對比。由圖18可以看出,材料燒蝕產(chǎn)物注入后流場中的NO+離子顯著增大,從而使得流場中的電子數(shù)密度均比純空氣狀態(tài)增大。

圖17 不同材料條件下對稱面電子數(shù)密度峰值對比Fig.17 Comparison of peak values of electron number density on symmetry planes of different materials

圖18 x=0.4 m位置流場NO+離子數(shù)密度對比曲線Fig.18 Comparison curves of NO+ion number density in flow field for x=0.4 m

3 結(jié)論

本文建立了考慮燒蝕產(chǎn)物影響的等離子體流場數(shù)值計算方法,針對阿波羅飛船采用的輕質(zhì)燒蝕材料,開展了升力體外形的燒蝕產(chǎn)物對等離子體流場的影響規(guī)律研究。得到如下主要結(jié)論:

1)對于升力體飛行器的等離子體流場,影響最重要的離子是NO+和N+。

2)熱解炭化材料最主要的熱解燒蝕產(chǎn)物是CO和H2;燒蝕產(chǎn)物進(jìn)入流場會使得激波脫體距離增加,等離子體厚度增加。

3)燒蝕產(chǎn)物進(jìn)入流場,對頭部的峰值電子數(shù)密度影響很小,隨著流動向下游發(fā)展影響逐漸增大;影響機(jī)理為,燒蝕產(chǎn)物中CH、CO組分會發(fā)生反應(yīng),產(chǎn)生O原子,O原子與N原子發(fā)生締合電離反應(yīng),產(chǎn)生更多NO+,使得流場中電子數(shù)密度增大。

4)成分存在差異的同一類防熱材料,燒蝕產(chǎn)物對流場電子數(shù)密度的影響增量存在差異;本文研究的一類防熱材料,高度60 km、馬赫數(shù)為25狀態(tài)下飛行器下游0.5 m位置對稱面峰值電子數(shù)密度,材料2的增量約為300%～320%。