A Fadai-Ghotbi

R Manceau
email: remi.manceau@lea.univ-poitiers.fr

J Borée

Revisiting URANS computations of the backward-facing step flow using second moment closures. Influence of the numerics

Keywords: Turbulence modelling, Unsteady flow, Near-wall flow, Elliptic Blending

Summary. Previous studies have shown that Unsteady Reynolds-Averaged Navier-Stokes (URANS) computations are able to reproduce the vortex shedding behind a backward-facing step. The aim of the present work is to investigate not only the quantitative predictions of the URANS methodology concerning the characteristic frequencies, but also the amplitude of the energy of the resolved eddies, by using the Elliptic Blending Reynolds Stress Model. This innovative low-Reynolds number second moment closure reproduces the non-viscous, non-local blocking effect of the wall on the Reynolds stresses, and it is compared to the standard k-ε and LRR models using wall-functions. Consistent with previous studies, in the 2D computations shown in the present article, the vortex shedding is captured with the correct Strouhal number, when second moment closures are used. To complete these previous analyses, we particularly focus here on the energy contained in the unsteady, resolved part and its dependency on the numerical method. This energy is less than 5% of the total energy and is strongly dependent on the mesh. Using a refined mesh, surprisingly, a steady solution is obtained. It is shown that this behaviour can be linked to the very small spatial oscillations at the step corner, produced by numerical dispersion, which act as perturbations that are sufficient to excite the natural mode of the shear layer, when the local Peclet number, comparing convection and diffusion effects, is high enough. This result suggests that URANS is not appropriate to quantitatively predict the amplitude of the large-scale structures developing in separated shear-layers, and that URANS results must be interpreted with care in terms of temporal variations of forces, temperatures, etc., in industrial applications using marginally fine meshes.

Introduction

Recently, a considerable interest has emerged in relatively low-cost (compared to Large Eddy Simulation) time-dependent computations of complex flows, in different fields of industrial activities, e.g., for thermal fatigue studies, fluid/structure interaction, noise prediction, etc. In particular, the URANS methodology has become quite popular, because of its successes in predicting the most energetic modes and their Strouhal number, in particular for vortex shedding behind bluff bodies [START_REF] Iaccarino | Reynolds averaged simulation of unsteady separated flow[END_REF][START_REF] Bosch | Simulation of vortex shedding past a square cylinder with different turbulence models[END_REF][START_REF] Johansson | Numerical simulation of vortex shedding past triangular cylinders at high Reynolds number using a k-ε turbulence model[END_REF][START_REF] Durbin | Separated flow computations with the k-ε-v 2 model[END_REF][START_REF] Lardeau | Unsteady RANS modelling of wakeblade interaction: computational requirements and limitations[END_REF]. The methodology is straightforward: it simply consists in solving standard RANS-model equations in a time-accurate mode. However, fundamental questions remain about the general definition of the URANS decomposition when the flow is stationary, since the ensemble average is independent of time, while URANS computations provide an unsteady solution. Moreover, even for non-stationary flows, like pulsed flows or flows around moving obstacles, for which the ensemble average is time-dependent, it is also worth pointing out that using standard RANS models in an unsteady computation is questionable. Indeed, although they sometimes provide satisfactory solutions, they are mainly formulated under the assumption of equilibrium, which can lead to a significant misrepresentations of the physical mechanisms present in unsteady flows: for instance, Carpy and Manceau [START_REF] Carpy | Turbulence modelling of statistically periodic flows: Synthetic jet into quiescent air[END_REF] have shown that linear eddy-viscosity models give a completely wrong dynamics in a pulsed jet, due to the lack of representation of the stress-strain misalignment. On the contrary, second-moment closures (SMCs) are able to reproduce the global dynamics of this flow without any case-specific modification. In the case of separated shear layers, it was shown by Lasher & Taulbee [START_REF] Lasher | On the computation of turbulent backstep flow[END_REF] that SMCs are able to reproduce the vortex shedding in a backward-facing step flow, as soon as sufficiently accurate numerical methods are used. Despite the fact that separation is fixed by the geometry, which avoids the appearance of an unsteady separation point, as can be observed on smooth surfaces, the backward-facing step flow is a challenging test case for unsteady computations, since it conjugates several fundamental mechanisms: vortex shedding, convection and pairing of these structures downstream, interaction between the vortices and the wall, flapping of the recirculation region at a low frequency. Numerous RANS studies (steady computations) of backwardfacing step flows (e.g., [START_REF] Abe | New turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows-1. Flow field calculations[END_REF][START_REF] Craft | Developments in a low-Reynolds-number second-moment closure and its application to separating and reattaching flows[END_REF][START_REF] Hanjalić | Contribution towards the second-moment closure modelling of separating turbulent flows[END_REF][START_REF] Parneix | A procedure for using DNS databases[END_REF]) have shown the difficulty in reproducing the main characteristics of the flow (recirculation length, backflow intensity, boundary-layer recovery), especially for small expansion ratios, where turbulence has a strong influence on the mean flow. On the other hand, only Lasher & Taulbee [START_REF] Lasher | On the computation of turbulent backstep flow[END_REF] have investigated this flow using a time-dependent computation (URANS) with a second-moment closure. However, their model was based on wall functions, and the initial aim of the present work was to evaluate, by performing 2D computations, the performance of a recent near-wall SMC, the Elliptic Blending Reynolds Stress Model (EB-RSM) [START_REF] Manceau | Elliptic Blending Model: a new near-wall Reynolds-stress turbulence closure[END_REF][START_REF] Manceau | An improved version of the Elliptic Blending Model. Application to non-rotating and rotating channel flows[END_REF], formulated in order to replicate the non-viscous, non-local blocking effect of the wall. This evaluation does not only concern the mean quantities, but also the quantitative prediction of the characteristic frequencies and the kinetic energy content of the resolved, unsteady large-scale structures. This study has actually revealed that the most important aspect of such computations is the influence of the numerics on the nature of the solution (steady/unsteady). Therefore, the main part of the present work is devoted to the analysis and understanding of this influence, and to the identification of the limitations of URANS for this type of flows.

Governing equations

The instantaneous flow is driven by the Navier-Stokes and continuity equations. The instantaneous velocity U * i is decomposed into a resolved part Ũi , including mean value and large-scale fluctuations, and a residual fluctuating part u ′′ i such that

U * i = Ũi + u ′′ i . (1)
The resolved velocity is obtained by the convolution product of a filter G with the instantaneous velocity as [START_REF] Gatski | Current trends in modeling research for turbulent aerodynamic flows[END_REF] Ũi (x, t)

= U * i = r∈V τ ∈[-∞,t] G(x -r, t -τ )U * i (r, τ ) drdτ, ( 2
)
where V is the fluid domain and the brackets denote the URANS filter, which can be defined as a phase average in the present case, corresponding to a filter function given by

G(x, t) = δ(x) lim N →∞ 1 N + 1 N n=0 δ(t + nT 0 ) , (3)
with T 0 the shedding time scale. The long-time average of U * i is denoted by U i , so that the large-scale fluctuation is u ′ i = Ũi -U i , and the total fluctuation is

u i = U * i -U i = u ′ i + u ′′ i .
The URANS continuity and Navier-Stokes equations are formally identical to the RANS equations: the unknown correlation u ′′ i u ′′ j , appearing in the momentum equation concerning Ũi , is the residual stress tensor, and has to be modelled. Using a second moment closure, the transport equations for u ′′ i u ′′ j and the kinetic energy dissipation rate ε can be written in a general form x/h

Figure 1: Wall-normal unit vector at the step corner, calculated from Eq. ( 14).

Ũi = 0; < u ′′ i u ′′ j >= 0; ε = 2ν lim y→0 k ′′ y 2 ; α = 0 Table 1: Boundary conditions. D Dt u ′′ i u ′′ j = Pij + φ ij -ε ij + ∂ l νδ lm + C µ T σ k u ′′ l u ′′ m ∂ m u ′′ i u ′′ j , (4)
Dε Dt = C ′ ε 1 P T -C ε 2 ε T + ∂ l νδ lm + C µ T σ ε u ′′ l u ′′ m ∂ m ε. (5)
where D/ Dt = ∂ t + Ũk ∂ k . The stress production term Pij and P = 1 2 Pkk are exact and do not need to be modelled; ε ij denotes the dissipation rate tensor, with ε = 1 2 ε kk ; φ ij is the velocity-pressure-gradient correlation, known as the pressure term. In standard models, this term is modelled using strong hypotheses like locality and local quasi-homogeneity, not valid in wall-bounded flows [START_REF] Manceau | Inhomogeneity and anisotropy effects on the redistribution term in Reynolds-Averaged Navier-Stokes modelling[END_REF]. In order to avoid the use of these hypotheses, Durbin [START_REF] Durbin | Near-wall turbulence closure modeling without "damping functions[END_REF][START_REF] Durbin | A Reynolds stress model for near-wall turbulence[END_REF] proposed the elliptic relaxation model, derived from the exact integral equation for φ ij [START_REF] Durbin | Near-wall turbulence closure modeling without "damping functions[END_REF][START_REF] Manceau | Inhomogeneity and anisotropy effects on the redistribution term in Reynolds-Averaged Navier-Stokes modelling[END_REF], which implies the resolution of six differential equations for

φ ij φ ij -L 2 ∇ 2 φ ij = φ h ij , (6)
where φ h ij is any "homogeneous" (valid away from the wall) model. In order to avoid the resolution of six additional equations with numerically stiff boundary conditions, Manceau [START_REF] Manceau | An improved version of the Elliptic Blending Model. Application to non-rotating and rotating channel flows[END_REF] suggests, following Manceau & Hanjalić [START_REF] Manceau | Elliptic Blending Model: a new near-wall Reynolds-stress turbulence closure[END_REF], to blend the homogeneous model and a near-wall model of φ ij as

φ ij = (1 -α 2 )φ w ij + α 2 φ h ij , (7)
where α is a blending coefficient which goes from zero at the wall to unity far from the wall. Contrary to Eq. ( 6), Eq. ( 7) is not a differential, elliptic equation, and consequently, in order to preserve the non-local influence of the wall, which yields the progressive transition from the homogeneous form to the near-wall form, Manceau [START_REF] Manceau | An improved version of the Elliptic Blending Model. Application to non-rotating and rotating channel flows[END_REF] proposes to use an elliptic equation for α, similar to the elliptic relaxation equations

α -L 2 ∇ 2 α = 1. ( 8
)
For the "homogeneous" part of the pressure term φ h ij , any existing model can be used. The SSG model [START_REF] Speziale | Modeling the pressure-strain correlation of turbulence: an invariant dynamical system approach[END_REF] is used here y/h = 1 with

φ h ij = -g 1 + g * 1 P ε εb ′′ ij + g 3 -g * 3 b ′′ kl b ′′ kl k ′′ Sij + g 4 k ′′ b ′′ ik Sjk + b ′′ jk Sik - 2 3 b ′′ lm Slm δ ij + g 5 k ′′ b ′′ ik Ωjk + b ′′ jk Ωik , (9)
C ε1 = 1.44; C ε2 = 1.83; C µ = 0.21; σ k = 1.0; σ ε = 1.15; A 1 = 0.03; C L = 0.161; C η = 80.0; C T = 6.0; g 1 = 3.4 ; g * 1 = 1.8 ; g 3 = 0.8 ; g * 3 = 1.3; g 4 = 1.25 ; g 5 = 0.4
y/h = 9 y/h = 0 x/h = -4 x/h = 0 x/h = 40 h x y U 0 Inlet Outlet
k ′′ = 1 2 u ′′ j u ′′ j , b ′′ ij = u ′′ i u ′′ j 2k ′′ - 1 3 δ ij , (10)
Sij = 1 2 ∂ j Ũi + ∂ i Ũj , Ωij = 1 2 ∂ j Ũi -∂ i Ũj . (11)
For the near-wall part of the pressure term, φ w ij has to satisfy the balance between the pressure term, the molecular diffusion and the dissipation term, in the vicinity of the wall, which can be written as

φ ij + ν∇ 2 u ′′ i u ′′ j -ε ij = 0. ( 12
)
Using the simple model of Rotta [START_REF] Rotta | Statistische Theorie nicht homogener Turbulenz[END_REF] for the near-wall part of the dissipation term (Eq. 16), Manceau & Hanjalić [START_REF] Manceau | Elliptic Blending Model: a new near-wall Reynolds-stress turbulence closure[END_REF] showed that φ w ij can be written in a general frame as

φ w ij = -5 ε k ′′ u ′′ i u ′′ k n j n k + u ′′ j u ′′ k n i n k - 1 2 u ′′ k u ′′ l n k n l (δ ij + n i n j ) . ( 13
)
Since α is zero at the wall, the wall-normal unit vector n can be evaluated from

n = ∇α ||∇α|| . ( 14
)
This definition enables the use of the model in complex geometries because the n vector continuously adapts itself to the shape of the wall, as shown in Fig. 1, and there is no need to define the nearest wall, which can be ill-defined, as in the step corner, for example. Similarly to φ ij , the standard, homogeneous model for the dissipation rate ε ij is not valid in the near-wall region. Manceau & Hanjalić [START_REF] Manceau | Elliptic Blending Model: a new near-wall Reynolds-stress turbulence closure[END_REF], using a DNS database, have shown that the dissipation tensor can be modelled in a way similar to φ ij

ε ij = (1 -α 2 )ε w ij + α 2 2 3 εδ ij . (15)
Away from the wall, the dissipation rate is assumed isotropic. The simple Rotta model [START_REF] Rotta | Statistische Theorie nicht homogener Turbulenz[END_REF] is used for the near-wall part, as previously mentioned

ε w ij = u ′′ i u ′′ j k ′′ ε. (16)
To take into account the increase of the production of dissipation in the nearwall zone, the coefficient of the generation term in the dissipation equation is taken as

C ′ ε 1 = C ε 1 1 + A 1 (1 -α 2 ) k ′′ u ′′ i u ′′ j n i n j . ( 17
)
This formulation gives the classical value C ε 1 far from the walls, and a value larger than C ε 1 in the near-wall zone. The length and time scales appearing in the model equations are the turbulent scales bounded by the Kolmogorov scales, to avoid singularities at the walls and to replicate the physical behaviour observed in DNS analyses [START_REF] Manceau | Inhomogeneity and anisotropy effects on the redistribution term in Reynolds-Averaged Navier-Stokes modelling[END_REF]

L = C L max k ′′ 3/2 ε ; C η ν 3/4 ε 1/4 , (18)
T = max k ′′ ε ; C T ν ε . (19)
The model constants are given in Tab. 2. Exact wall boundary conditions used for all variables are given in Tab. 1. Comparing to the original model of Manceau & Hanjalić [START_REF] Manceau | Elliptic Blending Model: a new near-wall Reynolds-stress turbulence closure[END_REF], the version of Manceau [START_REF] Manceau | An improved version of the Elliptic Blending Model. Application to non-rotating and rotating channel flows[END_REF] improves the robustness by getting rid of nonlinearities. The two-component limit of turbulence near the wall is well predicted in a channel flow [START_REF] Manceau | An improved version of the Elliptic Blending Model. Application to non-rotating and rotating channel flows[END_REF], and the model has been successfully applied to different configurations [START_REF] Thielen | Predictions of flow and heat transfer in multiple impinging jets with an elliptic-blending second-moment closure[END_REF][START_REF] Borello | Prediction of cascade flows with innovative second-moment closure[END_REF][START_REF] Choi | Predictions of turbulent heat transfer of super-critical carbon dioxide flow in a square duct with an elliptic-blending second moment closure[END_REF]. Another advantage of the model is that there is no explicit dependency on the distance to the wall, and it can therefore be used in complex geometries. Moreover, the derivation of the elliptic relaxation equations ( 6) from the Poisson equation for the fluctuating pressure [START_REF] Manceau | Inhomogeneity and anisotropy effects on the redistribution term in Reynolds-Averaged Navier-Stokes modelling[END_REF] is independent of the steady/unsteady character of the resolved field (velocity Ũi and pressure P), which indicates that the procedure leading to the elliptic blending model can still be applied in an unsteady context. The only possible influence of unsteadiness could be a modification of the near-wall balance of the Reynolds stress equations, from which the near-wall term φ w ij is derived. However, the analysis of the near-wall asymptotic behaviours given by Manceau & Hanjalić [START_REF] Manceau | Elliptic Blending Model: a new near-wall Reynolds-stress turbulence closure[END_REF] shows that the time-derivative is negligible compared to the other terms, meaning that the near-wall modelling is also valid in unsteady flows, contrary to usual methods based on damping functions.

Flow configuration and numerical methods

The 2D backward-facing step flow configuration is presented in Fig. 2. The test case of Driver & Seegmiller [START_REF] Driver | Features of a reattaching turbulent shear layer in divergent channel flow[END_REF] has been chosen for several reasons: the expansion ratio is low enough to ensure that the pressure effects do not overwhelm the turbulence effects [START_REF] Lasher | On the computation of turbulent backstep flow[END_REF], the experimental inlet conditions are well defined, and the aspect ratio is sufficiently large to neglect the effects of the lateral walls. The Reynolds number based on the inlet centerline velocity U 0 and the step height h is Re h = 37500, and the Reynolds number based on the inlet momentum thickness θ and the velocity U 0 is Re θ = 5000. The experiment [START_REF] Driver | Time dependent behavior of a reattaching shear layer[END_REF] underlines the existence of two dominant frequencies. A high frequency is associated with the vortex shedding due to the Kelvin-Helmholtz instability, with a Strouhal number of 0.20, based on the local vorticity thickness and shear velocity (U 0 /2), and a low frequency (St ≈ 0.06), characteristic of a flapping motion. In the present work, the EB-RSM is compared to the standard k-ε [START_REF] Launder | The numerical computation of turbulent flow[END_REF] and LRR (Rotta+IP version) [START_REF] Launder | Progress in the development of a Reynolds-stress turbulence closure[END_REF] models, using wall functions. Computations are performed with Code Saturne, a parallel, finite volume solver on unstructured grids, developed at EDF [START_REF] Archambeau | Code Saturne: A finite volume code for the computation of turbulent incompressible flows -Industrial applications[END_REF]. Space discretization is based on a collocation of all the variables at the centre of gravity of the cells. Velocity/pressure coupling is ensured by the SIMPLEC algorithm, with a Rhie & Chow interpolation in the pressure-correction step [START_REF] Ferziger | Computational methods for fluid dynamics[END_REF]. The Poisson equation is solved with a conjugate gradient method. Time advancement is based on an implicit Euler scheme, and time-step convergence has been checked carefully, by refining the time step up to a factor of 10. The convection terms of all the equations (for resolved momentum and turbulent quantities) are approximated by a second-order central-difference scheme (CDS). Some computations where also performed using a first-order upwind-difference scheme (UDS) for comparison. All the computations shown in the present article are 2D. This is justified by the fact that the purpose is to reproduce the primary instability of the shear layer which is 2D [START_REF] Wee | Self-sustained oscillations and vortex shedding in backward-facing step flows: simulation and linear stability analysis[END_REF]. The reproduction by URANS of the secondary instabilities leading to 3D structures is only possible if the primary instability is found. Six meshes, numbered from 1 to 6, are used for the high-Reynolds number models (k-ε and LRR, see Tab. 3). Six other meshes, numbered from 7 to 12, are also used for the EB-RSM low-Reynolds number model (see Tab. 4). The parameter f r , used in Tab. 3 and Tab. 4, is the factor of refinement in each direction compared to the reference mesh (mesh 1 for the high-Reynolds number models, and mesh 7 for the EB-RSM); N cell is the total number of cells and y + 1 is the distance to the wall, in wall units, of the centre of gravity of the first cell, evaluated at the inlet of the domain. As shown by Manceau & Hanjalić [START_REF] Manceau | Elliptic Blending Model: a new near-wall Reynolds-stress turbulence closure[END_REF], the EB-RSM does not require very fine meshes close to the wall, contrary to most of the other near-wall resolving second-moment closures. Mesh 7 corresponds to the most refined mesh used by Lasher & Taulbee [START_REF] Lasher | On the computation of turbulent backstep flow[END_REF].

A particular attention was paid to the inlet conditions to which the backwardfacing step flow is very sensitive. Since they do not correspond to a fully developed channel flow, separate computations of a developing channel flow were performed with each turbulence model to extract the boundary layer profiles, which are applied at the inlet of the domain in order to match the boundary layer thickness 4h upstream of the step corner observed in the experiments. Fig. 3 shows the mean velocity profile obtained with each turbulence model compared to the experimental data. Fig. 4 and Fig. 5 show respectively the Reynolds stresses and the total fluctuating kinetic energy profile, in comparison with the experiment, where it is assumed that w 2 ≃ 1/2(u 2 + v 2 ), the overbar denoting the long-time average.

Results

Mean flow

In comparison to the standard k-ε and LRR models with wall functions (denoted by WF), Figs. 6 and7 show that the EB-RSM significantly improves the skin friction coefficient C f in a region covering the whole recirculation bubble (0 < x/h < 8). The skin friction coefficients obtained with the EB-RSM, using the different meshes, are superimposed, showing that the results are not sensitive to a reasonable coarsening of the near-wall mesh. This point was also noticed by Manceau & Hanjalić [START_REF] Manceau | Elliptic Blending Model: a new near-wall Reynolds-stress turbulence closure[END_REF] in a channel flow and is very interesting for industrial applications where the meshes are marginally fine. Tab. 5 shows that the mean reattachment length l r is underestimated by 20% by the high-Reynolds models, as expected, and estimated within 5% error by the EB-RSM. However, Fig. 8 shows that the recirculation region is too thick: the shear strain is thus too weak, which implies an underestimation of the shear stress in the recirculation region, as can be seen in the figure. This is a shortcoming common to most of the RANS models (see, for instance, [START_REF] Leschziner | Proc. 10th ERCOFTAC (SIG-15)/IAHR/QNET-CFD Workshop on Refined Turbulence Modelling[END_REF]). The EB-RSM friction coefficient in Fig. 6 and mean velocity profiles in Fig. 8 moreover reveal that the recovery of the boundary layer downstream of reattachment is not satisfactory. This feature is common to all the turbulence models, SMCs as well as eddy-viscosity models [START_REF] Abe | New turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows-1. Flow field calculations[END_REF][START_REF] Craft | Developments in a low-Reynolds-number second-moment closure and its application to separating and reattaching flows[END_REF][START_REF] Durbin | Separated flow computations with the k-ε-v 2 model[END_REF][START_REF] Hanjalić | Contribution towards the second-moment closure modelling of separating turbulent flows[END_REF]. The prediction of the mean flow with the EB-RSM is very similar to previous results obtained with the same class of models (near-wall second-moment closures [START_REF] Craft | Developments in a low-Reynolds-number second-moment closure and its application to separating and reattaching flows[END_REF][START_REF] Hanjalić | Contribution towards the second-moment closure modelling of separating turbulent flows[END_REF]). Experiment k-ε LRR EB-RSM l r /h 6.26 5.14 4.72 6.51

Table 5: Mean reattachment length. Mesh 6 for the high-Reynolds models and mesh 12 for the EB-RSM.

Instantaneous flow

All the k-ε computations give steady solutions, whatever the mesh or the spatial discretization used. This model is based on the Boussinesq approximation which assumes alignment of the Reynolds stress anisotropy and the strain. As a consequence, the turbulent production is overestimated. For instance, in a 2D flow, as shown by Carpy & Manceau [START_REF] Carpy | Turbulence modelling of statistically periodic flows: Synthetic jet into quiescent air[END_REF], the exact production can be written as

P = k ′′ β(λ 1 -λ 2 ) cos(2ϕ), ( 20
)
where ϕ is the angle between the eigenvectors of the strain rate tensor and the eigenvectors of the anisotropy tensor in the x-y-plane, β is the positive eigenvalue of the strain rate tensor and λ 1 and λ 2 are the two eigenvalues of the anisotropy tensor in the x-y-plane (λ 1 > λ 2 ). The Boussinesq approximation implies that ϕ = 0, leading to a systematic overestimation of production when turbulence and strain are not aligned, which is the case if the strain is unsteady. This leads to an overestimated rate of transfer from resolved to modelled energy, leading to a rapid decrease of the former [START_REF] Carpy | Turbulence modelling of statistically periodic flows: Synthetic jet into quiescent air[END_REF]. Therefore, linear eddy-viscosity models should be avoided in unsteady computations.

On the contrary, the steady/unsteady character of the solution using the SMCs is strongly dependent on the numerics. Using an UDS, whatever the meshes and the model, the solution is steady due to the diffusive character of this scheme. Using a CDS and meshes similar to those used by Lasher & Taulbee [START_REF] Lasher | On the computation of turbulent backstep flow[END_REF] (mesh 7 for the EB-RSM), the solution is unsteady, exhibiting a single frequency, consistent with the results found by these authors. The correct order of magnitude of the Strouhal number is found with EB-RSM (St = 0.21) whereas St = 0.16 with LRR. This frequency is associated with the vortex shedding which can be visualized in Fig. 9, by isocontours of vertical resolved velocity fluctuations. Seven probes are positioned in the shear layer (y/h = 1) going from x/h = 1 to x/h = 7. The resolved velocity signal as a function of time, not shown here, is purely sinusoidal. Therefore, the measure of the convection velocity U c of the large-scale structures is easy and is given by the time delay between these probes. The measure gives a constant value U c ≃ 0.5U 0 independent of the streamwise location, as reported in the experiment [START_REF] Kiya | Structure of a turbulent seperation bubble[END_REF]. However, in the experiments, it is observed that some of the vortices created in the shear layer are incorporated into the recirculation zone rather than convected downstream, which is not observed in our computations: consequently, the flapping motion, which is linked to the incorporation and release of this vortices, is missed. However, the resolved, large-scale structures contain a very limited part of the total fluctuating energy: indeed, the ratio k res /k tot is less than 5% on the coarse mesh (mesh 7), where k res is the kinetic energy of the resolved structures and k tot = k res + k mod is the sum of the resolved and the modelled kinetic energies. These quantities are given by

k res = 1 2 u ′ i u ′ i = 1 2 Ũi Ũi -U i U i , (21)
k mod = k ′′ = 1 2 u ′′ i u ′′ i . ( 22
)
One could expect that refining the mesh would intensify the resolved energy, by reducing the numerical errors. Let us define the parameter M (x) by

M (x) = max y k res (x, y)/k tot (x, y) , ( 23
)
which is the maximum value of the ratio resolved energy/total energy on a vertical line at a given streamwise location. This ratio characterizes the energy contained in the resolved, large-scale eddies of the shear layer, coming from the vortex shedding. Fig. 10 shows the streamwise evolution of M (x) for each mesh used in the EB-RSM computation: a drastic reduction of the resolved energy and a steady solution for the finest mesh (mesh 12), meaning that M (x) = 0, are obtained. It is worth pointing out that this behaviour is not linked to the fact that the computations are 2D. Indeed, the primary instability of this flow being 2D [START_REF] Wee | Self-sustained oscillations and vortex shedding in backward-facing step flows: simulation and linear stability analysis[END_REF], 3D structures generated by secondary instabilities cannot be reproduces if the primary instability is missed. Therefore, 3D computations give exactly the same steady solution as 2D computations, which has been checked by using the same mesh [START_REF] Durbin | A Reynolds stress model for near-wall turbulence[END_REF] extended in the z-direction by adding 40 cell layers covering a width of 4h.

The results with the LRR model are qualitatively similar, and are therefore not shown here. The next section proposes an explanation for this surprising behaviour.

Investigation of the steady/unsteady behaviour

First, it is worth emphasizing that, contrary to LES, the URANS model is independent of the grid. Refining the grid enables the solution of the discretized equations to approach the solution of the continuous equations (the model), while in LES, where the filter is usually dependent on the local characteristics of the mesh, grid convergence cannot be obtained (except in the DNS limit).

Therefore, in URANS, a modification of the solution when refining the grid can only be due to discretization errors. Actually, the nature of the solution (steady/unsteady) can be traced to the excitation of the most unstable mode of the shear layer by weak spatial oscillations at the step corner, due to the CDS. Indeed, it is well known, at least for a 1D convection/diffusion equation, that central differencing can generate spatial numerical oscillations in a region of strong gradients in the streamwise direction, if the Peclet number exceeds 2 [START_REF] Ferziger | Computational methods for fluid dynamics[END_REF]. These oscillations do not appear when using a diffusive scheme, such as the UDS, which explains why the solution always remains steady with this scheme. In a one-dimensional situation, the Peclet number compares the time scales of diffusion τ D and convection τ C over the length of the cell. In the two-dimensional boundary layer just upstream the step corner, the convection and diffusion process can be considered as essentially in the streamwise and wall-normal directions, respectively, and associated time scales can be locally approximated by τ C = ∆x/U and τ D = ∆y 2 /ν eff . Thus, a local Peclet number can be defined as

P e(x, y) = τ D τ C = U ∆y 2 ν eff ∆x . (24)
In a turbulent boundary layer, ν eff accounts for the effective (molecular + turbulent) viscosity and is given by the Daly-Harlow model [START_REF] Daly | Transport equations in turbulence[END_REF], used in both SMCs:

ν eff = ν + C µ σ k T v ′′ 2 . ( 25
)
The wall-normal distance δ influenced by the total diffusion during the transport time τ C is equal to δ = ν eff ∆x/U , meaning that δ/∆y = P e -1/2 . As in a one-dimensional situation, we conjecture that for large values of the local Peclet number, δ is much smaller than the wall-normal cell size ∆y, which leads to spatial oscillations. Fig. 11 shows the streamwise evolution of the local Peclet number at y + = 15, which is the location of the streamline passing through the inflexion point of the shear layer just after the step, where the velocity profile is the most sensitive to perturbations, according to linear stability [START_REF] Wee | Self-sustained oscillations and vortex shedding in backward-facing step flows: simulation and linear stability analysis[END_REF]. It can be seen, comparing Fig. 10 and Fig. 11, that there is at least a qualitative link between the orders of magnitude of the local Peclet number and the energy contained in the resolved eddies. This comparison suggests that, for coarse meshes, the spatial numerical oscillations due to the high local Peclet number in a region of strong mean velocity gradient ∂ x U (very close to the step, say -0.1 < x/h < 0) act as perturbations that are sufficient to excite the natural mode of the shear layer. Fig. 12 attempts to highlight the existence of these very small numerical oscillations on the mean velocity. Since they are not directly visible on the mean velocity itself, the streamwise evolution of the long-time averaged velocity derivative ∂ x U is plotted instead, at the height y + = 15. It is seen that numerical oscillations exist and their amplitude decreases when the mesh is refined. These oscillations are not the footprint of the unsteadiness due to the vortex shedding, since what is plotted is the long-time averaged velocity. The steady solution was not observed by Lasher & Taulbee [START_REF] Lasher | On the computation of turbulent backstep flow[END_REF] because their most refined mesh corresponds to the coarsest mesh used in the present work (mesh 7). It is also important to note the existence of an unsteady solution, even on the most refined mesh (mesh 12), during a transient phase which can last for a long time, about twenty shedding time scales. Thus, computations must be performed on a sufficiently large duration to reach the steady solution.

When the solution is unsteady on the coarse meshes, the Strouhal number of the large-scale structures is reliable because it does not depend much on the mesh, since it corresponds to the most amplified mode of the shear layer, as predicted by linear stability [START_REF] Wee | Self-sustained oscillations and vortex shedding in backward-facing step flows: simulation and linear stability analysis[END_REF]. On the contrary, the amplitude of these structures is strongly dependent on the amplitude of the perturbations, which is driven by the numerical oscillations. Thus, obtaining an unsteady solution is a numerical artefact, which makes the use of URANS uncertain in such flows, in particular as concerns the prediction of the amplitude of the large-scale oscillations, which are of major interest for industrial applications.

Another important remark can be made by considering two meshes, denoted by the subscripts a and b. It is assumed that mesh b is finer than mesh a, with a factor of refinement in streamwise and wall normal directions denoted by f ab ≥ 1 and g ab ≥ 1 respectively. If the mean velocity and the turbulent viscosity are converged numerically, using definition [START_REF] Leschziner | Proc. 10th ERCOFTAC (SIG-15)/IAHR/QNET-CFD Workshop on Refined Turbulence Modelling[END_REF], it is easy to show that the ratio of the Peclet numbers associated with each mesh is

P e b P e a = ∆y b ∆y a 2 ∆x a ∆x b = 1 g ab 2 f ab . (26)
Considering, for example, that meshes a and b are mesh 7 and 12 (EB-RSM computation), we have in this case f ab = g ab = f r = 2 (see Tab. 4) and then P e 12 /P e 7 = 1/2, which is approximately verified in Fig. 11. Expression [START_REF] Manceau | Elliptic Blending Model: a new near-wall Reynolds-stress turbulence closure[END_REF] suggests that it is more favourable to numerical stability to refine only in the wall normal direction than in both directions, leading to elongated cells.

Conclusions

The 2D computations of the backward-facing-step flow have shown that obtaining an unsteady solution, whatever the model used (high or low-Reynolds number model), is due to the amplification of numerical oscillations in a region very close to the step, where the velocity gradients and local Peclet number are high. These numerical oscillations act as perturbations that are sufficient to excite the natural mode of the shear layer. A diffusive scheme, such as the UDS, do not exhibit such oscillations and, therefore, gives a steady solution whatever the mesh. The present work thus shows that the grid-converged solution of 2D URANS computations of the backward-facing step flow is a steady (RANS) solution. This indicates that the primary instability of this flow is not reproduced by URANS. Consequently, more complex structures, like those due to secondary instabilities cannot be obtained either, and 3D URANS computations also give a steady solution. The observations made using CDS and UDS schemes can be extended to more sophisticated schemes (upwind-biased second-order schemes, TVD schemes, etc.): the appearance of unsteady solutions can only be a numerical artefact due to the combination of a too coarse mesh and an oscillation-generating scheme.

As shown by previous studies [START_REF] Lasher | On the computation of turbulent backstep flow[END_REF], URANS is able to reproduce the correct Strouhal number of the vortex shedding, since it is imposed by the mean velocity profile, according to the linear stability [START_REF] Wee | Self-sustained oscillations and vortex shedding in backward-facing step flows: simulation and linear stability analysis[END_REF]. However, the amplitude of the vortices is strongly dependent on the numerics, and therefore, its prediction with the URANS methodology does not seem to be reliable. This result suggests that, whatever the convection scheme, in industrial applications using URANS, the reproduction of vortex shedding in separated shear layers will be driven by the competition between numerical diffusion and dispersion. This issue can be problematic in complex 3D applications for which the mesh is often only marginally fine. The interpretation of the results in terms of, e.g., temporal variations of forces or wall temperatures, is to be considered with care. This remark suggests that future work should be devoted to the development of hybrid RANS-LES methods, which enable to control the parameter resolved energy/total energy. The present work also showed that the Elliptic Blending Reynolds Stress Model [START_REF] Manceau | Elliptic Blending Model: a new near-wall Reynolds-stress turbulence closure[END_REF][START_REF] Manceau | An improved version of the Elliptic Blending Model. Application to non-rotating and rotating channel flows[END_REF], a recent near-wall second-moment closure, is sufficiently robust to be used in unsteady computations, although the grid-converged solution is steady in the present case. The steady results are very similar to previous results obtained with near-wall second-moment closures [START_REF] Craft | Developments in a low-Reynolds-number second-moment closure and its application to separating and reattaching flows[END_REF][START_REF] Hanjalić | Contribution towards the second-moment closure modelling of separating turbulent flows[END_REF], both in the recirculation and recovery regions. The mean skin friction in the recirculation region and the mean reattachment length are improved compared to the standard k-ε [START_REF] Launder | The numerical computation of turbulent flow[END_REF] or LRR [START_REF] Launder | Progress in the development of a Reynolds-stress turbulence closure[END_REF] models with wall-functions, and are not sensitive to a reasonable coarsening of the near-wall mesh. On the contrary, the reproduction of the recovery region after the reattachment point is not satisfactory. The EB-RSM is much more robust than previous models based on Elliptic Relaxation [START_REF] Durbin | A Reynolds stress model for near-wall turbulence[END_REF][START_REF] Manceau | Elliptic Blending Model: a new near-wall Reynolds-stress turbulence closure[END_REF] and can be easily implemented in RANS codes using second-moment closures.

Figure 2 :

 2 Figure 2: Geometry and computational domain of the backward-facing step flow.

Figure 3 :

 3 Figure 3: Mean velocity profile at the inlet (x/h = -4). Comparison between models and experiment [9].

Figure 4 :

 4 Figure 4: Reynolds stress profiles at the inlet (x/h = -4). Comparison between models and experiment [9]. • u 2 , 2 v 2 , ⋄ uv.

Figure 5 : 6 Figure 6 :Figure 7 :Figure 8 :

 56678 Figure5: Fluctuating kinetic energy profile at the inlet (x/h = -4). Comparison between models and experiment[START_REF] Driver | Features of a reattaching turbulent shear layer in divergent channel flow[END_REF].

Figure 9 :Figure 10 :Figure 11 :Figure 12 :

 9101112 Figure 9: Isocontours of vertical, resolved velocity fluctuations (mesh 10, EB-RSM computation). Dashed lines are negative values.

Table 2 :

 2 Values of the EB-RSM constants.

Table 3 :

 3 Mesh characteristics for the high-Reynolds number models (k-ε and LRR).

	Comparison between

Table 4 :

 4 Mesh characteristics for the low-Reynolds number model (EB-RSM).