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

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.

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 meanfield 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 twophase system is the important feature of the phasefield 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) where ρ = ρ A + ρ B is the total density and the order parameter ϕ = (ρ A + ρ B )/ρ describes the normalized density difference in the two fluids.
Vz 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) The interfacial tension is given by, (3) where ζ is a parameter proportional to the thickness of the interface, defined as, (4) 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): where M is the Cahn-Hilliard mobility and P Cs 2 t zn = + 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 r x zR W and the flow field the lattice speed c is defined by c = δ x /δ t , where δ x is the lattice distance, and δ t is the time step.The speed of sound c s can be related to c by c c 3 s = . τ f and τ g are independent relaxation parameters.We set to minimize numerical errors (Elzanfaly et al., 2015) and w i is the weight factor with w 0 =4/9, w 1-4 =1/9, w 5-8 =1/36.The equilibrium distribution functions are given by Liu and Zhang (2009;2011): where Γ is the tunable parameter regarding the Cahn-Hillard mobility, S X Considering the viscosities of the two phases, the following linear averaged viscosity is used as the mixture viscosity (Van Der Sman et al., 2006): where η A and η B are the viscosities of the two phases.
The local relaxation time τ f can be calculated from the local viscosity: (17) 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) 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 1 3 3 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 ).
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 gi 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 1 , e 5 and e 1 , g 5 and g 8 are unknown after the streaming step.In order to ensure a prescribed order parameter ϕ in at the inlet, these unknown PDFs must satisfy: Assuming that the g i (i=1,5,8) are distributed by their weight factors w i , we have: 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 θ S , the following formula is used to assign the order parameter ϕ S to the solid lattice sites next to the wall (Van der Graaf et al. , 2006).( 23)

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) 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 τ f = 1.0, σ = 0.01 and Γ = 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

S X
Pe UL Ma = g g g g g g g g g in 1 2 3 4 5 6 7 8 0 of the numerical results, which can be written as (Liu and Zhang, 2011): (25) where x 0 and y 0 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.

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 / ø 0 ) and the contact angle.Figure 3 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) 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): 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 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).where R is the initial radius of the droplet, U being the velocity of the moving plate, H is the channel height, ρ

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

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 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(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 breakup.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.
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.
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 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 breakup of the droplet is the result of competition between viscous and capillary forces.
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.
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 14 shows the dependence of the size of the droplet for different flow rate ratios in low Capillary

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.
As shown in Fig. 16, by decreasing the width of the lateral channel, longer droplets are formed inside the microchannel.

CONCLUSIONS
The multiphase lattice Boltzmann method was used to simulate the droplet formation and breakup process   • 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.
each lattice site, respectively.The macroscopic variables are related to the distribution functions by: 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.

Figure 2 .
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.

Figure 4 .
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 6 .
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 .
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 Capillary numbers for two flow rate ratios of 0.3 and 0.5 at fixed viscosity ratio of 0.5

Figure 9 .Figure 10 .
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

Figure 11 .
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 .
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.

Figure 13 .Figure 14 .Figure 15 .Figure 16 .
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 17 .
Figure 17.The effect of width ratio on droplet length for different Capillary numbers, R=0.3 and viscosity ratio of 0.5.