雙腔無閥式微幫浦流場理論分析

國 立 交 通 大 學

機 械 工 程 研 究 所

碩 士 論 文

雙腔無閥式微幫浦流場理論分析

Analysis of the flow in double-chamber valveless micropumps

研 究 生：吳 欣 恩

指 導 教 授：崔 燕 勇 博士

雙腔無閥式微幫浦流場理論分析

Analysis of the flow in double-chamber valveless micropumps

研 究 生：吳欣恩 Student：Shin-En Wu

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

國 立 交 通 大 學

機械工程研究所

碩士論文

A Thesis

Submitted to Institute of Mechanical Engineering

Collage of Engineering

National Chiao Tung University

In Partial Fulfillment of the Requirements

For the degree of

Master of Science

In

Mechanical Engineering

July 2008

雙腔無閥式微幫浦流場理論分析

研究生：吳欣恩 指導教授：崔燕勇 博士

國立交通大學機械工程研究所

摘要

本研究模擬雙腔體串聯式的無閥式微幫浦。在模擬中，考慮兩腔體振膜為同相位 的振動與加入相位差後，流體的流動情形，並且考慮在不同背壓下的情況，並估算出淨 流量。模擬所得到的淨流量約為單腔幫浦的 1.2 倍。 此外我們建立一套幫浦理論(Lumped system)來分析流場，可應用於單一腔體和雙腔 的幫浦中，能夠估算出不同背壓、相位差下，進出口的流量。在理論中加入慣性項以及 考慮振動腔體的材料性質來比較模擬與理論的流量變化。並且估算出在不同背壓以及不 同相位差下的淨流量並與模擬的結果相驗證。 其 Lumped system 理論中，加入慣性項的模式其流量與模擬的結果相比，在單腔幫 浦中相當接近，而在雙腔中不同背壓以及不同相位的情況下皆相當的類似。其淨流量在 單腔幫浦與模擬的結果較為接近，在雙腔中則稍有差異。不過其理論在單腔、雙腔或是 多腔體幫浦的應用上，具有較好的可用性。

Analysis of the flow in double-chamber valveless

micropumps

Student：Shin-En Wu Advisor：Prof. Yeng-Yung Tsui

Institute of Mechanical Engineering

National Chiao Tung University

ABSTRACT

In this thesis, A double chamber valveless micropump is presented. We

consider two chambers with in-phase motion or different phase motion, and the

net flow rates in different back pressures are calculated. The double chamber

valveless micropump is 1.2 times the net flow rate of the single chamber.

In addition, the flow rate is also analyzed by using the lumped system. The

system can be used in the single chamber and double chamber. Then, the inertia

term and blocking pressure is added into the system. The results of the lumped

system are similar to the At last, the composition of the net flow rate in the

lumped system and CFD simulation is presented.

In the lumped system analysis, the model of inertia term added is similar to

the CFD simulation. In the single chamber, the net flow rate of the lumped

system and CFD simulation are approximate. But it has a little difference in the

double chamber. Generally, the lumped system is useful to some kinds of pump

models.

致謝

在交大機械所這兩年中，首先我要感謝崔燕勇教授在學業上的指導，使我了解並學 習到解決問題的方法以及態度，而在論文的研究上也給予我很多的指導，使我能夠完成 學業。接下來要感謝我的父母在背後的支持，在求學的過程中，他們總是給予我鼓勵和 關心並且默默的付出，使我更能夠安心的在課業上努力。 此外也感謝胡育昌學長以及呂學霖、鄭東庭學長在研究與課業上的幫助，和林仕 文、李俊岩同學在實驗室中互相的幫忙與鼓勵，並且也感謝學弟俊佑、祚昌、孝修平時 生活上的幫助。 最後再次感謝兩年中陪伴我一起努力的家人、同學以及朋友，有你們的支持我才能 夠一路的向前走，真的很謝謝大家。

符號表

符號： 定義：

0 r 壓電薄膜的半徑 1 r 振動腔體半徑 MAX d 壓電薄膜中心最大位移量 f 壓電薄膜的振動頻率 c P 單腔腔體壓力 in P 入口端壓力 out P 出口端壓力 K 壓力損失係數 Q α 慣性項的權重係數 P α 阻壓項的權重係數 1 P 雙腔中左側腔體的壓力 2 P 雙腔中右側腔體的壓力 F 通量 m_& 質量流率 u X 方向之速度項 v Y 方向之速度項 w Z 方向之速度項 Vur 速度向量 ρ _水的密度 μ _{水的黏滯係數} in Q 進口端流量 c Q 中間段Nozzle/Diffuser 流量

out Q 出口端流量 b P 阻壓 net Q 淨流量 l 長度 t A 的喉部面積 l A Nozzle/Diffuser 的擴張部面積

目錄 摘要...i ABSTRACT ...ii 目錄...vi 圖目錄...viii 一、 簡介...1 1.1 前言： ... 1 1.2 無閥式幫浦的工作原理： ... 3 1.3 文獻回顧： ... 4 1.4 研究內容： ... 6 二、數學模式...8 2.1 基本假設 ... 8 2.2 統御方程式 ... 8 2.3 壓電薄膜振動近似曲線 ... 8 2.4 邊界條件 ... 9 2.4.1 出入口邊界條件 ... 9 2.4.2 壓電薄膜邊界條件 ... 9 2.4.3 牆的邊界條件 ... 9 三、離散方法...10

3.1 有限體積法(Finite Volume Method) ... 10

3.2 離散化 ... 10 3.2.1 非穩態項(Unsteady term) ... 10 3.2.2 對流項(Convection term) ... 11 3.2.3 擴散項(Diffusion term) ... 11 3.2.4 源項(Source term)...12 3.3 線性代數方程式...13 四、PISO法求解步驟...14 4.1 計算面上質量流率...14 4.2 壓力修正方程式 ... 15 4.2.1 第一步壓力修正 ... 15 4.2.2 第二步壓力修正 ... 16 4.2.3 壓力修正方程式之整理 ... 18 4.3 邊界條件的給定 ... 18 4.3.1 出入口邊界條件的給定 ... 18 4.3.2 固定牆的邊界條件 ... 19 4.4 PISO演算過程... 20 五、無閥式微幫浦之理論分析...22 5.1 薄膜振動假設方程式 ... 22

5.2 微幫浦簡易理論分析(Lumped System Analysis)... 23 5.2.1 單腔理論分析 ... 24 5.2.2 雙腔理論分析 ... 28 六、結果與討論...31 6.1 三維模擬的結果分析 ... 31 6.1.1 單腔三維模擬結果 ... 31 6.1.2 雙腔三維模擬結果 ... 32 6.2 Lump Model理論分析... 36 6.2.1 單腔理論之結果分析 ... 36 6.2.2 雙腔理論於同相位差不同背壓之結果分析 ... 38 6.2.3 雙腔理論於不同相位差不同背壓之結果分析 ... 40 七、結論...43 參考文獻...45

圖目錄 圖1.1 微流體系統的應用【15】 ...47 圖1.2 壓電式微幫浦【1】 ...48 圖1.3 熱氣式微幫浦【1】 ...48 圖1.4 靜電式微幫浦【1】 ...48 圖1.5 蠕動型幫浦【10】 ...48 圖1.6 無閥式微幫浦的工作原理【15】 ...49 圖1.7 Stemme, 第一個無閥式的微幫浦【15】 ...49 圖1.8 Gerlach, 大張角的無閥式微幫浦【15】 ...49 圖1.9 Olsson, 單腔與雙腔並連無閥式微幫浦【15】 ...49 圖1.10 Yang,【14】方型無閥式微幫浦 ...50 圖1.11 Olsson所實驗的單腔無閥式微幫浦【15】 ...50 圖1.12 本文中所模擬的微幫浦尺寸 ...51 圖3.1 非結構性網格 ...52 圖3.2 Over-relaxed approach ...52 圖3.3 壓力邊界示意圖 ...53 圖4.1 固定牆的邊界剪應力示意圖 ...53 圖5.1 三種不同的近似曲線於壓電薄膜 ...54 圖6.1 本文中模擬單腔無閥式幫浦之網格 ...54 圖6.2 背壓 0Pa時，五週期中進出口流量對時間的變化(單腔) ...55 圖6.3 背壓 5900Pa時，五週期中進出口流量對時間的變化(單腔) ...55 圖6.4 背壓 0Pa時，在 1/2 腔體厚度截面上，1/4 週期時流場與壓力變化...56 圖6.5 背壓 0Pa時，在 1/2 腔體厚度截面上，3/4 週期時流場與壓力變化...57 圖6.6 背壓 5900Pa時，在 1/2 腔體厚度截面上，1/4 週期時流場與壓力變化...58 圖6.7 背壓 5900Pa時，在 1/2 腔體厚度截面上，3/4 週期時流場與壓力變化...59 圖6.8 本文中模擬雙腔無閥式幫浦之網格 ...60 圖6.9 雙腔幫浦不同網格中，淨流量隨週期的變化 ...61 圖6.10 雙腔幫浦不同網格下，其淨流量與背壓的關係 ...62 圖6.11a 背壓 0Pa時，五個週期下流量相對時間的分佈圖(雙腔)...63 圖6.11b 背壓 2950Pa時，五個週期下流量相對時間的分佈圖(雙腔)...63 圖6.11c 背壓 5900Pa時，五個週期下流量相對時間的分佈圖(雙腔)...64 圖6.12a 背壓 0Pa下，0.25T時的壓力分布...65 圖6.12b 背壓 0Pa下，0.75T時的壓力分布...66 圖6.13a 背壓 2950Pa下，0.25T時的壓力分布...67 圖6.13b 背壓 2950Pa下，0.75T時的壓力分布...68 圖6.14a 背壓 5900Pa下，0.25T時的壓力分布...69 圖6.14b 背壓 5900Pa下，0.75T時的壓力分布...70 圖6.15a 背壓 0Pa下，在 1/2 腔體厚度的截面上，0.25T時的流場分布...71

圖6.15b 背壓 0Pa下，在 1/2 腔體厚度的截面上，0.75T時的流場分布...71

圖6.16a 背壓 2950Pa下，在 1/2 腔體厚度的截面上，0.25T時的流場分布...72

圖6.16b 背壓 2950Pa下，在 1/2 腔體厚度的截面上，0.75T時的流場分布...72

圖6.17a 背壓 5900Pa下，在 1/2 腔體厚度的截面上，0.25T時的流場分布...73

圖6.17b 背壓 5900Pa下，在 1/2 腔體厚度的截面上，0.75T時的流場分布...73

圖6.18a 背壓 0Pa中，0.2T與 0.4T下三個Nozzle/Diffuser元件的流線圖...74

圖6.18b 背壓 0Pa中，0.6T與 0.8T下三個Nozzle/Diffuser元件的流線圖...75

圖6.19 單腔與雙腔的幫浦淨流量在不同背壓下的比較 ...76

圖6.20a背壓 0Pa，相位差 60°時，其進出口與中間段Nozzle/Diffuser元件的流量變化..77

圖6.20b背壓 0Pa，相位差 90°時，其進出口與中間段Nozzle/Diffuser元件的流量變化..77 圖6.20c背壓 0Pa，相位差 150°時，其進出口與中間段Nozzle/Diffuser元件的流量變化 78 圖6.21 相位差 60°時，背壓為 0Pa，在 1/4 週期處，幫浦中的壓力與流場變化 ...79 圖6.22 相位差 60°時，背壓為 0Pa，在 3/4 週期處，幫浦中的壓力與流場變化 ...80 圖6.23 相位差 90°時，背壓為 0Pa，在 1/4 週期處，幫浦中的壓力與流場變化 ...81 圖6.24 相位差 90°時，背壓為 0Pa，在 3/4 週期處，幫浦中的壓力與流場變化 ...82 圖6.25 背壓 2950Pa，相位差 60°時，其進出口與中間段Nozzle/Diffuser元件的流量變化 ...83 圖6.26 背壓 2950Pa，相位差 90°時，其進出口與中間段Nozzle/Diffuser元件的流量變化 ...83 圖6.27 相位差 60°時，背壓為 2950Pa，在 1/4 週期處，幫浦中的壓力與流場變化 ...84 圖6.28 相位差 60°時，背壓為 2950Pa，在 3/4 週期處，幫浦中的壓力與流場變化 ...85 圖6.29 相位差 90°時，背壓為 2950Pa，在 1/4 週期處，幫浦中的壓力與流場變化 ...86 圖6.30 相位差 90°時，背壓為 2950Pa，在 3/4 週期處，幫浦中的壓力與流場變化 ...87 圖6.31 背壓 0Pa和 2950Pa下，不同相位下淨流量的變化 ...88 圖6.32 比較背壓 0Pa，相位差 60°與 300°時，Q 之間的相位關係... 89 _c 圖6.33 比較背壓 0Pa，相位差 90°與 270°時，Q 之間的相位關係... 89 _c 圖6.34 比較背壓 0Pa，相位差 150°與 210°時，Q 之間的相位關係... 90 _c 圖6.35 比較背壓 2950Pa，相位差 60°與 300°時，Q 之間的相位關係... 91 _c 圖6.36 比較背壓 2950Pa，相位差 90°與 270°時，Q 之間的相位關係... 91 _c 圖6.37 比較背壓 2950Pa，相位差 150°與 210°時，Q 之間的相位關係... 92 _c 圖6.38 單腔微幫浦理論中，不同週期下淨流量的變化圖 ...93 圖6.39a不同背壓下，ModelA與ModelB(α_Q = 、1 0.5)進出口流量對時間的變化(單腔) 94 圖6.39b不同背壓下，ModelC(α_P = 、1 0.0001)進出口流量對時間的變化(單腔) ...95 圖6.40a背壓 0Pa下，出口流量對時間的變化圖(單腔)...96 圖6. 40b背壓 5900Pa下，出口流量對時間的變化圖(單腔) ...96 圖6.41 不同背壓下，ModelA與ModelB(α_Q= 、1 0.5)腔體壓力對時間的變化(單腔)...97

圖6.42 Model A與Model B中，不同背壓下的淨流量變化(單腔) ...98 圖6.43 Model A與Model C中，不同背壓下的淨流量變化(單腔) ...98 圖6.44 不同背壓下，ModelA與ModelB(α_Q= 、1 0.5)進出口流量對時間的變化(雙腔).99 圖6.45a背壓 0Pa下，流量對時間的變化圖(雙腔)...100 圖6. 45b背壓 5900Pa下，流量對時間的變化圖(雙腔) ...100 圖6.46a同相位下，背壓 0Pa時進出口的流量與壓電薄膜所掃過的流量的關係圖 ...101 圖6.46b同相位下，背壓 5900Pa時進出口的流量與壓電薄膜所掃過的流量的關係圖 .101 圖6.47 不同背壓下，ModelA與ModelB(α_Q= 、1 0.5)腔體壓力對時間的變化(雙腔)...102 圖6.48 Model A與Model B中，不同背壓下的淨流量變化(雙腔) ...103 圖6.49 Model A與Model C中，不同背壓下的淨流量變化(雙腔) ...103

圖6.50a 0 背壓時，Model A(Steady Model)中，不同相位角下一週期流量變化 ...104

圖6.50b 0 背壓時，Model A(Steady Model)中，不同相位角下一週期流量變化 ...105

圖6.51a 0 背壓時，Model B(Unsteady Model)中，不同相位角下一週期流量變化 ....106

圖6.51b 0 背壓時，Model B(Unsteady Model)中，不同相位角下一週期流量變化....107

圖6.52 在背壓 0Pa下，相位角 60°時模擬與理論之流量比較 ...108 圖6.53 在背壓 0Pa下，相位角 90°時模擬與理論之流量比較 ...108 圖6.54 在背壓 0Pa下，相位角 150°時模擬與理論之流量比較 ...109 圖6.55 背壓 0Pa，相位差 60°時，薄膜振動所掃過的流量與進出口和中間段Nozzle/Diffuser 元件的流量比較... 110 圖6.56 背壓 0Pa，相位差 90°時，薄膜振動所掃過的流量與進出口和中間段Nozzle/Diffuser 元件的流量比較... 110 圖6.57 背壓 0Pa，相位差 150°時，薄膜振動所掃過的流量與進出口和中間段 Nozzle/Diffuser元件的流量比較 ... 111 圖6.58 在背壓 0Pa下，模擬與理論在不同相位之淨流量比較 ... 112

圖6.59a 2950 背壓時，Model A(Steady Model)中，不同相位角下一週期流量變化 .. 113

圖6.59b 2950 背壓時，Model A(Steady Model)中，不同相位角下一週期流量變化 .. 114

圖6.60a 2950 背壓時，Model B(Unsteady Model)中，不同相位角下一週期流量變化 ... 115

圖6.60b 2950 背壓時，Model B(Unsteady Model)中，不同相位角下一週期流量變化 ... 116

圖6.61 在背壓 2950Pa下，相位角 60°時模擬與理論之流量比較 ... 117

圖6.62 在背壓 2950Pa下，相位角 90°時模擬與理論之流量比較 ... 117

圖6.63 在背壓 2950Pa下，相位角 150°時模擬與理論之流量比較 ... 118

一、 簡介

1.1 前言：

近 年 來 ， 因 為 半 導 體 製 程 蓬 勃 的 發 展 ， 微 機 電 技 術 (Micro electro-mechanical

technology)也隨之進步，借由這些微小元件所組成的系統即稱之為微機電系統(Micro

Electro Mechanical Systems, MEMS )，包含了微致動器、微型馬達、微感測器等元件。

其應用的範圍也很廣，在半導體、通訊、電子、生物醫學、航太產業等地方都有相當的 發展，因此隨著微機電系統被運用在不同的領域上，各種相關的研究也紛紛開始去做深 入的探討。 本文所要探討的是一個微流體系統，因為微機電技術的發展，也讓微流體系統元件 製作更加的進步。微流體系統包含有微幫浦、微流道、微型閥門、微流速感測器、微致 動器等。因為體積小、成本低、精確度高、反應時間短等優點，所以能夠被廣泛的應用 在不同的領域上，如圖1.1 所示。在生物晶片上如聚合脢連鎖反應晶片、蛋白質晶片、 基因晶片等，利用微機電技術在微流道中加入樣本再以微幫浦驅動，可以使檢測分析的 時間縮短、分析系統變的簡便；在微量皮下藥劑輸送系統上，經由控制信號與幫浦振動 的頻率，將藥物定時定量的輸入到人體內的循環系統中，如此一來便不需要經常的去注 射藥物，可以減少病患的痛苦；而在機械電子電機方面，可應用於微流體冷卻系統上， 利用閥門式的壓電型微幫浦將冷卻液傳輸到熱交換區來達到冷卻的效果，溫度較高的液 體再由微幫浦送到液態冷卻器中冷卻，達成一個循環。 其中微幫浦在整個微流體系統中扮演一個重要的角色，主要是用來驅動整個流體的 運作，由驅動方式可分為機械式與非機械式兩種模式，機械式微幫浦主要可分為位移式 微幫浦跟動力型微幫浦兩種，下表為各種較常見的機械式微幫浦，

機械式微幫浦 位移式微幫浦 動力型微幫浦 閥門幫浦(Check-valve Pump) 蠕動式幫浦(Peristaltic Pump) 無閥式幫浦(Valve-less Pump) 轉動式幫浦(Rotary Pump) 超音波幫浦(Ultrasonic Pump) 離心式幫浦(Centrifugal Pump) 表1.1 機械式微幫浦的分類【1】 其驅動的方式有壓電式(Piezoelectric)、氣壓式(Pneumatic)、靜電式(Electrostatic)、熱氣 壓式(Thermopneumatic)、電磁式(Electromagnetic)等。 依照驅動方式來分類的話，較常見的有： z 壓電式微幫浦(圖 1.2)： 這種幫浦是採用壓電材料所製成的，是利用材料的壓電效應，使薄膜震動。而所謂 的壓電效應，是在材料的兩端施加正負電場，會讓材料產生變形，這種效應是可逆的， 如果反過來讓材料產生變形，在材料兩側便會有電荷產生。 z 熱氣式微幫浦(圖 1.3)： 其外型是在主要工作流體的腔體對面另外加上一個腔體，裡頭含有第二種工作流 體，這種幫浦的工作原理便是加熱或冷卻第二種工作流體，使其體積膨脹與收縮，如此 反覆的運作便會壓迫主腔體，達到驅動流體的效果。 z 靜電式微幫浦(圖 1.4)： 靜電式微幫浦是在振動的薄膜上施加電壓，利用靜電的吸引力讓薄膜往上移動，此 時入口閥門打開使流體流入，之後取消施加的電壓使薄膜回復原狀，此時入口閥門關閉 出口閥門打開，於是流體便從出口流出。 依照腔體的種類來區分的話，又可以分為以下幾種類型： z 單腔式微幫浦： 其幫浦只擁有一個腔室，並加上不同的致動器來驅使腔體產生振動，其作用的模式 較為單純，大致可區分為排水與吸水兩種模式。

z 串聯式微幫浦： 是由兩個以上的腔體串連而成的幫浦，如蠕動式幫浦(Peristaltic Pump)，即是由三個 (圖 1.5)或三個以上的腔體串接而成的，驅動流體的方式是利用各腔體間不同相位差的振 動造成體積的變化來輸送流體。 z 並聯式微幫浦： 將兩個(或多個)腔體的進口端與進口端、出口端與出口端相連接，即為並聯式的幫 浦，腔體間的振動多為同一相位(in-phase)或是反相位(anti-phase)兩種模式來做討論。 非機械式微幫浦則是利用感應電場來驅動流體，沒有可動閥門等元件，但是每種非 機械式的微幫浦幾乎都有相對應的工作流體，較常見的有電動流體動力式幫浦(EHD)、 電磁式幫浦(MHD)、電滲式幫浦(Electro-osmotic Pump)、電泳式幫浦(Electro-phoretic

Pump)、電濕式幫浦(Electro-wetting Pump)、氣泡式幫浦(Bubble Pump)等。

機械式微幫浦所具有的優點為能夠驅動不同種類的流體，傳統的機械式微幫浦所採 用的止回閥系統，在薄膜振動頻率高時，閥門也必須相對應的配合，再加上壓力的變化， 容易造成閥門的磨耗及損壞。因此便發展出了無閥式的微幫浦，利用漸縮/漸擴管 (diffuser/nozzle)來取代以往原有的閥門。 在本研究中，則是模擬機械式微幫浦中，雙腔串聯的無閥式微幫浦。 1.2 無閥式幫浦的工作原理： 無閥式微幫浦的工作原理，主要可分為兩種模式來討論，一個為吸水模式，一個則 為排水模式，如圖1.6。 z 吸水模式： 此時薄膜向上移動，腔體體積增大，流體經由兩端的出入口流進腔體由於從 diffuser 方向流進的壓力損失較小，從nozzle 方向流進的壓力損失較大，所以進口端(diffuser 方 向)流進的流量較出口端(nozzle 方向)來的大。 z 排水模式： 此時薄膜向下移動，腔體的體積減少，所以流體從兩端向外排出，這時入口端的流

流量較小，出口端的壓力損失小，流量較大。

1.3 文獻回顧：

在文獻回顧中，我們探討前人相關的文獻中，對微幫浦的實驗、模擬、以及理論和

Nozzle/Diffuser 元件相關的研究。

以下為微幫浦相關的實驗文獻：

1990 年，SHUICHI SHOJI et al.【2】, 利用矽板製作出有閥門的壓電式微幫浦，並將其串 聯與並聯，測得最大流量約為 40μl/min，最大幫浦壓力約為 1mH2O。

1993 年，Stemme et al.【3】 作出第一個無閥式的微幫浦(圖 1.7)，代替以往有閥門的微

幫浦，製造過程採用的材料為黃銅，其幫浦的直徑為19mm，實驗的結果中，最大的流

量為16ml/min，最大壓力為 2mH2O，頻率為 100Hz。

1995 年，Gerlach et al.【4】 本文中採用大角度的漸縮管(圖 1.8)，並且定義 nozzle 方向為

正向流動方向，流場在雷諾數小於 15 時為層流，大於 100 則變為紊流，其幫浦壓力最 大為 7kPa，振動頻率最高為 10kHz。 1995 年，Ollson et al.【5】 第一個用黃銅製造出圓形平板狀的無閥式微幫浦，將兩個腔 體並聯，每個腔體上下都有一個壓電薄膜，以期獲得較大的流量。每個腔體的厚度為 1mm，直徑為 13mm，並發現在相位差 180 度、頻率 540Hz 的時候有最大的流量為 16ml/min。 1997 年，Olsson et al.【6】 利用熱壓法製造平面式圓形腔體的無閥式微幫浦，腔體採用 單一腔體與雙腔並聯兩種來作比較(圖 1.9)，其中閥門(diffuser)張角為 7 度，實驗的結果 中，壓電薄膜在頻率 2200Hz 時的最大流量約為 1ml/min，最大的背壓達 5900kpa，而在 並聯的情況下，最大流量約為 1.9ml/min。 2001年，Nguyen et al. 【7】利用印刷電路板做為基材，製造出單腔的無閥式微幫浦以 及三個腔體的蠕動式幫浦，其最大的流量可達3ml/min，並利用CFD計算軟體CFD-ACE 模擬幫浦內部的流動。 2004年，Lee et al. 【8】利用矽板以及玻璃製作出蠕動式的微幫浦，並採用壓電材料來 驅動，其工作流體為水和空氣，在工作流體為水的情況下流量可達12μl/min，而工作流

體為空氣的情況下，流量可達60μl/min。對於供應電壓與壓電材料振幅之間的關係，採 用Coventor simulator與CoventorWare兩套軟體模擬，其結果和實驗有相當的準確性。 2006 年，關恕【9】對蠕動型的微幫浦的性能去做分析，對不同頻率以及不同的薄膜作 動的規律來探討其流量的輸出，並且改良電路的設計以提高傳輸的流量，可提升淨流量 1.9~2.8 倍。 2008 年，Jang et al. 【10】運用壓電材料作為致動器，來驅動蠕動式微幫浦，其輸入訊 號為步階函數，並且實驗三組不同的相位序列來比較其效能，最後發現，以4 階段相位 序列的訊號輸入並在頻率為100Hz 時有最佳的效能，流量為 17.6μl/min。 以下則是有關微幫浦的模擬：

1997 年，Athavale et al.【11】用來模擬的軟體為 CFDRC，採用 CFD-ACE 中的

Fluid-Thermo-Structure 模組來處理液固介面間的關係，而模擬的微幫浦為 1993 年 Stemme 所做的第一個無閥式微幫浦，其最大流量與實驗非常符合，而在最大背壓的地 方誤差較大。 2004 年，羅卓錚【12】利用商用軟體 CFDRC 模擬無閥式的微幫浦，研究中採用新的設 計，在前後進出口的流道中置入擋體，利用擋體的張角便可以產生類似於漸擴/漸縮管的 效果，結果發現最佳擋體的張角為5 度，其模擬的結果與實驗比較後，發現在最大背壓 下，最大淨流量有20%的誤差。 2005 年，曾裕博【13】利用 Fluent 來模擬方形的微幫浦，利用給定一個入口的振動曲 線來取代移動網格的建立，並且分析幫浦的效率，推導出的式子與漸擴管的效率有關， 並且比較不同張角的漸擴管下，效率的差異。

2006 年，Yang et al.【14】利用 CFD-ACE 模擬單腔、雙腔並聯、雙腔串聯(圖 1.10)，將

其振動週期分為四個區段。在雙腔並聯、無相位差的情況下，其流量約為單腔微幫浦的

兩倍，而在雙腔串聯的情況下，微幫浦的流量與兩個腔體振動的相位差有關係，並發現

在相位差為90 度時，有最佳的流量。

做比較，另外還建立了一套簡單的數學理論，借由找出噴嘴/擴散器的壓力係數與流量的 關係來分析不同背壓下的流量、腔體壓力以及效率。 下列為有關微幫浦理論的研究： 1998 年，Olsson et al.【16】提出利用塊狀模式來分析壓電式的無閥式微幫浦，利用 MATLAB 計算連續方程式與伯努力方程式所整合出的三個方程式組，其結果與實驗結 果非常相似，利用塊狀模式可以得出幫浦內的壓力與流量值。 1998 年，Ullmann【17】，本文討論微幫浦的數學模式，在單一腔體的微幫浦中，一週期 裡面採用三種模式來討論其流量、壓力的變化，也針對雙腔串聯的微幫浦作其簡單的分 析，並提出六種不同的連接方式。 2001 年，Pan et al.【18】分析薄膜震動的方程式與微幫浦內壓力的關係，並且從三組不 同的壓力關係，去求得微幫浦流量的解析解，並找出最佳的工作範圍。 2003 年，Pan et al.【19】將壓電材料的振動方程式與流體的流動結合，考慮其流體的慣 性並討論進出口的背壓與壓電薄膜的靜電力來求得其流量的解析解。 下列為Nozzle/Diffuser 元件相關的研究： 1996 年，Olsson et al.【20】在設計上將漸擴管分成三個區段來分析，並且找出各區段中 的壓力損失係數，將其合併後可以求出漸縮/漸擴管的工作效率，在文中並提出利用管壁 摩擦係數求出第二區段的壓力差和流量的關係與實驗的結果比較，實驗採用的漸擴管長 為1.45μm 至 3.95μm，平均喉部面積為 30μm×100μm。 2004 年，Yang et al.【21】本文主要在分析漸縮/漸擴管中，壓力損失係數與雷諾數間的 關係，並且和實驗的結果互相驗證，可發現壓力損失係數會隨雷諾數而遞減，而且漸縮 管的壓力損失係數大於漸擴管。由於漸擴管在張角 20 度時流場會分離，所以實驗與理 論的值會有一段差距。 1.4 研究內容： 在研究中，其模擬的外型為雙腔串聯的無閥式微幫浦，所參考的幾何外型以 1997 年Ollson et al.【6】 的實驗尺寸(圖 1.11)為主，並且以擴張管串接兩個幫浦來分析，其 腔體的半徑為 3mm，深度為 0.2mm。腔體上壓電薄膜的半徑為 2mm，在兩腔體連接處

與進出口共以三個Nozzle/Diffuser 元件來連接，其 Nozzle/Diffuser 元件的張角均為 7 度， 長 1mm，喉部的寬度為 0.1mm，如圖1.12 所示。 在模擬上，其壓電薄膜的振動採用實驗中的數據，為 2200Hz，而薄膜的振動情形 我們採用文獻【13】中，梯型的振動曲線來近似其振動情形，以簡化邊界條件的設定。 模擬中，考慮兩腔體振膜為同相位的振動與加入相位差後，流體的流動情形，並且考慮 在不同背壓下的情況，估算出淨流量與背壓的關係。 在理論的分析上，因為在雙腔體幫浦中，區分其工作模式較為複雜，故將文獻【17】 中，Ullmann 所提出的數學模式將其改寫，使其不僅適用於單腔的幫浦，也能夠應用於 雙腔幫浦中，並且在式中加入慣性項來分析，以符合物理現象，此外，也考慮振膜的材 料性質加入阻壓來測試。 最後比較單腔和雙腔模擬和理論中，淨流量與背壓的關係，與不同相位差下其淨流 量的關係，並且驗證其準確性。

二、數學模式

在本章節中將我們採用的數學模式，分成基本假設、統御方程式、壓電薄膜振動近 似曲線、邊界條件四個部份來做介紹。 2.1 基本假設 在本研究中，我們採用水做為微幫浦的工作流體，另外我們還做了以下幾種假設： z 非穩態流場：因為壓電薄膜的振動隨時間而改變，所以進出口的流量也隨時間的變 化而不同。 z 不可壓縮流：假設水的密度為定值。 z 忽略重力項 2.2 統御方程式 求解流場的方程式為連續方程式與動量方程式 連續方程式： ∇ ⋅

( )

ρvr =0 (2.1) 動量方程式：

( ) ( )

v vv P

( )

v t ρ ρ μ ∂ + ∇ ⋅ = −∇ + ∇ ∇ ∂ r rr r (2.2) 上面二式中， ρ：水的密度，μ：水的黏滯係數。 2.3 壓電薄膜振動近似曲線 對於壓電薄膜的振動曲線，是一個複雜的高階偏微分方程式，為了簡化求解的困 難，因此本文中採用梯形曲線來近似。 z 梯形曲線： 1 0 1 0 0 1 cos(2 ) , ( , ) cos(2 ) , MAX MAX d ft if r r S r t r r d ft if r r r r r π π − ⋅ < ⎧ ⎪ =_⎨− − ⋅ ⋅ < < ⎪ ₋ ⎩ (2.3)

上式中，d_MAX為振動薄膜從水平位置到上死點的最大位移量，在本文中d_MAX為1.0 mμ ， 0 r 為腔體的半徑，r 為壓電薄膜的半徑， f 為壓電薄膜的振動頻率，參照實驗數據為₁ 2200Hz。 2.4 邊界條件 2.4.1 出入口邊界條件 在微幫浦左右兩側的出入口皆為壓力的邊界條件，在模擬上皆將其假設為出口來做 計算，可以考慮不同壓力下的情況，所採用的是對流邊界條件(Convective Boundary Condition)。 2.4.2 壓電薄膜邊界條件 壓電薄膜的邊界因為薄膜的振動，必須採用移動邊界，網格因此也必須改成移動網 格的形式，所以在計算的過程會耗費較久的時間。而在文獻【13】中，曾等人發現，當 振幅低於 1000nm 且振動頻率低於 3000Hz 之下，可以將移動邊界的條件用一個入口速 度來取代，不僅能夠加快運算時間，而且誤差很小。因此在本研究中，壓電薄膜的振動 是採用一個入口的速度分佈曲線來去近似它，避免因為移動網格所增加的龐大的運算 量。 在本研究中壓電薄膜這部份的邊界，所給定的薄膜振動速度即為薄膜振動曲線對時 間微分後的速度分佈： z 梯形曲線的振動速度：

(

)

(

)

1 0 1 0 0 1 2 sin 2 , ( , ) ( , ) 2 sin 2 , MAX MAX fd ft if r r S r t V r t r r fd ft if r r r t r r π π π π ⋅ < ⎧ ∂ ⎪ = _∂ =⎨ _⋅ − _{⋅ < <} ⎪ ₋ ⎩ (2.4) 2.4.3 牆的邊界條件 牆上採用無滑移的邊界條件(No-slip Condition)，Vw=0 uur 。所以必須考慮剪應力在靠 近牆上的相鄰網格中所受的影響。

三、離散方法

現今在數值的求解上發展出了許多的方法，像是有限元素法、有限差分法、有限體 積法等，本文中求解統御方程式的方法，是採用有限體積法，顧名思義，是將整個計算 的模型分成一個個微小的體積來處理，相較於其他的數值方法，有限體積法在工程運用 上較多，因為求解的環境多為流動的性質，尤其是與通量有關的物理現象，對於較複雜 的模型採用有限體積法，能夠較容易去求解，且較具有物理意義，此種方法也能夠應用 在非結構性網格上(圖 3.1 )。

3.1 有限體積法(Finite Volume Method)

由第二章中的動量方程式(2.2 式)，我們可以得到下列通式：

( )

( )

v

(

)

q t ρφ ρ φ μ φ ∂ + ∇ ⋅ = ∇ ⋅ ∇ + ∂ r (3.1) 將上式各項作體積分，表示成積分形式：

( )

(

)

V V V V dV v dV dV qdV t ρφ ρ φ μ φ ∂ _{+ ∇ ⋅ = ∇ ∇ +} ∂

∫∫∫

∫∫∫

∫∫∫

∫∫∫

r (3.2)

右側第一項為非穩態項(Unsteady term)，第二項為對流項(Convection term)；左側第

一項為擴散項(Diffusion term)，第二項為源項(Source term)。其中源項的來源是壓力項，

所以 q= −∇ 。 P

利用高斯散度定理(Gauss divergence theorem)，我們將(3.2)式積分形式的動量方程式

轉換成面積分的形式： V S S V dV v d S d S qdV t ρφ ρ φ μ φ ∂ + ⋅ = ∇ ⋅ + ∂

∫∫∫

∫∫

∫∫

∫∫∫

r ur ur (3.3) 3.2 離散化 本節中，將詳述非穩態項、對流項、擴散項以及源項的離散方法。 3.2.1 非穩態項(Unsteady term) 在非穩態項中，當體積隨時間的變化而改變時，我們將此項離散成如下列所示： ( o) V V dV t t ρ ρφ φ φ ∂ ⋅ Δ = − ∂

∫∫∫

Δ (3.4) 上式中，變數φ是在這時間下的值，而φo 則是在上一個時間下的值，故等式右側第

二項為已知值，在離散後我們將其移入源項中。 3.2.2 對流項(Convection term) (3.3)式的對流項，是由高斯散度定理，將體積分轉換成面積分，我們再將以表示成 差分的形式： ( ) C f f f f f f f f S v d S v S m F ρ φ⋅ =

∑

ρ φ ⋅ =

∑

φ =

∑

∫∫

r ur r uur & (3.5) 其中，Sf uur ：面的法向量

m& ：面上的質量流率(mass flow rate) _f

下標 f ：在控制體積任一面上中點的值 C f F ：對流通量(flux) 對於面上的值，我們利用相鄰網格線性內插求得。 在求解對流項的過程中，我們採用一階上風法(UD)和中央差分法(CD)兩種的混合 型。因為單獨採用一階上風法的話，會造成很大的數值擴散；而採用中央差分法的話， 會有震盪的的產生以致於較不易收斂。混合式表示如下： C ( UD) ( CD UD) f f f f F = F +γ F −F (3.6) 其中 UD f F 為上風法產生的對流通量、 CD f F 為中央差分法產生的對流通量，γ 為一個可調 整的值，介於0 和 1 之間，其值為 0 時，完全採用上風法；其值為 1 時，即為採用中央 差分法，我們為其準確性，通常取較接近1 的值。本研究中採用的γ 值為0.9。 我們將(3.6)式整理成下列形式：

[

]

[max( ,0) min( ,0) ] (1 ) max( ,0) min( ,0) C f f P f C f P C f P f C F m m m w w m m φ φ γ φ φ φ φ = + ⎡ ⎤ + _⎣ − + − − _⎦ & &

& & & (3.7)

其中w為加權因子(weighting factor)，值介於 0~1 之間。

上式中上風法求得的第一項至於矩陣中，而將第二項中央差分法和上風法的差值移

將經由高斯散度定理，從體積分轉換成面積分的擴散項離散之後，得到下式： D ( ) D f f f f f f S F =

∫∫

μ φ∇ ⋅d Sur=

∑

μ ∇φ ⋅Suur=

∑

F (3.8) 其中，μ_f ：面上的黏滯係數， D f F ：擴散通量(flux) 對於擴散項的處理，我們採用 over-relaxed 法，首先我們令面法向量表示如下： Sf = +d (Sf −d) uur ur uur ur (3.9) dur定義為主格點P 至相鄰格點 C 的向量(圖 3.2 )，其大小將決定計算中擴散量的大小，dur 的定義表示如下： 2 f f d d s f S S d e e e δ S δ = = ⋅ ⋅ uur uur ur uur ur uur ur ur uur (3.10) 其中，δur：主格點P 至相鄰格點 C 的距離向量 ed uur ：主格點P 至相鄰格點 C 的單位向量 eur_s：沿著Suur_f 方向的單位向量 把(3.9)式代入(3.8)中，將其整理如下： 2 ( ) ( ) ( ) f f D f C P f f f f S F S d S μ φ φ μ φ δ = − + ∇ ⋅ − ⋅ uur uur ur ur uur (3.11) 在計算上，我們將等式右邊第一項置於係數矩陣中，而將第二項以前一個時間下的 值移入源項中。 3.2.4 源項(Source term) 利用高斯散度定理，我們將源項中的壓力梯度表示如下： 1 1 1 f f f V S P PdV Pd S P S V V V ∇ = ∇ = = Δ

∫∫∫

Δ

∫∫

Δ

∑

ur uur (3.12) 上式中，面上的壓力不包含邊界上的壓力，為了要求得邊界上的壓力梯度，邊界上的壓 力必須要另外去計算。 其邊界壓力可以由下式求得(圖 3.3 )：

P_b−P_P = ∇ ⋅P δur (3.13) 其中，下標b表示壁面上的點，δur則是主格點 P 到b的距離向量。 因此從(3.12)式，壓力梯度可以近似如下： 1 f f 1 b b f f f f b P P S P S P S V V ≠ ⎛ ⎞ ∇ = = _⎜ + _⎟ Δ

∑

_{Δ ⎝}

∑

_⎠

uur uur uur

(3.14) 我們將(3.14)代回(3.13)式，便可以得知邊界上的壓力P 。 _b 1 1 1 P f f f b b b P P S V P S V δ δ ≠ ⎛ ⎞ + ⋅ ⎜ _Δ ⎟ ⎝ ⎠ = ⎛ _{− ⋅} ⎞ ⎜ _Δ ⎟ ⎝ ⎠

∑

uur ur uur ur (3.15) 3.3 線性代數方程式 將前小節所導出的非穩態項、對流項、擴散項、源項合併後，可以歸納出一個線性 代數方程式： A_Pφ_P =

∑

A_{C C}φ +Q (3.16) 其中，A_P A_C V t ρφ⋅ Δ = + Δ

∑

(3.17) 2 max( ,0) f f C f f S A m S μ δ = + − ⋅ uur ur uur & (3.18)

[

(1 )

]

max( ,0) min( ,0) ( ) f P C f P f C o f f f f Q m w w m m V S d P V t γ φ φ φ φ ρφ μ φ ⎡ ⎤ = _⎣ − + − − _⎦ ⋅ Δ − ∇ ⋅ − + ∇ ⋅ Δ + Δ

∑

& & &

uur ur (3.19)

在(3.18)式中，A 的值是由擴散項與對流項分出來的；而在(3.19)式中，源項的值則是由C

四、PISO 法求解步驟

在處理速度跟壓力的偶合關係式上，現今發展出了很多種不同的方法，像是 SIMPLE、SIMPLEC、PISO 等。主要是因為我們只對動量方程式離散，但是在連續方程 式上卻不滿足，因此必須要透過修正壓力與速度來做修正。本研究中，所採用的方法為 PISO 演算法，此方法對於暫態的求解，不需要再去反覆疊代，能夠減少求解的時間。 在接下來的幾節，我們將詳細敘述其過程。 4.1 計算面上質量流率 在上一節的計算中，我們利用動量方程式求解出速度場，但是並未滿足連續方程 式，因此必須找出壓力與速度的關係式，並求解出面上的速度來計算面上的質量流率。 我們將(3.16)式中的壓力項自源項中提出，表示成下列所示： * * * p p P P _P V V H P A ⎛Δ ⎞ = −_{⎜ ⎟} ∇ ⎝ ⎠ uur uuur (4.1) 其中 * _' * _{C C} p P A V Q H A + =

∑

uur uuur ；Q 為不包含壓力項的源項。 (4.2) ' 而面上的速度我們也可以從(4.1)式來去近似，表示如下： * * * f f f P f V V H P A ⎛Δ ⎞ = −_{⎜ ⎟} ∇ ⎝ ⎠ uur uuur (4.3) 其中 _f* _f f P f V H V P A ⎛Δ ⎞ = +_{⎜ ⎟} ∇ ⎝ ⎠ uuur uur (4.4) 將其代入(4.3)式可得： * _{( )} * f f f f P f P f V V V V P P A A ⎛Δ ⎞ ⎛Δ ⎞ = +_{⎜ ⎟} ∇ −_{⎜ ⎟} ∇ ⎝ ⎠ ⎝ ⎠ uur uur (4.5) 上標 “ “表示是由主格點與相鄰格點內插所得，如下所示： ∇Pf = ∇ + − ∇w P_C (1 w) P_P (4.6)

Vf =wVC + −(1 w V) P

uur uur uur

(4.7) 而(4.5)式中 P f V A ⎛Δ ⎞ ⎜ ⎟ ⎝ ⎠ 的計算，可近似為下式： 1 2 P f P P P C V V V A A A ⎡ ⎤ ⎛Δ ⎞ ⎛Δ ⎞ ⎛Δ ⎞ = ⎢ + ⎥ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ _⎣⎝ ⎠ ⎝ ⎠ _⎦ (4.8) 所以質量流率可改寫為： * * * * ( ) ( ) f f f f f f f f f f f P _f f f f f f f P f V m V S V S P P S A V V S P P d A ρ ρ ρ ρ ρ ⎛Δ ⎞ = ⋅ = ⋅ − _{⎜ ⎟} ∇ − ∇ ⋅ ⎝ ⎠ ⎛Δ ⎞ ≈ ⋅ − _{⎜ ⎟} ∇ − ∇ ⋅ ⎝ ⎠

uur uur uur uur uur

& uur uur ur 2 f f f f ( C P) f P f V S V S P P P A S ρ ρ δ δ ⎛_Δ ⎞ ⎜ _{⎟ ⎡ ⎤} = ⋅ − _{⎜ ⎟ ⎣} − − ∇ ⋅ _⎦ ⋅ ⎜ ⎟ ⎝ ⎠ ur uur uur ur ur ur (4.9) 4.2 壓力修正方程式 前一節中以求得的壓力 * P 和速度Vur*仍沒有滿足連續方程式，所以接下來我們必須 要加上修正量去修正壓力與速度。 4.2.1 第一步壓力修正 我們定義修正後的速度為Vur**；修正後的壓力為P ，所以在第一次修正中，速度與** 壓力的關係式表示如下： ** * ** p p P P _P V V H P A ⎛Δ ⎞ = −_{⎜ ⎟} ∇ ⎝ ⎠ uur uuur (4.10) 其中Vur**=Vur*+Vur'；P**=P*+ ，上標 “**” 表示修正後的值； “*” 表示修正前的值； P' “'” 表示修正的值。 接下來我們將(4.10)式與(4.1)式相減，其速度修正量和壓力修正量表示成下列的關 係：

' ' P P P P V V P A ⎛Δ ⎞ = −_{⎜ ⎟} ∇ ⎝ ⎠ uur (4.11) 而主控面上的速度修正式也可以表示成： ' ** * ' f f f f P f V V V V P A ⎛Δ ⎞ = − = −_{⎜ ⎟} ∇ ⎝ ⎠

uur uur uur

(4.12) 因此由上式可得修正後的質量流率： ' ** * * ' * ' ' ( ) f f f f f f f f f P f f f f f f f P f P f V m m V S m P S A V V m P d P S d A A ρ ρ ρ ρ ⎛Δ ⎞ = + ⋅ = − _{⎜ ⎟} ∇ ⋅ ⎝ ⎠ ⎛Δ ⎞ ⎛Δ ⎞ = − _{⎜ ⎟} ∇ ⋅ − _{⎜ ⎟} ∇ ⋅ − ⎝ ⎠ ⎝ ⎠

uur uur uur

& & &

ur uur ur & 2 * ' ' ' _{f f} f ( _{C P}) _{f f} ( _f ) P f f P f S V V m P P P S d A S A ρ ρ δ ⎛Δ ⎞ ⎛Δ ⎞ = − _{⎜ ⎟} − − _{⎜ ⎟} ∇ ⋅ − ⋅ ⎝ ⎠ ⎝ ⎠ uur uur ur & ur uur (4.13) 為了滿足質量守恆連續方程式，對所有面上流量的總合為 0，其式子表示如下： ** ₀ f f m =

∑

& (4.14) 將(4.14)式代入(4.13)式中，可以得到壓力修正方程式： ' ' 1 2 P P C C P P f A P =

∑

A P +S +S (4.14) 其中， 2 f C f P f f S V A A S ρ δ ⎛Δ ⎞ = _{⎜ ⎟} ⋅ ⎝ ⎠ uur ur uur P C f A =

∑

A 而後兩項源項表示如下： * 1 P f f S = −

∑

m& (4.15) ' 2 ( ) P f f f f P f V S P S d A ρ ⎛Δ ⎞ = _{⎜ ⎟} ∇ ⋅ − ⎝ ⎠

∑

uur ur (4.16) 4.2.2 第二步壓力修正 在第二次壓力修正中，壓力與速度的關係式如同第一次壓力修正式的形式一樣，其

關係式表示如下： *** ** *** ' p p P P _P V V H P Q A ⎛Δ ⎞ = −_{⎜ ⎟} ∇ + ⎝ ⎠ uur uuur (4.17) 接下來，我們將(4.17)式與(4.10)式相減後，得到速度修正與壓力修正量的關係式，表示 如下： '' *** ** '' '' p p p p p P p V V V V H P A ⎛Δ ⎞ = − = −_{⎜ ⎟} ∇ ⎝ ⎠ uuur uur uur uur

(4.18)

然後，面上的質量流率我們將其表示為： *** ** ''

f f f

m& =m& +m& (4.19)

從(4.18)式，我們可以推得第二次修正後的質量流率並將其表示如下： '' *** ** '' C C f f f f f f f f P f P A V V m m P S S A A ρ ⎛Δ ⎞ ρ = − _{⎜ ⎟} ∇ ⋅ + ⋅ ⎝ ⎠

∑

uur uur uur & & '' ** '' '' ( ) C C f f f f f f f f f P f P f P A V V V m P d P S d S A A A ρ ⎛Δ ⎞ ρ ⎛Δ ⎞ ρ = − _{⎜ ⎟} ∇ ⋅ − _{⎜ ⎟} ∇ ⋅ − + ⋅ ⎝ ⎠ ⎝ ⎠

∑

uur ur uur ur uur & (4.20) 因為要滿足質量守恆連續方程式，所以所有面上的質量流率總和為 0，表示如下： *** ** '' '' ₀ f f f f f f f f m = m + m = m =

∑

&

∑

&

∑

&

∑

& (4.21) 我們將(4.21)式代入(4.20)式，最後將其整理成類似(4.14)式的線性代數方程式： '' '' 1 2 P P C C P P f A P =

∑

A P +S +S (4.22) 其中， _{P C} f A =

∑

A 2 f C f P f f S V A A S ρ δ ⎛Δ ⎞ = _{⎜ ⎟} ⋅ ⎝ ⎠ uur ur uur 後兩項源項表示如下： '' 1 C C f P f f P A V S S A ρ = −

∑

⋅ uur uur (4.23) '' 2 ( ) P f f f f P V S P S d A ρ ⎛Δ ⎞ = _{⎜ ⎟} ∇ ⋅ − ⎝ ⎠

∑

uur ur (4.24)

4.2.3 壓力修正方程式之整理 在兩次的壓力修正式中，SP₂為非正交的源項，因此如果是正交的網格的話，SP2的 值為0。 而在求解的過程中，我們對每次的壓力修正式中，又將它再區分成兩部份來計算， 以得到較準確的值。第一部份中我們忽略S ，源項中只包含_P2 S ；而在第二部份中，源_P1 項只包含S_P2，將其整理如下。 z 正交修正(Orthogonal Corrector)部份： '(1) '(1) 1 P P C C P f A P =

∑

A P +S z 非正交修正(Non-orthogonal Corrector)部份： '(2) '(2) '(1) 2( ) P P C C P f A P =

∑

A P +S P 於正交修正式中得到壓力修正 '(1) P ，以此計算S_P2中之梯度，再代入壓力修正求 '(2) P ，此非正交修正步驟可重複執行，但通常只執行一次便足夠。 4.3 邊界條件的給定 4.3.1 出入口邊界條件的給定 因為在本研究中，出入口的邊界條件為一壓力邊界，因此我們必須在出入口處設定 壓力值，再藉由修正出入口壓力來修正出入口的質量流率。 其步驟表示如下： 第一步：將壓力邊界相鄰的網格點上的壓力設定成給定的值。 第二步：在求解壓力修正式時，我們必須設法讓壓力修正量的值為 0，也就 是要滿足P'= (圖 3.3 )。 0 所以，壓力修正方程式的係數與源項我們做以下的調整： ' ' 1 2 P P C C P P f A P =

∑

A P +S +S 其中S_P1=S_P2 = ，0 A_C = 0 第三步：求解完壓力修正方程式後，邊界相鄰的網格需再讓其滿足質量守恆

定理，即 f 0 f

m =

∑

& 。

邊界上的質量流率便可表示成下式：

m&f₄ = −(m&f₁+m&f₂+m&f₃)

在計算壓力邊界的速度上，我們採用對流邊界條件(Convective Boundary Condition)，其數學表示式為： u_c 0 t z φ φ ∂ ₊ ∂ ₌ ∂ ∂ (4.25) 對上式作離散後表示如下： 0 n o n n b b b C c u t z φ −φ φ −φ + = Δ Δ (4.26) 其中上標 “n”、“o” 代表新的值與舊的值，下標 “b”、“c” 則代表邊界和其相鄰的網格， 將(4.26)式整理如下： 1 o n n b C b Cr Cr φ φ φ = + + (4.27) 其中_Cr uc t z Δ = Δ ，Cr 為 Courant number。 因此，邊界上的外插速度我們可由(4.27)式求出： u ib +v jb +w kb r r r 最後，壓力邊界上的速度可修正為： 4 ( ) ( ) f b b b b b x b y b z m V u i v j w k u S v S w S ρ = + + + + uur & r r r (4.28) 4.3.2 固定牆的邊界條件 在牆上，我們採用的是無滑移的邊界條件(No-slip Condition)。在壁面與相鄰格點 上，我們以 P 代表相鄰壁面的主格點、Sur為壁面上之面向量、格點上的速度為Vuur_P、壁面 的速度為Vuur_b、壁面朝格點內的法向量為nr，如圖4.1 所示。 其中法向量nr表示為：

S ix S jy S kz n S + + = − r r r r ur (4.29) 而相對速度Vur為相鄰壁面主格點上的速度向量與壁面速度的差值，表示為： Vur uur uur=V_P −V_b (4.30) 所以，垂直壁面的速度向量可由下式求得： ( ) Vuur_⊥ = V n nur r r⋅ (4.31) 平行壁面的速度向量可表示為： Vuur ur uur_// = −V V_⊥ (4.32) 故壁面剪應力為： _b f S V_// n μ τ δ = ur uur uur (4.33) 其中，δ ： P 點到壁面之垂直距離。 n 4.4 PISO 演算過程 以下為 PISO 整個演算的過程： 1. 給定邊界的初始值。 2. 修正進口邊界的速度和流量。 3. 求解動量方程式得到預測的速度場Vur*。 4. 求解第一次壓力修正值 ' P ，在進一步求得第一次修正過後的速度Vur**、壓力 ** P 、流量m&**_f 5. 求解第二次壓力修正值 '' P ，在進一步求得第二次修正過後的速度Vur***、壓力 *** P 、流量m&***_f 。 6. 更新壓力邊界上的流量。 7. 重複上述步驟 2 至 6 即可求得每個暫態下的值，直到所限定的時間，最後輸出 結果。

在求解代數方程式的係數矩陣上，為了減少求解矩陣所需的記憶體，我們採用疊代

的方式，避免用直接求解所需負荷的記憶體空間。在求解動量方程式上，因為其係數矩

陣不具對稱性質，所以採用BICG 法，壓力修正方程式中，其係數矩陣具有對稱性，故

五、無閥式微幫浦之理論分析

5.1 薄膜振動假設方程式 為了模擬壓電薄膜的振動情形，我們必須求得薄膜振動的模式，以便於模擬邊界的 給定以及理論的計算。首先 Pan【18】 提出了一個解析解，將流體力學和結構力學做 整合(Fluid-membrane Coupling)，結合了一些材料的係數與電場相關的係數，其數學式 如下： 2 4 2 e W D W h f P t ρ ∂ ∇ + = − ∂ (5.1) 其中 3 2 12(1 ) Eh D v = − ：為腔體材料係數 E ：基材的楊式係數。 v：基材的浦松比。 h：基材的厚度。 ρ：基材的密度。 f ：週期性的作用力，和壓電係數與電場有關。 e P ：壓力。 因為(5.1)式為二階且複雜的偏微分方程式，所以求解較為困難，因此我們試圖找尋 其它近似的方程式來簡化。因此，有研究將其簡化成流場的問題，不去考慮材料的特性， 單純用曲線去模擬壓電材料的振動。 在本研究中，壓電薄膜的振動則是採用下式中的梯形振動曲線來近似。以下列出幾 種不同形式的振動曲線(圖 5.1)。 z 四次曲線： Timoshenko 【22】 將圓形平板的振動假設成一個四次曲線，數學式可表示成 下列式子： 2 2 max 0 ( ) 1 r S r d r ⎛ _{⎛ ⎞} ⎞ ⎜ ⎟ = _{⋅ − ⎜ ⎟} ⎜ _{⎝ ⎠} ⎟ ⎝ ⎠ (5.2) z 二次曲線(Parabolic)：

在文獻 【22】 中，假設曲線為隨半徑變化的二項式： 2 0 1 2 ( ) S r =a +a r+a r (5.3) 並且加入軸對稱的條件： 0 , 1 r= S = ；r 0 , S 0 r ∂ = = ∂ ；r=r0 ,S = 0 將上列條件代入(5.3)式，可推得： 2 0 ( ) MAX 1 r S r d r ⎛ _{⎛ ⎞} ⎞ ⎜ ⎟ = _{− ⎜ ⎟} ⎜ _{⎝ ⎠} ⎟ ⎝ ⎠ (5.4) z 梯形曲線(Trapezoid)： 在文獻 【22】 中，也提到了另一種振動的曲線，其壓電薄膜的位移近似為梯 形： 1 0 1 0 0 1 , ( ) , MAX MAX d if r r S r r r d if r r r r r < ⎧ ⎪ =⎨ _⋅ − _{< <} ⎪ ₋ ⎩ (5.5) 其中，r 為腔體的半徑；₀ r 為壓電薄膜的半徑。 ₁ 因此本研究中，採用(5.5)式中的梯形振動曲線來近似。 薄膜由水平掃過的最大體積∀ ，我們可以將梯型曲線對面積分求出： ₀ z 梯形曲線 將(5.5)式積分後得出： 0 0 0 1 + MAX MAX r r d dA d dA r r − ∀ = ⋅ −

∫∫

∫∫

2 3 2 3 1 0 0 1 1 0 1 2 1 1 1 6 2 3 MAX MAX d d r r r r r r r π π ⎛ ⎞ = ⋅ + _⎜ − + _⎟ − ⎝ ⎠ (5.6)

5.2 微幫浦簡易理論分析(Lumped System Analysis)

在文獻 【17】 中，Ullmann 對無閥式的微幫浦提出一套簡易的數學模式，主要針

對單一腔體的無閥式微幫浦，能估算出進出口的流量與腔體的壓力值。而對於雙腔串聯

內壓力的變化來求得。 在 Ullmann 所提出的數學模式中，根據其流動的方向，分成三種不同的模式來求解， 分別是排水模式、吸水模式以及過渡模式。為了簡化求解的過程，避免選擇使用何種模 式的困擾，本研究試圖將各個模式合併，以適合單腔與雙腔不同的模型。 5.2.1 單腔理論分析 Ullmann 對微幫浦提出的簡易數學模式，對於單一腔體的微幫浦，定義流動方向往 Diffuser 的方向為正，其方向皆朝向右，其流量及壓力的圖示如圖 5.2。對於幫浦流體流 動的方式，分成三種模式，分別是排水模式(Pin <Pout <Pc)，吸水模式(Pc <Pin <Pout)， 以及過渡模式(Pin <Pc < Pout)，其中P 為腔體壓力、c P 為入口端壓力、in P 為出口端壓out 力，而出口端壓力永遠大於入口端壓力，其背壓定義為P_out −P_in。 在漸縮/漸擴管元件中，流體的速度U與壓力差 PΔ 可以表示成含壓力損失 K 的方程 式： 1 2 2K Uρ = Δ (5.7) P 我們可以將其寫成流量與壓力差的形式： 2 2 2 _t K Q P A ρ _{= Δ (5.8)} 其中A 為漸縮/漸擴管的喉部面積。並將連續方程式式加入上述三種模式中，表示如下： _t z 排水模式(P_in <P_out <P_c)：此時腔體壓力大於出入口端壓力，流體在入口與出 口端皆為流出。 2 2 2 n in c in t K Q P P A ρ = − (5.9a) 2 2 2 d out c out t K Q P P A ρ = − (5.9b) Q_in +Q_out =∀₀⋅2πf sin(2πft) (5.9c) z 吸水模式(Pc < Pin < Pout)：此時腔體壓力小於出入口端壓力，相對於排水模式， 流體在入口與出口端皆為流入。

2 2 2 d in in c t K Q P P A ρ _{= − (5.10a)} 2 2 2 n out out c t K Q P P A ρ _{= − (5.10b)} −Qin −Qout =∀₀⋅2πf sin(2πft) (5.10c) z 過渡模式(P_in <P_c <P_out)：在此模式中，入口端的流動方向為流出，出口端為 流入。 2 2 2 n in c in t K Q P P A ρ = − (5.11a) 2 2 2 n out out c t K Q P P A ρ = − (5.11b) Q_in −Q_out =∀₀ ⋅2πf sin(2πft) (5.11c) 上述方程式中壓力損失係數 K 的下標 “n” 為漸縮管方向、下標 “d” 為漸擴管方向、 0 ∀ 為薄膜振動掃過的最大體積、 f 為振動頻率。除了Q 、in Qout、P 以外皆為已知值，c 因此可由上述模式求出不同模式下的壓力以及流量。 由 Ullmann 【17】所提出的理論中，將微幫浦流動的方式分成不同模式來討論，而 壓力損失係數則分為流經漸縮管方向以及漸擴管方向，但隨著腔體的增加，要找到相對 應的流動模式便越來越複雜，為了避免採用不同模式所造成的繁瑣，因此我們試圖將不 同模式下的方程式合併，不僅可適用於單一腔體的理論上，亦可應用於雙腔體的幫浦中。 (1) Model A(Steady Model)：首先，在 5.9~5.11 式中，流量與壓力差的關係式中壓力損

失係數 K 與流動的方向有關，向右為漸擴管方向而向左為漸縮管方向，所以利用流 量的正負變化，我們將其合併為一式，再加上連續方程式，合併後的方程式表示如 下： max( ,0) min( ,0) ₂ 2 in in d n in in in c in in t Q Q K K Q Q P P Q Q A ρ ⎡ ⎤ − = − ⎢ ⎥ ⎣ ⎦ (5.12a) max( , 0) min( , 0) ₂ 2 out out

d n out out c out

out out t Q Q K K Q Q P P Q Q A ρ ⎡ ⎤ − = − ⎢ ⎥ ⎣ ⎦ (5.12b) Q_out −Q_in =∀₀2πfsin(2πft) (5.12c)

利用 5.12~5.14 式中，流量的判斷便能夠適用所有的模式，並不需要再區分成不同的流 動模式，其中流動方向的判斷上，Q 、_in Q_out採用前一時間下的流量，為已知值， o in o in in in Q Q t Q t Q = ) ( ) ( ， o out o out out out Q Q t Q t Q = ) ( ) ( 其中 o in Q 、Q_outo 為前一時間下出入口的流量 上述式中的壓力損失係數K 、_d K ，在不同流量下有不同的值，因此我們將上一個_n 時間下的流量代入： ( ) ( ( ))

( )

o d d d K t =K Q t =K Q ；K tn( )=K Q tn( ( ))=Kn

( )

Qo 其中 o Q 為前一時間下的流量。 流量和壓力損失係數的關係我們採用文獻【15】中，呂學霖利用模擬漸縮/漸擴管元 件，得到流量和壓力損失係數的曲線，並且Curve Fitting 出的關係式： 7 0.921 1.315 10 0.5981 d K = × − Q− + (5.13) 6 0.8112 1.173 10 1.204 n K = × − Q− + (5.14)

(2) Model B(Unsteady Model)：因為幫浦的流動為隨時間變化的模式，所以腔體內流體

的流動會受到之前時間的影響，故考慮流體的慣性，在Steady Model 中加入慣性項 【19】，使其更符合物理現象。其方程式表示如下：

(

)

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

(

)

2 max( , 0) min( , 0) 1 2 2

out out out

d n out out Q out

out out t t l Q Q l dQ K K Q Q P P Q Q A _{A A} dt ρ _α ρ ⎡ ⎤ − + = − ⎢ ⎥ ⎣ ⎦ + (5.15b) ) 2 sin( 2 0 f ft Q Qout − in =∀ π π (5.15c) 其中A 為漸縮/漸擴管的擴張部面積、_l α_Q為權重係數、l為漸縮/漸擴管的長度。α_Q的值 介於 0~1 間，以調整慣性項所造成的影響。在判斷流動方向中，Q 、in Qout也採用前一 時間的流量代入，壓力損失係數則如同Steady Model 採用(5.13)、(5.14)式，將前一時間

的流量帶入求得。

在求解 Steady Model 及 Unsteady Model 上，除了Q 、in Qout、P 以外皆為已知值，c

因此可由聯立方程式組求出其解，另外在求解Unsteady Model 中的慣性項上，我們採用

全隱式差分(Fully implicit scheme)的方法來取代，

o dQ Q Q dt t − = Δ 在求解的過程中，必須假設初始條件，我們假設第一個時間下的出入口流量分別為壓電 振膜振動流量的1/2，腔體壓力P 則假設為c 0Pa，在幫浦振動一週期內則分成400 個 Time

Steps 來計算，初始條件(Initial Condition)表示如下：

9 3 0 1 2 sin(2 ) 2.16 10 / 2 in Q = − ∀ πf π f dt⋅ = − × − m s； 9 3 0 1 2 sin(2 ) 2.16 10 / 2 out Q = ∀ π f πf dt⋅ = × − m s； 1 400 dt f = P=0Pa； (3) Model C：除了加入慣性項，在文獻【17】中，Ullmann 也提到了壓電振膜材料所 產生的阻壓P (Blocking Pressure)。其振動薄膜所掃過的總體積(_b V )由於考慮P 的_b 影響，可表示如下式， atm b atm c P P P P ft V − − ∀ − − ∀ = ₀(1 cos(2π )) ₀ (5.16) 上式中，等式右邊第一項為振動薄膜所掃過的體積，而第二項則是考慮其材料剛性 所加入的，此項為負值。所以振動薄膜所掃過的總體積在考慮振動腔體材料性質與不考 慮振動腔體材料性質的情況下相比，其體積是較為減少的 因此我們也將 5.16 式對時間微分並加入連續方程式中，探討阻壓(Blocking Pressure) 所造成的影響，其方程式表示如下： 2 max( ,0) min( ,0) 2 in in d n in in in c in in t Q Q K K Q Q P P Q Q A ρ ⎡ ⎤ − = − ⎢ ⎥ ⎣ ⎦ (5.17a) 2 max( , 0) min( , 0) 2 out out

d n out out c out

out out t Q Q K K Q Q P P Q Q A ρ ⎡ ⎤ − = − ⎢ ⎥ ⎣ ⎦ (5.17b)