張坤, 智小琦, 肖游, 王帥

(1.中北大學(xué) 機(jī)電工程學(xué)院, 山西 太原 030051;2.中國兵器裝備集團(tuán)自動(dòng)化研究所有限公司 智能制造事業(yè)部, 四川 綿陽 621000;3.湖北航天化學(xué)技術(shù)研究所 航天化學(xué)動(dòng)力技術(shù)重點(diǎn)實(shí)驗(yàn)室, 湖北 襄陽 441003)

0 引言

研究彈藥在熱環(huán)境中的響應(yīng)特性對(duì)彈藥在壽命期的安全使用具有重要意義。因?yàn)閺椝師o論是受熱刺激、強(qiáng)電磁輻射或二者的綜合作用,最終導(dǎo)致彈藥響應(yīng)的均是因含能材料受熱后的自持反應(yīng)所致。目前,模擬彈藥受熱刺激作用的研究方法主要是快速烤燃和慢速烤燃?？焖倏救寄M彈藥在大火作用下的響應(yīng)特性,慢速烤燃模擬彈藥在暗火作用下的響應(yīng)特性。研究手段主要采用試驗(yàn)和仿真技術(shù)。試驗(yàn)方法主要采用北大西洋公約組織標(biāo)準(zhǔn)化協(xié)定STANAG 4240(快速烤燃試驗(yàn))和STANAG 4382(慢速烤燃試驗(yàn)),試驗(yàn)研究內(nèi)容包括彈藥尺寸[1-2]、約束條件[3-4]、裝藥密度及間隙[5-6]、含能材料配方[7]、熔鑄炸藥點(diǎn)火機(jī)理[8]以及緩釋結(jié)構(gòu)特性[9]等對(duì)彈藥響應(yīng)溫度、響應(yīng)烈度的影響研究,已達(dá)到既深又廣的程度;仿真研究與試驗(yàn)研究既相輔相成又互相補(bǔ)存和促進(jìn),能獲得點(diǎn)火時(shí)刻的溫度及其梯度、點(diǎn)火點(diǎn)位置、點(diǎn)火前壓力[10]、內(nèi)部機(jī)理的分析[11-12]及點(diǎn)火后的瞬時(shí)響應(yīng)特性[13]等內(nèi)容。其中關(guān)于點(diǎn)火后的含能材料自持反應(yīng)特性的仿真研究軟件,目前只有美國ALE3D軟件和Uintah軟件。

從上述情況可見,研究彈藥的熱刺激,無論從試驗(yàn)還是仿真都取得了很大的成就。但是關(guān)于純理論研究點(diǎn)火點(diǎn)位置和點(diǎn)火溫度的文章,目前鮮有報(bào)道,而理論研究是彈藥熱刺激研究不可或缺的內(nèi)容,是對(duì)熱爆炸理論的完善和補(bǔ)存,具有重要的意義。本文擬采用疊加原理和分離變量方法,首次將炸藥非反應(yīng)性熱傳導(dǎo)與自熱反應(yīng)熱傳導(dǎo)拆分,以一維炸藥模型為基礎(chǔ),研究慢速烤燃的溫度分布及點(diǎn)火點(diǎn)位置的解析解,分析點(diǎn)火點(diǎn)位置的一維變化規(guī)律,以期為烤燃研究奠定理論基礎(chǔ)。

1 一維數(shù)學(xué)理論模型

1.1 理論推演的假設(shè)

為簡化理論計(jì)算,需對(duì)炸藥的烤燃過程作如下假設(shè):

1)炸藥為凝聚炸藥,炸藥在加熱過程中不發(fā)生相變;

2)彈體外壁均勻升溫,升溫速率為固定值;

3)忽略殼體與炸藥間的接觸熱阻;

4)炸藥的熱導(dǎo)率、比熱容、密度均為恒值;

5)炸藥活化能、指前因子、反應(yīng)熱均為常數(shù);

6)炸藥按照零階動(dòng)力學(xué)定律分解。

1.2 慢速烤燃一維模型

炸藥的熱傳導(dǎo)方程用式(1)[14]描述:

(1)

(2)

當(dāng)則炸藥厚度L=2δ,其中δ為厚度值的一半。則一維無限平板模型的熱傳導(dǎo)可用式(3)來描述:

(3)

式中:α為熱擴(kuò)散系數(shù),且α=λ/ρc;S(x,t)為熱傳導(dǎo)方程的源項(xiàng),用式(4)表示:

(4)

式(3)的定解條件如下:

1)初始條件:T(x,0)=Tr,其中Tr為初始溫度;

2)邊界條件:T(±δ,t)=κt+Tr,其中κ為升溫速率。

2 求解

2.1 溫度場的拆分

當(dāng)沒有化學(xué)反應(yīng)時(shí),炸藥熱傳導(dǎo)問題的解析解即非反應(yīng)性熱傳導(dǎo)比較容易求解。但是,具有化學(xué)反應(yīng)的熱傳導(dǎo)方程即反應(yīng)性熱傳導(dǎo)是極難獲得解析解的。為進(jìn)一步求解熱傳導(dǎo)方程式(3),需要用到疊加原理[15],將原問題T(x,t)分解為非反應(yīng)性熱傳導(dǎo)項(xiàng)Tt(x,t)和自熱反應(yīng)熱傳導(dǎo)項(xiàng)Ts(x,t)兩個(gè)子問題。其中非反應(yīng)熱傳導(dǎo)式為

(5)

式(5)相應(yīng)的定解條件如下:

1)初始條件:Tt(x,0)=Tr;

2)邊界條件:Tt(±δ,t)=κt+Tr。

自熱反應(yīng)熱傳導(dǎo)的方程如下:

(6)

式(6)的定解條件如下:

1)初始條件:Ts(x,0)=0;

2)邊界條件:Ts(±δ,t)=0。

則原方程(式(3))的解T(x,t)可通過T(x,t)=Tt(x,t)+Ts(x,t)求得。

2.2 非反應(yīng)性熱傳導(dǎo)問題的求解

式(5)的熱傳導(dǎo)子問題不包含自熱源項(xiàng),即屬于沒有化學(xué)反應(yīng)的惰性情況下的熱傳導(dǎo)問題。由于變量x、t相互獨(dú)立,使用分離變量法[16]將Tt(x,t)寫為

Tt(x,t)=X(x)U(t)

(7)

式中:X(x)為位置x的函數(shù);U(t)為時(shí)間t的函數(shù)。

將式(7)代入式(5),有如下等式:

(8)

(9)

式中:k為常數(shù)。

式(8)根據(jù)邊界條件Tt(±δ,t)=κt+Tr,通過Strum-Liouville理論求解有

Xn(x)=Cnsin[nπ(-x+δ)/(2δ)]

(10)

式中:Cn為常數(shù),X0(x)=C0;固有值λn=nπ,n=0,1,2,3,…。

對(duì)式(9)求解,有

Un(t)=Bne(-απ2n2t)/4δ2

(11)

式中:Bn為常數(shù)。

則非反應(yīng)性熱傳導(dǎo)方程式(5)的通解為

(12)

式中:An為常數(shù),且對(duì)于Tt(x,t)解函數(shù)首項(xiàng)Tt0(x,t)有Tt0(x,t)=A0。

代入初始條件Tt(x,0)=Tr,進(jìn)一步求得非反應(yīng)性熱傳導(dǎo)方程的精確解Tt(x,t)為

(13)

隨著時(shí)間t的不斷增加,無窮級(jí)數(shù)項(xiàng)將逐漸減小,非反應(yīng)性熱傳導(dǎo)所產(chǎn)生的溫度分布將逐漸穩(wěn)定。其中初始溫度Tr作為常數(shù)項(xiàng),不影響非反應(yīng)性熱傳導(dǎo)所產(chǎn)生的溫度分布。

2.3 自熱反應(yīng)熱傳導(dǎo)問題的求解

對(duì)于式(6)的自熱反應(yīng)熱傳導(dǎo)問題,非齊次項(xiàng)(即源項(xiàng))為Arrhenius方程,直接求解是相當(dāng)復(fù)雜的。通過進(jìn)一步將源項(xiàng)S(x,t)泰勒展開:

(14)

當(dāng)自熱反應(yīng)產(chǎn)生的溫度Ts小于非反應(yīng)性熱傳導(dǎo)溫度Tt的1%時(shí),源項(xiàng)S(x,t)可取泰勒展開式首項(xiàng)近似代替,作近似后的自熱反應(yīng)熱傳導(dǎo)方程如下:

(15)

式中:

(16)

對(duì)于式(15)的求解,采用分離變量法[17]。將Ts(x,t)、S(x,t)分別按固有函數(shù)展開為Dn(t)、Sn(t),于是有

(17)

式中:

(18)

(19)

再將其反變換后,則有

(20)

這樣,可解得自熱反應(yīng)熱傳導(dǎo)方程的解Ts(x,t)為

Ts(x,t)=

(21)

觀察解的結(jié)構(gòu),不難發(fā)現(xiàn)式(21)中Ts(x,t)包含Arrhenius溫度積分,直接求出其精確的結(jié)果極其困難[18]。但Ts(x,t)中常數(shù)顯然不影響對(duì)其最高溫度所在位置的求解,即指數(shù)前因子A、反應(yīng)熱Q不影響自熱反應(yīng)熱傳導(dǎo)溫度最高值所在位置,進(jìn)而不影響炸藥在烤燃中的點(diǎn)火點(diǎn)位置。隨著升溫速率κ的增加,炸藥兩側(cè)升溫加快,則點(diǎn)火點(diǎn)位置向邊界移動(dòng)。

3 慢速烤燃試驗(yàn)

3.1 試驗(yàn)設(shè)計(jì)

設(shè)計(jì)尺寸為φ75 mm×295 mm烤燃彈,殼體材料選用45號(hào)鋼,兩端帶螺紋的端蓋連接,端蓋和殼體壁厚均為5 mm。內(nèi)部裝藥為甘肅銀光化學(xué)工業(yè)集團(tuán)有限公司生產(chǎn)的RDX炸藥,單節(jié)藥柱直徑為57.4 mm,長度為47.5 mm,共6節(jié),藥柱總長285 mm,長徑比為5,裝藥密度為1 610 kg/m-3。

為獲得烤燃試驗(yàn)過程中,炸藥內(nèi)部的溫度的變化情況,在每節(jié)藥柱中間刻槽埋入φ1 mm的上海自動(dòng)化儀表股份有限公司生產(chǎn)的WRNK191型鎧裝電偶(精度為0.1 K,熱響應(yīng)時(shí)間≦3 s)。藥柱與藥柱之間以及藥柱與熱電偶之間均采用蟲膠漆粘接固定。熱電偶的溫度測點(diǎn)由上至下編號(hào)1～5號(hào),兩相鄰測點(diǎn)及邊緣測點(diǎn)至端面的間距均為47.5 mm,溫度測點(diǎn)位置如圖1所示。由于采用等長壓制藥柱試驗(yàn),無法獲得其他位置的溫度情況。

圖1 溫度測點(diǎn)的位置示意圖Fig.1 Location of temperature measuring points

試驗(yàn)彈試驗(yàn)前狀態(tài)如圖2所示。試驗(yàn)采用圓柱形保溫箱進(jìn)行加熱,試驗(yàn)裝置如圖3所示。

圖2 試驗(yàn)彈體Fig.2 Tested ammunition

圖3 慢速烤燃試驗(yàn)裝置Fig.3 Slow-cook test setup

使用日本島電公司生產(chǎn)的MR13溫控儀對(duì)烤燃彈外壁進(jìn)行升溫控制,并使用美國Fluke公司的1586A多路測溫儀采集各測點(diǎn)的溫度數(shù)據(jù),采樣周期為1 s。從室溫305 K開始,以3.3 ℃/h的固定速率升溫,直至烤燃彈發(fā)生響應(yīng)。

3.2 試驗(yàn)結(jié)果

試驗(yàn)藥柱中5支熱電偶和外壁1支熱電偶在點(diǎn)火時(shí)刻的測量結(jié)果如表1所示,試驗(yàn)中各測點(diǎn)溫度變化曲線如圖4所示。由于炸藥內(nèi)部的自熱反應(yīng)隨溫度的升高而逐漸加劇,臨近響應(yīng)時(shí),各測點(diǎn)的溫度迅速升高。從5條測溫曲線可見,藥柱的兩側(cè)1號(hào)和5號(hào)熱電偶測得的溫度最高,3號(hào)點(diǎn)測得的溫度最低,可見點(diǎn)火區(qū)域位于藥柱兩側(cè),點(diǎn)火時(shí)外壁溫度為468.9 K。

表1 點(diǎn)火時(shí)刻熱電偶的測量結(jié)果

圖4 溫度變化曲線Fig.4 Temperature history curves

烤燃彈響應(yīng)后,試驗(yàn)點(diǎn)產(chǎn)生較深的彈坑,加熱套筒與烤燃彈殼體完全碎裂。由于試驗(yàn)在野外溝壑中進(jìn)行,響應(yīng)后僅回收到少部分碎片,如圖5所示。根據(jù)響聲及以往同種炸藥的烤燃試驗(yàn)情況,綜合判斷烤燃彈的響應(yīng)等級(jí)為爆轟。

圖5 試驗(yàn)結(jié)果Fig.5 Test results

4 算例及驗(yàn)證

4.1 檢驗(yàn)理論與試驗(yàn)的一致性

為檢驗(yàn)理論與試驗(yàn)的一致性,以試驗(yàn)藥柱為例,利用表2中的物性參數(shù)和動(dòng)力學(xué)參數(shù),對(duì)式(21)的解析解進(jìn)行計(jì)算。

表2 RDX炸藥物性參數(shù)和動(dòng)力學(xué)參數(shù)

對(duì)于任意長圓柱,單從徑向或軸向考慮,傳熱均屬于一維問題,僅考慮徑向傳熱時(shí),因溫度沿縱軸均勻分布,沒有區(qū)別。而單考慮軸向傳熱時(shí),因溫度分布不均勻,能體現(xiàn)出位置的差別。同理,僅考慮自熱兩個(gè)方向也均屬于一維問題,故對(duì)藥柱軸向的溫度分布進(jìn)行計(jì)算[19]。

已知炸藥厚度 、初始溫度為305 K、升溫速率κ=3.3 ℃/h。使用Maple數(shù)學(xué)軟件[20]和Clenshaw-Curtis正交方法[21]求解,控制無窮級(jí)數(shù)的舍入誤差不大于0.1%。理論解的藥柱沿軸線方向的部分溫度值如表3所示。

表3 沿軸線方向的溫度值

外壁溫度為484.6 K時(shí),理論解與測點(diǎn)測量值的溫度分布如圖6所示。各測點(diǎn)測量得到的溫度分布趨勢與理論解的溫度分布趨勢相符,但是由于試驗(yàn)值是三維結(jié)果,而理論僅是一維結(jié)果,因此兩曲線的彎曲程度有所不同(即對(duì)于試驗(yàn)藥柱,沿軸向的傳熱屬于一維問題[19])。

圖6 溫度分布Fig.6 Temperature distribution

為直觀地比較理論解的相對(duì)溫度分布與試驗(yàn)測得的相對(duì)溫度分布,將圖6中試驗(yàn)測得的溫度值、理論解中截取與試驗(yàn)測量范圍(-0.095～0.095 m)對(duì)應(yīng)區(qū)間的溫度值,再根據(jù)各自的溫度幅值分別作如下歸一化處理:

(22)

式中:xN為歸一化后的溫度值;xi為歸一化前的溫度值;xmax為歸一化前數(shù)據(jù)中的最高溫度值;xmin為歸一化前數(shù)據(jù)中的最低溫度值。通過歸一化將試驗(yàn)測得的溫度值、理論解與試驗(yàn)對(duì)應(yīng)位置的溫度值映射到0～1之間,得到各個(gè)位置的相對(duì)溫度分布,如圖7所示。理論解各個(gè)位置的溫度之間的比例與試驗(yàn)結(jié)果一致,即理論解的相對(duì)溫度分布與試驗(yàn)測得的相對(duì)溫度分布相符,進(jìn)一步驗(yàn)證了理論解的溫度分布的正確性。

圖7 相對(duì)溫度分布Fig.7 Relative temperature distribution

4.2 數(shù)值驗(yàn)證

為進(jìn)一步驗(yàn)證理論解的正確性,使用數(shù)值方法對(duì)解析解式(21)進(jìn)行驗(yàn)證。Suceka[22]對(duì)熱傳導(dǎo)方程式(1)進(jìn)行了有限差分近似。對(duì)于無限大平板而言,關(guān)于位置x的偏導(dǎo)數(shù)可以通過以下差分來近似:

(23)

式中:Δx為將厚度為L的炸藥被均勻分為n份后的寬度;上下標(biāo)用于區(qū)分網(wǎng)格對(duì)應(yīng)的時(shí)間和位置,上標(biāo)j為時(shí)間序列,下標(biāo)i為位置序列。

對(duì)時(shí)間t的偏導(dǎo)數(shù)進(jìn)行同樣的有限差分近似:

(24)

式中:Δt為時(shí)間增量。

使用有限差分替代熱傳導(dǎo)方程式(1)中的偏導(dǎo)數(shù),即得到炸藥熱傳導(dǎo)離散方程:

(25)

同理,建立非反應(yīng)性熱傳導(dǎo)離散方程如下:

(26)

對(duì)于式(25)與式(26)具有相同定解的條件如下:

通過式(25)可以獲得任意時(shí)刻炸藥內(nèi)部的溫度分布;將式(25)與式(26)同時(shí)刻的溫度結(jié)果作差,即可獲得該時(shí)刻炸藥自熱反應(yīng)熱傳導(dǎo)的溫度分布。

以上兩種顯式方法的穩(wěn)定性由空間步長和時(shí)間步長的比值決定。對(duì)于熱傳導(dǎo)方程,如果滿足以下穩(wěn)定條件[23],則該方法是穩(wěn)定的:

(27)

穩(wěn)定性條件保證了有限差分解的穩(wěn)定性,該條件僅在炸藥分解緩慢的情況下,能保證解的精確性。對(duì)于慢速烤燃,點(diǎn)火時(shí)自熱反應(yīng)產(chǎn)生的熱量很低,故適用于該判據(jù)。

計(jì)算選取RDX炸藥為研究對(duì)象,計(jì)算中使用的物性參數(shù)和動(dòng)力學(xué)參數(shù)的值如表2所示,炸藥厚度、初始溫度Tr為305 K、升溫速率κ=3.3 ℃/h。代入以上參數(shù),對(duì)熱傳導(dǎo)離散方程式(25)、熱傳導(dǎo)離散方程式(26)賦予初始條件與邊界條件,以時(shí)間增量Δt=1 s循環(huán)迭代離散方程(即j+1時(shí)刻的溫度分布由j時(shí)刻的溫度分布計(jì)算得出),并將各個(gè)時(shí)刻的溫度分布進(jìn)行輸出。

點(diǎn)火點(diǎn)的溫度變化如圖8所示。在D點(diǎn)(加熱時(shí)間約為177 000 s時(shí))發(fā)生點(diǎn)火,自熱反應(yīng)逐漸加劇,溫度不再可控,進(jìn)而發(fā)生響應(yīng)。

圖8 點(diǎn)火點(diǎn)溫度的理論解與數(shù)值解Fig.8 Theoretical and numerical results of the temperature of ignition point

炸藥點(diǎn)火溫度的理論計(jì)算結(jié)果與采用文獻(xiàn)[22]方法的數(shù)值解具有良好的匹配度,理論計(jì)算值為468.8 K;數(shù)值計(jì)算值為468.6 K,誤差不超過1%。

點(diǎn)火時(shí)刻(外壁溫度為480.7 K時(shí)),炸藥總體溫度(即非反應(yīng)性熱傳導(dǎo)溫度與自熱反應(yīng)熱傳導(dǎo)溫度之和)分布與僅自熱反應(yīng)所產(chǎn)生的溫度分布如 圖9 所示,點(diǎn)火點(diǎn)位置位于炸藥兩側(cè),與試驗(yàn)結(jié)果相符。總體溫度分布主要受升溫邊界的影響,由于點(diǎn)火前達(dá)到分解能壘的分子數(shù)很少,自熱反應(yīng)所產(chǎn)生的溫度較低。自熱反應(yīng)所產(chǎn)生的溫度的理論計(jì)算結(jié)果與數(shù)值解相比,最大誤差為8.8%(誤差主要來自于式(14)對(duì)源項(xiàng)的近似處理);自熱反應(yīng)溫度最高值所在位置的計(jì)算結(jié)果與數(shù)值解相比,最大誤差不超過1%。

圖9 炸藥厚度L=0.285 m時(shí)的理論計(jì)算結(jié)果與數(shù)值計(jì)算結(jié)果Fig.9 Theoretical and numerical results with explosive thickness L=0.285 m

5 結(jié)果與討論

自熱反應(yīng)溫度最高值的位置是由炸藥自熱反應(yīng)與邊界的熱散失共同作用的結(jié)果。圖10為不同厚度條件下,炸藥自熱反應(yīng)溫度最高值隨時(shí)間的變化情況。由圖10可見,隨著時(shí)間的增加,自熱反應(yīng)溫度最高值所在的相對(duì)位置(即|x/δ|)將趨于穩(wěn)定,直至發(fā)生點(diǎn)火。

圖10 不同厚度炸藥自熱反應(yīng)溫度最高值所在的相對(duì)位置Fig.10 Relative location of the maximum value point of self-heating reaction temperature of explosiveswith different thicknesses

在炸藥厚度較薄的情況下,自熱反應(yīng)溫度最高位置將迅速移動(dòng)至中心;隨著炸藥厚度L的增加,位置移動(dòng)的相對(duì)量逐漸減小,當(dāng)炸藥厚度L>0.3 m時(shí),從50 000 s直至發(fā)生點(diǎn)火,自熱反應(yīng)溫度最高位置的相對(duì)移動(dòng)量不大于2%。顯然可以認(rèn)為在烤燃試驗(yàn)進(jìn)行到一定時(shí)間之后,溫度相對(duì)分布趨于穩(wěn)定,自熱反應(yīng)溫度最高值位置隨時(shí)間不再發(fā)生改變,直至點(diǎn)火。這樣,可以通過理論計(jì)算,預(yù)則最高溫度及其所在位置,為降低響應(yīng)等級(jí)的結(jié)構(gòu)設(shè)計(jì)提供便利。

點(diǎn)火時(shí),自熱反應(yīng)溫度最高值所在位置至邊界的距離變化如圖11所示。當(dāng)炸藥厚度L>0.3 m時(shí),隨著炸藥厚度的進(jìn)一步增加,自熱反應(yīng)溫度最高值所在位置至邊界的距離變化不大于2%,且最高值所在位置至邊界的距離逐漸趨于穩(wěn)定。即隨著炸藥厚度的增加,邊界的熱散失條件逐漸成為影響點(diǎn)火點(diǎn)位置的重要因素。

圖11 自熱反應(yīng)溫度最高值所在位置至邊界的距離Fig.11 Distance from the maximum value point of self-heating reaction temperature to boundary

當(dāng)炸藥厚度L>0.3 m時(shí),從邊界至中心的自熱反應(yīng)溫度梯度如圖12所示,隨著炸藥厚度的增加,點(diǎn)火點(diǎn)附近的自熱反應(yīng)溫度梯度趨于恒定。隨著炸藥厚度的增加,自熱反應(yīng)溫度分布的規(guī)律相似。

圖12 炸藥邊界到中心的自熱反應(yīng)溫度梯度Fig.12 Temperature gradient of self-heating reaction from boundary to center of explosives

6 結(jié)論

本文通過理論解析的方法研究一維凝聚炸藥慢速烤燃的點(diǎn)火點(diǎn)位置及點(diǎn)火溫度。得出以下主要結(jié)論:

1)將炸藥非反應(yīng)性熱傳導(dǎo)與反應(yīng)性熱傳導(dǎo)拆分,通過理論方程及計(jì)算可得出慢烤條件下一維凝聚炸藥內(nèi)部的溫度分布情況及點(diǎn)火點(diǎn)位置,并與試驗(yàn)及數(shù)值計(jì)算結(jié)果進(jìn)行對(duì)比,驗(yàn)證了理論方程及計(jì)算的正確性。

2)通過理論計(jì)算可得到一維炸藥點(diǎn)火時(shí)刻僅自熱反應(yīng)所產(chǎn)生的溫度分布和自熱反應(yīng)溫度最高值所在位置;獲得點(diǎn)火時(shí)刻一維凝聚炸藥溫度最高值所在位置及溫度梯度沿厚度的變化規(guī)律。為烤燃試驗(yàn)點(diǎn)火溫度及位置的預(yù)測提供了可靠的理論依據(jù)。