飞秒激光照射金纳米颗粒的分子动力学模拟

2019-06-06牛泽伟

上海理工大学学报 2019年2期

牛泽伟，李凌，关阳

（上海理工大学能源与动力工程学院，上海 200093）

激光烧结技术是近年来发展起来的先进加工制造技术，它的快速、高精度、节能等优点使得该项技术在各个领域得到广泛应用[1-2]，并且得到越来越多的关注和研究。连续介质学方法是一种在宏观尺度对激光烧结过程进行研究的数值模拟方法，双温度模型（TTM）是其中典型并且发展较快的一种。该模型能够模拟激光与金属颗粒的相互作用与导热过程，可以很好地描述颗粒内部两步加热过程及电子和晶格之间的不平衡传热的特点[3]。随着激光参数的进一步减小，如飞秒激光脉冲，它的能量密度很大，时间和空间上的尺度又很小，这使得烧结过程中的传热现象呈现出不一样的特性。分子动力学方法（MD）能够从微观层面描述每个原子、分子的运动轨迹，从而很好地描述微观原子在高度非平衡条件下的相变过程。同时，还能将过程中包括压力因素等一些常规方法忽略的因素考虑进去，能更加全面地观察激光与纳米颗粒相互作用的过程机理[4]。但经典的分子动力学是基于电子晶格热平衡的假设之下推导而来的，不能准确地描述激光与材料作用时电子晶格间的非平衡传热过程。本文兼顾了双温度模型与分子动力学模型各自的优点，采用在分子动力学方法的基础上引入双温度模型的方法，对颗粒内部电子和晶格不平衡相变传热过程进行了模拟分析。

1 计算模型和数学方法

1.1 计算模型

现主要以单个金纳米颗粒为对象，应用分子动力学方法和双温度模型耦合的方法研究其在飞秒激光照射下的传热及相变过程。考虑现实中激光照射到粉末床时，每一个颗粒表面都会发生反射现象，因而假设颗粒表面的热流密度均匀[5]，如图1（a）所示。颗粒初始温度为300 K，半径R=1.2 nm，原子的初始位置满足面心立方晶格结构，共包含457个原子。金的势函数采用EAM势[6]，颗粒周围均采用自由边界条件。

激光脉冲随时间满足高斯分布，如图1（b）所示，激光半高宽度tp=200 fs，模拟过程从-2tp时刻开始，时间步长为2.5 fs，整个模拟过程为25 ps，激光能量密度J=1 J/m2。

图1 计算模型Fig. 1 Computational model

1.2 MD和TTM的耦合

对于双温度模型（TTM）和经典分子动力学方法（MD）本文不作介绍，而是直接说明MD和TTM的耦合方法。加入激光作用之前使颗粒内部原子处于平衡状态[7]。

纳米粒子中，由于粒子尺寸小于平均碰撞自由程，热传导方程不适用，粒子间的能量传递由原子之间的相互碰撞完成。加入激光作用后，由于分子动力学中并没有电子的模型，因此，在与双温度方程结合时，不能直接地表现出电子温度的变化以及电子将能量传给晶格的过程。本文采用的办法是，先利用双温度模型计算出在每一层中电子传入晶格的能量，这部分能量就相当于在原子原来所受到的力上添加了一个额外的作用力Fi，即[8]

式中：mi和ri分别为第i个原子的质量和位置；′为该原子热运动速度，是原子的运动速度与该原子所在层运动速度的差值；t为时间；ξ为修正系数。

ξ的计算公式为[9]

式中：G为电子-晶格的耦合因子[3]；为每一层的体积； Tl为每一层中晶格的平均温度；Te为电子温度，i为原子的序数；n为模型中原子的总个数；ΔtMD为分子动力学方法中的时间步长；ΔtFD为双温度模型中的时间步长。

1.3 颗粒内原子固液状态的判定

使用由Morris提出的序参数法[10]来区分固态原子、液态原子，从而确定颗粒发生相变程度和分析其中的熔化特点。每一个原子的序参数

式中：rij为小于截断半径rcut时原子i和j之间的距离；h是同原子i最近邻原子总数；qk为一组倒易矢量，其取值不唯一，但必须满足i为虚数符号，文献[11]给出了一组最佳的qk。

在截断半径范围内，每个原子和其相邻的h个原子的序参数平均值

序参数平均值越大，说明原子越有序。理想固态晶格的序参数为1，颗粒完全发生相变的原子序参数接近0。取平均序参数0.04作为区分固态液态的标准[11]。

2 结果与讨论

2.1 算法的验证

首先对算法和程序的正确性进行验证。对激光（tp=200 fs，J=0.2 J/m2）照射半径 r=1.2 nm 的金纳米颗粒的整个过程进行模拟。未加入激光前，当颗粒内部平衡时，用能量均分定理计算得到颗粒温度T=300.5 K。此时给颗粒加入激光能量2.77×10-18J，到达新的平衡时，颗粒温度T=445.5 K。计算得到固态金在温度300.5～445.5 K之间的比热容是127.9 J/（kg·K），这与文献[12]中金的比热容 128 J/（kg·K）相吻合。

2.2 飞秒激光照射金纳米颗粒的相变过程

对tp=200 fs，J=1 J/m2的激光照射半径为1.2 nm颗粒过程进行了模拟研究。通过模拟发现，在2 ps的时候，颗粒内部开始发生相变现象，并且熔化的部分在颗粒内部所有区域均有分布，并没有出现一个明显的固液熔化界面。图2为激光照射时间t=5 ps时平均序参数随半径深度的分布情况。由图2可以看出，在任一颗粒半径深度均有平均序参数在[0，0.04]的原子，说明任一半径大小处均有熔化的原子，表明颗粒熔化发生在颗粒内部的所有区域，不存在明显的固液分界面，这和文献[13]中观察到的现象相一致。

图3为颗粒内部原子在不同时刻的微观投影图，红色表示固态原子，蓝色表示液态原子。由图3可以看出，在0 ps时刻，平衡状态的原子排列比较规则并遵守面心立方结构，原子在各自的原胞内作小范围的振动。随着激光能量的加入，颗粒的温度逐渐升高，整个颗粒中原子的晶格结构遭到破坏，在5 ps时刻，原子排列变得不规则，已经有部分原子发生熔化。之后随着激光的不断加入，在7.5 ps时，原子排列的更加无序，发生相变的原子更多。最后在12.5 ps，整个颗粒已经完全熔化。并且在熔化过程中，由于受到应力和压力的作用，颗粒已经有一定程度的变形。

图2 当t=5 ps时平均序参数随颗粒半径分布情况Fig.2 Distribution of average order parameter (t=5 ps) along with the change of particle radius

图3 不同时刻颗粒内部原子的微观投影图Fig.3 Microscopic projection of the gold nanoparticle at different moments

图4 为不同时刻颗粒内原子平均序参数的分布情况，图中虚线为固液分界线，虚线上半部分表示未发生相变的原子，虚线下半部分表示已经发生相变的原子。由图中可看出，在0 ps时所有原子都在固液分界线以上，随着激光的加入，颗粒内原子的平均序参数逐渐由固液分界线以上移至线以下，颗粒逐渐发生熔化，在12.5 ps时所有原子平均序参数都在固液分界线以下，颗粒完全熔化。

对烧结过程中的温度分布和熔化进程进行模拟分析，图5为电子和晶格温度随时间的变化情况。由图5可以看出，激光加入以后，电子温度Te在激光加入后大幅上升，在2.5 ps达到最高值，同时，由于从电子吸热，晶格温度Tl也随之上升，最后在5 ps附近稳定在900 K左右，10 ps附近电子和晶格温度达到平衡。图6为固态原子在颗粒所有原子中所占比例随时间的变化关系。由图6可以看出，颗粒开始发生熔化的时间为t=2 ps，而后随着激光的不断加入，发生熔化的颗粒逐渐增多，固态原子比例逐渐减少，在12 ps附近颗粒完全熔化。

图4 不同时刻颗粒内原子平均序参数分布情况Fig.4 Distribution of average order parameter at different moments

图5 晶格和电子温度随时间的变化Fig.5 Variations of lattice temperature Tl and electron temperature Te with time

图6 固态原子比例随时间的变化Fig.6 Variations of the fraction of solid atoms with time

2.3 不同激光能量大小的影响

研究激光参数对熔化过程的影响，保持颗粒半径和激光脉冲宽度不变，对激光强度分别为1.2 J/m2和0.6 J/m2的情况进行了模拟研究。图7为激光能量强度为1.2 J/m2和0.6 J/m2时纳米颗粒在不同时刻的原子微观投影图。同时对比图3可以看出，当激光能量增大时，在相同时刻颗粒内部熔化的比例较大，熔化速度更快，颗粒形状由于受到应力作用而发生的变化更明显，体积略微增大。当激光强度较小时，颗粒没有完全熔化，只有部分原子发生了熔化。当J=0.6 J/m2时，随着模拟过程的进行，当激光能量照射结束后，部分熔化的原子发生了再凝固现象，因而熔化的原子数又有所减少。

图7 当激光能量不同时颗粒内部原子不同时刻的微观投影图Fig.7 Microscopic projection of the gold nanoparticle at different moments and different intensities of laser energy

图8 和图9分别为不同激光能量下颗粒内部晶格的温度以及固态原子的比例随时间的变化曲线。当激光能量为0.6 J/m2时，晶格的平均温度稳定在820 K左右，虽然未达到颗粒的熔点835 K，但颗粒内部温度分布不均匀，已有部分原子达到熔点温度。因此，颗粒内部也有熔化现象的发生。随着激光能量的增大，晶格的平均温度也随之升高，当激光能量为1.2 J/m2时，晶格温度稳定在1 000 K左右，此时颗粒内部过热度已经很高。观察图8可以发现，当激光能量为1.2 J/m2和1.0 J/m2时，晶格温度达到最大值后会有下降的趋势。这是因为晶格从电子吸收能量后，温度急剧上升，上升至最大值时，熔化过程还未完全结束，而后随着内部原子的不平衡传热，部分能量转化为未熔化部分原子的熔化潜热，因而平均温度呈现一定的下降。当激光能量为1.2 J/m2时，熔化速度较快，在晶格温度上升至最大值时，未熔化原子的原子比例较少，相应的转化为熔化潜热的能量就低，所以，平均温度下降没有激光能量为1.0 J/m2时那么明显。从图9还可看到，激光能量较大时（J=1.0 J/m2和J=1.2 J/m2），颗粒最终完全熔化，并且激光能量越大，开始熔化的时刻越早，熔化的速度越快，总的熔化时间越短。当激光能量较小，等于0.6 J/m2时，颗粒没有完全熔化，最终固态原子的比例稳定在0.36左右。