Computational investigations of the mixing performance inside liquid slugs generated by a microfluidic T-junction

Since its emergence in the 1990s, microfluidics has received extensive research interests due to its unique characteristics. Compared to their traditional counterparts in lab-scale, microfluidic devices offer a number of advantages such as reduced sample and reagent consumptions, short analysis time, high throughput, and improved heat and mass transfer efficiencies;^1,2 therefore, they have been adopted by the biology and chemistry societies with numerous applications.^3–5 However, conventional single-phase microfluidics confronts several challenges. Since the flow in microfluidic systems is usually operated in the laminar flow regime, mixing between the fluid streams moving in the same direction is governed by diffusion across the fluid interface.⁶ In order to enhance the mixing process, passive microfluidic mixers⁷ or active mixing techniques such as electro-kinetic flow⁸ and surface acoustic waves⁹ are required by many applications to bring in convective mixing. Due to the parabolic velocity profile of the laminar flow, dispersion of the reagents occurs in conventional microfluidic devices and leads to cross-contamination and dilutions of reagent samples.⁴ In addition, the direct contact of reagents with channel walls may cause surface deposition, which is especially critical in the applications of biomaterial and nanomaterial syntheses. For example, Wagner and Kohler have reported that the synthesized nanoparticles deposited preferentially on the channel wall and led to reactor fouling.¹⁰ If the surface deposition is not controlled properly, the channels may be clogged eventually.

The droplet-based microfluidics has been developed to address the above mentioned challenges and to improve the performance of microfluidic devices for many applications.^11–13 This technique introduces an immiscible carrier fluid such as a gas or oil stream to encapsulate aqueous reagents inside discrete slugs/droplets with nanoliter sizes, and the slugs or droplets can be used for different reactions and analyses.^4,14 Compared to the conventional single-phase microfluidics, the droplet-based microfluidics offers several unique advantages.^4,14,15 For example, mixing inside small slugs is very efficient because of the secondary flow inside slugs. The shear force between the rigid channel wall and moving slugs generates internal recirculating flow inside both phases.¹⁶ The recirculating flow accelerates mixing in the transverse direction thus results in fairly uniform reagent concentrations inside slugs in a short time.^11,12 The encapsulation of the reagent solutions inside slugs can eliminate the transverse dispersion effectively.¹³ Between the channel wall and the travelling slugs, a thin film of the carrier fluid is developed that prevents the direct contact of reagent with channel wall,^13,17 thus the surface deposition can be avoided as well. In addition, since slugs are travelling with same speed, they have identical residence time when they leave the microfluidic channel; in other words, the reagents inside each slug experience identical reaction time. This characteristic enables the production of uniform particles in the applications of nanomaterial and biomaterial syntheses.^1,5 Because of such advantages, the droplet-based microfluidics has been utilized extensively in material synthesis, reaction kinetics measurement, protein crystallizations, and many other applications.

However, several challenges still remain in applying droplet-based microfluidics. First is to control the flow pattern properly in slug/droplet flow regime. The liquid–liquid two phase flow exhibits several types of flow patterns in microfluidic devices. The flow may appear as droplets/slug flow, jet flow, parallel flow, or annual flow based on operating conditions and the geometry of the devices.¹⁸ For example, slug flow shifts to jet flow or parallel flow in a microfluidic T-junction if the operating flow rate is increased above certain critical values.^19,20 The dispersed phase does not break into slugs in the junction region but forms a long segment; slugs with large size distributions are separated from the segment downstream the junction due to capillary instabilities.^19,21,22 In some circumstances, the segment may be very long.²² The jet flow is not desired for droplet-based microfluidics because the breakup process is unstable. The second challenge is to achieve rapid mixing inside slugs, which is the key issue in the applications of in situ kinetic studies and material syntheses. For example, Nichols et al. have recently utilized a droplet-based microfluidic system to measure the rate constant for interfacial mass transfer in the TALSPEAK (Trivalent Actinide–Lanthanide Separation by Phosphorous reagent Extraction from Aqueous Komplexes) process. They pointed out that the droplets must have sufficiently large mass transfer rates in order to conduct the measurement successfully.²³ In the applications of material syntheses, uniform reagent concentrations inside slugs are critical to obtain particles with minimum size distributions. For example, Zhao et al. have carried out the preparation of curcumin nanoparticles in a microfluidic reactor involving gas–liquid flow, and they reported that efficient mixing inside liquid slugs has led to nanoparticles with narrow size distributions.²⁴ Li et al. have also reported that the size distributions of the copper nanoparticles synthesized by liquid–liquid slug flow are influenced by the mixing performance inside the slugs.²⁵ 

Since 1961, Taylor has started the investigations of the internal flow inside slugs.²⁶ In the past decade, a number of experimental investigations have been reported to literatures which focus on the fluid dynamics of the droplet/slug flow in microfluidics devices. For example, Gunther et al. have exhibited the velocity field in a liquid slug between two adjacent gas bubbles by using confocal scanning microscopy and image processing.²⁷ By utilizing the high-speed confocal micro-particle image velocimetry (μ-PIV) technique, Kinoshita et al. have measured the three-dimensional velocity field of a moving slug travelling inside a microfluidic channel.¹⁶ Both these studies have revealed the recirculating flow patterns inside liquid slugs. From their experimental study, Tice et al. have observed the so-called “twirling” effect which occurred during the slug formation stage. The twirling effect re-orientates the flow direction of the dispersed phase and distributes the reagents in the transverse direction. Based on their observations, Tice et al. have concluded that the twirling effect played a critical role in the mixing performance inside slugs. For example, if the reagents were initially located on the path of the recirculating flow, the mixing was governed by convection resulted from the recirculating flow thus was very efficient; in contrast, if the reagents were initially distributed off the flow path, the mixing depended on the slow diffusion process.¹² Recently, Han et al. have reported their study about tracking the mixing process inside a single droplet moving along a microfluidic channel. Their results also demonstrated the dependence of mixing efficiency on the initial distributions of the reagents.²⁸ Nowadays, numerical simulations have provided an alternative way to understand the fundamental fluid dynamics. While the majority of the numerical studies have been carried out to understand the hydrodynamics of the slug formation process,^29,30 there are limited reports discussing about the mixing process inside slugs. Previous numerical studies carried out by Kashid et al.³¹ and Rhee and Burns³² have adopted the moving frame method with pre-defined concentration field to study the mixing process within one slug or between the adjacent slugs. Although their investigations have revealed the insights of the mixing process in some aspects, several key issues are still not answered. For example, the actual initial distribution of the reagents controlled by the twirling effect has not been considered in these investigations.

Since rapid mixing inside droplets/slugs is critical for many applications, the investigations on this topic are still continuing.^24,28 In this study, we have adopted numerical simulations to gain deep insights of the dynamic mixing process inside slugs. The aim of this study is to understand how the normal operating conditions such as operating velocity and flow ratios influence slug sizes as well as the mixing performance. The initial distribution which is controlled by the twirling effect is included in the numerical simulations by modeling the slug formation and mixing process simultaneously. The accuracy of the numerical model has been validated by comparing the predicted slug lengths and periods from simulations with the experimental data reported by Tice.¹² By adopting the mixing efficiency index, the mixing efficiency inside slugs is quantified with respect to residence time at various operating conditions. The simulation results may enhance the understanding of the mixing process inside slugs and provide guidance to the design and operation of droplet-based microfluidic systems.

Several numerical methods have been applied by researchers to study the slug formation process, i.e., volume of fluid (VOF) method,³³ level-set method,³⁴ lattice-Boltzmann method,³⁵ and phase field method.³⁶ Among them, VOF method has been applied extensively in modeling multiphase flow in microfluidic channels since it is efficient in tracking the topological changes of the fluid interface by using a relatively simple algorithm. This method also has other advantages such as reasonable accuracy, relative simplicity, and ability to solve complex free surface flow with minimum consumption of computational resource. In our previous work of studying the slug formation dynamics inside a millifluidic reactor, the slug lengths and formation frequency predicted by VOF method achieved good agreement with the experimental observations.²⁵ Therefore, VOF method was selected in this work to predict the slug formation inside a microfluidic T-junction. Coupled with VOF method, a species transport model was applied to the dispersed phase to study the mixing process of two solutions.

The VOF method relies on the fact that different phases are not interpenetrating, and the volume fractions of these phases are addable in computational cells. It solves only one set of Navier–Stokes equations for all the phases thus saves the computational cost significantly. In this study, the continuous and the dispersed phases are denoted as p and q, respectively. Each cell in the computational domain may contain one or two phases. A new variable α_q (dimensionless) is introduced to identify the underlying fluid volume of qth phase inside each computational cell. If α_q = 1, it implies that this cell is full of qth phase; in contrast, α_q = 0 represents that the cell is empty of qth phase. If α_q is within the range from 0 to 1, it implies that this cell contains the fluid interface between these two phases.

The governing equations—continuity (1) and momentum equations (2) are solved throughout the computational domain

in which $v→$ is the velocity vector (unit: m/s); ρ (unit: kg/m³) and μ (unit: Pa s) are the averaged density and viscosity based on the volume fraction of each phase. They are calculated by Eqs. (3) and (4), respectively,

The volumetric body force term $F→$ (unit: kg/(m² s²)) in Eq. (2) contains only the interfacial tension, which is computed by the continuum surface force (CSF) model³⁷ as shown in the following equation:

where σ (unit: N/m) is the surface tension coefficient; $n→$ (unit: 1/m) is interface normal vector; $n̂$ (dimensionless) is the unit interface normal vector, and κ (unit: 1/m) is the local surface curvature. They are estimated by

The dimensionless Bond number (defined as $Bo=(ρp−ρq)gL2/σ$⁠, of which g = 9.81 m/s² and the characteristic slug length L = 100 μm) indicates the relative importance of the buoyancy force to the interfacial tension force. A typical value of Bo in this study is 0.0063. Such a small number suggests that the buoyancy force is negligible compared to the interfacial tension force.³⁸ Therefore, the buoyancy force was not included in the volumetric body force term $F→$⁠.

The VOF method uses the geometric reconstruction scheme which adopts a piecewise linear interpolation³⁹ to obtain the reconstructed interface. The interface is tracked by solving an additional advection equation (9) to yield the volume fraction of each fluid

When the interface contacts the channel wall, the wall adhesion force exists in the boundary cells and affects the local interface curvature. The wall adhesion force is taken into account by specifying a contact angle θ_w, and the unit interface normal in the cells near the wall is calculated by

in which $n̂w$ and $t̂w$ are the unit vectors in the normal and tangential directions to the wall. The CSF model combines the contact angle with the normal direction of the surface calculated by one cell away from the wall to determine the local curvature κ. This curvature is then applied to Eq. (5) to calculate the body force term near wall. In this study, a constant contact angle is specified to all the wall boundaries, and its value is assumed to be identical at all flow rates. The effect of the contact angle values on slug sizes is discussed in Appendix  B.

In their experiments, Tice et al. adopted a perfluorodecalin (PFD) solution as the continuous phase. The dispersed phase contained two solutions: one was an aqueous dye solution; the other was a transparent solution. As shown in Fig. 1, the PFD solution was injected from the straight inlet while the dye and water solutions were injected from the side inlet. In the species transport method, a species transport equation, Eq. (11), is used to model the mixing process

in which $Ci$ is defined as the concentration of the dye solution, $v→$ is the velocity field solved by Eqs. (1) and (2), and D is the diffusion coefficient of dye in water. The dye solution is assumed to have a constant diffusivity as $10−11$ m²/s. In the side inlet, the dye solution has the concentration value of 1, while the water solution has 0. Since the solubility of the dye and water solutions in the PFD solution is negligible, zero flux condition is applied to the fluid–fluid interface while solving Eq. (11); in other words, the dye and water solutions cannot transport to the continuous phase, and the mixing process occurs only inside the dispersed phase.

The mixing performance is quantified by calculating the mixing efficiency index η (shortened as mixing index in the following text), which evaluates the deviation from the perfectly mixed state through

In this equation, N is the total number of cells inside slugs. The variable $C¯s$ represents the scaled concentration value in each cell. The variable $C¯s0$ is the scaled concentration values if the solutions are completely unmixed. One should notice that $C¯s0$ is either 0 or 1 depending on the initial fluid in this cell. If the cell is initially occupied by water, $C¯s0$ is 0; otherwise, it is 1 as the cell is initially filled with dye. The variable $C¯s*$ is the scaled concentration values if the solutions are perfectly mixed. When these variables are used to estimate Eq. (12), η is in the range from 0 to 1.

During the simulations, a user defined function was used to extract the position of each cell, the volume fraction of the dispersed phase (⁠$αq$⁠), and the concentration values (⁠$C¯s$⁠) out of the computational domain. The perfect mixing condition is then calculated by

in which V is the volume of each computational cell and $αq$ is the volume fraction of the dispersed phase.

In this work, the commercial CFD package ANSYS Fluent 13 was adopted to solve the governing equations shown in Secs. II A and II B. The PISO (pressure-implicit with the splitting of operators) scheme was selected for the pressure–velocity coupling of the momentum equations and the continuity equation. For the spatial discretization, Green–Gauss node based method was used for gradients; the PRESTO! (pressure staggering options) scheme was used for pressure interpolation; the second order upwind difference scheme was used for the momentum equations.

The three-dimensional (3D) geometry and the computational grids of the microfluidic T-junction were generated by ANSYS ICEM. The dimensions of the channels followed the experimental setup reported by Tice et al. in 2003, of which the width and depth of the channels were both 50 μm.^11,12 As shown in Fig. 1, the length of the straight channel after the junction region was set as 1000 μm. This length was verified to be sufficient to hold at least two slugs. The inlet channels were truncated in the simulations so as to reduce the computational cost; however, proper values were assigned to ensure that laminar flows had been fully developed prior to the junction. The two inlets were specified as the velocity inlet boundaries while the outlet was specified as the outflow boundary. The no-slip boundary condition was imposed on the channel walls associated with a constant contact angle. Proper values were assigned to the time step in order to ensure that the maximum Courant number (Cr) during the simulations was less than 0.25. The simulations were performed with 8 processors on the cluster Philip (the High Performance Computing facility located at Louisiana State University).

The fluidic properties used in the simulations were taken from the paper published by Tice et al. in 2003.¹² In their work, the continuous phase was made by a 10:1 mixture of PFD and 1H, 1H, 2H, 2H-perfluoroocatanol. The dye and water solutions used by them had properties very close to pure water. They have measured the interfacial tension through experiments of which the value was reported to remain constant for all their explored cases.¹² The values of the density and viscosity of the two phases as well as the interfacial tension are listed in Table I.

The accuracy of the VOF method on predicting the slug formation dynamics was validated by comparing the predicted slug lengths and formation frequencies from various operating conditions to the experimental data reported by Tice et al.¹² In their experiments, they adopted a length variable “period” to define the slug formation frequency, which avoided measuring the slug formation time directly. The period is defined as the center to center distance between two adjacent slugs. The operating conditions were controlled by two parameters: water fraction (wf, defined as the ratio of the volumetric flow rate of the dispersed phase to the total flow rate) and total superficial velocity (U_total). For a particular droplet-based microfluidic system, wf and U_total are the primary controlling parameters over slug formation processes; therefore, the effects of these two parameters on the slug lengths and periods are discussed in Sec. III B. After the VOF method has been validated, the mass transfer model is applied to the dispersed phase. The dynamic mixing process inside slugs is quantified by the mixing indices η with respect to the residence time inside the T-junction. The discussions about the mixing efficiencies at various operating conditions are provided in Sec. III D.

In order to eliminate the potential errors resulted from the inadequate grid resolutions, a mesh dependence study was carried out with three types of structured grids. The edge lengths of these grids were specified as 5 μm, 2 μm, and 1 μm, respectively, corresponding to the coarse, medium, and fine grid resolutions. These grids were tested with the operating conditions of U_total = 50 mm/s and wf = 0.40. The averaged slug lengths and periods predicted from the simulations by using these grids are shown in Table II. As the variations of the slug lengths and periods obtained from the medium and fine grid resolutions are within 5% of each other, the numerical model is deemed to reach mesh independence by using medium and fine grid resolutions. Thus, the medium grid resolution was adopted for the subsequent investigations. The entrance lengths of the straight and side inlet channels were verified to be sufficient to ensure fully developed laminar flow as well. The detailed discussions related to the mesh dependence study for large flow rates and the effect of entrance lengths are provided in Appendix  A.

The fluid flow in microfluidic channels is usually operated with very small Reynolds numbers. (Re is defined as Re = ρU_totalD/μ, in which ρ and μ are the volume-averaged density and viscosity of the continuous phase; D is hydraulic diameter of the microfluidic channel.) In this study, Re ranges approximately between 0.39 and 2.93. The small values of Re indicate that the inertial effect is negligible compared to the viscous effect. Therefore, the slug formation process is governed by three effects primarily. Namely, they are the interfacial tension, the viscous shear stress, and the pressure gradient across the fluid interface.^36,40 The interactions of these three effects can lead to three distinct breakup regimes, i.e., squeezing, dripping, and jetting regimes, inside a microfluidic T-junction. Moreover, the breakup process influences the mixing performance inside the slugs as well (quantitatively presented in Sec. III D). In this section, a typical slug breakup process inside the microfluidic T-junction is demonstrated with the help of simulations.

As seen in Fig. 2, the slug breakup processes in the squeezing and dripping regimes can be summarized as four steps:

• The tip of the dispersed phase penetrates into the continuous phase, and the incipient slug starts to grow in the junction region.

• As the dispersed phase accumulates in the junction region, it blocks the flowing channel of the continuous phase. To maintain its constant flow rate, the continuous phase accelerates around the tip of the dispersed phase. A high pressure is built up in the upstream of the continuous phase so as to drive the flow. Meanwhile, the viscous shear stress is elevated by the increased velocity as well.

• Because of the shear stress and the high pressure in the continuous phase, the fluid interface is pushed toward the upper corner of the junction, resulting in a visible neck in the junction region. If the shear force is sufficiently strong, the neck may break up before the junction is fully blocked by the dispersed phase. This phenomenon usually occurs in the dripping regime. If the shear force is relatively weak, the breakup process is dominated by the pressure gradient across the interface. The dispersed phase has to block the entire flowing area of the continuous phase so as to accumulate sufficiently high pressure. This breakup pattern is usually referred as in the squeezing regime.

• When the neck is squeezed to a critical width, it is snapped instantly by the interfacial tension force. A slug with a well-defined shape is detached from the dispersed phase. The tip of the dispersed phase retracts to the junction region, and a new cycle starts.

In droplet-based microfluidic devices, the slug breakup process is usually controlled by manipulating the flow rates of the continuous phase and the dispersed phase. In their experiments, Tice et al. studied the relationship of slug sizes versus wf when U_total was constant. They found that short slugs with a bullet shape were obtained when the wf was small; if the wf was large, the resultant slugs exhibited long and cylindrical shapes.¹² Similar phenomena were also observed from slug flow inside microfluidic channels using silicone–water⁴⁰ and kerosene–water⁴¹ solutions. All these studies suggest that slug sizes can be tuned easily through manipulating the flow rates of the two phases.

In order to validate the VOF model, simulations were carried out to reproduce the slug sizes and periods that reported by Tice.¹² The wf values were varied from 0.14 to 0.84, while U_total was fixed as 50 mm/s. A constant contact angle, θ_w = 140°, was specified to the wall boundaries. Although the wetting properties of the fluids with respect to the channel wall determine the flow patterns crucially,⁴² the analysis shown in Appendix  B indicates that the slug lengths and periods are marginally influenced by the value of contact angle in slug flow regime except for those cases with extremely large wf values. For the case with wf = 0.84, the simulation with the contact angle of 140° exaggerates the wall adhesion force and consequently generates extremely long slugs. After different values of the contact angle are tested by simulations, the slug length predicted by the contact angle of 150° provides the best match to the experimental results.

The obtained slug lengths and periods from simulations with respect to various wf values are plotted in Fig. 3(a), which agree well with the reported results by Tice.¹² Clearly, the relationship of slug lengths versus wf demonstrates two sections. In the section, where wf ranges from 0.14 to 0.4, the slug sizes increase slowly with wf. As shown in Fig. 3(b), the slugs break off from the dispersed phase almost immediately when their tips contact the lower side channel wall. This behavior indicates that the breakup process is governed by the dripping mechanism. Based on their numerical investigations on the slug breakup inside a microfluidic T-junction, De Menech et al. have reported that the breakup mechanism transits to dripping when the capillary number exceeds 0.01 (Ref. 36) (Ca = μ_cU_total (1 − wf)/σ, of which μ_c is the viscosity of the continuous phase). The Ca numbers of those cases with wf = 0.14–0.40 in our study conform to their findings.

In the other section with wf > 0.60, the slug sizes increase quickly with wf. As shown in Fig. 4(b), the dispersed phase forms a segment in the downstream of the junction before the slugs are separated from the dispersed phase, and the segment length increases with wf. The calculated Ca numbers indicate that the breakup mechanism in this section is governed by the squeezing mechanism. In the squeezing regime, the pressure difference between the two phases dominates the breakup process, and the dispersed phase has to block the entire cross-section of the straight channel so as to build up a sufficiently high pressure in the continuous phase. If wf increases while U_total remains constant, the viscous shear stress is weakened by the reduced flow rate of the continuous phase; therefore, the blocking time of the dispersed phase is elongated so as to build up the necessary pressure in the upstream of the continuous phase. Consequently, long slugs are produced in the squeezing regime.

In their experiments, Tice et al. have studied the slug breakup process with various total velocities (U_total = 20, 50, 100, and 150 mm/s), and they noticed that the obtained slug lengths and periods were almost identical for those explored U_total values when wf was fixed.¹² Following their experiments, numerical simulations were carried out at these U_total values while the values of wf and θ_w remained constant (wf = 0.52 and θ_w = 140°). The results are shown in Fig. 4.

As seen in Fig. 4(a), the predicted slug lengths and periods agree well with the experimental results when U_total is relatively low (U_total = 20 and 50 mm/s). The breakup processes shown in Fig. 4(b) indicate that these low velocities lead to stable slug flow. Slugs are separated from the dispersed phase in the junction region, which were observed in the experiments as well. However, if U_total is increased to 100 mm/s, the predicted breakup process by simulations deviates from the experimental observation. Compared to cases with low U_total, the slug lengths and periods obtained from this velocity are reduced noticeably. As shown in Fig. 5(b), the breaking point of the slugs propagates downstream of the junction, forming an aqueous jet from the junction to the breakup point. The simulation result indicates that the breakup process is governed by jetting mechanism, but the experimental one demonstrates that it is still in the dripping regime. If U_total is increased to 150 mm/s, the dispersed phase forms a longer jet than that of U_total = 100 mm/s, and the breakup point moves further downstream. But the corresponding breakup process in the experiments reported by Tice et al. was still in the dripping regime.¹² 

One possible reason for the deviation may be the failure of modeling of the wall adhesion force. As discussed in Appendix  B, constant contact angles ranging from 140° to 180° have been specified to the wall boundaries, and the resultant breakup modes are still governed by the jetting mechanism. The discrepancy between the simulation and experiments suggests that the numerical model using a constant contact angle is not appropriate to predict those breakup processes with high U_total values. According to the dynamic contact angle model proposed by Kalliadasis,⁴³ the value of the contact angle is a function of Ca. For the cases with high U_total, the values of the contact angle vary along the fluid interface with respect to the velocity magnitudes. To model such cases correctly, the VOF method should be coupled with the dynamic contact angle model, which specifies the value of the contact angle as a function of local velocity magnitudes.

In a short conclusion, the VOF model can predict the slug formation dynamics at low flow velocities (U_total = 20 mm/s and 50 mm/s) with good accuracy. It fails to capture the flow pattern at higher velocities (U_total = 100 mm/s and 150 mm/s). Since the focus of this study stays on mixing with low flow velocities, the dynamic contact angle is not adopted by this study. The following investigations are performed with constant contact angles at U_total = 20 and 50 mm/s.

When slugs are travelling along a microfluidic channel, their fluid elements close to the channel walls are moving slower than those in the center of the channel. Since the fluid inside slugs is immiscible with the one in the continuous phase, the fast moving elements cannot penetrate the interface thus are forced to move toward the channel wall. Therefore, secondary flow in the form of a three-dimensional toroidal vortex is formed inside slugs. As shown in Figs. 5(a) and 5(b), such vortices exist in both short and long slugs, and they exhibit two symmetric parts: upper and lower parts in the middle x–y plane. Similarly, toroidal vortices are developed inside the adjacent plugs of the continuous phase as well. As these two vortices in the continuous phase are rotating opposite to the one inside the slug, two small vortices are generated in the front and rear ends of the slug, which are indicated by the dash lines in Fig. 5(a). Similar vortex patterns have also been visualized through experiments^16,44 as well as numerical simulations⁴⁵ by other researchers. If the dye and water solutions are initially distributed in the front and back parts of the slug, the recirculating flow resulted from these vortices can mix the solutions efficiently. In contrast, if the solutions are initially distributed in the upper and lower parts, the mixing process is governed by diffusion between these two parts. The recirculating flow has no effect on accelerate mixing because the solutions are distributed off the flow path. Therefore, the twirling effect plays a critical role in mixing as it controls the initial distributions of the solutions during the slug formation stage.

Following the experimental conditions reported by Tice et al., simulations using species transport model coupled with VOF method were carried out at various wf values while U_total remained as 50 mm/s. As shown in Fig. 6(b), the dye concentration profiles predicted by simulations agree well with the experimental snapshots shown in Fig. 6(a). While the dispersed phase is growing in the junction region, the continuous phase exerts a viscous shear force on the fluid interface. The combined effect of the shear effect and the re-orientation of the flow induces a vortex in the dispersed phase. As a result, the dye solution is brought to the front tip of the slug by this vortex. This phenomenon was referred as “twirling” by Tice et al.¹² The twirling effect not only involves in the initial mixing of the dye and water solutions but also distributes these solutions transversely to the front and back parts of the slugs; therefore, after the slugs are separated from the dispersed phase, the recirculating flow can mix the solutions efficiently. When the dispersed phase contacts the channel wall, the flow path of the continuous phase is completely blocked, and the twirling effect is terminated. By comparing the dye concentration contours in Fig. 6(b), one may notice that the twirling effect is weakened as wf increases. When wf is small, especially for the cases of wf = 0.14, 0.2, and 0.3, the twirling effect is so strong that the dye–water interface in the dispersed phase is distorted toward the dye solution; after the slug is separated, the dye solution is mainly distributed on the front and rear parts of the slug. As the dye and water solutions are located on the path of the recirculating flows, the vortices inside the slug can mix them rapidly. When wf exceeds 0.4, the twirling effect is weakened, and the dye–water interface remains straight at the junction region. After the resultant slugs are separated from the dispersed phase, the majority of the dye solution locates in the lower part. The mixing process is slow since it depends mainly on the diffusion between the upper and lower parts of the slug.

In order to quantify the dynamic mixing process, the mixing indices of slugs are sampled with respect to the residence time. The sampling procedure starts at the moment when the slug is separated from the dispersed phase. It stops as the slug approaches the outlet of the T-junction. The obtained mixing indices are plotted as a function of the residence time, which yields the mixing efficiency curves shown in Fig. 6(c). As mentioned in the previous paragraph, the initial mixing at t = 0 is determined by the strength of the twirling effect. By comparing the initial mixing indices at t = 0 in Fig. 6(c), one may notice that their values decrease with wf, which provides another evidence that the twirling effect is weakened by using large wf. The mixing efficiency curves shown in Fig. 6(c) exhibit two regimes during the mixing process. They are referred as convection and diffusion dominated regimes, respectively. For the cases of wf = 0.20–0.40, the mixing process is dominated by convection initially. The mixing indices of the slugs increase quickly with respect to residence time because of the combined effects of the strong twirling effect and recirculating flow. After the dye and water concentration in the upper and lower parts become uniform individually, the mixing process is then governed by diffusion between these two parts; therefore, the mixing index increases slowly in this regime. For those cases having large wf (wf > 0.4), the entire mixing process is dominated by diffusion because the twirling effect is weak. From Fig. 6(c), one may conclude that wf plays a critical role in determining the mixing efficiency. When all the slugs reach the outlet of the T-junction, the mixing index for wf = 0.30 exceeds 90%, which is approximately three times of that for wf = 0.84.

Now a question is raised: how does wf influence the twirling effect? As discussed in Sec. III C, wf primarily controls the slug sizes, thus it determines the amount of fluid inside the slugs to be distributed by the twirling effect. Long slugs require strong twirling effects because they have more fluid to be distributed than that of short slugs. If the twirling effect of certain strength is sufficient to distribute all the fluid in short slugs, it may only influence a small portion of fluid in long slugs. Besides the amount of fluid, wf play roles in influencing the twirling time as well. If U_total remains constant, the velocity of the dispersed phase increases with the increase of wf, thus the dispersed phase takes less time to reach the lower-side channel wall. Consequently, the twirling effect has less time to distribute the fluid. Interestingly, longer twirling time does not essentially lead to better mixing for small slugs. From Fig. 6(c), one may notice that wf = 0.3 provides even better mixing than that of wf = 0.14 or 0.20. The slug sizes in these three cases are similar, but the long twirling time in the cases of wf = 0.14 and 0.20 leads to a so called “over-twirling” phenomenon. As shown in Fig. 6(b), the dye solution is now resembled more in the upper part instead of the lower parts. The corresponding mixing indices of wf = 0.14 indicate that the mixing process reaches to the diffusion dominated regime quickly because the majority of the fluid deviates from the recirculating flow path. Since “over-twirling” undermines the mixing inside slugs, it is undesired for those applications where efficient mixing is critical.

Based on their study of using a droplet-based microfluidic system to conduct the in situ kinetics measurement, Nichols et al. have pointed out that the mixing inside the slugs must be in “kinetic regime.”²³ The “kinetic regime” defined by them is similar as the convective dominated regime shown in this section. Our results indicate that such regime can only be achieved by using short slugs with wf < 0.4. In addition, excessively short slugs, i.e., with wf < 0.3, have relatively low efficiencies as well, which reminds the researchers that it is important to maintain a proper wf value or slug sizes in order to avoid the “over-twirling” phenomena.

As discussed in Sec. III D 2, the twirling effect is a vortex induced by the shear force between the two phases. Since the mixing of dye and water solutions inside long slugs is not efficient, one may ask if the mixing efficiencies in long slugs can be improved by using high total velocity U_total. When U_total is increased, the shear force exerted by the continuous phase is strengthened, which consequently enhances the twirling effect as well as the recirculating flow inside slugs. To understand the impact of U_total on mixing, additional simulations were carried out with U_total = 20 and 50 mm/s and wf = 0.20 and 0.52. As discussed in Sec. III C, the slug sizes are almost identical for a particular wf value; therefore, the effect of the amount of fluid on mixing is eliminated by comparing the slugs with same wf but different values of U_total.

The contour plots of the dye concentration shown in Fig. 7(a) indicate that the mixing efficiencies in small slugs (wf = 0.20) cannot be improved significantly by increasing U_total. When U_total = 20 mm/s is used, the long twirling time results in a good initial mixing but “over-twirling” phenomenon as well. As shown in Fig. 7(b), if U_total is increased to 50 mm/s, the initial mixing index decreases approximately 5% than that of U_total = 20 mm/s because of the reduced twirling time; however, as the twirling effect is strengthened by increasing U_total, the solutions are distributed more effectively in the front and rear parts of the slugs. After the slugs is separated, the enhanced recirculating flow can mix the solutions rapidly in the upper and lower parts, which is evident by the large slope of the mixing efficiency curve of U_total = 50 mm/s in the convection dominated regime. When the slugs approaches the outlet of the channel, the mixing indices increase about 3% higher than that of U_total = 20 mm/s. On the other hand, the enhanced twirling effect results in “over-twirling” as well. As seen in Fig. 7(a), the majority of the dye solution is distributed in the upper part of the slug; consequently, the mixing efficiency is undermined by enhanced twirling strength to some extent.

If wf is increased to 0.52, the obtained slug size is close to two times than that of wf = 0.20. As seen from Fig. 7(b), mixing is not improved by increasing U_total. Compared to the large amount of fluid inside the slug, the amount that can be distributed by twirling is still small even though the twirling effect is strengthened. The mixing indices corresponding to U_total = 20 mm/s and 50 mm/s are close to each other when the slugs reach the outlet.

Based on these results, one may conclude that mixing inside slugs cannot be improved significantly by increasing U_total. Since the short slugs already have efficient mixing performance because of their small sizes, increasing the shear force has limited effect on improving mixing. On the other hand, long slugs have relatively large amounts of fluid to be distributed through the twirling effect. Even though the twirling strength is enhanced by increasing U_total, the impact on mixing is still insufficient compared to the large amount of fluid inside slugs. Therefore, it is more effective to reduce the slug sizes than increasing the total velocities in terms of improving the mixing performance inside slugs.

In this work, we have adopted the VOF method coupled with the species transport model to quantitatively investigate the dynamic mixing process of dye and water solutions inside slugs generated by a microfluidic T-junction. After eliminating the potential errors resulting from the grid resolution of the computational mesh and the entrance lengths of the channel inlets, we have validated the numerical model by comparing the obtained slug lengths and periods with the experimental data reported by Tice et al.²² at various operating conditions. The simulation results have shown that the VOF method can capture the slug formation dynamics and the internal velocity fields inside the slugs successfully. The mixing indices of slugs have been employed to understand the evolutions of mixing efficiencies inside the microfluidic channels. When the dispersed phase enters the junction region, the shear effect exerted by the continuous phase induces the twirling effect, which distributes the solutions transversely to the front and rear parts of the slugs. After the slugs are detached from the dispersed phase, the recirculating flow inside the slugs can mix the solutions efficiently. The twirling effect is noticeably influenced by the slug sizes. For long slugs, the twirling effect cannot distribute effectively because of the large amount of fluid contained in the slugs. In contrast, the twirling effect appears to be very effective for short slugs as the amount of fluid is small, even though long twirling time in certain operating conditions may cause “over-twirling” phenomenon that undermines the mixing. Although one can enhance the twirling strength by increasing the total velocity, its impact on improving mixing is not as effective as reducing slug sizes. This study provides fluid mechanical insights to the design and operation of droplet-based microfluidics where efficient mixing inside slugs is essential in such cases as nanomaterial synthesis, studies of in situ reaction kinetics, cell activation in cryopreserved cells, and protein crystallization.

The authors would also like to thank the HPC and LONI facility at Louisiana State University for computational support. This material is based upon work supported as part of the Center for Atomic Level Catalyst Design, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, and Office of Basic Energy Sciences under Award No. DE-SC0001058.

In Sec. III A, we have demonstrated that the medium grid resolution is able to control the artificial diffusion caused by the numerical scheme with U_total = 50 mm/s. The capability of this grid to handle high flow rates was tested by the case with U_total = 150 mm/s and wf = 0.52. The corresponding slug lengths obtained from the medium and fine grid resolutions were 102 and 100 μm, while the periods were 144 μm and 142 μm, respectively. As the variations of the slug length and period obtained from the medium and fine grid resolutions are within 5% of each other, the simulations with the medium grid resolution are believed to be able to capture the breakup dynamics accurately.

The microfluidic channels used in the experiments have very long entrance lengths for the straight and side inlets. In the simulations, these inlet channels were truncated so as to save the computational cost. However, the entrance length must be sufficient to ensure fully developed flow before the streams reach the junction region. The empirical formula proposed by Dombrowski et al.⁴⁶ indicates that the required inlet length to provide fully developed flow can be calculated by

where L_e is the equivalent entrance length, d_h is the hydraulic diameter, and Re is the Reynolds number. For the highest flow rate used in this study (U_total = 150 mm/s), the required inlet lengths calculated by Eq. (A1) were Le_c = 33.21 μm for the continuous phase and Le_d = 39.23 for the dispersed phase. For the lowest flow rate (U_total = 20 mm/s), the required inlet lengths were Le_c = 31.96 μm and Le_d = 32.08 μm. These results suggest that the selected entrance lengths of 50 μm are sufficient to ensure fully developed flows prior to the T-junction. To further validate the sufficiency of the entrance lengths, two simulations were performed with the inlet lengths of 150 μm and 50 μm, respectively. The operating conditions were set as: U_total = 150 mm/s and wf = 0.2. As shown in Fig. 8, these two channel lengths provide almost identical simulation results: the averaged slug lengths are 65 μm and 64 μm; the corresponding periods are 218 μm and 220 μm. Therefore, the entrance length of 50 μm is believed to be sufficient to ensure fully developed flows prior to the T-junction. In this study, the entrance lengths were all set as 50 μm.

As mentioned in Sec. II, the wall adhesion force is taken into account by specifying a constant contact angle θ_w near wall. Several researchers have reported that the contact angle only plays roles in determining the flow patterns (i.e., stratified flow or slug flow), but the slug lengths are marginally influenced by the contact angle.^30,41,42,47 Since the continuous phase preferentially wets the channel wall, a thin film is developed between the moving slug and the rigid channel walls. Qian et al. recommended using a quasi-non-wetting condition that assigns a constant value to the contact angle instead of revealing the thin film.⁴⁸ Gupta et al. have used a very fine grid resolution to reveal the thin film, and they suggested using at least 5 grids inside the thin film to capture the correct film thickness.⁴⁹ Based on the operating conditions used by this study, a typical thin film thickness estimated from Bretherton correlation is about 3 μm,⁵⁰ thus a much finer grid resolution than the current one is required to capture the thin film. As this study focuses on the mixing inside slugs, revealing the thin film is not necessary for the aim of this work. Therefore, we hypothesized that the contact angle does not affect slug lengths. This hypothesis was validated by numerical simulations.

Virtually estimated from the experimental pictures, θ_w is close to 140°. To investigate the effect of θ_w on slug formation dynamics, numerical simulations were carried out with two θ_w values (θ_w = 140° and 180°) at U_total = 50 mm/s and various wf (wf = 0.14–0.84). As shown in Fig. 9, these two θ_w values result in almost identical slug lengths and periods for the explored wf range except for the cases with high wf (wf = 0.73 and 0.84). For low wf cases (wf = 0.14–0.6), the effect of wall adhesion force is small because the dispersed phase only contacts the channel wall for a very short time. If wf is large, the dispersed phase contacts the channel wall for a long time, which indicates that the adhesion force exerted by the wall has sufficient time to influence the slug breakup process. The simulation using θ_w = 140° exaggerates the wall adhesion force significantly for the cases with wf = 0.84. If θ_w is varied from 140° to 180° for the case with wf = 0.84, the obtained slug lengths and periods decrease with θ_w as shown in Fig. 10(c). No slugs are formed if the contact angle is 120° as shown in Fig. 10(a1). If θ_w is increased to 140°, the dispersed phase forms a long segment from the junction region to the breaking point where the slug is separated from the dispersed phase. Such breakup pattern is discrepant from that been observed in the experiments, thus it is believed to be resulted from the numerical artifact. If θ_w is further increased to 150°, the slug length and periods decrease significantly, and the breaking point moves back to the junction region. If θ_w is set as 180°, the slug length and period are reduced slightly, but the breaking point moves downstream again. Among these θ_w values, the slug length and period predicted by the contact angle of 150° match the experimental results best for this particular case.

When U_total was set over 100 mm/s, simulations using θ_w = 140° could not predict the correct slug breakup dynamics as shown in Fig. 4. The other values were assigned to θ_w to test if they were able to predict the correct breakup dynamics. As shown in Fig. 10(b), θ_w = 120° results in deformed slug shapes, which is discrepant from the experimental observations. If θ_w is increased from 140° to 180°, the slug lengths and periods are marginally affected by its values as seen in Fig. 10(c). The dispersed phase still forms a segment in the downstream of the junction; however, the lengths of these segments decrease with θ_w. The breaking point also moves towards the junction with the increment of θ_w values.

All these results indicate that the breakup dynamics cannot be predicted by using the constant contact angle model for the cases with large U_total. The high velocity gradient across the channel wall leads to significant different values of contact angle along the interface. An accurate model that takes account of the effect of three-phase contact line velocity is required for the cases with U_total > 100 mm/s. Since it is beyond the scope of this manuscript, we do not discuss it here.