CN110459269A - 一种核反应堆材料辐照损伤的多尺度耦合模拟方法 - Google Patents

Abstract

本发明涉及一种核反应堆材料辐照损伤的多尺度耦合模拟方法，该方法将MD、KMC、CD方法结合起来，发挥各自的优势，用MD模拟级联碰撞过程，模拟结果作为KMC的输入，用KMC模拟级联缺陷的退火过程，模拟结果作为CD的输入，用CD模拟辐照缺陷的长时演化过程，则可实现从缺陷产生到微观结构表征的大时空尺度模拟，从而实现服役条件下的材料性能预测。本发明突破了单一尺度的限制，实现核材料辐照损伤的原子尺度缺陷产生到微观结构表征的长时间、大规模模拟，可以深入理解核材料辐照损伤机理，预测辐照诱导材料微结构演化。

Description

一种核反应堆材料辐照损伤的多尺度耦合模拟方法

技术领域

本发明涉及核材料辐照损伤计算机模拟技术，具体涉及一种核反应堆材料辐照损伤的多尺度耦合模拟方法。

背景技术

核反应堆中关键部件材料(如：压力容器、堆内构件、燃料包壳等)的服役性能直接影响着核反应堆的安全性和经济性。结构材料在服役过程中受到高能中子辐照，撞击晶格点阵原子，产生级联碰撞引起初始损伤缺陷，大部分缺陷会在高温下复合湮灭，未被复合的缺陷在服役温度、应力等作用下，经过长时扩散、聚集形成缺陷团簇、析出物、孔洞、气泡、位错环等微观结构，他们与材料内位错、晶界等发生作用，在服役环境(高温、高压、高辐射)下，最终引起材料宏观性能变化，主要表现为辐照脆化、辐照硬化、辐照肿胀、辐照蠕变等辐照效应，严重影响核反应堆的安全运行。由于核材料辐照损伤的实验难度大、周期长且成本高昂，少量实验结合计算机模拟技术是研究材料性能降级的主要方法。

材料辐照损伤的时间尺度跨越飞秒到数年(10^-15s～years)，空间尺度跨越纳米到米量级(nm～m)，是一个从原子碰撞过程到宏观性能表征的多尺度物理过程。目前国际上针对核材料的研究主要采用多尺度模拟(Wirth B D,Odette G R,Marian J,etal.Multiscale modeling of radiation damage in Fe-based alloys in the fusionenvironment[J].2004.)结合高性能计算技术，探究核材料辐照效应的微观机理，其中，多尺度模拟主要集中在原子尺度和微观尺度的研究上。

分子动力学(MD，Molecular Dynamics)是一种有效的原子尺度模拟技术，能够研究级联碰撞过程和缺陷的原始形态、缺陷形成能和迁移能，是一种确定性的研究方法，能够精确追踪体系内粒子的运动过程。但由于算法、计算机内存和计算能力的限制，目前MD模拟在时间尺度上最长可达到秒量级，空间尺度上最大可达到10¹²个原子(～μm³)，其时空尺度不能同时达到。动力学蒙特卡洛(KMC，Kinetic Monte Carlo)方法将分子动力学方法与蒙特卡洛方法相结合，通过构造随机过程或概率模型，可以模拟大规模、长时间的体系演化过程，在材料辐照损伤的微观尺度模拟上得到了广泛应用。相比较MD和KMC这些模拟微观缺陷演化过程更为详细的方法来说，团簇动力学(CD，Cluster Dynamics)方法简化了底层的物理模型，具有计算速度快、计算效率高且无时空尺度限制等优点，其所能模拟的辐照剂量可达上百dpa，模拟结果可与实验直接进行对比。

为保证模拟结果的准确性，MD模拟的时间步长通常为飞秒量级，使得MD所能模拟的级联碰撞(即，初始缺陷的产生过程)非常短(约几十皮秒)，在如此短的时间内尚不能发现合金元素或杂质对辐照缺陷的影响。而KMC的迭代时间步长可取到纳秒量级，所能模拟的时间尺度可达到秒级别。但KMC模拟需要追踪体系中每个缺陷的演化过程，由于计算的复杂性，使得KMC能模拟的辐照剂量一般小于1dpa(dpa，displacement per atom)，很难模拟含有多种缺陷的复杂体系的演化过程。基于平均场近似的CD方法，通过简化底层的物理模型，可以模拟上百dpa下的辐照缺陷演化过程，所能模拟的时间尺度可达到小时或年。但CD的计算效率是以忽略了微观结构演化的波动性和相关性为代价的，而KMC模拟可以通过考虑原子或近原子细节的微观结构来克服CD的这种缺点。但KMC和CD均不能实现缺级联碰撞过程的模拟。

发明内容

本发明的目的在于针对目前核反应堆材料辐照损伤多尺度模拟中单一尺度模拟局限性的问题，提供一种将分子动力学、动力学蒙特卡洛和团簇动力学方法相结合的多尺度耦合模拟方法，通过将三个不同尺度的模拟方法串连耦合起来，借助并行计算技术，突破单一尺度的限制，实现核材料辐照损伤的原子尺度缺陷产生到微观结构表征的长时间、大规模模拟，深入理解核材料辐照损伤机理，精确预测核材料辐照性能。

本发明的技术方案如下：一种核反应堆材料辐照损伤的多尺度耦合模拟方法，包括如下步骤：

(1)初始化系统设置和原子信息，包括系综、温度、压力、第一个被中子撞击离位的点阵原子(PKA)的入射方向及能量、原子的初始坐标、初始速度、邻居原子；

(2)用分子动力学(MD)方法模拟一种PKA导致的级联碰撞过程，同一初始条件下多次重复模拟，对分子动力学模拟得到的原子信息进行初始缺陷统计分析，得到初始缺陷分布；

(3)将所述初始缺陷分布作为动力学蒙特卡洛(KMC)模拟的输入信息，用动力学蒙特卡洛方法模拟缺陷退火过程，对动力学蒙特卡洛模拟得到的间隙和空位分布进行缺陷统计，得到缺陷团簇的尺寸及数密度分布；

(4)对于每一种级联碰撞过程，重复步骤(2)～步骤(3)，得到不同PKA、温度条件下的缺陷团簇尺寸及数密度分布；

(5)将所述缺陷团簇的尺寸及数密度分布信息作为团簇动力学(CD)模拟的输入信息，由团簇动力学模拟缺陷团簇的长时演化过程，采用CD方法，兼顾考虑缺陷的空间信息，得到能够表征的微观结构类型及空间分布，为核材料的性能预测提供信息。

进一步，如上所述的核反应堆材料辐照损伤的多尺度耦合模拟方法，步骤(2)中所述的原子信息包括原子的位置、力、速度；所述的初始缺陷分布包括间隙缺陷和空位缺陷的分布。

进一步，如上所述的核反应堆材料辐照损伤的多尺度耦合模拟方法，步骤(2)中分子动力学模拟级联碰撞过程包括如下步骤：

(2-1)将模拟时间t₁和时间步数k₁初始化为0；

(2-2)判断t₁是否小于等于总的分子动力学模拟时间t_mdTime，k₁是否小于等于总的分子动力学模拟步数N_mdSteps；如果是，则进入步骤(2-3)，如果否，则模拟结束；

(2-3)遍历每个原子及邻居原子，计算作用于每个原子上的力；

(2-4)求解牛顿运动方程，计算每个原子的加速度、速度、位置；

(2-5)将体系前进时间Δt₁，同时时间步数k₁加1，与邻居进行通信，将边界原子信这句话不通顺，模拟时间和时间步数如何调整，文字描述建议与图3所示对应息存储在ghost区中，返回步骤(2-2)。

进一步，如上所述的核反应堆材料辐照损伤的多尺度耦合模拟方法，步骤(3)中动力学蒙特卡洛方法模拟缺陷退火过程包括如下步骤：

(3-1)初始化级联缺陷的类型、数量、及位置信息；

(3-2)计算每个缺陷往每个近邻方向跃迁的概率；

(3-3)计算所有缺陷的跃迁概率总和R，该总跃迁概率为每个进程上的每个缺陷的跃迁概率之和；

(3-4)将模拟时间t₂和时间步数k₂初始化为0；

(3-5)判断t₂是否小于等于总的动力学蒙特卡洛模拟时间t_kmcTime，k₂是否小于等于总的动力学蒙特卡洛模拟步数N_kmcSteps；如果是，则进入步骤(3-6)，如果否，则模拟结束；

(3-6)选择一个随机数r₁(r₁∈[0,1])，计算时间步选择另一个随机数r₂(r₂∈[0，1])，选择跃迁事件μ，其中μ满足R_i为第i个跃迁事件的跃迁概率，R为总跃迁概率；

(3-7)更新与跃迁事件μ相关的缺陷，若被更新的缺陷属于边界缺陷，则与邻居进程进行通信，将更新缺陷存储在ghost区中；

(3-8)更新缺陷的跃迁概率；

(3-9)更新体系的总跃迁概率；

(3-10)体系前进时间Δt₂,同时时间步数k₂加1,返回步骤(3-5)。

进一步，如上所述的核反应堆材料辐照损伤的多尺度耦合模拟方法，步骤(5)中所述缺陷团簇的长时演化过程包括：团簇的长大、复合、被晶界吸收、扩散。

进一步，如上所述的核反应堆材料辐照损伤的多尺度耦合模拟方法，步骤(5)中团簇动力学模拟缺陷团簇的长时演化过程包括如下步骤：

(5-1)将团簇动力学模拟区域均匀网格划分，采用立方体网格，每个网格内的缺陷视为均匀分布，缺陷的空间坐标为网格中心坐标，将每个网格根据其中心坐标分配到相应的进程中，并初始化每个网格6个方向的邻居网格；

(5-2)初始化每个网格中的缺陷列表和反应列表，对于每一种缺陷，包括其类型、数量、所在网格及进程编号，对于每一种可能发生的反应，包括其反应物类型和产物类型，以及对应的反应速率；

(5-3)根据进程数D将整个模拟区域均匀划分为D个子区域，计算每个子区域d的总反应速率A^d，D个子区域的最大总反应速率为A_max＝max{A^d}；

(5-4)将模拟时间t₃和时间步数k₃初始化为0；

(5-5)判断t₃是否小于等于总的速率理论模拟时间t_rtTime，k₂是否小于等于总的速率理论模拟步数N_rtSteps；如果是，则进入步骤(5-6)，如果否，则模拟结束；

(5-6)选择一个随机数r₁(r₁∈[0,1])，计算时间步选择另一个随机数选择反应μ，其中μ满足 为子区域d中的第i个反应的反应概率；

(5-7)更新与反应μ相关的缺陷，若被更新的缺陷所在网格是邻居进程的边界网格，则与邻居进程进行通信，将更新缺陷更新到各自的边界网格的缺陷列表中；

(5-8)更新缺陷的反应速率；

(5-9)更新体系的总反应概率；

(5-10)体系前进时间Δt₃,同时时间步数k₃加1,返回步骤(5-5)。

本发明的有益效果如下：

核材料的辐照损伤是一个跨越原子碰撞到力学性能失效的多尺度过程，传统的单一尺度的模拟技术存在算法和存储能力的限制，不能实现高辐照剂量下的缺陷长时演化模拟。本发明通过将不同尺度的模拟方法串连起来，借助并行计算技术和超级计算机的强大内存，能够实现缺陷产生到宏观性能表征的长时间、大规模模拟。

多尺度模拟的关键在于各尺度之间的耦合和参数的高效传递，各尺度所采用的物理模型的准确性以及参数传递的正确性是决定多尺度模拟结果正确性的关键。本发明通过多次迭代模拟消除不确定因素导致的结果不正确性，通过缺陷统计实现不同模拟对象之间的转换，以确保参数的高效传递。

附图说明

图1是MD、KMC、CD三种模拟方法可模拟的时间和空间尺度示意图。其中MD主要用来模拟缺陷产生过程(即，级联碰撞过程)，KMC主要用来模拟级联缺陷退火过程，CD主要用来模拟缺陷长时演化过程。

图2是本发明的模拟方法流程图。其中，“级联碰撞1～n”表示不同PKA及之后的级联碰撞过程，“缺陷退火1～n”对应“级联碰撞1～n”，表示不同级联碰撞过程产生的级联缺陷的退火模拟。

图3是MD模拟流程图。

图4是KMC模拟流程图

图5是MD和KMC模拟的二维通信示意图。其中实线为进程边界，P0～P8为进程编号，虚线与P4号进程边界所围区域为P4号进程的ghost区域，箭头表示发送/接收原子信息到邻居进程。

图6是KMC模拟区域到CD模拟区域的三维映射图，即初始化CD模拟区域。

图7是CD模拟流程图。

图8是CD模拟中每个体积元的6个方向的三维示意图。

图9是CD模拟的二维通信示意图。其中实线为进程边界，P0～P8为进程编号，P1、P3、P5、P7区域为P4号进程的邻居进程。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

本发明提供的核反应堆材料辐照损伤的多尺度耦合模拟方法，将分子动力学方法MD、动力学蒙特卡洛方法KMC、团簇动力学方法CD结合起来，发挥各自的优势，用MD模拟级联碰撞过程，模拟结果作为KMC的输入，用KMC模拟级联缺陷的退火过程，模拟结果作为CD的输入，用CD模拟辐照缺陷的长时演化过程，则可实现从缺陷产生到微观结构表征的大时空尺度模拟，从而实现服役条件下的材料性能预测。

具体的方法实现包括了如下步骤：

(1)初始化系统设置和原子信息，包括系综、温度、压力、第一个被中子撞击离位的点阵原子(PKA，Primary Knock-on Atom)的入射方向及能量、原子的初始坐标、初始速度、邻居原子。其中PKA由中子与点阵原子间的二体碰撞近似(Binary CollisionApproximation，BCA)计算得到，中子主要来源于裂变反应释放的快中子(能量在0.1～10MeV，平均能量为2MeV)。由于PKA是第一个被撞离位的原子，所以PKA与入射中子是1对1的关系。

(2)用分子动力学(MD)方法模拟某种PKA导致的级联碰撞过程，同一初始条件下重复模拟15～20次，对分子动力学模拟得到的原子信息(原子的位置、力、速度等)进行初始缺陷统计分析，得到初始缺陷分布(包括间隙缺陷和空位缺陷的分布)。MD模拟级联碰撞过程见附图3，主要步骤为：

A.将模拟时间t₁和时间步数k₁初始化为0；

B.判断t₁是否小于等于总的MD模拟时间t_mdTime，或k₁是否小于等于总的MD模拟步数N_mdSteps；

是：转C；

否：模拟结束；

C.遍历每个原子及邻居原子，计算作用于每个原子上的力；

D.求解牛顿运动方程，计算每个原子的加速度、速度、位置；

E.体系前进时间Δt₁，时间步数k₁加1，与邻居进行通信，将边界原子信息存储在ghost区中，见附图5，其中ghost区即为附图5中P4的边界与虚线所围区域。转B。

(3)将所述初始缺陷分布作为动力学蒙特卡洛(KMC)模拟的输入信息，用动力学蒙特卡洛方法模拟缺陷退火过程，对动力学蒙特卡洛模拟得到的间隙和空位分布进行缺陷统计，得到缺陷团簇的尺寸及数密度分布。KMC模拟缺陷退火过程见附图4，主要步骤为：

a.初始化级联缺陷的类型、数量、及位置信息；

b.计算每个缺陷往每个近邻方向跃迁的概率；

c.计算所有缺陷的跃迁概率总和R，该总跃迁概率为每个进程上的每个缺陷的跃迁概率之和；

d.将模拟时间t₂和时间步数k₂初始化为0；

e.判断t₂是否小于等于总的KMC模拟时间t_kmcTime，或k₂是否小于等于总的KMC模拟步数N_kmcSteps；

是：转f；

否：模拟结束；

f.选择一个随机数r₁(r₁∈[0,1])，计算时间步选择另一个随机数r₂(r₂∈[0,1])，选择跃迁事件μ，其中μ满足R_i为第i个跃迁事件的跃迁概率，R为总跃迁概率；

g.更新与跃迁事件μ相关的缺陷，若被更新的缺陷在属于边界缺陷，则与邻居进程进行通信，将更新缺陷存储在ghost区中，见附图5；

h.更新缺陷的跃迁概率；

i.更新体系的总跃迁概率；

j.体系前进时间Δt₂，同时时间步数k₂加1。转e。

(4)若有多个级联碰撞及随后的退火过程，则重复步骤(2)-(3)，得到多种不同缺陷分布的KMC模拟区域，则将所有这些KMC模拟区域初始化为RT的模拟区域；若只有一个级联碰撞及随后的退火过程，得到一种缺陷分布的KMC模拟区域，则将该KMC模拟区域复制多份，得到CD的模拟区域，见附图6。

(5)将所述缺陷团簇的尺寸及数密度分布信息作为团簇动力学(CD)模拟的输入信息，由团簇动力学模拟缺陷团簇的长时演化过程(例如：团簇的长大、复合、被晶界吸收、扩散等)，采用团簇动力学方法，兼顾考虑缺陷的空间信息，得到能够表征的微观结构类型及空间分布，为核材料的性能预测提供信息。CD模拟缺陷的长时演化过程见附图7，主要步骤为：

①将CD模拟区域均匀网格划分，采用立方体网格，每个网格内的缺陷视为均匀分布，缺陷的空间坐标为网格中心坐标，将每个网格根据其中心坐标分配到相应的进程中，并初始化每个网格6个方向的邻居网格，见附图8；

②初始化每个网格中的缺陷列表和反应列表，对于每一种缺陷，包括其类型、数量、所在网格及进程编号，对于每一种可能发生的但应，包括其反应物类型和产物类型，以及对应的反应速率；

③根据进程数D将整个模拟区域均匀划分为D个子区域，计算每个子区域d的总反应速率A^d，D个子区域的最大总反应速率为A_max＝max{A^d}；

④将模拟时间t₃和时间步数k₃初始化为0；

⑤判断t₃是否小于等于总的CD模拟时间t_rtTime，或k₂是否小于等于总的CD模拟步数N_rtSteps；

是：转⑥；

否：模拟结束；

⑥选择一个随机数r₁(r₁∈[0，1])，计算时间步选择另一个随机数选择反应μ，其中μ满足 为子区域d中的第i个反应的反应概率；

⑦更新与反应μ相关的缺陷，若被更新的缺陷所在网格是邻居进程的边界网格，则与邻居进程进行通信，将更新缺陷更新到各自的边界网格的缺陷列表中，见附图9；

⑧更新缺陷的反应速率；

⑨更新体系的总反应概率；

⑩体系前进时间Δt₃，同时时间步数k₃加1。转⑤。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若对本发明的这些修改和变型属于本发明权利要求及其同等技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (6)

1.一种核反应堆材料辐照损伤的多尺度耦合模拟方法，包括如下步骤：

(1)初始化系统设置和原子信息，包括系综、温度、压力、第一个被中子撞击离位的点阵原子PKA的入射方向及能量、原子的初始坐标、初始速度、邻居原子；

(2)用分子动力学方法模拟一种PKA导致的级联碰撞过程，同一初始条件下多次重复模拟，对分子动力学模拟得到的原子信息进行初始缺陷统计分析，得到初始缺陷分布；

(3)将所述初始缺陷分布作为动力学蒙特卡洛模拟的输入信息，用原子动力学蒙特卡洛方法模拟缺陷退火过程，对动力学蒙特卡洛模拟得到的间隙和空位分布进行缺陷统计，得到缺陷团簇的尺寸及数密度分布；

(4)对于每一种级联碰撞过程，重复步骤(2)～步骤(3)，得到不同PKA的缺陷团簇尺寸及数密度分布；

(5)将所述缺陷团簇的尺寸及数密度分布信息作为团簇动力学模拟的输入信息，由团簇动力学模拟缺陷团簇的长时演化过程，采用团簇动力学方法，兼顾考虑缺陷的空间信息，得到能够表征的微观结构类型及空间分布，为核材料的性能预测提供信息。

2.如权利要求1所述的核反应堆材料辐照损伤的多尺度耦合模拟方法，其特征在于：步骤(2)中所述的原子信息包括原子的位置、力、速度；所述的初始缺陷分布包括间隙缺陷和空位缺陷的分布。

3.如权利要求2所述的核反应堆材料辐照损伤的多尺度耦合模拟方法，其特征在于：步骤(2)中分子动力学模拟级联碰撞过程包括如下步骤：

(2-1)将模拟时间t₁和时间步数k₁初始化为0；

(2-2)判断t₁是否小于等于总的分子动力学模拟时间t_mdTime，k₁是否小于等于总的分子动力学模拟步数N_mdSteps；如果是，则进入步骤(2-3)，如果否，则模拟结束；

(2-3)遍历每个原子及邻居原子，计算作用于每个原子上的力；

(2-4)求解牛顿运动方程，计算每个原子的加速度、速度、位置；

(2-5)将体系前进时间Δt₁，同时时间步数k₁加1，与邻居进行通信，将边界原子信这句话不通顺，模拟时间和时间步数如何调整，文字描述建议与图3所示对应息存储在ghost区中，返回步骤(2-2)。

4.如权利要求1所述的核反应堆材料辐照损伤的多尺度耦合模拟方法，其特征在于：步骤(3)中动力学蒙特卡洛方法模拟缺陷退火过程包括如下步骤：

(3-1)初始化级联缺陷的类型、数量、及位置信息；

(3-2)计算每个缺陷往每个近邻方向跃迁的概率；

(3-3)计算所有缺陷的跃迁概率总和R，该总跃迁概率为每个进程上的每个缺陷的跃迁概率之和；

(3-4)将模拟时间t₂和时间步数k₂初始化为0；

(3-5)判断t₂是否小于等于总的动力学蒙特卡洛模拟时间t_kmcTime，k₂是否小于等于总的动力学蒙特卡洛模拟步数N_kmcSteps；如果是，则进入步骤(3-6)，如果否，则模拟结束；

(3-6)选择一个随机数r₁(r₁∈[0,1])，计算时间步选择另一个随机数r₂(r₂∈[0,1])，选择跃迁事件μ，其中μ满足R_i为第i个跃迁事件的跃迁概率，R为总跃迁概率；

(3-7)更新与跃迁事件μ相关的缺陷，若被更新的缺陷属于边界缺陷，则与邻居进程进行通信，将更新缺陷存储在ghost区中；

(3-8)更新缺陷的跃迁概率；

(3-9)更新体系的总跃迁概率；

(3-10)体系前进时间Δt₂,同时时间步数k₂加1,返回步骤(3-5)。

5.如权利要求1所述的核反应堆材料辐照损伤的多尺度耦合模拟方法，其特征在于：步骤(5)中所述缺陷团簇的长时演化过程包括：团簇的长大、复合、被晶界吸收、扩散。

6.如权利要求5所述的核反应堆材料辐照损伤的多尺度耦合模拟方法，其特征在于：步骤(5)中团簇动力学模拟缺陷团簇的长时演化过程包括如下步骤：

(5-1)将团簇动力学模拟区域均匀网格划分，采用立方体网格，每个网格内的缺陷视为均匀分布，缺陷的空间坐标为网格中心坐标，将每个网格根据其中心坐标分配到相应的进程中，并初始化每个网格6个方向的邻居网格；

(5-2)初始化每个网格中的缺陷列表和反应列表，对于每一种缺陷，包括其类型、数量、所在网格及进程编号，对于每一种可能发生的反应，包括其反应物类型和产物类型，以及对应的反应速率；

(5-3)根据进程数D将整个模拟区域均匀划分为D个子区域，计算每个子区域d的总反应速率A^d，D个子区域的最大总反应速率为A_max＝max{A^d}；

(5-4)将模拟时间t₃和时间步数k₃初始化为0；

(5-5)判断t₃是否小于等于总的速率理论模拟时间t_rtTime，k₂是否小于等于总的速率理论模拟步数N_rtSteps；如果是，则进入步骤(5-6)，如果否，则模拟结束；

(5-6)选择一个随机数r₁(r₁∈[0,1])，计算时间步选择另一个随机数选择反应μ，其中μ满足 为子区域d中的第i个反应的反应概率；

(5-7)更新与反应μ相关的缺陷，若被更新的缺陷所在网格是邻居进程的边界网格，则与邻居进程进行通信，将更新缺陷更新到各自的边界网格的缺陷列表中；

(5-8)更新缺陷的反应速率；

(5-9)更新体系的总反应概率；

(5-10)体系前进时间Δt₃，同时时间步数k₃加1，返回步骤(5-5)。