金属学报, 2020, 56(1): 112-118 DOI: 10.11900/0412.1961.2019.00257

研究论文

立方晶体弹性常数和EAM/FS势函数的关系

段灵杰, 刘永长,

天津大学材料科学与工程学院水利安全与仿真国家重点实验室  天津 300354

Relationships Between Elastic Constants and EAM/FS Potential Functions for Cubic Crystals

DUAN Lingjie, LIU Yongchang,

State Key Lab of Hydraulic Engineering Simulation and Safety, School of Materials Science and Engineering, Tianjin University, Tianjin 300354, China

通讯作者: 刘永长,ycliu@tju.edu.cn,主要从事金属结构材料研究

责任编辑: 肖素红

收稿日期: 2019-08-01   修回日期: 2019-09-11   网络出版日期: 2019-12-25

基金资助: 国家自然科学基金项目.  U1660201
国家磁约束核聚变能源研究专项课题项目.  2015GB119001

Corresponding authors: LIU Yongchang, professor, Tel: (022)85356410, E-mail:ycliu@tju.edu.cn

Received: 2019-08-01   Revised: 2019-09-11   Online: 2019-12-25

Fund supported: National Natural Science Foundation of China.  U1660201
National Magnetic Confinement Fusion Energy Research Project.  2015GB119001

作者简介 About authors

段灵杰,男,1992年生,博士生

摘要

为了确定bcc和fcc晶体的EAM/FS势函数参数,研究了EAM/FS势函数与弹性常数的关系。对bcc和fcc结构晶体,分别导出压强(P)和体积弹性模量(B)、弹性常数(C44)和剪切弹性模量(Cp=(C11-C12)/2)的利用嵌入函数、对势函数和电子密度分布函数的表达式。发现C44Cp的大小不仅取决于所考虑的原子和周围原子的距离,还受近邻原子排列状况的影响。将关于结合能(utot)和PBC44Cp的5个拟合方程转化为求最小值的优化模型,给出了5种典型bcc结构晶体(V、Mo、Nb、Ta、W)和3种典型fcc结构晶体(Cu、γ-Fe、Ni)的结合能中待定参数的值。采用上述参数计算得到的最小结合能与实验结果一致,此时所对应的原子间距与晶格常数相同,表明了方法的有效性。

关键词: 立方晶体 ; 弹性常数 ; EAM势 ; FS势

Abstract

Potential functions are extensively applied in molecular dynamics (MD) simulation of metals. Selection of them is a very important step in MD simulations due to its effects of the precision and reliability of the simulations. They are one of the most important reference data during the process of calculation. In order to cover the shortage of pairwise potentials for modelling transition metals, EAM/FS many-body potentials have been introduced since 80's of last century. For the sake of determining parameters in the EAM/FS potential functions of bcc and fcc crystals through macro mechanical properties, relations between the EAM/FS potential functions and elastic constants were investigated in this work. Expressions of the pressure (P) and the bulk modulus (B), elastic constant (C44) and shear elastic modulus (Cp=(C11-C12)/2) in terms of the embedding function, pair potential function and the electron density distribution function were deduced for bcc and fcc structures, respectively. It was found that the magnitude of the C44 and Cp depends on the distances between the considered atom and surrounding atoms, but also the configuration of surrounding atoms. Finally, by converting five fitting equations about the cohesive energy (utot) and P, B, C44, Cp into an optimization model of finding minimum value, the values of the six undetermined parameters in the cohesive energy were given for five typical bcc crystals (V, Mo, Nb, Ta and W) and three typical fcc crystals (Cu, γ-Fe, Ni), respectively. For each crystal, calculation errors show accuracy of parameter values. The obtained calculation results, for the minimum cohesive energy and the corresponding atomic distance, fit well with the reported experimental data, by adopting the above values of the parameters, which indicates the effectiveness for our method.

Keywords: cubic crystal ; elastic constant ; EAM potential ; FS potential

PDF (1397KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

段灵杰, 刘永长. 立方晶体弹性常数和EAM/FS势函数的关系. 金属学报[J], 2020, 56(1): 112-118 DOI:10.11900/0412.1961.2019.00257

DUAN Lingjie, LIU Yongchang. Relationships Between Elastic Constants and EAM/FS Potential Functions for Cubic Crystals. Acta Metallurgica Sinica[J], 2020, 56(1): 112-118 DOI:10.11900/0412.1961.2019.00257

势函数被广泛应用于金属体系的分子动力学模拟,是计算过程中重要的基本参数之一[1,2,3,4,5,6]。为了弥补对势在过渡金属建模时的不足,人们提出了一些多体势模型。Daw和Baskes[7]基于经验和量子力学模型提出了嵌入原子方法 (embedded atom method,EAM)。同时,Finnis和Sinclair[8]根据密度泛函理论提出了Finnis-Sinclair势(FS势)。这些模型得到了广泛的应用,Zhang等[9,10]将EAM用于原子尺度的材料设计理论,并提出了修正的EAM,用以研究bcc结构过渡金属。Johnson[11,12]研究了EAM参数和空位形成能的关系,并将EAM用于合金的研究。Wang等[13,14]应用EAM计算出液态镍基718合金的密度数据,并用分子动力学模拟了单相固溶合金的热膨胀和相变。Zou等[15]用修正的EAM研究了Ti-Ni合金密度的温度依赖性。EAM模型常用于计算bcc结构过渡金属二元合金的热、动力学性质[16],也有学者用FS势研究了Ni62Nb38合金的结构[17]和过冷Ni-Zr合金的热物理性质[18]

EAM势和FS势基于不同的假设,但在数学上有相同的形式[19,20,21,22,23,24]。通常,EAM/FS势由2部分组成:一部分是位于晶格点阵上的原子之间的相互作用对势;另一部分是镶嵌原子在电子云中所引起的多体相互作用的嵌入势能。根据经验量子论,先对对势函数、嵌入函数和电子密度分布函数进行基本假设,而具体参数则需要实验数据拟合来确定。要通过拟合宏观弹性模量的实验结果来确定EAM/FS势函数,必须先确定弹性模量与EAM/FS势的数学关系。

对于立方晶体,有3个独立的弹性常数C11C12C44。广义Hooke定律可写成[25]

σ11=C11ε11+C12ε22+C12ε33σ22=C12ε11+C11ε22+C12ε33σ33=C12ε11+C12ε22+C11ε33σ23=2C44ε23                               σ13=2C44ε13                              σ12=2C44ε12                              

式中,σij为应力分量,εij为应变分量(i, j=1, 2, 3)。C11C12可以通过体积弹性模量(B)和剪切弹性模量(Cp)表示:

B=13(C11+2C12),     Cp=12(C11-C12)

BCp是便于计算的[26,27]。对于C44Cp,考虑应变能密度W=12i,jσijεij,由广义Hook定律,有如下表达式:

W=12C11(ε112+ε222+ε332)+C12(ε11ε22+ 

ε11ε33+ε22ε33)+2C44(ε232+ε132+ε122)

应用应变张量ε=0ε120ε1200000W可简化为:W=2C44ε122。因此:

C44=142Wε122

应用ε=ε11000-ε110000,则W可简化为:W=(C11-C12)ε112。故有:

Cp=142Wε112

本工作先推导出弹性模量和EAM/FS势函数之间的数学关系,即确定了C44Cp与EAM/FS势的数学表达式,随后确定了bcc和fcc结构金属中与压强(P)、BC44Cp有关的嵌入函数以及对势函数和电子密度分布函数的表达式,最后分析讨论了势参数的拟合中应注意的问题,并通过拟合5种典型bcc结构金属和3种典型fcc结构金属的实验数据,确定出它们的结合能参数。

已往文献[8,20,28,29]中只列出势参数的最后拟合结果,很少具体给出势函数与弹性模量之间的数学表达式,导致势函数的可靠性很难验证,也缺乏普适性。

1 弹性常数和EAM/FS

考虑bcc或fcc结构的理想晶体,其中所有原子是等价的。因此,一个原子的结合能(utot)可写成:

utot=uP+uM                  uP=12V(ri)                uM=fρ,  ρ=ϕri

式中,uP是对势的贡献,uM是多体势的贡献,Σ表示对原点原子所有周围的相关原子求和,ri是所考虑的原点处的原子与周围的第i个原子间的距离,V(ri)是描述排斥的核-核作用的对势函数,ϕ(ri)是电子密度分布函数,f(ρ)是嵌入函数,ρ表示原点原子周围的相关原子在原点处产生的电子云密度之和。

a表示晶格常数。对于bcc结构晶体,第一层有8个原子,它们与原点处的原子有相同的距离ri=3a2 (i=1, 2, , 8);第二层有6个原子,与原点原子的距离是ri=a (i=9, 10, , 14)同理,对于fcc结构晶体,第一层有12个原子,ri=2a2 (i=1, 2, , 12);第二层有6个原子,ri=a (i=13, 14, , 18)

P定义为:

P=-dutotdΩ

B可写成:

B=-ΩdPdΩ

式中,Ω=a3/k是原子体积,对于bcc结构晶体,k=2;对于fcc结构晶体,k=4。

由式(6)得到W为:

W=utotΩ=WP+WM

这里引入了符号:

WP=k2a3Vri,   WM=ka3f(ρ)

由式(4)和(5)可导出C44Cp的表达式。为导出C44,引入位移:

δr=0ε120ε1200000xyz

这里xyz表示空间位置,ε12是小量,故有[30]

xε12=y,    yε12=x,    zε12=0

对于距离r=x2+y2+z2,偏导数可写成:

rε12=rxxε12+ryyε12=2xyr

对于对势项,计算得:

ε12V(ri)=V'(ri)2xiyiri
2ε122Vri=riVri-V'ri4xi2yi2ri3+
2V'rixi2+yi2ri

这里(xi, yi)是第i个原子的坐标。对于多体项,得到:

ε12f(ρ)=f'(ρ)ε12ϕ(ri)
2ε122f(ρ)=f(ρ)ε12ϕ(ri)2+
 f'(ρ)2ε122ϕ(ri)

用类似于式(14)和(15)的方法可得:

2ε122fρ=fρϕ'ri2xiyiri2+ 
 f'(ρ)riϕ(ri)-ϕ'(ri)4xi2yi2ri3+
2ϕ'(ri)xi2+yi2ri                   

利用式(4)中的C44,可导出对势和多体势的贡献(C44PC44M)分别为:

C44P=142WPε122=k8a3riV(ri)-V'(ri)4xi2yi2ri3+
2V'(ri)xi2+yi2ri                                        
C44M=142WMε122=k4a3f(ρ)ϕ'(ri)2xiyiri2+
f'(ρ)riϕ(ri)-ϕ'(ri)4xi2yi2ri3+
2ϕ'(ri)xi2+yi2ri                                  

这样C44为:

C44=C44P+C44M

对于Cp

δr=ε11000-ε110000xyz

这里ε11是小量,可得:

xε11=x,  yε11=-y,  zε11=0

以及

rε11=x2-y2r

类似于C44的导出,对势和多体势对Cp的贡献(CpPCpM)分别为:

CpP=142WPε112=k8a3riV(ri)-V'(ri)
(xi2-yi2)2ri3+2V'(ri)xi2+yi2ri

以及

CpM=142WMε112=k4a3f(ρ)ϕ'(ri)xi2-yi2ri2+
 f'(ρ)riϕ(ri)-ϕ'(ri)(xi2-yi2)2ri3+
2ϕ'(ri)xi2+yi2ri                                        

Cp为:

Cp=CpP+CpM

可见,CpC44依赖于周围原子的坐标。

2 bccfcc结构晶体的弹性模量

在前文讨论的基础上,分别对bcc和fcc结构晶体考虑周围的两层原子来推导弹性模量的表达式。

2.1 bcc结构晶体(k=2)

对于bcc结构晶体,第一层的8个原子有坐标:±a2, a2, a2±a2, -a2, a2±a2, a2, -a2±a2, -a2, -a2,与原点原子的距离为ri=3a2;第二层中的6个原子有坐标:(±a, 0, 0),(0, ±a, 0),(0, 0, ±a),与原点原子的距离为ri=aΩ=a3/2。在式(6)中,此时ρ=8ϕ(3 a2)+6ϕ(a)utot=4V3 a2+3Va+fρ

式(7)和(8)可写成:

P=-23a223V'(3a2)+3V'(a)+ 
f'(ρ)43ϕ'(3a2)+6ϕ'(a)
B=-29a243V'(3a2)+6V'(a)-3aV(3a2)-
-3aV(a)+f'(ρ)83ϕ'(3a2)+
12ϕ'(a)-6aϕ(3a2)-6aϕ(a)-
4af(ρ)23ϕ'(3a2)+3ϕ'(a)2

从式(19)~(21)以及(25)~(27),计算C44Cp,得:

C44=29a243V'(3 a2)+3aV(3 a2)+9V'(a)+ 
 49a243ϕ'(3a2)+3aϕ(3a2)+9ϕ'(a)f'(ρ)
Cp=13a243V'(3 a2)+3V'(a)+3aV(a)+ 
 23a243 ϕ'(3 a2)+3ϕ'(a)+3aϕ(a)f'(ρ)

2.2 fcc结构晶体(k=4)

对于fcc结构晶体,第一层中的12个原子有坐标:±a2, a2, 0±a2, -a2, 0±a2, 0 ,a2±a2, 0, -a20, ±a2, a20, ±a2, -a2,与原点的距离为ri=2a2;第二层中的6个原子有坐标:(±a, 0, 0),(0, ±a, 0),(0, 0, ±a),与原点的距离为ri=aΩ=a3/4。在式(6)中,此时ρ=12ϕ2a2+6ϕautot=6V(2a2)+3V(a)+fρ

由式(7)、(8)、(21)和(27)得到:

 P=-4a22 V'(2 a2)+V'(a)+2f'(ρ)          2 ϕ'(2 a2)+ϕ'(a)
B=-43a222V'2 a2+2V'a-aV2 a2-
aVa+4f'ρ2ϕ'2 a2+ϕ'a-
2af'(ρ)ϕ(2 a2)+ϕ(a)-12af(ρ)
2ϕ'(2a2)+ϕ'(a)2                           
C44=1a232V'(2 a2)+aV(2 a2)+4V'(a)+  2a232 ϕ'(2 a2)+aϕ(2 a2)+4ϕ'(a)f'(ρ)
Cp=12a272V'(2 a2)+aV(2 a2)+4V'(a)+4aV(a)+1a272ϕ'(2a2)+
 aϕ(2 a2)+4ϕ'(a)+4aϕ(a)f'(ρ)

3 势参数的数值拟合

为完成数值拟合,首先需要确定嵌入函数、对势函数和电子密度分布函数的形式,按照紧束缚理论的经典模型,取f(ρ)=-Aρ,这里A>0是待定常数。假设对势函数和电子密度分布函数为经验模型:

V(r)=(c-r)3(c0+c1r+c2r2)(rc)0(r>c)
ϕ(r)=(d-r)3(rd)0(r>d)

式中,常数cd分别表示位于第二层和第三层之间的待定的与原点原子的距离,c0c1c2为待定系数。对于bcc结构晶体,a<c<2aa<d<2a;而对于fcc结构晶体,a<c<1.5aa<d<1.5a。这样假设的对势函数和电子密度分布函数,在分界点cd处都有连续的二阶导数。

6个待定常数Acdc0c1c2的值,需要通过拟合utotBC44Cp的实验结果以及利用平衡条件P=0来确定。对bcc结构,在A>0,a<c<2aa<d<2a范围内;对fcc结构,在A>0,a<c<1.5aa<d<1.5a范围内,确定6个常数,优化下式以使拟合误差最小:

ER=(utot-ûtot)2+P2+(B-B̂)2+
(C44-Ĉ44)2+(Cp-Ĉp)2

式中,ûtotB̂Ĉ44Ĉp分别表示utotBC44Cp的实验值,ER表示拟合误差。这里仅考虑了5种典型bcc结构晶体和3种典型fcc结构晶体,利用符号计算软件MATHEMATICA解决优化问题。这些晶体的实验数据和势参数的拟合值列在表1,2,3,4中,其中在表2和4的最后一列给出了较为理想的拟合误差。表1[8]中除W之外的4种晶体,随晶格常数增大,结合能绝对值也增大。而表3[17,31]中的3种fcc晶体,随晶格常数增大,结合能绝对值减小。从表2和4的拟合误差可知,对势函数和电子密度分布函数模型(36)和(37)更适合bcc结构晶体。

表1   5种典型bcc结构晶体的实验数据(弹性常数的量纲为1011 Pa)[8]

Table 1  Experimental data for the five typical bcc structure crystals (elastic constants with dimension 1011 Pa)[8]

Crystala / nmutot / eVBC44Cp
V0.30399-5.311.5510.4260.546
Mo0.31472-6.822.6261.0891.516
Nb0.33008-7.571.7100.2810.567
Ta0.33058-8.101.9610.8240.524
W0.31652-8.903.1041.6061.590

Note: a—crystal lattice constant, utot—cohesive energy, B—bulk modulus, C44—elastic constant, Cp—shear elastic modulus

新窗口打开| 下载CSV


表2   5种典型bcc结构晶体的势参数的拟合值

Table 2  Fitted values of potential parameters for the five typical bcc structure crystals

CrystalcdAc0c1c2ER
V3.16283.87991.711018.83-16.253.9725.125×10-11
Mo3.23174.36251.3678262.85-200.7038.8809.408×10-10
Nb3.39274.09492.590926.30-22.685.4051.103×10-10
Ta3.48634.29682.049825.13-17.983.5849.312×10-12
W3.26014.47631.5395275.77-207.6739.6552.896×10-3

Note: A—constant in the embedding function; d—constant in the electron density distribution function; c, c0, c1 and c2—constants in the pair potential function; ER—fitting error

新窗口打开| 下载CSV


表3   3种典型fcc结构晶体的实验数据(弹性常数的量纲为1011 Pa)[17,31]

Table 3  Experimental data for the three typical fcc structure crystals (elastic constants with dimension 1011 Pa)[17,31]

Crystala / nmutot / eVBC44Cp
Cu0.3615-3.541.4200.8180.257
γ-Fe0.3591-4.291.6701.1700.475
Ni0.3518-4.451.8761.3170.552

新窗口打开| 下载CSV


表4   3种典型fcc结构晶体的势参数的拟合值

Table 4  Fitted values of potential parameters for the three typical fcc structure crystals

CrystalcdAc0c1c2ER
Cu4.02054.01860.753712.1130-1.44710.251400.001806
γ-Fe3.99213.99720.717783.1839-2.30860.415890.002837
Ni3.91073.91660.764073.6419-2.70790.500410.013935

新窗口打开| 下载CSV


将参数的拟合值以及嵌入函数、对势函数和电子密度分布函数代入2.1和2.2节中的utot的表达式,将a记作R。得到utot关于R的解析表达式。图1和2分别给出了5种典型bcc结构晶体和3种典型fcc结构晶体的utot-R曲线。这些曲线中utot的最小值与utot的实验值是一致的,这个最小值所对应的R值与晶格常数一致。

图1

图1   5种典型bcc结构晶体的结合能随原子间距(utot-R)的变化曲线

Fig.1   utot-R curves for five typical bcc structure crystals (R—measurement of atomic spacing)


图2

图2   3种典型fcc结构晶体的结合能随原子间距(utot-R)的变化曲线

Fig.2   utot-R curves for three typical fcc structure crystals


需要说明的是,本工作的主要研究内容仅涉及弹性模量用势函数的导出问题,由于对势函数和电子密度分布函数采用了分段函数形式,忽略了截断距离外的势能,从而不可避免地会引起R偏离a时,utot的数值误差加大。

4 结论

(1) 给出bcc和fcc结构晶体适用的EAM/FS势函数和弹性常数之间的数学关系,其中C44Cp=(C11-C12)/2不仅依赖于所考虑的原点原子和周围原子的距离,还取决于近邻原子的排列状况。

(2) 对bcc和fcc结构晶体,利用嵌入函数、对势函数和电子密度分布函数导出了PBC44Cp的表达式。推导过程中,未定义嵌入函数、对势函数和电子密度分布函数的具体函数形式,具有普适性。

(3) 使用本工作提出的方法,验证了5种典型bcc结构晶体(V、Mo、Nb、Ta、W)和3种典型fcc结构晶体(Cu、γ-Fe、Ni)的结合能。在假设含待定常数的嵌入函数、对势函数和电子密度分布函数具体表达式的基础上,将关于utotPBC44Cp的5个拟合方程转化为求最小值的优化问题,从而准确地确定了结合能中的待定参数。

/