炸藥沖塞點(diǎn)火的數(shù)值模擬

張韓宇，段卓平，孫文遂，黃風(fēng)雷

(北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京100081)

0 引言

炸藥沖塞試驗(yàn)(也稱作大型跌落試驗(yàn))是評價藥柱同時受到撞擊、剪切以及摩擦等力綜合作用下炸藥安全性的重要試驗(yàn)之一。國外從20 世紀(jì)80年代起開展此項(xiàng)工作的研究［1］，如美國洛斯阿拉莫斯實(shí)驗(yàn)室、海軍水面作戰(zhàn)中心等對PBX-9404、TNT、B-3 等炸藥進(jìn)行了沖塞試驗(yàn)。在國內(nèi)，申春迎等［2］對高聚物粘結(jié)炸藥進(jìn)行了沖塞試驗(yàn)，得到了一些炸藥發(fā)生點(diǎn)火的跌落高度以及對應(yīng)的跌落撞擊速度

數(shù)值計算方面，F(xiàn)rey［3］使用線粘塑性材料本構(gòu)方程和阿累尼烏斯方程計算了炸藥在剪切機(jī)械刺激下的點(diǎn)火反應(yīng)，忽略了應(yīng)力波的傳播，計算模型為準(zhǔn)穩(wěn)態(tài)模型，并未反映出材料的力學(xué)不穩(wěn)定性。隨著損傷力學(xué)的發(fā)展，損傷變量引入到炸藥非沖擊點(diǎn)火模型之中，Bennett 等［4］建立了粘彈性統(tǒng)計裂紋力學(xué)模型(viso-SCRAM)，Todd 等［5］建立了損傷點(diǎn)火反應(yīng)模型(DMGIR)，這兩個模型均能較好地預(yù)測炸藥在非沖擊作用下的點(diǎn)火反應(yīng)，但模型中損傷力學(xué)參數(shù)眾多，需要大量試驗(yàn)獲得相應(yīng)參數(shù)，其推廣使用尚需時日。炸藥沖塞屬于機(jī)械刺激范圍，其特點(diǎn)為低強(qiáng)度、長脈沖。本文采用節(jié)點(diǎn)約束-分離方法、聯(lián)合材料模型和化學(xué)動力學(xué)方程描述炸藥沖塞過程和點(diǎn)火反應(yīng)，考察單元尺寸對數(shù)值模擬的影響。

1 節(jié)點(diǎn)約束－分離方法

炸藥沖塞點(diǎn)火問題是熱-力-化耦合問題。炸藥沖塞過程中單元承受較大的變形，易出現(xiàn)嚴(yán)重的網(wǎng)格畸變，炸藥通常采用歐拉算法，它雖然可以解決網(wǎng)格的大變形問題，但不能夠精確體現(xiàn)炸藥破壞的物質(zhì)邊界，同時也影響炸藥破壞面溫度的描述。而傳統(tǒng)的拉格朗日算法計算大變形和破壞會產(chǎn)生網(wǎng)格畸變，一般采用侵蝕算法，即將炸藥設(shè)為連續(xù)介質(zhì)，當(dāng)炸藥受到拉壓后發(fā)生嚴(yán)重變形，達(dá)到指定數(shù)值后刪除單元網(wǎng)格，這樣可避免計算中出現(xiàn)網(wǎng)格畸變和負(fù)體積等問題，但這與實(shí)際炸藥裂紋生成、擴(kuò)展或者損傷斷裂等有很大差別，且在刪除單元后，炸藥的熱力學(xué)(溫度)計算會發(fā)生錯誤而終止。因此，上述方法均不能有效解決炸藥沖塞點(diǎn)火過程中的熱-力-化耦合問題。

本文數(shù)值模擬中炸藥材料采用節(jié)點(diǎn)約束-分離方法描述其力學(xué)破壞。節(jié)點(diǎn)約束-分離方法是指將節(jié)點(diǎn)按參與構(gòu)成單元的個數(shù)進(jìn)行復(fù)制，相鄰單元坐標(biāo)相同的節(jié)點(diǎn)相互獨(dú)立，通過定義單元應(yīng)變失效值對所有空間位置重疊的節(jié)點(diǎn)族進(jìn)行“捆綁”［6］。當(dāng)單元的應(yīng)變達(dá)到指定失效值時，節(jié)點(diǎn)之間的約束失效，相鄰單元坐標(biāo)相同的節(jié)點(diǎn)相互分離。本文采用自編程序進(jìn)行節(jié)點(diǎn)拆分，圖1給出了本文計算實(shí)例中某4 個單元的節(jié)點(diǎn)約束-分離示意圖。當(dāng)達(dá)到單元應(yīng)變失效值時，節(jié)點(diǎn)相互分離，單元演變成獨(dú)立單元。在計算中設(shè)置自接觸，單元之間不會發(fā)生穿透，單元碰撞時依然有相互作用。采用此方法時不需要刪除炸藥單元，材料仍采用拉格朗日算法，可以描述炸藥的大變形、斷裂破壞等現(xiàn)象以及宏觀裂紋的生成、擴(kuò)展，并支持炸藥的熱-力-化耦合計算。歐拉算法和傳統(tǒng)的拉格朗日算法均不能有效解決炸藥沖塞點(diǎn)火過程中的熱-力-化耦合問題。

圖1 4 個單元的節(jié)點(diǎn)約束-分離示意圖Fig.1 Sketch map of node tie-breaking of four elements

2 數(shù)值計算模型

2.1 計算模型

文獻(xiàn)［2］中試驗(yàn)裝置如圖2所示，炸藥粘結(jié)在代用材料內(nèi)，沖塞從鋼板中間的孔穿過，整個試驗(yàn)件從高空跌落撞擊靶板，沖塞鉆入炸藥，對炸藥主要產(chǎn)生剪切和撞擊作用。炸藥尺寸為φ152 mm×102 mm，沖塞與炸藥的接觸直徑為28 mm，沖塞長38 mm.試驗(yàn)測得PBX-2 炸藥臨界點(diǎn)火的跌落撞擊速度為27.7 m/s，等效自由跌落高度為39.1 m.

采用LS-DYNA 有限元軟件建立二維模型，計算模型采用拉格朗日算法。由圖2可以看出，試驗(yàn)件為對稱結(jié)構(gòu)，1/2 計算模型如圖3所示，采用軸對稱實(shí)體單元劃分網(wǎng)格，計算時直接將試件碰撞鋼靶的跌落撞擊速度加載于試驗(yàn)件上。根據(jù)最小單元尺寸由程序來設(shè)置結(jié)構(gòu)計算的初始時間步長。為了熱學(xué)計算的穩(wěn)定性，人為設(shè)置熱計算時間步長與結(jié)構(gòu)計算時間步長相同，并設(shè)置步長縮放因子為0.67，確保數(shù)值計算穩(wěn)定進(jìn)行［7-8］。

圖2 炸藥沖塞試驗(yàn)裝置示意圖Fig.2 Configuration of spigot test assembly

圖3 炸藥沖塞的1/2 計算模型Fig.3 1/2 calculation domain for spigot test

2.2 材料模型

炸藥熱-力-化耦合計算中，炸藥材料模型包括力學(xué)模型和熱學(xué)模型兩部分。

炸藥力學(xué)模型采用熱彈塑性本構(gòu)模型［8］描述。

式中:T 為由炸藥熱彈塑性材料和熱分解放熱反應(yīng)共同引起的炸藥溫度為總應(yīng)變率;為熱應(yīng)變率，α 為線膨脹系數(shù)，δij為克羅內(nèi)克符號;是與溫度相關(guān)的彈性矩陣。

式中E、ν 分別為為楊氏模量和泊松比，與溫度相關(guān)。

材料屈服函數(shù)為

式中:sij為應(yīng)力偏量;σy(T)=Ep(T)εpeff+σ0(T)為屈服應(yīng)力;εpeff為等效塑性應(yīng)變;σ0為初始屈服應(yīng)力;Ep為塑性硬化模量，σ0和Ep均為溫度的函數(shù)。

炸藥熱學(xué)模型采用各向同性熱傳導(dǎo)模型描述，熱傳導(dǎo)微分方程為

式中:ρ 為材料密度;cv為材料比熱;k 為導(dǎo)熱率;Qs為材料內(nèi)部的熱源，在此Qs=η σij+Qt，η 為功轉(zhuǎn)熱比率，Qt為放熱量，如(5)式所示。

對于炸藥材料而言，溫度上升會引起材料的熱分解放熱反應(yīng)，以阿累尼烏斯方程來表示炸藥的熱分解放熱量。炸藥熱分解放熱量為

式中:Q 為反應(yīng)熱;A 為指前因子;E 為活化能;R 為氣體常數(shù)。

代用材料、沖塞和靶板的熱學(xué)模型均采用各向同性熱傳導(dǎo)模型描述，代用材料、沖塞和靶板的熱學(xué)模型采用塑性隨動硬化模型描述。PBX-2 炸藥的熱力學(xué)參數(shù)如表1所示，炸藥失效應(yīng)變?yōu)?.688%［9］.

3 計算結(jié)果

3.1 數(shù)值計算方法和材料模型適應(yīng)性驗(yàn)證

圖4給出了不同時刻炸藥沖塞破壞和溫度分布，圖中的溫度為炸藥節(jié)點(diǎn)的溫度，炸藥單元的溫度為該單元各個節(jié)點(diǎn)的平均溫度。由圖4可以看出，溫度較高的區(qū)域基本位于沖塞破壞面，這是由于塑性功轉(zhuǎn)化為熱，熱量在非常短的時間內(nèi)聚集在材料局部區(qū)域。

表1 PBX-2 炸藥熱力學(xué)參數(shù)［10］Tab.1 Thermodynamic parameters of PBX-2［10］

圖4 不同時刻炸藥沖塞破壞及溫度(v=30 m/s)Fig.4 Evolutions of intrusion failure and temperature distributions of explosive (v=30 m/s)

圖5為炸藥的破壞情況，由圖可以看出，使用節(jié)點(diǎn)約束-分離方法可以有效地避免網(wǎng)格發(fā)生畸變，而且能較好地模擬出炸藥裂紋形成、破壞等力學(xué)行為。

圖5 炸藥的節(jié)點(diǎn)分離和網(wǎng)格變形以及局部放大Fig.5 Node breaking and mesh deformation of explosive and close up of local surface

沖塞過程中，隨著炸藥溫度的升高，熱分解速率將遵循阿累尼烏斯規(guī)律迅速增加［3，11］。圖6給出了跌落撞擊速度為30 m/s 時計與不計化學(xué)動力學(xué)方程計算的炸藥最高溫度曲線對比。由圖可知，計化學(xué)動力學(xué)方程計算炸藥的溫升時，在574 μs 時刻，炸藥單元的溫度曲線出現(xiàn)拐點(diǎn)，即dT/dt→∞，炸藥開始發(fā)生點(diǎn)火反應(yīng)，點(diǎn)火溫度為770 K.不計化學(xué)動力學(xué)方程計算炸藥的溫升時，炸藥溫度雖然有一定的升高，但并沒有出現(xiàn)溫度拐點(diǎn)，炸藥未發(fā)生點(diǎn)火反應(yīng)。因此，模擬炸藥點(diǎn)火反應(yīng)需要計及炸藥熱分解放熱反應(yīng)對溫升和點(diǎn)火的貢獻(xiàn)。從上述數(shù)值計算可以看出，節(jié)點(diǎn)約束-分離方法、聯(lián)合熱彈塑性材料模型和化學(xué)動力學(xué)方程適用于描述炸藥的沖塞破壞和點(diǎn)火反應(yīng)，清晰地顯示出炸藥破壞面的形成，得到了炸藥的點(diǎn)火溫升曲線。

3.2 單元尺寸對數(shù)值模擬影響及臨界速度確定

炸藥沖塞破壞屬于應(yīng)變局部化問題，數(shù)值計算結(jié)果依賴于有限元網(wǎng)格密度［12］。下面考察0.4 mm、0.5 mm 和0.8 mm 3 種不同單元尺寸對沖塞破壞和溫度分布等計算結(jié)果的影響。

圖7為相同跌落撞擊速度(v =36 m/s)條件下采用3 種不同單元尺寸計算的炸藥沖塞破壞和溫度分布(260 μs 時刻)。單元尺寸0.4 mm 時沖塞深度為7.4 mm;單元尺寸0.5 mm 時沖塞深度為7.5 mm;單元尺寸0.8 mm 時沖塞深度為7.6 mm.可以看出，破壞面形成和整體破壞效果比較一致，炸藥的力學(xué)破壞受單元尺寸影響較小。但是破壞面上單元網(wǎng)格變形不同，引起的塑性變形溫升也不同，單元尺寸對炸藥沖塞破壞達(dá)到的最高溫度有一定影響。

圖6 計及/不計化學(xué)動力學(xué)方程計算的炸藥最高溫度曲線對比(v=30 m/s)Fig.6 Comparison of peak temperature of explosive from the models with/without chemical kinetics equation (v=30 m/s)

圖8給出了相同跌落撞擊速度(v =36 m/s)條件下3 種單元尺寸計算的炸藥所能達(dá)到的最高溫度對比。單元尺寸0.4 mm 時點(diǎn)火時間308 μs，點(diǎn)火溫度為765 K;單元尺寸0.5 mm 時點(diǎn)火時間336 μs，點(diǎn)火溫度為758 K;單元尺寸0.8 mm 時點(diǎn)火時間為434 μs，點(diǎn)火溫度為731 K.可以看出，隨著網(wǎng)格加密，單元尺寸越小，溫升越快，炸藥點(diǎn)火時間越短。

從計算耗費(fèi)的時間來看，單元尺寸0.4 mm 耗費(fèi)機(jī)時約為131 min，單元尺寸0.5 mm 耗費(fèi)機(jī)時約為65 min，單元尺寸0.8 mm 耗費(fèi)機(jī)時約為22 min.單元尺寸越小，花費(fèi)的時間將成倍增加。

圖9給出了單元尺寸為0.4 mm 計算32 m/s、30 m/s、29 m/s 3 種跌落撞擊速度下炸藥最高溫度的曲線對比?？梢钥闯?，跌落撞擊速度越大，炸藥的溫升越高。當(dāng)?shù)渥矒羲俣葹?2 m/s 時，點(diǎn)火時間為406 μs，點(diǎn)火溫度為756 K;跌落撞擊速度為30 m/s時，點(diǎn)火時間574 μs，點(diǎn)火溫度為770 K;跌落撞擊速度為29 m/s 時，炸藥未發(fā)生點(diǎn)火反應(yīng)，因此，單元尺寸0.4 mm 下炸藥沖塞點(diǎn)火的最小跌落撞擊速度約為29.5 m/s.另外，單元尺寸0.5 mm 的最小跌落撞擊速度為31.5 m/s，單元尺寸0.8 mm 的最小跌落撞擊速度為35.5 m/s.單元尺寸0.4 mm 和0.5 mm 計算的點(diǎn)火反應(yīng)和最小跌落撞擊速度趨于一致。對比3 種單元尺寸下炸藥沖塞破壞、點(diǎn)火反應(yīng)、計算機(jī)時以及最小跌落撞擊速度來看，單元尺寸0.4 mm 計算結(jié)果已基本消除了網(wǎng)格尺寸帶來的計算誤差，能夠反應(yīng)PBX-2 炸藥沖塞過程中的力學(xué)破壞和點(diǎn)火反應(yīng)。因此得到炸藥沖塞點(diǎn)火的最小跌落撞擊速度為29.5 m/s，與試驗(yàn)測得的臨界跌落撞擊速度為27.7 m/s［3］吻合較好。

圖7 3 種單元尺寸計算的炸藥沖塞破壞和溫度分布對比(v=36 m/s，t=260 μs)Fig.7 Comparison of intrusion failures and temperature distributions of explosive with three different sizes of elements(v=36 m/s，t=260 μs)

計算結(jié)果與試驗(yàn)結(jié)果相比略高，這是由于數(shù)值模擬采用的材料模型僅考慮炸藥沖塞過程中塑性功轉(zhuǎn)熱引起的炸藥溫升，而沒有考慮材料裂紋之間摩擦生熱等對炸藥溫升的貢獻(xiàn)。

4 結(jié)論

圖8 3 種單元尺寸計算的炸藥最高溫度曲線對比(v=36 m/s)Fig.8 Comparison of peak temperature of explosive with three different sizes of elements (v=36 m/s)

圖9 3 種跌落撞擊速度條件下的炸藥最高溫度曲線對比Fig.9 Comparison of peak temperature of explosive calculated with three different drop velocities

1)采用節(jié)點(diǎn)約束-分離方法、聯(lián)合熱彈塑性材料模型和化學(xué)動力學(xué)方程對炸藥沖塞點(diǎn)火反應(yīng)進(jìn)行了數(shù)值模擬，該方法和模型能夠有效地描述炸藥沖塞破壞過程和點(diǎn)火反應(yīng)。

2)數(shù)值計算中，單元尺寸對炸藥沖塞過程中的力學(xué)破壞影響較小，但對炸藥溫升、點(diǎn)火反應(yīng)以及臨界跌落撞擊速度有較大影響。對比0.4 mm、0.5 mm和0.8 mm 3 種單元尺寸的計算結(jié)果，單元尺寸0.4 mm 的計算結(jié)果能夠反應(yīng)PBX-2 炸藥沖塞過程中的力學(xué)破壞和點(diǎn)火反應(yīng)，并與試驗(yàn)結(jié)果比較吻合。

References)

［1］ Lee P R.Hazard assessment of explosives and propellants［M］∥Zukas J A，Walters W P.Explosive effects and applications.New York:Springer Press，1997.

［2］ 申春迎，向永，代曉淦，等.高聚物黏結(jié)炸藥的沖塞試驗(yàn)研究［J］.火炸藥學(xué)報，2010，33(2):29 -32.SHEN Chun-ying，XIANG Yong，DAI Xiao-gan，et al.Study on the spigot test of polymer bonded explosives［J］.Chinese Journal of Explosives and Propellants，2010，33(2):29 -32.(in Chinese)

［3］ Frey R B.The initiation of explosive charges by rapid shear［C］∥Proceedings of the 7th Symposium (International)on Detonation.Annapolis，USA:Office of Naval Research，1981:36 -42.

［4］ Bennett J G，Haberman K S，Johnson J N，et al.A constitutive model for the non-shock ignition and mechanical response of high explosives［J］.Journal of the Mechanics and Physics of Solids，1998，46(12):2303 -2322.

［5］ Todd S N.Non-shock initiation model for plastic bonded explosive PBX-5:theoretical results［J］.Shock Compression of Condensed Matter，2007，955(1):1006 -1009.

［6］ 張曉天，賈光輝，黃海.基于節(jié)點(diǎn)分離Lagrange 有限元方法的超高速碰撞碎片云數(shù)值模擬［J］.爆炸與沖擊，2010，30(5):499 -504.ZHANG Xiao-tian，JIA Guang-hui，HUANG Hai.Simulation of hypervelocity-impact debris clouds using a Lagrange FEM with node separation［J］.Explosive and Shock Waves，2010，30(5):499 -504.(in Chinese)

［7］ 惲壽容，涂侯杰.爆炸力學(xué)計算方法［M］.北京:北京理工大學(xué)出版社，1995.YUN Shou-rong，TU Hou-jie.Explosion mechanics calculation method［M］.Beijing:Beijing Institute of Technology Press，1995.(in Chinese)

［8］ Hallquist J O.Ls-dyna theory manual［M］.Livermore，USA:Lawrence Livermore National Laboratory，2006.

［9］ 羅景潤.PBX 的損傷、斷裂及本構(gòu)關(guān)系研究［D］.綿陽:中國工程物理研究院，2001.LUO Jing-run.Study on damage，fracture and constitutive relation of PBX［D］.Mianyang:China Academy of Engineering Physics，2001.(in Chinese)

［10］ 董海山，周芬芬.高能炸藥及相關(guān)物性能［M］.北京:科學(xué)出版社，1989.DONG Hai-shan，ZHOU Fen-fen.Properties of high explosives and related materials［M］.Beijing:Science Press，1989.(in Chinese)

［11］ 孫寶平，段卓平，黃風(fēng)雷.炸藥摩擦點(diǎn)火評價實(shí)驗(yàn)數(shù)值模擬［J］.北京理工大學(xué)學(xué)報，2011，31(6):638 -642.SUN Bao-ping，DUAN Zhuo-ping，HUANG Feng-lei.Numerical simulation for assessment test of friction ignition in solid explosive［J］.Transaction of Beijing Institute of Technology，2011，31(6):638 -642.(in Chinese)

［12］ Needleman A.Material rate dependence and mesh sensitivity in localization problems［J］.Computer Methods in Applied Mechanics and Engineering，1988，67:69 -85.