δ相钚-镓合金自辐照损伤的分子动力学模拟

卢焘， 田新， 王继虎

(北京航空气象研究所， 北京 100085)

钚是一种重要的核材料，为了确保核设施的安全可靠，需要深入了解钚的结构与性质。钚能够自发α衰变，形成氦核(约5.14 MeV)和铀核(约86 keV)，在晶格中引起级联碰撞，从而产生过饱和的点缺陷(弗仑克尔缺陷对)。过饱和的点缺陷经过迁移、聚集能够演化为团簇、位错、位错环和孔洞等更为复杂的本征缺陷，这些本征缺陷能够捕获氦核形成氦泡。钚自发α衰变导致的材料内部损伤是钚固有的自辐照效应，将显著地影响材料的结构与性质。例如位错和位错环会影响材料的延展性，氦核会引起材料氦脆，氦泡和孔洞会造成材料肿胀。自辐照效应是钚材料老化的主要原因，也是评估核设施可靠性与安全性的重要依据[1-2]。

钚的晶体结构极其复杂，在从室温到熔点(约 640 ℃)的较小温度范围内存在6种同素异形体，各种同素异形体之间的性质差异极大，其膨胀、电阻和导热性不随温度线性变化，因此冶金工艺十分复杂。常温下，钚的稳定相为简单单斜结构的α钚，然而α钚性能脆，耐腐蚀性差，不易加工和保存。面心立方结构的δ钚具有良好的延展性，其稳定存在的温度范围为319～450 ℃。为了将δ钚稳定到室温，人们在δ钚中加入少量镓、铝等合金化元素，得到常温下亚稳态的δ相固溶体合金[3]。

由于钚性质异常独特，其放射性很强、毒性极大，同时化学性质活泼，具有自燃性，在空气中极易氧化腐蚀，因此钚的实验研究十分困难。分子动力学模拟是研究极端环境下材料结构与性质的有效工具。对金属和合金而言，表征原子间相互作用通常采用嵌入原子(embed atom method，EAM)势[4]，然而EAM势函数的一个假设是原子的电荷密度呈球对称分布，因此并不适用于钚。Baskes等保持EAM势原有的理论框架不变，针对原子的电荷密度球形对称的假设加以修正提出了修正嵌入原子(modified embed atom method，MEAM)势[5]。MEAM势在基体电荷密度求和中引入原子电荷密度分布的角度依赖因素，对s、p、d态电子的分布密度分别进行计算，但电子总密度仍然等于各种电子密度的线性叠加。Valone等基于Baskes等提出的MEAM势模型，拟合了钚-镓-氦势函数[6]，已得到了国内外研究者的广泛应用[7-10]。

MEAM势能够很好地描述Pu-Ga-He体系的热力学性质，但涉及自辐照损伤模拟的问题时，由于原子在级联碰撞的瞬间距离极近，MEAM势不再适用这种极端情况。针对级联碰撞问题，最常用的势函数为齐格勒-比尔扎克-利特马克(Ziegler-Biersack-Littmark，ZBL)势[11]。

本文采用分子动力学方法，结合MEAM势和ZBL势，模拟了δ相钚-镓合金中原子的级联碰撞，并探索了温度、镓、氦、位错等的影响，这些研究有助于从原子尺度理解钚的自辐照效应机理。

1 模型与方法

本文所有分子动力学模拟都通过大尺度原子/分子模拟软件(large-scale atomic/molecular massively parallel simulator，LAMMPS)进行[12-13]。模拟结果后处理采用开源可视化工具(open visualization tool ，OVITO)[14-15]。

势函数选择MEAM势，其函数形式为：

对势函数φ形式为：

Eu函数形式为：

式中：Ec、re、Ω、B分别为参考结构的原子内聚能、最近邻原子间距、原子体积、体积模量；δ为修正参数。

式中ρa(0)为原子电荷密度。

嵌入势F函数形式为：

式中：A为可调参数；ρ0为密度尺度参数。

式中：ρ(0)为球对称电荷密度；Γ为角修正参数，其函数形式为：

式中：t(l)为拟合常数;ρ(l)的函数形式为：

ρ(l)(R)=ρ(0)exp[-β(l)(R/re-1)]

式中β(l)为拟合常数。

对于原子间距小于0.1 nm的短程相互作用，采用ZBL势替代MEAM势的对势项。ZBL势的函数形式为：

式中：ε0为真空介电常数；Zi和Zj为原子序数；e为电子电荷。

MEAM势采用Valone等拟合的参数[6]，如表1～3。

表1 单一元素势函数参数Table 1 Parameters for pure elements

表2 不同元素相互作用势函数参数

表3 势函数角筛选参数Table 3 Parameters for angular screening

随机选取不同比例的钚原子置换为镓原子，获得不同镓浓度的δ相钚-镓合金模型。

为考察氦原子对δ相钚-镓合金自辐照损伤的影响，本文创建了2种模型，第1种是直接在上述模拟盒子中随机添加数密度为0.1%的氦原子，第2种是在上述模拟盒子中随机选取0.1%的钚原子或镓原子置换为氦原子。

在得到上述各种模型后，需要采用共轭梯度算法对体系结构进行优化，而后在等温等压系综(constant pressure-constant temperature ensemble， NPT)下，通过Nose-Hoover恒温器算法控制体系温度，通过Parrinello-Rahman恒压器算法[16]控制体系压强为零，采用速度Verlet算法[17]让模型在不同温度下运行20 ps(时间步长1 fs，运行20 000步)，使模型达到热力学平衡。

将模拟盒子分为2个部分，厚度为1 nm的边界层在正则系综(canonical ensemble， NVT)下通过Nose-Hoover恒温器算法控制温度，内部原子在微正则系综(micro-canonical ensemble， NVE)下控制体积和总能量不变，引入不同能量的初级离位原子(primary knock-on atom, PKA)，采用0.01～1 fs的变化时间步长，控制任一时间步长内任一原子的位移小于0.01 Å，通过速度Verlet算法让模型在不同温度下运行1 ns，每100步记录一次所有原子的坐标。

利用Wigner-Seitz缺陷分析方法寻找级联碰撞过程中的空位和间隙原子。具体方法为：根据初始构型划分Wigner-Seitz原胞作为参考结构，在后续模拟中，若某一初始原胞空间内出现2个或以上原子，则该原胞空间内的原子为间隙原子，若某一原胞空间内无原子，则该原胞空间为空位。

利用晶体近邻原子分析法得到晶体结构参数，通过密堆六方结构寻找位错。

2 温度、镓浓度对自辐照损伤的影响

为考察不同温度、不同镓浓度对δ相钚-镓合金自辐照损伤的影响，本文模拟了不同温度下(300 K和600 K)，不同镓浓度(3%和6%)的δ相钚-镓合金受到不同能量(1 keV和5 keV)初级离位原子辐照后的级联碰撞过程。如图1所示为300 K温度下，δ相钚-镓(3%)合金受到5 keV初级离位原子辐照后，不同阶段的点缺陷分布图。由图可知，级联碰撞分为2个阶段：第1阶段，初级离位原子将近邻原子撞离晶格位，被撞离晶格位的离位原子继续将近邻原子撞离晶格位，这一阶段级联碰撞占据主导地位，点缺陷迅速增加；随着离位原子能量的降低，级联碰撞效应逐渐减弱，大约3 ps，点缺陷达到峰值，此后进入第2阶段，即离位峰弛豫阶段，缺陷复合效应占据主导，点缺陷缓慢减少，直到30 ps后达到亚稳态。

图2(a)为1 keV初级离位原子辐照后，级联碰撞过程中产生的空位数随时间的变化，温度分别为300 K和600 K，镓浓度分别为3%和6%。由图可知，温度主要影响了峰值空位数，600 K温度下，峰值空位数较大，但30 ps后，不同温度条件下遗留的空位数几乎相同。另一方面，镓浓度对峰值空位数影响不大，但镓浓度越高，30 ps后遗留的空位数越大。

图1 300 K温度下，5 keV初级离位原子，级联碰撞过程中空位和间隙原子分布Fig.1 Distribution of vacancies and interstitial atoms during the collision cascades with 5 keV PKA at 300 K

图2 1 keV、5 keV初级离位原子级联碰撞过程中产生的空位数随时间的变化Fig.2 The vacancy number depending on the time during the collision cascades with 1 keV PKA and 5 keV PKA

图2(b)为5 keV初级离位原子辐照后，级联碰撞过程中产生的空位数随时间的变化，温度分别为300 K和600 K，镓浓度分别为3%和6%。由图可知，温度和镓浓度都显著地影响了峰值空位数和遗留空位数，温度越高，镓浓度越大，辐照后产生的峰值空位数和遗留空位数都越大。

由上述讨论可知，δ相钚-镓合金的自辐照损伤在高温环境下更为显著，而镓元素的存在会加剧自辐照损伤。

3 氦对自辐照损伤的影响

为考察氦对δ相钚-镓合金自辐照损伤的影响，本文模拟了300 K温度下，含有0.1%氦间隙原子或0.1%氦置换原子的δ相钚-镓合金受到不同能量(1 keV和5 keV)初级离位原子辐照后的级联碰撞过程。

由于氦原子溶解于δ相钚-镓合金的最稳定形态是占据晶格位，将钚原子挤到间隙位[18]。因此第1种模型等效于在晶格中加入0.1%的空位-氦对和0.1%的钚间隙原子，第2种模型等效于在晶格中加入0.1%的空位-氦对。级联碰撞过程中，增加的空位数随时间的变化如图3所示。

图3 初级离位原子能量1 keV、5 keV，在300 K温度下级联碰撞过程中增加的空位数随时间的变化Fig.3 The additional vacancy number depending on the time during the collision cascades with 1 keV PKA and 5 keV PKA at 300 K

比较无缺陷模型和2种模型的结果可知，空位-氦对的存在使峰值空位数增大，最终遗留空位数也增大。2种模型的区别在于是否预先存在钚间隙原子，因此比较2种模型的结果可知，钚间隙原子的存在使峰值空位数减小，最终遗留空位数增大。上述现象出现的原因可能是，氦原子与空位具有很强的结合能，因而在级联碰撞的过程中，产生的空位将与氦原子结合，从而抑制了空位与间隙原子的复合。需要注意的是，钚的自辐照本身就会产生大量氦，而氦又会加剧钚的自辐照损伤，这一机制是钚材料老化的主要原因。

4 位错对自辐照损伤的影响

图4展示的是4种不同入射方向情况下，空位数随时间的变化。由图可知，1/2〈110〉{111}刃位错的存在显著地降低了峰值空位数和最终遗留的空位数，即位错的存在显著地抑制了材料的自辐照损伤。

为了分析位错对材料自辐照损伤的抑制机制，考察了级联碰撞过程中不同阶段位错和点缺陷的分布形态。图5展示的是300 K温度下，δ相钚-镓(5%)合金内1/2〈110〉{111}刃位错附近，5 keV的初级离位原子沿-y方向入射后，不同阶段的位错、空位和间隙原子分布图。在级联碰撞过程中，位错能够吸收空位和间隙原子，使得峰值空位数和最终遗留的空位数大幅降低。在经过1 ns时间弛豫后，位错通过缓慢吸收空位和间隙原子，最终发生攀移，形成一个割阶。这种现象产生的机制可能是位错的存在引起晶格形变，在位错线两侧分别产生一个拉伸应力区和一个压缩应力区。拉伸引力区能够吸收间隙原子，压缩应力区能够吸收空位，从而使得位错向点缺陷方向发生攀移。

图4 在300 K温度下，δ相钚-镓(5%)合金内1/2〈110〉{111}刃位错附近，5 keV初级离位原子能量沿不同方向入射后，级联碰撞过程中产生的空位数随时间的变化Fig.4 The vacancy number depending on the time during the collision cascades with 5 keV PKA from different directions around the 1/2〈110〉{111} edge dislocation at 300 K

图5 在300 K温度下，δ相钚-镓(5%)合金内1/2〈110〉{111}刃位错附近，5 keV初级离位原子沿-y方向入射后，位错、空位、间隙原子分布图Fig.5 Distribution of the dislocation, vacancies and interstitial atoms during the collision cascades with 5 keV PKA along -y direction around the 1/2〈110〉{111} edge dislocation at 300 K

5 结论

1)钚的自辐照损伤在高温下更为显著，其中，1 keV辐照下温度主要影响级联碰撞过程中的峰值空位数，5 keV辐照下温度对级联碰撞过程中的峰值空位数和最终遗留的空位数都有显著影响；

2)镓和氦都会加剧钚的自辐照损伤，其中，镓和氦都以替代原子存在并显著增加级联碰撞过程中的峰值空位数和最终遗留的空位数，而被氦替代原子挤出的钚间隙原子减少了级联碰撞过程中的峰值空位数并增加了最终遗留的空位数；

3)位错能够抑制钚的自辐照损伤，其可能的机制是：位错的存在引起晶格形变，在位错线两侧分别产生一个拉伸应力区和一个压缩应力区。拉伸引力区能够吸收间隙原子，压缩应力区能够吸收空位，从而使得位错向点缺陷方向发生攀移，形成一个割阶。

这些结论有助于从原子尺度理解δ相钚-镓合金的自辐照损伤机理，并为将来更大尺度的计算模拟提供重要输入。需要注意的是，本文设置的初级离位原子最大能量仅为5 keV，而钚自辐照产生的初级离位原子能量理论上最高能够达到64 keV，下一步还需要对更高能量的自辐照损伤开展研究。