考慮非傅立葉效應的固液相變分子動力學模擬

關 陽， 李 凌， 牛澤偉

(上海理工大學，能源與動力工程學院，上海 200093)

1 引 言

飛秒激光自開始用于材料加工以來，憑借其獨特的優勢成為微細加工的理想工具，并被廣泛應用于不同的領域[1-4].飛秒激光與材料相互作用的中間過程對結果影響很大，尤其是加工的材料越來越小，加工時間也越來越短，所以研究加工過程中的傳熱機理有著重要意義. 目前對飛秒激光和金屬材料的相互作用機理已有大量的研究，很多都是采用連續介質模型，使用最廣泛的是基于傅立葉定律的拋物線兩步加熱模型[5-7].連續介質模型的缺點是無法獲得激光與材料相互作用過程的微觀特性.分子動力學方法是一種原子尺度的數值模擬方法，可以從微觀角度對材料的相變過程進行直觀描述.但經典的分子動力學方法不能描述激光與材料作用時電子晶格間的非平衡傳熱過程.并且在處理飛秒激光加工材料的過程中，由于激光脈寬接近或小于電子和晶格的弛豫時間，需要考慮非傅里葉效應，尤其是發生快速的相變過程時.

因此，為了兼顧連續介質模型與分子動力學方法的優點，本文分別將基于傅立葉定律的拋物線模型(parabolic model)和考慮非傅立葉效應的雙曲線模型(hyperbolic model)與分子動力學方法相結合來研究飛秒激光照射金箔過程中發生的固液相變現象，并對得到的結果進行對比分析，從而獲取了固液界面的變化規律，隨后采用考慮非傅立葉效應的模型探討了激光能量密度和脈沖寬度對相變過程的影響.

2 物理模型和數學方法

2.1 物理模型

激光照射金箔的物理模型如圖1所示.脈寬tp=200 fs、能量密度J=50 J/m2的激光從上邊均勻照射金箔，激光的作用時間為800 fs.金箔厚度L=12.24 nm，初始溫度為300 K，電子和晶格的弛豫時間分別是40 fs和800 fs[8].根據文獻[9]，當模擬的原子數超過500時，即具有統計學意義，本文計算的模型共包含1089個原子.金的勢函數采用嵌入原子(EAM)勢函數[10].長(x)、寬(y)方向均采用周期性邊界條件，金箔下表面設為彈性表面[11]，上表面設為自由表面.

圖1 物理模型Fig. 1 Physical model

2.2 耦合雙溫度模型的分子動力學模擬

基于傅立葉定律的拋物線模型下電子等效熱傳導方程為[12]：

(1)

非傅立葉效應下考慮電子和晶格弛豫效應的影響，雙曲線模型中電子等效熱傳導方程為[8]：

(2)

式中Ce為電子熱容[13]，t為時間，T為溫度，τ為弛豫時間，下腳標e和l分別代表電子和晶格，ke為等效電子導熱系數[14]，G為電子-晶格的耦合因子[15]，S為源項[12].

由于金箔中電子的熱熔比晶格的熱熔小，當激光作用金箔時，電子溫度在短時間內急劇升高，而此時晶格溫度的變化很小，所以忽略激光與晶格的直接作用.與此同時，電子聲子的耦合作用使金箔內原子產生了附加運動，隨著電子溫度的不斷升高原子運動逐漸加劇，隨即通過統計原子的運動狀態從微觀角度描述晶格系統.

加入激光后原子的受力方程為[16]：

(3)

式中Fi為第i個原子在激光作用前所受的力，mi為第i個原子的質量，ri(x，y，z，t)為某時刻原子所在位置，t為時間，viT為原子的運動速度vi與該原子所在層的運動速度vc的差值.ξ為修正系數，其計算公式為[17]：

(4)

n=ΔtMD/ΔtTTM，ΔtMD為分子動力學模型中的時間步長，ΔtTTM為雙溫度模型中的時間步長；Tl為每一個晶格的平均溫度，計算公式為：

(5)

式中KB為玻爾茲曼常數，Ncell為每一個晶格中的原子個數.

2.3 固液狀態的判定

根據文獻[17]中的序參數法來判定原子的固液狀態，每一個原子的序參數為：

(6)

式中兩原子i和j之間的距離rij=|rij|，且rij小于截斷半徑rcut；h為原子i截斷半徑內的原子總數；q為一組倒易矢量[16].

在截斷半徑內，與每個原子鄰近的h個原子的序參數的平均值為：

(7)

根據文獻[10]，當序參數在[0，0.04]時，原子為液態；當序參數在[0.04，1]時，原子為固態.序參數的平均值越大，原子越有序，對于理想晶格，原子的序參數為1；而對于液態物質，原子的序參數接近0.

3 結果與討論

3.1 算法的驗證

為了對算法進行驗證，采用耦合了兩種雙溫度模型的分子動力學方法對脈寬tp=200 fs、能量密度J=5 J/m2的激光照射12.24 nm厚金箔的傳熱過程進行模擬研究，模型初始溫度T1=250 K.

兩種模型中電子和晶格溫度均在12 ps左右達到平衡，在不考慮非傅立葉效應時，達到溫度平衡后金箔的溫度T2=313 K；考慮非傅立葉效應時，達到溫度平衡后金箔溫度T3=320 K.

金箔在整個過程中吸收的能量：

Q=(1-R)·J·S=(1-0.6)×

5×(3×4.08×10-10)2=2.996×10-18J

(8)

式中R為金箔表面的反射率[16]，S為激光照射到金箔的面積.

能量守恒方程為[18]：

Q=cp·m·ΔT

(9)

3.2 拋物線模型和雙曲線模型的結果對比

論文進而分別采用耦合了拋物線模型和雙曲線模型的分子動力學方法對脈寬為200 fs、能量密度為50 J/m2的激光照射12.24 nm厚金箔的固液相變過程進行了模擬研究.

圖2為兩種模型中金箔原子在60 ps時的平均序參數分布圖.當平均序參數在[0，0.04]時，金箔為液態；平均序參數在[0.04，1]之間時，金箔為固態.圖2(a)是拋物線模型的平均序參數分布圖，從圖中可以看出當0≤z≤105×10-10m時，金箔為固態；當105×10-10m≤z≤131×10-10m時，金箔為液態.圖2(b)是雙曲線模型的平均序參數分布圖，從圖中可以看出當0≤z≤71×10-10m時，金箔為固態；當71×10-10m≤z≤131×10-10m時，金箔為液態.在60 ps時，雙曲線模型的熔化深度大于拋物線模型.

圖3(a)給出了兩種模型下金箔固液界面位置隨時間的變化情況，可以看出熔化初始階段固液界面下降明顯，然后熔化速度減慢，而后熔化速度又再次加快，直到金箔完全熔化.金箔熔化速度的改變是由于金箔中心生成了激光引發的壓力[19]，該壓力向上下兩個表面傳播，中心壓力高并逐漸向兩側減小，從而引起了熔化速度的變化.圖3(b)所示為固液界面溫度隨時間的變化情況，初始階段，界面溫度隨時間增加而升高，而后又逐漸降低.這是由于模擬中激光脈沖的時間只有800 fs，激光先將能量迅速傳遞給金箔表面的電子，而后金箔表面的電子以緩慢的速度將能量傳遞給金箔表面的晶格，直至晶格溫度升高并達到熔點從而開始熔化.在這個過程中，晶格中的能量不斷向金箔深處傳遞，于是固液界面溫度呈現出先升高達到峰值后再下降的趨勢.

圖2 兩種模型中金箔原子的平均序參數隨金箔厚度的分布情況Fig. 2 Distribution of average order parameter versus gold thickness in two models

從圖3中還可以看出在同一時刻，雙曲線模型的熔化深度和界面溫度明顯大于拋物線模型.雙曲線模型下金箔在10 ps時開始熔化，在18 ps時固液界面溫度達到最高值1200 K；拋物線模型下金箔在19 ps時開始熔化，在21 ps時固液界面溫度達到最高值1153 K.在110 ps時，雙曲線模型下的金箔完全熔化；在250 ps時，拋物線模型下的金箔仍未完全熔化.造成這一結果的原因是(2)式中的最后一項，把S-G(Te-Tl)看成S’，

圖3 兩種模型下金箔固液界面位置(a)和溫度(b)隨時間的變化情況Fig. 3 Solid-liquid interface (a)location and (b)temperature versus time in two models

3.3 非傅立葉效應下的相變傳熱研究

在非傅立葉效應下對脈寬為200 fs、能量密度為50 J/m2的激光照射12.24 nm厚金箔的固液相變過程進行了模擬研究，圖4為120 ps內金箔中電子和晶格的平均溫度隨時間的變化情況.從圖中可以看出當激光作用于金箔后，金箔內的電子溫度急劇升高并達到最大值，然后隨時間逐漸降低；晶格溫度首先隨時間緩慢增加，當金箔在27 ps左右達到溫度平衡后晶格溫度產生隨時間產生周期性波動.

圖4 金箔中電子和晶格的平均溫度隨時間的變化情況Fig. 4 Average temperature of electrons and lattices in gold foil versus time

不同時刻金箔原子的微觀排列如圖5所示.從圖中可以看出，原子在初始平衡狀態下排列比較規則，在各自的原胞內小范圍振動.隨著激光能量的加入，金箔的溫度逐漸升高，同時上層原子的晶格結構最先遭到破壞，原子排列逐漸變得不規則，金箔開始熔化.隨著時間推移，這種不規則排列開始向金箔底部擴散，固液界面也隨之向金箔底部移動，最終在110 ps時，金箔全部熔化.

論文進而研究了非傅立葉效應下激光參數對金箔熔化過程的影響.圖6(a)所示為在非傅立葉效應下，不同的激光能量密度下金箔固液界面位置隨時間的變化情況.由圖可知，當激光能量密度較小(30 J/m2)時，金箔只有表面很淺一層發生熔化，其他部分未發生熔化，隨著激光能量密度的不斷增大，晶格系統吸收的熱量也不斷增加，從而金箔的熔化深度不斷加大，熔化速度也逐漸加快，當激光能量密度J=50 J/m2時，金箔在110 ps全部熔化.在同一時刻，能量密度越大，熔化深度越大.圖6(b)所示為不同能量密度下金箔固液界面溫度隨時間的變化情況.在同一時刻，激光能量密度越大，金箔中晶格系統吸收的熱量就越多，界面溫度達到的峰值也越高，并且達到峰值所需的時間也越短.

圖5 不同時刻金箔原子的微觀排列圖 (a)t=0 ps(b)t=20 ps(c)t=40 ps(d)t=60 ps(e)t=80 ps(f)t=100 ps(g)t=110 psFig. 5 Microscopic projections of gold foil at different moments (a)t=0 ps(b)t=20 ps(c)t=40 ps(d)t=60 ps(e)t=80 ps(f)t=100 ps(g)t=110 ps

圖6 不同能量密度下金箔固液界面位置(a)和溫度(b)隨時間的變化情況Fig. 6 (a)Location and (b)temperature of solid-liquid interface versus time at different laser flux densities

圖7 不同脈沖寬度下金箔固液界面位置(a)和溫度(b)隨時間的變化情況Fig. 7 (a)Location and (b)temperature of solid-liquid interface versus time at different laser pulse widths

激光的脈沖寬度也是影響激光燒結過程的一個重要因素.圖7所示為不同脈寬下金箔固液界面位置和溫度隨時間的變化情況.在同一時刻，隨著脈寬的增大，金箔熔化深度越小，總的熔化時間越長，而界面溫度也越小.

4 結 論

本文利用耦合雙溫度模型的分子動力學方法從微觀角度對飛秒激光照射金箔的固液相變過程進行了模擬研究，在本文研究的參數范圍內獲得以下結論：

(1)當飛秒激光照射金箔時，金箔的晶格結構隨著溫度的升高逐步遭到破壞，金箔開始熔化，固液界面隨時間逐漸向金箔底部移動，界面溫度隨能量在金箔中的傳遞先升高后降低.

(2)在相同條件下，考慮非傅立葉效應時，金箔熔化的時間較早，并且固液界面溫度和熔化深度比不考慮非傅里葉模型的計算結果大.

(3)較大的激光能量密度和較小的激光脈沖寬度會導致較高的固液界面溫度和較大的熔化深度，且熔化時間較短.