A Review on the Hydrodynamics of Taylor Flow in Microchannels: Experimental and Computational Studies

Taylor flow is a strategy-aimed flow to transfer conventional single-phase into a more efficient two-phase flow resulting in an enhanced momentum/heat/mass transfer rate, as well as a multitude of other advantages. To date, Taylor flow has focused on the processes involving gas–liquid and liquid–liquid two-phase systems in microchannels over a wide range of applications in biomedical, pharmaceutical, industrial, and commercial sectors. Appropriately micro-structured design is, therefore, a key consideration for equipment dealing with transport phenomena. This review paper highlights the hydrodynamic aspects of gas–liquid and liquid–liquid two-phase flows in microchannels. It covers state-of-the-art experimental and numerical methods in the literature for analyzing and simulating slug flows in circular and non-circular microchannels. The review’s main objective is to identify the considerable opportunity for further development of microflows and provide suggestions for researchers in the field. Available correlations proposed for the transition of flow patterns are presented. A review of the literature of flow regime, slug length, and pressure drop is also carried out.


Introduction
Multiphase flow in micro-sized structures refers to a microflow in which two or more distinct phases are recognizable, i.e., a carrying or continuous phase and one or more dispersed phases. A two-phase flow denotes a combination of two distinct phases, including gas, liquid, and solid particles. Gas-liquid (GL) and liquid-liquid (LL) are two common types of multiphase flows in microchannels encountered in various practical applications, such as biomedical, pharmacological, engineering, and commercial purposes. Immiscible LL two-phase flows can also be observed in many industrial applications, where the dispersive liquid flow is introduced as droplets into the carrying liquid flow. Intensified liquid-liquid extraction as a separation process is another application for twophase flow in oil and gas industries. The gas bubbles with equivalent diameters of the microchannel diameter in a GL Taylor flow are surrounded by a thin liquid film of the continuous flow and separated by the liquid slugs [1]. The interaction between the different phases, flow rates, thermophysical properties, and geometrical details of the channel characterizes the flow pattern, liquid film thickness, gas bubble or droplet shape, pressure drop, heat, and mass transfer rates. Knowledge of two-phase flow characteristics enables researchers and designers to optimize channel size and operating conditions. Minimizing the pressure drop can reduce corrosion and erosion while providing a high heat transfer rate [1,2]. An essential feature of Taylor flow is the flow pattern in the liquid slug, which can also be classified in recirculation flow and bypass flow [2,3]. Therefore, knowledge of hydrodynamic properties of Taylor flow is certainly required to understand transport phenomena characteristics for improving the performance of microchannels and enabling operational conditions to be optimized.
A review on the hydrodynamic characteristics of Taylor flow was conducted by Angeli and Gavriilidis [1] for circular and non-circular small channels. Correlations for film thickness measurements were summarized and the key effects of Capillary number (Ca) on film thickness were discussed in detail. Recently, Etminan et al. [4] presented all of the correlations proposed for film thickness measuring in the literature. They found the best-fitted correlation between some experimental data and a wide range of Ca. It can be concluded that there is still a lack of certainties for the gradient of surface tension at the interfacial region resulting in discrepancies between experimental predictions and numerical/theoretical findings. Moreover, more attempts should be devoted to finding specific flow and operational conditions for enhancing the performance of microchannels.
For multiphase flows, it is necessary to predict the transport phenomena in terms of flow parameters, i.e., superficial flow velocity components, slug length, liquid film thickness, and the flow rates of each phase. Multiphase flow is valid as long as the separation between phases is recognizable at a scale greater than the molecular level. The ability to describe the behavior of multiphase flow requires a set of basic definitions and assumptions. The following subsections have been designed to represent the importance of these basic definitions and the role of dimensionless groups in describing the two-phase flows' behavior. This paper is aimed to review two-phase flows in microchannels, with a particular focus on the hydrodynamic aspects. Fundamentals of multiphase flows, such as basic definitions and dimensionless parameters are studied. Identification of flow pattern and bubble/droplet formation are discussed in experimental investigations, numerical simulations, and flow regime transitions. Correlations of slug lengths and pressure drop are presented and classified in terms of cross-sectional areas of microchannels.

Basic Definitions in Two-Phase Flows
The capillary length has been used to identify the compaction in micro-sized applications, such as microchannels in heat exchangers. A micro-sized channel can be adequately assumed when the hydraulic diameter of the channel is less than the capillary length. In contrast, for large-diameter channels, the Laplace number (La) can be properly considered as a suitable length scale for calculations instead of the bubble or droplet diameter. Many other definitions have been developed to describe the characteristics of multiphase flows, including fluid flow, and heat and mass transfers, which are presented in Table 1. It should be noted that due to the presence of more than one phase and the interaction between all phases involved in multiphase flow, the thermophysical properties have been expressed differently compared to those of a single-phase flow [3][4][5]. The ratio between interfacial and gravitational (buoyancy) effects Slip ratio S U g /U l The ratio of average real velocity of the gas and liquid phases Average velocity of gas phase U g The ratio of volumetric flow rate of the gas phase to tube cross-sectional area occupied by the gas phase flow Table 1. Cont.

Name Symbol Definition Description
Average velocity of liquid phase U l The ratio of volumetric flow rate of the liquid phase to tube cross-sectional area occupied by the liquid phase flow The ratio of the tube cross-sectional area (or volume) occupied by the gas phase to the tube cross-sectional area (or volume) Volumetric quality (dynamic holdup) β Q g /Q t The ratio of the volumetric flow rate of the gas phase to the total volumetric flow rate Two-phase friction multiplier ϕ 1+ C X + 1 X 2 0.5 A function of the Lockhart-Martinelli parameter (X) and the Chisholm constant (C)

Dimensionless Parameters
A dimensionless number can represent the ratio of two different forces or physical quantities, which play a significant role in the flow pattern and interaction between phases (see Table 2). The length scale of a channel, which appears for the majority of dimensionless numbers, is often defined as the diameter of tubes or the hydraulic diameter of ducts. Table 2. Some of the most common non-dimensional numbers in multiphase flow. The ratio of the capillary and the gravitational (buoyancy) effects

Name
Ohnesorge Oh The ratio of the viscous force to the inertia and the surface tension forces Reynolds Re ρUd/µ The ratio between the inertia and the viscous forces Suratman Su The ratio of the surface tension to the viscous forces The ratio of the inertial forces to the interfacial forces The relation of gravitational to viscous forces is Archimedes number (Ar) that the formal definition is given in Table 2. This number frequently appears in tube-shaped chemical process reactor designs, which describe the motion of fluids caused by the Hua et al. [7] found that ellipsoidal gas bubbles were produced as the Ar and Bond (Bo) numbers were gradually increased. In solid particle multiphase flows, the Ar has also been related to the drag force coefficient and the Re [8][9][10] and the wind threshold velocities of particles in GL flow [11]. The Ar number was also correlated with the inverse viscosity number (N) by Quan [12] in an upward or downward co-current Taylor flow, N = Ar 0.5 =   ρ l ρ l − ρ g gd 3 µ l 2   0.5 (2) For a large N, the elongated tail of the bubbles oscillated for upward flow, and the shortened tail was only observed for downward flow.
The Bond number (Bo), also called Eötvös number (Eo), is the inverse of the Ar number and denotes the ratio between the gravitational and capillary effects. Since the Bond number varies with the square of length-scale, each change, even small, in channel diameter causes significant variation in the Bo. For air-water two-phase flow in channels of a diameter less than 1 mm, which are very common in micro-sized applications, the Bo number is approximately in the order of 0.1 and the importance of viscous effects is predominant. Therefore, the gravitational effects can be negligible in most GL flows in microchannels. Hua et al. [7] and Uno and Kintner [13] showed that the gas bubble velocity was significantly impacted by the wall effects for high Bond numbers. They showed that bubble shape was significantly affected by the Bo number: elongated bubbles with spherical caps at both ends appeared for Bo > 50. Some other specific regions of Bo and BoRe 0.5 l numbers have been proposed by investigators to categorize the behavior of multiphase flows, for example, dominant surface tension for Bo ≤ 1.5, non-negligible inertia, surface tension, and viscous effects for 1.5 < Bo ≤ 11, and negligible surface tension effects for Bo > 11 [14][15][16][17]. The absorption of CO 2 in micro-scale reactors where the Chisholm parameter was correlated with the Bo number was studied by Ganapathy et al. [18]. Prajapati and Bhandari [19] quantified the instability increase of boiling flow in a microchannel for the region of Bo < 1.
The Cahn number (Cn) is defined as the ratio of the interface width (ξ) and the tube diameter or any other length scale (d), Cahn and Hilliard [20,21] and Cahn [22] focused their studies on the free energy of the volume of an isotropic system, where the densities of components involved are nonuniform. The concept of the Cahn number was introduced by them when they found an increase in the thickness with temperature. Soon after, the efforts to obtain the governing equation in a non-equilibrium situation led to a well-known Cahn-Hilliard equation involving the Cahn number for the first time [23][24][25][26]. For example, Choi and Anderson [27] coupled the Cahn-Hilliard theory with the Extended Finite Volume Method (XFVM) for the dynamic modeling of suspended particles in two-phase flows.
The importance of viscous and capillary effects in two-phase flow in microchannels has been correlated by the capillary number based on the liquid slug velocity, The capillary number is a dimensionless group to explain how viscous and surface tension forces affect the interface between the gas and liquid phases, and also between two immiscible liquids. Although, the length-scale does not appear in the capillary number, surface tension forces become more significant relative to gravity as the cross-sectional area of the channel is decreased. For most microflow applications, the approximate values of superficial velocity range from 10 µms −1 to 1 ms −1 , when the viscosity is 10 −3 kg m −1 s −1 , and the surface tension can be reasonably considered at 0.05 N m −1 . These flow conditions cause the Ca to range from 2 × 10 −7 to 2 × 10 −2 . For high viscosity liquids, such as silicone oil at a mixture velocity of 1 ms −1 , the Ca becomes~1. The Ca frequently appears in correlations measuring liquid film thickness , pressure drop [30,35,44,48,[53][54][55][56][57], and slug length [58][59][60][61] in multiphase flows. Recently, the poor accuracy of classical pressure drop correlations with an increase of the Ca was emphasized by Ni et al. [51]. This was due to the presence of waves around the tail of the bubble, change of the semi-spherical head of the bubble, and powerful circulation inside the liquid slugs. Patel et al. [52] indicated that the film thickness at the corner of a square microchannel decreased with an increase of Ca, where a linear expression was derived for the film thickness in terms of Ca. Several other researchers showed an increasing trend in film thickness with Ca. These hydrodynamic features will be discussed in the following sections.
The Re expresses the ratio of inertial to viscous effects and predicts the flow pattern in different situations: where ρ, µ, and U are the density, dynamic viscosity of the continuous phase fluid, and the velocity of flow, respectively. The thermophysical properties of a multiphase flow are different from a single-phase flow, which has been discussed by Awad [5] in detail.
Regarding the value of the inertia forces, the flow regime remained laminar for the Re numbers less than~2300 [62,63]. Since the cross-sectional area of a microchannel often has a diameter of 1 mm or smaller and the liquid velocity is less than 1 m s −1 , the Re for water as the continuous phase is~1000 or less. This means that the flow regime in the microchannel applications is normally laminar and the predominant effects of the viscous forces must be taken into account. Because more than one phase is involved, different approaches have been used by researchers to determine the value of Re and the thermophysical properties in the multiphase flows. To realize the importance of the Re definition, consider two following cases: first, the inertia effects of the gas phase are insignificant, when the gas and the liquid superficial velocities are in the same order. For instance, the air to water density ratio in two-phase flow is around 1.225 × 10 −3 , which indicates the dominant role of the liquid phase density in defining the Re. This situation can be recognized in a bubbly flow pattern. Second, the inertia effects of both phases are important when the gas superficial velocity is relatively high compared with the liquid superficial velocity. The friction factor is necessarily calculated for the pressure drop calculation, while the Re appears in other correlations of hydrodynamic concepts, see Section 6 for further information. The Re is in the correlations for liquid film thickness and pressure drop obtained by Heiszwolf et al. [39] and Han and Shikazono [46]. The product of the Re and the We numbers was introduced as a transition criterion from slug to slug-bubbly regimes by Suo and Griffith [64]. Jayawardena et al. [65] proposed the ratio of liquid phase to gas phase the Re numbers to identify the transitions from bubbly to slug and from slug to annular patterns. Several other non-dimensional numbers in multiphase flows have been gathered and discussed by Awad [5].
As the diameter of the channel and the velocity of flow decrease, the inertial, and gravitational effects become negligible, while the surface tension becomes a predominant factor. The knowledge of the hydrodynamics of multiphase flow is strictly intertwined with the dimensionless numbers presented in Table 3. According to this table, an accurate description of the flow map, bubble or droplet formation, and transition from one flow regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, and Fr. While, Bo, Ca, and Re have appeared in the pressure drop correlations so far, the friction factor is only governed by the Re. The liquid film thickness has been correlated or described with Re, Ca, and We numbers. Finally, the slug length has only included the Re and the Ca numbers. regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, and Fr. While, Bo, Ca and Re have appeared in the pressure drop correlations so far, the friction factor is only governed by the Re. The liquid film thickness has been correlated or described with Re Ca, and We numbers. Finally, the slug length has only included the Re and the Ca num bers. regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, an and Re have appeared in the pressure drop correlations so far, the fr governed by the Re. The liquid film thickness has been correlated or Ca, and We numbers. Finally, the slug length has only included the R bers.  with the dimensionless numbers presented in Table 3. According to this table, an accurate description of the flow map, bubble or droplet formation, and transition from one flow regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, and Fr. While, Bo, Ca, and Re have appeared in the pressure drop correlations so far, the friction factor is only governed by the Re. The liquid film thickness has been correlated or described with Re, Ca, and We numbers. Finally, the slug length has only included the Re and the Ca numbers. with the dimensionless numbers presented in Table 3. According to this table, an accurate description of the flow map, bubble or droplet formation, and transition from one flow regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, and Fr. While, Bo, Ca and Re have appeared in the pressure drop correlations so far, the friction factor is only governed by the Re. The liquid film thickness has been correlated or described with Re Ca, and We numbers. Finally, the slug length has only included the Re and the Ca num bers. with the dimensionless numbers presented in Table 3. According to th description of the flow map, bubble or droplet formation, and transi regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, an and Re have appeared in the pressure drop correlations so far, the fr governed by the Re. The liquid film thickness has been correlated or Ca, and We numbers. Finally, the slug length has only included the R bers. factor. The knowledge of the hydrodynamics of multiphase flow is strictly intertwined with the dimensionless numbers presented in Table 3. According to this table, an accurate description of the flow map, bubble or droplet formation, and transition from one flow regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, and Fr. While, Bo, Ca and Re have appeared in the pressure drop correlations so far, the friction factor is only governed by the Re. The liquid film thickness has been correlated or described with Re Ca, and We numbers. Finally, the slug length has only included the Re and the Ca num bers. factor. The knowledge of the hydrodynamics of multiphase flow is with the dimensionless numbers presented in Table 3. According to th description of the flow map, bubble or droplet formation, and transi regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, an and Re have appeared in the pressure drop correlations so far, the fr governed by the Re. The liquid film thickness has been correlated or Ca, and We numbers. Finally, the slug length has only included the R bers. gravitational effects become negligible, while the surface tension beco factor. The knowledge of the hydrodynamics of multiphase flow is with the dimensionless numbers presented in Table 3. According to th description of the flow map, bubble or droplet formation, and transi regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, an and Re have appeared in the pressure drop correlations so far, the fr governed by the Re. The liquid film thickness has been correlated or Ca, and We numbers. Finally, the slug length has only included the R bers. been gathered and discussed by Awad [5].
As the diameter of the channel and the velocity of flow decrease, the inertial, and gravitational effects become negligible, while the surface tension becomes a predominant factor. The knowledge of the hydrodynamics of multiphase flow is strictly intertwined with the dimensionless numbers presented in Table 3. According to this table, an accurate description of the flow map, bubble or droplet formation, and transition from one flow regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, and Fr. While, Bo, Ca, and Re have appeared in the pressure drop correlations so far, the friction factor is only governed by the Re. The liquid film thickness has been correlated or described with Re, Ca, and We numbers. Finally, the slug length has only included the Re and the Ca numbers. been gathered and discussed by Awad [5].
As the diameter of the channel and the velocity of flow decrease, the inertial, and gravitational effects become negligible, while the surface tension becomes a predominan factor. The knowledge of the hydrodynamics of multiphase flow is strictly intertwined with the dimensionless numbers presented in Table 3. According to this table, an accurate description of the flow map, bubble or droplet formation, and transition from one flow regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, and Fr. While, Bo, Ca and Re have appeared in the pressure drop correlations so far, the friction factor is only governed by the Re. The liquid film thickness has been correlated or described with Re Ca, and We numbers. Finally, the slug length has only included the Re and the Ca num bers. been gathered and discussed by Awad [5].
As the diameter of the channel and the velocity of flow decrease, the ine gravitational effects become negligible, while the surface tension becomes a pre factor. The knowledge of the hydrodynamics of multiphase flow is strictly in with the dimensionless numbers presented in Table 3. According to this table, a description of the flow map, bubble or droplet formation, and transition from regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, and Fr. Whi and Re have appeared in the pressure drop correlations so far, the friction fac governed by the Re. The liquid film thickness has been correlated or described Ca, and We numbers. Finally, the slug length has only included the Re and the bers. been gathered and discussed by Awad [5].
As the diameter of the channel and the velocity of flow decrea gravitational effects become negligible, while the surface tension beco factor. The knowledge of the hydrodynamics of multiphase flow is with the dimensionless numbers presented in Table 3. According to th description of the flow map, bubble or droplet formation, and transi regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, an and Re have appeared in the pressure drop correlations so far, the fr governed by the Re. The liquid film thickness has been correlated or Ca, and We numbers. Finally, the slug length has only included the R bers. been gathered and discussed by Awad [5].
As the diameter of the channel and the velocity of flo gravitational effects become negligible, while the surface ten factor. The knowledge of the hydrodynamics of multiphas with the dimensionless numbers presented in Table 3. Accor description of the flow map, bubble or droplet formation, a regime to another are highly depended on Bo or Eo, Ca, Re, W and Re have appeared in the pressure drop correlations so governed by the Re. The liquid film thickness has been corr Ca, and We numbers. Finally, the slug length has only inclu bers. been gathered and discussed by Awad [5].
As the diameter of the channel and the veloc gravitational effects become negligible, while the su factor. The knowledge of the hydrodynamics of m with the dimensionless numbers presented in Table  description of the flow map, bubble or droplet for regime to another are highly depended on Bo or Eo and Re have appeared in the pressure drop correl governed by the Re. The liquid film thickness has Ca, and We numbers. Finally, the slug length has o bers. to annular patterns. Several other non-dimensional numbers in multiphase flows have been gathered and discussed by Awad [5].
As the diameter of the channel and the velocity of flow decrease, the inertial, and gravitational effects become negligible, while the surface tension becomes a predominant factor. The knowledge of the hydrodynamics of multiphase flow is strictly intertwined with the dimensionless numbers presented in Table 3. According to this table, an accurate description of the flow map, bubble or droplet formation, and transition from one flow regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, and Fr. While, Bo, Ca, and Re have appeared in the pressure drop correlations so far, the friction factor is only governed by the Re. The liquid film thickness has been correlated or described with Re, Ca, and We numbers. Finally, the slug length has only included the Re and the Ca numbers. to annular patterns. Several other non-dimensional numbers in multiphase flows have been gathered and discussed by Awad [5].
As the diameter of the channel and the velocity of flow decrease, the inertial, and gravitational effects become negligible, while the surface tension becomes a predominan factor. The knowledge of the hydrodynamics of multiphase flow is strictly intertwined with the dimensionless numbers presented in Table 3. According to this table, an accurate description of the flow map, bubble or droplet formation, and transition from one flow regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, and Fr. While, Bo, Ca and Re have appeared in the pressure drop correlations so far, the friction factor is only governed by the Re. The liquid film thickness has been correlated or described with Re Ca, and We numbers. Finally, the slug length has only included the Re and the Ca num bers. to annular patterns. Several other non-dimensional numbers in multiphase fl been gathered and discussed by Awad [5].
As the diameter of the channel and the velocity of flow decrease, the ine gravitational effects become negligible, while the surface tension becomes a pre factor. The knowledge of the hydrodynamics of multiphase flow is strictly in with the dimensionless numbers presented in Table 3. According to this table, a description of the flow map, bubble or droplet formation, and transition from regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, and Fr. Whi and Re have appeared in the pressure drop correlations so far, the friction fac governed by the Re. The liquid film thickness has been correlated or described Ca, and We numbers. Finally, the slug length has only included the Re and the bers. to annular patterns. Several other non-dimensional numbers in mul been gathered and discussed by Awad [5].
As the diameter of the channel and the velocity of flow decrea gravitational effects become negligible, while the surface tension beco factor. The knowledge of the hydrodynamics of multiphase flow is with the dimensionless numbers presented in Table 3. According to th description of the flow map, bubble or droplet formation, and transi regime to another are highly depended on Bo or Eo, Ca, Re, We, Su, an and Re have appeared in the pressure drop correlations so far, the fr governed by the Re. The liquid film thickness has been correlated or Ca, and We numbers. Finally, the slug length has only included the R bers. to annular patterns. Several other non-dimensional numbe been gathered and discussed by Awad [5].
As the diameter of the channel and the velocity of flo gravitational effects become negligible, while the surface ten factor. The knowledge of the hydrodynamics of multiphas with the dimensionless numbers presented in Table 3. Accor description of the flow map, bubble or droplet formation, a regime to another are highly depended on Bo or Eo, Ca, Re, W and Re have appeared in the pressure drop correlations so governed by the Re. The liquid film thickness has been corr Ca, and We numbers. Finally, the slug length has only inclu bers. to annular patterns. Several other non-dimension been gathered and discussed by Awad [5].
As the diameter of the channel and the veloc gravitational effects become negligible, while the su factor. The knowledge of the hydrodynamics of m with the dimensionless numbers presented in Table  description of the flow map, bubble or droplet for regime to another are highly depended on Bo or Eo and Re have appeared in the pressure drop correl governed by the Re. The liquid film thickness has Ca, and We numbers. Finally, the slug length has o bers. A graphical flow map, which displays the specific region for each flow regime along with their transition lines, uses a few important characteristics of multiphase flow, such as superficial velocity, Re, Ca, and We. Table 4 summarizes the non-dimensional numbers and other flow parameters that use the x-and y-axis coordinates for quantitative-based flow map diagrams. Table 4. Some non-dimensional numbers and thermophysical properties selected by authors for displaying flow maps.
As a guide to reading this table; the kinetic energy of flow (ρV 2 ) was used by the author's group of 7 which its row and column are colored red. The first effort to display a flow pattern diagram dates back to 1953 when Baker [84] computed the mass velocity of each phase and two other expressions involving the viscosity and densities (G/λ and Bλψ/G) of both phases (see '5' in Table 4). These dimensional terms were not chosen by others afterward. In contrast, it is the interfacial velocities of the phases that have frequently been selected by scholars for more than half a century (see '1' in Table 4). After velocity, it seems that Re, We, and Ca numbers have the highest use by researchers for showing the flow pattern clearly and more understandably.

Identification of Flow Patterns and Bubble/Slug Formation
The first step of transport phenomena analysis of multiphase flows is to explain the flow characteristics in a phase-mapped diagram. The transition lines in such maps specify the criterion of a transition from one flow regime to another. A two-phase flow usually consists of a train of dispersed bubbles or droplets carried by a continuous phase, which potentially enhances the rates of heat and mass transfers in industrial applications. Although the flow patterns in microchannels and large-diameter channels show similarities, there are several differences. As the channel diameter increases, the laminar flow regime starts to be unstable due to the inertial effects. Kawaji and Chung [71] and Akbar et al. [70,87] indicated flow patterns in micro and minichannels, where the stratified flow regime was not present for larger channels. The effects of increasing the Reynolds number on the bubble profile were numerically and experimentally discussed by Kreutzer et al. [88]. Their results indicated a more elongated nose and flattened tail due to the lack of large enough interfacial forces to keep the hemispherical shape of the bubble caps as Reynolds number increased. Figure 1 presents a schematic of GL flow patterns in capillaries. The geometrical details of the channel, configuration of the channel, properties of phases, superficial velocities of phases, and types of junctions are the predominant factors to determine flow patterns in microchannels. The wetting or drainage effects of each phase are responsible for at least part of the flow maps in multiphase flows. In a bubbly flow regime, the individual bubbles move through the continuous liquid phase at very low liquid superficial velocities. When the bubble size is small, the interaction between bubbles is negligible. Conversely, increasing the gas-to-liquid volumetric flow rate causes the bubbles to coalesce and a breakup occurs ( Figure 1). Bubble coalescence becomes stronger as the gas-to-liquid volumetric flow rate is increased, to make large long bubbles with a rounded front that spans the cross-sectional area of the channel. Characteristics of the bullet-shaped bubbles most often called Taylor or slug flow in which the Taylor bubbles are surrounded by a thin layer of the carrier liquid phase, Figure 1. Some small dispersed bubbles are between the Taylor bubbles as the gas-to-liquid volumetric flow rate is enhanced to make a transition from slug to churn patterns (Figure 1e). In a churn flow pattern, neither dispersed nor carrier phases are continuous and irregular plugs of gas flow through the liquid phase and a wide variety of bubble lengths can be observed. This chaotic flow pattern occurs when the velocity is increased (Figure 1f). In large-diameter tubes, an oscillating chaotic flow is much more remarkable than that in small-diameter tubes. This flow regime is also named a semi-annular or unstable-slug flow pattern. In an annular flow pattern, a very thin thickness of the liquid phase remains on the walls of the channel for two flow conditions: the gas-to-liquid volumetric flow rate increases and the velocities of phases rise at low liquid volume fraction (Figure 1g). While very small-diameter droplets of the liquid phase are in the core of the gas phase, the flow pattern is mist (Figure 1h).
An experimental study of a hydrophobic ionic/water two-phase flow in two T and Y inlet junctions on the flow patterns was carried out by Tsaoulidis et al. [89], resulting in a negligible impact of the inlet configuration. They also investigated the effects of capillaries' materials made of glass, Fluorinated Ethylene Propylene (FEP), and Tefzel on flow patterns illustrated in Figure 2. Although the boundaries of the FEP and Tefzel are similar, the plug flow regime occupies a larger area in the flow map compared to the FEP capillary, and an annular pattern is observed in higher mixture velocities and lower ionic volume fraction. They did not observe annular and drop patterns in the glass capillary when the continuous phase was water. An experimental study of a hydrophobic ionic/water two-phase flow in two T and Y inlet junctions on the flow patterns was carried out by Tsaoulidis et al. [89], resulting in a negligible impact of the inlet configuration. They also investigated the effects of capillaries' materials made of glass, Fluorinated Ethylene Propylene (FEP), and Tefzel on flow patterns illustrated in Figure 2. Although the boundaries of the FEP and Tefzel are similar, the plug flow regime occupies a larger area in the flow map compared to the FEP capillary, and an annular pattern is observed in higher mixture velocities and lower ionic volume fraction. They did not observe annular and drop patterns in the glass capillary when the continuous phase was water.

Experimental Investigations
Transport phenomena for multiphase flows in capillary passages involve the behavior of each phase individually, the mutual interactions of phases, and the interactions between phases and solid boundaries. The presence of gas bubbles in a multiphase flow can be considered as an obstacle to the liquid phase flow, which is known as the Jamin effect [90]. This phenomenon is prominent enough to decrease the productions of the petroleum industry due to the high-pressure drop and the flow pattern transition. At the interface between two components of a GL two-phase flow, the surface tension gradient is substantial, which affects the interaction between phases and other transport phenomena. Historically, this effect is known as the Marangoni effect, which pushes the liquid phase to flow away from low surface tension regions and was first explained by Gibbs [91].
A comprehensive experimental study on the flow maps by means of high-speed Xray photography and flash methods simultaneously was reported by Hewitt and Roberts

Experimental Investigations
Transport phenomena for multiphase flows in capillary passages involve the behavior of each phase individually, the mutual interactions of phases, and the interactions between phases and solid boundaries. The presence of gas bubbles in a multiphase flow can be considered as an obstacle to the liquid phase flow, which is known as the Jamin effect [90]. This phenomenon is prominent enough to decrease the productions of the petroleum industry due to the high-pressure drop and the flow pattern transition. At the interface between two components of a GL two-phase flow, the surface tension gradient is substantial, Processes 2021, 9, 870 9 of 40 which affects the interaction between phases and other transport phenomena. Historically, this effect is known as the Marangoni effect, which pushes the liquid phase to flow away from low surface tension regions and was first explained by Gibbs [91].
A comprehensive experimental study on the flow maps by means of high-speed X-ray photography and flash methods simultaneously was reported by Hewitt and Roberts [66], which recognized the previous correlations for each flow pattern proposed by Baker [84] and Wallis [92]. A rough sketch of streamlines in the slug region of a GL flow was experimentally obtained by Taylor [31], where a stagnation point on the vortex and a stagnation ring on the curvature of the bubble were observed for two extreme cases of velocity fraction, (V g − V m )/V g , greater and less than 0.5. The high-speed photography was compared with X-ray technology to provide clear images of air-water flow regimes through a capillary [66]. Irandoust and Andersson [34] experimentally depicted Taylor's bubble profile, where the film thickness remained constant in the middle of the gas bubble. They also confirmed that the gas bubble was elongated in the front and squeezed in the rear parts. A review study on the flow maps was conducted by Akbar et al. [70] to show the dominant effects in each flow regime using available experimental data. They proposed a hydraulic diameter of 1 mm as a criterion for classifying the previous studies into the five regions: bubbly, plug or slug (surface tension dominated), annular (inertia dominated), dispersed, and transition zones regarding the amounts of superficial velocities and the Weber numbers. Micro-particle image velocimetry (µ-PIV) and fluorescence microscopy techniques were used by Günther et al. [73] to recognize a micro-segmented GL flow in a rectangular cross-section area channel. They showed that the roughness of the inner side of the channel and compressibility of the gas bubbles induce an imbalance into the flow field and accelerate the mixing rate accordingly. The periodic switching of recirculation regions in the liquid slugs caused a higher mixing rate in a meandering channel compared with a straight channel ( Figure 3). Dore et al. [93] investigated the dynamics of water/ionic two-phase flow by exploring slug formation and recirculation flow in fluid segments at T-junction, straight and curved microchannel employing µ-PIV. As shown in Figure 4, the circulation patterns in water plug consisted of two main mirrored and counter rotating vortices for low mixture velocity. As the mixture velocity increased, two secondary vortices arose at the nose cap of water plug.   A combination of experimental observations and numerical simulations was employed by Meyer et al. [94] to show velocity profile, flow patterns in a 2 mm × 2 mm vertical channel, Figure 5. Both methods predicted a Couette velocity distribution within the liquid film region, a mirrored main recirculating region, and two small vortices in the nose and rear caps with a direction of opposite. They also realized some differences between the results of the two methods as shown previously by [95].
Experimental analysis by Zhao et al. [82] revealed that the volume of dispersed flow was significantly affected by interfacial forces, inertia force, and the volumetric ratio of dispersed phase, which correlated by the multivariable least squares method,  The flow patterns of an LL two-phase flow in a rectangular cross-section duct were experimentally investigated by Zhao et al. [82]. The flow maps and the flow transitions were correlated by the inertia forces of each phase and the interfacial forces. They realized six distinctive flow maps in terms of dispersed phase droplet formation process at a T-junction as follows: • Slug regime; when the interfacial tension is greater than inertial forces, and the Weber numbers are 7.61 × 10 −6 < We ws < 4.87 × 10 −2 and 5.94 × 10 −6 < We ks < 5.94 × 10 −4 . • Monodispersed droplet regime; when the carrier phase flow rate is increased, and the Weber numbers.

•
The droplet population regime; when 1.07 < We ws < 30.43 and 3.8 × 10 −4 < We ks < 2.38 × 10 −1 . This flow regime is observed at the center of the channel where the inertial effects of the continuous phase are significantly greater than the interfacial tension. • Chaotic thin striations regime; when both flow rates are increased the We numbers are in a range of 0.17 < We ws < 30.43 and 4.29 < We ks < 53.5. This flow map will be eventually changed into annular due to the instability of the flow.
A combination of experimental observations and numerical simulations was employed by Meyer et al. [94] to show velocity profile, flow patterns in a 2 mm × 2 mm vertical channel, Figure 5. Both methods predicted a Couette velocity distribution within the liquid film region, a mirrored main recirculating region, and two small vortices in the nose and rear caps with a direction of opposite. They also realized some differences between the results of the two methods as shown previously by [95].
The diameter of the channel has a major role in the measured void fraction and a flow pattern which is shown in Figure 6 for the small diameter tube of 38.1 mm (upper row) and large diameter tube of 101.6 mm (lower row). An increase in tube diameter resulted in more compressed gas bubbles toward the top wall and occupied a smaller region of the cross-sectional area [96]. The measurement also showed that a further increase in superficial velocity ratio makes a more uniform gas phase distribution realized by a smaller portion of gas volume in the top region of the tube. As a consequence, the bubbly to plug/slug transition occurs at lower superficial velocity as the channel's diameter increases (i.e., = 0.51 m/s for the smaller pipe and = 0.25 m/s for the larger pipe). A critical Reynolds number of ⁓300 was found by Butler et al. [97], where the timeaveraged gas phase volume fraction showed a significant difference by means of the PIV measurements in capillaries with a T-junction. As shown in Figure 7, a large mirrored recirculation region was observed between two consecutive O2 bubbles with a relative velocity of 2V tp -V g resulting in the most efficient convective transport (regimes 1, 2, 11, and 13 in Figure 7). They also revealed the key role of the recirculating motion in the liquid plug on the dynamics of mass transfer and its dependency on the bubble velocity.  Experimental analysis by Zhao et al. [82] revealed that the volume of dispersed flow was significantly affected by interfacial forces, inertia force, and the volumetric ratio of dispersed phase, which correlated by the multivariable least squares method, where the dynamic hold-up fraction (ε) and equivalent radius (R) in terms of the volume of dispersed phase (V d ) are, The diameter of the channel has a major role in the measured void fraction and a flow pattern which is shown in Figure 6 for the small diameter tube of 38.1 mm (upper row) and large diameter tube of 101.6 mm (lower row). An increase in tube diameter resulted in more compressed gas bubbles toward the top wall and occupied a smaller region of the cross-sectional area [96]. The measurement also showed that a further increase in superficial velocity ratio makes a more uniform gas phase distribution realized by a smaller portion of gas volume in the top region of the tube. As a consequence, the bubbly to plug/slug transition occurs at lower superficial velocity as the channel's diameter increases (i.e., j g = 0.51 m/s for the smaller pipe and j g = 0.25 m/s for the larger pipe).  A critical Reynolds number of~300 was found by Butler et al. [97], where the timeaveraged gas phase volume fraction showed a significant difference by means of the PIV measurements in capillaries with a T-junction. As shown in Figure 7, a large mirrored recirculation region was observed between two consecutive O 2 bubbles with a relative velocity of 2V tp -V g resulting in the most efficient convective transport (regimes 1, 2, 11, and 13 in Figure 7). They also revealed the key role of the recirculating motion in the liquid plug on the dynamics of mass transfer and its dependency on the bubble velocity.  Two different structures named high-concentrated and low-concentrated were observed by Yao et al. [98] to characterize the mass transfer and flow patterns of liquidliquid slug flow at the bend region of a rectangular meandering microchannel. As is shown in Figure 8a, smaller flow rates created a non-crossing evolution pattern where the red filaments occupied the outer and inner layers when moving through the bend. Conversely, when the flow rates were increased, Figure 8b, the outer red filament was shifted to the central region of the bend. Table 5 presents selected experimental studies on the flow patterns of two-phase flows in tubular and non-tubular microchannels and gas bubble or droplet formation. Two different structures named high-concentrated and low-concentrated were observed by Yao et al. [98] to characterize the mass transfer and flow patterns of liquid-liquid slug flow at the bend region of a rectangular meandering microchannel. As is shown in Figure 8a, smaller flow rates created a non-crossing evolution pattern where the red filaments occupied the outer and inner layers when moving through the bend. Conversely, when the flow rates were increased, Figure 8b Table 5 presents selected experimental studies on the flow patterns of two-phase flows in tubular and non-tubular microchannels and gas bubble or droplet formation. Table 5. The literature of selected experimental studies on the flow pattern maps and bubble or droplet formation.

Comment(s) Cross-Section Phases Reference
Several correlations as the function of velocity and pressure drop in each phase were suggested for different flow regimes Circular GL [84] Combination of X-ray and high speed flash photography technique High liquid flow rates Circular GL [66] Recognizing a novel flow map using Su number for microgravity two-phase flow Circular GL [65] Boiling heat transfer of R141b Heat transfer coefficient correlations Different flow regime in small-and large-diameter tubes Observation of local dry-out on the channel wall Circular GL (single phase) [99] Using particle image velocimetry Recirculation flow pattern with a high degree of mixing Counter-rotating vortices were observed inside the liquid slugs Velocity profile inside the slugs Circular Square GL [37] Increasing gas superficial velocity led to the development of the slug flow Low liquid superficial velocity made longer bubbles, shorter liquid slugs profile Slug and annular patterns At high liquid superficial velocity, churn flow was established Circular Semi-Triangular GL [68]

Comment(s) Cross-Section Phases Reference
Five flow regimes were observed including bubbly, wedging, slug, bubbly, and dry flow for moderate void fraction The classification of patterns regarding liquid fraction The effect of a sharp return of the channel Void fraction measurement The gravitational effects were taken into account for a bubbly pattern with a spherical gas bubble Liquid droplets may be observable on the walls in wedging pattern Square GL [72] Using µ-PIV and fluorescent microscopy imaging The liquid phase segments were attached at the corners of the cross-section Gas phase flow improved the mixing and the residence time features Rectangular GL [73] The hemispherical ends of gas bubbles were not maintained as the Ca increases The nose and the tail of the bubble elongated and flattened as the Re increased The Marangoni effect was observed in experiment efforts and was not taken into account in the numerical simulations Circular GL [44] Surface modifications Contact angle measurements Meandering microchannel When the length of the bubble was equal to the channel diameter, the flow map was annular At low void fraction when the surface energy was low, the isolated symmetric bubble patterns were observed At moderate void fractions, the flow pattern was asymmetric Square GL [100] Discussed in the text Rectangular LL [82] An increase in Ca changes the droplet profile Development of droplet formation with the increase in Ca The significant deviation between the measured film thickness and the Bretherton and Taylor predictions Circular LL [45] At relatively low flow rates, the slug pattern was established with the length equal to the inner diameter of the tube At high volumetric flow rate, the deformed interface regime was made with long water slugs and small cyclohexane droplets Circular LL [101,102] Using a µ-PIV technique to capture bubble formation in a T-junction About 25% of the liquid phase passed the gas bubble The gas bubble formation process was included in three steps: Gas bubble growing until occupied the junction Gas bubble developing was decreased as the liquid phase passes the bubbles and the bubble neck was tightened Gas bubble neck was rapidly decreased until approached one-quarter of the diameter before breaking up Square GL [103] An empirical correlation for a transition from Taylor to unstable slug flow regimes At lower superficial velocities, the Taylor regime was established At higher superficial velocities, the Taylor flow was transited into the dripping flow pattern Circular GL [74,76] Boiling heat transfer of Fluorinert FC-77 As the heat flux enhanced, the bubble generation rate increased along with an increase of bubble length At moderate heat flux, flow pattern went back and forth between churn and wispy-annular flow regimes Rectangular GL (single phase) [104]

Comment(s) Cross-Section Phases Reference
Bubbly flow regime was found at a low gas superficial velocity and a high liquid superficial velocity An increase in gas superficial velocity generated slug regime, where the gas-bubble length was longer than the diameter High liquid velocity and moderated gas superficial velocity produced churn flow map Low liquid superficial velocity and high gas superficial velocity transited the churn to the slug-annular Circular GL [105] PIV measurements were taken into account Backflow around the liquid slugs were observed The liquid slugs moved faster than that of average flow due to the lubricating effects Circular LL [106] Using a µ-PIV technique to capture flow structure Slug or Taylor pattern was observed for low superficial velocities As the superficial velocities increased, the lengths of bubbles, and slugs became more variable At high enough gas velocity, the gas phase penetrates the liquid and the length of slugs became shorter Circular GL [107] Using a laser signal system to provide flow pattern images The channel diameter affects the flow pattern and transition conditions In the bubbly flow regime, the size of the bubbles was approximately was the same as channel diameter Circular GL [108] A model of mass and momentum balance Stratified flow regime was observed for low oil to gas superficial velocity ratios An increase in the velocity of the gas phase makes the interface to be wavy Circular GL [109] The number of small bubbles in liquid plug/slug was significantly increased with increasing V g The size of small bubbles was decreased with increasing V g or increasing V l Increasing V g or decreasing V l slightly increased the depth of the plug/slug bubble Increasing pipe size enhanced the contribution from large bubbles to a total void fraction Large bubble was accelerated due to small bubble coalescence as flow developed, leading to α decreases Circular GL [96] Indicating of wave growth and coalescence during slug flow formation For low superficial velocities, the flow pattern remained stratified smooth By an increase of superficial gas velocity, the flow pattern was changed to wavy Circular GL [79]

Numerical Simulations
Numerical modeling methods can be categorized into four main groups; electronic, atomistic, mesoscale, and continuum [110]. In the following, we only review some of the widely-used numerical methods in the modeling of Taylor flows, and more in-depth reviews on the rest of numerical studies can be found in literature, i.e., [111].
Two different methods may be considered in numerical analysis when describing the motions of fluids in two-phase flows: continuous fluid (Eulerian) and discrete particle (Lagrangian). The first method solves the governing partial differential equations to predict the motion of the continuum-based fluid flow, while the second method follows the movement of fluid particles or molecules to predict fluid flow and heat transfer. The interaction between different phases is calculated by both methods to determine the coalescence and interference of interfacial forces [112]. Hashim et al. [113] simulated transport phenomena and biological components, such as proteins, DNA, and cells using the Eulerian method. The innovative aspect of their simulation was to describe the mixing concentration distribution at the interaction of biological components. Apte et al. [114] numerically investigated the liquid fuel spray from a nozzle using the Eulerian-Lagrangian approach with a pointparticle approximation for the liquid droplets. Ni et al. [51] investigated numerical study of the GL and LL Taylor flows on an Arbitrary Lagrangian-Eulerian (ALE) formulation to track interfacial and curvatures of slugs. They showed the predominant roles of the viscosity ratio and density ratio of continuous phase to the disperse phase for quite large Ca causing a non-stagnant liquid film region. The Lagrangian approach was also used by [115][116][117][118][119], but not limited to. Therefore, these methods can be widely employed in applications concerning fluid motion.
The Lattice Boltzmann Method (LBM) is a Computational Fluid Dynamics (CFD) method that simulates a fluid density on a lattice framework with streaming and collision (relaxation) processes instead of solving the Navier-Stokes equation directly [120]. This method is an efficient tool for modeling complex fluid systems that include complex boundaries, microscopic interactions, propagations, and the collusions of particles [121][122][123][124][125][126][127][128]. The LBM in a GL two-phase flow modeling was proposed by Seta and Kono [129], based on a particle velocity-dependent forcing scheme. Their results revealed that the effects of high potential caused instability in a three-dimensional model. A diffuse interface free energy LBM was utilized by Komrakova et al. [130] to study the behavior of a single liquid-droplet suspended in another liquid. They found that for droplets with a radius of less than 30 lattice units, a smaller interface thickness is required. Li et al. [131] comprehensively reviewed different uses of LBM proposed by others for simulating thermal boundary treatments, interaction forces between two different phases, mechanical stability conditions, liquidvapor phase changes, fuel cells, droplet collisions, and energy storage systems. Recently, Fei et al. [132] developed a three-dimensional, multiple-relaxation-time LBM based on a set of non-orthogonal vectors for modeling realistic multiphase flows. Some other uses of the LBM have also been developed by many investigators, e.g., [133][134][135].
Another widely used method for modeling two-phase flows is the Volume of Fluid (VoF), which is a numerical technique that determines the location of a free surface, GL, and LL interface by following the Eulerian approach. This numerical method is a powerful tool for modeling the interfaces of incompressible and immiscible two-phase flows because it is able to predict the interface between two phases even though the interface is becoming too weak to capture the curvature of the interfacial line [136]. Ketabdari [137] discussed VoF regarding the governing equations, as well as its advantages and capabilities predicting the interface line between two phases with large deformation. Osher and Sethian [138] devised a new numerical algorithm to realize and determine the propagation of curvaturedependent speed, which can be combined with VoF. This algorithm approximates the equations of motion by using hyperbolic conversation laws. Many investigators employed this method when simulating a two-phase flow's interface, e.g., [139][140][141][142]. The general applications of this method and a comprehensive review were conducted by Osher and Fedkiw [143], which discusses the variants and extension methods, such as fast methods for steady-state problems, diffusion generated motion and the variation level set approach. The growth of liquid film thickness was governed by one equation [144]. They assumed a constant pressure gradient in the channel direction and a linear velocity distribution over the liquid film thickness, emphasizing that the model was valid only for a film thickness less than near-wall grid size.
An effort to predict the bubble profile was analytically conducted by Bretherton [30], where he found that considering gravitation as a buoyancy force produced nonsensical results in the vertical configuration of a capillary tube. A critical value ρgR 2 /σ of 0.842 was introduced to show no upward movement of the bubble for the amount less than the critical value. Kolb and Cerro [145] found an axisymmetric bubble profile for the Ca greater than 0.1, where the lubrication's law predicts the interface profiles and the flow patterns adequately. The effects of bubble and slug lengths on the mass transfer in a GL flow was experimentally studied by Berčič and Pintar [146]. Their results revealed that the amount of nitrite transported was independent of the gas bubble length and primarily depended on the liquid slug length. The flow field predictions were numerically obtained by Giavedoni and Saita [36], who illustrate the streamlines for showing the bubble profiles and the slug region. They found that with an increase in the Ca, the recirculation flow disconnected the GL interface and moved downstream. Brauner et al. [147] followed an exact analytical solution to model the interface curvature of two-dimensional stratified flow. Heil [43] numerically simulated the influence of fluid inertia, Re, on the flow patterns and the propagation of an air bubble on two flexible and rigid walls of a GL flow in a twodimensional channel. Fujioka and Grotberg [148] numerically analyzed the propagation of liquid plugs inside a two-dimensional channel using a finite volume method. They showed that as the liquid plug grew, the frontal part swept the interfacial surfactant from the precursor liquid film. A pair of recirculation regions was found inside the liquid plug shortening as the magnitude of the velocity decreased. They also noted that micelle production and its transport process must be added to the numerical model when plug propagation occurs with a much higher concentration. Their results revealed that as the Re increases, a recirculation flow appears inside the plug core and is then skewed toward the rear part of the bubble. Kreutzer et al. [44] denoted that the inertia effects due to an increase in the Reynolds number elongated the nose and flattened the rear menisci of the gas bubble profile.
The influence of pressure gradient on the flow patterns and bubble profile were illustrated in Figure 9. According to this figure, as the pressure gradient increases, the elongated gas bubble is stretched along the tube axis and the impact of the no-slip conditions at the walls becomes weaker on the bubble shape. Strong, and clockwise circulation is predicted next to the head cap of the bubble resulted in bubble elongation when the pressure is intermediate. Falconi et al. [149] numerically and experimentally studied a vertical-upward three-dimensional Taylor flow in a square mili-channel using the VoF in-house finite volume code and µ-PIV measurements, respectively, Figure 10a. Instantaneous streamlines in Figure 10b show three vortices within the bubble; a large central vortex, two small toroidal vortices at the caps of the bubble. The main central vortex represents the highest z-component velocity next to the channel axis with the direction of rotation of opposite to the smaller vortices. Valizadeh et al. [150] conducted a numerical parametric study of non-Newtonian turbulent flow in a spiral duct by investigating the effects of geometry, consistency index, and power-law index values of viscosity. Since their study was dealt with a wide range of large Reynolds numbers in mili-channel, we do not discuss their findings anymore.
Both experimental and 3D numerical simulations were carried out by Abdollahi et al. [151] for liquid-liquid Taylor flow in a square channel with hydraulic diameters of 1 and 2 mm. Figure 11 displays the non-dimensional contour of the radial velocity component highlighting its key effects on heat transfer enhancement and recirculating flow within the slugs and the carrying phase. The highest radial velocity occurred at the nose and rear of the slug to make two vortices rotating in the opposite direction of the main recirculation zone in the middle of the slugs as observed by [139,149]. In a constant length unit cell, as the dispersed phase volume fraction increased, the liquid film thickness remained constant. The large recirculating zone in the middle of the slug show less axial flow and more radial flow as the droplet length was decreased. Their results also showed a huge increase of heat transfer rate up to 700% compared to single-phase flow indicating more effectiveness of short slugs on the heat transfer rate. Processes 2021, 9, x FOR PEER REVIEW 21 of 45    Xu et al. [152] numerically simulated 3D Taylor flow in a microchannel with a square cross-section under ultrasonic oscillation. A harmonic, vertical, and ultrasonic oscillating applied to microreactor as, d y = A cos (ω 0 t) (9) where the amplitude of the microreactor oscillation (A) is in the order of 2-8 microns and the angular frequency as a function of a constant vibrating frequency of f 0 = 20 kHz is They found the channel oscillation affected the flow pattern and hydrodynamic characteristics in both liquid slugs and bubbles as shown in Figure 12 for Taylor flow in bends and oscillating bubbly flows, for example. Twisted flow structure in the liquid slugs improved mixing rate and transport phenomena. Wavy interface of bubble enhanced the interfacial area accelerating the mass transfer rate across the channel and improved the performance of microreactor. The sub-harmonic bubble surface wave was also created by the pressure pulsation in the parallel direction of the oscillation.
A summary of analytical and numerical investigations on the flow patterns and the formations of the gas/liquid slugs is presented in Table 6. Xu et al. [152] numerically simulated 3D Taylor flow in a microchannel with a square cross-section under ultrasonic oscillation. A harmonic, vertical, and ultrasonic oscillating applied to microreactor as, d y = A cos (ω 0 t) (9) where the amplitude of the microreactor oscillation (A) is in the order of 2-8 microns and the angular frequency as a function of a constant vibrating frequency of f 0 = 20 kHz is They found the channel oscillation affected the flow pattern and hydrodynamic characteristics in both liquid slugs and bubbles as shown in Figure 12 for Taylor flow in bends and oscillating bubbly flows, for example. Twisted flow structure in the liquid slugs improved mixing rate and transport phenomena. Wavy interface of bubble enhanced the interfacial area accelerating the mass transfer rate across the channel and improved the performance of microreactor. The sub-harmonic bubble surface wave was also created by the pressure pulsation in the parallel direction of the oscillation.  Xu et al. [152] numerically simulated 3D Taylor flow in a microchannel with a square cross-section under ultrasonic oscillation. A harmonic, vertical, and ultrasonic oscillating applied to microreactor as, d y = A cos (ω 0 t) (9) where the amplitude of the microreactor oscillation (A) is in the order of 2-8 microns and the angular frequency as a function of a constant vibrating frequency of f 0 = 20 kHz is They found the channel oscillation affected the flow pattern and hydrodynamic characteristics in both liquid slugs and bubbles as shown in Figure 12 for Taylor flow in bends and oscillating bubbly flows, for example. Twisted flow structure in the liquid slugs improved mixing rate and transport phenomena. Wavy interface of bubble enhanced the interfacial area accelerating the mass transfer rate across the channel and improved the performance of microreactor. The sub-harmonic bubble surface wave was also created by the pressure pulsation in the parallel direction of the oscillation.
A summary of analytical and numerical investigations on the flow patterns and the formations of the gas/liquid slugs is presented in Table 6.  A summary of analytical and numerical investigations on the flow patterns and the formations of the gas/liquid slugs is presented in Table 6. Table 6. Selected numerical and analytical studies on the flow patterns and the bubble or the droplet formation.

Comment(s) Cross-Section Phases Reference
The film thickness at the advance of the rear meniscus oscillated The front and rear menisci were slightly different in curvatures Equilibrium bubble profile under surface tension and gravitational effects in a vertical tube Circular GL [30] The film thickness The profiles of front and rear film thickness The film thickness for very elongated gas bubbles Circular and Square GL [35] Bubble profile Flow fields around and between bubbles The motion of a bubble Square GL [145] The influence of velocity on the mass transport phenomenon The effects of liquid slug lengths on flow parameters Circular GL [146] Stronger recirculating flow region as the Ca decreased A liquid backflow was appeared at non-dimensional liquid thickness less than 1/3 Single stagnation point at the vertex of bubble curvature for no-recirculating flow conditions A decrease in the Ca made recirculating flow and moved the stagnation point further the vertex Parallel plates and Circular GL [36] The effects of flow rates on the interface curvature The influence of Eo and wall adhesion on the interface curvature Circular Two-phase (parametric) [147] The rear profile of bubble meniscus versus Ca and Re The effect of Re on the free surface undulations Gas bubble profile The curvature of the gas bubble Circular GL [38] The flow inertia effect was more significant even for deformable wall channels Flexible wall channels were more sensitive than rigid wall channels to the propagation of air bubbles into the walls at low Re and Re/Ca situations The pressure gradient in faraway positions of the tip of the air bubble was generated by Poiseuille flow 2D GL [43] The propagation of liquid plug The adsorption/desorption process of the surfactant was modeled The Marangoni stress results in nearly zero surface velocity at the front meniscus 2D GL [148] Refer to Table 5 Circular GL [44] Slug flow development with time Either uniform or parabolic inlet velocity profile made the same-sized slugs The inlet mixing level influence on the slug lengths 2D GL [58] Bubble shape in slug flow regime The effects of pressure gradient on the bubble profile The presence of gas bubbles made the circulating regions stronger causing a higher momentum transport 2D GL [139] Curved vertical microchannel Slug flow development with time The gas bubbles moved faster than liquid slugs The impact of inlet geometry on the slug development Circular GL [153]

Comment(s) Cross-Section Phases Reference
The VoF and level-set methods The lack of enough adhesion on wall deformed the rear interface of the slugs 2D LL [106] The VoF method Bubble shapes and formation The effects of superficial velocities on the bubble profiles The influence of nozzle diameter on the bubble formation and shapes 2D and 3D GL [154] Structured-square grid minimizes inaccuracies in the surface tension calculation The process of bubble formation was happened periodically 2D GL [155] The VoF and level set techniques Two commercial simulating software; Ansys and TransAT Bubble formation and development The effect of mixture velocity on the flow pattern and the bubble shapes 2D GL [156] Bubble curvature in slug flow The characteristics of flow were dominantly determined by adherent liquid film thickness An optimized model to minimize the pressure drop and maximize the heat transfer rate 2D GL [26]

Flow Regime Transitions
A transition from one flow regime to another often occurs when the flow or boundary conditions are changed. The actual flow maps and the transition conditions specifically depend on a group of effective features during an experiment or the assumptions in numerical simulations [157][158][159][160]. For example, Bottin et al. [161] experimentally obtained various flow regimes for a horizontal two-phase flow in a pipe with an inner diameter of 0.1 m. They attempted to describe different flow regimes regarding the values of superficial velocities of the gas and liquid phases in the earlier experiments [162][163][164][165][166][167]. Some of the transition criteria from a flow regime to another in microchannels are presented in Table 7.
Jayawardena et al. [65] showed the ratio of gas to liquid Reynolds numbers as a function of the Suratman number for recognizing the flow maps. An accurate transitional value of the Su occurs between 10 4 and 10 7 from a bubble to a slug transition. However, the transition from a slug regime to an annular regime has two different criteria regarding the value of Su, greater and less than 10 6 . The motion of the gas bubble through a microchannel with a square cross-sectional area was experimentally studied by Cubaud and Ho [72] and Cubaud et al. [100]. Regarding the liquid fraction, Q l / (Q l + Q g ), they classified the flow patterns into several categories: bubbly, wedging, slug, annular, and dry, which are illustrated in Figure 13. Their experiments specified some critical liquid fractions for a transition from one flow regime to another, which are~0.75 from bubbly to wedging, 0.20 from wedging to slug,~0.04 from slug to annular, and~0.005 from annular to dry transitions. They also found that at low liquid velocity, the gas bubbles clogged at the sharp corner of the channel bend before merging and eventually passing the bend. The corner of the channel bend trapped a gas bubble to merge into another upcoming bubble, producing a single larger bubble.  Table 7. Summary of proposed criteria for the transition of flow pattern.

Comment(s) Criterion Mechanism
Cross-Section Phases Reference Bubbly to slug transition The collision frequency was so low when the void fraction was smaller than 10% The collision frequency was rapidly increased at a void fraction above 25% so that transition to slug flow was rapid even in a strongly liquid When the rate of coalescence was much more than that of breakup and at a void fraction of 0.25 − 0.3

Bubble coalescence and void fraction
Circular GL [168] Churn to annular happens at a value of Minimum gas and liquid superficial velocities Circular GL [92] Slug to slug-bubbly transition No stratified flow pattern was observed in channels with diameter less than 1 mm (it was also verified by [67]) ReWe = 2.8 × 10 5 where, Re = ρ l dV g / 2µ l We = ρ l dV g 2 / 2 The production of Re and the We numbers Rectangular GL [64] The churn to annular transition occurred at lower in comparison to the flow reversal transition the same as [92] The same as [92] Circular GL [66] Bubbly to slug transition 0.25 Void fraction Circular GL [163] The   Table 7. Summary of proposed criteria for the transition of flow pattern.

Comment(s) Criterion Mechanism Cross-Section Phases Reference
Bubbly to slug transition The collision frequency was so low when the void fraction was smaller than 10% The collision frequency was rapidly increased at a void fraction above 25% so that transition to slug flow was rapid even in a strongly liquid When the rate of coalescence was much more than that of break-up and at a void fraction of 0.25 − 0.3

Bubble coalescence and void fraction
Circular GL [168] Churn to annular happens at a value of Minimum gas and liquid superficial velocities Circular GL [92] Slug to slug-bubbly transition No stratified flow pattern was observed in channels with diameter less than 1 mm (it was also verified by [67]) ReWe =2.8 × 10 5 where, Re = ρ l dV g / 2µ l We = ρ l dV 2 g / 2 The production of Re and the We numbers Rectangular GL [64] The churn to annular transition occurred at lower in comparison to the flow reversal transition the same as [92] The same as [92] Circular GL [66] Bubbly to slug transition 0.25 Void fraction Circular GL [163]  Re g = 2 × 10 −9 Su 2 Suratman number Circular GL [65] Bubble to slug transition Upward and vertical flow (N/A) Void fraction Circular GL [173] Bubbly to slug transition The formulation of interfacial transport equations The classification of bubble interactions The categorization of basic mechanisms of bubble coalescence and breakup (N/A) Interfacial area concentrations Circular GL [174] Upward and vertical flow for a bubble to slug transition Void fraction Rectangular GL [175] Discussed in the text Discussed in the text Liquid fraction Square GL [72,100] The bubble size distribution played a crucial role in the flow regime transition and flow patterns (N/A) Population Balance Model (PBM) Circular GL [176] Slug into parallel transition happened at small flow rates and large superficial flow ratio By increasing the flow rate, the chaotic flow was changed into annular Equivalent radius was correlated to the We numbers and holdup fraction

Comment(s) Criterion Mechanism Cross-Section Phases Reference
Bubbly to slug transition The influence of the bubble population at the inlet region on the regime transition An increase in bubble diameter transits the flow map from bubbly to slug A decrease in tube diameter changed the flow patterns from slug to bubbly regimes Bubble diameter PBM Circular GL [177,178] Stable slug regime in the square duct at Re g < 24 Stable slug regime in the circular duct at Re g < 36 Stable slug regime in the rectangular duct at Re g < 40 Stable slug regime in the trapezoidal duct at Re g < 72 Annular flow regime at Reynolds number of gas phase Square and Trapezoidal GL [77] Slug to churn transition and flow map identification α ≥ f α Void fraction Circular GL [179] A transition on the slug producing from squeezing to dripping occurred at the junction Q g > 20 Gas volumetric flow rate Rectangular GL [180] Above the transition value of Ca, an elliptic or oval cross-section area of the gas bubble was observed, and viscous shear overcame the interfacial tension 0.01 Capillary number Quasitrapezoidal GL [80]

Slug Length
The hydrodynamic characteristics of slug flows are highly dependent on the profile and length of the liquid slug. Bubble length measuring was experimentally conducted by Marchessault and Mason [29], where the change in the shape of the liquid film thickness around short bubbles was likely caused by end effects, buoyancy, and cylindrical model assumption. According to Schwartz et al. [32], the liquid film thickness left behind a gas bubble in a capillary was adequately close to the predictions of lubrication theory when the length of bubbles was short. They also found that for the length of bubbles shorter than the tube diameter, the thin film layer region was no longer observed. According to Kashid and Agar [101] and Kashid et al. [102], slug length is significantly affected by the capillary number and the junction dimensions. In a slug flow regime, the lengths of slugs of both phases are greater than their diameters, but at relatively low volumetric flow rates the slug pattern is established with the length equal to the inner diameter of the tube. In the case of a high volumetric flow rate and surrounded water slugs by cyclohexane, the deformed interface regime is made with long water slugs and small cyclohexane droplets. A few simple linear equations to compute the stable slug lengths were reported by [163,170,181,182], which indicated that the ratio of the liquid slug length to the tube diameter remained at a constant value of 16, 20.7, 7.5, and 21.5, respectively. Table 8 presents the correlations for measuring the length of liquid slugs of two-phase flow in a microchannel including the void fractions. Beyond the data provided in Table 8, two other recent studies, Han and Chen [183] and Qian et al. [184], are concerned with the slug length and its formation, but the correlations of slug length were not proposed by the authors. The first study numerically simulated the droplet formation at a T-junction in a microchannel. The effective diameter of droplets was studied in terms of contact angle, continuous phase viscosity, flow rate ratio, and interfacial tension. The latter experimental investigation considered the influence of flow rate on the liquid slug generation in a microchannel with a square cross-sectional area. They found the slug lengths vary from 0.923 to 1.186 mm based on different flow rates of the continuous and dispersed phases. The distance between subsequent water slugs was also measured.
The slug length was extensively depended on physical means of the experiment, operational conditions, and delivering system The squeezing mechanism affected the size of droplets or bubbles The effects of flow rates on the slug lengths w = 0.75 mm h = 0.28 mm 0.08 ms −1 ≤ V l ≤ 0.5 ms −1 0. 155

Pressure Drop
The interactive structure of different phases and between phases and solid walls represent the predominant parameters to describe the flow pattern of multiphase flow and pressure loss through the channel. Two homogenous and separated flow models have been employed by researchers to predict the frictional pressure drop [62]. The first model postulates the same velocity for both gas and liquid phases, which implies that the slip ratio at the interactive boundaries is equal to one. This model considers two or more different phases as a single phase. The values of the flow properties are dependent on the quality, while the frictional pressure drop can be obtained by single-phase theory [63]: where z indicates the axial direction of the channel. The Hagen-Poiseuille equation is a physical law to compute the pressure drop of a Newtonian and incompressible flow through a long and circular tube with a constant cross-section area [63]. The flow regime remained laminar due to the Reynolds number of the microflows, where the friction factor becomes f = 16/Re for round tubes and the equation above can be rearranged as: where L and d are the length and the diameter of the tube, respectively. The presence of the second phase in two-phase flow causes an additional term to the pressure loss of a single-phase flow and does not satisfy Poiseuille's law anymore. Therefore, the total pressure drop (∆p t ) is the sum of the single-phase pressure drop ∆p s and the interaction pressure drop (∆p int ) caused by introducing the secondary phase. These pressure drop components can be shown together as follows: The non-dimensional pressure drop and length scale in Equation (12) were proposed by Walsh et al. [54], as below: L * = L s d (15) where µ and L s denote the dynamic viscosity of the continuous phase and the length of liquid slugs, respectively. The pressure drop over a unit cell can also be described by three components [51,88,107]: Furthermore, Abiev [192] assumed the total pressure losses in the capillary tubes were a summation of two components, named the formation of a new surface during the motion of bubble and the arrangement of a velocity profile in liquid slugs as below (refer to Table 9 for further details): ∆p t = ∆p ∆F + ∆p trans (17) and Numerical [44,53]   Besides, Song et al. [180] expressed the total pressure drop considering the total number of water slugs (n), including the pinching off slug, ∆p t = ∆p g + ∆p cr + (n − 1)∆p s (18) where the first term indicates the pressure drop of single-phase flow of gas, the second term is the critical pressure drop to break up the water slugs, and the last term denotes the pressure drop due to the moving slugs through the capillary. According to the literature, a circular cross-sectional area has been assumed for the gas bubble, where the velocity gradient is limited to the confined liquid region surrounding the bubbles, not inside the bubbles [29]. However, the linear ratio between the pressure and velocity fields is valid for a straight tube, independent of the diameter of the channel. Introducing gas bubbles into the continuous phase flow changes the relation to be nonlinear. The moving bubbles in a tube leave a backfill of liquid left behind the bubble, where the volumetric flow rate of this backflow can be calculated by πrδV g . Marchessault and Mason [29] also showed the pressure gradient as a function of surface tension, film thickness, and the radius of the tube. They experimentally determined that at high velocities the liquid film thickness increases in the direction of the gas-bubble movement, when the profile of the bubble shows a tendency to be conic at the front. Table 9 presents the correlations to calculate the pressure drop created by bubble/droplet in a microflow. One of the first empirical correlations of pressure loss in a capillary tube for horizontal and vertical configurations was experimentally and analytically obtained by Bretherton [30]. He discovered that the pressure drop was enhanced from a coefficient of 14.894 at the leading to 18.804 at the trailing regions due to the slight difference between the frontal and rear menisci of the bubble. Ratulowski and Chang [35] proposed a correlation for pressure drop over a finite gas bubble approaching Bretherton's law for low Ca along with two other correlations for a square tube. Using a Finite Element Method (FVM), Giavedoni and Saita [36] determined the pressure gradient between two specific points; the first at the stagnation point of the gas bubble, and the second in the uniform film thickness region. Their results revealed an increasing trend in pressure gradient with the capillary number. They also found that the backflow region disappears when the capillary number is sufficiently high. Interfacial pressure along the free surface on the rear half of the gas bubble was found by them two years later, where the pressure was significantly decreased because of leaving liquid phase from the uniform film thickness region [38]. Kawaji and Chung [71] found that the shape of the ducts had an impact on the frictional pressure loss or the void fraction. They showed that the total pressure drop was a summation of pressure loss in the single-phase part and the two-phase part of a unit cell. The size impact of the cross-section area of a rectangular microchannel was experimentally carried out by Harirchian and Garimella [104] for a boiling heat transfer situation. A slight decrease in pressure drop was found with a heat flux increase caused by liquid viscosity reduction as the liquid phase temperature increased. They also observed an increase in pressure loss as the cross-sectional area decreased. Niu et al. [105] showed the total pressure drop in a two-phase microflow was the summation of the wall friction and the sudden expansion components showing good agreement when compared with homogenous and separated models.

Conclusions
In this review paper, the literature of GL and LL two-phase flows through microchannels with different cross-sectional areas were discussed. Key results can be summarized in the following way.

•
The maps of Taylor flow and the transition boundaries between each flow regime have been explained by x-y graphs, where the axes are defined by the superficial phase velocities, dimensionless numbers, and volumetric flow rates.

•
Most often the gas bubbles and dispersed liquid droplets are surrounded by a thin film of carrying phase, except in non-circular tubes with very high void fraction where the dry-out at the inner surface of the tube has been observed.

•
There is still no universal agreement to classify all flow regimes due to the experimental apparatus and capturing equipment, which force investigators to name the flow regimes differently.

•
As more efficient numerical solutions and supercomputers emerge, they are used to solve interesting microflow problems with acceptable accuracy. Since the hydrodynamic of the film thickness is crucial and informative, high-resolution images can be produced via CFD. Meanwhile, there is still room to pay more attention to this region because of the important effect of film thickness on the frictional pressure drop, bubble or droplet profile, and heat and mass transfer. • Contact angle can play a major role in the flow pattern as long as two phases are in contact with the inner surface of the tube, except for Taylor flow is, i.e., gas bubbles or dispersed liquid droplets surrounded by a liquid film.

•
The transition criterion from one flow regime into another one has been investigated by many researchers, resulting in numerous criteria based on the nature of the flow, experimental approach, and key parameters.

•
The quantitative attempts to correlate slug length, liquid film thickness (recently summarized by [4] in a comprehensive review paper), and pressure drop have been compared. These correlations include non-dimensional numbers, superficial velocity ratios, volumetric flow rate ratios, void fraction, and other thermophysical characteristics of phases.