Upwind Interface Source method
The finite volume method is possible for treating numerically hyperbolic systems of conservation laws, it is robust and presents the advantage to be conservative (we refer to  for a survey of its properties). We look at the semi-discrete scheme, called method of lines, where only a space discretization of equation (1.1) is performed.
We consider a mesh of R made up of cells Ci, with center xi, i ∈ Z and nonuniform length Δxi; we denote xi+ 1 the cell interfaces, so that the control 2 x +x volume can be identified as Ci = [xi− , xi+ ) and xi = I− 1 I+ 1 . Then, we construct a piecewise constant representation of the function z(x) on the mesh, whose coeﬃcients are zi = Δ1xI CI z(x) dx for example.
In this context, the discrete unknowns are expected to be approximations of the mean values of u on the mesh cells (the conservative quantities are cell centered), ui(t) ≈ Δxi CI u(t, x) dx, t ∈ R+, i ∈ Z, while the numerical fluxes are defined at the interfaces of the mesh.
To correctly treat the source term is more diﬃcult than it seems and centered schemes give unsatisfactory results, as it is well reported in the literature: a direct discretization of the source term by cell-averages, for instance, can not preserve the steady state solutions. A better approach is based on the Upwind Interface Source method (upwinding of external terms was originally formulated by Roe ), where the source term is also upwinded at the interfaces.
The general finite volume scheme for equation (1.1) can thus be written in the explicit form Δxi dui + (Ai+ 1 − Ai− 1 )+B+ + B− = 0, (2.1) dt 2 2 i− i+ 2 2 dropping the time dependence of the numerical functions for simplicity. We proceed to explain the notation in the previous formula. We first intro-duce the discrete fluxes Ai+ 1 = A(ui, ui+1), A ∈ C1, (2.2) where the numerical function A is chosen as a consistent approximation of the analytical flux, A(u, u) = A(u). (2.3)
Because of the choice of particular source terms (1.4), the function z(x) is defined up to a constant. Therefore, without loss of generality, we suppose the source term is discretized at the cell interfaces by means of functions B± 1 = B±(ui, ui+1, zi+1 − zi), B± ∈ C2, (2.4) i+ 2 and, in view of (1.4), it is natural to impose B+(u, v, 0) = B−(u, v, 0) = 0. (2.5)
The last condition refers to the interpretation of the numerical solver (2.1), applied to the model (1.1)-(1.4). According to the finite volume formalism, we can identify Bi− 21 + Bi+ 21 ≈ CI z (x)b(u)dx; (2.6) formally, this leads to deduce Bi−+ 21 + Bi++ 21 ≈ xI I+ 2 z (x)b(u)dx + x 1 z (x)b(u)dx, (2.7) x 1 xI+1 that is a way to perform an interfacial approximation of zero order terms. Such a discretization is also upwinded, in the sense that Bi−+ 12 represents the contribution of the waves coming from the left of the interface xi+ 12 and moving towards the cell Ci if they have a nonpositive velocity, while Bi++ 12 represents the waves moving forwards from the right of xi+ 1 and counted only if they have nonnegative velocity. Notice that, when the2 problem becomes homogeneous (for example, z (x) = 0 in equation (1.4), motivated by the analogy with the Saint-Venant model), this scheme reduces to the usual finite volume approximation for a scalar conservation law.
We observe that all what is stated in this section and in the following is also valid for a fully explicit scheme (obtained, for instance, using a standard forward Euler method for the time discretization), Δxi (un+1 − un) + (An − An 1 ) + Bn,+1 + Bn,−1 = 0, Δt i i i+ i− i− i+ where we introduce a time-step Δt and set tn = nΔt, n ∈ N. We then have to consider an additional restriction on the size of the ratio , the usual CFL condition, to guarantee numerical stability.
We define general conditions on the discretizations Bi±+ 12 so that the numerical scheme (2.1) preserves the steady state solutions. Note that all the methods developed in the references mentioned above are compatible with the formalism introduced in Section 2 and can be put in form (2.1), as we will do later for some particular cases.
By integrating the stationary equation associated with (1.1)-(1.4), we ob-tain the algebraic relation (1.6) for smooth steady state solutions. A discrete version is given by D(ui) + zi = Cst, ∀ i ∈ Z. (3.1)
We consider appropriate hypotheses on D, to ensure the existence of a unique Lipschitz continuous solution of that problem, namely that D is strictly monotonic. This assumption is restrictive and not always satisfied in realistic situations, but it is usually made (we are aware of that does use it in ).
For all initial data defined according to (3.1), a solver preserving steady states has to verify (Ai+ 1 − Ai− 1 ) + Bi− 1 + Bi+ 1 = 0.
This last statement can be formulated in terms of numerical functions, thanks to definition (2.2) and (2.4), so that it writes A(ui, ui+1) − A(ui−1, ui) + B+(ui−1, ui, zi − zi−1) + B−(ui, ui+1, zi+1 − zi) = 0, for all ui−1, ui, ui+1, such that D(uj )+zj = H, j = i−1, i, i+1. In particular, we choose alternatively ui−1 = ui and ui = ui+1 (then we deduce, respectively, zi−1 = zi and zi = zi+1), to obtain equivalent conditions at the interfaces, also exploiting properties (2.3) and (2.5),
A(ui, ui+1) − A(ui) + B−(ui, ui+1, zi+1 − zi) = 0,
A(ui+1) − A(ui, ui+1) + B+(ui, ui+1, zi+1 − zi) = 0.
We summarize the previous statements in the following proposition.
Lemma 3.1. A numerical scheme in form (2.1)-(2.5) is well-balanced, i.e. it preserves the steady state solutions (3.1), if and only if the equalities
A(u, v) − A(u) + B−(u, v, z+ − z−) = 0, (3.2)
A(v) − A(u, v) + B+(u, v, z+ − z−) = 0, (3.3)
hold true, for all u, v, z−, z+ such that
D(u) + z− = D(v) + z+. (3.4)
We call well-balanced or Steady State Preserving schemes the numerical solvers for the problem (1.1)-(1.2) which satisfy those conditions.
We present some numerical schemes to which Lemma 3.1 applies. We check these approaches enable to preserve the discrete steady state solutions, ac-cording to the result stated above.
B.P.V. method In , the authors introduce their solver in a compact form, taking into account the source term directly in the definition of the numerical fluxes, Δxi dudti + (A−i+ with A−i+ 12 = A(ui, u−i+1), 1 −A+ 1)=0,
The numerical flux used in  is given by a standard Engquist-Osher function, but one readily finds out that similar methods can be formulated with any consistent flux function A. The points u−i+1 and u+i−1 are defined by means of the relations
D(u−i+1) + zi = D(ui+1) + zi+1,
D(u+i−1) + zi = D(ui−1) + zi−1.
To reproduce this scheme in form (2.1), we identify Bi− 1 = A(ui−1, ui) − A(ui−1, ui), Bi+ 1 = A(ui, ui+1
For a steady state, we immediately deduce from (3.1), (3.7) and (3.8) that u−i+1 = ui, u+i−1 = ui, then resulting in the conditions (3.2) and (3.3).
This method extends to more general classes of function D, such as quadratic functions, and it also applies to hyperbolic systems of conservation laws endowed with a kinetic interpretation (see ). We remark that, combined with specific approximate Riemann solvers, the algorithm (3.5)-(3.6) can be interpreted as the well-balanced scheme de-rived by Greenberg and LeRoux  or by Gosse and LeRoux , for which the condition (3.2)-(3.4) is verified.
The quasi-steady wave-propagation algorithm The basic idea of the method developed by LeVˆeque  is to introduce a new Riemann problem in the center of each mesh cell, with values u−i on the left half of the cell and u+i on the right half, whose flux diﬀerence exactly cancels the eﬀect of the source term. These artificial states are defined so that the cell-average is preserved, ui− = ui − δi, ui+ = ui + δi, 1 (ui− + ui+) = ui; (3.9) moreover, if δi is chosen according to the in-cell balance condition A(ui+) − A(ui−) + zib(ui)Δxi = 0, (3.10) then the jumps occurring at the cell interfaces will correspond to perturba-tions from the steady states. Note that (3.10) represents a discrete version of the stationary problem associated with equation (1.1). The explicit formula for the scheme thus obtained looks like the classical Godunov solver,
Δxi dui + Δ+A(ui+−1, ui−) + Δ−A(ui+, ui−+1) = 0, dt
Δ+A(ui+−1, ui−) = A(ui−) − A(ui∗− 1 ), (3.11)
Δ−A(ui+, ui−+1) = A(ui∗+ 1 ) − A(ui+), (3.12)
where u∗i+ 12 now denotes the solution to the modified Riemann problem at the cell interfaces, between values u+i and u−i+1. If the solution we are looking for is quasi-steady then we deduce from (3.9) and (3.10) that u+i ≈ u−i+1, as δi tends to 0, so that the steady states are asymptotically preserved.
As for previous examples, this method can extend to any consistent nu-merical flux function A, by rewriting the flux diﬀerences (3.11)-(3.12) in the more general form Δ+A(u+ , u−) = A(u−) − A(u+ , u−), i−1 i i i−1 i
Jin’s formulas A simple scheme for handling hyperbolic systems of con-servation laws with source terms is proposed in , which preserves the steady state solutions exactly at the cell interfaces. For methods based on a generalized Riemann solver, it takes the form (2.1) and the source term is discretized by + Ai+ 1 − Ai− 1 Bi− 1 + Bi+ 1 = (zi+ 1 − zi− 1 ), (3.13) using interface values Ai+ 12 = A(ui+ 12 ) and Di+ 12 = D(ui+ 12 ), rather than the cell-averages. Consequently, if one gets a steady state at the interfaces, Di+ 1 + zi+ 1 = Cst, ∀i ∈ Z, a direct computation leads to verify that it is preserved, as we have Ai+ 1 − Ai− 1 (Ai+ − Ai− ) + (zi+ − zi− )=0.
A more generic scheme, again proposed in , is defined by dui bi− 1 + bi+ 1 Δxi + (Ai+ 1 − Ai− 1 ) + 2 2 (zi+ 1 − zi− 1 )=0. (3.14) dt 2
Although it is not possible to derive an explicit expression of D for a gen-eral flux function A, some applications considered by the author (shallow water equations, for instance) show this method yields formally second order approximations to the steady states at the interfaces, as suggested by an asymptotic expansion of (3.14).
The numerical discretizations formulated by Jin are called Steady State Capturing schemes, that is a weaker definition since only the interface values are preserved. According to the idea to process the source term by making use explicitly of relations on the steady states, a Steady State Preserving variation of (3.13), which agrees with the general formalism (2.1)-(2.4), is given by
1 = A(ui−1, ui) − A(ui)(zi − zi−1),
2 D(ui−1) − D(ui)
1 = A(ui, ui+1) − A(ui)(zi+1 − zi).
2 D(ui+1) − D(ui)
Again we can readily check the condition (3.2)-(3.4) for this method.
In order to investigate theoretical properties of the Upwind Interface Source method, a crucial question to discuss is that equation (2.1) verifies the consistency with the continuous equation (1.1)-(1.4).
We indicate a rigorous definition of consistency, which also results to be satisfied by well-balanced schemes.
Definition 4.1. A numerical scheme in form (2.1)-(2.5) is said to be con-sistent with (1.1) if the following limit is verified, locally uniformly in u,
B+(u, u, λ) + B−(u, u, λ)
lim = b(u). (4.1)
λ → 0 λ
We point out that the above definition of consistency for the source term, in finite volume sense, does not imply that the consistency error vanishes just as for the flux terms. Indeed, because of the choice of an arbitrary nonuniform spatial mesh, the space-step Δxi could be very diﬀerent from the length of an interfacial interval Δxi+ 12 = |xi+1 − xi| = Δxi/2 + Δxi+1/2; therefore, the corresponding interfacial discretizations (2.6) and (2.7) may have very diﬀerent values. The condition (4.1) we have established is closer to (2.7), which is the most appropriate interpretation of the discrete source term for the general method illustrated in this paper.
As it will be seen clearly in next section, consistency plays an important role to achieve convergence properties of a numerical solver, in particular to prove that the strong limit of discrete approximations (as the mesh is refined) is the suitable weak solution of the continuous problem.
The following result guarantees consistency for the numerical schemes described in Section 3.
Lemma 4.2. We assume D is monotonic. Let a numerical solver for the system (1.1)-(1.4) satisfy the conditions (3.2)-(3.4) in Lemma 3.1, then the property (4.1) is verified. In other words, all Steady State Preserving schemes are consistent.
Proof. We perform a Taylor expansion of the relation (3.4) and we deduce that, for some ξ ∈ (u, v),
D (ξ)(u − v) = z+ − z−. (4.2)
After adding equality (3.2) to (3.3), thanks to the definition (1.3) and (1.6), this leads to B+(u, v, z+ − z−) + B−(u, v, z+ − z−) = A(u) − A(v) = a(ζ) (z+ − z−), (4.3) D (ξ) for some ζ ∈ (u, v). We also note that the regularity assumed for the numer-ical functions enables to perform the general approximations B±(u, v, z+ − z−) = B±(u, u, z+ − z−)
Table of contents :
1 Présentation du problème
1.1 Les lois de conservation hyperboliques avec terme source g´eom´etrique
1.2 Le syst`eme des ´equations de Saint-Venant
1.3 La question num´erique des solutions stationnaires
2 L’approche numérique
2.1 La m´ethode ”Upwind Interface Source” pour les termes sources
2.2 L’extension aux discr´etisations d’ordre deux
2.3 Les r´esultats de convergence
3 La méthode cinétique
3.1 Interpr´etation cin´etique du syst`eme de Saint-Venant
3.2 La formulation g´en´erale du sch´ema cin´etique
3.3 Le sch´ema cin´etique pour la m´ethode ”Upwind Interface Source”
4 Conclusions et perspectives
4.1 Syst`emes modifi´es et comparaisons exp´erimentales
4.2 Introduction `a la page web
4.3 Une application aux ´ecoulements granulaires
1 La méthode ”Upwind Interface Source” pour les lois de conservation hyperboliques
1 Convergence of the Upwind Interface Source method for hyperbolic conservation laws
2 Upwind Interface Source method
5 ALax-Wendroff type convergence theorem
2 First and second order error estimates for the Upwind Interface Source method
1.1 Formalism of the Upwind Interface Source method
1.2 What is a second order scheme for the Upwind Interface Source method ?
1.3 Convergence and error estimates
2.1 Stability estimate
2.2 Consistency estimate
2.3 Proof of Theorem 1.2
3.1 Stability estimate
3.2 Consistency estimate
3.3 Proof of Theorem 1.3 and Theorem 1.4
4 Remarks and numerical evidence
2 L’approximation num´erique des ´equations de Saint- Venant et applications aux ´etudes exp´erimentales
3 A kinetic scheme for the Saint-Venant system with a source term
2 Preliminaries about the Saint-Venant equations
2.1 Properties of the system
2.2 Kinetic approach
3.1 The formulas
3.2 Properties of the numerical scheme
4 Numerical implementation
4.1 Computation of the integrals
4.2 Some numerical tests
4 Second order approximation of the viscous Saint-Venant system and comparison with experiments
2 Formalism of the numerical method
4 Experimental configuration and numerical results
5 Numerical modelling of avalanches based on Saint-Venant equations using kinetic scheme
3 Flow and friction law
3.1 Simple friction law
3.2 Flow variable friction law
4 Numerical model
4.1 Finite volume method
4.2 Kinetic formulation
6.1 Curvature effects
6.2 The Coulomb friction law
6.3 Pouliquen’s friction law
6.4 Mass stopping