蒙特卡罗法模拟压力容器钢中的团簇

1、蒙特卡罗法模拟压力容器钢中的团簇（材料科学与工程学院，上海200072）摘 要：反应堆压力容器钢的脆化可能源于团簇。由于在电镜下很难识别纳米团簇，所以通过蒙特卡罗与分子动力学，模拟团簇与位错的交互作用。模拟结果表明：空位簇阻碍作用最强，其余团簇差异微小；Cu-Ni-v簇的阻力与位错线的取向有关；小尺寸Cu，Cu-v或Cu-Ni-v簇对位错阻碍作用较弱； 将Ni替换为Cu原子不会导致强度的任何显著增加。关键词：压力容器；团簇；蒙特卡罗模拟；位错Monte Carlo Simulation of Clusters in Pressure Vessel Steel (School of Materi

2、al Science and Engineering, Shanghai 200072, China)Abstract: The embrittlement of the reactor pressure vessel steel may be due to clusters. Because it is difficult to identify nanoclusters under electron microscopy the interaction between clusters and dislocations is simulated by Monte Carlo and mol

3、ecular dynamics. The results show that the latent cluster has the strongest effect and the other clusters are small. The resistance of Cu-Ni-v cluster is related to the orientation of dislocation line. The Cu，Cu-v or Cu-Ni-v cluster The hindrance is weak; replacing Ni with Cu does not result in any

4、significant increase in strength.Key words：Pressure vessel; cluster; Monte Carlo simulation; dislocation1. 团簇团簇1一般指由几个甚至上万个原子或分子形成的相对稳定的聚集体，介于孤立原子分子与宏观物质之间，用无机分子来描述显得太大，用小块固体描述又显得太小，许多性质既不同于单个原子分子，又不同于固体和液体，也不能用两者性质的简单线性外延或内插得到。因此，人们把团簇看成是介于原子、分子与宏观固体物质之间的物质结构的新层次，是各种物质由原子分子向大块物质转变的过渡状态。团簇结构的形成是由组成团

5、簇的原子数目和原子间的相互作用势决定的2。利用势函数计算原子团簇的势能的最小值，来预测团簇的最低势能时的结构。 2. 第一原理分子动力学模拟团簇2.1第一原理分子动力学原理通过第一原理分子动力学研究团簇1，一般无法严谨求解，因而在实际计算中必须采用近似。一个基本的近似就是基于Born-Oppenheimer近似，将电子的运动和原子核的运动分开，以此描述体系电子结构和总能量。在团簇原子核坐标固定的情况下体系总能量的求解常采用两种方法，密度泛函理论（Density-Functional Theory，DFT）和紧束缚（Tight-Binding，TB）方法。密度泛函理论从第一原理出发，利用自洽的方

6、法计算出系统的电子结构，可以定量的反映原子间的价键强度。但是第一原理密度泛函方法的计算量非常大，不适合计算原子数上百的大团簇体系。利用第一原理计算团簇稳定构型过程中，并没有充分考虑团簇电子结构的性质。而紧束缚方法通过一套物理参数来描述原子间的价键关系。由这些物理参数组装体系电子部分的大哈密顿矩阵，然后在原子核固定的条件下由量子力学的方法求出体系的电子结构，其计算量比密度泛函小 2-3 个量级。根据第一原理（密度泛函和紧束缚）分子动力学，在团簇原子核坐标固定的情况下，计算出体系的电子结构和总能量。再利用 Hellmann-Feynman 定理对电子总能微分就可以求出每个原子核感受到的力。原子核在

7、该力作用下按照经典力学方程运动，从而可以进行结构优化、动力学模拟等计算。2.2第一原理分子动力学的优点当团簇包含的原子个数较少时，可结合基于密度泛函的第一原理分子动力学方法进行计算；而对于更大的团簇体系， 为了进一步减小计算量，还可以使用Tight-Binding 方法1。由于密度泛函方法的计算量非常大，不适合计算大团簇体系，因此在计算大团簇体系的时候我们采用基于第一原理的紧束缚方法。由于紧束缚方法的精度取决于所用物理参数的有效性，只要有合适的参数，紧束缚方法就能用较小的计算量研究非常大的团簇体系。基于第一原理的紧束缚方法和密度泛函方法相比较，不仅保持了第一原理的精度，而且可以以较小的计算量计

8、算较大的团簇体系。紧束缚方法在对原子数上百的碳纳米管的计算中取得了很好的效果。3. 蒙特卡洛法模拟团簇3.1蒙特卡罗法原理蒙特卡洛法又称随机抽样法，基本原理2：对影响可靠度的随机变量进行大量的随机抽样，然后把大量的抽样数据代入功能函数，去确定结构失效概率式。由概率定义知，某事件的概率可以用大量试验中该事件发生的频率来估算，当样本容量足够大时，可以认为该事件的发生频率即为其概率。因此，可以先对影响其可靠度的随机变量进行大量的随机抽样，然后把这些抽样值一组一组地代入功能函数式，确定功能函数是否失效，最后从中求得函数的失效概率。3.2 Metropolis蒙特卡罗法3.2.1 算法原理Metropo

9、lis蒙特卡罗算法是一种权重方法（或称随机抽样方法）。在热平衡情况下，对于用极限频率分布ens() 表征的特定统计系综，利用这种方法可生成系统状态3。在化学组分恒定的系综（正则、微正则和等温等压系综）中，其概率分布是系统哈密顿量的函数；而在化学组分变化的系综（巨正则系综）中，概率分布则是化学势的函数。对于由N个经典粒子组成的系统，其热力学状态可以用位置矢量（r1 ，rN ）的三个分量和动量矢量（p1 ，pN ）的三个分量描述。这 6N个分量可以看作矢量在6N维相空间中的分量。3.2.2.正则和微正则系综的Metropolis方法3 在正则和微正则系综中，蒙特卡罗方法通过任意地或系统地改变粒子位

10、置可以给出新的组态。例如，把原子的位置由i变到j就代表系统状态的变化，并记为i j，在相空间它表示系统由i点运动到j 点。根据所考虑粒子之间的相互作用，这一点位移将使系统能量由 H( i )变为H( j) 。然而，原子的这一位移并不影响系统的组成。计算出新组态相关的系综哈密顿量变化H(ij)，即H(ij)= H( j) - H( i ) 式中，H( j)为试探态的哈密顿量；H( i )为初始态的哈密顿量。然而，根据下面介绍的定则，可以计算出新的原子排列的情况。如果新的能量值比变化前的能量值小，则该位移使得系统处于较低能量状态。这时，变更即被接受且其位移原子保留在新的位置上。如果新的能量值大于前

11、一能量值，则该变更仅以一定的概率ens(j) /ens(i) 被接受，在正则系统中这个概率由玻耳兹曼因子给出，即：ensjensi=exp-Hij式中， =1/kB T。当H(ij) 取正值时，接受组态变化的概率为：Pij exp-Hij 根据Metropolis等人的理论，若假定在0到1之间由一个随机数 ，则新的组态可按照下列规定确定：当exp- H(ij)时，接受变更，否则不接受变更。如果新的组态被驳回，亦即没有被接受，则把原位置记为新位置并记数一次，然后采用其它随机选定的原子重复上述过程。3.2.3巨正则系综的 Metropolis方法3 巨正则系综中，其基本组分的浓度是没有限定的。这样

12、，就可任意选择一个原子并通过改变一次原子种类而得到一个新的组态j。显然，这一过程将影响到系统的化学组成。同样，组态的变更将以一定的概率被接受。在巨正则系综，因其组成变化而引起的能量变化U(ij)，也可由计算得到。对新的组态可按照下列标准进行判断：如果得到的新的能量小于原来的能量，则组态的这一变更被接受；如果能量大于原来的能量，则被一定的概率接受。当U(ij)为正值时，被接受的概率：Pij exp-Uij （2.16）式中，U 表示混合能及化学势该变量之和。判断的依据与正则和微正则的方式相同。3.3蒙特卡罗法优点MC方法，主要适用于原子尺度和显微尺度的模拟。MC 法用随机数来控制粒子运动，并使其

13、符合 Boltzmann 分布。MC方法并非只能用来解决包含随机过程的问题，例如:对计算定积分的问题可将其转换成相关的随机过程， 然后用Monte Carlo方法进行求解。MC与传统数学方法相比，具有直观性强，简便易行的优点，能处理一些其他方法无法解决的问题，并且容易在计算机上实现，在很大程度上代替许多大型、难以实现的复杂实验。3.4蒙特卡罗法在材料中的应用Monte Carlo方法在材料中的应用主要分为4：经典例子体系中的应用，如团簇结构、薄膜生长、合金偏析和高分子材料等；经典自旋体系，如Ising自旋玻璃、一维下的热平衡、二维下的二级相变。图1 MC法在经典例子体系中的应用3.5 蒙特卡罗

14、模拟的步骤： 根据欲研究的物理系统的性质，建立能够描述该系统特性的理论模型，导出该模型的某些特征量的概率密度函数（即构造或描述问题的概率过程）； 从概率密度函数出发进行随机抽样，实现从已知概率分布的抽样，得到特征量的一些模拟结果；有了明确的概率过程后，为了实现过程的数值模拟，必须实现从已知概率分布的随机数的抽样，进行大量的随机模拟实验，从中获得随机变量的大量试验值。产生已知概率分布的随机变量，是实现MC方法的关键步骤，其中最基本的是（0，1）均匀分布。 对模拟结果进行分析总结，预言物理系统的某些特性。 模拟结果的检验。对于MC方法模拟纳米团簇的具体算法4见图2。图2 MC法的具体算法流程图4.

15、 具体应用：刃型位错与bcc铁中的Cu-Ni空位簇的相互作用54.1案例背景反应堆压力容器钢脆化的硬化原因6：铜高度不溶于铁，在中子照射下导致团簇（2nm）形成，富铜的析出物引起辐照硬化。纳米空位簇和位错环造成的基质损伤。由于不能检测低于TEM分辨率的位错环，无法论证不同类型钢中的位错环的存在7，位错与特定合金元素的关联不能通过实验明确地建立。所以通过蒙特卡罗与分子动力学，考虑它们与移动位错的相互作用。研究刃型位错与包含铜原子，镍原子和空位的纳米团簇以及纯纳米空位簇和纯铜簇的相互作用6。bcc铁中的空位簇、Cu簇与位错的相互作用MD研究，已经表明：只有大的纯Cu团簇（> 3nm）作为位错

16、的强障碍物；可见的位错环（2nm）可能作为位错的强障碍；低于TEM分辨率的位错环在与位错相互作用时自发地改变它们的伯氏矢量并且被立即吸收；温度对位错与1nm以下团簇的相互作用几乎没有影响；空位簇是比纯Cu团簇更强的障碍物，小于1nm的空位簇容易被位错剪切。4.2模拟方法为了确定给定组成的每个簇中的溶质和空位的能量最有利的排列，我们使用Lattice Metropolis Monte Carlo（LMMC）采样技术获得复杂簇的原子构型。LMMC计算在包含根据bcc结构排列的2000个原子的立方体盒中进行，其中主轴沿着100，010和001方向，并且空位，Ni和Cu原子初始随机分布。不同晶格位置处

17、不同类型的原子交换导致总能量的变化，在晶体中产生的三元原子间势。模拟中，拒绝总能量增加的步骤（即原子位置交换），运行总能量减少的步骤，直到晶体的总能量收敛到恒定值。随后，将团簇模型输出到晶体模拟盒中的位错滑移面上。基于探索不同的溶质/空位比率以及用Ni原子替换部分Cu原子时的化学效应的想法，选择了许多团簇模型，见表1，引入纯Cu簇和空位簇，以便参考比较。用于确定簇中的原子和空位的空间分布的程序如下：使用原子离析分析，耦合结构最近邻分析，以及通过选择具有高势能的原子来监测模拟期间位错核、团簇的位置和结构。表1 Cu-Ni-v簇的构型及其总结合能序号空位数铜原子镍原子团簇大小能量差11020030

18、-13.4421040050-20.5331060070-27.1340151530-6.012510301050-18.59610303070-18.807700070-73.368500050-46.679300030-23.9210070070-21.4611050050-13.8212030030-7.5晶体(bcc)模拟盒的主轴x，y和z分别沿着111、2和10。滑移面为(110)，位错线与y方向平行，并且具有平行于x轴的伯氏矢量=12 111。 沿x和y方向应用周期边界条件。沿着z，盒子被分成三部分：上部和下部由几个原子平面组成，其中原子刚性地固定在其原始位置，而内部区域中的原子在

19、MD运行期间自由移动。 位错的动力通过在x方向上的块体刚性移动而产生8。变形引起的剪切应力为s= Fx/Axy，其中Fx是由于内部区域中的原子在下部刚性块上沿x方向的总力，Axy是x-y截面面积。为了为该弯曲提供足够的空间，MD箱的内部区域的尺寸为沿着x，y和z有100×3、30×6和25×2个原子平面。因此，盒体积为约25×21×10 nm3，并含有约50万个移动原子。首先放置直径为约1nm的Cu-Ni-v或Cu-v簇，其中心位于位错滑移平面上，离位错线约5nm，如图3所示。图3 MD模拟设置的示意图；刃型位错位于含有约50万个原子的晶体中用

20、MC获得的Cu-Ni-v簇放在滑移面上；团簇包含70个Cu原子和10个空位4.3结果分析4.3.1 互动机制和一般趋势最初，位错滑向团簇，产生交互作用，应力增加，位错线弯曲，最终从团簇断开，在团簇出口表面上留下台阶6。图4可以看出 ，空位簇阻碍作用最强，其余团簇差异微小。观察到位错和沉淀物之间的反应通过相同的机理发生，而与团簇的化学组成或尺寸无关。图4 刃型位错通过空位簇，Cu70，Cu60v10和Cu30Ni30v10簇的剪切应力图5中显示了纯纳米空隙和纯Cu沉淀的断裂处位错的关键构型。图b，与空位簇的相互作用导致位错偶极子的形成，并且临界角接近零；Cu团簇（cos77°= 0.2

21、1） 空位簇（cos18°= 0.95）。在与空位簇的相互作用中，观察到从团簇延出的两个螺旋臂爬升，因此位错运动中吸收了一些空位。空位的总数与纳米空隙直径成正比;纳米空隙越大，空位被吸收的越多。图5 位错与V70和Cu70簇在（110）平面交互后的位错核（a）显示了模拟框的实际尺寸；（b）障碍物附近的位错的构造被放大4.3.2 Cu-Ni-v簇的阻力发现剪切给定的Cu-Ni-v簇所需的临界应力取决于位错线取向。 图6中提供了取向对克服Cu60v10簇所需的临界应力的影响的实例。即使临界应力非常小，我们可以看出，一个特定的方向提供比其他两个明显更高的阻力。然而，该差异倾向于随着温度的升

22、高而消失。即使在与Cu-Ni-v簇的反应中没有观察到螺旋偶极子，存在于簇中的部分空位被刃型位错吸收。 发现实际吸附的空位数取决于簇的具体取向，其大小和模拟温度。 通常，温度越高，吸附的空位越多。 因此，这些团簇的有效尺寸在与刃型位错的相互作用之后减小，团簇发生原子重排。通过提高温度，Cu-Ni-v簇的临界应力也降低。然而，由于由热原子振动引起的分辨应力的波动，这些团簇的缺点导致这种减少难以发现。 图6，取向对克服Cu60v10簇所需的临界应力的影响。即使临界应力非常小，11 表现出更高的阻力。然而，该差异倾向于随着温度的升高而消失。图6 位错通过Cu60v10簇的剪切应力；三个曲线对应于微晶旋转后新的x轴分别与100,010或001方向一致时的情况图7a和b分别显示了在此研究的最小和最大团簇的数据。 在所有温度下，所有Cu-Ni-v簇的强度小于或等于相同尺寸的纯Cu簇的强度。 因此可以得出结论：所有小的Cu，或Cu-v或Cu-Ni-v簇对位错阻碍作用较弱， 将Ni替换为Cu原子不会导致强度的任何显着增加。图7 模拟中不同簇对应的剪切应力5. 结论使用LMMC技术获得低于电