Sediment transport modelling : Relaxation schemes for Saint-Venant – Exner and three layer models

In this note we are interested in the modelling of sediment transport phenomena. We mostly focus on bedload transport and we do not consider suspension sediment processes. We first propose a coupled numerical scheme for the classical Saint-Venant – Exner model. It is based on a relaxation approach and it works with all sediment flux function. We exhibit that this coupled approach is more stable than the splitting approach that is mostly used in industrial softwares. Then we derive an original three layers model in order to overcome the difficulties that are encountered when using the classical Exner approach and we write a related relaxation model. 1. Sediment transport modelling : A major issue In rivers, mean sediment discharge may represent several hundred cubic meters of gravels or silt per year. Therefore, the sediments must be taken into account in order to predict the river bed evolutions. For hydroelectricity managers sediment transport modelling is a major issue. Dams stop water and consequently impact sediment transport, as well as sediment deposition may disturb water intakes or bottom gate opening, for example see Fig. 1. Besides, in some rivers like the Loire river, water intakes of nuclear power station may be protected against the deposition of sand. In order to understand sediment transport dynamics and to suggest managing solutions, EDF has been working on sediment transport modelling tools for bedload and suspension sediment transport. This work focuses on the modelling of bedload transport which refers to gravel transport and pushes aside the transport of fine sediments by suspension. Up to now, one very classical approach is to approximate the solid phase equation by a simplified one : the well-known Exner equation [10]. The Exner equation is obtained by writing a mass conservation on the solid phase in interaction with the fluid. There is no dynamic effect in the solid phase 1 BANG Project, INRIA-Paris-Rocquencourt & LAGA, University Paris 13 Nord, France, e-mail : [email protected] 2 LJLL, University Paris Diderot, France, e-mail : [email protected] 3 Lab. J.A. Dieudonné & EPU Sofia Nice, University of Nice, France, e-mail : [email protected] 4 EDF R&D & Saint-Venant Lab., France, e-mail : [email protected] 5 EDF R&D, France, e-mail : [email protected] 6 BANG Project, INRIA-Paris-Rocquencourt & CETMEF, France, e-mail : [email protected] 7 IAANS, University of Stuttgart, Germany, e-mail : [email protected] 8 LAMFA, University of Picardie, France, e-mail : [email protected] c © EDP Sciences, SMAI 2012 ha l-0 06 74 36 3, v er si on 1 27 F eb 2 01 2 2 ESAIM: PROCEEDINGS Figure 1. Example of silt deposition pattern in reservoirs (water flows from the back to the front). ρs(1− p) ∂zb ∂t + ∂Qs ∂x = 0, where zb is the bed elevation, Qs the bed load, p the porosity of the gravel bed. The bed load may be expressed by empirical formulae of the form Qs = Ag(u)|u|u (1) where u is the velocity in the fluid (Grass [13] formula) or Qs = f(τb) where τb is the boundary shear stress (see for example Meyer-Peter and Muller [18], Einstein [8] or Engelund and Fredsoe [9] formulas). The classical fluid model is the shallow water equations. It is coupled with the solid phase by the bottom level evolution. The coupled model stands 1 : ∂H ∂t + ∂Q ∂x = 0, (2) ∂Q ∂t + ∂ ∂x ( Q H + g 2 H ) = −gH ∂zb ∂x , (3) ρs(1− p) ∂zb ∂t + ∂Qs ∂x = 0, (4) where Q = Hu is the water discharge and H the water height. There are two major ways to solve this three-equation system : the coupled or the uncoupled approaches. The softwares developed at EDF-RD are based on an uncoupled resolution of the system. Firstly, the hydraulic part of the system is solved by using the software MASCARET [12] and then the computed fluid quantities are sent to the software COURLIS [3] that solves the Exner equation. The fluid part and the solid one are coupled through the time evolution of the bottom level. In the best way, the shallow-water equations and the Exner equation are coupled at each time step without local iterations. EDF actual modelling tools are used in complex 1For simplicity’s sake, equations are written for the one dimension problem in rectangular channels ha l-0 06 74 36 3, v er si on 1 27 F eb 2 01 2 ESAIM: PROCEEDINGS 3 cases, for example the draw-down of a reservoir with a steep slope and deep zones. Engineers are confronted with the limitation of the codes. Fig 2 is an example of instabilities observed in the case of a dam break calculation using an uncoupled scheme (0.3s after dam break). This result is in agreement with [7] where the authors show that, with a splitting method, the instabilities cannot be always avoided in supercritical regions. Thus alternative approaches must be searched. Figure 2. Dam break over a moveable bottom with a splitting approach Free surface (top) and Bottom topography (bottom) In this note, we investigate two different approaches. First we consider the coupled approach that consists in solving the system (2)-(4) at once with three unknowns which are H, Q and zb. We then have a hyperbolic system to solve and some specific finite volume schemes have been recently developed, see [1, 6]. For stiff cases this approach seems to be more robust than the uncoupled one [7]. Here we propose a relaxation approach [4] of the model that allows us to deal with a large class of sediment fluxes and we construct the related relaxation solver. Second we go one step further and we derive a new three layers model that is able to propose a more accurate modelling of the phenomenon than the classical Exner approach. The obtained model is quite similar to some models studied in [20]. We also propose a relaxation approach for this new model. In Section 2, we use a relaxation approach to develop a stable coupled numerical solver for the classical SWExner model (2)-(4) and we propose some numerical test cases in Section 3. Then in Sections 4 and 5 we derive the new three layers model and we propose a related uncoupled relaxation model. 2. A relaxation solver for the Saint-Venant – Exner model In this section, we consider the classical SW-Exner model (2)-(4) and we propose a relaxation approach for this model. Then we construct the related numerical solver, see also [17]. This work has two main motivations : first, to construct a stable coupled solver for the numerical simulation of model (2)-(4) ; second, to introduce the main ideas that will be used to derive the relaxation model associated to the three layers model that is derived in the following sections. The hyperbolic nature of the SW-Exner model (2)-(4) strongly depends on the formula that is chosen for the sediment flux Qs [7]. Moreover the computation of the eigenvalues (that are needed if we want to compute an approximate Riemann solver) of the system is also related to this formula. Here we overcome this difficulty by ha l-0 06 74 36 3, v er si on 1 27 F eb 2 01 2 4 ESAIM: PROCEEDINGS considering a relaxation model where we relax the sediment flux. We also relax the hydrostatic pressure as it is usually done for shallow water system. This method was first introduced in [16]. It consists in the introduction of auxiliary variables and leads to a larger but more suitable system. The method is well-described in [4] where a large number of classical hyperbolic solver are interpreted as relaxation schemes. In our case we finally obtain the following model, see also [17] ∂H ∂t + ∂Hu ∂x = 0 (5) ∂Hu ∂t + ∂ ∂x ( Hu + Π ) = −gH ∂zb ∂x (6) ∂Π ∂t + u ∂Π ∂x + a H ∂u ∂x = 1 λ (Π− gH 2 2 ) (7) ∂zb ∂t + ∂Q̃s ∂x = 0 (8) ∂Q̃s ∂t + ( b H2 − u ) ∂zb ∂x + 2u ∂Q̃s ∂x = 1 λ (Q̃s −Qs) (9) where λ is a small parameter and with Π and Q̃s the relaxed quantities for the fluid pressure and for the sediment flux. The parameters a and b have to be chosen sufficiently large to ensure the stability of the model. This system obviously tends (at least formally) to the original SW-Exner model (2)-(4) when λ tends to zero. The main advantage of the relaxation system (5)-(9) is that all the eigenvalues are easy to compute whatever the sediment flux Qs is and they are ordered whatever the choice of a and b is, provided a 6= b. Moreover all the fields are linearly degenerated and then the solution of the Riemann problem is (quite) easy to compute. We can write system (5)-(9) on the quasi-linear form ∂tVf +A(Vf )∂xVf = Sz + Sλ