Comptes Rendus Mécanique

. We present a (partial) historical summary of the mathematical analysis of ﬁnite di ﬀ erence and ﬁnite volume methods, paying special attention to the Lax–Richtmyer and Lax–Wendro ﬀ theorems. We then state a Lax–Wendro ﬀ consistency result for convection operators on staggered grids (often used in ﬂuid ﬂow simulations), which illustrates a recent generalization of the ﬂux consistency notion designed to cope with general discrete functions

The finite volume method (FVM) has been used for over 60 years in computational fluid mechanics, and more than 30 years in solid mechanics.However, while the finite element method, also introduced in mechanics, had already been the object of several mathematical works at the turn of this millenium, the American Mathematical Society's (AMS) 2000 Classification only mentions the term "finite volume" in sections 74S10 (Mechanics of deformable solids) and 76M12 (Fluid mechanics); it was not until 2010 that this term appeared in section 65 "Numerical analysis" (65M08, 65N08).
In some respects, the FVM is close to the finite difference method (FDM): for instance it has often been called "conservative finite difference method" by the hyperbolic numerical community, and "finite difference method" in the oil industry.It might be for this reason that the first attempts to show the convergence of the FVM were to try and copy the FDM framework.We show in the first section below that the famous Lax-Richtmyer theorem developed for linear FD schemes fails to give an adequate answer, even for a linear finite volume (FV) scheme.We then turn to the Lax-Wendroff theorem that gives two fundamental tools for the analysis of FV schemes, and recapitulate the main steps of the proof of convergence of the schemes that were derived in the 90's.In Section 4, we extend the Lax-Wendroff theorem to the case of staggered meshes and show how it can be used for the celebrated MAC grid.

The finite difference method and the Lax-Richtmyer theorem
In classical numerical analysis textbooks, we are taught that, in order to show the convergence of a finite difference (FD) scheme, one should show its stability and its consistency.The founding result in this regard is the Lax-Richtmyer theorem [1] due to P. D. Lax, who exposed it in a seminar at NYU in 1954.The so-called Lax equivalence theorem can be summarized as follows (see for example [2, Theorem 1.5.1]):Theorem 1.1 (Lax-Richtmyer).Consider a linear partial differential equation for which the initial value problem is well posed, and a FD scheme consistent for its approximation; then this scheme is convergent if and only if it is stable.
In some articles and textbooks (see for instance [3, p. 142]), maybe because of the name "equivalence theorem", and because uniform meshes are considered, the Lax-Richtmyer theorem is replaced by the equivalence "consistency + stability ⇐⇒ convergence" However, Theorem 1.1 does not state the equivalence (1), which, in fact does not hold in the general case of a PDE which is discretized with a non-constant space step, even in the linear case.
As an example, let us consider the approximation on R×]0, T [, where T > 0 is the final time of the study, of the linear transport equation for some given a > 0, and initial data u ini ∈ C ∞ c (R, R), with ∂ t (resp.∂ x ) the time (resp.space) partial derivative.Let ū(x, t ) = u ini (x − at ) be the exact (unique) solution of this problem.The FDM applied to (3) is classically defined by choosing a strictly increasing sequence (x i ) i ∈Z of real numbers, such that h := max i ∈Z (x i +1 − x i ) < ∞ and h = min i ∈Z (x i +1 − x i ) > 0, and a time step δt = T /N , for N ∈ N with N > 1 (here a constant time step is considered for simplicity).The initial data is discretized by defining the following quantities that depend on h implicitly: Since some upwinding is indeed necessary for stability purposes, the approximation of ∂ x u at the point x i is upwinded, so that the FD scheme reads In this context, the definition of the terms consistency, stability and convergence is the following: • Consistency: it requires two conditions: -the consistency of the discretization of the initial condition, that is lim -the consistency of the discretization of the PDE (2); setting t n = nδt for n ∈ N, it reads: • Stability: There exists C depending only on u ini (and thus not on h nor on δt ) such that • Convergence: Clearly, condition (4) ensures the consistency condition (6) on the initial condition.Moreover, the consistency condition (7) on the PDE can be obtained by Taylor expansions.Finally, the stability condition (8) can be shown under the CFL (Courant-Friedrichs-Lewy) condition aδt ≤ h.The conditions for the Lax-Richtmyer to hold are therefore satisfied, and the convergence of the scheme (4)-( 5) is thus proven.
Consider now a variant of the scheme (5) obtained by keeping (4), but replacing h i −1/2 by h i = (x i +1 − x i −1 )/2 in (5): In the specific case where x 2k+1 − x 2k = h/2 and x 2k+2 − x 2k+1 = h for all k ∈ Z, we get that h i = 3/4h for all i ∈ Z and the consistency property no longer holds.Therefore, if the equivalence (1) were true, the scheme (4)- (9) would not be convergent in the above sense.However, let us show that this scheme is in fact convergent in the same sense.Let us write the finite difference scheme obtained at the points x i = (x i + x i +1 )/2 instead of the points x i , defining the values u n i : and Then the scheme (10)- (11) again satisfies the consistency condition (6), and satisfies the consistency condition (7) since h i = x i − x i −1 .The stability condition (8) can again be shown under the CFL condition aδt ≤ h ≤ 3/4h.Hence the Lax-Richtmyer theorem yields the convergence of the values , the maximum principle applied to the difference u n i − u n i , solution to the equation obtained by subtracting (9) to (11), shows that max i ∈Z,n∈ 0,N −1 leading to the convergence of the scheme in the same sense as above.This example shows that on the one hand the direction ⇐ of (1) is not true, and that on the other hand the Lax-Richtmyer theorem cannot be applied directly to obtain the convergence of such a scheme, since the scheme is not consistent in the above defined sense.In fact, the scheme (9) is the 1D upwind FV scheme with the control volumes ](x i −1 + x i )/2, (x i + x i +1 )/2[, as shown in the next section.In the scheme (5), the partial derivative ∂ x u is upwinded.In the scheme (9), the unknown u itself is upwinded.Note that the equivalence (1) also does not hold in the case of an elliptic operator, see the example of a non-uniform 1D mesh in [4] or [5,Section 5.2].For the analysis of such schemes, other notions must be introduced, and it seems that Peter D. Lax first identified them, in collaboration with B. Wendroff.

The finite volume method and the Lax-Wendroff theorem
In a very famous article of 1960, P. D. Lax and B. Wendroff [6] consider discretization schemes for nonlinear hyperbolic systems of conservation laws and show that if a conservative scheme with consistent fluxes, in the sense that they define and which is stated below (see Section 2.2), converges a.e. and boundedly towards a limit, then this limit is necessarily a weak solution of the system.We call this property Lax-Wendroff (LW) consistency.The notions of flux consistency and flux conservativity highlighted in [6] are truly fundamental for the convergence analysis of the FVM for the hyperbolic equations considered in [6], as are their extensions to elliptic and parabolic conservation equations.In order to explain these terms, we consider, again in the 1D case, a general differential form of a conservation law, written on the whole space R and on a time interval ]0, T [ where 0 < T < +∞ is the final time: stating the conservation of the quantity u at each point x ∈ R and each time t ∈ ]0, T [, with F a vector function depending only on x and t through the unknown u.In addition to this equation, an initial condition is assumed to be given on u.Simple examples of such a conservation law include • the transport equation F(x, t ) = au(x, t ) (linear hyperbolic equation, introduced in the previous section), The FVM consists in approximating the integral form of the conservation law, that is to say the balance on the time-space cuboid ]x, x + δx[×]t , t + δt [ (for given δx > 0 and δt > 0), rather than the PDE itself (this corresponds to the way such an equation is derived from physical conservation principles).Note that with a non-zero source term on the right-hand side of ( 12), the equation is then usually renamed "balance law", and this modification poses no problem for a FV discretization.The integral form relative to the differential form (12) of a conservation law reads Let (]x i −1/2 , x i +1/2 [) i ∈Z be a family of intervals of R (also called control volumes or grid cells), with , and let δt = T /N , N > 1 (the time step could be taken non-constant).Let x i be some chosen point in the cell ]x i −1/2 , x i +1/2 [ (the choice of this point is constrained by the flux consistency property in the case of elliptic or parabolic problems, but not in the hyperbolic case, contrarily to the FD consistency in the sense of the previous section).The discrete unknowns, {u n i , i ∈ Z, n ∈ {0, . . ., N }} are expected to be approximations of u(x i , t n ) with t n = nδt .The integral form ( 13) is written on each control volume ]x i −1/2 , x i +1/2 [ and time interval ]t n , t n+1 [, leading to the following FV scheme (with an explicit time scheme for the sake of simplicity): where , and F n i +1/2 is the numerical flux, which will be expressed in terms of the discrete unknowns to yield a numerical approximation of is the numerical flux outgoing ]x i −1/2 , x i +1/2 [ to the right and its opposite is the numerical flux outgoing ]x i +1/2 , x i +3/2 [ to the left: this copies the situation of the exact flux F(x i +1/2 , t n ).This is the well-known "local conservativity" or "flux conservativity" property, which is important in physical applications, but also fundamental in the mathematical analysis of the FVM.Indeed, it is thanks to this property that we may hope to prove some convergence properties of the method, both for elliptic or parabolic type and for hyperbolic equations despite the loss of the consistency in the FD sense, as presented in the first section of this paper.

Flux conservativity
Writing the FV scheme in one dimension naturally ensures the numerical flux conservativity, since only one flux is defined: F i +1/2 at the interface x i +1/2 .In a multi-dimensional framework, (d = 2 or 3) the PDE ( 12) is written as ∂ t u + divF = 0, where F is a vector function of x and t and div the space divergence operator.The scheme ( 14) is now written for a control volume K : where |K | (resp.|σ|) is the volume or surface of K (resp.the surface or length of |σ|) and |σ|F n K ,σ is the numerical outgoing flux from K through the face σ; it is an approximation of the ougoing normal flux σ F(x, t n )•n K ,σ (where n K ,σ is the normal vector to σ outward K ), which is expressed in terms of the discrete unknowns (u n M ) M ∈M : In 2 or 3D, the numerical flux is defined on either side of the interface σ. Suppose that the interface σ separates the control volumes K and L, which we write as σ = K |L, then the flux conservativity reads

Flux consistency
Let us turn back to the 1D case, for ease of notations.The numerical flux is said to be consistent if for a sequence of time and space discretizations, indexed by m and such that h (m) → 0 and δt (m) → 0 as m → +∞, one has where F n i +1/2 is the quantity obtained from F n i +1/2 when replacing the discrete unknowns by the values of a regular function u: this definition of consistency is equivalent to the usual notion of consistency introduced by Lax (written here for a two-point scheme): to be Lipschitz continuous, or at least, "lip-diag" see [7,Remark 5.2].In the context of the heat equation (F (u) = −∇u), the numerical flux is consistent in the above sense, since, for a regular function u, and [.We notice that in this case the flux depends on the choice of the points x i .

Stability
There are different notions of stability of a numerical scheme.The notion which is of interest in the context of the convergence of a numerical scheme for a general, possibly nonlinear, PDE is an estimate on the approximate solutions, independently of the mesh.For instance, the L ∞ stability of a linear FD scheme for a linear elliptic equation may be obtained by writing the scheme in matrix form and by obtaining a bound of the infinity norm of the inverse of this matrix.Even though the FV approximate solutions are piecewise constant, the estimates are obtained in a norm which is in close relation with the one used for the estimates on the solutions of the continuous problem, and which depends of course on the considered problem; the notion of stability in the FVM is therefore linked to the stability of the continuous problem.Let us take two examples: • the heat equation on [0, 1] with homogeneous Dirichlet conditions: the natural norm for the continuous problem is L 2 (H 1 0 ), and the associated discrete norm corresponding to the choice ( 17) is L 2 (H 1 0,d ) with , where M is the number of control volumes and u 0 = u M +1 = 0 and with h i +1/2 = x i +1 − x i .We denote by u the weak solution of the heat equation and u app the solution of the timeimplicit scheme (implicit schemes are a natural choice for parabolic equations to avoid a condition of type δt ≤ C h 2 ).The function u app is a piecewise constant function which is equal to u n i on the cuboid ]x i −1/2 , x i +1/2 [×]t n , t n+1 [.In the continuous PDE setting, an L 2 (H 1 0 ) estimate on u is obtained by taking u as a test function in the weak formulation of the heat equation, and integrating by parts.Similarly, the L 2 (H 1 0,d ) estimate on u app is obtained by multiplying the i -th discrete equation by δt u n i , summing over i and n and performing discrete integrations by parts (obtained by changes of indices in the sums).An L 2 (L 2 ) estimate on the approximate solutions then follows with a discrete Poincaré inequality [5, Lemma 9.1].
• the transport equation on R n , with initial condition u ini ∈ L ∞ : the natural norm for the continuous problem is L ∞ .It is classical and easy to show that the scheme (14), with the upwind flux F i +1/2 = a u i , is stable, see for example [5,Lemma 20.1], under a CFL [8] condition.

The linear case: stability + conservativity + flux consistency =⇒ convergence
Convergence for linear operators are often obtained through error estimates (the same technique as in the proof of stability is applied to the error between the approximate solution and the exact one); the compactness analysis of sequences of approximate solutions is another means, which also gives the existence of a solution (see [8] for a seminal paper on this type of proof ).
Let us detail this second means, which extends to nonlinear operators.Consider a sequence of approximate solutions, on meshes whose space and time steps tend to 0. The general principle of proof of convergence is then the following.Thanks to the stability, a uniform (with respect to the discretization parameters) estimate on the approximate solutions holds in a Lebesgue space, and so there exists a sub-sequence of this sequence which converges weakly (or -weakly) in that same Lebesgue space.Each equation of the scheme is then multiplied by the interpolation of a regular test function and by the time step; the summation of all these equations over the time index and space indices, and an integration by parts (using the conservativity of the fluxes) are performed so as to shift the discrete derivatives from the discrete unknowns to the regular test function (the flux consistency which is used on the test functions is thus that given by the dual discrete operator).Since the problem is linear, the weak convergence of the approximate solutions suffices to pass to the limit in all terms; however the parabolic and hyperbolic cases exhibit some different difficulties.

Elliptic or parabolic case
The passage to the limit can be performed thanks to the consistency of the flux, see e.g.[5, Theorem 8.1] for the 1D case.In the multidimensional case, the consistency of the flux for the heat equation (Laplace operator) is obtained for meshes which respect an orthogonality condition [5, Definition 9.1].In this case, the resulting matrix is symmetric and the flux consistency is identical to that given by the dual discrete operator.Anisotropic operators and general meshes have been the object of several different works in the last decades, we refer to [9] for a review of these methods, several of them leading to non-symmetric matrices.In this latter case, one can either prove convergence by an error estimate, using the flux consistency of the primal discrete operator, or a compactness method, using the flux consistency of the dual discrete operator, i.e. on the test functions.In both cases the main difficulty is to establish the stability of the scheme.

Hyperbolic case
For stability reasons, the numerical flux is upwinded and introduces an error term whose convergence to 0 must be shown.This term is a sum of products of differences of the values of the solutions in neighbouring cells by discrete derivatives of the test function.On 1D or multiD Cartesian uniform or non-uniform meshes, if u ini ∈ L ∞ ∩ BV (an additional argument is used to handle the case u ini ∈ L ∞ ) this term is shown to tend to 0 thanks to a uniform BV estimate on the approximate solutions (and on the continuous solution), see [10,11] who consider the nonlinear case; this proof uses the "TVD" character (total variation diminishing) of the monotone schemes [12].Unfortunately, in the case of an unstructured mesh in multiD, even for a linear equation, even if u ini ∈ L ∞ ∩ BV , it can be shown that the upwind scheme is not TVD; a counterexample is given in [13].A "weak BV " estimate based on the numerical diffusion of the scheme is established therein in order to prove convergence.
Let us notice that in both the parabolic and hyperbolic cases, the Lax-Wendroff theorem is not used directly: • in the parabolic case, because the continuous flux function cannot be applied to the approximate solutions, • in the hyperbolic case, because the approximate solutions only converge weakly, whereas the Lax-Wendroff theorem supposes a strong convergence.
Nevertheless in both cases, the fundamental notions introduced in [6] are used, namely: • flux conservativity: it is the property that leads to a weak form of the FV scheme by shifting the discrete derivatives of the discrete unknown to the discrete derivatives of the interpolant of the test function, • flux consistency: it is the property that is used to prove the fact that a limit of approximate solutions is a weak solution thanks to the convergence of the discrete derivatives of the interpolant of the test function to the exact derivatives of the same test function.

The nonlinear case: more compactness needed
In the nonlinear case, weak convergence is not sufficient; indeed, if a sequence (u n ) n∈N converges only weakly in a Lebesgue space to a limit u, there is no reason why the sequence ( f (u n )) n∈N should converge to f (u), even weakly.

Elliptic and parabolic equations
In the elliptic or parabolic setting, one of the fundamental tools to obtain more compactness is Kolmogorov's compactness theorem, a consequence of which is that any bounded sequence of L p , 1 ≤ p < +∞ which is "equicontinuous on average" admits a convergent sub-sequence (see for example [14,Theorem 8.16]).Equicontinuity in mean amounts to showing that the difference between the function and its translates in time and space converges to 0 in the L p norm, uniformly with respect to the time and space step, see for example [5,Lemma 18.3] in the case of a nonlinear parabolic equation of Stefan type.Once the compactness of the sequence of approximate solutions in L 2 (L 2 ) is proven, we can then exhibit a sub-sequence tending to ū in L 2 (L 2 ).
By passing to the limit in the "very weak" form of the scheme (that is with the discrete divergence of the discrete normal derivatives of the test functions) we can then show, as in the linear case, that each term (in time and space) converges to the corresponding term in the "very weak" formulation of the continuous problem.In the case of the Stefan problem, namely ∂ t u − ∆ϕ(u) = 0, this gives the convergence of the approximate solutions to the exact solution if ϕ is an increasing Lipschitz continuous fucntion.There is however an additional difficulty if ϕ is only non-decreasing.In this case, it is possible to prove compactness in L 2 (L 2 ) of ϕ(u app ) (where u app is the approximate solution) but not of u app for which only an L 2 (L 2 ) bound holds.It is possible however to conclude using the Minty trick (see [5,Chapter 4]).

Hyperbolic equations
Now consider a nonlinear hyperbolic conservation law of the form there exists a unique entropy solution of this problem [15].In order to show that a scheme approximates this entropy solution, it is first shown that it satisfies a discrete entropy equation.The case of Cartesian meshes has been studied independently by Kuznetsov [10] and Crandall and Majda [11].As in the linear case, if u ini ∈ BV , a BV estimate on the approximate solutions holds, uniformly with respect to the space and time step; Helly's lemma, which is itself a direct consequence of Kolmogorov's compactness theorem may then be invoked to obtain the convergence of a subsequence of the approximate solutions in L 1 (L 1 ), and one can then use the Lax-Wendroff theorem (which generalizes easily to the entropy formulation) [11], [5, Section 21.5]).It is also possible to handle the case u ini ∈ L ∞ , using a contraction principle in L 1 for the exact solution and for the approximate solution, see [11] for instance.
In the general case of a non-Cartesian mesh, a suitable BV estimate seems out of reach, and the proof of convergence is performed with the following steps.
• Consider a sequence of space discretizations, indexed by m, and of time steps δt (m) ; for any m, the mesh size is defined as the maximum diameter of the cells and denoted by h (m) .Assume that h (m) → 0 and δt (m) → 0 as m → +∞.• Owing to the L ∞ estimate on the sequence of approximate solutions (u (m) ) m∈N , there exists a sub-sequence which converges in a "nonlinear weak sense", i.e. there exists This notion of convergence is equivalent to the convergence to a Young measure [16]; it may seem a little simpler to handle in the sense that it involves a function, µ, rather than a measure; however this function depends on an additional parameter α ∈]0, 1[, which we will have to get rid of in order to reach convergence to an entropy weak solution.
• Using the numerical diffusion of the scheme, we get a uniform weak BV estimate on the sequence of approximate solutions; this estimate is called "weak" for two reasons: on the one hand it involves the differences and not the differences |u K − u L |, where (K , L) denotes a pair of control volumes sharing a common interface; on the other hand, it only requires that the sum of these differences does not blow up too fast: the difference between the discrete gradient of the interpolated test function and the gradient of the test function itself behaves like the mesh size, and to pass to the limit, we only need that the sum of the differences involving the discrete unknowns be bounded by a term in C /h 1−ε with ε > 0 and C > 0 independent of h.• Using the nonlinear weak convergence and the weak BV estimate, we pass to the limit on the weak form of the discrete entropy and obtain a so-called "process solution" which is an entropy solution up to an integral with respect to the additional parameter α. • Starting from the discrete entropy inequalities that are verified by the scheme, the process solutions are shown to satisfy an entropy inequality.A uniqueness result on the process solutions can then be obtained thanks to a variable doubling technique "à la Kruskov" [17]; this result differs from Di Perna's [16] in that it takes into account the initial condition in the weak entropy formulation, which allows to avoid the more restrictive conditions on the mesh [18].The uniqueness of the process solution entails that a process solution is the unique entropy weak solution.It also yields the (strong) convergence, in L p -spaces, of the approximate solution to the exact solution.
The proof of convergence of the FVM has been obtained for several other problems than the one considered here.However, for systems of PDEs, it is often difficult to obtain compactness results, and LW-consistency then seems an interesting way to make sure that an eventual limit of the scheme is indeed a weak solution of the system.In the next section, we show how this is feasible even on staggered grids, which are often used in the numerical simulation of fluid flows.

LW-consistency and staggered grids
Rectangular staggered grids have been used since the sixties in fluid mechanics (in the wellknown MAC method [19]) including environmental flows [20], see also [21].The mathematical analysis of the MAC scheme has been the object of several recent works, see e.g.[22,23].Systems of partial differential equations for which no existence or uniqueness is known have also been discretized on staggered grids.Let us mention in particular the compressible Euler equations: one of the advantages of using a staggered grid is to produce a scheme that can be mathematically Figure 1.Primal and dual meshes and associated notations for the MAC case.Left: the primal cells; the edges σ and σ belong to F (1) and the edges τ and τ to F (2) .Center: the dual cells associated to F (1) .Right: the dual cells associated to F (2) .
proven to be asymptotically stable to the incompressible limit, as shown for the isentropic case in [24].For such systems, a Lax-Wendroff type theorem is very useful; indeed, although true convergence cannot be proven for lack of compactness properties, such a theorem allows to state that if the scheme converges, and provided some bounds on the approximate solutions are satisfied (these bounds are generally not attainable mathematically but can be verified numerically), the limit of the scheme is an entropy weak solution of the system, see [25,26].The proof of this result may be obtained thanks to the generalization of the Lax-Wendroff theorem to general grids, which include staggered grids such as the MAC grid [27].One of the major additional difficulties for staggered grids is that the discrete unknowns are piecewise constant on different grids.
Consider a rectangular domain Ω ⊂ R 2 (the 3D case can be tackled in the same way), and a possibly non-uniform rectangular grid.We denote by F the set of edges of the mesh, and the internal edge separating the cells K and L is denoted by σ = K |L (see Figure 1).This mesh will be referred to in the following as the primal mesh, and denoted by P .Two dual meshes (three in 3D) are considered, each consisting in a partition of Ω indexed by the vertical and horizontal elements of F, i.e.Ω = ∪ σ∈F (i ) D σ , i = 1, 2, where F (1) (resp.F (2) ) denotes the set of vertical (resp.horizontal) edges.The cells (D σ ) σ∈F are referred to as the dual cells.A half dual cell D K ,σ is half of the rectangle K with side σ (see Figure 1).For an internal edge σ = K |L, the dual cell D σ is the subset of K ∪ L defined as D σ = D K ,σ ∪ D L,σ ; for an external edge σ of a cell K , D σ is the subset D K ,σ of K .
To illustrate the use of the generalized Lax-Wendroff theorem proven in [27], let us consider as a simple example the mass equation of, say, the compressible Euler equation: where ∂ t ρ denotes the time derivative of the density ρ, and div the space divergence.The scalar unknown ρ is associated to the primal cells: The unknowns associated to the i -th component of u are located at the center of the edges of the i -th dual mesh.The associated approximate vector function thus reads: u(x, t ) = (u 1 (x, t ), u 2 (x, t )) t where, for i = 1, 2, Let e (i ) denote the i -th unit vector; the discretization of (20) reads: and, for σ = K |L, ρ n σ stands for a convex combination of ρ n K and ρ n L (for instance the upwind or a MUSCL choice with respect to u n σ ).The initial value for the scalar unknown ρ is defined by We define the space step by h(P ) = max( h(1) , h(2) ), and the time step by δt = max n∈ 0,N −1 (t n+1 − t n ).A sequence of grids is said to be quasi-uniform if the quotients h(1) /h (2) and h(2) /h (1) are bounded by a constant, independent of the grid.
The proof of Lemma 4.1 is given in [7] and is quite simple, especially using the ad hoc tools developed in [7,27].However, the developed arguments may be extended to more complex operators, to deal for instance with the momentum and energy balance equations of the compressible Euler equations, even if the proofs are more tricky.

Conclusion
In this paper, we have presented some concepts for the mathematical analysis of FV schemes, paying special attention to the flux consistency issue and its most direct consequence, i.e. the LW-consistency of the scheme.We emphasize that we gave here a very partial picture of the mathematical world of finite volumes.In many problems (some of them evoked here), the convergence of the discrete solution to a limit follows by compactness arguments in norms strong enough to "feed" the consistency study; in fine, this yields a stronger result than just the scheme consistency, namely the convergence (up to a subsequence if the uniqueness of the solution of the continuous problem is not known) of the numerical solutions to the (a) continuous one (see e.g.[22]).Many parabolic equations enter this framework, including, focusing on fluid flow simulations, incompressible, possibly variable density, or steady barotropic Navier-Stokes equations.In some problems, the continuous solution may be reasonably supposed (or even proven) to be regular, and an error analysis is possible.
The LW-consistency issue is especially important for practical applications in fluid flow simulations.Indeed, in many cases of interest, stronger results are out of reach, and this property is the only one left to mathematically support the design of schemes.This is for instance the case for multi-dimensional flows governed by hyperbolic systems, as shallow-water equations, Euler equations or models for multi-phase flows.For instance, the study evoked in Section 4 is motivated by such a situation: in the last ten years, a class of staggered schemes has been designed for hyperbolic flow problems [25,26,28], and implemented in the open-source software CALIF 3 S developed at IRSN [29]; they are now routinely used for industrial safety applications like hydrogen explosion problems, supposing inviscid or at least vanishing viscosity flows.The accuracy of the numerical schemes involved here is essentially supported by LW-consistency studies [26].