冷卻速率對Lennard-Jones體系凝固過程中結構與動力學性質的影響

楊邦鍫，孫永麗，孫亞娟，楊 平，厲華明

(1. 太原理工大學物理與光電工程學院，太原 030024； 2. 天津工業大學理學院，天津 300160)

冷卻速率對Lennard-Jones體系凝固過程中結構與動力學性質的影響

楊邦鍫1，孫永麗1，孫亞娟2，楊 平1，厲華明1

(1. 太原理工大學物理與光電工程學院，太原 030024； 2. 天津工業大學理學院，天津 300160)

用分子動力學模擬方法研究了五種不同冷卻速率對Lennard-Jones體系凝固過程中結構與動力學性質的影響.采用兩種不同的方法來確定玻璃轉變溫度Tg，并且對結晶溫度Tc、徑向分布函數g(r)、均方位移函數MSD與擴散系數D、平均配位數進行比較分析.結果表明：冷卻速率影響Lennard-Jones體系凝固過程中的結構.當使用足夠高的冷卻速率冷卻時，體系發生玻璃化轉變，而且冷卻速率越快，玻璃轉變溫度越高；當冷卻速率較小時，體系形成晶體，而且冷卻速率越慢，結晶溫度越高，結晶程度也越高.同時發現，冷卻速率對擴散系數和平均配位數也有很大影響，二者在體系發生玻璃轉變時都有一個緩變的過程，表明了過冷液相區的存在.

Lennard-Jones體系； 冷卻速率； 結構； 動力學性質

1 引 言

分子動力學模擬自產生以來得到廣泛的應用，特別隨著計算機技術的高速發展，更是得到了絕大多數研究者的青睞.分子動力學模擬可以研究各種類型的材料，比如高分子材料[1-3]、金屬材料[4-6]、非晶體合金[7-9]等.近年來關于液態金屬凝固過程中微觀結構的演變特性引起了研究者的廣泛關注和重視，而冷卻速率對其微觀結構的影響也引起了激烈的討論[10,11].其中液態金屬Cu的快速凝固的模擬研究取得了許多重大的進展，這些研究所采用的模型勢有FS勢[12]、TB勢[13]、EAM勢[14,15]、由Wang等人[16]所發展的擴展非局域模型贗勢[17]、Quantum Sutton-Chen (QSC)[18,19]多體勢[20]等.由此可見，即使研究同一種金屬材料，也能采用不同的勢函數進行分析.

為此，本文采用Lennard-Jones模型，對液態Lennard-Jones體系在五個不同冷卻速率下的凝固過程進行模擬研究，分析其結晶溫度Tc與玻璃轉變溫度Tg、徑向分布函數g(r)、均方位移函數MSD與擴散系數D、平均配位數，以較為深入的研究冷卻速率對Lennard-Jones體系凝固過程中結構與動力學性質的影響.

2 模擬條件與方法

模擬計算的條件為：將4000個Lennard-Jones原子置于一立方體盒中, 然后讓體系在周期性邊界條件下運行.Lennard-Jones原子間的相互作用勢為

(1)

式中，r表示原子間的距離.本文采用對比單位(Reduced unit s) , 用*表示，并以原子半徑σ、原子間相互作用能量和原子質量m作為基本單位.

模擬計算均從T*=1.0開始，首先讓體系等溫等壓運行50000 步使之處于平衡態，再分別以9.672×1011K/s，4.836×1011K/s，1.209×1011K/s，3.023×1010K/s，1.209×1010K/s的冷卻速率冷卻至T*=0.1，其中溫度每隔0.05讓體系進行一次充分弛豫，以測量該體系的結構組態, 即記錄下每個原子的空間坐標.根據這些記錄，對五個不同冷卻速率下各溫度的徑向分布函數進行對比分析，并且分析其它相關量，以探討不同冷卻速率對Lennard-Jones體系快速凝固過程中結構與動力學性質的影響.

3 模擬計算結果與分析

3.1 結晶溫度Tc與玻璃轉變溫度Tg分析

在不同冷卻速率下，液體在凝固后可以得到晶體與非晶體，生成晶體的溫度被定義為結晶溫度Tc，生成非晶的溫度被定義為玻璃轉變溫度Tg.圖1 給出了Lennard-Jones體系在不同冷卻速率下平均原子體積隨溫度的變化關系.可以看出，在冷卻速率為9.672×1011K/s，4.836×1011K/s與1.209×1011K/s時，在降溫過程中平均原子體積沒有急劇的減少，說明在這三種冷卻速率下，Lennard-Jones體系形成了非晶結構.反之，在3.023×1010K/s與1.209×1010K/s的冷卻速率下，平均原子體積出現明顯的轉折，說明體系已經形成晶體.此外，從圖中還可以看出，在形成晶體的兩種冷卻速率下，平均原子體積出現的轉折點不一致，且速率較小的轉折點所對應的溫度較高，說明冷卻速率越慢，結晶溫度Tc越高.

圖1 Lennard-Jones體系在不同冷卻速率下平均原子體積隨溫度的變化關系Fig. 1 The average atomic volume as a function of temperature for Lennard-Jones system at different cooling rates

玻璃轉變溫度Tg是判定非晶形成能力的重要參數(Trg=Tg/Tl)[21]，因此對于玻璃轉變溫度的測量是一個關鍵問題.下面采用兩種常用的方法來確定玻璃轉變溫度Tg.

(1)體積變化判定法

這種方法是根據液體在快速冷卻過程中，發生玻璃化轉變時體積與溫度的變化關系，將曲線兩端的直線部分外推至相交[22,23]，交點即為Tg，如圖2(a)所示(紅色箭頭所指即為玻璃轉變溫度Tg).由此可見，采用這種方法來計算玻璃轉變溫度所需要的數據比較容易得到，而且在實驗測定和計算模擬都可以應用.

(2)徑向分布函數參量法

這種方法是 Wendt-Abraham等人[24]提出的一種根據非晶化轉變前后的結構變化來確定玻璃轉變溫度Tg的方法.這種方法首先定義一個參量R=gmin/gmax(g(r)曲線第一波谷的谷值與第一波峰的峰值之比)，然后通過不同溫度下的R值與溫度的關系，得到斜率不同的兩條直線，其中的斜率變化點，即是玻璃轉變溫度Tg，如圖2(b)所示(黑色箭頭所指即為玻璃轉變溫度Tg).

通過比較分析圖2(a)和圖2(b)，可以發現這兩種方法所確定的玻璃轉變溫度Tg相差不大.同樣的，我們可以利用這兩種方法分別對冷卻速率為4.836×1011K/s與1.209×1011K/s的情況進行分析，確定各自的玻璃轉變溫度Tg.結果表明，當使用足夠高的冷卻速率進行冷卻時，Lennard-Jones體系在快速凝固后形成非晶，而且冷卻速率越快，玻璃轉變溫度Tg越高.

圖2 冷卻速率為9.672×1011 K/s時Lennard-Jones體系的平均原子體積和徑向分布函數參量R隨溫度的變化關系Fig. 2 The average atomic volume and the radial distribution function parameter R as a function of temperature for Lennard-Jones system at the cooling rate of 9.672×1011 K/s

3.2 徑向分布函數g(r)分析

徑向分布函數g(r)被廣泛用來描述液態和非晶態的結構特征，它表示距離原子r位置處發現另一原子的幾率.圖3(a)、(b)分別給出了液態下Lennard-Jones體系在冷卻速率為9.672×1011K/s，1.209×1010K/s時的徑向分布函數g(r*)曲線隨溫度的變化情況.從圖中可以看出，隨著溫度的降低，第一個波峰由駝峰逐漸變高變銳.根據徑向分數函數的定義，這種情況說明了距離原子r位置處發現另一原子的幾率越來越大，表明了體系內相鄰原子成鍵的幾率也越來越大，有序度逐漸加強.在冷卻速率為9.672×1011K/s的快速冷卻下，當溫度降到T*=0.4時，g(r*) 曲線的第二峰開始出現劈裂, 而且隨著溫度的降低，第二峰劈裂的越來越明顯，并且劈裂的兩個峰是前高后低.這正是非晶態金屬結構的重要特征之一，這也表明快速冷卻到T*=0.4附近就開始形成非晶.

從圖3(b)可見，g(r*)曲線從T*=0.4開始，在第二峰與第三峰前后，明顯出現了多個無規律的小峰，與冷卻速率為3.023×1010K/s時的相比，這個速率下各小峰出現的位置和形狀更有規律.說明體系在冷卻速率為1.209×1010K/s時的結晶有序度比冷卻速率為3.023×1010K/s時的要高.

圖3 各冷卻速率下徑向分布函數g(r*)曲線隨溫度的變化關系 (a) 9.672×1011 K/s；(b) 1.209×1010K/sFig. 3 The radial distribution function g(r*) as a function of temperature at different cooling rates (a) 9.672×1011 K/s；(b) 1.209×1010 K/s

為了比較分析不同冷卻速率下體系形成的固體結構，圖4分別給出了不同冷卻速率下各個溫度的徑向分布函數的對比關系.由于文中采用了五種冷卻速率，有三種速率使系統玻璃化，另外兩種晶化，所以根據這種這種情況將冷卻速率分為兩組，如圖4(a)、(b)所示.由圖可見，在T*=0.5時徑向分布函數g(r*)在不同冷卻速率下幾乎完全一致，說明在液態時的g(r*)與冷卻速率無關.比較圖4(a)和圖4(b)，可以發現圖4(b)中的g(r*)相對于圖4(a)的較敏感，而且圖4(a)的三種冷卻速率下g(r*)曲線在體系玻璃化后幾乎完全重合，說明使系統玻璃化的不同冷卻速率對體系玻璃化所形成的徑向分布函數無多大影響.但對于使系統結晶的不同冷卻速率來說，由圖4(b)可見，溫度從T*=0.4開始，兩種冷卻速率下的g(r*)均有明顯的第一峰和第二峰.而且隨著冷卻速率的降低，第一峰和第二峰在強度上都有所升高，但第一峰的寬度有所減小，致使第一近鄰配位數減少,有序度增強.在第一峰和第二峰之間均出現了肩峰，說明在液態體系在冷卻過程中形成了一種中程有序結構(MRO)[25].

3.3 均方位移函數MSD與擴散系數D分析

均方位移(Mean square displacement，簡稱MSD)描述了體系中原子在t時間內的平均位移的平方.其表達式為

MSD=〈|ri(t+t0)-ri(t0)|2〉

(2)

式中，ri(t0)表示原子在t0時刻的位置，〈...〉表示系綜平均.根據Einstein關系[26]，它與擴散系數滿足如下關系式

(3)

圖4 不同冷卻速率下各個溫度的徑向分布函數的對比關系Fig. 4 Comparison of the radial distribution function of each temperature at different cooling rates

其中，c是常數，D是擴散系數.由此可知擴散系數D可以通過均方位移隨時間的變化曲線的斜率計算得到.

圖5給出了在冷卻速率為4.836×1011K/s與3.023×1010K/s時擴散系數D隨溫度的變化關系.由圖可見，溫度從T*=1下降到T*=0.45，兩種冷卻速率下lnD隨1/T*的變化曲線完全重合，而且呈線性關系，說明原子擴散系數D隨溫度T*的變化滿足Arrhenius關系，即lnD～1/T*.對于冷卻速率為4.836×1011K/s的曲線來說，溫度下降至T*=0.45之后，D隨T*的變化偏離Arrhenius關系，且有一個緩變的過程，說明體系進入了過冷液相區.待溫度下降至T*=0.4后，兩種冷卻速率下的擴散系數D基本不隨溫度的變化，說明此時體系流動性很小，已經凝固成晶體或非晶體.

圖5 在冷卻速率為4.836×1011 K/s與3.023×1010 K/s時擴散系數D隨溫度的變化關系Fig. 5 The diffusion coefficient D as a function of temperature at cooling rates of 4.836×1011 K/s and 3.023×1010 K/s

3.4 平均配位數的分析

配位數是指體系中任一原子周圍最近鄰且距離近于相等的原子數目[27]，是表征物質結構的基本參數.我們在模擬研究中發現，體系在冷卻過程中，其原子的平均配位數(第一近鄰)隨溫度的變化較為敏感.由此我們給出了兩種冷卻速率下(冷卻速率分別為9.672×1011K/s與1.209×1010K/s)平均配位數隨溫度的變化關系，如圖6所示.

由圖可以明顯看出，在兩種冷卻速率下，體系的原子平均配位數隨溫度的變化在高溫時差別很小，而且都隨著溫度的降低而增大，而當溫度降到T*=0.4以下時則有明顯的不同.以較快的冷卻速率(9.672×1011K/s)冷卻時，體系在凝固形成非晶的過程中，其平均配位數先是比較穩定而緩慢地增加，其后保持不變，說明了非晶態的形成；而以較慢的冷卻速率(1.209×1010K/s)冷卻時，體系的原子平均配位數在T*=0.4以下基本保持不變，說明晶體已經形成.由此可見，平均配位數可以在某種程度上反映出Lennard-Jones體系在凝固過程中微觀結構的演變特性，因而它可以作為一種方法來研究體系的結晶與玻璃化轉變.

4 結 論

根據上述五種不同冷卻速率對液態Lennard-Jones體系凝固過程中結構與動力學性質的影響的模擬研究結果和討論，可得如下結論：

(1)冷卻速率對結晶溫度與玻璃轉變溫度的影響.當冷卻速率為3.023×1010K/s與1.209×1010K/s時，體系形成晶體，而且冷卻速率越慢，結晶溫度Tc越高；當冷卻速率為9.672×1011K/s，4.836×1011K/s與1.209×1011K/s時，體系發生玻璃化，而且冷卻速率越快，玻璃轉變溫度Tg越高.

(2)冷卻速率對Lennard-Jones體系形成的最終固體結構的影響.當使用足夠高的冷卻速率進行冷卻時，冷卻速率對Lennard-Jones體系發生玻璃化后所形成的徑向分布函數無多大影響；反之，當使用的冷卻速率過慢致使Lennard-Jones體系結晶時，冷卻速率對Lennard-Jones體系發生晶化后所形成的徑向分布函數有很大影響，而且冷卻速率越慢，有序度越強，結晶程度也越高.

(3)冷卻速率對擴散系數的影響.當Lennard-Jones體系在高溫液體時，原子擴散系數D隨溫度T*的變化滿足Arrhenius關系，即lnD～1/T*，與冷卻速率無關；而當冷卻速率足夠高致使Lennard-Jones體系形成非晶時，原子擴散系數D不會隨溫度T*發生突變，而是有一個緩變的過程，說明了過冷液相區的存在.

(4)冷卻速率對平均配位數的影響.不同冷卻速率下Lennard-Jones體系的原子平均配位數在剛開始降溫階段都會隨溫度的降低而增大.而在高溫液體之后，若使用的冷卻速率過快致使Lennard-Jones體系形成非晶，那么其平均配位數先是比較穩定而緩慢地增加，而后保持不變；若使用的冷卻速率過慢致使Lennard-Jones體系結晶，那么在其結晶后平均配位數一直保持不變，沒有一個緩變的過程.

[1] Yao D G, Kim B Y. Simulation of the filling process in micro channels for polymeric materials [J].J.Micromech.Microeng., 2002, 12(5): 604.

[2] Qiao R, Brinson L C. Simulation of interphase percolation and gradients in polymer nanocomposites [J].Compos.Sci.Technol., 2009, 69(3): 491.

[3] Liu Z S, Hong W, Suo Z G,etal. Modeling and simulation of buckling of polymeric membrane thin film gel [J].Comp.Mater.Sci., 2010, 49(1): S60.

[4] Yoshidaa F. Material models for accurate simulation of sheet metal forming and springback [J].AIPConf.Proc., 2010, 1252: 71.

[5] Taherizadeh A, Green D E, Ghaei A,etal. A non-associated constitutive model with mixed iso-kinematic hardening for finite element simulation of sheet metal forming [J].Int.J.Plasticity, 2010, 26(2): 288.

[6] Toni M D, Coudert F X, Paranthaman S,etal. Molecular simulation of a Zn-triazamacrocyle metal-organic frameworks family with extraframework anions [J].J.Phys.Chem. C, 2012, 116(4): 2952.

[7] Schenk T, Holland M D, Simonet V,etal. Icosahedral short-range order in deeply undercooled metallic melts [J].Phys.Rev.Lett., 2002, 89(7): 075507.

[8] Luo W K, Sheng H W, Alamgir F M,etal. Icosahedral short-range order in amorphous alloys [J].Phys.Rev.Lett., 2004, 92(14): 145502.

[9] Huang L, Wang C Z, Hao S G,etal. Short- and medium-range order in amorphous Zr2Ni metallic alloy [J].Phys.Rev. B, 2010, 81(9): 094118.

[10] Yi X H, Liu R S, Tian Z A,etal. Formation and evolution properties of clusters in liquid metal copper during rapid cooling processes [J].Trans.NonferrousMet.Soc.China, 2008, 18(1): 33.

[11] Wang J F, Huang S, Guo S F,etal. Effects of cooling rate on microstructure, mechanical and corrosion properties of Mg-Zn-Ca alloy [J].Trans.NonferrousMet.Soc.China, 2013, 23(7): 1930.

[12] Zhang T, Zhang X R, Ding S L. The molecular dynamics simulation of cooled liquid Cu by FS potential model [J].J.At.Mol.Phys., 2003, 20 (3): 357 (in Chinese) [張弢, 張曉茹, 丁世良. 用FS多體勢模型模擬金屬銅的冷卻過程[J]. 原子與分子物理學報, 2003, 20(3): 357]

[13] Liu C S, Xia J C, Zhu Z G,etal. The cooling rate dependence of crystallization for liquid copper: A molecular dynamics study [J].J.Chem.Phys., 2001, 114(17): 7506.

[14] Chen F F, Zhang H F, Qin F X,etal. Molecular dynamics study of atomic transport properties in rapidly cooling liquid copper [J].J.Chem.Phys., 2004, 120(4): 1826.

[15] Pang H, Jin Z H, Lu K. Relaxation, nucleation, and glass transition in supercooled liquid Cu [J].Phys.Rev. B, 2003, 67(9): 094113.

[16] Wang S, Lai S K. Structure and electrical resistivities of liquid binary alloys [J].J.Phys. F:Met.Phys. , 1980, 10(12): 2717.

[17] Lin Y, Liu R S, Tian Z A,etal. Effect of cooling rates on microstructures during solidification process of liquid metal Zn [J].ActaPhys.Chim.Sin. , 2008, 24(2): 250 (in Chinese)

[18] Qi Y, Cagin T, Kimura Y,etal. Molecular-dynamics simulations of glass formation and crystallization in binary liquid metals: Cu-Ag and Cu-Ni [J].Phys.Rev. B, 1999, 59(5): 3527.

[19] Sutton A P, Chen J. Long-range finnis-sinclair potentials [J].Philos.Mag.Lett., 1990, 61(3): 139.

[20] Yi X H, Liu R S, Tian Z A,etal. Simulation study of effect of cooling rate on evolution of microstructures during solidification of liquid metal Cu [J].ActaPhys.Sin., 2006, 55(10): 5386 (in Chinese)

[21] Turnbull D. Under what conditions can a glass be formed? [J].Contemp.Phys., 1969, 10(5): 473.

[22] Berthier L, Biroli G. Theoretical perspective on the glass transition and amorphous materials [J].Rev.Mod.Phys., 2011, 83(2): 587.

[23] Biroli G, Garrahan J P. Perspective: The glass transition [J].J.Chem.Phys., 2013, 138(12): 12A301.

[24] Abraham F F. An isothermal-isobaric computer simulation of the supercooled liquid/glass transition region: Is the short-range order in the amorphous solid fcc? [J].J.Chem.Phys., 1980, 72(1): 359.

[25] Li H. Shoulder-peak formation in the process of quenching [J].Phys.Rev. B, 2003, 68(2): 024210.

[26] Barkai E, Fleurov V N. Generalized Einstein relation: A stochastic modeling approach [J].Phys.Rev. E, 1998, 58(2): 1296.

[27] Zhou L L, Liu R S, Hou Z Y,etal. Simulation study of effects of cooling rate on evolution of micro-cluster structures during solidification of liquid Pb [J].ActaPhys.Sin., 2008, 57(6): 3653 (in Chinese)

Effects of cooling rates on structures and kinetic properties during solidification process of the Lennard-Jones system

YANG Bang-Qiao1, SUN Yong-Li1, SUN Ya-Juan2, YANG Ping1, LI Hua-Ming1

(1.College of Physics and Optoelectronics, Taiyuan University of Technology, Taiyuan 030024, China;2. School of Science, Tianjin Polytechnic University, Tianjin 300160, China)

A molecular dynamics simulation study was performed for the effects of five different rates on structures and kinetic properties during solidification process of the Lennard-Jones system. Glass transition temperatureTgwas measured by using two different methods, and the crystallization temperatureTc, the radial distribution function, the mean square displacementMSD, the diffusion coefficientDand the mean coordination number were used to compare and analyze. The results show that the cooling rate plays a critical role on structures during solidification process of the Lennard-Jones system. As the cooling rate being high enough, the system will be transferred into glass, and the greater the cooling rate, the higher the glass transition temperature. As the cooling rate being smaller, the system will be crystallized, and the slower the cooling rate, the higher the crystallization temperature, and the higher the degree of crystallization. At the same time, it is found that the diffusion coefficient and the mean coordination number are largely influenced by cooling rates, and both have a slow graded process of the glass transition occurred in the system, which shows the existence of supercooled liquid region.

Lennard-Jones system; Cooling rate; Structure; Kinetic property

國家自然科學基金項目(11204200)

楊邦鍫(1989—),男，浙江溫州人，碩士研究生，研究方向為分子動力學模擬.E-mail: daxuewuli001@163.com

孫永麗. E-mail: sunyongli@tyut.edu.cn

103969/j.issn.1000-0364.2015.08.022

O751

A

1000-0364(2015)08-0647-06

投稿日期:2013-12-08