金属钨级联碰撞中势函数的影响

刘丽霞，陈阳春，邱荣阳，胡望宇,*，邓辉球

(1.湖南大学 材料科学与工程学院，湖南 长沙 410082；2.湖南大学 物理与微电子科学学院，湖南 长沙 410082)

金属钨(W)材料因其具有高熔点、高热导性和良好的抗溅射性等优良特性，已被认为是用于高温和辐照应用最具前景的第一壁候选材料[1-4]。在聚变反应堆中，14.1 MeV能量的中子辐照是面向等离子体材料所要承受的主要考验。当高能中子轰击到第一壁材料表面，与材料发生剧烈碰撞，会在一定程度上改变材料内部的微观结构，产生一系列辐照缺陷，如点缺陷、缺陷团簇、位错环、孔洞等，从而影响材料的宏观性能，如辐照硬化、脆化、肿胀、蠕变等[5-7]。

在研究辐照模拟时，由于高能强流中子难以直接实现辐照模拟，通常以具有一定动能的初级碰撞原子(PKA)，即被中子撞击后获得一部分能量的晶格原子来撞击材料的点阵原子模拟级联碰撞过程，通过分析碰撞级联之后材料发生的变化来分析其辐照损伤情况及其抗辐照性能。由于实验装置的限制及科技的发展，计算模拟在研究材料的抗辐照性能方面发挥着越来越重要的作用。目前被广泛应用的计算模拟方法通常有第一性原理(FP)、分子动力学(MD)、动力学蒙特卡罗(KMC)、团簇动力学(CD)和有限元方法(FEM)等，其中分子动力学方法因其可模拟整个碰撞级联过程中每个原子的动态演化过程，且所模拟的时间和空间尺度与中子辐照初期级联碰撞的尺度相一致，而成为模拟级联碰撞的首选方法。近年来，在缺乏高能中子装置进行辐照实验的情况下，尽管国内外在聚变反应堆材料模拟的中子辐照方面进行了许多工作，但仍缺乏系统且完善的模拟平台来实现聚变材料基于中子辐照下从微观到宏观尺度的变化。因此，实施分子动力学研究来建立完善的级联碰撞数据库，为后续的模拟提供输入参考，从而实现不同模拟尺度间信息的有效传递是非常有必要的，将会促进核聚变反应堆中钨材料的发展，具有重要的潜在应用价值。

Nordlund[8]回顾了近年来材料辐照效应计算机模拟的历史，详细介绍了辐照模拟的研究方法，充分阐明了分子动力学模拟的重要性。Fikar等[9-10]使用分子动力学方法研究了不同钨势函数的级联模拟，发现尽管不同势函数的离位阈能不同，各势函数级联碰撞之后会产生相似的结果，即稳定状态下的总缺陷数目相差不大，但缺陷的空间分布有一定的差异，这可能与间隙迁移能有关。Setyawan等[11]使用LAMMPS[12]代码通过分子动力学方法研究了钨中高能入射PKA在不同温度梯度下(300～2 050 K)的级联碰撞模拟，结果表明稳定状态幸存的缺陷数目对PKA能量有较强的依赖性，而对温度的依赖性很小；Yang等[13]研究了金属钨的级联碰撞模拟，结果表明在稳定状态幸存的缺陷数目虽然随辐照温度的升高而略有减少，但对温度的依赖效果并不显著；Warrier等[14]研究了面心立方铜和体心立方钨两种材料中的级联碰撞模拟，发现在1～5 keV的PKA低能区间内，当进行使离位原子数目稳定在平均值附近的模拟次数统计时，在钨体系中比在铜体系中所需统计的次数要少，即钨体系中更快达到平衡。Fellman等[15]利用分子动力学方法模拟了钨中含空位团簇和孔洞的级联碰撞模拟，结果表明当级联发生区域与空位团簇所在区域重叠时，最后产生的稳定状态下的缺陷数目会减少，这与Byggmästar等[16]研究的钨中含间隙缺陷团簇的级联碰撞模拟结果相似，Byggmästar等发现当级联发生区域与间隙缺陷团簇所在区域重叠时，级联产生的最终稳定状态下的缺陷数目有所下降，且当缺陷团簇的尺寸与级联区域体积的大小相当时，级联产生的新增缺陷对数量几乎为0。Fu等[17]使用一个新开发的WRe合金势函数研究了钨及钨铼合金中高能PKA的级联模拟，发现缺陷团簇的尺寸和数目均随PKA能量的增加而增加，且纯钨中间隙团簇和位错环的移动性比钨铼合金中要高。Zhang等[18]使用分子动力学方法进行了钨中在晶界附近不同温度下的级联碰撞模拟，结果表明空位缺陷的数目与温度无关。

目前已有大量使用分子动力学方法模拟钨基材料级联碰撞的研究，且所使用的原子间势函数大多各有优缺点，为了使其更好地用于辐照损伤模拟，钨势函数还在不断地进一步优化改进。随着势函数的不断发展，近几年来出现了一些新的优化版本的钨基势函数。不同原子间势函数在辐照级联模拟中的影响有待进一步研究，尤其是不同势函数在级联碰撞过程中辐照点缺陷、缺陷团簇和位错环等分布信息，可为钨基材料初级辐照损伤的理解及退火过程缺陷长时间的演化模拟提供基础，为用于辐照模拟钨势函数的选择和优化提供参考。

本文对金属钨的势函数进行比较测试，分析辐照过程中缺陷的产生、缺陷团簇和位错环分布信息。

1 不同势函数的中子辐照级联模拟

折衷考虑到势函数预测准确性和在大尺度计算时的效率，主要选取的势函数类型有Embedded-atom方法(EAM)[19]和Finnis-Sinclair(F-S)形式[20]两种类型。首先选取了两个目前为止使用得较频繁的势函数：一个是由Ackland等[21]开发的F-S形式的钨势函数(势函数AT)，另一个是Juslin等[22]在势函数AT的基础上进行了优化和修正的钨势函数(势函数JW)。然后选取了最近发表的几个各有优势的钨基合金势函数：Marinica等[23]发表了3个EAM势函数(势函数MV2、MV3和MV4)，其中的势函数MV2被Bonny等[24]进行了更新和优化，开发了一个新的WRe合金势函数，即势函数MV2-B；势函数MV4由Setyawan等[25]进行了优化，开发了另一个新的WRe合金势函数，即势函数MV4-S；Chen等[26]新发表了一个可准确预测钨中位错环形成能的F-S形式WRe合金势函数(势函数Chen)。本工作选取上述5个势函数进行动态级联碰撞模拟测试。

为更好地模拟在级联过程中原子间的短程相互作用，首先对各势函数进行Ziegler-Biersack-Littmark(ZBL)修正，并计算了其离位阈能，具体的连接参数和形式参考文献[27]。级联模拟采用三维周期性边界条件，模拟晶胞的最外层6层原子(3×a0，a0为0.316 52 nm)设置为边界区域，如图1所示Region-Ⅰ，模拟均在300 K温度下进行。在辐照级联模拟启动前，首先使用共轭梯度法将系统进行静态驰豫，然后整个模拟晶胞均通过Nose-Hoover热浴和控压方法在等温等压(NPT)系综下进行动态弛豫20 ps，时间步长设置为1 fs。弛豫充分之后，为避免沟道效应，且使得级联碰撞在盒子中心区域产生，参考文献[17，28-29]，在模拟晶胞的中心附近随机选取1个钨原子作为PKA原子，沿典型高指数〈135〉晶格方向启动级联碰撞模拟，一旦超出模拟盒子边界，则该模拟结果数据无效。在整个级联碰撞过程中，模拟晶胞外层边界区域(Region-Ⅰ)通过速度标定法来控制该区域晶格体系温度保持为300 K，内部区域(Region-Ⅱ)则在微正则(NVE)系综下自发进行。级联过程使用变时间步长，且控制在10-3～10-7ps之间，模拟的总时间根据PKA能量(EPKA)大小的不同选择为20～60 ps。具体的模拟时间、模拟盒子的边长和模拟事件总数列于表1。由于级联碰撞过程具有较大的随机性，为降低误差，本文采用随机选择不同PKA原子的方式对所有势函数和所模拟PKA能量下的每种情况都进行了15次模拟级联碰撞过程。本文所使用的分子动力学模拟软件为LAMMPS[12]，辐照级联后产生的点缺陷分析方法为Winger-Seitz(W-S)原胞方法[30]，而位错分析提取为DXA位错算法[31]。此外，这两种方法均可通过可视化和数据分析软件OVITO[32]来实现。

图1 模拟盒子示意图Fig.1 Schematic of simulation box

表1 势函数级联碰撞模拟参数Table 1 Collision cascade simulation parameter for potential

2 缺陷的产生

表2列出不同势函数在PKA能量为10 keV和50 keV级联过后稳定状态下弗兰克尔缺陷对(FPs)的数目。由表2可见，尽管各势函数的离位阈能有差别[27]，但在PKA能量为10 keV下，各势函数FPs数目的差别很小，这与之前Fikar等[9-10]获得的结果一致；在PKA能量为50 keV下，不同势函数之间的FPs数目略有一些差别，如果增加模拟次数其差异可能会缩小。此外，FPs的数目随着PKA能量的增加而显著增加。

图2示出各势函数FPs数目在PKA能量为50 keV下随时间的演化，表3列出各势函数达到热峰状态的时间和演化趋于平衡的时间。结果表明，50 keV的PKA能量下，各势函数到达级联热峰的时间(T热峰)和趋于平衡状态的时间(T平衡)均相差不大，均在约1 ps产生最多的缺陷，达到热峰状态，随着级联的演化，间隙与空位缺陷快速湮灭复合，且在10 ps内趋于平衡，缺陷状态基本稳定，各势函数间最后稳定状态下的FPs数目差别不大。

表2 不同势函数在稳定状态下的FPs数目Table 2 Number of FPs in stable state obtained by different potentials

图2 不同势函数下FPs数目随时间的演化Fig.2 Number of FPs as a function of simulation time obtained by different potentials

表3 不同势函数下级联到达热峰和平衡稳定状态的时间Table 3 Time for cascade simulation with different potentials to reach peak state and equilibrium stable state

红色和蓝色小球分别代表钨间隙原子和空位a，c——热峰状态；b，d——最终稳定状态图3 级联碰撞两种典型的缺陷分布Fig.3 Two typical defects distributionswith collision cascade

此外，注意到级联热峰阶段缺陷分布的空间形貌主要有两种：一种是集中型，其热峰状态构型如图3a所示，另外一种是较为分散的连续型次级联结构，如图3c所示。相对而言，热峰阶段为集中型时，随着级联的后续演化，间隙和空位的湮灭复合后，稳定状态级联中心可能会形成略大的空位核区域，由较大的间隙团簇和单间隙原子包围，如图3b所示，且最后稳定状态下的FPs数目会略多于分散型形貌，缺陷团簇的尺寸也会更大，总体上可能会使其拥有更高的团簇分数。热峰阶段为较分散的连续型次级联结构时，最后稳定状态的缺陷多以单间隙或单空位以及小尺寸的团簇形式存在，很少能观察到较大的缺陷团簇，如图3d所示。这与Fu等[17]在钨及WRe合金高能中子辐照中观察到的结果一致，且此类现象在高能PKA级联时会更显著。

图4示出各势函数在PKA能量为50 keV下的级联过后稳定状态的缺陷构型。由图4可见，间隙缺陷多以单间隙〈111〉哑铃(dumbbell)形式存在，缺陷团簇数目不多，总体上间隙型缺陷团簇比空位型多。且PKA能量为50 keV下观察到有位错环和位错线的出现，位错环大多是由间隙型团簇形成。而在PKA能量为10 keV情况下，5种势函数在稳定状态下的间隙缺陷构型大多数为单间隙〈111〉哑铃，团簇数目较少，可忽略。本工作中，定义包含2个及2个以上的净缺陷数(间隙原子/空位)为缺陷团簇。在已有文献[17，33-35]中，间隙为第3近邻、空位为第2近邻的缺陷团簇截断距离判据已被广泛使用。本文以势函数Chen为例，计算了其稳定状态下间隙与空位缺陷的径向分布函数g(r)，如图5所示。由图5可见，级联过后对于间隙原子，主峰在第3近邻(NN3)处，而对于空位，主峰在第2近邻(NN2)处。基于以上结果，本文中缺陷团簇截断距离判据选取间隙为第3近邻、空位为第2近邻。

总体而言，尽管各势函数的离位阈能略有差别[27]，但最后稳定状态下的FPs数目的差别很小，且各势函数FPs数目随PKA能量的增加均显著增加。

3 缺陷团簇分布

红色和蓝色小球分别代表钨间隙原子和空位，间隙以哑铃形式显示图4 势函数PKA能量为50 keV下稳定状态的缺陷构型Fig.4 Defect configuration in stable state obtained by different potentials at PKA energy of 50 keV

图5 稳定状态下间隙和空位的径向分布函数Fig.5 Radial distribution function of interstitial and vacancy in stable state

表4 不同势函数下缺陷团簇分布情况Table 4 Distribution of defect cluster obtained by different potentials

图6示出各势函数PKA能量为50 keV下的团簇分数。由图6可见：所有势函数的间隙团簇分数均大于空位团簇分数，这表明间隙更容易形成团簇；团簇分数较大的势函数AT和MV2-B，在15次级联模拟中热峰阶段缺陷主要呈现集中型形貌(所占比例分别为60%和80%)，与之相比，较为分散的次级联型形貌出现概率略小；其他势函数(Chen、JW、MV4-S)热峰阶段呈现集中型形貌的概率为30%～40%，致使这些势函数拥有略小的团簇分数，且在势函数MV4-S中，间隙和空位的团簇分数相差不大。此现象可能与各势函数的间隙迁移能有关[9-10]，势函数AT、JW、MV2-B、MV4-S和Chen所计算的〈111〉哑铃迁移能分别为0.03、0.01、0.04、0.03、0.13 eV。〈111〉单间隙拥有较小的迁移能会使得间隙缺陷具有更强的运动性，级联后它们能快速逃离级联中心区域，形成更广的缺陷分布范围，降低了与空位复合的概率和形成团簇(尤其是较大尺寸团簇)的可能性，最终使得缺陷团簇分数变小。势函数JW的间隙迁移能比其他势函数略小，其拥有最小的间隙团簇分数；而尽管势函数Chen所计算的〈111〉哑铃迁移能略大于其他势函数，但它在本文所模拟的结果中却未拥有最高的团簇分数，这与在目前总模拟次数中势函数Chen在热峰阶段呈现集中型形貌的概率约为40%有关(低于势函数AT和MV2-B)，继续增加模拟次数可能会有所变化。

图6 势函数PKA能量为50 keV下的缺陷团簇分数Fig.6 Defect cluster fraction obtained by different potentials at PKA energy of 50 keV

图7示出各势函数PKA能量为50 keV下的团簇信息分布。图7中5种势函数的间隙和空位团簇均以小尺寸(包含2～20个间隙/空位)团簇为主。由图7a可见，势函数AT和MV2-B出现稍大尺寸间隙团簇的概率比势函数Chen、JW和MV4-S的要大；图7b表明，所有势函数均以小尺寸的空位团簇为主，略大一些的空位团簇主要在势函数AT和JW及MV4-S中观察到。

从缺陷团簇分布方面来看，各势函数的缺陷团簇均以小尺寸(包含2～20个间隙/空位)的团簇为主；各势函数的间隙迁移能对其级联过程中缺陷的空间分布范围有影响，而缺陷的空间分布与团簇分数有关，热峰时缺陷为集中型形貌时，最终稳定状态下产生缺陷团簇的数目会更多且尺寸会更大，这种情况下获得一个较大团簇分数的可能性更高。

图7 势函数PKA能量为50 keV下的缺陷团簇分布Fig.7 Defect cluster distribution obtained by different potentials at PKA energy of 50 keV

4 位错信息分布

图8示出各势函数PKA能量为50 keV下的位错信息分布。在模拟过程中，5种势函数均观察到有位错线及位错环的产生。根据文献[26]所计算的各势函数位错环的形成能来看，势函数Chen与第一性原理(DFT)中1/2〈111〉位错环的形成能低于〈100〉位错环形成能的趋势一致，且结果最相符。势函数AT和JW也与DFT结果有相似的趋势，但这两个势函数1/2〈111〉和〈100〉位错环之间的形成能差值相对较小，而势函数MV2-B和MV4-S所计算的位错环形成能的结果与DFT结果不一致。势函数Chen、AT、JW和MV4-S 4种势函数级联后所产生的位错环以1/2〈111〉间隙型位错环为主，与实验上观察到的结果一致[36]。而势函数MV2-B以〈100〉位错环出现较多，1/2〈111〉位错环略少，这与势函数MV2-B所预测的位错环形成能为〈100〉位错环低于1/2〈111〉位错环有关。此外，势函数AT和JW分出现了1次〈100〉间隙型位错环和空位环，而势函数Chen和MV4-S目前未观察到〈100〉位错环，这可能与模拟次数的限制有关，也可能与位错环形成能有关。势函数Chen所预测的〈100〉位错环的形成能高于〈111〉位错环，所以势函数Chen的模拟结果观察到的位错环以〈111〉类型为主。在本文PKA能量为50 keV下模拟结果所观察到的位错环以小尺寸的为主，而MV4-S势函数在小尺寸下(原子数<20)亦是1/2〈111〉位错环的形成能低于〈100〉位错环，随位错环尺寸的增大，则为1/2〈111〉位错环的形成能高于〈100〉位错环，所以势函数MV4-S在小尺寸位错环下可能会以〈111〉类型为主，而大尺寸位错环更可能为〈100〉类型。

图9示出各势函数PKA能量为50 keV下位错环的尺寸和数量分布。由图9a可见，在势函数Chen、AT和JW中占主导地位的为1/2〈111〉小尺寸位错环，略大尺寸的1/2〈111〉位错环主要在势函数AT中出现。由图9b可见：〈100〉位错环主要在势函数MV2-B中出现，且以间隙型位错环为主；势函数AT中观察到1次小尺寸的间隙型〈100〉环，势函数JW中观察到1次大尺寸的空位〈100〉环。

从位错信息分布来看，势函数Chen、AT、JW和MV4-S 4种势函数级联后所产生的位错以伯格斯矢量为1/2〈111〉间隙型位错环为主，与实验结果一致；而势函数MV2-B以〈100〉位错环出现较多，1/2〈111〉位错环略少，与实验结果略有出入。这与各势函数预测的不同类型的位错环形成能大小有关，与在第一性原理计算中1/2〈111〉位错环的形成能低于〈100〉位错环的结果相比，在这5个势函数中，势函数Chen的计算结果与第一性原理计算结果最相符。

结果为15次重复模拟统计的总和a——1/2〈111〉位错环；b——〈100〉位错环图9 势函数PKA能量为50 keV下位错环的尺寸和数量分布Fig.9 Size and number distributions of dislocation loop obtained by different potentials at PKA energy of 50 keV

5 小结

本文对5个典型的钨势函数进行了PKA能量在10 keV和50 keV下的中子辐照级联碰撞模拟，系统分析和讨论了级联碰撞过程中缺陷的产生与分布、缺陷团簇和位错环的数目与结构等信息，获得的结果对于钨基材料初级辐照损伤的理解以及退火过程缺陷长时间的演化模拟提供了基础，为用于辐照模拟钨势函数的选择和优化提供了参考。得到的主要结论如下：

1) 对于级联碰撞到达热峰状态的时间及平衡稳定状态下FPs数目，不同势函数的模拟结果没有明显差别。

2) 在所模拟的PKA能量下，不同势函数产生的缺陷团簇均以小尺寸(净缺陷数<20)为主，对于势函数AT和MV2-B，缺陷团簇分数比其他势函数较高，且较易出现大尺寸缺陷团簇，这与其热峰阶段易呈现集中型缺陷形貌有关。

3) 从位错环的分布来看，势函数Chen、AT、JW和MV4-S级联碰撞模拟后所产生的位错环以1/2〈111〉间隙型位错环为主，与实验结果一致，MV2-B势函数出现的〈100〉间隙位错环比1/2〈111〉间隙位错环略多，与实验结果不符。