Lattice-Boltzmann Simulation and Experimental Validation of a Microfluidic T-Junction for Slug Flow Generation

We investigate the interaction of two immiscible fluids in a head-on device geometry, where both fluids are streaming opposite to each other. The simulations are based on the two-dimensional (2D) lattice Boltzmann method (LBM) using the Rothman and Keller (RK) model. We validate the LBM code with several benchmarks such as the bubble test, static contact angle, and layered flow. For the first time, we simulate a head-on device by forcing periodicity and a volume force to induce the flow. From low to high flow rates, three main flow patterns are observed in the head-on device, which are dripping-squeezing, jetting-shearing, and threading. In the squeezing regime, the flow is steady and the droplets are equal. The jetting-shearing flow is not as stable as dripping-squeezing. Moreover, the formation of droplets is shifted downstream into the main channel. The last flow form is threading, in which the immiscible fluids flow parallel downstream to the outlet. In contrast to other studies, we select larger microfluidic channels with 1-mm channel width to achieve relatively high volumetric fluxes as used in chemical synthesis reactors. Consequently, the capillary number of the flow regimes is smaller than 10−5. In conclusion, the simulation compares well to experimental data.


Introduction
In recent years, microfluidics became an attractive field of research, since it has a wide range of applications such as lab-on-a-chip devices, cosmetics, drug delivery, microfabrication, and chemical synthesis [1,2].In continuous-flow microfluidics, droplet flow can be used to act as an isolated reactor with the minimum consumption of reagents [3,4].Therefore, droplet formation in liquid/liquid or gas/liquid systems is an important area of study.One of the simplest but most common methods to generate droplets is the use of appropriate channel geometries [5][6][7].The T-junction is the simplest microfluidic geometry for producing droplets.
The scope of this study is droplet formation in a 1-mm channel geometry.In this range, droplet flow is successfully used in liquid/liquid systems for extraction [8][9][10] and synthesis of fine chemicals and pharmaceuticals [11][12][13][14].In gas/liquid systems, applications include membrane fuel cells [15] and chemical synthesis [16][17][18].For liquid-liquid extraction (LLE) processes, flow chemistry offers the benefit of excellent control over interfacial exchange parameters.In combination with chemical reactions in flow, LLE can be integrated as an in-line purification technique.Due to the combination of process steps, faster exchange of matter, and reduction in chemical consumption, LLE in microsystems leads to highly efficient processes [8][9][10].For chemical synthesis, the short diffusion times for mass and heat transfer lead to a new field of process intensification, which includes high temperatures, high pressures, and high concentrations (solvent-free) to explosive conditions [11][12][13][14].To achieve highly productive microfluidic systems, a parallelization of the microchannel is crucial [19].Furthermore, to build precise, uniform flow distributors, a good knowledge of the droplet formation process is necessary.
One configuration of a T-junction is a so-called head-on device.It consists of three channels to mix two immiscible fluids.In this device, the dispersed phase (DP) and the continuous phase (CP) flow toward each other, as one can see from Figure 1.The dispersed phase should not wet the walls.Shear stress and pressure drop play an important role in the droplet generation.The continuous phase creates a shear stress which detaches the dispersed phase and generates a droplet.The droplet length and frequency in T-junctions depends on various parameters such as flow rate, geometry, viscosity, surface tension, and contact angle.Over the last few years, numerous studies on the experimental and numerical level were performed to determine the effect of each parameter on the droplet formation.In the following sections, we discuss these parameters in detail.
ChemEngineering 2019, 3, x FOR PEER REVIEW 2 of 14 of process steps, faster exchange of matter, and reduction in chemical consumption, LLE in microsystems leads to highly efficient processes [8][9][10].For chemical synthesis, the short diffusion times for mass and heat transfer lead to a new field of process intensification, which includes high temperatures, high pressures, and high concentrations (solvent-free) to explosive conditions [11][12][13][14].
To achieve highly productive microfluidic systems, a parallelization of the microchannel is crucial [19].Furthermore, to build precise, uniform flow distributors, a good knowledge of the droplet formation process is necessary.
One configuration of a T-junction is a so-called head-on device.It consists of three channels to mix two immiscible fluids.In this device, the dispersed phase (DP) and the continuous phase (CP) flow toward each other, as one can see from Figure 1.The dispersed phase should not wet the walls.Shear stress and pressure drop play an important role in the droplet generation.The continuous phase creates a shear stress which detaches the dispersed phase and generates a droplet.The droplet length and frequency in T-junctions depends on various parameters such as flow rate, geometry, viscosity, surface tension, and contact angle.Over the last few years, numerous studies on the experimental and numerical level were performed to determine the effect of each parameter on the droplet formation.In the following sections, we discuss these parameters in detail.

Geometry
Shi et al. [20] studied numerically the droplet formation in a T-junction using a lattice Boltzmann (LB) free-energy model.They concluded that the ratio of the main channel width to the lateral channel width of the T-junction plays an important role in the droplet formation.They found that the droplet length increases for a fixed capillary number with an increase in the ratio of dispersed to continuous phase inlet.

Contact Angle
Shazia Bashir et al. [1] found that the surface properties of the main channel in T-junctions determine the droplet length.The contact angle describes the interaction between the fluid and channel wall.They found that the droplet break-up time decreases as the contact angle increases from 120 to 180 degrees.Moreover, the fluid velocity in the microchannel increases with incrementing contact angle.Finally, when the wetting property of the channel wall is hydrophilic, the advancing and receding contact angle changes [20].

Flow Rate
Haihu Liu et al. [21] studied the influence of the flow rate ratio on droplet formation.They found in both simulation and experiment that the droplet diameter decreases gradually upon increasing capillary number.They separated the capillary number using a critical value into large and small capillary numbers.Droplet size variation with increase in contact angle is higher in small capillary numbers; however, it decreases in high capillary numbers [21,22].Gupta et al. [23] found that the droplet volume grows when the dispersed flow rate increases.Moreover, a linear relationship

Geometry
Shi et al. [20] studied numerically the droplet formation in a T-junction using a lattice Boltzmann (LB) free-energy model.They concluded that the ratio of the main channel width to the lateral channel width of the T-junction plays an important role in the droplet formation.They found that the droplet length increases for a fixed capillary number with an increase in the ratio of dispersed to continuous phase inlet.

Contact Angle
Shazia Bashir et al. [1] found that the surface properties of the main channel in T-junctions determine the droplet length.The contact angle describes the interaction between the fluid and channel wall.They found that the droplet break-up time decreases as the contact angle increases from 120 to 180 degrees.Moreover, the fluid velocity in the microchannel increases with incrementing contact angle.Finally, when the wetting property of the channel wall is hydrophilic, the advancing and receding contact angle changes [20].

Flow Rate
Haihu Liu et al. [21] studied the influence of the flow rate ratio on droplet formation.They found in both simulation and experiment that the droplet diameter decreases gradually upon increasing capillary number.They separated the capillary number using a critical value into large and small capillary numbers.Droplet size variation with increase in contact angle is higher in small capillary numbers; however, it decreases in high capillary numbers [21,22].Gupta et al. [23] found that the droplet volume grows when the dispersed flow rate increases.Moreover, a linear relationship between the droplet volume and flow rate for low capillary number was reported by Demench et al. [22].

Viscosity
The viscosity of both the continuous phase and dispersed phase influences the droplet diameter.For a fixed continuous phase viscosity, an increase in the dispersed phase viscosity results in a decrease in droplet length [24].Similarly, the droplet volume grows as the viscosity increases [23].

Surface Tension
The droplet length increases for higher surface tension values and viscosity ratios of the fluids [1].In this paper, we present the simulation result of the droplet formation in a head-on device using a two-phase LB method (LBM).We compare the numerical results to experimental data.The size of the microfluidic channels of our system is 1 mm to reach relatively high volumetric fluxes as used in chemical synthesis reactors.Thus, the lateral dimension of our microchannel is large compared to the majority of other publications in this field [25].As a consequence, the capillary number of the flow regimes is smaller than 10 −5 , which is approximately two orders of magnitude smaller than other studies.Additionally, we use T-junctions made from poly(methyl methacrylate) (PMMA) and a water/heptane system which has no clear distinction between the dispersed and the continuous phases due to a contact angle of approximately 90 degrees.For this system, we show the influence of the volumetric flux on the droplet formation in the head-on device.

Multiphase LBM
The simulations presented here are based on the LBM using the Rothman and Keller (RK) model.The LBM multiphase method in this study is based on publications of Liu and Leclaire et al. [26][27][28].They presented a detailed derivation of the model in their publication.We decided to use the RK model, also called the color gradient model.The advantages of this method are that the surface tension, viscosity, and contact angle are independently adjustable.Moreover, this method separates the phases with a clear phase interface.

Rothman and Keller or Color Gradient Model
For two-phase flow, two distribution functions are defined.The macroscopic properties are defined using the distribution functions.As an example, the density can be written as: The velocity is calculated using the density and average weight of discrete velocities as: Finally, the pressure for the D2Q9 is calculated as: For further information about the parameters, please refer to Reference [29].
The evolution of the distribution function in time for each phase k can be expressed as: The collision Ω k i is calculated in three steps.

Single-Phase Collision
This equation is for the relaxation in the single-phase system.It assures that the system obeys the macroscopic values of the Navier-Stokes equations.
where τ k is the relaxation time in the collision process and f k,eq i is the local equilibrium distribution of the fluid velocities.τ k is related to kinematic viscosity by: In this project, we used a single relaxation time for both phases.The density average of kinematic viscosity is defined as: The local equilibrium f k,eq i is determined by the Maxwell-Boltzmann distribution.
where φ k i adjusts the compressibility of the fluid and w i is the weighting parameter for each node.

Collision or Perturbation
Operator The perturbation operator reproduces the surface tension in the model.To separate the phases, a color field ψ(x, t) can be defined.
Reis and Phillips [30] introduced: The interfacial properties adjusted by A k and B i can recover the Navier-Stokes equations.
2.5.Two-Phase Collision (Recoloring) The recoloring separates two fluids from each other.The RK model provides a clear interface between the phases in comparison to other methods.
ChemEngineering 2019, 3, 48 5 of 14 where ρ is the total density, f i is the total particle distribution function, and ϕ i is the angle between the color gradient and the lattice direction.f i and ϕ i take values between 0 and 1 [31,32].
The operator allows mixing at the interface between the phases.The interfacial thickness can be adjusted by β, but the interfacial tension is determined by A k .

Model Implementation and Hardware Requirements
The described LB model was implemented in Matlab (Version R2017a).The required memory was very moderate.Typically, a lattice of 512 × 512 = 262,144 nodes requires about 20 Mb of random-access memory (RAM).The bottleneck of the transient simulations was given by the time steps.As pointed out in the previous section in Equation (6), the time step ∆t is closely related to the parameters of the simulated fluids.Therefore, the generation of one flow pattern, as shown in Section 5, takes several hours on a modern 3-GHz central processing unit (CPU).The exact values of CPU time depend on the chosen volumetric flux.

Model validation
In order to verify our model, we performed several benchmark tests.

Laplace Test
There is a force balance at the interface between the two phases.The change in the interface form results in the pressure difference between the fluids.Surface tension can be described using the Laplace law [33].
where σ is the surface tension and ∆p is the pressure difference.This equation represents the linear relationship between surface tension and curvature 1/r.The surface tension can be estimated by simulating a series of drops with various radii (see Figure 2).
ChemEngineering 2019, 3, x FOR PEER REVIEW 5 of 14 The operator allows mixing at the interface between the phases.The interfacial thickness can be adjusted by , but the interfacial tension is determined by  .

Model Implementation and Hardware Requirements
The described LB model was implemented in Matlab (Version R2017a).The required memory was very moderate.Typically, a lattice of 512 × 512 = 262,144 nodes requires about 20 Mb of randomaccess memory (RAM).The bottleneck of the transient simulations was given by the time steps.As pointed out in the previous section in Equation (6), the time step ∆t is closely related to the parameters of the simulated fluids.Therefore, the generation of one flow pattern, as shown in Section 5, takes several hours on a modern 3-GHz central processing unit (CPU).The exact values of CPU time depend on the chosen volumetric flux.

Model validation
In order to verify our model, we performed several benchmark tests.

Laplace Test
There is a force balance at the interface between the two phases.The change in the interface form results in the pressure difference between the fluids.Surface tension can be described using the Laplace law [33].
where  is the surface tension and ∆ is the pressure difference.This equation represents the linear relationship between surface tension and curvature 1/r.The surface tension can be estimated by simulating a series of drops with various radii (see Figure 2).

Static Contact Angle
To validate the contact angle, it should be implemented in a capillary free model.Therefore, the system in Figure 3 can be used to assign two contact angles for each surface.The model size was 32 × 200 lattice units (lu) and, in the x-direction, a periodic boundary condition was applied.Figure 3 was calculated using  = 1 and viscosity = 1/6.In each simulation, a contact angle was assigned to the top and bottom (for instance, 85 and 95 degrees).

Static Contact Angle
To validate the contact angle, it should be implemented in a capillary free model.Therefore, the system in Figure 3 can be used to assign two contact angles for each surface.The model size was ChemEngineering 2019, 3, 48 6 of 14 32 × 200 lattice units (lu) and, in the x-direction, a periodic boundary condition was applied.Figure 3 was calculated using ρ = 1 and viscosity = 1/6.In each simulation, a contact angle was assigned to the top and bottom (for instance, 85 and 95 degrees).

Layered Flow for Immiscible Two-Phase Flow
In the Laplace test, we verified the surface tension, and the fluid interaction with the surface was verified with the contact angle.These methods are mostly steady; however, the fluid behavior in Tjunctions is dynamic.The last validation method verified the numerical method in a dynamic test.The analytical solution for the layered flow was calculated from the Navier-Stokes equation [34].
where is the pressure gradient, ℎ is the channel radius,  is the position in the channel, and  and  are the kinematic viscosity for red and blue phases.We simulated the layered flow for a 40-lu channel as shown in Figure 4.The viscosity ratio was chosen as 1:5 at a Reynolds numbers of Re = 25 with respect to the higher viscous fluid and Re = 125.At the solid walls, the half-way bounce-back boundary was used to reproduce the non-slip boundary condition.

Layered Flow for Immiscible Two-Phase Flow
In the Laplace test, we verified the surface tension, and the fluid interaction with the surface was verified with the contact angle.These methods are mostly steady; however, the fluid behavior in T-junctions is dynamic.The last validation method verified the numerical method in a dynamic test.The analytical solution for the layered flow was calculated from the Navier-Stokes equation [34].
where dp dx is the pressure gradient, h is the channel radius, x is the position in the channel, and µ b and µ r are the kinematic viscosity for red and blue phases.
We simulated the layered flow for a 40-lu channel as shown in Figure 4.The viscosity ratio was chosen as 1:5 at a Reynolds numbers of Re = 25 with respect to the higher viscous fluid and Re = 125.At the solid walls, the half-way bounce-back boundary was used to reproduce the non-slip boundary condition.
and  are the kinematic viscosity for red and blue phases.
We simulated the layered flow for a 40-lu channel as shown in Figure 4.The viscosity ratio was chosen as 1:5 at a Reynolds numbers of Re = 25 with respect to the higher viscous fluid and Re = 125.At the solid walls, the half-way bounce-back boundary was used to reproduce the non-slip boundary condition.

Simulation of the Droplet Formation in the Head-On Device
There are different methods to stimulate the flow in LBM.Velocity and pressure boundary conditions are two possible approaches.These boundary conditions operate by varying the densities at specified points.Another method is to initiate the flow by volume or external force.This method is convenient to implement especially in combination with periodic boundary conditions.In this study, all the mentioned boundary conditions were applied.In the following sections, we discuss each method in detail.

Velocity-Pressure Boundary Conditions
The head-on device geometry is shown in Figure 5a.The velocity boundary condition was defined at the inlets, and the pressure boundary condition was defined at the outflow.These boundaries were implemented based on Sukop et al. [35].Unfortunately, the application of these boundary conditions in LBM led to several problems.At first, reproducing the flow rates of the experiment at a Reynolds numbers of approximately 10 resulted in an unrealistically high pressure gradient in the main channel.Consequently, the droplet size expanded gradually in the main channel.In addition, at low flow rates, these boundary conditions corresponded to a compressible behavior of the fluids.Finally, the velocity boundary condition did not satisfy the conservation of mass in the system, which was also reported in References [35,36].To overcome these drawbacks, we used a volume force in combination with periodic boundary conditions as explained below.

Simulation of the Droplet Formation in the Head-On Device
There are different methods to stimulate the flow in LBM.Velocity and pressure boundary conditions are two possible approaches.These boundary conditions operate by varying the densities at specified points.Another method is to initiate the flow by volume or external force.This method is convenient to implement especially in combination with periodic boundary conditions.In this study, all the mentioned boundary conditions were applied.In the following sections, we discuss each method in detail.

Velocity-Pressure Boundary Conditions
The head-on device geometry is shown in Figure 5a.The velocity boundary condition was defined at the inlets, and the pressure boundary condition was defined at the outflow.These boundaries were implemented based on Sukop et al. [35].Unfortunately, the application of these boundary conditions in LBM led to several problems.At first, reproducing the flow rates of the experiment at a Reynolds numbers of approximately 10 resulted in an unrealistically high pressure gradient in the main channel.Consequently, the droplet size expanded gradually in the main channel.In addition, at low flow rates, these boundary conditions corresponded to a compressible behavior of the fluids.Finally, the velocity boundary condition did not satisfy the conservation of mass in the system, which was also reported in References [35,36].To overcome these drawbacks, we used a volume force in combination with periodic boundary conditions as explained below.

Periodic Boundary Conditions and Volume Force
As mentioned before, external or volume force is another method to initiate the flow field.The external force can be set for each fluid separately in the domain.We combined the volume force with periodic boundary conditions, which required a modified geometry as depicted in Figure 5b.As a result of the periodic boundary condition, we guaranteed the conservation of mass for both phases exactly [35].
As shown in Table 1, we wanted to simulate the same flow rate for both phases.This required a proportional integral (PI) controller.We placed sensors 1 and 2 close to the inlet and sensor 3 on the main channel.Sensors 1 and 2 controlled the flow for an equal flow rate simultaneously, while sensor

Periodic Boundary Conditions and Volume Force
As mentioned before, external or volume force is another method to initiate the flow field.The external force can be set for each fluid separately in the domain.We combined the volume force with periodic boundary conditions, which required a modified geometry as depicted in Figure 5b.As a result of the periodic boundary condition, we guaranteed the conservation of mass for both phases exactly [35].
As shown in Table 1, we wanted to simulate the same flow rate for both phases.This required a proportional integral (PI) controller.We placed sensors 1 and 2 close to the inlet and sensor 3 on the main channel.Sensors 1 and 2 controlled the flow for an equal flow rate simultaneously, while sensor 3 was used to maintain the defined flow rate in the main channel (Figure 6).The values of the sensors were constantly monitored.Based on these values, the controller changed the external force to achieve the desired flow rates.conditions.

Periodic Boundary Conditions and Volume Force
As mentioned before, external or volume force is another method to initiate the flow field.The external force can be set for each fluid separately in the domain.We combined the volume force with periodic boundary conditions, which required a modified geometry as depicted in Figure 5b.As a result of the periodic boundary condition, we guaranteed the conservation of mass for both phases exactly [35].
As shown in Table 1, we wanted to simulate the same flow rate for both phases.This required a proportional integral (PI) controller.We placed sensors 1 and 2 close to the inlet and sensor 3 on the main channel.Sensors 1 and 2 controlled the flow for an equal flow rate simultaneously, while sensor 3 was used to maintain the defined flow rate in the main channel (Figure 6).The values of the sensors were constantly monitored.Based on these values, the controller changed the external force to achieve the desired flow rates.

Flow Conditions
In addition to the described boundary conditions, we chose flow conditions as close as possible to the experiment described in Section 4. The simulation was performed on three different grid resolutions, i.e., channel widths of 10, 20, and 40 lu corresponding to a channel width of 1 mm.Additionally, we could match exactly the viscosity ratio of heptane/water with the experiment.There was a small difference between the contact angle in the simulation (80 • ) and the experiment (90 • ) in order to increase the numerical stability of the simulation.Moreover, we could perform stable simulations with a Reynolds number of 1-17, which was comparable to the range in the experiment.Finally, it was not possible to match the capillary number between the experiment and simulation.This was due to the relatively large channel width.The simulation of capillary numbers down to 10 −5 would require a grid resolution one hundred times higher and was not feasible.Nevertheless, the simulation results show that the simulation and the experiment were comparable, since both two-phase systems were in a similar regime.All flow parameters for simulation and experiment can be found in Table 1. a the 1-mm channel was simulated using various channel sizes.

Droplet Formation Process
Thorsen et al. [37] first introduced droplet generation using T-junctions, and it was followed by Shui et al. [38] in head-on devices.In our simulations and in the experiments, we observed five steps of droplet generation.This is depicted in Figure 7 as follows: (1) the disperse phase enters the junction; (2) it grows until the entire width of channel is filled; (3) the continuous phase squeezes the disperse phase at the junction; (4) the disperse phase begins to break at the junction; (5) finally the droplet dispatches into the continuous phase.These steps correspond to the data published by Shui [39]. of droplet generation.This is depicted in Figure 7 as follows: (1) the disperse phase enters the junction; (2) it grows until the entire width of channel is filled; (3) the continuous phase squeezes the disperse phase at the junction; (4) the disperse phase begins to break at the junction; (5) finally the droplet dispatches into the continuous phase.These steps correspond to the data published by Shui [39].Additionally, Shui et al. [39] reported three different flow regimes at the outlet of the head-on device.In their publication, they called these flow regimes as follows: (I) dripping-squeezing, (ii) Additionally, Shui et al. [39] reported three different flow regimes at the outlet of the head-on device.In their publication, they called these flow regimes as follows: (I) dripping-squeezing, (ii) jetting-shearing, and (iii) threading.In this study, we were able to reproduce the three flow types using LBM.In the following sections, each flow form is discussed in detail.

Dripping-Squeezing
In this regime, the droplets are mostly homogeneous and equal (Figure 8).Here, the surface tension forces are greater than shear forces.The flow is stable, and droplets form rapidly.The homogeneous droplets have practical use in microfluidics, as these droplets act as individual reactors.
ChemEngineering 2019, 3, x FOR PEER REVIEW 9 of 14 jetting-shearing, and (iii) threading.In this study, we were able to reproduce the three flow types using LBM.In the following sections, each flow form is discussed in detail.

Dripping-Squeezing
In this regime, the droplets are mostly homogeneous and equal (Figure 8).Here, the surface tension forces are greater than shear forces.The flow is stable, and droplets form rapidly.The homogeneous droplets have practical use in microfluidics, as these droplets act as individual reactors.

Jetting-Shearing
This range initiates with a transition from homogeneous droplets to varying droplet sizes.The increase in flow rate results in a new flow form.In this range, the viscosity is dominant at the junction and, therefore, the tip of the disperse flow enters along the continuous phase in the main channel and forms a parallel/layered flow (steps 1-3 in Fig. 9).The surface tension force creates a round form at the tip of the flow (step 4).The tip advances prior to the wall surface.As it contacts the wall, the continuous phase squeezes the disperse phase and a droplet detaches (step 5).The droplet generation cycle progresses and a continuous phase droplet forms (step 6).At this point in the simulation, this process ends as the main channel is filled with layered flow and the droplets detach downstream.In the experiment, we observed that the flow changes its behavior in a cycle as follows: the droplets are generated rapidly with various sizes; afterward, the flow changes to a layered flow and the cycle continues.

Jetting-Shearing
This range initiates with a transition from homogeneous droplets to varying droplet sizes.The increase in flow rate results in a new flow form.In this range, the viscosity is dominant at the junction and, therefore, the tip of the disperse flow enters along the continuous phase in the main channel and forms a parallel/layered flow (steps 1-3 in Figure 9).The surface tension force creates a round form at the tip of the flow (step 4).The tip advances prior to the wall surface.As it contacts the wall, the continuous phase squeezes the disperse phase and a droplet detaches (step 5).The droplet generation cycle progresses and a continuous phase droplet forms (step 6).At this point in the simulation, this process ends as the main channel is filled with layered flow and the droplets detach downstream.In the experiment, we observed that the flow changes its behavior in a cycle as follows: the droplets are generated rapidly with various sizes; afterward, the flow changes to a layered flow and the cycle continues.
cycle progresses and a continuous phase droplet forms (step 6).At this point in the simulation, this process ends as the main channel is filled with layered flow and the droplets detach downstream.In the experiment, we observed that the flow changes its behavior in a cycle as follows: the droplets are generated rapidly with various sizes; afterward, the flow changes to a layered flow and the cycle continues.

Threading
As the flow rate increases, the stable droplet flow turns into unstable flow with various droplet sizes.In the experiment, the droplets are generated rapidly but they are non-uniform.At this point, by increasing the flow rate, the viscous forces are dominant with neither surface tension nor shear stresses existing.At higher flow rates, the pressure drop is eliminated at the junction; therefore, the tip of the dispersed flow does not have the time to grow and block the channel and shape the layered flow [39].However, in this regime, the droplets are generated downstream of the T-junction.We observed this phenomenon in the experiment and simulation, as depicted in Figure 9.The transient flow was also reported by References [21,23,39].
Experiment and simulation showed a transition between jetting-shearing and threading at comparable Reynolds numbers of Re > 9 in the simulation and Re > 12 in the experiment.

Experimental Data
The flow experiments were carried out in a self-built micro-structured head-on device made from poly (methyl methacrylate) (PMMA) (Figure 10).The design consisted of a square cross-section of 1 × 1 mm which was manufactured by mechanical milling (MD24 CNC, GOLmatic Werkzeugmaschinen, Birkenau, Germany).The head-on device was cautiously bonded with a second PMMA plate (100 × 100 × 15 mm) using acetone.The second PMMA plate was made to be slightly thicker due to the fluid connections (1/4"-28 UNF-Flat-Bottom port) for the inlet and outlet streams.The inlet and outlet streams were maintained with a 1/16" FEP capillary (fluorinated ethylene propylene).The continuous feed stream of the liquid reactants was provided by a syringe pump (Pump 11 Elite, Harvard Apparatus, Cambridge, MA, USA) equipped with two 25-mL glass syringes (ILS GmbH, Stützerbach, Germany).The outlet stream was collected in a headspace vial.For the determination of the two-phase flow characteristics, the liquid stream was imaged with a digital camera (modelEOS 60D EF-S 15-85 mm, Canon, Osaka, Japan).The T-junction was illuminated with a light source (ACE Light Source, Polytec GmbH, Waldbronn, Germany) for improved visualization.The experimental set-up is displayed in Figure 10.
(ILS GmbH, Stützerbach, Germany).The outlet stream was collected in a headspace vial.For the determination of the two-phase flow characteristics, the liquid stream was imaged with a digital camera (modelEOS 60D EF-S 15-85 mm, Canon, Osaka, Japan).The T-junction was illuminated with a light source (ACE Light Source, Polytec GmbH, Waldbronn, Germany) for improved visualization.The experimental set-up is displayed in Figure 10.The droplet formation process was studied in a head-on device with heptane (VWR, United States of America) which was used without further purification and deionized water.No surfactants were used in this study to stabilize the flow pattern.We used water and heptane since they have similar densities.The water in contact with PMMA was hydrophobic and the heptane was hydrophilic.The experiments were carried out at room temperature and atmospheric pressure.The operating parameters of the flow experiments including channel dimensions, flow rates, and characteristic dimensionless numbers are displayed in Table 1.The droplet formation process was studied in a head-on device with heptane (VWR, United States of America) which was used without further purification and deionized water.No surfactants were used in this study to stabilize the flow pattern.We used water and heptane since they have similar densities.The water in contact with PMMA was hydrophobic and the heptane was hydrophilic.The experiments were carried out at room temperature and atmospheric pressure.The operating parameters of the flow experiments including channel dimensions, flow rates, and characteristic dimensionless numbers are displayed in Table 1.

Quantitative Comparison between Simulation and Experimental Data
In addition to the qualitative agreement between our LBM simulations and the flow experiments, we discuss here the size of the generated droplets as a quantitative measure.However, we must note that we are comparing two-dimensional (2D) simulations to three-dimensional (3D) experimental results.
As one can see from Figure 11, the length of the droplets in the experiments was more or less constant over the investigated range of flow rates.Additionally, the droplets of the continuous phase (water) were slightly larger than those of the dispersed phase (heptane).The LBM simulation results were comparable but showed some remarkable differences.
We performed a range of simulations to study the droplet diameter's relationship with the flow rate.In the simulations, a Reynolds number larger than nine led to the jetting-shearing regime, while, in the experiments, the droplet formation was stable up to a Reynolds number of approximately 12.
As seen in Figure 11, the droplet lengths between the simulations and experiments were comparable.Interestingly, both methods showed larger droplets for the continuous phase (water) and smaller ones for the dispersed phase (heptane).This was independent of the flow rate and, obviously, a result of the difference in the viscosity of both fluids.In the simulations, smaller flow rates corresponded to an increased variance in the droplet length.This variance was directly linked to our simulation approach where we used PI controllers to generate constant volumetric fluxes.Additionally, we saw a similar effect when we changed the lattice resolution, whereby we got larger variances in droplet length with higher resolution.
In summary, we found a good agreement between the simulations and the experimental data.To our belief, the observed differences were mainly an effect of the simulation being reduced to 2D and that the boundary conditions required the use of a PI controller.
to our simulation approach where we used PI controllers to generate constant volumetric fluxes.Additionally, we saw a similar effect when we changed the lattice resolution, whereby we got larger variances in droplet length with higher resolution.
In summary, we found a good agreement between the simulations and the experimental data.To our belief, the observed differences were mainly an effect of the simulation being reduced to 2D and that the boundary conditions required the use of a PI controller.

Conclusions
In this work, a multi-phase LBM with the RK model was used to simulate a head-on device with a channel width of 1 mm.Compared to other microfluidic devices, a channel width of 1 mm is relatively large, but is required for relevant volumetric fluxes in chemical synthesis reactors.The simulation on this scale was challenging due to the small capillary number of the flow regime.

Conclusions
In this work, a multi-phase LBM with the RK model was used to simulate a head-on device with a channel width of 1 mm.Compared to other microfluidic devices, a channel width of 1 mm is relatively large, but is required for relevant volumetric fluxes in chemical synthesis reactors.The simulation on this scale was challenging due to the small capillary number of the flow regime.Additionally, we encountered several problems using pressure and velocity boundary conditions with the LMB which were already reported in other publications.We overcame these problems by forcing periodicity of the geometry and the use of a volume force.The volume force acted differently on both fluids and equal volumetric fluxes for both fluids were regulated by a PI controller.
Although the simulations were only performed in 2D, we could reproduce the three main flow regimes-squeezing, dripping, and jetting.The transition between these regimes depended on the Reynolds numbers which were comparable between simulation and experiment.For the experiment, we used a self-built micro-structured head-on device with water and heptane as the continuous and dispersed phases, respectively.We used the observed droplet length as a quantitative measure to compare our simulation results to the experiments.There was relatively good agreement between the simulation and the experiment regarding the length of the droplets, while the variance of the length in the simulation was rather large.The variance was caused by the PI controller.To our understanding, significant differences between the simulations and experiments were due to the comparison of 2D LBM results with 3D experimental data.

Figure 2 .
Figure 2. Change in pressure drop across the droplet surface with bubble radius, in a simulation domain of 100 × 100 lattice units (lu).

Figure 2 .
Figure 2. Change in pressure drop across the droplet surface with bubble radius, in a simulation domain of 100 × 100 lattice units (lu).

Figure 3 .
Figure 3. Force balance between the two fluids and the solid phase.In the lattice Boltzmann method (LBM), the fluids are indicated as hydrophobic (index b) and hydrophilic (index r).Depending on the model parameters, the LBM can be used to simulate a contact angle between 0° and 180° (right).

Figure 4 .
Figure 4. Velocity profile in two-phase layered flow.

Figure 3 .
Figure 3. Force balance between the two fluids and the solid phase.In the lattice Boltzmann method (LBM), the fluids are indicated as hydrophobic (index b) and hydrophilic (index r).Depending on the model parameters, the LBM can be used to simulate a contact angle between 0 • and 180 • (right).

Figure 4 .
Figure 4. Velocity profile in two-phase layered flow.Figure 4. Velocity profile in two-phase layered flow.

Figure 4 .
Figure 4. Velocity profile in two-phase layered flow.Figure 4. Velocity profile in two-phase layered flow.

Figure 5 .
Figure 5. (a) Standard geometry of the T-junction; (b) modified geometry for periodic boundary conditions.

Figure 5 .
Figure 5. (a) Standard geometry of the T-junction; (b) modified geometry for periodic boundary conditions.

Figure 6 .
Figure 6.Sensor positioning and the flow direction to adjust the flow rate.Figure 6. Sensor positioning and the flow direction to adjust the flow rate.

Figure 6 .
Figure 6.Sensor positioning and the flow direction to adjust the flow rate.Figure 6. Sensor positioning and the flow direction to adjust the flow rate.

Figure 7 .
Figure 7. Five steps of droplet formation at the T-junction; the simulation is shown on the left (channel width 10 lu) and the experiment is shown on the right.

Figure 7 .
Figure 7. Five steps of droplet formation at the T-junction; the simulation is shown on the left (channel width 10 lu) and the experiment is shown on the right.

Figure 9 .
Figure 9. Droplet generation process in middle flow for 40 lu.

Figure 9 .
Figure 9. Droplet generation process in middle flow for 40 lu.

Figure 10 .
Figure 10.Experimental set-up for the flow experiments (left), and the head-on-device made from poly (methyl methacrylate) (PMMA) by mechanical milling (right).

Figure 10 .
Figure 10.Experimental set-up for the flow experiments (left), and the head-on-device made from poly (methyl methacrylate) (PMMA) by mechanical milling (right).

Figure 11 .
Figure 11.Flow rate relationship with droplet size in simulation and experiment.

Figure 11 .
Figure 11.Flow rate relationship with droplet size in simulation and experiment.

Table 1 .
Boundary conditions in simulation and experiment.