天津大学材料科学与工程学院水利安全与仿真国家重点实验室 天津 300354
为了确定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)的结合能中待定参数的值。采用上述参数计算得到的最小结合能与实验结果一致,此时所对应的原子间距与晶格常数相同,表明了方法的有效性。
关键词:
势函数被广泛应用于金属体系的分子动力学模拟,是计算过程中重要的基本参数之一[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个独立的弹性常数C11、C12和C44。广义Hooke定律可写成[25]:
式中,σij为应力分量,εij为应变分量(i,j=1, 2, 3)。C11和C12可以通过体积弹性模量(B)和剪切弹性模量(Cp)表示:
B和Cp是便于计算的[26,27]。对于C44和Cp,考虑应变能密度
应用应变张量
应用
本工作先推导出弹性模量和EAM/FS势函数之间的数学关系,即确定了C44和Cp与EAM/FS势的数学表达式,随后确定了bcc和fcc结构金属中与压强(P)、B、C44和Cp有关的嵌入函数以及对势函数和电子密度分布函数的表达式,最后分析讨论了势参数的拟合中应注意的问题,并通过拟合5种典型bcc结构金属和3种典型fcc结构金属的实验数据,确定出它们的结合能参数。
已往文献[8,20,28,29]中只列出势参数的最后拟合结果,很少具体给出势函数与弹性模量之间的数学表达式,导致势函数的可靠性很难验证,也缺乏普适性。
考虑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依赖于周围原子的坐标。
在前文讨论的基础上,分别对bcc和fcc结构晶体考虑周围的两层原子来推导弹性模量的表达式。
对于bcc结构晶体,第一层的8个原子有坐标:
式(7)和(8)可写成:
从式(19)~(21)以及(25)~(27),计算C44和Cp,得:
对于fcc结构晶体,第一层中的12个原子有坐标:
由式(7)、(8)、(21)和(27)得到:
为完成数值拟合,首先需要确定嵌入函数、对势函数和电子密度分布函数的形式,按照紧束缚理论的经典模型,取
式中,常数c和d分别表示位于第二层和第三层之间的待定的与原点原子的距离,c0、c1和c2为待定系数。对于bcc结构晶体,
6个待定常数A、c、d、c0、c1和c2的值,需要通过拟合utot、B、C44和Cp的实验结果以及利用平衡条件P=0来确定。对bcc结构,在A>0,
式中,
表15种典型bcc结构晶体的实验数据(弹性常数的量纲为1011Pa)[8]
Table 1
表43种典型fcc结构晶体的势参数的拟合值
Table 4
将参数的拟合值以及嵌入函数、对势函数和电子密度分布函数代入2.1和2.2节中的utot的表达式,将a记作R。得到utot关于R的解析表达式。图1和2分别给出了5种典型bcc结构晶体和3种典型fcc结构晶体的utot-R曲线。这些曲线中utot的最小值与utot的实验值是一致的,这个最小值所对应的R值与晶格常数一致。
图15种典型bcc结构晶体的结合能随原子间距(utot-R)的变化曲线
Fig.1utot-Rcurves for five typical bcc structure crystals (R—measurement of atomic spacing)
图23种典型fcc结构晶体的结合能随原子间距(utot-R)的变化曲线
Fig.2utot-Rcurves for three typical fcc structure crystals
需要说明的是,本工作的主要研究内容仅涉及弹性模量用势函数的导出问题,由于对势函数和电子密度分布函数采用了分段函数形式,忽略了截断距离外的势能,从而不可避免地会引起R偏离a时,utot的数值误差加大。
(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个拟合方程转化为求最小值的优化问题,从而准确地确定了结合能中的待定参数。
1弹性常数和EAM/FS势
2 bcc和fcc结构晶体的弹性模量
2.1 bcc结构晶体(k=2)
2.2 fcc结构晶体(k=4)
3势参数的数值拟合
图1
图2
4结论
来源--金属学报