高壓下β-HMX熱分解機理的ReaxFF反應分子動力學模擬

周婷婷 石一丁 黃風雷

(北京理工大學,爆炸科學與技術國家重點實驗室,北京100081)

高壓下β-HMX熱分解機理的ReaxFF反應分子動力學模擬

周婷婷 石一丁 黃風雷*

(北京理工大學,爆炸科學與技術國家重點實驗室,北京100081)

采用ReaxFF反應分子動力學方法研究了不同壓縮態β-HMX晶體(ρ=1.89,2.11,2.22,2.46,2.80, 3.20 g·cm-3)在T=2500 K時的熱分解機理,分析了壓力對初級和次級化學反應速率的影響、高壓與低壓下初始分解機理的區別以及造成反應機理發生變化的原因.發現HMX的初始分解機理與壓力(或密度)相關.低壓下(ρ＜2.80 g·cm-3)以分子內反應為主,即N－NO2鍵的斷裂、HONO的生成以及分子主環的斷裂(C－N鍵的斷裂).高壓下(ρ≥2.80 g·cm-3)分子內反應被顯著地抑制,而分子間反應得到促進,生成了較多的O2、HO等小分子和大分子團簇.初始分解機理隨壓力的變化導致不同密度下的反應速率和勢能也有所不同.本文在原子水平對高壓下HMX反應機理的深入研究對于認識含能材料在極端條件下的起爆、化學反應的發展以及爆轟等具有重要意義.

HMX;熱分解; 壓力;ReaxFF;分子動力學

1 引言

1,3,5,7-四硝基-1,3,5,7-四氮雜環辛烷(HMX)等硝胺類炸藥因優越的爆轟性能而被廣泛應用于武器裝備及民用領域,這些炸藥通過復雜的化學反應釋放出巨大的能量.對化學反應過程的深入認識,如初始反應機理、詳細的反應路徑、反應產物的結構等,有利于建立宏觀的燃燒和爆轟模型.Lewis等1應用密度泛函理論(DFT)對氣相α-HMX分子分解過程的研究表明N－N鍵的斷裂、HONO的離解、CN鍵的斷裂及主環的分解是主要的反應路徑.他們認為N－N鍵斷裂生成NO2在單分子的初始分解中起著主要作用,而HONO的離解及C－N鍵的斷裂在凝聚態HMX的分解過程中更有可能發生. Chakraborty等2通過DFT計算提出了單分子β-和δ-HMX分解的三種機理:N－N鍵斷裂、HONO的離解、主環斷裂生成CH2N2O2并進一步分解成CH2O和N2O.最近,Sharia等3應用DFT和過渡態理論(TST)研究了單分子和小體系HMX晶體(16個分子)的初始分解機理,得到了N－N鍵斷裂、HONO離解及主環斷裂的反應熱和活化能.姜富靈等4采用DFT在B3LYP/6-31G(d)水平上研究了氣相β-HMX分子N－NO2鍵斷裂的焓變及自由能變.

對常壓下凝聚態HMX熱分解機理的研究有較多的實驗報導.研究結果普遍認為在分解過程中生成了H2CO、N2O、HONO、HCN等中間產物,并進一步分解為N2、H2O、CO和CO2等最終產物.5-8Brill5提出了HMX分解的兩種競爭機理,即HMX→4(HONO+HCN)和HMX→4(H2CO+N2O).Tang等6認為用多步反應來解釋凝聚態HMX的熱分解過程更為合適.Tarver等7提出了HMX熱分解的三步模型:N－N鍵和C－N鍵的斷裂生成初始產物;初始產物進一步分解為H2CO、N2O、HONO和HCN等中間產物;中間產物反應生成H2O、N2、CO和CO2等最終產物.

在壓力作用下,凝聚態含能材料的分解過程更加復雜.因為壓力會對含能材料的微觀結構產生影響,進而影響化學反應.Gilman8認為,炸藥在受到沖擊作用時,強烈的壓縮使炸藥密度增大,共價鍵彎曲,當這種作用達到臨界值時炸藥發生局域金屬化,導致電子激發.對共價分子,化學鍵的伸縮和彎曲使最高被占據分子軌道(HOMO)和最低空分子軌道(LUMO)之間的能隙降低.在分子晶體中,HOMO-LUMO帶隙的閉合使得價帶電子發生離域,導致化學反應的發生.9Margetis等10應用電荷自洽密度泛函緊束縛方法(SCC-DFTB)對靜水和單軸壓縮下凝聚態硝基甲烷(NM)能帶結構的研究表明,C－N鍵的彎曲是造成帶隙減小的主要原因,同時壓縮作用使C－H鍵拉伸并最終造成氫質子的離解.Manaa11應用第一性原理研究了三氨基三硝基苯(TATB)帶隙與硝基彎曲的關系,結果表明硝基彎曲會使帶隙閉合,而硝基的彎曲可能是由剪切作用造成.Lu等12用DFT研究了靜水壓縮下β-HMX的結構和振動特性,發現壓力顯著地減小了N－N鍵長,說明N－N鍵在初始分解中起著重要作用. Kuklja等13,14應用DFT對剪切作用下FOX-7能帶結構與態密度的研究發現,剪切作用使得帶隙及化學反應的能壘降低,說明剪切作用可能是引發FOX-7熱分解的主要原因.

在極端條件下(如高溫高壓),對凝聚態含能材料化學反應機理的研究因材料本身的性質、使用環境的特殊性以及反應的快速性而比較困難,特別是實驗方面的研究.量子力學(QM)方法雖然能夠準確地預測化學反應的過渡態、反應能壘、反應產物等,但主要用于研究單分子和小體系的靜態性質.基于QM的分子動力學(MD)方法能夠研究原子及分子結構的動力學過程,但因較大的計算量使得其在凝聚態含能材料特別是大體系上的應用受到限制. Manaa等15,16采用QM-MD方法研究了常壓下單晶胞δ-HMX(6個分子)在溫度為3500 K時的熱分解過程,結果表明N－NO2鍵斷裂是初始分解機理.最近,Zhu等17采用從頭算MD方法研究了小體系β-HMX晶體(4個分子)在沖擊作用下的化學過程,他們認為N－O鍵的斷裂和主環的分解在初始分解中占主導地位.

基于第一性原理的ReaxFF反應力場18的發展為研究大體系的物理化學性質提供了有效的途徑. ReaxFF反應分子動力學(ReaxFF-RMD)方法能夠在較短的計算時間內模擬含有上百萬個原子的體系在不同外加載荷下的物理化學過程,從原子和分子尺度提供詳細和準確的信息.在含能材料領域, ReaxFF-RMD已經被廣泛地用于研究熱、壓縮、沖擊及剪切作用下的物理化學問題,為含能材料的化學反應機理提供了非常有價值的信息.采用這一方法,Rom等19研究了在一定溫度范圍內不同壓縮態的液態NM的熱分解機理,結果表明NM的初始分解機理與壓縮量有關.Zhang等20研究了不同溫度和密度下HMX和TATB晶體的熱分解過程,得到了常壓下初級和次級反應的活化能以及HMX與TATB熱分解機理的主要區別.Strachan等21研究了不同密度環三亞甲基三硝胺(RDX)晶體的熱分解過程,得到了各產物數量隨密度的變化規律以及低密度與高密度晶體分解機理的主要區別.Strachan等22研究了RDX在不同沖擊速度下的化學反應,結果表明沖擊速度對RDX反應機理有較大影響.Zybin23、An24和Zhou25等分別研究了季戊四醇四硝酸酯(PETN)、RDX和HMX在壓縮-剪切作用下材料的沖擊各向異性,得到了與實驗一致的研究結果.

在壓縮或沖擊作用下,炸藥晶體被強烈壓縮,分子間距顯著減小,分子發生明顯變形,使得化學反應過程比常壓下的反應更為復雜,與氣態單分子的分解更加不同.因此,研究高壓下炸藥的分解機理對于深入認識炸藥在極端條件下的化學反應、沖擊起爆以及爆轟等具有重要意義.盡管應用分子動力學方法對凝聚態HMX在壓力作用下的化學反應已有報導,20但是卻不夠詳細和深入,還有較多問題值得研究,如壓力如何影響初始分解機理、反應路徑及反應速率;常壓和高壓下反應機理的主要區別;反應機理發生變化的原因等.本文將應用ReaxFF-RMD方法研究不同壓縮態β-HMX晶體的化學反應機理,以期解決以上問題.

2 模擬方法及說明

2.1 ReaxFF勢函數

ReaxFF反應力場以鍵級(BO)為核心.18鍵級是原子間距離的函數,在分子動力學的每一次循環時進行計算.當化學鍵斷裂時,與價鍵相關的能量和力變為零,因而該力場能夠描述化學鍵斷裂和生成.

ReaxFF勢函數的表達式如式(1)所示.20它分成以下幾個部分:基于鍵級的共價相互作用(Ebond、Eval和Etors)、原子間庫侖力(ECoulomb)、分子間范德華作用力(EvdWaals)、氫鍵(EH-bond)以及各修正項(Elp、Eover、Eunder、Epen和Econj)用以正確描述不同化學環境下化學鍵的斷裂和生成.

2.2 初始構象的構建

β-HMX晶胞來源于中子衍射晶體數據,26其分子及晶體結構如圖1所示.首先構建4×2×3的超晶胞模型,共含有48個β-HMX分子、1344個原子.從低壓到高壓選取6個研究體系,晶胞參數來源于金剛石壓腔實驗27中壓力p=0,2.5,4.6,10.6,26.0,42.6 GPa所對應的值,對應的密度分別為1.89、2.11、2.22、2.46、2.80和3.20 g·cm-3.

2.3 模擬過程說明

首先對模型進行幾何優化,以獲得合理的初始構型;對優化后的體系進行升溫,并確保升溫過程中晶體不分解;升溫至目標溫度T=2500 K后進行等溫等容分子動力學模擬(NVT-MD).采用Berendsen方法對溫度進行調節,使溫度在設定值附近波動.本文所有的計算均采用周期性邊界條件.模擬的時間步長為0.1 fs,模擬時間為20 ps.在20 ps內,主要分解產物已經全部出現且足以描述初級及次級化學反應以及壓力對反應機理的影響.每隔50 fs記錄一次原子坐標和速度以及每個原子對的鍵級,用以對分解產物進行分析.以BO=0.3作為產物生成與否的判據,當BO≥0.3時,認為化學鍵形成,產物生成.

圖1 β-HMX的分子結構及晶胞結構Fig.1 Molecular and crystal structures for β-HMX

2.4 反應速率分析方法

對不同壓縮態HMX晶體在熱分解過程中的反應速率分階段進行分析.

(1)初始分解反應(吸熱階段):采用一階衰減指數式(2)19對HMX分子數量隨時間的變化進行擬合,得到不同壓縮態HMX晶體熱分解的初始反應速率.

式中,N0為HMX分子的初始數量;t0為HMX開始分解的時間;k1為初級反應速率.

(2)次級化學反應(放熱階段):對勢能隨時間的變化的分析發現,勢能在反應的初始階段存在一個最大值,對應的時間設為tmax.當t＜tmax時,勢能因體系吸熱反應而升高;當t＞tmax時,勢能因體系放熱反應而降低,因而tmax將反應分為吸熱和放熱兩個階段.采用一階衰減指數式(3)19對放熱階段的勢能隨時間的變化進行擬合,得到不同壓縮態HMX晶體熱分解的次級反應速率.

式中,ΔUexo=U(tmax)-U∞,U∞為當t趨近于無窮大時勢能的平衡值;U(tmax)為勢能最大值;k2為次級反應速率.

(3)最終產物的生成.對不同壓縮態HMX晶體熱分解的最終產物隨時間的變化采用式(4)進行擬合,得到最終產物的生成速率.在本文的模擬時間范圍內(20 ps),最終產物還沒有完全達到穩定值,因而本文沒有對擬合得到產物的穩定值C∞進行討論.

式中,C∞為t趨近于無窮時產物的穩定值;ti為產物開始形成的時間;k3為產物的生成速率.

3 模擬結果與討論

3.1 反應速率

3.1.1 初始反應速率

不同壓縮態晶體中HMX分子數量隨時間的變化如圖2所示,HMX分子數量隨時間逐漸減少,在同一時刻,分子數量隨著密度的增加先增加后減少,分界點為ρ=2.80 g·cm-3.采用式(2)對HMX分子數量隨時間的變化進行指數擬合,得到初始分解速率k1隨密度的變化,如圖3所示.對ρ=1.89,2.11, 2.22,2.46,2.80,3.20 g·cm-3的體系,k1分別為8.13、6.80、5.13、2.04、7.52、26.32 ps-1,即當密度ρ＜2.80 g· cm-3時,k1隨密度的增加而減小,當密度ρ≥2.80 g· cm-3時,k1隨密度的繼續增加而增大.這說明HMX的初始分解速率與壓力有關,低壓會減慢HMX的初始分解速率,高壓則會加快初始熱分解.對ρ= 2.80,3.20 g·cm-3的高壓縮態體系,對應的體積壓縮量ΔV為32%和41%,即當ΔV≥32%時,HMX的初始分解速率反而隨密度的增加而加快,這與Rom等19對不同壓縮量的NM在T=2500-3500 K的溫度范圍內的熱分解的研究結果相似.他們發現當ΔV＜40%時,NM的初始分解速率降低;當ΔV≥40%時,初始分解速率反而加快.

圖2 T=2500 K時不同壓縮態β-HMX晶體中的HMX分子數量隨時間的變化Fig.2 Evolution of HMX molecule number for β-HMX crystals with different densities at T=2500 K

圖3 T=2500 K時不同壓縮態β-HMX晶體的初始反應速率k1和次級反應速率k2Fig.3 Initial reaction rates k1and the second reaction rates k2for β-HMX crystals with different densities at T=2500 K

不同壓縮態HMX晶體在分解過程中勢能隨時間的變化如圖4所示,晶體的勢能隨密度的增加逐漸增大.隨著密度的增加,晶體體積減小,原子間距和分子間距減小.原子間距的減少,分子內成鍵距離縮短,使得分子內原子間相互作用增強,勢能增加.分子間距的減小,分子間的平衡距離被打破,分子間的排斥力占據主導地位,使得體系能量升高,勢能增大.不同壓縮態HMX的勢能均隨時間先升高后降低,最大值對應的時間為tmax.當t＜tmax時,體系因初始分解反應而吸收能量,勢能升高;當t＞tmax時,體系因次級反應生成中間產物和最終產物而放出能量,勢能降低.對ρ=1.89,2.11,2.22,2.46,2.80, 3.20 g·cm-3的體系,tmax分別為0.1、0.15、1.1、1.85、1.6和0.05 ps,說明低壓下(ρ＜2.80 g·cm-3)tmax隨密度的增加而增大,吸熱階段延長;高壓下(ρ≥2.80 g·cm-3) tmax隨密度的繼續增加而減小,吸熱階段縮短.tmax隨密度的變化表明,當壓力較小時,晶體密度較小(ρ＜2.80 g·cm-3),體積較大,分子間還存在一定的可壓縮空間,因此吸熱過程可以持續一段時間.隨著壓力的增加,當晶體達到臨界結構(ρ≥2.80 g·cm-3)時,分子間距很難再被壓縮,此時體系熱容過小,不再能夠通過改變分子間距而吸收熱量,而是直接發生分子間的放熱反應,tmax變短.

圖4 T=2500 K時不同壓縮態β-HMX晶體中平均每個HMX分子的勢能隨時間的變化Fig.4 Evolution of potential energy per HMX molecule for β-HMX crystals with different densities at T=2500 K

HMX的初始分解速率k1和勢能達到最大值的時間tmax隨密度的變化表明,HMX的初始分解機理在高壓下(ρ≥2.80 g·cm-3)與低壓下(ρ＜2.80 g·cm-3)有所不同.常壓下,分子內N－N鍵的斷裂是最主要的初始分解機理;低壓下,分子內N－N鍵的斷裂仍然處于主導地位,但由于原子間距減小,鍵能增強, N－N鍵的斷裂變得困難,導致初始反應變慢,吸熱階段延長;高壓下,分子嚴重變形,分子間距及原子間距顯著減小,使得分子間或分子內相鄰非成鍵原子之間更容易接觸而發生反應,導致較多小分子或大分子團簇的生成而迅速放出熱量,吸熱階段縮短.Rom等19對不同壓縮態NM的熱分解的研究結果表明,在T≤3500 K時,低壓下NM的初始分解機理以單分子分解為主,而高壓下以分子間反應為主.他們將受激發分子的周圍分子比作熱浴,其作用是降低受激發分子的溫度;對單分子分解而言,增加密度會增強激發分子與熱浴之間的作用,使溫度降低,從而降低受激發分子的分解速率;對雙分子反應而言,增加密度會減少分子之間的空間,從而加快反應速率.因此,當溫度不是非常高時,增加密度對不同的反應類型會有相反的影響,這與本文的研究結果相同.壓力對HMX初始分解機理和反應路徑的具體影響見3.2節.

3.1.2 次級反應速率

采用式(3)對不同壓縮態HMX晶體的勢能隨時間的變化進行擬合,得到HMX分解過程中次級反應速率k2隨密度的變化,如圖3所示.對ρ=1.89, 2.11,2.22,2.46,2.80,3.20 g·cm-3的體系,k2分別為0.06、0.086、0.097、0.104、0.117和0.128 ps-1,即k2隨密度的增加而增大,表明壓力會加快次級反應速率.但k2并不是無限地增大,而是隨著密度的增加趨于一個極值.Rom等19對不同壓縮量NM熱分解的研究也發現NM中間產物的反應速率隨壓縮率的增加而增加.壓力減小了自由空間,使得初始分解產物更容易相互碰撞發生次級反應而生成中間產物和最終產物,次級反應速率加快.

3.2 產物分析

3.2.1 初始及中間反應產物

不同壓縮態HMX晶體在T=2500 K時的主要分解產物(NO2、HONO、NO3、HNO3、HNO、NO、H2O、N2、CO2等)隨時間的變化如圖5所示.常壓下,在HMX的熱分解過程中觀察到三種主要的初始分解機理:N－NO2鍵的斷裂、HONO的離解和HMX分子主環的斷裂(C－N鍵的斷裂),這與前人的研究結果1-3一致.NO2是最主要的初始分解產物,即在HMX的初始反應中,N－NO2鍵斷裂最容易發生. NO2的數量在達到最大值后迅速減小,這是由于次級反應消耗掉大量的NO2,生成NO3、HNO3等中間產物.HONO是僅次于NO2的初始分解產物,并進一步分解生成NO、HO、HNO等次級產物.HMX分子主環斷裂的主要產物為C2H2O2N2,并迅速分解成HCN、HCON、H2CO、CON等中間產物.相比于NO2和HONO的數量,由主環斷裂所生成的產物很少,說明主環斷裂在HMX的熱分解過程中不占主導地位,這與Sharia等3的研究結果相同.

Chakraborty等2采用B3LYP/6-31G(d)計算得到的氣相HMX分子中N－NO2鍵斷裂和HONO離解的能壘分別為166.52和186.61 kJ·mol-1.Sharia等3采用平面波密度泛函理論和廣義梯度近似(GGA)得到的這兩個反應的能壘分別為179.08和181.59 kJ·mol-1.他們考查了三種主環斷裂模式,3即(1) C4H8N8O8→2C2H4N4O4,并進一步分解為兩個CH2N2O2;(2)C4H8N8O8→CH2N2O2+open RDX;(3) C4H8N8O8→4CH2N2O2,對應的能壘分別為406.68、293.72和347.27 kJ·mol-1.3因此,N－N鍵斷裂和HONO離解因較低的反應能壘比C－N主環斷裂更容易發生.對凝聚態HMX晶體,Sharia等3認為主環斷裂因較高的能壘及較慢的反應速率而不太可能出現,因而他們只計算了N－NO2鍵斷裂和HONO離解的能壘,分別為200.41和218.82 kJ·mol-1.3所以N－NO2鍵斷裂和HONO離解在HMX晶體的初始分解中占主導地位,且N－NO2鍵斷裂最容易發生.

圖5 T=2500 K時不同壓縮態β-HMX晶體中平均每個HMX分子的主要分解產物隨時間的變化Fig.5 Evolution of products per HMX molecule for β-HMX crystals with different densities at T=2500 K

壓力對HMX晶體的初始分解機理會產生影響.從圖5可以看出,初始分解產物NO2和HONO的數量隨體系密度的增加而減少;高壓下有較多的小分子生成,如O2、HO等,且數量與NO2和HONO相當.下面將詳細分析密度對主要反應路徑的影響.

3.2.1.1 N－NO2鍵斷裂

壓力對初始反應路徑N－NO2鍵斷裂的影響如圖6(a)所示,NO2的數量隨著體系密度的增加急劇減少.對ρ=1.89,2.11,2.22,2.46,2.80,3.20 g·cm-3的體系,NO2的最大值分別為1.6、1.25、1.2、0.8、0.5和0.3;在t=1 ps時,單個HMX分子中N－NO2鍵斷裂的數量隨密度的增加而減小,如圖6(b)所示,說明壓力顯著地抑制了N－NO2鍵的斷裂.Strachan等22對RDX在不同沖擊速度下的分解機理的研究結果也表明當沖擊速度較大時(vimp＞8 km·s-1),NO2的數量反而隨沖擊速度的繼續增加而降低.N－N鍵長隨著體系密度的增大逐漸減小(如Supporting Information中表S1所示),鍵能增強,導致壓力作用下NNO2鍵的斷裂比較困難,NO2數量減少.Lu等12采用DFT對靜水壓縮下β-HMX的結構和振動特性的研究表明壓力顯著地減小了N－N鍵長.對反應路徑的分析發現,壓力作用下HMX分子中N－NO2鍵的斷裂是一個可逆過程,即C4H8N8O8?C4H8N7O6+ NO2,N－NO2鍵斷裂后形成的NO2很容易在較小的空間里再次與C4H8N7O6結合生成HMX,導致NO2數量的減少.且C4H8N7O6進一步分解產生NO2也更加困難,進一步造成NO2數量的減少.

NO3和HNO3的數量隨著體系密度的增加而減少,如圖6(c,d)所示.NO3主要從HMX分子脫落以及游離態NO2與O原子結合而成.在壓力作用下, NO2數量急劇減少,導致它與O結合生成NO3的數量減小.HNO3主要為NO3與H原子結合以及NO2與HO結合而成,即NO3+H→HNO3和NO2+HO→HNO3.由于壓力抑制了NO2和NO3的產生,導致在壓力作用下HNO3的數量相應減少.

3.2.1.2 HONO離解

壓力對HONO離解的影響與N－NO2鍵斷裂類似,HONO的數量隨密度的增加逐漸減少,如圖7(a)所示.當ρ≥2.80 g·cm-3時,HONO的數量低于HNO3;對ρ=3.20 g·cm-3的體系,只有少量的HONO生成,其數量低于O2和HO.HONO主要由硝基(－NO2)上的O原子吸引－CH2上的H原子產生.2在壓力作用下,HMX分子中二面角θ(H2－C1－N2－N1)和θ(H3－C2－N3－N4)隨密度的增加逐漸增大(見Supporting Information中表S2所示),表明壓力使－NO2和－CH2向兩個相反的方向彎曲,增大了N－O…H－C之間的距離,從而降低了硝基與氫原子結合生成HONO的概率.特別是高壓下(ρ≥2.80 g·cm-3)二面角增大明顯,嚴重阻礙了硝基與氫原子的結合.另外,由于高壓抑制了游離態NO2的產生,導致其與H原子結合生成HONO的概率減小,進一步造成HONO數量的減少.這與Strachan等22對RDX在不同沖擊速度下的研究結果不同,他們發現HONO的數量隨著沖擊速度的增加而有所增加.

圖6 T=2500 K時不同壓縮態β-HMX晶體中平均每個HMX分子的分解產物NO2(a)、斷裂的N－NO2鍵(b)、NO3(c)和HNO3(d)的數量隨時間的變化Fig.6 Evolution of quantities of NO2(a),broken N－NO2bonds(b),NO3(c),and HNO3(d)per HMX molecule for β-HMX crystals with different densities at T=2500 K

圖7 T=2500 K時不同壓縮態β-HMX晶體中平均每個HMX分子的分解產物HONO(a)、NO(b)和HNO(c)的數量隨時間的變化Fig.7 Evolution of quantities of HONO(a),NO(b),and HNO(c)per HMX molecule for β-HMX crystals with different densities at T=2500 K

高壓對HONO離解的阻礙作用導致NO和HNO數量相應減少,如圖7(b,c)所示.NO和HNO主要由HONO的分解產生:HONO→HO+NO,NO+ H→HNO.由于壓力抑制了HONO的產生,造成NO減少,進而NO與H結合的概率減小,導致HNO的數量相應減少.對反應路徑的分析發現,HONO的分解是一個可逆反應,即HONO分解產生的HO與NO會再次結合生成HONO:HO+NO→HONO,進一步導致NO和HNO減少.

3.2.1.3 分子主環斷裂

HMX分子主環斷裂的主要產物是H2CO、HCON和HCN.這些產物的數量較少,且隨著體系密度的增加有所減少,特別是高壓縮態體系,但是這些因主環斷裂而生成的產物出現的時間卻隨著密度的增加有所提前,如Supporting Information中圖S1(a-c)所示.這說明壓力會加快HMX分子主環斷裂的速率,但高壓會抑制C－N鍵斷裂后產物的進一步分解.在壓力作用下自由空間減小,狹小的空間不利于環形大分子的存在,HMX分子受壓變形,主環中二面角發生扭轉,促進了C－N鍵斷裂. Goto等28認為由壓力引起的分子變形可以引發化學反應,因而在研究含能材料的點火起爆時不能忽視.HMX分子主環中二面角θ(N2－C2－N3－C1)和θ(N3－C1－N2－C2)隨密度的變化最為明顯(見Supporting Information中表S2),特別是對ρ≥2.80 g· cm-3的高壓縮態體系,說明高壓下分子主環變形嚴重使得N3－C1和N2－C2鍵容易斷裂.HMX分子主環上的C－N鍵長隨著密度的增大逐漸減小(見Supporting Information中表S1),但N3－C1和N2－C2鍵不同,對高壓縮態體系鍵長反而增加,進一步說明N3－C1和N2－C2鍵在高壓下容易斷裂.這與Zhu等17對β-HMX晶體在沖擊波作用下的化學過程的研究結果相同,他們認為主環的斷裂是主要的初始分解路徑之一.

3.2.1.4 大分子團簇的形成

與常壓下相比,壓縮態HMX體系在分解過程中容易形成大分子團簇,且大分子所含原子數隨體系密度的增加而增加,這與Strachan等21對不同密度RDX晶體的熱分解機理的研究結果相同.圖8顯示了在t=20 ps時生成的最大分子團簇中C、N原子數量隨密度的變化,當ρ≥2.80 g·cm-3時,C、N原子數顯著增加.常壓下HMX分解程度高,分解出的分子質量較小;隨著密度的增大,分子間距減小,自由空間減少,HMX分子相互擠壓變形,C－N鍵斷裂,顯著的空間位阻效應使得C－N鍵斷裂后易于發生分子間反應生成分子量很大的分子團簇,而不易于進一步分解成較小的分子,分解程度降低.

圖8 T=2500 K時不同壓縮態β-HMX晶體在t=20 ps時生成的最大分子團簇所含的C、N原子數Fig.8 Numbers of C and N atoms in the maximum molecule formed at t=20 ps for β-HMX crystals with different densities at T=2500 K

3.2.1.5 小分子的生成

壓力促進了一些分子量很小的分子的形成,如O、H、HO、O2.隨著體系密度的增加,O2和HO出現的時間提前,且數量增加,如Supporting Information中圖S2所示.對ρ≥2.80 g·cm-3的高壓縮態體系,O2和HO隨時間的變化與NO2和HONO的對比見圖9 (a-d).O2出現在HMX的初始分解階段,與N－NO2鍵斷裂的時間相同,早于HONO的生成;數量略少于NO2,與HONO相當.HO出現的時間略遲于NO2,與HONO相同;數量略少于NO2,與HONO相當甚至更多(ρ=3.20 g·cm-3).這說明在高壓下(ρ≥2.80 g· cm-3),O2、HO等小分子的形成在HMX的初始分解機理中起著重要的作用,這與常壓和低壓下HMX的分解機理不同.在低壓下(ρ＜2.80 g·cm-3),HMX的初始分解以N－NO2鍵斷裂和HONO離解等分子內反應為主,O2和HO等小分子很少出現.Strachan等22對不同沖擊速度下RDX晶體的熱分解機理的研究也表明由分子間反應生成的HO隨沖擊速度的增加而顯著增加.

對化學反應路徑的分析發現,高壓下O、H主要是從HMX分子中直接脫離而成.在高壓作用下, HMX分子受到擠壓嚴重變形,HMX分子主環以外的部分C－H鍵和N－O鍵伸長,使得H、O原子比較容易脫落.不同壓縮態HMX分子中的C－H鍵長和N－O鍵長隨著密度的增加而減小,但是對高壓縮態體系,C2－H4、N1－O1及N4－O3鍵反而伸長(如Supporting Information中表S1所示),說明這些鍵在高壓下變得比較容易斷裂,從而生成較多游離的H、O原子.Zhu等17采用從頭算分子動力學方法對β-HMX晶體在沖擊波作用下的化學過程的研究也表明,N－O鍵的斷裂在初始分解中起著重要作用.Margetis等10應用SCC-DFTB方法對凝聚態NM在靜水和單軸壓縮下結構的研究發現,壓縮作用使C－H鍵拉伸,最終造成氫質子的解離.

圖9 T=2500 K時高壓縮態β-HMX晶體中平均每個HMX分子的分解產物O2(a)、NO2(b)、HO(c)和HONO(d)的數量對比Fig.9 Comparisons of quantities of O2(a),NO2(b),HO(c),and HONO(d)per HMX molecule for β-HMX crystals with high densities at T=2500 K

這些游離小分子容易在較小的空間存在,且碰撞后進一步生成HO、O2等小分子.另外,高壓下分子間距顯著減小,分子主環外的相鄰原子比較容易結合而生成HO、O2等小分子.對ρ=1.89,2.11,2.22, 2.46,2.80,3.20 g·cm-3的體系,分子間最近的H原子與 O原子間距分別為 0.2674、0.2586、0.2535、0.2423、0.2296和0.2073 nm,即分子間N－O…HC間距隨著密度的增加而減小,有利于分子間H、O原子的結合.對反應路徑的分析表明,HO主要來源于HONO的分解(HONO→HO+NO)以及H、O原子的直接結合(O+H→HO).高壓抑制了HONO的產生,因此通過HONO分解產生的HO很少,大量的HO主要為后者產生.

3.2.2 最終產物的形成

HMX分解的最終產物主要是H2O、N2及CO2,如圖5所示.隨著體系密度的增加,H2O出現的時間逐漸推后,數量逐漸減少.對ρ=1.89,2.11,2.22,2.46, 2.80,3.20 g·cm-3的體系,H2O出現的時間分別為0.1、0.15、0.25、0.4、1.0和2.15 ps,最大值分別為2.0、1.85、1.9、1.6、1.25、0.75,如Supporting Information中圖S3(a)所示.對不同壓縮態體系生成的H2O隨時間的變化按式(4)進行指數擬合,得到H2O的反應速率k3分別為0.061、0.056、0.050、0.045、0.032和0.016 ps-1,說明H2O的反應速率隨著密度的增加逐漸減小,特別是高壓下反應速率顯著降低,說明壓力不利于H2O的產生.對反應路徑的分析發現H2O的產生與NO2的生成和HONO分解為HO和NO有關,即(1)C4H8N8O8→C4H8N7O6+NO2,(2)NO2+H→HONO,(3)HONO→HO+NO,(4)HO+H→H2O.NO2和HONO的數量隨密度的增加逐漸減少,進而分解產生的HO和NO也相應減少,最終導致H2O的減少.且通過對反應路徑的分析還發現,高壓下HO+ H→H2O是一個可逆反應,即壓力促進了H2O分解為HO和H,進一步造成H2O的減少.

N2的數量隨密度的變化比較復雜,并不是隨密度的增加而單調變化.當ρ＜2.46 g·cm-3,N2的數量隨著密度的增加逐漸增大;當ρ≥2.46 g·cm-3,N2的數量隨密度的繼續增大而減小,如Supporting Information中圖S3(b)所示.對ρ=1.89,2.11,2.22,2.46,2.80, 3.20 g·cm-3的體系,N2的數量在t=20 ps時分別達到0.85、1.25、1.3、1.1、0.75和0.5,說明低壓會促進N2的生成,而高壓會抑制N2的產生.按式(4)擬合得到的不同壓縮態體系生成N2的反應速率k3分別為0.027、0.040、0.042、0.035、0.021和0.018 ps-1,說明低壓下反應速率隨密度的增加而升高,高壓下反應速率隨密度的繼續增加而降低.對反應路徑的分析發現N2主要源于HMX分子主環的斷裂生成C2N2O2及CHN2O2,進一步分解成CN2O和HN2O,最終生成N2.由于壓力促進了主環的變形,加速了C－N鍵的斷裂,而高壓抑制了C－N鍵斷裂后產物的進一步分解,最終導致低壓下N2的數量隨密度的增加而增加,而高壓下隨密度的增加而減少.N－NO2鍵斷裂后產物的進一步反應也會生成N2.由于高壓顯著地抑制了N－N鍵的斷裂,因而通過這一路徑生成的N2也會減少.在不同沖擊速度下,RDX熱分解生成的N2隨著沖擊速度的增加而顯著增加,22這與本文的研究結果不完全相同.

壓力對CO2的影響與N2相似,但沒有N2顯著.當ρ＜2.46 g·cm-3,CO2的數量隨著密度的增加稍微有所增加;當ρ≥2.46 g·cm-3,CO2的數量隨著密度的繼續增加而有所減少,如Supporting Information中圖S3(c)所示.對ρ=3.20 g·cm-3的體系,CO2的數量反而有所增加,這主要是由于高壓下生成了較多的O和O2,使得C容易被氧化而生成較多的CO2.CO2是HMX主環斷裂的最終產物,其隨密度的變化說明壓力對主環斷裂的影響較為復雜,壓力會加速CN鍵斷裂的速率但高壓會抑制C－N鍵斷裂后產物的進一步分解,導致最終產物的減少.

4 結論

采用ReaxFF反應力場和分子動力學方法,研究了不同壓縮態的HMX晶體在T=2500 K時的熱分解過程.研究結果表明,HMX的初始分解機理與密度或壓力有關,高壓下(ρ≥2.80 g·cm-3)的初始熱分解機理與低壓(ρ＜2.80 g·cm-3)下有所不同.常壓下, HMX晶體熱分解過程中三種主要的初始分解路徑是N－NO2鍵的斷裂、HONO的離解和C－N主環的斷裂.低壓下,HMX的初始分解機理以分子內NNO2鍵斷裂和HONO離解為主,初始分解速率隨密度增加而降低,體系以吸熱反應開始,勢能隨時間先升高再降低.高壓下HMX晶體的熱分解程度降低,但分子間反應得到促進而生成了較多的HO、O2等小分子和大分子團簇,使得HMX單分子存在的時間變短,初始分解速率隨密度增加而升高,體系主要以放熱反應開始,勢能隨時間逐漸降低.高壓縮態下HMX的初始熱分解機理與低壓縮態下最主要的區別是前者除部分分子內反應以外還有較多的分子間反應,而后者以分子內反應為主.次級反應速率隨密度的增加有所增加,說明壓力會加快HMX的次級反應.本文的研究工作有助于深入認識炸藥在高溫高壓下的初始反應機理、反應速率、詳細的反應路徑,以及壓力對化學響應的影響規律,這對于炸藥的沖擊起爆及極端條件下的化學反應研究具有借鑒意義和參考價值.

壓力對HMX初始熱分解機理的影響是否受溫度的影響還有待研究.與RDX在不同沖擊速度下的分解產物相比,靜水壓縮與沖擊作用對含能材料的影響不完全相同;與HMX在一定沖擊作用下的分解機理相比，高壓縮態體系的分解機理與其有一些相似之處.因而有必要進一步研究凝聚態HMX在不同沖擊作用下的化學反應,以得到靜水壓縮與沖擊壓縮對材料化學響應的影響及區別.

ReaxFF力場的發展為研究大體系的物理化學性質提供了非常有效的途徑,且是目前主流的能用于研究凝聚態材料化學反應的力場.在含能材料領域,ReaxFF反應分子動力學方法得到了廣泛而成功的應用,為凝聚態含能材料的化學反應提供了非常有價值的信息.當然也應該看到ReaxFF力場存在的不足.由于ReaxFF勢函數比較復雜,力場參數較多,在通過擬合QM或實驗數據得到力場參數時會產生一定的誤差,且力場的參數化是一項具有較大難度和挑戰性的工作,這使得ReaxFF的計算精度不如QM.除了力場參數,力場的函數形式也至關重要,比如對不在一個數量級的鍵能與分子間作用能的處理,不完善的函數形式可能會使計算結果不夠準確.因此,進一步優化力場參數,發展函數形式是必要且有意義的工作.

Supporting Information Available: Evolution of quantities of H2CO,HCON,and HCN per HMX molecule,evolution of quantities of O2and HO per HMX molecule,and evolution of quantities of H2O,N2,and CO2per HMX molecule for β-HMX crystals with different densities at T=2500 K shown in Fig.S1, Fig.S2,and Fig.S3,respectively;bond length and dihedral an-gle in β-HMX molecule for crystals with different densities shown in Table S1 and Table S2,respectively,have been included.This information is available free of charge via the internet at http://www.whxb.pku.edu.cn.

(1) Lewis,J.P.;Glaesemann,K.R.;VanOpdorp,K.;Voth,G.A. J.Phys.Chem.A 2000,104,11384.doi:10.1021/jp002173g

(2) Chakraborty,D.;Muller,R.P.;Goddard,W.A.,III.J.Phys. Chem.A 2001,105,1302.doi:10.1021/jp0026181

(3) Sharia,O.;Kuklja,M.M.J.Phys.Chem.B 2011,115,12677.

(4) Jiang,F.L.;Zhai,G.H.;Ding,L.;Yue,K.F.;Liu,N.;Shi,Q. Z.;Wen,Z.Y.Acta Phys.-Chim.Sin.2010,26,409.[姜富靈,翟高紅,丁 黎,岳可芬,劉 妮,史啟禎,文振翼.物理化學學報,2010,26,409.]doi:10.3866/PKU.WHXB20100128

(5) Brill,T.B.J.Prop.Power 1995,11,740.doi:10.2514/3.23899

(6) Tang,C.J.;Lee,Y.J.;Litzinger,T.A.J.Prop.Power 1999,15, 296.doi:10.2514/2.5427

(7) Tarver,C.M.;Chidester,S.K.;Nichols,A.L.J.Phys.Chem. 1996,100,5794.doi:10.1021/jp953123s

(8) Gilman,J.J.Phil.Maga.B 1995,71(6),1057.doi:10.1080/ 01418639508241895

(9) Gilman,J.J.Phil.Maga.B 1993,67(2),207.doi:10.1080/ 13642819308207868

(10) Margetis,D.;Kaxiras,E.;Elstner,M.;Frauenheim,T.;Manaa, M.R.J.Chem.Phys.2002,117(2),788.doi:10.1063/ 1.1466830

(11) Manaa,M.R.Appl.Phys.Lett.2003,83(7),1352.doi:10.1063/ 1.1603351

(12) Lu,L.Y.;Wei,D.Q.;Chen,X.R.;Lian,D.;Ji,G.F.;Zhang,Q. M.;Gong,Z.Z.Mol.Phys.2008,106,2569.doi:10.1080/ 00268970802616343

(13) Kuklja,M.M.;Rashkeev,S.N.;Zerilli,F.J.Appl.Phys.Lett. 2006,89(7),71904.doi:10.1063/1.2335680

(14) Kuklja,M.M.;Rashkeev,S.N.Phys.Rev.B 2007,75(10), 104111.doi:10.1103/PhysRevB.75.104111

(15) Manaa,M.R.;Fried,L.E.;Melius,C.F.;Elstner,M.; Frauenheim,T.J.Phys.Chem.A 2002,106(39),9024.doi: 10.1021/jp025668+

(16)Manaa,M.R.;Fried,L.E.;Reed,E.J.J.Computer-Aided Materials Design 2003,10(2),75.doi:10.1023/B:JCAD. 0000036812.64349.15

(17)Zhu,W.H.;Huang,H.;Huang,H.J.;Xiao,H.M.J.Chem. Phys.2012,136,044516.doi:10.1063/1.3679384

(18) van Duin,A.C.T.;Dasgupta,S.;Lorant,F.J.Phys.Chem.A 2001,105(41),9396.

(19) Rom,N.;Zybin,S.V.;van Duin,A.C.T.;Goddard,W.A.,III; Zeiri,Y.;Katz,G.;Kosloff,R.J.Phys.Chem.A 2011,115, 10181.doi:10.1021/jp202059v

(20) Zhang,L.Z.;Sergey,V.Z.;van Duin,A.C.T.;Siddharth,D.; Goddard,W.A.,III.J.Phys.Chem.A 2009,113,10619.

(21) Strachan,A.;Kober,E.M.;van Duin,A.C.T.;Oxgaard,J.; Goddard,W.A.,III.J.Chem.Phys.2005,122(5),54502.doi: 10.1063/1.1831277

(22) Strachan,A.;van Duin,A.C.T.;Chakraborty,D.;Dasgupta,S.; Goddard,W.A.,III.Phys.Rev.Lett.2003,91(9),098301.doi: 10.1103/PhysRevLett.91.098301

(23) Zybin,S.V.;Goddard,W.A.,III;Xu,P.;van Duin,A.C.T.; Appl.Phys.Lett.2010,96,081918.doi:10.1063/1.3323103

(24)An,Q.;Liu,Y.;Zybin,S.V.;Kim,H.;Goddard,W.A.,III. J.Phys.Chem.C 2012,116(18),10198.doi:10.1021/ jp300711m

(25)Zhou,T.T.;Zybin,S.V.;Liu,Y.;Huang,F.L.;Goddard,W.A., III.J.Appl.Phys.2012,111(12),124904.doi:10.1063/ 1.4729114

(26) Choi,C.S.;Boutin,H.P.Acta Crystallogr.B 1970,26,1235.

(27) Yoo,C.S.;Cynn,H.J.Chem.Phys.1999,111(22),10229.doi: 10.1063/1.480341

(28)Goto,N.;Fujihisa,H.;Yamawaki,H.J.Phys.Chem.B 2006, 110,23655.doi:10.1021/jp0635359

June 4,2012;Revised:August 3,2012;Published on Web:August 3,2012.

Thermal Decomposition Mechanism of β-HMX under High Pressures via ReaxFF Reactive Molecular Dynamics Simulations

ZHOU Ting-Ting SHI Yi-Ding HUANG Feng-Lei*
(State Key Laboratory of Explosion Science and Technology,Beijing Institute of Technology,Beijing 100081,P.R.China)

The thermal decomposition mechanisms of condensed phase β-HMX at various densities(ρ= 1.89,2.11,2.22,2.46,2.80,3.20 g·cm-3)and at 2500 K were studied using ReaxFF reactive molecular dynamics simulations.The effects of pressure on the initial and secondary reaction rates,the main differences in the initial decomposition mechanisms between highly compressed and less compressed systems,as well as the reasons for these variations were analyzed.It was determined that the initial decomposition mechanisms of HMX were dependent on pressure(or density).At low densities(ρ＜2.80 g· cm-3),intramolecular reactions are dominant,these being N－NO2bond dissociation,HONO elimination, and concerted ring fission by C－N bond scission.At high densities(ρ≥2.80 g·cm-3),intramolecular reactions are well restrained,whereas intermolecular reactions are promoted,leading to the formation of small molecules,such as O2and HO,and large molecular clusters.These changes in the initial decomposition mechanisms lead to different kinetic and energetic behaviors,as well as variations in the distribution of products.These results obtained through this work are significant in that they assist in understanding the chemical reactions involved in the initiation,reaction development,and detonation of energetic materials under extreme conditions.

HMX;Thermal decomposition;Pressure;ReaxFF;Molecular dynamics

10.3866/PKU.WHXB201208031

?Corresponding author.Email:huangfl@bit.edu.cn;Tel:+86-10-68914518.

The project was supported by the National Natural Science Foundation of China(10832003).

國家自然科學基金(10832003)資助項目

O643;O642