葉 青，余永剛

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

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

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

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

1 理論模型

1.1 物理模型

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

(1)復合推進劑的自熱反應遵循與壓力相關的一階、二階Arrhenius 定律；

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

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

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

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

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

1.2 基本方程

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

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

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

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

式中：X、Y 分別為AP 和HTPB 的質量分數，Z 為AP 分解產物的質量分數，X=ρX/ρ，Y=ρY/ρ，Z=ρZ/ρ；β 為AP 和HTPB 的質量當量比。

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

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

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

連續(xù)型方程：

動量方程：

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

1.3 邊界條件及初始條件

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

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

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

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

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

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

2 數值結果與分析

2.1 數值模型驗證

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

表 1 AP/HTPB 推進劑化學反應動力學參數[16]Table 1Chemical reaction kinetic parameters ofAP/HTPB propellant[16]

圖 2 實驗裝置及試件結構示意圖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ā)動機慢速烤燃數值計算，固體推進劑自熱反應和火箭發(fā)動機溫度邊界條件通過用戶自定義函數(UDF)加載至軟件，并采用理想氣體方程模擬發(fā)動機內的自然對流。分離式求解方法選用隱式算子分割算法(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，推進劑內徑R2=49 mm，噴管長度L3=292 mm，推進劑外徑R3/R4=72/150。計算中分別對AP/HTPB 推進劑內點A (600 mm，100 mm)、絕熱層中部點B (600 mm，154 mm)推進劑肩部點C (885 mm，145 mm)、殼體內壁點D (1 100 mm，163 mm)及噴管喉部點E (1 270 mm，0 mm)進行監(jiān)測。固體火箭發(fā)動機烤燃數值模擬的材料物性[16-17]如表2 所示，ρ 為密度，cp為比熱容，λ 為導熱系數。

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

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

3 結果與討論

3.1 著火過程

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

表 2 材料物性參數Table 2 Parameters of materials

圖 6 不同網格下點C 處組分X 的質量分數隨時間變化的曲線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ā)動機殼體、絕熱層和推進劑溫度水平的不同，熱量是從徑向由外至內傳遞的。在重力的作用下（軸向方向），空腔內氣體因溫度差異形成自然對流，氣體流動速度較小，此時氣體流動速度最高為0.224 m/s，出現在軸線上400～600 mm 的位置。圖8(b)是30.60 h 時的溫度云圖，發(fā)現經過一段時間后，推進劑內部溫度升高但溫度分布依然不均勻。此時在推進劑的肩部出現一個高溫區(qū)（橢圓環(huán)標志處），溫度達到450 K，說明熱量在此處積累，推進劑反應速度加快。圖8(c)為推進劑發(fā)生著火時的溫度云圖，推進劑肩部高溫區(qū)域中心溫度不斷升高（圓環(huán)標志處），并最終發(fā)生著火如圖8(c)所示，在推進劑的1/2 軸對稱結構上有一橢圓形著火區(qū)域，該區(qū)域范圍為(868～880 mm，143～150 mm)，實際著火位置應為一圓環(huán)。著火延遲期為30.71 h，著火溫度為526.52 K，此時殼體壁面溫度為479.56 K。

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

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

3.2 升溫速率的影響

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

圖 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

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

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

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

4 結 論

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

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

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

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

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