Cu66Ti34非晶合金凝固过程的分子动力学模拟

邱克强, 厉 虹, 孙 晶, 尤俊华, 任英磊, 李庆丰

(沈阳工业大学 材料科学与工程学院，沈阳 110870)

块体Cu基非晶合金自从被报道以来[1]，因其具有优异的力学性能和相对较低的材料成本而受到关注。目前通过成分设计可以获得多种 Cu基块体非晶合金[2−6]，其中Cu50Zr50二元合金可以形成直径为2 mm的非晶合金[7−9]。考虑Ti具有较低的成本和与Zr类似的性质，因此，Cu-Ti二元非晶合金已成为研究热点。但是，目前的实验结果表明，Cu-Ti二元合金还不具有制备成块体非晶合金的能力[10−11]。

目前，基于分子动力学(MD)方法对Cu-Zr二元系非晶合金原子团簇结构演化的模拟比较多[12−13]，而从凝固过程出发，对Cu-Ti二元非晶合金性能的模拟还很少，本文作者基于 WADLEY等[14]和ZHOU等[15]提出的普适嵌入原子模型(GEAM 势)，从动力学角度模拟过冷液态合金在不同温度下原子的均方位移(MSD)、偶关联函数和焓随温度的变化关系，进一步给出定压比热容与温度的函数关系式，从而判定玻璃转变的发生。

1 计算方法

1.1 势函数模型

GEAM势的表达形式如下：

式中：Fni、Fi、Fe和n为模型参数，可由结合能和体模量计算得到；ρe为平衡时的电子密度。

两体作用势可写成如下形式：

式中：r为原子间的距离；re为平衡时最近邻原子间的距离；A、B、α和β为模型参数；κ、λ为有关截断半径的附加参数。当计算合金时，合金势采用以下形式[16]：

式中：φab(r)为合金的两体势；φaa(r)和φbb(r)分别为a和b组分的两体势；fa(r)和 fb(r)分别为a和b组分的电子密度函数。Cu66Ti34二元合金的势参数如表1所列。

1.2 模拟方法

模拟过程中采用二元合金Cu66Ti34，将2 000个原子，按照比例置于立方体盒子中，并让系统在三维周期边界条件下以4×1013K/s冷却速率进行冷却，时间步长选定为2.5×1015s。模拟过程中采用NTP系综，首先让系综在1 600 K等温运行80 000步以使系综达到平衡态，然后以给定冷却速率逐步冷却到300 K，最后在300 K弛豫40 000步，其中每100 K让系统等温运行2 000步，用于系统的数据统计，从而算出体系的各种性质参数。

表1 Cu66Ti34 二元合金的势参数Table 1 Potential parameters of binary Cu66Ti34 alloy

2 计算结果及讨论

2.1 偶关联函数

偶关联函数被广泛用来描述液态和非晶态的结构特征，是一个重要的形态谱参数，通常定义如下[17]：

式中：gi,j(r)表示以a原子为中心、在距离为r到r+Δr的球壳范围内发现原子的概率；L为模拟元胞的边长；Ni和Nj分别是原子i和j的个数；nαj是指在 i原子周围r到r+Δr范围内发现j原子的个数；Δr是计算的步长。

图1所示为Cu66Ti34在4×1013K/s冷却速率下的偏偶关联函数和总偶关联函数的图像。从图1可以看出，高温时偶关联函数都显示出典型的液态结构，而在低温300 K时，gCu-Cu(r)、gTi-Ti(r)、gCu-Ti(r)和gtot(r)第二峰均发生分裂，标志着非晶结构的形成。由于总的偏偶关联函数 gtot是所有原子偏偶关联函数三维信息的平均[18]，因此，总偶关联函数的第二峰分裂不如偏偶关联函数明显。与高温液态结构的 g(r)相比较，低温时第一峰的宽度变窄、高度增加、变锐，这表明随着温度的降低，原子配位数增加，短程有序性增强；另一方面，非晶结构的偶关联函数在第二峰过后和液态结构的无太大区别，偶关联函数都是比较平缓的曲线，说明仍然是长程无序的。从3种偏偶关联函数来看，gCu-Ti(r)的第一峰都高于 gCu-Cu(r)和 gTi-Ti(r)的，这说明合金中异类原子之间的相互作用较强，容易形成非晶[19]。在600 K时，总偶关联函数就开始出现比较微弱的第二峰分裂，随着温度的降低，此现象表现更加明显。

2.2 均方位移(MSD)

图1 Cu66Ti34二元合金不同温度下偶关联函数曲线Fig.1 Pair correction function curves of Cu66Ti34 binary alloy at different temperatures: (a)gCu-Cu; (b)gTi-Ti; (c)gCu-Ti; (d)gtot

粒子位移平方的平均值称为均方位移，均方位移随时间的变化表征了液态金属粒子的扩散行为。为了便于观察Cu66Ti34合金中Cu和Ti原子在急冷到各个温度下MSD的图像变化，以900 K为分界，将图像分开展示，图2(a)和(c)所示分别为Cu66Ti34合金中Cu和 Ti原子在急冷到 1 500、1 400、1 300、1 200、1 100、1 000和900 K时MSD(图中以D表示)的变化；图2(b)和(d)所示分别为Cu和Ti原子急冷到800、700、600、500、400和300 K时MSD随时间的变化情况。在800 K之前，它们都表现出液态的特点，尤其是在1 200 K之前，MSD曲线的斜率较大，说明原子具有较强的扩散能力。从图2(b)和(d)可以看出，800 K时，MSD曲线都变得比较复杂，从 MSD的数值看到，Cu和 Ti在800 K时MSD的最大值均已小于1×10−2nm2，合金熔体已经很粘稠。600 K时，在0.25 ps后，曲线的斜率降低，此时可认为扩散受阻，同时，也是动力学明显呈现冻结状态的起始温度。到500 K时，原子扩散能力很弱，但局域结构仍可以表现出随时间的波动性。到300 K时，原子几乎没有扩散，说明原子已经找到它们的平衡位置。通过图2(a)与(c)，图2(b)与(d)中 MSD的数值比较，还可以看出 Cu的扩散能力比Ti的强。

为了探讨原子的扩散性与玻璃转变温度之间的关系，图3给出了偶分布函数第一谷的最小值与第一峰最大值之比R[20−21]的分布规律，即R=g(r)min/g(r)max。在冷却速率为4×1013K/s条件下，测得玻璃化转变温度Tg约为600 K，略低于Cu50Ti50的玻璃化转变温度641 K[22]，考虑到Cu50Ti50具有较高的固相温度，因此预测 Cu66Ti34非晶合金玻璃化转变温度在600～641 K是符合实际的。

2.3 热力学性质分析

在冷却速率为4×1013K/s条件下，Cu66Ti34合金熔体凝固过程中焓随温度的变化如图4所示，总体上，合金的生成焓随温度的降低而降低。这说明液态金属过冷温度越低，合金中所含热量越少，越有利于非晶的形成。模拟得到的数据在400～900 K之间有波动，说明在这个温度范围内，此体系应处于过冷态，结构变化比较复杂，玻璃化转变在这个温度范围内发生。能量变化较为连续，只是在一定温度以下曲线斜率减小，这是非晶形成的标志[23]。当凝固过程中没有发生结晶时，这一过程表现在能量、体积曲线上， 也是一条连续变化的曲线，但在Tg处应有一个微小的斜率变化[24−25]。由于曲线具有连续性，因此，可以对数据进行整体拟合，得到如下关系式：

图2 在不同温度下Cu66Ti34合金中Cu原子与Ti原子的均方位移曲线Fig.2 MSD curves of Cu and Ti atoms at different temperatures: (a)MSD for Cu at 900−1 500 K; (b)MSD for Cu at 300−800 K;(c)MSD for Ti at 900−1 500 K; (d)MSD for Ti at 300−800 K

图3 Cu66Ti34合金偶关联函数第一谷的最小值与第一峰最大值之比Fig.3 Ratios of minimum to maximum of the first peak of PCF for Cu66Ti34 alloy

图4 Cu66Ti34合金凝固过程中焓随温度变化的关系Fig.4 Enthalpy of Cu66Ti34 vs temperature during solidification process

金属液体处于Fragile和Strong两种极限之间，比较接近Fragile极限，这个系统的特征是在Tg点时比热容有大的突变。计算定压比热容有两种方法[26]：一是利用能量的涨落来计算，二是用热力学求导的方法计算，即cp=dh/dT。由于第一种方法需要的原子数目较多，故采用第二种方法计算。即

当T的取值范围为300～1 600 K时，可得到图5所示曲线。从图5可以看出，函数在800～900 K之间存在最大值。这条曲线虽不能看出 cp的阶跃性变化[25]，但可以估计玻璃化转变的发生。BAILEY等[27]认为函数的最大值点即是 Tg点。对 cp求极值，得到Tg=892 K，这样求出的Tg值比较具体。这与用偶分布函数第一谷的最小值与第一峰最大值之比求得的 Tg点相比，处于温度更高的位置。一方面，BAILEY[27]所考虑的是极值点，而通过热分析获得的cp具有阶跃变化特征，这个点是偏向温度更低的范围。另一方面，是由于 cp的增量是由液态金属平移自由度所做的贡献[28]，这说明玻璃化转变与原子的运动有关，而原子的运动本身具有滞后性，从而导致了玻璃化转变相对于时间和温度具有滞后性，进一步说明Tg只是动力学参量。综合研究从动力学和热力学性质分别得到的玻璃化转变温度可知，用动力学性质计算的Tg与实验结果更接近。

图5 Cu66Ti34合金凝固过程比热容cp随温度的变化Fig.5 Specific heat capacity cp vs temperature during solidification process of Cu66Ti34 alloy

3 结论

1) GEAM势函数能够较好地反映Cu66Ti34合金玻璃化转变的动力学性质，通过偶分布函数第一谷的最小值与第一峰最大值之比获得玻璃化转变温度为 600 K，与相近成分Cu50Ti50的实验值接近。

2) 在模拟条件下根据均方位移数值，可以知道Cu的扩散能力要比Ti的强，且都在600 K开始呈现明显冻结状态。

3) 合金的生成焓随温度的降低而降低，比热容与温度具有二次函数关系。可以利用cP函数图像来估计玻璃化转变的发生。由于Tg的动力学性质，采用动力学方法计算的Tg与实验结果更接近。

[1]LIN X H, JOHNSON W L. Formation of Ti-Zr-Cu-Ni bulk metallic glasses[J]. Journal of Applied Physics, 1995, 78(11):6514−6519.

[2]XU D H, DUAN G, JOHNSON W L. Unusual glass-forming ability of bulk amorphous alloys based on ordinary metal copper[J]. Physical Review Letters, 2004, 92(24):245503−245506.

[3]LEE S W, HUH M Y, FLEURY E, LEE J C.Crystallization-induced plasticity of Cu-Zr containing bulk amorphous alloys[J]. Acta Mater, 2006, 54(2): 349−355.

[4]DAI C L, GUO H, SHEN Y, LI Y, MA E, XU J. A new centimeter-diameter Cu-based bulk metallic glass[J]. Scripta Materialia, 2006, 54: 1403−1408.

[5]SHEN Y, MA E, XU J. A group of Cu(Zr)-based BMGs with critical diameter in the range of 12 to 18 mm[J]. Journal of Materials Science and Technology, 2008, 24(2): 149−152.

[6]JIA P, GUO H, Li Y, XU J, MA E. A new Cu-Hf-Al ternary bulk metallic glass with high glass forming ability and ductility[J].Scripta Materialia, 2006, 54: 2165−2168.

[7]TANG M B, ZHAO D Q, PAN M X, WANG W H. Binary Cu-Zr bulk metallic glasses[J]. Chinese Physics Letters, 2004, 21(5):901−903.

[8]XU D H, LOHWONGWATANA B, DUAN G, JOHNSON W L,GARLAND C. Bulk metallic glass formation in binary Cu-rich alloy series-Cu100−x(x = 34, 36, 38.2, 40 at.% )and mechanical properties of bulk Cu64Zr36glass [J]. Acta Materialia, 2004, 52(9):2621−2624.

[9]WANG D, LI Y, SUN B B, SUI M L, LU K, MA E. Bulk metallic glass formation in the binary Cu-Zr system[J]. Applied Physics Letters, 2004, 84 (20): 4029−4031.

[10]SATO S, TAKAHASHI Y, SANADA T, WAGATSUMA K,SUZUKI S. Small-angle X-ray scattering characterization of precipitates in Cu-Ti alloys[J]. Journal of Alloys and Compounds,2009, 447(1/2): 846−850.

[11]GHOSH G. First-principles calculations of structural energetics of Cu-TM (TM =Ti, Zr, Hf)intermetallics[J]. Acta Materialia,2007, 55: 3347−3374.

[12]OGATA S, SHIMIZU F, LI J, WAKEDA M, SHIBUTANI Y.Atomistic simulation of shear localization in Cu-Zr bulk metallic glass[J]. Intermetallics, 2006, 14(8/9): 1033−1037.

[13]LO Y C, CHOU H S, CHENG Y T, HUANG J C, MORRIS J R,LIAW P K. Structural relaxation and self-repair behavior in nano-scaled Zr-Cu metallic glass under cyclic loading:Molecular dynamics simulation[J]. Intermetallics, 2010, 18(5):954−960.

[14]WADLEY H N G, ZHOU X W, JOHNSON R A, NEUROCK M.Mechanisms, model and methods of vapor deposition[J].Progress in Material Science, 2001, 46(5): 329−377.

[15]ZHOU X W, WADLEY H N G, JOHNSON R A, LARSON D J,TABAT N, CERZO A, PETFORD-LONG A K, SMITH G D W,CLIFTON P H, MARTENS R L, KELLY T F. Atomic scale structure of sputtered metal multilayers[J]. Acta Materialia, 2001,49(19): 4005−4015.

[16]JOHNSON R A. Analytic nearest-neighbor model for fcc metals[J]. Physical Review B, 1988, 37(8): 3924−3931.

[17]WASEDA Y. The structure of non-crystalline materials[M]. New York: Mcgraw-Hill, 1980: 292−300.

[18]梁永超, 刘让苏, 朱轩民, 周丽丽, 田泽安, 刘全慧. 液态Mg7Zn3合金快速凝固过程中微观结构演变机理的研究[J]. 物理学报, 2010, 59(11): 7930−7940.LIANG Yong-chao, LIU Rang-su, ZHU Xuan-min, ZHOU Li-li,TIAN Ze-an, LIU Quan-hui. Simulation study of evolution mechanisms of microstructures during rapid solidification of liquid Mg7Zn3alloy[J]. Acta Physica Sinica, 2010, 59(11):7930−7940.

[19]朱志雄, 张 鸿, 刘超峰, 齐卫宏, 易丹青, 李志成. Ni-Al合金凝固过程的分子动力学模拟[J]. 中国有色金属学报, 2009,19(8): 1409−1416.ZHU Zhi-xiong, ZHANG Hong, LIU Chao-feng, QI Wei-hong,YI Dang-qing, LI Zhi-cheng. Molecular dynamics simulation for solidification process of Ni-Al alloys[J]. The Chinese Journal of Nonferrous Metals, 2009, 19(8): 1409−1416.

[20]WENDT H R, ABRAHAM F F. Empirical criterion for the glass transition region based on Monte Carlo simulations[J]. Physical Review Letters, 1978, 41(18): 1244−1246.

[21]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]. Journal of Chemical Physics, 1980, 72(1): 359−365.

[22]NEUDECKER M, MAYR S G. Dynamics of shear localization and stress relaxation in amorphous Cu50Ti50[J]. Acta Materialia,2009, 57: 1437−1441.

[23]刘 俊, 赵九洲, 胡壮麒. Cu-50%Ni合金快速凝固过程中原子团簇演变的分子动力学模拟[J]. 金属学报, 2005, 41(2):219−224.LIU Jun, ZHAO Jiu-zhou, HU Zhuang-qi. Molecular dynamic simulation of atomic clusters evolution Cu-50%Ni in alloy during rapid solidification[J]. Acta Metallurgica Sinica, 2005,41(2): 219−224.

[24]CAHN R W, HAASEN P. Physical metallurgy[M]. North Holland: Elsevier Science B V, 1996.

[25]ÖZDEMIRKART S, TOMAK M, ULUDOĞAN M, ÇAĞIN T.Structural, thermodynamical, and transport properties of undercooled binary Pd-Ni alloys[J]. Materials Science and Engineering A, 2006, 435/436: 736−744.

[26]杨 弘, 陈 民. 深过冷液态Ni2TiAl合金热物理性质的分子动力学模拟[J]. 物理学报, 2006, 55(5): 2418−2421.YANG Hong, CHEN Min. A molecular dynamics simulation of thermodynamic properties of undercooled liquid Ni2TiAl alloy[J].Acta Physica Sinica, 2006, 55(5): 2418−2421.

[27]BAILEY N P, SCHIØTZ J, JACOBSEN K W. Simulation of Cu-Mg metallic glass: Thermodynamics and structure [J].Physical Review B, 2004, 69 (14): 144205−144215.

[28]陈 光, 傅恒志. 非平衡凝固新型金属材料[M]. 北京: 科学出版社, 2004: 99.CHEN Guang, FU Heng-zhi. New metal materials of nonequilibrium solidification [M]. Beijing: Science Press, 2004: 99.