Optimization and Extended Applicability of Simplified Slug Flow Model for Liquid-Gas Flow in Horizontal and Near Horizontal Pipes

The accurate prediction of pressure loss for two-phase slug flow in pipes with a simple and powerful methodology has been desired. The calculation of pressure loss has generally been performed by complicated mechanistic models, most of which require the iteration of many variables. The objective of this study is to optimize the previously proposed simplified slug flow model for horizontal pipes, extending the applicability to turbulent flow conditions, i.e., high mixture Reynolds number and near horizontal pipes. The velocity field previously measured by particle image velocimetry further supports the suggested slug flow model which neglects the pressure loss in the liquid film region. A suitable prediction of slug characteristics such as slug liquid holdup and translational velocity (or flow coefficient) is required to advance the accuracy of calculated pressure loss. Therefore, the proper correlations of slug liquid holdup, flow coefficient, and friction factor are identified and utilized to calculate the pressure gradient for horizontal and near horizontal pipes. The optimized model presents a fair agreement with 2191 existing experimental data (0.001 ≤ μL ≤ 0.995 Pa∙s, 7 ≤ ReM ≤ 227,007 and −9 ≤ θ ≤ 9), showing −3% and 0.991 as values of the average relative error and the coefficient of determination, respectively.


Introduction
Gas-liquid, two-phase slug flow in pipes is a commonly observed flow pattern in many industries such as petroleum, chemical, nuclear, ocean engineering, power plant, etc. The slug flow pattern has a repeating cycle of liquid slug body and liquid film region, coming with the fluctuation of pressure loss (see Figure 1). Based on the visual observations by Dukler and Hubbard [1], the slug has higher kinetic energy than that of the liquid film [2]. This intermittency can cause mechanical vibrations in the pipe with high structural loads threatening the stability of the system. Therefore, the understanding and prediction of pressure loss have great importance, while it is usually complicated since the calculation process requires the iteration of many variables. In slug flow modeling for horizontal pipes, several pioneering models such as Wallis [3], Dukler and Hubbard [1], as well as Taitel and Barnea [4] were developed. These models commonly adopt the unit-cell concept. The unit-cell is an approximation based on the concept of an ideal slug unit, with a sharp change between slug and bubble regions. After Wallis [3], who introduced the equivalent unit-cell concept, Dukler and Hubbard [1] proposed the unit-cell model for horizontal flow. The balance equations are written in a frame of reference moving with the cell, so the flow appears steady [5]. Fully developed flow is assumed in both bubble and slug regions within the unit cell [6]. The unit-cell slug flow model constitutes nine parameters, which describe a slug unit moving with the translational velocity, vT. The sub-models in each region are often based on the two-fluid model for the bubble region, and the mixture model for the liquid slug section. In addition to the closure relations required for each sub-model, expressions for the bubble front velocity, the void fraction in the slug region, and a length scale (slug length or slug frequency) are necessary [7].
As Shoham [8] enumerated, Taitel and Barnea [4] suggested two methods for the pressure loss calculation. The first method consists of a global force balance on the entire slug unit, as formulated by Equation (1).
where HL Avg. is the average liquid holdup of the slug unit. When the uniform equilibrium thickness in the liquid film is assumed, Equation (1) can be approximated and reformulated to Equation (3).
where HLLS and HLTB is the slug liquid holdup and the liquid film holdup, respectively. The second method neglects the pressure loss in the liquid film (and gas pocket) region to calculate pressure loss in the slug unit, as written by Equation (6).
where ∆pMIX is the pressure loss in the mixing region at the front of the liquid slug body. To be more exact, the term ∆pMIX includes variations of both kinetic and potential energy as the upstream liquid film is picked up by the following slug body. However, this term generally represents the effect of acceleration (kinetic). More recently, Brito et al. [9] proposed a model to predict pressure gradient, simplifying Taitel & Barnea [4]. One of the main achievements of their model is the simple, fast, and convenient calculation of dp/dL with accuracy. In general, many mechanistic models developed for several decades require complex numerical calculations, including iteration. Furthermore, those calculations necessitate computational techniques to proceed with each step. The utilization of commercial software may provide the most reliable prediction of dp/dL, while it needs significant economic availability to use it. Meanwhile, Brito et al. [9] only targeted medium and high viscosity liquids (0.039 Pa•s ≤ μL ≤ 0.601 Pa•s), neglecting the pressure loss in the liquid film region (vLTB ≈ 0 and τF,G ≈ 0). From the visual observation with a high-speed video, Brito et al. [9] argued that most of the liquid is transported in the liquid slug. Their model is verified by experimental data obtained by Gokcal [10] and Brito [11] for the laminar flow condition, namely mixture Reynolds number up to 2000, while applying to the turbulent flow condition is not presented. As the authors limited the applicability of the model to relatively highly viscous liquids with laminar flow conditions, they simplified the friction factor as fS = 16/ReM, the flow coefficient as Co = 2.0, and utilized the correlation of Kora [12] for the slug liquid holdup, HLLS.
It should be noted that the suitable prediction of slug liquid holdup, HLLS, is essential for the proper calculation of dp/dL. The slug liquid holdup means the entrained gas void fraction in the liquid slug body. More entrained bubbles can reduce both mixture density and viscosity of liquid slug body, possibly decreasing the friction between the slug body and the inner pipe wall.
The understanding of HLLS should be started from the comprehension of entrainment at the slug front and deformation in the slug body. As reported by Kim [13], the main forces affecting the shape of entering bubble at the slug front are the shear of the surrounding liquid phase, buoyancy force, and the surface tension force. The entrainment mechanisms were observed to follow the plunging jet effect described by Kiger and Duncan [14]. The entering bubbles become sharpened as an increase in the shear, making an easy pass through the liquid slug front, simultaneously encountering higher resistance at higher surface tension. Unfortunately, a single unified model to predict HLLS has not been proposed appropriately, forcing to select the best correlations for different experimental and industrial conditions. A single bubble rising in a quiescent (or stagnant) liquid pool has been characterized by a number of previous studies. For instance, Luther et al. [15], Cieslinski and Mosdorf [16], Ohta et al. [17], and Wu et al. [18] investigated the behavior of single bubble in the airlift reactor. Bunner and Tryggvason [19], Lu and Tryggvason [20], and Ziegenhein and Lucas [21] observed the characters of bubbles in vertical bubbly flow. As Bunner and Tryggvason [19] commented, these studies are valuable for many industrial processes, such as boiling heat transfer, cloud cavitation in hydraulic systems, stirring of reactors, aeration in water purification, bubble columns and centrifuges in the petrochemical industry, cooling devices of nuclear reactors, and scavenging of dissolved gases in separation process. However, the bubble entrainment phenomena in liquid-gas two-phase slug flow have not been analyzed correctly as this pattern has a transient behavior of repeating liquid slug body and film region.
In this study, the simplified model proposed by Brito et al. [9] is modified and optimized to extend its applicability to turbulent flow conditions for horizontal and near horizontal pipes. The main postulation of neglecting the pressure gradient in the liquid film region is further supported by velocity fields measured by particle image velocimetry (PIV) in the previous study. The calculation of the wall shear stress of liquid slug body is advanced by adopting a suitable correlation of friction factor. The optimal correlations of slug liquid holdup, HLLS, and flow coefficient, Co, are identified and utilized to calculate the pressure gradient. The proposed pressure gradient model is verified by 2191 previously obtained experimental data of Gokcal [10], Brito [11], Kim [27] and [13], Ekinci [28], Mukherjee [29], and Kokal [30] with 0.001 Pa•s ≤ μL ≤ 0.995 Pa•s, 0.024 mm ≤ ID ≤ 0.076 mm, 7 ≤ ReM ≤ 227,007, and −9° ≤ θ ≤ 9°.

Methodology
Brito et al. [9] proposed a simplified pressure gradient model derived from Equation (7) (or Equation (3)) developed by Taitel and Barnea [4] as follows: where the ratio LS/LU is determined from the mass balance for the case of equilibrium liquid film as follows: where WL is the input liquid mass-flow rate, vLLS is the average liquid velocity in the slug body, and vLTB is the velocity in the liquid film. The authors argued low velocity of the liquid film (vLTB ≈ 0 and τF,G ≈ 0) with negligible accelerational pressure gradient (∆pMIX ≈ 0). The inner velocity field measured by Kim et al. [2] utilizing particle image velocimetry (PIV) also supports the postulation of low velocity of the liquid film, neglecting the pressure loss of liquid film. Kim et al. [2] performed the experimental study with 0.250 Pa•s ≤ μL ≤ 0.960 Pa•s for horizontal 50.8-mm (2-in.) ID pipe and laminar flow conditions. The in-situ velocity profiles at the center of pipe cross-section were measured and analyzed in liquid phase only. The measured velocity fields illustrate that the velocity of the liquid film is close to zero (or much slower than the liquid slug body), especially near the inner pipe wall (see Figure 2). It should be noted that the consecutive images of slug bodies in Figure 2 have some discontinuities caused by the limited laser frequency of the PIV system.
The illustrated results in Figure 2 were obtained with highly viscous flow conditions. This may be posed the question of "Is the low velocity of liquid film also valid for low liquid viscosity?" The recent experimental results by Gokcal [10,31] Kora [12], Al-safran et al. [32], Brito [11], and Kim [27] elucidated that slug length increases and slug frequency decreases as the liquid viscosity decreases. Gokcal [10] described that an increase in liquid viscosity caused a decrease in Reynolds number and the turbulence in the mixing zone at the front of the slug. Accordingly, the slug length decreases with the decrease of mixing length for increasing liquid viscosity. Experimental results of Brito [11] with relatively medium and high liquid viscosity (0.039-0.166 Pa•s) also indicated shorter slug lengths as liquid viscosity increases. These results may indicate that the effects of liquid film region on the total pressure loss becomes more negligible as liquid viscosity decreases. This conclusion might still be a rough postulation to support the neglect of the film region in a simplified slug flow model even with low viscosity conditions. A more rigorous experimental measurement may be required in the future study, although the application of the simplified model embraces all of the fluid conditions in this study.
Substituting Equation (10) into Equation (7), the largest pressure loss within the slug unit in horizontal and near horizontal flow can be approximated by Equation (11). Unlike the model suggested by Brito et al. [9], the first term remains for the near horizontal flow condition. In this study, the inclination range of near horizontal condition is chosen as -9 ≤ θ ≤ 9, arbitrarily.
where vLLS can be formulated using the mass balance equation of the liquid and gas phase.
The entrained gas-bubble velocity in liquid slug body, vGLS, has different values depending on its location in the slug. As visually observed by Kim [13], the vortices generated at the mixing zone in the front of the liquid slug body induce shear on the air phase as the air is further entrained into the slug body. The entrained bubbles are squeezed owing to the induced shear (see Figure 3). From the middle of the liquid slug body, the velocity of entrained bubbles is assimilated to the one of the liquid slug. Therefore, vGLS is postulated to be the contribution of the mixture velocity, vGLS ≈ CovM in this study. Based on the above assumption, Equation (12) can be reformulated to calculate the average liquid velocity in the slug body, vLLS, as follows: As a result, the final form of the simplified slug flow model for the horizontal and near horizontal condition can be expressed by Equation (14).
There is no term of slug length, LS, in the suggested model as the possible effects of liquid film region on dp/dL are neglected. This might limit the applicability of the currently proposed model to the case, such as boiling of water. For instance, the non-equilibrium slug flow model (NESM) for vertically upward flow described by Barbosa and Hewitt [33] emphasizes the generation of large Taylor bubbles. This phenomenon might be caused by the abrupt vapor growth in the sub-cooled near zero quality regions. As a result, the slug length could be a key factor when it comes to the heat and phase exchange. In this study, as Kim [13] experimentally observed, the amount of entrained bubbles in the liquid slug body is in a steady-state condition. On the other hand, this is only possible when the pipe is fully insulated which indicates negligible heat flux as the fluid (particularly water) temperature is equal to the ambient or pipe wall temperature.
Recently, Kim [13] reported that shorter slug length might cause the pseudo-slug flow at relatively lower mixture velocity conditions as liquid viscosity increases. The author presented the computational fluid dynamic (CFD) result of Pineda et al. [34], indicating 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, generating unstable slug body or so-called the pseudo-slug flow (see Figure 4). In the current study, the data sources of pressure gradient include the pseudo-slug flow (typically when vM ≥ 5 m/s), and the proposed model shows suitable agreements with them as will be presented in Section 3. However, a more rigorous analysis and modeling study will be required in future studies to understand the pseudo-slug flow mechanistically, considering its length scale. It should be noted that the proposed model postulates that the system temperature is constant following the axial direction of pipe, so the viscosities of each phase are constant. The originally simplified model only focuses on the laminar flow condition, i.e., ReM < 2000. In Brito et al. [9], the shear stress in the liquid slug was simplified using f = 16/ReM and τS = 8 μSvS/d. In addition, as the only highly viscous liquid was targeted, the slug liquid holdup correlation of Kora [12] was utilized. The application of the original model to the turbulent flow condition was inappropriate owing to the limited ReM range of HLLS correlation of Kora [12] and Fanning friction factor, and constant value of flow coefficient, Co = 2.0.
In this study, the improved composite correlation for the two-phase slug flow friction factor developed by Garcia et al. [35] is adopted to overcome the previous limitation. The proposed correlation includes a wide range of ReM that is applicable to both laminar and turbulent flow conditions (see Figure 5 and Equation (15)  In the case of the flow coefficient, Co, which needs to be optimized to calculate vLLS in Equation (14), several different correlations of translational velocity are compared to the previously obtained experimental data. Tables 1 and 2 summarize the experimental dataset and the compared correlations, respectively, to optimize the proper value of Co. The model of Choi et al. [36] requires the proper correlation of gas void fraction. Woldesemayat and Ghajar [37] improved the correlation of gas void fraction by Dix (Coddington & Macian [38]), validating it with a very extensive number of data, i.e., 2845 experimental points, as given in Equation (16). This model shows a good match, | | < 8%, with the measured average liquid holdup (HL Avg.) data obtained by Brito [11], Mukherjee [29], and Kokal [30] (see Figure 6). These data were not verified by Woldesemayat and Ghajar [37], previously. Therefore, their correlation is utilized to calculate the average gas void fraction in this study.  Likewise, the experimental data of slug liquid holdup, HLLS, and some correlations are listed in Tables 3 and 4, respectively. The best models of vT and HLLS for each experimental dataset, evaluated by six statistical parameters, relative performance factor, Frp, and the coefficient of determination, R 2 , are utilized to calculate the pressure gradient with Equation (14).
Six statistical parameters are used to evaluate the performance of the models, namely, the average relative error (ε1), absolute average relative error (ε2), standard deviation of relative error (ε3), average actual error (ε4), absolute average actual error (ε5), and standard deviation of actual error (ε6). Actual error (ei) and relative error (ej) expressed in Equations (17) and (18) The average relative error, ε1, indicates how large the relative errors are on the average.
The absolute average relative error, ε2, indicates how large the absolute relative errors are on the average.
The standard deviation of relative error, ε3, indicates the degree of scattering of the relative errors around their average value.
The average actual error, ε4, indicates the overall trend of the measured values.
The absolute average actual error, ε5, indicates the magnitude of the average absolute error.
The standard deviation of actual error, ε6, indicates the dispersion of the results around their average.
Tables 5 and 6 summarize the evaluation results of the best correlations in different experimental datasets. Figures 7 and 8 illustrate the comparisons between experimental data and predicted values calculated by the correlations.
In addition to the six statistical parameters, Al-safran [44] suggested to use the relative performance factor (Frp) given in Equation (25) and the coefficient of determination (R 2 ) given in Equation (26), as listed in Tables 5 and 6.
where the numerator of Equation (26) is the sum of the squares of the deviations between the calculated values and the average of all the experimental data. The denominator is the sum of the squares of the deviations between the measured values and the average of all the measured values.
The correlation that has the lowest value of performance factor compared to other correlations with the coefficient of determination close to 1 represents the best agreement with the experimental data.

Results
This section presents the performance of a modified and optimized slug flow model. The dp/dL experimental data measured by Gokcal [10], Brito [11], Kim [27] and [13], Ekinci [28], Mukherjee [29], and Kokal [30] are utilized to evaluate the performance of the model. The suitable correlations of Co (vT) and HLLS that are statistically selected in Section 2 are used to calculate dp/dL with Equation (14). The slug liquid holdup data were not measured by Gokcal [10] and Kim [27]. Therefore, the correlation of Gregory et al. [22] which shows the best agreement with the data of Brito [11] and Kim [13] is adopted for Gokcal [10] and Kim [27]. In the case of Mukherjee [29], both measurements of vT and HLLS were missed. Alternatively, the correlations of Fabre [6] and Gomez et al. [43] that present a good agreement with the data of Marcano [26] are used to calculate dp/dL. The fluid properties of Mukherjee [29] and Marcano [26] are similar as both of them used Kerosene for the liquid phase. Table 7 summarizes the pressure gradient dataset with chosen correlations for Co and HLLS.
As visualized in Figure 9 and evaluated in Table 8, the suggested model shows suitable agreement with experimental data from various conditions (0.001 ≤ μL ≤ 0.995 Pa•s, 7 ≤ ReM ≤ 227,007 and −9 ≤ θ ≤ 9). All of experimental data coincide with the calculated values lower than 10% of average relative error, ε1, except the data obtained by Brito [11] (ε1 Brito [11] = −16%). The coefficient of determination, R 2 , presents appropriate performance with values higher than 0.97, except with the data from Kim [27] (R 2 Kim [27] = 0.880). Figure 9. Comparison between the current dp/dL model and experimental data.

Conclusions and Discussion
In this study, the previously simplified model is modified and optimized to predict the pressure gradient of slug flow in horizontal and near horizontal pipe conditions. The proper utilization of friction factor, slug liquid holdup, and flow coefficient, modifying the balance equation, enables the proposed model to be applicable to turbulent flow, i.e., high mixture Reynolds number, and near horizontal conditions. The suggested model presents suitable agreements with 2191 previously obtained experimental data of Gokcal [10], Brito [11], Kim [27] and [13], Ekinci [28], Mukherjee [29], and Kokal [30] including 0.001 Pa•s ≤ μL ≤ 0.995 Pa•s, 0.024 mm ≤ ID ≤ 0.076 mm, 7 ≤ ReM ≤ 227,007, and −9° ≤ θ ≤ 9°. Whereas, the suggested methodology still requires the knowledge or the experimental data to select the proper correlations of slug characteristics such as the slug liquid holdup and the flow coefficient. A numerous number of data have been obtained during several decades from various experimental conditions. The intelligent optimization methodology employing the existing data sets, such as machine learning, might be a suitable route to exclude the subjective utilization of correlation unless the experimental data for necessary slug characteristics are provided.