Mg熔体凝固过程中的分子动力学模拟*

孙 慧，坚增运，姜冰清，刘翠霞

(1.西安工业大学 材料与化工学院，西安 710021；2.西安工业大学 组织部，西安 710021)

diffusion coefficient

近年来，研究金属熔体中原子团簇与熔体结晶相在结构上的相关性对认识金属的凝固规律、控制金属凝固过程有着重要的意义[1-3].金属熔体的结构和性质直接影响到凝固过程的形核和晶体生长的特性[4-6]，但是由于金属熔体中原子团簇的结构无法进行实验测定，所以人们一直没有找到金属熔体中原子团簇与熔体中结晶相在结构上面的关系，目前的形核理论关于金属熔体中临界晶核、亚临界晶核与熔体中结晶相在结构上的关系仅仅建立在推测基础上，文献[4]曾指出研究非晶中亚临界晶核的结构及形成过程对完善和发展结晶动力学具有重要意义，但是对一般金属来说，目前实验所能获得的冷却速率远远达不到其发生玻璃化转变所需的冷速，计算机模拟技术的飞速发展为上述问题的解决提供了可能，分子动力学模拟是研究金属结构的一个有效方法，它能获得比目前实验大得多的冷却速率，可以从原子尺度直接模拟出熔体、非晶体及晶体的微观结构[7-10].

国内外涌现出了大量的关于分子动力学模拟金属熔体凝固过程中微观结构演变规律的报道[11-13]，由文献[14]采用分子动力学对液态金属的凝固过程进行了研究，拉开了以分子动力学模拟研究物质微观变化的序幕.文献[15]采用分子动力学方法对NaCl溶液进行了模拟研究，发现自扩散系数与溶液的温度和浓度有关，温度越高，自扩散系数越大；文献[16]采用EAM势对凝固过程进行了模拟研究，发现凝固以后的结构和冷速有关，当金属熔体在较慢的冷速下冷却，形成晶体，在较快的冷速下形成非晶，同时发现分析玻璃化转变温度对于分析晶体团簇结构有着重要的作用.文献[17]采用分子动力学方法对Pb-K合金的凝固过程进行了模拟研究，通过与实验数据的对比发现模拟所得结果与理论结果相符合，证明了分子动力学模拟结果的可靠性.其他一些研究者采用EAM势以及相应的微观结构分析方法对液态金属Al凝固以后的组织进行了研究.以镁为基体的材料是目前工程应用中最轻的金属材料，具有较高的比强度和比刚度，故镁金属是轻质金属材料的理想基体.随着各个学科的飞速发展，以镁为基体的合金材料的优势逐渐体现出来，然而由于镁金属的化学活性较高，并且容易燃烧同时成孔性差，所以对于金属镁熔体的性能的实验研究有一定局限性.分子动力学模拟方法不仅可以对实验难以达到和无法观察的细节进行研究，同时也能为传统实验提供科学的预判.目前采用分子动力学模拟的方法对镁熔体的团簇种类和团簇数量与冷速的关系以及亚临界晶核、临界晶核晶体在结构上的关系也并没有进行深入研究，而搞清熔体结构与结晶相结构之间的关系是正确认识形核生长特性和规律的前提.

本文采用分子动力学模拟方法研究了Mg熔体以不同冷却速率凝固后的微观结构的演变规律，并通过理论计算确定出Mg熔体凝固后所得的晶态结构为临界晶核的临界冷速，最后对不同初始温度和不同保温温度下Mg熔体的扩散系数进行了研究.

1 模拟过程及分析方法

模拟软件采用LAMMPS软件，模拟对象为13 500个镁原子，模拟采用三维周期性边界条件，选择这个边界条件主要是能够避免边界效应或表面效应对研究模型力学作用的影响，适合本模拟的研究模型，为了使模拟结果准确，故选择时间步长为2 fs，步长过大会使部分原子溢出，影响模拟结果的准确性.所选择的势函数为1990年ACKLAND G J等的势函数，由于模拟系统可以和外部的时刻进行热交换，因此模拟系统的温度能够保持恒定，采用Nose-Hoover控温控压算法来控制系统的温度压力.在常压下，将系统在980 K下运行100 000步使其达到液态平衡；然后将Mg熔体在1 400 K下运行100 000步使其转化为液态并弛豫保温500 000步，以不同冷速降温至50 K并保温弛豫，采用径向分函数、H-A键型分析和晶体团簇法分析Mg熔体的微观结构.扩散系数的模拟在1 400 K下运行50 000步使Mg熔体熔化并达到平衡，然后以1.0×1011.0K·s-1的冷速冷却至800 K，温度每降低100 K，弛豫保温5 000步测量均方位移.

2 模拟结果与讨论

2.1 熔体热历史对Mg熔体均质形核过冷度影响的分子动力学模拟

2.1.1 熔点的模拟

采用1990年ACKLAND G J等的势函数，以13 500个Mg原子作为模拟体系进行熔点的测量，绘制温度与时间的关系如图1所示.从图中可以看出，随着时间的推移，整个体系中的温度并不是一个定值，而是不断变化的.但是这种变化又有一定的趋势，总是在890 K附近波动，金属Mg理论熔点为921 K，模拟结果比实际值仅相差31 K，这在分子动力学模拟中基本可以忽略不计，说明此势函数下模拟结果可以达到稳态，因此认为金属Mg的模拟熔点Tm=890 K.

图1 采用ACKLAND G J等的势函数模拟Mg的熔点

2.1.2 初始液态平衡性验证

模拟体系的初始状态是否为液态平衡体系对于模拟结果的可靠性和准确性有着至关重要的影响，因此在开始进行一切模拟运算前都要先进行平衡性验证.图2是不同原子数的液态镁在980 K时的径向分函数.

图2 不同原子数体系中液态镁在980 K时的径向分布函数

采用原子数分别为2 048、13 500和50 000个Mg原子为研究对象，将其在980 K(处于研究对象熔沸点之间的任意温度值，为了方便后续实验的说明，这里取一个较低温度值)的初始温度下运行50 000步使其转化为液态并弛豫保温50 000步，所得体系中不同原子数的径向分布函数如图2所示.从图中可以看出，虽然体系中原子数不同，但是其径向分数的液态结构是一致的.因此可以认为Mg熔体经过100 000步的保温弛豫过程能够使液态Mg达到平衡.

2.1.3 初始温度对形核过冷度影响的模拟

将Mg熔体分别在七个不同的初始温度(980 K，1 080 K，1 200 K，1 280 K，1 350 K，1 380 K，1 400 K)下运行50 000步使其转化为液态并弛豫保温50 000步得到弛豫组织.再将弛豫组织以1.0×1011.0K·s-1的冷速快速冷却到200 K，绘制能量随温度变化曲线如图3所示.由图可知，Mg熔体以1.0×1011.0K·s-1的冷速冷却到200 K.一方面，初始温度越高，升温完成后，体系能量越高；另一方面，体系降温过程中能量均随着温度的降低逐渐减小，当温度分别为619 K、614 K、606 K、590 K、584 K、584 K、584 K时，能量出现突变，急剧减小，使原本平滑的曲线出现明显的拐点，说明当温度足够高时，能量曲线出现转折，故在此处系统开始凝固.根据热力学理论，在恒温恒压的条件下，系统总是自发地趋向低能量的稳定态，当过热温度低于固液两相的理论凝固温度时，体系中固相的自由能总是低于液相的自由能，此时体系就会发生凝固结晶这说明金属熔体在降温的过程中发生了相变，形成了晶态凝固组织，此时拐点所对应的温度为Mg熔体凝固过程中的实际相变温度.

图3 不同初始温度下Mg熔体凝固过程能量 随温度的变化曲线

为了研究过冷度随初始温度的变化关系，结合上述能量变化图获得初始温度对应的过冷度见表1.

表1 不同初始温度的Mg熔体冷却 到200 K时的过冷度

从表1可以看出，随着过热温度的升高，实际相变温度Tg大致上是逐渐减小的，其对应的过冷度ΔT逐渐增大，当温度达到1 350 K以后，Tg和ΔT都趋于稳定，不再变化，此时金属的结晶温度维持不变，相应的过冷度也将确定不变.从图4可以看出在过热温度较低时，体系内存在一些晶体质点，原子键对不易破坏，此时体系较为平衡，所以过冷度较小，但是随着过热温度的升高，金属原子的扩散能力增强，体系的平衡状态被逐渐打破，原子形核率增大，过冷度增大，曲线出现一次转折；当过热温度持续增大时，体系内原子完全挣脱束缚，在体系内任意运动，为了达到新的平衡状态，故整个系统的自由能应该在最小值附近，故而过冷度稳定在306 K.

图4 Mg熔体凝固过程中过冷度随 初始温度的变化曲线

2.1.4 冷速对均质形核过冷度的影响的分子动力学模拟及检验

将Mg熔体在初始温度为1 400 K下运行50 000步使其转化为液态并弛豫保温50 000步得到弛豫组织.再将弛豫组织分别以十个不同冷速：1.0×1011.0K·s-1,，1.0×1011.5K·s-1，1.0×1011.8K·s-1， 1.0×1012.0K·s-1，1.0×1012.3K·s-1，1.0×1012.5K·s-1，1.0×1013.0K·s-1，1.0×1013.5K·s-1，1.0×1014.0K·s-1，1.0×1015.0K·s-1(选择这些冷速主要是因为目前实验条件下所能获得的最大冷速为109K·s-1，在这个冷速条件下许多晶化能力强的金属和合金不能实现玻璃化转变，这限制了对非晶形成机理及结晶动力学的研究，故将冷速区间细化可以更好的找出临界冷速点)冷却到200 K(此温度点可以在室温下的任一温度，为了计算效率更高更快，所以选择温度为200 K)并弛豫保温后体系原子的势能的变化曲线如图5所示，在不同冷速下，随着温度的降低，体系势能总体呈现降低的趋势.当冷速降低到1.0×1012.5K·s-1后体系原子的势能变化曲线上会出现拐点，这标志着体系内部发生相变，凝固组织为晶态；当冷速大于1.0×1012.5K·s-1时，体系原子的势能随着温度的降低平缓下降不出现突变点，这表明体系内部没有发生相变，凝固组织为非晶态.因此，可以初步确定Mg熔体形成非晶的冷速应大于1.0×1012.5K·s-1.

图5 Mg熔体在不同冷速下原子势能随温度的变化曲线

为了研究过冷度随冷却速度的变化关系，结合上述能量变化图获得不同冷速对应的过冷度见表2.从表2可以看出，初始温度为1 400 K时Mg熔体以不同冷速冷却后，相变温度Tg大致上是逐渐减小的，其对应的过冷度ΔT逐渐增大.随着冷速的增大，熔体凝固后的结晶温度减小，说明过冷度受冷却速率Rc的影响较大，当冷速增加时，整个体系受到破环的温度就越低，由于受冷却速率的影响，产生了时间的滞后性，使得金属的过冷度越大.

表2 1 400 K时Mg熔体在不同冷速下的过冷度

2.1.5 模拟所得过冷度结果的可靠性验证

采用文献[18]中的理论模型计算Mg熔体凝固时的过冷度，当熔体在降温过程中发生均质形核时，液相结构中将会产生晶核，这一现象可用下面的公式来表示：

(1)

式中：NA为阿佛加德罗常数；r*为临界晶核半径；Vm为摩尔体积.

可通过下式计算出临界晶核的原子半径r*：

(2)

式中：σ为金属的固/液界面能；ΔSf为金属的熔化熵；ΔTv为凝固后的均质形核过冷度.

查出固/液界面能σ，而均质形核过冷度ΔTv可由下式确定：

(3)

其中Iv为

(4)

式中:Av为均质形核率固定值(Av=1 043±1 m-3·s-1)；ΔGA为界面扩散激活能；σ(T)为界面能；ΔHv为金属结晶潜热.如果固/液界面能知道，过冷度ΔTv就能计算出来.

对于公式中的参数值，可以通过查阅有关资料获得，并将结果绘制成表3.模拟结果是否可靠，需要同其他理论或者方法作为对比.凝固时的最大均质形核过冷度，需要将理论计算结果与模拟结果进行对比.选1.0×1011.0K·s-1和1.0×1011.5K·s-1两个冷速条件下的过冷度进行比较，其中，各项数据绘制见表3.

表3 镁在不同冷速下的均质形核过冷度

图6是模拟所得过冷度与理论对比的情况，通过固/液界面能理论模型获得冷速变化时对应的过冷度，在图6中，实线表示由界面能理论所得过冷度与冷速之间的变化关系，蓝色星状符号表示模拟所得结果.从图中对比可发现，模拟结果与理论吻合的很好，表明本次模拟的结果是准确可靠的.

2.2 熔体热历史对Mg熔体凝固组织微观结构影响的分子动力学模拟

2.2.1 径向分布函数分析结果

为了进一步确定Mg熔体凝固后形成非晶的冷速，采用径向分布函数进行分析.径向分布函数曲线如图7所示.

图6 模拟结果与界面能理论结果的对比曲线

图7 Mg熔体在1 400 K下以不同冷速降温 至200 K后的径向分布函数

从图7可知，曲线10峰型圆钝，线型平缓，具有液态的典型特征，曲线1峰型明显而尖锐，线型波动性大，说明熔体凝固后形成晶态结构.曲线2到曲线9，峰型逐渐变钝，劈裂成两个峰的第二峰逐渐消失，线型逐渐变得平缓，说明此时体系中既有非晶态的特征又具有晶态的特征.分析结果表明，冷速小于等于1.0×1011.0K·s-1时，熔体冷却降温后的结构为晶态组织；冷速大于等于1.0×1012.5K·s-1时，降温后的结构为非晶态组织，冷速在1.0×1011.0K·s-1到1.0×1012.5K·s-1之间时，熔体降温后的结构为晶态和非晶态混合组织.

2.2.2 键型指数分析

Mg熔体以不同冷速降温到200 K弛豫保温后的键对分析结果如图8所示，可以看出，当冷速为1.0×1011.0K·s-1时，表征 fcc 结构1421键型的分数为0.5，随着冷速的增加，1421键型的分数先增加后降低，当冷速超过1.0×1012.0K·s-1时降低的趋势逐渐平缓；表征 hcp结构1422 键型的分数随着冷速的增加总趋势也是降低的，但没有1 421键型降低得快，当冷速超过1.0×1012.0K·s-1时降低的趋势逐渐平缓.表征非晶结构的1551，1541，1431 键型的分数一开始为0，随着冷速的增加才开始出现并升高，表明 Mg 熔体以 1.0×1011.0K·s-1冷速降温凝固后完全形成晶体结构；冷速大于1.0×1011.0K·s-1时，表征非晶结构 1551，1541，1431 键对的分数随冷速的增加而增大，这说明随着冷速的增大形成非晶的分数增多.这与径向分布函数分析所得到的结论基本相符合.

图8 最终平衡构型中各键对分数随冷速的变化

2.2.3 最大晶体团簇中的原子数及其结构

表4为Mg 熔体以不同冷速冷却凝固后最大晶体团簇中的原子数.

表4 不同冷速下最大原子团簇中的原子数Nc及 fcc原子数Nfcc，bcc原子数Nbcc，hcp原子数Nhcp

从表4中可以看出，当冷速为1.0×1011.0K·s-1时，晶体团簇由fcc，hcp结构原子构成，其中hcp结构原子占绝大多数，fcc结构原子只有很少一部分.随着冷速的增大，最大原子团簇的原子总数，以及fcc，hcp结构原子数量都呈减少趋势.

为了能直观真实的显示最大晶体团簇的微观结构,采用周期性边界条件对研究原子进行了镜像处理,结果如图9所示.可以看出，Mg熔体以大于1.0×1011.0K·s-1冷速凝固后形成的最大晶体团簇均是由 hcp 和 fcc 晶态结构组成的混合偏聚体，其中绝大多数为hcp结构，混有少量的fcc 结构；当冷速处于1.0×1011.0K·s-1到1.0×1012.3K·s-1之间时，凝固后形成的最大晶体团簇中原子大小均匀排列紧密；以大于1.0×1012.3K·s-1的冷速冷却凝固后的最大晶体团簇，原子开始杂乱排布；随着冷速的增加,最大团簇的原子个数逐渐减少.

图9 不同冷速下最终平衡构型中最大晶体团簇结构(紫色表示hcp 晶态原子，黄色表示fcc 晶态原子)

2.3 Mg熔体扩散系数的分子动力学模拟

2.3.1 凝固过程中不同保温温度下的扩散系数

使用LAMMPS软件进行分子动力学模拟研究并不能直接得到我们想要的扩散系数，直接得到的是一个与扩散系数相关的量：均方位移MSD.根据菲克第一定律以及能斯特-爱因斯坦( Nernst-Einstein)方程，可以通过做MSD与时间的关系曲线的斜率来得出扩散系数.在1 400 K下对Mg熔体进行加热，运行50 000步，使之熔化并达到平衡；然后以1.0×1011.0K·s-1的冷速冷却至800 K，温度每降低100 K，弛豫保温500步，则可得到其均方位移.图10是模拟所得不同保温温度下的均方位移与时间的关系曲线.

由图10可以看出，一方面，在同一保温温度下，随着保温时间的延长均方位移的数值会增大，但是增大的速度会减缓，表现在图中就是均方位移曲线的斜率逐渐减小，即在同一保温温度下扩散系数会随着保温时间的延长而减小；另一方面，随着保温温度的降低，熔体中原子的均方位移逐渐减小，曲线越来越平缓，所以原子的扩散系数也就越来越小.

图10 不同保温时间下的均方位移

为了能直观的反映扩散系数与保温温度的变化规律，利用Einstein关系式，对不同温度下的扩散系数D0进行统计见表5，从表5中可以看出随着保温温度的降低，熔体原子的扩散系数逐渐降低.

2.3.2 不同初始温度下凝固过程中的扩散系数

在上面已经提到Mg熔体在980 K下就能完全熔化达到初始液态平衡结构.这里为了研究不同初始温度对Mg熔体扩散系数的影响，只需将熔体温度分别从1 000 K、1 200 K和1 400 K、以1.0×1011.0K·s-1的冷速降温到200 K，分别测定其凝固过程的MSD，从而进一步分析凝固过程中扩散系数的变化规律即可.图11为不同初始温度下的均方位移.

表5 不同保温温度下的均方位移及对应的扩散系数

图11 不同初始温度下凝固过程中的均方位移

从图11中可以看出，初始温度越高，开始降温时的均方位移增加越快，随着降温过程的进行，温度越低均方位移的变化幅度减小，扩散系数也随之减小.当时间足够长时，均方位移不再随着变化，出现一个稳定的值，此时，扩散系数也不在随着初始温度的改变再改变.

3 结 论

1) 当冷速小于等于1.0×1011.0K·s-1时，Mg熔体凝固后形成晶态组织，具有fcc结构；当冷速大于等于1.0×1012.5K·s-1时，Mg熔体凝固后形成非晶态组织；当冷速在1.0×1011.0K·s-1到1.0×1012.5K·s-1之间时，Mg熔体凝固后形成晶态结构和非晶态结构组成的混合组织，由fcc足够长时，均方位移不再变化，出现一个稳定值，此时，扩散系数不再随着初始温度的改变而改变.

2) 在一定的冷速范围内，Mg熔体凝固后得到的组织是晶体团簇和非晶组成的混合体，最大晶体团簇尺寸随着冷速的增加而减小.

3) Mg熔体中晶态结构原子团簇尺寸小于临界晶核尺寸时，用H-A键型指数法虽能确定出一定数量晶态结构原子键对的存在，但径向分布函数反映不出晶态结构的特征.

4) 在同一保温温度下扩散系数随着保温时间的延长而减小，随着保温温度的降低，熔体中原子的均方位移逐渐减小，原子的扩散系数越来越小；不同初始温度下随着初始温度的升高，开始降温时的均方位移增加较快，随着降温过程的进行，温度越低均方位移的变化幅度减小，扩散系数也随之减小，当时间位移不再随着变化，出现一个稳定的值，此时，扩散系数也不在随着初始温度的改变再改变.

参考文献：

[1] LEE B S,BURR G W,SHELBY R M,et al.Observation of the Role of Subcritical Nuclei in Crystallization of a Glassy Solid[J].Journal of Science,2009,326(3):980.

[2] 坚增运,相敏,朱满,等.Al凝固特性随冷却速率变化规律的分子动力学模拟[J].西安工业大学学报,2015,4(35):311.

JIAN Zengyun,XIANG Min,ZHU Man,et al.Molecular Dynamic Simulation of the Dependence of Solidification Behaviour of Aluminum on the Cooling Rate[J].Journal of Xi’an Technological University,2015,4(35):311.(in Chinese)

[3] 张爱龙,刘让苏,梁佳,等.冷却速率对液态凝固过程中微观结构演变影响的模拟研究[J].物理化学学报,2005,21(4):347.

ZHANG Ailong,LIU Rangsu,LIANG Jia,et al.A Simulation Study for Effects of Cooling Rate on Evalution of Microstructures During Solidification Process of Liquid Metal Ni[J].Acta Physico-Chimica Sinica,2005,21(4):347.(in Chinese)

[4] GIBSON J M.Viewing the Seeds of Crystallization[J].Journal of Science,2009,326(5955):942.

[5] 坚增运,钟亚男,许军锋,等.冷速对Au熔体凝固组织影响的分子动力学模拟[J].西安工业大学学报,2016,36(6):468.

JIAN Zengyun,ZHONG Yanan,XU Junfeng,et al.Molecular Dynamic Simulation of the Influence of Cooling Rate on Solidification Structure of Au Melt[J].Journal of Xi’an Technological University,2016,36(6):468.(in Chinese).

[6] QI D W,WANG S.Icosahedral Order and Defects in Metallic Liquids and Glasses[J].Physical Review B,1991;44(2):884.

[7] LI Y D,HAO Q H,CAO Q L,et al.Two-order-parameter Description of Liquid Al Under Five Different Pressures[J].Physical Review B,2008,78(17):174202.

[8] LIU H R,LIU R S,ZHANG A L,et al.A Simulation Study of Microstructure Evolution During Solidification Process of Liquid Metal Ni[J].Chinese Physics B,2007,16(12):3747.

[9] LI J,LIU R,ZHOU Z,et al.Microstructural Transitions During the Rapid Cooling Proess of Liquid Metal Al[J].Journal of Materials Science & Technology,1998(14):461.

[10] WOLDE P R T,DAAN F.Numerical Calculation of the Rate of Crystal Nucleation in a Lennard-Jones System at Moderate Undercooling[J].Journal of Chemical Physics,1996,104(24):9932.

[11] MONDAL K,MUTRY B S.Prediction of Maximum Homogeneous Nucleation Temperature for Crystallization of Metallic Glasses[J].Journal of Non-Crystalline Solids,2006,352(50):5257.

[12] SCHENK T,HOLLANDMORITZ D,SIMONET V,et al.Icosahedral Short-Range Order in Deeply Undercooled Metallic Melts[J].Physical Review Letters,2002,89(7):46.

[13] QIN J,GU T,YANG L.Structural and Dynamical Properties of Fe78Si9B13 Alloy During Rapid Quenching by First Principles Molecular Dynamic Simulation[J].Journal of Non-Crystalline Solids,2009,355(48):2333.

[14] ALDER B J,WAINWRIGHT T E.Phase Transition for a Hard Sphere System[J].Journal of Chemical Physics,1957,27:1208.

[15] RAHMAN A.Correlations in the Motion of Atoms in Liquid Argon[J].Physics Review,1964,136(2):405.

[16] PARK S,SCHULTEN K.Calculating Potentials of Mean Force from Steered Molecular Dynamics Simulations[J].Journal of Chemical Physics,2004,120(13):5946.

[17] NOSE S.A Unified Formulation of the Constant Temperature Molecular Dynamics Methods[J].Journal of Chemical Physics,1984,81(1):5.

[18] JIAN Z,KURIBAYASHI K,JIE W Q.Solid-liquid Interface Energy of Metals at Melting Point and Undercooled State[J].Materials Transactions,2002,43(4):721.