壁面蠕動具相位差之管道流場數值分析

(1)

國立交通大學

機械工程學系

碩士論文

壁面蠕動具相位差之管道流場數值分析

Numerical Analysis of the Flow in a Peristaltic Channel

with Phase Shift

研究生:鄭凱鴻

指導教授: 崔燕勇教授

(2)

壁面蠕動具相位差之管道流場數值分析

Numerical Analysis of the Flow in a Peristaltic Channel with Phase Shift

研究生：鄭凱鴻 Student：Kai Hung Cheng

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

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

碩士論文

A Thesis

Submitted 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 2012

(3)

壁面蠕動具相位差之管道流場數值分析研究生:鄭凱鴻指導教授:崔燕勇博士

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

摘要

本研究利用非結構性移動網格和有限體積法去求解二維蠕動式管道流場，之後探討上/下管壁進行一蠕動波以不同相位差、雷諾數、振幅比、波數與背壓造成流場中平均流量的影響。再將數值方法結果與潤滑模型中理論解比較，即得知在小波數模擬計算之結果與理論解差異甚小。根據數值結果得知隨著相位差增加則平均流量會遞減，且在相位差為零有最大平均流量。在逆壓與順壓的情況下，由數值結果得知提高振幅比與波數皆可提高平均流量的輸出。再經數值結果發現在較小波數與背壓為零情況下，雷諾數對平均流量的影響幾乎無關，然而高波數的情形下平均流量隨著雷諾數增加；另外在固定背壓與振幅比的情形下，發現波數與平均流量之間有遞增的關係。並且從數值結果得知隨著背壓增加則平均流量減少，兩者之間呈線性遞減的關係。最後運用徑線圖配合著流線與速度分佈，觀察管道中流體的瞬時流場與質點運動情況，使本研究更加清楚蠕動式管道內流體運動的情形。

(4)

Numerical Analysis of the Flow in a Peristaltic Channel with Phase Shift

student: Kai Hung Cheng advisor:Dr.Yeng-Yung Tsui

Department of Mechanical Engineering National chiao Tung University

ABSTRACT

In this study, the two-dimensional peristaltic channel flow is solved by using the unstructured moving grid and the Finite Volume method. We then

investigate the influence on the mean flow rate caused by the peristaltic waves in the upper/lower wall under different phase shift, Reynolds number, wave number, amplitude ratio and backpressure.Furthermore, when comparing the numerical results with the theory result of the lubrication model, we know their difference is very few when with smaller wave numbers. Based on the numerical results, we found that the average flow rate will decrease when the phase shift increasing, while the mean flow is maximized when the phase shift is at zero.Under the condition of adverse pressure and positive pressure, we can get the conclusion from the

calculated results that the output of mean flow rate can be enhanced by increasing the amplitude ratio and the wave number. Furthermore, the Reynold number has little effect on mean flow rate under the low wave number and zero back pressure condition, but the mean flow rate will increase in line with Reynold number under high wave number condition.When fixing the backpressure and the amplitude ratio, the wave number has liner-increasing relation with the mean flow rate.Based on the numerical result, we further know that the mean flow rate has linear- decreasing relation with backpressure: mean flow rate will decrease when the backpressure increasing.Finally, by using streamline chart as well as pathline and velocity profile, we can observe the transient flow and particle motion in the channel and get a clearer picture of how fluid moves in the peristaltic channel.

(5)

誌謝

非常感謝崔燕勇教授在求學路上細心的叮嚀和不吝嗇教導我們，不僅讓我們理解做研究最重要就是嚴謹的態度，還有在往後的日子如何用積極正確的態度去面對人生的各種問題。另外我還要感謝林仕文學長總是在課業上幫助我們解決問題，以及胡育昌和郭大慶學長總是在半夜的時候出現在實驗室，讓總是晝伏夜出的我在半夜上出現問題時教導我，還有在我們有任何問題時都能適時的給予幫助的賴胤男、黃義政、林子翔學長，以及吳奉起和黃裕堂和陳虹汝的陪伴才可以讓我在新竹的日子過的如此充實，最後感謝黃俊龍、葉保杰和丁崑哲學弟你們讓實驗室充滿了歡樂的氣氛，雖然丁崑哲學的總是在實驗室愛放著王識賢的經典台語老歌，並摧殘著我在課業上剩餘的腦細胞，但也讓我發現少女時代的魅力，使我在這所校園裡感受到異性的溫暖也鼓勵著我更有力氣應付於課業。還有江健銘總是分享自己真人真事的有趣經驗給我。還有替我解決英文知識上的問題的江珉磬，以及很多朋友的幫助和愛護這本論文才可以完成。最後感謝家人們的支持，讓我可以隨意的做自己想做的事情，雖然我在新竹的短短的日子裡發生很多的事情，睡覺睡到地板爆掉、騎車騎到手機噴掉、開個鐵捲門也可以讓門毀掉還有大年初一的好日子前女友提分手和剝皮妹的故事等..太多的事情，但是我相信這只是一個過渡期，讓我更可以在將來面對人生各種問題上有更高的抗壓性，總之太多需要感謝的人了才可以讓我完成學業。

(6)

目錄中文摘要...i 英文摘要...ii 致謝...iii 目錄...iv 表目錄...vii 圖目錄...viii 符號說明...x 第一章緒論...1 1.1 前言...1 1.2 蠕動式傳輸現象簡介...1 1.3 文獻回顧...2 1.4 研究目的...6 第二章數學分析與物理模型...7 2.1 簡介...7 2.2 流場中的假設與物理模型...7 2.3 統御方程式...7 2.4 無因次化... ...8 2.5 邊界條件...9

(7)

第三章空間守恆定理...11 3.1 簡介... ...11 3.2 空間守恆定理...11 第四章統御方程式的離散...13 4.1 離散動量方程式...13 4.1.1 非穩態項...13 4.1.2 對流項...13 4.1.3 擴散項...14 4.1.4 源項...15 4.2 代數方程式...16 第五章數值法法...17 5.1 預設步驟...17 5.2 第一次修正步驟... ...17 5.3 第二次修正步驟... .... ...20 5.4 入出口壓力邊界的流量計算...21 5.5 入出口壓力邊界的速度計算...22 5.6 PISO 演算法流程...22 第六章結果與討論...24 6.1 網格測試...24

(8)

6.2 週期性穩態之時間測試...25 6.3 潤滑模型的理論解...25 6.4 數值方法之計算結果...26 6.5 流場結構的影響... ...33 6.5.1 徑線與速度分佈...33 6.5.2 沿管道中心線之瞬時壓力與平均壓力分布...35 第七章結論...36 參考文獻...38 附錄...40

(9)

表附錄表 6.1(a)振幅比ε=0.1 平均流量網格測試...42 表 6.1(b)振幅比ε=0.3 平均流量網格測試...42 表 6.1(c)振幅比ε=0.6 平均流量網格測試...42 表 6.2(a)振幅比ε=0.1 之網格測試誤差百分比...43 表 6.2(b)振幅比ε=0.3 之網格測試誤差百分比...43 表 6.2(c)振幅比ε=0.6 之網格測試誤差百分比...43 表 6.3 週期性平均流量測試... ...44

(10)

圖目錄圖 2.1 蠕動管道示意圖...45 圖 6.1 在不同的雷諾數，相位差對平均流量的影響...46 圖 6.2 在不同的波數，相位差對平均流量的影響...46 圖 6.3 在不同的背壓，下相位差對平均流量的影響...46 圖 6.4 相位差=0、波數α=0.1 於流場 1/4 管道長之流線圖………47 圖 6.5 相位差 = /2、波數α=0.1 於流場 1/4 管道長之流線圖...48 圖 6.6 相位差 、波數α=0.1 於流場 1/4 管道長之流線圖...49 = 圖 6.7 相位差=0、波數α=0.1 於流場 1/4 管道長之流線圖...50 圖 6.8 相位差 = /2、波數α=0.1 於流場 1/4 管道長之流線圖...51 圖 6.9 相位差 、波數α=0.1 於流場 1/4 管道長之流線圖...52 = 圖 6.10 在不同的相位差下，雷諾數對平均流量的影響...53 圖 6.11 在不同的波數下，雷諾數對平均流量的影響...53 圖 6.12 在不同的背壓下，雷諾數對平均流量的影響...53 圖 6.13 在不同的雷諾數 Re=1、100 下，振幅比對平均流量的影響...54 圖 6.14 在不同的波數α=0.1、0.3125 下，振幅比對平均流量的影響....54 圖 6.15 在不同的背壓 Pb=-1.0、0、1.0 下，振幅比對平均流量的影響..54 圖 6.16 相位差 = /2、振幅比ε=0.1 於流場 1/4 管道長之流線圖...55 圖 6.17 相位差 = /2、振幅比ε=0.8 於流場 1/4 管道長之流線圖...56 圖 6.18 在不同雷諾數 Re=1、20、100 下，波數對平均流量的影響...57 圖 6.19 在不同的振幅比ε=0.1、0.6 下，波數對平均流量的影響...57 圖 6.20 在不同背壓 Pb=-1.0、0、1.0 下，波數對平均流量的影響...57 圖 6.21 相位差 = /2、波數α=0.1 於流場 1/4 管道長之流線圖...58 圖 6.22 相位差 = /2、波數α=0.3125 於流場 1/4 管道長之流線圖...59 圖 6.23 在不同相位差下，背壓對平均流量的影響...60 圖 6.24 在不同波數下，背壓對平均流量的影響...60 圖 6.25 在不同振幅比下，背壓對平均流量的影響...60 圖 6.26a 在相位差=0與背壓 Pb=0 下，流體質點的運動路徑...61 圖 6.26b 在相位差 = /2與背壓 Pb=0 下，流體質點的運動路徑...61 圖 6.26c 在相位差 = 與背壓 Pb=0 下，流體質點的運動路徑...61 圖 6.27a 在相位差=0與背壓 Pb=0 下，管道於流場的速度分佈...62 圖 6.27b 相位差 = /2與背壓 Pb=0 下，管道於流場的速度分佈...62 圖 6.27c 相位差 = 與背壓 Pb=0 下，管道於流場的速度分佈...62 圖 6.28a 在相位差=0與背壓 Pb=1.0 下，流體質點的運動路徑...63 圖 6.28b 在相位差 = / 2與背壓 P=1.0 下，流體質點的運動路徑...63

(11)

圖 6.28c 在相位差 = 與背壓 Pb=1.0，流體質點的運動路徑...63 圖 6.29a 相位差=0與背壓 Pb=1.0 下，管道於流場的速度分佈...64 圖 6.29b 相位差 = /2與背壓 Pb=1.0 下，管道於流場的速度分佈...64 圖 6.29c 在相位差 = 與背壓 Pb=1.0，管道於流場的速度分佈...64 圖 6.30a 在振幅比ε=0.1，流體質點的運動路徑...65 圖 6.30b 在振幅比ε=0.8，流體質點的運動路徑...65 圖 6.31a 振幅比ε=0.1 下，管道於流場的速度分佈...66 圖 6.31b 振幅比ε=0.8 下，管道於流場的速度分佈...66 圖 6.32a 在波數α=0.1，流體質點的運動路徑...67 圖 6.32b 在波數α=0.3125，流體質點的運動路徑...67 圖 6.33a 波數α=0.1 下，管道於流場的速度分佈...68 圖 6.33b 數α=0.3125 下，管道於流場的速度分佈...68 圖 6.34 不同相位差下，沿著中心線的壓力分佈...69 圖 6.35 不同振幅比下，沿著中心線的壓力分佈...69 圖 6.36 不同波數下，沿著中心線的壓力分佈...69

(12)

符號說明

v ：流體速度 g v ：網格速度 c v ：對流速度 wall v ：壁面速度  ：流體密度  ：流體黏滯係數 p ：流場壓力 t ：時間 1 y ：上壁面邊界移動方程式 2 y ：下壁面邊界移動方程式 h ：管道入口的一半高度 a ：振幅  ：波長 c_：波速 * v ：無因次流體速度 * g v ：無因次網格速度 * 1 y ：無因次上壁面邊界移動方程式 * 2 y ：無因次下壁面邊界移動方程式

(13)

* p ：無因次壓力 * t ：無因次時間  ：振幅比  ：波數 Re ：雷諾數 T ：週期

(14)

第一章緒論 1.1 前言近年來在醫學上或者工業上，蠕動式傳輸應用在幫浦或者微流道來輸送流體被廣泛的研究。在蠕動式幫浦上藉由蠕動式傳輸來達成計量精準且免受汙染上，而在微流道上蠕動式傳輸有益於在多種不同試劑上進行混合與相互作用藉此可以增加效率，不同以往在微流道時只依賴擴散作用來達成目的，而此舉節省在工業上的生產時間。在生物力學上，大多數的生物系統都藉由蠕動式傳輸來運送流體，舉例如尿從腎臟傳送到膀胱或者食糜在胃腸道的移動、精子在輸精管的輸出卵子在子宮頸的傳輸、淋巴液在淋巴管的傳輸或者血液在血管的舒張與收縮以及心肺機的製造也都是使用蠕動傳輸原理，也因此探討蠕動管道傳輸現象變得相對的重要。 1.2 蠕動式傳輸現象簡介蠕動式傳輸現象即為一行進波沿著韌性管壁在管道上運動，其管道內部充滿流體，藉由管壁連續收縮與擴張來輸送流體，且此傳遞波為連續且週期性的運動。一般我們探討二維的蠕動傳輸現象，在本研究中，我們將管壁分為具對稱性與非對稱性的蠕動管道流，即上/下壁面是否為對稱性。其中管壁面上之變形在計算上因牽扯流固耦和較為複雜，因此假設為剛性壁面。

(15)

1.3 文獻回顧: 在文獻回顧中分為潤滑理論和擾動理論與其他數值方法，並於下列介紹文獻之不同。 (1) 潤滑理論 Shapiro et al. [1]文中使用二維、對稱性蠕動管道，探討在極長的波長和低雷諾數下的蠕動傳輸現象，求解Ｎ-S 方程式，在求解中忽略對流項得出分析解，此種假設又稱為潤滑理論（lubrication theory）；相較於固定座標下，使用波座標為參考座標觀察，則可將流場視為穩態;即可發現陷流

(Trapping)和渦流(bolus)的現象，而此結果與 Weinberg et al. [6]的實驗相符合。

Jiménez-Lozano and Sen [2]探討二維、不可壓縮之牛頓流體，且管壁上/下為對稱管壁，依照潤滑理論的假設，求解簡化後的 N-S 方程式，文中使用波座標為參考座標，並觀察蠕動時

造成逆流（reflux）、陷流(Trapping)和增流的現象，然後討論

三種現象的分歧點穩定性。

Ramachandra Rao and Mishra [3]說明一個二維、不可壓縮之牛頓流體且非對稱性的蠕動式管道，即上/下管壁具有相位差時，依照潤滑理論的假設，求解簡化後的 N-S 方程式，文中

(16)

探討不同的管壁高度、波振幅與相位差造成平均流量的影響，並以波座標為參考座標觀察逆流（reflux）與陷流(trapping) 的現象，進而發現陷流現象只存在於管道截面具有差異時才發生。

Hayat and Ali [4]文中說明生物上的蠕動傳輸行為，若將黏滯性設為定值，則似乎不太接近實際案例；因此文中將黏滯性設為一對空間的的變數，藉此探討二維非對稱性的蠕動管道，依照潤滑理論的假設，使用 Mathlab 軟體去求解簡化後的 N-S 方程式，並討論不同管道高度、振幅、相位差與黏滯性對管道中壓力梯度造成的影響，並發現當黏滯性上升時，壓力梯度相對提高。

Eytan and Elad [5]文中為探討子官內流體運動，藉此利用一基本蠕動式管道進行分析，依照潤滑理論的假設，求解簡化後的 N-S 方程式，並討論在不同相位角時，造成截面上的速度分布與軸上壓力分佈的影響，然後以波座標為參考座標，分析在不同振幅比產生的逆流和陷流現象，其結果在上/下管壁為具對稱性時，此時可找出最大速度與流量；當相位差上升至π 時，此時管壁截面不改變發現最小速度與流量。

(17)

et al. [1]的理論解，利用實驗將其一端設為固定壁，另一端為行進 Sin 波的彈性壁，藉此觀察平均流量與壓力梯度;並運用影像技術找出產生陷流與回流現象的區域，其結果與 Shapiro 的理論相符合。

(2) 擾動理論

Fung and Yih.[7]文中探討一個二維、不可壓縮之牛頓流體的蠕動管道，且管道兩側設為定壓力下，利用振幅與管道高度的比值，以冪(power)進行展開，則振幅需很小，但相較於 Shapiro 限制下，此模型可得到較高雷諾數的分析解，而此法又稱為擾動理論，文中探討當壓力梯度達到一臨界值則會產生逆流的現象。

Yin and Fung [8]文中主要透過實驗方法去驗證 Fung and Yih.[7]的理論是否吻合，利用實驗儀器將上壁面設為固定牆，下壁為一行進 Sin 波，藉此使管道內的流體傳輸，經由修正 Fung and Yih. [7]的理論解與實驗數據比較，其結果發現在較高振幅時與理論解誤差甚大，並透過影像技術觀察一截面的速度分布情形。

Usha and Ramachandra Rao [9]文中藉由兩種不相容的液體與不同的黏滯係數比例進行二維的蠕動管道傳輸分析，並將擾動級數以管高和波長的比值進行展開去求解流線函數，使用波座標為參考座標討論波數、黏滯係數比例、雷諾數對陷流現象的影響，結果發現當渦

(18)

流（eddy）產生時，渦流的大小與黏滯係數比例和波數呈正相關。 Wilson and Panton [10]文中研究二維非對稱性的蠕動管道傳輸現象，探討一純（pure）彎曲的蠕動管道，即管道彎曲時，管道截面不隨彎曲而改變，利用擾動理論將波振幅與波長的比值去求解 N-S 方程式，進而得出一分析解，且此假設需較小波振福的前提下。 (3)其他數值方法解

Tong and Vawter [11]文中使用有限體積法去求解一軸對稱的蠕動管道傳輸，並討論波長與逆流的現象，結果發現當無淨流量時逆流發生於中心軸上。 Pozrikidis [12]在蠕動流場中，假設流場為二維且遲緩 (Creeping Motion)下，利用邊界積分法來求解簡化後之 N-S 方程式，文中探討管道高度、波振幅、相位差以及平均壓力梯度造成流場的影響。並比波座標分析流場情形，討論振幅比和相位差會影響陷流和渦流現象，結果發現渦流大小和振幅比與波數呈正相關。

Takabatake and Ayukawa [13]在蠕動流場中，使用有限差方去求解 N-S 方程式，並探討波振幅與波長、雷諾數對流場的影響。而在此篇中無限制雷諾數波長與波振幅，所以相對其他而言可到較廣泛的結果。

(19)

壁面為固定牆，透過非結構性網格與有限體積法來離散統御方程式，需滿足空間守恆定理，並使用 PISO 演算法求解 N-S 方程式。藉此來探討背壓、波數、振幅比等參數對平均流量的影響，與本研究不同是將引進相位差的變化對平均流量的影響，藉此探討各個不同參數對平均流量的影響，再與其 Shapiro et al. [1]的理論解做比較。並運用流線圖，搭配著徑線與速度分布觀察管道內之流場運動情形。 1.4 研究目的隨著工業與醫學上迅速地發展與需要，對於蠕動式傳輸的現象也相當重視，藉由蠕動式傳輸可增加微管道時的擴散效應與了解蠕動管壁的流場特性藉此運用在醫學與工業上。本文主要在使用數值模擬方法來了解蠕動式管道流場的特性。藉由非結構網格和有限體積法對統御方程式的離散以及滿足空間守恆定理條件下，使用以 PISO 演算法求解 N-S 方程式。在本文主要研究具相位差時蠕動之管道流場，並探討不同相位差時，造成蠕動傳輸管道在一個週期平均流量的影響並藉引進潤滑理論比較之。最後研究五種不同的參數:振幅比、波數、雷諾數、相位差、背壓，隨著五個不同的參數設定下，造成平均流量的變化，藉此探討對管壁內流場的情形。並引進徑線與管道內之速度分布配合著流線圖，則更清楚管道內流體運動情形。

(20)

第二章數學分析與物理模型 2.1 簡介本章主要介紹求解流場時所需要的物理模型與數學模式。2.2 節描述流場中的假設與物理模型；2.3 節裡描述流場所使用之統御方程式以及移動邊界方程式；2.4 節無因次化分析統御方程式；2.5 節邊界條件。 2.2 流場中的假設與物理模型本研究中，模擬蠕動管道流之物理模型，如圖(2.1)且管壁為一週期性之時變函數的波動，藉由管壁蠕動來驅動流體進/出管道而達到輸送的目的。故假設流體為: (1) 牛頓流體(Newton Fluid) (2) 不可壓縮流(Incompressible Flow) (3) 二維非穩態(Unsteady) (4) 忽略重力項 (5) 層流 2.3 統御方程式 2.3.1 ALE 座標系統數學模型 在計算流體力學時，若一般情形下網格靜止不動，則可使用 Eulerian 座標系統，因此利用 ALE 座標系統。在 ALE(Arbitary Lagrangian

(21)

-Eulerian)座標下若網格移動速度與流體相同時，則可視為 Lagrangian 座標系統，倘若網格靜止不動可視為 Eulerian 座標系統。所以 ALE 座標系統適用於分析移動邊界上的問題。其 ALE 座標系統統御方程式如下: 連續方程式: 0 v ∇ =i _(2.1) 動量方程式:

(

( g

)

( )

v v v v p v t ρ∂ + ∇ ρ − = −∇ + ∇ μ∇ ∂ i i _(2.2) 其中μ為流體黏滯係數、ρ為流體密度、∇p為壓力梯度、v為流體速度、vg為網格移動速度。比較 Eulerian 座標系統，即發現 ALE 座標系統在動量方程式中多了網格移動速度項vg。在本文中，移動邊界的壁面之運動方程式是隨時間變化且連續由左往右的固定波速行進的弦波，其移動邊界數學式為: 1 2 (x,t) cos ( ) y h a π x ct λ = + − _{；為上壁面 (2.3)}

( )

(

)

2 2 , cos[ ] y x t h a π x ct ϕ λ = − − − − _{；為下壁面 (2.4)} 其中λ為波長、c為波速，h 為管道的一半高度，a 為波動振幅，ϕ 為相位差角度。 2.4 無因次化為了易於分析，此節將御方程式做無因次化的處理，取管道高度的一

(22)

半長度a和移動方程式中的波速c為特徵長度、速度，則可定義無因次化參數如下: x x h ∗ ₌ _， y y h ∗ ₌ _， * u u c = ， * v v c = ， * g g v v c = ， * 2 p p c ρ = ， ct t h ∗ ₌ _， 1 = h ∗ ∇ ∇ 其無因次化統御方程式經過整理後可寫成下列通式: 連續方程式: * * ( )v 0 ∇ i = (2.5) 動量方程式:

(

)

( )

* * * * * * * * * * * 1 ( Re g v v v v p v t ∂ _{+ ∇} ₋ _{= −∇} ₊ _{∇ ∇} ∂ i i _(2.6) 其中無因次參數雷諾數Re定義為Re ρch μ ≡ 。移動邊界:

(

)

(

)

* * * * * 1 , 1 cos 2 y x t = +ε πα x −t _{；為上壁面 (2.7)}

(

)

(

)

* * * * * 2 , 1 cos (2 ) y x t = − −ε πα x −t −ϕ ；為下壁面 (2.8) 其無因參數振幅比 A h ε = 、波數α ≡ h λ 。 2.5 邊界條件 (A)出入口邊界條件: 在本文中，考慮出入口的兩端在設定為壓力邊界條件下，由於不同的壓力差造成淨流量變化，而出入口的速度則由對流邊界做計算:

(23)

壓力邊界條件 :

Δ =

P

定值 (B)牆邊界條件:

(24)

第三章空間守恆定理 3.1 簡介在本文中計算移動邊界與移動網格時，會造成網格體積隨時間變化，而探討此一學問為空間守恆定理。其表示式如下: 0 g s d d v d s dt

∫∫∫

∀ ∀ +

∫∫

i = (3.1) 其中vg代表網格速度 3.2 空間守恆定理首先由連續方程式表示:

( )

0 div v tρ ρ ∂ ₊ ₌ ∂ (3.2) 利用有限體積法進行體積分得如下:

( )

0 d div v d t ρ ρ ∀ ∀ ∂ ⎛ _{⎞ ∀+} _{∀ =} ⎜_∂ ⎟ ⎝ ⎠

∫∫∫

(3.3) 再經由 Gauss 散度定理與 Leibniz 的微分法則可將上式整理成: ( g) 0 s d d v v d s dt

∫∫∫

_∀ ρ ∀ −

∫∫

ρ − i = (3.4) 從上式可知，若假設流體速度趨近於 0 時，即是空間守恆定理(3.1) 式再將(3.4)式進行離散化可得:

( )

0 n o g f f v s t ∀ − ∀ − = Δ

∑

i (3.5) 在每個時階中，新與舊的整體體積差值，會等於控制體積各個面所掃過的體積加減。

(25)

f n o f t t Δ∀ ∀ − ∀ ₌ Δ Δ

∑

(3.6) 將(3.6)式代入(3.5)式可得:

( )

g f v s t Δ∀ = Δ i (3.7) 接著讓連續方程式滿足空間守恆定理，因此將(3.4)式離散

( ) ( )

₍

₎

0 n o g f f v v s t ρ ρ ρ ∀ − ∀ _⎡ _⎤ + _⎣ − _⎦ = Δ

∑

i (3.8) 再將(3.7)式代入上式:

( ) ( )

0 f f n o _n _o f f f v s t t ρ∀ − ∀ρ ₊ _ρ ₋_ρ ⎛∀ − ∀ ⎞₌ ⎜ ⎟ Δ

∑

i _⎝ Δ _⎠ (3.9) 即知下式: =0 f f f f v s ρ

∑

i (3.10) 上式即為不可壓縮流之連續方程式的差分型式，則可表示為: 0 v ρ ∇i = (3.11) 其上式為不可壓縮流之連續方程式的微分型式。當在不可壓縮條件下進一步得出: 0 v ∇ =i (3.12) 由此發現在滿足空間守恆定理條件下推導出的連續方程式，其與移動網格的速度無關的(3.12)式。

(26)

第四章統御方程式的離散 在本文流體力學計算中，使用有限體積法可求解通量的問題，將微分方程式對某個控制體積做積分，整理後可得出離散化的代數方程式。最後使用 PISO 演算法來求解流場，其 PISO 演算法於下一章討論。 4.1 離散動量方程式我們可將之前動量方程式(2.2)式進行積分後可得下式：

(

g

)

s s d vd v v v ds v ds pd dt

∫∫∫

_∀ ρ ∀ +

∫∫

ρ − i =

∫∫

μ∇ i −

∫∫∫

_∀ ∇ ∀ (4.1) (非穩態項) (對流項) (擴散項) (源項) 其中∀：控制體積積分範圍、s：控制表面積積分範圍。 4.1.1 非穩態項非穩態使用一階近似離散，可表示如下：

(

n o

)

p p d vd v v dt t ρ ρ ∀ Δ∀ ∀ ≈ − Δ

∫∫∫

(4.2) 其中:ΔΩ為體積變量、 Δt為時間變量。 vnp為新的速度、 o p v 為舊的速度。 4.1.2 對流項 將(4.1)式中的對流項，經過高斯散度定理轉換成面積分，再經離散後整理得下式:

(

)

r g _f f v v v ds m v ρ − ≈

∑

∫∫

i c f F =

∑

(4.3)

(27)

其中s代表面法向量、m_fr為面上的質量流率、F_fc乃對流通量、下標 f 為控制體積中任一面所代表的點，對應於其面上的值。並且定義 r f m 如下： r f m =vf =w vnb+ −(1 w v) p (4.4) 因為保持空間的守恆，所以移動網格速度vg所產生的相對流量將代入 (3.7)式來計算，則表示後如下: r f f m =m f t ρ Δ∀ − Δ (4.5) 求解對流項採用一階上風法與中央差分法的線性組合，則可表示為：

( )

n

(

)

o c UD CD UD f f f f F = F +⎡_⎣γ F −F ⎤_⎦ (4.6) 其中γ 的數值介於 0 到 1 之間的因子，當γ 為0為一階上風法、γ為1中央差分法。而 UD f F 為一階上風差分的對流通量、F_fCD為中央差分法的對流通量。接著中央差分法，vf =w vnb+ −(1 w v) p ，此式vf 為相鄰格點的內差值，其中ω為權重係數、p_{為主格點、下標}nb代表為位置p 的相鄰格點。將上述代入(4.6)整理可得

(

)

(

)

max , 0 max , 0 { (1 ) max( , 0) max( , 0) } n n c p nb f f f o o o o p nb p nb f f f F m v m v m w v v m v m v γ ω = + − + ⎡ ₋ ₊ ⎤₋ ₋ ₋ ⎣ ⎦ (4.7) 4.1.3 擴散項 上述中(4.1)式中的擴散項經高斯散度定理轉換為面積分，離散後表示如下：

(28)

(

)

= _fD f f f s v ds v s F μ∇ ≈

∑

μ∇

∑

∫∫

i i (4.8) 其中μ_f為黏滯係數、 D f F 為擴散通量。在處理非結構性網格時，我們在擴散項上利用 Over-Relax 方法去近似原來的面法向量，則可表示為:

(

)

f f s = +d s −d (4.9) 其中:d為主網格中心至鄰近網格中心的向量，為正交向量；且等號右邊第二項為非正交向量。其d定義如下： 2 f pnb pnb f s d s δ δ = i (4.10) 其中δpnb為由主網格中心點與鄰近網格中心點連線所構成之向量。。則(4.11)式可以整理成如下所示：

(

)

f 2

(

n n

) ( )

o

(

)

D nb p f _f f f pnb f s F v s v v v s d s μ μ μ δ = ∇ i = − + ∇ i − i (4.11) 其上標ｏ代表前一個時間的已知數，最後將擴散項已知數移動至源項。 4.1.4 源項 源項中的將壓力梯度為已知值，故源項即為壓力項，則離散如下：

( )

pd p ∀ −

_∫∫∫

∇ ∀ ≈ − ∇ ∀ (4.12) 由於壓力梯度計算乃採取前一個時間下的計算值，將本項全部移至源

(29)

項處理。 4.2 代數方程式使離散後的各項代回動量方程式(4.1)式可得到一線性代數的方程式，在離散的過程裡代數方程式內具有特定的幾項是因為採用上一次時間下的計算結果給定，則求解的代數方程式表示如下： p nb p nb f A v =

∑

A v +Q (4.13)

(

)

2 max , 0 f p f pnb f f s A m t s μ ρ δ ⎡ ⎤ Δ∀ _⎢ _⎥ = + + ⎢ ⎥ Δ ⎣ ⎦

∑

_i (4.14)

(

)

2 min , 0 f nb f pnb f s A m s μ δ = − i (4.15) 其中:

(

)

( )

(

)

{

_CD _UD o o

}

o

( )

p f f f f f Q F F v s d v p t ρ γ μ Δ∀ ⎡ ⎤ = −_⎣ − _⎦ + ∇ − + − ∇ ∀ Δ

∑

i (4.16) (4.16)式為動量方程式離散後的已知項組成的源項。且由(4.13)式即為數值解中所要求解的矩陣。

(30)

第五章數值方法 本章以數值方法中的 PISO 演算法來處理壓力與速度項的耦合關係。且在滿足連續方程式的質量守衡情形下，使用的 PISO 演算法來求解。 5.1 預測步驟 首先利用壓力 o p ，再由(4.14)式求出速度場v*： * * ( o ) p nb p nb p f A v =

∑

A v + − ∇ Δ∀S p (5.1) 其中v*為未修正之速度場 o p 為未修正之壓力場但此時的速度場v*並未滿足連續方程式，因此進行下列修正步驟中使其來修正壓力和速度值進而滿足連續方程式。 5.2 第一次修正步驟 將(5.1)式移項整理後，則可表示如下： v*p

( )

* _o p _p p H p A Δ∀ = − ∇ (5.2) 同樣地控制體積面上的速度亦可表示為：

( )

* * _o f f _f p _f v H p A ⎛ _∀ ⎞ = −⎜_⎜ ⎟_⎟ ∇ ⎝ ⎠ (5.3) 上式中 *f f

( )

_f p _f H v p A ⎛_Δ∀⎞ = +⎜_⎜ ⎟_⎟ ∇ ⎝ ⎠ ，則代回可得速度與壓力關係式：

(31)

( )

* _o f f _f _f p _f p _f v v p p A A ⎡ _⎛_Δ∀_⎞ ⎤ _⎛_Δ∀_⎞ ⎢ ⎥ = +⎜_⎜ ⎟_⎟ ∇ −⎜_⎜ ⎟_⎟ ∇ ⎢ _⎝ _⎠ ⎥ _⎝ _⎠ ⎣ ⎦ (5.4) 其中: ∇p_f 為主格點的內插值 vf 為鄰近格點的內插值：上標” ”為主格點p與控容面 f 相鄰格點內插，如下表示

(

1

)

f p nb p p p w p w p ∇ = ∇ + − ∇ (5.6)

(

1

)

f _{p nb} _p _p v =w v + −w v (5.7) 上式中wp為權重因子。未修正的面上質量流率如下所示： * * f f f f m =ρ v is

(

)

2 * * f f f pnb f f nb p f pnb p _f f s v s p p p A s ρ ρ δ δ ⎛_Δ∀⎞ ⎡ ⎤ = − ⎜_⎜ ⎟_⎟ _⎣ − − ∇ _⎦ ⎝ ⎠ i i i (5.8) 本次修正所產生之速度場用符號v**來表示之並滿足離散後的連續方程式，而對應之壓力場用符號 ** p 來表示，其修正後的速度與壓力關係如下：

( )

** * _* p p _p p v H p A Δ∀ = − ∇ (5.9) 其中v**p 為第一次修正後的速度、p*p為第一次修正後的壓力。將(5.9)式與(5.2)式相減後，可得速度和壓力修正的關係如下：

( )

' _' p _p p v p A −Δ∀ = ∇ (5.10)

(32)

同樣的，控制體積表面上的速度v'f 亦可表式為：

( )

' _' f _f p _f v p A ⎛_−Δ∀⎞ =⎜_⎜ ⎟_⎟ ∇ ⎝ ⎠ (5.11) 上式中，v'p =v**p −v*p、p'p = p*p −pp。則第一次修正的質量流率可表示如下： ** * ' f f f m =m +m (5.12) 其中 _' ' f f f f m =ρ v ⋅s 。修正後的質量流率如下：

(

)

(

)

2 ** * f ' ' ' f f f f nb p f f pnb p _f f p _f s m m p p p s d A s A ρ ρ δ ⎛_Δ∀⎞ ⎛_Δ∀⎞ = − ⎜_⎜ ⎟_⎟ − − ⎜_⎜ ⎟_⎟ ∇ ⋅ − ⋅ ⎝ ⎠ ⎝ ⎠ (5.13) 修正後的面上質量流率滿足質量守恆定理： ** 0 f f m =

∑

(5.14) 因此將(5.13)式代入上式可整理得： ' ' ' p p nb nb f a p =

∑

a p +Q (5.15) 2 f p f pnb f p f s a A s ρ δ ⎛_Δ∀⎞ = ⎜_⎜ ⎟_⎟ ⎝ ⎠

∑

_i (5.16) 2 f nb f pnb p f s a A s ρ δ ⎛_Δ∀⎞ = ⎜_⎜ ⎟_⎟ ⎝ ⎠ i (5.17)

(

)

* ' ' _f _f _f f f f p _f Q m p s d A ρ ⎛Δ∀⎞ = − + ⎜_⎜ ⎟_⎟ ∇ ⋅ − ⎝ ⎠

∑

(5.18) 由方程式(5.15)式獲得第一次的壓力修正因子的壓力場 ' p。

(33)

5.3 第二次修正步驟 此步驟之速度v***滿足下列式子：

(

)

*** ** _** p p _p p v H p A Δ∀ = − ∇ (5.19) 將(5.19)式與(5.9)式相減可得：

( )

" ' _" p p _p p v H p A Δ∀ = − ∇ (5.20) 其中，v"p =v***p −v**p、p"p = p**p − p*p。流量修正表示如下：

(

)

₍

₎

' 2 " " " " nb f nb f _nb _p _f f f f f f f pnb p _f f p _f p A v s p p m p s d s A s A A ρ ρ ρ δ ⎛ ⎞ − ⎛_Δ∀⎞ ⎛_Δ∀⎞ ⎜ ⎟ = − ⎜_⎜ ⎟_⎟ − ⎜_⎜ ⎟_⎟ ∇ ⋅ − +⎜ ⎟⋅ ⋅ _⎜ _⎟ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠

∑

(5.21) 第二次的流量修正如下： *** ** " f f f m =m +m (5.22) 上式滿足質量守恆定理，故可得： ** " 0 f f f f m + m =

∑

(5.23) 由於 ** f f m

∑

早已滿足了質量守恆，故 ** 0 f f m =

∑

，因此上式可化為： " 0 f f m =

∑

(5.24) 則由(5.21)式代入上式可整理為： " " " " " Q p p nb nb f a p =

∑

a p + (5.25) 2 " f p f pnb f p s a A s ρ δ Δ∀ =

∑

i (5.26)

(34)

2 " f nb f pnb p s a A s ρ δ Δ∀ = i (5.27) (5.25)式即第二次壓力修正方程式，則第二次修正後之壓力場 *** p 即求，同時速度亦可由(5.19)式求得。其中(5.25)式中等號右邊第二項B 又稱為源項(source term)，其乃由不遵守質量守衡的速度v'所產生，數學式表示如下：

( )

(

)

' " 1 " Q _f _nb nb f _f f f p f f p A v s p s d A A ρ ρ ⎛ ⎞ ⎡ ∀ ⎤ = − ⎜_⎜ ⎟_⎟ + ⎢ ∇ − ⎥ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦

∑

i

∑

i (5.28) 5.4 入出口壓力邊界條件的流量計算為了求解壓力邊界上的流量，需先將進出口邊界處設定壓力值，再利用壓力值來修正出入口的流量，在此分為三個步驟。步驟一:首先將壓力邊界相鄰的網格點上，設定為欲給定之壓力值。步驟二:求解兩次壓力修正方程式時，將壓力邊界相鄰的網格點上之壓力修正值 ' p 設為零。因此壓力修正方程式的係數和源項表示如下: ' ' ' p p nb nb f a p =

∑

a p +Q 其中a_nb =0、Q'=0 其中pnb為邊界壓力。步驟三:求解完壓力修正方程式，接著讓壓力邊界相鄰的網格滿足質量守恆。

(35)

0 f f m =

∑

(5.29) 5.5 入出口壓力邊界條件的速度計算壓力邊界出/口端的速度則是透過對流邊界來離散，其中數學式如下： c 0 v v v t x ∂ ₊ ∂ ₌ ∂ ∂ 離散後如下表示：

(

) (

)

1 0 n o c n n nb nb nb c v v v v v t − + x − = Δ Δ (5.30) 則上標n、o代表新、舊時間和下標nb、c代表邊界和其相鄰的網格，則整理後： 1 o n n b c b v Crv v Cr + = + (5.31) 其中 vc t Cr x Δ = Δ 。最後透過 5.2 節計算得出的壓力邊界上的流量來修正邊界上的速度: = ( )( ) f b _b _b b x b y m v u i v j u S v S ρ + + 5.6 PISO 演算法流程 1. 設定初始值。 2. 計算新的網格與更新網格中所有資料。 3. 求解動量方程式，得到速度。 4. 求解第一次壓力修正得到 p'，並得到滿足質量守衡的速度場與壓力場。

(36)

5. 第二次的壓力修正得到p'' ，並對壓力及速度再做一次修正。 6. 上敘步驟 2 至步驟 5 為一個完整的 Time Step 的流程。持續步驟達

(37)

第六章結果與討論本章首先進行網格測試，藉由網格測試挑選出最適合的網格密度，接著利用數值方法模擬對稱性與非對稱性蠕動式管道流場，透過影響蠕動式流場的五個主要無因次參數:振幅比ε= /a h、波數α= /h λ、雷諾數 Re ρch μ = 、背壓P = b₂ b p c ρ 、相位差φ。且文中敘述相位差即是上管壁與下管壁相位角差值，若相位差為零(φ=0)，則稱此管壁為對稱性蠕動管道；若相位差為φ π= ，則蠕動式管壁截面積為定值。再將其結果進行分析，經由本研究中的數值方法計算出分析解與 Shapiro et al.[1] 所提出潤滑模型中的理論解做比較。以及觀察蠕動管道內流動情形，更進一步了解在管壁內的質點運動狀況。 6.1 網格測試本節使用管道高度為 2、長度為 100，而蠕動方式在一個週期下的平均流量當作參考；來進行四種不同網格密度 10x100、20x200、 40x400、60x600 精準度分析。並對應網格密度可知分別為Δ = 0.1x 、 0.5、0.25、0.16，Δ =y 0.2、0.1、0.05、0.03和時間間距Δt為10−3，在上/下蠕動壁面無相位差φ，且管道兩端無壓力差P =0_b 、波數α =0.1 時，測試不同振幅比ε與雷諾數Re下造成平均流量的影響，則測試結果如表(6.1)所示；其以網格密度 60x600 的平均流量作為參考，則誤

(38)

差百分比如表(6.2)所示。由表(6.1)與表(6.2)的結果得知，低雷諾數下網格密度 20x200、40x400、60x600 之平均流量差異甚小；但較高雷諾數下網格密度 20x200 和 60x600 流量超過 5%的誤差百分比，因此在本章之後的計算上，皆使用網格密度 40x400 來做測試。 6.2 週期性穩態之時間測試如表 (6.3) 所示，在背壓為零的情形下，考慮不同波數 =0.1 0.3125 α 、，設定振幅比ε=0.6、相位差φ=π/2和雷諾數 Re=100 之出口平均流量。從表(6.3)可知，在考慮高振幅比、高雷諾數與高波數的情況下，蠕動管道的平均流量約在第二十五個週期內達到穩定的週期性平均流量，因此下列於高波數的平均流量皆使用第二十五個週期；而考慮低波數時，平均流量約莫於第二十個週期達到穩定。 6.3 潤滑模型的理論解此節中利用 Shapiro et al.[1]所提出潤滑模型，進而得到一組的理論解，且理論解的詳盡的描述於附錄(A)，此理論解為無因次的平均流量和背壓P_b、振幅比ε 、雷諾數Re、相位差φ的關係式: 2 5 2 2 * 2 2 Re P [4 2 (1 cos )] 4[4 2 (1 cos )] 2 [8 2 (1 cos )] 6 [8 2 (1 cos )] b Q L ε φ ε φ ε φ ε φ − + − + = − − + + + + (6.1) 其中L代表無因次管道長，當蠕動為具對稱性完全收縮(ε=1)管道時，

(39)

有最大平均流量Q*max =2;而當背壓P =0b 時，平均流量 * Q 只為振幅比與相位差的函數，即可表示為: 2 * 2 4[4 2 (1 cos )] 2 [8 2 (1 cos )] Q ε φ ε φ − + = − + + (6.2) 當上式中振幅比ε =0為零時，此平均流量為零(Q*=0)。 6.4 數值方法之計算結果在本節中藉由有限體積法的計算結果，考慮五種不同的參數振幅比（範圍ε =0~1）、波數（範圍α =0.1~0.3125）、雷諾數（範圍R e=1~100）、相位差（範圍φ=0~π）、背壓（範圍P_b=-1.0~1.0）之設定，在其參數下造成平均流量Q*的影響，並與 6.2 節所提出的理論解進行比較。相位差對平均流量的影響 (a)改變雷諾數由圖(6.1)所示，設定波數α=0.1、振幅比ε=0.6 與背壓P =0_b ，在三種不同雷諾數下，相位差與平均流量的變化。從圖中可觀察出在相位差φ=0有最大的平均流量，隨著相位差增加則平均流量減少，當到達相位差φ π= 流量最小。此時雷諾數對平均流量的影響不大，當相位差較小時，雷諾數 Re=100 之平均流量略低於雷諾數 Re=1;當相位差較大時，則雷諾數 Re=100 略高於雷諾數 Re=1。並於理論解 6-2 式

(40)

可得知在背壓為零時，雷諾數的影響與平均流量無關，但圖中可發現與雷諾數略有影響，此是因為理論解中假設波數幾近於零的關係。 (b)改變波數由圖(6.2)考慮在波數α =0.1與波數α =0.3125下，固定振幅比 = ε 0.6、雷諾數 Re=100 與背壓P =0_b 的情形下，測試在相位差對平均流量的影響。從圖中可觀察在波數α =0.3125的平均流量明顯高於波數 0.1 α = ，且隨著相位差改變時平均流量減少幅度較小。由圖中觀察理論解在相位差φ π= 平均流量為零；但在數值結果上，相位差φ π= 在波數α =0.1的平均流量略高於零，隨著波數增加至0.3125時，有更大的平均流量。此理論解的結論是因為解假設波數幾近於零，因此數值分析結果可適用於波數不為零的假設上。 (c)改變背壓由圖(6.3)中，固定波數α =0.1、振幅比ε 0.6、= 雷諾數 Re=100，測試在不同背壓情形下，相位差對平均流量的影響。圖中可知當背壓 P_b=0的平均流量高於背壓P_b=1.0。當背壓P_b=1.0會產生回流的現象，此時平均流量為負，隨著相位差增加則流體徑向速度增加但軸向速度卻相對減少，因此回流現象較易出現在相位差較大時。當相位差φ= 0.649π 平均流量為零;而理論解為平均流量為零處發生於相位差φ= 0.554π。並於相位差φ=0、π/2、π在不同背壓下的流線圖，如圖 6.4 至圖 6.9 所示。

(41)

雷諾數對平均流量的影響 (a)改變相位差由圖(6.10)所示，設定波數α=0.1、振幅比ε=0.3 與背壓P =0_b ，在三種不同相位差的管道下，雷諾數與平均流量的變化。從圖中可發現數值分析結果隨著雷諾數的上升，平均流量的影響卻很小，在此情形下平均流量幾乎為一直線；在此觀察出當理論解背壓為零時，平均流量為一定值與雷諾數無關。然而數值之結果中，當相位差φ=0時，平均流量大於其他兩者相位差，且隨著雷諾數增加，平均流量稍微減少，此是因為在管壁震盪頻率較小下，雷諾數增加則黏滯性減少，流體在流經收縮管壁時，速度略快於低雷諾數則造成管壁的徑向推力減少，因而影響流體的軸向速度變慢。 (b)改變波數由圖(6.11)所示，固定相位差φ=π/2、振幅比ε=0.3 與背壓 P =0_b ，考慮兩種不同波數下，探討隨著雷諾數改變與平均流量的影響。從圖中可觀察出在波數α=0.1 時，平均流量的改變與雷諾數無關，與理論解結果相同；但在波數α=0.3125 時，因為管壁震盪頻率較高，因此與波數α=0.1 相比下平均流量較高，隨著雷諾數上升，黏滯性越低，造成管道的平均流量呈線性上升。因此可知在 Shapiro et al.[1]所提出潤滑模型適用於波數較低的情形下。

(42)

(c)改變背壓由圖(6.12)所示，固定相位差φ=π/2、波數α=0.1 與振幅比ε =0.3，考慮三種不同背壓下，探討隨著雷諾數改變與平均流量的影響。從圖中可知，當背壓正值時，壓力會阻止流體向出口前進，當雷諾數越大時，黏滯性越弱，則流體越容易受到壓力小的方向移動，因此當背壓P =1.0_b 時，雷諾數增加則平均流量呈線性減少。當背壓P =-1.0_b 時，壓力會造成流體往出口方向前進，因此當雷諾數越大則平均流量越大。最後探討背壓為零時，平均流量幾乎不受雷諾數的影響，此結果也符合 Shapiro et al.[1]的理論解。振幅比對平均流量的影響 (a)改變雷諾數如圖(6.13)所示，在固定相位差φ=π/2、波數α =0.1與背壓P_b=0 下，考慮兩種不同雷諾數，探討隨著振幅比改變對平均流量的影響。從圖中可觀察出隨著振幅比增加則平均流量變大，在較小的振幅比(ε =0.3)時，不同雷諾數的平均流量差異很小，當振幅比越大時則平均流量在雷諾數 Re=100 時比雷諾數 Re=1 大。 (b)改變波數如圖(6.14)所示，在固定相位差φ=π/2、雷諾數 Re=100 與背壓P_b=0 下，考慮兩種不同波數，探討隨著振幅比改變對平均流量的影響。由

(43)

圖中波數較高的曲線，隨著一開始振幅比上升，流量有明顯的提高，但之後振幅比在提高，流量的上升趨於平緩，而波數較低的曲線，隨著振幅比提高，曲線斜率緩緩增加，並推測振幅比在提高，最終兩條曲線會交會。 (c)改變背壓如圖(6.15)所示，在固定相位差φ=π/2、波數α =0.1與雷諾數 Re=100下，考慮三種不同背壓，探討隨著振幅比改變對平均流量的影響。當背壓P =-1.0_b 時，在較小的振幅比情形，可發現流體在壓力的幫助下，直接往出口流動;當振幅比逐漸變大的時候，因管壁的變形較大，促使管壁給流體的徑向推力較大則提升流體的軸向速度，近而使平均流量增加。而背壓P =1.0_b 時，壓力在流場中為阻力，因此在管道中阻止流體向出口流動則平均流量減少。而圖中可推測，當振幅比越來越大，三種不同背壓的曲線慢慢地接近，推測在振幅比ε＝1 時，流場將不受背壓的影響。並於振幅比ε=0.1、0.8在不同背壓下的流線圖，如圖 6.16 至圖 6.17 所示。波數對平均流量的影響 (a) 改變雷諾數如圖(6.18)所示，在固定相位差φ=π/2、振幅比ε=0.3 與背壓

(44)

P =0_b 下，考慮三種不同雷諾數，探討隨著波數改變對平均流量的影響。可發現在波數 0.1 時，三種不雷諾數幾乎在同一點，隨著波數增加則平均流量較大，且雷諾數高者明顯比雷諾低平均流量大，即可發現在波數較高時，雷諾數對流量的影響越來越大；此是當波數越大時，管壁震盪的頻率變高，則使管道流體加速造成流量增加。並由理論解可發現隨著波數增加，雷諾數與平均流量無關。 (b) 改變振幅比如圖(6.19)所示，在固定相位差φ=π/2、雷諾數 Re=100 與背壓P =0_b 下，考慮兩種不同振幅比，探討隨著波數改變對平均流量的影響。隨著波數增加下，振幅比越大則平均流量越高，此是當振幅比增加時則管壁上的波振幅越大，使管道壓力變化變大，近而使流體加速造成平均流量增加。並發現理論解下固定振幅比，波數的增加與平均流量無關，且振幅比越大平均流量越大。 (c) 改變背壓如圖(6.20)所示，在固定相位差φ=π/2、振幅比ε=0.3 與雷諾數 Re=100 下，考慮三種不同背壓，探討隨著波數改變對平均流量的影響。由圖中可發現三條曲線隨著波數的增加，明顯使平均流量也跟著增加，且彼此曲線的斜率相當接近。並於波數α =0.1、0.3125在不同背壓下的流線圖，如圖 6.21 至圖 6.22 所示。

(45)

背壓對平均流量的影響 (a)改變相位差如圖(6.23)所示，在固定波數α =0.1、振幅比ε=0.3 與雷諾數 Re=100，考慮三種不同相位差下，探討隨著背壓改變對平均流量的影響。從圖中可發現在固定相位差時，背壓與平均流量的關係呈線性關係，隨著背壓增加則平均流量減少，並發現具較小相位差的管道中，其曲線較為其他平緩，隨著相位差增加則曲線越傾斜。 (b)改變波數如圖(6.24)所示，在固定相位差φ=π/2、振幅比ε=0.3 與雷諾數 Re=100，考慮兩種不同波數下，探討隨著背壓改變對平均流量的影響。從圖中發現平均流量和背壓的曲線一開始幾乎呈線性下降，當背壓增加時曲線斜率的會緩緩的改變，並發現波數越大斜率改變越明顯。即可推測當波數越小時，則曲線越接近一條斜直線，則如同理論解的曲線。 (c)改變振幅比如圖(6.25)所示，在固定相位差φ=π/2、波數α =0.1與雷諾數 Re=100，考慮兩種不同振幅比下，探討隨著背壓改變對平均流量的影響。從圖中可發現在固定振幅比時，平均流量與背壓為線性的關係，隨著背壓增加則平均流量下降。當振幅比提高則曲線的斜率將改變，

(46)

在較高振幅比下斜率較趨於水平。然而平均流量為零時，即可將 6-1 式改寫成: 2 2 5 2 72 (1 cos ) P Re[4 2 (1 cos )] b Lε φ ε φ + = − + 經由上式，可知在振幅比ε=0.3下，理論解平均流量為零時，背壓為 0.227，而數值結果為 0.296。 6.5 流場結構的影響本節探討流場中質點移動的路徑和一週期內的速度分佈與沿著管道中心線的瞬時壓力。 6.5.1 徑線（pathline）與速度分佈透過徑線（pathline）可以觀察出流場中的質點運動情況，並比較管道中的某一截面速度分佈和流線。 a.不同相位差如圖 6.26 至圖 6.29 所示，在固定波數α =0.1、振幅比ε=0.6、背壓P =0_b 、1.0下，考慮不同相位差φ=0、π/2、π，探討管道中質點徑線與速度分佈的情形。在不同相位差時，以管道中心線為基準，上方的流體質點以順時針方向運動;下方的流體質點以逆時針方向運動，且流體質點重複著相同移動軌跡往入\(出)口前進。並由速度分布圖發現於管道收縮時，管道中截面積變化，當管道截面積收縮至最小時，有最大的負軸向速度。當相位差φ=0 時，流線圖(圖 6.4)與徑

(47)

線圖(6.26a)所示，則可發現上/下的流體質點移動路徑以中心線為基準呈對稱性，並可觀察出中心處為直線方向往出口前進，並從速度分佈（圖 6.27a）可發現當管壁截面積最小（上/下管壁同時達到波谷）時，速度有回流的趨勢。當相位差φ=π時，流線圖(6.6)、徑線圖 (6.26c)與速度分佈圖(6.27c)所示，由流線圖和徑線圖發現當管道舒縮時，管道截面積未改變，管道內流體幾乎有上/下移動方向，並從速度分佈得知軸向速度極小。流線圖(6.9)與徑線(6.28c)所示，當背壓P =1.0_b 和相位差φ=π時，由圖中觀察逆壓情形下，流體有回流現象產生，並對照流線圖與速度分布圖發現只有在接近波峰處有軸向速度分佈。 b.不同振幅比如圖 6.30 與圖 6.31 所示，在固定相位差φ=π/2、背壓P =1.0_b 、波數α =0.1下，考慮不同振幅比ε=0.１、0.8，探討管道中質點徑線與速度分佈的情形。由圖中可發現當振幅比較小時，因流場為逆壓阻止流體往前，則發現流場呈回流現象，並對照其流線圖(6.16)。當振幅比較大時，則流場因管壁振幅較大，雖然逆壓阻止流體前進的，但因管壁變形較大，而推向流體往前的速度較大，因此流體流向出口，並對照其流線圖(6.17)。 c.不同波數

(48)

如圖 6.32 與圖 6.33 所示，在固定相位差φ=π/2、背壓P =1.0_b 、振幅比ε=0.１下，考慮不同波數α =0.1、0.3125，探討管道中質點徑線與速度分佈的情形。即發現當波數越大時，因為壁面振盪的頻率較高，因此讓流體往前。但於靠近中心線處受逆壓影響而有回流的現象。 6.5.2 沿著管道中心線之瞬時壓力與平均壓力分布由圖(6.34)至圖(6.36)所示，可知沿著中心線的壓力分佈，當管道截面積最大時，則壓力分佈有局部最大值，當管道截面積最小時有壓力局部最小，並發現壓力隨著管道位置逐漸上升至 1.0。在不同相位時如(圖 6.34)，於相位差φ=π壓力振盪頻率較高，但振幅卻相對小於其他相位差者。在不同振幅比時的壓力分佈於圖(6.35)所示，可發現當振幅比越大時，則壓力分佈上的振幅越大。於圖(6.36)所示在不同波數下，因為波數增加，波長縮短，因此震盪頻率增加，藉此發現波數較大則震盪頻率較多。

(49)

第七章結論 1.平均流量隨著相位差呈遞減關係，當相位差越大則流量越小。在背壓為零時，固定相位差下，不同雷諾數對平均流量的影響甚小，而不同於理論解中雷諾數與平均流量無關的結果。 2.設置背壓為零，在低波數的情況下，隨著雷諾數提升則平均流量卻無改變，其與理論解之結果相同;但在高波數的情況下，隨著雷諾數提升則平均流量有明顯的上升。然而設置逆壓的情形，雷諾數與平均流量為線性遞減；若順壓的情形，雷諾數與平均流量為線性遞增的關係。 3.平均流量隨著振幅比呈非線性曲線上升，振幅比越大則流量越高; 另外在背壓為零時，並且較小振幅比的情況下，雷諾數幾乎無影響流量，但在較高振幅比的情況下雷諾數越大流量越高。 4.平均流量會隨著波數增加，並且在隨著波數增加下，雷諾數逐漸提高對流量的影響。 5.平均流量與背壓呈線性遞減的關係，並且在提高波數與振幅比的情況下，皆無法改善流量隨著背壓遞減的速率。 6.透過振幅比波數皆可以提升平均流量，然而雷諾數卻要在較高波數的情形下才可以明顯提升平均流量。其中振幅比相對於波數下，更可以大幅提升平均流量；隨著背壓提升下，流量下降速率的情形

(50)

並未改變，更得知當波數越小其數值結果越接近理論解。 7.由徑線可發現，以管道中心為基準上方質點皆以順時針方向移動，下方皆以逆時針方向移動到入/出口。在設置背壓為零下，相位差為零的情形下流體質點上方與下方運動情形具對稱性;而當相位差為π的情形下，流體質點幾乎為上下來回運動。 8.透過速度分布可以發現當管道收縮置截面積最小時則流體會有回流的趨勢。並且於壓力分布可知管道於舒縮時會造成截面積的改變，且當管道舒張至最大時有最大的局部壓力，而收縮至截面積最小時有最小的局部壓力。

(51)

參考文獻

[1] Shapiro AH, Jaffrin MY, Weinberg SL.”Peristaltic pumping with long wavelengths at low Reynolds number”, Journal of Fluid Mechanics , 37(4), pp.799‐825, 1969. [2] Jiménez‐Lozano Joel, Sen Mihir. “Streamline topologies of two‐dimensional

peristaltic flow and their bifurcations”, Chemical Engineering and Processing: Process Intensification , 49, pp.704‐715,2010.

[3] Mishra Mishra, Ramachandra Rao Adabala. “Peristaltic transport of a Newtonian fluid in an asymmetric channel”, Zeitschrift für Angewandte Mathematik und Physik ,54,pp.532‐550, 2003.

[4] Hayat Tasawar, and Ali* Nasir. “Effect of variable viscosity on the peristaltic transport of a Newtonian fluid in an asymmetric channel”, Applied Mathematical Modelling,32,pp.761‐774 ,2008.

[5] Osnat Eytan, and David Elad.”Analysis of Intra‐uterine Fluid Motion Induced by Uterine Contractions”, Bulletin of Mathematical Biology ,61,pp.221‐238, 1999. [6] Weinberg SL, Eckstein EC, Shapiro AH. “An experimental study of peristaltic

pumping”, Journal of Fluid Mechanics , 49(3),pp.461‐479,1971.

[7] Fung YC, Yih CS. “Peristaltic transport”, ASME Journal of Applied Mechanics , 45,pp.669‐675, 1968.

[8] Yin FCP, Fung YC. “Comparison of theory and experiment in peristaltic transport”, Journal of Fluid Mechanics ,47(1),pp.93‐112, 1971.

[9] Usha Srinivasan, Ramachandra Rao Adabala.” Effects of curvature and inertia on the peristaltic transport in a two‐fluid system”, International journal of Engineering Science, 38,pp.1355‐1375, 2000.

[10] Wilson DE, Panton RL. “Peristaltic transport due to finite amplitude bending and contraction waves”, Journal of Fluid Mechanics ,90(1),pp.145‐159, 1979.

[11] Tong P, Vawter D. “An analysis of peristaltic pumping”, ASME Journal of Applied Mechanics, 39(4),pp.857‐863, 1972.

[12] Pozrikidis C. “A study of peristaltic flow”, Journal of Fluid Mechanics, 180,pp.515‐527, 1987.

[13] Takabatake S, Ayukawa K. “Numerical study of two‐dimensional peristaltic flows”, Journal of Fluid Mechanics , 122,pp.439‐465, 1982

(52)

(53)

附錄 A(潤滑理論之流量方程式) 如圖(2.1)在一蠕動式管道流場，依據潤滑模型假設管壁為無限長的波長和自由慣性（inertia-free）下，在固定座標下波速為ｃ，而在相對固定座標下，使用波座標可將流體視為穩態，且波座標與固定座標兩者關係以下表示之: x=X −ct , u= −U c (A1) y=Y , v=V_(A2) 其中u與v是在波座標下的 x 與 y 軸方向速度則在潤滑模型的假設下，Navier-Stokes Equations 在 x 軸可以簡化成下式: 2 2 dp u dx μ y ∂ = ∂ (A3) 而在管壁在波座標下之移動方程式為: 1 2 ( ) cos x y x h a π λ = + (A4) 2 2 ( ) - -a cos( x- ) y x h π φ λ = (A5) 此邊界狀況如下表示之: 1 ( ) u y= y = − , c u y( = y₂)= − c (A6) 則可獲得在 u 方向的速度解 2 1 2 1 2 1 [ ( ) ] 2 dp u y y y y y y c dx μ = − + + − (A7) 波座標下的流量為:

(54)

1 2 3 1 2 1 2 1 ( ) ( ) 12 y y dp q udy y y c y y dx μ − =

∫

= − − − (A8) 透過波座標下，可將流體視為穩態，而在固定座標下的流量為 1

(

)

2 1 2 ( ) y y Q=

∫

u+c dy= +q c y −y (A9) 而在一個週期下平均流量，如下表示之：. 0 1 2 T Q Qdt q hc T =

_∫

= + (A10) 根據上式中的（A7）即可將找出壓力梯度為: ₃ ₂ 1 2 1 2 12 ( ) ( ) dp q c dx μ y y y y ⎛ ⎞ = − _⎜ + _⎟ − − ⎝ ⎠ (A11) 而在(A11)式中，需利用下列積分公式解之.

(

)

₍

₎

2 2 3 2 0 2 2 2 cos d π θ πα α β+ θ = α −β

∫

(A11a)

(

)

₍

(

₎

)

2 2 2 3 5 2 0 2 2 2 cos d π θ π α β α β θ _α _β + = + −

∫

(A11b) 將上式中(4)、(5)的 y1與 y2合併，並對此式進行一個波長的積分，意即壓力在每單位波長上升的多寡，其式可表示如下: 2 2 2 2 5 2 2 2 3 2 6 [8 2 (1 cos )] 24 [4 2 (1 cos )] [4 2 (1 cos )] q a A ac p a A a A λ μλ φ μλ φ φ − + + Δ = − − + − + (A12) 將式(A9)合併，並引入無因次參數，則產生無因次的平均流量Q*如下 * 2 5 2 2 * 2 2 Re [4 2 (1 cos )] 4[4 2 (1 cos )] 2 [8 2 (1 cos )] 6 [8 2 (1 cos )] p Q ε φ α λ ε φ ε φ ε φ Δ − + − + = − − + + + + (A13) 其中 Q* =Q hc/ , * 2 / p_λ p_λ ρc Δ = Δ ε =a h/ α =h/λ

(55)

表 6.1(a).波數α=0.1、振幅比ε=0.1、相位差φ=0、背壓 P=0 平均流量測試網格密度 10x100 20x200 40x400 60x600 Re=1 3.157x10-2 3.150x10-2 3.143x10-2 3.140x10-2 Re=10 3.186x10-2 3.142x10-2 3.139x10-2 3.133x10-2 Re=20 3.214x10-2 3.115x10-2 3.099x10-2 3.082x10-2 Re=100 3.528x10-2 3.069x10-2 2.894x10-2 2.840x10-2 表 6.1(b).波數α=0.1、振幅比ε=0.3、相位差φ=0、背壓 P=0 平均流量測試網格密度 10x100 20x200 40x400 60x600 Re=1 2.746x10-1 2.689x10-1 2.686x10-1 2.681x10-1 Re=10 2.750x10-1 2.695x10-1 2.684x10-1 2.674x10-1 Re=20 2.800x10-1 2.688x10-1 2.660x10-1 2.630x10-1 Re=100 3.243x10-1 2.759x10-1 2.559x10-1 2.449x10-1 表 6.1(c).波數α=0.1、振幅比ε=0.6、相位差φ=0、背壓 P=0 平均流量測試網格密度 10x100 20x200 40x400 60x600 Re=1 9.380x10-1 9.273x10-1 9.251x10-1 9.239x10-1 Re=10 9.490x10-1 9.285x10-1 9.249x10-1 9.234x10-1 Re=20 9.684x10-1 9.315x10-1 9.233x10-1 9.198x10-1 Re=100 1.121 9.454x10-1 8.992x10-1 8.860x10-1

(56)

表 6.2(a).波數α=0.1、振幅比ε=0.1、相位差φ=0、背壓 P=0 之誤差百分比網格密度 10x100 20x200 40x400 Re=1 0.54 0.31 0.09 Re=10 1.69 0.28 0.19 Re=20 4.28 1.07 0.51 Re=100 24.22 8.06 1.90 表 6.2(b).波數α=0.1、振幅比ε=0.3、相位差φ=0、背壓 P=0 之誤差百分比網格密度 10x100 20x200 40x400 Re=1 0.02 0.002 0.001 Re=10 2.84 0.78 0.37 Re=20 6.46 2.20 1.14 Re=100 29.77 10.40 2.40 表 6.2(c).波數α=0.1、振幅比ε=0.6、相位差φ=0、背壓 P=0 之誤差百分比網格密度 10x100 20x200 40x400 Re=1 0.01 0.003 0.001 Re=10 2.77 0.55 0.06 Re=20 5.28 1.26 0.38 Re=100 26.52 6.70 1.48

(57)

表 6.3.a 週期性平均流量測試波數α=0.1(相位差=π/2、振幅比ε=0.6、雷諾數 Re=100、背壓 Pb=0.0) 週期 1 3 5 7 9 11 13 15 17 19 平均流量 7.527 x10-2 4.510 x10-1 5.064 x10-1 5.253 x10-1 5.318 x10-1 5.340 x10-1 5.348 x10-1 5.350 x10-1 5.351 x10-1 5.352 x10-1 週期 21 23 25 27 29 30 平均流量 5.352 x10-1 5.352 x10-1 5.352 x10-1 5.352 x10-1 5.352 x10-1 5.352 x10-1 表 6.3.b 週期性平均流量測試於波數α=0.3125(相位差=π/2、振幅比ε=0.6、雷諾數 Re=100、背壓 Pb=0.0) 週期 1 3 5 7 9 11 13 15 17 19 平均流量 9.326 x10-2 7.608 x10-1 8.634 x10-1 9.269 x10-1 9.718 x10-1 1.003 1.041 1.052 1.059 1.065 週期 21 23 25 27 29 30 平均流量 1.070 1.072 1.073 1.073 1.073 1.073

(58)

壁面蠕動具相位差之管道流場數值分析

國 立 交 通 大 學

機械工程學系

碩 士 論 文