铁在液态铅铋中团聚和扩散的分子动力学模拟

Molecular dynamics simulation of the aggregation and diffusion of iron in liquid lead-bismuth

摘要

液态铅铋被视为第4代铅冷快堆的堆芯冷却剂，其与结构钢材的相容性问题是目前研究的重点之一。采用分子动力学方法，对铁原子在液态铅铋中的团聚过程和铁团簇的扩散行为进行模拟，分析体系温度、铁团簇尺寸和质量浓度对铁团簇之间聚合过程，以及铁团簇扩散行为的影响规律和微观机理。结果表明，铁原子在液态铅铋中会倾向于形成团簇，并且高温能够加快团聚过程。研究还发现在液态铅铋中铁团簇之间的聚合需要克服能垒，在473 K下半径为10 Å的2个铁团簇聚合需要克服的能垒约为2.14 eV，并且体系温度的升高和铁团簇尺寸的增大有利于铁团簇之间发生聚合。另外，通过计算可知，铁团簇的扩散系数随着其质量浓度的降低而增大，随着其尺寸的增大而减小，随着体系温度的升高而增大。铁团簇尺寸对扩散激活能的影响并不显著，半径为10、12.5和15 Å的铁团簇的扩散激活能分别为0.128、0.140和0.134 eV。

Abstract

Liquid lead-bismuth is considered as the core coolant of the fourth-generation lead-cooled fast reactors， and its compatibility with structural steels is one of the focuses of current research. In this paper， the aggregation process of iron atoms in liquid lead-bismuth and the diffusion behavior of iron clusters were simulated by using the molecular dynamics method. The influence laws and microscopic mechanisms of the temperature of the system， the size and mass concentration of iron clusters on the coalescence process between iron clusters and the diffusion behavior of iron clusters were analyzed. The results show that iron atoms tend to form clusters in liquid lead-bismuth and that high temperatures can accelerate the aggregation process. In addition， it is found that the coalescence between iron clusters in liquid lead-bismuth needs to overcome energy barriers， the energy barrier to be overcome for the coalescence of two iron clusters of radii 10 Å at 473 K is about 2.14 eV， and the increase of system temperature and the increase of iron clusters size are favorable for the occurrence of coalescence between iron clusters. Furthermore， the diffusion coefficients of iron clusters were calculated， and found to increase with the decrease of their mass concentration， decrease with the increase of their size， and increase with the increase of the temperature of the system. The effect of iron clusters size on the diffusion activation energies was not significant， and the diffusion activation energies of iron clusters with radii of 10， 12.5， and 15 Å were 0.128， 0.140， and 0.134 eV， respectively.

Graphical abstract

关键词

分子动力学 / 铁团簇 / 液态铅铋 / 扩散系数

Key words

molecular dynamics / iron clusters / liquid lead-bismuth / diffusion coefficient

铅铋共晶合金因其低熔点、高沸点、高热导率等热力学特性，被视为第4代铅冷快堆的堆芯冷却剂^［1-2］。然而，在高温下液态铅铋与结构钢材的相容性问题是发展铅冷快堆过程中亟待解决的难题之一^［3］。腐蚀会造成结构钢材的服役性能下降，腐蚀产物的沉积更会引起严重的热工水力问题^［4］。

一般来说，结构钢材在液态铅铋中的腐蚀可以分为溶解腐蚀和氧化腐蚀，即结构钢材中铁、铬和镍等合金元素在液态铅铋中的溶解，以及结构钢材表面发生的氧化^［5］。当液态铅铋中溶解氧浓度较低时，结构钢材主要发生溶解，随着溶解氧浓度的升高，结构钢材与氧发生反应并在表面形成氧化层，通过控制溶解氧浓度，可以有效控制氧化层厚度，从而在保证液态铅铋传热效率的前提下，减缓对结构钢材的进一步腐蚀^［6］。

科研人员对液态铅铋腐蚀不同组分的结构钢材开展了大量的实验研究。例如，Barbier和Rusanov^［7］对不同类型结构钢材置于流动液态铅铋中的腐蚀行为进行了对比分析，液态铅铋中溶解氧质量浓度维持在10^-6 %量级，发现实验所用的4种类型结构钢材均形成了保护性氧化层，却未表现出溶解腐蚀的特征。其中，奥氏体钢的氧化层厚度最薄（小于1 μm），表现出最佳的抗氧化腐蚀性，对于另外3种类型的马氏体钢，氧化层厚度介于10～20 μm，且抗氧化腐蚀的能力由强到弱依次为：EP823，T91，Optifer Ⅳ马氏体钢。Hosemann等^［8］对D9奥氏体钢在流动液态铅铋中的氧化行为进行了分析，发现结构钢材表面形成的氧化层并非简单的双层结构，而是至少包含3层的复杂结构。Zhu等^［9］将T91铁素体/马氏体钢置于不同溶解氧浓度的静止液态铅铋中，对形成氧化层的结构、形貌和成分进行详细分析，发现在饱和溶解氧浓度条件下，氧化层由铁铅复合氧化物层、外氧化层（主要成分为Fe₃O₄）、内氧化层（主要成分为铁铬尖晶石）和内氧化区（FeCr₂O₄与结构钢材基体的混合物）组成，随着溶解氧浓度的降低，氧化层厚度逐渐变薄。

液态铅铋腐蚀结构钢材的实验通常长达数千至上万小时，并且难以捕捉腐蚀过程的微观动态变化，计算机技术的发展为从原子尺度揭示腐蚀的本质特征提供了可能。例如，部分科研人员采用基于密度泛函理论的第一性原理方法，通过计算铅铋原子在铁表面的吸附能和表面铁原子的逃逸能等物理量，发现相较于铅原子，铋原子更容易吸附在铁表面，并且铅铋原子的吸附能够降低铁原子从表面溶解所需的能量，显著增强了表面铁原子的溶解倾向^［10-11］。同时，还发现结构钢材中其他元素的存在，比如硅等，对铅铋原子在铁表面的吸附具有抑制作用，从而减缓了表面铁原子的溶解趋势^［12］。虽然第一性原理的计算结果精度很高，但是计算复杂度也很大，计算规模通常限制在数百原子以内，时间尺度也局限于ps量级。分子动力学（molecular dynamics， MD）模拟能够在保证一定精度的前提下，实现更大规模（数万至数百万原子）和更长时间（ns至ms量级）的模拟。Zhou等^［13］对铁与液态铅铋固液界面处的溶解渗透现象进行了MD研究，通过统计不同温度下铋原子的渗透数量和铁原子的溶解数量，发现当温度由823 K升高至873 K时，铋原子渗透数量和铁原子溶解数量会显著增加，溶解渗透行为明显加剧。Chen等^［14］研究了液态铅铋对含有空位缺陷的纯铁的腐蚀行为，通过计算铅铋原子替代铁原子以及与空位结合所需的能量，发现空位缺陷的存在能够增强腐蚀作用，铅铋原子更倾向于通过空位缺陷渗透到铁中。

截至目前，液态铅铋腐蚀结构钢材的微观尺度模拟研究主要关注固液界面的相容性问题，对溶解腐蚀产物在液态铅铋中运动过程的研究尚不充分。鉴于此，本文首先对铁原子在液态铅铋中的团聚过程进行MD模拟，接着通过在模拟区域中引入铁团簇，依次分析铁团簇间的相互作用，以及铁团簇的质量浓度和尺寸对其扩散系数的影响规律和微观机理。

1 模型与方法

为正确描述原子间的相互作用，本文的MD模拟采用由Gao等^［15］基于嵌入原子方法^［16］发展的铁铅铋三元合金势函数，其中，铅铋的参数取自Belashchenko^［17］的结果，铁的参数取自Zhou等^［18］的结果，通过电子密度缩放因子s=0.17^［15］使2组参数兼容协调，截断半径的取值为9.05 Å（1 Å=0.1 nm，下同）^［15］。该势函数的准确性通过计算铁铅之间和铁铋之间的径向分布函数，液态铅铋的黏度，以及铁在液态铅铋中的扩散系数等物理量得到了较好的验证^［19］。本文的模拟采用1 fs的时间步长和三维周期性边界条件，Nosé-Hoover方法^［20-21］控制体系温度，Parrinello-Rahman方法^［22］控制体系压力。

在铁原子团聚过程的模拟中，模拟区域的尺寸设置为20a\times20a\times20a，其中，a=4.91 Å为FCC铅的晶格常数，铅铋原子的质量分数分别为44.5%和55.5%，并将其中的650个铅铋原子随机替换为铁原子。这里，引入的铁原子浓度高于实际情况下铁在液态铅铋中的饱和溶解度，目的是加速铁原子的团聚过程，以便在有限的模拟时间内统计铁团簇的形成规律。体系温度范围介于473～1 500 K，覆盖真实工况下的工作温度区间。考虑到原子运动具有一定的随机性，在每种工况下均进行5次独立的模拟，并对统计结果进行平均处理。

在计算铁团簇扩散系数的模拟中，模拟区域的尺寸根据铁团簇的质量浓度进行设置，并将一个BCC铁团簇置于区域中心以模拟铁原子团聚形成的团簇。为了避免界面处铁原子与铅铋原子之间的强排斥作用，在铁团簇和其邻近的铅铋原子之间设置了0.1 nm的间隙。先在等温等压系综下弛豫0.5 ns使体系达到平衡状态。由于模拟区域中只含有一个铁团簇，其运动具有一定的随机性，然而，扩散系数是大量粒子统计平均的结果，因此，每种工况均在正则系综下进行60～100次独立的单个铁团簇模拟，并将多次统计结果进行平均处理，以模拟含有大量铁团簇的体系，从而计算扩散系数，具体公式如下

D=\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\mathrm{M}\mathrm{S}\mathrm{D}}{\mathrm{6}t}.（1）

其中，D为粒子的扩散系数，t为模拟时间，\mathrm{M}\mathrm{S}\mathrm{D}为均方位移，其定义如下

\mathrm{M}\mathrm{S}\mathrm{D}=〈|\mathit{r}(t)-\mathit{r}(\mathrm{0}){|}^{\mathrm{2}}〉.（2）

这里，|\mathit{r}(t)-\mathit{r}(\mathrm{0})|为单个粒子在时间间隔t内的位移，〈〉表示对组内所有粒子的位移平方取均值。此外，液态铅铋的黏度采用反扰动非平衡分子动力学方法^［23-24］，利用动量通量与速度梯度的线性关系进行计算，具体做法为：将模拟区域沿某一方向（如z方向）划分为若干层，周期性交换中间层和边缘层粒子沿x方向的动量，形成速度梯度从而计算黏度。

本文所有模拟均使用LAMMPS^［25］完成，并使用OVITO^［26］进行可视化后处理。

2 结果与讨论

2.1 铁原子的团聚

在液态铅铋中随机分散着一定数量的铁原子，观察在一定温度下铁原子的团聚过程。图1展示了在473 K下铁原子的运动演化过程，为清楚地显示铁原子的分布，对铅铋原子的演示尺寸进行了缩小处理。从图中可以看出，在模拟的初始阶段（见图1（a）、1（b）），距离较近的铁原子相互靠近并团聚成小团簇，此后，小团簇之间相互聚集，同时吸收附近孤立的铁原子，团簇尺寸逐渐增大（见图1（c））。当模拟进行至100 ns时（见图1（d）），可以看到模拟区域中几乎没有孤立的铁原子，这是因为在473 K下铁在液态铅铋中的溶解度很小^［27］，几乎所有铁原子都已形成团簇。综上所述，团聚过程可以概括为先由分散的、离得近的铁原子团聚成小团簇，再由小团簇聚集成大团簇。其他温度下的团聚过程与此类似，不再赘述。

图2所示为铁团簇数量和孤立的铁原子数量随时间的变化曲线。这里，判断一个铁原子是否属于一个团簇的标准是其与团簇中任意一个铁原子的距离是否在临界值之内，该临界值取在铁晶体径向分布函数第1谷的位置，位于第1近邻和第2近邻原子之间^［28］，本文取为0.27 nm。从图中可以看出，在所有温度下，铁团簇数量均快速增加并达到峰值，孤立的铁原子数量则急剧减少，对应上述过程的第一阶段，即距离较近的铁原子相互靠近并团聚成小团簇。而后，铁团簇数量缓慢减少，孤立的铁原子数量也缓慢减少并趋于稳定，对应上述过程的第二阶段，即小团簇之间相互聚集，同时吸收少量附近孤立的铁原子，团簇尺寸逐渐增大。

铁原子的团聚过程受多种因素的影响，考虑到温度越高，铁原子动能越大，可以更容易地克服周围铅铋原子的束缚作用，从而发生团聚。因此，本文首先研究温度对团聚过程的影响，模拟了不同温度下铁原子的团聚过程。一方面，不同温度下形成的铁团簇构型有所差异，图3所示为在873 K下铁原子的运动演化过程，通过对比图1（d）和图3（d），发现在473 K下铁团簇大多呈细长的棒状构型，而在873 K下铁团簇呈近球状或近椭球状构型，这是由于随着温度的升高，原子热运动加剧，铁团簇融合得更加充分，从而导致形成的铁团簇构型有所差异。另一方面，不同温度下铁团簇数量和孤立的铁原子数量随时间的变化曲线也有区别，如图2所示。在473 K下铁团簇数量在0.5 ns达到峰值，随着温度的升高，达到峰值所需的时间逐渐减少，比如，在873 K和1 500 K下铁团簇数量的峰值分别出现在0.1 ns和0.05 ns（见图2（a）），并且温度越高，在模拟的初始阶段孤立的铁原子数量减少越快（见图2（b）），这说明高温能够加速铁原子团聚成小团簇的过程。此外，通过对比不同温度下铁团簇数量达到峰值后曲线的下降阶段，可以看出温度越高，铁团簇数量减少越快（见图2（a）），这说明高温能够加速小团簇聚集成大团簇的过程。因此，温度升高，会同时加快铁原子团聚成小团簇和小团簇聚集成大团簇2个阶段。

2.2 铁团簇的聚合

为分析铁团簇间的相互作用，首先将一个铁团簇置于模拟区域中，并保持其质心位置固定不变，接着将另一个铁团簇从距离较远处逐渐靠近直至两者接触并聚合。以半径为10 Å的铁团簇为例，图4所示为在473 K下铁原子平均势能随2个铁团簇中心距离的变化曲线。其中，图4（a）为铁团簇在真空中时铁原子平均势能的变化情况。从图中可以看出，当2个铁团簇距离较远时，势能曲线保持水平，说明2个铁团簇之间不存在相互作用力。当2个铁团簇中心距离为21.2 Å时恰好发生聚合，之后，随着中心距离的减小，势能曲线呈下降趋势，说明2个铁团簇之间表现为吸引作用，并且曲线斜率逐渐减小，吸引作用逐渐减弱。当中心距离为18.6 Å时，势能达到局部最小值，2个铁团簇达到此时融合程度下的稳态构型。然而，当接着减小中心距离时，势能曲线转而呈上升趋势，说明2个铁团簇之间表现为排斥作用。随着中心距离的进一步减小，势能发生突降，聚合后的铁团簇内部原子重新排列，2个铁团簇加深融合，达到融合程度加深后的稳态构型。

图4（b）为铁团簇在液态铅铋中时铁原子平均势能的变化情况，与在真空中时的变化情况略有差异。从图中可以看出，当2个铁团簇距离较远时，势能曲线近似水平，说明彼此之间不存在相互作用力，但受到液态铅铋的影响曲线会出现些许波动。值得注意的是，2个铁团簇发生相互作用的临界距离有所增加，当中心距离为23.0 Å时，势能曲线开始呈上升趋势，说明在液态铅铋中2个铁团簇相互靠近时先表现为排斥作用，直至20.8 Å时势能达到局部最大值，受到周围铅铋原子的阻碍作用，铁团簇之间的聚合需要克服约2.14 eV的能垒。之后，与真空中类似，随着中心距离的减小，2个铁团簇发生聚合，势能曲线呈下降趋势，说明2个铁团簇之间表现为吸引作用，当中心距离为18.6 Å时，势能达到局部最小值，2个铁团簇达到此时融合程度下的稳态构型。

通过对比图4（a）和4（b）不难看出，真空中铁原子平均势能高于液态铅铋中铁原子平均势能，这是由于铁团簇在真空中和液态铅铋中的表面能并不相同，表面能E_{\mathrm{s}}可以采用如下公式进行计算

E_{\mathrm{s}}=\frac{E_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}-n{E}_{\mathrm{c}}}{\mathrm{4}\mathrm{\pi }{R}^{\mathrm{2}}}.（3）

其中，E_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}为铁团簇的总势能，E_{\mathrm{c}}为无限大块体中单个铁原子的平均势能，n为铁团簇含有的原子数，R为铁团簇的等效半径。由式（3）可得，在473 K下半径为10 Å的铁团簇置于真空中和液态铅铋中时的表面能分别为1.56和0.84 J/m²。显然，铁团簇在真空中的表面能更高，从而使得铁团簇中的铁原子平均势能更高。

由于液态铅铋的存在，铁团簇间的相互作用将发生一定变化，铁团簇之间的聚合需要克服能垒，这是由于在液态铅铋中，铁团簇的表面会吸附铅铋原子，形成一定厚度的吸附层，处于吸附层的铅铋原子以类固态的形式存在并伴随铁团簇运动。实际上，纳米团簇表面吸附液体原子的现象在很多研究中都有发现^［28-29］。在液态铅铋中2个铁团簇相互靠近时，首先发生相互作用的是表面吸附的铅铋原子，因此临界距离有所增加，由真空中的21.2 Å增加至23.0 Å。然而，由于2个铁团簇之间充满液态铅铋，表面吸附的铅铋原子之间的吸引作用受到屏蔽，当这些铅铋原子之间的距离小于平衡距离时，将发生排斥作用，体现在图4（b）中，即铁原子平均势能升高。随着距离的进一步减小，处于吸附层的铅铋原子不断发生重排，并向2个铁团簇接触点周围扩散，直至2个铁团簇将要发生聚合，势能达到局部最大值。

由以上分析可知，在液态铅铋中铁团簇之间的聚合需要首先克服一定能垒。为进一步验证该结论，模拟了不同尺寸的铁团簇在不同温度下的随机运动过程，体系温度范围介于473～1 500 K，铁团簇的半径分别为6、8、10和12 Å，观察2个铁团簇在运动过程中是否发生聚合。

图5展示了不同温度下半径为12 Å的铁团簇在液态铅铋中的运动过程。从图中可以看出，当温度较低时，2个铁团簇相互靠近后又分开，并未发生聚合，如图5（a）、5（b）所示，在673 K和873 K下，2个铁团簇分别在10.02 ns和0.92 ns运动至最近，而后又分离，这是因为铁团簇的动能不足以克服聚合所需的能垒。当温度升高时，2个铁团簇发生聚合，且温度越高，铁团簇融合得越充分，在1 073 K下（见图5（c）），2个铁团簇在18.14 ns接触并开始聚合，聚合后的铁团簇保持哑铃状构型，在1 273 K下（见图5（d）），2个铁团簇在0.26 ns接触并开始聚合，且它们融合得更加充分，在后续的模拟时间内保持稳定的棒状构型。

除体系温度会影响铁团簇的聚合过程，铁团簇的尺寸也有一定的影响。图6展示了在973 K下不同尺寸的铁团簇在液态铅铋中的运动过程。从图中可以看出，当尺寸较小时，2个铁团簇不易发生聚合，如图6（a）和6（b）所示，半径为6 Å和8 Å的铁团簇分别在8.90 ns和16.06 ns相互靠近，而后又分离。当尺寸增大时，2个铁团簇则更容易发生聚合，如图6（c）和6（d）所示。通过分析聚合过程中铁团簇的构型，发现对于尺寸较大的铁团簇，彼此之间的接触面积更大，会有更多的铁原子相互作用，从而提高了聚合的可能性。值得注意的是，由于尺寸较大的铁团簇内部原子排列更接近理想晶格，在结构上的稳定性使其对外界扰动的抵抗能力更强，聚合后仅在2个铁团簇接触处融合充分，如图6（d）所示，半径为12 Å的铁团簇在1.44 ns接触并开始聚合，在后续的模拟时间内呈哑铃状构型。相比之下，如图6（c）所示，半径为10 Å尺寸较小的铁团簇在1.56 ns接触并开始聚合，当模拟进行至6.84 ns时整体上融合得更加充分，呈棒状构型。以上现象说明，体系温度的升高和铁团簇尺寸的增大有利于2个铁团簇发生聚合。

2.3 铁团簇的扩散

铁团簇在液态铅铋中的扩散行为对腐蚀产物的沉积过程有着重要影响，本节通过计算铁团簇的扩散系数，并分析质量浓度和尺寸的影响，以帮助更好地理解铁团簇在液态铅铋中的扩散行为。

2.3.1 质量浓度的影响

为分析铁团簇的质量浓度对其扩散系数的影响规律，将一个半径为15 Å的球形铁团簇置于液态铅铋中，通过改变模拟区域的尺寸来控制铁团簇的质量浓度。当模拟区域越大时，铅铋原子数量越多，相应地，铁团簇质量浓度越低，图7所示为铁团簇的扩散系数随其质量浓度的变化曲线。从图中可以看出，随着铁团簇质量浓度的降低，其扩散系数迅速增大并趋于稳定。实际上，质量浓度反映的是铁团簇之间的相互约束作用。当质量浓度较高时，相邻铁团簇之间存在相互约束，因此运动扩散能力受到限制；随着质量浓度的降低，相邻铁团簇的距离逐渐增大，彼此之间约束作用减弱，从而使运动扩散能力得到提升；当质量浓度降低到一定程度时，铁团簇之间相互作用消失，扩散系数不再发生变化（见图7红色虚线）。

铁团簇扩散系数的变化规律也可以从液态铅铋黏度的变化规律进行解释。已有的研究结果表明，在液体中加入异质纳米团簇可以改变其黏度^［30-31］。为此，计算了含有不同质量浓度铁团簇的液态铅铋的黏度，如图8所示。从图中可以看出，铁团簇的存在会导致液态铅铋的黏度增大，随着铁团簇质量浓度的降低，这种影响逐渐减弱，液态铅铋的黏度逐渐减小并最终趋于纯液态铅铋的黏度（见图8红色虚线）。以上说明液态铅铋黏度的变化规律与铁团簇扩散系数的变化规律完全对应。

2.3.2 尺寸的影响

根据斯托克斯-爱因斯坦关系式^［32］

C=\frac{k_{B}T}{D\eta R}.（4）

其中，D为铁团簇的扩散系数，\eta为液态铅铋的黏度，R为铁团簇的等效半径，C为常数。由此可见，铁团簇的扩散系数不仅与液态铅铋的黏度有关，还与铁团簇的尺寸有关。为此，保持模拟区域尺寸为22a\times22a\times22a，将一个球形铁团簇置于液态铅铋中，通过改变其半径，分析铁团簇的尺寸对其扩散系数的影响规律。体系的原子数量如表1所示。

首先考虑铁团簇的尺寸对其扩散系数的影响，计算结果如图9所示。为了对比验证，图中同时展示Gao等^［15］的计算结果，包括单个铁原子以及含有83个原子的小尺寸铁团簇的扩散系数。从图中可以看出，本文对于单个铁原子扩散系数的计算结果与Gao等^［15］的计算结果完全一致，证明了本文计算方法的正确性。还可以看出在相同温度下，铁团簇的扩散系数总是小于单个铁原子，并且随着铁团簇尺寸的增大，其扩散系数逐渐减小，这一规律与式（4）相符合。

另外，从图9可以看出，随着体系温度的升高，单个铁原子和铁团簇的扩散系数均会逐渐增大，这也符合普遍的热力学规律。扩散系数与体系温度的关系可以根据阿伦尼乌斯公式进行拟合

D=D_{\mathrm{0}}{\mathrm{e}}^{-\frac{Q}{k_{\mathrm{B}}T}}.（5）

其中，D_{\mathrm{0}}为扩散系数因子，k_{\mathrm{B}}为玻尔兹曼常数，T为体系温度，Q为扩散激活能。通过拟合可得半径为10、12.5和15 Å的铁团簇的扩散激活能分别为0.128、0.140和0.134 eV，铁团簇的尺寸虽然会影响扩散系数的大小，却对扩散激活能的影响并不显著。

3 结论

本文采用MD方法模拟铁原子在液态铅铋中的团聚过程，并分析铁团簇间的相互作用，以及铁团簇的扩散行为，旨在更好地理解铅冷快堆冷却回路中腐蚀产物的沉积过程，主要结论如下：

1）在液态铅铋中铁原子会倾向于团聚成小团簇，小团簇会接着聚集形成尺寸更大的团簇，温度升高可以加快以上过程。

2）在真空中，铁团簇间的相互作用表现为先吸引后排斥，而在液态铅铋中，受到表面吸附的铅铋原子的阻碍，2个铁团簇在相互靠近时先表现为排斥作用，因此铁团簇之间的聚合需要克服一定能垒，并且体系温度的升高和铁团簇尺寸的增大有利于2个铁团簇发生聚合。

3）铁团簇的运动扩散能力随着其质量浓度的降低而提升。铁团簇的扩散系数还与其尺寸和体系温度有关，铁团簇尺寸增大，扩散系数减小，体系温度升高，扩散系数增大。

参考文献

[1] Zhang J S， Li N. Review of the studies on fundamental issues in LBE corrosion［J］. Journal of Nuclear Materials， 2008， 373（1/2/3）： 351-377. DOI： 10.1016/j.jnucmat.2007.06.019 .

[2] Zhang J S. A review of steel corrosion by liquid lead and lead-bismuth［J］. Corrosion Science， 2009， 51（6）： 1207-1227. DOI： 10.1016/j.corsci.2009.03.013 .

[3] Ballinger R G， Lim J. An overview of corrosion issues for the design and operation of high-temperature lead- and lead-bismuth-cooled reactor systems［J］. Nuclear Technology， 2004， 147（3）： 418-435. DOI： 10.13182/NT04-A3540 .

[4] Wang W T， Yang C X， You Y H， et al. A review of corrosion behavior of structural steel in liquid lead-bismuth eutectic［J］. Crystals， 2023， 13（6）： 968. DOI： 10.3390/cryst13060968 .

[5] Zhang J S. Lead-bismuth eutectic （LBE）： a coolant candidate for gen. IV advanced nuclear reactor concepts［J］. Advanced Engineering Materials， 2014， 16（4）： 349-356. DOI： 10.1002/adem.201300296 .

[6] Gorynin I V， Karzov G P， Markov V G， et al. Structural materials for atomic reactors with liquid metal heat-transfer agents in the form of lead or lead： Bismuth alloy［J］. Metal Science and Heat Treatment， 1999， 41（9）： 384-388. DOI： 10.1007/BF02469876 .

[7] Barbier F， Rusanov A. Corrosion behavior of steels in flowing lead-bismuth［J］. Journal of Nuclear Materials， 2001， 296（1/2/3）： 231-236. DOI： 10.1016/S0022-3115（01）00521-9 .

[8] Hosemann P， Dickerson R， Dickerson P， et al. Transmission electron microscopy （TEM） on oxide layers formed on D9 stainless steel in lead bismuth eutectic （LBE）［J］. Corrosion Science， 2013， 66： 196-202. DOI： 10.1016/j.corsci.2012.09.019 .

[9] Zhu Z G， Zhang Q， Tan J B， et al. Corrosion behavior of T91 steel in liquid lead-bismuth eutectic at 550 °C： effects of exposure time and dissolved oxygen concentration［J］. Corrosion Science， 2022， 204： 110405. DOI： 10.1016/j.corsci.2022.110405 .

[10] Song C， Li D D， Xu Y C， et al. Corrosion related properties of iron （100） surface in liquid lead and bismuth environments： A first-principles study［J］. Chinese Physics B， 2014， 23（5）： 056801. DOI： 10.1088/1674-1056/23/5/056801 .

[11] Xu Y C， Song C， Zhang Y G， et al. An energetic evaluation of dissolution corrosion capabilities of liquid metals on iron surface［J］. Physical Chemistry Chemical Physics， 2014， 16（31）： 16837-16845. DOI： 10.1039/c4cp01224k .

[12] Long X L， Shi J L， Zhu J， et al. First-principles calculation of the resistance to lead-bismuth eutectic corrosion on Fe （111） surface of austenitic stainless steel［J］. Surface Science， 2022， 725： 122132. DOI： 10.1016/j.susc.2022.122132 .

[13] Zhou T， Gao X， Ma Z W， et al. Atomistic simulation of α-Fe（100）-lead-bismuth eutectic （LBE） solid-liquid interface［J］. Journal of Nuclear Materials， 2021， 555： 153107. DOI： 10.1016/j.jnucmat.2021.153107 .

[14] Chen L M， Xu S， He X X， et al. Molecular dynamics study of corrosion behavior of iron with vacancies exposed to lead-bismuth eutectic［J］. Materials and Corrosion， 2023， 74（5）： 793-802. DOI： 10.1002/maco.202213375 .

[15] Gao Y， Takahashi M， Cavallotti C， et al. Molecular dynamics simulation of metallic impurity diffusion in liquid lead-bismuth eutectic （LBE）［J］. Journal of Nuclear Materials， 2018， 501： 253-260. DOI： 10.1016/j.jnucmat.2018.01.044 .

[16] Daw M S， Baskes M I. Embedded-atom method： derivation and application to impurities， surfaces， and other defects in metals［J］. Physical Review B， 1984， 29（12）： 6443-6453. DOI： 10.1103/physrevb.29.6443 .

[17] Belashchenko D K. Computer simulation of the properties of liquid metals： gallium， lead， and bismuth［J］. Russian Journal of Physical Chemistry A， 2012， 86（5）： 779-790. DOI： 10.1134/S0036024412050056 .

[18] Zhou X W， Johnson R A， Wadley H N G. Misfit-energy-increasing dislocations in vapor-deposited CoFe/NiFe multilayers［J］. Physical Review B， 2004， 69（14）： 144113. DOI： 10.1103/physrevb.69.144113 .

[19] Gao Y， Raos G， Cavallotti C， et al. Molecular dynamics simulation on physical properties of liquid lead， bismuth and lead-bismuth eutectic （LBE）［J］. Procedia Engineering， 2016， 157： 214-221. DOI： 10.1016/j.proeng.2016.08.359 .

[20] Nosé S. A unified formulation of the constant temperature molecular dynamics methods［J］. The Journal of Chemical Physics， 1984， 81（1）： 511-519. DOI： 10.1063/1.447334 .

[21] Hoover W G. Canonical dynamics： equilibrium phase-space distributions［J］. Physical Review. A， General Physics， 1985， 31（3）： 1695-1697. DOI： 10.1103/physreva.31.1695 .

[22] Parrinello M， Rahman A. Crystal structure and pair potentials： A molecular-dynamics study［J］. Physical Review Letters， 1980， 45（14）： 1196-1199. DOI： 10.1103/physrevlett.45.1196 .

[23] Müller-Plathe F. Reversing the perturbation in nonequilibrium molecular dynamics： an easy way to calculate the shear viscosity of fluids［J］. Physical Review. E， Statistical Physics， Plasmas， Fluids， and Related Interdisciplinary Topics， 1999， 59（5 Pt A）： 4894-4898. DOI： 10.1103/physreve.59.4894 .

[24] Xu B， Liu S， Yu X. Molecular simulation study of the effect of iron clusters on the viscosity of liquid lithium［J］. Journal of University of Chinese Academy of Sciences， 2023， 40（2）： 165-172. DOI： 10.7523/j.ucas.2021.0062 .

[25] Thompson A P， Aktulga H M， Berger R， et al. LAMMPS-a flexible simulation tool for particle-based materials modeling at the atomic， meso， and continuum scales［J］. Computer Physics Communications， 2022， 271： 108171. DOI： 10.1016/j.cpc.2021.108171 .

[26] Stukowski A. Visualization and analysis of atomistic simulation data with OVITO-the open visualization tool［J］. Modelling and Simulation in Materials Science and Engineering， 2010， 18（1）： 015012. DOI： 10.1088/0965-0393/18/1/015012 .

[27] Weeks J R， Romano A J. Liquidus curves and corrosion of Fe， Ti， Zr， and Cu in liquid Bi-Pb alloys［J］. Corrosion， 1969， 25（3）： 131-136. DOI： 10.5006/0010-9312-25.3.131 .

[28] Xu C， Li Z Y， Chen L G， et al. Atomic-scale insights into the precipitation behaviors of copper atoms in liquid lithium［J］. Fusion Engineering and Design， 2023， 194： 113901. DOI： 10.1016/j.fusengdes.2023.113901 .

[29] Cui W Z， Shen Z J， Yang J G， et al. Effect of chaotic movements of nanoparticles for nanofluid heat transfer augmentation by molecular dynamics simulation［J］. Applied Thermal Engineering， 2015， 76： 261-271. DOI： 10.1016/j.applthermaleng.2014.11.030 .

[30] Gan X L， Xiao S F， Deng H Q， et al. Clustering of Fe atoms in liquid Li and its effect on the viscosity of liquid Li［J］. Nuclear Fusion， 2016， 56（4）： 046004. DOI： 10.1088/0029-5515/56/4/046004 .

[31] He Y R， Jin Y， Chen H S， et al. Heat transfer and flow behaviour of aqueous suspensions of TiO₂ nanoparticles （nanofluids） flowing upward through a vertical pipe［J］. International Journal of Heat and Mass Transfer， 2007， 50（11/12）： 2272-2281. DOI： 10.1016/j.ijheatmasstransfer.2006.10.024 .

[32] Zwanzig R. On the relation between self-diffusion and viscosity of liquids［J］. The Journal of Chemical Physics， 1983， 79（9）： 4507-4508. DOI： 10.1063/1.446338 .

基金资助

国家自然科学基金(12172356)

国家自然科学基金(U23B20110)