RDX基PBX藥柱烤燃過程的數(shù)值模擬

蒲翰濤, 王 星, 趙寒月, 陳科全, 蔣道建, 路中華

(1. 中國工程物理研究院化工材料研究所, 四川 綿陽 621999; 2. 中國工程物理研究院安全彈藥研發(fā)中心, 四川 綿陽 621999; 3. 中物院高性能數(shù)值模擬軟件中心, 北京 100080; 4. 北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所, 北京 100080)

1 引 言

在儲(chǔ)存、運(yùn)輸和使用過程中,彈藥有遭受意外火災(zāi),發(fā)生燃燒、爆炸甚至殉爆的風(fēng)險(xiǎn), 因此合理評估炸藥熱安全性非常重要。目前對彈藥的熱安全性研究主要是烤燃試驗(yàn)和數(shù)值模擬,烤燃試驗(yàn)是針對裝藥的宏觀響應(yīng)特性進(jìn)行直接觀測,數(shù)值模擬則針對炸藥熱分解過程進(jìn)行綜合分析,具有成本低、周期短、數(shù)據(jù)全面等優(yōu)點(diǎn),與烤燃試驗(yàn)相輔相成,是目前研究的熱點(diǎn)。

20世紀(jì)90年代中期,烤燃數(shù)值模擬主要是編譯程序,對炸藥點(diǎn)火特征信息進(jìn)行計(jì)算。Jones[1]基于Fortran語言開發(fā)了HEAT軟件,對烤燃過程中的炸藥溫度分布進(jìn)行了研究。而90年代末至今,相關(guān)專業(yè)軟件的發(fā)展為烤燃過程的研究提供了更好的平臺,Howard[2],McClelland[3]等利用ALE3D軟件對多種炸藥烤燃試驗(yàn)過程中的反應(yīng)動(dòng)力學(xué)進(jìn)行了研究。國內(nèi)專家學(xué)者也對炸藥烤燃試驗(yàn)做了大量的研究工作。由于國外專業(yè)軟件的禁運(yùn),我國的研究工作主要是借助現(xiàn)行商業(yè)軟件。王沛[4],陳朗[5]張曉立[6]陳科全[7]等人分別利用FLUENT、LS-DYNA3D、ABAQUS軟件對烤燃試驗(yàn)全過程進(jìn)行了分析。限于商業(yè)軟件的功能,這些研究并未針對炸藥的反應(yīng)動(dòng)力學(xué)進(jìn)行細(xì)致表征。

為此,本課題組自主編寫了“含能材料動(dòng)態(tài)響應(yīng)數(shù)值模擬軟件”(Dynamic Response Simulation package for Energetic Materials),并使用該軟件中的反應(yīng)動(dòng)力學(xué)模塊,針對黑索今(RDX)基高聚物粘結(jié)炸藥(PBX)的烤燃過程,進(jìn)行了數(shù)值模擬,研究其在不同升溫率下相關(guān)組分的變化規(guī)律。

2 計(jì)算方法及模型校驗(yàn)

2.1 計(jì)算方法

在加熱過程中,熱量以熱傳導(dǎo)形式通過殼體向炸藥內(nèi)部傳遞,受熱后,炸藥會(huì)發(fā)生復(fù)雜的多步熱分解反應(yīng),先將熱能轉(zhuǎn)化為化學(xué)能,再將化學(xué)能轉(zhuǎn)化為熱能釋放出來,隨著熱量的積累,反應(yīng)速率會(huì)逐漸加快,最終導(dǎo)致炸藥點(diǎn)火。

對于以RDX為基的PBX炸藥,其受熱分解過程可用(1)式、(2)和(3)式[8],相關(guān)物理化學(xué)參數(shù)見表1。

A→Br1=Z1exp-E1/RTA

(1)

B→Cr2=Z2exp-E2/RTB

(2)

C→Dr3=Z3exp-E3/RTC2

(3)

式中，r為反應(yīng)速率,mol·L-1·s-1;Z為指前因子,s-1;E為活化能,J·mol-1;m為反應(yīng)序號,m=1,2,3;R為普適氣體常數(shù),8.314 J·mol-1·K-1;A為RDX炸藥;B為固體中間產(chǎn)物;C為氣體中間產(chǎn)物;D為氣體終產(chǎn)物; [A],[B],[C],[D]分別為其對應(yīng)的質(zhì)量分?jǐn)?shù),kg·m-3。反應(yīng)(1)是固體吸熱一級化學(xué)反應(yīng),反應(yīng)(2)是固體中間產(chǎn)物向氣體中間產(chǎn)物轉(zhuǎn)變的放熱一級化學(xué)反應(yīng),反應(yīng)(3)是氣體中間產(chǎn)物分解為氣體終產(chǎn)物的強(qiáng)放熱二級化學(xué)反應(yīng),也是炸藥點(diǎn)火的主導(dǎo)反應(yīng)。

表1 炸藥和殼體的物性參數(shù)及反應(yīng)動(dòng)力學(xué)參數(shù)[8-9]

Table 1 Physical parameters and chemical kinetic parameters for RDX-based PBX and steel

materialdensity/kg·m-3specificheatcapacity/J·kg-1·K-1thermalconductivity/W·m-1·K-1mreactionheat/J·kg-1activationenergy/J·mol-1pre-exponentialconstant/1014s-1RDX164013300.2512.68×1061.97×1056.39843×10172-8.03×1061.85×1054.74008×10173-6.560×1071.43×1059.54×1014steel8030502.4816.27----

因此其動(dòng)力學(xué)方程可表示為:

(4)

(5)

(6)

(7)

在二維柱坐標(biāo)系中,熱量傳遞方程的形式為[10]:

(8)

S=Q1Z1exp-E1/RTA+Q2Z2exp-E2/RTB+

Q3Z3exp-E3/RTC2

(9)

Tr,z,t=Ts

(10)

式中,ρ是密度,kg·m-3;c是比熱,J·kg-1·K-1;T是溫度,K;λ是熱導(dǎo)率,W·m-1·K-1;S為炸藥自熱反應(yīng)源項(xiàng),Q為反應(yīng)的反應(yīng)熱。方程(10)表示第一類邊界條件,其意義是在邊界上給定溫度值。

2.2 計(jì)算流程

計(jì)算流程如圖1所示,主要包括多步化學(xué)反應(yīng)動(dòng)力學(xué)計(jì)算和和熱傳導(dǎo)兩部分。多步化學(xué)反應(yīng)動(dòng)力學(xué)部分包括計(jì)算化學(xué)組份變化及熱源項(xiàng),其難點(diǎn)是: (1)化學(xué)反應(yīng)動(dòng)力學(xué)計(jì)算與宏觀的烤燃過程存在較強(qiáng)的時(shí)間多尺度問題,并且化學(xué)反應(yīng)控制方程也存在較強(qiáng)的剛性。為了削弱這些問題,程序根據(jù)溫度梯度來自適應(yīng)調(diào)節(jié)子程序中的時(shí)間步長,當(dāng)溫度梯度超過一定限制時(shí)將時(shí)間步長縮減為原來的1/10。(2)化學(xué)反應(yīng)初期反應(yīng)物消耗較低導(dǎo)致數(shù)值計(jì)算存在誤差問題,程序以產(chǎn)物的生成量作為反應(yīng)進(jìn)度的標(biāo)準(zhǔn)以減小數(shù)值誤差。

熱傳導(dǎo)部分主要是利用有限元的基本思想進(jìn)行計(jì)算,包括建立左端矩陣、源項(xiàng),以及邊界條件加載。對于時(shí)間積分部分,主要是考慮到在烤燃前期溫度變化緩慢,故采用隱式方法進(jìn)行時(shí)間離散,取較大的時(shí)間步長進(jìn)行計(jì)算; 在烤燃后期,由于溫度急劇升高,計(jì)算中所取時(shí)間步長必須足夠小,故采用顯式算法進(jìn)行時(shí)間離散。顯式算法和隱式算法之間的耦合切換通過溫度變化率調(diào)整設(shè)定。

圖1 計(jì)算流程圖

Fig.1 Sketch diagram of calculation method

2.3 模型校驗(yàn)

采用文獻(xiàn)[9]中的試驗(yàn)結(jié)果校驗(yàn)本文計(jì)算方法。試驗(yàn)裝置由加熱烘箱、殼體、藥柱、熱電偶等組成[9],如圖2所示。初始溫度為300 K,在升溫速率為1 K·min-1,殼體厚度為3 mm,藥柱尺寸為Φ40 mm×50 mm的條件下測量藥柱中心(M點(diǎn))和殼體外壁(N點(diǎn))溫度,直至裝藥結(jié)構(gòu)點(diǎn)火,并記錄點(diǎn)火時(shí)間于表2中。

圖2 加熱烘箱和烤燃試件[9]

1—炸藥, 2—?dú)んw, 3，5—熱電偶, 4—加熱烘箱

Fig.2 Schematic diagram of heater and cook-off test specimen

1—explosive, 2—shell, 3,5—thermocouple, 4—heater

從表2中可以看出,本文計(jì)算的點(diǎn)火時(shí)間與試驗(yàn)值存在-6.8%的偏差,計(jì)算的藥柱中心和殼體外壁溫度與試驗(yàn)值分別存在-3.55%和-0.86%的偏差,分析原因可能有以下三點(diǎn):

(1)試驗(yàn)中,下端蓋與凸起臺階緊密接觸,流入熱量略小于殼體外壁。本研究為了簡化計(jì)算,將試件外側(cè)以及上下表面全部施加第一類邊界條件,導(dǎo)致其總熱量的積累要高于試驗(yàn)環(huán)境,因此點(diǎn)火時(shí)間縮短,中心點(diǎn)和邊界點(diǎn)的溫度也低于試驗(yàn)值。

表2 本文計(jì)算結(jié)果與試驗(yàn)結(jié)果[9]

Table 2 Comparison of the numerical simulation results and experimental ones[9]

itemignitiontime/sdeviation/%temperatureatthepointM/Kdeviation/%temperatureatthepointN/Kdeviation/%experiment[9]10560-458.35-467.95-thiswork9842-6.8442.1-3.55463.91-0.86

(2)試驗(yàn)件的殼體與炸藥在局部存在少量間隙,熱量在此處傳遞以熱對流和熱輻射為主,這與數(shù)值模擬中的無縫接觸有一定區(qū)別,也將導(dǎo)致偏差產(chǎn)生。

(3)數(shù)值模擬中沒有考慮RDX在高溫下會(huì)發(fā)生一定程度的相變,這也會(huì)導(dǎo)致偏差產(chǎn)生。

但在實(shí)際烤燃過程中,殼體的導(dǎo)熱系數(shù)遠(yuǎn)大于炸藥,少量間隙對總體溫度影響較小,并且僅有極少量炸藥在溫度急劇升高的區(qū)域可能發(fā)生相變,表2對比結(jié)果表明本研究建立的計(jì)算模型較為合理,故以本程序?yàn)榛A(chǔ)進(jìn)行了RDX基PBX炸藥烤燃試驗(yàn)數(shù)值模擬研究。

3 計(jì)算結(jié)果及分析

3.1 RDX基PBX炸藥烤燃試驗(yàn)數(shù)值模擬

針對裝藥尺寸為Φ40 mm×50 mm,殼體厚度為3 mm的烤燃模型,在升溫速率為1 K·min-1條件下進(jìn)行計(jì)算,起爆點(diǎn)、中心點(diǎn)、邊界點(diǎn)溫度歷程曲線如圖3所示。在加熱初期,熱量由殼體向炸藥內(nèi)部傳遞,殼體溫度上升較快,炸藥內(nèi)部溫度變化較小,在殼體和炸藥內(nèi)部存在明顯的溫度梯度,其原因是殼體的熱導(dǎo)率遠(yuǎn)大于炸藥的熱導(dǎo)率。隨著時(shí)間推移,熱量逐漸積累,起爆點(diǎn)溫度在9790 s時(shí)刻反超邊界點(diǎn)溫度,之后迅速升高,直至9842 s時(shí)溫度梯度無窮大,計(jì)算終止,此時(shí)炸藥處于熱失控狀態(tài),即認(rèn)為發(fā)生了點(diǎn)火。

烤燃過程中,點(diǎn)火區(qū)域組分含量曲線如圖4所示。在加熱初期,各物質(zhì)的含量變化極小,400 K時(shí),初級反應(yīng)開始加快,至457 K時(shí),有1%RDX發(fā)生初級分解反應(yīng); 與此同時(shí),固體中間產(chǎn)物B的含量逐漸增加,在457 K時(shí),其含量為0.9%; 次級反應(yīng)由于固體中間產(chǎn)物B的含量較低,故反應(yīng)時(shí)間拖后,在450 K左右開始才開始有明顯變化; 460 K時(shí),第三步反應(yīng)開始加快; 當(dāng)氣體終產(chǎn)物D的質(zhì)量分?jǐn)?shù)為0.006時(shí),在第三步反應(yīng)強(qiáng)烈熱效應(yīng)的作用下,引發(fā)點(diǎn)火。因此氣體終產(chǎn)物的含量可以作為炸藥能否點(diǎn)火的判據(jù)之一。

圖3 三個(gè)特征點(diǎn)溫度隨時(shí)間的變化曲線

Fig.3 Temperature as a function of time at three characteristic points

圖4 1 K·min-1時(shí)炸藥點(diǎn)火區(qū)域各組分的質(zhì)量分?jǐn)?shù)

Fig.4 Mass fraction of each component in the ignition point of explosive at a heating rate of 1 K·min-1

3.2 升溫速率對炸藥烤燃響應(yīng)特性的影響

為了研究在烤燃過程中特征點(diǎn)溫度和組分含量對升溫速率的響應(yīng)規(guī)律,針對藥柱尺寸Φ40 mm×50 mm,殼體厚度為3 mm的模型,計(jì)算了0.5,1,5,10 K·min-1共4種升溫速率下RDX基PBX炸藥烤燃過程,為了簡化計(jì)算,忽略升溫速率對炸藥熱導(dǎo)率、比熱容等參數(shù)的影響,結(jié)果如圖5,圖6所示。

圖5 點(diǎn)火時(shí)間隨升溫速率的變化曲線

Fig.5 Ignition time as a function of heating rate

圖6 點(diǎn)火時(shí)刻的特征點(diǎn)溫度隨升溫速率的變化曲線

Fig.6 Temperature of characteristic points at ignition time as a function of heating rate

圖5為點(diǎn)火時(shí)間隨升溫速率變化的曲線。由圖5可見,隨著升溫速率提高,點(diǎn)火時(shí)間急劇衰減。圖6為點(diǎn)火時(shí)刻,模型各特征點(diǎn)溫度隨升溫速率變化的曲線。由圖6可見,隨著升溫速率提高,中心點(diǎn)的溫度近似線性下降,其原因是點(diǎn)火時(shí)間縮短,中心點(diǎn)接受的熱量隨之減少; 起爆點(diǎn)溫度基本不變,其原因是該溫度主要由熱分解反應(yīng)峰值決定,與升溫速率無關(guān); 邊界點(diǎn)溫度增大,其原因是升溫速率的提高在邊界點(diǎn)的作用最明顯。由此表明: 升溫速率越低,點(diǎn)火時(shí)環(huán)境溫度越低,炸藥危險(xiǎn)性越高。

圖7和圖8反映了在點(diǎn)火時(shí)刻,起爆區(qū)域各分組份含量與升溫速率的關(guān)系。由圖7可見,隨著升溫速率減小,RDX炸藥A和固體中間產(chǎn)物B含量變化較小。其原因是初級和次級反應(yīng)對升溫速率不敏感。由圖8可見,隨著升溫速率減小,氣體中間產(chǎn)物C和氣體終產(chǎn)物D明顯增大。其原因是第三步強(qiáng)放熱化學(xué)反應(yīng)的反應(yīng)程度增大。由此表明: 升溫速率越低,氣體終產(chǎn)物D的含量就越高,炸藥反應(yīng)就越劇烈,炸藥的危險(xiǎn)性就越高。

圖7 點(diǎn)火時(shí)刻起爆區(qū)域A,B組分含量隨升溫速率的變化曲線

Fig.7 Mass fraction of components A and B in the initiation area at ignition time as a function of heating rate

圖8 點(diǎn)火時(shí)刻起爆區(qū)域C,D組分含量隨升溫速率的變化曲線

Fig.8 Mass fraction of components C and D in the initiation area at ignition time as a function of heating rate

4 結(jié) 論

(1) RDX基PBX炸藥在烤燃過程中,400 K時(shí)初級分解反應(yīng)開始加速,450 K時(shí)次級分解反應(yīng)明顯發(fā)生,進(jìn)而誘導(dǎo)強(qiáng)放熱的第三步化學(xué)反應(yīng),表現(xiàn)為氣體終產(chǎn)物的釋放,最終發(fā)生點(diǎn)火。

(2) 點(diǎn)火時(shí)間隨升溫速率的增加而急劇衰減,中心點(diǎn)溫度降低,起爆點(diǎn)溫度基本不變,邊界點(diǎn)溫度增大。因此升溫速率越低,炸藥點(diǎn)火的環(huán)境溫度越低,危險(xiǎn)性越高。

(3) 升溫速率越低,在點(diǎn)火時(shí)刻,起爆區(qū)域的氣體終產(chǎn)物D的含量就越高,第三步強(qiáng)放熱反應(yīng)就越劇烈,后續(xù)的爆轟過程中,釋放的能量就更大,炸藥的危險(xiǎn)性就越高。

本研究主要針對熱傳導(dǎo)機(jī)制作用下的烤燃過程數(shù)值模擬,關(guān)于復(fù)雜問題,如熱輻射和熱對流邊界,殼體力學(xué)響應(yīng)等,還需進(jìn)一步探究。

參考文獻(xiàn):

[1] Jones D A, Parker R P. Heat flow calculations for the small-scale cook-off bomb test, AD-A236829[R]. US: DTIC, 1991.

[2] Howard W M, Nichols A L, Mcclelland M A. ALE3D Simulations of Gap Closure and Surface Ignition for Cookoff Modeling[M]. United States. Department of Energy, 2006.

[3] Mcclelland M A, Maienschein J L, Howard W M, et al. ALE3D simulation of heating and violence in a fast cookoff experiment with LX-10[R]. Lawrence Livermore National Laboratory (LLNL), Livermore, CA, 2006.

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

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。

[5] 陳朗, 王沛, 馮長根. 考慮相變的炸藥烤燃數(shù)值模擬計(jì)算[J]. 含能材料, 2009, 17(5): 568-573.

CHEN Lang, WANG Pei, FENG Chang-gen. Numerical study of explosive cook-off with phase change[J].ChineseJournalofEnergeticMaterials(HannengCailiao),2009, 17(5): 568-573.

[6] 張曉立, 洪滔, 秦承森, 等. 烤燃彈熱點(diǎn)火的LS-DYNA數(shù)值模擬研究[J]. 彈箭與制導(dǎo)學(xué)報(bào), 2011, 31(5): 110-112.

ZHANG Xiao-li, HONG Tao, QIN Cheng-cen, et al. LS-DYNA numerical simulation for thermal ignition of cook-off bomb[J].JournalofProjectiles,Rockets,MissilesandGuidance, 2011, 31(5): 110-112.

[7] 陳科全, 黃亨建, 路中華, 等. 基于ABAQUS的PBX炸藥烤燃試驗(yàn)數(shù)值計(jì)算[J]. 火炸藥學(xué)報(bào), 2014, 37(2): 31-36.

CHEN Ke-quan, HUANG Heng-jian, LU Zhong-hua, et al. Numerical calculation of cook-off test for PBX explosive based on ABAQUS[J].ChineseJournalofExplosivesandPropellants, 2014, 37(2): 31-36.

[8] McClelland M A, Tran T D, Cunningham B J, et al. Cookoff response of PBXN-109: material characterization and ALE3D model[C]∥Proceedings of JANNAF 19th Propulsion Systems Hazards Subcommittee Meeting. 2000: 191.

[9] 牛余雷, 南海, 馮曉軍, 等. RDX基PBX炸藥烤燃試驗(yàn)與數(shù)值計(jì)算[J]. 火炸藥學(xué)報(bào), 2011, 34(1): 32-37.

NIU Yu-lei, NAN Hai, FENG Xiao-jun, et al. Cook-off test and its numerical calculation of RDX-based PBX explosive[J].ChineseJournalofExplosivesandPropellants, 2011, 34(1): 32-37.

[10] 馮長根. 熱爆炸理論[M]. 北京: 科學(xué)出版社,1988: 97-117.

FENG Chang-gen. Theory of Thermal Explosion[M]. Beijing: Science Press, 1988: 97-117.