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

周 源,齊 強(qiáng),陳志剛,徐 明

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

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

周 源1,齊 強(qiáng)1,陳志剛1,徐 明2

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

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

絕熱層;熱燒蝕;時(shí)動(dòng)邊界

0 引言

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

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

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

1 計(jì)算模型

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

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

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

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

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

1.2 傳熱計(jì)算模型

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

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

(1)

2)熱解面上能量守恒

(2)

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

圖1 傳熱計(jì)算模型

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時(shí),絕熱層表面出現(xiàn)碳化層,根據(jù)邊界條件方程,有

(9)

(10)

7)初始條件

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

1.3 數(shù)值離散方法

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

(11)

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

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

在計(jì)算過(guò)程中,使用變空間步長(zhǎng)的差分方法,節(jié)點(diǎn)坐標(biāo)不斷更新。在每個(gè)時(shí)間步長(zhǎng)上,根據(jù)燒蝕速率來(lái)確定移動(dòng)邊界節(jié)點(diǎn)的位置,然后重新進(jìn)行離散,再計(jì)算溫度場(chǎng)。

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

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

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

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

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

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

圖4 燒蝕線和熱解線的相對(duì)位置

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

圖5 碳化層表面線燒蝕率

3 結(jié)論

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

2)預(yù)測(cè)-校正格式可以用于求解時(shí)動(dòng)邊界的瞬態(tài)傳熱問(wèn)題。

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

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

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

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

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

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

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

[8] 南京大學(xué)數(shù)學(xué)系計(jì)算數(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

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

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

V435

A