Rémy Garcia

Claude Michel

Marie Pelleau

Michel Rueher

student Garcia Rémy

Towards a

Keywords: floating-point numbers, round-off error, constraints over floating-point numbers, domain of errors

come

Introduction

Floating-point computations induce errors due to rounding operations required to close the set of floating-point numbers. These errors are symptomatic of the distance between the computation over the floats and the computation over the reals. Moreover, they are behind many problems, such as the precision or the numerical stability of floating-point computation. Especially when users omit to take into account the nature of floating-point arithmetic and use it directly like real number arithmetic. Identifying, quantifying and localizing those errors are tedious tasks that are difficult to achieve without tools automating it. A well-known example of computations deviation due to errors on floating-point numbers is Rump's polynomial [START_REF] Rump | Algorithms for verified inclusions: Theory and practice[END_REF]:

333.75b 6 + a 2 (11a 2 b 2 -b 6 -121b 4 -2) + 5.5b 8 + a 2b
where a = 77617 and b = 33096. The exact value of this expression, computed using the GMP library, is -54767 66192 ≈ -0.827396056. However, when this expression is evaluated on simple floats with a rounding mode set to the nearest even, the computed result is ≈ -6.3382530011411×10 29 , which is far apart from the real value. The difference between these two results, about -6.3382530011411×10 29 , emphasizes the need for round-off error analysis tools.

Floating-point computation errors have been the subject of many works based on an overestimation of actual errors. Let us mention the abstract interpreter Fluctuat [START_REF] Goubault | Static analysis of numerical algorithms[END_REF][START_REF] Ghorbal | A logical product approach to zonotope intersection[END_REF] that combines affine arithmetic and zonotopes to analyze the robustness of programs over floating-point numbers. Another more recent work, PRECiSA [START_REF] Titolo | An abstract interpretation framework for the round-off error analysis of floating-point programs[END_REF][START_REF] Moscato | Automatic estimation of verified floating-point round-off errors via static analysis[END_REF], relies on static analysis to evaluate round-off errors in a program. Nasrine Damouche [START_REF] Damouche | Improving the numerical accuracy of programs by automatic transformation[END_REF] and Eva Darulova [START_REF] Darulova | Sound compilation of reals[END_REF] have developed techniques for the automatic enhancement of numerical code. Their approach is based on an evaluation of round-off errors to estimate the distance between the expression over floats and the expression over reals. These approaches compute an error estimation which can be refined by splitting the search space into subdomains. However, it is not possible to directly reason on those errors, for example, by computing input values that exercise a given error. In order to overcome this lack of reasoning capabilities and to enhance the analysis of errors, we propose to incorporate in a constraint solver over floats [START_REF] Zitoun | Search strategies for floating point constraint systems[END_REF][START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF][START_REF] Botella | Symbolic execution of floating-point computations[END_REF][START_REF] Michel | Exact projection functions for floating point number constraints[END_REF][START_REF] Michel | Solving constraints over floating-point numbers[END_REF], the domain of errors which is dual to the domain of values. Both domains are associated with each variable of the problem. Additionally, we introduce projection functions that filter those domains as well as mechanisms required for the analysis of errors. More precisely, we focus on the analysis of deviation between computations over the floats and computations over the reals.

Our approach is based on interval arithmetic for approximating the domains of errors. In addition, a search applied on reduced domains computes input values that satisfy constraints on errors. We deliberately ignore the possibility of an initial observational error on input data even though initial computational errors are handled. Thus input data are assumed with an initial error of zero. For the sake of simplicity, we restrain ourselves to the four classical arithmetic operations. This simplification permits exact computation of values over reals1 for both values of expressions and computation of errors. Finally, the rounding mode is left to default, i.e. a rounding to the nearest even.

Notation and definitions

Floating-point numbers

The set of floating-point numbers is a finite subset of the rationals that has been introduced to approximate real numbers on a computer. The IEEE standard for floating-point arithmetic [9] defines the format of the different types of floatingpoint numbers as well as the behavior of arithmetic operations on these floatingpoint numbers. In the sequel, floating-point numbers are restricted to the more common one i.e. simple binary floating-point numbers represented using 32 bits and double binary floating-point numbers represented using 64 bits.

A binary floating-point number v is represented by a triple (s, e, m) where s is the sign of v, e, its exponent and m, its mantissa. When e > 0, v is normalized and its value is given by :

(-1) s × 1.m × 2 e-bias
where the bias allows us to represent negative values of the exponent. For instance, for 32 bits floating-point numbers, the size of s is 1 bit, the size of e is 8 bits, the size of m is 23 bits and the bias is equal to 127.

x + denotes the smallest floating-point number strictly larger than x while x -denotes the largest floating-point number strictly smaller than x. In other words, x + is the successor of x while x -is its predecessor.

An U lp, which stands for unit in the last place, is the distance which separate two consecutive floating-point numbers. However, this definition is ambiguous for floats that are a power of 2 like 1.0: in such a case, and if x > 0, then x + -x = 2 * (x -x -). To make this point clear, an explicit formulation of this distance is used whenever required.

Quantification of computation deviations

Computation on floating-point numbers is different from computation over real numbers due to rounding operations. Since the set of floating-point numbers is a finite subset of the real, in general, the result of an operation on the floats is not a float. In order to close the set of floating-point numbers for those operations, the result should be rounded to the nearest float according to a direction chosen beforehand.

The IEEE 754 norm [9] defines the behavior of floating-point arithmetic. For the four basic operations, it requires correct rounding, i.e. the result of an operation over the floats must be equal to the rounding of the result of the equivalent operation over the reals. More formally, z = x y = round(x • y) where z, x and y are floating-point numbers, is one of the four basic arithmetic operations on the floats, namely, ⊕, , ⊗, and , • being the equivalent operation on the reals, and round being the rounding function. This property bounds the error introduced by an operation over floats to ± 1 2 ulp(z) for correctly rounded operations with a rounding mode set to round to the nearest even float, which is the most frequent rounding mode.

When the result of an operation on the floats is rounded, it is different from the one expected on the reals. Moreover, each operation that belongs to a complex expression is likely to introduce a difference between the expected result on the reals and the one computed on the floats. Whereas for a given operation, the computed float is optimal in terms of rounding, the accumulation of these approximations can lead to significant deviations, like the one observed with Rump's polynomial.

Compared to its equivalent over the reals, the deviation of a computation over the floats takes root in each elementary operation. Therefore, it is possible to rebuild it from the composition of each elementary operation behavior. Input variables can come with errors attached due to previous computations. For example, for the variable x, the deviation on the computation of x, e x , is given by e x = x R -x F where x R and x F denote the expected results on the reals and on the floats respectively. Contrary to an observational error, e x is signed. This choice is required to capture correctly specific behaviors of floating-point computations, such as error compensations.

As such, computation deviation due to a subtraction can be formulated as follows: for z = x y, the error on z, e z , is equal to (x R -y R ) -(x F y F ). As e x = x R -x F and e y = y R -y F , we have

e z = ((x F + e x ) -(y F + e y )) -(x F y F )
The deviation between the result on the reals and the result on the floats for a subtraction can then be computed by the following formula:

e z = e x -e y + ((x F -y F ) -(x F y F ))
In this formula, the last term, ((x F -y F ) -(x F y F )), characterizes the error produced by the subtraction operation itself. Let e denotes this subtraction operation error. The formula can then be simplified to:

e z = e x -e y + e
The formula comprises two elements: firstly the combination of deviations from input values and secondly, the deviation introduced by the elementary operation.

Addition : z = x ⊕ y → ez = ex + ey + e⊕ Subtraction : z = x y → ez = ex -ey + e M ultiplication : z = x ⊗ y → ez = x F ey + y F ex + exey + e⊗ Division : z = x y → ez = y F ex -x F ey y F (y F + ey) + e Fig. 1

. Computation of deviation for basic operations

Figure 1 formulates computation deviations for all four basic operations. For each of these formulae, the error computation combines deviations from input values and the error introduced by the current operation. Notice that, for multiplication and division, this deviation is proportional to input values.

All these formulae compute the difference between the expected result on the reals and the actual one on the floats for a basic operation. Our constraint solver over the errors on the floats relies on these formulae.

Domain of errors

In a classical CSP, to each variable x is associated x its domain of values. It denotes the set of possible values that this variable can take. When the variable takes values in F, its domain of values is represented by an interval of floats:

x F = [x F , x F ] = {x F ∈ F, x F ≤ x F ≤ x F } where x F ∈ F and x F ∈ F.
Computation errors form a new dimension to consider. They require a specific domain in view of the distinct nature of elements to represent, but also, due to the possible values of errors which belong to the set of reals. Therefore, we introduce a domain of errors, which is associated with each variable of a problem. Since all arithmetic constraints processed here are reduced to the four basic operations, and since those four operations are applied over floats, i.e. a finite subset of rationals, this domain can be defined as an interval of rationals with bounds in Q:

e x = [e x , e x ] = {e x ∈ Q, e x ≤ e x ≤ e x }
where e x ∈ Q and e x ∈ Q.

Another domain of errors is required for the smooth running of our system: it is the domain of errors on operations, denoted by e , that appears in the computation of the deviations (see Figure 1). Contrary to previous domains, it is not attached to each variable of a problem but to each instance of an arithmetic operation of a problem.

Like the domain of errors attached to a variable, it takes values in the set of rationals. Thus, we have: This triple, composed of the domain of values, the domain of errors, and the domain of errors on operations, is required to represent the set of phenomena due firstly to possible values on variables and secondly, to different errors that come into play in computation over floats.

Projection functions

The filtering process of our solver is based on classical projection functions to reduce the domains of variables. Domains of values can be computed by projection functions defined in [START_REF] Michel | Exact projection functions for floating point number constraints[END_REF] and extended in [START_REF] Botella | Symbolic execution of floating-point computations[END_REF] and [START_REF] Marre | Improving the floating point addition and subtraction constraints[END_REF] but new ones are required for the domains of errors.

Those projections on the domains of errors are made through an extension over intervals of formulae from Figure 1. Since these formulae are written over reals, they can naturally be extended to intervals. For example, in the case of the subtraction, we get the four projection functions below: where e x , e y , and e z are the domains of errors of variables x, y, and z respectively and e is the domain of errors on the subtraction.

Figure 2 gives projection functions for the three other arithmetic operations. Note that the projection on y F for the division requires solving a quadratic equation and requires the computation of a square root. Thanks to outward roundings, correct computation of such a square root on rationals is obtained using floating-point square root at the price of an over-approximation. Note also that these projections handle the general case. For the sake of simplicity, special cases like a division by zero are not exposed here.

Addition : ez ← ez ∩ (ex + ey + e⊕) ex ← ex ∩ (ez -ey -e⊕) ey ← ey ∩ (ez -ex -e⊕) e⊕ ← e⊕ ∩ (ez -ex -ey) M ultiplication : ez ← ez ∩ (x F ey + y F ex + exey + e⊗) ex ← ex ∩ ez -x F ey -e⊗ y F + ey ey ← ey ∩ ez -y F ex -e⊗ x F + ex e⊗ ← e⊗ ∩ (ez -x F ey -y F ex -exey) x F ← x F ∩ ez -y F ex -exey -e⊗ ey y F ← y F ∩ ez -x F ey -exey -e⊗ ex Division : ez ← ez ∩ y F ex -x F ey y F (y F + ey) + e ex ← ex ∩ (ez -e )(y F + ey) + x F ey y F ey ← ey ∩ ex -ezy F + e y F ez -e + x F y F e ← e ∩ ez - y F ex -x F ey y F (y F + ey) x F ← x F ∩ (e -ez)y F (y F + ey) + y F ex ey y F ← y F ∩ [min(δ 1 , δ 2 ), max(δ1, δ2)] with δ1 ← ex -(ez -e )ey - √ ∆ 2(ez -e ) δ2 ← ex -(ez -e )ey + √ ∆ 2(ez -e ) ∆ ← [0, +∞) ∩ ((ez -e )ey -ex) 2 + 4(ez -e )eyx F

Fig. 2. Projection functions of arithmetic operation

Projection functions, on the domain of errors, support only arithmetic operations and assignment, where the computation error from the expression is transmitted to the assigned variable. Since the error is not involved in comparison operators, their projection functions only manage domains of values.

The set of those projection functions is used to reduce all variables' domains until a fixed point is reached. For the sake of efficiency, but also to get around potential slow convergence, the fixed point computation is stopped when no domain reduction is greater than 5%.

Links between the domain of values and the domain of errors

In order to take advantage of domain reductions of one domain in another domain, clear and strong links must be established between the domain of values and the domain of errors. What naturally occurs for domains of values thanks to constraints on values requires more attention when it comes to the relations between dual domains. A first relation between the domain of values and the domain of errors on operations is based upon the IEEE 754 norm, which guarantees that basic arithmetic operations are correctly rounded. Since the four basic operations are correctly rounded to the nearest even float, we have

(x y) - (x y) -(x y) - 2 ≤ (x • y) ≤ (x y) + (x y) + -(x y) 2
where x -and x + denote respectively, the greatest floating-point number strictly smaller than x and the smallest floating-point number strictly larger than x. In other words, the result over floats is, at a half-ulp, the distance between two successive floats from the result over reals. Thus, the error on an operation is contained in this ulp:

- (x y) -(x y) - 2 ≤ e ≤ + (x y) + -(x y) 2
This equation sets a relation between the domain of values and the domain of errors on operations: operation errors can never be greater than the greatest half-ulp of the domain of values on the operation result. The projection function for the domain of errors on operations is obtained by extending this formula to intervals:

e ← e ∩ - min((z -z -), (z -z -)) 2 , + max((z + -z), (z + -z)) 2
Finally, these links are refined by means of other well-known properties of floating-point arithmetic like the Sterbenz property of the subtraction [START_REF] Sterbenz | Floating Point Computation[END_REF] or the Hauser property on the addition [START_REF] Hauser | Handling floating-point exceptions in numeric programs[END_REF]. Both properties give conditions under which these operations produce exact results. As is the well-known property that 2 k * x is exactly computed provided that no overflow occurs.

Constraints over errors

Usually, constraints available in a solver establish relations between variables of a problem. The duality of domains available in our solver requires introducing a distinction between the domain of values and the domain of errors. In order to preserve the current semantic of expressions, variables keep on representing possible values. A dedicated function, err(x), makes it possible to express constraints over errors. For example, abs(err(x)) ≥ , denote a constraint which demands that the error on variable x be, in absolute value, greater or equal to . It should be noted that since errors are taking their values in Q, the constraint is over rationals.

When a constraint involves errors and variables, the latter, with domains over floats, are promoted to rationals. Therefore, the constraint is converted to a constraint over rationals.

Preliminary experiments

Projection functions and constraints over errors are being evaluated in a prototype based on Objective-CP [START_REF] Hentenryck | The objective-cp optimization system[END_REF], which already handles constraints over floats thanks to the projection functions of FPCS [START_REF] Michel | Solving constraints over floating-point numbers[END_REF]. All experiments are carried out on a MacBook Pro i7 2,8GHz with 16GB of memory.

Predator prey

Predator prey [START_REF] Darulova | Sound compilation of reals[END_REF] has been extracted from the FPBench test suite2 : double predatorPrey(double x) { double r = 4.0; double K = 1.11; double z = (((r * x) * x) / (1.0 + ((x / K) * (x / K)))); return z; } When x ∈ [ 1 10 , 3 10 ], Objective-CP using a simple and single filtering process reduces the domain of z to [3.72770587e-02, 3.57101682e-01] and its error domain to [-1.04160431e-16, 1.04160431e -16]. Fluctuat reduces the domain of z to [3.72770601e-02, 3.44238968e-01] and its error to [-1.33729201e-16, 1.33729201e-16]. Thus, while Fluctuat provides better bounds for the domain of values, Objective-CP provides better bounds for the domain of error. Moreover, our solver can use a search procedure to compute error values that are reachable. For example, we can search for input values such that the error on z will be strictly greater than zero.

With this constraint the solver output z = 3.354935286988540155128646e-01 with an error of 3.096612314293906301816432e-17.

Since rational numbers are used for computations of errors it is crucial to take solving time into account. Table 1 shows times in seconds for the generation of round-off error bounds on some benchmarks from FPBench. Solving time for Objective-CP are correct, especially as a search procedure is done in addition of filtering. Those times comfort us in the use of rational numbers for representing errors.

Conclusion

In this paper, we introduced a constraint solver capable of reasoning over computation errors on floating-point numbers. It is built over a system of dual domains, the first one characterizing possible values that a variable of the problem can take and the second one defining errors committed during computations. Moreover, there are particular domains, bind to instances of arithmetic operations in numerical expressions of the constraints, which represent errors in those operations. Our solver, enhanced with projection functions and constraints over errors, offers unique possibilities to reason on computation errors. Preliminary experiments are promising and will naturally be reinforced with more benchmarks.

Such a solver might appear limited by the use of rational numbers and multiprecision integers. However, a thorough examination of the solver behavior has shown that its main limit lies in its approximation of round-off errors. Firstly, round-off errors are not uniformly distributed across input values. As a result, finding input values that satisfy some error constraints has often to resort to an enumeration of possible values. Secondly, round-off error of operations are overestimated. Such an overestimation does not perform the fine domain reductions that would allow an efficient search. Therefore, a deeper understanding and a tighter representation of the round-off error on an operation basis is a must to actually improve the behavior of our solver.

Further work include extending support for a wider set of arithmetic functions, improving the search to quickly find solutions in the presence of constraints over errors and to add global optimization capabilities, for example, by using a branch-and-bound method. Then, by formulating a problem as an optimization problem, we should be in a position to determine for which input values the error is maximal.

Another direction of improvement is the combination of CSP with other tools dedicated to round-off error like abstract interpreter in an approach similar to what has already been done for domains of values [START_REF] Ponsini | Verifying floating-point programs with constraint programming and abstract interpretation techniques[END_REF].

e

  = [e , e ] = {e ∈ Q, e ≤ e ≤ e } where e ∈ Q and e ∈ Q.

  e z ← e z ∩ (e x -e y + e ) e x ← e x ∩ (e z + e y -e ) e y ← e y ∩ (-e z + e x + e ) e ← e ∩ (e z -e x + e y )

Table 1 .

 1 Times in seconds for the generation of round-off error bounds. For Objective-CP a searh is also used. (bold indicates the best approximation and italic indicates the second best)

	Gappa Fluctuat Real2Float FPTaylor PRECiSA Objective-CP
	carbonGas 0.152 0.025	0.815	1.209	3.830	0.060
	verhulst 0.034 0.043	0.465	0.812	0.789	0.032
	predPrey 0.052 0.031	0.735	0.916	0.477	0.050
	turbine1 0.165 0.028	67.960	2.906	110.272	0.232

By using a rational arithmetic library for computation on the reals and down to the memory limit.

See fpbench.org.


Caterina Calgaro
email: caterina.calgaro@univ-lille.fr

Claire Colin

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

A combined Finite Volumes -Finite

Keywords: Low-Mach model, Finite Volume method, Finite Element method, Projection scheme

HAL is

Introduction

Variable density -low Mach numbers flows have been widely studied in the recent literature because of their applicability in various phenomena such as flows in hightemperature gas reactors, meteorological flows, flows with convective and/or conductive heat transfert, combustion processes and many others. In such cases, the resolution of the full compressible Navier-Stokes system is not adapted, because of the sound waves speed which is much faster than the entropy or the vorticity ones. Consequently, this choice would impose a too strong time-step limitation in the framework of explicit solvers, leading to unreachable numerical simulations. The Boussinesq incompressible model is not a better alternative for such low-speed phenomena. Indeed, the compressibility effects can not be totally canceled because of large variations of temperature and density, even if pressure ones are much smaller. Consequently, some models have been formally derived, leading to the filtering of the acoustic waves by the use of some formal asymptotic expansions [START_REF] Majda | The derivation and numerical solution of the equations for zero Mach number combustion[END_REF][START_REF] Lions | Mathematical topics in fluid mechanics[END_REF][START_REF] Munz | The extension of incompressible flow solvers to the weakly compressible regime[END_REF][START_REF] Principe | Mathematical models for thermally coupled low speed flows[END_REF].

We recall that there exists two families of methods to compute flows at low-Mach number regime. On the one hand, there are the so-called density-based solvers, corresponding to methods used for the simulation of supersonic and transonic flows, which have been adapted to make them efficient and robust in the case of a low-Mach flow (see e.g. [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF][START_REF] Noelle | A weakly asymptotic preserving low Mach number scheme for the Euler equations of gas dynamics[END_REF][START_REF] Penel | Coupling strategies for compressible-low Mach number flows[END_REF][START_REF] Dellacherie | Construction of modified Godunov-type schemes accurate at any Mach number for the compressible Euler system[END_REF][START_REF] Herbin | Low Mach number limit of a pressure correction MAC scheme for compressible barotropic flows[END_REF]) and references therein, using for example some preconditionning techniques [START_REF] Turkel | Preconditioned methods for solving the incompressible and low speed compressible equations[END_REF][START_REF] Luo | Extension of harten-lax-van leer scheme for flows at all speeds[END_REF][START_REF] Vigneron | An implicit finite volume method for the solution of 3D low Mach number viscous flows using a local preconditioning technique[END_REF][START_REF] Tyliszczak | Application of time preconditioning and highorder compact discretization method for low Mach number flows[END_REF][START_REF] Volkov | Preconditioning of gas dynamics equations in compressible gas flow computations at low Mach numbers[END_REF]. On the other hand, there are the so-called pressure-based solvers, coming from the incompressible case. The pressure variations become independent from the state equation, and are coupled to the divergence condition on the velocity [START_REF] Dellacherie | On a diphasic low Mach number system[END_REF][START_REF] Dellacherie | On a low Mach nuclear core model[END_REF]. In that case, considering the non-linearity of the coupled system, a fixed-point iterations process is usually performed, see e.g. [START_REF] Beccantini | Numerical simulations of a transient injection flow at low Mach number regime[END_REF][START_REF] Gravemeier | An algebraic variational multiscale-multigrid method for large-eddy simulation of turbulent variable-density flow at low Mach number[END_REF][START_REF] Avila | A finite element dynamical nonlinear subscale approximation for the low Mach number flow equations[END_REF]. A fractional step method, initially developed in [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Chorin | On the convergence of discrete approximations to the Navier-Stokes equations[END_REF][START_REF] Temam | Navier-Stokes equations. Theory and numerical analysis[END_REF] and progressively improved in [START_REF] Guermond | Calculation of incompressible viscous flows by an unconditionally stable projection FEM[END_REF][START_REF] Guermond | On stability and convergence of projection methods based on pressure Poisson equation[END_REF][START_REF] Guermond | An overview of projection methods for 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] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF] is used most of the time for the momentum equation. The pressure field comes from the resolution of a Poisson equation, for which the right-hand side contains some derivatives of the density or the temperature. Concerning the space discretization, lots of papers deal with Finite Element methods [START_REF] Heuveline | On higher-order mixed FEM for low Mach number flows: application to a natural convection benchmark problem[END_REF][START_REF] Martinez | A finite element method for low-speed compressible flows[END_REF], which can also been stabilized in the case of convection dominated regimes [START_REF] Principe | A stabilized finite element approximation of low speed thermally coupled flows[END_REF][START_REF] Liu | An implicit finite element solution of thermal flows at low Mach number[END_REF][START_REF] Gravemeier | An algebraic variational multiscale-multigrid method for large-eddy simulation of turbulent variable-density flow at low Mach number[END_REF][START_REF] Gravemeier | Residual-based variational multiscale methods for laminar, transitional and turbulent variable-density flow at low Mach number[END_REF][START_REF] Avila | A finite element dynamical nonlinear subscale approximation for the low Mach number flow equations[END_REF]. Others are also devoted to Finite Differences [START_REF] Nicoud | Conservative high-order finite-difference schemes for low-Mach number flows[END_REF][START_REF] Dellacherie | On a diphasic low Mach number system[END_REF][START_REF] Lessani | Time-accurate calculation of variable density flows with strong temperature gradients and combustion[END_REF].

In this paper, we propose a combined Finite Volume -Finite Element method, which was initially developed for the simulation of incompressible and variable density flows [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF]. This method is based on a time splitting allowing to solve the mass conservation equation by a Finite Volume method, and the momentum equation associated with the free divergence constraint on the velocity by a Finite Element one. It allows in particular to preserve the constant density states and to ensure the discrete maximum principle [START_REF] Calgaro | L ∞ -stability of vertexbased MUSCL finite volume schemes on unstructured grids: simulation of incompressible flows with high density ratios[END_REF]. It also has been used to simulate some mixture flows such as avalanches [START_REF] Calgaro | Modeling and simulation of mixture flows: application to powder-snow avalanches[END_REF]. Following the same philosophy, we propose to adapt this method to the case of a low-Mach model, providing a new pressure-based solver. The originality of our approach lies in the fact that the density is computed from the mass equation with a Finite Volume method, the other variables of the problem being approximated by a Finite Element method. In our work, the equation of state in not explicitly imposed. Moreover, the scheme recovers the properties of our previous proposed scheme at the incompressible limit, namely the preservation of the constant states, as well as the discrete maximum principle for the density.

The paper is organized as follows. In Section 2 the governing equations are recalled, and the choice of the final system to be solved is justified among several equivalent sets of equations. In Section 3, the combined Finite Volumes -Finite Elements is carefully described. Section 4 is devoted to some numerical simulations to exhibit the ability of the code. First, some analytical benchmarks are proposed and underline the accuracy of the scheme. Then, a transient injection flow [START_REF] Beccantini | Numerical simulations of a transient injection flow at low Mach number regime[END_REF][START_REF] Gravemeier | Residual-based variational multiscale methods for laminar, transitional and turbulent variable-density flow at low Mach number[END_REF][START_REF] Avila | A finite element dynamical nonlinear subscale approximation for the low Mach number flow equations[END_REF] as well as the natural convection of a flow in a cavity [START_REF] Le Quéré | A Chebyshev collocation algorithm for 2d non-Boussinesq convection[END_REF][START_REF] Heuveline | On higher-order mixed FEM for low Mach number flows: application to a natural convection benchmark problem[END_REF][START_REF] Gravemeier | Residual-based variational multiscale methods for laminar, transitional and turbulent variable-density flow at low Mach number[END_REF][START_REF] Avila | A finite element dynamical nonlinear subscale approximation for the low Mach number flow equations[END_REF] are simulated.

Governing equations 2.1 Choice of the system

The equations modelling low Mach number flows are derived by inserting the asymptotic expansions of the variables with respect to the Mach number M in the Navier-Stokes compressible equations [START_REF] Majda | The derivation and numerical solution of the equations for zero Mach number combustion[END_REF][START_REF] Lions | Mathematical topics in fluid mechanics[END_REF][START_REF] Dellacherie | On a diphasic low Mach number system[END_REF]. One of the characteristics of the process is that the pressure splits into two terms. Denoting x ∈ R d the space variable and t ∈ R + * the time one, we write: p(x, t) = P (t) + π(x, t),

where P is called the thermodynamic pressure and π the dynamic pressure. P only depends on t, and π is in the order of M 2 . The other variables considered are the velocity u(x, t), the density ρ(x, t) and the temperature T (x, t). The non-dimensioned characteristic numbers of the flow are given by: Moreover, assuming that the state law of the fluid is close to the one of a perfect case (see assumption 2.5. of [START_REF] Dellacherie | On a diphasic low Mach number system[END_REF]), a reference temperature of the flow is also deduced by

T ref = P ref C P ref ρ ref .
Finally, given α ref a reference value for the compressibility coefficient at constant pressure, a non-dimensioned number β ref is defined by:

β ref = α ref T ref .
According to single time scale and single space scale asymptotics, the continuity, momentum and temperature equations in the non-dimensioned formulation in an open polygonal domain Ω ⊂ R d are given by (see [START_REF] Lions | Mathematical topics in fluid mechanics[END_REF][START_REF] Dellacherie | On a low Mach nuclear core model[END_REF], see also [START_REF] Dellacherie | On a diphasic low Mach number system[END_REF] adapted in case of a single-phase flow):

∂ρ ∂t + ∇ • (ρu) = 0, (1)
ρ ∂u ∂t + u • ∇u + ∇π - 1 Re ∇ • τ = - 1 Fr 2 ρ e y , (2)
ρ C P ∂T ∂t + u • ∇T -β ref α T dP dt - 1 Re Pr ∇ • (λ∇T ) = 0, (3)
where the density is obtained with the following equation of state:

ρ = ρ(T, P ). ( 4
)
Here, τ is the viscosity stress tensor defined by

τ = µ ∇u + ∇u T - 2 3 ∇ • u I ,
where µ(T, P ) is the fluid viscosity, λ(T, P ) the thermal conductivity and e y = (0, 1) T . The quantities α(T, P ) and C P (T, P ) are respectively the compressibility coefficient at constant pressure and the calorific capacity at constant pressure, in their nondimensioned formulation. They are given by:

α(T, P ) = - 1 β ref ρ ∂ ρ ∂T (T, P ), (5)
and

C P (T, P ) = ∂h ∂T (T, P ), (6)
where h = e + P ρ is the enthalpy of the fluid, and e = ẽ(T, P )

is the internal energy. We assume that the equations of state ( 4) and ( 7) satisfy the thermodynamic hypothesis 3.2. given in [START_REF] Dellacherie | On a diphasic low Mach number system[END_REF], namely that there exists a strictly convex function s( 1 ρ , e), defined as the entropy of the fluid, verifying :

-T ds = de + P d 1 ρ .

This last assumption allows us to derive the equation (3) from the energy equation. 4) in (1) and then [START_REF] Calgaro | Modeling and simulation of mixture flows: application to powder-snow avalanches[END_REF], we obtain:

Now, with the notation D

t = ∂ ∂t + u • ∇, using (
∇ • u = β ref α D t T - 1 ρ ∂ ρ ∂P dP dt .
Now the equation (3) takes the formulation:

∇ • u = - 1 ΓP dP dt + β ref β Re Pr P ∇ • (λ∇T ), ( 8
)
where Γ is defined by

Γ(T, P ) = ρ c 2 P (T, P ), (9)
c = ∂ ρ ∂P - β 2 ref α 2 T C P -1
being the sound velocity ; and β is defined by : β(T, P ) = α P ρ C P (T, P ). [START_REF] Dellacherie | On a diphasic low Mach number system[END_REF] Finally, the thermodynamic pressure time evolution is determined by integrating (8) over the domain Ω ⊂ R d :

dP dt Ω 1 Γ β ref β = -P Ω ∇ • u β ref β + 1 Re Pr ∂Ω λ ∇T • n, (11)
where n is the outer unit normal to Ω on ∂Ω.

In this work, we propose a new way to proceed, using the fact that the system (1)-( 2)-( 3)-( 4)-( 11) is equivalent to the system (1)-( 2)-( 3)-( 8) and [START_REF] Dellacherie | On a low Mach nuclear core model[END_REF]. In particular, we verify that the equation of state ( 4) is implicitly imposed. Indeed, starting from (1) and using [START_REF] Chorin | On the convergence of discrete approximations to the Navier-Stokes equations[END_REF] to subsitute the term ∇ • u, we get:

D t ρ = -ρ - 1 Γ P dP dt + β ref β Re Pr P ∇ • (λ∇T ) .
Now, using (3) to substitute the diffusive term gives :

D t ρ = ρ Γ P + β 2 ref ρ β α T P D t P - β ref ρ 2 β C P P D t T.
Finally, thanks to the expression of α, Γ and β respectively given in (5), ( 9) et [START_REF] Dellacherie | On a diphasic low Mach number system[END_REF], we obtain

D t ρ = ∂ ρ ∂P D t P + ∂ ρ ∂T D t T.
Consequently, supposing that the equation of state is initially fulfilled, we can deduce (4).

Perfect gases. In most of applications (see for exemple [START_REF] Heuveline | On higher-order mixed FEM for low Mach number flows: application to a natural convection benchmark problem[END_REF][START_REF] Le Quéré | Modelling of natural convection flows with large temperature differences: a benchmark problem for low Mach number solvers[END_REF][START_REF] Dellacherie | On a diphasic low Mach number system[END_REF][START_REF] Beccantini | Numerical simulations of a transient injection flow at low Mach number regime[END_REF][START_REF] Avila | A finite element dynamical nonlinear subscale approximation for the low Mach number flow equations[END_REF][START_REF] Gravemeier | Residual-based variational multiscale methods for laminar, transitional and turbulent variable-density flow at low Mach number[END_REF]), calorifically perfect gases are considered and we also restrict the rest of this work to the study of such fluids. For these gases, the equation of state in the dimensionless form 1 is given by

ρ = γ P (γ -1)T . ( 12
)
Consequently,

α ref = 1 T ref
is a natural choice, and we get β ref = 1. Thus we obtain:

C P = 1, β = γ -1 γ , c = (γ -1)T , Γ = γ.
For the sake of simplicity, we also assume that the viscosity µ and the thermal conductivity λ are constant, namely λ = µ = 1. We thus obtain the Navier-Stokes equations at low Mach number regime for a calorifically perfect gas:

∂ρ ∂t + ∇ • (ρu) = 0, (13)
ρ ∂u ∂t + u • ∇u + ∇π - 1 Re ∇ • τ = - 1 Fr 2 ρ e y , (14)
ρ ∂T ∂t + u • ∇T - dP dt - 1 Re Pr ∆T = 0, (15)
∇ • u = - 1 γP dP dt + γ -1 γ Re Pr P ∆T, (16)
dP dt = - γ P |Ω| ∂Ω u • n + γ -1 |Ω| Re Pr ∂Ω ∇T • n, (17)
with

τ = ∇u + ∇u T - 2 3 ∇ • u I.
1 Denoting with star superscripts the dimensioned variables, the equations of states are defined with

ρ * = P * R T * and e * = R T * γ -1
, where γ = 1.4 and R = 287Jkg -1 K -1 are the gas specific heat ratio and the gas constant, respectively. Then, we obtain

C * P = γR γ -1 and α * (T * ) = 1 T * .

Initial and boundary conditions

The initial conditions for the system ( 13)-( 14)-( 15)-( 16) and ( 17) are given in Ω by: u(x, 0) = u 0 (x), T (x, 0) = T 0 (x), P (0) = P 0 , ρ(x, 0) = γ P 0 (γ -1)T 0 (x) .

The boundary conditions on the velocity are given by:

u = u D on ∂Ω,
where u D is a function to be specified. We set

∂Ω = ∂Ω D ∪ ∂Ω N with ∂Ω D ∩ ∂Ω N = ∅.
We define n as the outward unit normal of Ω and ∂Ω in the part of ∂Ω such that u D • n < 0. The boundary conditions on the temperature and density are respectively given by:

T = T D on ∂Ω D , ∇T • n = F N on ∂Ω N and ρ = ρ in on ∂Ω in ,
with T D , F N and ρ in some functions to be specified.

The combined Finite Volumes -Finite Elements method

The combined Finite Volumes -Finite Elements scheme (C-FV-FE scheme) is based on a time splitting, allowing in particular to solve equation ( 13) by a finite volumes solver, and equations ( 14), ( 15), ( 16) by a finite elements one, using the same strategy as the one developed in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF] for the variable density incompressible Navier-Stokes system. We first describe the splitting used in time in order to reach the globally optimal accuracy in time according to the accuracy in time of each part of the scheme.

The time splitting

Let ∆t be the time step and t n = n∆t. Functions approximated at time t n will be identified with superscript n. We assume that T n-1 , ρ n-1 and u n-1 as well as T n , ρ n and u n are known approximated values.

1. Since we want at least the global mass conservation in the domain, the thermodynamic pressure P n+1 is computed by integrating the equation of state over Ω, using an extrapolation of the density and the temperature:

P n+1 = γ -1 γ Ω 2ρ n -ρ n-1 Ω 1 2T n -T n-1 -1 . ( 18
)
2. The time derivative of the thermodynamic pressure dP dt n+1 is obtained through the use of equation ( 17) and an extrapolation of the temperature:

dP dt n+1 + γP n+1 |Ω| ∂Ω u n+1 D .n = γ -1 |Ω|Re Pr ∂Ω D ∇(2T n -T n-1 ) • n + γ -1 |Ω|Re Pr ∂Ω N F n+1 N . (19)
3. In order to ensure a local mass conservation in the domain, the density ρ n+1 at time t n+1 is computed by solving the continuity equation ( 13) using a second order Runge-Kutta scheme in time:

ρn+1 -ρ n ∆t + ∇ • (ρ n u n+ 1 2 ) = 0, ( 20
)
ρ n+1 -ρ n ∆t + 1 2 ∇ • (ρ n u n+ 1 2 ) + ∇ • (ρ n+1 u n+ 1 2 ) = 0, (21)
with

u n+ 1 2 = (2u n -u n-1 ) + u n 2 = 3u n -u n-1 2 (22)
and

ρn+1 = ρ n+1 = ρ n+1 in on ∂Ω in . (23)
4. The temperature T n+1 is computed by solving the temperature equation ( 15) using a BDF2-scheme in time:

ρ n+1 3T n+1 -4T n + T n-1 2 ∆t + (2u n -u n-1 ) • ∇T n+1 - 3P n+1 -4P n + P n-1 2 ∆t - 1 Re Pr ∆T n+1 = 0 (24)
and

T n+1 = T n+1 D on ∂Ω D , (25)
∇T n+1 • n = F n+1 N on ∂Ω N . ( 26
)
5. The velocity u n+1 and the pressure π n+1 are computed by solving the momentum equation ( 14) associated with the compressibility constraint ( 16) using a BDF2scheme in time:

ρ n+1 3u n+1 -4u n + u n-1 2 ∆t + (2u n -u n-1 ) • ∇u n+1 + ∇π n+1 (27)
- 1 Re ∇ • τ n+1 = - 1 Fr 2 ρ n+1 e y , ∇ • u n+1 = - 1 γ P n+1 dP dt n+1 + γ -1 γ Re Pr P n+1 ∆(2T n -T n-1 ), (28)
with

∇ • τ n+1 = ∆u n+1 + 1 3 ∇ ∇ • u n+1 (29)
and

u n+1 = u n+1 D on ∂Ω.
Similarly to other references [START_REF] Guermond | Un résultat de convergence d'ordre deux en temps pour l'approximation des équations de Navier-Stokes par une technique de projection incrémentale[END_REF][START_REF] Guermond | Error analysis of a fractional time-stepping technique for incompressible flows with variable density[END_REF] or to our previous contributions for incompressible fluids [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF][START_REF] Calgaro | Modeling and simulation of mixture flows: application to powder-snow avalanches[END_REF][START_REF] Calgaro | Simulations of non homogeneous viscous flows with incompressibility constraints[END_REF], an extrapolation was used to approximate the convective velocity at time t n+1 by 2u nu n-1 in ( 24) and [START_REF] Lessani | Time-accurate calculation of variable density flows with strong temperature gradients and combustion[END_REF]. Moreover, the velocity used in ( 20) and ( 21) is extrapolated at time (t n + t n+1 )/2 as indicated in [START_REF] Guillard | On the behaviour of upwind schemes in the low Mach number limit[END_REF], which is necessary to reach the second order accuracy. Note that in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF], alternatively to this extrapolation process, the second order in time was obtained thanks to a Strang splitting. Finally, an extrapolation of the temperature (and eventually of the density) was done in ( 18), ( 19) and [START_REF] Lions | Mathematical topics in fluid mechanics[END_REF]. The choice 2T n -T n-1 in ( 19) and ( 28) has to be the same for compatibility reasons, [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF] being obtained by integrating [START_REF] Lions | Mathematical topics in fluid mechanics[END_REF] over Ω. We underline that with these choices, we did not observe any numerical instabilities. Following the same ideas as in [START_REF] Beccantini | Numerical simulations of a transient injection flow at low Mach number regime[END_REF][START_REF] Gravemeier | An algebraic variational multiscale-multigrid method for large-eddy simulation of turbulent variable-density flow at low Mach number[END_REF][START_REF] Avila | A finite element dynamical nonlinear subscale approximation for the low Mach number flow equations[END_REF], we also implemented a fixed-point method for the C-FV-FE scheme (see Appendix A for the description of the corresponding fixed-point algorithm), to evaluate its impact on some configurations (see section 4). Remark 3.1. As it will be explained in the next section, such a way to proceed allows in particular to solve the mass equation by a FV method. Consequently, it ensures a local mass conservation, while preserving the maximum principle on the density at the incompressible limit of the fluid. In that sense, the proposed numerical scheme can be seen as a generalization of the one previously developed for the variable density incompressible model [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF].

Remark 3.2. Instead of solving equations ( 18) and (19), we could compute P n+1 by solving the ordinary differential equation (17) using a BDF2-scheme in time, and approximate dP dt n+1 in (28) as made in [START_REF] Heuveline | On higher-order mixed FEM for low Mach number flows: application to a natural convection benchmark problem[END_REF]. However, although this strategy can be used for the simulation of evolutive smooth solutions, according to [START_REF] Le Quéré | A Chebyshev collocation algorithm for 2d non-Boussinesq convection[END_REF] we observe that for steady states this scheme can prevent the algorithm from reaching any stationary solution.

Remark 3.3. Instead of solving the temperature equation (24), we could have in mind to use the state equation (12) to compute T n+1 . Nevertheless, this would lead to a lack of regularity of T n+1 , which is necessary in the constraint (28) in order to derive the values of u n+1 and π n+1 .

Space discretization

The discretization in space is based on a triangulation of the domain Ω ⊂ R 2 by a set of triangles defining a regular mesh τ h in the Ciarlet sense [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. Each component of the velocity u h is discretized by some P 2 -Lagrange finite elements, and the pressure π h by some P 1 -Lagrange finite elements, leading to H 1 (Ω) conforming approximations fulfilling the usual discrete LBB condition and consequently ensuring the stability of the discrete problem. The temperature T h is also discretized by some P 2 -Lagrange finite elements, leading to a H 1 (Ω) conforming approximation. The density ρ h is discretized by piecewise constant values on a dual mesh τ * h associated with τ h , allowing to consider a vertex-based finite-volume schemes for the resolution of the mass equation. The density field can also be interpreted as a P 1 -Lagrange finite elements field, since a value of the density is naturally associated with each node of any triangles. The degrees of freedom of each variable corresponding to the associated space discretization are displayed in Figure 1. For further details, we refer to [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF]. are simply computed by using the discrete versions of ( 18) and ( 19), namely:

Primal mesh τ h

Dual mesh τ *

P n+1 = γ -1 γ Ω 2ρ n h -ρ n-1 h Ω 1 2T n h -T n-1 h -1 (30)
and

dP dt n+1 + γP n+1 |Ω| ∂Ω u n+1 D .n = γ -1 |Ω|Re Pr ∂Ω D ∇(2T n h -T n-1 h ) • n + γ -1 |Ω|Re Pr ∂Ω N F n+1 N . (31)
3.2.2 Solving the temperature by a FE method

Assuming that ρ n+1 h , T n h , T n-1 h
, ũn h , ũn-1 h , P n+1 , P n and P n-1 are known, the value of T n+1 h is computed by the resolution of equation ( 24), namely:

ρ n+1 h 3T n+1 h -4T n h + T n-1 h 2 ∆t + (2ũ n h -ũn-1 h ) • ∇T n+1 h - 3P n+1 -4P n + P n-1 2 ∆t - 1 Re Pr ∆T n+1 h = 0. ( 32
)
It is performed considering its weak finite elements formulation, associated with the boundary conditions specified in ( 25)-( 26).

Solving the velocity by a FE method

We detail here the projection method used to derive u n+1 h and π n+1 h , contrary to [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF] in which an Uzawa solver was considered. It constitutes a natural adaptation from the section 4 of [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF] to the low-Mach model case. We assume that ρ n+1 h , π n h , P n+1 , P n , P n-1 and T n+1 h are known, as well as the auxiliary variables ũn h , ũn-1 h , φ n h , and φ n-1 h specific to the projection method (where (ũ 0 h , φ 0 h ) and (ũ 1 h , φ 1 h ) are initialized in the same way as in [START_REF] Guermond | A splitting method for incompressible flows with variable density based on a pressure Poisson equation[END_REF]). First, the velocity field ũn+1 h which does not fulfill the constraint ( 16) is computed by solving the weak finite elements formulation of the parabolic equation, based on a BDF2 integration scheme, similarly to [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]:

         ρ n+1 h 3ũ n+1 h -4ũ n h + ũn-1 h 2∆t + (2ũ n h -ũn-1 h ) • ∇ũ n+1 h + ∇ π n h + 4 3 φ n h - 1 3 φ n-1 h - 1 Re ∇ • τh n+1 = - 1 Fr 2 ρ n+1 h e y , ũn+1 h | ∂Ω = u n+1 D , ( 33
) with ∇ • τ n+1 h = ∆ũ n+1 h + 1 3 ∇ ∇ • (2ũ n h -ũn-1 h ) . (34)
Comparing ( 29) and [START_REF] Nicoud | Conservative high-order finite-difference schemes for low-Mach number flows[END_REF], we notice that an extrapolation is used because of the projection method, since the two velocity components are computed successively. Then, the pressure π h n+1 is defined by:

π n+1 h = π n h + φ n+1 h .
Here, φ n+1 h is the solution of the weak finite elements formulation of the elliptic equation given by:

               ∇ • 1 ρ n+1 h ∇φ n+1 h = 3 2∆t ∇ • ũn+1 h + 1 γP n+1 dP dt n+1 - (γ -1) γ Re Pr P n+1 ∆(2T n h -T n-1 h ) , ∇φ n+1 h • n| ∂Ω = 0, (35)
where for any triangle K ∈ τ h we define:

1 ρ n+1 h | K = 1 3
A i ∈K 1 ρ n+1 h (A i ) , (36)
with A i the three vertices belonging to the triangle K.

Remark 3.4. Note that the choice of

ρ n+1 h | K given in (36) is not unique. Another possible definition for ρ n+1 h | K is 1 ρ n+1 h | K = 3 A i ∈K ρ n+1 h (A i )
and leads to similar results. The essential point is to define a constant density per triangle identically in [START_REF] Noelle | A weakly asymptotic preserving low Mach number scheme for the Euler equations of gas dynamics[END_REF], [START_REF] Principe | A stabilized finite element approximation of low speed thermally coupled flows[END_REF] and [START_REF] Turkel | Preconditioned methods for solving the incompressible and low speed compressible equations[END_REF], see remark 3.6.

Finally u n+1

h can be defined by:

u n+1 h = ũn+1 h - 2 ∆t 3 ρ n+1 h ∇φ n+1 h , (37)
even if there is no need to evaluate it from the practical point of view, as mentioned in [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF].

Solving the density with a FV method

Assuming that ρ n h , ρ n-1 h , ũn h , ũn-1 h , φ n h and φ n-1 h are known, the value of ρ n+1 h is computed by the resolution of the discrete version of equations ( 20)-( 21), namely:

ρn+1 h -ρ n h ∆t + ∇ • (ρ n h u * ,n+ 1 2 h ) = 0, ρ n+1 h -ρ n h ∆t + 1 2 ∇ • (ρ n h u * ,n+ 1 2 h ) + ∇ • (ρ n+1 h u * ,n+ 1 2 h ) = 0, (38)
where

u * ,n+ 1 2 h | K = 3u * ,n h | K -u * ,n-1 h | K 2 ( 39
)
is defined by

u * ,n h | K = 1 |K| K ũn h - 2 ∆t 3 ρ n h ∇φ n h . (40)
The Finite-Volume method is carefully detailed in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF], and its generalization to ensure the L ∞ -stability in the case of incompressible flows with the use of the so-called τ -limiters is given in [START_REF] Calgaro | L ∞ -stability of vertexbased MUSCL finite volume schemes on unstructured grids: simulation of incompressible flows with high density ratios[END_REF]. Here, we point out the fact that from the values of the velocity obtained by the Finite Element scheme, we need to deduce the values of an auxiliary velocity u * ,n+ 1 2 h at the interfaces of the density control volumes surrounding each node of the triangulation. These interfaces correspond to the dotted lines of the dual mesh τ * h displayed in Figure 1. Following the same strategy as in the incompressible case [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF], this value has to be piecewise constant on each triangle K of the mesh τ h . Remark 3.5. Definition (39) of u * ,n+ 1 2 h | K allows to ensure that in the case of a constant flow density in space (and consequently also constant in temperature), the scheme preserves the constant states imposed in the continuous model by the incompressibility constraint. In other words, the weak divergence property in the sense of the Finite Element projection method is transferred to the Finite Volume method, as it was proved in the context of a direct resolution by an Uzawa solver in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF].

Remark 3.6. In the incompressible case, let us note that the developed scheme allows to preserve the constant states. Indeed, let assume the following properties: -

Ω A 1 ρ n h ∇φ n h • ∇ψ A = 3 2∆t Ω A ∇ • ũn h ψ A + 1 γP n dP dt n Ω A ψ A + γ -1 γ Re Pr P n Ω A ∇(2T n-1 h -T n-2 h ) • ∇ψ A ,
where ψ A is the P 1 basis function associated to an internal node A and Ω A is the support of ψ A (see Figure 1). Assumptions (H1), (H2) and (H4) lead to:

Ω A ũn h - 2∆t 3 ρ n h ∇φ n h • ∇ψ A = 0,
and definition (37) of u n h gives:

Ω A u n h • ∇ψ A = 0.
It is the analogous of relation ( 23) in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF], which allows to obtain:

∂C A u * ,n h • n = 0,
where C A is the FV control volume associated to node A. Then, (H3) leads to ρ n+1 h = ρ n h . Finally, equation [START_REF] Martinez | A finite element method for low-speed compressible flows[END_REF], assumptions (H1) and (H2) associated to the previous results give

T n+1 h = T n h .
Note that ( 32) and ( 33) are both implicit schemes in time for the resolution of parabolic equations, so that they are not constrained by a CFL condition. The only explicit scheme in time is applied to the continuity equation [START_REF] Principe | Mathematical models for thermally coupled low speed flows[END_REF], so that the time step is constrained by a CFL condition and proportional to h/||u|| L ∞ , see [START_REF] Calgaro | L ∞ -stability of vertexbased MUSCL finite volume schemes on unstructured grids: simulation of incompressible flows with high density ratios[END_REF]. In practice, if ∆t does not satisfy this CFL condition, a smaller time step ∆t F V satisfying the CFL condition is defined and sub-steps are made in the finite volume scheme to compute the density.

Numerical simulations 4.1 Analytical Benchmarks

Constant states

We first want to illustrate the fact that for incompressible flows, if the density and the temperature are initially homogeneous, then they remain constant (see Remark 3.6). It was already done for the density in the case of the variable density incompressible system using an Uzawa solver (see [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF]). To do this, we consider the following analytical solution:

                   u ex (x, y) = 4 -y(x -1) 2 (x + 1) 2 (y -1)(y + 1) x(y -1) 2 (y + 1) 2 (x -1)(x + 1) , ρ ex = 1, T ex = γ γ -1 P ex = 1, π ex = - y F r 2 , (41)
in the square domain Ω = [-1, 1] 2 . A source term in the right-hand-side of equation ( 14) is consequently added to the gravity term. Dirichlet boundary conditions are prescribed for the temperature on the whole boundary (i.e. ∂Ω = ∂Ω D ). Since ∂Ω in = ∅, there is no need to specify the value of ρ at the boundary. Simulations are performed on unstructured meshes at Re = 1, Pr = γ and 1/Fr 2 = 9.81 up to time t f = 1 = N ∆t with ∆t = h max = 0.0625 and h max the maximum space step of the mesh. We give in Table 1 the values of max

0≤n≤N ||ρ n h -ρ ex || L ∞ (Ω) and max 0≤n≤N ||T n h -T ex || L ∞ (Ω) . max 0≤n≤N ||ρ n h -ρ ex || L ∞ (Ω) max 0≤n≤N ||T n h -T ex || L ∞ (Ω) 2.91e-13
1.02e-13 We observe on Table 1 that the density and the temperature remain constant nearly to the machine error during the whole simulation time, as it can be seen at the final time t f in Figure 2 on the mesh corresponding to h max = 0.0625. The density local maximum error is located in the vicinity of the boundaries. Consequently, the scheme ensures the preservation of the constant values of density and temperature in the presence of a gravity field, in the case where the velocity field is divergence free.

Density

Temperature Pressure

Horizontal velocity Vertical velocity

Analytical solution

Now, in order to investigate the accuracy of the scheme, convergence tests are performed for a smooth solution and a non-solenoidal velocity field. The analytical solution is given by:

                       u ex (t, x, y) = - 1 γ(2 + sin(2πt)) 2(γ -1)(2 + cos(2πt)) (1 + x 2 + y 2 ) 2 + π cos(2πt) x y , ρ ex (t, x, y) = γ(2 + sin(2πt)) (γ -1)(2 + cos(2πt)) (1 + x 2 + y 2 ),
T ex (t, x, y) = 2 + cos(2πt) 1 + x 2 + y 2 , P ex (t) = 2 + sin(2πt), π ex (t, x, y) = sin(x) sin(y) sin(2πt), ( 42) in the square Ω = [-1, 1] 2 . The Reynolds and Prandtl numbers are both equal to 1, the gravity term is not considered and appropriated source terms are added in the right-hand-sides of equations ( 13), ( 14) and [START_REF] Guermond | Un résultat de convergence d'ordre deux en temps pour l'approximation des équations de Navier-Stokes par une technique de projection incrémentale[END_REF]. Non-homogeneous Neumann boundary conditions on the temperature are prescribed on the whole boundary of Ω (i.e. ∂Ω = ∂Ω N ), so that F N = ∇T ex • n on ∂Ω ; and ρ in = ρ ex on ∂Ω in = ∂Ω.

The simulations are performed on unstructured grids in the range 1/150 ≤ h max ≤ 1/50, up to the final time t f = 0.2, using ∆t = h max . In order to evaluate the performances of the splitting proposed in section 3.1, the C-FV-FE scheme results are compared with the ones obtained using some fixed-point iterations, see Appendix A. In that case, the fixed-point iterations are performed until the L 2 -norm of two successive iterates is smaller than 10 -10 for all variables. From the practical point of view, we observe that the fixed point converges in 6 or 7 iterations for the worst cases. Finally, these results are also compared to those obtained with one of the schemes proposed in [START_REF] Beccantini | Numerical simulations of a transient injection flow at low Mach number regime[END_REF] (pressure-based solver, asymptotic approach 1), for which some fixed-point iterations are required.

We observe, whatever the considered scheme, that the thermodynamic pressure P converges at order 2 in the L ∞ (0, t f ) norm. We plot in Figure 3 the L ∞ (0, t f ; L 2 (Ω)) norm of the errors on the density ρ h , the temperature T h , the velocity u h and the dynamic pressure π h as a function of h max in a log/log scale. On the one hand, all schemes provide a convergence rate at order 2 for the density, the temperature and the velocity. Concerning the dynamic pressure, a rate slightly larger than 1.5 is obtained. More precisely, for the C-FV-FE scheme, a rate between 1.65 and 1.85 is observed, as we can see in Table 2. These results are in good agreement with the incompressible constant density case, for which it was proved in [START_REF] Guermond | On the error estimates for the rotational pressurecorrection projection methods[END_REF] that the errors in time in the L 2 (Ω)-norm for the velocity and the pressure are of order 2 and 3/2 respectively. In fact, the numerical convergence rate obtained for the dynamic pressure is slightly better than the theorical expected one. On the other hand, we can see that the results obtained using the C-FV-FE scheme with or without a fixed-point iterations procedure correspond to the same orders of convergence. In particular, the errors are quite the same, except a small difference in the density error. In conclusion, the fixed-point iterations are not necessary for the C-FV-FE scheme on this analytical solution benchmark to obtain optimal convergence orders. We also point out that the fixed-point iterations are crucial for the scheme proposed in [START_REF] Beccantini | Numerical simulations of a transient injection flow at low Mach number regime[END_REF]. Indeed, without the fixed-point iterations, the scheme does not converge for the dynamic pressure, and orders of convergence are smaller than one for the velocity and the temperature.

In Figure 4, we plot the L ∞ (0, t f ; L 2 (Ω)) error on the discrete state equation, defined by:

err h = max 0≤n≤N ρ n h - γ P n (γ -1)T n h 2 L 2 (Ω),
obtained with the C-FV-FE scheme. As explained in section 2.1, the state equation is imposed implicitly. As we can see, err h converges towards zero at order 2, which corresponds to the expected behaviour because of the previous convergence rates obtained in ρ h , T h and P . This benchmark was initially proposed in [START_REF] Beccantini | Numerical simulations of a transient injection flow at low Mach number regime[END_REF] and also considered in [START_REF] Gravemeier | Residual-based variational multiscale methods for laminar, transitional and turbulent variable-density flow at low Mach number[END_REF] and [START_REF] Avila | A finite element dynamical nonlinear subscale approximation for the low Mach number flow equations[END_REF]. Considering the non-dimensioned equations, the domain is defined by a rectangle Ω = [-1.5 ; 1.5] × [0 ; 7] (see Figure 5(a)), defining a cavity in which a calorifically perfect gas is initially at rest. The inital temperature and thermodynamic pressure values are given by: T 0 = 300 and P 0 = 1.

(

) 43
The Reynolds, the Prandtl and the Froude numbers of the fluid are respectively equal to Re = 40, Pr = 0.71 and Fr = 0.042. Zero Dirichlet boundary conditions for the velocity and zero Neumann ones for the temperature are specified on all boundaries, except for a small hole in the bottom wall defined by

∂Ω D = ∂Ω in = - l 2 ; l 2 × {0}, with l = 0.2: u(t, x, y)| ∂Ω N = 0, ∇T • n| ∂Ω N = 0.
On ∂Ω D , the fluid is injected at the temperature T D = 600, subject to a parabolic inflow profile. The momentum is imposed on ∂Ω in from t = 0 up to t f = 6 by:

(ρ u) i n(x, y) = 0 ; 6 × 2.87 • 10 -3 γ ṁi n l 2 (γ -1) l 2 4 -x 2 T for 0 ≤ t ≤ 6, (44)
where ṁi n = 1 is the average momentum. Because of this kind of boundary condition, after the computation of P n+1 we update the velocity as follows:

u n+1 D = 0, 2.87 • 10 -3 T D P n+1 6 ṁi n l 2 l 2 4 -x 2 T = 0, 258.3 P n+1 l 2 4 -x 2 T .
Note that in [START_REF] Gravemeier | Residual-based variational multiscale methods for laminar, transitional and turbulent variable-density flow at low Mach number[END_REF] and [START_REF] Avila | A finite element dynamical nonlinear subscale approximation for the low Mach number flow equations[END_REF], different boundary conditions are considered. Consequently, the time evolution of the temperature is relatively similar but not exactly comparable with the results presented in [START_REF] Beccantini | Numerical simulations of a transient injection flow at low Mach number regime[END_REF]. If the boundary conditions given in [START_REF] Gravemeier | Residual-based variational multiscale methods for laminar, transitional and turbulent variable-density flow at low Mach number[END_REF] are enforced in the C-FV-FE scheme, we can observe analogous temperature distribution and velocity field as those presented in [START_REF] Gravemeier | Residual-based variational multiscale methods for laminar, transitional and turbulent variable-density flow at low Mach number[END_REF].

First of all, we check the grid convergence property. We use some structured meshes like the one displayed in Figure 5(b). We plot in Figure 6 the isovalues of the temperature at t f = 6. Results are obtained on three meshes corresponding respectively to grids 60 × 60, 120 × 120 and 180 × 180, and using ∆t = h max . Even if the jet obtained with the 60 × 60 grid seems to be a little delayed, we can see by comparing the solutions obtained for grids 120 × 120 and 180 × 180 that they are close to each other. Also, we plot in Figure 7 the evolution of the velocity components u x and u y and of the temperature T along the vertical lines x = 0, x = -L/4 and horizontal ones y = H/4, y = H/2 and y = 3H/4 at t = 6. The temperature distribution and the velocity field are nearly the same in both cases, so that the grid convergence can be considered as achieved. Let us note that although the shape of the cross-sections showed in Figure 7 are comparable to the profiles plotted in Figures 8 and9 in [START_REF] Beccantini | Numerical simulations of a transient injection flow at low Mach number regime[END_REF], their amplitudes are slightly greater. These differences can be explained by the ability of our scheme to preserve locally the density of the fluid. In order to confirm the observations made in section 4.1.2, we compare in Figure 8 the temperature distributions obtained with the C-FV-FE scheme and with the scheme from [START_REF] Beccantini | Numerical simulations of a transient injection flow at low Mach number regime[END_REF] (pressure-based solver, asymptotic approach 1) at time t f = 6 on a 120 × 120 grid, without and with the fixed-point process. In that second case, we require the relative error to be smaller than 10 -10 for all variables, leading to about 15 iterations. We can first observe that with some fixed-point iterations, the propagation speed of the jet is very similar for the two schemes (compare (b) and (d)). Nevertheless, without the fixed-point iterations, this propagation speed is nearly the same for the C-FV-FE scheme (compare (a) and (b)), whereas the phenomenon is delayed for the other scheme (compare (c) and (d)), like observed in Figure 10 in [START_REF] Beccantini | Numerical simulations of a transient injection flow at low Mach number regime[END_REF].

Since the law state is imposed implicitly in the C-FV-FE scheme, we check in Figure 9(a) the convergence towards zero of the discrete state equation error in the L 2 (0, t f ; L 2 (Ω)) norm. In any case (with or without fixed-point iterations), it goes towards zero, even if the convergence is faster with some fixed-point iterations. Finally, we want to investigate the thermodynamic pressure evolution. Similarly to [START_REF] Beccantini | Numerical simulations of a transient injection flow at low Mach number regime[END_REF], neglecting the diffusive term in [START_REF] Guermond | Calculation of incompressible viscous flows by an unconditionally stable projection FEM[END_REF], we get:

dP dt + γP |Ω| ∂Ω u.n 0. ( 45
)
Using the value of (ρ u) i n given in (44), and the state equation [START_REF] Dellacherie | Construction of modified Godunov-type schemes accurate at any Mach number for the compressible Euler system[END_REF] to evaluate ρ i n, we can explicitly solve (45) to obtain an approximation of the thermodynamic pressure:

P (t) P 0 + 258.3 γ 750 |Ω| t. (46)
Figure 9(b) displays the evolution of the approximate thermodynamic pressure given by ( 46), the one computed with the C-FV-FE scheme and also by the scheme from [START_REF] Beccantini | Numerical simulations of a transient injection flow at low Mach number regime[END_REF] on a 120 × 120 grid. Once again, results are close to each other. with the C-FV-FE scheme, and obtained the values of 2.26e -02 and 4.62e -04 respectively. Consequently, as the second term is two orders of magnitude below the first one, here it can be neglected. This fact can also explain that the pressure given by (46) in slightly lower that the one computed using the C-FV-FE scheme, since the second term is not only very small, but also positive.

4.3

The natural convection in a cavity

The original benchmark

This flow example has been proposed in [START_REF] Le Quéré | Modelling of natural convection flows with large temperature differences: a benchmark problem for low Mach number solvers[END_REF] as a benchmark problem for natural convection flows with large temperature gradients. It was also considered in e.g. [START_REF] Heuveline | On higher-order mixed FEM for low Mach number flows: application to a natural convection benchmark problem[END_REF][START_REF] Avila | A finite element dynamical nonlinear subscale approximation for the low Mach number flow equations[END_REF][START_REF] Gravemeier | Residual-based variational multiscale methods for laminar, transitional and turbulent variable-density flow at low Mach number[END_REF]. We consider a square cavity Ω = [0, 1] 2 containing a calorifically perfect gas, see Figure 10(a). The gas is initially at rest with uniform temperature and pressure : u 0 = 0 m.s -1 , T 0 = 600 K and P 0 = 101325 Pa.

g T 0 T h ∇T • n = 0 ∇T • n = 0 Tc (a) (b)
On all walls, the no-slip condition is imposed for the velocity, namely : similarly to [START_REF] Heuveline | On higher-order mixed FEM for low Mach number flows: application to a natural convection benchmark problem[END_REF]. We are interested in the solution when the steady-state is reached.

u| ∂Ω = 0, so that ∂Ω in = ∅. A temperature of T h = T 0 (1 + ε) (respectively T c = T 0 (1 -ε)) is imposed
A quantity of interest for this benchmark is the Nusselt number, which represents the heat transfer from the hot to the cold wall. It is defined as:

Nu(x, y) = l ref T h -T c ∇T (x, y) • n, ∀(x, y) ∈ ∂Ω.
An average Nusselt number is computed on both hot and cold walls, and is respectively denoted by Nu c and Nu h . We also define

Nu av = Nu c + Nu h 2 .
Finally, we are interested in the value of the ratio of thermodynamic pressure P over the initial pressure P 0 , and in the difference between its value and the reference one.

Like in [START_REF] Heuveline | On higher-order mixed FEM for low Mach number flows: application to a natural convection benchmark problem[END_REF], meshes are defined by the number of segment N on each of the four boundaries of the domain. Moreover, they are refined in the vicinity of the boundaries using a geometrical progression and ensuring an aspect ratio (defined as the ratio of the length of the largest segment in the mesh over the length of the smallest one) equal to 15.1, see Figure 10(b) for an example in the case N = 8. At each time step, few fixed-point iterations are performed: we require the relative error on all the variables to decrease of four order of magnitude or, except in a very short interval of time [0, t s ] (for instance t s = 10∆t), we enforce at most two fixed-point iterations. Time iterations are performed until the relative residual on all the variables is less that 10 -5 , leading to the numerical steady state. Let us note that the fixed-point iterations are here necessary to reach this steady state and to avoid some small oscillations on the velocity field.

The grid convergence study is given in Table 3, and the temperature and velocity fields are displayed for the mesh N = 256 in Figure 11. We can observe that each of the Nusselts numbers Nu c and Nu h goes towards its reference value when the mesh is being refined. The value of Nu av goes towards zero and is divided by a factor four when the value of N is twice, as well as the value of P/P 0 -(P/P 0 ) ref . We also check in Figure 12 the convergence towards zero of the discrete state equation in the L 2 (0, t f ; L 2 (Ω)) norm, and observe a convergence rate around 0.5 when the mesh is being refined.

The case without gravity

Finally, we investigate a similar case to the one developped in section 4.3.1, but without any gravity field (Ra = 0). In such a simple configuration, the exact solution is known and given by: u ex (x, y) = 0, T ex = T h + (T h -T c ) x, P ex = P 0 , ρ ex = P 0 R T ex and π ex = 0.

(47)

We aim to investigate the behaviour of the scheme from a small perturbation around this steady state. As expected, we observe that if the computation is initialized using the exact solution (47), the fluid remains at rest and the exact solution is preserved, as it was the case for the constant states test in subsection 4.1.1. Then, we start the computation by using ρ η init and T η init instead of ρ ex and T ex respectively given by: ρ η init = ρ ex + η sin(2kπx) sin(2kπy) and T η init = P 0 R ρ η init , with η = 0.01 and k = 3. Time iterations are performed until the relative residual on all the variables is less than 10 -10 , leading to the numerical steady state u h . If we compare the solution obtained with the unperturbed steady state (47), we observe that the L 2 (Ω) norms of the difference between the unperturbed state (47) and the numerical steady one for the velocity, the temperature and the dynamic pressure are less than 10 -10 whatever the mesh used, concluding that the state (47) is reached for these variables. Since the finite volume method ensures the mass preservation, the quadrature error made at initial time to evaluate the mass implies the mesh to be refined in order to ensure the convergence of ρ (s) h and P (s) towards ρ ex and P ex respectively. We display in Figure 13 the relative errors on the density, the thermodynamic pressure and the state equation. It can be observed that all these quantities tend towards zero as expected, and that the convergence is faster for P (s) than for ρ (s) .

Conclusion

In this paper, a combined Finite Volume -Finite Element method based on a time splitting has been developed. The main ingredients are:

• to solve the mass conservation equation by a FV method instead of exploiting the equation of state;

• to solve the temperature and the momentum equations by a FE method, using a projection method in order to fulfill the constraint ( 16);

• to keep a particular definition of the velocity (see [START_REF] Temam | Navier-Stokes equations. Theory and numerical analysis[END_REF]) in order to verify the divergence constraint also in the FV scheme.

We compared the numerical results of the proposed method with those obtained implementing one of the schemes proposed in [START_REF] Beccantini | Numerical simulations of a transient injection flow at low Mach number regime[END_REF]. We checked that the C-FV-FE scheme allows to preserve the constant states as well as the maximum principle on the density at the incompressible limit of the fluid, and gives the optimal rates of convergence on a smooth benchmark using an analytical solution. We then investigated a problem of injection of hot gas into a cavity filled with the same gas, obtaining results very similar to those reported in the literature, and observing that a fixed-point iterations procedure is not necessary, even if it can increase the convergence process. Finally, we simulated a fluid subject to the natural convection in a cavity, for which high temperature gradients are involved in the vicinity of the boundaries. Once again, the developped scheme provided satisfactory qualitative results as well as convergence rates.

AFigure 1 :

 1 Figure 1: Space discretization: patch Ω A . Meshes and Degrees of Freedom (DoF) for each variable.

3. 2 . 1 1 h , ρ n h and ρ n- 1 h

 2111 Solving the thermodynamic pressure Assuming that the approximated values of T n h , T n-are known, the values of P n+1 and dP dt n+1

(ρ n h = ρ n- 1 h= ρ n- 2 h

 12 H1) P n = P n-1 and dP dt n are constant in space, (H4) ∂Ω u n+1 D • n = 0. We want to prove that T n+1 h = T n h and ρ n+1 h = ρ n h . First, (30) and (31) with assumptions (H2), (H3) and (H4) lead to P n+1 = P n and dP dt n+1 = 0. Then, the weak formulation of (35) at time t n writes:

Figure 2 :

 2 Figure 2: Density, temperature, pressure and the two velocity components at t f = 1 for the mesh hmax = 0.0625.

Figure 3 : 2 Figure 4 :

 324 Figure 3: Errors in L ∞ (0, t f ; L 2 (Ω)) norm versus hmax, log/log scale.

Figure 5 :

 5 Figure 5: (a) : The cavity Ω. (b) : The grid mesh 4 × 4.

Figure 6 :

 6 Figure 6: C-FV-FE scheme, temperature at t = 6. From the left to the right: Grids 60 × 60, 120 × 120 and 180 × 180. 40 isovalues from 330 to 600, uniform distribution.

Figure 7 :

 7 Figure 7: Evolution of ux, uy and T along the lines x = 0, x = -L/4, y = H/4, y = H/2 and y = 3H/4 at t = 6, for the grids 120 × 120 and 180 × 180.

Figure 8 :Figure 9 :

 89 Figure 8: Temperature at t = 6. Comparison between C-FV-FE scheme and scheme in [2]. (a) : C-FV-FE scheme. (b) : C-FV-FE scheme with fixed-point iterations. (c) : scheme in [2] without fixed-point iterations. (d) : scheme in [2].

Remark 4 . 1 .

 41 In order to give a justification to the approximation used in the derivation of (45), we computed the average values in time of γP |Ω| ∂Ω u.n and γ -1 |Ω|Re Pr ∂Ω ∇T • n

Figure 10 :

 10 Figure 10: (a) : The differentially heated cavity. (b) : The refined mesh.

2 εg l 2 ref ν 2 = 10 6

 226 on the left (respectively right) wall with ε = 0.6. The horizontal walls are insulated. Denoting by ∂Ω N = [0, 1] × {0, 1}, we thus have :∇T • n| ∂Ω N = 0.The Prandtl and Rayleigh numbers of the flow are respectively given by Pr = 0.71 and Ra = Pr ,

Figure 11 :

 11 Figure 11: Temperature and velocity field for the mesh N = 256.

2 Figure 12 :

 212 Figure 12: Error on the state equation as a function of hmax.

  Figure 13: ||ρ(s) h -ρex|| L 2 (Ω) ||ρex|| L 2 (Ω) , |P (s) -Pex| |Pex| and ρ (s) h -P (s) R T (s) h L 2 (Ω)versus hmax.

Table 1 :

 1 Constant density case, errors in density and temperature.

  Error P Rate Error ρ Rate Error T Rate Error u Rate Error π Rate

	h max
	2.00e-2 5.12e-3	-	2.84e-3	-	3.31e-3	-	8.49e-4	-	6.44e-2	-
	1.56e-2 3.31e-3	1.76 1.76e-3 1.94	2.13e-3	1.79 4.72e-4 2.38 4.29e-2 1.65
	1.11e-2 1.74e-3	1.89 9.04e-4 1.95	1.06e-3	2.05 2.10e-4 2.38 2.38e-2 1.73
	7.81e-3 8.95e-4	1.89 4.54e-4 1.96	5.48e-4	1.87 1.03e-4 2.02 1.26e-2 1.80
	6.67e-3 6.53e-4	1.99 3.32e-4 1.97	3.87e-4	2.19 7.48e-5 2.02 9.39e-3 1.85

Table 2 :

 2 C-FV-FE scheme. Errors in L ∞ (0, t f ; L 2 (Ω)) norm and corresponding convergence rates.

  Nu h Nu av P/P 0 P/P 0 -(P/P 0 ) ref Ref. values [24] 8.85978 -8.

		Nu c
			85978	0	0.85633	0
	N = 64	8.87358	-8.98999 -0.1164 0.84769	-0.00864
	N = 96	8.83670	-8.93469 -0.0980 0.84819	-0.00814
	N = 128	8.83551	-8.86779 -0.0323 0.85375	-0.00258
	N = 192	8.84402	-8.85730 -0.0133 0.85512	-0.00121
	N = 256	8.85452	-8.86184 -0.0073 0.85564	-0.00069

Table 3 :

 3 Heated cavity, ε = 0.6 and Ra = 10 6 .

Acknowledgements

This work was supported in part by the Labex CEMPI (ANR-11-LABX-0007-01).

A Appendix

In this section, we will detail the fixed-point iterations algorithm. We perform a loop of index l. Knowing P n+1,l , T n+1,l , u n+1,l and π n+1,l , we show how to compute P n+1,l+1 , T n+1,l+1 , u n+1,l+1 and π n+1,l+1 . Note that for a generic variable a, a n+1,0 = a n . 1. The thermodynamic pressure P n+1, l+1 is computed by the global mass conservation:

3. The new density ρ n+1, l+1 is computed by solving:

with:

4. The temperature T n+1, l+1 is computed by solving:

∆T n+1, l+1 = 0, 5. The velocity u n+1, l+1 and the pressure π n+1, l+1 are computed by solving:

dP dt n+1, l+1 + γ -1 γ Re Pr P n+1, l+1 ∆T n+1, l , with ∇ • τ n+1, l+1 = ∆u n+1, l+1 + 1 3 ∇ ∇ • u n+1, l .