杜昌隆, 夏威豪, 楊嘉杰, 李 捷

(武漢理工大學 船海與能源動力工程學院,武漢 430063)

0 引 言

微米通道流體系統(tǒng)在工程、生物和醫(yī)學領(lǐng)域的應(yīng)用不斷增加,使微納米領(lǐng)域成為近十年來最具吸引力的流體力學研究領(lǐng)域之一．早在1809年,Reuss在對多孔滲水層實驗研究時發(fā)現(xiàn)了動電現(xiàn)象[1]．Rice等[2]研究了細長圓柱形毛細管中的動電流,討論了在雙電層(EDL)效應(yīng)的影響下,外加電場和壓力產(chǎn)生的電黏、流動電勢和電滲等現(xiàn)象．劉浩等[3]采用格子Boltzmann模型研究了T型微通道中不同冪律指數(shù)和黏度比對液滴形成的影響．Zhao等[4]分析了冪律流體在狹縫微通道中的電滲流(EOF),并考察了流動行為指數(shù)、雙電層厚度和電場的影響,用更通用的卡勞非Newton模型對電滲透遷移率進行了數(shù)值研究．王爽等[5]在微通道具有調(diào)制zeta電勢時,在電場力與Lorentz力混合驅(qū)動中得出流函數(shù)以及速度分布的解析解,并且調(diào)制的壁面電勢會產(chǎn)生一個垂直速度分量,從而導(dǎo)致渦旋的形成．上述研究僅限于簡單的非Newton流體模型,沒有表現(xiàn)出彈性特性．Afonso等[6]報道了使用Phan-Thien-Tanner(PTT)模型和有限可擴展非線性彈性FENE-P模型對微通道內(nèi)兩個同心圓柱體之間的黏彈性流體EOF的解析解．在早期工作的基礎(chǔ)上,Afonso等[7]也研究了不對稱zeta電位微通道中黏彈性流體的EOF．Sousa等[8]推導(dǎo)了考慮壁面耗盡層的解析解．對于PTT流體的純EOF,Park等[9]首先推導(dǎo)出Helmholtz-Smoluchowski速度,為計算微通道內(nèi)的體積流量提供了一種簡單的方法,同時也對微通道中壓力-電滲混合驅(qū)動力的黏彈性流體EOF進行了研究．目前,已有不少研究人員將重點放在微通道中壓力驅(qū)動和電滲透流動的結(jié)合上[10-13]．Dutta等[14]研究了直流道中電滲與壓力混合流的一維速度解析解,研究發(fā)現(xiàn)在雙電層效應(yīng)的影響下,可以在流道壁面上施加H-S 速度作為滑移邊界條件．Ferrás等[15]使用簡化的PTT模型提出了在微通道中,電滲和壓力驅(qū)動的組合流動下的黏彈性流體的解決方案．高峰等[16]與羅艷等[17]數(shù)值模擬了二維矩形流道的壓力-電滲驅(qū)動流．由于黏彈性流體的彈性存在,非Newton流體的 EOF具有時間依賴性,即使在低Reynolds數(shù)(Re)下也表現(xiàn)出不穩(wěn)定性．Bryce和Freeman[18]首次報道了當施加的電場超過閾值時,通過2∶1∶2的微尺度收縮/膨脹,聚丙烯酰胺(PAA)溶液的電彈性不穩(wěn)定性．Song等[19]通過T形微通道對黏彈性聚環(huán)氧乙烷(PEO)溶液的EOF中的彈性不穩(wěn)定性進行了實驗和數(shù)值研究,結(jié)果表明,不穩(wěn)定開始的閾值電場高度依賴于 PEO 濃度．

上述關(guān)于EOF不穩(wěn)定性的研究大多考慮純EOF,雖然電滲是一種合適的驅(qū)動方式,但由于流場中存在機械閥或由于入口和出口之間的水頭差,壓力梯度可能隨著電滲沿微通道出現(xiàn),同時生物體內(nèi)物質(zhì)運輸過程中也存在壓力梯度．此外在某些特殊應(yīng)用中,可以采用電滲和壓力混合驅(qū)動來實現(xiàn)更加精確的流體控制,以實現(xiàn)工程中的實際要求,例如采用電滲壓力混合驅(qū)動提高石油采收率,分析由不穩(wěn)性產(chǎn)生的渦流對微流體混合的作用．因此,采用電場、壓力等驅(qū)動與控制的方式,研究生物流體在微通道內(nèi)流動與混合的問題,對于生物流體芯片的基礎(chǔ)研究、設(shè)計開發(fā)與其應(yīng)用具有十分重要的理論和實際意義．目前關(guān)于電滲壓力混合流的彈性不穩(wěn)定性相關(guān)研究還較少．針對EOF的研究主要是雙電層厚度、雙電層電勢對EOF流型的影響,關(guān)于壓力梯度對電滲-壓力驅(qū)動流速度分布影響的研究也少有見到．本文將基于Poisson-Boltzmann(PB)模型,采用電滲壓力混合驅(qū)動條件,分析了壓力和聚合物濃度對流場速度與渦流不穩(wěn)定性的影響．

1 模型、控制方程與驗證

1.1 幾何模型

在微通道內(nèi),帶負電的壁面從溶液中吸引帶正電的離子,從而形成雙電層．施加電壓時,雙電層擴散區(qū)域中的陽離子向陰極方向遷移,同時攜帶配位的水分子流動．從而溶液向負極凈流,即所謂的EOF．

假設(shè)流動發(fā)生在二維微通道中,并且假定:1) 流體為Oldroyd-B流體,為假象流體,是非Newton不可壓縮的; 2) 熱物性在整個區(qū)域內(nèi)恒定; 3) 流動完全發(fā)展; 4) 壁面zeta電位恒定; 5) 電解液是對稱的; 6) 離子是點電荷; 7) 流體的介電常數(shù)恒定,不受總場強的影響．微通道幾何形狀如圖1所示．

圖1 幾何結(jié)構(gòu)示意圖Fig.1 The structure diagram of the microchannel

1.2 控制方程

根據(jù)質(zhì)量守恒和動量守恒原理,在流場中任意地方的流體都一定要滿足連續(xù)性方程,通過連續(xù)性不可壓縮流體的假設(shè)條件,可以得到黏彈性流體流動的基本控制方程,包括質(zhì)量守恒方程(式(1))和動量守恒方程(式(2)):

?·u=0,

(1)

(2)

式中u為速度矢量,ρ和p分別表示流體密度和壓力,t為時間,f為任何外部體積力,τ′為額外應(yīng)力張量．為了模擬黏彈性流體流動,常用的方法是將額外應(yīng)力張量分離為溶劑貢獻張量(τs)和聚合物貢獻張量(τ),且τ′=τs+τ．為了得到一組閉合的方程,每個貢獻張量都需要一個本構(gòu)方程,本構(gòu)方程通常為

(3)

(4)

本文中f代表電體力,存在于雙電層中,雙電層的厚度用Debye長度表示:

(5)

式中,ε,k,T,C0分別為介電常數(shù)、Boltzmann常數(shù)、絕對溫度、體離子濃度,z,e分別為離子價態(tài)的絕對值和元電荷．

黏彈性流體選用Oldoryd-B本構(gòu)模型,其方程為

(6)

式中λ為黏彈性流體的弛豫時間．

本研究中所使用的流體是受電場作用的弱電解質(zhì)．在這種情況下,動量方程(2)中的f應(yīng)包含電場力:

(7)

其中E是電場強度,I為單位矩陣,ε=ε0εr表示溶液的介電常數(shù),ρe為電荷密度．

假設(shè)離子遵循 Boltzmann平衡,PB模型為

(8)

式中,F是Faraday常數(shù),ci,0是物種i的參考濃度,zi是物種i的電荷價,其中內(nèi)在電勢是ψ0．在不失一般性的情況下,其中固有電勢ψ0=0．式(8)的右側(cè)表示 PB模型的電荷密度．對于這個模型,唯一與電相關(guān)的未知數(shù)是兩個電勢ψ和ψext．此外,從式(8)可以看出,PB模型中的流量變量沒有影響(單向耦合)．為了增加式(8)的隱式性,其源項可以通過Taylor級數(shù)展開直至一階導(dǎo)數(shù),將方程轉(zhuǎn)化為

(9)

本文的具體參數(shù)設(shè)置見表1．

表1 模型部分參數(shù)Table 1 Numerical simulation parameters

本研究采用的邊界條件如下:

入口:?u·n=0,p=pin,E=E0,?ψ·n=0;

出口:?u·n=0,p=pout,E=0,?ψ·n=0;

通道壁面:u=0,?p·n=0,?E·n=0,ψ=ζ．

1.3 模型驗證

控制方程由RheoTool使用有限體積法進行數(shù)值求解,RheoTool是開源平臺OpenFOAM中實現(xiàn)的開源黏彈性EOF求解器．

(10)

H為通道半高,?p為通道x方向的壓力梯度,ε和ψ0分別為介電常數(shù)和通道壁面電勢．采用無量綱速度u/ush對比,其中

(11)

Ex為流動方向上的電場,η為溶液黏度．

數(shù)值仿真分析Γ=2的結(jié)果對比如圖2所示,由通道中間(x=0)沿y軸速度分布與Afonso等[7]的結(jié)果對比,可以發(fā)現(xiàn)數(shù)值結(jié)果與文獻結(jié)果吻合,驗證了采用的數(shù)值模型的正確性．

圖2 Γ=2速度分布驗證Fig.2 For Γ=2,the velocity distribution verification

2 結(jié) 果 分 析

2.1 渦流的不穩(wěn)定性

由于黏彈性流體在電勢超過一定閥值后會產(chǎn)生彈性不穩(wěn)定性[18],在縮放管的入口處產(chǎn)生渦流,并隨著時間波動．本節(jié)分析了濃度為300 ppm的聚合物溶液在外加電勢E=20 V時不同壓力下的流動情況．此時流體的聚合物黏度ηp=6.66×10-3kg/(m·s),弛豫時間λ=0.024 2 s．圖3和圖4分別展示了純EOF和擁有反向壓力梯度的混合驅(qū)動流下流線隨時間的變化趨勢．在施加反向壓力之后,渦流在任何時刻都比純EOF下大．但左右兩端儲池上下壁面附近的流線隨時間變化不大．在圖3(a)中,我們觀察到收縮通道入口處下方的渦已經(jīng)形成,上方的渦開始形成．隨后,入口處上方渦曲線彎曲加大(如圖3(b)所示),最終形成圖3(c)中入口處上下對渦一起存在,并不斷產(chǎn)生和消失．漩渦的大小和形狀在不同時期表現(xiàn)出顯著差異．在存在反向壓力時,這種情況出現(xiàn)轉(zhuǎn)變,渦流不再發(fā)生產(chǎn)生和消失的現(xiàn)象,渦流隨時間波動,且上下對渦一直波動存在．

(a) 0.60 s (d) 0.66 s

(b) 0.62 s (e) 0.68 s

(c) 0.64 s (f) 0.70 s圖3 p=0 Pa,Cp=300 ppm,E=20 V,不同時刻的流線Fig.3 For p=0 Pa,Cp=300 ppm,E=20 V,the streamlines at different moments注 為了解釋圖中的顏色,讀者可以參考本文的電子網(wǎng)頁版本,后同．

(a) 0.60 s (d) 0.66 s

(b) 0.62 s (e) 0.68 s

(c) 0.64 s (f) 0.70 s圖4 pout=3 Pa,Cp=300 ppm,E=20 V,不同時刻的流線Fig.4 For pout=3 Pa,Cp=300 ppm,E=20 V,the streamlines at different moments

如圖3所示,對于E=20 V,Cp=300 ppm的純電滲情況,渦流的寬度和長度與收縮微通道的高度幾乎相同．如圖4所示,當采用電滲壓力混合驅(qū)動時,渦流在反向壓力存在的情況下,渦流的高度大約是通道高度的兩倍,可以發(fā)現(xiàn)反向壓力讓渦流增大了,這是因為反向壓力的存在,縮管入口處的流體流入受阻,進而導(dǎo)致流體在入口處產(chǎn)生聚集．

圖5顯示了在兩種不同壓力下,縮放管入口處、中間、出口處3點的速度．兩種壓力下,3個位置在不同時刻的速度圍繞一個值上下波動．對于純EOF,即圖5(a),收縮微通道入口的平均速度為0.05 mm/s,標準偏差為0.02 mm/s;通道中心的平均速度為 0.626 mm/s,標準偏差為0.015 mm/s;收縮微通道出口的平均速度為0.35 mm/s,標準偏差為0.025 mm/s．對于pout=3 Pa,即圖5(b),收縮微通道入口、中心和出口的平均速度分別為0.103 mm/s,0.326 mm/s和0.160 mm/s,標準偏差分別為0.024 mm/s,0.022 mm/s和0.021 mm/s．在同一聚合物濃度下,壓力的存在使入口和中間處的速度波動都變大,而出口處的速度波動變化不明顯．入口處的平均速度在壓力存在的情況下是增大的,而中間和出口處的平均速度是減小的,通過圖3和圖4中入口處的渦流可以發(fā)現(xiàn),入口處的渦流在反向壓力存在下是增大的,入口處平均速度增大可能是渦流增大導(dǎo)致的．

(a) E=20 V,pout=0 Pa,Cp=300 ppm (b) E=20 V,pout=3 Pa,Cp=300 ppm圖5 3個不同位置的速度隨時間的變化Fig.5 The change of the velocity with time at 3 different locations

Newton流體和300 ppm的PAA溶液沿y軸的速度分布如圖6所示．使Newton流體的黏度與300 ppm的PAA溶液的總黏度相同,速度分布如圖6(a)所示．在混合驅(qū)動的情況下,對于Newton流體,混合驅(qū)動的速度等于純電滲驅(qū)動與純壓力驅(qū)動之和,這說明疊加原理在縮放管條件下是符合的．但是,在300 ppm的PAA溶液中,由于黏彈性流體彈性應(yīng)力的存在,疊加原理不再適用,如圖6(b)所示．純EOF與壓力驅(qū)動流在通道中間的速度相加和為0.958 9 mm/s,而混合驅(qū)動產(chǎn)生的速度為0.919 7 mm/s,相加的結(jié)果較混合驅(qū)動的速度更大,增大程度約為4.3%．

(a) 總黏度與300 ppm PAA溶液相同的Newton流體 (b) 300 ppm PAA溶液(a) The Newtonian fluid with the same total viscosity as that of the 300 ppm PAA solution (b) The 300 ppm PAA solution圖6 x=0,沿y軸的速度分布Fig.6 The velocity distributions along the y axis

2.2 渦流隨壓力的變化

由上一小節(jié)可以發(fā)現(xiàn),不同的壓力會使渦流產(chǎn)生變化．為了進一步觀察渦流隨壓力的變化,設(shè)出口壓力為0 Pa,采用不同的進口壓力(5 Pa,3 Pa,1 Pa,0 Pa,-1 Pa,-3 Pa,-5 Pa,-10 Pa)．此處為便于理解,壓力大小為正數(shù)時代表壓力梯度與電勢相同,反之則相反．圖7展示了在聚合物濃度為300 ppm,入口電勢大小E=20 V,t=0.8 s時刻下的流線圖．在10∶1∶10的縮放通道中,當壓力為0 Pa時(純EOF),縮放管的入口處產(chǎn)生渦流．

(a) 5 Pa (e) -1 Pa

(b) 3 Pa (f) -3 Pa

(c) 1 Pa (g) -5 Pa

(d) 0 Pa (h) -10 Pa圖7 Cp=300 ppm,入口處電勢E=20 V,t=0.8 s時刻下的流線,pout=0 Pa,采用不同的進口壓力Fig.7 For Cp=300 ppm,entrance potential E=20 V,the streamlines at t=0.8 s,pout=0 Pa,with different inlet pressures

如圖7(a)—7(d)所示,當壓力梯度和電勢梯度方向相同時,壓力越大,入口處的渦流會逐漸變小直至消失,這是因為正向壓力梯度的存在,流體在入口處的流速加快,使得由彈性不穩(wěn)定引起的渦流減?。蓤D7(c)可以發(fā)現(xiàn),當正向壓力為1 Pa時,渦流近似于消失狀態(tài),但入口處仍保持紊流．由圖7(b)可以看出當正向壓力足夠大時,縮放管入口處流線的紊亂情況減弱,渦流徹底消失．反之,當壓力梯度與電勢梯度相反時,如圖7(d)—7(h),反向壓力會阻礙流體在縮管入口的流入,流體在入口處聚集,從而使原本由EOF不穩(wěn)定性產(chǎn)生的渦流變大．反向壓力越大,入口處的渦流尺寸越大．同時,足夠大的反向壓力會使縮管中流速以壓力驅(qū)動為主導(dǎo),如圖7(h)所示．當反向壓力達到10 Pa時,壓力驅(qū)動在流動中占主導(dǎo)地位,渦流大小擴大到儲池最左端,縮放管中流體流向的轉(zhuǎn)變,渦流旋轉(zhuǎn)方向也倒轉(zhuǎn),大小和位置都發(fā)生巨大變化,此時渦流的邊界已經(jīng)超出儲池邊界．

為分析壓力對渦流大小的影響,采用如圖8所示的測量方法測量渦流大小[21],渦流大小XR以渦流最左端到渦流中心的距離為準．通過觀察發(fā)現(xiàn)在渦流不對稱時,上半部分渦流和下半部分渦流整體波動趨勢一致,故只測量下半部分渦流．

圖8 渦流大小示意圖Fig.8 Schematic diagram of the vortex size

圖9 E=20 V,Cp=300 ppm:渦流大小隨時間大小波動和不同出口壓力下渦流的平均大小Fig.9 For E=20 V,Cp=300 ppm:the vortex current size changing with time and the vortex current sizes under different outlet pressures

2.3 渦流隨聚合物濃度的變化

此部分分析不同聚合物濃度對縮放管入口渦流大小的影響．固定儲池進出口壓力,即給定一個反向壓力pout．分析從Newton流體到600 ppm聚合物濃度下的渦流大?。攭毫Ρ容^小時,如圖10中的1 Pa和3 Pa．此時與純電滲類似,渦流的大小隨著聚合物濃度的增大而增大．當壓力為3 Pa時,渦流大小在聚合物濃度達到500 ppm后趨于平緩．圖11為通道中點(0,0)處速度的離散系數(shù),離散系數(shù)的計算公式為cv=σ/μ,其中σ為標準差,μ為平均值．從圖中可以看出,隨著壓力的增大,同一聚合物濃度下離散系數(shù)均增大,速度的不穩(wěn)定性均增強,反向壓力存在可以增大黏彈性流體的不穩(wěn)定性．

圖10 E=20 V,渦流大小隨聚合物濃度的變化 圖11 通道中點速度的離散系數(shù)Fig.10 For E=20 V,the change of the vortex size with the polymer concentration Fig.11 The dispersion coefficient of the midpoint velocity of the channel

當壓力增大到5 Pa時,Newton流體在縮放管入口處已有渦流的產(chǎn)生,這與純電滲和較小反向壓力時不同,在純電滲或較小反向壓力驅(qū)動時,縮放管入口處沒有渦流產(chǎn)生．隨著聚合物濃度的增大,渦流在100 ppm達到一個最大值,然后隨著聚合物濃度的增大而減小,在聚合物濃度達到300 ppm后下降趨勢變緩．這是因為足夠大的反向壓力會使入口處的流速進一步降低,EOF驅(qū)動的流體在入口處無法進入,即使Newton流體也在入口處聚集,產(chǎn)生漩渦．當聚合物濃度達到100 ppm后,流體的彈性增大對渦流的影響大于速度減小對渦流的削弱影響,導(dǎo)致渦流進一步增大．但隨著聚合物濃度繼續(xù)增大,流體黏度增大,入口處的渦流速度降低,進而渦流減小并趨于穩(wěn)定．出口處壓力為5 Pa時,縮放管中點的速度變化規(guī)律也不再和低壓力時混合流或者純EOF下相同．從圖12可以看出,隨著聚合物濃度增大,流體黏度增加,速度受壓力影響減小,速度分布曲線的凹陷程度降低．聚合物濃度越小Helmholtz-Smoluchowsko速度越大[22],但反向壓力導(dǎo)致的速度減小程度也越大,進而出現(xiàn)通道中心速度在高聚合物濃度時比低聚合物濃度時的速度大．

圖12 pout=5 Pa,沿y軸的速度分布Fig.12 The velocity distributions along the y-axis,pout=5 Pa

3 結(jié) 論

本文采用PB模型與Oldroyd-B流體模型,建立了10∶1∶10的縮放微通道,對黏彈性流體的電滲壓力混合驅(qū)動流進行了數(shù)值研究,提出了壓力對黏彈性流體EOF不穩(wěn)定性的影響．在一定的電勢梯度下,當聚合物濃度為300 ppm時,分析了純EOF與具有反向壓力(pout=3 Pa)的電滲壓力混合驅(qū)動流．反向壓力的存在使得縮放管入口處由于黏彈性流體彈性不穩(wěn)定性產(chǎn)生的渦流增大,壓力每增大1 Pa渦流變大25 μm,壓力增大縮管進出口速度與中點的速度方差都在一定程度上增大,速度不穩(wěn)定性現(xiàn)象增強．

對于Newton流體與300 ppm的Oldroyd-B流體,疊加原理在Newton流體中適用,而在聚合物濃度為300 ppm的Oldroyd-B流體中由于黏彈性應(yīng)力存在,疊加速度比混合驅(qū)動的速度要大,產(chǎn)生了約4.3%的誤差．

在反向壓力較小時,渦流隨著聚合物濃度的增大而增大,并逐漸趨于穩(wěn)定．而在反向壓力較大時,例如pout=5 Pa,渦流的大小有一個先降后升的情況．這是因為當反向壓力足夠大時,在Newton流體下在入口處已經(jīng)產(chǎn)生了渦流,而后隨著聚合物濃度的增大,彈性不穩(wěn)定性增強渦流變大,繼續(xù)增大聚合物濃度導(dǎo)致黏度增大,由于流體黏度越大速度越低,導(dǎo)致渦流大小趨于穩(wěn)定．