拉伸加載下PBX 炸藥力學(xué)性能的分子動(dòng)力學(xué)模擬

高飛艷，劉 睿，陳鵬萬(wàn)，龍 瑤，陳 軍,3

（1. 北京理工大學(xué)爆炸科學(xué)與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 100081；2. 北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所, 北京 100088；3. 北京大學(xué)應(yīng)用物理與技術(shù)中心, 北京 100871）

聚合物黏結(jié)炸藥（polymer bonded explosive，PBX）是一種高度填充的顆粒復(fù)合材料，主要由含能單質(zhì)炸藥晶粒和聚合物黏結(jié)劑以及其他增塑劑、鈍感劑等混合而成[1-3]。PBX 具有低感度、高爆速、良好的物理安定性等特點(diǎn)，廣泛應(yīng)用于火箭推進(jìn)劑、彈頭和采礦等軍事和民用領(lǐng)域。相對(duì)于單質(zhì)炸藥晶粒，黏結(jié)劑的彈性模量相對(duì)較低且易于變形，在低速撞擊、拉伸和壓縮等外部刺激下，復(fù)雜的熱-力-化學(xué)耦合反應(yīng)可能沿著炸藥晶體與黏結(jié)劑之間的界面產(chǎn)生，從而影響PBX 的穩(wěn)定性[4-5]。因此，研究炸藥晶體與黏結(jié)劑之間的界面行為對(duì)于PBX 的設(shè)計(jì)和應(yīng)用具有重要意義。

目前，人們對(duì)PBX 的力學(xué)行為、結(jié)構(gòu)變形和斷裂機(jī)理開(kāi)展了大量研究。例如：Yeager 等[6]通過(guò)中子反射儀觀察到在HMX 晶體與黏結(jié)劑Estane 之間存在一個(gè)明顯的界面混合層（厚度約6 nm），該混合層可能影響?zhàn)そY(jié)劑與炸藥之間的界面附著力；Rae 等[7]通過(guò)巴西實(shí)驗(yàn)研究了PBX 的拉伸斷裂行為，發(fā)現(xiàn)微裂紋通常在大顆粒的邊界形成，并垂直于加載方向擴(kuò)展，導(dǎo)致PBX 斷裂，在裂紋處黏結(jié)劑被拉成纖維細(xì)絲，這被認(rèn)為會(huì)影響炸藥的拉伸強(qiáng)度；Xiao 等[8]通過(guò)動(dòng)態(tài)拉伸實(shí)驗(yàn)研究了RDX 基PBX 的拉伸性能，發(fā)現(xiàn)其拉伸行為表現(xiàn)出明顯的高應(yīng)變率依賴性，炸藥顆粒與黏結(jié)劑之間的界面脫粘是主要的失效模式。然而，僅通過(guò)實(shí)驗(yàn)很難完全揭示PBX 在拉伸加載下的力學(xué)響應(yīng)，如應(yīng)力局部化、微觀結(jié)構(gòu)演變等。

數(shù)值模擬方法是研究PBX 力學(xué)行為的有效手段，近年來(lái)，人們針對(duì)含能晶體、聚合物黏結(jié)劑及損傷界面提出了不同的本構(gòu)模型。Wang 等[9]考慮了由于黏結(jié)劑破裂和脫粘引起的速率相關(guān)損傷，為黏結(jié)劑與PBX 界面建立了損傷黏彈性模型，并模擬了PBX9501 的動(dòng)態(tài)力學(xué)響應(yīng)，發(fā)現(xiàn)PBX9501 的拉伸強(qiáng)度與應(yīng)變速率有關(guān)，且斷裂方式為黏結(jié)劑脫粘。Dai 等[10]使用內(nèi)聚力模型（cohesive zone modeling，CZM）研究了PBX9501 在單軸準(zhǔn)靜態(tài)拉伸下的失效演化，發(fā)現(xiàn)在整個(gè)拉伸過(guò)程中，拉伸應(yīng)變起源于較軟的黏結(jié)劑，拉伸應(yīng)變的持續(xù)累積導(dǎo)致黏結(jié)劑脫粘。分子動(dòng)力學(xué)（molecular dynamics，MD）方法作為一種廣泛使用的數(shù)值模擬技術(shù)，可以揭示界面局部結(jié)構(gòu)演化和微觀力學(xué)行為[11-14]。Xiao 等[15-17]應(yīng)用COMPASS 力場(chǎng)計(jì)算了一系列HMX 基PBX 炸藥的彈性力學(xué)性能，提出了多種黏結(jié)劑（Estane5703、PEG、PVDF、PCTFE 等）對(duì)PBX 結(jié)構(gòu)和力學(xué)性能的影響。Wang 等[18]通過(guò)MD 模擬研究了應(yīng)變率對(duì)RDX 基PBX 的動(dòng)態(tài)力學(xué)響應(yīng)和失效機(jī)制的影響，發(fā)現(xiàn)PBX 的初始損傷主要發(fā)生在黏結(jié)劑區(qū)域。Lv 等[19]使用MD 模擬研究了溫度和應(yīng)變率對(duì)TATB-F2314 界面力學(xué)性能及微觀結(jié)構(gòu)的影響，發(fā)現(xiàn)黏結(jié)劑F2314 會(huì)擴(kuò)散進(jìn)入TATB 層，導(dǎo)致TATB-F2314 界面處出現(xiàn)“內(nèi)嵌狀”界面混合層，并且隨著溫度的升高，界面混合層厚度逐漸增大。目前，通過(guò)MD 方法，在原子尺度上對(duì)含能晶體-黏結(jié)劑界面的微觀變形機(jī)理和力學(xué)性能的研究較少，主要原因是難以建立高顆粒填充且能真實(shí)反映PBX 細(xì)觀特性的結(jié)構(gòu)模型。

本研究通過(guò)MATLAB 編程建立復(fù)雜的PBX 炸藥模型，即含氟聚合物F2311 包覆在多個(gè)以HMX為基的炸藥顆粒外形成HMX-F2311 界面模型，采用MD 方法模擬單軸拉伸下HMX-F2311 炸藥的力學(xué)性能，通過(guò)分析HMX-F2311 的拉伸應(yīng)力-應(yīng)變曲線、界面微觀結(jié)構(gòu)演變和能量分布，探討應(yīng)變率對(duì)PBX 力學(xué)性能的影響。

1 力場(chǎng)和模型的建立

1.1 力 場(chǎng)

HMX-F2311 的力場(chǎng)由HMX 勢(shì)、F2311 勢(shì)以及HMX 與F2311 之間的界面相互作用勢(shì)組成，因此HMX-F2311 系統(tǒng)的總能量可寫(xiě)為

1.2 初始模型

首先，使用MATLAB 程序建立F2311 包覆HMX 的PBX 界面初始模型，如圖1(a)所示。模型中x、y方向各有4 個(gè)HMX 晶體顆粒，z方向有2 個(gè)HMX 晶體顆粒，黏結(jié)劑F2311 均勻包裹在顆粒的邊緣。模型尺寸約為200 ?×200 ?×100 ?，原子數(shù)約為280 000。

HMX-F2311 界面模型的優(yōu)化過(guò)程分為4 個(gè)步驟：能量最小化、優(yōu)化原子坐標(biāo)、退火消除殘余應(yīng)力以及獲得平衡構(gòu)象。首先，對(duì)HMX-F2311 界面模型進(jìn)行幾何弛豫，采用力收斂公差為10-8kcal/(mol·?)的共軛梯度法進(jìn)行能量最小化。其次，在NPT系綜下進(jìn)行MD 模擬，以優(yōu)化原子坐標(biāo)，壓力為標(biāo)準(zhǔn)大氣壓，溫度為300 K，步長(zhǎng)為0.1 fs，模擬時(shí)間為40 ps。接著，為降低系統(tǒng)的能量，消除殘余應(yīng)力，對(duì)模型進(jìn)行退火處理。將體系溫度迅速上升至350 K，弛豫50 ps 后，降至300 K，繼續(xù)弛豫50 ps。退火過(guò)程在NPT系綜下進(jìn)行，模型尺寸將發(fā)生明顯變化。最后，采用NVT系綜對(duì)體系進(jìn)行結(jié)構(gòu)優(yōu)化，模擬時(shí)間為20 ps。優(yōu)化后的HMX-F2311 界面結(jié)構(gòu)見(jiàn)圖1(b)。NVT優(yōu)化過(guò)程中的壓力和能量隨時(shí)間的變化曲線如圖2 所示。系統(tǒng)的壓力和能量最終趨于穩(wěn)定，上下波動(dòng)小于5%，表明此時(shí)HMX-F2311 結(jié)構(gòu)達(dá)到平衡。將模擬得到的平衡結(jié)構(gòu)作為初始構(gòu)型，進(jìn)行不同應(yīng)變率下單軸拉伸的MD 模擬。

圖1 HMX-F2311 界面的初始模型(a)以及優(yōu)化后HMX-F2311 界面的穩(wěn)定結(jié)構(gòu)(b)Fig. 1 Initial model of HMX-F2311 interface (a) and the stable structure of HMX-F2311 interface after optimization (b)

圖2 優(yōu)化過(guò)程中壓力和能量隨時(shí)間的變化Fig. 2 Pressure and energy fluctuations vs. time during optimization

1.3 拉伸加載

2 動(dòng)態(tài)拉伸下應(yīng)變率對(duì)PBX 界面的影響

為了研究應(yīng)變率對(duì)HMX-F2311 界面拉伸變形的影響，在應(yīng)變率為0.5×1010、1.0×1010、2.5×1010和5.0×1010s-1下進(jìn)行了一系列動(dòng)態(tài)拉伸模擬，應(yīng)力-應(yīng)變曲線如圖3(a)所示。不同應(yīng)變率下的拉伸應(yīng)力-應(yīng)變曲線表現(xiàn)出相似的特征，均可分為3 個(gè)階段：彈性階段（Ⅰ）、塑性階段（Ⅱ）和破壞階段（Ⅲ）。在拉伸加載初期，應(yīng)力隨應(yīng)變的增加而線性增加，即彈性階段（ ε=0.025）；在應(yīng)力-應(yīng)變曲線中沒(méi)有明顯的屈服點(diǎn)，應(yīng)力持續(xù)非線性增加至PBX 的拉伸強(qiáng)度（ ε=0.080），即塑性階段；達(dá)到拉伸強(qiáng)度后，隨著應(yīng)變的增大，應(yīng)力持續(xù)減小，PBX 的力學(xué)性能逐步劣化失效，為破壞階段。在模擬中，采用應(yīng)力-應(yīng)變曲線中的極值應(yīng)力作為PBX 的拉伸強(qiáng)度，并將彈性階段的應(yīng)力-應(yīng)變曲線斜率作為PBX 的彈性模量。如圖3(b)和圖3(c)所示，拉伸強(qiáng)度和彈性模量均隨拉伸應(yīng)變率的增加而增大，隨對(duì)數(shù)應(yīng)變率的增加而線性增大。

圖3 HMX-F2311 在不同應(yīng)變率下的應(yīng)力-應(yīng)變曲線 (a)以及拉伸強(qiáng)度 (b) 和彈性模量 (c) 與對(duì)數(shù)應(yīng)變率的關(guān)系Fig. 3 (a) Stress-strain curves of HMX-F2311 under different tension loading; (b) tension strength vs.logarithmic strain rate; (c) elastic modulus vs. logarithmic strain rate

為了比較不同應(yīng)變率（0.5×1010、1.0×1010、2.5×1010和5.0×1010s-1）下PBX 的微觀結(jié)構(gòu)演變，通過(guò)OVITO 可視化程序計(jì)算了拉伸過(guò)程中每個(gè)原子的局部應(yīng)變分布，如圖4 所示。在低應(yīng)變率（0.5×1010、1.0×1010s-1）拉伸加載下：當(dāng) ε =0.05 時(shí)，靠近黏結(jié)劑區(qū)域的應(yīng)變明顯高于炸藥晶體的內(nèi)部應(yīng)變，表明在拉伸加載下原子應(yīng)變主要集中在較軟的黏結(jié)劑F2311 上；當(dāng)應(yīng)力達(dá)到拉伸強(qiáng)度后（ ε=0.10），炸藥晶體與黏結(jié)劑F2311 界面處出現(xiàn)了部分孔洞和微裂紋；當(dāng) ε =0.15 時(shí)，靠近加載側(cè)的界面逐漸開(kāi)裂，出現(xiàn)黏結(jié)劑脫粘現(xiàn)象；隨著拉伸的持續(xù)，黏結(jié)劑的脫粘降低了PBX 的承載能力，導(dǎo)致材料的力學(xué)性能逐漸劣化；當(dāng) ε=0.25時(shí)，PBX 的拉應(yīng)力降到最低，在靠近加載側(cè)的界面處沿垂直于加載方向產(chǎn)生了一條主裂紋。PBX 的失效模式為黏結(jié)劑脫粘。

圖4 不同拉伸加載應(yīng)變率下HMX-F2311 界面在x 方向的應(yīng)變分布Fig. 4 Distribution of strain field in the x direction of HMX-F2311 interface under different tension strain rates

如圖5 所示，相比于低應(yīng)變率拉伸加載，在高應(yīng)變率（2.5×1010和5.0×1010s-1）拉伸加載下，PBX 微觀結(jié)構(gòu)演變整體上差異不大。加載初期，PBX 的拉伸原子形變主要集中在軟黏結(jié)劑附近。隨著拉伸的持續(xù)，相對(duì)較弱的HMX-F2311 界面發(fā)生脫粘。不同的是，在低應(yīng)變率拉伸下形成了一條大致垂直于加載方向的主裂紋，且斷裂位置靠近加載面一側(cè)。隨著拉伸應(yīng)變率的增加，在遠(yuǎn)離加載側(cè)的界面處也出現(xiàn)了大量微裂紋，破壞路徑分布在整個(gè)模型上。當(dāng) ε =0.20 時(shí)，低應(yīng)變率拉伸加載下HMX-F2311 靠近加載側(cè)的界面處已經(jīng)產(chǎn)生了一條裂紋，而高應(yīng)變率加載下HMX-F2311 界面處只存在孔洞，裂紋并未開(kāi)始擴(kuò)展，表明隨著應(yīng)變率的增大，HMX-F2311 的抗拉性能增強(qiáng)。HMX-F2311 在單軸動(dòng)態(tài)拉伸下的斷裂失效是由于黏結(jié)劑脫粘，與Rae 等[7]通過(guò)顯微鏡觀察到的斷裂方式相似。

圖5 不同拉伸應(yīng)變率下HMX-F2311 界面在x 方向的應(yīng)變分布Fig. 5 Distribution of strain field in the x direction of HMX-F2311 interface under different tension strain rates

為了從能量角度理解PBX 在單軸拉伸過(guò)程中的形變機(jī)制，分析了不同拉伸應(yīng)變率（0.5×1010和5.0×1010s-1）下PBX 界面的勢(shì)能隨應(yīng)變的變化情況，如圖6 所示。在低應(yīng)變率（0.5×1010s-1）拉伸加載下：當(dāng)PBX 拉伸處于彈性階段和塑性階段時(shí)，勢(shì)能隨著應(yīng)變的增大逐漸增大；當(dāng)達(dá)到HMX-F2311 的拉伸強(qiáng)度時(shí)，勢(shì)能也達(dá)到極值（0.68×104kcal/mol），隨后迅速下降。然而，在高應(yīng)變率（5.0×1010s-1）拉伸加載下，當(dāng)達(dá)到HMX-F2311 的拉伸強(qiáng)度后，勢(shì)能仍然持續(xù)升高至極值（3.26×104kcal/mol），隨后下降。隨著拉伸應(yīng)變率的增加，HMX-F2311 的勢(shì)能增加了4 倍。

圖6 不同應(yīng)變率下總勢(shì)能隨應(yīng)變的變化Fig. 6 Variations of potential energy with stain under different strain rates

模擬體系的勢(shì)能一般分為成鍵相互作用和非鍵相互作用，其中：成鍵相互作用主要包括鍵能、角能、二面角能以及非正常二面角扭轉(zhuǎn)能等，而非鍵相互作用包括范德華力作用和庫(kù)侖相互作用。拉伸應(yīng)變率為0.5×1010和5.0×1010s-1時(shí)PBX 界面各勢(shì)能分項(xiàng)隨應(yīng)變的變化曲線如圖7 所示。為了便于分析，勢(shì)能以及各能量分項(xiàng)均減去未拉伸時(shí)的初始值，其中：ΔEpotential、ΔEbond、ΔEangle、ΔEdihed、ΔEvdwl、ΔEimp分別為勢(shì)能、鍵能、角能、二面角能、范德華力相互作用能、非正常二面角扭轉(zhuǎn)能的變化量，ΔEvdwl、ΔEbond、ΔEangle、ΔEdihed以及ΔEimp均由LAMMPS 直接輸出。從圖7 可以看出：?jiǎn)屋S拉伸過(guò)程中影響勢(shì)能變化的主要因素是范德華力相互作用，尤其在5×1010s-1的高應(yīng)變率下，范德華力相互作用對(duì)勢(shì)能的變化起決定性作用；成鍵能對(duì)勢(shì)能的貢獻(xiàn)較小，其值隨應(yīng)變率的變化不明顯；鍵能和角能在拉伸過(guò)程中對(duì)勢(shì)能有一定的影響，尤其是彈性和塑性階段；二面角能和非正常二面角扭轉(zhuǎn)能在整個(gè)拉伸過(guò)程中幾乎沒(méi)有變化。

圖7 不同應(yīng)變率下各個(gè)勢(shì)能分項(xiàng)隨應(yīng)變的變化Fig. 7 Variations of each energy sub-item of potential energy with strain under different strain rates

PBX 的勢(shì)能隨拉伸應(yīng)變率的增加而快速增大的原因可能與黏結(jié)劑F2311 的構(gòu)象有關(guān)。在較低的拉伸應(yīng)變率下，F(xiàn)2311 鏈有更多的時(shí)間進(jìn)行平衡，并在每一步拉伸后進(jìn)行滑動(dòng)、旋轉(zhuǎn)和展開(kāi)以達(dá)到平衡，從而減少界面變形的約束。然而，在高應(yīng)變率下，F(xiàn)2311 鏈改變?cè)尤∠虻臅r(shí)間較短，導(dǎo)致HMXF2311 界面處的非平衡狀態(tài)和更大的變形阻力，從而需要較大的應(yīng)力才能使PBX 界面進(jìn)一步形變。由此可知，低應(yīng)變率更有利于黏結(jié)劑F2311 達(dá)到構(gòu)象平衡。

3 結(jié) 論

采用MD 方法模擬了HMX-F2311 在不同拉伸速率下的力學(xué)性能和斷裂機(jī)理，得到如下結(jié)論。

(1) HMX-F2311 的力學(xué)性能表現(xiàn)出高應(yīng)變率依賴性。拉伸強(qiáng)度和彈性模量隨著應(yīng)變率的增加而增大，并與對(duì)數(shù)應(yīng)變率呈線性關(guān)系。

(2) 在單軸拉伸加載下，應(yīng)變主要集中在黏結(jié)劑F2311 附近，并逐漸在界面處形成孔洞和裂紋。PBX 的斷裂是由于黏結(jié)劑F2311 的脫粘。在低應(yīng)變率下形成了一條垂直于加載方向的主裂紋，隨著拉伸應(yīng)變率的增加，破壞路徑將分布在整個(gè)模型上。

(3) 在單軸拉伸加載過(guò)程中，HMX-F2311 的勢(shì)能隨拉伸應(yīng)變率的增加而迅速增大。在高應(yīng)變率下，范德華力相互作用對(duì)拉伸過(guò)程中勢(shì)能的演變起決定作用。

研究結(jié)果對(duì)于理解PBX 界面的微觀結(jié)構(gòu)演變及失效機(jī)理具有重要意義，也有助于通過(guò)建立材料微觀結(jié)構(gòu)與性能之間的關(guān)系設(shè)計(jì)和改進(jìn)PBX 復(fù)合材料。