凹槽結構管道微混合器流場分析

全文

(1)國立交通大學機械工程學系碩士論文. 凹槽結構管道微混合器流場分析 Analysis of the mixed flow in microchannels with grooved structure. 研究生：謝崇民指導教授：崔燕勇. 教授. 中華民國九十四年七月.

(2) 凹槽結構管道微混合器流場分析 Analysis of the mixed flow in microchannels with grooved structure. 研究生：謝崇民. Student：Chung-Ming Hsieh. 指導教授：崔燕勇. Advisor：Yeng-Yung Tsui. 國立交通大學機械工程學系碩士論文. A Thesis Submitted to Department of Mechanical Engineering College of Engineering National Chiao Tung University in partial Fulfillment of the Requirements for the Degree of Master of Science in Mechanical Engineering July 2005 Hsinchu, Taiwan, Republic of China. 中華民國九十四年七月.

(3) 凹槽結構管道微混合器流場分析學生：謝崇民. 指導教授：崔燕勇博士. 國立交通大學機械工程學學系碩士班. 摘. 要. 兩種流體的混合在許多的微流體系統中是必要的過程，由於微流道的尺寸非常小，因此流體以層流方式流動，混合過程是以分子擴散方式進行，為了在合理的時間內達成有效的混合，必須控制流體使使兩流體界面增加，藉此才能降低擴散路徑，加強分子擴散行為以完成整個混合過程。具凹槽結構之微管道可產生二次流場，其可延展並且折疊流體介面，降低流體間之擴散路徑。在本文中，討論幾何外型參數對直型溝槽混合器混合性能的影響，並比較人字形混合器及具傾斜阻塊混合器之混合效率，數值計算結果顯示，溝槽較深的溝槽及較寬的溝槽有較佳的混合效果，雙傾斜溝槽比單傾斜溝槽有更佳的混合品質，於 SHM 管道中置入傾斜阻塊，可藉由影響流場型態及產生橫向質傳效應方式提升混合效率，因此，在有限的尺寸下，OBSHM 概念應可提供微混合器設計有用的資訊. i.

(4) Analysis of the mixed flow in microchannels with grooved structure student：Chung-Ming Hsieh. Advisors：Dr.Yeng-Yung Tsui. Department of Mechanical Engineering National Chiao Tung University. ABSTRACT. Two-fluids mixing is an essential processes in many microfluidics systems, Due to the small dimensions of the micro channels, the flow is predominantly laminar and mixing is limited by molecular diffusion. In order to efficiently mix in a reasonable time, fluids must be controlled to increase the interfacial area between the fluids so that the diffusion path is decreased to enhance the molecular diffusion in mixing processes. micro channels with grooved structure can create the secondary flow which stretched and folded material line to reduce the diffusion path between fluid streams. in the thesis, a study of the influence of various geometrical parameters of straight ridges micromixer on mixing performance is reported, then, comparing the efficiency between the staggered herringbone micromixer(SHM) and new type micromixer which we call OBSHM. The results indicate that a large groove depth and width has better performance then others. the mixing quality of microchannls with double grooved structure is superior to the the single structure, the oblique blocks placed in the SHM can improve mixing efficiency by affecting the flow pattern and alter the flow to create lateral mass transport, therefore, in the limits of dimensions, the ideal of OBSHM can provide useful information in the design of micomixers.. ii.

(5) 誌. 謝. 誠摯感謝指導教授崔燕勇博士在這兩年求學生涯中，無論在課業及論文研究上，殷切指導與教誨，使學生得以順利完成學業，並感謝中山科學研究院第二研究所提供此次在職進修的機會，同時也要感謝吳添成學長、胡育昌學長在學業上的指導，另在口試過程中，謝謝林振德教授、傅武雄教授、楊文美教授在論文上的建議及指導，在此獻上最深的敬意。此外，特別感謝父母多年養育之恩，以及兄長、女友韶文在這段時間的不斷的支持與鼓舞，沒有他們的辛苦付出與鼓勵，就無法順利完成學業。最後感謝在我生涯當中，曾經共事的長官、同事，因為他們的教誨與惕勵，使我日益成長; 學校師長們的認真教導，更使我在學業專業領域獲益匪淺，感謝求學階段，陪伴一起走過的學長、朋友、同學及學弟們，謝謝你們！！謝崇民 2005.7.14 于交大機械系計算流體實驗室. iii.

(6) 目錄中文摘要英文摘要誌謝目錄表目錄圖目錄符號說明第一章 1.1 1.2 1.3 1.4 1.5 1.6 1.7 第二章 2.1 2.2 2.3 2.4 第三章 3.1 3.2 3.3 3.4 3.5 第四章 4.1 4.2 4.3 4.4 第五章. ………………………………………………………………… ………………………………………………………………… ………………………………………………………………… ………………………………………………………………… ………………………………………………………………… ………………………………………………………………… ………………………………………………………………… 緒論…………………………………………………………… 前言…………………………………………………………… 生物晶片之分類……………………………………………… 微流系統……………………………………………………… 混合機制……………………………………………………… 微混合器之分類……………………………………………… 文獻回顧……………………………………………………… 研究動機與方向……………………………………………… 混合器幾何外型介紹………………………………………… 直型傾斜溝槽形混合器……………………………………… 直型雙傾斜溝槽混合器……………………………………… 人字形混合器………………………………………………… 具傾斜阻塊人字形混合器…………………………………… 數學模式……………………………………………………… 基本假設……………………………………………………… 統御方程式…………………………………………………… 邊界條件……………………………………………………… 渦度計算……………………………………………………… 混合效率的定義……………………………………………… 數值方法……………………………………………………… 傳輸方程式離散化…………………………………………… 解線性代數方程式…………………………………………… 壓力速度耦合關係式………………………………………… 解題步驟……………………………………………………… 結果與討論……………………………………………………. iv. i ii iii iv vi vii x 1 1 1 2 2 3 4 10 12 12 13 15 16 18 18 18 20 21 21 23 23 26 26 30 32.

(7) 5.1 5.2 5.3 5.4 5.5 5.6 第六章參考文獻. 網格測試與實驗值比較……………………………………… SRM 混合器混合效能分析與探討…………………………… DSRM 混合器混合效能分析與探討………………………… 直型傾斜溝槽混合器最佳化外型…………………………… 人字形、具傾斜阻塊人字形混合器混合效能分析………… SRM、DSRM、SHM 及 OBSRM 流線軌跡探討…………… 結論…………………………………………………………… …………………………………………………………………. v. 32 33 37 38 39 41 42 43.

(8) 表目錄表 5.1. 中心截面在 X 方向渦度強度…………………………………. vi. 45.

(9) 圖目錄圖 1.1 圖 1.2 圖 1.3 圖 1.4 圖 1.5 圖 1.6 圖 1.7 圖 1.8 圖 1.9 圖 1.10 圖 1.11 圖 1.12 圖 3.1 圖 3.2 圖 5.1 圖 5.2 圖 5.3 圖 5.4 圖 5.5 圖 5.6 圖 5.7 圖 5.8 圖 5.9 圖 5.10 圖 5.11 圖 5.12 圖 5.13 圖 5.14 圖 5.15 圖 5.16 圖 5.17 圖 5.18 圖 5.19 圖 5.20 圖 5.21 圖 5.22 圖 5.23. 參考文獻[3]混合器結構外型與分析結果參考文獻[4]混合器外型與模擬結果參考文獻[5]混合器外型與流場、質點分佈圖參考文獻[6]混合器外型圖參考文獻[7]不同雷諾數流場分佈參考文獻[8] 混合器外型與分析結果參考文獻[9] 混合器外型與濃度分佈圖參考文獻[10]外型示意圖參考文獻[12]外型示意圖參考文獻[13]混合器外型與原理參考文獻[14]雙層微結構參考文獻[15]微噴嘴第一階段數值計算區域示意圖……………………………… 第二階段數值計算區域示意圖……………………………… SRM 截面格點數與 Stroock 實驗混合效率比較圖………… SHM 截面格點數與 Stroock 實驗混合效率比較圖………… SRM 網格結構………………………………………………… SRM 週期網格結構…………………………………………… SRM 第一週期截面位置圖…………………………………… SRM 第一週期截面速度、流線及濃度分佈圖……………… SRM 溝槽深度與混合效率關係圖…………………………… SRM 溝槽寬度與混合效率關係圖…………………………… SRM 溝槽傾斜角度與混合效率關係圖……………………… SRM 各週期截面位置示意圖………………………………… DSRM 各週期截面位置示意圖……………………………… SHM 各週期截面位置示意圖………………………………… OBSHM 各週期截面位置示意圖…………………………… SRM(Gd=0.2d)各週期截面濃度分佈圖……………………… SRM(Gd=0.42d)各週期截面濃度分佈圖…………………… SRM(Gd=0.6d)各週期截面濃度分佈圖……………………… SRM(Gd=0.2d)各週期截面質點分佈圖……………………… SRM(Gd=0.42d)各週期截面質點分佈圖…………………… SRM(Gd=0.6d)各週期截面質點分佈圖……………………… 傾斜截面示意圖……………………………………………… 傾斜截面中心點橫向速度位置示意圖……………………… 不同溝槽深度在傾斜中心截面中點橫向速度分佈圖……… Gd=0.2d中心截面A’在y-z平面速度及渦度分佈圖………… vii. 46 46 46 47 47 47 48 48 49 49 49 50 51 51 52 52 53 53 54 54 57 57 57 58 58 58 58 59 59 60 60 61 61 62 62 63 64.

(10) 圖 5.24 圖 5.25 圖 5.26 圖 5.27 圖 5.28 圖 5.29 圖 5.30 圖 5.31 圖 5.32 圖 5.33 圖 5.34 圖 5.35 圖 5.36 圖 5.37 圖 5.38 圖 5.39 圖 5.40 圖 5.41 圖 5.42 圖 5.43 圖 5.44 圖 5.45 圖 5.46 圖 5.47 圖 5.48 圖 5.49 圖 5.50 圖 5.51 圖 5.52 圖 5.53 圖 5.54 圖 5.55 圖 5.56 圖 5.57 圖 5.58 圖 5.59 圖 5.60 圖 5.61 圖 5.62. Gd=0.2d中心截面B’在y-z平面速度及渦度分佈圖………… Gd=0.42d中心截面A’在y-z平面速度及渦度分佈圖……… Gd=0.42d中心截面B’在y-z平面速度及渦度分佈圖……… Gd=0.6d中心截面A’在y-z平面速度及渦度分佈圖………… Gd=0.6d中心截面B’在y-z平面速度及渦度分佈圖………… SRM(Gw=0.25d)各週期截面濃度分佈圖…………………….. SRM(Gw=1d)各週期截面濃度分佈圖……………………… SRM(Gw=0.25d)各週期截面質點分佈圖…………………….. SRM(Gw=1d)各週期截面質點分佈圖……………………… 不同溝槽寬度在傾斜中心截面中點橫向速度分佈圖……… Gw=0.25d中心截面A’在y-z平面速度及渦度分佈圖……… Gw=0.25d中心截面B’在y-z平面速度及渦度分佈圖……… Gw=1d中心截面A’在y-z平面速度及渦度分佈圖…………. Gw=1d中心截面B’在y-z平面速度及渦度分佈圖…………... SRM(θ＝60°)各週期截面濃度分佈圖……………………….. SRM(θ＝75°)各週期截面濃度分佈圖………………………... SRM(θ＝60°)各週期截面質點分佈圖………………………... SRM(θ＝75°)各週期截面質點分佈圖……………………….. 不同溝槽傾斜角度在傾斜中心截面中點速度分佈圖……… (θ＝60°)中心截面 A’在 y-z 平面速度及渦度分佈圖……… (θ＝60°)中心截面 B’在 y-z 平面速度及渦度分佈圖……… (θ＝75°)中心截面 A’在 y-z 平面速度及渦度分佈圖……… (θ＝75°)中心截面 B’在 y-z 平面速度及渦度分佈圖……… DSRM 網格結構……………………………………………… DSRM 週期性網格結構……………………………………… DSRM 第一週期橫截面位置圖……………………………… DSRM 第一週期橫截面速度、流線及濃度分佈圖…………… DSRM 各週期截面濃度分佈圖………………………………. DSRM 各週期截面質點分佈圖………………………….…… DSRM V.S SRM 混合效率比較圖…………………………… SRM、DSRM、最佳化 SRM 混合效率比較圖……………… SHM 網格結構………………………………………………… SHM 週期性網格結構………………………………………… SHM 第一週期橫截面位置圖……………………………….. SHM 第一週期截面速度、流線及濃度分佈圖……………… SHM 各週期截面濃度分佈圖………………………………… SHM 各週期截面質點分佈圖………………………………… OBSHM 網格結構…………………………………………….. OBSHM 週期性網格結構…………………………………….. viii. 64 64 64 64 64 65 65 66 66 67 68 68 68 68 69 69 70 70 71 72 72 72 72 73 73 74 74 77 77 78 79 80 80 81 81 84 84 85 85.

(11) 圖 5.63 圖 5.64 圖 5.65 圖 5.66 圖 5.67 圖 5.68 圖 5.69 圖 5.70 圖 5.71 圖 5.72. OBSRM 第一週期截面位置圖………………………………... OBSRM 第一週期截面速度、流線及濃度分佈圖…………… 第一週期 xy 平面速度及濃度分佈圖………………………… OBSHM 各週期截面濃度分佈圖……………………………... OBSHM 各週期截面質點分佈圖……………………………... SHM V.S OBSHM 混合效率比較圖………………………… SRM 流線軌跡圖……………………………………………… DSRM 流線軌跡圖…………………………………………… SHM 流線軌跡圖……………………………………………… OBSHM 流線軌跡圖…………………………………………. ix. 86 86 89 90 90 91 92 92 92 92.

(12) 符號說明. Ai. 微小面積. Bh 阻塊高度 Bl 阻塊長度 Bw 阻塊寬度 C0 溶質初始容度 Ch 管道高度 CW 管道寬度. c c c∗. 濃度. D. 擴散係數. F. 通量. 平均濃度無因次化濃度. Gd 溝槽深度 Gp 溝槽間距 Gw 溝槽寬度 L. 管道長度. LC 週期性幾何外型長度 Li 入口至第一週期性幾何外型距離 LO 出口至第最後一週期性幾何外型距離. m. 質量流率. M 混合指數 N. 溝槽總數. P. 壓力. Pe 裴立數 Re 雷諾數 K S f 控容面之面積向量. u v JK V V. X 方向速度項 Y 方向速度項速度向量. 入口平均速度 x.

(13) w. Z 方向速度項. w f 權重因子希臘符號. Ω Ω. 渦度. 平均渦度. μ 黏滯係數 ρ 密度. Φ 變數 θ 溝槽傾斜角度 ψ 溝槽傾斜角度. Γ. 係數參數. 下標. b c f in out p. 壁面上中點控制體積周圍網格之中心格點控制體積任一面上之中點入口出口主格點. 上標. c. 對流項. d. 擴散項. CD 中央差分法 P. 週期. UD 上風差分法 *,** 預測值、修正值. xi.

(14) 第一章. 緒論. 1.1 前言近年來，微機電系統（micro electromechanical systems）發展已逐漸邁入成熟階段，藉由半導體製程技術，結合機械工程、材料工程、電子及生物科技等，使得 MEMS 有了更多元的發展與應用，其中，透過生物科技或醫學應用與 MEMS 微製造技術的結合，其所發展之領域，即所謂的生醫微機電(BioMEMS)。生醫微機電系統的發展，在近程主要是診斷技術的應用，隨著基因體時代來臨，在結構基因體、功能性基因體及蛋白質基因體的發展皆面臨高靈敏度及低成本需求，生物晶片（Biochip）即是應此需求，作為結合生醫微機電與生物分子控制而發展的技術。 1.2 生物晶片之分類廣義而言，依其特性可分為兩類：感測晶片（Sensing Chip）與處理晶片（Processing Chip）兩種。感測型晶片是以微陣列技術(microarray) 為主之晶片，包含蛋白質晶片（Protein-chip）、基因晶片（Gene-chip）與生物感測晶片（Biosensor-chip）等三種，而蛋白質晶片與基因晶片又可統稱微陣列晶片(microarray chip)，其乃將不同排列組合核甘酸短序列做成探針，探針端附有發光酵素，藉由色彩訊號以表示基因序列。處理型晶片是藉由微機電微小化技術，使生物分析儀器微小化，處理分析樣品功能尺寸介於釐米與微米間，其優點可加速反應，減少試劑與檢體用量，降低製造成本，簡單的說，就是運用縮小結構技術，將實驗室進行之實驗流程縮小至晶片上進行，故通常又稱微流體晶片、實驗室晶片（lab-on-a-chip)或微型全分析系統，此微型全分析系統概念被全世（Micro-Total-analysis-System, µ-TAS）界許多研究機構看好，紛紛投入資金研發，雖然其挑戰性高，但其所 1.

(15) 需分析樣本量少、攜帶方便及操作簡單之優點，可以廣泛應用於不具實驗室設備單位，因此，在這些誘人的因素下，加入這個研究領域的組織或學校研究機構愈來愈多，一旦技術成熟，對於未來生物分析檢測技術將貢獻良多。 1.3 微流系統（Microfluidics System）實驗室晶片可透過微小之晶片，完成傳統實驗室所需的分析流程。一般此種處理型晶片包含前處理、雜交與檢測三種過程，此過程必須藉由所謂微流系統，以達到驅動、傳輸、混合、反應、分離等目的。因此，微流系統包含了各種之微元件（Microdevices），如微幫浦（Micropump）、微流道（Microchannels）、微分離器（Microseparator）、微混合器（ Micromixers ）、微反應器（ Microreactors ）、微感應器（Microsensors）等，藉由這些元件的結合，可精確控制流體流量，具有體積小、反應時間快速及適合量產之優點，而微混合器，控制檢體與試劑混合效率良窳的關鍵元件，對於分析結果影響甚巨，因此為本篇討論之主題。 1.4 混合機制討論微混合器前，必先瞭解其混合機制與一般巨觀之差異。從巨觀角度，流體混合可以藉由雷諾數的增加所產生的紊流機制，透過渦流擴散達到混合，因此往往只需利用外力，如加入攪拌裝置或改變流道設計於流道中置入阻礙物即可產生混合，然而，當尺寸微小化後，流體於微流道中的流體特性，可能已與巨觀流體有所不同，因此無法以一般流體混合的觀念思考，在微流道中，因流體流速緩慢，流體流動完全以層流機制進行，因雷諾數可以視為流體流動慣性力與黏滯力相對大小的常數，當雷諾數很小時，整個流場將由黏滯力所主導，因此流場中無法透過紊流機制達到混合的目的，其只能完全透過分子擴 2.

(16) 散方式進行，然而，以分子擴散機制進行混合是一種非常緩慢耗時的過程，在微流體晶片中是非常不符合需求的，因此，如何設計一高效率的混合器至今仍持續被各界研究中。 1.5 微混合器之分類一般微混合器可劃分為兩類：主動式混合 (active mixing)與被動式混合（passive mixing)。主動式混合利用外加元件產生驅動力，如超音波、電磁力、壓力梯度增加混合效果;被動式混合則是藉由結構外型的變化，使流體產生移動及渾沌對流的效果。一般而言，主動式混合器的混合較佳，但是其繁雜的製程及需考量外加驅動元件與晶片之整合，反觀被動式混合器的簡單結構，成為各界設計微流系統中最佳的選擇方法之一。被動式混合器之設計約有下列幾種原理[1][2]： 1.5.1 薄層化（Lamination）根據費克第二定律(Fick’s Second Law) t=. ∂c ∂ 2c = D 2 ，進而可得知 ∂t ∂x. x2 ，混合時間（ t ）與混合路徑 x 的平方成正比關係，因此藉由降低 2D. 混合路徑可提升混合效率，依據此理論，利用薄層化或將流道分為數個次流道的技術，藉由分離（spliting）與再結合(recombining)，增加流體接觸面積，加速流體混合效果。 1.5.2 注射法（Injection）與薄層化剛好相反，其乃將流體分成好幾個次流體並將其注入另一流體當中，此種方法也可以增加接觸面積，提升混合效率。 1.5.3 混沌對流（Chaotic advection）由動態系統理論得知，當流場為二維暫態或三維(暫態或穩態) ，混沌粒子運動將會發生。因此，混沌對流是一種藉由在空間、時間完全規律之速度流場內的一種混合過程，這樣的流場使得粒子剛開始非常接近，然後將分離得非常遠，最後得到很好的混合效果。 3.

(17) 1.6 文獻回顧數值模擬與實驗部分 Stroock 等人(2002)[3]說明在具平滑壁面流道之微混合器中，因流體速度緩慢，試劑進行充分效混合非常困難，因流場中，流體混合完全以分子擴散方式進行，其混合長度（ ∇Ym ）隨 Peclet number(Pe) 呈線性成長（ ∇Ym = Pe × l ），因此設計底部具溝槽結構之流道，使流場產生與主流方向垂直的橫向混沌流場，實驗結果證明其可降低混合長度（混合長度將隨 Peclet number 的對數成線性關係 ∇Ym = ln( Pe) × l ）。並比較平滑壁面微混合器(simple microchannel Mixer)、直條形溝槽結構微混合器(striaight ridges micromixer，簡稱 SRM)、與人字形結構微混合器（staggered herringbone mixer ，簡稱 SHM).在不同 Peclet number 時混合的狀況，其結果發現：SHM 具有較佳的混合效果，且其混合長度隨 Peclet number 增加而增長，如圖 1.1 所示。 Liu 等人(2004)[4]探討兩種不同性質(水/甘油水溶液)之流體在（蜿蜒形）方波形微混合器（square wave mixer）、三維蜿蜒形微混合器（3D serpentine micromixer）與人字形微混合器（staggered herringbone micromixer）中，觀察藉由濃度的變化，對混合行為之影響。首先以數值方法比較三種混合效率，結果發現 serpentine micromixer 與 staggered herringbone micromixer 遠優於 square wave mixer，後續再比較 serpentine micromixer 與 staggered herringbone micromixer 兩種，觀察濃度變化與混合效率之關係。選擇甘油溶液濃度 ψ=0、0.2 及 0.4，雷諾數在 Re=1 及 10 的情況下，在 serpentine mixer 混合器中，低雷諾數(Re=1)時,因分子擴散作用，混合效率隨著甘油溶液濃度之增加而變差，但高雷諾數 (Re=10)時，因流體對流效應的增加，隨著雷諾數之增加引起混沌對流， 4.

(18) 混合效率隨著甘油溶液濃度之增加而變佳，但在 staggered herringbone mixer(SHM)混合器中，混合效率隨甘油溶液濃度之增加而變差，且雷諾數之變化（Re<10）對混合效果無明顯的改變，SHM 混合效果優於 serpentine micromixer。如圖 1.2 所示。 Wang 等人(2003)[5]利用 CFD 套裝軟體(MemCFD)及質點追蹤技術，分析流體進入具直型溝槽之Ｔ形微流道現象，並根據溝槽結構產生之週期性速度流場特性，模擬此週期流場狀況，並利用 Poincare map 之觀念，判斷流場中粒子之紊亂程度。其發現雷諾數(Reynold Number) Re=5 時，當溝槽振幅 α 值（aspect ratio groove）為 0.05 時，流場形態並沒有發生不規則流動的情形，當 α=0.30 時，流體型態變得較為不規則，但是並沒有明顯的證據說明此現象是 chaotic 所引起;其原因推測為流道底面非等向性之溝槽結構，使流體產生與主流方向相垂直之二次流體，增加兩流體間接觸面積，提升混合效果，並得知管道長度與 α 的指數值成正比關係，與流量或雷諾數大小無關，如圖 1.3 所示。 Gobby 等人(2003)[6]利用計算流體力學數值模擬方法，研究在氣體在 T 形微混合器中混合特性，其假設流體為層流、等溫，入口分別為甲醇與氧兩種氣體，考慮氣體的可壓縮性及滑移邊界性質，進口溫度為 400K，藉由流速、高寬比、進口角度及節流孔變化，探討對混合之影響。模擬結果顯示，混合長度隨著流體速度增加而增加，並受到混合器高寬比之影響，在流道寬度不變時，高寬比增加時，混合長度與其成反比;進口角度的改變方面，對混合效率與壓力變化沒有特別的影響;利用節流閥設計概念，能夠有效增加混合效果，當其節流孔越小，混合效果也越好，但相對在擴張區有較高之壓力降。如圖 1.4 所示。 Engler 等人(2004)[7]藉由數值模擬與實驗兩種方法，探討僅靠對流 5.

(19) 效應進行混合之微混合器中，其渦流現象對混合效果之影響。在 T 形靜微混合器（static micromixer）中，隨著雷諾數增加，產生渦度依大小可分成三種型態，由小到大分別為 stratified flow regime,vortex flow regime 與 engulfment flow 等三種，而 engulfment flow 產生之非對稱流場正是導致混合效能提升的原因。雷諾數趨近 200 時 engulfment flow 將產生，如圖 1.5 所示。 Liu 等人(2000)[8]設計一三維 S(serpentine)形流道，觀察其雷諾數在 6-70 之間流場混合情形，發現混合機制除擴散外，還包括因幾何形狀產生的渾沌對流，加強了混合的效果，混合效率隨雷諾數增加而提升，同時並與另方波管及直形管做比較，混合性能優於兩者。另發現直行管反而隨雷諾數增加，混合效果反而降低的現象。如圖 1.6 所示。 Schönfeld(2004)[9]等人利用數值計算軟體分析具有淺浮雕溝槽結構管道及方形截面彎管之混合效果。本文分淺浮雕溝槽結構及方形截面彎管結構兩部分討論。在淺浮雕溝槽結構方面，比較具有單壁面及雙壁面淺浮雕溝槽結構之混合效率，結果發現具雙壁面淺浮雕結構管道具有較佳的混合效果，且雷諾數的大小對相對橫向速度的影響非常微小。另外在方形截面彎管模擬方面，因離心力作用使的管道有二次流產生，雷諾數越高，混合效率越好。其重要之無因數參數共有兩項： Re(Reynold Number)與 k (Dean Number（ k = Re ( d R) 2 )；臨界值 kc =200， 1. 當 k < kc ，二次流由兩個旋轉方向相反之渦流構成，當 k > kc ，靠近管道外壁區域將再產生另一組渦流。 k =300～400 時，渦度有明顯的增強。並比較淺浮雕溝槽結構及方形截面結構彎管之混合效率，藉由相對橫向速度之大小可以發現，彎曲方形截面結構管道具較佳之混合效果。如圖 1.7 所示。 6.

(20) Kim 等人(2004)[10]提出新的被動式混合器（barrier embrdded chaotic micromixer）概念，簡稱 BEM，參考圖 1.8。藉由在具傾斜溝槽流道中，週期性置入阻礙物，使其產生交替作用之速度場，造成流體產生延伸與折疊，以增加混合效果。利用實驗與數值模擬方法，與平滑壁面及底部傾斜溝槽 slanted groove（SGM）等兩種微混合器做比較，明顯發現三者之混合效率（強度）皆隨著雷諾數增加而下降；並定義所需的混合長度，BEM 之混合長度隨雷數增加成對數成長，與 Stroock 等人得結果相符合。而 BEM 的混合長度約比 SGM 短近 2 倍；BEM 所需之混合長度分別在最低及最高流率時，比平滑壁面流道少約 4～10 倍左右，由此可知 BEM 提供較佳之混合結果。 Tasi 等人(2002)[11]利用氣泡幫浦於微管道內產生氣泡，造成兩液體發生混合，並借其爆炸時使流體接觸面積產生波浪形狀，增加流體接觸面積，進而增加混合效果。在氣泡產生頻率為 200HZ 時，此波浪流之波長小於理論之擴散距離，因此會有較佳之混合效果。 Tang等人(2004)[12]由數值模擬與實驗方式,探討具有溝槽底部之Y 形微混合器（參考圖 1.9）其幾合尺寸與流體旋轉之關係，在數值模擬方面:兩進口流量相等，範圍從 5mm/s~50mm/s,流體擴散係數D(diffusion coefficient)為 1.26×10-5cm2/s，管道長度 10mm,寬度變化由 100~600mm 深度變化由 100~400mm，溝槽寬度與深度變化由 50~400mm；在實驗方面：管道之長寬深分別為 30mm,600 mm,200 mm三種，溝槽之寬、深為 200 mm、225 mm兩種，溝槽總數共 56 個。因底部溝槽結構可使得流體產生扭曲，增加流體接觸面積，增加混合效果。其結果發現:較淺的流道產生較強的流體旋轉，但是混合效率並不一定是隨著流道寬度增加而增加，因為擴散距離增加亦同時影響混合效率。另外，溝槽高寬比值愈大(但須小於 1.5) ，對混合效果較有幫住；太寬或太窄的溝 7.

(21) 槽結構,會降低混合效率. Park 等人(2004)[13]使用流體分離(breakup)方法在低雷諾數時加強產生流體旋轉的效能，藉由微流道的螺旋外型產生旋轉流體，共設計兩種混合器（參考圖 1.10），分別研究管道區段數目與分離(breakup)方法對混合效率的影響，從第一種混合器可以發現，在高雷諾數值 Re=50) 彎曲區域愈多，混合效果愈好，但是在低雷諾數之情況下並沒有相同的結果，因為流體的扭曲在低雷諾數(Re<1)時並不容易發生，另以第二種混合器，利用 breakup 方法，在低雷諾數時流體產生連續破裂碎塊，因此增加接觸面積，經實驗結果，在雷諾數分別為 1、10、50 時，通過 4mm 長之管道，可達到近 70%的混合效率。 Schwesinger 等人(1996)[14]設計雙層微結構（如圖1.11），利用將兩種不同液體相疊加的觀念，增加液體的接觸面積，提升擴散效能，他們的作法是將AB二種液體由左右兩側合成一股之後，再將之分成上下兩股；如此，上下兩股中各有一半A液及B液。被上下分開的兩股，再導引至左右兩側，然後再次合流之後，再將之分成上下兩股流。至此，左右兩股中各有A、B液交互相疊加。經多次的合分，A、B流體以層狀交互相疊而幫助A、B兩液間的混合。. Knight等人(1998)[15]的研究中，利用10µm寬的微噴嘴（如圖 1.12），將液體A 注入由其兩側垂直方向流入的B液中，由於受兩側Ｂ液擠壓的關係，而形成如光學聚焦的作用，A液的厚度約在2µm以下，而且被B液包夾，如此透過二液間的交互擴散作用使得兩液迅速混合。. Kang等人(2004)[16]指出在微流道混合數值分析中，以傳統解動量及質傳方程式時，因流體介面變形甚巨或管道長度過長，在流體介 8.

(22) 面間之計算時，因數值擴散效應(numerical diffusion)或網格不夠緊密，無法獲得較為準確的答案。並礙於電腦處理效能，乃藉由流場之週期性變化，只求解某週期間之流場並獲得濃度分佈。在此，其利用新的數值方法，也就是”colored particle tracking method” 簡稱CPTM，分析 Slanted groove micromixer(SGM)、staggered herringbone micromixer(SHM)及barrier embedded micromixer(BEM)等三種流場速度及濃度分佈情形。其分析步驟有三，第一先以Galerkin/least-aquares(GLS) 數值方法解動量方程式並獲得週期性流場，再以質點追蹤方法獲得最後週期之質點分佈，最後利用所謂Mixing entropy(S)，獲得截面上濃度分佈，結果發現，混合效率優劣為SHM>BEM>SGM。. Chang等人(2004)[17]以數值方法模擬電滲流於微管道流體流道情形。在具有電滲流之微管道中加入方形阻塊，當流體通過阻塊時，因受阻塊作用使其流道受限縮減，因此增加擴散效應，當阻塊增加，混合效率提升，但因過多的阻塊，勢必需增加流道長度才可得到理想的混合效果，利用電動力學的觀念，於阻塊上方加入異質性物質，使其電性與管道表面相異，流體在阻塊上方將產生回流，迫使流過阻塊上方區域更窄，扭曲也更嚴重，可以改善混合效率。. 混合效果評估部分微混合器的研究，除了混合器的設計外，如何針對混合效能進行評估也是相當重要的。目前評估的方法不外乎是利用顯微鏡與CCD 達成偵測影像的目的，並藉由染色程度或螢光強度的分析達到量化的結果。螢光強度分析 Kamholz等人(2001)[18]利用具有螢光性質的物質螢光黃（fluorescein） 9.

(23) 進入流體，當兩種流體混合時，螢光強度的改變可在CCD顯微鏡下觀察，並藉由螢光強度的大小，透過程式量化獲得混合成效。Mengeaud 等人(2002)[19]在Z型微流道中，於一端進口處注入10mM磷酸鹽 (phosphate)，另一端注入螢光黃（fluorescein），亦透過混合後之螢光強度判斷混合成效。色度分析其原理利用兩種流體以不同顏色染色，藉其混合時之顏色差異判斷混合效果。Miyake等(1993)[20]人觀察透明水溶液與加入黃色染料（螢光素鈉）之混合情形，其可透過顯微鏡觀察黃色染料之分佈狀況，以決定混合效率。. 酸鹼度分析 Liu等人(2000)[8]在兩分岔流道中，一端注入濃度0.31mol/L之酚酞與乙基酒精之混合液，另一端注入0.33mol/L之氫氧化納與乙基酒精之混合液，其PH值為13，混合過程藉由酚太酸鹼指示計顏色變化，觀察出混合成效，此指示劑可能反應出色度的改變或螢光強度的改變。. 1.7研究動機與方向雖然主動式(active)微混合器之混合效果優於被動式(passive)混合器，但其需另外提供額外能量與外力之形式，不易整合或製作在晶片上，產生之電場或熱能亦可能傷害分析的生物分子，因此，被動式混合器簡易之製作方式及原理，成為混合器設計之主流，但在低雷諾數之平滑管道流場內，根據文獻[21]得知，擴散時間、擴散距離及擴散係. d2 ，再由Peclet number之定義， P e = Vd D ，可以數之關係如下: t = 2D 推得混合長度( Lm ) Lm =Pe×d，當Peclet No.值愈大，所需之管道長度愈 10.

(24) 長，一般微流道中 Peclet No.範圍在102~106間，由此可知，在實驗室晶片之設計上是非常不切實際，因此，必須結由改變流道外型結構，如增加阻塊、蜿蜒管道設計、底部溝槽設計等，藉此增加接觸面積，加速混合效果。本篇採非結構性網格，利用計算流體力學方法，模擬在溝槽式結構管道，藉由幾何外型參數之改變，探討對混合效率之影響，並討論新型態之外型對混合效果之改善情形。. 11.

(25) 第二章混合器幾何外型介紹本次討論混合器之幾何外型共有四種，分別為：直型傾斜溝槽混合器（Striaght ridge micromixer，簡稱SRM）、直型雙傾斜溝槽混合器（ Double Striaght ridge micromixer，簡稱 DSRM）、人字形溝槽混合器（staggered herringbone mixer，簡稱SHM）、具傾斜阻塊人字形混合器（簡稱OBSHM），其中SRM與SHM兩種外型為Stroock等人(2002)[1] 提出，在此藉由改變其幾何外型，探討溝槽深度（Gd）、溝槽寬度（Gw）、溝槽傾斜角度（θ）對混合效率之影響，另發展出兩種新型態混合器，其一將SRM單溝槽混合器於上壁面增加溝槽結構形成雙傾斜溝槽形混合器（DSRM），另將SHM管道中加入傾斜阻塊，發展成為具傾斜阻塊人字形混合器（OBSRM），透過數值模擬，藉此瞭解混合器之混合效率。相關之幾何外型及參數如下所示。. 2.1 直型傾斜溝槽形混合器（Striaght ridge micromixer；簡稱SRM）. Z Y Lc. X. θ. Y X. SRM示意圖. 其中： Cw:管道寬度；Ch: 管道高度；Gd：溝槽深度；Gw：溝槽寬度；θ：溝 12.

(26) 槽傾斜角度； Gp:溝槽間距(2*Gw)；Li：入口到第一個週期性幾何外型之距離； Lc：週期性幾何外型長度；Lo：出口到第最後一週期性幾何外型之距離 L：管道全長；幾何外型 A B C D E F G. 管道高管道寬溝槽深溝槽寬溝槽傾溝槽間度(Ch) 度(Cw) 度(Gd) 度(Gw) 斜角(θ) 距(Gp) 1d 2d 0.42d 0.5d 45 1d 參考 (Gd=0.2d) 1d 2d 0.2d 0.5d 45 1d (Gd=0.6d) 1d 2d 0.6d 0.5d 45 1d (Gw=0.25d) 1d 2d 0.42d 0.25d 45 0.5d (Gw=1d) 1d 2d 0.42d 1d 45 1d (θ=60) 1d 2d 0.42d 0.5d 60 1d (θ=75) 1d 2d 0.42d 0.5d 75 1d 變動參數. 幾何進口段長單週期長出口段長週期數單週期溝槽 (P) 數外型度（Li) 度(Lc) 度（Lo) A 2d 10d 3d 12 10 B 2d 10d 3d 12 10 C 2d 10d 3d 12 10 D 2d 10d 3d 12 24 E 2d 10d 3d 12 5 F 2d 10d 3d 12 10 G 2d 10d 3d 12 10 流道總長（L）＝Li+（Lc×P）+Lo. 總溝槽流道總數(N) 長（L） 120 125d 120 125d 120 125d 240 125d 60 125d 120 125d 120 125d. 2.2直型雙傾斜溝槽混合器（Double Striaght ridge micromixer，簡稱 DSRM）於SRM混合器上壁面增加溝槽，溝槽位置與下溝槽交錯排列，溝槽幾何參數與下溝槽相同，如下圖所示. 13.

(27) Z Y. Lc. X. θ. Y X. DSRM示意圖其中：. Cw:管道寬度；Ch: 管道高度；Gd：溝槽深度；Gw：溝槽深度；θ、ψ：溝槽傾斜角度；Gp:溝槽間距(2*Gw)；Li：入口到第一個週期性幾何外型之距離；. Lc：週期性幾何外型長度；Lo：出口到第最後一週期性幾何外型之距離 L:管道長度. 管道高度(Ch) 1d. 管道寬度(Cw) 2d. 溝槽深度(Gd) 0.42d. 溝槽寬度(Gw) 0.5d. 14. 溝槽傾斜角(θ) 45. 溝槽間距(Gp) 1d.

(28) 進口段長單週期出口段長週期數 (P) 度（Li) 長度(Lc) 度（Lo) 2d 10d 3d 12. 單週期上單週期下總溝槽流道總數(N) 長（L）溝槽數溝槽數 10 9 228 125d. 2.3人字形混合器(Staggered herringbone mixer；簡稱SHM). Z Y. Lc. X. θ Gl φ. Y X. SHM示意圖其中： Cw:管道寬度；Ch: 管道高度；Gd：溝槽深度；Gw：溝槽深度；θ、φ：溝槽傾斜角度；Gp:溝槽間距(2*Gw)；Li：入口到第一個週期性幾何外型之距離； Lc：週期性幾何外型長度；Lo：出口到第最後一週期性幾何外型之距離 L：管道長度長溝槽短溝槽管道高管道寬溝槽深溝槽寬溝槽間對稱性傾斜角傾斜角度(Ch) 度(Cw) 度(Gd) 度(Gw) 距(Gp) 比(Gl/Gs) (θ) (φ) 1d 2d 0.42d 0.5d 45 -45 1d 2. 15.

(29) 進口段長單週期出口段長週期數 (P) 度（Li) 長度(Lc) 度（Lo) 1d 13.34d 2.28d 8. 單週期溝槽數 12. 總溝槽流道總數(N) 長（L） 96 110d. 2.4.具傾斜阻塊人字形混合器(staggered herringbone micromixer with oblique block；簡稱OBSHM). Z. Lc. Y X. θ. ψ. Y X. OBSHM示意圖其中： Cw：管道寬度；Ch：管道高度；Gd：溝槽深度；Gw：溝槽深度；θ、 ψ：溝槽（阻塊）傾斜角度；Gp：溝槽間距(2*Gw)；Li：入口到第一個週期性幾何外型之距離；Lc：週期性幾何外型長度；Lo：出口到第最後一週期性幾何外型之距離 L：管道長度；Bl：阻塊長度；Bw：阻塊寬度；Bh：阻塊高度；. 管道高度(Ch). 管道寬度(Cw). 溝槽深度(Gd). 1d. 2d. 0.42d. 溝槽寬長溝槽傾短溝槽傾溝槽間對稱性度(Gw) 斜角(θ) 斜角(φ) 距(Gp) 比(Gl/Gs) 0.5d. 45 16. -45. 1d. 2.

(30) 阻塊長度(Bl). 阻塊寬度(Bw). 阻塊高度(Bd). 0.5d. 0.67d. 1d. 阻塊傾阻塊傾斜角(θ) 斜角(φ) 45 -45. 進口段長單週期出口段長週期單週期總溝單週期阻塊流道總度（Li) 長度(Lc) 度（Lo) 數(P) 溝槽數槽數阻塊數總數長（L） 1d 13.34d 2.28d 8 12 96 11 88 110d. 17.

(31) 第三章. 數學模式. 3.1 基本假設(Basic Assumption) 本篇模擬之混合器具有週期性幾何外型之特性，為減少建立之網格數量過多造成電腦計算過於耗時的困擾，因此必須獲得單一週期之流場分佈及特性，藉此重複運用此單一週期速度場分佈，計算其下游週期之濃度分佈。我們將採兩階段進行，第一階段先計算具有進出口之管道（如圖 3.1 所示），獲得單一週期之流場分佈，第二階段則是週，利用第一階段得到之速度場計算濃度期性流場區域（如圖 3.2 所示）分佈。本研究流場的基本假設如下：. 1. 三維(3-Dimensional)、穩態(Steady)且不可壓縮(Incompressible)之流場。. 2. 工作流體(Working Fluids)為牛頓流體(Newtonian Fluids)。 3. 層流(Laminar Flow)、恆溫流場(Isothermal Flow)。 4. 不考慮物體力(Body Force)。 5. 混合過程不考慮化學反應產生。、密度（density）及黏滯力(Viscosity) 6. 擴散係數（Diffusion coefficient）不隨濃度變化而改變，為一定值。. 7. 壁面與流體流體之間無滑移現象存在。 8. 流場具有週期性速度場特性。 3.2 統御方程式(Governing Equation) 由以上的假設，我們可以推得所需的統御方程式如下：連續方程式：(Continuity Equation). div( ρV ) = 0. (3.1). 動量方程式 ( Navier-Stokes Equation ) 18.

(32) div ( ρVV ) = − gradP + div ( µgradV ). (3.2). 質傳方程式 (Mass Transfer Equation) div ( ρVc ) = div ( Dgradc ). (3.3). 計算質傳方程式時，因假設流體性質不隨濃度變化而改變，因此，不與連續方程式、動量方程式聯立求解。其中， V 為速度， ρ 為密度， P 為壓力， µ 為黏滯係數， D 為溶質在溶液中的擴散係數，C 為溶質的濃度（表示單位體積內所含之質量）。將上述之統御方程式無因次化，相關無因次化參數的定義如下：. p ∗ x∗ = x 、 y∗ = y 、 z∗ = z 、 u ∗ = u 、 v∗ = v 、 w∗ = w 、 p = d d d V V ρV 2 、 V Re = ρVd µ 、 c. ∗. = c / c0 、 P e = Vd D. 其中 d ：管道高度. V ：入口處的平均速度值。. c0 ：溶質的初始濃度經無因次化後的統御方程式如下：. ∂u∗ ∂v∗ ∂w∗ + + = 0 連續方程式： ∂x∗ ∂y∗ ∂z∗. (3.4). X 方向動量方程式：. u∗. ∗ ∂u∗ ∗ ∂u∗ ∂p∗ 1 ⎛ ∂2u∗ ∂2u∗ ∂2u∗ ⎞ ∗ ∂u v w + + = − + ⎜ + + ⎟ ∂x∗ ∂y∗ ∂z∗ ∂x∗ Re ⎝ ∂x∗2 ∂y∗2 ∂z∗2 ⎠. (3.5). Y 方向動量方程式： ∗ ∂v∗ ∗ ∂v∗ ∂p∗ 1 ⎛ ∂2v∗ ∂2v∗ ∂2v∗ ⎞ ∗ ∂v u ∗ + v ∗ + w ∗ = − ∗ + ⎜ ∗2 + ∗2 + ∗2 ⎟ ∂x ∂y ∂z ∂y Re ⎝ ∂x ∂y ∂z ⎠ ∗. Z 方向動量方程式：. 19. (3.6).

(33) ∗ ∂w∗ ∗ ∂w∗ ∂p∗ 1 ⎛ ∂2 w∗ ∂2 w∗ ∂2 w∗ ⎞ ∗ ∂w u +v + w ∗ = − ∗ + ⎜ ∗ 2 + ∗2 + ∗2 ⎟ ∂x∗ ∂y∗ ∂z ∂z Re ⎝ ∂x ∂y ∂z ⎠ ∗. (3.7). 質傳方程式： ∗ 1 ⎛ ∂2c∗ ∂2c∗ ∂2c∗ ⎞ ∂c∗ ∗ ∂c∗ ∗ ∂c u ∗ + v ∗ + w ∗ = ⎜ ∗2 + ∗2 + ∗ 2 ⎟ ∂x ∂y ∂z Pe ⎝ ∂x ∂y ∂z ⎠ ∗. (3.8). 3.3 邊界條件(Boundary Conditions) 3.3.1 第一階段區域邊界條件 1.速度邊界條件： (1)入口處為一均勻流場， u∗ = 1、 v∗ = 0 、 w∗ = 0. ∂u∗ ∂v∗ ∂w∗ (2)出口速度梯度為零， ∗ = ∗ = ∗ = 0 ∂x ∂x ∂x (3)壁面之邊界條件，採用無滑動(No Slip)的邊界條件，u∗ = v∗ = w∗ = 0 2.濃度邊界條件： (1)入口處半邊為 1，另一半邊為 0 Cw ≤ y ≤ Cw ) 、 c∗ 2. 濃度： c∗ = 1 (. = 0 (0 ≤ y <. Cw ) 2. (2)壁面為不可滲透表面，故壁面上之通量為零，也就是在壁面垂直 * 方向通量為零。 ∂c ∂n = 0 （ n 為與壁面垂直之方向）. (3) 出口濃度梯度為零 . ∂c* ∂x = 0 3.3.2 第二階段區域邊界條件以第一階段計算獲得單一週期速度場計算濃度分佈. 1.濃度邊界條件： * P * P −1 (P：週期) (1)入口: c in = c out. (2)壁面為不可滲透表面，故壁面上之通量為零，也就是在壁面垂直 * 方向通量為零。 ∂c ∂n = 0 （ n 為與壁面垂直之方向）. 20.

(34) (3) 出口濃度梯度為零 . ∂c* ∂x = 0 3.4 渦度計算: 為瞭解因底部溝槽結構管道其產生之渦度強弱與幾何外型之關係，因此計算在傾斜溝槽中心截面之渦度強度，其定義如下：. ⎛ ∂w ∂v ⎞ ⎛ ∂u ∂w ⎞ ⎛ ∂v ∂u ⎞ 渦度( Ω )= ∇×V = ⎜ − ⎟ i + ⎜ − ⎟ j + ⎜ − ⎟ k ⎝ ∂y ∂z ⎠ ⎝ ∂z ∂x ⎠ ⎝ ∂x ∂y ⎠ 平面法向量 n 沿 X 方向（ i ）之渦度強度如下：. Ωx = ⎡⎣(∇×V )in⎤⎦ii = (. ∂w ∂v − ) ∂y ∂z. 沿截面渦度總和 (Ωtol ) ：. Ωtol = ∑Ωxi × Ai. (3.9). 其中. Ω xi ：橫截面上第 i 個元素之渦度； Ai ：橫截面上第 i 個元素之面積截面渦度平均值 (Ω) ： Ω =. ∑Ω × A ∑A xi. i. (3.10). i. 3.5 混合指數的定義: 雖然從混合濃度之分佈狀態可大約判別混合效果之好壞，但仍需定義一明確之混合效率，將濃度分佈加以量化，以獲得明確指標，作為混合之參考依據。在此採用 Stroock(2002)[3]標準差之定義:. M=. ∑(C − C) ∑A i. 2. Ai. i. n. 其中. ∑C A C= ∑A i =1. i. i. i. 21. (3.11).

(35) M ：混合指數（mixing index）， M =0 表示完全混合； M =0.5 表示完全. 不混合 Ci ：截面上第 i 個元素之濃度 C ：截面上平均濃度. Ai ：截面上第 i 個元素之面積. 22.

(36) 第四章. 數值方法. 4.1 傳輸方程式離散化 (Discretization) 傳輸方程式數學模式如下：. div ( ρV Φ ) = div (Γgrad Φ ) + q (4.1) 等號左邊為對流項，右邊為擴散項及源項，利用高斯散度定理(Gauss. Divergence Theorem)與中點定理(Midpoint Rule)，將方程式離散化，便於轉化成數值方式計算。. 4.1.1 對流項 (Convection Term)離散化. ∫∫∫ (∇ ⋅ ρVΦ )dV = ∫∫ ( ρVΦ ) ⋅ dS ⇒ ∑ (ρVΦ ) ∆v. f. S. f. ⋅ S f = ∑ F fc. (4.2). f. C 而 Ff ≡ mf Φ f. (4.3). 其中. F c ：面上的對流通量. (Convection Flux)。. m ：質量流率(mass flow rate)。 S f ：面上之法向量。（下標 f ：控制體積之任一面上之中點）. 在計算對流項的過程，為求得準確度高且易收斂的解，採用中央差分和上風差分混合方式，表示方式如下： F C = F UD + γ (F CD − F UD ). (4.4). F CD ：中央差分法所得之對流通量。. F UD ：上風差分法所得之對流通量。 γ：範圍在 0<γ<1 ， (γ 值接近於 1 值時，可以確保高階準確) 在疊代過程中，為確保係數矩陣能具有對角優勢(Diagonally 23.

(37) Dominant)，故採 Deffered Correction 的方法[22]，即將上式中以上風法求得之第一項置於係數矩陣，另將中央差分法與上風差分法之差所得的第二項置於源項。. 4.1.2 擴散項 (Diffusion Term)離散化. ∫∫∫ ∇(Γ∇Φ)dV = ∫∫ Γ∇Φ ⋅ dS ⇒ ∑ ( Γ∇Φ ) ∆V. f. f. S. ⋅ S f = ∑ F fd. (4.5). f. 而 F fd ≡ ( Γ∇Φ ) f ⋅ S f. (4.6). 其中， F d ：面上的擴散通量 (Diffusion Flux)。. Γ ：係數參數，用於動量及質傳方程式有下列兩種情況：動量方程式→ Γ ＝µ（ µ ：黏滯係數）質傳方程式→ Γ ＝D（D：擴散係數）在求解擴散項時，另外定義向量 d 為沿主格點 P 至相鄰格點 C 方向的向量，因此 S f 可表示如下： S. f. = d +( S f - d ). (4.7). 將其代入(4.6)式後， F fD = Γ f ∇Φ f ⋅ S f = Γ f ∇Φ f ⋅ d + Γ f ∇Φ f ⋅ ( S f − d ). (4.8). 因為 d 的大小決定數值計算時之擴散量，影響計算時的穩定性，. 因此利用 over-relaxed approach[23]方法處理 d 值，方法如下： d ≡. S. f. ed ⋅ es. ed =. S. δ. 2. δ （ e s：沿 S f 之單位向量）. f. ⋅S. (4.9). f. 將(4.9)代入(4.8) FfD =. Af S f. δr ⋅ S f. 2. (Φ C − Φ P ) + Af ∇Φ f ⋅ ( S f − d ). (4.10). (4.10)式採 over-relaxed 方法獲得之第一項因具有較大的擴散係數，將其置於係數矩陣，獲得之第二項置於源項。 24.

(38) 4.1.3 源項(Source Term)離散化壓力梯度同樣採高斯散度定理及中點定理方式計算，因此 ∇P =. 1 ∆V. ∫∫∫ ∇ PdV. =. ∆V. 1 ∆V. ∫∫ Pd S. ⇒. s. 1 ∆V. ∑. Pf S. (4.11). f. f. （ ∆V ：單一計算網格的體積）. e i 方向之分量表示為: ∂P 1 = ∇ P ⋅ ei = ∂xi ∆V. ∑ P (S f. f. ⋅ ei ). (4.12). 4.1.4 邊界壓力之計算因計算壓力梯度項時，需利用邊界上之壓力值，推導邊界壓力，方式如下：. Pb − P p = ∇ P ⋅ δ. (4.13). b：壁面上的中點； δ ：P 到 b 之距離向量下列方式近似壓力梯度： ∇P =. 1 ∆V. ∑ f. Pf S. f. =. 1 ∆V. ⎛ ⎜ Pb S b + ⎜ ⎝. ∑. f ≠b. Pf S. f. ⎞ ⎟ ⎟ ⎠. (4.14). （下標 f ：除 b 以外之其餘各面）因此推導出邊界上之壓力為： ⎛ 1 ⎜ Pp + ∑ Pf S f ⋅δ ⎜ ∆V f ≠b ⎝ Pb = 1 ⎛ ⎞ Sb ⋅δ ⎟ ⎜1 − ∆V ⎝ ⎠. ⎞ ⎟ ⎟ ⎠. (4.15). 25.

(39) 4.2 解線性代數方程式將對流項、擴散項及源項合併後，以下列方法得到傳輸方程式之代數方程式： ApΦ. p. ∑. =. AC Φ C + Q. (4.16). C. 其中 A p =. ∑A. C. C. AC =. Γf Sf. 2. δ ⋅Sf. + max(−mf ,0). (4.17). 而源項 Q 表示如下：. {. }. Q = ∑ −γ ⎡⎣mf (wf Φp + (1− wf )ΦC ) −⟨max(mf ,0)ΦP + min(mf ,0)ΦC ⟩⎤⎦ c. (. (4.18). ). +∑ Γf ∇Φf Sf − d + qp∆V f. 下標 C：指控制體積周圍網格之中心格點。下標 f ：指控制體積任一面上之中點。另考量(4.16)式疊代的穩定性，加入一所謂鬆弛因子 β. (under-relaxation factor，0< β <1 )[22]，故修正後之動量代數方程式如下所示： AP. β. Φ P (κ +1) = ∑ AC Φ C ( κ +1) + QΦ + C. 1− β. β. AP Φ P (κ ). (4.19). 其中上標 k+1 代表新值，k 則代表前次疊代之值。. 4.3 壓力速度耦合關係式根據 Patankar 之 SIMPLE 法則[24]，從動量方程式計算獲得之速度、壓力修正關係式，將其代入連續方程式後，即可獲得壓力修正方程式，然後用其計算出之壓力修正值，重新修正速度及質量流率，其方式如下：. 4.3.1 計算面上的質量流率 26.

(40) 在滿足連續方程式之條件下，需先求得面上之質量流率，而要獲得面上之流量必須獲得面上之速度，因此先要瞭解速度與壓力之關係式為何，再利用內差法差分計算面上的速度與壓力關係式，最後求得面上質量流率。在此以類似 Rhie & Chow[25]之線性內插法來計算面上速度。方式如下：將(4.16)式中之壓力項自源項提出，可得到主格點之速度與壓力關係式： V. ⎞ ⎟ ⎟ ⎠. ⎛ ∆V = H − ⎜ ⎜ A p ⎝. p. ∑. 其中 H =. ∇ Pp. (4.20). p. A C V C +Q. '. Q'：不含壓力項之源項. C. Ap. (4.21). 以此類推，格子面上之速度可表示為： V. f. = H. 其中. H. ⎛∆ V − ⎜⎜⎜ ⎜⎝ A p. f. = V. f. f. ⎞⎟ ⎟⎟ ∇ P f ⎠⎟ f. ⎛∆ V + ⎜⎜⎜ ⎜⎝ A p. (4.22). ⎞⎟ ⎟⎟ ∇ P ⎠⎟ f. (4.23). f. 將(3.23)代入(3.22)式得: ⎛∆V V f = (V f + ⎜⎜⎜ ⎜⎝ A p. ⎞⎟ ⎛∆V ⎜ ⎟⎟⎟ ∇ P f ) − ⎜⎜ ⎜ ⎠ ⎝ Ap. ⎟⎞ ⎟⎟⎟ ∇ P f ⎠f. (4.24). 上標”—” 表示由控制體積之中心格點 P 及相鄰共面 f 之 C 格點內插而得，其表示如下： V f = w p V C + (1 − w p )V p. (4.25). ∇P f = w p ∇PC + (1 − w p )∇Pp. (4.26). ⎛ ∆V ⎜ ⎜ A p ⎝ ⎛ ∆V ⎜⎜ ⎝ Ap. ⎞ ⎟ ⎟ ⎠. 可由下式方法近似： f. ⎞ 1 ⎡⎛ ∆V ⎟⎟ = ⎢ ⎜⎜ 2 ⎢⎝ Ap ⎠f ⎣. ⎞ ⎛ ∆V ⎟⎟ + ⎜⎜ ⎠C ⎝ Ap. ⎞ ⎤ ⎟⎟ ⎥ ⎠ P ⎥⎦ 27.

(41) (4.27) 則面上的質量流率可表示為： m. f. = ρV f ⋅ S. f. = ρ fV f ⋅ S − ρ. f. ⎛ ∆V ⎜⎜ ⎝ Ap. ⎞ ⎟⎟ ⎠f. (∇ P. f. ⋅S. f. − ∇P. f. ⋅S. f. ). (4.28). 4.3.2 壓力修正方程式依據SIMPLE法則，求解動量線性代數方程式，以前一次疊代求解之格點中心 P 點速度場V. ∗. 與壓力P*仍未滿足連續方程式，因此需再次. 修正，因此先利用動量方程式所得之速度及壓力關係式，再代入連續方程式，可獲得壓力修正方程式，藉此獲得壓力修正量p’，並加以修正速度與質量流率，獲得修正後之速度及壓力分別為V ⎛ ∆V V p′ = − ⎜ ⎜ A p ⎝. **. 、 P **. ⎞ ⎟ ∇ P p′ ⎟ ⎠. (4.29). ** * 其中格點中心之壓力修正量 Pp′ = Pp − Pp. (4.30). 同理，網格面上速度之速度修正式： V f′ = V. ** f. − V. * f. ⎛ ∆V = −⎜ ⎜ A p ⎝. ⎞ ⎟ ∇ P f′ ⎟ ⎠f. (4.31). 由上式可繼續推得質量流率： m *f* = m *f + ρ f ν 'f ⋅ S f = m *f − ρ f ( ⎛ ∆V = m − ρf⎜ ⎜ A ⎝ p * f. ∆V ) f ∇ Pf' ⋅ S f Ap. ⎞ ⎛ ⎞ ⎟ ∇ Pf' ⋅ d − ρ f ⎜ ∆ V ⎟ ∇ Pf' ( S f − d ) ⎜ A ⎟ ⎟ ⎝ p ⎠f ⎠f. (4.32). 2. ⎛ ∆V ⎞ ⎛ ∆V ⎞ S f ⎟ (PC' − Pp' ) − ρ f ⎜ ∇P′ ⎟ ⋅ (S f − d ) = m*f − ρ f ⎜ ⎜A ⎟ ⎜ A ⎟ δ ⋅S f ⎝ p ⎠f ⎝ p ⎠f. (4.33). 令修正後之速度場滿足連續方程式：. ∑. m. ** f. = 0. (4.34). f. 將(4.33)代入(4.34)則可得到壓力修正方程式： 28.

(42) A p P p′ =. ∑. AC PC′ + S p 1 + S p 2. (4.35). C. 其中： AC. ⎛ ∆V = ρf⎜ ⎜ A ⎝ p. AP =. ∑. Sf ⎞ ⎟ ⎟ δ ⋅S ⎠f. 2. f. AC. c. S. p1. = −. ∑. m. * f. (4.36). f. ⎞ ⎛ ∆V S p2 = ∑ ρ f ⎜ ∇P ′ ⎟ ⋅ ( S f − d ) ⎟ ⎜ A f ⎠f ⎝ p. (4.37). 4.3.3 解壓力修正方程式求解壓力修正方程式時，因 S P 2 項中包含 P ′ ，可以採用下列兩步驟連續修正 successive correction[22]近似︰步驟一：只考慮含 S. A p P p′ ( 1 ) =. p1. 部分，計算獲得第一次壓力修正量 P ′ ( 1 ) ，表示方式如下：. ∑. A C PC′ ( 1 ) + S. (4.38). p1. C. 步驟二：再利用求得之壓力修正量 P ′. (1 ). 計算 S. p2. 部分，並再次求解壓力修正. 方程式，獲得第二次壓力修正量 P ′ ( 2 ) A p P p′ ( 2 ) =. ∑. A C PC′ ( 2 ) + S p( 12). (4.39). C. ⎛ ∆V. ⎞ ∇P ′ (1) ⎟ ⋅ ( S f − d ) ⎟ ⎝ Ap ⎠f. (1) 其中 S p 2 = ∑ ρ f ⎜⎜ f. (4.40). 可重複疊代上述第二步驟以獲得更高階之修正量，本文使用連續修正兩次，待求解壓力修正方程式後，即可得到修正後之壓力。 29.

(43) Pp ** = Pp * + Pp '(1) + Pp ′(2) ∗∗. 修正後的新速度 V 如下： ⎛ ∆V V f∗∗ = V f∗ − ⎜⎜ ⎝ AP. ⎞ ⎟⎟ ∇Pf′ ⎠f. (4.41). 由(4.33)式，經兩步驟修正獲得新的面上質量流率： ∗∗ ∗ 步驟一：m f = m f. 2 ⎛ S V ∆ ⎜ −ρf⎜ ⎜ AP δ r ⋅ S ⎝. ⎛ ∆V ∗∗ ∗ 步驟二： m f = m f − ρ f ⎜⎜ ⎝ AP. ⎞ ⎟ ⎟⎟ ( PC′ − PP′ ) ⎠f. ⎞ ⎟⎟ ∇ Pf′ ⋅ S f − d ⎠f. (. (4.42). ). (4.43). 4.4 解題步驟 0. 1.先給定初始速度、壓力值V 0 、 P 。 2.再解動量方程式獲得預測速度場V * 。 3.計算質量流率 mi * ,再由壓力修正方程式推得 P ′ 。 i **. 4.修正速度、壓力及質量流率，得到V ** 、 P ** 、 m 。 5.將修正後的新值作為初始值，重覆步驟 2~4，直至獲得收斂解。 6.將獲得收斂解之V. **. 代入質傳方程式，經多次疊代獲得濃度 C 之收斂. 解。進行求解代數方程式之係數矩陣，一般結構性網格具有數個對角線的. banded matrix，可以利用 line solver 及 Stone’s SIP 方法求解，但使用非結構性網格時，其並不具有 banded matrix 之特性，雖可以使用 point. slover 求解 unbanded matrix，但其收斂速度太慢，在此使用 preconditioned conjugate gradient method;而動量方程式不具對稱性質，採用 BICG(Biconjugate Gradient)法;壓力修正方程式其係數矩陣具對稱性質，故採用 ICCG(Incomplete Cholesky Conjugate)法，另外 30.

(44) 對質傳方程式則是採用 GS(Gauss-Seidal)法。. 31.

(45) 五、結果與討論 5.1 網格測試與實驗值比較網格密度影響模擬結果之正確性，網格太少，模擬結果將與實驗誤差很大，網格過密，將耗費太多之電腦運算時間，因此挑選適當網格數目非常重要。本文與 Stroock(2002)[3]等人實驗結果比較，分別測試三種網格密度，比較直型傾斜溝槽混合器（Striaght ridge micromixer）與人字形混合器(Staggered herringbone mixer)在雷諾數(Re)~0.01(網格測試定義之雷諾數根據實驗定義，與本文定義有所不同)，裴立數. (Pe)~2000 時，網格疏密對實驗值的影響，在數值計算中，採中央差分法（γ=1）求解速度及壓力場，濃度場則是結合中央差分及上風法之混合法（γ=0.8）求解，速度、壓力及濃度使用之鬆弛因子分別為 0.8、0.7、. 0.8、0.8 及 0.9。 5.1.1 直型傾斜溝槽形混合器(SRM)網格測試分別測試網格主流道截面格點數為 24x4 , 35x10 及 48x17 三種，如下表所示：. 主流道橫截溝槽橫截週期性幾何總網格面網格數面網格數外型網格數數 24x8 24x4 19200 測試網格1 27648 35x10 35x5 70000 測試網格2 108000 48x17 48x7 147600 208800 測試網格3. 備考 ∆x:0.125d∆y：0.08d;∆z：0.08d ∆x:0.06d;∆y：0.06d; ∆z：0.07d ∆x:0.06d;∆y：0.04d;∆z：0.04d. 測試結果如圖 5.1 所示，在橫向網格數為 48x17 時，模擬值與實驗值較為接近，而網格愈粗，差距則愈大。. 5.1.2 人字形混合器(SHM)網格測試主流道橫溝槽橫截面週期性幾何總網格數截面網格網格數外型網格數 24x8 24x4 24576 測試網格1 28800 36x5 92160 測試網格2 36x10 108000 48x7 201984 237312 測試網格3 48x17 32. 備考 ∆x:0.125d∆y：0.08d;∆z：0.08d ∆x:0.06d;∆y：0.06d; ∆z：0.07d ∆x:0.06d;∆y：0.04d;∆z：0.04d.

(46) 測試結果如圖 5.2 所示，與傾斜溝槽形混合器網格測試結果相同，在橫向網格數為 48x17 時，模擬值與實驗值較為接近。由上述結果，後續模擬網格，將採用橫向截面為 48x17 的網格密度，也就是各方向之單位 Cell 長度為：∆x：0.06d；∆y：0.04d；∆z：0.04d，而雷諾數（Reynolds Number）與裴立數(Peclet Number)與實驗值相同，分別以 Re=0.01；Pe=2000 進行數值計算。. 5.2 SRM 混合器混合效能分析與探討 5.2.1 計算網格共建立兩個計算網格，第一個計算網格（如圖 5.3）主要求得單一週期內之速度場分佈，第二個計算網格（如圖 5.4）主要目的在計算下游各週期之濃度場分佈情形，格點總數因幾何外型不同而有些差異，詳細數據如下表所示. 幾何外型 A(Ref.) B C D E F G. 週期性網格格點數 216000 198000 234000 216000 216000 216000 216000. 總格點數 306000 288000 324000 294300 306000 306000 306000. 5.2.2 流場特性. 為了解速度場分佈情形對混合效率之影響，因此觀察第一週期內橫向(Y-Z)截面之流場分佈狀態。選取之橫向截面方向為由流體出口處往入口處觀察所得之截面。分別截取 13 個橫向截面（A~M）位置（如，便於觀察流場對濃度分佈之影響。由橫向截面速度、流線及圖 5.5） 33.

(47) 濃度分佈圖（如圖 5.6）明顯發現，一逆時針方向旋轉之流場形成，其原因乃由底部傾斜溝槽結構造成，傾斜溝槽使流體產生沿橫向(Y 軸) 方向之速度分量，在管道下方溝槽區域流體沿溝槽向右流動，撞擊壁面後上翻進入管道後向左向下流動，在逆時針旋轉之流場帶動下，原本兩流體介面將慢慢向左偏移（上壁面比靠近溝槽區快），而管道右下角區域則是緩緩往上壁面移動，使二流體介面拉伸變形，截面 A~D 剛好是一個週期，到截面 E 時又回復到與截面 A 之流場型態。由此可之，底部溝槽結構引起之逆時針旋轉之二次流場，可使流體介面扭曲，增加流體接觸面積，為增加混合效果之主因之一。. 5.2.3 SRM 混合器幾何外型對混合效率之影響為了瞭解混合器幾何外型對混合效率之影響，因此變化其溝槽深度、溝槽寬度及溝槽傾斜角度，其餘幾何外型參數不變情況下觀察其對混合之影響，由圖 5.7~5.9 模擬結果發現，三組不同溝槽深度（0.2d、. 0.42d及 0.6d）中，溝槽深度Gd=0.6d時混合效果最好，溝槽深度Gd=0.2d 時混合最差；不同溝槽寬度（0.25d、0.5d及 1d）中，溝槽寬度Gw=1d 時混合效果最好，溝槽寬度Gd=0.25d時混合最差；當溝槽傾斜角度為. 45°、60°及 75°三種情況時，溝槽傾斜角度θ=45°時混合最好，θ=75°最差。後續將從質傳方程式計算之濃度分佈、質點追蹤方法獲取各截面分佈狀態、及渦度強度三方面探討不同幾何參數影響混合效率之主因。. 1.溝槽深度(Gd)對混合效果之影響（1）質傳方程式濃度計算結果分析：觀察各週期位置（如圖 5.10 所示）橫向截面濃度分佈情形，由圖. 5.14~5.16 可發現，兩不同流體介面將隨著下游週期位置增加而扭曲，其中又以溝槽深度(Gd)為 0.6d時所造成之介面拉伸最長，0.25d最短，當到達最後一週期位置，流體旋轉之角位移亦以 0.6d最多，約兩圈 34.

(48) (720°)，0.42d其次（約 630°），而 0.2d最差（約 360°），在最後一週期濃度分佈，可清楚分辨出在Gd=0.6d時，混合效果最佳，而Gd=0.2d最差。（2）質點追蹤方法於入口右側初始位置放置 800 個質點，可以獲得不同質點所行走之路徑軌跡，在穩態流場中也就是所謂的流線，藉此可以獲得不同截面各質點位置分佈情形。此種方法可以避免因進行質傳方程式計算造成之數值擴散，透過質點分佈情形可判斷混合之好壞。由圖 5.17~圖 5.19，可以清楚發現，因為管道底部傾斜溝槽結構，使得靠近底部的質點沿橫向方向拉伸，逆時針方向慢慢旋轉向渦度中心移動，藉由介面增加，可以縮短擴散所需之時間，我們以質點旋轉一圈為參考，溝槽深度Gd=0.2d（如圖 5.17）時需在第 12 個週期才可完成，. Gd=0.42d（如圖 5.18）時約在第 10 週期完成，而Gd=0.6d（如圖 5.19）時，只需到第 6 週期即可形成，明顯看出溝槽愈深，混合效果較佳。（3）渦度強度分析為瞭解溝槽深度與渦度大小之關係，依據(3.9)式(3.10)式計算A’與. B’傾斜中心截面位置（圖 5.20）沿X方向渦度大小，在此之前，先觀察其傾斜截面中心點A、B位置（圖 5.21）之橫向速度（V）分佈，從圖. 5.22 中心點A可以看出，隨著高度增加，在溝槽最高位置附近有最大的，之後慢慢速度遞減並轉為Y-方向，並隨高度增橫向速度（往Y+方向），但其最大速度遠小於Y+方向，整個流場為加有大速度值（往Y-方向）逆時針旋轉型態。從最大橫向速度比較，溝槽深度最淺者（Gd=0.2d），但在Y-方向的橫最小，而 Gd=0.42d及Gd=0.6d最大（兩者幾乎接近），較深的溝槽（Gd=0.6d）向速度以Gd=0.6d最大。在中心點B（不含溝槽）有有較大的橫向速度（沿Y-方向）。由速度及渦度場分佈（圖 5.23~圖. 5.28），我們可以發現在A’截面上之速度場皆呈現逆時針旋轉，至於渦度分佈，靠近壁面及溝槽位置之渦度皆為負值，而在中間區域則為正 35.

(49) 值，計算結果顯示（參考表 5-1），渦度大小以Gd=0.6d最強，Gd=0.42d 次之，而Gd=0.2d最小，另在B’截面上，流場則由右下角流至左下角，渦度在靠近上、左、右壁面為負值，其餘為正，渦度計算結果（參考，亦是以Gd=0.6d最強，而Gd=0.2d最小，因此可以表 5-1(a)~ 表 5-1(b)）得知，溝槽愈深，可以獲得較強之渦度，增加混合效果。. 2.溝槽寬度(Gw)對混合效果之影響（1）質傳方程式濃度計算結果分析：同樣觀察各週期位置（如圖 5.29~圖 5.30 所示）橫向截面濃度分佈情形，並與參考之幾何外型比較，整個介面扭曲的情形以溝槽寬度最寬（Gw=1d）時最為激烈，寬度愈狹窄（Gw=0.25d）效果較差，從流體旋轉角度看來，在最後一週期，最寬的溝槽（Gw=1d）時，流體旋轉將近 720°以上，而最窄的溝槽（Gw=0.25d）流體旋轉約 360°，明顯不利於兩流體介面間持續的加速變化。（2）質點追蹤方法從質點分佈現象圖 5.31、圖 5.32 可以看出，當溝槽寬度較窄時，質點形成之條紋寬帶分佈仍以片狀區域為主，但隨溝槽（Gw=0.25d），流體因旋轉使得由塊狀區慢慢形成帶狀連續區域，寬度增加（Gw=1d）擴散距離縮短，利於兩流體間分子擴散進行。（3）渦度強度分析由圖 5.33，流體的中心速度隨溝槽愈寬而增加，因此得知，窄型溝槽不利於流體在橫向方面速度之發展；至於A’與B’截面之流場及渦度場（圖 5.34~圖 5.37）分佈情形大致與前面相同，在最寬的溝槽中有最大的渦度值及面積。另從渦度強度計算所獲得之結果（如表 5-1(c) 及表 5-2(d)），當溝槽寬度Gw=1d時，渦度強度遠大於其他兩者，而 Gw=0.25d，獲得之渦度強度最小。. 3.溝槽傾斜角度（θ）對混合效果之影響（1）質傳方程式濃度計算結果分析： 36.