Volume of Fluid (VOF) Method as a Suitable Method for Studying Droplet Formation in a Microchannel

Abstract

Microfluidics is a rapidly advancing field focused on optimizing microdevices for applications such as organ-on-a-chip systems and enhancing laboratory analyses. Understanding the physical parameters of droplet generation is crucial for these devices. Computational fluid dynamics (CFD) techniques are essential for providing insights into the limitations and efficiency of numerical methods for studying fluid dynamics and improving our understanding of various application conditions. However, the influence of different numerical methods on the analysis of physical parameters in problems involving droplet generation in microchannels remains an area of ongoing research. This study implements the Volume of Fluid (VOF) method to investigate key physical parameters, including droplet size and the effect of the capillary number on fluid regimes, in droplet generation within a microchannel featuring a T-junction geometry. We compare the VOF method with the widely used Level Set Method (LSM) to evaluate its suitability for this context. The results show that the VOF method agrees with the LSM in fundamental outcomes, such as the reduction in droplet diameter as the flow rate ratio increases and the identification of the capillary number’s influence on fluid regime classification. The VOF method provides a clearer understanding of transitions between fluid regimes by detecting stages of non-uniformity in droplet size. It identifies a transition region between regimes with variations in droplet size, proving to be effective at mapping fluid flow regimes. This study highlights the potential of the VOF method in offering more detailed insights into instabilities and transitions between fluid regimes at the microscale.

1. Introduction

The development of microdevices has played a significant role in biomedical applications, particularly those related to the optimization of lab-on-a-chip devices, which have become essential for adapting classical laboratory methods to the microscale volumes of biological phenomena [1].

Systems that integrate microdevices are founded on the principles of microfluidics, a technology focused on the manipulation of small volumes of fluids within microchannels typically sized between a few tens and several hundreds of micrometers [2]. The study of microfluidics applied to microdevice-based systems for biotechnological applications can focus on analyzing the physical parameters of the dynamics of biological fluids flowing through microchannels. This analysis enables the understanding of how fluid behavior, dictated by physical parameters, relates to and influences the biological parameters essential for optimizing a properly functional system.

Droplet microfluidics is classified as a subcategory of microfluidics that is of great interest, as it allows for the creation, manipulation, and analysis of tiny droplets, which can serve as independent microreactors for various chemical and biological reactions [3]. In this context, the ability to produce monodisperse microdroplets plays a significant role in devices used for organ-on-a-chip applications [4].

The implementation of computational tools emerges as an important means to facilitate the study of microfluidics in microdevices. Computational fluid dynamics (CFD) is characterized by the application of robust computational methods to model fluid flow conditions, predicting aspects of heat, mass, and momentum transfer [5] by employing numerical analyses and data structures to identify and solve fluid flow problems [6].

Computational fluid dynamics (CFD) techniques enable the determination of critical parameters such as fluid velocity, pressure, surface tension, and density, which are challenging to ascertain through experimental techniques. The analysis provided by computational modeling of microfluidic phenomena under a wide range of experimental conditions allows for the extraction of relevant information for the development and optimization of microdevices, as well as their operating conditions and functionality [7].

For droplet-based microfluidic systems, analyzing critical parameters becomes a fundamental task, as it is necessary to understand and identify their influence on the dimensions and shapes of the generated droplets, considering that monodisperse droplets with specific sizes are crucial for different applications [8]. Thus, precise control over the dimensions and shapes of the microdroplets generated in microchannels is important, as it influences critical parameters for the functioning of microdevices [9]. As for the analysis of microdroplet generation, mapping the flow regime is also an important task, as it helps to determine the best conditions for forming microdroplets with the desired properties [10].

Many studies have applied computational and experimental methods to analyze physical parameters in the droplet generation process within microdevices [11]. Experimental investigations have examined the influence of fluid elasticity on droplet formation in T-junction geometries, indicating that elastic properties in low-viscosity fluids promote the generation of secondary droplets with varying sizes [12]. Complementary numerical analyses have studied the formation of both single and compound droplets in microfluidic flow-focusing devices, utilizing finite element methods with adaptive mesh refinement to evaluate viscoelastic effects on droplet behavior [13]. Furthermore, the dynamics of droplet formation in flow-focusing configurations were predicted using the Volume of Fluid (VOF) method, revealing qualitative agreement between simulations and experimental observations [14].

It is observed in the literature that many studies have implemented different computational methods to evaluate the physical properties of droplet formation in microchannels. However, few efforts have been made to comparatively assess the efficiency of different methods and to understand the similarities and discrepancies between them. Thus, the present work aims to understand the main differences between two well-established computational methods within the context of microfluidics through the study of droplet formation in a microchannel with T-junction geometry.

Building upon this foundation, we aim to analyze the effects of fluid velocity and capillary number ($Ca$) on the formation of microdroplets in a T-junction geometry channel using the Volume of Fluid (VOF) method. This includes identifying the influence of these parameters on droplet dimensions and their shapes, as well as examining the fluid regimes that dominate microdroplet formation. Alongside this, we seek to compare the VOF method with the Level Set Method (LSM) employed in the numerical–experimental study that served as our validation benchmark [11].

The distinct methods analyzed in this study exhibit unique potentialities and characteristics. While many studies highlight the Level Set Method as an excellent tool for analyzing droplet generation in microchannels due to its high-order accuracy and good mass conservation, we aim to explore how well the Volume of Fluid method adapts to microfluidic problems. By focusing on its potential for mass conservation, we hope to better understand phenomena and disturbances in the droplet generation process, particularly under conditions of fluid flow instability and variations in the capillary number ($Ca$).

2. Mechanisms of Droplet Generation

The formation of droplets in microchannels occurs due to the instability of the fluids present. Droplet generation typically occurs within three main regimes, known as squeezing, dripping, and jetting [15].

Microdroplets can be formed in a geometry that features a T-junction (Figure 1), where two perpendicular channels intersect. The main channel contains a continuous phase fluid, and the orthogonal channel is the inlet of a internal phase. At the beginning of droplet formation, the dispersed phase fluid flows into the main channel. In the following stage, this fluid forms a plug whose length is approximately the same as the width of the main channel. At this stage, there is an interaction between shear forces and interfacial forces.

The effect of interfacial tension and viscosity acting between two considered immiscible fluids is represented by the capillary number ($Ca$), which thus influences the conditions that determine the dominance of each regime in the droplet formation process, making it a determining factor in the characterization and control of droplet formation dynamics. This parameter can then be determined as

(1)$Ca={\displaystyle \frac{\mu U_c}{\sigma }}$

where U (m/s) is the average velocity of the carrier fluid, $\mu$ ($\mathrm{Pa}·\mathrm{s}$) is the dynamic viscosity, and $\sigma$ (N/m) is the surface tension. The squeezing and dripping mechanisms typically occur due to the influence of the physical confinement exerted by the microchannel walls. In the squeezing regime, observed for values of $Ca<{10}^{-2}$, surface forces dominate over shear stresses, causing the drop length to depend on the flow rate ratio of the fluids and the dimensions of the microchannel. The transition between the two regimes occurs for values of $Ca>{10}^{-2}$, where, in this case, the dripping regime becomes dominant, characterized by the shear forces overtaking the surface forces [16].

In the squeezing regime, the flow rate of the fluid that constitutes the dispersed phase ($Q_d$) is greater than the flow rate of the continuous phase fluid ($Q_c$), and thus, the droplet elongates along the main microchannel, presenting a part of its shape that tapers. The breakup of the droplet occurs with the detachment of the dispersed phase from the channel junction due to the pressure drop over the droplet in the main channel. In the dripping regime, the flow rate of the dispersed phase fluid ($Q_d$) is lower than the flow rate of the continuous phase fluid ($Q_c$). Thus, the breakup of the droplet occurs immediately after a portion of the dispersed fluid is formed that extends completely over the cross-section of the main channel, creating a kind of plug. For high $Ca$ values, shear forces become even more crucial. The continuous phase applies a shear stress on the growing droplet, causing it to detach near the junction at the upper wall, while still allowing the main channel’s cross-section to remain unblocked. This phenomenon is referred to as a shearing regime. Once detached, the droplet travels downstream within the main channel. Concurrently, the tip of the dispersed phase retracts to the end of the inlet channel and the cycle recommences. Therefore, the droplet size formed in the microchannel is governed by the balance between shear force and interfacial tension [11].

In this study, using computational fluid dynamics techniques, we employed a T-junction geometry to lead an analysis of properties such as the size and shape of microdroplets generated by the interaction between two fluids, oil (the continuous phase) and water (the dispersed phase), as well as a mapping of the flow regimes in which the microdroplets were generated. We considered $U_d$ as the velocity of the dispersed phase fluid and $U_c$ as the velocity of the continuous phase fluid. The flows (oil as continuous phase and water as dispersed phase) are considered two-dimensional and laminar. Both phases are considered to have constant fluid density and viscosity, and the capillary number ($Ca$) was determined based on the properties of the continuous phase.

3. Computational Method and Numerical Solution

Different methods can be used for numerical simulations of droplet formation under conditions of strong surface tension effects [17]. These methods can be classified as interface tracking or interface capture methods. Among the main methods applied are the Level Set Method, marker methods, the Phase Field method, and the VOF method, which can be combined with spatial and temporal discretization approaches, such as boundary integral methods, finite element methods, and finite difference schemes [18].

In general, computational methods for solving numerical simulations in fluid dynamics can be categorized into explicit interface tracking methods and implicit interface capture methods. The first group of methods involve tracking the interface using a moving grid or through the use of massless markers positioned along the interface. In addition, a moving grid coupled with a fixed grid can also be used. Among these methods, the Front Tracking method stands out [19].

Explicit methods reconstruct the interface by relying on a set of discrete points stored along the interface, with the resolution directly depending on the number of such points. Conversely, implicit methods employ a scalar function defined on a fixed grid to approximate the interface, where its position is inferred through an indicator function. As a result, these approaches do not provides the exact location of the interface.

In problems of dimensions n, explicit methods utilize a grid of dimensions $(n-1)$ to reconstruct the interface, whereas implicit methods require solving an equation of dimensions n throughout the computational domain. A notable advantage of the implicit approach is its ability to address more complex problems, including those involving topological changes. Prominent examples of implicit methods include the Volume of Fluid (VOF) method, the Phase Field method, the Lattice Boltzmann method, and the Level Set Method.

A brief comparison of the main methods used in multiphase fluid dynamics simulations is presented below. Particular attention is given to the Volume of Fluid (VOF) method, as it was chosen for use in this study.

3.1. Front Tracking Method

The Front Tracking method, classified as an explicit interface tracking method, was developed by the Glimm and Tryggvason groups, who developed a three-dimensional Front Tracking algorithm to address Rayleigh–Taylor instability issues [20]. In its development, Unverdi and Tryggvason also investigated the rising of single and multiple bubbles [21]. This method employs interconnected markers to dynamically reconstruct the interface grid, establishing a relationship where the front grid moves relative to the fixed global grid. This approach provides an accurate representation of the location and geometry of the interface. In this context, the interfacial tension is calculated on the front grid and subsequently transferred to the fixed global grid.

Assuming that the physical properties of the fluid are constant within each phase, the Navier–Stokes equations are computed using a fixed grid. A discrete delta function, applied along the smoothed interface, is used to incorporate the interfacial jump condition into the momentum equations. This multigrid scheme demonstrates the complexity inherent in the Front Tracking method. Furthermore, it lacks the ability to resolve the formation of thin films resulting from droplet rupture [19].

The Front Tracking method has been improved and developed to address volume conservation and multiphase problems [22]. Some studies have implemented the Front Tracking method to demonstrate that discontinuity curves in two-dimensional interfaces can be tracked in a one-dimensional manner [23]. Furthermore, the algorithm has been uniformly applied to n-dimensional problems. Researchers have sought to develop Front Tracking methods to handle unconnected points in meshes [24].

3.2. Level Set Method

The Level Set Method, classified as an implicit interface capturing method, was demonstrated through a novel algorithm by Osher and Sethian, who developed a specific methodology for tracking a moving interface within a fixed grid system [25]. The Level Set function $\varphi (\overrightarrow{x},t)$ is defined throughout the domain near the interface and represents the signed minimum distance from any point to the interface, measured along the normal direction. The interface between two fluids is implicitly represented by the function $\varphi (\overrightarrow{x},t)=0$. This means that the interface is described as the set of points where the Level Set function equals zero, allowing for a smoother and more continuous representation of interfaces, even in the presence of deformation. For values of $\varphi (\overrightarrow{x},t)>0$, the computational domain is referenced as being inside the interface—this is, on one side of the fluid. Conversely, for values of $\varphi (\overrightarrow{x},t)<0$, the domain is considered to be outside the interface [26].

The method has been developed over the years, with significant contributions made by Olsson and Kreiss [27], who utilized it for interface capturing. One advantage of this method is that numerical calculations can be performed on a fixed Cartesian grid without the neeed for parameterizing the interface [7]. Due to its strong performance in accurately delineating phase boundaries, the Level Set Method has been widely applied to problems involving moving interfaces and free boundaries within the context of microfluidics [11]. The method has also been implemented in two-dimensional and three-dimensional problems, such as particle encapsulation in droplets transported by fluid flow [28]; droplet solidification and spreading [29]; droplet falling down [30]; breakup into smaller droplets [31]; and topological changes in the interface [32].

3.3. Volume of Fluid Method

The Volume of Fluid (VOF) method, classified as an implicit interface capturing method, was developed and extensively discussed by researchers such as Scardovelli and Zaleski [33], as well as Hirt and Nichols [34]. These authors presented techniques for employing the VOF approach to capture interfaces within non-uniform mashes. The method is based on a volume fraction function, which is then defined as a volume fraction field across the entire computational domain. The VOF method has been widely used in areas beyond multiphase fluid problems, such as heat transfer [35], thermocapillary motion [36], mass transfer [37], and particle-laden flows [38].

In this study, we guided the computational fluid dynamics simulations using the VOF method. This method characterizes the interface by using the volume fraction of one fluid phase or component, represented as C. In this method, the volume fraction function is introduced and defined within a volume fraction field across the entire computational domain [39]. In the bulk phase, within multi-fluid computational cells, the volume fraction C satisfies $0<C<1$. Generally, the VOF model utilizes phase averaging to estimate the proportion of continuous and dispersed phases within each computational cell [19].

A variable $\alpha$ is introduced such that $\alpha =1$ indicates a cell fully occupied by the continuous phase; $\alpha =0$ corresponds to a cell entirely filled with the dispersed phase; and values of $\alpha$ between 0 and 1 denote cells containing the interface between both phases. The function of the volume fraction is

(2)${\displaystyle \frac{\partial \alpha }{\partial t}}+\overrightarrow{U}·\nabla \alpha =0$

and can be solved on a fixed grid. In order to solve the function, the approximate position of the interface needs to be found. Additionally, the interface reconstruction must be carried out to determine the weighted viscosity and density for the computational cells, thus making it possible to calculate the volume flux for the convective terms present in the governing equations [19].

The density $\rho$ and viscosity $\mu$ for the water and oil phases are determined based on a linear dependence [34]. For the continuous liquid (primary) phase (3) and dispersed phase (droplet) (4), it is described as

(3)$\rho ={\rho }_1\alpha +{\rho }_2(1-\alpha )$

(4)$\mu ={\mu }_1\alpha +{\mu }_2(1-\alpha )$

3.4. Governing Equations

The governing equations stem from the Navier–Stokes equation, which represents a continuous perspective of Newton’s second law per unit of volume. For a laminar and incompressive flow, it can be written as

(5)$\rho \left({\displaystyle \frac{\partial \overrightarrow{U}}{\partial t}}+\overrightarrow{U}·\nabla \overrightarrow{U}\right)=-\nabla p+\mu {\nabla }^2\overrightarrow{U}+\overrightarrow{f}$

Additionally, the continuity equation is also solved:

(6)$\nabla ·\overrightarrow{U}=0$

In the equations above, $\rho$ is the density (${\mathrm{kg}/\mathrm{m}}^3$), $\overrightarrow{U}$ is the velocity (m/s), p is the pressure (Pa), so $\nabla p$ represents its gradient, $\mu$ is the viscosity ($\mathrm{Pa}·\mathrm{s}$), and $\overrightarrow{f}$ represents the external forces applied to the fluid (${\mathrm{N}/\mathrm{m}}^3$).

3.5. Numerical Calculation of Surface Tension Force

Considering the Navier–Stokes formulation (Equation (5)), the term $\overrightarrow{f}$ incorporates the surface tension forces, which we denote as ${\overrightarrow{f}}^{\sigma }$. The surface tension force can be considered per unit volume. Numerical approximations for calculating surface tension forces have been developed over the years [40]. For the purposes of this article, we adopt the Continuous Surface Force (CSF) model, implemented through the resolution software used for the simulations.

Considering an elementary volume ($\mathrm{\Omega }$) intersected by two points A and B, the surface tension force is given by

(7)${\int }_{\mathrm{\Omega }}{\overrightarrow{f}}^{\sigma }dS={\oint }_A^B\sigma d\overrightarrow{t}$

where $\sigma$ is the surface tension and $\overrightarrow{t}$ is the unit tangent vector. Using the Frenet parameterization for curves, we have $d\overrightarrow{t}=\kappa \overrightarrow{n}ds$, where $\kappa$ is the curvature, $\overrightarrow{n}$ is the unit normal vector, and s is the curvilinear coordinate. Thus, we obtain the following:

(8)${\int }_{\mathrm{\Omega }}{\overrightarrow{f}}^{\sigma }dS={\oint }_A^B\sigma \kappa \overrightarrow{n}ds={\int }_{\mathrm{\Omega }}\sigma \kappa \overrightarrow{n}{\delta }_s$

where ${\delta }_S$ is the surface Dirac delta function, which is nonzero only at the interface.

Using the Continuous Surface Force model by Brackbill et al. [41], we can write the following relation:

(9)$\sigma \kappa \overrightarrow{n}{\delta }_s=\sigma \kappa \nabla H(\overrightarrow{x}-\overrightarrow{x_s})$

where H denotes the Heaviside function and $x_s$ represents the interface location, offering a smoothed Dirac delta approximation of the surface tension force. In the CSF model, a void fraction field $H_{ϵ}=C$ is employed to approximate the Heaviside function, with $ϵ=\Delta$ corresponding to the grid resolution, such that

(10)$\underset{ϵ\to 0}{lim}{H}_{ϵ}=H$

The curvature of the interface can be determined using the height function approach, proposed by Cummins et al. [42] and Popinet [43], which is based on the concept that an interface can be described and defined as a graph of a function in a local coordinate system. The process consists of three steps: determining the local orientation of the coordinates through the vector component $\overrightarrow{n}$; estimating the height function by summing the void fraction values; and finally, calculating the curvature using central differencing [40].

Starting from the principle that the present work seeks its analyses through two-dimensional simulations, some considerations must be taken into account. The surface tension should be defined along a line, rather than a surface, with the normal vector $\overrightarrow{n}$ and the tangent vector $\overrightarrow{t}$ adjusted in a two-dimensional plane. Furthermore, the surface tension force can be considered along the interface line, which can be represented as

(11)${\int }_{\Gamma }\sigma dL$

where $\Gamma$ is the interface line. The height function then describes the position of the interface line in one direction, and the curvature estimation can be made based on the orientation of the interface in a single plane. The Dirac delta function $\delta$ will also act along the interface line.

4. Computational Procedure and Simulational Setup

For this study, the CFD software Fluent was employed, which utilizes control-volume-based methods to discretize the governing equations into algebraic forms suitable for numerical solution. The geometry modeling used in the simulations was carried out using Ansys SpaceClaim (2024 R2) software and the mesh generation was carried out afterward, as shown in Figure 2. The quadrilateral mesh was used, providing more accurate calculations of surface tension effects compared to triangular or tetrahedral meshes [44]. The dimensions of the channels can be seen in Figure 3.

4.1. Mesh Independece Study

To ensure the adequacy of the number of nodes in the droplets and computational performance, we conducted a study on mesh independence. To this end, we performed a series of simulations with different mesh element sizes, analyzing how the mesh dimensions would alter or influence the droplet size parameters. The mesh independence study aims to identify the mesh element size from which variations in physical parameters (such as droplet size, in this case) become minimized despite further changes in element size. This makes it possible to employ a computational mesh that provides good accuracy while maintaining a reasonable computational cost. Thus, keeping all other variables fixed, we conducted simulations using different mesh element sizes.

To this end, we chose to test meshes with element sizes of 2.5 $\mathsf{\mu }$m, 5.0 $\mathsf{\mu }$m, 7.5 $\mathsf{\mu }$m, and 10 $\mathsf{\mu }$m in order to measure the computational effort. Figure 4 shows the different results for the droplet diameter for the various mesh element sizes. The analysis allowed us to consider the mesh size of 5 $\mathsf{\mu }$m as the most efficient, as it provided considerably good accuracy while maintaining a reasonable computational effort.

The analysis of the results presented in the graph allows us to conclude that the mesh element size can directly influence the studied physical parameter—in this case, the droplet size. Although smaller mesh sizes, such as 2.5 $\mathsf{\mu }$m, provide excellent precision in defining the interface, allowing an accurate estimation of droplet size, such mesh configurations require significant computational effort. However, larger meshes, such as those beyond 7.5 $\mathsf{\mu }$m, reduce computational cost at the expense of the quality and definition of the fluid interface, thus reducing the accuracy of the simulation results. Therefore, we selected the 5 $\mathsf{\mu }$m mesh, as it demonstrated an optimal balance, offering reliable simulation results with negligible differences compared to more refined meshes while maintaining a reasonable computational effort.

4.2. Simulational Setup

The validity of our model was established through comparison with the results of Bashir et al. (2011) [11]. For this validation, the computational setup parameters and fluid physical properties were kept consistent with those utilized in that study. To perform the simulation, we defined the initial conditions of phase interactions based on anti-diffusion interfacial modeling, assuming the surface tension coefficient to be constant and equal to $0.005$ N/m. The inlet conditions were derived from the initial velocities of the fluids. A no-slip condition was applied on the walls of the geometry, and the contact angle between the phases, which constitutes the wall adhesion, was set to 135°. The solution method was characterized by the coupled with volume fractions scheme for the pressure–velocity coupling scheme. The spatial discretization was carried out using the Green–Gauss Cell-Based Gradient method, while PRESTO! was used for the pressure discretization scheme. The second-order upwind scheme was used for the discretization of the momentum equation. The simulations reached convergence when the time step accuracy was within ${10}^{-5}$, allowing for the creation of a database.

For all the simulations performed, we characterized the fluids with the physical properties presented in Table 1.

In order to understand the behavior of the fluids with respect to the generation and formation of microdroplets, we conducted simulations by varying the velocity of the continuous phase fluid ($U_c$) so that the ratio $U_c/U_d$ ranged from 1 to 12, considering the velocity of the dispersed phase fluid ($U_d$) fixed at $U_d=0.012$ m/s. Thus, we also varied the capillary number within this range. These values were chosen to ensure the compatibility of the present study with the work of Bashir et al., which served as a reference for data comparison [11].

5. Model Validation

As presented, considering the existence of different methods and simulation implementation approaches for the analysis of fluid dynamics behavior in microdevices, in this study, we aim to compare and validate our analyses through the work of Bashir et al., who implemented the Level Set Method (LSM), which is distinct method from the one used in the present study. Thus, we seek to identify the discrepancies resulting from the differences between the approaches, aiming for a more specific analysis of the Volume of Fluid (VOF) method’s implementation. In this sense, we propose an analysis of how the method employed in this study prioritizes physical parameters and under which conditions it differs from the Level Set Method, which was experimentally validated through the work of Bashir et al. [11].

Bashir et al. [11] presented an experimental investigation of droplet generation in water and oil phases using a T-junction geometry, with the same dimensions for the lateral and main channels as those used in our study. They identified a certain divergence between the sizes of droplets generated experimentally and those analyzed through simulations at low velocity ratio values ($U_c/U_d$). However, for higher flow rates, this difference diminishes. This discrepancy may be due to variations in the experimental physical properties. Moreover, performing 2D simulations tends to fail in accurately capturing the forces of tension acting on the droplet. Figure 5 shows a graphical comparison of the relationship between the dimensionless droplet length, given by the ratio of its diameter ($L_d$) to the channel width ($W_c$) in which the microdroplets are generated, and the velocity ratio ($U_c/U_d$).

In this study, we conducted computational simulations considering a viscosity ratio ${\eta }_d/{\eta }_c$ close to 0.8 so that the results obtained could be compared with the experimental model provided by the literature, thus validating the analyses carried out. It can be observed that, for low fluid flow rate values, there is a certain divergence in the results (Figure 3), a phenomenon also reported by Bashir et al. As the velocity ratio values increase, the numerical data obtained through the simulations in this study converge, to some extent, with the experimental data of the model on which this study is based. We can attribute the observed differences in droplet size behavior to the numerical methods used. While the Level Set Method, employed by Bashir et al., uses a distance function to define the interface, the Volume of Fluid (VOF) method uses a volume fraction function to represent the interface between phases. In this way, VOF is considered a more robust method in terms of mass conservation, while the Level Set Method may offer better accuracy in the shape of the interface but introduces errors with respect to mass conservation [45].

6. Results and Discussion

The analyses and subsequent discussions aim to emphasize the physical aspects of the droplet generation process in the studied microchannel by comparing the implementation of the Volume of Fluid (VOF) method with the numerical results obtained by Bashir et al. using the Level Set Method (LSM) [11]. As will be identified, the method implemented in this study enabled a specific analysis of regions of unstable droplet formation, which allowed the development of a more complete phase diagram of the fluid regime, also considering the transition zones.

For a T-junction microfluidic channel geometry, a pressure gradient produced by the shear stress forces of the continuous phase is responsible for the behavior of the dispersed fluid, which tends to break at the junction of the geometry and eventually forms a droplet [46]. The ways in which the droplet is generated, as well as its behavior along the channel, can be influenced by parameters such as fluid velocities, influencing the fluid flow rate ratio, as well as the viscosity ratio of the fluids. The following results show the effects on droplet formation, with respect to their sizes and fluid flow regimes, by altering the fluid flow rates.

6.1. Effects of Flow Velocity Ratio on Droplet Generation

A series of twelve simulations were conducted to validate our analyses on the formation of microdroplets in terms of their size through the presented experimental model in order to analyze the results obtained from the implementation in our simulation model. Thus, we varied the fluid velocity ratio from 1 to 12, keeping the velocity of the dispersed phase fluid (water) constant at $U_d=0.012\mathrm{m}/\mathrm{s}$ and altering the velocity of the continuous phase fluid (n-dodecane oil). The correlation between the formed droplet size and the velocity ratio is shown in Figure 6.

As can be observed, in general, the droplet size, regarding its diameter dimensions, decreases as the velocity ratio of the fluids increases. Similarly to the results presented by Bashir et al., the range of difference between the analyzed sizes decreases at high velocity ratios, whereas for low values, the differences become more significant. It is noted that the discrepancies observed when comparing our results with those from the numerical validation model occur during moments of fluid regime transition, as will be explained in the analysis of the capillary number. This can be explained by the fact that the method used in our analyses tends to show better performance in conserving the mass of the simulated system [45].

The discrepancies found, particularly under low-velocity-ratio conditions, are related to fluid regime transition moments, which explain a certain irregularity in droplet size. During transitions between regimes, a significant variation in droplet size may occur, as there is greater instability in the control over droplet separation, contributing to a broader distribution of droplet sizes. In these transitions, a combination of forces begins to act in a variable manner, further corroborating the irregularity [47].

Figure 7 shows the dependence of fluid regimes on the fluid velocities. We observe three regimes being formed. The first regime, observed for relatively low velocity ratios, is the squeezing regime, which has the pressure force as the cause of droplet breakup; the generated droplets are large and often irregular in size due to the longer detachment time and higher interfacial tension values [48]. The second regime observed is the dripping regime, in which a droplet is formed near the junction between the channels of the geometry; the velocity of the dispersed phase fluid is fast enough that the viscous stress exerted by the continuous phase fluid exceeds the interfacial tension. Thus, to comply with the mass conservation principle, the liquid thread begins to thin at the point of highest velocity. In this situation, a high capillary pressure is generated at the end of the fluid stretching, breaking the dispersed phase and, consequently, forming a droplet. The third regime observed, more subtly, is the jetting regime, which occurs with the increase in the velocity of the continuous phase fluid. The droplet tends to be formed at the end of a jet of the dispersed phase fluid due to the high flow inertia of this phase, causing the breakup to occur due to the instability of the jet [39]. A transition region from the squeezing to the dripping regime can also be observed. In this region, the droplets are not formed uniformly and do not tend to display a regular breakup mode.

6.2. Effects of Capillary Number on Droplet Size

By altering the fluid flow ratio through changes in the continuous fluid velocity ($U_c$), we also modify the capillary number for each simulation. From this, we could observe the behavior of the fluid phases with respect to droplet formation. Initially, as shown in Figure 8, it can be observed that the droplet size decreases as the capillary number increases, similar to what is observed with the fluid velocity ratios, since, under our conditions, the only factor that alters the capillary number is the variation in the velocity of the continuous fluid.

To conduct a more detailed analysis, we identified the moment of breakup that generated the droplet in the main microchannel The visual analysis of droplet breakup and formation enables a classification of the identified regimes according to the associated capillary number. Figure 9 shows the breakup process and the subsequent formation of droplets in the main channel concerning the capillary number ($Ca$) for each of the twelve simulations performed.

For relatively low capillary number values ($Ca$≤ 0.012), the presence of the squeezing regime is identified, which, as explained, results in the formation of larger and wider droplets, often without a well-defined contour. The time for the detachment of droplets is longer and there are irregularities between the interface of fluids, often producing droplets with poorly and non-uniform sizes [49]. For capillary number values within the range of $0.0012<Ca\le 0.023$, a transition stage between the regimes is observed. The droplets are formed irregularly, and the moments of droplet breakup exhibit different characteristics. In these transition regions, the pressure drop reaches a critical value, causing the transition from the squeezing regime to the dripping regime. It is also observed that, in this region, the droplet size does not continue to show a local decreasing trend; there is a slight increase in its diameter before returning to the trend of decrease. For capillary number values $0.023<Ca\le 0.032$, we observe droplet formation through the dripping regime, where droplets are generated near the junction of the geometry. In this regime, the role of shear force is significant for breakup. The droplets are typically more monodisperse and homogeneous compared to the squeezing regime and the transition phase [50]. Finally, for relatively high capillary number values ($Ca>0.032$), we identify the jetting regime. In this regime, droplets can be produced due to a Rayleigh–Plateau instability (i.e., capillary instability). The droplets are produced at the end of an extension of the dispersed phase fluid [51,52].

The mechanisms of droplet breakup and formation in the dripping and jetting regimes are different, influencing the shape and size of the generated droplets. In the dripping regime, the dispersed phase is influenced by several forces, including pressure, interfacial tension, viscous forces, and inertia. At the junction, the interfacial tension opposes droplet detachment, whereas pressure and viscous forces drive the dispersed phase downstream. These forces act on three distinct levels: initially, the dispersed phase enters the junction, inducing a pressure drop and velocity gradient in the continuous phase; subsequently, the droplet elongates within the main channel, but its growth is restrained by interfacial tension; and finally, when the dispersed phase obstructs the junction, the pressure drop reaches a critical threshold, leading to thinning of the elongated droplet. The breakup then occurs after the detachment of the dispersed phase from the junction. Thus, droplets are formed near the junction and attain their stable shape [53]. The jetting regime produces even smaller droplets due to the narrower elongation of the dispersed phase and lower values of interfacial tension and shear [54].

Figure 10 allows the visualization of different flow regimes in which fluids act during droplet generation, relating them to the variation in capillary number. It is possible to visually detect and classify, according to the capillary number and fluid flow rate of each regime, the regions where each type of breakup and droplet generation predominates. It can be observed that in the squeezing regime, the droplet size decreases as the capillary number increases. After reaching a certain value of $Ca$, there is a transition between regimes, and at this stage, the droplet size no longer follows a uniform trend of decrease or increase. Upon entering the dripping regime, the droplet size decreases as the capillary number increases. This decreasing trend continues into the jetting regime, but with a progressively slower rate of decrease.

Based on the analyses enabled by the results obtained in this study, it is possible to comparatively establish qualitative characteristics of the different evaluated methods, namely the VOF method, which was the primary approach employed in this research, and the LSM method, widely used in other studies and serving as numerical validation and comparison for this work. Table 2 provides a comparative description of the characteristics of the analyzed methods.

7. Conclusions

In this study, we aimed to understand how physical parameters, such as fluid velocity ratio and capillary number, affected the formation and generation of microdroplets in microchannels. To this end, we employed a numerical and simulation-based method that, once validated against literature data, provided crucial information for understanding the conditions of droplet formation. Initially, we analyzed how the velocity ratio between the dispersed and continuous phases ($U_c/U_d$) influenced the droplet size. We observed an inverse relationship, where droplet diameter decreased as the velocity ratio increased. Additionally, we examined the influence of the capillary number on the characteristics of different formation regimes. The employed simulation method, the Volume of Fluid (VOF) method, proved effective in identifying a transition region between formation regimes, characterized by instability in droplet sizes. The approach enabled a more detailed description of the mapping of fluid regimes.

Despite not offering the same level of accuracy in resolving the fluid interface as the Level Set Method, the VOF method demonstrated efficiency and proved to be suitable for analyzing microfluidic problems. It emerged as a viable alternative for numerical simulations in this field, balancing computational cost and accuracy. Moreover, the method was efficient in presenting key physical parameters, such as droplet size and formation regime characteristics, which are crucial for understanding the droplet generation process. The information obtained from this study contributes to a better understanding of the implementation of the VOF method in the numerical and simulation-based analysis of droplet formation in microchannels, with potential future applications in the development of microdevices.

Author Contributions

Conceptualization, F.S.P.d.S. and P.N.L.-F.; methodology, F.S.P.d.S. and P.N.L.-F.; validation, F.S.P.d.S. and P.N.L.-F.; formal analysis, F.S.P.d.S. and P.N.L.-F.; investigation, F.S.P.d.S. and P.N.L.-F.; writing, F.S.P.d.S.; writing—review and editing, F.S.P.d.S. and P.N.L.-F. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Schematic of the t-junction microchannel. A lateral channel is used for the inlet of the dispersed phase, while the continuous phase flows through the main channel. At the junction, the two phases converge, initiating the droplet formation process. The arrows in the schematic indicate the direction of the flows.

Figure 2 Meshing applied to the modeled geometry.

Figure 3 Dimensions of the channels in the T-junction geometry used in the simulations. The inlet channel for the dispersed phase has a width of $60\mathsf{\mu }\mathrm{m}$, while the inlet channel for the continuous phase has a width of $100\mathsf{\mu }\mathrm{m}$.

Figure 4 Mesh independence study. A 5 $\mathsf{\mu }$m mesh was selected as it offered a good balance between accuracy and computational cost, with no significant differences in the analyzed parameter compared to finer meshes.

Figure 5 Comparison of droplet length for flow velocity ratios between the present study and the validation model (Bashir et al., 2011 [11]), which includes experimental and numerical data.

Figure 6 Dimensionless droplet length as a function of the velocity ratio ($U_c/U_d$). A general decreasing trend in droplet size is observed as $U_c/U_d$ increases, except near $U_c/U_d\approx 4$, where an irregular pattern emerges with an unexpected increase in droplet length.

Figure 7 Phase diagram illustrating the different droplet formation regimes as a function of the continuous phase velocity ($U_c$) for a fixed dispersed phase velocity ($U_d=0.012\mathrm{m}/\mathrm{s}$). The regimes are defined over the following intervals of $U_c$: squeezing regime ($0.015\le U_c<0.060\mathrm{m}/\mathrm{s}$), transition regime ($0.060\le U_c<0.085\mathrm{m}/\mathrm{s}$), dripping regime ($0.085\le U_c<0.115\mathrm{m}/\mathrm{s}$), and jetting regime ($0.115\le U_c\le 0.140\mathrm{m}/\mathrm{s}$).

Figure 8 Dimensionless droplet length ($L_d/W_c$) as a function of the capillary number (Ca). The results indicate a general decrease in droplet length as the capillary number increases. A slight increase in droplet length is observed around $Ca\approx 0.012$, representing a local deviation from the overall decreasing behavior.

Figure 9 Droplet break up for different capillary numbers. The breakup moments of the dispersed phase at the junction were captured under different capillary number ($Ca$) conditions. These observations provide valuable insights into the droplet formation regimes as a function of the capillary number.

Figure 10 Phase diagram of fluid regimes depending on the capillary number. The squeezing regime is observed in the range of $0.003<Ca<0.012$. A transition region, associated with unstable droplet formation, occurs in the range of $0.012<Ca<0.023$. The dripping regime appears within $0.023<Ca<0.032$, and for $Ca>0.032$, the system enters the jetting regime.

Table 1

Physical properties of the fluids used in the simulations.

Fluid	Density ($\mathit{\rho }$, kg/$m^3$)	Viscosity ($\mathit{\eta }$, mPa · s)
Water	1000	1.0
n-Dodecane Oil	750	1.34

Table 2

Qualitative comparison between the VOF method and LSM based on the present study and the numerical validation by Bashir et al. [11].

Criterion	VOF (Present Study)	LSM (Bashir et al.) [11]
Interface Representation	Volume fraction function	Distance function
Mass Conservation	High (better conservation)	Limited (mass loss observed)
Droplet Size Prediction at Low $U_c/U_d$	Underestimation compared to experimental data	Closer to experimental data
Droplet Size Prediction at High $U_c/U_d$	Improved agreement with experiments	Good agreement with experiments
Identification of Regime Transitions	Effective in identifying transition and instability regions	Less detailed in transition mapping
Suitability for Microfluidics	Suitable and efficient for analyzing physical parameters	Higher accuracy in interface shape

References

1. Castillo-León, J. Microfluidics and lab-on-a-chip devices: History and challenges. Lab-on-a-Chip Devices and Micro-Total Analysis Systems: A Practical Guide; Springer: Cham, Switzerland, 2014; pp. 1-15.

2. Whitesides, G.M. The origins and the future of microfluidics. Nature; 2006; 442, pp. 368-373. [DOI: https://dx.doi.org/10.1038/nature05058] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16871203]

3. Taly, V.; Kelly, B.T.; Griffiths, A.D. Droplets as microreactors for high-throughput biology. ChemBioChem; 2007; 8, pp. 263-272. [DOI: https://dx.doi.org/10.1002/cbic.200600425] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17226878]

4. Yap, Y.F.; Tan, S.H.; Nguyen, N.T.; Murshed, S.S.; Wong, T.N.; Yobas, L. Thermally mediated control of liquid microdroplets at a bifurcation. J. Phys. D Appl. Phys.; 2009; 42, 065503. [DOI: https://dx.doi.org/10.1088/0022-3727/42/6/065503]

5. Malekjani, N.; Jafari, S.M. Simulation of food drying processes by Computational Fluid Dynamics (CFD); recent advances and approaches. Trends Food Sci. Technol.; 2018; 78, pp. 206-223. [DOI: https://dx.doi.org/10.1016/j.tifs.2018.06.006]

6. Minniti, F.; Fiamma, V.; Barbaro, G.; Puntorieri, P. Numerical Simulation Of Velocity And Pressure Field Around A Cylinder At Moderate Reynolds Number. Technology; 2019; 10, pp. 47-54. [DOI: https://dx.doi.org/10.34218/IJARET.10.4.2019.006]

7. Hernández-Cid, D.; Pérez-González, V.; Gallo-Villanueva, R.; González-Valdez, J.; Mata-Gómez, M. Modeling droplet formation in microfluidic flow-focusing devices using the two-phases level set method. Mater. Today Proc.; 2022; 48, pp. 30-40. [DOI: https://dx.doi.org/10.1016/j.matpr.2020.09.417]

8. Li, Y.; Wang, K.; Xu, J.; Luo, G. A capillary-assembled micro-device for monodispersed small bubble and droplet generation. Chem. Eng. J.; 2016; 293, pp. 182-188. [DOI: https://dx.doi.org/10.1016/j.cej.2016.02.074]

9. Teh, S.Y.; Lin, R.; Hung, L.H.; Lee, A.P. Droplet microfluidics. Lab A Chip; 2008; 8, pp. 198-220. [DOI: https://dx.doi.org/10.1039/b715524g]

10. Mastiani, M.; Seo, S.; Jimenez, S.M.; Petrozzi, N.; Kim, M.M. Flow regime mapping of aqueous two-phase system droplets in flow-focusing geometries. Colloids Surfaces A Physicochem. Eng. Asp.; 2017; 531, pp. 111-120. [DOI: https://dx.doi.org/10.1016/j.colsurfa.2017.07.083]

11. Bashir, S.; Rees, J.M.; Zimmerman, W.B. Simulations of microfluidic droplet formation using the two-phase level set method. Chem. Eng. Sci.; 2011; 66, pp. 4733-4741. [DOI: https://dx.doi.org/10.1016/j.ces.2011.06.034]

12. Husny, J.; Cooper-White, J.J. The effect of elasticity on drop creation in T-shaped microchannels. J. Non-Newton. Fluid Mech.; 2006; 137, pp. 121-136. [DOI: https://dx.doi.org/10.1016/j.jnnfm.2006.03.007]

13. Zhou, C.; Yue, P.; Feng, J.J. Formation of simple and compound drops in microfluidic devices. Phys. Fluids; 2006; 18, 092105. [DOI: https://dx.doi.org/10.1063/1.2353116]

14. Davidson, M.R.; Harvie, D.J.; Cooper-White, J.J. Flow focusing in microchannels. ANZIAM J.; 2004; 46, pp. C47-C58. [DOI: https://dx.doi.org/10.21914/anziamj.v46i0.946]

15. De Menech, M.; Garstecki, P.; Jousse, F.; Stone, H.A. Transition from squeezing to dripping in a microfluidic T-shaped junction. J. Fluid Mech.; 2008; 595, pp. 141-161. [DOI: https://dx.doi.org/10.1017/S002211200700910X]

16. Zhu, P.; Wang, L. Passive and active droplet generation with microfluidics: A review. Lab A Chip; 2017; 17, pp. 34-75. [DOI: https://dx.doi.org/10.1039/C6LC01018K] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27841886]

17. Liovic, P.; Francois, M.; Rudman, M.; Manasseh, R. Efficient simulation of surface tension-dominated flows through enhanced interface geometry interrogation. J. Comput. Phys.; 2010; 229, pp. 7520-7544. [DOI: https://dx.doi.org/10.1016/j.jcp.2010.06.034]

18. Hooper, R.; Cristini, V.; Shakya, S.; Lowengrub, J.S.; Derby, J.J.; Macosko, C.W. Modeling multiphase flows using a novel 3D adaptive remeshing algorithm. Trans. Eng. Sci.; 2001; 30, 34.

19. Liu, J.; Trung Nguyen, N. Numerical simulation of droplet-based microfluidics–A review. Micro Nanosyst.; 2010; 2, pp. 193-201. [DOI: https://dx.doi.org/10.2174/1876402911002030193]

20. Glimm, J.; Grove, J.W.; Li, X.L.; Shyue, K.M.; Zeng, Y.; Zhang, Q. Three-dimensional front tracking. SIAM J. Sci. Comput.; 1998; 19, pp. 703-727. [DOI: https://dx.doi.org/10.1137/S1064827595293600]

21. Unverdi, S.O.; Tryggvason, G. A front-tracking method for viscous, incompressible, multi-fluid flows. J. Comput. Phys.; 1992; 100, pp. 25-37. [DOI: https://dx.doi.org/10.1016/0021-9991(92)90307-K]

22. Dijkhuizen, W.; Roghair, I.; Annaland, M.V.S.; Kuipers, J. DNS of gas bubbles behaviour using an improved 3D front tracking model—Model development. Chem. Eng. Sci.; 2010; 65, pp. 1427-1437. [DOI: https://dx.doi.org/10.1016/j.ces.2009.10.022]

23. Mao, D.k. Towards front-tracking based on conservation in two space dimensions II, tracking discontinuities in capturing fashion. J. Comput. Phys.; 2007; 226, pp. 1550-1588. [DOI: https://dx.doi.org/10.1016/j.jcp.2007.06.004]

24. Gois, J.P.; Nakano, A.; Nonato, L.G.; Buscaglia, G.C. Front tracking with moving-least-squares surfaces. J. Comput. Phys.; 2008; 227, pp. 9643-9669. [DOI: https://dx.doi.org/10.1016/j.jcp.2008.07.013]

25. Osher, S.; Sethian, J.A. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. J. Comput. Phys.; 1988; 79, pp. 12-49. [DOI: https://dx.doi.org/10.1016/0021-9991(88)90002-2]

26. Peng, D.; Merriman, B.; Osher, S.; Zhao, H.; Kang, M. A PDE-based fast local level set method. J. Comput. Phys.; 1999; 155, pp. 410-438. [DOI: https://dx.doi.org/10.1006/jcph.1999.6345]

27. Olsson, E.; Kreiss, G. A conservative level set method for two phase flow. J. Comput. Phys.; 2005; 210, pp. 225-246. [DOI: https://dx.doi.org/10.1016/j.jcp.2005.04.007]

28. Yap, Y.F.; Chai, J.; Wong, T.N.; Nguyen, N.T.; Toh, K.C.; Zhang, H.; Yobas, L. A procedure for the motion of particle-encapsulated droplets in microchannels. Numer. Heat Transf. Part B Fundam.; 2008; 53, pp. 59-74. [DOI: https://dx.doi.org/10.1080/10407790701632485]

29. Zheng, L.; Zhang, H. An adaptive level set method for moving-boundary problems: Application to droplet spreading and solidification. Numer. Heat Transf. Part B Fundam.; 2000; 37, pp. 437-454. [DOI: https://dx.doi.org/10.1080/10407790050051137]

30. Ni, M.J.; Abdou, M.; Komori, S. A variable-density projection method for interfacial flows. Numer. Heat Transf. Part B Fundam.; 2003; 44, pp. 553-574. [DOI: https://dx.doi.org/10.1080/716100497]

31. Croce, R.; Griebel, M.; Schweitzer, M.A. Numerical simulation of bubble and droplet deformation by a level set approach with surface tension in three dimensions. Int. J. Numer. Methods Fluids; 2010; 62, pp. 963-993. [DOI: https://dx.doi.org/10.1002/fld.2051]

32. Nochetto, R.H.; Walker, S.W. A hybrid variational front tracking-level set mesh generator for problems exhibiting large deformations and topological changes. J. Comput. Phys.; 2010; 229, pp. 6243-6269. [DOI: https://dx.doi.org/10.1016/j.jcp.2010.04.035]

33. Scardovelli, R.; Zaleski, S. Direct numerical simulation of free-surface and interfacial flow. Annu. Rev. Fluid Mech.; 1999; 31, pp. 567-603. [DOI: https://dx.doi.org/10.1146/annurev.fluid.31.1.567]

34. Hirt, C.W.; Nichols, B.D. Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys.; 1981; 39, pp. 201-225. [DOI: https://dx.doi.org/10.1016/0021-9991(81)90145-5]

35. Davidson, M.R.; Rudman, M. Volume-of-fluid calculation of heat or mass transfer across deforming interfaces in two-fluid flow. Numer. Heat Transf. Part B Fundam.; 2002; 41, pp. 291-308. [DOI: https://dx.doi.org/10.1080/104077902753541023]

36. Ma, C.; Bothe, D. Direct numerical simulation of thermocapillary flow based on the volume of fluid method. Int. J. Multiph. Flow; 2011; 37, pp. 1045-1058. [DOI: https://dx.doi.org/10.1016/j.ijmultiphaseflow.2011.06.005]

37. Bothe, D.; Fleckenstein, S. A volume-of-fluid-based method for mass transfer processes at fluid particles. Chem. Eng. Sci.; 2013; 101, pp. 283-302. [DOI: https://dx.doi.org/10.1016/j.ces.2013.05.029]

38. Vincent, S.; De Motta, J.C.B.; Sarthou, A.; Estivalezes, J.L.; Simonin, O.; Climent, E. A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows. J. Comput. Phys.; 2014; 256, pp. 582-614. [DOI: https://dx.doi.org/10.1016/j.jcp.2013.08.023]

39. Chekifi, T. Computational study of droplet breakup in a trapped channel configuration using volume of fluid method. Flow Meas. Instrum.; 2018; 59, pp. 118-125. [DOI: https://dx.doi.org/10.1016/j.flowmeasinst.2017.11.013]

40. Mohan, A.; Tomar, G. Volume of fluid method: A brief review. J. Indian Inst. Sci.; 2024; 104, pp. 229-248. [DOI: https://dx.doi.org/10.1007/s41745-024-00424-w]

41. Brackbill, J.U.; Kothe, D.B.; Zemach, C. A continuum method for modeling surface tension. J. Comput. Phys.; 1992; 100, pp. 335-354. [DOI: https://dx.doi.org/10.1016/0021-9991(92)90240-Y]

42. Cummins, S.J.; Francois, M.M.; Kothe, D.B. Estimating curvature from volume fractions. Comput. Struct.; 2005; 83, pp. 425-434. [DOI: https://dx.doi.org/10.1016/j.compstruc.2004.08.017]

43. Popinet, S. An accurate adaptive solver for surface-tension-driven interfacial flows. J. Comput. Phys.; 2009; 228, pp. 5838-5866. [DOI: https://dx.doi.org/10.1016/j.jcp.2009.04.042]

44. Garimella, R.V.; Shashkov, M.J.; Knupp, P.M. Triangular and quadrilateral surface mesh quality optimization using local parametrization. Comput. Methods Appl. Mech. Eng.; 2004; 193, pp. 913-928. [DOI: https://dx.doi.org/10.1016/j.cma.2003.08.004]

45. Tryggvason, G.; Bunner, B.; Esmaeeli, A.; Juric, D.; Al-Rawahi, N.; Tauber, W.; Han, J.; Nas, S.; Jan, Y.J. A front-tracking method for the computations of multiphase flow. J. Comput. Phys.; 2001; 169, pp. 708-759. [DOI: https://dx.doi.org/10.1006/jcph.2001.6726]

46. Chao, Y.; Shum, H.C. Emerging aqueous two-phase systems: From fundamentals of interfaces to biomedical applications. Chem. Soc. Rev.; 2020; 49, pp. 114-142. [DOI: https://dx.doi.org/10.1039/C9CS00466A]

47. van Steijn, V.; Kleijn, C.R.; Kreutzer, M.T. Flows around confined bubbles and their importance in triggering pinch-off. Phys. Rev. Lett.; 2009; 103, 214501. [DOI: https://dx.doi.org/10.1103/PhysRevLett.103.214501]

48. Mehraji, S.; Saadatmand, M. Flow regime mapping for a two-phase system of aqueous alginate and water droplets in T-junction geometry. Phys. Fluids; 2021; 33, 072009. [DOI: https://dx.doi.org/10.1063/5.0051789]

49. Montanero, J.M.; Gañán-Calvo, A.M. Dripping, jetting and tip streaming. Rep. Prog. Phys.; 2020; 83, 097001. [DOI: https://dx.doi.org/10.1088/1361-6633/aba482]

50. Nunes, J.; Tsai, S.; Wan, J.; Stone, H.A. Dripping and jetting in microfluidic multiphase flows applied to particle and fibre synthesis. J. Phys. D Appl. Phys.; 2013; 46, 114002. [DOI: https://dx.doi.org/10.1088/0022-3727/46/11/114002]

51. Moon, B.U.; Abbasi, N.; Jones, S.G.; Hwang, D.K.; Tsai, S.S. Water-in-water droplets by passive microfluidic flow focusing. Anal. Chem.; 2016; 88, pp. 3982-3989. [DOI: https://dx.doi.org/10.1021/acs.analchem.6b00225]

52. Seemann, R.; Brinkmann, M.; Pfohl, T.; Herminghaus, S. Droplet based microfluidics. Rep. Prog. Phys.; 2011; 75, 016601. [DOI: https://dx.doi.org/10.1088/0034-4885/75/1/016601] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22790308]

53. Lai, D.; Frampton, J.P.; Sriram, H.; Takayama, S. Rounded multi-level microchannels with orifices made in one exposure enable aqueous two-phase system droplet microfluidics. Lab A Chip; 2011; 11, pp. 3551-3554. [DOI: https://dx.doi.org/10.1039/c1lc20560a] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21892481]

54. Gupta, A.; Kumar, R. Flow regime transition at high capillary numbers in a microfluidic T-junction: Viscosity contrast and geometry effect. Phys. Fluids; 2010; 22, 122001. [DOI: https://dx.doi.org/10.1063/1.3523483]