渠道壁面流體滑移現象對微流道熱沈熱傳的影響

國立交通大學

機械工程研究所

碩士論文

渠道壁面流體滑移現象對微流道熱沈熱傳的影

響

The Effects of Fluid Slip at Walls on Heat

Transfer of Microchannels Heat Sink

指導教授：林振德教授

研 究 生：游本權

渠道壁面流體滑移現象對微流道熱沈熱傳的影

響

研究生：游本權 指導教授：林振德

摘要

隨著流體滑移現象漸漸在微米跟奈米尺度的微流道中被觀測出 來，實驗也證實當流體通過疏水性的微流道壁面會在壁面產生滑移的 現象。本文以計算流體力學的方法，針對流體通過壁面產生滑移現象 和未產生滑移現象進行模擬，也藉由紐賽數跟雷諾數的關係了解在不 同水力直徑下，不同形狀的微流道之熱傳行為及產生滑移後的影響。 分析結果顯示具滑移壁面的流道熱傳效果較佳。

The Effects of Fluid Slip at Walls on Heat

Transfer of Microchannels Heat Sink

Student: Ben-Chain Yu Advisor: Prof.Jenn-Der Lin

Abstract

Fluid slip has been observed experimentally in micro and nanoscale flow device by several investigators. Some

researcher examines a possible mechanism for the measured fluid slip, for water flowing over a hydrophobic surface. In this study, the computational fluid dynamics approach is employed to simulate and analyze the effect of fluid slip at walls on heat transfer. Results are presented in terms of Nusselt number at varies values of Reynolds number for microchannel heat sinks of different shape. Numerical results show that more effective of fluid slip at walls on heat transfer of

誌謝

這篇論文能夠完成，首先我要感謝林振德老師在這段期間的指 導，老師對於觀念上適切的引導，幫助我在進行論文研究時能夠有清 晰的思路，在此謝謝老師。 我非常珍惜在交大的每一個日子，感謝學長陳志堅、賴志銘、黃 永賢、張智通與潘柏霖對於我的照顧，還有昱宏、俞任、自皓在平時 的互相打氣，還有學弟懋勳、豐庭你們的幫忙，雖然相處的時間短暫， 但沒有你們就沒有這樣一段值得回味的實驗室時光，謝謝大家。 最後感謝我的家人及女友的支持與鼓勵，讓我能夠無後顧之憂地專心 於學業；最後，僅以本論文獻給所有關心我的人。

目錄

中文摘要...i 英文摘要...ii 誌謝 ...iii 表目錄 ...iv 圖目錄 ...v 符號表 ...x 第一章 緒論...1 1-1、前言...1

1-2、文獻回顧...3

1-3、研究目的...5 第二章 物理模式與數值方法... 7 2-1、物理模式...8 2-2、數值方法...12 2-2-1 離散方程式...13 2-2-2 數值模擬...14 2-2-2 SAMPLE Algorithm ...15 2-3 熱傳與 Nusselt Number...16 第三章 結果與討論...18

3-1 軟體數值碼的驗證與格點測試...18 3-2 網格測試...19 3-3 無滑移邊界條件與滑移邊界條件的熱傳效益...19 第四章 結論...23 參考文獻...24 表...30 圖...33

表目錄

表 4-1 圓形的幾何圖形與尺寸...30 表 4-2 矩形的幾何圖形與尺寸...31 表 4-3 矩形的幾何圖形與尺寸...32 表 4-4 熱沈與水物理性質...32 表 4-5 網格測試的格點數與總網格數...33

圖目錄

圖 1-1 IC 元件在封裝型態上的發展與演進 ...33 圖 1-2 IC 元件在引腳的發展與演進...33 圖 1-3 引起電子元件損壞的主要因素[1] ...34 圖 1-4 矽晶片表面上行成一氣體薄層...35 圖 1-5 流體經過表面有粒子的壁面情形...37 圖 2-1 物理模型...37 圖 2-2 矩形流道模擬剖面圖(流道長為 5mm) ...38 圖 2-3 梯形流道模擬剖面圖(流道長為 5mm)...38 圖 2-4 圓形流道模擬剖面圖(流道長為 5mm) ...39 圖 2-5 邊界條件示意圖...40 圖 2-6 氣體分子與管壁碰撞示意圖...41 圖 2-7 分子碰撞前後示意圖...42 圖 2-8 CFD-RC 求解過程...43 圖 2-9 二維三角格點...44 圖 3-1 文獻裡 y 方向(300μm)和 z 方向(30μm)的速度...45 圖 3-2 流道滑移模擬網格測試...46 圖 3-3 圓形流道入口流動情形...46 圖 3-4 矩形流道入口流動情形...47

圖 3-5 梯形流道入口流動情形...47 圖 3-6 不同雷諾數與水力直徑下圓形流道 Nu 的變化...48 圖 3-7 水力直徑為 50

μm

矩形流道在不同雷諾數下 Nu 的變化..48 圖 3-8 水力直徑為 100

μm

矩形流道在不同雷諾數下 Nu 的變化.49 圖 3-9 水力直徑為 200

μm

矩形流道在不同雷諾數下 Nu 的變化.49 圖 3-10 水力直徑為 100

μm

矩形流道在較低雷諾數下 Nu 的變化.50 圖 3-11 水力直徑為 50

μm

梯形流道在不同雷諾數下 Nu 的變化..50 圖 3-12 水力直徑為 100

μm

梯形流道在不同雷諾數下 Nu 的變化 51 圖 3-13 水力直徑為 200

μm

梯形流道在不同雷諾數下 Nu 的變化 51 圖 3-14 水力直徑為 100

μm

矩形流道相同截面機在不同雷諾數下 Nu 的變化...52 圖 3-151 水力直徑為 100

μm

矩梯形流道在不同滑移長度下(a)Zhu

and Granick (b)Tretheway and Meihart (c)Tyrell and Attard 不同雷諾數下 Nu 的變化...53

符號表

A 面積 C 未知數 f 流體 h 高度 q 分子速度 P 壓力 S 固體 T 溫度 U 速度 V 速度 u x 方向速度 v y 方向速度 w z 方向速度 ∗ u x 軸方向速度預測 u＇ x 軸方向速度修正值 ∗ v y 軸方向速度預測 v＇ y 軸方向速度修正值 ∗ w z 軸方向速度預測 w＇ z 軸方向速度修正值 β 滑動長度 μ 黏滯係數 α accommodation coefficient λ 平均自由路徑 Γ 界面

ρ 密度 Φ 任意相關變數 nb Φ 鄰近格點的任意相關變數

( )

∇Φ nb ∇Φ垂直到面的量 上標 * 預測值 , 修正值 下標 f 控制體積的面 g 氣體 w 水 t 紊流 Φ 與相關變數有關的值

第一章 緒論

1-1 前言

隨著半導體產業的高度發展與技術不斷的精進，現今電腦產業和 半導體製程的技術也有快速改良發展，特徵尺寸(Feature size)由 1980 年的 2 微米製程到 1995 年的 0.25 微米製程一直到現今的 0.18 微米製程，並且科學家預測在 2008 年時會發展為 0.07 微米的奈米製 程。由於特徵尺寸的極速下降，電子產品必朝向外型短小輕薄(如圖 1-1、1-2 所示)、性能提昇與價格降低的趨勢發展。由於裝置與元件 的發熱密度也因此不斷的增加，過熱的問題日益嚴重，於是造成了產 品的可靠度降低以及壽命的減少，所以散熱問題也就成為半導體產業 技術上的一大挑戰。根據美國空軍１９９５年所發表的研究報告中指 出，造成超過５０％的半導體元件損壞或缺陷的主要的原因來自於熱 的問題，如圖 1-3 所示。許多有關熱所造成的損壞與缺陷的研究文 獻，於各知名的期刊與國際研討會之文獻中歷歷可見[1-3]。

微機電系統(micro-electro-mechanical systems , MEMS)在歐洲 被稱為微系統科技(micro system technology)，著名的費曼博士 (1965 年諾貝爾物理獎得主)在 1959 年美國物理學年會上發表 「There＇s Plenty of Room at the Bottom」的專題演講中，首先 提到把機器微型化的概念，而「微機器(micromachines)」此一名詞

在 1978 年首次正式出現在國際學術研討會的名稱中；接著彼得森博 士 在 1982 年 發 表 了 著 名 的 「 以 矽 為 機 械 材 料 (silicon as a mechanical material)」研究報告，1989 年猶他州鹽湖城(Salt Lake City)的一場研討會(Micro-Tele-Operated Robotics Workshop)中， 則具體提出「微機電系統」此一名稱。而所謂微機電系統[4]是一種 結合機械、電子、材料、控制、物理、生醫、化學、光學等多重技術 整合的研發領域，其特色為兼具微小化、可量產之新製造技術，不僅 可以有效協助製造業改善製程，使產品提高品質、性能與可靠度，進 而提高附加價值，同時可降低製造成本及能源消耗。 就 微 流 體 系 統 (Micro-Fluidics system) 而 言 ， 以 微 流 道 (micro-channel)作為新型的熱沈(heat sink)與熱管之構型，時有所 聞。近幾年許多文獻均指出，層流流體在微流道熱沈中的熱傳效果優 於傳統大尺寸管路中的紊流態，而相變化過程所帶走的熱量更高，足 可應付未來高功率電子產品所產生的廢熱。而近年來，由於微機電 (MEMS)及奈米科技所興起的熱潮，故可製造出水力直徑(hydraulic diameter)在數十至數百微米的各種型態的流道，透過液體強制對流 的使用，將高功率電子產品產生的高熱降低，使之能於容許的溫度下 使用，來試圖解決電子產品過熱問題的損壞並增加其可靠度。

1-2 文獻回顧

如何將高功率電子元件所產生的廢熱排出外界是一項重要的課 題。近幾年由於微機電技術昌盛，製造數十至數百微米不同幾何形狀 的微流道熱沈已不是難事。最早提出結合微機電系統與微流道熱沈技 術概念的是 1981 年的Tuckerman and Pease[5]，其中探討微流道熱 沈的熱傳性能，並對微流道熱沈施以最佳化設計。實驗結果顯示最大

熱通量達到 790 W/cm2，可應用於高功率密度的超大型積體電路上

(VLSI)。Pfahler[6]、Choi[7]、Mohiuddin-Mala[8]等人亦提到微流 道熱沈不僅擁有超佳的散熱效果，其流體行為或熱傳特性與傳統大尺 寸管道的實驗結果大不相同。

Peng and Wang[9-13]則從事了一連串以相變化為主的微流體研 究。其發現當流道壁施以較高的熱通量或過熱溫度時，根據現有的曲 線預測應該產生核沸騰，結果卻沒有氣泡微流道中產生。換言之，流 道的幾何尺寸大小會嚴重影響流體的相變化行為，此對於兩相流之研 究而言，是非常有趣且重要的發現。近幾年對於微流道熱沈兩相流的 研究仍持續進行[14-16] 。

除了兩相流的研究外，Wu and Chang[17]針對微流道不同表面粗 糙度的情況下，探討流道的熱傳現象，其結果顯示出，粗糙度越大的 流道具有較高的 Nusselt number 與 Apparent friction，且此物理

性質會隨著雷諾數的增大而增大，但增大的比例有趨緩的現象。微流 道熱沈熱傳現象除實驗的驗證外，亦有一些科學家，運用計算流體力 學的方式與現有的實驗結果相互比較[18-19]，但在其壁面邊界的設 定上仍設為無滑動的邊界條件(no slip condition)，此邊界條件於 大尺寸的解析上並不會產生失真的現象，但對於尺寸下降至微米等級 時，是否合適尚待考量。 近一百多年來的科學家及工程師都以無滑移現象用在流體通過 固體表面的邊界條件上，但事實上，早在兩百多年前 Nevier 對於邊 界滑移已經提出一個假說，即滑移的速度與邊界表面的剪應力成正 比，亦即邊界滑移是成立的。邊界為何會滑移，就如同我們嘴巴吹氣 要將一個佈滿粒子的表面清乾淨，而較大的粒子比較容易被氣體帶走 而較小的粒子則緊緊的黏在表面上，如圖 1-1，這是因為越接近壁面 速度越接近零。近幾年，有人提出在微米或奈米等級的尺度下，在流 道壁面上會產成滑移的現象。Ruckenstein 和 Rajora[20]等人使用毛 細玻璃做為流道表面，讓液體流經疏水性的表面。而他們在實驗中發 現壓力降對邊界滑移所造成的影響遠比化學勢能理論(chemical potential theory)來的大，這個結果提供了邊界滑移可能是產生在 流道表面距離一個空隙的位置。Barrat 和 Bocquet[21] 則假設在有 奈米孔徑的流道下 non-wetting 所帶給流體滑移的環境。Pit[22]等

人根據了假設，在流道的表面塗佈了一層疏水性(hydrophobic)的

perfluorodecanetrichlorosilane 和 octadecyltrichlorosilane 而分別產生 170 跟 400nm 的滑移長度(slip length)。而 Derek and carl[23]更用高解析度的 Micro-PIV(Micron-resolution particle image velocimetry)拍攝進行邊界滑移的流體，觀測到微米等級下， 在管牆上會產生滑動的現象，其造成邊牆滑動的主因乃有一層氣體吸 附在邊牆上，因此產生液體與氣體界面的滑動。此氣體層的存在已由 James W. G. Tyrrell and Phil[24]在原子力顯微鏡(Atomic force microscopy)下觀測到(如圖 1-4)，Trethway 和 Meinhart[25]的實驗 結果也顯示在疏水性(hydrophobic)的微流道會產生滑移，相反的親 水性(hydrophilic)的微流道則不會產生滑移現象。在熱傳方面,Yu 和 Ameel[26]等人也對矩形流道中空氣流體的滑移現象所造的影響

1-3 本文研究目的

發熱量與散熱量能力為電子元件的可靠度及壽命的主要因素，根 據「10℃理論」[27]，當電子元件每升高 10℃，其壽命則相對減少 一半，可見溫度對 IC 的重要性，這也是熱管理(Thermal Management) 技術日益受到重視的原因。所謂的「熱管理」指的是電子元件中熱的

產生與熱控制的一門技術，它主要的目的是使電子元件晶片的接合處 溫度(Junction Temperature)維持在容許的範圍內，使得電子元件的 可靠度與壽命提升。因此電子元件熱的控制是相當重要的一大課題。 電子產品熱的生成主要是由於晶片中百萬個電晶體運作時所產生 的，這些問題雖然可由降低電壓的方式來減少發熱量，但是仍然不能 解決因功能提升及體積縮小所產生發熱密度增加的問題。 本文研究之目的，係以比較不同形狀與水力直徑微流道熱沈的熱 傳現象。此外，管牆的邊界條件以滑移與不滑移的條件加以探討熱傳 現象是否有其差異性。

第二章 物理模式與數值方法

2-1 物理模式

本研究主要採用 CFD-ACE(U)軟體進行模擬分析，微流道熱沈中 流道的設計上以矩形、梯形與圓形為其流道形狀，流道長 5mm，如圖 2-1 以矩型流道為例。研究中針對不同流道形狀、水力直徑與管道邊 牆是否滑動進行探討，滑移長度為 0.92μm[29]，探討其不同雷諾數 下 Nusselt Number 的變化情形，其模擬之三種不同流道結構剖面圖 如圖 2-2、2-3、2-4。 對於本文欲模擬結構的統御方程式與邊界條件如下所示: 連續方程式: 0 ) ( = ⋅ ∇ + ∂ ∂ v t r ρ ρ (2-1) 動量方程式: g p v v v t r r r r _{ρ τ ρ} ρ +∇⋅ =−∇⋅ +∇⋅ + ∂ ∂ ) ( ) ( ) ( (2-2) 能量方程式: ] ) [( ) ( ) ( h hv k k T t +∇⋅ =∇⋅ + t ∇ ∂ ∂ _{ρ ρ r (2-3)} 其中 g_{表物體重力(Body Force)} p_{表靜態壓力(Static Pressure)}

τ _{表應力張量(Stress Tensor)；} ] 3 2 ) [(∇vr +∇vrT − ∇⋅vrI = µ τ I表單位張量(Unit Tensor) h_表焓 本研究為穩定狀態，則統御方程式可改為: 連續方程式: 0 ) ( = ⋅ ∇ vr _(2-4) 動量方程式: g p v vrr τ ρr ρ =−∇⋅ +∇⋅ + ⋅ ∇ ( ) ( ) _(2-5) 能量方程式: ) ( ) ( hv =∇⋅ k∇T ⋅ ∇ ρ r (2-6) 對於固體區內，因為無流動現象，統御方程式為: ∇ ⋅(k∇T ) = 0 (2-7) 邊界條件: 1.在固體區: 有三個邊面設為絕熱的狀態，另外一邊面有一固定熱通量進入，如 圖 2-5 所示。在固體與液體的界面處其物理量連續，即溫度與能量 連續，其數學表示式如下所示: Γ Γ = , , f S T T (2-8)

Γ Γ _∂ ∂ − = ∂ ∂ − n T k n T k S _f f S (2-9) 其中 S 表固體 f 表流體 T 表溫度 Γ表界面 2.在液體區: 流體入口處:

0

0

293

=

=

=

=

w

v

u

u

K

T

in in 其中 u、v、w 表 x、y、z 方向速度，入口流數隨雷諾數而改變 流體出口處: 流體出口排放至標準狀況下，所以位於流道終端處會由軟體 程式自動計算。 在固體與液體界面處: 本論文考慮固液間無滑移與有滑移兩種情況，其邊界條件分訴如下 (Ⅰ)無滑動邊界: 固體與液體的界面處其物理量連續，即溫度與能量連續其數

學表示式如下所示: Γ Γ = , , f S T T (2-10) Γ Γ _∂ ∂ − = ∂ ∂ − n T k n T k S _f f S (2-11) (Ⅱ)滑動邊界: 滑移速度 Uslip 與滑移長度 β 之間的關係: 液體與物體間界面任何物理皆為連續，即速度連續、剪應力、 溫度與能量連續。 假設 Us 為滑移速度，Uw 為邊界流道壁面速度(固定壁面速度， Uw=0) w w slip dn dU U U 2 λ( ) α α − − = (2-12) 其中 α 表容納係數(accommodation coefficient) λ 表平均自由路徑(mean free path)

_w dn dU ) ( λ 表壁面剪應力 溫度滑移 w w s dn dT T T 22 λ( ) α α − − = (2-13) 在多年前，Navier 對於邊界滑移已提出一個假說，即滑移的速度與

剪應力成正比，如下式所示: ) ( ) 0 ( dy dV y V_x = = β x (2-14) 因為β的單位為長度，因此稱為滑動長度或滑動系數(slip length or

slip coefficient)，以下將根據 Navier 假說來探討 2-D 全展流下β

與壁面速度 Uslip 的關係： 圖 2-6 為流體層與管壁的放大圖，在圖中 與 分別表流體分子與 管壁碰撞前後 y 方向的切線速度，但對於流體分子與管壁的碰撞應有 兩種機制。第一種機制為彈性碰撞(Elastic collision)，即 等於 (分子與管壁碰撞的入射角與反射角相等)。第二種機制為擴散反 射(Diffuse reflection)，即分子會被吸附於管壁上，分子會跟隨著 管壁的速度運動，隨後可能再放射出去，因此 的速度會等於管壁 的速度。基於碰撞可能會有兩種機制，因此引進容納系數 (Accommodation coefficient) 1 q q₂ 1 q 2 q 2 q α，以便於將兩種不的機制加以整合，

而α的定義為分子經由擴散反射的比例(The fraction of molecules

undergoing diffuse reflection)。藉由容納系數的引進 可表為

下列的形式: 2

q

1 2 U (1 )q q = α wall + − α (2-15) wall U :管壁的速度 流體在管壁的上的邊界速度，取其碰撞前後的切線速度的平均值即

2 2 1 q q U_wall = + _(2-16) wall

U

:流體在管壁上的速度 將(2-15)式代入(2-16)式中可得: 1 2 2 2U q U_wall =α _wall + −α _(2-17) 上式仍無法解得流體在管壁上的邊界速度，為了求得其速度，可進 一步引進流體平均自由路徑的方法[18]，即分子碰撞前後的距離等於 λ 3

2 _{(λ:mean free path)，如圖 2-7 所示，因此 與牆的距離為}

1 q λ 3 2 _。 因此利用平均自由路徑的方法展開切線速度的函數 可得下式 (取到一階): ) ( y U_wall 0 1 ] ) ( [ 3 2 = + = y wall wall dy y dU U q λ (2-18) 將(2-18)式代入(2-17)式中可得下式 0 ] [ 3 2 2 = − + = y wall wall wall dy dU U U λ α α (2-19) 又管壁的速度為零，所以流體在管壁的速度為 0 ] [ 3 2 2 = − = y wall wall dy dU U λ α α (2-20) 對照(2-14)式我們可以得到滑移長度 λ α α β 3 2 2− = (2-21)

2-2 數值方法

CFD-ACE(U)解的步驟是，運用對控制體積積分的技巧，先將統

御 方 程 式 離 散 成 代 數 方 程 式 ， 也 就 是 一 般 所 稱 的 離 散 方 程 式 (Discretization Equations) ， 然 後 以 SIMPLEC (Semi-Implicit Method for Pressure-Linked Equation)演算法求得其正確壓力與速 度場，計算過程如 2-8 圖所示。較詳細的求解方式將在下面各節討論。 3-2-1 離散方程式 為了將統御方程式轉換為代數方程式，也就是為了穫得離散方程 式。離散方程式是由統御方程式對每個控制體積作積分而得。以一穏 態的統御方程式為例，對一任意的控制體積 V 積分的型式可表示如 下:

∫

ρφ

v

⋅

d

A

=

∫

Γ

φ

∇

φ

⋅

d

A

+

∫

V

S

φ

dV

r

r

r

(2-22) 其 中φ 表 任 意 相 關 變 數 ， vr _{表 速 度 向 量 ，} Ar _{表 面 積 向 量 以 及} 表擴散係數又， 表源項 φ Γ S_φ 式(2-22)可運用到計算域每一個控制體積。以一二維三角格點為 例，如圖 2-9 所示，式(2-23)可離散成下式: V S A A v _{n f} N f f f N f f faces faces φ φ

φ

φ

ρ

=

∑

Γ ∇ +

∑

r r ( ) r (2-23) 其 中 Nfaces 表 格 點 的 面 數 ， φf 表 透 過 面 f 轉 換 的 φ 值 ， 表每個面的質量通量， f f A vr ⋅ r ρ Af r 表面 的面積向量 ，f (∇φ)n為∇φ垂 直到面 的量 f

對於(2-23)式的對流項φf 值須進一步根據格點中心的值，運用 upwind scheme 作修改。 在式(2-23)的離散方程中，包含在格點中心未知變量φ與在鄰近 格點的未知值。一般而言，(2-23)式仍然可能還是非線性，因此進一 步線性化的行式可表示為下式:

∑

+ = nb nb nb p a b a

φ

φ

(2-24) 其中下標nb指的是鄰近的格點， 及ap anb是對於φ與φnb線性化的係數， 為源項。 b 3-2-2 數值模擬 數值模擬的計算乃藉由商用軟體CFDRC（Version 2003，CFD Research Corporation，Huntsville，AL）的CFD-ACE+TM模組，以有限 體積法進行三維的微流道混合數值模擬；此軟體專為模擬微流體流 動、熱傳等現象所設計，計算過程如 2-8 圖所示，操作步驟如下： 1. 軟體前處理器 CFD-GEOM：以內部之繪圖軟體繪製出微流體 裝置之形狀及構造，在此選用結構性網格(Structured grid)，由 於網格數(cell)大小直接影響到模擬的正確性及運算時間，為 了減少計算時間及電腦硬體本身對於運算網格數數目的限制 下，因此控制混合器網格數目數量約在 860000~930000 個。 接著定義邊界性質，如出入口、壁面等，完成後匯出檔案給

CFD-ACE-GUI 處理器讀取。 2. CFD-ACE-GUI處理器：給定流體性質如密度、黏滯係數及滑 移長度，選用水為模擬之流體。邊界條件我們依雷諾數給定 不同的流速，收斂標準為 10-5，進行穩態(Steady-State)之運算， 最後再將初始條件、疊代次數及其他參數設定完成後即開始 進行模擬。 CFD-VIEW 後處理器：將模擬結果溫度與速度場加以呈現。在求解溫 度與速度場時，則利用 SIMPLEC 法。 3-2-3 SIMPLEC Algorithm 壓力與速度修正主要是使在疊代的過程中，能滿足所有的物理 現象，計算出較符合實際狀況的值，其修正的過程如下: 1. 首先對壓力場 P*作預測。 2. 將 P*代入動量方程式求得相關的速度場 u*、v*與 w*。 3. 求得的 P*、u*、v*與 w*並不滿足連續方程式，因此運用連續方程 式建構壓力修正 P＇，相關的速度修正 u＇、v＇、w＇進而穫得， 因此可得到正確的壓力與速度場，如下所示: P =P*+ P＇ (2-25) u =u*+ u＇ (2-26)

v =v*+ v＇ (2-27) w =w*+ w＇ (2-28) 4. 假如計算出的值會因流體性質與源項影響流場，必須先解其它的 變數(如溫度、濃度、紊流量)，若流場不被影響，其餘變數值待 其收斂後再進行求解。 5. 將計算出的壓力 P 當為新的預測值，重複步驟 2 到 4 直到其解收 斂為止。

3-3 熱傳與紐賽數(Nusselt Number)

本文主要探討紐賽數 Nusselt Number(Nu)在不同流道形狀、流 體進口流速及固液邊界滑移情況，紐賽數的定義為對流熱傳係數與每 單位特徵長度流體熱傳導係數之比。為求得此值，主要需獲得對流熱 傳係數，此值可由下式計算而得: m W T A Q h ∆ = (2-29) ) ( out in p T T Mc Q = − (2-30) ) ( 2 1 out in m 流道壁面的平均度 T T T = − + ∆ (2-31) 其中 Q表傳入工作流體的熱量

W表整個流道的表面積

m

out in

A

∆T 表平均溫度差(mean temperature difference)

T 與T 表微流道進出口的溫度 將(3-28)與(3-30)代入 Nessult Number 定義式中，可獲得下式 m W f in out h p T A k T T D Mc Nu ∆ − = ( ) (3-32) 上式為本文計算 Nuessult Number 的主要公式。

第三章 結果與討論

本文採用水為工作流體，針對圓形、矩形與梯形與三種不同型態 的流道在不同水力直徑下比較其熱傳效益。其幾何圖形與尺寸如表 3-1 到 3-3 所示。相關熱沈與水物理性質如表 3-4 所示。在邊界條件 上，運用滑移與無滑移條件再作進一步的探討。

3-1 軟體數值碼的驗證與格點測試

為了判定 CFD-ACE(U)能運用於微米等級模擬與如何分佈格點的 疏密補捉正確的物理現象，本文引用文獻[29]的例子加以比較，其尺 寸為 30μm 高 300μm 寬 4cm 長的微流道，流速為 0.08m/s。在此的 比對，主要是與文獻所求得流體經過壁面產成滑動現象，求得流道 2cm 位置上 y 方向跟 z 方向的速度分佈，運用的格點數為 20x200x400 (z,y,x)。由圖 3-１的速度場分佈，我們可看出相同位置上，y 方向 上 y=o 的速度及文獻的誤差值為 0.221%，z 方向上 z=0 的速度及文獻 的誤差值為 0.25%，由此可證明， CFD-ACE(U)能運用於微米等級模 擬流體的滑動現象。另外，我們也在圖中觀察到，流體在 z 方向已經 達到完全發展流(fully develop)，而 y 方向尚未達到完全發展流。

3-2 網格測試

為了使內部網格區分大小、格點數，不至於影響到最後結果 的正確性。因此在數值模擬之前，格點的測試的工作有其必要性。通 常而言，網格區分越細密、則所計算出來的值也越精密，相對花費時 間也較多。但在某些特殊情況下，網格區分越細密也會導致結果不正 確及數值不穩定現象。因此可藉由格點測試的工作，有限的電腦資源 中尋求網格分部的最佳化，本數值模擬以方形流道為例，邊長為 50μm，在雷諾數等於 50 下，探討網格疏密與結果的影響。在滑移邊 界條件下測試 A、B、C、D 四種網格(表 3-5)測試結果如圖 3-2 所示， 因 85000 與 95000 網格數的出口溫度之間的誤差為 0.001%，故我們 以此為設定網格的基準，而流道格點最密的區域主要位於流道內與流 體進出口區、固體區靠近流道處與熱源處，如此的加密方式，主要是 此三處為物理量梯度變化最顯著的區域。

3-3 無滑移邊界與滑移邊界的熱傳效益

微流道熱沈擁有超佳的散熱效果，已經由許多文獻所證實。在 本節中，首先探討無滑移邊界下，不同流道形態與水力直徑的差異 性，和滑移邊界下與無滑移邊界流道之間的差異性。在以往無滑移的

邊界條件運用於大尺寸的流道中，其所獲得的結果並不會產生失真的 現象，但當流道的特徵尺寸(如水力直徑)小到微米或奈米等級的情況 下，再將邊界條件視為無滑移的現象，是否能與真實情況吻合需再進 一步加以確認。 圖3-3 到 3-5 為圓形、矩形和梯形在流道入口處的流速分 佈，由圖形可以看出，流道的入口 x=0 的位置，流速分佈相當平均， 差別在於邊界滑移的流道壁面有速度，隨著邊界層(boundary layer) 逐漸變小最後到達完全發展流。我們可以發現從入口到完全發展流的 這段入口區域 (entrance region)皆以邊界滑移時較小，也就是在 滑移邊界下，流體較快達到完全發展流，有邊界滑移的流道跟無邊界 滑移的流道的入口區域 距離差距約為20μm。 L X L X 圖 3-6 為圓形流道在不同水力直徑與雷諾數(Re)下，紐賽數 Nusslet Number (Nu)的變化情形。由此圖可知，在圓形流道中，紐 賽數會隨雷諾數與水力直徑的變大而變大，考慮邊界滑移的情況下， 其紐賽數較無滑移情況下大，且在滑移邊界下，紐賽數隨著雷諾數的 變化比較大。此外，當水力直徑越大增加的量也會越大。 圖 3-7 到 3-9 為矩形流道在水力直徑為 50、100、200μm 在不同 雷諾數(Re)下，紐賽數(Nu)的變化情形，圖中顯現三種不同高寬比的 流道，其幾何尺寸列於表 3-2。由圖觀察到，當在不同高寬比及相同

的水力直徑下會有不同的紐賽數。在相同水力直徑下，高寬比越大時 矩形流道擁有較大的截面積，所以造成不同的熱傳效果，也就是截面 積越大質量流率越大，紐賽數也越大。而由三個圖比較水力直徑由 50、100、200 其變化也有逐漸邊大的趨勢。而由圖 3-10 我們也可以 看出，當雷諾數低於 100 時，紐賽數的變化量加劇，而滑移邊界與無 滑移邊界的紐賽數差距也較小。此外紐賽與雷諾數的變化與圓形流道 相同；但在滑移邊界下，矩形流道的變化量較趨於平緩，其變化趨勢 大致跟無滑移邊界條件相同。 圖 3-11 到 3-13 梯形流道在水力直徑為 50、100、200μm 在不同 雷諾數(Re)下，紐賽數(Nu)的變化情形。將梯形流道與矩形流道相互 比較，發現在不同邊長的比例相同水力直徑下，同一雷諾數 Nu 大約 相同。造成此種現象的主因我們認為，雖然邊長的比例不同水力直徑 相同，但梯形在不同邊長的比例下截面積較為接近，通過流道的流量 大約相同也使熱傳的效益差距不大，才會造成較為類似的結果。 圖 3-14 為矩形在水力直徑為 100μm 下，高寬比為倒數，比較其 不同雷諾數(Re)下與紐賽數(Nu)的變化情形。在此我們可以發現，相 同水力直徑及截面積的流道，因固體區的溫度分佈相當平均，也就是 流道各個壁面的熱通量接近，使其紐賽數結果相同。圖 3-15 我們以 水力直徑為 100μm ，typeⅠ的矩形流道為例，比較在文獻中(a) Zhu

和 Granick 所計算出的滑移長度 β=2.5、(b) Trethway 和 Meihart 藉由 μ-PIV 量測出滑移長度 β=0.92 及(c) Tyrell 和 Attard 量測 奈米氣泡(nanobubbles)測出的滑移長度 β=0.5，三種不同滑移長度 紐賽數隨雷諾數的變化情形，其結果顯示出滑移長度越大，散熱效果 越佳。

第四章 結論

本研究主要探討微流道邊界滑移對微流道熱沈的熱傳現象。在微 流道熱傳方面，主要探討圓形、矩形與梯形三種微流道紐賽數隨雷諾 數隨雷諾數的變化。當中我們發現，矩形流道在相同的水力直徑下，， 若矩形流道的高寬比越大，其熱傳效果會越好，但對於梯形流道與圓 形流道而言，不同的邊長的比例相同的水力直徑下，截面積影響力最 大。 在邊界滑移部份，具邊界滑移的流道形成完全發展流較快。在熱 傳效益上，結果顯示出具滑移邊界的流道熱傳效益較佳，且當滑移長 度越大時，流道有更好的熱傳效果。 因此，經過本文之探討，能讓我們對於渠道壁面流體滑移現象對 微流道熱沈熱傳有更深一層的視野，並期望能在壁面滑移的研究上能 帶給人們更廣的視野及提供日後研究流道壁面滑移者有效之參考資 料。

參考文獻

[1] L. T. Yeh, “Review of Heat Transfer Technologies in Electronic Equipment,＂ J. of Electronic Packaging, Vol. 117, pp. 333-339, 1995.

[2] C. Chapman, “The Basics of Package/Device Cooling,＂ Electronic Packaging and Production, pp. 57-60, May 1998. [3] Tummala, “Fundamentals of Microsystems Packaging,＂ Mc

Graw-Hill, 2002.

[4] T. R. Hsu,MEMS and Microsystems：design and manufacture McGraw-Hill,2002

[5] D. B. Tuckerman, R. F. W. Pease, ＂High-Performance Heat Sinking for VLSI＂, IEEE Electronic Device Letter , Vol. EDL-2(5), No. 4, pp. 126-129, 1981.

[6] J. Pfahler, J. Harley, H. H. Bau, and J. Zemel, “Liquid Transport in Micron and Submicron Channels＂, sensors and Actuators, A21-23, pp. 181-314, 1990.

[7] S. B. choi, R. F. Barron, R. O/ Warrington, “Liquid Flow and Heat Transfer in Microtubes＂, ASME Micromechanical Sensors, Actuators and Systems, DSC-Vol. 32, pp. 123-128,

1991.

[8] Weilin Qu, Mohiuddin Mala, Dongqing Li, “Presure-Driven Water Flows in Trapezoidal Silicon Microchannels＂, Int. J. Heat and Mass Transfer, Vol. 43, pp. 553-364, 1999. [9] X. F. Pen, B. X. Wang, “Forced Convection and Flow Boiling

Heat Transfer for Liquid Flowing through Microchannels＂, Int. J. Heat Mass Transfer, Vol. 36, No.14, pp. 3421-3427, 1993.

[10] B. X. Wang, X. F. Peng, “Experimental Investigation on Liquid Forced Convection Heat Transfer through

Microchannels＂, Int. J. Heat and Mass Transfer, Vol. 37 Suppl. 1, pp. 73-82, 1994.

[11] X. F. Peng, G. P. Peterson, B. X. Wang, “Frictional Flow Characteristics of Water Flowing through Rectangular Microchannels＂, Experimental Heat Transfer, Vol. 7, pp. 249-264, 1994

[12] X. F. Peng, G. P. Peterson, B. X. Wang, “Heat Transfer Characteristics of Water Flowing through Microchannels＂, Experimental Heat Transfer, Vol. 7, pp. 2659-283, 1994.

[13] X. F. Peng, G. P. Peterson, B. X. Wang and H. B. Ma, “Experimental Investigation of Heat Transfer in Flat Plates With Rectangular Microchannels＂, Int. J. Heat and Mass Transfer, Vol. 38, No. 1, pp. 127-137, 1995. [14] D. Brutin, F. Topin, L. Tadrist, “Experimental Study of

Unsteady Convective Boiling in Heated minichannels＂, Int J. of Heat and Mass Transfer, Vol. 46, pp. 2957-2965, 2003.

[15] Linan Jiang, Man Wong, “Phase Chang in Microchannel Heat Sinks with Integrated Temperatuere Sensors＂, J. of

Microelectromechanical Systems, Vol. 8, No. 4, pp. 358-365, 1999.

[16] Lian Zhang, Jae-Mo Koo, Linan Jiang, “Measurements and Modeling of Tow-Phase Flow in Microchannels with Nearly Constant Heat Flux Boundary Conditions＂. J. of Microelectromechanical Systems, Vol. 11, No. 1, pp. 12-19, 2002..

[17] H. Y. Wu, Ping Chang, “An Experimental Study of Convection Heat Transfer in Silicon Microchannels with Different Surface Condition＂, Int. J. of Heat and Mass Transfer , Vol.

46, pp. 2547-2556, 2003

[18] K. C. Toh, X. Y. Chen, J. C. Chai, “Numerical Coputation of Fluid Flow and Heat Transfer in Microchannels＂, Int. J. of Heat and Mass Transfer, Vol. 45, pp. 5133-5141, 2002 [19]Andrie G. Fedorov, Raymond Viskanta, “Three-Dimensional

Conjugate Heat Transfer in the Microchannel Heat Sink for Electronic Packaging＂, Int. J. of Heat and Mass Transfer, Vol. 43, pp. 399-415, 2000

[20]Ruckenstein,E. and Rajora, P.“On the No-slip Boundary Condition of Hydrodynamics＂.J. of Colloid and Interface Science,96,pp. 488-491,1983

[21]Barrat,J.and Bocqet,L.“Large Slip Effect at a Nonwetting Fluid-Solid Interface＂.Physical Review Letters,82,pp. 4671-4674,1999

[22]Pit, R., Hervet, H.,and Leger, L.“Direct Experimental Evidence of Slip in Hexadecane: Solid Interface＂.Physical Review Letters,85,980-983,2000

[23]Derek C. Tretheway and Carl D. Meinhart “Apparent fluid slip at hydrophobic microchannels walls,＂ Physics of Fluids,

Vol. 14, N0.3, 2002.

[24]James W. G. Tyrrell and Phil Attard “Images of Nanobubbles on Hydrophobic Surfaces and Their Interaction,＂Physical Review Letters, Vol. 87, NO. 17, 2002.

[25]Derek C. Tretheway, Xiaojun Liu, and Carl D. Meinhart “Analysis of Slip Flow in Microchannel＂Department of Mechanical Engineering University of California, Santa Barbara. CA 93106, 2003.

[26]Ship Yu, Timothy A.Ameel,“Slip-flow heat transfer in rectangular microchannels＂ Int. J. of Heat and Mass Transfer,Vol.44,pp.4225-4234,2001

[27] 鐘文仁，“IC 封裝製程與 CAE 應用＂全華出板社，2003。 [28] Weilin Qu, Issam Mudawar, “Analysis of Three-dimensional

Heat Transfer in Micro-channel Heat Sinks＂, Int. J. of Heat and Mass Transfer, Vol. 45, pp. 3973-3985, 2002. [29] Derek C. Tretheway, Xiaojun Liu, and Carl D. Meinhart,

“Examination of the Slip Boundary Condition By μ-PIV and Lattice Boltzmann Simulations＂,ASME,2002

questions about complex fluids flowing past solids＂ nature materials, VOL 2, 2003

表 4-1 圓形的幾何圖形與尺寸

表 4-3 梯形的幾何圖形與尺寸 表 4-4 熱沈與水物理性質 Tin (K) Heat Flux (W/m2) Kwater (W/m-K) Ksilicon (W/m-K) Cpwater (J/kg-K) υ (m2/s) 293 360000 0.61 148 4179 1*E-6

表 4-5 網格測試的格點數與總網格數 (y,z)格點數 總網格數 網格測試 A 75x36 662904 網格測試 B 79X39 763344 網格測試 C 82x42 859599 網格測試 D 87x44 960876

圖 1-1 IC 元件在封裝型態上的發展與演進

Temperature 55% Humidity 19% Dust 6% Vibration 20% 圖 1-3 引起電子元件損壞的主要因素[1]

圖 2-2 矩形流道模擬剖面圖(流道長為 5mm)

圖 3-1 文獻與 CFD-RC 所作出 y 方向(300μm)和 z 方向(30μm)的速 度

圖 3-2 流道滑移模擬網格測試

圖 3-4 矩形流道入口流動情形

圖 3-6 不同雷諾數與水力直徑下圓形流道 Nu 的變化

圖 3-8 水力直徑為 100

μm

矩形流道在不同雷諾數下 Nu 的變化

圖 3-10 水力直徑為 100

μm

矩形流道在較低雷諾數下 Nu 的變化

圖 3-12 水力直徑為 100

μm

梯形流道在不同雷諾數下 Nu 的變化

圖 3-14 水力直徑為 100

μm

矩形流道相同截面積在不同雷諾數下 Nu 的變化

圖 3-15 水力直徑為 100

μm

矩梯形流道在不同滑移長度下(a)Zhu and Granick (b)Tretheway and Meihart (c)Tyrell and Attard