Numerical study of hydrodynamic characteristics of gas liquid slug flow in vertical

Multiphase lows occur in wide applications including; nuclear, chemical, and petroleum industries. One of the most important low regime encountered in multiphase low is the slug low which is often encountered in oil and gas production systems. The slugging problems may cause looding of downstream processing facilities, severe pipe corrosion and the structural instability of pipeline and further induce the reservoir low oscillations, and a poor reservoir management. In the present study, computational luid dynamics simulation is used to investigate two phase slug low in vertical pipe using the volume of luid (VOF) methodology implemented in the commercial code ANSYS Fluent. The viscous, inertial, and interfacial forces have signiicant effect on the hydrodynamic characteristics of two-phase slug low. These forces can have investigated by introducing a set of dimensionless numbers, namely; inverse viscosity number, N f , Eotvos number, E o , and Froude number, Fr TB . The simulation accounts for the hydrodynamic features of two phase slug low including; the shape of Taylor bubble, bubble proile, velocity and thickness of the falling ilm, wake low pattern, and wall shear stress distribution. The CFD simulation results are in good agreement with previous experimental data and models available in literature.


Introductıon
Multiphase flow is a flow field in which more than one phase exists.It could be two-phase or three-phase or even four-phase flow depending on the number of phases encountered in the flow domain.Different combinations between the phases such as; gas bubbles in liquid, liquid droplets in gas, and solid particles in either gas or liquid are all different types of multiphase flow.Multiphase flow is complex in nature compared to single phase flow.The difficulty of multiphase flow arises from; the need of good knowledge of the chemical, physical and thermodynamics properties of each phase, the variation in the fluid composition with time as it flows in multiple locations, slugging problem, corrosion and erosion, and the variation in pressure and temperature leading to deposit formation.
Multiphase flow occurs in a wide range of applications including; natural processes as rivers and clouds formation, nuclear systems as nuclear waste processing, fluidized beds, boiling and condensation in nuclear reactors, chemical processes as food mineral processing and transportation, transport of cement, grains, and metal powders and petroleum industries as oil processing, oil and gas transport in pipelines, sloshing in offshore separator devices.
Multiphase flow is classified according to the distribution of the phases, which is known as flow regime, or flow pattern.Identifying the encountered flow regime in multiphase flow field is an important aspect prior performing any investigation on the flow.For gas-liquid flow in pipes, one of the common and complex pattern encountered in multiphase flow is known as slug flow.Slug flow is an intermitted flow between stratified and annular flow.It is characterized by an elongated bullet shaped bubble, known as Taylor bubble, and a liquid slug, which is a liquid film flowing downwards between the bubble interface and pipe wall.
Computational fluid dynamics (CFD) has been proven to be a powerful practical tool for the analysis and simulation of the hydrodynamic characteristics of slug flow in pipes.The main complex feature of gas-liquid slug flow is the deformable interface [1].The volume of fluid (VOF) method originally developed by Hirt CW [2] is often used to simulate complex multiphase flows including slug flow, and is powerful in tracking the interface between fluids [3][4][5][6][7].
The hydrodynamic characteristics of gas-liquid vertical slug flow is governed by three main forces namely; viscous, inertial, and interfacial forces.A dimensionless analysis for two-phase gas-liquid slug flow can be done by applying Pi-Buckingham theorem.This simplifies the problem to three dimensionless groups, which are; Eötvös number, , Froude number,   , and inverse viscosity number,   , defined as follows: Eötvös number, , is the ratio between gravitational forces, and surface tension forces, and given by;  =  *  *  2 ⁄.
Where; ,   ,  , , , and   are the density of the liquid phase (kg/m 3 ), dynamic viscosity of the liquid phase (Pa.s), gravitational acceleration (m/s 2 ), pipe diameter (m), surface tension coefficient Froude number,   , and inverse viscosity number.A validation of the present numerical code is done by performing comparison with both experimental data, and with correlations available in literature.

Model Development
Following the concept of slug unit cell, two-phase vertical slug flow can be characterized by an elongated bullet shaped bubble, known as Taylor bubble, and a liquid slug which is a liquid film flowing downwards between the bubble interface and pipe wall.The computational time is an important issue that need to be considered while performing any multiphase flow simulation.Thus, in order to reduce the computational time, the present simulation is performed in 2D coordinate system, with axial symmetry around the centerline of the pipe.The simulation domain is a vertical pipe of diameter, D, and length, L. The flow domain is constructed and solved using the VOF methodology implemented in the computational dynamic software package, ANSYS Fluent (Release 15.0).In all simulated cases, a uniform grid of quadrilateral control elements is applied.Different grids were tested to check solution convergence.The present simulation has been performed for 2D coordinate system, unsteady, and laminar flow with constant fluid properties.The two phases were assumed as incompressible, viscous, immiscible, and not penetrating each other.

Governing equations
In the present model, the fluids share a well-defined interface and hence the volume of fluid (VOF) method for two phase flow has been selected.The VOF model is a surface-tracking technique applied to a fixed Eulerian mesh.This model is designed for two or more immiscible fluids to track the interface between them.This model solves a single set of momentum equation that is shared by the two fluids, and the volume fraction of each of the fluids in each computational cell is followed throughout the domain.Details of the governing equations and the treatment of the interface can be obtained from the Fluent user guide [12].The equations being solved in the VOF model models are as given in this section.
One set of continuity and momentum equations are solved for the two-phase system.Firstly, the Reynolds average continuity equation in a VOF model for n number of phases can be expressed as follows; ( ) Where; , , ,  and   and are the volume-fraction-averaged density, time, time average velocity vector of the flowing fluid, number of the phases (for the present two phase flow = 2), and mass source respectively.The mass source,   , is set to zero in the present case.In addition, a single Reynolds average momentum equation is solved throughout the computational domain, and all phases share the same resulting velocity field.The general momentum equation can be written as follows; Where; , , and  and are the pressure in the domain, volumefraction-averaged viscosity of the flowing fluid, and the additional forces as gravitational term and surface tension, respectively.In other words, the left hand side of equation ( 2) represents the unsteady term and convection terms.While, the right hand side represents the pressure term, the diffusion term, the body force and other external forces that might act on the system.The continuum surface force (CSF) of Brackbill [13] is used to account for the surface tension effects.
(N/m), and Taylor bubble rise velocity (m/s), respectively [8].The main hydrodynamic characteristics of gas-liquid vertical slug flow are; the shape of the Taylor bubble, Taylor bubble rise velocity, liquid film velocity, liquid film thickness, and wall shear stress distribution.Despite the fact that hard effort has been done in the modelling process of gasliquid slug flow, a needs for closure correlations based on experimental data is still required.These relationships included slug characteristics such as; slug frequency, slug length, slug liquid hold up, and slug unit velocity.This is in addition to other relationships as shear force, wall friction, and gas entrainment rate.
In literature, many experimental and numerical studies were done on the rise of single Taylor bubble in both stagnant and flowing fluid in vertical pipes.Taha T [9] investigated the motion of a single Taylor bubble in vertical pipes using the volume of fluid (VOF) method implemented in the CFD software FLUENT.The simulation accounted for the Taylor bubble velocity, bubble shape, wake flow pattern, and wall shear stress distribution for the rise of single Taylor bubble in both stagnant and flowing liquid including both laminar and turbulent flow conditions.Their study involved wide range of Eötvös number, Froude number, and Morton number, and the results showed that the bubble shape is dependent on the liquid viscosity, the surface tension, and independent on bubble length.In addition, the results showed that the inverse viscosity number,   , greatly influence the wake structure, and increases the wall shear stress, and finally reduce the liquid film thickness.Similarly, [1] investigated the hydrodynamic characteristics of single Taylor bubble rising in stagnant liquid using the VOF method.The study further included the near wall treatment, and hence accounting for the near wall mass transfer and wall shear stress.The results revealed that the key parameters governing the twophase flow induced corrosion are; the mass transfer coefficient and wall shear stress.Furthermore, (8) simulated the rise of single Taylor bubble through stagnant Newtonian liquid using the VOF method, with particular focus to the liquid wake structure, and analyzing the results in terms of Taylor bubble velocity, flow around the bubble nose, liquid film and wake region.The study accounted for laminar flow regime only.The authors developed a map gathering all information reached while studying the wake structure.This maps could show the existence or absence of wake structure with the knowledge of Morton number, Eotvos numbers, and the type of concavity of the bubble rear.Moreover, Yan K [10] investigated the hydrodynamic characteristics of single Taylor bubble rising in stagnant liquid with further consideration of the small dispersed bubble in the liquid slug zone.The study accounted for the effect of small dispersed gas bubbles in liquid slug zone on the flow hydrodynamics features and CO 2 corrosion rate.It was concluded that the small dispersed gas bubbles result in higher fluctuations in the liquid slug zone, which subsequently increase the mass transfer and wall shear stress.Lastly, Araújo J [11] investigated the rising of two consecutive Taylor bubbles through vertical stagnant Newtonian liquids under laminar regime using the VOF method.The results accounted for bubble-bubble interaction, and showed the dependency of the wake on the separation distance between the bubbles.
The main aim of the present study is to investigate the hydrodynamic characteristics of gas-liquid slug flow in vertical pipe by applying computational fluid dynamics CFD simulation using the volume of fluid (VOF) methodology implemented in the commercial software ANSYS Fluent.The flow consists of single Taylor bubble rising in stagnant fluid.The study includes the investigation of some of the main hydrodynamic features of slug flow, which are; Taylor bubble shape, Taylor bubble rise velocity, velocity fields, liquid film thickness, and liquid film velocity for a range of Eötvös number, , The VOF formulation relies on the fact that two or more fluids (or phases) are not interpenetrating.For each additional phase added to the model, a variable is introduced which is the volume fraction of the phase in the computational cell.In each control volume, the volume fractions of all phases sum to unity.If the volume fraction of the q th fluid in a cell is given by , thus the following relationship is valid for each computational cell; The fields for all variables and properties are shared by the phases and represent volume-averaged values, as long as the volume fraction of each of the phases is known at each location.Thus the variables and properties in any given cell are either purely representative of one of the phases, or representative of a mixture of the phases, depending upon the volume fraction values [12].Hence, there are three possible conditions: 1.If  = 0; the cell is empty of the  h fluid.The Reynolds average continuity equation ( 1) and the momentum equation ( 2) are thus dependent on the volume fractions of all phases through the volume-fraction-averaged properties;  and .Hence, depending on the local value of as discussed above, the volumefraction-averaged density and viscosity are calculated as follows; Where; the subscripts, ,  and  refer to phase, liquid phase (primary phase) and the gas phase (secondary phase).Tracking the interface between the two phases is achieved by the treatment of the volume fraction of the  h fluid,  through solving a separate continuity equation, given by the AFU Guide [12]; .
Where; , °p q and °  are the mass source term, mass transfer from phase  to phase  and the mass transfer from phase to phase and the mass transfer from phase  to phase , respectively.According to ANSYS Fluent user guide (12) the source term on the right-hand side of equation ( 6) is by default set to zero.The volume fraction equation will only be used to solve the volume fraction of the  h fluid, , and not the primary phase (liquid phase).The gas phase is computed according to the constraint given in equation (3).

Model geometry and boundary conditions
As mentioned earlier, the solution domain is a vertical pipe with diameter, D, and length, L, with symmetry along the centerline.The length of domain is eleven times larger than pipe diameter.Figure 1 shows the boundary conditions and the initial Taylor bubble shape used in the simulation.The initial bubble shape is a cylinder connected to a hemisphere with the same radius giving an overall bullet shape of Taylor bubble.The length and radius of the Taylor bubble are given by;   , and   respectively.This initial shape is simulated until a steady bubble shape is reached.Different bubble shapes were tested and final steady shape of bubble was similar but this only affects the solution convergence.
The simulation is performed by attaching a reference frame to the rising Taylor bubble.Enabling moving reference frame (MRF) in the simulation, causes the rising Taylor bubble to be stationary and the pipe wall moves downwards with velocity equal to that of the bubble [14].The initial guess of Taylor bubble velocity,   , is estimated according to the general correlation of Wallis GB [15], which is given by, 0 0.01 3.37 0.345 Once the Taylor bubble ceases moving up or down in the axial direction, and hence pseudo-steady solution is reached, the velocity is then adjusted.
The initial guess of the liquid film thickness,   , is estimated using the following equation [16]; Referring to Figure 1, the boundary conditions could be summarized as follow; At the top of domain, the inlet flow boundary condition is applied with liquid entering at average uniform velocity equal to velocity of Taylor bubble,   =   ,   = 0.This uniform profile assumption is due to the investigation of the rise of Taylor bubble in stagnant liquid.
At the bottom of the domain, the outflow boundary condition is applied as the liquid phase is the only phase available.The symmetry boundary condition is applied at the pipe centerline; / = 0,   = 0.At the wall, the no-slip condition is applied with wall moving downwards with the following velocities;   = ,   = 0.
The pressure variation in the gas phase is assumed to be constant.The boundary conditions at the gas-liquid interface are given by; 1.The kinematic condition assuming full slip at the interface is applied; u. n = 0.
2. The dynamic boundary condition, which is also known by stress jump condition can be divided into two separate boundary conditions; the tangential stress balance assuming zero interfacial shear stress along the interface; ˆ( ).= 0, n s τ and the normal stress balance;   +  = constant, where; τ , n , ŝ , ,   and  and are the shear stress, unit normal vector at the interface, unit tangential vector at the interface, surface tension, liquid phase side pressure, and curvature of the interface, respectively.
According to Mao [14], the curvature of the interface is expressed in terms of radii of the curvature of the bubble surface, as follows;

Solution strategy and convergence criterion
As mentioned earlier, a transient simulation is carried out in the present case in order to consider the transient behavior of two phase flow.The simulation was carried out using the explicit VOF model.The PISO pressure-velocity was selected.The spatial discretization scheme used are as follows; Green Gauss Node Based for gradient, PRESTO for pressure, Compressive for volume fraction, Quick scheme for momentum, and first order implicit for transient formulation.The scaled absolute values of the residual of the calculated values of mass, velocity in x, and y directions were monitored including a convergence criterion of 10 -3 for each time step.

Results and Discussion
As mentioned earlier, the aim of the present investigation is to study the hydrodynamics characteristics of single Taylor bubble rising in a stagnant Newtonian liquid using set of dimensionless number, which are; Eötvös number,   , Froude number,   , inverse viscosity number,   , and Morton number, M. Figure 2 shows a schematic of a single Taylor bubble rising in stagnant liquid with some of the main hydrodynamic features of slug flow, including; Taylor bubble shape, Taylor bubble velocity,   , liquid film thickness,   , liquid film velocity,   , and wall shear stress distribution,   .
Figure 3 provide the velocity fields and streamlines of the liquid phase in two phase slug flow in vertical pipe.As shown in the figure, the flow could be divided into three zones a, b, and c namely; Taylor bubble nose zone, falling liquid film zone, and Taylor bubble wake zone (liquid slug zone).In the Taylor bubble nose zone, the Taylor bubble moves up with velocity, due to buoyancy, pushing the liquid sideways where liquid film zone starts to develop.In the falling liquid film zone, the liquid moves downwards with velocity,   , and decreasing liquid film thickness,   .Once a balanced between the gravitational and friction forces is reached, a constant terminal liquid film velocity and thickness is developed.In the Taylor bubble wake zone, the falling liquid film starts to plugs into the liquid slug ending with highly mixing zone in the wake structure of the bubble.

Taylor bubble shape
This section introduces the effect of inverse viscosity number,   , on Taylor bubble shape and its profile.Four cases are simulated according to the experimental work by Campos J [17], and their relevant properties are given in Table 1.In all simulation cases, the gas phase used is air with density =1.225 kg/m 3 , and viscosity =1.7894E-05 Pa.s.The liquid phase is aqueous glycerol solution with varying viscosity range.The simulation domain is a vertical pipe of 19 mm diameter, and 209 mm (11D) length.Figure 4 shows the effect of on Taylor bubble shape in a glycerol solution.The increase in gradually increases the concave shape   of the Taylor bubble bottom surface.The simulated results are in good agreement with the experimental observations of Goldsmith H [18]; in highly viscous flow (viscosity dominated flow) the Taylor bubble has spheroid shape where the top end of bubble is prolate, and the bottom end is oblate.While, in low viscosity flow the flattening or concaving shape of Taylor bubble bottom end is observed.There simulation results are as well in good agreement with the numerical work [8,9,19].
To validate the current simulation results, Figure 5 gives a comparison of the bubble shape profile for different values of inverse viscosity number, (cases 1 and 2) with the results of (9).As indicated in the figure, there is a good agreement between both results.The small deviations in the results are due to the different methods used to estimate the initial guess of the liquid film thickness prior simulation.
Figure 6 shows the effect of inverse viscosity number, on the bubble shape profile for wide range of   .The findings suggest that the increase in increases the bluntness of the bubble nose, and the bubble bottom edge is more flattened.In addition, the increase in decreases the liquid film thickness,   .The results are in good agreement with the work of Zheng et al. [19].

Taylor bubble rise velocity
The Taylor bubble rise velocity,   , is one of the main hydrodynamic features frequently used in the description of two phase slug flow.In literature, many experiments were done to estimate the Taylor bubble rise velocity [15,[20][21][22].Table 2 gives the numerical and experimental results for Taylor bubble velocity for all simulation cases with the deviation between results.All deviations obtained are below 10%, which is an accepted range.

Liquid film
The flow in the liquid film is investigated by studying the liquid film thickness,   , and the liquid film axial velocity,   .Figure 7 shows the effect of inverse viscosity number,   , on the liquid film thickness, for the four cases even in Table 1.At constant value of the dimensionless thickness of the liquid film,   /, decreases with the increase in the dimensionless distance from bubble nose, /.It then remains constant at around / =1.
Increasing the inverse viscosity number,   , decreases the liquid film thickness.On the other hand, the increase in increases the liquid film axial velocity, as shown in Figure 8. Due to gravity, the falling liquid film is accelerated along the Taylor bubble with a decrease in the liquid film thickness.Thus, at film axial velocity,   , increases with the increase in the dimensionless distance from bubble nose, /.In conclusion, once a balanced between the gravitational and friction     forces is reached, a constant liquid film thickness, and velocity is established.The results are in good agreement with the work of Zheng et al. [19].

Wake flow structure
The wake structure is one of the vital hydrodynamic characteristics of vertical slug flow especially in describing the interaction and coalescences between successive Taylor bubbles (8).Prior studding the wake structure, it is important to identify the flow behavior in the liquid film, whether it is laminar or turbulent.This is done by estimating the value of critical Reynolds number,   =   *   * ⁄  , where;   is the average liquid velocity in the liquid film.
To ensure laminar flow regime in the liquid film in the present study, the selected simulation cases are based on the experimental work of Campos et al. [17], where the authors concluded that for laminar flow regime  should be less than 500.
Figure 9 shows the wake flow structure for different value of .In all cases the flow is characterized by closed axisymmetric wake as indicated by Campos et al. [17], and no gas entrainment in the liquid slug zone is noticed.At low values of , the liquid from the liquid film zone expands directly and smoothly in to bubble wake zone, which is noticed by parallel streamlines in case 1. Increasing the values of  leads to the development of circulatory vortex in the bubble wake as liquid plugs into the Taylor bubble bottom.The scale and intensity of the vortex increases with higher values of .The wake flow pattern agrees as well with the work of [8,9,19].

Wall shear stress
If the two phase slug flow problem is involved in heat or mass transfer, then the wall shear stress becomes a main significant hydrodynamic parameter.This processes are often referred to as slug flow induced corrosion [1,8,19,[20][21][22][23][24].The main problems results from slug flow corrosion are; pipeline damage, decrease pipeline lifetime and may possibly lead to shut down of the pipeline.
Figure 10 illustrates the wall shear stress distribution, , around a slug unit for different values of   .It is obvious that the wall shear stress has the same distribution for all values of which starts with an increase in the wall shear stress near the bubble nose then it reaches a maximum positive value with the formation of a constant liquid film characteristics (thickness and velocity).The wall shear stress then starts to decrease once again till it reaches zero value in the bubble tail or wake zone.The results show that the wall shear stress is dependent on the distance between the Taylor bubble and the pipe wall, known as the liquid film thickness,   , which decreases with the increase in value of   .This is shown in Figure 10 where the maximum peak of the wall shear stress profile is decreased from   = 84 to higher values of   = 325.In addition, Figure 11 shows that there is a good agreement between the present results and the work of Taha [9].

Conclusions
In this work, a detailed numerical simulation of the hydrodynamic characteristics of gas-liquid slug flow in vertical pipe is developed using the volume of fluid (VOF) methodology implemented in the commercial software ANSYS Fluent.The main hydrodynamic characteristic studied are; the Taylor bubble shape and profile, the Taylor bubble rise velocity, the liquid film thickness and velocity, the wake flow structure, and the wall shear stress distribution.The simulation was done using air as the gas phase, and aqueous glycerol solution as the liquid phase.The following conclusions can be pointed out;    1.The Taylor bubble bottom depends on the liquid viscosity where the increase in the inverse viscosity number,   , increases the concave shape of the Taylor bubble bottom surface.
2. The calculated Taylor bubble rise velocity,   , is in an acceptable range when compared with experimental values and commonly used correlations in the literature.
3. The liquid film zone can be describes using the liquid film thickness,   , and the liquid film axial velocity, that are both directly affected by the inverse viscosity number   .
4. The wake flow structure has a closed axisymmetric nature for all the simulation cases with the development of circulatory vortex in the bubble wake with the increase in   .
5. The wall shear stress,   is mainly dependent on the liquid film thickness,   , and has a peak positive value in the stabilized liquid film zone.

Recommendations for Future Work
This study set out to develop a basic simulation model for gas-liquid slug flow in vertical pipe under laminar flow regime.It is recommended that further research be undertaken in the following areas; investigating the hydrodynamic characteristics of slug flow for different fluid system including the effect of viscosity and density ratios, investigating the hydrodynamic characteristics of slug flow under turbulent regime with > 500, exploring the hydrodynamic characteristics of slug flow including the flow of two consecutive Taylor bubbles in vertical pipe, and studying the wake flow pattern of single Taylor bubble or two consecutive Taylor bubbles under turbulent flow regime in terms of wake volume and length.
2. If  = 1; the cell is occupied mainly by the  h fluid.3.If 0 <  < 1; the cell contains the interface between the  h fluid and the other fluid.

Figure 1 :
Figure 1: Schematic of solution domain showing the boundary conditions.

Figure 2 :
Figure 2: Schematic of the main hydrodynamic features of two phase slug flow in vertical pipe.

Figure 3 :
Figure 3: Velocity fields (left) and streamlines (right) of gas-liquid slug flow in vertical pipe of 19 mm diameter and 209 mm length, N f = 84, Fr = 0.289, and Eo = 66.3 using moving reference frame.

Figure 5 :
Figure 5: Validation of simulation results for Taylor bubble shape profile for cases 1, and 2 with the work of (9) -is axial distance from bubble nose.

Figure 6 :
Figure 6: Effect of on Taylor bubble shape profile for cases 2, 3 and 4 -is axial distance from bubble nose.

Figure 7 :
Figure 7: Effect of   on liquid film thickens - is axial distance from bubble nose.

Figure 8 :
Figure 8: Effect of   on liquid film axial velocity - is axial distance from bubble nose.

Figure 10 :
Figure 10: Effect of N f on the wall shear stress distribution around a slug unit-x is axial distance from bubble nose.

Figure 11 :Volume 5 •
Figure 11: Validation of simulation results for the wall shear stress distribution around a slug unit with the work of (9) - is axial distance from bubble nose.

Table 1 :
Properties of the fluid systems used in the present study.

Table 2 :
Numerical and experimental values of Taylor bubble velocity, for all simulation cases given in Table1with corresponding deviations.