超快激光加熱薄膜的格子玻爾茲曼模擬

毛煜東呂慧麗孫浩森陳明九于明志

(山東建筑大學(xué)熱能工程學(xué)院，山東濟南250101)

0 引言

近年來，激光加熱技術(shù)已展現(xiàn)出強大優(yōu)勢，并廣泛應(yīng)用于激光焊接、激光鉆孔、激光熔覆、激光照明、三維印刷技術(shù)以及熔化動力學(xué)等諸多領(lǐng)域[1-5]。由于該技術(shù)具有高效率、高功率密度、最小的附帶材料損壞、較低的燒蝕閾值等優(yōu)點而受到了越來越多的關(guān)注。目前該技術(shù)在材料的熱加工過程中具有很好的應(yīng)用，如激光清洗、激光圖案化以及超音速激光沉積等[6-8]?，F(xiàn)代微加工技術(shù)要求微結(jié)構(gòu)的高精度加工生產(chǎn)能力，以及材料熱處理中加熱時間迅速和位置的高精度控制，探索微結(jié)構(gòu)的傳熱過程是十分重要的。1822年，傅立葉經(jīng)實驗研究導(dǎo)熱過程并對其歸納總結(jié)，得到了描述宏觀導(dǎo)熱過程的基本定律—傅立葉導(dǎo)熱定律，并逐步地應(yīng)用到了機械、冶金、建筑以及電氣等領(lǐng)域[9]。然而，傅立葉定律暗含熱以擴散的方式傳播且熱傳播速度無限大，即如果對導(dǎo)熱體內(nèi)某一點施加熱擾動，則其他部分會同時感受到此熱擾動帶來的溫度變化，即傅立葉定律忽略了熱傳導(dǎo)過程中溫度梯度和熱流矢量之間的弛豫時間。

玻爾茲曼輸運方程BTE(Boltzmann Transport Equation)用于解決微尺度系統(tǒng)中的熱輸運問題[10]是十分有效的手段。BTE是統(tǒng)計力學(xué)中用以描述不平衡態(tài)分布函數(shù)演化規(guī)律的方程，但BTE中含有的碰撞算子項是一個非常復(fù)雜的非線性微分積分，極難以求解。格子玻爾茲曼方法 LBM(Lattice Boltzmann Method)是BTE常用的數(shù)值解法之一，LBM依據(jù)非平衡統(tǒng)計力學(xué)的運動論，在解決邊界，實現(xiàn)編程和處理復(fù)雜的多尺度耦合問題方面具有明顯的優(yōu)勢，是求解BTE的有效手段，受到了廣泛的關(guān)注[11-12]?；贚BM而來的格子BGK模型LBGK(Lattice Bhatnagar-Gross-Krook)簡化了BTE中的碰撞算子項，用弛豫時間近似表達(dá)了不同粒子間的碰撞達(dá)到的復(fù)雜機制，并且不考慮外力的影響。由于計算效率高，能精準(zhǔn)地推導(dǎo)出納維—斯托克斯方程N-S(Navier-Stokes equations)，常用于研究不可壓縮流體和微尺度的導(dǎo)熱問題，也是LBM主要的研究模型。

近年來，LBM廣泛應(yīng)于描述微納米尺度傳熱的問題中。Escobar等[13]應(yīng)用LBM研究了一維半導(dǎo)體材料中的瞬態(tài)多空間尺度和時間尺度的導(dǎo)熱過程，與基于擴散-彈道輸運方程、傅立葉導(dǎo)熱方程和CV模型的結(jié)果做了對比，并對主要的導(dǎo)熱模型的使用范圍做了簡要總結(jié)。華鈺超等[14]用熱質(zhì)理論的觀點分析了彈道擴散導(dǎo)熱機理，基于聲子玻爾茲曼方程推導(dǎo)了修正邊界條件模型，數(shù)值求解了修正的普適導(dǎo)熱定律，并與蒙特卡羅模擬進行對比。張珂等[15]基于LBM和雙曲兩步模型結(jié)建立了一個LBM兩步方程，并利用此方法模擬研究了納米薄膜在超快激光照射過程中的熱響應(yīng)特性，分析了照射過程中薄膜內(nèi)溫度隨時間及空間的變化規(guī)律，探討了激光強度以及薄膜厚度對金屬薄膜熱響應(yīng)的影響。文章應(yīng)用LBM方法對超快速激光加熱硅薄膜的一維導(dǎo)熱問題進行研究，來描述這種時間超短，尺寸超小的超快速納米傳熱問題，展示薄膜內(nèi)部能量密度分布情況，并對結(jié)果進行分析。

1 理論分析

LBGK模型的表達(dá)式由方程(1)表示為

式中:f(r，ξ，t)為分子速度分布函數(shù)；r為空間位置矢量；ξ為分子速度矢量；t為時間參數(shù)；feq為熱力學(xué)平衡態(tài)時的速度分布函數(shù)；v為熱載子的群速度；τ0為松弛時間(也叫弛豫時間)是粒子發(fā)生2次碰撞的時間間隔，表征粒子碰撞后達(dá)到平衡態(tài)的快慢程度。

1.1 LBM方法驗證

傅立葉導(dǎo)熱定律可以很好地描述宏觀問題，文章采用LBM方法解決微觀—宏觀問題，應(yīng)先驗證LBM在解決宏觀熱傳導(dǎo)問題時的有效性。建立一個簡單的薄膜導(dǎo)熱問題，即在薄膜左側(cè)施加一個無量綱數(shù)為1的熱擾動，其控制方程由式(2)表示為

式中:Ui為能量密度；x為位置參數(shù)；U0i為平衡能量密度；vi為聲子群速度；τ為聲子弛豫時間。

考慮常溫初始條件，將式(3)～(6)所示的無量綱變換為

式中:x?為無量綱位置參數(shù)；L為薄膜的特征長度；Cv為硅的體積熱容；t?為無量綱時間參數(shù)；U?為無量綱能量密度參數(shù)；U1?為向x軸正方向傳輸?shù)臒o量綱能量密度；為向x軸負(fù)方向傳輸?shù)臒o量綱能量密度；T0為初始時刻溫度；E1為初始溫度T0下的能量。

在做與式(3)～(6)相同的無量綱變換后，其邊界條件由式(7)表示為

1.2 應(yīng)用LBM模擬超快激光加熱薄膜問題

建立基于LBM的一維熱傳導(dǎo)模型來模擬超快激光在硅薄膜中的能量分布。在硅薄膜中，聲子是主要的能量載體?；贐TE，考慮超快激光加熱薄膜問題的一維導(dǎo)熱模型由式(8)表示為

式中:i=1，2；S為能量吸收率，可由式(9)[16]表示為

式中:J為激光能量發(fā)射密度；R為表面反射率；tp為激光脈沖的持續(xù)時間；δ為激光穿透深度。

使用LBM離散式求解方程(8)，忽略高階項，并做與式(3)～(6)相同的無量綱變換后，則主導(dǎo)方程可由式(10)和(11)表示為

式中:W=Δt?=Δx?/Kn；ξ=L/δ；γ=τ/tp；W、ξ、γ均為中間參數(shù)，無特殊含義Kn=l/L為克努森數(shù)；l為聲子平均自由程。

根據(jù)傅立葉定律，將上述的超快激光加熱問題由經(jīng)典熱傳導(dǎo)方程式進行描述，由式(12)表示為

式中:kF為硅的導(dǎo)熱系數(shù)。

基于相同邊界條件，可將函數(shù)的解展開，由式(13)～(15)表示為

1.3 應(yīng)用LBM模擬雙側(cè)超快激光加熱薄膜問題

利用相同的超快激光脈沖同時從硅薄膜的兩側(cè)表面進行加熱，并設(shè)置與單側(cè)加熱相同的鏡面反射邊界條件和常溫初始條件。此時描述該問題的BTE控制方程由式(16)～(18)表示為

將式(17)和(18)代入式(16)中，可將控制方程變?yōu)槭?19)形式

利用LBM對式(19)進行離散求解，方法與式(10)和(11)相同，不再贅述。

2 數(shù)值結(jié)果與分析

為便于分析，超快脈沖激光加熱的硅膜采用以下參數(shù):體積熱容Cv為1.66×106J/(m3·K)，硅的聲子平均自由程l為41 nm，聲子弛豫時間τ為6.53 ps，持續(xù)時間tp為 0.65 ps，激光通量J為312 J/m2，穿透的光學(xué)深度δ為15.3 nm，以及激光的表面反射率R為0.93。

2.1 驗證LBM在微觀—宏觀問題的有效性

分析式(2)～(7)描述的薄膜導(dǎo)熱問題，分別利用LBM和經(jīng)典傅里葉熱傳導(dǎo)理論對該問題進行求解和模擬，觀察薄膜內(nèi)部的能量傳遞現(xiàn)象，如圖1所示?？梢钥闯?，通過LBM獲得的結(jié)果清楚地表現(xiàn)出擴散行為，并顯示出良好的擴散性能，這與傅立葉導(dǎo)熱理論所給出的結(jié)果一致。因此，LBM可用于模擬宏觀擴散狀態(tài)下的熱傳導(dǎo)。

圖1 利用LBM和傅立葉定律得到的Kn=0.001的無量綱能量密度分布圖

改變克努森數(shù)，將熱傳導(dǎo)過程引入到過渡區(qū)，并觀察LBM和傅立葉定律得到的結(jié)果之間的差異，如圖2所示。當(dāng)克努森數(shù)為1時，此時薄膜內(nèi)熱傳遞處于過渡區(qū)，基于LBM和傅立葉定律的解決方案表現(xiàn)出截然不同的行為?？梢钥闯?，當(dāng)無量綱時間為0.1時，通過LBM模擬獲得的能量密度分布具有顯著的特征:能量密度在受加熱側(cè)(左側(cè))邊界有明顯的跳躍；隨著時間的流逝，由LBM獲得的熱能以熱波的形式在薄膜內(nèi)部傳輸；當(dāng)t?達(dá)到1時，能量密度又在右側(cè)邊界發(fā)生了的躍變。然而，傳統(tǒng)的傅立葉定律卻不能捕捉這一現(xiàn)象，由于其內(nèi)含熱傳遞速度無限大的假設(shè)，在過渡區(qū)，薄膜內(nèi)任何一處都可以立即感受到受加熱側(cè)施加的熱擾動。

圖2 利用LBM和傅立葉定律得到的Kn=1的無量綱能量密度分布圖

2.2 單側(cè)激光加熱納米薄膜導(dǎo)熱問題

分析式(8)和(9)所描述的利用超快激光加熱厚度為82nm薄膜的導(dǎo)熱問題。通過傅立葉定律和LBM方法得到的無量綱能量密度分布如圖3所示。從圖3(a)中可以看出，基于傅立葉定律的熱傳導(dǎo)是一個擴散過程。由于強大的激光脈沖加熱，短時間內(nèi)能量密度就在薄膜的左側(cè)區(qū)域迅速向右側(cè)擴展。隨著時間的流逝，能量密度逐漸降低，并最終在薄膜內(nèi)趨于穩(wěn)定。

圖3 激光加熱薄膜單側(cè)時薄膜內(nèi)無量綱能量密度分布圖

然而，對于相同的加熱問題，LBM卻得到了不同的結(jié)果，如圖3(b)所示。在初始階段，同傅立葉定律獲得的結(jié)果相比，受到脈沖激光擾動的左邊界附近產(chǎn)生了更高的能量密度。該無量綱密度峰值達(dá)到了1.18，較傅立葉結(jié)果的0.8相比，超出約47.5%。這表明，傳統(tǒng)的傅立葉熱定律在描述這此類微納米尺度下的超快加熱問題時會嚴(yán)重低估熱載子傳輸?shù)淖畲竽芰棵芏取?/p>

隨著時間的流逝，薄膜內(nèi)部能量是以波的形式從左端傳遞到右端，且隨著能量波的傳遞，其能量密度峰值逐漸下降，這一現(xiàn)象并未出現(xiàn)在傅立葉定律得出的結(jié)果中。LBM的結(jié)果還發(fā)現(xiàn)，在鏡面反射的邊界條件下，能量密度會經(jīng)過一個較傅立葉結(jié)果更為漫長的時間后才能在薄膜內(nèi)趨于穩(wěn)定。

2.3 雙側(cè)激光加熱納米薄膜導(dǎo)熱問題

分析式(16)～(19)所述利用相同的超快激光脈沖同時從硅薄膜的兩側(cè)表面進行加熱的問題。在克努森數(shù)為0.5時，利用傅里葉定律和LBM方法得到的薄膜內(nèi)部的能量密度分布情況如圖4所示。如同預(yù)期的那樣，傅立葉的結(jié)果顯示了熱傳導(dǎo)的擴散行為，也并未發(fā)現(xiàn)能量波狀傳遞的行為，如圖4(a)所示。初始時刻，左右兩側(cè)的無量綱能量密度峰值均為約0.8，同單側(cè)激光加熱時相同，經(jīng)過無量綱時間2后，薄膜內(nèi)部能量密度達(dá)到穩(wěn)定，此時無量綱能量密度值為0.5。

由圖4(b)可以看到，在初始時刻受激光脈沖擾動的兩個邊界附近的區(qū)域能量密度的快速升高(即溫度快速提升)；隨后，從左右邊界發(fā)起的2組熱波從不同方向同時向硅膜的中間區(qū)域傳輸，在此期間，能量密度峰值逐漸降低。當(dāng)無量綱時間達(dá)到1時，2組熱波峰在薄膜的中間區(qū)域匯合，相互的碰撞作用導(dǎo)致了薄膜中心區(qū)域的能量密度顯著升高。此時，在中心點處的能量密度的峰值為1.095，約為圖3(b)所示的單側(cè)激光加熱情況下所預(yù)測的峰值的兩倍，這一現(xiàn)象不僅是傅立葉結(jié)果無法捕捉的，而且在單側(cè)激光加熱時也不會出現(xiàn)。最后，當(dāng)無量綱時間達(dá)到16時，薄膜內(nèi)能量密度趨于穩(wěn)定；相比于傅立葉結(jié)果，在鏡面反射的邊界條件下，趨于穩(wěn)定所需時間要更為漫長。

圖4 激光加熱薄膜兩側(cè)時薄膜內(nèi)無量綱能量密度分布圖

3 結(jié)論

文章應(yīng)用LBM方法和傅立葉定律對超快激光加熱硅薄膜的一維導(dǎo)熱問題進行研究，分析了薄膜內(nèi)部能量密度分布情況。主要結(jié)論如下:

(1)發(fā)現(xiàn)在克努森數(shù)較小的連續(xù)介質(zhì)區(qū)，LBM得到的結(jié)果符合傅立葉定律給出的結(jié)果；然而在系統(tǒng)特征尺度和分子平均自由程相當(dāng)?shù)倪^渡區(qū)，LBM方法可以捕獲到能量以波動的方式進行傳輸?shù)臒岵ㄌ匦?，但傳統(tǒng)的傅立葉定律無法得到。

(2)利用LBM數(shù)值模擬激光加熱硅薄膜的一維導(dǎo)熱問題，發(fā)現(xiàn)克努森數(shù)對薄膜內(nèi)無量綱能量密度的分布具有重要的影響。

(3)如果利用超快激光分別在薄膜的兩側(cè)進行加熱，薄膜內(nèi)部的能量密度分布情況將有所改變，尤其是當(dāng)沿相反方向傳播的兩個熱波交匯在一起時，會發(fā)生碰撞并使得能量密度有顯著的提升，這些現(xiàn)象均在傅立葉定律結(jié)果中不能呈現(xiàn)。此外，研究發(fā)現(xiàn)傳統(tǒng)的熱傳導(dǎo)理論可能嚴(yán)重低估了過渡區(qū)聲子輸運的最大能量密度。