CFD study of droplet formation in a cross-Junction microfluidic device: investigating the impact of outflow channel design and viscosity ratio

Abstract

This study employed computational fluid dynamics (CFD) to investigate various regimes of droplet formation in a square microfluidic cross-junction device. The first part of the study focused on investigating the impact of viscosity ratios ranging from 0.1 to 5 on droplet characteristics under different flow rate ratios. The transition regime between the dripping and jetting exhibited particularly interesting behavior, with droplet volume being proportional to viscosity ratio. However, when conditions were far from the transition regime, the volume of droplets did not vary significantly with respect to viscosity ratios greater than one. The results of the study indicate the feasibility of reducing droplet volume to the appropriate level while minimizing the impact on the continuous phase's flow rate and the risk of transitioning to jetting. Moreover, in the final part of the study, we investigated the potential of shifting the transition boundary between the squeezing and jetting regimes. Our findings indicate that modifying the outflow channel's design can increase the frequency of droplet generation without altering the operating conditions for both the continuous and dispersed phases.

Keywords:

• CFD
• droplet microfluidics
• viscosity ratio
• frequency
• instability

1. Introduction

When miniaturization became a requirement for bio-system research decades ago, microfluidics was proposed as a method for generating micro-scale systems that facilitate highly precise physical and chemical analysis (Mark et al., 2010). For instance, when cells or magnetic beads are kept in a droplet or protein crystallization happens in nano-liter droplets, particles are often encased. This happens when droplets merge, split, or form Zheng et al. (2003). Inertial microfluidics is effective for circulating tumor cell (CTC) enrichment using platforms with controllable inlet/outlet conditions (Mashhadian & Shamloo, 2019). Several numerical models have been proposed to optimize and evaluate the performance of microfluidic platforms, facilitating the development of more efficient systems (Shi, 2023).

Fluids behave differently in microchannels than in macrochannels, due to dominant effects of wall, friction, electro-kinetic forces and surface tension. These effects, which are usually negligible in macrochannels, now play a significant role (Chen & Ren, 2017; Wu et al., 2008).

Droplet generation occurs when viscous and surface forces dominate over inertial forces at low Reynolds numbers (Squires & Quake, 2005). Droplet generation and manipulation have diverse applications, such as drug delivery (Zhao, 2013), tissue engineering (Agarwal et al., 2015) and Polymerase chain reaction(PCR) (Wang & Burns, 2009). However, high-throughput droplet generation presents challenges, including polydispersity of droplet volume and robustness in relation to droplet spacing. The study by Glawdel et al. (2011) revealed that even small variations in droplet spacing or initial conditions could result in issues with droplet sorting downstream.

Droplet formation in microfluidic systems can be achieved through various techniques that involve mixing two immiscible fluids. These techniques can be classified into two broad categories: passive and active. Active techniques involve the use of external fields, such as electrical or thermo-mechanical fields (Zhu & Wang, 2017) to create local actuation and form droplets. In contrast, passive methods modify hydrodynamical factors, such as the flow rates of the continuous and dispersed phases (Tan et al., 2006), physio-chemical properties and geometrical attributes of the device (Sugiura et al., 2002). Droplets are typically generated in configurations such as T-junctions, cross-junctions, and flow-focusing devices. Christopher and Anna (2007). Passive methods are simple, low-cost, and easy to fabricate, while active methods offer greater control but can be complex and expensive.

The viscosity of the fluids, interfacial tension, and inlet/outlet boundary conditions are among the factors that affect droplet formation, leading to non-linear behavior at the interface of two fluids. This behavior can result in the formation of droplets in various regimes, including squeezing, dripping, jetting, tubing, threading, and viscous displacement (Cubaud & Mason, 2008). Experimental and numerical studies, however, often focus on the first three regimes.

At low flow rates of the continuous and dispersed phases, droplets are formed in the squeezing regime, where the length of the dispersed droplets is equal to or greater than the channel width. In the dripping regime, the length of the droplets is reduced and is less than the width of the exit channel. As the flow rates increase, the likelihood of the jetting regime occurring increases. Within this regime, the formation of droplets becomes challenging as a cylindrical jet of dispersed phase enters the exit channel, and the tail trailing the developing droplet may grow longer. Capillary instability is observed in threads generated from all immiscible liquid pairings in the jetting regime. Ultimately, the expansion of asymmetric sinusoidal deformations causes the thread to disintegrate into multiple droplets (Eggers, 1997).

To address the complexities of droplet interface behavior in various regimes, researchers have introduced dimensionless numbers such as the Capillary number($Ca=\frac{\mu v}{\sigma }$) and flow rate ratio($Q_r=\frac{Q_d}{Q_c}$) in the literature. Here, µ, v, $Q_d$, $Q_c$, w, and σ represent the dynamic viscosity, average velocity, flow rate of the dispersed phase, flow rate of the continuous phase, characteristic length of the device, and interfacial tension, respectively. The geometry or configuration of the device can also play a critical role in droplet generation, as demonstrated by previous studies. Chekifi (2018) numerically investigated the behaviour of droplet generation in each regime of droplet formation based on configuration of the device and physical parameters such as flow rate and interfacial tension on the process.

By investigating the impact of dispersed phase inertia and viscosity ratio, recent studies aim to establish a more comprehensive guideline for predicting droplet breakup. In a recent numerical study conducted by Nekouei and Vanapalli (2017), the effect of viscosity ratio on droplet volume in capillaries of dripping in a T-junction was investigated. The absence of an established connection between viscosity ratio and droplet volume in fixed capillaries resulted in a lack of discussion on how droplet volume varies with respect to viscosity ratio. Additionally, the impact of viscosity ratio on droplet generation in the transition regime between dripping and jetting was not investigated. Moreira, Campos, and Miranda (2022) conducted a comprehensive experimental study on extracting droplet volumes using a flow-focusing device to create gelatin microdroplets. Although their study was limited to qualitative categories of viscosity ratio(low and high), they were able to closely observe droplet generation in both jetting and dripping regimes.

The purpose of this study is to create a more general method of predicting droplet's characteristics. In order to achieve this, we need to consider using more dimensionless numbers. Weber number($We=\frac{\rho v^{2}w}{\sigma }$) reflects the balance between the opposing forces of inertial and interfacial forces acting on a droplet during break-up.

Also viscosity ratio of the dispersed and continuous phases(${\mu }_r=\frac{{\mu }_d}{{\mu }_c}$) is one of the effective parameters on the dynamics of droplet formation, and it can play a crucial role in the efficacy of encapsulation process Cubaud et al. (2012). We examined the use of dimensionless numbers for classifying different regimes of droplet generation and predicting droplet volume and frequency of formation across a wide range of viscosity ratios in different regimes.

Additionally, we proposed a novel approach by modifying the device's geometry to prevent jetting and increase formation frequency in the squeezing to jetting regime, which occurs when the dispersed phase has high inertia.

2. CFD model description

2.1. Initial geometry

The schematic of the initial design of the square cross-junction used in this study is shown in Figure 1. The cross-section has equal width and depth, with $w=h=100\mu \mathrm{m}$. The dispersed phase inlet channel has a width of $w_d=w=100\mu \mathrm{m}$ and a length of $L_d=120\mu \mathrm{m}$. Each of the continuous phase inlet channels has a width of $w_c=w=100\mu \mathrm{m}$ and a length of $L_c=100\mu \mathrm{m}$. In the initial design, the outflow channel has a uniform width of $w_{exit}=w=500\mu \mathrm{m}$ and a length of $L_{exit}=15w=1500\mu \mathrm{m}$. Although this geometry serves as a starting point for assessing the impact of factors such as surface tension and viscosity, it is important to note that we will also examine the effect of different levels of local expansion in the outflow channel on the flow regimes and droplet characteristics in this work.

Figure 1. Schematic view of (a) CFD model's domain and (b) Configuration of the dispersed and continuous phases' inlets and main outlet(front view).

2.2. Governing equations

2.2.1. Hydrodynamics equations

In our approach, a set of equations governing the conservation of mass and momentum are solved for the two immiscible fluids, as follows:

Conservation of mass: (1) $\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\mathrm{\nabla }.\left(\rho \mathbf{u}\right)=0$(1) Conservative form of momentum equation: (2) $\begin{array}{rl}& \frac{\mathrm{\partial }\mathbf{\left(}\rho \mathbf{u}\mathbf{\right)}}{\mathrm{\partial }t}+(\mathbf{u}\cdot \mathrm{\nabla })\mathbf{\left(}\rho \mathbf{u}\mathbf{\right)}\\ & =-\mathrm{\nabla }p+\mathrm{\nabla }\cdot \mu (\mathrm{\nabla }\mathbf{u}+\mathrm{\nabla }{\mathbf{u}}^T)+{\mathbf{f}}_{\mathbf{s}\mathbf{t}}+{\mathbf{f}}_{\mathbf{e}\mathbf{x}\mathbf{t}}\end{array}$(2) where $\mathbf{u}$, ρ, µ, p, ${\mathbf{f}}_{\mathbf{s}\mathbf{t}}$ and ${\mathbf{f}}_{\mathbf{e}\mathbf{x}\mathbf{t}}$ denote velocity vector, density, dynamic viscosity, pressure, surface tension force and external forces respectively. Before proceeding with the remaining equations, certain assumptions must be established. Both fluids are considered to be incompressible. Additionally, due to the small scales at which we are working, it is assumed that the effects caused by gravity are negligible.

2.2.2. Phase field method

The phase field method utilizes a diffuse-interface approach to represent a finite mixing energy stored at the interface. To define the boundaries of the interface, a phase field variable is introduced, such that the concentration of the components is represented as $\frac{1+\varphi }{2}$ and $\frac{1-\varphi }{2}$, respectively. The phase field variable is defined such that $\varphi =-1$ for the dispersed phase and $\varphi =1$ for the continuous phase.

The Ginzburg-Landau-Wilson free energy is the fundamental form for representing the free energy density function as a function of the phase field variable and its gradient: (3) $f_m(\varphi ,\mathrm{\nabla }\varphi )=\frac{\lambda }{2}{\left|\mathrm{\nabla }\varphi \right|}^2+f_0\left(\varphi \right)$(3) Here, λ denotes the free energy density, and the bulk energy function $f_0\left(\varphi \right)$ is defined as a double-well potential function: (4) $f_0\left(\varphi \right)=\frac{\lambda }{4{ϵ}^2}({\varphi }^2-1{)}^2$(4) The bulk energy attempts to divide the interface into two domains and make it as sharp as possible, whilst the former part of the free energy density function represents gradient-free energy resulting from interfacial interactions seeking to make interface more diffuse and maximize mixing. The equilibrium profile of the interface is determined by the trade-off between these two components.

To determine the surface tension as a function of the field variable, we first consider a one-dimensional interface located at h = 0, and relate the mixing energy obtained from the free-energy density function to the surface energy: (5) $\gamma =\lambda {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\{\frac{1}{2}{\left(\frac{\mathrm{d}\varphi }{\mathrm{d}z}\right)}^2+f_0\}\mathrm{d}z$(5) To obtain the field variable, we have to consider a mixing equilibrium by setting $G=\frac{\delta F}{\delta \varphi }=0$. Using Equation (3(3) $f_m(\varphi ,\mathrm{\nabla }\varphi )=\frac{\lambda }{2}{\left|\mathrm{\nabla }\varphi \right|}^2+f_0\left(\varphi \right)$(3) ), we have: (6) $\frac{\delta F}{\delta \varphi }=\lambda \{-\frac{{\mathrm{d}}^2\varphi }{\mathrm{d}{z}^2}+f_0^{\prime }(\varphi \left)\right\}=0$(6) Considering the boundary conditions $f_0(\pm \mathrm{\infty })=0$, $\frac{\mathrm{d}\varphi }{\mathrm{d}z}=0$ at $z=\pm \mathrm{\infty }$, and $\varphi \left(0\right)=0$, we obtain: (7) $\frac{1}{2}{\left(\frac{\mathrm{d}\varphi }{\mathrm{d}z}\right)}^2=f_0\left(\varphi \right)$(7) By solving the above equation, we obtain the equilibrium profile of the interface: (8) $\varphi \left(z\right)=tanh\left(\frac{z}{\sqrt{2}\delta }\right)$(8) Where δ denotes the interface thickness. Substituting the preceding equation into Equation (5(5) $\gamma =\lambda {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\{\frac{1}{2}{\left(\frac{\mathrm{d}\varphi }{\mathrm{d}z}\right)}^2+f_0\}\mathrm{d}z$(5) ) yields the following relationship between interfacial tension, capillary width of the interface, and free energy density parameter: (9) $\gamma =\frac{2\sqrt{2}}{3}\frac{\lambda }{\delta }$(9) The equation above shows that if the interface width approaches the sharp-interface limit, the free energy density parameter λ should also decrease. Liu and Shen (2003) demonstrated that a diffuse-interface with an extremely narrow width would converge to a classical Navier-Stokes system with a sharp interface.

2.2.3. Dynamic evolution of interface

For the interface to behave naturally, the free energy of the system must decrease over time. This can be expressed mathematically as follows: (10) $\frac{\mathrm{D}F}{\mathrm{D}t}\le 0$(10) By substituting the free-energy density from Equation (3(3) $f_m(\varphi ,\mathrm{\nabla }\varphi )=\frac{\lambda }{2}{\left|\mathrm{\nabla }\varphi \right|}^2+f_0\left(\varphi \right)$(3) ) into the above equation and simplifying, we have: (11) $G\frac{\mathrm{D}\varphi }{\mathrm{D}t}\le 0$(11)

2.2.4. Cahn-Hillard equation

Based on the criteria given in Equation (11(11) $G\frac{\mathrm{D}\varphi }{\mathrm{D}t}\le 0$(11) ), Cahn-Hilliard equation was suggested as an evolution Equation (Cahn & Hilliard, 1959) that incorporates Fick's Law, which states that mass diffusion is proportional to the gradient of the chemical potential across the interface. The Cahn-Hilliard equation satisfies phase and mass conservation by definition: (12) $\frac{\mathrm{D}\varphi }{\mathrm{D}t}=\mathbf{-}\mathbf{\nabla }\mathbf{.}\mathit{J}$(12) Here, $\mathbf{J}$ is the mass flux, and to satisfy the criterion of Equation (11(11) $G\frac{\mathrm{D}\varphi }{\mathrm{D}t}\le 0$(11) ), the flux can be considered as $\mathbf{J}=-M\mathrm{\nabla }G$, where M is the mobility of the order parameter. By substituting the chemical potential, we arrive at the following transport equation: (13) $\begin{array}{rl}& \frac{\mathrm{\partial }\varphi }{\mathrm{\partial }t}+\mathit{u}\mathbf{.}\mathbf{\nabla }\mathit{\varphi }=\tau {\mathrm{\nabla }}^2[-{\mathrm{\nabla }}^2\varphi +\frac{\varphi ({\varphi }^2-1)}{{\delta }^2}]\end{array}$(13) (14) $\begin{array}{rl}& n.M\mathrm{\nabla }G=0\end{array}$(14) $\tau =M\times \lambda$ determines the interface's relaxation time; if it is too large, it can excessively dampen the flow. Too little diffusion will result in a drastic change in interface thickness (Jacqmin, 1999). Equation (14(14) $\begin{array}{rl}& n.M\mathrm{\nabla }G=0\end{array}$(14) ) represents the boundary condition along the walls, where n denotes the normal vector to the surface of the wall.

The last step to solve Equations (1(1) $\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\mathrm{\nabla }.\left(\rho \mathbf{u}\right)=0$(1) ) & (2(2) $\begin{array}{rl}& \frac{\mathrm{\partial }\mathbf{\left(}\rho \mathbf{u}\mathbf{\right)}}{\mathrm{\partial }t}+(\mathbf{u}\cdot \mathrm{\nabla })\mathbf{\left(}\rho \mathbf{u}\mathbf{\right)}\\ & =-\mathrm{\nabla }p+\mathrm{\nabla }\cdot \mu (\mathrm{\nabla }\mathbf{u}+\mathrm{\nabla }{\mathbf{u}}^T)+{\mathbf{f}}_{\mathbf{s}\mathbf{t}}+{\mathbf{f}}_{\mathbf{e}\mathbf{x}\mathbf{t}}\end{array}$(2) ) requires calculating the viscosity and density of the interface, which must be solved in conjunction with the phase field model equations (i.e. Equations (13(13) $\begin{array}{rl}& \frac{\mathrm{\partial }\varphi }{\mathrm{\partial }t}+\mathit{u}\mathbf{.}\mathbf{\nabla }\mathit{\varphi }=\tau {\mathrm{\nabla }}^2[-{\mathrm{\nabla }}^2\varphi +\frac{\varphi ({\varphi }^2-1)}{{\delta }^2}]\end{array}$(13) ) & (14(14) $\begin{array}{rl}& n.M\mathrm{\nabla }G=0\end{array}$(14) )). There are a number of suggestions that have been made in recent studies to calculate the density and viscosity in the interface of the two fluids. Volume average and Heaviside average are the more popular ones. In this study, we used volume average to calculate the density and viscosity in the interface, as mentioned in the following equations: (15) $\begin{array}{rl}& \rho ={\rho }_1+({\rho }_2-{\rho }_1)\varphi \end{array}$(15) (16) $\begin{array}{rl}& \mu ={\mu }_1+({\mu }_2-{\mu }_1)\varphi \end{array}$(16) Where ρ denotes density, ϕ denotes phase field parameter and µ denotes viscosity.

3. Results and discussion

3.1. Validation

Experiments from Chen and Ren (2017) are utilized to validate our model by demonstrating its accuracy and utility. To recreate the experiment in two dimensions, certain assumptions were made.

The primary assumption was the insignificance of the velocity and its gradients in the z-direction with respect to x, y, and z in the governing hydrodynamic equations. Additionally, the shallow channel depth causes the droplet to be squeezed between the lateral walls(see Figure 2), resulting in zero curvature in the plane of the lateral walls. So according to the definition of interfacial tension($F_s=\kappa \gamma$) this results in a zero interfacial tension in the z-direction. So just solving for hydrodynamic and phase field equations in the 2D would seem to be plausible. For a further comparison of 2D and 3D simulations of droplet formation and breakup, refer to Bedram et al. (2015), where the authors numerically compared the results and reported a maximum difference of 7 percent between 2D and 3D simulations.

Figure 2. (a) Front plane and (b) Side plane view of the dispersed droplet in the exit channel.

To precisely recreate the boundary and initial conditions, we have considered fully developed laminar flow at the inputs and zero pressure at the outflow. Although the exit in experiments is not immediately exposed to atmospheric pressure, this assumption is not far from the truth due to the short length of the microchannel remaining.

In the experiment, microchannels were fabricated using polydimethylsiloxane (PDMS), which is a hydrophobic material. On top of that, channels were initially filled with silicone oil, and since silicone oil swells PDMS, the interior surface of the channels is coated with oil and won't have direct contact with water. So, the wall contact angle between phases can be set to a number very close to zero (relative to the oil phase and 180 degrees relative to the water phase). Any number between zero and 8 degrees would show similar behavior of wall contact angle to the experiment according to our trial and error.

The next step is concerned with determining the parameters involved in the phase field method. Those are the mobility tuning parameter (χ) and interface thickness (δ). Interface thickness (δ) is assumed to be half the mesh element size. It is recommended by the literature (Jafari & Okutucu-Özyurt, 2016) to use this default setting and instead of changing the interface thickness manually, change the mesh size until similar results as the experiment are reached.

As it was discussed in the Section 2.2.3, the mobility parameter relates to relaxation time of interface which has a physical meaning, however its magnitude which is $M=\chi {\delta }^2$ cannot be perceived from experimental data and the magnitude of χ must be searched in conventional values used in the literature which is between 0.01 and 10 (Anandan et al., 2015; Ma et al., 2021). After performing multiple simulations with different χ values, it was observed that for $\chi <0.3$, the interface would dissipate and collapse, and the results wouldn't converge. On the other hand, for $\chi >0.3$, the simulations converged, and the results were compared for different values of χ(see Figure 3). To obtain the most similar results to the experiment, we chose $\chi =1$.

Figure 3. Frame by frame Comparison of (a) numerical results from our model and (b) experimental results obtained from Chen and Ren (2017) (Reprinted from Chen and Ren (2017),Copyright (2023), with permission from Elsevier), (c) The impact of varying χ values on droplet length.

The continuous phase in the experiment was a water solution with Glycerol. The dispersed phase was silicon oil, and the surface tension between the two phases was reported to be $37\mathrm{m}\mathrm{N}/\mathrm{m}$.

To adequately evaluate the numerical findings of our model with the experiment of Chen and Ren (2017), we can adopt their definition of a dimensionless parameter ($V_d^{\ast }$) for characterizing droplet volume under different operating conditions, such as viscosity ratio and flow rate ratios. This parameter is calculated using: (17) $V_d^{\ast }=\frac{V_d}{w^2\times h}$(17) Here, $V_d$, w, and h denote droplet volume, channel width, and channel depth, respectively. Our numerical model was run for various operating conditions, and the time spent between two droplet generation ($\mathrm{\Delta }t$) was calculated. Using the general relation between flow rate, time, and droplet volume ($V_d=\mathrm{\Delta }t\times Q_d$), we were able to determine the volume of droplets generated for each flow rate ratio. Using Equation (17(17) $V_d^{\ast }=\frac{V_d}{w^2\times h}$(17) ), we calculated $V_d^{\ast }$ to compare with the experimental results obtained by Chen and Ren (2017). In Figure 3, a frame-by-frame comparison between numerical results from simulation and experimental results are shown. For a more detailed comparison see Table 1.

Table 1. Physical properties of dispersed and continuous phases used in the study of Chen and Ren (2017) and our study.

Fluid(Glycerol-%wt conc.)	Dynamic viscosity, µ(mPa.s)	Density, $\rho \left(\mathrm{k}\mathrm{g}/{\mathrm{m}}^3\right)$	Interfacial tension with Silicon oil,γ(mN/m)
10% Glycerol	1.37	1024.3	37.1
80% Glycerol	44.1	1186.6	34.8
Silicon oil	0.93	930	−

Table 2. Validation of our numerical model with experimental results of Chen and Ren (2017) (Reprinted from Chen and Ren (2017), Copyright (2023), with permission from Elsevier) in different flow rate ratios.

Flow rate ratio($Q_r=Q_d/Q_c$)	$V_{d,simulation}^{\ast }$	$V_{d,experiment}^{\ast }$	Δt(s)	Error(%)
0.34	1.76	1.85	2.19E-02	4.67
0.45	1.83	1.9	1.71E-02	3.90
0.54	1.86	1.96	1.45E-02	5.01
0.59	1.96	2.05	1.40E-02	3.03
0.65	2.00	2.10	1.30E-02	4.52

3.2. Meshing

A quadrilateral mesh was used to mesh the domain, which was then partitioned into two main components as shown in Figure 4(b). Domain B, centered at the cross-junction, required a uniform mesh size that matched the interface thickness, as mentioned in the previous section. After multiple simulations, we determined that setting the mesh size in domain B to 1.77 µm provided an accurate capture of the interface. Furthermore, when the number of mesh elements in domain A was increased while the number of elements in domain B remained constant, we achieved mesh grid independence and closely matched the experimental data (refer to Figure 4(a)).

Figure 4. (a) Plot showing the the mesh grid independence of results for the case of $Q_r=0.34$, (b) Configuration of the two components of the mesh.

3.3. Mathematical model

In our study, the microfluidics cross-junction, as depicted in Figure 1(b), consists of two flows from separate inflow channels intersecting at a cross-shaped junction, with characteristic dimensions such as the width of the inflow channels and the exit channel, depicted in Figure 1(a).

The droplet-forming liquid enters the junction with a flow rate of $Q_d$, and the cross-flow of the second, continuous phase flowing at the rate of $Q_c$, deforms the expanding interface, either forcing it to break up before it reaches the exit channel(the dripping regime) or necking it at the junction(the squeezing regime) with the dispersed droplet penetrating inside the exit channel.

Further deformation of the interface after entering the exit channel may result from pressure drops across the developing droplet or shear stress at the interface's neck. A schematic showing influential forces is shown in Figure 5.

Figure 5. (a) Schematic of forces acting on droplet and curvature of droplet profile at the neck, (b) Important dimensions of forming droplet.

Research on the analytical prediction of droplet volume (Christopher & Anna, 2007; Garstecki et al., 2006; Husny & Cooper-White, 2006; Jena et al., 2021; Thorsen et al., 2001; Xu & Nakajima, 2004) agrees on three forces–the interfacial force, the viscous force, and the pressure drop force–as being primary contributors to droplet detachment. In the following sections, we will define the essential dimensions for characterizing the forces acting on a droplet and derive each force. Then, we will investigate numerically the influence of viscosity ratio on each portion of the analytical relation predicted for droplet volume in the squeezing and dripping regimes.

3.4. Effects of viscosity ratio on droplet volume and formation frequency

3.4.1. I: squeezing regime

In the squeezing regime, droplets extend as far as $L_{droplet}$ in the exit channel with a width of $w_{exit}$, completely compressed between the out-of-plane walls. $w_g$ indicates the width of the gap between walls and droplets where the continuous phase is flowing. Using these parameters and the relations originally provided by Garstecki et al. (2006), we can calculate the forces acting on a forming droplet.

The difference in Laplace pressure between the droplet's tip and neck provides the capillary force, which opposes droplet breakup. As shown in Figure 5(a), the droplet's tip has a curvature of $\frac{2}{w_{exit}}$ in the cross-channel direction and a curvature of $\frac{2}{h}$ in the out-of-plane direction. The neck of the droplet has a curvature of $\frac{2}{h}$ in the depth direction and an in-plane curvature of $\frac{1}{R_{neck}}$. Thus, we can write the capillary force as follows: (18) $F_{\sigma }\approx \sigma [(\frac{2}{h}+\frac{1}{w_{exit}})-(\frac{2}{h}+\frac{2}{R_{neck}})]A$(18) Here, σ, h, $w_{exit}$, and A denote interfacial tension, exit channel's depth and width, and projected area of the droplet in the channel direction, respectively. If we disregard the rounded corners of the droplet, the projected area of a droplet, as seen in the Figure 2(b), is essentially a rectangle with an area of $A=w_{exit}h$. Furthermore, to simplify the equation, Figure 5(a) reveals that $\frac{R_{neck}}{w_{exit}}\sim O\left({10}^0\right)$, thus $R_{neck}$ may be substituted with $w_{exit}$ in the equation. Therefore, the equation can be written as: (19) $F_{\sigma }\approx \sigma \left(\frac{-1}{w_{exit}}\right)w_{exit}h\approx -\sigma h$(19) The continuous phase that penetrates through the gap between the droplet and channel walls induces a shear viscous stress on the droplet. This stress can be approximated as $\tau =\frac{{\mu }_c{u}_g}{\sigma }$, where ${\mu }_c$, $u_g$, and σ denote the dynamic viscosity of the continuous phase, the average velocity of the continuous phase in the gap, and the interfacial tension, respectively. If we rewrite the equation to replace $u_g$ with $\frac{Q_c}{w_{g}h}$, we obtain the following equation: (20) $F_{\tau }\approx \frac{{\mu }_c{Q}_c}{w_g^2}A\approx \frac{{\mu }_c{Q}_c{L}_{droplet}{w}_{exit}}{w_g^2}$(20) The last force is the squeezing force, which is created by the squeezing pressure pressing on the droplet due to channel blockage. For an approximation of this pressure loss, a complete derivation can be found in previous studies (Garstecki et al., 2006; Stone, 2005), which states that: (21) $F_{\mathrm{\Delta }p}\approx \mathrm{\Delta }{p}_{c}A\approx \frac{{\mu }_c{Q}_c}{w_g^2}\frac{L_{droplet}}{w_g}\approx \frac{{\mu }_c{Q}_c{L}_{droplet}{w}_{exit}^2}{w_g^3}$(21) The net force acting on a droplet as a control mass in the direction of the exit channel, as stated by Newton's law, is equal to its momentum change in the same direction. Therefore, we can write: (22) $\mathrm{\Sigma }F=F_{\mathrm{\Delta }p}+F_{\tau }+F_{\sigma }=\frac{\mathrm{d}\left(mv\right)}{\mathrm{d}t}=u_{cm}\frac{\mathrm{d}m}{\mathrm{d}t}+m\frac{\mathrm{d}v}{\mathrm{d}t}$(22) m and $u_{cm}$ denote the mass and velocity of the control mass, respectively. Observing the process of droplet generation in experimental studies (Chen & Ren, 2017), we can conclude that the change in velocity of the forming droplet as a control mass is negligible. Therefore, we can assume that $\frac{\mathrm{d}v}{\mathrm{d}t}\approx 0$. According to the assumptions stated by Jena et al. (2021), the shallow depth of the channel will provide minimal velocity gradients, and the velocity of the control mass, $u_{cm}$, will be affected by the velocity of the larger velocity, whether $u_d$ or $u_c$. However, a plausible assumption is that $u_{cm}=\frac{u_c+u_d}{2}$.

Since every inlet and outlet in our geometry has a width of w, for a more generic form, we can substitute $l_{droplet}=\frac{L_{droplet}}{w}$ for the non-dimensional length of the droplet and $Q_r=\frac{u_d}{u_c}$ for the flow rate ratio. By substituting the above terms into Equation (22(22) $\mathrm{\Sigma }F=F_{\mathrm{\Delta }p}+F_{\tau }+F_{\sigma }=\frac{\mathrm{d}\left(mv\right)}{\mathrm{d}t}=u_{cm}\frac{\mathrm{d}m}{\mathrm{d}t}+m\frac{\mathrm{d}v}{\mathrm{d}t}$(22) ), we obtain: (23) $\begin{array}{rl}& {\mu }_c{u}_c{L}_{droplet}\frac{w_{exit}}{w_g^2}h+{\mu }_c{u}_c{L}_{droplet}\frac{w_{exit}^2}{w_g^3}h-\sigma h\\ & =u_{cm}\frac{\mathrm{d}m}{\mathrm{d}t}=\left(\frac{u_d+u_c}{2}\right){\rho }_d{u}_{d}wh\end{array}$(23) (24) $\begin{array}{rl}& \frac{{\mu }_c{u}_c{l}_{droplet}}{\sigma }({\left(\frac{w}{w_g}\right)}^2+{\left(\frac{w}{w_g}\right)}^3)-1\\ & =\left(\frac{Q_r+1}{2Q_r}\right)\frac{{\rho }_d{u}_d^2{w}_d}{\sigma }\end{array}$(24) In the non-dimensional form, the Capillary number of the continuous phase, $C{a}_c=\frac{{\mu }_c{u}_c}{\sigma }$, and the Weber number of the dispersed phase, $W{e}_d=\frac{{\rho }_d{u}_d^{2}w}{\sigma }$, are clearly present. These numbers indicate how much each force contributes to droplet breakup.

For instance, if the Weber number of the dispersed phase is relatively large, the inertia of the dispersed phase will have a significant effect on droplet volume. However, in the case of a squeezing regime with low inertia and the droplet being squeezed between the walls, Equation (25(25) $l_{droplet}=\frac{1+\left(\frac{Q_r+1}{2Q_r}\right)W{e}_d}{C{a}_c\left(\right(\frac{w}{w_g}{)}^2+\left(\frac{w}{w_g}{)}^3\right)}$(25) ) demonstrates that $\frac{F_{\mathrm{\Delta }p}}{F_{\tau }}=\frac{w}{w_g}\gg 1$, suggesting that the pressure drop force has a much stronger influence on droplet volume than the viscous force and inertia of the dispersed phase.

Noting that the viscosity ratio does not appear in Equation (24(24) $\begin{array}{rl}& \frac{{\mu }_c{u}_c{l}_{droplet}}{\sigma }({\left(\frac{w}{w_g}\right)}^2+{\left(\frac{w}{w_g}\right)}^3)-1\\ & =\left(\frac{Q_r+1}{2Q_r}\right)\frac{{\rho }_d{u}_d^2{w}_d}{\sigma }\end{array}$(24) ) but has been shown to have an effect on droplet breakup, it is reasonable to assume that the viscosity ratio impacts the gap width. Therefore, a useful experiment would include measuring the relationship between viscosity ratio and gap width under different operating conditions. Now that we have investigated the effects of every important parameter influencing droplet breakup, we can predict each parameter's effect on droplet volume based on the equation: (25) $l_{droplet}=\frac{1+\left(\frac{Q_r+1}{2Q_r}\right)W{e}_d}{C{a}_c\left(\right(\frac{w}{w_g}{)}^2+\left(\frac{w}{w_g}{)}^3\right)}$(25) And lastly, to derive the volume of the droplet, based on our assumption of modeling the droplet as a spherocylinder, the droplet's volume can be computed as: $V_d=L_{droplet}{w}^{2}h+\frac{\pi w^{2}h}{4}$.

In this study, we used the physical properties of the continuous phase as reported by Chen and Ren (2017) and the properties of the dispersed phase as reported by Jena et al. (2021) for viscosity. Numerical simulations were conducted to explore the impact of viscosity ratio and flow rate ratio on droplet volume in the squeezing regime. To ensure that the system remained in the squeezing regime, we maintained the Capillary number of the continuous phase close to $C{a}_c=0.001$. The viscosity ratio was varied between 0.01 and 5, and the flow rate ratio was varied between 0.25 and 1.375. The results indicate that the gap width does not appear to be affected by the viscosity ratio when the flow rate ratio is varied from 0.25 to 1.375. Instead, the gap width appears to depend solely on the flow rate ratio(see Figure 6(a)). This finding suggests that the droplet volume in the squeezing regime is predominantly determined by the flow rate ratio, as previously reported by Garstecki et al. (2006). Therefore, by adjusting the flow rate ratio, we can tune the droplet volume to the desired value in the squeezing regime.

Figure 6. Resulted graphs for (a) Gap width in different flow rate ratios and viscosity ratio, (b) Droplet volume for different flow rate ratios and viscosity ratios.

Furthermore, as the flow rate ratio increases, we approach the transition regime from squeezing to dripping or jetting, and our numerical results deviate from the theoretical prediction (see Figure 6(b)). This deviation occurs because the reduction in droplet length makes the spherocylindrical droplet model less accurate, and the simplifications used in the derivation of forces are no longer valid. Therefore, we should exercise caution when interpreting our results near the transition regime and consider using more accurate models and methods.

3.4.2. II: dripping and the transition regimes

The dripping regime is characterized by the formation of small droplets and occurs when the capillary number of the continuous phase surpasses a critical threshold, which previous research has shown to be $C{a}_c=0.01$ (Christopher & Anna, 2007; Cubaud & Mason, 2008).

In the dripping regime, we expect the Capillary number to play a significant role, as high Capillary numbers can result in significant inertia in the droplet's neck. The inertia of the continuous phase can cause the droplet to break up before entering more than w(width of channel) in the exit channel, or in the case of high dispersed phase flow rates, cause the velocity of the control mass(droplet) to significantly accelerate, resulting in the droplet formation process entering the jetting regime. Therefore, controlling the Capillary number and the flow rate ratio is crucial in determining the droplet size and preventing the formation of multiple small droplets in the dripping regime.

The transition from the squeezing regime to the dripping regime is not sudden. As we increase the Capillary number from 0.001 to 0.01, the features of the squeezing regime are gradually replaced by those of the dripping regime. In the dripping regime, viscous shear stresses are primarily responsible for droplet breakup, but in the squeezing regime, pressure forces dominate. Additionally, the gap between the droplet and the walls widens, such that $\frac{w_g}{w}\approx O\left({10}^0\right)$, which is much larger than the corresponding value in the squeezing regime.

According to Christopher and Anna (2007), droplet breakup in the dripping regime is often divided into two distinct stages: filling and necking. During the filling stage, the droplet continues to spread as a cylinder-shaped parallel flow alongside the continuous phase, as shown in Figure 7(a). As the droplet approaches the entrance of the exit channel, the viscous force at the interface increases until the total of the viscous and pressure drop forces is equal to the capillary force and inertia of the dispersed phase. This moment marks the beginning of the necking stage.

Figure 7. (a) Filling stage and (b) Necking stage in the dripping regime.

By calculating the balance of forces according to Newton's law, Christopher et al. determined the equation for the droplet's size during the filling stage: $(1-\overline{b}{)}^3=\overline{b}\times Ca$, where $\overline{b}$ is the dimensionless droplet length. Christopher et al. adopted the same assumption as Garstecki for the duration of necking, which assumes that the rate of necking of the interface is equal to the average velocity of the continuous phase in its inflow channel (i.e. $u_{neck}=u_c=\frac{Q_c}{wh}$).

Regarding this assumption, the volume that is injected into the forming droplet after the necking stage is equal to $V_{neck}=\frac{Q_{d}w}{u_{neck}}$. Eventually dimensionless volume of droplet can be found as: $\overline{V_d}={\overline{b}}^2+Q_r$, which doesn't account for viscosity ratio in necking or filling stages. Other studies also haven't been successful in deriving the accurate prediction for droplet volume, since the shape of droplet in the necking stage is complex and transforms so quickly that deriving the accurate formulation for each force is impossible. Similarly, we cannot distinguish the forces operating on the droplet in the dripping and the transition regimes as easily as we could for droplet production in the squeezing regime in our experiment. However, dimensionless numbers may still be used to evaluate the dynamics of droplet break-up in the dripping regime.

We conducted simulations by varying the viscosity ratio from 0.05 to 5 and recording the observations for three different flow rates ($Q_r=0.50,1.75$) and three different Capillary numbers of the continuous phase ($C{a}_c=0.005,0.01,0.05$). At a constant flow rate ratio, Our results show that as the Capillary number increases, the droplet volume decreases, as seen in Figure 8(a). Similarly, the volume of droplets decreases as the viscosity ratio increases, although for ${\mu }_r>1$, the decrease in droplet size with respect to an increase in viscosity ratio becomes less significant.

Figure 8. (a) Droplet's volume against viscosity ratio for $Q_r=0.25$ in different capillaries, (b) Droplet's volume versus viscosity ratio in different flow rate ratios and $C{a}_c=0.01$.

In the transition regime between dripping and jetting ($C{a}_c=0.05$), the relationship between droplet volume and viscosity ratio becomes linear, and the effect of viscosity ratio becomes increasingly significant. Another experiment examined the effect of viscosity ratio at different flow rates, and Figure 8(b) demonstrates that for high flow rates ($Q_r>1$), where the droplet generation regime falls in the transition between dripping and jetting, the slope of the droplet volume versus viscosity ratio becomes more decreasing and does not follow the power-law trend observed in regimes where no transition is occurring.

To accurately predict droplet volume, we derived formulas for two transition scenarios: transitions from squeezing to dripping and from dripping to jetting. In the transition from squeezing to dripping, we fit the data using a power-law formulation. For example, for $Q_r=0.5$ and $C{a}_c=0.01$, we fit the data for different viscosity ratios using a power-law formula for the dimensionless droplet volume as follows: (26) $V_d^{\ast }=\frac{V_d}{w}=f\left(Q_r\right)C{a}_c^m{\mu }_r^n$(26) By fitting the acquired data on a power-law curve, we can conclude that m = −0.35 and n = −0.063. Furthermore, from Figure 8(a), we can observe that the power coefficient of ${\mu }_r$ appears to be a function of $C{a}_c$ at low flow rate ratios in the transition regime between squeezing and dripping. Therefore, for a more general formula for predicting droplet volume, we apply linear proportionality to the flow rate, as used by previous studies. We can write this as: (27) $V_d^{\ast }=(A{Q}_r+B)C{a}_c^m{\mu }_r^{g\left(C{a}_c\right)}$(27) Here, $g\left(C{a}_c\right)$ denotes a function of the Capillary number of the continuous phase. However, as the flow rate ratio or Capillary number of the continuous phase is increased, we approach the second transition regime, where the proportionality between droplet volume and viscosity ratio becomes linear: (28) $V_d^{\ast }=A(Q_r,C{a}_c)+B(Q_r,C{a}_c){\mu }_r$(28) Both A and B are functions of the flow rate and Capillary number of the continuous phase. The graphs in Figure 8 illustrate the linear functions A and B for various flow rate ratios and Capillary numbers. The decreasing trend of droplet volume versus viscosity ratio suggests that, since the frequency of droplet formation is given by $f=\frac{Q_d}{V}$ for a given droplet volume, a higher frequency of droplet formation can be achieved with a larger flow rate of the dispersed phase by increasing the viscosity ratio. It is important to note that in all the experiments conducted for different Capillary numbers in the dripping and transition regimes, the product of the droplet Weber number and $\frac{Q_r+1}{2Q_r}$ was kept constant at $6.6\times {10}^{-4}$ to solely observe the effect of the Capillary number of the continuous phase. Furthermore, it is noteworthy that when the viscosity ratio was increased, there was no threading observed. This could be a great alternative to increasing the flow rate of the continuous phase to reduce droplet volume without transitioning to the threading regime.

The variation of droplet volume in a fixed Capillary and Weber number system can be an indicator of another operating parameter, namely the viscosity ratio, which can be utilized to control the droplet's frequency. For example, Figure 8(b) shows that two different operating conditions with the same droplet volume and fixed capillary number of 0.01, but different flow rates, result in different frequencies of formation. Specifically, the first case with $(Q_r,{\mu }_r)=(1,2.5)$ leads to a higher frequency of formation due to the higher flow rate, as determined by the relationship $f=\frac{Q_d}{V_d}$. The second case with $(Q_r,{\mu }_r)=(0.5,2.5)$ leads to a lower frequency of formation with the same droplet volume as the first case due to the lower flow rate.

Furthermore, this study proposes a novel approach for determining the viscosity of an unknown fluid in the dispersed phase. By exploiting the known relationship between droplet volume and viscosity ratio at different Capillary numbers, the unknown viscosity can be determined by injecting the unknown fluid as the dispersed phase and manipulating the operating conditions, such as the flow rate. The droplet volume and viscosity ratio can then be measured to solve for the unknown viscosity. This technique has significant applications in various fields, including microfluidic devices and industrial processes.

3.5. Controlling jet streams

A phenomenon commonly encountered when creating droplets in high flow rates is jetting instability, also known as Rayleigh-Plateau instability. This phenomenon results in the extension of the thread that forms between droplets as they grow, leading to longer thread length.

A mathematical examination of the term $u_{cm}\frac{\mathrm{d}m}{\mathrm{d}t}$ in Equation (25(25) $l_{droplet}=\frac{1+\left(\frac{Q_r+1}{2Q_r}\right)W{e}_d}{C{a}_c\left(\right(\frac{w}{w_g}{)}^2+\left(\frac{w}{w_g}{)}^3\right)}$(25) ) reveals that it can be decomposed into two distinct components: $u_{cm}$, which accounts for the inertia effects of both the dispersed and continuous phases, and $\frac{\mathrm{d}m}{\mathrm{d}t}$, which is solely determined by the flow rate of the dispersed phase. In the context of jetting, the velocity of the control mass ($u_{cm}$) is significantly influenced by the velocity of the continuous phase, causing the inertia of the control mass to dominate over the pressure drop force and viscous shear stress. This can delay or impede droplet breakup and suggests that droplet generation may not occur at all flow rate ratios or capillaries of the continuous phase.

If this barrier were to be overcome, it would suggest that higher flow rates, previously inaccessible for continuous droplet generation due to jetting, could now be achieved. As per the relationship $f=\frac{Q_d}{V_d}$, this would result in an increase in the frequency of droplet generation.

To overcome jetting, we conducted experiments by introducing sudden expansions along the path of droplet formation in the jetting or transition to jetting regime. Figure 9(a) illustrates the characteristic dimensions of the sudden expansion in the exit channel.

Figure 9. (a) Characteristic dimensions of a sudden expansion along the exit channel, (b) Critical points in a normal channel, (c) Critical points in exit channel with a sudden expansion.

Our initial simulation involved observing the jet stream in the transition regime between the dripping and jetting regimes, both in a normal exit channel and an exit channel with a sudden expansion, as shown in Figure 11(a-b), respectively. Despite the fact that the expansion slowed the jet tip, it continued to grow in the exit channel after passing through the expanded section depicted in Figure 11. Shear stress contours, as depicted in Figure 10, also revealed only minor differences between the two scenarios. This suggests that, despite the channel expansion, the increase in viscous shear or pressure drop forces was insufficient to overcome the combined effects of inertial and surface tension forces.

Figure 10. Shear stress contours in (a) Normal exit channel, (b) Schematics of droplet generation in a normal exit channel, (c) Shear stress contours in the vicinity of expanded part of the exit channel, (d) Schematics of droplet generation in an exit channel with sudden expansion.

Figure 11. Sequence of frames showing that a sudden expansion won't slow down the jet's tip and resolve jetting in the transition from the dripping to the jetting regime.

In the next experiment, we investigated jet streaming in the transition regime between squeezing and jetting. In the squeezing regime, the pressure force is the primary factor responsible for fragmentation. As shown in Figure 12, a sequence of frames demonstrates how the expansion in the exit channel prevented the jet's streaming tip from developing. As a result, the pressure in that region increased, contributing to the pressure drop caused by the channel blockage at the tip.

Figure 12. Schematics of (a) Sequence of frames for normal exit channel droplet generation with $Q_c=8\frac{\mu L}{min}$ and $Q_d=16\frac{\mu L}{min}$, (b) Sequence of frames for droplet generation with $Q_c=8\frac{\mu L}{min}$ and $Q_d=16\frac{\mu L}{min}$ with an expansion of $H_{\exp }=1.5w$, $W_{\exp }=w$ and $d_{\exp }=2.5w$.

Capillary pressure also appears to play a significant role in this phenomenon. The tip of a blocked stream that enters an expansion experiences capillary pressure due to the sudden change in free interfacial energy, as illustrated in Figure 13.

Figure 13. A droplet expanding into a sudden expansion along the channel, (a) Before entering, (b) After entering, (c) Configuration of interfacial tensions when three surfaces meet at the wall.

To better comprehend the physics and thermodynamics of the capillary pressure barrier, we can refer to the concept of free interfacial energy (Kim & Whitesides, 1997; Man et al., 1998). The free interfacial energy describes the total energy stored at the interface between the phases of a multiphase flow or static state, including solid-liquid-gas phases. It can be expressed as: (29) $U_T=A_{os}{\gamma }_{os}+A_{ws}{\gamma }_{ws}+A_{ow}{\gamma }_{ow}$(29) where $A_{os}$, $A_{ws}$, and $A_{ow}$ denote the interfacial surface areas between solid-liquid, solid-gas, and liquid-gas phases, respectively. ${\gamma }_{os}$, ${\gamma }_{ow}$, and ${\gamma }_{ws}$ are the corresponding surface tensions. Young's equation, which relates the surface tensions between the solid-liquid, solid-gas, and liquid-gas interfaces, can be incorporated by considering the static contact angle, ${\theta }_c$, resulting in ${\gamma }_{os}+{\gamma }_{ow}\cos \left({\theta }_c\right)={\gamma }_{ws}$. By differentiating this equation with respect to the volume of injected dispersed phase, we obtain the capillary pressure: (30) $P_{cap}=-\frac{\mathrm{d}{U}_T}{\mathrm{d}{V}_P}={\gamma }_{ow}(\cos {\theta }_c\frac{\mathrm{d}{A}_{os}}{\mathrm{d}{V}_P}-\frac{\mathrm{d}{A}_{ow}}{\mathrm{d}{V}_P})$(30) Figure 13(a-b) illustrate the expansion of a droplet. By varying θ from zero to 90 degrees, we can calculate the capillary pressure using Equation (31(31) $R\left(\theta \right)=\frac{w+2r_f(1-\cos (\theta \left)\right)}{-2\cos (\theta +{\theta }_c)}$(31) ). To accomplish this, we first need to determine the geometrical parameters $A_{os}$, $A_{ow}$, and $A_{sw}$. Additionally, dv represents the infinitesimal volume of the dispersed droplet that has entered the expanded channel. From the figure, we can write $A_{os}=r_f\theta h$, $A_{ow}=R\left(w\right)(\pi -2\text{β})h$, and $\text{β}=\pi -({\theta }_c+\theta )$ can be determined from the geometric relations. Similarly, for $R\left(\theta \right)$, we can write: (31) $R\left(\theta \right)=\frac{w+2r_f(1-\cos (\theta \left)\right)}{-2\cos (\theta +{\theta }_c)}$(31) In the study conducted by Man et al. (1998), the pressure of the fluid upon entering the channel was measured. They observed that when approximately 2 picoliters of droplet entered the channel, a negative pressure of 2 kPa was generated to prevent the fluid from growing into the expanded part. In our numerical model, we selected several critical points within the geometry to measure the pressure quantitatively. The results are presented in Figure 14.

Figure 14. Pressure fluctuation in critical points for (a) Normal exit channel, (b) Expanded exit channel.

We modified the width and position of the expansion along the exit channel and observed that a sudden expansion with a width 1.5 times that of the channel and a length equal to the channel width could prevent jetting, as shown in Figure 12.

When we investigated the process of producing jet streams with lower dispersed phase flow rates, the importance of capillary pressure became apparent. For example, at flow rates of $Q_d=8\frac{\mu L}{min}$ and $Q_c=10\frac{\mu L}{min}$, even though the process occurred theoretically in the transition regime between squeezing and jetting, no configuration of sudden expansion in the exit channel induced the droplet to break up. This observation demonstrates that even significant sudden expansions in the channel, which would increase the pressure at the droplet's tip, cannot cause the droplet to break apart.

Therefore, we propose that when a droplet experiences a sudden expansion, its velocity abruptly decreases and the tip, due to the high flow rate, begins to grow and fill the gap width. This causes the droplet to encounter a sudden expansion, similar to the one shown in Figure 9, and ultimately experience the capillary pressure.

The results presented in Figure 12(a-b) demonstrate that an expansion with a width of 2 times the channel width and a length equal to the channel width can prevent jetting and facilitate droplet generation at the expansion location. Figure 14 illustrates the pressure fluctuations measured at critical points CP1 through CP6 (see Figure 9(b-c)). The measurements indicate that creating a local expansion in the exit channel can significantly increase the pressure at the droplet thread's neck, thereby contributing to droplet detachment. In a normal exit channel, the pressure in the droplet's neck drops after approximately 15 ms, but in an exit channel with a sudden expansion, the jet's tip is slowed down, allowing pressure to build up in critical locations. This demonstrates the success of droplet generation.

Since flow rates are low in the squeezing regime, the rate of jet tip displacement is also low. Therefore, an expansion can help resolve jet instability by preventing the growth of the jet's tip. From this perspective, another criterion is to locate the optimal site for a sudden expansion.

Having successfully resolved jetting, we can now produce droplets of a size that was previously only achievable with a lower dispersed phase flow rate, in the flow rates of the transition regime between squeezing and jetting. Using the relationship $f=\frac{Q_d}{V_d}$, we can conclude that the frequency of droplet formation has increased.

3.6. Conclusion

In our previous studies, we employed both numerical and experimental methods to analyze microfluidic devices (Besanjideh et al., 2022; Ghazimirsaeed et al., 2021; Naghdloo et al., 2019; Nasiri et al., 2021; Shamloo et al., 2016, 2020). In this study, we investigated the effect of viscosity ratio on droplet volume for varying flow rate ratios using a cross-junction device as our first observation, employing numerical methods. In the squeezing regime, we observed a weak relationship between droplet volume and both the viscosity ratio and the Capillary number of the continuous phase.

As we transitioned from the squeezing to the dripping regime, it was observed that the volume of the forming droplet became more dependent on the Capillary number of the continuous phase. We found that the relationship between droplet volume and Capillary number followed a power-law with coefficients dependent on the Capillary number of the continuous phase, with a weak dependency towards the flow rate ratio.

Furthermore, we observed that the relationship between droplet volume and viscosity ratio in the transition regime between dripping and jetting undergoes a complete change in behavior. It becomes linear, with coefficients that are dependent on the flow rate ratio and Capillary number of the continuous phase.

Subsequently, when operating conditions are away from transition regimes, we observed that the relationship between droplet volume and viscosity ratio weakens for viscosity ratios greater than one. This is supported by the findings of Nekouei and Vanapalli (2017).

Based on our results, it is possible to reduce droplet volume by increasing the viscosity of the continuous phase without altering the inlet flow rates or channel geometry. Therefore, by adding shear-thinning fluids to the continuous phase or modifying the dispersed phase, it is possible to increase or decrease the viscosity ratio, which can alter the droplet volume. This enables us to regulate the frequency of droplet production without affecting the flow rates.

Subsequently, the jetting regime was investigated to identify possible methods to modify the geometry of the exit channel and prevent jetting. Our findings revealed that capillary pressure plays a crucial role in eliminating jetting. For specific flow rates where the transition from squeezing to jetting occurs, we observed that an abrupt expansion in the exit channel can halt the droplet, causing its tip to extend and fill the exit wall. This results in a sudden change in the droplet's contact area with the walls, which generates a negative pressure and helps in droplet break-up.

3.7. Future work

Although our study investigated a broad range of viscosity ratios, we had to limit ourselves to capillaries less than 0.05. This is because for capillaries greater than that, the threading regime would be observed, which requires numerical methods to accurately capture the interface at the thread. Experimental studies could complement our work by studying different capillaries in the threading regime and investigating the possibility of creating satellite droplets.

In addition to investigating the flow rates at which jet resolution is achievable using this method, future numerical or experimental studies should also explore several geometry alterations to engage capillary pressure as an effective force contributing to the breakup of the droplet before the creation of jet streams. Moreover, polydispersity of droplets after the resolution of jetting is an important factor to consider and should be investigated in future studies.

Disclosure statement

No potential conflict of interest was reported by the author(s).