立方晶体弹性常数和EAM/FS势函数的关系
天津大学材料科学与工程学院水利安全与仿真国家重点实验室 天津 300354
Relationships Between Elastic Constants and EAM/FS Potential Functions for Cubic Crystals
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
基金资助: |
|
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: |
|
作者简介 About authors
段灵杰,男,1992年生,博士生
为了确定bcc和fcc晶体的EAM/FS势函数参数,研究了EAM/FS势函数与弹性常数的关系。对bcc和fcc结构晶体,分别导出压强(P)和体积弹性模量(B)、弹性常数(C44)和剪切弹性模量(Cp=(C11-C12)/2)的利用嵌入函数、对势函数和电子密度分布函数的表达式。发现C44和Cp的大小不仅取决于所考虑的原子和周围原子的距离,还受近邻原子排列状况的影响。将关于结合能(utot)和P、B、C44、Cp的5个拟合方程转化为求最小值的优化模型,给出了5种典型bcc结构晶体(V、Mo、Nb、Ta、W)和3种典型fcc结构晶体(Cu、γ-Fe、Ni)的结合能中待定参数的值。采用上述参数计算得到的最小结合能与实验结果一致,此时所对应的原子间距与晶格常数相同,表明了方法的有效性。
关键词:
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:
本文引用格式
段灵杰, 刘永长.
DUAN Lingjie, LIU Yongchang.
势函数被广泛应用于金属体系的分子动力学模拟,是计算过程中重要的基本参数之一[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]。
对于立方晶体,有3个独立的弹性常数C11、C12和C44。广义Hooke定律可写成[25]:
式中,σij为应力分量,εij为应变分量(i, j=1, 2, 3)。C11和C12可以通过体积弹性模量(B)和剪切弹性模量(Cp)表示:
应用应变张量
应用
本工作先推导出弹性模量和EAM/FS势函数之间的数学关系,即确定了C44和Cp与EAM/FS势的数学表达式,随后确定了bcc和fcc结构金属中与压强(P)、B、C44和Cp有关的嵌入函数以及对势函数和电子密度分布函数的表达式,最后分析讨论了势参数的拟合中应注意的问题,并通过拟合5种典型bcc结构金属和3种典型fcc结构金属的实验数据,确定出它们的结合能参数。
1 弹性常数和EAM/FS势
考虑bcc或fcc结构的理想晶体,其中所有原子是等价的。因此,一个原子的结合能(utot)可写成:
式中,uP是对势的贡献,uM是多体势的贡献,Σ表示对原点原子所有周围的相关原子求和,ri是所考虑的原点处的原子与周围的第i个原子间的距离,V(ri)是描述排斥的核-核作用的对势函数,ϕ(ri)是电子密度分布函数,f(ρ)是嵌入函数,ρ表示原点原子周围的相关原子在原点处产生的电子云密度之和。
用a表示晶格常数。对于bcc结构晶体,第一层有8个原子,它们与原点处的原子有相同的距离
P定义为:
而B可写成:
式中,Ω=a3/k是原子体积,对于bcc结构晶体,k=2;对于fcc结构晶体,k=4。
由式(6)得到W为:
这里引入了符号:
由式(4)和(5)可导出C44和Cp的表达式。为导出C44,引入位移:
这里x、y、z表示空间位置,ε12是小量,故有[30]:
对于距离
对于对势项,计算得:
这里(xi, yi)是第i个原子的坐标。对于多体项,得到:
用类似于式(14)和(15)的方法可得:
利用式(
这样C44为:
对于Cp,
这里ε11是小量,可得:
以及
类似于C44的导出,对势和多体势对Cp的贡献(CpP和CpM)分别为:
以及
Cp为:
可见,Cp和C44依赖于周围原子的坐标。
2 bcc和fcc结构晶体的弹性模量
在前文讨论的基础上,分别对bcc和fcc结构晶体考虑周围的两层原子来推导弹性模量的表达式。
2.1 bcc结构晶体(k=2)
对于bcc结构晶体,第一层的8个原子有坐标:
式(7)和(8)可写成:
从式(19)~(21)以及(25)~(27),计算C44和Cp,得:
2.2 fcc结构晶体(k=4)
对于fcc结构晶体,第一层中的12个原子有坐标:
由式(7)、(8)、(21)和(27)得到:
3 势参数的数值拟合
为完成数值拟合,首先需要确定嵌入函数、对势函数和电子密度分布函数的形式,按照紧束缚理论的经典模型,取
式中,常数c和d分别表示位于第二层和第三层之间的待定的与原点原子的距离,c0、c1和c2为待定系数。对于bcc结构晶体,
6个待定常数A、c、d、c0、c1和c2的值,需要通过拟合utot、B、C44和Cp的实验结果以及利用平衡条件P=0来确定。对bcc结构,在A>0,
表1 5种典型bcc结构晶体的实验数据(弹性常数的量纲为1011 Pa)[8]
Table 1
Crystal | a / nm | utot / eV | B | C44 | Cp |
---|---|---|---|---|---|
V | 0.30399 | -5.31 | 1.551 | 0.426 | 0.546 |
Mo | 0.31472 | -6.82 | 2.626 | 1.089 | 1.516 |
Nb | 0.33008 | -7.57 | 1.710 | 0.281 | 0.567 |
Ta | 0.33058 | -8.10 | 1.961 | 0.824 | 0.524 |
W | 0.31652 | -8.90 | 3.104 | 1.606 | 1.590 |
表2 5种典型bcc结构晶体的势参数的拟合值
Table 2
Crystal | c | d | A | c0 | c1 | c2 | ER |
---|---|---|---|---|---|---|---|
V | 3.1628 | 3.8799 | 1.7110 | 18.83 | -16.25 | 3.972 | 5.125×10-11 |
Mo | 3.2317 | 4.3625 | 1.3678 | 262.85 | -200.70 | 38.880 | 9.408×10-10 |
Nb | 3.3927 | 4.0949 | 2.5909 | 26.30 | -22.68 | 5.405 | 1.103×10-10 |
Ta | 3.4863 | 4.2968 | 2.0498 | 25.13 | -17.98 | 3.584 | 9.312×10-12 |
W | 3.2601 | 4.4763 | 1.5395 | 275.77 | -207.67 | 39.655 | 2.896×10-3 |
表3 3种典型fcc结构晶体的实验数据(弹性常数的量纲为1011 Pa)[17,31]
Table 3
Crystal | a / nm | utot / eV | B | C44 | Cp |
---|---|---|---|---|---|
Cu | 0.3615 | -3.54 | 1.420 | 0.818 | 0.257 |
γ-Fe | 0.3591 | -4.29 | 1.670 | 1.170 | 0.475 |
Ni | 0.3518 | -4.45 | 1.876 | 1.317 | 0.552 |
表4 3种典型fcc结构晶体的势参数的拟合值
Table 4
Crystal | c | d | A | c0 | c1 | c2 | ER |
---|---|---|---|---|---|---|---|
Cu | 4.0205 | 4.0186 | 0.75371 | 2.1130 | -1.4471 | 0.25140 | 0.001806 |
γ-Fe | 3.9921 | 3.9972 | 0.71778 | 3.1839 | -2.3086 | 0.41589 | 0.002837 |
Ni | 3.9107 | 3.9166 | 0.76407 | 3.6419 | -2.7079 | 0.50041 | 0.013935 |
将参数的拟合值以及嵌入函数、对势函数和电子密度分布函数代入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势函数和弹性常数之间的数学关系,其中C44和Cp=(C11-C12)/2不仅依赖于所考虑的原点原子和周围原子的距离,还取决于近邻原子的排列状况。
(2) 对bcc和fcc结构晶体,利用嵌入函数、对势函数和电子密度分布函数导出了P、B、C44和Cp的表达式。推导过程中,未定义嵌入函数、对势函数和电子密度分布函数的具体函数形式,具有普适性。
(3) 使用本工作提出的方法,验证了5种典型bcc结构晶体(V、Mo、Nb、Ta、W)和3种典型fcc结构晶体(Cu、γ-Fe、Ni)的结合能。在假设含待定常数的嵌入函数、对势函数和电子密度分布函数具体表达式的基础上,将关于utot、P、B、C44和Cp的5个拟合方程转化为求最小值的优化问题,从而准确地确定了结合能中的待定参数。