核电材料辐照损伤的多尺度高通量计算模拟

2022-02-18薛飞王忆刘向兵赖文生季骅刘剑波柳百新

装备环境工程 2022年1期

薛飞，王忆,，刘向兵，赖文生，季骅,，刘剑波，柳百新

（1.苏州热工研究院有限公司，江苏苏州 215004；2.清华大学，北京 100084）

随着核电的发展，核电材料在服役期内的安全性与可靠性已成为人们普遍关心的问题。核电机组堆内具有强中子辐照环境，在运行过程中，包容和防止放射性物质外溢的部件材料与中子相互作用将发生辐照损伤，即产生点缺陷和缺陷团及其演化的离位峰、位错环、原子偏聚、微空洞以及析出相等缺陷。这些辐照缺陷的演化会引起材料性能的变化。例如，对反应堆压力容器用低合金钢来说，中子辐照缺陷的产生会导致辐照脆化；对堆内构件用奥氏体不锈钢来说，中子辐照缺陷的产生会导致辐照加速应力腐蚀开裂。因此，对核电安全来讲，掌握核电材料服役过程中辐照缺陷结构演变的规律及其对材料性能的影响效应至关重要，这也是国内外研究一直关注的焦点问题。目前，在辐照缺陷结构演变实验表征方面取得了一些进展，如已有研究表明(见图1)，反应堆压力容器用低合金钢经中子辐照后损伤缺陷主要有三种：1）空位、空位-溶质原子对等基体缺陷（MD-Matrix Defects）；2）富Cu团簇（CEC-Cu Enriched Clusters）；3）P在晶界的偏聚。这些材料组织结构内部纳米尺度的缺陷，阻碍位错的滑移，从而引起材料的强度、硬度提高，并最终导致材料脆化。

图1 反应堆压力容器中子辐照损伤缺陷[3]Fig.1 The neutron irradiation induced defects in reactor pressure vessel[3]

然而，核电材料寿命预测和性能评价仍面临较大困难。核电材料服役行为评价难度大、费用高、周期长、数据匮乏。而常用的试验堆或离子辐照属于加速辐照环境，与商用堆服役环境相比，其辐照注量率差异大、缺陷累积和材料损伤进程不一致。如何建立不同辐照环境之间的映射关系、确立等效辐照缺陷结构，是当前推动各种辐照损伤研究方法的融合与互补的重要课题。多尺度高通量计算模拟可以预测不同辐照环境中的缺陷累积进程，揭示缺陷演化的机理，在辐照损伤研究领域起关键作用。文中对多尺度高通量计算模拟方法研究核电材料辐照损伤与辐照缺陷演化的发展和现状进行了综述。

1 辐照缺陷结构的演化

核电材料辐照损伤的微观机制是裂变产生的中子注入材料后，通过散射将动能传递给材料原子，造成原子级联损伤、形成辐照缺陷，在长时间的累积作用下，高密度的辐照缺陷将显著影响材料的性能。核电材料在服役期间，中子注入是长期且持续发生的。尽管如此，商用核电的堆内中子注量率并不高，以反应堆压力容器为例，其中子注量率大约处在1×10cms量级。所以在相邻两次的中子注入的间隔期内，辐照缺陷有充足的时间进行新生、分解、迁移、湮灭、复合和聚集等各种反应过程，辐照缺陷的演化不仅是单纯的数量或尺寸上的累积，还包括缺陷结构的变化。辐照缺陷演化过程中存在哪些缺陷结构，它们演化的热力学、动力学规律是什么，这些都是多尺度计算模拟需要研究的关键问题。

1.1 辐照缺陷结构的多尺度层级

中子注入材料后，初级撞出原子在中子散射过程中获得较大动能，被撞离晶格点位、随后与其他原子发生碰撞并耗散掉过剩能量。在初级辐照损伤过程中，材料内部会产生大量的空位和自间隙原子，即Frenkel缺陷。空位和自间隙原子仅占据单晶格点位，引发的畸变较小，具有高度可动性，既能在初级辐照损伤中引入，也有可能在其他大型缺陷分解的过程中被释放出来，是辐照缺陷演化进程中最活跃的缺陷结构类型。

随着空位/自间隙原子的扩散迁移，同种类的辐照缺陷相遇时，它们复合长大，形成空位团/自间隙原子团，而空位与自间隙原子相遇则会发生湮灭效应、部分修复辐照损伤。空位团/自间隙原子团等复合点缺陷具有1 nm左右的尺寸大小，占据多个晶格点位。由于尺寸的成长，其可动性相较于单点缺陷的空位/自间隙原子下降，运动方式也会发生改变，一般较大的自间隙原子团只能进行一维移动。

当复合缺陷成长到一定的尺寸时，缺陷引入的晶格畸变将驱使其坍塌形成位错环。位错环可以进行一维迁移运动，迁移方向由其伯格斯矢量决定。位错环发生分解和复合反应时具有伯格斯矢量的守恒性。在-Fe基的材料中，单个位错环也可能自主地发生伯格斯矢量不守恒的转化反应。空位型位错环具有亚稳特征，它们最终会转化为球形的空洞，从而丧失可动性。位错环及其可能转化成的空洞等辐照缺陷具有1×10~1×10nm的尺度量级，进入了常规显微分析可观测的范围，它们对位错运动造成显著阻碍，弱化了材料的韧性，是核电材料性能降低的主要来源。

上述辐照缺陷结构及其性质见表1，可以看出各种辐照缺陷在结构上具有不同的空间特征，进而导致其可动性及迁移路径各不相同，形成多尺度层级。根据辐照缺陷结构扩散反应机制的不同，必须采用多尺度高通量计算模拟分别进行研究。

表1 辐照缺陷结构的多尺度层级Tab.1 Multi-scale hierarchy of radiation defects

1.2 辐照缺陷结构的特征能量

核电材料内的辐照缺陷在演化过程中可分解成新生、分解、湮灭、复合和迁移等各种微观反应，这些反应都属于热激活过程，因而其反应难易、快慢受控于反应前后的能量变化。依据具体反应类别，可将反应的能量变化分类为缺陷的生成能、结合能和迁移能三大类。

当一个辐照缺陷新生于材料中时，其引入的能量变化称为缺陷的生成能，如下式：

式中：是这个缺陷占据的晶格点位的数量；［+()］是该缺陷内嵌于个基体原子时的系统能量，是所有原子处在原本的理想晶格点位时的参考能量。

当两个或多个辐照缺陷相遇发生复合反应时，其能量变化为辐照缺陷的结合能，如下式：

式中：(D)是第个缺陷D的生成能；(++…+D)是个缺陷复合形成的复合缺陷的生成能。分解反应作为复合反应的逆过程，同样可使用结合能计算其能量变化。

和分别决定了辐照缺陷在新生与复合/分解过程发生的难易程度，反应概率为：

辐照缺陷的迁移在原子尺度上是通过一个或一系列连续的原子跃迁行为链式组合而成的，单原子的跃迁的难易程度由原子的迁移能决定，该能量可计算为：

式中：是跃迁原子位于跃迁路径鞍点位置时的系统能量；是跃迁原子位于跃迁路径初始亚稳定位置时的系统能量。根据过渡态理论，式（4）的迁移能决定了原子跃迁频率：

式中：是原子尝试跳跃频率。

占据多个晶格点位的辐照缺陷的整体迁移过程由一系列原子跃迁的子过程链式组合而成。因此，需要采用一个表观的迁移能来表征迁移的快慢，具体数值可基于扩散理论计算。

上述的生成能、结合能决定了辐照缺陷的热力学性质，而迁移能决定了辐照缺陷的动力学性质。简而言之，只要构建了这几个辐照缺陷结构特征能量的数据库，就可以计算分析辐照缺陷新生、分解、湮灭、复合和迁移等各种微观反应过程的难易、快慢，从而阐明辐照缺陷的演化规律。下文将结合具体的多尺度高通量计算模拟方法介绍相关的发展。

2 辐照缺陷结构演化的多尺度高通量计算模拟

如前所述，在核电材料中，辐照缺陷演化可分解为新生、分解、湮灭、复合和迁移等各种微观反应过程。为了分析这些反应过程对材料结构的损伤效果，评估它们对材料性能的影响，一方面需要分而治之，阐明个体的反应过程的难易/快慢；另一方面，需要统而合之，阐明总体上各反应的耦合作用。因此，多尺度高通量计算模拟的基本框架和研究内容可以归纳至图2。具体包括以下几个主要内容：1）采用第一性原理计算和分子动力学计算来预测分析辐照缺陷结构的热力学稳定性和动力学可动性，形成辐照缺陷的热力学/动力学数据库；2）利用分子动力学计算进行高通量计算，获取辐照的初级辐照损伤数据库；3）在缺陷扩散反应动力学中综合利用各底层尺度计算模拟方法构建的数据库，依据中子的注量率导入相应的初级辐照损伤数据，基于辐照缺陷的热力学/动力学数据评估新生、分解、湮灭、复合和迁移等各种反应过程的耦合作用；4）最终，根据缺陷扩散反应动力学的综合分析得到辐照缺陷的密度和分布等统计性结果，以便预测材料的服役行为。

图2 辐照损伤高通量多尺度计算模拟的研究内容和框架Fig.2 The high-throughput multi-scale simulation framework of radiation damage

2.1 第一性原理计算

第一性原理计算是多尺度高通量计算模拟的基础，具有不依赖人为假设、可靠性强的优点。但是第一性原理计算的模拟成本高、计算效率随尺寸增长呈指数率下降，因而主要适用于小型单点缺陷的研究。表2列出了-Fe中空位（Vac）和自间隙原子（SIA）的生成能，自间隙原子的生成能显著高于空位，因而更倾向于被位错、晶界等吸收以使系统能量降低。此外，当自间隙原子具有<110>型哑铃对结构时，其生成能最低，因此其更容易出现。而在其他的bcc金属中，自间隙原子主要以<111>型挤列子结构存在，只能进行一维迁移；由于特殊的哑铃对结构，-Fe中的自间隙原子可进行三维迁移运动，更容易与其他辐照缺陷相遇并反应，有利于降低空洞生成的概率，但也会增加位错环形成的可能性。

表2 α-Fe中空位/自间隙原子的生成能[11]Tab.2 Formation energy of vacancies/self-interstitial-atoms in α-Fe[11]

核电用钢中普遍含有多种合金溶质，它们对钢材的性能有重要贡献。对于辐照缺陷，溶质的存在会改变原子间相互作用，影响辐照缺陷的生成和迁移。由于溶质种类繁多、作用复杂，溶质对辐照缺陷的影响规律实质上已经成为当前第一性原理计算的主要研究课题。研究发现，在第一近邻的距离溶质普遍与空位存在吸引作用，因而空位更容易出现在溶质附近。而溶质原子一般倾向于排斥自间隙原子，Mn元素是唯一已知的有利于自间隙原子生成的特例。研究表明，溶质对空位及自间隙原子所起到吸引/排斥的强弱与溶质的有效体积因子具有较强的关联性，一般溶质有效体积越大，其吸引/排斥效果越强，如图3所示。从图3a中可以看出以Fe/Ru/Os作为分界点，两侧的元素随着偏离分界点原子序数越远，溶质体积因子变得越大。图3b绘制了溶质-缺陷结合能与溶质体积因子的关系，从中可见，对于同类型缺陷而言，溶质体积因子与结合能之间近似地符合一定的线性关系，不同类型缺陷的线性关系的斜率互不相同。

图3 溶质对空位和自间隙原子缺陷能量影响的统计性规律[13]Fig.3 The statistical law of the solute effects on the formation energies of vacancy and self-interstitials[13]: a) the relationship between the effective volume and the element of solute; b) the relationship between formation energies and the effective volume

在辐照缺陷的可动性方面，第一性原理计算表明自间隙原子的可动性远远大于空位的，在-Fe中，自间隙原子的迁移能一般介于0.3~0.4 eV之间，显著小于空位的迁移能（~0.6 eV）。溶质对辐照缺陷迁移有显著作用，第一性原理计算发现，溶质普遍易于与空位相互拖曳，进行同向的耦合迁移扩散，溶质原子也可能因辐照缺陷与缺陷阱的相互作用而发生辐照诱导偏析。溶质可能降低自间隙原子的迁移速率，这仍有待进一步验证。

2.2 分子动力学计算

基于多体势的传统分子动力学具有计算效率高、时间成本与体系大小几乎呈线性变化的特征，因而适合开展大规模计算模拟。对于空位团/自间隙原子团、位错环等自身尺寸大、影响区也较大的大型辐照缺陷，传统分子动力学计算可以克服第一性原理计算的高成本障碍，是更为现实的计算模拟方法。

式（6）的多体势方法是初级辐照损伤模拟和辐照缺陷热力学/动力学性质研究的主要方法。多体势研究表明，空位团的热稳定性比自间隙原子团差，而空位团/自间隙原子团的迁移能与单空位/单自间隙原子基本相当；溶质易于偏聚在辐照缺陷引起的畸变区附近。这些数据都可传递给缺陷扩散反应动力学以进行缺陷密度的计算模拟。

随着信息技术发展，机器学习与多体势方法的交叉融合有望解决传统多体势难以解决的一些问题。赖文生课题组提出了一种新型机器学习多体势，通过原子构型描述子和人工神经网络相结合来构建新型势函数。其描述子定义如下：

式中：G是基于距离的原子对描述子；G是基于键角的原子三体描述子。通过前述的两类描述子对材料体系内的原子坐标进行处理，转化为人工神经网络的输入层节点，然后导入如图4a所示的多层神经网络来进行自动的学习以获取机器学习多体势。在经过数万组的学习数据的训练后，新型机器学习多体势可以准确重现第一性原理计算的参考数据，如图4b所示。图4c展示了这种机器学习多体势所预测的1200 K温度条件下，由fcc γ-Fe向bcc α-Fe转变的晶体结构，其中大部分原子已转变为理想的bcc结构(深色原子)，但小部分原子仍没转变为bcc结构(浅色原子)。这个预测结果与实验吻合，解决了传统势函数在Fe的相转变问题上面对的困难。由于式（7）定义的描述子可以对任意种类原子的坐标进行转化，新型机器学习多体势也可以扩展到合金中。

图4 机器学习多体势的构建与应用[29]Fig.4 The construction and application of machine-learning potential: a) the neural network adopted by the machine-learning potential; b) comparison of atomic forces predicted by the machine-learning potential and the first-principles calculations; c) the predicted atomic structure for the transformation from γ-Fe to α-Fe

分子动力学有着很高的计算效率，在辐照缺陷的生成和迁移方面有重要意义。但是分子动力学的时间尺度一般限于纳秒量级，而大型辐照缺陷在进行迁移时可能花费秒量级或更长时间。然而，在进行辐照缺陷的迁移的计算时，没有必要对一个连续的时间段进行完全的牛顿力学求解，其原因是原子迁移过程仅仅是原子运动的偶发稀有概率事件，在大部分的时间范围内，原子都在其亚平衡的势阱附近作热振动。加速分子动力学可以避免计算资源耗费在对辐照缺陷迁移没有贡献的热振动上，而集中资源计算原子的跃迁行为。

空位团、自间隙原子团等辐照缺陷的迁移过程由一系列原子跃迁链式组合而成，每次原子跃迁后都会形成新的局域势阱，这些局域势阱互相耦合连接。传统的分子动力学加速方案，如Hyper动力学和Meta动力学方法等只能处理简单势阱体系，不适用于辐照缺陷迁移的研究。为了解决这个问题，近两年来，高宁等发展了自适应加速分子动力学，其原理是通过在势能面上叠加一个根据局域势阱自适应变化的加速势场，从而使辐照缺陷可以跨越连续的势阱实现有效迁移运动，如图5a所示。采用该方法，高宁等人计算了空位-氦气团的扩散系数，其温度依赖关系如图5b所示。

图5 自适应加速分子动力学示意[33]Fig.5 The scheme of the adaptive self-accelerated molecular dynamics: a) theoretical scheme of the self-adaptive accelerated molecular dynamics; b) its application on the He-V migration

辐照缺陷的迁移和反应都在材料内部发生，完全服从核电材料的热力学、动力学规律。辐照缺陷最初的引入过程即初级辐照损伤过程受到外部中子注入过程的影响，不完全服从材料本身的规律，在核电材料缺陷演化的研究中必须要考虑建立初级辐照损伤的数据库。

初级辐照损伤的效果不仅与核电材料有关，还与注入中子的能量和角度有关，在正常的核电运营条件下，注入中子的能谱比较稳定，而注入角度则完全随机、不受核电材料控制，对初级辐照损伤影响更大。第一性原理计算表明，在不同的注入角度上，中子撞击原子使其离位的离位阈能差异非常显著，在<100>和<111>方向上非常容易撞离原子，而在<110>方向则相对困难，如图6a所示。很显然，单独的一次或几次计算模拟不能全面地描述初级辐照损伤现象。对此，赖文生课题组提出了如图6b所示的辐照损伤高通量计算方法，扫描遍历各种可能的中子注入角度，并考虑注入微区的尺寸因素，自动且并发地生成注入区结构模型并产生分子动力学计算的参数配置，随后进行批量损伤模拟并识别初级辐照缺陷。通过该模型开展高通量计算，然后依据注入微区三维尺寸和注入中子能量及注入角度的球坐标（）进行编码和入库，最终可获取比较完备的初级辐照损伤数据库，如图6c所示。

图6 离位阈能数据库的构建Fig.6 Construction of the displacement threshold energy database: a) the angular dependence of the threshold displacement energy[34; b) the high-throughput calculation flowchart of primary damage; c) coding and storage of the simulation data of primary damage data

2.3 缺陷扩散反应动力学

在多尺度高通量计算模拟框架中，第一性原理计算和分子动力学主要处理的是个体的辐照缺陷在演化过程中的热力学和动力学的物理机制及其数据库建设。缺陷扩散反应动力学将综合利用所有这些基础数据，计算分析它们的耦合作用，最终获取缺陷密度和分布等与辐照损伤及材料性能直接相关的物理量。按照对空间的划分，缺陷扩散反应动力学可分为离散空间模型的动力学蒙特卡罗和连续介质模型的速率理论。

动力学蒙特卡罗对空间做出离散的划分，认为每个辐照缺陷都是一个整体的对象，不对其内部的原子结构细节作任何区分。在这种模拟方法中，通过人为定义的列表对辐照缺陷的种类、大小和取向等性质加以限定，根据第一性原理计算和分子动力学计算设置扩散反应类型以及生成能和迁移能数据库。结合1.2节中的过渡态理论（式5）计算当前所有辐照缺陷对象发生分解、复合、迁移等反应事件的频率，从而促使系统进行缺陷演化，实现多尺度的链接。

由于对空间的离散化处理，动力学蒙特卡罗更适合需要有空间分辨率的研究。比如研究初级辐照损伤的退火现象，或研究位错/晶界等高维晶体缺陷吸收辐照缺陷的强度。动力学蒙特卡罗研究发现，在初级辐照损伤的退火过程中，存在多个回复阶段，在回复末期一般仍存有大量的单空位，而自间隙原子全部复合成自间隙原子团，如图7所示。

图7 不同温度条件下初级辐照损伤退火后的辐照缺陷结构[15]Fig.7 The defect structures of annealed samples after primary damage under different temperatures[15]

动力学蒙特卡罗的研究体系已经可以达到数十或数百纳米量级，但是距离核电材料的器件量级依然较小。速率理论方法比动力学蒙特卡罗更加粗粒化，它将辐照缺陷体系视作一个具有浓度性质的物理场参量，从而高效地分析辐照缺陷的浓度演化规律。

速率理论的基本控制方程可归纳如下:

式中：C是某个感兴趣的辐照缺陷的浓度；G是中子注入引入的该种缺陷的浓度流；其后的三组描述了缺陷反应所导致的浓度流，其中J描述了大小为的缺陷长大为的流量，其他可依此类推。式（10）根据缺陷结合能 E计算了缺陷分解/复合的流量，而式（11）根据缺陷迁移速率计算了缺陷相遇发生反应的频率。以上的特征能量均可从第一性原理计算或分子动力学计算获得。

速率理论计算得到的304不锈钢在经受离子辐照后自间隙性位错环（Frank loops）及空洞（cavities）的密度随深度的变化如图8所示。一般而言，速率理论计算将综合利用第一性原理计算、分子动力学计算以及动力学蒙特卡罗得到的数据库进行耦合作用分析，得到具有统计意义的缺陷密度及分布。因而，这种多尺度计算框架对核电材料服役行为具有直接的指导意义。

图8 304不锈钢离子辐照后缺陷密度随深度变化曲线Fig.8 The relationship between the defect density and the depth of post-irradiation 304 stainless steels: a) the self-interstitial type dislocation loops; b) the voids

3 辐照缺陷演变与材料性能的关联性

尽管不同辐照条件下材料性能退化的规律不完全一致，但材料缺陷结构与材料性能之间具有一定的构效关系，明晰了材料的缺陷就可以指导材料性能的预测。具体而言，根据经典的金属材料硬化模型，核电材料的硬化效果来自于辐照损伤后的微观组织内各种不同种类、结构的缺陷所贡献的硬化因子的综合作用。

Russel-Brown机制中：是施密特系数；是析出相切变模量；是基体切变模量；是位错的柏氏矢量；是位错的内切割半径；是位错的外切割半径；是析出相半径；是析出相平均间距。Orowan机制中：是强化因子；是泰勒系数；是切变模量；是位错柏氏矢量；是基体缺陷数密度；是基体缺陷平均直径。

屈服强度的增加Δ 与材料脆化评价指标韧脆温度转变增量 Δ的关系可整理为，

式中：是关联系数，取决于实验变量和材料的原始性能。

可以看出，核电材料辐照后缺陷结构的种类、尺寸、密度和分布决定了材料的硬化性能，当以上几个缺陷结构因素具有相近数值时，材料应具有等效的辐照损伤性能。各种辐照环境下的缺陷演变规律可以采用多尺度高通量计算模拟进行预测，从而获取等效缺陷结构及其产生条件，其中缺陷结构演变导致的特征能量变化是辐照缺陷演化的关键，也是多尺度高通量计算模拟的基础。因此，在材料基因组思想的指导下，笔者把缺陷结构的特征能量统称为“材料基因组结构能”，提出高通量模拟计算与高通量实验相结合的核电材料服役评价方法：通过高通量计算模拟获得材料缺陷结构数据库，同时通过高通量实验结合数据挖掘技术得到材料“缺陷结构-性能”的构效映射关系。这种方法有望改变核电材料服役评价的传统方法，加快核电材料评价进程。在材料基因组结构能的基础上已经可以实现第一性原理计算、分子动力学、缺陷扩散反应动力学等多尺度模拟方法有效耦合，构成了核电材料中缺陷结构模拟预测的辐照损伤多尺度计算模拟链条。把这些缺陷结构信息导入到构效关系的研究中，并采用相场、有限元等方法开展多物理场耦合模拟，将可以较为完整地研究核电材料的辐照损伤。在后续高通量实验数据验证方面，将基于“结构-性能”构效映射关系建立以材料基因组结构能为参数传递的工程模型，并推进模型在核电材料服役行为的评价和寿命预测中的工程应用，切实提升核电材料服役评价的预测效率及可靠性。