Anal. Chem. 2000, 72, 4767-4777
Accelerated Articles
Electroosmotic Fluid Motion and Late-Time Solute Transport for Large Zeta Potentials Stewart K. Griffiths* and Robert H. Nilson
Sandia National Laboratories, P.O. Box 969, Livermore, California 94551-0969
Analytical and numerical methods are employed to determine the electric potential, fluid velocity, and late-time solute distribution for electroosmotic flow in a tube and channel at ζ potentials that are not necessarily small. The electric potential and fluid velocity are in general obtained by numerical means. In addition, new analytical solutions are presented for the velocity in a tube and channel in the extremes of large and small Debye layer thickness. The electroosmotic fluid velocity is used to analyze latetime transport of a neutral nonreacting solute. Zero- and first-order solutions describing axial variation of the solute concentration are determined analytically. The resulting expressions contain eigenvalues representing the dispersion and skewness of the axial concentration profiles. These eigenvalues and the functions describing transverse variation of the concentration field are determined numerically using a shooting technique. Results are presented for both tube and channel geometries over a wide range of the normalized Debye layer thickness and ζ potential. Simple analytical approximations to the eigenvalues are also provided for the limiting cases of large and small values of the Debye layer thickness. Integrated microchannel devices offer the potential for increased speed and reduced cost in the analysis and synthesis of chemical and biological species.1-5 Such devices typically employ transverse channel dimensions from a few micrometers to ∼1 mm (1) Manz, A.; Fettinger, J. C.; Verpoorte, E.; Ludi, H.; Widmer, H. M.; Harrison, D. J. Trends Anal. Chem. 1991, 10, 144-149. (2) Manz, A.; Harrison, D. J.; Verpoorte, E. M. J.; Fettinger, J. C.; Paulus, A.; Ludi, H.; Widmer, H. M. J. Chromatogr. 1992, 593, 253-258. (3) Harrison, D. J.; Manz, A.; Fan, Z. H.; Ludi, H.; Widmer, H. M. Anal. Chem. 1992, 64, 1926-1932. (4) Effenhauser, C.; Manz, A.; Widmer, H. M. Anal. Chem. 1993, 65, 26372642. (5) Jacobson, S. C.; Hergenroder, R.; Koutny, L. B.; Ramsey, J. M. Anal. Chem. 1994, 66, 1114-1118. 10.1021/ac000539f CCC: $19.00 Published on Web 09/12/2000
© 2000 American Chemical Society
and have longitudinal dimensions of only a few centimeters. As a result, they may also enable compact portable systems for analysis and synthesis having process capabilities comparable to those now available only in stationary equipment. Electroosmotic flows6,7 are widely used in microchannel systems for both separation processes and routine sample handling.8-15 One benefit of such flows is that species samples may be transported over long ranges with little dispersion16 due to nonuniform fluid speeds provided that the Debye layer thickness is small compared to the tube or channel width. Under this restriction, the velocity profile across a channel is very flat, and all variations in the fluid speed occur close to the channel walls. However, as the Debye layer grows or the channel width is reduced, the velocity profile approaches a parabola and the accompanying dispersion increases to that of a pressure-driven flow. In contrast, increasing ζ potentials yield a more uniform velocity17 and so tend to counteract this influence of increased Debye layer thickness. Dispersion is thus important when either channel geometry or the operating conditions are being optimized. Most previous studies of transport and dispersion in electroosmotic flow have employed the Debye-Hu¨ckel approximation (6) Probstein, R. F. Physicochemical Hydrodynamics; John Wiley & Sons: New York, 1995. (7) Rice, C. L.; Whitehead, R. J. Phys. Chem. 1965, 69, 4017-4024. (8) Chiem, N.; Harrison, D. J. Anal. Chem. 1997, 69, 373-378. (9) Effenhauser, C. S.; Bruin, G. J. M.; Paulus, A.; Ehrat, M. Anal. Chem. 1997, 69, 3451-3457. (10) Salimi-Moosavi, H.; Tang, T.; Harrison, D. J. J. Am. Chem Soc. 1997, 119, 8716-8717. (11) Dasgupta, P. K.; Liu, S. Anal. Chem. 1994, 66, 3060-3065. (12) Salimi-Moosavi, H.; Tang, T.; Harrison, D. J. Electrophoresis 2000, 21, 107115. (13) Jacobson, S. C.; Hergenroder, R.; Koutny, L. B.; Ramsey, J. M. Anal. Chem. 1994, 66, 2369-2373. (14) Kutter, J. P.; Jacobson, S. J.; Ramsey, J. M. Anal. Chem. 1997, 69, 51655171. (15) Kutter, J. P.; Jacobson, S. J.; Matsubara, N.; Ramsey, J. M. Anal. Chem. 1998, 70, 3291-3297. (16) Taylor, G. I. Proc. R. Soc. London 1953, A 219, 186-203. (17) Gross, R. J.; Osterle, J. F. J. Chem. Phys. 1968, 49, 228-234.
Analytical Chemistry, Vol. 72, No. 20, October 15, 2000 4767
that the ζ potential is small.18 This limitation permits linearization of the Poisson-Boltzmann equation governing the electric potential, enabling closed-form solutions describing the fluid velocity7 and, thereby, admits very general closed-form or series solutions for the coefficient of dispersion. Such studies have addressed combined pressure-driven and electrokinetic motion in a tube,19-21 adsorption of reactive species on tube walls,20,21 nonuniform ζ potentials,22,23 species distribution along and across a tube or channel,24 and the influence of a rectangular channel geometry.25-27 While analytical solutions are valuable in understanding the phenomena of dispersion, some practical applications may not satisfy the Debye-Hu¨ckel condition. Organic solvents in contact with glass may, for example, exhibit ζ potentials in excess of 200 mV.28 Only a few previous studies have addressed transport and dispersion when ζ potentials are large and the Debye-Hu¨ckel approximation is not valid. Gas et al. analyzed flow and electrokinetic transport in a tube at large ζ potentials by assuming that the effective thickness of the Debye layer is small compared to the tube radius.29 They employed integral formulas to obtain the coefficient of dispersion by numerical means. Their tabulated results were presented in dimensional form for several values of the ζ potential. Andreev and Lisin30,31 similarly employed numerical methods along with integral equations to obtain the coefficient of dispersion for a tube without relying on the Debye-Hu¨ckel approximation. Their dimensional results, expressed as a plate height, were plotted as a function of the applied electric field for a range of ion concentrations, ζ potentials, tube radii, and fluid transport properties. These studies also addressed species interactions with the tube wall, electrophoretic species motion, radial temperature variations, and column resolution in separation processes. These previous studies used numerical methods to obtain the coefficient of dispersion for electrokinetic transport when the ζ potential is large. Numerical methods are generally required in these cases, but the dimensional results are then applicable only to specific conditions. These previous studies thus lack the broad applicability provided by analytical solutions. In the present study, we use combined analytical and numerical methods to obtain the fluid velocity and a two-term series describing the late-time distribution of a neutral nonreacting species carried by electroosmotic flow in a tube and a channel. No a piori assumptions regarding either the Debye layer thickness or ζ potential are employed. Numerical solutions for the coefficient of dispersion and coefficient of skewness are first presented in terms of the Peclet number over a wide range of the normalized Debye layer thickness and ζ potential. These results are plotted (18) (19) (20) (21) (22) (23) (24) (25) (26) (27) (28)
Gas, B.; Stedry, M.; Kenndler, E. Electrophoresis 1997, 18, 2123-2133. Datta, R.; Kotamarthi, V. R. AIChE J. 1990, 36, 916-926. Martin, M.; Guiochon, G. Anal. Chem. 1984, 56, 614-620. McEldoon, J. P.; Datta, B. Anal. Chem. 1992, 64, 227-230. Keely, C. A.; van de Goor, T. A. A. M.; McManigill, D. Anal. Chem. 1994, 66, 4236-4242. Potocek, B.; Gas, B.; Kenndler, E.; Stedry, M. J. Chromatogr., A 1995, 709, 51-62. Griffiths, S. K.; Nilson, R. H. Anal. Chem. 1999, 71, 5522-5529. Andreev, V. P.; Dubrovsky, S. G.; Stepanov, Y. V. J. Microcolumn Sep. 1997, 9, 443-450. Zhang, X.; Regnier, F. E. J. Chromatogr., A 2000, 869, 319-328. McGuffin, V. L.; Tavares, M. F. M. Anal. Chem. 1997, 69, 152-164. Wright, P. B.; Lister, A. S.; Dorsey, J. G. Anal. Chem. 1997, 69, 32513259.
4768
Analytical Chemistry, Vol. 72, No. 20, October 15, 2000
in dimensionless form, making them generally applicable to any combination of fluid properties and the applied electric field. In addition, closed-form expressions for the fluid velocity and the coefficients of dispersion and skewness are given for the limiting cases of a small and large Debye layer thicknesses. These simple expressions are applicable to all ζ potentials and are asymptotically exact for small and large values. GOVERNING EQUATIONS Consider the electroosmotic fluid motion and neutral solute transport in a tube or a channel of infinite width. Assuming that the fluid is incompressible and that transport properties are constant, the time-dependent concentration field is governed by
∂c/∂t + u‚∇c ) D∇2c
(1)
where c is the local solute concentration, t is time, u ) ui + vj is the local fluid velocity, and D is the coefficient of diffusion. Further assuming that the flow is steady, that there are no applied pressure gradients, and that inertial effects are small, the momentum equation may be written as
µ∇2u ) Fe∇φ
(2)
where µ is the fluid viscosity, Fe is the net local charge density, and φ is the electric potential. Finally, for a dielectric constant, , that does not vary with position, the Poisson equation governing the electric potential is
∇2φ ) -Fe
(3)
and the charge density for equivalent ions may be related to the electric potential through the Boltzmann distribution, Fe ) -2Fzce sinh(zFφ/RT), where F is the Faraday constant, z is the ion charge number, ce is the bulk-fluid ion concentration far removed from any surface, R is the universal gas constant, and T is the temperature. We now introduce a set of dimensionless variables. The new normalized dependent variables are taken as c* ) c/co, u* ) u/U, and φ* ) φ/ζ, where co is some reference concentration yet to be specified, U is the mean axial fluid speed spatially averaged across the tube or channel, and ζ is the electric potential at the tube or channel wall. The new independent variables are x* ) (x - Ut)/a, y* ) y/a, and t* ) Dt/a2, where x and y are the axial and transverse coordinates and a is the tube radius or channel half-height. This normalization leads to three new parameters, the normalized Debye length,
λ*2 )
2
(aλ) ) 2FRT zca 2 2
2
(4)
e
the normalized wall potential, ζ* ) zFζ/RT, and the Peclet number, Pe ) Ua/D indicating the relative magnitudes of advective and diffusive transport rates. Introducing these normalized variables into the primitive governing equations and rearranging slightly yields
(
)
∂c* ∂c* + Pe u*‚∇c* ) ∇2c* ∂t* ∂x*
(5)
β∇2u* ) -∇2φ* E*
(6)
and
λ*2∇2φ* )
1 sinh(ζ*φ*) ζ*
(7)
for eqs 1-3, respectively. The new dependent variable, E* ) -∇φ/Ex is the electric field vector normalized by the applied axial electric field, Ex. Note that the operators ∇ and ∇2 above implicitly involve derivatives with respect to the normalized independent variables x* and y* when applied to any normalized dependent variable. The normalization additionally introduces one dimensionless unknown constant, β ) -µU/ζEx. This constant is the mean axial fluid speed normalized by the Helmholtz-Smoluchowski speed for flow past a plane charged surface. Its value is given by the condition
(n + 1)
∫
1
0
n
u*y* dy* ) 1
exists along the center of the tube or channel when λ* , 1, but no such interior region exists when λ* is order one or larger. In the latter case, Debye layers from the opposing walls overlap, and this implicit assumption would appear to break down. This is not always the case, however, and eq 11 is widely used in modeling electrokinetic phenomena for large λ*.17,30-34 When λ* is large, the electrically neutral region exists outside the tube or channel, say in an upstream reservoir, so eq 11 is still rigorously justified for large λ* provided that the Peclet number based on the Debye layer thickness, Peλ ) λ*Pe, is small.6 By analogy to the Graetz problem35 in convective transport, eq 11 should also remain valid for cases in which the tube or channel length is large compared to Peλ. Under either constraint, eq 11 yields the correct result that the potential and ion concentrations become uniform across the tube or channel when λ* becomes large. Given the governing equation (10) and parallels between the u* and φ* boundary conditions, we see that the normalized fluid velocity can be expressed directly in terms of the normalized electric potential as
u* ) (1 - φ*)/β
(12)
(8)
in accordance with the definition of the mean fluid speed, U. The parameter n in eq 8 is used to describe either the planar or axisymmetric geometries by taking n ) 0 or 1, respectively. Recognizing that the incompressible flow in a long tube or channel must be one-dimensional and that the radial component of the fluid velocity is therefore everywhere zero, the diffusionadvection eq 5 may be rewritten as
∂c* ∂c* ∂2c* ∂c* 1 ∂ + Pe(u* - 1) ) + y*n ∂t* ∂x* ∂x*2 y*n ∂y* ∂y*
(
)
(9)
Substituting this result into eq 8, the normalized mean fluid speed may be written explicitly as
β)-
µU ) (n + 1) ζEx
∫ (1 - φ*)y* 1
0
n
dy*
(13)
Note that φ* here implicitly represents only the intrinsic portion of the electric potential, as described by eq 11. The fluid velocity thus possesses no axial variation due to the applied electric field. METHOD OF SOLUTION
where u* is the normalized fluid speed in the axial direction. Now recognizing that the second derivative of the electric potential in the axial direction is small compared to that in the transverse direction and that the normalized electric field is defined such that E*‚i ) 1, the axial component of the momentum equation (6) may be written as
β d du* ∂φ* 1 ∂ y*n )- n y*n n dy* dy* ∂y* ∂y* y* y*
(
)
(
)
(10)
and the Poisson-Boltzmann equation becomes 2
∂φ* λ* ∂ 1 sinh(ζ*φ*) y*n ) ∂y* ζ* y*n ∂y*
(
)
(11)
Boundary conditions for the fluid velocities are du*/dy* ) 0 at y* ) 0 and u* ) 0 at y* ) 1. Those for the electric potential are ∂φ*/∂y* ) 0 at y* ) 0 and φ* ) 1 at y* ) 1. Note that this form of the Poisson-Boltzmann equation relies on an implicit assumption that the fluid is electrically neutral in some region of zero potential far from any surface. Such a region
To solve eq 9, we now introduce one more transformation of the axial coordinate, η ) x*/2(t*1/2). Introducing this new variable into the diffusion-advection equation yields
∂2c* ∂c* ∂c* + 2η ) 2xt*Pe(u* - 1) + 2 ∂η ∂η ∂η 1 ∂ ∂c* ∂c* 4t* y*n ∂t* y*n ∂y* ∂y*
[
(
)] (14)
In a previous analysis,24 we showed that eq 14 possesses a latetime solution of the form (29) (30) (31) (32)
Gas, B.; Stedry, M.; Kenndler, E. J. Chromatogr., A 1995, 709, 63-68. Andreev, V. P.; Lisin, E. E. Electrophoresis 1992, 13, 832-837. Andreev, V. P.; Lisin, E. E. Chromatographia 1993, 37, 202-210. Levin, S.; Marriott, J. R.; Neale, G.; Epstein, N. J. Colloid Interface Sci. 1975, 52, 136-149. (33) Oldham, I. B.; Young, F. J.; Osterle, J. F. J. Colloid Interface Sci. 1963, 18, 328-336. (34) Morrison, F. A.; Osterle, J. F. J. Chem. Phys. 1965, 43, 2111-2115. (35) Kays, W. F. Convective Heat and Mass Transfer; McGraw-Hill: New York, 1966.
Analytical Chemistry, Vol. 72, No. 20, October 15, 2000
4769
c*t*m/2 ) f0 +
(
)
df0 1 f1 + Pe g + dη 1 2xt*
(
or R1. For the higher-order terms, however, the conditions g′2(1) ) 0 and g′3(1) ) 0 require that
)
df1 d 2f 0 1 2 g + Pe f + Pe g2 + ... (15) 4t* 2 dη 1 dη2
R0 ) -(n + 1)
∫ [(u* - 1)g ]y* 1
1
0
n
dy*
(23)
and where the functions gj depend only on the transverse position, y*, and the functions fk depend only on the transformed axial position, η. Values of the parameter m ) 0 or m ) 1 identify the cases of a translating interface or instantaneous source, respectively. Substituting eq 15 into eq 14, grouping like powers of time, and separating the results into functions of the axial and transverse directions yields governing equations for the functions fk and gj. For the first few terms these are
(1 + R0Pe2) f ′′0 + 2ηf ′0 + 2mf0 ) 0
(16)
(1 + R0Pe2) f ′′1 + 2ηf ′1 + 2(m + 1)f1 ) -R1Pe3f 0′′′ (17) and
g′1 g′′1 + n - (u* - 1) ) 0 y*
(18)
g′2 g′′2 + n - (u* - 1)g1 ) R0 y*
(19)
g′3 g′′3 + n - (u* - 1)g2 - R0g1 ) R1 y*
(20)
Note that primes applied to the functions fk denote differentiation with respect to η, while those applied to gj denote differentiation with respect to y*. The separation constants or eigenvalues R0 and R1 are yet to be determined. Based on symmetry about the center line and the condition of zero solute flux through the tube or channel walls, boundary conditions for the functions g1 and g2 are g′1(0) ) g′1(1) ) g′2(0) ) g′2(1) ) 0. In addition, these functions must satisfy
gjj ) (n + 1)
∫
1
0
gj y*n dy* ) 0
(21)
for all j in order that all axial variation of the concentration field is carried by the functions fk alone. Because the governing equations (18-20) are second order, only one constant of integration is available to satisfy the two boundary conditions on the derivative g′j(y*). As a result, one of the two boundary conditions can be satisfied only by appropriate choices for the eigenvalues, Rk. From the governing equation (18) and condition g′1(0) ) 0, the derivative g′1(1) at the tube or channel wall can be written as
g′1(1) )
∫ (u* -1)y* 1
0
n
dy* ≡ 0
(22)
By eq 8, the first-order solution therefore automatically satisfies all required conditions and so imposes no constraint on either R0 4770
Analytical Chemistry, Vol. 72, No. 20, October 15, 2000
R1 ) -(n + 1)
∫ [(u* - 1)g 1
0
2
+ R0g1]y*n dy*
(24)
respectively. These eigenvalues are thus uniquely determined by the two boundary conditions g′j(0) ) g′j(1) ) 0. The two constants of integration for each function gj are then determined by the integral constraint (21) and the condition g′j(0) ) 0. Note that the term R0g1 in eq 24 arises naturally in the derivation of this expression from eq 20, but by eq 21 it does not contribute to the integral. NUMERICAL PROCEDURE As discussed later, eqs 16 and 17 describing the axial variation of the solute concentration possess closed-form analytical solutions. Equation 11, however, possesses no such solution for the electric potential in a tube or channel of finite width and nonnegligible values of ζ*. In general, eq 11 and the dependent equations describing the fluid velocity and transverse variation of the solute concentration must be solved numerically. For this we employ a very accurate shooting technique based on the following procedure. Given values of λ* and ζ*, eq 11 is first integrated from y* ) 0 to y* ) 1 using the boundary condition φ*′(0) ) 0 and an initial guess for φ*(0). The resulting value of φ*(1) does not necessarily satisfy the required boundary condition φ*(1) ) 1, so an improved value of φ*(0) is computed based on the discrepancy between the observed and required boundary values. This process in repeated to convergence, each time using the improved estimate of φ*(0) to begin the integration. These integrations are performed using a double-precision routine, DDERKF,36 with specified relative and absolute error tolerances of 10-12 and 10-20, respectively. New values of φ*(0) are computed automatically using the nonlinear algebraic equation solver DNSQE.36 The convergence criterion for the results shown here is a relative error of 10-8 in the value of φ*(0). This criterion yields an error below 10-6 in satisfying the boundary condition φ*(1) ) 1. Once the electric potential is known, the corresponding value of β is computed from eq 13. The required quadrature is performed by integrating eqs 11 and 13 together as a coupled pair, using the correct value of φ*(0) previously determined. From eq 12, the velocity profile is then known. Next, eqs 18, 19, 21, 23, and 24 are added to the system of equations and integrated along with eqs 11 and 13 to yield g1(y), g2(y), R0, and R1. This final system of equations consists of 10 coupled, first-order ordinary differential equations. These are again integrated by means of DDERKF, using the relative and absolute error tolerances described above. To satisfy the integral constraints imposed by eq 21, the system of equations is simply integrated multiple times. During the first (36) Haskell, K. H.; Vandevender, W. H.; Walton, E. L. The SLATEC Mathematical Subprogram Library: SNL Implementation; SAND80-2992, Sandia National Laboratories, Albuquerque, NM, 1980.
integration, the unknown value g1(0) is taken as zero and the equations are integrated from y* ) 0 to y* ) 1, yielding a tentative value of gj1 from eq 21. The integration is then repeated, this time setting g1(0) ) -gj1. By this method eq 21 is satisfied exactly. Once g1(y) is known, the system of equations is integrated again, with g2(0) ) 0, to yield R0 and a tentative value of gj2. One more integration, setting g2(0) ) -gj2, gives g2(y) satisfying eq 21. A final integration yields R1. This numerical procedure was checked against previously published solutions describing the electric potential in a tube for ζ* ) 2.7 and several values of λ* between 0.01 and 10.17 The two results agree within the accuracy of reading the earlier plot. The present results show similarly good agreement with previous analytical solutions based on inner and outer expansions for ζ* between 0 and 10 and λ* between 0.01 and 1.32 The procedure was further benchmarked against known analytical solutions describing the electric potential for negligibly small ζ* in both tubes and channels.6,7 In those cases, the two results agree within a relative error of 10-6 for all values of y* and values of λ* between 10-3 and 103. Finally, computed values of R0, R1, g1(y), and g2(y) were compared with known analytical solutions for a parabolic profile of the fluid speed.16,37,38 Such a profile arises in the limit λ* f ∞. The numerical solutions reproduce these analytical results within relative errors of 10-6 for both tube and channel geometries. The shooting procedure described above works well when λ* is greater than ∼10-3. For smaller values, the shooting target becomes difficult to hit since very small changes in φ*(0) then yield large changes in φ*(1). This sensitivity results from parasitic solutions to eq 11 that grow exponentially in y*. The rate of this growth increases with decreasing λ*. Still worse, such exponential growth can lead to a numerical overflow before reaching the target location, y* ) 1, unless φ*(0) is very close to the correct value. Integrating eq 11 backward from y* ) 1 helps to solve this problem but introduces new problems of comparable severity. Because of these difficulties, an alternate method is employed when λ* is very small. When λ* is small, the electric potential varies only in a region close to the tube or channel wall. The boundary condition φ*′(0) ) 0 becomes nearly superfluous in this limit, and the electric potential can be well approximated by34
φ* ∼
ζ* 4 arctanh e(y*-1)/λ* tanh ζ* 4
[
( )]
(25)
as λ* f 0. This is an exact solution to eq 11 for all ζ* and satisfies the boundary condition φ*(1) ) 1. It is approximate only in the sense that the condition φ*′(0) ) 0 is not satisfied exactly except in the limit λ* f 0. In practice, however, this solution is very accurate for all λ* less than ∼0.01. Shooting is thus unnecessary when λ* is small, and the electric potential can be described instead by this closed-form result. The values of g1(y), g2(y), R0, and R1 in this alternate method are still computed by the balance of the numerical procedure described above since all of the required integrals of φ*, as expressed by eq 25, are not known in analytical form. FLUID MOTION IN A TUBE The computed normalized electric potential and axial fluid speed are shown in Figure 1 for electroosmotic flow in a tube.
Figure 1. Normalized electric potential (φ*) and normalized fluid speed (u*) for electroosmotic flow in a tube at various ζ potentials.
Figure 2. Normalized mean fluid speed for electroosmotic flow in a tube (solid) and channel (dashed) at various ζ potentials.
These sample results are for λ* ) 0.3 and a range of values of ζ*. Here we see that increasing ζ* is somewhat analogous to decreasing λ*. In both cases, the electric potential and fluid velocity develop a steeper gradient near the tube wall, but then exhibit an increased region of nearly uniform values near the center line. We thus expect that increased ζ* will yield more uniform concentration fields and reduced late-time solute dispersion. Figure 2 shows computed values of the normalized mean fluid speed, β, corresponding to the conditions of Figure 1. Here we see that increasing ζ* always increases β. This is again because larger ζ* yields an effectively smaller Debye layer thickness and so yields a normalized electric potential closer to zero over a broader portion of the tube cross section. By eq 13, β is thus increased. We also see that β approaches unity as λ* becomes small and falls as 1/λ*2 for all ζ* as λ* grows large. In addition to numerical results such as those of Figures 1 and 2, we have obtained analytical solutions for the fluid velocity in the extremes of small and large λ*. Such solutions serve to quantify the general observations above and further provide valuable benchmark results for testing more general finiteAnalytical Chemistry, Vol. 72, No. 20, October 15, 2000
4771
difference and finite-element methods. Recall from eq 12 that the local fluid speed under all conditions is given by u* ) (1 - φ*)/ β. Thus, local axial fluid speeds, u*, are uniquely determined by specifying only the normalized electric potential, φ*, and the normalized mean axial fluid speed, β. These are given below in analytical form for small and large λ*. In the important practical limit of small λ*, the electric potential is given by eq 25 and the corresponding normalized mean fluid speed is
β∼1-
4λ* [Li (ξ) - Li2(-ξ)] ζ* 2
as
λ* f 0
(26)
computed using eq 31 agree with those of eq 26 within 3% for ζ* g 10, and within 1% for ζ* g 14. Analytical solutions to eqs 11 and 13 can also be obtained in the limit of large λ*. In this extreme, the potential and normalized mean fluid speed are
φ* ∼ 1 -
sinhζ* (1 - y*2) 4ζ*λ*2
as
λ* f ∞
(32)
and
β∼
1 sinhζ* 8λ*2 ζ*
as
λ* f ∞
(33)
where ξ ) tanh(ζ*/4) and Li2(ξ) is the dilogarithm function,
Li2(ξ) ) -
∫
ξ
0
ln(1 - s)
ds s
∞
)
∑ k)1
ξk
(27)
k2
This result, correct to order λ*, was obtained by integrating eq 25 in accordance with eq 13 defining the normalized mean fluid speed. It agrees with the value of β computed numerically within 5% for all ζ* when λ* is less than 0.5; when λ* is less than 0.2, the two agree within 1%. The asymptotic behavior of eq 25 in the limit of small ζ* is
(
φ* ∼ e(y*-1)/λ* 1 -
)
ζ*2 [1 - e2(y*-1)/λ*] 48
(28)
as λ* f 0 and ζ* f 0. The corresponding normalized axial speed in this double limit is
(
β ∼ 1 - 2λ* 1 -
)
ζ*2 72
as
(
)
(
)]
(30)
as λ* f 0 and ζ* f ∞. This yields a corresponding normalized mean fluid speed of
π2 β ∼ 1 - λ* ζ*
as
λ* f 0, ζ* f ∞
(31)
This expression for β can be obtained directly from eq 26 since ξ f 1 as ζ* f ∞ and Li2(1) - Li2(-1) ) π2/4. Values of 1 - β 4772
Analytical Chemistry, Vol. 72, No. 20, October 15, 2000
(
)
1 ζ*2 1 + 6 8λ*2
as
λ* f ∞, ζ* f 0
(34)
Again this is consistent with the result β ) 1-2λ*I1(1/λ*)/I0(1/ λ*) ∼ 1/8λ*2 for λ* f ∞ and ζ* ) 0. The expansion of eq 33 for large ζ* is
β∼
2 y* - 1 -ζ*/2 y* - 1 ln sech2 e - tanh ζ* 2λ* 2λ*
[
β∼
λ* f 0, ζ* f 0 (29)
This expression for β is consistent with our previous result24 of β ) 1 - 2λ*I1(1/λ*)/I0(1/λ*) ∼ 1 - 2λ* for the case of vanishingly small ζ*. The values of 1 - β computed using eq 29 agree with those of eq 26 within 3% for ζ* e 2 and within 15% for ζ* e 4. Note that the mean normalized fluid speed varies weakly with ζ* when ζ* is small. The term ζ*2/72 yields a correction to 1 - β of only ∼5% for ζ* ) 2. Thus, by eq 13, the dimensional mean fluid speed should grow about linearly with ζ for small ζ* when λ* is small. In the alternate extreme of large ζ*, the asymptotic behavior of eq 25 is
φ* ∼ -
These expressions were obtained by solving eq 11 for large λ* using a perturbation expansion in 1/λ* and then using the resultant electric potential to compute β from eq 13. As a result, they carry the additional restriction that sinhζ*/ζ* , λ*2. Equation 33 reproduces the values of β determined numerically within ∼2% for all λ* g 3, provided that sinhζ*/ζ* < λ*2/50. Note that the fluid velocity in this limit of large λ* recovers the expected parabolic profile u* ) (1 - φ*)/β ) 2(1 - y*2) for pressure-driven flow in a tube. The asymptotic behavior of eq 33 for small values of ζ* is
1 eζ* 8λ*2 2ζ*
as
λ* f ∞, ζ* f ∞
(35)
We note that the effect of ζ* on the normalized mean fluid speed at large λ* is much stronger than it is when λ* is small. This is especially the case as ζ* becomes large. In this limit, the dimensional fluid speed increases exponentially with increasing ζ so long as eζ*/2ζ* remains much smaller than λ*2. FLUID MOTION IN A CHANNEL Results paralleling those of Figure 1 for a tube are shown for a channel in Figure 3. Normalized mean fluid speeds for flow in a channel appear in the previous Figure 2 as dashed curves. These results are very similar to those for a tube, except that β for a channel is everywhere a bit higher than that for a tube at the same values of λ* and ζ*. As for a tube, local fluid speeds in a channel are given by u* ) (1 - φ*)/β. The flow field is thus again determined uniquely by the normalized electric potential, φ*, and normalized mean speed, β. In the limit λ* f 0, the electric potential given in eq 25 applies equally to a tube and a channel since the Debye layer thickness is then very small compared to the tube diameter or channel width. Despite this, the value of β differs between the two geometries owing to the difference in eq 13 for the two cases. For a channel,
β∼
1 sinhζ* 3λ*2 ζ*
λ* f ∞
as
(40)
As with eqs 32 and 33, these results were obtained by solving eq 11 for large λ* using a perturbation expansion and so again carry the restriction that sinhζ*/ζ* , λ*2. Equation 40 reproduces the numerically determined values of β within ∼4% for all λ* g 4 provided that sinhζ*/ζ* < λ*2/100. Note that these expressions recover the parabolic profile u* ) (1 - φ*)/β ) 3(1 - y*2)/2 for pressure-driven flow in a channel. The asymptotic behavior of eq 40 for small values of ζ* is
β≈
Figure 3. Normalized electric potential (φ*) and normalized fluid velocity (u*) for electroosmotic flow in a channel.
(
2λ* [Li (ξ) - Li2(-ξ)] ζ* 2
β∼ as
λ* f 0
λ* f ∞, ζ* f 0
as
(41)
Again, this is consistent with β ) 1 - λ* tanh(1/λ*) ∼ 1/3λ*2 for λ* f ∞ and ζ* f 0. The corresponding expansion of eq 40 for large ζ* is
the normalized mean fluid speed for small λ* is
β∼1-
)
1 ζ*2 1+ 2 6 3λ*
(36)
1 eζ* 3λ*2 2ζ*
λ* f ∞, ζ* f ∞
as
(42)
The additional restriction that sinhζ*/ζ* , λ*2 also applies here. where again ξ ) tanh(ζ*/4) and Li2(ξ) is the dilogarithm function given by eq 27. As before, this result was obtained using the analytical expression for the electric potential given by eq 25. The accuracy and applicable range of this result are comparable to the parallel result for a tube, given by eq 26. In the limit of small λ* and small ζ*, the electric potential in a channel is again the same as that in a tube and is given by eq 28. The corresponding asymptotic behavior of eq 36 in this small ζ* limit is
(
β ∼ 1 - λ* 1 -
)
ζ*2 72
as
λ* f 0, ζ* f 0 (37)
AXIAL VARIATION OF THE CONCENTRATION FIELD The axial variation of the mean solute concentration, spatially averaged across the tube or channel cross section, is governed by eqs 16 and 17. These equations were previously solved for initial and boundary conditions describing two special cases of practical importance.24 The first case is the transport in the vicinity of a traveling solute interface, as indicated by m ) 0. The second is the transport of a solute peak following injection of a instantaneous planar source into the fluid stream. This problem is indicated by m ) 1. For the translating interface, the results are
f0 ) This expression for a channel is again consistent with our previous solution β ) 1 - λ* tanh(1/λ*) ∼ 1 - λ* for vanishing ζ*.24 The asymptotic behavior of eq 36 at large ζ* is
β ∼ 1 - λ*
π2 2ζ*
as
λ* f 0, ζ* f ∞
φ* ∼ 1 -
sinhζ* (1 - y*2) 2 2λ* ζ*
as
λ* f ∞
(39)
f1 )
R1Pe3 2xπ(1 + R0Pe2)3
(1 - 2η′2)e-η′
2
(44)
where
η′ )
η
(45)
x1 + R0Pe2
For the case of an instantaneous planar source of unit strength, the zero- and first-order functions are
f0 ) and
(43)
and
(38)
The corresponding potential in this limit is again given by eq 30, which applies to both the tube and channel geometries when λ* is small. The accuracy and applicable range of eqs 37 and 38 are comparable to those of the corresponding results for a tube. Analytical solutions to eqs 11 and 13 have also been obtained for flow in a channel at large λ*. In this case, the results are
1 erfc η′ 2
1
e-η′
2
2xπ(1 + R0Pe )
(46)
2
Analytical Chemistry, Vol. 72, No. 20, October 15, 2000
4773
Figure 4. First-order correction describing the transverse variation in solute concentration for electroosmotic flow in a tube.
Figure 5. Eigenvalues for flow in a tube. Zero-order eigenvalue R0 is the coefficient of dispersion; R1 is coefficient of skewness.
and
midplane of the solute peak but are antisymmetrically increased on the center line ahead. Concentrations at the tube walls exhibit just the opposite behavior. We also see in Figure 4 that increasing values of the ζ potential tend to reduce the transverse variation of the concentration field. This is consistent with the previous observation that increasing ζ* produces a more uniform fluid velocity over a large central portion of the tube cross section. Computed values of R0 and R1 are shown in Figure 5. These values along with eqs 43-47 provide the complete two-term solution describing axial variation of the mean concentration field in the late-time limit. Here we see that R0 is positive for all values of λ* and ζ* and that increasing ζ* always decreases its value. Since R0 is the same as the coefficient of dispersion, solute dispersion is thus reduced when ζ* is increased. This is again due to the fact that increasing ζ* produces a steeper gradient of the fluid velocity at the tube wall but also gives a more uniform velocity near the centerline. Increasing ζ* also decreases the magnitude of R1 for all λ*, which by eqs 44 and 47 reduces the late-time skewness of the mean solute concentration. Note that values of R1 are negative on the right of Figure 5 but become positive below some critical value of λ* that depends on ζ*. For the problem of a translating solute interface, the centroid of the mean concentration field is thus shifted forward of its asymptotic position for small values of λ* but is shifted rearward of this position for large λ*. The centroid of the concentration field is similarly shifted from its asymptotic position for the case of an instantaneous planar source. For the special case in which R1 vanishes, the late-time solutions exhibit no asymmetry except that due to still higher-order terms. Although the eigenvalues in Figure 5 exhibit a complex dependence on λ* and ζ*, they follow well-defined asymptotes in the limits of both small and large λ*. When λ* is large, the charge distribution and electric potential are nearly uniform across the tube. In this case, the applied electric field produces a nearly uniform body force that acts on the fluid in manner analogous to the pressure gradient in pressure-driven flows. The result in both cases is a parabolic profile of the fluid velocity, having a maximum speed u* ) 2 that occurs on the center line of the tube. In this limit of large λ*, our present solutions recover the well-known
f1 )
R1Pe3 2xπ(1 + R0Pe )
(3 - 2η′2)η′e-η′
2 2
2
(47)
where η′ is again given by eq 45. Noting the definition of η given above, we see that eqs 43 and 46 describe a diffusion process in which the molecular diffusivity, D, is replaced by an effective diffusivity D′ ) D(1 + R0Pe2). The eigenvalue R0 is thus equivalent to the more familiar coefficient of dispersion. The eigenvalue R1 is also equivalent to the coefficient of skewness. The solutions for f0 and f1 describing the axial variation of the concentration field are applicable to both the tube and channel geometries and to all profiles of the fluid velocity. The geometry and velocity distribution influence these solutions only through the eigenvalues R0 and R1, which are determined by solving the equations describing the transverse variation of the solute concentration. TRANSVERSE VARIATION OF THE CONCENTRATION FIELD Figure 4 shows computed values of the function g1(y*) describing the first-order term of the transverse variation of the concentration field. These results are useful in understanding the effects of a nonuniform fluid velocity on the concentration field. Recall from eq 15 that the late-time concentration field is given by c*t*m/2 ) f0 + (f1 + Pe f ′0g1)/2xt* + O[1/t*]. Thus, the function g1 in the concentration field is multiplied by both the Peclet number and the derivative of f0. For the case of a moving interface having a high concentration on the trailing side, the derivative of f0 is everywhere negative. Thus, from Figure 4 we see that the electroosmotic velocity profile increases the concentration on the tube center line, y* ) 0, at all axial positions. Likewise, the concentration at the tube wall is everywhere reduced. For the case of an instantaneous source, however, the derivative of f0 is positive for η < 0 and negative for η > 0. In this case, concentrations on the center line are reduced behind the 4774 Analytical Chemistry, Vol. 72, No. 20, October 15, 2000
results of Taylor16 and Chatwin37 for pressure-driven flow in a tube. These are
R0 ∼
1 48
R1 ∼ -
and
1 2880
as
λ* f ∞
(48)
The applicable range of these results is λ* > 2 and ζ*2 < λ*. Over this range, the discrepancy between eq 48 and the numerical results of Figure 5 is less than 2% for R0 and less than 3% for R1. In many practical applications of electroosmotic flow, the normalized Debye length is very small. In the limit, the results of Figure 5 can be expressed as
R0 ∼
1 2 λ* h0 2
as
λ* f 0
(49)
where the factor λ*2/2 is the asymptotic value of R0 for λ* f 0 and negligible ζ*.24 The function h0 depends only on ζ* and is well-approximated by
h0 ≈
2592 + 24ζ*2 2592 + 96ζ*2 + ζ*4
(50)
This function is an empirical fit to the numerical results of Figure 5, taking into account the asymptotic behavior at small and large values of ζ*. These limiting behaviors are h0 ) 1 - ζ*2/36 as ζ* f 0 and h0 ) π4/4ζ*2 ≈ 24/ζ*2 as ζ* f ∞. Equation 49, with h0 given by eq 50, agrees with the numerical results within 4% for values of λ* less than 0.01 and all values of ζ*. Equations 49 and 50 are also in good agreement with the numerical results of Gas et al. for transport in a tube at very small λ*.29 In the present notation, they obtained the result R0 ) λ*2(B/λ)2/2 for λ* , 1. Here, B is the effective thickness of the boundary layer, given by β in their text. Thus, in terms of their result, h0 ) (B/λ)2. Computed values of B (and λ) are given in Table 1 of their paper for a range of ζ potentials between 20 and 200 mV (0.8 e ζ* e 8.0). All these values are consistent with h0 ) (B/λ)2 to within ∼5% when h0 is given by eq 50. Note that the h0 correction to eq 49 is over a factor of 3 at ζ* ) 8. The behavior of R1 in the limit of small λ* similarly can be expressed as the product of the solution for ζ* f 0 and a function accounting for the effect of nonnegligible ζ*. In this case, the result is
R1 ∼
1 λ*3h1 12
as
λ* f 0
(51)
where
h1 ≈
5760 + 120ζ*2 5760 + 360ζ*2 + 12ζ*3 + ζ*4 + ζ*5
(52)
As with h0, the function h1 represents an empirical fit to the numerical results. Here the asymptotic behaviors are h1 ) 1 ζ*2/24 as ζ* f 0, and h1 ) π6/8ζ*3 ≈ 120/ζ*3 as ζ* f ∞. The applicable range of eq 51 is again roughly λ* < 0.01. For all smaller (37) Chatwin, P. C. J. Fluid Mech. 1970, 43, 321-352.
Figure 6. First-order correction describing the transverse variation in solute concentration for electroosmotic flow in a channel.
λ* and all ζ*, eqs 51 and 52 yield values for R1 that are within 20% of the values computed numerically. This relative error falls to 2% for λ* < 0.001. In this latter range, all inaccuracy in eq 51 is attributable to the approximation π6/8 ≈ 120 in the expression for h1. Results describing transverse variation of the concentration field in a channel closely parallel those for a tube presented above. Figure 6 shows computed values of the function g1(y*) for λ* ) 0.3 and values of ζ* between 0 and 12. We see that the behavior of this function is qualitatively very similar to that shown in Figure 4 for a tube. The main difference is that the intercept g1(y*) ) 0, indicating a neutral point at which the solute concentration is neither increased nor reduced, is shifted closer to the midpoint between the center line and wall for the case of a channel. This occurs because of the constraint gj1 ) 0 given by eq 21. The radial geometry of the tube shifts this neutral point toward the wall, close to the median of the cross-section area, while the planar geometryof the channel requires that the neutral point occur somewhat near the midpoint for any function that is reasonably antisymmetric about the intercept. The previous comments concerning regions of solute enrichment and reduction apply also to the channel geometry. The eigenvalues R0 and R1 for a channel are shown in Figure 7. These were computed using the same numerical procedure previously described for a tube. We see that the two geometries yield very similar results for R0 but that R1 differs significantly between the two. In the case of the tube, this coefficient of skewness changes signs from positive to negative in the vicinity of λ* ) 1. The results for a channel show no such behavior; R1 is everywhere positive. In the limit λ* f ∞, the values of R0 shown in Figure 7 recover the known coefficient of dispersion for pressure-driven flow in a channel of infinite width.38 We find no published record of R1, the coefficient of skewness, for this geometry. These asymptotic (38) Aris, R. Proc. R. Soc. London 1959, A 252, 538-550.
Analytical Chemistry, Vol. 72, No. 20, October 15, 2000
4775
Peak spreading by diffusion and dispersion is often expressed in terms of a theoretical plate height, H, given by H ) σ2/L, where L ) Ut is the column length or length of run of the species band. In terms of the present nomenclature, this can be expressed as
H)
σ2 2a ) (1 + R0Pe2) L Pe
(57)
Again, a is the channel half-height or tube radius. Finally, the extent of spreading is also sometimes described by a plate number, N, where N ) L/H. This is equivalent to
N)
Figure 7. Zero- and first-order eigenvalues for flow in a channel. Increasing ζ* always reduces these coefficients of dispersion and skewness.
values for the channel are
R0 ∼
2 105
and
R1 ∼
4 17325
as
λ* f ∞
(53)
In the opposing limit of small λ*, we find that the eigenvalues for a channel are
R0 ∼
1 2 λ* h0 3
as
λ* f 0
(54)
R1 ∼
4 λ*3h1 45
as
λ* f 0
(55)
and
where h0 and h1 are the same as those given above for a tube. The functions h0 and h1 are the same for both the tube and channel geometries because they apply only in the limit of very small λ*. In this limit, both geometries appear to be planar over the small region in which the electric potential varies. The fluid velocities in the two geometries are thus the same in this limit, at least within a constant, and the values of R0 and R1 must therefore share the same dependence on ζ* for the two geometries. VARIANCE, PLATE HEIGHT, AND PLATE NUMBER The eigenvalues presented above determine the axial variation of a concentration peak in accordance with eqs 15, 46. and 47. At sufficiently late times, the mean species distribution along the channel is thus Gaussian, and the variance of this distribution grows linearly in time. This is given by
σ2 ) 2D(1 + R0Pe2)t
(56)
The term 2Dt in eq 56 accounts for spreading due to diffusion alone, while the term 2R0Pe2Dt accounts for dispersion. 4776
Analytical Chemistry, Vol. 72, No. 20, October 15, 2000
Pe L2 L L ) ) H σ 2a 1 + R Pe2
()
(58)
0
Spreading over a fixed length of run is minimized when the plate number is maximum. Based on eq 58, this occurs for the condition Pe ) x1/R0 for a fixed length L and fixed tube or channel dimension a. SUMMARY We have examined the electroosmotic fluid motion and latetime solute transport in a tube and in a channel for cases in which the ζ potential is not necessarily small. Using both numerical and analytical methods, the transverse variation of the electric potential and fluid speed were computed over a broad range of the normalized Debye layer thickness, λ*, and the normalized ζ potential, ζ*. These fluid velocities were then used to compute the late-time distribution of a neutral nonreacting solute carried in the flow. The numerical procedure used here is based on a shooting method. This procedure was checked against previously published analytical and numerical solutions for both large and small values of the ζ potential. In comparisons with analytical results for a small ζ potential, the two agree within a relative error of 10-6 for all values of the transverse position and values of λ* between 10-3 and 103. Analytical solutions were also obtained for the fluid velocity in a tube and channel in the asymptotic limits of small and large λ*. The first of these applies when the Debye layer thickness is small compared to the transverse tube or channel dimension and is valid for all values of the normalized ζ potential. The second is applicable to large values of the Debye layer thickness. These solutions are again valid for all ζ* but carry the restriction that sinhζ*/ζ* < λ*2. The accuracy and applicable range of each solution is discussed, and expansions of the solutions for large and small ζ* are provided. Based on these fluid velocities and a series describing the full late-time concentration field, the mean axial variation of the field was determined in closed form. The zero- and first-order solutions describing this variation contain two unknown constants that arise as eigenvalues in the series solution. These eigenvalues are the coefficients of axial dispersion and skewness. They are determined numerically along with the solutions describing transverse variation of the concentration field. The numerical procedure reduces the governing equations to a system of 10 coupled first-order ordinary differential equations, which are integrated in the transverse direction using a standard integration routine.
The results presented here recover the well-known solutions for dispersion in pressure-driven flows when the Debye length is sufficiently large. In this limit, the axial dispersion is proportional to the square of the Peclet number based on the transverse dimension of the tube or channel and is independent of the ζ potential. The skewness in this limit is proportional to the cube of this Peclet number and is also independent of the ζ potential. In the limit of a small Debye layer thickness, we find that the dispersion varies as the square of the Peclet number based on the Debye length. The skewness varies as the cube of this Peclet number. In this limit, both the dispersion and skewness exhibit a first-order dependence on the ζ potential. Simple approximations to the dependence are presented. We find that increasing values of the ζ potential always reduce both late-time dispersion and skewness. ACKNOWLEDGMENT The authors thank Dr. R. S. Larson for his very careful review of this paper. This work was funded by a Sandia Phenomenological Modeling and Engineering Simulations LDRD. Sandia National Laboratories is operated by Sandia Corp. for the United States Department of Energy. NOMENCLATURE
u
axial fluid speed
u
local fluid velocity
U
mean fluid speed
x
axial position
y
transverse position
z
charge number
Rk
eigenvalues
β
normalized mean fluid speed (β ) -µU/ζEx)
dielectric constant
µ
viscosity
λ
Debye length
Fe
charge density
φ
electric potential
ζ
surface electric potential
E*
electric field: E* ) -∇φ/Ex
t*
time: t* ) Dt/a2
u*
local fluid speed: u* ) u/U
NORMALIZED VARIABLES
x*
axial position: x* ) (x - Ut)/a
a
tube radius or channel half-height
y*
transverse position: y* ) y/a
c
solute concentration
η
axial position: η ) x*/2xt*
ce
ion concentration
λ*
Debye length: λ* ) λ/a
D
diffusivity
ζ*
surface potential: ζ* ) zFζ/RT
D′
effective diffusivity inclusive of dispersion
Ex
applied axial electric field: Ex ) -dφ/dx
fk
axial concentration functions
SUBSCRIPTS AND SUPERSCRIPTS j,k
order of solution
F
Faraday constant
gj
transverse concentration functions
*
asterisk denotes normalized variable
L
distance of travel (L ) Ut)
h
bar denotes spatial average
m
interface (m ) 0) or plane source (m ) 1)
ˆ
hat denotes initial distribution
n
channel (n ) 0) or tube (n ) 1)
Pe
Peclet number: Pe ) Ua/D
R
ideal gas constant
t
time
Received for review May 11, 2000. Accepted July 31, 2000.
T
temperature
AC000539F
Analytical Chemistry, Vol. 72, No. 20, October 15, 2000
4777
