Evaluation of drift-velocity closure relationships for highly viscous liquid-air slug flow in horizontal pipes through 3D CFD modelling

Journal Pre-proofs Evaluation of drift-velocity closure relationships for highly viscous liquid-air slug flow in horizontal pipes through 3D CFD modelling Juan Pablo Valdés, Paula Pico, Eduardo Pereyra, Nicolás Ratkovich PII: DOI: Reference:

S0009-2509(20)30069-5 https://doi.org/10.1016/j.ces.2020.115537 CES 115537

To appear in:

Chemical Engineering Science

Received Date: Revised Date: Accepted Date:

23 September 2019 12 January 2020 1 February 2020

Please cite this article as: J. Pablo Valdés, P. Pico, E. Pereyra, N. Ratkovich, Evaluation of drift-velocity closure relationships for highly viscous liquid-air slug flow in horizontal pipes through 3D CFD modelling, Chemical Engineering Science (2020), doi: https://doi.org/10.1016/j.ces.2020.115537

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier Ltd.

Evaluation of drift-velocity closure relationships for highly viscous liquidair slug flow in horizontal pipes through 3D CFD modelling Juan Pablo Valdés1, Paula Pico1*, Eduardo Pereyra2, Nicolás Ratkovich1 1 Department

2 McDougall

of Chemical Engineering, Universidad de los Andes, Carrera 1 # 18a-12 Bogota, Colombia

School of Petroleum Engineering, The University of Tulsa, Tulsa, OK 74104, United States *Corresponding Author: [email protected]

ABSTRACT (Formatted) This study assesses the accuracy and validity of seven common drift velocity correlations for horizontal pipes through a 3D-Computational Fluid Dynamics (CFD) model. For this assessment, 13 experimental measurements from the literature were used to validate the model, and 11 case-studies were proposed in order to expand the range of analysis to high-viscosity fluids (μ>0.3 Pa·s). The results indicated that the CFD model offers an excellent drift velocity prediction with an average absolute relative error of 6%. The correlations with the highest predicting capabilities were Jeyachandra’s (2012) and Livinus’ (2018). These studies considered the combined effect of the Viscosity and Eötvös Numbers. In addition, a new correlation was developed by fitting the CFD results for all case-studies proposed. This new correlation improves the prediction achieved by other authors as the overall average error obtained (18%) is 4% lower than the one found for the correlations of Livinus (2018) and Jeyachandra (2012). ABSTRACT In the past 50 years, several correlations have been developed in order to predict the drift velocity of a liquidair slug-flow system inside horizontal pipelines. However, a great portion of these correlations is based on either very limited experimental data in terms of fluid properties and operating conditions or on often invalid assumptions. Considering this, the present study assesses the accuracy and validity of seven of the most common non-complex drift velocity correlations through the development of a 3D-Computational Fluid Dynamics (CFD) model using the commercial software Star-CCM+. For this assessment, 13 experimental measurements from the literature were used to validate the CFD model, and 11 case-studies were proposed in order to expand the range of analysis to high-viscosity fluids (μ>0.3 Pa·s) often found in industrial settings. The two-phase system was modelled with an Eulerian-Eulerian approach, coupled with the Volume of Fluid (VOF) method. Three VOF parameters were calibrated in order to obtain the best configurations for different viscosity ranges. The results indicated that the CFD model offers an excellent prediction of the drift velocity given than the average absolute relative error obtained was 6%, even considering highly viscous fluids (μ ~6.86 Pa·s). The model allowed to confirm that the slug units and gas bubbles do not behave symmetrically in the axial direction and therefore, a 3-D approach substantially improves the accuracy of the model when compared with some 2-D models developed previously. The correlations that showed the highest predicting capabilities were Jeyachandra’s 2012 correlation and Livinus’ 2018 correlation. These two studies considered the combined effect of the Viscosity and Eötvös Numbers, indicating that both parameters significantly influence the drift velocity. In addition, a new correlation was developed by the fitting of the CFD results for all the case-studies proposed. This new correlation slightly improves the prediction achieved by other authors as the overall average error obtained (18%) is 4% lower than the one found for the correlations of Livinus (2018) and Jeyachandra (2012) and around 11% and 6% lower for the high-viscosity range (Nvis > 0.1) when compared to those two correlations, respectively. Keywords: Drift Velocity, CFD, Closure Relationships, Slug flow, Two-phase, High-Viscous Liquid.

1

NOMENCLATURE R

Buoyancy Reynolds number

AF C

Fitting parameters for the proposed correlation Angle Factor Courant Number


Re

Reynolds Number

SF

Sharpening Factor

CFL

Courant-Friedrichs-Lewy

T

Fluid temperature, ºC

CFD

Computational Fluid Dynamics Distribution parameter or flow Coefficient Pearson correlation coefficient

t

Physical Time, s

∆t

Timestep, s

VOF

Volume of fluid

𝐯

Velocity vector, m·s-1

vd

Drift Velocity, m·s-1

vdi,j ext

Extrapolated solution of grid 𝑖 based on grid 𝑗, m·s-1

a,b,c,e,f

Co Cp d eri,j a eri,j ext

Distance travelled by the bubble, m Relative difference between the solution of grid 𝑖 and 𝑗, % Relative difference between the solution of grid 𝑖 and the extrapolated solution of grid 𝑖 based on grid 𝑗, %

vd, CFD i vd, exp i

Drift velocity of case i predicted by CFD, m·s-1 Experimental drift velocity of case i, m·s-

𝐅

Eötvös Number Eötvös Number defined by each correlation Source term at the interphase, kg·m-2·s-2

Fr

Froude Number


Greek symbols αk

kth-phase volume fraction

GIT

Grid Convergence Index of grid i in terms of grid j, % Grid Independence Test

ρ

Density, kg·m-3

g

Gravity, m·s-2

μ

Dynamic viscosity, Pa·s

h

Average size of the cell, cm

μart

Artificial viscosity

IMD

Interface Momentum Dissipation


Σ

Dimensionless surface tension

ID

Pipeline’s Internal Diameter, m

σ

L

Length, m

N

Total number of Cells

Nvis

p

Viscosity number Viscosity number defined by each correlation Apparent order of discretization

P

Pressure, Pa

Eo Eo

GCIi,j

Nμ

1

Average ∆x

Length interval, m

θ

Surface tension between air and the liquid, N·m-1 Inclination angle of the pipeline, °

χ2 GOF p-value

p-value calculated by the Chi-squared Goodness of fit test

Subscripts L

Liquid

g

Gas

1. INTRODUCTION The presence of multiphase flow, particularly gas-liquid, through pipeline systems is a commonly encountered phenomenon in the operation of chemical, geothermal, nuclear, and Oil & Gas (O&G) industries [1]. However, due to its complex behaviour, the understanding of multiphase flows is minimal in comparison to single-phase flows, especially when dealing with high-viscosity liquid phases. Most of the current multiphase models and correlations are heavily dependent on experimental data gathered for low viscosity 2

fluids (viscosities commonly under 0.02 Pa·s) [2]. Additionally, these models do not consider properly the combined effect of viscosity, surface tension, and operational conditions of the pipeline, which results in low predictive capabilities and a narrow range of applicability [3]. Moreover, extensive use of computational models for multiphase transport on pipelines has been made, but mostly one dimensional (1D) models have been considered [4]. Therefore, there is a great interest to extend the range of applicability of the existing knowledge for medium and high viscosity fluids, given the recent focus of the O&G industry on heavy crude oils as they account for 70% of the available reserves nowadays [5]. There is also increasing interest in developing three-dimensional (3D) Computational Fluid Dynamics (CFD) models for pipeline design, which are benchmarked against representative experimental flow cases. Several design and sizing parameters for multiphase flows in pipelines strongly depend on the spatial distribution of the phases present or the ‘flow pattern’. One of the most common flow patterns observed for gas-liquid flow in pipelines is the slug flow, which describes the motion of a slug unit [1]. The slug unit is exemplified in Figure 1, where its main characteristics, including slug length, bubble length, and the slug parts, are highlighted. One of the critical closure relationships in two-phase slug flow modelling is the translational velocity, which primarily depends on the drift velocity [2]. The drift velocity parameter refers to the velocity at which the gaseous phase penetrates and displaces the stagnant liquid phase within the tube, as depicted in Figure 1.

Figure 1. Illustration of the slug unit taken and adapted from [1]

The drift velocity and motion of large gas bubbles in slug flow in pipelines has been extensively analyzed in the past decades. Several studies have been conducted on vertical or inclined pipes experimentally (Dumitrescu (1943) [6], Davies and Taylor (1950) [7], Zukoski (1966) [8], Bendiksen (1984) [9], Weber et al. (1986) [10], Alves (1993) [11], Jeyachandra (2012) [2], Moreiras et al. (2014) [5]) and computationally (Taha and Cui (2006) [12], Ramdin (2011) [13], Ramdin and Henkes (2012) [4]). However, the scope of this work will center entirely in a horizontal arrangement. Some of these works provide unified correlations for horizonal, vertical, and inclined pipelines, for which they will be included in the present work. The first analytical study to estimate the drift velocity in horizontal pipes was carried out by Benjamin (1986) [14], following a potential flow theory, which neglects the effects of surface tension and viscosity. The drift velocity therefore only included the effects of gravity and pipeline geometry. Benjamin (1968) demonstrated analytically that the drift velocity component is not zero, even for a horizontal arrangement of the pipeline. Other studies, such as Wallis (1969) [15], Dukler & Hubbard (1975) [16] and Heywood & Richardson (1979) [17] initially claimed a null drift velocity in the estimation of the translational velocity for horizontal flow, considering that gravity cannot act in the horizontal direction. However, it was shown experimentally ([18],[9]) that the drift velocity component exists for the horizontal flow and that it can even exceed the vertical drift velocity value. Experimentally, Zukoski (1966) [8] was one of the first to investigate the effects of liquid viscosity, surface tension, and pipe inclination on the behaviour of elongated gas bubbles travelling through a stagnant liquid column with varying pipe diameters. He determined that the effects of liquid viscosity become negligible on the drift velocity component for the elongated gas bubble at Reynolds (Re) numbers larger than 200. On the other hand, it was observed that for Re number lower than 4, the drift velocity becomes inversely proportional 3

to the liquid viscosity. Gardner and Crow (1970) [19] experimentally investigated the motion of a large air bubble displacing stationary water on a horizontal channel and analyzed the influence of surface tension on the bubble’s motion. It was found that, even for a dimensionless surface tension (Σ~10 ―4), the effects are still present and must be considered in any theoretical treatment. Additionally, a hypothesis on the constant value of the radius of curvature of the two-phase interface against channel depth was developed to explain the experimental influence observed for the surface tension. Later on, based on Zukoski’s study, Weber (1981) [18] developed a correlation that only accounts for surface tension forces acting on the nose of the bubble. This correlation was based on a balance of forces between the hydrostatic pressure and the capillary pressure. Weber (1981) found a critical Eötvös number (Eö) value at which no drift velocity on horizontal pipelines could exist. Simultaneously, Wilkinson (1982) [20] carried out experiments and proposed a similar dimensionless model to Weber (1981) based on Benjamin’s 1968 analysis and conservation of hydrostatic and capillary forces between the tail and the tip of the bubble. Bendiksen (1984) [9] studied experimentally the influence of the angle of inclination and pipe diameter, and found that the drift velocity can be expressed as the sum of a horizontal and vertical component, weighted by the cosine and sine of the angle of inclination, respectively. The study of Weber et al. (1986) [10] showed that the drift velocity was very sensitive in the near horizontal region for sufficiently large Eo numbers, which correspond to bigger pipe diameters and/or fluids with smaller surface tension. Later on, Gockal (2008) [21] focused his experimental study on the effect of high-viscosity oils on two-phase slug flow in horizontal and upward inclined pipelines. He concluded that, for decreasing Reynolds numbers, the bubble velocity decreases. Intermittent slug flow and elongated bubbles were the dominant flow pattern generated. These studies led to the formulation of initial correlations that considered the effects of viscosity on the calculation of pressure drop, translational velocity, drift velocity, slug length, and frequency. Jeyachandra (2012) [2] studied the drift velocity experimentally considering different viscosity and pipe diameter ranges to those used in the study of Gockal (2008). A new empirical correlation for the drift velocity on horizontal and inclined pipelines was proposed. Moreiras et al. (2014) [5] proposed a closure relationship for the drift velocity on horizontal flow as a function of the viscosity number Nvis, based on experimental data with large Eo numbers for which the surface tension effects could be neglected. More recently, Livinus (2018) [3] carried out an extensive work analyzing most of the previously-mentioned models and closure relationships for the drift velocity and proposed a generalized correlation based on previous and own experimental data. Although a very well-developed correlation was proposed, the range of applicability in terms of viscosity for the horizontal drift velocity was not tested for highly-viscous fluids over 6 Pa·s. The experimental data gathered and measured for horizontal pipelines commonly considered viscosities under 1 Pa·s. Various computational studies have corroborated and extended the conclusions obtained by previous experimental results. Andreussi and Bonizzi (2009) [22], as well as Ramdin and Henkes (2012) [4] performed simulations on the commercial CFD software FLUENT with a Volume of Fluid (VOF) approach. They found that the drift velocity of the gas bubble in a horizontal configuration not only decreases with decreasing Reynolds number, but also decreases as it travels along the pipeline. This trend was observed particularly for high-viscosity fluids. Andreussi and Bonizzi used a 2D approach and a half-section 3D symmetry plane. Ramdin (2011) [13], using 2D and 3D simulations, found a dependency of the Froude number (Fr) on the Re and Eö number, which accounts for the effect of viscosity and surface tension, respectively. Kroes & Henkes (2014) [23] studied the drift velocity of elongated gas bubbles in liquid-filled horizontal pipelines using FLUENT (ANSYS 14). The authors obtained a reasonably-good agreement with the inviscid analytical solution and the experimental data. It was found that the drift velocity decreases with increasing liquid viscosity and over time as the bubble travels. 2D and half-section 3D approaches were considered. However, the 3D model approximated the half-section with a symmetry plane. Lu (2015) [24] carried out an experimental and CFD investigation on horizontal gas-liquid slug-flow, in order to describe the slug initiation, growth, and collapse along the pipeline. Six CFD codes were employed (STARCD, CFX, FLUENT, LedaFlow, TRIOMPH, TransAT), to compare convergence methods in terms of accuracy and fast processing. All simulations carried out considered 3D half-section pipelines with symmetry 4

planes. Similarly, Sanderse et al. (2015) [25] simulated gas bubbles travelling through a channel using the CFD software FLUENT and considering a 2D approach. The bubble velocity and pressure drop were predicted along the channel and compared against Benjamin’s 1968 analytical solution with good agreement. As can be noted from the state of the art previously presented, the gas bubble dynamics for slug flow have been extensively investigated, both through experimental and CFD approaches. Nevertheless, most of the experimental focus has been given to low and medium viscosities (<1 Pa·s), and most of the CFD studies have mostly considered 2D or half-section 3D symmetric approaches, given their simple and inexpensive yet accurate results for some flow cases. However, it has been demonstrated in previous studies that the flow distribution for two-phase systems can either be symmetric or asymmetric depending on the flow conditions, fluid properties, and pipe inclination [26]. Therefore, a more robust 3D CFD model, which models appropriately high viscosity two-phase flows and fits with previous analytical/empirical correlations and relationships, is yet to be done. Moreover, an evaluation of the existing correlations for a more comprehensive set of flow cases with high viscosities and different flow conditions has not been carried out. Although the study of Livinus (2018) [3] proposed a simplified correlation based on previous closure relationships, high viscosity cases associated to mature oil wells were not considered. This analysis could give valuable insights into the applicability of the current closure relationships for heavy crude oils, which correspond to the majority of the reserves available nowadays [5]. Accordingly, the main objective of this work is to generate a 3D CFD model that replicates and extrapolates flow cases for very high viscosities and different pipe diameters and to compare its results against correlations and data presented in previous works. This 3D CFD model is intended to capture non-symmetric flow phenomena of the slug units and to predict flow cases for a wide range of Nvis and Eo numbers, in order to assess the applicability of the current closure relationships and correlations proposed in the literature. To achieve this, a Eulerian multiphase approach with a VOF model was be used. A sensitivity analysis of various VOF model configurations was carried out to further validate the CFD model constructed. Finally, based on the validating and extrapolated cases calculated through CFD, a new correlation for horizontal drift velocity will be proposed with the purpose of considering the broader range of scenarios studied in this work.

2. MATERIALS AND METHODS This section contains a detailed description of the case-studies considered for the analysis of the drift velocity and the existing closure relationships. This section also details the 3D-CFD model constructed, including the spatial and temporal discretization methods used, the physical models selected and the post-processing of the results.

2.1. Case-Studies and drift velocity correlations evaluated As mentioned in the introduction section, this study aims at analyzing the validity range and assumptions of specific drift velocity correlations proposed by other authors using CFD. Considering this, the fluid properties and operating conditions of several experimental works were taken as validation sources for the CFD model developed, as well as a benchmark for extrapolated cases proposed in this study. Table 1 shows the density, viscosity, surface tension, and operating conditions of the 24 cases considered (13 real cases from other authors and 11 extrapolated cases).

5

Table 1. Case-studies contemplated in this study. All cases considered air at atmospheric conditions as the gaseous phase

Case

Fluid

ρL kg ( m3)

μL (Pa·s)

σ (N/m)

𝑇𝐿 (°𝐶)

ID (mm)

d (m)

Re

Nvis

Eo

Source

1

Oil 1

856.704

0.039

0.028

120.000

50.800

0.500

324.154

0.001

774.585

[5]

2

Oil 2

884.000

0.342

0.029

28.100

50.800

0.500

20.747

0.011

771.704

[21]

3

Oil 3

882.000

0.298

0.029

30.000

50.800

0.500

26.613

0.009

769.958

[21]

4

Oil 4

876.000

0.165

0.029

39.400

50.800

0.500

72.280

0.005

764.720

[21]

5

Oil 5

878.000

0.201

0.029

36.100

50.800

0.500

46.600

0.006

766.466

[21]

6

Oil 6

881.000

0.248

0.029

32.800

50.800

0.500

35.190

0.008

[21]

7

Oil 7

889.000

0.692

0.029

19.200

50.800

0.500

6.918

0.022

769.085 776.069

8

Oil 7*

889.000

0.692

0.029

19.200

44.000

0.500

-

0.027

582.208

9

Oil 7*

889.000

0.692

0.029

19.200

37.300

0.500

-

0.035

418.399

10

Oil 7*

889.000

0.692

0.029

19.200

24.000

0.500

-

0.067

173.219

11

Oil 8

888.000

0.593

0.029

21.000

50.800

0.500

9.357

0.019

775.196

[21] [21] (E.C.) [21] (E.C.) [21] (E.C.) [21]

12

Oil 9

876.000

0.150

0.029

40.000

152.400

0.914

453.908

0.001

6882.481

[2]

13

Oil 10 Corn Syrup 95%* Corn Syrup 95%*

888.000

0.574

0.029

-

152.400

0.914

103.738

0.003

6976.762

1340.000

1.830

0.081

29.000

50.800

0.600

-

0.038

418.809

1340.000

1.830

0.081

29.000

44.000

0.600

-

0.047

314.191

[2] [10] (E.C.) [10] (E.C.)

14 15 16

Corn Syrup 95%

1340.000

1.830

0.081

29.000

37.300

0.600

1.006

0.061

225.791

[10]

17

Corn Syrup 95%*

1340.000

1.830

0.081

29.000

24.000

0.600

-

0.117

93.4784

[10] (E.C.)

18

Corn Syrup

1410.000

6.120

0.087

29.000

37.300

0.600

0.085

0.192

221.201

19

Honey*

1468.687

6.857

0.054

25.000

50.800

0.500

-

0.130

691.054

20

Honey

1468.687

6.857

0.054

25.000

44.000

0.500

0.198

0.162

518.430

21

Honey*

1468.687

6.857

0.054

25.000

37.300

0.500

-

0.207

372.565

22

Honey*

1468.687

6.857

0.054

25.000

24.000

0.500

-

0.401

154.244

23

Athabasca Bitumen*

1044.888

57.000

0.030

40.000

44.000

0.500

-

1.887

661.489

24

Athabasca Bitumen*

1044.888

168.000

0.033

31.000

44.000

0.500

-

5.562

601.354

*These additional extrapolated cases (E.C.) were generated using the properties of the fluid of each specific reference and varying the diameter of the pipeline.

The cases highlighted with the * symbol correspond to those extrapolated cases that have been included in this study to expand the analysis range of Nvis and Eo numbers, but that have not been tested experimentally by any of the studies mentioned. Therefore, as will be seen in the Results and Discussion section, only the cases not highlighted were used as validation for the 3D-CFD model. From Table 1 it can be noted that the cases selected cover a wide range of Nvis numbers (1E-2 to ~10) but are mainly focused on high-viscosity 6

[10] [27] (E.C.) [27] [27] (E.C.) [27] (E.C.) [28] (E.C.) [28] (E.C.)

fluids (high Nvis numbers) as well as on various Eo numbers. Most closure relationships that have been proposed up to date only include the influence of one of these dimensionless numbers or limit their applicability to very narrow ranges. Considering this, 24 particular cases were selected to enable a more thorough and comprehensive understanding of the behaviour of the drift velocity under different operating conditions. It is relevant to mention that recent studies, such as Andreussi (2009) [22] and Pico et al. (2018) [27], have demonstrated that the drift velocity exhibits an inverse relationship with the distance travelled by the bubble, and that this relationship becomes more prominent as the viscosity of the liquid phase increases. Hence, the drift velocity cannot be considered a constant parameter during the development of an individual test. As can be noticed in Table 1, the drift velocity of each case-study is taken at a fixed distance travelled by the bubble (d), which implies that all conclusions drawn from this study are only applicable for fixed distances. A total of seven of the most commonly known correlations developed for the prediction of the drift velocity considering two-phase flow in horizontal, circular pipelines were evaluated in this study. A summary of the mathematical expression of each of these correlations in terms of specific dimensionless numbers is shown in Table 2. Table 2. Horizontal drift velocity correlations evaluated in this study

Correlation

Reference

vd = 0.542 gD

Benjamin (1968) [14]

vd = (0.54 ― 1.76Eo―0.56) gD Where, Eo = gD2(ρL ― ρg)σ ―1 vd = (0.5 ― 1.36Eo―0.5) gD Where, Eo = ρLgD2σ ―1

Weber (1981) [18]

Wilkinson (1982) [20]

vd = vhdcosθ + vvdsinθ Where, vhd = 0.542 gD

vvd = 0.351 gD

0.46

vd = (0.53e ―13.7Nμ Where, Nμ = μρL―1(gD3)

[

vd = 0.54 ―

Eo―0.1

)

―0.5

]

Nμ

Bendiksen (1984) [9]

gD

Eo = gD2ρLσ ―1

Jeyachandra (2012) [2]

gD(ρL ― ρg))0.5

―0.5(

ρ 0.01443Nμ + 1.886 L

Moreiras (2014) [5]

―0.5 Where, Nvis = μ(gD ρL ― ρg)ρL) 3(

7

FrH = 10 ―m

m = ( ―0.02861x3) + (0.5987𝑥2) + ( ―4.139𝑥) + 9.843 Where, 𝑥 = 𝐿𝑜𝑔10(𝑅 ∙ 𝐸𝑜) 0.5

R = (D3g(ρL ― ρg)ρL) μ ―1

Livinus (2018) [3] [29]

Eo = gD2(ρL ― ρg)σ ―1

The correlations evaluated vary significantly from each other in terms of their considerations regarding the effects of different parameters (such as the surface tension and viscosity), their applicability range, and their assumptions. As can be seen in Table 2, the theoretical correlation developed by Benjamin (1968) [14] predicts a horizontal drift velocity which depends only on the diameter of the pipeline, ignoring the effects of viscosity and surface tension This occurs similarly to the correlation developed by Bendiksen (1984) [9] for horizontal pipes. The semi-empirical correlations of Weber (1981) [18] and Wilkinson (1982) [20] include the effects of the surface tension through the Eo Number at low Fr numbers but do not consider the effects of the fluid’s viscosity. On the contrary, the study of Moreiras et al. (2014) [5], based on the observations made by Zukoski (1966) [8], neglects the effects of the surface tension given that, in most practical applications, the Eo overcomes the critical value and does not influence the drift velocity. The correlations of Jeyachandra (2012) [2] and Livinus (2018) [3] [29] include the influence of both, surface tension and viscosity on the drift velocity at different diameters. It is essential to highlight that the correlation published by Livinus (2018) [3] had a miscalculation on the expression formulated, for which the revised and corrected model was the one considered and analyzed in this study [29].

2.2. Numerical Simulations The simulations were performed using the commercial software STAR-CCM+ v.13.04, where the main steps of numerical modelling (geometry generation/preparation, spatial discretization, physical models’ selection, and solving of the conservation equations) were made. These CFD simulations provided predictions of the drift velocity at different operating conditions in order to evaluate the existing closure relationships. In addition, they gave valuable insights in regards to the behaviour of the two-phase system itself that would otherwise require sophisticated experimental techniques to obtain. 2.2.1. Pre-Processing 2.2.1.1. Specification of the physical models The physical models used to emulate the behaviour of the system were selected considering the robustness of the available options, as well as the computational requirements and limitations of this study. Firstly, as can be seen in Table 1 for the cases used for validation, even the fluid with the lowest viscosity (Oil 1) would have a drift velocity low enough to consider the flow as fully laminar (Re<500). Therefore, the Navier-Stokes equations were solved without the need for turbulence models. Secondly, it was assumed that the system is fully isothermal, and no chemical reactions occur between the phases, for which the energy equation and source terms related to chemical reactions in the mass conservation equation were neglected. Thirdly, as a result of the low velocities that both phases exhibit in this application, it can be stated that the Mach number will unquestionably be much lower than 0.3 in all the case-studies, and therefore, the flow is incompressible [30]. The Eulerian-Eulerian approach, coupled with the VOF model, was used to simulate the two-phase system, as well as the interaction between the phases. The use of the VOF model, which was first proposed and introduced by Hirt and Nichols [31], has previously been successful in various applications of the O&G industry, such as the modelling of an Electrical Submersible Pump [32], prediction of two-phase slug flow pattern in pipelines [33], cyclones [34], oil spilling [35], among others. This method assumes that all phases present within the system share the same velocity and pressure fields. Therefore, all fluid properties at each cell are averaged by the volume fraction of the phases, resulting in only one set of conservation equations to be solved. The interface is then tracked through the use of an indicator function, which in this case

8

corresponds to the volume fraction [36]. A summary of the essential features of the formulation of this method is shown in Eq.(1)-(5) [37]. Sum of the volume fraction of all the phases within a cell:

∑α = 1;

(1)

k

k

Indicator function (volume fraction):

{

1→𝑇ℎ𝑒 𝑐𝑒𝑙𝑙 𝑖𝑠 𝑐𝑜𝑛𝑡𝑎𝑖𝑛𝑒𝑑 𝑤𝑖𝑡ℎ𝑖𝑛 𝑡ℎ𝑒 𝑔𝑎𝑠 ― 𝑝ℎ𝑎𝑠𝑒 𝑟𝑒𝑔𝑖𝑜𝑛 α𝑔 = 0→𝑇ℎ𝑒 𝑐𝑒𝑙𝑙 𝑖𝑠 𝑐𝑜𝑛𝑡𝑎𝑖𝑛𝑒𝑑 𝑤𝑖𝑡ℎ𝑖𝑛 𝑡ℎ𝑒 𝑙𝑖𝑞𝑢𝑖𝑑 ― 𝑝ℎ𝑎𝑠𝑒 𝑟𝑒𝑔𝑖𝑜𝑛 0 < α𝑔 < 0 →𝑇ℎ𝑒 𝑐𝑒𝑙𝑙 𝑖𝑠 𝑐𝑜𝑛𝑡𝑎𝑖𝑛𝑒𝑑 𝑤𝑖𝑡ℎ𝑖𝑛 𝑡ℎ𝑒 𝑖𝑛𝑡𝑒𝑟𝑝ℎ𝑎𝑠𝑒

(2)

Conservation equations with volume-averaged properties: ∂ρ𝐯

(3)

+ ∇.(ρ 𝐯 𝐯 ) = ―∇P + ∇.[μ(∇𝐯 + ∇𝐯T)] + ρ 𝐠 + 𝐅 ∂t ∂ρ + ∇.(ρ 𝐯 ) = 0 ∂t Advection equation for the volume fraction: ∂αk + v∇ αk = 0 ∂t

(4)

(5)

In this case, the source term in Eq.(3) was taken as the surface tension force between the gas and liquid phases. Eq.(5) represents the evolution of the volume fraction of the secondary phase (air in this case). The discretization of this particular equation is of critical importance in the implementation of the VOF method, given that the sharpness of the interphase depends directly on the discretization scheme used [36]. In this study, a second-order upwind method, combined with a High-Resolution Interface Capturing (HRIC) model, were used as discretization schemes. The upwind scheme can provide a robust and stable solution; however, it can also introduce high numerical diffusion and smear the interphase, which are problems that are corrected by the less robust HRIC model [38]. Previous CFD studies have found that the results obtained with the standard settings of the VOF method implemented in STAR-CCM+ tend to be significantly less accurate when handling liquid phases with viscosities over a value of ~0.3 Pa·s [27]. Based on the VOF calibration made in the study by Pineda et al. [33], three different VOF parameters related to the HRIC discretization scheme or the phase interaction were tested in order to obtain the best configuration for the simulations. These parameters were the Angle factor (AF), the Sharpening factor (SF), and the use of the Interphase Momentum Dissipation (IMD) model through the specification of an artificial viscosity. A total of six configurations for these VOF parameters were tested with two representative fluids in terms of viscosity (Oil 2 and Oil 4, cases 2 and 4, respectively). These fluids were chosen in order to represent the limit values for two ranges of viscosities (>0.3 Pa·s and <0.3 Pa·s). The time estimated on the CFD simulations was compared against the experimental data reported in the study of Moreiras et al. [21]. Table 3 shows the configurations established. The results obtained from these tests and the selected VOF configuration for the rest of the case studies are shown in section 3.1. Table 3. VOF configurations tested. The value of 0 for the IMD indicates that the model is deactivated.

Configuration

SF

AF

μart

1 2 3 4 5 6

0 1 1 0 1 0

0.05 0.2 0.05 0.2 0.05 0.2

0 10 10 0 0 10 9

In applications where the surface tension forces are significant, numerical diffusion problems tend to arise, causing the interphase to be smeared. The SF includes an additional term to the VOF transport equation, reducing the numerical diffusion and allowing sharp interphase to be obtained [39]. Higher values of the AF are mostly applied in simulations where the downwind discretization scheme is used along with the HRIC model, and the interphase is parallel to the flow. The reason being that the AF adds a correction to the wrinkled interphases that are obtained as a consequence of this configuration, forcing the interphase to align with the grid [1]. In the system of interest of this study, the interphase is indeed parallel to the flow; however, an upwind discretization scheme is used and, therefore, and as will be seen further on, this parameter will not have a significant influence on the results. The Continuum Surface Force (CSF) model was used to calculate the surface tension forces between the phases. This model tends to generate parasitic currents in the regions close to the interphase as a result of parameter discontinuities, numerical approximation flaws, and discretization errors. The IMD model adds an additional momentum dissipation term to the momentum equation at the regions close to the interphase in order reduce these parasitic currents that cannot be diminished by finer grid parameters [40]. 2.2.1.2. Boundary and initial conditions Three main boundaries were set to the fluid region of the geometry: the air inlet/liquid outlet, the lateral walls of the pipeline, and the wall located on the opposite side of the inlet/outlet. The two latter boundaries were modelled as no-slip walls with a 0º contact angle between the fluid and the boundary, whereas the air inlet/liquid outlet was modelled as a Pressure Outlet with an air volume fraction of 1 (free flow inlet and outlet). As for the initial conditions, the volume fraction of the liquid phase was set to 1 in the fluid region and the velocity of the system was set to 0 m/s in order to obtain an initially-stagnant system. This combination of boundary and initial conditions ensures the correct behaviour of the flow at the inlet/outlet boundary, i.e., the gas phase enters the pipeline while the liquid phase initially contained within the pipeline exits the system (refer to Figure 2).

(a)

(b)

10

Figure 2. Geometrical domain simulated. (a) Boundary regions of the geometry. (b) Axial view of the velocity vectors of each phase at the inlet/outlet boundary.

2.2.1.3. Spatial discretization, grid independence test, and grid convergence index The spatial resolution of the conservation equations was made through the finite-volumes method. Considering this, the geometry was discretized using an unstructured, non-conformal grid that included polyhedral-shaped cells in the core regions of the geometry, as well as five prism layers in the zones near the lateral and posterior walls. While previous studies of flow through pipelines have mostly used some variation of orthogonal or tetrahedral grids [33], [41], [42], some other authors have recommended the use of polyhedral cells given that the number of neighboring cells tends to be higher, resulting in a better approximation of the gradients, lower skewness angles between cells (less sensitive to cell stretching) and lower numerical diffusion [27], [43]. Additionally, these studies have found that, in general, a lower number of polyhedral volumes is needed to obtain that same level of accuracy when compared to tetrahedral volumes. It is relevant to highlight some features that were added to the grid, which aided to the overall quality and the convergence of the solution. For instance, an additional refinement was placed at the zones reaching the prism layers (Figure 3(a)) in order to correctly capture the high-velocity gradients that arise in the zones near the walls. Additionally, the maximum growth rate from cell to cell was set to 1.3 to obtain an average change in cell spacing that does not overcome 30% in most of the geometry.

(a)

(b)

Figure 3. Grid constructed for pipeline geometry. (a) Frontal view with radial refinement. (b) Lateral cross-section view

The average size of the polyhedral cells in the core regions of the geometry (h) was determined through a Grid Independence Test (GIT) and certain accuracy considerations. The GIT performed used one of the casestudies found in the work of Pico et al. [27] as a source for comparison. This case used hydraulic oil as the operating fluid and a pipeline of 44 mm, where the drift velocity was measured when the air bubble had travelled a distance of 0.5 m. The characteristics of the four grids evaluated can be found in Table 4. These cell sizes were specified in order to meet the spacing criteria between grids that is recommended in the literature: ℎ𝑓𝑖𝑛𝑒 ℎ𝑐𝑜𝑎𝑟𝑠𝑒 > 1.3 [44]. Figure 4 shows the drift velocity values obtained with each grid size tested. Table 4. Grids evaluated in the GIT

Grid name

Grid number

h (cm)

Number of cells

Ultrafine

0

1.0

3’985,728

Fine

1

1.4

1’005,014

Normal

2

2.0

440,803

Coarse

3

3.5

262,022

11

Figure 4. Result of the application of the GIT for hydraulic oil, according to the testing conditions and fluid properties used by Pico et al. (2018) [27]

As portrayed in Figure 4, the number of cells (or average cell size) tends to influence the predicted drift velocity up to the ‘fine’ grid. From this grid onwards, additional refinement does not seem to alter the solution significantly (difference between the vd obtained with the ‘fine’ and ‘ultrafine’ grids of 0.08%, as opposed to the 0.28% difference found between the ‘normal’ and ‘fine’ grids). Accordingly, it can be stated that the solution is approximately independent from the spatial discretization from ~1 million cells onwards, and therefore, the configuration of the ‘fine’ grid was selected to perform all simulations. The uncertainty due to the spatial discretization of the selected grid was estimated through the Grid Convergence Index (GCI) methodology proposed by Celik et al. (2008) [44]. This methodology was applied to three different cases: hydraulic oil (μL = 0.092 Pa ∙ s) and mineral oil (μL = 0.034 Pa ∙ s ) according to the operating conditions of Pico et al. [27], and Case 2 of the present study (μL = 0.342 Pa ∙ s, refer to Table 1). The ‘ultrafine’ grid was not included in the calculations of this method given that it is assumed that, from the ‘fine’ grid and onwards, the solution is grid-independent. Eq.(6)-(9) show the mathematical expression of the most relevant parameters used to estimate the GCI. These include the apparent order of discretization (p), and the drift velocity that would be obtained in a completely converged solution (vdext or extrapolated solution). The subscripts i and j represent the number assigned to each grid evaluated (1 for ‘fine’, 2 for ‘normal’, and 3 for ‘coarse’). ϵgri,j = vdi ― vdj

(6)

rgri,j = hi hj

(7)

p=

|

1

rgr2,1

vdi,j ext

=

|

ln|ϵgr3,2 ϵgr2,1| + ln

rpgri,j vdj ― vdi rpgri,j ― 1

|

(

rpgr2,1 ― sgn(ϵgr3,2 ϵgr2,1) rpgr3,2 ― sgn(ϵgr3,2 ϵgr2,1)

)|

(8)

(9)

From the calculation of p and vdi,j ext, the quality of the discretization can be evaluated by various parameters i,j (Eq.(10)-(12)). The erext parameter quantifies the deviation of the solution obtained by each grid with respect to the extrapolated (fully-converged) solution, the eri,j a is percentage difference between the solution of the different grids, and the grid convergence index (GCIi,j) is a measure of the contribution of the former two parameters and can be understood as the error introduced to the solution by the spatial discretization [44]. 12

Figure 5 contains a comparison between the extrapolated solution and the solution predicted by each grid for each fluid, as well as the uncertainty to the solution associated to the selected grid for the three fluids. Table 5 contains the detailed results of the GCI method. eri,j a =

|

eri,j ext =

i,j

GCI =

|

vdj ― vdi vdj

|

vdi,j ext ― vdj vdi,j ext

(10)

|

(11)

1.25eri,j a

(12)

rpgri,j ― 1

(a)

(b)

Figure 5. Results of the GCI method for the three grids tested and the three fluids. (a) Comparison between the predicted value by the simulations and the extrapolated solution. (b) Error bars associated to the solution obtained with the ‘fine’ grid for each fluid Table 5. Details of the results obtained from the application of the GCI method for the selected grid

Parameter

Mineral Oil

Hydraulic Oil

Oil 2

ϵgr2,1 (𝑚 s)

0.0006

0.0007

-0.0004

ϵgr3,2 (𝑚 s)

-0.0019

0.0008

-0.0023

p (𝑚 s)

2.1683

1.4658

2.8711

0.2785

0.2530

0.1753

vd2,1 ext

𝑒𝑟2,1 𝑎

(%)

0.2237

0.2802

0.2096

𝑒𝑟2,1 𝑒𝑥𝑡

(%)

0.1921

0.4097

0.1173

2,1

𝐺𝐶𝐼 (%) Uncertainty introduced by the ‘fine’ mesh (m/s)

0.2396 6.5857 ∙ 10

0.5100 ―4

1.2954 ∙ 10

0.1469 ―3

2.5710 ∙ 10 ―4

Figure 5(a) indicates that, as expected, an increase in the number of cells brings the solution closer to the extrapolated solution. As mentioned previously, the formal discretization scheme selected for the convective terms was second-order upwind. This order is in reasonable agreement with the apparent order estimated through the GCI method (p) shown in Table 5 for the three testing fluids, which is a general indication that 13

the estimation of the extrapolated solution is satisfactory. From Table 5 it can also be noticed that oscillatory convergence (𝜖𝑔𝑟3,2 𝜖𝑔𝑟2,1 < 1) is only observed for mineral oil [45], [46], but not for higher-viscosity fluids. The GCI values of the ‘fine’ grid indicate that the spatial discretization uncertainties of the simulations are expected to be of the order of 0.1-0.5%, corresponding to ± 2.5 ∙ 10 ―4 to ± 1.3 ∙ 10 ―4 m/s (plotted in Figure 5(b)). As will be seen in the results section, the magnitudes of these uncertainties represent a very low portion of the total error obtained by the simulations. 2.2.2. Numerical Solving As a result of the transient nature of this application, it is of significant importance to select an adequate time discretization scheme. A fully implicit first-order time discretization scheme was selected for all the simulations given that it is unconditionally bounded and, therefore, does not have an upper limit for the appropriate time-step [30]. The value for this parameter was determined using the Courant–Friedrichs–Lewy (CFL) condition (Eq.(13)), which recommends the use of a C ≤ 1. However, it is also recommended the use of a Courant number (C) below 0.1 for the regions near the interphase [33], [47]. Considering this, the timestep was calculated aiming for a 𝐶 ≤ 0.1 and using the minimum cell size encountered as well as the drift velocity values found in the studies that were used as a benchmark (refer to section 2.1). For most of the cases analyzed, the average C number over the entire geometrical domain did not exceed the value of 0.01, and the maximum C number was always below 0.6. The number of iterations per time-step was set as a variable value depending on a threshold set for the residuals of some fundamental transport equations at each iteration. These equations were the continuity, x-momentum, y-momentum, z-momentum, and the air mass balance equations and the threshold values needed to proceed to the following time-step were 1 ∙ 10 ―12, 1 ∙ 10 ―9, 1 ∙ 10 ―8, 1 ∙ 10 ―8, and 1 ∙ 10 ―12, respectively. The number of iterations per time-step generally oscillated between 15 and 20. The convective terms in the momentum equations were discretized using a second-order scheme. C=

U∞Δt

(13)

Δx

2.2.3. Post-Processing The stopping criterion established for the simulations was set with a cross-section plane located at a specific measuring distance, depending on the case study. The air bubble tip was assumed to have arrived at the set plane when the sum of the air volumetric fraction over all the cells contained within the plane reached the threshold value of 0.001. The drift velocity obtained through CFD was calculated as the average velocity of the bubble travelling the specific distance given that this is the methodology followed by most of the experimental studies used in this work to validate the CFD model (see Section 2.1).

3. RESULTS AND DISCUSSION The following section is divided into three main parts. The first one contains the selection of the most suitable numerical configuration for the CFD model, focused on several VOF settings. In this section, a sensitivity analysis of three main parameters is carried out for two flow cases taken from the literature. Secondly, a comparison of the CFD model against experimental data from several flow cases found in the literature, as summarized in Table 1, is presented to validate the further applicability of the numerical model. Moreover, in this subsection, the captured 3D asymmetric characteristic of the slug flow is shown for some representative cases. Finally, extrapolated simulation results are shown and used to validate the fittingness of various closure relationships previously proposed to estimate the drift velocity on horizontal pipes for two-phase flow. 3.1. Sensitivity Analysis for the VOF Model Configuration As previously discussed in subsection 2.2, a calibration of the VOF settings was performed for case-studies 2 and 4, adopting as reference the time taken for the air bubble to reach a given set point, depending on the case considered. The time estimated on the CFD simulations was compared against the experimental data reported 14

in the literature. Following the configurations given in Table 3, the results obtained for the calibration are presented as follows in Table 6. Table 6. Summary of the calibration results obtained for the different VOF configurations tested.

Configuration

Oil 2 time (s)

Oil 4 time (s)

Oil 2 error (%)

Oil 4 error (%)

1 2 3 4 5 6

2.3390 2.8595 2.8560 2.3390 2.3290 2.8660

1.9055 2.2650 2.2650 1.9045 1.9015 2.2675

26.09 9.64 9.75 26.09 26.49 9.43

2.13 21.40 21.40 2.08 1.92 21.54

Before selecting the most suitable configuration for the VOF settings, the effect of each parameter on the elongated bubble front was studied, as shown in Figure 6 and Figure 7. The AF was not analyzed separately, given that the profiles and results obtained between configurations 1 and 4 were nearly the same, as noticed from Table 6. Figure 6 considers solely the effect of the SF, leaving the AF and IMD in their default settings (configurations 1 and 5 of Table 6). As seen from this figure, the effect of the SF is not very notorious on the penetrating bubble and its velocity field, which results consistent with the very similar times calculated for configurations 1 and 5. Nonetheless, a slight decrease in the swirling patterns close to the interface can be noticed, mainly for Oil 4. This is achieved thanks to the anti-diffusion effect, which sharpens the interface and prevents a high dilution of the gas phase onto the liquid phase [39]. Overall, this effect slightly speeds-up the velocity of the bubble, which benefits low viscosity fluids, as it reduces unphysical currents (swirls) generated at the interphase due to numerical diffusion. However, for higher viscosity fluids, it has the opposite effect, as it diminishes the viscous effect of the liquid phase on the penetrating bubble.

Oil 2

SF = 0

Oil 4

SF = 0

Oil 2

SF = 1

Oil 4

SF = 1

(a)

(b)

Figure 6. Elongated bubble tip and velocity vectors for (a) Oil 2 and (b) Oil 4 considering the high and low level of the SF. The top bubbles correspond to configuration 1 and the bottom ones to configuration 5.

On the other hand, Figure 7 considers the effect of the IMD parameter combined with the SF on the bubble velocity field. It can be perceived straight-away that the inclusion of the IMD has a significant effect on the velocity vectors generated at the front of the bubble. This response is also clearly observed in Table 6, where the inclusion of this parameter drastically alters the experimental error obtained for both fluids (configuration 1 vs. 2, 3, and 6). When the IMD is deactivated, high-velocity vectors at the tip of the bubble can be detected, especially for Oil 2. In contrast, when the IMD is activated, most of these vectors are substantially reduced or completely removed. These high-velocity vectors correspond to parasitic currents generated at the interphase due to dominant surface tension forces, high curvatures, and a slow movement of the gas bubble [39]. The 15

bubble slow-motion pushes the system towards an equilibrium between the pressure and surface tension forces. However, the numerical approximations for this conservation of forces on the VOF model do not guarantee this equilibrium due to non-linearities in the discretization [39]. The IMD adds extra dissipation near the interphase (more on the gas phase) in the form of an augmented viscosity in order to reduce these unphysical currents and mitigate this numerical flaw.

Oil 2

SF = 0 IMD OFF

Oil 4

SF = 0 IMD OFF

Oil 2

SF = 1 IMD ON

Oil 4

SF = 1 IMD ON

(a)

(b)

Figure 7. Elongated bubble tip and velocity vectors for (a) Oil 2 and (b) Oil 4 considering low and high levels of SF and IMD (configurations 1 and 3).

As noticed from Table 6, the implementation of the IMD improves the accuracy of Oil 2 (high viscosity range) but deteriorates the results obtained for Oil 4 (lower viscosity range). This observation may be related to the physical properties of both oils and their influence on the penetration of the gas bubble. Higher viscosities of the liquid phase will imply a lower bubble velocity, which will induce more parasitic currents near the tip of the bubble, as already explained. On the other hand, as the liquid viscosity drops the bubble’s velocity rises, diminishing the impact of the flaws associated with the VOF numerical approximations and hence, reducing the formation of these unphysical currents. For this reason, the implementation of the IMD on low-viscosity fluids will add an extra artificial dissipation, generating a viscous effect that will slow down the bubble’s penetration in the liquid phase. In contrast, when the liquid phase exhibits higher viscosities, the artificial dissipation added will mostly remove parasitical currents and will not intensify the liquid’s viscous effect on the gas bubble. These effects can also be appreciated in the bubble’s shape, especially for Oil 4, as it loses part of its curvature when the IMD is activated. As observed in previous studies [33], the IMD tends to represent a smoother interface due to the reduction of the parasitic currents. Consequently, in low viscosity cases, the IMD can modify the hydrodynamics of the bubble by altering its shape, as observed for Oil 4 in Figure 7, which results in a significant modelling error and accounts for the deviation obtained in Table 6. Based on all the parameters analyzed above, a configuration was selected for each of the viscosity ranges considered. For lower viscosity fluids, represented by Oil 4, configuration 5 was selected, with the IMD deactivated to avoid the alteration of the bubble’s hydrodynamics. The SF selected at its highest setting to reduce even further the numerical dissipation, which may slow down the bubble’s motion. For the higher viscosities, configuration 6 was selected, with the IMD activated to reduce the unphysical currents accelerating the bubble’s velocity and with the SF deactivated to avoid an unnecessary reduction of the viscous effect near the interphase. The AF did not have a significant influence on the solution obtained from the CFD model, as seen in Table 6 and, therefore, was not considered as a decisive factor in the calibration performed. 3.2. Validation of the CFD model The selected VOF configurations were applied to all simulations depending on their liquid viscosity, and the accuracy/reliability of the general model was verified by comparing the results of the 13 validating cases 16

simulated against the experimental results taken from previous studies (refer to Section 2.1). Figure 8 shows the results of this comparison. The percentage deviation from the experimental drift velocity of each case was calculated with Eq.(14), and the average absolute relative error was calculated with Eq.(15). vd, error i =

(

vd, CFD i ― vd,exp i

[

vd,exp i

1 % Error = N

N

∑

)

∗ 100%

]

‖vd, CFD i ― vd,exp i‖ vd,exp i

i=1

(14)

∗ 100%

(15)

0.6

CFD drift velocity (m/s)

0.5 0.4

CFD drift velocity (m/s)

Data points X=Y -10% 10%

0.3 0.2 0.1

0.14

Data points X=Y

0.12

-10% 10%

0.1 0.08 0.06 0.04 0.02 0

0 0

0.1

0.2

0.3

0.4

Experimental drift velocity (m/s)

0.5

0.6

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

Experimental drift velocity (m/s)

(a) (b) Figure 8. CFD vs. Experimental drift velocity for the 13 cases stablished for validation. (a) Complete range. (b) Zoom on the high-viscosity cases (low drift velocities)

Figure 8 suggests that, overall, the CFD model achieves a very reasonable prediction of the drift velocity for all the 13 cases evaluated. The maximum and minimum deviations from the experimental data points were found to be around 16.6% for Case 13 (Nvis = 0.003 and Eo = 6967.138) and -0.11% for Case 5 (Nvis = 0.006 and Eo = 765.397), respectively. Furthermore, the average absolute relative error calculated was 6.48%, which confirms that, on average, the drift velocities obtained by this model would lie within the 10% error region, even considering highly viscous fluids commonly found in industrial scenarios. These particular results represent a significant improvement from other CFD models proposed previously given that the accuracy of these models tends to diminish very sharply as the viscosity of the fluid increases [4], [22], [27]. In fact, in the study of Pico et al. [27], the values of the deviation between the experimental and CFD points can reach up to ~30% when operating with fluids of similar viscosities to those of Cases 18 and 20 (the deviations obtained with these cases were 0.68% and -8.54%, respectively). From the results shown in Figure 8, it can be noted that the CFD model over-predicts the drift velocity in 9 out of the 13 validating cases. These 9 cases generally correspond to the lower viscosity and lower Nvis number ranges. On the contrary, the cases where the CFD model under-predicts the results tend to be on the higher viscosity range. However, the results also suggest that the Eo number plays a vital role in the predictive capabilities of the model. In general terms, at higher Eo numbers, the model will tend to overestimate the drift velocity, whereas, at lower values, the model will under-estimate the results. Considering this, it can be stated that, as the surface tension forces become more dominant over the gravitational forces (decreasing the Eo number), the model will over-calculate this dominance, resulting in lower drift velocities. The opposite statement is valid for the case where the gravitational forces dominate over the surface tension forces. This behaviour can be seen clearly by noticing that the highest under-prediction (-10.66% deviation) is

17

achieved with Case 16, one of the cases with the lowest Eo number and the highest over-prediction (+16.65% deviation) is achieved with Case 13, the case with the highest Eo number. Besides the selection of the best configuration for the VOF parameters, the accuracy of the results can be attributed in part to the successful modelling of the complete 3-D system. As mentioned previously, most computational studies on the matter limit their reach to 2-D simulations or half-section 3-D simulations with an axial symmetry plane. As can be observed in Figure 9 for selected cases with different Nvis and Eo numbers, the velocity fields and velocity magnitudes of both phases, as well as the volume occupied by the air bubble, differ from one section of the pipeline to the other (separated by the red dotted line). This difference appears to be more prominent at high Nvis or high Eo numbers, such as the ones found in Cases 18, 20, and 12 (Figure 9(e), (f), and (b), respectively). This suggests that, even in confined flows that exhibit a fully-laminar behaviour, the assumption of axial symmetry in a two-phase system is not entirely valid as it neglects certain physical phenomena that occur as a result of the presence of an additional phase in the system. The improvement of the prediction achieved in this study compared to the literature, especially for cases with highly-viscous liquids, is an indication that these non-symmetric phenomena are relevant for modelling and simulation purposes.

(a)

(b)

(c)

(d)

18

(e)

(f)

Figure 9. Cross-section view of the velocity magnitude and velocity contour maps of certain representative cases when the bubble had travelled a distance of 20 cm. (a) Case 1 (b) Case 12 (c) Case 2 (d) Case 7 (e) Case 18 (f) Case 20. The vertical dotted line represents the geometrical center of the pipeline. The continuous red line represents the spatial volume occupied by the gas phase.

3.3. Evaluation of the existing correlations The drift velocity correlations summarized in Table 2 were analyzed and evaluated, considering the casestudies presented in Table 1. For this subsection, the dimensionless numbers employed were defined as shown in Eq.(16)-(19). It must be mentioned that the original definition of these dimensionless numbers, as given by Moreiras (2014) [5], was simplified by assuming that ρL ― ρg ≈ ρL, given that the difference between the densities is approximately three orders of magnitude for all fluids considered. Therefore, it is important to highlight that the analysis of the existing correlations was strictly limited to experimental and 𝑘𝑔 simulated cases where the gas phase corresponds to atmospheric gas with a density around ρg ≈ 1𝑚3. Fr = vd(gD) ―0.5

(16)

Eo = gD2ρLσ ―1

(17)

Nvis = μ(ρ2LgD3) Re = vdρLDμL―1

―0.5

(18) (19)

The general behaviour of the Fr number for the 24 cases is shown in Figure 10. It can be seen that the proposed extrapolated cases expand the scope of the analysis as they include Nvis numbers that have not been evaluated previously (Nvis > 0.1), and that could be commonly found in industrial applications. These extrapolated cases also include a more diverse Eo number span, with generally lower values than the ones found in the literature.

19

0.45

0.45 0.40

0.40

Validating

Extrapolated

0.35

Extrapolated

0.30

0.30

0.25

0.25

Fr

Fr

0.35

Validating

0.20

0.20

0.15

0.15

0.10

0.10

0.05

0.05

0.00 0.0001

0.00 0.001

0.01

0.1

1

10

10

100

Nvis

1,000

10,000

Eo

(a)

(b)

Figure 10. CFD-predicted behaviour of the 𝐹𝑟 number for the 24 cases. (a) 𝐹𝑟 vs. 𝑁𝑣𝑖𝑠 and (b) 𝐹𝑟 vs. 𝐸ö

3.3.1. Benjamin (1968) & Bendiksen (1984) Correlations

0.60

0.60

0.50

0.50

0.40

0.40

0.30

0.30

Fr

Fr

Both Benjamin (1968) and Bendiksen (1984) proposed drift velocity correlations based on perfect-fluid and potential flow theory on horizontal pipelines, neglecting the contribution of frictional losses due to viscosity and assuming a stagnant liquid phase [9], [14]. Only energetic losses and dissipation associated with geometrical factors of the channel (height-to-diameter ratio) were included. As mentioned by Livinus (2018), Bendiksen’s model was developed and tested for a liquid viscosity of 0.001 Pa·s, which restricts its application to very low viscosities [3]. Furthermore, Benjamin’s analytical approach neglected the effects of viscosity and the mixing of the fluids at the interphase. As a result, drift velocities measured for fluids with high viscosities will present significant deviations from the predicted values by these correlations. As observed in Figure 11, all the CFD data points are scattered far away from the models’ prediction, especially the cases with very high Nvis or very low Eo values. The only two cases that approach the 30% deviation value correspond to the lightest oils (Oil 1 and 9), which incidentally also exhibit high Eo values. This implies that the frictional viscous losses are lower in comparison to the other cases, and the effects of the surface tension forces can be neglected, approaching the system towards a perfect-fluid scenario.

0.20

0.20 0.10 0.00 0.0001

Benjamin (1968) / Bendiksen (1984)

Benjamin (1968) / Bendiksen (1984)

0.10

CFD Data

CFD Data

0.00

0.001

0.01

0.1

1

10

Nvis

1

10

100

1000

10000

Eo

(a) (b) Figure 11. CFD data and Benjamin & Bendiksen drift velocity correlation plotted for (a) 𝐹𝑟 vs. 𝑁𝑣𝑖𝑠 and (b) 𝐹𝑟 vs. 𝐸ö

3.3.2. Weber (1981) & Wilkinson (1982) Correlations

20

Weber (1981) and Wilkinson (1982) formulated drift velocity correlations in which only capillary forces were considered through the Eo number. Weber’s correlation was developed and fitted for the low viscosity liquids flowing through horizontal pipelines of varying diameters considered in the study of Zukoski (1966) [18]. The interest of this model was to account for the effects of the surface tension and the pipe diameter on the pressure and gravitational gradients driving the motion of the bubble, neglecting the viscous contribution. A balance of forces was carried out between the capillary pressure at the tip of the bubble, given by Laplace’s equation, and the pressure exerted by the displaced liquid slug (hydrostatic pressure). This analysis led to a critical Eo value (<8.3) from which the expected drift velocity (and hence the Fr number) is zero due to dominant surface tension forces over gravitational and inertial ones, as seen in Figure 12. On the other hand, Wilkinson (1982) simultaneously proposed a very similar model with basically the same phenomenological concept. This correlation was based on the conservation of flow forces between two reference points: the nose of the bubble and a stagnation point. Between these points, a Bernoulli balance was formulated, in which the pressure term at the stagnation point was replaced with the capillary pressure inside the air cavity generated due to surface tension [20]. This balance was reformulated with dimensionless parameters and left explicitly in terms of the Fr number. Nevertheless, the specific correlation examined in Figure 12 was established based on the observation that, when the shape and dimensions of the bubble at the immediate vicinity of the stagnation point remain unchanged, the gravitational and surface tension forces are dominant on the system. 0.60 Weber (1981)

0.50

CFD Data Wilkinson (1982)

Fr

0.40 0.30 0.20 0.10 0.00 1

10

100

1000

10000

Eo

(a) (b) Figure 12. (a) CFD data and Weber/Wilkinson drift velocity correlation plotted for 𝐹𝑟 vs. 𝐸ö and (b) predicted (Correlation) vs. calculated (CFD) 𝐹𝑟 number for the case studies considered

Both correlations have a poor performance for fluids with very high viscosities, which correspond to low Fr numbers on the x-axis of Figure 12(b), as expected from their main assumptions previously discussed. None of these correlations consider the viscous effects and both highly rely on systems with dominant gravitational and surface tension forces. Therefore, fluids that exhibit both strong capillary and viscous forces, such as case 17, will not be accurately described by either of these correlations. This highlights a significant flaw in these models, demonstrating that viscous forces are fundamental in the estimation of the drift velocity on horizontal pipes. The study of Livinus (2018) also reported high errors between the predicted values with Weber’s model when compared to experimental data reported in the literature, ranging between -20% to 60% [3]. The bestfitting cases obtained in the present study (10% - 47% deviation) correspond to the smallest Nvis, either due to low viscosities or very large diameters (cases 1, 4, 11, and 12). These are the cases in which the viscous forces will not be significant in the bubble’s penetration. Wilkinson offers a slightly better prediction than Weber (~10% average improvement) since the formulation of his model does not directly neglect the viscous effects but considers a limit case in which other forces will be more significant. 3.3.3. Jeyachandra (2012) Correlation

21

Jeyachandra (2012) was the first study to point out that all existing correlations are based on Benjamin’s analysis for inviscid fluid flow or fitted with experimental data for low viscosity fluids (mostly water). Therefore, the main intent of this study was to extend the experimental scope given by Gockal et al. (2009) to different diameters and viscosities to propose an empirical closure relationship [2]. The correlation developed was the first to include the effects of viscosity and surface tension simultaneously, but its approach was entirely based on a mathematical fitting of the tendency observed for the data measured and gathered from the literature. 0.60

0.60 Jeyachandra (2012) Eo=5

0.50

CFD Data

0.40

Jeyachandra (2012) Eo=100

0.40

0.30

Jeyachandra (2012) Eo=1000

0.30

Fr

Fr

0.50

Jeychandra (2012) Eo=10000

0.20

CFD Data Jeyachandra (2012) Nvis=1E-3

0.20

0.10 0.00 0.0001

Jeyachandra (2012) Nvis=1E-4

Jeyachandra (2012) Nvis=1E-2 Jeyachandra (2012) Nvis=0.1

0.10 0.00

0.001

0.01

0.1

1

1

10

10

Nvis

100

1000

10000

Eo

(a)

(b) 0,08

Predicted Fr (Jeyachandra)

0,06

0,04 Data Points 0,02

X=Y 30% -30%

0 0

0,02

0,04

0,06

0,08

CFD Calculated Fr

(c) (d) Figure 13. CFD data and Jeyachandra drift velocity correlation plotted for (a) 𝐹𝑟 vs. 𝑁𝑣𝑖𝑠, (b) 𝐹𝑟 vs. 𝐸ö and predicted (Correlation) vs. calculated (CFD) 𝐹𝑟 number for (c) all test cases considered and (d) high viscosity fluids with low 𝐹𝑟 numbers.

Nevertheless, this correlation is in excellent agreement with the CFD data of the present study and with the experimental data collected by both Jeyachandra and Livinus. Figure 13(c) shows that most of the data points lie comfortably within the ± 30% deviation zone, especially for high Fr numbers. Excluding all highviscosity fluids (Nvis > 0.1), the absolute relative error obtained when compared to the CFD data was 11%, which is consistent with the error span for the experimental data reported in Jeyachandra (2012) (+8% to 13%) and Livinus (2018) (~20%) [2], [3]. Despite that, the average error depicted for highly viscous fluids in Figure 13(d) is significantly higher, around 45%. Some significant deviations can be observed, particularly for the most viscous fluids (cases 17, 22, 23, and 24), which coincide with the experimental conditions that were not tested by Jeyachandra (high Nvis and small diameters). This weakness in the correlation is a consequence of its mathematical construction, in which the exponential relationship was fitted considering an asymptotic region determined only by the fluids with the highest 22

viscosities tested in that study. Hence, fluids with higher Nvis than the ones explored by Jeyachandra will tend to be incorrectly predicted by this model, as well as cases with strong capillary and surface tension forces (smaller diameters). Despite this flaw, the relationship proposed performs successfully, with an average absolute error of 23% for all cases tested.

3.3.4. Moreiras (2014) Correlation Moreiras (2014) followed a similar procedure to that of the past correlation analyzed. In this study, a unified closure relationship for horizontal pipes was formulated based on the mathematical behaviour observed between the Fr and the Nvis number for his data and experimental data taken from the literature. The developed correlation follows the potential flow solution for Nvis values approaching zero and tends to a zerodrift velocity for increasing Nvis numbers [5]. However, this model was constructed on two main assumptions: i) only cases with Eö numbers above the critical value, as proposed in previous investigations [8], [15], were considered (surface tension effects were neglected), and ii) only pipe diameters equal or above to 0.0373 m are valid to be used with the correlation given. The second assumption comes as a consequence of the first one, given that at small diameters, the Eö number drops and the surface tension forces become dominant. 0.60 0.50

Fr

0.40 0.30 0.20 Moreiras (2014)

0.10 CFD Data

0.00 0.0001

0.001

0.01 Nvis (a)

0.1

1

10

0,08

Predicted Fr (Moreiras)

0,06

0,04 Data Points X=Y

0,02

30% -30% 0 0

0,02

0,04 0,06 CFD Calculated Fr

0,08

(b) (c) Figure 14. (a) CFD data and Moreiras drift velocity correlation plotted for 𝐹𝑟 vs. 𝑁𝑣𝑖𝑠 and predicted (Correlation) vs. calculated (CFD) 𝐹𝑟 number for (b) all test cases considered and (c) high viscosity fluids with low 𝐹𝑟 numbers.

23

In Figure 14(b), similar behaviour can be observed to that of Jeyachandra’s model, in which most of the data points with high Fr numbers are well situated inside the ± 30% deviation region. These data points usually correspond to the study cases with medium-ranged viscosities (0.04 Pa·s to 0.6 Pa·s) and large diameters (ID>37.3 mm) (cases 1-13). These cases are similar, or the same, as the ones used by Moreiras to adjust the correlation parameters, which explains the excellent agreement observed. If all high-viscosity cases (Nvis > 0.1) are excluded, the average absolute relative error obtained for this study’s data is 18%, which is marginally higher than the one obtained for Jeyachandra, but is consistent with Livinus (2018) reported performance for this model (ranging between -20% and 40%) [3]. However, as noticed from Figure 14(c), and in contrast to Figure 13(d), Moreiras’ proposed correlation does not provide an accurate prediction for most of the high viscosity cases analyzed in this study. The average absolute error dramatically rises over 200%, over predicting by one order of magnitude or more the drift velocity estimated through CFD. This gap occurs similarly as for Jeyachandra’s model, in which an asymptotic region for very high Nvis values is reached. Although both models share the same shortcoming, Jeyachandra’s inclusion of the Eö parameter benefits the model as it generates a steeper descent of the Fr number for increasing Nvis numbers. This characteristic allows Jeyachandra’s model to give a more accurate prediction of the drift velocity for highly viscous cases. Moreiras correlation does not reach low Fr numbers sufficiently fast for high Nvis numbers. The average absolute error obtained is 88% overall. 3.3.5. Livinus (2018) Correlation Livinus’ 2018 simplified relationship seeked to tackle all the shortcomings mentioned previously and the limitations of the most common non-complex drift velocity correlations. This study pointed out that all these correlations have always tested a narrow range of experimental data, which consequently generated low predictive capabilities for a broader scope of operating conditions [3] [29]. The corrected model was based on the fitting of a third-degree polynomial, considering the log-log relationship between the Fr number and a combination of the Eö and the buoyancy Re number (inverse of the Nvis number) [29]. The consideration of the log-log relationship generates a much better adjustment of the experimental data, given that it covers an extensive range of orders of magnitude for the Fr and the Nvis/Eö numbers. Likewise, combining the Nvis and Eö numbers in one single independent variable is very convenient, as it accounts for the effect of both, the surface tension and viscous forces on the motion of the bubble without giving a preference to either one. Thus, almost every possible experimental scenario can be well represented by this model, regardless of the dominant force acting on the penetrating bubble or the fluid properties and geometrical conditions of the pipeline. 2.00 1.80

Livinus (2018)

1.60

CFD Data

1.40

-Log10(Fr)

1.20 1.00 0.80 0.60 0.40 0.20 0.00 0

2

4

6 Log10(REo)

8

10

12

(a)

24

0,08

Predicted Fr (Livinus)

0,06

0,04

Data Points 0,02

X=Y 30% -30%

0 0

0,02

0,04

0,06

0,08

CFD Calculated Fr

(b) (c) Figure 15. (a) CFD data and Livinus drift velocity correlation plotted for –𝐿𝑜𝑔10(𝐹𝑟) vs. 𝐿𝑜𝑔10(𝑅 ∙ 𝐸ö) and predicted (Correlation) vs. calculated (CFD) 𝐹𝑟 number for (b) all test cases considered and (c) high viscosity fluids with low 𝐹𝑟 numbers.

As expected, Livinus’ relationship achieved an excellent agreement against the case-studies considered in this work, both for low to medium viscosities (Nvis < 0.1), and very high viscosities (Nvis > 0.1). On average, a relative error of 23% against all case-studies of this work was found, which is the same as Jeyachandra’s model and almost the same as the precision observed for Moreiras’ relationship, excluding the high-viscosity cases. This result is also consistent with the assessment of the model carried out by the authors of this correlation, where an average absolute relative error of 16% is reported for all the experimental data collected from multiple previous studies [29]. Very few cases (cases 17, 19, 22 and 24) exhibited higher deviations. In particular, these are cases with very dominant viscous forces and important surface tension forces, given by high Nvis (Nvis > 0.1) and low Eö numbers (Eö < 600). This limitation was already pointed out by Livinus (2018), where the performance of this correlation was analyzed for several experimental flow cases taken from the literature. It was noticed that the proposed relationship presented significant deviations for non-Newtonian fluids in small pipes (ID<0.012m), data taken from Shosho and Ryan (2001) [48], and for small diameters at a horizontal arrangement (ID<0.057m), data taken from Losi and Poesio (2016) [49]. These cases correspond to a similar range of Eö and Nvis numbers as the one mentioned previously, which accounts for the inconsistencies observed. Even though this correlation slightly out-performs Jeyachandra’s model for high-viscosity cases (5% lower average absolute relative error), it somewhat loses accuracy for other general flow cases with lower Nvis numbers (4% increase in the average absolute relative error). 3.3.6. Proposed Correlation Taking advantage of the analysis presented previously for all the existing non-complex closure relationships, a new correlation is proposed with particular emphasis on the flow cases studied in this work. Given that the complete Nvis and Eö number ranges evaluated in this study have not been considered in previous studies, the new model developed should out-perform the existing relationships’ predictions for horizontal drift velocity on highly viscous two-phase flow, a common scenario encountered nowadays in the O&G industry. Furthermore, bearing in mind that the CFD model has been thoroughly revised and adequately validated with experimental cases from the literature, the data used to test the new correlation is reliable and provides a realistic result of its expected performance. It is important to emphasize that the proposed correlation was developed exclusively for gas-liquid flows where the gas phase corresponds to atmospheric gas with a density 𝑘𝑔 around ρg ≈ 1𝑚3. Consequently, this correlation has not been developed or tested considering other heavier gases or different two-phase flow scenarios, such as liquid-liquid flows.

25

The formulation of the new model was based on a similar principle as the last three correlations developed for the drift velocity of horizontal flows (Jeyachandra (2012), Moreiras (2014) and Livinus (2018)), which is a mathematical fitting based on the relationship observed between the Fr and the Nvis/Eö numbers. As discussed throughout this section, an accurate prediction of the drift velocity must always include the influence of both, the viscous and surface tension forces. Considering that, the contemplated model accounted for both of these forces with the inclusion of the Nvis and Eö numbers. In contrast to the previous correlations, this model incorporated these numbers in separate terms, in order to attempt to consider the contribution of each set of forces separately. The expression suggested is given on Eq.(20), where the parameters were fitted by a simple method of least squares against the CFD calculated values. These parameters are given in Table 7. Fr = a ―

(

Nvis

)

b + cNvis

― (𝑒exp (E𝑓o))

(20)

Table 7. Fitted parameters for the proposed horizontal drift velocity correlation.

a

b

c

e

f

0.76419

0.02418

2.15042

0.28879

-0.50932

The first two terms of Eq.(20) are based on the relationship of Moreiras (2014), given that this study provided a better fit for most of the low-to-medium range viscosity cases in comparison to Livinus (2018) and Jeyachandra (2012). Moreover, Moreiras’ model describes more accurately experimental cases where the viscous forces are dominant, and surface tension ones can be neglected. When the Eö number increases, the exponential term of Eq.(20) (third term) tends to e, and the correlation will resemble the behaviour described by Moreiras, as a ― e ≈ 0.48, which is close to the potential flow value given by Moreiras (0.54). This third term accounts for the surface tension forces and has an exponential nature in order to correct the flaw discussed on Moreiras’ correlation regarding the asymptotic region and slow descent of the Fr number with increasing Nvis numbers. Regardless of the Eö value considered, the third term will never become zero, which causes the Fr to have a steeper decline that the one found in Moreiras’ model. This feature improves the prediction of very low Fr numbers at high Nvis numbers, but will cause a slight under-prediction for low Nvis and high Fr. 0.50

0.50

Proposed Eo=5

0.45

CFD Data

0.45

0.40

Proposed Eo=100

0.40

Proposed Eo=1000

0.35

Proposed Nvis=1E-4

0.30

0.30

CFD Data

0.25

0.25

Proposed Nvis=1E-3

0.20

0.20

Proposed Nvis=1E-2

0.15

0.15

Proposed Nvis=0.1

0.10

0.10

0.05

0.05

0.00 0.0001

Proposed Eo=10000

Fr

Fr

0.35

0.001

0.01 Nvis

(a)

0.1

1

10

0.00 1

10

100

1000

10000

Eo

(b)

26

0,08

Predicted Fr (Porposed)

0,06

0,04 Data Points X=Y

0,02

30% -30% 0 0

0,02

0,04

0,06

0,08

CFD Calculated Fr

(c) (d) Figure 16. CFD data and proposed drift velocity correlation plotted for (a) 𝐹𝑟 vs. 𝑁𝑣𝑖𝑠, (b) 𝐹𝑟 vs. 𝐸ö and predicted (Correlation) vs. calculated (CFD) 𝐹𝑟 number for (c) all test cases considered and (d) high viscosity fluids with low 𝐹𝑟.

In comparison to the previous correlations, a very good fit of the CFD data produced in this work can be observed for all cases studied. All data points lie very close to the X=Y line, with very few cases lying close to the 30% deviation lines. Only two significant deviations can be observed for this correlation, which correspond to cases 22 and 24. These two cases have been mispredicted by all models considered in this work, given that they represent very dominant viscous and surface tension forces simultaneously, for which a simple empirical relationship is not sufficient. On average, for all test cases studied, this model obtains an absolute relative error of 18%, which is 4% lower than Jeyachandra’s and Livinus’ correlations, and the same as the one calculated for Moreiras (excluding the high viscosity cases). Furthermore, the average absolute relative error obtained for this model considering high Nvis numbers (Nvis > 0.1) is 34%, which is 11% and 6% lower than the one found in Jeyachandra’s and Livinus’s studies, respectively. Likewise, for cases with lower Nvis (under 0.1), the correlation proposed does not compromise its accuracy as it achieves an average absolute relative error of 11%, which is equal to Jeyachandra’s, 4% and 7% lower than Livinus’ and Moreiras’, respectively. Finally, Figure 17 indicates that the proposed correlation predicts the highest number of cases within the 0% to 10% relative error range, which suggests a general improvement from the previous existing correlations. 18

Number of data points

16

Jeyachandra (2012) Proposed

14

Moreiras (2014)

12

Livinus (2018)

10 8 6 4 2

90 10 0

80

70

60

50

40

30

20

0 10

-1 00 -9 0 -8 0 -7 0 -6 0 -5 0 -4 0 -3 0 -2 0 -1 0

0

Vd, error (%) Figure 17. Distribution of the drift velocity errors obtained by each correlation

The performance of the different correlations was further evaluated through the Pearson linear correlation test and the Chi-squared Goodness of Fit Test (χ2 GOF). When comparing the Pearson correlation coefficients depicted in Figure 12 - Figure 16 (Cp), it can clearly be noted that Jeyachandra, Moreiras and Livinus’ models provide a substantially higher value (Cp~0.98) when compared to Wilkinson/Weber’s correlations (Cp 27

~0.66). However, the correlation proposed in this study shows a slightly higher coefficient in comparison to the former models mentioned (0.994 for the proposed model compared to the 0.989, 0.987 and 0.985 for Jeyachandra, Moreiras and Livinus, respectively). Pearson coefficients closer to 1 imply that the two variables in the correlation plots have a stronger linear correlation. In this particular application, higher Cp values imply that the predictions of the models are closer to the CFD data. This can be considered as another demonstration that the developed correlation provides a slightly better fit against the CFD data when compared against the other models. The p-values given for the Chi-Square Goodness of Fit (GOF) test show, as expected, that only Weber and Wilkinson models display statistically significant differences between the CFD data and the predicted values, with p-values lower than a significance level of 0.01.

4. CONCLUSIONS This study had the objective of developing a 3D-CFD model able to capture the complex behaviour of a twophase system flowing through a horizontal pipeline. This model was utilized to evaluate the existing correlations for the prediction of the drift velocity and to propose an improved correlation, placing particular focus on high-viscosity fluids commonly found in real world situations such as unconventional reservoirs. For these purposes, experimental data points were taken from previous studies as a means of validation for the CFD model. Additionally, various extrapolated cases were proposed in order to expand the range of analysis of the Nvis and Eö numbers. In general terms, the CFD model developed achieved high predicting capabilities given that the average absolute relative error for all 13 validating cases was 6.48%, which included cases with high-viscosity fluids (Nvis>0.1) that, before this study, had not been modelled successfully in the context of drift flux in pipelines. The improvement in the accuracy of the model was partly attributed to the correct selection of the two-phase physical models and VOF parameter calibration. This calibration indicated that, when operating with fluids of viscosities lower than 0.3 Pa·s, sharper interphases achieved through the addition of terms to the VOF equation (Sharpening Factor) are beneficial to the solution. In the case of higher-viscosity fluids, the modelling of the surface tension forces through the Continuum Surface Force model (CSF) creates parasitic currents and unrealistic high-velocity vectors at the interphase. Therefore, an additional dissipation term in the momentum equation has to be added in order to correct the behaviour and obtain more reliable results. The CFD model also allowed understanding some key features of the two-phase system, notably, that the phenomena that occur within the pipeline do not have a symmetric nature in terms of the interphase shape, bubble front and velocity fields. As a consequence, any efforts at modelling and simulating a liquid-air slugflow in pipelines with a high degree of accuracy should consider a complete 3-D system. The analysis of the seven drift velocity correlations through the 24 case-studies proposed suggested that the influence of the surface tension, the viscosity, and the capillary forces is indeed very significant in the overall behaviour of the gas bubble. Hence, the correlations that provided better predictions are mainly those that consider both, the Nvis and Eö numbers, such as Jeyachandra’s and Livinus’ correlations. In contrast, correlations such as the ones constructed by Benjamin (1968), Weber (1981), Wilkinson (1982), Bendiksen (1984) and even Moreiras (2014) will always fail at predicting the drift velocity for a wide range of cases. The results also indicated that the proposed correlation provides a simple, yet accurate-enough relationship between the drift velocity and the most relevant operating parameters and fluid properties. The proposed equation out-performs all seven previous correlations, even with high-viscosity fluids, which makes it more suitable for real-life scenarios at the industrial level. However, it should be noted that the behaviour of the system differs significantly depending on the viscosity of the fluid and a single simplified correlation would always be insufficient to cover all possible ranges. Finally, it can also be stated that the prediction by either CFD means or empirical correlations of cases in which the Nvis is high and the Eö is low (high viscosity and high surface tension) is particularly problematic given that the combined effect of these two parameters has to be taken into account. Considering this, the authors suggest that any future studies on the matter should place particular focus on these types of cases. The 3D CFD model developed and tested throughout this work represents a valuable tool to analyze comprehensively real-world situations encountered in the O&G sector. Through the computational model, it is possible to perform accurate predictions on the detailed dynamics of multiphase flows handled in upstream facilities without the need to invest large amounts of resources on in-situ experimental tests and studies. 28

Furthermore, the CFD model may provide solutions to existing problems through a robust assessment of a wide range of optimized methodologies and technologies. This represents a major advantage in O&G applications nowadays, since most of the currently available knowledge relies on empirical models and experimental data which are limited to a narrow spectrum of operational conditions and fluid properties.

References [1]

H. Pineda, T.-W. Kim, N. Ratkovich, and E. Pereyra, “CFD Modelling of High Viscosity Liquid-Gas Two-Phase Slug Flow,” Universidad de los Andes, 2016.

[2]

B. Gokcal, B. C. Jeyachandra, A. Al-Sarkhi, C. Sarica, and A. Sharma, “Drift-Velocity Closure Relationships for Slug Two-Phase High-Viscosity Oil Flow in Pipes,” SPE J., vol. 17, no. 02, pp. 593–601, 2012.

[3]

A. Livinus, P. Verdin, L. Lao, J. Nossen, M. Langsholt, and H. G. Sleipnæs, “Simplified generalised drift velocity correlation for elongated bubbles in liquid in pipes,” J. Pet. Sci. Eng., 2018.

[4]

M. Ramdin and R. Henkes, “Computational Fluid Dynamics Modeling of Benjamin and Taylor Bubbles in Two-Phase Flow in Pipes,” J. Fluids Eng., vol. 134, no. 4, p. 041303, 2012.

[5]

J. Moreiras, E. Pereyra, C. Sarica, and C. F. Torres, “Unified drift velocity closure relationship for large bubbles rising in stagnant viscous fluids in pipes,” J. Pet. Sci. Eng., vol. 124, pp. 359–366, 2014.

[6]

D. T. Dumitrescu, “Strömung an einer Luftblase im senkrechten Rohr,” ZAMM J. Appl. Math. Mech. / Zeitschrift für Angew. Math. und Mech., vol. 23, no. 3, pp. 139–149, 1943.

[7]

R. M. Davies and G. Taylor, “The Mechanics of Large Bubbles Rising through Extended Liquids and through Liquids in Tubes,” Proc. R. Soc. A Math. Phys. Eng. Sci., vol. 200, no. 1062, pp. 375–390, 1950.

[8]

E. E. Zukoski, “Influence of viscosity, surface tension, and inclination angle on motion of long bubbles in closed tubes,” J. Fluid Mech., vol. 25, no. 4, pp. 821–837, 1966.

[9]

K. Bendiksen, “An Experimental Investigation of the Motion of Long Bubbles in Inclined Tubes,” Int. J. Multiph. Flow, vol. 10, no. 4, pp. 467–483, 1984.

[10]

M. E. Weber, A. Alarie, and M. E. Ryan, “Velocities of extended bubbles in inclined tubes,” Chem. Eng. Sci., vol. 41, no. 9, pp. 2235–2240, 1986.

[11]

I. N. Alves, O. Shoham, and Y. Taitel, “Drift velocity of elongated bubbles in inclined pipes,” Chem. Eng. Sci., vol. 48, no. 17, pp. 3063–3070, 1993.

[12]

T. Taha and Z. F. Cui, “CFD modelling of slug flow in vertical tubes,” Chem. Eng. Sci., vol. 61, no. 2, pp. 676–687, 2006.

[13]

M. Ramdin and R. A. W. M. Henkes, “CFD for Multiphase Flow Transport in pipelines,” in International Conference on Ocean, Offshore and Arctic Engineering, 2011, pp. 1–11.

[14]

T. B. Benjamin, “Gravity currents and related phenomena.,” J. Fluid. Mech., vol. 31, no. 02, pp. 209– 248, 1968.

[15]

G. B. Wallis, One-dimensional two-phase flow. New York: McGraw-Hill, 1969.

[16]

A. E. Dukler and M. G. Hubbard, “A Model for Gas-Liquid Slug Flow in Horizontal and Near Horizontal Tubes,” Ind. Eng. Chem. Fundam., vol. 14, no. 4, pp. 337–347, 1975.

[17]

N. I. Heywood and J. F. Richardson, “Slug flow of air-water mixtures in a horizontal pipe: Determination of liquid holdup by γ-ray absorption,” Chem. Eng. Sci., vol. 34, no. 1, pp. 17–30, 1979.

[18]

M. E. Weber, “Drift in intermittent two phase flow in horizontal pipes,” The Canadian Journal of 29

Chemical Engineering, vol. 59, no. 3. pp. 398–399, 1981. [19]

G. C. Gardner and I. G. Crow, “The motion of large bubbles in horizontal channels,” J. Fluid Mech., vol. 43, no. 2, pp. 247–255, 1970.

[20]

D. L. Wilkinson, “Motion of air cavities in long horizontal ducts,” J. Fluid Mech., vol. 118, pp. 109– 122, 1982.

[21]

B. Gockal, Q. Wang, H.-Q. Zhang, and C. Sarica, “Effects of High Oil Viscosity on Oil/Gas Flow Behaviour in Horizontal Pipes,” Soc. Pet. Eng. Proj. Facil. Constr., pp. 1–11, 2008.

[22]

P. Andreussi, M. Bonizzi, and A. Vignali, “Motion of elongated gas bubbles over a horizontal liquid layer,” 14th Int. Conf. Multiph. Prod. Technol., pp. 309–318, 2009.

[23]

R. F. Kroes and R. A. W. M. Henkes, “CFD for the motion of elongated gas bubbles in viscous liquid,” BHR Gr. - 9th North Am. Conf. Multiph. Technol. 2014, no. 1, pp. 283–298, 2014.

[24]

M. Lu, “Experimental and computational study of two-phase slug flow,” Imperial College of London, 2015.

[25]

B. Sanderse, M. Haspels, and R. A. W. M. Henkes, “Simulation of Elongated Bubbles in a Channel Using the Two-Fluid Model,” J. Dispers. Sci. Technol., vol. 36, no. 10, pp. 1407–1418, 2015.

[26]

M. Tshuva, D. Barnea, and Y. Taitel, “Two-phase flow in inclined parallel pipes,” Int. J. Multiph. Flow, vol. 25, no. 6–7, pp. 1491–1503, 1999.

[27]

P. D. Pico, J. P. Valdés, N. Ratkovich, and E. Pereyra, “Analysis of the Drift Flux in Two-Phase GasLiquid Slug-Flow along Horizontal and Inclined Pipelines through Experimental (Non)-Newtonian and CFD Newtonian Approaches,” J. Fluid Flow, Heat Mass Transf., vol. 5, pp. 53–70, 2018.

[28]

M. A. B. Khan, A. K. Mehrotra, and W. Y. Svrcek, “VISCOSITY MODELS FOR GAS-FREE ATHABASCA BITUMEN.,” J. Can. Pet. Technol., vol. 23, no. 3, pp. 47–53, 1984.

[29]

A. Livinus, “A STUDY OF A SINGLE ELONGATED BUBBLE IN LIQUID IN PIPES,” Cranfield University, 2018.

[30]

B. Andersson, R. Andersson, L. Håkansson, M. Mortensen, R. Sudiyo, and B. Van Wachem, Computational fluid dynamics for engineers. 2011.

[31]

C. W. Hirt and B. D. Nichols, “Volume of fluid (VOF) method for the dynamics of free boundaries,” J. Comput. Phys., vol. 39, no. 1, pp. 201–225, 1981.

[32]

H. Pineda et al., “Phase distribution analysis in an Electrical Submersible Pump (ESP) inlet handling water-air two-phase flow using Computational Fluid Dynamics (CFD),” J. Pet. Sci. Eng., vol. 139, pp. 49–61, 2016.

[33]

H. Pineda-Pérez, T. Kim, E. Pereyra, and N. Ratkovich, “CFD modeling of air and highly viscous liquid two-phase slug flow in horizontal pipes,” Chem. Eng. Res. Des., 2018.

[34]

M. A. Reyes-Gutiérrez, L. R. Rojas-Solórzano, J. C. Marín-Moreno, A. J. Meléndez-Ramírez, and J. Colmenares, “Eulerian-Eulerian Modeling of Disperse Two-Phase Flow in a Gas-Liquid Cylindrical Cyclone,” J. Fluids Eng., vol. 128, no. 4, p. 832, 2006.

[35]

H. Yang, J. Lu, and S. Yan, “Preliminary Numerical Study on Oil Spilling from a DHT,” TwentyFourth Int. Ocean Polar Eng. Conf., vol. 3, no. 1994, pp. 610–617, 2014.

[36]

V. R. Gopala and B. G. M. van Wachem, “Volume of fluid methods for immiscible-fluid and freesurface flows,” Chem. Eng. J., 2008.

[37]

S. Jeon, S. Kim, and G. Park, “CFD simulation of condensing vapor bubble using VOF model,” Proc. World Acad. Sci. …, vol. 3, no. 12, pp. 209–215, 2009.

30

[38]

Y. Tanabe, “Assessment of Volume of fluid Method for high-pressure gas injection into liquid,” Royal Institute of Technology. Stockholm, Sweden, 2015.

[39]

S. Muzaferija and M. Peric, “Simulation of Free Surface Flows With STAR-CCM+,” in STAR Japanese Conference 2013, 2013.

[40]

Y. Tanabe, “Assessment of Volume of fluid Method for high-pressure gas injection into liquid,” Royal Institute of Technology. Stockholm, Sweden, 2015.

[41]

N. M. C. Martins, N. J. G. Carriço, H. M. Ramos, and D. I. C. Covas, “Velocity-distribution in pressurized pipe flow using CFD: Accuracy and mesh analysis,” Comput. Fluids, 2014.

[42]

N. M. C. Martins, A. K. Soares, H. M. Ramos, and D. I. C. Covas, “CFD modeling of transient flow in pressurized pipes,” Comput. Fluids, 2016.

[43]

M. Spiegel et al., “Tetrahedral vs. polyhedral mesh size evaluation on flow velocity and wall shear stress for cerebral hemodynamic simulation,” Comput. Methods Biomech. Biomed. Engin., 2011.

[44]

I. B. Celik, U. Ghia, P. J. Roache, C. J. Freitas, and H. C. Raad, “Procedure for Estimation and Reporting of Uncertainty Due to Discretization in CFD Application,” J. Fluids Eng., vol. 130, 2008.

[45]

L. Kwasniewski, “Application of grid convergence index in FE computation,” Bull. Polish Acad. Sci. Tech. Sci., 2013.

[46]

C. J. Roy, “Grid convergence error analysis for mixed-order numerical schemes,” AIAA J., 2003.

[47]

N. Ratkovich, S. K. Majumder, and T. R. Bentzen, “Empirical correlations and CFD simulations of vertical two-phase gas-liquid (Newtonian and non-Newtonian) slug flow compared against experimental data of void fraction,” Chem. Eng. Res. Des., 2013.

[48]

C. E. Shosho and M. E. Ryan, “An experimental study of the motion of long bubbles in inclined tubes,” Chem. Eng. Sci., vol. 56, no. 6, pp. 2191–2204, 2001.

[49]

G. Losi and P. Poesio, “An experimental investigation on the effect of viscosity on bubbles moving in horizontal and slightly inclined pipes,” Exp. Therm. Fluid Sci., vol. 75, pp. 77–88, 2016.

31

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

32

Evaluation of the existing drift-velocity closure relationships considering highly viscous gas-liquid two-phase slug flow in horizontal pipelines through 3-D CFD modelling Juan Pablo Valdés1, Paula Pico1*, Eduardo Pereyra2, Nicolás Ratkovich1 1 Department

2 McDougall

of Chemical Engineering, Universidad de los Andes, Carrera 1 # 18a-12 Bogota, Colombia

School of Petroleum Engineering, The University of Tulsa, Tulsa, OK 74104, United States *Corresponding Author: [email protected]

Highlights     

A 3D-CFD model was developed to predict the drift velocity of highly-viscous fluids. Seven drift velocity correlations were evaluated at different dimensionless numbers. High predicting capabilities found for the CFD model (6.48% average error). Results show that both, Viscosity and Eötvös numbers, influence the drift velocity. A new correlation was proposed, improving the prediction compared to other works.

33

Contributor Roles Taxonomy (CRediT) author statement The individual contributions of each author to the development of the manuscript are as following: Juan Pablo Valdés: Methodology, Formal analysis, Investigation, Writing – original draft and review, and Project Administration. Paula Pico: Methodology, Formal analysis, Investigation, Writing – original draft and review, and Project Administration. Eduardo Pereyra: Conceptualization, Supervision, and Resources. Nicolás Ratkovich: Conceptualization, Supervision, and Resources.

34