• 沒有找到結果。

噴嘴/擴散器式微幫浦之數值模擬及理論分析

N/A
N/A
Protected

Academic year: 2021

Share "噴嘴/擴散器式微幫浦之數值模擬及理論分析"

Copied!
113
0
0

加載中.... (立即查看全文)

全文

(1)

國 立 交 通 大 學

機械工程研究所

碩士論文

噴嘴/擴散器式微幫浦之數值模擬及理

論分析

Numerical Simulation and Theoretical

Analysis of Nozzle/Diffuser

Micropumps

研 究 生:呂學霖

指導教授:崔燕勇

(2)

噴嘴/擴散器式微幫浦之數值模擬及理論分析

Numerical Simulation and Theoretical Analysis of

Nozzle/Diffuser Micropumps

研 究 生:呂學霖 Student:Shiue-Lin Lu 指導教授:崔燕勇 Advisor:Yeng-Yung Tsui 國 立 交 通 大 學 機械工程研究所 碩 士 論 文 A Thesis

Submitted to Institute 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 2007

Hsinchu, Taiwan, Republic of China

(3)
(4)

噴嘴/擴散器式微幫浦之數值模擬及理論分析

研究生: 呂學霖

指導教授: 崔燕勇 博士

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

摘 要

本文是以計算流體力學的方法在分析噴嘴/擴散器式微幫浦內的流場與效 率。計算時使用非結構性和非交錯式網格;對於速度與壓力的偶合則採取 PISO 演算法則並以有限體積法將統御方程式離散化。而模擬此流場時是採用不同出口 壓力為邊界條件,並計算幫浦之淨流量。 本文針對 Ollson【5】的單腔式微幫浦進行分析,主要是探討不同背壓下幫 浦的性能。在研究中測試了不同的壓電薄膜的近似曲線對於流量的影響,結果顯 示較佳的近似曲線會使得模擬與實驗結果較為穩合,不同近似曲線會影響微幫浦 的淨流量。而當增加出口背壓時,流量與效率會隨之減少。 此外,本文採用了兩數學種模式來分析噴嘴/擴散器元件,找出了噴嘴/擴散 器元件的壓力係數Kn,Kd 與流量的關係,並且將其應用在 lumped mass 理論。透 過此簡易的理論,便可分析微幫浦不同背壓下的流量、效率及腔體中心壓力、。 本研究中比較了模擬和理論的結果,發現此理論估算的結果與模擬稍有差異,為 了解決此問題,本文發展出了不同的方式估算效率,以使得此理論在效率的估算 的準確性提高,而具有較高的可用度。

(5)

Numerical Simulation and Theoretical analysis of

Nozzle/Diffuser Micropumps

Student:Shiue-Lin Lu

Adviser:Dr. Yeng-Yung Tsui

Institute of Mechanical Engineering

National Chiao Tung University

ABSTRACT

In this thesis, CFD is used to analyze viscous flow and efficiency of micropumps . Unstructured and nonstaggered meshes are adopted in computations. The PISO algorithm is used to take the velocity-pressure coupling and the governing equations are discretized by the finite volume method. Pressure boundary condition was apply in our simulation.

In this thesis, a single chamber Micropump pump [5] is analyzed and we explore the pumping performance. The PZT is approximated by different mathematical curves. Testing indicates that the approximated curve need to be optimized to achieve better results. Different approximated curves will affect the pumping flow. When the back pressure is increased,the flow rate and pump efficiency get lower.

In order to find the pressure coefficient Kn&Kd of the Nozzle/Diffuser, we analysis the Nozzle/Diffuser element by using two different models.Furthermore,we use the lumped mass model to analyze the net flow and umping efficiency of the micorpump in different back pressures.We compare the results of the 3D simulation and the lump model.It was shown that there exits a little difference certain degree of difference between the two models.To solve this problem,we develop a different method to calculated the pump efficiency. This method improve the accuracy of the theorical solution and make Lumped Mass Model more useful .

(6)

誌謝 在交大這兩年來的日子,覺得成長了許多,其中最感謝崔燕勇教授在課業上 的指導與栽培,讓我除了對於知識成長之外,也學到了解決問題的方法與觀念, 在此特地對老師的教導獻上感激。 令外也感謝育昌學長和 George 在研究過程中的教導,以及同學飛鴉,彼此 相依為命與鼓勵,在平常互相嘴炮的過程中,無形地在課業中也成長了許多。也 感謝實驗室學弟欣恩、仕文、俊岩平常在生活的相處與幫助。 此外,感謝 TAMAGOYA 的伙伴們,感謝你們陪我渡過在新竹這兩年,大家 一起努力一起歡笑,讓我的生活更加豐富,沒有你們也許我也撐不過這兩年。最 後感謝我的父母與妹妹,在求學過程中一路支持,給我自主的空間,真的謝謝大 家。

(7)

符號表 符號: 定義: α 不同近似曲線的混合比 β 流量比 β1 排水模式下流量比的均值 β2 吸水模式下流量比的均值 0 r 壓電薄膜的半徑 1 r 振動腔體半徑 MAX d 壓電薄膜中心點最大位移量 f 壓電薄膜的振動頻率 ηη1η2η3 幫浦效率 pump P 背壓 i P 左端入口壓力 o P 右端出口壓力 c P 腔體中心點壓力 F 通量 m& 質量流率 P 壓力項 Re 雷諾數

u

X 方向之速度項

v

Y 方向之速度項 W Z 方向之速度項

V

v

速度向量 ρ 密度 μ 分子黏滯係數 Φ 相關變數

φ

流量比 θ Nozzle/Diffuser 流道夾角

(8)

W Nozzle/Diffuser 喉部寬度 L Nozzle/Diffuser 長度

H 微幫浦、Nozzle/Diffuser 深度

H

(9)

目錄 第一章、微幫浦簡介... 1 1.1 微流體系統簡介:...1 1.2 現存微幫浦的簡介...1 1.3 無閥式微幫浦的原理簡介...3 1.4 文獻回顧...3 1.5 研究內容...6 第二章、數值方法... 8 2.1 基本假設...8 2.2 統御方程式...8 2.3 壓電薄膜振動的近似...8 2.4 邊界條件(Boundary Condition)...9 2.4.1 壓電薄膜的邊界條件...9 2.4.2 出口的邊界條件...10 2.4.3 固定牆的邊界條件...10 第三章、動量方程式的離散 ... 11

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

3.2 有限體積法的離散...11 3.3 離散化...12 3.3.1 非穩態項 (Unsteady Term)...12 3.3.2 對流項 (Convection Term)...12 3.3.3 擴散項 (Diffusion Term)...13 3.3.4 源項 (Source Term)...14 3.4 線性代數式的整理...14 3.4.1 代數方程式...14 3.4.2 矩陣的求解...15 第四章、PISO演算法 ... 16 4.1 PISO演算法的簡介...16 4.2 速度與壓力之間的偶合...16 4.2.1 預測步驟(Predictor Step)...16 4.2.2 第 1 次修正步驟 (1st Corrector Step ) ...17

(10)

4.2.3 第 2 次修正步驟 (2nd Corrector Step)...19 4.2.4 壓力修正方程式的整理...20 4.3 進出口邊界條件的給定...20 4.3.1 出口壓力邊界的流量計算...21 4.3.2 出口壓力邊界的速度計算...21 4.4 固定牆邊界條件...22 4.5 PISO演算流程...22 第五章、無閥型微幫浦簡易理論... 24 5.1 壓電式幫浦的薄膜振動模擬...24 5.2 壓電式微幫浦的效率計算...26 5.2.1 真實效率ηreal的計算:...26 5.2.2 效率η 的計算:...28 1 5.2.3 效率η 的計算: ...29 2 5.2.4 效率η3的計算:...29 5.2.5 效率計算的比較:...30 5.3 Nozzle/Diffuser 元件的探討...30 5.3.1 穩態方形截面Nozzle/Diffuser探討...30

5.3.2 Region 2 的計算-Model A:...31

5.3.3 Region 2 的計算-Model B:...32 5.4 無閥型微幫浦的簡易理論(Lump Model )...33 第六章、結果分析與討論 ... 36 6.1 壓力邊界的驗證...36 6.1.1 T型管的測試...36 6.1.2 Y型管的測試...36 6.2 不同背壓下的 3D模擬結果...37 6.2.1 網格測試...37 6.2.2 考慮不同薄膜振動假設的流量結果...37 6.2.3 流場結構分析...38 6.2.4 幫浦效率探討...39 6.3 單腔微幫浦Lump-Model的理論分析:...40 6.3.1 Model A – 透過CFD 的計算...40 6.3.2 Model B – Yang的數學模式...40

(11)

6.3.3 Lump Model的結果分析與比較...41

第七章、結論 ... 43

Reference... 45

(12)

圖目錄 《圖 1.1》微幫浦在工程上的應用 ... 49 《圖 1.2》微幫浦的致動來源【1】 ... 49 《圖 1.3》壓電式微幫浦的運作方式... 50 《圖 1.4》無閥型壓電式微幫浦的發展史 ... 50 《圖 1.5》Ullman【6】所提的六種串聯方式 ... 51 《圖 1.6》前人對於微幫浦的模擬 ... 51 《圖 1.7》K-S Yang【14】所模擬多腔體方形微幫浦... 52 《圖 1.8》1998 年Ollson【5】所製微幫浦搆造示意圖 ... 52 《圖 1.9》1998 年Ollson【5】所製微幫浦之Nozzle/Diffuser元件尺寸 ... 52 《圖 1.10》1998 年Ollson【5】所製微幫浦尺寸... 53 《圖 2.1》本文所欲分析的微幫浦的幾外形。圖為 1998 年,Ollson【5】所製單腔 體圓形微幫浦示意圖... 54 《圖 2.2》薄膜中心點位移和速度隨時間的變化 ... 54 《圖 2.3》圓形壓電薄膜於不同半徑下的位移量(在上死點時)... 55 《圖 2.4》圓形壓電薄膜於不同半徑下的位移量(在下死點時)... 55 《圖 2.5》採用混合曲線下,在不同背壓下的壓電薄膜形狀 ... 56 《圖 2.6》本文模擬之幾何與邊界條件 (A)3D網格 (B) 邊界條件... 56 《圖 3.1》非結構性網格(Unstructured Grid) ... 57 《圖 3.2》Over-Relax Approach... 57 《圖 4.1》出口壓力邊界示意圖 ... 58 《圖 4.2》出口壓力邊界流量的計算... 58 《圖 4.3》邊界剪應力的計算示意圖... 59 《圖 4.4》PISO演算法流程圖 ... 60 《圖 5.1》壓電薄膜和幫浦內部的作用... 61 《圖 5.2》比較三種不同近似曲線於壓電薄膜的振動... 61 《圖 5.3》兩出口端流量隨時間的變化 (背壓=0) ... 62 《圖 5.4》方程式(5.18)式的積分範圍... 62 《圖 5.5》Diffuser-Element 示意圖,將Nozzle/Diffuser元件的壓力係數分成三個區 間探討:Region 1-入口區間,Region 2-漸縮/漸擴區間,Region 3-出口區間。 ... 63

(13)

《圖 5.7》(A)方管摩擦係數f與不同深寬比的關係 (B)以文獻【18】方式所求得管 道內部的壓力損失... 63 《圖 5.8》以文獻【18】方式計算Kn,Kd隨不同雷諾數(Re)的變化 ... 64 《圖 5.9》以文獻【18】方式計算於不同雷諾數(Re)下的ηnd... 64 《圖 5.10》Diffuser流向與Nozzle流向示意圖 ... 65 《圖 6.1》壓力邊界測試(1) - T型分岐管的幾何與網格... 66 《圖 6.2》壓力邊界測試(1)的結果–流量比( in m m Fi & &1 = )隨雷諾數的變化。... 66 《圖 6.3》壓力邊界測試(1)的結果–雷諾數為 10、200、400 的流場流線圖... 67 《圖 6.4》壓力邊界測試(2) – Y型分岐管幾何 ... 68 《圖 6.5》壓力邊界測試(2) – Y型分岐管網格 ... 69 《圖 6.6》壓力邊界測試(2)的結果 –流量比( in m m Fi & &2 = )隨出口壓力差dp的變化。 ... 69 《圖 6.7》壓力邊界測試(2)的結果-不同壓差下的流 ... 70 《圖 6.8》本文欲模擬的網格示意圖... 71 《圖 6.9》本文欲模擬的網格上視圖... 72 《圖 6.10》本文欲模擬的網格側視圖... 72 《圖 6.11》背壓 = 0 Pa,兩端出口流量Q2,Q1 隨時間變化 ... 73 《圖 6.12》背壓 = 2950 Pa,兩端出口流量Q2,Q1 隨時間變化... 73 《圖 6.13》背壓 = 5310 Pa,兩端出口流量Q2,Q1 隨時間變化... 74 《圖 6.14》背壓 = 0 Pa,每個週期下的淨流量(mL/min) ... 74 《圖 6.15》採用不同振動假設(2.3)、(2.4)式模擬所得淨流量與背壓關係圖... 75 《圖 6.16》比較不同振動假設(2.3)、(2.4)、(2.8)式模擬所得淨流量與背壓關係圖 ... 75 《圖 6.17》比較不同網格數量下(65288、117232、210258),淨流量與背壓關係76 《圖 6.18》所取流場截面示意圖 ... 76 《圖 6.19》背壓=2950 Pa,t=0.25T ,壓力場分佈圖 ... 77 《圖 6.20》背壓=2950 Pa,t=0.75T,壓力場分佈圖 ... 77 《圖 6.21》背壓=0,排水模式(t =0.25T ),截面圖(z=0.0001 mm) ... 78 《圖 6.22》背壓=0,排水模式(t =0.25T ),流場與壓力場截面圖(z=0.0001 mm) 78 《圖 6.23》背壓=0,吸水模式(t =0.75T),流場截面圖(z=0.0001 mm)... 79 《圖 6.24》背壓=0,吸水模式(t =0.75T),壓力場截面圖(z=0.0001 mm) ... 79

(14)

《圖 6.25》背壓=5310 Pa,排水模式(t=0.25T ),流場截面圖(z=0.0001 mm) ... 80 《圖 6.26》背壓=5310 Pa,吸水模式(t=0.75T),流場截面圖(z=0.0001 mm)... 80 《圖 6.27》背壓=0 Pa,排水模式,不同時間下Nozzle/Diffuser流場圖(z=0.0001 mm) ... 81 《圖 6.28》背壓=0 Pa,吸水模式,不同時間下Nozzle/Diffuser流場圖(z=0.0001 mm) ... 82 《圖 6.29》比較不同網格下,真實效率η 與背壓的關係...83 r 《圖 6.30》β(Q2/Q1)於一個週期內的變化... 83 《圖 6.31》Pump Model下,流量比β1 (Q2/Q1) 與背壓的關係(網格 = 117232).. 84 《圖 6.32》Supply Model下,流量比β2 (Q2/Q1) 與背壓的關係(網格 = 117232) 84 《圖 6.33》(5.19)式計算的效率η2與真實效率(η)在不同背壓下的關係 ... 85 (網格 = 117232) ... 85 《圖 6.34》(5.23)式計算的效率η3與真實效率(ηr)在不同背壓下的關係 ... 85 (網格 = 117232) ... 85 《圖 6.35》Kn,Kd隨流量的變化 (Model B)... 86 《圖 6.36》欲研究的Nozzle/Diffuser幾何尺寸 ... 86 《圖 6.37》Nozzle/Diffuser CFD模擬的網格上視圖與與測視圖... 87 《圖 6.38》Nozzle/Diffuser 模擬的邊界條件给定 ... 87

《圖 6.39》Nozzle/Diffuser模擬結果Ⅰ- 於不同進口流量(Q)下,Diffuser Direction 之壓力分佈圖 ... 88

《圖 6.40》Nozzle/Diffuser模擬結果Ⅱ - 於不同進口流量下(Q),Nozzle Direction 之壓力分佈圖 ... 89 《圖 6.41》壓降係數Cp隨不同流量下的變化 ... 90 《圖 6.42》Kn,Kd隨流量的變化 (Model A) ... 90 《圖 6.43》不同背壓下,腔體中心壓力Pc隨時間變化(一個週期內) ... 91 《圖 6.44》腔體中心點最大壓力Pc (T* = 0. 25) 與不同背壓下的關係 ... 92 《圖 6.45》腔體中心點最小壓力Pc (T* = 0. 75) 與不同背壓下的關係 ... 92 《圖 6.46》背壓= 0 Pa,右端出口流量(Q2)隨時間的變化(一個週期) ... 93 《圖 6.47》背壓=2950 Pa,右端出口流量(Q2)隨時間的變化(一個週期) ... 93 《圖 6.48》背壓=5310 Pa,右端出口流量(Q2)隨時間的變化(一個週期) ... 94 《圖 6.49》比較實驗、模擬、理論,來分析淨流量在不同背壓下的變化。 ... 94 《圖 6.50》以模擬、理論的流量結果,比較在不同背壓下的真實效率。 ... 95

(15)

《圖 6.51》背壓 = 0 Pa,比較Model A、Model B與真實模擬結果中,一個週期內 β(Q2/Q1)隨時間變化... 95 《圖 6.52》背壓 = 2950 Pa,比較Model A、Model B與真實模擬結果中,一個週期 內β(Q2/Q1)隨時間變化... 96 《圖 6.53》背壓 = 5310 Pa,比較Model A、Model B與真實模擬結果中,一個週期 內β(Q2/Q1)隨時間變化... 96 《圖 6.54》比較Model A、Model B與 3D模擬結果中,β (Q2/Q1)與背壓的關係。1 ... 97 《圖 6.55》比較Model A、Model B與真 3D擬結果中,β (Q2/Q1)與背壓的關係。2 ... 97 《圖 6.56》比較Model A、Model B與 3D模擬結果中,效率η 與背壓的關係。.98 2 《圖 6.57》比較Model A、Model B與 3D模擬結果中,效率η3與背壓的關係。. 98

(16)

第一章、微幫浦簡介

1.1 微流體系統簡介: 微流體元件及系統技術主要發展控制、感測、反應及分析微量流體。其元件 包括微閥、微幫浦、微流量計、微噴嘴微流道、微混合器等,並可整合為不同功 能之智慧化微型流體系統晶片。 在微流體系統,微幫浦為最關鍵的元件之一。在這幾年內發展出了多種型式 的微幫浦,而其也廣泛的應用在各個不同的領域,如圖 1.1 所示: 燃料電池的應用:閥式壓電微幫浦可應用在直接甲醇燃料電池(DMFC)裝置 燃料的傳遞上。燃料傳遞系統必須提供足夠的燃料,且又必須能排除同時產 生的二氧化碳。此系統可透過無閥式壓電微幫浦交換燃料和二氧化碳。 生物晶片上的應用:生物晶片是指在玻璃片、塑膠晶片等載體上,利用微機 電加工技術,放上特定的生物性材料,例如蛋白質或核酸,再經由一連串的 偵測反應,便可得到想要獲得的一些資訊以用來進行計算、比對等檢驗分 析。而其中推動流體的技術便需要用到微幫浦。先利用微機電技術在微流道 上放置樣本,在以微幫浦驅動,使液體在管道上進行混合等化學反應。 在支持冠狀動脈繞道移植上的應用:微幫浦技術可應用再心臟手術上,經由 實驗得知,由微幫浦支持可使在連續冠狀動脈阻塞時,維持在較好的血液動 力學狀態和心肌收縮性。

在藥物傳輸系統上的應用:藥物傳輸系統 Drug Delivery System (DDS)是結合 微電子系統將藥物直接植入體內。將藥物帶入體內循環系統,到達病變產生 的地方,再將藥物釋放出來。而壓電材料貼附在矽薄片上,便形成幫浦結構, 而微幫浦的流量控制再由頻率與控制信號來達成。 在熱傳上的應用---用於冷卻:微幫浦可應用在微縮實驗室晶片當晶片作用溫 度上升時,冷卻液走的是一條小循環路線,閥門式壓電微幫浦將冷卻液導入 熱交換區以達到冷卻效果。通過熱交換區溫度較高的液體再經由微幫浦送到 液態冷卻器冷卻,使得進入熱交換區的液體維持在低溫,完成一個循環,達 到冷卻效果。 1.2 現存微幫浦的簡介 在微機電的領域中,微幫浦被應用來當做為流體的輸送和控制。因此,如何 提高機械的驅動效率和精確的控制流體仍然是微流體系統中主要的研究目標。根 據微幫浦驅動方式,一般將其分為機械式和非機械式此兩類。機械式微幫浦的分

(17)

類如表 1.1 所示,以其作動方式可分好幾種型,如蠕動型(peristaltic pump)、轉軸 型幫浦(rotary pump)、擴散型幫浦(diffuser pump)和位移型幫浦 (displacement pump) 數種。機械式微幫浦原理是透過振動薄膜或者是其它可動元件週期性的運動以造 成流體輸送。因此機械式微幫浦包含可動式的元件(moving part)和微閥門系統 (micro-valve system)。 機械式微幫浦的分類 致動來源 腔室 閥門 Piezoelectric electrostatic electromagnetic pneumatic thermopneumatic Thermomechanic single chamber moving volumes peristaltic check valves valveless 表 1.1 機械式微幫浦的分類【1】 而由上表的分類,可知根據其可動元件的動力來源,較常見的數種如圖 1.2 所示: z 靜電式微幫浦(Electrostatic miropump): 一般的機械式微幫浦的運作可分為兩個模式(mode),對於靜電式微幫 浦,在吸水模式(supply mode)時,會在振動薄膜上施加電壓訊號,此時會造 成靜電吸引力而使薄膜向上移動,在此模式時,會將出口閥門關閉而進口閥 門打開,會造成流體會由進口流進腔室裡頭;在排水模式(pump mode)時, 則不在振動薄膜上通電,因此薄膜會恢復原狀而向下移動,在此模式時,會 將進口閥門關閉而出口閥門打開,因此流體會由出口流出去。 z 熱氣式微幫浦(thermopneumatic): 其原理是透過加熱驅動腔室內體積增大,而冷卻後,腔室內體積又回到 原來的大小,如此反覆運作達到驅動流體的效果。 z 壓電式微幫浦(PZT): 壓電式微幫浦是利用材料的壓電效應。壓電效應最早是在1880年由法國 居禮(Curie)兄弟於石英中發現。所謂的壓電效應簡單的說,就是在材料某兩 端施以壓力產生變形,則在材料的某特定兩端會產生正負電荷的分布,此效 應為正電壓效應。同樣的,在材料某兩端施以正負電場,則在材料的某特定 兩端將會產生變形,此即為負電壓效應。因此當我們在振動薄膜中使用壓電

(18)

材料,便可利用負電壓效應,施加電場使其上下振動,進而達到推動流體的 效果。本文所要討的微幫浦就是屬於此種類型。

令外一大類為非機械式微幫浦,這類的微幫浦運用感應電場來驅動液體,其 特色是是以固定的電極構造,藉由施加電壓產生電場來推動液體,而不需要任何 之可動元件。常見的有:電液體力學驅動幫浦(EHD)、電磁式微幫浦(MHD)、電 滲式幫浦(electro-osmotic pump)、電泳式幫浦(electrophoretic pump)與電濕式幫浦 (electro-wetting pump)等數種。 一般而言,機械式微幫浦具有構造簡單和容易驅動多數流體等特色。但傳統 的機械式微幫浦都需要止回閥系統,當微幫浦的振動薄膜做週期性振動時,微閥 門也必須相對應的開關以相互配合,因此在很高頻率的作動下,容易造成微閥門 的磨損和疲勞,而使得整個微幫浦系統損壞。因此,這幾年便發展出無閥式微幫 浦(valve-less Micropump),利用噴嘴/擴散器來取代原有閥門。 1.3 無閥式微幫浦的原理簡介 無閥式微幫浦與傳統的幫浦不同點在於利用漸縮/漸擴管(Nozzle/Diffuser)來 取代原本需要使用的微閥門。常見的無閥式微幫浦使用一上下振動的壓電薄膜, 而根據其運作方式,在一個週期內可分為兩個模式討論,如圖 1.3 所示: z 吸水模式(Supply Mode): 如圖 1.3A,在此模式時,壓電薄膜向上移動,此時腔內的體積增大,流體會 經由兩端的 Nozzle/Diffuser 元件流進腔體,進口端的流體經由 Diffuser 方向流進, 壓力損失較小;而出口端是經由 Nozzle 方向流進,壓力損失(Pressure Loss)大,因 此從進口端進入的流量會大於出口端。 z 排水模式(Pump Mode): 如圖 1.3B,在此模式下,壓電薄膜向下移動,此時腔體內被壓縮而體積變小, 便將流體往兩端同時排出,此時進口端是走 Nozzle 方向流出,壓力損失大;而 出口端則反之。因此在此模式下,從出口端排出的流量會大於進口端。因此,如 此反覆作動,便可達到輸送流體之功能。 1.4 文獻回顧 在本文中,文獻回顧主要分為無閥型微幫浦的製造實驗與幫浦理論的建立、 模擬與 Nozzle/Diffuser 的研究。 1989 年,荷蘭人 Van de pol【1】最早提出無閥門的構想。其主要概念為在出 口端和入口端可利用噴嘴/擴散器,可造成不同的壓降,便可取代傳統止回閥, 達成無需閥門便可控制流體流向的構想。

(19)

1993 年,Stemme et al.【2】以黃銅加工出了第一個無閥型微幫浦。幫浦腔直 徑為 19mm,而進出口均 Diffuser 來取代閥門,而 diffuser 角度均不超 20 度,振 動薄膜則是用圓形的薄膜片(圖 1.4A)。當時做出的流量最大可達 16ml/min。 1995 年,Gerlach et al.【3】以半導體製程的方式,在矽板蝕刻出大角度的漸 縮管。其作動方式都先前的一樣,不同點在流體的流動方向恰和先前的相反,是 以 nozzle 為正流動方向而不是以 Diffuser 為流動方向(圖 1.4B)。他們也發現當大 角度時會有回流,且在 Re>100 時,亦會產生。 1995 年,Ollson et al.【4】首度嘗試製造出平面式圓型腔體無閥微幫浦,並 且將兩個不同腔室並聯,其腔體直徑均為 6mm,而在兩個腔體上面均有壓電薄 膜振動。在此篇研究中也針對無閥型微幫浦的效率與流量做分析。 1997 年,Ollson et al.【5】以熱壓法製造圓形單一腔體的微幫浦,並且以 7 度張角的 Diffuser 元件連接腔體與進出口(圖 1.4C)。其結果顯示壓電薄膜在頻率 2200Hz 時最大流量可達 1ml/min,而最大背壓則有 5900kpa。令外,其實驗中也 並聯了兩個相同腔體比較(圖 1.4D),且最大流量可高達 1.9ml/min。 而在 Ollson 提出一系列微幫浦的製作之後,關於無閥型微幫浦製造與實驗的 相關文獻研究如雨後春筍般的發展。然而更進一步,提出對於微型幫浦數學模式。 1998 年,Ullmann【6】建立了微幫浦數學模式,他將微幫浦腔內的振動,在 一個週期之內切成三個區間來探討,在這三個區間內,透過壓力係數 K 來估算 其流量。令外,此篇研究中也提出了微幫浦可能的六種串接方式(圖 1.5):

1997 年,Ollson et al.【7】提出了塊狀模式來分析 Nozzle/diffuser 效果。在其研 究中針對壓電形的無閥式微幫浦,應用流體力學中的連續方程式和伯努力方程式 方式及整合出了三個方程組,並以 MATLAB 求解此方程組得到出入口流量、入 口流量和振動腔的壓力值。 2001 年,Pan et al.【8】 提出了微幫浦內振動薄膜與振動腔內壓力之間的偶 合關係式。並在三種不同的壓力差的情況下,估算出微幫浦的流量的解析解。 而關於微幫浦的模擬研究的相關文獻則較少。最早起步是在 1997 年開始: 1997 年,Athavale et al.【9】針對 1993 年 Stemme【2】所製做的微幫浦幾何 做 流 場 的 模 擬 ( 圖 1.6A) , 在 壓 電 薄 膜 振 動 的 模 擬 上 利 用 CFD-ACE 中 的 Fluid-Thermo-Structure (FSI)模組,此模組是用來處理液固界面的問題。其模擬出 的最大流量與實驗相當穩合,但最大背壓則差了一點。

2000 年,新加坡 et al.【10】等人以 CFD-ACE 模擬圓型腔體無閥式微幫浦 (圖 1.6B),其幾何類似於 1997 年 Ollson【5】所製的單腔幫浦。其壓電薄膜的位

(20)

移曲線模擬採用 Timoshenko 的圓板振動假設,並以方波訊號的 FOURIER 級數做 為振動函數的假設。 2004 年,羅卓錚【11】以 CFD-ACE 針對 1997 年 Ollson【5】的單腔幫浦做 模擬,其論文中壓電薄膜位移以兩種不同曲線假設,並以三角函數做來回振動, 模擬出的最大流量與實驗約有 30%的誤差。 2005 年,Fan et al.【12】利用有限體積法軟體 CFDRC 針對壓電無閥微幫浦 做分析(圖 1.6C),其包含了流場分析、壓電薄膜振動和電流這三個元素,並整合 此三者關係的統御分程式:Navier-Stoke 方程式、質量守恆方程式及薄膜振動方 程式【8】。而結果顯示出薄膜的振動頻率在 7500Hz 之內,流量會隨著頻率而上 升;而超過此頻率後,則流量反而會隨著頻率上升而下降。 2005 年,曾裕博【13】以 Fluent 模擬方形腔體的微幫浦,其論文建立了關 於 Nozzle/Diffuser 效率與幫浦效率之間的數學關係,並對不同張角的分析其效能。

2006 年,Yang et al.【14】以 CFD-ACE 分別模擬單一腔體、雙腔體串聯與單 腔體並聯的微幫浦(圖 1.7),並分析多腔體微幫浦在不同相位角下的流量結果。其 結果顯示,串聯式的微幫浦在兩個腔體振動的相位差達 100 度時,有較佳的流量。

2007 年,韓國 Jeong et al.【15】等人比較 FSI model 與 Timoshenko 假設,在 低頻的情況下以此兩種不同壓電薄膜的位移的方式的差異,而幫浦的幾何形狀則 類似於 1993 年 Stemme 所製微幫浦。其結果可知此兩種假設位移形狀不太類似, 所得出的流量也不盡相同。 僅管針對微幫浦全幾何的流場模擬上,可透過 CFD 的計算以得到估計的淨 流量,但在計算上則需耗費大量時間。因此根據 Stemme【1】所提出的幫浦理論 中可知,影響幫浦效率的主要關鍵在於 Nozzle/Diffuser 元件,因此亦有文獻針對 於 Nozzle/Diffuser 元件效率的探討:

1995 年,Ollson et al.【16】針對 Nozzle/Diffuser 元件提出了壓力係數的計算 方式。其原理是將此元件分為三個區間: z 入口區:因流體流進 Nozzle/Diffuser 元件的壓力損失 z 漸縮/漸擴區:因流過 Nozzle/Diffuser 而造成的壓力損失 z 出口區:當流體流出 Nozzle/Diffuser 元件的壓力損失 並利用壓降係數進而求出其工作效率,而結果顯示出採用漸縮管為通道的效 果較佳。

2000 年,Ollson et al.【17】以有限元素法的軟體 ANSYS,以數值模擬的方式 (CFD)計算 Diffuser 的流場,並比較先前實驗的結果。他們先前實驗的結果是在

(21)

Diffuser 以不同的角度:7、9.8、13,其效率會隨著其角度的增加,或者是長度的 拉長而有較佳的表現。在模擬中,分別以層流模式和紊流模式方式去模擬,但結 果與先前實驗有著一定的誤差。

2003 年,Yang et al.【18】針對不同張角的 Nozzle/Diffuser 實驗,比較在不同 雷 諾 數 下 的 Diffuser Efficiency , 並 以 方 管 的 摩 擦 力 損 失 的 數 學 計 算 來 預 測 Nozzle/Diffuser 的效率,其結果顯示在 5 度張角的時後,實驗結果與其預測相當 穩合,但 20 度張角的結果卻有相當的誤差。

2004 年,Singhal et al.【19】於低雷諾數情況下,以套裝軟體 fluent 對 nozzle/diffuser 進行分析。其所做的幾何外型為平板式漸擴管和圓錐式的漸擴管。 並在入口邊界給定 Fully Develop 條件,其結果顯示出在低雷諾數下,壓力係數和 Re 數有很大關係,尤其在小角度時。令外也發現圓錐式和平板式的漸擴管結果 相似。 1.5 研究內容 本文研究中,主要是以 1997 年 Ollson et al.【5】所實驗的微幫浦為參考幾 何,如圖 1.8,利用計算流體力學的方式來分析此型微幫浦的效果。其幾何幾形 狀為一個深度 0.2mm 的單一圓形腔體的微幫浦,腔體直徑為 3mm。腔體上放置 直徑 2mm 的壓電振動薄膜片。在腔體與進出口間以兩個 Nozzle/Diffuser 元件連 結,Nozzld/Diffuser 元件的張角為 7 度,長度為 1mm,而喉部的寬度為 0.1mm, 如圖 1.9 所示。 本研究主要是以不同的曲線來模擬幫浦作動時的流場現象,並探討流量與幫 浦背壓之間的關係,其中兩端背壓的壓差變化從 0~5900Pa,振動頻率 f =2200Hz 並且更進一步找出最佳的近似曲線來模擬幫浦的流場。 而為了簡化時間時,本文模擬中採用對稱邊界條件,僅模擬一半的微幫浦, 其幾何尺寸如圖 1.10 所示。令外,本文研究中以三種不同數量網格的來探討網 格對於流量結果的影響。 而透過對流量的分析,便可更進一步分析幫浦的效率。本文研究中亦探討在 不同背壓下,幫浦效率的變化。本文中以流量比參數β取代前人的ηnd,並參考 前人的方式,建立了數種不的效率的估算方式,以期能透過對 Nozzle/Diffuser 的 分析便可得知幫浦的效能。 令外,我們也就 Nozzle/Diffuser 元件做探討,本文中採用兩種不同的數學模 式來分析 Nozzle/Diffuser 元件,並建立壓力係數K ,n Kd和流量之間的關係。 最 後 , 我 們 透 過 Ullmann 【 6 】 所 提 的 數 學 模 式 , 並 整 合 先 前 所 建 立

(22)

Nozzle/Diffuser 數學模式,建立起一套微幫浦理論 Lump Model,以理論解求得β 值,並將其代入先前建立的效率估算式子中,更進一步得到幫浦效率,並以此與 模擬的結果相比較來驗證其準確性。

(23)

第二章、數值方法

本章的內容主要敘述我們求解流場所使用的數學模式。在 2.1 節我們先說明 對於本次模擬的基本假設。在 2.2 裡敘述我們求解流場所使用的統御方程式。2.3 節裡頭敘述對於壓電材料薄膜的振動位移所採用的近似曲線函數。最後,對於邊 界條件的給定則會在 2.4 節內敘述 2.1 基本假設 本研究中,我們模擬的壓電幫浦的幾何形狀如圖 2.1,其所使用的工作流體 為水,故我們可假設流體為: z 不可壓縮流( Incompressible Flow):我們假設水的密度為常數。 z 忽略重力項 (Neglect body force)

z 非穩態流場 (Unsteady): 在本次模擬中,腔體底部的壓電材料的振動會隨時間變化,因此流場的 變化亦隨時間而改變,當壓電薄膜向上振動時時,流體會經由水管向外流 出;而當往下時,流場的方向則反之。固我們假設流場為非穩態的情況。 2.2 統御方程式 z 連續方程式(Continuity Equation):

( )

=0 ⋅ ∇ ρV (2.1) z 動量方程式(Momentum Equation):

( ) ( )

VV P

( )

V t V ∇ ∇ + −∇ = ⋅ ∇ + ∂ ∂ρ ρ μ (2.2) 其中μ: 黏滯係數,ρ:流體密度 ,∇P:壓力梯度。 2.3 壓電薄膜振動的近似 在本文研究中,對於壓電薄膜的振動位移是採用兩種不同的曲線【11】,分 別以梯形與二次曲線來近似,而薄膜隨時間的來回振動,則以 Euler Beam 的假設 為三角函數(圖 2.2): z 二次曲線(Parabolic Curve): ) 2 ( 1 ) , ( 2 0 2 1 d Cos f t r r t r Z MAX ⋅ ⋅ ⎦ ⎤ ⎢ ⎣ ⎡ − − = π (2.3) 圖 2.3,圖 2.4 所示的分別為壓電薄膜在上死點與下死點時的位移。 z 梯形曲線( Trapezoid Curve ):

(24)

= ) , ( 2 r t Z 0 1 1 0 0 1 , ) 2 ( , ) 2 ( r r r IF t f Cos d r r r r r r IF t f Cos d MAX MAX < < ⋅ ⋅ ⋅ − − − < ⋅ ⋅ − π π (2.4) 其中r 為腔體半徑,0 r1為壓電薄膜的半徑。dMAX為中心點由水平至上死點的 最大位移量,本文假設dMAX=1.0μm。 f 為振動頻率,依其實驗設定為 f =2200Hz。 另外在模擬中,我們也考慮兩個曲線做線性混合,以期能在不同背壓的情況 下都能得到較佳的結果。在不同背壓下,以α 值來線性混合兩種曲線,其中α 值 為 0~1 之間,而位移變化如圖 2.5 所示。 z 混合曲線( Blending Curve ): ) , ( ) 1 ( ) , ( ) , ( 1 2 3 r t Z r t Z r t Z =α⋅ + −α ⋅ (2.5) 當 ) ( C , 1 ) ( , 0 二次曲線 梯形 urve Parabolic Trapezoid = = α α 而幫浦背壓則為Ppump =α⋅PMAX。 2.4 邊界條件(Boundary Condition) 在本文模擬中,採用對稱邊界條件,因此僅模擬真實微幫浦幾何的一半,其 網格如圖 2.6A 所示。而邊界條件如圖 2.6B,圖中 Inlet 與 Outlet 的部份均採用壓 力邊界條件,壓電薄膜邊界條件的給定在接下來下一節中說明。 2.4.1 壓電薄膜的邊界條件 由於在邊界上有壓電薄膜的振動,因此必需將其邊界設為移動邊界,而流體 的計算上也必需藉由移動網格的計算。在文獻中【13】中,曾等人對於壓電薄膜 的振動模擬,是以進口速度取代。而其結果顯示當振幅小於 1000nm 且振動頻率 低於 3000Hz 之下,以進口速度來取代薄膜振動其與採用移動邊界做動態網格的 計算相比,誤差在可接受範圍之內。故在本文研究中對於壓電薄膜的振動是以一 進口速度來取代。 給定的進口速度即為壓電薄膜的振動速度,因此將方程式(2.3)-(2.5)式對時間 微分後可得到進口壓電薄膜的振動速度: „ 二次曲線( Parabolic Curve): ) 2 ( 1 2 ) , ( ) , ( 2 0 2 1 1 d Sin f t r r f t t r Z t r V MAX ⋅ ⋅ ⎦ ⎤ ⎢ ⎣ ⎡ − ⋅ = ∂ ∂ =

π

π

(2.6) „ 梯形曲線( Trapezoid Curve ):

(25)

= ) , ( 2 r t V 0 1 1 0 0 1 , ) 2 ( 2 , ) 2 ( 2 r r r IF t f Sin d r r r r f r r IF t f Sin d f MAX MAX < < ⋅ ⋅ ⋅ − − ⋅ < ⋅ ⋅ ⋅

π

π

π

π

(2.7) „ 混合曲線( Blending Curve ):

(

1

)

( , ) ) , ( ) , ( 1 2 3 r t V r t V r t V =α⋅ + −α ⋅ (2.8) 2.4.2 出口的邊界條件 出口邊界設定我們採用的是壓力邊界,本文所研究的幾何共有兩個出口,並 考慮兩端出口在不同壓力差下流量的變化。壓力邊界的數學式子和流量的計算會 在稍後的 4.3 節中介紹。 2.4.3 固定牆的邊界條件 固定牆的物理意義為 No-Slip Condition,因此邊界為固定牆鄰近網格需考慮 剪應力的影響。其中詳細的數學式子在稍後的 4.4 節中介紹。

(26)

第三章、動量方程式的離散

在本章裡,我們會離散前一章所提到的統御方程式及採用非結構性網格所使 用的數值方法。3.1 節會先介紹有限體積法的概念,3.2 節裡應用有限體積法來處 理我們所使用的動量方程式。3.3 節會進一步詳細的描述每一個項的離散化。最 後在 3.4 節裡,動量方程式將被離散成代數方程式,面對此線性代數的問題,我 們以 BICG 方法來求解此矩陣。

3.1 有限體積法 (Finite Volume Method)

過去 30 年來,隨著電腦的發展,越來越多人應用數值方法來求解工程問題。 而一些求解偏微分方程式的方法也被發展出來,如有限元素法(FEM)、有限體積 法(FVM)、有限差分法(FDM)…等。現今,有限元素法是最廣泛應用在求解偏微 分方程式的方法。但面對特定的工程問題,尤其是在對於需要處理通量的問題, 有限體積法仍是目前較為可靠的方法。在本文研究中,將採用有限體積法。對於 計算流體力學而言,以有限體積法來解 Navier-Stokes 方程式比其他數值方法來 說,更較有物理的意義,而其亦可應用於非結構性網格(圖 3.1)。在接下來的文章 中,將敘述如何將有限體積法應用於非結構性網格上。 3.2 有限體積法的離散 在前一章裡已敘述所使用的統御方程式(2.2 式),我們可將其寫成下列通式:

( )

( )

(

)

q V t +∇⋅ =∇⋅ ∇ + ∂ ∂

ρφ

ρ

φ

μ

φ

(3.1) 將其積分後,其積分式可表示如下:

∫∫∫

∫∫∫

∫∫∫

∫∫∫

∀ ∀ ∀ ∀ ∀ + ∀ ∇ ∇ = ∀ ⋅ ∇ + ∀ ∂ ∂ qd d d V d t

ρφ

(

ρ

φ

) (

μ

φ

) (3.2) 【非穩態項】 【對流項】 【擴散項】 【源項】 其中φ為我們所要計算的速度向量,而 q 為源項,其來源為壓力項−∇1P。 透過高斯散度定理的轉換後,對流項和擴散項將由體積分轉換為面積分,方程式 將被整理如下 :

∫∫∫

∫∫

∫∫

∫∫∫

∀ ∀

+

=

+

qd

S

d

S

d

V

d

t

ρφ

S

ρ

φ

S

μ

φ

(3.3) 其中 q 為源項其來源包含壓力項−∇P,故q=−∇P,。 :μ 為流體黏滯係數。 在本文的求解中,我們將壓力梯度視為已知數,故將其放入源項。在下一節

(27)

裡,將會進一步的敘述每一個項的離散方法。 3.3 離散化 3.3.1 非穩態項 (Unsteady Term) 對於有移動邊界的不可壓縮流,其網格的體積會隨著時間而改變。因此非穩 態項將被離散成如下式 :

( )

(

o

)

t d t φ φ ρ ρφ − Δ ∀ Δ ⋅ = ∀ ∂ ∂

∫∫∫

∀ (3.4) 其中 φ,φo 分別為我們要求解的變數其在這一個時間和前一個時間的當下 值。因此右手邊的第二項即為前一個時間的值,即為已知數,故在計算上可將其 移進源項。 3.3.2 對流項 (Convection Term) 經由高斯散度定理的轉換後,原本對流項從體積分轉為面積分,而對流項可 被離散為下式:

(

)

∫∫

⋅ = ⋅ = = f C f f f f f f f f f S f V S m F S d V

φ

ρ

φ

φ

ρ

v v v v & (3.5) 其中

S

v

f :面之法向量。m& 代表流經面上的質量流率。 f 而下標 f,代表的是 ρ,V,φ 在面上的值。對於其在面上的值,我們可透過 其相鄰網格以線性內差的方式求之。 在本文的研究中,求解對流項所採用的方式為一階上風法和中央差分法的混 合式。因此,對流項可整如下:

( ) (

UD

)

f CD f UD f C f

F

F

F

F

=

+

γ

(3.6)

【New】 【Old time level】 其中γ 的值介於 0~1 之間,其表示之意義為: scheme difference central scheme upwind , 1 , 0 = = γ γ 在本文研究中採用的γ =0.9。而右手邊的第二項我們將以前一個時間項的值, 故將其移至源項中計算。 z 在一階上風法裡,φf 的值取決於其上游:

(

)

(

f

)

f c f p f f f m m m m φ φ φ φ = < = > , 0 , 0 & & & & 由外流向主格點 由主格點流向外

(28)

z 在中央差分方式裡,φf 的值則為其相鄰網格的內差,其表示如下:

(

)

p c f

w

φ

w

φ

φ

= 1

+

(3.7) 其中 w 為權重係數,w 值在 0~1 之間。 因此,最後(3.6)式可寫成如下列式的型式:

[

]

[

]

[

f p C f P f C

]

C f P f c f m m w w m r m m F

φ

φ

φ

φ

φ

φ

) 0 , max( ) 0 , max( ) 1 ( ) 0 , max( ) 0 , max( & & & & & − − − + − + − + = (3.8) 3.3.3 擴散項 (Diffusion Term) 經由高斯散度定理的轉換,擴散項由體積分轉為面積分,將其離散後可得下 列式:

∫∫

=

=

=

f D f f f f f s D

F

S

S

d

F

μ

φ

v

μ

(

φ

)

v

(3.9) 其中μf 為於此面上的黏滯係數,計算上是以其所相鄰的網格做線性內插。 針對非結構性網格而言,對於擴散項的處理,我們利用 Over-Relax 的方法 來近似其原來的面法向量。 因此,首先我們將原來的面法向量 S 分成兩部份, 其表示如下:

(

S

d

)

d

S

f

=

+

f

(3.10) 上式右手邊的第一項我們稱其為正交向量,而第二項為非正交向量。而正交向量 d 的訂義如下:

δ

δ

v

v

v

v

v

f f

S

S

d

=

2 (3.11) 其中δv定義如圖 3.2 所示。將 (3.10)代入(3.9),擴散項可整理成下式:

)

(

)

(

2

d

S

S

S

F

C P f f f f f f D f

r

r

r

v

r

Γ

+

=

φ

φ

φ

δ

μ

(3.12) 【New time level】 【Old time level】

在上式中,右手邊的第二項為非正交項,在計算上,∇φf 以前一個時間項所 計算的值近似,因此可其將其放至源項。

(29)

3.3.4 源項 (Source Term) 在計算中,我們將壓力梯度視為已知的值,因此可將放至源項。因此源項的 來源即為壓力梯度−∇P。對於非結性網格而言,在求解壓力梯度項時,會需要 利用邊界上的壓力值,為了得到邊界上的壓力,我們可做以下的推導:

δ

v

=

P

P

P

b p (3.13) 其中,b 表示壁面上的中點;而

δ

v

代表 P 到 b 之距離向量。 如前述所示,壓力梯度可由下方的式子去近似之,

⎟⎟

⎜⎜

+

Δ

=

Δ

=

≠ b f f f b b f f f

S

P

S

P

S

P

P

1

v

1

v

v

(3.14) 此處下標

f

表示除了

b

之外的其餘各面。 因此我們得到了邊界上之壓力為: ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∀ Δ − ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⋅ ∀ Δ + =

≠ δ δ v v v v b b f f f p b S S P P P 1 1 1 (3.15) 3.4 線性代數式的整理 3.4.1 代數方程式 透過先前幾小節的敘述,我們將非穩態項、擴散項、對流項及壓力項合併之 後,可表示成下列的線性代數式:

+

=

A

Q

A

P

φ

pn C

φ

C (3.16)

t

A

A

P C

Δ

Δ

+

=

ρφ

(3.17) 其中

φ

P

,

φ

C 為此線性代數式的變數,而AC, 為: Q

)

0

,

max(

2 f f f f C

m

S

r

S

A

r

&

r

+

=

δ

μ

(3.18)

[

]

[

]

t 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 ( r r & & &        (3.19) 其中AC項的值僅來自於對流項和擴散項。而源項 Q 則包含從非穩態項、對 流項、擴散項和壓力項貢獻得來。

(30)

3.4.2 矩陣的求解 透過(3.15)式的表示,動量方程式被表示成求解線性代數的問題。在工程上, 對於求解矩陣的方式,大致可分為直接方式和疊代方式這兩種類型的方法。直接 方式來求解矩陣需要較大的記憶體空間來儲存矩陣。而對於非結性網格系統,其 產生的係數矩陣為稀疏矩陣(Sparse Matrix),而網格數量的大小會決定矩陣行列 數,當網格數量到達好幾萬時,電腦記憶體便無法負荷這麼龐大的矩陣。而疊代 式的求解則不需要如此大的記憶體空間。因此,在本文研究中對於求解矩陣,採 用疊代方式的 BICG 和 ICCG 來求解。

(31)

第四章、PISO 演算法

在前面一個章節裡,敘述了動量方程式的求解方式。但在物理上,其並未滿連續 方程式。為滿足質量守恆,可以從離散過後的動量方程式裡頭找出速度與壓力的 關係式,並藉此修正壓力和速度值。在今日,已發展出幾種壓力和速度的偶合方 式,如 SIMPLE、SIMPLEC、PISO 等。在本文的研究中,我們採用 PISO 此種方 式。本章節主要是簡介 PISO 及壓力修正方程式的離散。在 4.1 節我們會簡介 PISO 的特性。4.2 節會進一步敘述其離散過後的方程式。4.3 節中會敘述本文所使用的 定壓邊界條件。4.4 節則會描述對固定牆邊界的數學式子。最後,整個 PISO 演算 的流程將會整理在 4.5 節。 4.1 PISO 演算法的簡介 在本文中,採用 PISO 演算法使我們能得到流場的至暫態解。PISO 的優點在 於此演算法對於壓力和速度之間的偶合關係,不需經由反覆的疊代。因此對於流 場的暫態解,採用 PISO 可節省許多時間。

PISO 演算法最早是由 Issa【24】於 1986 年提出。PISO 四個字代表意思為 (Pressure Implicit with Splitting of Operator),其求解流程包含兩個步驟,我們會在 接下來的文章中敘述其如何使用於非結構性網格的運算上。 4.2 速度與壓力之間的偶合 4.2.1 預測步驟(Predictor Step) 在前一個章節裡,透過求解動量方程式而得出了速度場V*,而在求解過程 中其壓力是採用一個現有(guessed)的壓力值p ,但其並未滿足連續方程式,因此* 我們必須透過修正壓力值以使其滿足。為了滿足連續方程式,首先必需計算面上 速度。我們可以先將動量方程式中壓力項自源項提出,表示如下列型式: * * * * P P P c p p

P

A

V

H

V

⎟⎟

⎜⎜

⎛ ∀

Δ

=

(4.1) 其中 P C C c p A Q V A V H ⎟=

+ ⎠ ⎞ ⎜ ⎝ ⎛ * * * (4.2) 同樣地,控容面上之速度與壓力關係式可透過(4.1)式的近似表示為: * * * * f f P f f P A V H V ⎟⎟ ∇ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = r (4.3 )

(32)

其中 f f P f f

P

A

V

H

⎟⎟

⎜⎜

⎛ ∀

Δ

+

=

*

r

,故將其代入上式可得: * ) ( f f P f f P f f P A P A V V ⎟⎟ ∇ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ + = (4.4) 而在(4.4)式中 f P A ⎟ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ 的計算,是以主控面上相鄰的兩網格取平均值,即: ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ + ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ C P P P f P A A A 2 1 (4.5) 因此面上的質量流率可由下面的方式改寫為: d P P A S V S P P A 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 ⋅ ∇ − ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ⋅ ⋅ ≈ ⋅ ∇ − ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ⋅ ⋅ = ⋅ ⋅ = ) ( ) ( * * * * ρ ρ ρ ρ ρ r r r r &

(

)

[

P P P r

]

S r S A S V C P f f P f f f f δ δ ρ ρ − −∇ ⋅ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ⋅ ∀ Δ − ⋅ ⋅ = r r v 2 (4.6) 4.2.2 第 1 次修正步驟 (1st Corrector Step ) 在前一個步骤中已得到的速度V*和壓力P 在接下來的文章裡我們會敘* 述如何計算這些修正量。於 PISO 演算法中,第一次修正步驟中主格點與相鄰格 點速度與壓力關係式可表示為: * * * * * * p P P c P P

P

A

V

H

V

⎟⎟

⎜⎜

⎛ ∀

Δ

=

(4.7) 其中

V

**

=

V

*

+

V

'

,

P

**

=

P

*

+

P

'

。其中的上標**和*表示的即為修正後與修正 前的值,上標*的值是由求解動量方程式得到,而上標 '即為我們所要修正的修 正量。 將方程式(4.7)與(4.1)相減後,主格點中修正速度V 和修正壓力梯度p' ∇ 可表P' 示為下列的關係式:

(33)

f P P P p

P

a

V

⎟⎟

⎜⎜

⎛ ∀

=

'

(4.8) 其中的∇PP' =∇Pp** −∇Pp*,上述的方程式是在網格主格點上的推出,其亦可使 用在每個主控面上。因此在每個主控面上的新流量的修正值可表示如下:

'

* * * f f f

m

m

m

&

=

&

+

&

(4.9)

其中 m& :未修正前的流量,其計算上是透過方程式(4.9)式。 f* ' f m& :要修正的流量值。即m&f = − fVpSrf ⋅ ' ' ρ 。流量修正量可表示成下列式: .

(

S d

)

P A d P A S P A m f f f P f f f P f f f f P f f − ⋅ ′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ⋅ ′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − = ⋅′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − = ⋅ r r &

ρ

ρ

ρ

' (4.10) 其中 f f

S

S

d

=

δ

2 ,其定義在前一章節中我們已敘述過.。將(4.10)式代入(4.9) 式中,方程式可改寫為:

(

S

d

)

P

A

d

P

A

m

m

f f P f f f P f f f

⎟⎟

⎜⎜

⎛ ∀

Δ

⎟⎟

⎜⎜

⎛ ∀

Δ

=

*

ρ

ρ

* *

&

&

(4.11) 在網格主格點 P 上,我們可將方程式(4.11)式用於在網格的每個控制面上 (Control Surface)。為了滿足質量守恆,因此對於網格的控制面上,其修正過後總 合流量為應為 0,固可表示為:

=

+

=

f f f f f f

m

m

m

&

**

&

*

&

'

0

(4.12) 最後,將方程式(4.11)代入(4.12)中,方程式可整理成下列的線性代數式:

(

)

⎟∇ ′⋅ − − ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∀Δ + ′ = ′ f f f f p f f c c P p P S d m A P A P A

ρ

& * (4.13)

=

f C P

A

A

(4.14)

(34)

其中 f f f p f c S S A A ⋅ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∀Δ =

δ

ρ

2 。 4.2.3 第 2 次修正步驟 (2nd Corrector Step) 在修正步驟 1 時,我們透過求解第一次壓力修正方程式(4.13),而得到了第 一次修正過後速度和壓力值。在《第 2 次修正步驟》裡。類似於《第 1 次修正步 驟》先前的方法,在此次修正中,於主格點速度與壓力的關係如下式:.

( )

*** * * * * * p P P P P P

P

A

V

H

V

⎟⎟

⎜⎜

⎛ ∀

Δ

=

(4.15) 因此,我們將(4.15)式與(4.7)式相減後,可得速度修正量與壓力修正量的關係 式如下:

( )

'' '' ' ' P P P p P

P

a

V

H

V

⎟⎟

⎜⎜

⎛ ∀

=

(4.16) 接下來,在面上的流量即可表示為: ' ' * * * * * f f f m m

m& = & + & (4.17)

其中m&f**在《第 1 次修正步驟》中已求得。而流量修正量 ' ' f m& 可表示為:

(

)

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

⋅ ' ' '' '' ' ' '' ' ' ρ ρ ρ ρ ρ & (4.18) 將 (4.18)式代入(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 d P A m m ⎟⎟ ∇ ⋅ − + ⋅ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ⋅ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − =

' ' ' ' ' ' * * ' ' ρ ρ ρ & & (4.19) 同等於先前的方法,其必需滿足連續方程式,因此可表示成: 0 ' ' * * *

=

= f f f f m m& & (4.20)

(35)

最後,將(4.19)式代入(4.20)式中,方程式可表示成下線性代數式:

(

)

f p f C C f f f f p f f C C p p S A V A d S P A P A P A ⎟ ∇ ⋅ − − ⋅ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∀Δ + =

' ' '' '' ''

ρ

ρ

(4.21)

=

f c P

A

A

(4.22) 其中 f f f p f c

S

S

A

a

⎛ ∀

Δ

=

δ

ρ

2 4.2.4 壓力修正方程式的整理 在前兩節裡,在兩次的修正步驟中分別導出了兩個壓力修正方程式,我們可 將其整理成一個通式,如下表所示(表 4-1): 2 1 ' ' p p f C C P P

P

A

P

S

S

A

=

+

+

《修正步驟1》 《修正步驟2》 1 p S

f f m& * f P f C C f S A V A

' ' ρ 2 p S

(

)

d S P Ap f f f ⎟⎟∇ ⋅ − ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∀Δ ' ρ P

(

S d

)

Ap f f f ⎟⎟∇ ⋅ − ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∀Δ '' ρ 表 4.1 兩次壓力修正方程式的源項 其中Sp2為非正交的源項,對於正交的網格中,Sp2為 0。 在求解過程,為了使結果更為可靠,因此對於非結構網格系統中,於每個修 正步驟(Correct Step)中,在求解壓力修正方程式的過程裡,我們再將其分成兩部 份運算:第一部份只源項僅包含Sp1,而第二部份則源項只包含Sp2 „ 正交修正(non-orthogonal corrector) : ' ' p1 f c c p pP A P S A =

+ „ 非正交修正(non-orthogonal corrector) : ' ' 2 p f c c p pP A P S A =

+ 對於幾何外形較為複雜的網格,可重覆比較多次的非正交修正。 4.3 進出口邊界條件的給定

(36)

4.3.1 出口壓力邊界的流量計算 在本文研究中出口以壓力邊界為出口條件,其給定方式為在出口面設定一壓 力值,利用此壓力以修正出口部份的流量。定壓力邊界條件的作法如下: 1. 首先,先在壓力邊界相鄰的網格壓力值設定為欲給定的壓力值。 2. 在求解兩次壓力修正方程式之時,需設法將邊界相鄰網格的壓力修正值 設為 0,如圖 4.1 所示,必需在求解矩陣的之前將矩陣調整以滿足 ' =0 p 。 因此,需將壓力修正方程式的係數和源項做調整,即: 2 1 ' ' p p f C C P P

P

A

P

S

S

A

=

+

+

Sp1 =Sp2 =0 , AC =0 3. 在求解完壓力修正之後,需讓邊界相鄰的網格滿足質量守恆定理:

∫∫

⋅ =

=0 f f S m S d V & ρ 如圖 4.2 所示,便可透過上式求得邊界上流量為: m&f,4 =−

(

m&f,1+m&f,2 +m&f,3

)

4.3.2 出口壓力邊界的速度計算 令外,在壓力邊界出口上的速度,則是透過對流邊界條件(Convective Boundary Condition)得來【20】,其數學式子為: 0 = ∂ ∂ + ∂ ∂ z u t c

φ

φ

(4.23) 經離散後可得下列式: 0 = Δ − + Δ − z u t o C n B c o B n B

φ

φ

φ

φ

(4.24) 上式中的上標n,o分別表示新與舊的值,而 下標B,C則代表的是邊界與其相鄰的 網格,故可將方程式整理成下式: Cr Cr cn o B o B + + = 1

φ

φ

φ

(4.25) 其中的 Cr 為 Courant Number,定義為

z

t

u

Cr

c

Δ

Δ

=

(37)

接著,透過前一節計算得到的壓力邊界上的流量來修正邊界上的速度:

(

)

(

u i w k

)

S w S v S u m V b b z b y b x b f b = + + +vb j+ 4 ,

ρ

& (4.26) 4.4 固定牆邊界條件

在牆的部份,其物理條件為無滑移邊界條件(No Slip condition)。為了計算在 邊界旁邊網格的剪應力,我們必需先計算出平行於壁面上的速度V '' 。但在此首 先,我們需先找出平行壁面的向量S''和相對速度V 。(圖 4.3) (1) 平行壁面的向量 '' S : 我們在計算時,垂直壁面的法向量我們定義為: j Sy x + =S i S 。因此平行壁面向量可表示為 '' = y +Sx j i S S (2) 相對速度V : 相鄰邊界網格其速度向量和邊界上速度向量的差值,即

(

u u

) (

i v v

)

j V V V = pb = pb + pb =(Δ iu) +(Δv) j 其中 0 , : : = b b p V V V 在靜止的牆則 壁面速度 相鄰壁面網格的速度 而平行於壁面的速度V ''即為VS''上的投影,其可表示如下: S S S V ' ' 2 '' '' '' =V (4.27) 透過上式,我們可將其整理為如下式:

V ''=(ΔuSySy−ΔvSxSy)i+(ΔvSxSx−ΔuSxSy)j (4.28) 4.5 PISO 演算流程 對於本文中所採用的整個 PISO 演算流程,其表述如下: 步驟 1:给定初始邊界值。 步驟 2:修正進口邊界上速度、流量。 步驟 3:求解動量方程式,得到速度V*。 步驟 4:求解第一次壓力修正方程以得到p ,並進而求出第一次修正過後的 '

(38)

速度、壓力和流量 • * * * * * * , ,P mf V 。 步驟 5:求解第一次壓力修正方程以得到p ,並進而求出第一次修正過後的 '' 速度、壓力和流量 • * * * * * * * * * , ,P mf V 。 步驟 6:計算壓力邊界上的流量。 步驟 7:上述的步驟 2 至步驟 6 為一個完整的 Time step 的流程,因此重複步 驟 2 至 6 直到達到所求解的時間,最後將結果輸出,整個 PISO 流程在圖 4.4 中 描述。

(39)

第五章、無閥型微幫浦簡易理論

本章的內容主要是簡介無閥型微幫浦理論。在 5.1 節我們先簡介關於壓電薄膜位 移的計算方式。在 5.2 節裡我們提出了三種方式針對本文所模擬的無閥式微幫浦 的效率估算方式。在 5.3 節裡則是探討 Nozzle/Diffuser 元件的壓力係數估算方式。 最後,關於單腔微幫浦的數學理論則在 5.4 節裡描述。 5.1 壓電式幫浦的薄膜振動模擬

Fan【8】將流體與結構力學整合 fluid-membrane Coupling 並提出了解析解, 其研究中指出壓電型微幫浦的壓電薄膜振動的數學模式是一個複雜的偏微分方 程式,其包含了腔體壓力、壓電材料造成的週期性作用力和腔體(Diaphragm)本身 的振動(圖 5.1),如下所示: P e f t W h W D = − ∂ ∂ + ∇ 2 2 4

ρ

(5.1) 其中 ) 1 ( 12 2 3 v Eh D − = 為腔體(Diaphragm)的材料係數。 E Diaphragm 基材的楊式係數(Elastic modulus) V Diaphragm 壓電片基材浦松比(Poisson ratio) H Diaphragm 壓電片基材的厚度(Thickness of membrane)

ρ Diaphragm 壓電片基材的密度(Density of membrane)

而等號右手邊的 f 為一週期性的作用力,其計算需透過壓電材料的壓電係e 數和電場以求之。而右手邊第二項的 P 為壓力,壓力的計算必須在每個時間點 (Time Step)下,透過求解流場得到。但求解(5.1)式太過複雜,因此一般在模擬中, 對於壓電薄膜的位移近似提出了不同的方式: ¾ fluid-membrane Coupling Fan【12】等人以 CFD-ACE 模擬方形壓電型微幫浦,將(5.1)式簡化為一維的 問題,將電場、壓電片的振動和流體力學做整合,所使用數學模式包含了壓電薄 膜振動方程式、Navier-Stokes 方程式和連續方程式,並依此求解計算出了在不同 頻率下的薄膜振動曲線。僅管如此,在求解圓形壓電薄膜的位移仍是一個難題。

¾ FSI Module(Fluid-Thermo-Structure Module)

FSI 模組是一個應用來處理液固界面的問題,1997 年 Athavale【9】以及 2007 年 Jeong【15】均曾使用 CFD-ACE 的 FSI 模組模擬 Stemme 型的微幫浦。FSI 的原 理是求解完流場之後,接著求解固體振動方程式並算出位移量,接著重建網格,

參考文獻

相關文件

A factorization method for reconstructing an impenetrable obstacle in a homogeneous medium (Helmholtz equation) using the spectral data of the far-eld operator was developed

(In Section 7.5 we will be able to use Newton's Law of Cooling to find an equation for T as a function of time.) By measuring the slope of the tangent, estimate the rate of change

We would like to point out that unlike the pure potential case considered in [RW19], here, in order to guarantee the bulk decay of ˜u, we also need the boundary decay of ∇u due to

Al atoms are larger than N atoms because as you trace the path between N and Al on the periodic table, you move down a column (atomic size increases) and then to the left across

Due to the important role played by σ -term in the analysis of second-order cone, before developing the sufficient conditions for the existence of local saddle points, we shall

To complete the “plumbing” of associating our vertex data with variables in our shader programs, you need to tell WebGL where in our buffer object to find the vertex data, and

For the data sets used in this thesis we find that F-score performs well when the number of features is large, and for small data the two methods using the gradient of the

In order to identify the best nanoparticle synthesis method, we compared the UV-vis spectroscopy spectrums of silver nanoparticles synthesized in four different green