Numerical Investigation of Flow Patterns and Plug Hydrodynamics in a 3D T-junction Microchannel

A three-dimensional (3D) numerical simulation of a two-phase flow liquid/liquid is performed in a rectangular microchannel with a T-junction. The volume of fluid (VOF) method was used under ANSYS Fluent to capture the interface between the two phases. The dynamic mesh adaptation technique together with the assumption of symmetry plane helps us to reduce the computational cost. The study focuses on the flow patterns and hydrodynamics of plugs. So, the influence of the flow rate ratio q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q$$\end{document}, the capillary number Ca\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Ca$$\end{document}, and the viscosity ratio on the liquid film, the plug/droplet shape, and velocity are examined here. Particularly, the plug/droplet lengths predicted by the simulation show good agreement with the experimental and correlation available in the literature. The results revealed six distinct flow patterns by dispersing water in a continuous phase of silicone oil. By decreasing the flow rate ratio as well as the viscosity ratio, the liquid film thickness increases in the corners and side planes. In turn, this greatly impacts the liquid film velocity and the plug velocity. Furthermore, capillary number (based on two-phase flow velocity) is also shown to have a greater impact than viscosity ratio and flow rate ratio on plug shape, with the curvature radii of the tail always larger than the front one.


Introduction
In recent years, researchers from various areas of science have become interested in two-phase flow, especially the flow of immiscible liquids in microchannels. The high surface-to-volume ratio, highly uniform size, and highprecision control of flow rate make microchannels an attractive tool for many technological processes (Kashid et al. 2011) (Murphy et al. 2018;Qian et al. 2019a, b).
Liquid-liquid two-phase microfluidics is particularly the most important subject of research in this field (Oishi et al. 2009;Qiao et al. 2012;Song et al. 2006). Furthermore, it is widely used in a variety of applications, including extraction, nitration, and nanoparticle synthesis. These processes benefit from either their favorable mass transfer characteristics or their droplet formation that allow the encapsulation of fluids.
For such applications, the precise and efficient generation and control of droplets are essential. For this purpose, droplet generation methods can be classified as either passive or active.
Passive techniques are generally more economical, easier to implement, and generate droplets without external actuation. The purpose of these techniques is to identify the mechanisms and characteristics of droplet breakup modes that can occur in microfluidic cross-flow (e.g. T-junction or Y-junction) (Muijlwijk et al. 2016;Nisisako and Hatsuzawa 2010;Paquin et al. 2015), co-flow , flow-focusing (Fu et al. 2012), and step emulsification configurations (Chakraborty et al. 2017).
In contrast, active techniques use additional energy input by promoting interfacial instabilities for droplet generation. Generally, these methods either use external forces from electric (Gompper and Schick 2007), magnetic , and centrifugal fields, or methods that modify intrinsic properties of flows or fluids.
A T-junction is one of the most commonly used geometries for creating droplets. Thorsen et al. (Thorsen et al. 2001) presented it for the first time as a way to produce oil droplets in water. In this configuration, the continuous phase enters directly into the main channel, whereas the dispersed phase is fed perpendicularly through the secondary channel. The drop generation mechanism was described earlier by Garstecki et al. (2006) and different regimes can be observed, depending on the capillary number Ca c .
In most studies, the effect of the fluid viscosity of the continuous phase is examined. Tice and al (Tice et al. 2004) investigated the formation of plugs in multiphase flows of fluids more viscous than water in a rectangular microchannel with a hydraulic diameter D h = 50 m . They observed that the capillarity number, Ca, can be used to predict the conditions required for plugs to form reliably when the fluids in the flow have approximately the same viscosity. In contrast, when there is difference in viscosity between the two fluids (continuous phase and dispersed phase), threshold values of Ca differed by a factor of 1.5 − 3.4 from those of fluids with matched viscosities.
Based on a 2D microchannel T-junction with the same entry width ((W c = W d = 100 m) for the two phases), Liu and al (Liu and Zhang 2009) investigate numerically the role of viscosity ratio in the droplet break-up process. They found that the viscosity ratio λ has a more significant effect at high Ca ( Ca > Ca c ) where the large viscosity ratio forces the droplet detachment point to move downstream. They also found that the predicted drop diameter in the squeezing regime is practically independent of the viscosity ratio, however, in the dripping regime, the influence of the viscosity ratio becomes more obvious as Ca increases.
At the same period, Gupta and al (Gupta and Kumar 2010) in a three-dimensional simulation of two immiscible fluids through a T-junction revealed that the transition to stable jets can be delayed and droplets can be formed even at very high flow rate ratios Q by decreasing significantly the viscosity ratio . Furthermore, for a fixed Ca and Q , they mentioned that the droplet formation frequency can be increased by increasing the viscosity ratio . The same results were highlighted by Wehking and al (Wehking et al. 2014) who reported specific characteristics during droplet generation for the different regimes (Droplet regime in a T-junction "DTJ", Droplet regime in the channel "DC", and Parallel flow "PF"). Particularly, the authors confirm the earlier observation of Gupta and al (Gupta and Kumar 2010) by finding that the transition region from DC to PF shifts to a high flow rate ratio due to the high shear stress of the continuous phase on the dispersed phase at low viscosity ratio, which facilitates droplet formation.
Droplet formation involving highly viscous fluids is studied experimentally in a square microchannel T-junction by Bai and al (Bai et al. 2016). They noticed that in the dripping regime, the effect of high viscosity c of the continuous phase on the droplet size is so significant that Ca is not enough to correlate droplet size. Besides, they noticed that the critical Ca number increases with the continuous phase viscosity, which shows that there is no single critical Ca value to identify the various flow patterns when various combinations of fluids with various viscosity ratios λ are involved. Characteristics of liquid-liquid slug flow in a microchannel are investigated experimentally by Yao and al (Yao et al. 2018) in a square microchannel with a hydraulic diameter of D h = 600um . Their results reveal that both the continuous and dispersed phase viscosities affect the formation of droplets, though the latter has a more profound impact. Additionally, low channel wetting leads to an earlier transition from squeezing to a dripping regime (i.e., a smaller Ca number).
The effect of viscosity ratio on the size of droplets is considered in the 3D numerical simulations by Mahdi and al (Nekouei and Vanapalli 2017). They used the volumeof-fluid (VOF) method to perform droplet formation of Newtonian fluids in a microfluidic rectangular T-junction device with width and height ratios set to W = 0.5 and H = 0.33 respectively. In the range of the viscosity ratio 0.01 < < 15 , the researchers found that increasing the Ca number decreases the droplet size. However, the reduction in the droplet size with the capillary number is stronger for < 1 than for > 1 . Furthermore, they announced that at a particular capillary number Ca, the droplet size is not significantly influenced for < 1 , but it increases when > 1.
Li and Angeli (Li and Angeli 2017) explored the hydrodynamics of the liquid-liquid plug through an experimental and numerical study in circular microchannels of 0.2 and 0.5 mm in diameter. According to their results, the liquid film thickness does not increase linearly with capillary, and for the same Ca, it is slightly lower in the 0.2 mm channel than in the 0.5 mm one. Additionally, the researchers noticed a change in shape for the front and rear menisci of the plug with high Ca . The size of the channel, however, has no significant effect on the shape of the plug.
Recently, Kovalev et al. (2018) present an experimental study of the liquid-liquid system with an extremely low viscosity ratio ( = 10 −3 ) in a T-shaped microchannel. The researchers observed new flow patterns including plug, droplet, plug with suspended droplets, slug, throatannular and parallel flow and that the transition between the flow regimes is mainly described by dispersed phase Reynolds number. The authors used a typical parameter (Q d ∕Q c × Ca bulk ) to describe the plug interface geometry and flow structure. An increase of this parameter leads to a transition from the typical Taylor bubble shape of the plug to a dumbbell-like plug shape. Furthermore, the velocity distributions, measured by Particle Tracing Velocimetry (PTV) technique, inside plugs are essentially threedimensional in contrast to other previous studies reported in the literature. Finally, they found that the curvature radii of the plug at the rear decreases linearly with increasing Ca , whereas the curvature radii at the front is almost constant in contrast with Li and Angeli (Li and Angeli 2017) results.
More recently, Kovalev and al (Kovalev et al. 2021) conducted a more in-depth study of the influence of viscosity ratio on liquid-liquid two-phase flow in a T-shaped microchannel with 200 × 200 m inlets and a 400 × 200 m outlet channel where the two fluids enter in opposition. The castor oil was used as a continuous phase while distilled water and water-glycerine solutions as the dispersed phase. In the range of these parameters, the researchers observed three basic flow patterns and pointed out that the boundary between segmented (plug and droplet) and continuous (parallel) types of flow was found to shift depending on . Particularly in the plug flow case, at low viscosity ratios, = 0.0013 and = 0.0067 the plug interface shape deforms downstream of the flow which requires a channel length downstream of the junction of about 38 × w out . Finally, the researchers confirm the observation mentioned earlier by Liu et al. (Liu and Zhang 2009) about the effect of viscosity ratio in the squeezing and dripping regimes.
In the experimental investigation of Ma and al , aqueous glycerol solutions and a mixture of toluene and dimethyl silicone oil were chosen to vary the viscosity of the continuous and dispersed phases respectively within a rectangular microchannel. The researchers confirm the overall observations made by Kovalev and al (Kovalev et al. 2021) and mentioned the shift of the slug-drop flow transition to be at a higher capillary number of the continuous phase (Ca c ) while the slug-parallel flow transition to a higher Reynolds number of the dispersed phase (Re d ) . The authors indicated that for a lower Ca number based on the plug velocity ( Ca D < 0.03 ), the film is thinner in liquid-liquid flow than in gas-liquid flow. There is a distinct plug shape characteristic under high Ca D that is marked by both sharpened caps with a more sharpened cap at the rear end. Besides, the researchers revealed that the shape of the plug deforms slowly under squeezing regimes when compared to dripping regimes. Wang et al. (2021) studied the liquid-liquid two-phase flow patterns and hydrodynamics of plugs in a circular channel of 2mm diameter. Silicone oil and water are used as working fluids. Their results present five distinct flow patterns. By increasing the capillarity number based on the two-phase flow, Ca TP , the curvature radius of the plug head becomes smaller, while the tail plug one increases gradually with Ca TP . Besides, the researchers revealed, at low Ca TP , the appearance of a concave surface on the lateral side of the plug which slowly relaxes until it completely disappears for high Ca TP . For predicting the plug length, the researchers proposed a correlation based on the Reynolds number Re TP (based on the velocity of the two phases) and the flow rate ratio. Satisfactory results were reported using this correlation compared to the experimental data as well as to the data available in the literature. From what has been reviewed in the literature, most of the available numerical studies on liquid-liquid two-phase flow have been carried out in microchannels of width or diameter less than 300 m . Furthermore, there were no numerical studies that address the flow patterns and hydrodynamics of liquid-liquid plugs in rectangular microchannels with fluid combinations that provide a low viscosity ratio. To this end, the present work aims to investigate the effect of parameters such as Ca TP , flow rate ratio Q , and viscosity ratio λ on the hydrodynamics of plugs in a rectangular microchannel with a cross-section of 600 × 300 m . The liquid film, plug shape, and plug velocity will be discussed as well as the flow patterns.

Geometrical Model and Boundary Conditions
A 3D T-junction rectangular microchannel was chosen to study the flow patterns and hydrodynamics of liquid-liquid plugs in various conditions. The schematic diagram of the microchannel geometry with boundary conditions is shown in Fig. 1. The geometry and dimensions are the same as those of Ma and al . It is a microchannel T-junction with a width W c = W d = 600 m and a height of h = 300 m . To minimize computational costs, we took a smaller computational domain than the whole channel. However, the length of inlet channels is fixed to be sufficiently long to ensure fully developed flow (Ookawara et al. 2006;Reddy Cherlo et al. 2010). This length is determined by the correlation introduced by Dombrowski and al (Dombrowski et al. 1993).
D h is the hydraulic diameter ( D h = 400 m in this study) and the entrance length L e of the inlet channels is consequently calculated for different Reynolds.
A uniform velocity is established at both inlets, while at the outlet the outflow condition is used. With this boundary condition combination, mass conservation is provided for liquid systems (Reddy Cherlo et al. 2010). Notice that Raj and al. (Raj et al. 2010) found that using pressure outlet or outflow as the outlet boundary condition did not affect predicted slug length. In order to eliminate the effects of wall wettability on plug formation, the contact angle is set to 180 • . A non-slip boundary condition is used on all the walls.
To further reduce the computational cost, a symmetry along the XY plane (see Fig. 1) is considered.

Mathematical Model
To capture the different positions of the interface between phases, two main approaches are commonly used: Lagrangian and Eulerian methods. If the interface deformations are very significant, Eulerian methods are generally the most efficient. Three models including VOF (Malekzadeh and Roohi 2015;Mastiani et al. 2017;Khan et al. 2018;Qian et al. 2019a, b), Lattice Boltzmann (Ba et al. 2015), Level set (Wong et al. 2017) and Coupled Level set and volume of fluid (CLSVOF) (Chakraborty et al. 2019), are mainly adopted. The volume of fluid method (VOF) introduced by Nichols and Hirt (Hirt and Nichols 1981) is one of the multiphase models implemented in ANSYS Fluent and used in this simulation. It is based on the fact that two or more fluids (or phases) are immiscible. The total volume fraction of all phases in each control volume is equal to one. In addition, all variables shared the volume fraction of each phase and represent volume-average value, as long as the volume fraction of each phase is known within a calculation cell (Canonsburg 2012).

Governing Equations
The isothermal two-phase flow of Newtonian fluids of constant density was described using the one-fluid approach.
The continuity and momentum equations are as follows: The volume fraction φ allows tracking the phase: it equals zero in the continuous phase, and one in the dispersed phase. (1)

Fig. 1 Schematic diagram and meshing of the T-junction microchannel
So, the physical properties such as density and dynamic viscosity are expressed as the volume-fraction weighted average between the two fluids.
Where the indices c and d denote the continuous and dispersed phase, respectively. Finally, the continuity equation for the volume fraction function is expressed as follows: Surface tension is treated as a source term in the momentum equation using the continuum surface force CSF model (Brackbill et al. 1992). The surface tension force per unit volume F in Eq. (3) was calculated as: Where is the surface tension, n is the interface normal vector defined as the gradient of the volume fraction, i is the Dirac delta function specifying the interface, is the interface curvature evaluated from the divergence of the unit normal surface n.

Working Fluids
In this investigation, a variety of fluid combinations were employed with various viscosity ratios λ. We used silicone oil or glycerol/water mixture as continuous phase, and toluene or water as a dispersed phase. Tables 1 and 2 summarized details of the combinations.
These fluid properties are identical to those used in previous studies by Ma et al. (Ma et al. 2021). and Wang et al. .

Parameter Settings
The simulation parameters used in this investigation are given in Table 3.

Mesh Grid Independence
A structured mesh was generated using ICEM CFD and different mesh sizes were used in order to achieve grid independence. For the third combination (C 3 ), where the superficial velocity of the dispersed phase was set to 0.0185m • s −1 and the continuous phase velocity was fixed at 0.02763m • s −1 , four mesh size was chosen, 2 × 10 −5 m, 1.7143 × 10 −5 m , 1.5 × 10 −5 m and 1.33 × 10 −5 m respectively. In order to further reduce the computational cost, a dynamic mesh adaption based on the gradient of the volume fraction was used. The mesh adaption parameter, which includes the refinement and coarsening thresholds, is fixed equal to 0.15 and 0.1, similar to Mehdizadeh and al (Mehdizadeh et al. 2011) and Alireza and al (Alireza et al. 2018) and the mesh is generated every five-time step.
The results showed that the four basic meshes used gave almost identical results in terms of droplet length, detachment time (Table 4), and static pressure in the centerline along the channel (Fig. 2), and that differences could be attributed to the sharpness of the interface between the two phases. Finally, we chose the mesh that gave us the appropriate results, a sharp interface, and an acceptable computational time.

Validation
The results of the simulations are validated by comparing the dimensionless droplet length with the experimental  results and correlation performed by Ma et al. (2021), in a range of flow rate ratios varying from 0.5 to 2. Throughout the validation, the flow rate of the dispersed phase was kept constant Q d = 0.2ml∕min while that of the continuous phase was changed. Figure 3 shows the present simulations for two contact angles ( 158 • ;180 • ) in terms of dimensionless plug/droplet length (L∕W) compared with the experimental results of Ma et al. (2021) and their correlation.
As can be seen, the present simulation results are in good agreement with the experimental data as well as the correlation proposed by Ma et al. (2021).
For a flow rate ratio q ≤ 1 , there is no clear difference in plug/droplet length when moving from a contact angle of 158 • to180 • . Furthermore, the simulation results deviate slightly from the experimental data.
On the other hand, for q > 1 , the contact angle starts to affect the plug/droplet length progressively with increasing q . Besides, the deviation between the simulation and experimental results gradually decreases until it becomes negligible for q = 2.
Generally, and for a partial wetting of the dispersed phase with the channel walls same as , The average relative error between the simulation and the experimental data does not exceed 7% . This indicates that the VOF method can be used to correctly simulate the liquid-liquid two-phase flow, even with a high viscosity ratio λ.

Flow Patterns
The fourth fluid combination (C 4 ) ( Table 1) with a viscosity ratio λ = 0.0086 is considered to analyse the flow patterns for a liquid/liquid system. A contact angle of 180 • was used to overcome wettability effects.
According to Tsaoulidis et al. (2013), the flow patterns are not affected by the channel configuration (T or Y junction). However, flow parameters such as the flow rate of the two phases, the properties of the two fluids, or even the wettability of the walls can strongly influence the flow patterns.
As can be seen in Fig. 4, six flow patterns were observed in this study: jet flow, droplet flow, plug flow, throat annular flow, annular flow, and stratified flow.
The appearance of flow patterns results from competition between the forces present at the junction. These forces are   . In contrast, for gravitational effects, the Bond number is very small ( Bo = 1.55 × 10 −3 ), which means that gravity effects are neglected in this study. When the surface tension force is the dominant force compared to the viscous and inertia of both phases, regular flow patterns (Plug and droplet) are generated (Darekar et al. 2017;Guan et al. 2019).
Plug flow (Fig. 4c) is observed at low and comparable flow rates of both phases. In this flow pattern, plugs are usually produced under the squeezing or shearing regime in a microchannel, with a length usually equal to or greater than the width. Under certain conditions, a liquid film of the continuous phase appears, separating the plug from the channel walls. Besides, the plugs preserve a defined shape after the flow has fully developed.
By increasing the flow rate Q c of the continuous phase while keeping the dispersed phase flow rate fixed, the droplets pattern was observed (Fig. 4b). The main difference between droplet flow and plug flow is the size of the droplets and the detachment mechanism.
In this flow pattern, the Reynolds number of the dispersed phase Re d remains small and the dispersed phase was easily broken into droplets due to the high shear exerted by the continuous phase (Re c > 0.4) . This leads to a short detachment time and droplets with a length smaller than the channel width.
By further increasing the flow rate of the continuous phase Q c while keeping that of the dispersed phase fixed, jet flow (Fig. 4a) was formed. This pattern is dominated by the strong shear and drag forces of the continuous phase (Lei et al. 2021).
In this pattern, the dispersed phase remains stuck to the channel wall and pulled away from the junction due to the drag force exerted by the continuous phase. Then, the high shear of the continuous phase causes the dispersed phase to break up rapidly into small droplets (Lei et al. 2021). These droplets are usually much smaller than the width of the microchannel.
When the flow rate ratio becomes large with an increase in the flow of the dispersed phase, the throat annular flow was observed. In this flow pattern, the inertial force of the dispersed phase becomes so great ( Re d = 66.4 ) that the continuous phase is no longer able to break it up into plugs.
According to Wang et al. (Wang et al. 2021), the interface shape and periodic waves observed in this flow pattern are mainly caused by surface tension effects in the junction region.
The same flow pattern was observed in the study of Liu et al. ) with a large viscosity of the continuous phase. Yagodnitsyna et al. (Yagodnitsyna et al. 2016) observed the throat annular flow with a partial wetting of both fluids with the channel walls. However, a contact angle of 180 • is set in this study to eliminate the effects of wettability.
If the dispersed phase flow rate continues to increase while the flow rate of the continuous phase remains low, the flow pattern changes to an annular flow. In this flow pattern, the inertial force of the dispersed phase ( Re d = 132.8 ) can no longer be ignored and has an important role in the appearance of the annular flow by overcoming the surface tension force.
As seen in (Fig. 4e), the dispersed phase (water) occupies the middle of the channel continuously, while the continuous phase (Silicone oil) is spread out at the sides of the channel, forming a thin liquid film that wets the channel walls. It should be noted that this flow pattern was observed far from the junction. Therefore, a 33 mm channel was used to observe it. At high flow rates for both phases, a stratified flow can be observed. From the Re numbers of both phases ( Re c = 0.9;Re d = 121 ), it is clear that the inertial force of the dispersed phase overcomes the viscous and surface tension forces.
From Fig. 4f, it can be seen that the dispersed phase moves in the middle of the channel with a non-regular wavy interface. This wavy interface is caused by shear instability as a result of the difference in velocities between the two phases .
This stratified flow is similar to that of Wang et al. . However, the one of Wang et al. remains attached to the channel wall, which is not the case in our simulations.
This may be due to the effects of the contact angle which is one of the important parameters affecting the flow patterns in small channels (Barajas and Panton 1993;Tsaoulidis et al. 2013). According to (Kovalev et al. 2018;Yagodnitsyna et al. 2016), a high contact angle of the dispersed phase with the channel walls favours the development of a stratified flow. In such a case, the dispersed phase is easily transported by the high shear of the continuous phase .
The competition between viscous and inertial effects at the junction is a major factor in the formation of the plug/ droplet regime. In our simulations, the plug/droplet was formed in a transition regime (between the squeezing and dripping regime). The pressure built up upstream of the junction along with the viscous shear exerted by the continuous phase on the emerging plug/droplet contributes to the formation mechanism. On the other hand, the inertial force of the dispersed phase resists detachment, thus giving the plugs/droplets a certain length ).
Different flow rates of the two phases were used in our simulations, leading to a change in the period of plug/droplet formation. This is mainly due to the change in the competitive relationship between the viscous and inertial forces at the junction. Figure 5 shows the time period of plugs/droplets T P measured and estimated as a function of flow rate ratio q . The flow rate of the dispersed phase was constant in this section (Q d = 0.2ml∕min) . In addition, two combinations of the fluid were used for two different viscosity ratios λ (see Table 1).
The time period of plugs/droplets T p was estimated by measuring the time interval required for a unit cell (plug and slug) to pass through a cross-section. Thus, our results were compared to the time period calculated using the following equation: As can be seen, the time period T P increases with increasing q . This indicates the resistance of the dispersed phase to (10) T P = L D + L S U D detachment as q increases, resulting in longer plugs. Besides, the estimated and measured T P are in good agreement for both viscosity ratios λ, with an average relative error of no more than 2% for both λ.

Effect of Ca Number
A series of simulations were performed using silicone oil as the continuous phase and water as the dispersed phase (see Table 1). Along this section, the viscosity ratio is set to = 0.0086 and a fixed flow rate ( q = 1 ). The capillary number is calculated based on the mixing velocity U TP . So, we chose Ca TP = 0.054, 0.11, 0.21 and 0.43 to study the effect of Ca number on particularly: the liquid film in the corners Corner and at the side walls y , the plug shape, the plug velocity.

Film Thickness
The liquid film separating the dispersed phase from the microchannel wall and the liquid film on the corners are two important parameters affecting heat and mass transfer in the microfluidic systems .
Currently, most studies have focused on the liquid film in gas/liquid systems in circular, rectangular and square microchannels (Aussillous and Quere 2000;Han and Shikazono 2009;Kreutzer et al. 2005;Yao et al. 2019). However, few experimental and numerical studies are available in the literature that concern the liquid film in liquid/liquid systems with square or rectangular microchannels.
It should be noted that the calculus results, for plug lengths greater than 1.57 ( L P ∕W > 1.57 ) (see Fig. 6), show a uniform liquid film area. For this purpose, we considered a cross-section in the middle of the plug to calculate the liquid film thickness in the sidewall ( y ) and in the corner ( Corner ). However, for plug lengths less than 1.57 ( L P ∕W < 1.57 ) (see Fig. 6), the liquid film varies along the plug (non-uniform). Therefore, the liquid film thickness is obtained by averaging the film thickness in 5 sections along the plug. The same observation was noted by Li et Angeli (2017).
The Bond number Bo is less than 0.022 in this study (Bo ≪ 1) which means that the effects of gravity are neglected at this scale. Therefore, viscosity, surface tension, and inertia are the only forces governing the flow in general and the liquid film thickness in particular.
The evolution of the dimensionless liquid film ( ∕W ) as a function of Ca TP is represented in Figs. 7 and 8.
From Fig. 7, it can be seen that increasing the mixture velocity U TP is traduced in an increase in the liquid film in the corners and sidewalls of the microchannel. Increasing the mixture velocity leads to an increase in the viscous stress of the continuous phase and will also increase the inertial force of the dispersed phase. It's the competition between these two forces which define the thickness and shape of the liquid film in the corners Corner and side walls y .
As can be seen in Fig. 8, for U TP = 0.0185m∕s , a thin and uniform liquid film surrounds the plug. As the mixture velocity increases, viscous and inertial stress start to affect the shape and thickness of the liquid film ( y , Corner ). So, at U TP = 0.037m∕s , the liquid film becomes thick, non-uniform and a non-negligible flow of the continuous phase around the plug appears.

Plug Shape
In this section, we discuss the effect of U TP on the interface shape at the front and tail of the plug as well as the sidewall. As previously stated, the flow ratio is kept constant, q = 1.
The U TP and the flow rate ratio q are the most important parameters affecting the shape of the plug (Kovalev et al. 2018).
The dimensionless curvature radii R tail and R front are shown in Fig. 9. Globally, we notice that R tail and R front decrease with increasing Ca TP under a fixed flow ratio (q = 1).
Furthermore, we can see that the curvature radius of the plug tail is always greater than that of the front. The shape of the plugs is affected by the viscous and inertial forces as the only factor that we change here is U TP . This result agrees with

Ma et al. (2021) explanation that increasing U c increases both viscous stresses (Ca c ) and inertial stresses (Re c ).
For Ca TP < 0.2 , the back of the plugs is more affected than the front. However, for Ca TP > 0.2 , there is a minor difference between the front and the back of the plug.
As like mentioned in Fig. 8, for a low U TP (low Ca TP ), there is a concave surface on the lateral side of the plug surface. Wang et al. (2021) explain this by the fact that the interfacial force causes the interface to shrink which leads to the formation of a concave surface on the lateral side of the plug.
By increasing U TP ( Ca TP increases) the concave surface disappears and the plugs/droplets take on an oval shape at high Ca TP .

Plug Velocity
The plug velocity U D versus mixture velocity U TP is shown in Fig. 10 for a constant flow ratio (q = 1).
Two methods are used to determine the plug/droplet velocity ( U D ) in this section.
The first method entails the displacement of a plug/droplet through a cross-section. This cross-section is chosen away from the junction, thus ensuring a fully developed flow. Hence, the plug velocity is derived by dividing the plug/droplet length by the time required for it to pass through the cross-section.
The second method consists of averaging the bulk velocity obtained inside the plugs/droplets for different planes (seven) along the z-direction.
The bulk velocity in each plane is determined by dividing the flow rate obtained inside the plug/droplet by its width.
As can be seen, the U D (based on the first method) increases linearly with U TP . Moreover, U D is always higher than the U TP . This is normal as the film surrounding the plug increases with increasing U TP as can be seen in Fig. 8. This leads to a gap between U D and U TP . The same observations were reported by Li et al. (2017).
The present simulation results are compared to the literature correlations listed in Table 5.
Globally the results are in good agreement with the correlation of Liu et al. (2005) with a mean error of 14% . For the second correlation proposed by Kovalev et al. (Kovalev et al. 2018), for low viscosity ratio (λ = 0.001) liquid/liquid systems in a rectangular channel, an average relative error of 11% is observed for U TP < 0.12 range. Figure 10 also presents the plug velocity U D determined from the bulk velocity inside the plug (second method). It can be seen that the two methods have a very good agreement with an average relative error that doesn't exceed 6%.
In the following sections, the U D obtained from the first method is used to calculate the Capillary Number Ca D .

Effect of Flow Rate Ratio q
In this section, we analyse the effect of flow rate ratio q on the liquid film, droplet/plug shape, and velocity. We fixed the flow rate of the dispersed phase (Q d = 0.2ml∕min) , while   (Kovalev et al. 2018) that of the continuous phase is varied. Also, two viscosity ratios λ ( = 0.0127 and = 0.0086) are considered in this section. For each λ, different flow rate ratios are used.

Film Thickness
Figures 11 and 12 present the evolution of the dimensionless liquid film thickness ( Corner , y ) as a function of the Ca D number (based on the plug velocity) for two viscosity ratios λ. The figures also present the predictions obtained from the correlations available in the literature (see Tables 6 and 7). For Ca D < 0.12 , the liquid film thickness obtained in the corner ( Corner ) with our simulation seems to be in good agreement with the correlation proposed by Kreutzer et al. (Kreutzer et al. 2005) for a gas/liquid system in square microchannels. As Ca D increases, the gap between our simulation and the correlation also increases with a deviation between 4 and 26% ( 15% on average) for λ = 0.0086 and between 1 and 19% ( 6% on average) for λ = 0.0127.
However, when comparing our results with the correlation of Yao et al. (Yao et al. 2019) and the correlation of Ma et al. (2021), Corner is overestimated by both correlations. The deviation between our results and the Yao et al. correlation doesn't exceed 15% for (λ = 0.0086) and 20% (λ = 0.0127) and 28% for (λ = 0.0086) and 27% (λ = 0.0127) for the Ma et al. correlation. This discrepancy between our results and the correlation of Ma et al. may be attributed to the partial wetting of both phases with the channel walls used in the Ma et al. study (2021), which may affect the liquid film thickness ( Corner and y ).
Aussillous and Quéré (Aussillous and Quere 2000) suggest that the liquid film at high Ca (high velocity) is governed by the competition between inertial and viscous forces. This implies the presence of the dimensionless numbers as Re, We, Ca in the correlations that predict the liquid film. To date, no correlation has been able to predict Corner accurately for high Ca D and low viscosity ratios .
The liquid film thickness y evolution versus λ, as can be seen in Figs. 11 and 12, agrees quite well with the correlation proposed by Yao et al. (Yao et al. 2019) with a deviation between 4 and 15% ( 10% on average) for λ = 0.0086 and between 2 and 34% ( 9% on average) for λ = 0.0127.
In contrast, we can see that y deviate from the correlations proposed by Aussillous et Quéré (Aussillous and Quere 2000) and Irandoust et Andersson (Irandoust and Andersson 1989) as Ca D increases.
It's worth mentioning that both correlations are based on experimental results for gas/liquid systems in a circular microchannel. Furthermore, the y predicted by these two models ignores the inertia effects (Re or We) at high velocity and is based only on the competition between the surface tension force and the viscous force (Ca) . This would explain the discrepancy with our simulations and indicates that both models are adequate for predicting the liquid film at low Ca (Li and Angeli 2017). Han and Shikazono (2009) in a gas/liquid system propose a correlation that takes into account the effects of inertia, viscosity, and surface tension (see Table 7). However, this correlation underestimates the liquid film at high capillary (high velocity) by 28% → 48% ( 36% on average) and by 8% → 70%(37% on average) for the two cases λ = 0.0086 and λ = 0.0127 respectively. This is probably due to the viscosity of the two phases (Li and Angeli 2017). Moreover, the shape of the channel along with the competition between forces (interfacial, viscous, inertia) at the junction may also contribute to this discrepancy. Using experimental results obtained with a liquid/liquid system in circular microchannels, Li et Angeli (2017) propose a correlation to predict the liquid film. This correlation is based on the Ca number to highlight the effects of the interfacial tension and viscosity, and on the We number to incorporate inertia effects. However, there is poor agreement between our results and y predicted by this correlation. The deviation is between 9 and 90% ( 52% on average) for λ = 0.0086 while it is about 9% and 95%(57% on average) for λ = 0.0127. One possible reason for this discrepancy is the channel size and shape which plays an important role in the behaviour of the liquid film in the liquid/liquid system. Li and Angeli (2017) noted an average difference of 13.66% in liquid film thickness between the 0.2 mm and 0.5 mm microchannel. Figure 13 shows the dimensionless liquid film as a function of flow rate ratio q . The flow rate of the dispersed phase Q d is kept constant.
As can be seen in this Fig. 13, the dimensionless liquid film decreases as q increases. At low q (high Ca c ) the liquid film ( Cr , y ) is the result of competition between the viscous effects of the continuous phase and the inertia effects of the dispersed phase. By increasing q ( Ca c decreases), the viscous effects of the continuous phase decrease and the liquid film result from the competition between the surface tension force of the continuous phase and the inertia force of the dispersed phase.

Plug Shape
For two viscosity ratios and different flow rate ratios, the shape of the plug/droplet is shown in Figs. 14 and 15.
Under the influence of the viscous effects of the continuous phase, a typical plug/droplet is observed (i.e. R tail of the plug/ droplet greater than the R front ). Kovalev et al. (2018) and Ma et al. 2021) reported a slender plug/droplet under high Ca with a sharper tail of the plug/droplet than the front. It is noted that both studies used partial wetting of the dispersed phase with the microchannel walls, which is probably the main reason for the observed trend.
Furthermore, for λ = 0.0127, a concave surface in the lateral sides of the interface is observed for a high flow ratio (low Ca ). The interfacial tension force is behind this shrinkage of the interface . By decreasing q ( Ca increases) the interface in the lateral sides becomes flat.
The front and rear curvature radii of of the plug/droplet ( R tail andR front ) dimensionalized by the channel width, are shown in Fig. 16.
As can be seen, the tail curvature radii of the plug (R tail ) increases with increasing flow rate ratio for both viscosity ratios. Similarly, the curvature radii of the front (R front ) also increases with increasing q . However, for = 0.0086 , there is not much difference in the front curvature radii when we  move from q = 1.5 to q = 2 . According to these observations, it is clear that not only Ca TP affects the front and back curvature radii of the plug but also the flow rate. The same conclusion is made by Kovalev et al. (2018). The plot of R tail -R front scaled by the channel width is also given in Fig. 16. It shows that R tail is always greater than R front and there is no significant change in the R tail -R front plot.
This indicates that the change in q affects nearly the same way the front and the back of the plug.
For, = 0.0086 , q = 1 (Q c = Q d = 0.2ml∕min) , the plug transverse interface is shown for different locations along the microchannel in Fig. 17.
As we can see, there is not much difference in the plug interface along the microchannel. It implies that the plug shape stabilises very quickly downstream of the junction. With a smaller viscosity ratio ( < 0.0086), Kovalev et al. (Kovalev et al. 2021) observed deformation of the plug interface downstream of the junction.

Plug Velocity and Velocity Profile
In this section, the plug/droplet velocity and velocity profiles inside the plug and in the liquid film region are investigated. We fixed the viscosity ratio (λ = 0.0086) and the flow rate of both phases to Q c = Q d = 0.2ml∕min;(q = 1).
The velocity vector field in the plane of symmetry (z = 0) is shown in Fig. 18.
As can be seen, the velocity vectors are oriented from the rear interface of the plug to the centre, and from the centre to the front interface of the plug. It should be noted that the velocity vectors at the plug centre are more important compared to the plug side interfaces. Furthermore, in the region of the liquid film, the low velocity vectors indicate a low flow rate through this region compared to the total flow rate.
In the case where the continuous phase is more viscous than the dispersed phase (which is the case in our study), the continuity of the viscous stresses in the region of the liquid film drags the liquid of the dispersed phase (water in our case) which is near the interface. The same topology was observed in similar conditions by Sarrazin et al. (2006). At seven locations distanced by 150 m between them inside the plug, the longitudinal velocity profiles dimensionalized by U max , are presented in Fig. 19. The maximum velocity U max is located close to the plug centre, specifically at the location ( X C − Δx).
As we can see, most of the profiles show a parabolic shape, with a peak in the middle of the channel ( y∕W = 0 ). The maximum velocity is always located in the central region of the plug. As it moves towards the side interfaces of the plug, the velocity gradually decreases until it is no longer present at the wall. Besides, in the area near the front and tail cap interface, the velocities are relatively low with a flattened velocity profile. Figure 20 shows the longitudinal velocity profile normalised by the maximum velocity at the middle of the plug for different positions along the channel.
The longitudinal velocity profiles are almost superposed for all positions. This indicates that the plug develops rapidly downstream of the junction, which makes the longitudinal velocity profile independent of the longitudinal position. The same result was obtained by Sarrazin et al. (2006).
As can be seen in Fig. 21, even if the inlet flow for both phases is the same, the velocity values in the liquid film are lower than the velocities inside the plug. This confirms what was mentioned earlier about the low flow rate in this region.
Under the no-slip condition, the velocity is zero at the wall and then increases progressively as one approaches the interface between the two phases. To predict the plug/droplet's velocity, the following expression is used.
Where A D , A f , A , U f and U D are respectively the droplet area, the film area, channel cross section area, liquid film velocity and droplet velocity. In our study, the plugs/droplets are surrounded by a liquid film with a low flow rate compared to the total flow rate. A flow rate of less than 5% of the total flow rate is observed in the liquid film area for λ = 0.0086, while for λ = 0.0127, a flow rate of less than 4% was observed.
This confirms the literature results where it is reported that the flow rate in the liquid film region is less than 5% of the total flow rate (Morris 1995;. So, the previous expression becomes: This linear relationship between U D and U TP was observed in the literature (Gupta et al. 2013;Taha and Cui 2004). However, a power function was employed by Kovalev et al. (Kovalev et al. 2018) which they claim gives better results compared to a linear function.
In the present study, for λ = 0.0086 and 0.08 ≤ Ca ≤ 0.16 , the constant C was estimated as 1.65 , while for λ = 0.0127 and 0.053 ≤ Ca ≤ 0.107 the constant is 1.44. Figure 22 shows the predicted and measured plug/droplet velocity as a function of the mixing velocity for two viscosity ratios. Two correlations from the literature are used in this section (Table 5).
From Fig. 22 it can be seen that U D is always higher than the U TP . This is due to the presence of the liquid film around the plug/droplet. As the film increases, the gap between U D and U TP also increases (Li and Angeli 2017). Finally, for both λ, the present calculus results approved the linear tendency.
Moreover, U D obtained from simulations agrees well with the U D predicted by both expressions. The average relative error for both cases (λ = 0.0086; λ = 0.0127) is less than 5%.
Compared with Kovalev's correlation (Kovalev et al. 2018), the correlation correctly predicts U D for both viscosity ratios. In general, there is an average error of 8% and 7% respectively for λ = 0.0086 and λ = 0.0127. It is worth noting that Kovalev's correlation was originally proposed for liquid-liquid systems with low viscosity ratios (λ = 0.001) in a rectangular microchannel. When compared to Liu et al. 's correlation (Liu et al. 2005), the U D obtained from simulations follows the same trend as predicted by the correlation for both λ. Our results are slightly different from those predicted with an average error of 17% for λ = 0.0127 while for λ = 0.0086, an average error of 7% is noticed.

Effect of Viscosity Ratio λ
In this section, with a flow rate of 0.2ml∕min for both phases, the effect of the viscosity is analysed by considering four viscosity ratios λ. In this sense, we varied the viscosity of the continuous phase and kept the dispersed phase one fixed. However, for = 0.0086 the viscosity of both phases was changed using the same combination of fluids reported in Wang et al. (2021) (see Table 2). Figure 23 shows the evolution of the dimensionless liquid film thickness versus the viscosity ratio λ. As we can see, the dimensionless liquid film ( Corner ; y ) decreases non-linearly as the λ increases.

Film Thickness
The thickness of the liquid film results from the competition between the interfacial tension and the viscous force of the continuous phase (Ca number), as the Re number of the dispersed phase is practically constant throughout this part. Therefore, a decrease in λ ( C increases) leads to an increase in the viscous stress of the continuous phase ( Ca c increases). This means that the film thickness is governed by the viscous force of the continuous phase.
As λ increases ( C decreases), the viscosity of the dispersed phase becomes comparable to the continuous phase one. This induces a strong shear force at the interface between the two phases, resulting in a decrease in the liquid film thickness (Li and Angeli 2017;Tsaoulidis and Angeli 2016).
For the low viscosity ratio ( = 0.0086; = 0.0127) considered in this section, a non-uniform liquid film region is   (Fig. 24). A non-negligible leakage flow of the continuous phase is observed in this region as mentioned early by Liu et al. (2005).
The liquid film region becomes uniform and thinner between the sidewalls and the interface separating the two phases by increasing .

Plug Shape
The effect of the viscosity ratio on the shape of the front and tail interface as well as on the interface along the lateral sides is discussed below.
Globally, as we can see in Fig. 25, an increase in the λ ( c decreases) leads to an increase in the R front and R tail interface of the plug.
For λ = 0.0086, R tail and R front are strongly influenced by the viscous stress of the continuous phase. This effect gradually decreases as λ increases.
Indeed, at low Ca c ( Ca c < 0.007 ), the plug/droplet shape becomes independent since R tail and R front are not affected by the operating parameters ( U c , c ), and only the plug length changes in this range (Gupta et al. 2013;Lac and Sherwood 2009). However, for Ca c ( 0.006 < Ca c < 0.054 ) as is the case here, R tail and R front are strongly influenced by the change in λ.

Plug Velocity and Velocity Profile inside Plug
The velocity profile inside the plug and in the film region as a result of the λ is shown in Fig. 26. The profiles were taken in a section at the middle of the plug, normalized by the maximum velocity derived from the velocity profile of λ = 0.0086.
From Fig. 26, we can see that all profiles show a parabolic shape independently of the viscosity ratio λ. Moreover, a peak is located in the midline of the channel (y/W = 0), and this is for all λ.
By increasing λ, the maximum velocity value decreases but remains at the centre of the plug. Furthermore, a flattened velocity profile is observed moving from λ = 0.0086 to λ = 0.0656. This is mainly due to the decrease in the liquid film thickness.
The velocity profile in the liquid film region, on the two sides of the plug, is shown in Fig. 27 for different viscosity ratios λ.
As can be seen, the velocity profiles in the region of the liquid film are almost similar on both sides of the plug's lateral interface. However, for λ = 0.0127, The flow appears to be disturbed in the targeted section, resulting in a slightly different velocity profile on both sides in the liquid film region.
It should be noted that the velocity in the liquid film region is generally low compared to the velocity inside the plug.
For λ = 0.0086, due to the leakage flow surrounding the plug, the velocity in the film is large compared to all other viscosity ratios λ. By increasing λ ( y decreases), the velocity in the film decreases remarkably with a value that does not overcome 4% of the maximum velocity. Figure 28 shows the evolution of the plug velocity as a function of λ. As we can see, for all values of λ, U D remains higher than the mixture velocity U TP due to the presence of the liquid film around the plug (Li and Angeli 2017).
However, by increasing λ, U D decreases rapidly and approaches U TP . The decrease in U D is a result of the decrease of the liquid film in the lateral sides of the channel y .

Conclusion
A 3D simulation of the hydrodynamics of liquid-liquid twophase flow and flow patterns in a rectangular microchannel was carried out with the commercial software ANSYS Fluent using the volume of fluid (VOF) method to capture the interface. On the other hand, a dynamic mesh adaptation based on the volume fraction is used to increase the accuracy and make the interface as sharp as possible.
Different combinations of fluids were used to investigate the effect of parameters such as Ca TP , flow ratio q and viscosity ratio λ on the liquid film ( Corner ; y ), shape and velocity of plugs/droplets.
The main obtained results are summarised below: 1. The present numerical results are consistent and agree well with the experimental results of Ma et al. (2021). It was shown that for q ≤ 1 , there was a minor effect of the contact angle on the plug/droplet length. However, for q > 1 , the contact angle starts to affect the plug/droplet length progressively. 2. Varying the flow rates of both phases (silicone oil-water), we observed six flow patterns. These include jet flow, droplet flow, plug flow, throat annular flow, annular flow, and stratified flow. The annular flow was observed far from the junction whereas for the stratified flow, the flow does not attach to the channel wall and moves towards the centre of the channel due to wettability effects. 3. By increasing the Ca TP through the U TP , the leakage rate around the plugs increases. For U TP > 0.0185m∕s , the interface shrinkage in the lateral sides disappears, indicating the dominance of viscous and inertial effects on the formation mechanism. Otherwise, the droplet velocity U D shows a linear relationship with U TP . Furthermore, in all conditions, U D remains higher than U TP due to the presence of the liquid film around the plugs/droplets. 4. The liquid film thickness ( Corner ; y ) was compared to correlations from the literature. There is a lack of correlations that predict the liquid film for liquid-liquid systems in rectangular microchannels with a low viscosity ratio λ. For the liquid film evolution versus flow rate ratio q , we found that Corner and y decrease as q increases. But for the plug shape, R front and R tail increase slightly with increasing q . Moreover, same as Li and Angeli (2017), Wang et al. , and contrarily to Ma et al. (2021), Kovalev et al. (2018), R tail is found always higher than R front . The longitudinal velocity profiles inside the plug show parabolic shape with a maximum in their centre. Furthermore, in the region of the liquid film, the velocity becomes lower compared to the velocity inside the plug. 5. Regarding the effect of λ, the liquid film thickness ( Corner ; y ) decreases while increasing the viscosity ratio λ. Therefore, for the plug shape, R front and R tail increase with increasing λ. However, by further increasing λ, their effect on R front and R tail decreases remarkably. The longitudinal velocity profiles inside the plug all exhibit a parabolic pattern with a maximum in the middle of the plug whatever λ..
Finally, this study validates the use of an adaption technique for two-phase liquid-liquid flow systems in microchannels that exceed 500 m . Physically, it could help to better understand the viscous effects and other parameters on liquid-liquid two-phase flow systems and thus increase the control within the microreactors.

Authors' Contributions
The author Mohammed SAID developed the numerical model, conducted the simulations, wrote the manuscript, and prepared all the figures. The Co-author N. NAIT BOUDA supervised