辛之夼，聂宁明，贺新福，王彦棡，吴 石，王 珏

(1.中国科学院 计算机网络信息中心，北京 100190；2.中国科学院大学 计算机科学与技术学院，北京 100049；3.中国原子能科学研究院，北京 102413)

材料辐照损伤一直是材料计算重要的研究领域。材料在反应堆中经受高能中子轰击产生过饱和缺陷，辐照缺陷在温度场、应力场等作用下经过长时间演化，导致材料微结构演化乃至宏观性能退化，是影响反应堆安全性和经济性的关键因素之一。随着物理理论研究的深入和高性能计算技术的发展，研究人员采用实验加先进建模再加计算模拟的手段对材料辐照损伤的复杂机理开展研究。

速率理论[1](RT)是一种用于模拟材料辐照损伤微观结构长时间演化的重要方法，相比较分子动力学(MD)、动力学蒙特卡罗(KMC)[2]这类直接模拟演化轨迹的方法，RT方法基于平均场近似，忽略了缺陷演化的空间信息，在均匀介质中采用反应速率表征缺陷间相互作用的快慢，从而实现对微观结构演化行为的高效模拟。RT方法模拟相较于模拟金属中辐照缺陷随时间演化的常用方法——实体动力学蒙特卡罗方法[3](OKMC)而言，其能模拟约100 dpa下各种类型缺陷的演化过程，远强于OKMC方法，因此在模拟高损伤剂量材料微观结构演化时，能明显体现出优势。而且RT在介观层面上进行模拟，是建立材料宏观性能与微观组织之间关联关系的重要手段。近30年来，国际上采用RT方法成功模拟了反应堆结构材料内辐照诱导各类微观缺陷随损伤剂量的演化行为，如位错环[4-5]、析出物[6]等，并用于预测辐照脆化[7]、辐照肿胀[8]等宏观性能。

RT方法的应用难点主要在于对刚性极强的大规模RT方程组进行数值求解。受限于计算能力，RT方程求解方法大致可归为确定性方法和随机性方法两类，这两类方法的基本思想均是采用近似方法对方程组进行方程数量缩减，以减少计算量。对于确定性方法，一种是按团簇尺寸进行分组求解的分组方法(Group-Method)[9]：用“平均方程”来代替方程组，能有效地减少欲求解方程的数量，降低了计算复杂度，但分组需要技巧。另一种为 Fokker-Planck方法[10]：对于小尺寸团簇用主方程求解获得离散的解，对于大尺寸团簇将其转变为Fokker-Planck方程进行连续化近似求解。该方法对于求解二元体系(间隙-空位)很有效，但对于有3种或以上类型缺陷体系的求解很受限制。对于随机性方法，其中比较经典的是随机团簇动力学方法(SCD)[11]。SCD方法是随机模拟算法(stochastic simulation arithmetic,SSA)[12]在核材料模拟中的一个扩展，利用随机方法求解方程组，减少了欲求解方程的数量，从而降低了求解难度。两种方法的目的均是为了在保证准确度的同时，减少计算量，提高求解效率。

本文引入指数时间差分(ETD)方法对大规模的RT方程组进行并行求解。ETD方法[13]是一类用于求解刚性常微分方程组的高精度数值算法，与传统的半隐式欧拉法相比，ETD方法均能达到对应的最优时间收敛率且表现优于传统的半隐式欧拉法[14]。ETD方法具有良好的稳定性[15]，其对非线性方程的线性部分进行精确求解，对非线性部分进行插值近似，在保证格式精度的同时，可很好地解决由于稳定条件的限制导致时间步长过小的问题。利用ETD方法的数值计算优势对RT主方程直接求解，可得到更符合RT方程的数值解，从而对材料辐照损伤模拟有更好的模拟效果。

本文针对模拟辐照损伤的RT方法和ETD方法进行介绍，分析ETD方法在求解RT主方程中使用的可能性，并介绍ETD方法在空位团簇演化模型和点缺陷可动模型两类RT方程求解中的不同计算实现。

1 RT模型和数值求解算法介绍

辐照导致材料内微观缺陷的形核、长大行为是一个跨越多个时间(从缺陷产生约10-15s到微结构演化约s/h/a)、空间尺度(点缺陷约0.1 nm到微结构约μm)的演化过程，描述微观结构演化的MD、KMC等方法受计算资源的限制，往往无法模拟微结构的长时间演化行为，而RT忽略了体系的随机效应和空间的关联性，因此在模拟缺陷演化的过程中不会受到时空尺度的限制，其模拟结果能与微观表征结果(如TEM、SANS等)进行直接对比。因此RT在辐照导致材料内微观缺陷演化过程的研究中得到了非常广泛的应用。

1.1 RT方法

RT方法类似于化学反应过程的求解，将缺陷之间的相互作用简化为一组微分方程来描述。在辐照条件下，材料中同时存在多种缺陷，包括点缺陷、缺陷团簇以及位错、晶界等缺陷阱。缺陷的产生、聚集、释放以及阱吸收等相互作用均会在缺陷演化RT主方程中予以考虑。RT的一般形式如式(1)所示。

CNDEF(k,t)-∑w(j,l)×CNDEF(j,t)-

Stot×CNDEF(j,t)

(1)

其中：CNDEF(j,t)为t时刻尺寸为j的NDEF类型缺陷的浓度；GNDEF(j)为级联碰撞过程产生尺寸为j的NDEF类型缺陷的速率；w(k,j)为缺陷尺寸由k变为j的转化速率；Stot为材料内固有缺陷对尺寸为j的缺陷的吸收速率；主方程右侧的第2项表示由尺寸k的缺陷团簇生成尺寸j的团簇的反应速率，右侧第3项表示由尺寸j的团簇生成其他尺寸的团簇的生成速率，右侧第4项表示材料内部固有缺陷对尺寸j的团簇的吸收项。

RT方法通过RT主方程(式(1))描述团簇数密度随时间变化的演化过程。方程维数与缺陷数量相关。伴随损伤剂量增加，位错环、空洞尺寸逐渐增大，预测高剂量下位错环等缺陷的演化行为，通常需要百万量级的方程数量。伴随着合金成分的增加，方程数量呈几何型增长。且这一类方程通常是刚性病态的微分方程组，对数值求解提出了较高要求。因此，对该大规模微分方程组的高效精确的数值求解是其主要应用难点。

1.2 ETD算法

ETD算法是由Cox和Matthews提出的一类用于数值求解刚性微分方程的时间积分方法[16-17]。由于隐式算法格式通常具有较高的计算稳定性，刚性微分方程的数值求解通常采用隐式求解方法。但隐式算法格式在计算过程中需对方程进行迭代求解，计算量较大，对于长时间模拟不友好。在ETD方法的求解过程中，首先将刚性方程组分解成线性和非线性部分，对于被积函数中的线性部分直接进行精确计算，对非线性部分采用插值逼近的方式进行求解，这样可得到非线性部分的高精度积分。所以ETD方法可在保证精度的情况下，对刚性微分方程进行显式地稳定求解，同时相较于隐式算法格式大幅降低了计算量。Cox和Matthews给出了ETD-RK格式[16]，应用到具有快速衰变特征的刚性方程的求解上。ETD还被应用于相场方程求解[18-19]，结果表明与传统相场求解方法相比，求解效果与使用半隐式欧拉法的求解效果相比具有更好的收敛精度。

2 ETD方法应用于RT模型求解的典型算例

基于上述优势，采用ETD方法对Ni中空位团簇演化以及奥氏体不锈钢中微结构(位错环以及空洞)演化过程建立RT模型并进行数值求解，并详细介绍整个应用过程。

2.1 空位团簇演化模型

1) 模型建模与算法格式推导

当金属材料内某个原子逃离原有晶体点阵位置，则在该位置出现空缺，并导致空位形成，在一定温度条件下，金属材料内的空位处于热平衡状态[20]。长时间热老化服役条件下，材料内的空位会在热力学驱动下发生扩散，并在迁移过程中不断发生聚集、释放反应。空位间的聚集、俘获反应导致空位团簇不断长大，而空位的释放则引起空位团缩小，伴随服役时间延长，材料内的空位型缺陷不断发生相互作用，从而引起各尺寸空位缺陷浓度随时间发生变化，其演化方程形式如式(2)所示。

(2)

其中：Cvac为单个空位的浓度；Cn为尺寸为n的空位团簇的数密度，n=2,…,N；βn为吸收速率；αn为发射速率，可通过下式计算：

(3)

(4)

由于模拟过程中没有新的空位点缺陷产生，所以要求保持物质总量Qtot守恒：

(5)

初始条件通过下式计算：

(6)

为后续ETD格式的推导，上述模型可表示为如下矩阵形式：

(7)

其中：C(t)为t时刻的缺陷浓度，C(t)=(Cvac,C2,…,Cn)；右边第1项为线性部分，第2项F为非线性部分，表示非线性算子对C(t)的作用。

首先将上式乘以积分因子e-tL，然后在1个时间步(从t=tk到tk+1=tk+Δt)上对方程进行积分，得到如下形式：

C(tk+1)=eΔt·LC(tk)+

(8)

其中：k为迭代的时间步数；Δt为时间步长。采用数值积分求解式中积分项。右端的积分项采用不同数值近似算法可得到不同的ETD格式。

ETD1格式：

C(t+Δt)=C(t)·eΔt·L+

(9)

其中，I为单位对角阵。在空位团簇演化模型求解中，L由下式的对角阵给出：

(10)

F(t)由下式计算：

(11)

其中：Cvac为单尺寸空位原子在t时刻的浓度；Cn为尺寸为n的团簇在t时刻的浓度；n为团簇尺寸；N为团簇的总体尺寸；F(t)矩阵的1～N行分别对应非线性算子对尺寸为1～N的团簇的作用。

ETD4RK格式：

(12)

L由下式计算：

(13)

poly1、poly2、poly3由下式计算：

(14)

(15)

(16)

F(t)由下式计算：

(17)

a(t)、b(t)、c(t)由下式计算：

a(t)=φ1Ct+φ2F(C(t),t)

(18)

(19)

C(t)=φ1a(t)+

(20)

φ1、φ2、φ3由下式计算：

(21)

(22)

(23)

其中，γ为调整系数，用于数值求解时的精度控制。

2) 空位团簇演化模型算例结果分析

基于文献[9]中的物理模型与参数，本文尝试采用上一节中给出的ETD算法对速率方程进行求解验证。模型参数列于表1。

采用ETD算法对空位及其团簇演化的RT主方程进行求解，并将1阶的ETD1和4阶ETD4RK格式计算结果与Group-Method方法的计算结果[9]进行对比，两种方法的对比结果如图1所示，图中示出了金属镍及镍基合金热老化103s条件下的空位团簇尺寸分布。采用1阶、4阶ETD算法的计算结果与Group-Method方法的计算结果基本相近，证明了ETD算法在求解空位团簇演化RT模型方程中的可行性与有效性。

为验证算法稳定性，对基于缺陷数密度的守恒量进行数值验证。在热老化条件下，虽然空位能聚集形成空位团簇，但没有外来缺陷的影响，总的空位数量始终保持守恒，其守恒量的形式可表达为：

表1 金属镍及镍基合金热老化速率方程模型参数Table 1 Model parameters of thermal aging rate equation for nickel metal and nickel-based alloys

图1 金属镍及镍基合金823 K热老化103 s时空位团簇尺寸分布采用不同算法的计算模拟结果Fig.1 Simulated results of vacancy cluster size distribution of nickel metal and nickel-based alloys during thermal aging at 823 K for 103 s using different algorithms

(24)

其中：Cn为缺陷数密度；n为缺陷尺寸；t0为初始时间；tm为第m个时间步，在演化过程中，空位数量值保持稳定，且与初始状态空位数相同。

经计算，初始浓度设置为1.0×10-7的条件下，守恒量为1.054 7×10-7。在106次迭代后，守恒量如图2所示，仍然保持在这个值附近。在迭代计算的前100步左右，守恒量未立即达到1.054 7×10-7，这是由于初期计算缺陷数密度时，缺陷数密度计算结果基本均小于机器精度。随着计算的进行，缺陷数密度计算结果逐渐增大，最后守恒量数值逐渐增加最后趋于稳定。

2.2 点缺陷可动模型

1) 模型建模与算法格式推导

点缺陷可动模型采用Poker于2004年针对中子辐照奥氏体不锈钢位错环演化行为建立的团簇动力学模型。在压水堆服役环境下，堆内构件用奥氏体不锈钢内的微观缺陷以间隙型位错环为主，这也成为其辐照促进腐蚀开裂的1个主要诱因。为模拟预测奥氏体不锈钢内位错环的演化行为，国内外先后建立了描述辐照诱导奥氏体不锈钢内位错环演化行为的团簇动力学模型，中国原子能科学研究院开发的团簇动力学RADIEFF软件[21]，已广泛应用于模拟预测不同辐照条件下奥氏体不锈钢内位错环的演化行为。与上述空位团簇演化模型不同，点缺陷可动模型需考虑间隙型缺陷的演化行为以及间隙、空位型缺陷之间的相互作用，其具体方程形式如式(25)所示。

图2 守恒量验证结果Fig.2 Verification result of conserved quantity

(25)

其中：

(26)

方程中所涉及的缺陷反应类型包括缺陷的产生、缺陷之间的复合、吸收、释放等，点缺陷可动速率方程模型涉及的缺陷反应列于表2。

表2 点缺陷可动速率方程模型涉及的缺陷反应Table 2 Defect reaction involved in point defect moveable rate theory model

对于点缺陷模型方程(式(25))这样复杂的非线性方程组，采用ETD方法进行数值求解，对方程组进行时间上的离散处理以获得显式的时间步进迭代公式。首先对方程(25)进行处理。为表述简便，采用矩阵形式(式(27))进行ETD算法格式的描述。

(27)

将方程组(25)的右端项划分为非线性项和线性项两个部分。线性项矩阵L选取形式如式(28)所示，对于其余的部分组合为非线性项部分Fc。

L矩阵仅选取了方程组(25)右端项线性部分系数的主对角线元素和副对角线元素，而未选取其他线性项系数。是因为经过实验发现，不论是在L矩阵中添加其他线性项的系数还是在L矩阵直接取为方程的Jacobi矩阵，均未对计算结果精度有更明显的改进。综合考虑计算精度和计算效率，本文选取了形式如式(28)的L矩阵进行ETD算法格式的构造。

(28)

对方程(27)的求解等价于求解下面的方程：

(29)

对式(29)进行积分后再进行数值逼近，推导过程与空位团簇模型中1阶算法格式推导类似，得到1阶的ETD格式的迭代式(式(9))。

为提高计算结果的准确度，采用预估校正法对每步时间步所得的结果进行修正，设计2阶的预估校正的ETD算法格式。使用式(9)求得初步的近似值Fc0，被称为预测值。此时预测的精度可能很差，所以再使用1次式(9)，对其进行1次矫正，得到校正值Fc1。

Cn+1=Cn·eΔt·L+(Fc1-Fc0)·

(30)

由于团簇尺寸在百万量级，即方程组数量N=106，直接的数值求解在单进程上计算效率很低甚至无法计算。本文对上述求解格式进行并行求解，将整个方程分解为M个部分，利用M个进程，每个进程负责1个子区域的计算，达到子区域之间并行计算的目的，则有N=n1,n2,n3,…,nM，如图3所示。

并行划分后，模拟区域的边界是原始方程对应的边界，而子区域划分产生的相邻子区域之间的边界为内部边界。考虑边界条件的计算，求解时需分为0号进程边界和其他进程边界。

图3 点缺陷可动模型的并行划分Fig.3 Parallel partitioning of point defect movable model

将划分到每个进程上的部分方程分为控制域和计算域两个部分，控制域组合成原始方程的解，计算域中通过重叠计算的方式通信数据，重叠方式如图4所示。由于并行划分，在划分处会导致数值的丢失。实现过程中采用重叠计算的方式来降低并行划分时数值丢失对结果的影响，通过通信边界值来尽可能地减少数值的丢失。

图4 重叠方式示意Fig.4 Overlapping method illustration

2)点缺陷可动模型算例结果分析

对点缺陷可动模型进行ETD算法求解的算例计算。针对辐照诱导奥氏体不锈钢内位错环演化行为建立团簇动力学模型涉及大量模型参数，包括辐照参数、材料微观参数等，这些参数通常来源于低尺度的模拟计算或由实验数据获得的结果。为对比ETD算法的有效性，本模型所选参数与Poker等[22]所给参数一致(表3)。

采用预估校正ETD方法的计算结果、团簇动力学程序RADIEFF模拟结果以及实验结果对比如图5、6所示。

表3 中子辐照SA-304 SS位错环演化的团簇动力学模型所需参数Table 3 Parameter required for cluster dynamics model of neutron-irradiated SA-304 SS dislocation ring evolution

图5 损伤剂量与直径之间的关系Fig.5 Relationship between irradiation damage dose and diameter

由图5、6可知，采用ETD方法获取的计算结果与经典团簇动力学方法[22]的模拟结果相吻合，且与实验观察数据的变化趋势相符，在高损伤剂量下，模拟结果略高于实验结果，但二者没有量级的差别。利用ETD方法求解的计算数密度结果与304钢的实验值较为接近，数值求解格式解得的半径分布与实验结果非常接近，从而验证了ETD方法在点缺陷可动模型求解中的可行性与有效性。

图6 位错环数密度随辐照损伤剂量的变化 Fig.6 Variation of dislocation ring number density with irradiation damage dose

3 总结与展望

本文对RT方法的原理和建模以及ETD方法进行了介绍。通过对空位团簇演化模型和点缺陷可动模型两类RT模型求解的计算实现，对ETD方法在求解模拟辐照损伤模型中的应用进行了详细介绍，计算结果显示了ETD方法在求解RT方程中的有效性。总之，使用ETD方法求解RT模型的优势，主要体现在，ETD方法能对RT主方程进行高精度的数值求解，减少了额外简化方程所带来的误差，而并行求解的引入使得对于大规模RT方程组的直接数值求解成为可能。可通过直接对主方程进行求解的方式，减少额外简化方程所带来的误差，从而更好地对辐照损伤进行模拟。

但ETD方法在应用实现中仍存在不足。ETD计算过程中时间步长和不同精度算法的选择需随模型方程刚性变化而变化。为更加高效地实现ETD算法在速率方程高精度稳定求解，自适应的变阶变步长ETD算法的实现需在下一步工作中继续开展研究。