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

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

引 言

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

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

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

1 計算模型與計算方法

1.1 單個分子的建立

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

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

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

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

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

用2個RDX分子替換“完美”型晶體中的2個HMX分子(圖3(a)中標記為黃色)，得到摻雜率為1.39%的缺陷晶體模型，如圖3(b)所示，將所得缺陷模型標記為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的基礎上，繼續(xù)用2個RDX分子取代模型中的2個HMX分子，得到摻雜率為2.78%的缺陷晶體模型，標記為Model-3。

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

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

1.4 計算條件設置

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

2 結果與討論

2.1 平衡判別與平衡體系

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

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

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

2.2 結合能

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

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

(1)

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

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

2.3 感度

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

2.3.1 引發(fā)鍵鍵長

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

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

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

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

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

(2)

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

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

2.3.3 內聚能密度

內聚能密度(CED)的計算公式如下：

(3)

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

通過分子動力學模擬，得到不同模型的內聚能密度、范德華力與靜電力，結果如表2所示。

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

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

2.4 爆轟性能

本研究采用修正氮當量法[19]計算炸藥的爆轟參數并預測其能量特性。

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

(4)

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

對于混合炸藥，氧平衡系數的計算公式如下：

OB=∑wiOBi

(5)

式中：wi為混合炸藥中第i種組分所占的質量分數，%；OBi為第i種組分對應的氧平衡系數。

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

(6)

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

根據修正氮當量理論，通過計算得到不同模型的爆轟參數，結果如表3所示，其中炸藥的密度可以直接從分子動力學結果中提取得到。

表3 不同模型的密度與爆轟參數

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

2.5 力學性能

力學參數可由彈性系數求得，滿足廣義胡克定律，其表達式為[21-22]：

σi=Cijεj

(7)

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

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

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)

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

力學參數之間存在如下關系：

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

(10)

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

(11)

(12)

通過計算，得到不同模型的彈性系數與力學性能參數，結果如表4與圖6所示。

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

圖6 不同模型的力學參數Fig.6 Mechanical parameters of different models

從表4與圖6可看出，“完美”型模型(Model-1)的拉伸模量、體積模量與剪切模量的值最大，分別為13.389、8.225和5.490GPa，而柯西壓(C12-C44)的值最小(0.209GPa)，表明炸藥的剛性與硬度最強，而延展性較差。由于晶體缺陷的影響，E、K、G的值減小，而柯西壓增大。在缺陷晶體中，Model-2(摻雜率為1.39%)的E、K、G最大，分別為13.283、8.139和5.408GPa，柯西壓為0.317GPa；而Model-7(摻雜率為16.67%)對應模量最小，分別為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。模量減小，表明炸藥的剛性與硬度降低；柯西壓增大，表明炸藥的延展性增強。在缺陷晶體中，模量減小的原因可能是由于晶體結構遭到破壞，分子排列方式發(fā)生了變化，從而使晶體的對稱性破壞，抵抗變形的能力減弱，在外界作用下，晶體更容易產生形變。

3 結 論

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

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

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

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