廖志強, 龍芋宏, 江 威, 童友群, 馮唐高

(桂林電子科技大學(xué)機電工程學(xué)院, 桂林 541004)

高能量短脈沖激光作用水分子的動力學(xué)模擬

廖志強, 龍芋宏, 江 威, 童友群, 馮唐高

(桂林電子科技大學(xué)機電工程學(xué)院, 桂林 541004)

基于高能量短脈沖激光作用水分子實驗的復(fù)雜性,需用采取一種比較準(zhǔn)確且較為方便的方法來預(yù)測實驗結(jié)果是必要的.本文采用分子動力學(xué)方法(Molecular Dynamics)對高能量密度短脈沖激光加載下的超臨界水分子進行熱力學(xué)分析及結(jié)構(gòu)研究.結(jié)果表明,隨著能量快速加入加載區(qū)域,水分子系統(tǒng)的溫度迅速提高.同時,由于分子迅速向四周擴散,在加載區(qū)域產(chǎn)生空泡.其中,只有25.5%的能量用以提高水分子系統(tǒng)的動能(溫度),其余的能量都用于增大水分子系統(tǒng)的勢能(水分子間距).伴隨著溫度的提高,水分子熱運動加快,有序程度逐漸減弱,O—H間距增加,水分子間的氫鍵作用減弱,分子極性降低.分子動力學(xué)模擬方法研究短脈沖高能量激光對水分子的作用,對水下激光微加工研究有一定指導(dǎo)意義.

短脈沖; 高能量; 激光; SPC水分子; 分子動力學(xué)

1 引 言

水是地球上最常見的物質(zhì)之一,是包括人類在內(nèi)所有生命生存的重要資源,也是生物體最重要的組成部分.水同時是地球上最珍貴的資源之一,她能在地球表面的熱力學(xué)條件下自然存在為汽、液、固三種物理形態(tài),在很多學(xué)科及領(lǐng)域中,或多或少都涉及到水的問題[1].水是我們?nèi)粘？茖W(xué)研究中用得最多的溶劑、冷卻劑、介質(zhì)之一.水在激光進行水下精密微加工時是一種良好的冷卻劑,同時水也會與加工件產(chǎn)生相互影響.因此,研究水分子的微觀結(jié)構(gòu)對于生產(chǎn)研究是迫切需要的[2].在實驗室中對水下金屬進行短脈沖高能量激光加載實驗是一項復(fù)雜而且耗時的工作.此類實驗激光的脈寬是ns或者ps級,加工半徑也是μm級甚至nm的,這些對于儀器在時間和空間上的精確性要求非常高,因而,這些會很大程度上的提高了實驗的難度,同時需求昂貴測試手段.計算機模擬現(xiàn)已成為和理論分析及實驗研究并列的三大研究手段之一,它普遍適用于解決多體問題.計算機模擬處理問題的優(yōu)勢在于：處理較為復(fù)雜的系統(tǒng)時,能減少繁重的實驗工作,節(jié)省大量的人力、物力和財力[3].隨著計算機技術(shù)的飛速發(fā)展,采用計算機模擬方法進行仿真實驗和獲得產(chǎn)品性能參數(shù)已經(jīng)日益普及了[4].分子動力學(xué)模擬作為一種常見的材料結(jié)構(gòu)分析手段而備受關(guān)注,它對材料的平衡、結(jié)構(gòu)和動態(tài)性質(zhì)方面有著良好的分析能力[5].前人在實驗或仿真中發(fā)現(xiàn),水在使其急速升溫的高熱流下發(fā)生爆發(fā)沸騰,在該區(qū)域形成空泡,并產(chǎn)生沖擊波[6].本文采用分子動力學(xué)模擬手段對水在短脈沖高能量激光作用條件下的結(jié)構(gòu)及熱力學(xué)性質(zhì)進行分析.

2 勢能模型

首先本文要進行水分子模型的建立,在對水分子進行模型建立時要選擇適合的勢能參數(shù).常見的勢能模型有：MCY量子力學(xué)模型[7]和Simple Point Charge(SPC)、Transferable Intermolecular Potential Function(TIP2、TIP3、TIP4)、Central Force(CF)、ST2模型(由Stillinger F H提出)等.SPC模型的優(yōu)勢在其形式簡單,利用SPC力場及分子動力學(xué)計算可以得到合理的水分子徑向分布函數(shù)與各種動力學(xué)性質(zhì).所以,本文選擇應(yīng)用SPC模型,其相互作用勢如下：

(1)

式(1)中,前半部分是靜電相互作用,第二部分為短程Lennard-Jones作用,SPC模型假定只有氧原子間存在Lennard-Jones 作用.其中,u(ri,rj)是水分子之間的相互作用勢能,rij是原子i與j之間的距離,qi和qj分別表示第i,j個原子上所帶的電荷.ε,σ為氧原子之間的Lennard-Jones作用參數(shù).SPC水分子的具體參數(shù)如表1.

對于庫倫長程力的處理,我們選擇PPPM(Particle-ParticleParticle-Mesh)算法,這是由于PPPM算法的速度比Ewald算法的速度快.而Ewald方法為各種計算靜電作用的方法中準(zhǔn)確性最高最為可行的方法.

表1 SPC水的勢能參數(shù)

3 分子模擬過程

在進行水分子動力學(xué)模擬時,本文要同時考慮分子靜電作用與水分子復(fù)雜的剛體運動,即同時考慮LJ勢能和庫倫力.本文模擬所用的水分子總數(shù)為2000個.在建立好分子模型后,我們首先使用NPT系綜讓系統(tǒng)逐漸平衡,使用NPT系統(tǒng)即意味著在平衡階段保持系統(tǒng)中的原子數(shù)量N,系統(tǒng)壓力P和系統(tǒng)溫度T不變.為了保證運動方程求解的穩(wěn)定性和求解的準(zhǔn)確性,模擬中設(shè)置的步長為0.01fs.模擬時應(yīng)用兩體有效勢能近似、最小鏡像準(zhǔn)則等基本假設(shè),模擬過程中采用周期性邊界條件以防止原子的丟失及控制模擬盒子的大小.勢能截斷采用球形截斷法,截斷半徑為0.9nm.采用SHAKE算法[8]固定每個水分子的鍵長和鍵角,此種算法的好處是增加計算效率.牛頓運動方程積分方法采用經(jīng)典的運動方程“l(fā)eap-frog”,這種方法是在最早的Velert解法上改良而來,它的優(yōu)點是不容易出現(xiàn)誤差,其表達式為：

(2)

(3)

式中,r是各水分子的位置,v是各原子的速度,m是各原子的質(zhì)量,Fi是各原子所受力在各方向上的分量,δt是計算時的積分步長.

系統(tǒng)經(jīng)過1.2ps達到平衡態(tài),初始溫度為298K,初始壓力為1atm.溫度控制采用Nose-Hoover[9]溫度控制技術(shù).壓力控制技術(shù)采用Parrinllo-Rahman[10]算法.系統(tǒng)到達平衡態(tài)后,使用NVE系綜對其進行弛豫,時間為0.1ps,之后系統(tǒng)達到穩(wěn)定狀態(tài).

在系統(tǒng)達到穩(wěn)定狀態(tài)后,進行非平衡模擬.在非平衡模擬階段,在設(shè)定加載區(qū)域添加一個1ns的脈沖能量(由密集的小脈沖能量構(gòu)成),其能量密度為1Kcal/mol,并撤銷Nose-Hoover溫度控制,只保留Parrinllo-Rahman壓力控制,即NPH系綜.因此,系統(tǒng)的大小可以隨著壓力的變化而自由膨脹或者收縮.

圖 1 脈沖能量加載區(qū)域水分子隨時間的位置圖Fig.1 The distribution of water molecules in heating area

4 模擬結(jié)果及分析

4.1 熱力學(xué)分析

在仿真計算完后,本文用VMD來實現(xiàn)水分子的可視化.圖1為脈沖能量加載區(qū)域水分子隨時間的位置圖.從圖1可以看出隨著時間的前進,該區(qū)域內(nèi)的分子數(shù)量急劇減少.由于在水分子模擬體系氣泡形成階段的模擬中只對模擬體系施加了壓力控制,同時又有高能量進入加載區(qū)域,該區(qū)域的分子逐漸運動到其他區(qū)域,使得加載區(qū)域形成空泡,從而產(chǎn)生沖擊波.陳笑等人[6]提出如果在水下金屬表面加工時這種脈沖式的沖擊波碰撞到金屬表面會加劇金屬表面的蝕除,并且產(chǎn)生反向的沖擊回波,減弱下個脈沖產(chǎn)生的沖擊波.

其后,本文對仿真獲得的數(shù)據(jù)進行分析.分析其熱力學(xué)特性及其徑向分布函數(shù).仿真的相關(guān)熱力學(xué)數(shù)據(jù)如下：

圖2和圖3反映了水分子體系在給定的能量密度下1ns時間內(nèi)的溫度及氣泡形成過程體系中水分子總能量、勢能和動能隨時間的變化的過程.如圖2所示,在1Kcal/mol的能量密度下,1ns的能量加載會使水分子的溫升率達到2×1011K/s,在這樣高的溫升率下,水分子會在短時間能氣化,甚至產(chǎn)生離子化.如圖3所示,約有25.5%的能量轉(zhuǎn)變?yōu)榱讼到y(tǒng)的動能,用于增加系統(tǒng)的溫度.其余的能量轉(zhuǎn)變?yōu)榱讼到y(tǒng)的勢能,而分子的勢能增大使受熱區(qū)域分子間距增大,轉(zhuǎn)化為用來形成汽泡的潛熱.這部分能量是無法用來提升系統(tǒng)溫度的.

4.2 水分子結(jié)構(gòu)分析

徑向分布函數(shù)可以解釋為系統(tǒng)的“區(qū)域密度”,徑向分布函數(shù)反映出液體中分子聚集的特性,可藉此了解液體的“結(jié)構(gòu)”[11].通過水分子的徑向分布函數(shù)圖可以分析高能短脈沖能量下水的微觀結(jié)構(gòu)特征.

圖 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的徑向分布函數(shù).在能量加載之前,O—O徑向分布函數(shù)的峰值分別在0.275nm和0.446nm出現(xiàn).然而在能量加載之后,O—O徑向分布函數(shù)的峰值分別為0.293nm和0.455nm,且峰值都稍微降低.這表明了溫度迅速升高,分子熱運動加快,水分子的有序程度逐漸減弱.而從O—H徑向分布函數(shù)圖可以得出第一和第二峰峰值降低,峰位稍微右移,而峰谷卻略微提高.這與BruniF[12]等人之前做的中子衍射實驗結(jié)果一致.這是由于水分子系統(tǒng)溫度迅速升高,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 結(jié) 論

本文采用分子動力學(xué)方法對高能量密度短脈沖激光加載下的超臨界水分子的熱力學(xué)和結(jié)構(gòu)進行了模擬分析.模擬結(jié)果表明系統(tǒng)溫度迅速升高后,在加載區(qū)域產(chǎn)生空泡.其中,只有25.5%的能力用以提高水分子系統(tǒng)的動能,其余的能量都用于增大水分子系統(tǒng)的勢能.伴隨著溫度的提高,水分子熱運動加快,有序程度逐漸減弱,O—H間距增加,水分子間的氫鍵作用減弱,分子極性降低.本文采用分子動力學(xué)研究短脈沖高能量激光對水分子的作用,對水下激光微加工有一定指導(dǎo)作用.分子動力學(xué)模擬方法能較方便地分析物質(zhì)微觀領(lǐng)域的變化,能部分取代實驗研究手段,有一定的實用性.下一步工作是進一步深入研究高能短脈沖激光作用水下金屬的模型和分析.

[1]SunW,ChenZ,HuangS.Moleculardynamicssimulationofwaterambienttosupercriticalconditions[J].ComputersandAppliedChemistry, 2005, 22(05)： 406(in Chinese) [孫煒, 等.從常溫到超臨界條件下水的分子動力學(xué)模擬[J]. 計算機與應(yīng)用化學(xué), 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) [周健, 等.不同溫度下水的分子動力學(xué)模擬[J]. 計算機與應(yīng)用化學(xué), 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) [鄒玉, 等. 爆發(fā)沸騰形核分子動力學(xué)模擬[J]. 工程熱物理學(xué)報, 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) [徐榮青, 等. 激光燒蝕水下金屬產(chǎn)生沖擊波和空泡效應(yīng)的研究[J]. 光學(xué)學(xué)報, 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]. 北京:化學(xué)工業(yè)出版社, 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

國家自然科學(xué)基金(51065007, 61366009)；廣西自然科學(xué)基金(2012GXNSFAA053202)；廣西制造系統(tǒng)與先進制造技術(shù)重點實驗室主任項目(13-051-09-004Z)；桂林電子科技大學(xué)廣西信息科學(xué)實驗中心(20130313)

廖志強(1988—),男,廣西人,碩士研究生,主要研究方向為激光加工機理研究.E-mail： fengshaxia21@gmail.com

O561

A

1000-0364(2015)02-0264-05