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

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

(桂林電子科技大學機電工程學院, 桂林 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