低壓縮比蠕動式微幫浦分析

國 立 交 通 大 學

機 械 工 程 學 系

碩士論文

低壓縮比蠕動式微幫浦分析

Analysis of low-compression ratio peristaltic micropump

研 究 生

： 張 祚 昌

指 導 教 授

： 崔 燕 勇

博士

低壓縮比蠕動式微幫浦分析

Analysis of low-compression ratio peristaltic micropump

研 究 生：張祚昌 Student：Tso-Chang Chang

指導教授：崔燕勇 Advisor：Yeng-Yung Tsui

國 立 交 通 大 學

機 械 工 程 學 系

碩士論文

A Thesis

Submitted to Department of Mechanical Engineering College of Engineering

National Chiao Tung University in partial Fulfillment of the Requirements

for the Degree of Master of Science

In

Mechanical Engineering July 2009

Hsinchu, Taiwan, Republic of China

低壓縮比蠕動式微幫浦分析

研究生：張祚昌 指導教授：崔燕勇 博士

國立交通大學機械工程學系

摘 要

本研究中探討二腔及三腔式低壓縮比的蠕動式微幫浦。壓縮比約為

0.00706。測試不同的蠕動相位及背壓下，對微幫浦淨流量的影響。此外我

們利用一套塊狀理論，分析微幫浦的流場，在非穩態的模式下可預測微幫

浦內壓力及流量的變化。由模擬的結果可知，逆向三相位的蠕動形式，在

低壓縮比的蠕動式微幫浦中能產生最大的淨流量，另外微幫浦的腔體數目

增加後，可承受較大的背壓。

Analysis of low-compression ratio peristaltic micropump

Student：Tso-Chang Chang Advisor：Dr. Yeng-Yung Tsui

Department of Mechanical Engineering National Chiao Tung University

ABSTRACT

A double and triple chambers low-compression ratio peristaltic micropump is presented in this thesis. The compression ratio is almost 0.00706. We calculate the net flow rate of micropump with different peristaltic motion and different back pressures. Furthermore, the flow of micropump is also analyzed by using lump system. The unsteady model is capable of predicting the pressure and flow rate of micropump inside. By the simulation result, the peristaltic motion of reverse three-phase achieved a maximum flow rate. Besides, the micropump can sustain maximum backpressure when increase the chambers of micropump.

誌謝

在這碩士的兩年生活之中，非常感謝崔燕勇教授在研究上給我的教導與督促，當我 遇到困難時，老師都能適時的給我一個指點，使我能順利完成研究所的學業，真的是由 衷的感謝。在研究中遇到困難時，也很感謝吳添成、胡育昌、吳欣恩、林仕文、呂學霖 等學長的幫助及指導，使我更能融入研究之中並能順利的解決問題。在研究室的生活之 中，則非常感謝同學俊佑及孝修在課業及生活上的鼓勵及幫忙，還有感謝學弟大慶、光 桓及信宏在生活上的協助，讓我有個難忘的兩年新竹生活。也感謝隔壁實驗室同學仁 鈞、俊傑等，及在新竹的大學朋友在平時能互相鼓勵，使我能在繁忙的研究生活之中， 能獲得片刻的舒緩。 最後特別感謝我的家人，尤其是我的爸媽全力支持我讀研究所，讓我無後顧之憂的 專心在課業上，還有我外婆、妹妹及家人平日的關心問候，在此說聲深深的感謝，謝謝 你們的支持，謝謝。

目錄

中文摘要...i 英文摘要...ii 誌謝...iii 目錄...iv 表目錄...vi 圖目錄...vii 符號說明...x 第一章 簡介 ...1 1.1 前言...1 1.2 微幫浦簡介...2 1.2.1 無閥式微幫浦簡介 ...3 1.2.2 蠕動式微幫浦簡介 ...3 1.3 文獻回顧...5 1.4 研究目的...9 第二章 數學模式 ...10 2.1 基本假設...10 2.2 統御方程式...10 2.4 邊界條件...10 2.4.1 對稱邊界條件 ...10 2.4.2 壓電薄膜邊界條件 ... 11 2.4.3 出入口端邊界條件 ...12 2.4.4 壁面邊界條件 ...12 第三章 方程式離散化 ...13 3.1 有限體積法...13 3.2 離散動量方程式...13 3.2.1 非穩態項 ...14 3.2.2 對流項 ...14 3.2.3 擴散項 ...15 3.2.4 源項...16 3.3 代數方程式...16

第四章 PISO 演算法 ...18 4.1 壓力與速度修正...18 4.1.1 預測步驟 ...18 4.1.2 第一次修正 ...18 4.1.3 第二次修正 ...20 4.1.4 修正後方程式整理 ...22 4.2 邊界條件設定...22 4.2.1 進出口邊界條件之流量計算 ...22 4.2.2 進出口邊界條件之速度計算 ...23 4.2.3 固定牆邊界條件 ...23 4.3 PISO 演算步驟...23 第五章 微幫浦塊狀理論分析 ...25 5.1 塊狀系統理論分析(LUMP-SYSTEM ANALYSIS)...25 5.1.1 單腔無閥式 ...26 5.1.2 二腔蠕動式 ...27 5.1.3 三腔蠕動式 ...30 第六章 結果與討論 ...33 6.1 網格點測試 ...33 6.2 各式蠕動相位的特性分析...33 6.2.1 壓縮比(ε ) ...34 6.2.2 蠕動相位介紹 ...34 6.2.3 順向三相位與逆向三相位 ...34 6.2.4 對稱性的蠕動相位 ...35 6.3 塊狀理論分析...36 6.3.1 流向分析 ...36 6.3.2 淨流量分析 ...38 6.4 CFD 模擬結果...38 6.5 腔體數的影響...41 第七章 結論 ...42 參考文獻...43

表目錄

表6.1 三腔(蠕動形式：順向三相位) ...46 表6.2 三腔(蠕動形式：逆向三相位) ...46 表6.3 三腔(蠕動形式：四相位) ...46 表6.4 三腔(蠕動形式：六相位) ...46 表6.5 二腔(蠕動形式：順向三相位) ...47 表6.6 二腔(蠕動形式：逆向三相位) ...47 表6.7 二腔(蠕動形式：四相位) ...47

圖目錄

圖1.1 無閥式微幫浦一個週期的運作【28】 ...48 圖1.2 二腔蠕動式微幫浦一個週期的運作：(A)三相位；(B)四相位【17】 ...48 圖1.3 三腔蠕動式微幫浦一個週期的運作：(A)三相位；(B)四相位；(C)六相位【21】 ...49 圖1.4 壓縮比示意圖【29】 ...49 圖1.5 高壓縮比蠕動相位示意圖 (A)二腔三相位；(B)二腔四相位 ...50 圖1.6 高壓縮比蠕動相位示意圖 (A)三腔三相位；(B)三腔四相位；(C)三腔六相位....51 圖1.7 低壓縮比蠕動相位示意圖 (A)二腔順向三相位；(B)二腔逆向三相位；(C)二腔四 相位...52 圖1.8 低壓縮比蠕動相位示意圖 (A)三腔順向三相位；(B)三腔逆向三相位(C)三腔四相 位；(D)三腔六相位...53 圖1.9 三腔蠕動式微幫浦尺寸示意圖 ...54 圖2.1 模擬壓電薄膜振動的形式-梯形曲線...55 圖3.1 非結構性網格示意圖 ...55 圖3.2 OVER-RELAXED APPROACH ...55 圖4.1 壓力邊界設定示意圖 ...56 圖6.1 相異網格之中，三腔逆向三相位每一週期的淨流量值(漸縮/漸擴管) ...58 圖6.2 塊狀理論分析中，三腔逆向三相位每一週期的淨流量值(漸縮/漸擴管) ...58 圖6.3 三腔逆向三相位在不同背壓下的淨流量(漸縮/漸擴管) ...59 圖6.4 二腔順向三相位(背壓=0)流量隨時間變化圖(漸縮/漸擴管) ...60 圖6.5 二腔順向三相位(背壓=0)壓力隨時間變化圖(漸縮/漸擴管) ...61 圖6.6 二腔逆向三相位(背壓=0)流量隨時間變化圖(漸縮/漸擴管) ...62 圖6.7 二腔逆向三相位(背壓=0)壓力隨時間變化圖(漸縮/漸擴管) ...63 圖6.8 三腔順向三相位(背壓=0)流量隨時間變化圖(漸縮/漸擴管) ...64 圖6.9 三腔順向三相位(背壓=0)壓力隨時間變化圖(漸縮/漸擴管) ...65 圖6.10 三腔逆向三相位(背壓=0)流量隨時間變化圖(漸縮/漸擴管) ...66 圖6.11 三腔逆向三相位(背壓=0)壓力隨時間變化圖(漸縮/漸擴管)...67

圖6.12 二腔四相位(背壓=0)流量隨時間變化圖(漸縮/漸擴管) ...68 圖6.13 二腔四相位(背壓=0)壓力隨時間變化圖(漸縮/漸擴管) ...69 圖6.14 三腔四相位(背壓=0)流量隨時間變化圖(漸縮/漸擴管) ...70 圖6.15 三腔四相位(背壓=0)壓力隨時間變化圖(漸縮/漸擴管) ...71 圖6.16 三腔六相位(背壓=0)流量隨時間變化圖(漸縮/漸擴管) ...72 圖6.17 三腔六相位(背壓=0)壓力隨時間變化圖(漸縮/漸擴管) ...73 圖6.18 二腔四相位(背壓=0)流量隨時間變化圖(直管)...74 圖6.19 二腔四相位(背壓=0)壓力隨時間變化圖(直管)...75 圖6.20 三腔四相位(背壓=0)流量隨時間變化圖(直管)...76 圖6.21 三腔四相位(背壓=0)壓力隨時間變化圖(直管)...77 圖6.22 三腔六相位(背壓=0)流量隨時間變化圖(直管)...78 圖6.23 三腔六相位(背壓=0)壓力隨時間變化圖(直管)...79 圖6.24 不同背壓下，三腔逆向三相位的淨流量比較(漸縮/漸擴管) ...80 圖6.25 不同薄膜最大振幅之下，藉由塊狀分析及 CFD 計算的淨流量值(漸縮/漸擴管) ...80 圖6.26 不同薄膜最大振幅及背壓下，藉由塊狀分析及 CFD 計算的淨流量值(漸縮/漸擴 管)...81 圖6.27(A) 二腔順向三相位 1/6 週期時，腔內壓力與流線圖 ...82 圖6.27(B) 二腔順向三相位 3/6 週期時，腔內壓力與流線圖...83 圖6.27(C) 二腔順向三相位 5/6 週期時，腔內壓力與流線圖...84 圖6.28(A) 二腔逆向三相位 1/6 週期時，腔內壓力與流線圖 ...85 圖6.28(B) 二腔逆向三相位 3/6 週期時，腔內壓力與流線圖...86 圖6.28(C) 二腔逆向三相位 5/6 週期時，腔內壓力與流線圖...87 圖6.29(A) 二腔四相位 1/8 週期時，腔內壓力與流線圖 ...88 圖6.29(B) 二腔四相位 3/8 週期時，腔內壓力與流線圖...89 圖6.29(C) 二腔四相位 5/8 週期時，腔內壓力與流線圖...90 圖6.29(D) 二腔四相位 7/8 週期時，腔內壓力與流線圖 ...91 圖6.30(A) 三腔順向三相位 1/6 週期時，腔內壓力與流線圖 ...92 圖6.30(B) 三腔順向三相位 3/6 週期時，腔內壓力與流線圖...93

圖6.30(C) 三腔順向三相位 5/6 週期時，腔內壓力與流線圖...94 圖6.31(A) 三腔逆向三相位 1/6 週期時，腔內壓力與流線圖 ...95 圖6.31(B) 三腔逆向三相位 3/6 週期時，腔內壓力與流線圖...96 圖6.31(C) 三腔逆向三相位 5/6 週期時，腔內壓力與流線圖...97 圖6.32(A) 三腔四相位 1/8 週期時，腔內壓力與流線圖 ...98 圖6.32(B) 三腔四相位 3/8 週期時，腔內壓力與流線圖...99 圖6.32(C) 三腔四相位 5/8 週期時，腔內壓力與流線圖...100 圖6.32(D) 三腔四相位 7/8 週期時，腔內壓力與流線圖 ...101 圖6.33(A) 三腔六相位 1/12 週期時，腔內壓力與流線圖 ...102 圖6.33(B) 三腔六相位 3/12 週期時，腔內壓力與流線圖...103 圖6.33(C) 三腔六相位 5/12 週期時，腔內壓力與流線圖...104 圖6.33(D) 三腔六相位 7/12 週期時，腔內壓力與流線圖 ...105 圖6.33(E) 三腔六相位 9/12 週期時，腔內壓力與流線圖...106 圖6.33(F) 三腔六相位 11/12 週期時，腔內壓力與流線圖 ...107 圖6.34 單腔及雙腔微幫浦在不同背壓下淨流量比較(漸縮/漸擴管)【26】 ...108 圖6.35 二腔及三腔蠕動式微幫浦在逆向三相位時，淨流量隨背壓的變化(漸縮/漸擴管) ...108 圖6.36 二腔及三腔蠕動式微幫浦在四相位時，淨流量隨背壓的變化(漸縮/漸擴管) ..109

符號說明

A_l 漸縮/漸擴管右側擴散端面積 A t 漸縮/漸擴管左側喉部面積 a 蠕動相位數 d_MAX 壓電薄膜最大振動幅度 f 壓電薄膜振動頻率 K 漸縮/漸擴管造成的壓力損失係數 l 漸縮/漸擴管的長度 m 質量流率 n 壓電薄膜作一簡諧運動所劃分的 timestep 數 P 壓力 Q 體積流率 QΔ 腔體因壓電薄膜位移造成的體積變化率 r ₀ 腔體半徑 r₁ 壓電薄膜半徑 T 無閥式微幫浦運作一個週期所需的時間 VK 流體速度 α_Q 塊狀理論非穩態的權重係數 ε 壓縮比 μ _{流體的黏滯係數} ρ _{流體的密度} τ 蠕動式微幫浦運作一個週期所需的時間 ∀MAX 薄膜掃過的最大體積

第一章 簡介

1.1 前言

近年來半導體逐漸成為台灣最重要的產業，使得相關的微機電系統(Micro Electro Mechanical System，MEMS)製程技術也隨之蓬勃發展。機械微型化的概念最早由費曼博 士(R. P. Feynman)在 1959 年所提出，接著 1982 年彼得森博士(K. E. Petersen)發表了一篇 「Silicon as a mechanical material，以矽作為機械材料」研究報告，至 1989 年確立了「微

機電系統」此一名稱【1】。微機電系統技術是一種結合光學、機械、電子、材料、控制、

物理、化學、生醫等多重領域的整合型微小化系統製造技術，其具有較小的體積、重量 及耗能等優點，同時亦能降低其製造成本。利用微機電系統技術可製造出微制動器(micro actuator)、微型馬達(micro motor)，微感測器(micro sensor)等微型元件，可應用在包括資 訊、通訊、電子、生醫及半導體等相關產業。

本文即探討在微機電系統中關於微流體系統的領域。微流體系統包含了微幫浦 (micro pump)、微流道(micro channel)、微閥門(micro valve)等，可整合成各種功能不同 的智慧型微流體晶片。在醫學的應用方面，可藉由微型的幫浦作精準微量的藥物輸送來 取代傳統的針孔注射，如此一來對於需接受化學療程的患者便能減輕其所承受的痛苦， 若能再配合微電子系統設備，更可達成定時定量藥物傳遞的功能，對於很多醫學治療方 法來講可完美的敘述微量藥物的流量且能程序化的控制是非常重要的。在冷卻應用方 面，需要靠微幫浦從低溫冷凍室抽取氦(He)或液態氮(N2)到熱交換區，以維持低溫的工 作環境，達成冷卻的效果，例如超導體只能在攝氏負兩百度以下的極低溫環境下工作， 此時就可將冷卻液和超導體合併在同一個矽晶片上，一面為超導體另外一面就為冷卻 液，透過微機電製程技術在矽晶片上製作出微小流道，再配合微幫浦去抽取冷卻液便能 有冷卻的功能，使超導體可穩定的運作。在能源應用方面，閥式壓電微幫浦可應用在甲 醇燃料電池上，用來傳遞燃料也可排除其化學變化後所產生的生成物。由此可知微流體 系統可應用的層面相當廣泛。

1.2 微幫浦簡介

在微流體系統的領域中，微幫浦常被應用在流體的輸送控制上，可算是系統中最重 要的一環，因此如何提高微幫浦的機械效率和精準控制微幫浦的流量是微流體系統中最 主要的研究目標。

一般根據微幫浦【2】的運作方式可分為機械式和非機械式兩類。常見的機械式幫 浦有閥門式幫浦(check-valve pump)、蠕動式幫浦(peristaltic pump)、無閥式幫浦(valve-less pump) 、 轉 動 式 幫 浦 (rotary pump) 、 超 音 波 式 幫 浦 (ultrasonic pump) 與 離 心 式 幫 浦 (centrifugal pump)。非機械式微幫浦有電液體力學驅動幫浦(EHD)、電磁式微幫浦 (MHD)、電滲式幫浦(electro-osmotic pump)、電泳式幫浦(electro-phoretic pump)、電濕式 幫浦(electro-wetting pump)與氣泡式幫浦(bubble pump)。

由上述可發現，非機械式的微幫浦主要是利用感應電場來牽引流體內的離子流動， 來達到推動流體的目的，雖然不需外加任何的可動元件，但因受到流體內離子多寡的影 響，因此無法廣泛的應用在各種流場而受到限制。另一大類為機械式微幫浦，它主要包 含了致動腔室和微閥門系統。腔室中有一振動薄膜，依據此一薄膜振動的動力來源可分 為氣動式(pneumatic)、熱氣動式(thermopneumatic)、熱驅動式(thermomechanic)、電磁式 (electromagnetic)、形狀記憶合金式(shape memory alloy)、靜電式(electrostatic)、壓電式 (piezoelectric)等數種。以下將依據薄膜振動的動力源，介紹三種常見的微幫浦。 — 熱氣動式微幫浦(thermopneumatic micropump)： 其原理為輸入電壓使加熱器因電阻而產生熱能，驅使致動腔室內的空氣膨脹造成薄 膜的變形；當不再供應電壓時，致動腔室內的空氣因不再受熱而自然冷卻，薄膜變會自 行回復原狀，如此反覆運作便可使薄膜產生振動，達到驅動主腔室內流體流動的效果。 — 靜電式微幫浦(electrostatic micropump)： 此微幫浦是在振動薄膜上施加電壓，造成靜電吸引力來使薄膜產生移動。當薄膜移 動使腔體內體積增加時，此時進口閥門打開出口閥門關閉，使流體流入腔體內，此為吸 水模式。當切斷薄膜上的電壓時，薄膜便會回覆原始狀態，此時就打開出口閥門關閉進

口閥門，使流體從出口流出，為排水模式。結合進水模式和排水模式即完成一個循環。 — 壓電式微幫浦(piezoelectric micropump)： 壓電式微幫浦是指振動薄膜是使用壓電材料製成。當壓電材料因體積變化而產生電 壓稱為正壓電效應；而若施加電壓於壓電材料使其產生變形則稱為負壓電效應。利用壓 電片的負壓電效應，便能藉由控制電壓而使壓電薄膜產生預期的振動，進而達到推動流 體的效果。 1.2.1 無閥式微幫浦簡介 無閥式微幫浦最大的特點在於利用漸縮/漸擴管(nozzle/diffuser)來取代傳統微幫浦 中閥門的設計，如此便能簡化微幫浦的構造與製程，減少機械元件耗損的可能性。其工 作原理可分為吸水模式和排水模式如(圖 1.1)來討論。 — 吸水模式(supply model) 振動薄膜向上移動使腔體體積增加，造成腔體内壓力小於進出口兩端，使流體分別

從進口由 diffuser 方向和出口由 nozzle 方向流進腔體，因 diffuser 方向的壓力損失小於

nozzle 方向的壓力損失，所以流進進口端的流量大於流進出口端的流量。 — 排水模式(pump model)

振動薄膜向下移動使腔體體積減少，造成腔體内壓力大於進出口兩端，使流體分別

從進口由 nozzle 方向和出口由 diffuser 方向流出腔體，因 diffuser 方向的壓力損失小於

nozzle 方向的壓力損失，所以從出口端流出的流量會大於從進口端流出的流量。經由此 週期性的運作，便會有一淨流量從出口端流出。 1.2.2 蠕動式微幫浦簡介 蠕動式微幫浦是由流道將多個腔體與輸入、輸出端串連起來所組成，且每個腔體內 都有一個致動薄膜，藉由控制腔內振動薄膜的運作順序，來達成運送流體往特定方向流 動的功能。二腔式的蠕動式微幫浦有二種薄膜振動相位，分別為三相位及四相位如(圖 1.2)；最常見的蠕動式微幫浦由三個腔體所組成，依薄膜振動順序，則主要可分為三相 位、四相位和六相位三種如(圖 1.3)。影響蠕動式微幫浦效能最重要的因素即為壓縮比ε

如(圖 1.4)【29】。 stroke volume ( ) dead volume ε = 薄膜掃過的最大體積 腔體體積 未被薄膜掃過 (1.1) 以下將對高壓縮比及低壓縮比的各式蠕動相位作比較。 — 高壓縮比： 在文獻中提到的蠕動式微幫浦，主要是把從入口處吸取的流體，藉由各腔體中薄膜 的蠕動行為，來達成傳遞流體的效果，其中因高壓縮比，振動薄膜有較大的振動幅度， 所以也扮演著被動閥門的角色，避免流體在腔體間，有回流的現象發生。若是薄膜的振 幅等同於腔體的深度，將類似一個開關，可完全阻隔各腔體間的流體，並可穩定地將定 量的流體傳遞在各個腔體之間，其中二腔三相位(圖 1.5(A))及三腔三相位(圖 1.6(A))即為 最好的例證。三腔六相位(圖 1.6(C))蠕動形式的發展，是希望能夠從入口端吸取更多的 流體進入蠕動式微幫浦的腔體內部，因此六相位比三相位多了步驟二，流進腔體內部的 流體量進而增加。六相位的步驟三及步驟六看似多餘且耗費時間，所以即發展出刪除步 驟三及步驟六的四相位(圖 1.6(B))蠕動形式，是故在高壓縮比的情況下，四相位所產生 的淨流量較大，六相位次之，三相位的淨流量最小。 — 低壓縮比： 在本文之中，因假設的薄膜振動頻率較高(2200 Hz)，故考量實際應用時，其產生的 振動幅度較小，有較低的壓縮比，無法扮演被動閥門的角色，因此在薄膜上下擺動時， 流體會在連接腔體與出入口端的漸縮/漸擴管間自由的進出，如(圖 1.7)、(圖 1.8)，藉由 漸縮/漸擴管的設計，使流體在每一週期出入腔體內部時，能有一向右往出口端的淨流量 產生。文獻中提出的三相位-本文中稱為順向三相位蠕動形式(圖 1.7(A))、(圖 1.8(A))， 在薄膜振動較小的情況下，每個漸縮/漸擴管間的流體，有三分之二週期的時間是往左側 入口端流動，僅有三分之一週期的時間是往右側出口端流動，造成流體往反方向輸送， 因此文中設計一逆向三相位的蠕動形式，其運作方式與順向三相位剛好相反，是故腔內 流體有三分之二週期的時間是往右側出口端流動，使能達到預期較大的流量值。而三腔 四相位(圖 1.8(B))及六相位(圖 1.8(C))為對稱的蠕動形式，各有半個週期的時間是分別流

往左側入口端及右側出口端，因此其產生的淨流量也較小。

1.3 文獻回顧

本章節主要以無閥式微幫浦的實驗與模擬、蠕動式微幫浦相關實驗及微幫浦理論分 析做相關文獻回顧。 — 無閥式微幫浦的實驗與模擬： 1990 年 Shoji et al.【3】在矽晶片上製作了兩種有閥門形式的微幫浦，其一為並聯 式的雙腔微幫浦，另一種為有類似緩衝存儲槽的串聯式雙腔微幫浦。實驗結果發現其在 背壓為1 mH2O 時，並聯式微幫浦的最大流量為 40 μl/min。

1993 年 Stemme E. and Stemme G.【4】第一個提出無閥門式幫浦的設計，其原理是

使用兩個圓錐狀的nozzle/diffuser 元件連接一個直徑為 19 mm 的腔體內含振動薄膜，實 驗後發現在頻率為100 Hz 下，最大流量為 16 ml/min，最大背壓可達 2 mH2O。 1995 年 Gerlach et al.【5】探討連結大角度漸縮管的壓電無閥式微幫浦流場情形， 其定義nozzle 流動方向為正向，diffuser 流動方向為反向。發現流場在雷諾數小於 15 時 為層流，當雷諾數大於100 時為紊流，並在 nozzle/diffuser 流場區域會有回流的產生。 1995 年 Olsson et al.【6】提出一個平面狀的無閥式微幫浦設計。此微幫浦並聯了兩 個直徑為13 mm 深度為 1 mm 的圓形腔體，每個腔體都各與兩個 nozzle/diffuser 連接。 經測試發現當此兩腔軆相位差為180°時，有最大流量 16 ml/min 及最大背壓 1.7 mH2O。 此篇研究也針對無閥式微幫浦作流量分析，且測量出nozzle 流向與 diffuser 流向在不同 壓力降下的流量差異。

1996 年 Olsson et al.【7】針對無閥式幫浦中 diffuser 元件的長度與張角，對幫浦效

率的影響作研究。文中把 diffuser 分為三段作討論，分別求出每段的壓力損失係數，並

與實驗結果相互比較與分析。經測量結果得出在紊流中 diffuser 元件的流向性與幫浦效

率比在層流中好。實驗採用diffuser 元件長度介於 1.45 到 3.95 mm 之間，而 diffuser 張

角介於1.9°到 6.8°之間。

幫浦，並在模擬薄膜振動時利用CFD-ACE 中的 Fluid-Thermo-Structure(FSI)模組，來處 理有關工作流體和振動薄膜間液固介面的問題。其模擬出的流量與實驗相當吻合，但在 最大背壓上則有明顯的誤差。 1998 年 Olsson et al.【9】利用熱壓法製造出無閥式微幫浦，其採用壓電材料作為薄 膜振動的動力源，且使用張角為7°的 diffuser 元件作為連接出入口端與腔體的流道。實 驗結果發現頻率在2200 Hz 下，單一腔體的無閥式微幫浦最大流量為 1 ml/min，最大背 壓為5.9 kPa。若是並聯兩個腔體的無閥式微幫浦，其最大流量可達 1.9 ml/min，最大背 壓為7.7 kPa。 2004 年 Yang et al.【10】製作漸縮/漸擴管(nozzle/diffuser)元件並測試與分析其壓力 損失係數與雷諾數的關係。結果發現壓力損失係數會隨雷諾數增加而遞減，在相同雷諾 數時，漸縮管的壓力損失大於漸擴管。漸縮/漸擴管的長度對於壓力損失的影響很小。當 漸擴管的張角大於 20°時，會有分離現象的產生，因此理論分析與實驗結果會有較大的 誤差。 2004 年羅卓錚【11】在其研究中利用 CFDRC 軟體模擬 Olsson 提出的無閥式微幫浦 設計。文中提到擋體式無閥門微幫浦，是指在靠近進出口兩端的流道中，設計一個梯形 擋體以代替nozzle/diffuser 的功能。模擬結果發現其最佳的擋體張角為 5°，但最大流量 與Olsson 實驗比較有 30%的誤差。 2005 年曾裕博【12】利用 Fluent 軟體模擬無閥式微幫浦，提出可用速度進口條件 代替薄膜振動時移動網格的設定，以簡化運算的難度與時間。研究中並測試 diffuser 元 件在不同雷諾數與張角下，其效率的差異。 — 蠕動式微幫浦相關實驗：

1968 年 Fung and Yih【13】最早提出蠕動行為對流場造成的影響。文中指出完全依 靠管壁上蠕動形式的移動，分析其對二維直管內部流場產生的效果。因蠕動行為而產生 的壓力梯度，如大於某一極限值，將會有回流產生，並探討各種雷諾數下，回流產生的 條件。

流場的相關性，文中利用無因次化的理論分析二維流場，並用實驗的結果作驗證。 1990 年 Smits【15】提出一個六相位蠕動式微幫浦的蠕動形式，並詳述其微幫浦腔 體是由矽晶片及玻璃薄膜並搭配壓電材料所組成。其微幫浦的最大流量可達 100 μl/min，最大背壓為 60 cmH2O。 2001 年 Cao et al.【16】提出一個能藉由幫浦內部體積變化，而達到傳遞流體功能 的蠕動式微幫浦，蠕動即為在其運作過程中，每個腔體內部的壓電薄膜都是各自在活動 的。文中並利用有限元素分析(FEA)的 ANSYS 軟體，模擬不同厚度及直徑下，壓電薄 膜所產生的最大振幅，與 MEMS 製程技術製作出的實際實驗結果作比對。其得到當壓 電薄膜厚度為40 μm，直徑為 12 μm 時有最大的振幅，與模擬結果相符。 2003 年 Berg et al.【17】製作出兩腔體的蠕動式微幫浦，並提出兩種蠕動形式做實 驗比較，且用簡易的理論分析預測其流量大小。微幫浦中每個腔體的體積為1.17 μl，有 較小的 dead volume。實驗結果得知類似三腔體中三相位的三相位蠕動形式其最大流量 為0.098 μl/s，兩階段式四相位蠕動形式最大流量為 0.107 μl/s。

2005 年 Teymoori and Abbaspour-Sani【18】設計一個靜電式的蠕動式微幫浦，其有 諸多優點如，與其他蠕動式微幫浦相比有較小的尺寸、採用被動式閥門取代傳統機械式 閥門較不易損壞、流量易控制及較小的能源消耗等。文中採用六相位的蠕動形式，壓縮 比約為0.8，微幫浦的尺寸為 7 mm × 4 mm × 1 mm。其淨流量可達 9.1 μl/min，與 ANSYS 模擬結果相符合。 2005 年 Goldschmidtboing et al.【19】建立一套塊狀理論模型，分析壓電型蠕動式微 幫浦的流場，依照其模擬後的結果，製造理想的微幫浦，微幫浦中的腔體深度為80 μm， 薄膜所能達到的最大振幅逾10 μm。文中並提到兩種腔體的蠕動形式，分別為四相位及 六相位，再根據不同的振動頻率及背壓，實際量測出微幫浦的淨流量與模擬作比對，結 果相符。而此壓電型蠕動式微幫浦有最大淨流量為1.4 ml/min，最大可承受的背壓為 40 kPa。 2006 年關恕【20】以 MEMS 技術製作出三腔壓電蠕動式微幫浦，腔體的深度為 30 μm，並實驗比較不同的驅動相位，當腔內壓電薄膜隨頻率增加時，其最大位移量變化

的情形，其薄膜的振幅介於0.02 μm~1.8 μm 之間。亦改良電壓輸入的迴路設計，並有效 提升微幫浦的淨流量與效率。

2008 年 Jang and Yu【21】製作了串聯三個腔體的蠕動式微幫浦，實驗比較各腔體 中薄膜振動的順序對微幫浦流量的影響。三腔蠕動式微幫浦依其作動的方式，可分為三 相位、四相位和六相位共三種，在供應100 伏特的電壓下，四相位蠕動形式的微幫浦有 最大流量為17.6 μl/min，其微幫浦內腔體深度為 10 μm，壓電薄膜有最大位移量為 2.91 μm。 — 微幫浦理論分析： 1998 年 Ullmann【22】利用流量與壓力的關係式，分析單腔壓電無閥式微幫浦的流 場，並把微幫浦的流場情形分為三種模式作運算，分別求出每個時間點下的流量與壓 力。文中也提及串聯雙腔無閥式微幫浦的分析，其所需考慮的模式較多。 1999 年 Olsson et al.【23】提出塊狀數學模型來分析無閥式微幫浦。其理論是把微 幫浦分為數個區塊作討論，並利用連續方程式和伯努力方程式整合出三個方程式系統， 再用MATLAB 解出其流量與壓力，模擬結果與實驗數據相符。 2001 年 Pan et al.【24】找出壓電薄膜振動的統御方程式，且探討薄膜振動最適合的 頻率範圍。在研究中把微幫浦的流場情形分為三個模式作計算，再藉由漸縮/漸擴管的壓 力損失係數代入壓力的關係式中，估算出微幫浦的流量。 2003 年 Pan et al.【25】將流體的流動與壓電薄膜振動的統御方程式作偶合，並考慮 流體慣性對流場造成的影響，以求得流量。 — 微幫浦計算模擬與理論分析比較： 2008 年吳欣恩【26】在其微幫浦的理論分析中，提出合併型的壓力與流量關係式， 在計算單腔微幫浦流場時不必再分為吸水模式、排水模式與過度區等多種模式，簡化數 學公式，且解決探討多個腔體時將遇到畫分模式的問題。文中在lump system 理論中加 入流體的慣性項作計算後，與有限體積法的CFD 運算模擬結果相似。

2008 年 Tsui and Lu【27】藉由有限體積法的 CFD 計算模擬單腔無閥式微幫浦的流 量變化情形，其在薄膜振動上使用往復式的速度邊界條件，取代移動網格的設定，薄膜

振動的位移則融合二次曲線與梯形曲線，計算在不同背壓下淨流量的數值與實驗作比 較。文中並利用壓力與速度的關係式，提出簡易的數學理論，分析微帮浦的流量與壓力 變化，且提出三種計算微幫浦的效率公式，並繪出不同背壓時微幫浦的效率變化。

1.4 研究目的

隨著微機電系統(MEMS)科技迅速的發展，人們可利用其微加工的技術，製作出功 能齊備但尺寸更小的各項機械元件，利用這些微米化的元件便可作出高效率、低耗能且 小型化的電子產品，並廣泛地應用在週邊的生活之中。其中具有傳輸流體功能的幫浦， 亦得利於微機電系統技術的進步，使其體積縮小便於攜帶，在內部安裝電子化的微致動 器後，亦可藉由電腦作程序化的控制及準確的傳遞流體，使得微幫浦可應用在醫學藥物 傳送、生物化學試劑檢測等需要精確控制流量的應用層面，因此近年來吸引很多學者研 究探討關於微幫浦的製造與分析。 微幫浦效能的優劣，可藉由兩方面作分析比較，其一是微幫浦運作一個週期內所能 傳遞流體的淨流量大小；另一方面是在微幫浦的出口端施加一背壓，觀察期週期內淨流 量所受的影響，當淨流量值為零時，其出口端的背壓值，即為該微幫浦最大所能承受的

背壓。本文中串聯二至三個單腔無閥式微幫浦，並施加Berg et al.【17】與 Jang and Yu

【21】整理提出的各種二腔及三腔蠕動相位形式，藉由 Tsui and Lu【27】與吳欣恩【26】 提出的塊狀理論分析，可快速地計算出微幫浦一週期內的淨流量，及腔內壓力與流量的 變化關係，並利用有限體積法模擬計算結果作比對。期望找出一蠕動相位形式，能提升 低壓縮比微幫浦的工作效率，包含有較高的淨流量與可承受的背壓，並探討增加腔體數 目時，在相同的蠕動相位下，微幫浦的效能是否能隨之提升。

第二章 數學模式

本章主要介紹求解流場時所需用到的數學模式。首先對模擬的流場作基本假設，接 著敘述求解流場時所需用到的統御方程式，以及近似壓電薄膜振動的曲線，最後對於邊 界條件的設定作探討。

2.1 基本假設

本文中模擬壓電二腔及三腔微幫浦的幾何形狀，及詳細的尺寸可參考(圖 1.9)。使用 水作為微幫浦的工作流體。以下對模擬的流場作基本假設。

z 三維流場(three-dimensional flow)、忽略重力項(neglect body force)。 z 不可壓縮流(incompressible flow)：工作流體-水的密度為定值。 z 非穩態流場(unsteady flow)：腔體中的壓電薄膜振動隨時間變化，因此腔內流量也 為時間的函數。 z 層流(laminar flow)：本文模擬的微幫浦尺寸較小，流速常介於 10-5_到10-6 _{m/s 之間} 較為緩慢，故雷諾數小，假設幫浦內部的流場為層流。

2.2 統御方程式

連續方程式(continuity equation)：

∇

⋅

(

ρ

V

G

)

=

0

(2.1) 動量方程式(momentum equation)：

( ) ( ) ( )

V VV V P t ρ ρ μ ∂ + ∇ ⋅ = ∇ ⋅ ∇ − ∇ ∂ JG JGJG JG (2.2) ρ 為水的密度；μ 為水的黏滯係數。

2.4 邊界條件

2.4.1 對稱邊界條件 由於對稱的幾何形狀，僅模擬真實微幫浦的一半，故幫浦數值計算時採用對稱邊界。

2.4.2 壓電薄膜邊界條件 實際壓電薄膜振動的曲線為一高階的偏微分方程式，為了簡化計算的困難度，因此 薄膜的位移假設為一梯形曲線如圖(2.1)，隨時間呈現一簡諧運動(harmonic motion)。 梯形曲線(trapezoidal curve)位移函數： = ) , ( tr Z 0 0 1 1, cos(2 ) MAX r r d Min ft r r π ⎛ − ⎞ − ⋅ _{⎜ ⎟} − ⎝ ⎠ (2.3) 上式中r 為腔體半徑(₀ r = 3 mm)，₀ r 為壓電薄膜的半徑(₁ r = 2 mm)，₁ d_MAX為壓電薄膜從水 平位置到上死點的最大位移量(d_MAX = 1 μm)。壓電薄膜振動的頻率假設為 f = 2200 Hz。 壓電薄膜會隨時間的變化而產生位移，因此邊界需設定為移動邊界，並利用移動網 格作計算，但會增加計算量耗費時間。由於薄膜振動幅度最大僅為 1 μm，相較於幫浦 腔體高度0.2 mm，其比值為 1：200，因此可於數值模擬中假設薄膜位置固定，取而代 之給定薄膜運動之速度作為邊界條件。將(2.3)式微分得 梯形曲線(trapezoidal curve)運動函數： 0 0 1 ( , ) 2 _MAX 1, r r sin(2 ) V r t f d Min ft r r

π

⎛ − ⎞

π

= ⋅ ⋅ _{⎜ ⎟}⋅ − ⎝ ⎠ (2.4) 壓電薄膜振動為梯形曲線的形式如(2.3)式，水平向上掃過的最大體積∀MAX可計算如下： 1 0 1 0 max 0 1 2 2 0 0 0 0 0 1 MAX MAX r r MAX MAX r r r d dA d dA r r r r d rdrd d rdrd r r π π θ θ − ∀ = + ⋅ − − = ⋅ + ⋅ −

∫∫

∫∫

∫ ∫

∫ ∫

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ _{− +} − ⋅ + ⋅ = 3 1 2 1 0 3 0 1 0 2 1 3 1 2 1 6 1 2 r r r r r r d r d MAX MAX π π (2.5) 代入各項數值後，可計算出薄膜由水平向上掃過的最大體積為1.9897×10-11 _m3_，而薄膜 由最高點到最低點掃過的體積是其兩倍為3.9794×10-11 m3。又每個腔體的體積為 5.6549×10-9 _m3_{，藉由計算壓縮比的(1.1)式，可得知壓縮比ε = 0.00706，故本研究探討} 的蠕動式微幫浦屬於低壓縮比。

2.4.3 出入口端邊界條件

本文中微幫浦的進出口兩端設定為壓力邊界條件，並考慮在進出口端有壓力差時， 微幫浦的流量變化。

2.4.4 壁面邊界條件

第三章 方程式離散化

3.1 有限體積法

早期流體力學的發展以實驗觀察與理論分析為主，隨著數值方法的進步和計算機運 算處理速度的精進，人們開始試著以數值模擬搭配實驗數據求得最佳化的成果。用來求 解工程應用上問題的數值方法有很多種，如有限體積法(Finite Volume Method)、有限元 素法(Finite Element Method)、有限差分法(Finite Difference Method)等。其中有限體積法 是把微分方程式對某個控制體積作積分，經過整理後得出一組離散化的代數方程式，此 種方法適合求解通量等物理問題，因此本文數值模擬採用有限體積法求解動量方程式 (momentum equation)，並應用於非結構性網格(如圖 3.1)。

3.2 離散動量方程式

由動量方程式(2.2)式改寫成下列通式：

( )

( )

V

( )

q t +∇⋅ = ∇∇ + ∂ ∂ ρφ _{ρ φ μ φ} (3.1) 將上式對一控制體積作積分後，表示如(3.2)式。其分別由非穩態項(unsteady term)、對流 項(convection term)、擴散項(diffusion term)及源項(source term)組成。

∫∫∫

+

∫∫∫

∇⋅ =

∫∫∫

∇⋅ ∇ +

∫∫∫

∂ ∂ V V V V dV q dV dV V dV t

ρφ

(

ρ

φ

) (

μ

φ

) (3.2) 上式中φ為欲求的流場性質，在本文中即為速度向量，源項因忽略了重力影響只剩壓力 項，故q=−∇P。

經由高斯散度定理(Gauss Divergence Theorem)，將(3.2)式中對流項與擴散項的體積分轉 換為面積分，並整理如下：

_∂

∂

∫∫∫

+

∫∫

⋅

=

∫∫

∇

⋅

+

∫∫∫

V S S V

dV

q

S

d

S

d

V

dV

t

ρφ

ρ

φ

μ

φ

(3.3)

3.2.1 非穩態項 在非穩態項中，使用一階近似作離散，可表示如下：

(

o

)

V t V dV t φ φ ρ ρφ − Δ Δ ⋅ = ∂ ∂

∫∫∫

(3.4) 上式中φ為此刻時間點上的流場性質，φo 為前一個時間點上的流場性質，故φo 為已知 可移到等號右側的源項作計算。 3.2.2 對流項 (3.3)式中的對流項，經過高斯散度定理轉換為面積分，再經離散化整理後可表示為 下式：

∫∫

⋅ =

∑

(

⋅

)

=

∑

=

∑

f C f f f f f f f f f S F m S V S d V

φ

ρ

φ

φ

ρ

G ) G G G  ( (3.5) 上式中，Sf G 為面之法向量，m 為通過任一面的質量流率，f C f F 為對流通量。下標符號 f 為控制體積中任一面上所代表的點，對於在其面上的值，可以藉由相鄰網格以線性內插 的數學式求得。 求解對流項時採用一階上風法(Upwind scheme)與中央差分法的混合型：

( ) (

fUD

)

CD f UD f C f

F

F

F

F

=

+

γ

−

(3.6) 上式中，γ 的數值介於 0 到 1 之間，γ =0為一階上風法，γ =1為中央差分法。 UD f F 為 一階上風法求得的對流通量， CD f F 為中央差分法求得的對流通量。如果完全採用一階上 風法，會產生很大的數值擴散誤差較大，只採用中央差分法會造成震盪不易收斂，故本 文令γ =0.9，使數值易收斂又較為精確。 — 一階上風法，φf 為其上游的數值：

(

)

(

)

0 0 f P f f C f m m φ φ φ φ ⎧ = > ⎪ ⎨ = < ⎪⎩   ， ， (3.7) 上式中，下標 P 為主格點，下標C為外部的相鄰格點。

— 中央差分法，φf 為相鄰網格的內差值： C P f wφ wφ φ =(1− ) + (3.8) 上式中，w為權重係數，其w值介於0 到 1 之間。 將上述(3.6)式、(3.7)式、(3.8)式整理後可得：

[

]

(

)

[

]

[

]

{

f P C f P f C

}

C f P f C f m m w w m m m F φ φ φ φ γ φ φ ) 0 , max( ) 0 , max( 1 ) 0 , max( ) 0 , max(      − + − + − + − + = (3.9) 上式中，與γ 相乘大括弧內的算式為前一個時間點運算後的數值，故可移到源項作計算。 3.2.3 擴散項 (3.3)式中右側的擴散項，經過高斯散度定理轉換為面積分，再經離散化整理後可表 示為下式：

∫∫

∇

⋅

=

∑

∇

⋅

=

∑

f D f f f f f S

F

S

S

d

G

μ

(

φ

)

G

φ

μ

(3.10) 上式中，μ_f 為面上的黏滯係數， D f F 為擴散通量。 對於非結構性網格，擴散項中面的法向量採用 over-relaxed 的方法作處理。原先面 的法向量SG_f 可表示成下式：

S

_f

=

d

+

(

S

_f

−

d

)

(3.11) 上式中dG為主格點 P 與相鄰格點C相連的向量，定義如下： PC f PC f

S

S

d

δ

δ

G

G

G

G

G

⋅

=

2 (3.12) 上式中，δG_PC(如圖 3.2)為主格點 P 到相鄰格點C的距離向量。 將上述(3.10)式、(3.11)式、(3.12)式整理後可得：

(

)

(

)

2

d

S

S

S

F

_{C P f f f} f PC f f D f

G

G

G

G

G

−

⋅

∇

+

−

⋅

=

φ

φ

μ

φ

δ

μ

(3.13) 上式中，等號右側第二項為前一個時間點運算後的數值，故可移到源項作計算。

3.2.4 源項

(3.3)式中右側第二項的源項，經過高斯散度定理轉換為面積分，再經離散化整理後 可表示為下式：

∫∫∫

∫∫

∑

Δ = Δ = ∇ Δ f f f S V S P V S d P V dV P V G G ₁ 1 1 (3.14) 在數值模擬中，壓力梯度為已知，故源項即為壓力梯度。但是邊界上的壓力梯度需 經由下列式子計算：

P

_b

−

P

_P

=

∇

P

⋅

δ

K

_PC (3.15) 上式中，P_b為壁面上的壓力，δG為主格點 P 與壁面b點相連的向量。 面上的所有壓力梯度可由下式表示： _⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + Δ = Δ = ∇

∑

∑

≠ b f f f b b f f f P S P S V S P V P 1 G 1 G G (3.16) 整理後可得邊界上的壓力P_b： ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ _⋅ Δ − ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⋅ Δ + =

∑

≠ PC b b f PC f f P b S V S P V P P δ δ G G G G 1 1 1 (3.17)

3.3 代數方程式

將 3.2 節中離散化後所得到的非穩態項、對流項、擴散項和源項整理之後，可求得 下列之線性代數方程式：

A

_P

φ

_Pn

=

∑

A

_C

φ

_C

+

Q

(3.18) 上式中，A_P 、 A_C 和Q 分別為：

t

V

A

A

_{P C}

Δ

Δ

⋅

+

=

∑

ρφ

(3.19) max( ,0) 2 f f PC f f C m S S A G G  G − + ⋅ =

δ

μ

(3.20)

(

)

[

]

[

]

{

}

t V V P d S m m w w m Q o f f f C f P f C P f Δ Δ ⋅ + Δ ⋅ ∇ + − ⋅ ∇ − − + − + − =

∑

μ φ ρ φ φ φ φ φ γ ) ( ) 0 , max( ) 0 , max( 1 G G    　　　　　　 (3.21) (3.20)式中 AC 項由部分的對流項與部分的擴散項所組成。(3.21)式中Q 包含非穩態 項、對流項、擴散項和源項前一個時間點的疊代值。

第四章 PISO 演算法

在上個章節中，利用已知的壓力梯度求解動量方程式得到速度場，但因為未計算連 續方程式，故不滿足質量守恆。為了滿足質量守恆定律，進而發展了各種壓力與速度偶

合的關係式，並對壓力與速度作修正。現今常見的壓力與速度偶合演算法有SIMPLE、

SIMPLEC、PISO 等，本文採用 PISO(Pressure Implicit with Splitting of Operators)，此為

一種先預測後再修正之演算法則，以下為對PISO 演算步驟之描述。

4.1 壓力與速度修正

4.1.1 預測步驟 首先將(3.18)式源項 Q 中的壓力項提出，並將速度場代入φ可得： A_PVG_P* =

∑

A_C VG_C*+

(

Q′−∇P_P*⋅ΔV

)

(4.1) 上式中， Q′ 為不包含壓力項的源項。 利用現有之壓力 * P ，求解動量方程式可求得速度VG*。 4.1.2 第一次修正 上節中求得了速度 * VG 與壓力 P ，接下來敘述如何修正此速度及壓力。 將(4.1)式整理簡化如下所示： * * * P P P P P P A V H V _⎟⎟ ∇ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ − = G G (4.2) 上式中， P C C P A Q V A H =

∑

+ ′ * * G G 。 控制體積面上的速度依(4.2)式的形式可近似如下： * * * f f f p _f V V H P A ⎛_Δ ⎞ = −⎜_⎜ ⎟_⎟ ∇ ⎝ ⎠ G G (4.3)

上式中， f f P f f P A V V H _⎟⎟ ∇ ⎠ ⎞ ⎜⎜ ⎝ ⎛Δ + =G G_* 代回可得： * * f f f f P _f p _f V V V V P P A A ⎡ _{⎛Δ ⎞ ⎤ ⎛Δ ⎞} =⎢ +_{⎜ ⎟} ∇ ⎥ ⎜−_⎜ ⎟_⎟ ∇ ⎢ ⎝ ⎠ ⎥ ⎝ ⎠ ⎣ ⎦ G G (4.4) 上式中， f f A V ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Δ 為主格點與相鄰格點的平均值，即下列所示： ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ + ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ C P P P f P A V A V A V 2 1 _(4.5) 故質量流率可寫為： d P P A V S V S P P A V S V S V m f f f P f f f f f f f f P f f f f f f f f ⋅ ∇ − ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ − ⋅ ≈ ⋅ ∇ − ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ − ⋅ = ⋅ = ) ( ) ( * * * * ρ ρ ρ ρ ρ G G G G G G G 

(

)

2 f f f f f f C P PC P PC f S V V S P P P A S ρ ρ δ δ ⎛ ⎞ Δ ⎜ _{⎟ ⎡ ⎤} = ⋅ − _{⎜ ⎟ ⎣} − − ∇ ⋅ _⎦ ⋅ ⎜ ⎟ ⎝ ⎠ G G G G G G (4.6) 依照(4.2)式的形式，修正後的速度與壓力的關係可表示為下式： ** * ** P P P P P P A V H V _⎟⎟ ∇ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ − = G G (4.7) 上式中，VP =VP +VP′ G G G_{** *} ；PP =PP +PP′ * * * _{，上標 * 為未修正的值，上標 ** 為第一次修正} 後的值，上標 ' 為第一次修正的修正量。 將(4.7)式與(4.2)式相減可得： _P P P P P A V V _⎟⎟ ∇ ′ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ − = ′ G (4.8) 控制體積面上的速度與壓力關係式可表示如下：

f f P f P A V V _⎟⎟ ∇ ′ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ − = ′ G (4.9) 第一次修正的質量流率可表示如下：

m



*_f*

=

m



*_f

+

m



′

_f (4.10) 上式中，m′_f =

ρ

_fVG_f′⋅SG_f 。 修正後的質量流率可計算如下：

(

S

d

)

P

A

V

d

P

A

V

m

S

P

A

V

m

S

V

m

m

f f f P f f f P f f f f f P f f f f f f f

G

G

G



G



G

G





−

⋅

′

∇

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ Δ

−

⋅

′

∇

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ Δ

−

=

⋅

′

∇

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ Δ

−

=

⋅

′

+

=

ρ

ρ

ρ

ρ

* * * * *

(

)

P

(

S

d

)

A

V

P

P

S

S

A

V

m

_{f f} f P f P C f f f P f f

G

G

G

G

G



_⎟⎟

∇

′

⋅

−

⎠

⎞

⎜⎜

⎝

⎛ Δ

−

′

−

′

⋅

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ Δ

−

=

ρ

δ

ρ

2 * (4.11) 為了滿足質量守恆，故控制體積面上的淨流量為0，表示如下：

∑

**

=

∑

*

+

∑

′

f

=

0

f f f f

m

m

m







(4.12) 整合(4.11)式與(4.12)式，並用代數形式表示，可得修正式如下： ′ =

∑

′ + + f P P C C P PP A P S S A _{1 2} (4.13) 上式中， f f f P f C S S A V A G G G ⋅ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ = δ ρ 2 ； = −

∑

f f P m S ₁  ； *

∑

_⎟⎟

∇

′

⋅

−

⎠

⎞

⎜⎜

⎝

⎛ Δ

=

f f f f P f P

P

S

d

A

V

S

₂

(

)

G

G

ρ

。 4.1.3 第二次修正 利用第一次修正求得壓力 ** P 與速度VG**，做第二次的壓力與速度修正，其方式和第

一次壓力與速度修正的形式類似，故第二次修正，主格點的速度與壓力關係式可表示如 下： * * * * * * * * P P P P P P A V H V _⎟⎟ ∇ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ − = G G (4.14) 將(4.14)式與(4.7)式相減可得： P P P P P P A V H V _⎟⎟ ∇ ′′ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ − ′′ = ′′ G G (4.15) 第二次修正的質量流率表示如下： f f f m m m*** = **+ ′′ (4.16) f

P f C C f f f f P f f

S

A

V

A

S

P

A

V

m

G

G

G



⋅

′′

+

⋅′′

∇

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ Δ

−

=

′′

∑

⋅

ρ

ρ

(4.17) 結合(4.16)式與(4.17)式的修正後質量流率可表示如下： *** **

(

)

f P f C C f f f f P f f f P f f f S A V A d S P A V d P A V m m G G G G G   ⋅ ′′ + − ⋅′ ′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ − ⋅′ ′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ − =

∑

⋅ ⋅ ρ ρ ρ (4.18) 為了滿足質量守恆，故控制體積面上的淨流量為0，表示如下：

∑

***

=

∑

**

+

∑

′′

f

=

0

f f f f

m

m

m







(4.19) 結合(4.12)式與(4.19)式可得：

∑

***

=

∑

′′

f

=

0

f f

m

m





(4.20) 整合(4.18)式與(4.20)式，並用代數形式表示，可得修正式如下： ′′=

∑

′′+ + f P P C C P PP A P S S A _{1 2} (4.21) 上式中， f f f P f C S S A V A G G G ⋅ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ Δ = δ ρ 2 ； f P f C C f P S A V A S G G ⋅ ′′ − = ρ

∑

1 ；

∑

_⎟⎟

∇

′′

⋅

−

⎠

⎞

⎜⎜

⎝

⎛ Δ

=

f f f f P f P

P

S

d

A

V

S

₂

(

)

G

G

ρ

。 4.1.4 修正後方程式整理 為了使計算結果更為準確，在求解壓力修正式的過程中，將源項分為兩部份作運 算，第一部份源項只有 S ；第二部份源項只有_P1 S_P2。表示如下： 正交修正(orthogonal corrector)

∑

′+ = ′ f P C C P PP A P S A 1 (4.22) 非正交修正(non-orthogonal corrector)

∑

′ + = ′ f P C C P PP A P S A ₂ (4.23)

4.2 邊界條件設定

4.2.1 進出口邊界條件之流量計算 本文進出口邊界設定為壓力邊界，故在進出口端給定一壓力值，並以此修正進口與 出口的質量流率，其步驟如下所示： 步驟一：在壓力邊界相鄰的網格點上設定一給定的壓力值。 步驟二：求解兩次修正式的過程中，需設法將邊界相鄰網格點的壓力修正量為 0(如 圖4.1)，以滿足P′=0。故將壓力修正式的係數與源項作調整，表示如下： ′ =

∑

′+ + f P P C C P PP A P S S A 1 2 (4.24) 上式中，S_P1 = S_P2 = A_C =0。 步驟三：求解完壓力修正式後，需使相鄰的網格滿足質量守恆：

∫∫

⋅ =

∑

=0 f f S m S d V  ρ (4.25) 故邊界上的質量流率可表示為如(圖 4.1)： mf_,B =−

(

mf_,1+mf_,2 +mf_,3

)

(4.26)

4.2.2 進出口邊界條件之速度計算 壓力邊界出口端的速度，是藉由對流邊界條件(Convective Boundary)計算求得，其 數學形式為下列所示： =0 ∂ ∂ + ∂ ∂ z u t c

φ

φ

(4.27) 離散化(4.27)式可得： =0 Δ − + Δ − z u t o C n B c o B n B

φ

φ

φ

φ

(4.28) 上式中，上標n 為此刻時間點上的值，上標 o 為前一個時間點的值。下標 B 為邊界的網 格，下標C 為與邊界相鄰的網格。 將(4.28)式整理成下列所示：

Cr

Cr

_co o B n B

₊

+

=

1

φ

φ

φ

(4.29) 上式中，Cr 為 Courant number， z t u Cr c Δ Δ = 。 4.2.3 固定牆邊界條件 固定牆上採用無滑移邊界條件(no-slip condition)，VG =0。

4.3 PISO 演算步驟

第一步：給定邊界的初始值。 第二步：修正進口邊界上的速度與流量。 第三步：解動量方程式，求得速度 * VG 。 第四步：求解第一次壓力修正得到 P′ ，並利用 P′ 求出第一次修正後的速度 ** VG 、壓力P 、質量流率** m*f*。 第五步：求解第二次壓力修正得到 P ′′ ，並利用 P ′′ 求出第一次修正後的速度 *** VG 、壓力P***、質量流率m*_f**。

第六步：計算壓力邊界上的流量

第五章 微幫浦塊狀理論分析

5.1 塊狀系統理論分析(Lump-system analysis)

1998 年 Ullmann【22】提出一個數學模式，分析無閥式微幫浦的流場，並可估算出 腔體內的壓力變化及進出口的流量大小。其數學模式是依照流體流動的方向，將單腔無 閥式微幫浦一個週期的運作，分為三種模式作計算，分別為排水模式、吸水模式與過渡 模式，並定義流體流動方向向右為正。以下對其數學模式做介紹： 流體流經漸縮/漸擴管時，會因其幾何外形而造成壓力損失，其結合速度V 與壓力差 P Δ 可表示成： K V2 =ΔP 2 1 _ρ (5.1) 將其改寫成流量 Q 與壓力差 PΔ 的形式為： Q P A K t Δ = 2 2 2 ρ (5.2) 上式中，A 為漸縮/漸擴管的喉部面積。K 為壓力損失係數，流體流向漸縮管(nozzle)時，_t n K K = ；流體流向為漸擴管時，K =Kd。 — 排水模式(P_c > P_out >P_in)：此模式中腔體內壓力大於進出口端壓力，流體經漸縮管 (nozzle)方向由腔體流到進口端；經漸擴管(diffuser)方向由腔體流到出口端。 in c in t n Q P P A K ₂ 2 = − 2 ρ (5.3A) out c out t d Q P P A K ₂ 2 = − 2 ρ (5.3B)

(

ft

)

f Q

Q_in + _out =∀_MAX ⋅2π sin 2π

− (5.3C) — 吸水模式(Pout >Pin >Pc)：此模式中進出口端壓力大於腔體內壓力，流體經漸擴管 (diffuser)方向由進口端流進腔體內部；經漸縮管(nozzle)方向由出口端流進腔體內 部。 in in c t d Q P P A K ₂ 2 = − 2 ρ (5.4A)

c out out t n Q P P A K ₂ 2 = − 2 ρ (5.4B)

(

ft

)

f Q

Qin − out =∀MAX ⋅2π sin 2π (5.4C)

— 過渡模式(P_out >P_c >P_in)：此模式中流體經漸縮管(nozzle)方向由腔體流到進口端； 經漸縮管(nozzle)方向由出口端流進腔體內部。 in c in t n Q P P A K ₂ 2 = − 2 ρ (5.5A) c out out t n Q P P A K ₂ 2 = − 2 ρ (5.5B)

(

ft

)

f Q

Q_in − _out =∀_MAX ⋅2π sin 2π

− (5.5C) 上述方程式中，∀MAX為壓電薄膜水平向上掃過的最大體積。 f 為振動頻率。 5.1.1 單腔無閥式 從上述的方程組中可發現，Ullmann【22】提出的微幫浦理論，需依照流體流動的 方向，將方程式劃分為數種模式作討論，當微幫浦腔體數目增加時，就會增加劃分模式 的困難度，計算上也較為複雜。因此吳欣恩【26】提出一種方式，可將不同模式下的方 程式作合併，使其理論可套用在多腔體的微幫浦中。 吳欣恩【26】發現在(5.9)式到(5.11)式中，流體流動的方向會影響壓力損失係數K 與n d K 的選取，因此藉由流量正負變化的關係與 max 和 min 的指令，使其自動選取適合的K_n 或K 值作運算，達成將三個方程組合併成為一個方程組的目的，其中d Q 與in Qout為前一 個時間點的數值。 壓力損失係數與流量的關係，採用呂學霖【28】藉由模擬漸縮/漸擴管元件，Curve Fitting 出的關係式，表示如下： ₌1.315_×10−7 _⋅ −0.921₊0.5981 Q K_d (5.6) ₌1.173_×10−6_⋅ −0.8112 ₊1.204 Q K_n (5.7) 上列式子中， Q 為前一個時間點的數值。 文獻【26】提出微幫浦中流體流動的情形隨時間而變化，為非穩態的形式，故考慮

流體慣性項的影響，在穩態形式中加入慣性項，使分析後的流場情形能更貼近實際的物 理現象。其方程組表示如下：

(

)

dt P P dQ A A l Q Q A K Q Q K Q Q in in l t Q in in t n in in d in in = − + + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − 2 1 2 ) 0 , min( ) 0 , max( ρ₂ α ρ (5.8A)

(

)

out out l t Q out out t n out out d out out P P dt dQ A A l Q Q A K Q Q K Q Q − = + + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − 2 1 2 ) 0 , min( ) 0 , max( ρ₂ α ρ (5.8B)

(

ft

)

f Q

Qout − in =∀MAX ⋅2π sin 2π (5.8C)

上列式子中，A 為漸縮/漸擴管的擴張部面積；_l l為漸縮/漸擴管的長度。當α_Q =0為穩

態；當αQ =1為非穩態。

求解慣性項時採用全隱式差分法(fully implicit scheme)，表示如下： t Q Q dt dQ o Δ − = (5.9) 上式中， o Q 為前一個時間點的數值。 在無背壓時，初始條件(initial condition)的設定如下： s m t f f Q_{in MAX}2 sin(2 ) 2.16 10 9 3 2 1_{∀ Δ = ×} − = π π s m t f f

Q_{out MAX}2 sin(2 ) 2.16 10 9 3

2 1_{∀ Δ = ×} − = π π Pa P=0 上式中， f n t × = Δ 1 。n為壓電薄膜上下振動一個週期中劃分的 time step。在本文中設 400 = n ， f =2200為振動頻率。故單腔無閥式微幫浦運作一個週期所需的時間為 1 T f = 。

藉由上述的聯立方程組，求解在Steady Model 與 Unsteady Model 下，Q 、in Qout、P

在不同time step 的數值，並計算出微幫浦一個週期的淨流量。

5.1.2 二腔蠕動式

而Q 代表流經中間段漸縮/漸擴管(nozzle/diffuser)或直管(straight tube)的流量，依照單腔c 無閥式的形式將求解的方程組表示如下：

(

)

1 2 2 1 2 ) 0 , min( ) 0 , max( P P dt dQ A A l Q Q A K Q Q K Q Q in in l t Q in in t n in in d in in = − + + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − ρ α ρ (5.10A)

(

)

1 2 2 2 1 2 ) 0 , min( ) 0 , max( P P dt dQ A A l Q Q A K Q Q K Q Q _c l t Q c c t n c c d c c = − + + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − ρ α ρ (5.10B)

(

)

out out l t Q out out t n out out d out out P P dt dQ A A l Q Q A K Q Q K Q Q − = + + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − 2 2 2 1 2 ) 0 , min( ) 0 , max( ρ α ρ (5.10C) 1 Q Q Qc − in =Δ (5.10D) 2 Q Q Q_out − _c =Δ (5.10E) 上述的方程組中，下標1代表第一個腔體；下標 2 代表第二個腔體。 QΔ 為腔體因壓 電薄膜振動造成內部體積改變的變化率，其變化的情形因蠕動形式的不同而有差異。文 獻【17】提出二種雙腔蠕動式微幫浦的蠕動形式如(圖 1.2)分別為三相位與四相位。本文 中提出逆向三相位的蠕動形式作比較，原始的三相位在文中稱為順向三相位，蠕動形式 如(圖 1.7)。蠕動式微幫浦運作一個週期所需的時間為τ ， f a n a nf 2 2 1 _{× =} = τ ，a為相位 數。其不同的蠕動形式所造成的體積變化率表示如下： z 2 腔-順向三相位

(

)

(

)

⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ < < ⋅ ⋅ ∀ − < < < < ⋅ ⋅ ∀ = Δ τ τ π π τ τ τ π π t t f f t t t f f Q MAX MAX 3 2 if , 2 sin 2 3 2 3 if , 0 3 0 if , 2 sin 2 1

(

)

⎪⎩ ⎪ ⎨ ⎧ < < < < ⋅ ⋅ ∀ − = Δ τ τ τ π π t t t f f Q MAX 3 2 if , 0 3 2 0 if , 2 sin 2 2 3 3 2f 2T τ = = 。 z 2 腔-逆向三相位