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

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

《锆级联的分子动力学模拟.docx》由会员分享，可在线阅读，更多相关《锆级联的分子动力学模拟.docx（25页珍藏版）》请在冰豆网上搜索。

锆级联的分子动力学模拟

锆辐照级联的分子动力学模拟

研究小组成员：

韩维强（主持人）、司宇、张锦文

指导老师：

侯怀宇

2012.3~2013.9

1绪论

1.1研究背景和意义

能源,信息,材料,通常被称为新科技革命的三大支柱。

人类文明的进程和能源的发展革命息息相关。

自工业革命以来，以化石燃料作为主要能源的时代已经持续了两百多年。

当今世界，人类对能源的需求与日俱增，而地球上探明的可使用的化石燃料如煤炭，石油，天然气和其他可燃气体迅速减少，因而寻找一种可替代化石能源的其他能源是迫切而必然的选择。

半个多世纪以来，核电，火电和风力发电一直作为当今世界电力的三大支柱发挥着不可代替的作用。

核能，以其高效，清洁，安全的优点备受青睐。

在我国，目前已建成秦山核电站、广东大亚核电站、江苏田湾核电站和岭澳核电站四座核电站，为我国电力事业发挥巨大作用[1]。

核材料应具备以下条件：

a.符合核条件。

材料的吸收中子截面要小；杂质含量满足核级纯；低活性，感生放射性小，剩余发热量小，半衰期短。

b.运行性能。

应具备良好的物理性能；高抗辐射性能；与冷却剂相容性好。

c.经济适用性。

材料具备工业生产能力，可批量生产；低成本；焊接性好。

基于各项性能考虑，工业应用的核材料主要有铝及铝合金、镁合金、不锈钢、锆合金。

锆合金比不锈钢的熔点高300-400℃,热胀系数小2/3,热导率高18%,热中子吸收截面小一个量级;力学性能、加工性能好,同UO2相容性好,尤其是对300-400℃的高温水、高温蒸汽也具有良好的抗腐蚀性能和足够的热强性。

基于此,锆被选作反应堆燃料。

元件的包壳材料和结构材料,如压力管、容器管和定位格架等,是“核安全最坚韧防线”。

核反应堆的必要条件除了等离子物理外，最为重要的就是壁材料。

所以说，核反应堆的诞生基于锆的发展、成熟[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.5