脉冲激光烧蚀氩晶体颗粒的分子动力学模拟

梁小龙, 李 凌

(上海理工大学能源与动力工程学院,上海 200093)

脉冲激光烧蚀氩晶体颗粒的分子动力学模拟

梁小龙, 李 凌

(上海理工大学能源与动力工程学院,上海 200093)

采用分子动力学方法对氩晶体颗粒在皮秒脉冲激光照射下内部发生的传热过程以及相变现象进行模拟研究.通过记录氩原子速度和原子位置随时间的变化情况分析了颗粒内热量的传递过程,并计算了单位体积内氩原子数目的空间分布情况以及随时间的变化规律,从而分析了颗粒内部的相变过程.研究结果表明:当激光强度较低时,晶体内部仅仅发生传热过程而没有熔化发生,加入的激光能量随时间由外向内传递;当激光强度增加到足够晶体熔化的时候会发生相变,此时固液相态之间并没有明显的界面,而是存在一个纳米级别的过渡区域;当照射的激光强度增加时,过渡区域的移动速度和移动深度都将增加.

分子动力学;脉冲激光;相变;温度分布

近年来由于激光技术的飞速发展,激光烧结技术在产品加工制造业内的应用日益广泛.激光烧结中多以粉末作为成型材料,激光束按照计算好的扫描路径对粉末床进行层层烧结成型,每层成型后的粉末层层叠加最终形成所要设计的三维产品.激光与材料相互作用的过程是一个十分复杂的物理过程,涉及很多学科领域,虽然已有大量文献对此领域进行研究,但是激光与物质相互作用的机制还依然不是十分清楚[1].目前国际上关于这方面的研究主要有两方面:一方面是实验研究;另一方面是数值模拟.激光烧结实验研究是随着激光技术的产生与发展展开的,这些实验主要能够获得材料的烧结形貌和烧结阈值等一些信息.而数值模拟则主要侧重于对烧结过程中传热机理的研究.传统的方法是连续介质力学,这种方法是一种宏观尺度的模拟方法,在处理有些问题时会遇到一些困难,如状态方程的确定、固液界面的产生、材料熔化的飞溅和高应变率下材料的选择等[2].而分子动力学模拟是一种在原子尺度上的确定性模拟方法,能够有效避免连续介质假设的局限.它将经典的牛顿运动力学应用到原子体系中,通过对牛顿方程积分来得到每个时刻原子的位置、速度和受力等微观信息,再用经典统计物理的相关理论来描述系统的宏观性质[3],这样就能避免连续介质假设的局限性.Herrmann等[4]首次使用分子动力学的方法模拟了超短波脉冲激光烧蚀硅的过程.Ohmura等[5]采用分子动力学的方法得到了脉冲激光的烧蚀深度和烧蚀的粒子数等.国内,王志军等[6]利用双温模型模拟了飞秒激光烧蚀金属镍的热影响区.白明泽等[7]采用耦合一维双温模型的分子动力学方法研究了纳米级的铝膜在飞秒激光辐照下的熔化机制.早期的分子动力学模拟技术受到计算机硬件水平的限制,其研究的时空尺度不够理想.近些年来,随着计算机技术和势函数理论的发展,分子动力学模拟技术得以推广.本文运用分子动力学的方法来研究激光烧蚀氩晶体颗粒过程,通过对每个氩原子进行受力分析,追踪其运动轨迹,以统计的方法[8]来模拟激光能量被吸收后颗粒内部的传热过程、发生相变时固液界面的产生与移动,以及激光强度对熔化过程的影响.

1 理论模型

1.1 物理问题

本文主要以单个固态氩颗粒为对象,应用分子动力学方法研究其在皮秒脉冲激光照射下的传热过程.考虑实际中激光照射到粉末床时,每一个颗粒表面都会发生反射现象,因而假设颗粒表面的热流密度均匀[9],如图1所示.颗粒的初始温度设为70 K,计算的氩原子个数为16 754.以球心为坐标原点建立坐标系,原子的初始位置由面心立方晶格结构确定.原子的初始速度根据初始温度在麦克斯韦玻尔兹曼速度分布中随机选取[10].

图1 物理模型Fig.1 Physical model

模拟过程中的时间步长为10 fs,初始时刻为5 ps,结束时刻为25 ps,其中包括了整个脉冲宽度.入射激光在空间中均匀分布,脉冲半高宽度为5 ps.激光脉冲随时间满足高斯分布,在5 ps时开始加载,15 ps时加载结束.

1.2 数学描述

分子之间势函数的正确确定是应用分子动力学方法的关键,论文采用普遍应用的较简单的Lennard-Jones 12-6势函数.在进行模拟的过程中,每一个原子都要受到它周围的每一个原子的Lennard-Jones作用力,其受到的合外力应该满足牛顿定理,即

式中,mi为原子i的质量;ri为原子i的位置;t为时间;Fij为原子i与原子j的相互作用力,其大小可由Lennard-Jones势求得,即Fij=-∂φij/∂rij.其中Lennard-Jones势φij为

式中,ε为势能阱的深度;σ是相互作用的原子i与原子j的势能为0时的距离;rij为原子i与原子j的距离.因此作用力Fij可以整理得到

在计算速度时采用Verlet算法[11],即

式中,vi为原子i的速度;δt为计算过程中的时间步长.

在计算过程中,为了减少计算时间,当两个原子的距离大于0.851 5 nm[12]时忽略其相互作用力.氩晶体颗粒对激光能量的吸收深度为2.5 nm[13],在此深度内分层进行.在激光照射过程中,激光能量的强度沿着颗粒表面向球心呈指数衰减[14],即

式中,I0为入射激光的初始光强;R为氩晶体表面反射率;α为激光在氩晶体中的吸收系数;z为距离晶体表面的距离.沉积在每一层上的能量被平均分配到该层的每一个原子上,原子吸收的激光能量全部转换为原子的动能,使得温度得以升高[15].

当没有发生相变时,非平衡状态下不同位置处的温度可由该位置所在区域一定厚度内的所有原子平均平动动能v—计算得到[16]

式中,kB是玻尔兹曼常数.

2 模拟计算及结果分析

首先对算法和程序的正确性进行验证.当系统达到平衡后,由式(6)得到颗粒温度为69.70 K.此时给颗粒加入能量8.02×10-18J,并当其达到新的平衡状态后,根据比热容c=q/(M·Δt)计算得到,固态氩在71.36 K时的比热容为0.935 kJ/(kg·K),这与文献[16]的0.895 kJ/(kg·K)基本吻合.

2.1 入射激光强度为0.01 J/m2

首先模拟了入射激光强度为0.01 J/m2时的情况,通过模拟发现颗粒内部没有相变发生.图2分别显示了氩颗粒在5,10,15,20和25 ps时内部温度的分布情况,l为距颗粒球心的距离.

从图2可以看出,在5 ps时,颗粒内外温度处于平衡状态,此时激光开始照射.到第10 ps时,激光脉冲达到最大值,在此过程中,靠近颗粒表面位置处首先升温,而靠近颗粒球心处温度基本维持不变,这时热量还没有传递到内部去.随着激光的继续照射,颗粒表面的温度进一步升高,颗粒内部由于分子之间的碰撞,热量逐渐被传递到更深处.到第15 ps时,激光照射停止,而后外部热源消失,表面温度逐渐降低,但由于内部温度的不平衡,热量继续向球心传递,靠近颗粒球心部分的温度继续升高.在第25 ps时,颗粒内部温度逐渐趋于一致,并且因为激光能量的加入,此时新的平衡温度将会高于初始的平衡温度.

图2 温度分布(I=0.01 J/m2)Fig.2 Temperature distribution(I=0.01 J/m2)

2.2 入射激光强度为0.1 J/m2

当激光强度增加为0.1 J/m2时,模拟结果发现颗粒温度在升高的同时表面局部区域发生了熔化现象.颗粒内原子在不同时刻的位置分布如图3所示.

图3 原子的位置分布(I=0.1 J/m2)Fig.3 Positions of atoms(I=0.1 J/m2)

从图3可以看到,在t=5 ps时,颗粒处于平衡状态,此时激光开始照射,颗粒吸收能量温度开始上升.在t=10 ps时,由于此时接受激光照射的时间只有5 ps,能量的传递时间还不够长,升温引起的形状变化不够明显,所以颗粒形状基本保持不变.而后随着激光的继续照射,颗粒的体积开始逐渐增大,靠近颗粒表面的原子密度有所降低,到t=15 ps时,激光停止照射.但由于之前吸收的热量使得靠近表面区域的温度还依然很高,热量继续向内部传递,并且颗粒的体积继续增大,靠近表面的原子分布更加稀疏.在t=20 ps和t=25 ps时,颗粒的形状进一步发生改变,靠近表面处的变化最为明显,能看到颗粒表面外部大量飞散的原子.由于同一种物质的相态不同时,其内部所含原子数的疏密程度也不同.为了对此现象进一步分析,计算出了不同位置处的原子数密度n,即单位体积的原子个数.图4为颗粒内部不同位置处的原子数密度及其随时间的变化情况.

图4 原子数密度分布(I=0.1 J/m2)Fig.4 Distribution of the number density of atoms(I=0.1 J/m2)

观察图4,在t=5 ps时,颗粒处于固态,其原子数密度在球内半径方向上均匀分布,此时加入激光照射.在t=10 ps时,靠近颗粒表面的位置处,原子数密度相比于平衡时有所降低,这是因为表面部分吸收了能量,温度有所上升,因而密度有所下降.在t=15 ps时,靠近球心处的原子数密度与初始状态基本一致,而靠近颗粒表面处原子数密度进一步降低,在距离半径5 nm左右的一段范围之内基本维持不变.经过计算,此段区域的平均原子数密度为1.04×1028个/m-3,小于固态氩的平均原子数密度,基本接近液态氩的原子数密度的最大值1.0× 1028个/m-3[16],因而可以判断此处发生了熔化.在t=20 ps和t=25 ps时,虽然激光照射已经停止,但表面温度依然很高,熔化继续向颗粒内部移动,并且还有一部分原子从表面处溢出形成气态的氩,因而使得靠近表面处的原子数密度继续减小.通过计算还发现当熔化发生时,在固态与液态之间没有明显的分界线,而是一段过渡区域,此过渡区域内,分子数密度沿径向逐渐减小.为了对相变发生时的固液区域进行描述,确定了一个判定标准,当原子数密度大于1.80×1028个/m-3时,为固态氩;当原子数密度小于1.0×1028个/m-3时,不考虑气态的情况下,为液态氩;而原子数密度在两个数值之间的为固液过渡区域[16-17].根据此判断标准,在t=5 ps时,原子数密度分布均在固态范围内.t= 10 ps时,在半径3.5 nm处开始出现固液过渡区域,但均未小于液态的判断标准,说明相变刚开始发生,还未形成液态.继续经过5 ps的激光照射,在t=15 ps时,靠近表面处有液态区域出现,固液过渡区域的范围为距离球心2.8～3.4 nm,固液过渡区域的前端相比于10 ps时,在激光的作用下向颗粒中心移动了0.7 nm.激光停止照射后到t= 20 ps时,固液过渡区域范围为2.6～3.5 nm.到t=25 ps时,固液过渡区域范围为2.8～3.4 nm,说明激光照射停止后的一段时间内过渡区域的移动也基本停止.

2.3 入射激光强度为0.12 J/m2

为了研究激光强度对固液过渡区域的影响,进一步增大激光强度到0.12 J/m2,不同位置处原子数密度随时间的变化情况如图5所示.

由图5可以看到,随着时间的推移,颗粒表面也发生了熔化现象,并且熔化部分也随着激光的照射向颗粒中心移动.与激光能量0.1 J/m2相比,当激光强度为0.12 J/m2时,在相同时刻、相同位置的原子数密度更小.这说明当激光强度升高时,能量的传递和温度的升高更加剧烈,使得相变过程更加迅速.同时,通过定量计算可知,t=10 ps时,距离半径3.2 nm处开始出现固液过渡区域.在t=15 ps时,固液过渡区域的范围为2.0～3.4 nm,说明在5 ps的激光作用下,固液过渡区域向颗粒内部移动了1.2 nm,与在能量密度为0.1 J/m2的激光照射下,相同的时间段内移动的0.7 nm相比,发现更强的激光能量密度使得固液过渡区域移动更为快速.激光停止照射后,在t= 20 ps时,固液过渡区域的范围为2.0～3.2 nm,而后直到在t=25 ps时,固液过渡区域位置基本没有变化.此阶段中固液过渡区域前端到达的最深处距离为距颗粒球心2.0 nm,比能量密度为0.1 J/m2的2.6 nm更深入,说明激光能量越强,熔化发生的区域越深.

图5 原子数密度分布(I=0.12 J/m2)Fig.5 Distribution of the number density of atoms(I=0.12 J/m2)

3 结 论

通过对氩颗粒在不同强度激光照射下传热过程的模拟研究,得到以下结论:

a.当激光能量比较弱时,氩晶体颗粒吸收激光能量而温度升高,但不会发生相变.靠近颗粒表面处的原子吸收能量,温度首先升高,而后随着时间的推移,通过原子之间的碰撞,能量被逐渐传递到颗粒更深的区域.

b.当激光能量足够强时,氩晶体颗粒会发生相变现象,此时固液相态之间并没有明显的界面,而是一段纳米级别的固液过渡区域,并且在激光的照射下向颗粒内部逐渐移动.

c.当相变发生时,激光能量越强,固液过渡区域向颗粒中心移动的速度越快,同时最终达到的位置也越深,即熔化深度也越大.

[1] 王丽梅,曾新吾.266 nm飞秒激光烧蚀单晶硅的分子动力学模拟[J].强激光与离子束,2008,20(8): 1360-1364.

[2] 李玉华,马法军,戴能利,等.超短脉冲激光对无机硅材料的损伤[J].中国激光,2007,34(7):1009-1013.

[3] 王馨,张江,吴金闪.统计物理学基础研究新进展[J].上海理工大学学报,2012,34(3):205-220.

[4] Herrmann R F W,Gerlach J,Campbell E E B. Molecular dynamics simulation of laser ablation of silicon[J].Nuclear Instruments and Methods in Physics Research Section B,1997,122(3):401-404.

[5] Ohmura E,Miyamoto I.Computer simulation on fusion and vaporization of metal due to laser irradiation[J]. Journal of High Temperature Society,1994,20(6): 228-335.

[6] 王志军,贾威,倪晓昌,等.飞秒激光烧蚀金属镍热影响区的数值模拟[J].激光技术,2007,31(6): 578-580.

[7] 白明泽,程丽,唐红,等.激光诱导铝纳米膜融化机理的分子动力学模拟[J].物理化学学报,2010,26(12): 3143-3149.

[8] 毕桥,方锦清,刘杰.独树一帜的统计物理理论——子动力学[J].上海理工大学学报,2012,34(2): 118-137.

[9] Zhang Y W,Chen J K.Ultrafast melting and resolidification of gold particle irradiated by pico-to femtosecond lasers[J].Journal of Applied Physics, 2008,104(5):054910.

[10] Allen M P,Tildesley D J.Computer simulation of liquid[M].London:Oxford University Press,1987.

[11] Swope W C,Anderon H C,Beren P H,et al.A computer simulation method for calculation of equilibrium constants for the formation of physical clusters of molecules:application to small water cluster[J].Journal of Chemical Physics,1982,76(1): 637-649.

[12] Rapaport D C.The art of molecular dynamics simulation[M].London:Cambridge University Press,2004.

[13] Wang X W,Xu X F.Molecular dynamics simulation of thermal and thermomechanical phenomena in picoseconds laser material interaction[J]. International Journal of Heat and Mass Transfer, 2003,46(1):45-53.

[14] 刘璇,王杨,赵丽杰.飞秒激光蚀除金属的分子动力学模拟[J].微细加工技术,2004(4):56-63.

[15] 辛建婷,祝文军,刘仓理.飞秒激光辐照铝材料的分子动力学数值模拟[J].爆炸与冲击,2004,24(3): 207-211.

[16] Wang X W.Thermal and thermomechanical phenomena in laser material interaction[D].Indiana: ETD Collection for Purdue University,2001.

[17] Wang X W.Thermal and thermomechanical phenomena in picosecond laser copper interaction[J]. Heat Transfer,2004,126(3):355-364.

(编辑:丁红艺)

Molecular Dynamics Simulation of Argon Crystal Particle Irradiated by Pulsed Laser

LIANGXiaolong, LILing
(School of Energy and Power Engineering,University of Shanghai for Science and Technology,Shanghai 200093,China)

The phenomena of heat transfer and phase change of solid argon crystal irradiated by pulsed laser was studied through the molecular dynamics simulation.In order to analyze the process of heat transfer inside the particle,the changes of velocities and locations of atoms were recorded at different time.Then the process of phases change was discussed through analyzing the variation of the number of argon molecules per unit volume.The calculation results show that when the laser intensity is low,no melting will happen.When the laser energy is increased enough,the phase change occurs and a transitional area of nanometer level instead of a clear interface is revealed between solid and liquid phase.The moving speed and moved depth of the transition zone increase with the increasing of laser intensity.

molecular dynamics;pulsed laser;phase change;temperature distribution

O 532+.25

A

1007-6735(2015)03-0245-06

10.13255/j.cnki.jusst.2015.03.008

2014-05-03

国家自然科学基金资助项目(51476102)

梁小龙(1989-),男,硕士研究生.研究方向:微尺度传热.E-mail:liangxiaolongeu@163.com

李 凌(1976-),女,副教授.研究方向:微尺度传热.E-mail:liling@usst.edu.cn