CFD Simulation of Droplet Formation in Microchannels by a Modified

5 Mar 2014 - ABSTRACT: A level set method embedded in a Computational Fluidic Dynamics (CFD) simulation is a useful approach in the study of ...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/IECR

CFD Simulation of Droplet Formation in Microchannels by a Modified Level Set Method Wenjie Lan,† Shaowei Li,*,‡ Yujun Wang,§ and Guangsheng Luo*,§ †

State Key Laboratory of Heavy Oil Processing, China University of Petroleum, Beijing 102249, China Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China § Department of Chemical Engineering, Tsinghua University, Beijing 100084, China ‡

ABSTRACT: A level set method embedded in a Computational Fluidic Dynamics (CFD) simulation is a useful approach in the study of microfluidic two-phase flow. In this work, we established a new governing equation for the level set function to overcome the “mass loss” problem. The feasibility of the modified method was proved by the simulation of a droplet formation in a co-tube microchannel. The modified method was then used to simulate the droplet flow in a focusing shape microchannel. The simulation results showed that the droplet size would decrease with an increase in the continuous phase flow rate. In the early stage of the droplet formation process, there was strong inner circulation flow that enhanced the mass transfer between the two phases and generated a concentration notch at the head of the droplet. In the later stage of the droplet formation process, the droplet moving velocity became higher, and there was no longer inner circulation flow in the droplet. The mass transfer speed slowed, and there was an enrichment of the material at the head of the droplet. The total mass transfer ratio in the droplet formation process reached 50%. The simulation results were consistent with both our previous experiment results and the results reported in the literature, which further proved the feasibility of the modified method.

1. INTRODUCTION Recently, the rapid development of microreaction and micromixing technology has led to a considerable variety of microfluidic devices. These devices have been used in many fields for their high efficiency, safety, repeatability, and facile controllability.1,2 In particular, significant advances have been made in the use of microfluidic devices for controlling twophase flow patterns. Many investigations have been carried out to analyze how the operation conditions and fluid properties, as well as channel structures, affect multiphase flow.2−7 Well controlled two-phase flows in microdevices have the advantages of a large interfacial area, short transfer distance, and fast mixing performance. This can reduce mass transfer limitations to achieve better performances relative to conventional scale systems.2 Thus, the formation of microdroplet in microdevices has a very important role in applications of chemical reaction,8−10 liquid−liquid extraction,11−13 biological analysis,14 crystallization,15 polymer synthesis,16−19 structural material preparation,20−22 and nanoparticle synthesis.6,23−26 However, many underlying mechanisms for droplet formation are still poorly understood, even though the droplet size and droplet size distribution are generally regarded as important product properties.27 In order to realize the precise control of the droplet flow, a better understanding of the complex droplet formation behavior in microfluidic devices is needed. It is difficult to obtain a deep understanding just relying on experimental investigations because there are too many variables, such as flow speed, viscosity, capillary force, geometry, and the wetting property of the microchannel wall. Moreover, measurement in microfluidic systems is also a difficult task. It is hard to obtain the detail of the flow field in a microfluidic system. Because of these two reasons, the current investigation on microfluidic system is often empirical. © 2014 American Chemical Society

Therefore, it is imperative to develop a numerical simulation method that provides an optimal alternative solution to such a complex problem. Compared to the experimental studies, there are relatively few investigations that use numerical techniques to dissect twophase flows in microfluidic devices. Existing numerical methods used to solve two-phase flow problems include the fronttracking method,28−30 boundary integral method,31,32 finite element method,33,34 volume of fluid method,35−38 phase field method,39,40 Lattice Boltzmann method,41−43 and level set method.44−51 All of the above methods have their advantages and disadvantages. For example, the volume of fluid method does not require special procedures to model topological changes of the front. However, it is difficult to calculate the curvature of the front from volume fractions. A high-order scheme needs to be introduced for reconstructing a highly curved surface.52,53 The level set method is one of the most powerful techniques available to determine implicit surfaces within a fixed grid system.52 The strength of the level set method lies in its ability to accurately compute two-phase flows with complex topological changes.45−49 Moreover, this method is excellent in accurately computing problems with surface tension.53 However, with the classical governing equation of the level set function, the level set function may be distorted by the flow field after some steps of calculation. Therefore, the gradient of the level set function may become too large or too small. This inadequacy can lead to an inaccurate approximation of the Received: Revised: Accepted: Published: 4913

September 16, 2013 February 27, 2014 March 5, 2014 March 5, 2014 dx.doi.org/10.1021/ie403060w | Ind. Eng. Chem. Res. 2014, 53, 4913−4921

Industrial & Engineering Chemistry Research

Article

In the level set method, the level set function ϕ is usually defined as the signed distance function to the interface. Then n and curvature κ are expressed as the unit normal vector ⇀ ⇀ n = ∇ϕ (3)

normal, the mean curvature, and the Gaussian curvature of the interface; therefore, a “mass loss” problem may occur. This is an inherent problem of the classical level set method, and increasing the computing accuracy will not solve this issue. A commonly used method to overcome this problem is the reinitialization of the level set function.54−60 Generally, the level set function is reinitialized in order to maintain it as a true signed distance function, and so the method is also known as redistancing. Besides that, a lot of new methods and technologies to deal with this problem have been reported in recent years. For example, the Discontinuous Galerkin (DG) method promotes good mass conservation by a DG implementation of the level set transport equation.61−64 Olsson et al.65,66 reduced the mass conservation errors by replacing the usual signed distance function with a hyperbolic tangent function. Some other researchers reduced the mass loss of the level set method by coupling another method with a good mass conservation property, such as the volume of fluid method,67−72 front-tracking method,73,74 or marker particles method.75 A gradient-augmented level set method, which augmented the level set function values by gradient information, was also developed, and the mass loss issue was addressed successfully.76,77 Though approaches to control or limit the mass loss of the level set method have been studied for decades, the research is ongoing in this topic, and improvements and alternatives are much needed. To solve the mass loss problem and simultaneously maintain the simplicity of the level set method, we present a modified method to avoid the distortion of the level set function by introducing a new governing equation in this work. In this method, the reinitializing process that requires extra iterative computing is no longer needed. The new governing equation itself can keep the level set function a signed distance function. The feasibility of the modified method was first proved by the simulation of droplet formation in a cotube microchannel. Then, the two-phase flow in a focusing shape microchannel, which is one of the most important microchannel types for droplet formation,78−80 was simulated with this modified method. The “mass loss” problem was indeed overcome. The simulation results with the modified method were compared with the experimental results, and a good agreement was observed.

and κ = −∇·⇀ n = −∇·(∇ϕ)

The force term in eq 2 is normal to the interface and in proportion to the interface curvature and is expressed as λδ ⇀ fS = − ∇·(∇ϕ)∇ϕ (5) ε Herein, ε is the thickness of the interface, and δ is the Dirac delta function, which is defined as

⎧1, |ϕ| ≤ ε δ=⎨ ⎩ 0, |ϕ| > ε ⎪



ρ

D⇀ u ⇀ FM − ∇p + μ∇2 ⇀ u + fS = ρ⇀ Dt

(6)

The governing equation for the level set function in the classical level set simulation is written as

Dϕ =0 Dt

(7)

Here is the so-called material derivative or Lagrangian derivative, which means that the ϕ value for a fluid group does not change with the flow development. Equation 7 is tenable at the interface because the ϕ value at the interface is always equal to zero. However, the equation may not be tenable out of the interface. As shown in Figure 1, the distance from the

Figure 1. Distortion of the level set function.

interface (ϕ = 0) to the ϕ = ϕ0 surface is ϕ0 at time t = 0. After a duration of t1, the ϕ = ϕ0 surface is no longer an isometric surface of the interface because of the distortion of the level set function. In other words, ϕ is no longer the distance function to the interface at t = t1. This will create an error to the calculation ⇀ of fS from eq 5 and will cause the “mass loss” problem. We will show the influence of “mass loss” on the simulation results in Section 4. We try to overcome the “mass loss” problem by introducing a new governing equation that is tenable for fluid groups both at and outside of the interface. With this new governing equation, ϕ should be a real distance function from the interface. For this purpose, moving coordinates, based on the level set function, are established. Geodesic coordinates on the two principal directions of the isosurface for the level set function is first constructed. As shown in Figure 2, we denote the unit vectors of these two directions as ⇀ e1 and ⇀ e 2 . Then, we set the third coordinate on the gradient direction of the level set function with the coordinate value x3 = ϕ.

2. SIMULATION METHOD 2.1. Theory Model. As a base of the computational fluidic dynamics simulation, the governing equations for conservation laws of mass and momentum should be considered, which could be written in the following forms with the assumption that the fluid is incompressible ∇·⇀ u =0

(4)

(1)

(2)

In these equations, the density and viscosity of the fluid are considered to be variables that allow the simulation of multiphase systems and give rise to sharp changes in density ⇀ and viscosity at the interface. fS is a force term introduced by the interface, and its expression is dependent on the interface describing method. 4914

dx.doi.org/10.1021/ie403060w | Ind. Eng. Chem. Res. 2014, 53, 4913−4921

Industrial & Engineering Chemistry Research

Article

⎛ ∂C ⎞ ⎛ ∂C ⎞ ⎜DA ⎟ = ⎜DA ⎟ (flux continuity at the interface) ⎝ ∂ϕ ⎠c ⎝ ∂ϕ ⎠d (13)

(C)c = K (C)d (interfacial equilibrium)

Here, K is the equilibrium constant. The subscripts c and d denote the continuous phase and the dispersed phase respectively. 2.2. Geometry and Modeling Method. To prove the feasibility of the modified level set method, the method was used to simulate two-phase flow in a co-tube microchannel as shown in Figure 3(a). The simulation by the original level set

Figure 2. Moving coordinates based on the level set function.

The velocity components on the three coordinate directions are u1, u2, and u3, respectively. Consider a fluid group FG1, of which the level set function value is ϕ̃ at time t. Then the third coordinate value of FG1 is x3 = ϕ ̃ at time t. After an infinitesimal time dt, the third coordinate value of FG1 becomes x3 = ϕ ̃ + u3dt . At the same time, the third coordinate value of the corresponding fluid group on the interface (has the same value of x1 and x2 with FG1) becomes u3dt |ϕ= 0 . Therefore, the distance of FG1 to the interface at time t + dt is ϕ|t + dt = ϕ ̃ + u3dt − (u3dt )|ϕ = 0 . On the basis of the

Figure 3. Geometry of the microchannels: (a) Co-tube microchannel. (b) Focusing shape microchannel.

method was also implemented for comparison. The size of the inner tube is 200 μm diameter × 35 μm thick. The diameter of the outer tube is 400 μm. The dispersed phase flows into the inner tube while the continuous phase flows into the annular space. The dispersed phase is sheared by the continuous phase to form droplets at the outlet of the inner tube. The modified level set method was then applied to the simulation of the flow and mass transfer in a focusing shape microchannel, which has been used in our previous experiment for investigating the flow pattern and mass transfer.81 The geometry of the focusing shape channel is shown in Figure 3(b). The branch channels on both sides are the inlets of the continuous phase, while that in the middle is the inlet of the dispersed phase. The size of all the branch channels is 150 μm wide × 100 μm high. After the dispersed phase is sheared by the continuous phase to form droplets, the droplets will flow into a bigger channel that is 300 μm wide × 250 μm high. The simulation results were compared with our previous experimental results.81 The mesh for the geometry shown in Figure 3 was created by CFX Build. As shown in Figure 4, the mesh consisted of 0.528

definition of material derivative, we can get the following equation ϕ| −ϕ|t Dϕ = t + dt = u3 − u3|ϕ = 0 Dt dt

(8)

This equation is based on the moving coordinates and will result in complex calculation. To simplify calculation complexity, eq 8 should be translated to a more universal form. We can easily get the following equations u3 = ⇀ u ·⇀ e3

⇀ e 3 = ∇ϕ

(9) (10)

By submitting these two equations to eq 8, we have Dϕ ⇀ = u ·∇ϕ − (⇀ u ·∇ϕ)|ϕ= 0 Dt

(11)

This is an expression independent of coordinate selection. Though we obtain eq 11 in a specific coordinate system, it can be used in any coordinate system. Equation 11 is just the governing equation of the level set function, which is tenable for the whole fluid field. Comparing eqs 11 and 7, we can see there are two additional terms on the right-hand side of eq 11. The two additional terms can be explained as a compensation for the distortion of the level set function. Unlike eq 7, eq 11 will not bring error to the level set function if the computing accuracy is high enough. Furthermore, to simulate the mass transfer between the two phases, a mass transfer equation is needed, which is expressed as follows

DC = DA ∇2 C Dt

(14)

Figure 4. Tetrahedral mesh used for CFX simulation.

M unstructured tetrahedral cells for the co-tube microchannel and 0.414 M cells for the focusing shape microchannel. The mesh size was 12 μm for both of the microchannels. Mesh dependency was not explored but can be assumed to be similar to other comparable methods. A transient-state simulation was carried out by using the commercial software CFX 5.6. A time step dependence study was performed to determine the time step size. A time step size of the order of the physical time scale (∼10−4 s) was initially used. Then double and half this time step size were used to rerun the simulation for a same period.

(12)

where C is the concentration of a tested material, and DA is the diffuse coefficient of the material. Two interfacial conditions are needed to calculate the mass transfer process. These two conditions can be expressed as 4915

dx.doi.org/10.1021/ie403060w | Ind. Eng. Chem. Res. 2014, 53, 4913−4921

Industrial & Engineering Chemistry Research

Article

Figure 5. Comparison between the modified and original level set method. (a) and (b): Simulation results of the original level set method. (c), (d), and (e): Results of the modified level set method.

Figure 6. Simulation result for the droplet formation process, uc = 0.088 m/s and ud = 0.04 m/s.

Repeated tests indicated that a time step of 10−5 s is an appropriate value in our simulation, in which good convergence was achieved. The maximum number of coefficient loops for each time step was set as 10. The convergence criterion was set as that the normalized root-mean-square of the equation residuals reaches 0.00001. CFX’s high resolution discretization for both the advection and turbulence terms were used for the calculation.

broken into droplet. Figure 5(b) is the final station of the simulation, and it is also a steady state because the “mass loss” balances the dispersed phase feed. Therefore, we conclude that the original level set method cannot be directly used to simulate droplet formation in the microchannels. In comparison, Figure 5(c), (d), and (e) show a successful simulation of the droplet formation by using the modified level set method. We can also observe that the level set function maintains a distance function to the interface in the duration of the droplet formation. Though some small deviation is caused by the calculating error, the droplet formation process is successfully simulated. The results prove that the modified level set method is feasible for the simulation of the droplet formation in a microchannel. 4.2. Droplet Formation and Mass Transfer in Focusing Shape Microchannel. Figures 6 and 7 show simulation results for the droplet formation and mass transfer processes when the flow velocity of the continuous phase is 0.088 m/s and the flow velocity of the dispersed phase is 0.04 m/s. In comparison, Figures 8 and 9 show the results when the flow velocity of the continuous phase is changed to 0.132 m/s. The flow velocities here are all apparent velocities, which can be

4. RESULTS AND DISCUSSION 4.1. Comparison between Modified and Original Level Set Methods. Figure 5 is the cross section of the simulation result in the co-tube microchannel. The level set function ϕ changes from 0 (blue) to 0.04 mm (red) in this figure. Figure 5(a) and (b) are the simulation results of the original level set method, while panels (c), (d), and (e) are the results of the modified level set method. We can see from Figure 5(a) and (b) that the distance from the ϕ = 0.04 mm surface to the ϕ = 0 surface is not 0.04 mm at the front of the forming droplet. The level set function cannot maintain its feature (a distance function from the interface) in the original level set method. As a result, the dispersed phase cannot be 4916

dx.doi.org/10.1021/ie403060w | Ind. Eng. Chem. Res. 2014, 53, 4913−4921

Industrial & Engineering Chemistry Research

Article

Figure 7. Simulation result for the mass transfer in the droplet formation process, uc = 0.088 m/s and ud = 0.04 m/s.

Figure 8. Simulation result for the droplet formation process, uc = 0.132 m/s and ud = 0.04 m/s.

obtained by calculating the ratio of the flow rate and the channel cross area. Comparing Figures 6 and 8, we see that the droplet size decreases with an increase in the continuous phase flow rate. An increase in the continuous phase flow rate can also result in an increase in the distance from the break point to the entrance. High flow rate of the continuous phase leads to a long dispersed phase neck. This simulation result is consistent with our previous experiment.81 From Figures 7 and 9, we see that there is a concentration notch at the head of the droplet during the droplet formation. The reason is that the continuous phase flows around the droplet in high speed and takes away the solute transferred out from the droplet continuously, which causes a high mass transfer speed at the droplet surface and creates the lowest concentration at the head of the droplet. Marry et al.82 observed a similar phenomenon in their experiment characterizing mass transfer with the fluorescent labeling technique. The simulation results are consistent with

the experiment observation. Otherwise, we see from Figures 7(d) and 9(d) that the concentration distribution in the droplet experiences a translation when the droplet is broken off. The concentration at the head of the droplet changes from the lowest to the highest. The material in the droplet moves from the rear to the head after the droplet is separated. This phenomenon is caused by the changes in the inner circulation flow when the droplet is disintegrated. In Marry’s experiment82 as well as our experiment,81 a similar phenomenon was observed. 4.3. Effect of Inner Circulation on Mass Transfer. To further investigate the mass transfer ratio in the droplet formation process, we calculated the integral of the concentration distribution in the droplet. The result was compared with the total mass of the material in the droplet, and the mass transfer ratio was obtained, which was denoted as y. Figure 10 shows the mass transfer ratio variation with time. We see from the figure that a higher continuous phase flow rate 4917

dx.doi.org/10.1021/ie403060w | Ind. Eng. Chem. Res. 2014, 53, 4913−4921

Industrial & Engineering Chemistry Research

Article

Figure 9. Simulation result for the mass transfer in the droplet formation process, uc = 0.132 m/s and ud = 0.04 m/s.

comparison, the highest mass transfer speed is reached at the side flanks of the droplet because of the high convection effect of the continuous phase. The flow field also explained why the mass transfer is high at the early stage of the droplet formation process. The direct reason is the high speed of the continuous phase flow velocity at the droplet surface and the inner circulation flow in the droplet. Figure 11(c) and (d) show the flow field image for the later stage of the droplet formation process and its magnification, corresponding to Figure 6(c) and 8(c). We see that the moving velocity of the droplet increase, and there is no vortex in the droplet anymore. The reason is that the shearing force of the continuous phase on the droplet is bigger than the interfacial force. Thus, the droplet moves in a relative velocity similar to the continuous phase, and the inner vortex of the droplet disappears. Moreover, we see that the flow velocity at the rear of the droplet is higher than that at the front of the droplet. This is why the material in the droplet moves from the rear to the head after the droplet is broken off. In fact, this moving effect starts before the droplet is broken off. Because there is little convection effect of the continuous phase and no inner circulation flow in the droplet, the mass transfer in the later stage of the droplet formation process is mainly contributed by diffusion and thus slows down.

Figure 10. Mass transfer ratio varying with time in the droplet formation process.

leads to a higher mass transfer ratio at the droplet formation stage. The mass transfer ratio at the droplet formation stage can reach 50%, which is consistent with our previous experiment.81 Moreover, we see a new result that the mass transfer speed is high at the early stage of the droplet formation and decreases at the later stage of the droplet formation. This phenomenon has not yet been reported so far because of the measurement technique limitation. This can be explained by the inner circulation flow, which we will present later. Figure 11 shows the simulated flow field of the system, which will help us to further understand the mass transfer process in the droplet formation process. Figure 11(a) and (b) are the flow field image for the early stage of the droplet formation process and its magnification, corresponding to Figures 6(a) and 8(a). We see that the flow velocity of the continuous phase at the droplet surface is very high, and the continuous phase has an obvious shearing effect on the droplet. Thus, there are two obvious vortexes in the droplet. The vortexes take the low concentration fluid to the front of the droplet and cause the formation of a notch in the concentration profile. This proved that the formation of the concentration notch is not only because of the high mass transfer speed at the front. In

5. CONCLUSIONS We established a new governing equation for the level set method by using moving coordinates and differential analysis on a curved surface. Compared to the original governing equation, the new governing equation is valid in the whole flow field. The modified method can overcome the “mass loss” problem and avoid the reinitialization process. Comparison between the original method and the modified method is carried out in the simulation of two-phase flow in a co-tube microchannel. The droplet formation process cannot be simulated by the original method because of the “mass loss” problem, while it can be successfully simulated by the modified method. 4918

dx.doi.org/10.1021/ie403060w | Ind. Eng. Chem. Res. 2014, 53, 4913−4921

Industrial & Engineering Chemistry Research

Article

Figure 11. Simulation result for the flow field of the system in the droplet formation process, uc = 0.132 m/s and ud = 0.04 m/s. (a) Flow field of the system at time t = 0.5 ms. (b) Flow field in the droplet. (c) Fow field of the system at time t = 1.5 ms. (d) Flow field in the droplet.

The modified method was used to simulate the droplet flow in the focusing shape microchannel. The simulation results showed that an increase in the continuous phase flow rate led to a decrease in the droplet size and an increase in the distance from the droplet break point to the entrance. In the early stage of the droplet formation process, there was strong inner circulation flow because of the shearing effect of the continuous phase. This enhanced the mass transfer between the two phases and caused the formation of a concentration notch at the head of the droplet. In the later stage of the droplet formation process, the droplet moving velocity significantly increased, and there was no longer inner circulation flow in the droplet. The mass transfer speed slowed, and there was an enrichment of the material at the head of the droplet because the velocity at the rear of the droplet was higher than that in the front. The total mass transfer ratio in the droplet formation process could reach 50%. The simulation results were consistent with our previous experiment results and the results reported in the literature, which further proved the feasibility of the modified method.



(2) Günther, A.; Jensen, K. F. Multiphase microfluidics: From flow characteristics to chemical and materials synthesis. Lab Chip 2006, 6, 1487. (3) Xu, J. H.; Luo, G. S.; Li, S. W.; Chen, G. G. Shear force induced monodisperse droplet formation in a microfluidic device by controlling wetting properties. Lab Chip 2006, 6, 131. (4) Xu, J. H.; Li, S. W.; Wang, Y. J.; Luo, G. S. Controllable gas-liquid phase flow patterns and monodisperse microbubbles in a microfluidic T-junction device. Appl. Phys. Lett. 2006, 88, 133506. (5) Guillot, P.; Colin, A. Stability of parallel flows in a microchannel after a T junction. Phys. Rev. E 2005, 72, 066301. (6) Chan, E. M.; Alivisatos, A. P.; Mathies, R. A. High-temperature microfluidic synthesis of CdSe nanocrystals in nanoliter droplets. J. Am. Chem. Soc. 2005, 127, 13854. (7) Zheng, B.; Roach, L. S.; Ismagilov, R. F. Screening of protein crystallization conditions on a microfluidic chip using nanoliter-size droplets. J. Am. Chem. Soc. 2003, 125, 11170. (8) Song, H.; Chen, D. L.; Ismagilov, R. F. Reactions in droplets in microfluidic channels. Angew. Chem., Int. Ed. 2006, 45, 7336. (9) Burns, J. R.; Ramshaw, C. The intensification of rapid reactions in multiphase systems using slug flow in capillaries. Lab Chip 2001, 1, 10. (10) Burns, J. R.; Ramshaw, C. A microreactor for the nitration of benzene and toluene. Chem. Eng. Commun. 2002, 189, 1611. (11) Dittrich, P. S.; Tachikawa, K.; Manz, A. Micro toal analysis systems: Latest advancements and trends. Anal. Chem. 2006, 78, 3887. (12) Xu, J. H.; Luo, G. S.; Chen, G. G.; Tan, B. Mass transfer performance and two-phase flow characteristic in membrane dispersion mini-extractor. J. Membr. Sci. 2005, 249, 75. (13) Kumemura, M.; Korenaga, T. Quantitative extraction using flowing nano-liter droplet in microfluidic system. Anal. Chim. Acta 2006, 558, 75. (14) Grodrian, A.; Metze, J.; Henkel, T.; Martin, K.; Roth, M.; Kohler, J. M. Segmented flow generation by chip reactors for highly parallelized cell cultivation. Biosens. Bioelectron. 2004, 19, 1421. (15) Zheng, B.; Tice, J. D.; Ismagilov, R. F. Formation of droplets of in microfluidic channels alternating composition and applications to indexing of concentrations in droplet-based assays. Anal. Chem. 2004, 76, 4977. (16) Steinbacher, J. L.; McQuade, D. T. Polymer chemistry in flow: New polymers, beads, capsules, and fibers. J. Polym. Sci., Part A: Polym. Chem. 2006, 44, 6505.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (S.L.). *E-mail: [email protected] (G.L.). Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS



REFERENCES

This work was supported by the National Natural Science Foundation of China (21106077, 21376132) and the Science Foundation of China University of Petroleum, Beijing (2462013YJRC025).

(1) Ehrfeld, W.; Hessel, V.; Löwe, H. Microreactors: New Technology for Modern Chemistry.; Wiley-VCH, Weinheim, Germany, 2000. 4919

dx.doi.org/10.1021/ie403060w | Ind. Eng. Chem. Res. 2014, 53, 4913−4921

Industrial & Engineering Chemistry Research

Article

(17) Zourob, M.; Mohr, S.; Mayes, A. G.; Macaskill, A.; Pérez-Moral, N.; Fielden, P. R.; Goddard, N. J. A micro-reactor for preparing uniform molecularly imprinted polymer beads. Lab Chip 2006, 6, 296. (18) Nie, Z.; Xu, S.; Seo, M.; Lewis, P. C.; Kumacheva, E. Polymer particles with various shapes and morphologies produced in continuous microfluidic reactors. J. Am. Chem. Soc. 2005, 127, 8058. (19) Quevedo, E.; Steinbacher, J.; McQuade, D. T. Interfacial polymerization within a simplified microfluidic device: capturing capsules. J. Am. Chem. Soc. 2005, 127, 10498. (20) Dendukuri, D.; Tsoi, K.; Hatton, T. A.; Doyle, P. S. Controlled synthesis of nonspherical microparticles using microfluidics. Langmuir 2005, 21, 2113. (21) Dendukuri, D.; Pregibon, D. C.; Collins, J.; Hatton, T. A.; Doyle, P. S. Continuous-flow lithography for high-throughput microparticle synthesis. Nat. Mater. 2006, 5, 365. (22) Khan, S. A.; Günther, A.; Schmidt, M. A.; Jensen, K. F. Microfluidic synthesis of colloidal silica. Langmuir 2004, 20, 8604. (23) Yen, B. K. H.; Günther, A.; Schmidt, M. A.; Jensen, K. F.; Bawendi, M. G. A microfabricated gas-liquid segmented flow reactor for high-temperature synthesis: The case of CdSe quantum dots. Angew. Chem., Int. Ed. 2005, 44, 5447. (24) Hung, L. H.; Tseng, W. Y.; Choi, K.; Tan, Y. C.; Shea, K. J.; Lee, A. P. Controlled Droplet Fusion in Microfluidic Devices. In 8th International Conference on MicroTAS, Malmo, Sweden, 2004, Vol. 2, p 539. (25) Sotowa, K. I.; Irie, K.; Fukumori, T.; Kusakabe, K.; Sugiyama, S. Droplet formation by the collision of two aqueous solutions in a microchannel and application to particle synthesis. Chem. Eng. Technol. 2007, 30, 383. (26) Shestopalov, I.; Tice, J. D.; Ismagilov, R. F. Multi-step synthesis of nanoparticles performed on millisecond time scale in a microfluidic droplet-based system. Lab Chip 2004, 4, 316. (27) McClements, D. J.; Chanamai, R. Physicochemical properties of monodisperse oil-in-water emulsions. J. Dispersion Sci. Technol. 2002, 23, 125. (28) Unverdi, S. O.; Tryggvason, G. A front-tracking method for viscous, incompressible, multi-fluid flows. J. Comput. Phys. 1992, 100, 25. (29) Hua, J.; Lou, J. Numerical simulation of bubble rising in viscous liquid. J. Comput. Phys. 2007, 222, 769. (30) Hua, J.; Zhang, B.; Lou, J. Numerical simulation of microdroplet formation in coflowing immiscible liquids. AIChE J. 2007, 53, 2534. (31) Yeo, L. Y.; Matar, O. K.; de Ortiz, E. S. P.; Hewitt, G. E. Film drainage between two surfactant-coated drops colliding at constant approach velocity. J. Colloid Interface Sci. 2003, 257, 93. (32) Pozrikidis, C. Expansion of a two-dimensional foam. Eng. Anal. Boundary Elem. 2002, 26, 495. (33) Notz, P. K.; Chen, A. U.; Basaran, O. A. Satellite drops: Unexpected dynamics and change of scaling during pinch-off. Phys. Fluids 2001, 13, 549. (34) Zhou, C.; Yue, P.; Feng, J. J. Formation of simple and compound drops in microfluidic devices. Phys. Fluids 2006, 18, 092105. (35) Weber, M. W.; Shandas, R. Computational fluid dynamics analysis of microbubble formation in microfluidic flow-focusing devices. Microfluid. Nanofluid. 2007, 3, 195. (36) Davidson, M. R.; Harvie, D. J. E.; Cooper-White, J. J. Flow focusing in microchannels. ANZIAM J. 2005, 46(E), C47. (37) Abrahamse, A. J.; van der Padt, A.; Boom, R. M.; de Heij, W. B. C. Process fundamentals of membrane emulsification: Simulation with CFD. AIChE J. 2001, 47, 1285. (38) Ohta, M.; Yamamoto, M.; Suzuki, M. Numerical-analysis of a single drop formation process under pressure pulse condition. Chem. Eng. Sci. 1995, 50, 2923. (39) De Menech, M. Modeling of droplet breakup in a microfluidic T-shaped junction with a phase-field model. Phys. Rev. E 2006, 73, 031505.

(40) Kuksenok, O.; Jasnow, D.; Yeomans, J.; Balazs, A. C. Periodic droplet formation in chemically patterned microchannels. Phys. Rev. Lett. 2003, 91, 108303. (41) Dupin, M. M.; Halliday, I.; Care, C. M. Simulation of a microfluidic flow-focusing device. Phys. Rev. E 2006, 73, 055701(R). (42) van der Graaf, S.; Nisisako, T.; Schroen, C. G. P. H.; van der Sman, R. G. M.; Boom, R. M. Lattice Boltzmann simulations of droplet formation in a T-shaped microchannel. Langmuir 2006, 22, 4144. (43) Yu, Z.; Heraminger, O.; Fan, L. S. Experiment and lattice Boltzmann simulation of two-phase gas-liquid flows in microchannels. Chem. Eng. Sci. 2007, 62, 7172. (44) Osher, S.; Sethian, J. A. Fronts propagating with curvaturedependent speed: Algorithms based on Hamilton-Jacobi formulations. J. Comput. Phys. 1988, 79, 12. (45) Sussman, M.; Smereka, P.; Osher, S. A level set approach for computing solutions to incompressible two-phase flow. J. Comput. Phys. 1994, 114, 146. (46) Chang, Y. C.; Hou, T. Y.; Merriman, B.; Osher, S. A level set formulation of Eulerian interface capturing methods for incompressible fluid flows. J. Comput. Phys. 1996, 124, 449. (47) Fedkiw, R. P.; Aslam, T.; Merriman, B.; Osher, S. A nonoscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method). J. Comput. Phys. 1999, 152, 457. (48) Kang, M.; Fedkiw, R. P.; Liu, X. D. A boundary condition capturing method for multiphase incompressible flow. J. Sci. Comput. 2000, 15, 323. (49) Sussman, M.; Puckett, E. G. A coupled level set and volume-offluid method for computing 3D and axisymmetric incompressible twophase flows. J. Comput. Phys. 2000, 162, 301. (50) Tatineni, M.; Zhong, X. L. Numerical simulation of unsteady low-Reynolds-number separated flows over airfoils. AIAA J. 2000, 38, 1295. (51) Cubaud, T.; Tatineni, M.; Zhong, X. L.; Ho, C. M. Bubble dispenser in microfluidic devices. Phys. Rev. E 2005, 72, 037302. (52) Watanabe, Y.; Saruwatari, A.; Ingram, D. M. Free-surface flows under impacting droplets. J. Comput. Phys. 2008, 227, 2344. (53) Sussman, M.; Almgren, A. S.; Bell, J. B.; Colella, P.; Howell, L. H.; Welcome, M. L. An adaptive level set approach for incompressible two-phase flows. J. Comput. Phys. 1999, 148, 81. (54) Adalsteinsson, D.; Sethian, J. A. A fast level set method for propagating interfaces. J. Comput. Phys. 1995, 118, 269. (55) Dervieux, A.; Thomasset, F. A finite element method for the simulation of Rayleigh−Taylor instability. Lect. Notes Math. 1979, 771, 145. (56) Dervieux, A.; Thomasset, F. Multifluid incompressible flows by a finite element method. Lect. Notes Phys. 1981, 11, 158. (57) Helmsen, J.; Puckett, E.; Colella, P.; Dorr, M. Two new methods for simulating photolithography development in 3D. Proc. SPIE 1996, 2726, 253−261. (58) Tsai, Y. H.; Cheng, L. T.; Osher, S.; Zhao, H. K. Fast sweeping algorithms for a class of Hamilton−Jacobi equations. SIAM J. Numer. Anal. 2003, 41, 673. (59) Zhao, H. K. A fast sweeping method for eikonal equations. Math. Comput. 2005, 74, 603. (60) Cheng, L. T.; Tsai, Y. H. Redistancing by flow of time dependent eikonal equation. J. Comput. Phys. 2008, 227, 4002. (61) Marchandise, E.; Remacle, J.-F. A stabilized finite element method using a discontinuous level set approach for solving two phase incompressible flows. J. Comput. Phys. 2006, 219, 780. (62) Marchandise, E.; Remacle, J.-F.; Cheveaugeon, N. A quadraturefree discontinuous Galerkin method for the level set equation. J. Comput. Phys. 2006, 212, 338. (63) Pochet, F.; Hillewaert, K.; Geuzaine, P.; Remacle, J.-F.; Marchandise, É . A 3D strongly coupled implicit discontinuous Galerkin level set-based method for modeling two-phase flows. Comput. Fluids 2013, 87, 144. 4920

dx.doi.org/10.1021/ie403060w | Ind. Eng. Chem. Res. 2014, 53, 4913−4921

Industrial & Engineering Chemistry Research

Article

(64) Owkes, M.; Desjardins, O. A discontinuous Galerkin conservative level set scheme for interface capturing in multiphase flows. J. Comput. Phys. 2013, 249, 275. (65) Olsson, E.; Kreiss, G. A conservative level set method for two phase flow. J. Comput. Phys. 2005, 210, 225. (66) Olsson, E.; Kreiss, G.; Zahedi, S. A conservative level set method for two phase flow II. J. Comput. Phys. 2007, 225, 785. (67) Sussman, M.; Puckett, E. G. A coupled level set and volume-offluid method for computing 3D and axisymmetric incompressible twophase flows. J. Comput. Phys. 2000, 162, 301. (68) Sussman, M.; Smith, K. M.; Hussaini, M. Y.; Ohta, M.; Zhi-Wei, R. A sharp interface method for incompressible two-phase flows. J. Comput. Phys. 2007, 221, 469. (69) Sussman, M. A second order coupled level set and volume-offluid method for computing growth and collapse of vapor bubbles. J. Comput. Phys. 2003, 187, 110. (70) Van der Pijl, S. P.; Segal, A.; Vuik, C.; Wesseling, P. A massconserving level set method for modeling of multi-phase flows. Int. J. Numer. Methods Fluids 2005, 47, 339. (71) Sun, D. L.; Tao, W. Q. A coupled volume-of-fluid and level set (VOSET) method for computing incompressible two-phase flows. Int. J. Heat Mass Transfer 2010, 53, 645. (72) Wang, T.; Li, H.; Feng, Y.; Shi, D. A coupled volume-of-fluid and level set (VOSET) method on dynamically adaptive quadtree grids. Int. J. Heat Mass Transfer 2013, 67, 70. (73) Shin, S.; Yoon, I.; Juric, D. The local front reconstruction method for direct simulation of two- and three- dimensional multiphase flows. J. Comput. Phys. 2011, 230, 6605. (74) Basting, S.; Weismann, M. A hybrid level set−front tracking finite element approach for fluid−structure interaction and two-phase flow applications. J. Comput. Phys. 2013, 255, 228. (75) Enright, D.; Fedkiw, R.; Ferziger, J.; Mitchell, I. A hybrid particle level set method for improved interface capturing. J. Comput. Phys. 2002, 183, 83. (76) Nave, J. C.; Rosales, R. R.; Seibold, B. A gradient-augmented level set method with an optimally local, coherent advection scheme. J. Comput. Phys. 2010, 229, 3802. (77) Anumolu, L.; Trujillo, M. F. Gradient augmented reinitialization scheme for the level set method. Int. J. Numer Methods Fluids 2013, 73, 1011. (78) Gañań -Calvo, A. M. Generation of steady liquid microthreads and micron-sized monodisperse sprays in gas streams. Phys. Rev. Lett. 1998, 80, 285. (79) Anna, S. L.; Bontoux, N.; Stone, H. A. Formation of dispersions using “flow focusing” in microchannels. App. Phys. Lett. 2003, 82, 364. (80) Ward, T.; Faivre, M.; Abkarian, M.; Stone, H. A. Microfluidic flow focusing: Drop size and scaling in pressure versus flow-rate-driven pumping. Electrophoresis 2005, 26, 3716. (81) Xu, J. H.; Tan, J.; Li, S. W.; Luo, G. S. Enhancement of mass transfer performance of liquid−liquid system by droplet flow in microchannels. Chem. Eng. J 2008, 141, 242. (82) Mary, P.; Studer, V.; Tabeling, P. Microfluidic droplet-based liquid-liquid extraction. Anal. Chem. 2008, 80, 2680.

4921

dx.doi.org/10.1021/ie403060w | Ind. Eng. Chem. Res. 2014, 53, 4913−4921