O Angelini

C Chavant

E Chénier
email: eric.chenier@univ-mlv.fr

R Eymard

S Granet

Finite volume approximation of a diffusion-dissolution model and application to nuclear waste storage 1

Keywords: Two phase Darcy flow in porous media, finite volume method, nuclear waste storage

The study of two phase flow in porous media under high capillary pressures, in the case where one phase is incompressible and the other phase is gaseous, shows complex phenomena. We present in this paper a numerical approximation method, based on a two pressures formulation in the case where both phases are miscible, which is shown to also handle the limit case of immiscible phases. The space discretization is performed using a finite volume method, which can handle general grids. The efficiency of the formulation is shown on three numerical examples related to underground waste disposal situations.

Introduction

A large community of scientists is concerned by understanding the mechanical and hydraulic behaviour of deep repository radioactive waste, in reason of its large impact on environment and human safety. This implies to be able to model and simulate complex phenomena such as the de-saturation and re-saturation of geological media, gas production induced by the corrosion of steel containers, within complex 3D heterogeneous and anisotropic domains including singular zones such as galleries and cells intersections. Moreover, materials with highly contrasted physical properties are involved in long time phenomena (from thousand to millions of years).

Hence the simulation of these physical features happens to be a complex task, and their validation is a major concern for the safety improvement of the industrial devices. Computational benchmarks, such as the Couplex Gaz benchmark [START_REF] Talandier | Synthèse du benchmark couplex-gaz[END_REF], are useful for the definition of relevant physical models and numerical methods. Indeed, the Couplex Gaz benchmark has shown that the Darcy flow of two immiscible phases, the first one being an incompressible liquid phase and the second one the gaseous phase, can lead, in presence of high capillary pressures, to unphysical situations and to drastic numerical difficulties. This model implies the displacement of a free boundary between zones where the two fluid phases are simultaneously present (called in this paper the "under-saturated" zone) and zones where the only present fluid phase is the liquid phase (called in this paper the "saturated" zone). In an unexpected way, it can lead to the existence of zones where the liquid phase has been removed because of capillary forces and cannot be instantaneously replaced by the gaseous phase, hence creating a vacuum volume. It is interesting to exhibit the origin of such behaviour from the immiscible model equations, that we now describe. Throughout all this paper, we shall denote by l and g the liquid and gaseous phases, and by w and h the water and gas components (the choice of h for the gas component is linked with the fact that this model should apply to hydrogen generation due to acid attack of metallic containers, which may arise during long term storage). In the immiscible model, the liquid phase only contains the water component and the gaseous phase only contains the gas component. This leads to the following conservation equations:

                   ∂(φ(x)ρ l S l )
∂t div k(x) ρ l k rl (S l ) µ l (∇P l -ρ l g) = 0, ∂(φ(x)ρ g (P g )(1 -S l )) ∂t div k(x) ρ g (P g )k rg (S l ) µ g (∇P g -ρ g (P g )g) = 0, P g -P l = P c (S l ) and ρ g (P g ) = M h P g /RT.

(

In [START_REF] Abadpour | Method of negative saturations for modeling two-phase compositional flow with oversaturated zones[END_REF], the density of the liquid phase ρ l is assumed to be constant, that of the gaseous phase ρ g (P g ) is given by the ideal gas law, the phase viscosities µ l and µ g are constant, the relative permeabilities k rl (S l ), k rg (S l ) and capillary pressure P c (S l ) functions are given for example by Van Genuchten laws. In the Couplex Gaz benchmark [START_REF] Talandier | Synthèse du benchmark couplex-gaz[END_REF], some materials happen to have a very low permeability ( k(x) varies between 10 -14 and 10 -21 m 2 ) and the intensity of the capillary pressures is very high, about 10 7 Pa. The porosity φ(x) varies within the range [0.05, 0.3]. A consequence of such data is that the initial gas pressure remains negligible (with an order of magnitude about 10 5 Pa), compared to the liquid pressure, whose order of magnitude (which is negative, as it is classical for such high capillary pressures) is about -10 7 Pa. Hence, to a first approximation, P g -P l -P l , and the liquid saturation can be approximately solved by the relation

S l = C(-P l ),
where the capacity function C(P) is, as in the case of Richards equation, the reciprocal function of the capillary pressure for P ≥ 0 and equal to 1 for P ≤ 0 (see Figure 1). Substituting this expression into the liquid component conservation equation, we get that P l is solution of the following equation, which appears to be (again, to a first approximation) decoupled from the gas conservation equation:

∂(φ(x)ρ l C(-P l )) ∂t -div k(x) ρ l k rl (C(-P l )) µ l (∇P l -ρ l g) = 0. (2)
Then Equation ( 2) is the classical Richards equation (which is of elliptic/parabolic type). The solution of this wellposed problem, denoted by P l (x, t), allows for expressing the liquid saturation by S l (x, t) = C(-P l (x, t)). Let us now consider the gas conservation equation, divided by M h /RT , in which we now assume that the liquid saturation is the imposed function given by the previous expression:

∂(φ(x)P g (1 -S l (x, t))) ∂t -div k(x) k rg (S l (x, t)) µ g ( 1 2 ∇P 2 g - M h P 2 g RT g) = 0. (3)
Then Equation (3), whose unknown is the function P g , behaves in the region S l (x, t) < 1 as the solution of the degenerate parabolic equation ∂u ∂t -div(∇u m ) = 0 with m = 2. Such an equation, sometimes called the "porous medium equation" (in reference to a fractured medium filled by a fluid: u is the fracture width, assumed to be proportional to the pressure of the fluid, whereas the permeability behaves as an increasing function of the fracture width, say mu m-1 ) is degenerate in the sense that the value u = 0 propagates with a finite velocity. This phenomenon arises for the solution P g of Equation (3) (this is shown in numerical simulations provided at the end of this paper). As a consequence, the gas phase cannot instantaneously fill the volume released by the removal of the water phase, since the value P g = 0 can only enter into this domain with a finite velocity. Let us emphasise that such a behaviour of the solution of Equations ( 1) does not occur in the case where ρ g is assumed to be constant, instead of ρ g = M h P g /RT . It is no longer observed in case of low capillary pressure.

These mathematical features introduce some complexity in the numerical procedures. For example, the numerical formulation, issued from the oil reservoir engineering framework, consists in selecting P l and S l as primary unknowns for the approximate problem. This formulation is shown to be very efficient in the case of the drainage of a vertical sand column, despite of the fact that the diffusion, due to the capillary pressure, vanishes at the moving boundary between the under-saturated and the saturated zones [START_REF] Liakopoulos | Transient flow through unsaturated porous media[END_REF]. Unfortunately, this formulation does not succeed in presence of the non-physical situation presented above, that is the apparition of vacuum zones. In order to overcome this difficulty, it is natural to modify Equations (1) in order to get rid of this nonphysical situation. Two new physical features have then to be introduced. Firstly, it is natural to let the water component, constituting the liquid phase pass to the gaseous state, hence providing strictly positive gas pressures. Secondly, the possible dissolution of the gas component into the liquid phase has to be considered. Hence, in this paper, we now consider this miscible model, where each phase can contain both components. Although this more complete model avoids the occurrence of the nonphysical situation precised above, the numerical approximation of this model gives rise to new difficulties, since the composition of both phases has to be prescribed by the model. The determination of these compositions is a consequence of equilibrium laws, which only apply in the under-saturated zone. Unfortunately, if one assumes that, initially, some part of the porous domain only contains pure water, the gas component conservation equation degenerates in this zone. Then the choice of the variables is essential in order to avoid divisions by zero. After elimination of the local equations, there remain only two non-local conservation equations, with two unknown functions of x, t, called the primary unknowns. Different choices for the primary unknowns are available in the literature, some of them being different in the saturated and the under-saturated zones. The choice P l , S l , analysed in [START_REF] Panfilova | Ecoulement diphasiques en milieu poreux : modèle de ménisque[END_REF][START_REF] Gagneux | Analyse mathématique de modèles non linéaires de l'ingénierie pétrolière[END_REF], does not allow for solving the gas conservation equation in the saturated zone. A solution to this problem can be, as done in [START_REF] Bourgeat | Two-phase, partially miscible flow and transport modeling in porous media; application to gas migration in a nuclear waste repository[END_REF], to choose the liquid pressure, P l and the total gaseous mass density X = ρ h g S g + ρ h l S l as main variables. Focusing on the mathematical treatment of the capillary pressure curve, we observe that, in the case of Van Genuchten-Mualem laws, the derivative of the capillary pressure with respect to the saturation tends to infinity as the saturation tends to 1; the change of variable S l = S l (P l , X) does not allow to get rid of this singularity.

Another possibility, presented in [START_REF] Abadpour | Method of negative saturations for modeling two-phase compositional flow with oversaturated zones[END_REF], is to extend the notion of phase saturation, allowing for negative values, and values greater than one. In another method, developed in [START_REF] Wua | On the selection of primary variables in numerical formulation for modeling multiphase flow in porous media[END_REF], primary variables are changed, depending on the value of the saturation: they are switched when, during a time step, the gas phase appears in the saturated zone. These methods remain sensitive to capillary pressure singularities.

This paper is devoted to a new numerical procedure for approximating the solution of liquid-gas two phase miscible flow in porous media problem. This method is based on the choice of the two pressures, P l and P g , as primary unknowns, whatever be the value of the saturation. We show in the paper that it leads to remarkably smooth and stable numerical approximations. Note also that such a choice is entirely compatible with the choice of finite element techniques for the space discretization, which are also shown to be efficient for the poromechanical coupling between porous medium flow and soil deformation [START_REF] Chavant | Thermo-hydro-mechanical numerical modeling: Application to a geological nuclear waste disposal[END_REF][START_REF] Gerard | Numerical modelling of coupled mechanics and gas transfer around radioactive waste in long-term storage[END_REF]. On the other hand, finite volume methods on structured grids have been used for decades in the oil engineering framework, where multi-phase multi-component flows in porous media have to be computed. This is justified by the fact that finite volume methods are natural in the framework of nonlinear hyperbolic equations. Indeed, such a framework arises in the case where the diffusion due to the capillary pressure is not large enough when, during oil production, the displacement of the fluids is mainly due to the fact that the fluids are injected and produced. In this case, upstream weighting techniques are necessary for the stability of the numerical schemes. Surprisingly, in the physical cases considered in this paper, such upstream weighting numerical techniques remain necessary due to vanishing diffusion effects within the saturated zone [START_REF] Coussy | A vanishing diffusion process in unsaturated soils[END_REF]. In order to obtain a method which can gather the advantages of the finite element and finite volume methods, we propose here a finite volume formulation suited to unstructured general meshes. This method is based on a technique developed for heterogeneous anisotropic diffusion equation [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes. SUSHI: a scheme using stabilisation and hybrid interfaces[END_REF].

This paper is structured as follows. In section 2, we present the model for diffusion-dissolution liquid-gas flow in a porous medium. Then, in section 3, we present the discretization method. Finally, in section 4, we show the efficiency and the accuracy of the formulation, on typical test case used in order to simulate the storage of nuclear waste. In particular, we consider the case where it is used as a regularisation method for the model (1).

Physical and mathematical model

We now present a standard model for miscible two phase flow in a porous medium. The density ρ c p of each component c = w, h in each phase p = l, g is linked to the density of the phase p = l, g by the relation:

ρ p = ρ h p + ρ w p . (4)
The molar concentration C p of phase p = l, g, the molar concentration C c p and the molar fraction X c p of component c = w, h in phase p = l, g are respectively defined by

C c p = ρ c p M c , C p = C h p + C w p and X c p = C c p C p , (5)
where M c denotes the molar mass of component c. As done in the introduction of this paper, the liquid saturation S l is chosen as variable for the relative permeability k rp and capillary pressure P c curves, whereas the gas saturation S g satisfies

S l + S g = 1. (6)
The liquid phase and gas phase pressures are linked with S l by the relation

P g -P l = P c (S l ). (7)
The mass conservation equations of the two components c = w, h read:

               ∂m w l ∂t + ∂m w g ∂t + div ρ w l F l ρ l + ρ w g F g ρ g + J w l + J w g = 0 ∂m h l ∂t + ∂m h g ∂t + div ρ h l F l ρ l + ρ h g F g ρ g + J h l + J h g = 0, (8)
with m c p = φρ c p S p . In ( 8), we denote by F p the flux given by Darcy's law:

F p ρ p = - kk rp µ p (∇P p -ρ p g), (9)
where we denote by g the gravity acceleration. We define J c p as the diffusive flux given by Fick's law [START_REF] Coussy | [END_REF][START_REF] Mainguy | Modèles de diffusion non-linéaires en milieux poreux. applications à la dissolution et au séchage des matériaux cimentaires[END_REF]:

J w g = -φM w S g D w g C g ∇X w g , J h g = -φM h S g D h g C g ∇X h g , J w l = -φM w S l D w l C l ∇X w l , J h l = -φM h S l D h l C l ∇X h l .
Since, in the applications under consideration in this paper, Fick's flow of water in the liquid phase is negligible in front of the liquid phase Darcy flow, we set in (8):

J w g = -φM w S g D w g C g ∇X w g , J h g = -φM h S g D h g C g ∇X h g , J w l = 0, J h l = -φM h S l D h l C l ∇X h l . (10)
In order to close the system, equilibrium laws for components c = w, h are required. Since our concern is to focus on gas injection or on the regularisation of immiscible two-phase flow, we now assume that ρ w g = 0, and we are no longer accounting for the water component within the gaseous phase. We assume that the equilibrium of the gas component h between the two phases is prescribed by Henry's law, which can be expressed by:

ρ h l = HM h P g and S l < 1 or ρ h l ≤ HM h P l and S l = 1 , (11)
denoting by H Henry's constant at the temperature of the domain (supposed to be constant in space and time) and assuming that P c (1) = 0.

Let us now rewrite the system of all preceding equations, remarking that, in the saturated zone S l = 1, all terms of System (8) which depend on P g are multiplied by S g = 0 or k rg (1) = 0. Hence P g , which is linked to ρ h l in the region S l < 1 by [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes. SUSHI: a scheme using stabilisation and hybrid interfaces[END_REF], can be freely defined in the region S l = 1. Hence we choose to set

P g = ρ h l HM h in both cases S l = 1, and S l < 1. ( 12
)
Then, from [START_REF] Gagneux | Analyse mathématique de modèles non linéaires de l'ingénierie pétrolière[END_REF], the relation C h l = HP g holds in both saturated and under-saturated zones and from [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes. SUSHI: a scheme using stabilisation and hybrid interfaces[END_REF], the inequality P g ≤ P l holds in the saturated zone S l = 1. Therefore, P g is continuous across the transition area between the saturated and the under-saturated zones, and we can write the relation

S l = C(P g -P l ), (13)
where, as in the introduction of this paper, the capacity function C(P) is defined as the reciprocal function to the capillary pressure for P ≥ 0 and equal to 1 for P ≤ 0. This function C(P) is then in some cases differentiable at point P = 0 (this holds for example in the case of Van Genuchten's functions, see Figure 1). Hence the system of all equations, from (4) to [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes. SUSHI: a scheme using stabilisation and hybrid interfaces[END_REF], resumes to the following system of equations, with respect to the two unknowns P l and P g :

             ∂φA w (P l , P g ) ∂t -div kB w l (P l , P g )(∇P l -ρ l g) = 0
∂φA h (P l , P g ) ∂t div kB h l (P l , P g )(∇P l -ρ l g) + kB h g (P l , P g )(∇P g -ρ g (P g )g) + φD h l (P l , P g )∇P g = 0,

where A w (P l , P g ) = ρ l C(P g -P l ), B w l (P l , P g ) = ρ l k rl (C(P g -P l )) µ l , ρ g (P g ) = M h P g RT , A h (P l , P g ) = ρ g (P g )(1 -C(P g -P l )) + HM h P g C(P g -P l ),

B h l (P l , P g ) = HM h P g k rl (C(P g -P l )) µ l , B h g (P l , P g ) = ρ g (P g ) k rg (C(P g -P l )) µ g , D h l (P l , P g ) = M h C(P g -P l )D h l ρ w l M w H. (15)
An advantage of this formulation, from the point of view of numerical implementation, is that the coefficient ∂A h (P l , P g )/∂P g does no longer vanish in the saturated zone, where it is equal to HM h . If we assume that D h l > 0, we also get that the coefficient of ∇P g does not vanish in the saturated zone, since it is equal to D h l HM h > 0. Nevertheless, if we approximate the immiscible problem (1) by System ( 14), using a small value for H and letting D h l = 0, the coefficient of ∇P g vanishes in the saturated zone (this is the approximation method used in the third example in section 4).

The numerical scheme

Our aim is to extend to the system (14), the SUSHI scheme (Scheme Using Stabilisation and Hybrid Interfaces), presented in [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes. SUSHI: a scheme using stabilisation and hybrid interfaces[END_REF][START_REF] Eymard | A new finite volume scheme for anisotropic diffusion problems on general grids: convergence analysis[END_REF][START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes. SUSHI: a scheme using stabilisation and hybrid interfaces[END_REF][START_REF] Eymard | A cell-centered finite-volume approximation for anisotropic diffusion operators on unstructured meshes in any space dimension[END_REF] in the case of a heterogeneous and anisotropic pure diffusion problem. This scheme leads to the approximation, for a scalar field u defined on the domain, of the fluxσ k∇u.n K,σ dγ, in the case where the domain Ω has been partitioned into a mesh M (for example, triangles, quadrangles,. . . ), K is one element of the mesh M and σ is an edge of K (see figure 2). The continuous field u is approximated by discrete values u K , approximating u at given points x K ∈ K (in the case of triangles, we choose here the barycentre of the triangle, but this is not mandatory, and we could as well choose the circumcentre for example), and u σ , approximating u at the barycentre x σ of σ.

x K n K,σ |σ| d K,σ x σ K σ K Figure 2: Geometrical description of control volume K
First, an approximate gradient of u in K is reconstructed, thanks to the relation

∇ K U K = 1 |K| σ∈E K |σ|(u σ -u K )n K,σ ,
where

U K = (u K , (u σ ) σ∈E K ).
Then an approximate gradient is reconstructed in the cone K σ with vertex x K and basis σ,

∇ K,σ U K = ∇ K U K + β K d K,σ (u σ -u K -∇ K U K • (x σ -x K )) n K,σ ,
where β K is a parameter which is equal to √ d for getting back the 5-point scheme in the case of rectangular meshes. In case of highly distorted meshes, this parameter can be adjusted in order to prevent from oscillations [START_REF] Angelini | A finite volume scheme for diffusion problems on general meshes applying monotony constraints[END_REF]. We assume that the permeability k takes the constant value k K in K. Then (F K,σ (U K )) σ∈E K is defined as the unique family of values such that the following relation holds

σ∈E K F K,σ (U K )(v σ -v K ) = σ∈E K |K σ | k K ∇ K,σ U K • ∇ K,σ V K , ∀V K = (v K , (v σ ) σ∈E K ) ∈ R 1+#E K . ( 16
)
Note that, in [START_REF] Panfilova | Ecoulement diphasiques en milieu poreux : modèle de ménisque[END_REF], the value |K σ |, measure of the cone K σ , is given by

|σ|d K,σ
d , where d is the space dimension. It results from [START_REF] Panfilova | Ecoulement diphasiques en milieu poreux : modèle de ménisque[END_REF] an expression of F K,σ under the form:

F K,σ (U K ) = σ ∈E K M σ,σ K (u K -u σ ) (17)
where the local matrices (M σ,σ K ) σ,σ ∈E K are symmetric and positive. In the case of a monophasic diffusion problem, the balance of the diffusion terms in the control volume K is equal to σ∈E K F K,σ (U K ), and the continuity equation F K,σ (U K ) + F L,σ (U L ) = 0 is requested for all common interface σ to two neighbouring control volumes K and L. Since the scheme ( 14) also contains pure diffusion terms, we also need to have an approximate ofσ ∇u.n K,σ dγ (this term no longer includes the permeability matrix k). Hence we similarly define the family (G K,σ (U K )) σ∈E K as the unique family of values such that

σ∈E K G K,σ (U K )(v σ -v K ) = σ∈E K |K σ | ∇ K,σ U K • ∇ K,σ V K , ∀V K = (v K , (v σ ) σ∈E K ) ∈ R 1+#E K , (18)
which provides, for the expression of G K,σ (U K ), an expression similar to [START_REF] Smaï | Apparition-disparition de phase dans un écoulement diphasique eau/hydrogène en milieu poreux : Injection de gaz dans un milieu saturé en eau pure[END_REF].

We may now use these functions F K,σ and G K,σ for defining the approximation of ( 14). We first use, for the time discretization, an Euler backward (implicit) scheme, using the time step ∆t > 0. The discrete unknowns are now the values P l,K , P g,K approximating at the end of the time step the liquid and gas pressures at point x K , and the values P l,σ , P g,σ approximating at the end of the time step the liquid and gas pressures at point x σ . The approximation of the term

F w l,K,σ = - σ kB w l (P l , P g )(∇P l -ρ l g).n K,σ dγ
is then given, in the case where σ is a common interface between two neighbouring control volumes K and L, by

F w l,K,σ = B w l (P l , P g ) ups K,L (F K,σ (P l,K , (P l,σ ) σ∈E K ) -|σ|ρ l k K g • n K,σ ), (19)
where the upstream weighting expression B w l (P l , P g ) ups K,L is defined by

B w l (P l , P g ) ups K,L = B w l (P l,K , P g,K ) if (F K,σ (P l,K , (P l,σ ) σ∈E K ) -k K |σ|ρ l g • n K,σ ) ≥ 0 B w l (P l , P g ) ups K,L = B w l (P l,L , P g,L )
otherwise.

The continuity equation

F K,σ (P l,K , (P l,σ ) σ∈E K ) -|σ|ρ l k K g • n K,σ + F L,σ (P l,L , (P l,σ ) σ∈E L ) -|σ|ρ l k L g • n L,σ = 0, (20)
is also solved. It could be possible to consider different expressions for the continuity equation (20), involving for example nonlinear terms at the interface. The advantage of ( 20) is that it remains linear and does not degenerate if a phase disappear. Note also that a consequence of (20) and of the definition ( 19) of F w l,K,σ is the continuity equation F w l,K,σ + F w l,L,σ = 0. Similar expressions are derived for the terms involved in the conservation of the gas component. The upstream weighting definition for B h l (P l , P g ) is again depending on the sign of the expression F K,σ (P l,K , (P l,σ

) σ∈E K ) -k K |σ|ρ l g • n K,σ
, whereas the upstream weighting definition for B h g (P l , P g ) is depending on the sign of the expression F K,σ (P g,K , (P g,σ

) σ∈E K ) -k K |σ|ρ g (P g,K )g • n K,σ . Finally, the diffusion term D h l,K,σ = - σ φD h l (P l , P g )∇P g .n K,σ dγ is approximated by D h l,K,σ = 1 2 (φ K D h l (P l,K , P g,K )G K,σ (P g,K , (P g,σ ) σ∈E K ) -φ L D h l (P l,L , P g,L )G L,σ (P g,L , (P g,σ ) σ∈E L )) .
Note that all terms A c (P l , P g ), B c p (P l , P g ), D h l (P l , P g ) involve some functions which depend on the rock type (capacity function, relative permeabilities). Hence a dependence with respect to the control volume has not been written for the sake of simpler notations. The scheme can then be written, using the superscript "-" for expressing values at the beginning of the time step:

|K|φ K ∆t (A w (P l,K , P g,K ) -A w (P - l,K , P - g,K )) + σ∈E K F w l,K,σ = 0, ∀K ∈ M |K|φ K ∆t (A h (P l,K , P g,K ) -A h (P - l,K , P - g,K )) + σ∈E K (F h l,K,σ + F h g,K,σ + D h l,K,σ ) = 0, ∀K ∈ M F K,σ (P l,K , (P l,σ ) σ∈E K ) -|σ|ρ l k K g • n K,σ + F L,σ (P l,L , (P l,σ ) σ∈E L ) -|σ|ρ l k L g • n L,σ = 0, ∀σ ∈ E int F K,σ (P g,K , (P g,σ ) σ∈E K ) -|σ|ρ g (P g,K )k K g • n K,σ + F L,σ (P g,L , (P g,σ ) σ∈E L ) -|σ|ρ g (P g,L )k L g • n L,σ = 0, ∀σ ∈ E int . (21)
Hence we get a system of nonlinear discrete equations. We approximate its solution, using the exact Newton's method (taking advantage of the fact that all functions are given by an analytical expression in the numerical examples considered in this paper). Since our aim is to focus on the schemes, we only used direct Gaussian elimination solvers. The convergence criterion is tuned such that an average number of 2 or 3 iterations is performed at each time step. This numerical scheme has been implemented within Code_Aster [8], which is an open source finite element code used for research/development purposes in EDF.

Numerical results

We present in this paper three applications of the use of Scheme (21). The first one is a benchmark problem proposed by F. Smaï et al. [START_REF] Smaï | Apparition-disparition de phase dans un écoulement diphasique eau/hydrogène en milieu poreux : Injection de gaz dans un milieu saturé en eau pure[END_REF]; the second one is a two-dimensional extension of this benchmark problem. The third one corresponds to the use of Scheme (21) for regularising the immiscible problem described in the introduction of this paper.

Test 1: Gas injection in a saturated porous media in one dimension

The aim of the test proposed in [START_REF] Smaï | Apparition-disparition de phase dans un écoulement diphasique eau/hydrogène en milieu poreux : Injection de gaz dans un milieu saturé en eau pure[END_REF] is to check the numerical schemes in the case of production of hydrogen due to acid attack of metallic containers. The production of hydrogen is represented by a mass injection rate Q during the first 500 000 years, equal to 0 after, at the left boundary of a one dimension domain Ω =]0m; 200m[ (the section of the domain is assumed to be equal to 1 m 2 ), initially saturated with pure water. We consider two cases:

1. the first one is that provided in [START_REF] Smaï | Apparition-disparition de phase dans un écoulement diphasique eau/hydrogène en milieu poreux : Injection de gaz dans un milieu saturé en eau pure[END_REF], where Q = 1.766 10 -13 kg.m -2 .s -1 = 5.57 10 -6 kg.m -2 .year -1 , 2. the second case, motivated by additional interesting physical phenomena, is given by Q = 4.76 10 -13 kg.m -2 .s -1 = 1.50 10 -5 kg.m -2 .year -1 .

The water flux is assumed to be null at the left boundary. At the right boundary, the value 10 6 Pa is imposed for the liquid pressure and it is assumed that there is no free nor dissolved gas.

For this test we use the Van Genuchten-Mualem model, defined by the values of the parameters P r , S lr , S gr , m and n, for expressing the capillary pressure and the relative permeability:

S le = S l -S lr 1 -S lr -S gr , P c = P r (S -1 m le -1) 1 n , k rl = S le (1 -(1 -S 1 m le ) m ) 2 , k rg = (1 -S le )(1 -S 1 m le ) 2m . ( 22
)
The parameters involved by ( 22), as well as the physical data and fluid characteristics are resumed in Table 1. The time variable t belongs to the range ]0, 10 6 [ years, discretized in 556 time steps whose size is between 1 day at the beginning of the simulation and 16 666 years at its end. Figure 3 presents the profiles of the gas pressure at different times, in the two cases. We observe that, in Case 1, the gas pressure mainly behaves as expected, in the sense that, at each location, it increases with respect to the time, except for the small values of x at t = 500 000 years, where this pressure is smaller than that obtained for t = 100 000 years. This behaviour is more clearly observed in Case 2: the gas pressure is lower for t = 500 000 years than that obtained for t = 100 000 years, in almost the whole region S g > 0. This is due to the fact that the main part of the volume needed for the gas phase is occupied at t = 100 000 years. The increase of this volume imposes an outward water rate at the right boundary of the domain, and therefore a strictly negative gradient of the liquid pressure along the domain. When the total volume of the gaseous phase has reached its maximum value, there is no further need for a large liquid flow from the left to the right of the domain. Then no negative gradient for the liquid pressure is necessary, and the liquid pressure tends to the constant value imposed at the right boundary, whereas a small gas pressure gradient remains necessary in the zone S l < 1 for transporting the injected gas to the right boundary. We present in Figure 4 the values of the pressures of both phases at the left boundary of the domain, with respect to the time variable, only for Case 1 (the results for Case 2 show the same behavior). After the end of the gas injection (i.e., for t > 500 000 years), the gas component is progressively eliminated thanks to the value 0 imposed at the right boundary. Hence the gas phase disappears in the domain, which imposes a flux of water, from the right to the left, in order to replace the volume occupied by the gas phase. This induces a positive gradient to the liquid pressure. Hence the liquid pressure must become smaller at the left boundary than that imposed at the right boundary, in order to get a positive gradient. This behaviour finishes when the liquid saturation is again equal to 1 in all the domain (in Figure 4, the indications S = 1 or S < 1 respectively mean that S l (0, t) = 1 or S l (0, t) < 1). We recall that the values of the gas pressure, which are lower than that of the liquid pressure must be interpreted as a gas concentration in the liquid phase.

Porous medium

Figure 5 presents the liquid saturation obtained at different times, only for Case 1 (the results for Case 2 are again similar), computed using the capacity function from the pressure values of both phases. We observe that the location of the free boundary is an easy consequence of the discrete conservation equations. Figure 6 presents the outward rates (liquid and gas) at the right boundary, with respect to time. The inversion of the sense of the water flow after the end of injection is clearly put in evidence.

Test 2: Gas injection in a saturated porous media in two dimensions

This test is a two dimensional adaptation of the benchmark problem presented in section 4.1. We assume that the domain is horizontal with thickness 1 m (Figure 7). The porous medium and fluid data are again given by ( 22) and Table 1. The domain, depicted in left part of Figure 7, is initially saturated with pure water. The value 10 6 Pa is imposed to the liquid pressure in the upper right corner of the domain, which is assumed to be free of gas component. The total gas rate, injected in the lower left corner of the domain, is equal to 1.76 10 -11 kg.m -2 .s -1 during the first 500 000 years. The other boundaries are assumed to be impervious. We use a triangular mesh with 1632 elements (right part of Figure 7). The time interval is discretized in 466 time steps, whose size varies between 1 day and 16 666 years. Figure 8 shows the gas pressure and saturation at different times into the domain. It shows a similar behaviour in 2D to that observed in 1D. The red colour in the saturation field shows the saturated zone, which corresponds to a gas pressure lower than the liquid pressure. The boundary between the saturated and the under-saturated zone exactly corresponds to the isovalue line P g = 10 6 Pa. This test confirms that the numerical scheme is suited for handling this problem in 2D (further tests will be performed in 3D).

Test 3: approximation of immiscible two-phase flow in one dimension

This test is firstly dedicated to numerically validate the analysis of the behaviour of the solution, proposed in the introduction of this paper, and secondly to show that the numerical method presented in this paper leads to an easy numerical approximation of the problem. The data are those provided by [START_REF] Chavant | Cas test diphasique[END_REF]. We consider the 1D domain ]0, 1[ m (assumed to have a section equal to 1 m 2 ). It is composed of two materials, M1 and M2, as shown in Figure 9.

The capillary pressure and the relative permeabilities are defined by

S l (P c ) = 1 + P c A 1/(1-B) -B , k rl (S l ) = 1 + a(S -b l -1) c -d , k rg (S l ) = (1 -S l ) 2 (1 -S 5/3 l ), (23)
where the parameters A, B, a, b, c, d and the porous medium data are provided by Table 2. We assume no flow boundary conditions at the left and right parts of the domain. We consider the following two cases for the initial conditions:

Case 1: S l = 0.77, P g = 10 5 Pa and P l = P g -P c (S l ) in M1, S l = 1, P g = 0 Pa and P l = 10 5 Pa in M2.

Case 2: S l = 0.77, P g = 10 5 Pa and P l = P g -P c (S l ) in M1, S l = 1, P g = 10 5 Pa and P l = 10 5 Pa in M2.

A 1D mesh with 200 elements has been used, and the time period [0, 10 11 ] s (3169 years) is considered. A number of 3730 time steps is required, for both initial data cases, whereas the size of the time steps varies from .1 s to 5. 10 6 s (arbitrary above limit) during the computation. The gas pressure is shown in Figure 10 for Case 1 and in Figure 11 for Case 2. Since, in Case 2, the amount of dissolved gas in water is very small accounting for the value of H provided by Table 2, one could expect that both cases would lead to close results. This is indeed what we observed, since, for the other variables (liquid pressure and saturation, total mass of gas per unit of porous medium volume), both cases lead to identical results up to the plotting precision due to the fact that the amount of dissolved gas is very small. So only one figure is proposed for the liquid pressure (Figure 12), the liquid saturation (Figure 13) and the global density of gas (defined by the quantity m h l + m h g of the balance equation ( 8)) is plotted in Figure 14. We observe that the numerical results follow the analysis developped in the introduction of this paper. In a first period (t < 10 5 s), water flows from region M2 to region M1 (see Figure 13). This leads to a reduction of the porous volume occupied by the gas phase in M1, and provides an increase in the value of the gas pressure (Figures 10 and11). The volume which was previously occupied by the water phase in M2 cannot be instantaneously filled by the gas phase present in M1. In case 1 (Figure 11), an exact vacuum is observed in M2, whereas in case 2, the main part of the dissolved gas pass to the gas phase. But, since in case 2, the amound of dissolved gas is very small, the pressures obtained in the gas phase are very low and close to 0 in the case where 1-S l is large enough. For 1-S l → 0 as x → 1, we naturally observe in Figure 11 that P g → 10 5 , which is equal to the initial equilibrium pressure of the dissolved gas. Let us notice that the location in M2 of the boundary between the regions S l < 1 to S l = 1 is more precisely assessed in Figure 11 by the location at which the gas pressure tends to its initial value, about 0.6 m for t = 100 s and 0.95s for t = 1000 s, than by Figure 13. We observe that this location moves with a finite velocity, which is expected for the solution of (2), similar to Richards equation.

At time t = 20000 s, the gas pressures obtained in both cases are still different but already very close. For larger times, the pressures obtained in both cases become undistinguishable. We then observe in M2, for t = 10 6 , 10 7 and 10 8 s, the characteristic shape corresponding to the porous medium equation ∂u/∂t -∂ 2 (u m )/∂x 2 = 0 with m = 2, that is a finite slope at u = 0. At large times, the gas pressure becomes constant, which is expected since the gas phase is mobile in the whole domain. Since the volume occupied by the gas phase at the end of the simulation in M1 and M2 is the same as the one which was initially occupied by the gas phase in M1, the final pressure is about 10 5 Pa, which is the same as the initial one (again, even in case 2, the amount of dissolved gas does not significantly modify this quantity).

Let us now comment the results obtained for the liquid pressure in Figure 12 and the liquid saturation in Figure 13. This pressure is highly negative, which corresponds to traction cases observed for materials such as clays and cement concrete. We must recall that the pressures are linked to the free energy of the water phase, which accounts for all types of interactions with the solid phase (such interactions being particularly complex in the case of clays), and Darcy's law provides a dissipation mechanism compatible with the second principle of thermodynamics [START_REF] Coussy | [END_REF]. We observe that, at short times, the variation of the saturation is greater in region M2 than in region M1, due to the contrast between the porosity values. At all times, the saturation presents a discontinuity in x = 0.5, due to the fact that, since both phases flow at this location, the gas and the liquid pressures must remain continuous. Therefore, the saturation must respect the equation P (1) c (S (1) l (0.5, t)) = P (2) c (S (2) l (0.5, t)), where upper indices (1) and ( 2) denote the values and functions respectively available in regions M1 and M2, at all times t > 0, which leads to different left and right limits of the saturation in x = 0.5. At large times, the liquid pressure becomes constant, which is expected since the water phase is mobile in the whole domain, and therefore the liquid saturation, resulting from the capillary curves and the difference between the gas and the liquid pressures, becomes constant in M1 and M2. Let us write the equations satisfied by the asymptotic state at large times:

           P g (∞) = P g (0), φ (1) S (1) l (∞) + φ (2) S (2) l (∞) = φ (1) S (1) l (0) + φ (2) S (2) l (0) P (1) c (S (1) l (∞)) = P (2) c (S (2) l (∞)), P l (∞) = P g (∞) -P (1) c (S (1) l (∞)) = P g (∞) -P (2) c (S (2) l (∞)). (26)
These equations are satisfied for P l (∞) -2 10 7 Pa, S (1) l (∞) 0.844 and S (2) l (∞) 0.548. Let us finally comment Figure 14, showing the global density of gas. The integral of this quantity with respect to x should remain constant, and equal to the total amount of gas present in the domain. We observe that this property is graphically respected, in particular at the final time, where the area between the curve in M2 is equal to the difference between the initial and final curves in M1. Figure 14 confirms the lack of displacement of the gas phase in the short times (the curves remain nearly horizontal in domain M1 for t < 20000 s ). One can also observe that the penetration of the gas phase into the place previously occupied by the water phase in M2 is mainly driven by the porous medium type equation on the pressure.

Conclusion

The numerical examples show the good behaviour of the formulation introduced in this paper for two-phase flow modelling. Many advantages result from this formulation: it can handle gas injection case, and it can be used for the approximation of immiscible compressible flows. The numerical implementation of this pressure-pressure formulation, thanks to the extension of the SUSHI scheme to this two-phase flow problem, presents stability properties which make this scheme suited for simulations of industrial problems in which the capillary pressure phenomena are dominating (this is the case in the studies concerning nuclear waste storage, but not necessarily the case in oil recovery studies for example).

There now remains to show the properties of the scheme on 3D examples, and on various data with highly contrasted properties.

Nomenclature

Figure 1 :

 1 Figure 1: Capacity function

Figure 3 :Figure 4 :

 34 Figure 3: Gas pressure at different times (left: case 1, right: case 2).

Figure 5 :Figure 6 :

 56 Figure 5: Saturation at different times in Case 1.

Figure 7 :

 7 Figure 7: Dimensions of the domain (left), mesh used for the simulation (right) Porous medium M1 Porous medium M2 Characteristics of fluids Parameter Value Parameter Value Parameter Value k (m 2 ) 10 -20 k (m 2 ) 10 -19 D h l (m 2 .s -1 ) 0 φ 0.3 φ 0.05 H (mol.Pa -1 .m -3 ) 3.8 10 -10 A (Pa) 1.5 10 6 A (Pa) 10 10 6 µ l (Pa.s) 10 -3 B 0.06 B 0.412 µ g (Pa.s) 1.8 10 -5 a 0.25 a 1 ρ l (kg.m -3 ) 10 3 b 16.67 b 2.429 M w (kg.mol -1 ) 18 10 -3 c 1.88 c 1.176 M h (kg.mol -1 ) 28.96 10 -3 d 0.5 d 1

Figure 8 :

 8 Figure 8: Gas pressure (left) and liquid saturation (right) at different times: top 50 years, middle 1 000 years, bottom 10 000 years

Figure 9 :

 9 Figure 9: The geometry.

Figure 10 :

 10 Figure 10: Gas pressure (Pa) with respect to position (m) in Case 1 (no initial gas constituant in M2)

Figure 11 :

 11 Figure 11: Gas pressure (Pa) with respect to position (m) in Case 2 (water initially saturated with gas in M2).

Figure 12 :

 12 Figure 12: Liquid pressure (Pa) with respect to position (m).

Figure 14 :

 14 Figure 14: Global density of gas (kg/m 3 ) with respect to position (m).

pppp

  Mass flux for phase p = l, g g Gravity acceleration J c Mass diffusive flux for component c = h, w in phase p k Absolute permeability C(P) Capacity function µ p Viscosity of phase p = l, g n K,σ Unit vector normal to σ outward to K φ Porosity ρ p Density of phase p = l, g σ Face of a control volume x σ Barycentre of face σ x K Centre of control volume K C p Molar concentration of phase p = l, g C c Molar concentration of the component c = h, w in phase p = l, g D c Diffusion coefficient of component c = h, w in phase p = l, g d K,σ Orthogonal distance between point x K and face σ D h l,K,σ Approximation of diffusion flux of component h in phase l outward to K at face σ F K,σ Approximation of normal gradient outward to K integrated over face σ F c p,K,σ Approximation of Darcy flux of component c = h, w in phase p = l, g outward to K at face σ g Subscript for gaseous phase H Henry's law constant at the temperature of the domain h Superscript for gaseous component K Control volume, element of M K σ Cone with vertex x K and basis σ k rp Relative permeability of phase p = l, g L Control volume, element of M l Subscript for liquid phase m Parameter of Van Genuchten-Mualem law M c Molar mass of component c = h, w m c Volumic mass of component c = h, w in phase p = l, g n Parameter of Van Genuchten-Mualem law P c Capillary pressure P p Pressure of phase p = l, g P r Parameter of Van Genuchten-Mualem law R Ideal gas constant

Table 1 :

 1 Data of benchmark problem[START_REF] Smaï | Apparition-disparition de phase dans un écoulement diphasique eau/hydrogène en milieu poreux : Injection de gaz dans un milieu saturé en eau pure[END_REF]

			Characteristics of fluids
	Parameter	Value	Parameter	Value
	k(m 2 ) φ	5 10 -20 0.15	D h l (m 2 .s -1 ) µ l (Pa.s)	3 10 -9 10 -3
	P r (Pa)	2 10 6	µ g (Pa.s)	9 10 -6
	n	1.49	H ( mol.Pa -1 .m -3 ) 7.65 10 -6
	S lr	0.4	M w (kg.mol -1 )	18 10 -3
	S gr	0	M h (kg.mol -1 )	2 10 -3
	m	1 -1 n	ρ l (kg.m -3 )	10 3

Table 2 :

 2 Physical data

  Figure 13: Liquid saturation with respect to position (m).

	1										1
									0 s											1.e6 s
									100 s 1000 s		0.95									1.e7 s 1.e8 s
									5000 s											1.e9 s
	0.95								20000 s		0.9									1.e10 s
																				1.e11 s
											0.85
	0.9										0.8
											0.75
	0.85										0.7
											0.65
	0.8										0.6
											0.55
	0.75										0.5
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
	0.09										0.09
									0 s											1.e6 s
	0.08								100 s 1000 s		0.08									1.e7 s 1.e8 s
									5000 s											1.e9 s
	0.07								20000 s		0.07									1.e10 s 1.e11 s
	0.06										0.06
	0.05										0.05
	0.04										0.04
	0.03										0.03
	0.02										0.02
	0.01										0.01
	0										0
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1


O Touazi

E Chénier
email: eric.chenier@univ-mlv.fr

R Eymard

Simulation of natural convection with the Collocated Clustered Finite Volume scheme

This paper presents numerical results obtained in the case of natural convection within non constant fluid density, using the Collocated Clustered Finite Volume (CCFV) scheme. The continuous equations are first given in a dimensionless form. Then we present the finite volume scheme with the principles and the spatial discretization used. Analytical tests illustrate the numerical behavior of this scheme according to the type of grid, of the pressure stabilization method and check the robustness of this scheme. Next, the results obtained on the square thermally driven cavity under large temperature differences show that the CCFV scheme accurately fits the reference results.

Introduction

Various fluid flows are successfully modeled by the Navier-Stokes equations within the Oberbeck-Boussinesq approximation: the fluid properties are assumed to be constant in the mass, momentum and energy conservation equations, except in the buoyancy term where a linear dependence of the density with respect to the temperature is taken into account. This approximation is justified when the range of relative temperature variations is small enough (the mathematical validation of this approximation is examined in [START_REF] Gray | The validity of the Boussinesq approximation for liquids and gases[END_REF][START_REF] Mihaldjan | A rigourous exposition of the Boussinesq approximations applicable to a thin layer of fluid[END_REF]), and then leads to a simpler mathematical formulation than the complete equations. However, when the temperature gap increases significantly, the variations of the fluid properties cannot be neglected.

The Navier-Stokes and energy equations with variable fluid properties involve many non-linear terms, which have to be handled by numerical methods. Adaptations of the collocated finite volume schemes [START_REF] Archambeau | Code Saturne: A Finite Volume Code for Turbulent flows[END_REF][START_REF] Chénier | Numerical results using a colocated finite volume scheme on unstructured grids for incompressible fluid flows[END_REF][START_REF] Mathur | Pressure boundary conditions for incompressible flow using unstructured meshes[END_REF] to this case are often necessary in the engineering framework. In the Collocated Clustered Finite Volume method, the discrete unknowns are all located at the center of the cells of a unique mesh. Since the equivalent of the well-known "inf-sup" condition in the finite element setting (see [START_REF] Pironneau | Finite element methods for fluids[END_REF]), or also called LBB condition (see [START_REF] Brezzi | On the existence uniqueness and approximation of saddlepoint problems arising from Lagragian multipliers[END_REF]), is not satisfied by a collocated scheme, stabilization methods have to be implemented in order to avoid unrealistic pressure oscillations, sometimes called "checkerboard modes". Among the pressure stabilization techniques, let us cite the Brezzi and Pitkaranta approach [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF], first introduced in the framework of the finite element method, which consists in adding in the mass conservation equation a stabilization term proportional to the laplacian of the pressure. Another common pressure stabilization technique is the "Momentum interpolation Method", presented by Rhie and Chow in [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF] and usually practiced within the use of time stepping methods in finite difference or volume schemes. We have introduced in [START_REF] Chénier | Numerical results using a colocated finite volume scheme on unstructured grids for incompressible fluid flows[END_REF][START_REF] Eymard | On colocated clustered finite volume schemes for incompressible flow problems[END_REF] the Collocated Clustered Finite Volume (CCFV) scheme, presenting the following characteristics:

(1) This scheme holds for a large class of grids, the only geometric requirement being that in every control volume, one can select a particular point, called the center of the control volume, such that each interface between two adjacent control volumes is orthogonal to the line joining their centers. This condition applies for rectangular or parallelepipedic meshes, triangular (2D) or tetrahedral (3D) Delaunay meshes, and Voronoï meshes. (2) The scheme is fully coupled, and the discrete gradient operator is defined as the adjoint of the discrete divergence operator, thus ensuring the decrease of the kinetic energy for incompressible and isothermal fluid flows.

The discrete non-linear mass, momentum and energy equations are simultaneously solved using an under-relaxed Newton's method, and the linear systems are solved using an iterative linear solver, namely the BiCGStab method preconditioned by an incomplete LU factorization. (3) The pressure stabilization method is performed by selecting clusters of control volumes, i.e. small sets of control volumes, such that a stabilization term is imposed in the equation of the mass conservation, only on the edges between neighboring control volumes belonging to the same cluster. This procedure allows keeping stable discrete solutions, while the approximate solution is shown to converge to a solution of the Navier-Stokes equations.

We present in this paper the extension of the CCFV scheme to non constant fluid properties, thus permitting to handle the low Mach number equations for large temperature variations. Note that this extension leads to several new difficulties:

(1) New non-linear terms have to be handled.

(2) Since the density is not constant, one has to account for cross derivative terms involved by the viscous stress tensor. (3) The thermodynamic pressure must be computed in agreement with the temperature field.

We take advantage in this paper of the fact that the CCFV is natively a non-linear coupled scheme, the extension of which the new non-linear terms appears to be straightforward. New developments, based of recent research results on finite volume methods [START_REF] Eymard | A cell-centred finite volume approximation for anisotropic diffusion operators on unstructured meshes in any space dimension[END_REF] have allowed for accounting for cross derivative terms resulting from variable density.

The outlines of this paper are the following. We give in section 2 the fluid flow equations with non constant physical properties. In section 3 are presented the numerical scheme and the main features of its practical implementation. Numerical results are then described in section 4 showing some convergence properties, the effect of the pressure stabilization method and comparisons with a benchmark, and some conclusions are drawn in section 5.

Low Mach number approximation

The low Mach number model, suggested by Paolucci in the early 1980s [START_REF] Paolucci | Filtering of sound from the Navier-Stokes equations[END_REF], is used in studying very low velocity flows such as natural convection in a large temperature difference [START_REF] Paolucci | Transition to chaos in a differentially heated vertical cavity[END_REF][START_REF] Quéré | A Chebyshev Collocation Algorithm for 2D Non-Boussinesq Convection[END_REF]. In this approach, the total pressure is the sum of a mean pressure P * (t) in the fluid domain Ω, called the thermodynamic pressure, and fluctuations of the pressure p * ( x * , t * ), supposed to be such that Ω p * ( x * , t * )dv( x * ) = 0 for all time t * . From now on, the dimensional quantities are labbeled with * .

The dimensionless low Mach number equations for an ideal gas in an enclosure with impervious walls at the initial constant state (P * 0 , T * 0 ) and under the gravitational force per mass unit -g * e z are:

∂ ∂t ρ + ∇.(ρ v) = 0 (1) ∂ ∂t (ρ v) + ∇.(ρ v ⊗ v) = -∇p - RaP r ε T ρ e z + P r ∇. τ (2) ∂ ∂t (ρT ) + ∇.(ρT v) = γ -1 ε T γ dP dt + ∇.(k ∇T ) (3) Ω ρ dv( x) = Ω dv( x), Ω p( x, t)dv( x) = 0 (4)
τ represents the viscous stress tensor, that we assume to be given by the Stokes relation with I, the identity tensor:

τ = µ ∇ v + ( ∇ v) t - 2 3 ( ∇. v)I (5)
Notice that radiative fluxes are neglected in the energy equation and the diffusive heat flux is expressed using the Fourier law q = -k ∇T . As ρ is supposed to verify the ideal gas law, we obtain:

ρ = P 1 + ε T T (6)
where ρ, v, p, P , and T are respectively the non-dimensional density, velocity, spatially pressure fluctuations (often called dynamic pressure), thermodynamic pressure and temperature where ρ

= ρ * /ρ * 0 , v = v * /v * 0 , p = p * /(ρ * 0 v * 0 2
), P = P * /P * 0 and T = (T * -T * 0 )/∆T * . The space variable x * and velocity have been normalized by the reference height H * , and by the ratio between the thermal diffusivity at the reference state, a * 0 , and H * . The dynamic viscosity and thermal conductivity are normalized by their reference values at (P * 0 , T * 0 ). Three dimensionless parameters appear in these equations: the Rayleigh number Ra = (g * H * 3 ε T )/(a * 0 ν * 0 ), the Prandtl number P r = ν * 0 /a * 0 and the parameter ε T = ∆T * /T * 0 .

The next section of the paper is devoted to the description of the discrete approximations of terms involving in the low Mach number model.

The Collocated Clustered Finite Volume scheme

We now intend to present the scheme for the approximation of a solution (P, p, T, v) to system (1-3) with ( 4)- [START_REF] Pironneau | Finite element methods for fluids[END_REF]. Only the space discretization is presented with centered approximations. Upwind schemes and time discretization can be found in [START_REF] Chénier | Numerical results using a colocated finite volume scheme on unstructured grids for incompressible fluid flows[END_REF] and convergence proofs in [START_REF] Eymard | On colocated clustered finite volume schemes for incompressible flow problems[END_REF].

In the following subsections, we first give definitions related to meshes, then we present the space discretization of system (1-3) and associated principles for this finite volume scheme.

Admissible meshes

The domain Ω is completely recovered by a family M of unstructured convex polyhedra, called control volumes. The detailed notations for the geometric quantities involved by the discretization are illustrated in figure 1 in the 2dcase.

Let us introduce some notations: we denote by M the set of all the control volumes. For each control volume K, we denote by x K the coordinates of its center, by m(K) its area (if d = 2) or its volume (if d = 3), by N K the set composed by the neighbors of K. We denote by E K,int the set of the internal edges and by E K,ext , the set of the external edges. If the control volume L is a neighbor of the control volume K, we denote by m(σ KL ) the length (d = 2) or the area (d = 3) of the interface σ KL between K and L (σ KL ∈ E K,int ).If the control volume K has an edge σ at the boundary domain (σ ∈ E K,ext ), we denote by m(σ) the length (d = 2) or the area (d = 3) of the interface σ at the boundary. We denote by x σ the orthogonal projection of x K on σ, by n K,σ the unit vector normal to σ, outward of K, and by d K,σ , the distance between x K and x σ . We then denote by x Kσ the center of gravity of the interface σ.

The main objective to elaborate this numerical scheme was its simplicity, in particular for the expression of diffusion fluxes. So the approximation of the flux normal to the edge of a control volume σ KL depends only on two adjacent cells what permits to systematically satisfy monotonic solutions for diffusive equations. To allow it, the requirements to mesh the domain are:

σ KL ⊥ ---→ x K x L (7) ---→ x K x L • n KL > 0 (8)
From ( 7), we get the orthogonality property whereas relation [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF] defines the order of the cell center occurence in the direction of n KL . These conditions are checked for rectangular or parallelepiped and dual Delaunay-Voronoï grids.

Recall that a Voronoï mesh is such as the control volume K is as the set of points closer to x K than to any other x L . Therefore, the edge separating cells K and L belongs to the mid-plan between x K and x L what implies that requirements [START_REF] Brezzi | On the existence uniqueness and approximation of saddlepoint problems arising from Lagragian multipliers[END_REF] and ( 8) are automatically checked. For triangular (2d) or tetrahedral (3d) grids, the center x K of the cell K is the circumcenter (condition ( 7)); the constraint ( 8) is automatically verified once the Delaunay criterion is satisfied.

Spatial discretization of the Navier-Stokes and energy equations

The objective in the development of this numerical scheme is that solutions of the discrete Navier-Stokes and energy equations p K , v K , T K , accurately approximate p, v, T at the point x K of each control volume K ∈ M. Our discrete spatial operators are also constructed so that, by multiplying the momentum and the thermal equations expressed on cell K with the velocity v K and the temperature T K , and by summing over all the cells, we recover discrete approximations of the kinetic and thermal energy balances. For example, the works of the pressure forces and of the non-linear inertial forces must vanish under Boussinesq approximation.

For the sake of clarity, the discrete spatial operators are presented for enclosures with impervious walls and only for Dirichlet boundary conditions. The following notations are used in the next sections. We denote by (A( x)) D the discrete approximation of the continuous form A( x) and by (A( x)) D,K the discrete approximation evaluated at the cell center x K ; v (i) is the i-th component of the velocity vector v and ∂ i marks out the partial derivative into the i-th direction.

Mass equation

Let us first give the approximation of the steady mass equation integrated over the cell K, K ∇.(ρ v) dv( x) that is:

m(K) ∇.(ρ v) D,K = L∈N K F KL (9)
where F KL is the approximation of the mass flux

F KL = σ KL ρ v( x)• n KL ds( x)
between two neighboring control volumes K and L. Notice that no boundary contribution is introduced because of the impervious walls.

To preserve a good coupling between pressure and velocity, the approximation of the velocity crossing the edges σ KL must be second order accurate at the gravity center of the face [START_REF] Touazi | Etude numérique d'un schéma de volumes finis sur des maillages non structurés pour des écoulements de fluide visqueux à densité constante ou faiblement variable[END_REF]. This condition is easily fulfilled by using a linear interpolation of the velocity v( x KL ) = I∈M α I v( x I ). Practically, nodes x K and x L are involved plus, eventually, one or two other nodes close to x KL .

Without any loss of generality, we now suppose that only centers x K and x L play a part in the interpolation. This case corresponds to a mesh where the intersection point between the line ( x K , x L ) and the edge σ KL is exactly x KL , the gravity center of the face σ KL . Parallelepipeds in 3d [START_REF] Chénier | Numerical results using a colocated finite volume scheme on unstructured grids for incompressible fluid flows[END_REF] or rectangles and Delaunay triangles in 2d satisfy these requirements. The discrete mass flux then writes:

F KL = ρ KL m(σ KL )(a KL v K + a LK v L ) • n KL -λ KL m(σ KL ) d KL (p L -p K ) (10)
with

a KL = ----→ x KL x L . n KL d KL , a KL + a LK = 1
The density ρ KL on interface σ KL is also evaluated with a linear interpolation.

A stabilization method is necessary to prevent from the apparition of pressure oscillations in collocated schemes. We use the method presented in [START_REF] Chénier | Numerical results using a colocated finite volume scheme on unstructured grids for incompressible fluid flows[END_REF] that consists in modifying the mass equation by introducing a stabilization contribution which corresponds to the second term of the right-hand-side of equation [START_REF] Eymard | On colocated clustered finite volume schemes for incompressible flow problems[END_REF], where λ KL is a stabilization parameter defined on each interface between two control volumes K and L.

Recall briefly the principle of pressure stabilization method. The family of control volumes M is partitioned into disjoint subsets including groups of neighboring control volumes, with any algorithm ensuring that two control volumes belonging to the same subset are not too distant one from each other (Fig. 2). An example of such an algorithm can be found in [START_REF] Chénier | Numerical results using a colocated finite volume scheme on unstructured grids for incompressible fluid flows[END_REF].

Then λ KL is chosen according to the following rule:

• λ KL = λ for any pair of neighboring control volumes K and L belonging to the same subset, • λ KL = 0 otherwise. The value λ is a positive value, chosen large enough for preventing from the apparition of instabilities. The main advantage of gathering the cells by clusters is that the mass redistribution due to the stabilization term is limited to each cluster.

Transport terms

For φ ∈ {v (i) , T }, the discrete expression of the non-linear term integrated over the cell K is:

m(K) ∇.(ρ vφ) D,K = L∈N K G KL ( 11
)
where G KL is the approximation of the flux

G KL = σ KL (ρ v( x)• n KL )φ( x)ds( x).
To satisfy the discrete kinetic and thermal energy balances, the transport velocity must be equal to the discrete mass flux crossing each edge (Eq. ( 10)) and the centered approximation of φ must be used:

G KL = F KL φ K + φ L 2
By taking into account the mass equation for a steady flow, L∈N K F KL = 0, relation ( 11) also writes:

m(K) ∇.(ρ vφ) D,K = L∈N K F KL φ L -φ K 2

Pressure

The discrete approximation of the pressure gradient is based on its duality with the discrete velocity divergence:

K∈M m(K) v K • ( ∇p) D,K = - K∈M     L∈N K m(σ KL )(a KL v K + a LK v L ) n KL   p K  
Thus we obtain

m(K)( ∇p) D,K = L∈N K m(σ KL )a KL (p L -p K ) n KL (12)
By using the geometric property L∈N K m(σ KL ) n KL + σ∈E K,ext m(σ) n K,σ = 0, another form for the pressure gradient can be established:

m(K)( ∇p) D,K = L∈N K m(σ KL )(a KL p L + a LK p K ) n KL + σ∈E K,ext m(σ)p K n K,σ (13)
Notice that relations [START_REF] Paolucci | Filtering of sound from the Navier-Stokes equations[END_REF] and ( 13) do not provide consistent approximations of the equivalent forms K ∇p dv( x) and δK p nds( x), except when a KL = 0.5.

Viscous diffusion

In the low Mach number model, the viscous term integrated over the cell K writes:

K ∇. τ dv( x) = K ∇.(µ ∇ v) dv( x) + K ∇.(µ ∇ t v) dv( x) - 2 3 K ∇(µ ∇. v) dv( x) = K a(µ, v) dv( x) + K b(µ, v) dv( x) - 2 3 K c(µ, v) dv( x)
Let us now give the approximations of these three integrals.

• The approximation of a(µ, v) is:

m(K)( a(µ, v)) D,K = L∈N K A KL + σ∈E K,ext A K,σ (14)
where A KL and A K,σ are the approximations of

A KL = σ KL ∇ v( x)• n KL ds( x) and A K,σ = σ ∇ v( x) • n K,σ ds( x).
Thanks to the mesh properties ( 7) and ( 8), we simply obtain:

A KL = µ KL m(σ KL ) d KL ( v L -v K ) and A K,σ = µ σ m(σ) d K,σ ( v σ -v K ) (15)
where v σ is the approximation of v( x σ ), µ KL and µ σ are the dynamic viscosity computed on edges. • The discrete contribution K b(µ, v) dv( x) is achieved by considering the weak form:

d i=1 Ω u (j) ∂ i (µ∂ j v (i) )dv( x) = - d i=1 Ω ∂ i u (j) µ∂ j v (i) dv( x) (16)
with u (j) the j-th component of a vectoriel test function u, u (j) being null on Dirichlet boundary conditions. By using the i-th component of the discrete strong gradient expression (∂ i φ) D,K [START_REF] Eymard | A cell-centred finite volume approximation for anisotropic diffusion operators on unstructured meshes in any space dimension[END_REF],

(∂ i φ) D,K = 1 m(K) L∈N K m(σ KL ) d KL (φ L -φ K )(x (i) KL -x (i) K )+ σ∈Eext m(σ) d K,σ (φ σ -φ K )(x (i) Kσ -x (i) K ) ( 17
)
the discretization of the right-hand-side of (Eq. 16) writes:

d i=1 Ω µ ∂ i u (j) ∂ j v (i) dv( x) D = d i=1 K∈M m(K)µ K (∂ i u (j) ) D,K (∂ j v (i) ) D,K
The k-th component of the discrete contribution of K b(µ, v)dv( x) is obtained by substituting (∂ i u (j) ) D,K and (∂ j v (i) ) D,K by the expressions [START_REF] Eymard | On a stabilized colocated Finite Volume scheme for the Stokes problem[END_REF] for φ = u (j) and φ = v (i) and by choosing the test function u

(j) L = δ KL δ jk : m(K)(b (k) (µ, v)) D,K = L∈N K B (k) KL + σ∈E K,ext B (k) K,σ (18)
with

B (k) KL = m(σ KL ) d KL d i=1 µ K (∂ k v (i) ) D,K (x (i) KL -x (i) K ) + µ L (∂ k v (i) ) D,L (x (i) L -x (i) KL ) B (k) K,σ = m(σ) d K,σ d i=1 µ K (∂ k v (i) ) D,K (x (i) Kσ -x (i) K ) + µ L (∂ k v (i) ) D,L (x (i) L -x (i)
Lσ )

• The k-th component of discrete contribution of K c(µ, v)dv( x) writes:

m(K)(c (k) (µ, v)) D,K = L∈N K C (k) KL + σ∈E K,ext C (k) K,σ (19)
with

C (k) KL = m(σ KL ) d KL d i=1 µ K (∂ i v (i) ) D,K (x (k) KL -x (k) K ) + m(σ KL ) d KL d i=1 µ L (∂ i v (i) ) D,L (x (k) L -x (k) KL ) C (k) K,σ = m(σ) d K,σ d i=1 µ K (∂ i v (i) ) D,K (x (k) Kσ -x (k) K ) + m(σ) d K,σ d i=1 µ L (∂ i v (i) ) D,L (x (k) L -x (k) Lσ )
Note that the use of the strong gradient [START_REF] Eymard | On a stabilized colocated Finite Volume scheme for the Stokes problem[END_REF] in K a(µ, v)dv( x) may give rise to numerical instabilities. Therefore, to prevent from these oscillations, relation [START_REF] Quéré | A Chebyshev Collocation Algorithm for 2D Non-Boussinesq Convection[END_REF] must absolutly be applied instead.

Thermal diffusion

In the energy equation, the discrete contribution (h(k, T )) D,K provided by the approximation of the diffusion term K h(k, T ) dv( x) = K ∇.(k ∇T )) dv( x) is similar to relations ( 14) with (15):

m(K)(h(k, T )) D,K = L∈N K H KL + σ∈E K,ext H K,σ
The approximations H KL and H K,σ of the conducted heat fluxes on σ KL ,

H KL = σ KL k ∇T ( x)
• n KL ds( x), and on σ, H K,σ = σ k ∇T ( x) • n K,σ ds( x), write:

H KL = k KL m(σ KL ) d KL (T L -T K ) and H K,σ = k σ m(σ) d K,σ (T σ -T K )
where T σ is the approximation of T ( x σ ), k KL and k σ are thermal conductivities on σ KL and σ.

Discrete equations and solvers

We get a system of non-linear equations, which depends on the unknowns of the control volume K itself, on the unknowns of the neighbors L ∈ N K of K or on boundary σ ∈ E K,ext , and also, because of ( b(µ, v)) D and ( c(µ, v)) D exist, on unknowns of the neighbors M ∈ N L of L. Note that under the Boussinesq approximation, the mass equation ∇. v = 0 implies c(µ, v) = 0 and if furthermore the viscosity is constant, b(µ, v) = 0. In that case, the equation on cell K only depends on its neighbors L ∈ N K or on boundary σ ∈ E K,ext .

In our implementation, the non-linear system is solved using an under-relaxed Newton method. The resulting linear equations are then solved by an iterative linear solver using the BiCGStab method preconditioned by an incomplete LU factorization, with variable fill-in levels (a detailed description can be found in [START_REF] Chénier | Numerical results using a colocated finite volume scheme on unstructured grids for incompressible fluid flows[END_REF]). The thermodynamic pressure P is then updated at the end of each Newton iteration, solving the discrete form of the total mass conservation equation given in (4).

Natural convection results

The computational domain is meshed either by rectangular cells or by Delaunay triangles thanks to a free software developed by Saltel and Hecht [START_REF] Saltel | EM C 2 Un logiciel d'édition de maillages et de contours bidimensionnels[END_REF].

Analytical tests

We consider the stationary fluid flow and heat transfer in a square cavity with a low Mach number model. The dimensionless velocity, pressures, density and temperature fields are solutions of equations [START_REF] Gray | The validity of the Boussinesq approximation for liquids and gases[END_REF][START_REF] Mihaldjan | A rigourous exposition of the Boussinesq approximations applicable to a thin layer of fluid[END_REF][START_REF] Archambeau | Code Saturne: A Finite Volume Code for Turbulent flows[END_REF][START_REF] Chénier | Numerical results using a colocated finite volume scheme on unstructured grids for incompressible fluid flows[END_REF][START_REF] Mathur | Pressure boundary conditions for incompressible flow using unstructured meshes[END_REF][START_REF] Pironneau | Finite element methods for fluids[END_REF] where the time derivatives are zero. No-slip boundary conditions are applied on the isothermal walls.

To study the accuracy of the numerical method in the non-Boussinesq regime, forcing terms are added in the momentum and energy equations in order to converge towards the prescribed solution: (1) ) sin 2 (πx (2) ) (1) ) sin 2 (πx (2) ) v (1) = ∂Ψ ∂x (2) = 2π cos(πx (2) ) sin(πx (2) ) sin 2 (πx (1) )

P = 1 T = sin 2 (πx
ρ = P 1 + ε T T p = sin 2 (πx
v (2) = -∂Ψ ∂x (1) = -2π cos(πx (1) ) sin(πx (1) ) sin 2 (πx (2) )

The velocity and the density (i.e. the temperature) are chosen so that ∇.(ρ v) = 0. Sutherland laws are used to evaluate the dimensionless dynamical viscosity as a function of the temperature:

µ = (ε T T + 1) 3/2 1 + S * /T * 0 1 + ε T T + S * /T * 0 ( 20
)
where S * = 110.5 K is the Sutherland constant. The variation laws of the thermal conductivity k and dynamical viscosity µ, as a function of the temperature, are supposed identical, and computations are carried out with Ra = 1 and P r = 1. The deviation from the analytical solution is measured by the discrete L 2 -norm:

||e(φ)|| 2 2 = K∈M m(K)(φ K -φ( x K )) 2 , φ ∈ {v (i) , T, p, ρ}

Spatial accuracy

Figures 3a-3d show the L 2 -norm of errors between the analytical and numerical solutions as a function of mean mesh size h = 1/ √ N , N being the number of triangles, and with three reduced temperature differences, ε T ∈ {0, 1, 2}. Remark that the Boussinesq approximation is recovered when ε T = 0. For all ε T , the velocity, temperature and density are second-order accurate whereas the pressure is, like most of collocated schemes, first-order accurate only.

Comparison of pressure stabilization methods

The simulations are performed for the Boussinesq model, namely with ε T = 0, and for Delaunay triangular meshes. Our pressure stabilization method, described in section 3.2.1, is compared with two other methods.

• The Brezzi and Pitkaranta method (BP) [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF] consists in introducing a stabilization contribution λh α (∆p) D in the mass equation, where α ∈]0, 2[, h is the characteristic mesh size, α ∈ is a positive small value and ∆ D p ≡ a D (1, p) is the discrete laplacian operator applied on pressure and defined by ( 14):

( ∇. v) D = λh α (∆p) D

• The Rhie and Chow interpolation method (RC) [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF] can also be viewed as a penalization technique of the mass equation:

( ∇. v) D = λ (∆p) D -∇.( ∇p) D D
where the second term of the right-hand-side is constructed with the strong gradient [START_REF] Eymard | On a stabilized colocated Finite Volume scheme for the Stokes problem[END_REF].

In figures 4a-4b are drawn the velocity and pressure errors as a function of λh α , for three pressure stabilization techniques and for a mesh size h = 0.019. The accuracy of the all solutions decreases if the coefficient λh α is too small: the pressure stabilization is insufficient and oscillations appear. For large coefficients, the errors for the velocity and pressure increase except when the stabilization is based on clusters. Our clusterization technique is very few sensitive to λh α but the (RC)-method seems to provide more accurate pressure solutions for a wide gap of the coefficient.

Figures 5a-5d illustrate the convergence properties of these methods as a function of h and for two coefficients, λ ∈ {1, 0.0001}. For λ = 1, the slopes of the curve for the (BP)-method become worse, less than 2 and 1 for the velocity and pressure respectively. No significant difference distinguishes the convergence behavior of the cluster and (RC)-methods. The advantage of the cluster stabilization method with respect to the Rhie and Chow technique stems from the stability property of the resulting scheme. While the clusters and (BP)-methods introduce stable contributions in the discrete scheme and thus provide a L 2 stability for unsteady isothermal solutions [START_REF] Chénier | Numerical results using a colocated finite volume scheme on unstructured grids for incompressible fluid flows[END_REF][START_REF] Eymard | On colocated clustered finite volume schemes for incompressible flow problems[END_REF][START_REF] Eymard | On a stabilized colocated Finite Volume scheme for the Stokes problem[END_REF][START_REF] Eymard | Convergence Analysis of a Colocated Finite Volume Scheme for the Incompressible Navier-Stokes Equations on General 2D or 3D Meshes[END_REF], no such property can easily be proved for the (RC)-method because of the negative contribution -∇.( ∇p) D D in the mass equation.

Benchmark comparisons

We consider the problem handled in the benchmark [START_REF] Paillère | Modelling of natural convection flows with large temperature differences: A benchmark problem for low Mach number solvers ESAIM[END_REF], stated in the above dimensionless form, which concerns the natural convection flow in a differentially heated square cavity. The vertical walls are isothermal at the uniform temperatures T (0, y) = T h = 1/2 and T (1, y) = T c = -1/2 with 0 ≤ y ≤ 1, and the bottom and top boundaries are adiabatic. The ratio ε T is given and equal to 1.2, and therefore the Boussinesq approximation no longer holds. The other dimensionless numbers are Pr = 0.71, γ = 1.4, Ra = 10 6 or Ra = 10 7 .

As in [START_REF] Paillère | Modelling of natural convection flows with large temperature differences: A benchmark problem for low Mach number solvers ESAIM[END_REF], we use the Sutherland law (Eq. 20) for the viscosity and assume µ(T )/k(T ) constant, which leads to a similar law for the thermal conductivity. We recall the definition of the local, Nu( x), and average, Nu x , Nusselt numbers:

Nu( x) = k(T ( x)) ∇T ( x) • n ∂Ω ( x) and Nu x = 1 0 Nu( x)dy,
The computations were performed either on a non-uniform distribution of Delaunay triangles or on different orthogonal grids where the sides of the rectangular control volumes follow a sinusoidal law in each direction to accurately represent the boundary layers. The stabilization parameter is equal to λ = 10 -6 and the iLU preconditioning is based on a three levels of fill-in.

Rayleigh number 10 6

In Table 1 are indicated the relative errors, with respect to solutions obtained with our scheme with a 304 × 304 rectangular grid, for the results computed on Delaunay triangles and drawn from [START_REF] Paillère | Modelling of natural convection flows with large temperature differences: A benchmark problem for low Mach number solvers ESAIM[END_REF]. Relative differences of benchmark solutions being smaller than 1%, our reference result is in good agreement with the literature. Bold underlines values smaller than 0.005% and show the good agreement with Vierendeels's results. The accuracy of results obtained with Delaunay triangular cells are good, except for the Nusselt number along the cold wall, max y (Nu(1; y)) and min y (Nu(1; y)). For non-Boussinesq flow model, the boundary layer along the cold wall is thicker than it is along the hot wall. The mesh size used being identical along the two vertical walls, it was therefore obvious that the accuracy of the heat transfer will be worst on the cold surface. The minima of the Nusselt number value are also larger than their corresponding maxima. The explanation is due to the fact that the coordinates, where min y (Nu(0; y)) and min y (Nu(1; y)) occur, are very close to the horizontal walls where a finer mesh should be desirable too.

The effect of the mesh size was also studied but it is not presented in this paper. Quadratic convergence has successfully been checked for the set of results except for minimal values of the Nusselt number when triangular meshes were considered.

Rayleigh number Ra = 10 7

For Ra = 10 7 , only results obtained with rectangular grid are shown, the computational cost to solve the Navier-Stokes and energy equations on a fine enough triangular grid was judged too expensive. As shown in Table 2, our results are in excellent agreement with benchmark solutions [START_REF] Paillère | Modelling of natural convection flows with large temperature differences: A benchmark problem for low Mach number solvers ESAIM[END_REF]. Relative differences less than 0.005% are anew exhibited in bold. As for Ra = 10 6 , the second order accuracy of the scheme was successfully checked for the set of results.

Conclusion

Extention of the Collocated Clustered Finite Volume scheme to the low Mach number model was presented in detail. Convergence property of the new scheme was checked by considering an analytical solution. Effect of the pressure stabilization method was studied with respect to other stabilization techniques. Applications to natural convection problems with large temperature differences have led to results which are in good agreement with the literature.

List of Figures

x K n KL σ KL K L x KL x L d KL d K,σ n K,σ
x Kσ x σ σ Fig. 1. Two neighboring control volumes K and L. One assumes that for all control volume K there exists a point denoted by x K in K such that, for any neighbor L of K, the straight line ( x K , x L ) is orthogonal to the common edge σ KL . Then the length between x K and x L is denoted by d KL , and the unit vector, normal to σ KL and outward to K, is denoted by n KL . We then denote the center of gravity of the interface σ KL by x KL . Similar notations are used when an edge σ belongs to δΩ, the boundary of the domain Ω. Comparison with results given in [START_REF] Paillère | Modelling of natural convection flows with large temperature differences: A benchmark problem for low Mach number solvers ESAIM[END_REF] for Ra = 10 7

ε T =0 ε T =1 ε T =2 (c) ||e(T )|| 2 -4.5 -4 -3.5 -3 -2.5 -2 -2.4 -2.2 -2 -1.8 -1.6 -1.4 -1.2 log 10 (||e(ρ)|| 2 ) log 10 (h) ε T =1 ε T =2 (d) ||e(ρ)|| 2

22

Fig. 2 .-

 2 Fig. 2. Delaunay grid and associated clusters

Fig. 3 .Fig. 4 .

 34 Fig. 3. Spatial accuracy for velocity, pressure, temperature and density

1

  Two neighboring control volumes K and L. One assumes that for all control volume K there exists a point denoted by x K in K such that, for any neighbor L of K, the straight line ( x K , x L ) is orthogonal to the common edge σ KL . Then the length between x K and x L is denoted by d KL , and the unit vector, normal to σ KL and outward to K, is denoted by n KL . We then denote the center of gravity of the interface σ KL by x KL . Similar notations are used when an edge σ belongs to δΩ, the boundary of the domain Ω.

			19
	2	Delaunay grid and associated clusters	20
	3	Spatial accuracy for velocity, pressure, temperature and
		density	21
	4	Effect of the stabilization parameter λh

α on the L 2 -norm of the velocity and pressure, α = 0 for the cluster method and the Rhie and Chow technique, (RC), whereas α = 1 for the Brezzi and Pitkaranta method, (BP). 22 5 Velocity and pressure errors as a function of h and λ for the cluster, Rhie and Chow (RC) and Brezzi and Pitkaranta (BP) methods. 23

Table 1

 1 Velocity and pressure errors as a function of h and λ for the cluster, Rhie and Chow (RC) and Brezzi and Pitkaranta (BP) methods. Comparison with results given in[START_REF] Paillère | Modelling of natural convection flows with large temperature differences: A benchmark problem for low Mach number solvers ESAIM[END_REF] for Ra = 10 6

		-1.2								-1.2
		-1.6								-1.6
	log 10 (||e(u)||) 2	-2.8 -2.4 -2							log 10 (||e(u)||) 2	-2.8 -2.4 -2
		-3.2					Clusters (BP), α = 1 (RC)		-3.2					Clusters (BP), α = 1 (RC)
		-3.6								-3.6
		-2.1	-2	-1.9	-1.8	-1.7	-1.6	-1.5	-1.4	-2.1	-2	-1.9	-1.8	-1.7	-1.6	-1.5	-1.4
					log 10 (h)							log 10 (h)
			(a) Velocity, λ = 1				(b) Velocity, λ = 0.0001
		1.2
		0.8								1.2
		0.4								0.8
		0								0.4
	log 2 10 (||e(p)||)	-1.6 -1.2 -0.8 -0.4							2 log 10 (||e(p)||)	0 -1.6 -1.2 -0.8 -0.4					Clusters (BP), α = 1 (RC)
		-2.4 -2				Clusters (RC) (BP), α = 1			-2.4 -2
		-2.8								-2.8
		-2.1	-2	-1.9	-1.8	-1.7	-1.6	-1.5	-1.4	-2.1	-2	-1.9	-1.8	-1.7	-1.6	-1.5	-1.4
					log 10 (h)							log 10 (h)
			(c) Pressure, λ = 1				(d) Pressure, λ = 0.0001
	Fig. 5.

Table 2

 2