2.1 2.1 熱傳導方程式 熱傳導方程式 熱傳導方程式 熱傳導方程式

第二章 第二章 第二章 第二章

熱傳導的數值計算方法 熱傳導的數值計算方法 熱傳導的數值計算方法 熱傳導的數值計算方法

在這一章節裡，我們將描述有限差分方法[1](Finite Difference Method，

簡稱，FDM)來計算熱流模型，並且介紹基本熱傳導方程式架構以及數值方法 (numerical method)。

熱傳學是研究熱能的傳遞方式及預測在特定情況下，能量傳遞率的科學；也 是處理熱能、動能及質量傳遞的科學。在熱力學中，由於系統和外界之間溫度差 所引起的通過邊界之能量稱為熱；熱是與過程有關的，本身無法直接測量或觀 察，但熱所產生的效果是可以測量並且觀察的。因此我們利用熱傳導是討論因溫 度的改變所引起的能量轉換和傳遞，且可以了解熱能與時間和位置的關係，並且 求得系統的溫度分佈。光碟的讀取、寫入以及擦拭的情況，都是利用雷射光照射 到光碟上所產生的能量與熱的轉換來達成，因此可以利用熱傳導來計算空間中不 同的時間裡溫度的分佈，進而調整雷射光功率來達到讀取、寫入、擦拭的情況。

因此光與熱的交互作用後所產生的溫度分佈，就可以利用熱傳導的計算來求得。

熱的應用對於相變化光碟是非常重要的，相變型光碟材料非結晶相的熱穩定 性相當重要，以必免在高溫的氣候下讓資料遺失，因此在研究相變型光碟時若能 有效控制雷射光所提供的能量使材料溫度可以提升的多寡，就可以有效控制材料 的熱穩定性或是藉由熱產生的溫度高低可以適當的添加一些別種元素來控制材 料的融點或是升溫的速度，因此這些溫度上的量測都需要利用熱傳導學來加以研 究才能真正達到有效的成果。

2.1 2.1

2.1 2.1 熱傳導方程式 熱傳導方程式 熱傳導方程式 熱傳導方程式

當一個雷射光聚焦在光碟片上時部份的入射光會被記錄層(dye/phase change)和其他的吸收層(metallic)所吸收，這些被吸收的光能會被轉換成熱, 而熱就由 Q(x,y,z,t)[W/cm³]來表示熱源，熱源是位置與時間的函數。溫度所滿 足的熱擴散方程式[2]為：

) , , , ( ]

[ T Q x y z t t

C

_p

T =∇• ∇ +

∂

∂

κ

ρ

(2.1)

] / / [W m °C

κ

為熱傳導係數(heat conductivity)，為物質的傳輸性質，是熱以傳 導方式通過材料中快或慢的指標，

ρ

[Kg/m

³

]為材料的密度(density of the material)，其物理意義為物體內含物質的疏密程度，C

_p

[J/Kg/° 為比熱C]

(specific heat)，其物理意義為每單位質量升高一度所需要的能量。在光碟片 中，隨著每一層材料的不同因此各有不同的材料係數。

2.2 2.2

2.2 2.2 數值方法 數值方法 數值方法 數值方法

當解析的方法在計算上很困難或是根本不存在的時候，我們就可以利用一些 數值方法來求得我們要的解，在這一小節裡我們將介紹一些數值方法並且利用有 限差分法來計算熱傳導方程式。數值分法常用來解析複雜幾何形狀、非均勻邊界 條件，或與溫度有關之邊界和性質變化問題。對於偏微分方程的數值分析解大略 可以分為有限元素法與有限差分法兩種。在工程問題的處理上利用有限元素法 (Finite Element Method，簡稱，FEM)是一種實際而有效的方法，尤其在力學問 題的處理上更有其獨到之處，解決問題的重點是在於元素形狀的選取，此法也比 有限差分法來得複雜。因此我們所採用的數值計算方法為有限差分法。

2.2.1 2.2.1 2.2.1

2.2.1 有限元素法 有限元素法 有限元素法 有限元素法

有限元素法[3]的基本觀念是，任何連續量(continuous quantity)，像是 溫度、壓力、或位移等，均可用一不連續函數的型式作近似的表示。此型式乃為 有限區域的集合分段連續函數所組成。因此將所要計算的空間分割成許多不互相 重疊之次空間(元素)，並將元素內之物理量用節點之物理量內差代表之。由於有 限元素法基本上為非結構性網格(unstructured grid)，因此適合幾何形狀複雜 的物件或計算空間。隨著元素數目的增加，以及元素尺寸的減小，將會得到更精 確的結果，如果元素本身滿足收斂性要求，近似解也將會是個收斂的結果。

舉例而言，當遇到一個材料性質不均勻的物質時，可以將該材料作有限度的 分割，分割出有限的小單元，由於單元很小可以假設每一小單元內部其性質是相 同的，而單元與單元的性質不同，這樣就化繁為簡，可以進行分析，元素切割越 小，所得結果越精細，此稱為有限元素分析。有限元素法建立在不連續的型式上，

在區域內的有限點均屬同一性質，這些點稱為節點(node)，每個節點連續量的值 表 示 所 需 求 出 的 一 個 變 數 。 由 區 域 分 割 而 成 的 有 限 個 小 區 域 稱 為 元 素 (element)，將這些元素的共同節點連接起來並集合估計區域的形狀。

有限元素法的優點為鄰近元素的材料性質不一定要相同。如此，能使這種方 法應用到由不同材料所組成物體上，且在不規則形狀的邊界，能用直線的元素作 近似或用曲線邊界作正確的配合。而缺點為運算複雜,必須要使用大量的記憶裝 置，而比起有限差分法較不直觀。

有限元素法與有限差分法最大的差別在於後者是以微分的型態去得到解，而 前者卻是以積分的型態來解決問題。有限元素法中最重要的數學觀念是變分法 (variational calculus)觀念，這方法中要將一積分式之類得函數給予極小化。

以下列的式子為例:

∫ =

=

L

x

dx x u x u x F I

0

' ( ))

),

(

,

(

在上式中，需變動 u(x)使積分 I 達到最小值，這是變分法最主要的目的。

2.2.2 2.2.2 2.2.2

2.2.2 有限差分法 有限差分法 有限差分法 有限差分法

有限差分法是用來求解(2.1)式的數值解，也就是把連續的偏微分方程式 (PDE)離散化，進而離散(discrete)近似成我們所要求得的數值解，這裡離散是 指用有限的點數來求得數值解，一般來說增加空間離散點數會同時增加數值解的 精確度。因此把偏微分方程式離散化後，我們就稱此新的方程式為離散方程式。

考慮一個一維平板，長度為 L，在座標上其長度的方向為 x 方向，我們利用 離散的方法將空間0≤x ≤ L離散化[1]，

x

_i

=(i−1)∆x, i=1,2,...N

N 為空間中離散的總點數，而我們知道整個離散區域長度 L 與離散總點數 N，就 可求得點數間的距離，

−1

=

∆ N x L

在時間上，我們可以將時間 t 離散化，

t

_j

=(j−1)∆t, j=1,2,....J M 為時間間隔的總數，因此我們可以求得時間的間隔為，

1

max

= −

∆ J t t

有限差分法的好處在於方法不複雜，把方程式離散化的過程也很容易，可以很直 觀的利用一些數學計算來解決我們所遇到的問題。

2.2.3 2.2.3 2.2.3

2.2.3 一維 一維 一維熱傳導方程式 一維 熱傳導方程式 熱傳導方程式 熱傳導方程式

考慮一維與時間有關的熱傳導方程式[4]：

0 , 0

) , , ( )

, (

2 2

≥

≤

∂ ≤

= ∂

∂

∂ x L t

x t x T C

K t

t x T

ρ

(2.2) (2.2)式為一個典型的一維平板熱傳導方程式，其平板長度為 L，而(2.2)式 的解必須要在 x=0 與 x=L 給與邊界條件，其邊界條件假設為：

T(0,t)= ， T

₀ T ( L , t ) = T _L

， T(x,0)= f

₀

(x) (2.3)

2.2.4 2.2.4 2.2.4

2.2.4 差分方程式 差分方程式 差分方程式 差分方程式

我們要離散(2.2)式偏微分方程式，(2.2)式等號左邊的一階(first order) 微分式可由中央差分來近似，其型式為[1]：

t

t x T t t x T t

t x T

∆

−

∆

≈ +

∂

∂ ( , ) ( , ) ( , )

(2.4) (2.2)式等號右邊二階(second order)微分式可利用中央差分，其近似的型式為 [1]：

_{2 2}

2

) (

) , ( ) , ( 2 ) , (

) , (

x

t x x T t x T t x x T x

t x T

∆

∆

− +

−

∆

≈ +

∂

∂

(2.5)

把(2.4)式與(2.5)式取代(2.2)式，可得到一維熱傳導方程式離散化後的差分方 程式[1]，

)

2

(

) , ( ) , ( 2 ) , ( )

, ( ) , (

x

t x x T t x T t x x T C

K t

t x T t t x T

∆

∆

− +

−

∆

≈ +

∆

−

∆ +

ρ

(2.6)

由(2.6)式可以得知欲求這一刻的溫度，可以由前一刻各位置的溫度做簡單的運 算來求得。將方程式做簡單的處理，我們可得到其離散的型式為[1]：

( 2 ) )

(

^{2 1 1}

1 j

i j i j i j

i

T T T

x C

t

T

⁺

k

₊

− +

₋

∆

= ∆

ρ

(2.7) (2.7)式中，上標符號 j 表示時間的增量數，t=jΔt，我們可以看出在時間t

^{j + 1}

時 的溫度T

_i ^{j + 1}

是由前一個時間t 中的溫度組合。因此所以未知的溫度就可由前一時

^j

刻的溫度分布來求得。圖(2.1)描述(2.7)式的演算法，在點x=i∆x,t =(j+1)∆t的 溫度可由在較早時刻的三個點(x=(i−1)∆x,x=i∆x,x=(i+1)∆x,t= j∆t)的溫度求 得。圖中的 j 代表著在這一刻的時間t ，而 j+1 代表著下一刻我們要運算的時間

_j

+ 1

t

j

，而 j-1 代表著前一刻運算的時間t

_{j − 1}

。

圖中實心的方塊代表已知的初始值，空心的方塊代表邊界條件，空心的圓圈代表 由差分近似所計算的數值解。

2.3 2.3 2.3

2.3 二維熱傳導方程式與差分方程式 二維熱傳導方程式與差分方程式 二維熱傳導方程式與差分方程式 二維熱傳導方程式與差分方程式

在我們了解了一維的熱傳導方程式與差分方程式後，接下來我們將討論二 維熱傳導的情況，二維的情況較接近於真實的情況，例如一個二維平板，我們若 要了解在這個平板面上各個位置的溫度分佈，就可利用二維熱傳導方程式來計 算。

2.3 2.3 2.3

2.3.1 .1 .1 .1 二維 二維 二維熱傳導方程式 二維 熱傳導方程式 熱傳導方程式 熱傳導方程式

二維熱傳導方程式的型式為[5]：

( , , ) ( , , ) ) ) (

, , (

2 2 2

2

y t y x T x

t y x T C

K t

t y x T

∂ + ∂

∂

= ∂

∂

∂

ρ

(2.8) x，y 分別代表在直角座標系統上的兩個方向，T(x

_i

,y

_j

,t

_k

)代表在某一時間t 在空

_k

間座標(x

_i

,y

_j

)上的溫度，而(2.8)式中的係數與一維熱傳導方程式係數的定義相 同。

t=0

x=0

1 i-1 i i+1 t

N x=L j+1

j j-1

圖 2.1 有限差分法空間格子點的分佈與計算

2.3 2.3 2.3

2.3.2 .2 .2 .2 二維熱傳導差分方程式 二維熱傳導差分方程式 二維熱傳導差分方程式 二維熱傳導差分方程式

(2.8)式等號左邊偏微分式的差分型式為[1]：

t

t y x T t t y x T t

t y x T

∆

−

∆

≈ +

∂

∂ ( , , ) ( , , ) ( , , )

(2.9) (2.8)式等號右邊偏微分式的差分型式為[1]：

_{2 2}

2

) (

) , , ( ) , , ( 2 ) , , (

) , , (

x

t y x x T t y x T t y x x T x

t y x T

∆

∆

− +

−

∆

≈ +

∂

∂

(2.10)

_{2 2}

2

) (

) , ,

( ) , , ( 2 ) , ,

( ) , , (

y

t y y x T t y x T t y y x T y

t y x T

∆

∆

− +

−

∆

≈ +

∂

∂

(2.11)

將(2.9) ，(2.10) ，(2.11)式代入(2.8)式，且取∆x=∆y，我們可以得到二 微熱傳導差分方程式：

∆ +

∆

− +

∆

= +

∆

−

∆ +

)

2

(

) , , ( ) , , ( )

, , ( ) , , (

x

t y x x T t y x x T C

K t

t y x T t t y x T

ρ

₂

) (

) , , ( 4 ) , ,

( ) , ,

(

x

t y x T t y y x T t y y x T

∆

−

∆

− +

∆

+

(2.12)

將(2.12)式經過簡單的整理，我們可得到二微熱傳導的離散型式：

( 4 )

)

(

^{2 1 , 1 , , 1 , 1 ,}

, 1 ,

k j i k

j i k

j i k

j i k

j i k

j i k

j

i

T T T T T

x C

t T K

T + + + −

∆ + ∆

=

_{+ − + −}

+

ρ

(2.13)

(2.13)式為若要計算下一刻空間位置的溫度T(x

_i

,y

_j

,t

_{k + 1}

)，可以利用前一刻空間

位置的溫度T(x

_i

,y

_j

,t

_k

)以及此位置周圍四個位置的溫度

) , , ( ), , , ( ), , , ( ), , ,

(x

_{i 1}

y

_j

t

_k

T x

_{i 1}

y

_j

t

_k

T x

_i

y

_{j 1}

t

_k

T x

_i

y

_{j 1}

t

_k

T

_{+ − + −}

來求得。

2.3 2.3 2.3

2.3.3 .3 .3 .3 穩定性條件 穩定性條件 穩定性條件 穩定性條件

(2.13)式中比值

₂

) ( x C

t K

∆

∆

ρ

在二維的系統裡也有其限制條件，其原理與二維 相同，其比值必須小於或等於 1/4 的情況，

4 1 ) (

²

≤

∆

∆ x C

t K

ρ

(2.14) 才會使整個計算達到穩定而不發散的情況。

舉例說明，若假設初始空間間隔設定為 x∆ ，時間間隔為 t∆ ，當我們改變空

間間隔為2∆x的時候,時間間隔必須變為原先的四倍,也就是4∆t，這樣才會達到 收斂的條件，符合(2.14)式的結果。

2.4 2.4

2.4 2.4 二維平板熱傳導模擬 二維平板熱傳導模擬 二維平板熱傳導模擬 二維平板熱傳導模擬

考慮一個二維平板熱傳導的例子，此平板長與寬分別為 60cm 與 30cm，平板 的其中一邊保持 100℃，其他三邊則維持 0℃，求當保持穩定狀態時平板上的溫 度分佈。我們選擇這個例子來做為熱傳導模擬的分析主要是因為本篇論文是模擬 光碟片的熱傳導現象，我們的系統以二維的架構為主，而光碟是由許多的薄膜所 組成，每層薄膜都類似於一個二維平板的結構，因此我們選擇一個簡單的二維平 板結構模擬其熱傳導的現象，就可以進一步了解薄膜受熱的情況，因此二維結構 的光碟熱模擬就類似於一個二維的平板，所以我們可以由最簡單的二維平板例子 出發，進而推廣到光碟結構的情況。

2.4 2.4 2.4

2.4.1 .1 .1 .1 邊界條件 邊界條件 邊界條件 邊界條件

上述的邊界條件我們可以用以下的數學式來表示，而溫度為位置的函數，

因此其表示式為[5]:

C y

T

C y

T

C x

T

C x

T

°

=

°

=

°

=

°

=

100 ) , 60 (

0 ) , 0 (

0 ) 20 , (

0 ) 0 , (

(2.15)

圖(2.2)為二維平板結構的示意圖以及所對應的邊界條件。

2.4.2 2.4.2 2.4.2

2.4.2 熱傳導方程式與模擬結果 熱傳導方程式與模擬結果 熱傳導方程式與模擬結果 熱傳導方程式與模擬結果

由於此例子為不隨時變且無熱源無熱散的情況，因此熱傳導方程式就變 為 Laplace 方程式:

₂ 0

2 2

2 =

∂ + ∂

∂

∂

y T x

T

(2.16)

熱傳導方程式轉換成差分式後，其型式為[5]:

) 4(

1

0 4 ,

1 , 1 , , 1 , 1 ,

, 1 , 1 , , 1 , 1

− +

− +

− +

− +

+ + +

=

=

− + +

+

j i j i j i j i j

i

j i j i j i j i j i

T T T T T

T T

T T T

(2.17)

經由反覆迭代後，當溫度達到收斂時，就可以得到在這樣的邊界條件下空間的溫 度分佈，而反覆迭代的做法就是由(2.17)式裡計算出T

_{i , j}

，而T

_{i , j}

是由周圍四個點

的溫度平均來求得的，因此我們可以利用這樣的方法來求得整個空間平面上每個 點的溫度,然後在利用這些計算出來的溫度分佈做為下一次計算的起始值，於是 這樣反覆的迭代到所求得的溫度已經不再變化，最後的值就是我們所要求的解。

圖(2.3)為空間中溫度的分佈圖，從圖我們可以知道在靠近邊界條件 C

y

T(60, )= 100° 的地方,其溫度比靠近其他地區的位置要來的高，這符合了熱流 的物理觀念，也就是熱是由高溫流向低溫。

T=0℃

x y

T=0℃

60 10

T=0℃

T=100℃

圖 2.2 二維平板結構示意圖以及其邊界條件

0 20

40 60

0 10

20 30

0 20 40 60 80 100

y(cm) Temperature( ℃ )

x(cm)

圖 2.3 二維平板空間溫度分佈圖

2.4.3 2.4.3 2.4.3

2.4.3 誤差分析 誤差分析 誤差分析 誤差分析

在上一小節，我們計算了在二維熱傳導情況下的數值解，而這個例子是存 在著解析解，其解析解的型式為[5]:

∑ ^∞

=

+ +

= −

1

1

0 sin( )

) / sinh(

) / sinh(

1 ) 1 ( 2

n

n

L x n L w n

L y n T n

T π

π π

π

(2.18)

T 為在邊界

0

T(60,y)= 100°C的溫度，L 與 W 分別為二為平板的長與寬，因此我們 可以利用解析解與數值解所計算在相同位置上的溫度值來比較，可以知道數值計 算的誤差。而造成誤差的原因有兩種，一種是在利用有限差分法數值計算時，反 複迭代的次數不夠，也就是整個系統沒有達到收斂的情況，而另一種造成數值解

誤差很大的原因則為格點的間隔太大，因為在理論的情況下，若所取的格點大小

∆ 越小，其所得到的值會越精確，會越接近實際連續的情況。 x

在誤差分析上，我們做了 6 種情況，隨著 x∆ 的兩倍而減小也就是點數兩倍 的增加，然後每一種不同的情況下與解析解做誤差分析，於是把整個空間每個點 的誤差取平方後相加起來，此為

χ ²

(chi-squared)[1]。因此把每個情況的

χ ²

分 別與每個情況的 x∆ 做 log-log 圖，(2.16)式為與時間無關的二次偏微分熱傳導 方程式，因此當我們在討論格點大小的收斂以及誤差時，其誤差與∆ 有關，因 x

²

為我們利用中央差分來近似，而中央差分與格點大小平方有關，所以其近似誤差 為Ο(h

²

)項，所以對格點大小以及

χ ²

分別取 log，以格點大小為橫軸，

χ ²

為縱軸 做 log-log 圖，其斜率要為 2 才能表示誤差已達到收斂，圖 2.4 為誤差的 log-log 圖，我們共做了 6 種不同 x∆ 的情況，因此種共有 5 個斜率，其斜率分別為 2.4、

2.3、2.2、2.1、2.0，由斜率我們可以知道當我們所取的 x∆ 越小的時候，誤差 斜率會越接近 2，表示當 x∆ 越小值會越精確。另外要提到的就是當 x∆ 越小點數 越多的時候，所需要迭代的次數也必須要跟著變多，整體的數值才會達到收斂。

最後，

χ ²

的值越小表示越接近精確的值，因為我們是把數值解與解析解比 較，因此若解析解與數值解差值越小，表示越精確，而

χ ²

值越小表示整體的誤 差越精確。

2.5 2.5 2.5

2.5 穩定性與收斂性 穩定性與收斂性 穩定性與收斂性 穩定性與收斂性

在數值方法的運算中，我們必須要了解此方法的穩定性以及收斂性，這一章 節我們將針對有限差分法來分析其收斂性與穩定性。

2.5 2.5 2.5

2.5.1 .1 .1 .1 收斂與穩定的定義 收斂與穩定的定義 收斂與穩定的定義 收斂與穩定的定義

收斂(convergence)的意義為當 x∆ 與 t∆ 的值趨近於零的時候，由數值分析 所計算出的值必須與理論值所計算出的值一樣。穩定(stable)的意義是當運算步 驟上產生誤差時，此誤差不隨著運算步驟的增加而增大，相反的此誤差應隨著運 算而減小或是保持不變。

2.5 2.5 2.5

2.5.2 .2 .2 .2 熱傳導數值計算穩定條件 熱傳導數值計算穩定條件 熱傳導數值計算穩定條件 熱傳導數值計算穩定條件

在(2.7)式中比值

₂

) ( x C

t K

∆

∆

ρ

是一個影響數值計算穩定性的參數，在計算一維 熱傳導的情況時，比較必須要小於或等於 1/2 的情況

圖 2.4 斜率誤差圖,

χ ²

與格點大小的 log-log 圖

2 1 ) (

²

≤

∆

∆ x C

t K

ρ

(2.19) 才會使得整個數值計算達到穩定，當然隨著計算維度的不同,所需要達到穩定的 條件也會跟著改變。

我們以解析的方法來做穩定性之分析[5]，

X

AX =

λ

(2.20) 在上面的矩陣中，λ為矩陣 A 之固有值(eigenvalue)，而 X 為對應的固有向量 (eigenvector)。若矩陣 A 為-N×N 之矩陣，而其 N 個固有值均不相同的話，則其 N 個固有向量必定為線性獨立。考慮一個熱傳導的問題，其邊界條件固定，而其

內 部 有 N 個 未 知 的 溫 度 點 。 可 以 矩 陣 的 型 式 表 示 為 ，

 

 

 







 

 

 







=

 

 

 





−

 

 

 



 −

=

 

 

 







 

 

 







+ + +

p N

p p

p N p p

T T T F

F F

F

T T T

. .

. . .

) 2 1 ( ...

. . .

) 2 1 (

. .

.

2 1 0

0 0

0

1 1 2

1 1

(2.21)

以T

^{p +1}

= AT

^p

表示，A 表示上式之係數矩陣，T

^{p + 1}

與T 分別為在時間

^P

t

^{p + 1}

與t 時

^p

的溫度分佈。我們想要使得運算穩定，就必須控制誤差的變化，使誤差在運算後 小於運算前，也就必須使矩陣 A 之固有值之最大絕對值小於或等於 1。

矩陣 A 之固有值為，

) 1 ( sin 2 4 1

₀ ²

− +

N F n

π

n=1,2,...N (2.22) 其絕對值小於或等於 1:

) 1 1 ( sin 2 4 1

1

₀ ²

≤

− +

≤

− N

F n

π

(2.23) 因為F

₀

≥0，所以 2.16 式不等式變為，

) 1 ( sin 2 4 1

1

₀ ²

− +

≤

− N

F n

π

在改寫為，

) 1 ( sin 2

2 / 1

2 0

+

≤

N F n

π

(2.24)

由 2.24 式可知不等式右邊數值的最小值為(1/2) ，而

_{0 2}

) ( x C

t F K

∆

= ∆

ρ

所以一維穩定條件

2 1 ) (

²

≤

∆

∆ x C

t K

ρ

,必須要成立才會使得整體達到收斂。

Reference

1. R. H. Landau, M. J. P. Mejia,“Computational physics problem solving with computer”,(John Wiley and Sons,Inc.,New York ,1997) 2. H. S. Carslaw and J. C. Jaeger,“Conduction of heat in

solids”,(Oxford U. P.,London,1954)

3. 董建良,“應用有限元素法(Applied Finite Element Analysis)”,(科 技圖書股份有限公司,二版,1983)

4. M. D. Mikhailov, M. N. Ozisik,“Unified analysis and solutions of heat and mass diffusion”,(John Wiley and Sons,Inc.,1984) 5. 黃文雄,“熱傳學”,(中央圖書出版社,初版,1985)