高能量短脉冲激光作用水分子的动力学模拟

廖志强, 龙芋宏, 江 威, 童友群, 冯唐高

(桂林电子科技大学机电工程学院, 桂林 541004)

高能量短脉冲激光作用水分子的动力学模拟

廖志强, 龙芋宏, 江 威, 童友群, 冯唐高

(桂林电子科技大学机电工程学院, 桂林 541004)

基于高能量短脉冲激光作用水分子实验的复杂性,需用采取一种比较准确且较为方便的方法来预测实验结果是必要的.本文采用分子动力学方法(Molecular Dynamics)对高能量密度短脉冲激光加载下的超临界水分子进行热力学分析及结构研究.结果表明,随着能量快速加入加载区域,水分子系统的温度迅速提高.同时,由于分子迅速向四周扩散,在加载区域产生空泡.其中,只有25.5%的能量用以提高水分子系统的动能(温度),其余的能量都用于增大水分子系统的势能(水分子间距).伴随着温度的提高,水分子热运动加快,有序程度逐渐减弱,O—H间距增加,水分子间的氢键作用减弱,分子极性降低.分子动力学模拟方法研究短脉冲高能量激光对水分子的作用,对水下激光微加工研究有一定指导意义.

短脉冲; 高能量; 激光; SPC水分子; 分子动力学

1 引 言

水是地球上最常见的物质之一,是包括人类在内所有生命生存的重要资源,也是生物体最重要的组成部分.水同时是地球上最珍贵的资源之一,她能在地球表面的热力学条件下自然存在为汽、液、固三种物理形态,在很多学科及领域中,或多或少都涉及到水的问题[1].水是我们日常科学研究中用得最多的溶剂、冷却剂、介质之一.水在激光进行水下精密微加工时是一种良好的冷却剂,同时水也会与加工件产生相互影响.因此,研究水分子的微观结构对于生产研究是迫切需要的[2].在实验室中对水下金属进行短脉冲高能量激光加载实验是一项复杂而且耗时的工作.此类实验激光的脉宽是ns或者ps级,加工半径也是μm级甚至nm的,这些对于仪器在时间和空间上的精确性要求非常高,因而,这些会很大程度上的提高了实验的难度,同时需求昂贵测试手段.计算机模拟现已成为和理论分析及实验研究并列的三大研究手段之一,它普遍适用于解决多体问题.计算机模拟处理问题的优势在于：处理较为复杂的系统时,能减少繁重的实验工作,节省大量的人力、物力和财力[3].随着计算机技术的飞速发展,采用计算机模拟方法进行仿真实验和获得产品性能参数已经日益普及了[4].分子动力学模拟作为一种常见的材料结构分析手段而备受关注,它对材料的平衡、结构和动态性质方面有着良好的分析能力[5].前人在实验或仿真中发现,水在使其急速升温的高热流下发生爆发沸腾,在该区域形成空泡,并产生冲击波[6].本文采用分子动力学模拟手段对水在短脉冲高能量激光作用条件下的结构及热力学性质进行分析.

2 势能模型

首先本文要进行水分子模型的建立,在对水分子进行模型建立时要选择适合的势能参数.常见的势能模型有：MCY量子力学模型[7]和Simple Point Charge(SPC)、Transferable Intermolecular Potential Function(TIP2、TIP3、TIP4)、Central Force(CF)、ST2模型(由Stillinger F H提出)等.SPC模型的优势在其形式简单,利用SPC力场及分子动力学计算可以得到合理的水分子径向分布函数与各种动力学性质.所以,本文选择应用SPC模型,其相互作用势如下：

(1)

式(1)中,前半部分是静电相互作用,第二部分为短程Lennard-Jones作用,SPC模型假定只有氧原子间存在Lennard-Jones 作用.其中,u(ri,rj)是水分子之间的相互作用势能,rij是原子i与j之间的距离,qi和qj分别表示第i,j个原子上所带的电荷.ε,σ为氧原子之间的Lennard-Jones作用参数.SPC水分子的具体参数如表1.

对于库伦长程力的处理,我们选择PPPM(Particle-ParticleParticle-Mesh)算法,这是由于PPPM算法的速度比Ewald算法的速度快.而Ewald方法为各种计算静电作用的方法中准确性最高最为可行的方法.

表1 SPC水的势能参数

3 分子模拟过程

在进行水分子动力学模拟时,本文要同时考虑分子静电作用与水分子复杂的刚体运动,即同时考虑LJ势能和库伦力.本文模拟所用的水分子总数为2000个.在建立好分子模型后,我们首先使用NPT系综让系统逐渐平衡,使用NPT系统即意味着在平衡阶段保持系统中的原子数量N,系统压力P和系统温度T不变.为了保证运动方程求解的稳定性和求解的准确性,模拟中设置的步长为0.01fs.模拟时应用两体有效势能近似、最小镜像准则等基本假设,模拟过程中采用周期性边界条件以防止原子的丢失及控制模拟盒子的大小.势能截断采用球形截断法,截断半径为0.9nm.采用SHAKE算法[8]固定每个水分子的键长和键角,此种算法的好处是增加计算效率.牛顿运动方程积分方法采用经典的运动方程“leap-frog”,这种方法是在最早的Velert解法上改良而来,它的优点是不容易出现误差,其表达式为：

(2)

(3)

式中,r是各水分子的位置,v是各原子的速度,m是各原子的质量,Fi是各原子所受力在各方向上的分量,δt是计算时的积分步长.

系统经过1.2ps达到平衡态,初始温度为298K,初始压力为1atm.温度控制采用Nose-Hoover[9]温度控制技术.压力控制技术采用Parrinllo-Rahman[10]算法.系统到达平衡态后,使用NVE系综对其进行弛豫,时间为0.1ps,之后系统达到稳定状态.

在系统达到稳定状态后,进行非平衡模拟.在非平衡模拟阶段,在设定加载区域添加一个1ns的脉冲能量(由密集的小脉冲能量构成),其能量密度为1Kcal/mol,并撤销Nose-Hoover温度控制,只保留Parrinllo-Rahman压力控制,即NPH系综.因此,系统的大小可以随着压力的变化而自由膨胀或者收缩.

图 1 脉冲能量加载区域水分子随时间的位置图Fig.1 The distribution of water molecules in heating area

4 模拟结果及分析

4.1 热力学分析

在仿真计算完后,本文用VMD来实现水分子的可视化.图1为脉冲能量加载区域水分子随时间的位置图.从图1可以看出随着时间的前进,该区域内的分子数量急剧减少.由于在水分子模拟体系气泡形成阶段的模拟中只对模拟体系施加了压力控制,同时又有高能量进入加载区域,该区域的分子逐渐运动到其他区域,使得加载区域形成空泡,从而产生冲击波.陈笑等人[6]提出如果在水下金属表面加工时这种脉冲式的冲击波碰撞到金属表面会加剧金属表面的蚀除,并且产生反向的冲击回波,减弱下个脉冲产生的冲击波.

其后,本文对仿真获得的数据进行分析.分析其热力学特性及其径向分布函数.仿真的相关热力学数据如下：

图2和图3反映了水分子体系在给定的能量密度下1ns时间内的温度及气泡形成过程体系中水分子总能量、势能和动能随时间的变化的过程.如图2所示,在1Kcal/mol的能量密度下,1ns的能量加载会使水分子的温升率达到2×1011K/s,在这样高的温升率下,水分子会在短时间能气化,甚至产生离子化.如图3所示,约有25.5%的能量转变为了系统的动能,用于增加系统的温度.其余的能量转变为了系统的势能,而分子的势能增大使受热区域分子间距增大,转化为用来形成汽泡的潜热.这部分能量是无法用来提升系统温度的.

4.2 水分子结构分析

径向分布函数可以解释为系统的“区域密度”,径向分布函数反映出液体中分子聚集的特性,可藉此了解液体的“结构”[11].通过水分子的径向分布函数图可以分析高能短脉冲能量下水的微观结构特征.

图 2 加载中的温度曲线Fig.2 Temperature change during energy adding in

图 3 加载的能量曲线图Fig.3 Energy change during energy adding in

图4,5,6分别为水分子加载能量前后的O—O、O—H、H—H的径向分布函数.在能量加载之前,O—O径向分布函数的峰值分别在0.275nm和0.446nm出现.然而在能量加载之后,O—O径向分布函数的峰值分别为0.293nm和0.455nm,且峰值都稍微降低.这表明了温度迅速升高,分子热运动加快,水分子的有序程度逐渐减弱.而从O—H径向分布函数图可以得出第一和第二峰峰值降低,峰位稍微右移,而峰谷却略微提高.这与BruniF[12]等人之前做的中子衍射实验结果一致.这是由于水分子系统温度迅速升高,O—H间距增加,水分子间的氢键作用减弱,分子极性降低.

图 4 加载能量前后g(r)O—O对比Fig.4 g(r)O-O comparison

图 5 加载能量前后g(r)O—H对比Fig.5 g(r)O-H comparison

图 6 加载能量前后g(r)H—H对比Fig.6 g(r)H-H comparison

5 结 论

本文采用分子动力学方法对高能量密度短脉冲激光加载下的超临界水分子的热力学和结构进行了模拟分析.模拟结果表明系统温度迅速升高后,在加载区域产生空泡.其中,只有25.5%的能力用以提高水分子系统的动能,其余的能量都用于增大水分子系统的势能.伴随着温度的提高,水分子热运动加快,有序程度逐渐减弱,O—H间距增加,水分子间的氢键作用减弱,分子极性降低.本文采用分子动力学研究短脉冲高能量激光对水分子的作用,对水下激光微加工有一定指导作用.分子动力学模拟方法能较方便地分析物质微观领域的变化,能部分取代实验研究手段,有一定的实用性.下一步工作是进一步深入研究高能短脉冲激光作用水下金属的模型和分析.

[1]SunW,ChenZ,HuangS.Moleculardynamicssimulationofwaterambienttosupercriticalconditions[J].ComputersandAppliedChemistry, 2005, 22(05)： 406(in Chinese) [孙炜, 等.从常温到超临界条件下水的分子动力学模拟[J]. 计算机与应用化学, 2005, 22(05)： 406]

[2] Zhou J, Lu X, Wang Y,etal. Molecular dynamics simulation of water at different temperatures [J].ComputersandAppliedChemistry, 1999, 16(04)： 241(in Chinese) [周健, 等.不同温度下水的分子动力学模拟[J]. 计算机与应用化学, 1999, 16(04)： 241]

[3] Zhou Y, Huai X, Liang S Q. Molecular dynamics simulation of bubble nucleation in explosive boiling [J].JournalofEngineeringThermophysics, 2009, 39(06): 992(in Chinese) [邹玉, 等. 爆发沸腾形核分子动力学模拟[J]. 工程热物理学报, 2009, 39(06): 992]

[4] Allen M P, Tildesley D J.ComputerSimulationofLiquids[M]. Clarendon Press, Oxford, 1997: 230-234.

[5] Haile J M.MolecularDynamicsSimulation[M]. New York: Wiley, 1992: 145-146.

[6] Xu R, Chen X, Chen J,etal. Shock wave and cavitation effects by laser ablation of metal in water [J].ActcOpticaSinica, 2004, 24(12): 1643(in Chinese) [徐荣青, 等. 激光烧蚀水下金属产生冲击波和空泡效应的研究[J]. 光学学报, 2004, 24(12): 1643]

[7] Matsuoka O, Clementi E, Yoshimine M. CI study of the water dimer potential surface [J].Chem.Phys., 1976, 64: 1351.

[8] Andersen H C. Rattle: A “velocity” version of the shake algorithm for molecular dynamics calculations [J].JournalofComputationalPhysics, 1983, 52(1): 24.

[9] Hoover W G. Canonical dynamics: equilibrium phase-space distributions [J].PhysicalReviewA, 1985, 31(3): 1695.

[10] Parrinello M, Rahman A. Crystal structure and pair potentials: a molecular dynamics study [J].Phys.Rev.Lett., 1980, 14: 1196.

[11] Chen Z L.Theoryandpracticeofmolecularsimuletion[M]. Beijing: Chemical industry, 2007: 88-103(in Chinese)[陈正隆. 分子模拟的理论与实践 [M]. 北京:化学工业出版社, 2007: 88-103]

[12] Bruni F, Ricci M A, Soper A K. Unpredicted density dependence of hydrogen bonding in water found by neutron diffraction [J].Phys.Rev. B, 1996, 54: 11876.

Molecular dynamics simulation of water heated by short pulse and high-energy density laser

LIAO Zhi-Qiang, LONG Yu-Hong, JIANG Wei, TONG You-Qun, FENG Tang-Gao

(School of Mechanical and Electrical Engineering, GuiLin University of Electronic Technology, Guilin 541004, China)

Due to the complexity of the experiments for water heated by short pulse and high energy density laser, an accurate and simple method is needed to predict the result. Molecular Dynamics is applied to the thermodynamic analysis and structure research of water heated by short pulse and high energy density laser. As a result, the temperature of water system increases rapidly with energy added in. Meanwhile, bubbles generate due to the quick diffusion of water. Only 25.5 percent of the energy is used to increase the kinetic energy of water, the rest improves the potential energy. With the increasing of temperature, the thermal motion of water increases, the order degree of water reduces continue continuously, the distance of bonding O—H increases, the hydrogen-bond of water weakens and the molecular polarity reduces.

Short pulse; High energy; Laser; SPC water; Molecular dynamics

103969/j.issn.1000-0364.2015.02.015

2013-10-31

国家自然科学基金(51065007, 61366009)；广西自然科学基金(2012GXNSFAA053202)；广西制造系统与先进制造技术重点实验室主任项目(13-051-09-004Z)；桂林电子科技大学广西信息科学实验中心(20130313)

廖志强(1988—),男,广西人,硕士研究生,主要研究方向为激光加工机理研究.E-mail： fengshaxia21@gmail.com

O561

A

1000-0364(2015)02-0264-05