Prediction of Horizontal Gas–Liquid Segregated Flow Regimes with an All Flow Regime Multifluid Model

Abstract

The generalized multifluid modelling approach (GEMMA) has been developed to handle the multiplicity of flow regimes and the coexistence of interfaces of largely different scales in multiphase flows. The solver, based on the OpenFOAM reactingEulerFoam family of solvers, adds interface resolving-like capabilities to the multifluid solver in the cells occupied by large interfaces. In this paper, GEMMA is further developed to predict stratified and slug flow regimes in horizontal ducts. The suppression of the turbulence and the wall-like behaviour of large interfaces is modelled with an additional dissipation source. This enables an accurate prediction of the velocity and of the turbulence kinetic energy in a stratified channel flow and the capturing of the formation and the travel of liquid slugs in an annulus. Large interfaces are identified and tracked, not only in the smooth and wavy stratified regimes but also in the much more perturbed interfaces of liquid slugs. The present work confirms GEMMA to be a reliable approach to provide all flow regime modelling capabilities. Further development will be focused on large interface momentum-transfer modelling, responsible for the overestimation of the interfacial shear and the limited liquid excursion during slugs, and the extension to interface break-up and the entrainment of bubbles and droplets, to handle the entire range of regimes encountered in horizontal flows.

1. Introduction

The thermo-fluid dynamics of multiphase gas–liquid flows is often driven by the transport of mass, momentum and energy across the interface. The rate and physics of these transport mechanisms is proportional to the interfacial area, the density of which strongly depends on the interface morphology and the spatial arrangement between the phases involved. From the phase arrangement, multiphase gas–liquid flows have been historically classified into different flow regimes [1,2]. In dispersed regimes, a large number of relatively small gas bubbles or liquid droplets are dispersed inside a continuous flow of liquid/gas. On the opposite side of the spectrum, in segregated regimes, large continuous structures of both phases coexist, separated by a smaller number of large-scale interfaces.

A typical example of the latter are stratified flows where, driven by gravity, gas and liquid flows are separated by a large horizontal free surface. The stratified flow regime is relevant in the thermal hydraulics of nuclear reactors and multiphase processes in multiple engineering sectors, such as oil and gas. In water-cooled nuclear reactors, the development of co-current and counter-current stratified flow conditions during loss of coolant accidents can threaten the cooling capability of the system and the structural integrity of the reactor vessel and of the cooling system [3,4]. Stratified flows are a key feature in reactor safety problems such as pressurized thermal shock, condensation-induced water hammer and counter-current flow limitation. Therefore, the accurate modelling of their thermo-fluid dynamics will positively impact the accuracy and reliability of safety analyses on the loss of coolant accidents progression. In nuclear thermal hydraulics, as in the majority of other applications, smooth stratified flows rarely exist, but the interface is often perturbed by waves of different length scales, the entrainment of bubbles and the displacement and deposition of liquid droplets [5,6]. Larger perturbations can lead to flow regime transition with the formation of intermittent liquid slugs or a thin continuous liquid film around the periphery of any conduit. The formation of slugs is a well-known issue for oil and gas production and transport, causing loss of productivity and efficiency and extra fatigue on equipment, which have led to significant investments and efforts into reducing or eliminating the occurrence of slug flow [7,8]. Clearly, there are many advantages to be gained from the accurate numerical prediction of stratified flows, specifically if a model can be extended to predict the perturbation of the interface and the transition from stratified to other flow regimes, including detecting and potentially preventing the occurrence of slug flow conditions.

However, the modelling of flow regime transition and the broad range of interface scales that are found, especially in transition and mixed conditions, are still a challenge even for computational fluid dynamics (CFD) [9,10]. In stratified flows, the large, free-surface interface can be solved with relative ease by interface-resolving methods, such as the volume of fluid (VOF) or level set approaches, where all the interfacial scales and transport processes are resolved on the computational grid [7]. However, as soon as the grid needs to be refined enough to capture the perturbations of the interface, or small-scale bubbles that are entrained in the liquid phase or droplets that are displaced by the gas flow, the computational cost becomes prohibitive [11]. Averaged Eulerian–Eulerian multifluid models, where interfacial transport is entirely modelled, are better equipped to model these processes, and in more general dispersed flow regimes [12,13], but they are instead limited in the modelling of large, intermittent interfaces.

To resolve this scale-induced problem, different authors have recently coupled interface resolving-like capabilities to averaged multifluid models [4,14,15,16,17,18,19]. These interface resolving capabilities are only used with large interfaces, while small details and dispersed regimes are modelled with the multifluid approach, reducing significantly the overall computational load. Most models include compression forces that counteract the diffusion of sharp interfaces [14,16,19,20], which is a major drawback in multifluid models, and the capability to select an appropriate closure framework, usually including surface tension and a large-scale drag model when large interfaces are present. Most models identify the latter from an a priori function of the local volume fraction, used to distinguish between dispersed and continuous regions, but without accounting for the local morphology of the interface [15,16,20]. In contrast, the recently published generalized multifluid modelling approach (GEMMA) [21] avoids using a priori functions of the volume fraction and directly detects and tracks large interfaces from the local mesh resolution and the morphology of the interface.

In this work, the GEMMA formulation is further developed and applied to horizontal and annular smooth and perturbed stratified flows. The accuracy of the solver is evaluated, with a focus on two specific aspects. One is the drag modelling framework employed and the model used with large-scale interfaces [22], which can be expected to provide the vast majority of the interfacial drag in stratified flows. The second aspect is a well-known issue found in stratified flow, which is the modelling of turbulence near the large interface. This interface has an influence on the flow similar to that of a solid wall, with the additional complication that the interface is a moving surface subject to deformation, instabilities and changes of shape. Turbulent fluctuations decrease near the interface, especially on the gas side, and if the density difference between the phases is large [23,24,25].

When Reynolds-averaged Navier–Stokes (RANS) models are employed to model turbulence, the large velocity gradients near the interface produce unphysically high levels of turbulence that can only be avoided with dedicated damping mechanisms [26]. Among the methods employed, a few researchers have relied on wall-function type treatments in the near-interface computational cells, which, however, need to be spatially located first with targeted algorithms [4,27]. A more popular option has been to rely on source terms that force the turbulence dissipation to obey to a wall-like behaviour. This kind of damping was introduced first in the ω equation of the k-ω model by Egorov et al. [26] and multiple authors have followed this approach [15,28,29]. Instead, Frederix et al. [23] adapted the ω source to the k-ε model and introduced a new interface related length scale that is independent of the computational grid.

The GEMMA model is first used to simulate the stratified air–water flows in a horizontal channel studied by Fabre et al. [30]. The turbulence damping model from Frederix et al. [23] is adapted to the GEMMA formulation and tested. The model is included in a generalized interfacial turbulence formulation that is already available for additional physics to be included, such as bubble-induced turbulence in bubbly flows or turbulence generation from waves and perturbations at large interfaces. Two experiments at different Reynolds numbers from Fabre et al. [30] have been simulated. Successively, GEMMA is used to predict the occurrence of slug flow in an annular channel configuration that has been investigated numerically by Friedemann et al. [7] with an interface resolving method and experimentally at the Institute for Energy Technology (IET) in Norway. The annular flow configuration is relevant to the oil and gas sector, where it can be found in multiple circumstances such as between the drill pipe and casing during drilling or in multistage horizontal fracking systems [7]. Results in stratified flows aim at cementing the GEMMA solver as a promising all flow regime multiphase computational tool, capable in this context of handling the transition to wavy and slug flow regimes for which different future developments are identified.

2. GEMMA Model

The GEMMA solver is built on top of the native reactingMultiphaseEulerFoam solver of the OpenFOAM 7.0 CFD package [31]. reactingMultiphaseEulerFoam is a multifluid solver capable of handling a system of n compressible phases and equipped to solve mass, momentum and energy conservation equations for each phase. GEMMA adds to the multifluid solver a large-scale formulation to handle large interfaces (in this context, large is a relative term and needs to be related to the local resolution of the computational grid, as will be discussed in more detail later). The formulation includes mechanisms to detect large interfaces and a compression term, similar to that used in other models of this kind [14,19,20]. The compression is entirely numerical, and maintains the interface sharp by counteracting the numerical diffusion that affects multifluid solvers. In addition, an extended blending method ensures that proper closures, still necessary to model transfer processes at the interface, are used as a function of the local flow regime. In this work, only adiabatic flows are considered and therefore, only mass and momentum conservation equations are solved:

$$\frac{\partial }{\partial t}\left({\alpha }_k{\rho }_k\right)+\nabla ·\left({\alpha }_k{\rho }_k{\mathit{U}}_k\right)+\nabla ·\left[{\alpha }_k\left(1-{\alpha }_k\right){\rho }_k{\mathit{U}}_c\right]=0.$$

(1)

$$\frac{\partial }{\partial t}\left({\alpha }_k{\rho }_k{\mathit{U}}_k\right)+\nabla ·\left({\alpha }_k{\rho }_k{\mathit{U}}_k{\mathit{U}}_k\right)=-{\alpha }_k\nabla p+\nabla ·\left[{\alpha }_k\left(\mathbf{\tau }+{\mathbf{\tau }}^{Re}\right)\right]+{\alpha }_k{\rho }_k\mathit{g}+{\mathit{M}}_k.$$

(2)

In the previous equations, α, ρ and U are the volume fraction, density and velocity of each phase, which is identified by the subscript k. p is the pressure, g the gravitational acceleration and τ and τ^Re the viscous and turbulent stress tensors. In Equation (1), U_c is the compression velocity that counteracts numerical diffusion in the presence of sharp gradients, such as that of the phase volume fraction across an interface. The compression velocity is directed normal to the interface and it is proportional to the relative velocity between the phases U_r:

$${\mathit{U}}_c=C_{\alpha }\left|{\mathit{U}}_r\right|\frac{\nabla \alpha }{\left|\nabla \alpha \right|}.$$

(3)

The activation of the compression term for large interfaces only is obtained through the C_α scalar field, which is equal to 1 in numerical cells where large interfaces are detected and 0 in dispersed regions or when continuous regions of the same phase exist. The value of C_α is established by a large interface detection algorithm, discussed later in Section 2.1.

In Equation (2), M_k is the interfacial momentum transport term and models the dynamic interaction between the phases as a sum of forces, each accounting for a different effect. In total, 6 forces are considered, those being the drag (F_d), lift (F_l), wall lubrication (F_w), virtual mass (F_vm), turbulent dispersion (F_td) and surface tension (F_st) forces:

$${\mathit{M}}_k={\mathit{F}}_d+{\mathit{F}}_l+{\mathit{F}}_w+{\mathit{F}}_{td}+{\mathit{F}}_{vm}+{\mathit{F}}_{st}.$$

(4)

In dispersed regions, all forces except surface tension are included. In large interface regions, instead, only drag, using a specific large interface model, and surface tension remain active. This is achieved by extending the blending capabilities existing in OpenFOAM. For a gas–liquid system, linear and hyperbolic functions allow the blending of the closures as a function of the volume fraction, from gas dispersed in liquid (e.g., bubbly flows) to liquid dispersed in gas (e.g., droplet flows). These functions are extended by joining them with the C_α field, effectively allowing the solver to locally apply the desired closures where large interfaces are present. For the drag force, which is active in all regimes, the blended formulation reads:

$${\mathit{F}}_d=\left[1-\left(1-C_{\alpha }\right)f_{gl}-\left(1-C_{\alpha }\right)f_{lg}\right]{\mathit{F}}_{d,LI}+\left(1-C_{\alpha }\right)f_{gl}{\mathit{F}}_{d,gl}+\left(1-C_{\alpha }\right)f_{lg}{\mathit{F}}_{d,lg}.$$

(5)

In Equation (5), LI identifies the large interface drag closure and gl and lg the gas-in-liquid and liquid-in-gas closures and blending functions f. In cells where a large interface is present, C_α is equal to 1 and the large interface closure is employed. In the other cells, C_α is 0 and the gas-in-liquid or liquid-in-gas closure is employed depending on the local volume fraction (more details on how this practically works are provided in Section 3.1.1, when the results for the channel flow from Fabre et al. [30] are discussed).

The surface tension is active only for large interfaces (C_α = 1):

$${\mathit{F}}_{st}=C_{\alpha }{\mathit{F}}_{st,LI}.$$

(6)

All the other forces are only active in gas-in-liquid and liquid-in-gas regions, and they are deactivated in the presence of large interfaces using C_α:

$${\mathit{F}}_{\left(l,w,vm,td\right)}=\left(1-C_{\alpha }\right)f_{gl}{\mathit{F}}_{gl}+\left(1-C_{\alpha }\right)f_{lg}{\mathit{F}}_{lg}.$$

(7)

For the modelling of the flow regime addressed in this work, only drag and surface tension are relevant and the specific closures are discussed in more detail in Section 2.2.

2.1. Large Interface Detection

The GEMMA solver’s local formulation is controlled by the value of C_α, which governs the interface compression and the blending of the interfacial closures. While in the native OpenFOAM formulation, and in most similar models available in the literature, closure blending is an a priori function of the volume fraction, in GEMMA it depends on the local morphology of the interface. This is achieved with a specific methodology to identify large interfaces in the multifluid field from the local length scale of the interface and the local resolution of the computational grid and to set the value of C_α accordingly. As mentioned before, in this modelling context, any reference to a “small” or a “large” interface needs to be related to the computational grid used to resolve it (e.g., a small bubble can be treated as a large interface if finely resolved, as performed in [21]). An additional advantage of the GEMMA approach is the capability of choosing which interface scales are resolved.

A first criterion verifies that a sufficient mesh resolution is available locally for the interface length scale to be resolved by calculating the local interface resolution quality (IRQ), as originally proposed in [19]:

$$IRQ=\frac{2}{\Delta \kappa }>IR{Q}_{crit}.$$

(8)

In Equation (8), Δ is the local mesh size and κ the interface curvature, and IRQ increases with the resolution of the interface. Resolution is acceptable if higher than a critical value, assumed in this work to be equal to 2. The same value has been already used in the validation of GEMMA against several test cases [21], and it is maintained in the present work. If the IRQ condition is met, two conditions can be used to evaluate the local length scale of the interface and effectively activate the large scale formulation of GEMMA. The first is based on the local volume fraction gradient being higher than a critical value γ [4,18]:

$$\frac{\left|\nabla \alpha \right|}{max\left(\left|\nabla \alpha \right|\right)}>\gamma .$$

(9)

Alternatively, when a population balance model is coupled with the multiphase solver, the large interface mode is activated from the local dispersed phase size:

$$d_{SM}>\Gamma \Delta ,$$

(10)

where d_SM is the Sauter-mean diameter, provided by the population balance model, and it needs to be larger than the local mesh resolution by at least a factor Γ. Finally, a last check is made to deactivate the large interface model in regions that are already continuous and where the void fraction is higher than 0.99 or lower than 0.01.

At the present time, GEMMA is already equipped to work with the class-method population balance model [32] available in the reactingEulerFOAM family of solvers. For the stratified flows modelled in this paper, however, a population balance model is not necessary and the volume fraction gradient criterion in Equation (9) is used, with a value of γ equal to 0.1.

2.2. Interfacial Momentum Transfer

From Equation (4), only drag and surface tension forces play a role in the flows modelled in this study. In dispersed regimes, drag is modelled with the Ishii and Zuber [33] correlation when gas is dispersed in the liquid flow (e.g., bubbly flow) and with the Schiller and Naumann [34] correlation when liquid is dispersed in the gas phase (e.g., droplet flow). Instead, for large interfaces, the segregated flow regimes drag model from Marschall [22], formulated as an interfacial force density from unbalanced pressures and viscous stresses, is used:

$${\mathit{F}}_{d,LI}=\left[0.5\frac{{\rho }_m\delta \left|{\mathit{U}}_g-{\mathit{U}}_l\right|}{\frac{{\alpha }_g{\alpha }_l{\mu }_g{\mu }_l}{{\mu }_g+{\mu }_l}}+8\frac{\frac{{\alpha }_g{\alpha }_l{\mu }_g{\mu }_l}{{\alpha }_g{\mu }_l+{\alpha }_l{\mu }_g}}{\frac{{\mu }_g{\mu }_l}{{\mu }_g+{\mu }_l}}\right]\frac{\left|\nabla \alpha \right|}{\delta }\frac{{\mu }_g{\mu }_l}{{\mu }_g+{\mu }_l}\left({\mathit{U}}_g-{\mathit{U}}_l\right).$$

(11)

In Equation (11), μ is the dynamic viscosity, ρ_m the volume-fraction-weighted mixture density and δ the width of the interface, modelled as $\delta =\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\left|\nabla \alpha \right|$}\right.$ [22].

For modelling the surface tension force, the continuum surface force method of Brackbill et al. [35] is employed, where the force per unit volume is given by the surface tension, the interface curvature and the gradient of the volume fraction. A density correction is included, given the large density ratio existing between the phases considered [21,36], and the volume fraction field used in the calculation of the curvature is smoothed by successive interpolation from cell centres to faces (smoothed continuum surface force from Ubbink [37]) to limit the appearance of parasitic currents. Finally, the surface tension force is calculated in each cell for each phase k by summation over all the n_k phases that share an interface with k in the cell and weighted using the cell phase volume fraction:

$${\mathit{F}}_{st,k}={\alpha }_k{\displaystyle \sum }_{i=1}^{n_k}{C}_{\alpha ,ki}{\sigma }_{ki}\kappa \nabla \alpha \frac{2{\rho }_m}{\Delta {\rho }_{k,i}}.$$

(12)

2.3. Turbulence Modelling

The last term in Equation (2) that requires modelling is the turbulent stress tensor. Here, a mixture k-ε RANS approach is used, which solves balance equations for the mixture turbulence kinetic energy k_m and the mixture turbulence kinetic energy dissipation rate ε_m [38]:

$$\frac{\partial }{\partial t}\left({\rho }_m{k}_m\right)+\nabla ·\left({\rho }_m{k}_m{\mathit{U}}_m\right)=\nabla ·\left({\mu }_m^{tot}\nabla k_m\right)+P_m-{\rho }_m{\epsilon }_m+S_{k,m}.$$

(13)

$$\frac{\partial }{\partial t}\left({\rho }_m{\epsilon }_m\right)+\nabla ·\left({\rho }_m{\epsilon }_m{\mathit{U}}_m\right)=\nabla ·\left({\mu }_m^{tot}\nabla {\epsilon }_m\right)+C_1{P}_m\frac{{\epsilon }_m}{k_m}-C_2{\rho }_m\frac{{\epsilon }_m}{k_m}{\epsilon }_m+S_{\epsilon ,m}.$$

(14)

In Equations (13) and (14), U_m is the mixture velocity, ρ_m is the mixture density, obtained from the weighted average of the phase densities, and P_m is the mixture turbulence production due to shear, obtained from the mass-weighted average of the phase values. The use of a mixture model improves the stability of the model by avoiding solving for phase-related turbulence quantities in regions where the phase volume fraction is negligible (i.e., where the other phase is continuous). In these regions, local extremes and strong oscillations in the production and dissipation of turbulence are often found, which, however, become negligible in the mixture balance with respect to the local continuous phase. In addition, in stratified flows, the mixture turbulence model essentially converges to the local phase turbulence model on each side of the interface.

In Equations (13) and (14), interfacial turbulence production and dissipation mechanisms are modelled with sources S_k,m and S_ε,m. The same blending method used for interface transport closures is extended to the modelling of interfacial turbulence, which can therefore include both dispersed regime and large interface mechanisms. In this work, a mechanism of turbulence suppression near large interfaces is included through an additional dissipation source. This is needed to model the wall-like behaviour of the interface and to avoid the turbulence increasing to unphysical levels. Specifically, the suppression model developed by Frederix et al. [23] is adapted to the mixture turbulence formulation, obtaining a single source for the mixture from the sources of the individual phases:

$$S_{\epsilon ,m}=C_{\alpha }{\displaystyle \sum }_k\left[C_2{\alpha }_k{\rho }_k{\left(\frac{{\nu }_k}{{\delta }^2}\right)}^2{k}_k\right].$$

(15)

C_α guarantees that the suppression, assumed symmetric between the phases, is only active near large interfaces. The parameter δ is fixed at 10⁻⁴, as suggested in [23]. As mentioned, while only turbulence suppression near large interfaces is considered, the model is built in a generalized fashion and is equipped to account, depending on the local flow regime, for additional physical effects such as bubble- or droplet-induced turbulence or the generation of turbulence by interfacial wakes.

3. Results and Discussion

3.1. The Experiment of Fabre et al.

In this section, the GEMMA solver is used to predict the stratified flow experiment in a rectangular channel of Fabre et al. [30]. In the experiment, air and water were injected in a 12 m long channel having a 0.1 m × 0.2 m cross-section (height × width). Two different cases, 250 and 400, were predicted, for which values of water and air flow rates and liquid height measured in the channel are reported in

Table 1. Summary of the conditions in the experiments of Fabre et al. [30] selected for comparison.

Experiment	Ql [ls−1]	Ul [ms−1]	Qg [ls−1]	Ug [ms−1]	hl [m]
250	3.0	0.395	45.4	3.66	0.038
400	3.0	0.476	75.4	5.50	0.0315

. For both cases, a generally smooth stratification was recorded, with only a mild wavy behaviour, more evident for the 400 case which was in the wavy flow regime.

For the 250 case, where three-dimensional effects were not expected, a two-dimensional geometry was employed in the predictions [23], with results recorded at the same 9.1 m distance from the inlet used in the experiment. The computational grid had 60 × 3600 elements, sufficient for obtaining mesh-independent results. Results from the mesh sensitivity study are provided in Appendix A. Inlet velocities match the experimental flow rates, except for the gas velocity for the 250 experiment, which was increased to 4.0 ms⁻¹ following previous authors [20,23], and the inlet phase distribution was fixed to match the measured liquid level. The no-slip condition was imposed on the upper and lower walls, where the law of the wall was enforced in the first near-wall cell by using the usual single-phase high-Reynolds-number wall function. A constant pressure condition was imposed on the outlet section. The use of a constant pressure profile was found to perturb the flow and the phase distribution upstream in the vicinity of the outlet section. However, the measurement point was at a sufficient distance upstream from the outlet section to avoid any distortion on the results.

For the 400 experiment, instead, a three-dimensional 40 × 60 × 1250 grid was used, for a total of 3,000,000 cells in the domain. The computational effort was limited by reducing the domain length to 7 m, using results from the 250 case that showed how a 3 m distance from the inlet section was sufficient to obtain, on average, fully developed conditions. For this case, results were recorded 4 m from the inlet section.

3.1.1. 250 Experiment

The two-dimensional phase distribution for the 250 experiment, showing a smooth stratified flow of air on top of water and only minimal perturbations at the interface, is shown in Figure 1a. The water level 9.1 m away from the inlet is shown in Figure 1b as a function of time. The value is obtained by averaging at each time step the liquid volume fraction over the height of the channel. The water level decreases for the first 10 s, after which a very mild increase is observed, with good agreement found with the experimental value of 0.038 m, shown in the plot as a constant line.

The two-dimensional distribution of C_α is shown in Figure 1c. The figure demonstrates how the GEMMA solver detects the presence of the large interface in the stratified regime, which leads to the activation of the large interface mode for the entire length of the channel, with C_α = 1. In Figure 1d, the value of the three components of the drag force in Equation (5) with the height in the channel (aligned with the y-axis in this horizontal channel as well as in the annular channel in Section 3.2) is shown to illustrate the blending closure mechanism and the drag model behaviour. In the continuous region, the gas-in-liquid and liquid-in-gas components are active, although their contribution is negligible given the negligible amount of dispersed phase. Near the interface (where C_α = 1), they are deactivated and the large interface model is the only contributor to the interfacial drag.

Figure 2 shows profiles of velocity and turbulence quantities along the height of the channel compared against the experimental measurements of Fabre et al. [30]. The subscript m identifies a mixture weighted average of phase values for velocity, turbulence kinetic energy and turbulence energy dissipation rate. In Figure 2a,b, results without the suppression term in Equation (15) are also included to illustrate the impact this has on the predictive capabilities of the model. When Equation (15) is included, the air phase exhibits an almost parabolic profile, except for a slight shift towards the wall with respect to the measurements. Without turbulence suppression, in contrast, due to the turbulence kinetic energy being largely overestimated near the interface, the air velocity profile is drastically different. The much higher turbulent mixing reduces the velocity gradient near the interface but also greatly extends the region of intense momentum transport above it, where the velocity increases at an almost linear rate. Due to this, and the consequent higher flow resistance in the area, the air velocity is reduced in the air domain above the interface and the velocity peak is shifted significantly towards the upper wall.

With turbulence suppression, the velocity remains in good agreement with the experiment in both water and air regions, except for the mentioned slight shift of the air velocity profile towards the upper wall. Prediction of turbulence kinetic energy is also greatly improved in Figure 2b. The large difference in values in the figure prevents a detailed understanding, therefore a closer look at the results obtained with turbulence suppression in the water and air regions is given in Figure 2c,d. In the water, the turbulence kinetic energy decreases from the lower wall towards the interface (Figure 2c). In the air, the turbulence kinetic energy also decreases from the lower wall towards the interface, goes through a minimum, after which it increases again as the interface is approached. These trends are correctly predicted by the CFD model, although discrepancies still remain. At the interface on the air side, and despite the turbulence suppression model, the turbulence kinetic energy remains overestimated. On the water side of the interface, instead, the kinetic energy drops to a negligible value, much lower than in the experiments. The same behaviour has been recorded by other authors, and some solutions [39,40], such as the use of an asymmetric suppression model [40] limited to the air phase, are being investigated and will provide useful input to future work aimed at further improving the GEMMA approach. Turbulence also remains underpredicted near the lower wall, but this is an expected consequence of the use of wall functions.

Figure 2e,f display the turbulence energy dissipation rate and the turbulent viscosity, which again demonstrate how much lower values are obtained when turbulence suppression is included in the model. The turbulent viscosity, as shown in Figure 2f, is reduced by more than an order of magnitude. An order of magnitude difference is also seen between values in air and water, which are visualized better in the window on the top right of the figure, because of the much larger velocity and turbulence kinetic energy that air possesses. The much higher turbulent viscosity in the absence of suppression promotes turbulent mixing near the interface. The perturbation induced by the presence of the interface penetrates much further in the air region, where the velocity increases almost linearly away from the interface and the velocity peak moves much closer to the upper wall to enforce mass conservation (Figure 2a). When suppression is included in the model formulation, on the other hand, the turbulent viscosity and the turbulent momentum transport are drastically reduced, allowing the model to recover the wall-like behaviour of the interface and its role in suppressing the wall normal turbulent fluctuations [24,25].

From Fabre et al. [30], measurements are also available for pressure drop, wall shear stresses for the liquid and gas and the interfacial shear stress. These are reported in

Table 2. Pressure drop and wall and interfacial shear stresses from the CFD simulations and the experiments of Fabre et al. [30], with percentage relative errors from GEMMA in brackets.

	τw,l [Pa]	τw,g [Pa]	τi [Pa]	dp/dx [Pam−1]
GEMMA 250	0.392 (−12.7%)	0.055 (−5.2%)	0.112 (93%)	−2.70 (22.2%)
Fabre et al. [30] 250	0.449	0.058	0.058	−2.10
GEMMA 400	0.634 (−21.9%)	0.162 (23.7%)	0.343 (28.9%)	−9.0 (34.3%)
Fabre et al. [30] 400	0.812	0.131	0.266	−6.7

together with the values predicted with GEMMA. The wall shear stresses are calculated from the turbulent viscosity and the derivative of the velocity in the first near-wall cell:

$${\tau }_w=-{\left(\nu +{\nu }_t\right)}_k{\left(\frac{dU}{dy}\right)}_w.$$

(16)

The interfacial shear can be estimated from the momentum equation, which, after assuming the flow to be one-dimensional, steady and fully developed, can be written for the liquid or the gas as [41]:

$$-A_{l-g}{\left(\frac{dp}{dx}\right)}_{l-g}-{\tau }_{w,l-g}{S}_{w,l-g}\pm {\tau }_i{S}_i=0,$$

(17)

Here, A is the flow area of each phase and S the surface area per unit length of the wall and the interface. Under the assumptions made for the flow, the pressure gradient in the two phases is equal and we can focus on the gas phase only. Using the pressure gradient and the wall shear stress from the CFD predictions, the interfacial shear can be estimated from:

$${\tau }_i=\frac{-A_g{\left(\frac{dp}{dx}\right)}_g-{\tau }_{w,g}{S}_{w,g}}{S_i}$$

(18)

For the estimation of the area of the interface, a flat smooth surface is assumed. From Table 2, while for the 250 experiment the wall shear stress shows good agreement with the experimental value, the interfacial shear is instead largely overestimated. Even so, the pressure drop remains closer to its experimental counterpart. Additional and more detailed studies will be required, even though the larger pressure gradient obtained with a lower wall shear stress supports the overestimation in the interfacial shear in the CFD results. Going back to the momentum balance in Equation (2), the interfacial shear is modelled with the interfacial momentum transfer source and this is essentially provided by the large interface drag, which is therefore the model that will require additional research moving forward. Overall, the pressure drops are more difficult to predict due to the presence of the interface and its influence on the velocity profiles in the gas and liquid regions, which inevitably also affects the prediction of wall shear stresses in single-phase regions near the walls.

3.1.2. 400 Experiment

The 400 experiment case has the same liquid flow rate and a higher gas flow rate, therefore a slightly stronger wave pattern is expected at the interface, although still mild and not sufficient to cause any entrainment of liquid into the gas or vice versa. In the predicted three-dimensional void distribution along the channel in Figure 3a, a periodic wavy behaviour at the interface is apparent. Consistently with the experiments, this is only mild and does not cause any significant perturbation of the continuous interface structure. In Figure 3b, the C_α field is displayed, confirming the proper detection of the interface region from the GEMMA solver. Essentially, C_α tracks the interface surface. The vertical location of the interface, and consequently the liquid height, remains almost constant along the channel. This is confirmed in Figure 4 by the liquid height as a function of time taken 4 m from the inlet. The predicted liquid level remains almost constant at the inlet value and it is in very good agreement with the experiment. The figure further confirms the very mild wavy behaviour observed, also in agreement with Fabre et al. [30].

In Figure 5, results similar to those provided for case 250 are shown on a vertical line through the channel and 4 m from the inlet. Figure 5a shows the velocity profile. Qualitatively, agreement is found and the correct shape of the velocity profile is predicted in the air region. In the water region, good quantitative agreement is also found, although in the air region the velocity is overestimated. Given that the correct shape of the velocity profile is recovered, and the liquid height is in good agreement although the CFD result is everywhere higher than the experimental value, the overestimation is likely related to the provided experimental air flow rate that has been used to calculate the inlet velocity.

The turbulence kinetic energy is reported in Figure 5b for the entire channel, with details of the water and air regions given in Figure 5c,d. Results are similar to those of the 250 case, although the turbulence is underestimated in the air region. The shape of the kinetic energy profile is well predicted and the suppression sufficient to avoid overestimating the turbulence near the interface. At the walls, the underestimation of the turbulence kinetic energy can be attributed to the use of wall functions. On the water side, the marked drop near the interface occurs again, although it is milder with respect to the 250 case. The value of delta equal to 10⁻⁴, derived by [23], is again used based on the 250 case. Therefore, further investigation should aim at establishing if prediction on the air side can be improved by using a different value or a more physically based model.

3.2. Slug Flow in an Annular Channel

After the horizontal channel of Fabre et al. [30], GEMMA is here applied to a different configuration and flow regime, a horizontal annulus where a slug flow regime is established. The condition is of interest in the oil and gas sector, where the annular configuration is found in multiple equipment types where the periodic formation of slugs can lead to operational issues and deterioration of equipment [7,8]. For this reason, it has been the focus of multiple experimental and numerical studies. Friedemann et al. [7] predicted with the interFoam VOF solver of OpenFOAM the experiments made at the IET in Norway. Here, we predict the same experiment with GEMMA, with the aim of evaluating its capability of predicting the slug flow with respect to an interface resolving modelling technique. If accurate, GEMMA may potentially enable the use of less refined meshes, provided it is able to efficiently model in the multifluid field bubbles that are entrained in the liquid, and liquid droplets that are dispersed in the gas, around the boundaries of the liquid slugs near the upper wall of the pipe. In contrast, the modelling of these same structures would require an immense computational effort with interface-resolving methods, where all interfacial details need to be resolved.

In the experiment at IET, a stratified flow of oil and gas was studied in a very long horizontal annulus, with 0.05 m inner diameter and 0.1 m outer diameter. Measurements were taken at a very long distance from the inlet, from 36 m to 39 m. Therefore, as also performed in Friedemann et al. [7], we only simulated a small portion of the annulus and apply periodic conditions. Velocities and turbulence were mapped from the outlet to the inlet and the gas volume fraction is maintained at its initial inlet value, taken as 0.47 from the experiment. The mixture velocity was 3.6 ms⁻¹, the density and kinematic viscosity were 801 kg m⁻³ and 2·10⁻⁶ m²s⁻¹ for the oil and 23.82 kg m⁻³ and 7.56·10⁻⁷ m²s⁻¹ for the air, and the surface tension as 0.0285 N m⁻¹. Friedemann et al. [7] tested different pipe lengths, finding both 3 m and 7 m pipe sections to provide good results. Therefore, we adopted the 3 m test section in this work. Only a half section of the pipe was considered, and symmetry conditions were applied at the upper and lower lateral surfaces. The flow was initialized as stratified, with the volume of air occupying 47% of the volume, and the simulation was run for 60 s. For comparison, two meshes were built with a number of elements in line with the “extremely coarse” and the “very fine” meshes from [7]. The coarse mesh had 116,100 elements (30 × 18 × 215, circumference × radius × length) and the fine mesh 430,080 (48 × 28 × 320).

Figure 6 shows the gas volume fraction at the beginning of the simulation and after 15 and 45 s, with the coarse mesh. At the same three times, the distribution of C_α is also displayed. At the beginning of the simulation, the flow is stratified and, as the simulation progresses, instability develops and liquid slugs, almost reaching the upper portion of the annulus wall, appear and travel along the channel. The flow is characterized by the periodic build-up and dissipation of slug structures (two of these are present at t = 15 s in Figure 6b, but only one at t = 45 s in Figure 6c). In the build-up phase, disturbances on the interface grow in size, travel faster than the surrounding liquid and accumulate further liquid by incorporating downstream waves and disturbances. Liquid height can reach the top of the channel, and the upper portion of the slug is generally perturbed by gas entrained in the liquid slug or small liquid structures dispersed from the liquid main body. The breaking of the liquid front on top of the slug can lead to significant gas entrainment and the eventual dissipation of the slug structures (two of these breaking structures are visible in Figure 6b. The flow direction in the figure is from right to left) [7]. The distribution of C_α highlights the interface between the oil and the gas. This is correctly identified when it is almost flat at the beginning of the simulation but also when it is much more perturbed in the presence of slugs. In both regimes, the large interface mode gets activated in the interface region. Similar qualitative findings are observed in results obtained with the fine mesh.

Figure 7 shows the time behaviour of the oil volume fraction averaged over the annulus cross-section, which allows for the detection of the slug frequency in the channel and comparison with the results of Friedemann et al. [7]. The figure reports the coarse mesh on the left and the fine mesh results on the right. Qualitatively, the behaviour is similar to the results in [7]. The number of slugs passing through the channel is similar, although they are characterized by a generally lower oil volume fraction, and values higher than 0.8 or lower than 0.4 are rarely observed. In Friedemann et al. [7], higher and lower thresholds for the detection of slugs are 0.7 and 0.4, and a slug frequency of 1.3–1.5 Hz is predicted (depending on the channel length and mesh), which is close to the experimental value of 1.43. Here, comparable values around 1.3 can only be obtained by reducing the higher and increasing the lower threshold for slug detection to 0.65 and 0.5, due to the more limited excursion of the oil volume fraction associated with the slugs. Therefore, the mechanism leading to slug formation is predicted by GEMMA, the slug interface is detected and the large interface mode applied when a slug travels in the annulus. However, slugs have a smaller liquid volume fraction with respect to the VOF results in Friedemann et al. [7] and the oil essentially never occupies the entire channel, with no significant differences between the two meshes (with the same issue, although more limited, also noticed in [7]). Probably, this is associated with the accuracy of the large interface drag model and the overall large interface momentum transport modelling, which will be targeted by specific developments moving forward, as also suggested by the interfacial shear stress results in the stratified channel flow. Other possible explanations for the discrepancies in the liquid level excursion can be found in the modelling of the dispersed regions and entrained gas in the upper portion of the slugs and the impact of the length of the computational domain, which will be also investigated in future studies.

4. Conclusions

In this work, we have further developed our all flow regime approach GEMMA to predict horizontal segregated flow regimes, specifically the stratified channel flow of Fabre et al. [30] and the slug flow in a horizontal annulus studied experimentally at IET in Norway and computationally by Friedemann et al. [7]. The stratified flow confirms the capability of the model to detect the large interface and activate the large interface mode that maintains it sharp, essentially realizing a tracking of the interfacial surface in the multifluid field. The flow regime is well predicted and only a mild wavy behaviour perturbs the interface at the highest flow rate (400 experiment). Essential for the correct prediction of the fluid dynamics of stratified flows is the damping of the turbulence at the interface, similar to that near a solid wall. It is modelled with a specific source of turbulence energy dissipation, active only in large interface regions, that allows for the correct prediction of velocity and turbulence energy profiles and the liquid level, in both the 250 and the 400 experiments. Further improvements are envisaged, mostly in terms of correcting the excessive drop in turbulence kinetic energy on the liquid side of the interface. Additional model validation will be also pursued using available datasets for co-current and counter-current stratified flows in horizontal channels [5,28]. In the slug flow regime, the formation of liquid slugs is well predicted by GEMMA, and the large interfaces are again detected and followed as they are deformed and perturbed by the slugs travelling along the annular channel. The slug frequency measured in the experiments can be reproduced by GEMMA, but with a more limited excursion of the liquid volume fraction between its maximum and minimum values with respect to both the experiments and the computational modelling with VOF [7].

An analysis of the pressure drop in the stratified flows has highlighted that further improvement is possible to the large interface drag, which is responsible for the modelling of the interfacial shear in the flow. Although the pressure drop predicted is not largely different to the experimental measurements, because of the large contribution from the wall shear stresses that are well predicted, the analysis highlights a significant overestimation of the interfacial shear stress. Therefore, focused development will target the large interface drag model, which should also help to improve the quantitative prediction of the liquid volume fraction excursion during slug flow. Overall, the results confirm that the GEMMA solver is a credible all flow regime modelling concept, aiming at providing improved modelling capabilities of the thermal hydraulics of the loss of coolant accidents in nuclear reactors and the detection of slug flow formation in the transport of multiphase mixtures. This work has focused on horizontal segregated flow regimes, and further developments, to take full advantage of the approach, will extend the modelling to the break-up of the large interface with entrainment of droplets and bubbles and multiscale transfer between the averaged and resolved interfaces.