Voltammetric Monitoring of Transient Hydrodynamic Flow Profiles in

Adnane Kara , Arnaud Reitz , Jessy Mathault , Syllia Mehou-Loko , Mehran ... Michelle Rogers , Chi Leong , Xize Niu , Andrew de Mello , Kim H. Parker ...
0 downloads 0 Views 374KB Size
Download PDF
Anal. Chem. 2007, 79, 626-631

Voltammetric Monitoring of Transient Hydrodynamic Flow Profiles in Microfluidic Flow Cells Mary Thompson and Richard G. Compton*

Physical and Theoretical Chemistry Laboratory, Oxford University, South Parks Road, Oxford OX1 3QZ, United Kingdom

We consider the transition to steady-state flow in the inlet region of a hydrodynamic channel cell and show that a microelectrode positioned within this inlet region allows chronoamperometric results to be recorded, from which information about the extent of the development of the flow profile may be deduced as well as information about the precise dimensions of the microfluidic channel. The use of microfluidics in the areas of analysis and synthesis has provoked great interest in recent years. Much work has been reported on the manufacture of microfluidic devices and cells1-3 as well as on their numerous uses. These include diverse applications in the fields of biological and environmental analysis,4-7 flow characterization,8-12 and sensing devices for organic and inorganic analysis,13-17 as well as uses in chemical and biological synthesis.18-22 A large proportion of these applications involve * To whom correspondence should be addressed. Fax: +44 (0) 1865 275410. Phone: +44 (0) 1865 275413. E-mail: [email protected]. (1) Daniel, D.; Gutz, I. G. R. Talanta 2005, 68, 429. (2) Mitrovski, S. M.; Nuzzo, R. G. Lab Chip 2005, 5, 634. (3) Pumera, M.; Merkoci, A.; Alegret, S. Trends Anal. Chem. 2006, 25, 219. (4) Wheeler, A. R.; Moon, H.; Kim, C.-J.; Loo, J. A.; Garrell, R. L. Anal. Chem. 2004, 76, 4833. (5) Kwakye, S.; Goral, V. N.; Baeumner, A. J. Biosens. Bioelectron. 2006, 21, 2217. (6) Marle, L.; Greenway, G. M. Trends Anal. Chem. 2005, 24, 795. (7) Morier, P.; Vollet, C.; Michel, P. E.; Reymond, F.; Rossier, J. S. Electrophoresis 2004, 25, 3761. (8) Amatore, C.; Belotti, M.; Chen, Y.; Roy, E.; Sella, C.; Thouin, L. J. Electroanal. Chem. 2004, 573, 333. (9) Bayraktar, T.; Pidugu, S. B. Int. J. Heat Mass Transfer 2006, 49, 815. (10) Amatore, C.; Oleinick, A.; Klymenko, O. V.; Svir, I. ChemPhysChem 2005, 6, 1581. (11) Amatore, C.; Klymenko, O. V.; Svir, I. ChemPhysChem 2006, 7, 482. (12) Amatore, C.; Oleinick, A.; Svir, I. Electrochem. Commun. 2004, 6, 1123. (13) Clegg, A. D.; Rees, N. V.; Klymenko, O. V.;Coles, B. A.; Compton, R. G. J. Am. Chem. Soc. 2004, 126, 6185. (14) Dayon, L.; Josserand, J.; Girault, H. H. Phys. Chem. Chem. Phys. 2005, 7, 4054. (15) Goral, V. N.; Zaytseva, N. V.; Baeumner, A. J. Lab Chip 2006, 6, 414. (16) Wang, J.; Escarpa, A.; Pumera, M.; Feldman, J. J. Chromatogr., A 2002, 952, 249. (17) Evenhuis, C. J.; Guijt, R. M.; Macka, M. Haddad, P. R. Electrophoresis 2004, 25, 3602. (18) Lee, S. H.; Lee, H. J.; Oh, D.; Lee, S. W.; Goto, H.; Buckmaster, R.; Yasukawa, T.; Matsue, T.; Hong, S.-K.; Ko, H.; Cho, M.-W.; Yao, T. J. Phys. Chem. B 2006, 110, 3856. (19) Baxendale, I. R.; Deeley, J.;Griffiths-Jones, C. M.; Ley, S. V.; Saaby, S.; Tranmer, G. K. Chem. Commun. 2006, 2566. (20) Zhou, X.; Cai, S.; Hong, A.; You, Q.; Yu, P.; Sheng, N.; Srivannavit, O.; Muranjan, S.; Rouillard, J. M.; Xia, Y.;Zhang, X.; Xiang, Q.; Ganesh, R.; Zhu, Q.; Matejko, A.; Gulari, E.; Gao, X. Nucleic Acids Res. 2004, 32, 5409.

626 Analytical Chemistry, Vol. 79, No. 2, January 15, 2007

electrochemistry as the monitoring technique,1,2,5,8,15,16,23-28 a good example of which is the high-speed channel electrode29,30 which has been used for various analyses.13,31 The modeling and interpretation of microfluidic cell data as well as the optimization of the design of the cells themselves typically presumes a steady-state flow has developed between the inlet to the cell and the spatial zone of interest. Depending on the Reynolds number employed, the distance before such a steadystate flow is achieved may not be insignificant. For example, in a rectangular channel of width d much greater than its height 2h (see Figure 1a) so that the flow profile can be regarded as twodimensional, the flow pattern is characterized by the Reynolds number, Re:

Re ) hv0/ν

(1)

where v0 is the fluid velocity at the center of the channel and ν is the kinematic viscosity (cm2 s-1) of the fluid. A laminar flow results if the Reynolds number remains below ∼2000.32,33 A lead-in length, le, given by

le ) (1 + 0.08(Re))h

(2)

is required to establish the steady-state flow profile,34 which is (21) Kartalov, E. P.; Quake, S. R. Nucleic Acids Res. 2004, 32, 2873. (22) Hung, P. J.; Lee, P. J.; Sabounchi, P.; Aghdam, N.; Lin, R.; Lee, L. P. Lab Chip 2005, 5, 44. (23) Matthews, S. M.; Du, G. Q.; Fisher, A. C. J. Solid State Electrochem. 2006, 10, 817. (24) Yi, C.; Zhang, Q.; Li, C.-W.; Yang, J.; Zhao, J.; Yang, M. Anal. Bioanal. Chem. 2006, 384, 1259. (25) Zhan, W.; Alvarez, J.; Crooks, R. M. J. Am. Chem. Soc. 2002, 124, 13265. (26) Kim, S. K.; Lim, H.; Chung, T. D.; Kim, H. C. Sens. Actuators, B 2006, 115, 212. (27) Sullivan, S. P.; Johns, M. L.; Matthews, S. M.; Fisher, A. C. Electrochem. Commun. 2005, 7, 1323. (28) Henley, I.; Fisher, A. Electroanalysis 2005, 17, 255. (29) Rees, N. V.; Alden, J. A.; Dryfe, R. A. W.; Coles, B. A.; Compton, R. G. J. Phys. Chem. 1995, 99, 14813. (30) Rees, N. V.; Dryfe, R. A. W.; Cooper, J. A.; Coles, B. A.; Compton, R. G.; Davies, S. G.; McCarthy, T. D. J. Phys. Chem. 1995, 99, 7096. (31) Rees, N. V.; Klymenko, O. V.; Coles, B. A.; Compton, R. G. J. Electroanal. Chem. 2002, 534, 151. (32) Levich, V. G. Physicochemical Hydrodynamics; Prentice Hall: New York, 1962. (33) Schlichting, H. Boundary-layer Theory; McGraw-Hill: New York, 1979. (34) White, F. M. Viscous Fluid Flow; McGraw-Hill: New York, 1974. 10.1021/ac0612022 CCC: $37.00

© 2007 American Chemical Society Published on Web 12/14/2006

Figure 2. A schematic diagram of the developing parabolic flow in the inlet region of a channel flow cell.

Figure 1. (a) Channel flow cell, showing the dimensions and coordinate system. (b) Uniform grid used for solution of the concentration profiles via the ADI method.

parabolic in nature:

[

vx ) v0 1 -

]

(h - y)2 h2

(3)

where y is defined in Figure 1a and vx is the fluid velocity in the x direction. The subject of microfluidics is concerned with small values of h and often large values of v0. Accordingly, the distance le is typically not necessarily insignificant in normal experimental practice. The purpose of the work described in this paper is to demonstrate computationally the proof-of-principle that the extent to which steady-state flow has been realized can be conveniently determined experimentally by means of a microelectrode located in the wall of the flow cell. Specifically we consider the evolution of the parabolic flow profile following entry into a rectangular channel flow cell (Figure 2) and calculate the current transient on a microband electrode located in one wall of the cell resulting from a potential step from a value at which no current flows to one relating to the transport-controlled discharge of an electroactive species dissolved in the fluid.35 In experimental practice dissolved oxygen is likely to find use in this respect. Currenttime transients are reported for different electrode locations, and their analysis is shown to be a suitable measure of the velocity profile in the cell. Moreover, transient analysis is shown to permit the measurement of the channel geometry. THEORY In this work we are concerned with the development of the parabolic flow profile within the inlet region of a channel flow cell and the electrochemical response obtained at an electrode positioned flush with one wall of the cell within the inlet region, when an electroactive species is dissolved in the flowing solution. In this section we describe how the flow profile is modeled numerically and how in turn the mass transport equation describing the electroactive species is solved using this profile. The flow is in a stage of development in the inlet region of a channel cell so that an electrochemical response, unlike that obtained in a region of fully developed flow, will be dependent on (35) Bard, A. J.; Faulkner, L. R. Electrochemical Methods, Fundamentals and Applications; John Wiley and Sons Inc.: New York, 2001.

Figure 3. Developing parabolic flow in the inlet region of a channel flow cell, showing the flow velocity ω as a function of the dimensionless coordinates Y and the stretched X*.

the axial position of the electrode. It follows that it may be possible at least in principle to deduce the flow profile and hence the cell dimensions from analysis of electrochemical data. It is this application that is addressed in this work. Here we are concerned with the development of parabolic flow in the “inlet” region32,33 of a duct (see Figure 2). Above the inlet there may be a reservoir of bulk solution so that at the entrance to the duct the flow profile is planar and then gradually develops to form the fully developed or parabolic flow at a certain distance from the duct entrance.33,34 In most realistic situations, the length of the inlet region is such that to detect a significant difference in response from that in a fully developed flow region, electrodes of microdimensions must be employed. Therefore, it is also necessary to consider the effect of axial diffusion in this region, and as such the model developed and described below uses the full mass transport equation relevant to the type of duct considered. For the channel electrode the mass transport equation used is

(

)

∂c ∂2c ∂2c ∂c ) D 2 + 2 - vx ∂t ∂x ∂x ∂y

(4)

where D is the diffusion coefficient of species c, vx is the axial flow velocity, and x and y are the coordinates in the axial direction and the direction of the channel height, respectively. Numerical Description of the Flow Profile. Modeling of the flow profile, vx, as a function of x and y was achieved using the work of Sparrow et al.36 This work on flow development in (36) Sparrow, E. M.; Lin, S. H.; Lundgren, T. S. Phys. Fluids 1964, 7, 338.

Analytical Chemistry, Vol. 79, No. 2, January 15, 2007

627

Figure 4. Concentration profiles showing species B at the steady state at an electrode positioned within the inlet region of a channel cell: (a) xgap ) 0 cm, (b) xgap ) 0.002 cm, and (c) xgap ) 0.005 cm. Other parameters: xe ) 0.001 cm, ν ) 1 × 10-3 cm2 s-1, d ) 0.3 cm, h ) 0.1 cm, D ) 1 × 10-5 cm2 s-1, Vf ) 5 × 10-4 cm3 s-1. The X* axis is nonuniform to allow easier observation of the concentration profiles obtained.

Figure 5. Concentration profiles showing species B at the steady state at an electrode positioned within the inlet region of a channel cell: (a) Vf ) 1 × 10-2 cm3 s-1, (b) Vf ) 5 × 10-3 cm3 s-1, and (c) Vf ) 5 × 10-4 cm3 s-1. Other parameters: xe ) 0.001 cm, ν ) 1 × 10-3 cm2 s-1, d ) 0.3 cm, h ) 0.1 cm, D ) 1 × 10-5 cm2 s-1, xgap ) 0.002 cm. The X* axis is nonuniform to allow easier observation of the concentration profiles obtained. Table 1. Cell Dimensions and Other Parameters Used in the Simulations param

value

param

value

xe d Vf

0.001 cm 0.2 cm 1 × 10-3 cm3 s-1

v D

1 × 10-3 cm2 s-1 1 × 10-5 cm2 s-1

entrance regions provides expressions for the flow velocity as a function of the duct cross-section coordinates at any value of the stretched axial coordinate, denoted by X*. Here X* is given by

X*/h X* ) Uh/v

u (5)

and X* is obtained from the relationship dx )  dX*, where  is the “stretching factor”.36 The expression given by Sparrow et al. can be evaluated numerically, and the stretched axial coordinate can be converted to the physical coordinate using the expression for the stretching factor , where dX )  dX*. The method of Sparrow et al. is in excellent agreement with other analyses, providing a much better approximation than the so-called “Targ-type” solution36 where  is assumed to be unity such that X ) X* throughout the duct length. The full derivation of the flow velocity solution can be found in Sparrow’s paper,36 with only a brief description given here, together with the specific application to the parallel-plate channel. Simulation of the chronoamperometric response of an electroactive species dissolved in the fluid flowing through the cell 628

Analytical Chemistry, Vol. 79, No. 2, January 15, 2007

was carried out in several stages. First the flow velocity profile was evaluated numerically for the region of interest. Then the relevant convection-diffusion equation was solved using the calculated flow velocities to obtain a concentration profile. The axial coordinate was converted from X* either to the physical dimensionless coordinate X or to the real coordinate x (cm), and finally the current response required was evaluated using the concentration profile obtained. Sparrow’s work derives the velocity solution for the parallelplate channel to be

U

∞

) ω ) 1.5(1 - Y2) +

∑

{

2 cos RiY

i)1 R i

2

cos Ri

}

- 1 e-Ri X* 2

(6)

where X* ) (x*/h)/(Uh/ν) and Y ) y/h. This equation corresponds to a sum of the fully developed flow velocity solution, ufd (first term), and an extra term denoted by u* (second term). In this result the first 25 values of Ri are listed in Sparrow’s paper36 and are used to reconstruct the flow velocity profile in this work, and ω is the flow velocity normalized with respect to the average flow velocity, U. For a channel cell, U is given by U ) Vf/2hd, where Vf is the volume flow rate and h and d are the channel half-height and width, respectively, as shown in Figure 1a. Figure 3 shows the development of the velocity profile within the inlet region of a channel, with ω shown as a function of Y and X*, from the inlet at X* ) 0 to X* ) 0.2, where the flow is parabolic. It is seen that initially there is a planar portion to the

Figure 6. Graph showing the values of Nulim as a function of xgap for different values of h. Other parameters: xe ) 0.001 cm, ν ) 1 × 10-3 cm2 s-1, d ) 0.2 cm, D ) 1 × 10-5 cm3 s-1.

flow profile in the center of the channel height coordinate, Y, which gradually diminishes along the X* coordinate to give a maximum velocity in the center, until the fully developed parabolic flow profile is obtained. Conversion to the physical dimensionless coordinate X (or the real coordinate x) can be computed by numerical evaluation of Sparrow’s expression for the stretching factor, , where dX )  dX*:

∂ω ∫ (2ω - 1.5ω )(∂X* ) ∂Y ) ∂ω + ∫ ( ) ∂Y (∂ω ∂Y ) ∂Y 1

2

0

1

Y)1

(7)

2

0

Electrochemical Response. For simplicity the mass transport equation was normalized in the same way as required for the solution of the velocity profile using the method of Sparrow et al.:36

∂c ∂2c ∂2c ∂c ) 2 + p1(i) - p2(i,j) ∂τ ∂Y ∂X* ∂X*2

(8)

p1(i) ) v2/U2h2i2

(9)

p2(i,j) ) vωi,j/DUi

(10)

where

Figure 7. Transient results for i/ilim as a function of time (s) for a range of xgap (cm) values: (a) h ) 0.1 cm, (b) h ) 0.25 cm, and (c) h ) 0.4 cm. The discontinuity seen in (c) for xgap ) 0 cm is an artifact arising from the magnification of the data obtained to highlight the difference between curves for different values of xgap.

for the half step implicit in Y, where

λy )

Equation 8 is discretized and solved using the alternating direction implicit (ADI) method over a uniform grid, so the discretized equations for solution are

λx1(i) )

k+1 k+1 k+1 ci-1,j (-λx1(i) - λx2(i,j)) + ci,j (1 + 2λx1(i) + λx2(i,j)) + ci+1,j

λx2(i,j) )

k k k (-λx1(i)) ) -ci,j-1 (λy) + ci,j (1 - 2λy) + ci,j+1 (λy) (11)

for the half time step implicit in X* and k k k ci-1,j (λx1(i) + λx2(i,j)) + ci,j (1 - 2λx1(i) - λx2(i,j)) + ci+1,j k+1 k+1 k+1 (λx1(i)) ) ci,j-1 (-λy) + ci,j (1 + 2λy) + ci,j+1 (-λy) (12)

∆τ/2 (∆Y)2

(∆τ/2)p1(i) (∆X*)2 (∆τ/2)p2(i,j) ∆X*

(13)

(14)

(15)

These sets of linear equations are then solved using the Thomas algorithm.37 With reference to Figure 1b, the normalized current, Nu, obtained at the electrode can be computed from the concentration (37) Fisher, A. C.; Compton, R. G. J. Phys. Chem. 1991, 95, 7538.

Analytical Chemistry, Vol. 79, No. 2, January 15, 2007

629

Figure 8. Contour plots showing the time taken (s) to reach nilim as a function of h (cm) and xgap (cm) for (a) n ) 2 and (b) n ) 1.2. Table 2. Pairs of Values for h and xgap for Which the Time To Reach nilim Matches the Experimental Dataa h xgap

h xgap

0.15 0.0061

0.20 0.0030

2ilim 0.25 0.0021

0.30 0.0019

0.35 0.0014

0.40 0.0012

0.15 0.0066

0.20 0.0030

1.2ilim 0.25 0.0019

0.30 0.0014

0.35 0.0010

0.40 0.0009

a Comparison of the two sets of data reveals one pair (shown in bold font) satisfying both contour plots.

profile obtained at each time step using Nxd c

Nu )

∑ Nxu

- ci,Ny ∆X*i ∆Y

i,Ny-1

(16)

Computational Details. Programs were all written using Delphi Pascal and were executed on a PC with a 2 GHz Pentium 4 processor and 1 GB of RAM. A typical simulation from X* ) 0 to X* ) 0.2 required a uniform grid of 4000 × 2000 (Ny × Nx) nodes for convergence and would run in several minutes. RESULTS AND DISCUSSION The results presented in this section contain some general results such as the analysis of concentration profiles and limiting current responses and the effect of varying the various parameters such as the flow rate and the cell and electrode dimensions. We then present a potentially useful practical method, which by using the results obtained from the simulations to provide proof-ofprinciple allows accurate measurement of the cell dimensions of a microfluidic channel electrode using experimental chronoamperometric data. Concentration Profiles in the Inlet Region. Figure 4 shows the concentration profiles once the steady-state response is achieved for various values of xgap, the distance from the inlet of the channel to the upstream edge of the electrode (see Figure 1a). The profiles show the concentration of species B at the electrode once the steady state is reached, in a potential step experiment where A ( e- f B such that the boundary condition at the electrode surface is a ) 0. We observe that, as xgap increases 630

Analytical Chemistry, Vol. 79, No. 2, January 15, 2007

from 0 cm in Figure 4a to 0.005 cm in Figure 4c, so too does the extent of the diffusion layer of B, which is consistent with the decrease in convective flux close to the electrode surface as we move along the axial coordinate. Figure 5 depicts the steady-state concentration profiles for a range of flow rates, Vf. Here we see that a decreasing flow rate means a growth in the diffusion layer and a greater influence of axial diffusion, exactly as would be expected. Effect of Varying Channel Dimensions. To see more clearly the effect of variation in the real parameters, we present results with some variables fixed as given in Table 1 and vary the real dimensions xgap and h. Sparrow’s results show that the full parabolic flow is achieved when X* ≈ 0.2, and from eq 5 we see that the value of the stretched real coordinate x* for which X* ) 0.2, at fixed d and ν, increases as h increases. A similar effect is observed for increasing Vf. As expected from the extension of x in the nonparabolic regime for increasing h, we find that steady state ilim values at higher h show a more steady decrease, while those at smaller h level off very quickly to the value expected for parabolic flow (Figure 6). The difference is observed more clearly if we examine the transient responses of i/ilim against t for each h value used; these graphs are shown in Figure 7, where we see that the transient responses become more distinct from one another as h increases. Evaluating Cell Dimensions from Chronoamperometry. The results presented above mean that, in practice, accurate information could be gained about the precise channel and electrode dimensions by analyzing simple chronoamperometric data. Here we demonstrate this procedure for the case where the unknowns of interest are the channel half-height, h, and the length of the channel upstream of the electrode, xgap. However, we note that since the chronoamperometry result depends on two dimensionless variables, p1 and p2, the analysis procedure could be used to determine other unknown real parameters, notably the flow rate Vf or the diffusion coefficient D. The chronoamperometry simulations can be analyzed by finding the time taken for the current to reach a value of nilim, where ilim is the steady-state limiting current. Simulations were carried out using the real parameters as given in Table 1 and for a range of xgap and h values. Thus, surfaces can be drawn for the time in seconds to reach nilim as a function of xgap and h (Figure 8). Subsequent analysis of simulated “experimental” data (for

which h and xgap are unknown) would then yield a set of pairs of xgap-h values that would be consistent with such a result. Repeating the process for a different value of n then reveals a unique solution for the xgap-h pair that satisfies the time to reach nilim on both xgap-h surfaces. Alternatively, having evaluated the pairs of values from one surface, a best estimate of the unknown dimensions can be obtained by comparing experimental and simulated values for the absolute limiting current. To demonstrate this process, an “experimental” transient was simulated, using the values xgap ) 0.003 cm and h ) 0.2 cm. The values of the time taken to reach nilim were evaluated from the current-time transient obtained and were found to be 0.0064 s for n ) 2 and 0.027 s for n ) 1.2. The contour plots shown in Figure 8 were then used to find the pairs of xgap-h values that satisfied these conditions: the value of h was fixed at 0.15, 0.2, 0.25, and so on, and the value of xgap that gave the appropriate time was evaluated. This process was repeated for both n ) 2 and n ) 1.2, and the results are shown in Table 2. We see from the table that there is a unique pair of values that satisfies the time to reach nilim for both n ) 2 and n ) 1.2; thus, this method provides a suitable means for determining these cell dimensions from a single chronoamperometry experiment. CONCLUSIONS We have demonstrated that information on the dimensions of a channel or other parameters describing an electroactive species may be obtained by means of a simple chronoamperometric

experiment using a microelectrode positioned within the inlet region of the channel. We have reconstructed the flow profile numerically using the analysis of Sparrow et al.36 to give an accurate description of the flow velocities as the flow changes from a linear profile to a fully developed parabolic profile in the inlet length. We then used the flow profile obtained in the mass transport equation describing the electroactive species within this region; solving the mass transport equation allowed us to simulate the chronoamperometric response of an electroactive species dissolved in the fluid. We have noted that the normalized form of the diffusionconvection equation gives results dependent on the two dimensionless parameters p1 and p2. To give a useful, practical tool for experimentalists, we have demonstrated how the simulated results could be used to elucidate precisely the two unknown dimensions xgap and h. This method of analysis could provide a basis for accurate measurements in microfluidic channels, particularly for determining precise cell dimensions not immediately obtainable by other methods. ACKNOWLEDGMENT M.T. thanks the EPSRC for financial support.

Received for review July 3, 2006. Accepted November 8, 2006. AC0612022

Analytical Chemistry, Vol. 79, No. 2, January 15, 2007

631