A model for unsteady mixed flows in non uniform closed water pipes and a wellbalanced finite volume scheme
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 PFSmodel 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 Roelike 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, wellbalanced 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 nonconservative 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 PFSmodel 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 FSmodel. In Section 3, we adapt the derivation of the FSmodel to derive the pressurized model, called Pmodel, from the D compressible Euler equations. Writing the source terms of these two models, P and FSmodel, 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 PFSmodel . 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 wellbalanced 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

: SerretFrenet 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 Bcoordinate of any particle.

: width of the pipe at altitude with (resp. ) is the Ycoordinate of right (resp. left) boundary point at altitude

: angle

: crosssection area

: radius of the crosssection

: 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 Bnormal 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
2 Formal derivation of the Fsmodel for free surface flows
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 NavierStokes 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 SerretFrenet 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:
(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 noleak 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:
where is an arbitrary abscissa. We set and we denote by the Bcoordinate of any fluid particle in the SerretFrenet 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:
where denotes the radius, the physical water height at section . We denote (respectively Ycoordinate 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:
2.2 Incompressible Euler model in the curvilinear coordinates
Following the work in [3, 4], we write System (1) in the SerretFrenet 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:
and, for any scalar function , one has:
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
where we denote by the angle in the plane.
Using Lemma 2.1, the incompressible Euler system in the variables reads:
(2) 
where is the determinant of the transformation and
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 noleak condition reads:
(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.
2.3 Formal derivation of the Fsmodel for free surface flows
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 NavierStokes 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:
In the sequel, we set and (i.e. we only consider laminar flows).
Under these hypotheses, we have . Thus, the rescaled System (2) reads:
(4) 
where
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:
(5)  
(6)  
(7) 
Let us introduce the conservative variables and representing respectively the wet area and the discharge defined as:
where is the mean value of the velocity :
We integrate the preceding system (567) along the crosssection with the approximation and . Then, returning to the physical variables, the free surface model, that we call FSmodel, reads:
(8) 
where and are respectively the classical term of hydrostatic pressure and the pressure source term defined by:
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 crosssection at position and at height . The additional term is defined by . It is the coordinate of the center of mass:
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 ManningStrickler law (see e.g. [29]):
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 Pmodel 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 FSmodel (8). Starting from the D compressible Euler equations in cartesian coordinates,
(9) 
(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 coneshaped pipe. Notice that we consider only the case of infinitely rigid pipes, thus the sections are only dependent.
We then write Equations (910) 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:
(11) 
(12) 
We choose the linearized pressure law:
(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 sectionaveraging method performed in Subsection 2.3, we integrate System (12) over the crosssection . Noting the averaged values over by the overlined letters (except ), and using the approximations the shallow water like equations read:
(14)  
where is the velocity in the plane. We denote by the area of the crosssection 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 nonpenetration condition (see Fig. 4):
Omitting the overlined letters (except ), we introduce the conservative variables
the FS equivalent wet area  (16)  
(17) 
and dividing Equations (14)(3) by we get:
(18) 
As introduced previously for the FSmodel in Section (2.3), we may introduce the friction term given by the ManningStrickler law (see e.g. [29]):
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 FSmodel (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 Pfsmodel
The formulations of the FSmodel (8) and Pmodel (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:
(19) 
Next, we define the physical wet area S by:
(20) 
and a modified pressure law (see Fig. 5) which ensures its continuity through the change of state by:
(21) 
Remark 4.2

Indeed, when a change of state occurs we have:
which ensures the continuity of the pressure.

The flux gradient is discontinuous through the change of state since
Finally, from the Pmodel (18), the FSmodel (8), the definition of (19), the definition of S (20) and the pressure law (21), the PFSmodel for unsteady mixed flows can be simply expressed into a single formulation as:
(22) 
where , , and denotes respectively the friction, the pressure source and the geometry source term defined as follows:
and stands for . represents the coordinate of the water level:
(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:
which is not exactly the pressure law of the Pmodel (18). Indeed, the derivation of the Pmodel 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 PFSmodel (22) satisfies the following properties:
Theorem 4.1

The right eigenvalues of System (22) are given by:
with where is the width of the free surface (see Fig. 3).
Then, System (22) is strictly hyperbolic on the set:

For smooth solutions, the mean velocity satisfies
(24) The quantity is called the total head.

The still water steady state reads:
(25) 
It admits a mathematical entropy
(26) which satisfies the entropy relation for smooth solutions
(27)
Notice that the total head and are defined continuously through the transition points.
Remark 4.4
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:
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 :
while in the pressurized case, it is:
Finally, the entropy for the PFSmodel reads:
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 PFSmodel 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 PFSmodel 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.