葉 青，余永剛

（南京理工大學能源與動力工程學院，江蘇 南京 210094）

隨著火箭導(dǎo)彈等武器在戰(zhàn)場上的大量應(yīng)用和使用環(huán)境的日漸苛刻，固體火箭發(fā)動機的生存能力面臨更大的挑戰(zhàn)。導(dǎo)彈武器在受到敵方打擊或由于自身彈藥著火時而引起燃燒或爆炸，將會對工作人員和武器發(fā)射平臺造成巨大危害。因此，裝填復(fù)合推進劑的固體火箭發(fā)動機熱安全性問題引起了高度重視?？救紝嶒灪涂救紨?shù)值仿真計算是研究和評估彈藥和含能材料熱安全性的常用方法。

目前針對炸藥和推進劑的烤燃特性已進行了廣泛研究。學者們通過烤燃實驗對炸藥烤燃過程進行了研究，以裝藥尺寸、裝藥密度、裝藥孔隙率等烤燃條件為變量進行實驗，分析其對炸藥烤燃過程和響應(yīng)程度的影響[1-4]?？救紝嶒?zāi)軌蛑苯佑行У卦u價彈藥熱安全性，但成本高、周期長、測量數(shù)據(jù)不全且危險大，而針對烤燃實驗進行數(shù)值仿真計算，可以直觀地改變升溫速率、裝藥尺寸和約束等烤燃條件，預(yù)測熱反應(yīng)過程并進行綜合分析，其結(jié)果具有前瞻性。在此基礎(chǔ)上，烤燃特性研究由實驗研究深化為烤燃實驗與數(shù)值計算相結(jié)合。陳朗等[5]通過烤燃實驗和數(shù)值仿真研究炸藥熱反應(yīng)規(guī)律，提出了按照炸藥內(nèi)部熱量傳遞方向，把炸藥烤燃分為慢速、中速和快速烤燃。高峰等[6]利用自行設(shè)計的慢烤實驗結(jié)合數(shù)值模擬計算研究不同物理界面和不同界面厚度對黑索金(RDX)基高能炸藥烤燃特性的影響。牛余雷等[7]在不同升溫速率下，進行了不同尺寸的固黑鋁炸藥裝藥烤燃實驗和數(shù)值仿真計算，分析了裝藥尺寸對炸藥烤燃臨界環(huán)境溫度和響應(yīng)程度的影響。Aydemir 等[8]針對含能材料開發(fā)二維烤燃數(shù)值模型以預(yù)測彈藥的瞬態(tài)溫度分布、著火時間和著火位置。Ho[9]針對HTPB/AP 和HTPB/RDX 復(fù)合推進劑進行小規(guī)模烤燃彈實驗，研究推進劑的熱力學性質(zhì)和加熱速率對烤燃過程的影響。Komai 等[10]對縮水甘油疊氮聚醚(GAP)/高氯酸銨(AP)推進劑和HTPB/AP 復(fù)合推進劑進行慢速烤燃實驗，發(fā)現(xiàn)GAP/AP 推進劑的熱反應(yīng)比HTPB/AP 復(fù)合推進劑更溫和，后者烤燃裝置的破壞更嚴重。陳中娥等[11]利用同步差示掃描-熱重聯(lián)用儀、掃描電鏡和慢速烤燃實驗，對比分析了硝酸酯增塑聚醚(NEPE)推進劑和HTPB 推進劑的熱分解特性與慢速烤燃行為的關(guān)系。楊后文等[12]、Yang 等[13]針對某固體火箭發(fā)動機建立了二維烤燃簡化模型研究熱載荷作用下AP/HTPB 推進劑的熱安全性。Li 等[14]建立底排藥柱烤燃計算模型研究裝藥尺寸對底排藥烤燃響應(yīng)特性的影響，發(fā)現(xiàn)裝藥內(nèi)徑和長度對烤燃響應(yīng)時間有影響。

綜上所述，目前對含能材料的熱安全性研究以小型烤燃實驗結(jié)合數(shù)值模擬為主。由于固體火箭發(fā)動機加熱體積大，烤燃溫場精確控制困難，危險性大，特別是大尺寸固體火箭發(fā)動機烤燃實驗的難度更大，危險性也更大。目前針對大尺寸的固體火箭發(fā)動機的熱安全性分析研究報道較少。本文中以此為背景，擬對大尺寸固體火箭發(fā)動機的慢速烤燃特性進行數(shù)值分析。針對某固體火箭發(fā)動機建立二維烤燃模型，首先根據(jù)實驗結(jié)果驗證數(shù)值模型的合理性，在此基礎(chǔ)上分別計算在3 種慢速加熱速率下固體火箭發(fā)動機慢速烤燃時的著火溫度、著火位置和延遲時間。以期研究結(jié)果可為固體火箭發(fā)動機熱安全性分析提供參考。

1 理論模型

1.1 物理模型

某固體火箭發(fā)動機結(jié)構(gòu)[15]如圖1 所示，由殼體、絕熱層、復(fù)合推進劑、環(huán)氧樹脂擋板和噴管組成。本文中針對該固體火箭發(fā)動機的尺寸建立二維軸對稱烤燃模型，采用如下假設(shè)：

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

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

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

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

(5)烤燃條件下發(fā)動機內(nèi)存在氣體流動，假設(shè)氮氣為理想可壓氣體，密度隨溫度變化。

圖 1 固體火箭發(fā)動機結(jié)構(gòu)簡圖Fig. 1 Schematic drawing of solid rocket motor

1.2 基本方程

針對AP/HTPB 固體推進劑，AP 的熱分解反應(yīng)和最終放熱反應(yīng)采用兩步總包反應(yīng)[16]描述：

其中，β 為AP 和HTPB 的質(zhì)量當量比，反應(yīng)(1)和(2)的化學反應(yīng)速率R1和R2遵循與壓力相關(guān)的一階、二階Arrhenius 定律，分別為：

式中：A1、A2為指前因子；E1、E2為反應(yīng)活化能；R 為通用氣體常數(shù)；ρX、ρY、ρZ分別為AP、HTPB 和AP 分解產(chǎn)物Z 的密度；p 為壓力，按照理想狀態(tài)方程p=ρRT/M 計算，ρ 為密度，T 為溫度，M 為摩爾質(zhì)量。

推進劑中3 種組分，AP、HTPB、AP 分解產(chǎn)物（X、Y、Z）的組分方程如下：

式中：X、Y 分別為AP 和HTPB 的質(zhì)量分數(shù)，Z 為AP 分解產(chǎn)物的質(zhì)量分數(shù)，X=ρX/ρ，Y=ρY/ρ，Z=ρZ/ρ；β 為AP 和HTPB 的質(zhì)量當量比。

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

式中：ρi為密度，kg/m3；ci為比熱容，J/(kg·K)；T 為溫度，K；t 為時間，s；λi為導(dǎo)熱率，W/(m·K)；r、x 分別為徑向坐標和軸向坐標，i=1，2，3，4 分別表示殼體、絕熱劑、固體推進劑和氣體空腔；qi為內(nèi)熱源，q1=q2=q4=0，q3=R1Q1+R2Q2為AP/HTPB 的自熱反應(yīng)放熱率，Q1和Q2分別為反應(yīng)(1)和反應(yīng)(2)的反應(yīng)熱，kJ/kg。

烤燃條件下發(fā)動機豎直放置，噴口朝上，發(fā)動機內(nèi)存在氣體流動，假設(shè)氮氣為理想可壓氣體，密度隨溫度變化?？涨粌?nèi)氮氣的守恒型運動方程組如下。

連續(xù)型方程：

動量方程：

式中：u、v 為空腔內(nèi)氮氣的軸向、徑向流動速度，m/s；g 為重力加速度，m/s2；μ 為氮氣的黏度，kg/(m·s)。

1.3 邊界條件及初始條件

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

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

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

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

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

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

2 數(shù)值結(jié)果與分析

2.1 數(shù)值模型驗證

為驗證烤燃模型的正確性，根據(jù)文獻[10]中的實驗工況對AP/HTPB 推進劑的慢速烤燃情況進行數(shù)值模擬，并與實驗數(shù)據(jù)進行比較。AP/HTPB 推進劑的化學反應(yīng)動力學參數(shù)[13]如表1 所示, A 為指前因子，E 為活化能，Q 為反應(yīng)熱。實驗裝置及試件結(jié)構(gòu)如圖2 所示，鋼筒內(nèi)徑為25 mm，兩端用鋁制剪切板密封。實驗中初始溫度為428 K，采用3.3 K/h 的升溫速率加熱試件，并在推進劑中心安裝熱電偶，實驗記錄了1 000～1 350 min 期間溫度變化曲線，著火時樣品溫度為504.55 K。根據(jù)實驗工況進行數(shù)值計算，計算結(jié)果如圖3 所示，著火溫度為503.43 K，著火延遲期為1 350.5 min，著火溫度誤差為0.22%，著火位置有兩處，位于推進劑的軸線上，數(shù)值計算得到的推進劑中心溫度變化曲線與實驗記錄結(jié)果的對比如圖4 所示，橫坐標表示對烤燃裝置的加熱時間，縱坐標為推進劑中心的溫度，可以發(fā)現(xiàn)數(shù)值模擬結(jié)果與實驗測量結(jié)果吻合較好。由此可見，兩步總包反應(yīng)模型能夠較好反映AP/HTPB 推進劑的烤燃特性，可用于裝填A(yù)P/HTPB 推進劑火箭發(fā)動機的烤燃特性數(shù)值預(yù)測。

表 1 AP/HTPB 推進劑化學反應(yīng)動力學參數(shù)[16]Table 1Chemical reaction kinetic parameters ofAP/HTPB propellant[16]

圖 2 實驗裝置及試件結(jié)構(gòu)示意圖Fig. 2 Sketch map of experimental device andspecimen structure

圖 3 實驗裝置著火時刻t=1350.5 min 溫度云圖Fig. 3 Temperature distribution at the ignition timeof the test device (t=1 350.5 min)

圖 4 推進劑中心溫度隨對烤燃裝置加熱時間的變化Fig. 4 The temperature in the center of the propellant varyingwith the heating time of the cook-off device

2.2 計算模型及計算方法

采用基于單元格心有限體積法的FLUENT 軟件進行固體火箭發(fā)動機慢速烤燃數(shù)值計算，固體推進劑自熱反應(yīng)和火箭發(fā)動機溫度邊界條件通過用戶自定義函數(shù)(UDF)加載至軟件，并采用理想氣體方程模擬發(fā)動機內(nèi)的自然對流。分離式求解方法選用隱式算子分割算法(PISO)，壓力插值格式選用PRESTO！方法，密度、能量和組分方程對流項均采用二階迎風格式離散。以0.05 K/s 將邊界升溫到400 K 并保持8 h 后分別以3.6、7.2 及10.8 K/h 加熱速率對發(fā)動機殼體進行加熱，直至固體推進劑著火。固體火箭發(fā)動機尺寸如圖5 所示，發(fā)動機總長L1=1 400 mm，發(fā)動機外徑R1=348 mm，噴管直徑R5=127.5 mm，推進劑長度L2=820 mm，推進劑內(nèi)徑R2=49 mm，噴管長度L3=292 mm，推進劑外徑R3/R4=72/150。計算中分別對AP/HTPB 推進劑內(nèi)點A (600 mm，100 mm)、絕熱層中部點B (600 mm，154 mm)推進劑肩部點C (885 mm，145 mm)、殼體內(nèi)壁點D (1 100 mm，163 mm)及噴管喉部點E (1 270 mm，0 mm)進行監(jiān)測。固體火箭發(fā)動機烤燃數(shù)值模擬的材料物性[16-17]如表2 所示，ρ 為密度，cp為比熱容，λ 為導(dǎo)熱系數(shù)。

圖 5 固體火箭發(fā)動機尺寸及監(jiān)測點位置Fig. 5 Sizes of the solid rocket motor and locations of monitoring points

針對該固體火箭發(fā)動機軸對稱結(jié)構(gòu)，采用1/2結(jié)構(gòu)模型。采用四邊形非結(jié)構(gòu)網(wǎng)格進行網(wǎng)格劃分，并取3 套網(wǎng)格Mesh 1、Mesh 2、Mesh 3 進行網(wǎng)格獨立性驗證，網(wǎng)格數(shù)分別為41 507、112 559、199 424，對升溫速率為10.8 K/h 的烤燃工況進行數(shù)值分析。圖6 為監(jiān)測點C (885 mm，145 mm)組分X 的質(zhì)量分數(shù)隨時間變化的曲線，網(wǎng)格Mesh 3 在32 000 s 時，組分X 的質(zhì)量分數(shù)為0.167 8，Mesh 1 和Mesh 2 的X 組分的質(zhì)量分數(shù)分別為0.174 03、0.169 3，與Mesh 3 的誤差分別為3.71%、0.89%，可以發(fā)現(xiàn)Mesh 2 的數(shù)值結(jié)果與網(wǎng)格加密后的數(shù)值結(jié)果一致，最終選取網(wǎng)格Mesh 2 進行數(shù)值計算。

3 結(jié)果與討論

3.1 著火過程

圖7～8 是加熱速率為3.6 K/h 工況下有無考慮發(fā)動機空腔內(nèi)自然對流的數(shù)值計算結(jié)果。對比圖7～8 可以發(fā)現(xiàn)，在22.22 h，是否考慮自然對流對推進劑內(nèi)的溫度分布影響并不明顯，對空腔內(nèi)的溫度分布影響顯著。同一時刻未考慮自然對流的空腔溫度明顯低于考慮自然對流情況，這表示未考慮自然對流情況下推進劑內(nèi)部熱傳導(dǎo)方向依然為徑向由外至內(nèi)，而考慮自然對流的熱傳導(dǎo)方向是推進劑周邊向推進劑內(nèi)傳遞。觀察30.60 h 時推進劑內(nèi)壁面低溫區(qū)域分布(溫度低于415 K)，發(fā)現(xiàn)未考慮空腔自然對流的低溫區(qū)域面積大于考慮自然對流的，且推進劑肩部沒有出現(xiàn)高溫區(qū)，最終發(fā)生著火的時刻為30.96 h、著火溫度為540.44 K，而考慮自然對流的著火時刻為30.71 h，著火溫度526.52 K，分別相差0.81%和2.64%。對于本文這種尺寸較大的固體火箭發(fā)動機而言，空腔內(nèi)的自然對流通過影響慢烤過程中空腔內(nèi)的溫度分布，進而影響推進劑溫度分布并最終導(dǎo)致著火延遲期的改變，因此空腔內(nèi)自然對流在熱安全性精確分析中不可忽略。

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

圖 6 不同網(wǎng)格下點C 處組分X 的質(zhì)量分數(shù)隨時間變化的曲線Fig. 6 Variation of mass fraction of componentX at point C in different grids with time

圖 7 不考慮自然對流不同時刻溫度云圖Fig. 7 Temperature distribution at different times without natural convection

圖 8 考慮自然對流不同時刻溫度云圖Fig. 8 Temperature distribution at different times with natural convection

分析固體火箭發(fā)動機慢速烤燃情況下AP/HTPB 推進劑的著火過程，以升溫速率為3.6 K/h 為例，圖8(a)是22.22 h 時的溫度云圖和流線圖，溫度云圖清晰地反映了發(fā)動機殼體、絕熱層和推進劑溫度水平的不同，熱量是從徑向由外至內(nèi)傳遞的。在重力的作用下（軸向方向），空腔內(nèi)氣體因溫度差異形成自然對流，氣體流動速度較小，此時氣體流動速度最高為0.224 m/s，出現(xiàn)在軸線上400～600 mm 的位置。圖8(b)是30.60 h 時的溫度云圖，發(fā)現(xiàn)經(jīng)過一段時間后，推進劑內(nèi)部溫度升高但溫度分布依然不均勻。此時在推進劑的肩部出現(xiàn)一個高溫區(qū)（橢圓環(huán)標志處），溫度達到450 K，說明熱量在此處積累，推進劑反應(yīng)速度加快。圖8(c)為推進劑發(fā)生著火時的溫度云圖，推進劑肩部高溫區(qū)域中心溫度不斷升高（圓環(huán)標志處），并最終發(fā)生著火如圖8(c)所示，在推進劑的1/2 軸對稱結(jié)構(gòu)上有一橢圓形著火區(qū)域，該區(qū)域范圍為(868～880 mm，143～150 mm)，實際著火位置應(yīng)為一圓環(huán)。著火延遲期為30.71 h，著火溫度為526.52 K，此時殼體壁面溫度為479.56 K。

上文設(shè)置的5 個監(jiān)測點是溫度變化明顯且具有代表性的位置，溫度監(jiān)測結(jié)果如圖9(a)所示，殼體(D)溫度和噴管喉部處(E)的溫度升溫趨勢基本一致，而絕熱層內(nèi)(B)、推進劑肩部(C)和推進劑內(nèi)部(A)處的升溫速率在約17.2 h 時出現(xiàn)了明顯的下降，而在圖9(b)中監(jiān)測點C 處的推進劑在此時反應(yīng)速率加快。原因是17.2 h 之前推進劑中的溫度變化是由外向內(nèi)的熱傳導(dǎo)主導(dǎo)的，推進劑導(dǎo)熱能力弱，導(dǎo)致推進劑內(nèi)部溫度分布不均，熱量在推進劑外表面堆聚。在此之后推進劑外壁處AP 的分解反應(yīng)和最終產(chǎn)物的生成反應(yīng)速率加快，需要吸收一定的熱量，導(dǎo)致該處推進劑溫升變慢，而推進劑內(nèi)部和周圍絕熱層受其影響升溫速率下降。當溫度上升到著火點附近時，推進劑內(nèi)AP 分解反應(yīng)基本完成，最終產(chǎn)物生成反應(yīng)釋放大量熱量，使得推進劑及周圍的絕熱層溫度上升加速，推進劑肩部的溫度變化最為明顯，最終在肩部位置著火。

圖 9 不同監(jiān)測點的溫度和監(jiān)測點C 處組分的質(zhì)量分數(shù)的變化曲線Fig. 9 Temperature curves of different monitoring points and mass fraction curves of components at point C

3.2 升溫速率的影響

針對3 種不同加熱速率工況進行固體火箭發(fā)動機慢速烤燃數(shù)值計算，加熱速率分別為3.6、7.2 和10.8 K/h。圖10、11 為7.2 和10.8 K/h 加熱速率下著火時刻的溫度云圖，由圖8(c)、10 和11 可發(fā)現(xiàn)，不同升溫速率下著火位置均出現(xiàn)在推進劑的肩部，原因是推進劑徑向尺寸較大，烤燃過程中熱量不斷向推進劑內(nèi)部傳遞，而推進劑肩部被絕熱層呈直角包圍，兩者導(dǎo)熱率的差異使得熱量更易在此處積累。3 種工況的著火位置局部放大圖顯示隨著升溫速率的提高，著火時刻殼體溫度明顯升高，著火位置向推進劑與絕熱層交界處移動，且著火區(qū)域的二維截面由橢圓形變?yōu)榘霗E圓形。這是因為升溫速率的提高使得絕熱層溫升加快，導(dǎo)致推進劑內(nèi)的高溫點位置更靠近絕熱層，且推進劑即將著火時釋放熱量來不及向絕熱層傳遞，因此著火位置向兩者交界處移動。

圖 10 在7.2 K/h 的加熱速率下，著火時刻(t=79 409.5 s)的溫度云圖Fig. 10 Temperature distribution of propellant at ignition time (t=79 409.5 s) under the heating rate of 7.2 K/h

圖 11 在10.8 K/h 的加熱速率下，著火時刻(t=67 257.5 s)的溫度云圖Fig. 11 Temperature distribution of propellant at ignition time (t=67 257.5 s) under the heating rate of 10.8 K/h

根據(jù)計算結(jié)果獲得了不同加熱速率下的著火特征參數(shù)，如表3 所示。以0.05 K/s 將殼體外表面快速升溫到400 K，并保持8 h 的情況下，慢速加熱速率在3.6～10.8 K/h 范圍內(nèi)，著火延遲期td隨著加熱速率k 增加而縮短，兩者滿足非線性關(guān)系td=65.666 38/(65.666 38+0.357 83 k-0.011 68 k2)，如圖12 所示。著火溫度和殼體溫度均隨著升溫速率增大而升高。隨著加熱速率的增大，著火中心逐漸向推進劑與絕熱層的交界處移動，且著火區(qū)域面積縮小。

表 3 不同加熱速率下的著火特征參數(shù)Table 3 Ignition characteristic parameters at different heating rates

圖 12 著火延遲期隨升溫速率的變化Fig. 12 Ignition delay vaying with heating rate

4 結(jié) 論

針對較大尺寸的固體火箭發(fā)動機在3 種不同升溫速率慢速烤燃過程進行數(shù)值分析，得到以下結(jié)論：

（1）基于AP/HTPB 推進劑的慢速烤燃特性建立二維軸對稱烤燃模型，針對小尺度慢速烤燃實驗進行數(shù)值模擬驗證，兩者數(shù)據(jù)結(jié)果吻合良好，證明所建模型是合理的；

（2）建立考慮固體火箭發(fā)動機內(nèi)空腔自然對流的慢速烤燃模型，并分析有無考慮空腔內(nèi)自然對流慢速烤燃數(shù)值結(jié)果。有、無考慮發(fā)動機內(nèi)空腔自然對流情況下發(fā)動機內(nèi)熱傳導(dǎo)方向不同，發(fā)動機溫度分布差異明顯，且著火時刻和著火溫度分別相差0.81%、2.64%，在固體火箭發(fā)動機的熱安全性精確分析中不可忽略；

（3）固體火箭發(fā)動機慢速烤燃過程中，前期主導(dǎo)推進劑溫度變化的是熱傳導(dǎo)，后期則是由推進劑的自熱反應(yīng)主導(dǎo)，當推進劑溫度上升到著火點附近時，AP 已經(jīng)基本分解完全，此時最終產(chǎn)物生成反應(yīng)釋放大量熱量，使得推進劑溫度迅速上升，并在肩部位置著火；

（4）3.6、7.2 和10.8 K/h 加熱速率下，固體火箭發(fā)動機的著火點均出現(xiàn)在推進劑的肩部，但隨著加熱速率的升高，著火位置向推進劑與絕熱層交界處移動，且著火區(qū)域的二維截面由橢圓形變?yōu)榘霗E圓形。3 種加熱速率對應(yīng)的著火延遲期、著火溫度及著火時殼體溫度分別為30.71、22.06、18.68 h，526.52、528.10、530.64 K，和479.56、496.82、508.77 K。