Study of Two Layered Immiscible Fluids Flow in a Channel with Obstacle by Using Lattice Boltzmann RK Color Gradient Model

. Lattice Boltzmann method (LBM) is employed in the current work to simulate two-phase flows of immiscible fluids over a square obstacle in a 2D computational domain using the Rothman-Keller color gradient model. This model is based on the multiphase Rothman-Keller description, it is used to separate two fluids in flow and to assess its efficacy when treating two fluids in flow over a square obstacle with the objective of reducing turbulence by adjusting the viscosities of the two fluids. This turbulence can cause major problems such as interface tracking techniques in gas-liquid flow and upward or downward co-current flows in pipes. So, the purpose of the study is to replace a single fluid with two fluids of different viscosities by varying these viscosities in order to reduce or completely eliminate the turbulence. The results show that to have stable, parallel and non-overlapping flows behind the obstacle, it is necessary that the difference between the viscosities of the fluids be significant. Also, showing that the increase in the viscosity ratio decreases the time corresponding to the disappearance of the vortices behind the obstacle. The results presented in this work have some general conclusions: For M≥2, the increase in the viscosity differen ce leads to an increasing of friction between fluids, reducing of average velocity of flow and decreasing the time corresponding to the disappearance of the vortices behind the obstacle. However, for M≤1/2, the opposite occurs.


Introduction
Immiscible fluids flow in a 2D channel are widely encountered in many scientific and industrial applications and they remain one of the challenges in computational fluid mechanics.These involve filtration systems, heat exchangers, recovery of hydrocarbons from oil and gas reservoirs, capillary networks, battery slurries, highly conductive pastes for solar cells and many other processing systems and many other that are a fundamental problem in fluid mechanics.Immiscible fluids e.g. oil and water can be described using chemical concept, when the fluids mix together, their bonds are not broken in order to form new bonds, and the two fluids will not form one cohesive solution.Numerically, several studies have been carried out (Bitsch et al., 2014, Bitsch et al., 2016, Schneider et al., 2016, Schneider et al., 2017), among these methods, we find the lattice Boltzmann method (LBM) which has achieved considerable success to simulate hydro-thermodynamics problems.It is especially suitable for the direct numerical simulation of immiscible multiphase fluid flows and also a numerical approach of the computational fluid dynamics (CFD) in engineering.It focuses on constructing simplified kinetic models that incorporate the physics of microscopic processes from which the macroscopic flow characteristics are computed on the basis of the Navier-Stokes (NS) macroscopic equations.LBM has been used to study several fluid flows problems such as viscous two layered fluids flow through a channel (Leclaire et al., 2012), wetting and spreading phenomena (Huang et al., 2007), collision and bubble rising phenomena (Huang et al., 2014), Rayleigh-Taylor instability (He et al., 1999).
Multiphase LBM models have been developed in past years, including the color-gradient model originally proposed by Rothman and Keller in 1988 and is labeled as Rothman-Keller (RK) model (Rothman et al., 1988).It is developed by Gunstensen et al. (Gunstensen et al., 1991) and Grunau et al. (Grunau et al., 1993) for binary fluids with different densities and viscosity ratios.Then, the pseudo-potential proposed by Shan and Chen (Shan et al., 1993) which is the very popular model in the LBM community for reason of conceptual simplicity and computational efficiency.Well after, Inamuro et al. (Inamuro et al., 2004) improved the free-energy model of Swift et al. (Swift et al., 1995) to achieve higher density ratios and make it more satisfactory.The fourth model is the interface tracking proposed by He et al. (He et al., 2000).Flow control is one of the main objectives of the fluid mechanics community and it is important in a wide range of systems and devices.For this, we are interested in this phenomenon in the present study.This paper focuses on 2D RK color-gradient model which is developed for immiscible two-phase fluids, the basic idea behind this model lies on the introduction of two distribution functions labeled 'red' and 'blue' to represent two different fluids, through adding an extra binary fluid collision operator (perturbation term) to the lattice Boltzmann equation to generate the interfacial tension.In the RK model the surface tension, the density and viscosity ratios can be adjusted independently.For this reason, we have chosen to use this model.However, RK model gives poor results for two-phase parallel flows with different densities in a channel.As a result, to overcome this deficiency, Huang et al. (Huang et al., 2015) found some extra terms in the recovered momentum equation.These extra terms may affect the numerical results significantly.Recently, several works are conducted using the RK color-gradiant as Lafarge et al. (Lafarge et al., 2021) who improved color-gradient method for lattice Boltzmann modeling of two-phase flows, Sadeghi et al. (Sadeghi et al. 2021) who studied immiscible displacement mechanisms in pore doublets and Mora et al. (Mora et al., 2021a(Mora et al., , 2021b) ) who investigated the optimization isotropy of the gradient at small radius of curvature interfaces such as those that occur in flow through a porous medium for the choice of the interfacial thickness parameter =0.5 and the optimal surface-tension isotropy in the RK color-gradient lattice Boltzmann method for multiphase flow.
In this paper, we study two fluids having an immiscible character to determine their behaviors through a square obstacle by varying their viscosity ratio.This study is organized as follows: in Section 2, we describe the color-gradient LBM model.In section 3, we present the validation of our numerical code.In section 4, we discuss the numerical results obtained.The conclusion is presented in the last section.

𝑘
. Each colored distribution function undergoes the collision and streaming steps.Thus, the evolution of the distribution of particles    (, ) is expressed by the following equation: where   + (, ) represents the distribution function after the recoloring step and  is the lattice spacing time which is equal to 1 in the simulation.The collision step can be described as (Mora et al., 2021): where    * (, ) is the post-collision state, both terms (Ω   ) 1 and (Ω   ) 2 are used to describe two collision steps.The Bhatnagar, Gross and Krook (BGK) approximation is adopted for the first term (Ω   ) 1 as (Behrend et al., 1994): where   = 1/√3 is the speed of sound and   represents the weight factors defined by: We can obtain the density of the fluid   (, ) and the velocity components (, ) by the following equations: Theoretically, it is true for   ≥   and this constraint 0 ≤   ≤   < 1 must be respected in order to avoid negative pressures which are presented as: The relaxation parameter  can be given by the interpolation scheme constructed by Grunau et al. (Grunau et al., 1993) is a free positive parameter that affects interface thickness and is usually set as 0 <  ≤ 1.  is a function that takes the value 1 or −1, depending on whether it is evaluated at a position that contains only the red fluid or only the blue fluid.
To get the average viscosity of the two fluids components, we approximate the viscosity as an extensive variable and it can be expressed by: In the RK model, the surface tension is modeled by means of the perturbation operator which is defined in the literature by (Mora et al., 2021): where   is a parameter that affects the interfacial tension, (, ) is the color-gradient which is calculated by: The interfacial tension in the N-S equations can be recovered by adding the correct term   presented in Eq. ( 17) which is defined for the  2  9 model as follows: The recoloring step for each component is used to achieve separation of the two fluids and to maximize the amount of red fluid at the interface sent into the red fluid region and the amount of blue fluid at the interface sent into the blue fluid region, while respecting the principles of conservation of mass and total momentum.It was explained in (Latva-Kokko et al., 2005) and it presented by the following equations: ,  is the parameter which adjusts the interfacial thickness, its values are between 0 and 1 to ensure positive particle distribution functions and cos(  ) =   . |  |. || ⁄ is the cosinus of the angle between the color gradient (, ) and the discrete velocity   .   (,  = 0) represents the equilibrium distribution function, which are evaluated using zero speed.We can note that when two components have identical densities, it is not necessary to calculate both collision step Eqs. ( 4) and ( 17) separately for each component.The two collision steps become: where  = ∑   /2  and   = ∑ . A is a parameter, which determines the interfacial tension.

Elimination of the unwanted extra term
The RK model gives poor results for two-phase parallel immiscible fluid flows with different densities in a channel.Therefore, Huang et al. (Huang et al., 2013) found some extra terms which added to the lattice Boltzmann transport equation.These terms are defined by:  24).This extra term can be considered as a forced term in N-S equations and can affect numerical results significantly.It can be defined as (Huang et al., 2015): Where The derivative terms appear in the extra unwanted term are evaluated over y-direction and we adopt the velocity   .Because, the vertical velocity   is assumed to be zero everywhere inside the computational domain and the flow is steady over x-direction.Therefore, the flow also depends on the x-direction and the derivate of parameter over x-direction is equal to zero    = 0, where  denotes density, velocity, and pressure.The central finite difference method is used to calculate the density gradient and derivatives terms   (  ) and   (    (  )) as (Huang et al., 2015): The N-S Eq. ( 24) can be simplified as (Huang et al., 2015)       2 (  ) + G is a body force.In our simulation, the body force is added to lattice Boltzmann equation as (Nie et al., 2015):

Drag coefficient expression
The average drag coefficient   is a critical parameter that is calculated in the channel to quantify the drag resistance of the obstacle, here we define it by: where     and     are the average drag coefficients of fluid 1 and fluid 2, respectively.They are defined by: (, ) is the drag force which is related to the velocities   and distribution function    (, ) by: where  −  is the reverse direction of the distribution function    .D is the solid body boundary cells.

Model validation
The validation of our study is done in two steps, the first is to compare our results using the color gradient model with the case of two-phase immiscible fluids flow carried out by Huang et al. (Huang et al., 2013).The second step, we investigated the single-phase flow around a square obstacle using the color gradient model by giving the same density and viscosity values i.e.   =   and   =   by comparing our results with the work reported by Breuer et al. (Breuer et al., 2000).

Two immiscible layered fluids flow in a 2D channel
In this part, we simulate immiscible layered two-phase flow between two stationary parallel plates as shown in Fig. 2. The periodic boundary was applied for both left and right sides of the channel (Periodic flow).The bounce-back boundary conditions were selected for the upper and lower plates.The viscosities for both fluids are calculated as   =   2 ( − 0.5) and  =     ⁄ is the viscosity ratio.The wetting phase (Fluid 1) is placed in the region  < || ≤  and the non-wetting phase (Fluid 2) in the central region 0 ≤ || ≤ .At the inlet the velocity parabolic profil is considered.

Viscosity effect
The purpose of this test case is to verify the efficiency of our numerical code to simulate the immiscible layered fluids.The velocity profiles for identical fluid densities are presented in Fig. 3.The viscosity ratio is the only parameter that changed in this calculation.For the same values of  (1/ 50 and 50) and mesh size (10 × 100) used in the reference (Huang et al., 2013), we find the same curves of velocities.We note that our simulation errors between the analytical and numerical LBM model are: Fig. 3 (Huang et al., 2013).For both cases, the maximum velocity is localized in the center region of the channel.As the flow progresses inside this channel, it can also be shown that the shape of the velocity profile changes with the viscosity ratio.The error between numerical and analytical solutions is defined as: Here the summation is over the lattice nodes j and  0 is the analytical solution.The convergence criterion is:

Density effect and elimination of the extra unwanted term
Table 1 presents different cases used in this study with identical viscosities and different density ratios.Further, this test is performed to show the effect of the unwanted extra term that appears in equation ( 24) and the parameter  which adjusts the interfacial thickness.Fig. 4 shows the different velocity profiles obtained by the RK-Color gradient model in comparison with the analytical solution and the work of Huang et al. (Huang et al., 2013).In our study, the unwanted term plays an important role and its effect is clearer for different densities cases.We can observe that when the extra unwanted term is eliminated, the numerical solution agrees with the analytical solution with an estimated errors of 2.64% and 1.11% for both cases (b) and (d), respectively.Again, it agrees with the numerical solution of the reference with an estimated errors of 0.029% and 0.0098% for both cases (b) and (d), respectively.Hence, the parameter  is one of the important and free parameters in this simulation, it is used to adjust the thick of the interface between the two fluids which gives a good precision for an important difference density as illustrated in the reference (Huang et al., 2013).

Single phase flow in a channel with obstacle
In this study, as shown in Fig. 5, we consider a rectangular channel of length L and height H in which a square obstacle of dimension D = 40 (l.u) is placed in the middle of the horizontal walls of the channel of dimension 2000 × 320(l.u) 2 .The upper and lower walls of the channel are assumed to be solid (fixed plate) in which Bounce-back conditions have been adopted, the sides of the channel are assumed to be open and the conditions of Zou and He (Zou et al, 1997) are implemented.The flow is assumed to have maximum velocity   equal to 0.057.
In this part, our numerical code is used to simulate a single-phase flow through a channel containing an obstacle.This problem was studied in detail by (Breuer et al., 2000).Here, we present a comparison between our results and those of these authors.The Reynolds number   is the main parameter which influences the flow behavior, the numerical calculations are performed for   = 1, 30, 60 and 100.The computed results comparison is based on streamlines, drag coefficient and velocity profiles at different positions.It's noted that the different regimes flow is captured by our numerical code.Indeed, for   < 60, the regime is steady.In this case, the flow is unidirectional and uniform as shown in Fig. 6(a).Two counter rotating vortices appear symmetrically about the flow axis behind the square obstacle as shown in Fig. 6(c).From   = 60 (Fig. 6(e)), the flow becomes unsteady with the well-known von Karman vortex street with periodic vortex shedding from the cylinder (Admi et al. 2022(a), Admi et al. 2022(b)).For this flow regime, we are interested in the velocity following the main flow direction (  ) and therefore, we compared the velocity profiles to those of the reference.
Figures 7(a) and 7(b) illustrate the velocity components profiles   () and   () along a centerline  = /2 over x−direction for   = 100 and at a time instant at which the cross-stream velocity u at an axial position of 10D behind the cylinder changes its sign from minus to plus.We can show that before the obstacle, the component of the velocity   () is uniform along the center line until it is equal to zero at the solid.Just after this obstacle, it becomes negative, this is due to the disturbance caused by the obstacle and far from it, the flow tends to become regular accompanied by oscillations until it reaches its value at the entrance of the channel.However, for the orthoradial component   (), the curve shows the same behavior noted for   () before the obstacle.Whereas, after the obstacle,   () follows a wave variation and takes an amortissement from maximum to minimum values.
The last validated result is the drag coefficient   which is calculated on the obstacle and it depicted in Fig. 8.This coefficient gives information about the resisting force of a body immersed in a fluid environment that depends on flow direction, size, shape and placement of the body.We note that the drag coefficient decreases with increasing .This behavior is exactly similar to that noted in the reference (Breuer et al., 2000).Once again, we can conclude that for all the results found, graphical comparisons are made and show a good agreement between our results and those of the reference.Thus, we can say that our code is reliable and able to simulate single phase flows in channels with obstacles.This prompted us to study the same geometry with multi-component flow.In what follows, the calculations are carried out for  = 100.This value corresponds to a drag coefficient   = 1.378.

Physical problem studied
In the case of a single-phase flow in a channel with block, the increase in  produces Von Karman vortices behind this block (Moussaoui et al, 2009a(Moussaoui et al, , 2010b)).The amplitude of these oscillations increases with Re and consequently, the friction coefficient and the drag force increase on the walls of the channel.In some cases, this can cause technical problems.Thus, the aim of our study is to show how these vortices can be reduced.For this, we propose to replace the single-phase flow by two fluids flow.Indeed, we simulate immiscible layered fluids flows between two plates past a square obstacle in 2D channel with a wetting phase (high density) located between 0 ≤  ≤  and a non-wetting phase (low density) located between a ≤  ≤ .Fluids 1 and 2 are defined in red and blue, respectively.At the channel inlet and outlet, Zou and He (Zou et al, 1997) boundary conditions are applied.No-slip (bounce-back) boundary conditions are used at the solid boundaries.Midway the plates, a square obstacle of length D = H/8 is placed at distance  = /3 from inlet of the channel.The inlet velocity is parabolic with a maximum value   in the middle as shown in Fig. 9.
Around the obstacle, the boundary conditions used are bounce-back.For both distribution functions, the values of    (, ) = 0 inside the obstacle.The non-wetting condition is set by the parameter  = 10 −4 .The ratio of densities is  = 3 because these densities are adjusted by the parameters α1 = 0.8 and α2 = 0.4.It should be noted that the obtained results are validated in the case of a single phase by giving the same values of the densities and viscosities.

Mesh size sensitivity
The grid size is varied to find the best compromise between computational costs and accuracy.For this, the average drag coefficient around the obstacle (   ) is calculated for different meshes and fluid viscosities (  ).Table 2 illustrates these results and those of references (Breuer et al., 2000, UI-Islam et al., 2009, Sohanker et al., 1998, Benamour et al., 2015).It can be seen from this table that the mesh (1000 × 160) gives results which are very close to those of the literature and as   = 0.0887 is not considered here, this mesh is chosen to carry out the present study.

Density behavior over time t*
Fluid flow motion control is one of the major concerns of the mechanical fluids engineering.For this reason, we have proposed a method that allows us to control two immiscible fluids by adjusting their viscosities.At the beginning, the effect of viscosity on fluids flow after being distorted by the square obstacle is evaluated.Reynolds number is set to be  = 100 (  =   /  ).In a first calculation, the viscosity of fluid 1 is chosen to be twice the viscosity of fluid 2. Fig. 10 shows the hydrodynamic behavior of fluids from their densities for a double viscosity rate noted  =  1  2 ⁄ = 2 over time t* which describes the characteristic time, it is defined by t* = t/1000.As for a single phase and before reaching the obstacle, the flows are juxtaposed, stable and parallel to the two plates.However, after the obstacle, the fluids mix and vortices appear.Contrary to the case of a single fluid, these vortices decrease in frequency and amplitude over time and disappear around t* = 400 and the flow becomes steady.In a second calculation, the viscosity of fluid 1 is chosen such that it is half of the viscosity of the fluid 2. This case is presented in Fig. 11, from this figure, we note a very large difference compared to  = 2.In fact, the oscillations are strongly damped and disappear completely around t*= 200.The fluids flow becomes stable and steady for a very short time (t* = 120).Thus, it can be concluded that the fluids behavior and the disappearance of vortices strongly depend on the viscosity ratio M.

Effect of viscosity ratio on fluids flow
A short time (t* = 50) corresponding to strong oscillations is chosen to present the behaviors of fluids due to the change of the viscosity ratio .These results are illustrated in Fig. 12 and 13.It is noted that in the case where fluid 1 is less viscous than fluid 2 ( ≤ 1/3), the fluids circulate in parallel without overlapping over the entire channel.Thus, in this case, the obstacle does not disturb the flow.However, slight oscillations appear behind and further from this obstacle for  = 1/2.As  continues to increase and becomes greater than 2 i.e. fluid 1 becomes more viscous than fluid 2, these oscillations change shape and become larger indicating a large overlap between the fluids.
An unexpected result shows that the oscillations are damped with the further increase of M. This is due to the increase in viscous effects between the fluids.We notice that for the present value of t*, the oscillations disappear for  = 5, whereas for other values of t*, the behavior of the fluids changes.For this, we have defined the parameter h* = h/H which expresses the ratio between the height of the first oscillation of fluid 1 after the obstacle and the height of the channel.Then we have plotted its variations as a function of time in Fig. 14.This figure shows that this parameter decreases in amplitude until the zero-value indicating the total disappearance of the vibrations in a time which decreases with the increase of M for  > 2. Indeed, the oscillations disappear around 385, 250, 175 and 100 for  = 2, 3, 4 and 5, respectively.For  < 1/2 the opposite fluids behavior is noted.Thus, to have a more stable flow behind the obstacle, the difference between the two viscosities of the fluids must be very large.In fact, in the case of a single phase, the effect of viscosity decreases from the wall towards the central axis of the channel.This allows the Karman vortex to form and move a bit freely behind the obstacle.However, in the case of two fluids with different viscosities, at the interface which generally exists in the central zone of the channel, there is a braking phenomenon between the two fluids.This prevents vortex from forming and moving freely compared to the single fluid case.

Variation of velocity components as a function of the viscosity ratio 𝑀 at time t* = 50
The velocities   () and   () along a centerline  = /2 over x-direction for   = 0.1328 and for different values of M are depicted in Figs. 15 and 16, respectively.We can clearly notice that for all viscosity ratio M, the velocity component   () varies in the same way as in a single phase with a small perturbation at the canal inlet due to the friction between the two fluids which caused a reducing of their average velocity.For  ≤ 1/2, the velocity keeps almost the same shape for all values of M.However, for  > 2, The velocity component   ()decreases with increasing viscosity and tends to become linear as it moves away from the obstacle.Contrary to the single-phase case and at the channel outlet, the velocity not regains its initial value due to the viscous effect.From Fig. 16, it is very interesting to note that the oscillatory behavior of   () after the obstacle noted for a single fluid disappears in the case of two fluids of very high or very low viscosity ratio.

Effect of the viscosity on the drag coefficient
Figure 17 shows the evolution of average drag coefficient   versus the viscosity ratio  for different times spaced by a step of 100 except for the first two values.It should be noted that for all values of  ( 1 and  1),   decreases with time.For  1,   becomes almost constant for a time which exceeds that which corresponds to the disappearance of the oscillations.Also, when the viscosity of the fluid 1 increases ( increases), the drag force exerted by this fluid on the obstacle increases and consequently   increases.However, for  1, this last remark remains valid when the viscosity of the fluid 2 increases ( decreases).

Conclusion
Numerical computations were performed to simulate the two layered immiscible fluids flow in 2D channel past a square cylinder obstacle.The RK color-gradient model is used.The main objective of this article was to investigate the dynamic behaviors of the two fluids through the obstacle by varying their viscosity ratio in unsteady regime ( = 100), and to verify the important effect of viscosity on flow control.The results were presented in terms of figures showing the behavior of density and components of velocity over time.From these results, we can draw the following conclusions: -To have stable, parallel and non-overlapping flows behind the obstacle, it is necessary that the difference between the viscosities of the fluids be significant i.e.  = 5 or  = 1/3.-For  ≥ 2, the increase in the viscosity difference leads to an increasing of friction between fluids and a reducing of average velocity of flow.-For  ≥ 2, the increase in the viscosity ratio decreases the time corresponding to the disappearance of the vortices behind the obstacle.However, for  ≤ 1/2, the opposite occurs.

-
The average drag coefficient increases with viscosity ratio for  ≥ 2. These results show that the effect of the obstacle on the fluids flow can be eliminated by varying the viscosity ratio.Indeed, the vortices caused by the obstacle can be reduced or eliminated.At the same time, a good compromise must be sought between this objective and the drag coefficient which increases with the viscosity ratio.In addition, the numerical results obtained for all cases show that this method can be used for other flows with multiple phases and components in 3D.This will be the subject our future work.

Fig. 1
Fig. 1 Illustration of a lattice node of the D2Q9 Model.

Fig. 2
Fig. 2 Illustration of the computational physical problem.
(a): err = 2.26%, Fig. 3(b): err = 3.9% and the errors between our numerical results and those of the reference are: Fig. 3(a): err = 1.89%,Fig. 3(b): err = 1.28 %.Then, we can conclude that our results are in good agreement with the analytical solution and with the work of Huang et al.

Fig. 3
Fig. 3 Velocity profiles for the flow of immiscible layered fluids of identical densities   =   .

Fig. 4 Fig 5
Fig. 4 Compared results with those of Huang et al. (Huang et al., 2013) and the analytical solution ((a), (b), (c) and (d) are the cases of Table1).
Fig. 9 Illustration of the computational physical problem.

Fig 12 .
Fig 12. Behaviours of fluids for a viscosity ratio ranging from 2 to 5 at t* = 50.
Fig. 14 Variation of characteristic height h* as a function of t* and .

Table 2
Average drag coefficient as a function of mesh size for  = 100.