## Abstract

A lattice Boltzmann equation method based on the Cahn–Hilliard diffuse interface theory is developed to investigate the bubble formation process in a microchannel with T-junction mixing geometry. The bubble formation process has different regimes, namely, squeezing, dripping and jetting regimes, which correspond to the primary forces acting on the system. Transition from regime to regime is generally dictated by the capillary number *Ca*, volumetric flow ratio *Q* and viscosity ratio *λ*. A systematic analysis is performed to evaluate these effects. The computations are performed in the range of 10^{−4}<*Ca*<1, 1<*Q*<20 and 10^{−2}<*λ*<1, with the equilibrium contact angle varying from 30^{°} to 150^{°}.

## 1. Introduction

The use of two immiscible fluids in microchannel geometries has become a growing topic of interest in many industrials applications. Mixing geometry and fluid properties at inlet conditions are the most important criteria to determine the type of flow that will be generated in a microchannel. T-junction is one of the most used mixing geometries to form two-fluid flows in microchannels. Many experimental and numerical studies have been performed to establish a relationship among geometry, volume flow rate and the dynamics of the flow for gas–liquid systems.

Garstecki *et al.* [1] showed experimentally that the pinch-off process at a low capillary number *Ca* was governed by the balance between the hydrostatic pressures in both fluids. They presented a scaling law that related the dimension of the channel and the volumetric flow ratio *Q* with the final length of the bubble. Van Steijn *et al.* [2] described the bubble formation process using high-resolution, temporally resolved μ-PIV. They were able to predict the final bubble length without any fitted parameters by taking into account the leakage past the bubble. Later, the same authors [3] confirmed that the collapse and pinch-off processes were triggered by a sudden reversed flow of liquid from the tip of the gas thread to the neck. Velocity of the fluid at the corners depends on the size of the thread, which is a measure of the space occupied by the dispersed fluid in the channel. Velocity increases with the length of the bubble.

Numerically, Qian & Lawal [4] presented a study of gas–liquid Taylor flow, and evaluated the influence of inlet conditions, fluid properties, and wetting properties. The study showed that the flow pattern became bubblier and the void fraction increased as the velocity ratio of two fluids was decreased. In addition, variation of the contact angle considerably changed the shape of the slug, making it either concave or convex. Santos & Kawaji [5] identified two different flows in their experiment: slug and stratified flows. The slug flow was further classified into three categories: snapping, breaking, and jetting flows. They also carried out a two-dimensional modelling to capture the effect of surface tension. In addition, void fraction was correlated linearly with homogeneous void fraction. Lastly, it was shown that the velocity slip depended on how the dispersed fluid covered the cross-sectional area of the channel and the presence of liquid at the edges. Numerical studies using lattice Boltzmann method (LBM) are mainly focused on liquid–liquid systems. Van der Graaf *et al.* [6] performed a three-dimensional study, which investigated the effects of *Ca* and contact angle on the volume of the droplets. It was determined that size decreased as *Ca* and contact angle increased. Gupta *et al.* [7] performed a liquid–liquid three-dimensional study for high *Ca*. As a result, the dependence between shear, surface tension and viscosity ratio on the transition between dripping and jetting regimes was established.

The aim of this paper is to present the factors that affect the formation process in a gas–liquid system using an LBM that is based on the Cahn–Hilliard diffuse interface approach [8]. This study will investigate the effect of *Ca*, *Q*, density ratio, viscosity ratio and equilibrium contact angle on the formation process of a bubble in a T-junction microchannel. The paper is organized as follows: §2 introduces the numerical simulation method used in this paper. Section 3 presents the validation of the model using experimental results. Section 4 presents the factors that affect the bubble formation process. Lastly, §5 provides a summary of the key results obtained in this study.

## 2. Two-phase lattice Boltzmann method

LBM is a mesoscopic approach to solve fluid flows by evaluating the evolution of particle distribution functions owing to advection and collision over a discrete lattice mesh. Two-distribution lattice Boltzmann equations (LBEs) are used: one for the pressure evolution and momentum equations, and the other for the advective Cahn–Hilliard equation. The following are LBEs for each distribution function:
2.1and
2.2In the above equations, *g*_{α} and *h*_{α} are the particle distribution functions, **e**_{α} is the α-direction microscopic particle velocity, *ρ* is the mixture density, *C* is the composition of one species, *p* is the dynamic pressure, **u** is the volume-averaged velocity [9], *c*_{s} is the speed of sound and *λ* is the relaxation time. The equilibrium distribution functions are, respectively
2.3and
2.4where *t*_{α} is the weight. *μ* is the chemical potential, and *M*>0 is the mobility in the Cahn–Hilliard equation. *Γ*_{α} is given as . The macroscopic equations recovered from equations (2.1) and (2.2) are
2.5
2.6
and
2.7
where *η* is the dynamic viscosity. *ρ* and *C* are related linearly by *ρ*=*ρ*_{l}*C*+*ρ*_{g}(1−*C*), where *ρ*_{l} and *ρ*_{g} are the density of the liquid and gas, respectively. The mixing energy density of an isothermal system is related to *C* by *E*_{mix}(*C*,∇*C*)=*E*_{0}(*C*)+*κ*/2|∇*C*|^{2} [10], where *E*_{0}=*βC*^{2}(1−*C*)^{2} is the bulk energy, *β* is a constant and *κ* is the gradient parameter. The classical part of the chemical potential (*μ*_{0}) and the thermodynamic pressure (*p*_{0}) are related to *E*_{0} by *μ*_{0}=∂_{C}*E*_{0} and *p*_{0}=*C*∂_{C}*E*_{0}−*E*_{0}. Equilibrium profile is obtained when the mixing energy is minimized and reads *μ*=*μ*_{0}−*κ*∇^{2}*C*=const. in one dimension. The interface profile in a plane at equilibrium is then given by , where *z* is the coordinate normal to the plane interface and *D* is the numerical interface thickness, which is chosen based on accuracy and stability. This study uses *D*=4, a resolution required to resolve the interface region using the isotropic discretization suggested by Lee & Fischer [11], to eliminate the parasitic currents in the system. Having *D* and *β*, one can compute the gradient parameter *κ*=*βD*^{2}/8. Surface tension *σ* is a measure of free energy per unit area at constant temperature, and defined as . The solution of equation (2.7) requires two boundary conditions. The boundary condition for ∇^{2}*μ* ensures no mass flux owing to the chemical potential gradient in the direction normal to the solid boundary, **n**⋅∇*μ*|_{s}=0, where **n** is the unit normal vector. The boundary condition for ∇^{2}*C* can be established by minimizing the total free energy subject to the specified wall free energy. This study uses a cubic boundary condition given as , where *C*_{s} is the composition at a solid surface that generally differs from the bulk or equilibrium composition and *ϕ*_{c} is a constant related to the equilibrium contact angle *θ*_{eq} by . *θ*_{eq} is defined as the angle between the wall and the gas. For detailed derivation of these boundary conditions, see Lee & Liu [9].

## 3. Validation over experimental results

In the T-junction, gas (labelled with subscript g) penetrates the main channel at a constant velocity *v*_{g}, which is already filled with liquid (labelled with subscript l) from a vertical direction and a plug starts to grow. A distortion of the gas in the downstream direction is created by the pressure gradient and the flow in the main channel. The gas stream starts necking and finally breaks, forming a plug. This plug continues downstream in the main channel and the process starts all over again. Simulations are carried out to determine the effect of two numerical parameters: the equilibrium contact angle *θ*_{eq} and density ratio *ρ*_{r}=*ρ*_{g}/*ρ*_{l} in the formation process. Evaluation of each parameter is performed by comparing the final shape of the bubble. Other non-dimensional numbers are defined as capillary number *Ca*=*η*_{l}*v*_{l}/*σ*, viscosity ratio *λ*=*η*_{g}/*η*_{l} and volumetric flow ratio *Q*=*v*_{g}/*v*_{l}, where *v*_{l} is the liquid velocity at the entrance. Results of simulations are presented using non-dimensional time defined as *T*=*t*/*t*_{n}, where *t*_{n}=*H*/*v*_{g}, with *H* being the height of the channel. The values of the non-dimensional numbers are fixed at *Ca*=0.005, *Q*=0.8, *λ*=1 and *T*=12, unless otherwise stated.

### (a) Effect of contact angle θ_{eq}

The effect of contact angle is evaluated by comparing four different cases of *θ*_{eq}=30^{°}, 60^{°}, 90^{°} and 150^{°} for *Ca*=0.005, *Q*=0.8, *λ*=1 and *ρ*_{r}=0.01. Figure 1 shows the final shapes of the bubble, where the contact angle is measured in the liquid phase. As the contact angle decreases, the wall attracts the liquid and rejects the gas, and thus the shape of the slug becomes convex; however, for large contact angle, it tends to become concave.

### (b) Effect of inertia

Liquid–liquid systems have density ratio of the order of unity in magnitude. However, gas–liquid systems can easily reach low density ratio of the order of 10^{−3}. Most previous numerical studies assumed that both gas and liquid have the same density, claiming that at the microscale the buoyancy has little or no effect on the evolution of the flow. To the best of our knowledge, however, no numerical comparison has been yet presented verifying this assumption. Four different density ratios are considered: 1, 0.1, 0.01 and 0.001 at *Q*=0.8 and *Q*=5.0. *Ca* is fixed at 0.005. T-junction microchannel flow with these parameters corresponds to the squeezing regime, where surface tension is the dominant force and the size of the bubble is proportional to *Q*. Figure 1*a*,*b* shows the final shape of the bubble for *Q*=0.8 and *Q*=5.0, respectively. In the case with *Q*=0.8, the pinch-off process is triggered when the density ratio is the smallest, while in the case with *Q*=5.0, the final shapes at the front and middle of the bubble are almost overlapping each other. From these results, it can be seen that the main difference in these cases is concentrated in the area close to the tee, and the major source of this dissimilarity is the inertia. Using a non-dimensional number for the degree of inertia *I*=*Re*/*Ca*, where *Re*=*ρ*_{l}*v*_{l}*H*/*η*_{l}, the inertia of the liquid for the case with *Q*=0.8 is more than six times higher than that for the case with *Q*=5. Wagner *et al.* [12] showed that the increase in inertia produced larger deformation, which in this case leads to bubble pinch-off. Therefore, comparing the given cases, inertia plays an increasingly important role as the value of *Q* is reduced. As the value of *Q* is reduced, the velocity of both fluids becomes more homogeneous. At this point, liquid has to create a higher force to trigger pinch-off, if *ρ*_{r}∼1. However, as *ρ*_{r} decreases, the force will be lower and pinch-off will take place with more ease. Wagner *et al.* presented the evidence which suggested that inertia was a necessary component for the instability that triggers two-dimensional pinch-off. Furthermore, they questioned whether the simulations performed in three dimensions would also undergo the sudden transition that was observed in two dimensions. This implies that setting the densities of two fluids to be the same when simulating a gas–liquid system may not be an acceptable condition unless the ratio of inertia in the two fluids is very low. A density ratio of 0.001 is used in this study.

### (c) Comparison with experiments

Our LBM is assessed by simulating the benchmark cases selected from the experimental study of van Steijn *et al.* [2]. The working fluids are ethanol and air, and the surface of the channel is considered fully wetting. *H*=800 μ*m* and the length of the channel is *L*≈200 *H*. Comparison between the experimental and simulation results is performed using the silhouette of the bubble at different stages of the formation process. Experimental results are presented using the time for a formation cycle, which corresponds to *T*_{c}=170 ms. Figure 2 presents the formation process and the in-plane streamlines near the bubble at four different times. In general, the numerical results for the evolution of the gas at each time are in good agreement with the experimental results. *T*/*T*_{c}=0.03 and *T*/*T*_{c}=0.3 correspond to the filling stage, where the bubble enters and fills the main channel distorting the fully developed pattern. *T*/*T*_{c}=0.75 illustrates the squeezing stage, where the gas has completely blocked the main channel and is being pushed along. *T*/*T*_{c}=1.0 presents the final stage of the process, the pinch-off, where the neck of the thread is reduced to a minimum and the bubble detaches. The case in the experiment corresponds to the squeezing regime, which indicates that the non-dimensional size of the bubble is related to the volumetric flow ratio using a scaling law developed by Garstecki *et al.* [1] given as *L*_{bu}/*H*=*α*_{1}+*α*_{2}*Q*, where *L*_{bu} is the length of the bubble, the constants *α*_{1,2} are of the order of unity and depend on the geometry of the channel. Using *Q*=1.5 and *α*_{1}=*α*_{2}=1.5 as given by van Steijn *et al.* [2], *L*_{bu}/*H*=3.75. Simulated value is *L*_{bu}/*H*=3.94, yielding a per cent difference around 5 per cent.

## 4. Bubble formation process

The influence of the interfacial and viscous forces can be adjusted by varying the values of *Ca*, *Q* and *λ*. The following sections will explore the effect that each of these non-dimensional numbers has on the formation process of a bubble. All these results are performed with *θ*_{eq}=30^{°}, *ρ*_{r}=0.001, and presented at *T*=18.0.

### (a) Effect of capillary number *Ca*

Four values of *Ca* were considered, 0.005, 0.01, 0.1 and 1.0 at *Q*=1.0 and *λ*=0.02. Figure 3 shows changing *Ca* yields three different regimes: squeezing, dripping and jetting. The squeezing regime corresponds to *Ca*<0.01 [1]. In this regime, the length of the bubble decreases as *Ca* increases. The dripping and jetting regimes are reached for higher *Ca*; however, there is no established critical *Ca* in the literature. In this study, *Ca*∼0.1 is used as a transition value. This value is selected by comparing the features of the final shapes. For *Ca*=0.1, the gas stream is thicker than that for *Ca*>0.5. In addition, the tip of the stream for *Ca*=0.1 has a shape of a bullet, juxtaposed to the tear-like shape for other cases. The cases presented in the dripping or jetting regime do not present any pinch-off. This implies that the shear force is not high enough to overcome surface tension.

### (b) Effect of volumetric flow rate *Q*

Figure 4*a* shows the effect of *Q* on the terminal shape for the squeezing regime at *Ca*=0.005 by comparing *Q*=0.5, *Q*=1.5, *Q*=3.0 and *Q*=4.5 at *λ*=0.02. As described previously and as shown in the plot, this regime is mainly dominated by the value of *Q*. As *Q* decreases, the bubble becomes smaller and the flow becomes bubblier. In addition, the non-dimensional length of the bubble as a function of *Q* is compared with the scaling law introduced in §3 [1]. This comparison is presented in figure 4*b* with very good agreement. Figure 5 illustrates the effect of *Q* for dripping and jetting regimes at *Ca*=0.05 and *Ca*=0.1, respectively, at *λ*=0.02, using *Q*=0.3, *Q*=0.5, *Q*=1.0 and *Q*=3.0. The results indicate that the reduction of *Q* decreases the size of the bubble and increases void fraction. Also, pinch-off tends to take place when *Q*<0.5.

### (c) Effect of viscosity ratio *λ*

In order to evaluate the effect of *λ*, four values *λ*=0.02, 0.05, 0.1 and 1.0 are considered with the three values of *Ca* adopted in the previous section. Garstecki *et al.* [1] showed that in the squeezing regime, the size of the bubble was solely affected by *Q*. Therefore, change in *λ* is expected to have little or no effect in the formation process. This assumption is validated by the results presented in figure 6, where *Q*=0.5. Evaluation of *λ* in the dripping and jetting regime is performed using *Q*=0.3 and *Q*=0.5. These values of *Q* were selected to see how *λ* affects flows with and without pinch-off. Figure 7*a*,*c* presents the effect of *λ* for the dripping and jetting regimes, respectively, with *Q*=0.3. As can be seen in the figure, the changing value of *λ* does not have a significant effect on the frequency and the void fraction of the flow. However, the size of the bubble tends to increase as *λ* decreases. Figure 7*b*,*d* shows that the value of *λ* for *Q*=0.5 does not play an important role in the final shape of the flows.

## 5. Conclusion

Bubble formation process in a T-junction microchannel is investigated by systematically changing *Ca*, yielding three distinct regimes. Transition from squeezing to dripping takes place around *Ca*∼0.01. Transition between dripping and jetting is around *Ca*∼0.1. In the squeezing regime, *L*_{bu}/*H* was successfully compared with the experimental scaling law. In addition, it was shown that in the squeezing regime, the viscosity ratio *λ* has very little effect on the frequency, size or location of the bubble formation. Formation process in the dripping and jetting regimes is influenced more by *Q* rather than *λ*. It has been found that the popular computational condition that the density ratio for a gas–liquid system is unity is not valid unless the inertia ratio between the fluids is considerably low. It is suggested that inertia is a necessary component for the instability that triggers two-dimensional pinch-off. Further analysis is required to establish this instability in three-dimensional cases. Surface wettability also plays a very important role in the formation process. As the surface becomes more non-wetting, gas is repelled from the wall, promoting the pinch-off process.

## Footnotes

One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.

- This journal is © 2011 The Royal Society