Langmuir 2007, 23, 12421-12428
12421
Generalized Model for Time Periodic Electroosmotic Flows with Overlapping Electrical Double Layers Suman Chakraborty* and Amit Kumar Srivastava Department of Mechanical Engineering, Indian Institute of Technology, Kharagpur-721302, India ReceiVed July 13, 2007. In Final Form: August 27, 2007 In this article, an analytical model is devised for analyzing time periodic electroosmotic flows through nanochannels within the continuum regime, without presuming the validity of the Boltzmann distribution of ionic charges. The charge density distributions are obtained from the conservation considerations of the individual ionic species and other thermochemical constraints and are subsequently utilized to derive the potential distribution within the electrical double layer (EDL). This, coupled with the Navier-Stokes equation, yields a closed-form expression of the timedependent velocity field that is valid under overlapped EDL conditions. This expression is first validated in asymptotic limits of thin EDLs, for which closed form expressions have been benchmarked in the literature. Further analyses are carried out to bring out the influences of the frequency of the electrical field on the electroosmotic flow features in the presence of overlapped EDLs.
Introduction Miniaturization and integration of biomedical and biotechnological protocols into microfabricated lab-on-a-chip devices,1-5 commonly referred to as “micro total analysis” (µTAS) systems, are becoming increasingly important in the fields of genomics, proteomics, and clinical analysis. Such bioanalytical microdevices quite often rely on the electroosmotic transport mechanisms for sample mixing, manipulation, and separation processes, primarily motivated by the advantageous features6-11 of scaled-down analysis times, the capability of operating without moving components, reduced sample volumes, and minimized sample dispersion effects. The fundamental origin of electroosmotic transport lies in the fact that solid substrates in contact with aqueous solutions of polar solvents are likely to acquire a net surface charge12 through the action of several possible mechanisms such as ion adsorption, isomorphous substitution, and the ionization of surface groups. As an example, for glass- or silicon-fabricated microfluidic devices, the surface silanol groups might deprotonate, giving rise to a negatively charged surface and an associated diffuse charged layer of hydronium ions adhering to the substrate. The formation of this charged layer, also termed the electric double layer (EDL), results in a redistribution of charge density close to the substrate-liquid interface.12 The oppositely charged ions are attracted to the counter ions of the surface, and the co-ions are repelled. These attractive and repulsive forces, combined with the random thermal motion of ions (entropic interactions), * Author to whom correspondence should be addressed. E-mail:
[email protected]. (1) Chakraborty, S.; Mittal, R. J. Appl. Phys. 2007, 01, 104901. (2) Das, S.; Chakraborty, S. Electrophoresis, in press, 2007 (DOI 10.1002/ elps.200700265). (3) Chakraborty, S. J. Phys. D 2006, 39, 5356-5363. (4) Chakraborty, S.; Paul, D. J. Phys. D 2006, 39, 5364-5371. (5) Chakraborty, S. Lab Chip 2005, 5, 421-430. (6) Das, S.; Chakraborty, S. Anal. Chim. Acta 2006, 559, 15-24. (7) Das, S.; Das, T.; Chakraborty, S. Microfluid. Nanofluid. 2006, 2, 37-49. (8) Das, S.; Das, T.; Chakraborty, S. Sens. Actuators, B 2006, 114, 957-963. (9) Das, S.; Chakraborty, S. J. Appl. Phys. 2006, 100, 014098. (10) Das, S.; Chakraborty, S. AIChE J. 2007, 53, 1086-1099. (11) Das, S.; Subramanian, S.; Chakraborty, S. Colloids Surf., B 2007, 58, 203-217. (12) Hunter, R. J. Zeta Potential in Colloid Science: Principles and Applications; Academic: New York, 1981.
result in a diffuse arrangement of charged species within the EDL. An electric potential gradient, applied tangentially to the EDL, therefore causes the mobile portion of the EDL to migrate toward the cathode or anode, depending on the polarity of the EDL. This electromigration of ions constituting the EDL causes a viscous shearing of the adjacent bulk-liquid molecules, thereby resulting in an effective bulk-liquid motion that is commonly known as electroosmotic flow (EOF). The typical features of electroosmotic transport, as described above, are likely to be substantially modified for several reasons. For instance, under the action of periodic electrical fields, the electroosmotic flow becomes transient and frequency-dependent. Commonly termed “ac electroosmosis”, this phenomenon has been attracting growing attention from the research community over the past few years, leading to several pioneering experimental10-17 and theoretical18-24 investigations. Dutta and Beskok18 were among the early researchers that analytically investigated the time periodic EOF between two parallel plates, illustrating interesting similarities/dissimilarities with the Stokes second problem. Using the Green function approach, Kang et al.19 developed a theoretical model for ac electroosmosis in cylindrical microchannels, accommodating the possibilities of incorporating an arbitrarily large wall (zeta) potential in the mathematical analysis framework. Oscillating laminar EOFs in (13) Minor, M.; van der Linde, A. J.; van Leeuwen, H. P.; Lyklema, J. J. Colloid Interface Sci. 1997, 189, 370-375. (14) Ramos, A.; Morgan, H.; Green, N. G.; Castellanos, A. J. Phys. D: Appl. Phys. 1998, 31, 2338-2353. (15) Ramos, A.; Morgan, H.; Green, N. G.; Castellanos, A. J. Colloid Interface Sci. 1999, 217, 420-422. (16) Barraga´n, V. M.; Bauza´, C. R. J. Colloid Interface Sci. 2000, 230, 359366. (17) Oddy, M. H.; Santiago, J. G.; Mikkelsen, J. C. Anal. Chem. 2001, 73, 5822-5832. (18) Dutta, P.; Beskok, A. Anal. Chem. 2001, 73, 5097-5102. (19) Kang, Y. J.; Yang, C.; Huang, X. Y. Int. J. Eng. Sci. 2002, 40, 22032221. (20) Bhattacharyya, A.; Masliyah, J. H.; Yang, J. J. Colloid Interface Sci. 2003, 261, 12-20. (21) Yang, J.; Bhattacharyya, A.; Masliyah, J. H.; Kwok D. Y. J. Colloid Interface Sci. 2003, 261, 21-31. (22) Erickson, D.; Li, D. Langmuir 2003, 19, 5421-5430. (23) Reppert, P. M.; Morgan, F. D. J. Colloid Interface Sci. 2002, 254, 372383. (24) Kang, Y.; Yang, C.; Huang, X. J. Micromech. Microeng. 2004, 14, 12491257.
10.1021/la702109c CCC: $37.00 © 2007 American Chemical Society Published on Web 10/20/2007
12422 Langmuir, Vol. 23, No. 24, 2007
infinitely extended microchannels have been subsequently investigated by several researchers, including Bhattacharyya et al.,20 Yang et al.,21 Erickson and Li,22 Reppert and Morgan,23 Kang et al.,24 Macros et al.,25 Luo et al.,26 Xuan and Li,27 and Olesen et al.28 Most of the continuum-based mathematical models of ac electroosmosis, as mentioned above, are essentially based on the presumption that the celebrated Boltzmann distribution of ionic charge density remains a valid proposition for the theoretical analysis. However, in many critical situations typically encountered in nanochannel flows, the characteristic EDL thickness may be of comparable order to that of the characteristic system length scales. In these situations, strong EDL interactions, eventually leading to overlapped EDL conditions, become a realistic possibility. As a consequence of such an EDL overlapping phenomenon, the far-field boundary conditions of a predefined free-stream ionic concentration and zero electrical potential can no longer be applied to the central plane of symmetry of the fluidic channel. This in turn implies that the Boltzmann equation for ionic charge distribution might no longer be applicable.29-31 Under these circumstances, the electrical charge density field needs to be more fundamentally calculated by invoking the pertinent chemical equilibrium conditions, consistent with the ionic species conservation equations that dictate the number concentration of cations and anions in the central planes of the channel. Despite the well-understood limitations of the Boltzmann distribution of ionic charges under overlapped EDL conditions mentioned above, no comprehensive study has yet been reported in the literature on the analysis of ac electroosmotic flows through nanochannels with possible EDL overlaps. The aim of the present work, therefore, is to develop a comprehensive mathematical model for ac electroosmotic flows in parallel plate channels, incorporating the possibility of the channel hydraulic radius to be substantially smaller than the characteristic EDL thickness, despite being operative well within the continuum regime. For ready comparisons with benchmark solutions on ac electroosmotic flows in microchannels with thin EDLs, the mathematical model is derived for time periodic fully developed electroosmotic flows in straight channels, without precluding the possibilities of overlapped EDL conditions. Striking dissimilarities with the cases characterized with thin EDL conditions are also pinpointed, with the later condition being obtained as an asymptotic limiting condition of the generalized mathematical model postulated in this work. The vorticity transport in the diffuse EDL is also analyzed, and the contrasting features with the PoissonBoltzmann-based theories of ac electroosmosis are emphatically addressed.
Generalized Model for Time Periodic Electroosmotic Flows The present mathematical analysis considers an electroosmotic flow between two parallel plates of large width, separated by a distance of 2h, as shown in Figure 1. The confinement is considered to be narrow enough that the thicknesses of the EDLs formed in the vicinity of the two plates may turn out to be larger (25) Marcos, C.; Yang, K.T.; Ooi, T. N.; Wong, J. Micromech. Microeng. 2005, 15, 301-312. (26) Luo, W-J.; Pan, Y-J.; Yang, R-J. J. Micromech. Microeng. 2005, 15, 463-473. (27) Xuan, X.; Li, D. J. Colloid Interface Sci. 2005, 289, 291-303. (28) Olesen, L. H.; Bruus, H.; Ajdari, A. Phys. ReV. E 2006, 73, 056313. (29) Qu, W.; Li, D. J. Colloid Interface Sci. 2000, 224, 397-407. (30) Ren, C. L.; Li, D. Anal. Chim. Acta 2005, 531, 15-23. (31) Talapatra, S.; Chakraborty, S. Eur. J. Mech. B/Fluids, in press, 2007 (doi:10.1016/j.euromechflu.2007.06.005).
Chakraborty and SriVastaVa
Figure 1. Schematic diagram depicting the overlapped EDL field between two parallel flat plates in the presence of a time-varying axial electric field. The distance h is assumed to be larger than the characteristic EDL thickness, λ.
than the microchannel depth parameter, h. A time periodic axial electric field, Ex ) E0 sin(Ωt) is applied to the system, which generates the necessary driving force for electroosmotic transport of the incompressible polar fluid contained between the two plates. It is important to mention in this context that in the continuum fluid flow theory governed by the Navier-Stokes equations the fluid properties (such as density and temperature) are assumed not to vary appreciably over length and time scales comparable to the molecular mean free path and the molecular relaxation time, respectively. Although local near-wall density fluctuations may not necessarily mean the breakdown of the continuum theory altogether, special treatment of the interfacial conditions may be necessary to capture the underlying implications in the scaledup continuum limits. Such analysis may necessitate the consideration of certain noncontinuum interaction mechanisms, such as short-range forces due to the discrete structure of the liquid molecules as well as hydrophobic or hydrophilic interactions between the liquid molecules and the atoms of the solid surface (and the corresponding hydration forces). These interactions give rise to a molecular packing close to the fluid-solid interfaces and the resulting distribution of solvation or structural forces. Another ultra-short-range interaction potential may also come into play because of the straining of the chemical bonds between the liquid molecules and the solvation-induced layer structure close to the walls, rendering the continuum hypothesis to be inappropriate over local scales. The characteristic length scales over which these noncontinuum effects become operative is typically on the order of 1 nm. Thus, for channel hydraulic diameters to be less than 10 times the characteristic dimension of typical fluid molecules, continuum simulation techniques tend to become vulnerable. However, in the present study, typical nanochannel dimensions that are at least 1 order of magnitude larger than the above-mentioned critical limit are considered, so the classical continuum conservation equations remain applicable. One may begin the analysis by referring to the Navier-Stokes equation for incompressible fluids, subjected to electroosmotic body forces, given as
F
(∂V∂tB + (VB‚∇)VB) ) -∇P + µ∇ BV + F BE 2
e
(1)
where P is the pressure, F is the fluid density, µ is the dynamic viscosity, and B V ) (u, V) is a continuity-satisfying (divergence
Model for Time Periodic Electroosmotic Flows
Langmuir, Vol. 23, No. 24, 2007 12423
free) velocity field. The final term on the right-hand side of eq 1 represents the electroosmotic body force, where Fe is the electric charge density, and B E is the externally applied electric field. For any binary fluid consisting of two kinds of ions of equal and opposite charge z+ and z- , the charge density is given as Fe ) (n+ - n-)ze, where n+ and n- are the number density of the positively and negatively charged ions, respectively. In the classical EDL theory, the distributions of n+ and n- are given by the Boltzmann distribution,12 which is explicitly expressible in the form of exponential functions of the local EDL potential. However, the validity of this Boltzmann distribution is based on the conjecture that the EDLs formed in vicinity of the microchannel walls do not overlap and protrude into the centerline, so a far-stream boundary condition for the ionic charge densities can be imposed on the channel centerline. In the present model, the above assumptions are not necessarily presumed to be valid. Accordingly, a more general analysis is executed, starting from the fundamental equations of thermodynamic equilibrium, following the work of Qu and Li.29 To begin with, we consider the aqueous phase in the vicinity of the bottom plate of the channel and note that for thermodynamic equilibrium the electrochemical potential of the ions needs to be constant everywhere, which implies that
dµ ji )0 dy
(2)
where the subscript “i” indicates type-i ions and the electrochemical potential, µi, is defined as
µ j i ) µi + zieψ
(3)
Here, µi and zi are the chemical potential and valence of type-i ions, and e is the charge of a proton. Furthermore, from thermodynamic considerations, one may write
µi ) µ0i + kBT ln ni
(4)
where µ0i is a constant for type-i ions, kB is the Boltzmann constant, T is the absolute temperature of the solution, and ni is the number concentration of the type-i ions. Differentiating eq 4 with respect to y, one obtains
dµi 1 dni ) kBT dy ni dy
(5)
Combining eq 5 with eqs 2 and 3, one gets
The net charge density per unit volume, Fe, in that case, can be described as
Fe ) ze(n+ - n-)
The fundamental equation for the electrostatic potential distribution within the system can now be described by the Poisson equation (not necessarily by its simplified form, which, coupled with the Boltzmann distribution, is known as the PoissonBoltzmann equation), which states that
Fe d2ψ )2 dy
where is the permittivity of the medium. By substituting eqs 8a and 8b into eq 10, one gets
[
(
nh- exp
n+ )
(
(
zie (ψ - ψh) kBT
ze (ψ - ψh) kBT
ze (ψ - ψh) kBT
(11)
d 2ψ ze2 h (n+ + nh-)ψ ) dy2 kBT ze h ze h (n + nh+)ψh (12) (n - nh+) kBT -
[
]
With the boundary conditions of ψ ) ψ0 at y ) 0 (where ψ0 is approximated as ζ, the zeta potential, in most electrokinetic models, although there is a minor difference between these two in a strict sense) and ψ ) ψh at y ) h, one gets, on solving eq 12,
ψ ) ψh + A1 exp(xK1y) + A2 exp(-xK1 y) - K2 (13)
(6)
)
)
ze (ψ - ψh) exp kBT
n- ) nh- exp
)]
(
For large characteristic values of zeψ/kBT, one needs to integrate eq 11 numerically to obtain the potential distribution as a function of y. However, for small values of zeψ/kBT, the Debye-Hu¨ckel linearization principle12 can be invoked so that eq 11 simplifies to
)
K1 )
ze2 h (n + nh-) kBT +
(13a)
K2 )
kBT (nh- - nh+) ze (nh + nh )
(13b)
-
(7)
A1 )
Thus, for a simple symmetric (z:z) electrolyte, one may write
nh+
)
ze ze d 2ψ ) - nh+ exp (ψ - ψh) 2 k dy BT
Equation 6 can be solved by imposing the following boundary conditions at y ) h, ψ ) ψh, and ni ) nhi to yield
(
(10)
where
zie dni dψ )ni kBT
ni ) nhi exp -
(9)
(8a) (8b)
A2 )
ψ0 - ψ h 1 - exp(2 xK1h)
ψ0 - ψh exp(2xK1h) - 1
+
+
K2 1 + exp(xK1 h)
exp(2xK1h) +
K2 exp(2xK1h) 1 + exp(xK1h)
(13c)
(13d)
The concentrations of cations and anions can be obtained by substituting eq 13 into eqs 8a and 8b to yield
12424 Langmuir, Vol. 23, No. 24, 2007
[
n+(y) ) nh+ -
n-(y) ) nh-
[
Chakraborty and SriVastaVa
K3 exp(xK1(y - h)) 2 K3 exp(xK1(h - y)) + K3 + 1 (14a) 2
]
K3 exp(xK1(y - h)) + 2 K3 exp(xK1 (h - y)) - K3 + 1 (14b) 2
]
nondimensional velocity in the form U(χ, θ) ) F(χ) G(θ) ) Im[exp(iθ) F(χ)] so that eq 17 becomes
[
exp(iθ) d2F 1 × 2 2 2n κ dχ oR ze ze (ψ - ψ) - nh- exp (ψ - ψ) nh+ exp kBT h kBT h
i exp(iθ) F(χ) )
(
(
)
(
))]
(18)
Because exp(iθ) cannot be zero at all times, one can obtain from eq 18
where
K3 )
(
nh- - nh+ nh-
+
nh+
For the solution of the fluid flow equation (eq 1), consistent with the charge density distribution given by eq 9, we assume hydrodynamically fully developed conditions for which the convective term in the Navier-Stokes equation becomes identically equal to zero. Furthermore, we consider a sinusoidal time periodic pure electroosmotic flow in the absence of any pressure gradients so that eq 1 simplifies to the following form
∂2u ∂u ) µ 2 + FeE0 sin(Ωt) ∂t ∂y
F
(15)
where is B E0 the magnitude and Ω is the frequency of the unsteady external electric field, B E, and Fe is given by eq 9. Substituting eq 10 in eq 15, one may obtain
[
(
nhze exp (ψ - ψ) E0 sin(Ωt) (16) 2n0 kBT h where n0 is the nominal number density of positive or negative ions in the buffer. Using the definition of the Debye-Hu¨ckel parameter, given by ω ) 1/λ ) x2noe2z2/kBT (with λ being the characteristic EDL thickness, also known as the Debye length), and R ) ezζ/kBT, eq 16 may be cast in the following nondimensional form
[
(
(
)
(
(
(17)
where θ ) Ωt is the nondimensional time, χ ) y/λ is the nondimensional distance, κ ) (Ωλ2/ν)1/2 is a normalized frequency, and U ) u/uHS is the nondimensional velocity, with uHS ) -ζE0/µ being the Helmholtz-Smoluchowski velocity (the characteristic electrosmotic flow velocity for nonoverlapped EDL fields). The κ parameter can be interpreted as the ratio of the Debye length to a diffusion length scale (lD), based on the kinematic viscosity and the excitation frequency where lD ≈ (V/Ω)1/2. To obtain a time periodic solution of eq 17, one may employ the method of separation of variables to cast the
(19)
CF ) A1 exp(xiκχ) + B1 exp(-xiκχ)
(20)
where A1 and B1 are the coefficients of the homogeneous solution, which depend on the κ parameter. The particular integral of eq 19 can be obtained as
PI )
1 2κRxi
[
(
exp(-xiκχ)∫0 exp(xiκχ) χ
)
1
(
2κRxi
(
[
exp(xiκχ)
nh+ exp × 2no
)
(
nh
))
]
dχ -
∫0χ exp(-xiκχ) -2n+o exp ×
(
nh-
))
ze ze (ψ - ψ) + (ψ - ψ) exp kBT h 2no kBT h
]
dχ (21)
The complete solution for the velocity field accordingly becomes
U(χ, θ) ) Im[exp(iθ)A1 exp(xi κχ) + exp(iθ)B1 exp × (-xiκχ)] + Im
))]
))
The complimentary function for the solution of the above differential equation can be written as
(
) )]
ze ∂U 1 ∂2U sin θ h ) (ψ - ψ) n exp ∂θ κ2 ∂χ2 2noR + kBT h ze nh- exp (ψ - ψ) kBT h
)
nhze ze exp (ψh - ψ) + (ψ - ψ) kBT 2no kBT h
nh+ ze ∂ 2u ∂u (ψh - ψ) exp F ) µ 2 + 2noez ∂t 2n k ∂y o BT
(
(
ze d 2F 1 (ψh - ψ) - iκ2F(χ) ) nh+ exp 2 2n R k dχ o BT ze nh- exp (ψ - ψ) kBT h
[ [ exp(iθ) 2κRxi
( ( ) [ [
exp(-xiκχ)
(
∫0χ exp(xiκχ) ×
nh+ nhze ze exp exp (ψh - ψ) (ψ - ψ) 2no kBT 2no kBT h Im
(
-
exp(iθ) 2κRxi
exp(xi κχ)
)
∫0
nh+
χ
(
]]
))
dχ
-
nhexp(-xiκχ) exp × 2no
(
ze ze exp (ψ - ψ) (ψ - ψ) kBT h 2no kBT h
))
]]
dχ
(22)
Coefficients of the homogeneous solution are found by employing the following boundary conditions: (i) the no-slip wall boundary condition (i.e., F(χ ) 0) ) 0) and (ii) the centerline symmetry condition (i.e., F/(χ ) χ0) ) 0, where χ0() h/λ) is the value of χ at the central plane of symmetry). These yield
Model for Time Periodic Electroosmotic Flows
[
∫0χ
exp(-xiκχ)
A1 )
0
Langmuir, Vol. 23, No. 24, 2007 12425
( ( [ ∫
exp(xiκχ)
)
(
nh+ nhze ze exp - (ψh - ψ) exp (ψ - ψ) 2no kBT 2no kBT h 4xiκR cosh(xiκχ)
exp(xiκχ)
χ0
0
(
))
(
dχ
)
]
+
(
nh+ nhze ze (ψ - ψ) (ψ - ψ) exp(-xiκχ) exp exp 2no kBT h 2no kBT h 4 xiκR cosh(xiκχ) B1 ) -A1
It is important to mention here that because of predominant hydrophobic interactions the boundary conditions at the walls may exhibit “apparent slip” behavior instead of following the classical “no-slip” conjecture. Such deviations, not being the focal point of concern in the present study, are not considered here. For detailed modeling of the wall boundary conditions in the case of predominant hydrophobic interactions, one may refer to the work of Chakraborty.32,33
Asymptotic Limits With the limiting conditions ψh f 0 and nh- ) nh+ ) n0, as h/λ f ∞, the solution given by eq 22 assumes the following simplified form
U(χ, θ) ) Im[exp(iθ)A2 exp(xiκχ) + exp(iθ)B2 exp × exp(iθ) χ exp(-xiκχ) 0 exp(xiκχ) × (-xiκχ)] + Im x 2κR i exp(iθ) χ ze sinh exp(xiκχ) 0 exp × ψ dχ - Im kBT x 2κR i ze (-xiκχ) sinh ψ dχ (24) kBT
[
[
∫
( ( )) ]] [
∫ [ ( ( )) ]]
where
[
exp(-xiκχ)
A2 )
[
∫0
χ0
( ( )) ]
ze ψ exp(xiκχ) sinh kBT
exp(xiκχ)
∫0χ
0
+
( ( )) ]
exp(-xiκχ) sinh
ze ψ kBT
4xiκR cosh(xiκχ)
(23a) (23b)
overlapped EDL conditions, are not independent constants but rather are dependent on certain additional constraints. These compatibility conditions can be exploited in the manner outlined in Qu and Li29 to determine the values of the above-mentioned input parameters. It is interesting that to obtain the solution of the velocity field it appears that the numerical value of the ψ0 parameter becomes inconsequential because it gets nullified in the expression for (ψh - ψ) that needs to be evaluated to obtain the electroosmotic body forces. However, the values of nh- and nh+ are themselves implicitly dependent on ψ0, for which the determination of all of these parameters becomes necessary. To illustrate the underlying principles, one may consider that the solid surface is composed of inorganic oxides (such as silica, zirconia, titania, etc.). The solid surface can become either positively or negatively charged, depending on whether H3O+ or OH- ions (formed as a consequence of the dissociation of the aqueous medium carrying the electrolyte) can be adsorbed to the solid surface, consistent with the pH of the solution. The resultant surface charge density is given as34
σ0 )
eNsδ sinh(ψN* - R) 1 + δ cosh(ψN* - R)
ψN* ) 2.303(pHz - pH)
dχ
(26)
(27)
(25)
and
B2 ) -A2 The above-mentioned special case, in essence, mimics the classical theory, which presumes that the characteristic EDL thickness (λ) is taken to be much less than h, so the location of the central plane can be treated mathematically as a far-field one. In that case, because the charge density gradients are confined only within the EDLs adjacent to the solid surfaces, the effects of the near-wall potential distribution cannot effectively penetrate into the channel centerline.
Results and Discussions To implement the mathematical model, it is important to recognize here that parameters nh-, nh+, and ψ0 (or equivalently, R), necessary for the solution of the potential field under (32) Chakraborty, S. Appl. Phys. Lett. 2007, 90, 034108. (33) Chakraborty, S. Phys. ReV. Lett. 2007, 99, 094504 (1-4).
]
dχ
where Ns is the site density on the oxide surface, δ ) 2 × 10-∆pK/2 (with ∆pK being the dissociation constant difference; typically ∆pK ∼10) and ψN* is the nondimensional Nernst potential, given by
dχ
4xiκR cosh(xiκχ)
))
In eq 27, pHz is the pH value when the solid surface reaches the point of zero charge. (For most oxide surfaces, this value is between 5 and 6.) For a given value of pH, the ionic concentrations [H3O+] and [OH-], consistent with the 2H2O S H3O+ + OHdissociation reaction, can be described as [H3O+ ] ) 10-pH and [OH-] ) 10pH - 14. Accordingly, the number density of H3O+ and OH- ions in the central plane of symmetry can be obtained h 3 as nHh 3O+ ) N × 103 × [H3O+] and nOH - ) N × 10 × [OH ], where N is Avogadro’s number and the ionic concentrations are expressed in M (mol/L) units. The average values of these ionic concentrations may also be assumed to obey the reaction equilibrium constraints, so one may write
njH3O+ × njOH- ) N2 × 106 × Kw
(28)
where nji ) 1/h ∫h0 ni dy. In eq 28, Kw is the dissociation constant, taken as Kw ) 10-14 at 25 °C,24 provided that the ionic concentrations are taken in units of M (mol/L). Furthermore, (34) Healy, T. W.; White, L. R. AdV. Colloid Interface Sci. 1978, 9, 303-345.
12426 Langmuir, Vol. 23, No. 24, 2007
Chakraborty and SriVastaVa
Table 1. Values of Coefficient A1, as Function of K, for r ) 0.4278 at Different Values of χ0 κ
χ0 ) 100
χ0 ) 2
χ0 ) 0.5
0.1 0.01 0.001
-3.0843 + 3.5173i -35.0691 + 35.6862i -354.7529 + 354.8226i
-3.4592 + 3.5309i -34.9529 + 34.9602i -349.5655 + 349.5662i
-3.4117 + 3.4172i -34.1444 + 34.1449i -341.4471 + 341.4472i
from the mass conservation of individual ionic types, one can write29
njH3O+h ) nHh 3O+ + Nd h njOH- h ) nOH - -
(29a)
σ0 + Nd e
(29b)
where Nd is the number of H2O molecules dissociated into H3O+ and OH- ions. It is also important to mention here that to achieve overall ionic charge balance one also needs to include the dissociated ions of the buffer salt (we take KCl as an example here) and write
∫0h n+ dy ) ∫0h nK
+
∫0h n- dy ) ∫0h nCl
dy +
∫0h nH O
dy +
∫0h nOH
-
3
+
dy
(30a)
-
dy
(30b)
Noting that ∫h0 nK+ dy ) ∫h0 nCl- dy ) n0h ∫h0nH3O+ dy ) njH3O+h and ∫h0 nOH- dy ) njOH-h, the above equations can be rewritten as
∫0h n+ dy ) n0h + njH O h
(31a)
∫0h n- dy ) n0h + njOH h
(31b)
3
+
-
where n+ and n- are given by eqs 8a and 8b as functions of nhand nh+. Additionally, one may use the compatibility condition that
h ψh f 0 and nh- ) nh+ ) n0 as f ∞ λ
(32)
which represents the physical situation with nonoverlapped EDL conditions (where the infinity condition is practically reached for h/λ ≈ 10). With the inputs of the n0, pH, pHz, ∆pK, and Ns parameters, eqs 26, 28, 29a, 29b, 31a, 31b, and 32 represent seven independent equations with seven unknowns σ0, R, nh-, nh+, njH3O+, njOH-, and Nd. These equations can be simultaneously solved to obtain the relevant model parameters consistent with the present formalism. In the present study, the input parameters for obtaining nh-, nh+, and R are taken as follows:24 n0 ) 10-5 M, Ns ) 5 × 1014 sites/cm2, pHz 5.8, pH 7, and ∆pK ) 10. With the corresponding values of nh- and nh+, the velocity profile given by eq 22 is uniquely determined. All integrals necessary for this evaluation are obtained by employing the trapeziodal rule of numerical integration. The values of the A1 parameter, corresponding to a chosen set of problem data, are shown in Table 1 for illustration. Results from the present simulation are first validated for the nonoverlapped EDL case, for which comprehensive semianalytical solutions are presented in the work of Dutta and Beskok.18 The time periodic velocity profiles across the channel section, corresponding to this case, are presented with the choices of parameters mentioned earlier for different values of the nondimensional frequency, κ, by taking χ0 (≡ h/λ) ) 100, as depicted in Figure 2, with θ ) π/2. Fundamentally, the fluid flow originates from the stress produced by the interaction between the externally
applied electric field and the charge density distribution due to the induced potentials, which is balanced by viscous stresses. For very low values of κ (κ ≈ 0.001), the diffusion length scale (ld ≈ λ/κ) is 1 order of magnitude larger than the order of magnitude of the channel height (∼100λ), so a plug-shaped quasisteady velocity profile is obtained in Figure 2, corresponding to that case. For κ ≈ 0.01, the channel hydraulic radius and the diffusional length scales become comparable in order, characterizing the onset of deviation from a plug-shaped velocity distribution. For nondimensional frequencies of even higher magnitudes (κ ≈ 0.1), appreciable reductions in velocity are observed to occur outside the EDL. Such physical behavior is in accordance with the fact that at lower frequencies almost the entire potential drop takes place within the EDL whereas at higher frequencies the potential drop across the EDL goes to zero because of the polarization of the substrate/solution interface. The above can be explained from the consideration that at high frequencies the capacitance of the EDL tends to a constant limit, which is equivalent to that of a parallel plate capacitor with the permittivity of the bulk electrolyte. At high frequencies, the resistance also tends to be a constant, in inverse proportionality to the conductivity of the solution. These limiting values of capacitance and resistance primarily correspond to the impedance of the bulk electrolyte because the EDL impedance is negligible at high frequencies. The surface charge density, being proportional to the EDL potential drop, is very small at high frequencies. At lower frequencies, however, the net resistance and capacitance include the EDL impedance effects as well. Because the potential across the EDL is very small at low frequencies, the tangential field at the EDL is also accordingly small. Furthermore, it is interesting to observe from Figure 2 that in all cases the results predicted by using the overlapped and nonoverlapped EDL models effectively coincide. This can be attributed to the fact that the ionic concentration distribution tends to be the same as that predicted by the nonoverlapped EDL theory for cases in which the opposite plates are kept a sufficient distance apart so that the respective EDLs do not interfere with each other, as applicable to the situations depicted in Figure 2. To illustrate the implications of EDL overlap on the timedependent electroosmotic flow solutions, Figure 3a,b is plotted, representing the velocity profiles corresponding to χ0 ) 0.5 for κ ) 0.001 and 0.1, respectively, for different values of θ. In general, it is observed that the magnitudes of the velocities predicted with an overlapped EDL model are found to be higher than that predicted by classical theory. This is fundamentally attributed to the fact that the overlapped EDL theory predicts a stronger local EDL potential gradient than the classical theory. Under overlapped EDL conditions, the ionic concentration distribution differs considerably from that obtained with a thin EDL approximation. The key factor here is the variation in the local charge density, which is directly proportional to the difference between the anion and cation concentrations. With EDL overlap in the vicinity of negatively charged substrates, the number concentration of anions increases whereas the number concentration of cations decreases with increasing distance from the solid surface. Moreover, the local concentration difference between the anions and the cations turns out to be more for the overlapped EDL case (because of stronger concentration gradients prevailing near the solid boundaries). Thus, for channels with
Model for Time Periodic Electroosmotic Flows
Langmuir, Vol. 23, No. 24, 2007 12427
Figure 2. Velocity profiles for time periodic electroosmotic flows at various nondimensional frequencies, for θ ) π/2, with χ0 ) 100.
actual EDL overlaps (as for this case), appreciable velocity gradients persist throughout the entire channel section, which results in stronger electroosmotic effects and hence higher values of the peak velocity as compared to that predicted by the nonoverlapped EDL theory. Regarding the implications of the frequency of the applied electrical field on the velocity variations, it is interesting to observe from Figure 3a,b that the extent of the frequency dependence is not found to be as conspicuous as that observed for channels with much thinner EDLs. In general, it is intuitive that for a thin EDL, the fluid elements located outside the same can remain rather insensitive to the charge density gradients imposed on the system near the electrode/solution interface. This insensitiveness is further aggravated by the polarization of the electrode/fluid interface at relatively higher frequencies. In fact, at higher frequencies, the time scale of establishment of the charge distribution, or the so-called relaxation time, turns out to be so large that changes in the charge density distribution cannot effectively propagate until reaching the channel centerline. However, for overlapped EDLs, the message of the presence of the EDL potential gradient propagates effectively into the channel centerline, so that the above situation does not feature with prominence. The higher velocities predicted with overlapped EDL theory, under these conditions, can be attributed to the effective penetration of the EDL effects into the channel centerline, despite the flow suppression tendencies induced by the high-frequency effects. At low frequencies, however, the entire potential drop virtually takes place within the EDL. In such cases, the strengthened velocity field predicted by the overlapped EDL theory is primarily attributable to the higher charge density gradients obtained under overlapped EDL conditions. The frequency sensitivity of the velocity variations obtained in Figure 3a,b is further suppressed by the fact that both plots are normalized with respect to the distance χ/χ0, which always scales between 0 and 1, irrespective of the actual channel dimension. This enables one to obtain a kind of universal velocity characteristic that is only weakly sensitive to the imposed frequencies for strongly overlapped EDL conditions. The key issue that determines this factor is the effective penetration of the EDL effects into the bulk, rather than being dictated by the frequency alone to determine the extent of this penetration. However, it is interesting to observe that a strong θ dependence of velocity profiles is still exhibited despite the overlapped EDL conditions and the relatively weaker frequency dependence exhibited as a consequence. The most significant deviations in the velocity distributions, as predicted by the modified and the classical theories, can be observed for θ ) π/2, for which the
Figure 3. Velocity profiles for time periodic electroosmotic flows at various nondimensional times, with χ0 ) 0.5, for (a) κ ) 0.001 and (b) κ ) 0.1.
applied electric field operates with its greatest possible strength so that higher values of charge densities predicted with overlapped EDL condition can amplify the effects of the applied field to the greatest possible extent. Figure 4a,b depicts the normalized velocity gradients as a function of the nondimensional distances from either of the two plates, corresponding to χ0 ) 2 and 0.5, respectively. With the transverse components of velocity being identically equal to zero for fully developed flows, these plots are also reminiscent of the vorticity distribution along the channel cross section. For either of the two cases, this normalized vorticity distribution is only marginally sensitive to frequency variations, which is in sharp contrast to the situation observed for thinner EDLs. It can be observed by comparing Figure 4a,b that the higher value of χ0 is characterized with higher magnitudes of wall vorticity, expressed in terms of normalized coordinates. Because the vorticity is zero at the channel centerline, there is a gradual reduction of vorticity with increasing distance from the solid boundaries. For very narrow channels with strong EDL interactions (such as the case with χ0 ) 0.5), the relatively short distance over which the vorticity values transit from the wall condition to the centerline condition renders the corresponding profiles to be approximately linear. These profiles, however, are predominantly transient in nature because they exhibit strong dependences on the temporal parameter, θ. In sharp contrast to the thin EDL conditions, no maxima or minima are observed in the vorticity
12428 Langmuir, Vol. 23, No. 24, 2007
Chakraborty and SriVastaVa
in the flow field, thereby disallowing them to grow any further. It is also important to mention here that because with strongly interacting EDLs the characteristic velocity scales are somewhat smaller than the Helmholtz-Smoluchowski slip velocity (UHS), analogies with a time periodic flow with a prescribed wall-slip condition (the well-known Stokes second problem) cannot be effectively drawn. In fact, with EDLs of thicknesses comparable to the nanochannel hydraulic radii, substantial velocity gradients (and hence vorticities) do exist across the entire channel section, rendering the Helmholtz-Smoluchowski velocity to be a grossly inappropriate electroosmotic slip condition for such flows. Rather, the velocity and vorticity distributions qualitatively resemble pressure-driven time-oscillating flows under these conditions, leading to substantial effects of dispersion depending on the operating temporal regime. Nevertheless, because pressure-driven transport necessitates huge pumping power to transmit liquid samples through such narrow slits, electroosmosis remains the preferred choice of flow actuation under strongly overlapped EDL conditions.
Conclusions
Figure 4. Normalized cross velocity gradient profiles for time periodic electroosmotic flows at various nondimensional times with (a) χ0 ) 2 and (b) χ0 ) 0.5. The curves for κ ) 0.001 and 0.1 almost coincide for both cases, with a closer coincidence for case b.
distributions with strong EDL interactions13 because of the lack of any inflection points in the corresponding velocity profiles. Therefore, any flow instabilities as a consequence of vorticity inversion are virtually ruled out for strongly interacting EDLs, as compared to the classical situations marked with ultrathin EDLs. Physically, with strongly overlapped EDL conditions, the diffusive effects are strong enough to dissipate any perturbations
Analytical solutions have been obtained for time periodic electroosmotic flows under overlapped EDL conditions, which can be of potential importance in transmitting and manipulating biological and biochemical samples through nanoscale confinements. However, in sharp contrast to the corresponding phenomenon with thin EDL conditions, the following interesting behavior has been carefully noted for this regime: (i) Because of strong ionic concentration gradients prevailing over the entire channel section, no analogies can be drawn with time periodic fluid flows with oscillating wall velocity conditions (Stokes second problem). (ii) The transience in velocity variation is further amplified with overlapped EDL considerations because of the above effects. (iii) The frequency dependence of normalized flow velocity happens to be rather weak for overlapped EDL conditions, as compared to that for nonoverlapped EDL conditions. (iv) Although nonzero values of vorticity do exist over almost the entire channel section (except for the channel centerline) for channels with EDL thicknesses of comparable order to the characteristic system length scale, no extremum conditions of vorticity can be obtained in between under overlapped EDL conditions, thereby precluding any possibility of the onset of a resulting flow instability. Physically, vorticities diffuse deeper into the channel for thicker EDLs, suppressing the further growth of any perturbations. LA702109C
