李楊井敏武吉偉劉立強

(1.山東建筑大學(xué) 材料科學(xué)與工程學(xué)院，山東 濟南250101；2.國家包裝產(chǎn)品質(zhì)量監(jiān)督檢驗中心，山東 濟南250100)

0 引言

微晶玻璃作為一種新型無機材料，具備優(yōu)良的化學(xué)、力學(xué)、機械、光電性能，有著廣闊的應(yīng)用前景，相關(guān)研究一直是關(guān)注熱點[1]。 傳統(tǒng)微晶玻璃的制備通常從改變原料配比和熱處理方式等方面開展。為得到一組可行有效的配方，通常需要反復(fù)的實驗，不僅消耗了大量的人力、物力，而且成效低并存在產(chǎn)品易開裂、成品率低、結(jié)構(gòu)和性能關(guān)系尚不明細等問題[2]。 隨著微晶玻璃研究的深入，要求從微觀和介觀角度能夠準確、定量地描述微晶玻璃結(jié)構(gòu)以及析晶過程規(guī)律，如研究微晶玻璃高溫熔體液體結(jié)構(gòu)和性能的關(guān)系、堿金屬氧化物添加量對高溫熔體的黏度的影響、晶核劑的形核機理等。 分子動力學(xué)模擬作為凝聚態(tài)物理學(xué)常用的計算機模擬技術(shù)，從原子、分子層次認識物質(zhì)的組成，利用計算機直觀展示和量化晶體結(jié)構(gòu)，通過構(gòu)建物質(zhì)的微觀結(jié)構(gòu)和數(shù)值模擬熱力學(xué)運動過程，得到每個粒子運動規(guī)律、能量波動等信息[3-4]。 微晶玻璃體系一般性模擬集中在三元或四元體系，完全可以借助基于分子動力學(xué)開發(fā)的模擬軟件開展工作，通過調(diào)整密度大小、粒子的數(shù)目、力場參數(shù)、施加的溫度、壓力，構(gòu)建需要的微晶玻璃模型，從微觀層次上對微晶玻璃進行結(jié)構(gòu)優(yōu)化，從而表征宏觀性能，對于改進微晶玻璃的性能從而指導(dǎo)生產(chǎn)優(yōu)良性能的微晶玻璃具有社會價值。

1 分子動力學(xué)模擬概述

1.1 分子動力學(xué)的發(fā)展歷史

20 世紀50 年代初，ALDER 等[5]用分子動力學(xué)方法設(shè)計了一個具有周期性邊界條件的粒子系統(tǒng)解決有關(guān)動量與壓力的問題，第一次真正意義上從微觀角度模擬了物質(zhì)宏觀性能。 20 世紀70 至80 年代，ANDERSEN[6]對傳統(tǒng)的動力學(xué)方法加以約束條件，開創(chuàng)了恒壓條件下分子動力學(xué)方法。 GILLAN等[7]提出非平衡狀態(tài)動力學(xué)方法。 NOSE[8]提出恒溫條件下分子動力學(xué)方法。 CAR[9]解決了半導(dǎo)體和金屬勢函數(shù)難以模型化的問題。 CAGIN 等[10]提出了基于巨正則系統(tǒng)的分子動力學(xué)方法。 這些方法的提出都極大促進了分子動力學(xué)的的發(fā)展。 進入21世紀，隨著基礎(chǔ)學(xué)科基本理論的完善、計算機硬件的更新?lián)Q代以及考慮更多作用力的多體勢函數(shù)的開發(fā)，使緩慢發(fā)展的分子動力學(xué)得到迅速發(fā)展，已廣泛應(yīng)用于材料、醫(yī)藥、機械、化學(xué)、生物等多個學(xué)科，包含了晶體、非晶體、液態(tài)溶液、復(fù)合結(jié)構(gòu)等方面，尤其對于材料極限條件(超臨界、深過冷、高溫、高壓)的實驗，分子動力學(xué)模擬優(yōu)越性顯著[11-16]。

1.2 分子動力學(xué)模擬的原理

第三種科學(xué)研究手段是除理論分析和實驗觀察之外的分子動力學(xué)方法，也被稱為“計算機實驗”手段[17]。 分子動力學(xué)(Molecular Dynamics，MD)方法和蒙特卡洛(Monte Carlo，MC)方法既有相同點又有所區(qū)別，最初MD 和MC 方法的產(chǎn)生的原因一樣，都是為了計算積分，不同點在于MD 方法在微正則系綜算積分，而MC 方法在正則系綜算積分；與MC 方法相比，MD 方法優(yōu)勢在于計算的內(nèi)容和性質(zhì)更多，得到的信息更多，而且MD 方法可使用的軟件要遠多于MC 方法，適應(yīng)于大多數(shù)體系。 分子動力學(xué)模擬是借助MD 模擬軟件[18](如Moldy、Materials Studio、LAMMPS、DL_POLY 等)在原子及分子水平對多體系統(tǒng)進行求解的計算機模擬方法，在模擬過程遵循兩個假設(shè)[19]:(1) 體系中粒子運動遵循牛頓運動定律；(2) 粒子之間的相互作用滿足疊加定理。忽略量子效應(yīng)和多體作用的模擬研究對象，微觀粒子的運動定義為體系原子核的運動，因此忽略了原子內(nèi)電子的分布情況，這與實際存在的物理系統(tǒng)有一定的差別。

分子動力學(xué)模擬體系內(nèi)的所有粒子的運動都遵循牛頓運動方程，即Fi(t)＝ miai(t) ，F(xiàn)i(t)為粒子i在系統(tǒng)中受到的力；mi為粒子i的質(zhì)量；ai(t)為粒子i的加速度。 勢能函數(shù)U對坐標ri求導(dǎo)，可得對于體系任何一個粒子，都由式(1)和(2)表示為

式中v為速度矢量，對式(1)和(2)求解，可以求得體系中任一粒子的位置和速度。

1.3 有限差分算法

對于多元體系，牛頓運動方程求解有一定的局限性，一般利用有限差分法完成對運動方程的求解，常用的方法包括Verlet 算法、Velocity-Verlet 算法、Beeman 算法、Gear 算法、Rahman 算法、蛙跳算法(Leap-frog algorithm)、Nordsieck 算法。

1.3.1 Verlet 算法

其算法計算過程較簡單，對計算機性能要求不高，每次積分只計算一次力而且時間可逆，因此運用最為廣泛，而其缺點是計算結(jié)果精度不高、軌跡與速度無關(guān)[20]。

1.3.2 Velocity-Verlet 算法

Velocity-Verlet 算法，其特點是可以同時計算出粒子一段時間內(nèi)的位置、速度、加速度，并且準確性較高、計算量適中，已廣泛應(yīng)用于分子模擬中。 其缺點是計算過程相較Verlet 算法復(fù)雜，用時較長[21]。

1.3.3 Beeman 算法

相較于Verlet 算法，Beeman 算法更為復(fù)雜，是一種更為精確的速度表達式，其動能是由速度計算直接得到的，并盡可能得到能量守恒結(jié)果，但是實際應(yīng)用中對模擬條件要求較高，而且計算量較大[22]。

1.3.4 Gear 算法

Gear 算法由泰勒式展開，能夠預(yù)測每一個粒子新位置、速度、加速度[23]。

試驗地位于合作市卡加曼鄉(xiāng)新集村的甘南州農(nóng)科所綜合試驗站，海拔2 721 m，年平均氣溫3.0℃，年降水量639.8 mm左右，無霜期93 d左右，耕種亞高山草甸草原土，旱川地，地力中等，前茬作物為油菜。播種時施尿素150 kg/hm2、磷酸二銨225 kg/hm2作基肥一次性施入，人工犁開溝溜種條播。2017年4月3日播種，6月6日中耕除草，田間管理略高于大田。

1.3.5 Rahman 算法

其算法表達式很復(fù)雜、計算量較大，但是計算結(jié)果精確，不適用于中小型體系模擬。 實際中很少采用Rahman 算法[24]。

1.3.6 蛙跳算法(Leap-frog algorithm)

蛙跳算法由Verlet 算法發(fā)展而來，相較Verlet式帶有t時刻的速度，不計算2y(t) 和y(t - δt) ，減少了計算量，提高了精確度，軌跡與速度無關(guān)。 其缺點是位置項和速度項不能同時計算[25]。 其表達式由式(3)表示為

1.3.7 Nordsieck 算法

Nordsieck 算法適用于解決大系統(tǒng)粒子牛頓方程的求解[26]。

1.4 勢函數(shù)

用來描述系統(tǒng)原子間或者分子間作用力的函數(shù)稱為勢函數(shù)[27]，也稱為力場。 勢函數(shù)不僅決定模擬能否順利進行，還將影響計算結(jié)果的精度。 精準且簡練的勢函數(shù)往往會使模擬效率提高，同時模擬結(jié)果符合物理規(guī)律。 20 世紀中后期，ALDER 等[5]借助數(shù)學(xué)方法，將粒子形象化硬球，一定距離之內(nèi)會發(fā)生彈性碰撞，這種模型雖有缺陷，但是處理問題的思想促進了勢函數(shù)的發(fā)展，勢函數(shù)深刻體現(xiàn)了力場的關(guān)系。 在廣大科研人員共同推動下，出現(xiàn)了多種形式的勢函數(shù)，一般來說，根據(jù)相互作用力包含的粒子數(shù)分為二體勢、三體勢及多體勢，目前對勢、無方向性多體對泛函數(shù)、考慮角度效應(yīng)多體勢是最常見的種類[28]。 常見的勢函數(shù)見表1。

表1 常見的勢函數(shù)表

1.5 系綜

2 微晶玻璃的分子動力學(xué)模擬綜述

2.1 分子動力學(xué)模擬對微觀網(wǎng)絡(luò)結(jié)構(gòu)的分析

微晶玻璃的制備中存在高溫熔融階段，微晶玻璃高溫熔體的微觀網(wǎng)絡(luò)結(jié)構(gòu)特征一直是模擬的重點問題。 需要結(jié)合均方位移函數(shù)(Mean Square Displacement， MSD)、 徑 向 分 布 函 數(shù)( Radial Distribution Function，RDF)、鍵長鍵角分布、配位數(shù)和橋氧數(shù)量等多個角度進行綜合分析，粒子在固定時間段內(nèi)位移平方的平均值為均方位移[31]。 均方位移的量與原子的擴散系數(shù)存在著對應(yīng)關(guān)系，均方位移數(shù)值越大，粒子的擴散系數(shù)越大。 對于研究微晶玻璃網(wǎng)絡(luò)結(jié)構(gòu)，由MSD 可以推測擴散系數(shù)小的粒子很可能為網(wǎng)絡(luò)形成子或者網(wǎng)絡(luò)中間子，反之，則為網(wǎng)絡(luò)修飾子(不參與網(wǎng)絡(luò)結(jié)構(gòu)的構(gòu)建)。 RDF 描述的是某個原子周圍其他原子的分布特征，反映出給定某一粒子的位置坐標，在其半徑r的空間范圍內(nèi)發(fā)現(xiàn)另一原子的概率[32]，用于表征材料的有序程度，利用這一特性，可以探究微晶玻璃體系粒子之間的相互結(jié)合能力大小。 鍵長、鍵角分布規(guī)律能反映出粒子鍵能強弱。 配位數(shù)是原子的第一近鄰原子的個數(shù)，橋氧定義為連接兩個網(wǎng)絡(luò)四面體且同時被兩個四面體所共用的氧粒子[33]。 橋氧數(shù)量的多少反映了硅氧網(wǎng)絡(luò)的完整性、致密性。 在微晶玻璃網(wǎng)絡(luò)結(jié)構(gòu)中，四面體是構(gòu)成網(wǎng)絡(luò)骨架的主要部分，與橋氧數(shù)量類似其也能反映網(wǎng)絡(luò)的致密性[34]。

2.2 分子動力學(xué)在微晶玻璃研究中的應(yīng)用進展

分子動力學(xué)模擬由于獨特的實驗操作環(huán)境，在國內(nèi)發(fā)展較快，受到越來越多的科研人員的青睞，為了更好地了解微晶玻璃微觀結(jié)構(gòu)性質(zhì)，指導(dǎo)特定性能微晶玻璃的生產(chǎn)，開展分子動力學(xué)模擬研究微晶玻璃的性能十分必要，但是目前在微晶玻璃領(lǐng)域的應(yīng)用還較少。 ROSSANO 等[35]利用分子動力學(xué)模擬了CaO-FeO-2SiO2系玻璃，得到了鐵周圍氧原子的徑向分布情況，測量了四配位和五配位鐵原子的氧原子平均距離分別為1.99 和2.15 nm，這些基礎(chǔ)的有關(guān)原子距離的模擬為以后的諸多模擬正確性評價提供了一定的參考價值。 FUMIYA[36]選取NPT 系綜，采用Ewald 求和法模擬了Na2O·3SiO2熔體，發(fā)現(xiàn)Na—O 和Na—Na 的距離隨壓力的增加而縮短，Si—O—Si 鍵角隨壓力的增加而減小，說明了—Si—O—網(wǎng)絡(luò)的崩潰，表明了—Si—O—網(wǎng)絡(luò)的變形自由度增加，網(wǎng)絡(luò)在高壓下趨于崩潰，正由于這些結(jié)構(gòu)的變化，Na2O·3SiO2熔體產(chǎn)生了致密化。 相較大多數(shù)模擬以單一組分對其他基礎(chǔ)玻璃成分的影響，研究壓力變化與微觀網(wǎng)絡(luò)結(jié)構(gòu)的關(guān)系打開了另一種思路，基礎(chǔ)組分固然影響宏觀性能，但是某些特定溫度、壓力變化等外界因素也應(yīng)被重視。 朱才鎮(zhèn)等[37-38]采用了一種二體和三體的多體相互作用勢探究了CaO-Al2O3-SiO2系微晶玻璃成分比例組成、微觀結(jié)構(gòu)之間的關(guān)系，發(fā)現(xiàn)Ca/Al＝1/2 時，CaO-Al2O3-SiO2系微晶玻璃并不是傳統(tǒng)理論中認為的完整的網(wǎng)絡(luò)結(jié)構(gòu)，而是存在一定的非橋氧，Ca/Al＜1/2 時，Si比Al 更容易形成網(wǎng)絡(luò)中間體。 通過模擬了不同Ca、Al 元素比情況下網(wǎng)絡(luò)結(jié)構(gòu)的特征變化，發(fā)現(xiàn)了橋氧的存在，這是對傳統(tǒng)硅酸鹽網(wǎng)絡(luò)結(jié)構(gòu)完整的觀點的一次否定，體現(xiàn)了計算模擬用來解決現(xiàn)實問題、提供理論創(chuàng)新的作用。 吳永全等[39-40]采用了BMH勢函數(shù)的基礎(chǔ)上，利用經(jīng)典分子動力學(xué)模擬xCaO-(1-x)Al2O3(x為成分的變化)高溫熔體，重點探討了有關(guān)Al 的配位數(shù)及其網(wǎng)絡(luò)中的構(gòu)成、微結(jié)構(gòu)單元分布等結(jié)構(gòu)性質(zhì)，發(fā)現(xiàn)Al 的微觀結(jié)構(gòu)單元是四面體的形式，類似硅酸鹽網(wǎng)絡(luò)中Si—O 鍵的結(jié)構(gòu)特點，證明了Al 的配位數(shù)為4，Al 在網(wǎng)絡(luò)架構(gòu)中起到了網(wǎng)絡(luò)形成子的作用。 可以看出，大多數(shù)的模擬是以基礎(chǔ)玻璃、三元體系為主，這是因為元素越多，需要考慮的作用力越多，開發(fā)的勢函數(shù)越復(fù)雜，因此分子動力學(xué)研究領(lǐng)域中，開發(fā)簡潔而且精確的勢函數(shù)是重要的方向。 王艷偉等[41]采用EAM 原子勢，在NPT 系綜條件下模擬了Ni 的含量變化對Ni-Zr 非晶玻璃微觀結(jié)構(gòu)的影響，發(fā)現(xiàn)ico＋other(ico 為體系的十二面結(jié)構(gòu)，other 為未知的配位結(jié)構(gòu))結(jié)構(gòu)隨著Ni 數(shù)量的增加先增加后減少，ico＋other 結(jié)構(gòu)反映了非晶的形成能力，當Ni 含量為50%時，結(jié)構(gòu)數(shù)量最多，即說明了一定數(shù)量的Ni 有利于促進玻璃的形成。 趙亞賢等[42]利用經(jīng)典分子動力學(xué)模擬對比了Na-Al-Si系微晶玻璃和高堿高鋁酸鹽玻璃中堿金屬的擴散行為，統(tǒng)計了兩種玻璃的Si—O、Al—O 鍵長，擴散系數(shù)，發(fā)現(xiàn)高堿高鋁酸鹽玻璃中堿金屬離子的擴散系數(shù)更大，這是因為[AlO4]四面體的存在，相較于[SiO4]四面體體積更大，會使網(wǎng)絡(luò)結(jié)構(gòu)疏松，有利于堿離子擴散，證明了高含量的Al 為堿金屬離子的擴散提供了更多的擴散通道。 王海龍等[43]采用Mishin 嵌入原子勢模擬了金屬玻璃中的Cu 在急速冷條件下的形成過程，利用Verlet 算法求解運動方程，證明了隨著應(yīng)變率的提高，Cu 的塑性流動應(yīng)力增大，晶化程度增加，應(yīng)力效應(yīng)和溫度效應(yīng)共同作用導(dǎo)致了金屬玻璃的晶化，對于解釋金屬的晶化起到了一定的促進作用。 肖成[44]構(gòu)建了CaO-SiO2-Al2O3-Na2O四元體系的微觀玻璃模型，建立了1 873 K下體系中溫度和黏度的函數(shù)關(guān)系，發(fā)現(xiàn)Al2O3濃度高的體系傾向于形成更加復(fù)雜的網(wǎng)絡(luò)結(jié)構(gòu)。 BINGHUI 等[45]巧妙設(shè)計了一種自上而下的模擬方法，探究影響Al2O3-SiO2系微晶玻璃斷裂韌性的因素，發(fā)現(xiàn)并證明了增強程度隨著納米晶體的尺寸、長徑比和排列角度的變化而變化。 王亞文[46]模擬了形核劑在熔融態(tài)高爐渣中的擴散行為，高爐渣成分與微晶玻璃成分相似，非均勻形核過程需要加入晶核劑降低析晶活化能。 通過數(shù)值與物理模擬發(fā)現(xiàn)形核劑的混勻時間與形核劑密度、等效直徑成正比，與熔體黏度成反比。 可以看出，有關(guān)微晶玻璃的模擬已經(jīng)取得了一些進展，從鍵長、鍵角的分布，溫度、壓力等外界因素對網(wǎng)絡(luò)結(jié)構(gòu)的影響，以及各元素在網(wǎng)絡(luò)中的位置、作用等都做出了一定的解釋，這些研究為制備性能優(yōu)異的微晶玻璃提供了一定的幫助。

井敏課題組在微晶玻璃的制備方面經(jīng)過多年的積累，已有一套成熟制備微晶玻璃的工藝。 對于微晶玻璃微觀結(jié)構(gòu)的探究，引進了先進的材料工作室(Materials Studio，MS)軟件開展工作，取得了一系列成果。 張國瑩等[47]基于MS 軟件利用分子動力學(xué)模擬了Na2O-Al2O3-SiO2微晶玻璃體系熔體結(jié)構(gòu)，發(fā)現(xiàn)了作為網(wǎng)絡(luò)游離子的Na＋能夠促進Al3＋從網(wǎng)絡(luò)修飾子變成網(wǎng)絡(luò)形成子，且在一定范圍內(nèi)，隨著Na＋的增加，Si—O—Si 和Al—O—Al 鍵角的分布范圍變小。 王健健等[48-49]探究了Fe3＋在CaO-Al2O3-SiO2系微晶玻璃中的行為，發(fā)現(xiàn)隨著Fe3＋含量的增加Al—O 鍵長分布更廣，峰型有單峰向雙峰轉(zhuǎn)化的趨勢。 王正[50]針對RO/R2O-Al2O3-SiO2熔體黏度特性進行了相關(guān)模擬工作，采用BMH 勢，選用Shear模塊計算了Na-Al-Si 系和Ca-Al-Si 系基礎(chǔ)玻璃的剪切黏度。 發(fā)現(xiàn)隨著堿金屬元素的減少，Na-Al-Si系黏度增加要大于Ca-Al-Si 系，整體結(jié)構(gòu)上，Na-Al-Si 系相較Ca-Al-Si 系網(wǎng)絡(luò)結(jié)構(gòu)較松散。 張永豪等[51]結(jié)合BMH 勢，建立了MgO-Al2O3-SiO2-TiO2系微晶玻璃的微觀結(jié)構(gòu)模型，分別探究了TiO2晶核劑含量和SiO2含量對微觀結(jié)構(gòu)的影響，發(fā)現(xiàn)隨著SiO2添加量的增大，體系中Mg2＋的均方位移逐漸減小，當SiO2的含量為57.1%時，橋氧總量最多，硅酸鹽網(wǎng)絡(luò)結(jié)構(gòu)最完善。 TiO2含量增加會使Al—O 鍵長的分布更廣，峰形變寬。

3 展望

目前，采用分子動力學(xué)模擬技術(shù)對于微晶玻璃體系的模擬已經(jīng)取得很大的進展。 從微觀角度出發(fā)，通過原子間相互作用力、邊界條件的選取、熱力學(xué)性質(zhì)、運動方程的計算對微晶玻璃宏觀的性質(zhì)(黏度、致密度、硬度等)做出了解釋，改進了微晶玻璃的性能。 但是仍存在一些問題，如初始模型過于理想化，而實際晶體材料存在各種缺陷，應(yīng)當深入研究構(gòu)建有缺陷的建模方法；MD 方法有觀測時間和系統(tǒng)大小的限制，仿真的規(guī)模較小，原子數(shù)較少，缺乏普遍性，應(yīng)當發(fā)展較大規(guī)模的模擬；微晶玻璃通常使用經(jīng)典分子動力學(xué)模擬，忽略了電子極化效應(yīng)，對于電荷的相關(guān)信息無法獲取，勢函數(shù)過于依賴力場參數(shù)，數(shù)量較少，更加精確和簡潔的勢函數(shù)是發(fā)展分子動力學(xué)的重點；微晶玻璃模擬的模型也較為簡單，原子種類和數(shù)量較少，大多數(shù)模擬忽略了含量稀少的元素，或許會影響宏觀性能。 因此，多元體系的微晶玻璃的分子動力學(xué)模擬是發(fā)展的方向之一。