陈丹丹，贺新福，杨 文，储根深，白 鹤，胡长军,*

(1.北京科技大学，北京 100083；2.中国原子能科学研究院，北京 102413)

核反应堆材料的服役性能主要取决于辐照产生的缺陷的动力学行为，如缺陷的扩散、形核、长大、湮灭等，这些长时间尺度的演化行为会引起缺陷尺寸和数密度分布的变化，从而导致核材料服役性能的退化，如辐照脆化、辐照肿胀等。实体动力学蒙特卡罗(object kinetic Monte Carlo, OKMC)方法[1-4]和基于平均场速率理论发展来的团簇动力学(cluster dynamics, CD)方法[2,5-8]是两种广泛使用的研究上述问题的方法。前者通过追踪每个缺陷的随机扩散和相互作用模拟缺陷的演化；后者则假设缺陷处于各向同性的均匀介质中，通过缺陷反应的速率方程跟踪缺陷浓度的变化。然而，二者均有其优点和局限性。OKMC可模拟多种缺陷间的复杂行为，且可捕获缺陷间的空间相关性，但常受限于缺陷间复杂行为模拟所需的计算时间和计算量。CD方法一般具有很高的计算效率，可模拟很高的辐照剂量和时间尺度。然而，这种方法通常局限于模拟含少量可动缺陷的体系，很难处理复杂体系的缺陷演化。

空间分辨随机团簇动力学(spatially resolved stochastic cluster dynamics, SRSCD)[9-12]是近年来发展的一种模拟核材料辐照缺陷行为的新方法。它基于随机团簇动力学(stochastic cluster dynamics, SCD)[13]发展而来，将模拟体积划分为多个体积元，假设每个体积元内部缺陷均匀分布，可聚集和分解，体积元之间存在浓度差及其导致的扩散。缺陷间各反应的反应速率由经典的CD方法推导而来，反应的选择和时间增量则由经典的KMC算法确定。SRSCD方法一方面避免了CD在缺陷种类和行为复杂性方面的限制，另一方面减少了与OKMC相比的计算需求。

为扩大SRSCD的模拟体积，并解决扩大体积后带来的计算量，最有效的方式就是并行处理。在确定性方法中，并行区域按相同的时间步长向前推进，而在KMC方法中，由于每步的时间增量是基于系统中所允许的反应的加权随机选择的，可能会导致并行区域的异步推进。为解决这个问题，Martinez等[14]提出了一种同步并行KMC算法。该算法通过在经典KMC算法中增加选择空事件的可能性，为所有MPI(message passing interface)进程选择同一个时间增量，从而将并行区域同步向前推进。Dunn等将这种同步并行KMC算法应用在SRSCD的并行实现中，提出了同步并行SRSCD方法[12]。然而，由于其一维的进程拓扑方式(实际模拟区域为三维)以及使用了大量的阻塞式点到点通信，使得程序的并行效率较低，很难进行大规模的SRSCD模拟。

本文采用同步并行SRSCD方法，开发用于模拟核材料辐照损伤的大规模并行SRSCD程序——MISA-SCD1.0，并将其用于模拟反应堆压力容器(reactor pressure vessel, RPV)钢模型合金中富Cu团簇的析出，以再现富Cu团簇在辐照条件下的演化过程，验证程序的正确性并分析程序的并行性能。

1 SRSCD方法

SCD是CD的一种随机变体[15]，通过将平均场CD限制在有限体积V内来将缺陷浓度的演化方程转变为在有限体积V中演化整数值的缺陷数量Ni。

(1)

其中：i、j、k分别为不同的缺陷类型；G、R、K分别为体积V内的0阶反应(缺陷产生)、1阶反应(分解、被阱吸收而湮灭)、2阶反应(聚集、复合)的反应速率。

SCD通过随机抽样模拟有限体积内的缺陷演化过程，从而避免大量微分方程的组合展开，非常适合于处理复杂缺陷团簇的演化问题。若考虑多个这样的体积元，再加上不同体积元间的缺陷浓度差异以及扩散，则随机团簇动力学可具有空间分辨率，即SRSCD。

在SRSCD中，认为每个体积元内的缺陷均匀分布，缺陷可聚集和分解，可动缺陷则可在体积元之间扩散，缺陷间的各反应的反应速率由经典的CD方法推导而来，时间增量Δt和反应μ的选择则由经典的KMC算法确定。结合Marian等[13]和Dunn等[9]的工作，给出SRSCD方法中各类反应速率的计算如下。

1) 0阶反应速率

在辐照损伤的情况下，0阶反应被用来表示初始缺陷或其他注入物的来源。通常其反应速率用辐照损伤速率或其他注入物的注入速率表示。

2) 1阶反应速率

1阶反应包括分解(即发射点缺陷)和缺陷被阱吸收而消失。其中，分解反应的反应速率R为：

(2)

被阱吸收而消失的反应速率为：

R=SDiNi

(3)

其中：S为阱强度；Di为缺陷i的扩散率。

3) 2阶反应速率

2阶反应项解释了两种缺陷i和j相互碰撞的各种机制，如间隙团簇(空位团簇)对自间隙原子(空位)的吸收、间隙团簇与空位团簇的复合等，其反应速率R的计算如下：

(4)

对于2阶反应，考虑3种情况：缺陷i和缺陷j均是三维扩散(如空位团簇)；缺陷i和缺陷j中，1个是三维扩散，1个是一维扩散(如位错环)；缺陷i和缺陷j均是一维扩散。

4) 扩散速率

在SRSCD中，缺陷的扩散由体积元之间的浓度梯度驱动。假设扩散均匀、恒定，则缺陷i在体积元p、q之间扩散的反应速率R为：

(5)

2 MISA-SCD1.0实现

2.1 MISA-SCD1.0概述

MISA-SCD1.0的计算流程如图1所示，主要包括区域分解及预处理、KMC循环、后处理3部分。

图1 MISA-SCD1.0的流程图Fig.1 Flow chart of MISA-SCD1.0

区域分解及预处理部分首先根据进程数将模拟区域均匀划分，利用MPI笛卡尔拓扑建立三维进程拓扑结构，并将相应的子模拟区域及其网格映射到各进程上，并为每个进程构建Ghost区域，用于存储邻居进程上与当前进程邻接的网格，如图2所示，其中的网格即为上述体积元。对于每个网格(或进程)，还需建立其上、下、左、右、前、后6个方向的关联关系，即存储其6个邻接网格(或进程)编号，这6个方向也即缺陷扩散的6个方向，如图3所示。

图2 MISA-SCD1.0的区域分解示意图Fig.2 Domain decomposition diagram of MISA-SCD1.0

图3 网格6个邻接方向Fig.3 Six adjacent directions of mesh

KMC循环为整个程序的核心计算部分，采用同步并行KMC算法实现，各进程在各自负责的子区域块上独立选择反应并发生，并在需要时进行通信，以同步缺陷信息。

后处理部分用于统计并输出模拟的中间结果和最终结果。在后处理中，统计各子区域上每个网格中的缺陷类型及其数量并输出。

2.2 同步并行KMC算法

(6)

图4 时间异步推进示意图Fig.4 Schematic diagram of time asynchronous advance

D个区域中的最大总反应速率为：

(7)

(8)

(9)

(10)

2.3 同步通信策略

MISA-SCD1.0中的通信主要发生在更新缺陷过程，对于1个反应，当反应物(参与反应的缺陷)和/或产物(反应生成的缺陷)所处的网格为Boundary区域(位于当前进程且与邻居进程相邻)的网格和Ghost区域的网格时，需与邻居进程进行点对点通信，以同步反应物和/或产物的信息，图5示出了6种需通信的情况。

在一般操作中，将Boundary区域和Ghost区域的缺陷分开发送/接收，且发送/接收前，需通信待更新的缺陷个数，以确定接收缓冲区的大小，则需进行2×(6×2+6×2)=48次点到点通信(即分别发送和接收6个方向的Boundary和Ghost两个区域的缺陷信息)。在MISA-SCD1.0中，采用计算与通信重叠、通信合并的方式进行通信优化，如图6所示。由于1个进程在1个时间步内只发生1个反应，且允许的反应类型是已知的，因此1个时间步内待更新缺陷数目的最大值已知，直接将接收缺陷的缓冲区设置为该值，则省去缺陷个数的通信操作，通信次数减为2×(6×2)=24次。进一步将Boundary区域和Ghost区域的缺陷合并发送，通信次数最终减为2×(6×1)=12次。

图5 缺陷更新过程中的通信情况Fig.5 Communication during defect update

图6 通信优化过程Fig.6 Communication optimization process

2.4 数据结构设计

在进行数据结构设计时，链表能充分利用内存中的碎片空间，在进行插入和删除操作时，比数组更高效，适用于需频繁进行插入/删除操作的应用。由于KMC算法中，模拟的每步均需经过选择反应、更新缺陷、更新反应速率这3步，需频繁更新网格中的缺陷及其可能发生的反应，且缺陷和反应的数量是动态变化的。因此，在MISA-SCD1.0中采用链表来存储频繁更新的缺陷及其反应等信息，而对于网格、进程、输入参数等需频繁查找的信息则采用数组存储。MISA-SCD1.0中，为每个网格创建1个缺陷列表(defectList)和1个反应列表(reactionList)，如图7所示。缺陷列表中存储每种缺陷的属性(包括缺陷类型、数量等)，反应列表中则存储该网格内可能发生的反应及其速率(包括发生反应的反应物数量及类型、产物的数量及类型、反应物及产物所在的网格局部编号、反应物及产物所在的进程编号、反应速率等)。

图7 每个网格的缺陷列表和反应列表Fig.7 Defect list and reaction list for each mesh

3 实验结果及分析

3.1 正确性评估

RPV是压水堆中堆寿期唯一不可更换的核心设备，其服役性能直接影响反应堆的安全性。辐照脆化是RPV钢面临的主要性能问题，一直是核材料领域关注的热点。目前大量研究结果表明，富Cu团簇的析出是导致RPV钢脆化的主要原因之一。近年来，研究者们选取Fe-Cu合金作为RPV钢的模型合金，利用模拟和实验的手段，广泛研究了富Cu团簇的析出机理[6,16-21]。

为验证MISA-SCD1.0的正确性，并评估其模拟RPV钢中缺陷演化的能力，本节选取了不同Cu含量的Fe-Cu合金，分别模拟了电子辐照和中子辐照下的富Cu团簇析出过程，并与实验结果和其他类似的模拟结果进行了对比分析。

1) RPV钢模型合金中富Cu团簇析出模拟

MISA-SCD1.0进行RPV钢模型合金(Fe-Cu合金)中富Cu团簇析出模拟时，允许的缺陷类型及反应如下。

(1) 点缺陷，包括自间隙原子(I)、空位(V)、Cu原子(Cu)。这类缺陷被视为球形，可在三维空间中扩散，这类缺陷可形成团簇，也可与其反类型的缺陷相互湮灭。

(2) 空位团簇，尺寸大于等于2的空位团簇为不可动的球形团簇。空位团簇可吸收单空位、自间隙原子，也可发射单空位。空位团簇也可与Cu原子结合形成Cu_Vac团簇。

(3) 位错环，尺寸大于等于2的自间隙团簇被视为不可动的位错环(环形)。位错环可吸收单空位和自间隙原子，也可发射单空位。位错环不与Cu原子结合。

(4) Cu团簇，包括只含Cu原子的纯Cu团簇，以及同时含有Cu原子和空位的Cu_Vac团簇。这类团簇视为不可动的球形团簇，可发射、吸收单空位及Cu原子。

所有允许的反应类型前面已描述，计算反应速率所需的迁移能Em、扩散前置因子D0、结合能Eb，列于表1。其中，尺寸为3～16的空位团簇的结合能参数以及Cu_Vac团簇与Cu原子或空位的结合能参数来自Kulikov等的BCC Fe-Cu研究[16]。对于大尺寸的空位团簇(n>16)和自间隙团簇(n>2)，其结合能Eb(n)由下式给出：

(11)

其中：n为团簇尺寸；Ef为点缺陷的形成能；Eb(2)为双空位或双间隙团簇的结合能。

对于大尺寸(n>2)的纯Cu团簇的结合能则由下式给出：

Eb(n)=Ω-TΔS-

(12)

其中：Ω为Cu在Fe中溶解时的焓变，Ω=6 255kB[6]；ΔS为非构型熵，ΔS=0.866kB[6]；σ为表面能，为0.37 J/m2[20]。

表1 空位、自间隙和Cu的迁移参数和结合能参数Table 1 Migration and binding parameters for vacancies, self-interstitials and Cu

2) 电子辐照Fe-1.34at.%Cu中Cu析出模拟

本实验中，使用MISA-SCD1.0模拟电子辐照(2.5 MeV)下Fe-1.34at.%Cu中的Cu析出过程。辐照温度为290 ℃，剂量率为2×10-9dpa/s，模拟体积为400 nm×400 nm×200 nm，网格分辨率为10 nm，其他模拟参数列于表2，与Christien等[6]采用的参数相同。

表2 290 ℃电子辐照下Fe-1.34at.%Cu合金中Cu析出模拟的参数Table 2 Parameters for simulating Cu precipitation in Fe-1.34at.%Cu alloy under electron irradiation at 290 ℃

MISA-SCD1.0模拟获得的Cu团簇的数密度及平均半径随辐照剂量的变化，与类似的模拟结果以及实验结果的对比如图8所示，图中，红色点线为MISA-SCD1.0的模拟结果，蓝色实线是Christien等[6]采用传统CD方法获得的模拟结果，黑色圆圈是Mathon等[22]用小角中子散射(small-angle neutron scattering, SANS)的表征技术测量的实验结果。由于SANS测量团簇半径的检测极限约0.5 nm，因此，本文仅将含有10个以上Cu原子(半径约0.3 nm)的团簇(忽略Cu_Vac团簇中的空位)记入数密度和平均半径的统计中。虽然选取的统计阈值略小于SANS的检测极限值，但Bai等[20]的工作表明，此阈值的选取并不会影响Cu团簇后期结果的评估。

如图8a所示，在辐照的初始阶段，Cu团簇的总数密度迅速增加，在4 s左右达峰值，约5×1025m-3。这期间，Cu团簇的平均半径缓慢增加(图8b)。表明此期间主要是Cu团簇的形核和长大阶段。随后Cu团簇的演化进入粗化阶段，平均半径不断增大，数密度逐渐减小。整体上，MISA-SCD1.0的结果与文献中传统CD方法的模拟结果以及实验结果非常吻合。

3) 中子辐照Fe-0.3at.%Cu中的Cu析出模拟

为验证MISA-SCD1.0模拟中子辐照下的缺陷演化能力，进行中子辐照(能量大于1 MeV)Fe-0.3at.%Cu中的Cu析出过程模拟。辐照温度为300 ℃，剂量率为1.4×10-7dpa/s，模拟体积同样为400 nm×400 nm×200 nm，网格分辨率为10 nm，位错密度为5×10-5nm-2，辐照增强因子为7.7×105，其他模拟参数列于表2。

在中子辐照下，初始缺陷不再像电子辐照那样只有点缺陷，而是包含点缺陷和它们的小团簇。在使用传统CD方法的中子辐照模拟中，普遍对这种级联缺陷(即包含点缺陷和小团簇的初始缺陷)进行均匀化处理，以级联效率的方式引入，而在SRSCD中，则可直接将级联缺陷引入到模拟体系中。Bai等[20]的模拟结果表明，级联效率为0.4、空位迁移能为1.0 eV时，获得的模拟结果与实验测量结果最接近。而曹晗等[23]的分子动力学模拟结果表明，级联效率为0.4时对应的初级离位原子(primary knock-on atom, PKA)的能量约为20 keV。因此，在使用MISA-SCD1.0进行本次模拟时，引入20 keV的PKA产生的级联缺陷，作为中子辐照的模拟的缺陷产生项。

对于模拟获得Cu团簇，同样只统计了含有10个以上Cu原子的团簇，获得的中子辐照下的Cu团簇的数密度及平均半径随辐照剂量的变化如图9所示。其中，红色点线为MISA-SCD1.0的模拟结果，蓝色实线是Bai等[20]采用传统CD获得的模拟结果，黑色圆圈是Meslin等[21]的SANS实验测量结果。如图9a所示，相较于电子辐照，中子辐照下的Cu团簇粗化阶段发生在较高的辐照剂量范围内，且Cu团簇的平均半径较电子辐照下的半径小(图9b)。整体可看出，MISA-SCD1.0模拟获得Cu团簇数密度及平均半径，与CD的模拟结果和实验结果均高度一致。与电子辐照类似，中子辐照下MISA-SCD1.0模拟得到的Cu团簇数密度同样较CD模拟得到的结果略大，这主要是因为与平均场相比，MISA-SCD1.0的模拟限制在1个较小的体积内。

3.2 性能分析

对MISA-SCD1.0的性能测试在天河2号超算平台上进行，仅使用了CPU核，相关配置参数列于表3。

图8 290 ℃电子辐照下Fe-1.34at.%Cu合金中Cu团簇的数密度(a)和平均半径(b)随辐照剂量的变化Fig.8 Irradiation dose-dependent evolution of number density (a) and mean radius (b)of Cu clusters in Fe-1.34at.%Cu during electron irradiation at 290 ℃

图9 300 ℃中子辐照下Fe-0.3at.%Cu合金中Cu团簇的数密度(a)和平均半径(b)随辐照剂量的变化Fig.9 Irradiation dose-dependent evolution of number density (a) and mean radius (b) of Cu clusters in Fe-0.3at.%Cu during neutron irradiation at 300 ℃

表3 天河2号超算平台的相关参数Table 3 Parameters of Tianhe-2 supercomputer

为测试MISA-SCD1.0采用2.3节中的通信策略的性能提升效果，选取3.1节中的电子辐照Fe-1.34at.%Cu进行实验。固定每核的模拟体积为100 nm×100 nm×100 nm，网格分辨率为10 nm，模拟缺陷演化0.01 s。在24、48、96、192核上分别测试了通信优化前后的程序运行时间，如图10所示。可看到，优化后的计算速度有较明显提升，总体性能较优化前提升约30%，表明2.3节中的通信策略在进行大规模模拟时能带来良好的计算加速。

图10 MISA-SCD1.0通信优化前后的运行时间对比Fig.10 Running time comparison before and after communication optimization in MISA-SCD1.0

为测试MISA-SCD1.0大规模并行模拟的性能，同样选取Fe-1.34at.%Cu在290 ℃下的电子辐照模拟，在天河2号超级计算机上，从强、弱可扩展性两方面对MISA-SCD1.0进行测试分析。其中，强、弱可扩展的并行效率由下式给出：

(13)

其中：ηstrong和ηweak分别为强、弱可扩展性的并行效率；T1和Tn分别为程序在基准核数和n核上的运行时间；N为n核相对于基准核数的倍数。

在强可扩展性测试中，固定总模拟体积为4 800 nm×2 400 nm×2 400 nm，网格分辨率为10 nm，分别测试了MISA-SCD1.0在1 200、2 400、4 800、7 200、9 600核上，模拟缺陷演化0.001 s的程序运行时间，并以1 200核为基准，计算程序的并行效率，如图11所示。可看到，随着核数的增加，程序的整体运行时间逐渐减少，从1 200核扩展到4 800核的运行时间下降较快，表现出较好的并行效率，达70%以上。当核数进一步增加时，运行时间下降放缓，并行效率下滑。这主要是因为随着核数的增加，通信开销逐渐增加，且单核的计算规模减小，通信占比增大，使得整体的并行效率呈下降趋势。但整体上，从1 200核扩展到9 600核，并行效率维持在59%以上，表现出良好的强可扩展性。

图11 MISA-SCD1.0的强可扩展性测试Fig.11 Strong scalability test of MISA-SCD1.0

在弱可扩展性测试中，采取与强扩展测试相同的模拟对象和模拟时间，固定每核的模拟体积为200 nm×200 nm×200 nm，网格分辨率同为10 nm。随着核数的增加，成倍地扩大总模拟体积。分别测试了MISA-SCD1.0在240、480、960、1 920、3 840核上的运行时间，并以240核为基准，计算程序的并行效率，获得的结果如图12所示。可看到，随着总模拟体积及核数的增加，程序运行时间逐渐增加。这主要是因为，在并行KMC算法中，模拟的每个时间步均需进行全局通信以获得最大总反应速率，核数的增加使得这种全局通信的开销时间逐渐增加，进而导致并行效率逐渐下降。但总体上而言，扩展到3 840核时，并行效率仍在64%以上。图11、12的测试表明，MISA-SCD1.0具有良好的扩展性，能开展大规模的缺陷演化SRSCD模拟。

图12 MISA-SCD1.0的弱可扩展性测试Fig.12 Weak scalability test of MISA-SCD1.0

4 结论

空间分辨随机团簇动力学已发展成为模拟核材料辐照损伤的一种有效方法，且能近似地考虑缺陷的空间信息。本文基于空间分辨随机团簇动力学方法和同步并行KMC方法开发了用于模拟核材料辐照损伤的大规模并行程序MISA-SCD1.0，并介绍了其实现细节，包括程序的计算流程、通信策略以及数据结构设计。通过将MISA-SCD1.0用于模拟不同比例的Fe-Cu合金在电子辐照和中子辐照下的富Cu团簇析出过程，并与传统CD模拟结果以及实验结果进行对比，验证了程序的正确性。在此基础上，进行了程序的强、弱可扩展性测试。在强扩展性测试中，从1 200核扩展到9 600核，并行效率保持在59%以上，在弱扩展性测试中，从240核扩展到3 840核，并行效率保持在64%以上，表明MISA-SCD1.0具有良好的扩展性，能很好进行核材料辐照损伤的大规模并行SRSCD模拟，为材料的宏观性能预测提供有力的支撑。