A model for unsteady mixed flows in non uniform closed water pipes and a well-balanced finite volume scheme

# A model for unsteady mixed flows in non uniform closed water pipes and a well-balanced finite volume scheme

Christian Bourdarias, Mehmet Ersoy and Stéphane Gerbi
Laboratoire de Mathématiques, Université de Savoie
73376 Le Bourget du Lac, France
email: Christian.Bourdarias@univ-savoie.fr,email: Mehmet.Ersoy@univ-savoie.fr,email: Stephane.Gerbi@univ-savoie.fr
###### Abstract

We present the derivation of a new unidirectional model for unsteady mixed flows in non uniform closed water pipes. We introduce a local reference frame to take into account the local perturbation caused by the changes of section and slope. Then an asymptotic analysis is performed to obtain a model for free surface flows and another one for pressurized flows. By coupling these models through the transition points by the use of a common set of variables and a suitable pressure law, we obtain a simple formulation called PFS-model close to the shallow water equations with source terms. It takes into account the changes of section and the slope variation in a continuous way through transition points. Transition point between the two types of flows is treated as a free boundary associated to a discontinuity of the gradient of pressure. The numerical simulation is performed by making use of a Roe-like finite volume scheme that we adapted to take into account geometrical source terms in the convection matrix. Finally some numerical tests are presented.

Keywords : Shallow water, mixed flows, free surface flows, pressurized flows, curvilinear transformation, asymptotic analysis, VFRoe scheme, well-balanced finite volume scheme, hyperbolic system with source terms.

## 1 Introduction

The presented work takes place in a more general framework: the modelling of unsteady mixed flows in any kind of closed pipe taking into account the cavitation problem and air entrapment. We are interested in flows occurring in closed pipes with non uniform sections, where some parts of the flow can be free surface (it means that only a part of the pipe is filled) and other parts are pressurized (it means that the pipe is full). The transition phenomenon between the two types of flows occurs in many situations such as storm sewers, waste or supply pipes in hydroelectric installations. It can be induced by sudden changes in the boundary conditions as failure pumping. During this process, the pressure can reach severe values and may cause damages. The simulation of such a phenomenon is thus a major challenge and a great amount of works was devoted to it these last years (see [13],[14],[25],[29], for instance).

The classical shallow water equations are commonly used to describe free surface flows in open channels. They are also used in the study of mixed flows using the Preissman slot artefact (see for example [13, 29]). However, this technic does not take into account the depressurisation phenomenon which occurs during a water hammer except in recent works [21, 20, 22]. On the other hand the Allievi equations, commonly used to describe pressurized flows, are written in a non-conservative form which is not well adapted to a natural coupling with the shallow water equations.

A model for the unsteady mixed water flows in closed pipes and a finite volume discretisation have been previously studied by two of the authors [7] and a kinetic formulation has been proposed in [9]. We propose here the PFS-model which tends to extend naturally the work in [7] in the case of a closed pipe with non uniform section. For the sake of simplicity, we do not deal with the deformation of the domain induced by the change of pressure. We will consider only an infinitely rigid pipe.

The paper is organized as follows. Section 2 is devoted to the derivation of the free surface model from the D incompressible Euler equations which are written in a suitable local reference frame (following [3, 4]) in order to take into account the local effects produced by the changes of section and the slope variation. The construction of the free surface model is done by a formal asymptotic analysis. Seeking for an approximation at first order gives the model called FS-model. In Section 3, we adapt the derivation of the FS-model to derive the pressurized model, called P-model, from the D compressible Euler equations. Writing the source terms of these two models, P and FS-model, into a unified form and using the same couple of conservative unknowns as in [8], we propose in Section 4 a model for mixed flows, that we call PFS-model . We state some mathematical properties of this model. Section 5 is devoted to the extension of the VFRoe scheme described in [12, 16, 7] that was used for the case of uniform pipes. In Section 6, we show how to construct a convection matrix in order to get an exactly well-balanced scheme. Several numerical tests are presented in Section 7.

#### Notations concerning geometrical variables

• : cartesian reference frame

• : parametrization in the reference frame of the plane curve which corresponds to the main flow axis

• : Serret-Frenet reference frame attached to with T the tangent vector, N the normal vector and B the binormal vector

• : local variable in the Serret Frenet reference frame with the curvilinear abscissa, the width of pipe, the B-coordinate of any particle.

• : width of the pipe at altitude with (resp. ) is the Y-coordinate of right (resp. left) boundary point at altitude

• : angle

• : cross-section area

• : radius of the cross-section

• : outward normal vector to the wet part of the pipe

• n: outward normal vector at the boundary point in the -plane defined below

#### Notations concerning the free surface (FS) part

• : wet area

• : discharge

• : free surface cross section

• : physical water height

• : -coordinate of the water level, : width of the free surface

• : outward B-normal vector to the free surface

• : density of the water at atmospheric pressure

#### Notations concerning the pressurized part

• : pressurized cross section

• : density of the water

• : water compressibility coefficient

• : sonic speed

• : FS equivalent wet area

• : FS equivalent discharge

#### Notations concerning the Pfs model

• S: the physical wet area: if the state is free surface, otherwise

• : the coordinate of the water level: if the state is free surface, otherwise

#### Other notations

• Bold characters are used for vectors except for S

The classical shallow water equations are used to describe physical situations like rivers, coastal domains, oceans and sedimentation problems. These equations are obtained from the incompressible Euler system (see e.g. [2, 23]) or from the incompressible Navier-Stokes system (see for instance [10, 11, 17, 24]) by several techniques (e.g. by direct integration or asymptotic analysis). We adapt here the derivation in [3, 4] to get a new unidirectional shallow water model. We start from the D incompressible Euler equations where we neglect the acceleration following the -axis supposing the existence of a privileged main flow axis. We write then the Euler equations in the local Serret-Frenet reference frame in order to take into account the local effects produced by the changes of section and the slope variation. Then we derive a shallow water model by a formal asymptotic analysis (done in Subsection 2.3).

### 2.1 Incompressible Euler equations and framework

Let us consider the cartesian reference frame . In the corresponding coordinate system , the D incompressible Euler system writes:

 {div(ρ0U)=0∂t(ρ0U)+ρ0U⋅∇(ρ0U)+∇P=F (1)

where denotes the velocity with components , is the isotropic pressure tensor, the density of the fluid at atmospheric pressure and F is the exterior strength of gravity.

We close classically System (1) using a kinematic law for the evolution of the free surface: any free surface particle is advected by the fluid velocity U and on the wet boundary, we assume the no-leak condition where is the outward unit normal vector to the wet part of the pipe (see Fig. 2). We set the pressure to at the free surface.

We define the domain of the flow at time as the union of sections (assumed to be simply connected compact sets) orthogonal to some plane curve lying in to follow the privileged main flow axis. We choose the parametrization in the cartesian reference frame where k follows the vertical direction; is then the elevation of the point over the plane (see Fig. 1).

We define a local reference frame as follows: we introduce the curvilinear variable defined by:

 X=∫xx0√1+(b′(ξ))2dξ

where is an arbitrary abscissa. We set and we denote by the B-coordinate of any fluid particle in the Serret-Frenet reference frame at point with T the tangent vector N, the normal and B the binormal vector (see Fig. 1 and Fig. 3 for the notations). B is normal to in the vertical plane .

Then, at each point , is defined by the set:

 {(y,Z)∈R2;Z∈[−R(X),−R(X)+H(t,X)],y∈[α(X,Z),β(X,Z)]}

where denotes the radius, the physical water height at section . We denote (respectively Y-coordinate of the left (respectively right) boundary point of the domain at altitude , (see Fig. 3). We denote also by which is the -coordinate of the water level.

In the sequel, we will use a curvilinear map which will be an admissible transformation under the geometrical hypothesis on the domain:

Let be the algebraic curvature radius of the plane curve . We assume that:

 ∀x∈ΩF,|R(x)|>R(x).

### 2.2 Incompressible Euler model in the curvilinear coordinates

Following the work in [3, 4], we write System (1) in the Serret-Frenet reference frame at point by the transformation using the divergence chain rule lemma that we recall here:

###### Lemma 2.1

Let be a diffeomorphism and
the jacobian matrix of the transformation with determinant .

Then, for any vector field , one has:

 Jdiv(x,y,z)Φ=div(X,Y,Z)(JAΦ),

and, for any scalar function , one has:

 ∇(x,y,z)f=At∇(X,Y,Z)f,

where stands for the transpose of the matrix .

Let be the components of the velocity vector in the coordinates defined as where is the matrix

 \Uptheta=⎛⎜⎝cosθ0sinθ010−sinθ0cosθ⎞⎟⎠,

where we denote by the angle in the plane.

Using Lemma 2.1, the incompressible Euler system in the variables reads:

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩∂X(ρ0U)+∂Y(Jρ0V)+∂Z(Jρ0W)=0∂t(Jρ0U)+∂X(ρ0U2)+∂Y(Jρ0UV)+∂Z(Jρ0UW)+∂Xp=G1∂t(Jρ0V)+∂X(ρ0UV)+∂Y(Jρ0V2)+∂Z(Jρ0VW)+∂Y(Jp)=0∂t(Jρ0W)+∂X(ρ0UW)+∂Y(Jρ0VW)+∂Z(Jρ0W2)+J∂Z(p)=G2 (2)

where is the determinant of the transformation and

 G1=ρ0UWθ′(X)−Jgρ0sinθ,G2=−ρ0U2θ′(X)−Jgρ0cosθ.

The interested reader can find the details of the calculus in [3]. We have denoted by the derivative with respect to the space variable of any function .

On the wet boundary, the no-leak condition reads:

 (U,V,W)t.nwb=0. (3)
###### Remark 2.1

Notice that is the algebraic curvature of the axis at point and the function depends only on the variables . Moreover, under the hypothesis , we have in . Consequently, defines a diffeomorphism and thus the performed transformation is admissible.

In this section, we perform a formal asymptotic analysis on System (2). According to the work in [3, 17, 24], the shallow water equations can be obtained from the incompressible Navier-Stokes equations with particular boundary conditions. Here, we perform this analysis directly on the incompressible Euler system in order to get for some small parameter .

Let us introduce the usual small parameter where (the height) and (the length) are two characteristics dimensions along the B and T axis respectively. Moreover, we assume that the characteristic dimension along the j axis is the same as for the k axis. We introduce the other characteristics dimensions for time, pressure and velocity respectively and the dimensionless quantities as follows:

 ˜U=U/¯¯¯¯U,˜V=ϵV/¯¯¯¯U,˜W=ϵW/¯¯¯¯U,
 ˜X=X/L,˜Y=Y/H,˜Z=Z/H,˜p=p/P,˜θ=θ,˜ρ=ρ0.

In the sequel, we set and (i.e. we only consider laminar flows).

Under these hypotheses, we have . Thus, the rescaled System (2) reads:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂˜X˜U+∂˜Y(˜J˜V)+∂˜Z(˜J˜W)=0∂˜t(˜J˜U)+∂˜X(˜U2)+∂˜Y(˜J˜U˜V)+∂˜Z(˜J˜U˜W)+∂˜X˜p=G1ϵ2(∂˜t(˜J˜V)+∂˜X(˜U˜V)+∂˜Y(˜J˜V2)+∂˜Z(˜J˜V˜W))+∂˜Y(˜J˜p)=0ϵ2(∂˜t(˜J˜W)+∂˜X(˜U˜W)+∂˜Y(˜J˜V˜W)+∂˜Z(˜J˜W2))+˜J∂˜Z(˜p)=G2 (4)

where

 G1=ϵ˜U˜W˜κ(˜X)−sin˜θFr,L2−˜ZFr,H2(cos˜θ)′,
 G2=−ϵ˜U2˜ρ(˜X)−cos˜θFr,H2+ϵκ(X)˜Z˜Jcos˜θFr,H2,

is the Froude number along the T axis and the B or N axis where is any generic variable equal to or .

Formally, when vanishes, System (4) reduces to:

 ∂˜X˜U+∂˜Y(˜V)+∂˜Z(˜W) = 0 (5) ∂˜t(˜U)+∂˜X(˜U2)+∂˜Y(˜U˜V)+∂˜Z(˜U˜W)+∂˜X˜p = −sin˜θFr,L2 (6) −˜ZFr,H2(cos˜θ)′ ∂˜Z(˜p) = −cos˜θFr,H2 (7)

Let us introduce the conservative variables and representing respectively the wet area and the discharge defined as:

 A(t,X)=∫Ω(t,X)dYdZ,Q(t,X)=A(t,X)¯¯¯¯U

where is the mean value of the velocity :

 ¯¯¯¯U(t,X)=1A(t,X)∫Ω(t,X)U(t,X)dYdZ.

We integrate the preceding system (5-6-7) along the cross-section with the approximation and . Then, returning to the physical variables, the free surface model, that we call FS-model, reads:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩∂tA+∂XQ=0∂tQ+∂X(Q2A+gI1(X,A)cosθ)=gI2(X,A)cosθ−gAsinθ−gA¯¯¯¯Z(X,A)(cosθ)′ (8)

where and are respectively the classical term of hydrostatic pressure and the pressure source term defined by:

 I1(X,A)=∫h−R(h−Z)σdZ and I2(X,A)=∫h−R(h−Z)∂XσdZ

which are obtained from the integration of the pressure term in Equation (5) with (obtained from equation (7)).

In these formulas is the width of the cross-section at position and at height . The additional term is defined by . It is the -coordinate of the center of mass:

 ¯¯¯¯Z=∫Ω(t,X)ZdYdZ=∫h(t,X)−R(X)∫β(X,Z)α(X,Z)ZdYdZ=∫h(t,X)−R(X)Zσ(X,Z)dZ.

In System (8), we may add a friction term to take into account the dissipation of energy. We have chosen this term as the one given by the Manning-Strickler law (see e.g. [29]):

 Sf(A,U)=K(A)U|U|.

The term is defined by: , is the Strickler coefficient of roughness depending on the material, is the hydraulic radius and is the perimeter of the wet surface area (length of the part of the channel’s section in contact with the water).

## 3 Formal derivation of the P-model for pressurized flows

In this section, we present a new set of unidirectional shallow water like equations to describe pressurized flows in closed non uniform water pipes. This model is constructed to be coupled in natural way with the FS-model (8). Starting from the D compressible Euler equations in cartesian coordinates,

 ∂tρ+div(ρU)=0, (9)
 ∂t(ρU)+div(ρU⊗U)+∇p=F, (10)

where and denotes the velocity with components and the density respectively. is the scalar pressure and F the exterior strength of gravity.

We define the pressurized domain of the flow as the continuous extension of (see Subsection 2.1) defined by some plane curve with parametrization in the cartesian reference frame ; we recall that is then the elevation of the point over the plane (see Fig. 1). The curve may be, for instance, the axis spanned by the center of mass of each orthogonal section to the main mean flow axis, particularly in the case of a piecewise cone-shaped pipe. Notice that we consider only the case of infinitely rigid pipes, thus the sections are only -dependent.

We then write Equations (9-10) in the coordinates introduced in Subsection 2.1. As we want a unidirectional model, we suppose that the mean flow follows the -axis. To this end, we neglect the second and third equation for the conservation of the momentum.
By a straightforward computation, the mass and the first momentum conservation equation in the coordinates reads:

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩∂t(Jρ)+∂X(ρU)+∂Y(ρJV)+∂Z(ρJW)=0∂t(JρU)+∂X(ρU2)+∂Y(ρJUV2)+∂Z(ρJUW)+∂Xp=−ρJgsinθ+ρUW(cosθ)′ (11)

Applying the same asymptotic analysis developed in Subsection 2.3, Equations (9-10) read:

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩∂t(ρ)+∂X(ρU)+∂Y(ρV)+∂Z(ρW)=0∂t(ρU)+∂X(ρU2)+∂Y(ρUV)+∂Z(ρUW)+∂Xp=−ρgsinθ−gZ(cosθ)′ (12)

We choose the linearized pressure law:

 p=pa+ρ−ρ0βρ0 (13)

(see e.g. [29, 30]) in which represents the density of the fluid at atmospheric pressure , is some function set to zero and the water compressibility coefficient (equal to in practice). The sonic speed is then given by and thus .
For , is the outward unit vector at the point in the -plane and m stands for the vector (as displayed on Fig. 3).

Following the section-averaging method performed in Subsection 2.3, we integrate System (12) over the cross-section . Noting the averaged values over by the overlined letters (except ), and using the approximations the shallow water like equations read:

 ∂t(¯¯¯ρS)+∂X(¯¯¯ρS¯¯¯¯U) = ∫∂Ωρ(U∂X%m−V).nds (14) ∂t(¯¯¯ρS¯¯¯¯U)+∂X(¯¯¯ρS¯¯¯¯U+c2¯¯¯ρS) = −g¯¯¯ρSsinθ+c2¯¯¯ρS′ − g¯¯¯ρS¯¯¯¯Z(cosθ)′ + ∫∂ΩρU(U∂Xm−V).nds

where is the velocity in the -plane. We denote by the area of the cross-section of the pipe at position .

The integral terms appearing in (14) and (3) vanish, as the pipe is infinitely rigid, i.e. (see [8] for the dilatable case). It follows the non-penetration condition (see Fig. 4):

 ⎛⎜⎝UVW⎞⎟⎠.nwb=0.

Omitting the overlined letters (except ), we introduce the conservative variables

 A=ρρ0S the FS equivalent wet area (16) Q=AU the FS \emph{equivalent discharge}. (17)

and dividing Equations (14)-(3) by we get:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂t(A)+∂X(Q)=0∂t(Q)+∂X(Q2A+c2A)=−gAsinθ−gA¯¯¯¯Z(X,S)(cosθ)′+c2AS′S (18)

As introduced previously for the FS-model in Section (2.3), we may introduce the friction term given by the Manning-Strickler law (see e.g. [29]):

 Sf(S,U)=K(S)U|U|

where is defined by: , is the Strickler coefficient of roughness depending on the material and is the hydraulic radius where is the perimeter of the wet surface area (length of the part of the channel’s section in contact with the water, equal to in the case of circular pipe).

This choice of variables is motivated by the fact that this system is formally close to the FS-model (8) where the terms , , are respectively the counterparts of , , in System (18). Let us remark that the term is continuous through the change of state (pressurized to free surface or free surface to pressurized state) when the same curve plane is chosen (in practice, the main axis of the pipe). Then, we are motivated to connect “continuously” System (8) and (18) through transition points (through the change of state) by defining a continuous pressure law. It leads to a “natural” coupling between the pressurized and free surface model as we will see in Section 4.

## 4 The Pfs-model

The formulations of the FS-model (8) and P-model (18) are very close to each other. The main difference comes from the pressure law. In order to build a coupling between the two models, we have to define a pressure that ensures its continuity through transition points in the same spirit of [7]. As pointed out in the previous section, we will use the common couple of unknowns and the same plane curve (see Remark 4.1) to get a continuous model for mixed flows.

###### Remark 4.1

The plane curve with parametrization is chosen as the main pipe axis in the axisymmetric case. Actually this choice is the more convenient for pressurized flows while the bottom line is adapted to free surface flows. Thus we must assume small variations of the section ( small) or equivalently small angle as displayed on Fig. 4.

We introduce a state indicator (see Fig. 5) such that:

 E={1 if the state is pressurized: (ρ≠ρ0)0 if the state is free surface: (ρ=ρ0). (19)

Next, we define the physical wet area S by:

 S=S(A,E)={S if E=1A if E=0 (20)

and a modified pressure law (see Fig. 5) which ensures its continuity through the change of state by:

 p(X,A,E)=c2(A−S)+gI1(X,S)cosθ. (21)
###### Remark 4.2
• Indeed, when a change of state occurs we have:

 lim\lx@stackrelA→SASp(X,A,E)=gI1(X,S)cos(θ)

which ensures the continuity of the pressure.

• The flux gradient is discontinuous through the change of state since

 ∂F∂A(A,Q,0)=g∂∂AI1(X,A)cosθ≠c2=∂F∂A(A,Q,1).

Finally, from the P-model (18), the FS-model (8), the definition of (19), the definition of S (20) and the pressure law (21), the PFS-model for unsteady mixed flows can be simply expressed into a single formulation as:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂t(A)+∂X(Q)=0∂t(Q)+∂X(Q2A+p(X,A,E))=−gAb′+Pr(X,A,E)−G(X,A,E)−K(X,A,E)Q|Q|A (22)

where , , and denotes respectively the friction, the pressure source and the geometry source term defined as follows:

 Pr(X,A,E)=c2(AS−1)S′+gI2(X,S)cosθ with I2(X,S)=∫H(S)−R(X)(H(S)−Z)∂Xσ(X,Z)dZ,G(X,A,E)=gA¯¯¯¯Z(X,S)(cosθ)′,K(X,A,E)=1K2sRh(S)4/3

and stands for . represents the -coordinate of the water level:

 H=H(S)={h(A) if E=0R(X) if E=1. (23)
###### Remark 4.3 (Both models are recovered)

Setting in System (22), we obtain obviously the free surface model (8). For all pressurized states, when , the pressure law (21) reads, for instance, in the case of circular pipe:

 c2(A−S)+gI1(X,S)cosθ=c2(A−S)+gπR3cosθ

which is not exactly the pressure law of the P-model (18). Indeed, the derivation of the P-model is done with the linearized pressure law (13) (see Section 3) with . Thus, the property of the continuity of models (18)-(8) through a change of state is obtained if and only if is chosen as which is the hydrostatic pressure corresponding to a full section.

The PFS-model (22) satisfies the following properties:

###### Theorem 4.1
1. The right eigenvalues of System (22) are given by:

 λ−=U−c(A,E),λ+=U+c(A,E)

with where is the width of the free surface (see Fig. 3).

Then, System (22) is strictly hyperbolic on the set:

 {A(t,X)>0}.
2. For smooth solutions, the mean velocity satisfies

 ∂tU+∂X(U22+c2ln(A/S)+gH(S)cosθ+gb)=−gK(X,A,E)U|U|⩽0. (24)

The quantity is called the total head.

 u=0 and c2ln(A/S)+gH(S)cosθ+gb=0. (25)
4. It admits a mathematical entropy

 E(A,Q,E)=Q22A+c2Aln(A/% S)+c2S+gA¯¯¯¯Z(X,S)cosθ+gAb (26)

which satisfies the entropy relation for smooth solutions

 ∂tE+∂X((E+p(X,A,E))U)=−gAK(X,A,E)U2|U|⩽0. (27)

Notice that the total head and are defined continuously through the transition points.

###### Remark 4.4

The term is also called “corrective term” since it allows to write the Equations (24) and (27) with (26).

Proof of Theorem 4.1: the results (24) and (27) are obtained in a classical way. Indeed, Equation (24) is obtained by subtracting the result of the multiplication of the mass equation by to the momentum equation. Then multiplying the mass equation by and adding the result of the multiplication of Equation (24) by , we get:

 ∂t(Q22A+c2Aln(A/S)+c2S+gA¯¯¯¯Z(X,S)cosθ+gAb)+∂X((Q22A+c2Aln(A/S)+c2S+gA¯¯¯¯Z(X,S)cosθ+gAb+p(X,A,E))U)+c2(AS−1)∂tS=−gAK(X,A,E)U2|U|⩽0.

We see that the term is identically since we have when the flow is free surface whereas when the flow is pressurized. Moreover, from the last inequality, when , we have the classical entropy inequality (see [7, 8]) with :

 E(A,Q,E)=Q22A+gA¯¯¯¯Z(X,A)cosθ+gAb

while in the pressurized case, it is:

 E(A,Q,E)=Q22A+c2Aln(A/S)+c2S+gAb.

Finally, the entropy for the PFS-model reads:

 E(A,Q,E)=Q22A+c2Aln(A/% S)+c2S+gA¯¯¯¯Z(X,S)cosθ+gAb.

Let us remark that the term makes continuous through transition points and it permits also to write the entropy flux under the classical form .

## 5 Finite volume discretisation

In this section, we adapt the VFRoe scheme described in [12, 16, 7]. The new terms appearing in the PFS-model related to the curvature and the section variation are upwinded in the same spirit of [7]. The numerical scheme is adapted to discontinuities of the flux gradient occurring in the treatment of the transitions between free surface and pressurized states.

### 5.1 Discretisation of the space domain

The spatial domain is a pipe of length . The main axis of the pipe is divided in cells . denotes the timestep and we set .
The discrete unknowns are . For the sake of simplicity, the boundary conditions are not treated (the interested reader can find this treatment in details in [7]).

### 5.2 Explicit first order VFRoe scheme

We propose to extend the finite volume discretisation [7] to the PFS-model using the upwinding of the new source terms: the curvature and section variation of the pipe. In what follows, we do not write the dependency.

First, following Leroux et al. [18, 26] we use piecewise constant functions to approximate as well as the term and the cross section area . Adding the equations , and , the PFS-model can be written under a non-conservative form with the variable :

 ∂tW+∂XF(X,W)+B(X,W)∂XW=TS(W) (28)

where

 F(X,W)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝000QQ2A+p(X,A)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, TS(W)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝0000−gK(X,S)Q|Q|A