钨中自间隙原子团簇扩散行为的分子动力学模拟

李小椿，潘新东,2，周海山,*，罗广南,2

(1.中国科学院 等离子体物理研究所，安徽 合肥 230031；2.中国科学技术大学，安徽 合肥 230026)

核聚变能因其对环境友好和不产生放射性废料的高安全性等优势成为一种清洁的安全能源，其燃料——氘在海水中几乎取之不尽、用之不竭，是解决能源危机最有效的方法之一。然而，可控热核聚变技术现在还处于一个相当不成熟的阶段。因为库仑排斥作用使核聚变反应非常困难，须使反应气体电离成为等离子体并加热至极端的高温。当前一般使用磁约束热核聚变的方法，托卡马克装置是迄今为止最有发展前途的磁约束热核聚变装置[1-2]。国际上已启动国际热核聚变实验反应堆(ITER)计划[3-4]，为将来进一步发展聚变反应堆奠定基础。中国也正在进行聚变试验堆(CFETR)的工程设计[5-7]。

聚变能源应用的最终实现除要解决可控热核聚变这一物理问题外，在很大程度上取决于可控热核聚变装置托卡马克以及未来反应堆中关键材料问题的解决。核能界公认聚变堆材料是开发核聚变能的最关键技术之一，其中面向等离子体材料(PFM)的选择尤为关键[8]。在聚变堆严酷的工况环境下，PFM受三重辐照：等离子体辐照、热辐照、中子辐照。高能聚变中子辐照导致PFM中产生大量缺陷，这些缺陷的迁移和聚集会引起材料的肿胀和变形，从而严重影响材料的力学性能，降低部件的使用寿命，最终必将导致材料在服役条件下的失效。如此极端的服役条件对PFM的开发提出了巨大的挑战[9]。

钨(W)材料被视为未来聚变堆中最可能全面使用的PFM，其主要优势有高熔点、高热导率、不与氢发生化学刻蚀以及相对低的氢滞留等[10-12]。中国科学院等离子体研究所的EAST托卡马克将逐步过渡到全W的PFM[13]，同时W也是ITER和CFETR的重要候选材料[5,14]。

在未来聚变堆真实环境下，氘氚聚变产生的14 MeV高能中子辐照将在材料中产生严重的原子离位损伤和各种缺陷积累[15]。其中自间隙原子(self-interstitial atom, SIA)及其团簇是中子辐照损伤中最常见的缺陷种类，因此研究SIA团簇聚集及其动力学扩散行为具有重要意义。由于纯中子源辐照均存在样品活化带来的后续问题，因此国际上一直试图利用离子束辐照产生缺陷，以避免中子活化带来的巨大麻烦[16-17]。材料模拟是研究极端条件下材料行为的有效工具，其中分子动力学(MD)模拟是研究材料辐照损伤效应的重要工具[18]，是材料研究领域中对微观结构和微观过程研究的重要手段。本文采用MD模拟研究不同尺寸SIA团簇的扩散行为，拟为更大尺度的动力学蒙特卡罗和团簇动力学模拟提供准确的输入参数，为正确掌握和评价W中子辐照行为提供依据。

1 研究方法

在分子静力学和MD模拟方面采用开源并行程序LAMMPS[19]，LAMMPS是经典的免费开源MD模拟软件。本文所涉及的模拟均通过LAMMPS软件实现。MD模拟的准确性极大依赖于模拟选用的势函数，为正确掌握SIA团簇的性质，针对不同的性质，本文选用两套势函数开展相关模拟工作。

Ef(SIAn)=Etot(SIAn)-NWEc(W)

(1)

其中：Etot(SIAn)为包含n个SIA的体系总能量；NW为总的W原子数；Ec(W)为单个W原子在完美体心立方晶格中的内聚能。单个SIA与SIA团簇的结合能定义为：

Eb(SIA+SIAn-1)=

Ef(SIA)+Ef(SIAn-1)-Ef(SIAn)

(2)

v=v0exp(-Em/kBT)

(3)

其中：v0为跃迁频率指前因子；Em为跃迁激活能；kB为玻尔兹曼常数；T为温度。对式(3)两边分别取对数，可得到：

lnv=lnv0-Em/kBT

(4)

可见跃迁频率的对数与温度的倒数满足线性关系，线性拟合后就可得到跃迁频率指前因子和跃迁激活能。

2 结果与讨论

2.1 SIA团簇的形成能与结合能

首先采用分子静力学模拟获得了数量为1～30个的1/2〈111〉和〈100〉 SIA团簇的稳定结构、形成能和结合能，这部分模拟采用Chen势完成。图1示出了1/2〈111〉和〈100〉 SIA团簇的部分稳定构型。表1列出了这些团簇的形成能和结合能，作为参考也列出了第一原理计算的结果[24]。由表1可看出，1/2〈111〉 SIA团簇的形成能要低于〈100〉 SIA团簇，说明1/2〈111〉 SIA团簇较〈100〉 SIA团簇稳定。SIA团簇的形成能随团簇尺寸的增长而增加，且SIA与SIA团簇的结合能非常高，说明SIA团簇聚集后会稳定存在。

图1 1/2〈111〉和〈100〉 SIA团簇的最稳定结构Fig.1 The most stable structure of 1/2〈111〉 and 〈100〉 SIA clusters

表1 1/2〈111〉和〈100〉 SIA团簇的形成能与结合能Table 1 Formation energy and binding energy of 1/2〈111〉 and 〈100〉 SIA clusters

为给后续大尺度模拟提供更准确的输入，研究了1/2〈111〉和〈100〉 SIA团簇尺寸大于30时的形成能，结果如图2所示。受计算量限制，无法给出每个尺寸的SIA团簇的形成能，可采用以下经验公式拟合[24]：

(5)

采用团簇尺寸大于30之后的形成能拟合，得到了1/2〈111〉 SIA团簇的拟合参数为：a0=6.801 17、a1=-11.021 26、a2=60.138，而〈100〉 SIA团簇的拟合参数为：a0=8.845 46、a1=-20.636 8、a2=93.665 1。这些参数可为后续的大尺度模拟提供每个尺寸的SIA团簇的形成能。相应地，各团簇之间的结合能也可根据形成能轻易获得。据此，可为后续的大尺度模拟提供W中最常见的两种SIA团簇的形成能和结合能的完备数据库。

图2 大尺寸SIA团簇的形成能与经验拟合Fig.2 Formation energy of large SIA clusters and empirical fitting

2.2 SIA团簇的扩散行为

采用Ackland势模拟研究不同尺寸的1/2〈111〉 SIA团簇在不同温度下的扩散行为。SIA团簇的尺寸变化范围为1～109个，温度范围为300～900 K。对于每个SIA团簇在每个温度下的模拟，均通过随机产生初始速度的方法运行20次，因此有共20 ns的模拟数据，随后采用OVITO软件分析SIA团簇的运动轨迹。

通过对SIA运动轨迹的分析发现，在300～600 K温度条件下，单个SIA在MD模拟时间尺度内未观察到转向，表现为一维方向的快速运动，如图3a所示，SIA只在x方向(〈111〉方向)移动，而在另外两个方向几乎不动。只有当温度大于700 K时，单个SIA才可转向，如图3b所示，SIA在0.37 ns发生转向后，易在0.7 ns时再次转回稳定的〈111〉方向。随着温度的升高，单个SIA越易转向。在700 K的模拟下，20个模拟里只有2个模拟观察到了SIA的转向。而在800 K的模拟下，有13个模拟观察到SIA的转向。

图3 单个SIA在600 K(a)和800 K(b)下的扩散坐标与时间的关系Fig.3 Time dependence of diffusion coordinates of single SIA at 600 K (a) and 800 K (b)

而SIA团簇尺寸超过2个后，在300～900 K时，这些SIA团簇均非常稳定，表现为一维方向的运动。需说明的是，虽然模拟的20 ns的时间内未观察到2个以上SIA团簇的转向，并不代表这些尺寸的SIA团簇不会发生转向，受限于MD的局限，可能需更长的时间才能观察到。将模拟温度升高到1 200 K可发现，2个SIA团簇在1 000～1 200 K、3个SIA团簇在1 200 K时可观察到SIA团簇的转向，而4个以上的SIA团簇在1 200 K下也未发现转向。

2.3 SIA团簇的跃迁频率

为掌握SIA团簇的扩散行为，进一步分析了SIA团簇在不同温度下的跃迁频率，结果如图4所示。SIA团簇尺寸越大，跃迁频率越低，扩散越缓慢，而温度越高，SIA团簇扩散越快。根据式(4)对SIA团簇在不同温度下的跃迁频率的对数与温度的倒数进行线性拟合，可得到跃迁频率指前因子和跃迁激活能，拟合曲线如图5所示，在300～900 K下，SIA团簇的跃迁频率vn的对数与温度的倒数有着良好的线性关系，说明可采用式(3)定义跃迁频率与温度的关系。根据拟合结果，在表2列出了各种尺寸SIA团簇的跃迁能垒Em和跃迁频率指前因子v0。跃迁频率指前因子随SIA团簇尺寸的增加而递减，而跃迁能垒却几乎与SIA团簇尺寸无关。通过分析SIA团簇的运动轨迹可发现，SIA团簇里的每个SIA均在一条独立的〈111〉方向里移动，因此整个团簇的跃迁能垒与单个SIA的跃迁能垒类似。

图4 SIA团簇在不同温度下的跃迁频率Fig.4 Vibration frequencies of SIA clusters at different temperatures

图5 SIA团簇的跃迁频率与温度倒数的关系Fig.5 Relationship between vibration frequency and reciprocal temperature of SIA clusters

表2 SIA团簇的跃迁频率指前因子和跃迁能垒Table 2 Pre-exponential factor of vibration frequency and diffusion barrier of SIA clusters

在真实的环境中，SIA团簇有各种不同的尺寸，因此需根据经验公式，外推出各种尺寸SIA团簇的跃迁频率指前因子和跃迁能垒。SIA团簇的跃迁频率与团簇尺寸满足以下关系：

vn=v0n-Sexp(-〈Em〉/kBT)

(6)

其中：S为可拟合的参数；〈Em〉为平均跃迁激活能。对表2中间隙团簇尺寸大于7的数据进行拟合，可得到1/2〈111〉 SIA团簇的相关参数：v0=2.226×1012Hz、S=0.498和〈Em〉=0.021 7 eV。

图6为SIA团簇的跃迁频率指前因子与团簇尺寸的关系。从图6可看出，拟合参数可很好描述1～109个SIA团簇的跃迁频率指前因子的关系。据此，就可为后续的大尺度模拟提供每个SIA团簇的跃迁频率和跃迁能垒等参数。

图6 SIA团簇的跃迁频率指前因子与团簇尺寸的关系Fig.6 Relationship between pre-exponential factor of vibration frequency and SIA cluster size

3 结论

采用MD模拟研究了W中1/2〈111〉和〈100〉 SIA团簇的稳定结构，发现W中SIA团簇最稳定结构是1/2〈111〉 SIA团簇结构，且SIA与SIA团簇的结合能非常高，说明SIA团簇聚集后会稳定存在。给出了1/2〈111〉和〈100〉 SIA团簇的形成能和结合能，并给出了计算大尺寸SIA团簇形成能的经验参数。随后研究了1/2〈111〉 SIA团簇的动力学扩散行为，单个SIA在温度高于700 K时易扩散和转向。而两个以上的SIA团簇在300～900 K时主要表现为一维方向的运动，未观察到转向。根据不同温度下的跃迁频率，得到各种尺寸SIA团簇的跃迁能垒和跃迁频率指前因子。为准确描述各种尺寸SIA团簇的动力学行为，给出了一套计算SIA团簇跃迁频率的经验参数。目前的模拟结果将有助于理解W中SIA的团簇行为，从微观角度理解和解释实验，并为大尺度的模拟提供准确和完备的输入参数。