Franck Nicoud
email: nicoud@cerfacs.fr

Conservative High-Order Finite Difference Schemes for Low-Mach Number Flows

Keywords: 65M06, 76M20, 76N99 Low-Mach number, finite difference, staggered mesh, density gradient, conservative scheme, Poisson equation

Three finite-difference algorithms are proposed to solve a low-Mach number approximation for the Navier-Stokes equations. These algorithms exhibit fourth-order spatial and second-order temporal accuracy. They are dissipation-free, and thus well suited for DNS and LES of turbulent flows.

The key ingredient common to each of the methods presented is a Poisson equation with variable coefficient that is solved for the hydrodynamic pressure. This feature ensures that the velocity field is constrainted correctly.

It is shown that this approach is needed to avoid violation of the conservation of kinetic energy in the inviscid limit which would otherwise arise through the pressure term in the momentum equation. An existing set of finite-difference formulae for incompressible flow is generalized to handle arbitrary large density fluctuations with no violation of conservation through the non-linear convective terms. An algorithm which conserves mass, momentum and kinetic energy fully is obtained when an approximate equation of state is used instead of the exact one. Results from a model problem are used to show both spatial and temporal convergence rates and several test cases are presented to illustrate the performance of the algorithms.

Introduction

Large Eddy Simulation (LES) and/or Direct Numerical Simulation (DNS) can provide detailed information about turbulent flows that may be difficult to obtain experimentally. However, for the particular class of flow with low Mach number M and strong density variation, the classical compressible Navier-Stokes equations are not well suited for computation. The small time step limitation dictated by numerical stability requirements of explicit methods would require excessive computing times to solve practical flow problems. Indeed the sound waves move much faster than entropy or vorticity waves when M << 1. At the same time, in flows where the dominating mechanism is free, forced or mixed convection, the acoustical mode of energy carries only a small fraction of the energy present in the fluctuating part of the flow. These observations led several authors [START_REF] Majda | The derivation and numerical solution of the equations for zero mach number combustion[END_REF][START_REF] Ghoniem | Twenty-first symposium (international) on combustion[END_REF][START_REF] Mcmurthry | Direct numerical simulations of a reacting mixing layer with chemical heat release[END_REF][START_REF] Knio | Twenty-sixth symposium (international) on combustion[END_REF] to propose a set of low-Mach number equations which do not contain acoustic waves but can still describe the entropy and vorticity modes as well as compressibility due to exothermicity of chemical reactions. A fractional-step method is used most often, the pressure field being obtained by solving a Poisson equation with the time derivative of the density field as part of the source term [START_REF] Mcmurthry | Direct numerical simulations of a reacting mixing layer with chemical heat release[END_REF][START_REF] Mahalingam | Full numerical simulations of coflowing, axisymmetric jet diffusion flames[END_REF].

The projection step was found to be the most destabilizing part of the algorithm [START_REF] Cook | Direct numerical simulation of a turbulent reactive plume on a parallel computer[END_REF]. Even-ordered finite difference approximations to this derivative were found to be more stable but density ratios larger than 3 are difficult to compute. Sandoval (reported in [START_REF] Cook | Direct numerical simulation of a turbulent reactive plume on a parallel computer[END_REF]) found that by decreasing the Reynolds number, larger variations in density could be achieved. Larger density ratios seem computable by using a predictor-corrector time-stepping algorithm in which the predictor uses a second-order Adams-Bashforth time integration scheme and the corrector relies on a quasi-Crank-Nicolson integration with the inversion of a pressure Poisson equation at each step [START_REF] Rutland | Simulations of flame-vortex interactions[END_REF][START_REF] Najm | A semi-implicit numerical scheme for reacting flow[END_REF].

As far as incompressible Navier-Stokes equations are concerned, experience has shown that the kinetic energy must be conserved if a stable and dissipation-free numerical method is sought. Indeed, such property ensures that the sum of the square of the velocities cannot grow, even through nonlinear interactions between modes. As a consequence, a numerical scheme which conserves kinetic energy cannot be unstable. Moreover, it makes unnecessary the use of numerical stabilization through up-winding which is known to introduce too much artificial damping in DNS/LES computations.

Morinishi et al. [START_REF] Morinishi | Fully conservative higher order finite difference schemes for incompressible flow[END_REF] developed a set of fully conservative (mass, momentum and kinetic energy) high order schemes for incompressible flow. However, none of the numerical studies on low-Mach number flows cited above addressed this issue. In the present study it is shown that the global conservation of kinetic energy is a common feature of incompressible and low-Mach number flows in the inviscid limit. A nearly conservative fourth-order finite difference scheme is proposed in which one solves a Poisson equation with variable coefficient for pressure. Also, this algorithm makes use of a generalization of the Morinishi's finite difference formulae for variable density flow. These two ingredients lead to algorithms which are well suited for LES and/or DNS computations. In particular:

• no numerical dissipation from the spatial discretization is used to stabilize the computation; and,

• they can handle density ratios much larger than 3.

The low Mach number approximation and the numerical method are discussed in sections 2 and 3. An error analysis is conducted in section 4 and some numerical experiments are discussed in section 5 to show the potential of the algorithm.

Low Mach number approximation

In compressible flows, a natural parameter to measure the effects of compressibility is the ratio of the dynamic to the thermodynamic pressure, viz. γM 2 . To derive the low Mach number equations, one expands the dependent variables as a power series in ǫ = γM 2 , which is a small parameter (see [START_REF] Majda | The derivation and numerical solution of the equations for zero mach number combustion[END_REF][START_REF] Mcmurthry | Direct numerical simulations of a reacting mixing layer with chemical heat release[END_REF][START_REF] Cook | Direct numerical simulation of a turbulent reactive plume on a parallel computer[END_REF] for details on the derivation of the low Mach number equations).

Substituting these expansions into the compressible Navier-Stokes equations and collecting the lowest order terms in ǫ yields:

∂ρ ∂t + ∂ρu j ∂x j = 0 (1)
ρ ∂u i ∂t + ρu j ∂u i ∂x j = - ∂P ∂x i + 1 R e ∂τ ij ∂x j ( 2
)
ρC p ∂T ∂t + ρC p u j ∂T ∂x j = 1 R e P r ∂q j ∂x j + γ -1 γ dP o dt (3)
As body forces have no impact on the numerical methods presented in section 3, they have not been included in Eq. ( 2). All the variables are normalized using the reference state

ρ ref , u ref , T ref = P ref o /ρ ref , C ref p = C * p (T ref ), µ ref = µ * (T ref ) and k ref = k * (T ref ) where the superscript * represent dimen- sional quantities. Also R e = ρ ref u ref L ref /µ ref and P r = µ ref C ref p /k ref are the
Reynolds and the Prandtl number while γ is the ratio of specific heats at the reference state. u i , ρ, T , k, µ and C p stand for the non-dimensionalized velocity vector, density, temperature, thermal conductivity, dynamic viscosity and specific heat.

τ ij = µ ∂u i ∂x j + ∂u j ∂x i -2 3 δ ij ∂u k ∂x k
and q j = k ∂T ∂x j are the viscous stress tensor and the heat flux vector respectively. Moreover, P may be interpreted as the hydrodynamic pressure. In the low-Mach number approximation, the thermodynamic pressure P o only depends on time and must be computed if it is not constant. The equation of state is simply:

P o = ρT ( 4
)
Since density is uniquely determined by the temperature (and the thermodynamic pressure which is constant in space), the energy equation acts as a constraint which is enforced by the hydrodynamic pressure. Combining the equations ( 1), ( 3) and ( 4), this constraint is:

∂u i ∂x i = 1 P o (t)C p 1 R e P r ∂ ∂x j k ∂T ∂x j + γ -1 γ -C p dP o dt (5)
Integrating over the flow domain V leads to the following ODE for the thermodynamic pressure in a closed system:

dP o dt = 1 V γ-1 γ -C p dV 1 R e P r V ∂ ∂x j k ∂T ∂x j dV +P o (t) V u i ∂C p ∂x i dV (6)
Since V ∂ ∂x j k ∂T ∂x j dV = S k ∂T ∂x j dS j , this relation expresses how the rate of change of the mean pressure is affected by the heat flux through the surface S of the domain V and the gradients of heat capacity of the gas. In many practical applications such as piston engines and internal ballistics the fluid may be considered as a calorifically perfect gas so that C P = 1 and the time derivative of P o is simply:

dP o dt = γ V 1 R e P r V ∂ ∂x j k ∂T ∂x j dV = 1 V 1 R e P r S k ∂T ∂x j dS j (7)
Thus the constraint on the velocity field becomes:

∂u i ∂x i = 1 P o (t)R e P r ∂ ∂x j k ∂T ∂x j - 1 V V ∂ ∂x j k ∂T ∂x j dV (8)
Note that the numerical methods presented in section 3 remains usable even if the fluid is not considered as being a calorifically perfect gas. If the system considered is open, then the thermodynamic pressure is set by atmospheric conditions. If it is closed then the amount of mass in it, M 0 , is constant over time so that by integrating the equation of state over the whole domain one obtains the following expression for the thermodynamic pressure:

P o (t) = M 0 V 1 T dV (9)
Note that in the limit of an inviscid flow of a calorifically perfect gas the thermodymanic pressure remains constant over time (from Eq. ( 7)) and the velocity field is divergence-free (from Eq. ( 8)).

The solution (ρ, u i , T , P , P o ) is completely described by Eqs. ( 1)-( 4)

and [START_REF] Rutland | Simulations of flame-vortex interactions[END_REF]. The constraint Eq. ( 8) should also be satisfied since it is a linear combination of Eqs. ( 1), ( 3) and (4).

Numerical method

The numerical method chosen for solving the variable density momentum and temperature equations is a generalization of a fully conservative fourth order spatial scheme developed for incompressible flows on staggered grids by Morinishi et al. [START_REF] Morinishi | Fully conservative higher order finite difference schemes for incompressible flow[END_REF]. A scheme to solve the momentum equations in nonconservative form is described in the following subsection 3.1. A scheme with 'good' conservative properties is discussed in section 3.2. Both of these algorithms involve a variable coefficient Poisson equation for the pressure.

For completeness, section 3.3 presents an alternative formulation in which the pressure is obtained, as proposed in most of the previous studies, through a Poisson equation with constant coefficient and approximate source term.

Scheme in non-conservative form: AdvSC

For a uniform mesh, the advective term in the momentum equation, Eq.

(2), is discretized as:

ρu j ∂u i ∂x j ≡ 9 8 9 8 ρ (4j) u j 1x i - 1 8 ρ (4j) u j 3x i δ 1 u i δ 1 x j 1x j - 1 8 9 8 ρ (4j) u j 1x i - 1 8 ρ (4j) u j 3x i δ 3 u i δ 3 x j 3x j (10)
where the finite-difference operator with stencil n acting on φ with respect to x i is defined as

δ n φ δ n x i = φ(x i + nh i /2) -φ(x i -nh i /2) nh i (11)
and the interpolation operator with stencil n acting on φ in the x i direction is

φ nx i = φ(x i + nh i /2) + φ(x i -nh i /2) 2 . ( 12
)
ρ (4j) is a fourth order interpolation of ρ in the x j direction:

ρ (4j) = 9 8 ρ 1x j - 1 8 ρ 3x j (13)
When the density is constant, Eq. ( 10) reduces to the advective form (Adv.-S4) in Morinishi et al. [START_REF] Morinishi | Fully conservative higher order finite difference schemes for incompressible flow[END_REF]. The pressure term in Eq. ( 2) in discretized by:

∂P ∂x i ≡ (P res.) i = (∇ d P ) i = 9 8
δ 1 P δ 1 x i - 1 8
δ 3 P δ 3 x i ( 14
)
and the discrete divergence operator is defined consistently:

∂u i ∂x i ≡ ∇ d • (u i ) = 9 8
δ 1 u i δ 1 x i - 1 8
δ 3 u i δ 3 x i (15)
The viscous terms in Eq. ( 2) are written using the following generic form:

∂ ∂x j µ ∂u i ∂x j ≡ 9 8 δ 1 δ 1 x j µ (4i),(4j) 9 8 δ 1 u i δ 1 x j - 1 8
δ 3 u i δ 3 x j - 1 8
δ 3 δ 3 x j µ (4i),(4j) 9 8 δ 1 u i δ 1 x j - 1 8
δ 3 u i δ 3 x j (16)
Note that µ is successively interpolated in the i and j-direction (applying Eq. ( 13) to µ once in i, once in j) to give the fourth-order interpolation µ (4i),(4j) . The advective term for the temperature is discretized as:

ρu j ∂T ∂x j ≡ 9 8 ρ (4j) u j δ 1 T δ 1 x j 1x j - 1 8 ρ (4j) u j δ 3 T δ 3 x j 3x j (17)
A semi-implicit time marching algorithm is used in which the diffusion terms in the wall normal direction are treated implicitly with a Crank-Nicolson scheme, while a third order Runge-Kutta or second order Adams-Bashforth scheme is used for all other terms. The temperature equation is advanced first so that ρ n+1 is known via the state equation ρ = P o /T , where P o is first assessed using Eq. ( 9) written at time n + 1. Then a fractional step method is used to solve the momentum equation.

ρ (4i),n+1 u n+1 i -u n i ∆t = ρ (4i),n+1 u n+1 i -ûi ∆t +ρ (4i),n+1 ûi -u n i ∆t 10
Sub-step =

β k γ k ζ k RK-CN k=1
β k I n+1 + I n +γ k E n +ζ k E n-1 -2β k ∇ d P n -2β k ∇ d δP n+1 (18)
where I and E represent all the spatial implicit and explicit terms except for the pressure at n and the pressure update δP n+1 = P n+1 -P n .

The parameters β k , γ k and ζ k (k = 1, 3) can be chosen so that the mixed Runge-Kutta/Crank-Nicolson (RK-CN) time stepping is recovered after the third sub-step (see Spalart [START_REF] Spalart | Hybrid rkw3 + crank-nicolson scheme[END_REF]). One can also choose their values so that the mixed Adams-Bashforth/Crank-Nicolson (AB-CN) time stepping is obtained at each sub-step. The values for the coefficients β k , γ k and ζ k are given in table 1.

Equation ( 18) is then split into a decoupled set which is a second-order approximation in time to the original equation:

ρ (4i),n+1 ûi -u n i ∆t = β k I n+1 + I n + γ k E n + ζ k E n-1 -2β k ∇ d P n (19)
ρ (4i),n+1 u n+1 i -ûi ∆t = -2β k ∇ d δP n+1 (20)
Equation ( 19) is solved for ûi by using the discretizations (10), ( 14) and [START_REF] Suslov | Stability of mixed-convection flow in a tall vertical channel under no-boussinesq conditions[END_REF]. Then [START_REF] Kim | Turbulence statistics in fully developed channel flow at low reynolds number[END_REF] is divided by ρ (4i),n+1 before its discrete divergence is taken to obtain:

∇ d • 1 ρ (4i),n+1 ∇ d δP = 1 2β k ∆t ∇ d • ûi -∇ d • u n+1 i (21)
A similar Poisson equation with variable coefficient was solved in Bell &

Marcus [START_REF] Bell | A second-order projection method for variable-density flows[END_REF] to impose the divergence-free constraint for variable-density flows -see also [START_REF] Lai | A projection method for reacting flow in the zero Mach number limit[END_REF][START_REF] Lai | A projection method for combustion in the zero Mach number limit[END_REF].

Since the transport equation for T has been advanced prior to the momentum equation, the last term in the equation for the pressure variation is known from Eq. ( 8), written at time n + 1. The variable coefficient Poisson equation Eq. ( 21) for the pressure is solved using the (pre-conditioned) conjugate gradient algorithm. In the case where homogeneous directions exist, it is worth pre-conditioning Eq. ( 21) by the elliptic operator:

∇ d • 1 ρ (4i),n+1 ∇ d (22)
where <> denote a spatial averaging in the homogeneous directions. In this case one can make use of FFT-based Fast Poisson Solver at each iteration of the conjugate gradient algorithm. In the more general case it is worth solving for the modified unknown δP/ √ ρ [START_REF] Concus | Use of fast direct methods for the efficient numerical solution of nonseparable elliptic equations[END_REF]. The advantage of solving Eq. ( 21) to update the pressure is that the divergence-free constraint is recovered in the inviscid limit, as it has to be from Eq. ( 8). This is not the case when a backward approximation of ∂ρ ∂t is used to compute the source term of a linear Poisson equation for δP as proposed earlier [START_REF] Mcmurthry | Direct numerical simulations of a reacting mixing layer with chemical heat release[END_REF], [START_REF] Cook | Direct numerical simulation of a turbulent reactive plume on a parallel computer[END_REF] (see also section 3.3). The other advantage is that the pressure terms remain energy conserving in the high-R e number limit, as discussed in the following subsection.

Fully conservative scheme: DivSC, DivSC approx

Although the previous scheme was found to be accurate, it only conserves momentum and kinetic energy to its own order of accuracy. Experience has shown that the latter quantity must be conserved exactly if a stable and dissipation-free numerical method is sought. Morinishi et al. [START_REF] Morinishi | Fully conservative higher order finite difference schemes for incompressible flow[END_REF] developed a set of fully conservative (mass, momentum and kinetic energy) high order schemes for incompressible flow. In the general case of the Navier-Stokes equations without body force, the transport equation for the kinetic energy ρk per unit volume reads:

∂ρk ∂t + ∂ρu j k ∂x j = P S jj - ∂P u j ∂x j + ∂τ ij u i ∂x j -τ ij S ij (23)
Let us consider a periodic (or infinite) domain so that, after Eq. ( 23) is integrated over the domain, the flux terms ∂ρu j k ∂x j and ∂P u j ∂x j make no contribution. Due to the dissipation term τ ij S ij , the question of conservation of the kinetic energy is only relevant in the inviscid limit where τ ij = 0. We know from Eq. (8) that in this limit, the velocity field is divergence-free, that is S jj = 0. Thus global conservation of kinetic energy is a common feature of incompressible and low Mach number flows. The purpose of this section is to investigate how this property can be extended in discrete space.

Let us define the following discrete approximations of the possible forms for the non-linear term in the momentum equation:

(Adv.) i = 9 8 9 8 ρ (4j) u j 1x i - 1 8 ρ (4j) u j 3x i δ 1 u i δ 1 x j 1x j - 1 8 9 8 ρ (4j) u j 1x i - 1 8 ρ (4j) u j 3x i δ 3 u i δ 3 x j 3x j (24) (Div.) i = 9 8 δ 1 δ 1 x j 9 8 ρ (4j) u j 1x i - 1 8 ρ (4j) u j 3x i u i 1x j - 1 8
δ 3 δ 3 x j 9 8 ρ (4j) u j 1x i - 1 8 ρ (4j) u j 3x i u i 3x j (25) (Skew.) i = 1 2 ((Adv.) i + (Div.) i ) (26)
The forms (Adv.) i , (Div.) i and (Skew.) i are the discrete equivalent to the advective ρu j ∂u i ∂x j , conservative

∂ρu i u j ∂x j and skew-symmetric 1 2 ρu j ∂u i ∂x j + ∂ρu i u j ∂x j
form of the convective term. Note that the discrete operator in Eq. ( 24) is the same as that in Eq. ( 10). The following relations holds between these three discrete forms:

(Div.) i = (Adv.) i + u i 9 8 (Cont.) 1x i - 1 8 (Cont.
)

3x i (27) (Skew.) i = (Adv.) i + 1 2 u i 9 8 (Cont.) 1x i - 1 8 (Cont.
)

3x i (28) (Skew.) i = (Div.) i - 1 2 u i 9 8 (Cont.) 1x i - 1 8 (Cont.) 3x i (29)
where

(Cont.) = 9 8 δ 1 ρ (4j) u j δ 1 x j - 1 8
δ 3 ρ (4j) u j δ 3 x j (30)
is the discrete form of the divergence of the momentum vector ∂ρu j ∂x j . A key assumption in the semi-discrete analysis proposed in Morinishi et al. [START_REF] Morinishi | Fully conservative higher order finite difference schemes for incompressible flow[END_REF] for incompressible flow is that the operator (Cont.) is identically zero so that the three forms (Div.) i , (Adv.) i and (Skew.) i are equivalent. Since (Div.) i is conservative a priori for the momentum equation and (Skew.) i is conservative a priori in the kinetic energy equation, a fully conservative scheme is obtained as soon as the velocity constraint ∂u j ∂x j = 0 is imposed properly through the pressure correction step. In the present case where the density is not constant, the velocity constraint ∂u j ∂x j = 0 (in the high-R e limit) does not imply that ∂ρu j ∂x j is zero. Thus the discrete operators (Div.) i , (Adv.) i and (Skew.) i are not equivalent in the low Mach number case, meaning that a fully discrete analysis (including the time discretization) must be conducted to achieve conservation of kinetic energy.

A conservative scheme for the momentum can be derived by considering the momentum equation in its divergence form. The first guess for the velocity is obtained by solving:

ρ(4i) ûi -ρ (4i),n u n i ∆t = -γ k (Div.) n i -ζ k (Div.) n-1 i -2β k (P res.) n i ( 31
)
where ρ can be either ρ n or ρ n+1 or any intermediate value. Then the projection step is written as:

u n+1 i = ρ(4i) ρ (4i),n+1 ûi -2β k 1 ρ (4i),n+1 ∇ d δP ∆t ( 32
)
where the Poisson equation for δP must be:

∇ d • 1 ρ (4i),n+1 ∇ d δP = 1 2β k ∆t ∇ d • ρ(4i) ρ (4i),n+1 ûi -∇ d • u n+1 i (33)
Obviously, Eqs. (31), (32) and (33) constitute a scheme which is momentum conserving. To investigate whether it also conserves kinetic energy, let us multiply Eq. (31) by ûi + u n i and integrate over the whole domain. The overall contribution of the pressure term involving u n i in the kinetic energy equation behaves like:

V u n i (P res.) n i dV = V 9 8 u n i δ 1 P δ 1 x i 1x i - 1 8 u n i δ 3 P δ 3 x i 3x i dV (34) where 9 8 u n i δ 1 P δ 1 x i 1x i - 1 8 u n i δ 3 P δ 3 x i 3x i = 9 8 δ 1 u n i P 1x i δ 1 x i - 1 8 δ 3 u n i P 3x i δ 3 x i -P 9 8 δ 1 u n i δ 1 x i - 1 8
δ 3 u n i δ 3 x i (35)
The first two terms of the RHS of Eq. ( 35) do not contribute because they are in divergence form. The last one is identically zero because the variable coefficient Poisson equation, Eq. ( 33), is solved with

∇ d • u n+1 i = 9 8 δ 1 u n+1 i δ 1 x i -1 8 δ 3 u n+1 i δ 3 x i
= 0 imposed in the source term (i.e. the projection step Eq. ( 32) imposes the divergence-free constraint exactly). The other pressure term is ûi (P res.) n i and its overall contribution is of order ∆t because ûi = u n i + O(∆t) and the integral of u n i (P res.) n i is zero. Using Eq. ( 26), the overall contribution of the convective terms in the RHS of Eq. ( 31) may be written as

V (û i + u n i ) -γ k (Skew.) n i -ζ k (Skew.) n-1 i dV -γ k V ûi + u n i 2 u n i 9 8 (Cont.) n 1x i - 1 8 (Cont.) n 3x i dV -ζ k V ûi + u n i 2 u n-1 i 9 8 (Cont.) n-1 1x i - 1 8 (Cont.) n-1 3x i dV (36)
The first integral in Eq. (36) contributes to the order ∆t because (Skew.) i is kinetic energy conserving in nature and because ûi , u n i and u n-1 i are equal to the order ∆t. Thus an approximation to the order ∆t of the overall contribution of the full RHS of Eq. ( 31) is:

-γ k V ûi + u n i 2 u n i 9 8 (Cont.) n 1x i - 1 8 (Cont.) n 3x i dV -ζ k V ûi + u n i 2 u n-1 i 9 8 (Cont.) n-1 1x i - 1 8 (Cont.) n-1 3x i dV (37)
On the other hand, the contribution of the LHS of Eq. (31) may be written:

V ρ(4i) (û i ) 2 -ρ (4i),n (u n i ) 2 ∆t dV + V u n i ûi ρ(4i) -ρ (4i),n ∆t dV (38)
Comparing Eqs. (37) and (38) it appears that the discrete rate of change of the kinetic energy (the first integral in Eq. ( 38)) is at most of order ∆t if one defines the intermediate density as:

ρ -ρ n ∆t = -γ k (Cont.) n -ζ k (Cont.) n-1 (39)
In the context of a second order scheme, the same definition of ρ was adopted (C. Pierce, private communication) to achieve approximate conservation of kinetic energy.

Multiplying the projection step Eq. (32) by ûi +u n+1 i and integrating over the whole domain, the following expression can be derived:

V ρ (4i),n+1 u n+1 i 2 -ρ(4i) (û i ) 2 ∆t dV = V ûi u n+1 i ρ(4i) -ρ (4i),n+1 ∆t dV +O(∆t) (40)
This shows that the global rate of change of the kinetic energy is of order

∆t only if ρ(4i) -ρ (4i),n+1 = O(∆t n ), n ≥ 2.
Unfortunately, n is only 1 in the most general case.

A conservative scheme (DivSC approx ) is obtained if one accepts that the state equation, Eq. ( 4), is verified to the order ∆t only, viz:

ρ n+1 = ρ = P o T n+1 + O(∆t) (41)
In this case, the error in the global kinetic energy conservation is at most of order ∆t for each sub-step. If one accounts for the cancellation of error in the full third-order Runge-Kutta procedure, it can be shown that the error is actually of order ∆t 3 . The same order was achieved for incompressible flows in Morinishi et al. [START_REF] Morinishi | Fully conservative higher order finite difference schemes for incompressible flow[END_REF]. The Adams-Bashforth time integration is not considered here since it is unstable for inviscid flows (see section 4).

Scheme with approximate Poisson equation: DivSC ρ

The third algorithm is similar to the one in section 3.2 except for the Poisson equation for the pressure. It is briefly described in this section for completeness. One keeps the same scheme, Eq. ( 31), for the first guess of the velocity.

The projection step is now written as:

ρ (4i),n+1 u n+1 i = ρ(4i) ûi -2β k ∇ d δP ∆t (42)
where the Poisson equation for δP must be:

∇ d • ∇ d δP = 1 2β k ∆t ∇ d • ρ(4i) ûi -∇ d • ρ (4i),n+1 u n+1 i (43)
The divergence of momentum at time level n + 1 is unknown but it can be assessed by using the continuity equation:

∇ d • ρ (4i),n+1 u n+1 i = - δρ δt n+1 (44)
Since the temperature field is advanced first, the equation of state, Eq. ( 4), can be used to compute the density at level n + 1. Then a backward discretization for the time derivative of density can be used to assess the source term in Eq. ( 43). Since δP is of order ∆t, the source term in the Poisson equation is of order ∆t too. Since in this source term the time derivative of density is divided by the time step (see Eq. ( 43)), the backward discretization must be at least second-order accurate for consistency. A possible choice is:

δρ δt n+1 = (∆t n + ∆t n-1 ) 2 -∆t 2 n ρ n+1 -(∆t n + ∆t n-1 ) 2 ρ n + ∆t 2 n ρ n-1 ∆t n ∆t n-1 (∆t n + ∆t n-1 ) (45)
Obviously, Eqs. (31), ( 42) and (43) constitute a scheme which is momentum conserving. However it is not kinetic energy conserving because the divergence-free constraint is not recovered in the inviscid limit. Thus the overall contribution of the pressure term does not vanish as in section 3.2.

The test cases in section 5 show that this algorithm is less stable and accurate than DivSC and AdvSC. However, the Poisson equation, Eq. ( 43), is with constant coefficient and can be solved very efficiently. Subsequently, it seems that this approach was used in all the previous studies dealing with low-Mach number flows [1, 2,

. The approach involving a Poisson equation with variable coefficient as in sections 3.1 and 3.2 was preferred in [START_REF] Bell | A second-order projection method for variable-density flows[END_REF][START_REF] Lai | A projection method for reacting flow in the zero Mach number limit[END_REF][START_REF] Lai | A projection method for combustion in the zero Mach number limit[END_REF].

The conservative properties of the schemes described in section From the channel flow computations performed in the course of this study, (see section 5.5), the CPU time required per iteration for DivSC is between 1.2 and 2 times higher than what needed for DivSC ρ , the higher the temperature ratio, the higher the cost. The cost difference arises mainly through the resolution of the variable coefficient Poisson equation ( 33). This difference is expected to decrease in cases with complex geometry which require the use of an iterative method for Eq. ( 43) of algorithm DivSC ρ . Recall that a FFT-based Fast Poisson Solver has been used in the present study.

The cost for algorithm AdvSC is slightly smaller than for DivSC since the continuity equation is not advanced in time.

Error Analysis

In order to investigate the stability limit of the previous algorithms in the 1D linear case, a Fourier analysis was performed which is presented in section 4.1. Convergence tests are then described in sections 4.2 and 4.3 to show both time and space accuracy of the methods developed in this paper.

Stability

A Von Newmann analysis of the linear 1D convection-diffusion equation

∂u ∂t + c ∂u ∂x = ν ∂ 2 u ∂x 2 (46)
was performed to assess the stability of the first step of the algorithms described in the previous section. Note that in this simple case the forms (Adv.) i and (Div.) i in Eqs. ( 24) and (25) reduce to (convection velocity c omitted):

∂u ∂x ≡ 27 (u i+1 -u i-1 ) -(u i+3 -u i-3 ) 48 h (47)
Also, in the case of constant viscosity, the spatial differencing Eq. 16 for the diffusion term is proportional to:

∂ 2 u ∂x 2 ≡ -1460 u i + 783 (u i+1 + u i-1 ) -54 (u i+2 + u i-2 ) + (u i+3 + u i-3 ) 576 h 2 (48)
Denoting the semi-discrete operator A 4 such that ∂u/∂t = A 4 u n i , it is straightforward to demonstrate using Eqs. ( 47 Figure 1a shows the same quantities for the classical second-order semi-discretization where ∆t Â2 = -j CF L sin (ωh)+2F o(cos (ωh)-1). As expected, the fourth-order scheme appears to be less stable than the second-order one. Since both semidiscretizations Â2 and Â4 are anti-symmetric in the convection part, the Adams-Bashforth time differencing is always unstable for F o = 0. For finite viscosity, both spatial schemes are stable, although the domain of stability is smaller for the present formulation. For pure diffusion, viz. CF L = 0, the critical Fourier number for Â2 is 0.25 whereas for Â4 it is 0.197.

The same analysis was performed for the third-order Runge-Kunta time differencing (not shown). In the absence of viscosity (F o = 0), the critical CF L number is 1.54, compared to the equivalent value for the classical second-order scheme which is 1.8.

Time accuracy

Time accuracy of the numerical methods described in section 3 was assessed through a numerical experiment in the 2D periodic domain [0, 2π] × [0, 2π]. The initial condition was set to u(x, y, t = 0) = ∂ψ/∂y + sin (x) × cos (y)/(R e P r ), v(x, y, t = 0) = -∂ψ/∂x+cos (x)×sin (y)/(R e P r ), P (x, y, t = 0) = 0 and T (x, y, t = 0) = 1cos (x) × cos (y) so that the constraint on the velocity field, Eq. 8, is satisfied initially (assuming constant viscosity).

The stream function is equal to ψ = cos (x) × cos (y). Several simulations were performed with R e = 2000 and P r = 1 and using a uniform 16x16 grid. The number of time-steps (Adams-Bashforth time integration) used to cover a given time interval was increased from 4 to 32, corresponding to Courant numbers in the range [0.011,0.088]. One more computation was performed for which 256 iterations were necessary in order to generate a highly accurate reference solution. The error in the solutions obtained at larger ∆t was computed as the L 2 -norm of the difference in the solution for a given time increment when compared with the reference solution. Regarding the u-velocity, second-order accuracy is demonstrated in figure 2 for all the considered algorithms. The same convergence rate was obtained for the v-velocity and for the temperature. An interesting point to note is that the convergence rate for the pressure is only unity. This is not related to the density variations considered here since the same phenomenon is true for the incompressible Morinishi's algorithm [START_REF] Morinishi | Fully conservative higher order finite difference schemes for incompressible flow[END_REF]. As discussed by Perot [START_REF] Perot | An analysis of the fractional step method[END_REF],

the order of accuracy of the pressure update does not affect the order of accuracy of the velocity field.

Spatial accuracy

Other numerical experiments were conducted to study the spatial accuracy of the algorithms. Only the results for AdvSC and DicSC are discussed in this section since the spatial differencing for all the previous algorithms reduces to one or the other of these two approaches. The initial condition was the same as that used in section 4.2 and four different uniform grids were considered, namely ones containing 16 2 , 48 2 , 144 2 and 432 2 cells. The solutions obtained after a given integration time are denoted S 16 , S 48 , S 144 and S 432 respectively (S may denote either one of the velocity components, the pressure or the temperature). This grid sampling was chosen so that all the velocity and pressure nodes of the 16 2 staggered grid are common to all the grids. The L 2 -norm of the difference between two solutions can be computed on the coarser mesh without further interpolation. It is believed that in doing so one obtains a more reliable estimation of the spatial convergence rate. For all the computations the Courant number was kept to a small value (CF L ≈ 5 × 10 -3 ) so that the errors associated with the time differencing are small with respect to the spatial differentiation. The number of iterations was 16 for the coarsest grid, 432 for the finest one. The convergence rate σ was computed in two different ways. First the solution obtained with the finest grid (S 432 ) was interpreted as the 'exact' solution.

Then assuming a spatial error of the form C.h σ , the convergence rate was estimated by:

σ ≈ ln [L 2 (S n -S 432 )] -ln [L 2 (S m -S 432 )] ln 3 (50)
where (n,m) is either [START_REF] Suslov | Stability of mixed-convection flow in a tall vertical channel under no-boussinesq conditions[END_REF]48) or (48,144). The other way of computing σ does not require knowledge of the exact solution [START_REF] Bell | A second-order projection method for variable-density flows[END_REF] and assumes that the L 2 -norm of the difference between two successive grid levels is proportional to the grid spacing in the coarser mesh. It is then straightforward to show that:

σ ≈ ln [L 2 (S n -S m )] -ln [L 2 (S m -S k )] ln 3 (51)
where (n,m,k) is either [START_REF] Suslov | Stability of mixed-convection flow in a tall vertical channel under no-boussinesq conditions[END_REF]48,144) or (48,144,432). Both estimate Eqs.

(50) and (51) were applied to the velocity components, the pressure and the temperature fields. The results are shown in table 3. The fourth-order accuracy of the methods presented in section 3 is demonstrated.

Numerical tests

A few test cases were designed to illustrate the performance of the numerical methods described in the previous section, including the computation of a The time integration is AB-CN in sections 5.2 and 5.5. The results obtained are discussed in more detail in the following subsections.

1D Euler convection

If the Peclet number is infinite, the velocity field must be divergence-free, that is u must be constant in 1D. Also, the pressure should remain constant.

To test the ability of the different formulations to reproduce this feature of Eqs. ( 1)-( 3), consider the domain 0 ≤ x ≤ 1, periodic in x. The initial condition is u = u 0 = 1, P = 0 and T = 1 + A exp -x-x 0 a 2 with A = 1,

x 0 = 0.5 and a = 0.05. When the grid contains N x = 24 points, only 6 points are used to describe a Gaussian perturbation. Figure 4 shows P rms /ρ 0 u 2 0 as a function of the grid spacing (this quantity is also the error in P rms since the exact solution is P rms = 0). Three grid levels were considered: 24, 48 and 96 points in x. The rms of pressure is assessed for the time t = 20a/u 0 .

The CFL number is of order 0.5 in all cases. Since AdvSC renders the exact solution for this particular test case, the corresponding results are not included in the figure. The divergence is zero for DivSC because it is explicitly enforced through the Poisson equation, Eq. ( 33). On the other hand, the divergence-free constraint is not enforced in the algorithm DivSC ρ and the velocity does not remain constant. Since the error in the backward difference equation, Eq. ( 45), is proportional to a high-order time derivative of density, one expects the error in the source term of the Poisson equation, Eq. ( 43), to increase with the amplitude of the density perturbation. In other words, the error in the divergence of velocity should be proportional to the amplitude A of the initial temperature perturbation, at least for small values of A. This is illustrated in figure 5.

Small 1D perturbations

In the case where the Reynolds number is finite but where the perturbation in temperature is small (A << 1), the analytical resolution of the equations ( 1)-( 3) can be conducted and the structure of the perturbation which propagates is given by:

ρ ′ = - ρ 0 T 0 T ′ (52) u ′ = 1 R e P r ∂T ′ ∂x (53) p ′ = 4 3R 2 e P 2 r ∂ 2 T ′ ∂x 2 P r - 3 4 (54)
An interesting feature is that the pressure fluctuation should vanish in the limit P r = 3 4 . Figure 6 shows the error in Eq. ( 54) in the case N x = 24, a = 0.05, A = 0.01 and R e = 50. The initial condition is uniform for u and P and the physical time simulated is large enough (t ≈ 160a/u 0 ) so that the values reported in the figure are asymptotic values. Both DivSC and AdvSC are in complete agreement with the theory. In order to derive the analytical relationship Eq. ( 54) between p ′ and T ′ , one makes use of the constraint on the divergence of the velocity field which reduces to Eq. ( 53) in the present case. Since this constraint is not properly imposed by the algorithm DivSC ρ , there is no reason to expect that this method gives results in agreement with Eq. ( 54). Instead, the pressure fluctuation should only depend on the initial temperature fluctuation (shape and amplitude).

Accordingly, figure 7 shows that the rms of pressure behaves like P 2 r when scaled by 4

3R 2 e P 2 r ∂ 2 T ′ ∂x 2 rms
as in figure 6.

2D Random perturbations

To validate the results of section 3.2 with numerical tests, inviscid flow simulations were performed on a 2D periodic domain. The analytical solution dictates that the total momentum in each direction ρu i and total kinetic energy K = 1 2 ρu 2 i should be conserved in time. The domain is 0 ≤ x ≤ L, 0 ≤ y ≤ L and a 24x24 mesh is used. Solenoidal velocity fields are used as the initial condition together with random temperature fluctuations.

The initial mean kinetic energy is of order 1.5 while T rms ≈ 0.15 < T > at t = 0. Figure 8 shows the relative error for the total kinetic energy K 0 -K K 0 after an integration time of t = 0.125L/ K 0 . As expected from section 3.2, the error for the scheme DivSC does not behaves like ∆t 3 because of the violation in the conservation of kinetic energy related to the nonlinear equation of state Eq. ( 4). On the other hand, it appears that the same scheme with the approximate equation of state (Eq. ( 41)) conserves the global kinetic energy to the order ∆t 3 as predicted in section 3.2. The scheme DivSC ρ violates the conservation of kinetic energy even if the approximate equation of state is used. This is because the divergence-free constraint is not recovered in the inviscid limit when the Poisson equation, Eq. ( 43), is used. Thus the pressure term is not kinetic energy conserving for DivSC ρ (see section 3.2, Eq. ( 35)) which was found unstable for the present test case. The convective terms in the algorithm AdvSC violate the conservation of kinetic energy.

Linear stability in a channel

To check the accuracy of the algorithm in the case where the physical properties vary in space and time through the temperature, the evolution of low amplitude eigenmodes in laminar channel flow is simulated. The linear stability problem in a channel flow between two isothermal walls with temperature T 1 = 1 -δT 2 and T 2 = 1 + δT 2 was studied by Suslov & Paolucci [START_REF] Suslov | Stability of mixed-convection flow in a tall vertical channel under no-boussinesq conditions[END_REF] under the low Mach number assumption. They found that the critical Reynolds number increases with the parameter δT 2 . It is of order 40, 000 for δT 2 = 0.4, compared to 5, 772 in the isothermal case (δT = 0). In their analysis the dimensionless thermal conductivity and dynamic viscosity are given by Sutherland's law:

k(T ) = T 3/2 1 + S k T + S k µ(T ) = T 3/2 1 + S µ T + S µ (55)
where S k = 0.648 and S µ = 0.368 for air at T ref = 300K and normal pressure. The molecular Prandtl number is 0.76. In the computation, the length of the periodic domain in x is L = 2π/α, where α is the wave number of the mode of interest. The initial condition consists of a small amplitude (0.01%) random noise on u, v superimposed to the laminar solution of the problem (Suslov & Paolucci [START_REF] Suslov | Stability of mixed-convection flow in a tall vertical channel under no-boussinesq conditions[END_REF]). A stretched grid is used in the normal direction in order to capture the eigenvector accurately near the walls. The wall normal velocity points are distributed according to a hyperbolic tangent function (j = 0, 1, 2, ..., N ):

y v (j) = y j+ 1 2 = tanh γ 2j N -1 tanh(γ) (56)
A typical result is shown in figure 9. In this case the resolution is 24x100 with γ = 2 for the stretching parameter. The CFL number is fixed at 1.

The length of the domain is L = 2.4πh (α = 5 6 1 h ) and the Reynolds number is 45, 000, based on the maximum velocity and the channel half-height h.

schemes AdvSC, DivSC, DivSC approx and DivSC ρ give very similar results (AdvSC shown). The success of DivSC ρ in this test case is due to the fact that only small temperature (density) perturbations are considered for comparison with the linear theory. Thus the error in the time derivative of density is negligible compared to the other space discretization errors (see section 5.1) and the algorithm provides a reasonable answer.

The algorithms AdvSC and DivSC have been found to give consistent results over a wide range of flow and numerical conditions (temperature ratio from 1.01 to 2.33, grid size from 16x100 to 56x350). In order to establish the qualitative nature of the algorithms, typical 2D plots of velocity, pressure and temperature fluctuations are shown in figure 10 in which the complete domain is shown. Figure 11 shows the iso-lines for u and T for the region near the cold (bottom) wall of the channel. The conditions are the same as for figure 9. The modal structure shown in figures 10 and 11 is in very good agreement with previous theoretical findings [START_REF] Vasilyev | Stability of unstable stratified shear flow in a channel under non-boussinesq conditions[END_REF].

Turbulent flow in a channel

Several computations of a turbulent flow in a channel with isothermal walls have been performed. Results from DNS's are presented first. Results from coarse grid computations are then discussed in terms of stability.

Direct Numerical Simulations

Three Direct Numerical Simulations of a channel flow between two isothermal walls with temperatures T 1 and T 2 were performed. A few results are shown in this section to illustrate the performance of the method proposed in sections 3.1 and 3.2. A more detailed analysis can be found in [START_REF] Nicoud | Numerical study of a channel flow with variable properties[END_REF][START_REF] Nicoud | Dns of a channel flow with variable properties[END_REF].

Case T 2 /T 1 R τ 1 -R τ 2 R b 10 3 C f 1 10 3 C f 2 B q1 B q2 A 1.
The computations presented here were performed with AdvSC. Details of the test cases adopted are given in table 4.

The indices '1' and '2' denote the cold and hot wall respectively. The bulk Reynolds number R b is based on the bulk velocity and the values of density and dynamic viscosity corresponding to the bulk temperature.

The friction coefficients C f are based on the mean density in the channel and the maximum velocity while the heat flux parameter is defined as

B q = q w /ρ w C p u τ T w
, where u τ is the friction velocity √ τ w ρ w , q w the heat flux, C p the constant-pressure specific heat and T w the temperature at the wall considered. The friction Reynolds number R τ is based on the friction velocity. In Case A the temperature is almost uniform and the results may be compared to a previous incompressible DNS performed by Kim et al. [START_REF] Kim | Turbulence statistics in fully developed channel flow at low reynolds number[END_REF][START_REF] Kim | Transport of passive scalars in a turbulent channel flow[END_REF]. In Cases B-C one expects the temperature (density) variations to be large enough to modify the momentum balance through both viscous and inviscid effects. In each case the domain size is (4πh,2h,4πh/3) and the grid contains 120x100x120 cells. The mean flow is along the xdirection. The statistics were obtained over a time period of order 10h/u τ ,

where

u τ = u τ 1 +u τ 2 2
is the mean friction velocity. The wall normal velocity points are distributed according to a hyperbolic tangent function, Eq. (56).

The grid spacing is equivalent for all cases with ∆x + ≈ 20, ∆y + wall ≈ 0.3, ∆y + max ≈ 9 and ∆z + ≈ 6, where the superscript ' + ' denotes the usual wall scaling U + = U/u τ and y + = ρ w u τ y/µ w . The molecular Prandtl number is 0.76 and the dimensionless thermal conductivity and dynamic viscosity are chosen to be proportional to 1/ √ T so that the Reynolds number near the hot wall is not (too) small in comparison with its value near the cold wall [START_REF] Nicoud | Numerical study of a channel flow with variable properties[END_REF] -see table 4. It was shown in [START_REF] Nicoud | Dns of a channel flow with variable properties[END_REF] that the ratio G r /R 2 τ (G r is the Grashof number) is small compared to unity as long as the characteristic length scale of the channel is of order 1 cm or less. This would be a reasonable range for a true experiment with temperature differences typical of a laboratory combustion chamber. Consequently, the buoyancy effects are neglected in the present test cases.

Many experimental data support the validity of the Van Driest [START_REF] Van Driest | Turbulent boundary layer in compressible fluids[END_REF] transformation for wall-bounded turbulent flows with variable density. This transformation reads:

U + V D = u + o ρ ρ w 1/2 du + = 1 κ ln y + + C (57)
Figure 12 shows that the present calculations also give support to the transformation (57). When the classical incompressible coordinates are used instead, a logarithmic region can hardly be observed. Note also that Case A is in very good agreement with the law-of-the wall u + = 2.5 ln y + + 5.5, as proposed in [START_REF] Kim | Turbulence statistics in fully developed channel flow at low reynolds number[END_REF]. Regarding the velocity fluctuations and the Reynolds stress -u ′ v ′ , this computation is also in very good agreement with the incompressible DNS data from Kim et al. [START_REF] Kim | Turbulence statistics in fully developed channel flow at low reynolds number[END_REF] -see figure 13. Figure 14 shows a very good qualitative agreement for the temperature fluctuations in Case A when compared to a passive scalar calculation [START_REF] Kim | Transport of passive scalars in a turbulent channel flow[END_REF]. The slight quantitative difference is most likely due to the difference in Prandtl numbers (P r = 0.76 for Case A, P r = 0.71 in [START_REF] Kim | Transport of passive scalars in a turbulent channel flow[END_REF]). In each of the cases A-C, the streamwise velocity-temperature correlation coefficient is very high near the walls (of order 0.95, as in [START_REF] Kim | Transport of passive scalars in a turbulent channel flow[END_REF]). Subsequently the temperature field reveals the elongated streaky structure of turbulence near the wall. This is illustrated for Case A in figure 15 which shows iso-lines of u-velocity and temperature fluctuations in a plane located approximately 5 wall units above the cold wall. As already observed in previous studies, the temperature field is a good marker for the bursting events, as shown in figure 16 for Case E. Of interest also in this figure is the position of the iso-line

T = (T 1 + T 2 )/2.
It is clearly skewed towards the hot wall. Indeed, the Reynolds number in the hot region being smaller than in the cold region -see table 4-, the linear sub-layer close to the top wall is thicker. As a consequence, the mean temperature decreases faster near the top wall (y = +1) than it increases near the bottom one (y = -1). The iso-lines of the w-velocity in figure 16 reveals the streamwise vortices responsible for the streaky structure shown in figure 15.

Coarse grid computations

In algorithms of the type DivSC ρ where a constant coefficient Poisson equation (43) is solved for the pressure field, the projection step was found to be the most destabilizing part of the method and density ratios larger than 3 are difficult to compute [START_REF] Cook | Direct numerical simulation of a turbulent reactive plume on a parallel computer[END_REF]. In the present algorithms AdvSC, DivSC or DivSC approx , the exact constraint on the velocity field is enforced. Further computations of the channel flow were performed to see whether the methods proposed are more stable. The numerical and physical parameters are the same as before, except for the grid which is now 32x80x32. This mesh is too coarse to reproduce all the scales that are dynamically important. However, its resolution makes it computationally inexpensive to use and it was therefore considered suitable for investigating the stability characteristics of the different algorithms. The temperature ratios considered were

  ) and (48) that its Fourier transform Â4 (ω) may be expressed as: ∆t Â4 = -j CF L 24 (27 sin (ωh)sin (3ωh)) + F o 288 (-730 + 783 cos (ωh) -54 cos (2ωh) + cos (3ωh)) (49) where CF L = c ∆t/h and F o = ν ∆t/h 2 denote the Courant and the Fourier numbers respectively while j 2 = -1. For the full discretization to be stable, it is necessary that ∆t Â4 (ω) be contained in the stability region of the time differencing for all values of ω. The locus of ∆t Â4 (ω) is shown for different values of CF L and F o together with the stability region of the Adams-Bashforth algorithm in figure 1b.

  pulse of density in a viscous flow, 2D random fluctuations of velocity and temperature, 2D small amplitude fluctuations in a channel, and 3D large amplitude fluctuations in a channel. A brief description of the different test cases is provided in figure3which also summarize the assumptions and the objectives for each test case. The time integration must be RK-CN for the inviscid test cases while either RK-CN or AB-CN can be used for the viscous ones. RK-CN has been used for the test cases in sections 5.1, 5.3 and 5.4.

1 h- 1

 11 Unit of time is h/u τ . . . . . . . . . . . . . . . . . . . . . . . . 51 10 Fluctuations associated with the most unstable linear mode for α = 5 6 , R e = 45, 000 and δT 2 = 0.4. From top to bottom: u-velocity, v-velocity, pressure and temperature fluctuations. Dashed iso-lines are negative values. The mean flow is from left to right and the bottom wall is the cold one. The mesh size is 24x100. . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 11 Fluctuations of u-velocity (Top) and Temperature (Bottom). See the previous figure for the conditions. Zoom of region near the cold wall. . . . . . . . . . . . . . . . . . . . . . . . . 53 12 Mean velocity profiles. Top figure without the Van Driest transformation Eq. (57), bottom figure with it. : u + = y + and u + = 2.5 ln y + + 5.5. . . . . . . . . . . . . . . . . . . . . . . . 54 13 Root-mean-square velocities and Reynolds stress from the DNS (Case A) in wall units. Symbols from Kim et al. [20]. • : u + rms ; △ : v + rms ; : w + rms ; ⋄ : u ′ v ′ + . . . . . . . . . . . . . . . 55 14 Root-mean-square of the temperature fluctuations from Cases A, B and C, non-dimensionalized by the temperature difference T 2 -T 1 . : Case A; : Case B; : Case C; • : Passive scalar case with P r = 0.71, Kim and Moin [21]. . . 56 15 DNS of a channel flow with T 2 /T 1 = 1.01 (Case A). The position of x -z plane is y/h ≈ -0.975, i.e. 5 wall units above the cold wall. The flow is from left to right. Top: Contours of u + -velocity fluctuations: the range is [-2,+6], the increment is 1, dashed contours represent negative values. Bottom: Contours of temperature fluctuations (nondimensionalized by T 2 -T 1 ): the range is [-0.04,0.04], the increment is 0.02, dashed contours represent negative values. 57

Figure 1 :

 1 Figure 1: Von Newmann stability analysis: (a) second-order centered scheme, (b) present scheme. : Adams-Bashforth stability region; : locus of ∆t Â. Upper half-planes: F o = 0.1 and CF L varies (increment 0.2). Bottom half-planes: CF L = 0.5 and F o varies (increment 0.05).

Figure 2 :Figure 3 :

 23 Figure 2: L 2 -norm of the difference between the solution (u-velocity or pressure) and the 'exact' solution as a function of the number of iterations performed to compute the given time interval. • : DivSC; • : DivSC approx ; : AdvSC

Figure 4 :Figure 5 :Figure 8 :Figure 10 : 5 6 1 h

 4581051 Figure 4: Root-mean-square of pressure as a function of the grid spacing at time t = 20a/u 0 . : P rms ∝ ∆x 4 ; • : DivSC; :AdvSC; △ :DivSC ρ

Figure 16 :

 16 Figure 16: DNS of a channel flow with T 2 /T 1 = 4 (Case C). Typical z -y plane. Top: Contours of w + -velocity: the range is [-4,+4], the increment is 0.4, dashed contours represent negative values. Bottom: Contours of temperature (non-dimensionalized by T 1 ): the range is [1,4], the increment is 0.15, dashed contours represent values smaller that (T 2 + T 1 )/2.

Table 1 :

 1 Numerical parameters for time stepping RK-CN and AB-CN.

		4 15	8 15	0
	RK-CN k=2	1 15	5 12	-17 60
	RK-CN k=3	1 6	3 4	-5 12
	AB-CN k=1,3 1 2	3 2	-1 2

Table 2 :

 2 Properties of the different algorithms. A scheme has a given property if the entry is checked.

	tion and pressure terms (in the momentum equation) and the projection
	step respectively. The table also specifies whether or not the exact state
	equation (4) is enforced.

3 are summarized in table 2. Recall that AdvSC, DivSC and DivSC ρ stand for the schemes discussed in sections 3.1, 3.2 and 3.3 respectively. DivSC approx denotes the case where the state equation Eq. (

4

) is not enforced exactly -see section 3.2. The particular conservation properties of each of the schemes considered are denoted by cross. The columns 'convective', 'pressure' and 'projection' refer to kinetic energy conservation with respect to the convec-

Table 3 :

 3 Spatial convergence rates for AdvSC and DivSC. For each entry,

	1D high-amplitude pulse of density in an inviscid flow, a 1D small amplitude

Table 4 :

 4 Description of the DNS's.

		01	182 -185 2855	6.1	6.1	≈ 0	≈ 0
	B	2	195 -164 2810	7.0	5.0	-0.018 +0.016
	C	4	211 -151 2818	8.2	4.2	-0.041 +0.029

  1.01, 2, 4, 6, 8 and 10. The algorithm DivSC ρ is found to be unstable for T 2 /T 1 greater than 4. However, the other algorithms AdvSC, DivSC and DivSC approx remain stable up to T 2 /T 1 = 10. Larger temperature gradients were not considered since most of the low-Mach number flows have temperature ratio less than 10 (piston engine, rocket engines, burners, etc ...). For cases T 2 /T 1 ≥ 6, the temperature field was clipped to values greater than the temperature imposed at the cold wall T 1 . Indeed, with this coarse grid, a large amount of From the fine grid computations discussed earlier in this section, the proposed algorithms are accurate enough to represent the near wall complexity of turbulent flows with strong heat transfer. From the coarse grid simulations, these algorithms can handle larger temperature ratios than the methods used in most previous studies.

	6 Conclusion Maximum of divergence du dx as a function of the amplitude A
	of the temperature perturbation. Conservation of kinetic energy was specified as an analytical requirement for : du 47 dx ∝ A; △ :DivSC ρ 6 Root-mean-square of pressure as a function of the Prandtl
	a proper set of discrete equations in the zero-Mach number limit. The pro-number. P rms is non-dimensionalized by its theoretical value posed algorithms are fourth-order accurate in space and dissipation-free. A as a function of T ′ , viz. 4 3R 2 e P 2 r ∂ 2 T ′ . Note that the ab-∂x 2 rms key ingredient is a variable coefficient Poisson equation to solve for pressure. scissa is P r -3/4. : exact solution Eq. (54); • :DivSC; This elliptic operator ensures that the proper constraint is applied to the ve-:AdvSC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
	locity field at each (sub-)step of the time integration procedure. Noticeably, 7 Root-mean-square of pressure as a function of the Prandtl
	the divergence-free constraint is recovered in the inviscid limit as required number. P rms is non-dimensionalized by its theoretical value from the low-Mach number approximation to the Navier-Stokes equations. It is shown that this feature is necessary to avoid violating conservation of 4 ∂ 2 T ′ as a function of T ′ , viz. 3R 2 e P 2 r ∂x 2 rms
	kinetic energy. A proper discretization of the non-linear terms is also pro-
	posed that can handle any density variation and still remains conservative in
	kinetic energy. The only violation in the kinetic energy conservation appears
	in the projection step and is related to the non-linear state equation. A fully
	conservative scheme is obtained if instead an approximate state equation is
	used.

energy is present at the highest wave numbers (a common feature with more realistic LES's). For these scales the dispersion phenomenon is important and may generate negative values of temperature if a non-positive scheme is used. The clipping to values greater than the cold wall temperature is justified since there is no physical mechanism which can decrease the temperature below this point. It emulates a TVD or ENO scheme which should be used for temperature. Since such techniques are well known but beyond the scope of this study, we preferred to use the simple clipping approach T > T 1 when necessary. . The abscissa is P r . :P rms ∝ P 2 r ; △ :DivSC ρ . . . . . . . . . . . . . . . . . . 49 8 Kinetic energy conservation error ( K 0 -K / K 0 ) as a function of the Courant number. : ∆t 3 behavior; : ∆t behavior; • : DivSC; • : DivSC approx . . . . . . . . . . . . . . 50 9 Time evolution of the global energy of the fluctuations in the computational domain. : linear stability theory (Suslov & Paolucci); • : < u ′2 >; △ : < v ′2 >; : < T ′2 >.

Acknowledgments

The author would like to thank Dr. Oleg Vasilyev for many helpful discussions during the course of this study. C. Pierce and C. McNeil are acknowledged for their remarks on a draft of this manuscript. 16 DNS of a channel flow with T 2 /T 1 = 4 (Case C). Typical z -y plane. Top: Contours of w + -velocity: the range is [-4,+4], the increment is 0.4, dashed contours represent negative values. Bottom: Contours of temperature (nondimensionalized by T 1 ): the range is [1,4], the increment is 0.15, dashed contours represent values smaller that (T 2 + T 1 )/2. 58

The temperature ratio is T 2 T 1 = 2.33, i.e. δT 2 = 0.4. For these conditions, the flow is linearly unstable (see Suslov & Paolucci [16]). The code predicts a reasonable (5 % error) energy growth rate even if the number of grid points is rather small in the direction normal to the wall. Note that a fairly long time (10h/u τ ) is needed for the mode to settle in. Once the transition phase is finished, the temperature and the two velocity component develop with exactly the same rate, as dictated by linear stability theory. The four P rms is non-dimensionalized by its theoretical value as a function of T ′ , viz. P rms is non-dimensionalized by its theoretical value as a function of T ′ , viz. : linear stability theory (Suslov & Paolucci);

List of Tables

List of Figures