楊后文， 余永剛， 葉 銳

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

1 引 言

高氯酸銨(AP)/端羥基聚丁二烯(HTPB)推進劑由于其燃燒的穩(wěn)定性和燃燒速率的可控性，廣泛用于固體火箭發(fā)動機中。但在運輸、儲存、使用過程中，由于環(huán)境溫度變化或周圍熱源的影響，可能引起推進劑內(nèi)部溫度升高，發(fā)生失控的化學(xué)反應(yīng)和能量釋放，造成嚴重的災(zāi)難性事故。近年，含能材料的熱安全性問題日益引起了人們的高度關(guān)注，國內(nèi)外學(xué)者針對炸藥和推進劑的烤燃特性進行了廣泛研究。馮曉軍等[1]利用自行研制的烤燃實驗裝置，研究了JB-B、TNT和R852三種炸藥尺寸對慢速烤燃響應(yīng)特性的影響。Howard等[2]對LX-10炸藥慢速烤燃(2 ℃·h-1)過程進行了數(shù)值模擬，并在計算中考慮了炸藥和殼體間空氣層對烤燃過程的影響。王沛等[3]研究了四種邊界條件下固黑鋁炸藥(GHL)的熱反應(yīng)規(guī)律，結(jié)果表明，慢速烤燃下炸藥點火時的環(huán)境溫度比快速烤燃低。向梅等[4]利用有限元程序LS-DYNA3D對不同升溫速率下復(fù)合炸藥的烤燃過程進行了數(shù)值模擬，他們認為升溫速率較小時，復(fù)合藥柱的熱安全性取決于內(nèi)部高能炸藥的特性。關(guān)于AP/HTPB推進劑的慢烤燃行為已有大量的實驗研究。Ferschl等[5]利用超小型烤燃彈(SSCB)裝置研究了AP/HTPB推進劑的烤燃特性與熱化學(xué)性能之間的關(guān)系。他們認為推進劑烤燃響應(yīng)的劇烈程度與它的高爆熱性質(zhì)有關(guān)，當提高推進劑中增塑劑的含量時，AP/HTPB推進劑的烤燃響應(yīng)程度減弱。Komai等[6]利用差示掃描量熱分析儀(DSC)測量了AP/HTPB推進劑在不同升溫速率下的著火溫度。結(jié)果表明，著火溫度隨著升溫速率的增大而提高。Rodrigo等[7]利用慢烤裝置(SCO)重點研究了在6 ℃·h-1升溫速率下，AP/HTPB推進劑的慢烤燃特性。陳中娥等[8-9]則聯(lián)合使用差示掃描-熱重聯(lián)用儀(DSC-TG)、掃描電鏡(SEM)和慢烤試驗，研究了AP/HTPB推進劑的熱分解特性與慢速烤燃行為的關(guān)系。結(jié)果表明，AP低溫分解形成的大量孔隙，是導(dǎo)致AP/HTPB推進劑慢速烤燃響應(yīng)劇烈的重要原因，并在此結(jié)論基礎(chǔ)上提出了改善AP/HTPB推進劑慢速烤燃響應(yīng)的技術(shù)途徑。廖林泉等[10]通過真空安全性測試，研究了3種HTPB推進劑配方在75℃下的熱安全性，結(jié)果表明，AP/HTPB推進劑的熱安全性與AP粒度有關(guān)，粒度越小，發(fā)生分解反應(yīng)的活性越高。趙孝彬等[11]則采用慢速烤燃試驗裝置結(jié)合熱電偶測溫及傳感器測壓技術(shù)，研究了配方組成、燃速、升溫速率、約束條件、自由體積等對HTPB推進劑慢速烤燃特性的影響。綜上所述，目前國內(nèi)外對含能材料的烤燃計算主要以炸藥為主，而針對AP/HTPB復(fù)合固體推進劑慢烤燃特性的研究主要集中在實驗分析上，相應(yīng)的數(shù)值模擬研究還很少。

本研究針對某固體火箭發(fā)動機建立了一種二維簡化烤燃模型，采用基于Beckstead-Derr-Price(BDP)多火焰結(jié)構(gòu)[12]的兩步總包反應(yīng)機理來描述AP/HTPB推進劑的慢烤燃過程，分別對AP/HTPB推進劑在升溫速率為1.8，3.6，7.2 K·h-1下的慢烤燃行為進行數(shù)值模擬，得到了AP/HTPB推進劑最初的著火位置、著火延遲期以及著火溫度。研究結(jié)果對分析相關(guān)固體推進劑的烤燃特性和熱安全性具有一定的參考價值。

2 理論模型

2.1 基本假設(shè)

固體火箭發(fā)動機一般由殼體、絕熱層和推進劑裝藥三部分組成，本研究根據(jù)某發(fā)動機尺寸建立一種二維簡化烤燃模型，并采用如下簡化假設(shè)：

(1) AP/HTPB推進劑為均質(zhì)、各向同性的致密材料;

(2) 推進劑在整個模擬過程中為固態(tài)，不考慮相變影響;

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

(4) 推進劑的自熱反應(yīng)遵循與壓力相關(guān)的一階、二階Arrhenius定律。

2.2 基本方程

當發(fā)動機四周壁面受到熱作用時，溫度逐漸升高，熱量向系統(tǒng)內(nèi)部傳遞，使推進劑溫度不斷升高，最終發(fā)生著火。整個過程中，系統(tǒng)內(nèi)部熱傳遞、熱交換和熱分解過程可用如下方程[13]描述：

(1)

式中，i=1,2,3，分別表示殼體、絕熱層和固體推進劑。ρi、ci、λi、qi分別為對應(yīng)材料的密度, kg·m-3; 比熱容, J·kg-1·K-1; 導(dǎo)熱率, W·m-1·K-1和內(nèi)熱源，q1=q2=0,q3為固體推進劑自熱反應(yīng)放熱率。針對AP/HTPB推進劑的烤燃特性，將文獻[14-16]中提出的兩步總包反應(yīng)引入到烤燃理論分析中，包括AP分解吸熱反應(yīng)和粘合劑與AP分解產(chǎn)物的最終放熱反應(yīng)：

AP(X)→decomposition product(Z)

(A)

βZ+binder(Y)→final product

(B)

式中，反應(yīng)(A)和(B)的化學(xué)反應(yīng)速率R1和R2分別為

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

(2)

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

(3)

式中，A1、A2為指前因子，s-1;E1、E2為反應(yīng)活化能,kJ·mol-1;R為通用氣體常數(shù), J·mol-1·K-1;ρX、ρY、ρZ分別為AP、HTPB和AP分解產(chǎn)物Z的密度，kg·m-3;p為壓力，Pa，按照理想狀態(tài)方程pV=nRT計算。根據(jù)以上兩步反應(yīng)，AP/HTPB推進劑自熱反應(yīng)的能量方程和組分方程[17]如下：

(4)

(5)

(6)

(7)

式中，ρ為AP/HTPB推進劑的密度，kg·m-3;Q1和Q2分別為反應(yīng)(A)和反應(yīng)(B)的反應(yīng)熱, kJ·kg-1;m1、m2分別為AP和HTPB的質(zhì)量分數(shù)，m3為AP分解產(chǎn)物的質(zhì)量分數(shù)，m1=ρX/ρ，m2=ρY/ρ，m3=ρZ/ρ;β為AP和HTPB的質(zhì)量當量比(β=88/12)。

2.3 邊界條件與初始條件

四周壁面為溫度邊界條件：

t>0,Ts=T0+kt

(8)

式中，t為時間, s;Ts為壁面溫度, K;T0為初始壁溫, K;k為升溫速率, K·n-1。

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

Ta=Tb

(9)

λaTa=λbTb

(10)

式中，Ta、Tb與λa、λb分別為相接觸兩種材料的溫度與導(dǎo)熱系數(shù)。

初始條件為：

T0=300 K,m1=0.88,m2=0.12,m3=0

(11)

3 計算結(jié)果與分析

3.1 計算方法及基本參數(shù)

采用基于有限體積法的計算流體力學(xué)軟件FLUENT[18]對固體發(fā)動機的慢烤燃模型進行數(shù)值模擬。AP/HTPB推進劑的自熱反應(yīng)和邊界條件通過用戶自定義標量(UDS)和自定義函數(shù)(UDF)加載到FLUENT中進行計算。采用PISO算法[18]，即壓力的隱式算子分割算法。密度、能量和組分方程的離散采用二階迎風(fēng)格式。計算所需的物性參數(shù)與動力學(xué)參數(shù)[17,19-20]見表1和表2。

表1材料的物性參數(shù)

Table1Physical property parameters of materials

materialρ/kg·m-3cp/J·kg-1·K-1λ/W·m-1·K-1shell785046242heat?insulatinglayer14508400．168epoxyresinboard180012000．15AP/HTPBpropellant182612550．389

Note:ρis density,cpis specific heat capacity,λis thermal conductivity.

表2AP/HTPB化學(xué)動力學(xué)參數(shù)

Table2Chemical reaction kinetic parameters of AP/HTPB

reactionstepA/s-1E/kJ·mol-1Q/kJ·kg-1first?step800137．18-297second?step1100178．759643．2

Note:Ais pre-exponential constant,Eis activation energy,Qis reaction heat.

3.2 計算結(jié)果與實驗數(shù)據(jù)的對比

為了驗證本文所建模型的正確性，根據(jù)文獻[6]中對AP/HTPB推進劑的熱烤裝置進行建模，得出不同升溫速率k下的著火溫度Tc，并與實驗數(shù)據(jù)[6]相比較，如圖1所示。由圖1可知，數(shù)值模擬結(jié)果與實驗測量結(jié)果基本吻合，證明本研究采用的兩步反應(yīng)機理模型能夠較好地反映烤燃過程特性。

圖1logk與Tc的關(guān)系

Fig.1The relation of logkandTc

圖2AP/HTPB推進劑內(nèi)部的溫度時程曲線

Fig.2History curves of internal temperature of AP/HTPB propellant

另外，本研究還針對文獻[7,17]中的慢烤試驗裝置建立模型，模擬在升溫速率為6 K·h-1下AP/HTPB推進劑的慢烤燃過程，計算推進劑內(nèi)部溫度隨時間的變化，將計算結(jié)果與實驗結(jié)果[7,17]進行對比，如圖2所示。數(shù)值計算得到的著火溫度為520K，著火延遲期為52426 s，而實驗測量得到的著火溫度為532 K，著火延遲期為54420 s，著火溫度與著火延遲期的計算誤差分別為2.26%和3.66%。由此可見，本研究所采用的模型是合理可行的，可用于發(fā)動機慢烤燃特性的數(shù)值預(yù)測。

3.3 固體火箭發(fā)動機慢烤燃特性的數(shù)值預(yù)測

圖3為某模擬發(fā)動機[21]結(jié)構(gòu)簡圖，外徑Φ1=60 mm，殼體厚度δ1=2 mm，推進劑裝藥外徑Φ2=54 mm，內(nèi)徑Φ3=8 mm，長450 mm，殼體與推進劑之間有1 mm厚絕熱層，發(fā)動機內(nèi)密封氣體假設(shè)為空氣。不考慮噴管部分對計算的影響，由于藥柱中心的空氣腔直徑很小，僅考慮導(dǎo)熱效應(yīng)。燃燒室出口用5 mm厚的環(huán)氧樹脂板密封。四周壁面為溫度邊界條件，初始溫度為300 K。計算時，首先以0.05 K·s-1的升溫速率將壁面溫度升至450 K，然后保持該溫度一段時間，再分別以1.8，3.6，7.2 K·h-1的升溫速率加熱壁面，直至發(fā)生著火反應(yīng)。

圖3固體火箭發(fā)動機結(jié)構(gòu)簡圖

Fig.3Schematic drawing of solid rocket motor

圖4為升溫速率為1.8 K·h-1條件下，發(fā)動機橫向截面在不同時刻的溫度分布云圖?？梢姡谕饨鐭嵩醋饔孟?，由于殼體的導(dǎo)熱系數(shù)與絕熱層和推進劑的導(dǎo)熱系數(shù)相差較大，壁面溫度向內(nèi)部傳導(dǎo)緩慢，4500 s時，殼體溫度高于AP/HTPB推進劑內(nèi)部溫度，系統(tǒng)最大溫差接近40 K。56620 s時，由于升溫速率較慢，外壁熱量有足夠的時間傳向內(nèi)部，系統(tǒng)整體溫差減小。隨著AP/HTPB推進劑溫度的升高，推進劑緩慢的自熱反應(yīng)速率逐漸增加，產(chǎn)生的熱量來不及向周圍釋放，使得內(nèi)部藥柱溫度升高較快。102210 s時，AP/HTPB推進劑內(nèi)部溫度已經(jīng)高于殼體溫度。113335 s時，推進劑內(nèi)壁面中心區(qū)域附近溫度最高，即達到著火溫度，最初著火位置發(fā)生在藥柱內(nèi)壁縱向1 mm肉厚的環(huán)形區(qū)域內(nèi)。

a. 4500 sb. 56620 s

c. 102210 sd. 113335 s

圖4升溫速率為1.8 K·h-1時發(fā)動機橫向截面的溫度分布云圖

Fig.4Temperature distribution on transverse section of motor at a heating rate of 1.8 K·h-1

圖5為升溫速率3.6 K·h-1條件下，發(fā)動機橫向截面的溫度分布云圖。該升溫速率下，系統(tǒng)溫度分布情況與升溫速率為1.8 K·h-1時相似，只是最初著火位置與著火延遲期有所差別。5400～45800 s，由于加熱速率較慢，系統(tǒng)溫差逐漸減小。69130 s時，藥柱內(nèi)部由于緩慢的自熱反應(yīng)，溫度已經(jīng)高于殼體溫度。73106 s時，AP/HTPB推進劑高溫區(qū)域集中在藥柱內(nèi)壁縱向1 mm肉厚的環(huán)形區(qū)域內(nèi)，該區(qū)域距離藥柱內(nèi)壁面中心100 mm附近，這一區(qū)域即為該升溫速率下的最初著火位置。

a. 5400 sb. 45800 s

c. 69130 sd. 73106 s

圖5升溫速率為3.6 K·h-1時發(fā)動機橫向截面的溫度分布云圖

Fig.5Temperature distribution on transverse section of motor at a heating rate of 3.6 K·h-1

圖6是升溫速率為7.2 K·h-1時，不同時刻發(fā)動機橫向截面的溫度分布云圖。5000 s時，外壁溫度高于藥柱內(nèi)部溫度。隨著傳熱的進行，47000 s時，推進劑內(nèi)部溫度與殼體之間溫差縮小。由于殼體、絕熱層和環(huán)氧樹脂板之間物性的差別以及較大的長徑比，52860 s時，AP/HTPB推進劑在距離藥柱內(nèi)壁中心200 mm附近已經(jīng)發(fā)生了緩慢的自熱反應(yīng)，溫度逐漸高于周圍溫度。53547 s時，最高溫度區(qū)域首先出現(xiàn)在偏向于左端殼體一側(cè)的藥柱內(nèi)壁面上，縱向肉厚約為1 mm，橫向距離藥柱內(nèi)壁中心180 mm附近，這一區(qū)域即為7.2 K·h-1升溫速率下的最初著火位置。

a. 5000 sb. 47000 s

c. 52860 sd. 53547 s

圖6升溫速率為7.2 K·h-1時發(fā)動機橫向截面的溫度分布云圖

Fig.6Temperature distribution on transverse section of motor at a heating rate of 7.2 K·h-1

圖7給出了不同升溫速率下，空氣腔中心(0,0)、藥柱中心(0,4)、藥柱外表面(0,27)、絕熱層外表面(0,28)四個特征點的溫度時程曲線。由圖7可知，起初壁面升溫速率相對較快，由于殼體的導(dǎo)熱系數(shù)與絕熱層和推進劑的導(dǎo)熱系數(shù)相差較大，壁面溫度向內(nèi)部傳導(dǎo)緩慢，因而，殼體與絕熱層和推進劑交界面存在較大的溫度梯度。在緩慢升溫速率下，推進劑內(nèi)部溫度與壁面溫度差距減小，隨著壁面溫度的緩慢升高，AP/HTPB推進劑內(nèi)部開始發(fā)生緩慢的自熱反應(yīng)，導(dǎo)致推進劑內(nèi)部溫度逐漸高于壁面溫度，最終在內(nèi)部發(fā)生著火。當升溫速率分別為1.8，3.6，7.2 K·h-1時，AP/HTPB推進劑的著火延遲期分別為113335 s(31.48 h)，73106 s(20.31 h)，53547 s(14.87 h)，相應(yīng)的著火溫度分別為592，594，595 K。可見，在1.8～7.2 K·h-1范圍內(nèi)，升溫速率對點火溫度影響較小，但著火延遲期隨著升溫速率的增大而大幅縮短。

a. 1.8 K·h-1b. 3.6 K·h-1c. 7.2 K·h-1

圖7不同升溫速率下各特征點溫度隨時間的變化關(guān)系

Fig.7Curves of change in temperature with time for various feature points with different heating rates

4 結(jié) 論

根據(jù)AP/HTPB推進劑在3種升溫速率下慢烤燃過程的數(shù)值模擬結(jié)果，可得出以下結(jié)論：

(1) 根據(jù)文獻[6]和文獻[7,17]中AP/HTPB推進劑的烤燃實驗裝置建立模型，分別對相應(yīng)邊界條件下AP/HTPB推進劑的烤燃過程進行數(shù)值模擬，計算結(jié)果與實驗數(shù)據(jù)吻合較好，證明了所建的慢烤燃模型是合理可行的。

(2) 在1.8，3.6，7.2 K·h-1三種升溫速率下，AP/HTPB推進劑的最初著火位置均發(fā)生在藥柱內(nèi)壁縱向1 mm肉厚的環(huán)形區(qū)域內(nèi)。當加熱速率較小時(1.8 K·h-1)，AP/HTPB推進劑最初的著火位置為藥柱內(nèi)壁中心處附近。隨著升溫速率的增大，最初著火位置由中心區(qū)向藥柱殼體端面移動。

(3) 當升溫速率分別為1.8，3.6，7.2 K·h-1時，AP/HTPB推進劑的著火延遲期分別為31.48，20.31，14.87 h，相應(yīng)的著火溫度分別為592，594，595 K?？梢姡?.8～7.2 K·h-1范圍內(nèi)，升溫速率對著火溫度影響較小，但著火延遲期隨著升溫速率的增大而大幅縮短。

參考文獻:

[1] 馮曉軍，王曉峰，韓助龍. 炸藥裝藥尺寸對慢速烤燃響應(yīng)的研究[J]. 爆炸與沖擊, 2005, 25(3): 285-288.

FENG Xiao-jun, WANG Xiao-feng, HAN Zhu-long. The study of charging size influence on the response of explosives in slow cook-off test[J].ExplosionandShockWaves,2005,25(3): 285-288.

[2] Howard W M, McClelland M A, Nichols A L. ALE3D simulations of gap closure and surface ignition for cookoff Modeling[C]∥13th International Detonation Symposium, Norfolk, VA, United States, 2006.

[3] 王沛，陳朗，馮長根. 不同升溫速率下炸藥烤燃模擬計算分析[J]. 含能材料, 2009, 17(1): 46-49, 54.

WANG Pei, CHEN Lang, FENG Chang-gen. Numerical simulation of cook-off for explosive at different heating rates[J].ChineseJournalofEnergeticMaterials(HannengCailiao),2009,17(1): 46-49,54.

[4] 向梅，黃毅民，饒國寧，等. 不同升溫速率下復(fù)合藥柱烤燃實驗與數(shù)值模擬研究[J]. 爆炸與沖擊, 2013, 33(4): 394-400.

XIANG Mei, HUANG Yi-min, RAO Guo-ning, et al. Cook-off test and numerical simulation for composite charge at different heating rates[J].ExplosionandShockWaves, 2013,33(4): 394-400.

[5] Ho S Y, Ferschl T, Foureur J. Correlation of cook-off behavior of rocket propellants with thermo-mechanical and thermochemical properties[R],ADA274983, 1993.

[6] Komai I, Sato W. Reaction mechanism in slow cook-off Test of GAP-AP propellants[C]∥Insensitive Munitions and Energetic Materials Symposium (IMEMTS),Bristol,UK,24-28 April,2006.

[7] Rodrigo I.Caro, John M.Bellerby. Behavior of hydroxyl-terminated polyether (Htpe) composite rocket propellants in slow cook-off[J].InternationJournalofEnergeticMaterialsandChemicalPropulsion, 2008, 7(3): 171-185.

[8] 陳中娥，唐承志，趙孝彬. HTPB/AP推進劑的慢速烤燃特征[J]. 含能材料, 2006, 14(2): 155-157.

CHEN Zhong-e, TANG Cheng-zhi, ZHAO Xiao-bin. Characteristics of HTPB/AP propellants in slow cook-off[J].ChineseJournalofEnergeticMaterials(HannengCailiao), 2006,14(2): 155-157.

[9] 陳中娥，唐承志，趙孝彬. 固體推進劑的慢速烤燃行為與熱分解特性的關(guān)系研究[J]. 含能材料, 2005, 13(6): 393-396.

CHEN Zhong-e, TANG Cheng-zhi, ZHAO Xiao-bin. Relationship between slow cook-off behaviour and thermal decomposition characteristics of solid propellant[J].ChineseJournalofEnergeticMaterials(HannengCailiao), 2005,13(6): 393-396.

[10] 廖林泉，胥會祥，李勇宏，等. HTPB推進劑危險性實驗研究[J]. 火炸藥學(xué)報, 2010, 33(4): 28-31.

LIAO Lin-quan, XU Hui-xiang, LI Yong-hong, et al. Experimental study on hazard of HTPB propellants[J].ChineseJournalofExplosivesandPropellants, 2010,33(4): 28-31.

[11] 趙孝彬，李軍，程立國，等. 固體推進劑慢速烤燃特性的影響因素研究[J]. 含能材料, 2011, 19(6): 669-672.

ZHAO Xiao-bin, LI Jun, CHENG Li-guo, et al. Influence factors of slow cook-off characteristic for solid propellant[J].ChineseJournalofEnergeticMaterials(HannengCailiao), 2011,19(6): 669-672.

[12] Beckstead M W, Derr R L, Price C F. A Model of Composite Solid-Propellant Combustion Based On Multiple Flames[J].AIAAJournal, 1970, 8(12): 2200-2207.

[13] 陳廣南，張為華. 固體火箭發(fā)動機撞擊與熱安全性分析[M]. 北京: 國防工業(yè)出版社, 2008: 189-191.

CHEN Guang-nan,ZHANG Wei-hua.Safety analyses for solid rocket motors under insults of impact and heat [M]. Beijing: National Defence Industry Press,2008: 189-191.

[14] ZHOU X U, Jackson T L, Buckmaster J. Oscillations in propellant flames with edges[J].CombustionandFlame, 2003, 133(1-2): 157-168.

[15] Hegab A, Jackson T L, Buckmaster J, et al. Nonsteady burning of periodic sandwich propellants with complete coupling between the solid and gas phases[J].CombustionandFlame, 2001, 125(1): 1055-1070.

[16] Wang X, Jackson T L, Massa L. Numerical simulation of heterogeneous propellant combustion by a level set method[J].CombustionTheoryandModelling, 2004, 8(2): 227-254.

[17] Ki-hong Kim, Chang-kee Kim, Ji-chang Yoo. Test-based thermal decomposition simulation of AP/HTPB and AP/HTPE propellants[J].JournalofPropulsionandPower, 2011, 27(4): 822-827.

[18] 江帆，黃鵬. Fluent 高級應(yīng)用與實例分析[M].北京：清華大學(xué)出版社, 2008: 42-45.

JIANG Fan, HUANG Peng. Fluent advanced application and example analysis[M].Beijing: Tsinghua University Press, 2008: 42-45.

[19] Gwak M, Jung T, Yoh J J. Friction-induced ignition modeling of energetic materials[J].JournalofMechanicalScienceandTechnology, 2009, 23(7): 1779-1787.

[20] 原渭蘭，潘浪. 一種艦載導(dǎo)彈固體火箭發(fā)動機烤燃過程的數(shù)值計算方法[J]. 艦船科學(xué)技術(shù), 2009, 31(7): 129-132.

YUAN Wei-lan, PAN Lang. An numerical calculation method on cook-off of solid rocket motor of ship-based missiles[J].ShipScienceandTechnology, 2009, 31(7): 129-132

[21] 武曉松，陳軍，王棟，等. 固體火箭發(fā)動機工作過程數(shù)值仿真[M]. 北京: 高等教育出版社, 2006: 93-97.

WU Xiao-song, CHEN Jun, WANG Dong, et al. Numerical simulations of working progress of solid rocket engine[M].Beijing: Higher Education Press,2006: 93-97.