Simulation of the Generation of Solution Gradients in Microfluidic

Nov 16, 2011 - Zhongbin Xu , Xing Huang , Pengfei Wang , Huanan Wang , David A. Weitz ... Konstantinos Zografos , Robert W. Barber , David R. Emerson ...
0 downloads 0 Views 2MB Size
ARTICLE pubs.acs.org/IECR

Simulation of the Generation of Solution Gradients in Microfluidic Systems Using the Lattice Boltzmann Method Yangxu Hu,† Xianren Zhang,*,† and Wenchuan Wang† †

Division of Molecular and Materials Simulation, State Key Laboratory of OrganicInorganic Composites, Beijing University of Chemical Technology, Beijing 100029, China ABSTRACT: In this work, we used the lattice Boltzmann (LB) method to simulate the generation of concentration gradients in microfluidic networks. We first developed a model of microfluidic networks, in which the flow velocity of the laminar flow in the microchannel and the diffusion of solute molecules govern the solute gradients generated, and the obtained results are comparable with experimental data. Our simulation results indicated that the relative positions of the branching plates in different levels, at which a microchannel is split into several daughter microchannels, exert significant effects on the shape of concentration gradients. We also performed extensive simulations to study the dependence of the shape of the concentration gradient on the velocity of the flow u, diffusion coefficients of solutes D, and the length of the microchannel L. A dimensionless parameter, (L/u)/(H2/D), in which H is the hydraulic diameter of the main microchannel, was proposed in this work. It is found that for geometrically similar microfluidic networks, the parameter alone determines the shape of the generated concentration gradient. Therefore, the proposed parameter allows one to perform experiments or simulations with reduced models in microchannels and correlate the data to the other flows or geometry sizes.

1. INTRODUCTION Gradients of diffusible chemicals play an important role in some biological and chemical processes. For a variety of chemical processes, such as reactions,1 nucleation, and the growth of crystals,2 solution gradients are needed to be generated for different purposes. Well-defined gradients on surfaces and in solutions have also been extensively used to study chemotaxis,3 cell-signaling,4 and other biological processes. Technology of the generation of gradients on surfaces or in solutions in microfluidic systems with different fabrications has recently achieved rapid progresses. A universal gradientgeneration device has been proposed recently to generate gradients through controlling the diffusive mixing of species on surfaces or in solutions.5,6 The gradient-making network and its variations are suitable to a variety of fluids and solutions, and it can be used to produce gradients of different shapes, types, and kinds.5 Irimia et al.7 proposed a systematic approach to generate stable gradients of any profiles by the controlled mixing of two starting solutions. This model shows its flexibility in producing nonlinear gradient shapes. The microfluidic device designed by Li et al.8 can produce linear and nonlinear gradients along a microfluidic channel by microtunnel controlled stepwise addition of a sample solution. To obtain fixed concentration gradients, Holden et al.9,10 designed and built a laminar microfluidic diffusion diluter, in which the main channel is divided into a parallel array of microchannels. In addition, several different approaches were also proposed to obtain different gradient shapes.1113 There are two main advantages of the microfluidic-based method for the generation of concentration gradients: (i) Flexible gradients can be obtained for various gradient shapes by using different channel network designs. (ii) Stable gradients can be maintained for a long period of time. However, a limitation of this method is that the gradient shape is fixed for each device, r 2011 American Chemical Society

and each device can only be used in a static experiment. To solve this problem, a single microfluidic device was designed to generate dynamic temporal and spatial concentration gradients, 14 and a modular approach was proposed to construct microfluidic systems for the generation of gradients of arbitrary profiles.15 Recent development of integrated microfluidic devices in different areas can refer to the review work by Erickson.16 In this work, we develop a model of microfluidic networks and use the lattice Boltzmann method to simulate the generation of concentration gradients in microfluidic networks. Branching plates are placed at relative positions in different levels to get concentration gradients. Besides, we also investigate the effects of several factors, including fluid velocity, solute diffusion coefficient, the size of microchannel, and the length of the microchannel, on the generation of concentration gradients. More importantly, a parameter is proposed here to describe the concentration gradient at the outlet.

2. MODEL AND METHOD 2.1. Model. Joen et al.5 and Dertinger et al.6 proposed a

microfluidic network to generate gradients on surface or in solutions through repeatedly mixing and splitting of two solutions. Emerson et al.17 and Barber et al.18 also designed a microfluidic network to study biomimetics. The two networks adopted similar mechanisms to generate concentration gradients. In this work, we propose a model to simplify these Received: May 26, 2011 Accepted: November 16, 2011 Revised: November 14, 2011 Published: November 16, 2011 13932

dx.doi.org/10.1021/ie201136r | Ind. Eng. Chem. Res. 2011, 50, 13932–13939

Industrial & Engineering Chemistry Research

ARTICLE

Figure 1. Model of microfluidic networks. (a) The structure of microfluidic networks. (b) The simplified model of microfuidic networks. l0, l1, l2, and l3 are the length of the branch 0, branch 1, branch 2, and branch 3, respectively. C1C8 are the concentrations at the outlet of the pathway, respectively. (c) A schematic diagram of baffles and plates in the mcirochannel.

Figure 2. (a) Our simulation model. (b) Gradient profile at the outlet of the microchannel. The black square (0) is our simulation result, while the red line is the experimental result in ref 6. The simulation is carried out at the condition of L = 600, u = 0.0278, d = 0.01 (the box size is 600  60  60).

two networks for computer simulations, which makes the lattice Boltzmann method applicable, as is shown in Figure 1. In our

model, two streams enter into the fluid inlet of a square microchannel (Figure 1c) with the same velocity separately. Fluids 1 and 2 contain components B and A, respectively, and both of them are combined with the same solvent (e.g., A is salt, B is sucrose, and the solvent is water). Through the interface between the two fluids, components A and B can diffuse into fluids 2 and 1, respectively. Since the fluids are in laminar flow moving slowly, the mixing process that component A transfers from fluid 2 to fluid 1 and B from fluid 1 to fluid 2 is determined by diffusion rather than convection. After a certain distance of the diffusion process, the two fluids are separated by a plane to prevent further diffusion of component A from fluid 2 to fluid 1 and that of component B from fluid 1 to fluid 2. A baffle exists in the front of the plane to make the model more reasonable (Figure 1b and c). Thus, fluid 1 is of a high concentration of B and a low concentration of A. Since the diffusion is rather slow, component A in fluid 1 is not uniform, which leads to the generation of the concentration gradient of A, while the concentration gradient of B can be generated. Through repeatedly diffusing and separating of the two fluids, concentration gradients of A and B would be built up at the outlet. In general, the simplified model proposed here (Figure 1) is not only easy to use in computer simulations but also flexible enough to model different microfluidic networks. It can generate a series of gradients with different shapes and kinds. More importantly, our model catches the essential factors in processes of generation of concentration gradients: repeating the processes of mixing flows with different solvent concentrations and then splitting. Furthermore, through controlling the velocity of fluid in different inlets, our model can also be used to generate dynamic concentration gradients. 13933

dx.doi.org/10.1021/ie201136r |Ind. Eng. Chem. Res. 2011, 50, 13932–13939

Industrial & Engineering Chemistry Research

ARTICLE

Figure 3. Velocity field of the xy section at the middle width of the microchannel. The velocity field shows that the flow in the microchannel is laminar flow, and the mass transfer of the solute is controlled cooperatively by the diffusion in the solution and the flow.

2.2. Method. Here, we used the lattice Boltzmann (LB) method with a 3-dimensional D3Q19 (19 velocities on a simple cubic lattice of dimension 3) lattice Bhatnagar-Cross-Krook (BGK) model19,20 to simulate the generation process of a concentration gradient. In comparison with traditional computational fluid dynamics (CFD) methods starting from the continuum equations, the LB method has advantages in dealing with complex boundaries and incorporating of microscopic interactions, which are crucial in many microfluidic circumstances. The origin of the LB method can be traced back to the lattice-gas cellular automata (LGCA).21 In the LB method, the evolution of the particle distribution function fi satisfies the following equation19 eq

fi ðx þ ei , t þ 1Þ ¼ fi ðx, tÞ  ðfi ðx, tÞ  fi ðx, tÞÞ=τ

ð1Þ

where x is a node position in the lattice, ei is the i velocity allowed, e.g. i = 0, 1, 2,...,18 in the D3Q19 model, t is the time step, feq i is the equilibrium distribution, and τ is the BGK relaxation parameter, which determines the viscosity and diffusion coefficients. So, in our simulation the time evolution of Fi is governed by the

single relaxation time lattice Boltzmann equation eq

Fi ðx þ ei , t þ 1Þ ¼ Fi ðx, tÞ  ðFi ðx, tÞ  Fi ðx, tÞÞ=τ

ð2Þ

where Fi is the local density along the ei direction at a certain node position, x. The equilibrium distribution feq i in the D3Q19 model is given 22 by " # ci 3 u ðci 3 uÞ2 u2 eq fi ¼ ωi F 1 þ 2 þ  2 ð3Þ cs 2c4s 2cs Furthermore, the equilibrium distribution must satisfy the following constraint22

∑i fi eq ¼ F, ∑i ci fi eq ¼ Fu

ð4Þ

where ωi is a set of weight parameters, F is the total local density, ci is the discrete set of velocity, u is the fluid bulk velocity, and cs = (RT)1/2 is the sound velocity. 13934

dx.doi.org/10.1021/ie201136r |Ind. Eng. Chem. Res. 2011, 50, 13932–13939

Industrial & Engineering Chemistry Research

ARTICLE

Figure 4. Gradients profile established at the outlet of the microchannel. (a) shows the concentration gradients profile generated at the outlet of the microchannel. (b) The red square (0) indicates the concentration variation of component A along the y-direction at the plane of z = 30.

Under the condition of the low Mach number (Ma , 1 with Ma = u/cs), the Chapman-Enskog method is used to obtain the macroscopic diffusion eqution22 ∂F þ ∇ 3 ðFuÞ ¼ 0 ∂t

ð5Þ

∂ðFuÞ þ ∇ 3 ðFuuÞ ¼  ∇p ∂t þ ∇ 3 ½Fνð∇u þ ð∇uÞT Þ

ð6Þ

where ν is the viscosity, a parameter related to relaxation time, v = c2s (τ  1/2)δt. The Reynolds number of the flow is given by Re = (Fude)/(μ), where μ is the viscosity of the fluid, and de is the characteristic length, namely the hydraulic diameter of the microchannel. The hydraulic diameter can be expressed as de = (4Ae)/(P), where Ae and P are the area and the perimeter of the microchannel cross-section, respectively.

3. RESULTS AND DISCUSSION 3.1. Validation of the Simulation Method. In order to validate our algorithm used in this work, we designed a model, corresponding to the experimental network proposed by Dertinger et al.,6 as is shown in Figure 2a. In our model, two streams enter into the fluid inlet of a square microchannel with the same velocity separately. Fluids 1 and 2 contain components A and B, respectively, and both of them are combined with the

same solvent (e.g., A is salt, B is sucrose, and the solvent is water). First, two fluids are separated at branch 0 by a plane to prevent the diffusion of component A from fluid 1 to fluid 2 and component B from fluid 2 to fluid 1. Then, the two fluids flow into branch 1, a part of fluids 1 and 2 flows into one channel, and diffusion takes place thereafter. Branch 1 here corresponds to the inlet of the experimental model.6 The fluids flow through branches 2, 3, and 4, then. Finally, the concentration gradient profile (see Figure 2b) is generated at the outlet. We can find our simulation data are in good agreement with the experimental results, which proves that the algorithm used in this work is reasonable and reliable. 3.2. Generation of Gradients. In this work, a series of simulations were performed to investigate the behavior of the fluid in the microfluidic networks. To generate gradient profiles, the microchannel was characterized by a three-dimensional box with Lx = 480, Ly = 64, and Lz = 64, which correspond to the length (L), the height (H), and the width (W) of the microchannel, respectively. The size of the baffle plate (lx = 1, ly = 4, and lz = 64, see Figure 1c) along the y-direction (b = ly) is set to 4. The solvent density is set to 3.0, while both the initial local densities for components A and B are 0.15. Both fluids flow at a velocity of 0.03, with the solvent viscosity of 0.98. The diffusion coefficients for components A and B are set to 0.01 (all the above parameters are reduced with respect to the lattice constants). Note that the Reynolds number of the flow in the microchannel is about 5.877, which is low enough to induce a laminar flow. In the LB method, the concentration of A or B can be given by the local density of A or B. To investigate the process of the gradients generation, the velocity field of the xy plane with z = 30 (the middle of the channel) is given in Figure 3. The figure shows that the fluid in the microchannel is laminar, and the mass transfer of the solutes along the y-direction is controlled by diffusion, since there is mainly no velocity along the y-direction (Figure 3). The fluid velocity near the channel surfaces is smaller than that in the middle of the microchannel, which indicates the formation of stagnant layers near the surfaces (see Figure 3 for details). At the outlet (x = 479) of the microchannel, the concentration gradients are established, and a typical concentration profile is shown in Figure 4. Figure 4a shows the profile of the concentration of component A, and it demonstrates that the flow at the outlet presents an apparent concentration gradient. Since component B shows a similar behavior as component A, only the results of component A are discussed here. We also define a dimensionless parameter, c = C/C0, where C0 is the initial concentration of A, and C is the local concentration of A. Near the boundary (i.e., z = 0 and z = 60), the local velocity is smaller than that of the bulk solution because of the effect of the surfaces. The concentration distribution of the bulk solution in the channel shows approximately a linear change as a function of the distance from the center of the channel. Thus, along the z-direction the concentrations next to the surfaces are somewhat different from those in the middle of the channel. Figure 4b shows the concentration profile along the y-direction at z = 30. Furthermore, we performed extensive simulation runs to study the dependence of the shape of concentration gradients on different parameters, such as the velocity of the flow, properties of fluids, and the size of the microchannel. 3.3. Effects of the Microchannel Structure and Fluid Property on Gradient Generation. First, the effects of the microchannel’s structure, including the position of the plates and 13935

dx.doi.org/10.1021/ie201136r |Ind. Eng. Chem. Res. 2011, 50, 13932–13939

Industrial & Engineering Chemistry Research

ARTICLE

Figure 5. (a) Effects of baffle width on the generation of the gradients. The concentration gradients at the outlet of the microchannel are roughly the same when the widths of baffle correspond to b = 0, b = 2, b = 4, and b = 6 (all simulations were performed at the condition of u = 0.03, L = 480, D = 0.004). (b-d) Effect of the position of different plates on the generation of gradient. Here, some parameters are fixed, such as velocity (u = 0.030), the length of the microchannel (L = 480), and the diffusion of solute (D = 0.004). (b) shows the effect of plate I: l0 = 40, l0 = 120, (For clarity, the concentrations are shifted to the right with a distance of 10.), l0 = 200 (shifted to the right with a distance of 20), and l0 = 240 (shift to the right with a distance of 30). (c) shows the effect of the positions of plates of plate II: l0+l1 = 160, l0+l1 = 240 (shifted to the right with a distance of 10), l0+l1 = 320 (shifted to the right with a distance of 20), l0+l1 = 360 (shifted to the right with a distance of 30). (d) shows the effect of the positions of plates of plate III: l0+l1+l2 = 264, 300, 360, 420. (Note that the some concentrations are also shifted for clarity.)

the size of the baffle, on the gradient generation were investigated. The simulations with different baffle sizes (b = 0, 2, 4, 6, see Figure 1) were performed, and the results corresponding to the different baffle sizes are shown in Figure 5a. When the fluid flows slowly, the baffle shows negligible effects on the behavior of the fluid, especially on the diffusion of solutes. In general, all the chosen velocities (from u = 0.010 to u = 0.040) in this work are in the range, where the effect of the baffle plate can be neglected. Hence, the width of the baffle plate is set to 0 below. Plates at different positions in this model correspond to the different length of branched channels in the microfluidic networks. The position of a plate is characterized by the distance between the plate and the inlet of the microchannel. For example, the position of plate I is denoted by l0, plate II is l0+l1, and plate III is l0+l1+l2, as shown in Figure 1b (l0 = l1 = l2 = l3 = 0.25 L, where L is total length of the microchannel). By changing the position of plates of the same type along the x-direction while fixing those of the other types, the concentration gradients at the outlet of the microchannel would be changed accordingly, as are shown in Figure 5b-d. Figure 5b shows the dependence of concentration gradient on the position of plate I, which varies from x = 0 (l0 = 0) to x = 240 (l0 = 240) in the x-direction, while the positions of the other plates remain unchanged. The results show that the increase of l0 leads to a more smooth shape of concentration gradients. This is because the increase of l0 causes the diffusion of solute more sufficiently; as a result, the difference of solute concentrations between two branched channels becomes smaller. Figure 5c shows the effect of the position of plate II, and the figure indicates that the position of plate II determines the

concentration difference between neighboring branch II, δC2, which is equal to C2C3 or C6C7. Since δC2 increases with the decrease of l1, in order to get a smoother concentration gradient, l1 should decrease correspondingly, which means that the positions of plate II should be moved away from the inlet of the microchannel. We also studied the effects of the position of plate III (Figure 5d). The figure shows that the decrease of l2 tends to generate a larger concentration difference between neighboring branch III, δC3, which is equal to C1C2, C3C4, C5C6, or C7C8. So, in order to get a relatively smooth gradient profile, we fix a box with l0 = 1/3L, l1 = l2 = 1/6L, l3 = 1/3L to carry out simulations for them. Our simulations also indicate that the morphology of networks significantly affects the shape of the concentration gradients generated. For example, when we changed the positions of daughter plates along the y-direction from the microchannel (Figure 1b), a new microfluidic network (Figure 6a) formed, and the obtained profile of concentration gradients (Figure 6b) was also significantly changed. The velocity of the flow can also affect the generation of concentration gradients. Figure 7a shows the concentration gradients at several fluid velocities in the microchannel. The flow velocity affects the gradients profile through changing the diffusion time of solute molecules in the microchannel. The high fluid velocity leads to the decrease of diffusion time, which induces a large concentration difference between two neighboring channels. In contrast, if the fluid flows slowly, the concentration difference of the solutes along the y-direction would decrease. A similar trend was observed previously by experimental and theoretical studies. 9,10 Hence, the flow velocity plays an 13936

dx.doi.org/10.1021/ie201136r |Ind. Eng. Chem. Res. 2011, 50, 13932–13939

Industrial & Engineering Chemistry Research

ARTICLE

Figure 6. (a) Another simplified model of microfluidic networks, in which l0, H, and the corresponding position for daughter branches along the x-direction are the same as in Figure 1b. (b) Gradient profile of component A at the outlet of the microchannel. The simulation is carried out at the condition of L = 480, u = 0.025, and d = D = 0.01.

important role in the generation of the solute gradient, and its effect depends on the relative length of the splitting plate. The diffusion coefficient of solute molecules is found to exert a similar effect on the generated gradient profile (Figure 7b). This is because, to generate a solute gradient, a high diffusion coefficient works the same as the slow fluid flow. Finally, the diffusion time for solute increase with the length of the microchannel. Hence, the increase of the length of the microchannel causes a similar variation of the concentration profile (Figure 7c) as that is induced by the decrease of the flow velocity. Therefore, for different networks the generation of concentration gradients follows the same mechanism, namely, it is the competition between the flow velocity of the laminar flow in microchannels and the diffusion of solute molecules that govern the solute gradients generated. The competition will be discussed in the following section in detail. 3.4. Dimensionless Parameter for Characterization of the Concentration Gradient. The time required for fluid flowing through the microchannel can be roughly estimated as L/u, while a typical diffusion time for solutes reaching the surfaces can be written as t∼H2/D, where H is the characteristic width of the microchannel and D is the fluid diffusion coefficient. As a result, the generated concentration gradient at the outlet of the microchannel depends on the completion between the solute diffusion and fluid flow. A dimensionless parameter is defined here as (L/u)/(H2/D) to characterize the completion between the solute diffusion and

Figure 7. (a) Effects of the velocity on the generated concentration gradient. The concentration gradients at large flow velocities are shifted to the right for clarity with different distances: 10 for u = 0.015, 20 for u = 0.020, 30 for u = 0.025, 40 for u = 0.030, 50 for u = 0.035, and 60 for u = 0.040 (All simulations were performed at the condition of L = 480, D = 0.004.). (b) Effects of the diffusion coefficients on the generated concentration gradient. Four concentration profiles corresponding to large diffusion coefficients are shifted to the right for clarity with different distances: 10 for D = 0.007, 20 for D = 0.010, 30 for D = 0.013, and 40 for D = 0.016 (All simulations were performed at the condition of u = 0.03, L = 480.). (c) Effects of the length of the microchannel on the generated concentration gradient. The concentration gradients at large microchannel lengths are shifted to the right for clarity with different distances: 10 for L = 420, 20 for L = 480, 30 for L = 540, 40 for L = 600, 50 for L = 660, and 60 for L = 720 (All simulations were performed at the condition of u = 0.03, D = 0.004.).

fluid flow, and thus the shape of the generated concentration gradient at the outlet of the microchannel is found to be described by the parameter. Generated concentration gradients are found to collapse into a universal master curve as a function of the dimensionless parameter (see Figure 8), which can be approximately used to get the concentration at the outlet of the 13937

dx.doi.org/10.1021/ie201136r |Ind. Eng. Chem. Res. 2011, 50, 13932–13939

Industrial & Engineering Chemistry Research

ARTICLE

relative importance of molecule diffusion and fluid flow through the network, which dominates the generation of a concentration gradient in microfluidic networks. This parameter allows one to perform experiments or simulations with reduced models in microchannels and correlate the data to the other flows or geometry sizes.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

Figure 8. Generated concentration gradient at the outlet are found to merge into a correlation curve as a function of the dimensionless parameter, (L/u)/(H2/D), where symbols 0 H = 64, L = 480, u = 0.03; O (red) H = 64, D = 0.004, u = 0.03; 4 (blue) is H = 64, L = 480, D = 0.004. We can find at the outlet of the microchannel, the concentration changes with the change of microchannel’s structure and fluid properties. A reduced parameter can approximately present the concentration change; C2, C4, and C5 are 3 of the pathways at the outlet of the microchannel (see Figure 1b for details).

microchannel of a simulation with the parameter value between 0.012 and 0.063. Therefore, for geometrically similar microchannels, the shape of the generated concentration gradients is only determined by the value of the dimensionless parameter, although these microfluidic systems perhaps having different flow velocities, solute diffusion coefficients, and hydraulic diameters of the microchannels. This is because the dimensionless parameter quantifies the relative importance of the two factors controlling the fluid flowing through the microchannel: molecule diffusion and fluid flow.

4. CONCLUSIONS In this work, we used lattice Boltzmann (LB) method on a 3-dimensional D3Q19 lattice to investigate which factor affects the generation of a concentration gradient in microfluidic networks. Several factors, such as structure characteristics of the branching microfluidic network, flow velocity, solute diffusion coefficient, and the hydraulic diameter of microchannel, were included here to assess their effects on the shape of the generated concentration gradients. We first developed a model of microfluidic networks with hierarchical branches. Our simulation results indicate that the generated solute gradient is controlled by the completion between molecule diffusion and fluid flow through the network. The positions at which micrchannels are split into several daughter microchannels (also their relative positions) are found to exert huge effects on the shape of concentration gradients. We also performed extensive simulation runs to study the dependence of the shape of concentration gradients on the velocity of the flow, properties of fluids, and the length of the microchannel. The shape of the generated concentration gradient at the outlet of the microchannel is found to be controlled by a dimensionless parameter, (L/u)/(H2/D), with L being the length of the microchannel, u being the flow velocity of the flow, D being the solute diffusion coefficient, and H being the hydraulic diameter of the main microchannel. For geometrically similar microchannels, the shape of the generated concentration gradients is only determined by the value of the dimensionless parameter. This is because the dimensionless parameter quantifies the

’ ACKNOWLEDGMENT This work was financially supported by National Natural Science Foundation of China (No. 20736005 and No. 20876004). ’ REFERENCES (1) Seong, G. H.; Crooks, R. M. 2002. Efficient Mixing and Reactions within Microfluidic Channels Using Microbead-Supported Catalysts. J. Am. Chem. Soc. 2002, 124, 13360–13361. (2) Aizenberg, J.; Black, A. J.; Whitesides, G. M. Control of crystal nucleation by patterned self-assembled monolayers. Nature 1999, 398, 495–498. (3) Parent, C. A.; Devreotes, P. N. A Cell’s Sense of Direction. Science 1999, 284, 765–777. (4) Abhyankar, V. V.; Lokuta, M. A.; Huttenlocher, A.; Beebe, D. J. Characterization of a membrane-based gradient generator for use in cellsignaling studies. Lab Chip 2006, 6, 389–393. (5) Jeon, N. L.; Dertinger, S. K. W.; Chiu, D. T.; Choi, I. S.; Stroock, A. D.; Whitesides, G. M. Generation of Solution and Surface Gradients Using Microfluidic Systems. Langmuir 2000, 16, 8311–8316. (6) Dertinger, S. K. W.; Chiu, D. T.; Jeon, N. L.; Whitesides, G. M. Generation of Gradients Having Complex Shapes Using Microfluidic Networks. Anal. Chem. 2001, 73, 1240–1246. (7) Irimia, D.; Geba, D. A.; Toner, M. Universal Microfluidic Gradient Generator. Anal. Chem. 2006, 78, 3472–3477. (8) Li, C.; Chen, R.; Yang, M. Generation of linear and non-linear concentration gradients along microfluidic channel by microtunnel controlled stepwise addition of sample solution. Lab Chip 2007, 7, 1371–1373. (9) Holden, M. A.; Kumar, S.; Castellana, E. T.; Beskok, A.; Cremer, P. S. Generating fixed concentration arrays in a microfluidic device. Sens. Actuators, B 2003, 92, 199–207. (10) Holden, M. A.; Kumar, S.; Beskok, A.; Cremer, P. S. Microfluidic diffusion diluter: bulging of PDMS microchannels under pressure-driven flow. J. Micromech. Microeng. 2003, 13, 412–418. (11) Yamada, M.; Hirano, T.; Yasuda, M.; Seki, M. A microfluidic flow distributor generating stepwise concentrations for high throughput biochemical processing. Lab Chip 2006, 6, 179–184. (12) Campbell, K.; Groisman, A. Generation of complex concentration profiles in microchannels in a logarithmically small number of steps. Lab Chip 2007, 7, 264–272. (13) Kim, C.; Lee, K.; Kim, J. H.; Shin, K. S.; Lee, K.; Kim, T. S.; Kang, J. Y. A serial dilution microfluidic device using a ladder network generating logarithmic or linear concentrations. Lab Chip 2008, 8, 473–479. (14) Lin, F.; Saadi, W.; Rhee, S. W.; Wang, S.; Mittal, S.; Jeon, N. L. Generation of dynamic temporal and spatial concentration gradients using micr ofluidic devices. Lab Chip 2004, 4, 164–167. (15) Sun, K.; Wang, Z.; Jiang, X. Modular microfluidics for gradient generation. Lab Chip 2008, 8, 1536–1543. (16) Erickson, D. Towards numerical prototyping of labs-on-chip: modeling for integrated microfluidic devices. Microfluid. Nanofluid. 2005, 1, 301–318. 13938

dx.doi.org/10.1021/ie201136r |Ind. Eng. Chem. Res. 2011, 50, 13932–13939

Industrial & Engineering Chemistry Research

ARTICLE

(17) Emerson, D. R.; Cieslicki, K.; Gu, X.; Barber, R. W. Biomimetic design of microfluidic manifolds based on a generalized Murray’s law. Lab Chip 2006, 6, 447–454. (18) Barber, R. W.; Emerson, D. R. Optimal design of microfluidic networks using biologically inspired principles. Microfluid. Nanofluid. 2008, 4, 179–191. (19) Qian, Y. H.; D’Humieres, D.; Lattemand, P. Lattice BGK models for Navier-Stokes equation. Europhys. Lett. 1992, 17, 479–484. (20) Bhatnagar, P. L.; Cross, E. P.; Krook, M. A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems. Phys. Rev. 1954, 94, 511–525. (21) Frisch, U.; Hasslacher, B.; Pomeau, Y. Lattice-Gas Automata for Navier-Stokes Equation. Phys. Rev. Lett. 1986, 56, 1505–1508. (22) Guo, Z.; Zhao, T. S. Lattice Boltzmann model for incompressible flows through porous media. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2002, 66, 036304.

13939

dx.doi.org/10.1021/ie201136r |Ind. Eng. Chem. Res. 2011, 50, 13932–13939