周 源,齊 強,陳志剛,徐 明

(1 海軍航空工程學(xué)院,山東煙臺 264001;2 91515部隊,海南三亞 572000)

燃?xì)獍l(fā)生器絕熱層燒蝕數(shù)值仿真*

周 源1,齊 強1,陳志剛1,徐 明2

(1 海軍航空工程學(xué)院,山東煙臺 264001;2 91515部隊,海南三亞 572000)

為了研究燃?xì)獍l(fā)生器燃燒室絕熱層的傳熱燒蝕過程,建立了熱化學(xué)燒蝕計算模型和絕熱層傳熱計算模型?？紤]到燒蝕過程中絕熱層邊界的移動,根據(jù)預(yù)測-校正格式對模型進行離散,并采用擬牛頓法對絕熱層的燒蝕過程及傳熱過程進行了耦合計算。計算得到了絕熱層表面溫度、燒蝕率和燒蝕厚度等參數(shù)的變化規(guī)律。計算結(jié)果表明,預(yù)測-校正格式和擬牛頓法可以用于求解時動邊界的瞬態(tài)傳熱模型。

絕熱層;熱燒蝕;時動邊界

0 引言

燃?xì)獍l(fā)生器是一種小型固體火箭發(fā)動機,是彈射發(fā)射系統(tǒng)的工質(zhì)源和能源。藥柱在燃燒室內(nèi)點火燃燒,產(chǎn)生高溫高壓燃?xì)馐沟萌紵颐媾R極其惡劣的內(nèi)部熱環(huán)境。因此,研究燃燒室內(nèi)壁絕熱層的傳熱燒蝕過程對確保燃?xì)獍l(fā)生器安全工作具有重要意義。

對絕熱層的傳熱燒蝕問題已有很多的研究,文獻[1]和文獻[2]對熱解型絕熱材料的燒蝕過程建立了物理結(jié)構(gòu)模型,文獻[3]和文獻[4]采用有限差分法對絕熱層的瞬態(tài)傳熱進行計算,文獻[5]對于時動邊界上的熱傳導(dǎo)問題的求解提出了人工邊界的方法。

文中對燃燒室絕熱層的燒蝕傳熱過程建立了碳化層—熱解面—基體層的結(jié)構(gòu)模型,采用預(yù)測-校正格式對時動邊界模型進行離散求解,為絕熱層的燒蝕傳熱研究提供了理論依據(jù)。

1 計算模型

1.1 熱化學(xué)燒蝕計算模型

在燃?xì)獍l(fā)生器內(nèi),燃?xì)庵袇⒓訜g反應(yīng)的氣體組分主要是CO2、H2O和H2,將其余不參加燒蝕反應(yīng)的惰性氣體組分折合成N2,如表1所示。那么,熱化學(xué)反應(yīng)燒蝕率的計算采用以下3個化學(xué)反應(yīng)方程:

表1 燃?xì)庵械闹饕磻?yīng)氣體組分

根據(jù)質(zhì)量守恒定理,燒蝕反應(yīng)氣體組分來自于燃?xì)庵髁鲾U散和絕熱層熱解氣體,在燒蝕表面發(fā)生反應(yīng)后被燃?xì)庵髁鲙ё摺峄瘜W(xué)燒蝕程序方程組[6]表示如下:

(1+F)K6,w=K6,g+Bf6

1.2 傳熱計算模型

絕熱層在工作過程中形成碳化層—熱解面—原始材料層的結(jié)構(gòu),如圖1所示。將絕熱層內(nèi)部的導(dǎo)熱過程簡化為一維大平板的瞬態(tài)導(dǎo)熱[1],不考慮絕熱材料的熱膨脹及其引起的熱應(yīng)力,且材料參數(shù)隨溫度的變化忽略不計[7],由此建立傳熱計算模型。

1)碳化層內(nèi)的瞬態(tài)導(dǎo)熱

(1)

2)熱解面上能量守恒

(2)

式中:ΔHp表示絕熱材料的熱解潛熱;λ1表示原始材料層的導(dǎo)熱系數(shù);Tp表示絕熱材料的熱解溫度。

圖1 傳熱計算模型

3)原始材料層內(nèi)的瞬態(tài)導(dǎo)熱

(3)

式中:ρ1表示原始材料層的密度;c1表示原始材料層的比熱;λ1表示原始材料層的導(dǎo)熱系數(shù)。

4)原始材料層與殼體層之間的界面上能量守恒

(4)

式中,λ2表示殼體層的導(dǎo)熱系數(shù)。

5)殼體層內(nèi)的瞬態(tài)導(dǎo)熱

(5)

式中:ρ2表示殼體層的密度;c2表示殼體層的比熱;λ2表示殼體層的導(dǎo)熱系數(shù)。

6)邊界條件

①根據(jù)假設(shè),燃?xì)獍l(fā)生器工作期間與外界是絕熱的。殼體層與外部環(huán)境之間的界面上能量守恒,有

(6)

②當(dāng)Tw

(7)

Qin=Qrad+Qcon

(8)

③當(dāng)Tw≥Tp時,絕熱層表面出現(xiàn)碳化層,根據(jù)邊界條件方程,有

(9)

(10)

7)初始條件

絕熱層內(nèi)部初始溫度是常數(shù),取環(huán)境溫度。

1.3 數(shù)值離散方法

在燒蝕過程中絕熱層的厚度不斷減少,對于具有時動邊界的瞬態(tài)傳熱問題,采用預(yù)測-校正格式[8]進行離散求解。如圖1所示,對每一層結(jié)構(gòu)進行等分離散,則有:

(11)

通過與其它數(shù)值離散格式對比,該格式具有二階精度且無條件穩(wěn)定。

2 計算結(jié)果與分析

在計算過程中,使用變空間步長的差分方法,節(jié)點坐標(biāo)不斷更新。在每個時間步長上,根據(jù)燒蝕速率來確定移動邊界節(jié)點的位置,然后重新進行離散,再計算溫度場。

考慮到計算時間步長較小(取0.001 s),i時刻和i+1時刻初值變化很小,將i時刻程序的解作為i+1時刻程序初值[6],再采用擬牛頓法求解,以縮短計算時間。

圖2給出了絕熱層單位面積上熱解氣體質(zhì)量流率隨時間的變化規(guī)律。在燃?xì)獍l(fā)生器開始工作的極短時間內(nèi),原始材料層不發(fā)生熱解,表面未形成碳化層。隨著表面溫度的升高,絕熱材料開始熱解并產(chǎn)生極薄的碳化層。這時,通過熱解面?zhèn)魅虢^熱材料內(nèi)部的熱流很大,熱解速率較大。在0.009 s時,熱解氣體質(zhì)量流率曲線出現(xiàn)了一個很大的峰值,熱解氣體質(zhì)量流率達到6.697 1 kg/(m2·s);隨著絕熱材料的不斷熱解,碳化層厚度增加并逐漸穩(wěn)定,熱解速率減小,最后逐漸趨于定值。熱解氣體平均質(zhì)量流率為0.519 4 kg/(m2·s)。

圖2 熱解氣體質(zhì)量流率

圖3給出了絕熱層內(nèi)表面的溫度變化規(guī)律。在開始階段,溫度迅速升高,溫度越高,變化速率越小。在達到約2 500 K時,溫度趨于平衡。

圖3 絕熱層內(nèi)表面溫度

圖4給出了燃?xì)獍l(fā)生器工作過程中燒蝕線與熱解線的相對位置。在開始階段,絕熱層不發(fā)生熱解和碳化。直到絕熱層溫度達到熱解溫度,碳化層厚度迅速增大,并隨著溫度的升高逐漸減緩。從圖中可以看出,熱解速率明顯比燒蝕速率快很多。

圖4 燒蝕線和熱解線的相對位置

圖5給出了絕熱層內(nèi)表面碳化層的線燒蝕率隨時間變化曲線。在絕熱層表面未形成碳化層的時候,燒蝕率為零。在絕熱層內(nèi)表面溫度的升高,以及燃燒室壓力的增加,根據(jù)Arrhenius定律,碳化層燒蝕率逐漸增大。在后效段,燃燒室壓力急劇下降,線燒蝕率也急劇下降。根據(jù)曲線求得平均線燒蝕率為0.069 3 mm/s。

圖5 碳化層表面線燒蝕率

3 結(jié)論

1)文中建立的熱化學(xué)燒蝕模型和傳熱計算模型以及所采用的離散格式和計算方法實現(xiàn)了燒蝕變形與傳熱的雙向耦合。

2)預(yù)測-校正格式可以用于求解時動邊界的瞬態(tài)傳熱問題。

[1] 徐善瑋, 侯曉, 張宏安. 固體火箭發(fā)動機內(nèi)絕熱層燒蝕質(zhì)量損失計算 [J]. 固體火箭技術(shù), 2003, 26(3): 28-31.

[2] 張濤, 孫冰. 三維燒蝕內(nèi)部熱響應(yīng)數(shù)值計算研究 [J]. 宇航學(xué)報, 2012, 33(3): 298-304.

[3] 李宛珊, 王文洽. 二維熱傳導(dǎo)方程的有限差分區(qū)域分解算法 [J]. 山東大學(xué)學(xué)報: 理學(xué)版, 2011, 46(12): 1-5.

[4] 馮立偉. 熱傳導(dǎo)方程幾種差分格式的MATLAB數(shù)值解法比較 [J]. 沈陽化工大學(xué)學(xué)報, 2011, 25(2): 179-182.

[5] 張雪艷, 周愛霞. 時動邊界上熱傳導(dǎo)問題的求解方法 [J]. 廊坊師范學(xué)院學(xué)報: 自然科學(xué)版, 2010, 10(2): 9-12.

[6] 張斌, 劉宇, 王長輝，等. 長時間工作固體火箭發(fā)動機燃燒室熱防護層燒蝕計算 [J]. 固體火箭技術(shù), 2011, 34(2): 189-192.

[7] 楊春杰. 固體火箭發(fā)動機后效沖量研究 [D]. 長沙: 國防科技大學(xué), 2011.

[8] 南京大學(xué)數(shù)學(xué)系計算數(shù)學(xué)專業(yè). 偏微分方程數(shù)值解法 [M]. 北京: 科學(xué)出版社, 1979: 145-182.

Numerical Simulation for Ablation of Insulator in Gas-generator

ZHOU Yuan1, QI Qiang1, CHEN Zhigang1, XU Ming2

(1 Naval Aeronautical and Astronautical University, Shandong Yantai 264001, China; 2 No.91515 Unit, Hainan Sanya 572000, China)

In order to study thermal behavior and ablation of insulator in the combustor of gas-generator, the thermo-chemistry ablation model and heat transfer calculation model were established. Considering the fact that the boundary of insulator is moving during the ablation, the models were discredited with predicted-correction method. Coupling calculations were carried out for insulator ablation and temperature field with quasi-Newton method. Temperature of inner wall, ablation velocity and ablation thickness were calculated. The results show that the predicted-correction method and the quasi-Newton method can be applied in unsteady heat transfer and ablation problem with moving boundary.

insulator; ablation; moving boundary

2014-03-14

國家自然科學(xué)基金(51005242)資助

周源(1979-),男,安徽濉溪人,講師,博士研究生,研究方向:兵器發(fā)射理論與技術(shù)。

V435

A