Numerical Analysis of Electroosmotic Flow in Dense Regular and

Langmuir 2005, 21, 6097-6112

6097

Numerical Analysis of Electroosmotic Flow in Dense Regular and Random Arrays of Impermeable, Nonconducting Spheres Dzmitry Hlushkou,† Andreas Seidel-Morgenstern,†,‡ and Ulrich Tallarek*,† Institut fu¨ r Verfahrenstechnik, Otto-von-Guericke-Universita¨ t Magdeburg, Universita¨ tsplatz 2, 39106 Magdeburg, Germany, and Max-Planck-Institut fu¨ r Dynamik komplexer technischer Systeme, Sandtorstrasse 1, 39106 Magdeburg, Germany Received January 27, 2005. In Final Form: April 11, 2005 We present a numerical scheme for analyzing steady-state isothermal electroosmotic flow (EOF) in three-dimensional random porous media, involving solution of the coupled Poisson, Nernst-Planck, and Navier-Stokes equations. While traditional finite-difference methods were used to resolve the PoissonNernst-Planck problem, the (electro)hydrodynamics has been addressed with high efficiency using the lattice-Boltzmann method. The developed model allows simulation of electrokinetic transport under most general conditions, including arbitrary value and distribution of electrokinetic potential at the solidliquid interface, electrolyte composition, and pore space morphology. The approach provides quantitative information on a spatial distribution of simulated velocities. This feature was utilized to characterize EOF fields in regular and random, confined and bulk packings of hard (i.e., impermeable, nonconducting) spheres. Important aspects of pore space morphology (sphere size distribution), surface heterogeneity (mismatch in electrokinetic potentials at confining wall and sphere surface), and fluid phase properties (electrical double layer thickness) were investigated with respect to their influence on the EOF dynamics over microscopic and macroscopic spatial domains. Most important is the observation of a generally nonuniform pore-level EOF velocity profile in the sphere packings (even in the thin double layer limit) which is caused by pore space morphology and which is in contrast to the pluglike velocity distribution in a single, straight capillary under the same conditions.

1. Introduction Electrokinetic transport of bulk liquid and solutes through high-surface-area materials plays a central role in many analytical, technological, and environmental processes including dewatering of waste sludge and soil remediation, electrophoresis, electrochromatographic separations in particulate and monolithic beds, and microfluidic devices.1-15 Mass transport is realized by electromigration of ions (background electrolyte), electro* To whom correspondence should be addressed. Fax: +49-(0)391-67-12028. E-mail: [email protected]. † Otto-von-Guericke-Universita ¨ t. ‡ Max-Planck-Institut. (1) Dittmann, M. M.; Wienand, K.; Bek, F.; Rozing, G. P. LC•GC 1995, 13, 800-814. (2) Acar, Y. B.; Ozsu, E. E.; Alshawabkeh, A. N.; Rabbi, M. F.; Gale, R. J. CHEMTECH 1996, 26, 40-44. (3) Crego, A. L.; Gonza´lez, A.; Marina, M. L. Crit. Rev. Anal. Chem. 1996, 26, 261-304. (4) Coletta, T. F.; Brunell, C. J.; Ryan, D. K.; Inyang, H. I. J. Environ. Eng. 1997, 123, 1227-1233. (5) Colo´n, L. A.; Burgos, G.; Maloney, T. D.; Cintro´n, J. M.; Rodrı´guez, R. L. Electrophoresis 2000, 21, 3965-3993. (6) Whitesides, G. M.; Stroock, A. D. Phys. Today 2001, 54, 42-48. (7) Tanaka, N.; Kobayashi, H.; Nakanishi, K.; Minakuchi, H.; Ishizuka, N. Anal. Chem. 2001, 73, 420A-429A. (8) Stone, H. A.; Kim, S. AIChE J. 2001, 47, 1250-1254. (9) Tallarek, U.; Rapp, E.; Seidel-Morgenstern, A.; Van As, H. J. Phys. Chem. B 2002, 106, 12709-12721. (10) Delgado, A. V., Ed. Interfacial Electrokinetics and Electrophoresis; Marcel Dekker: New York, 2002. (11) Legido-Quigley, C.; Marlin, N. D.; Melin, V.; Manz, A.; Smith, N. W. Electrophoresis 2003, 24, 917-944. (12) Jemere, A. B.; Oleschuk, R. D.; Harrison, D. J. Electrophoresis 2003, 24, 3018-3025. (13) Broyles, B. S.; Jacobson, S. C.; Ramsey, J. M. Anal. Chem. 2003, 75, 2761-2767. (14) Chen, G.; Tallarek, U. Langmuir 2003, 19, 10901-10908. (15) Stone, H. A.; Stroock, A. D.; Ajdari, A. Annu. Rev. Fluid Mech. 2004, 36, 381-411.

phoresis (charged solutes or colloidal particles), and electroosmosis (bulk liquid) driven by shear stresses concentrated in the electrical double layer (EDL) close to a solid-liquid interface.10 The microscopic and macroscopic transient behavior and the magnitude, stability, and uniformity of the electroosmotic flow (EOF) in porous media are inherently related to the physicochemical nature of the surface, the pore space morphology, and properties of the background electrolyte. An analysis of these parameters on relevant spatiotemporal scales is important because it critically guides the performance and design strategy of an electrokinetic process with respect to alternative diffusive-convective transport schemes. Concerning bulk transport of liquid through fixed beds of fine particles and other granular or fibrous media as well as microfluidic devices, electroosmosis can offer distinct advantages over pressure-driven (hydraulic) flow. The EOF is generated by interaction of an externally applied electrical field with that region of the liquid electrolyte which has become locally charged at the interface to the stationary (and oppositely charged) solid surface of the confining medium.10 The extension into the bulk solution of the fluid-side domain of the quasiequilibrium EDL (usually characterized by the Debye screening length λD) can be as small as a few nanometers as compared to channel sizes of micrometer dimension. This has important consequences for the EOF dynamics under these conditions. First, from a macroscopic point of view, the liquid in this thin-double-layer (TDL) limit apparently slips at the wall,16,17 being in contrast to the parabolic velocity profile in Poiseuille flow which results directly from the radial distribution of shear stress over the whole channel. Second, because the ratio of electroos(16) Rice, C. L.; Whitehead, R. J. Phys. Chem. 1965, 69, 4017-4024. (17) Gross, R. J.; Osterle, J. F. J. Chem. Phys. 1968, 49, 228-234.

10.1021/la050239z CCC: $30.25 © 2005 American Chemical Society Published on Web 05/18/2005

6098

Langmuir, Vol. 21, No. 13, 2005

Hlushkou et al.

motic to hydraulic volumetric flow rates (at fixed electrical potential and pressure gradients) is inversely proportional to the squared channel diameter, EOF becomes increasingly effective in liquid transport through finer channels as their size is reduced. Thus, a benefit of using EOF is that chemical and biological species can be transported over relatively long distances, while the hydrodynamic dispersion can be limited almost to that by axial diffusion alone which has been verified theoretically18 and also shown experimentally.19,20 Applications of these EOF features are found in liquid and solute transport through porous media with low hydraulic permeability, in cases where still a high volumetric throughput (not only by the preferential channelling via larger pores) and little flow field dispersion are required, as for the efficient separation of target molecules in complex bioanalytical and pharmaceutical samples. In this respect, the use of EOF in capillary electrochromatography (CEC) allows achievement of much higher separation efficiencies than can be realized by conventional liquid chromatography (LC) employing hydraulic flow.1,3 In particular, the relatively high electroosmotic permeabilities allow the packing of small porous particles as adsorbent stationary phases which improve column efficiency and alleviates the use of prohibitively high operating pressures in liquid chromatography.21 In CEC, 50-150 µm i.d. fused-silica capillaries are usually packed with 3-10 µm spherical, porous particles and electrical field strengths up to 105 V/m are applied to drive the liquid and solutes through the fixed beds. The mobile phase is a 0.5-10 mM buffer solution giving λD of the order of 1-10 nm10 being often so much smaller than the dimensions in the void space between particles that the TDL limit is realized. For the design of porous materials (the CEC and microchip columns, in particular) with improved performance, it is important to evaluate the EOF velocity distributions a priori and reach optimization with respect to hydrodynamic dispersion by adjusting operating conditions (e.g., applied electrical field or ionic strength and pH of the mobile phase), column properties such as the column diameter, particle size distribution, and pore space morphology, or surface properties. Thus, theories and experimental approaches need to resolve pore-level EOF in three-dimensional random porous media to address microscopic and mesoscopic details of the velocity field underlying flow and dispersion on the macroscopic scale. A number of theoretical studies have been reported and used for investigation of EOF in open microchannel systems22-34 as well as

packed beds.35-37 However, since the EOF problem has no general solution, it can be solved analytically only for a limited set of simple configurations such as a flat solidliquid interface or an open, homogeneous channel. Most of the theoretical approaches used to study the EOF in random porous media like sphere packings are based on a representation of the materials as an assembly of individual channels which only allows finding the average velocity through a material. Rathore and Horva´th35 have extended the model of Overbeek and Wijga38 for EOF through capillary columns in CEC (which are usually packed with spherical-shaped adsorbent particles) to account for the effect of differences in electrokinetic potential (shear plane or ζ-potential) at the capillary inner wall and the particles external surface. They derived an expression for the EOF velocity in dependence of radial position with respect to the capillary center. According to that model, if the ζ-potential of the capillary wall is identical to that of the support material, the EOF velocity becomes constant over the whole column cross section, while differences in ζ-potential result in a decrease or increase, depending on the actual ζ-potential ratio, of local axial EOF velocity in a near-wall region of the bed (electrokinetic wall effect). It is noteworthy that this theoretical model does not consider the fluctuations in interstitial porosity and tortuosity of a packed bed caused by a more ordered (layered) geometrical structure in vicinity of the confining wall. At the same time, detailed experimental measurements of pressure-driven flow in packed cylinders39,40 demonstrate a radial dependence of axial flow velocity on local packing porosity. Thus, it can be expected that a similar correlation also prevails for the EOF as indicated by dynamic NMR microscopy measurements.41 Besides theoretical approaches the EOF problem has been extensively studied numerically.42-50 The majority of numerical studies is concerned with EOF in open microchannel systems. Until now, direct numerical simulation of EOF in porous media is a challenging task. Apart from the necessity to resolve, in general, the coupled hydrodynamic, electrostatic, and mass transport problems subjected to complex geometrical boundary conditions represented by the solid-liquid interface in random porous media, simulations usually must be accomplished at different length scales ranging from the Debye length of typically 1-10 nm to the characteristic dimension of the whole system (for instance, the complete column diameter) of hundreds of micrometers. These widely disparate length scales require huge computational resources. Coelho et

(18) Griffiths, S. K.; Nilson, R. H. Anal. Chem. 1999, 71, 5522-5529. (19) Paul, P. H.; Garguilo, M. G.; Rakestraw, D. J. Anal. Chem. 1998, 70, 2459-2467. (20) Tallarek, U.; Rapp, E.; Scheenen, T.; Bayer, E.; Van As, H. Anal. Chem. 2000, 72, 2292-2301. (21) Chen, G.; Pacˇes, M.; Marek, M.; Zhang, Y.; Seidel-Morgenstern, A.; Tallarek, U. Chem. Eng. Technol. 2004, 27, 417-428. (22) Ajdari, A. Phys. Rev. E 1996, 53, 4996-5005. (23) Griffiths, S. K.; Nilson, R. H. Anal. Chem. 2000, 72, 5473-5482. (24) Tsao, H.-K. J. Colloid Interface Sci. 2000, 225, 247-250. (25) Griffiths, S. K.; Nilson, R. H. Anal. Chem. 2000, 72, 4767-4777. Cummings, E. B.; Griffiths, S. K.; Nilson, R. H.; Paul, P. H. Anal. Chem. 2000, 72, 2526-2532. (26) Santiago, J. G. Anal. Chem. 2001, 73, 2353-2365. (27) Dutta, P.; Beskok, A. Anal. Chem. 2001, 73, 1979-1986. (28) Li, D. Colloids Surf., A 2001, 195, 35-57. (29) Ghosal, S. J. Fluid Mech. 2002, 459, 103-128. (30) Qian, S.; Bau, H. H. Anal. Chem. 2002, 74, 3616-3625. (31) Gleeson, J. P. J. Colloid Interface Sci. 2002, 249, 217-226. (32) Hsu, J.-P.; Kao, C.-Y.; Tseng, S.; Chen, C.-J. J. Colloid Interface Sci. 2002, 248, 176-184. (33) Yang, J.; Masliyah, J. H.; Kwok, D. Y. Langmuir 2004, 20, 38633871. (34) Xuan, X.; Li, D. J. Micromech. Microeng. 2004, 14, 290-298.

(35) Rathore, A. S.; Horva´th, Cs. J. Chromatogr. A 1997, 781, 185195. (36) Vallano, P. T.; Remcho, V. T. Anal. Chem. 2000, 72, 4255-4265. (37) Liapis, A. I.; Grimes, B. A. J. Chromatogr. A 2000, 877, 181215. (38) Overbeek, J. Th. G.; Wijga, P. W. O. Recl. Trav. Chim. Pays-Bas 1946, 65, 556-563. (39) Bey, O.; Eigenberger, G. Chem. Eng. Sci. 1997, 52, 1365-1376. (40) Giese, M.; Rottscha¨fer, K.; Vortmeyer, D. AIChE J. 1998, 44, 484-491. (41) Tallarek, U.; Scheenen, T. W. J.; Van As, H. J. Phys. Chem. B 2001, 105, 8591-8599. (42) Coelho, D.; Shapiro, M.; Thovert, J.-F.; Adler, P. M. J. Colloid Interface Sci. 1996, 181, 169-190. (43) Patankar, N. A.; Hu, H. H. Anal. Chem. 1998, 70, 1870-1881. (44) Arulanandam, S.; Li, D. Colloids Surf., A 2000, 161, 89-102. (45) Ren, L.; Li, D. J. Colloid Interface Sci. 2001, 243, 255-261. (46) MacInnes, J. M. Chem. Eng. Sci. 2002, 57, 4539-4558. (47) Dutta, P.; Beskok, A.; Warburton, T. C. J. Microelectromech. Syst. 2002, 11, 36-44. (48) Erickson, D.; Li, D. J. Phys. Chem. B 2003, 107, 12212-12220. (49) Fu, L.-M.; Lin, J.-Y.; Yang, R.-J. J. Colloid Interface Sci. 2003, 258, 266-275. (50) Li, B. M.; Kwok, D. Y. J. Chem. Phys. 2004, 120, 947-953.

Electroosmotic Flow in Porous Media

al.42 developed a numerical solution to the problem of EOF through porous media, in particular, a random packing of impermeable, nonconducting spherical particles. The authors used a relatively coarse spatial discretization step (8.54 grid points per sphere diameter) to reach reasonable computation times. Since their model could provide a good accuracy only when the discretization step became less than the EDL thickness, it restricts the application to the case of a thick EDL with respect to the sphere size. In addition, the Debye-Hu¨ckel (low ζ-potential) approximation and the assumption of small perturbations from equilibrium have been introduced to that model. Recently, we have developed a numerical scheme for simulating the EOF in spatially regular sphere systems.51 The spatial periodicity allowed to reduce the computational domain to a unit cell and employ fine computational grids (with up to 200 grid points per sphere diameter). This numerical model was used to investigate the dependence of EOF through a simple cubic (SC) array of hard spheres on the applied electrical field strength, value of the (uniform) ζ-potential, and ionic strength of the electrolyte solution. It was found that if the absolute value of the ζ-potential is higher than 50 mV and the aspect ratio ds/λD of the sphere diameter to the EDL thickness (approximated by the Debye length) exceeds 50, a high computational accuracy can be achieved only when the spatial discretization step is less than 0.01ds. Naturally, the use of very fine grids results in a substantial increase of computation times. One possibility for avoiding this problem is to assume that the EDL thickness is negligibly small compared to the sphere diameter (TDL limit with ds/λD f ∞) which allows decoupling the hydrodynamic and electrostatic problems. The purpose of this contribution is to present a numerical model and results of simulations for EOF through dense regular and random arrays of hard (nondeformable, impermeable, and nonconducting) spheres. Important aspects of pore space morphology (sphere size distribution), surface heterogeneity (mismatch in electrokinetic potentials at the confining wall and sphere surface), and fluid phase properties (electrical double layer thickness) are investigated concerning their influence on the EOF dynamics over microscopic and macroscopic spatial domains. This article is organized as follows. In section 2 we formulate the governing equations along with corresponding boundary conditions for the general case and the TDL approximation. In section 3 we outline briefly the employed numerical methods, in particular, the packing algorithm, general computational scheme, and lattice-Boltzmann approach to the simulation of the hydrodynamics in the EOF problem. In section 4 we present the results of simulations addressing the influence of the velocity boundary condition at the solid-liquid interface (slip vs no-slip), aspect ratio ds/λD (relative sphere diameter), global and local order in a sphere packing (regular vs random and confined vs bulk), and the distribution of ζ-potential and sphere size on the EOF dynamics over the involved length scales. 2. Theoretical Background A. General Description. The velocity field of an incompressible Newtonian fluid is governed by the Navier-Stokes equation. The steady flow velocity v of the electrolyte solution representing a divergence-free velocity field (∇‚v ) 0) is given by (51) Hlushkou, D.; Apanasovich, V.; Seidel-Morgenstern, A.; Tallarek, U. Chem. Eng. Commun., accepted for publication.

Langmuir, Vol. 21, No. 13, 2005 6099

Fv‚∇v ) -∇p + µ∇2v - Fe∇ψ + Fa

(1)

where p is hydrostatic pressure, µ and F are the dynamic viscosity and mass density of the fluid, respectively, Fe is charge density, ψ local electrical potential, and a denotes the acceleration due to gravitational forces. The local electrical potential is governed by the Poisson equation

Fe

2

∇ ψ)-

∑i nizi

qe )-

0r

(2)

0r

where 0 and r are the vacuum permittivity and dielectric constant of the fluid, qe is elementary charge, and ni and zi denote the number concentration and algebraic valence of the ith ionic species in the electrolyte solution, respectively. In turn, individual number flux ji is determined by the Nernst-Planck equation

ji ) -Di∇ni - Di

qezini ∇ψ + niv kBT

(3)

where Di is the diffusion coefficient, kB is the Boltzmann constant, and T is the temperature of the fluid. Assuming absence of chemical reactions number species concentrations obey the conservation principle expressed in the form of a continuity equation, ∇ji ) 0, which combines with eq 3 to yield

(

)

qezini ∇ψ + niv ) 0 ∇‚ -Di∇ni - Di kBT

(4)

Thus, for the determination of the velocity field the relevant set of coupled partial differential equations (eqs 1, 2, and 4) has to be solved subjected to the corresponding boundary conditions formulated at the solid-liquid interface Γ (inner boundary conditions) and at the edges Γout of the simulated spatial domain (outer boundary conditions). B. Inner Boundary Conditions. Numerical simulations presented in this paper were carried out with the zero-normal-flux boundary condition at the solid-liquid interface assuming hard (impermeable) spheres and no species flux due to surface reactions, adsorption, or dissociation of surface groups. The velocity at the solidliquid interface was either set as zero (no-slip velocity boundary condition), v|Γ ) 0, or defined by the HelmholtzSmoluchowski equation as a slip velocity10

v|Γ ) -

0rζE ) vslip µ

(5)

where E refers to the local strength of the applied electrical field. The latter boundary condition was realized in simulations involving the TDL approximation described below. Finally, the imposed electrical boundary condition was either of the Neumann or Dirichlet type

n‚∇ψ|Γ )

σ or ψ|Γ ) ζ + φ 0r

(6a,b)

where n is the local outer normal to Γ and σ is the surface charge density. The local electrical potential ψ at the solidliquid interface (or, more properly, at the shear plane) is assumed to be composed of two independent components, i.e., the ζ-potential and potential φ due to the externally applied electrical field (Eext).

6100

Langmuir, Vol. 21, No. 13, 2005

Hlushkou et al.

C. Outer Boundary Conditions. All sphere packings in this paper are assumed as spatially periodic. While the spatial periodicity of a regular sphere array such as the SC packing is an inherent property, spatial periodicity with respect to random sphere packings needs clarification. In our simulations the representative unit for random sphere arrays contained a large number of spheres, up to more than 1000. Realization of periodic boundaries assumes that the position of a sphere on one side of the packing (within the representative unit) influences the position of spheres at the opposite side. As a result, space is filled regularly at the macroscale, while reproducing the representative unit, but randomness prevails locally at the microscale, within this unit. Periodic boundaries are often encountered in modeling of flow and hydrodynamic dispersion in microscopically disordered, macroscopically homogeneous media.52,53 The application of periodic boundary conditions assumes that values of a generic (vector or scalar) physical field Θ like flow velocity, pressure, species concentration, or electrical potential at the edges of a representative unit are given by [|Θ|] ) const, where [|...|] is the difference between local values of that quantity at opposite points lying on the corresponding edges of the unit. It should be pointed out that periodic boundary conditions cannot be applied to turbulent flows, while for most electrokinetic systems the characteristic length scales (∼10-6-10-4 m), flow velocities (∼10-3 m/s), and viscosities (∼10-6 m2/s) are such that the Reynolds numbers are small (Re ∼ 10-3-10-1) and flows remain laminar. D. Thin Double Layer (TDL) Approximation. Although the analytical or numerical solution of the EOF problem in systems with complex pore space morphology is a difficult (or unresolvable) task, many systems meet requirements for application of the TDL approximation. The Debye length is given by

λD )

(

0rkBT

qe

2

)

∑ zi ni,∞ 2

1/2

(7)

Thus, λD is less than 100 nm in a 1:1 electrolyte solution with molar concentrations higher than 10-5 M. With the TDL approximation the shear plane (characterized by ζ) and the slipping plane (characterized by vslip) actually coincide, and the no-slip condition at the solid-liquid interface is replaced by the velocity boundary condition formulated at the slipping plane according to eq 5. Assuming that the packing contains spheres with large relative diameter, this replacement of the velocity boundary condition causes only insignificant error since the distance between slipping plane and solid-liquid interface is negligibly small compared to pore space dimensions between the spheres.26 In this picture, EOF slips past the spheres surface with a velocity given by eq 5. Thus, the TDL approximation allows rewriting the original problem (eqs 1, 2, and 4) in the form of the Stokes equation

µ∇2v ) ∇p

(8)

which is solved subjected to the slip-velocity boundary condition. Since the solution of the EOF problem is sought only outside the EDL, i.e., in the electroneutral solution, the applied electrical field is governed by the Laplace equation (52) Lowe, C. P.; Frenkel, D. Phys. Rev. Lett. 1996, 77, 4552-4555. (53) Maier, R. S.; Kroll, D. M.; Bernard, R. S.; Howington, S. E.; Peters, J. F.; Davis, H. T. Phys. Fluids 2000, 12, 2065-2079.

∇E ) ∇2ψ ) 0

(9)

Assuming nonzero fluid conductivity the above equation is solved subjected to the insulation-type boundary condition (eq 6a) with σ ) 0. 3. Numerical Methods A. Generation of Random Sphere Packings. Prior to the description of the approach followed in resolving numerically the EOF problem, we briefly discuss some aspects of modeling a dense packing of spherical particles. Random-close hard sphere packings have been generated using an algorithm based on the technique reported by Jodrey and Tory54 (swelling-up algorithm) which has been further developed to allow for an arbitrary shape of the confining container and distribution of the sphere size, as well as the possibility to continue the generation until a predefined packing porosity is realized. More details of the employed packing algorithm can be found elsewhere.55 In the present work a collection of 10 packings with mean interparticle porosity of inter ) 0.4 but different macroscopic shape (confined vs bulk packing) and size distribution of the hard spheres was employed for analysis (see Table 1). Packings 1-9 were generated in a cylindrical container (confined packings) with periodic boundaries in axial direction, while packing 10 was generated in a cubic container using periodic boundaries in all dimensions (unconfined or bulk packing). Packings 1-3 consist of monosized spheres, but they reveal slight differences in morphology due to different initial seeds used for examining the statistical reproducibility of the simulated results. Figure 1 demonstrates that despite identical macroscopic properties (such as the mean porosity, geometrical dimensions, and sphere size distribution), packings 1-3 show slightly different radial porosity distributions indicating differences in microstructure. At the same time, the distribution functions have common features, in particular, a similar location and magnitude of maxima and minima, as well as a similar decay of oscillations at increasing distance from the wall. It means that these packings have the same macroscopic geometrical properties and radial variations in porosity related to a more ordered (layered) structure near the confining wall. Packing 10 (which also consists of monosized spheres, but without a confining container) was used to analyze this effect resulting in the distortion of packing randomness close to a rigid wall. Radial distribution functions of interparticle porosity in the computer-generated randomclose hard sphere packings are very similar to those found experimentally in real packings of spherical particles56,57 (Figure 1) which indicates that the simulated packings represent adequately the real ones. Packings 4-7 consist of hard spheres with a Gaussian size distribution. They all have the same mean diameter, but standard deviations of 0.1, 0.2, 0.3, and 0.4, respectively (Table 1). Packings 8 and 9 contain bisized spheres with diameters ds and 0.5ds. Ratios between number fractions of the smaller and larger spheres are 1:1 and 1:8 for packings 8 and 9, respectively, which translates to ratios between volume fractions of the smaller and larger (54) Jodrey, W. S.; Tory, E. M. Phys. Rev. A 1985, 32, 2347-2351. (55) Hlushkou, D. Numerische Simulation von Stro¨mung und Massentransport in (elektro-) chromatographischen Systemen. Ph.D. Thesis, Otto-von-Guericke-Universita¨t Magdeburg, Germany, 2004. (http:// diglib/uni-magdeburg.de/Dissertationen/2004/dzmhlushkou.htm). (56) Benenati, R. F.; Brosilow, C. B. AIChE J. 1962, 8, 359-361. (57) Giese, M. Untersuchung der Stro¨mung in poro¨sen Medien unter Beru¨cksichtigung effektiver Viskosita¨ten. Ph.D. Thesis, TU Mu¨nchen, Germany, 1997.

Electroosmotic Flow in Porous Media

Langmuir, Vol. 21, No. 13, 2005 6101

Table 1. Employed Sphere Packings and Average Interparticle EOF Velocitiesa packing no.

distribution of sphere size

packing geometry

vinter (mm/s)

1 2 3 4 5 6 7 8 9 10

monosized (σ ) 0) monosized (σ ) 0) monosized (σ ) 0) Gaussian (σ ) 0.1) Gaussian (σ ) 0.2) Gaussian (σ ) 0.3) Gaussian (σ ) 0.4) bisized (d1 ) 2d2) bisized (d1 ) 2d2) monosized (σ ) 0)

cylindrical, confined cylindrical, confined cylindrical, confined cylindrical, confined cylindrical, confined cylindrical, confined cylindrical, confined cylindrical, confined cylindrical, confined cubic, bulk

1.314 1.317 1.315 1.342 1.337 1.322 1.324 1.303 1.282 1.175

a All sphere packings have the same mean interparticle porosity inter ) 0.4. Simulations were carried out with the slip-velocity boundary condition at the solid-liquid interface using Eext ) 50 kV/m and ζ ) ζw ) ζp ) -50 mV (packings 1-9). The liquid is a strong 1:1 electrolyte at T ) 298 K characterized by r ) 80, µ ) 8.9 × 10-4 kg m-1 s-1, and F ) 1 × 103 kg m-3.

Figure 1. Radial distribution of interparticle porosity against the normalized distance from the container inner surface for random-close hard sphere packings 1-3. Symbols correspond to the experimental data from Benenati and Brosilow56 (circles) and Giese57 (squares).

spheres of 1:8 and 1:1. The cylindrical packings 1-9 have a diameter of 10ds and length 5ds, while packing 10 has the dimension 10ds × 10ds × 10ds. For illustration, sphere packings 1, 4, 8, and 10 are shown in Figure 2. B. General Computational Scheme. The next step consisted of space discretization which was performed by introducing a uniform spatial grid. Each elementary cell of the grid was associated with either a solid (S-cell) or liquid phase (L-cell) depending on the location of the cell center. The set of coupled Navier-Stokes, Poisson, and Nernst-Planck equations subjected to the appropriate boundary conditions was solved following an iterative numerical scheme. At each iteration the three governing equations were solved by turns using the same uniform computational grid. During computations the electrical potential, species number concentration, and fluid velocity were determined at the centers of only the L-cells. If an L-cell was adjacent to an S-cell, the corresponding boundary conditions were involved to determine values of the physical quantities. In the case of the TDL approximation the solution was reduced to a single iteration for the Laplace and Stokes equations subjected to the insulation-type and slip-velocity boundary conditions, respectively. Solutions for the Nernst-Planck and Poisson (Laplace) equations were obtained by conventional finite-difference methods (explicit second-order schemes), while the lattice-Boltzmann approach was applied for the

solution of the hydrodynamic problem. Simulations stopped when the divergence rate (relative difference between values calculated at two subsequent iterations) of the mean axial EOF velocity became less than a predefined value (10-5 in all simulations of this work). Computational grid resolutions providing sufficient numerical accuracy were estimated by examining the dependence of mean axial EOF velocity51 and packing hydraulic permeability58 on the grid resolution. Specifically, we used a spatial resolution of 40 grid nodes per sphere diameter (ds) for simulations involving the TDL approximation (considering them as particular hydraulic flow problem with slip velocity boundary condition at the solid surface) and 200 grid nodes per sphere diameter for simulations without the TDL approximation. Further technical details concerning the employed numerical methods and developed computational scheme can be found in a previous article59 which presents the results of EOF simulations in microfluidic channels. C. Lattice-Boltzmann Method. The EOF velocity field was obtained using the lattice-Boltzmann method (LBM) which is a simplified kinetic approach with discrete space and time.60-64 The macroscopic fluid dynamics is approximated by synchronous interactions between fictitious fluid particles on a regular lattice invoking the idea that the fluid flow is determined mainly by collective behavior of many molecules and not by detailed molecular interactions. It has recently been shown62 that the LBM corresponds to the truncation of the Boltzmann equation in a Hermite velocity spectrum space or, in other words, that it can be viewed as special finite-difference approximation of the Boltzmann equation. Among the advantages of the LBM are its inherent parallelism (in view of computational efficiency) and the capability to deal with geometrically complex solid-liquid interfaces such as those encountered in random sphere packings. The simplest, most popular form of the lattice-Boltzmann equation is the lattice-BGK (Bhatnagar-GrossKrook) equation.61 The starting point in its derivation is the continuum Boltzmann-BGK equation

1 ∂F + e‚∇F ) - (F - F eq) ∂t τ

(10)

τ is a relaxation parameter (in this work we assumed τ ) 1) and F ) F (r,e,t) is the one-particle probability distribution function so defined that F (r,e,t)d3rd3e represents the number of fluid particles which, at time t, are located in a phase-space element d3rd3e about a particles coordinate r and velocity e. The equilibrium probability distribution function F eq is given by

F

eq

)

[

(e - v) F exp 3/2 2RT (2πRT)

]

2

(11)

where R is the gas constant. Macroscopic quantities such as the density F and momentum Fv are calculated as (58) Kandhai, D.; Tallarek, U.; Hlushkou, D.; Hoekstra, A. G.; Sloot, P. M. A.; Van As, H. Philos. Trans. R. Soc. London, Ser. A 2002, 360, 521-534. (59) Hlushkou, D.; Kandhai, D.; Tallarek, U. Int. J. Numer. Methods Fluids 2004, 46, 507-532. (60) Higuera, F. J.; Jime´nez, J. Europhys. Lett. 1989, 9, 663-668. (61) Chen, S.; Doolen, G. D. Annu. Rev. Fluid Mech. 1998, 30, 329364. (62) He, X.; Luo, L.-S. Phys. Rev. E 1997, 56, 6811-6817. (63) Qian, Y. H.; d’Humie`res, D.; Lallemand, P. Europhys. Lett. 1992, 17, 479-484. (64) Gou, Z.; Zheng, C.; Shi, B. Phys. Rev. E 2002, 65, art. no. 046308.

6102

Langmuir, Vol. 21, No. 13, 2005

Hlushkou et al.

Figure 2. Random-close hard sphere packings 1, 4, 8, and 10 (see also Table 1). For packing 10 the periodic image of spheres which intersect the face boundaries is added.

hydrodynamic moments of the distribution function F

F)

∫F de

and Fv )

∫eF de

(12a,b)

Assuming a small macroscopic velocity with respect to the speed of sound (low-Mach number approximation) F eq can be expanded and written in the following form up to the v2 term

F eq ≈

( )[

]

(e‚v)2 (e‚v) F e2 v2 + exp 1 + 3/2 2 2RT RT 2RT (2πRT) 2(RT) (13) The continuous particle velocity-space e can be reduced to a finite set of discrete velocities eR, and the hydrodynamic moments of F , as well as its fluxes, can be preserved exactly because the moment integrals (eqs 12a,b) can be replaced by quadratures exactly62 up to a certain order in e. The coherent discretization with step ∆rR ) eR∆t allows rewriting eq 10 in the form of the lattice-BoltzmannBGK equation

F R(r + ∆rR, t + ∆t) ) 1 F R(r,t) - [F R(r,t) - F Req(r,t)] (14) τ where F R is the distribution function of the Rth discrete

velocity and F Req can be expressed as

F Req ≈

wRF

× (2πRT)3/2 eR2 (eR‚v) (eR‚v)2 v2 + 1+ (15) exp 2RT RT 2(RT)2 2RT

( )[

]

The weight constants wR depend on lattice geometry (we employed the D3Q19 model63) and can be obtained by approximating the moment integrals (eqs 12a,b) with a third-order Hermite formula.62 Thus, eq 14 is iteratively solved in two steps.

collision: 1 F ˜ R(r,t) ) F R(r,t) - [F R(r,t) - F Req(r,t)] (16) τ streaming: ˜ R(r,t) F R(r + ∆rR, t + ∆t) ) F

(17)

Finally, the macroscopic quantities are evaluated by

F(r,t) )

∑R F R(r,t)

(18a)

and

v(r,t) )

∑ReRF R(r,t) F(r,t)

(18b)

Electroosmotic Flow in Porous Media

Langmuir, Vol. 21, No. 13, 2005 6103

Figure 3. Dependence of EOF velocity through the SC array of hard spheres on applied electrical field strength for different values of relative sphere diameter. The fluid phase is a strong 1:1 electrolyte at T ) 298 K characterized by r ) 80, µ ) 8.9 × 10-4 kg m-1 s-1, and F ) 1 × 103 kg m-3. Simulations were carried out using the no-slip velocity boundary condition at the solid-liquid interface (ζ ) -50 mV).

The body force associated with electrostatic interaction between the applied electrical field and net charge in the liquid was incorporated in the above scheme by adding a force term in the lattice-Boltzmann-BGK equation.64 To satisfy the no-slip velocity boundary condition at the solidliquid interface, we employed the conventional bounceback scheme where momentum from an incoming fluid particle is bounced back in the opposite direction as it hits the wall. In turn, imposition of the slip-velocity boundary condition assumes that the fluid velocity at nodes adjacent to the solid-liquid interface for each iteration must be determined by eq 5 (instead of eq 18b). It should be mentioned that while the LBM approach has been widely used for the simulation of hydraulic flow through porous materials, an extension of this formalism to the electrokinetic microfluidics in dense sphere packings until now has been an open problem. 4. Results and Discussion A. Simple Cubic (SC) Array of Spheres. We first simulated EOF through an SC array of hard, contacting spheres. The interparticle porosity of this system is inter ) (1 - π/6) ) 0.476, and its spatial periodicity allowed reduction of the simulated domain to a cube circumscribing a single sphere as representative unit. This allowed performing the simulations in much shorter computation times or, where needed, with an increased accuracy by using finer computational grids. We introduce the average interparticle axial EOF velocity in a sphere packing as

vinter )

1 Vinter

∫vx dVinter

(19)

where Vinter denotes the interparticle void space and vx is the local axial velocity. Axial refers to a vector component parallel to the applied electrical field. The volumeaveraged axial velocity vvol (superficial or open-tubular velocity) for a packing of impermeable spheres is given by vvol ) intervinter and is less than the average axial interparticle flow velocity (inter < 1). We begin the simulation of EOF through an SC array of hard spheres by imposing the no-slip velocity boundary condition at the solid-liquid interface. Figure 3 shows the computed values of vvol for selected relative sphere

Figure 4. Dependence of electroosmotic mobility (eq 20) in the SC array of hard spheres on relative sphere diameter (cf. Figure 3). Simulations were carried out with the no-slip velocity boundary condition at the solid-liquid interface, ζ ) -50 mV and Eext ) 50 kV/m.

diameters (ds/λD) in dependence of the externally applied electrical field strength (Eext). Under the given set of conditions the linear dependence observed in Figure 3 can be characterized by a constant electroosmotic mobility β defined as

β)

vvol Eext

(20)

which we use, in turn, to analyze the EOF dynamics as a function of sphere diameter, EDL thickness, and the fluid phase properties. To validate the numerical approach described in section 3 we compared calculated values of β with those reported by Coelho et al.42 These authors presented a model for EOF in sphere packings assuming a low ζ-potential, weak electrical field, and small applied pressure gradient. Further, they showed results only for small relative sphere diameters (ds/λD e 22.72) due to limitations in the computational resources. Figure 4 demonstrates a good agreement between the results obtained by either approach regardless of the underlying numerical methods, implying that both describe adequately the electrokinetic microfluidics in an SC array of impermeable, nonconducting spheres. In addition, our approach allows extention of the work of Coelho et al.42 to (i) higher electrical field strengths because it is not restricted to the smallperturbation limit, (ii) higher ζ-potentials as it does not involve the Debye-Hu¨ckel approximation, and (iii) larger relative sphere diameters due to higher computational efficiency. Concerning the last aspect we analyzed the dependence of electroosmotic mobility on the relative sphere diameter for 2.272 e ds/λD e 100 (Figure 4). Lower values of β at smaller relative sphere diameters (ds/λD < 50, Figure 4) can be explained by the fact that, as the EDL thickness approaches the sphere size, the EOF velocity locally cannot realize the maximum possible value due to noticeable EDL overlap. It first becomes pronounced close to the contact points of the neighboring spheres, but for ds/λD f 0 it will cover the whole interparticle pore space. This effect, which also reduces the average EOF velocity through the sphere packing, is less important for a thinner EDL or larger sphere diameter (ds/λD g 50) when the EDL overlap essentially disappears, as indicated by our simulations for ds/λD ) 50 and 100 (Figure 4). Finally, for ds/λD f ∞ the distance over which the fluid velocity rises from zero

6104

Langmuir, Vol. 21, No. 13, 2005

at the surface to a limiting value beyond the EDL has become infinitesimally short with respect to the EDL thickness. Then, the no-slip velocity boundary condition at the solid-liquid interface may be replaced using the slip velocity determined by eq 5. Results in Figure 4 demonstrate that the mean EOF velocity at increasing ds/λD approaches asymptotically a value obtained by imposition of the slip-velocity boundary condition (dashed line). Naturally, the study of EOF without the TDL approximation in sphere packings at large relative sphere diameter requires great computational efforts. It is caused by the need to retain a fine computational grid relative to the EDL thickness which, at the same time, has to cover comparatively large system dimensions. With the slip-velocity boundary condition it is possible to significantly reduce computation times. For instance, in our simulations they could be reduced by a factor of the order of 103. To evaluate the actual error introduced by the TDL approximation when simulating EOF in the SC array at a large, but not infinite, relative sphere diameter, we analyzed velocity fields for ds/λD ) 100 obtained by modeling with either the no-slip or the slip-velocity boundary condition. The former simulation was run with a computational grid for the representative unit of 200 × 200 × 200 and the total computation time was about 6 × 105 s, while the latter was performed with a 40 × 40 × 40 grid and took about 600 s. Other conditions, including Eext ) 50 kV/m and ζ ) -50 mV as well as fluid phase properties, were identical. The simulated superficial flow velocities (0.567 and 0.648 mm/s for no-slip and slipvelocity boundary conditions) differed by 12.5%. In comparison, the velocities obtained in both cases with a 40 × 40 × 40 grid already differed by ca. 50% because the grid resolution becomes too low compared to the EDL thickness. A difference of 12.5% still may be caused by remaining EDL overlap, even at ds/λD ) 100, which reduces the velocity in the simulation with the no-slip boundary condition. In this context, the reduction of EOF velocity due to EDL overlap in an open capillary with diameter dc corresponding to the smallest size of a throat in the SC array of spheres (dc ≈ 0.414ds and dc/λD ≈ 20) estimated with the model of Rice and Whitehead16 is about 6.5%. In addition, the results presented in Figure 4 indicate that the superficial EOF velocity at ds/λD ) 100 has not yet reached the asymptotic value corresponding to ds/λD f ∞. We were not able to perform simulations with the no-slip velocity boundary condition for ds/λD > 100 due to restrictions in computational resources. The trend in Figure 4 allows the anticipation that by a further increase of the relative sphere diameter a closer, if not quantitative, agreement with the TDL approximation is achieved. For a complementary analysis of the TDL approach, we compared probability functions (histograms) of the local axial interparticle velocity (vx) obtained from simulations with the apparent slip and no-slip velocity boundary conditions, i.e., in an SC sphere array for ds/λD f ∞ and ds/λD ) 100, respectively. Figure 5 indicates a structural similarity of both simulated velocity fields. (To facilitate visual comparison of the histograms, we generally provide curves generated by interconnecting center values of histogram bars.) Although the histogram related to the slip-velocity boundary condition is statistically noisier due to a smaller sample volume (inter × 403 ≈ 3.0 × 104 elements vs inter × 2003 ≈ 3.8 × 106 elements for the other one), it is obvious that both histograms display similar multipeak shapes. Distinctions in peak location and global velocity distribution can originate, in addition to the differences in velocity boundary condition, also in a

Hlushkou et al.

Figure 5. Probability distributions of axial interparticle EOF velocity in an SC array of hard spheres obtained with the slip (ds/λD f ∞) and no-slip (ds/λD ) 100) velocity boundary conditions. Fluid phase properties as in Figure 3, ζ ) -50 mV and Eext ) 50 kV/m.

remaining inequality of relative sphere diameters (ds/λD f ∞ vs ds/λD ) 100). At this point it should be mentioned that multiple peaks in the probability distribution for the local axial velocity were observed as well in the numerical simulation of hydraulic flow through an SC packing of hard spheres.65 A possible explanation for this feature can be found in the geometrical periodicity of this structure. For verification, we detected the location of computational grid cells corresponding to three selected peaks in the probability distribution (peaks 1, 2, and 3 in Figure 5), as well as to regions with a negative axial velocity. Figure 6 shows the three-dimensional distribution of computational grid cells within a representative unit of the SC array of spheres, along with projections onto the coordinate planes. As expected, the distributions are azimuthally symmetric with respect to the axial (x-)direction, but they show pronounced differences in the up- and downstream halves of the representative cube. This analysis indicates that the selected axial flow velocities and the range of negative velocities occur in specific regions of the pore space. For instance, the highest velocity (peak 3, Figure 5) dominates in throats formed by every four touching spheres in the yz-plane, i.e., perpendicular to the axial direction (Figure 6c). The region close to the contact points of these spheres is characterized by a stronger local electrical field strength due to the reduced cross-sectional area of the pore space which, in turn, results in a higher local EOF velocity. The grid cells corresponding to vx ≈ 2.6 mm/s (peak 3, Figure 5) are concentrated around these contact points. By contrast, the next lower velocity (peak 2, vx ≈ 1.3 mm/s) corresponds to regions distributed symmetrically around straight lines passing in the axial direction through the centers of the above-mentioned throats. These lines can be considered as the axes of open channels formed between spheres of the SC array. Finally, the lowest axial interparticle velocity (peak 1 in Figure 5, vx ≈ 0.64 mm/s) mainly corresponds to regions around the contact points of spheres neighbored in the axial direction. The higher probability compared to the other two velocities can be explained by the relatively large total volume of fluid regions behind peak 1. Compared to grid cells associated with the three selected peaks of the velocity probability distribution (Figure 5), fluid regions with a negative axial flow velocity are located (65) Schure, R. M.; Maier, R. S.; Kroll, D. M.; Davis, H. T. J. Chromatogr., A 2004, 1031, 79-86.

Electroosmotic Flow in Porous Media

Langmuir, Vol. 21, No. 13, 2005 6105

Figure 6. Three-dimensional distributions and projections onto Cartesian planes of computational grid cells containing axial interparticle flow velocities corresponding to (a) peak 1 (vx ≈ 0.64 mm/s), (b) peak 2 (vx ≈ 1.3 mm/s), and (c) peak 3 (vx ≈ 2.6 mm/s), as well as (d) negative values in the velocity probability distribution obtained for the simulation with no-slip velocity boundary condition of EOF through the SC array of hard spheres (see Figure 5, dashed line). Distributions and their projections are based on 13417, 5877, 4167, and 30525 points (peaks 1, 2, and 3 and negative velocities, respectively) which denote the corresponding cell centers. The external electrical field (Eext ) 50 kV/m) is applied in the positive x-direction.

in a thin shell close to the solid-liquid interface (Figure 6d). In this respect it should be noted that the numerical simulation of hydraulic flow through regular arrays of hard spheres (simple, body-centered, or face-centered cubic)65 has also detected regions with negative velocities. An analysis of their location in the body-centered cubic

array showed that they originate in vortices which develop around contact points of the particles.66 Figure 6d reveals that grid cells corresponding to negative axial velocities are mainly concentrated in cusp regions near the up- and downstream poles of the spheres. As with the hydraulic flow they can be related to recirculating flow behind a

6106

Langmuir, Vol. 21, No. 13, 2005

Hlushkou et al.

sphere. At the same time, a small population of negative velocities in grid cells near the equator of a sphere may also indicate the presence of spurious mass and momentum fluxes caused by the discrete approximation of continuous-space gradient operators. Further simulations using lattice-space gradient operators which avoid the occurrence of these spurious fluxes at thermodynamic equilibrium will help to clarify this aspect.67 It should be realized that the mean value of the simulated negative axial velocities is much smaller than that for positive velocities irrespective of whether the former originate in spurious fluxes, real recirculation cells, or a combination of both. By contrast, the numerical simulation with the slip-velocity boundary condition did not reveal negative axial velocities. It is in agreement with the TDL formalism assuming similitude for the irrotational electrical field and EOF velocity field under the given set of conditions.26 During further validation we compared results for simulation of EOF in the SC array of hard spheres (with the slip-velocity boundary condition) to those obtained by the macroscopic model of Overbeek and Wijga.38 According to their work which also involves the TDL approximation the volume-averaged EOF velocity through a packing of impermeable, nonconducting particles with uniform distribution of the ζ-potential is calculated by

vvol ) -

( )

0rζEext K* µ K∞

(21)

where K∞ is the conductivity of equilibrium electrolyte beyond the EDL and K* is the conductivity of the packing saturated with the same liquid. The conductivity ratio K*/K∞ for the SC array can be determined with an expression developed by Johnson and Sen68

K* )1K∞

3(1 - inter)

[

2 + (1 - inter) - 12

]

3(1 - inter) 4π

10/3

(U40)2 (22)

where the average interparticle porosity for the SC array of touching spheres is inter ) 1 - π/6. U40 is a lattice sum and about 3.1083 for the SC sphere array.69 Equations 20 and 21 allow expression of the electroosmotic mobility in the following form

β)-

( )

0rζ K* µ K∞

(23)

With eqs 22 and 23 we obtain β ) 1.383 × 10-8 m2 V-1 s-1 (using r ) 80, µ ) 8.9 × 10-4 kg m-1 s-1, and ζ ) -50 mV), while the corresponding value for β obtained in our simulation with slip-velocity boundary condition using eq 20 is 1.300 × 10-8 m2 V-1 s-1, differing from the analytical result by ca. 6.0%. This good agreement together with the earlier results indicates that EOF in hard sphere packings at a large relative sphere diameter can be reproduced by simulations which employ the TDL approximation and slip-velocity boundary condition at the solid-liquid interface. B. Random-Close Sphere Packings. As the following step we simulated the EOF in random-close packings of (66) Maier, R. S.; Kroll, D. M.; Kutsovsky, Y. E.; Davis, H. T.; Bernard, R. S. Phys. Fluids 1998, 10, 60-74. (67) Capuani, F.; Pagonabarraga, I.; Frenkel, D. J. Chem. Phys. 2004, 121, 973-986. (68) Johnson, D. L.; Sen, P. N. Phys. Rev. B 1988, 37, 3502-3510. (69) Venema, P.; Struis, R. P. W. J.; Leyte, J. C.; Bedeaux, D. J. Colloid Interface Sci. 1991, 141, 360-373.

Figure 7. Probability distributions of the axial EOF velocity in random-close hard sphere packings 1-3. Simulations employed the slip-velocity boundary condition. Fluid phase properties as in Figure 3, ζw ) ζp ) -50 mV and Eext ) 50 kV/m.

hard spheres confined by a cylindrical container. Arrays of monosized spheres were studied first (packings 1-3, see Table 1). To justify use of the TDL approximation, we selected a relative diameter ds/λD ) 1000 which corresponds, e.g., to a sphere diameter of 10 µm and Debye screening length of 10 nm. We compared the simulated volume-averaged axial EOF velocity to that obtained by the model of Overbeek and Wijga38 which assumes identical values of ζ-potential at the cylinder wall (ζw) and the particle surface (ζp). Further, we used K*/K∞ ) 0.26 ( 0.02, a value that has been obtained experimentally for a fixed bed of particles with an interparticle porosity and column-to-particle diameter ratio of 0.4 and 10.7, respectively.70 The electroosmotic mobility calculated by eq 23 for an aqueous phase at 298 K is β ) (1.035 ( 0.0796) × 10-8 m2 V-1 s-1 using r ) 80, µ ) 8.9 × 10-4 kg m-1 s-1, and ζ ) ζw ) ζp ) -50 mV. Under the actual conditions simulations of EOF in confined packings of monosized spheres with different microstructure (cf. Figure 1) resulted in close values of β of 1.051, 1.054, and 1.052 × 10-8 m2 V-1 s-1 for packings 1, 2, and 3, respectively. These values differ from the above analytical one by just 1.55, 1.84, and 1.64%. In addition, probability distributions for interparticle axial velocities (Figure 7) are also very similar, which allows us to conclude that differences in microstructure of packings 1-3 caused by different initial seeds do not noticeably affect the macroscopic EOF behavior. In contrast to EOF in the SC array (Figure 5) the velocity probability distributions for random-close sphere packings reveal a unimodal character (Figure 7). The most probable axial interparticle EOF velocities are slightly smaller than mean values due to the skewn shape of these distributions toward higher velocities. The main conclusion to be drawn from the analysis of EOF velocity fields in confined random sphere packings is that the local axial velocity is nonuniform (Figure 7). At first, this finding may be surprising as it is in contrast to the ideal pluglike velocity profile (anticipated for the TDL-limit along with a δ-like velocity probability distribution) in a single straight, open, i.e., unpacked capillary with uniform surface properties.19,20 On the other hand, it does not contradict the similitude for electrical and EOF velocity fields under the given set of conditions, as we will discuss below in more detail. Concerning a general nonuniformity of local axial EOF velocities, it may be expected that the radial velocity distribution in confined (70) Vallano, P. T.; Remcho, V. T. J. Phys. Chem. B 2001, 105, 32233228.

Electroosmotic Flow in Porous Media

Langmuir, Vol. 21, No. 13, 2005 6107

Figure 8. Radial distribution of the axial EOF velocity in sphere packing 1. Simulations employed the slip-velocity boundary condition. Fluid phase properties as in Figure 3, ζw ) ζp ) -50 mV and Eext ) 50 kV/m.

sphere packings bears a close relation with the radial porosity distribution, similar to hydraulic flow.71 It should be realized that, apart from the variations in hydrodynamic resistance, local differences in porosity due to a varying packing order close to the confining wall influence the local electrical field strength and, in turn, the local slip-velocity of EOF at the solid-liquid interface in view of eq 5. To correlate explicitly radial variations in axial EOF velocity with the packing porosity, we modified eq 21 which includes only mean values of the EOF velocity and conductivity ratio. The conductivity of a fixed bed of nonconducting spheres (saturated with electrolyte solution) relative to that of the bulk solution results from the72 (i) reduction of void space and total cross-sectional area, (ii) decrease of local electrical field strength and increase in migration distance for ions due to the tortuous nature of the interparticle channels, and (iii) dilatation and constriction of the channels. In practice, the following empirical relationship for conductivity ratio and interparticle porosity is often used73

K* ) inter1.5 ∞ K

(24)

Hence, eq 21 can be rewritten as

vvol(r) ) -

0rζEext inter1.5(r) µ

(25)

where vvol(r) denotes the volume-averaged axial EOF velocity in a thin annulus of radius r. In Figure 8 we compare the simulated vvol(r) in packing 1 to that obtained with eq 25. These data allow the conclusion that radial variations of the axial EOF velocity in confined cylindrical sphere packings strongly correlate with the radial porosity distribution. Further, the good agreement between the results of our simulation and predictions based on eq 25 implies that this correlation can be well reproduced using the analytical result by employing the actual radial porosity distribution. Equation 25 is applicable, however, only to systems with a uniform ζ-potential. (71) Maier, R. S.; Kroll, D. M.; Bernard, R. S.; Howington, S. E.; Peters, J. F.; Davis, H. T. Phys. Fluids 2003, 15, 3795-3815. (72) Choudhary, G.; Horva´th, Cs. J. Chromatogr., A 1997, 781, 161183. (73) De la Rue, R. E.; Tobias, C. W. J. Electrochem. Soc. 1959, 106, 827-833.

Figure 9. Radial distributions of axial EOF velocity in sphere packings 1 and 4-7 (Table 1). Simulations employed the slipvelocity boundary condition. Fluid phase properties as in Figure 3, ζw ) ζp ) -50 mV and Eext ) 50 kV/m.

With respect to hydrodynamic dispersion, a central issue for efficiency in (electro)chromatographic separations of, e.g., bioanalytical or pharmaceutical target molecules, the analysis and optimization of the underlying interparticle flow field are most important. The radial distribution of interparticle axial EOF velocity vinter(r) in packing 1 (Figure 8) is similar to that of vvol(r) and demonstrates the same features, in particular, damped oscillations with highest amplitude near the wall. This observation suggests that the use of spheres with nonuniform size distribution may reduce the heterogeneity of interparticle axial EOF velocities in confined packings caused by the wall effect. Smaller spheres better fill voids between larger ones, especially between larger ones and the hard wall, thereby reducing the amplitude of local porosity fluctuations. As a consequence, the EOF velocity field should become more uniform on a macroscopic scale over which the wall effect prevails. To verify this hypothesis we simulated EOF in confined sphere packings with different Gaussian and discretely bimodal size distributions (packings 4-9, see Table 1). Values of the simulated electroosmotic mobility β for packings 4-9 range from 1.026 × 10-8 m2 V-1 s-1 (packing 9) to 1.0736 × 10-8 m2 V-1 s-1 (packing 7) being close to those for the packings of monosized spheres of about 1.05 × 10-8 m2 V-1 s-1. Radial distributions of interparticle flow velocity obtained for packings 4-7 (Gaussian size distributions) are shown in Figure 9, together with vinter(r) for packing 1 (monosized spheres). In the core region all velocity distributions oscillate stochastically with similar amplitudes. The analysis further indicates that, starting from monosized spheres (σ ) 0, packing 1), the increasing standard deviation of the Gaussian size distributions (packings 4-7) leads to a decrease of EOF heterogeneity. It is most pronounced for packings 6 and 7 with broader sphere size distributions (σ ) 0.3 and 0.4). Thus, variations in particle size can reduce the nonuniformity of local axial EOF velocities, particularly in the critical wall region of a confined sphere packing. This will decrease hydrodynamic dispersion which, in turn, improves the long-distance transport by EOF, e.g., during chromatographic separation of chemical and biological species in adsorbent media consisting of random-close sphere packings. The radial distribution of vinter in packings 8 and 9 (Figure 10) with bisized spheres allows analysis of the influence of discrete particle fractions on radial EOF heterogeneity. Packing 8 which contains equal numbers of small (d ) 0.5ds) and large (d ) ds) particles (the

6108

Langmuir, Vol. 21, No. 13, 2005

Figure 10. Radial distributions of the axial EOF velocity in sphere packings 1, 8, and 9 (cf. Table 1). Simulations employed the slip-velocity boundary condition. Fluid phase properties as in Figure 3, ζw ) ζp ) -50 mV and Eext ) 50 kV/m.

resulting relative volume fractions are 1:8) reveals locally increased axial flow velocities compared to packing 1 about a distance r/ds ) 0.5 from the wall, corresponding to one diameter of the small particles. This can be explained by the multiplicity of oscillation periods related to the individual radial porosity distribution functions for each particle size. Because the first minimum and maximum of vinter(r) are located at distances of, respectively, 0.5ds and ds from the wall (Figure 8), it results in an antiphase superposition of the two oscillation harmonics. The influence of discrete particle fractions on the overall velocity distribution depends on relative weight. In packing 8 the small particles dominate in the peripheric annular region extending inward from the wall by a distance of 0.5ds. Therefore, their contribution to the behavior of vinter(r) in this immediate wall region is relatively strong. Further inward, both kinds of particles are found in almost equal amounts and the presence of the small particles in the core region leads to a just slightly smaller oscillation amplitude of vinter(r) compared to packing 1 containing only large particles. In contrast, packing 9 is characterized by a higher number fraction of the small particles (8:1 as compared to packing 8 with number fraction of 1:1) which results in equal volume fractions of both kinds of particles. Then, the small sphere size (d ) 0.5ds) dominates the behavior of vinter(r), particularly its oscillation period, also in the core region of packing 9 as seen in Figure 10. The above results demonstrate that the average interparticle axial velocity and volumetric rate of EOF in random-close sphere packings hardly depend on sphere size distribution for a given electrical field strength and constant macroscopic parameters such as mean porosity and geometrical dimensions of the packing. At thesame time, the sphere size distribution can affect the homogeneity and, consequently, the hydrodynamic dispersion characteristics of the underlying EOF velocity field. It is an aspect relevant to most microscale and analytical electrical-field related separation techniques because of the generally nonuniform size distribution of fine materials employed for the preparation of fixed beds used, e.g., in CEC and chip electrochromatography.1,3,5,12 To examine further the effect of packing confinement on the EOF velocity distributions, we continued with the simulation of EOF in a bulk (unconfined) random-close packing of monosized spheres. For estimating the volumeaveraged EOF velocity according to the model of Overbeek and Wijga,38 we used K*/K∞ ) 0.24 determined experi-

Hlushkou et al.

mentally by Zeng et al.74 for a packed capillary with porosity inter ) 0.37 and a column-to-particle diameter ratio dc/dp ) 150. Because the wall effect manifests itself only in a near-wall annulus of ca. 4-5 particle diameters,56 this aspect ratio appears high enough to assume the wall effect as being macroscopically negligible. Calculation by eq 23 yields β ) 0.955 × 10-8 m2 V-1 s-1 using r ) 80, µ ) 8.9 × 10-4 kg m-1 s-1, and ζ ) ζp ) -50 mV. The value obtained in our simulation for the bulk packing 10 is β ) 0.940 × 10-8 m2 V-1 s-1 resulting in a relative difference of only 1.57%. The simulated electroosmotic mobility and average axial velocity in the bulk packing are lower than for the confined packings 1-3 (see Table 1) despite the same average porosity and other parameters such as ζ-potential, Eext, and fluid phase properties. A possible explanation is related to the fact that confined sphere packings have a more ordered structure in the near-wall region (wall annulus) than in the bulk (core) region. As mentioned earlier, the local structure near the wall results in a reduced tortuosity of the flow paths and, consequently, it causes a local increase in both the conductivity ratio and flow velocity. This explanation is confirmed by Figure 11 which illustrates the distribution of the interparticle axial velocity obtained by averaging along the whole bed for packings 1 (confined) and 10 (bulk). The light-yellow regions near the wall of the confined packing (cf. Figure 8) indicate higher local axial EOF velocities, which also explain the higher average velocity under otherwise identical conditions. In addition to EOF heterogeneity introduced on a macroscopic scale by the local morphology in the nearwall region of confined sphere packings, a microscopic (pore-level) nonuniformity of interparticle axial EOF velocity seems to prevail in the core region. Figures 12 and 13 allow comparison of the velocity profiles in nearwall and core regions of confined packing 1. It is obvious that a nonuniform distribution of axial EOF velocity also exists in the interparticle pore space of the core region which originates in a very similar nonuniformity of the local electrical field strength (Figure 13). This interesting behavior should be interpreted in the context of an expected similitude for the electrical field and EOF velocity field in the TDL limit. On the basis of these data in Figures 12 and 13 it can be concluded that even the imposition of the slip-velocity boundary condition at the solid-liquid interface does not result in a locally flat velocity profile in the interstitial pore space of a sphere packing. By contrast, this boundary condition leads to a pluglike velocity distribution in a single-straight, open, homogeneous channel. As indicated, the velocity profile in a sphere packing is caused by the variations in cross-sectional area of the interstitial channels formed by neighboring spheres. This pore space morphology results in a nonuniform distribution of the local axial component of applied electrical field (Figure 13, bottom row), in particular, at the solid-liquid interface which leads to local variations in slip velocity (eq 5) used to formulate the velocity boundary condition at this interface. It should be pointed out that this interrelation responsible for the pore-level velocity profile in a sphere packing is different from the mechanism related to variations in tangential velocity within the EDL which cannot become active in the present study due to the TDL approximation. As indicated by Figure 13 the similitude between irrotational electrical and EOF velocity fields, expected to (74) Zeng, S.; Chen, C.-H.; Mikkelsen, J. C., Jr.; Santiago, J. G. Sens. Actuators, A 2001, 79, 107-114.

Electroosmotic Flow in Porous Media

Langmuir, Vol. 21, No. 13, 2005 6109

Figure 11. Distribution of axial EOF velocity (averaged over the complete length of a sphere packing) obtained by imposition of the slip-velocity boundary condition: left, packing 1 (cylindrical-confined); right, packing 10 (cubic-unconfined). Fluid phase properties as in Figure 3, ζw ) ζp ) -50 mV (packing 1) and ζp ) -50 mV (packing 10), Eext ) 50 kV/m.

Figure 12. Column cross-sectional distribution of axial EOF velocity in sphere packing 1. Fluid phase properties as in Figure 3, ζw ) ζp ) -50 mV and Eext ) 50 kV/m.

hold under the given set of conditions (thin EDL, uniform ζ-potential and fluid phase properties, low Reynolds number, etc.26), is indeed observed for both computed fields. The deviation from the ideal pluglike EOF velocity profile suggests that hydrodynamic dispersion with EOF through a sphere packing differs significantly from that in a single, open channel even for high aspect ratio of the transverse channel (and interstitial pore space) dimensions to the local EDL thickness. As a consequence, the use of capillary models for the simulation of dispersion phenomena in EOF

through random porous media, in particular, sphere packings can result in essential discrepancy with respect to more advanced models which account for the microscopic nonuniformity of the velocity profiles in individual channels, as well as a macroscopic nonuniformity due to the packing confinement. So far we were concerned with EOF in confined packings characterized by a uniform distribution of ζ-potential assuming, in particular, identical values at the sphere surface (ζp) and the confining wall (ζw). However, fixed beds in chromatographic columns, catalytic reactors, or electroosmotic micropumps are characterized by not only a different chemical but also a different physical nature of the involved surfaces, including magnitude and even sign of the ζ-potentials. For instance, electrokinetic potentials at the inner surface of quartz capillaries up to -100 mV and above are not unusual, although for buffers and ionic strengths encountered in electrochromatographic applications this value may be regarded as an upper limit and ζw typically ranges between -50 and -100 mV.75-78 By contrast, ζp of many commercial, e.g., cationexchange particles, used in CEC is significantly lower than ζw of the fused-silica capillaries.79 Then, although these surfaces with sulfonic acid (cation exchange particles) and silanol groups (quartz capillary) both carry negative charge density, the differences between ζp and ζw can become noticeable. Apart from an influence on the volumetric EOF rates, it will engender hydrodynamic dispersion which is aggravated at low column-to-particle diameter ratio.41 To analyze the effect of differences in the ζ-potential at the container wall and particle surface on the EOF in confined sphere packings, we carried out simulations for (75) Overbeek, J. Th. G. In Colloid Science; Kruyt, H. R., Ed.; Elsevier: Amsterdam, 1952; Chapter 5. (76) Hunter, R. J.; Wright, H. J. L. J. Colloid Interface Sci. 1971, 37, 564-580. (77) Vindevogel, J.; Sandra, P. J. Chromatogr. 1991, 541, 483-488. (78) Tsuda, T. In Electric Field Applications in Chromatography, Industrial and Chemical Processes; Tsuda, T., Ed.; VCH: Weinheim, 1995; Chapter 3. (79) Rathore, A. S.; Wen, E.; Horva´th, Cs. Anal. Chem. 1999, 71, 2633-2641.

6110

Langmuir, Vol. 21, No. 13, 2005

Hlushkou et al.

Figure 13. Profiles of the local EOF velocity (top) and electrical field strength (bottom) in the selected near-wall (A) and core (B) regions of packing 1 (see Figure 12). Fluid phase properties as in Figure 3, ζw ) ζp ) -50 mV and Eext ) 50 kV/m.

in sphere packing 1. A linear dependence vvol ) vvol(ζw) revealed by Figure 14 qualitatively agrees with that obtained by the model of Rathore and Horva´th.35 However, their model predicts a weaker dependence of the volumeaveraged axial flow velocity on ζw (under otherwise the same conditions) calculated by

vvol ) -

0rEextζp K* 8 µ K∞ d 2 c

∫0r

c

(

[

1+

]

)

ζw I0(2κr/ds) -1 r dr (26) ζp I0(κdc/ds)

Figure 14. Dependence of EOF velocity through sphere packing 1 on ζ-potential at the container wall simulated using the slip-velocity boundary condition. Fluid phase properties as in Figure 3, ζp ) -50 mV and Eext ) 50 kV/m. Velocity is normalized with respect to the uniform ζ-potential distribution (ζw ) ζp ) -50 mV). The solid line represents results obtained with the model of Rathore and Horva´th.35

packing 1 assuming constant ζp ) -50 mV, but varying ζw in the range ζp/2 e ζw e 2ζp. Figure 14 shows the volumeaveraged axial velocity as a function of ζw (normalized by the corresponding value for ζw ) ζp). While ζw is increased from -25 to -100 mV, we observe an only insignificant change of the normalized velocity from around 0.95 to 1.1. This weak dependence on ζw can be explained by the relatively small contribution of the column inner surface (with respect to the particle surface) to the average EOF

where I0 is the zero-order modified Bessel function of the first kind, dc ) 2rc is the diameter of the packed column or capillary, and κ is evaluated by

κ)3

(

)

ϑ(1 - inter) 2

1/2

(27)

ϑ is a parameter depending on packing structure (porosity, tortuosity, and the particle shape) which corrects the differences in drag Stokes force offered by a single (free) hard sphere, multiplied by the number of spheres in a given volume, compared to the actual drag force on the same number of spheres packed in that volume. (It may be noted that the corresponding expressions in the paper of Rathore and Horva´th35 erroneously contain the radius of the particles instead of their diameter.) The first term of the integrand in eq 26 represents the EOF component generated at the particle surface and does not vary with radial position. The second term is related to EOF

Electroosmotic Flow in Porous Media

Langmuir, Vol. 21, No. 13, 2005 6111

however, the inherent dispersive behavior of EOF in confined sphere packings depends, in particular, on the actual differences in ζ-potential at particle (ζp) and wall (ζw) surfaces, as well as on the column-to-particle diameter ratio. An optimization of fixed bed morphology in view of a reduced heterogeneity of the axial EOF velocity can be achieved by adjustment of ζw at low column-to-particle diameter ratio (Figure 15, ζw < ζp) or by a sufficient increase of the latter. 5. Summary

Figure 15. Radial distribution of axial EOF velocity in sphere packing 1 (for different values of the ζ-potential at the inner surface of the confining container, ζw as indicated) simulated using the slip-velocity boundary condition. Fluid phase properties as in Figure 3, ζp ) -50 mV and Eext ) 50 kV/m.

generated at the confining wall and varies with radial position. It has a maximum absolute value at the wall determined by the excess ζ-potential (ζw - ζp) and decays with distance toward the core region of the packing. A discrepancy between the results of our simulations and those obtained with eqs 26-27 (Figure 14) can be explained by the circumstance that the model of Rathore and Horva´th35 considers the random sphere packing as a homogeneous medium characterized by constant mean porosity and tortuosity. As a consequence, also the packing parameter ϑ is inherently assumed as radially constant. On the other hand, ϑ strongly depends on the density of a sphere packing. For example, values of ϑ corresponding to inter ) 0.4 and 0.8 calculated by the model of Happel80

9 9 3 - γ + γ5 - 3γ6 2 2 ϑ) 3 + 2γ5

(28)

with γ ) (1 - inter)1/3, are 85.12 and 5.64, respectively. Hence, radial variations of inter in confined sphere packings affect ϑ and are expected to result in a strong dependence of EOF velocity on radial position. This has been confirmed by our simulations for sphere packings with uniform distribution of the ζ-potential (cf. Figure 8). Equations 26-28 indicate that a local increase or decrease in porosity causes a corresponding change in velocity. Thus, an underestimation of the interparticle porosity in the nearwall region by use of a mean porosity in the model of Rathore and Horva´th35 (and also that adapted by Liapis and Grimes37) results in reduced values of the EOF component close to the confining wall. Although mean porosity generally overestimates local porosity in the core region of a sphere packing, it cannot fully compensate the reduction of EOF velocity in the near-wall region, especially at low column-to-particle diameter ratio (dc/dp). Figure 15 allows an estimation of the effect of differences in the ζ-potential at particle and wall surfaces on an EOF nonuniformity. This figure shows the radial distribution of interparticle axial velocity in packing 1 for different values of ζw (-25, -50, and -100 mV) at constant ζp ) -50 mV. The distributions coincide nearly perfectly as the distance from the wall becomes larger than 0.5ds indicating that the value of ζw affects the axial EOF velocity only in direct vicinity to the confining wall. In general, (80) Happel, J. AIChE J. 1958, 4, 197-201.

Electroosmosis in dense regular and random fixed arrays (packings) of impermeable, nonconducting spherical particles has been numerically studied employing relative sphere diameters, i.e., ratios of the sphere diameter to the thickness of equilibrium EDL in a range 2.272 e ds/λD e 100. Simulation of EOF in an SC array of hard spheres for increasing ds/λD has shown that the average velocity asymptotically approaches a value corresponding to the TDL limit (i.e., ds/λD f ∞). This value has been obtained by simulations using the slip velocity boundary condition at the solid-liquid interface, and it agrees well with the value provided by the macroscopic model for the EOF through a porous medium of Overbeek and Wijga.38 The presented numerical approach allows complete three-dimensional information to be obtained on EOF velocity fields in random porous media. In particular, our analysis of EOF in the SC array and random-close sphere packings indicates that the local (pore-level) axial velocity profile in the interparticle pore space is generally nonuniform, from a fundamental point of view, even for an infinitesimal EDL thickness and a uniform distribution of the ζ-potential. The deviation from the flat, pluglike EOF velocity profile observed in a single and straight, homogeneous channel (in the TDL limit) is caused by a nonuniform distribution of the local electrical field strength which, in turn, can be related directly to local variations in cross-sectional area of the interparticle void space. This dynamics of local EOF in a sphere packing as compared to the much simpler, more readily anticipated singlechannel porous medium realized, e.g., in capillary electrophoresis has a purely morphological origin. Further, for a given porosity of the sphere packing, ζ-potential, and applied electrical field, the heterogeneity of local axial EOF can be influenced by the particle size distribution with respect to the use of monosized spheres. It is related to better filling of void space and results in more homogeneous simulated velocity fields. At the same time, volumetric EOF does not depend noticeably on the particle size distribution. The more ordered structure in a region close to the column wall yet constitutes another origin of EOF velocity heterogeneity in confined sphere packings, resulting in less sinuous flow paths and higher local axial velocities. This is a characteristic caused by the inherent geometry of confined sphere packings and cannot be reduced significantly by choice of a sphere size distribution alone. We have analyzed the local axial EOF velocity in dependence of differences in ζ-potential at the column wall and particle surface. Simulations show that the influence of packing confinement on the EOF heterogeneity, particularly at a low column-to-particle diameter ratio, can be minimized by adjusting the wall ζ-potential. Besides physical insight and understanding gained with respect to the EOF velocity profiles in random sphere

6112

Langmuir, Vol. 21, No. 13, 2005

packings, the results obtained in this work provide implications for a reduction of hydrodynamic dispersion in capillary and chip electrochromatography, or transport of bulk fluid phase through porous media by electroosmosis, in general. Important future studies based on this work are related to the scaling of hydrodynamic dispersion with the average EOF velocity through sphere packings, especially compared to hydraulic flow, as well as the use of porous (permeable, conducting) spheres. In contrast to the hard spheres, these allow intraparticle EOF, buts due to EDL interaction in sufficiently small intraparticle poressdevelop ion-permselective transport behavior as

Hlushkou et al.

well which leads to concentration polarization around the charge-selective particles.81 Acknowledgment. This work was supported by the Deutsche Forschungsgemeinschaft (Bonn, Germany) under Grant TA 268/1-1 and by Fonds der Chemischen Industrie (Frankfurt a.M., Germany). LA050239Z (81) Tallarek, U.; Leinweber, F. C.; Nischang, I. Electrophoresis 2005, 26, 391-404.