Modulo Didático Eletroosmotico

Prévia do material em texto

Escoamento Eletro-osmótico em Micro-canais 
Área: Microfluidica e Micro-sistemas Eletromecânicos
Carolina Palma Naveira Cotta
Introdução
O presente módulo didático da disciplina Microfluidica e Micro-sistemas Eletromecânicos trata da solução de um problema exemplo,
usando a plataforma Mathematica, que lida com o escoamento eletro-osmótico em um microcanal com paredes de silicio e vidro, a
partir de uma solução eletrolitica que se deseja escoar. O objetivo é ilustrar o escoamento eletro-osmótico, como uma opção
interessante de bombeamento sem partes móveis em micro-sistemas com fluidos condutivos.
Escoamento Eletro-osmótico
Ÿ Problema
Considere um micro-canal de placas planas paralelas formado por uma parede de silicio e outra de vidro pyrex. O espaçamento
entre as paredes, h, é de 10 microns e a largura do canal, w, é de 200 microns, o que resulta em uma razão de aspecto
suficientemente grande para que possamos aproximar a geometria por um canal de placas paralelas. O canal tem comprimento
L de 2 cm. Deseja-se estabelecer um escoamento eletro-osmótico completamente desenvolvido no canal acima descrito, em
uma solução eletrolítica de KCl (cloreto de potássio) em água à temperatura de 25 C, com constante dielétrica Ε=80, concen-
tração iônica igual a n0=6.022 1020 m-3 (número médio de ions positivos e negativos por unidade de volume), viscosidade absoluta
Μ=0.9 10-3kgm.s, massa especifica Ρ=998 kg/m3. Os íons formados tem valencia z=1 (K+ e Cl -). O potencial zeta então assume
os valores de -100 mV na parede de silicio e de -59 mV na parede de vidro. Utilize a aproximação de Debye-Huckel, para
linearizar o problema, e encontre em forma analitica o campo de velocidades estabelecido para um campo elétrico externo
aplicado de 20 V/mm. Depois, compute e plote o campo de velocidades.
Figura 1. Fotografia do canal silicio-vidro fornecido pelo ESPCI (Prof. Patrick Tabeling), França.
Ÿ Solução
Primeiramente, preparamos uma lista com todos os dados numéricos do problema.
2 Modulo Didatico Eletroosmotico.nb
In[1]:= Dados2 = 8T ® 25 + 273.15, h ® 10 ´ 10^-6, w ® 200 ´ 10^-6,
L ® 2*10^-2, Ε ® 80., Μ ® 0.9 ´ 10^-3, n0 -> 6.022*10^20,
Ρ ® 998., z ® 1, Ζ1 ® -100*10^-3, Ζ2 ® -59*10^-3,
kb ® 1.3806*10^-23, q ® 1.60219*10^-19, Ε0 ® 8.854*10^-12,
Ex ® 20H10^-3L, tfinal ® 110, dpdx ® -10^5<
Out[1]= :T ® 298.15, h ®
1
100000
, w ®
1
5000
, L ®
1
50
, Ε ® 80., Μ ® 0.0009,
n0 ® 6.022´1020, Ρ ® 998., z ® 1, Ζ1 ® -
1
10
, Ζ2 ® -
59
1000
, kb ® 1.3806´10-23,
q ® 1.60219´10-19, Ε0 ® 8.854´10-12, Ex ® 20 000, tfinal ®
1
10
, dpdx ® -100 000>
O balanço de quantidade de movimento, para escoamento completamente desenvolvido, é então escrito em
função do potencial elétrico, com as condições de não-deslizamento nas duas paredes:
In[2]:= equ = Μ*D@u@yD, y, yD Š Ε*Ε0*Ex*D@Ψ@yD, y, yD
Out[2]= Μ u¢¢@yD Š Ex Ε Ε0 Ψ¢¢@yD
In[3]:= ccu1 = u@0D Š 0
ccu2 = u@hD Š 0
Out[3]= u@0D Š 0
Out[4]= u@hD Š 0
A equação que governa a distribuição espacial do potencial elétrico (sem aproximação de Debye-
Huckell) é dada como:
In[5]:= eqΨ = D@Ψ@yD, y, yD Š 2*n0*z*q *Sinh@z*q *Ψ@yDHkb*TLDHΕ*Ε0L
Out[5]= Ψ¢¢@yD Š
2 n0 q z SinhB q z Ψ@yD
kb T
F
Ε Ε0
Admitindo-se válida a aproximação de Debye-Huckell, a equação é linearizada na forma abaixo, com as condições
de contorno referentes aos potenciais zeta em cada parede:
In[6]:= eqΨnew = D@Ψ@yD, y, yD Š H2*n0*z*q HΕ*Ε0LL*z*q *Ψ@yDHkb*TL
Out[6]= Ψ¢¢@yD Š
2 n0 q2 z2 Ψ@yD
kb T Ε Ε0
In[7]:= ccΨ1 = Ψ@0D Š Ζ1
ccΨ2 = Ψ@hD Š Ζ2
Out[7]= Ψ@0D Š Ζ1
Out[8]= Ψ@hD Š Ζ2
A tentativa de determinação analitica do perfil de potencial elétrico (sem aproximação de Debye-Huck-
ell) através da 
função DSolve do Mathematica não é bem sucedida :
Modulo Didatico Eletroosmotico.nb 3
A tentativa de determinação analitica do perfil de potencial elétrico (sem aproximação de Debye-Huck-
ell) através da 
função DSolve do Mathematica não é bem sucedida :
In[9]:= DSolve@8eqΨ, ccΨ1, ccΨ2<, Ψ@yD, yD
Solve::ifun :
Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete
solution information.‡
DSolve::bvfail : For some branches of the general solution, unable to solve the conditions.‡
DSolve::bvfail : For some branches of the general solution, unable to solve the conditions.‡
Out[9]= 8<
Por outro lado, a determinação analitica do perfil de potencial elétrico com a aproximação de Debye-
Huckell é facilmente obtida:
In[10]:= Ψa@y_D = FullSimplify@
ExpToTrig@DSolve@8eqΨnew, ccΨ1, ccΨ2<, Ψ@yD, yD@@1, 1, 2DDD,
Assumptions ® 8n0 > 0, T > 0, kb > 0, Ε > 0, Ε0 > 0<D
Out[10]= CschB
2 h n0 q z
kb n0 T Ε Ε0
F Ζ1 SinhB
2 q Hh - yL z
kb T Ε Ε0
n0
F + Ζ2 SinhB
2 n0 q y z
kb n0 T Ε Ε0
F
A partir dessa solução, pode-se resolver diretamente o problema de escoamento, também utilizando a 
função DSolve para o campo de velocidades:
In[11]:= systemu =
Simplify@8equ, ccu1, ccu2< . D@Ψ@yD, y, yD ® D@Ψa@yD, y, yD,
Assumptions ® 8n0 > 0, T > 0, kb > 0, Ε > 0, Ε0 > 0<D
Out[11]= :2 Ex n0 q2 z2 CschB
2 h n0 q z
kb n0 T Ε Ε0
F Ζ1 SinhB
2 q Hh - yL z
kb T Ε Ε0
n0
F + Ζ2 SinhB
2 n0 q y z
kb n0 T Ε Ε0
F Š
kb T Μ u¢¢@yD, u@0D Š 0, u@hD Š 0>
In[12]:= ua@y_D = FullSimplify@DSolve@systemu, u@yD, yD@@1, 1, 2DD,
Assumptions ® 8n0 > 0, T > 0, kb > 0, Ε > 0, Ε0 > 0<D
Out[12]=
Ex Ε Ε0 -h Ζ1 + y HΖ1 - Ζ2L + h CschB 2 h n0 q z
kb n0 T Ε Ε0
F Ζ1 SinhB 2 q Hh-yL z
kb T Ε Ε0
n0
F + Ζ2 SinhB 2 n0 q y z
kb n0 T Ε Ε0
F
h Μ
Determina-se o comprimento de Debye em nanômetros:
4 Modulo Didatico Eletroosmotico.nb
In[13]:= Κ = Sqrt@2*n0*z^2*q^2HΕ*Ε0*kb*TLD
Out[13]= 2
n0 q2 z2
kb T Ε Ε0
In[14]:= Λ = 1Κ
Out[14]=
1
2
n0 q2 z2
kb T Ε Ε0
In[15]:= Λnm = Λ*10^9 . Dados2
Out[15]= 307.091
In[16]:= Z = hΛ . Dados2
Out[16]= 32.5637
Assim, ilustramos o perfil de potencial elétrico através do canal, em Volts (Região Neutra em cerca de 3 a 5 Ld):
In[17]:= Ψanum@y_D = Ψa@yD . Dados2  Simplify;
In[18]:= gΨa = Plot@Ψanum@y1000000D, 8y, 0, h*1000000 . Dados2<,
AxesLabel ® 8"y@ΜmD", "Ψ@VD"<, PlotRange ® 8-0.1, 0<D
Out[18]=
2 4 6 8 10
y@ΜmD
-0.10
-0.08
-0.06
-0.04
-0.02
Ψ@VD
Obtém-se também o perfil de velocidades e a posição do seu máximo na coordenada transversal do canal,
mostrado em [mm/s]:
Modulo Didatico Eletroosmotico.nb 5
In[19]:= uanum@y_D = ua@yD . Dados2  Simplify;
In[20]:= du@y_D = D@uanum@yD, yD;
In[21]:= eqmax = Hdu@yD  SimplifyL Š 0;
In[22]:= ymax = Solve@eqmax, yD@@1, 1, 2DD
Solve::ifun :
Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete
solution information.‡
Out[22]= 1.34346´10-6
In[23]:= umax = uanum@ymaxD
Out[23]= 0.00146752
In[24]:= umax*1000H*mms*L
Out[24]= 1.46752
Finalmente mostra-se o campo de velocidades para o caso em análise, onde claramente se observa as maiores
velocidades mais próximo à parede com maior potencial zeta.
In[25]:= gua = Plot@uanum@y1000000D*1000, 8y, 0, h*1000000 . Dados2<,
AxesLabel ® 8"y @ΜmD", "u @mmsD"<,
PlotRange ® 80, 1.2*umax*1000<D
Out[25]=
0 2 4 6 8 10
y @ΜmD
0.5
1.0
1.5
u @mmsD
Ÿ Análise Paramétrica
A seguir, com auxilio do comando Manipulate, podemos realizar uma análise paramétrica do comportamento do campo de veloci-
dades, variando-se os potenciais zeta nas duas paredes, o campo elétrico externo imposto, e a concentração iônica da solução.
6 Modulo Didatico Eletroosmotico.nb
In[26]:= Dados = 8T ® 25 + 273.15, h ® 10 ´ 10^-6, w ® 200 ´ 10^-6,
L ® 2*10^-2, Ε ® 80., Μ ® 0.9 ´ 10^-3, Ρ ® 998., z ® 1,
kb ® 1.3806*10^-23, q ® 1.60219*10^-19, Ε0 ® 8.854*10^-12<;
In[27]:= uamanip@y_D = ua@yD . 8Ζ1 ® x1, Ζ2 ® x2< . Dados;
In[28]:= Manipulate@Plot@8uamanip@y1000000D*1000 .
8x1 -> Ζ1*10^-3,x2 -> Ζ2*10^-3, Ex ® V *10^3, n0 ® n*10^20<,
uanum@y1000000D*1000<, 8y, 0, h*1000000 . Dados<,
AxesLabel ® 8"y @ΜmD", "u @mmsD"<, PlotRange ® 80, 5<,
PlotStyle ® 8RGBColor@1, 0, 0D, RGBColor@0, 0, 1D<D,
88Ζ1, -100, "Potencial Zeta 1 @mVD"<, 0, -250<,
88Ζ2, -59, "Potencial Zeta 2 @mVD"<, 0, -250<,
88V, 20, "Campo elétrico @VmmD"<, 5, 25<,
88n, 6.022, "Concentração ionica @10^20 ionsm3D"<, 0.1, 20<D
Out[28]=
Potencial Zeta 1 @mVD
Potencial Zeta 2 @mVD
Campo elétrico @VmmD
Concentração ionica @10^20 ionsm3D
0 2 4 6 8 10
y @ΜmD
1
2
3
4
5
u @mmsD
Ÿ Superposição de escoamentos eletro-osmótico e por diferença de pressão
Solução do escoamento só com o gradiente de pressão:
In[29]:= equpp = Μ*D@u@yD, y, yD Š dpdx
Out[29]= Μ u¢¢@yD Š dpdx
In[30]:= systemupp =
Simplify@8equpp, ccu1, ccu2< . D@Ψ@yD, y, yD ® D@Ψa@yD, y, yD,
Assumptions ® 8n0 > 0, T > 0, kb > 0, Ε > 0, Ε0 > 0<D
Modulo Didatico Eletroosmotico.nb 7
Out[30]= 8dpdx Š Μ u¢¢@yD, u@0D Š 0, u@hD Š 0<
In[31]:= uapp@y_D = FullSimplify@DSolve@systemupp, u@yD, yD@@1, 1, 2DD,
Assumptions ® 8n0 > 0, T > 0, kb > 0, Ε > 0, Ε0 > 0<D
Out[31]=
dpdx y H-h + yL
2 Μ
In[32]:= uappnum@y_D = uapp@yD . Dados2  Simplify;
Solução do escoamento devido ao campo elétrico e ao gradiente de pressão:
In[33]:= equp = Μ*D@u@yD, y, yD Š Ε*Ε0*Ex*D@Ψ@yD, y, yD + dpdx
Out[33]= Μ u¢¢@yD Š dpdx + Ex Ε Ε0 Ψ¢¢@yD
In[34]:= systemup =
Simplify@8equp, ccu1, ccu2< . D@Ψ@yD, y, yD ® D@Ψa@yD, y, yD,
Assumptions ® 8n0 > 0, T > 0, kb > 0, Ε > 0, Ε0 > 0<D
Out[34]= :dpdx kb T + 2 Ex n0 q2 z2 CschB
2 h n0 q z
kb n0 T Ε Ε0
F Ζ1 SinhB
2 q Hh - yL z
kb T Ε Ε0
n0
F + Ζ2 SinhB
2 n0 q y z
kb n0 T Ε Ε0
F Š
kb T Μ u¢¢@yD, u@0D Š 0, u@hD Š 0>
In[35]:= uap@y_D = FullSimplify@DSolve@systemup, u@yD, yD@@1, 1, 2DD,
Assumptions ® 8n0 > 0, T > 0, kb > 0, Ε > 0, Ε0 > 0<D
Out[35]=
1
2 h Μ
dpdx h y H-h + yL - 2 Ex Ε Ε0 Hh Ζ1 - y Ζ1 + y Ζ2L +
2 Ex h Ε Ε0 CschB
2 h n0 q z
kb n0 T Ε Ε0
F Ζ1 SinhB
2 q Hh - yL z
kb T Ε Ε0
n0
F + Ζ2 SinhB
2 n0 q y z
kb n0 T Ε Ε0
F
In[36]:= uapnum@y_D = uap@yD . Dados2  Simplify;
Perfil de velocidade com atuações opostas do campo elétrico e do gradiente de pressão:
In[37]:= uapbacknum@y_D = uap@yD . dpdx ® -dpdx . Dados2  Simplify;
Analise do Perfil de velocidade 
8 Modulo Didatico Eletroosmotico.nb
In[38]:= guap = Plot@8uanum@y1000000D*1000, uapnum@y1000000D*1000,
uappnum@y1000000D*1000, uapbacknum@y1000000D*1000<,
8y, 0, h*1000000 . Dados2<, AxesLabel ®
8"y @ΜmD", "u @mmsD"<, PlotRange ® 8-0.2, 2*umax*1000<D
Out[38]=
2 4 6 8 10
y @ΜmD
0.5
1.0
1.5
2.0
2.5
u @mmsD
Ÿ Verificação da Aproximação de Debye-Huckell
O balanço de quantidade de movimento, para escoamento completamente desenvolvido, (sem aproxi-
mação de Debye-Huckel)
Pseudo-Transiente : Como condição inicial foi dada a solução aproximada com Debye-Huckel
In[39]:= equ2 = Μ*D@ut@y, tD, y, yD Š
D@ut@y, tD, tD + Ex*2*n0*z*q *Sinh@z*q *Ψt@y, tDHkb*TLD
Out[39]= Μ utH2,0L@y, tD Š 2 Ex n0 q z SinhB
q z Ψt@y, tD
kb T
F + utH0,1L@y, tD
 Como condição inicial foi dada a solução aproximada com Debye-Huckel
In[40]:= ccui = ut@y, 0D Š uanum@yD
Out[40]= ut@y, 0D Š
0.00157404 - 64.5358 y - 0.00157404 CoshA3.25637´106 yE + 0.00157404 SinhA3.25637´106 yE
In[41]:= ccu1t = ut@0, tD Š uanum@0D
ccu2t = ut@h, tD Š uanum@hD
Out[41]= ut@0, tD Š 0.
Out[42]= ut@h, tD Š
0.00157404 - 64.5358 h - 0.00157404 CoshA3.25637´106 hE + 0.00157404 SinhA3.25637´106 hE
Potencial elétrico (sem aproximação de Debye - Huckell) 
Pseudo-Transiente : Como condição inicial foi dada a solução aproximada com Debye-Huckel
In[43]:= eqΨ2 = D@Ψt@y, tD, y, yD Š
D@Ψt@y, tD, tD + 2*n0*z*q *Sinh@z*q *Ψt@y, tDHkb*TLDHΕ*Ε0L
Modulo Didatico Eletroosmotico.nb 9
Out[43]= ΨtH2,0L@y, tD Š
2 n0 q z SinhB q z Ψt@y,tD
kb T
F
Ε Ε0
+ ΨtH0,1L@y, tD
In[44]:= ccΨi = Ψt@y, 0D Š Ψanum@yD
ccΨ1t = Ψt@0, tD Š Ψanum@0D
ccΨ2t = Ψt@h, tD Š Ψanum@hD
Out[44]= Ψt@y, 0D Š -0.1 CoshA3.25637´106 yE + 0.1 SinhA3.25637´106 yE
Out[45]= Ψt@0, tD Š -0.1
Out[46]= Ψt@h, tD Š -0.1 CoshA3.25637´106 hE + 0.1 SinhA3.25637´106 hE
Determinação númerica do problema acoplado não-linear
In[47]:= syscoupled =
8equ2, ccui, ccu1t, ccu2t, eqΨ2, ccΨi, ccΨ1t, ccΨ2t< . Dados2;
In[48]:= solcoupled = NDSolve@syscoupled, 8ut@y, tD, Ψt@y, tD<,
8y, 0, h . Dados2<, 8t, 0, tfinal . Dados2<,
MaxStepSize ® 8h1000 . Dados2, tfinal100 . Dados2<,
PrecisionGoal ® 5D@@1DD
NDSolve::mxsst :
Using maximum number of grid points 10000 allowed by the MaxPoints or MinStepSize options
for independent variable y.‡
NDSolve::eerr :
Warning: Scaled local spatial error estimate of 9.010419809730595`*^8 at t = 0.1` in the direction of
independent variable y is much greater than prescribed error tolerance. Grid
spacing with 10001 points may be too large to achieve the desired accuracy
or precision. A singularity may have formed or you may want to specify a
smaller grid spacing using the MaxStepSize or MinPoints method options.‡
Out[48]= 8ut@y, tD ® InterpolatingFunction@880., 0.00001<, 80., 0.1<<, <>D@y, tD,
Ψt@y, tD ® InterpolatingFunction@880., 0.00001<, 80., 0.1<<, <>D@y, tD<
Solução numérica da Velocidade eletrosmotica 
In[49]:= utnum@y_, t_D = solcoupled@@1, 2DD
Out[49]= InterpolatingFunction@880., 0.00001<, 80., 0.1<<, <>D@y, tD
Solução numérica do Potencial elétrico 
In[50]:= Ψtnum@y_, t_D = solcoupled@@2, 2DD
Out[50]= InterpolatingFunction@880., 0.00001<, 80., 0.1<<, <>D@y, tD
10 Modulo Didatico Eletroosmotico.nb
In[51]:= Plot3D@utnum@y, tD, 8y, 0, h . Dados2<, 8t, 0, tfinal . Dados2<D
Out[51]=
0
5.´10-6
0.00001
0.00
0.05
0.10
0.0005
0.0010
In[52]:= Plot3D@Ψtnum@y, tD, 8y, 0, h . Dados2<, 8t, 0, tfinal . Dados2<D
Out[52]=
0
5.´10-6
0.00001
0.00
0.05
0.10
-0.010
-0.005
0.000
In[53]:= Plot@8Ψanum@y1000000D, Ψtnum@y1 000000, tfinal . Dados2D<,
8y, 0, h*1000000 . Dados2<,
AxesLabel ® 8"y@ΜmD", "Ψ@VD"<, PlotRange ® 8-0.1, 0<D
Out[53]=
2 4 6 8 10
y@ΜmD
-0.10
-0.08
-0.06
-0.04
-0.02
Ψ@VD
Modulo Didatico Eletroosmotico.nb 11
In[54]:= Plot@8uanum@y1000000D*1000,
utnum@y1000000, tfinal . Dados2D*1000<,
8y, 0, h*1000000 . Dados2<,
AxesLabel ® 8"y @ΜmD", "u @mmsD"<,
PlotRange ® 80, 1.2*umax*1000<D
Out[54]=
0 2 4 6 8 10
y @ΜmD
0.5
1.0
1.5
u @mmsD
In[55]:= Plot@8Ψanum@y1000000D, Ψtnum@y1 000000, tfinal . Dados2D<,
8y, 0, h10*1000000 . Dados2<,
AxesLabel ® 8"y@ΜmD", "Ψ@VD"<, PlotRange ® 8-0.1, 0<D
Out[55]=
0.2 0.4 0.6 0.8 1.0
y@ΜmD
-0.10
-0.08
-0.06
-0.04
-0.02
Ψ@VD
In[56]:= Plot@8uanum@y1000000D*1000,
utnum@y1000000, tfinal . Dados2D*1000<,
8y, 0, h10*1000000 . Dados2<,
AxesLabel ® 8"y @ΜmD", "u @mmsD"<,
PlotRange ® 80, 1.2*umax*1000<D
12 Modulo Didatico Eletroosmotico.nb
Out[56]=
0.0 0.2 0.4 0.6 0.8 1.0
y @ΜmD
0.5
1.0
1.5
u @mmsD
Referências
1. Karniadakis, G., Beskok, A., e Aluru, N.(2005) Microflows and Nanoflows : Fundamentals and Simulation, Springer, NY.
2. Nguyen, N.T. e Werely, S.T. (2006) Fundamentals and Applications of Microfluidics, Artech House, 2nd ed..
3. Mikhailov, M.D., and R.M. Cotta, “Mixed Symbolic-Numerical Computation of Convective Heat Transfer with Slip Flow in Microchan-
nels”, Int. Comm. Heat & Mass Transfer, Vol. 32, Issues 3-4 , pp. 341-348, February 2005.
4. Cotta, R.M., M.D. Mikhailov, and S. Kakaç, “Steady and Periodic Forced Convection in Microchannels”, NATO ASI - Advanced
Study Institute on Micro-Scale Heat Transfer: Fundamentals and Applications in Biological and Microelectromechanical systems,
Çesme, Turkey, July 18-30, 2004; also, NATO Science Series II: Mathematics, Physics and Chemistry, Microscale Heat Transfer:
Fundamentalsand Applications, vol.193, S. Kakaç et al. (eds.), pp.49-74, 2005.
5. Soares, P.O., R.M. Cotta, and A.J. Silva Neto, “Convection in Microchannels with Electroosmotic Flow : Integral Transforms in
Pseudo - Transient Formulation”, 12 th Brazilian Congress of Thermal Sciences and Engineering, ENCIT 2008, Belo Horizonte,
Brazil, November 10 - 14, 2008.
Modulo Didatico Eletroosmotico.nb 13