Eulerian–Eulerian Modeling of Multiphase Flow in Horizontal Annuli: Current Limitations and Challenges

Abstract

Multiphase flows are present in many natural phenomena, processing technologies, and industries. In the petroleum industry, the multiphase flow is highly relevant, and special attention is paid to the development of predictive tools that determine flow conditions to guarantee safe and economic hydrocarbon extraction and transportation. Hydrodynamic aspects such as pressure drop and holdup are of primary relevance for the field engineer in daily operations like pumping power calculation and equipment selection and control. Multiphase flow associated with oil production is usually a mixture of liquids and gas. The hydrodynamic behavior has been studied in different pipeline configurations (i.e., vertical ascending/descending and horizontal/inclined pipelines). However, the available information about flow patterns as well as the general conditions present in horizontal annuli is incomplete, even if they are of fundamental relevance in today’s horizontal drilling, production, and well intervention in many oil wells around the world. This review aims to present an in-depth revision of the existing models developed to predict two-phase flow patterns and hydrodynamic conditions in annuli flow, focusing mainly on, but not limited to, horizontal configuration. Key flow parameters and effects caused by annuli geometry and the physical properties of fluids are extensively discussed in the present paper. Different empirical correlations and mechanistic and numerical models on two-phase flow through horizontal/inclined pipelines and in both concentric and eccentric annuli are analyzed. Some of these models partially agree with experimental results and show acceptable predictions of frictional pressure loss and flow patterns. Limitations in current models and challenges to be faced in the next generation of models are also discussed.

1. Introduction

An annulus is a space between two pipes, and this space can be found between the production tubing and casing, pipe string and the borehole wall, a gas injector and the production tubing, or between coiled tubing (inside the well) and the production tubing [1]. Multiphase annular flow is commonly found in various applications related to hydrocarbon recovery processes such as drilling, production, and well intervention. Although multiphase flow in annuli is common in the petroleum industry, the understanding of this flow is limited compared to the multiphase flow in circular pipes [2].

In recent decades, the drilling of horizontal wells has become a proven technology for oil and gas recovery. It is particularly useful in production from thin and naturally fractured reservoirs, reservoirs with gas and water conning issues, offshore environments, and various enhanced oil recovery projects [3]. Due to additional challenges associated with the annular geometry, the rigorous prediction of two-phase gas–liquid and liquid–liquid flow behavior is essential to provide the continuous and safe transportation of the production fluid from the downhole to processing facilities. Moreover, in recent decades, researchers have attempted to create correlations and develop models for predicting key-flow parameters such as flow patterns, liquid holdup, and pressure drop in these complex flow conditions with relative success, but there is a significant space for improvements in the models [1,4].

With the computation technologies’ progress, the numerical modeling approach has been in continuous development, and very realistic numerical results can be obtained in complex multiphase flows under appropriate assumptions [5,6].

This paper aims to review the current state of research on multiphase annular flow with a focus on liquid–liquid mixtures through horizontal annuli.

1.1. Empirical and Mechanistic Models in Multiphase Flows: An Overview

Many theoretical and experimental studies have been conducted on two-phase flow in horizontal pipelines and have developed several models. These models can be divided mainly into two groups: empirical models or mechanistic models. Early empirical models disregarded the flow patterns and treated the two-phase flow as it was a one-phase flow [7]. Lockhart and Martinelli [8], Duns and Ros [9], Hagedorn and Brown [10], Wallis [11], and Beggs and Brill [12] studies are among the most important empirical models. In this type of model, the average absolute error in the multiphase-flow pressure-gradient prediction is up to 50.51% [13].

After developing these initial correlations, the focus has been directed to flow pattern identification [7]. Moreover, significant effort has been placed in developing flow pattern maps for two-phase flows in pipelines [14]. Most of the studies on flow patterns transition definitions were conducted using conservation equations, as proposed by Taitel and Dukler [15] and Barnea [16]. Examples from Xiao et al. [17] used the Taitel and Dukler [15] model with some modifications and developed a comprehensive mechanistic model for the two-phase flow in near-horizontal and horizontal pipelines. Petalas and Aziz [18] introduced a mechanistic model applicable to all conditions commonly encountered in the petroleum industry. Moreover, Petalas and Aziz [18] also created empirical correlations applicable to all pipe layout configurations, fluid properties, and flow directions. Their correlations were developed using a large amount of field and experimental data to estimate interfacial friction factors and frictional pressure losses [19]. Caetano et al. [20] developed mechanistic models for bubble, dispersed bubble, slug and annular flows. The models based on physical phenomena of the flow include annulus features such as diameters of casing and tubing as well as the degree of eccentricity. The performance of the models in terms of the average percent error in predicting the pressure gradient in most flow patterns tend to overestimate the gradient in concentric cases and a tendency to underestimate in eccentric cases. In general, the prediction of the pressure drop in mechanistic models has an average absolute error between 30 and 43%.

1.2. Gas–Liquid Two-Phase Annular Flow

Aerated drilling fluid uses a gas–liquid mixture that helps prevent possible reservoir damage, and it is usually used in low-pressure or depleted reservoirs [21]. During the production stage, gas–liquid is introduced through the annulus in flowing wells for high production rates or when implementing gas-lift techniques. Furthermore, extensive studies have been conducted on the gas–liquid annular flow finding different correlations and developing mechanistic models.

Caetano et al. [22] conducted a theoretical and experimental study of upward two-phase flow through vertical concentric and fully eccentric annuli and developed one of the earliest models on the flow-pattern maps for the flow-through annuli. In the experiments, both air–water and air–kerosene mixtures were used. They developed flow pattern maps based on visual observations and mechanistic models for predicting the pressure drop and average liquid holdup for each flow regime in both concentric and fully eccentric annular passage. Moreover, Hasan and Kabir [23] recognized four major flow regimes in vertical and inclined annular flow: bubbly, slug, churn, and annular patterns. They developed the flow pattern approach that predicts the void fraction in the slug, bubbly, and churn flow regimes for air–water systems. Das et al. [24] and Sun et al. [25] experimentally studied air–water upflow through concentric annulus. They investigated the distribution of the two phases in the bubbly, slug, and churn flow regimes using the signals obtained from several conductivity probes and impedance void metters, respectively, located at different places in the annular space. Lage et al. [26] formulated a mechanistic model to predict the two-phase flow behavior in horizontal or slightly inclined fully-eccentric annuli. Their model was validated against the experimental data collected for air–water and diesel oil–N2 mixtures, and the absolute average percentage error of the friction loss was 14.1%, with the average percent error being around −2.3%, indicating a slight tendency to underestimate pressure losses. Ozbayoglu and Yuksel [27] performed an experimental study of the two-phase eccentric horizontal annulus. They used artificial intelligence techniques (like the decision tree, nearest neighbor, and back-propagation neural networks (NN)) to estimate the flow pattern and liquid holdup for air–water flow. The best performance was found with the backpropagation NN with over 90% of correct flow pattern determination, and the nearest neighbor NN technique had the highest accuracy in the estimation of the liquid holdup with 17% of the average absolute error. Osgouei et al. [7] developed a mechanistic model that was used to estimate the frictional pressure losses for a two-phased liquid–gas flow through a horizontal eccentric annular configuration. They estimated the frictional pressure losses with three different methods based on (a) experimental data; (b) the modified Lockhart–Martinelli pressure loss correction factor; and (c) computational intelligence models (regression trees, nearest neighbor, multilayer perceptron, and support vector machines). According to the results, the computational intelligence models had better results than the mechanistic model in predicting the pressure drop with an average absolute percent error of 9.79%, which is a major improvement compared to the modified Lockhart–Martinelli model that can estimate frictional pressure loss with an error of 20%.

1.3. Liquid–Liquid Two-Phase Flow in Pressurized Conduits

The simultaneous flow of oil and water is standard in hydrocarbon recovery. The amount of water in the extracted crude oil increases as wells age [28]; moreover, the importance of this water ratio increase is growing because there are many mature oil wells around the world. Additionally, the liquid–liquid flow is also introduced during the secondary recovery stage, when water is injected to support reservoir pressure. The flow characterization generally includes the following parameters: flow pattern identification and transition, effective viscosity estimation, pressure drop prediction, droplet formation, phase-inversion identification, and droplet size(s) distribution [29]. Trallero et al. [30] noted a large liquid–liquid momentum transfer, small buoyancy of the oil–water flow, and low free energy at the interface, which leads to a small droplet size of the dispersed phase. They also highlighted the diversity of oil properties, e.g., the oil can be more or less dense than water; it can be Newtonian or non-Newtonian; and the viscosity ratio in the oil–water flow can vary by several orders of magnitude. Moreover, Trallero et al. [30] also emphasized that understanding the flow behavior is critical in measuring the amount of free water that is in contact with the pipe, which can promote corrosion or erosion issues. Besides, a correct interpretation of the production logging tools response is dependent on the flow behavior.

Russell and Charles [31] conducted a series of early experimental studies on the liquid–liquid two-phase flow which resulted in the equations for flow rates calculation of two immiscible incompressible Newtonian fluids flowing in laminar regime. Charles et al. [32] studied horizontal oil–water flow with different oil densities, oil and water velocities and oil–water ratios. They identified five different flow patterns; found the dependence of holdup ratios on the contact with the pipe wall of either oil or water; in addition, established a relation between oil–water ratio and pressure drop. Charles et al. served as a reference for subsequent developments in the field [33]. However, despite a large amount of in-depth research of liquid–liquid flow in pipelines [34,35,36,37,38,39,40], fewer efforts have been made on understanding the liquid–liquid flow in the annulus configuration. Nevertheless, the next section is devoted to the relevant investigation in the annulus configuration focusing on its characterization in pressure drop, flow pattern, and liquid holdup.

2. Flow Pattern Characterization in Annular Flow

The flow pattern identification is a critical factor in estimating the pressure drop in multiphase flows accurately [41]. Many studies have developed flow pattern maps for two-phase flow in pipes, but little work has been done for the flow in the annulus [42]. Osgouei et al. [42] noted that the common approach for developing a flow pattern map in annular geometry is the researcher’s interpretation, as there is no standard flow pattern classification method and, consequently, a variation among recognized flow regimes may exist. The description of two-phase flow patterns in annuli is more complicated than for pipe flow since the annulus section has additional characteristics such as the eccentricity and inner/outer diameters along with the parameters that are considered for circular cross-sections [22]. Rouhani and Sohal [43] made one of the first reviews on various aspects of two-phase flows. They discussed experimental methods for the direct and indirect determination of flow regimes, including direct viewing, photography, X-ray and the statistical analysis of fluctuations. In addition, the impact of wall roughness, heat flow and flow acceleration on the transitions between flow regimes are considered in the review. Sadatomi and Sato [44] and Furukawa and Sekoguchi [45] performed experiments on air–water flow in vertical annular passages, which were one of the first studies using two-phase flow patterns in annuli. Sadatomi and Sato [44] found that channel geometry does not significantly affect the flow pattern transitions when the hydraulic diameter (Dh) is greater than 10 mm. On the other hand, Furukawa and Sekoguchi [45] analyzed the void fraction distribution for the air–water flow in the annulus with three different radius ratios and a derived equation for evaluating cross-sectional average void fraction. Kelessidis and Dukler [46] distinguished five major flow patterns in an experiment with air–water two-phase flow regimes associated with the co-current vertical upflow in an annulus: slug, bubbly, churn, annular, and annular-with-lumps.

Besides, Osamasali and Chang [47] conducted experimental studies on air–water two-phase flow in horizontal annuli; furthermore, they inserted rods of different diameters along the central axis of an acrylic pipe and identified stratified-smooth, stratified-wavy, plug, slug, and annular flow patterns. Trallero et al. [30] studied oil–water flow regimes in horizontal pipes and characterized six different flow patterns analyzing the patterns description, holdup, pictures, and pressure loss measurements of the existing works. These flow patterns can be divided into two main categories: stratified and dispersed. The stratified flow patterns consist of a stratified flow with a mixing layer at the interface, and the dispersed flow patterns are divided into the dispersion of oil in water and water in oil, water in oil emulsion, and oil in water emulsion. Wongwises and Pipathattakul [48] identified the following flow regimes for air–water flow in a horizontal concentric annulus: annular, churn, dispersed-bubble, plug, slug, annular/slug, bubble/plug, bubbly/slug-plug, and slug/bubble flow. Their experiments [48] were conducted in a concentric annular test section at the length of 880 mm with an outer diameter of 1.25 cm and an inner diameter of 0.8 cm. Jeong et al. [49] discovered that small bubbles tend to move toward the wall, and large bubbles toward the center during the experimental study of vertical, air–water upward flows in an annular passage. Wu et al. [50] extensively examined 3947 published experimental data points for gas–liquid flow maps in vertical pipes and annuli. They examined existing measurement techniques used to identify bubble, slug, churn, and annular flow regimes. During experimental investigations of gas–liquid flows at the horizontal concentric and fully eccentric annulus, Ibarra et al. [1] identified four main flow regimes together with transitional regions (flow with features of more than one regime) that depend on the range of gas and liquid flow rates, and those flow regimes are dispersed-bubble flow, churn flow, slug flow, and rolling waves (see Figure 1). This classification has been made according to visual observations from high-speed imaging and the cross-sectional holdup data’s probability density function. Ibarra et al. [1] found that, unlike the dispersed-bubble flow regime, the other three flow patterns are commonly more common in the eccentric annuli than in the concentric geometry.

3. Pressure-Drop in Annular Flow

The pressure drop prediction is among the most critical factors in oil and gas pipeline transportation, and it helps to predict the continuity and safety of the hydrocarbon’s conveyance. Due to the cross-section geometry complexity, the pressure loss estimation in the annulus flow is more challenging to the model compared with the circular pipe flow. The common practice to represent the annular geometry is to calculate an effective diameter. Anifowoshe and Osisanya [51] reviewed seven different definitions of equivalent diameter to calculate the frictional pressure losses in an annulus, and they theoretically investigated their effect on wellbore hydraulics. These authors distinguished four widely used equivalent diameter definitions: the hydraulic diameter, slot approximation, Crittendon criteria, and the Lamb approach (see

Table 1. Equivalent diameter definitions (adapted from [51]).

Definition	Formula
The hydraulic diameter	$d_H=d_o-d_i$
Slot approximation	$d_e=0.816\left(d_o-d_i\right)$
Crittendon criteria	$d_e=\frac{1}{2}\left[\sqrt[4]{d_o^4}-d_i^4-\frac{{\left(d_o^2-d_i^2\right)}^2}{\ln \left(\frac{d_o}{d_i}\right)}\right]+\frac{1}{2}\sqrt{d_o^2}-d_i^2$
Lamb approach	$d_e=\sqrt{d_o^2}-d_i^2-\frac{\left(d_o^2-d_i^2\right)}{ln\left(\frac{d_o}{d_i}\right)}$

). However, the hydraulic diameter definition was recognized as the most used in the industry due to its simplicity.

On the other hand, Metin and Ozbayoglu [52] formulated a mechanistic model to predict the gas–liquid flow pattern and frictional pressure losses in both circular and annular geometries. In their model, a “representative diameter” term (see Equation (1)) is introduced that gives up to six times more accurate results for intermittent flow and 4.5 times better accuracy for stratified flow than the hydraulic diameter method for annular geometries:

$$d_r=\sqrt{(d_o^2}-d_i^2)$$

(1)

Additionally, it was found that the liquid phase viscosity, along with the geometry of the pipeline, significantly affects the frictional pressure losses and flow pattern transition, while the liquid density had a small effect on those two parameters [19].

Arirachakaran et al. [34] noted that variations in the pressure gradient could be significantly related to phase inversion or flow pattern transitions. Phase inversion occurs when two immiscible liquids change their continuity in the flow. It is a broad area of study for flows in the dispersed regime. Despite the critical effects associated with phase inversion, its mechanism of generation is the least understood phenomena in dispersions [53]. Ngan et al. [53] also noted that the phase inversion changes the mixture properties like its rheology and drop size. Therefore, the phase inversion point, which is the phase fraction where inversion occurs, needs to be known to process the dispersions. Moreover, Arirachakaran et al. [34] also found that a dramatic change in pressure drop occurs at the phase inversion point and that the magnitude of this change depends on the mixture velocity and oil viscosity. Furthermore, Plasencia and Nydal [54] noted the influence of temperature and chemical composition of the phases on the dispersion’s stability that also affects the phase inversion point. In addition, it has been found that increasing the pipe wall wettability expands the range of phase fractions for the continuous fluid, meaning that the wetting properties of the pipe wall are a relevant factor [55].

Moises and Shah [56] conducted an experimental study of the eccentricity effect on pressure drop calculations for the drilling and the stimulation fluids that flow through annuli. They found that when the annular pressure losses in a fully eccentric annulus are compared to the pressure losses in the concentric annulus, it can be reduced by 60%. Moreover, the effect of eccentricity becomes more evident with the increase in fluid concentration. Pina and Carvalho [57] proposed a model to study the effect of the eccentricity along the well and developed a more precise description of the geometric configuration of directional wells. Moreover, Pina and Carvalho simplified and turned the governing equations into a two-dimensional differential equation using the lubrication theory. Their results show the effect of the maximum and minimum axial velocity in each cross-section, eccentricity on the friction factor, and the presence of azimuthal flow even when the inner cylinder is not rotating.

4. Liquid Holdup in Annular Flow

Experimental studies of two-phase flow in full-pipe systems for various fluid properties, inclinations, and pipe diameters have led to the development of several correlations and mechanistic models; some of them are summarized in

Table 2. Existing holdup correlations for two-phase flow in horizontal and inclined pipes.

Author	Equation
Gregory et al. [58]	$H_{LS}=\frac{1}{1+{\left(\frac{U_M}{\alpha }\right)}^{\beta }}$
Malnes [59]	$H_{LS}=1-\frac{U_M}{83{\left(\frac{g\sigma }{{\rho }_L}\right)}^{0.25}+U_M}$
Marcano et al. [60]	$H_{LS}=\frac{1}{1.001+0.0587U_M+0.0118U_M^2}$
Abdul-Majeed [61]	$H_{LS}=\left(1.009-C{U}_M\right)A$ where the parameter A accounts for the effect of the inclination angle: $C=0.006+1.3377\frac{{\mu }_G}{{\mu }_L}$ A = 1.0          if γ ≤ 0 A = 1.0 − sin φ     if γ > 0
Angeli and Hewitt [35]	${\epsilon }_0=\frac{1}{A}{{\displaystyle \int }}_0^A{\epsilon }_0^{\text{'}}dA= \frac{1}{A}{\displaystyle \sum }_i{\epsilon }_{oi}^{\text{'}}{A}_i$, and ${\displaystyle \sum }_i{A}_i=A$ where Ai is the area of the pipe cross-section surrounding point i: $X^2=\frac{{\left(\frac{dp}{dz}\right)}_0}{{\left(\frac{dp}{dz}\right)}_w}$
Gomez et al. [62]	$R_S=1.0e^{-\left(0.45{\theta }_R+2.48\ast {10}^{-6}R{e}_{LS}\right)}$   0 ≤ γR ≤ 1.57 where the slug Reynolds number is: $R{e}_{LS}=\frac{{\rho }_L{V}_{M}D}{{\mu }_L}$
Zhang et al. [63]	$H_{LS}=\frac{1}{1+\frac{T_{SM}}{3.16{\left[\left({\rho }_L-{\rho }_g\right)g\sigma \right]}^{1/2} }}$ where Tsm has the same units as the shear stress, and it includes both the wall shear stress and the contribution from the momentum exchange between the liquid slug and the liquid film in a slug unit: $T_{sm}=\frac{1}{C_e}[\frac{f_s}{2}{\rho }_s{U}_M^2+\frac{d}{4}\frac{{\rho }_L{H}_{Lf}\left(U_T-U_F\right)\left(U_M-U_F\right)}{l_s}]$ where Ce is the coefficient dependent on pipe inclination angle, and it is calculated with the equation below: $C_e=\frac{2.5-\left|\sin \left(\gamma \right)\right|}{2}$
Al-Safran et al. [64]	$H_{LS}=0.85-0.075\phi +0.057\sqrt{{\phi }^2}+2.27$ where φ is the independent correlation parameter and is calculated with the equation below: $\phi =\left(N_{Fr}{N}_{\mu }^{0.2}-0.89\right)$

.

A liquid holdup is a volume fraction of a fluid present in a given interval of a pipeline; this value cannot be estimated analytically, and it is a function of pipe geometry, flow pattern, and liquid properties [33]. Only a few studies can be found on liquid holdup estimation for the two-phase flow in a horizontal annular configuration [27]. For example, the liquid holdup can be estimated using mechanistic model measurements or calculations through empirical correlations. In addition, Ozbayoglu and Yuksel [27] were one of the first who used an image-processing technique to determine the liquid holdup. Their model measured a holdup in all flow patterns with an average absolute error rate of 0.078 (a fraction of the real value) in the best case. Zhang et al. [65] compared the liquid holdup in high-viscosity oil–gas two-phase flow to low-viscosity liquid–gas flow. They found the liquid holdup in the flow with high-viscosity oil to be considerably higher due to the slowdown of the oil phase. Besides, Ibarra et al. [1] revealed slightly higher holdup values in the eccentric annulus flow compared to concentric cases, which they associated with low velocity due to small clearance between the two pipes. Their experimental data also showed minor differences in liquid holdup between gas–oil and gas–water two-phase flows, except when liquid velocities are low. This difference may be due to the development of different flow regimes affected by geometric parameters and fluid physical properties.

In the experimental study of gas–water and gas–oil slug flow in pipes with concentric and fully eccentric annuli, Ibarra et al. [66] found that the eccentricity has a limited effect on the slug liquid holdup. Nevertheless, in the case of gas–water flow in concentric annular space, a more effective gas encapsulation in the liquid slug is detected due to the chaotic flow behavior. The reason for this behavior may lie in the wetting characteristics of the pipe, as oil tends to create a continuous thin film around the pipe in the gas region, while water does not fully wet the pipe in the gas region for high gas velocities, leading to chaotic behavior. In addition, the slug frequency was investigated and found to be similar for both geometries; however, for gas–oil flows at low gas velocities and medium to high liquid velocities, the slug frequency is considerably higher in the concentric annular space. They also concluded that the experimental data presented in the study supports the effectiveness of full tube models when estimating the gas–liquid flow in the annuli.

5. Numerical Modeling

The numerical modeling of multiphase pipe flow is usually done with one-dimensional (1D) simulation tools in the oil industry [6]. Large-scale models with a lower resolution based on these tools can model the whole system, from the reservoir to the production facilities. One-dimensional tools are practical for the industry as they provide a fast response to technical problems; however, current 1D simulators cannot capture hydrodynamic instabilities of the flow and transient phenomena, which is critical for defining slug flow regimes [6].

In contrast to 1D simulation tools, the computational fluid dynamics (CFD) technique is a great tool that allows us to investigate and understand flow-field characteristics deeply. This approach is widely used in different energy-related sectors as it reduces the uncertainty at the development stage and optimizes engineering designs; moreover, it provides a high resolution of velocity profiles and phase distribution [5,28]. The application of CFD technologies, together with the fast development of computer hardware, has sharply accelerated the investigation of a vast range of industrial problems in recent decades [67].

On the other hand, Eulerian–Lagrangian (E–L) and Eulerian–Eulerian (E–E) approaches are typically used in multiphase flow simulations. The E–L approach can model dispersed phases spread in a continuous phase-like droplets, bubbles, or sand particles. The continuous phase is treated via the Eulerian framework, and the disperse phase is treated via the Lagrangian framework, and the collisions of particles are treated statistically.

The advantage of the E–L approach is that it can provide information on the residence time and behavior of particles for a relatively wide range of particle sizes and velocities. However, for high particle concentration, the E–L approach’s advantages might be neutralized by a prohibitive computational cost [6]. Therefore, the E–L approach is frequently employed in conditions when the number of particles is relatively small, and when it is necessary to move an individual particle [68].

On the other hand, the E–E approach treats both phases as a continuum, and the conservation equations for mass, momentum, and energy are phase-averaged within each control cell, resulting in separate conservation equation sets coupled by phase interaction terms [28]. The velocity and temperature field are calculated for each phase, but both phases share the pressure field. Even though the E–E method can only lead to the spatial distribution of flow parameters, there are fewer calculation requirements; therefore, there are no restrictions on the volume or the number of fractions of dispersed phases, which makes this method more effective in large dispersions [68]. However, the E–E approach becomes computationally expensive when simulating many particle sizes/classes, as it is required to solve conservation equations for each of the sizes/classes [6].

Several numerical studies have been conducted on multiphase flow in horizontal annular geometry. For example, Sorgun [69] used CFD simulations with experimental validations to investigate the effects of pipe eccentricity on axial and tangential velocity, frictional pressure loss, and the effective viscosity of water and non-Newtonian fluids; moreover, Sorgun [69] found that when the pipe eccentricity rises, the frictional pressure loss inside the annulus is reduced. Furthermore, Sorgun [69] found that the tangential velocity for both Newtonian and non-Newtonian fluids increases dramatically. It was also noted that the effect of pipe eccentricity on tangential velocity became more severe with increased fluid viscosity. In addition, results suggest that the axial velocity and effective viscosity inside the annulus increase with the increase in the pipe eccentricity. On the other hand, McCaslin and Desjardins [14] conducted a numerical examination of the gravitational effects in horizontal liquid–gas flows from first principles, and they presented the governing parameters through dimensional analysis and described the behavior of three simulations using those parameters. They deduced that the gravitational effect has more impact than inertial effects. Busahmin et al. [70] proposed a numerical model for hole cleaning in the horizontal well-bores. When numerically solving the cutting transport equation to determine the cutting concentrations it was found that the concentration of cutting increases with an increasing rate of penetration and decreases with increasing drilling mud flow rate and time. Erge et al. [71] applied experimental, analytical, and numerical approaches to study the dependence of the annular pressure loss on drill-pipe eccentricity during the circulation of non-Newtonian drilling fluids, and they found that the effect of eccentricity on pressure loss has been under-estimated by all the models. Moreover, Sultan et al. [72] developed a numerical model of multiphase (gas–liquid–solid) flow and validated it by comparing CFD simulation results with previously collected experimental data and empirical correlations in both pipeline and annular geometries. Their model predicted that the average difference of pressure loss for the different cross-sectional areas of the pipelines and annuli is usually less than 15%, with a maximum value of 30%. Friedemann et al. [73] conducted numerical simulations of gas–liquid flow in a horizontal concentric annulus in 3, 5, and 7 m long domains. Moreover, they used the experimental results gathered at the Institute for Energy Technology in Norway as a comparison, and, for all cases, the maximum error considering any of the extreme values or the mean value was within 30%.

Friedemann et al. [74] ran two-phase flow simulations in an eccentric annulus with the OpenFoam CFD platform using periodic boundary conditions. They compared numerical data with the results of experiments conducted in a fully concentric/eccentric geometries, and they found that the numerical results of the holdup pattern corresponded within reasonable accuracy to the experimental data. However, the change of eccentricity led to a significant discrepancy in the pressure gradient results. The annulus eccentricity in simulations was set to 0.5, while the experiments were at E = 0 and E = 1.0.

5.1. Modeling of Turbulence

The CFD’s accuracy in modeling multiphase flow strongly depends on the approach used to model the turbulence. One of the most precise approaches for turbulent multiphase flow modeling is the direct numerical simulation (DNS) that solves the full three-dimensional, time-dependent Navier–Stokes equations without using a turbulence model. Furthermore, the absence of a turbulence model assumes a numerical solution of all the spatial and temporal scales of turbulence, resulting in its unrivaled accuracy compared to other methods. Although the DNS is great for the study of the enhancement and development of turbulence models, turbulence mechanisms, and the evaluation of point-to-point closure theories, the method is too expensive for practical flow conditions. Furthermore, DNS usage is limited to flows of low to moderate Reynolds numbers even with modern supercomputers [75].

Reynolds-averaged Navier–Stokes (RANS) models are time-averaged equations of fluid flow motion. These equations govern the averaged flow quantities’ transport, within the full spectrum of the turbulence scales, and significantly reduce the computational effort and resources required. It is the most practical and widely used approach for industrial flows. The RANS equations are frequently applied to compute time-dependent flows. The unsteadiness of these flows may be imposed externally or self-sustained [76].

The traditional linear-eddy-viscosity RANS models can be categorized into the following four main groups:

• Algebraic (zero-equation) models: it is the oldest class, and it uses partial differential equations to calculate only the mean fields, and for the turbulence values, it only uses algebraic expressions. Currently, the model is applied only for the initial prediction of the flow field before more accurate models are explored [77].
• Half-equation model: it is a two-layer model for studying separated flows controlled by pressure (also known as the Johnson and King model) [75].
• One-equation models: these models can be distinguished because they formulate one extra transport equation to compute the magnitude of the turbulence, usually the kinetic energy (k). For all of the one-equation models, there is still a need to assign a distribution along the length scale (L), which is algebraically determined and is commonly based on available experimental data. The one-equation models have been developed to enhance turbulent flow predictions by solving one additional transport equation [75].
• Two-equation models: besides the mean-flow Navier–Stokes equations, the models use two transport equations for two turbulence properties. The first property is usually the turbulence kinetic energy (k), while the second varies according to the model, and among the properties that can be used are the dissipation rate of turbulent kinetic energy (ε), the length scale (l), and the specific dissipation rate (ω). These models are preferred by the industry, being the first choice for general CFD calculations, with the most widely used standard models being k–ε (Launder and Sharma [78]) and k–ω (Wilcox [79]) [70].

Moreover, the turbulent k–ω shear stress transport-model (SST) introduced by Menter [80] merges the advantages of k–ε and k–ω models. Thus, in free turbulent flow, the Menter model behaves like k–ε, thus avoiding the common issues associated with the k–ω model, which is too sensitive to the initial conditions of turbulent flow [81].

Table 3. Turbulence models used in two-phase flow computational fluid dynamics (CFD).

Turbulence Model	Comment
k–ε model	It is the most widely tested and used two-equation model; Inability to accurately predict flows with unfavorable pressure gradients and additional deformations. It can lead to poor results for separated flows.
k–ω model	Achieves higher accuracy for boundary layers; More accurate for free shear flows and separated flows.
Shear stress transport (SST)	Near-wall regions are evaluated with k–ω, and free-streams are calculated with k–ε.

presents a summary of the most popular CFD turbulence models used in two-phase flows.

Large eddy simulation (LES) is another modeling approach that fully resolves large scales of the turbulent flow while the small sub-grid scales are modeled. The LES approach is especially practical for studying turbulence at high Reynolds numbers, the prediction of complex flows, development, and the evaluation of new turbulence models. The LES’s impact on industrial flow modeling is moderate and mostly limited to free shear flows, even if this method has been around for many decades. However, the LES model’s use is limited as most engineering applications use wall boundary layers, and the cost of applying LES in flow-modeling in these regions is too high [82].

5.2. Governing Equations

When simulating two-phase flow by applying the E–E approach, the following governing equations are solved:

The mass continuity equation:

$$\nabla •\left(U_{\alpha }\right)=0$$

(2)

Conservation of momentum equation:

$$\nabla •\left(r_{\alpha }\left({\rho }_{\alpha }{U}_{\alpha }\otimes U_{\alpha }\right)\right)=-r_{\alpha }\nabla p+\nabla •\left(r_{\alpha }{\mu }_{\alpha }\left(\nabla U_{\alpha }+{\left(\nabla U_{\alpha }\right)}^T\right)\right)+r_{\alpha }{\rho }_{\alpha }g+M_{\alpha }$$

(3)

where p is the static pressure of either phase, and M_α is the sum of interfacial forces like the lift force, drag force, buoyancy, and wall lubrication force [83].

The primary factor that includes the interaction between disperse and continuum phases is an interfacial drag force, which is part of M_α (Equation (3)). The interfacial drag force is calculated using Equation (4):

$$F_D=\frac{3}{4}{C}_D\frac{r_{\beta }{\rho }_{\alpha }}{d_{\beta }}\left|U_{\alpha }-U_{\beta }\right|\left(U_{\alpha }-U_{\beta }\right)$$

(4)

where C_D is the drag force, which is correlated in terms of the non-dimensional drag coefficient (Equation (5)):

$$C_D=\frac{D}{\frac{1}{2}{\rho }_{\alpha }{\left(U_{\alpha }-U_{\beta }\right)}^{2}A}$$

(5)

where D is the drag force magnitude, α is the continuous phase, β is the dispersed phase, (U_α − U_β) is the relative speed, and A is the projected area of the body [76]. Due to the deformable shape of bubbles/droplets, several semi-empirical equations can be used to adjust the shapes of an ellipse, sphere, and spherical caps that can be applied to calculate the drag coefficient.

In addition to the drag acting in the direction of flow, the so-called non-drag forces act mainly perpendicular to the direction of the flow and should also be considered; for example, the lifting force, the wall force, and the turbulent dispersion force [84]. The lifting force accounts for the interfacial momentum exchange between the bubbles and the liquid field while the wall force considers the hydrodynamic impulse acting on the bubble moving in the immediate vicinity of the solid wall [85]. The lifting force effect is the movement of bubbles towards or away from the pipe wall and the wall force effect is generally a sharp spike in the void fraction of the bubble near the wall [85].

Chuang and Hibiki [86] reviewed the progress made in the development of an interphase force model and the recent advancement in an interphase transportation model. Moreover, Chuang and Hibiki extensively reviewed the progress in modeling the interfacial force, presenting the recent advance of the interfacial area transport equation. They pointed out several concerns with the accuracy of the current two-phase flow CFD models. For example, they mentioned that there were limited data for validation, noticing that existing two-phase flow data mostly have been taken in relatively simple flow channels with low-pressure conditions (near atmospheric). Moreover, they emphasized the importance of creating a code benchmark process. Additionally, Chuang and Hibiki [86] discussed both the framework of the bubble-wall collision force and the structure of the bubble collision force while considering the effect of bubble coalescence on the bubble collision frequency.

5.3. Coalescence and Breakup

In recent decades, great effort was given to study the phenomenon of fluid droplets/bubbles coalescence and breakup [87]. Moreover, the resulting theoretical and experimental studies seek to determine the collision and coalescence efficiency [88,89,90], as well as the breakup frequency and daughter size distribution function [91,92,93,94].

The treatment of the coalescence and breakup mechanisms and their effects on the energy, species, and momentum transport processes is one of the most significant challenges in modeling multiphase flow [85]. Understanding the interaction between phases is important, as there is a strong relation of the mean motion of two-phase flow with the small-scale motion. These interactions are commonly dependent on the flow conditions, fluid properties, and geometry of the flow conduit [85]. Liao and Lucas [95] stated that the breakup mechanism, on the whole, can be stated as the balance between external stresses from the continuous phase and which, as the external stresses try to destroy the fluid particle, the surface stress of the particle, along with the viscous stress of the fluid, restores its form. Consequently, the hydrodynamic conditions in the surrounding liquid and the fluid characteristic determine the breakup of a particle. Liao and Lucas [95] divided the breakup mechanisms into four categories: interfacial instability, viscous shear stress, shearing-off process, and turbulent fluctuation and collision.

There are several ways of interaction recognized between dispersed droplets/bubbles and turbulent eddies within the continuous phase depending on their relative size [35], and those interactions can change the key parameters of multiphase flow.

Montoya et al. [85] categorized the breakup and coalescence mechanisms expected to affect a two-phase flow regime as follows:

• Fluid particle coalescence caused by random collisions due to turbulent eddies;
• Fluid particle coalescence caused by wake entrainment;
• Collision due to velocity gradient near the wall region;
• Collision due to different rise velocity of fluid particles with different sizes;
• Large-cap and churn-turbulent fluid particles’ breakup driven by flow instabilities;
• Fluid particles breakup due to the turbulent eddy impact;
• Small fluid particles shearing-off at the rim of large-cap, churn-turbulent, or slug fluid particles;
• Breakup caused by a laminar viscous force.

The most recognized and commonly used models that can be found in the literature are the Prince and Blanch [88] model for bubble breakup and the Luo and Svendson [92] model for coalescence. These models have been employed either with their original formulation or with certain modifications in several important studies [96]. However, these modifications were ineffective in the prediction of high void fraction regimes due to the insufficient consideration of the whole range of mechanisms for breakup and coalescence [85].

5.4. Models for Breakup Frequency

Coulaloglou and Tavlarides [97] developed the pioneer breakup kernel that is applicable for liquid–liquid dispersions. The main idea is that a deformed oscillating drop will break if the drop surface energy is lower than the turbulent kinetic energy transmitted from drop-eddy collisions. The breakup rate is defined as

$${\mathsf{\Omega }}_B\left(V\right)=\left(\frac{1}{breakage time}\right)\left(\frac{fraction of}{drops breaking}\right)= \frac{1}{t_B}\frac{\mathsf{∆}N\left(V\right)}{N\left(V\right)}$$

(6)

where N is the total number of drops, and t_B is the breakup time.

The drops’ kinetic energy distribution is presumed to be the same as that of the turbulent eddies. Moreover, assuming that the distribution of kinetic energy is a normal function, the equation for breakage frequency of a droplet is expressed as

$${\mathsf{\Omega }}_B\left(V\right)= C_I{V}^{-2/9}{\epsilon }^{1/3}\exp \left[-\frac{C_{II}\sigma }{{\rho }_d{V}^{\frac{5}{9}}{\epsilon }^{\frac{2}{3}}}\right]$$

(7)

where C_I and C_II are breakup constants, σ is the interfacial tension, and ρ_d is the dispersed phase density.

In the resulting equation, the breakup rate is presented as a function of the dispersed phase density. Prince and Blanch [88] compared Equation (7) with the experimental data for gas–liquid systems and revealed that the breakup rate values were several orders of magnitude lower than experimental observations. Moreover, Prince and Blanch [88] established a bubble breakup model based on the continuous phase density. They drew an analogy to collisions in ideal gases and studied the bubbles and turbulent eddies interaction to define the breakup frequency. In addition, they calculated the breakup rate similarly to the coalescence rate by multiplying the collision frequency with a collision efficiency, which is defined as a function of the time bubbles remain in contact and the time required for bubbles to coalesce [95]. The total breakup rate for all bubbles is determined with Equation (8):

$${\mathsf{\Omega }}_{BT}= {\displaystyle \sum }_i{\displaystyle \sum }_e{\mathsf{\Gamma }}_{ie}\exp \left(-\frac{u_{ci}^2}{u_{ie}^2}\right)$$

(8)

where Ω_BT is the total breakup rate, Γ_ie is the collision rate, u_ci is the critical eddy velocity (i.e., eddy velocity required to cause a bubble’s break-up), and u_te is the turbulent eddy velocity.

Prince and Blanch [88] also noted the need to obtain data about the bubble concentration and bubble size distribution to use the model:

$$n_i= \frac{\Phi R_T^2{H}_T}{\frac{4}{3}{\overline{r}}_{bi}^3}{x}_i$$

(9)

where x_i is the fraction of bubbles, and r_bi is the bubble’s radius.

Similar to Prince and Blanch [88], Luo and Svendsen [91] proposed a novel theoretical model based on the kinetic gas theory to calculate breakup frequency (Equation (10)). They pointed out the cost of the experiments required to determine unknown parameters, thus criticizing all previous models [95]. To develop the model, they assumed the following simplifications [92]:

• The turbulence to be isotropic;
• Consider only the binary breakage of fluid particles in a turbulent dispersion;
• The breakage volume fraction to be a stochastic variable;
• Determine the occurrence of breakup by the energy level of the arriving eddy. Moreover, the frequency of the particle oscillation to be greater than the eddies’ arrival frequency;
• The oscillations of a particle can only be caused by eddies of a length scale smaller than or equal to the particle diameter:

  $${\mathsf{\Omega }}_B\left(V:V{f}_{BV}\right)={{\displaystyle \int }}_{{\lambda }_{min}}^d{P}_B\left(V:V{f}_{BV},\lambda \right){\mathsf{ώ}}_{B,\lambda }\left(V\right)d\lambda$$

  (10)

  where ώ_B,λ(V) is the arrival frequency of eddies, λ is the eddy size (length-scale), V is the particle size (volume), P_B(V:V_I,λ) is the probability that particles of size V break into two when an arriving eddy hits the particle, and f_BV is an independent variable with a range of (0.1) [92].
  Given the random nature of the eddies motion, the collision frequency can be expressed as

  $${\mathsf{ώ}}_{B,\lambda }\left(d\right)=\frac{\pi }{4}{\left(d+\lambda \right)}^2{\overline{u}}_{\lambda }{\mathsf{ṅ}}_{\lambda }n$$

  (11)

  where ṅ_λ is the number of eddies per unit reactor volume, d is the particle diameter, and $\overline{v}$_λ is the turbulent velocity of eddies.
  The conditional probability for a particle to break up (P_B(V:Vf_BV,λ)) can be expressed as

  $$P_B\left(V:V{f}_{BV},\lambda \right)=1-{{\displaystyle \int }}_0^{X_c}\exp \left(-\chi \right)d\chi =\exp \left(-{\chi }_c\right)$$

  (12)
  The total breakup rate obtained by integrating the resulted equation for the breakup rate of particles over the whole interval can be expressed as

  $${\mathsf{\Omega }}_B\left(V\right)=\frac{1}{2}{{\displaystyle \int }}_0^1{\mathsf{\Omega }}_B\left(V:V{f}_{BV}\right)d{f}_{BV}$$

  (13)

  where 1/2 factor considers the effective range of f_BV = 0–0.5 or 0.5–1.

Luo and Svendsen’s model is widely used in later applications. It does not include any unknown or empirical parameters; moreover, it also normalizes the partial breakup frequency by the total breakup frequency, which allows the calculation of the daughter distribution directly [95].

5.5. Models for Coalescence Frequency

The coalescence is commonly divided into three steps: (1) two bubbles collide, and there is a small quantity of liquid trapped in between; (2) the bubbles make contact until the liquid film reduces to a critical thickness; (3) the liquid film breaks resulting in bubble coalescence. Coalescence only occurs if the bubbles interact for enough time for the intervening liquid film to reach the critical rupture thickness.

There are three main theories/criteria proposed for the coalescence process: the critical value of the relative velocity, the film drainage model, and the critical approach velocity model [87]. In any case, the conditions for coalescence are contact and collision.

Over the years, several models have been used to calculate the collision frequency and coalescence efficiency; some used an empirical approach, while some are based on physical quantities to determine each collision’s mechanisms.

5.6. Collision Frequency

Liao and Lucas [87] presented a comprehensive literature review on the mechanisms and models of coalescence of fluid particles. They have identified five primary sources of motion leading to the collision in turbulent bubble flow: buoyancy, eddy capture, global velocity gradients, turbulent fluctuations, and wake effect. Moreover, one of the primary reasons for bubble collisions are the fluctuations in the turbulent velocity of the liquid phase, and the collision is commonly assumed to be alike to particle collisions in an ideal gas (Prince and Blanch [88]). However, compared to the turbulent fluctuations, the other significant sources of motion leading to collisions are usually neglected in the research.

Nevertheless, Prince and Blanch [88] developed a phenomenological model for the bubble coalescence rate in turbulent gas–liquid dispersions. The model considers bubble collisions caused by the buoyancy, laminar shear, and turbulence.

5.7. Turbulent Collision Rate

The collision rate was derived following Kennard [98] and expressed as a function of bubble size, concentration, and velocity as

$${\mathsf{\theta }}_{ij}^{\mathrm{T}}=n_i{n}_j{S}_{ij}{\left({\overline{u}}_{ti}^2+{\overline{u}}_{tj}^2\right)}^{1/2}$$

(14)

where n_i and n_j are concentrations of bubbles of radius r_b, and r_bj, respectively, ū_t is the average turbulent fluctuating velocity of the bubble, and S_ij is the collision cross-sectional area of the bubbles.

There are also assumptions required to determine the turbulent velocity of a bubble: (1) the turbulence is isotropic, and (2) the bubble size lies in the inertial subrange.

Prince and Blanch [82] expressed the turbulent velocity as (Rotta [99]):

$$u_t=1.4{ϵ}^{\left(\frac{1}{3}\right)}{d}_b^{\left(\frac{1}{3}\right)}$$

(15)

The energy dissipation per unit mass was derived by dividing the expression for the power input to gas sparged vessels developed by Bhavaraju et al. [100] by the mass of the system:

$$\epsilon =\frac{Qg}{\pi R_T^2}\frac{P_2\ln \left(\frac{P_1}{P_2}\right)}{\left(P_1-P_2\right)}: P_1=P_2+{\rho }_{l}g{H}_T$$

(16)

Subscript 1 and 2 correspond to the position in the column. During the experiment, bubble samples were taken at two different points in the column, 80 and 180 cm above the gas bubblers.

Substituting all these values in Equation (14) allows us to calculate the turbulent collision rate:

$${\mathsf{\theta }}_{ij}^{\mathrm{T}}=0.089\pi n_i{n}_j{\left(d_{bi}+d_{bj}\right)}^2{\epsilon }^{1/3}{\left(d_{bi}^{2/3}+d_{bj}^{2/3}\right)}^{1/2}$$

(17)

5.8. Buoyancy-Driven Collision Rates

The buoyant collision rate is given by Friedlander [101]:

$${\theta }_{ij}^B=n_i{n}_j{S}_{ij} \left(u_{ri}-u_{rj}\right)$$

(18)

The rise velocity can be expressed as a function of size (Clift et al. [102]):

$$U_r={\left[\left(2.14\sigma \right)/{\rho }_l{d}_b+0.505g{d}_b\right]}^{1/2}$$

(19)

5.9. Laminar Shear Collision Rate

The functional form of the collision rate due to laminar shear is given by Friedlander [101] as

$${\mathsf{\theta }}_{ij}^{\mathrm{LS}}=n_i{n}_j\frac{4}{3}{\left(r_{bi}+r_{bj}\right)}^3\left(\frac{\overline{d{U}_l}}{dR}\right)$$

(20)

where U is the liquid circulation velocity, and R is the radial coordinate of the column.

The velocity profile for this condition is:

$$U_l=U_{l max}\left[\frac{1-R^2}{{\left(\alpha R_T\right)}^2}\right]$$

(21)

where αR_T is the transition point, U_lmax is the velocity in the center of the column.

The maximum circulating velocity is predicted using the Miyauchi and Shyu [103] model:

$$U_{l max}=\left(\frac{1-0.75\Phi }{1-\Phi }\right)\frac{\Phi g{D}_T^2}{48{\upsilon }_t}$$

(22)

where Φ is the gas holdup, and υ is the turbulent kinetic viscosity.

5.10. Collision Efficiency

Existing experimental studies indicate that the majority of collisions do not lead to coalescence and colliding bubbles separate from each other after the collision. Thus, a collision efficiency parameter must be introduced, which is a function of the contact time that passes while bubbles touch and the time required for bubbles to coalesce. The collision efficiency is given by Coulaloglou and Tavlarides [97]:

$${\lambda }_{ij}=\exp \left(-\frac{t_{ij}}{{\tau }_{ij}}\right)$$

(23)

where t_ij is the time required for bubbles’ coalescence, and τ_ij is the contact time of two bubbles.

To avoid the model’s cumbersomeness, Prince and Blanch [88] used a simplified model for calculating the thinning of the film without considering the Hamaker contribution since it is only noticeable at very low film thicknesses before the bubbles’ rupture. The coalescence time was determined by applying the approach velocity between the bubbles.

Thus, the coalescence rate of bubbles is expressed by the total collision frequency multiplied by the efficiency:

$${\mathsf{\Gamma }}_{ij}=\left\{{\theta }_{ij}^T+{\theta }_{ij}^B+{\theta }_{ij}^{LS}\right\}\times \exp \left\{-\frac{t_{ij}}{{\tau }_{ij}}\right\}$$

(24)

The overall coalescence rate is then:

$${\mathsf{\Gamma }}_T=\frac{1}{2}{\displaystyle \sum }_i{\displaystyle \sum }_j\{\left({\theta }_{ij}^T+{\theta }_{ij}^B+{\theta }_{ij}^{LS}\right)\times \exp (-\frac{t_{ij}}{{\tau }_{ij}})\}$$

(25)

where the 1/2 factor is included to avoid the double counting of coalescence between bubble pairs.

The modeling of bubbles’/droplets’ coalescence and breakup requires the solution of a population balance equation. There are many model concepts in the literature based on the two-fluid model approach, and one of the most widely used is the discretized population balance method; nevertheless, it requires significant extra calculations [84]. Wang et al. [90] performed a comparison of the predicted bubble size distributions using four common bubble breakup and coalescence models (Lehr et al., Luo and Svendsen, Prince and Blanch, and Wang et al.), and they found that the use of different models resulted in considerable distinction in the prediction of both the bubble size distributions and the flow regimes transition [68].

5.11. Multiple-Size Group Approach

Since the modeling of breakup and coalescence mechanisms are the weaknesses in simulations of poly-dispersed flows [15,84], a multiphase flow modeling method called the multiple size group (MUSIG) approach was implemented. In the MUSIG approach, the dispersed phase has a wide variety of bubbles/droplet sizes, and it is based on the Eulerian–Eulerian framework. It provides an efficient method for solving the conservation equations, closure equations for the interfacial transport of mass, energy, momentum, population balance equations, and bubble/droplet coalescence and breakup [104].

The MUSIG model was developed to process polydisperse multiphase flows, providing a framework for the inclusion of the population balance model into three-dimensional calculations [105]. The model was first proposed by Lo [106], with only one common momentum equation solved for all size groups and the diameter range of the dispersed phase divided into M size (homogeneous MUSIG). Mathematically, the MUSIG model is based on the two-fluid modeling approach and the population balance method. Furthermore, the population balance equation is applied to account for the mass conservation of the size fractions while considering the inter-fraction mass transfer resulting from bubble coalescence and breakup [84]. Due to the accurate calculation of coalescence and breakup for big-sized groups, the multiple size group model has been applied successfully on several occasions. However, the homogeneous model’s limitation was soon discovered, as it assumes that the slip velocities of bubbles are independent of bubble size [107].

Moreover, it is necessary to consider different velocity fields for bubbles of different sizes, and one of the ways to do it is using a multi bubble-size class test solver developed by Lucas et al. [84,108]. As a result of these studies, the concept of an inhomogeneous MUSIG model based on the Eulerian modeling framework proposed by Krepper et al. [109] was incorporated into CFD analyses, which is an enhancement of the existing homogeneous MUSIG model. Within this inhomogeneous MUSIG model, the dispersed gaseous phase is divided into N inhomogeneous velocity groups (phases) with a particular velocity field for each of the velocity groups, and the bubble diameter range within each of the velocity groups being subdivided into Mj bubble size classes. The corresponding models consider the processes of bubbles’ breakup and coalescence between all classes of bubble sizes Mj [84].

Huang et al. [110] performed a numerical study for three bubble-column ascending configurations occurring in the churn turbulent flow regime in a vertical circular column. They examined the effect of various bubble size models on numerical predictions. The results showed no improvement in homogenous nor inhomogeneous MUSIG models compared to the single bubble size models. Thus, Huang et al. [110] stated the need to improve breakup and coalescence closure descriptions and have experimental data to validate such models.

On the other hand, Lucas et al. [111] proposed an extension of the inhomogeneous MUSIG model, including the phase transfer effects. They stated that additional source and sink terms related to evaporation or condensation must be added to the continuity equation to simulate poly-dispersed flows with phase transfer. In addition, Lucas et al. [111] specified three phenomena to consider:

(a)

Shrinkage or expansion of bubbles caused by temporal and spatial pressure variation.

(b)

Mass transfer between the bubble size group and liquid phase.

(c)

Shrinkage or expansion of bubbles caused by mass transfer.

6. Concluding Remarks

Even though extensive theoretical and experimental studies have been conducted on multiphase annulus flow, currently there are still questions and challenges that need to be addressed. Some of these limitations and challenges in the study of annular flow are listed below:

-

Despite several available studies on two-phase flow at circular pipes, there is still a lack of studies on the two-phase flow through the annulus [7].

-

The need to develop correlations that work with several flow characteristics and fluid properties, as well as a wide range of velocities [112].

-

More attention should be paid to interfacial friction coefficient, entrainment fraction, effective viscosity, oil–water dispersion inversion, gas trapping in liquid [65].

-

The lack of general agreement in available flow pattern maps due to the diversity of oil properties [28].

-

Research on the slug characteristics in two-phase gas–liquid flows in annuli is limited, particularly in nearly horizontal geometries [66].

-

The lack of physics-based models in the commercial dynamic multiphase flow simulators and lack of industry-relevant high-quality multiphase flow data [1].

-

The effect of using different closure laws needs to be further studied. Moreover, the investigation of the CFD model capability to predict more complex flow patterns is also required [28].

-

The investigation of the eccentricity effect and presence of the secondary phase’s variable concentration upon the effect of the transition needs a proper consideration [113].

-

It is still challenging to consider and simulate a realistic number of particle size classes, which would increase simulation accuracy, using either Eulerian–Eulerian or Eulerian–Lagrangian frameworks [6].

-

Computational limitations still prevent an in-depth study that combines a fully turbulent gas phase with the effects of a deformable interface and non-unity density and viscosity ratios of a liquid–gas two-phase flow inside a horizontal pipe [28].

-

Improvements in the modeling of breakup closures must be made [95]. Moreover, most of the existing breakup closure equations are at their initial development stage and need improvement.

-

More fundamental and consistent coalescence models are required, which can be used in practice in the vast flow conditions and will consider all significant collisions [87]. Moreover, extensive experimental work needs to be conducted for a better understanding of the microscopic coalescence process.

-

The closure model for breakup and coalescence is the weakest point, limiting the application of the MUSIG approach [107].