DIFFERENT STAGES OF LIQUID FILM GROWTH IN A MICROCHANNEL: TWO-PHASE LATTICE BOLTZMANN STUDY

Abstract

A free energy model is used to describe the droplet formation and break-up process in a T-junction bio-microchannel. Droplets are created as a result of interaction of two immiscible liquids. Different stages for the droplet formation process are analyzed which are: a) growing in the x and y directions, b) growing in the x direction and c) detachment process. The effects of capillary number and flow rate ratio on the droplet formation stages are also studied. The influences of Capillary number, flow rate ratio, viscosity ratio and geometrical parameters on droplet break up, droplet size and detachment time are systematically studied. By increasing the flow rate ratio; the duration of droplet formation and the length of the x-growth stage are decreased for small capillary numbers. For larger capillary numbers; the droplet penetrates toward the downstream; therefore, the length of the x-growth stage is increased. The start of detachment process in the microchannel is also reported, which is related to narrowing of the neck of the liquid film. The results show that the detachment time is increased by decreasing the Capillary number. For Ca>0.02, the detachment time is independent of the flow rate ratios. Moreover; the effects of viscosity ratios on detachment time are not significant in comparison to the effects of capillary number. For Ca<0.04, the size of the droplet is independent of the viscosity ratio, but after the critical Capillary number (i.e., Ca=0.04), the size of the droplet is varied by the viscosity ratio. The time between two consecutive drops is also decreased by increasing the Capillary number. Moreover, this time is decreased by increasing the flow rate ratio until Ca=0.04. After this Capillary number, the flow rate ratios have no significant effect on the time between two consecutive droplets. An exhaustive validation study is performed including (a) the Laplace equation in the stationary droplet; (b) a contact angle test; (c) Taylor deformation test in shear flow and (d) comparison of the droplet length as a function of flow rate ratio between the present work and other studies.

Keywords
Stages of Liquid Film Growth; Microchannel; Lattice Boltzmann Method; Critical Capillary number

INTRODUCTION

The dynamics of droplet formation in a second fluid inside microchannels have important applications in drug delivery (Richter et al., 1997), food industries (Cal, 2005), mixing in microfluidic multiphase systems (Losey et al., 2002), and a wide range of applications in biology and chemistry. Monodisperse droplets can be generated via different methods in microfluidic devices, including cross-junction devices (Yasuno et al., 2004; Sugiura et al., 2004; Elzanfaly et al., 2015), T-junctions (Cheng et al., 2016; Liu et al., 2015; Sbragaglia et al., 2016; Thorsen et al., 2001; Nisisako et al., 2002; Xu et al., 2006; Garstecki et al., 2006; Van der Graaf et al., 2006; Christopher et al., 2008; Azarmanesh and Farhadi, 2015; Malekzadeh and Roohi, 2015), flow-focusing devices (Anna et al., 2003; Cubaud et al., 2005; Garstecki et al., 2005; Fu et al., 2009) and co-flowing devices (Umbanhowar et al., 2000; Hua et al., 2007; Alizadeh et al., 2015).

In recent years, both experimental and numerical studies have investigated on the mechanisms of droplet formation in microchannels (Cramer et al., 2004; Fu and Pan, 2005; Günther and Jensen, 2006; Waelchli and Rohr, 2006; Wu et al., 2008; Yu et al., 2007). Thorsen et al. (2001) performed experiments to investigate the formation of water droplets in oil. They reported that both shear forces and surface tension control the formation of droplets. Garstecki et al. (2006) analyzed two-phase flow in a T-junction microchannel experimentally. They suggested that the formation of the droplet is related to the pressure difference across the droplet neck at low Capillary numbers. On the other hand, when the Capillary number is large, the droplets are formed by capillary instability (Günther and Jensen, 2006). Liu and Zhang (2011) investigated the influences of viscosity ratio of two phases, flow rate ratio and Capillary number on water droplet generation in oil in a cross-junction. Hua et al. (2007) used the front-tracking/finite volume method to investigate the mechanism of droplet formation in a co-flowing microchannel. The effects of flow rate of the continuous phase, viscosity, and the surface tension on the size of droplets were investigated. The correlations of droplet size, Reynolds number, Weber number, Capillary number and viscosity ratio were also reported by the authors. In recent years, the Lattice Boltzmann method(LBM) has emerged as a powerful tool for many technical applications involving complex fluid dynamics, in particular, in the breakup and formation of droplets in micro devices. By using meso and microscopic Boltzmann's kinetic equation in LBM, the macroscopic fluid dynamics can be obtained (Sukop and Thorne, 2006). The objective of this study is to achieve a deeper understanding of droplet formation in a T-junction microchannel. The influence of Capillary number, flow rate ratio, viscosity ratio, contact angle and geometrical parameters on droplet break up, droplet size and the detachment process are systematically studied. It is important to note that the different stages of droplet growth in the microchannel are analyzed. In other word, the duration of droplet formation in the microchannel (related to the length of the x-growth stage along the microchannel) is investigated based on the capillary numbers.

There are several popular multiphase models in the LBM, including the color model, the Shan-Chen model, the free-energy model and the mean-field model (Liu et al., 2012; Liu et al., 2016; Chen et al., 1999; Gunstensen et al., 1991; Shan and Chen, 1993; Gouyet et al., 2003). The color model and the Shan-Chen model are widely used in two phase flows because of their easy implementation. However, the disadvantages of the color model and Chan-Chen model are large spurious velocities at the interface. Moreover; the local momentum is not conserved in these methods. The spurious velocities are also high at the interface for the interparticle potential model. The computing time is quite long for the mean-field model. Considering both the strangeness and weakness of each model, the present simulation is based on the free energy model in the Lattice Boltzmann framework. Local conservation of momentum, low spurious velocities and a controllable interface for two immiscible phases are the main advantages of the free energy model.

NUMERICAL METHOD

The present Lattice Boltzmann simulations were conducted by using the scheme developed by Liu and Zhang (2009), which is a variant of Swift et al. (1996). The existence of a free-energy functional which controls the equilibrium properties of the two-phase system is the important feature of the phase-field models. The dynamics of a multiphase system, including local momentum conservation, low spurious velocities and thin interface of two immiscible fluids can be captured strongly by the phase-field models.

Phase field theory

In an incompressible two-phase system consisting of A and B phases with densities of ρ_A and ρ_B, the free energy functional is employed to describe the dynamics of two-phase flow (Kendon et al.,2001):

(1) $\mathit{F}\left(\mathit{\rho },\mathit{\varphi },▿\mathit{\varphi }\right)=\int \left[\Psi \left(\mathit{\varphi }\right)+\frac{1}{2}\mathit{k}{\left|▿\mathit{\varphi }\right|}^2+\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^2\ln \mathit{\rho }\right]\mathit{dV}$

where $\mathit{\rho }={\mathit{\rho }}_{\mathit{A}}+{\mathit{\rho }}_{\mathit{B}}$ is the total density and the order parameter $\varphi =\raisebox{1ex}{$\left({\mathit{\rho }}_{\mathit{A}}-{\mathit{\rho }}_{\mathit{B}}\right)$}\!\left/ \!\raisebox{-1ex}{$\mathit{\rho }$}\right.$ describes the normalized density difference in the two fluids. $\mathit{\psi }\left(\varphi \right)=\frac{1}{4}\mathit{a}{\left({\varphi }^2-1\right)}^2$ is the bulk free-energy density. The term $\frac{1}{2}\mathit{K}{\left|▽\varphi \right|}^2$ is the interfacial energy density with k related to the interfacial tension. The final term in the free energy is introduced to enforce incompressibility. The chemical potential µ is given by Kendon et al. (2001),

(2) $\mathit{\mu }=\raisebox{1ex}{$\partial \mathit{F}$}\!\left/ \!\raisebox{-1ex}{$\partial \mathit{\varphi }={\mathit{\psi }}^{\prime }\left(\mathit{\varphi }\right)-\mathit{k}{▿}^2\mathit{\varphi }$}\right.=\mathit{a}\mathit{\varphi }\left({\mathit{\varphi }}^2-1\right)-\mathit{k}{▿}^2\mathit{\varphi }$

The interfacial tension is given by,

(3) $\mathit{\sigma }=\frac{4\mathit{k}}{3\mathit{\zeta }}$

where ζ is a parameter proportional to the thickness of the interface, defined as,

(4) $\mathit{\zeta }=\sqrt{\raisebox{1ex}{$2\mathit{k}$}\!\left/ \!\raisebox{-1ex}{$\mathit{a}$}\right.}$

From Eqs. (3) and (4); the interfacial tension and the interface thickness can be adjusted by the parameters k and a. In the numerical solution, the interfacial thickness parameter ζ is a free parameter which should be set small to keep a sharp interface. However, ζ cannot be set too small in order to prevent instability. The dynamics of two phase flow can be described by the continuity equation, the N-S equation and the Cahn-Hilliard equation as follows (Zhou et al., 2010):

(5) $▿.\mathit{u}=0$

(6) $\mathit{\rho }\left(\frac{\partial \mathit{u}}{\partial \mathit{t}}+\mathit{u}.▿\mathit{u}\right)=-▿\mathit{P}+▿.\left[\mathit{\eta }\left(▿\mathit{u}+▿{\mathit{u}}^{\mathit{T}}\right)\right]+{\mathit{F}}_{\mathit{s}}$

(7) $\frac{\mathit{\delta }\mathit{\varphi }}{\mathit{\delta }\mathit{t}}+\mathit{u}.▿\mathit{\varphi }=▿.\left(\mathit{M}▿\mathit{\mu }\right)$

where M is the Cahn-Hilliard mobility and $\mathit{P}=\mathit{\rho }{\mathit{C}}_{\mathit{s}}^2+\varphi \mathit{\mu }$ is the modified pressure. In the Lattice Boltzmann model, the interfacial force F_s can be simply considered as a forcing term.

Lattice Boltzmann method

Two distribution functions f_i and g_i are used to calculate the order parameter $\varphi \left(\stackrel{r}{\mathit{x}}\right)$ and the flow field $\stackrel{r}{\mathit{u}}\left(\mathit{x}\right)$ at each lattice site, respectively. The macroscopic variables are related to the distribution functions by:

(8) $\mathit{\varphi }\left(\stackrel{r}{\mathit{x}},\mathit{t}\right)=\sum _{\mathit{i}}{\mathit{g}}_{\mathit{i}}\left(\stackrel{r}{\mathit{x}},\mathit{t}\right)$

(9) $\mathit{\rho }\left(\stackrel{r}{\mathit{x}},\mathit{t}\right)=\sum _{\mathit{i}}{\mathit{f}}_{\mathit{i}}\left(\stackrel{r}{\mathit{x}},\mathit{t}\right)$

(10) $\mathit{\rho }\stackrel{r}{\mathit{u}}\left(\stackrel{r}{\mathit{x}},\mathit{t}\right)=\sum _{\mathit{i}}{\mathit{f}}_{\mathit{i}}\left(\stackrel{r}{\mathit{x}},\mathit{t}\right)\stackrel{r}{{\mathit{e}}_{\mathit{i}}}+\raisebox{1ex}{$\mathit{\mu }▿{\mathit{\varphi }\mathit{\delta }}_{\mathit{t}}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$

where $\mathit{F}=\mathit{\mu }▽\varphi$ is the interfacial force. The evolution equations are:

(11) ${\mathit{f}}_{\mathit{a}}^{\mathit{\sigma }}\left(\mathit{x}+{\mathit{e}}_{\mathit{a}}\Delta \mathit{t},\mathit{t}+\Delta \mathit{t}\right)={\mathit{f}}_{\mathit{a}}^{\mathit{\sigma }}\left(\mathit{x},\mathit{t}\right)-\frac{\left[{\mathit{f}}_{\mathit{a}}^{\mathit{\sigma }}\left(\mathit{x},\mathit{t}\right)-{\mathit{f}}_{\mathit{a}}^{\mathit{\sigma },\mathit{eq}}\left(\mathit{x},\mathit{t}\right)\right]}{{\mathit{\tau }}_{\mathit{f}}}+{\mathit{F}}_{\mathit{a}}$

(12) ${\mathit{g}}_{\mathit{a}}^{\mathit{\sigma }}\left(\mathit{x}+{\mathit{e}}_{\mathit{a}}\Delta \mathit{t},\mathit{t}+\Delta \mathit{t}\right)={\mathit{g}}_{\mathit{a}}\left(\mathit{x},\mathit{t}\right)-\frac{\left[{\mathit{g}}_{\mathit{a}}\left(\mathit{x},\mathit{t}\right)-{\mathit{g}}_{\mathit{a}}^{\mathit{eq}}\left(\mathit{x},\mathit{t}\right)\right]}{{\mathit{\tau }}_{\mathit{g}}}$

where

(13) ${\mathit{F}}_{\mathit{a}}=\left(1-\frac{1}{2{\mathit{\tau }}_{\mathit{f}}}\right){\mathit{w}}_{\mathit{a}}\left[\frac{\stackrel{r}{{\mathit{e}}_{\mathit{a}}}-\stackrel{r}{\mathit{u}}}{{\mathit{c}}_{\mathit{s}}^2}+\frac{\stackrel{r}{{\mathit{e}}_{\mathit{a}}}.\stackrel{r}{\mathit{u}}}{{\mathit{c}}_{\mathit{s}}^4}\stackrel{r}{{\mathit{e}}_{\mathit{a}}}\right].\stackrel{r}{{\mathbf{F}}_{\mathit{s}}}{\mathit{\delta }}_{\mathit{t}}$

For a two-dimensional model (D₂Q₉), the lattice velocities are chosen to be ${\mathit{c}}_0=0$ and ${\mathit{c}}_{\mathit{i}}=\left[\cos \left(\left(\mathit{i}-1\right)\frac{\mathit{\pi }}{2}\right),\sin \left(\left(\mathit{i}-1\right)\frac{\mathit{\pi }}{2}\right)\right],\mathit{i}=1,2,3,4$, the lattice speed c is defined by $\mathit{c}=\raisebox{1ex}{${\mathit{\delta }}_{\mathit{x}}$}\!\left/ \!\raisebox{-1ex}{${\mathit{\delta }}_{\mathit{t}}$}\right.$, where δ_x is the lattice distance, and δ_t is the time step. The speed of sound c_s can be related to c by ${\mathit{c}}_{\mathit{s}}=\raisebox{1ex}{$\mathit{c}$}\!\left/ \!\raisebox{-1ex}{$\sqrt{3}$}\right.$. τ_f and τ_g are independent relaxation parameters. We set ${\mathit{\tau }}_{\mathit{g}}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\left(3-\sqrt{3}\right)$}\right.$ to minimize numerical errors (Elzanfaly et al., 2015) and w_i is the weight factor with w₀ =4/9, w_1-4 =1/9, w_5-8 =1/36. The equilibrium distribution functions are given by Liu and Zhang (2009; 2011):

(14a) ${\mathit{f}}_{\mathit{a}}^{\mathit{eq}}={\mathit{w}}_{\mathit{a}}\left[{\mathit{A}}_{\mathit{a}}+\mathit{\rho }\left(3\stackrel{r}{{\mathit{e}}_{\mathit{a}}}.\stackrel{r}{\mathit{u}}+\frac{9}{2}{\left(\stackrel{r}{{\mathit{e}}_{\mathit{i}}}.\stackrel{r}{{\mathit{u}}_{\mathit{\sigma }}}\right)}^2-\frac{3}{2}\stackrel{r}{\mathit{u}}.\stackrel{r}{\mathit{u}}\right)\right]$

(14b) ${\mathit{g}}_{\mathit{a}}^{\mathit{eq}}={\mathit{w}}_{\mathit{a}}\left[{\mathit{B}}_{\mathit{a}}+\mathit{\varphi }\left(3\stackrel{r}{{\mathit{e}}_{\mathit{a}}}.\stackrel{r}{\mathit{u}}+\frac{9}{2}{\left(\stackrel{r}{{\mathit{e}}_{\mathit{i}}}.\stackrel{r}{{\mathit{u}}_{\mathit{\sigma }}}\right)}^2-\frac{3}{2}\stackrel{r}{\mathit{u}}.\stackrel{r}{\mathit{u}}\right)\right]$

(15) $\begin{array}{c}{\mathit{A}}_0=\left\{\begin{array}{c}\frac{\mathit{P}}{{\mathit{c}}_{\mathit{s}}^2}\left(\mathit{i}>0\right)\hfill \\ \frac{\left[\mathit{\rho }-\frac{\left(1-{\mathit{w}}_0\right)\mathit{P}}{{\mathit{c}}_{\mathit{s}}^2}\right]}{{\mathit{w}}_0}\left(\mathit{i}=0\right)\hfill \end{array}\right.,{\mathit{A}}_{1-8}=\frac{\mathit{P}}{{\mathit{C}}_{\mathit{s}}^2}\hfill \\ \hfill \\ {\mathit{B}}_{\mathit{a}}=\left\{\begin{array}{c}\frac{\Gamma \mathit{\mu }}{{\mathit{c}}_{\mathit{s}}^2}\left(\mathit{i}>0\right)\hfill \\ \frac{\left[\mathit{\varphi }-\frac{\left(1-{\mathit{w}}_0\right)\Gamma \mathit{\mu }}{{\mathit{c}}_{\mathit{s}}^2}\right]}{{\mathit{w}}_0}\left(\mathit{i}=0\right)\hfill \end{array}\right.,{\mathit{B}}_{1-8}=\frac{\Gamma \mathit{\mu }}{{\mathit{C}}_{\mathit{s}}^2}\hfill \end{array}$

where Γ is the tunable parameter regarding the Cahn-Hillard mobility,

$\mathit{M}=\Gamma {\mathit{\delta }}_{\mathit{t}}\left({\mathit{\tau }}_{\mathit{g}}-\frac{1}{2}\right)$

Considering the viscosities of the two phases, the following linear averaged viscosity is used as the mixture viscosity (Van Der Sman et al., 2006):

(16) $\mathit{\eta }\left(\mathit{\varphi }\right)=\frac{1-\mathit{\varphi }}{2}{\mathit{\eta }}_{\mathit{A}}+\frac{1+\mathit{\varphi }}{2}{\mathit{\eta }}_{\mathit{B}}$

where ${\mathit{\eta }}_{\mathit{A}}$ and ${\mathit{\eta }}_{\mathit{B}}$ are the viscosities of the two phases. The local relaxation time τ_f can be calculated from the local viscosity:

(17) $\mathit{\eta }=\mathit{\rho }{\mathit{c}}_{\mathit{s}}^2{\mathit{\delta }}_{\mathit{t}}\left({\mathit{\tau }}_{\mathit{f}}-\frac{1}{2}\right)$

As mentioned in Liu and Zhang (2009), the Navier-Stokes equations converge to the classical sharp interface behavior as the interface thickness reduces toward zero along with the diffusivity 1/Pe, where Pe is the Peclet number:

(18) $\mathit{Pe}=\raisebox{1ex}{$\mathit{UL}$}\!\left/ \!\raisebox{-1ex}{$\mathit{Ma}$}\right.$

where U and L are the characteristic velocity and length of the system. A suitable Peclet number should be chosen in the simulation to balance the effects of convection and diffusion at the interface. Therefore, in the present simulation, ζ =(1.5 * lattice size) is selected. In this condition, the thickness of the interface is about 4~5 lattices. With this amount of ζ and having a surface tension σ, using Eqs. (3) and (4), the value of a is achieved. M is the Cahn-Hilliard mobility and can be obtained with ${\mathit{\tau }}_{\mathit{g}}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\left(3-\sqrt{3}\right)$}\right.$ and Γ,. By using M;a; U and L; the Peclet number is O (10)~O (100).

DROPLET FORMATION IN A T-JUNCTION MICROCHANNEL

The purpose of this study is to investigate the formation mechanism of the droplet in microchannels. A T-shaped microchannel is selected, which includes two inputs and one output (for disperse and continuous phases) (see Fig. 1). As shown in Fig. 1; The continuous phase flows through the main channel (with width W_c) and the dispersed phase is injected through the lateral channel (with width W_d).

Figure 1
(a) Geometry of the T-junction microchannel; (b) Droplet formation in the T-junction; d and S are related to neck and front of the liquid film; respectively.

The no-slip boundary condition is applied at the solid walls using halfway bounceback, which can prevent the boundary “mass leakage", especially for a flow with small velocity (Succi, 2001). The constant flow rate or pressure boundary condition can be imposed based on the work of Zou and He (1997). It is assumed that the inlet and outlet boundaries contain only a single-component of the phases, where the unknown g_i can be obtained by using the method proposed by Hao and Cheng (2009). For example, if the inlet boundary is assumed to be perpendicular to the y -direction with the lattice velocities e₁, e₅ and e₁, g₅ and g₈ are unknown after the streaming step. In order to ensure a prescribed order parameter ${\varphi }_{\mathit{in}}$ at the inlet, these unknown PDFs must satisfy:

(19) ${\mathit{\varphi }}_{\mathit{in}}={\mathit{g}}_1+{\mathit{g}}_2+{\mathit{g}}_3+{\mathit{g}}_4+{\mathit{g}}_5+{\mathit{g}}_6+{\mathit{g}}_7+{\mathit{g}}_8+{\mathit{g}}_0$

(20) ${\mathit{g}}_1+{\mathit{g}}_5+{\mathit{g}}_8={\mathit{\varphi }}_{\mathit{in}}-\left({\mathit{g}}_4+{\mathit{g}}_2+{\mathit{g}}_6+{\mathit{g}}_7+{\mathit{g}}_3+{\mathit{g}}_0\right)$

(21) ${\mathit{\varphi }}^{*}={\mathit{\varphi }}_{\mathit{in}}-\left({\mathit{g}}_4+{\mathit{g}}_2+{\mathit{g}}_6+{\mathit{g}}_7+{\mathit{g}}_3+{\mathit{g}}_0\right)$

Assuming that the g_i (i = 1,5,8) are distributed by their weight factors w_i, we have:

(22) ${\mathit{g}}_{\mathit{i}}=\frac{{\mathit{w}}_{\mathit{i}}{\mathit{\varphi }}^{*}}{{\mathit{w}}_1+{\mathit{w}}_5+{\mathit{w}}_8}\mathit{i}=1,5,8$

In the present numerical simulation, a constant velocity is used for the inlet boundary as proposed by Zho and He (1997) and the open boundary condition (Hao and Cheng, 2009) is selected for the outlet of the microchannel. The bounce back boundary condition (Sukop, 2006) is also used for the walls. Wetting properties are usually characterized by the contact angle on a surface. For a desired static contact angle ${\mathit{\theta }}_{\mathit{s}}$, the following formula is used to assign the order parameter ${\varphi }_{\mathit{s}}$ to the solid lattice sites next to the wall (Van der Graaf et al., 2006).

(23) $\cos \left({\mathit{\theta }}_{\mathit{s}}\right)=\frac{1}{2}{\mathit{\varphi }}_{\mathit{s}}\left(3-{\mathit{\varphi }}_{\mathit{s}}^2\right)$

MODEL VALIDATIONS

Stationary droplet: test case I

According to the Laplace equation, the pressure difference Δp across the interface of a droplet (with radius R_d) in equilibrium can be explained as:

(24) $\Delta \mathit{p}=\frac{2\mathit{\sigma }}{\mathit{R}}$

where Δp denote the pressure difference across the inside and outside of the droplet, and σ is the surface tension between the two fluids. Simulations were conducted for different sizes of droplet. The droplet is located at the center of the lattice domain with 120(120 lattices, for a fixed interfacial tension parameter. The pressure difference at equilibrium is compared to Laplace's law. The obtained results are shown in Fig. 2(a) by using ${\mathit{\tau }}_{\mathit{f}}=1.0$, $\mathit{\sigma }=0.01$ and $\Gamma =4$. Figure 2(a) shows that the results can be fitted with a straight line with a slope of 0.01 (LB lattice), which is the surface tension between the two fluids. The validation was verified for viscosity ratios of 1 and 10. Equation (25) can also be used to verify the accuracy of the numerical results, which can be written as (Liu and Zhang, 2011):

(25) $\mathit{\varphi }\left(\mathit{x},\mathit{y}\right)=\tanh \frac{\mathit{R}-\sqrt{{\left(\mathit{x}-{\mathit{x}}_0\right)}^2+{\left(\mathit{y}-{\mathit{y}}_0\right)}^2}}{\mathit{\zeta }}$

where x₀ and y₀ are the coordinates of the center of the droplet. Figure 2(b) displays the order parameter as a function of the distance from the droplet center, which is in good agreement with the theoretical profile given by Eq. (25). This shows that our method can capture the interface correctly.

Figure 2
(a) Validation of Laplace equation for a droplet in an infinite domain; (b) Comparison of the profile of order parameter along the cross section of a droplet between the present model and the theoretical profile presented in Liu and Zhang (2011). The solid line is the theoretical profile and discrete symbols are the present model.

Evaluation of the contact angle: test case II

Different contact angles can be taken into consideration by adjusting the interaction strength of the two immiscible fluids with the walls of the microchannel. For this purpose, the present code (free energy LBM code) is employed for a single droplet in equilibrium inside the channel by choosing different static contact angles. Figure 3(a) shows droplet shapes for three different contact angles. The contact angles are also measured from the numerical contours and are compared with the theoretical solution (i.e., Eq. (26) of Van der Graaf et al. (2006)). The comparison is shown in Fig. 3(b). The results exhibit a linear relationship between the order parameter (ø_wall / ø₀) and the contact angle. Figure 3(b) shows that the numerical results are in good agreement with the analytical solution (Van der Graaf et al. ,2006).

(26) $\cos \left(\mathit{\theta }\right)=\frac{3}{2}\frac{{\mathit{\varphi }}_{\mathit{wall}}}{{\mathit{\varphi }}_0}\left(1-\frac{1}{3}{\left(\frac{{\mathit{\varphi }}_{\mathit{wall}}}{{\mathit{\varphi }}_0}\right)}^2\right)$

Figure 3
(a) Droplet contour for different contact angles. (b) Contact angle as a function of φ_wall / φ₀

Single droplet in shear flow: test case III

Taylor deformation is considered to investigate the droplet formation behavior of the multiphase model. A droplet is placed between two shearing plates. The plates are moving in the opposite direction to obtain linear shear in The Stokes regime. The droplet deformation is a function of the shear rate (expressed as the Capillary number) at a constant Peclet number. The definition of the Reynolds number (Re), Capillary number (Ca), Peclet number (Pe) and shear rate (γ) are as follows:

(27) $\mathit{Re}=\frac{\mathit{\gamma }{\mathit{R}}^2\mathit{\rho }}{\mathit{\eta }}$

(28) $\mathit{Ca}=\frac{\mathit{\gamma }\mathit{R}\mathit{\eta }}{\mathit{\sigma }}$

(29) $\mathit{\gamma }=\raisebox{1ex}{$2\mathit{U}$}\!\left/ \!\raisebox{-1ex}{$\mathit{H}$}\right.$

(30) $\mathit{Pe}=\frac{\mathit{\gamma }\mathit{R}\mathit{\zeta }}{\mathit{aM}}$

where R is the initial radius of the droplet, U being the velocity of the moving plate, H is the channel height, ρ is the density and η is the dynamic viscosity of droplet. λ Is the viscosity ratio between droplet and fluid. The simulations are performed at Re=0.1 and for a droplet with a radius of 32 lattice cells in a computational domain of 256*128 lattices. In the steady condition the shape of the droplet is elliptic, which is usually characterized by the deformation parameter (Df):

(31) $\mathit{Df}=\frac{\mathit{L}-\mathit{B}}{\mathit{L}+\mathit{B}}$

where L is the major and B is the minor axis of the ellipse. For a drop in the Stokes regime and for low Capillary numbers, it is expected that Df follows the Taylor relation as (Roths et al., 2002):

(32) $\mathit{Df}=\mathit{f}\left(\mathit{\lambda }\right)\mathit{Ca}$

Figure 4(a) shows the steady state droplet shapes for various Capillary numbers. All profiles cross two points that are in good agreement with the numerical results of the boundary integral method (Zhou and Pozrikidis, 1993). Figure 4(b) gives the deformation parameter as a function of the Capillary number, where Df =1.5Ca is obtained based on the present simulation results. A comparison with the VOF model (Zhou and Pozrikidis, 1993), the Swift et al. (1996) model and the results of Liu and Zhang (2011) is also reported here.

Figure 4
(a) Stable droplet shapes at different Capillary numbers Ca={0.05,0.1,0.2,0.3}. By increasing the Capillary number the ellipse is stretched and the length of the large-diameter is increased. (b)Taylor deformation parameter Df as a function of Capillary number.

Figure 5 shows the comparison of the droplet shape under simple shear flow for (a) Ca=0.9 and (b) Ca=1.0. In Fig. 5, at Ca=0.9, the maximal deformation is not enough to “pinch-off" the droplet. The break up is observed at the critical Capillary number (i.e., Ca=1.0) and two small droplets are detached from the bulk. These profiles agree well with the available results in Zhou and Pozrikidis (1993).

Figure 5
Droplet deformation under shear flow for (a) Ca=0.9 and (b) Ca=1.0.

RESULTS AND DISCUSSION

The droplet formation in a T-shape microchannel can be defined by several dimensionless numbers, the interfacial tension γ, the inlet volumetric flow rates (Q_d and Q_c), fluid viscosities (η_c and η_d) and fluid densities (ρ_c and ρ_d) where the subscripts 'c' and 'd' refer the continuous and dispersed phases, respectively. One of the important dimensionless numbers is the Capillary number (Ca) of the continuous phase, $\mathit{Ca}=\frac{{\mathit{U}}_{\mathit{c}}{\mathit{\eta }}_{\mathit{c}}}{\mathit{\gamma }}$, which is the ratio of the viscous force and the surface tension force, where U_c is the average inlet velocity, η_c is the dynamic viscosity of the continuous phase, and γ is the interfacial tension. Re is the Reynolds number describing the ratio of inertia to viscous stresses, $\mathit{Re}=\frac{\mathit{\rho }{\mathit{U}}_{\mathit{c}}{\mathit{w}}_{\mathit{c}}}{{\mathit{\eta }}_{\mathit{c}}}$. The ratio of flow rates ($\mathit{R}=\frac{{\mathit{Q}}_{\mathit{d}}}{{\mathit{Q}}_{\mathit{c}}}$) and the viscosity ratio of the two immiscible fluids ($\mathit{\lambda }=\frac{{\mathit{\eta }}_{\mathit{d}}}{{\mathit{\eta }}_{\mathit{c}}}$) are other key dimensionless parameters. For typical microfluidic systems, the Bond number is small due to the small density difference of the two liquids.

Effects of capillary number, viscosity and flow rate ratios

In this section, the dynamics of droplet break-up in the microchannels are studied. Moreover, the interaction of dominant forces and their effects on the breakup process are investigated. These forces are surface tension forces, shear stresses and pressure. The droplet break-up is a result of interaction of these forces in different conditions. While the continuous phase is flowing in the main channel, the droplet penetrates towards downstream due to the applied tension and pressure. In the early steps of droplet formation, the dispersed phase enters the main channel and moves towards the downstream due to the continuous phase flow direction. When it occupies the channel's width, the droplet's bottleneck starts to narrow. Finally, the neck is disconnected and the droplet is formed. This process contains four steps;

1. Dispersed phase enters the main channel;

2. The droplet grows in the y direction and fully occupies the channel width;

3. The droplet penetrates towards down streams (the bottleneck is narrowing);

4. The droplet is disconnected and moves towards downstream;

Fig. 6 shows these four stages from the entering of the dispersed phase into the main channel until droplet detachment for two different Capillary numbers of 0.005 and 0.02. The above mentioned stages can be observed for all Capillary numbers, but for small Capillary numbers (see Fig. 6(a)), the edge of the T-junction is the detachment point. By increasing the Capillary number (see Fig. 6(b)), the detachment point shifts from the edge of the T-junction towards the downstream of the flow far from the edge of the T-junction.

Figure 6
Droplet formation from penetration until detachment for two Capillary numbers; Ca=0.005 and Ca =0.02; at a fixed viscosity ratio of 0.5 and considering a fixed flow rate ratio of R =0.3

Figure 7 shows temporal variations of S (S is shown in Fig. 1). The droplet formation contains three levels: (a) growing in both the x and y directions, (b) growing in only the x direction and finally (c) detachment. In growing in the x and y directions, the plug shape grows in both the x and y directions. It advances in the main channel until the tip of the drop connects the lower wall of the microchannel. At this point, droplet growth in the y direction is stopped. In this condition the formed droplet continues growing in the x direction. As a result of continuous phase pressure, the neck (i.e. =d) is stretched towards downstream and starts to narrow until the droplet starts to disconnect. The time required from droplet formation to disconnection is called the ''detachment time''. As shown in Fig. 7, the longest time is related to the growing process in both the x and y directions. The first increasing section of the curves in Fig. 7 is related to the growing in both the x and y directions. The flat section of the curves in this figure reflects the growth in only the x direction. The second increasing section shows the detachment stage. In other words, the neck of the liquid film is narrowing in the third section of the curves.

Figure 7
The process of droplet growth in the microchannel; (a) Small capillary numbers; (b) Larger capillary numbers; (c) Duration of step1 for different Capillary numbers for two flow rate ratios of 0.3 and 0.5 at fixed viscosity ratio of 0.5 and (d) Duration of step2 for different Capillary numbers for two flow rate ratios of 0.3 and 0.5 at fixed viscosity ratio of 0.5

Figure 7(a) shows the variations of S with time for different Capillary numbers. The Capillary numbers selected are small values in this figure, i.e., Ca<0.02. Two flow rate ratios are also used in the simulation. As can be seen, by increasing the flow rate ratio (R) and also increasing the Capillary number, the duration of droplet formation and the length of the x-growth stage are obviously decreased. Figure 7(b) shows temporal variations of S for larger Capillary numbers (Ca>0.02) at a fixed flow rate ratio (i.e., R =0.5). As can be seen, the droplet penetrates toward the downstream; therefore the length of the x-growth stage is increased. Figures 7(c) and 7(d) show the duration of the x- and y- growth stage (step1) and the duration of the x-growth stage (step2) in the microchannel, respectively. The duration of both steps is obviously decreased by increasing the Capillary number. For small Capillary numbers, i.e., Ca<0.02, the lengths of step1 and step2 are decreased by increasing the flow rate ratio (R). In contrast, for larger Capillary numbers (Ca>0.02) the curves are independent to flow rate ratio.

Capillary number for the continuous phase and flow rate ratio are important parameters to distinguish the mechanism of break up in the microchannels. Three flow regimes are categorized in the two phase flow in the microchannels: (a) squeezing, (b) dripping and (c) jetting. In the squeezing regime, the dispersed phase is likely to occupy the whole width of the channel and the edge of the T-junction is the detachment point. The pressure applied on the dispersed phase from the continuous phase plays the main role in droplet break-up. In the dripping regime, the viscous stresses are important as mentioned by Garstecki et al. (2006). In the jetting regime, the detachment point shifts from the edge of the T-junction towards downstream of the flow far from the edge of the T-junction. Moreover, the detachment point shifts towards downstream at larger Capillary numbers, flow rate ratios and viscosity ratios and forms a stable jet.

The process of detachment for different Capillary numbers (Ca=0.005-0.07) and for two viscosity ratios of 0.5 and 1.0 is shown in Fig. 8.

Figure 8
Detachment point for different Capillary numbers (Ca=0.005, 0.02, 0.04, 0.05, 0.07) and for two viscosity ratios of 0.5 (left frame) and 1.0 (right frame). The flow rate ratio is equal to 0.5.

As can be seen in Fig. 9, by increasing the Capillary number, the detachment point shifts towards downstream and the detachment length is increased. Figure 9 also illustrates that the detachment length is increased by increasing the flow rate ratio and viscosity ratio and forms a stable jet.

Figure 9
The droplet detachment length: (a) different Capillary numbers considering two flow rate ratios of 0.3 and 0.5 at a fixed viscosity ratio of 0.5; (b) different Capillary numbers considering two viscosity ratios of 0.5 and 1.0 at a flow rate ratio of 0.5

According to Fig. 10(a), the detachment time is increased by decreasing the Capillary number. For Ca ≥ 0.04, the detachment time is independent of the flow rate ratios. For Capillary numbers smaller than Ca=0.04; the detachment time is decreased by increasing the flow rate ratio. As shown in Fig. 10(b), decreasing the Capillary number leas to an increase in the detachment time; nonetheless, the viscosity ratios do not have a significant effect on the detachment time.

Figure 10
The droplet detachment time; (a) different Capillary numbers considering different flow rate ratios of 0.15, 0.3 and 0.5 at a fixed viscosity ratio of 0.5 (b) different Capillary numbers considering viscosity ratios 0.25, 0.5 and 1.0 at a flow rate ratio of 0.5

Figure 11 illustrates the formation of droplets for two different viscosity ratios of 0.5 and 1.0 for a constant flow rate of 0.3. According to Fig.11, in low Capillary numbers, the droplet fills all the width of the channel and detachment happens in the sharp part of the T-shaped junction for different viscosity ratios; so the viscosity ratio has no effect on the size of the droplet. By increasing the Capillary number, a larger viscosity ratio leads to pushing the position of detachment to the downstream. In this case, the break-up of the droplet is the result of competition between viscous and capillary forces.

Figure 11
The droplet generation for R =0.3 for different Capillary numbers of Ca=0.005, 0.01, 0.04, 0.07 for two viscosity ratios of 0.5 and 1.0

Figure 12(a) shows the variations of the size of the droplet at different capillary numbers for viscosity ratios of 1.0, 0.5 and 0.25. For Ca ≤ 0.04, the size of the droplet is independent of the viscosity ratio, but after the critical Capillary number (i.e., Ca=0.04), the droplet size is completely affected by viscosity ratio. In other words, for larger viscosity ratios, a smaller droplet is observed in the microchannel. According to Fig. 12(b), for Capillary numbers smaller than 0.02, the size of the droplet is completely influenced by the flow rates of the input phases. In other words, a large flow rate ratio leads to an increase in the droplet size. But for Capillary numbers larger than the critical value (i.e. Ca=0.02), the size of droplet is independent of flow rate ratios.

Figure 12
Droplet length as a function of the Capillary numbers (a) for viscosity ratios of 0.25, 0.5 and 1.0 at a fixed flow rate ratio of R =0.3; (b) for flow rate ratios of 0.5, 0.3 and 0.15 at a fixed viscosity ratio of 0.5.

According to Fig. 13, the time between two consecutive drops is decreased by increasing the Capillary number. Moreover, this time is decreased by increasing the flow rate ratio until Ca=0.04. After this Capillary number (i.e., Ca=0.04), the flow rate ratios have no effect on the time between two consecutive droplets.

Figure 13
The time between two consecutive drops for different Capillary numbers by considering three flow rate ratios of R =0.15, 0.3, 0.5 at a fixed viscosity ratio of 0.5

Figure 14 shows the dependence of the size of the droplet for different flow rate ratios in low Capillary numbers. In low Capillary numbers, typically less than 0.01, the droplet formation is a result of the pressure force which is applied by the continuous phase on the dispersed phase. Shear stress and the viscous forces at the interface can be ignored in comparison with the surface tension. The pressure difference between two phases plays a main role in the break-up process. To overcome the forces at the interface of two phases, the droplet blocks all the T-shaped part of the main channel, hence it blocks flow of the continuous phase. At this time, the length of droplet is almost equal to the width of the channel, which results in an increase of pressure on the upstream side of the forming droplet. The neck of the droplet (=d) is pressed and this growth in the droplet continuous in flow rate ratio (Q_d). The tip of the dispersed phase is also elongated until the neck of the dispersed phase is squeezed. The rate of reducing the thickness of the neck of the dispersed phase is almost equal to the mean speed of the continuous phase ${\mathit{U}}_{\mathit{squeeze}}\approx \raisebox{1ex}{${\mathit{Q}}_{\mathit{c}}$}\!\left/ \!\raisebox{-1ex}{$\mathit{hw}$}\right.$. During this process, the droplet is elongated with a rate of ${\mathit{U}}_{\mathit{growth}}\approx \raisebox{1ex}{${\mathit{Q}}_{\mathit{c}}$}\!\left/ \!\raisebox{-1ex}{$\mathit{hw}$}\right.$. Therefore, the final length of the droplet can be expressed as $\mathit{L}\approx \mathit{w}+\raisebox{1ex}{${\mathit{dQ}}_{\mathit{d}}$}\!\left/ \!\raisebox{-1ex}{${\mathit{Q}}_{\mathit{c}}$}\right.$. In other words, one can write the following relation for the length of the droplet:

(33) $\frac{\mathit{L}}{{\mathit{W}}_{\mathit{c}}}=1+\mathit{\alpha }\frac{{\mathit{Q}}_{\mathit{d}}}{{\mathit{Q}}_{\mathit{c}}}$

which has a suitable consistency with the results of Garstecki et al. (2006), the relation of Christopher et al. (2008) and the results presented by experimental data. L is the length of droplet, W_c is the width of main channel and α is a constant. In this part, in order to recheck the dynamics of two phase flow, the proposed LBM code is validated by the Comsol software (COMSOL Multiphysics 5.0). The geometry of the microchannel is simulated by Comsol software in order to check the length of the droplet for different flow rate ratios. The results of Tang et al. (2014) and Liu and Zhang (2009) are also presented in Fig. 14.

Figure 14
Comparison of normalized droplet length as a function of flow rate ratio (Q_d/Q_c) between the current model, Tang et al. (2014), Liu and Zhang (2009) and Comsol simulation (Finite Element).

Effect of width ratio of the two inlets

The ratio of the main channel to the lateral channel width can change the length of the droplet. In this section, the effects of different width ratios (W_d/W_c) on the length of the droplet are examined. The contact angle between fluid and walls is considered to be nearly 180ᵒ. Densities of the two phases are the same with the value of 1 (in lattice units) and the viscosity ratio between the two phases is considered to be 0.5. Figure 15 shows the variations of the size of the droplet for different width ratios. Different flow rate ratios of 0.15, 0.3 and 0.5 and in two Capillary numbers of Ca=0.005 and Ca=0.01 are selected for this simulation. Moreover, the width of the main channel is kept constant.

Figure 15
The droplet length at different width ratios considering different flow rate ratios at Capillary numbers: (a) Ca=0.01 and (b) Ca=0.005

As shown in Fig. 16, by decreasing the width of the lateral channel, longer droplets are formed inside the microchannel.

Figure 16
Flow patterns for various width ratios at Ca=0.01, R= 0.3 and viscosity ratio of 0.5.

Figure 17
The effect of width ratio on droplet length for different Capillary numbers, R=0.3 and viscosity ratio of 0.5.

CONCLUSIONS

The multiphase lattice Boltzmann method was used to simulate the droplet formation and breakup process of immiscible fluids in T-junction microchannels. The steps of droplet formation in the microchannel were analyzed in the paper. The effect of capillary number and flow rate ratio on the different steps of droplet formation was investigated. The influences of Capillary number, flow rate ratio, viscosity ratio and geometrical parameters on droplet break up, droplet size and detachment time are systematically studied.

• By increasing the flow rate ratio; the duration of droplet formation and the length of the x-growth stage were decreased for small capillary numbers. For larger capillary numbers, the droplet penetrated toward the downstream; therefore, the length of the x-growth stage was increased.

• The start of the detachment process in the microchannel was also reported.

• The detachment time was increased by decreasing the capillary number. For Ca>0.02, the detachment time was independent of the flow rate ratios.

• The effects of viscosity ratios on the detachment time were not significant in comparison to the capillary number.

• For Ca<0.04, the size of the droplet is independent of the viscosity ratio, but after this critical Capillary number, the size of the droplet was varied by the viscosity ratio.

• The time between two consecutive drops was also decreased by increasing the Capillary number. Moreover, this time was decreased by increasing the flow rate ratio until Ca=0.04. After this Capillary number, the flow rate ratios had no significant effect on the time between two consecutive droplets.

• An exhaustive validation study was performed in this two phase study.

NOMENCLATURE
• Ca  Capillary number
• C_s  Speed of sound
• e_i  Lattice speed of particles moving in direction i
• f, g  Particle distribution function
• f_a  Particle equilibrium distribution function
• g_a(eq)  Particle equilibrium distribution function
• p  Pressure
• Pe  Peclet number
• R  Flow rate ratio
• Re  Reynolds number

Greek symbols
• σ  Surface tension
• τ  Relaxation time
• ρ  Density
• µ, η  Dynamic viscosity
• λ  Viscosity ratio
• ϒ  Shear rate

ACKNOWLEDGMENTS

This work has been supported by the Center for International Scientific Studies & Collaboration (CISSC).

REFERENCES

• Alizadeh, M.; Rahni, M. T. and Yazdi, M. E. Numerical investigation of emulsion process in microchannels, using index-function lattice Boltzmann method. Modares Mechanical Engineering, 15(4), 13-22 (2015).
• Anna, S. L., Bontoux, N., and Stone, H. A. Formation of Dispersions Using "flow Focusing" in Microchannels. Applied Physics Letters 82(3), 364-366 (2003). doi: 10.1063/1.1537519.
  » https://doi.org/10.1063/1.1537519
• Azarmanesh, M., and Farhadi, M. The Effect of Weak-inertia on Droplet Formation Phenomena in T-junction Microchannel. Meccanica 51(4), 819-834 (2015). doi: 10.1007/s11012-015-0245-6.
  » https://doi.org/10.1007/s11012-015-0245-6
• Bai, L., Fu, Y., Zhao, S. and Cheng, Y. Droplet formation in a microfluidic T-junction involving highly viscous fluid systems. Chemical Engineering Science 145, 141-148 (2016).‏ doi: 10.1016/j.ces.2016.02.013
  » https://doi.org/10.1016/j.ces.2016.02.013
• Ba Y, Liu H, Li Q, Kang Q, Sun J. Multiple-relaxation-time color-gradient lattice Boltzmann model for simulating two-phase flows with high density ratio. Physical Review E 94: 023310 (2016). doi: 10.1103/PhysRevE.94.023310.
  » https://doi.org/10.1103/PhysRevE.94.023310
• Ba Y, Liu H, Sun J and Zheng R. Three dimensional simulations of droplet formation in symmetric and asymmetric T-junctions using the color-gradient lattice Boltzmann model. International Journal of Heat and Mass Transfer 90: 931-947 (2015).‏ doi: 10.1016/j.ijheatmasstransfer.2015.07.023
  » https://doi.org/10.1016/j.ijheatmasstransfer.2015.07.023
• Cal K. How Does the Type of Vehicle Influence the in Vitro Skin Absorption and Elimination Kinetics of Terpenes?. Archives of Dermatological Research 297(7): 311-315 (2005). doi: 10.1007/s00403-005-0622-4.
  » https://doi.org/10.1007/s00403-005-0622-4
• Christopher G F, Noharuddin N N, Taylor J A, and Anna S L. Experimental Observations of the Squeezing-to-dripping Transition in T-shaped Microfluidic Junctions. Physical Review E 78(3) (2008). doi: 10.1103/physreve.78.036317.
  » https://doi.org/10.1103/physreve.78.036317
• Cramer C, Fischer P, and Windhab E J. Drop Formation in a Co-flowing Ambient Fluid. Chemical Engineering Science 59(15), 3045-3058 (2004). doi: 10.1016/j.ces.2004.04.006.
  » https://doi.org/10.1016/j.ces.2004.04.006
• Cubaud T, Tatineni M, Zhong X, and Ho C M. Bubble Dispenser in Microfluidic Devices. Physical Review E 72(3) (2005). doi: 10.1103/physreve.72.037302.
  » https://doi.org/10.1103/physreve.72.037302
• Elzanfaly E S, Hassan S A, Salem M Y, and El-Zeany B A. Continuous Wavelet Transform, a Powerful Alternative to Derivative Spectrophotometry in Analysis of Binary and Ternary Mixtures: A Comparative Study. Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 151: 945-955 (2015). doi: 10.1016/j.saa.2015.06.100.
  » https://doi.org/10.1016/j.saa.2015.06.100
• Fu B R, and Pan C. Flow Pattern Transition Instability in a Microchannel with CO2 Bubbles Produced by Chemical Reactions. International Journal of Heat and Mass Transfer 48(21-22): 4397-4409 (2005). doi: 10.1016/j.ijheatmasstransfer.2005.05.010.
  » https://doi.org/10.1016/j.ijheatmasstransfer.2005.05.010
• Fu T, Ma Y, Funfschilling D, and Li H Z. Bubble Formation and Breakup Mechanism in a Microfluidic Flow-focusing Device. Chemical Engineering Science 64(10), 2392-2400 (2009). doi: 10.1016/j.ces.2009.02.022.
  » https://doi.org/10.1016/j.ces.2009.02.022
• Garstecki P, Fuerstman M J, Stone H A, and Whitesides G M. Formation of Droplets and Bubbles in a Microfluidic T-junction—scaling and Mechanism of Break-up. Lab on a Chip 6, (3): 437(2006). doi: 10.1039/b510841a.
  » https://doi.org/10.1039/b510841a
• Garstecki P, Stone H A, and Whitesides G M. Mechanism for Flow-Rate Controlled Breakup in Confined Geometries: A Route to Monodisperse Emulsions. Physical Review Letters 94(16) (2005). doi: 10.1103/physrevlett.94.164501.
  » https://doi.org/10.1103/physrevlett.94.164501
• Gouyet J F, Plapp M, Dieterich W, and Maass P. Description of Far-from-equilibrium Processes by Mean-field Lattice Gas Models. Advances in Physics 52(6), 523-638 (2003). doi: 10.1080/00018730310001615932.
  » https://doi.org/10.1080/00018730310001615932
• Gupta A, Sbragaglia M, Belardinelli D and Sugiyama K. Lattice Boltzmann simulations of droplet formation in confined channels with thermo capillary flows. Physical Review E 94, (6): 063302(2016).‏ doi: 10.1103/PhysRevE.94.063302.
  » https://doi.org/10.1103/PhysRevE.94.063302
• Gunstensen AK, Rothman DH, Zaleski S, and Zanetti G. Lattice Boltzmann Model of Immiscible Fluids. Physical Review A43(8), 4320-1327 (1991). Doi: 10.1103/physreva.43.4320.
  » https://doi.org/10.1103/physreva.43.4320
• Günther A, and Jensen KF. Multiphase Microfluidics: From Flow Characteristics to Chemical and Materials Synthesis. Lab on a Chip 6(12), 1487-1503 (2006). doi: 10.1039/b609851g.
  » https://doi.org/10.1039/b609851g
• Hao L, and Cheng P. Lattice Boltzmann Simulations of Liquid Droplet Dynamic Behavior on a Hydrophobic Surface of a Gas Flow Channel. Journal of Power Sources 190(2), 435-446 (2009). doi: 10.1016/j.jpowsour.2009.01.029.
  » https://doi.org/10.1016/j.jpowsour.2009.01.029
• He X, Chen S, and Zhang R. A Lattice Boltzmann Scheme for Incompressible Multiphase Flow and Its Application in Simulation of Rayleigh-Taylor Instability. Journal of Computational Physics152(2): 642-663 (1999). doi: http://dx.doi.org/10.1006/jcph.1999.6257.
  » http://dx.doi.org/10.1006/jcph.1999.6257.
• Hua J, Zhang B, and Lou J. Numerical Simulation of Microdroplet Formation in Coflowing Immiscible Liquids. AIChE Journal 53(10), 2534-2548 (2007). doi: 10.1002/aic.11287.
  » https://doi.org/10.1002/aic.11287
• Kendon VM, Cates ME, Pagonabarraga I, Despla JC, and Bladon P. Inertial effects in three-dimensional spinodal decomposition of a symmetric binary fluid mixture: a lattice Boltzmann study. Journal of Fluid Mechanics 440: 147-203 (2001). doi: 10.1017/S0022112001004682.
  » https://doi.org/10.1017/S0022112001004682
• Liu H, and Zhang Y. Lattice Boltzmann Simulation of Droplet Generation in a Microfluidic Cross-junction. Communications in Computational Physics 9(05), 1235-1256 (2011). doi: 10.4208/cicp.231009.101110s.
  » https://doi.org/10.4208/cicp.231009.101110s
• Liu H, and Zhang Y. Droplet Formation in a T-shaped Microfluidic Junction. Journal of Applied Physics 106(3) (2009). doi: 10.1063/1.3187831.
  » https://doi.org/10.1063/1.3187831
• Liu H, and Zhang Y. Droplet Formation in Microfluidic Cross-junctions. Physics of Fluids 23(8) (2011). doi: 10.1063/1.3615643.
  » https://doi.org/10.1063/1.3615643
• Liu H, Valocchi J.A, Kang Q. Three-dimensional lattice Boltzmann model for immiscible two-phase flow simulations. Physical Review E 85, (4): 046309(2012). doi: 10.1103/PhysRevE.85.046309.
  » https://doi.org/10.1103/PhysRevE.85.046309
• Losey M W, Jackman RJ, Firebaugh S L, Schmidt M A, and Jensen K F. Design and Fabrication of Microfluidic Devices for Multiphase Mixing and Reaction. Journal of Micro electromechanical Systems 11(6), 709-717 (2002). doi: 10.1109/jmems.2002.803416.
  » https://doi.org/10.1109/jmems.2002.803416
• Malekzadeh S, and Roohi E. Investigation of Different Droplet Formation Regimes in a T-junction Microchannel Using the VOF Technique in Open FOAM. Microgravity Science and Technology 27(3), 231-243 (2015). doi: 10.1007/s12217-015-9440-2.
  » https://doi.org/10.1007/s12217-015-9440-2
• Nisisako T, Torii T, and Higuchi T. Droplet Formation in a Microchannel Network. Lab on a Chip 2(1), 24-26 (2002). doi: 10.1039/b108740c.
  » https://doi.org/10.1039/b108740c
• Richter M, Woias P, and Weiβ D. Microchannels for Applications in Liquid Dosing and Flow-rate Measurement. Sensors and Actuators A: Physical 62(1), 480-483 (1997). doi: 10.1016/s0924-4247(97)01486-6.
  » https://doi.org/10.1016/s0924-4247(97)01486-6
• Roths T, Friedrich C, Marth M, and Honerkamp J. Dynamics and Rheology of the Morphology of Immiscible Polymer Blends on Modeling and Simulation. Rheologica Acta 41(3), 211-222 (2002). doi: 10.1007/s003970100189.
  » https://doi.org/10.1007/s003970100189
• Shan X, and Chen H. Lattice Boltzmann Model for Simulating Flows with Multiple Phases and Components. Physical Review E 47(3), 1815-1819 (1993). doi: 10.1103/physreve.47.1815.
  » https://doi.org/10.1103/physreve.47.1815
• Shan X, and Chen H. Simulation of Nonideal Gases and Liquid-gas Phase Transitions by the Lattice Boltzmann Equation. Physical Review E 49(4), 2941-2948(1994). doi: 10.1103/physreve.49.2941.
  » https://doi.org/10.1103/physreve.49.2941
• Shi Y, Tang Gh, and Xia Hh. Lattice Boltzmann Simulation of Droplet Formation in T-junction and Flow Focusing Devices. Computers & Fluids 90: 155-163 (2014). doi: 10.1016/j.compfluid.2013.11.025.
  » https://doi.org/10.1016/j.compfluid.2013.11.025
• Succi S. The Lattice Boltzmann Equation for Fluid Dynamics and beyond. Oxford: Clarendon Press (2001).
• Sugiura S, Nakajima M, and Seki M. Prediction of Droplet Diameter for Microchannel Emulsification: Prediction Model for Complicated Microchannel Geometries. Industrial & Engineering Chemistry Research 43(26), 8233-8238 (2004). doi: 10.1021/ie0494770.
  » https://doi.org/10.1021/ie0494770
• Sukop M C, Thorne Jr. Lattice Boltzmann Modeling Lattice Boltzmann Modeling. (2006). doi: 10.1007/3-540-27982-2.
  » https://doi.org/10.1007/3-540-27982-2
• Swift MR, Orlandini E, Osborn WR, and Yeomans JM. Lattice Boltzmann Simulations of Liquid-gas and Binary Fluid Systems. Physical Review E 54(5), 5041-5052 (1996). doi: 10.1103/physreve.54.5041.
  » https://doi.org/10.1103/physreve.54.5041
• Thorsen T, Roberts R W, Arnold F H, and Quake S R. Dynamic Pattern Formation in a Vesicle-Generating Microfluidic Device. Physical Review Letters 86(18), 4163-4166 (2001). doi: 10.1103/physrevlett.86.4163.
  » https://doi.org/10.1103/physrevlett.86.4163
• Umbanhowar PB, Prasad V, and Weitz DA. Mono disperse Emulsion Generation via Drop Break Off in a Coflowing Stream. Langmuir 16(2), 347-351(2000). doi: 10.1021/la990101e.
  » https://doi.org/10.1021/la990101e
• Van Der Sman R G, Galilean M. Invariant Lattice Boltzmann Scheme for Natural Convection on Square and Rectangular Lattices. Physical Review E 74(2) (2006). doi: 10.1103/physreve.74.026705.
  » https://doi.org/10.1103/physreve.74.026705
• Van der Graaf S, Nisisako T, Schroen C G P H, Van Der Sman R G M, and Boom R M. Lattice Boltzmann Simulations of Droplet Formation in a T-Shaped Microchannel. Langmuir 22(9), 4144-4152 (2006). doi: 10.1021/la052682f.
  » https://doi.org/10.1021/la052682f
• Waelchli S, and Rohr PR. Two-phase Flow Characteristics in Gas-liquid Micro reactors. International Journal of Multiphase Flow 32(7), 791-806 (2006). doi: 10.1016/j.ijmultiphaseflow.2006.02.014.
  » https://doi.org/10.1016/j.ijmultiphaseflow.2006.02.014
• Wu L, Tsutahara M, Kim LS, and Ha M. Three-dimensional Lattice Boltzmann Simulations of Droplet Formation in a Cross-junction Microchannel. International Journal of Multiphase Flow 34(9), 852-864 (2008). doi: 10.1016/j.ijmultiphaseflow.2008.02.009.
  » https://doi.org/10.1016/j.ijmultiphaseflow.2008.02.009
• Xu JH, Luo GS, Li SW, and Chen G. Shear Force Induced Monodisperse Droplet Formation in a Microfluidic Device by Controlling Wetting Properties. Lab on a Chip 6(1), 131-136 (2006). doi: 10.1039/b509939k.
  » https://doi.org/10.1039/b509939k
• Yasuno M, Sugiura S, Iwamoto S, Nakajima M, Shono A, and Satoh K. Mono dispersed Micro bubble Formation Using Microchannel Technique. AIChE Journal 50(12), 3227-3233 (2004). doi: 10.1002/aic.10276.
  » https://doi.org/10.1002/aic.10276
• Yu Z, Hemminger O, and Fan LS. Experiment and Lattice Boltzmann Simulation of Two-phase Gas-liquid Flows in Microchannels. Chemical Engineering Science 62(24), 7172-7183 (2007). doi: 10.1016/j.ces.2007.08.075.
  » https://doi.org/10.1016/j.ces.2007.08.075
• Zhou H, and Pozrikidis C. The Flow of Suspensions in Channels: Single Files of Drops. Physics of Fluids A: Fluid Dynamics 5, (2): 311(1993). doi: 10.1063/1.858893.
  » https://doi.org/10.1063/1.858893
• Zhou C, Yue P, Feng J, Ollivier-Gooch CF, and Hu H. 3D Phase-field Simulations of Interfacial Dynamics in Newtonian and Viscoelastic Fluids. Journal of Computational Physics 229(2), 498-511(2010). doi: 10.1016/j.jcp.2009.09.039.
  » https://doi.org/10.1016/j.jcp.2009.09.039
• Zou Q, and He X. On Pressure and Velocity Boundary Conditions for the Lattice Boltzmann BGK Model. Physics of Fluids 9, (6): 1591(1997). doi: 10.1063/1.869307.
  » https://doi.org/10.1063/1.869307

Publication Dates

• Publication in this collection
  Jul-Sep 2018

History

• Received
  17 Dec 2016
• Reviewed
  13 May 2017
• Accepted
  30 Sept 2017