中華大學碩士論文

(1)

中華大學碩士論文

題目：極快速脈衝雷射加熱過程之二維熱傳遞研究 An Analysis on Two Dimensional Heat Transfer During Ultrashort Laser Heating of Thin Films

系所別：機械與航太工程研究所學號姓名：M09208024 張毓正指導教授：許隆結博士

中華民國九十四年七月

(2)

(3)

(4)

(5)

摘要

本論文針對材料受極短脈衝雷射加熱過程中二維熱傳現象，以雙相延遲模式模擬加熱過程，對溫度分佈情形進行詳細的研究。並分析包括熱通量延遲時間、溫度梯度延遲時間、雷射脈衝半徑等參數對熱傳遞的影響。文中定義一無因次比值 B，其物理意義為溫度梯度延遲時間除上 2 倍的熱通量延遲時間。由徑向分析的結果可得知，若以傳統傅立業熱傳模型，或熱波模型來模擬極快速雷射脈衝加熱過程，其徑向溫度分佈與能量傳遞都異於雙相延遲模型。軸向分析的結果指出二維效應對溫度預測的影響。其結果並顯示若以 95nm 的雷射半徑對金的薄膜加工，相對於一維模型的溫度預測，在時間為一倍的脈衝時間時，已有 29%的相對誤差。而在時間等於十倍的雷射脈衝時，其相對誤差更高達 70%。針對固定位置溫度隨時間的變化，結果顯示，一維的溫度預測將高於二維所模擬之情形，且隨時間增加而差距逐步擴大。

(6)

Abstract

This study is intent to investigate the effect of two dimensional heat transport phenomena of a thin film subjected to femtosecond laser heating.

The dual phase lag model (DPL) was applied to simulate the temperature distribution. The effects of radius of the pulse laser and ratio of delay times on the temperature distributions were discussed in detail.A nondiemensional parameter B was defined as the ratio of time delays of temperature gradient over twice of time delay of heat flux. The radial directional analyses showed that both temperature and energy distributions predicted by DPL model were totally different from those predicted by diffusion and CV wave model. The axial directional analyses also showed a significant two dimensional effect on the temperature distribution. When gold film subjected to pulse laser heating with 95 nm beam radius, results also indicated that one dimensional analyses cause the relative errors on temperature distribution to be about 29%. The error could increase to 70% as time increases to 10 times of the laser pulse duration.

(7)

誌謝

動筆的開始，也正象徵著研究所生活的結束，難免有些不捨。想起當初學長的一席話，令我至今都奉為圭臬。研究所的深造，並不僅於大學教育的延續，更進一步在培養獨立思考與解決問題的能力。回顧兩年，我不斷期許自己朝這目標邁進。畢竟知識代表個人的內涵，

而態度才是左右人生的關鍵。故我引村上春樹的話來勉勵後進學弟。

『問，是一時之恥；不問，是一生之恥。』

得以完成這本論文，要感謝的人實在太多，首要感謝恩師許隆結博士。老師不但有計畫地指導，不厭其煩地督促與提點，更適時地在學生遇到重重難關時給予幫助。對於這份恩情，學生永遠銘記在心，

並竭誠希望學生的表現能讓老師滿意。再來要感謝口試委員邱弘興博士、陳鎮憲博士與朱朝煌博士。感謝他們對我論文的細心批閱，並指出總總缺失，使得這本論文更加完備。

接著要感謝學長霍呈翔、孟憲虹及學姊廖詩茵，他們的親切與熱心讓我能很快地融入研究所的一切。同學范牧樹、郭子瑋與陳亭棋更讓我有兩年精彩難忘的研究所生活點滴。學弟葉祖銘、廖乙安、李森正與劉佳原，謝謝你們容忍學長的所有無理取鬧的行為。最後要感謝我最親愛的家人們這二十多年來的辛勞與無怨無悔的付出，沒有你們就沒有今天的我，僅以此文獻給你們。

(8)

目錄

頁次

摘要...i

英文摘要...ii

誌謝...iii

目錄...iv

圖目錄...vi

表目錄...vii

符號說明...viii

第 1 章緒論...1

1.1 研究動機與目的...1

1.2 文獻回顧...2

1.3 研究方法...4

第 2 章巨觀與微觀熱傳簡介...5

2.1 傳統傅立葉熱傳模型...5

2.2 聲子-電子交相反應模型...6

2.3 聲子散射模型...8

2.4 熱波模型...9

2.5 雙相延遲模型...10

第 3 章脈衝雷射加熱薄膜理論模式...12

3.1 熱傳模型與熱源項...12

3.2 積分轉換法求解...14

3.3 數值方法...18

第 4 章結果與討論...23

(9)

4.1 程式驗證...23

4.2 徑向溫度分佈...24

4.3 軸向溫度分佈...29

4.4 時間軸溫度分佈...35

第 5 章結論...39

參考文獻...41

(10)

圖目錄

頁次

【圖 1】物理模式及座標圖...12

【圖 2】溫度分佈與文獻[19]比較結果...24

【圖 3】γ₀=1.0，β =0.01，B 值改變下溫度分佈圖...25

【圖 4】γ₀=1.0，β=0.1，B 值改變下溫度分佈圖...25

【圖 5】γ₀=1.0，β=0.5，B 值改變下溫度分佈圖...26

【圖 6】γ₀=1.0，β=1.0，B 值改變下溫度分佈圖...26

【圖 7】γ₀=3.0，β=0.5，B 值改變下溫度分佈圖...28

【圖 8】γ₀=10.0，β =0.5，B 值改變下溫度分佈圖...28

【圖 9】γ₀=20.0，β =0.5，B 值改變下溫度分佈圖...29

【圖 10】β =0.1，雷射半徑改變下溫度分佈圖...30

【圖 11】β=0.5，雷射半徑改變下溫度分佈圖...30

【圖 14】B=0.2，β=0.5，雷射半徑改變下溫度分佈圖...32

【圖 17】B=0.0，δ=1.0，γ₀改變下溫度分佈圖...36

【圖 18】B=0.2，δ=1.0，γ₀改變下溫度分佈圖...36

【圖 19】B=0.5，δ=1.0，γ₀改變下溫度分佈圖...37

【圖 20】B=1.0，δ=1.0，γ₀改變下溫度分佈圖...37

【圖 21】B=5.0，δ=1.0，γ₀改變下溫度分佈圖...38

【圖 22】B=10.0，δ =1.0，γ₀改變下溫度分佈圖...38

(11)

表目錄

頁次

【表 1】金屬之聲子-電子偶合因子...7

【表 2】Legendre 多項式之根及係數...21

【表 3】四種金屬之 B 值...27

【表 4】

β

=1β_P、5

β

_P、10

β

_P，與一維比較之相對誤差...34

(12)

符號說明

B 溫度延遲時間與熱通量延遲時間之比值 [

2

T q

τ τ ] qK

熱通量 [W/m²] Cp 定壓比熱 [ ^J₃

m K ] k 熱傳導係數 [W/mK]

τ

T 熱通量延遲時間 [s]

τq 溫度梯度延遲時間 [s]

R 單層介面反射率 L 拉式轉換

s 拉式轉換參數 H Hankel 轉換 σ Hankel 轉換參數

Q 熱源項 [W/m³] t 時間[s]

z 軸向[m]

r 徑向[m]

r0 雷射半徑[m]

d0 雷射直徑[m]

T 溫度[°k]

(13)

T0 室溫[°k]

η 無因次熱源項 [^Q^/

(

^kT⁰^/^ατ^q

)

^]

θ 無因次溫度 [ ⁰

0

T T

T

− ]

β 無因次時間 [

2 _q t τ ]

δ

無因次 z 軸 [

2 _q

z α τ ]

γ 無因次 r 軸 [

2 _q

r α τ ]

0

δp 無因次穿透深度 [

2

p q

δ α τ ] βp 無因次雷射脈衝時間 [

2

p q

t τ ] γ0 無因次雷射半徑 [ 0

2 _q

r α τ ] A 無因次薄膜厚度 [

2 _q

L

α τ ]

上標：

~ 經拉普拉斯轉換經 Hankel 轉換下標：

r 隨徑向變化 z 隨軸向變化 t 隨時間項變化

(14)

第一章緒論

1.1 研究動機與目的

隨著時代進步，科技不斷創新，精密加工技術從微米(10^-6m)到如今的奈米(10^-9m)，使得科技產品關鍵性零組件的精度或尺寸均已達微米或奈米等級，說明了高科技產業已朝向微小化的趨勢邁進。微小化有許多優點，如節省空間、節省材料、低污染或節約能源。然而也正因如此，伴隨微小化而衍生的加工問題也隨之而來，所以精密加工技術的研究便成為高科技產業發展的重要關鍵。

雷射加工技術是利用雷射光束與物質相互作用的特性對材料進行切割、焊接、表面處理、打孔等的一門加工技術，被廣泛應用於金屬、塑膠及其他材料的微加工製程。當材料的幾何尺寸到達微米或奈米等級時，加工的精度常常因為雷射的能量無法準確的控制而產生限制。近年來發展的超短脈衝雷射其雷射脈衝長度(Pulse Duration)可達到飛秒(Femtosecond, 10^-15s) [1] 甚至於原秒(Attosecond, 10^-18s) [2,3]

的等級，經由此一技術可精確控制雷射加工時的能量，避免上述的問題產生。然而在超短脈衝雷射加工過程中，因雷射加熱能量而發生於加熱區/或加工區的材料消融(Ablation) [4,5]，勢必會破壞材料晶格，

以至於影響尺寸與形狀的精度，甚至影響到一些物理性質。此外，利

(15)

用脈衝雷射來進行非破壞性檢測亦是常用的檢測方法，因此脈衝雷射加熱過程的熱傳遞研究變成是很重要的問題。就作者所知，由於統御方程式的複雜度，目前針對薄膜加熱的研究均受限於一維熱傳分析。

故本文運用雙相延遲模型研究皮秒(Picosecond, 10^-12s)等級內的熱傳遞現象。比對以傳統傅立葉模型及熱波模型所模擬的結果，來探討此三種模型間的異同。並詳細探討二維效應對薄膜溫度分佈的影響，指出雷射半徑在微米到奈米等級的加工上使用一維熱傳分析的誤差。

1.2 文獻回顧

由於雷射的脈衝時間在飛秒的等級，因此單一雷射脈衝加熱影響的時間約在幾個皮秒左右，在這麼短的時間內，傳統的傅立葉熱傳模型並無法準確的描述其熱傳遞現象，其原因在於傅立葉熱傳模型假定傳熱的載子(電子或聲子)經過多次的碰撞後，處於局部平衡狀態。但在金屬內傳熱載子碰撞所需的時間約在幾個皮秒左右，在非金屬內則約十個皮秒，因此在這麼短時間內傳熱載子無充裕時間與其他載子做碰撞而達平衡狀態[6-9]。欲描述微小時間內的熱傳遞現象，Kaganov 等人[10]及 Anisimov[11]等人提出兩階段加熱模型(Two-Step Heating Model)來描述金屬的微觀熱傳遞現象。此模型假設在材料受熱後，由於電子的熱容遠比晶格的熱容小，脈衝雷射加熱後，電子會馬上吸收

(16)

此熱量達到非常高的溫度，而晶格仍處於低溫狀態，電子透過電子- 聲子偶合因子(Electron-Phonon Coupling Factor)加熱晶格，經過幾個皮秒的時間後兩者才達到相同的溫度。此模型經過Qiu 等人[12]給予嚴謹的理論推導，並經由飛秒脈衝雷射加熱金屬的實驗量測加以證明其正確性[13-17]。Guyer 及 Krumhansl[18]考慮非金屬內的熱傳遞，

經由波茲曼傳遞方程式推導，忽略電子的傳遞效應討論聲子與聲子之間的碰撞，得到描述非金屬的微觀熱傳方程式，其方程式再經過推導可得到與兩階段加熱模型所得的方程式有相同的形式。Tzou[19]提出了雙相延遲模型(Dual Phase Lag Model)來描述微小時間內的微觀熱傳遞現象。此理論假設材料受熱後溫度梯度場與熱通量及熱傳遞分別在不同時間形成，熱傳遞在(t)時間開始，經過τT秒(時間 t+τT)及τq秒(時

間t+τq)後分別形成溫度梯度場及熱通量，其中τT及τq分別為溫度梯度場及熱通量場形成的延遲時間(Time Delay)。Tzou 經由此假設推導描述微觀熱傳遞的統御方程式，該方程式可得到與描述金屬內微觀熱傳的兩階段加熱模型以及描述非金屬的微觀熱傳理論有相同的結果。同時當不考慮溫度梯度的時間延遲下(τT=0)，雙相延遲模型可得到波動熱傳理論(C.V. Wave Model)[20,21]，而當不考慮熱通量及溫度梯度的時間延遲下(τT=τq=0)，雙相延遲模型可推導出傅立葉熱傳方程式。

(17)

1.3 研究方法

由上面的文獻回顧可知，對於皮秒及飛秒級的脈衝雷射加熱過程，其熱傳遞現象目前仍處於一維傳遞的探討。本論文針對材料受極短脈衝雷射加熱過程中二維熱傳現象進行詳細的研究，經由雙相延遲模型建立適當的統御方程式，使用拉式轉換(Laplace Transform)及 Hankel 轉換將統御方程式轉成常微分形式來求解，再配合黎曼和近似法(Riemann-Sum Approximation Method) 與高斯格點法 (Gaussian Quadrature Method)進行 Fortran 程式編寫，以求得拉式逆轉換及 Hankel 逆轉換的數值解。其求解過程將詳述於第三章。而第四章則主要以變更溫度場與熱通量場延遲時間之比值(B 值)，連結三種熱傳遞模型，即傅立葉模型、熱波模型與雙相延遲模型來針對所得之數值模擬結果進行分析探討。

(18)

第二章巨觀與微觀熱傳簡介

由微觀的角度而言，金屬的熱傳主要藉由聲子(Phonons)與電子的交相反應來傳遞能量。至於介電體、絕緣體和半導體，則單純由聲子與聲子間的碰撞所造成。傳統的熱傳理論,即傅立葉定理(Fourier’s Law)，所呈現的是巨觀下的行為，數以千萬計的聲子有著足夠的時間來達成所需的平衡。相較於時間等級在皮秒甚至飛秒下的熱傳行為，

則勢必需要其他的熱傳模型來加以描述。本節針對數種已被提出的熱傳模型來加以說明，介紹其不同基本假設及數學模型。

2.1

^{傳統傅立葉熱傳模型}

當一物體內有溫度梯度(Temperature Gradient)存在時，經驗顯示能量會從高溫區傳到低溫區。其中，熱傳導率(Heat Transfer Rate)和法線溫度梯度(Normal Temperature Gradient)成正比，代入比例常數後表示如下：

q= − ∇ (2.1) k T 而從能量守恆定理可知：

t C T

q _P

∂

= ∂

•

∇

− (2.2) 上式中 q 是熱傳遞速率，C 為定壓比熱， T_p ∇ 是熱流流動方向的溫度 梯度。正值常數k稱為材料的熱傳遞係數 (Thermal Conductivity)，加

(19)

上負號是為了滿足熱力學第二定律；亦即熱必須由高溫傳至低溫。方程式(2.1)稱為熱傳導的傅立葉定律。

上述假設配合能量方程式作一階近似展開後，傳統傅立葉熱傳方程式統御方程式即可表示成：

t T T

∂

= ∂

∇ α

2 1 (2.3)

根據傳統傅立葉熱傳導理論可得知，在擴散模型(Diffusion Model) 中，假設在固體中之熱傳遞速度為無窮大，因此使邊界條件或起始條件瞬間作用於物體中，無任何延遲時間，但熱傳遞速度為無窮大已不符合實際之物理現象，運用在微觀熱傳上則必須加以修正。

2.2

聲子-電子交相反應模型

金屬以聲子和自由電子來進行熱傳導，但是以自由電子為主，因為金屬中有大量的自由電子。其熱傳行為是互相影響的，故本節所介紹的聲子-電子交相反應模型(Phonon-Electron Interaction Model)是利用兩階段的過程來描述金屬的微觀熱傳行為，此模型又稱之兩階段加熱模型(Two Step Heating Model)。電子受熱後溫度升高，經由聲子- 電子偶合因子(Phonon-Electron Coupling Factor)加熱聲子，數學化此過程如下：

(20)

( ) ( )

e

e e e l

C T K T G T T

t

∂ = ∇ ∇ − −

∂ i (2.4) ( )

l

l e l

C T G T T t

∂ = −

∂ (2.5) 其中，C 為熱容，K 為電子的熱傳導係數，下標 e 與 l 則分別代表電子及聲子，G 則是聲子-電子偶合因子。表[1]列出數種金屬材質的 G 值：

【表1】金屬之聲子-電子偶合因子[19]

(2.5)式忽略了聲子的熱傳效應。從數學的觀點來看，(2.4)與(2.5)式裡有兩個未知數，即聲子熱容(C_l)與電子熱容(C_e)。利用(2.5)式聲子與電子的關係，將之代入(2.4)式，則可得到分別就聲子與電子角度下描述之兩階段加熱模型統御方程式。表示如下：

2

2 2

2

l l e l l e l

l l

C C C T C C T

T T

G t K t KG t

∂ + ∂ ∂

⎛ ⎞ ⎛ ⎞ ⎛ ⎞

∇ +⎜⎝ ⎟⎠∂ ∇ =⎜⎝ ⎟⎠ ∂ +⎜⎝ ⎟⎠ ∂ (2.6)

2

2 2

2

l l e e l e e

e e

C C C T C C T

T T

G t K t KG t

∂ + ∂ ∂

⎛ ⎞ ⎛ ⎞ ⎛ ⎞

∇ +⎜⎝ ⎟⎠∂ ∇ =⎜⎝ ⎟⎠ ∂ +⎜⎝ ⎟⎠ ∂ (2.7)

(21)

2.3

聲子散射模型

對於非金屬的材料而言，因沒有自由電子，故熱傳主要以聲子間的碰撞為主。而聲子散射模型(Phonon-Scattering Model)藉由解決波茲曼傳遞方程式，忽略其中電子傳遞效應的貢獻，可得到其統御方程式如下：

p 0

C T q

t

∂ + ∇ =

∂ i (2.8)

2 2

( )

1 2

3 5 2

p N

R

q c C c

T q q q

t

τ τ

∂∂ + ∇ + = ⎡⎣∇ + ∇ ∇i ⎤⎦ (2.9)

c 平均聲子速，τ_R為Umklapp Process 的遲緩時間(Relaxation Time)，τ_N 則為Normal Process 的遲緩時間。

對於非金屬的熱傳，高溫區的聲子數量與能量密度都高於低溫區，因此整個熱傳過程皆為聲子逐步由高溫區擴散至低溫區。由於聲子的流動主導整個熱傳過程，故當聲子碰撞後的速度方向與碰撞前的方向相反，便將其稱之為 Umklapp Process；反之，則稱為 Normal Process。Umklapp Process 雖然是動量不守衡的狀態，然其存在表達了一個重要事實，即熱傳導係數不會無限制的擴大。假如不存在此機制，即便沒有溫度梯度造成聲子流動，能量的傳遞依舊會保持定值。

連結(2.9)與(2.10)後，可得到聲子散射模型的統御方程式：

( )

²

2 9τ ∂ 2 3 ∂T 3 ∂ T

∇ + ∇ = + (2.10)

(22)

2.4

熱波模型

熱波模型(Thermal Wave Model)以傅立葉定律為基礎，提出熱通量與溫度梯度間存在時間延遲，而非瞬間傳遞。以聲子碰撞之延遲時間為背景，將熱傳現象中之熱通量與物體內之溫度梯度視為非同時發生，此即為所謂的CV Wave Model。將熱通量之時間項加入一延遲時間τ ，則方程式表示如下：

T t k

q q =− ∇

∂

+τ ∂ (2.11)

其中τ 為延遲時間。配合能量方程式展開可得熱波模型之統御方程式：

2 2 2

2 1 1

t T t C

T T

∂ + ∂

∂

= ∂

∇ α (2.12) 其中α 為熱擴散係數(

CP

= k

α )、C為熱波速度(

τ

= α

C )。溫度對時間

二次微分使熱傳遞形式從擴散變成以波的形式做熱傳遞。從上式可知傳統傅立葉熱傳模型暗示熱波速度為無限大。然相對微觀的熱傳行為，傅立葉熱傳模型的假設顯然不合理。熱波模型僅假設熱通量的延

遲，Tzou[19]則進一步提出，溫度場、熱通量與溫度梯度三者皆有延

遲時間存在，從而提出雙相延遲模型。

(23)

2.5

雙相延遲模型

延遲行為指出，熱通量與溫度梯度並非同時產生，此間關係互成

因果，這明顯地違背了傳統的熱傳假設。本節比較各種不同熱傳模型，以延遲行為將之連貫。

雙相延遲模型的理論基礎在於假設熱傳遞、溫度梯度與熱通量在

不同時間發生。當熱傳遞於時間(t)形成時，經過τT _與τ_q_{秒後分別形} 成溫度梯度場(t+τT ₎_{與熱通量場}_(t+τ_q₎_，其中

τT與τ_q分別為溫度場與熱通量場的延遲時間(Time Delay)，其表示如下：

(

r t q

)

k T(r t T)

q , +τ =− ∇ , +τ (2.13) 將上述方程式展開後可得：

( ) ( ) ( )

[

( )

]

⎭⎬

⎫

⎩⎨

⎧ ∇

∂ + ∂

∇

−

∂ ≅

+ ∂ T r t

t t r T k t t r t q

r

q , τ_q , , τ_T , (2.14) 再對方程式(2.15)取Divergence 後可得：

[ ] [ ]

^T

k t T k t q

q _q ² _T ∇²

∂

− ∂

∇

−

=

•

∂ ∇ + ∂

•

∇ τ τ (2.15) 配合能量方程式展開便可得到統御方程式如下：

[ ]

² ²₂

2 1 1

t T t

T t

Q Q T k

T _T t _q ^q

∂ + ∂

∂

= ∂

⎥⎦⎤

⎢⎣⎡

∂ + ∂ +

∂ ∇ + ∂

∇ α

τ τ α

τ (2.16)

式(2.16)與傳統傅立葉模型比較後，其中除了有波形項

2

q T

t τ

α

∂

∂ 外，還有溫度梯度造成的

[ ]

^T

T t

∇2

∂

τ ∂ ，此項會增強熱傳造成短時間內熱

傳遞影響的空間範圍增大。

(24)

將式(2.16)與(2.10)式的聲子散射模型或是(2.6)式與(2.7)式的兩階段加熱模型相比可知，倘若不考慮其係數，彼此仍有著相同的形

式。故只要藉由變動兩個延遲時間，便可使雙相延遲模型的統御方程式吻合其他熱傳模型的假設。若假設τT ₌τ_q₌₀_，則_(2.16)_{可簡化成傳} 統傅立業熱傳模型；忽略溫度梯度延遲時間，τT ₌₀_{，則簡化為熱波}

模型；若要化成兩階段加熱模型，則需使用下列假設：

e l

k

C C

α =

+ ^，

l T

C τ = G ，

1 1 1 1 q

e l

G C C

τ

⎡ ⎤−

= ⎢ + ⎥

⎣ ⎦ (2.17) 同樣的，以下列定義，則可清楚地模擬聲子散射模型。

2

3

Rc

α =τ ， 9

5

N T

τ = τ ，τ_q =τ_R (2.18)

如此，根據不同的延遲時間假設，雙相延遲模型便可完美地連結前述的所有熱傳模型，提供了新型態的熱傳方程式。本論文採用雙相延遲模型，並在對統御方程式無因次化後，藉改變兩個延遲時間比值

2

T

q

τ

τ ^，^即^B 值，來模擬傳統傅立業熱傳模型(Diffusion Model)，熱波(CV

Wave Model)，並比較此兩種熱傳模型與雙相延遲模型之差異處。

(25)

第三章脈衝雷射加熱薄膜理論模式

本論文針對脈衝雷射加熱過程之二維熱傳遞進行研究。現敘述過程如下：

3.1 熱傳模型與熱源項

考慮一半無窮平板，受到一脈衝雷射光加熱，在此採用圓柱座

標，其物理模式及座標系統如【圖一】所示，

假設材料為等向性(Isotropic)，材料受熱後的溫度分布為軸對稱，且滿足下面微觀熱傳方程式[13]：

2 2 2

2

1 1

[ ] [ ] ^q

T q

Q T T

T T Q

t k t t t

τ τ τ

α α

∂ ∂ ∂ ∂

∇ + ∇ + + = +

∂ ∂ ∂ ∂ (3.1)

方程式中 T 代表加熱後溫度的上升量、α為材料熱擴散率、

【圖1】物理模式及座標圖

(26)

2 2

2

1 (r )

r r r z

∂ ∂ ∂

∇ = +

∂ ∂ ∂ ，而Q為雷射光強度，表示如下：

2

0

( , , ) ( ) ( _p) exp

p

z r

Q r z t Q H t H t t

δ r

⎡ ⎛ ⎞ ⎤

⎡ ⎤ ⎢ ⎥

= ⎣ − − ⎦ ⎢⎣− − ⎜ ⎟⎝ ⎠ ⎥⎦

(3.2a)

0

0.94 1

p p

Q J R

t δ

⎛ − ⎞

= ⎜⎜⎝ ⎟⎟⎠ (3.2b)

其中Q₀為雷射吸收強度(Intensity of Laser Absorption),J為雷射能量 (Laser Fluence)， t_p 為脈衝時程， δ_p 為輻射穿透深度(Radiation Penetration Depth)。

方程式(3.1)比傅立葉熱傳方程式(2.3)多出幾項。 ^q ²^T₂ t τ

α

∂

∂ 為熱通量

延遲時間造成，此項的效應使熱由擴散形式變成熱波傳遞，熱波的傳播速度為有限值。另外，熱通量延遲時間讓方程式中除了正常的雷射

熱源項(Q(r,z,t))外，還增加了一額外熱源項 ^q ^Q

k t τ _∂

∂ 。 _T [ ²T] τ ^∂t ∇

∂ 為溫度梯

度延遲造成，為方程式中的最高次項，溫度對空間為二次微分，對時間為一次微分，此會增強熱傳讓熱影響範圍在短時間內增大。

假設材料在室溫下被加熱且無其他能量施加，薄膜厚度L，則起

始條件可寫成：

(27)

T = at T0 t =0,0≤ < ∞ ≤ ≤r ,0 z L (3.3a) T 0

t

∂ =

∂ at t =0,0≤ < ∞ ≤ ≤ (3.3b) r ,0 z L

由於加熱雷射脈衝非常短，所要考慮的暫態時間非常短，因此假設材料的表面及底部為絕熱狀態，而無窮遠處不受雷射影響。故邊界條件可寫成：

T 0 z

∂ =

∂ at z=0 or L , t> (3.4a) 0

0 at r=0, T T at r=0

T r

∂ = → ∞

∂ (3.4b)

3.2 積分轉換法求解

方程式(3.1)無因次化後，表示如下：

( )

²

2 2

4 2 2 2

B η θ θ

θ θ η

β β β β

⎛ ⎞

∂ ∂ ∂ ∂

∇ + ∂ ∇ +⎜⎝ + ∂ ⎟⎠= ∂ +∂ (3.5)

其中，無因次化變數為：

0 0

T T

θ = T⁻ ，η =Q/

(

kT0/ατ_q

)

，

2 _q β t

= τ ，

2 _q

δ z

= ατ ，

2 _q

γ r

= ατ ，

2

T

q

B τ

= τ

而雷射熱源Q 無因次化後可得：

(28)

( )

²

0

0 0

0

/ _q exp ^p _p

Q H H

kT

δ γ

η β β β

δ γ

ατ

⎡ ⎛ ⎞ ⎤

⎡ ⎤ ⎢ ⎥

=⎡⎣ ⎤⎦ ⎣ − − ⎦⎢⎣− − ⎜ ⎟⎝ ⎠ ⎥⎦

(3.6)

式中無因次化變數為：

0 2

p p

q

δ δ

= ατ ，

2

p p

q

β t

= τ ， ₀ ⁰

2 _q

γ r

= ατ

整理得：

0 r z t

η η η η η= (3.7)

其中 ₀ ⁰

0/ _q

Q η kT

= ατ

⎡ ⎤

⎣ ⎦ ，^η^t ⁼^⎡_⎣^H( )^β ⁻^H

(

^{β β}⁻ ^p

)

^⎤_⎦^，

0 z exp

p

η δ

δ

⎡ ⎤

= ⎢− ⎥

⎢ ⎥

⎣ ⎦，

2

0 r exp

η γ

γ

⎛ ⎞

= − ⎜ ⎟

⎝ ⎠

無因次化起始條件可寫成：

θ θ= 0 at β =0, 0≤ < ∞ ≤ ≤γ , 0 δ (3.8a) θ 0

β

∂ =

∂ at β =0, 0≤ < ∞ ≤ ≤γ , 0 δ (3.8b)

無因次化邊界條件可寫成：

θ 0 δ

∂ =

∂ at δ =0 or ,β >0 (3.9a)

0 at =0, 0 at =

θ γ θ θ γ

γ

∂ = → ∞

∂ (3.9b)

(29)

其中， 2 _q L

= ατ

方程式(3.5)表示無因次化的極短脈衝雷射加熱過程之二維暫態熱傳遞方程式，而方程式(3.6)則是無因次化的脈衝雷射強度，配合起始條件方程式(3.8)及邊界條件方程式(3.9)，可解出其溫度分佈。論文中對於空間方向γ 微分採用Hankel 轉換(Hankel Transfom)[22]來解決，就本人所知，相關探討的問題大都為一維問題，忽略了其他方向的傳遞，因此目前仍無人使用此方法來解方程式(3.5)，但就方程式(3.5) 的型態來看，含有空間微分的兩項皆為∇²θ的形式，因此利用Hankel 轉換來處理γ 方向微分便可以順利進行。定義Hankel Transform 如下：

( , , )

[ ]

0 ( , , ) ( )J0 d θ σ δ β = Η θ =

∫

^∞θ γ δ β γ σγ γ

在此，σ 為轉換參數。則式(3.5)轉換如下：

2 2 2

2 2

2 B 2 4 2 2 2

θ σ θ θ σ θ η η θ θ

δ β δ β β β

⎡∂ − ⎤+ ∂ ∂⎡ − ⎤ ⎡+ + ∂ ⎤= ∂ +∂

⎢∂ ⎥ ∂ ⎢∂ ⎥ ⎢ ∂ ⎥ ∂ ∂

⎣ ⎦ ⎣ ⎦ ⎣ ⎦ (3.10)

而方程式(3.6)轉換如下：

0 r z t

η η η η η= (3.11)

其中，

2 0

2 2 2

0 0

0 0( ) exp

2 4

r e J d

γ

γ γ σ γ

η γ σγ γ

⎛ ⎞

∞ −⎜ ⎟⎝ ⎠ ⎛ ⎞

= = ⎜− ⎟

⎝ ⎠

∫

(30)

對於時間微分則採用拉式轉換(Laplace Transform)來解決，Laplace Transform 定義如下：

( ^{, ,}^s) ^L ₀ ( ^{, ,} )^e ^s^β^d

θ σ δ = ^{⎡ ⎤}⎣ ⎦θ =

∫

^∞θ σ δ β ⁻ β

在此，s 為轉換參數。則式(3.10)與(3.11)經轉換整理如下：

2

d c f

d

θ θ

δ ⁻ = (3.12)

其中

(

² ² ² ²

)

_, (^{4 2} )

1 1

BS S S S

c f

BS BS

σ ⁺ σ ⁺ ⁺ − + η

= =

+ +

0 r z t

η η η η η= (3.13)

其中 ¹ ^s

s t

e s η _{= ⎜}^⎛ ⁻ ⁻β ^⎞_⎟

⎝ ⎠

經過拉式轉換及Hankel 轉換處理後，方程式(3.5)可轉換成對δ 微分的常微分方程式，配合邊界條件可將此常微分方程式解出，得到θ經過轉換的表示式θ σ δ( , , )s 。

0

1 2 3

F F p

A e A e A e

δ

δ δ δ

θ = + ⁻ + ⁻ (3.14)

其中，F = c

(31)

而

0

3 1

0 F p

F F

p

A e e

A Fe Fe

δ

− −

−

⎛ ⎞

⎜− + ⎟

= ⎜⎜ − ⎟⎟

⎝ ⎠

，

1 3

0 2

p

A F A

A F

−δ

=

( )

⁽ ⁾

3 0 2

0

4 2 _t _r / 1 1

p

A S η η η BS c

δ

⎡ ⎛ ⎞⎤

= − + ⎢⎢⎣ + ⎜⎜⎝ − ⎟⎟⎠⎥⎥⎦

3.3 數值方法

(3.14)式是一複雜的數學形式，無法藉由逆轉換表得到解析解 (Exact solution)，因此本論文以數值方法來求解逆轉換的結果。關於拉式逆轉換的數值方法，本文以黎曼和近似法 (Riemann-Sum Approximation Method)轉換， Hankel 轉換的逆轉換，其表示式是一積分式，積分範圍為零到無窮大，可以用數值方法來解此積分，本論文以高斯格點法(Gaussian Quadrature Method)來進行此數值積分，從而得到二維暫態溫度分布。

3.3.1 Laplace Transform 逆運算-黎曼和

Laplace 逆轉換定義如下：

( ^{, ,} ) ¹ ( ^{, ,} )

2

i s

i s e ds

i

γ β

θ σ δ β γ θ σ δ π

+ ∞

=

∫

− ∞ (3.15) 此形式為一複數積分，定義式中的Laplace 轉換參數 s 如下：

(32)

則(3.15)可改寫為：

( ^{, ,} ) ( ^{, ,} )

2 e i

s i e d

υβ ωβ

θ σ δ β θ σ δ υ ω ω

π

∞

=

∫

−∞ = + (3.17) 積分域由負無窮大到正無窮大，配合黎曼和近似法，定義

ω

為頻率，

τ

為半週期時程，即ω ⁿ^π

= τ ，將(3.17)式展開成下式：

( ^{, ,} ) ^{, ,}

2

i n

e in

s e

υβ βπ

π τ

θ σ δ β θ σ δ υ

τ τ

⎛ ⎞

∞ ⎜⎝ ⎟⎠

−∞

⎛ ⎞

=

∑

⎜⎝ = + ⎟⎠ (3.18) 此波形正負值對稱出現，且有著下式的關係，

( ) ( )

, , in ⁱⁿ , , in ⁱⁿ

e e

βπ βπ

τ τ

π π

θ σ δ υ θ σ δ υ

τ τ

⎛ + ⎞ + ⎛ − ⎞ −

⎜ ⎟ ⎜ ⎟

⎝ ⎠ ⎝ ⎠

_{2 Re} _{, ,} ⁱⁿ ⁱ( )ⁿ

e

βπτ

θ σ δ υ π τ

⎡ ⎛ ⎞ ⎤

= ⎢⎣ ⎜⎝ + ⎟⎠ ⎥⎦ (3.19) 所以，(3.18)式可改寫為下式：

( ) ( ) ( )

1

, , 1 , , Re , ,

2

in

n

e in

e

υβ βπ

π τ

θ σ δ β θ σ δ υ θ σ δ υ

τ τ

∞

=

⎡ ⎛ ⎞ ⎤

= ⎢⎣ +

∑

⎜⎝ + ⎟⎠ ⎥⎦ (3.20)

其中，Re 為和(Summation)的實數部。因為函數 ⁱ( )ⁿ

e

βπτ 有固定的週期時

程2π ，故

β

必定落於0≤ ≤β 2τ 的關係。令β τ= ，進一步將(3.20)式改寫如下：

( ) ( ) ( )

1

, , 1 , , Re , , 1

2

N n

n

e^υβ inπ

θ σ δ β θ σ δ υ θ σ δ υ

β ₌ τ

⎡ ⎛ ⎞ ⎤

≅ ⎢⎣ +

∑

⎜⎝ + ⎟⎠ − ⎥⎦ (3.21)

如此一來，上述的式子便可進行數值運算。然而，所得之結果仍需加以進行Hankel Transform 逆轉換，本文採取高斯格點法進行，以下對使用方法進行介紹。

(33)

3.3.2 Hankel Transform 逆運算-高斯格點法

若函數 ( , , )θ γ δ β 的 Hankel Transform 為 ( , , )θ σ δ β ，則根據定義 ( , , )

θ γ δ β 可表示為：

(

, ,

)

H ¹

(

, ,

)

0 J0

( )

d

θ γ δ β = ⁻ ⎡⎣θ σ δ β ⎤⎦=

∫

^∞θσ σγ σ (3.22) 對一任意函數，高斯格點法(Gaussian Quadrature)選擇了[a,b]間

x1，x₂，…，x_n n 個點以及c₁，c₂，…，c_n n 個常數使得下列積分近似之誤差為最少：

1

( ) ⁿ ( )

b

i i

a i

f x c f x

=

≈

∑

∫

(3.23)

對於任一區間均可用簡單的線性轉換^t⁼^⎡_⎣¹ (^{b a}⁻ )^⎤_⎦⁽²^{x a b}^{− −} ⁾^對應至

[-1,1]，故可得到下式之結果：

( ) ( ) ( )

1

( ) 1

2 2

b

a

b a t b a b a

f x dx f dt

−

− + + −

⎛ ⎞

= ⎜ ⎟

⎝ ⎠

∫ ∫

(3.24)

假設x₁，x₂，…，x_n為第 n 次的Legendre 多項式 p_n之根，則其係數c_i 定義如下：

1 1 1

n

j i

j i j

j i

x x

c dx

x x

− =

≠

= −

∏

−

∫

(3.25)

如此，可得到：

( ) ( )

1

1 1

n

i i

i

P x dx c P x

− =

=

∑

∫

(3.26)

(34)

【表2】Legendre 多項式之根及係數[23]

n Legendre Polynomials Roots x_i Coefficients c_i

1 2 3

4

5

這些根及係數大多已有列表[24]，表[2]列出數個 Legendre 多項式之根及係數關係。在此使用高斯格點積分，繼續進行Hankel 逆轉換。

由於(3.22)式積分域為零到無窮大，故本文定義一值以代替無窮大，

經反覆測試，確認500 已足夠大到不影響數值結果，並將積分區間 500 分成 50 等份，格點取 96 點。

關於零階Bessel Function，其形式也是一個積分式，故本論文使用高斯格點積分。由文獻[25]，可表示零階 Bessel Function 如下：

( )

^sin

( )

0 0

1 1

cos sin 2

J σγ ^π_πei^σγ ^θdθ ^π σγ θ θd

π ⁻ π

=

∫

=

∫

− (3.26)

(35)

經由高斯格點積分，轉換為：

( )

¹

0 1

1 cos sin

2 2

J σγ =

∫

₋ ^⎡⎢⎣−σγ ^⎛⎜⎝πt⁺π ^⎞⎟⎠^⎤⎥⎦dt

_, ^,

1

1 cos sin

2 2

i

n i n i

k

C σγ πR π

=

⎡ ⎛ + ⎞⎤

= ⎢− ⎜ ⎟⎥

⎝ ⎠

⎣ ⎦

∑

(3.27)

(36)

第四章結果與討論

利用上一章之數值方法模擬後，藉著改變溫度場與熱通量場的延遲時間之比值B來串連三種熱傳模型，分別為B=0 的熱波模型、

B=0.5的傳統傅立葉模型與 B等於其他值的雙相延遲模型。以下對

結果進行分析。本文先就所得結果與文獻Tzou [19]做比較，然後分別就徑向、軸向與固定點的時間變化加以討論。徑向比對旨在說明三種熱傳模型的差異，軸向分析則進一步探討 B值與二維效應產生的相對誤差為何。定點分析則顯示時間改變下的溫度分佈。

4-1 程式驗證

由於Tzou[19]利用雙相延遲模型模擬雷射加熱薄膜之一維熱傳

遞，經反覆測試，γ₀ =80已與一維熱傳情形所差無幾，故在程式中以極大雷射半徑γ₀ =80來模擬一維熱傳情形。並由於文獻未考慮穿透深度，所以在此假設無因次穿透深度極小δ_P₀ =0.001。其他參數設定分別為無因次雷射能量η₀ =500、薄膜與δ 軸深度δ = =1.0與雷射脈衝時程β_p =0.1。由四個不同 B值模擬三種熱傳行為。其對照結果如【圖 2】所示。由【圖2】可知與文獻之結果相當吻合，可以清楚看出雷射熱源穿過薄膜時，三種本論文所探討的熱傳遞模型彼此間之差異。

除了四個B值極度吻合外，連熱波模型碰到薄膜底板所產生的反彈

(37)

熱波也清楚的呈現其中。

【圖2】溫度分佈與文獻[19]比較結果

4-2

徑向溫度分佈

本節藉由找出溫度在徑向的分佈情形，比較三種模型間的異同，

其無因次參數為雷射能量η₀ =50，穿透深度δ_P₀ =0.001，雷射脈衝時程β_P=0.1以及薄膜厚度與δ 軸深度 =δ =1。【圖3】到【圖6】呈現了薄膜底面沿徑向的溫度分佈。四張圖分別預測了無因次雷射半徑γ₀=1 的情況下，6個不同B值(0，0.2，0.5，1.0，5.0，10.0)，在四個不同無因次時間β值(0.01，0.1，0.5，1.0)下的溫度分佈。

(38)

【圖3】γ₀=1.0，β =0.01，B值改變下溫度分佈圖

【圖 4】γ₀=1.0，β =0.1，B 值改變下溫度分佈圖

(39)

【圖 6】γ =1.0，β =1.0，B 值改變下溫度分佈圖

(40)

由圖中可看出，熱波(B=0)在時間β為0.01，0.1 與0.5，即【圖3】到【圖5】中，溫度都還沒上昇。反之，在時間為1.0的【圖 6】中則上昇的比其他模型還高出許多。這是因為預測處為薄膜底部溫度，

在前三個時間下，熱波的能量尚未來得及傳遞至底部，而在時間為 1.0 時，熱波能量傳達到底面，造成溫度上升。此情形吻合了熱波的基本假設，即能量的傳遞速度並非無窮快。上述的溫度分佈之預測，

對於極快速雷射加工的精度控制，必定有著極大誤差產生。

進一步探討B=0.5 的傅立葉熱傳模型，以及B=0.2，1.0，5.0 與 10.0的雙相延遲模形。可以發現在【圖3】到【圖6】下，B值不甚大的情況所預測的溫度差距並不會有太大的差別。然而隨B值的漸增，其差異也隨之逐步擴大，且溫度分佈的影響範圍也因而更廣。尤其對B 值較大的材料所產生的影響極巨。由表[3]可知，銅的B值約

【表3】四種金屬之 B值[19]

為81.5、而金則約為60，其他不同金屬材料的B 值也遠大於傅立葉模型的0.5 或熱波的0.0。如果以傅立葉熱傳模型或是熱波模型來進行暫態的溫度分佈模擬，則預測之溫度必然相差甚遠。

(41)

【圖 7】到【圖9】取與【圖3】到【圖 6】的相同參數設定下，

在β =0.5時，三個逐漸遞增的雷射半徑(γ₀=3.0，10.0，20.0)。

【圖8】γ =10.0，β=0.5，B值改變下溫度分佈圖

(42)

【圖9】γ₀=20.0，β=0.5，B值改變下溫度分佈圖

觀察其溫度分佈，幾乎可得到上述相同結論，即隨B值增加，

則所預測溫度與傅立葉熱傳模型之差距也隨之增大。

4-3

軸向溫度分佈

本節針對二維效應對薄膜溫度分佈之影響，討論不同雷射半徑γ₀

及B值之影響。在此先模擬熱波模型在不同的時間點所預測的溫度分佈，即β =0.1、0.5、1.0、1.5，其他無因次參數的設定分別為雷射能量50.0、穿透深度0.001 與七個雷射半徑來分析一維到二維間之差異。【圖10】到【圖 13】所討論的是B為 0.0的 CV Wave，從圖中可

(43)

清楚的看見熱波往前傳遞，以及碰到邊界後反彈的情況。

【圖 10】β=0.1，雷射半徑改變下溫度分佈圖

【圖11】β =0.5，雷射半徑改變下溫度分佈圖

(44)

(45)

從【圖10】到【圖 13】中亦可以看出，隨著時間改變，雷射半徑的改變造成的溫度預測差距越益明顯。在無因次雷射半徑5以上，

當半徑增加，所得的溫度分佈並沒太大改變。而在無因次雷射半徑5 以下的溫度分佈，隨著雷射半徑愈小，差距愈來愈大。由此可知，溫度隨著雷射半徑的減少而遞減。

【圖14】到【圖 16】則模擬另外兩種熱傳遞模型，即當B=0.5 時的傅立葉熱傳模型與B=0.2 和B=1.0時的雙相延遲模型。參數如同上述熱波模型所設定，無因次時間定為β =0.5。從圖中可看出，不論何種模型，當無因次雷射半徑大於5以上的溫度分佈幾乎沒太大差別。因此可推知在此參數設定下，二維效應在小於5的無因次雷射半徑下影響愈來愈大。

β

中 華 大 學 碩 士 論 文