In this study, detailed three-dimensional (3D) numerical simulations of intermittent multiphase flows were carried out to investigate the slug initiation process and various features of intermittent flows inside a horizontal pipe. Air and water are used as working fluids. The domain used for simulations is a 14.4 m long pipe with 54 mm inner diameter. The volume of fluid (VOF) model was used to capture the air/water interface and its temporal evolution. Using the developed computational fluid dynamics (CFD) model, the slug formation and propagation along horizontal circular pipe were successfully predicted and studied comprehensively. Slug length and the frequency of slug formation, as two main features of intermittent flow, were used to validate the model against experimental results and available correlations in the literature. Three-dimensional numerical simulation of intermittent flow proved to be a powerful tool in tackling limitations of experiments and providing detailed data about various features of the intermittent flow. The effect of gas and liquid superficial velocities on the liquid slug and elongated bubble length was explored. Moreover, the study revealed new findings related to the elongated bubble shape and velocity field in the slug unit.

## Introduction

Gas/liquid two-phase flows are encountered in many industrial applications such as condensation, evaporation, power generation, and crude transportation. One of the most important steps for understanding two-phase flows is the classification of flow distributions into various patterns which are known as flow patterns. Depending on the gas and liquid flow rates, pipe geometry, and pipe inclination angle, various flow patterns might appear in pipes which can be classified as stratified flow, wavy flow, intermittent flow, annular flow, and dispersed bubble flow.

The intermittent flow regime is one of the most common multiphase flow regimes that exist over a wide range of flow conditions in the horizontal two-phase flow [1]. An intermittent flow is generated when the liquid flow rate increases in a stratified flow regime. Therefore, waves grow until they occupy the entire pipe cross section and block the gas flow. Intermittent flows are characterized by the formation of large gas bubbles whose lengths are comparable to the pipe diameter, and they are known as “Taylor bubbles” or “elongated bubbles.”

These elongated bubbles tend to move along the top portion of the pipe due to buoyancy force. Two successive elongated bubbles are separated by a liquid body that may contain small gas bubbles and is called as “liquid slug.” Due to the transient nature of the intermittent two-phase flow, it is associated with large fluctuations in pressure and mass flow rate of liquid inside the pipes [2]. The accurate prediction of characteristics of intermittent flow is essential in the design of facilities where this type of flow occurs since these fluctuations may prevent sufficient delivery of liquid or damage the infrastructure. Therefore, it is important to identify and predict intermittent flow to prevent problematic issues. Although it has been studied experimentally and numerically for years, because of its complexity, a complete understanding of its nature has not been achieved yet.

where $A$ is the area of pipe cross section. $UL$ and $UG$ are the velocities of the liquid and gas phase, respectively. $\alpha $ indicates the volumetric fraction of each phase. Superficial velocity of a phase would be equal to the velocity of that phase if all pipe cross-sectional area was filled only by that phase.

Researchers are specifically interested in the prediction of the onset of the intermittent flow regime. Taitel and Dukler [3] proposed a mechanistic model to predict the horizontal flow patterns transitions in pipes. To anticipate initiation of an intermittent flow regime, they applied Kelvin–Helmholtz instability theorem on wave growth mechanism for inviscid fluids. Lin and Hanratty [4] improved Taitel and Dukler model by considering the effect of shear stress on the pipe wall using viscous long wavelength theorem. Currently, the slug stability theorem is used to describe the onset of intermittent flow. Based on this theorem, a volumetric balance between liquid attached rates at the slug front and liquid shedding rates at the slug tail imposes a criterion for the slug development [5]. Dinaryanto et al. [2] conducted an experimental study on slug initiation with air and water as working fluids and the inner pipe diameter of 26 mm. They concluded that slug initiation mechanism is different based on superficial velocities of liquid and gas. For high gas superficial velocity $(USG\u2009>\u20093\u2009m/s)\u2009$, wave coalescence was the main mechanism of slug initiation while for lower gas superficial velocity $(USG\u2009<\u20093\u2009m/s)\u2009$ wave growth was found to be the main mechanism. They also found that slug initiation position is closer to the pipe inlet at higher gas superficial velocity, but it is not very dependent on the liquid superficial velocity.

Considerable efforts have been devoted to predicting the slug hydrodynamic features such as slug length and frequency. Nydal et al. [6] experimentally evaluated liquid slug length with air and water as working fluids and pipe inner diameters of 53 and 90 mm. Mixture velocities were in a range of $7.5\u221223\u2009m/s$ through their experiments. They found out that slug length is fairly constant through all ranges of mixture velocities investigated. Rogero [7] also studied liquid slug length by conducting a series of experiments using $54\u2009mm$ inner diameter pipe. Based on Rogero's experimental results, slug length was about $15\u201320D$ for mixture velocities above 1.4 m/s. Generally, slug length is found to depend mostly on pipe diameter and working fluids rather than liquid and gas flow rates [8].

Slug frequency is defined as the number of slugs passing through a measuring device in 1 s. Gregory and Scott [9] proposed a correlation for slug frequency based on their experiments on carbon dioxide and water system. Their correlation was later developed for including the effects of inclination angle and the influence of mixture Froude number [10,11]. Heywood and Richardson [12] emphasized fluctuations in slug frequency and proposed a correlation to estimate average slug frequency. Dinaryanto et al. [2] found that slug frequency is not constant through the pipe and varies as intermittent flow develops. They also concluded that slug frequency will increase through all pipe length as liquid superficial velocity increases; however, it is not very dependent on the gas superficial velocity.

Although most studies on intermittent flow are experimental, numerical modeling of intermittent flow is becoming more popular in recent years. Despite the valuable information achieved through experimental investigations on intermittent flow, the complexity of the hydrodynamics requires methods to resolve some limitations of the experimental techniques. Computational tools have been proved to be very helpful to tackle limitations that are encountered in experimental investigations. The main obstacle for using numerical methods for intermittent flow modeling is its huge time expense [13]. Powerful computers are becoming more common today, which causes numerical simulations to be more practical in modeling large systems [14].

Different approaches can be used to model the transient nature of intermittent flows. Medina et al. [15] developed a slug tracking method to simulate the intermittent flow. The slug tracking method is based on the unit cell concept in which the numerical domain shrinks to an elongated bubble accompanied by a liquid slug. This method decreases the computational cost but requires information about the flow features (such as bubble and slug length, velocity, and volume fractions), usually obtained from empirical correlations or experimental tests, as an input to model or to validate it [16]. Issa and Kempf [17] studied slug initiation and development using a one-dimensional transient two-fluid model. The two-fluid model is based on solving continuity and momentum equations for each phase separately. They have developed a numerical code named TRIMOPH which was also used in the future studies of the group [18,19].

Due to the complexity of the intermittent flow, most recent publications use commercial CFD packages to simulate intermittent flows [13,20–23]. Intermittent flow simulations using fluent, cfx,transat and star-ccm+ can be found in the literature. The volume of fluid (VOF) and level-set are the mainly used methods for interface tracking of phases. Gupta et al. [24] compared the VOF method implemented in ansysfluent with the level-set method implemented in transat and concluded that both interface tracking techniques can be used to simulate intermittent flow within microchannels. Mehdizadeh et al. [23], Talimi et al. [13] and Khan et al. [20] employed VOF method in two-dimensional domains to simulate gas–liquid two-phase intermittent flow in microchannels.

Although several two-dimensional numerical simulations of intermittent flow were performed, few studies are found which performed simulations of intermittent flows in three-dimensional (3D) domains. Among the limited studies using 3D domains, Frank [25] employed ansyscfx 5.7 to simulate the slug formation and propagation along a circular pipe. He used a homogenous VOF method with shared velocity field assumption for the gas–liquid mixture. He used a transient liquid level at the inlet cross section to trigger intermittent flow in an 8 m pipe. A transient 3D air–oil intermittent flow simulation using VOF implemented in fluent was accomplished by Ban et al. [21]. Similar to Frank simulation, they also used transient liquid level at the inlet cross section to trigger intermittent flow formation. Although their simulation results are satisfactory, these simulations cannot be used to study the slug initiation process because gas and liquid flow rates supplied to the numerical domain are not steady.

Despite several experimental efforts were carried out to study the onset of intermittent flow, detailed numerical studies are rarely found on the subject. In the current work, 3D numerical simulations of intermittent flows were performed to study the slug initiation process and various features of intermittent flows inside a horizontal pipe. The slug length and frequency, as two main features of intermittent flow, are used to validate the model presented in this numerical study. This work is dedicated also to study intermittent flow in low superficial gas velocities $(USL<1\u2009m/s)$ which has not deeply studied. The computational tool ansysfluent is applied to predict less studied aspects of intermittent flow like elongated bubble shape and velocity field in and around the elongated bubble. During simulations, gas and liquid inlet conditions are fixed so slug formation and its frequency are not imposed as inputs to the numerical study.

## Governing Equations

$\alpha p$ is the volume fraction of phase $p$ and $|\u2207\alpha p|$ is the magnitude of the gradient of volume fraction of phase $p$. The detailed analysis of VOF model is well documented in Fluent theory guide [29].

## Computational Fluid Dynamics Model

The above-mentioned equations are solved using the finite volume method which is employed in the commercial fluid flow solver ansysfluent 16.1. The domain used for simulations is a 14.4 m long pipe with 54 mm inner diameter. At the inlet of the domain, the pipe cross section is divided into two semicircular sections. The upper semicircular section is air inlet and the lower semicircular section is water inlet. For water and air entering inlet boundaries, the velocity profile is assumed to be uniform. At the entrance of the pipe, there is a separator which extends by 20*D* along the pipe. The separator section is long enough to ensure a fully developed flow has been achieved. After the separator, the two phases are allowed to mix and form different two-phase flow patterns based on the values of superficial velocities. At the end of the computational domain, a pressure outlet boundary condition is applied. No slip velocity condition is imposed on pipe walls and the separator. The simulations are initialized as the upper part of separator is fully occupied with air, and the other parts of the pipe are occupied with water. In Fig. 1, the boundary and initial conditions are schematically displayed.

A 3D coordination system is used in order to yield the most realistic result of such an asymmetric flow pattern. In order to capture the flow development through the pipe, the length of the domain is 268 times larger than the pipe diameter. The grid used for this geometry consists of 621,000 hexahedral computational cells. Cross-sectional view of the grid used in simulations is illustrated in Fig. 2.

Finite volume method is used to discretize transport equations. The movement of the gas–liquid interface is tracked using VOF method. VOF method is the recommended multiphase modeling method when phases are immiscible, and the shape of the phases interface is of interest. Time discretization of phase continuity equations is through explicit scheme which yields more accurate results in time dependent simulations. Throughout simulations, geometric reconstruction and compressive schemes are applied as interface tracking schemes. Both of these methods are robust schemes to tackle numerical diffusions and capture sharp curvatures in phase interfaces. Geometric reconstruction scheme is more accurate compared to the compressive scheme in resolving interface curvatures but much more time-consuming. Additionally, implicit body force option is used, which improves the convergence for flows with large body forces like gravity.

Surface tension forces between phases are modeled through continuum surface force scheme which uses a nonconservative formulation. Turbulent flow modeled using $k\u2212\omega $ and turbulent damping option is enabled which adds source term in turbulence dissipation equation in interfacial cells. In intermittent flow modeling, turbulence damping helps resolving interfacial instability, which leads to the formation of slugs.

Within the CFD package, the pressure-based solver is applied. The pressure-implicit with splitting of operators is the chosen pressure–velocity coupling scheme. The PRESTO! (PREssure STaggering Option) interpolation scheme is used to compute the face pressure because of its higher accuracy especially near the interphases [13].

Courant number is defined as the ratio of particle moving distance during the assumed time-step to control volume dimension. High Courant number leads to unstable numerical approach and nonrealistic results. On the other hand, low Courant number leads to small time steps and consequently a large simulation time. In this study, the Courant number was fixed at 0.25 which provides a good balance between stable numerical approach and simulation time. The maximum number of iterations per each time-step was set to 200. Excluding the first time steps, each time-step was converged in less than 15 iterations. Scaled absolute values of residuals were set to $1\xd710\u22126\u2009$ as a convergence criterion for each time-step.

## Model Validation

In this study, the experimental measurements of Rogero [7] on intermittent two-phase flows were used to validate the numerical results. Scholars compare their models by available experimental results to validate them before presenting the final results of their model [14,33]. In this section, first, a grid size study is performed to show that our numerical results are mesh-size independent. Then, a validation test case is presented by comparing our numerical results with corresponding experimental measurements of Rogero [7] on intermittent two-phase flows. Moreover, some additional comparisons are made between our numerical results and available correlations in the open literature for slug flows.

### Grid Independency.

Grid dependency analysis is the first step for numerical model validation. In order to make sure that obtained results from our numerical flow solver are grid independent, the pipe geometry meshed with five different cell counts ranging from 120,000 to 2,160,000. A higher number of computational cells generally lead to a more accurate solution, but it is usually compensated with a higher computational cost [34]. In the cross section of the pipe, an O-grid mesh has been used as shown in Fig. 2. Elongated bubble velocity $(UEB)\u2009$ and minimum liquid film thickness below the elongated bubble $(\delta )\u2009$ are the criteria selected to check on the effect of numerical mesh size. Several simulations have been performed using different grid densities, and the results are summerized in Table 1. During the grid dependency analysis, air and water superficial velocities were $USG\u2009=0.8\u2009m/s$ and $USL\u2009=0.87\u2009m/s$, respectively. These are the highest velocities used for simulations which were selected because higher velocities are expected to be more sensitive to lower grid densities.

In Table 1, $UEB$ and $\delta $ are reported for different meshes, and corresponding deviations have been presented. Results are converged in the last three higher grid densities. The densest grid is used as a proper reference to calculate the deviations. It was concluded that 3D model with 621,000 hexahedral cells is fine enough to achieve reasonable mesh independent results. This mesh configuration is used for the simulations in this study.

### Validation of Computational Results.

A slug cell consists of an elongated gas bubble surrounded by a thin liquid film and accompanied by a liquid slug which separates two consecutive elongated gas bubbles. After the formation of a slug unit, it develops until the semi-equilibrium condition is reached in which its length becomes fairly constant. There is a good agreement between the length of liquid slug predicted by the numerical simulation, which is presented in Fig. 3, and previous experimental results presented by other authors, shown in Table 2. It should be noted that numerical simulations are identical to Rogero's [7] experimental investigations.

Based on Rogero's experiments, mean slug length is about 25–35*D* for mixture velocities less than $1\u2009m/s$. This mean length will decrease and become constant about 15–20*D* for mixture velocities above $1.4\u2009m/s$. Figure 3 corresponds to the simulation in which air and water superficial velocities were $USG\u2009=0.39\u2009m/s$ and $USL\u2009=0.87\u2009m/s$, respectively. Mixture velocity for this condition is 1.26 $m/s$. For the aforementioned simulation, as shown in Fig. 3, the size of the liquid slug is approximately 17*D*. This result is in good agreement with the experimental results presented in Table 2.

Recently, Hernandez-Perez et al. [39] proposed a correlation to predict the slug frequency in different pipe inclinations. Their formula is simplified to Gregory and Scott correlation for the horizontal pipes which confirm their correlation is still highly credited.

Two flow configurations were used to be compared with correlations. In case I, air and water superficial velocities were $USG\u2009=0.39\u2009m/s$ and $USL\u2009=0.87\u2009m/s$, respectively. In case II, air and water superficial velocities were $USG\u2009=0.6\u2009m/s$ and $USL\u2009=0.87\u2009m/s$. Table 3 presents the comparison between slug frequency predictions from most credited correlations and slug frequency results derived from the numerical simulations. Regarding that these correlations are not universal, there is a fairly good agreement between numerical results and correlations.

## Results and Discussion

A set of simulations were carried out in several gas and liquid superficial velocities. Figure 4 shows the flow pattern map for two-phase gas–liquid flows. As shown in Fig. 4, in a wide range of gas and liquid superficial velocities, the intermittent flow regime is the dominant two-phase flow pattern. Based on gas bubbles entrainment in liquid slug body, the intermittent flow regime can be divided into two sub regimes. In high gas superficial velocities, small gas bubbles are entrained in the liquid slug body due to the effect of inertial forces. While in low gas superficial velocities the liquid slug is free of small bubbles. Throughout this paper, the simulations were performed in low gas superficial velocity sub regime $(USG<1\u2009m/s)$. Gas and liquid velocities for each simulation case are selected based on the flow pattern map given in Fig. 4. Seven simulation cases were considered in which intermittent flow is expected to be the dominant two-phase flow regime and were named from A to G. Detailed data about these seven simulation cases are presented in Table 4.

### Slug Formation Process.

The slug formation process is well illustrated in Fig. 5. Figure 5(a) corresponds to $t=0(s)$ and shows that in the separator region of the pipe, air, and water are divided by a separator plate while the rest of the pipe is filled with water. In the first step, gas phase will advance in pipe along the top portion of the pipe as illustrated in Fig. 5(b). As a result of depletion and replenishment of liquid described by Taitel and Dukler [3], gas velocity will increase and form a bubble head at its front as shown in Fig. 5(c). When gas velocity increases, its pressure decreases, causing the liquid to come up more. Figure 5(d) illustrates the thinning process of the gas bridge which connects gas bubble to the main gas stream. This process continues until water completely fills the cross section of the pipe and forms an elongated bubble as depicted in Fig. 5(e). As slug unit develops along the pipe, the elongated gas bubble will become thicker and shorter. Meanwhile, liquid slug length increases until it reaches a semi-equilibrium condition. Figure 5(f) shows a slug unit in its semi-equilibrium condition. This process is common between all simulation cases performed in low gas superficial velocities ($USG<1\u2009m/s$).

Table 5 reports the equilibrium liquid slug length, the equilibrium elongated bubble length and the slug initiation position for each case of simulations. Throughout the simulations, the slug initiation position was in the range of 96*–*146*D* downstream of separator region. Taking into consideration the results reported in Table 5, one can conclude that as gas superficial velocity increases, the slug initiation position slightly moves toward the separator.

Formed slug units reached their semi-equilibrium condition at a distance between 30*D* and 50*D* downstream of the slug initiation position. Liquid slug length increased during slug unit development, but its equilibrium length was fairly independent of gas and liquid flow rates. This result is consistent with previous studies about liquid slug length. On the contrary, elongated bubble length increased as gas superficial velocity increased. Conte et al. [16] reached the same result about length of elongated bubbles while studying the features of intermittent flow experimentally in a 26 mm acrylic pipe.

### The Shape of Elongated Bubbles.

While several scholars have studied different aspects of liquid slugs, very little is known about elongated bubbles. For the simulations performed, elongated bubbles are 1.5–4.0 times larger than the liquid slugs (see Table 5). The contour plot of phase map at different cross sections of an elongated bubble is shown in Fig. 6. Due to the gravitational effect in horizontal pipes, elongated bubbles tend toward the upper portion of the pipe and form a liquid film underneath. The distance along pipe axis upstream of the elongated bubble nose is shown by *X*. For semi-equilibrium condition of case A, void fraction in pipe cross section sharply increases in the first 5*D* upstream of bubble nose. It reaches the highest quantity at a distance between 15*D* and 25*D* upstream of the bubble nose. Void fraction then decreases in the tail region of the elongated bubble.

Figure 7 depicts void fraction variation versus dimensionless distance upstream of the bubble nose for two simulation cases. It reveals more details about the elongated bubble shape. Figure 7 confirms that void fraction increases sharply in the first 5*D* upstream of bubble nose. Then void fraction grows slowly till it reaches its maximum quantity, which is about 70% for performed simulations. Finally, in the tail region, void fraction decreases monotonically till the end of the elongated bubble. By comparing the profiles of elongated bubbles as shown in Fig. 7, it can be concluded that in higher superficial gas velocities for which longer elongated bubbles are formed, the tail section is shorter and for shorter elongated bubbles a longer tail is expected. Conte et al.'s [16] experimental investigation reached the same conclusion while they did not provide a profile for the elongated bubble.

### Velocity Field in Slug Unit.

Figure 8 compares phase and velocity distribution at two cross sections of the elongated bubble. There is a meaningful relationship between phase and velocity distribution in each pipe cross section. The gas phase which forms elongated bubbles has a high velocity core. Along the walls, the velocity is reduced to satisfy no slip conditions on the walls. As peak velocity is located at the core of elongated bubble, higher velocity gradient is expected in the gas phase.

Axial velocity profiles at the bubble nose, bubble body, and bubble tail are depicted in Fig. 9. *X* and *r* represent the axial distance upstream of the bubble nose and the radial distance from the center of the pipe, respectively. Elongated bubble velocity has also been used as the reference velocity to achieve dimensionless axial velocity. As it appears from Fig. 9, velocity profile at the bubble nose (X = 0) is fairly similar to the single phase developed flow profile, except that the peak velocity is slightly shifted toward the upper portion of the pipe. It is reasonable as all previous research noticed the higher velocity of the elongated bubble compared to the liquid film. At 5*D* upstream of the bubble nose, this imbalance in velocity profile becomes more noticeable. Peak velocity occurs almost at elongated bubble core which is also well illustrated in Fig. 8. At the elongated bubble tail, the velocity profile obviously shows two peaks both in the upper half of the pipe cross section. While the first peak represents the high velocity gas flow, the second peak is the result of secondary circulating flow behind the elongated bubble where it becomes co-current with the main liquid stream.

## Conclusion

This study explores the ability of numerical modeling in the prediction of intermittent flow development and characteristics. The domain used for simulations is a 14.4 m long pipe with 54 mm inner diameter. As a result of grid independency analysis, it was concluded that 3D model with 621,000 hexahedral cells is acceptable to achieve reasonable converged results. Using VOF model, slug formation and propagation along horizontal circular pipe were successfully modeled and validated against experimental results. Liquid slug length and slug frequency were two intermittent flow features used to validate the numerical results.

Additionally, different features of intermittent flow were studied in detail such as the shape of elongated bubble and the velocity profile inside the slug unit. It is found that as gas superficial velocity increases, slugs tend to form closer to pipe inlet. But slug initiation position showed not to be dependent on the liquid superficial velocity. Formed slug units reached their semi-equilibrium condition at a distance between 30*D* and 50*D* downstream of the slug initiation position. During the development, the elongated gas bubble will become shorter and thicker. Meanwhile, liquid slug length increases to reach its final length. Elongated bubble length is a function of gas superficial velocity while liquid slug length is almost constant over a big range of gas and liquid superficial velocities. For larger gas superficial velocities that form longer elongated bubbles, the tail section of the bubble is shorter compared to smaller elongated bubbles. Velocity field in and around elongated bubble shows a high velocity core in gas flow. In tail cross section, circulation behind elongated bubble results in the creation of two velocity peaks.

## Nomenclature

- $A$ =
pipe cross section area $(m2)$

- $Ap$ =
interfacial area density of phase p $(m\u22121)$

- $B$ =
damping factor

- $Co$ =
Courant number

- $D$ =
pipe diameter $(m)$

- $f$ =
friction factor

- $fs$ =
slug frequency $(Hz)$

- $F$ =
surface tension force $(N)$

- $Fr$ =
Froude number

- $g$ =
acceleration of gravity $(m/s2)$

- $k$ =
turbulence kinetic energy $(J/kg)$

- $L$ =
length $(m)$

- $n\u0302$ =
surface normal vector

- $P$ =
pressure $(Pa)$

- $r$ =
distance in radial direction $(m)$

- $S$ =
source term in $\omega $ equation

- $t$ =
time $(s)$

- $U$ =
velocity $(m/s)$

- $V\u02d9$ =
volumetric flow rate $(m3/s)$

- $X$ =
distance from bubble nose $(m)$

### Greek Symbols

- $\alpha $ =
volume fraction

- $\beta $ =
closure coefficient for turbulence damping

- $\Delta n$ =
cell height normal to interface of phases $(m)$

- $\Delta t$ =
time step $(s)$

- $\Delta x$ =
time step $(s)$

- $\delta $ =
minimum liquid film thickness $(m)$

- $\kappa $ =
surface curvature $(m\u22122)$

- $\mu $ =
dynamic viscosity $Pa\xb7s$

- $\rho $ =
density $kg/m3$

- $\rho $ =
density $kg/m3$

- $\sigma $ =
surface tension $N/m$

- $\omega $ =
specific dissipation rate $(s\u22121)$