脈沖激光燒蝕氬晶體顆粒的分子動力學模擬

梁小龍, 李 凌

(上海理工大學能源與動力工程學院,上海 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