分子動(dòng)力學(xué)法研究摻雜缺陷對(duì)HMX/NQ 共晶炸藥性能的影響

杭貴云，余文力，王 濤，王金濤，苗 爽

(火箭軍工程大學(xué)，陜西 西安 710025)

引 言

HMX是一種綜合性能優(yōu)異的硝胺類高能炸藥，但因其感度較高，導(dǎo)致其應(yīng)用受到限制。硝基胍(NQ)是一種鈍感炸藥，安全性較好，但威力與能量密度較低。丁雄[1]研究了HMX/NQ共晶炸藥的性能，結(jié)果表明，共晶炸藥的感度低于HMX，且共晶具有較高的能量密度。因此，HMX/NQ共晶炸藥有望成為一種新型的高能炸藥。

HMX通常采用改進(jìn)乙酸酐法制備，在產(chǎn)物中會(huì)存在雜質(zhì)成分RDX。因此，在HMX/NQ共晶炸藥制備過程中，RDX也會(huì)進(jìn)入其中，導(dǎo)致晶體中存在摻雜缺陷，可能會(huì)影響炸藥的性能，如穩(wěn)定性、感度、能量特性和力學(xué)性能等，從而進(jìn)一步影響武器彈藥的安全性與作戰(zhàn)性能。因此，研究晶體缺陷對(duì)炸藥性能的影響具有重要意義。LaBarbera等[2]研究了缺陷對(duì)含能材料“熱點(diǎn)”形成的影響；Luscher等[3]研究了晶體缺陷對(duì)α-RDX感度的影響；Xue等[4]研究了位錯(cuò)缺陷對(duì)RDX沖擊感度的影響；彭亞晶等[5]研究了空位缺陷對(duì)RDX性能的影響。以往研究的晶體缺陷類型主要側(cè)重于空位與位錯(cuò)，而研究摻雜缺陷的報(bào)道相對(duì)較少。

本研究分別建立了“完美”型與含有不同摻雜率缺陷的HMX/NQ共晶炸藥模型。在Materials Studio 7.0(簡(jiǎn)稱MS)軟件中[6]，采用分子動(dòng)力學(xué)方法預(yù)測(cè)了模型的穩(wěn)定性、感度、能量特性和力學(xué)性能，綜合評(píng)估了晶體缺陷對(duì)炸藥性能的影響，以期為炸藥的性能評(píng)估提供理論參考。

1 計(jì)算模型與計(jì)算方法

1.1 單個(gè)分子的建立

HMX/NQ共晶炸藥由HMX與NQ組成，摩爾比為1∶1。在MS軟件中，分別建立HMX與NQ的單個(gè)分子模型，如圖1所示。

圖1 HMX與NQ的分子模型Fig.1 Molecule models of HMX and NQ

1.2 HMX/NQ共晶炸藥原始模型的建立

本研究中使用的HMX/NQ共晶炸藥結(jié)構(gòu)為文獻(xiàn)[1]給出的理論預(yù)測(cè)結(jié)構(gòu)。HMX/NQ共晶炸藥屬于單斜晶系，空間群為P21，晶格參數(shù)為a=2.264nm，b=0.420nm，c=0.798nm，α=90.00°，β=103.00°，γ=90.00°，單個(gè)晶胞中包含2個(gè)HMX與2個(gè)NQ分子[1]。HMX/NQ共晶炸藥的單個(gè)晶胞模型如圖2(a)所示，而后將單個(gè)晶胞模型擴(kuò)展為36(3×4×3)的超晶胞模型，其中包含72個(gè)HMX分子與72個(gè)NQ分子，一共144個(gè)分子，如圖2(b)所示。為了便于與含有晶體缺陷的模型進(jìn)行比較，將“完美”型晶體模型標(biāo)記為Model-1。

1.3 HMX/NQ共晶炸藥缺陷模型的建立

用2個(gè)RDX分子替換“完美”型晶體中的2個(gè)HMX分子(圖3(a)中標(biāo)記為黃色)，得到摻雜率為1.39%的缺陷晶體模型，如圖3(b)所示，將所得缺陷模型標(biāo)記為Model-2。

圖3 “完美”型晶體模型(Model-1)與摻雜率為1.39%的 缺陷晶體模型(Model-2)Fig.3 Perfect crystal model(Model-1) and defective crystal model with adulteration ratio of 1.39%(Model-2)

在Model-2的基礎(chǔ)上，繼續(xù)用2個(gè)RDX分子取代模型中的2個(gè)HMX分子，得到摻雜率為2.78%的缺陷晶體模型，標(biāo)記為Model-3。

采用同樣的方法，在Model-3的基礎(chǔ)上，繼續(xù)用4個(gè)RDX分子替換4個(gè)HMX分子，得到摻雜率為5.56%的缺陷晶體模型，標(biāo)記為Model-4。

類似地，分別用12個(gè)、16個(gè)、24個(gè)RDX分子替換初始模型中同等數(shù)量的HMX分子，得到摻雜率分別為8.34%、11.11%和16.67%的缺陷晶體模型，分別標(biāo)記為Model-5、Model-6和Model-7。

1.4 計(jì)算條件設(shè)置

分別對(duì)HMX/NQ共晶炸藥的“完美”型晶體模型和含有摻雜缺陷的晶體模型進(jìn)行能量最小化處理，對(duì)其結(jié)構(gòu)進(jìn)行優(yōu)化，而后進(jìn)行分子動(dòng)力學(xué)計(jì)算，其中溫度設(shè)置為295K，壓強(qiáng)設(shè)置為0.0001GPa，選擇NPT系綜與COMPASS力場(chǎng)[7-8]。選擇COMPASS力場(chǎng)是因?yàn)樵摿?chǎng)采用從頭算法，適用于凝聚相和不同類型物質(zhì)相互作用的研究，能在較大范圍內(nèi)對(duì)處于孤立體系和凝聚態(tài)體系的多種物質(zhì)的性能進(jìn)行準(zhǔn)確預(yù)測(cè)。溫度采用Andersen控溫方法[9]，壓強(qiáng)采用Parrinello控壓方法[10]，范德華力(vdW)的計(jì)算采用Atom-based方法，靜電作用的計(jì)算采用Ewald方法，時(shí)間步長(zhǎng)為1fs，總模擬步數(shù)為2×105步，其中前105步用于熱力學(xué)平衡計(jì)算，后105步用于統(tǒng)計(jì)分析。模擬過程中，每103fs保存一次軌跡，共得100幀軌跡文件。

2 結(jié)果與討論

2.1 平衡判別與平衡體系

在提取數(shù)值模擬計(jì)算結(jié)果時(shí)，需要讓混合體系達(dá)到熱力學(xué)平衡狀態(tài)，此時(shí)必須同時(shí)滿足溫度平衡和能量平衡。通常認(rèn)為當(dāng)溫度與能量波動(dòng)范圍為5%～10%時(shí)，體系已經(jīng)達(dá)到平衡狀態(tài)。以Model-4為例，模擬過程中混合體系的溫度變化曲線如圖4(a)所示，能量變化曲線如圖4(b)所示。

從圖4可以看出，模擬計(jì)算初期，溫度和能量均有所上升，且波動(dòng)幅度較大。隨著時(shí)間的推移，溫度和能量的波動(dòng)幅度逐漸減小，最終溫度波動(dòng)幅度約為±15K，能量波動(dòng)幅度約為±5%，偏差相對(duì)較小，表明混合體系已達(dá)到熱力學(xué)平衡狀態(tài)。對(duì)于其他的晶體模型，分子動(dòng)力學(xué)計(jì)算時(shí)，均以溫度平衡和能量平衡來判別混合體系是否達(dá)到平衡狀態(tài)。

圖4 混合體系的溫度和能量隨時(shí)間的變化曲線Fig.4 Change curves of temperature and energy of mixed system with time

2.2 結(jié)合能

對(duì)于HMX/NQ共晶及其摻雜缺陷晶體，結(jié)合能(Eb)的計(jì)算公式如下：

Eb=-Einter=-[Etotal-(ENQ+Eother)]

(1)

式中：Eb為結(jié)合能；Einter為分子之間的相互作用力；Etotal為混合體系達(dá)到平衡狀態(tài)時(shí)對(duì)應(yīng)的總能量；ENQ為去掉體系中的其他組分后，NQ分子對(duì)應(yīng)的總能量；Eother為去掉NQ分子后，體系中的HMX與RDX分子對(duì)應(yīng)的總能量。上述單位均為kJ/mol。

計(jì)算得到Model-1～Model-7的結(jié)合能分別為335.7、331.4、320.6、318.3、315.5、310.8、298.6kJ/mol。從以上數(shù)據(jù)可以看出，“完美”型晶體(Model-1)的結(jié)合能最大，表明HMX與NQ分子間的相互作用力最強(qiáng)，炸藥的穩(wěn)定性最好。缺陷晶體的結(jié)合能均有不同程度減小，且隨摻雜率的增加，結(jié)合能逐漸減小。當(dāng)摻雜率為1.39%(Model-2)時(shí)，結(jié)合能為331.4kJ/mol；當(dāng)摻雜率為16.67%(Model-7)時(shí)，結(jié)合能為298.6kJ/mol，結(jié)合能減小幅度為1.28%～11.05%。結(jié)合能減小，說明分子之間相互作用力的強(qiáng)度減弱，炸藥的穩(wěn)定性降低。結(jié)合能減小的原因可能是由于晶體缺陷的影響，炸藥的結(jié)構(gòu)遭到破壞，分子的排列方式發(fā)生變化，導(dǎo)致分子間的相互作用力減弱，結(jié)合能減小。

2.3 感度

依據(jù)“熱點(diǎn)”理論[11]與“引發(fā)鍵”思想[12]，同時(shí)參考以往的研究工作[13-16]，本研究選用引發(fā)鍵鍵長(zhǎng)、引發(fā)鍵鍵連雙原子作用能和內(nèi)聚能密度作為判據(jù)，從而預(yù)測(cè)炸藥的感度大小。

2.3.1 引發(fā)鍵鍵長(zhǎng)

在HMX/NQ共晶及其缺陷晶體中，HMX所占的比例最大，對(duì)外界的刺激最敏感。當(dāng)受到外界刺激時(shí)，HMX更容易發(fā)生化學(xué)反應(yīng)，使共晶炸藥發(fā)生分解或爆炸。HMX的引發(fā)鍵為分子中N-NO2基團(tuán)上的N—N鍵[12,17-18]，因此選擇HMX分子中的N—NO2鍵作為引發(fā)鍵來預(yù)測(cè)體系的感度。

以Model-2為例，圖5給出了平衡狀態(tài)時(shí)，體系中引發(fā)鍵(N—NO2鍵)的鍵長(zhǎng)分布情況。表1給出了平衡狀態(tài)下，不同模型中引發(fā)鍵的最可幾鍵長(zhǎng)(Lprob)、平均鍵長(zhǎng)(Lave)和最大鍵長(zhǎng)(Lmax)。

圖5 引發(fā)鍵的鍵長(zhǎng)分布(Model-2)Fig.5 Trigger bond length distribution of Model-2

ModelLprob/nmLave/nmLmax/nm10.13830.13820.153220.13820.13830.153930.13840.13840.154540.13830.13830.155950.13840.13830.157860.13830.13840.159470.13850.13840.1613

從圖5可以看出，當(dāng)混合體系達(dá)到平衡狀態(tài)時(shí)，引發(fā)鍵(N—NO2鍵)的鍵長(zhǎng)分布呈近似對(duì)稱的高斯分布。從表1可以看出，對(duì)于不同的模型，最可幾鍵長(zhǎng)與平均鍵長(zhǎng)近似相等，并且變化范圍很小，表明晶體缺陷對(duì)最可幾鍵長(zhǎng)與平均鍵長(zhǎng)的影響很小，而最大鍵長(zhǎng)變化很明顯，且不同模型之間的差異較大。對(duì)于“完美”型晶體模型(Model-1)，最大鍵長(zhǎng)為0.1532nm，對(duì)于摻雜率分別為1.39%(Model-2)、16.67%(Model-7)的晶體模型，最大鍵長(zhǎng)分別為0.1539、0.1613nm，增大幅度為0.46%～5.29%。最大鍵長(zhǎng)增大，說明引發(fā)鍵的強(qiáng)度減小，預(yù)示含能材料的感度升高，安全性降低。此外，從表1中還可以看出，隨著摻雜率的增加，最大鍵長(zhǎng)逐漸增大，表明炸藥的感度隨著缺陷數(shù)量的增加而逐漸升高。

2.3.2 引發(fā)鍵鍵連雙原子作用能

引發(fā)鍵鍵連雙原子作用能(EN-N)的計(jì)算公式如下[13-16]：

(2)

式中：E1為混合體系達(dá)到熱力學(xué)平衡狀態(tài)時(shí)對(duì)應(yīng)的總能量(kJ/mol)；E2為固定HMX中所有的N原子后，混合體系達(dá)到平衡狀態(tài)時(shí)的總能量(kJ/mol)；n為體系中HMX分子中包含的N—NO2鍵的數(shù)量。

通過計(jì)算，得到Model-1～Model-7的鍵連雙原子作用能分別為149.75、148.81、146.52、142.04、139.18、132.5、123.94kJ/mol。從以上數(shù)據(jù)可以看出，“完美”型晶體(Model-1)的鍵連雙原子作用能最大，而缺陷晶體的鍵連雙原子作用能均有不同程度的減小。其中，當(dāng)摻雜率為1.39%(Model-2)時(shí)，鍵連雙原子作用能為148.81kJ/mol；當(dāng)摻雜率為16.67%(Model-7)時(shí)，鍵連雙原子作用能為123.94kJ/mol；鍵連雙原子作用能減小的幅度為0.63%～17.24%。鍵連雙原子作用能減小，說明引發(fā)鍵的強(qiáng)度減弱，鍵斷裂時(shí)需要的能量降低，即含能材料的感度升高、安全性降低。此外，鍵連雙原子作用能的數(shù)據(jù)也表明，隨著摻雜率的增加，即缺陷數(shù)量的增加，鍵連雙原子作用能逐漸減小，預(yù)示炸藥的感度逐漸升高，安全性逐漸減弱。因此，晶體缺陷會(huì)對(duì)炸藥的安全性產(chǎn)生不利影響。在缺陷晶體中，鍵連雙原子作用能減小的原因可能是由于晶體結(jié)構(gòu)與分子排列方式發(fā)生改變，NQ的鈍感效應(yīng)降低，HMX的活性增強(qiáng)，從而使引發(fā)鍵鍵能減小、感度升高。

2.3.3 內(nèi)聚能密度

內(nèi)聚能密度(CED)的計(jì)算公式如下：

(3)

式中：CED為內(nèi)聚能密度，kJ/cm3；HV為摩爾蒸發(fā)熱，kJ/mol；RT為物質(zhì)氣化時(shí)所作的膨脹功，kJ/mol；Vm為摩爾體積，cm3/mol。

通過分子動(dòng)力學(xué)模擬，得到不同模型的內(nèi)聚能密度、范德華力與靜電力，結(jié)果如表2所示。

從表2可以看出，“完美”型晶體的內(nèi)聚能密度、范德華力與靜電力最大。在缺陷晶體中，內(nèi)聚能密度的最大值與最小值分別為0.832和0.748kJ/cm3；與“完美”型晶體相比，內(nèi)聚能密度減小幅度為0.83%～10.85%。內(nèi)聚能密度減小，表明炸藥由凝聚態(tài)變?yōu)闅鈶B(tài)時(shí)吸收的能量減小，預(yù)示炸藥的感度升高、安全性降低。此外，隨著摻雜率的增加，內(nèi)聚能密度逐漸減小，表明炸藥的感度隨著晶體缺陷的增加而逐漸升高。內(nèi)聚能密度減小的原因可能是晶體結(jié)構(gòu)遭到破壞，分子排列方式發(fā)生改變，從而導(dǎo)致分子間的范德華力與靜電力作用強(qiáng)度減弱，內(nèi)聚能密度減小、感度升高。

表2 不同模型的內(nèi)聚能密度、范德華力與靜電力

2.4 爆轟性能

本研究采用修正氮當(dāng)量法[19]計(jì)算炸藥的爆轟參數(shù)并預(yù)測(cè)其能量特性。

對(duì)于由C-H-O-N 4種元素組成、分子式為CaHbOcNd的炸藥，氧平衡系數(shù)(OB)的計(jì)算公式如下：

(4)

式中：a、b、c分別為炸藥分子中包含的C原子、H原子與O原子的數(shù)目；Mr為炸藥的相對(duì)分子質(zhì)量。

對(duì)于混合炸藥，氧平衡系數(shù)的計(jì)算公式如下：

OB=∑wiOBi

(5)

式中：wi為混合炸藥中第i種組分所占的質(zhì)量分?jǐn)?shù)，%；OBi為第i種組分對(duì)應(yīng)的氧平衡系數(shù)。

爆轟參數(shù)(D、p)的計(jì)算公式如下：

(6)

式中：D為炸藥的爆速，m/s；p為炸藥的爆壓，GPa；ρ為炸藥的密度，g/cm3；∑Nch為炸藥的修正氮當(dāng)量；pi為1 mol炸藥爆炸時(shí)生成第i種爆轟產(chǎn)物的摩爾數(shù)；Npi為第i種爆轟產(chǎn)物的氮當(dāng)量系數(shù)；BK為炸藥分子中第K種化學(xué)鍵出現(xiàn)的次數(shù)；NBK為炸藥分子中第K種化學(xué)鍵的氮當(dāng)量系數(shù)；Gj為炸藥分子中第j種基團(tuán)出現(xiàn)的次數(shù)；NGj為炸藥分子中第j種基團(tuán)的氮當(dāng)量系數(shù)。公式(6)中各變量的計(jì)算方法詳見參考文獻(xiàn)[20]。

根據(jù)修正氮當(dāng)量理論，通過計(jì)算得到不同模型的爆轟參數(shù)，結(jié)果如表3所示，其中炸藥的密度可以直接從分子動(dòng)力學(xué)結(jié)果中提取得到。

表3 不同模型的密度與爆轟參數(shù)

從表3可以看出，在不同的模型中，“完美”型晶體的密度、爆速與爆壓最大，分別為1.798g/cm3、8664m/s和34.00GPa，表明炸藥的能量密度最高。在缺陷晶體中，由于RDX的能量密度小于HMX，因此缺陷晶體的氧平衡系數(shù)、密度與爆轟參數(shù)均呈現(xiàn)出逐漸減小的變化趨勢(shì)且摻雜率越高，密度與爆轟參數(shù)越小。與“完美”型晶體相比，缺陷晶體的密度、爆速與爆壓減小幅度分別為0.89%～7.06%、0.68%～5.41%、1.85%～14.18%。密度與爆轟參數(shù)減小，表明炸藥的威力減小，能量密度降低，因此晶體缺陷會(huì)對(duì)炸藥的能量特性產(chǎn)生不利影響。

2.5 力學(xué)性能

力學(xué)參數(shù)可由彈性系數(shù)求得，滿足廣義胡克定律，其表達(dá)式為[21-22]：

σi=Cijεj

(7)

式中：σ為應(yīng)力；ε為應(yīng)變；Cij(i，j=1，2，…，6)為彈性系數(shù)，滿足Cij=Cji，因此獨(dú)立的彈性常數(shù)只有21個(gè)，對(duì)于完全的各向同性體，獨(dú)立的彈性常數(shù)只有2個(gè)(C11，C22)。

體積模量(K)與剪切模量(G)的計(jì)算公式如下：

KR=[S11+S22+S33+2(S12+S23+S31)]-1

(8)

GR= 15[4(S11+S22+S33)-4(S12+S23+S31)+

3(S44+S55+S66)]-1

(9)

式中：下標(biāo)R表示Reuss平均。柔量系數(shù)矩陣S=[Sij]等于彈性系數(shù)矩陣C的逆矩陣，即S=C-1。

力學(xué)參數(shù)之間存在如下關(guān)系：

E=2G(1+γ)=3K(1-2γ)

(10)

基于上述公式，可以求得拉伸模量(E)與泊松比(γ)的表達(dá)式：

(11)

(12)

通過計(jì)算，得到不同模型的彈性系數(shù)與力學(xué)性能參數(shù)，結(jié)果如表4與圖6所示。

表4 HMX/NQ共晶“完美”型模型與缺陷模型的力學(xué)性能

圖6 不同模型的力學(xué)參數(shù)Fig.6 Mechanical parameters of different models

從表4與圖6可看出，“完美”型模型(Model-1)的拉伸模量、體積模量與剪切模量的值最大，分別為13.389、8.225和5.490GPa，而柯西壓(C12-C44)的值最小(0.209GPa)，表明炸藥的剛性與硬度最強(qiáng)，而延展性較差。由于晶體缺陷的影響，E、K、G的值減小，而柯西壓增大。在缺陷晶體中，Model-2(摻雜率為1.39%)的E、K、G最大，分別為13.283、8.139和5.408GPa，柯西壓為0.317GPa；而Model-7(摻雜率為16.67%)對(duì)應(yīng)模量最小，分別為9.021、5.652和3.655GPa，柯西壓為1.996GPa。與“完美”型模型相比，E、K、G減小的幅度分別為0.106～4.368GPa、0.086～2.573GPa和0.082～1.835GPa，而柯西壓增大幅度為0.108～1.787GPa。模量減小，表明炸藥的剛性與硬度降低；柯西壓增大，表明炸藥的延展性增強(qiáng)。在缺陷晶體中，模量減小的原因可能是由于晶體結(jié)構(gòu)遭到破壞，分子排列方式發(fā)生了變化，從而使晶體的對(duì)稱性破壞，抵抗變形的能力減弱，在外界作用下，晶體更容易產(chǎn)生形變。

3 結(jié) 論

(1)與“完美”型晶體相比，缺陷晶體的結(jié)合能減小，分子之間的相互作用力減弱，炸藥的穩(wěn)定性降低。隨著缺陷數(shù)量的增加，結(jié)合能逐漸減小，炸藥的穩(wěn)定性逐漸減弱。

(2)缺陷晶體的引發(fā)鍵鍵長(zhǎng)增大，而鍵連雙原子作用能與內(nèi)聚能密度減小，表明炸藥的感度升高，安全性降低。

(3)由于RDX的能量密度小于HMX，因此缺陷晶體的密度與爆轟參數(shù)減小，威力與能量密度降低。隨著晶體缺陷數(shù)量的增加，炸藥的能量特性逐漸減弱。

(4)在缺陷晶體中，由于晶體結(jié)構(gòu)遭到破壞，炸藥的拉伸模量、體積模量與剪切模量減小，柯西壓增大，表明炸藥的剛性與硬度降低、延展性增強(qiáng)。在外界作用下，炸藥更容易產(chǎn)生形變。