文章编号：1004-0609(2011)09-2151-06

Cu₆₆Ti₃₄非晶合金凝固过程的分子动力学模拟

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

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

摘  要：采用分子动力学模拟了二元合金Cu₆₆Ti₃₄的凝固过程。原子间作用采用GEAM势，利用偶关联函数，均方位移(MSD)等分析方法，研究Cu₆₆Ti₃₄合金在4×10¹³ K/s冷却速度下的玻璃化转变温度、原子的扩散行为。结果表明，通过偶分布函数第一谷的最小值与第一峰最大值之比获得的玻璃转变温度为600 K，与相近成分Cu₅₀Ti₅₀的实验值接近；在800 K时，Cu和Ti的MSD最大值均小于1×10^-2 nm²，合金熔体很粘稠；在600 K时，曲线的斜率降低，在动力学上合金熔体已经凝固。定压比热容与温度成二次分布关系，存在一个峰值温度为892 K的热力学玻璃转变温度，证明了用动力学方法和用热力学方法获得的玻璃转变温度之间的差异。

关键词：Cu₆₆Ti₃₄合金；分子动力学模拟；玻璃转变温度；均方位移；热力学分析

中图分类号：O561.1　　     文献标志码：A

Solidification process of Cu₆₆Ti₃₄ amorphous alloy

 simulated by molecular dynamics

QIU Ke-qiang, LI Hong, SUN Jing, YOU Jun-hua, REN Ying-lei, LI Qing-feng

 (School of Materials Science and Engineering, Shenyang University of Technology, Shenyang 110870, China)

Abstract: The solidification process of Cu₆₆Ti₃₄ binary alloy was studied by molecular dynamics simulation. The general embedded atomic (GEAM) potential was selected as many-body interaction. The analytical methods, such as pair correction function (PCF) and mean square displacement (MSD), were employed to calculate the glass transition temperature (T_g) and the diffusion behaviors of atoms at cooling rate of 4×10¹³ K/s. The results show that the glass transition temperature obtained from inflection point of the curve is 600 K, which is the ratio of minimum to maximum of the first peak of PCF for Cu₆₆Ti₃₄ alloy, and is comparable to that of Cu₅₀Ti₅₀ alloy. The maximum of MSD is lower than 1×10^-2 nm² when the temperature is 800 K, showing that the melt is very sticky. However, when the temperature reaches 600 K, the slope of the curve becomes smaller, showing that the melt is frozen dynamically, which shows that the glass transition temperature varies possibly in certain range due to dynamical factor. Simultaneously, the calculated specific heat capacity (c_p) at constant pressure shows that it distributes in secondary order vs temperature. The peak temperature or the glass transition temperature is 892 K, which shows the difference between the glass transition temperatures obtained by dynamics and thermodynamics, respectively.

Key words: Cu₆₆Ti₃₄ alloy; molecular dynamic simulation; glass transition temperature; mean square displacement; thermodynamical analysis

块体Cu基非晶合金自从被报道以来^[1]，因其具有优异的力学性能和相对较低的材料成本而受到关注。目前通过成分设计可以获得多种Cu基块体非晶合  金^[2-6]，其中Cu₅₀Zr₅₀二元合金可以形成直径为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势的表达形式如下：

                    (1)

式中：代表间距为的两个原子i和j之间的对势作用；F_i(r_i)表示原子i的嵌入原子能；r_i可以表

示为，表示与原子i相距r_ij的原子j

在原子i所在的位置产生的电子密度，则嵌入原子能F_i(r_i)可以表示成：

    (2)

式中：F_ni、F_i、F_e和n为模型参数，可由结合能和体模量计算得到；r_e为平衡时的电子密度。

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

                          (3a)

                    (3b)

                    (3c)

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

               (4)

式中：为合金的两体势；和分别为a和b组分的两体势；和分别为a和b组分的电子密度函数。Cu₆₆Ti₃₄二元合金的势参数如表1所列。

1.2  模拟方法

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

表1  Cu₆₆Ti₃₄ 二元合金的势参数

Table 1  Potential parameters of binary Cu₆₆Ti₃₄ alloy

2  计算结果及讨论

2.1 偶关联函数

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

           (5)

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

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

2.2 均方位移(MSD)

粒子位移平方的平均值称为均方位移，均方位移随时间的变化表征了液态金属粒子的扩散行为。为了便于观察Cu₆₆Ti₃₄合金中Cu和Ti原子在急冷到各个温度下MSD的图像变化，以900 K为分界，将图像分开展示，图2(a)和(c)所示分别为Cu₆₆Ti₃₄合金中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^-2 nm²，合金熔体已经很粘稠。600 K时，在0.25 ps后，曲线的斜率降低，此时可认为扩散受阻，同时，也是动力学明显呈现冻结状态的起始温度。到500 K时，原子扩散能力很弱，但局域结构仍可以表现出随时间的波动性。到300 K时，原子几乎没有扩散，说明原子已经找到它们的平衡位置。通过图2(a)与(c)，图2(b)与(d)中MSD的数值比较，还可以看出Cu的扩散能力比Ti的强。

图1  Cu₆₆Ti₃₄二元合金不同温度下偶关联函数曲线

Fig.1  Pair correction function curves of Cu₆₆Ti₃₄ binary alloy at different temperatures: (a) g_Cu-Cu; (b) g_Ti-Ti; (c) g_Cu-Ti; (d) g_tot

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

2.3 热力学性质分析

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

h(T)=-0.000 010 214 T ³+0.027 335 3 T ²+

7.130 82 T -392 102 (J·mol^-1)

图2  在不同温度下Cu₆₆Ti₃₄合金中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  Cu₆₆Ti₃₄合金偶关联函数第一谷的最小值与第一峰最大值之比

Fig.3  Ratios of minimum to maximum of the first peak of PCF for Cu₆₆Ti₃₄ alloy

图4  Cu₆₆Ti₃₄合金凝固过程中焓随温度变化的关系

Fig.4  Enthalpy of Cu₆₆Ti₃₄ vs temperature during solidification process

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

c_p=dh/dT=-0.000 030 642 T ²+0.054 671 T+

7.130 82 (J·mol^-1·K^-1)

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

图5  Cu₆₆Ti₃₄合金凝固过程比热容c_p随温度的变化

Fig.5  Specific heat capacity c_p vs temperature during solidification process of Cu₆₆Ti₃₄ alloy

3  结论

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

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

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

REFERENCES

[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-Cu_100-x (x = 34, 36, 38.2, 40 at.% ) and mechanical properties of bulk Cu₆₄Zr₃₆ glass [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] 梁永超, 刘让苏, 朱轩民, 周丽丽, 田泽安, 刘全慧. 液态Mg₇Zn₃合金快速凝固过程中微观结构演变机理的研究[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 Mg₇Zn₃ alloy[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 Cu₅₀Ti₅₀[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] 杨 弘, 陈 民. 深过冷液态Ni₂TiAl合金热物理性质的分子动力学模拟[J]. 物理学报, 2006, 55(5): 2418-2421.
YANG Hong, CHEN Min. A molecular dynamics simulation of thermodynamic properties of undercooled liquid Ni₂TiAl 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 non- equilibrium solidification [M]. Beijing: Science Press, 2004: 99.

(编辑 何学锋)

基金项目：沈阳市科技局计划资助项目(1091177-1-00)

收稿日期：2010-09-13；修订日期：2011-03-11

通信作者：邱克强，教授，博士；电话：024-25499927；传真：024-25499928；E-mail：kqqiu@yahoo.com.cn