Xavier Nicolas
email: xavier.nicolas@univ-paris-est.fr

Marc Medale
email: marc.medale@polytech.univ-mrs.fr

Stéphane Gounand
email: stephane.gounand@cea.fr

Stéphane Glockner
email: glockner@enscbp.fr

Benchmark solution for a three-dimensional mixed convection flow -Detailed technical report

A solution to a benchmark problem for a three-dimensional mixed convection flow in a horizontal rectangular channel heated from below and cooled from above (Poiseuille-Rayleigh-Bénard flow) is proposed. This flow is a steady thermoconvective longitudinal roll flow in a large aspect ratio channel at moderate Reynolds and Rayleigh numbers (Re=50, Ra=5000) and Prandtl number Pr=0.7. The model is based on the Navier-Stokes equations with Boussinesq approximation. We propose a reference solution resulting from computations on large grids, Richardson extrapolation (RE) and cubic spline interpolations. The solutions obtained with finite difference, finite volume and finite element codes are in good agreement and reference values for the flow fields and the heat and momentum fluxes are given up to 4 to 5 significant digits. Some difficulties in the use of RE are highlighted due to the use of mixed Dirichlet and Neumann thermal boundary conditions on the same wall. The observed convergence orders of the numerical methods with RE are then discussed from the viewpoint of this singularity. A correction to the Taylor expansion involved in the RE formalism is proposed to take into account the singularity and to explain the majority of the RE behaviors observed. The results of the present study are published in two papers in Numerical Heat Transfer, Part B [1, 2].

Context and objectives

Mixed convection flows in channels of rectangular cross section are encountered in many industrial applications: thermal and chemical reactors, chimneys, solar collectors, thermal insulation of buildings, heat exchangers, etc. More specifically, the Poiseuille-Rayleigh-Bénard (PRB) configuration (i.e. mixed convection flows in horizontal rectangular channels heated from below) is representative of rectangular Chemical Vapor Deposition (CVD) reactors and of air flows in the cooling of printed electronic circuit boards, among others [START_REF] Nicolas | Numerical simulation of thermoconvective flows and more uniform depositions in a cold wall rectangular APCVD reactor[END_REF][START_REF] Wang | Instability and heat transfer in mixed-convection flow in a horizontal duct with discrete heat sources[END_REF][START_REF] Nicolas | Bibliographical review on the Poiseuille-Rayleigh-Bénard flows: the mixed convection flows in horizontal rectangular ducts heated from below[END_REF]. However, to the authors' best knowledge, three-dimensional benchmark numerical solutions of mixed convection flows in rectangular channels do not exist in the literature. In numerous numerical studies of PRB flows, for instance, numerical codes are only validated by comparisons with experimental data, particularly those of Chiu and Rosenberger [START_REF] Chiu | Mixed convection between horizontal plates -1. Entrance effects[END_REF] and Ouazzani et al. [START_REF] Ouazzani | Etude expérimentale de la convection mixte entre deux plans horizontaux à températures différentes -2[END_REF], in which laser Doppler anemometry measurements of local velocities are given. However, a precise agreement between the experimental and numerical results is hard to obtain because the thermal and dynamical initial and boundary conditions are not perfectly controlled experimentally and are very hard to accurately reproduce numerically [START_REF] Benzaoui | Efficient vectorized finite difference method to solve the incompressible Navier-Stokes equations for 3D mixed convection flows in high aspect ratio channels[END_REF]. The objective of the present paper is therefore to propose a first 3D numerical benchmark solution that can be used to validate numerical codes for the computation of thermoconvective instabilities in open rectangular channels.

This benchmark exercise was first proposed in the framework of the French Heat Transfer Society (SFT) by several laboratories involved in the numerical analysis of thermoconvective flows in closed cavities and open channels. A call for contributions was published in 2006 [START_REF] Medale | CALL FOR CONTRIBUTIONS: Towards numerical benchmark solutions for 3D mixed convection flows in rectangular channels heated from below[END_REF]. Initially, two configurations of PRB flows, covering two different flow ranges, were chosen. The first one concerned the present steady longitudinal roll flow in a large aspect ratio channel (A=L/H=50, B=W/H=10) at Reynolds number Re=50, Rayleigh number Ra=5000 and Prandtl number Pr=0.7. The second one concerned a fully-established space and time periodic transverse roll flow in a small aspect ratio channel (A=25, B=4) at small Reynolds number Re=0.1, Ra=2500 and Pr=7. In this paper, only the solution of the first test case is presented. The computation of 3D unsteady mixed convection flows in channels often requires computational domains of long and/or wide aspect ratios, fine space and time discretizations. Therefore, efficient numerical methods are needed to solve the conservation equations. The interest of this first test case is that its computational cost is quite reasonable and it is accessible with limited computational facilities. It is indeed steady and, the computational domain being extended up to the fully-established zone, a homogeneous Neumann boundary condition at the outflow accommodates the problem solution. It is therefore much more easy to compute than the second test case in which a steady state is never established since unsteady thermoconvective rolls are permanently transported by the flow.

In the present work, we solve the first benchmark problem using four different CFD research codes and three discretization methods: finite difference, finite volume and finite element methods. All contributors have mobilized a significant amount of computational resources to achieve reliable spatial convergence with each code. Approximate solutions have been obtained on successively refined grids so that Richardson extrapolation (RE) could be used to extent the results. This technique enables one to improve the accuracy of the discrete solutions when used in the asymptotic range of the numerical methods.

However, in the present test case, difficulties in the use of RE have appeared due to the mixed thermal boundary conditions on the channel bottom and top plates. Indeed, to try to reproduce the operating conditions of the PRB experiments by Pabiou et al. [START_REF] Pabiou | Wavy secondary instability of longitudinal rolls in Rayleigh-Bénard-Poiseuille flows[END_REF], adiabatic Neumann conditions are imposed near inlet while isothermal Dirichlet conditions are imposed downstream. This generates a thermal gradient discontinuity at the boundary condition junction. Because the method used to determine the reference solution is as important as the reference values themselves, the consequences of this singularity will be discussed in this article so that we can evaluate the degree of validity of the proposed reference solutions.

The geometry, the governing equations, the boundary conditions and the flow parameters of the simulated test case are described in section §2. The solvers of the different contributors are briefly presented in section §3 and references are given for more details. The methodology of RE adopted to compute the extrapolated reference solutions of the present benchmark problem is described in section §4. The fundamental assumptions for the validity of the RE technique are recalled in §4.1. The influence of the boundary condition singularity on the convergence order of RE is discussed in §4.2. Technical aspects to compute reference solutions of local values by RE are given in §4. [START_REF] Nicolas | Numerical simulation of thermoconvective flows and more uniform depositions in a cold wall rectangular APCVD reactor[END_REF]. The results are presented and analyzed in section §5. In §5.1, RE of integral values over the whole domain are discussed and used to determine the observed convergence order of our numerical methods when a singularity is present in the domain. Then, in §5.2, the streamwise and spanwise profiles of the velocity, temperature and wall Nusselt number are presented and selected extrapolated extrema on these profiles are given. The profiles of the convergence orders observed from RE of these quantities are also discussed from the viewpoint of the singularity. In §5.3, the extrapolated values of the heat and momentum fluxes through the channel boundaries are computed. Finally, in §5. [START_REF] Wang | Instability and heat transfer in mixed-convection flow in a horizontal duct with discrete heat sources[END_REF], we propose an explanation for the observed behaviors of RE and a correction to the Taylor expansion involved in the extrapolation formalism. The conclusions and the difficulties that have been raised during the study are summarized in section §6. The proposed benchmark is a PRB flow, made of ten steady longitudinal thermoconvective rolls, in the horizontal rectangular channel drawn in Figure 1. A Poiseuille flow is imposed at the channel entrance and the incoming fluid is cold, at temperature T c . After an entrance zone over which a zero heat flux is imposed on the four walls, the top horizontal wall is maintained at the cold temperature T c and the bottom wall is maintained at a higher temperature T h . The vertical lateral walls are adiabatic. A and B are the streamwise and spanwise aspect ratios of the computational domain and A e is the streamwise entrance aspect ratio (see Figure 1). The working fluid is Newtonian and the flow is governed by the 3D incompressible Navier-Stokes equations under the Boussinesq assumption. Using the channel height H, the mean flow velocity U mean , ρU 2 mean and H/U mean as reference quantities for lengths, velocities, pressure and time respectively, and using the reduced temperature θ = (T -T c )/(T h -T c ), the governing equations take the following dimensionless form:

Test case description

         ∇. -→ v = 0 ∂ -→ v ∂t + ( -→ v .∇) -→ v = -∇p + 1 Re ∇ 2-→ v + Ra P r Re 2 θ -→ k ∂θ ∂t + -→ v .∇θ = 1 P r Re ∇ 2 θ (1)
where x, y, z, t, -→ v = (u, v, w) and p are the dimensionless streamwise, spanwise and vertical coordinates, time, velocity vector and pressure, -→ k is the upward unit vector, P r is Prandtl number (= ν/α), Re is Reynolds number (= U mean H/ν) and Ra is Rayleigh number (= gβ(T h -T c )H 3 /(να)). The boundary and initial conditions for u, v, w and θ are: • at y = 0 and B, -→ v = -→ 0 and ∂θ/∂y = 0;

• at z = 0, -→ v = -→ 0
• at x = -A e , u = u P ois (y, z), v = w = 0 and θ = 0, where u P ois (y, z) is given either directly by an approximate solution of the Poisson equation = Re ∂p ∂x , with no-slip boundary conditions at y = 0 and B and at z = 0 and 1, or by the analytical solution of this equation computed in [START_REF] Nicolas | Linear stability of mixed convection flows in horizontal rectangular channels of finite transversal extension heated from below[END_REF] and given in appendix A.

• at x = A-A e , an outflow non-reflective boundary condition is imposed. The choice of this boundary condition was left free in [START_REF] Medale | CALL FOR CONTRIBUTIONS: Towards numerical benchmark solutions for 3D mixed convection flows in rectangular channels heated from below[END_REF]. Note however that the standard Neumann or Orlanski boundary conditions are appropriate for this test case since the flow is convection dominant.

• at t = 0, ∀x ∈ [-A e , A -A e ], u = u P ois (y, z), v = w = 0 and θ = 0.
The present test case is defined by: Re = 50, Ra = 5000, P r = 0.7, A = 50, B = 10 and A e = 2. The resulting flow pattern is the ten longitudinal roll steady flow which is presented in Figure 2. It is obtained by starting from the initial conditions given above and develops after a transient stage which will not be discussed here. It is symmetrical with respect to the median longitudinal vertical plane and can therefore be computed for y ∈ [0, B/2].

Contributors and solver description

Below, the solvers of the four contributors are presented in a few lines and references are given for more details. The numerical parameters for each of the four schemes are given in Table 1. In the table, we also indicate if the symmetry with respect to the median longitudinal vertical plane was used or not, the mesh size, N x × N y × N z , in each space direction, the time step value, ∆t, an estimation of the user time (restitution time) of each computation and the consistency orders1 , α°, of each space discretization method for each primitive variable. Note that, when symmetry is used, N y is the node number on the width B/2 of the computational domain.

Second order finite difference vectorized code: FD1

The test case solution "FD1" is computed using a finite difference method, optimized for vectorial computers. The time discretization scheme is a second-order Adams-Bashforth scheme, combining an explicit treatment of the convective term and an implicit treatment of the diffusive term. The equations are discretized in space on uniform, Cartesian and staggered grids using centered differences for the diffusive terms and a second order upwind or central scheme for the convective terms. With the two finest meshes used in this study (see Table 1), to avoid numerical scheme instabilities and save CPU time, the solution is computed first with the upwind scheme and ∆t = 0.01. Then, starting from this converged steady solution as initial condition, a new solution is computed with the central scheme and ∆t = 0.002. The time integration and the velocity-pressure coupling are computed by the projection method based on Goda's algorithm [START_REF] Goda | A multistep technique with implicit difference schemes for calculating twoor three-dimensional cavity flows[END_REF]. The Helmholtz equations for the temperature field and the components of the predicted velocity field are solved using an incremental factorization method of ADI type which preserves a second order time accuracy. The Poisson equation for the pressure increment is solved by a factorization method which consists in the partial diagonalizing of the mono-dimensional Laplace operators in the transverse directions y and z. The linear systems resulting from these two factorization methods are all tridiagonal and are solved by the TDMA algorithm. An Orlanski type boundary condition is used at Since this solver is highly vectorizable, the code is very efficient on vectorial supercomputers: for instance, it runs at 12.5 Giga Flops on average on the NEC-SX8 computer at IDRIS (the CNRS supercomputing center at Orsay, France), when the peak power of this computer is 16 Giga Flops. A detailed description of this code and of its validations and performances can be found in Benzaoui et al. [START_REF] Benzaoui | Efficient vectorized finite difference method to solve the incompressible Navier-Stokes equations for 3D mixed convection flows in high aspect ratio channels[END_REF].

Third order finite element parallelized code: FE2

The numerical model "FE2" is based on a segregated approach to build up separate integral forms associated with the set of coupled governing equations [START_REF] Nicolas | Benchmark solution for a threedimensional mixed convection flow -Part 1: reference solutions[END_REF]. The fluid flow problem is kept in primitive variable formulation and solved using an unconditionally stable projection algorithm [START_REF] Guermond | A new class of truly consistent splitting schemes for incompressible flows[END_REF]. As in most projection type algorithms the incompressibility constraint is enforced in the FE2 code through a pressure correction field computed from a pressure Poisson equation. The latter is obtained by taking the divergence of the momentum equation in equations [START_REF] Nicolas | Benchmark solution for a threedimensional mixed convection flow -Part 1: reference solutions[END_REF] and Neumann boundary conditions. Non homogeneous Neumann boundary conditions have been implemented for the pressure Poisson equation in a form derived from [START_REF] Orszag | Boundary conditions for incompressible flows[END_REF][START_REF] Karniadakis | High-order splitting methods for the incompressible Navier-Stokes equations[END_REF]:

∂p ∂ n = (- ∂ v ∂t -v∇ v - 1 Re ∇ × (∇ × v) + Ra P rRe 2 θ k) • n (2)
The mechanical stress and heat flux outlet boundary conditions arising at x = A for the momentum and energy equations have been treated with a formulation inspired from [START_REF] Papanastasiou | A new outflow boundary condition[END_REF] and adapted to the present framework combining a segregated approach for the mixed convection problem together with a projection algorithm. The spatial discretization of the three separate integral forms, associated with temperature, velocity and pressure unknowns, follows the standard finite element method, using tri-quadratic hexahedral finite elements for the velocity and temperature fields and tri-linear approximation for the pressure field. The non-linear algebraic system resulting from the discretization of the momentum equation is solved using a Newton-Raphson procedure, despite only partial convergence is required for solving this nonlinear system during the transient solution in the segregated procedure. The time integration is performed with a second order Backward Difference Formula scheme (BDF2) [START_REF] Guermond | A new class of truly consistent splitting schemes for incompressible flows[END_REF]. At each time step the three algebraic systems corresponding to the momentum, incompressible projection and energy conservation are solved with an iterative solver (Bi-Conjugate Gradient Stabilized, pre-conditioned with Additive Schwartz Method) provided in the PETSc toolkit [START_REF] Balay | Petsc Web pages[END_REF]. This implementation enables us to efficiently run high performance massively parallel computers (IBM SP4 and SP6 at IDRIS). Finally, the computational domain is discretized with three meshes uniformly spaced in the x, y and z directions and whose node numbers are given in Table 1, e.g., the finest mesh consists of 675 × 135 × 45 tri-quadratic hexahedra finite elements, built on 1351 × 271 × 91 nodes.

Second order finite volume parallelized code: FV3

The test case solution "FV3" is computed using the finite volume code Thétis developed at the I2M Institute. Time discretization of the Navier-Stokes and energy equations is implicit. A first order Euler scheme is used, with an implicit treatment of all the terms of the equations (after linearization of the inertial term of the Navier-Stokes equations and after uncoupling with the energy equation). The incompressibility constraint that couples the velocity and the pressure is solved using a pressure correction scheme [START_REF] Timmermans | An approximate projection scheme for incompressible flow using spectral elements[END_REF]. It consists of splitting the Navier-Stokes system into two stages, a velocity prediction stage and a pressure correction stage. The spatial discretization is based on the finite volume method on a velocity-pressure staggered grid of the Marker and Cells type. Pressure unknowns are located at the cell vertices whereas velocity components are face centered. A centered scheme of order 2 is used for the inertial and stress terms of the Navier-Stokes equations, as well as for the pressure correction step and the diffusive term of the energy equation. The convective term of the latter equation is discretized with the Quick scheme to avoid numerical instabilities [START_REF] Leonard | A stable and accurate convection modelling procedure based on quadratic upstream interpolation[END_REF]. A Neumann boundary condition is set on the outlet boundary for velocity and temperature. The code is parallelized in a distributed way and runs efficiently on hundreds of processors using the parallel solver library Hypre [START_REF] Falgout | The Design and Implementation of Hypre, a Library of Parallel High Performance Preconditioners, chapter in Numerical Solution of Partial Differential Equations on Parallel Computers[END_REF]. Among the different solvers and preconditioners available in this library, the most efficient for this problem are a GMRES solver for the prediction step and the energy equation with a point Jacobi preconditioner and a BiCGStab solver with a multigrid preconditioner for the correction step. Three meshes were used. The first one begins with 600 cells in direction x to avoid small oscillations observed in the temperature field with coarser meshes. Simulations are stopped when the stationary criteria of 10 -10 is reached and L 2 norm of the divergence is below 10 -10 . An SGI Altix 8200 cluster was used composed of 152 quadcore dual Intel Xeon processor blades.

Third order finite element stationary parallelized code: FE4

The numerical model "FE4" spatially discretizes equations (1) in stationary form (without the temporal derivative term) using an LBB-stable [START_REF] Ern | Eléments finis: théorie, applications, mise en oeuvre[END_REF] finite element method. No upwinding of the convective term is used. To find the solution of the discrete nonlinear stationary problem, we used a defect-correction solver [START_REF] Turek | Efficent solvers for Incompressible Flow Problems -An Algorithmic and Computational Approach[END_REF]. A certain linearization of the equations and additional regularization terms give an inexact tangent matrix. This tangent matrix is then approximately factorized as in the algebraic projection method [START_REF] Quarteroni | Factorization methods for the numerical approximation of Navier-Stokes equations[END_REF]. This leads to segregation of the linear systems to be solved for each scalar incremental unknown (3 velocities, 1 pressure, 1 temperature). The mesh is regular and consists in hexahedral elements, triquadratic (Q 2 ) for the velocity and temperature unknowns and linear discontinuous (P nc 1 ) for the pressure. Thus, the formal spatial discretization order of the method is 3 for the velocity and temperature unknowns and 2 for the pressure unknown. The total number of degrees of freedom is 16 million for the coarsest mesh and 73.3 millions for the finest one. The standard natural boundary condition on momentum for the discretization used [START_REF] Ern | Eléments finis: théorie, applications, mise en oeuvre[END_REF] is µ∂u n /∂n -P = 0. A boundary term in -P last (the last pressure estimation) is discretized and added to the right hand side of the boundary condition so that we get the desired µ∂u n /∂n = 0 when convergence is reached. The inexact tangent matrix is obtained from the following contributions: exact tangent matrix for the diffusion, pressure gradient and velocity divergence terms, fixed point linearization for the convective terms and a regularizing pseudo-time like mass term on the velocity and temperature. The linear systems are solved with BiCGSTAB preconditioned by an ILU(0) incomplete factorization [START_REF] Saad | Iterative methods for sparse linear systems[END_REF] for the velocity and temperature unknowns and FCG(1) preconditioned by algebraic multigrid for the pressure unknown. We used the algebraic multigrid method of Notay [START_REF] Notay | An aggregation-based algebraic multigrid method[END_REF] in sequential mode. To speed-up convergence towards the final steady state, a four-point acceleration method is used. The numerical model was implemented in Cast3m [START_REF]Cast3m web site[END_REF]: a freely available finite element code developed at CEA (French Atomic Energy Commission). The model was run on standard PC servers running Linux with up to 8 cores and 64 GB RAM. The most CPU intensive part of the model is the pressure linear system solves.

Richardson extrapolation method

General principle without singularity

When the approximate solutions of a continuous initial and boundary value problem are computed by discretization methods such as finite difference, finite volume or finite element methods, RE can be used to improve the precision of the discrete solutions. Indeed, provided that three main assumptions are satisfied (see below), it is possible to get an order of accuracy of at least O(h p+1 ) when the convergence order of the numerical method is O(h p ), where h is the mesh size. This technique then allows one to compute extrapolated primitive variables at any point of the computational domain as well as solution functionals such as differentiated or integrated quantities (heat and momentum fluxes, volume or surface averaged quantities, and so on). A concise and elegant presentation of RE to estimate a posteriori discretization errors in computational simulations can be found in [START_REF] Alvin | Methodology for characterizing modeling and discretization uncertainties in computational simulation[END_REF]. More details and deeper discussions on the theory are given in [START_REF] Roache | Perspective: a method for uniform reporting of grid refinement studies[END_REF][START_REF] Oberkampf | Verification and validation in computational fluid dynamics[END_REF][START_REF] Roy | Review of code and solution verification procedures for computational simulation[END_REF]. Examples of extrapolated solutions in natural and mixed convection problems can be found in [START_REF] Evans | The thermoconvective instability of plane Poiseuille flow heated from below: a proposed benchmark solution for open boundary flows[END_REF][START_REF] Le Quéré | Modelling of natural convection flows with large temperature differences: a benchmark problem for low Mach number solvers. Part 1. Reference solutions[END_REF].

RE first consists of computing the numerical solutions f h i (1 ≤ i ≤ N) of the discretized boundary value problem on N different nested uniform grids of size h i , with h 1 the coarsest grid and h N the finest one. If (assumption {A1}) the exact solution of the continuous problem, f exact , is sufficiently smooth to justify the use of Taylor expansion (at least up to the discretization order), then it can be written in the form:

f h i = f exact + C α h α i + O(h α+1 i ) (3)
where C α is a coefficient which is dependent on α but independent of h i . Then, the leading order α of the truncation error due to discretization, the coefficient C α and the exact solution f exact can be approximated from the discrete solutions, if two more assumptions are satisfied. The second assumption {A2} is that the mesh spacings h i used in the extrapolation must be small enough so that the discrete solutions f h i are located in the asymptotic convergence region that is the leading order term C α h α i of the truncation error truly dominates the total discretization error f exact -f h i . In this case, α will be considered as the observed convergence order from RE.

Thus, using three grids (N = 3) such as h 1 h 2 = h 2 h 3 , the approximations α, Cα and f ex of α, C α and f exact in equation ( 3) are given by [START_REF] Oberkampf | Verification and validation in computational fluid dynamics[END_REF][START_REF] Roy | Review of code and solution verification procedures for computational simulation[END_REF]:

α = ln f h 1 -f h 2 f h 2 -f h 3 ln h 1 h 2 Cα = f h 2 -f h 3 h α 2 -h α 3 (4)
and, using four grids (N = 4) such as h 1 h 2 = h 3 h 4 , they are given by:

α = ln f h 1 -f h 3 f h 2 -f h 4 ln h 1 h 2 Cα = f h 3 -f h 4 h α 3 -h α 4 (5)
f ex = f h 4 -Cα h α 4 with Cα = C α + O(h N -1 ) and f ex = f exact + O(h α+1 N ).
As a consequence, the approximation f ex of the asymptotic solution f exact will be better if h N is small and α is large. Thereafter, α and f ex will respectively be noted α and f ex .

The formal expression of the Taylor expansion ( 3) is valid for multidimensional problems, in any coordinates, including space and time, only if (assumption {A3}) the same grid refinement ratio is applied in all space and time directions. In our stationary problem, this means that the cell aspect ratios are kept constant from one grid to another. That is, if N uniform Cartesian grids of size ∆x i , ∆y i and ∆z i (i = 1, ..., N) are used for RE, the ratios ∆x i ∆z i must be equal whatever i, and the same holds for ∆y i ∆z i [START_REF] Roache | Perspective: a method for uniform reporting of grid refinement studies[END_REF][START_REF] Oberkampf | Verification and validation in computational fluid dynamics[END_REF][START_REF] Evans | The thermoconvective instability of plane Poiseuille flow heated from below: a proposed benchmark solution for open boundary flows[END_REF]. On smooth problems, the spatial convergence orders of the codes used to compute the FD1 and FV3 solutions have been shown to be equal to two (see [START_REF] Benzaoui | Efficient vectorized finite difference method to solve the incompressible Navier-Stokes equations for 3D mixed convection flows in high aspect ratio channels[END_REF] for the FD1 code) while, for the two finite element codes FE2 and FE4, they have been shown to be equal to three, for the temperature θ and the velocity components u, v, w. That is, for the four codes used in the present study, the spatial convergence order for u, v, w and θ is equal to the consistency order, α°, mentioned in Table 1. As a consequence, if the solution f exact of the problem is smooth enough (say f exact ∈ C 2 (Ω ∪ ∂Ω), where Ω ∈ R 3 is the computational domain and ∂Ω ∈ R 2 its boundary), the u, v, w and θ values that will be given as reference solutions from RE should only be obtained with an associated extrapolation coefficient α equal to two for the FD1 and FV3 solutions and to three for the FE2 and FE4 solutions. Otherwise, if α is very different from α°, this means that the discrete solutions used to compute the extrapolated solution are not in the asymptotic convergence region of RE (assumption {A2} is not satisfied) or that the solution of the continuous problem defined in §2 is not smooth enough (assumption {A1} is not satisfied). This is precisely what is observed in the present problem and what is developed in the following subsection.

Singularity influence on the RE convergence orders

The solution of the present problem is not smooth because the temperature gradient on the horizontal plates is discontinuous at x = 0. This is due to the use of homogeneous Neumann (adiabatic) conditions for -2 ≤ x < 0 and Dirichlet boundary conditions for x ≥ 02 . To understand the effect of this singularity on the convergence rate of the solutions and, as a consequence, on the convergence rates, α, associated with RE, let us refer to the finite element framework. Without any singularity, the theoretical finite element discretization error of an elliptic problem is given by the basic interpolation theory to be O(h α°) in the L 2 norm, where the consistency order of the discretization method is given by α°=k + 1, with k the polynomial approximation degree and h the characteristic mesh size (k=2 and α°=3 for θ, u, v and w in the FE2 and FE4 methods). But, as soon as a singularity enters the problem, the rate of convergence of the numerical model becomes O(h min(α°,r) ), where r measures the problem regularity influence on the actual convergence rate (see [START_REF] Zienkiewicz | The Finite Element Method[END_REF], section 14, p. 404). It seems therefore reasonable to consider that the convergence rates, α, obtained from RE of integral quantities based on a norm should be equal to min(α°, r), if the grids are located in the asymptotic convergence regions of the discrete solutions. However, what happens for RE of other quantities, such as local extrema for instance? Does α vary between min(α°, r) and max(α°, r)? Does the singularity only influence its neighborhood or the whole domain? These issues will be dealt with in §5.

Strang and Fix, in chapter 8 "Singularities" of reference [START_REF] Strang | An analysis of the Finite Elements Method[END_REF], propose a theoretical analysis to a priori determine the regularity r of a singular boundary value problem due to the discontinuity of its boundary conditions. They analyze the singularity influence on the convergence rate of finite element methods. Their analysis focuses on a test case corresponding to the displacement computation in a 2D domain with a crack. As shown in Figure 3, since the crack axis is a symmetrical axis, this test case can be viewed as a 2D pure diffusion problem (Poisson problem) with a mixed Dirichlet/Neumann condition on one of its boundaries. That is this test case presents the same singularity as in the present benchmark problem except it is 2D instead of being 3D and there is no convection. Strang and Fix [START_REF] Strang | An analysis of the Finite Elements Method[END_REF] show that the solution around the singularity behaves like x 1/2 , where x is here the distance to the singularity, and that the error of the finite element solutions of this elliptic problem is O(h) for any choice of element. Thus, the problem regularity is r = 1. But they also mentioned that away from the singularity a higher regularity can be expected due to a decrease of the singularity pollution. In the present benchmark exercise, the determination of the problem regularity will be done a posteriori in §5.1.2 by performing specific numerical spatial convergence studies. We will show, as in Strang and Fix [START_REF] Strang | An analysis of the Finite Elements Method[END_REF], that the regularity r of the solutions attached to the thermal field tends to 1 near the singularity and increases up to α°far away.

Boundary conditions : du/dn=0 A C A B C B D E G F D E 1 2 u=0
Figure 3: On the left, the square domain ( 1) with a crack is used by Strang and Fix [START_REF] Strang | An analysis of the Finite Elements Method[END_REF] to study the effect of the singularity at point B. On the right, the equivalent domain ( 2) is obtained using the symmetry of domain (1) through the line ABC. The domain (2) has the same boundary condition singularity at point B as the present benchmark problem for the temperature at x = 0 and z = 0 or 1 (see Figure 1).

Technical aspects of Richardson extrapolation of local values

Grid doubling or integer grid refinement ratios are not required for the validity of RE. Thus, in the general case, solutions of the coarsest grids are not computed at the nodes of the finest grid. However, if local values have to be evaluated by RE, it is necessary to have the values of f h i at the same points in equations ( 4) and [START_REF] Nicolas | Bibliographical review on the Poiseuille-Rayleigh-Bénard flows: the mixed convection flows in horizontal rectangular ducts heated from below[END_REF]. As suggested by Roache [START_REF] Roache | Perspective: a method for uniform reporting of grid refinement studies[END_REF], this is made possible by building an interpolation of the coarse solutions on the finest grid, the order of which is higher than the space discretization order of the used numerical method.

Since the consistency order, α°, of the numerical methods used in the present paper is equal to two or three (see Table 1), a cubic spline interpolation has been used to compute the solutions of the coarsest grids at the nodes of the finest one, before doing RE of the temperature, Nusselt number and velocity local extrema (see §5.2.3). Indeed, cubic spline interpolation is supposed to be third order if the solution is smooth enough and the points far enough from the boundaries. Other interpolation methods could have been used: for instance, an evaluation of the influence of Lagrangian interpolations on RE is proposed in [START_REF] Shyy | Evaluation of Richardson extrapolation in computational fluid dynamics[END_REF].

In the present paper, the values and the coordinates of the local extrema of the extrapolated solution are also computed using cubic spline interpolation between the finest grid nodes. This is illustrated in Figure 4 in which a zoom in the w streamwise profiles computed with the four grids of FD1, together with the extrapolated profiles, are presented. This figure allows to determine the FD1 value and the streamwise coordinate of the vertical velocity local extremum noted w 2 and x 2 in Table 3 1, together with the extrapolated solution. The latter can be considered as the asymptotic solution here because the α values monotonously varies between 2.01 and 2.16 when x varies between 24.7 and 25.3. In this figure, all the symbols correspond to the computational or extrapolation points and the curves to the cubic spline interpolation curves. The small window zooms in the local extremum of the extrapolated curve. The value w 2 = -0.472989 and the coordinate x 2 = 24.907 of this extremum are reported in Table 3.

Result presentation and analysis

Richardson extrapolation of integral quantities

Presentation and discussion of the results

The extrapolated values of integral quantities are computed first. These integrals are twice the mean kinetic energy, 2E c , on the whole domain of volume D, the mean pressure difference, ∆P io , between inlet and outlet and the mean temperature, T m , on the whole domain that are defined by:

2E c = 1 D D u 2 + v 2 + w 2 dx dy dz (6)
∆P io = 1 S i S i P dy dz - 1 S o So P dy dz (7)
T m = 1 D D θ dx dy dz (8)
They have been computed using either the middle point rule for the FD1 and FV3 solutions or 3 × 3 × 3 Gauss integration scheme for the FE2 and FE4 solutions. The advantage of these integrals is that they only depend on the primitive variables: no differentiation and no interpolation are needed to compute their values on each grid. Thus, in this case, the validity of RE and the values of the associated extrapolation coefficient, α Ec , α ∆P io and α Tm , only depend on the convergence orders of the numerical methods and on the satisfaction of the three assumptions {A1} to {A3} stated in §4.1.

The values of f = (2E c , ∆P io , T m ) on the finest grid (noted f f g ) and extrapolated by RE (noted f ex ) are given in Table 2 with the associated convergence order, α f , and with the relative distance between the finest grid and extrapolated solutions:

d f = (f f g - f ex )/|f f g |.
To check if the discrete solutions FD1 and FE4 are located in the asymptotic convergence region, the first and fourth contributors have also computed Log|f h i -f ex | = Log(h i ) on each of their grids and have compared the slopes of the linear interpolations of the values of this functions with the α f values. For FD1 solutions, each of these slopes is strictly equal to the associated α f value with four common figures. Referring to [START_REF] Oberkampf | Verification and validation in computational fluid dynamics[END_REF][START_REF] Roy | Review of code and solution verification procedures for computational simulation[END_REF], this proves that the grid asymptotic convergence region is reached by the FD1 contribution. On the other hand, for FE4 solution, the slopes of the linear interpolations and the α f values are nearly equal for T m , but a little bit different and smaller than 1 for ∆P io , and very different and negative for 2E c (see Table 2). Consequently, the asymptotic convergence region is not reached for ∆P io and 2E c and therefore RE can no longer be used for these FE4 quantities. Note that such a reasoning cannot be hold for the FE2 and FV3 solutions because RE is based on three grids only: in this case, the three points Log|f h i -f ex | = Log(h i ) are by construction automatically aligned whatever the f h i values and the slope is necessarily equal to α f . However, we note that 2E f g c and ∆P f g io values on the FE4 finest grid are very close to the extrapolated values 2E ex c and ∆P ex io of the three other contributors. Thus, for each quantity of Table 2, a reference value with its uncertainty margin and the precision of the reference value determination have been defined. The reference value, f ref , is equal to the arithmetic average of the extreme values of the FD1, FE2 and FV3 extrapolated values, plus the extrapolated or the finest grid FE4 value depending on whether RE succeeds or not. The uncertainty margin, f margin , is equal to the half difference between the two extreme values. The precision of the determination, f prec , is equal to the ratio of the uncertainty margin to the reference value. The reference values in Table 2 are given with a precision of the order of 10 -5 , with five common figures among the four extrapolated solutions for 2E c and ∆P io and four common figures for T m .

FD1 FE2 FV3 FE4 References f ref ±f margin f prec = fmargin f ref 2E f g c 2E ex c α Ec d Ec 1.
d f = f f g -f ex
|f f g | , between the extrapolated and finest grid values. FE4 column: the slopes of the linear interpolation of the functions Log|f h i -f h→0 | = Log(h i ) are noted into brackets in italic; the symbol °indicates an erroneous value due to the extrapolation failure (thus the FE4 finner grid value replaces the FE4 extrapolated value in the reference value determination). Right column: reference solutions with their tolerance margin and the precision of their determination.

It can also be noted that the relative distances between the finest grid solutions and the extrapolated solutions (when admissible) are higher for the second order methods (FD1 and FV3) than for the third order methods (FE2 and FE4). All these relative distances are higher than the precision of the determination of the reference values (except for 2E c in FE2 column), showing that RE improves the accuracy of the reference values, more substantially for the second order methods. Furthermore, for the FE4 contribution the maximum relative distance on 2E c (resp. ∆P io ) between the coarsest and finest grids are equal to 4.25 × 10 -6 (resp. 1.45 × 10 -5 ), which is smaller than the precision of the reference values given in Table 2.

It is interesting to analyze the Table 2 results going into details. Indeed, one can see that the coefficients α Ec and α ∆P io obtained when RE is admissible are respectively nearly equal to the consistency orders α u,v,w °and α p °of each numerical method (see Table 1). This corresponds to the expected behavior for a smooth problem without discontinuity (see §4.1). On the other hand, α Tm is nearly equal to one for the four contributions despite α θ °=2 for the FD1 and FV3 solutions and α θ °=3 for the FE2 and FE4 solutions. In §5.1.2 below, we are going to show that the lower than expected α Tm values are due to the singularity induced by the thermal boundary conditions.

Determination of the singularity influence

To analyze the influence zone of the singularity of the present benchmark problem, extrapolated values of the mean temperature T m (see equation [START_REF] Benzaoui | Efficient vectorized finite difference method to solve the incompressible Navier-Stokes equations for 3D mixed convection flows in high aspect ratio channels[END_REF]) and of the L 2 norms of temperature T L 2 and vertical velocity component W L 2 have been obtained by RE in several subdomains of the whole computational domain. Here the L 2 norm is defined by f L 2 = ( D f 2 dx dy dz) 1/2 . In all the subdomains considered, the extension is maximum in y and z transverse directions (that is 0 ≤ y ≤ B/2 or B depending on the contributor and 0 ≤ z ≤ 1) and it varies from x = -2 to x = x°for the upstream subdomains and from x = x°to x = A = 48 for the downstream subdomains, with -2 < x°<48. Then the convergence orders α Tm , α T L 2 and α W L 2 , observed from RE of T m , T L 2 and W L 2 respectively, were computed in the downstream and upstream subdomains. They are drawn as a function of x°in Figures 5(a-c). Note that, for the FD1 and FE4 solutions, the mean temperature computed in each subdomain is such that the points (Log(h i ), Log|T m h i -T ex m |) are perfectly aligned for the four grids considered: these solutions are therefore located in the asymptotic convergence region (see §5.1.1). However, this is not the case for T L 2 and W L 2 at the coordinates x°where α T L 2 and α W L 2 diverge in Figures 5(b,c).

Figure 5(a) shows that RE of T m succeeds except for the FD1 solution in the upstream subdomains near the entrance of the channel (divergence of α Tm ) and for the FV3 solution (negative values of α Tm for the downstream subdomains). On the other hand, Figure 5(b) shows that RE of T L 2 succeeds for the FV3 solution but fails for the FD1, FE2 and FE4 solutions in the downstream subdomains (divergence of α T L 2 ). These different behaviors are maybe due to the different discretization schemes used for the convective terms by the contributors: centered difference schemes for the FD1, FE2 and FE4 solutions and Quick scheme for the FV3 solution. If we only consider the cases with a successful RE, it appears that α Tm → 1 for the upstream subdomains located near the beginning of the heated plate (x°= 0), and α Tm increases when the subdomains are more and more located downstream: for the FD1 solution, the α Tm values on the downstream subdomains stabilize around 1.8 (that is around the consistency order α°= 2) for 20 ≤ x°≤ 48, while it increases beyond α°=3 for the two finite element solutions. The behavior of α T L 2 for the FV3 solution is similar to the behavior of α Tm for the FD1 solution: α T L 2 is nearly equal to one near x°= 0 for the upstream subdomains (α T L 2 ≈ 1.3 for x°=0 ) and it increases towards α°in the downstream subdomains (α T L 2 ≈ 1.7 for x°>20 ). It can therefore be considered that the influence zone of the thermal boundary condition singularity on the convergence orders of the quantities linked to temperature extends to x = 20 to 25 in the present PRB flow, regardless of the numerical method used. This zone precisely corresponds to the development zone of the longitudinal rolls as it can be seen in Figure 2(a). As a conclusion of this analysis, it appears that the effect of the singularity of the present problem spatially varies: it diminishes with the distance to the singularity. Furthermore the exponent r introduced in §4.2 to characterize the regularity of the solution tends to 1, regardless of the numerical method used, at least in the neighborhood of the upstream edge of the heated zone and for the quantities derived from the temperature field.

The dynamics fields should also be affected by this singularity through the buoyancy term in the momentum equation which makes a coupling of the velocity fields with the temperature field. To study this influence, RE of W L 2 in the different subdomains is presented in Figure 5(c). RE of W L 2 succeeds for the downstream subdomains (except for the FE2 solution for x°< 20) but fails for the FD1 and FV3 solutions in the upstream subdomains. It appears that α W L 2 ≈ 2 = α°for the FD1 and FV3 solutions in the downstream subdomains whatever the x°value, whereas α W L 2 ≈ 1.4 near inlet for the FE4 solution and α W L 2 tends to values that vary between 3 and 5 for x°>20 for the two FE solutions. Thus, the singularity of the thermal boundary conditions does not seem to affect the velocity field with the second order FD1 and FV3 methods for the used grids, whereas it influences the velocity field with the third order FE2 and FE4 methods. In this last case, the length of the influence zone of the singularity is the same as for the temperature integrals in Figures 5(a,b): it reaches x°≈25.

As a consequence, in the following, we have considered that the reference quantities defined in the present benchmark problem should be established from the extrapolated quantities only if the associated convergence rates α are such as 1 ≤ α ≤ α°. In practice, we take into account numerical errors by using superior tolerance margins on α to choose the conserved extrapolated values: in this paper, the reference solutions are defined from the extrapolated solutions with 1 ≤ α 2.5 for the FD1 and FV3 contributions and with 1 ≤ α 4 for the FE2 and FE4 contributions.

Richardson extrapolation of temperature, velocity and Nusselt number local extrema

Space profiles of the thermal and dynamical fields

In the following, we denote by Nu t and Nu b the local Nusselt numbers on the top and bottom walls respectively. They are defined by:

Nu t,b (x, y) = - H ∂T ∂Z Z=H,Z=0 T h -T c = - ∂θ ∂z z=1,z=0 (9)
In the variational context of finite element methods, it is possible to compute the Nusselt numbers in several ways. The "non consistent" way simply uses the definition (9), i.e. the z-derivative of the interpolation function for θ is computed. The "consistent" way exploits the duality between Dirichlet and Neumann boundary condition. The "consistent" flux at a Dirichlet boundary node is the one that would yield the same solution if prescribed instead of the Dirichlet condition. Details on how to compute such a flux in a finite element framework are given in references [START_REF] Gresho | The consistent Galerkin FEM for computing derived boundary quantities in thermal and/or fluids problems[END_REF][START_REF] Hughes | The Continuous Galerkin Method is locally conservative[END_REF]. A reported advantage of the "consistent" flux is that it is generally more precise than the non-consistent one. This is also what we have observed here (see §5. In Figure 6, the longitudinal profiles of the primitive variables θ, u, v and w are drawn along the lines (y, z) = (2, 0.2) and (5, 0.5) and the profiles of Nu t and Nu b are drawn along the lines y = 2 and y = 5. One can note that the longitudinal roll flow computed here is not fully established in space since all the longitudinal profiles go on slightly evolving in the streamwise direction at x = 48. The transverse profiles of θ, u, v and w are drawn in Figure 7 along the four lines at x = 10 and 30 and at z = 0.2 and 0.5. The transverse profiles of Nu t and Nu b are drawn in Figure 8 along the lines at x = 10 and x = 30. Only the first half of these transverse profiles is shown because the flow is symmetrical with respect to the median vertical plane (y = 5). The transverse profiles at x = 10 are located in the entrance region, more precisely at mid-length of the forced convection triangular zone, where only two longitudinal rolls are present along each vertical wall (see Figure 2(a)). On the other hand, the transverse profiles at x = 30 are located where ten well developed longitudinal rolls are present. 1, together with the extrapolated solution. The latter does not tend to the asymptotic solution near x = 22.35 because the α values (the black filled circles in the figure) diverge where the profiles intersect. All the symbols correspond to the computational or extrapolated points. The curves correspond to the cubic spline interpolation curves, except for the extrapolated solution where they are linear segments linking the extrapolated points.

All these profiles are computed from the FD1 solution on the finest mesh. The same profiles are obtained with the other numerical methods (FE2, FV3, FE4) if the comparisons are done at the same scales as those of Figures 6, 7 and 8. Note that it is not possible to draw the extrapolated profiles in the present problem because it is impossible to be located in the asymptotic convergence region along the whole profiles. In particular, RE diverges at points where the profiles computed on two distinct meshes intersect. Indeed, when f h i = f h j for h i = h j , α diverges in equations ( 4) or [START_REF] Nicolas | Bibliographical review on the Poiseuille-Rayleigh-Bénard flows: the mixed convection flows in horizontal rectangular ducts heated from below[END_REF]. This is illustrated in Figure 9 that focuses on a zone where a curve crossing is present. Such behaviors are also observed in [START_REF] Garbey | A least square extrapolation method for improving solution accuracy of PDE computations[END_REF].

Space profiles of the observed convergence rates from RE

To complete the preceding observations, a selection of streamwise and spanwise profiles of the convergence rates, α, observed from RE of the different primitive variables and Nusselt numbers are presented in Figures 10 and11. It can be shown that the α profiles are very erratic and that RE can even fail. This happens when the values of the studied quantity do not monotonously vary from one grid to the following finer one. This behavior is indicated by arbitrarily fixing α to zero in some profiles. As already seen in Figure 9, α profiles present several sharp overshoots and undershoots at the points where the field profiles on the different grids intersect [START_REF] Garbey | A least square extrapolation method for improving solution accuracy of PDE computations[END_REF]. This is the case for instance at x = 0 for all the variables of the four contributions, but also in nearly all the entrance zone for the FD1, FE2 and FE4 contributions. This is due to the probable conjunction of two causes. First the exact solution of the cold Poiseuille flow imposed as inlet boundary condition at x = -2 is nearly conserved untill the beginning of the heated plate at x = 0 on all the grids. Second, the FD1, FE2 and FE4 contributions use centered discretization schemes for the convective terms and very small oscillations are observed in their temperature and velocity streamwise profiles around x = 0 with their coarsest grids such as N x ≤ 601 (more precisely, no velocity oscillation is observed in the FD1 solution and very small velocity oscillations are observed on all the grids of the FE4 solution). These oscillations generally appear just around x = 0 because a streamwise acceleration of the flow due to the density variation near the bottom plate and high transverse thermal gradients are present at the same place. No oscillations are observed in the FV3 solutions because the Quick scheme is used (see §3.3).

In a general way, the FD1 and FV3 α profiles are much more regular than the FE2 and FE4 ones. The FD1 and FV3 α values for θ, u, v, w are nearly equal to α°=2 in the downstream zone for x > 20. On the other hand, the FE4 and FE2 RE can fail, even in the downstream zone, or can succeed but with associated α values very different from α°=3. This is very likely due to the higher precision of the finite element methods used here. Indeed, the values computed on each grid with these methods are very near one from the others: for instance, the maximum relative distance on the primitive variables The corresponding profiles of θ f g , u f g , v f g and w f g on the finest grid are also drawn on the figures.

between the coarsest and the finest grid solutions of the FE2 and FE4 contributions generally varies between 10 -4 and 10 -5 (or even less) when it varies between 10 -2 and 10 -3 for the FD1 and FV3 solutions. As a consequence, the finite element solutions are very sensitive to the numerical errors, to the entrance singularity and to the curve crossings. A way of limiting these negative effects on the RE with the high order methods would have been to increase the size ratio of the successive grids. However this solution has appeared impossible in the present case due to the computational costs on grids much finer than those already used. In Figures 10 and11, it can also be noted that the α profiles computed by the two finite element codes are very similar and that they diverge nearly at the same points. Furthermore, for the four contributions, the α values associated with θ and Nu t are generally smaller than α°for x < 20, except where over and undershoots are present. They tend to 1 for x < 10 due to the singularity influence (cf. §5.1.2). The α values FE4 References θ,u,w) local extrema along the line (y, z) = (5, 0.5) and of their streamwise coordinates, x; truncation error leading order, α, of RE and relative distance, d = f f g -f ex |f f g | , between the extrapolated and the finest grid values (when the extrapolated value is obtained). In the References column, the reference value and the margin on the primitive variables and on their coordinates are given, as well as the precision on the primitive variables.

f ref ± f margin x ref ± x margin f prec = fmargin f ref θ 1 x 1 α 1 d 1 0.
As it has just been discussed, extrapolated solutions of the present problem cannot be determined for the whole field but only for some local values, such as local extrema, or for integral quantities (see §5.1.1 and §5.3 for instance). The extrapolated values and the coordinates of forty local extrema, identified by small filled circles on the profiles of Figures 678along the lines at (y, z) = (5, 0.5) and (2, 0.2) and at (x, z) = (30, 0.5), have been computed using the method described in §4.3. The local extrema of the primitive variables, θ, u, v and w, are given in the Tables 3-5 and those of the Nusselt numbers Nu t and Nu b in the Tables 678.

FD1 FE2 FV3 FE4
References local extrema along the line (x, z) = (30, 0.5) and of their spanwise coordinates, y, and truncation error leading order, α, of RE. See the Table 3 caption for the description of the References column.

f ref ± f margin x ref ± x margin f prec = fmargin f ref θ 1 x 1 α 1 d 1 0.

FD1 FE2 FV3

FE4 References

f ref ± f margin x ref ± x margin f prec = fmargin f ref N u t1 x 1 α 1 d 1 N u cons t1
x local extrema along the line y = 5 and of their streamwise coordinates, x; truncation error leading order, α, of RE and relative distance, d = f f g -f ex |f f g | , between the extrapolated and finest grid values (when the extrapolated value is got). In the reference column, the reference value and the margin on the Nusselt numbers and on their coordinates are given, as well as the precision on the Nusselt numbers. For the FD1, FE2 and FV3 solutions, more than 70% of the whole extrema have been extrapolated with an associated extrapolation coefficient, α, whose value is equal to α = α°±20%, that is very close to the spatial consistency order, α°, of the numerical method used. For the other extrema, α values do not agree with the consistency order for the various reasons already listed above: intersection of the profiles computed on the different grids, influence zone of the thermal boundary condition singularity for x < 20 and, probably, influence of the boundary conditions on the cubic spline interpolations when the extrema are very close to the wall. For the FE4 solutions, the α values are very different of α°for the majority of the extrema. We consider that FE4 RE fails because its second assumption {A2} is not satisfied (the discrete solutions are not located in the asymptotic convergence region). Other explanations are proposed in §5.4.

FD1 FE2 FV3 FE4
References

f ref ± f margin x ref ± x margin f prec = fmargin f ref N u t1 x 1 α 1 d 1 N u cons t1
x 6 but along the line y = 2.

A part of the extrapolated values and coordinates of the local extrema for which 1.2 < α < 2.1 for the FD1 and FV3 solutions and 1.1 < α < 4.5 for the FE2 and FE4 solutions are given in Tables 345678, together with the associated α values and the relative distance d between the extrapolated and the finest grid values. In several cases (mainly concerning the FE4 solutions), only the values and the coordinates on the finest grid are provided because the associated α values are too high or too small compared with the α ranges given above. In these cases, the α values are replaced by * * * . References Table 8: Extrapolated values of f = (Nu t , Nu b ) local extrema along the line x = 30 and of their spanwise coordinates, y, and truncation error leading order, α, of RE. See the Table 6 caption for more details.

f ref ± f margin y ref ± y margin f prec = fmargin f ref N u t1 y 1 α 1 N u cons
Reference solutions for the local extrema and their coordinates have been determined in the same way as those of the integral values presented in §5.1.1. These reference solutions (denoted by f ref , x ref or y ref ), with their uncertainty margin (denoted by f margin , x margin or y margin ) and the precision of f ref (denoted by f prec ) are given in Tables 345678. Here, the reference values are equal to the arithmetic average of the minimum and maximum values of the extrapolated values of the four contributors, except when the FE2 or FE4 RE fails. In this case, the FE2 or FE4 extrapolated value is replaced by the FE2 or FE4 finest grid value. On the other hand, the FD1 and FV3 solutions on the finest grid are never used to define the reference solutions. Indeed, as it can be seen in Tables 3 and4, the relative distances d on θ, u, v, w for the FD1 and FV3 solutions are nearly always one or two orders higher than f prec whereas, for the FE2 and FE4 solutions, they are always on the same order or smaller than f prec (when RE is possible). This means that, with the definition and the precision of the reference values given here, RE of the discrete solutions obtained by the finite element methods are useless to determine the reference values of θ, u, v, w. RE is useful only to allow the second order FD1 and FV3 methods give solutions with a third order precision equal to the precision of the FE2 and FE4 methods.

The same observation as for the primitive variables can be done with the consistent Nusselt numbers Nu cons t and Nu cons b computed with the FE4 method. That is why only their values on the finest grid are proposed in the Tables 678. On the other hand, RE is useful to determine the reference values from the non consistent Nusselt numbers (compare d and f prec in Table 6 for instance). The extrapolated values of the non consistent Nusselt numbers and the finest grid values of the consistent Nusselt numbers are thus kept to define the reference values of the Nusselt numbers.

Following the methodology just described, the reference values of the primitive variable local extrema are given with four to five significant figures and those of the Nusselt number with three to four significant figures. Their coordinates are generally given with three significant figures in x direction and with four significant figures in y direction.

Heat and momentum fluxes through channel faces

Finally, we compare the dimensionless heat flux, Φ θ , and momentum fluxes, Φ u , Φ v and Φ w , through the boundary surfaces of the half channel obtained when the symmetry through the mid-plane at y = B/2 is taken into account. The flux definitions are given in Table 9. In this table, S i , S o , S f , S t and S b are respectively the inlet, outlet, front, top and bottom surfaces of the half channel, S s is the symmetry plane at y = B/2 and S tot = S i ∪ S o ∪ S f ∪ S s ∪ S t ∪ S b is the total surface of the half channel Ω/2. Note that, from the Navier-Stokes equation in [START_REF] Nicolas | Benchmark solution for a threedimensional mixed convection flow -Part 1: reference solutions[END_REF], the total momentum flux, Φ w , through S tot is equal to the integral of the buoyancy term, I buo = Ω/2

-Ra

Re²P r θ dV , on the half computational domain.

The diagonal terms of the momentum flux tensor depend on pressure. Since pressure is defined up to a constant, we decided to fix the value of this constant such that, for each grid, the momentum flux Φ u vanishes on the inlet surface S i . Furthermore, due to the symmetry conditions and our choice of boundary conditions, several other flux components are equal to zero. These are indicated in Table 9. We computed the fluxes in Table 9 on each grid then extrapolated these values by RE. The extrapolated values, or the values on the finest grid when RE fails, and the reference values with their tolerance margin are given in Table 10. The same criteria as those used in §5.1.1 and §5.2.3 to define the reference values and the tolerance margins of the integral quantities and local extrema are used here for the heat and momentum fluxes. In particular, the range of α kept to choose the extrapolated values used to define the references are 1 ≤ α ≤ 2.5 for the FD1 and FV3 solutions and 1 ≤ α ≤ 3.4 for the FE2 and FE4 solutions. Furthermore, the FD1 and FV3 solutions on the finest grid are not kept to define the reference. On the other hand, the FE2 and FE4 solutions on the finest grid are kept because, in this case, the relative distance between the extrapolated and finest grid solutions is still smaller than or on the same order as the precision of the reference solutions (not shown in Table 10). This methodology allows us to estimate the fluxes on the different surfaces with two to five significant digits, depending on the magnitude of the fluxes.

Φ θ Φ u Φ v Φ w S i S i -∂θ ∂x + RePr uθ dS S i p-
In Table 10, all the finite element fluxes are non consistent fluxes. The last part of the table gives the extrapolated total fluxes through S tot . The total flux balances are well satisfied since Φ θ,Stot , Φ u,Stot and Φ v,Stot are very near to zero and since Φ w,Stot is very near to the buoyancy term integral I buo3 . The maximum relative errors of the total flux balances, compared with the maximum elementary flux among the six elementary surfaces, are on the order of 10 -4 (the maximum relative error is equal to 10 -3 for Φ u through S tot computed with the FD1 method). Finally, one can note that the convergence

Φ θ Φ u Φ v Φ w FD1 S i -9.

Explanation of the observed behavior of Richardson Extrapolation

The preceding sections have shown a variety of behaviors when trying to apply Richardson extrapolation: (i) working behavior with an observed convergence order α equal to the consistency order α°of the discretization method; (ii) working behavior with an observed convergence order α located between 1 and α°; (iii) non-working behavior. In this section, we discuss these observations. The basic idea is to assume that the approximation error of a quantity f h can be written as two main terms within the Taylor expansion instead of one as in equation ( 3), section 4.1. Namely :

f h (h, C α°, α°, C r , r) = f exact + C α°h α°+ C r h r + O(h 1+max(α°,r) ) (10)
where C α°h α°i s the leading term of the approximation error of the regular part of the solution (same term as the one in equation ( 3)) and C r h r is the leading term of the approximation error of the singular part of the solution. As before °is the consistency order of the discretization and r measures the problem regularity influence on the actual convergence rate. Here f h → f exact when h → 0.

When h → 0, the term with largest exponent becomes negligible and Richardson extrapolation allows us to determine the smallest exponent and associated constant C as in section 4.1. However, in practice, we work with a fixed sequence of 3 (or 4) given h values, say {h 1 ; h

2 ; h 3 } = h 1 ; h 1 τ ; h 1 τ 2
with τ > 1. Scaling equation [START_REF] Pabiou | Wavy secondary instability of longitudinal rolls in Rayleigh-Bénard-Poiseuille flows[END_REF]

In the following, we use the scaled equation [START_REF] Nicolas | Linear stability of mixed convection flows in horizontal rectangular channels of finite transversal extension heated from below[END_REF], leaving out the tildes on f , h, Cα°a nd Cr for notation clarity. For example, our fixed sequence of (scaled) h, is now: 1; 1 τ ; 1 τ 2 . Then we choose typical numerical values C α°= 10 -4 << 1, τ = 2 and r = 1 and we plot the observed convergence order α as a function of the ratio ρ = Cr C α°w hen we apply the RE process (equation (4) of section 4.1) to our model function f h (equation ( 11)), neglecting the O(h 1+max(α°,r) ) term, in the four following cases : (i) ρ > 0, α°= 2; (ii) ρ > 0, α°= 3; (iii) ρ < 0, α°= 2; (iv) ρ < 0, α°= 3.

We also define a Richardson efficiency ratio σ as follows:

σ = log |f ex -f exact | max C α°h3 α° , |C r h 3 r | (12)
where f ex is the extrapolated function. If σ < 0, this means that RE has been successful in reducing the main component of the error compared to its value for the smallest h. Figure 12(a) (resp. Figure 12(b)) shows the profile of α and σ as a function of log |ρ| when ρ > 0 (resp. ρ < 0). On the two plots, we can distinguish three zones: Zone 1 log |ρ| -1 where the approximation error of the regular part of the solution dominates the approximation error of the singular part;

Zone 2 -1 log |ρ| 1 where the approximation error of the regular and singular parts of the solution have the same order of magnitude; Zone 3 1 log |ρ| where the approximation error of the singular part of the solution dominates. This zone corresponds to the asymptotic range in the present example since r = 1 and α°=2 or 3.

We can make the following remarks:

1. Zone 1 and 3 are the zones where RE is effective in reducing the error. This was expected for Zone 3 which is in the asymptotic range as defined in section 4.1, but not necessarily so for Zone 1;

2. In Zone 2, the behavior of RE depends on the sign of ρ: when ρ < 0, RE is not applicable, while if ρ > 0, RE still gives a result. However, as the profile of σ on Figure 12(a) shows, very little improvement in the reduction of the error is to be expected. We can conclude by saying that in Zone 2, RE is not very useful;

3. For the third-order methods, Zone 2 is larger than for the second-order methods (almost two decades instead of one); 4. As shown by the σ profiles, RE is less efficient at reducing the error for third-order methods than for second-order methods.

Remark 1 is consistent with our observations of the second order methods FD1 and FV3 for which RE seems to improve the results even though we are not in the asymptotic range. Therefore, most quantities seem to behave as if in Zone 1, with the notable exception of the mean temperature. Remark 3 is related to the fact that, for FE4 and FE2, RE was found to be much more difficult to apply than for the low order methods: this suggests that most quantities behave as if in Zone 2. This fact can be tempered with Remark 4 which suggests that less improvement in the error is to be expected for third order method than for second order method.

Applicability of RE was found to be higher for FE2 than for FE4: Remark 2 could provide an explanation for this suggesting that when close to or inside Zone 2, RE behavior can be found from non-working to almost working depending on the sign of ρ. Also, FE2 and FE4 not using the same finite elements, could be in different zones. In the test we have just presented, we have taken a mesh sequence with τ = 2, that is the cell length in each space direction is halved from one mesh to another. Remember that our fixed sequence of (scaled) h, is: 1; 1 τ ; 1 τ 2 . Now, we examine what happens when we take a smaller τ = 4 3 . This gives us figure 13a (resp. figure13b) which shows the profile of α and σ when ρ > 0 (resp. ρ < 0). The main thing to notice is that Zone 2 seems to be smaller for τ = 4 3 than for τ = 2. This is interesting because Zone 2 is the zone where RE doesn't work reliably. We explain this fact as follows: by taking a more condensed sequence of h, it is easier to be in the same convergence zone for all the hs.

To conclude this section, although further investigations would be necessary to assess that the proposed explanations are the right ones, we have found that the simple model in equation [START_REF] Nicolas | Linear stability of mixed convection flows in horizontal rectangular channels of finite transversal extension heated from below[END_REF] allowed us to reproduce most of the behavior we have observed in trying to apply RE to the benchmark of this article.

Conclusions

In this paper, a methodology has been presented in details to establish a first numerical benchmark solution of a three-dimensional mixed convection flow in a horizontal rectangular channel, partially heated from below and cooled from above. This methodology is based on the use of four different numerical methods (second order finite difference and finite volume methods, and third order finite element methods), Richardson extrapolations (RE) on very fine grids and cubic spline interpolations. Reference solutions are proposed for the dynamical and thermal fields, in the form of local, integral or differential quantities such as local extrema of the primitive variables and Nusselt numbers, surface heat and momentum fluxes, volume integrals of the temperature and kinetics energy, etc. These reference solutions are generally given up to four or five significant figures.

The difficulty of the present benchmark problem is that a discontinuity takes place in the thermal gradient over the bottom plate at x = 0, which not only significantly restricts the conditions of application of RE to establish reference solutions, but also complicates its analysis. Therefore the theoretical basis of RE are reminded and discussed from the viewpoint of this singularity. It is shown that the convergence order, α, observed from RE of the local and integral quantities is reduced to one in the neighborhood of the boundary condition discontinuity and tends to the consistency order, α°=2 or 3, of the used discretization methods far from the singularity. It is deduced from this result and other test cases that the problem regularity is close to r = 1 in the vicinity of the boundary condition discontinuity. Moreover, we have suggested in §5.4 a modified Taylor expansion to account for the problem singularity in the RE formalism. A simplified model problem has enabled us to reproduce most of the behaviors we have observed in the former benchmark problem and helped us to understand them.

The paper has also brought to the fore several practical difficulties in the proper usage and implementation of RE. It has been shown that the distance between the finest grid solutions and the extrapolated solutions is much smaller for the FE2 and FE4 methods than for the FD1 and FV3 ones. Furthermore the local behavior of α is much more oscillatory for the two third order methods than for the two second order ones. It was also shown that, for the used grids, RE cannot be applied locally on the whole fields due to the "crossing" of the computed quantities on the different grids. The FE2 and FE4 solutions have appeared very sensitive to these field variations and this behavior has been understood thanks to the model problem introduced in §5.4. Using larger grid size ratios (resulting in much finer grids) and therefore much greater computational costs than those required for the present paper would also have probably be another way to solve this problem.

Furthermore, this study has reminded us that the convergence order of a numerical model can be significantly deteriorated due to a loss of regularity of the solution and that the standard RE framework should not be used without taking precautions in this case. Finally, it is noteworthy that the four numerical models used for this benchmark have displayed their own sensibility to the various problem peculiarities (establishment zone, localized thermal gradient singularity, etc.) and, wherever the RE has been found to be applicable, the resulting convergence order could also depend on the quantity (primitive or derived variable) it is based on and its definition (L 2 norm, mean value, etc.).

Appendix B : about the mesh accuracy in the vertical direction and its consequences on the wall fluxes

In the cases of Rayleigh-Bénard free convection or PRB mixed convection, the thermal and velocity gradients are present on the whole thickness of the fluid layer due to the development of the thermoconvective rolls. This is clearly shown in Figure 14 in which the isolines of the temperature and velocity components are drawn in a half transversal plane at x = 45, with a constant step between two isolines in each sub-figure. This explains why uniform meshes are generally used in the vertical and spanwise directions when simulating such flows. As is shown in Figure 14, the temperature and velocity gradients are maximum in the z-direction near the top and bottom walls on a thickness equal to a quarter of the channel height. Therefore, in the present study, there is at minimum 10 nodes/cells/points in this zone since the coarsest grid of each contributor has 40 nodes/cells/points at minimum on the channel height. We have checked that this is enough to get accurate wall fluxes for the Reynolds and Rayleigh number values used in the present benchmark problem. For instance, in Figure 15, we present the space evolution of the local Nusselt number along two lines on the bottom plate of the channel, computed on the four grids (h 1 , h 2 , h 3 , h 4 ) of the FD1 solution. The extrapolation coefficient α is also given as well as the extrapolated Nusselt number values (when RE succeeds). One can see that the maximum relative discrepancy between the Nusselt number values on the coarsest (h 1 ) and finest (h 4 ) grids is at most 1%, while it is much less for the third order FE methods.

Figure 1 :

 1 Figure 1: Geometry and top and bottom thermal boundary conditions (the vertical lateral walls are adiabatic). The red dashed lines are some of the lines along which θ, u, v and w profiles and their extrema are calculated.

  and there is a Neumann thermal boundary condition, ∂θ/∂z = 0, for x ∈ [-A e , 0[, next to the Dirichlet condition, θ = 1, for x ∈ [0, A -A e ]; • at z = 1, -→ v = -→ 0 and there is a Neumann thermal boundary condition, ∂θ/∂z = 0, for x ∈ [-A e , 0[, next to the Dirichlet condition, θ = 0, for x ∈ [0, A -A e ];

Figure 2 :

 2 Figure 2: Temperature fields, θ, in the horizontal mid-plane and temperature field and velocity vector field in the transverse plane at x = 30 in the longitudinal roll flow of the present test case. The yellow dashed lines are some of the lines along which θ, u, v and w profiles and their extrema are calculated.

Figure 4 :

 4 Figure4: FD1 solutions of the vertical velocity component streamwise profiles along the line at (y, z) = (B/5, 0.5), computed on the four grids described in Table1, together with the extrapolated solution. The latter can be considered as the asymptotic solution here because the α values monotonously varies between 2.01 and 2.16 when x varies between 24.7 and 25.3. In this figure, all the symbols correspond to the computational or extrapolation points and the curves to the cubic spline interpolation curves. The small window zooms in the local extremum of the extrapolated curve. The value w 2 = -0.472989 and the coordinate x 2 = 24.907 of this extremum are reported in Table3.

  L 2 norm of vertical velocity component W L 2

Figure 5 :

 5 Figure 5: Space evolution of the convergence orders α observed from RE of the integral quantities T m , T L 2 and W L 2 in various upstream or downstream subdomains of the computational domain. 18

  2.3). In the present study, the FE2 Nusselt numbers are the non consistent ones while the FE4 contribution proposes the two Nusselt number types. The consistent Nusselt numbers will be denoted by Nu cons t and Nu cons b while the notations Nu t and Nu b will be kept to denote the non consistent Nusselt numbers and to denote the Nusselt numbers in a generic way.

Figure 6 : 5 .

 65 Figure 6: Longitudinal profiles of θ, u, v and w along the lines at (y, z) = (2, 0.2) and (y, z) = (5, 0.5) and longitudinal profiles of Nu t and Nu b along the lines at y = 2 and y = 5. The filled circles indicates the local extrema whose values and coordinates are given in Tables 3, 4, 6 and 7.

Figure 7 :

 7 Figure 7: Spanwise profiles of the primitive variables θ, u, v and w along the lines at (x, z) = (10, 0.2), (10, 0.5), (30, 0.2) and (30, 0.5). The filled circles indicates the local extrema whose values and coordinates are given in Table5.

  Figure 7: Spanwise profiles of the primitive variables θ, u, v and w along the lines at (x, z) = (10, 0.2), (10, 0.5), (30, 0.2) and (30, 0.5). The filled circles indicates the local extrema whose values and coordinates are given in Table5.

Figure 8 :

 8 Figure 8: Transverse profiles of the Nusselt numbers Nu t and Nu b along the lines at x = 10 and x = 30. The filled circles indicates the local extrema whose values and coordinates are given in Table8.

  Figure 8: Transverse profiles of the Nusselt numbers Nu t and Nu b along the lines at x = 10 and x = 30. The filled circles indicates the local extrema whose values and coordinates are given in Table8.

Figure 9 :

 9 Figure9: FD1 solutions of the temperature streamwise profiles along the line at (y, z) = (B/5, 0.5), computed on the four grids described in Table1, together with the extrapolated solution. The latter does not tend to the asymptotic solution near x = 22.35 because the α values (the black filled circles in the figure) diverge where the profiles intersect. All the symbols correspond to the computational or extrapolated points. The curves correspond to the cubic spline interpolation curves, except for the extrapolated solution where they are linear segments linking the extrapolated points.

Figure 10 :

 10 Figure 10: Comparison of the streamwise profiles of the extrapolation coefficients, α, computed by the four contributors for θ, u, v and w, along the line at (y, z) = (2, 0.2).The corresponding profiles of θ f g , u f g , v f g and w f g on the finest grid are also drawn on the figures.

Figure 11 :

 11 Figure 11: Comparison of streamwise and spanwise profiles of the extrapolation coefficients, α, for Nu t computed by the four contributors. The corresponding profiles of Nu t-f g on the finest grid are also drawn on the figures. 25

  For FE4 contribution, Nu t and Nu b are the extrapolated values of the non consistent Nusselt numbers and Nu cons t and Nu cons b are the consistent Nusselt number values on the finest grid. The symbol °indicates values that are excluded from the reference determination.

Table 9 :

 9 Definition of the heat and momentum fluxes through the channel faces

  α°= 2 or 3 and ρ < 0.

Figure 12 :

 12 Figure 12: Profiles of α and σ as a function of log |ρ| for τ = 2.

  α°= 2 or 3 and ρ < 0.

Figure 13 :

 13 Figure 13: Profiles of α and σ as a function of log |ρ| for τ = 4 3 .

Figure 14 :

 14 Figure 14: Isolines of the temperature, θ, and velocity components, u, v and w, in a half spanwise vertical plane at x = 45. A constant step is used between two isolines in each sub-figure.

  Streamwise Nu profiles along (y, z) = (5, 0) Spanwise Nu profiles along (x, y) =[START_REF] Roy | Review of code and solution verification procedures for computational simulation[END_REF] 0)

Figure 15 :

 15 Figure 15: FD1 streamwise and spanwise profiles of the local Nusselt number, along the lines (y, z) = (5, 0) and (x, z) = (30, 0), on the four different grids used (h 1 , h 2 , h 3 , h 4 ), as well as the extrapolated Nu profile (h ext ) and its extrapolation coefficients α.

Table 1 :

 1 Contributor N x × N y × N z [symmetry]

			∆t	User time	Consistency
				[computer type (organism/lab)]	orders α°M
	SME,	400 × 134 × 40	0.01	36 min on 1 processor	2 for θ, u, v, w;
	FD1	600 × 200 × 60	0.01	2 h 20 on 1 processor	2 for p
		800 × 268 × 80	0.002	25 h on 1 processor
		1200×400×120	0.002	100 h on 1 processor
		[no]		[NEC SX5 (IDRIS)]
	IUSTI,	601 × 121 × 41	0.01	19 min on 60 processors	3 for θ, u, v, w;
	FE2	901 × 181 × 61	0.01	1 h 40 on 150 processors	2 for p
		1351 × 271 × 91	0.005	43 h 15 on 225 processors
		[yes]		[IBM SP6 (IDRIS)]
	I2M	601 × 161 × 41	0.1	8 h on 152 processors	2 for θ, u, v, w;
	Institute,	901 × 241 × 61	0.1	12 h on 152 processors	2 for p
	FV3	1351 × 361 × 91	0.1	56 h on 152 processors
		[yes]		[ALTIX ICE 8200 (I2M Inst.)]
	CEA, FE4 601 × 121 × 49	0.5	200 h on 8 processors	3 for θ, u, v, w;
		751 × 151 × 61	0.5	400 h on 8 processors	2 for p
		801 × 161 × 65	0.5	450 h on 8 processors
		1001 × 201 × 81	0.5	1600 h on 8 processors
		[yes]		[PC 8 cores (CEA)]

Numerical parameters used by the different contributors the outflow boundary. Steady state solutions are obtained by integrating long enough.

  .

	-0.4705
		Nx=400
		Nx=600
	-0.471	extrapolated solution Nx=800 Nx=1200	-0.47298
			-0.472985
	-0.4715
	w		-0.47299	24.88	24.9	24.92 24.94
	-0.472
	-0.4725
	-0.473
	24.7	24.8	24.9	25 x		25.1	25.2	25.3

Table 2 :

 2 Left columns: finest grid (f f g ) and extrapolated (f ex ) values of the integral functions f = 2E c , ∆P io and T m , truncation error leading order, α f , from their RE and relative distance,

Table 3 :

 3 Extrapolated values or finest grid values (indicated by * * * ) of f = (

Table 4 :

 4

		FD1	FE2	FV3	FE4	References
						f ref ± f margin
						y ref ± y margin
						f prec =	fmargin f ref
	θ 1	0.24719	0.24716	0.24714	0.24715	0.24716 ± 0.00003
	y 1	1.0361	1.0368	1.0363	1.0365	1.0364 ± 0.0004
	α 1	1.90	2.91	1.96	3.38	1.2 × 10 -4
	θ 2	0.23816	0.23806	0.23816	0.23807	0.23811 ± 0.00005
	y 2	2.8991	2.9000	2.8993	2.8990	2.8995 ± 0.0005
	α 2	1.92	2.90	1.99	* * *	2.1 × 10 -4
	θ 3	0.77385	0.77384	0.77387	0.77383	0.77385 ± 0.00002
	y 3	3.9042	3.9039	3.9042	3.9041	3.9041 ± 0.0002
	α 3	1.94	2.86	2.09	* * *	2.6 × 10 -5
	u 1	1.06513	1.06513	1.06499	1.06509	1.06506 ± 0.00007
	y 1	1.0087	1.0087	1.0086	1.0086	1.0086 ± 0.0001
	α 1	1.95	2.87	2.01	* * *	6.6 × 10 -5
	u 2	1.15071	1.15059	1.15049	1.15062	1.15060 ± 0.00011
	y 2	3.9004	3.9005	3.9002	3.9002	3.9003 ± 0.0002
	α 2	1.94	2.96	1.93	* * *	9.6 × 10 -5
	u 3	1.74979	1.74962	1.74987	1.74978	1.74975 ± 0.00013
	y 3	4.4425	4.4425	4.4425	4.4425	4.4425 ± 0.0000
	α 3	1.94	2.94	2.02	* * *	7.4 × 10 -5
	v 1	0.0288474	0.0288520	0.0288523	0.0288521	0.0288499 ± 0.0000025
	y 1	0.4081	0.4076	0.4078	0.4088	0.4082 ± 0.0006
	α 1	1.93	* * *	2.09	* * *	8.7 × 10 -5
	v 2	0.035892	0.035912	0.035916	0.035900	0.035904 ± 0.000012
	y 2	0.7049	0.7043	0.7050	0.7051	0.7047 ± 0.0004
	α 2	1.97	3.15	2.03	* * *	3.3 × 10 -4
	v 3	0.032867	0.032907	0.032917	0.032878	0.032892 ± 0.000025
	y 3	4.7390	4.7395	4.7390	4.7388	4.7391 ± 0.0004
	α 3	1.96	* * *	1.67	3.62	7.6 × 10 -4
	w 1	0.372397	0.372372	0.372496	0.372397	0.37243 ± 0.00006
	y 1	0.2286	0.2285	0.2285	0.2285	0.2285 ± 0.0001
	α 1	2.00	2.91	1.98	* * *	1.6 × 10 -4
	w 2	-0.513013	-0.513049	-0.513221	-0.513050	-0.51311 ± 0.00011
	y 2	2.8998	2.8995	2.9000	2.8999	2.8998 ± 0.0003
	α 2	1.91	2.85	2.00	* * *	2.1 × 10 -4
	w 3	0.490347	0.490335	0.490478	0.490347	0.49041 ± 0.00007
	y 3	3.9028	3.9029	3.9027	3.9028	3.9028 ± 0.0001
	α 3	1.95	2.91	1.96	* * *	1.4 × 10 -4

Same as Table

3

but for f = (θ, u, v, w) along the line (y, z) = (2, 0.2).

Table 5 :

 5 Extrapolated values or finest grid values (indicated by * * * ) of f = (θ, u, v, w)

Table 6 :

 6 Extrapolated values or finest grid values (indicated by * * * ) of f = (Nu t , Nu b )

Table 7 :

 7 Same as Table

Table 10 :

 10 , α, of I buo are equal to those of T m in Table2since I buo and T m are proportional quantities.

			092 × 10 -9 [ * ]°0.0	-2.1364 × 10 -2 [1.07]	6.931 × 10 -5 [1.00]
	FE2	S i	-1.011 × 10 -8 [ * ]	0.0	-2.1344 × 10 -2 [2.72]	7.042 × 10 -5 [2.76]
	FV3	S i	-1.17 × 10 -8 [ * ]°0.0	-2.1348 × 10 -2 [2.18]	7.144 × 10 -5 [ * ]°F
	E4	S i	-1.031 × 10 -8 [ * ]	0.0	-2.1346 × 10 -2 [ * ]	7.058 × 10 -5 [ * ]
	Ref.	S i	-(1.021±0.010)×10 -8	0.0	-(2.1354±0.0010)×10 -2	(7.00 ± 0.07) × 10 -5
	FD1	So	-87.6290 [1.99]	72.1696 [2.03]	3.0294 × 10 -2 [0.88]	1.6804 × 10 -2 [1.02]
	FE2	So	-87.6274 [2.54]	72.1710 [2.80]	3.1236 × 10 -2 [2.37]	1.6576 × 10 -2 [2.89]
	FV3	So	-87.6328 [2.27]	72.1711 [1.99]	3.0356 × 10 -2 [1.56]	1.6822 × 10 -2 [1.86]
	FE4	So	-87.6292 [ * ]	72.1703 [0.93]	3.0262 × 10 -2 [1.04]	1.6828 × 10 -2 [ * ]
	Ref.	So	-87.630 ± 0.003	72.1704 ± 0.0008	(3.07 ± 0.05) × 10 -2	(1.670 ± 0.013) × 10 -2
	FD1	S f	∼ 0.0°-3.9727 [0.81]°-409.362 [2.00]	-1.7570 [0.95]°F
	E2	S f	-0.00303 [2.46]°-3.9804 [2.58]	-409.311 [2.56]	-1.7672 [3.38]
	FV3	S f	∼ 0.0°-3.9878 [1.77]	-409.366 [1.91]	-1.7683 [1.01]
	FE4	S f	0.04744 [1.69]°-3.9884 [1.65]	-409.381 [ * ]	-1.7675 [3.17]
	Ref.	S f	0.0	-3.984 ± 0.004	-409.35 ± 0.04	-1.7678 ± 0.0006
	FD1	Ss	∼ 0.0°∼ 0.0°409.307 [2.07]	∼ 0.0°F
	E2	Ss	-3.996 × 10 -3 [2.91]°-3.96 × 10 -4 [2.88]°409.301 [2.66]	1.027 × 10 -6 [3.39]°F
	V3	Ss	∼ 0.0°-7.928×10 -6 [2.90]°409.318 [2.03]	5.499 × 10 -7 [3.13]°F
	E4	Ss	4.547 × 10 -4 [2.88]°5.848 × 10 -5 [2.85]°409.316 [2.45]	1.989 × 10 -6 [3.10]°R
	ef.	Ss	0.0	0.0	409.31 ± 0.01	0.0
	FD1	S b	479.993 [2.07]	-35.4293 [1.16]	2.6366 [2.37]	-2249.602 [2.12]
	FE2	S b	480.018 [2.64]	-35.4043 [3.22]	2.6370 [2.33]	-2249.595 [2.26]
	FV3	S b	479.741 [0.76]°-35.4122 [1.08]	2.6358 [1.84]	-2249.678 [2.03]
	FE4	S b	479.928 [2.49]	-35.4029 [3.20]	2.6373 [2.28]	-2249.678 [2.85]
	Ref.	S b	479.97 ± 0.05	-35.416 ± 0.013	2.6366 ± 0.0008	-2249.64 ± 0.04
	FD1	St	-392.228 [2.17]	-32.7936 [1.08]	-2.5867 [2.55]	1930.952 [2.05]
	FE2	St	-392.368 [2.70]	-32.7810 [3.18]	-2.5864 [2.65]	1931.972 [2.26]°F
	V3	St	-392.382 [1.42]	-32.7861 [1.11]	-2.5868 [2.11]	1930.982 [1.99]
	FE4	St	-392.293 [2.68]	-32.7793 [3.15]	-2.5871 [2.50]	1930.979 [2.53]
	Ref.	St	-392.31 ± 0.08	-32.786 ± 0.008	-2.5868 ± 0.0004	1930.967 ± 0.015
	FD1	Stot	0.1140 [1.37]°0.07796 [0.83]°2.039 × 10 -4 [ * ]°-320.4334 [2.07]
	FE2	Stot	0.03801 [2.29]°0.00560 [3.86]°5.120 × 10 -2 [2.54]°-319.5350 [2.62]°F
	V3	Stot	-0.3576 [0.49]°0.007652 [ * ]°3.159 × 10 -3 [1.13]°-320.4303 [1.03]
	FE4	Stot	0.03687 [1.72]°-0.009371 [ * ]°2.666 × 10 -3 [1.61]°-320.4572 [ * ]
	Ref.	Stot	0.0	0.0	0.0	-320.444 ± 0.014
	FD1	I buo				-320.4245 [1.19]
	FE2	I buo				-320.4318 [1.18]
	FV3	I buo				-320.4330 [1.02]
	FE4	I buo				-320.4375 [1.18]
	Ref.	I buo				-320.431 ± 0.007

Extrapolated values of the heat and momentum fluxes through channel surfaces and buoyancy term integral. The values of the extrapolation coefficient α are noted by [ ]. Extrapolation failure is indicated by [ * ] and the flux value on the finest grid is given. When the flux value is strictly smaller than 10 -11 , it is noted ∼ 0.0. The symbol °indicates values that are excluded from the reference determination. The reference values 0.0 correspond to the theoretical values given in Table

9

. orders

  with fh = f h

fexact and h = h h 1 , one gets: fh ( h, Cα°, α°, Cr , r) = 1 + Cα°h α°+ Cr hr + O( h1+max(α°,r) )

the consistency order is the formal convergence order that is the leading order of the space discretization truncation error.

To explain the temperature gradient discontinuity, let's consider the isotherms near this singular point: the left-hand isotherms tend to be perpendicular to the wall due to the homogeneous Neumann condition while the right-hand isotherms tend to be parallel to the wall due to the constant Dirichlet condition. With a such change in the thermal boundary conditions, to avoid any singularity, a wall with a convex 90°step would be necessary.

For the FE2 solution, the term Φ w,Stot is slightly different from the term I buo because the value of Φ w on S t is overestimated. Indeed, the expected improvement associated with the consistent non homogeneous Neumann boundary condition for the pressure correction (see §3.2) turns out to become a weakness for the coarsest meshes considered (M1 and probably M2) to perform RE. Indeed, the normal projection of the second order derivative term associated with the viscous stress is quite stiff in the boundary layers and the M1 mesh resolution does not enable to accurately compute this term, which results in a poor pressure correction close to the horizontal walls where both thermal and hydrodynamic boundary layers develop. Further investigations are on the way to implement in a more subtle way this term for Q2 finite elements especially for coarse mesh problems. That is why the FE2 values of Φ w,St and Φ w,Stot are excluded from the reference value determination in Table10.

Acknowledgments

The first author gratefully acknowledges Pr. Shihe Xin, from CETHIL, UMR 5008 CNRS/Insa-Lyon, France, to have provided the finite difference code, FD1, that was developed by him when he was at LIMSI, CNRS, UPR 3251, Orsay, France. The authors acknowledge Donna Calhoun for proof reading of the article and providing many useful suggestions. This work was granted access to the HPC resources of IDRIS (NEC-SX5 vectorial supercomputer and IBM SP4 and SP6 parallel supercomputers) under the allocations 2006-1823 and 2007-1823 made by GENCI (Grand Equipement National de Calcul Intensif). The fourth author thanks the Aquitaine Regional Council for the financial support dedicated to a 256-processor cluster investment, located at I2M Institute.

Appendix A: about the computation of the analytical expression of the Poiseuille profile in a 3D channel

The analytical solution of the dimensionless Poiseuille profile is equal to u P ois (y, z) = U P ois (y,z)

Umean , where the dimensional Poiseuille profile U P ois (y, z) is given by [START_REF] Nicolas | Linear stability of mixed convection flows in horizontal rectangular channels of finite transversal extension heated from below[END_REF]:

where

∂P ∂X is the average velocity of the "two-dimensional" Poiseuille flow, i.e. in a two-dimensional channel or between two infinite plates, and where the dimensional average velocity U mean of the Poiseuille flow is given by:

Thus the inlet dimensionless Poiseuille profile u P ois (y, z) is given by the ratio of equations ( 13) and ( 14), where about 25 terms are kept in the series to obtain a sufficiently accurate entrance velocity profile. Note that in ( 13) and ( 14), the hyperbolic cosine at the denominator diverges when n is high. To avoid any problem, the two hyperbolic cosines of the series can be transformed in real exponentials via the Euler relations. Thus, by denoting N = 2n + 1, the ratio of the two hyperbolic cosines writes:

To avoid the divergence of exp(-2Nπ(y -B 2 ) when N is high, the Poiseuille profile must only be computed for B 2 ≤ y ≤ B. The symmetry with respect to the median vertical plane is used to compute the Poiseuille profile for 0 ≤ y ≤ B 2 .