CFD modeling of air and highly viscous liquid two-phase slug flow in horizontal pipes

Accepted Manuscript Title: CFD modeling of air and highly viscous liquid two-phase slug flow in horizontal pipes Authors: H. Pineda-P´erez, T. Kim, E. Pereyra, N. Ratkovich PII: DOI: Reference:

S0263-8762(18)30310-1 https://doi.org/10.1016/j.cherd.2018.06.023 CHERD 3233

To appear in: Received date: Revised date: Accepted date:

13-2-2018 30-5-2018 14-6-2018

Please cite this article as: Pineda-P´erez, H., Kim, T., Pereyra, E., Ratkovich, N., CFD modeling of air and highly viscous liquid two-phase slug flow in horizontal pipes.Chemical Engineering Research and Design https://doi.org/10.1016/j.cherd.2018.06.023 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.

CFD modeling of air and highly viscous liquid two-phase slug flow in horizontal pipes Pineda-Pérez, H.1, Kim, T.2,*, Pereyra, E.2, Ratkovich, N.1 1

Department of Chemical Engineering, Universidad de los Andes, Bogotá, Colombia McDougall School of Petroleum Engineering. The University of Tulsa, Tulsa, OK 74104, United States *Corresponding Author, [email protected] (Kim, T. W.), Tel: +1-918-934-0052 E-mail address: [email protected] (Pineda-Pérez, H.), [email protected] (Kim, T. W.), [email protected] (Pereyra, E.), [email protected] (Ratkovich, N.)

SC RI PT

2

Highlights

TE

D

M

A

N

U

Graphical abstract

CFD approach was used to simulate air-highly viscous oil two-phase slug flow.



Comparison of CFD with experimental PIV measurements was performed

EP



CC

satisfactorily.

A



Effect of sharpening factor in VOF model showed penetrating bubble into slug body.



Velocity profile is not fully developed near the bottom of slug body.



Velocity fields showed regions where gas bubbles are entrained in slug body.

1

Abstract

SC RI PT

Computational Fluid Dynamics (CFD) approach, encoded on STAR-CCM+ was used to simulate air and highly viscous liquid two-phase slug flow in a 50.8-mm (2-in.) ID and 18.9-m (62 ft.) long horizontal pipe. The Volume of Fluid (VOF) method, the Continuum Surface Force (CSF) model and the High Resolution Interface Capturing (HRIC) scheme were utilized. The liquid viscosity varied from 161 to 567 mPa·s based on the experimental data. A sensitivity analysis of the sharpening factor, the angle factor, and the Interface Momentum Dissipation (IMD) model was carried out to obtain the best combination of parameters. Besides, comparison with experimental measurements concerning slug frequency, average liquid holdup and velocity profiles in the liquid region was carried out. The comparison shows a fair correspondence in terms of behavior and magnitude with the experimental results obtained by Particle Image Velocimetry (PIV). It was observed that the velocity profile is not fully developed near the bottom of slug body. Finally, extrapolation to the extended velocity conditions was performed. The reduction of translational velocity at higher mixture velocity conditions was investigated. A certain amount of gas passing from the tail to the front of the liquid slug body was visualized by three-dimensional velocity profiles. This phenomenon supports the reduction of translational velocity above the certain level of mixture velocity.

U

Keywords

N

Two-phase slug flow; Horizontal pipe flow; High viscous liquid; Computational Fluid Dynamics; Liquid holdup; Velocity profile

TE

D

M

Cross-sectional area occupied by liquid-phase, Cross-sectional area occupied by gas-phase, Slug frequency, Liquid holdup, Average liquid holdup, Viscosity, Liquid viscosity, Mixture Reynolds number, Density, Surface tension, Wall shear stress, Volume occupied by gas-phase, Volume occupied by liquid-phase, Liquid velocity in 2-D windows of PIV and CFD, Mixture velocity, Superficial gas velocity, Superficial liquid velocity, Translational velocity,

EP

A

CC

AL AG fS HL HL Avg μ μL ReM ρ σ τW VG VL vL vM vSg vSL vT

A

Nomenclature

m2 m2 1/s mPa·s mPa·s kg/m3 N/m Pa m3 m3 m/s m/s m/s m/s m/s

Introduction Liquid-gas two-phase flows are common phenomena in the chemical, nuclear, and petroleum industries. They occur owing to the mixture or change of different phases. However, their understanding is limited compared to the single-phase liquid flow due to its complexity. Accordingly, many efforts have been made theoretically, experimentally and numerically to study two-phase flows. 2

U

SC RI PT

Despite all these efforts, a deeper understanding of inner flow field of two-phase flows is desired to achieve a better design of production and transportation system. When two-phase liquid and gas flow occurs inside a pipe, the phases can be distributed in diverse ways resulting in different flow patterns. Flow patterns depend on the operational conditions and fluids properties, and each flow pattern exhibits different flow characteristics. For horizontal pipes, the flow patterns can be divided into the dispersed bubble flow, stratified smooth- and stratified wavy flow, slug flow (also known as intermittent flow) and annular flow. An illustration of the phase distribution corresponding to each flow pattern is presented in Figure 1. Knowledge of the flow patterns would improve the design and sizing of two-phase flow systems (Pereyra et al., 2012)26.

N

Figure 1. Flow patterns presented in horizontal pipes: (a) Dispersed bubble flow (b) Slug flow (c) Stratified smooth flow, (d) Stratified wavy flow and (e) Annular flow.

CC

EP

TE

D

M

A

The slug flow is the most commonly observed flow patterns in two-phase flow. It is one of the most complex flow patterns due to the unsteady behavior of the phase distribution. This flow pattern is defined as alternating liquid slugs and gas bubbles flowing above liquid film, and a combination of them is a slug unit (Al-Safran, 2009)2. Figure 2 shows a representation of the slug unit, where the main features are; the slug length, the bubble length and the parts of the slug (front, body, and tail).

Figure 2. Illustration of the slug unit.

A

Several empirical and mechanistic models of liquid and gas two-phase slug flow have been developed for low viscosity liquids. Dukler and Hubbard (1975)9 proposed models to predict the characteristics of the slug flow based on flow rates, fluid properties, and pipe geometry. In addition, the comparisons of model predictions with experimental results for water and air two-phase flow were performed, showing fair agreements. Gregory et al. (1978)14 conducted experiments for light oil and air two-phase flow where slug characteristics were measured in 25.8-mm (≈ 1 in) and 51.2-mm (≈ 2 in) ID pipes. With the obtained data, they suggested a simple correlation for the prediction of slug liquid holdup with mixture velocity as the only variable. Taitel and Barnea (1990)31 suggested a method for pressure drop calculations in steady slug flow. They took into account the terms (such as the hydrodynamics in the liquid film) that were not considered by previous models. Zhang et al. (2003)36 presented a detailed mechanistic model for prediction of flow pattern transitions, liquid 3

A

CC

EP

TE

D

M

A

N

U

SC RI PT

holdups and slug characteristics in liquid and gas two-phase flow at all pipe inclinations. They validated the developed model for relatively low viscosity fluids (μL < 15 mPa·s). For high viscosity liquids, two studies are important to highlight, namely, Colmenares et al. (2001)7 and Al-Safran et al. (2015)3. Colmenares et al. (2001)7 modified the model proposed by Taitel and Barnea (1990)31 to estimate pressure drop. They compared the modified model with high viscosity (480 mPa·s) liquid and air two-phase flow, resulting in 6% of the average absolute error. Al-Safran et al. (2015)3 proposed an empirical non-linear regression correlation to predict slug liquid holdup for high viscosity (180-587 mPa·s) liquid and gas two-phase flow. They showed that the slug liquid holdup increases as the liquid viscosity increases especially beyond the critical aeration mixture velocity (approximately vM ≈ 2m/s). Recently, Computational Fluid Dynamics (CFD) studies of two-phase flow have been spotlighted utilizing various approaches. These approaches are the Eulerian-Lagrangian, Eulerian-Eulerian, and the Volume Of Fluid (VOF) methods. The Eulerian-Lagrangian approach is applied to the cases where there is a dispersed phase such as bubbly flows or liquid droplets. This approach does not allow the simulation of large interfaces. The Eulerian-Eulerian approach uses one set of equations for each fluid and its response is time averaged. Similar to the Eulerian-Lagrangian approach, this time averaging does not allow the application of the Eulerian-Eulerian approach when there is a large interface between the fluids. Lastly, the VOF approach uses only one set of equations only for the continuous phase, while the dispersed phase is modeled by using an additional transport equation for its volume fraction. This approach enables the track of the interface. Hence, it is suitable to analyze two-phase flow in pipes, since it is capable to predict the flow pattern and track the liquid holdup. For horizontal two-phase flow, De Schepper et al. (2008)8, Lo and Tomasello (2010)19 and Andersen (2012)4 have carried out CFD simulations. De Schepper et al. (2008)8 modeled all liquidgas and liquid-vapor flow patterns that were reported by the Baker Chart. They applied the VOF method and the Piecewise Linear Interface Calculation (PLIC) method. Their results presented a good agreement regarding the flow pattern. Lo and Tomasello (2010)19 simulated stratified flow using the VOF method, and showed the effects of the turbulence model on the CFD results. They obtained better results using the κ-ω turbulence model. Andersen (2012)4 performed 1-D, 2-D and quasi 3-D simulations using LedaFlow and Fluent with the VOF method, and compared the results from each other. On the other hand, the validation against experimental results was not carried out. On the experimental side, Gokcal et al. (2008)13 and Brito (2012)5 focused on the effect of liquid viscosity on slug flow characteristics in the horizontal condition. These two studies considered viscosities in the range of 180-587 mPa·s and 40-160 mPa·s, respectively. The relatively low mixture velocities and high viscosity liquids yielded laminar flows (or low mixture Reynolds numbers). They showed that the slug length reduces as the viscosity increases. This observation leads the question about the existence of fully developed velocity profile within the slug body owing to the shorten slug lengths. Especially, Brito (2012)5 presented that the velocity profile within the slug is not symmetric. The recordings by high-speed camera illustrated that lower velocities are observed near the bottom as compared to the upper part of the pipe. Furthermore, Brito (2012)5 observed that the ratio between slug translational velocity and mixture velocity reduces after about 2 m/s of mixture velocity. This observation can be related to the gas blowing through the top of the slug body, which can also affect the velocity profiles in liquid slug body. However, the in-situ velocity profiles both in liquid slug body and liquid film have remained doubtful. For the study of the velocity profile of the liquid phase in two-phase flow, Particle Image Velocimetry (PIV) and ultrasonic velocity profiling have been used before. Examples of these techniques are the studies of Nogueira et al. (2003)25 and Fan et al. (2013)11. Nogueira et al. (2003)25 used the PIV technique with the Pulsed Shadow Technique (PST). They obtained the velocity profiles and the bubble shape in a vertical pipe with an aqueous solution of glycerol and air two-phase flow. Fan et al. (2013)11 presented an approach using continuous wave ultrasonic Doppler technique to investigate the velocity characteristics of slug body and film in two-phase flow.

4

M

Experimental Programs

A

N

U

SC RI PT

Kim et al. (2018)18 performed experimental investigations of highly viscous liquid and air twophase slug flow in horizontal pipes. They adopted PIV as primary instrumentation. Their experiments helped to get insight into the velocity fields within the high viscous liquid slug body and film regions. The experimental results showed that the liquid slugs are not fully developed near the pipe bottom. They presented that a low momentum region at the slug front beneath the mixing zone towards the lower sections of the pipe leads to a continuous increase of velocity within slug body. This indicates that a new modeling approach may need to be proposed to modify the existing slug flow models for highly viscous conditions. Although there are many empirical and mechanistic studies about two-phase slug flow dynamics, in-depth study into the slug flow with high viscous liquids is still desired. As commented by many previous studies, the experimental measurement of the velocity profiles inside the two-phase flow is relatively difficult owing to the opacity of the fluid and limitations of apparatus. Nevertheless, Kim et al. (2018)18 has successfully applied the PIV technique to highly viscous oil. Whereas, the limited laser frequency made it difficult to acquire data for higher mixture velocity conditions (vM > 1.5). To overcome this limitation, CFD has been envisioned to extend experimental observations reported by Kim et al. (2018)18 to higher mixture velocity conditions. In this study, an analysis of air and highly viscous two-phase slug flow was performed by CFD with the VOF method in terms of the liquid holdup and the velocity profiles. A sensitivity analysis of different simulation parameters (i.e., sharpening factor, angle factor and interface momentum dissipation model) was carried out. In addition, a comparison with experimental measurements acquired at the Tulsa University Fluid Flow Project (TUFFP) using PIV and CPs was made. This comparison was performed under the low velocity conditions. After that, the CFD model was further applied to estimate the velocity profile in higher velocity conditions. This extrapolation was verified by the measured slug characteristics.

A

CC

EP

TE

D

The information about the experimental facility and instrumentation are described in this section. The experimental data used in this study were collected by the TUFFP as the subsidiary experiments of Kim et al. (2018)18. Two-phase slug flows were analyzed by using PIV for velocity profiles and using CPs for liquid holdups as a function of time. The experimental results allowed the validation of proposed CFD model with lower superficial velocity conditions. The experimental facility is schematically shown in Figure 3.

Figure 3. Schematic view of the experimental facility (Kim et al., 2018)18.

5

Table 1. Experimental range. Variable

Minimum

vSL [m/s]

0.1

vSg [m/s]

0.1

SC RI PT

The flow loop consists of 21.9-m (72 ft.) pipe and the test section is 18.9-m (62 ft.) long with a 50.8-mm (2-in.) ID pipe. The liquid is pumped by a 20-hp Kral screw pump from the main storage tank (3.03 m3) to the flow loop. Air is delivered to the test section using a dry rotary screw compressor, which has a capacity up to 1030 CFM at 100 psig. The liquid and the gas are mixed at a tee junction before the test section. The measurements of the mass flow of each phase and their densities are performed by using Micro Motion mass flow meters with ±0.05% uncertainty of measurement. In the test section, three pairs of Capacitance Probes (CPs) are arbitrarily located at 1.8, 10.3 and 14.7 m to observe the flow development and evolution through the pipe. The point where the PIV instruments are placed is located at 14 m from the mixing tee. After the mixture passes through the test section, it returns through a 76.2-mm (3-in.) ID pipe to the main storage tank, where the phases are separated by gravity. The experimental range is presented in Table 1 regarding the superficial liquid and gas velocities, vSL and vSg, respectively.

Maximum 0.8

U

3.5

A

N

The fluid properties of each experimental case were measured and reported. The average values are shown in Table 2. The liquid viscosity varied from 161 to 567 mPa·s, and this variation was achieved by heating the oil in a temperature range from 24 °C to 49 °C. Table 2. Average properties for the fluids used experimentally.

M

Property

3

Average density (ρ) [kg/m ] Average viscosity (μ) [mPa·s]

Oil

1.4 1.85×10

840.9 -2

161-567

0.033

D

Surface Tension (σ) [N/m]

Air

A

CC

EP

TE

The PIV technique was used to obtain the velocity profiles as a function of time. The tracer particles (silver-coated hollow glass sphere type) were mixed with the liquid phase. Images of the moving particles illuminated by laser were captured with a CCD camera. Two images in a pair have the micro-seconds time gaps, allowing to calculate the velocity magnitude and direction of the tracer particles. The applicability of these particles to represent the exact motion of liquid phase was verified by their settling velocity and frequency response. Further single-phase tests also support the validation of PIV technique for this study. More detailed descriptions of PIV can be found in Kim et al. (2018)18, and the method to estimate the uncertainty of the measured velocity field can be found in several literatures such as Raffel et al. (2007)27. After that, the velocity profiles obtained by PIV were compared with the results from CFD model. The CPs were placed in the test section to measure the liquid holdup in the pipe. The uncertainty of measured liquid holdup was estimated by the standard deviations of repeatedly conducted static calibration and two-phase experiments, which was between 15% and 30% of measured holdup depending on operational conditions.

CFD Modeling A set of reliable data were obtained by experiments, while the operational conditions were limited to lower superficial velocities (vM < 1.2 m/s). CFD simulations were then purposed to analyze and study two-phase slug flow with the extended range, namely, higher superficial velocities. The simulations were performed by the commercial software STAR-CCM+ v10.02.010 (Siemens) which is capable of doing all the CFD processes from the geometry generation to the post-processing of the results in a single interface. 6

Spatial Discretization

N

U

SC RI PT

The test section of the experimental loop was simulated to replicate the experimental facility of TUFFP. The geometry model consisted of a 50.8-mm (2-in.) ID and 18.9-m (62 ft.) long horizontal pipe as same as the experimental condition. Nevertheless, only half of the pipe cross section was simulated due to computing capabilities and the required number of simulations for this study (Figure 4). The geometry was divided into three sections; the first one was the development section which is 14-m long. The second section was 0.05-m long where the PIV instrumentation is placed. The last was a downstream section with 3.2-m long. The CPs were placed in the same position as the experiment. Figure 4 presents the part of modeled geometry and its surfaces.

A

Figure 4. Part of modeled geometry.

CC

EP

TE

D

M

It was necessary to discretize the domain in smaller volumes to solve the governing equations of the model. Multiple types of mesh could be exercised, but the orthogonal mesh has been recommended for the simulation of two-phase flow in pipes (Hernandez-Perez et al., 2011)16. Figure 5 illustrates an example of the cross-sectional and longitudinal discretization applied for this study.

Figure 5. Examples of the applied orthogonal mesh. (a) Cross-sectional and (b) Longitudinal discretization.

A

The effects of the discretization of cross-sectional and longitudinal directions were investigated to select the proper mesh size. Simulations were carried out by changing the number of divisions in those directions. Figure 6 shows the obtained average liquid holdup based on 863, 1725 and 3450 divisions in the longitudinal direction with 150 divisions in the cross section, for the total number of cells of 129450, 258750 and 517500, respectively. It was observed that a higher number of division in the longitudinal direction reduced the difference in liquid holdups between CFD and experimental (plotted as a line) results. As a result, 3450 divisions in the longitudinal section were selected to reduce errors.

7

HL CFD = 0.7659 (1725 divisions)

HL Avg , (-)

0.8

HL CFD = 0.7228 (3450 divisions)

0.75 0.7 0.65

HL CFD = 0.7660 (863 divisions)

0.6 1

2

3

4

5

6

7

8

SC RI PT

Number of cells x105, (-) Figure 6. Average liquid holdup as a function of the number of cells, changing only the longitudinal divisions.

HL CFD = 0.7009 (96 divisions)

N

HL CFD = 0.7228 (150 divisions)

HL CFD = 0.7200 (216 divisions)

A

0.75 0.7 0.65

M

HL Avg , (-)

0.8

U

To test the effect of the cross-sectional discretization, three cross section divisions were tested (96, 150 and 216) with 3450 division in the longitudinal direction, for the total number of cells of 331200, 517500 and 745200, respectively. Figure 7 shows the effect of the cross-sectional divisions on the average liquid holdup, showing the unclear trend among the results. A higher number of the division would produce a larger computational time, while the coarse divisions would reduce the resolution of the velocity profile in the test section. Therefore, 150 divisions were selected in the cross-sectional area based on the balance between computational time and resolution, with the total of 517500 cells.

0.6 1

2

3

4

5

6

7

8

TE

D

Number of cells x105, (-) Figure 7. Average liquid holdup as a function of the number of cells, changing only the cross-sectional divisions.

Physical Model Selection

EP

The main physical models that emulate the flow of high viscosity oil and air are described in this section.

Physical VOF Model

A

CC

The VOF method was used to model the two-phase flow phenomenon and capture its time-dependent behavior. In this method, a single set of momentum and turbulence equations are considered for the continuous phase. The dispersed phase is modeled by using a transport equation for its volumetric fraction (Hirt and Nichols, 1981)17. The HRIC scheme was used to capture a sharp interface. The VOF model has been used to simulate various cases in the oil and gas industry, such as vertical jets (Yuan, 2006)35, oil spilling (Yang et al., 2014)34, sloshing (Monroy and Sato, 2003)23, and cyclones (Reyes-Gutiérrez et al., 2005)29. The mentioned governing equations for momentum and transport conservation are described by Equations 1 and 2, respectively.   P   U i U j U j  U iU j      t xi x j xi  x j xi









    g j  Fj ,  

(1)

and  U i   0. t xi

8

(2)

where t is time, U is velocity, P is pressure, and g represents gravitational acceleration. F is the external force per unit volume, µ is the dynamic viscosity and ρ is density. Additionally, it should be noticed that the density and dynamic viscosity are volume average values that can be formulated by Equations 3 and 4, respectively. n

    q q ,

(3)

q 1

and n

    q q .

SC RI PT

(4)

q 1

Surface Tension Model

CC

EP

TE

D

M

A

N

U

The phase interaction should be specified for the VOF model. In this case, the surface tension effects are considered. The Continuum Surface Force (CSF) model was applied to calculate the surface tension force, regarding the oil as the primary phase and the air as the secondary phase. The problem with the CSF model is the generation of parasitic currents (Harvie et al., 2006)15. They are the unphysical currents generated by a discontinuity in the domain, such as an interface, and cannot be solved by a finer mesh. The Interface Momentum Dissipation (IMD) model is available to reduce those currents adding an extra dissipation near the free surface. This dissipation is proportional to the velocity gradients and the interface artificial viscosity. In this study, the use of the IMD model was tested with the value of 0.1 as the interface artificial viscosity (CD-adapcoTM, 2015)32. Another parameter tested in the CFD modeling is the sharpening factor. This factor is used to reduce the numerical diffusion in the domain adding an extra term in the VOF transport equation. This factor varies from 0.0 to 1.0, where 0.0 means that there is no reduction in numerical diffusion (CD-adapcoTM, 2015)32. The suitable results may be obtained by applying the value of 0.0 with the second order discretization scheme, while the sharper interface can be obtained with higher values of sharpening factor. Accordingly, 0.0 and 1.0 values were tested for this factor. In the HRIC scheme, the downwind scheme accurately captures interfaces perpendicular to the flow direction. However, if the interface is parallel to the flow direction the HRIC scheme tends to wrinkle it and aligns it with the mesh (Waclawczyk and Koronowicz, 2008)33. Final correction of the normalized face value was processed by using the angle factor, which affected the blending between the used downwind and the upwind scheme. The default value in STAR-CCM+ is 0.05, and two levels (0.05 and 0.2) were used to test. A segregated model was applied to solve the flow equations, using the second order upwind convection scheme. Since the two-phase phenomenon is time-dependent, an implicit unsteady study was performed. To reduce the computational effort and due to the low velocities presented inside the pipe, both liquid and gas phases were modeled with constant densities. The considered densities and viscosities corresponded to the experimental values. In a horizontal pipe, however, a large difference between the liquid and gas densities could affect the flow behavior. The gravity was considered to capture this phenomenon.

Turbulence Model

A

The κ-ω SST (Shear-Stress Transport) model was used concerning the turbulence. The κ-ω model has the advantage over the alternate κ-ε model, which can be applied to the viscous regions of the flow, e.g. the boundary layer near the pipe walls, without further modifications (Menter, 1992)21. The model, however, presents the problem that the boundary layer computation is very sensitive to the inlet boundary conditions. A modification of the model was implemented and then, the κ-ω SST was proposed, which does not have this shortcoming (Menter, 1992)22. This model blends the κ-ε and the κ-ω model, using the first one in the core of the domain and the second one near the walls. The high viscous liquid could yield a laminar flow regime in the liquid phase. Therefore, a γ (gamma) transition model was used for the transition between the laminar and turbulent regimes, the intermittency (γ) acts on the production term of the turbulent kinetic energy transport equation in the 9

SST model to simulate laminar and turbulence flows (Malan et al., 2009)20. It is important to highlight that the used turbulence model was shared by the phases since there was just one set of equations to solve. This model can be presented by Equations 5 to 7, where Equation 7 represents the effective viscosity of the fluid.  U U j   U j k  t  i   k    x j t x j xi 

 U i   k      *  k      t  ,  x j x j    k  x j         U U j  U i   U j  C1t  i    2      t  ,    k  x j xi  x j x j      x j 







(5)

SC RI PT

      t x j



and eff    t .

(6)

(7)

where k is the turbulent kinetic energy, ω is the specific dissipation rate and µt is the turbulent viscosity. σk, σω, C1, β and β* are constants.

Boundary and Initial Conditions Specification

CC

EP

TE

D

M

A

N

U

Twelve surfaces were presented in the domain. Four of those surfaces corresponded to the interfaces between the sections, and the remaining eight were divided into one inlet, one outlet, three pipe walls and three symmetry planes. A velocity inlet condition was specified for the inlet surface. In this type of boundary condition, both the velocity profile and the volumetric fraction of the phases must be defined. The pipe crosssectional area was divided for the volumetric fraction. The upper half of the surface had a liquid holdup of 0.0, while the lower half had a liquid holdup of 1.0 exemplifying the boundary condition (Figure 8). This division means that the gas enters the domain in the upper half of the pipe and the liquid in the lower half. The two times of the superficial velocities were placed for the liquid and the gas in the lower and upper half of the inlet surface, respectively, to match the volumetric flows with the experimental conditions. The outlet surface was modeled as a pressure outlet with a constant pressure of 0 Pa. The perimeter areas (or the pipe wall) were modeled by assuming no-slip conditions. The symmetry planes were applied since only half of the pipe was modeled.

A

Figure 8. Example of the liquid volume fraction boundary condition at the inlet.

Although the expected results must be invariant of the initial conditions, the right selection of these conditions could reduce the computational time required to obtain a developed two-phase flow. For this, the initial values of the holdup in the entire domain has a value of 0.0 in the upper half of the pipe and 1.0 in the lower half. In the same direction, two times of superficial liquid and gas velocities were specified in the lower half and upper half of the pipe, respectively.

Processing The simulations presented in this study are about transient flows. Therefore, a time step, inner iterations per time step and the maximum physical time must be specified. 10

SC RI PT

The Courant-Friedrichs-Lewy (CFL) condition was used to calculate the time step. It was calculated considering the length in the axis direction of the cells, which is 5 mm. In addition, the used velocity was twice the superficial velocity of the gas to ensure a low CFL number. For near the interface, it was recommended to use CFL numbers below 0.1 to capture the interface accurately. However, it was not possible to know the CFL number as it is a function of the velocity field. Then, a CFL number of 0.02 was defined to calculate the time step using the maximum velocity that was specified in the inlet condition. Based on this approach, an average CFL number of 0.025 was obtained in the simulations, oscillating between 0.0 and 0.3 in the entire domain and below 0.1 near the interface. More detailed descriptions of CFL number can be found in Ratkovich et al. (2013)28. The five times of inner iterations were used per time step. This number allowed the RMS (Root Mean Squared) residuals of the equations to go below 0.0001 for all time steps. The physical time depends on the mixture velocity, i.e. 3 seconds is required for the case of vSL = 0.2 m/s and vSg = 5 m/s, while 50 seconds is required for the case of vSL = 0.2 m/s and vSg = 0.2 m/s. After this time, additional 15 seconds was simulated to ensure that the flow was fully developed and the results were not influenced by the initial conditions. As a result, the former and the latter cases were simulated for 18 and 65 seconds, respectively.

Post-processing

EP

TE

D

M

A

N

U

The main results from the simulations were the velocity profiles and the volumetric fraction of the liquid (liquid holdup, HL). The volumetric fraction of liquid in the positions of the capacitance probes was retrieved for the liquid holdup. Images were retrieved in the PIV section for the velocity profile. Since the velocity profile was shared by both liquid and gas phases, a new function was written to assign a non-value to the gas phase. The selection of the threshold of 0.5 was selected to separate the gas liquid interphase in the VOF (Ratkovich et al., 2013)28. It is important to highlight, that the VOF model uses the transport of a volumetric fraction of each phase to obtain the phase distribution; the value is 0 (only gas) or 1 (only liquid) in the interrogation windows of the domain. However, the volumetric fraction can also take values between 0 and 1 owing to the mesh refinement. Therefore, to decide which region is occupied by which phase, a threshold of 0.5 was used. This means that the volumetric fraction below 0.5 was taken as liquid and above 0.5 was taken as gas. The translational velocity, vT, was obtained at 15 m and 14.73 m from the inlet. It was calculated using cross correlation approaches previously applied by Gokcal (2008) 12 and Brito (2012)5. They described the cross correlation function as a measurement of how the two signals correlate with each other as a function of the time displacement between them. If signals are identical, the cross correlation will be one and if they are completely dissimilar, the cross correlation will be zero. Then, vT is calculated as below.

vT 

LCS i  CS j



(8)

.

CC

where, ΔLCSi→CSj is the distance between capacitance sensor one and two in meters, τ is the temporal lag in seconds.

Tested Conditions

A

The range of superficial velocities that were compared to the experimental data can be found in Table 3. The effects of the liquid viscosity (µL) and the superficial gas velocity (vSg) on the flow fields and slug characteristics were studied in different cases. The experiments and simulations were carried out with a constant superficial liquid velocity (vSL) of 0.2 m/s. Table 3. Different cases to validate the model Case number

vSL [m/s]

vSg [m/s]

µL [mPa·s]

1 2 3

0.2 0.2 0.2

0.2 0.2 0.7

184 567 161

11

4 5 6

0.2 0.2 0.2

0.7 1.0 1.0

567 168 567

Different levels of the CFD parameters have been modified to summarize the variations in the tested models (Table 4). All the combinations of different levels were tested with the case 1, and the best combination was tested with the other cases.

Parameters and model

Level 1

Sharpening factor Angle factor IMD

0 0.05 Implemented

SC RI PT

Table 4. Summary of tested CFD model parameters.

Level 2

1 0.2 Not implemented

Extrapolated Conditions

A

N

U

After the model was validated, the extrapolated conditions were simulated by the developed model. The changes in the translational velocity and velocity profile with respect to the liquid viscosity and the superficial gas velocity were observed. Since the validation was performed with the superficial liquid velocity of 0.2 m/s, the same value was used for the extrapolated conditions. A total of 12 points were extrapolated and simulated (Table 5). The translational velocity values were obtained experimentally for cases up to 3 m/s of superficial gas velocity to have at least one comparable experimental parameter. Table 5. Levels of variables for extrapolated conditions.

D

TE

0.2 0.2 0.2 0.2 0.2 0.2

vSg [m/s]

µL [mPa·s]

1.3 1.6 2.0 3.0 4.0 5.0

180 and 570 180 and 570 180 and 570 180 and 570 180 and 570 180 and 570

M

vSL [m/s]

EP

Error Calculation

CC

Experimental data were compared to CFD numerical results by means of the absolute relative error:

 N x numi  xexpi 1 error (%)      N i 1  xexpi  

   100  

(9)

A

where x is the variable and N is the number of data points. The sub index num and exp stands for numerical results and experimental data, respectively.

Results and Discussions The results of this study can be divided into three sections. The first is the selection of the numerical parameters of the CFD model. The formation of the slug flow is described qualitatively and quantitatively based on the sensitivity analysis of parameters such as the sharpening factor and the angle factor. Secondly, the comparisons between the CFD and experimental results are presented,

12

validating the applicability of the developed numerical model. And last, the simulated results from the extrapolated velocity conditions are visualized.

Sensitivity Analysis and Selection of the Best Numerical Parameters

N

U

SC RI PT

The slug frequency and average liquid holdup were utilized to compare the experiments and the simulations for the selection of the best combination. It should be noticed again that the sensitivity analysis of CFD parameters was performed for the first experimental case (vSL = 0.2 m/s, vSg = 0.2 m/s and µL = 184 mPa·s; Table 3). Then, the selected model was applied to simulate remain cases. Before selecting the best model, an analysis of the effect of each CFD parameter was performed. The firstly tested parameter was the sharpening factor. The case with the sharpening factor value of 1.0 presented a higher penetration of the elongated gas bubble in the slug tail (Figure 9). The IMD model and an angle factor of 0.05 were used, but the effects of them on the combinations of other parameters were the same. The penetration allowed a better definition of slug flow, since the dispersed phase (gas phase) could break out and generate separated bubbles or elongated bubbles.

A

Figure 9. The effects of sharpening factor on the phase distribution (liquid holdup) at the slug tail.

A

CC

EP

TE

D

M

Figure 10 shows an example of the effects of the angle factor on the results in the slug front. It can be observed that a lower value of angle factor presented a smoother interface. The slug front was more perpendicular to the flow direction with a higher value of angle factor. In addition, a wave lip was presented by using an angle factor of 0.2. This characteristic was desired to compare with the experimental result to decide which value is adequate to use. In other words, a direct comparison with experimental results was required. Another effect of the angle factor was the curvature of the interface in the liquid film. In the case with an angle factor of 0.05 the interface in the liquid film was smooth and aligned with the mesh, while a higher value of angle factor presented a curvature in the interface.

Figure 10. The effects of angle factor on the phase distribution (liquid holdup) at the slug front.

The effects of the IMD model is presented in Figure 11. Two values of angle factor were used to have a better visualization of the effects of IMD model. It was observed that the implementation of the IMD model represented a smoother interface in the slug front. The longer wave lips could be observed without applying the IMD model, while the smooth slug front or shorter wave lips could be presented by using the IMD model. This might indicate that, in this case, the parasitic currents were larger at higher values of angle factor or that the artificial momentum dissipation applied in the IMD 13

U

SC RI PT

model was softening the phenomenon. Although the angle factor played a role in the IMD performance, the cases where wave occurred could have larger parasitic currents since the wave tube was scattered by the discretization.

N

Figure 11. The combined effects of the IMD model and angle factor on the phase distribution (liquid holdup) at the slug front.

D

M

A

A selection of the best model was performed by considering the effect of different CFD parameters analyzed above. The measured slug frequency was 0.565 slugs per second and the average liquid holdup was 0.677 the first experimental case (vSL = 0.2 m/s, vSg = 0.2 m/s and µL = 184 mPa·s; Table 3). Table 6 and Table 7 show the summaries of the slug frequency and the average liquid holdup, respectively, resulted from the simulations. For both slug characteristics, the best results were obtained with a sharpening factor of 1.0 and an angle factor of 0.05 using the IMD model. An angle factor of 0.2 showed the secondary best agreement with the experimental results.

TE

Table 6. Summary of slug frequency (slugs/s) results from the CFD simulations. With IMD model

EP

0.05 0.2

Angle factor

Without IMD model

Sharpening factor 0.0 1.0 0.2 0.6 0.35 0.6

Angle factor

0.05 0.2

Sharpening factor 0.0 1.0 0.3 0.9 0.3 0.5

A

CC

Table 7. Summary of the average liquid holdup results from the CFD simulations

Angle factor

With IMD model

0.05 0.2

Without IMD model

Sharpening factor 0.0 1.0 0.755 0.723 0.794 0.750

Angle factor

0.05 0.2

Sharpening factor 0.0 1.0 0.748 0.750 0.800 0.766

Comparison Between Experimental Measurements and the CFD Simulations Figure 12 describes the comparison of the velocity profile for cases 1 and 2. The experimental results show that a longer distance might be required to develop a quasi-parabolic velocity profile in lower viscosity case (case 1; 184 mPa·s). The same trend is presented in CFD results.

14

SC RI PT

The gas entrainment is presented experimentally as a blue zone, perturbing (or penetrating) the upper part of the pipe in higher viscosity case (case 2; 567 mPa·s). This is also presented in CFD results as a white zone in the same part of the pipe. However, this perturbation (or penetration) is not presented at a lower viscosity.

U

Figure 12. Velocity profiles at the slug front (cases 1 and 2).

A

CC

EP

TE

D

M

A

N

Figure 13 presents the slug body for cases 1 and 2. The amount of entrained gas in slug body increases as viscosity increases. The entrained gas is shown as bubbles in the slug body. These bubbles generate the slight fluctuations in the velocity profile at higher viscosity condition (case 2; 567 mPa·s). In this condition, the bubbles also generate more perturbations in the upper part of the velocity profiles. This can imply that gas phase may not infiltrate through the top of slug body at slug front at lower viscosity case (case 1; 184 mPa·s). On the other hand, the upward flow-stream or vorticities may lead the entrained bubbles to near the top of the pipe in slug bodies. In the CFD results, the entrained bubbles are presented as white zones in the upper part of the pipe, and a higher concentration of them is presented for higher viscosity condition (case 2; 567 mPa·s). The results indicate that the CFD simulation is possible to predict the regions where the bubbles are presented.

Figure 13. Velocity profiles in the slug body (cases 1 and 2).

Figure 14 shows the velocity profiles at the end of slug body (slug tail). A higher viscosity condition presents a lower liquid film height near the slug tail. Furthermore, the elongated bubble nose becomes sharper as viscosity increases, penetrating more length from the end of slug body. This behavior is also visualized by CFD model.

15

SC RI PT

Figure 14. Velocity profiles at the slug tail (cases 1 and 2).

A

CC

EP

TE

D

M

A

N

U

Since the entrained gas bubbles occupy more space in the slug body as viscosity increases, the liquid holdup in the slug body increases as liquid viscosity decreases. This statement agrees with Abdul-Majeed (2000)1 who proposed an inverse relationship between the slug liquid holdup and the liquid viscosity. Nevertheless, his study was performed with liquid viscosity below 6.7 mPa·s. Figure 15 illustrates the velocity profiles in the slug body for cases 3 and 4 (vSL = 0.2 m/s and vSg = 0.7 m/s). The experimental results show a larger gas entrance in the slug body as liquid viscosity increases. This behavior also can be resulted from the CFD model. The parabolic profiles still exist, while the maximum velocity is located at the upper part of the pipe. The velocity profiles show a slight convective acceleration in the vertical direction, especially near the pipe center.

Figure 15. Velocity profiles in the slug body (cases 3 and 4).

The higher superficial gas velocity that the PIV technique was able to analyze (with a constant superficial liquid velocity of 0.2 m/s) was 1.0 m/s. Figure 16 shows the velocity profile in for this condition (vSL = 0.2 m/s and vSg = 1.0 m/s) with 168 mPa·s and 567 mPa·s (cases 5 and 6, respectively). The target area at lower viscosity condition is in the slug body but near the slug front (case 5; 168 mPa·s). It presents the asymmetrical velocity profile with the higher velocities at the upper parts of the pipe. This trend also can be captured by CFD model. 16

A

N

U

SC RI PT

The target area at higher viscosity condition is in slug body but closer to the middle of it (case 6; 567 mPa·s). A relatively shorter slug length at higher viscosity condition resulted in the slightly asymmetrical velocity profiles in this region. In addition, the increased amount of entrained bubbles can be captured by CFD model.

Figure 16. Velocity profiles in the slug body (cases 5 and 6).

TE

D

M

A quantitative comparison between experimental and CFD results was performed. Three slug characteristics were compared, namely translational velocity (vT), average liquid holdup (HL Avg), and slug frequency (fS). Figure 17 shows the translational velocity as a function of the mixture velocity for the six cases listed in Table 3. It can be shown in both groups (μL ≈ 170 mPa·s (average) and 567 mPa·s) that the translational velocity simulated by CFD model presents a good agreement with the experimental data. The average relative error was 10 % for the cases with 170 mPa·s and 11.1% for the cases with 567 mPa·s.

1.5

2.5

vT, (m/s)

2

CC

vT, (m/s)

3

EP

3 2.5

1

0.5

A

1 0.5

(a) 170 mPa·s

(b) 567 mPa·s

0

0

0.4

2 1.5

0.6

0.8

1

0.4

1.2

0.6

CFD

0.8

1

1.2

vM, (m/s)

vM, (m/s) CFD

Experimental

Experimental

Figure 17. Measured translational velocity (PIV) vs. Simulated translational velocity (CFD).

Figure 18 presents a comparison of the average liquid holdup results between the experiments and CFD model. The CFD results give the suitable predictions compared to the experimental results. The average relative error values are 7.5% and 5.5% for the cases of 170 mPa·s and 567 mPa·s, respectively. 17

0.8

0.8

(b) 567 mPa·s

HL Avg, (-)

HL Avg, (-)

(a) 170 mPa·s 0.7 0.6 0.5

0.7 0.6 0.5

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

vSg, (m/s)

CFD

SC RI PT

vSg, (m/s) Experimental

CFD

Experimental

Figure 18. Measured average liquid holdup (CP) vs. Simulated average liquid holdup (CFD).

TE

D

M

A

N

U

The slug frequency as a function of the superficial gas velocity for the considered cases can be seen in Figure 19. The slug frequency slightly decreases as the superficial gas velocity increases with lower viscosity (170 mPa·s). This also can be resulted from CFD model. In general, the increase of superficial gas velocity leads to the increase of interfacial shear stress reducing the change of waves bridging the top of the pipe. On the other hand, the experimental results from the liquid viscosity of 567 mPa·s do not illustrate a clear trend. The fluctuations may be caused by the experimental method applied to obtain the slug frequency. As commented by Kim et al. (2018)18, a slug unit is converted from a capacitance sensor voltage to a vector of zeros and ones, defining a threshold value to characterize film and slug regions. Then, the slug frequency can be calculated by the number of slug bodies per data recording time. At higher viscosity, however, the increased amount and size of entrained bubbles in liquid slug body may lead to the unclear separation of film and slug regions. Therefore, the experimental results may count more slug bodies than real number of them, increasing the slug frequency. Similarly, the CFD simulation results in some fluctuations of the simulated slug frequency as increase in the superficial gas velocity at higher viscosity condition. Because the CFD model looks into the holdup fluctuation as a function of time and then counts the number of peaks in a certain time. As a result, the CFD simulation presents the average relative errors of 11.4% and 28% for low and high viscosity cases, respectively.

0.5 0.3

CC

0.2

0.4

fS, (slug/s)

0.9 0.7

1.1

(a) 170 mPa·s

EP

fS, (slug/s)

1.1

0.9 0.7 0.5

(b) 567 mPa·s

0.3 0.6

0.8

1

0.2

0.4

0.8

1

vSg, (m/s)

vSg, (m/s) CFD

0.6

Experimental

CFD

Experimental

A

Figure 19. Measured slug frequency (CP) vs. Simulated slug frequency (CFD).

The CFD simulation, moreover, can estimate the wall shear stress profile in the slug unit. Figure 20 shows the wall shear stress at the bottom of the pipe as a function of the distance from the inlet. It is observed that the wall shear stress repeats the cycles between 0 Pa and 12 Pa. The wall shear stress never has a plateau at the bottom of the pipe, continuously increasing in the liquid slug body. This phenomenon coincides with the experimental reports by Buitrago et al. (2017)6. They reported that the fully developed flow conditions are not reached at the bottom of the pipe in liquid slug body. This is because of a lower momentum region inside the liquid phase beneath the mixing zone and close to the lower sections of the pipe at slug front (Kim et al., 2018)18. 18

The wall shear stress decreases again from the end of slug body toward the upcoming slug front, namely in the liquid film region. Sharma et al. (1998)30 also reported this phenomenon. They observed a very rapid deceleration in the liquid layer from the end of slug body. They mentioned that the effect of interfacial shear exerted by the faster-moving gas bubble seems to be confined to the very vicinity of the interface.

10 5 0 10

11

12

13

14

SC RI PT

τW, (Pa)

15

15

16

17

Distance from inlet, (m)

U

Figure 20. Wall shear stress profile as a function of distance from the inlet (case 1; vSL = 0.2 m/s, vSg = 0.2 m/s and μL = 184 mPa·s).

N

Results from Extrapolated Conditions; Higher Velocities

Slug front

EP

Slug tail

CC

vSg (m/s)

TE

D

M

A

The chosen CFD model that shows the best performance was utilized to obtain the velocity profiles for the cases where the PIV technique was not available. Since there is no inner flow field information, a comparison of the visualized phase distribution was carried out. Figure 21 shows the phase distribution recorded by high-speed camera and the velocity profile simulated by CFD with the liquid viscosity of 180 mPa·s. A Photron FASTCAM SA3 60K-C2 high-speed camera with a Nikon AF-S VR Micro-Nikkor 105 mm f/2.8G IF-ED lense have been used. It is experimentally observed that the liquid slug completely fills the pipe cross section below 2 m/s of superficial gas velocities (with vSL = 0.2 m/s). Above that value, the gas passes in the upper part of the slug making the gaps near the pipe top. A similar trend is presented by CFD model. In addition, the sharp interfaces are also illustrated in CFD results. The velocity profiles show a higher velocity at the upper part of the pipe and in the sharp interfaces.

A

1.3

1.6

19

SC RI PT N

U

2.0

M

A

3.0

D

Figure 21. Phase distributions from high-speed camera and CFD model at slug front and tail (μL=180 mPa·s).

Slug tail

Slug front

A

CC

vSg (m/s)

EP

TE

Figure 22 shows the same comparison but with the liquid viscosity of 570 mPa·s. It is experimentally shown that a higher amount of bubbles entrained in the slugs. The liquid slug only fills the entire pipe when the superficial gas velocity is lower than 1.6 m/s. Similar to the results in Figure 21, the higher velocities are presented at the upper part of the pipe. The sharp interfaces are also presented.

1.3

20

SC RI PT

1.6

M

A

N

U

2.0

TE

D

3.0

EP

Figure 22. Phase distributions from high-speed camera and CFD model at slug front and tail (μL= 570 mPa·s).

A

CC

In Figure 21 and 22, the velocity profiles simulated by CFD model at slug front and tail seem to represent the stratified wavy flow. Nevertheless, the flow pattern was slug flow as exemplified in a three-dimensional view in Figure 23. It can be seen how the liquid is bridged to the pipe walls and also the upper part of the pipe is not entirely occupied by the liquid slug. Furthermore, Figure 23 shows the streamlines of the gas path. Those streamlines were generated in the long bubble behind the liquid slug, illustrating how the gas streamlines go across the slug. It indicates that the gas is blowing in the upper part of the slug body. Then, the liquid does not generate a hydraulic sealing between the tail and the front of the slug in relatively high gas velocity conditions.

21

SC RI PT

Figure 23. Velocity profile (in m/s) at the interface for the case with 5 m/s of superficial gas velocity with vSL = 0.2 m/s and μL = 570 mPa·s.

EP

TE

D

M

A

N

U

Figure 24 shows a three-dimensional representation of a slug velocity profile of case 1 (vSL = 0.2 m/s, vSg = 0.2 m/s and μL= 180 mPa·s). A smoother interface at the slug front is observed compared to higher superficial gas velocity condition in Figure 23. The gas streamlines were generated following the same way as performed in Figure 23. In contrast to the previous case, there is no gas passing through the upper part of slug body from the tail to the front.

CC

Figure 24. Velocity profile (in m/s) at the interface for a case with 0.2 m/s of superficial gas velocity with vSL = 0.2 m/s and μL = 180 mPa·s.

A

Figure 25 presents the translational velocity obtained by experiments and simulated by CFD model. An absolute average relative error shows 10% in the case of μL = 170 mPa·s and 14% in the case of μL = 570 mPa·s. The linear trends are observed up to 3.2 m/s and 2.3 m/s of mixture velocity with lower and higher viscosity cases, respectively. The similar trends can be resulted from CFD model. Above the certain level of mixture velocities, the simulated translational velocity from CFD model shows a decreasing slope. Moreover, even lower values of translational velocity are obtained as the mixture velocity increases. This reduction in the translational velocity can be explained by the incomplete hydraulic liquid sealing made by the gas path-line as shown before. This reduction in the translational velocity was reported by Gokcal et al. (2008)13 and Brito (2012)5. However, it is contrary to the model proposed by Nicklin et al. (1962)24 where the translational velocity presents a constantly linear trend with the mixture velocity. As a result, the gas passing through near the top of slug body and the continuous increase in wall shear stress profile near the bottom in the slug body are opposed to the assumptions of existing mechanistic models. 22

8

6

6

4 2

SC RI PT

8

vT , (m/s)

vT , (m/s)

The correlation suggested by Fabre (1994)10 was compared as can be seen in Figure 25. It results in the linear trend between the translational velocity and mixture velocity, while showing some discrepancies with the experimental data. The absolute average relative error between them is 34% in the case of μL = 170 mPa·s and 38% in the case of μL = 570 mPa·s, respectively. These error values are higher than the ones calculated by CFD model.

4 2

(b) 570 mPa·s

(a) 170 mPa·s 0

0 0.4

1.4

2.4

3.4

0.4

4.4

1.4

vM , (m/s) CFD

2.4

3.4

4.4

vM , (m/s)

Experimental

Fabre

CFD

Experimental

Fabre

Figure 25. Translational velocity vs. Mixture velocity for (a) μL = 170 mPa·s, (b) μL = 570 mPa·s.

Conclusions

A

CC

EP

TE

D

M

A

N

U

The selection of the best combination of CFD parameters to simulate air and highly viscous liquid two-phase slug flow in horizontal pipes was performed. The results present a fair correspondence with experimental results in terms of the average liquid holdup, slug frequency, translational velocity, and velocity profile. The effect of the sharpening factor in the VOF model was demonstrated, showing the penetrating bubble into the slug body from the end of it. In addition, the different angle factors in the HRIC scheme were analyzed, resulting in the dissimilar roughness of the interface. The velocity profiles in the slug body were simulated by the developed CFD model. They presented the suitable agreements with the experimental results based on both magnitude and behavior aspects. In addition, the velocity fields showed the regions where the gas bubbles are entrained in the slug body. The wall shear stress profiles at the bottom of the pipe were also calculated by CFD model. A continuous increase of wall shear stress from slug front to end was observed, while velocity profiles showed almost symmetrical parabolic shapes from certain points in slug bodies. This indicates that the liquid slug body is not fully developed at least near the pipe bottom area. This phenomenon coincides with the experimental observations of Kim et al. (2018)18 and Buitrago et al. (2017)6. The range of operational conditions was extrapolated to estimate the flow behaviors at higher velocity conditions. The derivative of the slug translational velocity with respect to the mixture velocity decreased above 2~3 m/s of mixture velocity. Furthermore, the reduction of translational velocity was shown by CFD model above 4.5 m/s. It was demonstrated that the decreased ratio (vT/vM) and reduced translational velocity were owing to the gas passing through the top of slug body. In other words, the liquid made an incomplete hydraulic sealing. Nevertheless, most of the existing mechanistic models do not catch this phenomenon and assumes a linear trend between translational velocity and mixture velocity. It is suggested to review and modify this assumption to avoid the inconsistencies with realistic evidence. CFD technique proved itself as a reliable tool in this study. It confirmed a lot of advantages against experiments such as complete information of the velocity and wall shear stress profiles inside the whole pipe. It presented the space-based results instead of the time-based ones.

23

Acknowledgement We gratefully acknowledge the support of TUFFP for the experimental study and the informative discussions. References

A

CC

EP

TE

D

M

A

N

U

SC RI PT

1. Abdul-Majeed, G. H., 2000. Liquid slug holdup in horizontal and slightly inclined two-phase slug flow. J. Pet. Sci. Eng. 27 (1-2), 27-32. 2. Al-Safran, E., 2009. Prediction of slug liquid holdup in horizontal pipes, J. Energy Resour. Technol. 131 (2), 023001 (1-8). 3. Al-Safran, E., Kora, C., Sarica, C., 2015. Prediction of slug liquid holdup in high viscosity liquid and gas two-phase flow in horizontal pipes. J. Pet. Sci. Eng. 133, 566-575. 4. Andersen, E., 2012. Multiple holdup solutions and the effect of interface level gradients. M.S. Thesis. Institutt for energi- og prosessteknikk, Norwegian University of Science and Technology, Trondheim, Norway. 5. Brito, R., 2012. Effect of medium oil viscosity on two-phase oil-gas flow behavior in horizontal pipes. M.S. Thesis. The University of Tulsa, Tulsa, Oklahoma, USA. 6. Buitrago, L., Kim, T., Pereyra, E., Sarica, C., 2017. Wall shear stress measurements of horizontal two-phase slug flow for high viscosity liquids using constant temperature anemometry. BHR Group 18th International Conference on Multiphase Technology, Cannes, France, June 7-9. 7. Colmenares, J., Ortega, P., Padrino, J., Trallero, J. L., 2001. Slug flow model for the prediction of pressure drop for high viscosity oils in a horizontal pipeline. SPE International Thermal Operations and Heavy Oil Symposium, Porlamar, Margarita Island, Venezuela, March 12-14. 8. De Schepper, S. C. K., Heynderickx, G. J., Marin, G. B., 2008. CFD modeling of all gas-liquid and vapor-liquid flow regimes predicted by the Baker chart. Chem. Eng. Journal 138 (1-3), 349357. 9. Dukler, A. E., Hubbard, M. G., 1975. A model for gas-liquid slug flow in horizontal and near horizontal tubes. Ind. Eng. Chem. Fundam. 14 (4), 337-347. 10. Fabre, J., 1994. Advancements of two-phase slug flow modeling. SPE paper 27961-MS presented at The University Centennial Petroleum Engineering Symposium, Tulsa, OK. 11. Fan, S., Yan, T., Yeung, H., Lao, L., 2013. Velocity characteristics of slug body and film for twophase gas-liquid slug flow using ultrasonic techniques, 8th Int. Conf. on Multiphase Flow ICMF 2013, Jeju, Korea, May 26-31. 12. Gokcal, B., 2008. An experimental and theoretical investigation of slug flow for high oil viscosity in horizontal pipes. Ph.D. Dissertation. The University of Tulsa, Tulsa, Oklahoma, USA. 13. Gokcal, B., Wang, Q., Zhang, H., Sarica, C., 2008. Effects of high oil viscosity on oil/gas flow behavior in horizontal pipes. SPE 102727, SPE Projects, Facilities & Construction vol. 3 (2). 14. Gregory, G. A., Nicholson, M. K., Aziz, K., 1978. Correlation of the liquid volume fraction in the slug for horizontal gas-liquid slug flow. Int. J. Multiph. Flow 4 (1), 33-39. 15. Harvie, D. J. E., Davidson, M. R., Rudman, M., 2006. An analysis of parasitic current generation in Volume of Fluid simulations. Appl. Math. Model. 30 (10), 1056-1066. 16. Hernandez-Perez, V., Abdulkadir, M., Azzopardi, B., 2011. Grid generation issues in the CFD modelling of two-phase flow in a pipe, J. Comput. Multiph. Flows 3 (1), 13-26. 17. Hirt, C. W., Nichols, B. D., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1), 201-225. 18. Kim, T. W., Aydin, T. B., Pereyra, E., Sarica, C., 2018. Detailed flow field measurements and analysis in highly viscous slug flow in horizontal pipes. Int. J. Multiph. Flow 106, 75-94. 19. Lo, S., Tomasello, A., 2010. Recent progress in CFD modelling of multiphase flow in horizontal and near-horizontal pipes. CD-adapco, pp. 209-219.

24

A

CC

EP

TE

D

M

A

N

U

SC RI PT

20. Malan, P., Suluksna, K., Juntasaro, E., 2009. Calibrating the γ-Reθ transition model for commercial CFD. 47th AIAA Aerospace Sciences Meeting Including The New Horizons Forum and Aerospace Exposition, Orlando, Florida, January 5-8. 21. Menter, F. R., 1992. Influence of freestream values on k-ω turbulence model predictions. AIAA Journal, 30 (6), 1657-1659. 22. Menter, F. R., 1992. Improved two-equation k-ω turbulence models for aerodynamic flows. NASA Technical Memorandum, Ames Research Center, Moffett Field, California. 23. Monroy, J., Sato, M., 2003. A VOF numerical model of storm waves on coral reefs. ISOPE-I-03305, 13th International Offshore and Polar Engineering Conference, Honolulu, Hawaii, USA, May 25-30. 24. Nicklin, D. J., Wilkes, J. O., Davidson, J. F., 1962. Two-phase flow in vertical tubes. Inst. Chem. Eng. 40, 61-68. 25. Nogueira, S., Sousa, R. G., Pinto, A. M. F. R., Riethmuller, M. L., Campos, J. B. L. M., 2003. Simultaneous PIV and pulsed shadow technique in slug flow: a solution for optical problems. Exp. Fluids 35, 598-609. 26. Pereyra, E., Torres, C., Mohan, R., Gomez, L., Kouba, G., Shoham, O., 2012. A methodology and database to quantify the confidence level of methods for gas-liquid two-phase flow pattern prediction. Chem. Eng. Res. Des. 90 (4), 507-513. 27. Raffel, M., Willert, C. E., Wereley, S., Kompenhans, J., 2007. Particle Image Velocimetry, Springer-Verlag Berlin Heidelberg. 28. Ratkovich, N., Majumder, S. K., Bentzen, T. R., 2013. 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. 91, 988-998. 29. Reyes-Gutiérrez, M. A., Rojas-Solórzano, L. R., Marín-Moreno, J. C., Meléndez-Ramírez, A. J., Colmenares, J., 2005. Eulerian-Eulerian Modeling of Disperse Two-Phase Flow in a Gas-Liquid Cylindrical Cyclone. J. Fluids Eng. 128 (4), 832-837. 30. Sharma, S., Lewis, S., Kojasoy, G., 1998. Local studies in horizontal gas-liquid slug flow. Nucl. Eng. Des. 184, 305-318. 31. Taitel. Y., Barnea, D., 1990. A consistent approach for calculating pressure drop in inclined slug flow. Chem. Eng. Sci. 45 (5), 1199-1206. 32. 2015 CD-adapcoTM, 2015. “STAR-CCM+ v.10.02.010”. 33. Waclawczyk, T., Koronowicz, T., 2008. Comparison of CICSAM and HRIC High-Resolution Schemes for Interface Capturing. J. Theor. Appl. Mech. 46 (2), 325-345. 34. Yang, H., Lu, J., Yan, S., 2014. Preliminary numerical study on oil spilling from a DHT. ISOPEI-14-138, 24th International Ocean and Polar Engineering Conference, Busan, Korea, June 1520. 35. Yuan, L., 2006. Numerical study of round vertical jets in cross flow with VOF method. ISBN 1880653-67-2, 7th ISOPE Pacific/Asia Offshore Mechanics Symposium, Dalian, China, September 17-21. 36. Zhang, H., Wang, Q., Sarica, C., Brill, J., 2003. A unified mechanistic model for slug liquid holdup and transition between slug and dispersed bubble flows. Int. J. Multiph. Flow 29 (1), 97107.

25