Hydrodynamic Dispersion due to Combined Pressure-Driven and

Jump to Dispersion due to the Combined Flow: General Expression - ... the local velocity of the combined flow, u, ... due to the purely pressure-drive...
1 downloads 0 Views 171KB Size
Download PDF
Anal. Chem. 2004, 76, 2708-2718

Hydrodynamic Dispersion due to Combined Pressure-Driven and Electroosmotic Flow Through Microchannels with a Thin Double Layer Emilij K. Zholkovskij

Institute of Bio-Colloid Chemistry of Ukrainian Academy of Sciences, Vernadskogo, 42, 03142, Kiev, Ukraine Jacob H. Masliyah*

Department of Chemical and Materials Engineering, 536 Chemical-Mineral Engineering Building, University of Alberta, Edmonton, Alberta, Canada T6N 2G6

The hydrodynamic dispersion of a nonadsorbed and nonelectrolyte solute is considered for the case of a flow driven through a straight microchannel by pressure and electric potential differences. The analysis is conducted using a thin double layer approximation developed in the previous paper (Zholkovskij, E. K.; Masliyah, J. H.; Czarnecki, J. Anal. Chem. 2003, 75, 901-909). On the basis of this approach, an expression is derived to address the dispersion coefficient for arbitrary electrokinetic potential, electrolyte type, and cross-section geometry. In the derived expression, the influence of cross-section geometry manifests itself through the channel hydrodynamic radius and through three dimensionless geometrical factors. The procedure for obtaining the geometrical factors is presented for an arbitrary cross-section geometry. The geometrical factors are evaluated for several examples of cross section: (i) unbounded parallel planes; (ii) circle; (iii) annulus; (iv) ellipse; (v) rectangle. The dependency of the dispersion coefficient on different parameters is discussed. It is shown that the dependencies are substantially affected by the cross-section geometry, electrolyte type, and electrokinetic potential. Electrically driven (electroosmosis) flow of a liquid through a capillary is widely employed in microfluidic devices as a mean for solute transport and separation.1-5 As compared with transport in pressure-driven (Poiseuille) flow, an important advantage of electroosmosis flow is associated with a substantially weaker longitudinal spreading (hydrodynamic dispersion) of the trans* Corresponding author. E-mail: [email protected]. (1) Harrison, D. J.; Fluri, K.; Seiler, K.; Fan, Z.; Effenhauser, C. S.; Manz, A. Science 1993, 261, 895-897. (2) Jacobson, S. C.; Hergenroder, R.; Kounty, L. B.; Ramsey, J. M. Anal. Chem. 1994, 66, 1114-1118. (3) Manz, A.; Effenhauser, C. S.; Burgraf, N.; Harrison, D. J.; Seiler, K.; Fluri, K. J. Micromech. Microeng. 1994, 4, 257-265. (4) Jacobson, S. C.; Ramsey, J. M. In High Performance Capillary Electrophoresis; Khaledi, M. G., Ed.; Join Wiley & Sons: New York, 1998; Vol. 146, pp 613633. (5) Gulbertson, C. T.; Jacobson, S. C.; Ramsey, J. M. Anal. Chem. 2000, 72, 5814-5819.

2708 Analytical Chemistry, Vol. 76, No. 10, May 15, 2004

ported molecules. Hydrodynamic dispersion occurs due to the nonuniformity of the velocity distribution within the channel cross section. Electroosmosis produces relatively weak dispersion since the electroosmotic velocity is nearly uniform within the cross section except the thin interfacial electric double layer region. At the same time, the parabolic-type Poiseuille flow is strongly nonuniform and, consequently, produces a much stronger dispersion. Thus, using electroosmosis instead of pressure-driven flow, the hydrodynamic dispersion is reduced thereby reducing peak broadening in microfluidic analytical methods.6,7 However, while dealing with microchannel networks, it is difficult to create a purely electroosmotic flow. Due to different mechanisms,7 electroosmosis is usually accompanied by a pressure-driven flow component. Although the parasitic Poiseuille flow is expected to be much weaker than the “main” electroosmotic flow, nevertheless, it gives a noticeable and, even, a dominant contribution to the hydrodynamic solute spreading. Theoretical analysis of the dispersion due to such a combined (electroosmotic and pressure-driven) flow is the objective of the present paper. Mathematical description of the hydrodynamic dispersion is often conducted using the Taylor-Aris theory.8-10 Taylor and Aris showed that, from a formal mathematical point of view, the hydrodynamic spreading is similar to longitudinal solute diffusion. Dynamics of such an apparent diffusion is characterized using an apparent diffusion (dispersion) coefficient, D*, which can be evaluated as

D* ) D +

〈[u - 〈u〉] [ω - 〈ω〉]〉 D

(1)

where u ) u(rb) is the velocity distribution within the cross section (Figure 1) (the radius vector, br ) bıxx + bıyy, belongs to the crosssection plane, x and y are the 2D Cartesian coordinates in the (6) Junno, K.; Sawada, H. Trends Anal. Chem. 2000, 19, 352-363. (7) Gas, B.; Stedry, M.; Kenndler, E. Electrophoresis 1997, 18, 2123-2133. (8) Taylor, G. Proc. R. Soc. London 1953, 219A, 186-203. (9) Aris, G. Proc. R. Soc. London 1956, 235A, 67-77. (10) Aris, G. Proc. R. Soc. London 1959, 252A, 538-550. 10.1021/ac0303160 CCC: $27.50

© 2004 American Chemical Society Published on Web 04/15/2004

interface, Pismen and Babchin,16 Martin and Guiochon,17 Martin et al.,18 Datta,19 McEldoon and Datta,20 and Griffiths and Nilson-24 evaluated the electroosmotic dispersion coefficient DE* for the circular and plane parallel channels. For high interfacial potentials, the electroosmotic dispersion was addressed by Andreev and Lizin,21,22 Gas et al.,23 Griffiths and Nilson,25 and Zholkovskij et al.26 In the above-mentioned papers,8-26 the obtained dispersion coefficient differs from the molecular diffusion coefficient, and the normalized value of this difference can be represented in the form

δDP,E )

cross-section plane as shown in Figure 1), D is the molecular diffusion coefficient of the transported nonelectrolyte solute, and 〈...〉 ) (∫s... ds)/S is the mean value over the cross-section area, S. The function ω ) ω(rb) satisfies the equation

(2)

For nonadsorbed solute, eq 2 is subject to the boundary condition

∇ B ω‚n b)0

)

2

(4)

where the superscripts P and E signify quantities attributed to the purely pressure-driven or electrically driven flow, respectively For the purely pressure-driven flow, we will express the lengthscale parameter dP in the form

Figure 1. Microchannel cross section.

∇2ω ) 〈u〉 - u

(

DP,E 〈uP,E〉dP,E / - D ) D D

at the cross-section perimeter line (3)

In (2) and (3), the gradient operator is given by ∇ B ) bıx∂/∂x + bıy∂/∂y; b n is the unity outward normal vector (Figure 1). Using the procedure given by (1)-(3) for a purely pressuredriven flow through channels, many authors analyzed the corresponding dispersion coefficients, DP* . During these studies, different cross-section geometries have been examined, for example, circle,8 ellipse,9 annulus,10 unbounded parallel planes,11 rectangle,12,13 etc. A general theory describing dispersion due to flow through a chromatographic column having arbitrary morphology was suggested by Giddings.14,15 Similarly to Aris,10 Giddings took into account redistribution of the solute between the mobile and stationary phases. This effect plays a mayor role in chromatography. A number of publications were concerned with the hydrodynamic dispersion in purely electrically driven (electroosmosis) flow.16-26 Assuming low electric potentials at the wall/solution (11) Dewey, R.; Sullivan, P. M. Z. Angew. Math. Phys. 1978, 30, 601-620. (12) Doshi, M. R.; Dalya, P. M.; Gill, W. N. Chem. Eng. Sci. 1978, 33, 795-804. (13) Chatwin, P. C.; Sullivan, P. J. J. Fluid Mech. 1982, 120, 347-358. (14) Giddings, J. C. J. Chromatogr. 1961, 5, 46-60. (15) Giddings, J. C. Dynamics of Chromatography; Edward Arnold Ltd.: London, 1984. (16) Pismen, L. M.; Babchin, A. J. J. Colloid Interface Sci. 1977, 62, 63-68. (17) Martin, M.; Guiochon, G. Anal. Chem. 1984, 56, 614-620. (18) Martin, M.; Guiochon, G.; Walbroehl, Y.; Jorgenson, J. W. Anal. Chem. 1985, 57, 559-561. (19) Datta, R. Biotechnol. Prog. 1990, 6, 485-493. (20) McEldoon, J. P.; Datta, R. Anal. Chem. 1992, 64, 227-230. (21) Andreev, V. P.; Lizin E. E. Electrophoresis 1992, 13, 832-837. (22) Andreev, V. P.; Lizin E. E. Chromatographia 1993, 37, 202-210. (23) Gas, B.; Stedry, M.; Kenndler, E. J. Chromatogr., A 1995, 709, 63-68.

(dP)2 ) KPh2

(5)

The parameter h ) S/l (where l is the cross-section perimeter) is referred to as the hydrodynamic radius of a channel. The dimensionless factor KP is defined by the cross-section geometry. Using (1), (4), and (5), the geometrical factor KP can be represented as

KP )

〈[

]

1 uP - 1 [ωP - 〈ωP〉] h2 〈uP〉

〉

(6)

For the purely electrically driven flow through the channel with thin double layer, Griffiths and Nilson25 suggested a simple relation for the length-scale parameter dE

(dE)2 )

g(ζ)KE κ2

(7)

where κ is the Debye parameter of the electrolyte

∑C Z

F2 κ2 )

2

i i

i

(8)

RT Here, F is the Faraday constant, Zi and Ci are, respectively, the ith ion charge (in Faraday units) and molar ionic concentration in an equilibrium solution,  is the solution dielectric permittivity, R is the gas constant, and T is the absolute temperature. In (7), g(ζ) is a geometry-independent function of the interfacial potential, ζ, which is also referred to as the electrokinetic or ζ-potential. For the case of symmetric electrolyte, Griffiths and Nilson25 suggested an interpolation for g(ζ). (24) Griffiths, S. K.; Nilson, R. H. Anal. Chem. 1999, 71, 5522-5529. (25) Griffiths, S. K.; Nilson, R. H. Anal. Chem. 2000, 72, 4767-4777. (26) Zholkovskij, E. K.; Masliyah, J. H.; Czarnecki, J. Anal. Chem. 2003, 75, 901-909.

Analytical Chemistry, Vol. 76, No. 10, May 15, 2004

2709

In our previous paper,26 using the approximation of a thin double layer, we obtained the following expression for the function g(ζ)

g(ζ˜ ) )

[

1

ζ˜

∫

Φ

ζ˜

0

x∑ 2

]

2

dΦ

νi[ exp(-ΦZi) - 1]

i

(9)

DISPERSION DUE TO THE COMBINED FLOW: GENERAL EXPRESSION Due to the linearity of the Stokes equations, the local velocity of the combined flow, u, can be represented as

u ) uP + u E

where ζ˜ ) Fζ/RT and νi ) Ci/∑nCnZn2. Expression 9 is valid for an arbitrary electrolyte system, and for symmetric electrolytes, it coincides with g(ζ) obtained by Griffiths and Nilson.25 The dimensionless factor, KE, in eq 7 is defined by the crosssection geometry. Griffiths and Nilson25 showed that KE ) 1/2 for the circular cross sections and KE ) 1/3 for the unbounded parallel planes. In the previous paper,26 we derived an expression that describes the geometrical factor KE for arbitrary cross-section geometry, namely

KE ) 〈W〉l - 〈W〉

geometry. Then, we will discuss how the dispersion coefficient depends on the geometry and on different parameters.

Combining (1)-(3) and (13) yields the following expression for the normalized variation of the dispersion coefficient

d×2

δD* ) δDp + δDE + 2〈up〉〈uE〉

d×2 )

∇ ˜ 2W ) 1

within the cross-section region

(11)

which is subject to the boundary condition

B ∇ ˜ W‚n b)1

at the perimeter line

(12)

B In (10) and (11), B r˜ ) b/h; r ∇ ˜ ) h∇ B and ∇ ˜ 2 ) h2∇2. Using (10)(12), the geometrical factors KE have been predicted for different cross-section geometries.26 In contrast with the above-mentioned theories,8-26 the present paper is concerned with the hydrodynamic dispersion for the case when both pressure-driven and electroosmotic flows are present simultaneously. The dispersion due to such a combined flow was analyzed by Datta and Kotamarthi,27 who considered a circular tube and a low interfacial potential. However, very often, the electrokinetic potentials are not low23,24 and the cross-section shapes substantially differ from a circle,1-5 which can substantially affect the dispersion. Therefore, the present paper is aimed at describing the dispersion coefficient for arbitrary interfacial potentials and for different cross-section geometries. Similarly to our previous paper,26 the present analysis will be based on the thin double layer approximation that drastically simplifies solution of the complex general problem. We will use the fact that the length-scale parameters of microchannel cross sections are much larger than the Debye length, κ-1. The present paper is organized according to the following structure. After obtaining some relationships, which are valid for arbitrary straight channels, we will predict the dispersion coefficient for the specific case of the channel having a thin double layer. Using the obtained general results, which are applicable for arbitrary cross-section geometry, we will consider several specific examples of the (27) Datta, R.; Kotamarthi, V. R. AIChE J. 1990, 36, 916-926.

2710

Analytical Chemistry, Vol. 76, No. 10, May 15, 2004

D2

(14)

where δDP and δDE are given by (4). The length-scale parameter d× can be expressed in a double manner

(10)

where 〈...〉 ) (Il... dl)l is an average over the cross-section perimeter, l. The function W ) W(rb) is obtained as a solution of the boundary problem containing the equation

(13)

1 〈[uP - 〈uP〉][ωE - 〈ωE〉]〉 ) 〈uP〉〈uE〉 1 〈[uE - 〈uE〉][ωP - 〈ωP〉]〉 (15) 〈uP〉〈uE〉

where ωP and ωE are obtained as solutions of the problem given by (2) and (3) after substituting there u ) uP and u ) uE, respectively. The latter equality in (15) is proved in Appendix A of the Supporting Information section. Thus, according to (14), for obtaining the dispersion coefficient due to the combined flow, one can use the results of all the previous publications,8-26 which give the variations of the dispersion coefficient due to the purely pressure-driven or electroosmotic flows (respectively, δDP and δDE). Additionally, one should obtain the cross term represented by the third term on the right-hand side of (14). The cross term can be obtained using any of equalities (15). To derive an expression for the cross term, we will represent the velocity distribution within the cross section in a widely used form24-31 that is applicable for arbitrary cross-section geometry26

uE(b) r )-

E [ζ - Ψ(b)] r η

(16)

where Ψ(rb) is the potential distribution within the cross section in the absence of an external electric field, E is the longitudinal electric field strength, and η is the liquid viscosity. Equation 16 is valid if the equilibrium potential takes a constant value at the interface

Ψ)ζ

at the cross-section perimeter line

(17)

where ζ is referred to as the electrokinetic potential. Substituting (28) Burgreen, D.; Nakache, F. R. J. Phys. Chem. 1964, 68, 1084-1091. (29) Rice, D.; Whitehead, R. J. Phys. Chem. 1965, 69, 4017-4024. (30) Sorensen, T. S.; Koefoed, J. J. Chem. Soc., Faraday Trans. 2 1974, 70, 665675. (31) Yang, C.; Li., D. J. Colloid Interface Sci. 1997, 194, 95-107.

(16) into the latter expression in (15) yields

d×2 ) 〈[〈ωP〉 - ωP]Ψ〉/[ζ - 〈Ψ〉]〈uP〉

(18)

Thus, to evaluate d× for given cross-section geometry, one should obtain the distributions ωp(rb), up(rb), and Ψ(rb) and substitute them into (18). APPROXIMATION OF THIN DOUBLE LAYER Using the approach of our previous paper,26 one can simplify the integrals associated with obtaining the cross-sectional mean values in (18). Consequently, we will consider a function Φ(ξ), which is set by a reversion of the following integral

∫

ξ ) sign(ζ)

dφ

ζ

Φ

x∑ 2

(19)

νk[exp(-φZk) - 1]

k

Ψ ) (RT/F)Φ(ξ)

ξ ) -κ(b r l - b)‚n r b

(21)

(iii)when obtaining the integrals of the structure ∫κ∆ b)Φdξ, to 0 f(r replace the upper limit, κ∆, by infinity. The above-described procedure yields the following expressions for the mean values present in (18)

〈(〈ωP〉 - ωP)Ψ〉 ≈

(

(24)

where the new dimensionless geometrical factor K× is expressed as

K× ) [〈ωP〉 - 〈ωP〉l]/〈uP〉h2

(25a)

where ωP(rb) is a solution of problem (2) and (3) formulated for the purely pressure-driven flow (u ) uP). It can be shown that K× g 0 regardless the sign of uP. Combining (25a) with the 2D version of the Greens theorem, one can derive another useful expression for K×

K× ) 〈[1 - uP/〈uP〉]W〉

(25b)

(20)

where the normalized coordinate ξ is defined with reference to a point at the cross-section perimeter line (Figure 1);

RT [ Fκh

ωP is considered to be independent of ξ. Combining (9) and (19)(24) one obtains

h d×2 ) [K×xg(ζ) + O(1/κh)] κ

The function Φ(ξ) given by (19) is a solution of the onedimensional (flat) version of the normalized Poisson-Boltzmann problem for the semi-infinite region. Under the condition κh . 1, Φ(ξ) can be used for approximating the function Ψ(rb) in the vicinity of the interface The approximate procedure employed in our previous paper26 enables a substantial reduction of the integrals having the structure ∫S f(rb)Ψ(rb) dS that occur while obtaining the crosssectional mean values. The essence of the procedure follows: (i) to replace the integration over the cross-section area, S, by the integration over the thin stripe with the thickness ∆ (κ-1 , ∆ , h) adjacent to the perimeter line (Figure 1); (ii) while integrating over the stripe, to approximate the equilibrium potential distribution, Ψ(rb), using the function Φ(ξ) given by (19)

〈Ψ〉 )

Figure 2. Illustration to definition of parameter τ.

∫ Φ(ξ) dξ + O(1/κh)] ∞

0

)∫

RT 1 〈ωP〉 - I ωP dl [ Fκh ll

(22)

∞

0

Φ(ξ)dξ +

O(1/κh)] (23) where ωP is a solution of the problem given by (2) and (3), where u ) uP. Hence, the length scale for the changes in ωP is h. The latter means that, within a stripe having thickness ∆ , h and adjacent to the interface (Figure 1), the function ωP is nearly constant along the normal ξ-axis. Therefore, while obtaining (23),

Details of this derivation are given in Appendix B (Supporting Iinformation). Thus, when W(rb/h) is a known function,26 using (25b) one can evaluate K× without preliminary obtaining ωP(rb). Combining (4)-(7), (13), (14), and (25), the ratio of the dispersion coefficient to the molecular diffusion coefficient is represented as

D*/D ) 1 + Pe2[τ2(κh)2KP + 2τ(1 - τ)κhK×xg(ζ) + g(ζ)KE + O(τ)] (26) where, for a given mean velocity, 〈u〉, the Peclet number is a geometry independent

Pe ) 〈u〉/Dκ

(27)

The parameter τ in (26) signifies a relative contribution of the Poiseuille into total flow

τ ) 〈uP〉/〈u〉

(28)

Illustrations given in Figure 2 show that τ can vary from -∞ to ∞. When the condition κh . 1 is satisfied, eqs (26)-(28) are Analytical Chemistry, Vol. 76, No. 10, May 15, 2004

2711

capable of describing the dispersion coefficient, D*, for any given electrolyte solution, channel cross-section geometry, and ζpotential. It should be noted that, since κh . 1, the omitted term, O(τ), on the right-hand side of (26) yields a negligible contribution for all values of τ. Thus, eq 26 yields an asymptotic expression that is valid for the whole range of variation of τ. Dependency on the electrokinetic potential is embedded in the function g(ζ) given by (9). As was discussed earlier,26 the behavior of g(ζ) strongly depends on the electrolyte type. The electrolyte type also affects the Debye parameter κ (eq 8). Three geometrical factors KP, KE, and K× and the hydrodynamic radius h are completely defined by the cross-section geometry. To evaluate the geometrical factor KP one should employ the Taylor-Aris procedure,8-13 which amounts to substituting into eq 6 the Poiseuille flow velocity distribution, uP, and a solution of problem (2) and (3), ωP. As stated earlier,24 the coefficient KE is obtained by substituting a solution of problem (11) and (12), W, into (10). The coefficient K× is evaluated using eqs 25a or 25b, where one should substitute, respectively, ωP or W and uP. EXAMPLES OF CROSS-SECTION GEOMETRY Studies of the hydrodynamic dispersion for different crosssectional geometries create a basis for microchannel optimization.33-35 We will analyze the geometrical factors KP, KE, and K× for five types of channel cross section: (i) unbounded parallel planes; (ii) circle; (iii) annulus; (iv) ellipse; (v) rectangle. For all these geometries, the purely hydrodynamic dispersion was described earlier.8-15 In the present section, combining these results with eq 5, we will obtain the corresponding values of KP. While obtaining KP from the Giddings results,14,15 it is necessary to substitute unity for the retention factor. Such a substitution corresponds to the case of a nonadsorbed solute. Except for ellipse, the coefficient KE was predicted for all the above geometries.26 Below, we will obtain KE for an ellipse and present its value for other geometries.26 As well, using eq 25a or 25b, we will obtain the third geometrical factor, K × , for all the above listed geometries Unbounded Parallel Planes. Using results of previous publications11,14,15,25,26 one can write

h ) b/2;

KP ) 2/105;

KE ) 1/3

(29)

where b is the distance between the parallel planes For the present geometry, the Poiseuille flow velocity distribution is well known

uP ) 3〈u〉[1 - (2x/b)2]/2

(30)

Combining (2) and (30) and solving problem (2) and (3) one obtains ωP

ωP ) 〈u〉(x/2)2[2(x/b)2 - 1]

(31)

Here, the X-axis is normal to the planes and its origin coincides (32) Landau, L. D.; Lifshitz, E. M. Fluid Mechanics, 2Nd ed. revised; Pergamon Press: Oxford, U.K., 1987.

2712

Analytical Chemistry, Vol. 76, No. 10, May 15, 2004

with the middle of the cross section. Combining (25a) and (31) one obtains

K× ) 1/15

(32)

Thus, eqs 29 and 32 yield h, KP, KE, and K× for unbounded parallel planes sections. Circle. For a circular cross section, the dispersion due to the purely Poiseuille flow was predicted in the classic work of Taylor.8 For thin double layer, the dispersion due to the electroosmotic flow was considered by Griffiths and Nilson25 and Zholkovskij et al.26 According to these references, the corresponding coefficients are

h ) a/2

KP ) 1/12

KE ) 1/2

(33)

where a is a circle radius. The function ωP(r) (r is the distance from the circle center) is obtained as a solution of (2) and (3), where the velocity distribution is expressed as

u ) 2〈u〉[1 - (r/a)2]

(34)

Consequently, ωP(r) takes the form

ωP ) 〈u〉r2[(r/a)2 - 2]/8

(35)

Combining (25a) and (35) one obtains

K× ) 1/6

(36)

Thus, eqs 33 and 36 yield coefficients h, KP, KE, and K× for a circular cross section. Let us now compare the above obtained result for the circular cross section with the expression derived by Datta and Kotamarthi27 for the same geometry and for the low interfacial potentials. In their study, the analysis was not restricted to large κa as in our case. The final result of Datta and Kotamarthi27 (eqs 37-42 from their paper) can be rewritten in the form that coincides with a particular version of eq 26 corresponding to the case of

g(ζ) ) 1

(37a)

The geometrical factors KE and K×, which are obtained using such a transformation of the Datta and Kotamarthi E result, become functions of the parameter κa, KD-K (κa) (33) Dutta, D.; Leighton T. D., Jr. Anal. Chem. 2001, 73, 504-513. (34) Dutta, D.; Leighton T. D., Jr. Anal. Chem. 2002, 74, 1007-1016. (35) Dutta, D.; Leighton T. D., Jr. Anal. Chem. 2003, 75, 57-70.

section geometry, the coefficient KE was obtained in our previous paper.24 Accordingly,

h ) (a1 - a2)/2 KE ) [γ2 - 4γ + 1 + 4γ2 ln(γ)/(γ2 - 1)]/2(γ - 1)2 (38) where γ ) a1/a2. As was shown,24 the expression for KE (38) yields the expected limiting cases: when γ f 1 (unbounded parallel planes), the geometrical factor KE f 1/3 (see eq 29); when γ f ∞ (circle), KE f 1/2 (see eq 33). For an annulus, the dispersion in purely pressure-driven flow was analyzed by Aris10 and Giddings.15 Using their results one can represent KP, as

KP ) [72γ4 ln3(γ) + 3(γ4 - 1)(γ4 + 8γ2 + 1) ln2(γ) 2(γ2 - 1)2(5γ4 + 32γ2 + 5) ln(γ) + 9(γ2 + 1)(γ2 - 1)3]/ [36(γ - 1)2(γ2 - 1)[γ2 - 1 - (γ2 + 1) ln(γ)]2] (39) Figure 3. Circular cross section: comparison with results of Datta and Kotamarthi.27

andKD-K × (κa), given as

KD-K ) × E ) KD-K

[

[

]

2κaI1(κa) 2 1 16 + κaI0(κa) - 2I1(κa) 12 (κa)2 (κa)3

2κaI1(κa)

κaI0(κa) - 2I1(κa)

(37b)

To obtain K× we will use (25a), where we will substitute the function ωP being a solution of (2) and (3). In the polar coordinate system with polar radius r, the velocity distribution is given by the well-known expression30

uP ) 2〈uP〉[(1 - F2)γ2 ln(γ) + (γ2 - 1) ln(F)]/ [(γ2 + 1) ln(γ) - γ2 + 1] (40)

]{ 2

[ ]}

I0(κa) I0(κa) 2 3 + 8 (κa)2 2κaI1(κa) 2I1(κa)

2

(37c) Here, I0(κa) and I1(κa) are the first kind modified Bessel functions of the zero and first order, respectively. According to the analysis conducted earlier,26 at ζ˜ f 0, the function g(ζ) f 1. Hence, for low interfacial potentials, condiE tion 37a is satisfied. Behavior of the functions KD-K (κa) and D-K K× (κa) given by (37b) and (37c) is illustrated in Figure 3. With increasing κa, the functions approach, respectively, 1/2 and 1/6, i.e., the values of the geometrical factors given by eqs 33 E and 36. Starting with κa = 70 (Figure 3, point A) KD-K (κa) deviates from 1/2 by less than 4%. The geometrical factor KD-K × (κa) reaches a value of 1/6 at κa = 11.3 (Figure 3, point B). It is interesting that, with further increase in κa, KD-K × (κa) also increases as far as it reaches a maximum at κa = 24. The KD-K × (κa) maximum slightly exceeds 1/6 (∼2.8%). When κa increases at κa > 24, the curve KD-K × (κa) asymptotically approaches 1/6. Thus, for low interfacial potentials, the present paper result given by (26)-(28), (33), and (36) is exactly transformed into the asymptotic form (κa . 1) of the Datta and Kotamarthi result.27 However, due to the properties of g(ζ),24 at large or moderate ζ˜ , the Datta and Kotamarthi result27 can substantially deviate from the predictions given by eqs 26-28, 33, and 36. Annulus. Now, we will consider a channel formed by coaxial circular cylinders with radii a1 and a2 (a1 > a2). For this cross-

where F ) r/a1. Combining (2), (3), and (40) one obtains ωP as a solution of (2) and (3)

ωP )

[

h2γ2〈u〉 2

2

(γ - 1) [γ - 1 - (γ2 + 1) ln(γ)]

× (γ2 - 1)(ln(γ) - 1)F2 -

F4 2 γ ln(γ) + 2

]

2(γ2 - 1)F2 ln(F) + 2 ln(γ)ln(F) (41)

Consequently, combining (25a) and (41) yields K×:

K× ) [24γ3 ln2(γ) - 2(γ3 + 1)(γ - 1)3 ln(γ) + 3(γ2 - 4γ + 1)(γ2 - 1)2]/ [12(γ - 1)2(γ2 - 1)[γ2 - 1 - (γ2 + 1) ln(γ)]] (42)

As was expected, at γ f 1, eqs 39 and 42 yield, respectively, KP f 2/105 and K× f 1/15; i.e., both coefficients approach the values inherent in the case of the unbounded parallel planes (eqs 29 and 32). Considering the limiting transition γ f ∞, one obtains from (39) and (42) KP f 1/12 and K× f 1/6, which corresponds to the case of the circular cross section (eqs 32 and 34). The behavior of the geometrical factors KP KE, and K× as functions of the parameter γ ) a1/a2 is illustrated in Figure 4 where all the curves display an asymptotic approach to the Analytical Chemistry, Vol. 76, No. 10, May 15, 2004

2713

h and geometrical factor KP (eq 5), the Aris result can be expressed as

h ) a/B0(λ) ) b/B0(1/λ) KP ) B02(λ)(5 + 14λ2 + 5λ4)/576λ2(1 + λ2) (44) Here, λ ) a/b is the ratio of the ellipse semi-axises, a and b. Expression for KE is derived by conducting the procedure given by (10)-(12) in the elliptic cylindrical coordinate system. Corresponding coordinates, µ and χ, are interrelated with the normalized 2D Cartesian coordinates (Figure 5a), as

x ˜ ) x/h ) B0x1 - λ-2 cosh(µ) cos(χ); ˜ ) y/h ) B0x1 - λ-2 sinh(µ) sin(χ) (45) y Consequently, a solution of the problem (11) and (12) for the elliptic cross sections, W(µ,χ), takes the form

{

W ) B0 1

(λ2 - 1) [cosh(2µ) + cos(2χ)] + 8λ2 (λ2 - 1)k

∞

∑

Bk cosh(2kµ) cos(2kχ)

λ k)1k[(λ + 1)2k - (λ - 1)2k]

}

(46)

Combining (10), (43), and (46), one obtains the expression for KE Figure 4. Geometrical factors for annulus as a function of γ ) a1/ a2: (a) KP(γ); (b) KE(γ) and K×(γ).

above-discussed limits. It is interesting that, with an increase in γ, the coefficient KP(γ) very gradually approaches the value, 1/12, corresponding to the circular cross section. In Figure 4b, even for γ ) 107 (point B), the coefficient KP is less than 1/12 by 10%. At smaller γ, one can observe a much sharper increase in KP with γ. For example, when γ ) 10 (point A), KP exceeds its value for the parallel planes, 2/105, more than by a factor of 2. Generally, KP(γ) can vary more than a factor of 4. Therefore, for any realistic γ > 10, the geometrical factor KP predicted for the annulus is expected to be noticeably different from the values attributed to both the parallel planes and circle. With increasing γ from the unity to infinity, each of the geometrical factors KE and K× increases by the factor of 3/2 (Figure 4b). One can observe a more rapid increase in KE which, at γ ) 100, practically coincides with 1/2 (circle). The coefficient K× increases slower than KE but much more rapidly than KP. Ellipse. To represent the geometrical factors KP, KE, and K× in a convenient form, we will introduce a function Bk(λ), as

Bk(λ) )

4 π

∫ x1 + (λ π/2

0

2

- 1) sin2(χ) cos(2kχ) dχ

(43)

For elliptic cross sections, the dispersion in the purely Poiseuille flow was described by Aris.9 In terms of the hydrodynamic radius 2714

Analytical Chemistry, Vol. 76, No. 10, May 15, 2004

E

K )

1 2λ

{

∞

Bk2 (λ + 1)2k + (λ - 1)2k

∑k k)1

+

(λ + 1)2k - (λ - 1)2k

B0 [B0(λ2 + 1) + 8λ 4B1(λ2 - 1)]

}

(47)

Combining (25b), (43), (45), (46), and the well-known expression for the purely pressure-driven velocity distribution within the elliptic cross section32

uP ) 2〈uP〉[1 - (x/a)2 - (y/b)2]

(48)

one obtains the geometrical factor K×, as

K× )

B0

[B0(λ2 + 1) + 2B1(λ2 - 1) - B2(λ2 - 1)2/ 48λ2 2(λ2 + 1)] (49)

Details of obtaining (46), (47), and (49) are given in Appendix C (Supporting Information). Thus, eqs 44, 47, and 49 yield the geometrical factors KP, KE, and K×, for elliptic cross section as functions of the ellipse aspect ratio, λ ) a/b. The behavior of these functions is illustrated in Figure 5b and c. Since KP,E(λ) ) KP,E(1/λ) and K×(λ) ) K×(1/λ), it is sufficient to discuss their variation within the range 0 < λ e 1. When λ f 0, all the geometrical factors, KP, KE, and K×, diverge

those obtained for circle (eqs 33 and 36) and, hence, differ from each other by a factor of 2. A plot of the function KE(λ), which is given by (47), is shown in Figure 5c. For a comparison, we additionally present the corresponding curves for two types of cross-section geometry, a rectangle and a rhombus, discussed in our previous paper.26 For the inscribed rhombus and circumscribed rectangle (Figure 5a), the coefficient KE can be represented as a function of the same aspect ratio λ. For a rectangle, the factor KE was shown26 to be independent of λ (KE(λ) ) 2/3). For a rhombus, the dependency obtained earlier26 can be rewritten as KE(λ) ) (1 + λ2)2/6λ2. When λ f 0, the latter expression yields KE f 1/6λ2 ≈ 0.167/λ2, whereas, using (47) one obtains the asymptotic expression for an ellipse: KE ) (5π2 - 49)/3π2λ2 ≈ 0.0117/λ2. Hence, when λ decreases toward zero, the geometrical factor, KE, predicted for an ellipse, diverges but remains more than in 10 times smaller than KE predicted for an inscribed rhombus (Figure 5a). For an ellipse with λ , 1, the coefficient KE substantially exceeds the constant value, 2/3, obtained for a rectangle.26 However, when λ is close to unity, KE obtained for an ellipse becomes smaller than 2/3, and, for λ ) 1, KE coincides with the value of 1/2 obtained for a circle.25,26 Rectangle. The purely electroosmotic dispersion in rectangular channels was considered in our previous paper26 where the coefficient KE was shown to be constant:

KE ) 2/3

h ) a/2(λ + 1)

(50)

Here, λ ) a/b is the ratio of the rectangle sides, a and b. Dispersion due to the purely pressure-driven flow through the rectangular channels was analyzed by Doshi et al.12 and Chatwin and Sullivan.13 In terms of the geometrical factor KP, the analytical result of Chatwin and Sullivan,13 can be represented as ∞

KP ) 18(λ + 1)2{

∑Q

∞

2

k

+

k)1

∑G

∞

2 m

∞

∑∑H

+8

m)1

2

km

/

k)1m)1

∞

[k2 + (mλ)2]}/π2[1 - (6/λ)

∑s ]

2

n

(51)

n)0

where

[

]

π sn ) tanh (2n + 1) λ /[(2n + 1)π/2]5 2 Figure 5. Elliptical cross section: (a) sketch of cross section (dotted line, inscribed rhombus; dashed line, circumscribed rectangle); (b) KP(λ) and K×(λ); (c) KE(λ).

as O(1/λ2). At λ ) 1, the geometrical factors take the minimum values which coincide with that obtained for a circular cross section. The curves in Figure 5b display the variation of KP and K×. At λ , 1, KP(λ) and K×(λ) take close values. Taking into account (43), the limiting transition λ f 0 in (44) and (49) yields KP f 5/36π2λ2 ≈ 0.0141/λ2 and K× f 11/80π2λ2 ≈ 0.0139/λ2. At the same time, when λ is close to unity, the coefficient KP and K× are noticeably different. When λ ) 1, the coefficients coincide with

Qk )

∞

2λ

∑ k

n)0λ

sn 2

(52)

;

+ [2k/(2n + 1)]2 Gm )

1

{

1

m3λ π2

+

2

}

m2sn

∞

∑

λn)01 - [2m/(2n + 1)]2

∞

Hkm ) λ

∑s /{1 - [2m/(2n + 1)] }{λ 2

n

2

+ [2k/(2n + 1)]2}

n)0

(53)

For a rectangular cross section, it is convenient to obtain K× with the help of (25b). Using the Cartesian coordinate system, whose Analytical Chemistry, Vol. 76, No. 10, May 15, 2004

2715

Table 1. Geometrical Factors for Different Cross-Section Geometries cross-section geometry parallel planes circle annulus (γ )a1/a2 g 1)

KP

2/105 1/12 γf1 2/105 γ ) 10 0.0402 γ ) 100 0.0580 γf∞ 1/12 ellipse (λ ) a/b) λ ) 0.01 140.8 λ ) 0.1 1.479 λ ) 0.5 0.1164 λ)1 1/12 rectangle (λ ) a/b) λ f 0 0.1515 λ ) 0.2 0.1596 λ)1 0.1340

Figure 6. Rectangular cross section: geometrical factors as functions of aspect ratio λ.

origin coincides with a vertex, the velocity distribution within the rectangular cross section can be represented in the form employed in the paper of Chatwin and Sallivan13 uP ) y*(1 - y*) -

∞

∑ π 8

3

n)0

6〈uP〉

cosh[(2n + 1)π(x* - λ/2)]

[

π (2n + 1)3 cosh (2n + 1) λ 2 1-

]

sin[(2n + 1)πy*]

∞

∑f

6

λ n)0

n

(54) where x* ) x/b, y* ) y/b, and the X-axis and Y-axis are parallel to the rectangle sides a and b, respectively. The earlier obtained solution of problem (11) and (12)26 can be represented as

W ) 2(1 + 1/λ)[(x* - λ/2)2/λ + (y* - 1/2)2]

(55)

Combining (25b), (52), (54), and (55) (see Appendix D of Supporting Information), one obtains

K× ) 1+λ ∞

15λ(λ - 6

∑

{

λ - 6 + 30

sn)

∞

∑s

[

n

n)0

λ+1+

]}

12(1 - λ) (2n + 1)2π2λ

n)0

(56)

Equations 50, 51, and 56 yield the geometrical factors KP, KE, and K×, for a rectangular cross section as a function of the rectangle aspect ratio, λ ) a/b. The behavior of these functions, for 0 e λ e 1, is illustrated in Figure 6. As shown earlier,26 the geometrical factor KE is constant. Two other coefficients, KP and K×, only 2716 Analytical Chemistry, Vol. 76, No. 10, May 15, 2004

KE

K×

KE - K×2/ KP

1/3 1/2 1/3 0.434 0.4907 1/2 118.3 1.741 0.5341 1/2 2/3 2/3 2/3

1/15 1/6 1/15 0.1120 0.1428 1/6 124.1 1.440 0.2008 1/6 0.2767 0.2720 0.2475

1/10 1/6 1/10 0.122 0.139 1/6 8.919 0.339 0.188 1/6 0.161 0.203 0.210

slightly vary with λ. At λ f 0, neither of the coefficients approaches the value corresponding to the unbounded parallel planes. At λ ≈ 0.2, the dependency of the geometrical factor KP on λ reaches a weakly expressed maximum, KP ≈ 0.1596. When λ approaches zero, the coefficient KP approaches 0.1515. For a square (λ ) 1), KP ) 0.134. Thus, the maximal deviation of KP, from its value corresponding λ ) 0, does not exceed 12%. The latter does not contradict the paper of Chatwin and Sullivan13 where the ratio of the dispersion coefficients attributed to λ ) 1 and λ ) 0 is ∼1/4. When the aspect ratio, λ, increases from 0 to 1, the hydrodynamic radius, h, decreases by a factor 1/2. Therefore, according to (4) and (5), for a purely Poiseuille flow, the dispersion coefficient decreases by a factor of ∼1/4 despite the slight variation of KP. When λ increases from 0 to 1, the geometrical factor K× ∞ monotonically decreases from 1/15 + 64∑n)0 1/[π(2n + 1)]5 = 0.2767 to, approximately, 0.2475. Thus, the variation of K× is ∼10.6%. For the above-discussed geometries, some values of the geometrical factors KP, KE and K× are collected in Table 1. DEPENDENCY OF THE DISPERSION COEFFICIENT ON DIFFERENT PARAMETERS In this section, we will discuss the behavior of the normalized dispersion coefficient, D*/D, which is given by eq 26, as a function of the squared Peclet number (Pe2), of the normalized hydrodynamic radius (κh), and of the flow regime parameter (τ). To discuss how the hydrodynamic dispersion depends on different parameters, we will consider the behavior of a normalized variation of the dispersion coefficient, namely, (D*/D - 1)/Pe2. According to eq 28, the parameter Pe2 in the denominator of the latter expression can vary in a wide range. To illustrate this, we will consider the case when the total mean velocity, 〈u〉, is of the same order as the electroosmotic mean velocity, 〈uE〉. For large organic molecules (D = 10-11m2/s) in distilled water (κ = 107 m-1;  = 7 × 10-10 J/mV2; η = 10-3 J‚s/m3), at a reasonably high electrokinetic potential (ζ = 50 mV) and electric field (E = 1 kV/ cm), the squared Peclet number takes a high value, Pe2 = 103, Under the same conditions, but for smaller molecules (D = 10-10 m2/s), eq 28 yields Pe2 = 10. Considering a millimolar electrolyte solution (κ = 108 m-1), one obtains Pe2 = 10 for D = 10-11 m2/s and Pe2 = 0.1 for D = 10-10 m2/s.

To display the behavior of (D*/D - 1)/Pe2 for the geometries discussed in the previous section, the curves presented in Figure 7 were plotted using eq 27 with the geometrical factors, KP, KE, and K×, given in Table 1. Similarly to the result of Datta and Kotamarthi,25 eq 27 represents (D*/D - 1)/Pe2 as a quadratic polynomial in terms of the parameter, τ. This polynomial is always positive and, at τ ) τmin < 0, takes a minimum value: (D*/D 1)min/Pe2. Applying a standard procedure, one obtains from (27)

τmin ) -

(

{

2 E P xg(ζ) 2K× - K K 1 K×xg(ζ) 1 + + O[1/(κh)2] κh κh KP KPK

(

}

(57)

)

D*min - 1 /Pe2 ) D g(ζ) KE -

×

){

K×2 K

P

1+

}

2 K× xg(ζ) + O[1/(κh)2] κh KP

(58)

The case of negative τ is illustrated in Figure 2c where, simultaneously, (i) the pressure-driven flow is directed against the electroosmotic flow and (ii) the mean velocity of the pressuredriven flow has a smaller magnitude than the mean electroosmotic velocity. For low electrokinetic potentials (g(ζ) f 1),26 the behavior of (D*/D - 1)/Pe2 as a function of τ is illustrated by the curve families of Figure 7. Each of the families was plotted for a given κh (κh ) 50 in Figure 7a and κh ) 100 in Figure 7b) and for different cross-section geometries. In both panels, the parameter τ varies within the range -0.1 e τ e 0.1. The latter inequality means that the pressure-driven flow velocity satisfies the following inequality: -〈uE〉/11 e 〈uP〉 e 〈uE〉/9 (Figure 2); i.e., |〈uP〉| , |〈uE〉|. Such regimes are of especial interest since, very often, in the microchannel networks, the Poiseuille flow is a parasitic effect that is much less than the electroosmotic flow. The behavior of the curves in Figure 7 demonstrates a strong influence of the cross-section geometry. All the curves pass through the minimum given by (57) and (58). The leading term on the right-hand side of eq 58 is proportional to the combination of the geometrical factors, KE - K2×/KP, which, for different geometries, can take different values (Table 1). Accordingly, for such geometries, the minimums of (D*/D - 1)/Pe2 substantially differ from each other. For example, (D*/D - 1)min/Pe2 predicted for an ellipse with the aspect ratio λ ) 0.1 exceeds that obtained for two parallel planes by nearly a factor of 4. For an ellipse with aspect ratio λ ) 0.01, the minimum exceeds the value obtained for the parallel planes by two orders. At the same time, both (KE - K2×/KP) and the minimum nearly coincide for a rectangle having the aspect ratio λ f 0 and for a circle (ellipse with aspect ratio λ ) 1). Note, that, in the case of a rectangle, the geometrical factors are either constant or nearly constant (Figure 6b). Therefore, we discuss here only the example of a rectangle with aspect ratio λ f 0. The curves plotted for an annulus (γ ) 10) and for parallel planes pass very close to the points corresponding to the above-mentioned minimums. Although, at τ ≈ 0.04 (Figure 7a) and τ ≈ 0.02 (Figure 7b), the

Figure 7. Normalized variation of dispersion coefficient, (D*/D 1)/Pe2, as a function of τ for ζ˜ , 1 and for different cross-section geometries: (a) κh ) 50; (b) κh ) 100.

four above-mentioned geometries are characterized by close values of (D*/D - 1)/Pe2, at other τ values, the different geometries result in substantially different values of (D*/D - 1)/ Pe2. As observed from eq 58 and in Figure 7, with a decrease in κh, the minimums for the curves become weaker. Such behavior is a manifestation of a general trend: with a decrease in κh, the minimums should gradually disappear, and at κh f 0, dispersion becomes independent of the pressure-driven flow contribution, τ. For the latter limiting case, κh f 0, the electric charge density and bulk force become uniform within the cross section. Consequently, for such a regime, the electroosmotic velocity distribution is similar to that for pressure-driven flow. Due to such a similarity, changes in τ do not affect the dispersion. According to (58), the behavior of (D*/D - 1)min/Pe2 as a function of the electrokinetic potential, ζ˜ , is defined by the properties of the function g(ζ˜ ) (eq 9). Since, the leading term on the right-hand side of (58) is simply proportional to g(ζ˜ ), the behavior of (D*/D - 1)min/Pe2 as function of ζ˜ will be qualitatively similar to the behavior of g(ζ˜ ) itself. In details, this function was analyzed in our previous paper26 where it was shown that the Analytical Chemistry, Vol. 76, No. 10, May 15, 2004

2717

behavior of g(ζ˜ ) strongly depends on the electrolyte type. For example, in the case of a 1:2 electrolyte with double charged coion, when ζ˜ increases from 0 to 2, one can expect that both g(ζ˜ ) and (D*/D - 1)min/Pe2 will decrease by a factor of ∼2. Under the same condition, but for a 2:1 electrolyte with double charged counterion, the functions are nearly constant. The curves illustrating the behavior of (D*/D - 1)min/Pe2 as function of ζ˜ are shown in Supporting Information (Figures SI1-SI3). CONCLUSIONS Using the approach of the previous study,26 we considered the hydrodynamic dispersion in a combined pressure-driven and electroosmotic flow. Being valid for microchannels with a thin double layer, the derived expression 26 addresses the dispersion coefficient for an arbitrary cross-section geometry, electrolyte system, and electrokinetic potential. According to (26), the dispersion coefficient depends on the following geometryindependent parameters: (i) Peclet number given by (27); (ii) the parameter τ, which characterizes the contribution of the Poiseuille flow and is defined by (28); and (iii) the function of electrokinetic potential g(ζ), which also depends on an electrolyte type and is defined by (9). In eq 26, the influence of the crosssection geometry manifests itself through the channel hydrodynamic radius, h, and three dimensionless geometrical factors: KP, KE, and K×. The procedures for obtaining the geometrical factors were presented for arbitrary cross-section geometry. The geometrical factors KP, KE, and K× have been evaluated for several examples of cross-section geometry: (i) unbounded parallel planes; (ii) circle; (iii) annulus; (iv) ellipse; (v) rectangle. For some geometries, values of KP, KE, and K× are collected in Table 1. According to the obtained results, the value of the dispersion coefficient is strongly affected by the cross-section geometry, electrolyte type, and electrokinetic potential of the wall. The derived expressions predict variation of the dispersion coefficients within a range of several orders. It was shown that dispersion in a combined flow depends on the parameter τ (contribution of the Poiseuille flow) with a deep minimum that is observed when a weak Poiseuille flow is directed against the electroosmotic flow (τ < 0 and |τ| , 1). The minimum dispersion coefficient takes a substantially smaller value than that for the case of purely Poiseuille or electroosmotic flow. ACKNOWLEDGMENT The authors gratefully acknowledge financial support from the Natural Sciences and Engineering Research Council of Canada.

D

molecular diffusion coefficient

D*

hydrodynamic dispersion coefficient

dP,E

length-scale parameters determined applying the Taylor-Aris procedure, respectively, to the sole Poiseuille and electroosmotic flow

E

longitudinal electric field strength magnitude

F

Faraday constant

f(ζ˜ ), g(ζ˜ )

functions of the electrokinetic potential

h

hydrodynamic radius (h ) S/l)

bıx,y

unity vectors

KP,

KE,

K×

geometrical factors describing combined flow

L

channel length

l

cross-section perimeter

b n

the unity vector normal to the channel wall

R

gas constant

br

radius-vector

S

cross-section area absolute temperature

T u,

uP,

uE

local velocity of, respectively, the total, Poiseuille and electroosmotic flows

x, y

2D Cartesian coordinates

Zi

electric charge of the ith ion expressed in the Faraday units;

Greek Letters 

dielectric permittivity

Φ

normalized solution of the Poisson-Boltzmann problem for semi-infinite region

γ

ratio of external to internal radius of annulus (γ ) a1/a2) liquid viscosity

η κ,

κ-1

Debye parameter and Debye length

λ

aspect ratio for elliptic and rectangular cross sections

F

normalized polar radius

τ

contribution of the Poiseuille into combined flow

ω, W

functions determined within framework of the Taylor-Aris procedure

ξ

normalized coordinate

Ψ

equilibrium electric potential

ζ

electrokinetic potential

Miscellaneous SUPPORTING INFORMATION AVAILABLE Additional information as noted in text. This material is available free of charge via the Internet at http://pubs.acs.org.

∇ B ∇2

2D gradient operator )∇ B ‚∇ B

2D Laplace operator

〈...〉

average over cross-section area

〈...〉l

average over cross-section perimeter line

a

GLOSSARY radius of a circle; an ellipse axis; a side of rectangle

a1,2

external and internal radiuses of annulus

b

an ellipse axis; a side of rectangle

Received for review August 28, 2003. Accepted January 29, 2004.

Ci

concentration of the ith ion in equilibrium electroneutral solution

AC0303160

2718

Analytical Chemistry, Vol. 76, No. 10, May 15, 2004