國 立 交 通 大 學
機 械 工 程 學 系
碩 士 論 文
雙壓電致動器幫浦之流固耦合分析
Analysis of the Fluid-Structure Coupling Flow in A
Bimorph Piezoelectric Micropump
研 究 生：陳 信 宏
指 導 教 授：崔 燕 勇 博士
雙壓電致動器幫浦之流固耦合分析
Analysis of the Fluid-Structure Coupling Flow in A Bimorph
Piezoelectric Micropump
研 究 生：陳信宏 Student：Shin-Hung Chen 指導教授：崔燕勇 Advisor：Yeng-Yung Tsui 國 立 交 通 大 學 機械工程學系 碩士論文 A ThesisSubmitted to Department 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 2010
Hsinchu, Taiwan, Republic of China
雙壓電致動器幫浦之流固耦合分析
研究生：陳信宏 指導教授：崔燕勇 博士 國立交通大學機械工程學系(研究所)碩士班 摘要 本文為研究雙壓電致動器幫浦之流固耦合分析，將 ANSYS 軟體與 CFD 程式 結合來分析此幫浦特性，利用 ANSYS 來計算固體力學部份，而 CFD 程式則負責 計算流體力學之部分。在耦合計算時，將幫浦內的流體壓力和流體剪力視為固體 力學之外力限制條件；將壓電振動平板之位移視為流體之牆邊界條件，來完成 ANSYS 軟體和 CFD 程式的結合。結果發現幫浦之流量同時受到電壓、擺動振幅 和振動頻率影響，且觀察出欲增加幫浦流量的最好之方法為增加電壓的大小。同 時發現將壓電振動平板放置於流體中並且受電壓之負載時，其共振頻率將會有明 顯的下降。 本文同時利用 ANSYS 之調和分析來計算此壓電致動器幫浦在不同外部負載 頻率下的振幅大小，並與 ANSYS 結合 CFD 的結果做比較，得到非常接近的結果， 因此可證明 ANSYS 與 CFD 程式的結合計算方法具有相當程度的可信度。而本程 式之優點為耦合時可以分析流場情況，例如：速度場、流線、壓力場，但 ANSYS 之調和分析無法分析流場的狀況，只能探討壓電振動平板受電壓後於水中之振幅 反應。Analysis of the Fluid-Structure Coupling Flow in A Bimorph Piezoelectric
Micropump
Student：Shin-Hung Chen Advisor：Dr. Yeng-Yung Tsui
Department of Mechanical Engineering National Chiao Tung University
ABSTRACT
This thesis is mainly aimed the calculation of the flow in oscillatory flow pumps using bimorph piezoelectric actuators, combining ANSYS finite element software with CFD program in order to solve Fluid-Structure interaction(FSI) problems. For computational simulations of the piezoelectric pump, the ANSYS is used; for computational simulations of the flow field, the CFD program is used. The result showed that flow rate of the flow pump is deeply affected by voltage, amplitude and frequency. In order to increase the flow rate of flow pumps is applied more high voltage. The resonance frequency of piezoelectric actuators is decreasing clear by setting piezoelectric actuators in fluids.
The amplitude of the piezoelectric actuator of this thesis had been compared with the calculation of ANSYS harmonic analysis, and we received good result. This coupled method can analysis the vector field, streamline, pressure of flow field, but ANSYS harmonic analysis can’t.
誌謝
能完成這篇論文必須要感謝崔燕勇教授，感謝崔老師兩年來的指導，從老師 的身上不只學習到學術知識，學到更多的是做事應有的態度以及分析和解決問題 的能力，還有增加了我的抗壓性，再次感謝崔老師，您的教導我相信一生都受用 無窮！同時也要感謝尹慶中老師和他的學生鄭博毅同學在壓電材料和 ANSYS 軟 體應用上給予我的幫助。 再來是感謝 307 實驗室全體，感謝吳添成學長總是能一通電話過去馬上就給 我協助。感謝胡育昌學長總是細心且熱情的幫我修正 Fortran 程式。感謝林仕文學 長從我一進實驗室那天起就對我特別的照顧，一起修課一同小酌、幫我解決困難， 課業和玩樂方面都照顧的非常周到。感謝陳俊佑學長，幫我解決所有研究上遇到 的問題，並且聽我吐苦水，充當我的心靈導師。感謝我的同窗郭大慶、鄭為陽和 何光桓，尤其是郭大慶你絕對是個好戰友。感謝學弟們林子翔、賴胤男、黃義政， 有你們在的實驗室變的熱鬧多了，招來兩個痞子房業勳和黃思齊，實驗室變的輕 鬆、壓力沒那麼大，希望未來實驗室有你們可以更蒸蒸日上。 接下來感謝幾個朋友，王派宣、王澤瑋、田希文、李其霖、吳建賢、周佳仁、 張國維、黃志瑋、鍾鼐駒、戴培祐、蘇彥甫，感謝你們在我碩士班生活的兩年， 甚至是我求學一路走來，陪我度過所有的喜怒哀樂，沒有你們我走不到這裡，謝 謝，未來還有更長的路，先謝過了兄弟們。 最後感謝我的家人，阿公、阿婆、外公、外婆、爸、媽、姊、兵叔、岱芸、 官俊，家人是最強大的後盾，謝謝你們對我的栽培和放縱，還有不知道歸類何處 的謝孟穎，謝啦。目錄
摘要 ... i ABSTRACT ... ii 誌謝 ... iii 目錄 ... iv 表目錄 ... viii 圖目錄 ... ix 符號表 ... xiv 一、緒論 ...1 1.1 前言(微幫浦介紹)...1 1.2 流固耦合介紹...2 1.3 文獻回顧...3 1.3.1 數值模擬相關的文獻 ...3 1.3.2 類魚幫浦相關的實驗文獻 ...5 1.4 研究內容...6 二、數學模式 ...7 2.1 基本假設...7 2.1.1 流體假設 ...7 2.1.2 固體假設 ...7 2.2 流體運動模式...7 2.2.1 ALE 座標系之統御方程式 ...7 2.2.2 流場邊界條件(Boundary Condition)...8 2.3 固體運動模式...8 2.3.1 結構運動方程式 ...8 2.3.2 壓電方程式 ...92.3.3 壓電參數矩陣極化方向轉換 ...9
2.3.4 振動平板的限制條件(Constraint)...11
三、數值方法 ...12
3.1 流體力學數值方法 ...12
3.1.1 有限體積法(Finite Volume Method) ...12
3.1.1.1 有限體積法的離散...12 3.1.1.2 離散化...12 3.1.1.2.1 非穩態項(Unsteady Term) ...12 3.1.1.2.2 對流項(Convection Term)...13 3.1.1.2.3 擴散項(Diffusion Term) ...14 3.1.1.2.4 源項(Source Term) ...15 3.1.1.3 線性代數式的整理...15
3.1.2 空間守恆定理(Space Conservation Law) ...15
3.1.2.1 空間守恆定理的簡介...15 3.1.2.2 質量源(Mass Sources)...16 3.1.2.3 空間守恆定理...17 3.1.2.4 整理 SCL 應用於 FVM 及 PISO ...18 3.1.3 PISO 演算法 ...18 3.1.3.1 PISO 演算法的簡介 ...18 3.1.3.2 速度與壓力之間的耦合...18 3.1.3.2.1 預測步驟(Predictor Step)...18 3.1.3.2.2 第一次修正步驟(1st Corrector Step) ...19 3.1.3.2.3 第二次修正步驟(2nd Corrector Step) ...21 3.1.3.2.4 壓力修正方程式的整理 ...22 3.1.3.3 邊界條件的給定...23 3.1.3.3.1 進出口壓力邊界的流量計算 ...23
3.1.3.3.2 進出口壓力邊界的速度計算 ...23
3.1.3.3.3 牆邊界條件的計算 ...24
3.2 固體力學數值方法 ...24
3.2.1 有限元素法(Finite Element Method)軟體基本架構 ...24
3.2.2 ANSYS 商用軟體 ...25 3.2.3 限制條件的給定 ...25 3.2.3.1 壁面剪力(Shear Force)...25 3.2.3.2 壁面壓力...26 3.3 網格的移動...27 3.4 流固耦合演算流程 ...27 四、結果與討論...28 4.1 ANSYS 分析...28 4.2 模態分析...28 4.2.1 ANSYS 建模 ...28 4.2.2 自然頻率 ...28 4.3 系統阻尼設定...29 4.4 調和分析...29 4.4.1 單純壓電振動平板的 ANSYS 建模...30 4.4.2 水中壓電振動平板的 ANSYS 建模...30 4.4.3 共振頻率(resonance frequency)和振幅 ...30 4.5 暫態分析...31 4.5.1 ANSYS 壓電振動平板建模 ...32 4.5.2 CFD 流場建模 ...32 4.5.3 流場分析 ...32 4.5.4 電壓對流量造成之影響 ...34 4.5.5 增加流道高度對流場之影響 ...34
五、結論 ...36 參考文獻 ...37 附錄 ...91
表目錄
表(一) PZT-4 的材料常數【17】 ...39 表(二) 黃銅的材料常數 ...39 表(三) 模態分析結果之前 10 組自然頻率 ...40 表(四) 阻尼比係數 ...40 表(五) 水的性質 ...40 表(六) 電壓 30V、60V、120V 在不同外部負載頻率下之流量結果 ...41 表(七) ANSYS 調和分析和 ANSYS 結合 CFD 結果比較，流道高度為 5mm 和 8mm 的 振動平板振幅 ...41圖目錄
《圖 2. 1》雙壓電致動器(a)擺動示意圖(b)平板材料剖析...42 《圖 2. 2》流場尺寸示意圖...42 《圖 2. 3》座標軸方向代號...43 《圖 2. 4》壓電致動器之電壓給定致動方式(parallel mode) ...43 《圖 3. 1》擴散項離散化 over-relaxed 法之標示圖...44 《圖 3. 2》計算邊界壓力示意圖...44 《圖 3. 3》利用網格速度計算體積改變量...45 《圖 3. 4》(a)顯式法和(b)隱式法的體積改變結果...45 《圖 3. 5》不同時階大小的隱式法結果(a)較大時階(b)較小時階...46 《圖 3. 6》利用各面移動前後的新舊點位置直接算出體積...46 《圖 3. 7》兩個時階中的網格位置變化(a)各面掃過的體積(b)顯示體積重疊 ...47 《圖 3. 8》壓力邊界示意圖...48 《圖 3. 9》壓力邊界之流量計算...48 《圖 3. 10》計算壁面剪力示意圖...49 《圖 3. 11》流場網格分割區塊 ...49 《圖 3. 12》一週期內的流場局部網格運動情形(a)1/4T(b)2/4T(c)3/4T(d)1T ...51 《圖 3. 13》流固耦合演算流程...52 《圖 4. 1》模態分析、調和分析之單純壓電振動平板、暫態分析的模型建立...53 《圖 4. 2》模態分析限制條件示意圖...53 《圖 4. 3》模態分析之模型完成圖...54 《圖 4. 4》第一模態之模態圖...54 《圖 4. 5》第二模態之模態圖...54 《圖 4. 6》第三模態之模態圖...55 《圖 4. 7》調和分析之單純壓電振動平板的限制條件示意圖...55《圖 4. 8》調和分析之單純壓電振動平板的模型完成圖...56 《圖 4. 9》調和分析之水中壓電振動平板的模型建立...56 《圖 4. 10》調和分析之水中壓電振動平板的限制條件示意圖...57 《圖 4. 11》調和分析之水中壓電振動平板的模型完成圖...57 《圖 4. 12》阻尼比為 0.01 時，ANSYS 調和分析在無水及加入有水元素(Fluid29)之比 較，振動平板尾端振幅隨頻率的變化...58 《圖 4. 13》阻尼比為 0.1 時，ANSYS 調和分析在無水及加入有水元素(Fluid29)之比較， 振動平板尾端振幅隨頻率的變化...58 《圖 4. 14》阻尼比為 1 時，ANSYS 調和分析在無水及加入有水元素(Fluid29)之比較， 振動平板尾端振幅隨頻率的變化...59 《圖 4. 15》阻尼比為 15 時，ANSYS 調和分析在無水及加入有水元素(Fluid29)之比較， 振動平板尾端振幅隨頻率的變化...59 《圖 4. 16》ANSYS 調和分析在無水及加入有水元素(Fluid29)之不同阻尼比振動平板尾 端振幅隨頻率的變化...60 《圖 4. 17》ANSYS 調和分析在無水及加入有水元素(Fluid29)之不同阻尼比振動平板尾 端振幅隨頻率的變化...60 《圖 4. 18》暫態分析的限制條件示意圖...61 《圖 4. 19》暫態分析之流體壓力給定...61 《圖 4. 20》暫態分析之流體剪力給定...62 《圖 4. 21》暫態分析之固力模型完成圖...62 《圖 4. 22》流場網格數目和尺寸圖(a)流場實際網格尺寸(b)網格 6976 格的局部放大圖 (c) 網格 13952 格的局部放大圖(d)網格 27904 格的局部放大圖...63 《圖 4. 23》三種網格測試的流量結果圖，阻尼比 15、電壓 120V、外部負載頻率 20hz ...63 《圖 4. 24》頻率 10，最後一週期振動平板尾端頂點的變化圖(a)位置隨時間的變化(b) 速度隨時間的變化...64
《圖 4. 25》頻率 10，振動平板振幅放大 50 倍運動圖(t=16/100T、41/100T、66/100T、 91/100T) ...65 《圖 4. 26》頻率 10，流場壓力隨時間變化圖(t=16/100T、41/100T、66/100T、91/100T) ...66 《圖 4. 27》頻率 10，流場流線隨時間變化圖(t=16/100T、41/100T、66/100T、91/100T) ...67 《圖 4. 28》頻率 20，最後一週期振動平板尾端頂點的變化圖(a)位置隨時間的變化(b) 速度隨時間的變化...68 《圖 4. 29》頻率 20，振動平板振幅放大 50 倍運動圖(t=16/100T、41/100T、66/100T、 91/100T) ...69 《圖 4. 30》頻率 20，流場壓力隨時間變化圖(t=18/100T、44/100T、68/100T、94/100T) ...70 《圖 4. 31》頻率 20，流場流線隨時間變化圖(t=18/100T、44/100T、68/100T、94/100T) ...71 《圖 4. 32》頻率 50，最後一週期振動平板尾端頂點的變化圖(a)位置隨時間的變化(b) 速度隨時間的變化...72 《圖 4. 33》頻率 50，振動平板振幅放大 50 倍運動圖(t=16/100T、41/100T、66/100T、 91/100T) ...73 《圖 4. 34》頻率 50，流場壓力隨時間變化圖(t=20/100T、45/100T、70/100T、95/100T) ...74 《圖 4. 35》頻率 50，流場流線隨時間變化圖(t=20/100T、45/100T、70/100T、95/100T) ...75 《圖 4. 36》頻率 100，最後一週期振動平板尾端頂點的變化圖(a)位置隨時間的變化(b) 速度隨時間的變化...76 《圖 4. 37》頻率 100，振動平板振幅放大 50 倍運動圖(t=16/100T、41/100T、66/100T、
《圖 4. 38》頻率 100，流場壓力隨時間變化圖(t=30/100T、50/100T、80/100T、1T) 78 《圖 4. 39》頻率 100，流場流線隨時間變化圖(t=30/100T、50/100T、80/100T、1T) 79 《圖 4. 40》阻尼比為 15 時，不同頻率於 0~5 秒的流量隨時間之變化 ...80 《圖 4. 41》ANSYS 調和分析和 CFD-ANSYS 耦合分析之結果：阻尼比為 15 時，振動 平板尾端的振幅大小...80 《圖 4. 42》頻率 10、20、50、100 在各個頻率的特定取樣時間的振動平板振幅放大 50 倍運動圖 ...81 《圖 4. 43》ANSYS 調和分析和 CFD-ANSYS 耦合分析之結果：阻尼比為 15 時，電壓 為 30V、60V、120V 振動平板尾端的振幅大小 ...82 《圖 4. 44》阻尼比為 15 時，電壓為 30V、60V、120V 時 CFD-ANSYS 耦合分析之流 量大小 ...82 《圖 4. 45》頻率 10、入口高度 5mm 和 8mm，最後一週期振動平板尾端頂點的變化圖 (a)位置隨時間的變化(b)速度隨時間的變化...83 《圖 4. 46》頻率 10，入口高度 8mm，流場壓力隨時間變化圖(t=16/100T、41/100T、 66/100T、91/100T) ...84 《圖 4. 47》頻率 10，入口高度 8mm，流場流線隨時間變化圖(t=16/100T、41/100T、 66/100T、91/100T) ...85 《圖 4. 48》頻率 20，入口高度 8mm，流場流線隨時間變化圖(t=18/100T、44/100T、 68/100T、94/100T) ...86 《圖 4. 49》頻率 50，入口高度 8mm，流場流線隨時間變化圖(t=20/100T、45/100T、 70/100T、95/100T) ...87 《圖 4. 50》頻率 100，入口高度 8mm，流場流線隨時間變化圖(t=30/100T、50/100T、 80/100T、1T)...88 《圖 4. 51》阻尼比為 15 時，入口高度 8mm，不同頻率於 0~5 秒的流量隨時間變化 ...89 《圖 4. 52》阻尼比為 15 時，入口高度為 5mm 和 8mm 的 CFD-ANSYS 耦合分析之流
量大小 ...89 《圖 4. 53》ANSYS 調和分析和 CFD-ANSYS 耦合分析之結果：阻尼比為 15 時，入口 高度為 5mm 和 8mm，振動平板尾端的振幅大小 ...90
符號表
符號 定義 c 剛性常數 C 阻尼係數 d 壓電應變常數 dK 正交向量 D 電位移 e 壓電應力常數 E 電場 C f F 面上的對流通量 D f F 面上的擴散通量 K 勁度(stiffness) m 質量流率 f m 面上的質量流率 r f m 流體相對於面上的質量流率 M 質量 P 流體壓力 S 應變 f SK 面的法向量 w SK 壁面之面法向量 t 時間 T 應力 VK 流體速度向量 g VK 網格速度向量 希臘符號 α 庫倫阻尼係數 β 結構阻尼係數 Pnb δK P 到 nb 距離向量 ε 介電常數 ζ 阻尼比 μ 黏滯係數 ρ 密度 w τK 壁面剪應力 ω 自然角頻率f Δ∀ 一個面掃過的體積 上標 n 此一時間的未知數 o 前一時間的已知數 ─ 主格點與鄰近格點之內差 ′ 修正量 下標 f 面上的值 nb 鄰近格點的值 P 主格點的值 w 壁面上的值
一、緒論
1.1 前言(微幫浦介紹) 利用微尺寸幫浦來傳輸液體或是氣體已被廣泛的研究，而微幫浦發展至今大 概分為兩大類：機械式的位移幫浦(displacement pumps)和非機械式的動力幫浦 (dynamics pumps)。位移幫浦主要有三種形式：蠕動式(peristaltic)、轉動式(rotary)、 往 復 式 (reciprocating) ； 其 動 力 來 源 包 含 ： 靜 電 式 (Electrostatic) 、 熱 氣 式 (Thermopneumatic)、壓電式(Piezoelectric)和電磁式(Electromagnetic)。其中往復式 幫浦是利用機械式的振動方式來推動流體，最著名的為薄膜幫浦(diaphragm pump) 。 非 機 械 式 的 動 力 幫 浦 則 是 利 用 離 心 力 (centrifugal) 、 電 水 力 學 (electrohydrodynamic)、聲學(acoustic streaming)和超音速(ultrasonic)等原理所運作 的幫浦。位移幫浦適合用來傳遞似氣體的流體和低雷諾數流體，而動力幫浦主要 用來分析流體的一些特性，例如流體的離子強度。 在眾多類型的微幫浦中，唯讀壓電微幫浦的研究不斷的推陳出新。結合壓電 致動器使微幫浦的型態更進一步，因壓電微幫浦的壓電致動器具有簡易的結構、 高輸出效率、高精確度和低耗電功率等優勢。因此可應用範圍極廣，例如生醫工 程裡的血液幫浦、微量皮下葯劑輸送系統：在傳遞藥物時是因為處在高壓環境， 進而產生的紊流和過大的剪應力會破壞生物組織，造成生物細胞的壞死，若經由 控制信號與微幫浦振動頻率，將藥物定時定量的輸入到人體內的循環系統中，如 此一來便不需要經常的去注射藥物，也可以減輕病患的痛苦。或是電子散熱元件： 電腦主機的冷卻風扇常見為旋轉葉片型，而此類型風扇通常伴隨著高電功率。雖 然在許多應用的成果令人滿意，但旋轉馬達所產生的噪音常常令人無法忍受。因 此發展出雙壓電致動型風扇，雙壓電致動型風扇的結構簡單且可精確的控制風 速。我們將旋轉葉片型風扇和雙壓電致動型風扇比較後可得知，雙壓電致動型風 扇沒有嘈雜的噪音、具有較長的壽命期(約 25000 小時)和低耗電功率(小於 100mW)。因此可以推斷雙壓電致動型風扇在許多應用層面上，將可取代傳統的旋轉型風扇。除此之外壓電致動器也可利用它的推進特性，應用於水中的推進裝置。 此種推進裝置是將兩片壓電片夾住一平板使它擺動，類似魚尾的擺動。這種類生 物的推進裝置具有高推進效率和無聲的特性，因此可運用於軍事單位上。本文即 是研究此類生物推進裝置，將其固定於水中，運用其推進原理成為新型的幫浦。 1.2 流固耦合介紹 流固耦合可應用於分析飛機的機翼振翅、大樓受風力的影響、流體流經彈性 管和血管問題。本文即是研究壓電振動平板和流體的交互作用，利用振動平板的 擺動來推動流體，然而振動平板在推動流體的同時，流體也會產生一反作用力去 影響振動平板的擺動，如此交互的作用稱為流固耦合。在解決耦合問題時，必須 同時解決流體力學和固體力學問題，而這兩種問題通常用不同的方法來求解，因 此兩種物理現象用一套程式來求解將是一個困難的挑戰。 工 業 應 用 和 科 學 研 究 的 許 多 問 題 中 ， 都 會 產 生 流 固 耦 合 (Fluid-Structure interaction，FSI)的物理現象。耦合的問題牽涉到多方面的物理現象，因此學者們 發展出一些方法來近似流固耦合。解決偶合問題共可分為完全統一法(monolithic approaches)和分段法(partitioned approaches)【3】。完全統一法在解流體問題和固體 問題是用同一套程式、且同時進行求解。分段法在解流體之統御方程式和解固體 之結構特性時需用不同的程式。在分段法中，通常以有限元素法(finite element method，FEM)來求解固體力學的問題；解決流體力學方面的問題時，較常見為利 用有限體積法(finite volume method，FVM)和有限元素法。而利用有限體積法來求 解固體力學問題也是可行，但較少見。分段法中另一項值得注意的是耦合的強度。 在每個時階(time step)中計算流體與固體的反覆交互作用時，決定耦合成功與否的 收斂條件稱作隱式分段法(implicit partitioned approach)。在一個時階中，隱式分段 法的耦合強度優於顯式法(explicit approach)。
解決 FSI 的流體問題時，目前已證實以 ALE(arbitrary Lagrangian-Eulerian)座標 系之統御方程式配合移動網格的方法來求解最合適。
1.3 文獻回顧 1.3.1 數值模擬相關的文獻 (A) 移動網格的研究 本文研究的流場中，壓電振動平板做大幅度的振動，流場中的網格會隨著振 動平板一起運動，因此我們必須探討移動邊界和移動網格的問題和方法。以下為 學者們對移動邊界和移動網格所做的相關研究。 1988 年，Demirdzic et al.【1】 研究移動邊界問題中，由質量守恆定理推導出 空間守恆定理，並將空間守恆定理結合有限體積法來求解問題。並解決空間守恆 可能衍生出人為的質量源誤差問題。 1998 年，Egelja et al.【2】 利用有限體積法求解自由表面流，邊界利用非正 交的固定邊界和移動網格來近似自由表面。運用 SIMPLE 壓力修正法結合移動網 格的空間守恆定理進行疊代，使流場滿足流體的連續性。
2008 年，Sternel et al.【3】 利用隱式分割(implicit partitioned)的 ALE 法來近 似流固耦合(FSI)的問題。並且深入研究求解過程中運用非線性多重網格法
(nonlinear multigrid)、適當的鬆弛因子(underrelaxation factor)和適當的移動網格對 耦合問題的準確性提升效果。最後得到線性配合橢圓的混合網格法和純橢圓網格 法，此二種網格法比較可以節省計算時間和提供較好的準確度。
2008 年，Zhang et al.【4】 研究包含移動邊界的三維非穩態不可壓縮流。在 計算移動網格時，空間差分利用有限體積法，時間利用 implicit dual time stepping scheme，則可滿足幾何守恆定律。選擇三角網格來建立類魚模型，給定魚尾的運 動方式，進行數值模擬。最後得到的數據和文獻比較，顯現出相同的趨勢。 2009 年，Srikanth et al.【5】 利用 CFD 來解決可壓縮的空氣流入典型吸入式 腔體的問題，並且研究移動網格和固定電極連結的方法。利用 Altair’s HyperMesh 軟體建立腔體模型和結構性網格。利用 ANSYS-CFX 給定邊界條件且進行流場的 模擬計算。結果顯示腔體移動的速度，確切影響吸入式腔體的壓力變化。並且得 到吸入式腔體的運動為定性的，可作為研究 electro-fluid dynamics 斷路器的基礎。
(B) 流固耦合的研究
2002 年，Slone et al.【6】 致力於只使用一套程式─PHYSICA─來解決流固 耦合的問題，且流場和固力結構體的網格處理也是直接一起分割，因此不會有兩 套軟體整合上的問題和繁雜的互傳資料。流固耦合的方式採用分段法，流力部分 在解 N-S 方程式是利用 FVM 來離散，利用 SIMPE 法則來作壓力修正。固力部分 是利用 FEM 來離散方程式，由 Newmark 預測修正演算法來解方程式。他們所分 析的流固耦合例子，為三維流場中置入一長方形棒，由流體對長棒作用，計算出 棒子受力後的彎曲變形和流場結果。
2006 年，Nakasone et al.【7】 利用 ANSYS 求解雙壓電晶片致動器和流體之 間的交互作用以及流量大小，並藉由改變致動器的振動頻率、振幅、平板厚度以 及流體的性質來探討整體的流量。此作者在模擬計算時並非是完全的流固耦合模 擬，是先使用 ANSYS 的調和分析計算流固耦合後壓電振動平板的隨外部負載頻 率改變的振幅，再利用不同外部負載頻率下的振幅大小，輸入到單純 CFD 程式進 行流場計算，最後經由實驗結果證實與模擬結果相似。
2007 年，Vatanabe et al.【8】 利用 ANSYS 求解雙壓電晶片雙致動器和流體 之間的交互作用以及流量大小。在管道中置入兩個壓電致動器，比較不同的擺放 方式對流場的出口流量的影響。此作者同樣在模擬計算時並非是完全的流固耦合 模擬，也是先使用 ANSYS 的調和分析計算流固耦合後壓電振動平板之隨外部負 載頻率改變的振幅頻率之間的關係，再利用不同外部負載頻率下的振幅大小，輸 入到單純 CFD 程式進行流場計算。在實驗和模擬比較數據的部份有較大誤差，但 由模擬分析可看出於管道中置入兩壓電致動器在無相位差時有較大的流量輸出。 2008 年，Xia and Lin【9】 利用非結構性網格和有限體積法來離散固體力學 方程式，並且評估有限體積法適合用在分析受流體影響的固體結構上。時間上採 用 implicit dual-time 的方法能得到較高的準確度，且作者採用之耦合方法為分段 法。最後作者和其他學者所研究的數值計算結果相比較得到相當好的結果。
2009 年，Lima et al.【10】 利用 ANSYS 求解雙壓電晶片致動器和流體之間的 交互作用以及流量大小，並藉由改變致動器的振動頻率、振幅、平板厚度以及流 體的性質來探討整體的流量。此作者在模擬計算時並非是完全的流固耦合模擬， 是先使用 ANSYS 的調和分析計算流固耦合後壓電振動平板的隨外部負載頻率改 變的振幅頻率之間的關係，再利用不同外部負載頻率下的振幅大小，輸入到單純 CFD 程式進行流場計算，求出平板第二模態震動比平板第一模態震動所產生的流 量大。未來並研究壓電電所產生的熱對流場的影響。 2009 年，陳俊佑【11】 研究壓電無閥式微幫浦流固耦合分析，利用 CFD 程 式結合商用程式 ANSYS 進行運算。ANSYS 負責固體力學部份之壓電致動器運 算，並將壓電致動器各位置變形的速度傳給 CFD 運算，CFD 則回傳壓力至 ANSYS 進行流固耦合的運算。並且探討在不同操作頻率下微幫浦之特性。 1.3.2 類魚幫浦相關的實驗文獻 2000 年，Yoo et al.【12】研究壓電型風扇來替代嘈雜的旋轉型葉片風扇。分 別以電壓 110 V、220 V 和頻率 60 Hz 來測試壓電風扇對電子裝置的冷卻效果，也 同時分析和測試不同的振動薄片所造成的影響。作者發現，當固定壓電片的長度 (L)，振動薄片的長度(l)增加反而造成振動薄片的共振頻率降低。最後測試出由 220 V 和 60 Hz 所驅動的連續型葉片(series-type fan)、振動薄片長度 l=31.8 mm、薄片 材質 S6 時，有最大位移 35.5 mm、最大風速 3.1 m/s。且振動薄片的長度影響共振 頻率較寬度來得顯著。 2007 年，Kosa et al.【13】研究一種新穎的推進法，原理是利用擺動的後尾使 之可以在低雷諾數的流場前進。先利用機械、電和流體的耦合模型─放大實驗的 模型─來分析出解析解。由解析分析出在後尾的長度等於 10 mm 時，可以使微小 的機器人以最佳的速度 5 cm/s 推進。由實驗證明此推進法在低雷諾數確實可發揮 效益。
1.4 研究內容 本文的研究，主要是利用計算流體力學和計算力學的方式來分析雙壓電振動 幫浦之流固耦合問題。將壓電材料施以正負交替的電壓驅動平板，使薄板產生週 期性的上下擺動。由於振動平板的擺動類似魚尾運動，我們將振動平板前端固定， 可造成流體向後方流動，達到流體傳輸的功能。 流體的數值模擬中，我們採用非結構性網格的數值方法，流體之統御方程式 採用 ALE 座標系統，以便於解決移動邊界的問題。利用有限體積法(Finite Volume Method)離散統御方程式，採用 PISO 壓力速度修正法和空間守恆定理(Space Conservation Law)來滿足連續方程式，可求得在每個時階流場中的情況。將流體所 計算出流體壓力和流體剪力施加於 ANSYS 中振動平板的邊界，ANSYS 之振動平 板再同時被施加電壓，來計算出下個時階的振動平板位移，觀察反覆的作動情況 下之流場情形。並比較施加不同外部負載電壓之大小和頻率，來分析幫浦的流量 大小，進而求得幫浦的最大流量。 並且探討壓電振動平板跟流體耦合後的振動情形，包含：共振頻率、擺動振 幅。比較不同外部負載電壓和外部負載頻率，壓電致動器之致動情形。
二、數學模式
本文模擬之雙壓電振動幫浦其物理模型共包含流體、壓電材料、金屬材料， 模型示意圖如《圖 2. 1》，利用壓電和金屬材料所組成之壓電振動平板的擺動來推 動流體，且流體將會回饋壓力和剪力給壓電和金屬，因此模擬計算必須同時利用 流體力學和固體力學的數值方法，並且將流體固體耦合。 本章之數學模型共分為流體運動模式和固體運動模式部分，因為兩者的運動 方程式之型式不同，以及不同的邊界條件關係，所以分別描述之。本文固體和流 體詳細的模型尺寸示意圖如《圖 2. 2》所示。 2.1 基本假設 2.1.1 流體假設 本文所假設的工作流體性質如下：不可壓縮流(incompressible Flow)，假設工 作流體的密度為常數。忽略重力項(neglect body force)。非穩態流場(unsteady)，在 本文模擬中，流場中的壓電振動平板會隨時間變化，平板的振動造成兩端進出口 有流體流進或流出，故假設流場為非穩態的情況。 2.1.2 固體假設 二維之振動平板假設為有限厚度之平板，壓電致動器型式採用棒型致動器。 2.2 流體運動模式 2.2.1 ALE 座標系之統御方程式 在計算流體力學問題時若網格為靜止不動，一般是使用 Eulerian 座標系統。 Eulerian 座標系之統御方程式： 連續方程式(Continuity Equation)：( )
V 0 t ρ ρ ∂ + ∇ ⋅ = ∂ K (2.1) 動量方程式(Momentum Equation)：( )
V(
VV)
P(
V)
t ρ ρ μ ∂ + ∇ ⋅ = −∇ + ∇ ⋅ ∇ ∂ K K K K (2.2) 但在探討流固耦合時，是無法單獨使用 Eulerian 座標系統求解的複雜流場，因此本文利用 ALE 座標系統。在 ALE(arbitrary Lagrangian-Eulerian)座標系統下之 網格，不管網格移動的速度分佈成線性或非線性，皆可容易的描述。當網格的移
動速度和流體速度相等時，可視為 Lagrangian 系統；當存在網格固定不動時，則 可視為 Eulerian 系統。所以 ALE 座標系適合用來分析流體的交界面、自由表面流 和最常見的移動邊界問題。本文研究的流場即是必須使用 ALE 座標系求解移動邊 界問題。 ALE 座標系之統御方程式： 連續方程式(Continuity Equation)：
(
V Vg)
0 t ρ ρ ∂ + ∇ ⋅⎡ − ⎤= ⎣ ⎦ ∂ K K (2.3) 動量方程式(Momentum Equation)：( )
(
g)
(
)
V V V V P V t ρ ρ μ ∂ ⎡ ⎤ + ∇ ⋅⎣ − ⎦= −∇ + ∇ ⋅ ∇ ∂ K K K K K (2.4) 其中μ 為黏滯係數、 ρ 為流體密度、∇P為壓力梯度、VK為流體速度、VKg為網格速 度、V −Vg K K 代表流體相對於網格的流速。比較 Eulerian 和 ALE 座標系之統御方程 式可明顯的發現，ALE 座標系之統御方程式只是在連續方程式和動量方程式中的 對流項多了相對速度的概念。 2.2.2 流場邊界條件(Boundary Condition) (A) 出入口邊界條件 本文模擬的進出口邊界為壓力邊界。 壓力邊界：考慮兩端出入口在不同壓力差的流量變化，出入口的速度計算同樣採 用對流邊界條件，其數學式子和詳細計算過程稍後在 3.1.3.3 節介紹。 (B) 牆邊界條件牆為無滑移邊界條件(No-Slip Condition)，牆和流體速度相等(u u− wall =0，
0 wall v v− = )。於本文研究中，在流體部分將壓電振動平板也設為牆邊界條件，因此 壓電振動平板為具有速度之牆邊界，且壓電振動平板運動時將會改變流場之幾何 形狀。 2.3 固體運動模式 2.3.1 結構運動方程式 結構受外力作用，其運動狀態可以用應力型式來表現，稱為 Navier equation。
結合適當的邊界條件和壓電方程式，便可求得棒型壓電致動器之位移、應變和應 力等性質。下式即為直角座標系之 Navier equation： 2 2 ij i i i i j u f u f x t σ ρ ρ ∂ ∂ = + = + ∂ ∂ (2.5) 其中σ 代表應力、ij ui代表位移、ρ 代表結構密度、 fi代表作用在結構上之外力。 2.3.2 壓電方程式 壓電材料是具有壓電效應和焦電效應的一種材料，這些效應構成了物理中的 力學、電學和熱學之間的交互關係。而早在 1880 年居禮兄弟(Piere Curie and Jacques Curie)即發現了壓電效應，至今仍然被廣泛的運用在生活之中，更表現出越簡單的 科技，實是越容易被廣泛的長久使用，【14】。 本文選用之壓電材料必須將電能(電學)轉換機械能(力學)，因此採用機電交互 作用為 e-form 型式的壓電方程式： E p pq q pj j T =c S −e E S i iq q ij j D =e S +ε E (2.6) 其中 E 代表電場、T 代表應力、D 代表電位移、S 代表應變。c 代表剛性常數、e 代表壓電常數、ε 代表介電常數。下標第一個符號代表應力所在平面的法線方向； 下標第二個符號代表應力的作用方向；其中 p,q=1、2、…、6；i,j=1、2、3；且由 《圖 2. 3》可知 1、2、3 各代表 x、y、z 軸方向，4、5、6 各代表 yz、zx、xy 方 向。因此可知(2.5)式和(2.6)式之應力應關係為σ11 =T1、σ22 =T2、σ33 =T3、σ23 =T4、 13 T5 σ = 、σ12 =T6。 2.3.3 壓電參數矩陣極化方向轉換 將 e-form 型式的壓電方程式完全以矩陣展開可得：
1 11 12 13 14 15 16 11 21 31 2 21 22 23 24 25 26 12 22 32 3 31 32 33 34 35 36 13 23 33 4 41 42 43 44 45 46 14 24 34 5 51 52 53 54 55 6 1 2 3 E E E E E E E E E E E E E E E E E E E E E E E E E E E E E T c c c c c c e e e T c c c c c c e e e T c c c c c c e e e T c c c c c c e e e T c c c c c c T D D D ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ 1 2 3 4 5 56 15 25 35 6 61 62 63 64 65 66 16 26 36 1 11 12 13 14 15 16 11 12 13 2 21 22 23 24 25 26 21 22 23 3 31 32 33 34 35 36 31 32 33 E E E E E E E S S S S S S S S S S S S S S e e e S c c c c c c e e e E e e e e e e E e e e e e e E e e e e e e ε ε ε ε ε ε ε ε ε ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎣ ⎦ ⎣ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ (2.7) 當壓電材料的極化方向為 Z 方向時，上方矩陣可表示如下： 1 11 12 13 31 2 12 11 13 31 3 13 13 33 33 4 44 15 5 44 15 6 66 1 15 11 15 11 2 31 31 33 33 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 E E E E E E E E E E E E S S S T c c c e T c c c e T c c c e T c e T c e T c D e e D e e e D ε ε ε ⎡ ⎤ ⎡ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎣ ⎦ 1 2 3 4 5 6 1 2 3 S S S S S S E E E ⎤ ⎡ ⎤ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎦ (2.8) 11 22 66 2 E E E c c c = − (2.9) 欲將 Z 方向極化的壓電材料參數轉換為 X 方向極化和 Y 方向極化，必須做座 標轉換的處理，因此 X 方向極化的壓電方程式表示如下： 1 33 13 13 33 2 13 11 12 31 3 13 12 11 31 4 66 5 44 15 6 44 15 1 33 31 31 33 15 11 2 15 11 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 E E E E E E E E E E E E S S S T c c c e T c c c e T c c c e T c T c e T c e D e e e e D e D ε ε ε ⎡ ⎤ ⎡ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎣ ⎦ 1 2 3 4 5 6 1 2 3 S S S S S S E E E ⎤ ⎡ ⎤ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎦ (2.10)
Y 方向極化的壓電方程式： 1 11 13 12 31 2 13 33 13 33 3 12 13 11 31 4 44 15 5 66 6 44 15 1 15 11 31 33 31 33 2 15 11 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 E E E E E E E E E E E E S S S T c c c e T c c c e T c c c e T c e T c T c e D e e e e D e D ε ε ε ⎡ ⎤ ⎡ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎣ ⎦ 1 2 3 4 5 6 1 2 3 S S S S S S E E E ⎤ ⎡ ⎤ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎦ (2.11) 本文之模型，定義水平方向為 X 方向、垂直方向為 Y 方向，因此選定 Y 方向 極化的壓電方程式，詳細的轉換座標方法可參考附錄 1，而本文之極化情況和致 動方式如《圖 2. 4》。 2.3.4 振動平板的限制條件(Constraint) 我們將壓電振動平板一端固定(clamped)，另一端為自由端，如《圖 2. 1》可做 往復式的擺動，其限制條件數學表示式如下：
( )
0, 0 , x( )
0, 0 u t = u t = (2.12) 其中u代表位移、ux代表應變(strain)。 壓電振動平板與流體接觸的周圍必須受外力的限制條件，外力即流體的壓力 和流體剪力。 壓電振動平板之壓電部分具有電壓的限制條件。 2.4 流固耦合 本文之模擬首先由電壓驅動振動平板移動，此時壓電振動平板的移動造成流 場的幾何形狀變形和流體的流動，流體產生壓力和剪力施加於振動平板上，如此 反覆的運動使流體運動速度逐漸變大並且到達週期性的變化，進而觀察中間轉換 過程和到達週期性變化時的流場情況以及固體結構運動情形。三、數值方法
本文之數值方法同樣也分為流體力學和固體力學部份，其中流體力學數值方 法運用於 CFD 程式共包含有限體積法、空間守恆定理和 PISO 演算法；固體力學 數值方法則是利用有限元素法之套裝商用軟體 ANSYS 來做分析。 3.1 流體力學數值方法 本文計算流體力學部份，是對第二章所提到的 ALE 座標系之統御方程式做處 理，利用有限體積法離散統御方程式，並且運用空間守恆定理和 PISO 演算法來滿 足流體之連續性，並進而求解出流場內各個位置的速度和壓力值。3.1.1 有限體積法(Finite Volume Method) 3.1.1.1 有限體積法的離散 我們將前一章裡所提到的統御方程式(2.4)式改成下列通式：
( )
(
)
(
)
g V V V V V q t ρ ρ μ ∂ ⎡ ⎤ + ∇ ⋅⎣ − ⎦= ∇ ⋅ ∇ + ∂ K K K K K (3.1) 將上式積分後可得下式：(
g)
(
)
Vd V V V d V d qd t ∀ ρ ∀ ρ ∀ μ ∀ ∂ ⎡ ⎤ ∀ + ∇ ⋅⎣ − ⎦ ∀ = ∇ ⋅ ∇ ∀ + ∀ ∂∫∫∫
∫∫∫
∫∫∫
∫∫∫
K K K K K (3.2) 【非穩態項】 【對流項】 【擴散項】 【源項】 透過高斯散度定理，將(3.2)式的對流項和擴散項由體積分轉換為面積分，可得下 式：(
g)
S S Vd V V V dS V dS qd t ∀ ρ ρ μ ∀ ∂ ∀ + − ⋅ = ∇ ⋅ + ∀ ∂∫∫∫
∫∫
∫∫
∫∫∫
K K K K K K Kw
w
(3.3) 3.1.1.2 離散化 3.1.1.2.1 非穩態項(Unsteady Term) 對於有移動邊界的不可壓縮流，其網格的體積會隨時間而改變。( )
V(
n o)
d V V t t ρ ρ ∀ ∂ ⋅ Δ∀ ∀ = − ∂ Δ∫∫∫
K K K (3.4) o VK 為前一個時間的已知數、VKn為這一個時間的未知數，因此在計算上可將VKo移進源項。 3.1.1.2.2 對流項(Convection Term) 將(3.3)式中的對流項離散化，可得下式：
(
)
C g f f f f f S V V V dS m V F ρ − ⋅ =∑
=∑
∫∫
K K K K Kw
(3.5)其中下標 f 代表面上的值、mf代表流經面上的質量流率(mass flow rate)、SKf代表
面之法向量、 C f F 代表面上的對流通量，而對於在面上的值，我們是透過相鄰網格 以線性內差的方式求之。流體相對於面上的質量流率定義如下：
(
)
r f f g f f m =ρ VK−VK ⋅SK (3.6) 上式為 ALE 座標系統御方程式的重點，因加進了網格速度的變化，網格的移動速 度會直接影響面上的質量流率。因此上式的計算方式必須配合空間守恆定理，詳 細過程稍後於 3.1.2.3 節作介紹。 另 外 ， 本 文 求 解 對 流 通 量 採 用 的 方 式 為 一 階 上 風 法 (first-order upwind difference)和中央差分法(second-order central difference)的混合方式，稱作混合差分 法(hybrid difference)，因此對流通量可表示如下：(
)
C UD CD UD f f f f F = F + γ F −F (3.7) 其中γ 的值介於 0~1 之間，當γ =0時為一階上風法、γ =1時為中央差分法。上式 等號右手邊第二項為前一個時間的值，故將其移至源項中計算。 在一階上風法裡，面上的值取決於其上游。例如：mf 由主格點流向外，面上 的值等於主格點的值；mf 由外流向主格點，面上的值等於鄰近格點的值。 在中央差分法裡，Vf K 的值取決於相鄰網格的內差：(
1)
f p nb VK = −w VK +wVK (3.8) 其中w為權重係數(weighting factor)、下標P代表主格點的值、下標nb代表鄰近格 點的值。最後將一階上風法和中央差分法的關係式代入(3.7)式可得：
(
)
(
)
max , 0 max , 0 C n n f f p f nb F = m VK + −m VK +γ{
mf ⎣⎡(
1−w V)
Kpo +wVKnbo⎦⎤−max(
mf, 0)
VKpo −max(
−mf, 0)
VKnbo}
(3.9) 其中上標n代表這個時間的未知數、上標o代表前一個時間的已知數。最後將面上 對流通量的已知數移進源項。 3.1.1.2.3 擴散項(Diffusion Term) 將(3.3)式中的擴散項離散化，可得下式：( )
D f f f f f f S V dS V S F μ∇ ⋅ ≈∑
μ ∇ ⋅ =∑
∫∫
K K K Kw
(3.10) 其中μ 代表面上的黏滯係數、f D f F 代表面上的擴散通量。針對非結構性網格問題， 在擴散項的部份我們利用 Over-Relax 的方法來近似原來的面法向量。首先將Sf K 分 成兩部分表示：(
)
f f SK = +dK SK −dK (3.11) 上式等號右手邊第一項為正交向量、等號右手邊第二項為非正交向量。而正交向 量dK的定義如下： 2 f Pnb Pnb f S d S δ δ = ⋅ K K K K K (3.12) 其中δPnb K 定義如《圖 3. 1》所示。將(3.11)式和(3.12)式代入(3.10)式可得：(
)
D f f f f f f f f f F =μ ∇ ⋅VK SK =μ ∇ ⋅ +VK dK μ ∇ ⋅VK SK −dK(
)
(
)
2 f f n n o nb P f f f Pnb f S V V V S d S μ μ δ = − + Δ ⋅ − ⋅ K K K K K K K K (3.13) 其中上標n代表這個時間的未知數、上標o代表前一個時間的已知數。且上式等號 右手邊第一項為正交項、第二項為非正交項。最後同樣地，將面上擴散通量的已 知數移進源項。3.1.1.2.4 源項(Source Term) 在計算中，我們將壓力梯度視為已知值，並放至源項。對於非結構性網格問 題，求解壓力梯度會需要用到邊界上的壓力，因此做以下推導： b P Pnb P −P = ∇ ⋅P δK (3.14) 其中下標b代表壁面上的中點、δKPnb代表P到b距離向量如《圖 3. 2》。 將壓力梯度利用(3.2)式和(3.3)式相同之方法近似，並離散化可得： 1 1 f f b b f f f f b P P S P S P S ≠ ⎛ ⎞ ∇ = = ⎜ + ⎟ ∀
∑
∀ ⎝∑
⎠ K K K (3.15) 其中 f ≠b代表不包含邊界的其餘各面。將(3.15)式代入(3.14)式可得邊界上之壓力： 1 1 1 P f f Pb f b b b Pb P P S P S δ δ ≠ + ⋅ ∀ = − ⋅ ∀∑
K K K K (3.16) 3.1.1.3 線性代數式的整理 將前一小節所離散的非穩態項、對流項、擴散項及源項代回(3.3)式，可整理 成下列的線性代數式： n n o P P nb nb nb A VK =∑
A VK + − ∇ Δ∀S P (3.17) 其中 P nb nb A A t ρΔ∀ = + Δ∑
(
)
2 max , 0 f f nb f Pnb f S A m S μ δ = + − ⋅ K K K(
)
(
)
{
o}
CD UD o o f f f f f f S F F V S d V t γ μ ρ Δ∀ ⎡ ⎤ = −⎣ − ⎦ + ∇ − + Δ∑
K K K K3.1.2 空間守恆定理(Space Conservation Law) 3.1.2.1 空間守恆定理的簡介
在計算移動邊界和移動網格的問題時，有一專門探討網格體積隨時間的變化 是否守恆之定理，稱作空間守恆定理(Space Conservation Law)，表示式如下：
0 g S d d V dS dt
∫∫∫
∀ ∀ +∫∫
⋅ = K Kw
(3.18) 其中Vg K 代表網格速度，詳細的推導過程稍後於 3.1.2.3 節中描述。 Demirdzic et al.【1】已研究空間守恆定理應用於有限體積法，只要能解決人 為造成的質量源(mass source)問題，在任意的網格移動速度和任意的時階皆能達到 精確的模擬結果。 3.1.2.2 質量源(Mass Sources) 處理移動網格的問題時，我們必須計算每個時階的網格體積改變量，《圖 3. 3》 即是利用網格速度來計算時階中的體積改變。若在計算過程中，體積改變量處理 的不完善，將可能產生質量源(mass source)，流場無法滿足空間守恆定理。因此必 須注意是否產生質量源【15】。 為了清楚的介紹質量源，接下來所舉的例子皆以簡單的二維四邊型網格做解 釋。以下即是產生質量源的例子。《圖 3. 4》為分別運用顯式法(explicit scheme)和 隱式法(implicit scheme)，利用網格面上的速度，來計算體積改變量。可發現顯式 法會少算一塊體積，即《圖 3. 4a》中的右上角空白區塊；而隱式法會多算一塊體 積，即《圖 3. 4b》中的右上角紅色區塊。少算或多算到的網格體積改變量，即稱 為質量源。質量源所造成的誤差，將會影響數值模擬的準確度和收斂性。 為了降低質量源造成的誤差，我們可將時階縮小，如《圖 3. 5》。可以發現《圖 3. 5b》的質量源小於《圖 3. 5a》的質量源，但此誤差還是會隨著時間累積。 因此我們直接利用，新舊網格各面的節點位置，算出網格體積改變量。如《圖 3. 6》，我們找出網格一個面移動前後的節點共四點，利用這四點的位置，算出此 面移動造成的體積改變量。每個面都運用此算法，如此一來可精確的算出網格體 積隨時間的改變量，以解決質量源的問題。 另外，在判斷網格體積的增加或漸少方面，我們是利用網格面的移動速度和 面向量內積，取其值之正負來確定網格各面所掃過的體積為增加或減少。最後將 網格各面所掃過的體積相加減，便可算出網格整體體積是增加或減小。3.1.2.3 空間守恆定理 將連續方程式(2.3)式整理如下：
(
g)
0 S d d V V dS dt∫∫∫
∀ ρ ∀ +∫∫
ρ − ⋅ = K K Kw
(3.19) 若將上式假想為流體速度趨近於 0，可得控制體積的形狀和位置隨時間的改變式， 即是空間守恆定理(3.18)式。再透過 3.1.1.2 節提到方式將(3.18)式離散化可得：(
)
1 0 n n g f f V S t + ∀ − ∀ − ⋅ = Δ∑
K K (3.20) 在每個時階中，控制體積各個面所掃過的體積相加減，會等於新跟舊的整體體積 的差值如下： 1 f n n f t t + Δ∀ ∀ − ∀ = Δ Δ∑
(3.21) 其中Δ∀f 為一個面掃過的體積，f =e w n s, , , ，如《圖 3. 7》所示。將(3.21)式代入(3.20) 式可得：(
)
f g f f V S t Δ∀ ⋅ = = ∀ Δ K K (3.22) 因此由(3.6)式所定義的 r f m 可表示如下：(
)
r f f g f f f f f f f m =ρ VK−VK ⋅SK =ρ VK ⋅SK − ∀ρ (3.23) 上式將Vg⋅Sf K K 提出變為∀f ，即是採用 3.1.2.2 節所提到的計算方法，直接算出網格 各面的體積改變量，因此 r f m (流體相對於面上的質量流率)將不會產生質量源。 接下來必須讓連續方程式滿足空間守恆定理，因此將(3.19)式離散： 1 0 n n r f f m t ρ∀ − ∀+ + = Δ∑
(3.24) 將(3.23)式代入上式： 1 1 0 n n n n f f f f f V S t t ρ∀ − ∀+ + ρ ⋅ −⎛⎜ρ ∀ − ∀+ ⎞⎟= Δ∑
⎝ Δ ⎠ K K0 f f f f V S ρ ⋅ =
∑
K K (3.25) 上式即為不可壓縮流之連續方程式的差分形式，且可表示如下：( )
ρV 0 ∇ ⋅ K = (3.26) 上式即為不可壓縮流之連續方程式的微分形式。因此本文的 ALE 統御方程式的連 續方程式(2.3)式將簡化為和網格速度無關的(3.26)式，並且再配合(2.4)式的動量方 程式來描述本文流體的運動模式。 3.1.2.4 整理 SCL 應用於 FVM 及 PISO 有限體積法中，在離散化對流項(3.6)式所定義出的 r f m (流體相對於面上的質量 流率)，必須運用(3.23)式的方式計算，避免造成質量源的累積。本章推導出(3.26) 式連續方程式將運用於稍後介紹的 PISO 壓力速度修正法當中的(3.41)式和(3.47) 式，詳細的過程稍後於 3.1.3.2 節中討論。 3.1.3 PISO 演算法 3.1.3.1 PISO 演算法的簡介 至今已發展出幾種壓力速度修正法，即包含可求解非穩態問題的 PISO 演算 法，其他還有專門求解穩態問題的 SIMPLE、SIMPLER 和 SIMPLEC 等。壓力速 度修正法應用於本研究，主要是為了滿足質量守恆以及修正第三章離散過後的統 御方程式所解出的速度場和壓力值。Issa【16】於 1986 年發展出非疊代的 PISO(Pressure Implicit with Splitting of Operators)壓力速度修正演算法，來計算非穩態可壓縮流的問題。而 Issa 證實 PISO 同時也適合求解需要反覆疊代的穩態流場問題。因此本文採用 PISO 演算法來求解 非穩態不可壓縮流流場。PISO 共包含一個預測步驟和兩次修正步驟，詳細過程將 在後面幾節介紹。 3.1.3.2 速度與壓力之間的耦合 3.1.3.2.1 預測步驟(Predictor Step) 將(3.17)式表示如下：
(
o)
P P nb nb P nb A VK∗=∑
A VK∗ + S− ∇ Δ∀P (3.27) 其中VK∗代表未修正的速度場、Po代表未修正的壓力值。利用現有(guessed)的壓力 值 o P 求解上式動量方程式，可得到速度場VK∗，但並不滿足連續方程式。因此我們 必須利用後面的兩次修正步驟，修正速度場和壓力值來滿足連續方程式和空間守 恆定理。 3.1.3.2.2 第一次修正步驟(1st Corrector Step) 將(3.25)式移項整理過後，主格點速度場可表示如下： nb nb o o nb P P P P p P P P P A V S V P H P A A A ∗ ∗ ∗ + ⎛ ⎞ ⎛ ⎞ Δ∀ Δ∀ = −⎜ ⎟ ∇ = −⎜ ⎟ ∇ ⎝ ⎠ ⎝ ⎠∑
K K K (3.28) 其中HP ∗ K 為Vnb ∗ K 的函數。 同樣的，控容面上的速度與壓力關係式可透過(3.28)式表示為： o f f f P f V H P A ∗ = ∗ −⎛Δ∀⎞ ∇ ⎜ ⎟ ⎝ ⎠ K K (3.29) 另外 f f f P f H V P A ∗ = +⎛Δ∀⎞ ∇ ⎜ ⎟ ⎝ ⎠ K K (3.30) 將(3.30)式代入(3.29)式，控容面上的速度與壓力關係式最後可表示成： o f f f f P f P f V V P P A A ∗ = +⎛Δ∀⎞ ∇ −⎛Δ∀⎞ ∇ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ K K (3.31) 上式括弧中的值是取主格點和鄰近格點的平均值，如下： 1 2 P f P P P nb A A A ⎡ ⎤ ⎛Δ∀⎞ ⎛Δ∀⎞ ⎛Δ∀⎞ = ⎢ + ⎥ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎣⎝ ⎠ ⎝ ⎠ ⎦ (3.32) 且(3.31)式中，上標”─”表示由主格點P與控容面 f 相鄰之格點nb內差而得，如下：(
1)
f P nb P P P w P w P ∇ = ∇ + − ∇ (3.33)(
1)
f P nb P P VK =w VK + −w VK (3.34) 其中wP為權重因子(weighting factor)。 未修正的面上質量流率mf ∗ ，可寫成下式： f f f f m∗ =ρ VK∗ ⋅SK(
)
2 f f f f f nb P f Pnb P f Pnb f S V S P P P A S ρ ρ δ δ ⎛ ⎞ ⎛Δ∀ ⎜⎞ ⎟ ⎡ ⎤ = ⋅ − ⎜ ⎟ ⎣ − − ∇ ⋅ ⎦ ⎜ ⋅ ⎟ ⎝ ⎠ ⎝ ⎠ K K K K K K (3.35) 而上式也必需避免產生質量源，因此利用(3.23)式的方法，可得：(
)
2 f f f f f f nb P f Pnb P f Pnb f S m V S P P P A S ρ ρ δ δ ∗ = ⋅ − ⎛Δ∀ ⎜⎞ ⎛ ⎞⎟ ⎡ − − ∇ ⋅ ⎤ ⎜ ⎟ ⎜ ⋅ ⎟⎣ ⎦ ⎝ ⎠ ⎝ ⎠ K K K K K K (3.36) 依照(3.28)式的形式，修正後的速度與壓力關係式可表示為： P P P P P V H P A ∗∗= ∗−⎛Δ∀⎞ ∇ ∗ ⎜ ⎟ ⎝ ⎠ K K (3.37) 其中VP ∗∗ K 為第一次修正後的速度、PP ∗為第一次修正後的壓力。 將(3.37)式與(3.28)式相減，可得速度修正量和壓力修正量的關係式如下： P P P P V P A ⎛Δ∀⎞ ′ = −⎜ ⎟ ∇ ′ ⎝ ⎠ K (3.38) 其中VK′ ≡VK∗∗−VK∗，∇ = ∇ − ∇PP' PP∗ PPo，上標” ' ”即為我們的修正量。上述的方程式 是在網格主格點上推出，其亦可使用在每個主控面上。 f f P f V P A ⎛Δ∀⎞ ′= −⎜ ⎟ ∇ ′ ⎝ ⎠ K (3.39) 第一次修正後的流量值可表示如下： f f f f f f f m m m m ρ V S ∗∗ ∗ ∗ ′ = + ′ = + ⋅ K K(
)
(
)
2 f f f nb P f f f P f Pnb f P f S m P P P S d A S A ρ ρ δ ∗ ⎛Δ∀⎞ ′ ′ ⎛Δ∀⎞ ′ = − ⎜ ⎟ − − ⎜ ⎟ ∇ ⋅ − ⋅ ⎝ ⎠ ⎝ ⎠ K K K K K (3.40) 為了滿足空間守恆定理，第一次修正後的流量值必須滿足(3.25)式，如下： 0 f f m∗∗ =∑
(3.41) 最後將(3.40)式代入(3.41)式，可整理成下列的線性代數式： 1 1 1 2 P P P P nb nb P P f A P′=∑
A P ′+S +S (3.42) 其中上標 P 代表壓力修正方程式、上標”1”代表第一次修正； 2 f P nb f P f Pnb f S A A S ρ δ ⎛Δ∀⎞ = ⎜ ⎟ ⋅ ⎝ ⎠ K K K P P P nb f A =∑
A 1 1 P f f S = −∑
m∗(
)
2 2 P f f f f P f S P S d A ρ ⎛Δ∀⎞ ′ = ⎜ ⎟ ∇ ⋅ − ⎝ ⎠∑
K K 3.1.3.2.3 第二次修正步驟(2nd Corrector Step) 透過求解(3.42)式，進而得到了第一次修正過後的速度VK∗∗和壓力值P∗。第二 次修正步驟中，主格點與相鄰格點速度與壓力關係式可表示為： P P P P P V H P A ∗∗∗= ∗∗−⎛Δ∀⎞ ∇ ∗∗ ⎜ ⎟ ⎝ ⎠ K K (3.43) 將(3.43)式和(3.37)式相減，速度修正量及壓力修正量可表示成下面關係： P P P P P V H P A ⎛Δ∀⎞ ′′= ′ −⎜ ⎟ ∇ ′′ ⎝ ⎠ K K (3.44) 其中VK′′ ≡VK∗∗∗−VK∗∗，VK′ ≡VK∗∗−VK∗，∇ = ∇PP′′ PP∗∗− ∇PP∗。 因此流量修正可表示成下式：(
)
(
)
2 nb nb f f f f nb P f f f f f P f Pnb f P f P A V S m P P P S d S A S A A ρ ρ ρ δ ′ ⎛Δ∀⎞ ⎛Δ∀⎞ ′′ = − ⎜ ⎟ ′′− ′′ − ⎜ ⎟ ∇ ′′⋅ − + ⋅ ⋅ ⎝ ⎠ ⎝ ⎠∑
K K K K K K K (3.45) 第二次修正後的流量值可表示如下： f f f m∗∗∗=m∗∗+m ′′ (3.46) 為了滿足空間守恆定理，第二次修正後的流量值也必須滿足(3.25)式，如下： 0 f f f f f f m∗∗∗ = m∗∗+ m′′ =∑
∑
∑
(3.47) 將上式結合(3.41)式可得： 0 f f m ′′ =∑
(3.48) 將(3.45)式代入(3.48)式，可整理成下列的線性代數式： 2 2 1 2 P P P P nb nb P P f A P′′ =∑
A P ′′+S +S (3.49) 其中係數 P P A 、AnbP和(3.42)式之係數相同，上標”2”代表第二次修正； 2 1 nb nb f P f f P f A V S S A ρ ⎛ ′⎞ ⎜ ⎟ = − ⎜ ⎟ ⋅ ⎜ ⎟ ⎝ ⎠∑
K K(
)
2 2 P f f f f P f S P S d A ρ ⎛Δ∀⎞ ′′ = ⎜ ⎟ ∇ ⋅ − ⎝ ⎠∑
K K 3.1.3.2.4 壓力修正方程式的整理 在兩次的修正步驟中，為了使結果更可靠，我們將每次的修正步驟再分成兩 部分運算：第一部分僅包含源項SP1，而第二部分僅包含SP2。求解壓力修正代數 方程式順序如下： 第一次修正 正交修正 ： 1 1 P P nb nb P f A P′=∑
A P ′+S 非正交修正 ： 1 2 P P nb nb P f A P′ =∑
A P ′+S第二次修正 正交修正 ： 2 1 P P nb nb P f A P′′=
∑
A P ′′+S 非正交修正 ： 2 2 P P nb nb P f A P′′=∑
A P ′′+S 其中SP2為非正交的源項，對於正交網格SP2 =0。另外，對於幾何外型較為複雜的 非結構性網格問題，可重複比較多次的非正交修正。 3.1.3.3 邊界條件的給定 3.1.3.3.1 進出口壓力邊界的流量計算 欲求解壓力邊界的流量，必須先在進出口邊界處設定壓力值，再利用此壓力 值修正出入口的流量。步驟如下： 第一步：將壓力邊界相鄰的網格點上之壓力，設定成給定的值。 第二步：再求解兩次壓力修正方程式中，將壓力邊界相鄰的網格點上之壓力修正 值設成 0，即P′=0如《圖 3. 8》。因此壓力修正方程式的係數和源項必須 做以下的調整： 1 2 P P P P nb nb P P f A P′=∑
A P ′+S +S (3.50) 其中SP1 =SP2 =0， 0 P nb A = 第三步：求解完壓力修正方程式，最後讓壓力邊界相鄰之網格滿足質量守恆定理： f 0 f m =∑
(3.51) 如《圖 3. 9》所示，壓力邊界上的質量流率可表示為下式： mf4 = −(
mf1+mf2+mf3)
(3.52) 3.1.3.3.2 進出口壓力邊界的速度計算 欲計算壓力邊界上的進出口速度，本文採用對流邊界條件(Convective Boundary Condition)，其數學式子為： 0 c u t x φ∗ φ∗ ∗ ∗ ∗ ∂ + ∂ = ∂ ∂ (3.53) 將上式離散後可得：0 n o n n b b b nb c u t x φ∗ φ∗ φ∗ φ∗ ∗ ∗ ∗ − + − = Δ Δ (3.54) 其中上標n、o代表新跟舊的值，下標b、nb代表壓力邊界和其相鄰網格的值，uc∗ 代表無因次之對流流速。將(3.54)式整理如下： 1 o n n b nb b Cr Cr φ φ φ∗ = ∗ + ∗ + (3.55) 其中Cr uc t x ∗ ∗ ∗ Δ = Δ ，Cr為 Courant number。 3.1.3.3.3 牆邊界條件的計算 流體力學計算中的流道上下壁面和壓電振動平板皆視為牆邊界條件。流道上 下壁面之牆是固定不動，因此牆速度Vwall和流體速度V 皆為 0，即(V =0、Vwall =0)。 而壓電振動平板的牆邊界條件中的牆是會隨著時間改變位置，因此數學式子為牆 的速度和流體速度相等，即(Vwall =V)。 3.2 固體力學數值方法
3.2.1 有限元素法(Finite Element Method)軟體基本架構
有限元素法即是將一個實體的物件，從點開始、線、面、體積將物體的實際 情況描述好後，將其分割為不同大小、不同形狀、不同材料和不同的區域，而每 個小區域稱之為元素，再進行處理。而一個實體物件不可能分割為無限個元素， 皆為有限的分割，因此稱為有限元素法。 分析者利用自己領域的需求，將每一個元素節點上的作用力方程式推導出 來，並將整個物體中的所有元素組合，構成最後的系統聯立方程式進行求解。 一般完整的有限元素法分析程式必須包含前處裡(Preprocessing)、求解器 (Solution)和後處理(Postprocessing)三大部分。 前處理：選用適當元素、定義材料特性矩陣、建立模型、元素切割。 求解器：選擇分析類型、定義限制條件、求解線性代數方程式。 後處理：將求解結果例如：位移、應變、應力、頻率等資料，利用動畫、等 位圖形、表格、數據歸納表…等不同方式將結果呈現出來。
