Fluid mixing in droplet-based microfluidics with a serpentine microchannel

Abstract

The hydrodynamics and mixing process within droplets travelling along a three dimensional serpentine microchannel are studied using a computational fluid dynamics simulation based on the volume-of-fluid approach. The fluid mixing within the droplet follows symmetric circulations in the straight section, which generates axial mixing. In the winding section, the asymmetric circulations lead to the reorientation of the fluids within the droplet, thus enhancing the mixing efficiency. The mixing performance is controlled by the spatial distribution of the mixing components and the circulation period within the droplet. The best mixing occurs when the droplet size is comparable with the channel width. When the droplet size is less than two times the channel width, the asymmetric circulations make it easy for the fluid to distribute in the axial direction, which leads to a fast mixing process. For larger droplets, the long circulation period becomes more significant, which causes lower mixing efficiency.

---

Introduction

Liquid inter-mixing in microfluidics has attracted wide interest in fundamental science and engineering practices because of the enhanced mixing efficiency and processing capability, which show great potential for applications in areas such as pharmaceutics, medical instruments, chemical analysis and biotechnology.^1–3 The most frequently used microfluidic system for mixing purposes is droplet flows because they eliminate the dispersion effect that is often associated with a single phase flow and the fluids are mixed faster as well.⁴ Compared to the conventional macro-scale mixing, fluid mixing in droplet flows is more uniform but takes a shorter time. They are especially suitable for processes that require precise control of the drop size, such as emulsification and crystallization.^5,6

It has been established that rapid mixing within a droplet is based on chaotic flow.^7–10 However, it is difficult for droplets travelling along a straight microchannel to develop chaotic diffusion due to the symmetrical circulations of the two flows in each half of the droplet,¹¹ even if the internal circulations can elongate and fold the contents of each half of the droplet. Chaotic advection is often formed in a serpentine microchannel, where stretching and folding of the entire droplet takes place in winding sections.¹² Song et al.^12,13 reported that during travelling in a serpentine microchannel, the droplet remained asymmetrical between the two inner circulating flows. Constrained by the channel geometry, the part of the droplet exposed to the inner arc of the serpentine channel generated a small circulation while the other part generates a big one. The asymmetric circulation allowed mass transfer along not only the axial but also the radial direction.

Kashid et al.¹⁴ used particle image velocity (PIV) technology and a computational fluid dynamics (CFD) numerical simulation method to investigate the fluid circulation within micro-droplets. The CFD particle tracking technique was also employed to investigate the changes in the internal circulation regions and stagnation regions. They indicated that the flow rate had little effect on the internal circulation at low flow velocities when a droplet had sufficient length.

Chandorkar et al.¹⁵ added a scalar equation to the volume-of-fluid methodology (VOF) for simulating liquid–liquid Taylor flow. Through tracing concentration changes in the micro-droplets, they examined the mixing process of multiphase flow under different operating conditions. Muradoglu et al.¹⁶ used the finite-volume (FV)/front-tracking (FT) method to study the effects of various non-dimensional parameters on the mixing quality during the moving process of droplets through a serpentine channel. They reported that the capillary number, the relative size of the droplet compared to the channel width, and the viscosity ratio of the dispersed fluid to the continuous fluid were the significant factors affecting the mixing process, however, the Reynolds number had little effect on the process.

Despite these studies, a detailed mixing mechanism of droplet-flow in serpentine micro-channels has not been established. There are no data for the influences of droplet formation, droplet size, and hydrodynamics on the mixing performance in a serpentine microchannel.

In this paper, we report the intermixing behavior of the droplet flow in a 3D serpentine channel-based microfluidic system. The volume-of-fluid methodology (VOF) was used to simulate the hydrodynamics and multiphase flow pattern. A user defined scalar (UDS) was added to the continuity equation of the dispersion phases to investigate the fluid mixing in droplets. Both the fluid distribution within the droplet and the droplet size seriously affect the mixing process. The best mixing occurs when the droplet size is comparable with the channel width.

Numerical model

Governing equations

The volume of fluid (VOF) methodology was selected because it could accurately describe the generation of the droplet, the multiphase interface and the wetting behavior of the fluids.¹⁷ For the VOF simulation, we used two types of incompressible Newtonian fluids (water and oil, with constant density) to calculate the dispersed phase and the continuous phase. Phase changes and viscous dissipation were neglected in the numerical model. So the governing equations for velocity (v) and pressure (p) can be described by the following equations: continuity, eqn (1):

∇·v = 0 (1)

and momentum, eqn (2):where ρ is the fluid density, μ is the fluid viscosity and F is the volumetric surface tension. p and μ were kept constant in the simulation.

In the present work, the surface tension and wall adhesion were considered to be the main parameters and the gravitational effect was neglected because of the small influence when compared with the surface tension, viscous and inertial effects. The external force vector can be calculated using eqn (2) where σ is the surface tension coefficient, α(x, t) is the volume fraction within the cell, i.e., α(x, t) = 1 in a cell filled with oil, α(x, t) = 0 with water, and κ is the free surface curvature which can be calculated using eqn (3):

where n is the perpendicular vector to the interface between the phases.

Once the primary phase and the secondary phase were defined, the interface between the two phases could be calculated using the convection, eqn (4):

The user defined scalar (UDS) was added to the continuity equation. The mixing process was tracked with a mixture fraction variable for the kth mixture, f_k, which was governed by the general transport, eqn (5).

where D represents the diffusion coefficient which is assumed to be constant for all mixture fractions.

Boundary conditions

A horizontal serpentine microchannel was used as the simulation object, as shown in Fig. 1. To investigate the mixing process within moving droplets, droplets were introduced from a cross-type inlet with fluids set in a radial direction initially. The cross section of the main channel was square with width w = 0.05 mm. The mesh was generated using the pre-processor software GAMBIT (ANSYS Inc.). A structured mesh was generated and the 3D geometry was considered, the mesh size of 0.003 mm was adopted in the simulation. Inlets 1, 2 and 3 were defined as velocity inlets. The outlet and the developed outflow were defined, and a no-slip boundary condition was set at the walls.

Fig. 1 Schematic of the microchannel with a cross inlet (channel depth of 0.05 mm).

Inlet 1 was used for feeding the continuous phase (oil, density of 1900 kg m⁻³, viscosity of 0.0051 Pas), while inlets 2 and 3 generated the dispersed phase (water, density of 1000 kg m⁻³, viscosity of 0.001 Pas). The surface tension (σ) between the two phases was set at 0.013 N m⁻¹, and the contact angle (θ) was 135°.

To accurately simulate the nonlinear interface movement of the multiphase flow, the PLIC interface reconstruction method¹⁸ was used to adjust the curvature of the moving interface during calculations. The CSF (continuum surface force) model¹⁹ was adopted to deal with the surface tension. The QUICK format for convection was used in the simulation to reduce the impact of numerical diffusion. PISO (pressure-implicit with splitting of operators) was adopted for pressure velocity coupling, which not only reduced the internal iteration per time step but also allowed larger under-relaxation factors. The simulations were carried out on a desktop computer (Intel Xeon 8 cores).

The time step size is an important parameter in the simulation of multiphase flow. The time step size can be calculated using the Co (Courant number), defined as:

where U, Δt and Δx are the average flow velocity, the time step size and the mesh size, respectively. A Co of 0.25 was used for the simulation originally to fix the step size under different mesh sizes.

Quantitative measurement of the mixing process

The mixing degree of a droplet along the channel at each time step was calculated using:where c is the mixing fraction in different cells, c_∞ is the mixing fraction in a well mixed cell (in our case c_∞ = 0.5), and c₀ is the initial mixing fraction (c₀ = 1.0).

Quantitative measurement of the circulation period

A discrete phase model (DPM) was used to investigate the circulation period within the droplet. 10 uniformly spaced tracer particles were placed along the radial centre line of the droplet. When each particle travelled along with the flow, the average time of 10 particles to finish one period was adopted as the circulation period. A third-order MUSCL scheme was used to reduce numerical diffusion, and a four-order Runge–Kutta algorithm was also used to minimize the time integration error.

Model validation

To validate this model, a T-junction microchannel was constructed with the same geometry dimensions as that reported by Song et al.¹³ The actual fluid properties reported were chosen for the simulation. The simulated phase pattern (Fig. 2) is almost the same as the experimental results.¹³ The ratio between the droplet length of the dispersed phase and the microchannel width is 1.4167 (the average flow velocity U = 50 mm s⁻¹, and the water fraction wf = 0.52), which is comparable with the experimental data. The simulation can effectively capture the droplet formation and predict the droplet size precisely.

Fig. 2 Simulated droplet formation for validation with experimental results¹³ (U = 50 mm s⁻¹, wf = 0.52, Re_d = 2.50, and Ca_d = 0.0036).

Results and discussion

Droplet formation

The formation of droplets can be characterized using the capillary number, Ca_d and the Reynolds number, Re_d.

Ca_d = Uμ_d/σ (8)

Re_d = wUρ_d/μ_d (9)

where ρ_d and μ_d are the density and viscosity of the dispersed phase, respectively. We operated our system at low values of Ca_d (0.0036) for the formation of droplets. The corresponding Reynolds numbers is 2.5.

Droplets of different sizes can be formed by changing the water fraction, wf = U_d/(U_d + U_c), where U_d is the flow rate of the dispersed phase and U_c is the flow rate of the continuous phase. The flow pattern was calculated by changing the water fraction in the range of 0.05–0.70, meanwhile keeping the fluid velocity in the main channel (U) constant at 50 mm s⁻¹. As shown in Fig. 3, with increasing the wf, the droplet size increased significantly.

Fig. 3 Droplet formation under different water fractions (U = 50 mm s⁻¹, Re_d = 3.10, and Ca_d = 0.0036).

The droplet size can be determined from the ratio of the flow rate between the dispersed phase and the continuous phase, wf/(1 − wf). As shown in Fig. 4, the dimensionless length of the droplet L_d/w was found to show a linear relationship with wf/(1 − wf), which can be fitted using:

L_d/w = 0.898 + 0.653wf/(1 − wf) (10)

Fig. 4 Dependence of the L_d/w ratio on the flow rate ratio (U = 50 mm s⁻¹, Re_d = 3.10, and Ca_d = 0.0036).

Mixing process

The user defined scalar was provided to show the mixing state, as shown in Fig. 5. Two dispersed phases (aqueous phase) entered the cross junction from inlets 2 and 3 and formed droplets. The two dispersed phases presented a radially symmetric distribution in the droplet, and each occupied one half of the droplet. The distribution evolved with the movement of the droplet along the microchannel. For example, the concentration of the tracer from inlet 2 changed from 1 (initial stage, red) to 0.5 (well mixed state, green).

Fig. 5 Fluid mixing in droplets with different water fractions (U = 50 mm s⁻¹, Re_d = 3.10, and Ca_d = 0.0036).

When the droplet travelled through the initial straight section of the channel, the fluids remained unmixed. Once the droplet entered the winding section, it became deformed significantly, and the fluids within the droplet were stretched and folded, which facilitated the diffusion process. Besides, a reorientation of the fluids occurred within the droplet, thus part of the fluids were distributed in the axial direction, which benefited mixing in the next straight section of the channel. With the increase of the water fraction, the mixing distance for good mixing increased significantly. The fluid mixing processes in droplets with different water fractions are given in the ESI† Movies for wf = 0.1, 0.3 and 0.7.

The mixing degree at different locations (from P1 to P15) through the microchannel was calculated and are shown in Fig. 6. It is obviously shown that regardless of the water fraction, the mixing degree within the droplet always reached a high level in the serpentine channel. At location P6, i.e. a droplet travelling distance of about 0.8 mm, the mixing degree reached a level of between 70% and 90%. However, the difference in droplet mixing with these water fractions was still obvious. Better mixing occurred when the water fraction was low, and increasing the water fraction required a longer travelling distance for the droplet to reach a mixing degree of 90%. Shorter droplets will be mixed in a shorter mixing distance.

Fig. 6 Effect of the droplet travelling distance in the microchannel on the mixing degree (U = 50 mm s⁻¹, Re_d = 3.10, and Ca_d = 0.0036).

Hydrodynamics within droplets

To gain a detailed understanding of the mixing mechanism, the hydrodynamics inside the droplet were studied. When a droplet moves through the channel, circulation flows are generated within the droplet due to the shear stress near the wall, and this contributes to the mixing of the liquid within the droplet.

To obtain the relative velocity field inside the droplet, the superficial velocity of the droplet (50 mm s⁻¹) was subtracted from the velocity field. The relative velocity field in the straight channel and the streamline of the circulations are shown in Fig. 7. The relative velocity field of the droplet in the straight section showed three pairs of symmetrical circulation patterns. The symmetric main circulation can be observed in the main part of the droplet, and two smaller circulations sat in the front and the back parts of the droplet. Here we refer to the main circulations as “major circulations” and the two side circulations as “minor circulations”. Each pair of the circulations remain symmetrical in the radial direction along the straight channel. Thus the two sides of the droplet will be mixed individually.

Fig. 7 Velocity field and streamline of the main circulations in a 2D plane in the straight section of the microchannel, (U = 50 mm s⁻¹, Re_d = 3.10, and Ca_d = 0.0036).

However, when the droplet travelled through the serpentine section, asymmetric circulations took place within the droplet. As shown in Fig. 8, the relative velocity field in the winding section shows three asymmetric circulations. The asymmetric circulations between the inner and outer sides of the arc can be seen. This led to fluid reorientation within the droplet. This result is consistent with those reported by other researchers.²⁰

Fig. 8 Velocity field in a 2D plane in the winding section of the microchannel (U = 50 mm s⁻¹, Re_d = 3.10, and Ca_d = 0.0036).

In addition, particles were added to investigate the circulation period in the center line of the droplet. As shown in Fig. 9, the circulation period of the major circulation has a linear relationship with the dimensionless droplet length, L_d/w, where L_d is the droplet length of the dispersed phase. It suggests that the time needed for effective circulation mainly depends on the droplet size. The fluid circulation is faster within smaller droplets.

Fig. 9 Relationship between the per-circulation time and the relative droplet length (L_d/w) (U = 50 mm s⁻¹, Re_d = 3.10, and Ca_d = 0.0036).

Mixing mechanism

Based on the mixing process and the relative velocity field within the droplet, the mixing of fluids in a serpentine microchannel has two different circulation states, (1) symmetric circulations in the straight section, and (2) asymmetric circulations in the winding section. Fig. 10 schematically illustrates mixing process while the droplets travel through the channel.

Fig. 10 Schematic diagram of the mixing process in a droplet.

In the straight section, symmetric major circulations within the droplets are the driving force of fluid mixing. The initial spatial distribution of the mixing components within the droplet is very important for fluid mixing in the droplet. Fluids remain unmixed when distributed in a radial direction (left and right). Meanwhile, the droplets with an axial distribution of the reagents (front and back) could achieve mixing easily.

When the droplet travelled in the winding section, the inner fluid reoriented, and the major and minor circulations caused the fluids to rotate in the opposite direction. Thus the inner fluid was distributed in the axial direction, which benefited the mixing process in the next straight section.

Besides, the droplet interface was stretched and folded during the periodic alternation of the straight and winding sections, thus providing chaotic mixing.

To understand the effect of the droplet size on the mixing degree, the relationship between the mixing distance and the droplet length was examined. Here, the point at which the mixing degree reached 90% was defined as complete mixing and the travelling distance required for the droplet to reach 90% mixing was referred to as the mixing distance. The mixing distance (L) and droplet length (L_d) were divided separately by the channel width (w = 0.05 mm) to become dimensionless. The result is shown in Fig. 11.

Fig. 11 Relationship between the relative mixing distance (L/w) and the relative droplet length (L_d/w) (U = 50 mm s⁻¹, Re_d = 3.10, and Ca_d = 0.0036).

With the increase in droplet size, a longer mixing distance is required for complete mixing. The mixing performance is controlled by two factors, i.e. the spatial distribution of the mixture and the circulation period within the droplet.

The best mixing was found to occur when the droplet size was comparable with the channel width, which was in accordance with the results reported by Bringer et al.²¹ and Tung et al.²² For smaller droplets (L_d/w < 1), a high circulation speed and a marked axial redistribution caused the best mixing performance. However, a much smaller droplet is unsteady and the flow pattern would transfer from slug flow to bubbly flow.

With the increase in droplet size (1 < L_d/w < 2.0), the mixing distance increased slowly. Although the asymmetric circulations caused the fluids to easily distribute in the axial direction, increasing the droplet length led to reduction in the circulation period. But when the droplet length was larger than 2.0, the asymmetric circulations had little effect on the reorientation, which enabled the radial distribution of the fluids within droplet. In addition, longer a circulation period also led to lower mixing efficiency.

Conclusions

The mixing of fluid in droplets follows two routes: symmetrical circulations in the straight section of the channel and asymmetrical circulations in the winding section. The asymmetrical circulations lead to the occurrence of fluid reorientation, which benefits the mixing process in the next straight section.

Our study has indicated that the droplet length shows an evident effect on the mixing degree. When the dimensionless droplet size is less than 2.0, asymmetric circulations cause the fluids to easily distribute in the axial direction, which leads to fast mixing. When the dimensionless droplet size is larger than 2.0, the reorientation generated by the asymmetric circulations becomes insignificant, thus the fluids within the droplet mainly distribute in the radial direction, observably increasing the mixing distance.

The best mixing performance occurred when the droplet size was comparable with the channel width. Under a slug multiphase flow regime, shorter droplets generated with a smaller water fraction are preferred for achieving fast mixing, but it is easy to transfer to a bubbly flow pattern and is not an economic method of operation due to the use of more of the continuous phase. These findings may be useful for designing efficient microfluidics for applications in the fields of emulsification, biochemical analysis and rapid chemical synthesis.

Acknowledgements

This work was supported financially by the National Natural Science Foundation of China (21276222), the National Basic Research Program of China (2011CB606001) and the State Key Laboratory of Chemical Engineering (SKL-ChE-13D01).

Notes and references

1. W. Ehrfeld, V. Hessel and V. Haverkamp, Microreactors, Wiley Online Library, 2000.
2. K. Jahnisch, V. Hessel, H. Lowe and M. Baerns, Angew. Chem., Int. Ed. Engl., 2004, 4, 406–446.
3. B. J. Burke and F. E. Regnier, Anal. Chem., 2003, 8, 1786–1791.
4. F. Sarrazin, K. Loubiere, L. Prat, C. Gourdon, T. Bonometti and J. Magnaudet, AIChE J., 2006, 12, 4061–4070.
5. A. Huebner, S. Sharma, M. Srisa-Art, F. Hollfelder, J. B. Edel and A. J. Demello, Lab Chip, 2008, 8, 1244–1254.
6. H. Song, D. L. Chen and R. F. Ismagilov, Angew. Chem., Int. Ed. Engl., 2006, 44, 7336–7356.
7. H. A. Stone, A. D. Stroock and A. Ajdari, Annu. Rev. Fluid Mech., 2004, 36, 381–411.
8. R. H. Liu, M. A. Stremler, K. V. Sharp, M. G. Olsen, J. G. Santiago, R. J. Adrian, H. Aref and D. J. Beebe, J. Microelectromech. Syst., 2000, 2, 190–197.
9. A. D. Stroock, S. K. W. Dertinger, A. Ajdari, I. Mezic, H. A. Stone and G. M. Whitesides, Science, 2002, 5555, 647–651.
10. K. Bajer and H. K. Moffatt, J. Fluid Mech., 1990, 212, 337–363.
11. K. Handique and M. A. Burns, J. Micromech. Microeng., 2001, 5, 548–554.
12. H. Song, J. D. Tice and R. F. Ismagilov, Angew. Chem., Int. Ed., 2003, 7, 768–772.
13. H. Song, M. R. Bringer, J. D. Tice, C. J. Gerdts and R. F. Ismagilov, Appl. Phys. Lett., 2003, 22, 4664–4666.
14. M. N. Kashid, I. Gerlach, S. Goetz, J. Franzke, J. F. Acker, F. Platte, D. W. Agar and S. Turek, Ind. Eng. Chem. Res., 2005, 14, 5003–5010.
15. A. Chandorkar and S. Palit, Sens. Transducers J., 2009, 7, 136–149.
16. M. Muradoglu and H. A. Stone, Phys. Fluids, 2005, 7, 073305.
17. M. van Sint Annaland, W. Dijkhuizen, N. Deen and J. Kuipers, AIChE J., 2006, 1, 99–110.
18. D. Youngs, Time-dependent multi-material flow with large fluid distortion[C], Numerical Methods for Fluid Dynamics, ed. K. W. Morton and M. J. Baines, Academic, New York, 1982, pp. 273–285.
19. J. U. Brackbill, D. B. Kothe and C. Zemach, J. Comput. Phys., 1992, 2, 335–354.
20. Z. Che, T. N. Wong and N. T. Nguyen, Int. J. Heat Mass Transfer, 2010, 9, 1977–1985.
21. M. R. Bringer, C. J. Gerdts, H. Song, J. D. Tice and R. F. Ismagilov, Philos. Trans. R. Soc. London, A, 2004, 1818, 1087–1104.
22. K. Y. Tung, C. C. Li and J. T. Yang, Microfluid. Nanofluid., 2009, 4, 545–557.

---

Footnote

† Electronic supplementary information (ESI) available: See DOI: 10.1039/c5ra21181f

---