锆级联的分子动力学模拟Word文档格式.docx
《锆级联的分子动力学模拟Word文档格式.docx》由会员分享，可在线阅读，更多相关《锆级联的分子动力学模拟Word文档格式.docx（28页珍藏版）》请在冰豆网上搜索。
基于此,锆被选作反应堆燃料。
元件的包壳材料和结构材料,如压力管、容器管和定位格架等,是“核安全最坚韧防线”。
核反应堆的必要条件除了等离子物理外，最为重要的就是壁材料。
所以说，核反应堆的诞生基于锆的发展、成熟[2]。
1.2国内外辐照级联模拟研究现状
目前，辐照损伤的模拟研究方法主要有两种：
实验模拟和计算机模拟。
在实验研究方面，高能电子辐照[3]、中子辐照[4]、离子辐照[5-7]和等离子辐照[8]等，是国内外模拟材料在核反应堆内的中子辐照，对核材料结构和性质的探究，欲发现其规律时常用的辐照模拟技术。
上述模拟方法各有特点，如离子辐照具有实验条件可控等独特优势、电子辐照不会引入杂质原子,因而可以利用它们的优点作综合分析研究。
在核材料的研究领域中,计算机模拟技术得到广泛的应用。
相对实验方法来说,计算机模拟方法实现起来比较容易,实验费用较低,且能够获得辐照条件下材料变化的微观细节,得到许多宏观实验中无法获得的信息,辅以相应的辐照模拟实验及理论分析,就能对辐照损伤进行更为深入细致的了解。
辐照损伤的计算机模拟的主要电对象是核反应材料（如SiC，不锈钢）和航天器工程材料（如石英玻璃）。
1.3本工作研究内容
本文针对核结构材料的辐照损伤过程,以金属Zr为研究对象,采用分子动力学（MD）方法对辐照损伤初始过程——级联碰撞进行计算机模拟,展现了辐照损伤在原子尺度的演变过程,并对此过程的温度影响因素进行了研究。
本文利用分子动力学方法，模拟计算了原子尺度下密排六方结构金属Zr（组成核结构材料Zr合金的主要元素）的辐照损伤过程,对其缺陷的产生演变机制进行了研究,重点讨论了PKA入射能量对材料辐照损伤过程的影响。
在一定温度下，本论文分别模拟了具有不同能量的初级离位原子PKA（1Kev,5Kev,10Kev）入射时的辐照过程。
模拟发现PKA能量增大时,辐照引起的级联碰撞过程更为剧烈,这说明更多原子得到了超过离位闭值Ed的能量,从而克服束缚力并从晶体点阵移出。
2分子动力学（MD）模拟方法
2.1分子动力学模拟的概念
分子动力学模拟（moleculardynamicssimulation,MD）是在评估和预测材料结构和性质方面模拟原子和分子的一种物质微观领域的重要模拟方法,通过计算机对原子核和电子所构成的多体体系中的微观粒子之间相互作用和运动进行模拟,在此期间把每一原子核视为在全部其他的原子核和电子所构成的经验势场的作用下按照牛顿定律进行运动,进而得到体系中粒子的运动轨迹,再按照统计物理的方法计算得出物质的结构和性质等宏观性能.。
简而言之即是应用力场及根据牛顿运动力学原理所发展的一种计算机模拟方法。
经典的分子动力学方法是由Alder[9]和wainwright[10]在1957年和1959年提出并应用于理想“硬球”模型的,他们发现了早在1939年就由Kirkwood根据统计力学预言的“刚性球组成的集合系统会发生由其液相到结晶相的相转变”，后来人们称这种相变为Alder相变,此举开创了用分子动力学模拟方法研究物质宏观性质的先例。
用分子动力学方法，有三点决定了其结果的可靠性：
势函数的优劣边界条件的合力配置以及计算机的运算速度。
势函数描述了原子间的相互作用，边界条件反映了模拟对象所处的真实环境，还有高性能的计算机是实现分子动力学模拟的必要保证。
2.2分子动力学方法的原理
分子动力学将连续介质看成是由N个原子或分子组成的粒子系统,体系中每个粒子都遵循牛顿第二运动定律。
MD模拟中体系随着时间的演变是通过对粒子的运动方程积分实现的"
通过数值求解得到粒子在相空间的运动轨迹,然后由统计物理学原理得出该系统相应的宏观动态和静态特性。
MD模拟是用来计算经典多体体系的平衡和传递性质的一种确定性方法。
经典是指体系中的粒子的运动都遵从经典力学定律。
简单的说,就是MD中处理的多体体系的粒子的运动遵从牛顿方程，即
Fi（t）=miai（t）（2.1）
式中，Fi（t）是粒子i所受的力，mi是粒子的质量，ai（t）是原子i的加速度。
原子i所受的力Fi（t）可以直接用势函数对坐标ri的一阶导数，即Fi（t）=-dU/dri,其中U为势函数。
不同于蒙特卡洛方法,分子动力学是确定性的方法,即,原子的初始位置和速度给出后,体系内任意时刻的坐标及速度都可以通过数值求解运动方程来确定。
D整个运行过程中的坐标和速度称为轨迹（trajeetory）,数值解普通微分方程的标准方法为有限差分法。
2.3周期性边界条件及作用势
2.3.1周期性边界条件
采用分子动力学对材料进行计算机模拟，必须首先对被计算的粒子系综给定适当边界条件。
这些条件一般分为四种：
1）自由边界条件。
常用于大型的自由分子的模拟。
2）固定边界。
在所有要计算到的粒子晶胞之外包裹几层结构相同的、位置固定的粒子,包层厚度大于粒子间相互作用力程范围。
包层部分代表了与运动粒子起作用的宏观晶体的那一部分。
这种边界条件常用于对点缺陷等性质的研究。
3）柔性边界。
它允许边界上粒子有微小移动以反映内层粒子对其施加作用力的情况。
常用于模拟缺陷延伸情况。
4）周期性边界条件模拟较大系统时,为了消除尺寸效应,常采用此类边界条件。
它使得利用较少粒子来模拟各种系综成为可能。
2.3.2原子间作用势
在MD模拟中，作用势函数的选取是个重头戏，它对模拟起着决定性的作用。
势函数是描述原子间作用的函数，原子间的力场相互作用控制着原子的行为，根本上是决定了材料的性质。
通俗的说，势函数是这样一个模型：
作用在粒子上的力F随着原子位置的变化而变化。
一般而言，主要有对势、多体势。
对势是指仅有两个原子的坐标决定的相互作用，这是仅考虑到原子两两作用，而忽略其他原子对它们的作用，半导体，金属以外其他的无机化合物的相互作用可以用这样的势函数较充分地描述。
多体势的研究对象通常是相互作用较为复杂的多粒子系综，也是实际研究的大多数情况。
2.4牛顿方程求解算法
我们需要用积分方法求解，这种情况下，牛顿运动方程可用有限差分方法来求解。
有限差分方法的思想是将积分分成很多小步，每一小步固定时间是δt，在时间t时刻,作用在每个粒子的力的总和等于它与其它所有粒子的相互作用力的矢量和。
根据此力,我们可以得到此粒子的加速度,结合t时刻的位置与速度,可以得到t+δt时刻的位置和速度。
这力在此时间间隔期间假定为常数。
作用在新位置上的粒子的力可以求出,然后可以导出t+2δt时刻的位置与速度等,以此类似。
MD模拟中,一个好的积分算法需要具有以下条件:
1,计算速度快;
2,需要较小的计算机内存;
3,允许使用较长的时间步长;
4,表现出较好的能量守恒。
计算中会涉及到截断误差（truncationerror）和计算机的舍入误差（round-oferror）.前者与选取的积分算法以及时间步长有关,后者与计算机如何实现这个算法以及存储的精度有关.通常,可以通过减小步长δt来减少这两个误差。
在模拟中,应该根据问题所需要的精度来选取合适的时间步长。
常见的算法有Verlet算法，跳蛙算法，Gear预测-校正法等等。
积分步长δt较大时，用Verlet算法比较合适，积分步长δt较小时，使用预估-校正法误差比较小。
本文使用的就是Gear预测-校正算法。
接下来，介绍一下这两种算法。
2.4.1Verlet算法
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-δt的位置（当前时刻下a=F/m）就可以得到下一时刻t+δt的位置.速度并未出现在verlet算法中,一个简单的方法是将2.1及2.2两式相减,即用t+δt时刻与t-δt时刻的位置差除以2δt,即可得到速度的
计算公式
（2.4）
速度的误差与时间步长δt的二次项相关,误差在O（δt3）的量级,而位置的误差仅与四次项相关。
Verlet算法虽然简单且存储要求适度,但有一些缺点:
计算过程容易造成精度损失;
方程中无显式速度项,则在下一步位置未确定前,速度项难以求出;
非自启动算法,新位置需由t和t-δt的位置得到。
针对Verlet算法的这些缺点,人们对Verlet算法做了一些改进,提出了新的算法,例如对速度的处理更好的办法还有速度Verlet算法和蛙跳算法。
2.4.2Gear预估-校正算法
Gear预估-校正是C.W.Gear在1971年提出的,是一类算法的总称。
其核心思想是通过当前时刻的位置及位置对时间的前几阶导数（一阶为速度，二阶为加速度）开始，推测下一个时刻的位置及速度。
通过预估位置得到的预估加速度与直接预估的加速度之间存在一个偏差；
然后通过加速度的差值来校正位置和速度等。
在校正的步骤里,需要事先确定对应好不同量的校正参数,因此即使对于同阶的算法,不同参数的选取会产生不同的具体算法。
这种方法可分为三步：
（2.5）
式中,v是速度（位置对时间的一阶导数）,a是加速度（位置对时间的二阶导数）,b是坐标对时间的三阶导数等。
第二步根据新预测的位置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是应用的最高阶微分数,N是原子数目。
2.5MD模拟的粒子系综
系综（Ensemble）即系统的综合，它是统计力学中的一个概念，是由1901年由吉布斯创立完成的。
在进行分子动力学模拟中，需要根据实际条件选定一定的系综。
分子动力学中运用系综，是为了研究系统的微观状态与宏观性质对应的规律。
分子动力学的研究对象是多粒子系综，其模拟的粒子数目有限，而且受到计算机性能的限制，但统计物理规律仍成立。
物质世界是由聚集在一起的大量粒子构成,而宏观物质的行为是由微观粒子间的相互作用决定的。
为了研究宏观体系的演化和结构,基本方法就是构建一个粒子系综。
系综是一个巨大的系统,它是指在一定的宏观条件下,大量组成、性质、尺寸和形状完全相同的、处于各种运动状态的、各自独立的系统的集合。
其中的每个系统各处在某一微观运动状态,而且各自独立。
微观运动状态在相空间中构成一个连续区域,与微观量相对应的宏观量是在一定的宏观条件下所有可能的运动状态的平均值。
MD模拟包括平衡态和非平衡态模拟。
主要系综包括：
1）能量和粒子数都固定的系综（N,V,E）：
体系与外界不能交换能量，系统原子数N，体积V，能量E都保持恒定。
一般而言，给定能量的精确初始条件无法知晓，为了把系统调节到给定的能量，需要先给出一个合理的初始条件，然后对能量进行调整，使系统达到所要状态。
一般是对速度V进行特别标度来调整能量，但这样会使系统偏离平衡，因此要给足时间来重新让系统建立平衡。
2）原子数N，体积V，和温度T都保持不变，总能量为零的系统（N,V,T）：
让系统与虚拟热裕保持热平衡状态来控制温度不变。
由于温度和动能有直接关系，所以要固定系统的动能。
3）原子数N，温度T，压力P都恒定的等温等压系综：
温度的恒定是通过调节系统速度或加一约束力来实现的。
P，V是相关量，是通过标定体积V来调节压力值的。
2.6MD运行基本步骤
虽然分子动力学核心思想比较简单,但实际上分子动力学的实现在很多方面都面临挑战,如模型的选择、初始条件的设定,为保证可靠性而进行的各种模拟方案,积分算法的选择,需要考虑到运动轨迹对初始条件及其他参数选择的敏感性,满足大计算量的要求,图形显示与数据分析等等[11]。
分子动力学的模拟流程大致如下：
1.选取势函数，即是建立一个计算模型。
势函数与物质描述密切相关。
当势函数模型确定时，可以根据物理学规律得到模拟中的守恒关系。
2.初始化坐标和速度。
MD过程中对运动方程的求解，必须知道每个粒子初始位置和速度。
（不同的积分方法要求不同的初始条件。
）
3.趋于平衡的计算。
MD模拟中需要设计一个使系统达到一个稳定的平衡状态的过程。
当系统的动能和势能的总能量在某个值附近波动时，认为已经达到平衡。
这段达到平衡状态的时间成为弛豫时间。
4.最后针对模拟对象进行物理量的计算。
参与计算的物理量是在模拟的最后进行的。
它是沿相空间轨迹求平均得到的。
2.7本文MD模拟采用的软件情况
本文分子动力学模拟采用的软件是Lammps。
Lammps是分子动力学模拟常用的一款软件，其源代码公开，免费下载，可以根据自己的需要修改lammps代码，重新编译；
lammps需要在linux环境下运行，可以串行和并行运算，高移植性的C++语言编写；
lammps是美国能源部下属的圣地亚国家实验室发布的，主要作者是StevePlimpton。
2.7.1Lammps软件优点
可以在Linux并行环境中运行，体现出强大的计算优势；
能模拟上百万的原子体系，气态、液态或者固态，在各种系综下；
提供了各种势函数可供选择（只是提供了各个函数的表达式，具体的参数是自己找的）。
2.7.2Lammps运行简要步骤
从文件读取，可包括拓扑信息→构造晶胞→复制命令→设置力场参数→模拟参数→输出文件。
3锆辐照损伤的MD模拟
分子动力学的方法虽然能很好地模拟在极端情况下的辐照现象，但是受到很多方面的影响。
如果进行一定程度的修改，会使得结果更贴近真实情况。
3.1模拟方法
3.1.1本文所用势函数
本文采用Mendelev在2007年提出的Zr的EAM势函数[12]。
使用的是文献[13]中的2号势函数。
EAM势函数的总能量：
，（3.1）
其下标i，j表示N原子系统中的每个原子，
是成对可能性，
是嵌入能量函数。
其中
，（3.2）
是另外一个成对可能：
“密度函数”为建立一个理想的EAM势函数需要给
寻找最佳函数。
本函数正确描述了HCP到bcc相的过度和液体结构数据。
作为第一个测试，其可能性被应用在研究HCP和bcc的自扩散，发现了异常。
伴随自我间迁移的高扩散性与bcc内在结构缺陷的形成有关。
3.1.2可变时间步长
辐照条件下，原子间的距离会变得非常小（~1À
甚至更小），从而原子间的作用力会非常大，标准的势函数很难描述原子间的作用势，单一的势函数会使得误差增大。
因而需要设置这样一种机制，即当原子受到高能量作用时，时间步长自动缩小，当原子受到低能量作用时，时间步长会自动增大。
这样的机制既保证了准确，又最大限度地控制了模拟的运行时间。
在辐照模拟中，常用的算法是依据计算出的时间步长，对原子在每一个MD时间步长内移动的距离加一个限制。
一般会给定一个时间步长（常用1fs,0.5fs…）在模拟中，我们采用可变时间步长的方法，步长不断调整。
设定一个最小时间步长
，最大时间步长为
。
每一个时间步长内移动的距离限制在一定范围，最大移动范围设定为
，使得所有原子基于当前原子的速度和力，移动的距离都不超过最大距离。
当从一个原子团簇开始模拟或运行一个或更多原子撞击固体引起级联损伤反应的模拟时可能都是有用的。
每步之后，下一步的时间步长dt按下列方式计算：
对于每个原子，时间步长计算，将导致它的位移
在接下来的集成步骤，为其当前的速度和力的作用。
由于执行此精密计算需要解决一个四次方程，产生一个更简单的估测。
估测是保守的，尽管它可能较小，但原子的位移保证不会
超过
3.1.3缺陷辨别方法及确定空位和离位原子
AB
图3-1完整晶体结构与有缺陷的晶体结构
在辐照模拟分析之前，需要说明缺陷的判别方法。
缺陷判别方法，就是对比完整晶体结构和有缺陷的晶体结构，根据两结构的不同来寻找缺陷。
如图3-1，完整晶体结构为A，有缺陷的晶体结构为B。
B中白圆圈表示空位，右下角的黑圆圈表示间隙原子。
在B中的某一原子如果不能在A中找到对应的原子，则认为它是间隙原子。
对格点进行循环检测：
1.每个格子点，确定其所在晶胞的坐标（x，y，z）；
2.确定其上下左右的9个晶胞，对每个晶胞内的原子进行原子-格点的距离计算；
根据
，确定格点的占据原子（可以多于一个）并确定离格点最近的原子，
3.如果离格点最近的原子距离大于
，则该点是空位；
4.格点占据者中，并非最近原子的，就是与该点相联系的间隙原子。
上法已确定一部分间隙原子，若原子并非任何一个格点的占据者，那它也是间隙原子。
3.1.4PKA速度的确定
本文使用MD对锆中辐照损伤过程进行了研究，模拟中，选定一原子赋予一定的初速度来模拟辐照过程。
原子被赋予的速度与PKA能量的转换关系为：
E=
mv2（3.3）
式中，
E为能量（单位：
J），
m为原子质量（单位：
kg），
v为原子速度（单位：
m/s）。
本模拟实验中，原子速度v单位为À
/ps，单位换算涉及常量有
1ev=1.6X10-19J;
C12=1.66X10-27kg;
1À
/ps=100m/s.
3.2模拟设置
3.2.1模拟体系范围
在模拟过程中,虽然缺陷不直接与热浴区域发生相互作用,但模拟体系中间活动区域的大小对缺陷的产生也有影响,这主要是因为系统的能量密度是由PKA能量的大小以及活动区域的体积来控制的[13]。
3.2.2平衡时间
MD模拟的优点在于计算量小、计算时间短。
模拟过程的计算时间长短由平衡时间来决定。
平衡时间是分子动力学模拟中一个重要的参数,如果过大则会增加不必要的计算量；
如果过小,模拟体系未达到平衡就输出模拟结果,使模拟情况不能反映真实情况。
所以选择合适的平衡时间对于计算时间节省和,实验结果的准确性是重要的。
本文对辐照模拟的实验时间均>
10ps。
模拟体系的大小对平衡平衡时间也有影响。
在MD模拟中,原子不断的运动,则原子与原子的之间的相互作用力会不断的改变,这就是能量产生波动的原因。
随着体系的变小,影响原子运动的条件增多,所以原子达到平衡态所需要的时间也就越长。
3.2.3弛豫时间
模拟中，需要初始化每个原子的初始速度和坐标，但初始化的矢量不一定就是对应势函数的最小值，需要最小化来弛豫应力。
体系的初始坐标和初始速度设置以后,在进行模拟体系的性质以前,首先必须使体系进行趋于平衡的过程。
在这个过程中,动能、势能相互转化,当动能、势能、总能量只在某平均值附近波动、趋于平稳时,就可认为体系就达到了平衡。
3.3模拟过程和结果
本文以HCP密排六方晶体结构的金属锆（Zr）为研究对象，使用分子动力学方法对其级联碰撞的过程进行了模拟研究。
首先对锆的结构和性质进行简要介绍，再结合辐照损伤理论和MD模拟方法，对锆级联过程的PKA入射能量的影响进行了模拟实验，并对结果进行了分析。
3.3.1Zr的基本性质
在1789年，德国化学家MartinKlaproth分析了一块锆石，并分离出了锆的“泥土”形式氧化锆，也就是氧化物ZrO2。
随着错合金在原子能工业上的应用,错工业有了迅速的发展，成为一种重要的战略材料,被誉为“原子时代的第一金属”[14]。
Zr在化学元素周期表中位于第5周期,其位置在IVB族的钦（Ti）之下,原子序数为40,平均相对原子质量为91.22。
纯锆在室温下具有HCP密排六方的晶体结构,称为α-Zr晶格常数为a=0.323nm,c=0.515nm。
-Zr在865℃时发生同素异晶转变,转变为体心立方β-Zr,相变时会有较大的体积收缩。
锆是钛的同族元素，也是同素异构体。
862℃以下呈密排六方晶格，称为α锆。
862℃以上呈体心立方晶格，称为β锆。
纯锆的熔点Tm是1850℃左右。
-Zr的结构如图3-2所示。
图3-2α-Zr的晶格结构
由于核性能优异、力学性能适中和加工性能良好,因而错已被普遍用作核动力水冷反应堆的燃料包壳管和结构材料,如压力管、容器管、孔道馆、导向管、定位格架、端塞和其他结构文件,这是错材的主要用途,占整个错加工材的80%[15]在改善锆合金机械性能和抗腐蚀性能，研究人员已经开发出了N18,N36等新锆合金，其焊接性能和力学性能也有所提高。
3.3.2级联碰撞的演化过程
100K温度下金属Zr中能量为1KeV的PKA原子入射模拟的结果表明：
Zr的级联碰撞过程非常短暂，大约是ps量级（1ps=10-12s）。
随着时间的推移，一开始缺陷迅速增多，随后渐渐下降，趋于平缓，最后达到一个稳定状态。
这是因为，当高于Zr原子离位阈值Ed的粒子入射时，晶体内会产生级联碰撞过程，产生间隙原子和空位并存的缺陷。
在碰撞末期，所有原子的能量都降到离位阈值以下，不再撞击晶格原子。
之后，通过间隙原子的弹性复合，出现“急增-衰减-平缓”的现象。
利用lammps软件对辐照损伤过程进行计算机计算，以便更直接地得到级联碰撞过程的演化情况。
观察的结果也和缺陷产生规律一致。
当PKA原子能量增大时，即辐照增强时，级联碰撞将更为剧烈。
这样，将有更多的原子被赋予大于离位阈值的能量，从而有更多的原子从点阵移出。
本文模拟了PKA能量分别是1KeV、5KeV、10KeV入射Zr基体。
结果发现，随着PKA能量的增大，Zr基体内产生的缺陷数目越多，这说明PKA能量的增加会使碰撞加剧，即辐照损伤加剧。
3.3.3MD元胞设置
如下图所示，-Zr为HCP结构，每个晶胞中有两个原子，位于（0，0，0）和（1/3，2/3，1/2）处。
把H
升级会员 