流體在平行板間受導電梯度與水平電場作用之動力穩定特性分析

國立臺灣大學工學院應用力學研究所 碩士論文

Graduate Institute of Applied Mechanics College of Engineering

National Taiwan University Master Thesis

流體在平行板間受導電梯度與水平電場 作用之動力穩定特性分析

Instability of electrokinetic flow in a thin fluid layer between two parallel plates with an electrical

conductivity gradients

蔡英廷 CAI,YING-TING

指導教授：陳發林 博士 Advisor : Falin Chen, Ph.D.

中華民國 98 年 6 月

June, 2009

致 謝

首先感謝的是指導老師陳發林教授，老師在校務百忙中亦能抽空指 導學生論文，並時常關心學生的論文進度與研究情緒。是否有其他外務 耽誤到論文進度或有其他生活上的困難，常讓學生對於老師貼心的問候 銘感五內。大同大學的張敏興教授，亦在每週關心學生論文進度，並給 予學生鼓勵，激勵學生對研究的熱情，使學生進度掌握得宜，口試順利 結束，且亦建議論文修改事項，使學生更明瞭論文缺漏與不足處，甚感 念老師用心處。亦非常感謝羅安成博士，有學長的幫助，使學生的理論 推導與程式的撰寫十分順利。遇到不懂的問題，學長總是殷殷期盼學弟 能自我學習，幫助學生能得到更多學習的機會，使學生更能深入了解已 身所不足處。且學長亦關心學生須注意論文撰寫與引用須小心，以免有 抄襲的問題發生。在在都非常感謝學長的教誨，非常謝謝。

博一的張森圳學長，更是良師益友，藉由跟學長的討論中，釐清學 習的漏洞，與未竟周全處，非常謝謝學長的抽空指導。而上屆的高鳴壕、

郭至益、任博川學長，陪伴學生苦悶的碩一生活，並給予心靈的支助，

時時幫學生打氣，非常謝謝。碩二的陳嘉倫、林佳慧、施經緯同學，感 謝你們，陪我度過這兩年的碩班生活，增添平日生活的趣味。碩一的陳 冠宏、陳立言、蔡之緯學弟，感謝你們豐富了我碩二的生活，使我每天 能以愉快的心情從事研究。實驗室的王姐、小藍、呂大哥，謝謝你們的 鼓勵，使實驗室的氣氛很融洽，謝謝你們。大學國樂社的朋友，志瑋、

開平、立中、東憲、賀婷、孟瑩、康康、球球、珮如、信如、冠伶、祥 和，謝謝你們的陪伴與邀約，讓我的碩士生活不再無聊，謝謝。還有千 緻也謝謝你對我的體貼與體諒，一直陪在我身邊，幫我加油。最後，謝 謝蔡明輝先生與陳美玉女士，感謝您讓我有個安心的學習後盾，謝謝你 們。

I

目 錄

目錄... I 圖目錄... IV 表目錄... VII 中文摘要... VIII 英文摘要... X 符號說明... XII

第一章 緒論... 1

1.1 研究背景... 1

1.2 文獻回顧... 2

1.3 研究動機與目的…... 4

1.4 研究方法…... 5

第二章 理論模式... 7

2.1 模型建立與假設... 7

2.2 統御方程式... 8

2.3 邊界條件... 10

2.4 基態解... 12

2.5 系統無因次化... 14

II

2.6 線性微擾化……... 19

2.7 正規模態展開... 20

2.8 忽略 EDL 效應... 24

第三章 數值分析方法... 28

3.1 謝比雪夫配置方法... 28

3.2 應用謝比雪夫配置方法... 30

3.3 QZ 演算法... 31

3.4 問題與討論... 31

第四章 計算結果... 33

4.1 概論... 33

4.2 忽略邊界 EDL 效應... 33

4.2.1 改變導電梯度的影響... 33

4.2.2 改變電施密特數的影響... 35

4.2.3 縱向模態... 36

4.3 考慮邊界 EDL 效應... 38

4.3.1 改變導電梯度的影響... 38

4.3.2 改變電施密特數的影響... 39

4.3.3 改變n的影響... 40

4.3.4 改變Q的影響... 40

III

4.3.5 改變

ψ

4.3.6 縱向模態... 42

4.4 綜合討論... 43

第五章 結論與展望... 64

參考文獻... 67

IV

圖 目 錄

圖 2.1 模型示意圖... 7

圖 2.2 EDL 示意圖... 11

圖 2.3 忽略 EDL 效應模型示意圖... 24

圖 4.1 橫縱向模態示意圖(a)平面；(b)立體。... 47

圖 4.2 當忽略 EDL，固定Sc_e =500的情形下，觀察H=1,3,10, 100 其(a)中性曲線；(b)波數與虛部成長率的關係。... 48

圖 4.3 當忽略 EDL，固定Sc_e =500的情形下，在不同導電梯度 的臨界流場形態(a)H=1,(k ,_c Ra )=(2.5,7474.3)；(b)H=10, _ec (k ,_c Ra )=(2.1,335.9)；(c)H=100, (_ec k ,_c Ra )=(2.3,199.9)。 49 _ec 圖 4.4 當忽略 EDL，討論其Sc_e ≥100 情況下，觀察 H=1 其(a) 中性曲線；(b)波數與虛部成長率的關係。... 50

圖 4.5 當忽略 EDL，討論其Sc_e ≥100 情況下，觀察 H=10 其(a) 中性曲線；(b)波數與虛部成長率的關係。... 51

圖 4.6 當忽略 EDL， Sc_e =500、k =2.1的情形下，觀察(a)H= 10,50,100 時，angle 與Ra 的中性曲線；(b)H=10, _e Ra_e = 10 時，angle 與7 s 的關係。... 52 _r 圖 4.7 當考慮 EDL，固定Sc_e =500、Q= 、1 n= −0.333、ψ =1 的情形下，觀察H=1,10,20,50,100 其(a)中性曲線；(b)波 數與虛部成長率的關係。... 53 圖 4.8 考慮 EDL 效應，固定固定Sc_e =500、Q = 、1 n= −0.333

、ψ =1的情形下，不同導電度的臨界流場形態。(a)H=1, (k ,_c Ra )=(2.5,5252.5)；(b)H=10, (_ec k ,_c Ra )=(2.1,339.9)； _ec

V

(c) H=100, (k ,_c Ra )=(2.3,203.1)。... 54 _ec 圖 4.9 當考慮 EDL，固定 H=10、Q= 、1 n= −0.333、ψ =1的 情形下，觀察Sc_e ≥100時，其(a)中性曲線變化；(b)波數 與虛部成長率的關係。... 55 圖 4.10 當考慮 EDL，固定 H=10、Q= 、1 n= −0.333、ψ =1的

情形下， 觀察n= −0.25, 0.333, 0.5, 1− − − 的(a)中性曲線變 化；(b)波數與虛部成長率的關係。…... 56 圖 4.11 當考慮 EDL，固定 H=10、n= −0.333、Sc_e =500、ψ =1 的情形下，觀察Q=1,10,100,10³時，其(a)中性曲線變化

；(b) 波數與虛部成長率的關係。... 57 圖 4.12 則塔電位效應圖。當固定 H=10、n= −0.333、Sc_e =500

、Q =1 的情形下，觀察ψ =0.1,1,10 的(a)中性曲線；(b)波 數與虛部成長率變化。... 58 圖 4.13 則塔電位效應圖。當固定 H=10、n= −0.333、Sc_e =500

、Q= 的情形下，觀察1 ψ =10,100,10 ,10^{3 4}與忽略EDL

效應的(a)中性曲線；(b)波數與虛部成長率變化比較。.. 59 圖 4.14 則塔電位影響圖。固定 H=10、Sc_e =500、n= −0.333、 Q= 的情形下，考慮1 ψ =0.1,1,10的臨界流場形態。(a) ψ =0.1, (k ,_c Ra )=(2.2,599.8)；(b)_ec ψ =1, (k ,_c Ra )=(2.1, _ec 339.9)；(c)ψ =10, (k ,_c Ra )=(2.1,336.1)；... 60 _ec 圖 4.15 當考慮 EDL 效應，Sc_e =500、k =2.1、n= −0.333、

Q= 、1 ψ =1的情形下，觀察(a) H=10,50,100 時，angle 與Ra 的中性曲線；(b)H=10,_e Ra_e =10⁷時，angle 與s 的 _r 關係圖。... 61 圖 4.16 固定Sc_e =500, 1Q= , 0.333n= − 的情形下，觀察ψ =0.1,

VI

1,10及忽略EDL 效應時，其(a)H 與Ra 的變化；(b)H _ec

與k 的變化。... 62 _c 圖 4.17 固定Sc_e =500, 1Q= , 0.333n= − 的情形下，觀察ψ =0.1, 1,10及忽略EDL 效應時，其(a)H 與Ra 的變化；(b)H 與 _ec

k 的變化。... 63 c

VII

表 目 錄

表 2.1 特徵尺度表... 15 表 4.1 物理系數表... 46

VIII

流體在平行板間受導電梯度與水平電場 作用之動力穩定特性分析

研究生 : 蔡英廷 指導教授 : 陳發林

國立台灣大學應用力學研究所

摘要

在微流體裝置(microfluidic device)中，以電力取代壓力來驅動流 體是最有效的方式。電流體力學(electrohydrodynamics, EHD)整合了電 學與流體力學，產生了相當豐富有趣的新現象。特別是當流體內部具 有導電度梯度(conductivity gradient)存在時，通電的結果導致自由電 荷的堆積，進而引發不穩定對流現象的發生。這種不穩定現象最常被 應用至微流體的混合器上，因此是目前微機電技術中相當重要的研究 課題之一。本論文主要是在利用線性穩定性分析方法探討具有導電梯 度分佈的電滲流場(electroosmotic flow, EOF)之穩定性特性。理論模型 則假設在填滿稀薄二元電解液(dilute binary electrolyte solution)的兩 無限平板間，通入一水平方向的電場，電雙層(electric double layer, EDL)內的流體受電場作用而產生邊界滑移速度，進而推動流體運動。

由於流場內部存在導電梯度，使得自由電荷與電力呈現不均勻的分佈。

當外加電場持續增加超過某一臨界值時，流場便會開始變得不穩定。

為了瞭解此滑移速度在此系統中所扮演的角色，我們分別計算有EDL 及無EDL兩個情況。研究結果發現，在低導電度梯度時，EDL所產生 邊界滑移速度的確增強了此流體系統的不穩定作用。然而在高導電梯

IX

度時，在EDL中的則塔電位(zeta potential)增強時，似乎對系統有更穩 定的影響。本研究藉由大量數值結果的呈現，嘗試解釋其中重要的物 理機制，成果將有助於電流體穩定現象的進一步了解。

關鍵辭:電滲透流，穩定性分析，水平電場，導電梯度，則塔電位

X

Instability of electrokinetic flow in a thin fluid layer between two parallel plates with an electrical

conductivity gradients

Advisee: CAI, YING-TING Advisor : Falin Chen, Ph.D.

Graduate Institute of Applied Mechanics College of Engineering

National Taiwan University

Abstract

In the microfluidic device, it is the most effective way that the fluid is driven by the electric force is better than pressure. The

electrohydrodynamics (EHD) is integrated electricity and the fluid mechanics have the interesting phenomenon. Specially the DC causes the

free charge to accumulate when the fluid has conductivity gradient, then occurs unstable convection phenomenon. This kind of stabilization is most

often applied to the micromixer. In the present the Micro Electro Mechanical Systems (MEMS) technology is quite important research subjects. This paper mainly considers electroosmotic flow (EOF) of the conductivity gradient distribution to use the linear stability analysis method.

The model use dilute binary electrolyte solution between two infinite plates passes over a horizontal electric field. The electric field makes the electric

double layer (EDL) in fluid to produce the boundary slip velocity then impetus fluid motion. As a result of the flow field has electric conduction

gradient inside to cause the free charge and the electric force are the

XI

non-uniform distribution. When the applied electric field continues to increase the critical value then the flow will become unstable. In order to

understand that this slip velocity which acts in this system, we calculate separately to have EDL and no two situations. we discover when low conductivity gradient , EDL had the boundary slip velocity indeed to strengthen system unstable. However, at high conductivity gradient when EDL of zeta potential is stronger that system is more stable. This paper use

massive value result to explain the important physical mechanism, the achievement will be helpful to the electrohydrodynamic stability a little

understanding.

Keywords: electroosmotic flow, stability analysis, horizontal electric field, conductivity gradient, zeta potential

XII

參數說明

英文字母

C 導電離子濃度(concentration of ions) mol m⋅ ^-3

D d dz/ ，微分運算子 －

d 兩平板間距離(depth of two plates) m E0 初始外加電場( imposed electric field) _{V m⋅} ^-1

E 電場向量(electric field vector) V m⋅ ^-1

fe 電物體力(electric body force) N

F 法拉第常數(Faraday constant) C mol⋅ ^-1 H Δσ σ/ 0，無因次導電梯度(dimensionless

conductivity gradients)

－

Jf σE+vρ_f，電流密度(electric current density vector) A m⋅ ^-2 keff 離子有效擴散係數(effective diffusivity) m s²⋅ ⁻¹

k z方向單位向量 －

k 無因次波數(dimensionless wave number) －

n 則塔電位(zeta potential)修正指數 －

n 法線單位向量(normal vector) －

p 壓力(pressure) _{N m⋅} ^-2

Q σ σ₀/ _r，無因次底板導電度 －

Rae εE²₀d²/μk_eff，電瑞利數(electric Rayleigh number) －

Re U d_ev /ν^，雷諾數(Reynolds number) －

Rv ψ_Rae^{1/ 2}，無因次電滲速度(dimensionless

electroosmotic velocity) －

s 擾動波成長率(growth rate) － Sce

ν

XIII

t 時間(time) s

t 切線向量(tangential vector) －

u x 方向流體速度 m s⋅ ⁻¹

Uev εE d₀² /μ ，電黏滯速度(electroviscous velocity)

v y 方向流體速度 m s⋅ ⁻¹

v 流體速度向量(velocity vector) _{m s⋅} ⁻¹

w z 方向流體速度 m s⋅ ⁻¹

Z 原子價數(valence number) －

希臘字母

ε 介電係數(permittivity) ^{C m⋅ -1⋅V-1}

ξ 則塔電位(zeta potential) V

θ 角度(angle) －

LD 德拜長度(Debye length) m

μ 動力黏滯係數(dynamic viscosity) _{N s m⋅ ⋅} ^-2 ν 運動黏滯係數(kinematic viscosity) _{m s}²_⋅ ^-1 ω 離子遷移率(mobility) mol s kg⋅ ⋅

ρ 密度(density) kg m⋅ ^-3

ρf 自由電荷密度(charge density) _{C m⋅} ^-3

σ 導電度(conductivity) _{S m⋅} ⁻¹

τ ^{σ0d2/εν ，}黏滯鬆弛時間(viscous relaxation time) 與電荷鬆弛時間(charge relaxation time)比值 ^－

ϕ 無因次電位(dimensionless potential) －

ψ

(

)

XIV

上標

￣ 基態解

* 無因次化

′ 微擾動量

^ 微擾動量之z方向分量

下標

o 位於下平板

d 位於上平板

r 參考值

x x方向分量

y y方向分量 z z方向分量

1

第一章 緒論

1.1 研究背景

電流體力學(Electrohydrodynamics)是一門整合流體力學與靜電 學的跨領域學科。而電流體的穩定性分析，近來受到微機電技術的需 求，如生物反應晶片或化工過程之微型化技術，凡舉化學分析、生醫 系統、微幫浦、注射器、成分分離器、混合器等等，其應用層面相當 廣泛。由於微機電製程的進步，現更可將多樣系統整合於單一晶片上，

即「System on a Chip」，得以降低製造成本，並縮小系統的整體尺寸，

以提升可攜度[1]。亦由於尺寸的縮小，許多在維度較大的系統中被 忽略的效應則有更重要的影響，包括電磁場效應、熱傳遞、流道平滑 度…等等。

在微流體裝置的應用上，主要有微混合器與微分離器兩種。首先，

在微流體混合器中，由於流場之雷諾數通常都相當低，流體多以層流 方式流動，致使流體之間的混合通常僅能藉分子擴散來達成。增進微 流體混合器效率所採用的方式一般可分為被動式與主動式兩種：被動 式混合係以在流道中設置擋塊，或改變流道幾何形狀等方式促使流體 混合，此種方法在混合器製造後即無法變動，於使用上較缺乏彈性。

而主動式混合係以加入電力、磁力、熱力…等作用力，或是利用旋轉、

震動等運動作為增進混合的機制。但需要廣泛的研究成果來加以應用，

且其中牽涉較高深學理，需要專才解決。目前以電力來引發流體不穩 定來製造對流與混合[2,3]為其主要發展技術且以電力控制亦有較容 易操控的優點。

除了微流體混合器的研究外，另一個重要應用是在所謂的『等電

2

位點聚焦電泳(isoelectric focusing, IEF)』技術。該技術源於生物技術 領域，主要被利用來分離蛋白質或高分子混合物之未知成份。在等電 位點聚焦電泳實驗操作中，若增加兩電極板間的電位差，通常可以加 速蛋白質分離的速度；但另一方面，電位差的增加亦會使得電泳力越 大而超越黏滯力的束縛，導致流體不穩定流動。此時乘載蛋白質的電 極緩衝液亦產生擾動，緩衝液黏性下降，分子自由擴散加劇，進而降 低PH 值分佈導致分離失效。因此，電壓施加的強弱，有其極限存在。

一般而言對於生物反應晶片來說，希望設計能增加系統的不穩定性混 合，來加速反應；但對於電泳分離技術而言，卻希望能提升系統的穩 定機制。故依據需求來設計流體的穩定性。

在微流體系統中，以驅動方式來說，主要有以下兩種方式：壓力 驅動與電動力驅動。在微流道中以壓力梯度驅動流體是相當沒有效率 的方式。且微流道的結構強度不高，無法承受太大的壓力梯度。因此，

電場驅動，是較合適的選擇。在微流體分析系統中，受處理的流體常 為可解離的溶液，或是具導電性的流體，若將流道的壁面設計成電極 板，便可藉由外加電場的方式驅動流體[4]。

1.2 文獻回顧

由前節提到，在生物晶片的研發上，為了加速物種混合效率，我 們希望能趨向設計系統產生不穩定；對IEF 則希望提升系統的穩定機 制。為了清楚了解這些物理機制的成因與其耦合作用，已有學者們先 後提出各種理論模型來分析與預測。

在 1965 年 Taylor & McEwan[5]研究在兩層不相容的液體界面做 過穩定性分析，1967 年 Melcher & Firebaugh [6] 率先研究弱導電的液

3

體在加入溫度梯度後所引起的導電梯度變化對穩定性的影響。

1970 年 Michael & O’Neill[7]研究一通入垂直方向的電場，並在 兩個導電流體中隔著一層非導電體平板的穩定性分析。以上的研究在 在都顯示，在沒有接觸面間的張力影響下，靜止的流體也會產生擾動 現象。Hoburg & Melcher [8] 與 Hoburg [9] 分別在 1976 年及 1977 年 則考慮兩個彼此可互相混合的液體介面存在一線性的導電性梯度 (conductivity gradient)分佈，並探討在分別受到切平面方向以及垂直 方向之電場作用下的穩定性行為。

Rhodes[10]在 1989 年也將兩可互融的流線噴流做分析，說明在

沒有受接觸面間張力的阻力下，原本圓柱狀的噴流會變成帶狀。

Melcher & Taylor[11]也做了許多加入電力驅動相關的研究。Baygents

& Baldessari[12]和 Hoburg & Melcher[13]有提到當電場與導電梯度相 互作用時，靜電荷也可在流體內部作用，而不僅限於表面上。

而 Baygents & Baldessari[12]，在假設不流動的流場，加入導電梯 度的變化，此導電梯度為離子濃度擴散造成。當黏滯力的反應時間遠 比電解液受電場作用的反應時間大，其導電度的統御方程式可以簡化 幾乎近似於擴散方程式。此項觀念亦可由Melcher[14]所提的觀念相 通。而此項相較於Hoburg & Melcher[15]有所不同。因在 Hoburg &

Melcher[15]的研究中所分析為共向性(colinear)的流場及導電梯度，雖 然其導電梯度一樣也是離子濃度擴散造成，但不同的是其導電梯度在 穩定性分析推導過程中會被消除。其消除的原因為當k_eff /ν <<1 時，

流體擴散所需的時間比黏滯力的反應時間長。

雖然k_eff /ν <<1，但對影響流體的運動還是有其重要性，當電流體

4

運動時，其離子的擴散效應會使得離子濃度均勻，衰減導電梯度的產 生，進而穩定流體。因為當離子的擴散效應減小，也就代表因離子擴 散所產生的導電梯度也減小。而當導電梯度減小，也就代表電場所能 對流體的作用力也減小，因此也就間接造成流體穩定。

再討論，如果擴散效應大於電泳效應，則流體將停止流動。因為 電場無法造成導電梯度讓帶電離子漂移。且流體的黏滯力在此扮演穩 定流體的功用，阻礙流體運動。假如k_eff /ν 假設趨於零，也就代表忽 略擴散效應的影響，等於只有討論外加電場對流體做穩定性分析。

在 Lin[16]文中討論在微流道中施加兩層均勻的電解液，上層導 電度較高，下層的導電度較低，其導電梯度僅存在於介面中。當施予 外加的水平電場時發現，所施予的電場越高，則其介面擾動的程度越 大。若固定電場，也發現當兩層的導電度差異越大，其擾動也越劇烈。

此一現象常用於，微流道的快速混合機制。

典型的微流道尺度大約在 10 microns 或更大，而杜拜長度(Debye lengths)大約在 10nm 等級，尺度相差 1000 倍。因此，將電雙層(electric double layer, EDL)在邊界的效應考慮為邊界滑動條件，其滑動速度與 局部電場和黏滯力有其比例關係。而Santiago[17]其物理模型中說明 如何考慮邊界EDL 效應，將 EDL 的電滲速度效應，考慮為邊界的滑 動速度效應而滑動邊界考慮在EDL 擴散層末端。Chen[18]更進一步 將EDL 中的電滲速度，其則塔電位(zeta potential)的部分，考慮與導 電度的影響關係。

1.3 研究動機與目的

5

可從實驗了解在厚長比甚小的微流道中施加與流動方向垂直的 導電梯度時，發現當施予流動方向的電場時，流體會隨著所施加的電 場大小產生顯著不穩定的擾動現象，且電場越大其不穩定現象越明顯 亦越快。當考慮導電梯度越高時，亦有相同的現象。因此，常將此現 象用於微混合(micromixer)裝置上。另外也發現若流道厚度尺度很小 時，反而會有使流體趨向穩定的現象。因此要如何控制厚度尺度才能 達到最有效的混合效率亦是許多學者研究的方向。

而由 Lin[16]的模型中，討論在微流道中考慮邊界 EDL 效應下，

兩層相異導電度的電解液予以水平電場作用，分析此系統的混合機制。

並說明當外加電場越大，則混合機制越明顯。

此系統考慮在常溫下，忽略溫度作用，討論其流體穩定性分析；

而其中較少討論其縱向模態的變化，且其所討論的導電梯度為在各層 中為均勻僅在交界面產生梯度變化。因此，本研究的重點放在將系統 考慮為平行板系統，並藉由給予上下板離子濃度經由擴散作用而產生 線性的導電梯度變化下。分別討論其橫向模態與縱向模態的線性穩定 性分析。

在本研究的物理模型，亦將討論以下狀況的線性穩定性分析: (1).

考慮EDL 在邊界的影響。 (2).考慮當 EDL 尺度遠小於板高，忽略 EDL 邊界效應的影響。

1.4 研究方法

本研究之物理模型可簡述為:兩具導電度差之平行板間充滿導電 溶液在常溫下，通入一平行於平板的均勻水平電場。考慮邊界上的電 雙層(electric double layer, EDL)影響，在本研究中 EDL 的影響主要考

6

慮在邊界上電滲流所造成的滑動速度，即電滲速度(electroosmotic velocity)所造成的邊界滑動[18]。

在本研究中，主要考慮則塔電位(zeta potential)、電場對電滲速度 的作用。在則塔電位中，依據Lin[16]考慮導電度對則塔電位有修正 關係。所以當給予上下板不同導電度時，將會引發不同的則塔電位，

進而在邊界上產生不同的電滲速度。而邊界導電度的給予，可以根據 Melcher[14]，當電荷鬆弛時間(charge relaxation time)<<黏滯擴散時間 (viscous diffusion time)時，可將電解質濃度直接線性正比於導電度分 佈。因此當給予上下板離子濃度，亦即給予邊界導電度。

本研究所使用的統御方程式除了連續方程式與動量方程式之外，

還使用電荷守恆方程式，考慮電場、導電度與自由電荷之間的關係，

導電度控制方程式。線性穩定性分析的流程，除了先確定物理模型，

使用的假設與統御方程式。之後由上述統御方程式內的各未知數配合 邊界條件解出基態解，再根據基態解選擇適合的特徵尺度將統御方程 式與基態解及邊界條件無因次化。利用微擾法對無因次化的方程式進 行線性穩定性分析，並以正規模態展開。運用數值分析方法求得未知 數，找出線性穩定性的中性曲線。

在數值方法中，本研究使用謝比雪夫配置方法，配合邊界條件，

求得在各種參數的組合下，要使流場發生不穩定現象所需的臨界波數 及臨界電場強度。

7

第二章 理論模式

2.1 模型建立與假設

本文所探討的模型如圖 2.1 所示。考慮兩平行電極板間充滿導電 性液體，及邊界上產生電雙層(EDL)效應。在本文 EDL 的效應主要討 論在邊界上產生電滲速度(electroosmotic velocity)即對系統產生邊界 滑動的影響。在忽略溫度效應的系統下，於兩平板邊界上施予相異的 離子濃度，則系統會因擴散作用，產生導電梯度效應。在此情況下，

加入均勻水平電場的作用時，分析其流場穩定性。

圖 2.1 模型示意圖

本問題會使用一些合理的假設。(1)在 x 及 y 方向為無限延伸平 板，則基態速度僅為與z 軸有關的 x 向速度 ( )u z ；(2)流體為不可壓縮，

表忽略流體密度在時間與空間的變化；(3)考慮等溫的稀薄電解溶液，

表系統所使用的電解液大部分為水所組成，且忽略溫度效應對系統的 變化，如電場作用後所產生的焦耳熱(joule heating)對流體的影響 亦予忽略；(4)流體為線性材料，介電係數、黏滯係數、離子有效擴

( )z σ

σ0

σd

E

x z

y d

剛性邊界 電雙層(EDL)

8

散係數視為常數；(5)忽略介電泳效應、電伸縮效應及黏滯消耗等假 設來簡化問題。

2.2 統御方程式

根據模型的簡化假設與參考相關文獻後，本研究所使用的統御方 程式共有下列4 條:

1. 連續方程式

在不可壓縮的狀態下，可由質量守恆關係將方程式簡化為下式。

∇⋅ =v 0 ...(2-1) 其中， v 含有三個方向的速度項v=( , , )u v w 。

2. 動量方程式

為慣性力、壓力、黏滯力、物體力(body force)的平衡方程式。因 系統忽略重力、磁力等物體力(body force)，且帶電離子受到電場作用 的影響產生轉移而造成離子間的碰撞。因此，Navier-Stokes equations 可修正成下式:

2 e

D p

ρ Dtv = −∇ + ∇ +μ

v f ...(2-2) 其中，f_e為電物體力(electric body force)，在下式中，右式第二項 介電泳力與第三項電伸縮因介電係數假設為常數而消去，故只考慮第 一項的庫倫力。

2 2

1 1

2 2

f

ρ ε ρ ε

ρ

⎡ ⎛∂ ⎞⎤

= + ∇ + ∇ ⎢⎣ ⎜⎝∂ ⎟⎠⎥⎦

fe E E E

...(2-3) 因此可以簡化動量方程式為

9

2 f

D p

ρ Dtv = −∇ + ∇ +μ ρ

v E

...(2-4)

3. 電荷守恆方程式

而由 Gauss’ law 的電荷守恆得知

⁰

f

t f

ρ

∂ +∇⋅ =

∂ J

...(2-5) 在電阻性物質中，(2-5)式的電流密度J_f 可以(2-6)式表示為導電度

σ 與電場作用產生的歐姆電流及電荷移動造成的電流 J^f =σE+vρ^f ...(2-6)

將(2-6)式帶入(2-5)式，並考慮連續方程式即可化簡成下式:

^{( )} D f

Dt

ρ = −∇⋅ Eσ

...(2-7)

4. 質傳方程式

因本研究的工作流體為離子流體，流體的導電性與離子濃度有關 係，而其導電性會隨不同位置與時間改變。導電度的變化與區域性的 帶電粒子密度有相當大的影響。

由 Melcher[14]可知，當

0

ε σ <<

0

d

wE <<^d2

ν 。在此前提下才可獲得 近乎擴散方程式關係的導電度的控制方程式:

^D K_eff ² Dt

σ = ∇σ

...(2-8) 其中，K_eff為由布朗運動影響所產生的離子有效擴散係數；

0

ε

σ 為

電荷反應時間(charge relaxation time)；

0

d

ωE 為粒子遷移所需時間 (particle migration time)；^d2

ν 為黏滯擴散所需時間(viscous diffusion

10

time)。即堆積電荷的反應時間要比帶電離子受電場作用的漂移時間要 短。而帶電離子受電場作用的漂移時間也要比黏滯力作用在流體粒子 的反應時間短。此時才可假設質傳方程式如(2-8)所示。

2.3 邊界條件

1. 導電度

當給予離子濃度可由關係式獲得導電度[16]。

2 2

Z F C

σ ω= ...(2-9) 其中，^Z為離子的原子價(valence number)；F為法拉第常數(Faraday constant)，一莫耳電子所帶的電荷；C為離子的莫耳濃度(concentration)；

ω 為離子遷移率(mobility)。

因此給予上板濃度C0，下板濃度^Cd，即給予邊界導電度。

⁰ (0) ( )d _d

σ σ

σ σ

⎧ =

⎨ =

⎩ ...(2-10)

2. 電場

在電場的邊界條件中，考慮電場不可穿透條件 ϕ 0

∇ ⋅ n = ...(2-11) EDL 即電雙層(electric double layer)為緊密層(stern layer)與擴散 層(diffuse layer)的簡稱。由圖 2.2 將 EDL 放大來看，緊密層因距離帶 電邊界最近，對離子的吸引力最強，產生密不可分的不動層。超過緊 密層之後，即是擴散層。在擴散層中因帶電粒子逐漸吸附中和而使其

11

吸引力減弱，使帶電粒子可脫離束縛移動。距離緊密層越遠越易脫 離。

而緊密層與擴散層中間的界面為剪力面(shear plane)。在擴散層 末端因電性中和的作用，趨於電中性，故在擴散層末端的電勢能趨於 電中性。而擴散層末端與剪力面的電勢能差即為則塔電位(zeta

potential )；擴散層的厚度亦稱為德拜長度(Debye length)。

系統考慮的邊界即位於擴散層末端的邊界 0

z LD

z ϕ

=

∂ ≈

∂ 。

圖 2.2 EDL 示意圖

3. 速度

考慮邊界不可穿透的條件。可知法線方向速度為 0

w= ...(2-12) 而切線方向的速度中，在邊界上產生電滲速度(electroosmotic velocity )。在此使用 Helmholtz-Smoluchowski formulation 來計算電滲 速度。

shear plane

zeta potentail

Z

stern layer diffuse layer

ξ

ϕ

LD

Solid

Debye length

12

εξ μ

⋅ =− ⋅

v t E t

...(2-13) 其中，(2-13)式中的則塔電位ξ ，考慮局部離子濃度的影響(即導 電度的影響)並作(−n)次方修正[16]，即下式所示:

^{( )}

n

r r

ξ σ

ξ σ

= −

...(2-14) 而考慮外加水平電場的作用，我們令電場為

^{E=E0ex− ∇ϕ}...(2-15) 上式代表，在 x 方向施加外加電場(E )₀ 與系統內部反映電場

(−∇ϕ)的疊加關係。將(2-14)式帶入(2-13)式可獲得邊界切線速度為:

( ) ⁿ

r r

εξ σ μ σ ⁻

⋅ =− ⋅

v t E t

...(2-16)

2.4 基態解

可由上述所說明的假設、統御方程式、邊界條件中，求出變數的 解，即為此模型的基態解。此為系統在穩態且忽略時間成長，亦不考 慮擾動變化時的狀態。以下在變數上標加上[－]代表基態解。

計算基態解時，所使用的合理假設如下:

(1)定常(steady state)，即 ^{( )} 0 t

∂ =

∂ ...(2-17) (2)導電度無 x 方向的變化，即 0

x

∂σ =

∂ ...(2-18) (3)電荷密度無 x 方向變化，即 ^f 0

x ρ

∂ =

∂ ...(2-19) (4)完全發展流

(5)因為 x-y 無限延伸的平行板，故令v=[ ( ),0,0]u z ...(2-20)

13

依據 2.1 節、2.2 節、2.3 節，所列出的假設、統御方程式及邊界 條件。求出系統中各未知變數的基態解。

1. 導電度

當給予上下板不同離子濃度時，因考慮黏滯擴散時間>>電荷鬆弛 時間。故電場對離子濃度的變化相當快，使得由離子濃度擴散所形成 的導電梯度為幾近線性變化。

        z ₀ d

σ = ^Δσ +σ

...(2-21) 2. 電場

由電場不可穿透的邊界條件(2-11)帶入電荷守衡方程式(2-7)可知，

電位基態解ϕ 為常數，則基態電場僅為外加電場E₀。

E = E₀e_x ...(2-22)

3. 自由電荷密度

在(2-22)式中得知，基態電場僅為外加電場E₀，系統內部並無反 應電位ϕ，由高斯方程式可知，將電場取散度即可得到自由電荷密度。

代表自由電荷密度在基態電場為電中性的狀態，系統內部並無電泳力 的發生。

f 0

ρ = ...(2-23) 4. 壓力

將上述(2-22)、(2-23)電場及電荷密度的基態解代入(2-4)式動量方 程式中，可解出本系統的基態壓力均為等壓，但在此設為零。即本系 統並不給予外加壓力驅動流體，或不考慮壓力對基態流場的影響。

14

P=0...(2-24) 5. 速度

將邊界速度條件(2-12)式、(2-16)式帶入動量方程式(2-4)式中，並 將上述電場、自由電荷密度及壓力的基態解(2-22)式、(2-23)式、(2-24) 式同時帶入動量方程式(2-4)中即可求出基態流場。我們已可由(2-23) 可知基態的電荷密度為電中性。外加電場並無法對系統產生電泳力的 作用。因此系統的基態流場是由邊界的滑動速度所造成的流場，近似 庫第流(Couette flow)因邊界移動所造成流場。不同的是，在本系統 的邊界滑動速度，為因導電度差異所造成邊界有不同的電滲速度，使 上下板因有不同的滑動速度來驅動流場。亦因系統無電泳力的作用，

使得流場為線性變化。

(0) u d^{( )} u⁽⁰⁾

u u z

d

= + −

0

0 0

( ) ( )

E ( )

n n

d

r n r r

r

d z

σ σ

ε ξ σ σ σ

μ σ

− −

−

⎡ − ⎤

⎢ ⎥

⎢ ⎥

= − +

⎢ ⎥

⎢ ⎥

⎣ ⎦...(2-25)

2.5 系統無因次化

我們通常會將統御方程式、邊界條件、基態解等方程式無因次化。

即去除單位絕對大小的影響，僅考慮階數(order)的大小以便合理的看 出變數的變化。

在本節，我們由 2.4 節解出的基態解選出系統的各個特徵量，將 變數上標加上＊來代表無因次變數。

15

表 2.1 特徵尺度表

特 徵 尺 度 說 明 符 號

長度 兩板間高度 d

速度 電黏滯速度 Uev

時間 －

ev

d U

電位 － E d0

導電度 兩板的導電度梯度 Δσ

壓力 － ^Uev

μ d

自由電荷密度 考慮導電度梯度的尺度 ⁰

0

E d

ε σ

σ Δ

其中，採用平行板高作為長度尺度。由 Hoburg & Melcher[15]及 Lin[16]，速度特徵尺度選用電黏滯速度U_ev。

U_ev ε^E20d

= μ

...(2-26) 此電黏滯速度為當電荷受電場作用移動時，在下游會累積電荷。

造成一反電位的產生，阻礙流體的流動，(2-26)式為考慮黏滯力與電 力平衡所求出。而自由電荷密度的特徵尺度採用高斯定理，對電場取 散度並考慮導電梯度的比例影響。

下列各式為無因次化的變數關係:

E0

∗ = E E

...(2-27)

16

d

∗ = x x

...(2-28) σ ^{σ σ0}

σ

∗ −

= Δ ...(2-29)

ev

t t d U

∗ =

...(2-30)

0 0

E

f f

d

ρ ε ρ σ

σ

∗ = Δ ...(2-31)

由高斯定理，ρ_f = ∇⋅ E(ε )

...(2-32) 把(2-32)帶入(2-31)中，將電荷密度化為電位的關係式:

ρ _f ^{σ0 σ0 2}ϕ

σ σ

∗ = ∇ ⋅∗ ∗ = − ∇∗ ∗

Δ E Δ

...(2-33)

ev

p p

U μ d

∗ =

...(2-34)

Uev

∗ = v v

...(2-35) 將 2.4 節(2-21)~(2-25)的基態解依(2-27)~(2-35)的關係式無因次化 成下列各式:

σ^∗ =z^∗...(2-36)

∗ = _x

E e ...(2-37)

f 0

ρ^∗ = ...(2-38)

0

p^∗ = ...(2-39)

0 0

0

( ) ( ) ( )

E

n d n n

r

r r r

u z

d

σ σ σ

ξ

σ σ σ

∗= − ⎧⎪⎨ − +⎡⎢ − − − ⎤⎥ ∗⎫⎪⎬

⎪ ⎣ ⎦ ⎪

⎩ ⎭

{ ( ) }

v

1 (H 1)

R

n n n

Q⁻ ⎡ Q ⁻ Q⁻ ⎤z^∗

= +⎣ + − ⎦

...(2-40)

17

其中，R_v為無因次電滲速度與電黏滯速度比；Q為無因次底板導 電度；H 為無因次導電度差。

R_v ^{ev E0}

eo r

U d

U ξ

= =

− ...(2-41) ⁰

r

Q σ

=σ

...(2-42)

0

H σ

σ

=Δ

...(2-43) 將 2.2 節(2-1)、(2-4)、(2-7)、(2-8)式的統御方程式依(2-27)~(2-35) 式的關係式無因次化成下列各式。(2-44)式為連續方程式，(2-45)式為 動量方程式無因次化。

∇ ⋅^∗ v^∗ =0...(2-44)

²

0

1 ( )

Re p ^f

t σ ρ

σ

∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗

∗

∂ + ⋅∇ = −∇ + ∇ +Δ

∂

v v v v E ¹ ( ² ( ² )( ))

Re ^{∗ ∗}p ^{∗ ∗ ∗}ϕ^{∗ ∗ ∗}ϕ

= −∇ + ∇ v + ∇ − + ∇e_x

...(2-45) 其中，Re 為雷諾數(Reynolds number)，為黏滯力與慣性力的比值。

2 2

E0

Re U d^ev ε d ν μν

= =

...(2-46) 由下列關係式，將電場(2-15)式與自由電荷密度(2-33)式以電勢能 函數取代，並代入電荷守衡控制方程式(2-7)式。

2 0 2

0 0

( ) ( ( )) (1 ) ( )

ev

D d

Dt U

σ σ σ

ϕ σ ϕ σ ϕ

ε σ σ

∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗

∗

⎡ Δ Δ ⎤

∇ = − ⎢ ∇ ⋅∇ − + ∇ + + ∇ − + ∇ ⎥

⎣ e^x e^x ⎦

(H ( )) (1 H ) 2

Re

τ ^{⎡ ∗}σ^{∗ ∗ ∗}ϕ σ^{∗ ∗}ϕ^∗⎤

= − ⎣ ∇ ⋅ − + ∇e^x + + ∇ ⎦...(2-47) 其中，τ τ= _viscous/τ_charge =σ₀d²/εν 為電荷反應時間與黏滯擴散反應時

間比。而τ_viscous =d²/ν ，τ_charge =ε σ/ ₀，分別為 2.2 節所述的黏滯反應時

18

間與電荷鬆弛時間。當σ₀ ≈10 S/m⁻⁴ ，d ≈10 m⁻³ ，ε ≈10 C/m V⁻⁹ ⋅ ，ν ≈10 m /s^{−6 2} ，

則τ ≈10 ~ 10^{5 6}。表示(2-47)中，左式電荷密度隨時間的變化項可近乎可

以忽略其影響。

而質傳方程式亦可無因次化如下式所示:

^{1 2}

eRe D

Dt Sc

σ_∗^∗ = ∇^∗σ^∗

^{1 2} Rae ^∗ σ^∗

= ∇

...(2-48) 其中，Sc_e為電施密特數(electric Schmidt number)，為擴散係數與 粘滯係數的比值，在本研究為代表離子擴散的效應。Ra_e為電瑞利數 (electric Rayleigh number)，在本研究為代表電場的效應。

_e

eff

Sc k

= ν

...(2-49) _e ^{ev E20 2}

eff eff

U d d

Ra k k

ε

= = μ

...(2-50) 本文亦將 2.3 節(2-10)、(2-11)、(2-12)、(2-16)式的邊界條件依(2-27) 式~(2-35)式的關係式無因次化成下列各式:

^{(0) 0} (1) 1 σ

σ

∗

∗

=

= ...(2-51) _{z 0,1} 0

z ϕ

∗

∗

∗ =

∂ =

∂ ...(2-52) ^{w z 0,1 0}

∗

= =

...(2-53)

^0,1 _{v 0,1} 1 ( )

R

n

z r z

σ σ

∗

∗

∗ ∗ − ∗

= =

⎡ ⋅ ⎤ = ⋅

⎣v t⎦ E t

...(2-54)

最後可以由Ra_e&R_v的關係得到一無因次參數

ψ

19

1 1

2 2

v

R 1 ^eff _{e e}

r

k Ra Ra

μ ψ

ξ ε

⎛ ⎞

= − ⎜ ⎟ =

⎝ ⎠ ...(2-55) 其中，

ψ

ψ

1

1 eff 2 r

μk

ψ ξ ε

⎛ ⎞

= − ⎜ ⎟

⎝ ⎠ ...(2-56) 將(2-55)式的關係帶入(2-25)後得到:

( )

{ }

( )

{ }

v

1 2

1 (H 1)

R

1 (H 1)

n n n

n n n

e

u Q Q Q z

Q Q Q z

ψRa

− − −

− − −

⎡ ⎤

= +⎣ + − ⎦

⎡ ⎤

= + ⎣ + − ⎦

...(2-57)

2.6 線性微擾化

為了瞭解基態流場的穩定特性，將 2.4 節中的基態解加上微小擾 動量(infinitesimal disturbance)帶回統御方程式與邊界條件，並忽略微 擾量的乘積項即高階項，可得線性擾動方程式(linear perturbation equation)。且在此後章節皆省去上標(＊)來代表無因次參數，上標(′) 代表微擾量。並設定如下變數:

( )^{,t = ( )z + ′}( )^,t

v x v v x ...(2-58)

( )z ( , )t

σ σ= +σ′ x ...(2-59) ( )z ( , )t

ϕ ϕ= +ϕ′x ...(2-60) ( ) ( , )

P=P z +P′ x t ...(2-61) 將上述擾動變數，帶入無因次化統御方程式(2-44)、(2-45)、(2-47)