Zr纳米粒子内HCP-BCC堆积结构转变的分子动力学模拟

王之涵， 于 宁， 张 林，2，

(1. 东北大学 材料各向异性与织构教育部重点实验室， 沈阳 110819； 2. 东北大学 轧制技术及连轧自动化国家重点实验室， 沈阳 110819； 3. 东北大学 材料科学与工程学院， 沈阳 110819)

1 引 言

纳米粒子随温度升高所发生的结构转变与块体相比存在明显的不同[1]. 自1909年Pawlow提出纳米粒子的熔点与其尺寸成反比关系以来[2]，研究者们对于具有自由表面的纳米粒子的熔化问题提出了多种热力学模型[3-5]并对多种金属纳米粒子进行了研究[6-13]. 这些研究发现，纳米粒子的尺寸显著影响着粒子的熔化行为. 小尺寸的纳米粒子倾向于堆积为稳定的二十面体(Icosahedron, Ih)构型，大尺寸粒子内的原子则堆积为与其体相相同的结构. 对于较大尺寸或大尺寸的纳米粒子，在它们发生结构转变前，粒子的表面在较低温度下就已率先发生了原子堆积结构的重排，且随着温度的升高，这种结构的重排持续进行，并缓慢向粒子内部扩展[14-17]. 这种发生在熔化之前的结构转变极大地影响着纳米粒子的熔化过程. 阐明这些粒子在发生固-液相变前所出现的结构转变对于理解纳米粒子熔化的发生过程，制备具有特定结构与功能的纳米材料极为重要. 然而与体相具有FCC结构和BCC结构的金属纳米粒子相比，对体相为HCP结构的粒子的研究则少得多. Ti和Zr是两种典型的六角结构金属，它们的体相在升温过程中会发生由HCP结构向BCC结构的转变. 对Ti纳米粒子升温过程的原子模拟表明，当粒子内所包含原子数目增加到400个以上时，在升温过程中该粒子会发生由HCP向BCC堆积结构的转变，进而出现HCP和BCC结构共存的现象[18-20]. 但目前对于Zr纳米粒子结构转变过程的模拟尚不多见.

位于周期表第5周期ⅣB族的过渡金属元素Zr，在室温下为密排六方晶体相，具有耐腐蚀、密度低和热中子吸收截面小等优良性质，是核反应堆堆芯的关键结构材料[21,22]. 近年来随着核反应堆技术的发展，研发满足新型堆芯高燃耗、长寿命要求的新型锆基合金成为核电用材料研究的重点[23,24]. 这样就需要对Zr纳米粒子的结构与性质随温度和尺寸的变化提供原子尺度上的理解，这种理解可以架起从微观的小尺寸原子聚集体到宏观体相材料结构转变的桥梁，从而为制备高性能锆基合金提供理论参考.

基于嵌入原子势(Embedded Atom Method, EAM)的分子动力学方法适合于研究金属体系内难以进行实验表征的微观结构转变和演化[25,26]. 本文采用分子动力学方法模拟研究升温过程中Zr纳米粒子内原子堆积结构的变化. 在模拟中，粒子形状的变化与其内部原子位置所决定的转动惯量有关，由于温度及粒子尺寸不同所导致的粒子内原子堆积结构变化的差异则由对分析技术以及原子堆积结构二维投影图显示. 通过对势能、对分布函数、势能差分、粒子形状、粒子内部堆积结构以及堆积图像的综合分析，能够确定不同粒径纳米粒子的HCP-BCC转变路径及转变温度.

2 模拟方法

本文中Zr原子间的相互作用势采用了Zhou等人提出的EAM势的形式[27]. 在EAM所包含的两体项和多体项中，两体项用于预测在合适的Cauchy压力下具有准确晶格参数的可均匀承受应力的晶格，多体项则主要影响缺陷附近处的弛豫结构. 作用在每个原子上的作用力由势能梯度的负值给出，通过求解牛顿运动方程，得到每个原子的位置与速度. 在本文的正则系综(NVT)模拟中使用Andersen热浴，并运用预报-校正算法对运动方程进行积分.

模拟系统的总势能与原子对间距密切相关. 在模拟过程中原子对间距会随着原子堆积结构的变化而变化，这将导致总势能在模拟过程中会随着时间步的增加不断改变. 图1显示了在300K的温度下，Zr13、Zr57、Zr401、Zr763、Zr1477这五个粒子的原子总势能随模拟时间的变化. 由图中可见，Zr13和Zr57的总势能均呈现出起始相对较高并随模拟时间的增加而快速降低，最后在较低值振荡变化；而大粒径的Zr401、Zr763、Zr1477的总势能变化较为平缓. 这是由于所构建的粒子初始结构为理想状态，粒子内的原子占据着标准HCP结构的格点位置，当所构建的粒子处于某一温度进行弛豫时，对于大尺寸粒子，仅少量表面原子会发生位置变化，内部的大部分原子仍堆积成标准HCP结构，当粒子由低温升温到室温时，这些粒子内的原子堆积结构已经经过了较为充分的结构弛豫，故大尺寸粒子在模拟开始后的势能变化没有出现大幅度的振荡. 而对于小尺寸粒子，温度变化对原子间距的影响则明显的多，这时尽管粒子内的原子堆积结构没有发生明显改变，但位于表面的原子，它们之间以及和其近邻的次表面原子之间的间距会较低温时明显增加，这就使得室温下模拟开始后的总势能出现明显的降低. 图2显示了这五个粒子在300 K的原子堆积结构图，可以看出粒径较小的Zr13和Zr57粒子在室温下堆积为Ih结构，而粒径较大的Zr401、Zr763和Zr1477粒子内的原子则堆积为HCP结构. 由于小尺寸粒子的Ih结构具有很高的稳定性，其升温过程的结构变化较为单一，所以下面主要对粒径较大的Zr401、Zr763和Zr1477粒子所发生的堆积结构转变进行详细的计算研究.

图1 部分Zr纳米粒子的原子总势能随模拟时间的变化Fig. 1 The potential energy per atom of Zr nanoparticle with the simulation time.

图2 上述五个粒子在300 K时的堆积结构二维投影图Fig. 2 Two-dimensional projection views of the atomic packing structure of the above five particles at 300 K.

由于总势能需经过一定时间的驰豫后才能稳定在一定范围，本文所提及的平均势能是选取每个温度点模拟的最后50000个时间步的总势能求平均而得到的，这个时间步数的选取要大于所计算的粒子系统的自由度数3N，这里N为粒子所包含的原子个数.

对分布函数g(r)给出了在粒子中找到距离为r的原子对的概率:

(1)

式中< > 表示系综平均，当两原子间距等于rij时，δ为1,否则为0.

为了分析纳米粒子的形状随温度的变化，需要计算粒子的形状因子. 首先，以纳米粒子的质心为原点建立坐标轴，然后计算粒子的转动惯量张量Iij：

(2)

(3)

式中m为Zr原子质量，δij为单位张量，xl为任意原子到粒子质心的距离，当i或j取1、2、3时对应的xi或xj分别表示任意原子到x、y、z轴的距离，式中的求和遍及所有原子. 在得到Zr纳米粒子的惯量张量以后，将其对角化，得到三个主转动惯量值. 然后将主转动惯量中的最大值和最小值分别记为I3和I1，这样便可以定义形状因子：

(4)

显然形状因子越接近1，纳米粒子越接近球形；形状因子越接近0，纳米粒子越接近线状长条形. 通过计算各个温度下的形状因子，可以得出Zr纳米粒子外形的变化.

图3 各特征键对的示意图，白色圆圈表示成键原子, 黑色圆圈表示两个成键原子的共有近邻原子Fig. 3 Schematic diagrams of various characteristic atomic pairs. The white circle represents the bonding atom, and the black circle represents the neighboring atoms shared by two bonding atoms.

对结构转变温度下各粒子内的堆积结构随模拟时间变化的分析可以通过键对分析(Pair Analysis,PA)来进行[28,29]. 它使用四个整数i，j，k，l来描述存在于各种相中的不同类型的粒子对，通过统计各特征键对的数目就能够表征粒子内不同的堆积结构.i=1表示这两个原子成键；i=2表示两个原子不成键；j表示这两个原子周围共有成键原子的数目；k代表这两个原子的共有成键原子间的成键数目；l是为了确定唯一的某一类键对而另外指定的参数. 其中二十面体结构以1551键对为特征键对；FCC结构以1421为特征键对；等量的1421和1422键对表征了HCP结构；等量的1441与1661键对用于表征BCC结构.

3 结果与讨论

可以把纳米粒子在两个相邻模拟温度下的原子平均势能差值除以这两个温度点的温度差后所得的数值定义为原子平均势能的差分值，它可以描述原子平均势能在不同温度范围内的变化趋势，从而给出纳米粒子结构转变温度的确切范围. 图5给出了Zr401，Zr763和Zr1477这三个粒子在升温过程中的原子平均势能差分值随温度变化的曲线. 如图中所示，当温度低于400 K时，这三个粒子的平均势能差分值具有相近的数值，随着温度的升高，平均势能差分值出现振荡变化，并各自在某一温度点处出现势能差分值的陡增. 其中，Zr401和Zr763粒子的势能差分值在发生陡增前的变化更为显著，这是由于这两个粒子相较于Zr1477具有很高的表面积/体积比，其表面原子的结构重排对于势能变化的影响也更大. 这里，Zr401粒子的势能差分值变化尤为剧烈，这意味着该粒子在更大的温度区间内存在着多种结构转变. 当温度分别为1250 K、1350 K和1450 K时，Zr401，Zr763和Zr1477粒子的势能差分值出现了剧烈变化的峰值，这表明，在1250-1300 K、1350-1400 K和1450-1500 K温度范围内，这三个粒子内发生了明显的原子堆积结构转变. 进一步通过二分法寻找结构转变温度的计算表明，这三个粒子的结构转变温度分别为1253 K 、1361 K和1475 K. 对于具有较小粒径的Zr401粒子，在725 K和828 K两个温度下也能观察到明显的结构转变. 得到这些温度点后便可以对各粒子在结构转变温度下的驰豫过程中的结构演变进行计算研究.

图4 部分粒子在300 K时的对分布函数和对应的堆积结构Fig. 4 The pair distribution functions and the atomic packing structures of some particles at 300 K.

图5 原子平均势能差分值随温度变化曲线Fig. 5 Atomic average potential energy difference value varying with temperature.

图6显示了Zr401粒子在这三个转变温度下的原子平均势能和形状因子随模拟时间的变化. 由图中可见，在725 K和828 K的较低温度下，尽管该粒子的势能没有出现跃变，但其形状因子却出现明显变化. 当温度为725 K，模拟时间超过300 ps后，该粒子的形状因子振荡下降，这表明该粒子由近球形被拉长；当温度升高到828 K时，随着模拟时间的增加，粒子的形状则由300 ps时的杆状逐渐恢复为近球形，并在600 ps后，形状因子呈现出小幅振荡变化. 势能和形状因子变化趋势的不匹配表明，在这两个较低的温度下，粒子并没有在整体上发生结构转变，但在粒子某些区域内的原子堆积结构却发生了变化. 这种转变虽然不会明显改变粒子的势能，但会显著影响粒子的外形. 当温度达到较高的1253 K，模拟时间少于1040 ps时，势能在低势能值附近振荡变化，随后，在较短的时间内，势能出现由低势能值向高势能值的跃升，并在1120 ps-1392 ps的时间内维持着这个高势能值. 之后，势能值再次出现明显的骤降. 这是由于Zr401粒子粒径较小，有很大比例的原子位于粒子表面，而这些表面原子比内部原子更易于运动，对于具有较小粒径的粒子而言，温度足够高时表面结构的变化极易引起整个粒子内原子堆积结构的转变. 在1253 K温度下，随着驰豫过程的进行，Zr401粒子内的原子堆积结构发生了两次明显的转变. 与此相对应，粒子的形状在这两次堆积结构转变的时间范围内也发生了明显的变化.

随着Zr纳米粒子所包含的原子数增加到763，由图7的势能变化曲线可以看出，尽管势能没有出现明显的跃变，但在650 ps之后势能值高于650 ps之前的势能值，且整体上呈现出振荡变化，与此同时，该粒子的形状因子也发生了剧烈的变化. 这表明，在该温度下，粒子内的原子堆积结构正不断发生着转变. 当粒子内所包含的原子数目增加到1000个原子以上时，随模拟时间的增加，会出现势能值的明显跃升. 对于Zr1477粒子，结构的转变发生在480 ps-560 ps的短暂时间内，在这个时间间隔内，粒子的形状因子出现了显著的振荡变化，随后形状因子的快速下降预示着粒子被拉长. 显然，这些粒子在结构转变温度点处的驰豫过程存在着多种结构转变，这种结构的转变下面通过键对分析来进一步观察.

图6 Zr401在各结构转变温度下的原子平均势能随模拟时间的变化和对应的形状因子随模拟时间的变化. (a)势能；(b)形状因子；Fig. 6 The potential energies of Zr401 nanoparticles with simulation time at various structural transition temperatures and the corresponding shape factors with simulation time. (a) The potential energy; (b)The shape factor

图7 Zr763和Zr1477在结构转变温度下的势能随模拟时间步的变化以及对应的形状因子随模拟时间步的变化(a)总势能；(b)形状因子；Fig. 7 The potential energies of Zr763, Zr1477 nanoparticles with the simulation time at the structural transition temperature and their shape factors with the simulation time. (a) The potential energy； (b)The shape factor

图8 Zr401、Zr763、Zr1477在各形状转变点处的键对随模拟时间的变化以及相应的堆积结构二维投影图；Fig. 8 The pair fractions of Zr401, Zr763, and Zr1477 at some shape transition points with simulation time and the corresponding two-dimensional projection views of the atomic packing structure.

图8显示了Zr401粒子在725 K、828 K和1253 K，Zr763粒子在1361 K，Zr1477粒子在1475 K下的各键对所占比例随模拟时间的变化，图的右侧给出了这三个粒子在不同温度、不同模拟时间下的原子堆积结构图. 如图中所示，对于Zr401粒子，当温度为725 K时，粒子在320 ps时主要存在着等量的1421和1422键对，之后1421键对的比例上升而1422的比例下降，达到640 ps时两键对的比例又趋于相等，这意味着320 ps-640 ps间粒子内部分原子存在着HCP-FCC堆积结构间的相互转化，随着模拟时间延长至1440 ps后，1422键对比例不断降低，此时粒子中出现了以HCP为主的HCP-FCC共存结构；随着温度升高到828 K，在320 ps下，1421键对比例略高于1422键对，当模拟时间达到480 ps后，1421键对比例上升并维持着一个较高值，1422键对比例下降并维持着一个较低值，同时由图中可见，在该温度下还存在着一定比例的1441、1661键对以及1551键对，并随着模拟时间的增加，1441和1661键对比例趋于相等. 这表明828 K下粒子内大量原子由原堆积成的HCP结构转化为FCC结构，并最终形成以FCC结构为主的FCC-HCP-BCC-Ih多结构共存状态；当温度升高到1253 K，模拟时间为640 ps时，粒子内相当数目的原子已经堆积成BCC结构，同时还存在着一定数量的Ih结构，这种结构的共存一直持续到1040 ps，随后1551键对所占比例快速增加的同时，1441和1661键对所占比例出现明显下降. 这表明随着模拟时间的增加，粒子内部分堆积成BCC结构的原子不断转化为Ih结构，最终形成以Ih结构为主，并和一定数量的BCC结构共存的混合结构. 对于较大粒径的Zr763粒子，在原子堆积结构转变为BCC前，粒子内的原子主要堆积为HCP结构和FCC结构，同时也存在着少量的Ih结构和BCC结构，并在480 ps-592 ps的时间内出现FCC结构和Ih结构的增长. 随后，在FCC结构出现大幅降低的同时，BCC结构出现显著增长，最后形成BCC结构和少量HCP、 Ih结构的共存.

5 结 论

本文采用基于嵌入原子势的分子动力学方法，研究了不同粒径Zr纳米粒子在升温过程中的结构转变，并着重讨论了在结构转变温度点处驰豫过程中的结构转变路径. 模拟结果表明，粒子内原子堆积结构的转变并不是两个结构间的简单转化，而且该过程还会受到粒径的影响. 对于较小粒径的Zr401粒子，在发生固-液转变前，粒子内的原子堆积结构发生了三次转变，即低温下首先从HCP结构转变为HCP-FCC共存，然后由HCP-FCC共存转变为FCC-BCC-HCP-Ih共存，随着温度升高逐渐转变为BCC-Ih共存，并最终形成Ih结构和少量的BCC结构共存的样式. 随着粒子粒径的增大，多结构的共存仅在较高的温度下出现，而且BCC结构的出现主要来自于粒子内原堆积为FCC结构的原子的位置变化. 大粒径粒子在完成BCC结构转变后，粒子内会出现明显的界面区.