基于有限元-離散元結(jié)合方法的Steven 實(shí)驗(yàn)三維數(shù)值模擬

黃彬彬，傅 華，喻 寅，劉倉(cāng)理

（中國(guó)工程物理研究院流體物理研究所沖擊波物理與爆轟物理重點(diǎn)實(shí)驗(yàn)室，四川 綿陽 621900）

1 引言

炸藥在跌落、翻倒等低速撞擊事故中易經(jīng)受低幅值（數(shù)十到數(shù)百兆帕）、長(zhǎng)脈沖（數(shù)百微秒）的非沖擊載荷，可能造成炸藥發(fā)生意外點(diǎn)火，引發(fā)安全性事故。但炸藥非沖擊點(diǎn)火的實(shí)驗(yàn)精細(xì)觀測(cè)存在多方面的困難與限制，不易獲得點(diǎn)火演化的細(xì)節(jié)與機(jī)制。因此，開展非沖擊點(diǎn)火的數(shù)值模擬研究，對(duì)深入揭示炸藥在低速撞擊下的響應(yīng)規(guī)律具有重要意義。

Chidester 等［1-3］設(shè) 計(jì) 了Steven 實(shí) 驗(yàn) 用 于 研 究 帶殼體約束的炸藥經(jīng)受低速撞擊下的非沖擊點(diǎn)火問題，建立了撞擊點(diǎn)火模型并采用實(shí)驗(yàn)獲得的臨界撞擊速度對(duì)模型參數(shù)進(jìn)行標(biāo)定。此后，眾多國(guó)內(nèi)外學(xué)者開展了相應(yīng)的數(shù)值模擬工作，分析了不同形狀彈頭加載下的點(diǎn)火速度閾值和點(diǎn)火位置等［4-7］。目前相關(guān)的數(shù)值模擬工作主要采用了有限元方法與離散元方法。其中有限元方法（Finite Element Method，F(xiàn)EM）具有高精度、高效率等優(yōu)點(diǎn)，是目前最成熟且應(yīng)用最廣泛的力學(xué)響應(yīng)計(jì)算方法，但是依賴網(wǎng)格的計(jì)算方法普遍存在不適合模擬大變形和損傷斷裂現(xiàn)象的缺點(diǎn)，雖然基于微裂紋界面的摩擦生熱的點(diǎn)火模型（Visco-SCRAM 模型）在炸藥低速撞擊點(diǎn)火數(shù)值模擬中也有較多應(yīng)用［8-11］，但在處理撞擊加載下炸藥局域點(diǎn)火前的壓剪斷裂與斷裂之后的裂紋摩擦等非連續(xù)性問題時(shí)還是面臨較大的困難。離散元方法（Discrete Element Method，DEM）的基本思想是把不連續(xù)體或連續(xù)體離散為具有一定物理意義的獨(dú)立“微元”或“粒子”，粒子間通過虛彈簧保持連接，粒子間的相互作用力通過彈簧的拉伸和壓縮來實(shí)現(xiàn)，因此在處理斷裂和破碎問題時(shí)具有天然優(yōu)勢(shì)，目前也被廣泛應(yīng)用于炸藥力學(xué)響應(yīng)機(jī)理的模擬研究［12-18］，但是由于離散元不擅長(zhǎng)描述材料的塑性變形，同時(shí)出于計(jì)算尺度和效率的考量，在應(yīng)用于炸藥材料模擬時(shí)，通常模擬中未考慮炸藥外部的金屬殼體對(duì)其非沖擊點(diǎn)火行為的影響。

本研究中建立了有限元-離散元結(jié)合方法，利用有限元模擬裝藥結(jié)構(gòu)響應(yīng)，利用離散元描述炸藥局部斷裂損傷，針對(duì)脆性炸藥PBX-2 的Steven 實(shí)驗(yàn)開展了三維數(shù)值模擬，獲得了裝藥受載應(yīng)力歷史、裂紋分布、溫升曲線及溫度分布。以炸藥裂紋摩擦以及自反應(yīng)放熱為主導(dǎo)點(diǎn)火機(jī)制，分析了由炸藥材料非均質(zhì)特性引起的的概率點(diǎn)火現(xiàn)象，并討論了鈍頭彈與平頭彈撞擊下點(diǎn)火位置的分布特性。

2 計(jì)算方法及模型

2.1 FEM-DEM 結(jié)合方法

FEM-DEM 結(jié)合方法的核心問題在于有限元單元與離散元粒子間接觸作用力的計(jì)算。由于離散元單元本身具有一定尺寸，區(qū)別于有限元節(jié)點(diǎn)僅僅是一個(gè)點(diǎn)，因此離散元粒子-有限元面單元的接觸可以根據(jù)其接觸位置的不同分為三類，即：粒子-表面（Paticles-to-Facets，簡(jiǎn)稱PTF）、粒子-棱邊（Paticles-to-Edges，簡(jiǎn)稱PTE）、粒子-頂點(diǎn)（Paticles-to-Vertices，簡(jiǎn)稱PTV），如圖1所示。

三種接觸方式下的離散元單元所受到的作用力為；

式中，κ為自定義接觸系數(shù)；δn為侵徹距離，m；n為單位法向量；E*為等效楊氏模量，MPa；r*為等效半徑，m；k法向剛度，N·m-1。

圖1 離散元單元與有限元面單元接觸模型Fig.1 Contact model between a discrete-element particle and a finite-element facet unit

2.2 DEM 斷裂模型

PBX-2 各組分的質(zhì)量百分比為HMX/TATB/黏結(jié)劑+鈍感劑=87/7/6，屬于典型的脆性炸藥，因此采用彈脆性本構(gòu)對(duì)其進(jìn)行描述。假設(shè)材料為各向同性，建立的離散元模型如圖2 所示，圖2a 為球形離散元結(jié)構(gòu)示意圖，圖2b 為相鄰單元間作用力示意圖。

圖2 離散元模型示意圖Fig.2 Schematic diagram of discrete-element model

當(dāng)兩個(gè)處于連接狀態(tài)的離散元單元間法向作用力fn與切向作用力fs滿足關(guān)系式（5）中關(guān)系時(shí)，則判定其連接斷裂，兩個(gè)單元由連接狀態(tài)變?yōu)榻佑|狀態(tài)：

式中，fn0和fs0分別為定義的最大法向與切向作用力，N。

為驗(yàn)證該計(jì)算方法以及所選取模型的有效性，對(duì)脆性炸藥PBX-2 的分離式霍普金森桿實(shí)驗(yàn)（SHPB）進(jìn)行了數(shù)值模擬。模型包含炸藥試樣、子彈、整形器、入射桿和透射桿，炸藥試樣尺寸為φ12 mm×6 mm，子彈尺寸為φ20 mm×300 mm，入射桿尺寸為φ20 mm×2000 mm，透射桿尺寸為φ20 mm×1000 mm。炸藥試樣采用離散元?jiǎng)澐志W(wǎng)格，半徑為125 μm，其余部分采用有限元?jiǎng)澐志W(wǎng)格，網(wǎng)格尺寸為2 mm，如圖3 所示。分別于A、B、C三點(diǎn)處記錄入射波、反射波和透射波信息。

“民營(yíng)企業(yè)反映突出的痛點(diǎn)、難點(diǎn)、堵點(diǎn)，給出了6方面35條政策,讓我們真正從政策中增強(qiáng)獲得感。”記者在采訪謝清森時(shí)，他用這句話表達(dá)自己內(nèi)心的最大感受。

圖3 SHPB 實(shí)驗(yàn)計(jì)算模型Fig.3 Calculated model of SHPB test

計(jì)算模型中采用彈脆性本構(gòu)，彈性段彈性模量E=5.5 GPa，平均斷裂應(yīng)變?chǔ)?0.01。模擬結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比如圖4 所示。從圖4a中可以看出，計(jì)算得到的三波（入射波、透射波、反射波）波形與實(shí)驗(yàn)結(jié)果符合較好，表明該計(jì)算方法能較好地處理FEM-DEM 接觸界面間的物理量傳遞。圖4b中應(yīng)力-應(yīng)變本構(gòu)關(guān)系模擬結(jié)果與實(shí)驗(yàn)結(jié)果大致吻合，呈現(xiàn)明顯的彈脆性特征，表明所采用的斷裂模型適用于模擬脆性炸藥材料的力學(xué)響應(yīng)。

圖4 SHPB 實(shí)驗(yàn)?zāi)M結(jié)果Fig.4 Simulation results of SHPB experiment

2.3 溫升點(diǎn)火模型

考慮到相比炸藥塑性功產(chǎn)生的熱量，顆粒間的摩擦溫升更容易產(chǎn)生較高的局部溫升導(dǎo)致點(diǎn)火，因此主要計(jì)算炸藥斷裂后的摩擦生熱以及炸藥自反應(yīng)放熱作為熱源。

兩個(gè)接觸的離散元粒子i與j之間的摩擦力可表示為：

式中，eij為耗散能量，J；Δt為時(shí)間步長(zhǎng)，s。

由于加載過程為瞬態(tài)過程，一個(gè)時(shí)間步長(zhǎng)內(nèi)單元間相互作用可視為絕熱過程，則粒子i的摩擦溫升為：

采用Arrhenius 方程［17］描述炸藥自反應(yīng)放熱

式中，S為源項(xiàng)，J?s-1；ρ為炸藥密度，g?cm-3；Q為反應(yīng)熱，J?kg-1；Z為指前因子，s-1；E為活化能，J?mol-1；R為普適氣體常數(shù)，8.314 J?mol-1?K-1；T為溫度，K；α 為已反應(yīng)炸藥分?jǐn)?shù)。初始溫度設(shè)為298 K，炸藥反應(yīng)動(dòng)力學(xué)參數(shù)見表1［18］。

表1 PBX-2 炸藥反應(yīng)動(dòng)力學(xué)參數(shù)［18］Table 1 Kinetic parameters of explosive reaction for PBX-2［18］

當(dāng)某個(gè)炸藥單元的溫度超過某個(gè)臨界溫度，根據(jù)Arrhenius 方程，將大量放熱導(dǎo)致溫度近乎直線上升超過3000 K，這種情況下我們可以判定炸藥發(fā)生了點(diǎn)火。因此，可以將溫度出現(xiàn)快速上升的拐點(diǎn)作為點(diǎn)火臨界點(diǎn)，其對(duì)應(yīng)的溫度為點(diǎn)火臨界溫度。

2.4 Steven 實(shí)驗(yàn)計(jì)算模型

模型包含炸藥試樣、蓋板、約束環(huán)、樣品盒和彈丸。炸藥試樣尺寸為Φ98 mm×13 mm，壓板尺寸為Φ120 mm×3.5 mm，樣品盒底部高度19 mm，彈丸主體部分半徑為31 mm。炸藥試樣采用離散元?jiǎng)澐志W(wǎng)格（半徑400 μm），其余部分采用有限元?jiǎng)澐志W(wǎng)格（網(wǎng)格尺寸4 mm），建立的FEM-DEM 模型如圖5 所示。

圖5 不同形狀彈頭Steven 實(shí)驗(yàn)計(jì)算模型Fig.5 Calculated model of Steven impact test for different projectile heads

模型中所用樣品為PBX-2 炸藥，彈丸、蓋板和樣品盒的材料為Q235 鋼，約束環(huán)材料為聚四氟乙烯，模型中各材料性能參數(shù)見表2［18］。

實(shí)際加載中，炸藥材料的孔洞坍縮、炸藥-粘結(jié)劑脫粘等行為會(huì)進(jìn)一步導(dǎo)致響應(yīng)局域化，因此模型中需要考慮炸藥材料的非均勻性。為模擬炸藥材料的非均勻性，選取隨機(jī)分布的占比為1%的離散元單元，設(shè)置其最大法向與切向作用力為其余單元的一半。

表2 模型中各材料性能參數(shù)Table 2 Performance parameters for various materials in the model

3 結(jié)果與討論

3.1 力學(xué)響應(yīng)與損傷裂紋

采用有限元/離散元結(jié)合方法得到鈍頭彈在撞擊速度為45 m·s-1時(shí)炸藥樣品底部中心處壓力的模擬結(jié)果，將其與文獻(xiàn)［15］實(shí)測(cè)曲線進(jìn)行對(duì)比，如圖6 所示。由圖6 可以看出，試樣底部中心處的壓力峰值約為0.10 GPa。彈丸撞擊后壓力迅速升高，在約170 μs 處達(dá)到峰值，之后隨彈丸反彈壓力持續(xù)下降，進(jìn)入卸載階段，卸載階段持續(xù)380 μs 左右。數(shù)值模擬的壓力峰值及壓力脈沖寬度都與實(shí)驗(yàn)測(cè)試曲線符合較好，表明該方法能較準(zhǔn)確地模擬Steven 實(shí)驗(yàn)中炸藥的受力過程。

圖6 炸藥底部中心實(shí)測(cè)和計(jì)算壓力曲線（鈍頭彈彈速45 m·s-1）Fig.6 Experimental and simulated pressures at the bottom center of explosives（the velocity of blunt projectile is 45 m·s-1）

圖7 和圖8 為不同彈丸撞擊下炸藥試樣內(nèi)部產(chǎn)生的力學(xué)響應(yīng)及損傷裂紋分布情況。在鈍頭彈實(shí)驗(yàn)中，裂紋由彈頭軸線處沿徑向發(fā)散，形成多段貫穿裂紋（圖7a 和圖7b）；在平頭彈實(shí)驗(yàn)中，裂紋沿彈頭圓周附近呈環(huán)狀分布（圖7c 和圖7d）。對(duì)比圖8 中的應(yīng)力分布情況可知，裝藥試樣受到鈍頭彈撞擊后，撞擊中心的炸藥受到彈丸擠壓（圖7a），該區(qū)域的炸藥向遠(yuǎn)離撞擊中心的方向運(yùn)動(dòng)，應(yīng)力最大值始終集中于軸線附近（圖7b）；裝藥試樣受到平頭彈撞擊后，在彈頭邊緣與蓋板接觸處發(fā)生強(qiáng)烈的剪切作用（圖7c），從而形成繞圓周分布的損傷裂紋，而彈頭下方的炸藥在軸向被壓實(shí)（圖7d）的同時(shí)也會(huì)產(chǎn)生較分散的損傷裂紋。

3.2 點(diǎn)火閾值與點(diǎn)火位置

圖7 不同形狀彈頭撞擊下PBX-2 炸藥的裂紋分布Fig7 Crack distribution of PBX-2 under the impact of different projectile heads

圖8 不同形狀彈頭撞擊下PBX-2 炸藥的應(yīng)力分布Fig.8 Stress distribution of PBX-2 under the impact of different projectile heads

圖9 為鈍頭彈與平頭彈以不同速度撞擊時(shí)炸藥內(nèi)部最高溫度單元的溫度-時(shí)間曲線。由圖9 可以看出，當(dāng)鈍頭彈與平頭彈速度低于40 m·s-1時(shí)，炸藥溫度上升比較緩慢，當(dāng)達(dá)到一定溫度后便不再升高，未發(fā)生點(diǎn)火。當(dāng)鈍頭彈與平頭彈的速度分別大于45 m·s-1和50 m·s-1時(shí)，炸藥溫度以更快的速率上升，達(dá)到點(diǎn)火溫度后溫度急劇上升發(fā)生點(diǎn)火，且點(diǎn)火時(shí)間隨彈速的增加而減小。分析可知，撞擊速度較低時(shí)，炸藥裂紋出現(xiàn)的時(shí)間較晚，摩擦溫升比較緩慢；隨著撞擊速度增加，裂紋間摩擦加劇，自反應(yīng)放熱加快，形成局部熱點(diǎn)導(dǎo)致點(diǎn)火。

圖9 不同速度彈丸撞擊下的溫升曲線Fig.9 Temperature rise curve under the impact of different projectile velocities

圖10 為炸藥達(dá)到點(diǎn)火臨界點(diǎn)時(shí)的溫度分布。鈍頭彈加載下，達(dá)到點(diǎn)火溫度的炸藥單元集中于彈頭軸線下方附近；平頭彈加載下，達(dá)到點(diǎn)火溫度的炸藥單元主要位于彈頭圓周下方，呈環(huán)狀分布。結(jié)合圖7 中的裂紋分布與圖8 中的應(yīng)力分布可知，鈍頭彈加載下，軸線下方炸藥單元所承受的壓力更大，發(fā)生斷裂后單元間的摩擦作用更強(qiáng)，導(dǎo)致更高的摩擦溫升，加速自反應(yīng)放熱，從而更快到達(dá)點(diǎn)火溫度；平頭彈加載下，彈頭邊緣與蓋板接觸處發(fā)生強(qiáng)烈的剪切作用，發(fā)生斷裂后此處的炸藥單元間相對(duì)位移速度更大，導(dǎo)致更高的摩擦溫升，加速自反應(yīng)放熱，從而更快到達(dá)點(diǎn)火溫度。

圖10 不同形狀彈頭撞擊下炸藥的溫度分布圖Fig.10 Temperature distribution of explosives under the impact of different projectile heads

表3 鈍頭彈不同撞擊速度下的點(diǎn)火情況Table 3 Ignition situation of blunt projectile under the impact of different velocities

重復(fù)隨機(jī)生成炸藥非均勻模型，對(duì)每速度撞擊下的點(diǎn)火情況各進(jìn)行25 次模擬，模擬結(jié)果見表3。借鑒落錘實(shí)驗(yàn)炸藥感度判定標(biāo)準(zhǔn)［16］（一定落高下重復(fù)實(shí)驗(yàn)點(diǎn)火概率超過50%），可得數(shù)值模擬中使炸藥點(diǎn)火的速度閾值分別約為43 m·s-1（鈍頭彈）和47 m·s-1（平頭彈），與文獻(xiàn)中結(jié)果（鈍頭彈43 m·s-1，平頭彈47 m·s-1）對(duì)比大致吻合，表明本文中建立的計(jì)算模型可以較好地模擬Steven 實(shí)驗(yàn)。

圖11 為表3 中點(diǎn)火概率100%時(shí)的25 次模擬中點(diǎn)火位置的平面分布統(tǒng)計(jì)：鈍頭彈撞擊后，點(diǎn)火位置集中分布于距離軸心4～11 mm 處，即點(diǎn)火點(diǎn)圍繞撞擊軸產(chǎn)生；對(duì)于平頭彈，點(diǎn)火位置集中分布于距離軸心27～32 mm 處，點(diǎn)火點(diǎn)呈環(huán)狀分布，與彈頭直徑大致相當(dāng)。圖11 中的點(diǎn)火位置分布與文獻(xiàn)［5］中實(shí)驗(yàn)觀測(cè)到的結(jié)果大致相符。

圖11 不同形狀彈頭撞擊下炸藥點(diǎn)火位置平面分布統(tǒng)計(jì)Fig.11 Plane distribution statistics of explosive ignition positions under the impact of different projectile heads

4 結(jié)論

本研究采用FEM-DEM 結(jié)合方法對(duì)Steven 實(shí)驗(yàn)中脆性炸藥PBX-2 的低速撞擊破壞和點(diǎn)火行為進(jìn)行了研究。

（1）建立了基于FEM-DEM 結(jié)合方法的三維非沖擊點(diǎn)火計(jì)算模型并編寫了Fortran 計(jì)算程序，通過對(duì)PBX-2 炸藥的SHPB 實(shí)驗(yàn)?zāi)M，驗(yàn)證了程序處理FEM-DEM 界面物理量傳遞的可靠性，以及所采用的DEM 斷裂模型模擬脆性炸藥斷裂破壞過程的合理性。

（2）考慮了炸藥斷裂破碎、裂紋摩擦升溫與炸藥熱分解，模擬了Steven 實(shí)驗(yàn)中在鈍頭彈與平頭彈加載下PBX-2 炸藥的受載應(yīng)力歷史、裂紋分布、溫升曲線等結(jié)果。鈍頭彈實(shí)驗(yàn)中，裂紋由彈頭軸線處沿徑向發(fā)散，形成多段貫穿裂紋，中心底部形成明顯凹陷；在平頭彈實(shí)驗(yàn)中，裂紋沿彈頭圓周附近呈環(huán)狀分布。

（3）基于考慮了非均質(zhì)特性的足夠數(shù)量PBX 炸藥樣本的模擬，給出了考慮隨機(jī)分散特性的Steven 實(shí)驗(yàn)中的炸藥點(diǎn)火閾值速度與點(diǎn)火位置的統(tǒng)計(jì)結(jié)果。模擬獲得的鈍頭彈與平頭彈撞擊下引起點(diǎn)火的閾值速度分別為43 m·s-1和47 m·s-1。在鈍頭彈加載下點(diǎn)火位置集中分布在撞擊中軸線附近；在平頭彈加載下形成點(diǎn)火環(huán)狀帶，且環(huán)狀直徑與彈頭直徑相當(dāng)。本研究從炸藥裂紋摩擦以及自反應(yīng)放熱的角度解釋脆性炸藥的非沖擊點(diǎn)火現(xiàn)象，模擬結(jié)果與實(shí)驗(yàn)報(bào)道定性相符。