應用D2Q9 格子玻爾茲曼模型研究超快激光傳熱過程

趙國晨孫浩森王先征于明志毛煜東

(1.山東建筑大學熱能工程學院，山東 濟南 250101；2.山東建筑大學學報編輯部，山東 濟南 250101)

0 引言

超快激光技術因具有加熱速度快、能量密度和精準性高等優點而廣泛應用于激光焊接[1]、激光加工[2-3]等諸多領域。 硅薄膜是一種比較常見的半導體材料[4]，在芯片[5]、太陽能電池[6-7]和各種電子系統[8]等領域都具有十分廣泛的應用。 借助超快激光加熱硅薄膜探索微/納米結構的傳熱過程，對于材料加工、熱處理[9]以及微電子系統的散熱[10]等技術具有重要的現實意義。

傅里葉定律可用來描述宏觀上的傳熱過程，但其假設熱的傳播速度無限大且是以擴散的方式在材料中傳播，對于超快速加熱過程并不適用。 格子玻爾茲曼方法(Lattice Boltzmann Method， LBM)是玻爾茲曼輸運方程(Boltzmann Transport Equation， BTE)常用的數值解法，具有計算工作量較小、邊界處理簡單、易于編程等優點，非常適合于求解微納米尺度下的熱輸運問題。 BTE 方程具有非常復雜的碰撞項，無法對其直接求解。 1954 年，BHATNAGAR、GROSS 和KROOK在不影響求解結果的前提下引入了一個簡單的碰撞算子模型簡化了BTE 模型，稱為LBGK 模型，成為后來LBM 的主要研究模型[11]。

硅薄膜作為半導體材料，其內部熱量傳遞的熱載子為聲子[12]。 德拜模型認為所有聲子都以單一、固定的速度傳播，而與聲子的頻率無關，將所有聲子分支壓縮成一個線性色散關系[13]。 在遵循德拜模型假設的前提下推導出了灰色LBM 模型(Gray LBM Model)，忽略了聲子特性的頻率依賴性，大大簡化了求解過程，因此也成為了模擬聲子輸運的主要模型。 文章應用灰色LBM 模型研究了超快激光加熱硅薄膜的二維導熱問題，探究了這種超短時間、超小尺寸的傳熱問題，展示了二維薄膜內部能量密度的分布情況，并對比分析了一維結果[14]。

1 理論分析

1.1 晶格模型

LBM 主要由玻爾茲曼輸運方程、晶格、格點至格點的傳輸限制等部分組成[15]。 晶格可以描述為對研究區域的空間離散化，進而可以擴展為對速度域進行離散化，因為粒子或粒子群在LBM 中只允許沿著連接兩個連續相鄰節點的方向移動。 兩個常用的二維晶格示意圖如圖1 所示。

圖1 二維晶格模型示意圖

D2Q5 是二維晶格中最簡單的模型，其每個格點只有4 個相鄰點，所以僅有4 個與之相關的離散傳播速度，這還不足以代表粒子可以傳播的所有可能方向。 為了取得更加精確的結果，文章選用傳播速度更多的D2Q9 模型，其每個格點有8 個相鄰點交換能量，可以很好地模擬二維傳熱問題[16]。

1.2 格點到格點的傳輸限制

格點到格點的傳輸限制迫使粒子在離散時間步長內從一個晶格位置跳到相鄰晶格位置，從而有效地將格點到格點的距離d與離散時間步長dt的大小相關聯，如圖2 所示，其中c為離散速度。 這一限制有助于簡化玻爾茲曼方程的離散推導。

圖2 格點到格點的傳輸

1.3 控制方程的推導和離散

Boltzmann-BGK 模型的表達式由方程(1)表示為

式中f為速度分布函數；v為粒子的群速度；τ為弛豫時間(粒子發生兩次碰撞的時間間隔)；?eq為平衡狀態下的速度分布函數。

由聲子能量密度表示的BTE 可由式(2)表示為

式中e為能量密度；e0為平衡能量密度；x和y為位置參數；τ0為聲子弛豫時間。

文章采用的D2Q9 模型，共有9 個離散速度，控制方程可由式(3)表示為

式中i＝0，1，2，…，8 為不同的速度方向。 當i＝0 時，；當i＝1，2，3，4 時，；當i＝5，6，7，8時，

為方便計算進行下列無量綱變換，由式(4)～(6)表示為

式中x*、y*為無量綱位置參數；L、M分別為薄膜厚度方向和長度方向的特征尺寸；t*為無量綱時間參數；e*為無量綱能量密度參數；bi為不同方向的比重系數，相加和為1；E1為初始溫度下的能量，E1＝e(T0)＝CvT0，其中Cv為硅的體積熱容；T0為初始溫度。

基于LBM 建立超快激光加熱薄膜問題的二維導熱模型，由式(7)表示為

式中S(x，y，t) 為激光能量的吸收率，可由式(8)表示為

式中J為激光能量的發射密度；R為表面反射率；tp為激光脈沖的持續時間；δ為激光的穿透深度；r0為激光熱源影響半徑。

將式(4)～(6)的無量綱變化帶入式(7)，得到無量綱化后的控制方程，由式(9)～(17)表示為

離散化后的無量綱化的控制方程可以由式(18)～(26)表示為

2 方法驗證

尺寸效應會對熱導率產生影響，是納米尺度傳熱的顯著特征之一。 ALVAREZ 等[17]基于擴展的不可逆熱力學，提出了一個尺寸相關的有效導熱系數表達式，由式(27)表示為

式中kF為硅的體導熱系數。 從表達式中可以清楚地看出有效導熱系數取決于克努森數，而克努森數反映了納米級導熱系數的尺寸效應。 此外，當平均自由程遠小于器件的特征長度時，有效熱導率近似等于體熱導率。 相反，當器件的尺寸趨于零時，克努森數趨于無窮大，有效導熱系數近似線性地取決于器件的尺寸。 該模型的優點為在較大克努森數范圍內的結果與實驗結果一致[17]。 AMON 等[18]運用LBM 研究了納米尺度硅薄膜中的熱傳輸，表明硅薄膜的面外有效導熱系數取決于克努森數和表面調節系數，其表達式可由式(28)表示為

式中表面調節系數α為擴散反射的熱載體分數，文章取α＝1。 該模型不僅反映了尺寸效應，而且還反映了熱載體在硅薄膜表面的擴散和鏡面散射的影響。

將文章所研方法獲得的熱導率結果與ALVAREZ 等[17]、AMON 等[18]的結果進行比較，如圖3 所示，在0 ～150 nm 厚度范圍內，文章方法所獲得的熱導率與兩者的結果較為吻合，尤其與ALVAREZ 等[17]的結果吻合度較高，而與AMON等[18]的結果誤差最大值≤20 W/(m·K)。

圖3 不同方法獲得的有效導熱系數對比圖

3 結果與分析

3.1 超快激光加熱硅薄膜厚度方向傳熱問題

采用絕熱邊界條件，涉及到的參數見表1。

表1 硅薄膜與超快激光基本參數表

為了方便與一維結果進行對比，將被激光照射后薄膜的溫度分布轉換為無量綱能量密度分布，結果如圖4 所示，其中圖4(a)～(c)為激光開始照射后不同時間的薄膜被加熱中心位置的厚度方向的無量綱能量密度分布圖。 超快激光在無量綱時間t*＝0 時照射薄膜一側的中心位置，此時被加熱區域的無量綱能量密度開始快速上升并且向薄膜內部傳遞。 在激光持續加熱階段，被加熱面的無量綱能量密度快速升高，而激光對薄膜內部的影響隨著厚度方向的深入逐漸減弱。 這是由于超快激光在很短的時間內將大量的能量傳遞到薄膜，使得薄膜靠近被加熱面區域的能量密度急劇升高，此時能量在薄膜內部傳遞的速度遠遠小于激光傳遞給薄膜能量的速度。

圖4 硅薄膜被激光加熱時厚度方向無量綱能量密度分布圖

Kn＝0.5 時無量綱能量密度的分布圖如圖4(a)所示，此時薄膜的厚度為82 nm，固定薄膜的縱向長度為薄膜厚度的10 倍。 當無量綱時間t*為0.04、0.10時，二維LBM 與一維LBM 模擬的結果基本吻合，薄膜左側被加熱面附近區域的能量密度明顯高于其他區域，且薄膜靠近右側位置的無量綱能量密度為0，說明此時激光加熱還未影響到薄膜右側位置。 在這兩個時間點時激光還處在加熱薄膜階段，被加熱位置的能量密度快速上升，此時由于時間較短，激光主要作用在薄膜厚度方向，而在縱向上的能量傳遞較小。 激光對薄膜的穿透深度有限，加熱主要作用在被穿透區域，所以在時間較短時激光對薄膜右側不會有明顯的影響。 當t*為0.60 和1.20時，可以看到能量在薄膜內部是以波的形式傳遞，且隨著能量的傳遞無量綱能量密度的峰值是逐漸減小的。 此時的一維與二維模擬的結果雖然在相同時刻時在薄膜內部傳遞的規律及波峰出現的位置是基本相同的，但是二維LBM 模擬的無量綱密度的數值已經明顯小于一維LBM 模擬的結果。 當激光開始照射薄膜后引起薄膜左側的能量急劇升高并開始沿厚度方向和縱向傳遞，相對于能量的升高速度，能量的傳遞會有一定的遲滯，這一點是被傅里葉模型忽略的，所以能量在傅里葉模型下都是以擴散形式傳遞的。 由于模擬采用的是絕熱邊界條件，所以當激光照射結束后，能量在薄膜內部的傳遞是守恒的。 一維LBM 忽略了能量在薄膜縱向上的傳遞，能量僅在橫向傳遞，隨著時間增加，二維LBM 會有越來越多的能量傳遞到縱向，直到整個薄膜達到平衡。

Kn＝1.0 和2.0 時無量綱能量密度的分布分別如圖4(b)和(c)所示，此時的薄膜厚度分別為41.0、20.5 nm。 隨著薄膜厚度的減小，激光對薄膜厚度方向上的影響也越來越大，薄膜內部能量傳遞的波動性減弱，能量密度達到穩定所需的時間也越來越少。對比圖4(a)～(c)發現薄膜最后達到穩定時的能量密度隨著克努森數的增大而逐漸增大。 為了進一步探究克努森數對薄膜內能量密度分布的影響，模擬了t*＝0.20 時不同克努森數的薄膜被加熱中心位置厚度方向的無量綱能量密度分布，如圖4(d)所示。 可以看出，克努森數越大其能量密度峰值也就越大，峰值達到的無量綱位置越遠。

一維LBM 與二維LBM 在不同克努森數情況下模擬的能量傳遞規律相同，二維LBM 考慮了能量的縱向傳遞，導致橫向傳遞的能量出現大幅度的減小，說明在研究薄膜中能量的傳遞時不能將縱向傳遞忽略。

3.2 超快激光加熱硅薄膜縱向傳熱問題

在不同t*時刻，薄膜被加熱面縱向能量密度分布如圖5(a)～(c)所示。 可以觀察到，在激光加熱階段被加熱的位置能量密度快速上升，而離加熱位置較遠的區域能量密度并沒有明顯變化，加熱結束后能量隨著時間逐漸向外傳遞，并展現出了與厚度方向上類似的波動傳遞的現象。 由于薄膜的縱向尺寸隨厚度減小而減小，而激光加熱半徑不變，所以在圖中可以觀察到隨著克努森數的增大，被加熱區域占總區域的比例也逐漸增大。 無量綱能量密度在被加熱的中心位置最高，并向兩側逐漸減小。 隨著時間的推移能量不斷向兩側傳遞，這一點與在厚度方向上的傳遞規律類似。 而相同時刻，薄膜厚度方向上的能量波動要大于縱向，在縱向上傳遞的能量是遠小于厚度方向的，所以薄膜縱向上堆積的能量會小于厚度方向，進而造成能量波動小于厚度方向。薄膜的縱向尺寸大于厚度方向，能量密度在縱向上達到穩定需要的時間要長于厚度方向。

圖5 硅薄膜被激光加熱時縱向無量綱能量密度分布圖

為了研究在某一時刻薄膜內部縱向的能量密度分布，模擬了當Kn＝1.0、t*＝0.10 時薄膜內部縱向的能量密度分布，結果如圖5(d)所示。 可以看出，薄膜在這一時刻，縱向上的能量密度隨薄膜厚度的深入而減小，影響范圍也逐漸變小，這是因為激光對薄膜的穿透深度會沿照射中心位置向兩側逐漸變淺，攜帶的能量也逐漸變小。 由此可見薄膜內部縱向與厚度方向上的能量傳遞是會相互影響的。

4 結論

基于玻爾茲曼輸運方程，采用D2Q9 的二維LBM 模型模擬了超快激光加熱硅薄膜的能量傳遞，主要結論如下:

(1) 薄膜被激光照射后，靠近加熱位置的區域能量密度快速升高，隨著時間的推移，影響的區域逐漸增大。 在薄膜被加熱區域附近，激光照射為能量來源，而薄膜內部的能量來源主要來自熱傳導。

(2) 能量在薄膜內部以波的形式傳播，在厚度方向和縱向上的傳遞都展現出了明顯的波動性，導致能量的峰值出現在了薄膜的內部，其位置會隨著時間而發生改變。

(3) 雖然當激光做熱源時在厚度方向上攜帶的能量要多于縱向，但厚度方向與縱向上的能量在薄膜內部傳遞時會相互影響、相互傳遞，會與一維的模擬結果有較大差距，所以在研究薄膜傳熱問題時應同時考慮厚度方向和縱向上的傳遞，不能簡單近似為一維傳熱問題。