鄭子玙,史金光,徐 旸,趙 渭

(1.南京理工大學 能源與動力工程學院,江蘇 南京 210094;2.中國兵器工業(yè)試驗測試研究院,陜西 渭南 714200)

隨著社會的進步及科技的發(fā)展,現(xiàn)代戰(zhàn)爭對于武器有著越來越高的要求,具有高速遠程打擊能力的高超聲速武器應(yīng)運而生。高超聲速一般是指飛行馬赫數(shù)大于5,高超聲速武器依靠其自身的高動能產(chǎn)生足夠大的破壞力,其自身的穿甲、防空、反導能力相對于普通武器有著很大的提高,因此世界各軍事強國近年來均在積極開展高超聲速武器的研究。

高超聲速彈丸頭部母線形狀有球錐形、圓弧形、拋物型以及指數(shù)型等,壓心是決定飛行器穩(wěn)定配平飛行的關(guān)鍵因素,球錐旋成體外形的顯著特點是其壓心在高超聲速條件下十分穩(wěn)定,因此本文選取球錐體頭部外形進行數(shù)值模擬分析。由于此類彈丸速度大,?？蛇_到2 000 m/s及以上,在彈丸表面存在非常嚴重的氣動加熱現(xiàn)象,氣動加熱會導致彈體表面溫度升高,在到達一定溫度后彈體表面會產(chǎn)生燒蝕,其幾何外形變化不僅會增加阻力,還可能引起影響穩(wěn)定性方面的非設(shè)計力和力矩,進而對彈丸的精確打擊產(chǎn)生一定影響。因此高超聲速彈丸表面的氣動熱預測與熱防護技術(shù)已成為該類彈藥研制的關(guān)鍵技術(shù)之一。文獻[3]通過耦合計算氣動力、外彈道和氣動燒蝕,研究了超高速彈箭在飛行彈道上的氣動燒蝕對外彈道的影響;文獻[4]采用有限體積TVD差分格式求解三維歐拉方程,數(shù)值模擬了來流馬赫數(shù)>5時的彈丸表面流場,利用數(shù)值計算求得相關(guān)參量,并代入求解熱流密度的經(jīng)驗公式,得出彈丸頭部的熱流密度和溫度分布情況;文獻[5]數(shù)值模擬了高超聲速彈丸的流場與結(jié)構(gòu)傳熱,及彈丸的燒蝕變形過程,并對彈丸燒蝕過程特性進行了分析;文獻[6]對某彈丸進行了氣動加熱計算,模擬了彈丸表面的氣動熱及對內(nèi)部裝藥傳熱的影響;文獻[7]采用流固耦合的方法,對不同飛行條件下彈丸頭部的氣動燒蝕進行了模擬,編寫了UDF仿真程序,實現(xiàn)了流固耦合面共同移動,得到了與工程算法較吻合的計算結(jié)果;文獻[8]對于不同湍流模型公式進行了數(shù)值比較,以轉(zhuǎn)捩點作為判據(jù)得出了最符合真實邊界氣體流動情況的湍流模型。文獻[9]通過數(shù)值算例驗證了基于參考焓法的氣動加熱工程算法的可行性,實現(xiàn)了高超聲速三維翼面溫度的準確預測,該方法可用于進行高超聲速飛行器氣動熱分析及熱防護設(shè)計。文獻[10]用不同熱力學模型進行分析,說明了熱力學模型對熱流密度的預測有一定的影響。文獻[11]利用不同的數(shù)值方法對高超聲速飛行器的壁面熱流密度進行計算并與風洞實驗數(shù)據(jù)做對比,研究了各數(shù)值計算程序?qū)Ω叱曀偌げǖ念A測能力。

目前國內(nèi)對于彈丸表面氣動熱的研究主要基于氣流加熱使得壁面升溫,很少考慮彈丸表面因氧化化學反應(yīng)產(chǎn)生的吸放熱現(xiàn)象,這就使得準確預測彈丸表面的氣動熱以及彈丸內(nèi)部的傳熱導熱問題變得困難。本文基于流固耦合,在考慮熱化學反應(yīng)的同時,建立了熱化學氧化燒蝕模型,通過對某典型結(jié)構(gòu)的高超聲速彈丸進行氣動熱計算,分析有、無化學反應(yīng)對于彈丸表面氣動熱的影響,為高超聲速彈丸的氣動熱預測與熱防護提供參考。

1 計算模型

1.1 物理模型

在彈丸高速飛行過程中,彈丸頭部所承受的壓力高,氣動熱環(huán)境惡劣,頭部處氧化反應(yīng)劇烈,因此截取彈丸頭部進行數(shù)值模擬。其物理模型如圖1所示,該彈丸頭部為鈍體模型,其頭部半徑為5 mm,頭部半錐角為5.3°,截取頭部錐體長度為380 mm。

圖1 彈丸頭部錐體模型

1.2 網(wǎng)格劃分

為使計算更加精確,采用結(jié)構(gòu)四邊形網(wǎng)格對模型進行網(wǎng)格劃分,整體及局部網(wǎng)格劃分如圖2和圖3所示,在近壁面處對網(wǎng)格進行邊界層加密,與離壁面最近的一層網(wǎng)格高度為0.001 mm,其附近的網(wǎng)格增長比為1.1,保證了壁面附近的+值小于1。

圖2 彈丸整體網(wǎng)格劃分

圖3 彈丸局部網(wǎng)格劃分

1.3 邊界條件

外邊界條件設(shè)置為壓力遠場邊界條件,外部壁面條件設(shè)置為內(nèi)部面邊界,內(nèi)部壁面條件設(shè)置為固體無滑移邊界,與影子壁面通過耦合邊界條件實現(xiàn)流固耦合之間的熱交換。

2 計算方法

2.1 基本假設(shè)

一般情況下,彈丸頭部為一旋成體,材料分布均勻,且具有各向同性。因此為了簡化計算,本文作以下假設(shè):①彈丸質(zhì)量分布均勻;②彈丸各個部位材料的熱化學性質(zhì)不變,加熱方式為空氣摩擦及表面化學反應(yīng);③主要進行短時間的氣動熱預測,并未考慮物體外形的變化。

2.2 數(shù)學模型

1)控制方程。

數(shù)值模擬采用N-S方程,在直角坐標系下,其中質(zhì)量守恒方程為

(1)

式中:為密度,為時間,為速度矢量。

動量守恒方程為

(2)

(3)

式中:為壓力;,分別為重力加速度在方向和方向的分量;,分別為速度在方向和方向的分量;為分子黏性系數(shù);為應(yīng)力張量。

能量方程為

(4)

式中:為流體比焓,為邊界層向流體的傳熱,為力對流體微元體所做的功。

2)湍流模型。

本文選擇SST-模型作為基本的湍流模型。其流動方程為

(5)

(6)

式中:,分別為和的湍流動能產(chǎn)生項;,分別為和的發(fā)散項;和分別為和的擴散率;和為擴散產(chǎn)生的湍流;為正交發(fā)散項。

3)輻射模型。

為準確描述彈體表面?zhèn)鳠釥顩r,本文采用離散坐標輻射(DO)模型來模擬彈體及邊界層中的輻射傳熱。

4)氧化傳熱模型。

假設(shè)彈丸的材料為金屬鎢,彈丸在高速飛行過程中,加熱的鎢與空氣中的氧氣發(fā)生反應(yīng),在彈體表面會形成三氧化鎢,其基本的化學反應(yīng)方程式為

2W+3O→2WO

(7)

文獻[14]研究發(fā)現(xiàn),高溫下鎢的氧化反應(yīng)速率可用下述方程描述:

(8)

(9)

其反應(yīng)的化學生成熱為

=-ΔΔ

(10)

式中:Δ為參與表面反應(yīng)顆粒的質(zhì)量,Δ為表面反應(yīng)釋放的熱量。

化學反應(yīng)生成熱一部分被彈體吸收,用表示;另一部分傳導至附面層中,用表示。

彈體表面化學反應(yīng)熱平衡方程為

(11)

式中:為單位時間內(nèi)參與反應(yīng)的鎢顆粒質(zhì)量,為定壓比熱容,為溫度,為傳至彈體中的化學反應(yīng)熱在總化學反應(yīng)生成熱中所占的百分比,為對流換熱系數(shù),為單位參與反應(yīng)鎢顆粒表面積,為外界溫度,為輻射給熱。

基于已有的氧化模型,在彈體表面鎢的氧化區(qū)域建立氧化傳熱區(qū)域。氧化傳熱的有效熱傳導系數(shù)由區(qū)域中流體熱傳導率和固體熱傳導率通過體積平均得到,表示為

=+(1-)

(12)

式中:為氧化傳熱區(qū)域內(nèi)部單位體積內(nèi)氧氣體積所占的百分比。

5)熱傳導模型。

在高溫邊界層中向物體表面的傳熱主要由導熱項、擴散傳熱項及輻射傳熱項組成。規(guī)定邊界層向固體導熱方向為正。

氣體導熱用傅里葉公式表示為

=-·grad

(13)

擴散傳熱是指多組元的混合氣體由于組元的濃度不均勻形成擴散流動而傳遞的熱量。令為組元的質(zhì)量濃度,為組元的焓,它包括熱焓和化學焓,表示為

(14)

令為組元的擴散速度,即組元對質(zhì)點速度的相對速度,,為組元的比熱容,為組元的擴散系數(shù),如果忽略由溫度引起的熱擴散影響,則擴散速度只與濃度有關(guān),表示為

(15)

僅由質(zhì)量擴散所傳遞的熱能為

=∑

(16)

組元的密度為

=

(17)

在加入氧化反應(yīng)后,壁面?zhèn)鳠峥捎?部分構(gòu)成,如圖4所示。圖中,為化學反應(yīng)熱中傳導至固體的部分,為化學反應(yīng)熱中傳導至流體的部分,為導熱,為擴散傳熱,為輻射傳熱至固體的部分,為輻射傳熱至流體的部分。

圖4 熱傳導示意圖

由上面表述可知,在高溫邊界層內(nèi)向物體表面?zhèn)鳠岬姆匠虨?/p>

(18)

式(18)將包含化學反應(yīng)放熱的能量源項添加到能量方程中,通過迭代進行熱量從流體傳遞到固體的熱傳導耦合求解。

3 計算結(jié)果及分析

3.1 數(shù)值方法驗證

為驗證本文采用傳熱計算方法的可靠性,以二維圓管高超聲速繞流耦合換熱為例進行驗證計算,并與實驗結(jié)果進行對比。

驗證模型如圖5所示,圓管內(nèi)徑為0.025 4 m,外徑為0.038 1 m。圓管材料為標準鋼,密度=8 030 kg/m,比熱=502.48 J/(kg·K),熱傳導系數(shù)=16.27 W/(m·K)?？諝鈦砹黢R赫數(shù)=647,壓強=648.1 Pa。驗證模型網(wǎng)格如圖6所示,對管壁附近網(wǎng)格進行加密處理。計算結(jié)果如圖7所示,將管壁附近熱流密度無量綱化后與實驗數(shù)據(jù)進行比較,為二維圓管外壁某一點和圓管中心連線與水平的夾角,從圖中可以看出二者吻合較好,表明本文采用的傳熱計算方法可行。

圖5 二維圓管模型示意圖

圖6 二維圓管模型網(wǎng)格示意圖

圖7 本文熱流密度結(jié)果與實驗數(shù)據(jù)對比圖

3.2 網(wǎng)格無關(guān)性驗證

為對某外形鎢制高超聲速炮彈進行彈丸表面氣動熱對比分析,首先給出鎢及三氧化鎢材料特性,如表1所示。

表1 鎢及三氧化鎢材料特性

本文對于算例開展了網(wǎng)格無關(guān)性驗證,在海平面標準條件下,設(shè)計了4套網(wǎng)格計算駐點溫度,對于網(wǎng)格無關(guān)性進行了驗證分析,網(wǎng)格總數(shù)分別為12萬、18萬、27萬、40萬,對應(yīng)駐點熱流分別為19.3 MW/m,19.7 MW/m,20.4 MW/m和20.6 MW/m。對比27萬網(wǎng)格數(shù)量的計算結(jié)果,12萬、18萬、40萬網(wǎng)格計算結(jié)果誤差分別為5.39%,3.43%,0.98%,綜合計算精度和計算速度考慮,本文擬采用27萬網(wǎng)格精度進行計算。

3.3 考慮有、無化學反應(yīng)對比

本文采用參考來流條件為海平面標準氣象條件,大氣壓力為101 325 Pa,彈丸表面初始溫度為288.15 K,飛行馬赫數(shù)為5.5。

圖8和圖9分別表示了彈丸表面不考慮化學反應(yīng)和考慮化學反應(yīng)的溫度云圖,從圖中可以看到,整個彈體溫度最高出現(xiàn)在頭部,說明此處為氣動熱最為嚴重的區(qū)域。有、無考慮化學反應(yīng)發(fā)生的彈體內(nèi)部溫度等值線均呈環(huán)狀分布,沿中心軸線溫度逐漸降低,溫度梯度逐漸增大,這是由于在彈丸飛行過程中,外部氣流溫度很高,對彈丸表面產(chǎn)生加熱現(xiàn)象,其溫升很快,彈身部分由于固體材料傳熱系數(shù)的限制溫升較慢。相較于不考慮化學反應(yīng),考慮化學反應(yīng)發(fā)生的彈頭部分及錐身表面處溫度更高,且彈身傳熱現(xiàn)象更為明顯。

圖8 彈丸表面未考慮化學反應(yīng)時的溫度云圖

圖9 彈丸表面考慮化學反應(yīng)時的溫度云圖

飛行1.5 s后,彈丸表面有、無考慮化學反應(yīng)時的溫度對比如圖10所示。從圖10可以看出,當考慮彈丸表面發(fā)生的氧化反應(yīng)后,彈丸駐點處溫度相較于沒有考慮化學反應(yīng)時有了明顯提高。在不考慮氧化反應(yīng)時,彈丸駐點處最高溫度為1 326.58 K,在考慮了彈體表面氧化反應(yīng)后駐點處最高溫度為1 549.20 K。究其原因為:彈丸在空氣中飛行時,彈體表面和空氣中的氧氣在高溫下發(fā)生劇烈的氧化反應(yīng),是一個放熱過程,其放出的熱量與高溫邊界層中氣體的擴散傳熱以及導熱累加,向壁面?zhèn)鬟f使得壁面溫度進一步升高。由于考慮氧化反應(yīng),相較于不考慮化學反應(yīng)的情況,經(jīng)歷同樣的飛行時間后彈丸表面溫度更高,會更快出現(xiàn)燒蝕。彈丸表面的熱流分布如圖11所示。

圖10 彈丸表面有、無考慮化學反應(yīng)溫度對比圖

圖11 彈丸表面有、無考慮化學反應(yīng)熱流分布對比圖

有、無考慮氧化反應(yīng)發(fā)生的彈體表面熱流密度分布規(guī)律大體相同,二者均為駐點處熱流密度最大,錐體部分熱流逐漸減少。在駐點后,由于轉(zhuǎn)捩區(qū)的存在,熱流值下降的過程中有一短暫升高?；瘜W反應(yīng)的放熱現(xiàn)象使得空氣中的溫度較不考慮氧化反應(yīng)時更高,其與彈體之間的溫差更大,彈體表面單位時間單位面積的傳熱量更多,導致熱流密度更大。在考慮氧化反應(yīng)后,錐體部分熱流密度有了很大的提高,氧化反應(yīng)的進行有顯著增加壁面熱流效應(yīng)的作用。

圖12為氧化傳熱區(qū)域內(nèi)三氧化鎢的質(zhì)量分數(shù)分布圖。從圖中可以看出,在彈丸頭部三氧化鎢質(zhì)量分數(shù)較高,說明氧化反應(yīng)在頭部較為劇烈,這也就是在彈丸表面有、無考慮化學反應(yīng)熱流對比圖中,頭部駐點區(qū)域彈丸熱流密度差異出現(xiàn)較大的原因。

圖12 氧化傳熱區(qū)域內(nèi)三氧化鎢質(zhì)量分數(shù)分布圖

3.4 不同飛行高度對比

為了對比分析在不同飛行高度下,有、無考慮化學反應(yīng)時對彈丸頭部氣動熱的影響,設(shè)計以下算例進行數(shù)值模擬分析,如表2所示,飛行馬赫數(shù)均為5.5(算例1),飛行時間為1.5 s。

表2 不同飛行高度條件下計算參數(shù)

通過對上述算例進行數(shù)值模擬,可以得到11組不同飛行高度下有、無考慮化學反應(yīng)時彈丸駐點處的溫度,如表3所示。為了量化駐點處考慮化學反應(yīng)對氣動熱的影響,表3中列出了不同高度下有、無考慮化學反應(yīng)時駐點處的溫度差。溫度差越大,說明不考慮化學反應(yīng)對于彈丸表面氣動熱的準確預測影響越大。

表3 不同飛行高度駐點處的溫度及溫差對比

從表3可以看出,隨著飛行高度的增加,彈丸駐點處及彈身表面的溫度會降低,飛行1.5 s后各高度上考慮化學反應(yīng)彈丸駐點處溫度較未考慮化學反應(yīng)駐點處溫度高約200 K。這是因為隨著飛行高度的增加,空氣密度等大氣參數(shù)降低,彈丸表面周圍氣體與彈體表面交換的熱量降低;考慮化學反應(yīng)后由于增加了一部分化學反應(yīng)產(chǎn)生的熱量,使得彈身表面及駐點處的溫度較未考慮化學反應(yīng)時有明顯增加。

3.5 不同來流速度對比

為了對比分析在不同飛行馬赫數(shù)下,有、無考慮化學反應(yīng)對彈丸頭部氣動熱的影響,設(shè)計3組不同飛行馬赫數(shù)的對比算例進行數(shù)值模擬分析,其外界大氣壓強均為101 325 Pa,彈丸表面初始溫度均為288.15 K,飛行馬赫數(shù)分別為5,6,6.5。

將上述3組算例與算例1結(jié)合起來,可以得到4組不同馬赫數(shù)下彈體表面溫度的對比結(jié)果。表4給出不同馬赫數(shù)下彈丸駐點處考慮有、無化學反應(yīng)時的溫度及溫度差對比。

表4 不同飛行馬赫數(shù)駐點處的溫度及溫差對比

從表4可以看出,隨著來流馬赫數(shù)的增加,彈丸駐點處的溫度增加,且溫度差值變大,氧化反應(yīng)變得更加劇烈,由此產(chǎn)生的氣動加熱現(xiàn)象更為嚴重,彈丸表面將更快到達熔點溫度并產(chǎn)生燒蝕現(xiàn)象。究其原因:隨著來流馬赫數(shù)的增加,彈丸表面氣動加熱現(xiàn)象更為劇烈,表面溫度更高;彈丸表面溫度升高后,反應(yīng)物的分子獲得能量,使得原本不會發(fā)生化學反應(yīng)的一部分分子獲得了更高的活化能而發(fā)生反應(yīng),加快了化學反應(yīng)速率,釋放了更多的熱量,因此隨著來流馬赫數(shù)的增加,駐點處的溫度相對提升更大。彈丸飛行馬赫數(shù)的增大導致彈丸表面化學反應(yīng)更為劇烈,因此在進行高馬赫數(shù)彈丸表面氣動熱預測時,應(yīng)充分考慮化學反應(yīng)對其表面溫度的影響。

4 結(jié)束語

本文基于流固耦合,在考慮熱化學反應(yīng)的同時,建立了熱化學氧化傳熱模型,對某外形高超聲速彈丸頭部進行了氣動熱計算,分析了有、無考慮化學反應(yīng)對于彈丸表面氣動熱的影響,通過對比不同飛行高度、不同飛行馬赫數(shù)等多種條件下的彈丸表面氣動熱計算結(jié)果,得到以下結(jié)論:

①在考慮化學反應(yīng)后,彈丸表面溫度及其表面熱流密度均有明顯提高,其中由于彈丸頭部駐點附近化學反應(yīng)劇烈,駐點處的溫度差異最大?？紤]化學反應(yīng)發(fā)生的彈丸頭部及錐身表面處溫度更高,且彈身傳熱現(xiàn)象更為明顯。相較于不考慮化學反應(yīng),考慮化學反應(yīng)彈丸表面熱流密度有較大提高,因此為準確預測高超聲速彈丸的氣動熱,需要考慮彈體表面的化學反應(yīng)。

②在飛行馬赫數(shù)為5.5,飛行時間為1.5 s的情況下,隨著飛行高度的增加,彈丸駐點處及彈身表面的溫度會降低,但各高度上考慮化學反應(yīng)彈丸駐點處溫度較未考慮化學反應(yīng)駐點處溫度高約200 K。為減小高超聲速彈丸長時間飛行時彈身表面的氣動熱,可設(shè)計高超聲速彈丸在較大高度上長時間巡航飛行,有利其抗氣動燒蝕與熱防護。

③隨著來流馬赫數(shù)的增加,化學反應(yīng)產(chǎn)生的熱量越多,有、無考慮化學反應(yīng)駐點處的溫度差變大。在標準海平面條件下,彈丸飛行馬赫數(shù)為6時,不考慮化學反應(yīng)駐點處溫度為1 561.63 K,考慮化學反應(yīng)駐點處溫度為1 800.49 K,后者已經(jīng)到達材料熔點溫度而產(chǎn)生燒蝕現(xiàn)象。因此,對于高超聲速這類彈箭,其飛行馬赫數(shù)較大,具有飛行距離遠與長時間滯空飛行等特點,為了有效降低其表面的氣動燒蝕,準確的氣動熱預測和熱防護對該類彈箭的設(shè)計具有重要的意義。