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

廖志強(qiáng), 龍芋宏, 江 威, 童友群, 馮唐高

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

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

廖志強(qiáng), 龍芋宏, 江 威, 童友群, 馮唐高

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

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

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

1 引 言

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

2 勢(shì)能模型

首先本文要進(jìn)行水分子模型的建立,在對(duì)水分子進(jìn)行模型建立時(shí)要選擇適合的勢(shì)能參數(shù).常見(jiàn)的勢(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)勢(shì)在其形式簡(jiǎn)單,利用SPC力場(chǎng)及分子動(dòng)力學(xué)計(jì)算可以得到合理的水分子徑向分布函數(shù)與各種動(dòng)力學(xué)性質(zhì).所以,本文選擇應(yīng)用SPC模型,其相互作用勢(shì)如下：

(1)

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

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

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

3 分子模擬過(guò)程

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

(2)

(3)

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

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

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

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

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

4.1 熱力學(xué)分析

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

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

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

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

徑向分布函數(shù)可以解釋為系統(tǒng)的“區(qū)域密度”,徑向分布函數(shù)反映出液體中分子聚集的特性,可藉此了解液體的“結(jié)構(gòu)”[11].通過(guò)水分子的徑向分布函數(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,且峰值都稍微降低.這表明了溫度迅速升高,分子熱運(yùn)動(dòng)加快,水分子的有序程度逐漸減弱.而從O—H徑向分布函數(shù)圖可以得出第一和第二峰峰值降低,峰位稍微右移,而峰谷卻略微提高.這與BruniF[12]等人之前做的中子衍射實(shí)驗(yàn)結(jié)果一致.這是由于水分子系統(tǒng)溫度迅速升高,O—H間距增加,水分子間的氫鍵作用減弱,分子極性降低.

圖 4 加載能量前后g(r)O—O對(duì)比Fig.4 g(r)O-O comparison

圖 5 加載能量前后g(r)O—H對(duì)比Fig.5 g(r)O-H comparison

圖 6 加載能量前后g(r)H—H對(duì)比Fig.6 g(r)H-H comparison

5 結(jié) 論

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

[1]SunW,ChenZ,HuangS.Moleculardynamicssimulationofwaterambienttosupercriticalconditions[J].ComputersandAppliedChemistry, 2005, 22(05)： 406(in Chinese) [孫煒, 等.從常溫到超臨界條件下水的分子動(dòng)力學(xué)模擬[J]. 計(jì)算機(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) [周健, 等.不同溫度下水的分子動(dòng)力學(xué)模擬[J]. 計(jì)算機(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ā)沸騰形核分子動(dòng)力學(xué)模擬[J]. 工程熱物理學(xué)報(bào), 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é)報(bào), 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)[陳正隆. 分子模擬的理論與實(shí)踐 [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

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

廖志強(qiáng)(1988—),男,廣西人,碩士研究生,主要研究方向?yàn)榧す饧庸C(jī)理研究.E-mail： fengshaxia21@gmail.com

O561

A

1000-0364(2015)02-0264-05