分子动力学模拟位错和界面的相互作用

1、分子动力学模拟位错和界面的 相互作用学校代码 10530学 号 200910081121分类号 密 级硕士学位论文分子动力学模拟位错和界面的相互作用学位申请人周银库指导教师陈尚达 副教授学院名称材料与光电物理学院学 科 专业材料科学与工程研 究 方 向金属薄膜的力学性能二零一二年五月Molecular dyn amics simulati ons of in teractio n betwee n dislocati ons and in terfacesCan didateYinku ZhouSupervisor Shangda Chen (AssociateProfessor)Colleg

2、e Faculty of Materials, Optoelectronicsand PhysicsProgramMaterial Science andEngin eeri ngSpecializatio nMechani cal Properties ofmetal filmDegreeEngin eeri ng MasterUn iversityXiangtan Uni versityDateMay, 2012湘潭大学学位论文原创性声明本人郑重声明：所呈交的论文是本人在导师的指导下独立进行研究所取得的研究成果。除了文中特别加以标注引用的内容外，本论文不包含任何其他个人或集体已经发表或撰写

3、的成果作品。对本文的研究做出重要贡献的个人和集体，均已在文中以明确方式标明。本人完全意识到本声明的法律后果由本人承担。作者签名：日期： 年 月 日学位论文版权使用授权书本学位论文作者完全了解学校有关保留、使用学位论文的规定，同意学校保 留并向国家有关部门或机构送交论文的复印件和电子版，允许论文被查阅和借 阅。本人授权湘潭大学可以将本学位论文的全部或部分内容编入有关数据库进行 检索，可以采用影印、缩印或扫描等复制手段保存和汇编本学位论文。涉密论文按学校规定处理。导师签名：作者签名：日期： 年 月 日日期： 年 月摘要纳米尺度的金属多层膜在屈服应力、塑性、 抗腐蚀性能等方面具有特殊的性能。 目 前

4、它已被广泛应用于航空航天、机械制造、 电子技术、 光学工程及计算机工程等各个领 域。而在薄膜材料的应用过程中， 薄膜的使用寿命和可靠性是人们普遍关注的焦点问题。 界面的结合性能是影响多层膜寿命和可靠性的关键指标， 而位错和界面的相互作用机理 决定着界面的结合性能， 即位错和界面的相互作用机理在薄膜的使用寿命和可靠性方面 扮演着关键角色。因此对位错和界面的相互作用机理的研究就显得特别有价值和意义。 随着高性能计算机的发展，原子模拟已成为材料性能预测与设计方面一种有效的方法。 本文用三维分子动力学方法研究了位错和界面的相互作用机理，具体如下：首先，用分子动力学方法研究了侧向拉伸载荷下位错从 bcc

5、-Fe/Ni 界面的形核和发 射过程。弛豫后，在 Fe(0 0 1)/Ni(0 0 1)和 Fe(0 0 1)/Ni(1 1 1)界面观察到无序的失配位错 网络， Fe(0 0 1)/Ni(1 1 0) 界面观察到长方形的失配位错网络。研究了晶体取向对 Fe/Ni 双层膜拉伸性能的影响。不同取向的对比发现 Fe(0 0 1)/Ni(1 1 0)系统的屈服强度最低。 和Fe薄膜进行了对比，发现Fe/Ni双层膜系统的塑性高于Fe薄膜的，而屈服强度低于 Fe 薄膜的。模拟结果显示，界面是位错的发射源，滑移位错从界面的失配位错线形核 和发射。同时界面也会阻碍位错运动，随着拉伸的进行， Fe 层中越来越

6、多的位错被塞 积在界面处，当到达到临界值时，迫使位错穿过 Fe/Ni 界面，从 Fe 层到 Ni 层。在 Fe 基体中位错主要在 1 0 1面滑移，而在 Ni 中主要在 1 1 1面滑移。其次，用分子动力学模拟了单轴拉伸载荷下不同扭转角的 Cu(001)/Ni(001) 界面的结 合性能。模拟结果显示，当扭转角小于 15.124度时，界面形成方格状的失配位错网络， 界面失配位错网络的密度随着扭转角的增加而增加。当扭转角大于 15.124 度时，在界 面形成面缺陷。 模拟发现界面构型对 Cu/Ni 系统的界面强度有着非常显著的影响。 随着 扭转角的增加屈服应力首先减小，直到扭转角为 5.906度

7、的最小值，然后增加，当其达 到扭转角为 15.124 度的最大值后，又开始减小 ，最后当扭转角约大于 20 度，屈服应 力几乎趋于一稳定的值。关键词： 分子动力学；界面；滑移位错；失配位错ABSTRACTNanoscale multilayered composites often possess extraordinary mechanical properties in terms of yield stress, ductility, and wear resistant. Now it has been widely used in aerospace, mechanical manu

8、facturing, electronics, optical engineering and computer engineering fields. In the applied process of the thin film materials, people commonly focus on the reliability and service life of the thin films. The binding property between the thin film and the substrate is a key indicator of metallic mul

9、tilayers 'the reliability and service life, the interaction mechanisms between dislocations and interfaces dominate the binding properties of the interfaces. Therefore, the interaction mechanism between dislocations and interfaces plays a vital role in the field of the reliability and service li

10、fe of the thin films.So it is an interesting and valuable thing to understand the interaction mechanisms between dislocations and interfaces. With the development of high-performance computer, atomic simulations have become an effective method in the field of material properties forecast and design.

11、In the presentwork, we have studied the interaction between dislocations and interfaces with 3D Molecular Dynamic Simulations.Firstly, molecular dynamics simulations were carried out to investigate the nucleation and emission of dislocations from an interface in a bcc-Fe/Ni bilayer subjected to tran

12、sverse loading. After relaxation, disordered types of dislocations were observed at both Fe(0 0 1)/Ni(0 0 1) and Fe(0 0 1)/Ni(1 1 1) interfaces, and rectangular dislocations types at Fe(0 0 1)/Ni(1 1 0) interface. The orientation effect on the mechanical properties of a Fe/Ni bilayer system was inve

13、stigated. The yield stress of the Fe(0 0 1)/Ni(1 1 0) system abtained is lowest. We also found that the yield stress of pure iron nanofilm was higher than that of a Fe/Ni bilayer system, and the ductility was lower than that of a Fe/Ni bilayer system for given temperature and strain rate. The simula

14、tion results obtained also show that the misfit dislocations at Fe/Ni interface acted as a source to nucleation and emission of glide dislocations. Glide dislocations nucleation and emission from misfit dislocation line at Fe/Ni interface. The existence of misfit dislocations and the lattice mismatc

15、h can also act as barriers to dislocation motion and transmission across the interface. More dislocations in Fe have been arrested at the Fe/Ni interface, which provides sufficient stress for dislocations to transmit from Fe toNi. Glide dislocations mainly occurred on 1 0 1 plane in Fe layer of FeNi

16、 bilayer, and 1 1 1 plane in Ni layer.Secondly, molecular dynamics simulations were carried out to study the mechanical properties of Cu(001)/Ni(001) interface boundaries with different twist angles subjected to uniaxial loading. The results obtained revealed that square misfit dislocations networks

17、 can be observed when the twist angle was lower than 15.124 , and the density of misfit dislocations increased with increasing twist angle. Face defects were formed when the twist angle was higher than 15.124 . It has been found that the interface configuration had a significant effect on the interf

18、ace strength of the Cu/Ni system. The yield stress was found to decrease first with increasing twist angle and it reached its lowest value at 5.906 twist angle. Subsequently, it increased with increasing twist angle till it reached its highest value at 15.124 of the latter; it then decreased again a

19、nd finally became almost constant when the twist angle was larger than approximately 20 .Key Words: Molecular Dynamics; Interface; Glide dislocation; Misfit dislocation5第1章引言11.1薄膜概述 11.2界面概述 11.3晶体位错相关理论概述 41.4晶体界面的分子动力学研究现状 61.5本文的研究思路、目的及意义 7第2章分子动力学方法92.1引言 92.2基本原理102.2.1积分方法102.2.2原子间的相互作用势 11

20、2.3边界条件 132.4温度、压力控制方法 142.4.1控温方法152.4.2控压方法162.5分析方法 162.5.1 中心对称参数(centrosymmetry parameter ) 162.5.2径向分布函数 172.6模拟软件与可视化 182.7小结 18第3章Fe/Ni双层膜系统中位错和界面的相互作用 203.1模拟方法 203.2结果及讨论 213.2.1弛豫后沿厚度方向的势能分布 223.2.2初始界面失配位错结构 223.2.3晶体取向的影响 233.2.4 bcc-Fe薄膜的变形机理 253.2.5 bcc-Fe层的滑移位错 263.2.6 Fe/Ni双层膜的变形机理

21、283.2.7 Ni层厚度的影响 31328温度的影响32329应变率的影响333.3小结 34第4章Cu(001)/Ni(001) 扭转界面的结合强度 354.1模拟模型和方法 354.2结果与讨论 374.2.1初始失配位错构型374.2.2 Cu(001)/Ni(001)系统的变形机理 384.2.3 Cu层中滑移位错的变形机理 404.3小结 39第5章 总结与展望 435.1工作总结 435.2工作展望 44参考文献 45致谢55个人简历、攻读硕士学位期间发表的论文 56#湘潭大学硕士毕业论文第 1 章 引言1.1 薄膜概述薄膜材料通常是指在二维方向上的尺度比另一维方向大很多，在其基

22、体材 料上通过化学或物理等方法制备另外一层材料，以达到某种特定功能，如提高 断裂韧性、增加塑性、增加使用寿命 , 以及美观度等等。薄膜材料与块体材料相 比，有比较高的比表面积，比较少的配位数，由于原子之间的化学键在表面处 突然断开，故使其表面原子相对内部原子有较高的表面势能和表面应力1-3 。有很多种薄膜的制备方法 4 ，如物理方法中的真空蒸发、溅射、离子束和离 子助、外延膜沉积技术等，以及化学方法中的电镀、阳极反应沉积法、化学气 相沉积等。各种制备方法有各自的优缺点和适用范围，实际应用时以具体情况 而定。薄膜的分类方式有多种 5。按性质来分，有合成膜与天然膜等 ; 按层数 ,有 单层膜和多层

23、膜等。由于具有独特的微观结构以及物理和化学性能，薄膜材料在最近几十年来 已被广泛应用于航空航天、机械制造、电子技术、光学工程及计算机工程等各 个领域。目前许多国家都把薄膜材料的研究设为大型研究项目。在各种类型薄膜中，金属多层膜在屈服应力、塑性、抗腐蚀性能等方面具 有特殊的性能 6-11。既而金属 /金属薄膜材料被广泛应用于各种结构功能材料中， 薄膜材料应用的可靠性和使用寿命很大程度上依赖于薄膜与基体的结合性能， 据统计，日、美、欧共体等国每年国民生产总值的 6%8% 12 都因材料的疲劳、 断裂、腐蚀、磨损等破坏而损失。而界面的结合性能是影响多层膜质量的关键 指标，这就要求我们对薄膜界面微观结

24、构有很好的了解。对界面结合性能的研究，已经是当今具有困惑性和挑战性的难点问题。目 前已有大量的实验研究方法，如纳米压痕法、鼓包法、拉伸法等。另外，随着 近年来高性能计算机的发展，计算机模拟方法也被越来越广泛的应用于这一领 域。1.2 界面概述界面一般包含表面、晶界和相界面等，为一种二维缺陷，即面缺陷。我们 称晶体与空气或液体接触的界面为表面。在表面上的原子，其相邻原子数比晶体内部要少，相当于一部分结合键被拆断，因而有较高的能量，产生了表面能， 表面能通常要比晶界能更大些。通常应用中的金属材料，大部分为多晶材料， 而非单晶材料。在多晶材料中，在晶粒和晶粒之间存在着界面，即晶界。晶界 处晶粒从一个

25、位相过渡到另一个位相，晶界把结构相同位相不同的亚晶粒隔开。 由于晶界处同时受到两侧晶粒的影响，处于两侧晶粒的过度状态，从而使晶界 具有特殊的构型，进而对晶体的结合性能产生着巨大的影响。如图1.1所示为Cu晶体中艺9（221） 9 =141.1的晶界结构模型13。上下两部分晶粒的晶向不同， 交叠在一起，在界面处形成 E结构，晶界处为两边晶粒的过度状态。晶界对多 晶材料的物理、化学、力学等性质有着非常显著的影响。材料的强度和断裂等力学行为，如偏聚、晶界扩散，以及晶界处初始滑移位错的形核等，都受到晶界结构的显著影响图1.1 Cu晶体中工9（221） 9 =141.1的晶界结构模型13晶界按晶粒间取向

26、差的大小可分为小角度晶界（取向差小于5 度）和大角度晶界（向差大于10度）。小角度晶界通常有可分为扭转晶界和倾侧晶界两种类型， 如图1.2和图1.3所示。重合位置点阵模型（CSL） 14:即在一些特殊位相的晶界处，有一些原子 同属于两边晶粒的格点，且自身形成了超晶格点阵模型。图1.2为面心立方结构 中的（001 ）面重位扭转晶界15，即当旋转角9=6.9度时所得的扭转界面，其中 重合点阵的格点（即图1.2中的大黑点）的数目相当于总格点数目的1/5，即为重合密度。设 工为重合密度的倒数。工越小，即界面处重合的密度越大，重合24的原子数也就越多。图1.3为<100>38度时，重位倾侧界

27、面示意图15。图中AC为阶，阶高为BC,阶长为AB。显然，阶越小，即图1.3中AC越短，重合密度就越高，工也 越小，阶中不接触的原子也越少，即晶界能也较小。图1.3 <100>38时，倾侧界面示意图15相界为不同两相所形成的界面，相邻两相不仅取向不同，而且结构、成分也 不同。按照原子在相界上排列不同，可把相界分为三种形式 16:FCC(a) 共格相界面：界面两边金属具有同种类型的晶格结构，如都为结构，且其晶格参数相差不大；(b) 半共格界面：界面两边金属，其晶格类型相同，但其晶格不匹配程度 稍微大点；(c) 非共格界面：界面两边金属具有不同的晶格结构，如 FCC/BCC界面。通常将

28、失配度定义为一个量= 2|d1d2|/(d1 + d2)，这里d1, d2分别晶体1和晶体2中原子间的距离。两近邻平行失配位错间的距离 L=b/ §这里b为失 配位错的 Burgers 矢量。在共格和半共格晶界或是相界面上， 由于晶格不匹配， 经弛豫后， 在界面处 会形成失配位错网络，而这些界面失配位错网络对界面的初始滑移位错的形核 和发射起着决定性作用，从而决定着界面的结合强度。金属多层膜之所以具有其都特的性能， 主要归因于存在着界面。 界面既可以 阻碍位错运动， 也是位错的发射源 17-21。位错和界面相互作用决定着金属多层膜 机械性能。关于位错与界面的相互作用机理的研究就变得非

29、常有价值和意义。1.3 晶体位错相关理论概述位错是指晶体中滑移部分与未滑移部分的分界线， 是晶体中原子的局部不规 则排列，为一维缺陷或线缺陷。晶体的塑性变形通过位错滑移来实现，若没有 位错的运动，变形就不能发生，位错滑移通常在特定的滑移系上进行。其存在 对材料的物理性能、化学性能，以及机械性能产生着巨大的影响。二十世纪初 期， Volterra 等人 22在连续弹性力学中提出位错的概念。实际应用中的材料， 其内部通常含有这样或那样的缺陷。 缺陷的存在对材料 各方面性能会产生巨大的影响。如会降低材料的使用寿命、力学性能等等。但 是有时候材料中存在缺陷又会增加其某方面的性能，如在半导体中掺进微量元

30、 素来提高其性能，加工增加位错来使金属材料强化等。因此，为了使材料更好 的为人类服务，对材料中缺陷的研究就显得非常重要。晶体力学性质的微观理论和位错有着非常密切的联系。 晶体的塑性变形通过 滑移来完成，在不同的滑移系下进行，而晶体滑移借助滑移位错来实现。近年 来，随着晶体力学性质的微观理论和实验的发展，进而促进了晶体缺陷相关理 论的发展。根据局部滑移的方式不同， 位错可分为刃型位错、 螺型位错，以及混合位错，如图1.4所示。图1.4(a )为刃型位错模型，从图中可以看出，有一额外的半原子 面终止在晶体中形成一条线缺陷，即位错线，图中用符号壬表示。由图1.4(a )可以看出位错线上部分原子受到压

31、缩，而下部分受到拉伸，在位错线处形成局部应变场，导致位错线处晶格的对称性遭到破坏，该处相对完整晶体，势能会增 加。当额外的半原子面在晶体的下半部分时，用 T标示。其多出的半原子面如同 刀刃一样插入晶体中，故其形成的位错称为刃型位错。(a)刃型位错23(b)螺型位错24(c)混合位错24图1.4三种类型位错示意图图1.4(b)为晶体中螺型位错的晶体外观。晶体上下两部分沿滑移面相互滑移 了一个原子距离，滑移终止在晶体中。螺型位错也是一种线缺陷，如图1.4(b)中的AB,即为位错线，是已滑移区与未滑移区的分界。大部分晶体材料中的位 错既不是纯粹的刃型位错也不是纯粹的螺型位错，他们往往是这两种类型都有

32、， 即为混合位错，如图1.4(c)所示。通常用柏氏矢量来表示晶体滑移的量。下面以刃型位错为例介绍一下柏氏矢 量。在图1.5(a)中封闭回路MNOPQ，而在图1.5(b)中相应的线路是不封闭的， 始点M，终点Q。我们定义：将终点联到始点的矢量(图中的 QM )，称为柏氏矢量。在图1.4三种类型位错示意图中，其中矢量b即代表柏氏矢量151.4晶体界面的分子动力学研究现状在过去十年里，对金属多层膜各方面性能的理解取得了巨大的进展。纳米 尺度的金属多层膜在屈服应力、塑性、抗腐蚀性能等方面具有特殊的性能6-11 0这种特殊性能主要归因于多层膜中存在着界面，界面既可以阻碍位错运动，也 是位错的发射源17-

33、21。获得精确工程界面的能力在当今新技术发展方面扮演着越 来越重要的角色25。位错和界面的相互作用对金属多层膜机械性能的影响依赖 于界面的类型。过去几十年里，分子动力学方法已被广泛应用于材料科学、物理学、化学、 生物学，以及医学等领域。随着高性能计算机的发展，原子模拟在材料性能预 测与设计方面已成为一种有效的方法。已经有大量的研究者用分子动力学方法 研究晶界和界面处位错的形核和发射过程26-34 0女口，Tschopp等人35用原子模拟 方法研究了纳米晶Cu中位错的形核和发射过程，表明晶界是位错形核和发射源；Spearot等人28用分子动力学模拟方法研究了单轴拉伸下位错从双晶界面的 形核和发射

34、过程；Derlet等人17用原子模拟方法模拟了位错和 Cu晶界的相互作 用，表明晶界既可以阻碍位错运动，也是新位错的发射源；Shao等人26用分子 动力学模拟研究了 Ni/Cu双层膜的纳米压痕过程，结果表明界面会阻碍滑移位 错位错穿过界面；Che ng等人36用分子动力学方法研究了Cu/Ni双层膜界面的初始失配位错网络，以及其对摩擦过程的影响，结果表明界面失配位错网络对 位错穿过界面有明显的阻碍作用；Zheng等人19用分子动力学方法模拟了 Ni/Al 双层膜的单轴拉伸过程，表明初始滑移位错从界面的失配位错线形核和发射等。然而，大部分这些研究主要关注于同种类型的材料和非扭转界面，如晶界、 FC

35、C/FCC 界面。不同类型的材料，如 FCC/BCC 界面，以及扭转界面被研究的 相对较少。对于 FCC/BCC 界面有如下： Hoagland 等人对非共格的 Cu/Nb 界面 进行了研究 18,37,其主要关注于在二维条件下位错和界面的相互作用；Wang 等人 38用类似方法更细致地研究了位错和不同类型界面的相互作用，认为Cu/Nb作为非共格界面系统，界面会阻碍滑移位错的发射；更近一点，Shao 等人16用分子动力学方法研究了 Cu/Nb 双层膜的纳米压痕过程，表明此非共格界面对滑 移位错的增殖起到强烈的阻碍作用。而对于扭转界面和扭转晶界的研究也非常 少，这其中包括在铜扭转晶界的晶界迁移和

36、晶界自扩散的研究39 ；未弛豫的Ag(001)/Ni(001)扭转界面能分析40;不同扭转角度下Cu纳米线机械性能的研究 41等。这些工作对 FCC/BCC 界面，以及扭转界面和扭转晶界的研究提供了非 常宝贵的信息。本文用分子动力学模拟研究了位错和 BCC-Fe/Ni 界面的相互作用，以及不 同扭转角下Cu(001)/Ni(001)双层膜的结合性能。1.5 本文的研究思路、目的及意义薄膜材料在当今社会中具有重要的应用价值， 目前它已被广泛应用于航空航 天、机械制造、电子技术、光学工程以及计算机工程等各个领域。而在薄膜材 料的应用过程中，薄膜的使用寿命和可靠性是人们普遍关注的焦点问题。而界 面的

37、结合性能是影响多层膜寿命和可靠性的关键指标。位错和界面的相互作用 机理决定着界面的结合性能。即位错和界面的相互作用机理在薄膜寿命和可靠 性方面扮演着关键角色。因此对位错和界面的相互作用机理的研究就显得特别 有价值和意义。位错和界面的相互作用机理虽然已经有很多科研工作者在研究， 但是大部分研究集中在位错和 FCC/FCC 界面，以及位错和晶界间的相互作用， 而位错和 FCC/BCC 界面，以及位错和扭转界面的相互作用机理研究的相对较 少。同时，近年来分子动力学模拟方法已成功应用于位错和界面相互相互作用 的研究中。Fe-Ni 合金作为一种重要的结构功能材料，由于具有较好的抗腐蚀性能、断 裂韧性，以

38、及较好的塑性 42,43，而被广泛应用于许多领域。下面简要介绍下本文的研究思路、研究目的和研究意义。研究思路是首先建立 Fe/Ni 双层膜模型和 Cu/Ni 扭转界面薄膜模型。 然后用 分子动力学方法模拟了这两模型的拉伸过程。 Fe/Ni 双层膜采用侧向拉伸，拉伸 应力仅应用于 Fe 薄膜一侧。 Cu/Ni 扭转薄膜模型拉伸过程中，沿轴向拉伸。研 究了弛豫后界面的初始失配位错网络，以及其对界面位错的形核和发射的影响， 进而对界面结合性能的影响，以揭示位错和界面的相互作用机理。研究目的是用计算机模拟方法研究位错从界面失配位错线形核和发射过程 的微观机理 , 以及如何设计薄膜材料使其界面结合强度最

39、高， 使材料更好地为人 类服务。本文研究的意义在于用分子模拟来研究金属 /金属薄膜材料界面的结合性 能，揭示界面位错的微观演化机理，为实验和理论研究金属 /金属薄膜材料界面 结合性能提供一定的支持和帮助，使之为人们更好的服务。第 2 章 分子动力学方法2.1 引言 如今，随着科学技术的不断发展，人们对自然世界的探索已不局限于宏观 领域，越来越趋于微观领域，如分子、原子尺度。由于宏观物质是由无可计数 的原子、分子构成，处理起来十分复杂，甚至超越了人的计算能力，以致处理 起来十分复杂。直至计算机模拟技术的出现 , 是这一问题处理起来方便许多 , 另 外也降低了成本，从而大大的促进了微观领域的发展。

40、计算模拟是化学、材料 学、医学、生物学等多个领域的交叉学科，是一门正在蓬勃发展的新兴学科， 是利用一定的算法结合相关计算机模拟软件对物质的组成、微观结构、性能以 及服役性能进行计算模拟与设计的学科。分子动力学模拟 44 是指用计算机模拟原子或分子的运动过程, 以揭示材料内在的本质规律。近年来，分子动力学模拟已成为纳米领域研究最常用的计算 方法。与蒙特卡洛方法、第一性原理，以及有限元一起成为计算机模拟技术的 主流技术，被称为计算机实验。计算机实验在当今科学研究领域扮演着重要角 色。在计算机实验中，通过建立理论模型，并借助计算机以及一定的算法和程 序来计算。其中，蒙特卡洛方法是研究随即性问题最有效

41、的方法，可以研究材 料介观尺度问题。第一性原理方法可以无需任何实验数据，完全从材料组成的 种类以及排列方式出发计算材料性能，其研究体系比分子动力学小得多。有限 元则是研究宏观尺寸的有效方法。这些计算机模拟方法被广泛的应用于科学研究，为理论研究和实验研究搭 建桥梁，成为除理论研究和实验研究外探索自然界的另一种有效方法。模拟有 时是一种理论，但有时又是一种实验。一方面，我们处理的是模型，而不是真 实的物体，这表明模拟属于一种理论方法；另一方面，通过计算机模拟证实一 个模型的程序又非常类似于实验：我们运行一个程序，然后分析结果，其方法 非常类似于实验过程。计算机模拟、理论研究和实验研究三者相互联系，

42、相互 补充，相互促进，共同推动着科学技术的发展。计算机模拟可以方便的观察到许多实验无法观察到的细节，甚至实验很难 达到或是无法达到的条件，揭示材料的本质规律。鉴于计算机模拟的优点，本 文使用分子动力学方法来研究位错和 Fe/Ni界面的相互作用机理，以及Cu/Ni界面 的结合强度。2.2基本原理分子动力学模拟通过用经典的牛顿运动方程来描述粒子的运动过程。221积分方法常用的积分方法有 Velet算法、Leap Frog算法、Velocity-Verlet算法和 Beeman算法等。 Verlet算法在分子动力学模拟中，最常用的时间积分方法为Verlet45算法。其对粒子的坐标进行泰勒展开得到式

43、2.1和式2.2：ri(tt) ri(t)tvi(t) 4 t2ai(t)+L(2.1)ri (tt) ri(t) tVi(t) 1 t2ai(t) L(2.2)式中r代表位移，v代表速度，a代表加速度，t代表时间。将式2.1和式2.2相 加并忽略高次项得到t+ S时刻的位置，即式2.3。ri(tt) 2ri (t) ri(tt) t2凹mi(2.3)式2.3即为Verlet算法的基本形式。将式 2.1和式2.2相减并忽略高次项可得粒 子运动的平均速度Vi(t) ri(t t) ri(t t)/2 t(2.4)或者Vi(t t/2) ri(tt) ri(t)/ t(2.5)Verlet算法具有

44、算法简单，存储要求适度等优点，但容易造成精度损失，且 非自启动算法。(2) Leap Frog 算法Leap -Frog算法46，即蛙跳算法。其表达式如下：Vi(tt/2)Vi(tt/2)/tFi(t)/mi(2.6)即得到Vi(tt/2)Vi(tt/2)tFi(t)/mi(2.7)由式2.5得ri(tt)ri (t)t Vi(tt/2)(2.8)Leap-frog算法比Verlet算法计算效率高、计算量小。但其缺点是位置和速 度不是同步的。(3) Velocity-Verlet 算法Velocity-Verlet算法47其优点是计算量适中，且给出了显式速度项，目前得 到了广泛的应用。其表达式

45、如下：斗(t)导h(2,h2miri (t h) Vi(t)F(t h) F(t)h2mi(2.10)在Velocity-Verlet算法中前一个时刻的速度、力和位置需要被储存。(4) Beeman 算法Beeman算法48也是在Verlet算法基础上发展的另一种计算精度更高算法, 但增加了其计算量。2.2.2原子间的相互作用势在分子动力学模拟中，势函数是影响计算结果的关键因素之一，其决定着 计算结果能否准确地反映实际材料的性质。原子间的相互作用力也可由势函数 求得，若势函数存在误差，贝U所计算的相互作用力也会受到影响。势函数的准 确性决定着分子动力学模拟的精度。(a)对势对势通常有Lenna

46、rd-Jones势、Morse势和Johnson势等势模型。 Lennard-Jones ( L-J )势49:(rj) 4 ( /汀2 ( /汀(2.11)式中 为能量参数，rij为原子i和原子j之间的距离，(为长度参数。当rj=o时， 势能为零；当口？(时，表现为引力势，由式中第二项来表示，其最小值为1.122(X当距离较远时，原子间的引力会越来越弱，趋近于零，通常为了计算方便，选 取截断半径rc,当rij?-c时，取势能为零，在Lennard-Jones势模型中，截断半径rc 通常取21/6 X；当口？(时，表现为斥力势，由式中第一项来表示，当原子间的距离 越近，部分电子云发生重叠，原子

47、间的斥力会越来越大。(2) Morse 势50:ij (ry) Ae2 (rij r0)2e ("j r0)(2.12)式中A、a ro分别是结合能，调节参数和平衡间距，这三个值都是取经验 参数。金属固体常用这种势函数来描述，常用于金属铜。 Johnson 势513(rij)An (rij Bn)5币 Dn(2.13)aFe常用Johnson势函数来描述。(b)多体相互作用势实际上，在多粒子系统中，就不能简单的应用以上对势模型，而必须考虑 多体间的相互作用。常见的多体势有嵌入原子法、Stillinger-Weber势、 Finnis-Sinclair 势等。(1) 嵌入原子势(EAM

48、 ) 52在嵌入原子势(Embedded Atom Method,简称EMA )中，材料中的每一个 原子被看作如同一个杂质原子镶嵌在包含所有其它原子的基体中。EMA势可以处理各种各样的金属体系，其中包括断裂、表面、杂质，以及合金等体系。系 统的总能量如下：1U ;ij (rj )+ Fi( i)(2.14)2 j iiEAM势是将晶体的总势能分成两部分组成，上式中右端第一项曲代表对势相互作用，口为原子i和原子j之间的距离。第二项Fi即嵌入能，是电子密度 p的函数，代表多体相互作用。式中1/2表示对势为两原子共有，式中 p可以表 达为：ij(rj(2.15)j i对于不同的金属，EAM势需通过拟

49、合金属的宏观参数而得到。(2) Stillinger-Weber 势53Stillinger-Weber势常用于描述半导体材料。1 1 2 V(rij)g(ij)g(rik) cos jik(2.16)2 ijijk3式中 加为ij键和ik键之间的夹角，g(r)为截断半径在第一近邻和第二近邻 之间的衰变函数。截断半径的选取是为计算方便，因为当距离非常远时，原子 间的相互作用势非常小，趋近于零，当截断半径取得越大，结果越精确，但计 算量也越大。当cosjk=-1/3时，这种结构类似金刚石四面体结构，其结构非常 稳定。(3) Finnis-Sinclair 势54以合金属能带紧束缚理论为基础而发展

50、的Finnis-Sinclair势，给出了详细的多体作用势。其表达式如下：2j (rij )Ak(Rk rij) H (Rk rij )(2.17)k 163ij (rj)ak rj) H(m)(2.18)k 1式中，当 x>0 时，H(x)=1 ;当 x<0 时，H(x)=0。 Ak、Rk、ak、rk 为常数，且 有R1>R2，r1>r2>->r6。它们的值随随具体应用的材料而定2.3边界条件选取合适的边界条件，对分子动力学模拟非常重要。但通常主要考虑以下 两点：首先，为了减小计算量，应尽可能减小模拟的单元，但是为了排除可能 的动力学扰动对结果造成影响，以

51、及为了满足统计学处理的可靠性要求，又要 求模拟原胞应该足够大；其次，要从物理角度考虑应变相容性、体积变化及环 境的应力平衡等实际耦合问题。当模拟大块体时，需使用三维周期性边界条件；当模拟薄膜材料时，需使 用二维周期性边界条件，如取 X、丫方向为周期性边界条件，Z方向为自由边界 条件；当模拟纳米线时，可选取一维周期性边界条件，如取Z轴方向为周期性边界条件，X、丫方向为自由边界条件；有时又需使用混合边界条件，即几种边 界条件的组合。实际应用时应根据实际情况选择合适的边界条件。由于实际的宏观物质由无可计数的原子、分子组成，而计算机的计算能力 又有限，所以为了进行模拟，就必须采用周期性边界条件。当使用

52、周期性边界 条件时，实际的粒子被假想在一个盒子里，这个盒子在空间三维方向被重复无 穷多次，填充满整个空间。即如果在盒子中位置r处有一粒子，那在式2.19所在位置处也能找到相同的粒子。r la mb nc, (l,m,n , )(2.19)这里l, m, n为任意整数，a, b,c为大小对应盒子边长的矢量。 所有这些镜像 粒子一起运动，但实际在计算程序中的只有一个粒子。这样大大的减小了计算 量，而且也不失精度。图2.1为分子动力学元胞三维周期性边界的二维图示。在基本分子动力学元胞中，当有原子离开分子动力学元胞时就有对应原子进入其中。OO O 00 000 0d 0 0L L L 00 000 0

53、00 000 0图2.1分子动力学元胞三维周期性边界的二维图示2.4温度、压力控制方法在分子动力学模拟中，我们往往要考虑模拟系统是否与外界发生作用，如温度、压力等。对于不同的模拟对象，往往要对原子数 N、体积V、压强P、温 度T等量中的一些进行控制，以实现不同的系综。如在NVT系综下需要控制原子数N、体积V、温度T保持不变。而在NPT系综下则需要对原子数 N、压强P、体积V进行控制241 控温方法在NVT、NPT系综下，甚至NVE系综下，我们期望模拟体系温度保持恒 定不变。而系统的温度T与动能K及速度有直接的关系，如下，N2Kmi (vj / 2 (3N 比山丁/2(2.20)i 1kB、U分

54、别为Boltzmann常数和第i个原子的速度，N为原子个数，Nc为约 束自由度数。温度控制方法通常有速度标度法、Berendsen热浴法和NosdHoover 控温方法等。(1)速度标度法55由式2.20可知温度和速度存在着一定的对应关系，故可通过对速度进行标 度进而达到控制温度的目的。若速度乘以标度因子入后，t时刻的温度的变化为(2.21)T ( 2 1)T(t)后(2'2)Treq，入分别为期望的参考温度和速度标度因子，给每一步乘上标度因子入来达到控制温度的目的。(2) Berendsen 热浴法56Bere ndsen热浴法是假想系统和一热浴相接触。对速度每一步都进行标度， 标度

55、因子为(2.23)式中，？t，tT分别为时间步长和耦合参数，Tbath, T(t)分别为t时刻系统的期望温 度和瞬时温度。(3) Nose-Hoover 控温方法57-58Nose-Hoover控温方法中，一恒定温度的热浴与体系相耦合。242 控压方法在NPT, NPH系综下，我们需要模拟定压情况下的分子动力学问题，下面主要介绍Berendsen压浴法和 Andersen方法。(1) Berendsen 压浴法59Bere ndsen压浴法是通过对体积乘以标度因子Cp，来通过控制系统体积来实 现对压力的控制。如同让一压力为期望值的压浴与系统接触。Cp 1 k-t Fbath P(t)(2.24

56、)tp式中，P(t)和Pbath分别为系统在t时刻的瞬时压力和期望压力；k和tP为耦合参数。 Andersen 方法60Andersen方法类似于图2.2中的活塞调节压力，通过活塞调节体积变化来控制压力。图2.2活塞法调节压力2.5分析方法2.5.1 中心对称参数(centrosymmetry parameter )为了可视化晶体缺陷，特别是位错的形核和发射过程，我们需要在晶体中可靠的确定这种缺陷结构的位置。基于不同的标准有许多不同的方法来研究晶 体缺陷，如原子的势能、配位数，以及中心对称参数等。本文采用Kelchner等人61介绍的中心对称参数方法。在固态物质中，中心对称参数对于区分原子为 完整晶体，局部缺陷(位错或堆垛层错)或表面原子非常有用。其定义为：PR R(2.25)i 1,这里a为近邻原子对数。对于 BCC结构，a等于4; FC