分子动力学模拟Pb₃₀Ag₇₀, Pb₄₀Ag₆₀和Pb₆₀Ag₄₀合金热力学性质的研究

贾国斌 杨斌 刘嘉

昆明理工大学真空冶金国家工程实验室

邢台钢铁有限责任公司

摘 要：

采用普适的嵌入原子模型模拟Pb30Ag70, Pb40Ag60和Pb60Ag403种合金体系的升温过程, 用于计算合金体系的微观及宏观性质, 首先模拟298¹498 K时体系的生成焓, 进而得到不同温度下合金体系生成自由能和过剩自由能, Pb30Ag70, Pb40Ag60和Pb60Ag403种合金的过剩自由能分别为1.26, 1.66和1.90 kJ, 与实验值符合较好, 3种合金的过剩自由能均为正值, 合金中原子间平均相互作用较小。同时计算合金的结合能及形成能等能量函数, 计算得到Pb-Ag合金的结合能随温度升高不断降低, 即合金中原子间平均相互作用不断降低。Pb-Ag合金的形成能为正值, 表明合金为正偏差体系, 随着温度的升高, 合金形成能的绝对值不断降低趋近于0, 可以判断合金与理想熔体间的偏离程度不断降低, 合金体系趋近于理想熔体, 合金中Pb含量越大, 合金与理想熔体的偏离程度越小, 通过对形成能的计算可以定量的判断合金与理想熔体的偏离程度。

关键词：

分子动力学;能量函数;普适的嵌入原子方法;

中图分类号： TG111.3

作者简介：杨斌, 通讯联系人, (E-mail:kgyb2005@126.com) ;

收稿日期：2009-12-25

基金：国家自然科学基金—云南联合基金 (U0837604) 资助项目;

Molecular Dynamics Simulation on Thermodynamic Properties of Pb₃₀Ag₇₀, Pb₄₀Ag₆₀ and Pb₆₀Ag₄₀ Alloy

Abstract：

The heating processes of Pb30Ag70, Pb40Ag60 and Pb60Ag40 alloys were simulated with generalized embedded atom method to calculate microscopic and macroscopic properties.The enthalpy of formation at 298¹498 K was simulated at first, then free energy of formation and excess free energy were obtained.The excess free energies of Pb30Ag70, Pb40Ag60 and Pb60Ag40 were 1.26, 1.66 and 1.90 kJ, respectively, this agreed well with the experimental values.The excess free energies of three alloys were all positive, so the average atomic interaction was slight.The energy functions, such as cohesive energy and formation energy were calculated, and cohesive energy of Pb-Ag alloy decreased with the increase of the temperature, indicating the decrease of average atomic interaction.Formation energy of Pb-Ag alloy was positive, so the alloy showed positive deviation.With the temperature increasing, absolute value of formation energy approached to zero, and the alloy approached to ideal melt.The more the lead content, the smaller the deviation degree.The calculation of formation energy could describe the deviation degree between the actual alloy and the ideal melt quantitatively.

Keyword：

molecular dynamics (MD) ;energy functions;generalized embedded atom method;

Received： 2009-12-25

混合效应使合金与理想熔体的热力学性质发生一定程度的偏离, 从而使合金的分离提纯过程变得复杂, 因此对合金体系热力学性质的研究, 变得尤为重要。 由于高温实验的复杂性, 目前仅有少量体系的热力学性质得到了研究, 难以满足科研工作者对合金热力学性质的需求。 几十年来许多研究者对溶液热力学性质的解析式提出了大量的理论与经验模型, 并取得了很大的进展 ^[1,2,3,4] , 但各种模型都有一定的适用条件和局限性。 随着计算机技术的日益发展和推广普及, 计算机模拟已成为材料和冶金领域中的一个可靠、 有效的研究方法, 其中分子动力学模拟方法由于其程序简单、 计算量小、 应用范围广, 得到了众多研究者的关注。

近年来各种经验和半经验多体势发展很快, 在精确研究材料的热力学和结构方面发挥很大作用。 适用于金属材料的多体势可以分为3类, 即嵌入原子法 (Embedded Atom Method-EAM) ^[5,6] , F-S势 (Finnis-Sinclair Model) ^[7] 和紧束缚TB势 (Tight Binding Formalism) ^[8] 。 由Daw和Baskes提出的EAM模型是应用最广泛的多体势模型, EAM模型是基于密度泛函理论的半经验模型, 认为原子系统的能量可由其电子密度函数精确给出。 EAM势能函数及改进式在研究金属单质及合金的晶体缺陷、 表面性质、 薄膜生长、 吸附性能及液态结构等方面取得较好的结果 ^[9,10,11,12] , 但在合金热力学性质的研究方面所做工作较少, 本文采用分子动力学方法模拟了不同成分下Pb-Ag合金的熔化过程, 计算298～1498 K下Pb-Ag合金的生成焓、 生成自由能、 过剩自由能、 结合能等能量函数, 同时从宏观和微观角度分析合金原子间平均相互作用的大小, 计算合金的形成能用来定量的描述实际合金与理想熔体的偏离程度。

1 分子动力学模拟

分子动力学模拟的关键在于选择准确的势能模型来描述分子间相互作用。 由Zhou提出的普适的嵌入原子模型 (GEAM势) ^[13,14] 可以用于若干金属和合金的分子动力学模拟, 是近年来常用于描述过渡族金属原子间相互作用的势能模型, 其基本方程如下:

E=12∑i,j,i≠jΦij(rij)+∑iFi(ρi)???(1)

式 (1) 中Φ_ij (r_ij) 代表了相距为r_ij的两个原子i, j之间的对势作用, F_i (ρ_i) 表示原子i的嵌入原子能, ρ_i可以表示为 ρi=∑j,j≠ifi(rij)?fi(rij) 表示与原子i相距r_ij的原子j在i原子所在的位置产生的电子密度, 可以表示为:

fi(rij)=fee-β(rre-1)1+(rre-λ)20???(2)

嵌入原子能F_i (ρ_i) 可以表示成:

Fi(ρi)=3∑i=0Fni(ρρn-1)i?ρ<ρnρn=0.85ρeFi(ρi)=3∑i=0Fi(ρρe-1)i?ρn<ρ<ρoρo=1.15ρeFi(ρi)=F0[1-ln(ρρe)η](ρρe)η?ρ0≤ρ???(3)

F_ni, F_i, F_o, η为模型参数, 可由结合能和体模量计算得到, ρ_e为平衡时的电子密度。

两体作用势为:

?(r)=A-αe(rre-1)1+(rre-κ)20-B-βe(rre-1)1+(rre-λ)20???(4)

r_e是最近邻原子间的平衡距离, A, B, α, β为4个模型参数, κ, λ, m, n是有关截断半径的附加参数, 一般情况下m=n=20。 λ=2κ, α/β=1.875

当计算合金时, 合金势采用 ^[15]

?ab(r)=12(fb(r)fa(r)?aa(r)+fa(r)fb(r)?bb(r))???(5)

其中?^ab (r) 为a-b合金的两体势, f^a (r) , f^b (r) 分别是a和b组分的电子密度函数, ?^aa (r) 和?^bb (r) 分别是a和b组分的两体势。

GEAM势函数中的参数见表1和2。

本文分子动力学模拟在容纳500个原子的立方晶胞中进行, 由于涉及的Pb, Ag两种元素均为面心立方结构, 因此模拟初期立方晶胞中Pb, Ag原子以面心立方结构随机排列, 模拟过程中采用NPT (等温等压) 系综, 使用周期性边界条件。 合金体系首先恒温在298 K, 驰豫1×10⁵个时间步长, 然后以8.5×10¹²K·s^-1的速率升温至1498 K, 最后在1498 K时再驰豫1×10⁵个时间步长, 时间步长选择2×10^-15s, 在模拟过程中压强保持0 Pa。

2 结果与讨论

2.1合金生成焓的计算

图1是计算得到3种合金体系不同温度下的生成焓。 从图1中可以看出计算得到Pb₃₀Ag₇₀, Pb₄₀Ag₆₀和Pb₆₀Ag₄₀合金的熔点温度分别为930, 925和830 K, 文献 [ 16] 中3种合金的实验测定的熔点分别为927, 893和815 K, 模拟相对误差分别为0.32%, 3.58%和1.81%, 模拟值与实验值符合较好。

表1GEAM势的势能参数 (两体势部分)

Table 1Potential parameters of GEAM (pair-potential)

Atom 1	Atom 2	re, i	fe, i	αi	βi	Ai/eV	Bi/eV	κi	Bj/eV
Ag	Ag	2.892	1.106	7.945	4.237	0.266	0.386	0.425	0.386
Pb	Pb	3.450	0.648	8.468	4.516	0.135	0.203	0.426	0.203
Ag	Pb	2.892	1.106	7.945	4.237	0.266	0.386	0.425	0.203
		κj	λj	λi	re, j	fe, j	αj	βj	Aj/eV
Ag	Ag	0.425	0.851	0.851	2.892	1.106	7.945	4.237	0.266
Pb	Pb	0.426	0.852	0.852	3.450	0.648	8.468	4.516	0.135
Ag	Pb	0.426	0.852	0.851	3.450	0.648	8.468	4.516	0.135

表2GEAM势的势能参数 (多体势部分)

Table 2Potential parameters of GEAM (many body potential)

Atom	ρe	Fno/eV	Fn1/eV	Fn2/eV	Fn3/eV	F0/eV	F1/eV	F2/eV	F3/eV	η	Fe/eV
Ag	15.539	1.730	-0.221	0.542	-0.967	1.750	0.000	0.984	0.521	1.149	1.751
Pb	8.907	1.420	-0.229	0.630	-0.561	1.440	0.000	0.921	0.109	1.172	1.440

图1 合金生成焓随温度的变化曲线

Fig.1 Temperature dependence of the enthalpy of formation

2.2合金过剩自由能的计算

根据合金体系生成焓, 可以方便地计算体系在升温过程中热力学性质。 合金体系生成自由能与温度的关系可以由公式 (6) 来表示。

ΔG (T) =ΔH (T) -TΔS (T) (6)

由于在熔点位置时, 体系的液态生成自由能等于固态生成自由能 ^[17] , 即: G_C (T_M) =G_L (T_M) , 因此在熔点位置有:

ΔS_LC (T_M) =ΔH (T_M) /T_M (7)

可以根据式 (7) , 计算任意温度下合金体系的生成熵。

ΔS(Τ)=Τ∫0ΚdΤ[?ΤΔΗ(Τ)]/Τ?Τ≤ΤΜΔS(Τ)=S(ΤΜ)+ΔSLC(ΤΜ)+Τ∫ΤΜdΤ[?ΤΔΗ(Τ)]/ΤΤ>ΤΜ???(8)

将 (8) 式代入 (6) 式即可计算出不同温度下合金体系的生成自由能。

在1273 K时, 模拟得到的Pb₃₀Ag₇₀, Pb₄₀Ag₆₀和Pb₆₀Ag₄₀ 3种合金体系的生成自由能分别为: -94.68, -97.36和-104.16 kJ, 文献 [ 16] 中3种合金的实验值为-94.72, -98.57和-105.33 kJ, 模拟相对误差分别为0.04%, 1.23%和1.12%, 可以看出模拟值与实验值符合较好。

二元Pb-Ag合金体系的过剩自由能可以由式 (8) 计算得到。

ΔG^E=ΔG_mix-ΔG idmix =ΔG_mix-x_PbRTlnx_Pb-x_AgRTlnx_AG (9)

其中ΔG^E是合金的过剩自由能, ΔG_mix是合金的混合自由能, ΔG idmix 是合金的理想混合自由能, x_i, x_j为组分i, j的摩尔分数, 合金的混合自由能可以由 (10) 式来表示。

ΔG_mix=ΔG-x_AgΔG_Ag-x_PbΔG_Pb (10)

其中ΔG_Ag, ΔG_Pb为纯Ag, Pb金属不同温度下的生成自由能, x_Ag, x_Pb分别为Ag, Pb组分的摩尔分数, ΔG为合金生成自由能。

图2为3种合金体系不同温度下的过剩自由能, 从图2中可以看出1273 K时计算得到的Pb₃₀Ag₇₀, Pb₄₀Ag₆₀和Pb₆₀Ag₄₀ 3种合金的过剩自由能分别为: 1.26, 1.66和1.90 kJ, 实验值 ^[16] 分别为: 0.92, 1.21和1.58 kJ。

从图2中可以看出3种合金的过剩自由能均为正值, 说明合金中原子间平均相互作用较小, 体系为正偏差体系; 随着温度的升高, 合金的过剩自由能增大, 即升温过程中合金中原子间平均相互作用不断降低, 升温过程有利于合金的分离。

2.3合金结合能的计算

合金体系内原子间相互作用还可以由合金的结合能来判断。 图3为计算得到的不同温度下金属单质及合金结合能。 在298 K时Pb, Ag两种金属单质的计算结合能为1.98, 2.81 eV, 实验值 ^[18] 分别为: 2.04, 2.85 eV, 模拟的相对误差为2.88%和1.45%。

从图3中可以看出:

(1) 相同温度下, 纯铅的结合能小于纯银的结合能, 说明铅原子间的相互作用比银原子间的相互作用小, 铅原子较银原子更容易形成孤立原子, 宏观上表现为铅的蒸气压比银高。

图2 合金过剩自由能随温度下的变化曲线

Fig.2 Temperature dependence of excess free energy

图3 不同温度下合金的结合能

Fig.3 Temperature dependence of the cohesive energy

(2) Pb-Ag合金的结合能随温度升高不断降低, 即合金中原子间平均相互作用不断降低, 说明温度越高, 合金液态表面越易形成孤立的金属原子, 使金属的蒸气压变大, 升温有利于金属的蒸发。

2.4合金形成能的计算

过剩自由能可以判断合金体系为正偏差或负偏差, 但不能表明合金体系与理想熔体的偏离程度, 从图2中可以看出, 合金体系的过剩自由能随温度的升高不断增大, 但不能说明合金与理想熔体的偏离程度不断增大, 相反温度越高合金体系与理想溶液的偏离程度越小。 合金的形成能可以从微观角度来分析合金熔体与理想熔体的偏离情况, 合金形成可以表示为:

E_f=-[E_{C (Pb-Ag)} - (n_PbE_{C (Pb)} +n_AgE_{C (Ag)} ) / (n_Pb+n_Ag) ] (11)

式 (11) 中n_Pb, n_Ag分别表示合金中铅、 银组分原子个数, E_{C (Pb-Ag)} , E_{c (Pb)} , E_{C (Ag)} 分别表示Pb-Ag合金体系、 纯铅、 纯银体系的结合能。 当合金体系的形成能为0时, 合金为理想熔体; 当合金体系的形成能为正时, 说明合金中原子间平均相互作用较小, 体系为正偏差; 当形成能为负值时, 合金中原子间平均相互作用较大, 体系为负偏差。 合金形成能的绝对值越大, 表明实际合金与理想熔体的偏离程度越大。 图4为合金体系的形成能随温度的变化曲线。

从图4中可以看出Pb-Ag合金体系的形成能为正值。 表明此温度下3种合金均为正偏差体系, 随着温度的升高, 合金形成能的绝对值不断降低趋近于0, 可以判断合金与理想熔体间的偏离程度不断降低, 合金体系不断趋近与理想熔体, 3种合金中Pb含量越大, 合金与理想熔体的偏离程度越小, 以上结果与文献 [ 16] 相一致, 即合金的形成能可以描述不同温度不同成分合金与理想熔体间的偏离程度。

图4 合金形成能随温度的变化曲线

Fig.4 Temperature dependence of formation energy

3 结 论

1. 计算了3种合金的能量函数, 可以看出计算值与实验值吻合, 说明GEAM势适用于Pb-Ag合金体系热力学性质的模拟。

2. 3种合金的过剩自由能均为正值, 说明合金中原子间平均相互作用较小, 体系为正偏差体系。 随着温度的升高, 合金的过剩自由能增大, 说明升温过程中合金中原子之间平均相互作用不断降低, 升温过程有利于合金的分离。

3. 体系的结合能随温度的升高不断降低, 表明合金中原子间平均相互作用不断降低, 因此升温过程有利于合金的分离。

4. 计算得到的合金形成能可以定量的描述实际合金与理想熔体的偏离程度。

参考文献

[1] Le Qichi, Zhang Xinjian, Cui Jianzhong, Lu Guimin.Progres-sion of research on thermodynamic model of metallic alloy melts[J].Acta Metallurgica Sinica, 2003, 39:35. (乐启炽, 张新建, 崔健忠, 路贵民.金属合金溶液热力学模型研究进展[J].金属学报, 2003, 39:35.)

[2] Cao Zhanmin, Song Xiaoyan, Qiao Zhiyu.Thermodynamicmodeling software FactSage and its application[J].ChineseJournal of Rare Metals, 2008, 32 (1) :216. (曹战民, 宋晓艳, 乔芝郁.热力学模拟计算软件FactSage及其应用[J].稀有金属, 2008, 32 (1) :216.)

[3] Chou K C.Anewsolution model for predicting ternary thermody-namic properties[J].Calphad, 1987, 11:293.

[4] Tao D P.Prediction of the thermodynamic properties of binarycontinuous solid solutions by infinite dilute activity coefficients[J].Thermochimica Acta, 2003, 408:67.

[5] DawMS, Baskes MI.Embedded atommethod:Derivation andapplication to impurities, surfaces, and other defects in metals[J].Phys.Rev.B, 1984, 29:6443.

[6] Foiles M S, BaskesS M I, Daw M S.Embedded-atom-methodfunctions for the fcc metals Cu, Ag, Au, Ni, Pd, Pt and their al-loys[J].Phys.Rev.B, 1986, 33:7983.

[7] Finnis M W, Sinclair J E.A simple empirical N-body potentialfor transition metals[J].Phil.Mag., 1984, 50:45.

[8] Tomanek D, Mukherjee S, Bennemann K H.Simple theory forthe electronic and atomic structure of small clusters[J].Phys.Rev.B, 1983, 28:665.

[9] Wang Li, Cong Rihong, Zhang Ningyan, Bian Fangxiu.Medi-um-range order of liquid metal in the quenched state[J].Physi-ca B, 2005, 355:140.

[10] Caro A, Turchi P E A, Caro M, Lopasso E M.Thermodynam-ics of Fe-Cu alloys as described by a classic potential[J].Jour-nal of Nuclear Materials, 2006, 349:317.

[11] Qi L, Zhang H F, Hu Z Q.Molecular dynamic simulation ofglass formation in binary liquid metal:Cu-Ag using EAM[J].Intermetallics, 2004, 12:1191.

[12] Kazanc S.Molecular dynamics study of pressure effect on crystal-lization behaviour of amorphous CuNi alloy during isothermal an-nealing[J].Physics Letters A, 2007, 365:473.

[13] Zhou X W, Wadley HNG, Johnson R A, Larson D J, Tabat N, Cerezo O A, Petford-Long A K, Smith G D W, Clifton P H, Martens R L, Kelly T F.Atomic scale structure of sputteredmetal multilayers[J].Acta Mater., 2001, 49:4005.

[14] Wadley HNG, Zhou X W, Johnson R A, Neurock M.Mecha-nisms, model and methods of vapor deposition[J].Progress inMaterial Science, 2001, 46:329.

[15] Johnson R A.Analytic nearest-neighbor modle for fcc metals[J].Phys.Rev.B, 1988, 37:3924.

[16] Hultgren R, Dessai P D, Hawkins D T, Gleiser M, Kelley K K.Selected Values of the Thermodynamic Properties of Binary Alloy[M].OHIO:ASMMetal Park, 1973.

[17] Teicher H.Melting transition in molecular-dynamics simulationsof the Ni0.5Zr0.5intermetallic compound[J].Phys.Rev.B, 1999, 59:8473.

[18] Smith C J.Metals Reference Book[M].London:Butterworrd, 1976.