Charles Demay
email: charles.demay@edf.fr

Jean-Marc Hérard
email: jean-marc.herard@edf.fr

A compressible two-layer model for transient gas-liquid flows in pipes

Keywords: Two-phase flow, Two-layer model, Pipe flow, Variable cross section, Hyperbolic system, Entropy

HAL is

Introduction

The present work focuses on the modeling of transient gas-liquid flows in pipes, especially air-water flows. This type of flow occurs in piping systems of several industrial areas such as nuclear power plants, petroleum industries or sewage pipelines. The presence of air in such facilities is usually unwanted as it may induce for instance pressure surges and traveling air pockets, leading to reduced efficiency and damages for pumping systems, see [START_REF] Murakami | Effects of entrained air on the performance of a horizontal axial-flow pump[END_REF][START_REF] Poullikkas | Effects of two-phase liquid-gas flow on the performance of nuclear reactor cooling pumps[END_REF] ; an accurate modeling is thus necessary to improve performances and reliability. Although being studied since many years, the macroscopic description of those flows is still complex to assess as they may display different regimes such as dispersed flow, stratified flow, pressurized flow (pipe full of water), slug flow and transitions between them, see [START_REF] Ishii | Thermo-fluid dynamic theory of two-phase flow[END_REF] for details. On the one hand, experimental studies have provided many interesting data on key points such as transition from stratified to pressurized regime [START_REF] Capart | Numerical and experimental water transients in sewer pipes[END_REF] and air entrainment [START_REF] Chosie | Experimental and numerical investigation on the motion of discrete air pockets in pressurized water flows[END_REF][START_REF] Escarameia | Investigating hydraulic removal of air from water pipelines[END_REF][START_REF] Pothof | Experimental study of air-water flow in downward sloping pipes[END_REF][START_REF] Pozos | Air entrapped in gravity pipeline systems[END_REF]. On the other hand, mathematical and numerical modeling still raise many challenges as pointed out in the recent literature review [START_REF] Bousso | Numerical modeling of mixed flows in storm water systems: critical review of literature[END_REF]. In particular, the free-surface regime is usually described by an incompressible flow with surface waves, while the pressurized regime is described by a compressible flow with acoustic waves. Thus, most of existing 1D models focus on this transition without computing the air phase which increases the modeling difficulties, see [START_REF] Bourdarias | A finite volume scheme for a model coupling free surface and pressurised flows in pipes[END_REF]. However, air-water interactions may greatly affect the flow behavior regarding hydraulic jumps or traveling air pockets. Thus, the purpose of the model presented herein is to account for air-water interactions in addition to regime transitions occurring in air-water pipe flows.

As a first step, we consider in this paper the 1D modeling of stratified air-water flows in pipes. Seeing this regime as a freesurface flow, a common way to deal with it is to use the Saint-Venant system [START_REF] De Saint-Venant | Théorie du mouvement non-permanent des eaux avec application aux crues des rivières et à l'introduction des marées dans leur lit[END_REF], also called shallow water equations, which results from a depth averaging process on the Euler set of equations and assumes a thin layer of incompressible fluid (liquid phase) with hydrostatic pressure law. It may also include friction and viscosity, see [START_REF] Gerbeau | Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation[END_REF]. Following this classical approach, the model proposed in [START_REF] Bourdarias | A finite volume scheme for a model coupling free surface and pressurised flows in pipes[END_REF] computes the water layer, and strict hyperbolicity is acquired with mathematical entropy. In addition, the transition to the pressurized regime is handled including a switching from an incompressible to a compressible description for water and a discontinuity of the pressure gradient. With the aim of computing the air phase, one may rather consider the extension of the Saint-Venant system to a multilayer system, one layer being associated with one height, one velocity and constant density, see [START_REF] Audusse | A multilayer Saint-Venant system : Derivation and numerical validation[END_REF]. The air layer is thus added in [START_REF] Bourdarias | Air entrainment in transient flows in closed water pipes: a two-layer approach[END_REF] with variable density and perfect gas pressure law, the compressibility of air playing a key role in our framework. This process results in a two-layer model with incompressible liquid and compressible gas, referred hereafter as the incompressible/compressible two-layer model. Nonetheless, the latter is hyperbolic only for small or large relative speed between both phases. In fact, this system inherits the difficulties from the common incompressible/incompressible two-layer systems, such as non-conservative terms, non-explicit eigenstructure and conditional hyperbolicity, see [START_REF] Abgrall | Two-layer shallow water system: a relaxation approach[END_REF][START_REF] Bouchut | An entropy satisfying scheme for two-layer shallow water equations with uncoupled treatment[END_REF]. Note that a multilayer Saint-Venant system which accounts for mass exchanges between the layers is proposed in [START_REF] Audusse | A multilayer Saint-Venant system with mass exchanges for shallow water flows. Derivation and numerical validation[END_REF] and is strictly hyperbolic in its two-layer version when the total water height is strictly positive. Dealing with two-fluid flows, the homogeneous model [START_REF] Chaudhry | Analysis of transient pressures in bubbly, homogeneous, gas-liquid mixtures[END_REF] or the drift-flux model [START_REF] Faille | A rough finite volume scheme for modeling two-phase flow in a pipeline[END_REF] may also be proposed as candidates. However, the first one considers one velocity field for the two phases, while the second one imposes velocity profiles for the relative motion. Thus, those models are unable to restore the complexities of the flow such as entrapped air pockets with unpredictable velocity. Another approach may be the two-velocity two-pressure models where compressibility is assumed for both phases with an associated barotropic pressure law. In that context, interesting mathematical properties are obtained such as hyperbolicity as well as entropy inequality. It was first introduced for separated flows with depth averaging in [START_REF] Ransom | Hyperbolic two-pressure models for two-phase flow[END_REF] and mainly used afterward in a statistical framework for bubbly or granular flows, see [START_REF] Baer | A two phase mixture theory for the deflagration to detonation (DDT) transition in reactive granular materials[END_REF][START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF][START_REF] Glimm | Two phase flow modelling of a fluid mixing layer[END_REF][START_REF] Kapila | Two-phase modeling of DDT: structure of the velocity-relaxation zone[END_REF]. This class of models can also been obtained from a variational approach as in [START_REF] Gavrilyuk | Mathematical and numerical modeling of two-phase compressible flows with micro-inertia[END_REF] and may be extended to multi-component fluids, see [START_REF] Müller | Closure conditions for non-equilibrium multi-component models[END_REF]. Nonetheless, the statistical framework is not relevant for stratified flows as the height and the relative position of each phase are not given by the model.

Combining the interesting properties of the two-layer and the two-velocity two-pressure frameworks, the model proposed herein is a compressible/compressible two-layer model as in [START_REF] Ransom | Hyperbolic two-pressure models for two-phase flow[END_REF]. Thus, water is considered compressible and does not follow the classical hydrostatic pressure law but a barotropic pressure law, such as stiffened gas law. The starting point is the isentropic Euler set of equations for both phases where the classical hydrostatic assumption is made for water (vertical acceleration is neglected). Adding a kinetic boundary condition at the interface between the liquid and gas layers, an averaging process results in a five-equation model in which the hydrostatic constraint is used to define the interfacial pressure. Following an entropy inequality, closure laws are provided for the interfacial velocity and source terms such as mass and momentum transfer. Contrary to the incompressible/incompressible or incompressible/compressible two-layer frameworks, the eigenstructure can be easily detailed and the hyperbolicity is acquired except for resonance conditions which may not occur for realistic air-water flows. Regarding the Riemann problem associated with the homogeneous problem (i.e., without any source terms), the nature of characteristic fields and the associated Riemann invariants can also be detailed. Furthermore, since non-conservative terms exist, the uniqueness of jump conditions is an important feature that is studied and acquired. Thus, one may build analytical solutions for the Riemann problem. In addition, positivity is guaranteed for heights and densities. Finally, note that when dealing with regime transitions, this framework includes a uniform description of water as a compressible flow without any switching on the pressure law.

The document is organized as follows. For the sake of clarity, the model is presented first in the 2D framework with a 1D averaging process where a vertical depth averaging process is used in the classical way of shallow water two-layer equations. Thus, plane channels with constant width are implicitly considered in Sect. 2 neglecting the spanwise variations of the flow. Closure laws for the interfacial velocity and source terms are provided in Sect. 3 and the resulting closed system is commented in Sect. 4 regarding its consistency with other well-known models. Its ability to deal with more complex pipes configurations such as sloping pipes or pressurized flows is also studied. section 5 details significant mathematical properties of the model. In a second step, the 3D framework with a 2D averaging process is handled in Sect. 6 with the aim of considering circular pipes with variable cross section in time and space. The derivation is done following the same approach as in Sect. 2. Air-water flows are considered throughout the paper but the model can applied in the general framework of gas-liquid flows and vapor-liquid flows.

Model development

Local governing equations

Throughout this section, a two-layer air-water flow through an horizontal pipe of height H is considered. Index k is attributed to each phase, 1 for water, 2 for air. Thus, as h k refers to the height of phase k, one has h 1 + h 2 = H. A 2D description of the flow is used such that all local variables depend on (x, z,t) excepting h k (x,t). In this framework, we implicitly deal with rectangular channels homogeneous in the spanwise direction. The extension to circular pipes with variable cross section will be investigated in Sect. 6. The geometric description is given in cartesian coordinates (0, x, z) in Fig. 2.1. The set of local governing equations is given by the 2D isentropic Euler set of equations for both phases. In that context, both phases are assumed compressible and the pressure P k of phase k depends on density ρ k . Denoting u k = (u k , w k ) the velocity vector of phase k, mass and momentum conservation equations write:

∂ ρ k ∂t + div(ρ k u k ) = M k , k = 1, 2, (2.1)
∂ ρ k u k ∂t + div(ρ k u k ⊗ u k + P k I d ) = D k + ρ k g, k = 1, 2, (2.2)
where M k and D k = D k e x stand, respectively, for mass and momentum transfer between phases, g = -ge z denotes the gravity field and I d is the identity matrix. Note that mass transfer terms are introduced here for generality, considering singlecomponent vapor-liquid flows for instance, but may be neglected for most of air-water flows. In addition, mass and momentum conservation of the mixture yields

M 1 + M 2 = 0 and D 1 + D 2 = 0.
Dealing with stratified flows in pipes, the classical hydrostatic assumption is made for water (i.e., vertical acceleration is neglected), so that (2.2) along z for k = 1 yields:

∂ P 1 (ρ 1 ) ∂ z = -ρ 1 g, (2.3)
which will be referred to as the hydrostatic constraint. Note that in the framework of barotropic pressure laws, this constraint implies that the density ρ 1 (x, z,t) is not homogeneous along the vertical coordinate z.

As the phases are compressible, state equations are required for pressures. For instance, perfect gas law may be used for air and isentropic stiffened gas law for water:

P 1 (ρ 1 ) = (P 1,ref + Π 1 ) ρ 1 ρ 1,ref γ 1 -Π 1 , P 2 (ρ 2 ) = P 2,ref ρ 2 ρ 2,ref γ 2
, with some reference density ρ k,ref and pressure P k,ref . Classical thermodynamic considerations impose γ k > 1, Π 1 ≥ 0 and P k (ρ k ) > 0 so that c k = P k (ρ k ) defines the celerity of acoustic waves. For air, γ 2 is set to 7/5 (diatomic gas) while γ 1 and Π 1 are fitted according to a reference state. Note that the development herein is independent of the chosen pressure laws. An additional interfacial kinetic boundary condition may be added on the interface separating air and water writing:

ϕ(x, z,t) = z -h 1 (x,t) = 0, so that: dϕ dt = 0,
and thus:

∂ ϕ ∂t + u I .∇ x,z ϕ = 0,
where u I = U I W I denotes the interfacial velocity. Using the expression of ϕ, it finally writes:

∂ h 1 ∂t +U I ∂ h 1 ∂ x = W I . (2.4)
Boundary conditions may be applied to both phases. On the walls, impermeability imposes:

w 1 (x, z = 0,t) = w 2 (x, z = H,t) = 0. (2.5)
On the interface, the continuity of the normal velocity is assumed:

(u 1 (x, z = h - 1 ,t) -u 2 (x, z = h + 1 ,t)).n I = 0, (2.6)
with n I = (-

∂ h 1 ∂ x 1) T 1+( ∂ h 1 ∂ x ) 2
, the unit normal vector at the interface. As a dynamic boundary condition, a pressure equilibrium is imposed on the interface neglecting surface tension effects:

P 1 (x, z = h - 1 ,t) = P 2 (x, z = h + 1 ,t) = P I (x,t), (2.7)
where P I denotes the interfacial pressure. Finally, we assume that the interfacial velocity satisfies:

u I (x,t) = β u 1 (x, z = h - 1 ,t) + (1 -β )u 2 (x, z = h + 1 ,t), β ∈ [0, 1]. (2.8)

Averaging process

Following the shallow water modeling approach, a vertical average across the layer depth is performed using the operator:

f k (x,t) = 1 h k h k +z k z k f k (x, z,t)dz with z 1 = 0, z 2 = h 1 , (2.9)
where f k is a function depending on the state variables. The density weighted averaging operator is also introduced:

f k (x,t) = ρ k f k ρ k . (2.10)
In addition, the interface value of each variable is denoted f * k such that:

f * 1 (x,t) = f 1 (x, z = h - 1 ,t), f * 2 (x,t) = f 2 (x, z = h + 1 ,t).
(2.11)

Mass conservation

Vertical integration of (2.1) gives:

h k +z k z k ∂ ρ k ∂t dz + h k +z k z k ∂ ρ k u k ∂ x dz + h k +z k z k ∂ ρ k w k ∂ z dz = h k +z k z k M k dz.
Using Leibniz's integral rule for non-constant bounds in the first two terms of the left-hand side and kinetic boundary condition (2.5) for the third one, one gets:

∂ h k ρ k ∂t + ∂ h k ρ k u k ∂ x -B k = h k M k ,
with:

B k = ρ * k ( ∂ h 1 ∂t + u * k ∂ h 1 ∂ x -w * k ) = ρ * k ((u * k -U I ) ∂ h 1 ∂ x + (W I -w * k )) according to (2.4) = ρ * k (u I -u * k ).(-∂ h 1 ∂ x 1) T = 0,
(2.12) using (2.6) and the colinearity between u I -u * k and u * 1 -u * 2 given by (2.8). Finally, using the definition ρ k u k = ρ k u k , the averaged mass conservation equations write:

∂ h k ρ k ∂t + ∂ h k ρ k u k ∂ x = h k M k . (2.13)

Momentum conservation

Dealing with the streamwise component of (2.2), vertical integration writes:

h k +z k z k ∂ ρ k u k ∂t dz + h k +z k z k ∂ ρ k u 2 k ∂ x dz + h k +z k z k ∂ ρ k u k w k ∂ z dz + h k +z k z k ∂ P k ∂ x dz = h k +z k z k D k dz.
Using Leibniz's integral rule and kinetic boundary condition (2.5), the first three terms of the left-hand side write:

∂ h k ρ k u k ∂t + ∂ h k ρ k u 2 k ∂ x -u * k B k ,
where B k = 0 according to (2.12). Concerning the pressure gradient term, it comes:

h k +z k z k ∂ P k ∂ x dz = ∂ h k P k ∂ x -P I ∂ h k ∂ x .
At first order, the turbulence of the flow is neglected, such that the closure law ρ k u 2 k = ρ k u 2 k is chosen. Thus, the averaged momentum equations considered are:

∂ h k ρ k u k ∂t + ∂ h k (ρ k u 2 k + P k ) ∂ x -P I ∂ h k ∂ x = h k D k . (2.14)

Hydrostatic constraint

In [START_REF] Ransom | Hyperbolic two-pressure models for two-phase flow[END_REF], the averaged system for 1D stratified flows gathers the averaged mass conservation equations (2.13), the averaged momentum equations (2.14) and the kinetic boundary condition (2.4). However, as gravitational effects have a leading role in the dynamics of stratified flows, it is proposed here to account for the so-called hydrostatic constraint (2.3) in the averaged system. Some approaches have already been studied in the literature, see [20,[START_REF] Morales De Luna | A Saint Venant model for gravity driven shallow water flows with variable density and compressibility effects[END_REF], but an original one is proposed below. The starting point is to integrate the hydrostatic constraint (2.3) between z and h 1 :

P 1 (ρ 1 (x, z,t)) = P I (x,t) + h 1 z ρ 1 (x, s,t)gds. (2.15)
Performing a second integration between 0 and h 1 it yields:

h 1 0 P 1 (ρ 1 (x, z,t))dz = h 1 P I (x,t) + g h 1 0 h 1 z ρ 1 (x, s,t)ds dz. (2.16)
Then, using an integration by parts of the double integral it comes:

h 1 0 P 1 (ρ 1 (x, z,t))dz = h 1 P I (x,t) + g z h 1 z ρ 1 (x, s,t)ds h 1 0 + g h 1 0 ρ 1 (x, z,t)zdz, = h 1 P I (x,t) + ρ 1 zgh 1 .
At this point, it is proposed to neglect the correlation between ρ 1 and z, such that the closure law for ρ 1 z writes:

ρ 1 z = ρ 1 z = ρ 1 h 1 2 .
(2.17)

Remark 2.1 As detailed in Appendix A, the validity of (2.17) can be clarified when choosing a linear pressure law for phase 1. Indeed, assuming that ε

= gh 1 c 2 1
1, one obtains the following estimate:

ρ 1 z -ρ 1 z ρ 1 z = ε 6 + O(ε 2 ),
which justifies (2.17) as in practice, c 1 ∼ 1500 m.s -1 , H ∼ 1 m and |ε| ≤ 10 -5 .

Thus, (2.16) yields:

P I = P 1 -ρ 1 g h 1 2 . (2.18)
In the isentropic framework, the above expression provides the closure law for P I regarding P 1 as a function of ρ 1 . Moreover, (2.18) may be read as an averaged hydrostatic pressure law for P 1 used with the averaged density ρ 1 . Therefore, the averaged momentum equations under hydrostatic constraint write:

∂ h k ρ k u k ∂t + ∂ h k (ρ k u 2 k + P k ) ∂ x -(P 1 -ρ 1 g h 1 2 ) ∂ h k ∂ x = h k D k . (2.19)

Resulting averaged system

Adding the interfacial kinetic boundary condition (2.4), a five-equation system corresponding to the five unknowns

(h 1 , ρ 1 , ρ 2 , u 1 , u 2 ) is obtained: ∂ h 1 ∂t +U I ∂ h 1 ∂ x = W I , (2.20a)
∂ h k ρ k ∂t + ∂ h k ρ k u k ∂ x = h k M k , (2.20b)
∂ h k ρ k u k ∂t + ∂ h k (ρ k u 2 k + P k ) ∂ x -(P 1 -ρ 1 g h 1 2 ) ∂ h k ∂ x = h k D k , (2.20c)
where k = 1, 2 and h 1 + h 2 = H. This system is relevant for two-layer gas-liquid or vapor-liquid flows in horizontal channels and relies on a vertical average across the layer depth. Compressibility is considered for both layers with barotropic pressure laws. Dealing with the averaged source terms, conservation of mass and momentum of the mixture impose:

h 1 D 1 + h 2 D 2 = 0, (2.21a
)

h 1 M 1 + h 2 M 2 = 0. (2.21b)
In practice, the averaged pressure laws P k (ρ k ) are classically expressed in terms of ρ k so that the closure P k (ρ k ) = P k (ρ k ) is chosen. Consequently, the celerity of acoustic waves in the averaged framework is defined by:

c k = P k (ρ k ).
(2.22)

The interfacial pressure P I is defined by the averaged hydrostatic constraint (2.18) whereas the interfacial velocity U I and source terms still need closure laws. In the following section, an entropy inequality is used to close the system.

Entropy inequality and closure laws

System (2.20) may develop discontinuous solutions even with continuous initial conditions. Therefore, an entropy inequality which allows to select the physically relevant solution is obtained in this section; closure laws will be proposed following with this inequality.

From now on, the operator notations are omitted until Sect. 6. Let us define:

m k = h k ρ k , E c,k = 1 2 m k u 2 k , E t,k = m k Ψ k (ρ k ), (3.1)
with Ψ k some function of ρ k . Notice that E c,k and E t,k respectively stand for the integrated kinetic energy and thermodynamic energy of phase k over the layer of height h k . Considering smooth solutions of (2.20) and combining (2.20a) with (2.20b), the equation for ρ k writes:

∂ ρ k ∂t + u k ∂ ρ k ∂ x + ρ k ∂ u k ∂ x + ρ k h k (u k -U I ) ∂ h k ∂ x = (-1) k ρ k h k W I + M k .
Multiplying the last equation by m k Ψ k , we get for E t,k :

∂ E t,k ∂t + ∂ u k E t,k ∂ x + m k ρ k Ψ k ∂ u k ∂ x + ρ 2 k (u k -U I )Ψ k ∂ h k ∂ x = (-1) k ρ 2 k Ψ k W I + (ρ k Ψ k +Ψ k )h k M k . (3.2)
Multiplying (2.20c) by u k and combining with (2.20b) we get for E c,k :

∂ E c,k ∂t + ∂ ∂ x (u k E c,k + u k h k P k ) -h k P k ∂ u k ∂ x -u k P 1 ∂ h k ∂ x + u k ρ 1 g h 1 2 ∂ h k ∂ x = u k h k D k - u 2 k 2 h k M k . (3.3)
For phase 1, the term

u 1 ρ 1 g h 1 2 ∂ h 1 ∂ x in (3.
3) may be used to introduce the gravitational potential energy:

E p,1 = ρ 1 g h 2 1 2 . (3.4)
After calculations, one obtains:

u 1 ρ 1 g h 1 2 ∂ h 1 ∂ x = ∂ E p,1 ∂t + ∂ u 1 E p,1 ∂ x +U I ρ 1 g h 1 2 ∂ h 1 ∂ x -ρ 1 g h 1 2 W I -g h 2 1 2 M k . (3.5)
The energy equation for the whole system is then obtained adding (3.2) with (3.3) for k = 1, 2. Using (3.5) and the following consistency relations:

∂ h 1 ∂ x + ∂ h 2 ∂ x = 0, h 1 D 1 + h 2 D 2 = 0, h 1 M 1 + h 2 M 2 = 0,
it yields:

∂ ∂t E c,1 + E t,1 + E p,1 + E c,2 + E t,2 + ∂ ∂ x u 1 (E c,1 + E p,1 + E t,1 ) + u 2 (E c,2 + E t,2 ) + u 1 h 1 P 1 + u 2 h 2 P 2 + κ u 1 ∂ u 1 ∂ x + κ u 2 ∂ u 2 ∂ x + κ h 1 ∂ h 1 ∂ x = κ w I W I + κ D 1 h 1 D 1 + κ M 1 h 1 M 1 , (3.6)
where:

κ u k = h k ρ 2 k Ψ k -h k P k , κ h 1 = (u 1 -U I )ρ 2 1 Ψ 1 -(u 2 -U I )(ρ 2 2 Ψ 2 -ρ 1 g h 1 2 ) + (u 2 -u 1 )P 1 , κ w I = ρ 2 2 Ψ 2 -ρ 2 1 Ψ 1 + ρ 1 g h 1 2 , κ D 1 = u 1 -u 2 , κ M 1 = (ρ 1 Ψ 1 +Ψ 1 + g h 1 2 ) -(ρ 2 Ψ 2 +Ψ 2 ) + 1 2 (u 2 2 -u 2 1 ).
In order to obtain a conservative equation, the contribution of non-conservative terms is canceled out. Thus, we define Ψ k (ρ k ) such that:

Ψ k (ρ k ) = P k (ρ k ) ρ 2 k . (3.8)
This closure yields κ u k = 0 and one obtains the consistency with the definition of thermodynamic energy when dealing with single-phase flows and the Euler system. Setting κ h 1 = 0 and using (3.8), one obtains:

(u 2 -U I )(P 1 -P 2 + ρ 1 g h 1 2 ) = 0, ∀(h 1 , ρ 1 , ρ 2 ),
which gives the closure for U I :

U I = u 2 . (3.9)
To conclude, one can state the following proposition:

Proposition 3.1 Smooth solutions of system (2.20) comply with the entropy inequality:

∂ E ∂t + ∂ G ∂ x ≤ 0, (3.10)
where the entropy E and the entropy flux G are defined by:

E = E c,1 + E p,1 + E t,1 + E c,2 + E t,2 , (3.11a) G = u 1 (E c,1 + E p,1 + E t,1 ) + u 2 (E c,2 + E t,2 ) + u 1 h 1 P 1 + u 2 h 2 P 2 , (3.11b)
with:

E c,k = 1 2 h k ρ k u 2 k , E t,k = h k ρ k Ψ k (ρ k ), E p,1 = ρ 1 g h 2 1 2 ,
and:

Ψ k (ρ k ) = P k (ρ k ) ρ 2 k ,
as soon as the following closure laws are used for the interfacial variables:

U I = u 2 , (3.12a
)

P I = P 1 -ρ 1 g h 1 2 , (3.12b)
and for the source terms:

W I = λ p (P I -P 2 ) = λ p (P 1 -ρ 1 g h 1 2 -P 2 ), (3.13a)
h k D k = (-1) k λ u (u 1 -u 2 ) + u 1 + u 2 2 h k M k , (3.13b
)

h k M k = (-1) k λ m ( P 1 + ρ 1 g h 1 2 ρ 1 +Ψ 1 ) -( P 2 ρ 2 +Ψ 2 ) , (3.13c)
where λ p , λ u and λ m are positive bounded functions which depend on the state variable

(h 1 , ρ 1 , ρ 2 , u 1 , u 2 ).
Proof (3.10) is obtained from (3.6) where the source terms are chosen to ensure that the inequality is verified. Indeed, the right-hand side of (3.6) writes:

S = h 1 D 1 (u 1 -u 2 ) + h 1 M 1 ( P 1 + ρ 1 g h 1 2 ρ 1 +Ψ 1 ) -( P 2 ρ 2 +Ψ 2 ) + u 2 2 -u 2 1 2 +W I (P 2 -P 1 + ρ 1 g h 1 2
),

so that choosing h k D k = h k D k + ( u 1 +u 2 2
)h k M k , one obtains:

S = h k D k (u 1 -u 2 ) + h 1 M 1 ( P 1 + ρ 1 g h 1 2 ρ 1 +Ψ 1 ) -( P 2 ρ 2 +Ψ 2 ) +W I (P 2 -P 1 + ρ 1 g h 1 2 ),
and the corresponding expressions for W I , h k D k and h k M k are chosen such that S ≤ 0.

Therefore, the closed system writes:

∂W ∂t + ∂ F(W ) ∂ x + B(W ) ∂W ∂ x = C(W ), (3.14)
where:

W = (h 1 , h 1 ρ 1 , h 2 ρ 2 , h 1 ρ 1 u 1 , h 2 ρ 2 u 2 ) T , F(W ) =       0 h 1 ρ 1 u 1 h 2 ρ 2 u 2 h 1 (ρ 1 u 2 1 + P 1 ) h 2 (ρ 2 u 2 2 + P 2 )       , B(W ) ∂W ∂ x =        u 2 ∂ h 1 ∂ x 0 0 -(P 1 -ρ 1 g h 1 2 ) ∂ h 1 ∂ x -(P 1 -ρ 1 g h 1 2 ) ∂ h 2 ∂ x       
, and:

C(W ) =             λ p (P 1 -ρ 1 g h 1 2 -P 2 ) -λ m ( P 1 +ρ 1 g h 1 2 ρ 1 +Ψ 1 ) -( P 2 ρ 2 +Ψ 2 ) λ m ( P 1 +ρ 1 g h 1 2 ρ 1 +Ψ 1 ) -( P 2 ρ 2 +Ψ 2 ) -λ u (u 1 -u 2 ) -( u 1 +u 2 2 )λ m ( P 1 +ρ 1 g h 1 2 ρ 1 +Ψ 1 ) -( P 2 ρ 2 +Ψ 2 ) λ u (u 1 -u 2 ) + ( u 1 +u 2 2 )λ m ( P 1 +ρ 1 g h 1 2 ρ 1 +Ψ 1 ) -( P 2 ρ 2 +Ψ 2 )            
, where λ p , λ u and λ m are positive bounded functions which depend on the state variable

(h 1 , ρ 1 , ρ 2 , u 1 , u 2 ).
4 Comments on the closed system

Consistency with the shallow water equations

When it comes to free-surface flows, the well-known (incompressible) shallow water equations are usually considered in the literature. It is thus interesting to check the consistency of the present model with that classical description. Regarding the closed system (3.14) without mass nor momentum transfers between the layers, averaged mass and momentum conservation equations for the water phase write:

∂ h 1 ρ 1 ∂t + ∂ h 1 ρ 1 u 1 ∂ x = 0, ∂ h 1 ρ 1 u 1 ∂t + ∂ h 1 ρ 1 u 2 1 ∂ x + ∂ h 1 P 1 ∂ x -P I ∂ h 1 ∂ x = 0,
where P 1 = P I + ρ 1 g h 1 2 (averaged hydrostatic constraint). Focusing on the pressure gradient and smooth solutions, it comes:

∂ h 1 P 1 ∂ x = ∂ ρ 1 g h 2 1 2 ∂ x + h 1 ∂ P I ∂ x + P I ∂ h 1 ∂ x ,
and the water layer system reads:

∂ h 1 ρ 1 ∂t + ∂ h 1 ρ 1 u 1 ∂ x = 0, (4.1a) ∂ h 1 ρ 1 u 1 ∂t + ∂ h 1 ρ 1 u 2 1 ∂ x + ∂ ρ 1 g h 2 1 2 ∂ x + h 1 ∂ P I ∂ x = 0. (4.1b)
Thus, if ρ 1 is considered as a constant in (4.1), one obtains formally the classical (incompressible) shallow water model with varying atmospheric pressure P I (x,t), see [START_REF] Agoshkov | Mathematical and numerical modelling of shallow water flow[END_REF] for instance. Furthermore, regarding the surface dynamic equation (3.14a), its interpretation is given in terms of pressure relaxation. Indeed, considering static fluids without mass transfer, one can write:

∂ h 1 ∂t = λ p (P I -P 2 ) = λ p Π (x,t), ∂ h k ρ k ∂t = 0,
where Π (x,t) = P I -

P 2 = P 1 -ρ 1 g h 1 2 -P 2 . The second equation gives h k ∂ ρ k ∂t = (-1) k ρ k ∂ h 1
∂t and Π (x,t) verifies:

∂ Π ∂t = - ∂ h 1 ∂t ( c 2 1 ρ 1 h 1 + c 2 2 ρ 2 h 2 ), = -λ p ( c 2 1 ρ 1 h 1 + c 2 2 ρ 2 h 2 )Π (x,t),
which yields:

Π (x,t) = Π (x, 0)exp(- t 0 λ p ( c 2 1 ρ 1 h 1 + c 2 2 ρ 2 h 2 )dt).
As λ p is a positive bounded function, one obtains the following asymptotic behavior:

P I -→ t→+∞ P 2 . (4.2)
In practice, dealing with water for phase 1 and pipe radii of about 1 m,

c 2 1 ρ 1 h 1 ∼ 10 9
Pa.m -1 and the relaxation is very fast. Consequently, the atmospheric pressure in (4.1b) quickly converges toward the pressure of the air phase, which makes sense regarding the layered configuration. Note that an explicit form for the pressure relaxation term λ p is proposed in [START_REF] Gavrilyuk | The structure of pressure relaxation terms: one-velocity case[END_REF] for bubbly flows and might be extended to our model.

The consistency of the proposed model with the incompressible shallow water model deserves also a deep investigation considering an asymptotic low Mach number development in the liquid, but this study lies beyond the scope of the work exposed herein.

Consistency with pressurized flows

The case where the pipe is full of phase 1 is referred to as pressurized flow. In practice, transitions from stratified to pressurized regime often occur in industrial facilities so that one may wonder if this configuration will be correctly handled by our model. Formally, considering h 1 = H, one obtains:

∂ ρ 1 ∂t + ∂ ρ 1 u 1 ∂ x = 0, ∂ ρ 1 u 1 ∂t + ∂ (ρ 1 u 2 1 + P 1 ) ∂ x = 0,
as soon as the source terms vanish when h 1 = H. This system gives the expected averaged equations to describe a rectangular channel full of phase 1 where compressibility effects provide the correct velocity of acoustic waves. Therefore, the model presented herein for stratified air-water flows degenerates correctly to the pressurized regime without any switching on the pressure law, contrary to the model presented in [START_REF] Bourdarias | A finite volume scheme for a model coupling free surface and pressurised flows in pipes[END_REF].

Comparison with the two-fluid two-pressure models

By construction, the model presented herein may be seen as a two-fluid two-pressure model. Indeed, looking at the system (3.14), one may recognize a two-velocity two-pressure model as developed in [START_REF] Baer | A two phase mixture theory for the deflagration to detonation (DDT) transition in reactive granular materials[END_REF][START_REF] Coquel | Closure laws for a two-fluid two-pressure model[END_REF][START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF][START_REF] Glimm | Two phase flow modelling of a fluid mixing layer[END_REF][START_REF] Kapila | Two-phase modeling of DDT: structure of the velocity-relaxation zone[END_REF] without the energy equations. The first difference relies on the averaging process; in the literature, one may find time averages, space averages or statistical averages, resulting in the same system structure. Another difference is given by the closure laws (3.12). In the overall framework, see [START_REF] Coquel | Closure laws for a two-fluid two-pressure model[END_REF], several closure laws are possible for (U I , P I ) to get both a linearly degenerated field for the 1-wave and an entropy inequality, whereas in our case, there is only one possibility satisfying the additional hydrostatic constraint (2.3). That being said, most of the good properties such as hyperbolicity (see next section) or entropy inequality have been inherited directly from the overall framework. One may also find similarities with the model developed in [START_REF] Gallouët | Hyperbolic relaxation models for granular flows[END_REF] where granular flows have been described in the statistical framework. Particularly, reading h k as the statistical fraction, similarities are found regarding the source terms: W I includes the granular stress contribution, h k D k the classical friction effects and h k M k the disequilibrium of Gibbs enthalpy. In the (O, x, y, z) frame of reference, the hydrostatic constraint writes:

∂ P 1 ∂ z = -ρ 1 g cos θ , (4.3)
and one may easily demonstrate that the closing relations for the interfacial variables are exactly the same as those presented in the previous section, replacing g by g cos θ :

U I = u 2 , P I = P 1 -ρ 1 g h 1 2 cos θ ,
and for the source terms:

W I = λ p (P 1 -ρ 1 g h 1 2 cos θ -P 2 ), h k D k = (-1) k λ u (u 1 -u 2 ) + u 1 + u 2 2 h k M k , h k M k = (-1) k λ m ( P 1 + ρ 1 g( h 1 2 cos θ + b) ρ 1 +Ψ 1 ) -( P 2 ρ 2 +Ψ 2 ) .
The entropy inequality being equivalent to (3.10) where E p,1 now writes:

E p,1 = ρ 1 g h 2 1 2 cos θ + ρ 1 gh 1 b,
with db dx = sinθ . Note that there is no approximation required on θ to get these relations. The only implicit assumption is that phases stay layered with phase 1 below. The validity of this system is thus questionable when θ gets close to π 2 . To conclude, one may easily deal with sloped pipes where θ is constant. The case with variable slopes is not considered here although it can a priori be handled considering the curvilinear formulation of Euler equations as in [START_REF] Bouchut | On new erosion models of Savage-Hutter type for avalanches[END_REF][START_REF] Bouchut | A new model of Saint Venant and Savage-Hutter type for gravity driven shallow water flows[END_REF][START_REF] Bourdarias | A mathematical model for unsteady mixed flows in closed water pipes[END_REF][START_REF] Morales De Luna | A Saint Venant model for gravity driven shallow water flows with variable density and compressibility effects[END_REF].

Mathematical properties

Dealing with the closed system (3.14), the homogeneous problem is studied through the eigenstructure and hyperbolicity analysis, the nature of characteristic fields and the associated Riemann invariants. Furthermore, as non-conservative terms exist, uniqueness of jump conditions is then studied carefully. Finally, positivity is obtained for the physical variables h k and ρ k , k = 1, 2.

Eigenstructure and hyperbolicity

Regarding the convective part of (3.14):

∂ h 1 ∂t + u 2 ∂ h 1 ∂ x = 0, (5.1a) ∂ h k ρ k ∂t + ∂ h k ρ k u k ∂ x = 0, (5.1b) ∂ h k ρ k u k ∂t + ∂ h k (ρ k u 2 k + P k ) ∂ x -(P 1 -ρ 1 g h 1 2 ) ∂ h k ∂ x = 0, (5.1c)
where h 1 + h 2 = H, and recalling that c k = P k (ρ k ) denotes the celerity of acoustic waves, the following proposition holds:

Proposition 5.1 The homogeneous problem (5.1) is hyperbolic under the condition:

|u 1 -u 2 | = c 1 . (5.2)
Its eigenvalues are unconditionally real and given by:

λ 1 = u 2 , λ 2 = u 1 -c 1 , λ 3 = u 1 + c 1 , λ 4 = u 2 -c 2 , λ 5 = u 2 + c 2 .
(5.3)

Defining the vector:

Y = (h 1 , ρ 1 , u 1 , ρ 2 , u 2 ) T , (5.4)
the corresponding right eigenvectors write in column:

R(Y ) =       1 0 0 0 0 η 1 ρ 1 ρ 1 0 0 η 2 -c 1 c 1 0 0 η 3 0 0 ρ 2 ρ 2 0 0 0 -c 2 c 2      
, where:

η 1 = ρ 1 c 2 1 -(u 1 -u 2 ) 2 (u 1 -u 2 ) 2 h 1 - g 2 , η 2 = u 1 -u 2 c 2 1 -(u 1 -u 2 ) 2 (- c 2 1 h 1 + g 2 ), η 3 = P 2 -P 1 + ρ 1 g h 1 2 h 2 c 2 2 .
Proof (5.1) may be written with respect to Y , which gives, for regular solutions:

∂Y ∂t +C(Y ) ∂Y ∂ x = 0, (5.5)
where

C(Y ) =         u 2 0 0 0 0 ρ 1 h 1 (u 1 -u 2 ) u 1 ρ 1 0 0 g 2 c 2 1 ρ 1 u 1 0 0 0 0 0 u 2 ρ 2 P 1 -ρ 1 g h 1 2 -P 2 ρ 2 h 2 0 0 c 2 2 ρ 2 u 2         .
The eigenvalues and the eigenvectors are then easily computed using the block structure of C(Y ). As c k = 0, the eigenvectors are linearly independent and span R 5 as soon as λ 1 = u 2 is different from the other eigenvalues, which may be rewritten under the condition (5.2).

Note that the hyperbolicity is not strictly verified everywhere since some eigenvalues may coincide if (5.2) is violated, leading to the so-called resonant behavior. However, in the context of air-water flows, c 1 ≈ 1500 m.s -1 , and the resonant situation is clearly out of the scope of interest and the model is not devoted to this unrealistic flow regime. Finally, looking at C(Y ) and R(Y ), phases 1 and 2 are only coupled by the first column which means that they evolve independently on each side of the 1-wave λ 1 = u 2 . Moreover, the system may be symmetrized using the Y variable, see [START_REF] Coquel | Two properties of two-velocity two-pressure models for two-phase flows[END_REF] for some counterpart.

Study of the Riemann problem

In this section, we focus on the Riemann problem associated with (5.1) in order to verify that the parametrization of waves λ k , k = 1, .., 5, is well defined. Using the conservative variable:

W = (h 1 , h 1 ρ 1 , h 2 ρ 2 , h 1 ρ 1 u 1 , h 2 ρ 2 u 2 ) T ,
the Riemann problem writes:

W (t = 0, x) = W L , if x < 0, W R , if x > 0,
where W L and W R are some constant states. After studying the nature of characteristic fields, one turns to Riemann invariants and jump conditions which need careful consideration regarding non-conservative terms.

Nature of characteristic fields and Riemann invariants

Considering the eigenstructure of (5.1) detailed in Proposition 5.1, the following proposition can be stated:

Proposition 5.2
The field associated with the 1-wave λ 1 is linearly degenerate while the fields associated with the waves λ k , k = 2, .., 5, are genuinely nonlinear. Moreover, denoting I k (W ) the vector of k-Riemann invariants associated with the k-wave, one obtains:

I 1 (W ) = u 2 , m 1 (u 1 -u 2 ), m 1 u 1 (u 1 -u 2 ) + h 1 P 1 + h 2 P 2 , P 1 + ρ 1 g h 1 2 ρ 1 +Ψ 1 + (u 1 -u 2 ) 2 2 T , I 2 (W ) = h 1 , ρ 2 , u 2 , u 1 + c 1 (ρ 1 ) ρ 1 dρ 1 T , I 3 (W ) = h 1 , ρ 2 , u 2 , u 1 - c 1 (ρ 1 ) ρ 1 dρ 1 T , I 4 (W ) = h 1 , ρ 1 , u 1 , u 2 + c 2 (ρ 2 ) ρ 2 dρ 2 T , I 5 (W ) = h 1 , ρ 1 , u 1 , u 2 - c 2 (ρ 2 ) ρ 2 dρ 2 T .
Proof One may readily see that:

∇ Y λ 1 (Y ).r 1 (Y ) = 0,
where r k (Y ) denotes the k th column of R(Y ). Thus, the 1-wave field is linearly degenerate. Moreover, a classical result derived from the Euler system with perfect gas law or stiffened gas law gives:

∇ Y λ k (Y ).r k (Y ) = 0, k = 2, ..., 5.
Thus, the associated fields are genuinely nonlinear.

Regarding the matrix R(Y ) in Proposition 5.1, the definition of the k-Riemann invariants, k = 2, .., 5, is also a classical result derived from the Euler system. Concerning the linearly degenerate field, λ 1 = u 2 is a 1-Riemann invariant by definition. The three remaining 1-Riemann invariants may be obtained using three conservation laws on which Rankine-Hugoniot jump conditions are applied with σ = λ 1 as speed of discontinuity. Mass conservation and momentum conservation summed on both phases provide the first two, while the last one is the entropy equality without source terms (see Proposition 3.1). Then, denoting

I 1 k the k th component of I 1 , checking that ∇ Y I 1 k (Y ).r 1 (Y ) = 0 is straightforward though tedious.
This result emphasizes the fact that the phases evolve independently on each side of the 1-wave. In those regions, h 1 is constant and one may notice that the system (5.1) reduces locally to two conservative Euler systems.

Jump conditions

System (5.1) contains two non-conservative products, u 2 ∂ x h 1 and (P 1ρ 1 g h 1 2 )∂ x h 1 . Regarding discontinuous solutions, one has to make sure that those products are well defined across genuinely nonlinear and linearly degenerate fields.

As h 1 is constant through the genuinely nonlinear fields, the system (5.1) can be considered locally conservative. Thus, one may state the following proposition: Proposition 5.3 For all genuine nonlinear fields corresponding to the k-waves, k = 2, ..., 5, the Rankine-Hugoniot jump conditions across a single discontinuity of speed σ write:

[h k ] = 0, [m k (u k -σ )] = 0, [m k u k (u k -σ ) + h k P k ] = 0,
where brackets [.] denote the difference between the states on both sides of the discontinuity.

Furthermore, as the field associated with the jump of h 1 is linearly degenerate, the non-conservative products u 2 ∂ x h 1 and

(P 1 -ρ 1 g h 1
2 )∂ x h 1 are well defined. Indeed, one may use the available 1-Riemann invariants detailed in Proposition 5.2 to write explicitly the 1-wave parametrization.

Finally, no ambiguity holds in the definition of jump relations and non-conservative products. In addition, one can build analytical solutions for system (5.1) which may be used to validate numerical schemes, see [19] for details. Nonetheless, the great complexity of (5.1) seems to prohibit the exact resolution of Riemann problems.

Positivity

Considering variables with physical meaning such as h k or ρ k , ensuring their positivity is a major requirement. Focusing on smooth solutions, the result is classical regarding the system (3.14) if one assumes that λ p and λ m may be written under the form λ p = m 1 m 2 λp and λ m = m 1 m 2 λm , where λp and λm are positive bounded functions depending on the state variable (see Proposition B.1 in Appendix B). In addition, one can demonstrate that the positivity requirements hold for discontinuous solutions of the Riemann problem associated with the homogeneous system (5.1). It relies on an explicit writing of elementary waves parametrizations.

Extension to circular pipes with variable cross section

As mentioned earlier, the 1D averaging operator proposed in (2.9) implicitly deals with the case of homogeneous channels. As most of pipes are actually circular in industrial facilities, the case of horizontal circular pipes is now considered. The analysis is also extended to variable cross section in space and in time corresponding to geometric constraint and pipe elasticity. To this end, a 2D averaging operator is applied on 3D Euler system following the same approach as in Sect. 2. Note that the variable cross section case is developed in [START_REF] Castro | Numerical simulation of two-layer shallow water flows through channels with irregular geometry[END_REF] considering open channels and two incompressible layers.

Local governing equations and geometric description

The framework is the set of 3D isentropic Euler equations with the same hypothesis as in Sect. 2, that is the classical hydrostatic assumption for water. Denoting u k = (u k , v k , w k ) the velocity vector of phase k and v k = (0, v k , w k ), the system along the longitudinal direction writes:

∂ ρ k ∂t + ∂ ρ k u k ∂ x + div y,z (ρ k v k ) = M k , (6.1a) ∂ ρ k u k ∂t + ∂ ρ k u 2 k ∂ x + div y,z (ρ k u k v k ) + ∂ P k (ρ k ) ∂ x = D k , (6.1b)
with the hydrostatic constraint for phase 1:

∂ P 1 (ρ 1 ) ∂ z = -ρ 1 g. (6.2)
It is assumed that the flow is homogeneous in the spanwise direction y, so that (ρ k ,u k ) are functions of (x, z,t) and the interface level η 1 (x,t) is constant along the y direction, see Fig. 6.1 for a sketch of the problem. The pipe is assumed to have a symmetry axis C which coincides with the x axis. In the cross section planes (Oyz), e z and e r denote, respectively, the unit normal vector to the interface and the outward unit normal vector to the wall, σ (x, z,t) is the width of the cross section and R(x,t) the radius of the pipe. A k (x,t) refers to the area filled by phase k such that A 1 (x,t) + A 2 (x,t) = S(x,t), the area of the cross section. The normal vector to the interface is given by n I , which differs from e z except when the height profile is horizontal.

σ (x, z,t) h 1 h 2 η 1 e z e y θ 1 • R(x,t)
In that framework, the position of the interface is determined from the symmetry axis C and given by η 1 (x,t) which verifies the following kinetic boundary condition:

∂ η 1 ∂t +U I ∂ η 1 ∂ x = W I , (6.3)
where u I = (U I ,V I ,W I ) T is the interfacial velocity. On the interface, the following boundary conditions still apply:

(u 1 (x, z = η - 1 ,t) -u 2 (x, z = η + 1 ,t)).n I = 0, (6.4a)
P 1 (x, z = η - 1 ,t) = P 2 (x, z = η + 1 ,t) = P I (x,t), (6.4b)
and we will assume that the interfacial velocity still satisfies:

u I (x,t) = β u 1 (x, z = η - 1 ,t) + (1 -β )u 2 (x, z = η + 1 ,t), β ∈ [0, 1]. (6.5)
On the walls, we will consider a contact boundary condition between phase k and the wall, which may be written as:

dOM b dt wall = u k,wall , (6.6)
where M b (x, y, z) is a material point belonging to the contour of the area filled by phase k in the cross section and defined as:

OM b = xe x + Re r
, on the wall, (6.7a) = xe x + ye y + η 1 (x,t)e z , on the interface. (

Note that the impermeability condition (2.5) used in Sect. 2 is consistent with (6.6).

Averaging process

Let us define Ω k (x,t) the integration domain for phase k,

Ω k (x,t) = {(y, z) ∈ R 2 ; - σ (x, z,t) 2 ≤ y ≤ σ (x, z,t) 2 , -R(x,t) + z k ≤ z ≤ -R(x,t) + z k + h k (x,t)},
with z 1 = 0, z 2 = h 1 (x,t), so that the averaging operator writes:

f k (x,t) = 1 A k (x,t) Ω k (x,t) f k (x, z,t)dΩ (6.8)
where A k (x,t) is the area filled by phase k in the cross section at abscissa x and time t:

A k (x,t) = Ω k (x,t)
dΩ . (6.9)

Applying the operator (6.8) to (6.1), one will use the Reynolds transport theorem in time and space given below:

Ω k (x,t) ∂ f k ∂t dΩ = ∂ A k f k ∂t - ∂ Ω k f k ∂ OM b ∂t .n k dl, (6.10a
)

Ω k (x,t) ∂ f k ∂ x dΩ = ∂ A k f k ∂ x - ∂ Ω k f k ∂ OM b ∂ x .n k dl, (6.10b)
and the divergence theorem:

Ω k (x,t) div y,z ( f k v k )dΩ = ∂ Ω k f k v k .n k dl, (6.11)
where OM b is defined in (6.7) and n k denotes the outward unit normal vector to ∂ Ω k .

∂ Ω k includes two different parts, a first part Γ I on the interface and a second part Γ w,k on the wall such that:

∂ Ω k = Γ I + Γ w,k . (6.12)
Thus, regarding the contour integral in (6.10b) and using the definition (6.7) of OM b , one writes for the wall part:

Γ w,k f k ∂ OM b ∂ x .n k dl = Γ w,k f k ∂ R ∂ x Rdθ = f k θ k 2π ∂ S ∂ x , (6.13)
where:

f k = 1 θ k Γ w,k f k dθ , θ k = Γ w,k dθ , (6.14)
and S(x,t) = πR(x,t) 2 is the area of the cross section. Note that f k corresponds to the linear average of f k along the wet (k = 1) or dry circular contour (k = 2).

In addition, note that the kinematic boundary condition (6.3) may be written with A 1 . Indeed, using (6.10b) and (6.13) with f k = 1, one obtains:

∂ A k ∂ x = σ I (-1) k+1 ∂ η 1 ∂ x + θ k 2π ∂ S ∂ x , (6.15)
where σ I is the width of the cross section at the interface. The same calculations with (6.10a) yield:

∂ A k ∂t = σ I (-1) k+1 ∂ η 1 ∂t + θ k 2π ∂ S ∂t ,
and it follows:

∂ A 1 ∂t +U I ∂ A 1 ∂ x = σ I W I + θ 1 2π ( ∂ S ∂t +U I ∂ S ∂ x ). (6.16)

Mass conservation

The 2D averaging of (6.1a) using (6.10) and (6.11) yields:

∂ A k ρ k ∂t + ∂ A k ρ k u k ∂ x - ∂ Ω k ρ k ∂ OM b ∂t + u k ∂ OM b ∂ x -v k .n k dl = A k M k .
On the interface, n k = (-1) k+1 e z and the definition (6.7) of OM b gives for the Γ I part:

Γ I ρ k (-1) k+1 ( ∂ η 1 ∂t + u k ∂ η 1 ∂ x -w k )dl = Γ I B k dl = 0,
the nullity of B k being provided by the kinetic boundary condition (6.3) combined with the continuity of the normal velocity on the interface (6.4a), see (2.12) for details. On the wall, n k = e r and the definition (6.7) of OM b gives for the Γ w,k part:

Γ w,k ρ k ∂ OM b ∂t + u k ∂ OM b ∂ x -v k .e r dl = Γ w,k ρ k dOM b dt wall -u k,wall .e r dl = 0, noticing that v k .e r = u k .e r , v k ∂ OM b ∂ y + w k ∂ OM b ∂ z
.e r = 0 on the wall, and using the contact boundary condition (6.6).

Finally, the contour integral on ∂ Ω k is zero and using the notation ρ k u k = ρ k u k , the averaged mass conservation equations write:

∂ A k ρ k ∂t + ∂ A k ρ k u k ∂ x = A k M k . (6.17)
Note that this equation is correct for all geometries of pipes, even for non-circular ones, as the nullity of the contour integral does not require a circular section.

Momentum conservation

The 2D averaging of (6.1b) yields:

∂ A k ρ k u k ∂t + ∂ A k (ρ k u 2 k + P k ) ∂ x - ∂ Ω k ρ k u k ∂ OM b ∂t + u k ∂ OM b ∂ x -v k .n k dl - ∂ Ω k P k ∂ OM b ∂ x .n k dl = A k D k .
The first contour integral is zero using the same decomposition as for the mass conservation. The contour integral on P k provides:

∂ Ω k P k ∂ OM b ∂ x .n k dl = Γ I (-1) k+1 P k ∂ OM b ∂ x .e z dl + Γ w,k P k ∂ OM b ∂ x
.e r dl,

= P I σ I (-1) k+1 ∂ η 1 ∂ x + P k θ k 2π ∂ S ∂ x , = P I ∂ A k ∂ x + ( P k -P I ) θ k 2π ∂ S ∂ x ,
using (6.15). Recall that P k corresponds to the linear average of P k along the wet (k = 1) or dry circular contour (k = 2).

Neglecting the turbulence of the flow at first order, the closure law ρ k u 2 k = ρ k u 2 k is chosen and one obtains:

∂ A k ρ k u k ∂t + ∂ A k (ρ k u 2 k + P k ) ∂ x -P I ∂ A k ∂ x = A k D k + ( P k -P I ) θ k 2π ∂ S ∂ x . (6.18)

Hydrostatic constraint

As in Sect. 2, one has to account for gravitational effects for phase 1 coupling (6.18) with the hydrostatic constraint (6.2). An original approach is proposed below. An integration of (6.2) between z and η 1 writes:

P 1 (ρ 1 (x, z,t)) = P I (x,t) + η 1 z ρ 1 (x, z,t)gds,
which gives:

A 1 P 1 = Ω 1 P I dΩ + Ω 1 η 1 z ρ 1 gds dΩ , (6.19a)
= A 1 P I + η 1 -R η 1 z ρ 1 gds σ dz, (6.19b)
as ρ 1 (x, z,t) and P I (x,t) do not depend on y. Let us define A (x, z,t) the primitive function of σ (x, z,t) along z which cancels out in -R:

A (x, z,t) = z -R σ (x, z,t)dz.
An integration by parts provides:

η 1 -R η 1 z ρ 1 gds σ dz = η 1 z ρ 1 gds A η 1 -R + η 1 -R ρ 1 gA dz, = η 1 -R ρ 1 gA dz, = gρ 1 A σ .
As in the depth-averaged case, it is proposed to neglect the correlation between ρ 1 and A σ (see Appendix A for details), such that the closure law for ρ 1 A σ writes:

ρ 1 A σ = ρ 1 A σ , with: A σ = 1 A 1 η 1 -R A dz, = 1 A 1 [(z -η 1 )A ] η 1 -R + η 1 -R (η 1 -z)σ dz , = 1 A 1 η 1 -R (η 1 -z)σ dz, = η 1 -z.
Thus, (6.19b) yields:

P I = P 1 -ρ 1 g 1 , (6.20)
where:

1 (x,t) = η 1 (x,t)z (6.21)

represents the distance between the interface and the center of mass of the wet section. (6.20) provides the closure law for P I regarding P 1 as a function of ρ 1 . Moreover, the latter may be read as a sectionaveraged hydrostatic pressure law for P 1 used with the averaged density ρ 1 . Therefore, the averaged momentum equations under hydrostatic constraint write:

∂ A k ρ k u k ∂t + ∂ A k (ρ k u 2 k + P k ) ∂ x -(P 1 -ρ 1 g 1 ) ∂ A k ∂ x = A k D k + P k -(P 1 -ρ 1 g 1 ) θ k 2π ∂ S ∂ x . (6.22)

Resulting averaged system

Adding the interfacial kinetic boundary condition (6.16), one obtains a five-equation system corresponding to the five unknowns

(A 1 , ρ 1 , ρ 2 , u 1 , u 2 ): ∂ A 1 ∂t +U I ∂ A 1 ∂ x = σ I W I + θ 1 2π ( ∂ S ∂t +U I ∂ S ∂ x ), (6.23a)
∂ A k ρ k ∂t + ∂ A k ρ k u k ∂ x = A k M k , (6.23b)
∂ A k ρ k u k ∂t + ∂ A k (ρ k u 2 k + P k ) ∂ x -(P 1 -ρ 1 g 1 ) ∂ A k ∂ x = A k D k + P k -(P 1 -ρ 1 g 1 ) θ k 2π ∂ S ∂ x , (6.23c)
where k = 1, 2 and A 1 + A 2 = S, θ 1 + θ 2 = 2π. Dealing with the source terms, conservation of mass and momentum of the mixture impose:

A 1 D 1 + A 2 D 2 = 0, (6.24a
)

A 1 M 1 + A 2 M 2 = 0. (6.24b)
As in subsect. 2.3, the closure P k (ρ k ) = P k (ρ k ) is chosen and the celerity of acoustic waves in the averaged framework is defined by:

c k = P k (ρ k ). (6.25)
Note that the structure of (6.23) is the same as (2.20) obtained with 1D averaging, replacing A k by h k . Moreover, considering rectangular pipes of constant width L, one has A k = Lh k , 1 = h 1 2 , and both systems are equivalent. Dealing with circular pipes, 1 is given by the formula (6.21) which may be detailed as:

1 = R 3 A 1 2 3 sin 3 θ 1 2 - 1 2 cos θ 1 2 (θ 1 -sin θ 1 ) , (6.26)
where θ 1 = 2arccos(1 -h 1 R ). S(x,t) is a given function which accounts for cross section variations and as in Sect. 2, P I is closed by the averaged hydrostatic constraint (6.20). The interfacial velocity U I and the source terms, including P k , still need closure laws. To this end, an entropy characterization is detailed in the following section.

Closure laws and system properties

From now on, the operator notations are omitted. As (6.23) and (2.20) are structurally identical, the calculations will not be detailed in this section. Proposition 6.1 System (6.23) admits the entropy inequality:

∂ E ∂t + ∂ G ∂ x + S t ∂ S ∂t ≤ 0, (6.27)
where the entropy E , the entropy flux G and S t are defined by:

E = E c,1 + E p,1 + E t,1 + E c,2 + E t,2 , (6.28a) G = u 1 (E c,1 + E p,1 + E t,1 ) + u 2 (E c,2 + E t,2 ) + u 1 A 1 P 1 + u 2 A 2 P 2 , (6.28b
)

S t = θ 1 2π (P 1 + ρ 1 g( η 1 -z -1 )) + θ 2 2π P 2 , (6.28c)
with:

E c,k = 1 2 A k ρ k u 2 k , E t,k = A k ρ k Ψ k (ρ k ), E p,1 = A 1 ρ 1 g(η 1 -1 ),
and:

Ψ k (ρ k ) = P k (ρ k ) ρ 2 k ,
as soon as the following closure laws are used for the interfacial variables:

U I = u 2 , (6.29a
)

P I = P 1 -ρ 1 g 1 , (6.29b)
and for the source terms:

P 1 = P 1 + ρ 1 g( η 1 -z -1 ), (6.30a
)

P 2 = P 2 , (6.30b)
σ I W I = λ p (P I -P 2 ) = λ p (P 1ρ 1 g 1 -P 2 ), (6.30c)

A k D k = (-1) k λ u (u 1 -u 2 ) + u 1 + u 2 2 A k M k , (6.30d
)

A k M k = (-1) k λ m ( P 1 + ρ 1 g(η 1 -1 ) ρ 1 +Ψ 1 ) -( P 2 ρ 2 +Ψ 2 ) , (6.30e)
where λ p , λ u and λ m are positive bounded functions which may depend on the state variable (A 1 , ρ 1 , ρ 2 , u 1 , u 2 ).

Proof Developing similar calculations as in Sect. 3 and using the closure laws (6.30c), (6.30d), (6.30e), one ends up with the following entropy inequality:

∂ E ∂t + ∂ G ∂ x + S t ∂ S ∂t + S x ∂ S ∂ x ≤ 0,
where E , G , S t are defined in (6.28) and:

S x = u 1 θ 1 2π (P 1 + ρ 1 g( η 1 -z -1 ) -P 1 ) + u 2 θ 2 2π (P 2 -P 2 ).
Thus, the suggested closure laws (6.30a) and (6.30b) for P 1 and P 2 cancel out the S x contribution and yield a decreasing entropy when dealing with section area constant in time. In addition, the latter are also consistent with the single-phase steady state at rest, see Remark 6.2.

Remark 6.1 In order to clarify (6.28c) and (6.30a), η 1z may be detailed as follows:

η 1 -z = 1 θ 1 Γ w,1 (η 1 -z)dθ = η 1 - 2R θ 1 sin( θ 1 2
).

Thus, when η 1 = R and A 1 = S, Rz = 1 = R and P 1 = P 1 . Moreover, note that A k = S yields S t

∂ S ∂t = P k ∂ S
∂t in (6.27), which is consistent with the single-phase case. Remark 6.2 The suggested closure laws (6.30a) and (6.30b) are consistent with the single-phase steady state at rest. Indeed, regarding the averaged momentum equation (6.23c) with u k = 0, one obtains:

∂ A k P k ∂ x -(P 1 -ρ 1 g 1 ) ∂ A k ∂ x = A k D k + P k -(P 1 -ρ 1 g 1 ) θ k 2π ∂ S ∂ x .
Considering uniform pressure and single-phase flow, that is ∂ P k ∂ x = 0, A k = S, θ k = 2π and D k = 0, the above equation reads:

(P k -P k ) A k =S ∂ S ∂ x = 0,
and the single-phase steady state at rest yields:

P k A k =S = P k A k =S ,
which is verified by (6.30b) and (6.30a), see Remark 6.1.

The closed system describing compressible two-layer flows in circular pipes with variable cross section writes:

∂W ∂t + ∂ F(W ) ∂ x + B(W ) ∂W ∂ x = C(W ) + D(W ), (6.31)
where:

W = (A 1 , A 1 ρ 1 , A 2 ρ 2 , A 1 ρ 1 u 1 , A 2 ρ 2 u 2 ) T , F(W ) =       0 A 1 ρ 1 u 1 A 2 ρ 2 u 2 A 1 (ρ 1 u 2 1 + P 1 ) A 2 (ρ 2 u 2 2 + P 2 )       , B(W ) ∂W ∂ x =        u 2 ∂ A 1 ∂ x 0 0 -(P 1 -ρ 1 g 1 ) ∂ A 1 ∂ x -(P 1 -ρ 1 g 1 ) ∂ A 2 ∂ x        , C(W ) =           λ p (P 1 -ρ 1 g 1 -P 2 ) -λ m ( P 1 +ρ 1 g(η 1 -1 ) ρ 1 +Ψ 1 ) -( P 2 ρ 2 +Ψ 2 ) λ m ( P 1 +ρ 1 g(η 1 -1 ) ρ 1 +Ψ 1 ) -( P 2 ρ 2 +Ψ 2 ) -λ u (u 1 -u 2 ) -( u 1 +u 2 2 )λ m ( P 1 +ρ 1 g(η 1 -1 ) ρ 1 +Ψ 1 ) -( P 2 ρ 2 +Ψ 2 ) λ u (u 1 -u 2 ) + ( u 1 +u 2 2 )λ m ( P 1 +ρ 1 g(η 1 -1 ) ρ 1 +Ψ 1 ) -( P 2 ρ 2 +Ψ 2 )          
, and:

D(W ) =       θ 1 2π ( ∂ S ∂t + u 2 ∂ S ∂ x ) 0 0 ρ 1 g η 1 -z θ 1 2π ∂ S ∂ x (P 2 -(P 1 -ρ 1 g 1 )) θ 2 2π ∂ S ∂ x      
, where λ p , λ u and λ m are positive bounded functions which depend on the state variable

(A 1 , ρ 1 , ρ 2 , u 1 , u 2 ).
The hyperbolicity of the homogeneous problem associated with (6.31) is readily obtained following the same approach as in Sect. 5, the eigenvalues being unchanged. Proposition 6.2 The homogeneous problem associated with (6.31) is hyperbolic under the condition:

|u 1 -u 2 | = c 1 . (6.32)
Its eigenvalues are unconditionally real, given by:

λ 1 = u 2 , λ 2 = u 1 -c 1 , λ 3 = u 1 + c 1 , λ 4 = u 2 -c 2 , λ 5 = u 2 + c 2 . (6.33)
Dealing with the Riemann problem, the same analysis as in Sect. 5 leads to same results for (6.31). Thus, the nature of characteristic fields associated with the k-waves, k = 1, .., 5, and Riemann invariants can be detailed, as well as positivity properties and definition of jump conditions. Remark 6.3 The comments of Sect. 4 apply to the system (6.31). Particularly, taking A 1 (x,t) = S(x,t) in (6.31), one obtains the following system:

∂ Sρ 1 ∂t + ∂ Sρ 1 u 1 ∂ x = 0, (6.34) ∂ Sρ 1 u 1 ∂t + ∂ S(ρ 1 u 2 1 + P 1 ) ∂ x = P 1 ∂ S ∂ x , (6.35)
as soon as the source terms vanish when A 1 = S. This resulting system correctly models a pressurized flow in a circular pipe with variable cross section where the pipe elasticity, taken into account here with the function S(x,t), has a great influence on the speed of acoustic waves.

Conclusion

A new model is proposed herein to deal with stratified gas-liquid or vapor-liquid flows in pipes with variable cross section.

It is a compressible two-layer model which results from an averaging process of the isentropic Euler set of equations with hydrostatic constraint on the liquid phase. The main difference with the two-layer models issued from the classical literature is that both phases are assumed compressible. Consequently, it more or less enters in the class of two-fluid two-pressure models and significant mathematical properties are obtained. The latter include hyperbolicity, entropy inequality and positivity. Regarding closure laws, the isentropic framework with hydrostatic constraint implies that the interfacial pressure is defined to satisfy this constraint although the interfacial velocity and source terms are provided by the entropy inequality. It is then observed that this closed system may correctly deal with pressurized flows without any switching on the pressure law. Indeed, this feature results from the compressible description of the stratified regime and opens the door, at least formally, to the modeling of transitions from stratified to pressurized flows and entrapped air pockets. In practice, note that industrial flows include low Mach velocities for both phases verifying |u k | c k . Thus, with the aim of performing numerical simulations with the model presented herein, the current work involves the development of a numerical scheme regarding the asymptotic low Mach number behavior of the system. This approach will be validated building analytical solutions thanks to the detailed Riemann invariants and jump conditions of the homogeneous part. The overall mathematical/numerical model will be then assessed using experimental data. Secondly, the robustness of the scheme will be studied regarding pressurized flows and vanishing phases, see [19] for details.

Regarding realistic configurations where H ∼ 1 m, one has ε 1 so that asymptotic expansions of (A.4) and (A.5) write:

ρ 1 z = r 1 h 1 2 1 - 2 3 ε + O(ε 2 ) , (A.6) ρ 1 z = r 1 h 1 2 1 - 1 2 ε + O(ε 2 ) . (A.7)
Therefore, one ends up with the following estimate:

ρ 1 z -ρ 1 z ρ 1 z = ε 6 + O(ε 2 ), (A.8)
which justifies the closure ρ 1 z = ρ 1 z. Note that this development is extended to the section-averaged case presented in Sect.

-Φ -∂ ∂t (Φ + -Φ -) -aΦ -∂ ∂ x (Φ + -Φ -) -Φ -(Φ + -Φ -) ∂ b ∂ x = -Φ -(Φ + -Φ -)m.
Defining the norm . = ( L 0 | . | 2 dx) 1/2 , one may obtain by integration over [0, L]:

∂ ∂t ( Φ -2 ) + L 0 a ∂ ∂ x (Φ -) 2 dx + 2 L 0 (Φ -) 2 ∂ b ∂ x dx = 2 L 0 (Φ -) 2 mdx.
Integrating by parts the second term of the left-hand side gives:

∂ ∂t ( Φ -2 ) + [a(Φ -) 2 ] L 0 = L 0 (Φ -) 2 (2m -2 ∂ b ∂ x + ∂ a ∂ x )dx.
The positive inlet boundary conditions give Φ -(x = 0,t) = Φ -(x = L,t) = 0, and thus, one can write:

∂ ∂t ( Φ -2 ) ≤ L 0 (Φ -) 2 |2m -2 ∂ b ∂ x + ∂ a ∂ x |dx ≤ Φ -2 sup x∈[0,L] |2m -2 ∂ b ∂ x + ∂ a ∂ x |.
Since the initial data on Φ is positive, the Gronwall's lemma gives for any time t in [0, T ]:

Φ -(t) = 0.
Therefore, Φ -is null and Φ remains positive on the whole domain

Fig. 2 . 1

 21 Fig. 2.1 Geometric description for horizontal channels

4. 4 Fig. 4 . 1

 441 Fig. 4.1 Geometric description for sloped pipes

Fig. 6 . 1

 61 Fig. 6.1 Geometric description for circular pipes with variable cross section

6 .

 6 Appendix B Positivity for heights and densities Proposition B.1 Let L and T be two positive and real constants. Assume that u k , ∂ x u k and the first two right-hand sides of (3.14) belong to L ∞ ([0, L]×[0, T ]) for k = 1, 2. Then, the latter equations associated with admissible inlet boundary conditions lead to:h k (x,t) ∈ [0, H], ∀(x,t) ∈ [0, L] × [0, T ], (B.1) ρ k (x,t) ≥ 0, ∀(x,t) ∈ [0, L] × [0, T ], (B.2) when restricting ourselves to regular solutions. Proof We place ourselves in the general case and we consider a function Φ from [0, L] × [0, T ] to R which verifies an equation of the form: b, m are smooth functions from [0, L] × [0, T ] to R and a, ∂ x a, b, ∂ x b, m belong to L ∞ ([0, L] × [0, T ]). Assume that Φ verifies positive inlet boundary conditions, that is Φ(x = 0,t) and Φ(x = L,t) positive for all t in [0,T]. Let introduce the decomposition Φ = Φ + -Φ -, with Φ + ≥ 0, Φ -≥ 0 and Φ + Φ -= 0. Multiplying (B.3) by -Φ - yields:

  [0, L] × [0, T ]. Finally, consider Φ = h k and Φ = Hh k to get (B.1), while Φ = m k provides (B.2).

Acknowledgements The authors would like to thank Benoît de Laage de Meux, Christian Bourdarias and Stéphane Gerbi for their helpful remarks. This work has been partially funded by ANRT and EDF through an EDF-CIFRE contract 749/2014.

Appendix A Error estimate for the closure ρ 1 z = ρ 1 z One considers the following linear pressure law for phase 1:

and yields:

where r 1 (x,t) = ρ 1 (x, 0,t). Using (A.3) and denoting ε

, one obtains: