# A Three-Dimensional Model of Turbulent Core Annular Flow Regime

Authors: Saliha Nouri; Zouhaier Hafsia; Salah Mahmoud Boulaaras; Ali Allahem; Salem Alkhalaf; Baowei Feng
Journal: Journal of Mathematics
 (2021)

---

## Abstract

In this study, three-dimensional (3D) turbulent core annular flow (CAF) regime is investigated numerically. The proposed model is based on the 3D Reynolds average Navier–Stokes (RANS) equations combined with a pure convective transport equation of the volume of fluid (VOF) to predict the interface between the oil and water phases. The k-ω turbulence model is adopted to better reproduce the oil and water flow characteristics. The two-phase (CAF) regime can be predicted by two inlet configurations: the T-junction (3D-T) and the straight pipe (3D-S). These two configurations are simulated and compared for pipe diameter D=0.026m and pipe length L=4m. For these two inlet configurations, the computed mixture velocity profile and the water volume fraction at a test section z=100D were compared to experimental measurements. The 3D-T configuration gives more appropriate results. The 3D-S slightly overestimates the maximum velocity at the test section and the lower and upper water layer of the (CAF) flow is shifted in the upward direction. For the 3D-T, the relative error in the pressure drop is 3.3%. However, for the 3D-S, this error is 13.0%.


---

## Body

## 1. Introduction

Two-phase oil-water flow is a particular case of multiphase flow and has a great interest especially in the petroleum, nuclear, and chemical industries. The distribution of each phase affects the pressure drop in a more complicated way than in a single-phase flow. The main tasks of the numerical and experimental investigations are the predetermination of the characteristics of the flow regimes, water holdup, and pressure drop along the transport cylindrical pipe. Analytical solutions were restricted to the formulations of the pressure drop (Brauner in [1]). Following the prescribed superficial velocities of the two phases, there are mainly different flow regimes: stratified, bubble, dispersed, plug, slug, and core annular flow (CAF). The slug and the stratified flow have been extensively studied experimentally and numerically (Cherlo et al., Desir et al., and Al-Yaari and Abu-Sharkh in [2–4]). Recently, the CAF regime has become the more interesting way for the transportation of oil.The experimental measurements are the primal method to study the distribution of the oil and water phases under different superficial velocities of the two phases. Charles et al., in [5], conducted the first experimental study of oil-water flow in a uniform cross section and horizontal pipe with equal densities and low oil viscosity. In these experiments, the CAF was observed. Elseth, in [6], studied experimentally the pressure drop and local phase fractions of the stratified and dispersed oil and water. Amundsen, in [7], investigated experimentally the effects of the pipe inclination and water cuts on the flow regime and the pressure drop for a constant mixture velocity. Balakhrisna et al., in [8], studied experimentally the flow patterns for oil-water and kerosene-water of high viscous oil through a sudden contraction and expansion of a horizontal pipe. The contraction and expansion two-phase flow coefficients are lower than the corresponding one for a single-phase water flow. Following the experimental study of Al-Awadi in [9], the water injection direction (horizontally or vertically) through the T-junction pipe does not affect the flow regime and pressure gradients.Numerical results and experimental investigations for a stratified oil-water flow in the horizontal pipe were conducted by Santos et al. in [10]. The considered pipe length in the computation is reduced compared to the experiments. Liu et al., in [11], validated the numerical model based on the volume of fraction (VOF) in the case of the oil-water flow through an elbow system for upward, downward, and horizontal configurations. Shi, in [12], conducted an experimental and numerical study of high-viscosity oil-water flow in horizontal pipes for two diameters. A T-junction cylindrical pipe was considered to predict the (CAF) regime. Due to the higher viscosity of the oil phase, a laminar parabolic profile was specified for the oil inlet at the branch pipe and a turbulent power profile was imposed for the water inlet at the main pipe. The SST k-ω turbulence model with turbulence damping at the interface leads to more accurate results. Al-Jadidi, in [13], conducted a series of numerical (CAF) flow using the large Eddy simulation (LES) model coupled with the VOF mixture model to simulate the heavy oil and water flow regimes. Konings, in [14], imposed the oil and the water inlet velocities at the same cross-section inlet to reproduce the CAF flow under periodic boundaries conditions at the inlet and the outlet. Jiang et al., in [15], studied numerically the CAF regime using a T-junction pipe. Ayuba et al., in [16], studied numerically the CAF flow using lateral and central injections of water and oil, respectively, through a horizontal pipe with a uniform cross section.In the previous numerical studies, the CAF is generated by different ways, but the comparison between the different approaches is not analyzed. The objective of the present study is to investigate the CAF by two methods: through a T-junction pipe (3D-T) and through a straight pipe (3D-S). The principal recommendations corresponding to the accuracy of the CAF regime generation is then deduced.


## 2. Problem Description

Figure1 represents the two geometric configurations usually adopted to generate the (CAF) regime. The T-junction configuration is formed by two circular pipes. The main and the branch pipe having the same diameter as indicated by the experimental conducted by Charles et al. in [5] which is D=0.026m. The length of the pipe upstream the junction is B=0.2m and downstream the junction is L=4m sufficient to obtain the developed flow. The length of the branch pipe is equal to B. In the configuration of a T-junction (3D-T), the oil enters at the inlet of the branch pipe (OI) and the water is injected at the inlet of the main pipe (WI). In the second configuration of a straight pipe (3D-S), the water is injected circumferentially at the same cross section as the core from the oil inlet. The test section is located at z=100andD=2.6m from the pipe inlet as considered by Shi in [12]. In all simulated cases, a uniform gravity field is imposed equal to gy=−9.81m/s2.Figure 1
Sketch of 3D geometries of the T-junction (3D-T) and a straight horizontal pipe (3D-S) for a (CAF) flow regime.In this case, the inner radius pipeR1 is given by the following equation (Konings) in [14]:(1)R1=R2Usw/Usoh+1,where R2 is the outer radius pipe R2=D/2=0.013m, Usw is the superficial velocity of water, Uso is the superficial velocity of oil, and h is the holdup ratio h=1.39.This givesR1=0.0084m for Usw=Uso=0.55m/s which correspond to a mixture velocity Um=Uso+Usw=1.1m/s and an inlet water cut, Cw=0.5.The fluid density and dynamic viscosity of each phase are(i)
For oil phase:ρo=997.2Kg/m3andμo=0.0168Pa.s(ii)
For water phase:ρw=998Kg/m3andμw=0.001Pa.sBased on these fluid characteristics, the corresponding Reynolds numbers for oil and water phases are, respectively,(2)Reso=ρoUsoDμo=848.8,Resw=ρoUswDμw=14271.4.Following these values, the turbulent flow model is needed for the two-phase mixture.


## 3. Mathematical Model

### 3.1. Governing Transport Equations

In the mixture model, the same Reynolds average Navier–Stokes (RANS) equations are used for the two immiscible phases constituted by water and oil. These equations can be written in the following form (Shi et al.) in [17]:The average mass conservation:(3)∂ρ∂t+∂∂xjρUj=0,j=1,2,3.The average momentum conservation equations”(4)∂∂tρUi+∂∂xjρUjUi=−∂P∂xi+∂∂xjμ+μt∂Ui∂xj+gi+Fi,i,j=1,2,3,where ρ is the mixture density and μ is the mixture dynamic molecular viscosity given by(5)ρ=ρ1α1+ρ2α2,μ=μ1α1+μ2α2,where ρo and ρw are, respectively, the density of the oil and water, αo and αw are the volume fraction of each phase such that αo+αw=1, Ui are the mean velocity components, P is mean pressure, μt is the dynamic turbulent viscosity, gi is the gravity acceleration field, and Fi are the components of external force per unit volume associated with the surface tension.Thek-ω model with shear stress transport is adopted as the turbulence model. Following this model, the turbulent viscosity is defined as (Wilcox) in [18]:(6)μt=fμρkω.The turbulent kinetic energyk transport equation:(7)∂∂tρk+∂∂xjρkUj=∂∂xjμ+μtσk∂k∂xj+ρPk−f2ε.The specific dissipation rateω transport equation:(8)∂∂tρω+∂∂xjρωUj=∂∂xjμ+μtσω∂ω∂xj+ρωf1C1ωPkk−C2ωω,where Pk is the volumetric production rate of k:(9)Pk=μtρ∂Ui∂xj+∂Uj∂xi∂Ui∂xj.fμ, f1, and f2 are low Reynolds damping functions.The turbulent dissipation rateε is given by(10)ε=CDωk.The following constants are commonly used:(11)σk=2.0,σω=2,CD=0.09,C1ω=59,C2ω=340.At high Reynolds numbers, the damping functions are set to unity and given by(12)fμ=1/40+Re,t/Rk1+Re,t/Rk,Rk=6.0,f1=1fμ0.1+Re,t/Rω1+Re,t/Rω,Rω=2.7,f2=5/18+Re,t/RB41+Re,t/RB4,RB=8.0,where Re,t is the turbulent Reynolds number:(13)Re,t=ρkμω.At high Reynolds numbers, the damping functionsfμ, f1, and f2 tend to unity.According to Brackbill et al., in [19], the surface tension at the interface between oil and water is transformed to a continuum surface force (CSF). Following this model, the surface tension force is proportional to the curvature κ:(14)F⇀=σκρ∇⇀α11/2ρ1+ρ2,where κ is the interface curvature defined in terms of the divergence of the unit normal as(15)κ=∇⇀.n⇀n⇀.The surface normal is defined as the gradient of the volume fraction of the oil phase:(16)n⇀=∇⇀α1.For the adopted validation case investigated experimentally by Charles et al. in [5], the surface tension between the oil and water phases is σ=0.045N/m. An academic version of a computational fluid dynamics code was used to solve these transport equations and the associated initial and boundaries conditions.


### 3.2. Initial and Boundary Conditions

At the initial time of simulationt=0s, the horizontal pipe with a uniform cross section is filled with water. At the pipe wall, a no-slip boundary condition is considered. At the outlet, zero pressure is specified. For the T-junction pipe, the oil inlet is located at the inlet of the branch pipe and a parabolic laminar profile is specified as follows:(17)Wo,in=2Ws,o1−ry2R22,where ry is the radial distance in the xz-plane: ry=x2+z2.For a single pipe, the oil inlet is located at the same location of the water inlet andry is replaced by rz=x2+y2.For water inlet, a turbulent power profile is imposed:(18)Ww,in=KWs,w1−rzR21/n,where ry is the radial distance in the xy-plane and K is given by the following relation (Schlichting) in [20]:(19)K=n+12n+12n2.The coefficientn depended on the water Reynolds number and is taken equal to n=6.6.The following equations are adopted to calculate the turbulent kinetic energy and the specific dissipation rate at the inlet:(20)k=IU2,ω=εCDk,ε=CD3/4k3/20.1H,where H is a characteristic inlet dimension as the hydraulic radius.


### 3.3. Numerical Methods

For the two inlet configurations (3D-T and 3D-S), the computational domain contains 1100 O-grids in the cross section of the pipe and 400 cells in the axial direction. Hence, a structured mesh was adopted to simulate the CAF of oil and water phases. The first cell layer in the cross section of the pipe is located at 0.2 mm in order to satisfy thaty+<5.The PISO algorithm (pressure implicit with splitting of the operator) is used to solve the velocity-pressure coupling. The same residual10−5 is adopted for all solved variables.As indicated by Shi in [12], the geometric reconstruction (piecewise-linear) approach permits to accurately capture the interface variation between oil and water phases.A fixed time step is chosenΔt=0.01s. For each inlet configuration, the fully developed flow (FDF) is attained if the pressure profile along the pipe axis is unchanged with time.


## 3.1. Governing Transport Equations

In the mixture model, the same Reynolds average Navier–Stokes (RANS) equations are used for the two immiscible phases constituted by water and oil. These equations can be written in the following form (Shi et al.) in [17]:The average mass conservation:(3)∂ρ∂t+∂∂xjρUj=0,j=1,2,3.The average momentum conservation equations”(4)∂∂tρUi+∂∂xjρUjUi=−∂P∂xi+∂∂xjμ+μt∂Ui∂xj+gi+Fi,i,j=1,2,3,where ρ is the mixture density and μ is the mixture dynamic molecular viscosity given by(5)ρ=ρ1α1+ρ2α2,μ=μ1α1+μ2α2,where ρo and ρw are, respectively, the density of the oil and water, αo and αw are the volume fraction of each phase such that αo+αw=1, Ui are the mean velocity components, P is mean pressure, μt is the dynamic turbulent viscosity, gi is the gravity acceleration field, and Fi are the components of external force per unit volume associated with the surface tension.Thek-ω model with shear stress transport is adopted as the turbulence model. Following this model, the turbulent viscosity is defined as (Wilcox) in [18]:(6)μt=fμρkω.The turbulent kinetic energyk transport equation:(7)∂∂tρk+∂∂xjρkUj=∂∂xjμ+μtσk∂k∂xj+ρPk−f2ε.The specific dissipation rateω transport equation:(8)∂∂tρω+∂∂xjρωUj=∂∂xjμ+μtσω∂ω∂xj+ρωf1C1ωPkk−C2ωω,where Pk is the volumetric production rate of k:(9)Pk=μtρ∂Ui∂xj+∂Uj∂xi∂Ui∂xj.fμ, f1, and f2 are low Reynolds damping functions.The turbulent dissipation rateε is given by(10)ε=CDωk.The following constants are commonly used:(11)σk=2.0,σω=2,CD=0.09,C1ω=59,C2ω=340.At high Reynolds numbers, the damping functions are set to unity and given by(12)fμ=1/40+Re,t/Rk1+Re,t/Rk,Rk=6.0,f1=1fμ0.1+Re,t/Rω1+Re,t/Rω,Rω=2.7,f2=5/18+Re,t/RB41+Re,t/RB4,RB=8.0,where Re,t is the turbulent Reynolds number:(13)Re,t=ρkμω.At high Reynolds numbers, the damping functionsfμ, f1, and f2 tend to unity.According to Brackbill et al., in [19], the surface tension at the interface between oil and water is transformed to a continuum surface force (CSF). Following this model, the surface tension force is proportional to the curvature κ:(14)F⇀=σκρ∇⇀α11/2ρ1+ρ2,where κ is the interface curvature defined in terms of the divergence of the unit normal as(15)κ=∇⇀.n⇀n⇀.The surface normal is defined as the gradient of the volume fraction of the oil phase:(16)n⇀=∇⇀α1.For the adopted validation case investigated experimentally by Charles et al. in [5], the surface tension between the oil and water phases is σ=0.045N/m. An academic version of a computational fluid dynamics code was used to solve these transport equations and the associated initial and boundaries conditions.


## 3.2. Initial and Boundary Conditions

At the initial time of simulationt=0s, the horizontal pipe with a uniform cross section is filled with water. At the pipe wall, a no-slip boundary condition is considered. At the outlet, zero pressure is specified. For the T-junction pipe, the oil inlet is located at the inlet of the branch pipe and a parabolic laminar profile is specified as follows:(17)Wo,in=2Ws,o1−ry2R22,where ry is the radial distance in the xz-plane: ry=x2+z2.For a single pipe, the oil inlet is located at the same location of the water inlet andry is replaced by rz=x2+y2.For water inlet, a turbulent power profile is imposed:(18)Ww,in=KWs,w1−rzR21/n,where ry is the radial distance in the xy-plane and K is given by the following relation (Schlichting) in [20]:(19)K=n+12n+12n2.The coefficientn depended on the water Reynolds number and is taken equal to n=6.6.The following equations are adopted to calculate the turbulent kinetic energy and the specific dissipation rate at the inlet:(20)k=IU2,ω=εCDk,ε=CD3/4k3/20.1H,where H is a characteristic inlet dimension as the hydraulic radius.


## 3.3. Numerical Methods

For the two inlet configurations (3D-T and 3D-S), the computational domain contains 1100 O-grids in the cross section of the pipe and 400 cells in the axial direction. Hence, a structured mesh was adopted to simulate the CAF of oil and water phases. The first cell layer in the cross section of the pipe is located at 0.2 mm in order to satisfy thaty+<5.The PISO algorithm (pressure implicit with splitting of the operator) is used to solve the velocity-pressure coupling. The same residual10−5 is adopted for all solved variables.As indicated by Shi in [12], the geometric reconstruction (piecewise-linear) approach permits to accurately capture the interface variation between oil and water phases.A fixed time step is chosenΔt=0.01s. For each inlet configuration, the fully developed flow (FDF) is attained if the pressure profile along the pipe axis is unchanged with time.


## 4. Validation of the Numerical Results

The numerical results are validated for low viscosity results in terms of the mixture velocity profile and the water volume fraction for the oil-water flow. The comparative three-dimensional study between the T-junction pipe and the straight pipe is made for the mixture oil-water velocity equal toUm=1.1m/s and a water cut Cw=0.5.The fully developed flow (FDF) is attained if the pressure profile along the pipe axis is unchanged with time. Figure2 represents the pressure along the pipe axis for different times after the water phase reaches the pipe outlet and for the T-junction inlet configuration denoted 3D-T. All these profiles have the same shape except the approaching outlet region. At the junction of the main and branch pipes, there is a sudden decrease of the pressure, and further with downstream, the pressure decreases gradually. The quasi-steady flow regime is attained from the time t=5s. The predicted pressure drop for this CAF is 498 Pa/m which is close to the experimental value measured by Charles et al. in [5] and equal to 482 Pa/m corresponding to a relative error equal to 3.3%. For the straight pipe (3D-S), the computational pressure profiles show that the FDF is attained at the same time and that the simulated pressure drop is 420 Pa/m and the relative error is around 13.0%.Figure 2
Pressure variation along the pipe axis for different time instants and for the 3D-T configuration.The simulated axial mixture velocity of CAF for the two inlet configurations denoted 3D-T and 3D-S is shown in Figure3. These mixture velocities are compared to the numerical results of Shi at the test line located at x=0andz=100andD=2.6m. The two velocity distributions give a small discrepancy with the Shi profile. The deformation of the axial velocity around the upper pipe wall is well reproduced. This deformation leads to a nonsymmetric character of these profiles and varies with the variation of the two-phase interface.Figure 3
Comparison between the simulated axial velocity for mixture velocity 1.1 m/s and the numerical results of Shi.For the 3D-T configuration, the computed maximum axial velocity isWmax=1.56m/s which is close to the Shi value in [12] equal to 1.52 m/s. The 3D-S slightly overestimates the maximum velocity which is equal to 1.63 m/s.Figure4 shows the simulated oil volume fraction distribution in the same line location x=0andz=100D. The CAF regime is more concentric for the 3D-S. However, the oil volume fraction for the 3D-T approaches the Shi profile. In all simulated cases, the film water film thickness at the upper pipe wall is slightly greater than at the lower pipe wall. Conformed to Shi in [12], the CAF regime is more concentric when the inertial force is dominant.Figure 4
Comparison between the simulated oil volume fraction for mixture velocity 1.1 m/s and the numerical results of Shi.The comparison of the oil volume fraction in theyz-plane for the two inlet configurations is represented in Figure 5. Each inlet configuration induces a different interaction between the oil and water phases. This results in different shapes of the CAF regime. The upper water film is thicker for the 3D-S, and this film is more irregular for 3D-T.Figure 5
Comparison of the oil volume fraction in theyz-plane for 3D-T and 3D-S configurations.These observations are confirmed by representing the oil volume fraction inxy-plane for different locations along the pipe axis, as shown in Figure 6. The CAF for 3D-T is more eccentric around the pipe axis especially in the downstream region of the horizontal pipe located away from z>2.6m.Figure 6
Comparison of the oil volume fraction in thexy-plane for 3D-T and 3D-S configurations.As confirmed by Shi simulations, the low turbulent kinetic energyk is located at the core of the CAF and high k located at the annular water region and close to the wall (Figure 7). The inlet configurations affect the distribution of k due to the eccentricity of the CAF in 3D-T configuration.Figure 7
Comparison of the turbulent kinetic energy in thexy-plane for 3D-T and 3D-S configurations.The cross-section distributions of the turbulent kinetic energy for the CAF regime are represented in Figure8. These distributions match the oil and water phase cross-section variations with a higher value in the water annular.Figure 8
Comparison of the turbulent viscosity in thexy-plane for 3D-T and 3D-S configurations.


## 5. Conclusions and Perspectives

The AF regime of the liquid-liquid two-phase flow is validated for two inlet configurations at the horizontal and circular pipe. Compared to the experimental results of Charles et al. in [5] and numerical simulations of Shi in [12] good agreement was found in terms of pressure drop, oil volume fraction distribution, and axial velocity. To predict the interface deformation due to the interaction of the oil and water phase, the present model is based on the VOF convective transport equation coupled with the RANS equations and k-omega turbulence model.This proposed model will be validated for the CAF regime with the high-viscosity oil-water flow, stratified, and slug flows. The large Eddy simulation (LES) can be used to better reproduce the small turbulence structure for this complex immiscible two-phase flow. As a future numerical investigation, we also propose a more complicated pipe configuration of the liquid-liquid flow such as those through an inclined or curved pipe. Finally, as a perspective of this study, we propose the development of a more sophisticated two-phase transport model taking into account the interface transfer between the liquid-liquid phases.



---

*Source: 5296700-2021-10-31.xml*