葉青， 余永剛

(南京理工大學 能源與動力工程學院, 江蘇 南京 210094)

0 引言

固體推進劑通常用作太空運載火箭、航天器、導彈的推進燃料。目前面世的推進劑種類多樣，可滿足各種特定火箭發(fā)動機的應用要求。復合推進劑是由氧化劑、粘合劑、固化劑、金屬燃料、燃燒催化劑和其他組分組成的固體推進劑。高氯酸銨(AP)是復合推進劑中使用最廣泛的氧化劑，基于AP的固體推進劑具有高能量和優(yōu)異的燃燒特性。固體推進劑在常溫常壓環(huán)境下相對安定，一旦在受到攻擊或由于自身彈藥著火而引起燃燒或爆炸，將會對固體火箭發(fā)動機的安全存儲和使用造成巨大危害，因此固體火箭發(fā)動機的熱易損性問題日益引起重視?？救荚囼灪涂救紨?shù)值仿真是研究和評估炸藥和固體推進劑熱易損性的常用方法，對彈藥的設計、制造、運輸、存儲及使用具有重要指導意義。

迄今為止，許多研究人員對炸藥和固體推進劑的熱分解特性和烤燃特性進行了研究。Hanson-Parr等[1]測定了固體火箭推進劑的氧化劑和粘結劑的熱參數(shù)。Atwood等[2]緩慢加熱AP基復合推進劑進行烤燃試驗，同時測量推進劑內部溫度，結果發(fā)現(xiàn)，在加熱過程中的某個時刻，推進劑內的溫度開始升高并超過了爐溫，并且放熱反應在著火前加速直到發(fā)生反應或自燃[3]。在Tran[4]和Kim等[5]所做的不同推進劑一維熱爆炸試驗中也觀察到這一現(xiàn)象，試驗結果表明，點火發(fā)生的溫度在503～608 K之間，具體著火溫度取決于加熱速率和推進劑配方。Yoh等[6]利用不同厚度的金屬套筒和蓋板組合成試驗設備進行炸藥的熱烤燃試驗，試驗結果表明金屬套筒與蓋板的厚度越小則烤燃響應越不劇烈，這一結果為烤燃試驗的設計提供了參考。而為了更好地指導烤燃試驗，Victor[7]針對小尺度烤燃試驗建立一維計算模型，能得到某升溫速率或某火焰溫度下烤燃響應的時間與溫度，為試驗設計提供了依據(jù)。周捷等[8]針對熔鑄炸藥慢速烤燃過程進行慢速烤燃試驗，并結合數(shù)值模擬觀察炸藥內部溫度場變化，發(fā)現(xiàn)炸藥內部溫度場為同心類橢圓狀分布，且對流是影響炸藥點火點位置分布的主要因素。此后針對含能材料的小型烤燃試驗研究越來越多，Komai等[9]對縮水甘油疊氮聚醚(GAP)/高氯酸銨(AP)推進劑和HTPB/AP復合推進劑進行慢速烤燃試驗，發(fā)現(xiàn)GAP/AP推進劑的烤燃響應比HTPB/AP復合推進劑更溫和，后者烤燃裝置的破壞程度更嚴重。Ho[10]針對HTPB/AP和HTPB/RDX復合推進劑進行小規(guī)?？救紡椩囼?，研究推進劑的熱力學性質和加熱速率對烤燃過程的影響。也有學者將烤燃試驗與數(shù)值模型結合起來，Baer等[11]針對炸藥進行小尺度烤燃試驗，結果表明炸藥著火之前的熱/化學/機械耦合過程對其本身造成了熱損傷，并建立二維模型試圖重現(xiàn)該過程，模型將著火前炸藥釋放的能量近似為炸藥中心點的巨大熱量，而這并不適用于其他工況下的烤燃響應過程，尤其是裝藥形狀較為復雜的情況。

除了小型烤燃試驗，科研人員針對裝載推進劑的固體火箭發(fā)動機也進行了烤燃試驗研究。Weigand等[12]針對裝填AP/HTPB推進劑的固體火箭發(fā)動機進行快速烤燃和慢速烤燃試驗研究，試驗記錄了相應的著火時間與著火溫度，但無法確定著火發(fā)生位置。由于固體火箭發(fā)動機加熱體積大，烤燃溫場精確控制困難，危險性大，在此基礎上，烤燃特性研究由試驗研究深化為烤燃試驗與數(shù)值計算相結合。Cocchiaro[13]為幫助建立適用于大型固體火箭發(fā)動機系統(tǒng)的烤燃試驗和烤燃模型的分類標準，針對當時大型固體火箭發(fā)動機的快速烤燃試驗和烤燃數(shù)值模型進行了分析與總結。一般通過定義火焰環(huán)境溫度[14]或者升溫速率[15]來模擬烤燃試驗條件進行數(shù)值仿真，Sahin等[16]認為火焰環(huán)境下周圍環(huán)境的通風量也能影響烤燃響應。Yang等[17]針對某固體火箭發(fā)動機建立了二維烤燃簡化模型，開展了火焰環(huán)境下AP/HTPB推進劑的熱安全性研究。Li等[18-19]建立底排藥柱烤燃計算模型，以升溫速率模擬試驗條件并研究裝藥尺寸對底排藥烤燃響應特性的影響，發(fā)現(xiàn)裝藥內徑和長度對烤燃響應時間有影響。宋柳芳等[20]則在大、中、小型試驗件的慢速和快速烤燃試驗基礎上，通過數(shù)值模擬研究了烤燃試驗中試驗件的尺寸對烤燃結果的影響。實際上，針對烤燃試驗進行數(shù)值仿真計算，可以直觀地改變升溫速率、裝藥尺寸和約束等烤燃條件，預測熱反應過程并進行綜合分析。

綜上所述，目前對含能材料的熱安全性研究以烤燃試驗結合數(shù)值模擬為主，針對裝藥結構復雜的固體火箭發(fā)動機的熱安全性分析研究鮮有報道。本文以此為背景，針對較大尺寸的星型裝藥固體火箭發(fā)動機建立三維烤燃模型，數(shù)值分析其烤燃特性。首先根據(jù)試驗結果驗證數(shù)值模型的合理性，在此基礎上分別計算了在快速、中速和慢速加熱速率下固體火箭發(fā)動機烤燃時的著火溫度、延遲時間和著火區(qū)域的形狀、大小、中心位置。研究結果可為固體火箭發(fā)動機熱安全性分析提供參考。

1 計算模型

1.1 物理模型

某固體火箭發(fā)動機結構[21]如圖1所示，由發(fā)動機殼體、絕熱層、復合推進劑和環(huán)氧樹脂端板組成。圖1中r為徑向距離(m)，z為軸向距離(m)，O為柱坐標系Orφz原點，φ為周向角(rad)。本文針對該固體火箭發(fā)動機建立三維烤燃模型，采用如下假設：

1) AP/HTPB復合推進劑的熱分解反應遵循與壓力相關的1階、2階Arrhenius定律；

2) 殼體與絕熱層以及絕熱層與推進劑之間無接觸熱阻；

3) AP/HTPB復合推進劑假設為擬均質、各向同性的致密材料，在整個模擬過程中均為固態(tài)，不考慮相變的影響[22]；

4) 各材料的物性參數(shù)及化學動力學參數(shù)取為常量；

5) 烤燃條件下發(fā)動機內氣體流動緩慢，忽略對流傳熱，僅考慮氣體和推進劑間的熱傳導。

圖1 固體火箭發(fā)動機結構簡圖Fig.1 Schematic diagram of solid rocket motor

1.2 基本方程

針對AP/HTPB固體推進劑，AP的熱分解反應和最終放熱反應采用兩步總包反應[19]描述:

(1)

(2)

式中：β為AP與HTPB的質量當量比；R1和R2分別為反應1((1)式)和反應2((2)式)的化學反應速率，

R1=A1ρAPp1.744exp(-E1/8.314T)，

(3)

R2=A2ρHTPBρDp1.75exp(-E2/8.314T)，

(4)

A1、A2為指前因子(s-1)，ρAP、ρHTPB、ρD分別為AP、HTPB和AP分解產物的密度(kg/m3)，p為壓力(Pa)，按照理想狀態(tài)方程p=ρgRgTg計算，ρg為混合氣體密度，Rg為混合氣體常數(shù)，Tg為混合氣體溫度，E1、E2分別為反應1和反應2的活化能(kJ/mol)，T為溫度(K)。

推進劑中AP、HTPB和AP分解產物的組分方程如下：

(5)

(6)

(7)

式中：X、Y分別為AP和HTPB的質量分數(shù)，Z為AP分解產物的質量分數(shù)，X=ρAP/ρ,Y=ρHTPB/ρ,Z=ρD/ρ，ρ為推進劑密度；t為時間(s)。

固體火箭發(fā)動機殼體壁面受熱，熱量向系統(tǒng)內部傳遞。殼體與絕熱層以及絕熱層與推進劑之間無接觸熱阻。殼體、絕熱層、固體推進劑和氣體空腔之間的熱傳遞、熱交換過程可以用如下三維非定常方程描述：

(8)

式中：i=m,l,p,a分別表示殼體、絕熱劑、固體推進劑和氣體空腔；ρi為密度(kg/m3)；ci為比熱容(J/(kg·K)；Ti為溫度(K)；λi為導熱率(W/(m·K))；qi為內熱源(kJ)，qm=ql=qa=0 kJ，qp為AP/HTPB的自熱反應放熱率，qp=R1Q1+R2Q2，Q1和Q2分別為反應1和反應2的反應熱(kJ/kg)；r、φ、z分別為柱坐標的3個變量，r為徑向距離(m)，φ為周向角(rad)，z為軸向距離(m)。

本研究以升溫速率描述中速和快速烤燃工況下的溫度邊界條件

Ts=T0+kt,kt>0 K，

(9)

式中：Ts為壁面溫度(K)；T0為初始壁溫(K)；k為升溫速率(K/s)。

根據(jù)慢速烤燃試驗中受熱表面的溫度變化情況，以0.05 K/s先將固體火箭發(fā)動機殼體外表面升溫到400 K并保持8 h，之后以慢速升溫速率加熱火箭發(fā)動機殼體外表面，因此慢速烤燃工況下殼體溫度邊界以分段函數(shù)表示：

(10)

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

Tm=Tn,

(11)

(12)

式中:λm、λn、Tm、Tn分別為交界面兩種材料的導熱系數(shù)和溫度。

殼體端面緊鄰絕熱層，該處傳熱對發(fā)動機溫度分布影響極小，將其設為絕熱邊界：

(13)

噴管端面材料為環(huán)氧樹脂擋板，其導熱系數(shù)很小，可近似視為絕熱邊界：

(14)

初始條件為

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

(15)

1.3 數(shù)值驗證

在圓環(huán)柱狀AP/HTPB推進劑烤燃試驗[21]的基礎上進行數(shù)值模擬，并與試驗數(shù)據(jù)進行比較，驗證模型的合理性。試驗裝置結構如圖2所示，樣品初始溫度為301.53 K，金屬殼體壁面平均升溫速率為0.83 K/s. 根據(jù)試驗工況和裝置尺寸進行數(shù)值計算，計算所用化學反應動力學參數(shù)[5]如表1所示，計算結果與試驗數(shù)據(jù)的對比如圖3所示。數(shù)值計算得到的著火溫度和著火延遲期為533.7 K、379 s，與試驗數(shù)據(jù)533.5 K、370 s的誤差分別為0.6%和2.43%，數(shù)值模擬結果與試驗測量結果吻合較好，說明該模型可以較好反映AP/HTPB推進劑的烤燃響應過程。

圖2 圓環(huán)柱狀烤燃裝置結構簡圖Fig.2 Structure of experimental device

表1 AP/HTPB推進劑化學反應動力學參數(shù)Tab.1 Chemical reaction kinetic parameters of AP/HTPB propellant

圖3 監(jiān)測點計算結果與試驗結果對比圖Fig.3 Comparison of experimental and numerical results of monitoring points in cook-off test

1.4 計算方法及參數(shù)設置

固體火箭發(fā)動機烤燃數(shù)值模擬的材料物性[23]如表2所示。根據(jù)快速、中速、慢速烤燃試驗條件下含能材料表面的升溫速率范圍[24]，以快速升溫速率0.05～1.45 K/s、中速升溫速率0.005～0.011 K/s和慢速升溫速率2.4～3.3 K/h分別對發(fā)動機殼體進行加熱，直至固體推進劑著火。

表2 材料物性參數(shù)Tab.2 Parameters of materials

火箭發(fā)動機為六翼星型孔裝藥，如圖4所示，推進劑外徑為R，內徑0.33R，星型孔溝槽深度為0.6R，肉厚0.4R，推進劑長度5.47R，發(fā)動機總長度為9.37R，半徑為1.16R. 由于固體火箭發(fā)動機為軸對稱結構，采用1/12結構進行計算，計算中設置了監(jiān)測點以便于記錄各部分溫度變化和推進劑內組分變化情況，點坐標以柱坐標表示，分別為殼體內點a(0.29R，π/12 rad，4R)、絕熱層中部點b(0.267R，π/12 rad，4R)及推進劑外壁點c(2.257R，π/12 rad，4R)和中部點d(0.232R，π/12 rad，4R)進行監(jiān)測，如圖4所示，其中圓圈區(qū)域e為軸向著火位置范圍。為驗證網格劃分的合理性，以升溫速率1.45 K/s烤燃工況的數(shù)值模擬為例，對其進行無關性驗證。表3為相同升溫速率下，不同網格數(shù)量情況下數(shù)值模擬結果。由表3可知，1號、2號網格計算結果與3號網格計算結果誤差為3.6%與0.36%，因此在保證計算結果準確的前提下節(jié)省計算時間，選取2號網格進行數(shù)值計算。

圖4 固體火箭發(fā)動機尺寸及監(jiān)測點Fig.4 Solid rocket motor size and monitoring point

表3 網格無關性驗證Tab.3 Grids-independent verification results

2 烤燃數(shù)值結果與分析

2.1 快速烤燃

圖5 1.15 K/s升溫速率下各監(jiān)測點溫升曲線及 著火中心組分變化Fig.5 Temperature curve at each monitoring point and composition change in ignition center at 1.15 K/s heating rate

圖6 快速升溫速率下推進劑著火時刻剖面溫度云圖Fig.6 Temperature profile at propellant ignition time at fast heating rate

固體火箭發(fā)動機在快速烤燃情況下，殼體溫度勻速上升，在絕熱層的作用下熱量緩慢而持續(xù)地傳給推進劑，三者的溫差始終隨時間增大。以升溫速率1.15 K/s工況為例(見圖5)，推進劑體積較大，有利于外界傳入熱量在推進劑與絕熱層相交壁面堆聚，推進劑外壁面溫度曲線(c點和e區(qū)域中心點溫度)幾乎重合，如圖5(a)所示，推進劑內部和發(fā)動機空腔溫度幾乎不變，著火時推進劑內部溫度仍不超過300 K. 根據(jù)著火前100 s著火位置組分變化情況(見圖5(b))，此過程放熱反應速率極小，著火時刻AP(X)含量較高為0.6左右，達到著火溫度后AP(X)和分解產物(Z)含量均銳減，說明此時吸熱反應和放熱反應速率都急劇增大，由于放熱量遠大于吸熱量，著火位置溫度驟升，發(fā)生著火?？焖倏救记闆r下推進劑的初始著火出現(xiàn)在z=5.92R剖面上，此剖面位于推進劑的右側端面，緊鄰絕熱層。圖6為4種快速烤燃工況下推進劑著火時刻剖面溫度云圖，可以發(fā)現(xiàn)隨著升溫速率的增加，著火區(qū)域逐漸擴大，升溫速率1.15 K/s和1.45 K/s工況下著火區(qū)域均為環(huán)狀且著火溫度也比升溫速率0.55 K/s和0.85 K/s工況的高，原因為升溫速率越大，著火時刻殼體與絕熱層溫度越高，有利于著火區(qū)域熱量的積累，因此著火溫度較高。表4為4種快速升溫速率下的烤燃響應特征，可以發(fā)現(xiàn)隨著升溫速率的升高，著火延遲期明顯縮短，著火時刻殼體壁面與著火點溫度均有明顯升高，著火位置的軸向坐標(z/R)增加，即向推進劑右側端面移動。

表4 快速烤燃響應特征Tab.4 Response characteristics of fast cook-off

2.2 中速烤燃

圖7 0.011 K/s升溫速率下各監(jiān)測點溫升曲線及 著火中心組分變化Fig.7 Temperature curve at each monitoring point and composition change in ignition center at 0.011 K/s heating rate

在中速烤燃情況下，殼體與絕熱層的溫差與快速烤燃工況下的溫差相比較小，傳入推進劑的熱量堆聚在推進劑與絕熱層相交壁面，以升溫速率0.011 K/s為例，查看發(fā)動機內各點溫升曲線(見圖7(a))，發(fā)現(xiàn)在4 h之后推進劑外壁面溫升速率變小，結合圖7(b)，可發(fā)現(xiàn)4 h之后AP(X)分解速率加快，需要吸收部分熱量。推進劑著火前，AP經過一段時間分解含量已減至0.15，分解產物含量達到0.36，此時吸熱反應和放熱反應迅速發(fā)生，釋放大量熱量，發(fā)生著火。著火位置出現(xiàn)在z=5.92R的剖面上，如圖8所示，著火區(qū)域呈不連續(xù)點狀圓環(huán)分布，且與快速烤燃工況一樣，著火區(qū)域徑向剖面面積隨升溫速率增大而增加。推進劑在翼槽位置徑向厚度較小，熱量傳遞更快導致推進劑溫度上升速率和推進劑反應速率也更大，因此著火位置出現(xiàn)在圖示位置。表5為4種中速升溫速率下的烤燃響應特點，可以發(fā)現(xiàn)隨著升溫速率的升高，著火延遲期明顯縮短，著火時刻殼體壁面溫度升高，著火中心坐標的變化很微小。

圖8 中速升溫速率下推進劑著火時刻 剖面溫度云圖Fig.8 Temperature profile at propellant ignition time at medium heating rate

2.3 慢速烤燃

在慢速烤燃情況下，先以3 K/min的升溫速率將發(fā)動機殼體加熱至400 K，并維持一段時間，再對殼體升溫。以升溫速率3.3 K/h為例，查看發(fā)動機內各點溫升曲線(見圖9(a))發(fā)現(xiàn)：在整個過程中發(fā)動機殼體、絕熱層和推進劑外壁面的溫差較小，而推進劑內壁面溫度從保溫階段結束到推進劑發(fā)生著火整個過程溫度變化不超過20 K；著火位置推進劑經過21 h的分解，AP分解反應已基本完成，AP含量減至0.1，分解產物含量達到0.63，此時吸熱反應和放熱反應迅速發(fā)生，釋放大量熱量，發(fā)生著火，著火位置出現(xiàn)在z=5.84R的剖面上。如圖10所示，與中速烤燃工況相同，慢速烤燃工況的著火區(qū)域呈不連續(xù)點狀圓環(huán)分布，著火區(qū)域徑向剖面面積隨升溫速率增大而增加。表6為4種慢速升溫速率下的烤燃響應特點，與快速烤燃和中速烤燃兩種工況類似，隨著升溫速率的升高，著火延遲期明顯縮短，著火時刻殼體壁面溫度和推進劑著火溫度升高，著火位置向推進劑右側端面移動，著火中心徑向位置變化不明顯，對稱于翼槽中線呈兩點分布，面積隨升溫速率升高而擴大，有合二為一的趨勢。

表5 中速烤燃響應特征Tab.5 Response characteristics of medium cook-off

圖9 3.3 K/h升溫速率下各監(jiān)測點溫升曲線 及著火中心組分變化Fig.9 Temperature curve at each monitoring point and composition change in ignition center at 3.3 K/h heating rate

圖10 慢速升溫速率下推進劑著火時刻剖面溫度云圖Fig.10 Temperature profile at propellant ignition time at slow heating rate

表6 中速烤燃響應特征Tab.6 Response characteristics of slow cook-off

綜合3種烤燃工況的預測結果可以發(fā)現(xiàn)，隨著升溫速率的增大，著火溫度與著火中心位置有所變化，著火時間明顯縮短。為進一步分析升溫速率與著火溫度之間的關系，根據(jù)表4、表5、表6繪制著火溫度隨升溫速率改變的變化曲線，如圖11所示。由圖11可知：著火溫度Ti與升溫速率k之間存在二次函數(shù)關系，擬合關系式為Ti=516.659 36-1.267 8k+7.479 4k2，相關系數(shù)為0.95，這個擬合關系可以表現(xiàn)慢速和中速烤燃工況下(升溫速率較小)，著火溫度由推進劑自身成分決定，隨升溫速度變化微小；在快速烤燃工況下(升溫速度較大)，著火溫度受熱傳導情況影響較大，呈明顯上升趨勢。

圖11 著火溫度與升溫速率關系Fig.11 Ignition temperature versus heating rate

3 結論

本文針對固體火箭發(fā)動機在3種烤燃過程(快烤、中烤、慢烤)進行數(shù)值分析。得到以下結論：

1) 基于圓環(huán)柱狀AP/HTPB推進劑的烤燃試驗進行三維軸對稱烤燃數(shù)值模擬，并與試驗數(shù)據(jù)進行比較，二者數(shù)據(jù)結果吻合良好，證明所建模型是合理的。

2) 在升溫速率0.55～1.45 K/s快速烤燃情況下，殼體溫度上升迅速，外界傳入的熱量使與絕熱層相接觸的推進劑溫度上升較快，在熱傳導作用下推進劑到達著火溫度時，AP分解量不超過一半。著火位置緊鄰推進劑右側端面，著火區(qū)域徑向剖面面積隨升溫速率升高而增大，升溫速率越高越靠近端面，升溫速率為1.45 K/s時，著火區(qū)域形狀接近環(huán)狀。

3) 在升溫速率0.005～0.011 K/s中速和2.4～3.3 K/h慢速烤燃情況下，著火區(qū)域均呈不連續(xù)點狀圓環(huán)分布(見圖8和圖10)，著火位置處的AP組分在著火前已基本分解完全。慢速烤燃工況下推進劑著火點以翼槽中線呈對稱兩點分布，著火區(qū)域徑向剖面面積隨升溫速率升高而增大，兩點逐漸靠近，最終呈中速烤燃工況下的著火區(qū)域形狀。推進劑著火位置在慢速烤燃工況下隨升溫速率升高向右側端面移動，而中速烤燃工況下著火位置在軸向并沒發(fā)生變化。

4) 在升溫速率2.4 K/h～1.45 K/s范圍內，AP/HTPB推進劑著火溫度Ti與升溫速率k之間存在二次函數(shù)關系，擬合關系式為Ti=516.659 36-1.267 8k+7.479 4k2，相關系數(shù)R=0.95，這一規(guī)律適用于裝填AP/HTPB推進劑的固體火箭發(fā)動機。