Simplified 1 D incompressible two-fluid model with artificial diffusion for slug flow capturing in horizontal and nearly horizontal pipes

The implemented numerical procedure has a little effect on the mass conservation since the threshold is very close to unity. In our model, the fluids are assumed to be incompressible, and, therefore, when a liquid slug exits the pipe a rapid change in the mass fluxes at the pipe outlet takes place. As a consequence, the integration error of the mass fluxes gives unphysical results. However, in every local section of the pipe the sum of superficial velocities (in the stratified region i.e. bubble) and the velocity of the liquid slug should be equal to the mixture velocity (to ensure the volume conservation). In Figure 1 (top) the ratio between the liquid velocity and the mixture velocity is plotted for a representative simulation showing that the ratio / is unitary in the liquid slug, i.e. = : thus, the threshold does not affect the imposed flow rates. Note that this condition is always satisfied in the stratified region.


Introduction
Slug flow is a highly intermittent flow regime (liquid slugs are followed by large gas bubbles at a randomly fluctuating frequency [1,2]) typical of many engineering applications, such as petroleum transport pipelines, chemical and nuclear industries and buoyancy-driven fermentation devices [3].This flow pattern may arise in horizontal and slightly inclined pipes from stratified flow because of the growth of instabilities: small perturbations appear on the liquid surface, and then, some may grow into larger waves to fill the pipe completely due to the well-known Kelvin-Helmholtz instability [4][5][6].Slug flow can also be established due to pipe slope changes, i.e., when the pipe inclination, initially horizontal or downward, turns upward: in this case, the liquid accumulates because of gravity, and the liquid volume fraction increases until the stratified-slug flow transition occurs [7].
In the past few decades, many numerical codes have been developed to numerically describe slug flow and to compute slugs' characteristics; one of the first methods was the 'unit-cell' approach, which enables a steady-state analysis of a control volume for gas bubbles and liquid slugs [8][9][10].Issa and Kempf [11] pointed out that "unit-cell" models cannot predict flow pattern transitions; moreover, steady-state models are not necessarily capable of predicting slug flow initiation in inclined pipes due to liquid accumulation.More recently, transient models were developed, and they can be classified mainly into three categories.The first category includes the empirical slug specification method, which provides slug characteristics (length, shape, speed, growth and decay) using some ad hoc correlations, by assuming a priori that slug flow is already developed [12][13][14].The second method is slug tracking: each single slug is tracked, i.e., the position of every slug tail and of every slug front is followed along the pipe using Lagrangian coordinates, and then, mass and momentum transport equations are solved at the slug front and tail.This approach was implemented in the commercial code OLGA, Oil and Gas simulator [15] and by Kjeldby et al. [16].These two methods are able to compute slug flow characteristics, but they cannot simulate the transition from stratified to slug flow, since the slug flow regime is a priori assumed.The third method is slug capturing, introduced in [11]: the transient mono-dimensional two-fluid model set of equations and closure laws, which describe the stratified flow regime, are solved numerically, and slug formation is an automatic outcome of the computed flow.In fact, slug formation, growth and decay arise naturally from the numerical solution of the two-fluid model for stratified flow, without introducing special correlations or other constraints.
Most of the exposed techniques and the commercial simulators are based on cross-sectionally-averaged one-dimensional (1D) models combined with coarse grids [15,17] to fulfil the strict requirements on computational speed.Unfortunately, as often pointed out [18][19][20], the 1D two-fluid model is ill-posed under certain conditions.To partially overcome this issue, Issa and Kempf [11] used grid diffusivity to dampen the short wavelengths' growth; this approach, however, does not ensure grid independence for every flow conditions since, as the grid is refined, the system is more sensitive to short wavelengths' instabilities.Fullmer et al. [21] propose to consider the viscous stress term (i.e., Reynolds stress) in the momentum equation and to account for the surface tension (see also [22]); the latter approach does not have a straightforward implementation because of the third-order derivative.Again, it would be possible to render the system hyperbolic with a physically-meaningful term such as the virtual mass, but this kind of description was developed mainly for dispersed flow; thus, in this work, this is not a viable strategy because we focus mainly on stratified and slug flow regimes.Another way to overcome the unbounded short waves' growth consists of adding some artificial diffusivity [23], introducing a second-order artificial diffusion term into the simplified two-fluid model, to obtain a formulation that provides convergent numerical solutions for flow conditions within the stratified and the stratified-wavy flow regime.Although the process of adding the artificial diffusivity in the numerical resolution of the simplified two-fluid model is not straightforward, we aim at exploring the possibility of this kind of regularization, which is based on sound mathematical principles: in fact, it has been demonstrated that if artificial axial diffusion is added to both mass and momentum equations, a cut-off wavelength is established below which all wavelengths are damped [23].
Afterwards, a transient roll-wave simulator based on a one-dimensional incompressible two-fluid model has been presented [24], together with a comparison of the numerical results with experimental data, showing good agreement; anyway, this approach does not take into account the slug flow regime because of the lack of a method to numerically describe the transition from two-phase to single-phase flow.In this work, instead, we developed a procedure to account for this phenomenon, which takes place when the liquid volume fraction tends to one as the slug body develops: this will extend the approach to slug flow.
The numerical description of the transition from two-phase to single-phase flow has been handled in very few papers, and at present, there is not fully consensual solution to cope with phase disappearance.Issa and Kempf [11] suggest to no longer solve the gas momentum equation and to set to zero the gas velocity when the gas volume fraction becomes negligible; however, they keep on solving the gas continuity equation, and when the gas volume fraction rises again above the prescribed threshold, the numerical solution of the momentum balance equation is re-introduced.Later, Nieckele et al. [25] followed this approach.Bonizzi et al. [26] introduced two thresholds to identify the transition from stratified to slug flow, one to observe when the gas volume fraction becomes small and one to observe the change of the volume fraction of the remaining gas, which must also be small.Zou et al. [27] give a special treatment of the momentum source terms if the volume fraction of the vanishing phase are lower than a certain value.Ferrari et al. [28] apply the approach by Issa and Kempf [11] to a five-equation two-fluid model, adopting a smoothing function between two thresholds and carefully handling the momentum source terms.Paillére et al. [29] describe an extension of the AUSM+scheme and propose to handle the velocity of the vanishing phase in a different way, if compared to [11]: they do not set the velocity arbitrarily to zero, but they make the vanishing phase velocity tend to the velocity of the remaining phase, thanks to a smoothing function.Cordier et al. [30] analyze the hyperbolicity of the two-fluid model and the loss of positivity of the numerical scheme in the presence of a vanishing phase.
The aforementioned transition criteria were developed for the four-equation two-fluid model [11,25,27,29], the five-equation two-fluid model [28], the six-equation two-fluid model [30] and the multi-field model [26]; as far as we know, no transition method has been developed for the case of a simplified two-fluid model, and thus, we develop here, for the first time, a completely new criterion for the two-equation model, on the basis of the previous literature.
Ansari and Shokri [31] developed a model to describe slug flow initiation, but its validity was restricted to the region of well-posedness of the two-fluid model: outside this region, the model is ill-posed, and the perturbation amplitude increases exponentially with decreasing mesh sizes.In the present paper, adopting a mathematical regularization, small wavelength instabilities are damped, and the model results in being well-posed.
In this paper, a simplified two-fluid model, regularized with artificial diffusion [23] is further developed introducing: • a criterion to properly describe the transition from two-phase to single-phase flow (which is required when a slug forms); • a method to choose the value of the artificial viscosity based on the fluid properties and flow rates, on the pipe geometry, on the space discretisation and on a desired cut-off wavelength; the latter is computed by solving the linear stability problem to dampen instabilities that possess a wavelength lower than a pipe diameter: this choice is coherent with the two-fluid model hypothesis, i.e., the resolution should not be lower than a representative length scale.
This paper is organized as follows: Section 2 presents the adopted simplified two-fluid model with artificial diffusion; Section 3 describes the numerical method, giving special care to the criterion developed to describe the transition from two-phase to single-phase flow.Section 4 shows the comparison of the numerical results with experimental data, when available, or with well-known empirical relations; finally, Section 5 reports some conclusions.

The Model
In this paper, we use a simplified version of the one-dimensional two-fluid model [23,32] for the derivation of the model.Similarly, in [33] and, more recently, in [34,35] and in [36], an analogous formulation of the two-equation two-fluid model is considered, which reduces the model from four to two equations.
The total mass conservation equation and the combined momentum equation for the stratified flow regime (see Figure 1) read: where the subscripts l, g and i stand for the liquid phase, the gas phase and the interface, respectively.
α represents the phase volume fraction, ρ the density, u the phase velocity, g the acceleration of gravity, β the pipe inclination, assumed to be constant along the entire pipe, and h the liquid height; A denotes the pipe cross-sectional area, τ the shear stresses and s the wetted perimeters.In this work, both fluids are considered as incompressible, and the effect of surface tension is neglected.The governing equations are completed by two algebraic relations: the volume fraction constraint, α l + α g = 1, and the relation between the superficial velocities and the mixture velocity, u ls + u gs = U m (t).The mixture velocity U m (t) is an input parameter, which can be expressed as a function of time, and it is constant in space due to the incompressibility hypothesis.We considered the value of U m as constant throughout the simulation for the sake of simplicity, but nothing prevents it from being a time-dependent function.
The two-fluid model is well-posed only in the region where the Inviscid Kelvin-Helmholtz (IKH) criterion [33] expressed by Equation (3) (see Appendix B) is satisfied: Since we are interested in the numerical description of slug flow, a regularization of the two-fluid model is required to prevent the unbounded growth of disturbances produced by small wavelengths' amplification beyond the IKH boundary.
A numerical regularization is embodied within the model introducing a second-order tensor diffusion (i.e., artificial viscosity) in the mass and in the combined momentum equations-Equations ( 1) and (2)-to neutralize the short wavelengths that are not correctly described by the two-fluid model [23].Even if the addition of artificial viscosity in both the mass and momentum equation may be physically questionable, it represents an acceptable compromise for engineering applications [37].
After including the artificial diffusion, the two-equation system becomes: ∂Ψ ∂t where F is the flux: and S is the source term: Closure relations for shear stresses are reported in Appendix A.
The vectors Q, Ψ and the matrix E contain the liquid flow variables, the independent variables and the constant artificial viscosity: This work focuses on the artificial viscosity to dampen instabilities at short wavelengths: thus, a method to set the element value of matrix E is described in Section 4.1.Holmås [23] designed the matrix E so that the initial value problem is well-posed for every flow condition (E must have real and positive eigenvalues), and all of the short wavelengths are stabilized.Coupling the artificial diffusion only to the liquid velocity in the momentum equation ( as the grid is refined, the short wavelengths grow bigger, and thus, numerical convergence is not achieved.The artificial diffusivity should be included also in the continuity equation and not only in the momentum equation [23], to obtain a cut-off wavelength under which all disturbances are damped.This approach aims at reproducing the regularization obtained by the surface tension description [18], and this can be achieved adopting a diagonal matrix, where both E 11 and E 22 have a positive value (namely, E 11 > 0 m 2 /s, E 22 > 0 m 2 /s and The method to choose the numerical viscosity values is discussed in Section 4.1: this reproducible method leads to a grid independent and regularized simplified two-fluid model for slug description; we aim at proposing a regularization that depends on flow and operative conditions and on physical considerations and that does not rely on an arbitrary choice of numerical viscosity values, as instead is done in [23].

Numerical Method and Two-Phase to Single-Phase Transition Modelization
In this section, we present the numerical method that has been embodied in an algorithm implemented in C language.In this work, the numerical method presented in [23] has been properly modified to account for the transition from two-phase to single-phase flow, which occurs during slug formation; we remark that, without this method, this approach would be applicable only in stratified and wavy flow and would be unable to describe slug flow; see Section 3.1.Moreover, boundary conditions have been consistently adapted to simulate a physical flow in pipelines.
First of all, the two-equation system, Equation (4), is split into an advection-source equation: and into a diffusion equation: to solve it with the second-order Strang operator splitting technique presented by [38].A finite volume method and an explicit time discretisation have been adopted to solve Equation ( 8): where the superscript n represents the time step, while the subscript j indicates the grid cell position.
The FORCE (First-Order Centered) scheme by [39] has been adopted to discretize the numerical fluxes F n j+1/2 and F n j−1/2 : Equation (11) states that each flux is computed as the mean of the fluxes obtained by the first-order Lax-Friedrichs method: and by the second-order two-step Lax-Wendroff method: where the predictor Ψ n+1/2 j+1/2 is: The implementation of the FLIC (Flux-Limited Centered) scheme [39], a high-order extension of the FORCE scheme, is discussed in the Supplementary Material.
The diffusion equation, Equation ( 9), is solved with the implicit Crank-Nicholson scheme discussed in [40]: this method consists of the numerical solution of a linear equations system to obtain Q n+1 j : This set of linear equations has been solved using the linear algebra package within the GNU Scientific Library.
The algorithm adopts an adaptive time step, respecting the Courant-Friedrichs-Levy condition: where λ n c,max is the maximum eigenvalue of the Jacobian of the advection equation, Equation (8), at time step n; the value of C has been set equal to 0.95 in all of the numerical results presented in Section 4, to obtain a large time step to optimize performances.When the liquid volume fraction is very close to one, the eigenvalues of the matrix A, which does not include the artificial viscosity, turn imaginary: here, the time step is computed on the basis of the absolute value of the complex number.This choice ensures that the time step is smaller (and therefore safer) than the one computed on the basis of the real or of the imaginary part.
Boundary conditions are imposed in the following manner: velocity and volume fraction are extrapolated from the last cell at the outlet, while constant liquid superficial velocity and volume fraction are imposed at the inlet; at the beginning of the simulation, the liquid volume fraction and the velocity are assumed to be uniform along the pipe, representing a stratified smooth flow regime.

Two-Phase to Single-Phase Transition Criterion
During the slug onset process, the transition from two-phase flow to single-phase flow occurs: when this happens, the liquid volume fraction grows and tends to unity.Concerning the numerical simulations, if the liquid volume fraction value is not bounded or its growth not prevented, it becomes bigger than one, which is not physically acceptable.Moreover, if a proper numerical treatment of this phenomenon is not introduced, the simulation stops because the numerical fluxes cannot be computed.We need, therefore, to introduce a criterion to handle these numerical issues; as detailed in Section 1, only a few works have developed a similar treatment, and an ad hoc criterion for the simplified two-fluid model has not been proposed yet.
During the model development, we observed that, when the liquid volume fraction tends to one, the gas velocity required in the numerical fluxes of Equation (11) increases: this happens because of the small value of the denominator in the equation employed to compute u g from the mixture velocity, the liquid velocity and the liquid volume fraction: Using this relation when α l → 1 leads to diverging values of u g , which are inserted in the numerical fluxes and, consequently, in the liquid velocity.This phenomenon has been observed by Zou et al. [27]: the velocity of the vanishing phase diverges and causes non-physical behaviors in the other variables.Our solution to prevent this phenomenon, due only to the numerics, consists of no longer computing gas velocity employing Equation ( 17) above the threshold α l > 0.999, but forcing it to zero.Together with this, when the liquid volume fraction rises above the prescribed threshold (which was set equal to 0.999), also the hydrostatic term ∂h/∂x is removed from the combined momentum equation, since it is suitable just for two-phase flow.In addition, we set the interfacial friction and the gas-wall friction to zero, since the slug region contains only liquid.These adjustments prevent the uncontrolled growth of the gas velocity and help the liquid volume fraction stay bounded.Anyway, after the third step of the Strang operator splitting, a check on the value of the liquid volume fraction has to be performed, and if it exceeds the threshold value of 0.999, it is set to this value (a similar limitation on the volume fraction value is exposed in [27]).Thus, the idea is that when the liquid volume fraction exceeds a prescribed threshold, the terms concerning the gas phase are properly switched off to let the code manage slug formation or breakage.Since the numerical fluxes regarding the mass equation have not been modified by our criterion, the value of α l is free to evolve, and at the slug tail, it drops under the threshold: when this happens, gas velocity is computed using Equation (17), and the terms regarding the gas phase are reintroduced into the calculations.
Further information (the effect of the transition criterion on mass conservation, its validation by the water-air separation test [41]) is reported in the Supplementary Material.

Numerical Results
This section presents the results of the previously-described code, aiming at simulating slug flow in practical application.As pointed out earlier, the two-fluid model is able to predict slug growth and development in an automatic way, starting from uniform stratified flow as the initial condition.Thus, in this work, no perturbations are needed to initiate slug flow, since the adopted two-equation model descends from the two-fluid model itself.
We chose to present the results referring to two different pipe geometries to show that the code can be used to obtain predictions on real systems.The first configuration, which will be referred to as Pipe1, represents an experimental setup available in the Fisica Tecnica Laboratory at the University of Brescia, where many experiments have been performed in order to measure slug characteristics [42,43].The second pipe configuration will be referred to as Pipe2 [11], whose geometrical properties and superficial velocity ranges are shown in Table 1.Both pipes are horizontal (i.e., β = 0 • ).The adopted fluids are water and air, and their physical properties (density and viscosity) are reported in Table 2.  Slug initiation and slug characteristics will be discussed in both cases, and numerical results will be compared either with empirical models or with experimental results.Finally, the algorithm's computational performance will be discussed.

On the Choice of the Artificial Viscosity Values
In this section, we will discuss the choice of the value of the elements E 11 and E 22 in the artificial viscosity matrix E, Equation ( 7), adopted to regularize the simplified two-fluid model.In previous work [23], the values where arbitrarily set, i.e., E 11 = 0.001 and E 22 = 0.01.Instead, we aim at providing a consistent method that can be applied for every flow condition; artificial viscosities are set case by case to obtain the desired cut-off at short wavelengths, and the choice is based on the analysis of the linear stability results.The two-fluid model is a one-dimensional averaged model: therefore, the assumption that the model resolution should not be lower than a representative length scale (i.e., the pipe diameter [37]) is consistent with the model hypothesis.The unbounded growth at short wavelengths is a consequence of the averaging process and of the absence of the surface tension description in the model; thus, numerical regularization aims at reproducing the small wavelengths' cut-off.
We used the results of the linear stability analysis to set the artificial viscosity values to obtain a cut-off wavelength of a representative pipe diameter, i.e., 0.022 mm in the case of Pipe1 and 0.078 mm in the case of Pipe2: the values of E 11 and E 22 are set so that the imaginary part of the amplification factor is lower than zero for wavelengths lower than a pipe diameter, i.e., Img(ω) < 0 for λ/D < 1.In Figure 2, the main steps of the algorithm are reported.In this way, we can obtain a well-posed model, and we set a lower cut-off to the size of perturbations allowed to grow in the case of unstable flow.Following the same approach proposed in [23], E 11 is set an order of magnitude lower than E 22 .The linear stability analysis is performed on the basis of the inlet superficial velocities, and the methodology is conceived of to provide an order of magnitude to the esteem of E 11 and E 22 , not a punctual value; thus, even if flow conditions along the pipe can differ from the inlet ones, the esteem provided by the proposed method is still acceptable.Figure 3 shows, for the investigated geometries and for different couples of superficial velocities, the amplification factor Img(ω) as a function of the wavelength, which is defined as the imaginary part of the propagation velocity [23,33].If Img(ω) > 0, the corresponding wavelengths are amplified; if Img(ω) < 0, they are damped.The values of E 11 and E 22 are the minima necessary to obtain a cut-off wavelength no higher than one pipe diameter, and they are specific for each couple of superficial velocities.The presence of the spatial discretisation introduces additional numerical diffusivity: in Figure 3, at the left-hand side of the black vertical line, all of the waves with wavelengths smaller than 2∆x, where ∆x = D/2, are damped.Anyway, as previously discussed, artificial diffusivity must be included.Once the cut-off wavelength has been set at one pipe diameter, the spatial discretisation must be at least half the diameter of the pipe or smaller, since a wavelength of 2∆x is the smallest representable by a ∆x grid.The independence of slug statistical characteristics from the spatial discretisation is shown in the Supplementary Material.

Flow Pattern Map and Code Validation
This section shows that the code is able to predict the transition between flow regimes (stratified and slug flow).Adopting the geometries described previously, several simulations were performed at different gas and liquid superficial velocities; for this purpose, the results of these simulations are plotted in Figure 4 and compared to flow transition boundary criteria.The Inviscid Kelvin-Helmholtz (IKH) boundary represents the transition between the well-posed and the ill-posed model (without including the artificial diffusivity) [33]: its points are the couples of superficial velocities that satisfy the equality in Equation (3).The Viscous Kelvin-Helmholtz (VKH) line represents the linear stability boundary under the long wave approximation; see the Supplementary Material for details.The artificial viscosity becomes negligible under the long wave approximation, but as shown in [23], the diffusion does not affect in a sensible way the VKH line.Thus, the stability diagram is divided into three main regions, and the flow regimes can be placed on the flow pattern map in the following manner: under the VKH line, the flow is stratified smooth, since the VKH curve represents the border under which the smooth stratified flow is stable; above the IKH line, the flow pattern is slug flow, since here, we observe the unbounded wave growth; between the VKH and the IKH line, the flow can be stratified wavy or slug, since this is a transition region between the other two flow pattern regimes.Moreover, the dashed slanting line represents the couple of superficial velocity that has correspondent equilibrium hold up of 0.5: if the hold up is higher than 0.5, in the wavy region, there exists a higher probability that some waves may grow, completely fill the pipe section and generate slugs [44].
With the only purpose of better discerning the emerging flow regime, between stratified smooth and stratified wavy flow, numerical simulations were performed imposing a sinusoidal perturbation on the inlet volume fraction, respecting the following criteria: • if the imposed perturbation was absorbed, the flow regime was classified as stratified smooth; • if the perturbation was neither adsorbed, nor amplified until slug development, the emerging flow regime was wavy stratified; • if the perturbation was amplified until the liquid completely filled the pipe, the slug flow regime was recognized.
Figure 4 shows that the code predicts the transition from stratified to wavy flow and from wavy to slug flow.It is possible to observe that, except for two points in Figure 4b, all of the smooth stratified results are under the VKH line, and then, before obtaining the slug flow regime, the transition to the wavy regime takes place; above the IKH line, the slug regime is obtained, as expected from theory.
In the Supplementary Material, a validation against the experimental flow pattern map reported in [45] is exposed.

Slug Characteristics
Several computations were carried out with different couples of superficial velocities for both geometries, with an end time of 300 s and a grid discretisation ∆x = 0.256D to estimate slug characteristics, such as slug mean velocity, slug mean frequency and slug length.Results have been compared to experimental data, in the case of Pipe1, or with well-known empirical models, in the case of Pipe2, in order to validate numerical predictions.An analysis of the effect of the grid discretisation and of the threshold adopted in the transition from two-phase to single-phase flow (see Section 3) and a comparison against the slug initiation phenomenology proposed in Ansari [46] and in Ansari and Shokri [31] are reported in the Supplementary Material.

Slug Mean Velocity
Figure 5a shows, on the left, the comparison of the slug mean velocities obtained numerically by the two-equation model and the measured ones, for the geometry Pipe1, while Figure 5b compares slug computed translational velocities, in the case of Pipe2, and the predictions given by the relation reported in [47]; their correlation (u slug = C 0 • U m + u d , where C 0 is the distribution parameter and u d is the drift velocity) is for bubble velocity, but we adopted it under the hypothesis that the slug front moves at the same speed of the bubble before it; according to [47], we consider C 0 = 1.2 for turbulent flow (the case of this study) and u d = 0 for horizontal flow.In both the cases, the agreement between the numerical results and the experimental data (or empirical prediction) is satisfactory, since they are within the ±20% lines.Moreover, at the right-hand side of Figure 5a, we show that the computed slug velocities (black rings) follow accurately the trend predicted in [47] (blue line).

Slug Mean Frequency
Slug mean frequencies computed on the Pipe1 geometry are compared in Figure 6a with the experimental results: the majority of the computed points are within or very close to the ±30% bound.For Pipe2, mean slug frequencies are reported in Figure 6b, where computed mean slug frequencies have been compared with the empirical correlation [48]: the slug frequencies computed with the two-equation model are in good agreement with this empirical correlation, and only a few points are out of the ±30% bounds.We can state that the numerical results compare quite well both with experimental data and with the proposed empirical correlation trend, since 87% of the analyzed numerical results are in the ±30% bounds, in the case of Pipe1, and 93%, in the case of Pipe2; however, the predictions are more accurate in the case of Pipe2, since in this case, 67% of the analyzed results are within the ±20% bounds, while only 27% in the case of Pipe1.

Slug Body Length Distribution
The computed slug body length distributions are reported in Figure 7: the typical log-normal distribution identified by experimental observation [49] is reproduced (the log-normal curves displayed in Figure 7 were fitted to the numerical data: the obtained mean and standard deviation are reported for each of the four cases).As concerns Pipe1, the experimental evidence indicates that the slug body lengths are typically in the range of 10-20-times the pipe diameter, while for Pipe2, the typical slug lengths are around 12-30 pipe diameters.In Figure 7a,b, we can see that the computed slugs are shorter than the ones observed experimentally, but that the statistical behavior of the distribution is preserved.The discrepancy between the computed slug lengths and the experimental ones could be due to the simple closure laws adopted in this work, because they do not describe faithfully the physical process that makes the slugs grow, that is the liquid slug deceleration given by the liquid-wall shear stresses.Moreover, in this model, the gas entrainment is not taken into account, and this also plays a crucial role in the slug development.

Computational Performances
In this paragraph, the computational performances of the presented algorithm will be briefly discussed, with the purpose of giving to the reader the possibility to compare the performance of the presented code with the ones of other numerical simulators.Several numerical tests have been performed with many couples of superficial velocities, all of them with and en time of 300 s and a grid discretisation ∆x = D/2.All of the simulations were performed serially on an Intel Xeon CPU, 2.30-GHz processor.To show the computational performances of the code on this machine, we define the ratio between the computational time and the simulated time: If the parameter Θ is lower than one, the simulation runs faster than real time.Figure 8a,b shows the dimensionless computational time Θ resulting from simulation respectively for Pipe1 and Pipe2; as we can see, keeping the liquid superficial velocity constant, the simulation time grows almost linearly, except for some points, as the superficial gas velocity increases: this velocity-dependent simulation time is due to the adaptability of the time step on the basis of the eigenvalues, which are strongly related to the phases velocities; as the phase velocities increases, also the eigenvalues grow, and this leads to smaller time steps; smaller time steps mean that more computational iterations are needed to reach the desired end time, and thus, a higher computational time is required.Although the computational time is very case specific, depending on the pipe length, on the grid discretisation and on the flow rates, the code shows promising computational performances if we evaluate the computational time of the presented example cases.Moreover, to further improve the computational efficiency, it would be possible to run the code in parallel: the conversion from a serial to a parallel architecture can be done, for example, using the domain decomposition method and the OpenMP API.The algorithm consists of dividing the pipe domain into N chunks, and the numerical resolution of each chunk is assigned to a thread.The chunk-thread assignment and the exchange of information among the threads (fundamental to send and receive boundary conditions) can be done using the features available in the OpenMP API.Once each thread has finalized a time iteration on the assigned chunk, boundary conditions are exchanged with the neighboring threads, and the subsequent iteration can take place.When the desired end time has been reached, the solution is reconstructed from the N solutions.As concerns practical applications, such as the preliminary design of the operative setup, where fast predictions are required, this code can be a a reliable tool for this kind of purpose.

Conclusions
In this work, a numerical code for slug capturing in a horizontal pipe was developed and implemented.We adopted a simplified two-fluid model regularized with artificial diffusion, to deal with the issue of ill-posedness.An original numerical procedure to describe the two-phase to single-phase flow transition was designed and implemented, thanks to which the numerical code is able to simulate the slug flow regime, as well as smooth and wavy stratified flow.A consistent choice of the coefficients of the artificial diffusion matrix was proposed, which depends on the pipe configuration, as well as inlet flow rates.
The numerical algorithm was tested on two different pipe geometries and validated against experimental results; it was shown that the obtained numerical results agree well with experimental data and with well-known empirical correlations.The developed simulator was proven able to describe in a satisfactory way the transitions between flow regimes (i.e., stratified smooth, stratified wavy and slug) and to compute well slug characteristics, especially mean frequency and mean velocity.As concerns slug lengths, the experimentally-observed log-normal distribution clearly emerges from the numerical results, but the computed mean length is underestimated: this could be because of the simple form of the adopted closure laws and the omission of gas entrainment; these issues should be, in our opinion, addressed in further works.
Finally, the computational performances of the code were discussed, showing that the simplicity of the model leads to a computationally-inexpensive resolution: this feature and the possibility of a straightforward conversion of the algorithm on a parallel architecture render the proposed simulator as a good candidate for applications where fast predictions are required.

Figure 2 .
Figure 2. The algorithm implemented to obtain the values of E 11 and E 22 .

Figure 3 .
Figure 3. Plot of the amplification factor for different couples of inlet superficial velocities.(a) Amplification factor for geometry Pipe1; and (b) amplification factor for geometry Pipe2.

Figure 6 .
Figure 6.Slug mean frequency: several numerical results are bounded by the ±30% confidence lines.(a) Pipe1: numerical results against experimental data; and (b) Pipe2: numerical results against empirical correlation.

Figure 7 .
Figure 7.Typical slug length distribution for the presented geometries.(a) Pipe1.Left: u sg = 1.1 m/s and u sl = 1.3 m/s.Right: u sg = 1.5 m/s and u sl = 1.6 m/s; and (b) Pipe2.Left: u sg = 2.0 m/s and u sl = 1.0 m/s.Right: u = 2.5 m/s and u sl = 1.0 m/s.

Figure 8 .
Figure 8. Dimensionless computational time Θ as a function of the superficial gas velocity: simulations with the lower superficial velocities are performed faster than real time.End time t = 300 s, ∆x = D/2, CFL = 0.95.(a) Dimensionless computational time Θ for the simulations performed on Pipe1; and (b) Dimensionless computational time Θ for the simulations performed on Pipe2.

Table 1 .
Geometrical characteristics, mesh sizes and superficial velocity ranges of the pipes analyzed in the numerical results.

Table 2 .
Physical properties of the adopted fluids.