Anal. Chem. 1998, 70, 4494-4504
Computer Simulations of Electrokinetic Transport in Microfabricated Channel Structures Sergey V. Ermakov, Stephen C. Jacobson, and J. Michael Ramsey*
Chemical and Analytical Sciences Division, Oak Ridge National Laboratory, P.O. Box 2008, Oak Ridge, Tennessee 37831-6142
A mathematical model describing electrokinetically driven mass transport phenomena in microfabricated chip devices is presented in this paper. The model accounts for principal material transport mechanisms such as electrokinetic migration (electrophoresis and electroosmosis) and diffusion. A computer code that implements the model is capable of simulating transport of materials during electrokinetic manipulation in 2-D channel structures. The computer code allows arbitrary channel geometries with various boundary conditions for the electric field and the sample concentration. Two fundamental microfluidic chip elements, a cross and a mixing tee, are of particular interest. An electrokinetic focusing experiment using a cross structure and mixing in a tee structure are simulated. Simulations revealed an optimum focusing voltage for which the ratio of sample concentration to sample width is maximized. They also verified that the mixing tee provides very accurate dilution/mixing characteristics for both charged and neutral samples. Good agreement between simulated and experimental data verified the accuracy of the mathematical model. Microfabricated fluidic (microchip) devices that perform chemical and biochemical analysis procedures have attracted much attention over the past few years. Early demonstrations of these devices involved chemical separations of various species including the separation of small ions in free solution by electrophoresis,1-8 capillary gel electrophoresis9-12 for the separation of nucleic acids, electrophoretic separation of proteins,13 and (1) Harrison, D. J.; Manz, A.; Fan, Z.; Lu ¨ di, H.; Widmer, H. M. Anal. Chem. 1992, 64, 1926. (2) Manz, A.; Harrison, D. J.; Verpoorte, E. M. J.; Fettinger, J. C.; Paulus, A.; Lu ¨ di, H.; Widmer, H. M. J. Chromatogr. 1992, 593, 253. (3) Effenhauser, C. S.; Manz, A.; Widmer, H. M. Anal. Chem. 1993, 65, 2637. (4) Seiler, K.; Harrison, D. J.; Manz, A. Anal. Chem. 1993, 65, 1481. (5) Harrison, D. J.; Fluri, K.; Seiler, K.; Fan, Z.; Effenhauser, C. S.; Manz, A. Science 1993, 261, 895. (6) Jacobson, S. C.; Hergenro¨der, R.; Koutny, L. B.; Warmack, R. J.; Ramsey, J. M. Anal. Chem. 1994, 66, 1107. (7) Jacobson, S. C.; Hergenro¨der, R.; Koutny, L. B.; Ramsey, J. M. Anal. Chem. 1994, 66, 1114. (8) Koutny, L. B.; Schmalzing, D.; Taylor, T. M.; Fuchs, M. Anal. Chem. 1996, 68, 18. (9) Ekstrom, B.; Jacobson, G.; Ohman, O.; Sjodin, H. International Patent WO 91/16966, 1990; Adv. Chromatogr. 1993, 33, 1. (10) Effenhauser, C. S.; Paulus, A.; Manz, A.; Widmer, H. M. Anal. Chem. 1994, 66, 2949. (11) Woolley, A. T.; Mathies, R. A. Proc. Natl. Acad. Sci. U.S.A. 1994, 91, 11348. (12) Woolley, A. T.; Mathies, R. A. Anal. Chem. 1995, 67, 3676.
4494 Analytical Chemistry, Vol. 70, No. 21, November 1, 1998
electrokinetically driven chromatographic separations of neutral species.14-16 The separative performance produced by these devices was similar to that of conventional capillary experiments, but results were obtained in time periods that are typically 2 orders of magnitude shorter using injection volumes in the nanoliter to picoliter range. The short analysis times produced by these devices are largely the result of the integrated electrokinetic valves which allow injection of small axial extent plugs of sample for separation that allows the use of separation distances of a few centimeters. The valves are actuated by controlling the electric potentials at the four ports of a “cross” channel structure. A potentially powerful extension of these microchip separation devices is the Lab-on-a-Chip concept, where sample processing and chemical reactions are integrated with chemical separations using electrokinetic material transport throughout the structures. Samples and reagents are guided through an appropriately designed microfluidic network in order to carry out an experimental protocol that includes mixing of reagents, control of reaction times, and injection and separation in an analysis channel. Demonstrated devices that monolithically integrated chemical and biochemical reactions with separations include derivatization reactions before17 and after18,19 electrophoretic separation and restriction enzyme digestion of nucleic acids followed by electrophoretic sizing.20 These latter devices include five ports which require electric potential control for precise fluidic manipulation and represent an advancement in measurement and fluidic complexity. Five-port microfluidic structures have also been recently demonstrated for the implementation of solvent-programmed micellar electrokinetic capillary chromatography21 and open-channel electrochromatography.22 Examples of more complicated microchips that include six or more ports include (13) Colyer, C. L.; Mangru, S. D.; Harrison, D. J. J. Chromatogr. A 1997, 781, 271. (14) Jacobson, S. C.; Hergenro¨der, R.; Koutny, L. B.; Ramsey, J. M. Anal. Chem. 1994, 66, 2369. (15) Moore, A. W.; Jacobson, S. C.; Ramsey, J. M. Anal. Chem. 1995, 67, 4184. (16) von Heeren, F.; Verpoorte, E.; Manz, A.; Thormann, W. Anal. Chem. 1996, 68, 2044. (17) Jacobson, S. C.; Hergenro¨der, R.; Moore, A. W.; Ramsey, J. M. Anal. Chem. 1994, 66, 4127. (18) Jacobson, S. C.; Koutny, L. B.; Hergenro¨der, R.; Moore, A. W.; Ramsey, J. M. Anal. Chem. 1994, 66, 3472. (19) Fluri, K.; Fitzpatrick, G.; Chiem, N.; Harrison, D. J. Anal. Chem. 1996, 68, 4285. (20) Jacobson, S. C.; Ramsey, J. M. Anal. Chem. 1996, 68, 720. (21) Kutter, J. P.; Jacobson, S. C.; Ramsey, J. M. Anal. Chem. 1997, 69, 5165. (22) Kutter, J. P.; Jacobson, S. C.; Matsubara, N.; Ramsey, J. M. Anal. Chem. 1998, 70, 3291. 10.1021/ac980551w CCC: $15.00
© 1998 American Chemical Society Published on Web 10/02/1998
structures designed to study enzyme kinetics,23-24 synchronized cyclic electrophoresis,25 and an integrated device for competitive immunoassay experiments.26 Lab-on-a-Chip devices for solving increasingly complex chemical and biochemical problems will require correspondingly complex fluidic designs and electrokinetic control strategies. Rather than relying solely on experimental means to design and test microfluidic structures, it seems prudent to develop computeraided design and analysis tools for prototyping and refining fluidic layouts. Ultimately these design and analysis tools should include the ability to model material transport and chemical and biochemical processes within microchips. Such tools will clearly reduce the time and expense for development of microchip devices designed to solve specific measurement problems. This paper describes our initial efforts to develop such design tools for electrokinetically driven microchip devices. Electrophoretic transport has been intensively investigated both theoretically and experimentally, and mathematical models have been developed to describe this phenomenon in detail27-29 (earlier theoretical works are reflected in the list of references in these books). These models are based on fundamental conservation laws expressed in a set of partial differential equations which describe electrophoresis to varying extents of complexity and detail. Attempts to resolve the equations for arbitrary initial and boundary conditions using analytical methods encounter enormous difficulties. Only a narrow class of problems can be solved in closed form (see, for example, refs 27 and 29). Most successful theoretical studies of sample transport in electrophoresis have been conducted using computer modeling which was initiated by Bier’s group.30 Computers have been applied to the modeling of most common electrophoretic techniques and classes of samples.31-34 More recently, computer modeling has been used to simulate capillary electrophoresis (CE).35-43 Both qualitative and quantitative agreement between theory and experiment has been (23) Hadd, A. G.; Raymond, D. E.; Halliwell, J. W.; Jacobson, S. C.; Ramsey, J. M. Anal. Chem. 1997, 69, 3407. (24) Hadd, A. G.; Jacobson, S. C.; Ramsey, J. M. Anal. Chem., to be submitted. (25) Burggraf, N.; Manz, A.; Effenhauser, C. S.; Verpoorte, E.; de Rooij, N. F.; Widmer, H. M. J. High Resolut. Chromatogr. 1993, 16, 594. (26) Chiem, N. H.; Harrison, D. J. Clin. Chem. 1998, 44, 591. (27) Babskii, V. G.; Zhukov, M. Yu.; Yudovich, V. I. Mathematical Theory of Electrophoresis; Naukova Dumka: Kiev, 1983 (in Russian; English translation by Consultants Bureau: New York, 1989). (28) Mosher, R. A.; Saville, D. A.; Thormann, W. The Dynamics of Electrophoresis; VCH: Weinheim, 1992. (29) Everaerts, F. M.; Beckers, J. L.; Verheggen, T. P. E. M. Isotachophoresis: Theory, Instrumentation and Application; Elsevier: Amsterdam, 1976. (30) Bier, M.; Palusinski, O. A.; Mosher, R. A.; Saville, D. A. Science 1983, 219, 1281. (31) Thormann, W.; Mosher, R. A. Adv. Electrophor. 1988, 2, 45. (32) Mosher, R. A.; Dewey, D.; Thormann, W.; Saville, D. A.; Bier, M. Anal. Chem. 1989, 61, 362. (33) Mosher, R. A.; Gebauer, P.; Caslavska, J.; Thormann, W. Anal. Chem. 1992, 64, 2991. (34) Mosher, R. A.; Gebauer, P.; Thormann, W. J. Chromatogr. A 1993, 638, 155. (35) Roberts G. O.; Rhodes P. H.; Snyder, R. S. J. Chromatogr. 1989, 480, 35. (36) Ermakov, S. V.; Righetti, P. G., J. Chromatogr. A 1994, 667, 257. (37) Ermakov, S. V.; Zhukov, M. Yu.; Capelli, L.; Righetti, P. G. Anal. Chem. 1994, 66, 4034. (38) Ermakov, S. V.; Zhukov, M. Yu.; Capelli, L.; Righetti, P. G. Anal. Chem. 1995, 67, 2957. (39) Ermakov, S. V.; Zhukov, M. Yu.; Capelli, L.; Righetti, P. G. Electrophoresis 1995, 16, 2149. (40) Ermakov, S. V.; Zhukov, M. Yu.; Capelli, L.; Righetti, P. G. Electrophoresis 1998, 19, 192.
obtained.36 Computer simulations have also been used to investigate and explain the unusual phenomena of sample peak splitting in CE,37,38 the existence of additional zones and peaks in isotachophoresis,39,40 and the impact of electroosmotic flow on CE separations.43 Most models are limited to 1-D cases. There are a few reports which solve 2-D and 3-D problems, such as modeling of a nonuniform electroosmotic flow due to nonuniform ζ potential in capillary electrophoresis44 and sample dispersion in free flow electrophoresis.45,46 These simulations were performed for simple separation column geometries (parallelepiped or cylinder). Electroosmotic flow in a more complicated geometry such as a rectangular channel intersection has been studied recently by Patankar and Hu.47 They made 3-D flow simulations of the focusing step which usually precedes the sample injection.6 This report studies sample mass transport in microchip electrophoresis using computer modeling. The 2-D mathematical model formulated below accounts for the principal physical phenomena affecting sample mass transport in microchip channels such as the electroosmotic flow, electrophoretic motion, and diffusion. Unlike ref 47, in which electroosmotic flow of the buffer was the target for simulations, we concentrate on the studies of the sample transport under various operating conditions. The simulations reproduce some features characteristic of sample manipulation on microchip, e.g., electrokinetic focusing and sample dilution. The simulation results are compared with experimental data. MATHEMATICAL MODEL Preliminary Considerations. Sample mass transport in electrophoresis is the result of three major mass transport mechanisms: convection, electrophoretic transport, and diffusion. Bulk convective flow of liquid arises from one or a combination of the following factors: (i) a pressure difference applied to the separation column ends, (ii) electroosmotic flow which has its origin at the column walls, and (iii) thermal convection due to Joule heating. Electrophoretic transport is determined by the electric field distribution, which depends on the boundary conditions and the conductivity of the solution at each point along the channel, which in turn depends on the chemical composition of the solution. Finally, diffusion takes place whenever a spatial nonuniformity in the composition of the solution exists. These three sample transport mechanisms depend on other physical phenomena accompanying electrophoresis: heat transport in the solution, chemical reactions within the buffer, and, of course, the electric field distribution. Mathematical models of electrophoresis formulated previously27-29 take into account all the phenomena mentioned above, but in many situations the contributions of (41) Schwer, C.; Gasˇ, B.; Lottspeich, F.; Kenndler, E. Anal. Chem. 1993, 65, 2108. (42) Schafer-Nielsen, C. Electrophoresis 1995, 16, 1369. (43) Thormann, W.; Zhang, C.-X.; Caslavska, J.; Gebauer, P.; Mosher, R. A. Anal. Chem. 1998, 70, 549. (44) Potocˇek, B.; Gasˇ, B.; Kenndler, E.; Sˇ teˇdry´, M. J. Chromatogr. A 1995, 709, 51. (45) Aksenov, A. A.; Gudzovsky, A. V.; Serebrov, A. A. Proc. 5th Int. Symp. On Computational Fluid Dynamics, Sendai, Vol. I, 1993 JSCF, 19. (46) Clifton, M. J.; Balmann Roux-de, H.; Sanchez, V. AIChE J. 1996, 42, 2069. (47) Patankar, N. A.; Hu, H. H. Anal. Chem. 1998, 70, 1870.
Analytical Chemistry, Vol. 70, No. 21, November 1, 1998
4495
different phenomena to the evolution of the electrophoretic system are not equally important. Basic Assumptions. First, transport processes have different time scales. Since we are interested in a macroscopic description of the electrophoretic system and its evolution as applied to microchip devices, we shall use the electrophoretic time scale as a reference point. Association-dissociation reactions for sample and buffer species, which in reality should also be considered as a transport mechanism, is the fastest process. Under normal circumstances, its time scale is many orders of magnitude less than that of other processes. Chemical composition of the solution immediately adjusts to any macroscopic changes in the system so that the acid-base equilibrium will exist locally at any given moment in time. Therefore, the equations of chemical kinetics can be substituted with equations of chemical equilibrium. In the simplest case, for a solution of univalent acids HAn and bases Bm, the dissociation scheme can be summarized by the following equations: + + HAn T H+ + An ; HBm T H + Bm
(1)
The equilibrium between the dissociated and the nondissociated forms is determined by the pH of the solution (pH ) -log[H+] ) and the equilibrium constants K na and K mb:
[A[H+Bm] n] , βm ) , an bm a Kn [H+] , β Rn ) a ) m Kn + [H+] Kbm + [H+] Rn )
(2)
+ where an ) [An ] + [HAn] and bm ) [B] + [H Bm] are the analytical concentrations of an acid and a base, Rn and βm are the corresponding dissociation degrees, and [H+] is the hydrogen ion concentration. We will further assume that the electric double layer arising at the chip wall-buffer solution interface is also adjusted instantly and that its primary characteristic, the local ζ potential, depends on the local pH value only. Convective transport and electrophoresis have similar time scales for cases in which convective flow is driven by electroosmosis, since electroosmotic and electrophoretic mobilities are usually of the same order. A characteristic time for electrophoresis can be defined as the time te that it takes for a particle to travel a distance equal to the channel width L, i.e., te ) L/Ve, where Ve is the electrophoretic velocity. The time scale for convection caused by pressure and temperature gradients depends on the particular conditions and can be either larger or smaller than the electrophoretic time. Diffusion is usually significantly slower than electrophoresis, and its characteristic time can be related to that of electrophoresis by the ratio Ped ) VeL/D, which sometimes is referred to as the diffusion Pe´clet number. For typical Ve , L, and D values, Ped lies between 100 and 5000, which means that diffusion is 100-5000 times slower than electrophoresis. Very often, however, the duration of the experiment is comparable with the transverse diffusion time scale, and, therefore, diffusion cannot be neglected. The next approximation concerns the calculation of the electric field. In an electrolyte solution, positively and negatively charged
4496 Analytical Chemistry, Vol. 70, No. 21, November 1, 1998
ions resulting from the reactions of dissociation give rise to the electric field which can be calculated using Poisson’s equation:
∇2φ )
1 F 0 e
(3)
where φ is the electric field potential, Fe is the electric charge density, is the dielectric constant of the solvent, and 0 ) 8.854 × 10-12 C V-1 m-1 is the dielectric permittivity of a vacuum. On the macroscopic level, the solution is electrically neutral. The electroneutrality is violated on scales comparable to the Debye radius rD around each ion:
rD )
x
kT0
∑
e2
(4)
nizi2
i
Here, k is Boltzmann’s constant, T is the absolute temperature, e is the electron charge, and ni is the concentration of the type i ions, while zi is the charge. Unless the solution is extremely dilute (concentration of electrolyte is less than 1 µM), the Debye radius is usually on the order of several nanometers; therefore, for “macroscopic” considerations with a length scale of tens of micrometers, the hypothesis of electroneutrality (Fe ) 0) is quite accurate for the bulk solution. Special attention has to be paid to the solution-wall interface, where the electrical double layer is formed. Unlike the bulk solution, in which ions are distributed quite randomly, here ions of the same charge are concentrated in close proximity to the wall having the same density of charge but of opposite sign. Such ion redistribution results in the existence of a wall potential relative to that of the solution far from the wall, known as the ζ potential. Although the thickness of the double layer is on the order of the Debye radius (see eq 4), its existence leads to electroosmotic flow when axial electric fields are present. The electroosmotic flow profile has been theoretically calculated for simple geometries, such as narrow capillary slits48 and cylindrical capillaries.49 It has been shown that it behaves very similarly to the potential profile φ(z) across the capillary or slit. The velocity component parallel to the wall has a flat profile; i.e., the velocity is constant almost everywhere, except in close proximity to the wall, where it drops sharply to zero. For an infinite rectangular channel and no external pressure applied, the following simple formula was obtained:48
(
V ) µeoE 1 -
)
φ(z) , ζ
φ(0) ) ζ, φ(z) f 0 at z . rD, µeo )
0ζ η
(5)
where E is the axial electric field, µeo is the electroosmotic mobility, η is the viscosity of the solution, and z is the coordinate normal to the channel wall. Thus, looking at the flow field on scales much larger than the double layer, it may be assumed that, along the wall, the liquid slips with a velocity V ) µeoE. This can be (48) Burgeen, D.; Nakache, F. R. J. Phys. Chem. 1964, 68, 1084. (49) Rice, C. L.; Whitehead, R. J. Phys. Chem. 1965, 69, 4017.
considered as the boundary condition for the equations of fluid motion. Below we will hypothesize that the above approximation is also valid in the case when the channels have finite length and a more complex geometry, e.g., when channels cross each other. The only requirement that remains is that the characteristic channel dimension must be much larger than the thickness of the double layer. In the case of the channel intersection, this hypothesis has been recently verified in the simulations made by Patankar and Hu.47 When dealing with high field strengths and current densities in chip electrophoresis, one should be aware of a potential increase in buffer and sample temperature and, as a result, increased axial dispersion. In capillary electrophoresis, the problem of Joule heating has been investigated in a number of papers.50-52 It has been found that the temperature profile is nearly parabolic for moderate applied power.51 At higher applied power, the temperature dependence of the buffer conductivity becomes important, and in this situation a parabolic profile underestimates the temperature gradient. Such positive feedback could eventually lead to a thermal runaway.51,52 Temperature inhomogeneities decrease the resolving power of CE experiments,50 and attempts to partially compensate for it by applying a Poiseuille flow to the capillary have been reported.52,53 Jansson et al.54 demonstrated that cylindrical microchannels (capillaries) are not as effective in heat dissipation because the ratio of the channel surface/channel volume is minimal. On the other hand, rectangular channels with a large aspect ratio can provide much better temperature stabilization. Such channel geometries are typical of microfabricated devices. Their advantage in maintaining a constant temperature for chip electrophoresis has been experimentally verified by Fan and Harrison.55 They showed that, for typical channel dimensions and low conductivity buffers, one can safely use field strengths up to 2.5 kV/cm, which exceeds the value tolerated in capillary electrophoresis by more than 2 times. Jacobson et al. successfully exploited even higher field strengths for ultrahigh-speed separations.56 They report that, for short separation times and specially designed chips, Joule heating did not significantly affect the resolution up to field strengths values of 50 kV/cm. In our model, we will assume that the running parameters do not exceed the critical values beyond which the temperature buildup in the channel due to Joule heating is more that 1-2 K; i.e., isothermal conditions are considered. An estimation of the temperature increase can be obtained from the following formula:54
∆T )
(
)
SσE2 π h + 8 λwall(ln 2 + (πw)/(4q)) Lλbuf
assuming that the channel outer wall has a fixed temperature. Here, S ) hL is the channel cross section area, q is the thickness of the channel wall, σ is the buffer conductivity, and λwall, λbuf are (50) Grushka, E.; McCormick, R. M.; Kirkland, J. J. Anal. Chem. 1989, 61, 241. (51) Jones, A. E.; Grushka, E. J. Chromatogr. 1989, 466, 219. (52) Gobie, W. A.; Ivory, C. F. J. Chromatogr. 1991, 516, 191. (53) Kutter, J.; Welsch, T. J. High Resolut. Chromatogr. 1995, 18, 741. (54) Jansson, M.; Emmer, A.; Roeraade, J. J. High Resolut. Chromatogr. 1989, 12, 797. (55) Fan, Z. H.; Harrison, D. J. Anal. Chem. 1994, 66, 177. (56) Jacobson, S. C.; Culbertson, C. T.; Daler, J. E.; Ramsey, J. M. Anal. Chem. 1998, 70, 3476.
the thermal conductivity of the substrate wall and buffer solution, respectively. For instance, for a 100 mM Tris-HCl buffer (pH 8.08) in a 10-µm × 50-µm channel with a quartz channel wall thickness of 100 µm, one can use field strengths (E) of up to 1.5 kV/cm while keeping the temperature difference between buffer and ambient below 1 K. When using less concentrated buffers, the field strength can be increased even more. Taking into account the above approximations, the set of equations describing the electric field distribution and sample mass transport is shown below:
σ)F
(
∑(z ) R a µ
a 2 a n n n n
+
n
∇×E B)0
(6)
∇‚Bj ) 0
(7)
Bj ) σE B
(8)
∑(z
b 2 b m) βmbmµm
+ µH[H+] +
m
)
µOH[OH-] (9)
∑z R a
a n n n
n
+
∑z
b mβmbm
+ [H+] - [OH-] ) 0
(10)
m
[H+][OH-] ) Kw
(11)
∂an + ∇((V B+B Vaep)an) ) Dan∇2an, B Vaep ) RnzanµanE B ∂t
(12)
∂bm + ∇((V B+B Vbep)bm) ) Dbm∇2bm, B Vbep ) βmzbmµbmE B ∂t
(13)
1 ∂V B V + (V B∇)V B ) - ∇p + ν∇2B ∂t F
(14)
∇‚V B)0
(15)
Here, ∇× and ∇‚ are the curl and divergence differential operators, and F is Faraday’s constant. Equation 6 states that the electric field is a potential (conservative) field, eq 7 stands for charge conservation, and eq 8 expresses Ohm’s law for the solution. Equation 9 describes the conductivity of the solution, assuming that the buffer and the sample compounds are monovalent acids and bases dissociating according to eqs 1 and 2. The current due to diffusion flux is supposed to be small, so we can exclude it from consideration. The electroneutrality condition is expressed by eq 10, which together with the ionic product of water (eq 11) allows the calculation of the concentration of hydrogen ions (pH value). Equations 12 and 13 describe buffer and sample component transport, while eq 14 is the Navier-Stokes equation describing bulk flow of the solution in conjunction with the continuity eq 15. Further Simplifications. The set of equations 2, 6-15 is rather complicated for immediate implementation in a computer code, and in many instances it can be further simplified while still preserving essential physical features. One such simplification is to decrease the dimensionality of the problem from 3-D to 2-D. We stipulate that all channels in the chip have a constant depth and their geometry depends only on x and y coordinates (Figure 1). Physicochemical properties of the chip substrate and cover Analytical Chemistry, Vol. 70, No. 21, November 1, 1998
4497
Figure 1. Elements of the microchip channel structure (cross) and the corresponding 2-D computational domain.
plate are assumed to be identical and homogeneous. Under these circumstances, the flow profile is assumed to be flat and independent of channel depth. 3-D simulation of the electroosmotic flow in the channel intersection performed by Patankar and Hu47 confirmed this. The next simplification concerns the calculation of the electric field and the sample and buffer composition. In the vast majority of experiments performed on microchips, zone electrophoresis is utilized as the separation mode. Usually during such experiments, the buffer concentration is many times (100 or more) higher than that of the sample species, which minimizes electromigration dispersion. In this situation, the conductivity of the solution is essentially constant across both buffer and sample zones (σ ) σbuffer ) const). A uniform buffer concentration assumes that the pH value is also uniform. With these assumptions, eqs 10 and 11 from which the concentration of H+ ion is calculated can be neglected. From eqs 12 and 13, only those describing the sample species should be considered since the buffer composition is constant. Equations 6-9 are reduced to one equation for the external electric potential Φ with the corresponding boundary conditions:
∇2Φ ) 0, E B ) -∇Φ
(16)
Electrophoretic and electroosmotic velocity will vary only with the local electric field strength, not with the pH value or concentration of the buffer/sample species. As a result, the flow field, electric field, and concentration field are “partially independent”. The electric field determines the flow field; together they determine the concentration distribution, but within the above assumptions there is no influence of the concentration field on the bulk flow 4498 Analytical Chemistry, Vol. 70, No. 21, November 1, 1998
or the electric field. The bulk flow in turn does not affect the electric field distribution. Since the electric field (eq 16) does not evolve in time, it needs to be calculated only once for given boundary conditions at uniform solution conductivity. Following the electric field, the bulk flow quickly reaches steady state. From hydrodynamics (see ref 57), it is known that one of the parameters determining the flow regime is the dimensionless Reynolds number, Re ) VL/ν, giving the ratio between inertial forces and viscous forces, where V is the velocity, L is the characteristic length, and ν is the kinematic viscosity. For a typical chip experiment, Re e 0.1 (V ≈ 10-1 cm/s, L ≈ 3 × 10-3 cm, ν ≈ 0.9 × 10-2 cm2/s); hence, viscous forces prevail and define a characteristic time scale at which the flow field reaches a steady state. The order of magnitude for the time scale can be estimated as tsteady ≈ L2/ν ≈ 10-3 s. Similar time scales have been calculated using more accurate simulations for cylindrical capillaries58 and 2-D straight channel.47 This time is usually small compared to all other characteristic times of sample manipulation, such as injection, separation, mixing, etc. Hence, when doing computer simulations with time steps similar or greater than tsteady, the fluid flow can be considered to be at steady-state conditions. This significantly reduces computational time since the flow needs to be calculated only when conditions change, such as the redistribution of electric potentials applied to reservoirs. Numerical Algorithm. The differential eqs 12-16 are solved using finite-difference algorithms similar to those described in ref 59. The simulation domain under consideration is covered by a rectangular grid (Figure 1), and then for each cell discrete versions of eqs 12-16 are derived. The most difficult problem involves solving the Navier-Stokes equations (eqs 14 and 15), since the discrete version of eq 15 cannot be used directly to find a solution. Instead, it must be used as a condition when solving another differential equation. The computational algorithm for which known velocity and pressure fields at time tn are used to calculate their values for the next moment, tn+1 ) tn + τ, is as follows:
1 V ˜ -Vn ˜ ) - ∇Fn + ν∇2V ˜ + (V n∇)V τ F ∇2(δp) )
∇V ˜ , δp ) pn+1 - pn τ
V n+1 - V ˜ ) -∇(δp) τ
(17) (18)
(19)
where the velocity V˜ is an intermediate value, and δp is the pressure increment in time. It is easy to check that a discrete version of eq 15 is satisfied in this algorithm once the solution for eqs 17-19 is found. By applying the divergence operator ∇‚ to both parts of eq 19 and taking into account eq 18, one gets ∇Vn+1/τ ) 0 or ∇Vn+1 ) 0. The discrete representation of eqs 12-16 results in a set of linear algebraic equations which are solved by the conjugate gradient method with incomplete Cholesky (57) Lamb, H. Hydrodynamics; Dover Publications: New York, 1945; p 738. (58) Dose, E.; Guiochon, G. J. Chromatogr. A 1993, 652, 263. (59) Fletcher, C. A. J. Computational techniques for fluid dynamics; SpringerVerlag: Berlin, 1988.
Figure 2. Experimental and simulated images of electrokinetic sample focusing. In each pair, the left image represents experimental data and the right image computational. Channel geometry as used in the experiment is indicated in panel a. (Reproduced with permission from ref 63.)
decomposition.60,61 These algorithms were implemented with inhouse written code. One of the major challenges was to make the code both flexible and computationally economic so that the channel geometry could be easily and interactively changed while keeping the computational overhead to a minimum. The computer code was tested by solving a simplified hydrodynamic problem: convective flow in a square cavity with a moving upper boundary, which is a classical, well-documented test case.62 The simulated results showed good agreement (less than 1% discrepancy in the values of the flow function) with the data from ref 62. More complex tests are limited by the lack of available analytical and numerical solutions. Computer simulations were performed for two basic chip elements: a cross used for sample focusing and injection and a channel tee used for sample mixing. The sample contained only one species. The simulation process contains several steps. First, it is necessary to design the channel structure using a simple drawing program. The cells belonging to the simulation domain are selected from a finite-difference grid. The set of such cells will form the desired channel geometry. Second, the boundary conditions for the fluid motion (V and p), the electric field (E and/ or Φ), and the concentration (c) at the channel walls and inlets/ outlets must be specified. Finally, other running parameters, such as analyte diffusion coefficient and electrophoretic mobility, and the electroosmotic mobility of the wall must be provided. The output data are the spatially dependent velocity, pressure, concentration, and potential fields. Other characteristics including analyte mass, dispersion, and flow rate can be calculated from the output data. RESULTS AND DISCUSSION Electrokinetic Focusing Using a Channel Cross Structure. In the first series of simulations, electrokinetic focusing experiments described previously in ref 63 were modeled. The actual (60) Saad, J.; Schultz, M. H. Math. Comput. 1985, 44, 417. (61) Kershaw, D. S. J. Comput. Phys. 1978, 26, 43. (62) Ghia, U.; Ghia, K. N.; Shin, C. T. J. Comput. Phys. 1982, 48, 387. (63) Jacobson, S. C.; Ramsey, J. M. Anal. Chem., 1997, 69, 3212.
Table 1. Field Strength (kV/cm) in the Chip Channels channel
run a
run b
run c
waste sample focusing (left and right)
2.82 1.12 0.85
2.63 0.61 1.01
2.43 0.11 1.16
sample focusing experiments were performed in a channel cross similar to the one drawn schematically in Figure 1. The simulation domain is shown in the lower part of the figure. For the simulations, the channel width was assumed to be 24 µm. This corresponded to the wider trapezoid base of the actual experimental channel. Three experiments with different focusing fields were simulated. Field strength values were taken from Figure 3b-d in ref 63 (see Table 1) and were used as boundary conditions for the entry ports of the channels in the simulations. Cationic sample (Rhodamine 6G) with an electrophoretic mobility of µep ) 1.4 × 10-4 cm2/(V‚s) and a diffusion coefficient D ) 3.0 × 10-6 cm2/s was used as the material to be spatially focused. It was assumed that in all experiments the electroosmotic mobility had the same value, µeo ) 4.0 × 10-4 cm2/(V‚s). Simulated concentration distributions are compared in Figure 2 with experimentally obtained CCD images (reproduced from ref 63). This figure shows three pairs of pictures corresponding to different focusing fields. In each pair, the left picture shows an experimental image, and the right one shows the simulated image. The channel edges have been highlighted in one of the experimentally obtained images, Figure 2a (left). The laserinduced fluorescence (LIF) is given in arbitrary units. The gray scale in the simulated images represents the concentration distribution normalized to the initial sample concentration. In general, the simulated and experimental images demonstrate good correlation, although a more detailed study reveals some differences. The experimental images have more diffused concentration boundaries, and the sample stream appears wider than in the simulations. A more quantitative comparison is presented in Figure 3, in which concentration profiles taken at four samplewaste channel cross sections are plotted. These cross sections Analytical Chemistry, Vol. 70, No. 21, November 1, 1998
4499
Figure 3. Simulated (solid line) and experimental (dashed line) concentration profiles across the sample-waste channel: (1) 3 µm upstream from the focusing chamber, (2) 9 µm into the focusing chamber, (3) 18 µm into the focusing chamber, and (4) 60 µm downstream from the focusing chamber
are shown in Figure 2c (left panel) by dashed lines. For comparison purposes, all concentration and LIF values are normalized. The simulated profiles (solid lines) in panels 2 and 3 are narrower than the experimental ones (dashed lines). In panel 1, the experimental profile gradually approaches zero at the channel edges, whereas in the simulation, the concentration on the wall has a value of ∼0.4. Among possible explanations for these discrepancies between experiment and theory are the simplifications adopted in the model and incomplete experimental data. In the first case, we refer to the fact that a 2-D model instead of a 3-D one was considered. The model, therefore, does not account for the actual trapezoidal shape of the channel cross section. The LIF signal observed in experiments is integrated over the depth of the channel. Since the channel depth close to the edges is variable (trapezoid shape), the “optical depth”, i.e., the thickness of layer emitting the light, is also variable. Therefore, close to the channel edges the intensity of light gradually diminishes, giving the impression that the sample stream is smeared. This can probably explain the observed results in Figure 3 panel 1. In addition, as the channel width changes with depth, the focusing field may also undergo slight variations along the z-coordinate. The electroosmotic and electrophoretic mobilities were not measured explicitly in these experiments, and, therefore, some typical values were taken for the simulation. The above experiments pose a question about how tightly the sample stream can be focused. For a fixed field strength in the waste channel, electrokinetic focusing has two limiting cases: in the first case, the focusing field, Ef, is absent (Ef ) 0) and the field in the sample channel, Es , is equal to that in the waste channel (Ew ) Es); in the second case, Es ) 0, while Ef ) 0.5Ew. In the first case the sample flow into the waste channel is close to maximum (small losses occur due to diffusion into focusing channels), while in the second it should be zero. However, if the sample channel has been previously filled with a sample, a small portion of it could still penetrate into the focusing chamber 4500
Analytical Chemistry, Vol. 70, No. 21, November 1, 1998
(channel intersection) by diffusion, where it will be captured by the field coming from the focusing channels and moved into the waste channel. Simulations were performed to confirm this effect. The channel geometry was the same as described above except that the channel width is assumed to be 18 µm. This value corresponds to the average channel width, not the largest one as before. The field strength in the waste channel was fixed at Ew ) 2 kV/cm. The values in the sample and focusing channel were varied under this constraint. Rhodamine B (neutral) was used as a sample (diffusion coefficient 3 × 10-6 cm2/s); the electroosmotic mobility was assumed to be 3 × 10-4 cm2/(V‚s). Simulated concentration profiles along the sample stream axis are plotted in Figure 4a. The profiles normal to the sample stream in the cross section at the exit of the focusing chamber are shown in Figure 4b. The different lines correspond to different field strengths in the sample/focusing channels. As expected, with the increase of Ef the sample stream becomes more tightly focused. Concurrently, the sample peak concentration in the focusing chamber drops quickly as Es approaches zero. At Es ) 0, the sample stream is still diffusively entering the focusing chamber (line 1, Figure 4), although its maximum concentration at the exit of the focusing chamber is ∼30 times lower than its initial value. At the same time, the stream width is only about 2 times smaller compared to the case in which Ef ) Es ) 1/3Ew (line 5, Figure 4b). As the focusing field is increased, the losses in sample concentration (Cmax) are much greater than the gain in focusing. Apparently, an optimum focusing scheme exists in which the ratio Cmax/χ (χ is the standard deviation calculated for the cross section concentration profile) reaches a maximum. In this particular situation, it is close to the case where Ef ) 0.4Ew, Es ) 0.2Ew (see table insert on Figure 4). From the simulations discussed above, two more conclusions can be drawn. First, electrokinetic focusing does not increase the sample concentration in the stream, but it rather confines the stream spatially. Second, by adjusting the field strength in the sample and focusing channels, it is possible to tune the amount of sample in the stream which is eventually going to be in the injected plug. The latter has been experimentally verified in an MEKC system16 and in a gel-filled microchannel64 with the electroosmosis suppressed. von Heeren et al. empirically optimized the focusing field so that the amount of sample injected is small enough to produce a symmetric sample peak, but on the other hand, the volume injected is sufficiently large not to hamper detection.16 In the next series of simulations, the impact of such factors as electrophoretic and electroosmotic mobilities on sample focusing was studied. Hypothetical cations, anions, and neutral molecules were used as samples. For charged ones, an absolute value for electrophoretic mobility was equal to 1.5 × 10-4 cm2/(V‚s), and the electroosmotic mobility was 3.0 × 10-4 cm2/(V‚s). The diffusion coefficient for each species was assumed to be 3 × 10-6 cm2/s. Results from this study are given in Figure 5, which shows the sample concentration profiles for the cross sections at the exit of the focusing chamber. Analysis of these data reveals that the cation sample stream (thin solid line with diamonds) is better focused than the neutral molecule stream (thick solid line), which (64) von Heeren, F.; Verpoorte, E.; Manz, A.; Thormann, W. J. Microcolumn Sep. 1996, 8, 373.
Figure 4. Concentration profiles along the sample stream axis (a) and normal to it at the exit of the focusing chamber (b). Electric field strength in the waste channel was held fixed at Ew ) 2 kV/cm. The sample and focusing channel electric field strengths are as follow: (1) Es ) 0, Ef ) 0.5Ew; (2) Es ) 0.04Ew, Ef ) 0.48Ew; (3) Es ) 0.1Ew, Ef ) 0.45Ew; (4) Es ) 0.2Ew, Ef ) 0.4Ew; (5) Es ) 1/3Ew, Ef ) 1/3Ew; (6) Es ) Ew, Ef ) 0.
Figure 5. Concentration profiles across the waste channel at the exit of the focusing chamber. The various profiles correspond to different samples and operating conditions: solid line with markers, cation µep ) 0.5µeo ) 1.5 × 10-4 cm2/(V‚s); free solid lines, neutral, thin µeo ) 1.5 × 10-4 cm2/(V‚s), bold µeo ) 3 × 10-4 cm2/(V‚s); other lines, anion, dashed µep ) 0.5µeo ) 1.5 × 10-4 cm2/(V‚s), dotted µep ) 1.5 × 10-4 cm2/(V‚s), µeo ) 0, dash-dot µep ) 0.1µeo ) 3 × 10-5 cm2/(V‚s).
in turn is better focused than the anion stream (dashed line). This is confirmed by experimental data.63 The maximum concentration values demonstrate the same dependence: highest for the cation and lowest for the anion. The cation stream is focused by the combined effect of electrophoresis and electroosmotic flow, while for a neutral molecule only the electroosmotic flow confines the
sample stream. In the case of an anion, the electrophoretic and electroosmotic effects act in opposite directions. However, if the electrophoretic mobility is small with respect to electroosmosis, the effect of electrophoretic focusing (or defocusing) is also small (see dash-dot line for the anion with µep ) 0.1µeo ) 0.3 × 10-4 cm2/(V‚s)). The dotted line represents the case when electroosmotic flow is absent and electrophoresis (µep ) 1.5 × 10-4 cm2/ (V‚s)) is the only focusing factor. Comparing this with purely electroosmotic flow (dashed line) and with the case of combined electrophoretic and electroosmotic effect (line with diamonds) shows the impact of electroosmosis and its interaction with electrophoresis. It is worth noting that, although the absolute value of the anion’s total electrokinetic mobility is the same in cases when the electroosmosis is present and when it is absent, (µep - µeo) ) µep ) 1/2 µeo, the corresponding concentration profiles are different. The combined effect of electroosmosis and electrophoresis (dashed line) produces a more focused profile with lower concentrations at the maximum compared to just electrophoretic focusing (dotted line). The same is true when electroosmotic focusing alone is compared to electrophoretic focusing and the absolute values of mobilities are equal, µep ) µeo ) 1.5 × 10-4 cm2/(V‚s) (thin solid line without markers versus dotted line). Sample Mixing in Tee-Channel Structure. The second simple channel geometry studied in simulations is the tee intersection. Such channel geometries are an important element in chip design since they can be used as mixers.5,21,23 Figure 6 shows an example of simulations in which the sample (Rhodamine B) is mixed with pure buffer in equal proportions. The channel width is 20 µm, the field strength in the downstream (mixing) channel is 100 V/cm, and in the side channels it is 50 V/cm. Because of the plug-shaped flow profile, i.e., uniform velocity over Analytical Chemistry, Vol. 70, No. 21, November 1, 1998
4501
Table 2. Simulated Sample Dilution for Several Field Strength Schemesa field coefficients Es/Eb
average minimum maximum standard deviation expected value a
0.9/0.1 N 0.9172 0.8595 0.9695 0.0415 0.9
0.6/0.4 N 0.59995 0.3992 0.7952 0.14943 0.6
0.5/0.5 N 0.49995 0.2921 0.7073 0.15668 0.5
0.4/0.6 N 0.4018 0.2229 0.6000 0.1467 0.4
0.2/0.8 N 0.19995 0.0802 0.3280 0.0935 0.2
0.5/0.5 C 0.5001 0.1983 0.8009 0.2277 0.5
0.4/0.6 C 0.4003 0.1232 0.6973 0.2170 0.4
0.2/0.8 C 0.2001 0.0380 0.3960 0.1355 0.2
0.5/0.5 A 0.4993 0.4325 0.5661 0.0504 0.5
0.2/0.8 A 0.1999 0.1590 0.2410 0.0309 0.2
0.1/0.9 A 0.1056 0.0824 0.1288 0.0175 0.1
N, neutral; C, cation; A, anion.
Figure 6. Sample (shown in white) mixing with buffer (black) in a tee channel structure. Vector lines show the sample particle motion, and their lengths reflect the relative velocities.
the channel cross section, the only source of sample mixing is diffusion. To obtain a uniform sample distribution across the channel, the mixing time should be long enough for the sample to diffuse across the channel. This can be achieved by using either a long mixing channel or a low field strength. The effect of the field strength on the sample concentration profile across the mixing channel is shown in Figure 7 (longer lines). The concentration profiles correspond to the cross section 70 µm downstream from the mixing tee. By decreasing the field strength in the mixing channel from 500 to 100 V/cm, the maximum concentration inhomogeneity drops from 90% of its initial value to about 30%. An alternative, which offers better and faster mixing, is a chip design in which the mixing channel is narrower than the side channels. Concentration profiles for such a chip with the mixing channel width one-half that of the side channels (10 µm) show much better homogeneity (Figure 7 shorter lines), especially when the field strength in the mixing channel is 200 or 100 V/cm. The improved mixing for narrower channels can be explained by analyzing the characteristic time scale for diffusion tdiff ) L2/D, which drops as the second power with the decrease of the length scale L (channel width). 4502 Analytical Chemistry, Vol. 70, No. 21, November 1, 1998
Figure 7. Sample concentration profiles across the mixing channel for different field strengths and channel widths. Initial sample concentration is 1. Data were collected at a cross section 70 µm downstream from the mixing tee. The following field strengths were used in the mixing channel: dash-dot line, 500; solid, 200; and dash, 100 V/cm. Longer lines are for 20-µm-wide mixing channel, and shorter lines correspond to a 10-µm-wide mixing channel.
An important feature of the tee mixer is the sample dilution characteristic, i.e., the dependence of the sample concentration in the mixing channel on the field distribution in the side channels. The possibility of manipulating the sample dilution by tuning the field was demonstrated experimentally by Kutter et al.21 They showed that, for neutral species, the temporal response of the dilution profile follows the temporal field profile in the side channels for different shapes of that profile. The question that remained was whether the dilution profile followed the quantitative characteristics of the field for charged species. For instance, if the ratio between the field magnitude in the sample and the buffer channel is 0.4, will the sample concentration in the mixing channel be 40% of its initial value? Such characteristics were studied by means of simulations for both neutral and charged species. In the simulations, the side channels were 20 µm wide, and the mixing channel was 10 µm wide. In the mixing channel, the field strength was fixed at E0 ) 500 V/cm, whereas the side channels were varied from 0.9E0/2 to 0.1E0/2. The factor of 1/2 appears since the side channels were twice as wide. The concentration was monitored in the mixing channel cross section 70 µm downstream from the intersection. The hypothetical substances and the electroosmosis characteristics used in these simulations were the same as described before. Table 2 gives the average
concentration value over the mixing channel cross section together with its minimum and maximum values and the standard deviation for several field schemes in the side channels. The field schemes are characterized by coefficients, which are multiplied by E0/2, giving the field strength in the sample, Es, and buffer, Eb, channels. Thus, 0.4/0.6 means Es ) 0.4E0/2 ) 100 V/cm and Eb ) 0.6E0/2 ) 150 V/cm. For convenience, the presumed sample concentration value is also given for each variant. As is evident from the data, the average sample concentration in the mixing channel is very close to the expected value for all species. This holds true despite a significant concentration inhomogeneity over the channel cross section. The major difference between different species is the magnitude of that inhomogeneity, which is largest for cations and smallest for anions. Obviously, cations spend the least amount of time in the mixing channel because the electroosmotic and electrophoretic velocities are summed. On the other hand, anions move slowest and, consequently, have more time to diffuse. As a result, a more uniform concentration profile is observed. Small deviations in the average sample concentration (∼5% from the expected value) appear only close to the extremes when Es/Eb ) 9 or Es/Eb ) 0.11.
k
Botzmann’s constant
K
association/dissociation equilibrium constant
L
length scale
n
concentration of charged species
p
pressure
Pe
Pe´clet number
rD
Debye radius
S
channel cross section area
t
time
T
absolute temperature
V
velocity
w
channel width
x, y, z
spatial coordinates
zi
electric charge (in electrons)
CONCLUSIONS The mathematical model and its computer implementation demonstrate that computer simulations can be a useful tool in the analysis of electrokinetic transport in microfabricated devices. Two basic electrokinetically driven microfluidic elements, the channel cross and the tee intersection, were targets for simulation. Sample transport phenomena specific to these elements were studied: electrokinetic focusing in the channel cross and sample mixing/dilution in the tee intersection. Good agreement between simulated and experimental data for electrokinetic focusing confirms that the mathematical model is accurate and can be used for the analysis and prediction of experimental results. It was found that the tightest sample focusing conditions may not be the most suitable for detection purposes, as dilution effects may outweigh the benefits of a narrower sample stream. Specific field combinations in the focusing and sample channels are shown to provide optimum focusing and detection conditions. Simulations performed for the tee intersection show that sample mixing can be enhanced by using a narrower mixing channel and a lower field strength. It was also demonstrated that the tee mixer provides very accurate dilution within a wide range of dilution ratios, regardless of the charge of the sample species. LIST OF SYMBOLS
Greek R
anion ionization degree
β
cation ionization degree
∇
nabla operator (divergence or gradient depending on a variable)
∇2
Laplacian operator
dielectric permittivity
0
dielectric permittivity of vacuum, 8.854 × 10-12 C V-1 m-1
φ, Φ
electric potential (internal, external field)
χ
standard deviation calculated for the cross section concentration profile
η
viscosity
λ
thermal conductivity
µ
electroosmotic/electrophoretic mobility
ν
kinematic viscosity
Fe
charge density
σ
specific conductivity
τ
time increment (in simulations)
ζ
wall zeta potential
Subscripts n
number of anionic species
m
number of cationic species
ep
electrophoretic
eo
electroosmotic
D
Debye
a
anion concentration
f
focusing
b
cation concentration
i
species number
c
capillary wall thickness
s
sample
D
diffusion coefficient charge of electron, 1.6 × 10-19 C
w
waste
e
b
buffer
E
electric field intensity
F
Faraday’s constant, 96 485 C/mol
h
depth of the channel
a
anion
H
hydrogen ion concentration
b
cation
j
current density
n
time step number
Superscripts
Analytical Chemistry, Vol. 70, No. 21, November 1, 1998
4503
ACKNOWLEDGMENT This research is sponsored by Oak Ridge National Laboratory Seed Money Program. It was supported in part by an appointment for S.V.E. to the Oak Ridge National Laboratory Postdoctoral Research Associates Program, administered jointly by the Oak Ridge National Laboratory and the Oak Ridge Institute for Science and Education. Oak Ridge National Laboratory is managed by
4504 Analytical Chemistry, Vol. 70, No. 21, November 1, 1998
Lockheed Martin Energy Research Corp. for the U.S. Department of Energy, under Contract DE-AC05-96OR22464.
Received for review May 20, 1998. Accepted August 20, 1998. AC980551W
