高 超 胡麗娟 張巧鳳 劉 珂 梁婉怡 謝耀平

(上海大學(xué)材料科學(xué)與工程學(xué)院材料研究所、微結(jié)構(gòu)重點(diǎn)實(shí)驗(yàn)室，上海 200072)

亞穩(wěn)態(tài)分子間復(fù)合物(metastable intermolecular composites, MIC)是一種由納米尺度的燃料(通常為Al)和氧化劑(MoO3、CuO、Fe2O3等)顆粒組成的新型含能材料。未加載環(huán)境下，MIC處于亞穩(wěn)態(tài)，在外界載荷(如激光、熱誘導(dǎo)等)下會(huì)發(fā)生反應(yīng)形成一種更穩(wěn)定的合金或復(fù)合物。MIC具有燃燒效率高、能量釋放快等特點(diǎn)。如Al和MoO3構(gòu)成的MIC含能材料發(fā)生反應(yīng)時(shí)，反應(yīng)傳播速率在600～1 000 m/s之間，能量釋放為1.12 kcal/g[1]，這可以與常見(jiàn)的硝胺基爆炸物相媲美。

為了探討MIC燃燒的反應(yīng)機(jī)制，許多學(xué)者對(duì)其做了研究。Asay等[2]試驗(yàn)研究了MIC材料Al/MoO3體系的燃燒過(guò)程，認(rèn)為對(duì)流機(jī)制是MIC材料燃燒的反應(yīng)機(jī)制。Dean等[3]以Al/CuO和Al/NiO作為試驗(yàn)材料，推測(cè)MIC的燃燒機(jī)制為熔化- 爆裂機(jī)制。但由于MIC燃燒反應(yīng)一般都在遠(yuǎn)離平衡態(tài)或者是極端條件下發(fā)生，這給試驗(yàn)研究帶來(lái)了一定困難，從而難以確定MIC材料燃燒的反應(yīng)機(jī)制。分子動(dòng)力學(xué)(molecular dynamics，MD)作為一種理論研究手段可以彌補(bǔ)試驗(yàn)方法在時(shí)間和空間尺度上的局限性，能在原子尺度上直觀地觀測(cè)到反應(yīng)過(guò)程中各個(gè)原子的運(yùn)動(dòng)狀況，進(jìn)而研究MIC燃燒的反應(yīng)機(jī)制。

MD采用經(jīng)驗(yàn)的勢(shì)函數(shù)描述原子的運(yùn)動(dòng)，因此找到合理的勢(shì)函數(shù)來(lái)準(zhǔn)確描述MIC的反應(yīng)過(guò)程是必要的。對(duì)勢(shì)可以較好地描述在平衡狀態(tài)的體系，嵌入原子勢(shì)[4](embedded- atom method，EAM)可以對(duì)金屬體系進(jìn)行較好的描述。但它們都不能描述鍵的斷裂和生成，不適用于描述MIC的燃燒反應(yīng)。由Van Duin等[5]開發(fā)的勢(shì)函數(shù)——反應(yīng)力場(chǎng)(reactive force filed，ReaxFF)可以精確描述復(fù)雜化學(xué)反應(yīng)過(guò)程中原子間的相互作用力，能很好地模擬鍵的斷裂與形成過(guò)程，是研究MIC燃燒反應(yīng)的有效方法。目前，運(yùn)用ReaxFF已經(jīng)成功描述了有機(jī)物[5]、金屬氧化態(tài)[6]、Al/Al2O3界面[7]等結(jié)構(gòu)的物理和化學(xué)性質(zhì)。由于必須獲取材料的ReaxFF參數(shù)才可以運(yùn)用ReaxFF對(duì)其進(jìn)行MD模擬，為了獲取適用于MIC材料的ReaxFF參數(shù)，需要對(duì)ReaxFF參數(shù)進(jìn)行擬合。通過(guò)將擬合獲得的ReaxFF參數(shù)引入ReaxFF，才可以將ReaxFF運(yùn)用于MIC體系的計(jì)算模擬。

本文選取典型的MIC材料Al/MoO3體系，擬合其ReaxFF參數(shù)，并驗(yàn)證ReaxFF對(duì)該體系的適用情況。通過(guò)分析ReaxFF勢(shì)函數(shù)和力場(chǎng)擬合方法，綜合運(yùn)用擬合好的ReaxFF勢(shì)函數(shù)計(jì)算Al/MoO3體系的一系列性質(zhì)，并將其與第一性原理計(jì)算值或試驗(yàn)值進(jìn)行對(duì)比，驗(yàn)證ReaxFF對(duì)Al/MoO3體系的適用情況。

1 ReaxFF反應(yīng)力場(chǎng)模型的建立

1.1 ReaxFF反應(yīng)力場(chǎng)

鍵級(jí)是ReaxFF的核心概念，它能很好地反映化學(xué)鍵的信息，從而能夠描述鍵的斷裂和形成過(guò)程。如式(1)～式(3)即為ReaxFF計(jì)算系統(tǒng)能量的表達(dá)式[5]。

Esystem=Ebond+Eangle+Etorsion+

Evdwaals+ECoulomb

(1)

(2)

(3)

1.2 ReaxFF擬合方法

ReaxFF參數(shù)決定著MIC的反應(yīng)過(guò)程分子動(dòng)力學(xué)模擬的準(zhǔn)確度及精度，因此擬合出合理可靠的ReaxFF參數(shù)可以說(shuō)是MIC的反應(yīng)過(guò)程分子動(dòng)力學(xué)模擬的一項(xiàng)核心，并具有一定挑戰(zhàn)性。ReaxFF參數(shù)的擬合主要分為三部分。首先，獲得合適的訓(xùn)練基(training sets)。訓(xùn)練基是應(yīng)用于ReaxFF參數(shù)擬合過(guò)程中的一組輸入數(shù)據(jù)，為第一性原理計(jì)算值或試驗(yàn)值，組成訓(xùn)練基的數(shù)據(jù)有鍵能、鍵角能、生成焓、晶格常數(shù)、狀態(tài)方程、空位形成能等。下一步是反應(yīng)力場(chǎng)的擬合，由于ReaxFF有大量的參數(shù)需要擬合，本文采用單參數(shù)連續(xù)優(yōu)化方法(successive one- parameter search method)[9]依次優(yōu)化一套力場(chǎng)參數(shù)中的每個(gè)參數(shù)，每次優(yōu)化使得ReaxFF計(jì)算值與訓(xùn)練基中標(biāo)準(zhǔn)值的偏差平方和最小，此時(shí)得到的參數(shù)即為該參數(shù)的最佳值。接著通過(guò)該自洽方法依次循環(huán)優(yōu)化每個(gè)參數(shù)，最終得到一套最佳的力場(chǎng)參數(shù)。最后，對(duì)該套力場(chǎng)參數(shù)進(jìn)行驗(yàn)證，確定其合理性。

2 結(jié)果與討論

為了擬合出適用于MIC材料Al/MoO3體系的ReaxFF參數(shù)，以前期研究開發(fā)的Al/Mo/O ReaxFF參數(shù)[10]為基礎(chǔ)，開發(fā)適用于Al/MoO3體系的Mo- O二元與Al- Mo- O三元參數(shù)，從而完善該ReaxFF參數(shù)，使其可以描述MIC材料Al/MoO3體系。

以Mo- O二元結(jié)構(gòu)MoO2、MoO3、Mo4O11、Mo8O23、Mo9O26和Al- Mo- O三元結(jié)構(gòu)Al8Mo12O48、Al16Mo24O96為研究對(duì)象，運(yùn)用第一性原理(density functional theory，DFT)計(jì)算這些結(jié)構(gòu)的幾何參數(shù)和能量作為訓(xùn)練基來(lái)擬合ReaxFF參數(shù)。DFT計(jì)算采用了VASP[11]程序包，利用投影綴加平面波方法(PAW)來(lái)描述電子與原子間作用，交換關(guān)聯(lián)泛函采用GGA- PW91處理交換關(guān)聯(lián)能，布里淵區(qū)內(nèi)則采用MP方法自動(dòng)產(chǎn)生K點(diǎn)。分子動(dòng)力學(xué)模擬則采用大型分子動(dòng)力軟件LAMMPS[12]進(jìn)行。

2.1 能量變化

由于DFT方法可以較精確地描述結(jié)構(gòu)的能量，可以視其狀態(tài)方程(equation of states，EOS)曲線為真實(shí)值。因此采用擬合完畢的Al/Mo/O ReaxFF分別計(jì)算了Mo- O二元結(jié)構(gòu)和Al- Mo- O三元結(jié)構(gòu)能量隨晶格常數(shù)變化的EOS，將其與DFT計(jì)算結(jié)果作對(duì)比，以驗(yàn)證ReaxFF對(duì)Al/MoO3體系的適用情況，如圖1所示。由圖1可以看出，ReaxFF勢(shì)對(duì)Al/MoO3體系的各種結(jié)構(gòu)的描述較好，能夠準(zhǔn)確區(qū)分出Al/MoO3體系的不同結(jié)構(gòu)，并較準(zhǔn)確地描述各個(gè)結(jié)構(gòu)的能量變化。需要指出的是，圖1(d)中，Al8Mo12O48和Al16Mo24O96的EOS曲線與DFT計(jì)算的結(jié)果有些許偏差，這是由于ReaxFF可以精確描述復(fù)雜化學(xué)反應(yīng)過(guò)程中原子間的相互作用力，但對(duì)金屬鍵的描述能力不如DFT那樣準(zhǔn)確。比如宋文雄等[13]使用ReaxFF對(duì)金屬Al以及Narayanan等[14]使用ReaxFF對(duì)Li/Al/Si/O體系含有金屬鍵結(jié)構(gòu)的描述都出現(xiàn)了一定的偏差。但這種情況并不影響力場(chǎng)的應(yīng)用，因?yàn)榱?chǎng)可以相對(duì)準(zhǔn)確描述出Al8Mo12O48和Al16- Mo24O96結(jié)構(gòu)的整體能量變化。

圖1 運(yùn)用DFT和ReaxFF勢(shì)函數(shù)計(jì)算的Al/MoO3體系結(jié)構(gòu)的狀態(tài)方程(EOS)曲線Fig.1 Equation of states of Al/MoO3 structures calculated by using DFT and ReaxFF

為了驗(yàn)證擬合好的ReaxFF對(duì)Al/MoO3體系表面的描述能力，以Al/MoO3體系晶體結(jié)構(gòu)中原子排布清晰的表面為研究對(duì)象，分別用DFT和ReaxFF計(jì)算了這些表面的表面能，如表1所示。由表1可以看出，ReaxFF計(jì)算結(jié)果與DFT計(jì)算結(jié)果一致，說(shuō)明擬合后的ReaxFF能夠描述Al/MoO3體系結(jié)構(gòu)的表面能。

2.2 晶體結(jié)構(gòu)

表1 運(yùn)用DFT與ReaxFF方法計(jì)算得到的材料表面能Table 1 Surface energy of some structures calculated by using DFT and ReaxFF eV

運(yùn)用DFT和ReaxFF計(jì)算了Mo- O二元結(jié)構(gòu)和Al- Mo- O三元結(jié)構(gòu)的晶格常數(shù)與結(jié)合能，如表2所示。運(yùn)用ReaxFF計(jì)算的不同結(jié)構(gòu)的晶格常數(shù)和結(jié)合能值與DFT計(jì)算值或試驗(yàn)值相差均小于5%，說(shuō)明擬合得到的力場(chǎng)參數(shù)較好地重現(xiàn)了Al/MoO3體系各個(gè)結(jié)構(gòu)的晶格參數(shù)和結(jié)合能，可以較準(zhǔn)確地描述MIC材料Al/MoO3體系的晶體結(jié)構(gòu)。ReaxFF計(jì)算得到的Al8Mo12O48和Al16- Mo24O96結(jié)構(gòu)的結(jié)合能分別為3.61和3.60 eV，這與DFT計(jì)算得到的2.89和2.87 eV相比有一定偏差，這是由于ReaxFF對(duì)金屬鍵的描述能力較弱造成的。另外，可以發(fā)現(xiàn)，DFT計(jì)算的晶格常數(shù)比試驗(yàn)值要大一些，這是由于DFT計(jì)算采用的是GGA勢(shì)函數(shù)，GGA勢(shì)函數(shù)在描述氧化物時(shí)會(huì)引起晶格常數(shù)偏大。Scanlon等[15]使用GGA勢(shì)函數(shù)對(duì)MoO2和MoO3結(jié)構(gòu)的晶格常數(shù)進(jìn)行了DFT計(jì)算，發(fā)現(xiàn)這些結(jié)構(gòu)的晶格常數(shù)均大于試驗(yàn)值，與本文結(jié)果一致。

表2 運(yùn)用DFT與ReaxFF方法計(jì)算得到的晶格常數(shù)和結(jié)合能及其試驗(yàn)值Table 2 Lattice constants and binding energies calculated by DFT and ReaxFF, and their experimental data

2.3 彈性性質(zhì)

晶體的狀態(tài)方程、熱容、德拜溫度、熔點(diǎn)等固態(tài)性質(zhì)都與其彈性性質(zhì)密切相關(guān)。Yun等[17]開發(fā)了Fe/Al/Ni合金的ReaxFF勢(shì)函數(shù)，并使用其計(jì)算了Fe/Al/Ni合金的彈性常數(shù)以驗(yàn)證ReaxFF適用性。Zhang等[7]運(yùn)用了ReaxFF計(jì)算了α- Al2O3的彈性常數(shù)，發(fā)現(xiàn)ReaxFF能夠準(zhǔn)確描述α- Al2O3的彈性性能?？梢?jiàn)，彈性性能已經(jīng)成為驗(yàn)證ReaxFF適用性的有效手段。通過(guò)計(jì)算彈性常數(shù)，可以驗(yàn)證擬合好的ReaxFF參數(shù)對(duì)晶體彈性性能的描述情況。根據(jù)彈性應(yīng)變能與應(yīng)變的關(guān)系來(lái)計(jì)算彈性常數(shù)，如式(4)所示。

(4)

對(duì)晶體施加不同的應(yīng)變ei、ej，根據(jù)晶體施加應(yīng)變前后應(yīng)變能的變化△E及晶體體積V即可求得晶體的彈性常數(shù)Cij。分別使用DFT與ReaxFF計(jì)算正交晶系MoO3、Mo4O11和Al8Mo12O48結(jié)構(gòu)的彈性常數(shù)，以驗(yàn)證擬合后的力場(chǎng)參數(shù)對(duì)晶體力學(xué)穩(wěn)定性的描述情況。正交晶系的彈性常數(shù)為9個(gè)，分別為C11、C12、C13、C23、C22、C33、C44、C55、C66。表3給出了用DFT和ReaxFF方法計(jì)算得到各個(gè)結(jié)構(gòu)的彈性常數(shù)。由表3可以看出，對(duì)于不同結(jié)構(gòu)，ReaxFF計(jì)算的彈性常數(shù)與DFT計(jì)算結(jié)果有一定偏差。由于ReaxFF添加了Morse函數(shù)來(lái)描述范德瓦爾斯力，而Morse函數(shù)考慮兩個(gè)原子間的相互作用，但不會(huì)考慮鍵角能、二面角能等多體相互作用。另一方面，ReaxFF又通過(guò)計(jì)算鍵級(jí)來(lái)描述多體相互作用。在擬合ReaxFF參數(shù)時(shí)，可能并未處理好Morse函數(shù)表示的范德華作用與鍵級(jí)描述的多體相互作用，導(dǎo)致ReaxFF計(jì)算的彈性常數(shù)偏大。但是ReaxFF計(jì)算的彈性常數(shù)與DFT計(jì)算結(jié)果的數(shù)量級(jí)一致，而且不同彈性常數(shù)間大小關(guān)系也幾乎與DFT計(jì)算結(jié)果一致。式(5)～式(7)給出了正交晶系的力學(xué)穩(wěn)定性判據(jù)，采用DFT和ReaxFF計(jì)算得到的彈性常數(shù)都滿足此判據(jù)，表明這兩種方法得到的晶體結(jié)構(gòu)是穩(wěn)定的，這與事實(shí)相符。因此ReaxFF預(yù)測(cè)Al/MoO3體系的彈性常數(shù)值雖然大于DFT計(jì)算結(jié)果，但是ReaxFF可以準(zhǔn)確地描述Al/MoO3體系的力學(xué)穩(wěn)定性。

C11>0;C22>0;C33>0;

C44>0;C55>0;C66>0

(5)

C11+C22-2C12>0;C11+C33-2C13>0;

C22+C33-2C23>0

(6)

C11+C22+C33+2(C12+C13+C23)>0

(7)

表3 運(yùn)用DFT與ReaxFF方法計(jì)算得到的彈性常數(shù)Table 3 Elastic constants of some structures calculated by DFT and ReaxFF kbar

3 結(jié)論

為了使用ReaxFF反應(yīng)力場(chǎng)對(duì)MIC體系燃燒反應(yīng)機(jī)制進(jìn)行理論研究，以常見(jiàn)MIC材料Al/MoO3體系為研究對(duì)象，對(duì)其ReaxFF參數(shù)的Mo- O二元與Al- Mo- O三元參數(shù)進(jìn)行了擬合優(yōu)化，并得到了一套較為完善的力場(chǎng)參數(shù)。運(yùn)用該力場(chǎng)參數(shù)計(jì)算了Mo- O和Al- Mo- O結(jié)構(gòu)的物理性質(zhì)，如狀態(tài)方程、晶格常數(shù)、表面能和彈性常數(shù)。并將計(jì)算結(jié)果與第一性原理(DFT)計(jì)算結(jié)果進(jìn)行對(duì)比，以驗(yàn)證該力場(chǎng)對(duì)Al/MoO3體系的適用性。

擬合好的Al/Mo/O ReaxFF勢(shì)對(duì)Mo- O和Al- Mo- O結(jié)構(gòu)的狀態(tài)方程(EOS)和晶格常數(shù)描述較好，能夠較準(zhǔn)確地描述各個(gè)結(jié)構(gòu)的能量變化、晶體結(jié)構(gòu)和表面能。在描述晶體彈性性能方面，擬合好的ReaxFF計(jì)算的Al/MoO3體系的彈性常數(shù)值，雖與DFT計(jì)算結(jié)果有一定偏差，但可以準(zhǔn)確描述Al/MoO3體系的力學(xué)穩(wěn)定性。這說(shuō)明擬合好的ReaxFF可以較準(zhǔn)確地描述Al/MoO3體系的物理性質(zhì)，可作為研究MIC材料Al/MoO3燃燒反應(yīng)機(jī)制的分子動(dòng)力學(xué)勢(shì)函數(shù)。

[1] BOCKMON B S, PANTOYA M L, SON S F, et al. Combustion velocities and propagation mechanisms of metastable interstitial composites [J]. Journal of Applied Physics, 2005,98(6): 064903.

[2] ASAY B W, SON S F, BUSSE J R, et al. Ignition characteristics of metastable intermolecular composites [J]. Propellants Explosives Pyrotechnics, 2004, 29(4): 216- 219.

[3] DEAN S W, PANTOYA M L, GASH A E, et al. Enhanced convective heat transfer in nongas generating nanoparticle thermites [J]. Journal of Heat Transfer, 2010, 132(11): 2201- 2210.

[4] DAW M S, BASKES M I. Embedded- atom method: Derivation and application to impurities, surfaces, and other defects in metals [J]. Physical Review B, 1984, 29(12): 6443- 6453.

[5] VAN DUIN A C T, DASGUPTA S, LORANT F, et al. ReaxFF: a reactive force field for hydrocarbons [J]. The Journal of Physical Chemistry A, 2001, 105(41): 9396- 9409.

[6] ARYANPOUR M, VAN DUIN A C T, KUBICKI J D. Development of a reactive force field for iron- oxyhydroxide systems[J]. The Journal of Physical Chemistry A,2010,114(21): 6298- 6307.

[7] ZHANG Q, CAGIN T, DUIN A V, et al. Adhesion and nonwetting- wetting transition in the Al/α- Al2O3interface [J]. Physical Review B, 2004, 69(4): 1129- 1133.

[8] MORTIER W J, GHOSH S K, SHANKAR S. Electronegativity- equalization method for the calculation of atomic charges in molecules [J]. Journal of the American Chemical Society, 1986, 108(15): 4315- 4320.

[9] VAN DUIN A C T, JAN M, DE GRAAF B. Delft molecular mechanics: a new approach to hydrocarbon force fields. Inclusion of a geometry- dependent charge calculation [J]. Journal of the Chemical Society, Faraday Transactions, 1994, 90(19): 2881- 2895.

[10]SONG W X, ZHAO S J. Development of the ReaxFF reactive force field for aluminum- molybdenum alloy [J]. Journal of Materials Research, 2013, 28(9): 1155- 1164.

[11] KRESSE G, FURTHMULLER J. Efficiency of ab- initio total energy calculations for metals and semiconductors using a plane- wave basis set [J]. Computational Materials Science, 1996, 6(1): 15- 50.

[12] PLIMPTON S. Fast parallel algorithms for short- range molecular dynamics [J]. Journal of Computational Physics, 1995, 117(1): 1- 19.

[13] 宋文雄,趙世金.ReaxFF反應(yīng)力場(chǎng)在金屬鋁中的適用性[J].含能材料,2012,20(5): 571- 574.

[14] NARAYANAN B, VAN DUIN A C T, KAPPES B B, et al. A reactive force field for lithium- aluminum silicates with applications to eucryptite phases [J]. Modelling and Simulation in Materials Science and Engineering, 2011, 20(1):386- 401.

[15] SCANLON D O, WATSON G W, PAYNE D, et al. Theoretical and experimental study of the electronic structures of MoO3and MoO2[J]. The Journal of Physical Chemistry C, 2010, 114(10): 4636- 4645.

[16] BREWER L, LAMOREAUX R H,FERRO R. The Mo- O System (Molybdenum- Oxygen) [J]. Journal of Phase Equilibria, 1980, 1(2): 85- 89.

[17] YUN K S, KWAK H, ZOU C, et al. Development and Validation of a ReaxFF Reactive Force Field for Fe/Al/Ni Alloys: Molecular Dynamics Study of Elastic Constants, Diffusion, and Segregation [J]. Journal of Physical Chemistry A, 2012, 116(49): 12163- 12174.