Caterina Calgaro
email: caterina.calgaro@math.univ-lille1.fr

Emile Chane-Kane
email: emile.chane-kane@inria.fr

Emmanuel Creusé
email: emmanuel.creuse@math.univ-lille1.fr

Thierry Goudon
email: thierry.goudon@inria.fr

L ∞ -stability of vertex-based muscl finite volume schemes on unstructured grids; simulation of incompressible flows with high density ratios

Keywords: Finite Volume method, maximum principle property, variable density flows, unstructured meshes

This work is devoted to the design of multi-dimensional finite volume schemes for solving transport equations on unstructured grids. In the framework of MUSCL vertex-based methods we construct numerical fluxes such that the local maximum property is guaranteed under an explicit Courant-Friedrichs-Levy condition. The method can be naturally completed by adaptive local mesh refinements and it turns out that the mesh generation is less constrained than when using the competitive cell-centered methods. We illustrate the effectiveness of the scheme by simulating variable density incompressible viscous flows. Numerical simulations underline the theoretical predictions and succeed in the computation of high density ratio phenomena such as a water bubble falling in air.

Introduction

We are concerned with the numerical simulation of the advection of a scalar quantity: given a velocity field u ∈ R N , depending on the time and space variables (t, x) ∈ R + × R N , we wish to 1 compute the solution ρ(t, x) ∈ R of ∂ t ρ + div x (ρu) = 0.

(1)

To this end, we shall use Finite Volume methods which are quite natural to this problem since they are based on the conservation relationship

d dt C ρ(t, x) dx = - ∂C ρ u • n(x) dσ(x)
which holds for any subdomain C. Therefore, the method is based on a suitable definition of the fluxes on the interfaces of the control volumes C that realize a tessellation of the whole computational domain. Anticipating on precise definitions, we distinguish among these methods between the vertex-based methods where variables are storred at the mesh vertices and cell-centered methods where variables are storred at the centroids of the cells. We refer to [START_REF] Mavriplis | Unstructured mesh discretizations and solvers for computational aerodynamics[END_REF] for a commented overview of the debate on the pros and cons of these methods. In this work, we find advantages in working in the vertex-based framework. The problem we address faces two difficulties:

• Firstly, as it is well known numerical approximation of (1) might generate diffusion and smooth out steep gradients. When dealing with multi-dimensional problems, as pointed in [DL01, DLL07], the origin of numerical diffusion is two-fold. Additionally to the usual diffusion in the direction of the velocity field, which already occurs in dimension one, diffusion phenomena also arise in transverse directions when the mesh is not aligned with the velocity.

• Secondly, we wish to work with unstructured grids. However, removing the Cartesian geometry of the mesh induces a loss of consistency: order one methods do not satisfy order one estimates on unstructured meshes. Fine analysis of this phenomenon can be found in [START_REF] Després | An explicit a priori estimate for a finite volume approximation of linear advection on non-cartesian grids[END_REF][START_REF] Després | Lax theorem and finite volume schemes[END_REF][START_REF] Merlet | Error estimate for finite volume scheme[END_REF][START_REF] Merlet | Error estimate for the finite volume scheme applied to the linear advection equation[END_REF][START_REF] Delarue | Probabilistic analysis of the upwind scheme for transport equations[END_REF].

The necessity to consider unstructured meshes can be motivated as follows. On the one hand, removing structure constraints make easier the tilling of complex geometries. On the other hand, bearing in mind complex and multi-physics problems, we realize that meshing uniformly on the finest scales is definitely non affordable. Then, we wish to couple the scheme with adaptive mesh refinements strategies which make the grids finer close to large gradients. The use of automated mesh generation produces unstructured grids. However, it might be difficult when using such refinements methods to respect the strong geometrical constraints which usually appear in the proof of stability statements. Finally, we address as a requirement for the scheme to preserve the local maximum principle satisfied by the solutions of (1). The motivation in preserving L ∞ estimates and the maximum principle is two fold. First of all the conservation of extrema can be considered as a fundamental physical property of the equation that a numerical scheme must reproduce (a mass or a temperature density should be positive). This is particularly crucial if one thinks of the transport equation as part of a more complex and coupled problem where the violation of the maximum principle might lead to definition troubles. Second of all spurious oscillations can lead, especially in coupled systems where ρ becomes a data for determining the evolution of the velocity, to severe stability issues. We shall illustrate all these aspects with the simulation of variable density incompressible viscous flows, by using the hybrid Finite Volume/Finite Element scheme introduced in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF]. The difficulty of such problems is embodied into two dimensionless parameters: the Reynolds number and the ratio between the densities of the heavy and the light fluids. The larger these numbers, the more difficult the problem. The method we discuss is well-suited to be used with the scheme of [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF] and we shall show that we can indeed consider high density ratios.

Clearly, a simple first order upwind method cannot produce satisfactory results and, anyway, useful engineering calculations generally require second-order spatial accuracy. Hence, a possible avenue consists in increasing the order of approximation by constructing piecewise linear approximations of the unknown. However, high order methods might have a bad behavior when the solution presents strong variations and it has been proved that a multidimensional TVD scheme could be at most of order one [START_REF] Goodman | On the accuracy of stable schemes for 2D conservation laws[END_REF]. In order to preserve the TVD property of the scheme [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF] and so to prevent the solution from non-physical oscillations near the discontinuities of the solution, a limiter (which can be interpreted as a "flux limiter" or a "slope limiter") has to be added. It avoids under and overshoots phenomena and it can lead in particular to a maximum principle on the discrete solution. This is the spirit of Monotone Upstream Scheme for Conservation Law (MUSCL) methods as introduced by Van Leer with the one-dimensional analysis of the resolution of hyperbolic conservation laws [START_REF] Van Leer | Towards the ultimate conservative difference scheme V. A second-order sequel to Godunov's method[END_REF]. Since then, the MUSCL technique has become a standard, used in many academic and industrial codes, see e. g. [START_REF] Mavriplis | Unstructured mesh discretizations and solvers for computational aerodynamics[END_REF]. A large bibliography database is devoted to the developpement of several one-dimensional limiters, and the study of the corresponding schemes convergence [START_REF] Roe | Generalized formulation of TVD Lax-Wendroff schemes[END_REF][START_REF] Sweby | High resolution schemes using flux limiters for hyperbolic conservation laws[END_REF][START_REF] Van Leer | On the relation between the upwind-differencing schemes of Godunov, Engquist-Osher and Roe[END_REF][START_REF] Osher | Convergence of generalized MUSCL schemes[END_REF][START_REF] Yee | Construction of explicit and implicit symmetric TVD schemes and their applications[END_REF]. Unfortunately, the generalization of this approach for multi-dimensional problems is not straightforward at all. Indeed, the way to determine the gradients involved in the reconstruction of the solution at the interfaces of a given control volume, having in mind to satisfy a maximum principle, is not a so easy task. Alternative approaches have been proposed. The Essentially Non Oscillatory scheme, see e.g. [Shu90, CJS + 98], is particularly adapted to Cartesian grids. It uses a large stencil and therefore it is not so convenient in a local mesh refinement context. Subject to recent developments, a promising approach is the Residual Distributed Method as discussed in [START_REF] Abgrall | High order fluctuation schemes on triangular meshes[END_REF]. It is also worth mentioning the work [START_REF] Després | The Vofire method for multicomponent flows on unstructured meshes[END_REF] which adopts a different viewpoint since, contrarily to the MUSCL method, the approximation remains discontinuous over the cells. In [START_REF] Després | The Vofire method for multicomponent flows on unstructured meshes[END_REF] the direction of the velocity is privileged for the construction of the fluxes which furthermore is based on anti-dissipative strategies. The method is shown to be L ∞ stable under a CFL condition comparable to the one that will be discussed below.

The first attempts to develop efficient MUSCL methods on multidimensional unstructured grids were concerned with the design of slope limiting procedures for cell-centered schemes dealing with grids made of triangles by the use of "monoslope" limiters [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF][START_REF] Durlofsky | Triangle based adaptive stencils for the solution of hyperbolic conservation laws[END_REF][START_REF] Liu | A maximum principle satisfying modification of triangle based adaptive stencils for the solution of scalar hyperbolic conservation laws[END_REF][START_REF] Batten | Positively conservative high-resolution convection schemes for unstructured elements[END_REF]. It consists in defining on each control volume a certain gradient, which is used to generate a linear reconstruction of the solution at the interfaces of this control volume, and then to limit it in order to ensure the local maximum principle. Monoslope limiters for vertex-based methods have been also proposed, see e. g. [START_REF] Dervieux | Compressible flow solvers using unstructured grids[END_REF], with applications in aerodynamics, but the most refined limiters strategies have been designed for cell-centered schemes. An elaborate improvement which takes more into account the multidimensional aspects of the problem is due to [START_REF] Hubbard | Multidimensional slope limiters for MUSCL-type finite volume schemes on unstructured grids[END_REF]. Very recently, a "multislope" variant has been proposed first for two-dimensional problems [START_REF] Buffard | Monoslope and multislope MUSCL methods for unstructured meshes[END_REF] and then extended to dimension three [START_REF] Clauzon | Analyse des schémas d'ordre élevé pour les écoulements compressibles. Application à la simulation numérique d'une torche à plasma[END_REF]. In these methods, the numerical gradient depends not only on the control volume under consideration, but also on the interfaces of the control volume. In the same vein, it is also possible to design a scheme which does not necessarily reduce to order one near extrema. In conterparts the scheme does not always satisfy a maximum principle any more. We refer to [START_REF] Christov | New nonoscillatory central schemes on unstructured triangulations for hyperbolic systems of conservation laws[END_REF] which proposes a new reconstruction inspired by ENO. However, the stability investigated in [START_REF] Christov | New nonoscillatory central schemes on unstructured triangulations for hyperbolic systems of conservation laws[END_REF] relies on geometrical constraint related to the finite speed of propagation and the L ∞ stability is not obtained. Nevertheless, this can be enough to obtain stable computations according to the physical problem under consideration and the numerical experiments are indeed very convincing. Note also that proving the L ∞ stability of the multislope cell-centered scheme needs geometrical constraints on the mesh, which can be difficult to verify by automated generation, especially in dimension three, see [START_REF] Ansanay-Alex | Un schéma éléments finis non-conformes/volumes finis pour l'approximation en maillages non structurés des écoulements à faible nombre de Mach[END_REF] where another variant is proposed from [START_REF] Buffard | Monoslope and multislope MUSCL methods for unstructured meshes[END_REF][START_REF] Clauzon | Analyse des schémas d'ordre élevé pour les écoulements compressibles. Application à la simulation numérique d'une torche à plasma[END_REF]. In the present paper, strongly inspired by [START_REF] Clain | L ∞ stability of the MUSCL methods[END_REF], a multidimensional flux limiter strategy is considered for the resolution of (1). However, by contrast to the mentioned references, we deal with vertex-based schemes instead of cell-centered schemes. In this context, it turns out that we are able to propose a gradient reconstruction for which the L ∞ stability can be justified provided a (explicit) Courant-Friedrichs-Levy condition is fulfilled. The proof works for divergence free velocity field, not necessarily constant. Furthermore, the stability is obtained without the restrictive geometrical requirements needed for the cell-centered methods. Numerical tests show the robustness and the accuracy of the numerical method we derive. In particular this treatment of the mass conservation equation allows to consider high density ratio when simulating multicomponent incompressible viscous flows.

The schedule of the paper is the following. In section 1, we explain our motivation that comes from the simulation of incompressible Navier-Stokes equations with inhomogeneous density. We briefly recall the governing equations as well as the numerical scheme introduced in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF]. Then, we detail the vertex-based Finite Volume scheme we propose for solving (1), with a multislope gradient reconstruction. In section 2, we establish the maximum principle property for the scheme under an explicitely given CFL condition. Finally, section 3 is devoted to numerical results and we make comparisons with other simulations in recent bibliography. Three classical benchmark test cases are addressed: the translated field, the rotational field and the falling droplet.

1 Motivation and description of the numerical scheme 1.1 Density dependent Navier-Stokes equations and time discretization Our motivation comes from the simulation of the variable density incompressible Navier-Stokes system. From now on, Ω stands for an open bounded polygonal subset on R 2 . We are interested in the following PDEs system

∂ t ρ + div x (ρu) = 0, (2) ∂ t (ρu) + Div x (ρu ⊗ u) + ∇ x p -µ∆ x u = f , (3) div x u = 0. ( 4
)
Here, ρ(t, x) ≥ 0 represents the density, p(t, x) ∈ R the pressure and u(t, x) ∈ R 2 the velocity field of the fluid. The description of the external force is embodied into the right hand side f (t, x) of (3) and µ > 0 stands for the (dynamic) viscosity. The unknowns depend on time t ≥ 0 and position x ∈ Ω ⊂ R 2 . Given vector fields u and v we set div x (u) = 2 i=1 ∂ x i u i , and u ⊗ v is the 2 × 2 matrix with components u i v j ; given a matrix valued function A we denote Div x A the vector having components 2 j=1 ∂ x j A ij . There are several difficulties for solving numerically the system. To start with, the divergence free constraint and the underlying definition of the pressure have to be considered carefully. In particular, it yields some constraints between the spaces of approximation for the velocity and the pressure (the so-called "inf-sup condition") and many Finite Elements methods have been designed for the treatment of the space homogeneous density incompressible Navier-Stokes equations. However, dealing with inhomogeneous densities leads to new numerical issues in order to treat the coupling with accuracy and the standard Finite Elements methods cannot be applied directly. The coupling also yields specific stability and conditioning questions. Classical illustrations consist in simulating the formation of Rayleigh-Taylor instabilities and the motion of droplets or bubbles in a fluid. Clearly, increasing the ratio of the extreme densities --that can be associated to the socalled Atwood number --makes the problem more challenging for numerics. We refer for instance to [Boy02, GQ00, CCG08] for a further description of these difficulties and discussion of techniques to treat the problem. In [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF], we have introduced an original hybrid Finite Volume/Finite Element scheme to compute the solution of (2)-( 3)-(4).

Let us denote ∆t the time step and t n = n ∆t, n ≥ 0. Let us assume that the numerical solution at time t n , namely (ρ n , u n , p n ), is known on the computational domain. The numerical scheme is based on a time splitting of the system (2)-(3)-(4) known as the "Strang splitting" [START_REF] Strang | On the construction and comparison of difference schemes[END_REF] :

1. The new density field, ρ n+ 1 2 , is computed by solving on the time interval (n∆t, (n + 1 2 )∆t) the transport equation

∂ t ρ n+ 1 2 + div x (ρ n+ 1 2 u n ) = 0, (5)
with suitable boundary conditions on ρ n+ 1 2 , using a vertex-based finite volume method.

2. The new velocity and pressure fields, u n+1 and p n+1 , are computed by the resolution on the time interval (n∆t, (n + 1)∆t) of the system

ρ n+ 1 2 ∂ t u n+1 + u n+1 • ∇ x u n+1 + ∇ x p n+1 -µ∆ x u n+1 = f n+1 , (6) div x u n+1 = 0, (7)
completed by the specification of boundary conditions on u n+1 , using a finite element method.

3. Finally, the new density field, ρ n+1 , is computed by solving on the time interval ((n+ 1 2 )∆t, (n+ 1)∆t) the transport equation

∂ t ρ n+1 + div x (ρ n+1 u n+1 ) = 0, (8)
with suitable boundary conditions on ρ n+1 , using once again the vertex-based finite volume method.

Then we go back to the first step (using n + 1 instead of n) to compute the solution at the following time step. A key point of the method introduced in [CCG08] is the "compatibility relation", which has to be carefully fulfilled so that the discrete divergence free constraint is ensured both for the finite element and the finite volume interpretations. We will not come back on this point here neither on the description of the global numerical scheme. We only consider that we have to solve (5) and ( 8), having at hand a given discrete velocity field u. In particular, in view of the application we have in mind, it turns out that

• it is natural to keep track of the interfaces and to appeal to mesh refinements strategies so that a fine grid is used only in regions of steep density gradients and a coarse grid elsewhere.

• as the density ratio increases (considering air and water this ratio is of order of 1000) the simulation becomes highly sensitive to violation of the maximum principle which might lead to severe numerical instabilities for the whole system.

Space discretization and notations

Let ∂Ω be the Lipschitz boundary of Ω, and n(x) the outer unit normal vector at x ∈ ∂Ω. Let T h be a partition of Ω (named primal mesh) composed of conforming and isotropic triangles

T k , with k ∈ [1, K].
For each element T k of the mesh, we denote B k the barycenter of the triangle and |T k | the area of T k . Moreover, we denote by A i , i ∈ [1, I], the vertices of T h . For a given node A i , let V(i) be the set of the indices of the neighboring nodes A j which are vertices of the elements T k sharing the common point A i . For a given triangle T i , let W(i) be the set of the indices of the neighboring triangles T j sharing at least one node with T i . We also denote by A ij the middle point of the edge [A i A j ]. For each vertex A i , we associate a polygon denoted by C i , which defines a second partition of Ω (named dual mesh). Here, this polygon C i is defined by joining the barycenters B k of the neighboring elements T k which share the node A i with the middle points A ij of the edges [A i A j ], for j ∈ V(i). This construction is referred to as the "CV1" control volume (see Figure 1), and we will see later that a variant can be proposed. We denote

Γ - ij (resp. Γ + ij ) the segment [B j-1 A ij ] (resp. [A ij B j ]) (here we visit the neighboring triangles around A i in the counterclockwise direction), A - ij (resp. A + ij ) the middle of Γ - ij (resp. Γ + ij ) and n - ij (resp. n + ij ) the unit outward normal to C i along Γ - ij (resp. Γ + ij ). Let |C i | stand for the area of C i and |Γ - ij | (resp. |Γ + ij |) for the length of Γ - ij (resp. Γ + ij ).
Then, the boundary of C i splits into several segments Γ - ij and Γ + ij and we write

∂C i = j∈V(i) (Γ - ij ∪ Γ + ij ).
We introduce now the following quantity C T h = max

A i A j A j+1 A j-1 M + ij M + ji N + ij A k + ij A k + ij +1 B j B j-1 T j T j-1 A ij n - ij n + ij A - ij A + ij Γ - ij Γ + ij
T j ∈T h max k∈W(i) h T j ρ T k , (9)
where h T j is the radius of the circombscribed circle to T j and ρ T k the smallest height in T k . The constant C T h measures the mesh size variation from an element to a neighboring one. It characterizes the mesh regularity, even if it does not exactly correspond to the usual elementwise aspect ratio in the Ciarlet sense [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF].

Finite volume scheme

As said above, let us assume that ρ n , the numerical solution at time t n , is known on the computational domain. The new density ρ n+1 is computed by solving on the time interval (n∆t, (n + 1)∆t) the transport equation

∂ t ρ n+1 + div x (ρ n+1 u) = 0. ( 10
)
Here and below, the velocity field is a given function u : R × R 2 → R 2 which fulfills the divergence free condition. We will make more precise this requirement at the discrete level later on. The computation of the solution is performed using a usual vertex-based Finite Volume (FV) scheme which provides the approximated value ρ n+1 i , i ∈ [1, I]:

ρ n+1 i ≈ 1 |C i | C i ρ(t n+1 , x) dx.
Let us remark that the vertex-based choice of the FV method (instead of the cell-centered one) allows to write a P 1 Lagrange interpolation of the density if needed, by using this constant value ρ n+1 i on C i as the value of ρ n+1 at node A i (see [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF]). Then, ρ n+1 i , i ∈ [1, I] is computed by using the formula

ρ n+1 i = ρ n i - ∆t |C i | j∈V(i) |Γ + ij |G + ij ((ρ + ij ) n , (ρ + ji ) n ) + |Γ - ij |G - ij ((ρ - ij ) n , (ρ - ji ) n ) , (11)
where the right hand side of (11) is evaluated using an upstream flux :

G ± ij (ρ 1 , ρ 2 ) = ρ 1 u ± ij • n ± ij if u ± ij • n ± ij ≥ 0, ρ 2 u ± ij • n ± ij if u ± ij • n ± ij < 0. ( 12
)
With this choice, the monotone property of the numerical flux is ensured, namely :

∂G ± ij ∂ρ 1 (ρ 1 , ρ 2 ) ≥ 0 and ∂G ± ij ∂ρ 2 (ρ 1 , ρ 2 ) ≤ 0. (13)
In (12), u + ij (resp. u - ij ) stands for the value of the velocity field at the point A + ij (resp. A - ij ). In fact, bearing in mind multi-physics coupled models, the velocity field is also known on a discrete way, and u + ij is rather an approximation of the continuous velocity associated to the point A + ij . (Coming back to [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF], it is defined from the finite element discretization of the momentum equation of the Navier-Stokes system.) Therefore, the discrete velocity field is defined so that it preserves the continuous divergence free property at the discrete level: for any given constant density ρ c on the whole domain, and whatever i ∈ [1, I],

j∈V(i) |Γ + ij |G + ij (ρ c , ρ c ) + |Γ - ij |G - ij (ρ c , ρ c ) = 0 (14)
holds. The consistency property ( 14) is a key point of the further analysis. A first order accuracy in space is obtained by choosing, as interface reconstruction involved in (11), the value (ρ

+ ij ) n = (ρ - ij ) n = ρ n i for all i ∈ [1, I] and j ∈ V(i).
It is well-established that in this case the L ∞ stability is satisfied under an appropriate CFL condition. The second order accuracy in space can be reached by using a MUSCL technique [VL79, VL84], with a reconstruction operator which gives new values ρ ± ij and ρ ± ji on both sides of Γ ± ij using the values ρ n i , ρ n j and ρ n k for some k ∈ V(i) ∪ V(j). In this work, we aim at developing such a reconstruction strategy with the goal of preserving the local maximum property :

∀i ∈ [1, I], min j∈V(i) (ρ n i , ρ n j ) ≤ ρ n+1 i ≤ max j∈V(i) (ρ n i , ρ n j ), ( 15
)
in a vertex-based finite-volume context. For further considerations it is convenient to introduce the quantity :

u i,∞ = max j∈V(i) max( u + ij l 2 (R 2 ) , u - ij l 2 (R 2 ) ). ( 16
)
1.4 Second order reconstruction 1.4.1 Monoslope procedure

A classical approach to reach the second order accuracy is based on a piecewise linear approximation of the density defined using a gradient vector of the unknown. The upstream or downstream gradient vectors can so be defined on each element T k or on each node A i of T h . Usually, (∇ρ

) | T k , k ∈ [1, K],
is the gradient vector computed from the three values of ρ n in the triangle T k , ∇ρ i = (

k∈Z(i) |T k |(∇ρ) | T k )/( k∈Z(i) |T k |), i ∈ [1, I],
with Z(i) the set of the indices of the triangles T k sharing ]. Using this approach on unstructured meshes, by considering a one-dimensional problem independently for each interface reconstruction, it is easy to see that the maximum principle will not be respected. To illustrate it, let us consider the following example where a density field is defined on the mesh given on figure 2 : ρ n is the density value at node A and ρ n i at node A i , 1 ≤ i ≤ 6 (see full circles). The density is identically equal to 1 (ρ n = ρ n 2 = ρ n 3 = ρ n 4 = ρ n 5 = ρ n 6 = 1), except at node A 1 where ρ n 1 = 1+ε, ε > 0. We are interested in the new value of the density at node A after one iteration in time, namely ρ n+1 , using (39) (see below) which is a variant of (11) in the particular case of such a structured mesh. The velocity u is supposed to be equal to u = (1, 0) T . We use the one-dimensional reconstruction with the β-scheme, β = 1 3 , and the Van-Leer limiter. Thus it is necessary to reconstruct the value of the density on each part of the interface of the control volume (see empty circles). Using one of the usual way to proceed [START_REF] Bruneau | A finite volume method with efficient limiters for solving conservation laws[END_REF], the reconstructed value ρn 1 of the density at interface Γ 1 located between A and A 1 is computed by :

A i [FLLO89, BR97
ρn 1 = ρ n + 1 2 βα(r 1 )∆ρ 1 + (1 -α(r 1 ))(1 -β)∆ρ 1 , where β = 1 3 , r 1 = β∆ρ 1 (1 -β)∆ρ 1 , α is the Van-Leer flux-limiter function (α(r) = r + |r| 1 + r ) , ∆ρ 1
is the "upstream" variation of the density, computed here using an averaged value of the density

A A 2 A 3 A 4 A 5 A 6 A 1 Γ 1 Γ 2 Γ 3 Γ 4 Γ 5 Γ 6 h h Figure 2: A particular mesh.
gradients in triangles surrounding the node A, and ∆ρ 1 is the "downstream" variation of the density, computed with the values of the density only at nodes A and A 1 . Using such a strategy, we naturally arrive at :

ρn 1 = 1 + 16 63 ε; ρn 2 = ρn 3 = ρn 4 = ρn 5 = ρn 6 = 1.
Consequently, a simple calculation using (39) and geometrical considerations of this particular control volume lead to :

h 2 ∆t (ρ n+1 -ρ n ) = - 32h 189 ε < 0,
and neither the maximum principle nor the local extremum diminushing (LED) properties of the scheme are fulfilled since ρ n+1 < 1. Such an approach can nevertheless ensure the maximum principle property in particular cases corresponding to Cartesian meshes on which a tensorial analysis can be performed, coming back under a given CFL to the one dimensional framework.

A first remedy is so to consider a monoslope method, as introduced in [Hub99, BC08] for the cell-centered framework. The idea of [START_REF] Hubbard | Multidimensional slope limiters for MUSCL-type finite volume schemes on unstructured grids[END_REF] can be adapted to the vertex-based context. It would consist in the computation of a gradient operator, the vector ∇ρ i for example, and to use only this gradient on the volume C i in order to obtain all the reconstructed values (ρ

+ ij ) n (resp. (ρ - ij ) n ) on the segments Γ + ij (resp. Γ - ij ), j ∈ V(i).
A limiting strategy would then be applied in order to avoid the creation of new extrema during the reconstruction process: it is indeed in this approach a necessary condition to ensure the maximum principle property (see the appendix of [START_REF] Batten | Positively conservative high-resolution convection schemes for unstructured elements[END_REF] for the proof). So we would have :

ρ ± ij = ρ i + α i ∇ρ i • A i A ± ij for any j ∈ V(i),
where α i ∈ [0, 1] has to be determined not to create any new extremum at the interface. Hence, the limiting procedure can reduce drastically the scheme locally to the order one in all the space directions around A i , while it is only needed for one of these directions. It constitutes a severe drawback, even for very regular solutions. Nevertheless, if this methodology is used, a sufficient condition to obtain the maximum principle remains that the linear approximation satisfies the following conservation property:

C i (ρ i + α i ∇ρ i • A i X) dX = |C i | ρ i ∀i ∈ [1, I].
It can be immediately seen that if a cell-centered FV approach automatically satisfies the conservation property, because of the triangular geometry of the control volume, it is not the case with the vertex-based FV approach for which the control volume is a given polygonal. This is the reason why a multislope multidimensional procedure is necessary in the vertex-based context, as confirmed by numerical experiments.

Multislope procedure

We need to intoduce a few notation. Given two points D and E in the plane, we denote (DE) the line passing through D and E, and (DE[ is the half-line starting from E and containing D (but not E). Moreover, we define

D i = ∪ A i ∈T k T k . Obviously, C i ⊂ D i .
This subsection is devoted to the description of the density reconstruction (ρ + ij ) n on the interface Γ + ij by a multislope strategy. Of course, a similar reconstruction will be derived for the value of (ρ - ij ) n , as well as for the values of (ρ + ji ) n and (ρ - ji ) n . In order to simplify the notations, the time indices are dropped; ρ + ij will so have to be understood as (ρ + ij ) n . Finally, the node A i is supposed not to belong to the boundary of Ω (in other words, it is an internal node). We introduce the point

M + ij = (A i A + ij ) ∩ [A j A j+1 ].
Definition 1 Let us define α + i,j ∈ [0, 1] such that :

A i M + ij = α + i,j A i A j + (1 -α + i,j ) A i A j+1 (17)
This definition, associated to a P 1 -Lagrange-piecewise interpretation of the density field on the primal mesh, allows to define a density value at nodes M + ij by setting :

ρ M + ij = ρ i + α + i,j (ρ j -ρ i ) + (1 -α + i,j )(ρ j+1 -ρ i ).
Now, we define the point

N + ij = (A i M + ij [ ∩ ∂D i .
Because of the definition of the monotone property of the flux (13), we clearly have A ± ij ≥ 0 and B ± ij ≥ 0.

On the one hand, lemma 2 allows to write :

∆ρ + ij = k∈V(i) w + ijk (ρ i -ρ k ) with w + ijk ≥ 0 and k∈V(i) w + ijk ≤ 7τ 12 C T h ,
and similarly :

∆ρ - ij = k∈V(i) w - ijk (ρ i -ρ k ) with w - ijk ≥ 0 and k∈V(i) w - ijk ≤ 7τ 12 C T h .
On the other hand, thanks to lemma 1, the convex property of the reconstruction leads to :

∆ρ + ij = θ + ji (ρ M + ji -ρ j ) + (ρ j -ρ i ) = θ + ji (ρ M + ji -ρ i ) + (1 -θ + ji )(ρ j -ρ i ). Since M + ji ∈ [A i A j+1 ], ρ M +
ji is defined by a linear interpolation between ρ i and ρ j+1 . So we get :

∆ρ + ij = θ + ji A i M + ji A i A j+1 (ρ j+1 -ρ i ) + (1 -θ + ji )(ρ j -ρ i ). ( 29
)
By identification, we can find some coefficients w+ ijk such that :

∆ρ + ij = k∈V(i)
w+ ijk (ρ k -ρ i ) with w+ ijk ≥ 0 and

k∈V(i) w+ ijk ≤ 2.
Similarly, we find wijk satisfying

∆ρ - ij = k∈V(i) w- ijk (ρ k -ρ i ) with w- ijk ≥ 0 and k∈V(i) w- ijk ≤ 2.
The numerical scheme (28) can now be recast as follows

ρ n+1 i = ρ i -∆t k,j∈V(i) (A + ij w + ijk (ρ i -ρ k ) -B + ij w+ ijk (ρ k -ρ i )) -∆t k,j∈V(i) (A - ij w - ijk (ρ i -ρ k ) -B - ij w- ijk (ρ k -ρ i )), (30)
which can be written as :

ρ n+1 i = c ii ρ i + k∈V(i) c ik ρ k , (31)
with :

           c ii = 1 -∆t k,j∈V(i) A + ij w + ijk + B + ij w+ ijk + A - ij w - ijk + B - ij w- ijk , c ik = ∆t j∈V(i) A + ij w + ijk + B + ij w+ ijk + A - ij w - ijk + B - ij w- ijk .
We obviously have

c ii + k∈V(i) c ik = 1. ( 32
)
Moreover, if ∆t is now chosen such that whatever

i ∈ [1, I], ∆t ≤   k,j∈V(i) A + ij w + ijk + B + ij w+ ijk + A - ij w - ijk + B - ij w- ijk   -1 , (33)
then we get :

0 ≤ c ii ≤ 1 and 0 ≤ c ik ≤ 1. ( 34
)
The scheme formulation (31) together with relations (32) and (34) on the coefficients leads to the local maximum preserving property, since ρ n+1 i is written as a convex combination of ρ i and ρ j , j ∈ V(i). Now, in order to recover (26), we just have to combine the fact that

A ± ij ≤ |Γ ± ij | |C i | u i,∞ and B ± ij ≤ |Γ ± ij | |C i | u i,∞
, associated to previous properties on the coefficients w ± ijk and w± ijk .

Remark 3 In the particular case of a constant velocity u associated to a uniform mesh like the one on figure 2 with hexagonal control volumes (C T h = 2) and with the use of the minmod limiter (τ = 1), the CFL condition (26) leads to :

∆t ≤ 9h 19( √ 2 + 2 √ 5) u ≈ h 12.4 u ( 35
)
At first sight, (35) can appear as a very severe CFL restriction, but it is very similar to other CFL-conditions on unstructured meshes used to ensure the maximum principle, see for instance [START_REF] Després | The Vofire method for multicomponent flows on unstructured meshes[END_REF][START_REF] Christov | New nonoscillatory central schemes on unstructured triangulations for hyperbolic systems of conservation laws[END_REF]. Moreover, we must bear in mind that all the proofs were established by considering sufficient conditions which always correspond to the worst case we could encounter. Consequently, it is observed in numerical computations that very good results are obtained even with less restrictive CFL numbers. Instead of considering the dual mesh made of CV1 control volumes displayed on figure 1, we could build another dual mesh using more simple ones, by joining the barycenters of all the triangles surronding A i leading to "CV2" control volumes (see figure 3). In that case, we define Γ ij = [B j-1 B j ] as the interface between the nodes A i and A j , and we so have to reconstruct the density value at the node

A variant of the numerical scheme

A i A j A j+1 A j-1 N ij A k ij +1 B j B j-1 T j T j-1 Q ij n ij A k ij Γ ij
Q ij = [B j-1 B j ] ∩ [A i A j ]. In general, Q ij is no more the middle of [A i A j ],
except in the case of very particular structured meshes. There are clearly several advantages with this construction. In particular, the downstream slope p down ij is now directly available without any linear interpolation process :

p down ij = ρ j -ρ i A i A j ,
even if the upstream gradient definition p up ij always needs the value of the density at node N ij = (A i A j [ ∩ ∂D i by a linear interpolation reconstruction :

p up ij = ρ i -ρ N ij A i N ij .
If we now introduce the mesh parameter α = min 1≤i≤I, j∈V(i)

A i A j A i Q ij ,
we have 1 ≤ α ≤ 2, with the property that the more regular the mesh, the more α close to 2. Reproducing the arguments used with the CV1 control volume, and assuming that the limiter verifies ∀r ∈ R + * , Ψ(r) ≤ min(αr, τ )

instead of (23), the scheme can be proved to satisfy the local maximum property under the CFL :

∆t ≤ min 1≤i≤I |C i | u i,∞ 1 + τ C T h α j∈V(i) |Γ ij | . (37)
Clearly ( 37) is less restrictive than (26).

Remark 4 In the particular case mentionned in remark 3, the CFL becomes

∆t ≤ 3h 4( √ 2 + 2 √ 5) u ≈ h 7.8 u . (38)
However there are also some drawbacks with the CV2 construction. Firstly, the node Q ij does not stand in the middle of [B j-1 B j ]. Consequently, the new density computation :

ρ n+1 i = ρ n i - ∆t |C i | j∈V(i) |Γ ij |G ij ((ρ ij ) n , (ρ ji ) n ), (39)
where the right hand side of (39) is as previously evaluated using an upstream flux :

G ij (ρ 1 , ρ 2 ) = ρ 1 u ij • n ij if u ij • n ij ≥ 0, ρ 2 u ij • n ij if u ij • n ij < 0, (40)
does not correspond to a second-order accuracy numerical integration scheme. Secondly, we have to define the quantity u ij , corresponding to the discrete value of the velocity field at Q ij from the previous values of u + ij and u - ij , corresponding respectively to the discrete values at nodes A + ij and A - ij (see figure 1), under the consistency constraint

j∈V(i) |Γ ij |G ij (ρ c , ρ c ) = 0. ( 41
)
It leads to further computational costs, as well as a loss of accuracy in the discrete velocity field which can affect the global scheme accuracy.

In section 3, the construction CV1 and CV2 will be tested to analyze and compare their behaviors, as well to confirm the theoretical considerations.

Numerical Results

In this section we present several numerical tests. We start by dealing with given velocity fields in order to discuss the quality of the approximation. Then, we consider the incompressible Navier-Stokes system and we show the ability of the scheme to deal with large density variation. The simulations are performed on unstructured meshes, possibly with refinement strategies. Throughout this section the tilling of the domain is made by using usual conforming and isotropic unstructured meshes composed of triangles, associated to either CV1 or CV2 control volumes. In particular we avoid with too small angles and all the geometric constants involved in the stability analysis remain controlled. For each mesh considered, the time step is fixed at ∆t = 0.2 h min / u ∞ , with h min the size of the smallest edge of the mesh (see in Table 1 the values of ∆t when u ∞ = 1). We first note that for all the following tests, this choice of CFL condition is sufficiently small to Mesh level numerically ensure the maximum principle, even if it is far less restrictive than the one expected by the theory. Indeed, it is possible to chose a larger ∆t than the one given by (26) or (37) according to the control volumes type CV1 or CV2 used, while preserving the maximum principle.

h

Translated field

The first numerical test consists in a convected profile on Ω =]-1, 1[ 2 , with a constant velocity field u = (λ, λ) T , where λ = 0.5. For the initial condition ρ(t = 0, x) = ρ 0 (x), two different functions are used. First, a regular function :

ρ 0 = 0.5(1 + cos(4πr)) if r ≤ 0.25, 0 otherwise,
and then a discontinuous one :

ρ 0 = 1 if r ≤ 0.25, 0 otherwise.
Setting r = (x 1 + 0.25) 2 + (x 2 + 0.25) 2 , the initial solutions are centered at the point (-0.25, -0.25).

The computations are performed until the final time t = 1.

Tables 2 and 3 list the order of the errors in the L 1 (Ω) norm and in the L ∞ (Ω) norm between the exact solution and the approximation at the final time t = 1, computed respectively by :

A i ∈T h |C i ||ρ i -ρ exact (A i )| and max A i ∈T h |ρ i -ρ exact (A i )|.
When a regular function is used, the L 1 (Ω) rate of convergence is very good (namely O(h 1.77 ) for CV1 and CV2 control volumes between the two most refined meshes), and seems to be better than in the cell-centered case [START_REF] Clauzon | Analyse des schémas d'ordre élevé pour les écoulements compressibles. Application à la simulation numérique d'une torche à plasma[END_REF]. Otherwise, a lower rate of convergence is achieved when a discontinuous function translates, leading to similar results than those obtained in [START_REF] Buffard | Monoslope and multislope MUSCL methods for unstructured meshes[END_REF][START_REF] Clauzon | Analyse des schémas d'ordre élevé pour les écoulements compressibles. Application à la simulation numérique d'une torche à plasma[END_REF].

Mesh level 1 to 2 2 to 3 3 to 4 4 to 5 L 1 in CV1 Qualitatively, 10 equidistributed isovalues of the density between ρ = 0.01 and ρ = 0.8 are displayed at initial and final times, for the regular density (Figure 4) as well as for the discontinuous one (Figure 5) on the mesh corresponding to mesh level 5 (see table 1). Since control volumes type CV1 and CV2 give very similar results, only results obtained with CV1 volumes are displayed. For the regular solution, the numerical diffusion is quite small, as expected with a second-order accuracy scheme. For the discontinuous one, the radial symmetry is well preserved.

Finally, we consider a structured mesh like the one considered in remarks 3 and 4, composed of 256 × 256 nodes. The initialisation is made with the discontinuous solution, and we are looking at the solution obtained after only one iteration in time, according to the time-step value ∆t which is no more equal to 0.2h min / u ∞ . We are particulary interested in the largest value of ∆t for which the maximum principle is still ensured. With the CV1 control volumes, this value is equal to 7.80E-03, compared to 8.3E-03 for the CV2 control volumes. On the one hand, it can be observed that these values are significantly larger than the one given by the theory (respectively 3.1E-04 from (35) for CV1 and 5.0 E-04 from (38) for CV2). On the other hand, it also shows that the use of the CV2 volume control leads to a slightly less restrictive CFL condition, as expected by the theory.

Rotational field

The second numerical test consists in a convected profile on Ω =] -1, 1[ 2 , with a stationary rotating velocity field u = (-x 2 , x 1 ) T . For the initial condition ρ 0 (x), the previous regular and discontinuous functions are used. In this case, the continuous as well as the discrete velocity fields verify the divergence free constraint div x u = 0 and, consequently, the numerical scheme uses the consistency property (14) or (41) according to the finite control volume used. As above, tables 4 and 5 list the order of the errors in the L 1 (Ω) norm (eventually in the L ∞ (Ω) norm too) between the exact solution and the approximation at the final time t = 1.50.

We observe that the rate of convergence is nearly the same using control volumes CV1 or CV2 for the regular function (around O(h 1.80 ) in the L 1 (Ω) norm between the two most refined meshes). Moreover, these convergence rates are as good as in the constant velocity field case. Concerning the discontinuous function, results are once again similar for both control volumes.

Mesh level 1 to 2 2 to 3 3 to 4 4 to 5 L 1 in CV1 The last test about the rotational field is devoted to the use of the local refinement mesh strategy. It takes advantage of the ability of the scheme in maintaining the local maximum property on fully unstructured meshes. This isotropic mesh-refinement process is based on the BAMG software [BAM], using the density value to detect the high gradient zones and to generate the new mesh at each time interval (∆t) mesh = 15∆t. Figure 8 displays the mesh as well as the isovalues of the density for the discontinuous function at several times of the simulation. Here the final time used is t = 2π. We obtain the same L 1 (Ω) error on the density than the one obtained in the non locally refined case by using only 8 000 points in the mesh compared to 77 000 points.

Falling Droplet

The goal of this test is to investigate the effectiveness of the proposed scheme by simulating variable density incompressible flows, in the context of very high density ratio such as a water bubble falling in air. This benchmark is already described in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF], but only with a moderate density ratio (namely 100) and a tensorial Cartesian mesh, for which the maximum principle property is simply ensured by the use of a one-dimensional limiter. Here, a higher density ratio is considered. An unstructured mesh coupled to a local mesh refinement strategy allows far less degrees of freedom for a better accuracy. In such a case, the use of the limiting process described above for the vertexbased finite volume scheme is absolutely essential to provide physical as well as stable simulations. A heavy "droplet" falls through a light fluid and impacts the flat surface of the heavy fluid in a cavity. The computational domain is (0, d) × (0, 2d), where d = 1 and at t = 0 the fluid is at rest with density:

ρ(x, y) = ρ M = 1000 if 0 ≤ y ≤ 1 or 0 ≤ r ≤ 0.2, ρ m = 1 if 1 < y ≤ 2 or 0.2 < r,
where r = (x -0.5) 2 + (y -1.75) 2 . As in [START_REF] Schneider | Extension of finite volume compressible flow solvers to multidimensional, variable density zero Mach number flows[END_REF], the equations are made dimensionless by using the following references: ρ m for density, d for length, d/G for time, so that the reference velocity is √ d G. In the dimensionless equations, the gravity term is f = (0, -ρ) and the Reynolds number is defined by Re = ρ m d 3/2 G 1/2 µ . In our test, the viscosity of the fluid is supposed to be constant in the whole domain and we have Re = 3132. The difficulty obviously comes from the very strong discontinuity at the interface. At initial time, the mesh is adapted to the initial condition using the BAMG mesh-refinement procedure [BAM], associated to the density gradient in order to detect the zones that need to be locally refined. Then, during the simulation, the time frequency of the remeshing process mainly depends on the quantity ||u|| ∞ h min , where ||u|| ∞ is the L ∞ (Ω) norm of the discrete velocity field, and h min the length of the smallest edge in the mesh. At each remeshing step, once the new mesh is obtained, the old velocity field u old is not only linearly interpolated on the new mesh, leading to a new discrete field u * , but also projected on the discrete divergence-free space in order to obtain the velocity field u new on the new mesh. This is a crucial point to ensure (14), and so to make the finite volume scheme maximum principle preserving. Numerically, it consists in the minimisation of ||u *u new || L 2 (Ω) under the constraint that u new is divergence-free from the finite element point of view. It leads to the resolution of a classical saddle-point problem.

For this computation, the smallest edge in the mesh is h min = 1/400, so that the number of nodes in the mesh evolves from 15000 to 21000 during the simulation. Several density isovalues are displayed on Figure 9 and the corresponding meshes on Figure 10. The maximum principle is perfectly ensured, and results are very consistent with those of [Boy02, PAB + 97, SBGK99] in slightly different contexts. Note also that these references make use of regular Cartesian grids.

After the splash of the droplet, some areas of lighter density contours appear within the heavy fluid, corresponding to a phenomenon for which some particles of lighter fluid are trapped in the heavy fluid after the splash. The change of topology can be clearly distinguished since the two independent regions are merging. The formation of water waves resulting from the splash moving toward the walls can also be remarked. Let us note that a difficulty, compared to smaller density ratios, lies in the bad condition number of the linear system (6)-(7) in the finite element step of the algorithm. This feature has to be taken into account in order to actually ensure a discrete divergence-free velocity filed, which is a fundamental asumption to provide the stability of the finite volume scheme. To obtain this satisfactory convergence process, the value of h min has so to be chosen sufficiently small according to the value of ρ M , and a suitable preconditionning strategy has to be used (we refer to [START_REF] Calgaro | Incremental incomplete LU factorizations with applications to time-dependent PDEs[END_REF] for progress in this direction).

Conclusion

In this paper, a L ∞ -stable MUSCL vertex-based finite volume scheme is developped using a multidimensional slope limiter. On the one hand, it is proved to preserve a local maximum principle property under a certain Courant-Friedrichs-Levy condition. By contrast to the cell-center finite volume case, the L ∞ stability can be obtained without any geometry restriction on the mesh. On the other hand, the proof of the maximum principle is obtained for any discrete divergence-free velocity field, generalizing the constant velocity field case. According to the chosen control volume, two variants are derived. Numerical tests on several classical benchmarks lead to very satisfactory convergence rates. The proposed Finite Volume scheme is incorporated in a more complete code to simulate the evolution of variable density flows with very high density ratios on unstructured meshes, combined to a local mesh refinement strategy.

Figure 1 :

 1 Figure 1: Construction of the dual mesh : control volume CV1 around the node A i .

Figure 3 :

 3 Figure 3: Construction of the dual mesh : control volume CV2 around the node A i .

Figure 4 :

 4 Figure 4: Density contours of a regular function at t = 0 (left) and t = 1 (right).

Figure 5 :

 5 Figure 5: Density contours of a discontinuous function at t = 0 (left) and t = 1 (right).

  of a discontinuous function: order of the errors for various control cell and meshes, t = 1.50.Once again, 10 equidistributed isovalues of the density between ρ = 0.01 and ρ = 0.8 are displayed at initial and final times in Figure6(resp. 7) when a regular density (resp. a discontinuous density) is convected by the rotating velocity field. The views of both control volumes CV1 and CV2 are very similar, then only the first one is displayed.

Figure 6 :

 6 Figure 6: Density contours of a regular function at t = 0 (left) and t = 1.50 (right).

Figure 7 :

 7 Figure 7: Density contours of a discontinuous function at t = 0 (left) and t = 1.50 (right).

Figure 8 :

 8 Figure 8: Mesh (top) and density contours (bottom) of a discontinuous function at t = 0, t = π/2, t = π and t = 2π with refinement strategy, CV1 control volumes.

  (a) t = 0.090 (b) t = 0.240 (c) t = 0.505 (d) t = 0.745 (e) t = 0.995 (f) t = 1.125 (g) t = 1.250 (h) t = 1.361

Figure 9 :

 9 Figure 9: ρ M = 1000, h min = 1/400, density contours = [200, 350, 400, 450, 500, 550, 600, 650, 800].

  Figure 10: ρ M = 1000, h min = 1/400, meshes.

Table 1 :

 1 Time step values with u ∞ = 1

		min	∆t
	1	4.42E-02 8.84E-03
	2	2.21E-02 4.42E-03
	3	1.10E-02 2.21E-03
	4	5.50E-03 1.10E-03
	5	2.76E-03 5.52E-04

Table 2 :

 2 Translation of a regular function: order of the errors for various control cell and meshes.

		1.62	1.52	1.69	1.77
	L 1 in CV2	1.65	1.60	1.76	1.77
	L ∞ in CV1 0.848	1.17	1.20	1.21
	L ∞ in CV2 0.878	1.20	1.20	1.22
	Mesh level 1 to 2 2 to 3 3 to 4 4 to 5
	L 1 in CV1 0.61	0.60	0.63	0.62
	L 1 in CV2 0.67	0.65	0.67	0.63

Table 3 :

 3 Translation of a discontinuous function: order of the errors for various control cell and meshes.

Table 4 :

 4 Rotation of a regular function: order of the errors for various control cell and meshes, t = 1.50.

		1.36	1.62	1.52	1.78
	L 1 in CV2	1.47	1.71	1.67	1.82
	L ∞ in CV1 0.89	1.15	1.20	1.26
	L ∞ in CV2 1.00	1.27	1.19	1.26
	Mesh level 1 to 2 2 to 3 3 to 4 4 to 5
	L 1 in CV1 0.61	0.63	0.65	0.62
	L 1 in CV2 0.67	0.64	0.66	0.64

Table 5 :

 5 Rotation

Definition 2 Let us introduce k + ij ∈ V(i) and α i,k + ij ∈ [0, 1] such that :

Once again, this definition allows to define a density value at nodes N + ij by setting :

To define the interface reconstruction ρ + ij involved in (11), we first set

where Ψ is the flux limiter function, which is assumed to satisfy :

Then, ρ + ij is finally defined as :

2 Local Maximum Principle property

Preliminary results

Definition 3 (Convexity property) The second-order reconstruction ρ + ij is said to be convex if there exists θ

Definition 4 (τ -limiter) Let τ > 0; Ψ is said to be a τ -limiter if it satisfies the property :

Remark 1 This definition slightly differs from the "Q-limiter" property arising in the cell-centered context [START_REF] Clauzon | Analyse des schémas d'ordre élevé pour les écoulements compressibles. Application à la simulation numérique d'une torche à plasma[END_REF]. For example, instead of 12/7, a coefficient depending on the mesh regularity has to be introduced. Remark that the minmod limiter is a τ -limiter with τ = 1, and that several other less-diffusive usual one-dimensional limiters can be adapted to easily become τ -limiters (for example with τ = 2 for the superbee or the Van-Leer limiters [START_REF] Leveque | Numerical methods for conservation laws[END_REF], see [START_REF] Clauzon | Analyse des schémas d'ordre élevé pour les écoulements compressibles. Application à la simulation numérique d'une torche à plasma[END_REF] p. 81).

Lemma 1 If Ψ is a τ -limiter, then the reconstruction ( 21) is a convex reconstruction.

Proof. The proof is very similar to the one developped in [START_REF] Clauzon | Analyse des schémas d'ordre élevé pour les écoulements compressibles. Application à la simulation numérique d'une torche à plasma[END_REF] in the cell-centered context. Without loss of generality, let suppose that ρ i < ρ M + ij , so that p down ij > 0. If p up ij < 0, then the scheme degenerates to order 1 and ( 22) is obviously fulfilled with θ + ij = 0. If p up ij > 0, (21) ensures that ρ + ij ≥ ρ i . On the other hand, since from geometrical considerations

A i M + ij and using the fact that Ψ is a τ -limiter, we get :

Lemma 2 There exists some coefficients w + ijk ≥ 0, k ∈ V(i), such that

holds, and furthermore, they satisfy

Remark 2 This results is an adaptation of the so-called "Chainais-Hillairet condition" , see [START_REF] Chainais-Hillairet | Second order finite volume schemes for a nonlinear hyperbolic equation : error estimate[END_REF], developped in the cell-centered approach. But in the cell-centered framework, it relies on additional geometrical assumptions on the mesh. Adopting a vertex-based viewpoint relaxes the geometrical constraints needed to prove such a statement.

Proof. Let us define : ∆ρ

By a simple identification, the coefficients w + ijk can obviously be derived to fulfill (24): we set

in any other case. Now, in order to prove (25), we just have to remark that Ψ(r ij ) ≤ τ and use definition (9).

L ∞ stability

Theorem 1 Let Ψ be a τ -limiter. The numerical scheme (11) associated to the numerical flux (12) and the second-order reconstruction at the interfaces described in subsection 1.4.2 respects the local maximum principle under the CFL condition:

Proof. Making use of the consistency property (14), the numerical scheme (11) can be written as follows (remind that the indices n are dropped in the right-hand side)

where we defined : ∆ρ

, and by using the finite variations theorem, we get the existence of 0 < ζ ± ij < 1 such that :

where we introduced :