Caterina Calgaro

Emmanuel Creusé

Thierry Goudon

Stella Krell

Simulations of non homogeneous viscous flows with incompressibility constraints

Keywords: Non homogeneous viscous flows, Navier-Stokes equations, Mixtures, Multifluid flows, Finite volume methods

come   L'archive ouverte pluridisciplinaire

Introduction

The motivation of this work comes from the mathematical modelling of mixtures or multifluid flows. In fact, there exists a large variety of partial differential equations (PDEs) systems intended to describe such flows, and, of course, we will focus on some specific models. The situations we are interested in can be characterized by the following features:

• There are strong density variations in the flows, and the numerical challenge is to capture and to follow with accuracy the strong gradients and the fronts of density variations.

• The set of equations involves a constraint on the divergence of the velocity field, hereafter denoted by u. The simplest of these constraints is the solenoidal condition div x (u) = 0, but we shall see more intricate situations.

Our objective consists in designing dedicated numerical methods for the simulations of these models. The difficulty comes from the fact that the system couples equations of different type, and the constraint mentioned above. The constraint relies on modelling assumptions. In order to set up performing methods, it could be helpful to understand the origin of this relation. We will give some hints in this direction. The discussion should be taken with full awareness that many aspects in the derivation of the equations are not that neat, and can be considered as questionable: mathematical models in this field remain under debate. According to prescriptions in [START_REF] Prosperetti | Computational Methods for Multiphase Flows[END_REF], our viewpoint is therefore fully pragmatic: let us say that we are just picking a set of equations, and we try to discuss some mathematical properties, and to set up specific numerical schemes that allow us to investigate the sensitivity of the model with respect to variations of the (physical and numerical) parameters. However, we should bear in mind that the robustness of the conclusion should be considered with caution: changing a "detail" in the modelling assumptions might dramatically affect the mathematical structure of the model, which thus would require another approach. Let us end this introduction with a few words about potential applications of our approaches. Powder-snow avalanches is a relevant example of the kind of fluid mixtures we are considering and variations about the Navier-Stokes system have been used successfully to reproduce certain features of laboratory avalanches [START_REF] Dutykh | Mathematical modeling of powder-snow avalanche flows[END_REF][START_REF] Etienne | Simulation numérique directe de nuages aérosols denses sur des pentes ; application aux avalanches de neige poudreuse[END_REF][START_REF] Etienne | Numerical simulations of dense clouds on steep slopes: Application to powder-snow avalanches[END_REF]. These complex models also naturally arise in combustion theory; nuclear safety engineering provides further relevant examples of applications, see for instance the works [START_REF] Dellacherie | On a low Mach nuclear core model[END_REF][START_REF] Bernard | Study of a low Mach nuclear core model for singlephase flows[END_REF][START_REF] Penel | Theoretical study of an abstract bubble vibration model[END_REF]] motivated by security computations for PWR reactors. Finally, we mention that the reasonings of mixture theory have been developed to derive models describing biofilms formation [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF].

The paper is organized as follows. We start by discussing a few relevant mathematical models that couple transport, diffusion and constraint equa-tions. In Section 3 we present the hybrid finite volume (FV)-finite element (FE) scheme developed in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF][START_REF] Calgaro | L ∞ stability of vertex-based MUSCL finite volume schemes on unstructured grids; simulation of incompressible flows with high density ratios[END_REF][START_REF] Calgaro | Modeling and simulation of mixture flows: Application to powder-snow avalanches[END_REF]. The method is designed to reach the second order accuracy, at least for smooth solutions, owing to a suitable adaptation of Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL) techniques for the transport equation. Here, we pay attention to the treatment of the boundary conditions. Proceeding naively might lead to spurious instabilities when the geometry of the mesh and of the computational domain is non trivial. An alternative numerical method based on the Discrete Duality Finite Volume (DDFV) framework is introduced in Section 4. The DDFV framework is an attempt to cope with the difficulty in defining diffusion fluxes on interfaces without imposing geometric restrictions on the mesh construction, see [START_REF] Hermeline | A finite volume method for the approximation of diffusion operators on distorted meshes[END_REF][START_REF] Hermeline | Approximation of diffusion operators with discontinuous tensor coefficients on distorted meshes[END_REF][START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF]. The method has been extended to deal with Stokes and homogeneous Navier-Stokes equations [START_REF] Delcourte | A discrete duality finite volume approach to Hodge decomposition and div-curl problems on almost arbitrary two-dimensional meshes[END_REF][START_REF] Krell | Stabilized DDFV schemes for Stokes problem with variable viscosity on general 2d meshes[END_REF][START_REF] Krell | Finite volume method for general multifluid flows governed by the interface Stokes problem[END_REF][START_REF] Krell | Stabilized DDFV schemes for the incompressible Navier-Stokes equations[END_REF]. Dealing with non homogeneous flows requires a specific attention to design a relevant treatment of the convection terms of the system [START_REF] Goudon | DDFV scheme for incompressible Navier-Stokes equations with variable density[END_REF]. Section 5 offers a set of numerical experiments where we compare the performances of the two methods. In particular, we bring out difficulties related to the treatment of the boundary terms and the sensitivity to the mesh construction.

Examples of non-homogeneous flows involving constraints on

div x (u)

Incompressible flows

Let us start by recalling a few facts about the simplest situation of incompressible flow which means that the velocity is required to satisfy div x (u) = 0.

(

Neglecting any difficulty associated to a possible lack of regularity of u : R × R N → R N , let us consider the characteristic curves, defined by the ODE system d dt X(t, x) = u(t, X(t, x)), X(0, x) = x.

The quantity X(t, x) ∈ R N is nothing but the position at time t ≥ 0 of a particle, driven by the velocity field u, which starts at time t = 0 from the position x ∈ R N . Let us consider a fixed domain D 0 , and consider its image at a latter time t > 0 (See Figure 1.) This is a simple exercice (see [START_REF] Goudon | Intégration; Intégrale de Lebesgue et Introduction à l'Analyse Fonctionnelle[END_REF]Chap. 4]) to check that (1) implies d dt

D(t) = {X(t, x 0 ), x 0 ∈ D 0 }. D 0 • x D(t) • X(t, x) u(t, X(t, x))
D(t) dx = 0.
In other words, volumes are conserved by the flow. Next, let ρ : R × R N → [0, ∞) stand for the mass density of the fluid. It obeys the mass conservation equation

∂ t ρ + div x (ρu) = 0.
Expanding the space derivative, and taking into account (1), we find

∂ t ρ + u • ∇ x ρ + ρ div x (u) = ∂ t ρ + u • ∇ x ρ = 0.
Accordingly, the chain rule leads to

d dt ρ(t, X(t, x)) = 0.
The density is conserved along the characteristic curves, and we deduce that

ρ(t, X(t, x)) = ρ 0 (x).
As a particular case, an homogeneous fluid, satisfying ρ 0 (x) = ρ > 0 for any x, remains homogeneous for ever. Coupled to the momentum equation, say

1 ∂ t (ρu) + Div x (ρu ⊗ u) + ∇ x p = µ∆u + ρf ,
where ρf is the density of applied forces and µ > 0 is the dynamic viscosity of the fluid, this is the case which is usually referred to as the Incompressible Navier-Stokes system. The third unknown, the pressure p : R × R N → R, can be seen as a Lagrange multiplier associated to the constraint (1). However, there is no reason to restrict to homogeneous flows and we can deal with variable initial density as well. Compared to the huge literature on the homogeneous Navier-Stokes system, there are only a few results on the analysis of the Incompressible Navier-Stokes system with variable density, see e. g. [START_REF] Boyer | Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models[END_REF][START_REF] Danchin | Local and global well-posedness results for flows of inhomogeneous viscous fluids[END_REF][START_REF] Lions | Incompressible Models[END_REF] and a few attempts of numerical simulations of such flows [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF][START_REF] Calgaro | L ∞ stability of vertex-based MUSCL finite volume schemes on unstructured grids; simulation of incompressible flows with high density ratios[END_REF][START_REF] Guermond | A projection FEM for variable density incompressible flows[END_REF][START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF][START_REF] Johnston | Accurate, stable and efficient Navier-Stokes solvers based on explicit treatment of the pressure term[END_REF][START_REF] Pyo | Gauge-Uzawa methods for incompressible flows with variable density[END_REF]. Let us now discuss further models where more intricate constraints appear.

Zero-Mach flows

The starting point of the discussion is the compressible Navier-Stokes system

∂ t ρ + div x (ρu) = 0, ∂ t (ρu) + Div x (ρu ⊗ u) + ∇ x P = Div x (τ ) + ρf , ∂ t (ρE) + div x ((ρE + P )u) = div x (k∇ x T ) + div x (τ u) + ρf • u,
with τ = µ(∇u + ∇u - 2 3 div x (u)I). Additionally to the mass density ρ and the velocity field u, the unknowns involve the total energy E ≥ 0, which splits into kinetic energy and internal energy E = |u| 2 2 + e with a state law relating e, P and ρ: P = RρT = (γ -1)ρe, R being the constant of perfect gases and γ > 1 an exponent that characterizes the fluid under consideration. In order to write the equations in dimensionless form, we introduce density, time, length and velocity units ρ , T, L, U = L/T respectively. Let E be the energy unit and we define the temperature unit as γ-1 R E and the pressure unit is P = ρ E . We can now introduce the Mach number

Ma = U √ E = ρ U 2 P .
1 Here we adopt different notations for the divergence of a vector field div x (u) = N j=1 ∂ xj u j and the divergence of a matrix-valued function [Div

x (A)] i = N j=1 ∂ xj A ij .
We assume that all the other dimensionless parameters -namely the Reynolds, Froude and Prandtl numbers 2 -that characterize the flow have the same order of magnitude and we consider the regime where the Mach number is small: Ma Re Fr Pr.

Keeping only Ma in the dimensionless form of the equation, we are led to

∂ t ρ Ma + div x (ρ Ma u Ma ) = 0, ∂ t (ρ Ma u Ma ) + Div x (ρ Ma u Ma ⊗ u Ma ) + 1 Ma 2 ∇ x P Ma = Div x (τ Ma ) + ρ Ma f , ∂ t (ρ Ma E Ma ) + div x ((ρ Ma E Ma + P Ma )u Ma ) = div x (k∇ x T Ma ) + Ma 2 (div x (τ Ma u Ma ) + ρ Ma f • u Ma ).
Furthermore, the Mach number also enters into the definition of the total energy

E Ma = Ma 2 γ -1 |u Ma | 2 2 + T Ma γ -1 , P Ma = ρ Ma T Ma . (2)
Let us assume that all unknowns admit limit ρ, u, T , etc... in a strong enough sense so that we can deal with the non linear terms. As Ma → 0 the momentum equation degenerates to ∇ x P = 0. Up to a convenient assumption on the boundary condition, we infer that the pressure does not depend on time and space, namely we have ρT (t, x) = P (constant). The mass conservation equation becomes ∂ t ρ + div x (ρu) = 0. Next, we multiply the momentum equation by a divergence free trial function -namely a vector valued function ϕ ∈ C ∞ c , with div x (ϕ) = 0 -in order to get rid of the stiff term (since

∇ x P Ma • ϕ dx = -P Ma div x (ϕ) dx). Letting Ma → 0 yields ∂ t (ρu) + Div x (ρu ⊗ u)) -ρf -Div x (τ ) • ϕ dx = 0, (3)
for any such trial function. Accordingly, we obtain

∂ t (ρu) + Div x (ρu ⊗ u) + ∇ x π = ρf + Div x (τ ), 2
We remind the reader that these numbers are defined by Re = ρ UL μ , Fr = U √ f L , Pr = μ ρ k respectively, with μ, k and f typical values for the dynamic viscosity, the thermal diffusivity and the force density (in most cases the latter is the gravity acceleration).

where the multiplier π appears as a consequence of the fact that (3) only holds for a restricted set of trial functions. We turn to the energy equation, bearing in mind [START_REF] Alazard | A minicourse on the low Mach number limit, Discrete Contin[END_REF]. In particular ρ Ma E Ma tends to 1 γ-1 ρT = P γ-1 . Therefore the limiting energy equation reads

γP γ -1 div x (u) = div x (k∇ x T ) = P div x (k∇ x (1/ρ)).
Namely, we obtain a constraint where the divergence of the velocity is related to derivatives of a function of the density. The analysis of such an asymptotic regime is quite challenging; we refer the reader to [START_REF] Alazard | A minicourse on the low Mach number limit, Discrete Contin[END_REF][START_REF] Alazard | Low Mach number flows and combustion[END_REF]. Note that it is a different question to search for a numerical scheme for the limiting system only or to construct a method for the compressible model which is able to handle correctly the regime Ma 1 and which avoids prohibitive stability conditions. On this aspect we refer the reader to [START_REF] Schneider | Extension of finite volume compressible flow solvers to multidimensional, variable density zero Mach number flows[END_REF][START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF][START_REF] Guillard | On the behavior of upwind schemes in the low Mach number limit[END_REF][START_REF] Haack | An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations[END_REF][START_REF] Herbin | Staggered schemes for all speed flows[END_REF][START_REF] Zaza | Contribution à la résolution numérique d'écoulements à tout nombre de Mach et au couplage fluide-poreux en vue de la simulation d'écoulements diphasiques homogénéisés dans les composants nucléaires[END_REF]. We also mention the derivation of further models corresponding to the Zero-Mach regime in [START_REF] Dellacherie | On a low Mach nuclear core model[END_REF][START_REF] Bernard | Study of a low Mach nuclear core model for singlephase flows[END_REF][START_REF] Penel | Theoretical study of an abstract bubble vibration model[END_REF]].

Modelling of mixtures: the Kazhikhov-Smagulov system

We are interested in flows that can be seen as a mixture of a disperse phase and a carrying phase, with respective mass densities ρ d , ρ f . Let u d , u f stand for their respective velocities. We start by writing the mass conservation equation for both phases:

∂ t ρ d,f + div x (ρ d,f u d,f ) = 0. ( 4
)
Adding these equations, we obtain

∂ t ρ + div x (ρu) = 0, (5)
where

ρ = ρ d + ρ f , ρu(t, x) = (ρ f u f + ρ d u d )(t, x),
define the mean mass density ρ and the mean mass velocity u respectively. Next, we write the following momentum balance equation for ρu

∂ t (ρu) + Div x (ρ u ⊗ u) + ∇ x p = ρg + Div x (µD(u)),
with D(u) = ∇ x u + ∇ x u . In the right side the force term ρg is proportional to the gravitational field g. In order to close the equations, we need to explain how the pressure field p is defined.

We bear in mind that we adopt an averaged description of the flow, where we consider that at each time t and each position x, the two phases can be present. The local amount of the disperse phase is evaluated through the volume fraction 0 ≤ φ(t, x) ≤ 1, which can be intuitively thought of through the formula φ(t, x) = lim r→0 Volume occupied by the disperse phase in B(x, r)

|B(x, r)| .

The system is closed by specifying the modelling assumption. Here and below we assume the incompressibility of both phase, in the sense that, if a single phase were occupying a given location, then its mass per unit volume is constant (given by ρd , ρf ). Accordingly, we simply have

ρ d (t, x) = ρd φ(t, x), ρ f (t, x) = ρf (1 -φ(t, x)).
Let us go back to (4) with this assumption. It yields

∂ t ρ f ρf + ρ d ρd = ∂ t (1 -φ + φ) = 0 = -div x ρ f u f ρf + ρ d u d ρd :=v .
It means that the mean volume velocity

v(t, x) = 1 -φ(t, x) u f (t, x) + φ(t, x)u d (t, x),
is divergence-free.

Remark 2.1. The distinction between the mean volume velocity and the mean mass velocity can be crucial in fluid mechanics: in a series of papers Howard Brenner [START_REF] Brenner | Unsolved problems in fluid mechanics: On the historical misconception of fluid velocity as mass motion, rather than volume motion, communication for the 100th anniversary of the Ohio State Chem[END_REF][START_REF] Brenner | Navier-Stokes revisited[END_REF] proposed to revisit the derivation of Navier-Stokes equations based on this distinction. It turns out that this approach is highly controversial, see e. g. [54]. However, we point out that Brenner's models have been showed to exhibit nice mathematical structure and stability properties [START_REF] Feireisl | New perspectives in fluid dynamics: Mathematical analysis of a model proposed by Howard Brenner[END_REF]. Moreover, similar corrections to the standard Navier-Stokes systems have been proposed recently motivated by numerical purposes [START_REF] Guermond | Viscous regularization of the Euler equations and entropy principles[END_REF].

Having observed this fact, the system of PDEs is closed by assuming a Fick's law that relates the mean mass velocity u, the mean volume velocity v and gradient of the density ρ:

u = v -κ∇ x ln(ρ), (6)
with a certain coefficient κ > 0. This assumption together with the fact that v is solenoidal yields div x (u) = -div x (κ∇ x ln(ρ)).

It is worth remarking that the closure relation can be expressed equivalently on the carrier fluid velocity u f = u -κ 1 ∇ x ln(ρ f /ρ) or as an evolution equation, of convection-diffusion type, for the disperse volume fraction

∂ t φ + div x (φu) = div x (κ 2 ∇ x ln(ρ f + (ρ d -ρf )φ),
with positive coefficients κ 1 or κ 2 . These observations lead to different options to choose the numerical unknowns. However, the equivalence between different formulations, which is clear for the continuous equations, might fail for the corresponding discrete systems, see [START_REF] Berthelin | Multifluid flows: a kinetic approach[END_REF].

For the introduction of such models for mixtures, and application to pollution spreading, we refer the reader to Kazhikhov-Smagulov [START_REF] Kazhikhov | The correctness of boundary value problems in a diffusion model in an inhomogeneous fluid[END_REF], see also [START_REF] Franchi | A comparison of Graffi and Kazhikov-Smagulov models for top heavy pollution instability[END_REF][START_REF] Franchi | Convection, diffusion and pollution: The model of Dario Graffi[END_REF][START_REF] Graffi | Il teorema di unicitá per i fluidi incompressibili, perfetti, eterogenei[END_REF][START_REF] Joseph | Fundamentals of Two-Fluid Dynamics. Part II: Lubricated Transport, Drops and Miscible Liquids[END_REF]. The analysis of this system dates back to [START_REF] Antontsev | Boundary Value Problems in Mechanics of Nonhomogeneous Fluids[END_REF][START_REF] Beirao | Diffusion on viscous fluids. Existence and asymptotic properties of solutions[END_REF][START_REF] Secchi | On the initial value problem for the equations of motion of viscous incompressible fluids in the presence of diffusion[END_REF]. Imposing a specific relation between the viscosity µ and the diffusion coefficient κ further estimates (roughly speaking that control ∇ x ρ) can be derived which are useful to justify the existence of solutions [START_REF] Bresch | Effect of density dependent viscosities on multiphasic incompressible fluid models[END_REF][START_REF] Liao | Quelques résultats mathématiques sur les gaz à faible nombre de Mach[END_REF][START_REF] Liao | A global existence result for a zero Mach number system[END_REF]. In [START_REF] Calgaro | Modeling and simulation of mixture flows: Application to powder-snow avalanches[END_REF][START_REF] Goudon | On a model for mixture flows: Derivation, dissipation and stability properties[END_REF], we reinterpret the derivation of the KS equations by means of hydrodynamic regimes, starting from a coupled fluid-kinetic system. This viewpoint allows us to derive a new model which can be interpreted as a generalized KS system and which has remarkable dissipation properties.

A hierarchy of models

It looks tempting to get rid of u and to make use of the relation div x (v) = 0, which is a more usual constraint. We plug (6) into [START_REF] Bell | A second-order projection method for variabledensity flows[END_REF]; it makes a diffusion term appear since

div x ρ κ∇ x ln(ρ) = κdiv x ρ ∇ x ρ ρ = κ∆ x ρ.
Similarly, we rewrite the momentum equation with the mean volume velocity v instead of u, see [START_REF] Calgaro | Modeling and simulation of mixture flows: Application to powder-snow avalanches[END_REF]. We obtain the following system for ρ, v and p

∂ t ρ + div x (ρv) = κ∆ x ρ, ρ ∂ t v + (v • ∇ x )v + ∇ x p = ρg + Div x (µD(v)) +κ(∇ x v -∇ x v )∇ x ρ +κ 2 ∇ x ∆ x ρ -Div x ∇ x ρ ⊗ ∇ x ρ ρ -κDiv x µD 2
x ln(ρ) .

Choosing the mean volume velocity v instead of the mean mass velocity u has, at least, two consequences. Firstly, it changes the mass conservation equation, which becomes a convection-diffusion equation instead of a mere transport equation. Secondly it makes several ugly terms appear in the momentum equation, that involve high-order derivatives of the density. Furthermore, in this form, we also realize that the "additional terms" induce intricate boundary conditions: in fact the Fick law hides further modelling issues related to the definition of the boundary conditions.

Let us point out that various simplifications are relevant:

• Assume a simple expression for µ, say constant or an affine function of ρ,

• Get rid of higher order terms by an asymptotic argument; as it appears in [START_REF] Kazhikhov | The correctness of boundary value problems in a diffusion model in an inhomogeneous fluid[END_REF],

• Assume an ad hoc relation between µ and κ: with µ = µ 0 + μρ and κ = μ, the higher order terms cancel out and, furthermore a remarkable energy identity can be established, as observed in [START_REF] Dutykh | Mathematical modeling of powder-snow avalanche flows[END_REF][START_REF] Bresch | Effect of density dependent viscosities on multiphasic incompressible fluid models[END_REF][START_REF] Danchin | On the well-posedness of the full low-Mach number limit system in general critical Besov spaces[END_REF][START_REF] Liao | A global existence result for a zero Mach number system[END_REF].

Example: simulations of powder-snow avalanches

Before entering into details of the numerical methods, let us discuss applications of Kazhikhov-Smagulov (KS) equations to the description of powdersnow avalanches. To this end, we rewrite the system in the following dimen-sionless form:

∂ t ρ + div x (ρv) = 1 Re Sc ∆ x ρ, ρ ∂ t v + (v • ∇ x )v + ∇ x p = 1 Fr 2 ρg + 1 Re Div x (µ(ρ)D(v)) + 1 Re Sc (∇ x v -∇ x v )∇ x ρ, div x (v) = 0,
where we make the following dimensionless parameters appear:

Reynolds number: Re = U L ν , Froude number: Fr = U(gL) -1/2 , Schmidt number: Sc = ν κ .
Here, L and U stand for length and velocity units, respectively, g gives the value of the gravitational acceleration, while ν is a typical value of the kinematic viscosity of the flow. Engineers also introduce the densimetric Froude number

Fr d = U gL ∆ρ ρd -1/2
, which also depends on the extreme densities within the flow ρd , with ∆ρ = ρf -ρd . Note that we have adopted the simplifying assumptions that allow us to get rid of terms with higher-order derivatives of the density. As a matter of fact, we remark that Sc = +∞ gives the usual Incompressible Navier-Stokes system. As said above, the system requires boundary conditions, on the density and the velocity; in what follows, simulations are presented with

v ∂Ω = 0, ∇ x ρ • ν ∂Ω = 0,
where ν is the unit outward normal of the domain.

The mathematical modelling of powder-snow avalanches, like many other mixture flows, is still widely open. Nevertheless, we have at hand experimental data and numerical simulations (based on different systems of PDEs though), which can be used to discuss the relevance of the results and the role of the different parameters. We use experimental data obtained at Cemagref-IRSTEA [START_REF] Naaim-Bouvet | Physical modelling of the interaction between powder avalanches and defence structures[END_REF][START_REF] Naaim-Bouvet | Numerical and physical modelling of the effect of a dam on powder avalanche motion: Comparison with previous approaches[END_REF][START_REF] Naaim-Bouvet | Approche macro-structurelles des écoulements biphasiques turbulents de neige et de leur interaction avec des obstacles[END_REF][START_REF] Clément-Rastello | Etude de la dynamique des avalanches de neige en aérosol[END_REF]. For real avalanches, we have Re 10 8 , and Fr d 1 which make this situation not affordable for simulation. But laboratory avalanches have been produced with less extreme parameters: Re 10 3 -10 4 , Fr d

1. Data are available but the comparison should still be considered with caution since many crucial parameters are not really well-known like the Schmidt number Sc, the mass densities of the constituents, the viscosity coefficient of the mixture, etc. Further discussions would be necessary on the boundary conditions to decide whether or not it is necessary to consider friction effects. Finally, the method and the quality of measurements could be also questionable, with several sources of uncertainties. The experimental results have been compared to simulations, performed on a k/ model, with tuned parameters, by using the Fluent software packages, on fixed Cartesian grids. Another set of simulations appeared in [START_REF] Etienne | Simulation numérique directe de nuages aérosols denses sur des pentes ; application aux avalanches de neige poudreuse[END_REF][START_REF] Etienne | Numerical simulations of dense clouds on steep slopes: Application to powder-snow avalanches[END_REF], based on a version of the KS system, with a Finite Element scheme coupled to characteristic methods for the convection terms; the simulations also use mesh refinements strategies in order to follow the advance of the front. Finally [START_REF] Dutykh | Mathematical modeling of powder-snow avalanche flows[END_REF] uses OpenFoam, on a fixed Cartesian mesh to produce numerical avalanches by using the KS system.

In [START_REF] Calgaro | Modeling and simulation of mixture flows: Application to powder-snow avalanches[END_REF], we simulate the KS system on the initial configuration represented in Figure 2. We use the hybrid FV-FE method described in Section 3, coupled to a mesh refinement method: see Figure 3 for a typical example of mesh adapted to the displacement of the avalanche. We deal with different sets of parameters:

• From numerical experiments in [START_REF] Dutykh | Mathematical modeling of powder-snow avalanche flows[END_REF] -Fr d = 1, Fr = 4.3589, Re = 2242, Sc = 1.0.

-Mesh refinements go up to 15.000 triangles, with smallest convex radius h min ≈ 7 × 10 -4 m (the domain size is 2.7 × 0.8 m and the initial avalanche size 0.3 × 0.3 m).

-We obtain very similar results (up to a slight change of units which need to be corrected in [START_REF] Dutykh | Mathematical modeling of powder-snow avalanche flows[END_REF]).

• From laboratory experiments in [START_REF] Naaim-Bouvet | Physical modelling of the interaction between powder avalanches and defence structures[END_REF][START_REF] Naaim-Bouvet | Numerical and physical modelling of the effect of a dam on powder avalanche motion: Comparison with previous approaches[END_REF][START_REF] Naaim-Bouvet | Approche macro-structurelles des écoulements biphasiques turbulents de neige et de leur interaction avec des obstacles[END_REF][START_REF] Clément-Rastello | Etude de la dynamique des avalanches de neige en aérosol[END_REF] -Fr d = 0.83, Fr = 0.3725, Re = 1000, Sc = 1.

-Mesh refinements go up to 50.000 triangles, with smallest convex radius h min ≈ 1.5 × 10 -4 m and maximum length of the edges h max ≈ 10 -2 m (domain size 2 × 0.5 m and initial height of the avalanche L= 0.09 m).

-The order of magnitude of the horizontal velocity agrees with the value recorded form the experiments (0.35 m s -1 ).

The simulations permit us to bring out the following qualitative results:

• Variations of Sc do not change significantly the velocity of the front. But the density profiles are modified: the KS model produces results noticeably different from the mere incompressible system.

• Numerical difficulty increases as Re increases or Fr decreases.

• Numerical results are qualitatively consistent with the experiments.

• The Reynolds number affects the details of the flow, but does not significantly influence the speed of the front.

• The Froude number greatly influences maximal speeds.

• The maximal velocity exceeds the front speed by 30% to 40%. In the following sections, we are going to describe numerical strategies to deal with such equations, restricting to the simpler case of incompressible, but non-homogeneous, flows.

h l 0 d s l l s h s h 0 θ g ρ + , ν + ρ -, ν - ρ, ν

A hybrid finite volume-finite element scheme

Principles of the method

Let us start by discussing the numerical simulation of the incompressible Navier-Stokes system in a two-dimensional domain Ω. We work with a (conformal) tessellation τ h made of triangles T , and we define Ω h = ∪ T ∈τ h T . Let us note that in the case where Ω is a polygonal domain, then τ h is defined such that Ω h = Ω whatever the tessellation τ h chosen. We associate to this primal mesh a dual mesh: the vertices of the primal mesh are considered as the centers of the cells of the dual mesh. As we shall see below, there are several options to construct the dual mesh: in Figure 4 it is obtained by joining the barycenters of the triangles to the middle of the edges.

Our scheme works by time-splitting the mass conservation and the momentum balance equation and by dealing with the equations on staggered grids defined on Ω h :

• The mass conservation

∂ t ρ + div x (ρu) = 0, (7)
is treated by a finite volume method, with cells defined by the dual mesh; • The momentum equation, together with the divergence free constraint

A 3 A 4 A 5 A 12 A 23 A 34 A 45 A 51 = A 15 A C A A 1 A 12 A 2 A 23 A 3 A 34 A 4 A 45 A 5 A 51 = A 15 Γ - 1 Γ + 1 Γ - 2 Γ + 2 Γ - 3 Γ + 3 Γ - 4 Γ + 4 Γ - 5 Γ + 5 n
∂ t (ρu) + Div x (ρu ⊗ u) -µ∆ x u + ∇ x p = ρf , div x (u) = 0, ( 8
)
is treated by a Finite Element method defined on the primal mesh (following the quite usual framework used for homogeneous flows). Hence (u, p) are approached by continuous and piecewise polynomials functions (say u is P 2 , and p is P 1 over the triangles).

Both parts use well-established techniques, but the difficulty consists precisely in the transmission between the two steps and the management of the unknowns defined on different grids and updated with different methods. This viewpoint allows us to make use of performing methods, that can be developed independently, for the two equations which have a different type. Therefore, the velocity and the pressure are updated by solving the (linear) system

ρ δt (u k+1 -u k ) + ρ u k • ∇ x u k+1 -µ∆ x u k+1 + ∇ x p k+1 = ρ f k+1 , div x (u k+1 ) = 0, (9)
where we are faced the question of how to define the density ρ to be used in the FE scheme to compute u. To this end, we have at hand a density field which is constant over the control volumes. Thus, given a triangle, we know 3 values of densities (ρ k+1 ); hence we can naturally define a P 1 reconstruction: this is the ρ that arises in [START_REF] Brenner | Unsolved problems in fluid mechanics: On the historical misconception of fluid velocity as mass motion, rather than volume motion, communication for the 100th anniversary of the Ohio State Chem[END_REF]. The update of (u, p) can be seen as a saddle point problem which, once discretized with respect to space, is expressed in the following matrix form

A B B T 0 u k+1 p k+1 = ρ f k+1 0 .
We point out that non homogeneities of the flows make the structure of the matrix A more complicated since ρ varies from place to place. If we think in terms of the underlying elliptic equation for the pressure (as usual when using projection methods), this equation is now with variable coefficients, which might lead to preconditioning issues, see [START_REF] Calgaro | Incremental incomplete LU factorizations with applications to time-dependent PDEs[END_REF].

The discussion of the mass conservation equation is more subtle. We consider a control volume C A with center A (see Figure 4). The numerical unknown

ρ k A is intended to approximate 1 |C A | C A ρ(t k , x) dx.
Integrating the the mass conservation equation over C A yields

ρ k+1 A -ρ k A δt + 1 |C A | i Γ i ρ k u • n dσ(x) = 0,
where

∂C A = ∪ nb i=1 Γ i , with Γ i = Γ + i ∪ Γ - i
, nb the number of neighbour nodes of A (respectively denoted A i , 1 ≤ i ≤ nb), and n the unit outward normal to ∂C A . Here, we address the question, given u piecewise polynomial on the triangles, to define a suitable u , piecewise constant on ∂C A , to obtain a relevant approximation of u • n on ∂C A .

To construct the numerical flux, the basic requirement is to preserve homogeneous states: if ρ is constant over the whole computational domain, it should remain constant forever. The naive construction which simply uses the evaluation of the FE function u, P 2 per triangle, on the interface ∂C A , fails to meet this objective, and we propose another construction. We bear in mind that the divergence free condition for the FE method actually reads

ω A div x (u)ψ A dx = 0 = i T i div x (u)ψ A dx,
for all P 1 basis function ψ A associated to the node A with supp(ψ A ) = ω A = ∪ nb i=1 T i . It turns out that this formula does not involve the value of u at the nodes A i , but only at the nodes A i and A ii+1 (1 ≤ i ≤ nb), corresponding to the degrees of freedom associated to the middle of the edges. Thus, we seek u + i =u - i+1 the velocity on Γ + i ∪ Γ - i+1 , as a convex combination of u i,i+1 , u i and u i+1 such that :

i |Γ - i |u - i • n |Γ - i + |Γ + i |u + i • n |Γ + i = 0,
this expression corresponding to the finite volume interpretation of the divergence free condition. It yields some algebraic relations and we end up with

u + i = u - i+1 = 1 3 u i + u i,i+1 + u i+1 .
As a matter of fact it is worth observing that

• u + i = u(A i i+1 ),
• The formula for u + i coincides with a linear interpolation at A i i+1 from u at the nodes A i , A i+1 and A i i+1 .

With the value of u + i at hand, a simple upwind finite volume scheme would naturally lead to

Γ + i ρ k u • n dσ(x) =    ρ k A |Γ + i | u + i • n if u + i • n > 0, ρ k A i |Γ + i | u + i • n if u + i • n < 0, (10)
and

Γ - i+1 ρ k u • n dσ(x) =    ρ k A |Γ - i+1 | u - i+1 • n if u - i+1 • n > 0, ρ k A i+1 |Γ - i+1 | u - i+1 • n if u - i+1 • n < 0. ( 11
)
3.2. Extension: MUSCL scheme on unstructured meshes, maximum principle and multislope methods Using upwind fluxes such as [START_REF] Brenner | Navier-Stokes revisited[END_REF] and [START_REF] Bresch | Effect of density dependent viscosities on multiphasic incompressible fluid models[END_REF] restricts the method to first order accuracy. To improve the accuracy, fluxes can be designed based on MUSCL strategies. While the method is completely clear in 1D or on Cartesian grids, difficulties arise when dealing with general unstructured meshes.

Following ideas introduced 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][START_REF] Clain | L ∞ stability of the MUSCL methods[END_REF] for Cell-Center methods, we design a multislope method for Vertex-Based methods, defining directional derivatives and limiters on each interfaces. We refer the reader to Figure 5, and we focus on the treatment to be done on Γ + 1 . A similar reasoning applies to all the other components of ∂C A . Let us denote n + 1 the unit outward normal of C A on Γ + 1 and let us suppose that u

+ 1 • n + 1 > 0. Let A + 1 be the mid-point of [A 1 A 12 ]. Our goal is to derive the value of ρ k A + 1
, to be used instead of ρ k A in the definition of the flux [START_REF] Brenner | Navier-Stokes revisited[END_REF], in order to increase the accuracy of the Finite Volume scheme. Let M + 1 and N + 1 be the two intersection points between (AA + 1 ) and ∂ω A (see Figure 5). We define

p up,+ 1 = ρ A -ρ N + 1 AN + 1 , and
p down,+ 1 = ρ M + 1 -ρ A AM + 1 .
Then the density is evaluated at node A + 1 by

ρ A + 1 = ρ A + p + 1 AA + 1 , with p + 1 = p up,+ 1 Lim p down,+ 1 p up,+ 1 ,
where Lim is a so-called "τ -limiter". For the convection equation, the L ∞ stability of the scheme can be established and improved accuracy is observed on numerical validations, see [START_REF] Calgaro | L ∞ stability of vertex-based MUSCL finite volume schemes on unstructured grids; simulation of incompressible flows with high density ratios[END_REF]. Eventually, the method can be coupled to mesh-refinement algorithms in order to follow strong density gradients. This scheme has been adapted to handle the additional terms of the KS system and to perform the computations of the numerical avalanches, see [START_REF] Calgaro | Modeling and simulation of mixture flows: Application to powder-snow avalanches[END_REF].

Boundary conditions

The treatment of boundary conditions is seldom addressed in details. However, it gives rise to important practical issues, in particular when we wish to maintain a high-order accuracy. We distinguish the cases of "Dirichlet" or "Neumann" boundary conditions to be imposed on the density. The former refers to the usual situation where the inflow density is imposed:

A A 1 A 2 A 5 M + 1 N + 1 A 3 A 4 T 1 T 5 n + 1 A 1 A 12 A + 1 Γ + 1
ρ(t, x) = ρ b (t, x) for (t, x) ∈ (0, T ) × ∂Ω such that u • ν < 0.
In this case, we just specify the prescribed value at the boundary ∂Ω h . As classically known, the approximation of the exact value of ρ b on ∂Ω given by an interpolated value ρ h on ∂Ω h introduces an additional error of order O(h 2 ). The latter refers to another common situation where either the boundary condition is intended to describe the reflection of the fluid by the boundary ∂Ω (wall boundary conditions) or it accounts for symmetry conditions, which are used to reduce the computational domain. This situation occurs for instance with the simulation of Rayleigh-Taylor instabilities, that will be presented in Section 5. We are going to discuss how the numerical fluxes should be adapted to this situation, in the MUSCL framework presented above.

We refer the reader to Figure 6: we focus on the cell with center A ∈ ∂Ω h , where the segments [A 2 A] and [AA 5 ] of the primal mesh belong to the boundary ∂Ω h . The flux through ∂Ω h ∩ ∂C A is evaluated as mentionned in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF]. In particular, the normal velocity at the boundary has to be carefully evaluated, see [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF]Section 2.4.3]. According to its sign, the normal flux of the density is evaluated from the numerical solution (outgoing flux) or by using the boundary condition (incoming flux). Obviously, in the degenerate case u • n A = 0, no flux is added. However, it remains to evaluate the flux density on the interior interfaces of ∂C A : for instance on Γ + 1 we need to define the value of ρ A + 1 . Going back to the construction in Section 3.2, the difficulty comes from the definition of the upstream gradient p up,+ 1 since the node N + 1 does not exist. A first solution is to let the scheme degenerate at order one, by simply setting ρ A + 1 = ρ A . We will see in Section 5 that this loss of accuracy can generate instabilities in the close vicinity of the boundary, leading for some configurations to unphysical phenomena. This effect is sensitive to the geometry of the domain and the shape of the mesh. As an alternative, we propose a reflection technique, which uses a suitable definition of "ghost cells", see Figure 6 which can be easily generalized for other mesh configurations. This method preserves the second order accuracy, as it will be shown by the numerical experiments. Let n A be the unit outward vector to ∂Ω at A, and let (N A ) be its orthogonal line at A. We define Ã2 , Ã1 and Ã5 as the points orthogonally symmetric with respect to (N A ) of A 2 , A 1 and A 5 respectively. Consequently this construction defines the ghost triangles T1 and T5 . In order to enforce homogeneous Neumann boundary conditions for ρ on ∂Ω, we set ρ Ã2 = ρ A 2 , ρ Ã1 = ρ A 1 and ρ Ã5 = ρ A 5 . Going back to the definition in Section 3.2, it allows us to define an upstream gradient, by using now the value of the density defined by linear interpolation at the fictitious node Ñ + 1 .

Remark 3.1. The geometry of the domain and the shape of the control volume are important in this discussion. In [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF] some simulations are performed by considering a mesh obtained by cutting through diagonals a preliminary tessellation made of squares, see Figure 12 below. On this structured mesh, control volumes of the dual mesh can be obtained by joining directly the mid- points of the vertices of the triangles, see [14, Fig. 2-(b)]. For this very specific mesh, the two options actually coincide. Indeed, using the ghost cells, we find that the upstream and downstream gradients, in the normal direction, are exactly the opposite of each other, and the limiter makes the scheme degenerate to first order.

A A 1 A 2 A 5 Ã2 Ã1 Ã5 ∂Ω ∂Ω M + 1 Ñ + 1 T 1 T 5 T1 T5 n + 1 A + 1 Γ + 1 n A (N A )

A DDFV method for non homogeneous viscous flows

Still having in mind the possible adaptations to handle complex models with non homogeneous constraints, let us describe another approach, where the mass conservation and the momentum balance equations are both treated by a finite volume method. It offers the advantage of a unified viewpoint in the discretization method. To this end we adopt the Discrete Duality Finite Volume (DDFV) framework, which, again, leads to consider unknowns on staggered grids. When approximating diffusion equations by a finite volume approach, we face the difficulty of defining the normal derivative at the interfaces of the control volumes. Except in very specific geometric situations, the information stored at the center of the volumes is not enough for that purpose. After a pioneering attempt of Y. Coudière, J.-P. Vila, P. Villedieu, the basis of the DDFV method has been set up by F. Hermeline [START_REF] Hermeline | A finite volume method for the approximation of diffusion operators on distorted meshes[END_REF][START_REF] Hermeline | Approximation of diffusion operators with discontinuous tensor coefficients on distorted meshes[END_REF] and K. Domelevo, P. Omnes [START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF]. The idea is two-fold: we increase the number of numerical unknowns so that a full discrete gradient can be defined, and we define discrete operators so that the usual duality formula that are obtained by integration by parts at the continuous level, are preserved. Since then, it has been the object of many extensions, including for the Stokes system [START_REF] Delcourte | A discrete duality finite volume approach to Hodge decomposition and div-curl problems on almost arbitrary two-dimensional meshes[END_REF][START_REF] Krell | Stabilized DDFV schemes for Stokes problem with variable viscosity on general 2d meshes[END_REF][START_REF] Krell | Finite volume method for general multifluid flows governed by the interface Stokes problem[END_REF][START_REF] Krell | Stabilized DDFV schemes for the incompressible Navier-Stokes equations[END_REF]. There exists also 3D versions of the method [START_REF] Krell | The discrete duality finite volume method for the Stokes equation on 3D polyhedral meshes[END_REF], which opens perspectives to deal with our scheme for non homogeneous flows in higher dimension as well.

Definition of the meshes and discrete operators

We refer the reader to [START_REF] Krell | Stabilized DDFV schemes for Stokes problem with variable viscosity on general 2d meshes[END_REF][START_REF] Krell | Stabilized DDFV schemes for the incompressible Navier-Stokes equations[END_REF] for an exhaustive description of the notations of the DDFV method. Let us explain the construction of the meshes. We start with a tessellation, the so-called primal mesh denoted by M. Its construction can be quite general: we can mix triangles and quadrangles, possibly with non-conformal elements, etc: see Figure 7-(a). With unknowns stored at x K and x L (the center of the control volume K and L) we can construct a discrete gradient, in the direction x K , x L ; but since x K , x L is not assumed to be orthogonal to the interface of the control volume, this is not enough to define the normal derivative at the interface. Therefore, we additionally consider numerical unknowns stored at the vertices x K * , see Figure 7-(b). Now, the subdomain defined by the four vertices x K , x K * , x L , x L * defines another control volume, the so-called diamond cell D. The tessellation obtained this way is referred to as the diamond mesh D; it is represented on Figure 7-(c). Finally, since we are dealing with unknowns stored at the vertices x K * of the primal mesh, we need to construct a control volume K * associated to these points. This is the dual mesh M * that can be obtained either by joining the centers x K of the primal mesh, see Figure 7-(d), or by joining the centers x K to the mid-point of the edges.

In what follows, we denote by T = (M, M * ) the sets of primal and dual cells. Given a cell C (resp. an edge), m C stands for the volume of this cell (resp. the length of the edge). The velocity unknowns will be stored on the primal u M = (u K ) K∈M and dual u M * = (u K * ) K * ∈M * meshes: we shall denote u T = (u M , u M * ). The pressure and density unknowns will be stored on the diamond mesh: ρ D = (ρ D ) D∈D and p D = (p D ) D∈D respectively. Next, we introduce discrete differential operators. Given a diamond cell D whose vertices are x K , x K * , x L , x L * , we can define discrete derivatives in the directions x K , x L and x K * , x L * :

⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ x K x L (a) Primal mesh M ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ x K x L x K * x L * (b) Diamond cell D ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ x K x L x K * x L * (c) Diamond mesh D ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ x K x L x K * x L * (d) Dual mesh M *
∇ D u T (x L -x K ) = u L -u K , ∇ D u T (x L * -x K * ) = u L * -u K * .
Since the vectors x K , x L and x K * , x L * are not colinear, we can reconstruct, from u T defined on the dual and primal meshes, a full discrete Jacobian matrix, defined on the diamond mesh:

∇ D u T = 1 2m D m σ (u L -u K ) ⊗ n σK + m σ * (u L * -u K * ) ⊗ n σ * K * ,
with m σ (resp. m σ * ) the length of (x K * , x L * ) (resp. (x K , x L )), n σK the unit vector orthogonal to the oriented edge σ = [x K * , x L * ], etc. The viscous term in the Navier-Stokes equation involves the symmetric part of this matrix, and the divergence can be defined as the trace:

D D u T = 1 2 (∇ D u T + (∇ D u T ) ) , div D u T = Tr ∇ D u T .
Next, in order to obtain a discrete analog of the Stokes formula we introduce the discrete dual operator

div K ξ D := 1 m K D∈D K m σ ξ D n σK ,
which picks a matrix-valued quantity defined on the diamond mesh, and returns a vector defined on the primal and dual cells.

Treatment of the convection equation

Let us start by explaining how we handle the mass conservation equation. We remind the reader that the discrete mass density is intended to be an approximation of the average over the diamond cells D ρ(t, x) dx. Therefore, we wish to mimic the flux

D div(ρu) = s∈∂D s ρu • n sD .
By using the simplest upwind fluxes, we introduce the operator

divc D : R D × R 2 T → R D , 24 ⊗ ⊗ ⊗ D D n sD s x L * x K * x K x L Figure 8:
The definition of the gradient on the diamond mesh uses the unknowns stored at the nodes x K , x L , x K * , x L * . which depends on a density field given on D and a velocity field given on T , defined by

m D divc D (ρ D , u T ) = s∈∂D F s,D F s,D = m s (u s,D ) + ρ D -(u s,D ) -ρ D where u s,D = u K + u K * 2 • n sD for s = [x K , x K * ] ∈ ∂D,
with the standard notation x + = max(x, 0), x -= -min(x, 0). Going back to the definition of the previous section, we observe that

div D u T = 1 m D s∈∂D m s u s,D .
We deduce that the following consistency relation holds

div D u T = divc D (1 D , u T ). ( 12
)
The scheme for the mass conservation equation then reads

ρ n+1 D -ρ n D δt + divc D (ρ n D , u n T ) = 0. ( 13
)
Let us bring out two crucial properties of the scheme:

• We can identify a stability condition in order to preserve the positivity of the discrete density:

if ρ n D ≥ 0 and δt ≤ u T ∞ 1 m D s∈∂D m s -1 then ρ n+1 D ≥ 0.
• As a consequence of ( 12), we notice that divergence free velocity fields preserve homogeneous states: if

ρ n D ≡ 1 and div D u n T = 0 then ρ n+1 D ≡ 1.
Having at hand the discrete density piecewise constant on the diamond mesh D, we can naturally define the density on the primal and dual meshes by projection:

ρ n+1 K = 1 m K D∈D K m(D ∩ K)ρ n+1 D , ρ n+1 K * = 1 m K * D∈D K * m(D ∩ K * )ρ n+1 D ,
for all K ∈ M, K * ∈ M * . These quantities will be useful to deal with the momentum equation.

Treating the convection term in the momentum equation by taking into account the constraint

The viscous term in the momentum equation is handled by standard definition of the DDFV framework for diffusion operator, see Section 4.1; similarly, for the pressure term, we refer the reader to [START_REF] Krell | Stabilized DDFV schemes for Stokes problem with variable viscosity on general 2d meshes[END_REF][START_REF] Krell | Stabilized DDFV schemes for the incompressible Navier-Stokes equations[END_REF]. What is original in this approach is the treatment of the convection term

K div(ρu ⊗ u) = σ∈∂K σ (ρu • n)u.
It looks like "the transport of u by ρu": this intuition guides the upwinding strategy. At the discrete level, bearing in mind upwinding principles, we are led to introduce b T :

R D × (R 2 ) T × (R 2 ) T → (R 2 ) T b K (ρ D , v T , u T ) = 1 m K σ∈∂K F K,σ (ρ D , v T ) + u K -F K,σ (ρ D , v T ) -u L
with the momentum flux We adopt the same definition on the dual cells.

F K,σ = - m(D ∩ L) m D s∈∂D,s⊂K F s,D + m(D ∩ K) m D s∈∂D,s⊂L F s,D . ⊗ ⊗ F K,σ F s,D x L * x K * x L x K
The definition of the momentum flux F K,σ is based on the following requirement: from (13), we wish a similar conservation relation on the primal cells

K ∈ M m K ρ n+1 K -ρ n K δt + σ∈∂K (F K,σ ) n = 0.
(and on the dual cells K * ∈ M * as well). The construction relies on the definition of the Nagtegaal Finite Element [START_REF] Nagtegaal | On numerically accurate finite element solution in the fully plastic range[END_REF], see [START_REF] Goudon | DDFV scheme for incompressible Navier-Stokes equations with variable density[END_REF]. Eventually, the scheme for the momentum equation reads

Find u n+1 T ∈ E 0 and p n+1 D ∈ R D such that ρ n+1 T u n+1 T -ρ n T u n T δt + b T (ρ n D , u n T , u n+1 T ) +div T (-2η D D D u n+1 T + p n+1 D Id) = ρ n+1 T f n+1 T , div D u n+1 T = 0, D∈D m D p n+1 D = 0, Here E 0 ⊂ (R 2 )
T is the space of discrete velocities that accounts for the boundary conditions imposed on the velocity (typically the no-slip condition). As far as we are concerned with the numerical analysis of the scheme, we can draw the following conclusions:

• For incompressible flows, the scheme is energy-stable in the sense that it satisfies the following analog of the energy inequality

1 2δt ρ n+1 T u n+1 T 2 T - 1 2δt ρ n T u n T 2 T + 1 2δt ρ n T (u n+1 T -u n T ) 2 T +C η |||∇ D u n+1 T ||| 2 2 ≤ ρ n+1 T f n+1 T , u n+1 T T .
• The scheme is semi-implicit, and to update the velocity and pressure requires to solve a linear system: the invertibility of the system can be established and the solution (u n+1 T , p n+1

D

) is indeed uniquely defined.

• The approach does not use the fact that u is solenoidal and it can be adapted to any constraint that prescribes the divergence of the velocity field.

Remark 4.1. Note that the treatment of the momentum equation differs from the viewpoint adopted with the hybrid VF-FE scheme: in (9), we dealt with the discretization of the non-conservative form of the equation, involving ρ(∂ t u + u • ∇u), while the DDFV scheme approaches the conservative term ∂ t (ρu) + Div x (ρu ⊗ u).

Numerical validation

Before starting the discussion on numerical grounds, it is worth comparing the degrees of freedom of the two methods, for 2D simulations. We consider tessellations made of triangles, with unknowns stored at the edges E, the vertices V and the barycenters T according to the following table

Hybrid DDFV Velocity V & E V & T Pressure V E Density V E
Asymptotically we have #E

Analytical test

We first evaluate the abilities of the schemes to recover an analytical solution and we check the rate of convergence for the two methods. Let Ω be the unit disk of R 2 . We consider the function defined in polar coordinates (r, θ) by

       ρ ex (t, r, θ) = 2 + (2r -r 2 )e 0.1(r 2 -2r) cos(θ -sin t), u ex (t, r, θ) = -r sin θ cos t r cos θ cos t , p ex (t, r, θ) = sin(r cos θ) sin(r sin θ) sin(t).
It is an exact solution of the variable density incompressible Navier-Stokes system ( 7)-( 8) on Ω (with source term chosen accordingly). Compared with the analytical solution used in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF] and [START_REF] Guermond | A projection FEM for variable density incompressible flows[END_REF], it should be noted that this solution obeys now ∂ρ ex ∂r (t, r, θ) = 0 on ∂Ω.

This manufactured solution permits us to evaluate the efficiency of the method developed in Section 3.3 in order to take into account the homogeneous Neumann boundary condition.

We consider the disk Ω approximated by a polygonal domain Ω h and discretized with two isotropic unstructured meshes :

(A) either the mesh has some symmetry properties (Figure 10, left), (B) or the unstructured mesh is more general (Figure 10, right).

For this solution, we have u • ν = 0 on ∂Ω where ν is the unit outward normal to Ω, but u • ν = 0 on ∂Ω h , and the normal flux of the density should be evaluated on the interfaces that intersect the boundary ∂Ω h . For the results discussed below, with the hybrid FV-FE method the dual mesh is constructed as in Figure 12 (barycentric mesh, joining barycenters of the primal mesh to midpoints of the edges). For the DDFV scheme we use as dual mesh either the barycentric mesh in Figure 12 or the classical mesh in Figure 15 (joining directly the barycenters of the primal cells); the results are not substantially different.

We plot the maximum error in time evaluated in L 2 (Ω) norm in space versus the number of cells (a half of the number of triangles) on the density, velocity and pressure. The top of Figure 11 corresponds to simulations made using the mesh (A) and the bottom of Figure 11 to those made using the mesh (B). The errors are plotted for the DDFV scheme and for the two different strategies on the boundary when we use the hybrid FV-FE method. The lines corresponding to a rate of convergence of order one (slope 1) and order two (slope 2) are also displayed. The computations are performed until the final time T = 1. In the FV-FE scheme, we choose δt F E = h max and we compute δt F V ≤ δt F E in order to verify the CFL condition, see [START_REF] Calgaro | L ∞ stability of vertex-based MUSCL finite volume schemes on unstructured grids; simulation of incompressible flows with high density ratios[END_REF]Section 3.2]. In the DDFV scheme, the time step is similarly driven by the stability condition of the convection operator.

In particular, for the FV-FE scheme, we compare the treatment of the boundary conditions using the "ghost cells" or letting the MUSCL scheme degenerate at order one. For the density, the convergence rates (obtained by linear regression) are respectively of order 1.69 for cases (A) and 1.76 for case (B) if we adopt the strategy presented in Section 3.3. With the naive approach, the convergence rates decrease at order 1.42 for case (A) and 1.53 for case (B). Despite this loss of accuracy on the density, we observe that the convergence rates for velocity and pressure remain optimal (order 2, and the error curves for the degenerate and the ghost cells methods coincide). The degradation of the rates of convergence for the density, which are always smaller than 2, is also due to the topology of the unstructured meshes and the phenomenon is well known in the finite volume literature. On the one hand, the orders remains very satisfactory when we introduce the "ghost cells" and they are very similar to those obtained in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF] for an analytical solution with Dirichlet boundary conditions. One the other hand, we note that a loss of accuracy in the close vicinity of the boundary eventually leads to a degraded performance of the scheme, even if the order remains higher than one. The DDFV scheme is of order 1 for both the density and the velocity, and a surprisingly higher accuracy for the pressure. In fact we note that the rates are slightly better on coarse grids.

Rayleigh-Taylor instability

This Section is concerned with the numerical simulation of a Rayleigh-Taylor instability. This problem has been considered in numerous papers (see [START_REF] Tryggvason | Numerical simulations of the Rayleigh-Taylor instability[END_REF] for the inviscid case or [START_REF] Bell | A second-order projection method for variabledensity flows[END_REF][START_REF] Guermond | A projection FEM for variable density incompressible flows[END_REF][START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF][START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF] for viscous fluids). We consider the Incompressible Navier-Stokes system

∂ t ρ + div x (ρu) = 0, ∂ t (ρu) + Div x (ρu ⊗ u) + ∇ x p -µ∆ x u = ρg, div x (u) = 0, in the two-dimensional domain Ω = (-d/2, d/2) × (-2d, 2d),
where µ > 0 is the dynamical viscosity, supposed to be constant. The fluid, initially at rest, is subjected to gravity and g = G(0, -1) (G = 9.81 ms -2 is the gravitational acceleration and for the length of the domain we set d = 1 m). The domain splits into two regions with noticeably different densities, the heavier fluid being initially above the light fluid. The interface is slightly smoothed since we set at time t = 0

ρ 0 (x, y) = ρ m + ρ M 2 + ρ M -ρ m 2 tanh y -η cos(2πx/d) 0.01d ,
with ρ M > ρ m > 0, and η > 0 the amplitude of the initial perturbation. Due to the perturbation of the flat interface, the heavy fluid falls down, the light fluid moves up, with the formation of a typical mushroom shape. Let us rewrite the equation in dimensionless form. To this end we choose d as the length scale, while the time unit is defined by

T = d G .
For the velocity unit we set U = d T . The Reynolds number is then given by

Re = ρ m d 3/2 G 1/2 µ .
Up to an obvious change of notations, we arrive at

∂ t ρ + div x (ρu) = 0, ∂ t (ρu) + Div x (ρu ⊗ u) + ∇ x p - 1 Re ∆ x u = -ρ 0 1 , div x (u) = 0,
in the two-dimensional domain Ω = (-1/2, 1/2) × (-2, 2). The problem is supplemented by no-slip boundary conditions on the horizontal boundaries. Actually, the solution has symmetries and we compute the solution on the half domain (0, 1/2) × (-2, 2) with the following boundary conditions for the velocity field: denoting u = (u, v), we impose On the horizontal boundaries: u = 0, v = 0, On the vertical boundaries:

u = 0, ∂ x v = 0.
The symmetry implicitly induces the Neumann boundary condition for the density. The difficulty of the problem essentially depends on

• the Reynolds number,

• the density ratio between the light and the heavy fluid, which is measured by the so-called Atwood number

At = ρ M -ρ m ρ M + ρ m .
The larger Re and At, the more difficult the simulation. As it is reported in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF], instabilities of numerical nature might appear, and the simulation is highly sensitive to mesh orientation and treatment of the boundary terms.

We go back to this delicate issue in the results presented here, for both the hybrid FV-FE method and the DDFV scheme, comparing the performances of the two schemes on this benchmark. We define the parameters as follows:

Re = 1000, ρ m = 1, ρ M = 7 (so that At = 0.75).

For the initial data, the regularization parameter is η = 10 -2 . The primal mesh is made of 50 × 400 squares which have been cut in two triangles, as displayed in Figure 12. Convergence grid tests was performed in order to ensure that the mesh as well as the FE time step are fine enough to reach the grid convergence.

Primal mesh

Dual mesh

∂Ω

Figure 12: Barycentric mesh: the dual mesh is constructed by joining the barycenters of the primal mesh to the midpoints of the edges.

The hybrid FV-FE method runs with the MUSCL scheme described in Section 3, where we use a modified Van-Leer τ -limiter (see [START_REF] Calgaro | L ∞ stability of vertex-based MUSCL finite volume schemes on unstructured grids; simulation of incompressible flows with high density ratios[END_REF]) and starshaped control volumes of Figure 12. We compare the results obtained by the scheme which degenerates at order 1 in the vicinity of the boundary (see Figure 13) with the "ghost cell" approach described in Section 3.3 (see Figure 14). We observe that the degenerate scheme produces spurious instabilities which emerge at the right vertical boundary, and then spread in the domain as time increases. With the "ghost cells" technique we take into account the homogeneous Neumann boundary condition on ρ and the unphysical instabilities are smoothed out. In this case, the simulation is far more faithful to other results obtained in the literature [START_REF] Bell | A second-order projection method for variabledensity flows[END_REF][START_REF] Guermond | A projection FEM for variable density incompressible flows[END_REF][START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF][START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF]. This effect relies on the shape of the control volumes: the simulation in [14, see Figure 14, based on a 40 × 320 mesh] used square-shaped control volumes and was free of these instabilities. This is due to the fact that, in this case, the homogeneous Neumann boundary condition on ρ is implicitly taken into account, owing to the symmetries of the mesh, see Remark 3.1. We treat the same problem by using the DDFV scheme (which is only first order accurate). The primal mesh is the same as above. We compare results obtained by using

• either the dual mesh represented in Figure 12, • or the dual mesh displayed in Figure 15, where the control volumes are simply built by joining directly the barycenters of each triangle.

The latter dual mesh contains square or octagonal control volumes. At first sight the mesh in Figure 15 has much more structure than the barycentric mesh in Figure 12, and one would naively bet it produces more stable simulations. The results can be found in Figures 16 and 17 we observe the sensitivity to the mesh construction: the dual mesh Figure 15 produces spurious instabilities that develop on the right vertical boundary, where the light fluid moves to the top. The role of the mesh construction can be understood by performing the same simulations on the entire domain [0, 1]: the instabilities disappear (not reproduced here). This is due to the fact that the mesh in Figure 15, when working on the half-domain [0, 1/2], adopts a specific treatment of the boundary cells that breaks the natural symmetries of the problem. Conversely, the barycentric mesh preserves the symmetry of the domain. It is remarkable that an a priori harmless variation in the mesh construction produces such a sensitive effect. We also compare the results with a first order version of the hybrid FV-FE scheme, where the MUSCL scheme for the transport equation is replaced by a mere upwind scheme, see [START_REF] Brenner | Navier-Stokes revisited[END_REF] and [START_REF] Bresch | Effect of density dependent viscosities on multiphasic incompressible fluid models[END_REF]; results are given in Figure 18 (the mesh is still given as in Figure 12). The upwind scheme is more diffusive and the interface spreads over a larger number of cells. The smallest structures in the domain do not appear. In turn, the diffusion slows down the fall of the mushroom, which is clearly delayed compared to the other results. Using a more refined mesh (composed by 70 × 560 squares, for instance) the interface is only slightly accelerated. The DDFV is less subjected to numerical diffusion; it is also robust in reducing the apparition of spurious instabilities in the foot of the mushroom. The Hybrid-MUSCL scheme has a better resolution of the interfaces. The experiments on the Rayleigh-Taylor instabilities are in full agreement with the analytical test-case. Finally, while we use different methods to address the same problem, we remark that the simulation is highly sensitive to mesh effects and treatment of the boundary conditions. As reported elsewhere, we point out that results for the larger times of simulation should always be considered with caution since it becomes difficult to distinguish between physical and numerical instabilities.

Primal mesh

Dual mesh ∂Ω

Figure 1 :

 1 Figure 1: Evolution of the domain D 0 .

Figure 2 :

 2 Domain and initial data configurations.

Figure 3 :

 3 Figure 3: An adapted mesh that follows the avalanche front.

Figure 4 :

 4 Figure 4: Primal and dual meshes for the FV-FE method.

Figure 5 :

 5 Figure 5: Construction of the upstream and downstream gradients for the definition of multislope limiters.

Figure 6 :

 6 Figure 6: Definition of the boundary condition.

Figure 7 :

 7 Figure 7: Construction of the primal, dual and diamond meshes

Figure 9 :

 9 Figure 9: Definition of the momentum fluxes from the mass fluxes known on the edges of the diamond mesh

Figure 10 :

 10 Figure 10: Unstructured meshes for a disk Ω : symmetric one (left) and non-symmetric one (right).

Figure 11 :

 11 Figure 11: Rates of convergence for hybrid FV-FE and DDFV schemes on the symmetric unstructured mesh (top) and the non-symmetric mesh (bottom).

Figure 13 :

 13 Figure 13: Rayleigh-Taylor instability. Simulation by the Hybrid FV-FE method, MUSCL scheme that degenerates at the boundary. Density contours 2 ≤ ρ ≤ 4 at times T = 1, 1.5, 2, 2.5, 3, 3.25, 3.5, 3.75 and 4 (in Tryggvason's scale).

Figure 14 :

 14 Figure 14: Rayleigh-Taylor instability. Simulation by the Hybrid FV-FE method, MUSCL scheme with "ghost cells" for the boundary. Density contours 2 ≤ ρ ≤ 4 at times T = 1, 1.5, 2, 2.5, 3, 3.25, 3.5, 3.75 and 4 (in Tryggvason's scale).

Figure 15 :

 15 Figure 15: Classical mesh: the dual mesh is constructed by joining directly the barycenters of the primal mesh.

Figure 16 :

 16 Figure 16: Rayleigh-Taylor instability. Simulation by the DDFV method. The dual mesh is constructed as in Figure 12. Density contours 2 ≤ ρ ≤ 4 at times T = 1, 1.5, 2, 2.5, 3, 3.25, 3.5, 3.75 and 4 (in Tryggvason's scale).

Figure 17 :

 17 Figure 17: Rayleigh-Taylor instability. Simulation by the DDFV method. The dual mesh is constructed as in Figure 15. Density contours 2 ≤ ρ ≤ 4 at times T = 1, 1.5, 2, 2.5, 3, 3.25, 3.5, 3.75 and 4 (in Tryggvason's scale).

Figure 18 :

 18 Figure 18: Rayleigh-Taylor instability. Simulation by the Hybrid FV-FE method, with the upwind scheme for the transport. Density contours 2 ≤ ρ ≤ 4 at times T = 1, 1.5, 2, 2.5, 3, 3.25, 3.5, 3.75 and 4 (in Tryggvason's scale).

2 #T and #V1 2 #T . Since we deal with two components for the velocity we get DoF Hybrid 5#T, DoF DDFV 6#T.

Acknowledgments

This work was partly supported by the Labex CEMPI (ANR-11-LABX-0007-01) and by INRIA Lille Nord-Europe (team RAPSODI).