脈沖激光燒蝕氬晶體顆粒的分子動(dòng)力學(xué)模擬

梁小龍, 李 凌

(上海理工大學(xué)能源與動(dòng)力工程學(xué)院,上海 200093)

脈沖激光燒蝕氬晶體顆粒的分子動(dòng)力學(xué)模擬

梁小龍, 李 凌

(上海理工大學(xué)能源與動(dòng)力工程學(xué)院,上海 200093)

采用分子動(dòng)力學(xué)方法對氬晶體顆粒在皮秒脈沖激光照射下內(nèi)部發(fā)生的傳熱過程以及相變現(xiàn)象進(jìn)行模擬研究.通過記錄氬原子速度和原子位置隨時(shí)間的變化情況分析了顆粒內(nèi)熱量的傳遞過程,并計(jì)算了單位體積內(nèi)氬原子數(shù)目的空間分布情況以及隨時(shí)間的變化規(guī)律,從而分析了顆粒內(nèi)部的相變過程.研究結(jié)果表明:當(dāng)激光強(qiáng)度較低時(shí),晶體內(nèi)部僅僅發(fā)生傳熱過程而沒有熔化發(fā)生,加入的激光能量隨時(shí)間由外向內(nèi)傳遞;當(dāng)激光強(qiáng)度增加到足夠晶體熔化的時(shí)候會發(fā)生相變,此時(shí)固液相態(tài)之間并沒有明顯的界面,而是存在一個(gè)納米級別的過渡區(qū)域;當(dāng)照射的激光強(qiáng)度增加時(shí),過渡區(qū)域的移動(dòng)速度和移動(dòng)深度都將增加.

分子動(dòng)力學(xué);脈沖激光;相變;溫度分布

近年來由于激光技術(shù)的飛速發(fā)展,激光燒結(jié)技術(shù)在產(chǎn)品加工制造業(yè)內(nèi)的應(yīng)用日益廣泛.激光燒結(jié)中多以粉末作為成型材料,激光束按照計(jì)算好的掃描路徑對粉末床進(jìn)行層層燒結(jié)成型,每層成型后的粉末層層疊加最終形成所要設(shè)計(jì)的三維產(chǎn)品.激光與材料相互作用的過程是一個(gè)十分復(fù)雜的物理過程,涉及很多學(xué)科領(lǐng)域,雖然已有大量文獻(xiàn)對此領(lǐng)域進(jìn)行研究,但是激光與物質(zhì)相互作用的機(jī)制還依然不是十分清楚[1].目前國際上關(guān)于這方面的研究主要有兩方面:一方面是實(shí)驗(yàn)研究;另一方面是數(shù)值模擬.激光燒結(jié)實(shí)驗(yàn)研究是隨著激光技術(shù)的產(chǎn)生與發(fā)展展開的,這些實(shí)驗(yàn)主要能夠獲得材料的燒結(jié)形貌和燒結(jié)閾值等一些信息.而數(shù)值模擬則主要側(cè)重于對燒結(jié)過程中傳熱機(jī)理的研究.傳統(tǒng)的方法是連續(xù)介質(zhì)力學(xué),這種方法是一種宏觀尺度的模擬方法,在處理有些問題時(shí)會遇到一些困難,如狀態(tài)方程的確定、固液界面的產(chǎn)生、材料熔化的飛濺和高應(yīng)變率下材料的選擇等[2].而分子動(dòng)力學(xué)模擬是一種在原子尺度上的確定性模擬方法,能夠有效避免連續(xù)介質(zhì)假設(shè)的局限.它將經(jīng)典的牛頓運(yùn)動(dòng)力學(xué)應(yīng)用到原子體系中,通過對牛頓方程積分來得到每個(gè)時(shí)刻原子的位置、速度和受力等微觀信息,再用經(jīng)典統(tǒng)計(jì)物理的相關(guān)理論來描述系統(tǒng)的宏觀性質(zhì)[3],這樣就能避免連續(xù)介質(zhì)假設(shè)的局限性.Herrmann等[4]首次使用分子動(dòng)力學(xué)的方法模擬了超短波脈沖激光燒蝕硅的過程.Ohmura等[5]采用分子動(dòng)力學(xué)的方法得到了脈沖激光的燒蝕深度和燒蝕的粒子數(shù)等.國內(nèi),王志軍等[6]利用雙溫模型模擬了飛秒激光燒蝕金屬鎳的熱影響區(qū).白明澤等[7]采用耦合一維雙溫模型的分子動(dòng)力學(xué)方法研究了納米級的鋁膜在飛秒激光輻照下的熔化機(jī)制.早期的分子動(dòng)力學(xué)模擬技術(shù)受到計(jì)算機(jī)硬件水平的限制,其研究的時(shí)空尺度不夠理想.近些年來,隨著計(jì)算機(jī)技術(shù)和勢函數(shù)理論的發(fā)展,分子動(dòng)力學(xué)模擬技術(shù)得以推廣.本文運(yùn)用分子動(dòng)力學(xué)的方法來研究激光燒蝕氬晶體顆粒過程,通過對每個(gè)氬原子進(jìn)行受力分析,追蹤其運(yùn)動(dòng)軌跡,以統(tǒng)計(jì)的方法[8]來模擬激光能量被吸收后顆粒內(nèi)部的傳熱過程、發(fā)生相變時(shí)固液界面的產(chǎn)生與移動(dòng),以及激光強(qiáng)度對熔化過程的影響.

1 理論模型

1.1 物理問題

本文主要以單個(gè)固態(tài)氬顆粒為對象,應(yīng)用分子動(dòng)力學(xué)方法研究其在皮秒脈沖激光照射下的傳熱過程.考慮實(shí)際中激光照射到粉末床時(shí),每一個(gè)顆粒表面都會發(fā)生反射現(xiàn)象,因而假設(shè)顆粒表面的熱流密度均勻[9],如圖1所示.顆粒的初始溫度設(shè)為70 K,計(jì)算的氬原子個(gè)數(shù)為16 754.以球心為坐標(biāo)原點(diǎn)建立坐標(biāo)系,原子的初始位置由面心立方晶格結(jié)構(gòu)確定.原子的初始速度根據(jù)初始溫度在麥克斯韋玻爾茲曼速度分布中隨機(jī)選取[10].

圖1 物理模型Fig.1 Physical model

模擬過程中的時(shí)間步長為10 fs,初始時(shí)刻為5 ps,結(jié)束時(shí)刻為25 ps,其中包括了整個(gè)脈沖寬度.入射激光在空間中均勻分布,脈沖半高寬度為5 ps.激光脈沖隨時(shí)間滿足高斯分布,在5 ps時(shí)開始加載,15 ps時(shí)加載結(jié)束.

1.2 數(shù)學(xué)描述

分子之間勢函數(shù)的正確確定是應(yīng)用分子動(dòng)力學(xué)方法的關(guān)鍵,論文采用普遍應(yīng)用的較簡單的Lennard-Jones 12-6勢函數(shù).在進(jìn)行模擬的過程中,每一個(gè)原子都要受到它周圍的每一個(gè)原子的Lennard-Jones作用力,其受到的合外力應(yīng)該滿足牛頓定理,即

式中,mi為原子i的質(zhì)量;ri為原子i的位置;t為時(shí)間;Fij為原子i與原子j的相互作用力,其大小可由Lennard-Jones勢求得,即Fij=-?φij/?rij.其中Lennard-Jones勢φij為

式中,ε為勢能阱的深度;σ是相互作用的原子i與原子j的勢能為0時(shí)的距離;rij為原子i與原子j的距離.因此作用力Fij可以整理得到

在計(jì)算速度時(shí)采用Verlet算法[11],即

式中,vi為原子i的速度;δt為計(jì)算過程中的時(shí)間步長.

在計(jì)算過程中,為了減少計(jì)算時(shí)間,當(dāng)兩個(gè)原子的距離大于0.851 5 nm[12]時(shí)忽略其相互作用力.氬晶體顆粒對激光能量的吸收深度為2.5 nm[13],在此深度內(nèi)分層進(jìn)行.在激光照射過程中,激光能量的強(qiáng)度沿著顆粒表面向球心呈指數(shù)衰減[14],即

式中,I0為入射激光的初始光強(qiáng);R為氬晶體表面反射率;α為激光在氬晶體中的吸收系數(shù);z為距離晶體表面的距離.沉積在每一層上的能量被平均分配到該層的每一個(gè)原子上,原子吸收的激光能量全部轉(zhuǎn)換為原子的動(dòng)能,使得溫度得以升高[15].

當(dāng)沒有發(fā)生相變時(shí),非平衡狀態(tài)下不同位置處的溫度可由該位置所在區(qū)域一定厚度內(nèi)的所有原子平均平動(dòng)動(dòng)能v—計(jì)算得到[16]

式中,kB是玻爾茲曼常數(shù).

2 模擬計(jì)算及結(jié)果分析

首先對算法和程序的正確性進(jìn)行驗(yàn)證.當(dāng)系統(tǒng)達(dá)到平衡后,由式(6)得到顆粒溫度為69.70 K.此時(shí)給顆粒加入能量8.02×10-18J,并當(dāng)其達(dá)到新的平衡狀態(tài)后,根據(jù)比熱容c=q/(M·Δt)計(jì)算得到,固態(tài)氬在71.36 K時(shí)的比熱容為0.935 kJ/(kg·K),這與文獻(xiàn)[16]的0.895 kJ/(kg·K)基本吻合.

2.1 入射激光強(qiáng)度為0.01 J/m2

首先模擬了入射激光強(qiáng)度為0.01 J/m2時(shí)的情況,通過模擬發(fā)現(xiàn)顆粒內(nèi)部沒有相變發(fā)生.圖2分別顯示了氬顆粒在5,10,15,20和25 ps時(shí)內(nèi)部溫度的分布情況,l為距顆粒球心的距離.

從圖2可以看出,在5 ps時(shí),顆粒內(nèi)外溫度處于平衡狀態(tài),此時(shí)激光開始照射.到第10 ps時(shí),激光脈沖達(dá)到最大值,在此過程中,靠近顆粒表面位置處首先升溫,而靠近顆粒球心處溫度基本維持不變,這時(shí)熱量還沒有傳遞到內(nèi)部去.隨著激光的繼續(xù)照射,顆粒表面的溫度進(jìn)一步升高,顆粒內(nèi)部由于分子之間的碰撞,熱量逐漸被傳遞到更深處.到第15 ps時(shí),激光照射停止,而后外部熱源消失,表面溫度逐漸降低,但由于內(nèi)部溫度的不平衡,熱量繼續(xù)向球心傳遞,靠近顆粒球心部分的溫度繼續(xù)升高.在第25 ps時(shí),顆粒內(nèi)部溫度逐漸趨于一致,并且因?yàn)榧す饽芰康募尤?此時(shí)新的平衡溫度將會高于初始的平衡溫度.

圖2 溫度分布(I=0.01 J/m2)Fig.2 Temperature distribution(I=0.01 J/m2)

2.2 入射激光強(qiáng)度為0.1 J/m2

當(dāng)激光強(qiáng)度增加為0.1 J/m2時(shí),模擬結(jié)果發(fā)現(xiàn)顆粒溫度在升高的同時(shí)表面局部區(qū)域發(fā)生了熔化現(xiàn)象.顆粒內(nèi)原子在不同時(shí)刻的位置分布如圖3所示.

圖3 原子的位置分布(I=0.1 J/m2)Fig.3 Positions of atoms(I=0.1 J/m2)

從圖3可以看到,在t=5 ps時(shí),顆粒處于平衡狀態(tài),此時(shí)激光開始照射,顆粒吸收能量溫度開始上升.在t=10 ps時(shí),由于此時(shí)接受激光照射的時(shí)間只有5 ps,能量的傳遞時(shí)間還不夠長,升溫引起的形狀變化不夠明顯,所以顆粒形狀基本保持不變.而后隨著激光的繼續(xù)照射,顆粒的體積開始逐漸增大,靠近顆粒表面的原子密度有所降低,到t=15 ps時(shí),激光停止照射.但由于之前吸收的熱量使得靠近表面區(qū)域的溫度還依然很高,熱量繼續(xù)向內(nèi)部傳遞,并且顆粒的體積繼續(xù)增大,靠近表面的原子分布更加稀疏.在t=20 ps和t=25 ps時(shí),顆粒的形狀進(jìn)一步發(fā)生改變,靠近表面處的變化最為明顯,能看到顆粒表面外部大量飛散的原子.由于同一種物質(zhì)的相態(tài)不同時(shí),其內(nèi)部所含原子數(shù)的疏密程度也不同.為了對此現(xiàn)象進(jìn)一步分析,計(jì)算出了不同位置處的原子數(shù)密度n,即單位體積的原子個(gè)數(shù).圖4為顆粒內(nèi)部不同位置處的原子數(shù)密度及其隨時(shí)間的變化情況.

圖4 原子數(shù)密度分布(I=0.1 J/m2)Fig.4 Distribution of the number density of atoms(I=0.1 J/m2)

觀察圖4,在t=5 ps時(shí),顆粒處于固態(tài),其原子數(shù)密度在球內(nèi)半徑方向上均勻分布,此時(shí)加入激光照射.在t=10 ps時(shí),靠近顆粒表面的位置處,原子數(shù)密度相比于平衡時(shí)有所降低,這是因?yàn)楸砻娌糠治樟四芰?溫度有所上升,因而密度有所下降.在t=15 ps時(shí),靠近球心處的原子數(shù)密度與初始狀態(tài)基本一致,而靠近顆粒表面處原子數(shù)密度進(jìn)一步降低,在距離半徑5 nm左右的一段范圍之內(nèi)基本維持不變.經(jīng)過計(jì)算,此段區(qū)域的平均原子數(shù)密度為1.04×1028個(gè)/m-3,小于固態(tài)氬的平均原子數(shù)密度,基本接近液態(tài)氬的原子數(shù)密度的最大值1.0× 1028個(gè)/m-3[16],因而可以判斷此處發(fā)生了熔化.在t=20 ps和t=25 ps時(shí),雖然激光照射已經(jīng)停止,但表面溫度依然很高,熔化繼續(xù)向顆粒內(nèi)部移動(dòng),并且還有一部分原子從表面處溢出形成氣態(tài)的氬,因而使得靠近表面處的原子數(shù)密度繼續(xù)減小.通過計(jì)算還發(fā)現(xiàn)當(dāng)熔化發(fā)生時(shí),在固態(tài)與液態(tài)之間沒有明顯的分界線,而是一段過渡區(qū)域,此過渡區(qū)域內(nèi),分子數(shù)密度沿徑向逐漸減小.為了對相變發(fā)生時(shí)的固液區(qū)域進(jìn)行描述,確定了一個(gè)判定標(biāo)準(zhǔn),當(dāng)原子數(shù)密度大于1.80×1028個(gè)/m-3時(shí),為固態(tài)氬;當(dāng)原子數(shù)密度小于1.0×1028個(gè)/m-3時(shí),不考慮氣態(tài)的情況下,為液態(tài)氬;而原子數(shù)密度在兩個(gè)數(shù)值之間的為固液過渡區(qū)域[16-17].根據(jù)此判斷標(biāo)準(zhǔn),在t=5 ps時(shí),原子數(shù)密度分布均在固態(tài)范圍內(nèi).t= 10 ps時(shí),在半徑3.5 nm處開始出現(xiàn)固液過渡區(qū)域,但均未小于液態(tài)的判斷標(biāo)準(zhǔn),說明相變剛開始發(fā)生,還未形成液態(tài).繼續(xù)經(jīng)過5 ps的激光照射,在t=15 ps時(shí),靠近表面處有液態(tài)區(qū)域出現(xiàn),固液過渡區(qū)域的范圍為距離球心2.8～3.4 nm,固液過渡區(qū)域的前端相比于10 ps時(shí),在激光的作用下向顆粒中心移動(dòng)了0.7 nm.激光停止照射后到t= 20 ps時(shí),固液過渡區(qū)域范圍為2.6～3.5 nm.到t=25 ps時(shí),固液過渡區(qū)域范圍為2.8～3.4 nm,說明激光照射停止后的一段時(shí)間內(nèi)過渡區(qū)域的移動(dòng)也基本停止.

2.3 入射激光強(qiáng)度為0.12 J/m2

為了研究激光強(qiáng)度對固液過渡區(qū)域的影響,進(jìn)一步增大激光強(qiáng)度到0.12 J/m2,不同位置處原子數(shù)密度隨時(shí)間的變化情況如圖5所示.

由圖5可以看到,隨著時(shí)間的推移,顆粒表面也發(fā)生了熔化現(xiàn)象,并且熔化部分也隨著激光的照射向顆粒中心移動(dòng).與激光能量0.1 J/m2相比,當(dāng)激光強(qiáng)度為0.12 J/m2時(shí),在相同時(shí)刻、相同位置的原子數(shù)密度更小.這說明當(dāng)激光強(qiáng)度升高時(shí),能量的傳遞和溫度的升高更加劇烈,使得相變過程更加迅速.同時(shí),通過定量計(jì)算可知,t=10 ps時(shí),距離半徑3.2 nm處開始出現(xiàn)固液過渡區(qū)域.在t=15 ps時(shí),固液過渡區(qū)域的范圍為2.0～3.4 nm,說明在5 ps的激光作用下,固液過渡區(qū)域向顆粒內(nèi)部移動(dòng)了1.2 nm,與在能量密度為0.1 J/m2的激光照射下,相同的時(shí)間段內(nèi)移動(dòng)的0.7 nm相比,發(fā)現(xiàn)更強(qiáng)的激光能量密度使得固液過渡區(qū)域移動(dòng)更為快速.激光停止照射后,在t= 20 ps時(shí),固液過渡區(qū)域的范圍為2.0～3.2 nm,而后直到在t=25 ps時(shí),固液過渡區(qū)域位置基本沒有變化.此階段中固液過渡區(qū)域前端到達(dá)的最深處距離為距顆粒球心2.0 nm,比能量密度為0.1 J/m2的2.6 nm更深入,說明激光能量越強(qiáng),熔化發(fā)生的區(qū)域越深.

圖5 原子數(shù)密度分布(I=0.12 J/m2)Fig.5 Distribution of the number density of atoms(I=0.12 J/m2)

3 結(jié) 論

通過對氬顆粒在不同強(qiáng)度激光照射下傳熱過程的模擬研究,得到以下結(jié)論:

a.當(dāng)激光能量比較弱時(shí),氬晶體顆粒吸收激光能量而溫度升高,但不會發(fā)生相變.靠近顆粒表面處的原子吸收能量,溫度首先升高,而后隨著時(shí)間的推移,通過原子之間的碰撞,能量被逐漸傳遞到顆粒更深的區(qū)域.

b.當(dāng)激光能量足夠強(qiáng)時(shí),氬晶體顆粒會發(fā)生相變現(xiàn)象,此時(shí)固液相態(tài)之間并沒有明顯的界面,而是一段納米級別的固液過渡區(qū)域,并且在激光的照射下向顆粒內(nèi)部逐漸移動(dòng).

c.當(dāng)相變發(fā)生時(shí),激光能量越強(qiáng),固液過渡區(qū)域向顆粒中心移動(dòng)的速度越快,同時(shí)最終達(dá)到的位置也越深,即熔化深度也越大.

[1] 王麗梅,曾新吾.266 nm飛秒激光燒蝕單晶硅的分子動(dòng)力學(xué)模擬[J].強(qiáng)激光與離子束,2008,20(8): 1360-1364.

[2] 李玉華,馬法軍,戴能利,等.超短脈沖激光對無機(jī)硅材料的損傷[J].中國激光,2007,34(7):1009-1013.

[3] 王馨,張江,吳金閃.統(tǒng)計(jì)物理學(xué)基礎(chǔ)研究新進(jìn)展[J].上海理工大學(xué)學(xué)報(bào),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] 王志軍,賈威,倪曉昌,等.飛秒激光燒蝕金屬鎳熱影響區(qū)的數(shù)值模擬[J].激光技術(shù),2007,31(6): 578-580.

[7] 白明澤,程麗,唐紅,等.激光誘導(dǎo)鋁納米膜融化機(jī)理的分子動(dòng)力學(xué)模擬[J].物理化學(xué)學(xué)報(bào),2010,26(12): 3143-3149.

[8] 畢橋,方錦清,劉杰.獨(dú)樹一幟的統(tǒng)計(jì)物理理論——子動(dòng)力學(xué)[J].上海理工大學(xué)學(xué)報(bào),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] 劉璇,王楊,趙麗杰.飛秒激光蝕除金屬的分子動(dòng)力學(xué)模擬[J].微細(xì)加工技術(shù),2004(4):56-63.

[15] 辛建婷,祝文軍,劉倉理.飛秒激光輻照鋁材料的分子動(dòng)力學(xué)數(shù)值模擬[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

國家自然科學(xué)基金資助項(xiàng)目(51476102)

梁小龍(1989-),男,碩士研究生.研究方向:微尺度傳熱.E-mail:liangxiaolongeu@163.com

李 凌(1976-),女,副教授.研究方向:微尺度傳熱.E-mail:liling@usst.edu.cn