锆级联的分子动力学模拟.

1、锆辐照级联的分子动力学模拟研究小组成员：韩维强（主持人）、司宇、张锦文指导老师：侯怀宇 2012.32013.91 绪论1.1研究背景和意义能源,信息,材料,通常被称为新科技革命的三大支柱。人类文明的进程和能源的发展革命息息相关。自工业革命以来，以化石燃料作为主要能源的时代已经持续了两百多年。当今世界，人类对能源的需求与日俱增，而地球上探明的可使用的化石燃料如煤炭，石油，天然气和其他可燃气体迅速减少，因而寻找一种可替代化石能源的其他能源是迫切而必然的选择。半个多世纪以来，核电，火电和风力发电一直作为当今世界电力的三大支柱发挥着不可代替的作用。核能，以其高效，清洁，安全的优点备受青睐。在我国，目

2、前已建成秦山核电站、广东大亚核电站、江苏田湾核电站和岭澳核电站四座核电站，为我国电力事业发挥巨大作用1。核材料应具备以下条件：a. 符合核条件。材料的吸收中子截面要小；杂质含量满足核级纯；低活性，感生放射性小，剩余发热量小，半衰期短。b. 运行性能。应具备良好的物理性能；高抗辐射性能；与冷却剂相容性好。c. 经济适用性。材料具备工业生产能力，可批量生产；低成本；焊接性好。基于各项性能考虑，工业应用的核材料主要有铝及铝合金、镁合金、不锈钢、锆合金。锆合金比不锈钢的熔点高300-400,热胀系数小2/3,热导率高18%,热中子吸收截面小一个量级;力学性能、加工性能好,同UO2相容性好,尤其是对30

3、0-400的高温水、高温蒸汽也具有良好的抗腐蚀性能和足够的热强性。基于此,锆被选作反应堆燃料。元件的包壳材料和结构材料,如压力管、容器管和定位格架等,是“核安全最坚韧防线”。核反应堆的必要条件除了等离子物理外，最为重要的就是壁材料。所以说，核反应堆的诞生基于锆的发展、成熟2。1.2国内外辐照级联模拟研究现状目前，辐照损伤的模拟研究方法主要有两种：实验模拟和计算机模拟。在实验研究方面，高能电子辐照3、中子辐照4、离子辐照5-7和等离子辐照8等，是国内外模拟材料在核反应堆内的中子辐照，对核材料结构和性质的探究，欲发现其规律时常用的辐照模拟技术。上述模拟方法各有特点，如离子辐照具有实验条件可控等独特

4、优势、电子辐照不会引入杂质原子,因而可以利用它们的优点作综合分析研究。在核材料的研究领域中,计算机模拟技术得到广泛的应用。相对实验方法来说,计算机模拟方法实现起来比较容易,实验费用较低,且能够获得辐照条件下材料变化的微观细节,得到许多宏观实验中无法获得的信息,辅以相应的辐照模拟实验及理论分析,就能对辐照损伤进行更为深入细致的了解。辐照损伤的计算机模拟的主要电对象是核反应材料（如SiC，不锈钢）和航天器工程材料（如石英玻璃）。1.3本工作研究内容本文针对核结构材料的辐照损伤过程,以金属Zr为研究对象,采用分子动力学(MD)方法对辐照损伤初始过程级联碰撞进行计算机模拟,展现了辐照损伤在原子尺度的演

5、变过程,并对此过程的温度影响因素进行了研究。本文利用分子动力学方法，模拟计算了原子尺度下密排六方结构金属Zr(组成核结构材料Zr合金的主要元素)的辐照损伤过程,对其缺陷的产生演变机制进行了研究,重点讨论了PKA入射能量对材料辐照损伤过程的影响。在一定温度下，本论文分别模拟了具有不同能量的初级离位原子PKA（1Kev,5Kev,10Kev）入射时的辐照过程。模拟发现PKA能量增大时,辐照引起的级联碰撞过程更为剧烈,这说明更多原子得到了超过离位闭值Ed的能量,从而克服束缚力并从晶体点阵移出。2 分子动力学（MD）模拟方法 2.1 分子动力学模拟的概念分子动力学模拟(moleculardynamic

6、s simulation , MD)是在评估和预测材料结构和性质方面模拟原子和分子的一种物质微观领域的重要模拟方法,通过计算机对原子核和电子所构成的多体体系中的微观粒子之间相互作用和运动进行模拟,在此期间把每一原子核视为在全部其他的原子核和电子所构成的经验势场的作用下按照牛顿定律进行运动,进而得到体系中粒子的运动轨迹,再按照统计物理的方法计算得出物质的结构和性质等宏观性能.。简而言之即是应用力场及根据牛顿运动力学原理所发展的一种计算机模拟方法。经典的分子动力学方法是由Alder9和wainwright10在1957年和1959年提出并应用于理想“硬球”模型的,他们发现了早在1939年就由Kir

7、kwood根据统计力学预言的“刚性球组成的集合系统会发生由其液相到结晶相的相转变”，后来人们称这种相变为Alder相变,此举开创了用分子动力学模拟方法研究物质宏观性质的先例。用分子动力学方法，有三点决定了其结果的可靠性：势函数的优劣边界条件的合力配置以及计算机的运算速度。势函数描述了原子间的相互作用，边界条件反映了模拟对象所处的真实环境，还有高性能的计算机是实现分子动力学模拟的必要保证。2.2 分子动力学方法的原理分子动力学将连续介质看成是由N个原子或分子组成的粒子系统,体系中每个粒子都遵循牛顿第二运动定律。MD模拟中体系随着时间的演变是通过对粒子的运动方程积分实现的"通过数值求解得

8、到粒子在相空间的运动轨迹,然后由统计物理学原理得出该系统相应的宏观动态和静态特性。MD模拟是用来计算经典多体体系的平衡和传递性质的一种确定性方法。经典是指体系中的粒子的运动都遵从经典力学定律。简单的说,就是MD中处理的多体体系的粒子的运动遵从牛顿方程，即 Fi(t)=miai(t) (2.1)式中，Fi(t)是粒子i所受的力，mi是粒子的质量，ai(t)是原子i的加速度。原子i所受的力Fi(t)可以直接用势函数对坐标ri的一阶导数，即Fi(t)=-dU/dri,其中U为势函数。不同于蒙特卡洛方法,分子动力学是确定性的方法,即,原子的初始位置和速度给出后,体系内任意时刻的坐标及速度都可以通过数值

9、求解运动方程来确定。D整个运行过程中的坐标和速度称为轨迹(trajeetory),数值解普通微分方程的标准方法为有限差分法。2.3周期性边界条件及作用势2.3.1 周期性边界条件采用分子动力学对材料进行计算机模拟，必须首先对被计算的粒子系综给定适当边界条件。这些条件一般分为四种：1）自由边界条件。常用于大型的自由分子的模拟。2）固定边界。在所有要计算到的粒子晶胞之外包裹几层结构相同的、位置固定的粒子,包层厚度大于粒子间相互作用力程范围。包层部分代表了与运动粒子起作用的宏观晶体的那一部分。这种边界条件常用于对点缺陷等性质的研究。3）柔性边界。它允许边界上粒子有微小移动以反映内层粒子对其施加作用力

10、的情况。常用于模拟缺陷延伸情况。4）周期性边界条件模拟较大系统时,为了消除尺寸效应,常采用此类边界条件。它使得利用较少粒子来模拟各种系综成为可能。2.3.2 原子间作用势在MD模拟中，作用势函数的选取是个重头戏，它对模拟起着决定性的作用。势函数是描述原子间作用的函数，原子间的力场相互作用控制着原子的行为，根本上是决定了材料的性质。通俗的说，势函数是这样一个模型：作用在粒子上的力F随着原子位置的变化而变化。一般而言，主要有对势、多体势。对势是指仅有两个原子的坐标决定的相互作用，这是仅考虑到原子两两作用，而忽略其他原子对它们的作用，半导体，金属以外其他的无机化合物的相互作用可以用这样的势函数较充分

11、地描述。多体势的研究对象通常是相互作用较为复杂的多粒子系综，也是实际研究的大多数情况。2.4 牛顿方程求解算法 我们需要用积分方法求解，这种情况下，牛顿运动方程可用有限差分方法来求解。有限差分方法的思想是将积分分成很多小步，每一小步固定时间是t，在时间t时刻,作用在每个粒子的力的总和等于它与其它所有粒子的相互作用力的矢量和。根据此力,我们可以得到此粒子的加速度,结合t时刻的位置与速度,可以得到t+t时刻的位置和速度。这力在此时间间隔期间假定为常数。作用在新位置上的粒子的力可以求出,然后可以导出t+2t时刻的位置与速度等,以此类似。MD模拟中,一个好的积分算法需要具有以下条件:1,计算速度快;2

12、,需要较小的计算机内存;3,允许使用较长的时间步长;4,表现出较好的能量守恒。计算中会涉及到截断误差(truncation error)和计算机的舍入误差(round-of error).前者与选取的积分算法以及时间步长有关,后者与计算机如何实现这个算法以及存储的精度有关.通常,可以通过减小步长t来减少这两个误差。在模拟中,应该根据问题所需要的精度来选取合适的时间步长。常见的算法有Verlet算法，跳蛙算法，Gear预测-校正法等等。积分步长t较大时，用Verlet算法比较合适，积分步长t较小时，使用预估-校正法误差比较小。本文使用的就是Gear预测-校正算法。接下来，介绍一下这两种算法。2.

13、4.1 Verlet 算法 Verlet算法于1967年被提出"在分子动力学中,积分运动方程运用最广泛的方法是Verlet算法。这种算法运用t时刻的位置和速度及t-t时刻的位置,计算出t+t时刻的位置r(t+t)。 verlet算法的推导可以通过对原子位置r(t)进行Taylor展开到三阶: (2.1) (2.2)式中,v是速度,a是加速度,b是位置r对时间t的三阶导数。两式相加,则得出 (2.3)这是verlet算法最基本的形式。式(2.2)是差分公式,其中(t)的奇次项全部被抵消,同时忽略(t4)及其高次项,差分误差为O(t4)的量级。在模拟中只需要提供原子当前时刻t及前一时刻t

14、-t的位置(当前时刻下a= F/m)就可以得到下一时刻t+t的位置.速度并未出现在verlet算法中,一个简单的方法是将2.1及2.2两式相减,即用t+t时刻与t-t时刻的位置差除以2t,即可得到速度的计算公式 (2.4)速度的误差与时间步长t的二次项相关,误差在O(t3)的量级,而位置的误差仅与四次项相关。Verlet算法虽然简单且存储要求适度,但有一些缺点:计算过程容易造成精度损失;方程中无显式速度项,则在下一步位置未确定前,速度项难以求出;非自启动算法,新位置需由t和t-t的位置得到。针对Verlet算法的这些缺点,人们对Verlet算法做了一些改进,提出了新的算法,例如对速度的处理更好

15、的办法还有速度Verlet算法和蛙跳算法。2.4.2 Gear预估-校正算法Gear预估-校正是C.W.Gear在1971年提出的,是一类算法的总称。其核心思想是通过当前时刻的位置及位置对时间的前几阶导数（一阶为速度，二阶为加速度）开始，推测下一个时刻的位置及速度。通过预估位置得到的预估加速度与直接预估的加速度之间存在一个偏差；然后通过加速度的差值来校正位置和速度等。在校正的步骤里,需要事先确定对应好不同量的校正参数,因此即使对于同阶的算法,不同参数的选取会产生不同的具体算法。这种方法可分为三步： （2.5）式中,v是速度(位置对时间的一阶导数), a是加速度(位置对时间的二阶导数),b是坐标

16、对时间的三阶导数等。第二步根据新预测的位置rip(t+t),计算t+t时刻的力F(t+t),然后计算加速度aic（t+t）。这加速度再与由Taylor级数展开式预测的加速度aip(t+t)进行对比,根据两者之差在校正步里用来校正位置与速度项。通过这种校正方法,可以估计预测的加速度误差为 (2.6)假定预测的量与校正后相差很小,我们可以说,它们互相成正比,这样校正后的量为 (2.7)Gear确定了一系列的系数c0,c1,.,展开式在三阶微分b(t)后被截断。采用的系数的近似值为c0=1/6,cl=5/6,c2=1和c3=1/3。Gear的预估一校正算法需要的存储量为3(O+l)N,O是应用的最高

17、阶微分数,N是原子数目。2.5 MD模拟的粒子系综系综（Ensemble）即系统的综合，它是统计力学中的一个概念，是由1901年由吉布斯创立完成的。在进行分子动力学模拟中，需要根据实际条件选定一定的系综。分子动力学中运用系综，是为了研究系统的微观状态与宏观性质对应的规律。分子动力学的研究对象是多粒子系综，其模拟的粒子数目有限，而且受到计算机性能的限制，但统计物理规律仍成立。物质世界是由聚集在一起的大量粒子构成,而宏观物质的行为是由微观粒子间的相互作用决定的。为了研究宏观体系的演化和结构,基本方法就是构建一个粒子系综。系综是一个巨大的系统,它是指在一定的宏观条件下,大量组成、性质、尺寸和形状完全

18、相同的、处于各种运动状态的、各自独立的系统的集合。其中的每个系统各处在某一微观运动状态,而且各自独立。微观运动状态在相空间中构成一个连续区域,与微观量相对应的宏观量是在一定的宏观条件下所有可能的运动状态的平均值。MD模拟包括平衡态和非平衡态模拟。主要系综包括：1）能量和粒子数都固定的系综（N,V,E）：体系与外界不能交换能量，系统原子数N，体积V，能量E都保持恒定。一般而言，给定能量的精确初始条件无法知晓，为了把系统调节到给定的能量，需要先给出一个合理的初始条件，然后对能量进行调整，使系统达到所要状态。一般是对速度V进行特别标度来调整能量，但这样会使系统偏离平衡，因此要给足时间来重新让系统建立

19、平衡。2）原子数N，体积V，和温度T都保持不变，总能量为零的系统（N,V,T）：让系统与虚拟热裕保持热平衡状态来控制温度不变。由于温度和动能有直接关系，所以要固定系统的动能。3）原子数N，温度T，压力P都恒定的等温等压系综：温度的恒定是通过调节系统速度或加一约束力来实现的。P，V是相关量，是通过标定体积V来调节压力值的。2.6 MD运行基本步骤 虽然分子动力学核心思想比较简单,但实际上分子动力学的实现在很多方面都面临挑战,如模型的选择、初始条件的设定,为保证可靠性而进行的各种模拟方案,积分算法的选择,需要考虑到运动轨迹对初始条件及其他参数选择的敏感性,满足大计算量的要求,图形显示与数据分析等等

20、11。分子动力学的模拟流程大致如下：1.选取势函数，即是建立一个计算模型。势函数与物质描述密切相关。当势函数模型确定时，可以根据物理学规律得到模拟中的守恒关系。2.初始化坐标和速度。MD过程中对运动方程的求解，必须知道每个粒子初始位置和速度。（不同的积分方法要求不同的初始条件。）3.趋于平衡的计算。MD模拟中需要设计一个使系统达到一个稳定的平衡状态的过程。当系统的动能和势能的总能量在某个值附近波动时，认为已经达到平衡。这段达到平衡状态的时间成为弛豫时间。4.最后针对模拟对象进行物理量的计算。参与计算的物理量是在模拟的最后进行的。它是沿相空间轨迹求平均得到的。 2.7本文MD模拟采用的软件情况

21、本文分子动力学模拟采用的软件是 Lammps。Lammps是分子动力学模拟常用的一款软件，其源代码公开，免费下载，可以根据自己的需要修改lammps代码，重新编译；lammps需要在linux环境下运行，可以串行和并行运算，高移植性的C+语言编写；lammps是美国能源部下属的圣地亚国家实验室发布的，主要作者是Steve Plimpton。2.7.1 Lammps软件优点可以在Linux并行环境中运行，体现出强大的计算优势；能模拟上百万的原子体系，气态、液态或者固态，在各种系综下；提供了各种势函数可供选择（只是提供了各个函数的表达式，具体的参数是自己找的）。2.7.2 Lammps 运行简要步

22、骤从文件读取，可包括拓扑信息构造晶胞复制命令设置力场参数模拟参数输出文件。3 锆辐照损伤的MD模拟分子动力学的方法虽然能很好地模拟在极端情况下的辐照现象，但是受到很多方面的影响。如果进行一定程度的修改，会使得结果更贴近真实情况。3.1模拟方法3.1.1 本文所用势函数 本文采用Mendelev在2007年提出的Zr的EAM势函数12。使用的是文献13中的2号势函数。EAM势函数的总能量：， （3.1）其下标i，j表示N原子系统中的每个原子， 是成对可能性，是嵌入能量函数。其中， （3.2）是另外一个成对可能：“密度函数”为建立一个理想的EAM势函数需要给寻找最佳函数。 本函数正确描述了HCP到

23、bcc相的过度和液体结构数据。作为第一个测试，其可能性被应用在研究HCP和bcc的自扩散，发现了异常。伴随自我间迁移的高扩散性与bcc内在结构缺陷的形成有关。 3.1.2可变时间步长辐照条件下，原子间的距离会变得非常小（1À甚至更小），从而原子间的作用力会非常大，标准的势函数很难描述原子间的作用势，单一的势函数会使得误差增大。因而需要设置这样一种机制，即当原子受到高能量作用时，时间步长自动缩小，当原子受到低能量作用时，时间步长会自动增大。这样的机制既保证了准确，又最大限度地控制了模拟的运行时间。在辐照模拟中，常用的算法是依据计算出的时间步长，对原子在每一个MD时间步长内移动的距离加一

24、个限制。一般会给定一个时间步长（常用1fs,0.5fs）在模拟中，我们采用可变时间步长的方法，步长不断调整。设定一个最小时间步长，最大时间步长为。每一个时间步长内移动的距离限制在一定范围，最大移动范围设定为，使得所有原子基于当前原子的速度和力，移动的距离都不超过最大距离。当从一个原子团簇开始模拟或运行一个或更多原子撞击固体引起级联损伤反应的模拟时可能都是有用的。每步之后，下一步的时间步长dt按下列方式计算：对于每个原子，时间步长计算，将导致它的位移在接下来的集成步骤，为其当前的速度和力的作用。由于执行此精密计算需要解决一个四次方程，产生一个更简单的估测。估测是保守的，尽管它可能较小，但原子的位

25、移保证不会超过。3.1.3 缺陷辨别方法及确定空位和离位原子122 A B 图3-1 完整晶体结构与有缺陷的晶体结构 在辐照模拟分析之前，需要说明缺陷的判别方法。缺陷判别方法，就是对比完整晶体结构和有缺陷的晶体结构，根据两结构的不同来寻找缺陷。如图3-1，完整晶体结构为A，有缺陷的晶体结构为B。B中白圆圈表示空位，右下角的黑圆圈表示间隙原子。在B中的某一原子如果不能在A中找到对应的原子，则认为它是间隙原子。对格点进行循环检测：1.每个格子点，确定其所在晶胞的坐标（x，y，z）；2.确定其上下左右的9个晶胞，对每个晶胞内的原子进行原子-格点的距离计算；根据，确定格点的占据原子（可以多于一个）并确

26、定离格点最近的原子，3.如果离格点最近的原子距离大于，则该点是空位；4.格点占据者中，并非最近原子的，就是与该点相联系的间隙原子。上法已确定一部分间隙原子，若原子并非任何一个格点的占据者，那它也是间隙原子。3.1.4 PKA速度的确定本文使用MD对锆中辐照损伤过程进行了研究，模拟中，选定一原子赋予一定的初速度来模拟辐照过程。原子被赋予的速度与PKA能量的转换关系为： E=mv2 (3.3)式中，E为能量（单位：J），m为原子质量（单位：kg），v为原子速度（单位：m/s）。本模拟实验中，原子速度v单位为À/ps，单位换算涉及常量有1ev=1.6X10-19J;C12=1.66X10-

27、27kg;1À/ps=100m/s.3.2模拟设置3.2.1 模拟体系范围 在模拟过程中,虽然缺陷不直接与热浴区域发生相互作用,但模拟体系中间活动区域的大小对缺陷的产生也有影响,这主要是因为系统的能量密度是由PKA能量的大小以及活动区域的体积来控制的13。3.2.2 平衡时间MD模拟的优点在于计算量小、计算时间短。模拟过程的计算时间长短由平衡时间来决定。平衡时间是分子动力学模拟中一个重要的参数,如果过大则会增加不必要的计算量；如果过小,模拟体系未达到平衡就输出模拟结果,使模拟情况不能反映真实情况。所以选择合适的平衡时间对于计算时间节省和,实验结果的准确性是重要的。本文对辐照模拟的实验

28、时间均>10ps。模拟体系的大小对平衡平衡时间也有影响。在MD模拟中,原子不断的运动,则原子与原子的之间的相互作用力会不断的改变,这就是能量产生波动的原因。随着体系的变小,影响原子运动的条件增多,所以原子达到平衡态所需要的时间也就越长。3.2.3 弛豫时间 模拟中，需要初始化每个原子的初始速度和坐标，但初始化的矢量不一定就是对应势函数的最小值，需要最小化来弛豫应力。体系的初始坐标和初始速度设置以后,在进行模拟体系的性质以前,首先必须使体系进行趋于平衡的过程。在这个过程中,动能、势能相互转化,当动能、势能、总能量只在某平均值附近波动、趋于平稳时,就可认为体系就达到了平衡。3.3模拟过程和结

29、果本文以HCP密排六方晶体结构的金属锆（Zr）为研究对象，使用分子动力学方法对其级联碰撞的过程进行了模拟研究。首先对锆的结构和性质进行简要介绍，再结合辐照损伤理论和MD模拟方法，对锆级联过程的PKA入射能量的影响进行了模拟实验，并对结果进行了分析。3.3.1 Zr的基本性质在1789年，德国化学家Martin Klaproth分析了一块锆石，并分离出了锆的“泥土”形式氧化锆，也就是氧化物ZrO2。随着错合金在原子能工业上的应用,错工业有了迅速的发展，成为一种重要的战略材料,被誉为“原子时代的第一金属”14。Zr在化学元素周期表中位于第5周期,其位置在IVB族的钦(Ti)之下,原子序数为40,平

30、均相对原子质量为91.22。纯锆在室温下具有HCP密排六方的晶体结构,称为-Zr晶格常数为a=0.323nm,c=0.515nm。a-Zr在865时发生同素异晶转变,转变为体心立方-Zr,相变时会有较大的体积收缩。锆是钛的同族元素，也是同素异构体。862以下呈密排六方晶格，称为锆。862以上呈体心立方晶格，称为锆。纯锆的熔点Tm是1850左右。a-Zr的结构如图3-2所示。图3-2 -Zr的晶格结构由于核性能优异、力学性能适中和加工性能良好,因而错已被普遍用作核动力水冷反应堆的燃料包壳管和结构材料,如压力管、容器管、孔道馆、导向管、定位格架、端塞和其他结构文件,这是错材的主要用途,占整个错加工

31、材的80%15在改善锆合金机械性能和抗腐蚀性能，研究人员已经开发出了N18,N36等新锆合金，其焊接性能和力学性能也有所提高。3.3.2 级联碰撞的演化过程 100K温度下金属Zr中能量为1KeV的PKA原子入射模拟的结果表明：Zr的级联碰撞过程非常短暂，大约是ps量级（1ps=10-12s）。随着时间的推移，一开始缺陷迅速增多，随后渐渐下降，趋于平缓，最后达到一个稳定状态。这是因为，当高于Zr原子离位阈值Ed的粒子入射时，晶体内会产生级联碰撞过程，产生间隙原子和空位并存的缺陷。在碰撞末期，所有原子的能量都降到离位阈值以下，不再撞击晶格原子。之后，通过间隙原子的弹性复合，出现“急增-衰减-平缓

32、”的现象。利用lammps软件对辐照损伤过程进行计算机计算，以便更直接地得到级联碰撞过程的演化情况。观察的结果也和缺陷产生规律一致。当PKA原子能量增大时，即辐照增强时，级联碰撞将更为剧烈。这样，将有更多的原子被赋予大于离位阈值的能量，从而有更多的原子从点阵移出。本文模拟了PKA能量分别是1KeV、5 KeV、10 KeV入射Zr基体。结果发现，随着PKA能量的增大，Zr基体内产生的缺陷数目越多，这说明PKA能量的增加会使碰撞加剧，即辐照损伤加剧。3.3.3 MD元胞设置如下图所示，a-Zr为HCP结构，每个晶胞中有两个原子，位于（0，0，0）和（1/3，2/3，1/2）处。把HCP晶格变成如

33、图中红线框所示的长方体晶格, 每个晶格中有4个原子，位于（0,0,0）、（0.5,5/6,0.5）、（0,1/3,0.5）和（0.5,0.5,0）。设原始坐标下晶格参数为a0和c0，则本文所采用的长方体晶格尺寸为： a0=0.3231nm, c0=5.148nm。a= a0,b=a0c=c0图3-3 HCP结构平面图3.3.4 模拟步骤（1）在温度设定为600K且不变的情况下，设定一个大小为203020总原子个数为16000个原子的的长方体元胞。初始晶胞大小为a0=0.3231nm, c0=5.148nm。在600K和零外压下运行30000个时间步长后，得到平衡状态下的元胞体积和晶胞的参数，结

34、果为 a = 0.3232 nm b = 0.5598 nm c = 0.5212 nm(2)在不同的PKA能量下，选取不同大小的MD盒子，其选取的MD盒子表3-1:表3-1 MD盒子设计PKA能量 MD元胞大小 总原子个数 1KeV 301520 36000 5KeV 502530 150000 10KeV 603040 288000先在NV条件下运行30000步，其时间步长为1fs.。运行结束经过驰豫过程后得到无缺陷状态下的平衡构型（3）选取MD胞中的一个原子，赋予它不同的初始动能（由 E=mv2 可知给予原子一定的能量它就拥有一定的初始速度）,速度的方向随机(可采用01中的随机数得到设定

35、两个随机数sin，sin ,从而确定速度的方向)。在NVE下运行，采用变时间步长法，在不同能量下重复多次，最终得到平均后的实验结果。能量和次数如表3-2所示:表3-2 不同PKA能量下的重复次数 PKA能量 重复次数 1KeV 5 5KeV 10 10KeV 153.4结果分析与讨论本文以密排六方晶体结构的Zr为研究对象，使用分子动力学计算模拟方法，对辐照级联进行了模拟研究。通过计算，对模拟实施过程中的时间步长、热裕、缺陷判别等方面进行了测试，建立了一套模拟参数。利用分子动力学计算软件lammps对级联过程进行了观察。3.4.1 模拟过程中的温升效应 NVE条件下，系综的粒子数和总能量保持恒定

36、，系综与外界没有能量交换。模拟设定的温度是600K。当PKA原子获得模拟给予的初始能量，原子即以高速度发射，获得动能。撞击过程必然导致其他晶格原子的移位，造成混乱，同时由于撞击PKA原子动能衰减。因为能量的转化，伴随着PKA原子动能的下降，原子的势能会上升，但动能并非全部转化为势能，一部分机械能的损失转化成了热量。微观世界里的剧烈撞击，短时间内转化了大量热能，这样体系温度迅速攀升。体系的平均温度可达到800K。由于PKA原子的撞击，离位原子产生了速度，离位原子的动能提高，再者，一部分动能的贡献来自于高温，即800K体系热量的传送。这样，体系温度又开始下降，最后，温度逐渐平稳下来。所以，整个过程

37、呈现一个温度快速上升继而逐渐下降最后趋于稳定状态。碰撞后体系中温度随时间的变化关系如图3-43-6。图3-4（1Kev）温度随时间的变化图3-5（5Kev）温度随时间的变化图3-6（10Kev）温度随时间的变化原子碰撞后动能转换为势能，动能减少，温度下降。后因势能弛豫，温度小幅上升，最后处于稳定状态。与此同时，动能和势能的变化如图3-73-9：图3-7 1Kev下动能和势能对时间的变化关系图图3-8 5Kev时动能和势能与时间的变化关系图3-9 10Kev下动能和势能与时间的变化关系3.4.2 级联过程中空位与离位原子在不同时刻的分布模型图3-10直观展示了整个级联过程中空位和离位原子分布的变

38、化。图中，红色圆点表示空位，绿色圆点表示离位原子，如图3-10（a）（e），在PKA能量为1Kev条件下，级联开始时，空位和离位原子都非常少，然后空位和离位原子都随着时间增多，大约进行到0.38ps时，空位和离位原子的数量都达到峰值，如图3-10（c）所示。之后空位和离位原子数量开始减少，最后达到一个稳定状态，如图3-10（e）所示。3-11和3-12的各图表现了相似的过程，不同之处在于PKA能量为5Kev和10Kev时产生空位和离位原子的峰值更大。 （a） （b） （c） （d） （e）图3-10 PKA能量为1Kev时各时刻的空位与离位原子分布(a) time=ps (b) time=0.

39、18ps (c) time=0.38ps (d) time=1.37ps (e) time=3.77ps （a） （b） （c） （d） (e)图3-11 PKA能量为5Kev时各时刻的空位与离位原子分布(a) time=ps (b) time=0.83ps (c) time=2.06ps (d) time=4.54ps (e) time=8.21ps （a） （b） （c） （d） (e)图3-12 PKA能量为10Kev时各时刻的空位与离位原子分布(a) time=ps (b) time=0.15ps (c) time=0.76ps (d) time=3.60ps (e) time=11.

40、01ps3.4.3 级联过程中产生的空位与离位原子统计基于对锆的基本性质的了解，用第2章介绍的MD模拟方法对锆级联的碰撞过程及此过程中设定不同PKA入射原子能量进行了模拟实验，并对结果进行了分析。模拟结果表明：材料在辐照过程中，随着时间的推移，晶体产生的缺陷数目一开始猛增，这是因为当PKA能量给出时，高于原子离位阈值的原子数目较多，碰撞概率很高。随后，因为间隙原子与空位的弹性复合，缺陷数又会慢慢下降，最后趋于平缓，达到一个稳定状态；随着PKA能量的增大，材料基体内级联碰撞产生的缺陷也会相应增多。本文的模拟结果较深入地了解了级联过程及其影响因素、改造材料结构对其抗腐蚀性，机械性能等有着重要作用。

41、模拟的结果以时间为横坐标，Flenkel 缺陷对数目为纵坐标，绘制了不同PKA入射原子能量下的关系图。如图3-133-15。 图3-13 time-Flenkel 缺陷对（1Kev） 图3-14 time-Flenkel 缺陷对（5Kev） 图3-15 time-Flenkel 缺陷对（10Kev）3.4.4 空位和离位原子的配位数统计级联发生开始前，体系中空位和离位原子很少，且都小于2，基本上都是离散的。说明发生级联碰撞前体系中的缺陷很少，而且没有原子团簇。大约在ps时，空位的配位数开始增大，同时有1配位的空位出现；大约在ps时，出现了2配位的空位；大约在ps时，出现了配位数大于2的空位。说

42、明此时已出现大量原子团簇。离位原子配位数的分布情况与空位相似，在此不再赘述。图3-16ae为在10Kev条件下空位配位数不同时刻的分布情况。（a）t=0.028ps, (b)t=0.061ps, (c)t=0.082ps, (d) t=0.101ps, (e)t=0.143ps(a) (b)(c)(d)(e) 图3-16 PKA能量为10Kev时各时刻的空位的分布4 总结本文利用分子动力学计算模拟方法，借助lammps计算机计算软件，以HCP密排六方Zr为研究对象，对材料辐照损伤过程进行了模拟，在原子尺度上研究了目前无法用实验的方法观察到的辐照损伤过程。重点研究了PKA入射能量对材料辐照损伤的

43、影响。此外，在本项目的研究过程中，设计了可变时间步长对温度、Flenkel缺陷对数目、空位和离位原子等变化关系的程序，一定程度上改进和发展了此程序。主要结论如下：1）PKA入射能量会对级联过程产生一定影响。在一定温度下，分别模拟了具有不同能量的PKA（1KeV、5 KeV、10 KeV）入射时的辐照过程。发现PKA能量增大时，辐照级联过程将更剧烈。2）材料辐照作用过程非常短暂，大约是ps量级；随时间的推进，缺陷数目迅速到达至高点，然后下降，最后趋于平缓；3）级联发生的一开始，由于碰撞原子对其他原子进行了碰撞，温度迅速升高。之后，动能急剧下降，离位原子与空位的逐渐产生使势能得到积累，同时温度下降，然后由于势能弛豫以及体系内部的弹性回复，部分离位原子与空位结合，势能减少，温度小幅上升，之后温度稳定在一定范围内波动。最后，温度稳定在720左右。过程中，温度呈现一种突升急降小幅上升平缓的态势。4）在PKA入射时体系内缺陷很少，随着时间的推移，