Rayleigh–Bénard Convection in a Near-Critical Fluid Using 3D Direct Numerical Simulation

Convection in a fluid close to its gas-liquid critical point (CP) has been a subject of growing interest since the exhibition of the piston-effect (PE), this thermo-acoustic effect responsible for the fast thermal equilibrium observed in such a fluid in the absence of convection. In 1987, under microgravity conditions, Nitsche and Straub observed a fast and homogeneous increase of the temperature inside a spherical cell containing SF6 slightly above the CP when it was subjected to a heating impulse. This phenomenon was then explained theoretically (Zappoli et al., 1990; Onuki et al., 1990; Boukari et al., 1990) by the well-known critical anomalies, more precisely by the divergence of the thermal expansion coefficient and the vanishing of its thermal diffusivity when approaching the CP. Indeed, the heating of a cell containing a supercritical fluid (SCF) induces along the heated wall a thin thermal boundary layer in which density shows large variations because of the divergence of the thermal expansion coefficient; this thermal layer expands compressing adiabatically the rest of the fluid leading by thermo-acoustic effects (the socalled PE) to a fast and homogeneous heating of the bulk of the cell. Several experiments were carried out subsequently, mainly in microgravity (Guenoun et al., 1993; Straub et al., 1995; Garrabos et al., 1998) but also on Earth (Kogan & Meyer, 1998), and confirmed the existence of the PE.


Introduction
Convection in a fluid close to its gas-liquid critical point (CP) has been a subject of growing interest since the exhibition of the piston-effect (PE), this thermo-acoustic effect responsible for the fast thermal equilibrium observed in such a fluid in the absence of convection. In 1987, under microgravity conditions, Nitsche and Straub observed a fast and homogeneous increase of the temperature inside a spherical cell containing SF 6 slightly above the CP when it was subjected to a heating impulse. This phenomenon was then explained theoretically (Zappoli et al., 1990;Onuki et al., 1990;Boukari et al., 1990) by the well-known critical anomalies, more precisely by the divergence of the thermal expansion coefficient and the vanishing of its thermal diffusivity when approaching the CP. Indeed, the heating of a cell containing a supercritical fluid (SCF) induces along the heated wall a thin thermal boundary layer in which density shows large variations because of the divergence of the thermal expansion coefficient; this thermal layer expands compressing adiabatically the rest of the fluid leading by thermo-acoustic effects (the socalled PE) to a fast and homogeneous heating of the bulk of the cell. Several experiments were carried out subsequently, mainly in microgravity (Guenoun et al., 1993;Straub et al., 1995;Garrabos et al., 1998) but also on Earth (Kogan & Meyer, 1998), and confirmed the existence of the PE.
Since 1996, many experimental and numerical studies were devoted to the interaction between the PE and natural convection. The Rayleigh-Bénard configuration (bottom heating) received a particular attention (Kogan et al., 1999;Amiroudine et al., 2001;Furukawa & Onuki, 2002) because the hydrodynamic stability of the SCF in that case is governed by an interesting and non-common criterion. Owing to the PE, the thermal field exhibits a very specific structure in the vertical direction. A very thin hot thermal boundary layer is formed at the bottom, then a homogeneously heated bulk settles in the core at a lower temperature, and at the top, a cooler boundary layer is formed in order to continuously match the bulk temperature with the colder temperature of the upper wall. The linear analysis, carried out by Gitterman and Steinberg in 1970, showed that the hydrodynamic stability of these thermal boundary layers, when subjected to a gravity field, depends on the interaction between two stability criteria which, for a normally compressible fluid, are separately available at very different space scales: on one hand, the classical Rayleigh criterion, derived from the Boussinesq approximation, hence available at small space scales, and on the other hand, the Schwarzschild criterion, usually encountered in atmospheric science, where the stabilizing effect of the hydrostatic pressure becomes appreciable. Indeed, because of the divergence of the isothermal compressibility of a SCF, the Schwarzschild criterion becomes available at small space scales; this was proven theoretically (Gitterman & Steinberg, 1970b;Carlès & Ugurtas, 1999), experimentally (Kogan & Meyer, 2001), and numerically (Amiroudine et al., 2001). Taking advantage of the interaction between those two stability criteria, a numerical study (Accary et al., 2005a) showed that, in spite of convection onset in the thermal boundary layers according to the classical Rayleigh criterion, a reverse transition to stability through the Schwarzschild line is possible without any external intervention. The hydrodynamic stability of the thermal boundary layers developed in this configuration has been exhaustively investigated (Accary et al., 2005b), and recently a numerical study (Accary et al., 2009) focused on the convective regime of the flow. Because of the particular physical properties of the fluid in the vicinity of the CP, the convective regime of the Rayleigh-Bénard problem is turbulent for unusually low intensities of heating (~mK). In this last study, 3D direct numerical simulations are carried out for Rayleigh numbers varying from 2.68×10 6 up to 160×10 6 . For a perfect gas (PG), this range of Rayleigh numbers corresponds to the transition between the soft and the hard turbulence; however, this is not always the case for the SCF because of its strong stratification induced by its high isothermal compressibility.
In § 2, the problem under consideration is presented. In § 3, the mathematical model is described together with the acoustic filtering of the Navier-Stokes equations. The details of the numerical method are presented in § 4 and the simulation conditions are mentioned. In § 5, several aspects of the Rayleigh-Bénard convection in a near-critical fluid are reported: the hydrodynamic stability of the thermal boundary layers, the convection onset and the beginning of the convective regime, the steady-state turbulent regime, details of the temperature and the dynamic fields, and the global thermal balance of the cavity. In § 6, a comparison with the case of the PG is presented at equal Rayleigh number. Finally, the chapter is concluded in § 7.

The problem under consideration
We consider a SCF in a cube-shaped cavity (of height H' = 10 mm) subjected to the earth gravity field g' = 9.81 m.s -2 (Fig. 1). The horizontal walls are isothermal while the sidewalls are insulated, and no-slip conditions are applied to all the walls. Initially the fluid is at rest, in thermodynamic equilibrium at a constant temperature T' i slightly above the critical temperature T' c , such that T' i = (1+)T' c , where  << 1 defines the non-dimensional proximity to the CP. Under the effect of its own weight, the fluid is stratified in density and in pressure, with a mean density equal to its critical value ' c . While maintaining the top wall at its initial temperature T' i , the temperature of the bottom wall is gradually increased (during one second) by T' (a few mK).
www.intechopen.com Fig. 1. Geometry of the cube-shaped cavity and the velocity and temperature conditions applied to the boundaries. The vertical axis z' is co-linear with the acceleration due to the earth gravity g'. Since the first seconds of heating, the temperature field is vertically stratified, divided in three distinct zones: two thermal boundary layers and the bulk of the cavity.

Equations governing near-critical fluid buoyant flows
The mathematical model for a SCF flow (Zappoli, 1992) is described by the Navier-Stokes and energy equations written for a Newtonian and highly conducting van der Waals fluid. Continuity: Van der Waals: Where P' is the pressure, T' is the temperature, and ' is the density. v'(u',v',w') is the velocity, g' = (0,0,-g'), e' is the internal energy, and ' is the viscous energy dissipation. ' is the dynamic viscosity, ' is the thermal conductivity, r' is the PG constant, a' and b' are respectively the energy parameter and the co-volume related to the critical coordinates T' c and ' c by: b' = 1/(3' c ) and a' = 9r'T' c /(8' c ).
In spite of its simplicity, the van der Waals' equation of state gathers the required conditions for the existence of the CP and yields to a critical divergence 1 as (T'/T' c -1) -1 of the thermal expansion coefficient ' P , of the isothermal compressibility ' T , and of the heat capacity at constant pressure C' P . The critical divergence of the thermal conductivity is described by: where  = 0.75 and ' 0 is the thermal conductivity for a PG. The heat capacity at constant volume C' V and the dynamic viscosity are supposed to be constant and equal to those of a PG, C' V0 and ' 0 respectively. With the van der Waals' equation of state, the expression of the internal energy is given by: In order to make the variables dimensionless, T' c , ' c , and r'' c T' c are used respectively as representative scales of the thermodynamic variables T', ', and P'. The independent variables of length x'(x',y',z') and time t' are scaled respectively by the height of the cavity H' and the PE time-scale 2 given by t' (Zappoli, 1992;Zappoli et al., 1999). Hence, the representative scale of velocity is V' PE = H'/t' PE . This scaling introduces the Reynolds numbers Re =, the Froude number Fr = (V' PE ) 2 /(g'H'), the Prandtl number based on the properties of the PG assumption (Pr 0 = ' 0 C' P0 /' 0 ), and the Mach Note that the PE time-scale obtained by (Onuki et al., 1990) is given by where D' T is the thermal diffusivity and  = C' P /C' V is the specific-heat ratio. Adapted for a van der Waals' gas 3 and in the assumption that  <<1, the PE time-scale

The acoustic filtering of the governing equations
Despite its high isothermal compressibility, the sound speed c' i n a S C F , d e f i n e d b y c' 2 = C' P / C' V ' T -1 , does not vanish at the CP according to the van der Waals' equation of state, indeed C' P /C' V and ' T diverge with the same critical exponent of -1, which allows the acoustic filtering of the equations. In the basic assumption that Ma<< 1, all the primary dimensionless variables of the problem  = t (v, T, P, ) can be expanded in series of Ma 2 (Paolucci, 1982) as follows: and O(Ma 2 ) parts of the governing equations resulting from this expansion that need to be solved are given by: O(1) continuity: 1 The real critical exponents (which are the same for all fluids) differ from those obtained from the van der Waals' equation of state that remains a good approximation to carry out qualitative studies. 2 The PE time-scale is the time necessary the PE to homogenize the temperature in the core of the cavity, this time-scale is between the acoustic time scale (H'/c') and the thermal diffusion one (H' 2 /D' T ) where c' is the sound velocity in the SCF and D' T is the thermal diffusivity.
O(1) van der Waals: In these equations, v(u,v,w), T, P, and  refer to the O(1) of the dimensionless variables, the superscript (0) has been omitted for conciseness. P th and P hyd are respectively the thermodynamic pressure (homogeneous in space but time varying according to the O(1) momentum equation) and the time independent hydrostatic pressure (P (0) = P th + P hyd ). a = 9/8 and b = 1/3 are the dimensionless parameters of the van der Waals' equation of state, e g = (0,0,-1). Before the heating begins, the initial dimensionless density distribution  i , the initial thermodynamic pressure P thi , and the hydrostatic pressure P hyd are obtained from the initial thermodynamic and static equilibrium with the constraint of a dimensionless mean density equal to 1. This results in: This low Mach number approximation (adapted to SCF buoyant flows) differs from the classical one where  i = 1 and consequently P hyd = 0. Indeed, owing to the divergence of the isothermal compressibility of the SCF, the hydrostatic pressure induces density variations ( i -1) comparable to those resulting from a weak heating. This has been done by keeping the buoyancy term ( 0 Ma 2 /Fr) (0) eg in the leading order O(Ma 2 ) of the momentum equation (Eq. 7) while in the classical low Mach number approximation, this term is shifted to the O(Ma 4 ) order. It has been shown (see § 5.2) that this modification is essential for a correct prediction of the convection onset in the thermal boundary layers (Accary et al., 2005c).

Numerical method
In describing the numerical method used in this analysis, it is assumed that the reader is familiar with the basis of the standard finite volume method and with the velocity-pressure coupling algorithms extensively reported in (Patankar, 1980). In this section, we will draw the outlines of the method the space and the time accuracies, the velocity-pressure coupling, the linear systems solvers, and the solver performance.

Space and time discretization
The computational domain is subdivided into a number of cells using a wall refined mesh for a better description of the solution in the boundary layers, the mesh is refined in the vicinity of the walls; as one moves away from the wall, the control volume size increases according to a geometric progression. If N is the number of cells in a direction (x for example), the dimensionless positions of the cells interfaces between the wall and the center of the computation domain (i.e. for 0 < x(i) ≤ 0.5) would be: For 0.5 ≤ x(i) < 1, cells interfaces position are obtained by symmetry with respect to x = 0.5. Depending on the value of q (q≥1), this mesh refinement is termed 'power q law' type; q=1 provides a uniform mesh, the simulations were carried out with q = 2. The variables location is staggered: the scalar variables are stored at the cells centers while the velocity components are defined at the midpoints of the cells faces perpendicular to the velocity direction. This staggering practice avoids the high-frequency noise in the solution resulting from the well-known problem of the zigzag pressure filed which would be made up of arbitrary values of pressure arranged in a checkerboard pattern (Patankar, 1980). In return, this staggering has no adverse consequences in the simple rectilinear domain considered here.
The convection-diffusion transport equation of a variable  to be solved in computational fluid dynamics can be written as: Where Γ  is the diffusion coefficient and S  is a source term. While integrating Eq. 12 over the control volume (CV) of a discrete variable  p and over a time step, the value of a variable  and its gradients are assumed to be constant on its CV faces; therefore the space accuracy of the method is already limited to the second order, depending on the interpolation scheme in the direction perpendicular to the considered face used to approach the value of  and its derivative.
The time integration is fully implicit providing the method a non-conditional stability as far as the time step is concerned; and to allow large time scale simulations, the unsteady terms are approached by a standard Euler time scheme with four time levels, leading to a third order truncation error in time. The time integration of Eq. 12 at the instant t n and for a uniform time step t is done as follows: Careful space discretization and integration of the transport equation is needed to reach the second order space accuracy ceiling of the method particularly for a non-uniform mesh. This concerns the integration of the source and the unsteady terms, and the evaluation of the variable  and its gradients at the faces of a CV. A linear interpolation is used to evaluate the density at the faces of a CV, while a harmonic mean is considered for the thermal conductivity as recommended in (Patankar, 1980). Figure 2 shows the grid structure and the notations in one dimension. For a second order space integration, the integrand should be localized in the centre of the CV; it is the case for the scalar variables (see Fig. 2), while for the velocity components a neighbor contribution must be considered when the grid is not uniform. For example, if  refers to a velocity component, its value * w at the centre of the velocity CV ( Fig. 2) can be written as: A central difference (CD) scheme approaches the diffusion terms. For the velocity components, since CV faces are localized at midway between two consecutive velocity nodes (see Fig. 2), the standard two-point formulation provides a second order approximation of the velocity gradients at the faces of a CV. It is not the case for a scalar variable on non-uniform mesh, where a three-point scheme is necessary to reach the second order accuracy. For example, the gradient of a scalar  at the face p can be written as: In this scheme, a correction is added to the standard CD formulation, a correction depending on (δ 2p -δ 1p ), hence that vanishes for a uniform mesh.
A second order hybrid scheme (SHYBRID) is used for the convection terms (Li & Rudman, 1995). It uses a four-point formulation to interpolate the values of a variable at the faces of its CV, and combines the QUICK, the second order upwind (SOU), and the CD schemes whose respective weights in the formulation depend on the local Peclet number which is the ratio of the convection flux to the diffusion one. The value  p of a variable  at a face p is determined as follows: where, The CD, SOU, and QUICK schemes fall into this formulation for appropriate choices of qp and q + p ; for the CD scheme q + p = qp = 0, for the SOU scheme qp =  and q + p = 1-, and for the QUICK one qp = δ 1p /(δ 1p +δ 2p +δ 3p ) and q + p = (1-)δ 2p /(δ 1p +δ 2p +δ 4p ). The values of qp and q + p are automatically adjusted during the simulation according to the local Peclet number (Pe p ) in order to minimize the potential oscillations by minimizing the remote-nodes contributions while maintaining positive neighbor-nodes coefficients (see Patankar, 1980, 2 nd basic rule). This minimization (Li & Rudman, 1995) A quick inspection of these choices shows that when the transport is diffusion-dominated (Pe p very small) SHYBRID scheme becomes the CD, when the transport is convectiondominated (Pe p very large or infinity) SHYBRID approaches the SOU scheme, and between these two limits it may go through the QUICK scheme for certain values of Pe p . Thus, SHYBRID is a stable second order accurate scheme for a wide range of the Peclet number. In order to suppress non-physical oscillations when predicting solutions with sharp gradients (Li & Rudman, 1995), a flux-correction transport was necessary; it was however useless for the gradient magnitudes encountered in our study.
Integrating Eq. 12 for  p , the discrete transport equation may be written after some manipulation (Patankar, 1980), as: where a P and a nb are discretization coefficients, S c is the discrete source term, and the subscript nb designates the six direct neighbors of the node P, any remote-node contribution resulting from second order space discretization is included in S c . The linear system (Eq. 18) is solved using iterative methods (Barrett et al., 1994); the pressure symmetric equation is solved using the Conjugate Gradient method with Jacobi preconditioning, while the Bi-Conjugate Gradient Stabilized with the same preconditioner is used for the other nonsymmetric transport equations. In addition, the pressure equation is preconditioned using a SCF equivalent of the artificial compressibility method proposed by (Chorin, 1997). In spite of the computational optimization of these solvers, most of the computation effort was spent on solving the linear systems and especially the pressure one.

The coupling algorithm
The pressure splitting in the low Mach number results into two more variables in the governing equations, the hydrostatic pressure P hyd and P th the thermodynamic one. P hyd is time-independent and is given by the initial stratification (Eq. 10), while P th (constant in space) can be determined at any moment using the conservation of the total mass whose dimensionless value is equal to 1. www.intechopen.com In order to determine P th , one must provide the function  =  (P th ,T); thus at each time step and at each iteration k of the velocity-pressure coupling algorithm, after computing the temperature field, the density is linearized using the van der Waals' equation as follows: The thermodynamic pressure at iteration k is computed from the conservation of the total mass of the fluid, since: Thus for each time step, a global iterative process (Fig. 3) consists first in solving the energy equation (since the flow is temperature driven), updating the thermodynamic pressure (Eq. 21) and the density (Eq. 20) using the linearized equation of state, and then solving the dynamic field. The global convergence is assumed obtained when the L  -norm of all governing equations residuals (momentum, energy, and state) reach an imposed stopping criterion. The velocity-pressure coupling is treated by a PISO algorithm (Jang et al., 1986); at each sweep, a velocity prediction, a pressure prediction, and a velocity correction are performed; at this stage, a second correction of the dynamic field is useless since the density will be severely perturbed in the next iteration. PISO algorithm determines the pressure field using a pressure equation and requires no pressure correction that introduces instability into the convergence process of unsteady solutions. More details about the numerical method can be found in (Accary et al., 2006) where the code has been thoroughly validated using an artificial analytical solution and on several benchmark problems of natural convection.

Hydrodynamic stability of the thermal boundary layers
As mentioned earlier, because of the PE, the temperature field is stratified vertically with three distinct zones since the first seconds of heating: the hot boundary layer, the cold boundary layer and the bulk of the cavity. Regardless the considered heating, as long as the flow is dominated by the diffusion and by the PE, the thermal boundary layers grow as  Figure 4(b) shows the corresponding density profiles; we notice that the density variations induced by the heating are comparable to those resulting from the hydrostatic pressure, which justifies the adaptation of the low Mach number approximation by including the fluid stratification in the model.  (Zappoli & Durand-Daubin, 1994) as a result of the PE action that increases the temperature of the core. For a SCF diffusing-layer, the local Rayleigh number based on h and T is given by (Gitterman & Steinberg, 1970b;Carlès & Ugurtas, 1999): To account for the high compressibility of the SCF, the classical expression of the Rayleigh number is modified in Eq. 22 by the adiabatic temperature gradient (T' a /H' = g'' P T' i /C' P ) obtained by moving a fluid particle along the hydrostatic pressure gradient. This term, that can be neglected for a normally compressible fluid, represents the stabilizing contribution of the Schwarzschild criterion commonly encountered for large air columns, and according to which the fluid layer is stable if: In the considered model, the adiabatic temperature gradient T' a /H' = 0.34 mK/cm and does not depend on the proximity  to the CP since ' P and C' P have the same critical exponent of -1. To better estimate the interaction between natural convection and stratification, the normalized intensity of heating of the bottom wall T is henceforth expressed in terms of T a = T' a /T' c .

Fig. 5(b) by the curve 180×t 3/2 ×e t ×erfc(t ½ ), and we can easily prove at long time scales that erfc(t ½ ) ~ e -t t -½ , which explains the linear time evolution of Ra corr (h,T). When the local
Rayleigh number exceeds the critical value of about 1100 (Chandrasekar, 1961), the hot boundary layer becomes unstable 4 . Convective cells start to get organized along the bottom plate; figure 6(b) enables the visualization of these vortical structures using the Q-criterion 5 (Dubief & Delcayre, 2000) along the bottom plate. Then, the intensity of these vortices rises exponentially with time; this can be easily seen on the time evolution of the mean enstrophy in the hot boundary layer shown in Fig. 5(b) and defined by: The intensity of these convective cells keeps rising until producing enough amount of convective transport to deform the isotherms causing the collapse of the thermal boundary layers. For T = 3T a , the collapse of the hot boundary layer occurs around t = 120 and corresponds to the symbol (□) in Figs. 5(a) and 5(b), afterwards the convective regime starts and the definition of the hot boundary layer holds no more. The cold boundary layer developed along the top plate is governed by the same mechanisms and its hydrodynamic stability depends on the same criterion (Accary et al., 2005b); the same scenario occurs for the cold boundary layer. In Fig. 7 figure 7 shows the evolution curves T(h) for the hot boundary layer (the solid blue lines) until the beginning of the convective regime, which corresponds to the symbol (□). Figure 7 shows also results obtained in a 2D approximation with periodic vertical boundaries (the dashed red lines) (Accary et al., 2005b). The boundary effects induced by the presence of the lateral walls in the 3D case accelerates the development of convection and, at equal intensity of heating, the convective regime is reached earlier in comparison with the 2D approximation with periodic vertical boundaries. For low intensities of heating, practically for T  0.72T a , once the hot boundary layer has become unstable, the intensity of the convective cells rises exponentially with time until deforming the isotherms. However, this deformation is not large enough to induce the collapse of the hot boundary layer that keeps growing and the curve T(h) crosses the Schwarzschild line back into the stable zone again and a reverse transition to stability obtained without any external intervention (Accary et al., 2005a). This phenomenon requires that the thermal boundary layers grow enough without reaching the centre of the cavity (in order to avoid their interaction); a height of 1.5H' at least is needed in this case. Fig. 7. Evolution of the temperature difference T across the hot boundary layer as a function of its thickness h. Time evolves in the arrows' direction and the symbols (□) correspond to the beginning of the convective regime. The neutral stability line was derived from Eq. 20. The 2D results were obtained using a cavity of height 1.5H' with periodic vertical boundaries; for T  0.72T a , a reverse transition to stability is obtained though the Schwarzschild line.
www.intechopen.com Figure 8 shows a comparison for the evolution curves T(h) between the low Mach number approximation adapted to the SCF flow (ALMN) described in § 3.2 and the classical one (CLMN) obtained from the model by setting  i = 1 in Eq. 10. The simulations are carried-out in a 2D approximation but with periodic vertical boundaries in order to suppress the disturbance resulting from the side walls. The straight dashed line with a slop of (-3) represents the classical Rayleigh criterion obtained from Eq. 22 by removing the adiabatic temperature gradient term. We notice that for a relatively strong heating (T > 3T a ) the flows predicted with both models are similar but not identical and the collapse of the thermal boundary layers occur at about the same time; this is due to the fact that a strong heating covers the effects of stratification, but the weight of this latter becomes more significant as the heating decreases. If we consider for example the case T=0.9T a , the collapse of the hot boundary layer is observed with about 90t' PE of time gap between the two models. For T=0.9T a , figure 9 shows the temperature fields obtained 100t' PE after the collapse of the hot boundary layer. Compared to Fig. 9(a) (ALMN), the thermal plumes are much more developed in Fig. 9(b) (CLMN) since their motion is not hindered by the stratification which, when taken into account, prevents the free growth of the plumes. However for the considered heating intensity (T=0.9T a ), the thermal plumes manage to deform the temperature field giving rise the slowly moving structures. As the heating gets weaker (T < 0.6T a for example), the buoyant force being not strong enough for pulling the fluid particles through the hydrostatic pressure gradient, the hot boundary layer predicted with the ALMN approximation remains stable. In return, not including the stratification, the CLMN model is unable to take account of this stabilizing effect and persists in predicting a convective instability (according to the classical Rayleigh criterion) provided that the height H' of the cavity allows enough growth of the thermal boundary layers.

The beginning of the convective regime
The convective regime starts with several plumes rising from within the thermal boundary layers as shown in Fig. 10(a). These plumes are encircled by donut-shaped structures shown by the Q-criterion in Fig. 10(b). Convection improves the heat transfer between the isothermal walls and the bulk of the cavity, resulting into a faster thermal balance in the whole fluid volume. For all the heating cases that we considered, the hot boundary layer has always become unstable before the cold one. As the heating increases, convection is triggered earlier since the instability criterion (Ra corr (h,T) > 1100) is satisfied earlier; consequently, the thickness of the thermal boundary layer is smaller when the convection arises and the size of the convective structures decreases as shown in Fig. 11. A detailed study of the size of the convective structures has been done in a 2D approximation in (Accary et al., 2005b).

Transition to turbulence
In the convective regime of the flow that follows the convection onset, the Rayleigh number, based on the total height H' of the cavity and on the temperature difference T' between the isothermal walls (Eq. 25), becomes a better indicator of the regime of the flow.
For T < T a , the Rayleigh number obtained from Eq. 25 is negative; this, however, does not prevent convection to arise in the thermal boundary layers when the local Rayleigh number (Eq. 22) exceeds 1100. But for T > T a , for example for T = 1.5T a , the term in front of the parentheses in Eq. 25, which diverges as  -1.5 , is very large and results in a Rayleigh number of 2.68×10 6 , while for a PG, the Rayleigh number is directly proportional to T.
The turbulent Rayleigh-Bénard convection is characterized by a statistically steady state of heat transfer. In the considered configuration, the settlement of the turbulent regime may be identified on the time evolution of the mean Nusselt numbers on the isothermal walls given by: For T = 7.5T a which corresponds to Ra corr = 80×10 6 (Eq. 25), figure 12 shows the time evolution of the mean Nusselt numbers on the bottom wall and of the top one. The convection onset is easily identified by the improvement of the heat transfer corresponding to the increase in the mean Nusselt numbers that stabilize afterwards around almost the same value, which indicates the settlement of the turbulent flow. Figure 13(a) shows the temperature field obtained in the turbulent regime. We notice first the appearance of crest-like patterns defining on the isothermal walls flat regions where the temperature is almost homogeneous in the (x,y) plan, we notice also the spreading of the isotherms along the adiabatic walls. Figure 13(b) shows the chaotic flow that takes place in the turbulent regime. The vortical structures have no particular shape; the tubular and toroïdal structures obtained at the beginning of the convective flow have completely disappeared. In order to better estimate the size of the vortical structures and its time evolution, a discrete Fourier transformation 6 of the vertical velocity component w has been carried out in both x and y directions. Along the line (y = y 0 , z = z 0 ) and for a wavelength H'/k associated to the mode k, the Fourier coefficient of w(x,y 0 ,z 0 ) is given by: Once the coefficients W x (k, y 0 , z 0 ) are computed for all (y 0 , z 0 ), the mean contribution of the mode k to the field of w is determined by: Figure 14 shows the contribution of the different modes to the spectrum of the component w in the x and y directions at the beginning of the convective regime and in the turbulent one. We notice an important contribution of small wavelengths ranging between H'/11 and H'/4 at the beginning of convection (at t = 89.25, see Fig. 12). But as time goes by, the spectra of w show a much higher contribution of large wavelengths exceeding sometimes half the of the cavity width. Similar results were obtained for the horizontal velocity components, u and v. Thus, the turbulent flow consists mainly of large vortical structures. Figure 15 shows cuts of the temperature field in the vertical median plans of the cavity with the corresponding velocity fields that confirm the presence of large convective structures in the steady-state turbulent regime. We notice that the temperature field consists mainly of two unstable thermal boundary layers exchanging heat and mass with the bulk of the cavity in  which the convective activity induces a quasi-homogeneous temperature. Figure 16 shows the time evolution, along the vertical axis of the cavity (x = y = 0.5), of the velocity magnitude and of the temperature at the free boundaries of the thermal layers (z = 0.05 and z = 0.95) 7 and at the centre of the cavity; the velocity components have the same order of magnitude. Figure  16(a) underlines the chaotic convection that takes place in the whole fluid volume; the velocity has been monitored at 25 different points of the cavity and confirms that chaotic behavior. In the steady-state turbulent regime, figure 16(b) shows a slight difference of the time-averaged temperature between positions z = 0.05 et z = 0.95, which reveals the existence of a temperature gradient in the bulk of the cavity that will be investigated in § 5.5.

The global thermal balance of the cavity
The steadiness of the mean Nusselt numbers on the isothermal walls (Fig. 12, turbulent regime) reflects the settlement of a statistically steady-state heat transfer across the cavity. However, figure 17 reveals the strong non-uniformity of the Nusselt numbers distributions on the isothermal walls. These patterns are directly related to those of the temperature field: the Nusselt number's minima are reached under the crest-like patterns shown in Fig. 13(a), while the maxima are obtained inside the cells determined by those patterns. These cells are thus characterized by very thin thermal boundary layers; for the temperature field shown in Fig. 13(a), the minimal normalized thicknesses of the thermal boundary layer were about 0.014 for the hot boundary layer and 0.012 for the cold one and were obtained where the distributions of the Nusselt numbers reach their maxima. Despite the strong non-uniform distributions of the Nusselt numbers, in the steady-state turbulent regime, the mean Nusselt numbers on both isothermal walls fluctuate around the same value.
For different intensities of heating and hence Rayleigh numbers, figure 18(a) reports the mean Nusselt numbers (the filled circles). For a PG, the experimental results (Poche et al., 2004) and those issued from a scaling theory (Siggia, 1994) show that Nusselt number behaves as Ra 2/7 . This behavior can be observed for the of the Nusselt number corrected by the adiabatic temperature gradient (Kogan & Meyer, 2001), given by: For T >> T a , Nu corr  Nu; but as the intensity of heating decreases, the corrected expression of the Nusselt number (the filled squares in Fig. 18(a)) enables the retrieval of the Ra 2/7 law. However, it should be reminded that the effective heat transfer is described by the classical expression of the Nusselt number given by Eq. 26, not by the corrected one.  At the global thermal balance of the cavity and for all intensities of heating, figure 18(b) reveals the existence of a mean temperature gradient in the bulk of the cavity equal to the adiabatic temperature one. This is a natural structure of the mean temperature field that ensures the minimal temperature gradients in the thermal boundary layers with the constraint of a globally stable bulk of the cavity. Indeed, if the mean temperature gradient in the bulk of the cavity were larger than the adiabatic temperature one, the bulk of the cavity would lose its hydrodynamic stability. In return, if the mean temperature gradient in the bulk of the cavity were smaller than the adiabatic temperature one, the bulk of the cavity would be 'too' stable, but this would increase the temperature gradients in the thermal boundary layers.

Comparison between a SCF and a PG, effects of stratification
The comparison between the Rayleigh-Bénard convection in a SCF and that in a PG is carried out here in the 3D case for a Rayleigh number of 2.6810 6 for which the density stratification of the SCF affects clearly the development of convection (T = 1.5T a ). The mathematical model described in § 3 was adapted to the PG case, mainly by setting a = b = 0,  = 0, and  = 0 in Eqs. 5 to 9, and by choosing reference values of temperature and density compatible with the PG assumption, these were set to 300 K and 1.8 Kg.m -3 respectively. An intensity of heating T' = 5 K was applied to bottom wall in the PG case and the height H' of the cavity, deduced from the classical expression of the Rayleigh number, is equal to 13.8 cm. A mesh of 100 3 and a time step of 0.125s have been used. Fig. 19. Time evolution of the mean Nusselt numbers on the isothermal walls obtained in the 3D case for a Rayleigh number of 2.6810 6 for a SCF (T = 1.5T a ) and for a PG. The curves were shifted by 2s for the PG to show prominently the first peak. For the SCF, the first peak of the mean Nusselt number on the bottom wall reaches the value of 370, and the beginning of convection at about 60 s is consistent with the result shown in Fig. 7, where the convective regime starts at t' = 120t PE = 58.9 s. Figure 19 shows the time 8 evolution of the mean Nusselt numbers for the SCF and for the PG 9 . The large temperature gradients obtained at the very first seconds of heating in the case of the SCF are responsible for the very high peak of the Nusselt number. For the SCF, figure 19 reports very similar evolutions of the mean Nusselt numbers on both isothermal walls. By contrast for the PG, while the mean Nusselt number on the bottom wall shows a similar behavior to that of the SCF during the diffusive regime, no heat transfer is detected on the top wall (Nu c = 0) until the beginning of convection. Because the PE is practically inexistent for the PG, the heat transfer is only activated on the top wall when the thermal plumes rising from the hot boundary layer reach it. Even though the Prandtl number is about 18 times smaller 10 (Verzicco & Camussi, 1999), convection in the PG is much more developed than in the SCF at the same Rayleigh number, as shown by Fig. 20. The fluctuating time evolution of the mean Nusselt numbers for the PG results from this intense convective activity. By contrast, the trace of the diffusion-dominated temperature field ( Fig.  20(a)) obtained for SCF due to its strong stratification is visible on the time evolution of the mean Nusselt numbers after the convection onset. Under these conditions, the global thermal balance of the cavity is mainly achieved by diffusion at long time scales because of the critical vanishing (as  ½ ) of the thermal diffusivity of the SCF. We notice finally that even though the temperature field of the SCF is diffusion dominated while it is convectiondominated for the PG, the corrected mean Nusselt number at the global thermal balance of the cavity is the same in both cases.

Conclusions
In this chapter, the mathematical model for SCF buoyant flows with the appropriate acoustic filtering has been recalled, then a description of the different stages of the SCF flow in a cubeshaped cavity heated from below were reported from the first seconds of heating until the settlement of a statistically steady-state of heat transfer, and this for Rayleigh numbers ranged from 2.68×10 6 up to 160×10 6 . While the scenarios of the convection onset and disappearance (reverse transition to stability) can be observed in a 2D approximation, the convective regime and the transition to turbulence requires 3D simulations. At the beginning of convection, tubular convective structures appear inside the thermal boundary layers while the thermal plumes are encircled by toroïdal vortical structures; the size of these structures decreases as the intensity of heating increases. In the turbulent regime, the convective structures grow until their size exceeds half of the cavity, and create on the isothermal walls several cells where an intense heat transfer takes place. Despite the non-homogeneous heat transfer on the isothermal walls, the steadiness of the mean Nusselt numbers around the same value reflects the global thermal balance of the cavity. The relation between that equilibrium Nusselt number and the Rayleigh number obtained for a PG (Nu ~ Ra 2/7 ) is applicable to the SCF, provided that the adiabatic temperature gradient is taken into account in the expressions of both numbers. In the turbulent regime, the temperature field consists mainly of two unstable thermal boundary layers and a bulk characterized by a mean temperature gradient equal to the adiabatic temperature one. For relatively high intensities of heating (T >> T a ), the global thermal balance of the cavity is achieved by a chaotic convection invading in the whole fluid volume. By contrast for weak intensities of heating (T  T a ), the strong density stratification, due to the high isothermal compressibility of the fluid, prevents the free development of convection whose penetrability is dramatically reduced; in this case, the thermal balance of the cavity is mainly achieved by diffusion and therefore on long time scales. Finally, the comparison between the SCF and the PG for the same Rayleigh number showed two major differences. The first, related to the PE, is the absence of heat transfer on the top wall for the PG until the beginning of convection; while for the SCF, the time evolutions of the mean Nusselt numbers on both isothermal walls are similar. The second, related to the stratification of the SCF and thus only encountered for T  T a , is the diffusion-dominated thermal balance of the cavity for the SCF, while it is convection-dominated for PG.