葉 青，余永剛

(南京理工大學(xué)能源與動(dòng)力工程學(xué)院，江蘇 南京 210094)

引 言

隨著火箭導(dǎo)彈等武器在戰(zhàn)場(chǎng)上使用環(huán)境的日漸苛刻，對(duì)固體火箭發(fā)動(dòng)機(jī)的生存能力提出了更高的要求。在現(xiàn)代戰(zhàn)爭(zhēng)中，彈藥在未投入使用前，存在遇襲或自身著火引起燃燒或發(fā)生爆炸的危險(xiǎn)。因此，降低彈藥易損性和提高彈藥熱安全性引起國(guó)內(nèi)外學(xué)者的高度重視。

烤燃試驗(yàn)和烤燃數(shù)值模擬是研究和評(píng)估彈藥及高能材料熱安全性的常用方法。國(guó)內(nèi)學(xué)者已針對(duì)炸藥進(jìn)行了大量烤燃試驗(yàn)研究，探討各因素對(duì)炸藥烤燃響應(yīng)的影響。高峰等[1]通過(guò)烤燃試驗(yàn)研究了物理界面及界面厚度對(duì)烤燃過(guò)程的影響。與炸藥的烤燃試驗(yàn)研究相比，針對(duì)固體推進(jìn)劑的烤燃試驗(yàn)研究起步較晚，但已有不少成果。陳中娥等[2]分析了HTPB推進(jìn)劑在慢速烤燃情況下的熱分解特性，認(rèn)為高氯酸銨(AP)熱分解過(guò)程形成的多孔性形貌是導(dǎo)致HTPB 推進(jìn)劑慢速烤燃響應(yīng)劇烈的主要因素；趙孝彬等[3]通過(guò)慢速烤燃試驗(yàn)研究了HTPE、GAP固體推進(jìn)劑慢速烤燃特性的影響因素，發(fā)現(xiàn)推進(jìn)劑成分、燃速、升溫速率及夾板約束對(duì)HTPE推進(jìn)劑和GAP推進(jìn)劑慢速烤燃的響應(yīng)程度均有影響；楊后文等[4]針對(duì)裝填A(yù)P/HTPB固體推進(jìn)劑的某小型固體火箭發(fā)動(dòng)機(jī)建立了二維烤燃簡(jiǎn)化模型，發(fā)現(xiàn)隨著火焰溫度的提高，著火延遲期縮短，著火溫度逐漸增大。烤燃試驗(yàn)耗資高、費(fèi)時(shí)長(zhǎng)、危險(xiǎn)性大且不確定因素多，而采用數(shù)值模擬方法可以更安全廉價(jià)地研究烤燃特性，且可以進(jìn)行預(yù)測(cè)性研究。

目前，國(guó)內(nèi)外對(duì)小型固體火箭發(fā)動(dòng)機(jī)的烤燃特性研究以烤燃試驗(yàn)結(jié)合數(shù)值模擬為主，而針對(duì)尺寸較大的固體火箭發(fā)動(dòng)機(jī)，由于加熱體積大，烤燃溫場(chǎng)精確控制困難，危險(xiǎn)性大，因此實(shí)驗(yàn)研究報(bào)道較少。相對(duì)于小型固體火箭發(fā)動(dòng)機(jī)而言，大尺寸固體火箭發(fā)動(dòng)機(jī)在運(yùn)輸、貯存和使用過(guò)程中的安全隱患更大。為研究大尺寸固體火箭發(fā)動(dòng)機(jī)的熱安全性，本研究針對(duì)某固體火箭發(fā)動(dòng)機(jī)建立二維烤燃模型，分別計(jì)算了在不同升溫速率下固體推進(jìn)劑快速烤燃時(shí)的著火溫度、著火位置和延遲時(shí)間。

1 理論模型

1.1 物理模型

某固體火箭發(fā)動(dòng)機(jī)結(jié)構(gòu)示意圖[5]如圖1所示，由殼體、絕熱層、復(fù)合推進(jìn)劑、環(huán)氧樹脂擋板和噴管組成。

針對(duì)該固體火箭發(fā)動(dòng)機(jī)的實(shí)際尺寸建立二維軸對(duì)稱烤燃模型，采用如下簡(jiǎn)化假設(shè)：

(1)復(fù)合推進(jìn)劑的自熱反應(yīng)遵循與壓力相關(guān)的一階、二階Arrhenius定律；

(2)殼體與絕熱層以及絕熱層與推進(jìn)劑之間無(wú)接觸熱阻；

(3)AP/HTPB推進(jìn)劑假設(shè)為擬均質(zhì)、各向同性的致密材料，在整個(gè)模擬過(guò)程中均為固態(tài)，不考慮相變的影響；

(4)各材料的物性參數(shù)及化學(xué)動(dòng)力學(xué)參數(shù)為常量，不隨溫度變化；

(5)由于烤燃條件下，發(fā)動(dòng)機(jī)內(nèi)氣體流動(dòng)速度很小，忽略對(duì)流傳熱，僅考慮氣體和推進(jìn)劑間的熱傳導(dǎo)。

1.2 基本方程

針對(duì)AP/HTPB固體推進(jìn)劑的烤燃特性，建立基于兩步總包反應(yīng)[6]的烤燃模型，包括AP的熱分解反應(yīng)和最終放熱反應(yīng):

AP(Z)→Decomposition product

(1)

βZ+Binder(Y)→Final product

(2)

其中，反應(yīng)(1)和(2)的化學(xué)反應(yīng)速率R1和R2分別為:

R1=A1exp(-E1/RT)ρXp1.744

(3)

R2=A2exp(-E2/RT)ρYρZp1.750

(4)

(5)

(6)

(7)

式中：A1、A2為指前因子，s-1；E1、E2為反應(yīng)活化能， kJ/mol；R為通用氣體常數(shù)， J/(mol·K)-1；ρX、ρY、ρZ分別為AP、HTPB和AP分解產(chǎn)物的密度， kg/m3；p為壓強(qiáng)， Pa，按照理想狀態(tài)方程p=ρRT/M計(jì)算，M為摩爾質(zhì)量。Q1和Q2分別為反應(yīng)(1)和反應(yīng)(2)的反應(yīng)熱， kJ/kg；X、Y分別為AP和HTPB的質(zhì)量分?jǐn)?shù)，Z為AP分解產(chǎn)物的質(zhì)量分?jǐn)?shù)，X=ρX/ρ,Y=ρY/ρ,Z=ρZ/ρ；β為AP和HTPB的質(zhì)量當(dāng)量比。

固體火箭發(fā)動(dòng)機(jī)周圍壁面受熱，熱量向系統(tǒng)內(nèi)部傳遞，使AP/HTPB推進(jìn)劑溫度不斷升高，最終發(fā)生著火。殼體、絕熱層和固體推進(jìn)劑之間的熱傳遞、熱交換過(guò)程可以用如下非定常二維軸對(duì)稱方程描述：

(8)

式中：i=1,2,3分別表示殼體、絕熱層和固體推進(jìn)劑。ρi為密度， kg/m3；ci為比熱容， J/(kg·K)；λi為導(dǎo)熱率， W/(m·K)；qi為內(nèi)熱源，q1=q2=0，q3=R1Q1+R2Q2為AP/HTPB推進(jìn)劑的自熱反應(yīng)放熱率；T為溫度， K；t為時(shí)間， s。

1.3 邊界條件及初始條件

快速烤燃情況模擬的是外界環(huán)境失火對(duì)固體火箭發(fā)動(dòng)機(jī)的影響。本研究以升溫速率描述該情況下發(fā)動(dòng)機(jī)四周壁面的溫度邊界條件：

Ts=T0+kt,t>0

(9)

式中：t為時(shí)間,s；Ts為壁面溫度，K；T0為初始壁溫，K；k為升溫速率，K/s。

殼體、絕熱層、推進(jìn)劑任意兩種固體材料之間的交界面滿足溫度連續(xù)和熱流連續(xù)性條件：

Ta=Tb

(10)

(11)

式中：λa、λb、Ta、Tb分別為交界面的兩種固體材料的導(dǎo)熱系數(shù)和溫度。

殼體端面和噴管端面為絕熱邊界：

(12)

(13)

式中：λ1、λ2分別為殼體及噴管端面的導(dǎo)熱系數(shù)；T1、T2分別為殼體及噴管端面的溫度。

初始條件為：

T0=290K,X=0.88,Y=0.12,Z=0

(14)

2 數(shù)值模擬

2.1 數(shù)值驗(yàn)證

為驗(yàn)證烤燃模型的合理性，根據(jù)文獻(xiàn)[7]中的實(shí)驗(yàn)工況對(duì)AP/HTPB推進(jìn)劑的快速烤燃情況進(jìn)行數(shù)值模擬，并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較。實(shí)驗(yàn)裝置及試件結(jié)構(gòu)如圖2所示，樣品初始溫度為284K，鋼管壁面的平均升溫速率約為1.95K/s。根據(jù)實(shí)驗(yàn)工況和裝置尺寸進(jìn)行數(shù)值計(jì)算，計(jì)算所用化學(xué)反應(yīng)動(dòng)力學(xué)參數(shù)[8]如表1所示，計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)的對(duì)比如圖3所示。

1 AP/HTPB
Table 1 Kinetic parameters of the chemical reaction of AP/HTPB propellant

Reaction stepA/s-1E/(kJ·mol-1)Q/ (kJ·kg-1)1800126.99-297.021100163.919643.2

由圖3可看出，數(shù)值計(jì)算得到的著火溫度和著火延遲期為508K、192s，與試驗(yàn)數(shù)據(jù)505K、195s的誤差分別為0.6%和1.72%，數(shù)值模擬結(jié)果與實(shí)驗(yàn)測(cè)量結(jié)果吻合較好。

由此可見，兩步總包反應(yīng)模型能夠較好地反映AP/HTPB推進(jìn)劑的烤燃過(guò)程特性，可用于裝填A(yù)P/HTPB推進(jìn)劑火箭發(fā)動(dòng)機(jī)的快速烤燃特性數(shù)值預(yù)測(cè)。

2.2 計(jì)算方法與基本參數(shù)

采用基于單元格心有限體積法的FLUENT軟件進(jìn)行固體火箭發(fā)動(dòng)機(jī)快速烤燃數(shù)值計(jì)算，固體推進(jìn)劑自熱反應(yīng)和火箭發(fā)動(dòng)機(jī)溫度邊界條件通過(guò)用戶自定義函數(shù)(UDF)加載至軟件，分離式求解方法選用隱式算子分割算法(PISO)，密度、能量和組分方程對(duì)流項(xiàng)均采用二階迎風(fēng)格式離散?？紤]到發(fā)動(dòng)機(jī)可能面臨的火焰環(huán)境，一般火焰環(huán)境可由堆積廢物燃燒、航空油料燃燒、柴油、汽油燃燒等產(chǎn)生，其升溫速率在0.55 ～7.70K/s范圍內(nèi)[9]。為此，選取相關(guān)火焰環(huán)境中的升溫速率，分別以1.45、1.75、1.95、2.25 及2.45K/s的升溫速率對(duì)發(fā)動(dòng)機(jī)殼體進(jìn)行加熱，直至固體推進(jìn)劑著火。固體火箭發(fā)動(dòng)機(jī)尺寸如圖4所示。計(jì)算中分別對(duì)殼體內(nèi)壁點(diǎn)A(1170，159)、絕熱層中部點(diǎn)B(600，159)、AP/HTPB推進(jìn)劑肩部點(diǎn)C(890，150)和推進(jìn)劑內(nèi)壁點(diǎn)D(600,49)進(jìn)行監(jiān)測(cè)。固體火箭發(fā)動(dòng)機(jī)烤燃數(shù)值模型的材料物性參數(shù)如表2所示。

表2 材料物性參數(shù)Table 2 Physical property parameters of materials

2.3 網(wǎng)格獨(dú)立性

針對(duì)該發(fā)動(dòng)機(jī)結(jié)構(gòu)，采用四邊形結(jié)構(gòu)網(wǎng)格進(jìn)行網(wǎng)格劃分，并取3套網(wǎng)格Mesh 1、Mesh 2、Mesh 3進(jìn)行網(wǎng)格獨(dú)立性驗(yàn)證，網(wǎng)格數(shù)分別為326390、634900、1291890，對(duì)加熱速率為2.45K/s的烤燃工況進(jìn)行數(shù)值計(jì)算。圖5為固體火箭發(fā)動(dòng)機(jī)內(nèi)點(diǎn)(1100,163)的溫度-時(shí)間曲線，網(wǎng)格Mesh 3在t=60s時(shí)溫度為369.9K，而Mesh 1和Mesh 2的溫度分別為 374.5、369.0K,誤差分別為1.24%、0.24%。

圖6為530s時(shí)噴管喉部x=1154mm處的溫度分布曲線，網(wǎng)格Mesh 3在r=0mm處溫度為369.5K，而Mesh 1和Mesh 2溫度分別為 365、369K，誤差分別為1.22%、0.13%，可以發(fā)現(xiàn)Mesh 2的數(shù)值結(jié)果與網(wǎng)格加密一倍的數(shù)值結(jié)果一致，最終選取網(wǎng)格Mesh 2進(jìn)行數(shù)值計(jì)算。

3 結(jié)果與討論

3.1 著火過(guò)程分析

為分析固體火箭發(fā)動(dòng)機(jī)快速烤燃情況下AP/HTPB推進(jìn)劑的傳熱特性，以升溫速率為1.45K/s時(shí)為例，圖7是固體推進(jìn)劑著火前后的溫度云圖。

由圖7可看出，外界熱量通過(guò)熱傳導(dǎo)傳到系統(tǒng)內(nèi)部，絕熱層和AP/HTPB推進(jìn)劑的導(dǎo)熱系數(shù)遠(yuǎn)小于殼體，熱量在殼體中積累導(dǎo)致殼體與絕熱層和固體推進(jìn)劑存在溫差，且溫差隨時(shí)間推移而增大。694s時(shí)，推進(jìn)劑肩部出現(xiàn)高溫區(qū)域，如圖7(b)所示，該處推進(jìn)劑因自熱反應(yīng)釋放熱量不能及時(shí)擴(kuò)散最終發(fā)生著火，著火區(qū)域中心位置為(888.5,148.5)mm，軸向?qū)挾扰c徑向厚度為1.5mm。

圖8為升溫速率為1.45K/s時(shí)徑向剖面x=885mm處不同時(shí)刻的溫度曲線。

由圖8可知，推進(jìn)劑溫度始終低于外界溫度，且在升溫過(guò)程中僅推進(jìn)劑外部(徑向135～150mm)區(qū)域溫度隨時(shí)間升高，推進(jìn)劑內(nèi)部溫度幾乎不變，兩者溫差隨時(shí)間推移增大，說(shuō)明外界傳入的熱量和推進(jìn)劑自熱反應(yīng)釋放的熱量均堆聚在外壁，這是因?yàn)橥七M(jìn)劑熱擴(kuò)散率小、導(dǎo)熱率低及徑向尺寸大，所以熱量難以傳遞到推進(jìn)劑內(nèi)部。而且內(nèi)、外部區(qū)域的溫差將使得推進(jìn)劑烤燃反應(yīng)速率不一致。推進(jìn)劑外部區(qū)域溫度高，反應(yīng)速率大，放熱量多，升溫速度也快，而內(nèi)部區(qū)域則相反，這一過(guò)程又加劇了內(nèi)外部區(qū)域間的溫差。

圖9為升溫速率為1.45K/s時(shí)4個(gè)監(jiān)測(cè)點(diǎn)處溫度和組分質(zhì)量分?jǐn)?shù)隨時(shí)間的變化曲線。

由圖9(a)可知，熱量由殼體傳導(dǎo)至絕熱層，兩者溫度上升趨勢(shì)相同；推進(jìn)劑肩部溫度上升十分緩慢，當(dāng)溫度超過(guò)320K之后可看出升溫加快；推進(jìn)劑內(nèi)壁溫度變化小。圖9(b)為C點(diǎn)的溫度和組分含量變化曲線。由圖9(b)可知，約520s時(shí)C點(diǎn)溫度約400K時(shí)，可明顯觀察到組分AP(X)開始分解產(chǎn)生組分分解產(chǎn)物(Z)；520～650s時(shí)，AP分解吸收熱量且分解速率隨溫度升高而加快，使得溫升曲線斜率略有下降，但在外界導(dǎo)熱作用下溫度依然不斷升高；約650s時(shí)，觀察到組分HTPB(Y)開始下降，表明HTPB與分解產(chǎn)物間的反應(yīng)開始，該反應(yīng)釋放熱量大于AP分解吸收的熱量，使得推進(jìn)劑溫度升高；694s時(shí)溫度達(dá)到一定，AP已分解了一半，HTPB、分解產(chǎn)物反應(yīng)釋放的熱量遠(yuǎn)大于AP分解吸收的熱量，因此溫度急劇升高，發(fā)生著火。

3.2 升溫速率的影響

圖10為升溫速率分別為1.75、1.95、2.25及2.45K/s工況下著火時(shí)的溫度云圖。

由圖10可知，與1.45K/s升溫速率的工況一樣，著火點(diǎn)均出現(xiàn)在推進(jìn)劑肩部環(huán)狀區(qū)域，說(shuō)明升溫速率對(duì)著火位置影響很小。5種升溫速率對(duì)應(yīng)的著火溫度、著火延遲期及著火中心位置見表3。

表3 不同升溫速率下著火溫度、著火延遲期及著火中心位置Table 3 The ignition temperature, ignition delay and ignition centrality position at different heating rates

圖11為5種升溫速率下監(jiān)測(cè)點(diǎn)C的溫升曲線。

由圖11可知，當(dāng)溫度到達(dá)450K時(shí)，因AP受熱分解吸熱，曲線斜率降低；溫度超過(guò)480K時(shí)，因HTPB及分解產(chǎn)物發(fā)生反應(yīng)放熱，曲線斜率升高。且隨著升溫速率的增大，曲線斜率的變化越不明顯，越接近升溫速率。由于升溫速率增大，外界傳入推進(jìn)劑的熱量增加，使得推進(jìn)劑溫度上升加快，從而該處反應(yīng)速率隨溫度加大，因此著火延遲期隨著升溫速率的增大而縮短。

綜合數(shù)值預(yù)測(cè)結(jié)果可以發(fā)現(xiàn)，隨著升溫速率的增大，著火溫度與著火中心位置變化不顯著，著火時(shí)間明顯縮短。為進(jìn)一步分析升溫速率與著火延遲期之間的關(guān)系，根據(jù)表3繪制著火延遲期隨升溫速率改變的變化曲線，如圖12所示，發(fā)現(xiàn)著火延遲期ti與升溫速率k之間存在線性關(guān)系，擬合關(guān)系式為：ti=823.96481-115.47452k，相關(guān)系數(shù)R=0.989。

4 結(jié) 論

(1)針對(duì)裝填A(yù)P/HTPB推進(jìn)劑固體火箭發(fā)動(dòng)機(jī)快速烤燃過(guò)程建立烤燃模型，并進(jìn)行小尺度實(shí)驗(yàn)驗(yàn)證，數(shù)值模擬結(jié)果與試驗(yàn)數(shù)據(jù)吻合良好，證明所建模型是合理的。

(2)固體火箭發(fā)動(dòng)機(jī)在快速烤燃過(guò)程中，由于導(dǎo)熱情況和徑向尺寸影響，推進(jìn)劑內(nèi)、外壁溫差較大，導(dǎo)致化學(xué)反應(yīng)速率不一致，使溫差進(jìn)一步擴(kuò)大，并最終在推進(jìn)劑肩部發(fā)生著火。

(3)AP/HTPB推進(jìn)劑的著火均發(fā)生在推進(jìn)劑肩部的環(huán)形區(qū)域中，不同升溫速率對(duì)著火溫度和著火中心位置基本沒影響，著火溫度約為530K，但對(duì)著火延遲期有較大影響，著火延遲期與升溫速率存在線性關(guān)系，通過(guò)最小二乘法擬合出關(guān)系式：ti=823.96481-115.47452k,R=0.989。