A modified forcing approach in the Rothman–Keller method for simulations of flow phenomena at low capillary numbers

The lattice‐Boltzmann method (LBM) is becoming increasingly popular for simulating multi‐phase flows on the microscale because of its advantages in terms of computational efficiency. Many applications of the method are restricted to relatively simple geometries. When a more complex geometry is considered—circular and inclined microchannels—some important physical phenomena may not be accurately captured, especially at low capillary numbers. A Y‐Y micro‐fluidic channel, widely used for a range of applications, is an example of a more complex geometry. This work aims to capture the various flow phenomena, with an emphasis on parallel flow and leakage, using the Rothman–Keller (RK) model of the LBM. To this purpose, we modify the forcing term to implement the surface tension for use at low capillary numbers. We compare the simulation results of the RK model with and without the force modification with experiments, Volume of Fluid and the phase field method and observe that the modified forcing term is an improvement over the current RK model at low capillary numbers, and it also captures parallel flow and leakage more accurately than the other simulation techniques.


INTRODUCTION
Microfluidic multiphase flow is a topic of increasing interest because of its applications and the possibilities it offers in various fields.It finds applications in a wide range of fields, be it in biology, 1 pharmaceuticals, 2 chemistry, 3 and radioisotope extraction. 46][7] For these reasons, a large number of experiments studying the flow pattern of these flows have been conducted in channels of various geometries such as Y-Y channels, 8,9 T-T channels, 8 and cross junction channels. 10he flow patterns commonly observed in literature include slug and parallel flow. 8In the case of slug flow, droplets of the dispersed phase form in the continuous phase, and parallel flow is defined as the regime where the continuous and dispersed phases flow parallel to each other in a microchannel.Slug flow generally occurs when interfacial forces are dominant, and parallel flow when inertial or viscous forces are dominant.In addition to these, many flow patterns intermediate to slug and parallel flow have been observed.These patterns vary for different geometries, wettabilities and fluids.Generally, in these intermediate regimes, parallel flow is observed for a certain length of the channel after which the interface deforms to form droplets. 8,11 In certain applications like radioisotope transfer, parallel flow is very important.When we observe parallel flow with the interface located exactly in the middle of the channel, the radioisotopes can be easily transferred from one fluid to another without necessitating a separation step as required in the conventional batch transfer or slug flow. 124][15] However, stable parallel flow with the interface located exactly in the middle is not easy to achieve, with one fluid generally leaking to the outlet of the other. 11This is generally attributed to multiple factors, such as the balance between pressure due to viscous drop and the Laplace pressure due to surface tension, 16 geometry 17 and wettability. 13o have a greater understanding of the flow phenomena, besides experimental evidence, numerical simulations are necessary.Numerous numerical studies have been conducted on microchannels to supplement the results obtained from the experiments.Some of the methods used for simulation include Volume of Fluid (VOF), 18 level-set (LS), 19 phase-field, 11 and lattice Boltzmann methods (LBM). 202][23] However, there are very few studies which numerically simulate the phenomenon of leakage, and the phase-field simulations of Liu 11 were unsuccessful in completely capturing the phenomenon.Many channel modifications have been proposed to minimize or eliminate leakage for LLE purposes. 24,17Compared to experimental fabrication, it is easier to test different designs numerically for leakage and stable parallel flow.For this, however, it is important to successfully simulate leakage in order to understand how geometry and wettability play a role in its development.
Two key advantages of the LBM which accounts for its increased usage in multiphase simulations 20,25,26 are its high degree of localization and lack of interfacial tracking. 27Its localization makes it easily amenable to parallelization. 28he lattice mesh can be extended to complex geometries and the mesoscale nature of LBM makes it very suitable for micro-fluidics. 29Since the phase-field method was unable to capture leakage successfully, the LBM is thus adopted for this study.
The LBM is based on the modelling of the Boltzmann equation. 28Here, fluids are modelled as probabilistic particle distributions, and the evolution of these particle distributions is related to the macroscopic fluid motions to retrieve the Navier-Stokes equations.Various multi-phase lattice Boltzmann (LB) models have been developed over the years, the most common being the Shan-Chen (SC) model, 30 the colour-gradient (CG) or Rothman-Keller (RK) method 31 and the free energy model. 32In the SC method, an interparticle potential is incorporated into the LBM framework to achieve the separation of different phases.The particle equilibrium distribution functions are modified in the free energy model to achieve thermodynamic consistency for flow separation.Each fluid has its individual particle distribution in the RK model and a perturbation operator (PO) is added for flow separation.Each model has its benefits and limitations, but we use the RK model in this work because of its accuracy, ease of application and flexibility as the surface tension is decoupled from the density.Moreover, the RK model has fewer numerical parameters than the free energy model. 33ultiphase channel flow simulations have been performed using the LBM for T channels 25,34 and cross junction channels. 20,35Most of these papers have considered mainly the case of droplet generation and development in a channel, with few like the paper by Liu et al. 20 showing parallel flow as a proof-of-concept, but not exploring further.The geometries considered in LB simulations are generally rectangular, and the outlet section of the channels is rarely modelled.The outlet sections of the channels are especially important when studying parallel flow, as leakage of one fluid into the outlet of another should be avoided in extraction applications.Some papers, like the one by Yu et al., 36 modified the cross channel inlets to make it inclined, but they simulated only droplet flow and operated at higher capillary numbers.
Considering the use of Y-Y micro-channels in a large number of micro-fluidic studies, especially those involving radioisotope extraction, 4,13,14 it is imperative to study the factors affecting the flow phenomena in a Y-Y channel numerically.Most LBM papers study such flows at higher capillary numbers ranging from 0.001 to 1, 20,[36][37][38] while many flows occur at lower capillary numbers. 9,11,39Thus, this article studies the capability of the LBM, specifically the RK model, to simulate the flow phenomena, especially parallel flow and leakage, across a range of low capillary numbers.The implementation of surface tension will be looked into to understand the workings of the model when dealing with such geometries and capillary numbers, and a solution is proposed to enhance the effectiveness of the RK model at such capillary numbers.VOF and phase-field method will be compared with the RK model to examine the capabilities of the simulation methods in each flow regime.
The article is organized as follows.The RK and VOF model is described in Section 2 along with the methods used for implementing the contact angle, and the in-house RK model code is subsequently validated in Section 3. The fluid properties and the channel geometry are given in Section 4. Section 5 details the simulation results which are compared to the experiments and simulations of Liu 11 in a Y-Y channel, and Section 6 discusses the numerical aspects which could influence the simulation results.Finally, the article is concluded in Section 7.

RK model for two-phase flow
The RK model used in this article is similar to the one used by Ba et al. 40 and Halliday et al. 41 The two immiscible fluids are indicated by the colours red and blue respectively.Each of these fluids has its own distribution function represented by f k i , where i is the lattice velocity direction and k is the colour of the fluid.The total distribution is the sum of these two distributions.The evolution of the distribution function is given by: where x and t are the position and time, e i is the lattice velocity in the ith direction, Δt is the time step and Ω k i is the collision operator.The collision operator consists of three parts: (Ω k i ) (1) is the standard single-phase collision operator, (Ω k i ) (2) is the PO responsible for interfacial tension and (Ω k i ) (3) is the redistribution operator which ensures that the two phases remain separated and the interface is maintained.
The Bhatnagar-Gross-Krook (BGK) operator for single-phase collision is given by: in which  is the relaxation time and f k,eq i is the equilibrium distribution function.The relaxation time can be related to the kinematic viscosity by  = c 2 s ( − Δt 2 ), with c s being the numerical speed of sound.It is related to the lattice speed c by c 2 s = c 2 ∕3.The macroscopic parameters can be obtained from the distribution functions: where  = ∑ k  k is the density of the fluid mixture and u is the fluid velocity.We use the modified equilibrium distribution function proposed by Ba et al. 40 to account for the error terms in the Chapman-Enskog expansion of the RK model. 33k,eq i .
Here,  k i is a parameter related to the density ratio of the fluid and is defined as: k is a free parameter related to the speed of sound in each fluid by 5 (1 −  k ), and to the pressure in each phase through the equation of state, 33 The density ratio is therefore related to  by: Since we are concerned with the two fluids of different viscosity, the relaxation times at the diffuse interface should be defined accordingly.We adopt the following scheme defined by Huang et al. 33

𝜏(
where . ≤ 1 is a free and positive parameter that affects the interface thickness, and is normally taken to be 0.98. N is the phase field/colour function indicating how much each fluid is located at a given node.It is defined as: The PO is similar to the one defined by Guo et al. 42 and reads as follows: where A k is the fraction of the interfacial tension contributed by the fluid k, and F s is the body force term which is used to implement the interfacial tension.Using the continuum surface force (CSF) model of Brackbill et al., 43 this term can be expressed as: Here,  is the surface tension coefficient which can be directly implemented in the model, and K is the local curvature of the interface K = −∇ S ⋅ n, with ∇ S = (I − nn) ⋅ ∇ being the surface gradient operator and n the outward-pointing unit normal vector of the interface which is given by: The curvature of the interface in 2D is given by: The local fluid velocity needs to be modified to accommodate the spatially varying body force. 42

𝜌u = ∑
Guo et al. 42 mentions that both the equilibrium (u * ) and local fluid velocities (u) need to be modified (Equation 13) for the model to correctly capture the Navier-Stokes equations and accommodate for spatial and temporal variations of the momentum with the CSF term (Equation 10).However, this modification, along with the CSF term, introduces spurious currents. 41These spurious velocities aren't a problem at higher capillary numbers but could significantly influence the result for cases when surface tension forces are dominant.The CSF term has been known to introduce spurious velocities not only for LBM, but also for Volume-of-Fluid and level-set methods. 44,45o accommodate lower capillary flows, we propose to neglect the velocity modification for u * and apply it only for u, that is, m = 0 for u * (Equation 13) This approach was used by Ladd and Verberg 46 for their forcing term.Their forcing term was shown to introduce additional terms in the momentum equation by Guo et al. 42 In summary, we apply Ladd and Verberg's velocity modification to the forcing term by Guo et al. 42 described in Equation ( 9) (Current approach, m = 0) and compare it with the standard approach proposed by Guo (Guo approach, m = 0.5).
Finally, we have the redistribution step, (Ω k i ) (3) , which is applied to enhance phase separation of the two fluids. 47 where  is a numerical parameter which influences numerical stability and interface thickness.A value of 0.7 is used for , which is the most common value mentioned in literature as it provides the best balance between accuracy, stability and spurious velocities. 40,48 i is the angle between the phase field gradient and the lattice direction, and it is given by: The partial derivatives for the phase-field gradient and curvature are calculated using a 9-point finite different stencil as evaluated for any variable  by Liu et al. 49

Boundary conditions
The halfway bounce-back method is applied at all channel walls to impose a no-slip boundary condition.At both inlets, the velocity is prescribed using the moving walls boundary condition method of Ladd. 50The modified pressure boundary condition method described by Douillet-Grellier et al. 51 is used at both outlets.

Contact angle at an inclined plane
This subsection details the contact angle implementation (CAI) at the boundary.A recent approach for applying contact angle in non-flat or inclined boundaries was developed by Leclaire et al. 52 and later modified by Xu et al. 48for greater accuracy and stability.According to this approach, the lattice nodes are divided into four categories: • C FB : the nodes which are located in the fluid domain and in contact with at least one node in the solid domain; • C FL : the nodes which are located in the fluid domain and not in contact with any node in the solid domain; • C SB : the nodes which are located in the solid domain and in contact with at least one node in the fluid domain; • C SL : the nodes which are located in the solid domain and not in contact with any node in the fluid domain.
The colour gradient ∇ N for the points at C FB requires the phase field at C SB according to the 9-point stencil in Equation (16).The contact angle is implemented at C FB through the modification of the colour gradient.The colour gradient is modified in such a way that the interface normal is at an angle  as shown in Figure 1.The value of the phase field at C SB is first obtained through a weighted average of its nearest fluid node neighbours.
The colour gradient at C FB can now be determined using Equation (16).Since the phase field determined at C SB is an estimate, the colour gradient at C FB is also an estimate.Thus, the normal at the interface is also an estimate.∇ N * (Figure 1).The predicted unit normal at the interface is given by: This predicted normal is compared to the two interface normals, n 1 and n 2 .The interface normals are determined from the normal perpendicular to the surface n s and contact angle.The surface normal is given by: where c l is the lth mesoscopic velocity associated with the isotropic discretization, s(x) is an indicator function that equals 1 for x ∈ C S and 0 elsewhere, and w is the weight function associated with the discretization.An eight-order discretization is used as it minimizes spurious velocities.The weights associated with these discretizations are the same as those specified in Sbragaglia et al. 53 Two interface normal vectors are possible for a given contact angle, and these can be determined from the surface normal.
The appropriate unit normal vector is chosen by calculating the Euclidean distances D 1 and D 2 from the predicted normal, and then selecting the correct interface normal.
The modified colour gradient is obtained by

Volume of Fluid
Ansys FLUENT 2022 R2 is used to simulate VOF.The Navier-Stokes equations for the individual fluids are simulated in VOF. 18The position of the interface is tracked using the volume fraction continuity equation and the density and viscosity of the two-phase mixture are described by: where  and  are the density and viscosity of the mixture respectively,  is the volume fraction of liquid 1.To include the surface tension, a source term is added to the Navier-Stokes equation based on the CSF model described by Brackbill et al. 43 This is given by: where  is the surface tension and  is the curvature.Equation ( 23) has a similar form to Equation (10), and this is expected as both are based on the CSF model.As mentioned earlier, the CSF term is known to induce spurious currents in VOF. 44,54The main reason for this is because the Laplace relation of surface tension and pressure difference is not exactly satisfied. 44For higher capillary numbers, this inaccuracy doesn't have a big impact on the solution, but for low capillary numbers, large spurious velocities are developed.Besides, the same gradient operator is used to calculate the smooth pressure term and the discontinuous surface tension term, which leads to an inaccurate estimation of gradients and curvature. 44One solution for this is to combine VOF with LS, as LS involves a continuous function for implementing surface tension. 55Instead of a volume fraction equation, an equation for the level-set function is solved.This option is available in FLUENT and the model is the same as that in the FLUENT theory guide. 56We run simulations with both methods for the Y-Y channel to understand the impact of the forcing term on simulations of low capillary number flows and leakage.

CODE VALIDATION
The LBM simulations are run using the methodologies described in Section 2. The model is first validated by running static contact angle simulations as described in the paper of Xu et al. 48Both the Guo and Current approaches are applied (Equation 13), and the same results are obtained.To validate a more dynamic case, the T-channel experiments of Graaf et al. 25 are simulated in 2D and the droplet diameters are compared for both approaches.The T-channel has two inlet junctions perpendicular to each other and an outlet located along the inlet of the continuous phase (Figure 2).The width of the channel is 100 μm.The densities of the continuous and dispersed phases are taken as  c = 1000 kg∕m 3 and  d = 1000 kg∕m 3 .The densities of the fluids used by Graaf et al. 25 are very similar to each other, so they F I G U R E 2 T channel setup for droplet simulations. 25[Colour figure can be viewed at wileyonlinelibrary.com] are made equal in the simulation for convenience.The dynamic viscosities are the same as those of the experiments, with  c = 1.95 × 10 −3 Pa s and  d = 6.71 × 10 −3 Pa s.A contact angle of 45 • with respect to the continuous phase is implemented and the surface tension is  = 5.0 × 10 −3 N/m.The capillary number of the dispersed phase is kept constant at Ca d = 0.011 while the continuous phase Ca c is varied from 0.0033 to 0.066.
As shown in Figure 3, the droplet forms the expected shape and follows the expected steps of expansion and necking during formation. 25The droplet diameters of the experiments and simulations are plotted in Figure 4. We observe an interface thickness of around 3-4 lu, and this is the expected thickness for  = 0.7. 57Since the experiments are in 3D, Graaf et al. 25 calculated an equivalent spherical diameter from the volume of the droplets.Following their procedure, we, therefore, calculate the equivalent circular diameter from the droplet area in our 2D simulations.In order to make a valid comparison, the Strouhal numbers (St) are equalized in 2D and 3D as this number is a measure independent of the dimensions and reflects the dynamics of droplet breakup.Hence, where f is the droplet frequency, L is the characteristic length of the flow (hydraulic diameter) and U is the flow velocity.The characteristic lengths in 2D and 3D are the same as the width and height are equal.The fluids flow only along the x direction in the main channel, so the velocities are also equal in both dimensions.This implies where Q is the flow rate, A is the area of the slug and V is the volume of the slug.Let d 2D be the equivalent circular diameter of the 2D slug and d 3D the equivalent spherical diameter of the 3D slug.Therefore, a 3D diameter is calculated from the 2D results using the results from Equation ( 25 From Figure 4, we can see that the droplet diameters closely match the experiments for the Guo approach for a wide range of Ca numbers, whereas the current approach is accurate for low Ca numbers (Ca < 0.016).Since the focus of this work is at very low values (Ca < 10 −3 ), the current approach is considered sufficient for this study.

NUMERICAL SETUP
The simulations are conducted in 2D and compared to experimental studies on the Y-Y channel by Liu. 11Liu used a glass microchip manufactured by the Institute of Micro-chemical Technology (IMT) as shown in Figure 5.The channel dimensions and the corresponding lattice dimensions are shown in Table 1.Water and toluene are the two fluids used for comparing parallel flow with Liu's results.The properties of these fluids are described in Table 2, along with the inlet velocities (u in ) at Ca = 10 −4 , as this would provide the reader with an idea of the velocity magnitude when talking about spurious velocities.The advancing and receding angles are unknown in the experiments and the implementation of the contact angle as described in Section 2.3 also doesn't allow for hysteresis, so hysteresis will not be included in this study.This table also shows the values of  k (Equation 6) and , the parameter related to interface thickness in Equation ( 14).The

TA B L E 2
Properties of the fluids used in the simulations.chosen relaxation time is close to the stability limit of 0.5. 28Such a low relaxation time is chosen so that the computational load is smaller, as a larger relaxation time would necessitate very large domains for the same capillary numbers.The flow is initialized in the simulation domain as shown in Figure 6.The red fluid corresponds to the aqueous phase and the blue fluid is the organic phase.Fully developed laminar flow is assumed at the inlet, so a shorter inlet length of 500 microns can be used.A channel length of 2 mm is used instead of the actual length of 2 cm.This assumption was justified by running the simulations for both channel lengths, 2 mm and 2 cm, with no difference in the results.

RESULTS
We first simulate the flow map obtained by Zheng Liu 11 for a water-toluene mixture.The capillary numbers applied are in the range of 3 ×10 −5 to 2 ×10 −3 .These values are relatively low compared to those used in literature, the lowest capillary number used in the literature for the RK model being 10 −3 by Wu et al. 35 and Liu et al., 49 both of whom simulated fluid flow in a cross junction channel.Both these papers compared their results with experiments, and in these experiments, the lowest capillary number was 10 −3 .The results, consisting of the flow pattern and leakage, are compared for both VOF and RK models with the experiments and phase-field simulations of Liu, 11 who observed three regimes-slug, transition and parallel flow.
The flow pattern map obtained from the experiments and simulations for water-toluene mixtures is shown in Figure 7.The Current approach and VOF are able to capture parallel and slug flow at the respective capillary numbers, but not the relatively narrow region of transition flow.The phase-field method, however, is able to successfully capture all the flow regimes-slug, parallel and transition flow.However, the benefits of the Current approach compared to the phase-field method is made clear when simulating leakage, and this will be discussed in Section 5.2.
In the case of the Guo approach, it was difficult to observe grid-independent results, especially when simulating parallel flow, but slug flow was successfully captured.Grid convergence could not be obtained for parallel flow at all the expected capillary numbers and this will be further discussed in Section 6.2.For this reason, the regimes obtained from the Guo approach are not plotted in Figure 7.

F I G U R E 6
Flow initialization in all simulations.The red fluid is the aqueous phase and the blue fluid is the organic phase.

Slug length
In this subsection, we look at the capabilities of each technique to accurately capture the slug length at various low capillary numbers.We once again use the experimental results of Liu, 11 but this time for water-heptane instead of water-toluene as no data on the slug length is available for the water-toluene system.The flow map of water-heptane is not presented in Figure 7 for clarity, as the same findings were also observed here.The slug lengths of all the techniques (VOF, phase-field, Guo and Current approaches) are compared with the experiments in Figure 8.
The slug lengths remain more or less constant in the case of the Guo approach and VOF with varying capillary numbers.The current approach, on the other hand, shows the right trend of slug lengths decreasing with increasing capillary number, but the lengths differ from the experiments by around 30%.The phase-field method is the only technique which is able to both capture the right trend and show accurate slug lengths, with the maximum error observed being around 10%.

Leakage
The outlet section of the Y-Y channel is now modelled in the case of parallel flow because we want to observe the phenomenon of leakage, where one fluid leaks into the outlet of the other. 11Perfect flow separation and stability are important when we consider the transfer of species from one fluid to another, 14 so it is imperative that simulation techniques (Current approach, VOF and phase field) are able to capture this phenomenon accurately.The Guo approach is not considered here as we couldn't obtain grid-independent results for parallel flow.The simulations are compared with the experiments of Liu. 11 Three different leakage regimes are observed when using the Current approach-leakage to the toluene outlet, leakage to the water outlet and leakage to both outlets.These regimes are shown in Figure 9.Only two leakage regimes are observed in the experiments, however, leakage to the toluene outlet and leakage to the water outlet, and the experimental regimes are shown in Figure 10.
Liu did not discuss the nature of the leakage regime near the outlet, that is, whether the fluid leaks in the form of droplets or parallel flow, so we restrict our discussions regarding the nature of the leakage regime to the simulations.In the case of leakage to the toluene outlet, the flow of water in the LBM simulations is always parallel.In the case of leakage to the water outlet, leakage occurs either in the form of droplets of toluene or as parallel flow.Droplets are observed in the case of leakage to both outlets.
The leakage regimes are plotted in Figure 11, where the capillary number of toluene and the ratio of flow rates have been varied.The figure shows that, while increasing the flow rate ratio, the leakage regime goes from leakage to the water outlet in the shape of droplets, to leakage to the water outlet as parallel flow, to leakage to both outlets and finally to leakage to the toluene outlet as parallel flow.
When comparing the results of the current approach with the experimental results of Liu 11 in Figure 11, the same trend is roughly visible, that is, leakage to the water outlet at low flow rate ratios and leakage to the toluene outlet at high ratios.The numerical simulations, however, show a transition regime with leakage to both outlets at flow rate ratios between 0.8 and 1, whereas the experiments don't.
In the case of VOF and VOF-LS, we observe droplets being generated near the intersection for all capillary numbers (Figure 12).Though leakage is observed, these results do not match the experiments of Liu.The phase-field simulations of Liu were also unable to capture leakage correctly, with the results showing either no leakage at all or leakage similar to that observed in VOF.
The Current approach, thus, produces the most accurate results when it comes to leakage compared to both VOF and phase-field.A possible explanation for the results observed using the Current approach can be gleaned from the findings of Latva-Kokko and Rohtman. 58They observed that the implementation of the RK model generates an artificial slip velocity, leading to a dynamic contact angle even though it wasn't specified.The role played by dynamic contact angles in leakage is less known, though it must be mentioned that Liu observed leakage, albeit inaccurately, in his simulations when he included an expression for a non-equilibrium dynamic contact angle, and no leakage when this expression was not included. 11Even when we consider the Guo approach at capillary numbers larger than 2 ×10 −3 , leakage was accurately captured for most flow rate ratios, except for 0.8-1 where no leakage was obtained.
Obviously, more study is needed to see how dynamic contact angle and slip are related to leakage.However, this section clearly shows the ability of the Current approach to capture leakage.

DISCUSSION
To understand the reasons for the results shown in Section 5, we dive a little deeper into the numerical aspects of each technique in this section, with a particular emphasis on spurious velocities.The current and Guo approaches, VOF and VOF-LS are looked into in the following subsections.The deficiencies of the phase field method are discussed in detail by Liu, 11 so we won't be discussing them here.

Current versus Guo
Figure 13 shows the capability of the Guo and Current approaches in capturing parallel flow in the microchannel.At a toluene capillary number of 10 −3 , the flow is parallel when using the Current approach (Figure 13C) which is the expected experimental result, whereas slug flow is obtained when applying the Guo approach (Figure 13A).A possible explanation for this difference can be found in the magnitude of the spurious velocities, which are created by the CSF term.Figure 14 shows the contours of the u x ∕u max (u max being the maximum expected fluid velocity at the centre of the main channel) for both cases at Ca = 10 −3 .While the Current approach shows the correct velocity profile and maximum velocity (Figure 14B), Figure 14A points to larger velocities at the interface when using the Guo approach (u x ∕u max =1.4), especially near the intersection of both inlets.These large spurious velocities are in line with the results of Xu et al, 48 where the largest velocities were observed at the contact points of the fluid-fluid and fluid-solid interfaces, especially at corners.The spurious velocities can be more prominently visualized at a lower capillary number of 10 −4 , as shown in Figures 15  and 16.At this range, we can see from Figure 15B that the velocities at the interface using the Current approach are only a little larger than the maximum fluid velocity (the ratio is around 1.2).On the other hand, the spurious velocities observed in the Guo approach are nearly an order of magnitude larger than the expected fluid u max (Figure 16B).At a capillary number of 10 −3 , the maximum fluid velocity is lesser than the interfacial velocities in the Guo approach, so these spurious velocities compete with the expected flow phenomena.If we operate at a higher capillary number, where the fluid velocity is larger than the spurious velocities, we should obtain parallel flow.At a capillary number of 2×10 −3 , the maximum fluid velocity is larger than the interfacial velocity of around 0.04, and, as seen in Figure 13B, parallel flow is indeed obtained.
To conclude, in LBM, the Current approach is more practical in predicting parallel flow in microchannel than the Guo approach at low capillary numbers (Ca < 10 −3 ).

Volume of Fluid
Since the VOF method uses the CSF term for implementing surface tension, one would expect similar problems as observed for the Guo approach.If we look at the u x contours at Ca = 10 −4 , we can see that the spurious velocities are very high for VOF (Figure 17B), and lower for VOF-LS (Figure 17C).The velocities are observed to be the largest at the intersection, with the spurious velocities being an order of magnitude larger for VOF and nearly 6 times larger for VOF-LS.Note that though the VOF-LS approach does reduce the spurious velocities present in VOF, it still requires more modifications to operate at very low capillary numbers.In summary, looking at Figures 14A and 17B,C, we can clearly see that the CSF term introduces large interfacial velocities regardless of whether we use VOF or LBM.This could be a possible reason why we observe constant slug lengths for VOF and the Guo approach in Figure 8, and why we observe the expected trend once we reduce the spurious velocities using the Current approach.

Grid convergence
One way of checking for grid independence is by plotting the slug lengths across different grid resolutions and observing the point of convergence, 11,21,25 that is, when the slug length remains constant with further refinement.It has been found that both the VOF and RK models show converged slug lengths at a mesh size of 8 μm.
A more strict criterion for grid convergence is proposed here, which is based on the ability of the scheme to capture parallel flow.Figure 18 shows the results for VOF as well as the RK approaches.This figure covers Ca numbers where parallel flow is to be expected (see Figure 7).
In Figure 18, we can see that the Current approach converges the quickest.Even for large mesh sizes, parallel flow is predicted for most Ca numbers.In the case of VOF and the Guo approach, however, the presence of large spurious velocities considerably influences the flow regime as indicated by Figure 13A, favouring slug instead of parallel flow.The Guo approach doesn't converge even at a resolution as low as 0.8 μm, hence finer domains are required for convergence.Both the VOF and VOF-LS models show slug flow in coarse grids, and capture parallel flow at a mesh size of 5 μm (Figure 12) after grid refinement, despite the large spurious velocities.
In conclusion, the Current approach predicts the transition from parallel flow to slug flow at relatively coarse meshes, which is favourable in terms of calculation time.

Comparison of the simulation techniques
The capabilities of the different approaches in capturing the flow phenomena at relatively low Ca numbers in a Y-Y channel are summarized in Table 3.We can clearly see that no technique is able to successfully and accurately capture all the flow phenomena, but overall, the phase-field and Current approaches are the most successful.The current approach is a better option when simulating parallel flow and leakage, while the phase-field approach is a good choice for capturing slug lengths and the overall flow map.It must be mentioned that considerable care has to be exercised in choosing the numerical parameters in the phase-field method, especially the numerical Peclet and Cahn numbers.Different Peclet numbers lead to changes in the TA B L E 3 Note: Based on the simulation results, each technique is accorded a value for each category.X corresponds to not capturing the phenomenon at all, 0 to qualitatively capturing the overall phenomenon but unable to capture the finer details (e.g., VOF showing constant slug lengths with varying capillary numbers and not the expected trend), + to successfully capturing the finer details, although inaccurately (e.g., inaccurate leakage regimes) and ++ to accurately capturing all the details.
observed simulation results when it comes to slug lengths and leakage, 11 and this limits the flexibility of the phase-field method.The current approach, on the other hand, doesn't have as stringent a requirement on the value of its numerical parameters, thus providing an added advantage when simulating multiphase flows at low capillary numbers.

CONCLUSION
This article examines the RK-LBM model's capability and accuracy in simulating flow regimes in a complex geometry-a Y-Y channel-under convective conditions and at low capillary numbers.The approach proposed by Xu et al. 48was used for implementing the contact angle.The forcing term proposed by Guo 42 was modified to reduce the magnitude of the spurious velocities and extend the application of the RK model to lower capillary numbers.The current approach is clearly an improvement over the existing Guo approach at low capillary numbers, showing quick grid convergence and successfully capturing leakage.Additionally, it also shows the expected trend of slug lengths decreasing with capillary number as opposed to VOF and the Guo approach.More accurate slug lengths are observed in the phase-field method, but the phase-field method isn't able to capture leakage accurately.Thus, the current approach is a viable option at lower capillary numbers especially when it comes to simulating parallel flow and leakage.
More studies are needed to come up with better solutions to capture intermediate phenomena such as transition and deformed interface flows.Modifications to the CSF term have been done in Volume-of-Fluid methods, 54,59 but these modifications have not been applied to complex cases such as micro-channel flows.Future models would have to work along the lines of accurate determination of the curvature (Equation 12) for improved results.

F I G U R E 1
Image of the wetting boundary condition.n s is the unit vector normal to the wall, n 1 and n 2 are the two interface normals, n * is the predicted normal and  is the contact angle.[Colour figure can be viewed at wileyonlinelibrary.com] ): F I G U R E 3 Droplet formation and development in a T channel at Ca c = 0.016.The red fluid is the continuous phase and the blue fluid is the dispersed phase, with the diffuse interface represented by the transitioning of colours from blue to red.The droplet develops in three stages: (A) Expansion, (B) necking (C) detachment. 25[Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 4 Comparison of droplet diameters obtained from the simulations using Guo (blue diamonds) and Current (black squares) approaches with the experiments (red circles) of Graaf et al. 25 for varying Ca C at a constant Ca D .[Colour figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 5 Y-Y microchannel used byLiu,11 with two inlets and outlets.TA B L E 1Channel and domain dimensions in physical and lattice units.
figure can be viewed at wileyonlinelibrary.com]

F I G U R E 9
Leakage regimes observed in LBM simulations of water (red)-toluene (blue) mixture.(A) Leakage to toluene outlet, (B) leakage to water outlet, (C) leakage to both outlets.The capillary number of toluene is fixed at 2 × 10 −3 and the flow rate of water is varied to obtain different regimes flow rate ratio is 1, 0.6 and 0.8 for (A)-(C) respectively.[Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 10 Leakage regimes as observed in the experiments 11 for water (dark)-toluene (light) mixture.(A) Leakage to toluene outlet (B) leakage to water outlet.[Colour figure can be viewed at wileyonlinelibrary.com]

F I G U R E 11
Comparison of leakage regimes of LBM simulations (current approach, black symbols) with experiments (regions above and below the horizontal line).The vertical axis corresponds to the ratio of the aqueous to the organic flow rate.A flow rate ratio of 0.66 (black horizontal line) points to the region where leakage switches from the water outlet to the toluene outlet. 11[Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 12 Parallel flow for water-toluene mixture simulated using VOF and VOF-LS.The red fluid is water and the blue fluid is toluene.[Colour figure can be viewed at wileyonlinelibrary.com]

F I G U R E 13
Flow pattern of water (red)-toluene (blue) system observed for the (A) Guo approach, Ca = 10 −3 , (B) Guo approach, Ca = 2 × 10 −3 and (C) Current approach, Ca = 10 −3 .[Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 14 Normalized u x contour for the (A) Guo and (B) Current approaches at a capillary number of 10 −3 .[Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 15 Slug flow for a capillary number of 10 −4 using the Current approach, with the same velocities for both fluids.(A) Density contour, (B) normalized u x contour.The red fluid is water and the blue fluid is toluene, with the diffuse interface represented by the transitioning of colours from blue to red.[Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 16 Slug flow for a capillary number of 10 −4 using the Guo approach, with the same velocities for both fluids.(A) Density contour, (B) normalized u x contour.The red fluid is water and the blue fluid is toluene.[Colour figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 17 Slug flow for a capillary number of 10 −4 using VOF and VOF-LS, with the same velocities for both fluids.(A) Density contour of VOF and VOF-LS, (B) normalized u x contour using VOF (C) normalized u x contour using VOF-LS.The red fluid is water and the blue fluid is toluene, with the diffuse interface represented by the transitioning of colours from blue to red.[Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 18 Grid convergence for the Current (circles) and Guo approaches (asterisks), and VOF (crosses).The vertical axis corresponds to the capillary number of toluene at which parallel flow is first observed in the simulation and the horizontal axis is the mesh size.Slug flow is observed for all the capillary numbers below the respective lines and parallel flow is above the lines.[Colour figure can be viewed at wileyonlinelibrary.com] Table comparing the capabilities of simulation techniques in capturing flow phenomena at low capillary numbers.