Nb具有优良的高温性能，既可以用来提高结构钢的强度，也可以增加高温合金的热稳定性，在超导体、焊接、电子、光学、机械领域得到了广泛的应用。Nb可以在非互溶性元素构建的多晶体中形成多层结构，这些多层结构的交叉点界面可以吸附辐射缺陷，具有可以减轻辐照损伤的能力[1]，因此Nb可以用作核反应堆结构材料的自愈保护层，是核聚变反应堆的重要候选材料之一[2,3,4]。从20世纪60年代早期的SNAP-50项目到90年代的SP-100项目，Nb都作为重要的核结构材料得到了应用[5]。辐照导致的材料损伤是引起辐照脆化的关键因素之一，在过去的实验研究[6,7]中，主要关注的是Nb及其合金的辐照肿胀和辐照后拉伸性能。

暴露在反应堆辐照环境中的材料会受到中子等粒子的辐照，使晶格原子获得能量，当获得的能量高于晶格原子的离位阈能时，这些原子会发生离位而成为初级离位原子(primary knock-on atom，PKA)，这些初级离位原子会在材料中形成级联碰撞而产生辐照位移损伤，这些位移损伤最终形成以Frenkel对(Frenkel pairs，FP)、空位团簇、间隙原子团簇以及位错环、空间四面体等形式存在的辐照缺陷[8]。发生级联碰撞过程的时间和空间尺度分别为皮秒和纳米量级，在实验上很难同时测量这种时空尺度的缺陷形成过程，然而分子动力学(molecular dynamics，MD)可以很好地模拟辐照损伤(或缺陷)产生和演化的动态过程，可以帮助人们了解辐照导致微观结构发生改变的物理机制，是一种有效的、被广泛应用于辐照损伤的研究方法[9]。关于金属Nb辐照损伤计算已有一些很有意义的成果，比如，Mat-thai等[10]利用Finnis-Sinclar多体经验势计算了单空位和双空位的形成能；Ackland等[11]计算了自间隙原子(self-interstitial atom，SIA)的形成能，发现<110>方向上形成哑铃结构时系统的能量最小；Lee等[12]研究得到修正嵌入原子(modified embedded atom method，MEAM)势函数的截断半径不会影响计算结果的结论，该小组还发现在<110>方向容易形成哑铃型结构。但是对于金属Nb的初级级联碰撞的模拟研究还鲜见报道。

本工作利用MD方法研究不同能量和初始运动方向的PKA在bcc结构金属Nb中所形成的级联碰撞过程，讨论了级联碰撞后的缺陷类型构型和分布情况、空位团簇、间隙原子团簇以及位错环等在微观尺度上的演化过程。

1 计算模型与方法

本工作利用分子动力学程序LAMMPS (large-scale atom/molecular massively parallel simulation)[13]计算了不同能量和初始运动方向的PKA在bcc结构Nb中产生级联碰撞过程以及缺陷形成、演化过程，并使用可视化软件OVITO (open visualization tool)[14]对该过程进行了直观地展现。模型构建时所用的晶格常数a=0.33205 nm，模拟所用模型的x、y、z方向分别对应晶体的[100]、[010]、[001]晶轴，3个方向均为周期性边界条件。根据测试计算，PKA能量越大，缺陷数量达到峰值时形成的缺陷区域越大，级联碰撞退火时间和系统达到稳定结构所需的模拟时间越长。在实施计算时，为了保证模拟的损伤区域保持在模拟体系内不穿越模拟边界，不同PKA能量构建了不同尺寸的模拟体系，对应的原子数和模拟时间具体如表1所示。根据能量守恒定律，一个具有一定动能的中子能产生相应能量的Nb PKA，离位原子动能的一部分以非弹性散射的形式耗散，不过MD不考虑电子自由度，不能模拟非弹性碰撞过程，由于辐照损伤主要是由弹性碰撞引起的，所以，这种近似对计算结果影响不大。在MD模拟计算辐照损伤中，通过两体碰撞能量传递公式和SRIM (the stopping and range of ions in matter)程序计算，模拟设置的PKA的能量5、10、20、30、40和50 keV，分别对应入射中子的能量为310、620、1240、1861、2482和3100 keV。

表1不同PKA能量级联碰撞模拟盒子尺寸和原子数

Table 1Cell size and number of atoms in the atomic displacement cascade simulation with different primary knock-on atom (PKA) energies (EPKA)

新窗口打开|下载CSV

在模拟级联碰撞之前，模拟系统采用粒子数(N)、压强(P)和温度(T)均不变的NPT系综，在0压强环境中对体系进行了60 ps弛豫，使模拟体系温度达到300 K，并形成稳定结构，将此结构作为级联碰撞t=0时刻的结构，在该结构中心附近随机地选取一个Nb原子作为PKA，按照设定的能量和初始运动方向赋予该原子初速度。PKA运动过程中会与其它原子发生相互作用，形成级联碰撞，从而使材料中产生辐照损伤。在级联碰撞的模拟过程中，考虑到nose-hoover热浴不仅更新原子的速度，还会更新原子的位置，由于该热浴不适合模拟级联碰撞，本工作采用了N、体积(V)和能量(E)均不变的NVE系综，同时使用Anderson热浴控制模拟系统温度保持300 K。另外，为了避免PKA运动过程中产生沟道效应，同时尽量使入射方向随机，每个能量的PKA选取了[122]、[135]、[235]、[4 11 9] 4个方向作为初始运动方向，沿每个方向进行5次计算以达到较好的统计性。模拟级联碰撞时，选用的时间步长越小，计算的结果越准确，但模拟所用计算机时会越长，经过测试计算，模拟过程使用改变时间步长方法：体系中运动速度最大的原子每个时间步长运动的距离不超过模拟盒子的1/4，同时每个步长的时间不大于0.01 ps，这样既兼顾计算效率又具有较高的计算准确性。引入PKA后，每100个时间步输出一次原子位置信息，用于后期的缺陷分析。缺陷分析采用Wigner-Seitz方法[9,15]，用t=0时刻晶格原子位置构造Wigner-Seitz原胞，级联碰撞某时刻的体系结构和Wigner-Seitz原胞对比，若Wigner-Seitz原胞内无原子，该位置就是空位，如果Wigner-Seitz原胞内有2个或者2个以上的原子，距离t=0时刻较远的原子为间隙原子。

2 势函数

本工作使用了Lin等[16]修改的Zr-Nb多体嵌入原子势(embedded atom method，EAM)势函数，该势函数已经被成功地应用于Zr-Nb体系的辐照损伤过程，能够很好地描述Zr-Zr、Zr-Nb、Nb-Nb原子间的相互作用[15]，并且取得了很好的计算结果。当原子间距离满足$?\text{<}{?}_{\mathrm{0}}$(${?}_{\mathrm{0}}$为EAM势函数截断半径)时，原子间相互作用使用式(1)所示的Ziegler-Biersack-Littmark (ZBL)短程相互作用的屏蔽势函数(${?}_{??}^{\mathrm{Z}\mathrm{B}\mathrm{L}}$)表示[17]：

式中，${?}_{?}$和${?}_{?}$是相互作用的第$?$个原子和第$?$个原子的原子序数；${?}_{??}$为原子$?$和原子$?$之间的距离；$?$是电子电量；l=$\frac{\mathrm{0.8854}{?}_{\mathrm{0}}}{{?}_{?}^{\mathrm{0.23}}\text{+}{?}_{?}^{\mathrm{0.23}}}$，其中，a0为Bohr半径，a0=${\mathrm{?}}^{\mathrm{2}}$/(${?}_{\mathrm{e}}{?}^{\mathrm{2}}$)=0.0529 nm (其中，${?}_{\mathrm{e}}$为电子质量，约化Planck常数$\mathrm{?}$=$?$/($\mathrm{2}\mathrm{\pi }$)，h=6.62607015×10-34J·s)；$?\text{(}?\text{)}$可表示为：

EAM势和ZBL短程相互作用势之间的结合使用平滑函数(${?}_{\mathrm{s}}\text{(}?\text{)}$)过渡[15]：

其中，${?}_{\mathrm{0}}$=0.23 nm[18]；${?}_{\mathrm{1}}$=0.1 nm[18]，为${?}_{\mathrm{s}}\text{(}?\text{)}$函数的截断半径；系数${?}_{\mathrm{0}}$、${?}_{\mathrm{1}}$、${?}_{\mathrm{2}}$和${?}_{\mathrm{3}}$由${?}_{\mathrm{s}}\text{(}?\text{)}$函数的连续性约束和一阶导数导出[15]。

3 结果与讨论

3.1 级联碰撞过程的Frenkel缺陷对随时间的演化及缺陷分布

一次完整的级联碰撞过程需要经历缺陷的产生、离位峰的出现、缺陷的湮灭以及级联退火而形成稳定缺陷4个阶段[19]。图1所示是能量为30 keV，初始运动方向为[235]的PKA在金属Nb中因级联碰撞产生的空位数(Frenkel缺陷对)随时间的演化曲线，图中还标出了几个主要时刻的级联中心附近的局部原子结构图。本工作模拟结果的曲线特征与其它bcc结构的材料级联碰撞研究结果[20,21,22,23]基本一致，在级联碰撞的初始阶段，空位数迅速增加并达到峰值(弹道阶段)，然后快速下降(恢复阶段)之后达到稳定阶段(10 ps左右以后)。图2是能量为30 keV，初始运动方向为[235]的PKA级联碰撞产生的缺陷在不同时刻的可视化截图，图中蓝色小球表示代表空位，黄色小球表示单间隙原子，红色小球表示双间隙原子。

图1

图1能量30 keV沿方向[235]的初级离位原子(PKA)产生的Frenkel缺陷对随时间的演化

Fig.1The number of Frenkel pairs induced by 30 keV PKA moving in [235] direction as a function of simulated time (Insets show the simulated snapshots of local structures near the cascade center,t—the time of simulation, black lines in insets show the size of the vacancy zone, PKA—primary knock-on atom)

图 2

图 2能量30 keV、沿方向[235]运动的PKA级联碰撞过程产生的缺陷随时间演化过程可视化截图

Fig.2Visual screenshots of defect evolution induced by 30 keV PKA moving in [235] direction during displacement cascades (The blue, yellow and red atoms represent vacancy, single interstitial atom and double interstitial atom, respectively)

(a)t=0.051 ps (b)t=0.11 ps (c)t=0.38 ps (d)t=1.49 ps

(e)t=5.09 ps (f)t=10.09 ps (g)t=15.09 ps (h)t=25.00 ps

结合图1可知，t=0.051 ps和t=0.11 ps属于弹道阶段，随着时间的推移，金属Nb中的空位数和间隙原子数都在增加，空位集中分布在级联中心处，间隙原子分布在空位周围，空位和间隙原子分布呈现为不规则的形状(图1中插图、图2a、图2b)。其它能量和方向PKA产生级联碰撞的弹道阶段时缺陷分布都有类似的结果，不同之处只是缺陷区域的尺寸。从该阶段缺陷演化过程中可以观察到，<111>、<110>方向出现了空位和间隙原子交替，并形成2个间隙原子与1个空位原子构成的哑铃状结构的换位挤塞子序列(图2a和b圈中的部分)，间隙原子通过这种换位挤塞子序列迁移到远离级联中心的位置。当t=0.38 ps时，Frenkel对数目接近峰值。结合图1和2可以看出，这时缺陷区域进一步扩大，级联中心的大量原子离位，在级联中心附近出现了密集的空位区，形成空位团簇，出现了比较明显的“空洞”(图1)，外围的间隙原子比较均匀地分布在“空洞”周围，这时空位和间隙原子分布大致呈球形(图2c)。当t=1.49 ps时，Frenkel对数量达到峰值，空位形成的“空洞”进一步增大(图1、图2d)，空位形状仍近似球形，间隙原子分布却呈现出各向异性，沿<111>方向的间隙原子多于其它方向。其原因是随着Frenkel对数量的增加，间隙原子会在晶体内部迁移，在bcc结构中，<111>方向的原子离位阈能较低[24,25,26,27]，致使该方向间隙原子更容易通过级联碰撞中的换位挤塞子序列迁移到离级联碰撞中心更远的位置。在弹道阶段出现的<111>挤塞子序列也是这个原因。由于碰撞级联而产生“热峰”效应，大量的间隙原子具有较高的能量，会与附近的空位发生复合，t=5.09 ps呈现的是复合阶段原子构型。由于空位与临近的间隙原子发生复合，使Frenkel对数目比峰值时刻明显减少。在远离级联碰撞中心的周围空位和间隙原子最先发生复合，而级联中心区域的空位和间隙原子相距较远，不容易发生复合(图1、图2e)。从t=10.09 ps开始，Frenkel对数目基本保持不变，复合阶段结束。从此时开始，体系中除了一些零星的点状缺陷(图2f圈中)外，还出现了空位成团和间隙原子成团的现象，这些空位团簇和间隙原子团簇之间相距较远，无法进一步复合，形成最终的稳定缺陷区。从t=15.09 ps到25.00 ps (模拟计算终止时刻)，缺陷数目和结构基本不发生变化，只有零星的缺陷移动，缺陷团簇基本保持不变(图2g和h)。

分析所有模拟结果发现，当EPKA较高(40和50 keV)时，碰撞级联分裂成多个行为相对独立的子级联。对级联碰撞后的稳定缺陷分析可知，弥散分布在研究空间内的间隙原子形成了很多哑铃状结构，对这些哑铃状结构位置坐标分析发现，它们大都沿<110>取向，其原因是沿该方向点缺陷的迁移能较低。在长时间的演化中，哑铃结构会在材料的晶界、位错等位置处发生湮灭[28,29]。另外，当EPKA为30、40和50 keV时，模拟计算15 ps左右后均出现大小不等的位错环，随着模拟演化时间的延长，这些位错环不断运动、融合，在计算结束时，最终形成了1/2<111>间隙型和<100>空位型2种不同类型的位错环。图2h所示为方向[235]、能量30 keV PKA的模拟中最终形成的3个1/2<111>间隙型位错环和1个<100>空位型位错环，在其它bcc结构材料中也发现了类似的结果[22]。在缺陷复合过程中，由于空位的迁移能比间隙原子高[26]，空位型缺陷迁移更困难，扩散能力更差，所以最终形成的空位型的点缺陷、空位团簇以及空位型的位错环更接近于级联中心区域[25]。

3.2 存活缺陷的数量

图3给出了不同EPKA引发级联碰撞产生的Frenkel对数量(NFP，每个能量PKA沿所有方向的平均数目)随时间演化的曲线。从图中可以看出，随着EPKA的增加，出现缺陷数量的峰值以及最终达到稳定状态所需的时间会变长。最终形成稳定的缺陷数目与EPKA相关，EPKA越大，缺陷数目越多，这些缺陷的数量可以用Norgett、Robinson和Torrens提出的模型(NRT模型)[27]来估算。但是Bacon等[30]经过多年的模拟研究发现，bcc结构的材料稳定阶段NFP和EPKA之间更符合幂函数关系：

图3

图3不同能量PKA级联碰撞过程中产生的缺陷数目随时间的演化

Fig.3The number of Frenkel pairs as a function of simulated time for different PKA energies (EPKA)

式中，b和m为幂函数参数，是受辐照材料的种类和温度影响的较小的常量。利用式(4)对本模拟稳定阶段的NFP和EPKA的关系进行了分段拟合，得到很好的拟合优度，EPKA为5~30 keV，拟合得到的b=2.699，m=0.895；EPKA为30~50 keV范围，拟合得到的b=0.786，m=1.261，拟合结果如图4所示。此结果与Allen等[31]的研究结果相符。对稳定阶段的缺陷构型分析发现，EPKA较低时，在级联损伤区形成的只有点缺陷，而EPKA较大时，会在模拟体系中引发二级级联碰撞甚至多级级联碰撞，从而形成缺陷的数目增大较多，导致幂函数拟合时的幂指数增大。如前所述，EPKA大于30 keV形成的缺陷中均发现了空位型位错环和间隙型位错环，这些位错环阻碍了空位和间隙原子的进一步迁移，从而使更多的缺陷存活，这是导致EPKA在不同区间对应的拟合系数不同的另一个原因。

图4

图4存活Frenkel对数目和级联碰撞效率随PKA能量的变化

Fig.4Curves of average number of surviving Frenkel pairs and cascade efficiencyvsEPKA

为了比较NRT模型和幂函数模型的不同，定义${?}_{\mathrm{F}\mathrm{P}}$/${?}_{\mathrm{N}\mathrm{R}\mathrm{T}}$为级联碰撞的效率(cascade efficiency)，其中${?}_{\mathrm{N}\mathrm{R}\mathrm{T}}$=$\mathrm{0.8}{?}_{\mathrm{P}\mathrm{K}\mathrm{A}}$/($\mathrm{2}{?}_{\mathrm{d}}\text{)}$(其中，${?}_{\mathrm{d}}$是有效平均离位阈能，在本工作中，使用ASTM标准E521-16推荐的60 eV)，为NRT模型计算的缺陷数量。计算得到不同能量的级联碰撞效率结果如图4所示。Bacon等[30]和Voskoboinikov等[32]和Gao等[33]计算非难熔bcc结构材料的级联碰撞效率的值在0.10~0.15之间，本工作的计算值在0.301~0.336之间，并且随着EPKA的增加，效率逐渐降低，这与Selby等[34]对难熔金属Mo以及Bj?rkas等[24]对难熔金属W和V等的计算结果类似，虽然工作的计算结果偏大，但效率随EPKA变化的趋势相同。

3.3 缺陷团簇分析

辐照损伤引起的缺陷团簇影响点缺陷的迁移、回复和湮灭，还会影响材料辐照后的热力学性能。本工作还对形成的空位和间隙原子缺陷团簇进行了分析。分析时空位团簇的截断半径为第二近邻 (second nearest neighbor，2NN)距离，间隙团簇截断半径为第三近邻(third nearest neighbor，3NN)距离，使用式(5)计算缺陷成团率(${?}_{\mathrm{c}\mathrm{l}}$)：

式中，${?}_{\mathrm{c}\mathrm{l}}$表示因级联碰撞形成团簇的点缺陷原子个数。在计算成团率时，用同一能量4个方向PKA级联碰撞的稳定阶段的最后5 ps的平均值，计算结果如图5所示。从分析结果可以看出，空位团簇分数在17%~35%之间，间隙原子团簇分数在23%~40%之间，间隙原子团簇率略高于空位团簇率。可见，金属Nb中级联碰撞形成的缺陷大部分还是以点缺陷的形成存在。根据文献[35]报道，Nb原子在<111>方向上的迁移能为0.008 eV，空位迁移能0.30 eV以上，间隙原子迁移能远低于空位迁移能，这就会使间隙原子更容易聚集形成团簇，如图2的缺陷结构所示。本工作计算结果还显示，随着EPKA的增加，空位团簇成团率有所增加，当EPKA大于20 keV，空位团簇形成率基本不再变化，这与Bj?rkas等[24]对V和Selby等[34]对Mo等难熔金属的模拟结果基本一致。Bj?rkas等[24]认为，对于熔点较高的难熔金属，由于“热峰”效应而形成的再结晶前缘运动较快，导致空位团簇分数较低。与空位相比，EPKA越高，间隙原子团簇率越高。

图 5

图 5间隙原子和空位的成团率随EPKA的变化

Fig.5Curves of vacancy cluster and interstitial cluster formation ratesvsEPKA

本工作对不同能量PKA形成的缺陷团簇尺寸(定义为团簇含有点缺陷个数)分布做了统计(对同一能量，不同初始运动方向PKA形成缺陷的统计平均)，结果如图6所示。从图中的缺陷分布可以看到，较高能量PKA更容易形成大的缺陷团簇。对不同类型的团簇尺寸比较可以看出，大部分间隙原子团簇为小型团簇，EPKA高达50 keV形成的最大团簇为7，而空位更容易形成较大的团簇，当EPKA为50 keV时，空位型缺陷形成了较大团簇，甚至出现了22个空位的大型团簇。

图6

图6不同能量PKA间隙原子团簇和空位团簇的尺寸和数量分布

Fig.6The number of interstitial clusters (a) and vacancy cluster (b) formed in each cascade as a function of the correspondingEPKA

4 结论

(1) PKA引发级联碰撞会在体系中产生辐照缺陷。这些缺陷大多数以点缺陷的形式存在，空位团簇率为17%~35%，间隙原子团簇率为23%~40%。

(2) 级联碰撞产生的间隙原子形成了大量的沿<111>和<110>方向的哑铃型结构。

(3) 当EPKA较大(大于30 keV)时，级联碰撞导致材料中产生了1/2<111>间隙型位错环和<100>空位型位错环。

(4) 对不同的EPKA区间(5~30 keV和30~50 keV)，Frenkel缺陷对数目满足不同的幂函数关系，拟合系数分别为b=2.699和m=0.895以及b=0.786和m=1.261。