Numerical simulation of collision between two droplets in the T-shaped microchannel with lattice Boltzmann method

In this study, the Lattice Boltzmann Method (LBM) is used to investigate the deformation of two droplets within microfluidic T-junctions (MFTD). In order to increase the accuracy the two immiscible fluids are modeled using the He-Chen-Zhang model. First, this model is applied to ensure that the surface tension effect existing between the droplets and the continuous fluid is properly implemented in the model. Then the collision and merging of the two droplets within the intersection of a T-shaped microchannel is investigated. For generating droplet formation the effects of relevant dimensionless parameters such as the Reynolds, the Weber numbers as well as a collision parameter affecting the two droplets during their motion and deformation are studied. It is found that by increasing the relative velocities of the inlet flows and droplet sizes, the deformation of the two droplets increases significantly. Our results also show that when the surface tension increases, it takes less time for the droplets to co...


List of Symbols Roman Letters
C α [m.s -1 ] Speed of sound D [lu] Initial diameter of the droplet e α [lu ts -1

I. INTRODUCTION
2][3] Microdroplets can be used within these devices to enhance the speed of biological and chemical processes.This is achieved by reducing the volume of the fluid used during the processes such as heating, cooling and separation.Another usage for microdroplets is in multiphase flows within microfluidic systems which have many practical applications such as diffusion time reduction of the species due to the smaller size of the microchannels.One of the challenges in such systems is difficulty in analyzing the motion and deformation of droplets within the microchannels.][6][7] There are some simulation techniques which are used to analyze the two-phase flow cases.Tan et al. 8 simulated the Ferro-fluid droplets within a microchannel with a T-shaped junction.Their work showed that the merging of these droplets can be improved by using the electrical forces to achieve a precise collision of the droplets.The formation of a single droplet for different interfacial tensions and various flow rates through a T-shaped microchannel has also been investigated by Van der Graaf et al.. 9 However, He et al. 10 studied the viscous Rayleigh-Taylor instability problem without considering the surface tension.
Ashgriz and poo 11 performed a two dimensional simulation of two droplets collision.The study was performed in coalescence and reflexive separation regions for a head on collision.They conducted extensive experiments to investigate water droplet collision dynamics with diameter ratios of 1, 0.5, and 0.75; and with Weber numbers and collision parameters ranging from 1-100 and 0-1, respectively.Pan and Suga 12 also reported four major regimes of binary collision: bouncing, coalescence, reflexive separation, and stretching separation.Acevedo-Malavé 13 conducted a study on the three dimensional simulation of off-center collisions of unequaled liquid droplet in a vacuum environment.In this work depending on Weber number two possible outcomes including permanent coalescence and flocculation were observed.
Morozumi et al. 14 also conducted a simulation to differentiate between coalescence and separation of the two droplets after collision.Coalescence of droplets using finite element method was conducted by Mashayek et al. 15 Hsu and Ashgriz 16 simulated the droplet collision phenomenon based on the volume of fluid method (VOF) considering theeffect of the surface tension.An axisymmetric numerical model for the collision of two equal size droplets was developed by Nikolopoulos et al 17 using the VOF method.Their study includes the detailed results about the size of a droplet formed from two colliding droplets.
Droplet motion in confined microchannels has been studied from many different aspects by researchers.For example, Chinnovet et al. 18 investigated a two-phase flow within a short horizontal rectangular channel with the height of 440 µm experimentally.Study of this two phase flows within the microchannel demonstrated an instability in the lateral parts of the channel that had a significant influence on the formation of different two phase flow regimes.In a different study, Liu et al 19 performed an investigation of the droplet formation within a microfluidic cross-junction.They used a numerical technique to estimate the influence of flow rate and viscosity ratio.Moreover, Amaya et al. 20 simulated the bubble formation within MFTD for three distinct regimes.The effects of the inertia, the size and the location of the bubbles on the flow behavior were discussed during the process of bubble formation.Ngo et.al 21 investigated droplet formations based on a two-dimensional numerical model in microfluidic double T-junctions (MFDTD).They predicted that the alternating droplet formation (ADF) region is shifted toward lower capillary numbers when channel width ratio is less than unity.
Other investigations also have studied coalescence of droplets using various numerical methods.Liu et al. 22 implemented a VOF method to track evolution of the gas-liquid free interface with emphasis on the structure of complicated fluid flow during the collision process.Ansari and Nimvari 23 studied the effect of the bubble viscosity on the internal circulation within a bubble using finite volume method and level-set method.Later, Hadidi et al. 24 also performed the effects of different strengths of uniform magnetic fields on the behavior and coalescence of the gas bubbles in a viscous liquid of dielectric fluids.They modeled the coalescence process using level set model to capture the interface of two droplets.Ma et al. 25 numerically investigate the single bubble formation of two coalescence bubbles which is accelerated by increasing superficial liquid velocity.They applied VOF model to characterize the behavior of gas-liquid interface movement.Deen et al. 26 performed a simulation model to predict the heat exchange rate between the liquid and a hot wall kept at a fixed temperature.Zhang et al. 27 analyzed the dynamics and mechanisms of binary collisions of two identical nano-droplets with various impact velocities both in vacuum and an ambient of nitrogen gas using MD simulation.
Studying dynamics of droplets gives useful information about how they behave inside the microchannels.Therefore, many researchers have been recently devoted to develop computational methods to study this problem.The Lattice Boltzmann method (LBM) is one of these computational techniques that its application is growing compared to other conventional simulation schemes.This method is based on mesoscopic kinetic equations and microscopic models that satisfy the Navier-Stokes equation.][35] In some coalescence researches carried out recently droplets are driven by surface tension in disregard of velocity field at the inlets of MFTD.For these cases the surface tension was added to the Navier-Stokes equations as the only force causing the coalescence phenomenon.In other words, small attention is paid on the effect of both the inlet velocity and surface tension simultaneously during the collision of two droplets.The goal of the present study is to examine the parameters that affect the deformation of two droplets during the collision and to investigate the time during which a unified droplet is formed within the MFTD.In this work, the energy budget variation based on kinetic and surface energy distribution of droplet before the collision and consequently during coalescence with respect to non-dimensional time at different surface tension is analyzed.The first section briefly introduces the parameters and the formulations implemented in the LBM.Afterwards, the coalescence collision zone is investigated in order to verify the accuracy of the simulation.Then, the coalescence collision of droplets is systematically studied based on the various parameters such as the droplets size and two inlet velocities of continuous flow and the microchannel fluid-droplet surface tension.By setting different values of these parameters, different positions and shape of unified droplet at the different dimensionless times through the MFTD can be obtained.Finally, in order to evaluate the deformation of droplets, some dimensionless numbers are discussed and their effects have been investigated on the droplets.

II. GOVERNING EQUATIONS FOR THE LBM
Two-dimensional LBM equations for the two-phase flow are implemented.Our primary equation is the Boltzmann equation which may be written for non-ideal fluids as below: 36

Df
In the above equation, f is the single particle distribution function, ξ and u represent the microscopic and macroscopic velocities respectively, f eq stands for the equilibrium distribution function which is the important parameter to characterize the effect of the molecules; what percentage of the molecules in a certain location of a container have velocities within a certain range, at a given instant of time, λ is the relaxation time, G denotes the gravitational force, R is the gas constant, T stands for the temperature, and F represents the effective interaction force that can be written as: where In the above equation, F s denotes the surface tension and the parameter k controls the strength of surface tension.The function ψ is the non-ideal part that accounts for the attractive force between the molecules calculated from the macroscopic pressure by In order to solve Eq. ( 1) without instability, the following transformation is used: where g is the pressure distribution function and Γ is related to the equilibrium distribution function defined as Substituting Eq. ( 5) into Eq.( 1) one gets the following equations: 36,37 It is obvious that parameter g has two additional terms compared to parameter f.In fact, as mentioned above, these two terms are related to the surface and volume forces.φ is the index function and its role is to track the interface between different phases.Using the Carnahan-Starling equation of state, ψ(φ) can be expressed as 38 where parameter a determines the strength of molecular interaction.To ensure phase separation, parameter a in Eq. ( 8) must be chosen to be more than10.601RT.In this research it is assumed that a = 12RT t, similar to the previous studies. 39,40n order to implement an explicit method, the following variable transformations are used: After substituting Eq. (9 and 10) into Eq.( 7), the following equations are obtained: The nine-velocity square lattice models are often referred to as D2Q9 model.Discretization of Eq. ( 7) is performed in a D2Q9 lattice structure, 41 with the microscopic velocities defined as: Where w α is the weighting factor given by: 41 In the LBM, the Courant-Friedrichs-Lewy (CFL) number is proportional to δ t /δ x ; where the discrete space interval δ X is defined as the reference length divided by the number of cells N used to discretize specific length.In the same way, δ t is defined as the reference time divided by the number of iteration steps N needed to reach certain time.Therefore, CFL number can be defined as 42 By choosing δ x = δ t = 1 and |e α | = 1, CFL number is equal to 1.The lattice speed of sound is c 2 s = RT which is constant depending on the LBM.In this work we set c 2 s = RT = 1/3.where Macroscopic variables are evaluated using the following equations: and The fluid density may be obtained by a simple interpolation using the following relation: where the fluid with low density and high density are noted as ρ l and ρ h .Parameters φ 1 and φ h are the minimum and maximum values of the index function, respectively.In this study, φ 1 = 0.022838 and φ h = 0.250291 are used as the limiting values of the index function.The kinematic viscosity is related to the dimensionless relaxation time as

III. NUMERICAL VALIDATION (LAPLACE'S LAW)
For studying the systemconsidered in the present study a visual C++ code was prepared.However, first, we carry out the static bubble test to predict Laplace's law for the surface tension.According to the Laplace's law, when the system reaches the equilibrium state, the pressure difference across the interface of a static droplet is related to its surface tension.That is where p in and p out are the pressures inside and outside of the droplet, respectively.Parameters R is the radius of the droplet at equilibrium, and σ is the surface tension coefficient.
In this case it is assumed that the viscosities are the same for both fluids (υ l = υ h = υ) and second order compact schemes used to evaluate the gradient and Laplacian terms are given as and In order to verify the Laplace's law, static droplets with the different radii are initially generated inside the computational domain with 121×121 lattice points and periodic boundary conditions.The densities of the light and heavy fluids are chosen to be 1 and 5, in turn.Figure 1 shows the pressure difference, ∆P versus 1/R using different values of the parameter k, along with the analytical solution shown by solid lines.The theoretical solution is given by 43 It should be stated that the above integration is calculated along the direction (n) normal to a flat interface.By referring to Fig. 1 it can be seen the present method can simulate and capture the phase interface accurately.

IV. MODEL DESCRIPTION
To study the dynamics of two phase flow, simulations are conducted inside the standard MFTD with two inlet ports and one outlet port as shown in Fig. 2. The width and the length of the horizontal and vertical channels are set to be 20 and 100 lattice cells, respectively.
In this section the effect of gravity force in the simulation of droplets is neglected.Thus, the Z direction will not be considered in the simulation.
The two droplets passing through the channel contain fluid A whereas fluid B has filled with the surrounding droplets.No-slip boundary condition is applied to the whole solid walls of the channel.According to the no-slip boundary condition, the bounce-back boundary condition should be used for the nodes on the solid surfaces. 43The bounce-back boundary condition ensures the conservation of mass and momentum.In D2Q9 model vectors are shown in Fig. 3.
One of the two common methods by which fluid actuation through MFTD achieved is pressure driven flow, in which the fluid is pumped through the device via positive displacement pumps.One of the basic laws of fluid mechanics for pressure driven laminar flow, the so-called no-slip boundary condition, states that the fluid velocity at the walls must be zero.This produces a parabolic velocity profile within the channel.To ensure the accuracy of the computed coalescence regime, the velocity profiles in the different meshes are compared.In Figure 4, the effect of the spatial discretization is shown for the numerical solution in the x = 20 lattice cells along with channel width.It is obvious that in this numerical analyze with increasing the number of meshes from 90×90 to 100×100, as can be seen in figure 4, there is almost no change in the velocity profile of flow streaming in MFTD which indicates the independency of our numerical method from number of meshes.
Figure 5 represents dimensionless distance between centroid of approaching droplet as a function of dimensionless time.Here, distance is normalized by channel width and time is normalized by h U drop where h denotes channel width and U drop is the approaching velocity of droplets.As can be seen in this figure the droplet trajectory is almost linear with time and our results are in a good agreement with experimental research given by Cristopher. 44ere, three cases were considered for the bridge formation at the early stages of the coalescence with ν = 0.2333 ,U ref = 0.2 at density ratio 5. Figure 6   In the above definitions D is the initial diameter of each droplet, U is the relative velocity of the droplets, T is the iteration parameter, ρ is the density, ν is the kinematic viscosity and µ is the dynamic viscosity.Each of the two droplets has an interface with its surrounding fluid.The closest distance between two droplets can be defined as X.The collision parameter is initially maximum.As the droplets move towards each other, this parameter decreases.Collision time is defined as the time when two droplets collide with each other.In the following sections the collision time is obtained for different cases.
Here, we investigate various issues related to the coalescence collision dynamics and merging of the two droplets.Investigation of the coalescence collision of two droplets includes the energy budget of collisions.The kinetic energy of the droplet pair, KE, is defined as: 45   where ρ, D and U are the density, initial diameter and velocity, respectively.Indexes r and l are related to right and left droplets calculated by U r = -U l for identical droplets.Total energy (TE) is defined as a normalized parameter as The surface energy of coalesced collision, SE, is given as

V. RESULTS AND DISCUSSION
A. Effect of different parameters

The effect of size
In this section by changing the droplets size, the dynamics in droplet generation is changed accordingly due to the correlation between two droplets. 46,47Figures 7 and 8 indicate the deformation of droplets with different sizes at various dimensionless times.Here, several parameters including ν = 0.2333, ρ l = 0.1, ρ h = 0.5, k = 0.15, U l = 0.1 and U r = -0.1 remain constant.As shown in Fig. 7, the first case indicates the deformation of two droplets with a size of 6 × 6 lattice cells.The sequence of merging droplets is also generated for the second case revealing the deformation of the two droplets with a size of 8 × 8 lattice cells as represented in Fig. 8. Reynolds number takes the values of 5.14 and 6.86 while Weber number increases from 60 to 80 in the Figs.7 and 8, respectively due to the fact that the Reynolds and Weber Number is proportional to droplets size.In addition, with the variation of droplets size Ohnesorge number decreases from 1.51 to 1.3.Comparing Figs.7.f and 8.f, it can be seen by increasing the Weber and Reynolds numbers, the deformation of the unified droplet is increased remarkably.
As to collision time, coalescence occurs at the dimensionless time rate of 0.73 by increasing size from 6 lattice cells to 8 lattice cells.Therefore, increasing the size of the droplets causes a decrease in their collision time values.Due to pressure rise between larger droplets, drainage of continuous flow occurs in the shorter time.Figure 9 displays the pressure variation between approaching droplets versus dimensionless time.As they approach together, the pressure field on the centerline between droplets increases.It has been seen that there is a higher pressure region between the large droplets compared to smaller ones.It can also be observed that as the droplets size becomes larger, the shape of the unified droplet departs from circle.This is mainly because increasing the size of droplets leads to larger surface deformation due to rise in Reynolds.
FIG. 9. Variation of average pressure between centroid of two droplets via non-dimensional time.

The effect of inlet velocity
It should be mentioned that in the many relevant researchers conducted previously there have been many restrictions.For instance, in the coalescence of two drops only the effect of the surface tension has been considered while in our research we also consider the effect of the velocity field on the deformation of these two droplets.In this case, numerical approaches can play an important role in studying the interfacial dynamics of fluid flows in small channels, where the evolution of droplet formation is monitored.To consider the effects of the continuous inlet velocity on the droplet formation, we analyze two cases by considering constant parameters such as: size = 6×6 lattice cells, ρ l = 0.1, ρ h = 0.5, ν = 0.2333, , and k = 0.15.In the first case, inlet velocities on the left and right sides of the channel are U l = 0.11 , U r = -0.1 respectively (Fig. 10).In the second case (Fig. 11), the left and the right velocities are set to the value of 0.12 and -0.1, respectively.In Figs. 10 and 11 Weber numbers are equal to 66.15 and 72.6, respectively.In these cases, Reynolds number also rises from 5.4 to 5.66.It was observed that when the difference between the inlet velocities is increased, the Weber number increases significantly compared to the Reynolds number due to variation of Weber number with quadratic relationship of velocity.Although, In Fig. 8.e, the resulting unified droplet is completely symmetric, In Fig. 10.e, at the dimensionless time t * = 24.5 there is no symmetry in the shape of the unified droplet.For high Weber number droplets considerably deform during the collision and turn to an asymmetric shape since the inertial force is small and surface tension force is dominant.
As the left inlet velocity rises from 0.1 to 0.11(Fig.10), the collision time reaches around 20.02.In this case, it is obvious that the less time is needed for the merging droplets to form a unified droplet compared to the case in which the inlet velocities are equal (Fig. 8).This is principally because when the velocities are low, the time it takes to drain the film is likely to be long.In the second case (Fig. 11), since the difference between the inlet velocities is bigger, the left droplet reaches the intersection of the channels more quickly, so there is no collision between the droplets within the microchannel.Accordingly, the collision time cannot be estimated for this case.Finally, each of the two droplets continues its motion in the vertical section of the channel.As the velocity of the left droplet exceeds over the value of 0.117, these droplets do not collide.A critical Reynolds number of 5.58 can be observed for this condition.When the Reynolds number is more than this critical Reynolds number, there is no contact between the two droplets.In the case of k equal zero, at the dimensionless time of 20, two droplets entering the vertical channel are more stretched along the junction.Comparing Figure 12 to Figure 13, in the case of larger surface tension, two droplets are more accelerated to each other and cause to drainage of gap between two droplets in shorter time.

The effect of the surface tension
Also, decreasing Weber number leads to a less deformation of the merging droplets (Fig. 13).By comparing Figs. 12 and 13 at the t * = 26.6, it is seen that by increasing parameter k, the combined droplet becomes more circular in shape.This is mainly because increasing the surface tension keeps droplets from propagation and hold them in the circular shape in which droplets have minimum potential of energy.

B. Energy budget
This section highlights the energy of two droplets during approaching together.As can be seen in Fig. 14 the less the distance of droplets is, the less the velocity of droplets becomes.An examination of Fig. 14  case with k = 0.In Fig. 15, the normalized kinetic, surface and total energy versus dimensionless time are given.The surface energy can be obtained by subtraction of the kinetic energy from the total energy.It is obvious that the kinetic energy and the surface energy graphs behave exactly opposite to each other.When droplets collide to each other, the kinetic energy attains its minimum value which will be larger for the case with k = 0.35 compared to that of case with k = 0.Then, the velocity of the combined droplet starts increasing through the vertical channel and, therefore, kinetic energy increases.Also, it has been proved that after merging of two droplets, the coalesced droplet firstly takes a non-equilibrium shape and then relaxes towards a circular configuration in order to minimize the surface energy.

VI. CONCLUSION
The coalescence of two droplets through the MFTD was studied using the Lattice Boltzman Method.The two droplets are driven by the velocity field of the continuous phase.Here, discussed AIP Advances 6, 115307 (2016) parameters were droplets size, inlet velocity and surface tension existing between two droplets and continuous flow.By changing each of the above parameters, coalescence of droplets were analyzed at different dimensionless times and it was found that deformation of these droplets varies significantly depending on these parameters.As the size of two droplets increases, it was observed that formation of the merged droplet occurs in a shorter time.This leads to further deformation of the droplets entering the vertical channel.Increasing the left inlet velocity (Fig. 7) results in further changes in the momentum of the left droplet that leads to an asymmetric unified droplet.When the inlet velocity increases beyond a critical value, there is no collision between two droplets (Fig. 11).Also, when the surface tension increases, the pressure difference between the inside and outside of the droplet increases.In this case, there is less resistance against merging of the two droplets; thus, coalescence process occurs faster (Fig. 13).It was also found that as the surface tension increases the shape of the unified droplet becomes more circular.By increasing the Weber and Reynolds numbers or decreasing the Ohnesorge number, deformation of two droplets increases.As the Reynolds and Weber numbers increase, the inertial forces prevail over the viscous and the surface tension forces.Once the inertial forces are directly proportional to the momentum of the droplets, an increase in these forces means more momentum leading to more deformation.Finally, kinetic and surface energies for different surface tension values were studied.

FIG. 2 .
FIG. 2. Problem description: Initial position of fluid A and continuous fluid B.

FIG. 7 .
FIG.4.Velocity profile in MFTD under conditions of pressure driven flow.