Xavier Gloerfelt

Paola Cinnella
email: paola.cinnella@ensam.eu

Large Eddy Simulation Requirements for the Flow over Periodic Hills

Keywords: Large eddy simulation, High-accuracy method, Periodic hill flow

Large eddy simulations are carried out for flows in a channel with streamwise-periodic constrictions, a well-documented benchmark case to study turbulent flow separation from a curved surface. Resolution criteria such as wall units are restricted to attached flows and enhanced criteria, such as energy spectra or two-point correlations, are used to evaluate the effective scale separation in the present large eddy simulations. A detailed analysis of the separation above the hill crest and of the early shear layer development shows that the delicate flow details in this region may be hardly resolved on coarse grids already at Re = 10 595, possibly leading to a non monotonic convergence with mesh refinement. The intricate coupling between numerical and modeling errors is studied by means of various discretization schemes and subgrid models. It is shown that numerical schemes maximizing the resolution capabilities are a key ingredient for obtaining high-quality solutions while using a reduced number of grid points. On this respect, the introduction of a sharp enough filter is an essential condition for separating accurately the resolved scales from the subfilter scales and for removing ill-resolved structures. The high-resolution approach is seen to provide solutions in very good overall agreement with the available experimental data for a range of Reynolds numbers (up to 37 000) without need for significant grid refinement.

Introduction

In large eddy simulations (LES), only the dynamics of the large scales is computed and the effects of smaller scales are modeled. Scale separation is however difficult to establish since the low-pass filtering arises from a complex combination of implicit filtering by the grid and the discretization schemes. Even when explicit filters are applied, the approximations introduced by the discretization methods modify the actual shape of the filter function. The question of the intricate interactions between numerical errors, leading to this badly-defined filter, and subgrid-scale (SGS) modeling errors has retained the attention of numerous authors since the studies of Ghosal [START_REF] Ghosal | An analysis of numerical errors in large-Eddy Simulations of turbulence[END_REF], Vreman et al. [START_REF] Vreman | Comparison of numerical schemes in large-eddy simulation of the temporal mixing layer[END_REF] or Kravchenko and Moin [START_REF] Kravchenko | On the effect of numerical errors in large eddy simulations of turbulent flows[END_REF]. Even using high-order schemes the subgrid contribution can be dominated by numerical errors in the high-wavenumber range. Coupling of numerical and modeling errors may lead to a non-monotonic dependency on grid spacing, making comparisons with DNS (Direct Numerical Simulation) or experimental results difficult. For instance, Meyers et al. [START_REF] Meyers | Optimal model parameters for multi-objective large-eddy simulations[END_REF] noted that better results for wall shear stress in a turbulent channel flow can be obtained by coarsening the grid. In [START_REF] Meyers | Optimal model parameters for multi-objective large-eddy simulations[END_REF][START_REF] Meyers | A computational error-assessment of central finite-volume discretizations in large-eddy simulation using a Smagorinsky model[END_REF] error-landscapes are used to provide a more systematic assessment of computational errors and their coupling with SGS modeling for homogeneous isotropic turbulence or channel flow. It is however difficult to generalize this kind of analysis to complex flows. The interactions are even more complex for compressible flow solvers, which always require some form of numerical dissipation or filtering to ensure the stability of the simulation. Well-established resolution criteria for free flows or attached wall flows exist, such as grid sizes expressed in wall units. However, these criteria do not account for the numerical scheme resolution capabilities on a given grid. For instance, a low-accurate scheme may lead to more restrictive grid requirements in order to rule out numerical errors; conversely, high-resolution schemes may allow for less severe grid refinement criteria, with respect to the recommended ones. A more in-depth level of validation, such as the analysis of flow unsteadiness, is sometimes required to have a better understanding of where and why resolution criteria are not satisfied. Some authors, such as Davidson [START_REF] Davidson | Large eddy simulations: how to evaluate resolution[END_REF][START_REF] Davidson | How to estimate resolution of an LES of recirculating flow[END_REF], have proposed guidelines to quantify the quality of LES independently, notably in presence of separated flows. Despite this, most LES studies still focus only on averaged quantities, such as mean velocity profiles or turbulent intensities, to judge of the quality of the simulations by comparison with some reference experimental or numerical data.

In the present study, the configuration of the flow in a channel with streamwise-periodic constrictions is used to gain insight on the LES requirements for the prediction of separated turbulent flows. This particularly challenging configuration has become a popular benchmark case for LES in the past years, with experimental and numerical references available as an ERCOFTAC database [START_REF] Rapp | 2d periodic hill flow[END_REF]. The periodic hill flow has been widely used in the past to assess a large variety of numerical methods, meshing strategies, subgrid-scale (SGS) models, and wall modeling assumptions (based either on RANS/LES coupling or on wall functions) [START_REF] Balakumar | DNS, enstrophy balance, and the dissipation equation in a separated turbulent channel flow[END_REF][START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF][START_REF] Davidson | Hybrid LES-RANS modelling: a one-equation SGS model combined with a k-ω model for predicting recirculating flows[END_REF][START_REF] Duprat | A wall-layer model for large-eddy simulations of turbulent flows with/out pressure gradient[END_REF][START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF][START_REF] Fröhlich | Hybrid LES-RANS methods for the simulation of turbulent flows[END_REF][START_REF] Temmerman | A hybrid two-layer URANS-LES approach for large eddy simulation at high Reynolds numbers[END_REF][START_REF] Temmerman | Investigation of wall-function approximations and subgrid-scale models in large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF][START_REF] Tessicini | Approximate near-wall treatments based on zonal and hybrid RANS-LES methods for LES at high Reynolds numbers[END_REF][START_REF] Xia | Constrained large-eddy simulation of separated flow in a channel with streamwise-periodic constrictions[END_REF][START_REF] Ziefle | Large-eddy simulation of separated flow in a channel with streamwiseperiodic constrictions[END_REF]. The level of validation is most of the time restricted to mean velocity profiles, skin friction and turbulent intensities. More in-depth flow analyses, such as those of Fröhlich et al. [START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF] or Breuer et al. [START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF], have however clearly shown that the flow exhibits a complex dynamics, involving separation and reattachment from a curved boundary, generation of an unsteady shear layer surrounding a large recirculation bubble, strong acceleration on the convex windward wall of the next hill, attached boundary layer on the ceiling subject to intense fluctuating pressure gradients, and the turbulence recycling due to the periodicity assumption. Previous studies have demonstrated that the boundary layer flow on the bottom wall is always far from equilibrium, which is very challenging for wall functions or RANS modeling of near-wall effects. Additionally, the flow is characterized by significant irregularity and intermittency, which makes the identification of reproducible coherent structures driving the flow dynamics a very diffult task. Solution quality often depends on the mean locations of the separation and reattachment points, which are highly unsteady due to the strong chaotic motions of the separation and reattachment lines. Thus, examining the mean behavior is not sufficient to discriminate about the quality of a numerical or modeling strategy.

The paper aims at contributing to the analysis of the numerical and physical ingredients that are required to realize LES of such a separated flow. After a brief presentation of the problem setup, LES of periodic hill flow at Re = 10 595 are carried out using a high-resolution numerical method, relying on an optimized finite difference scheme supplemented with a sharp-cutoff selective filter. The latter is tailored to ensure both numerical damping of ill-resolved solution modes and the necessary regularization of subgrid flow scales. Such a numerical approach was selected based on a previous study, published in this journal [START_REF] Aubard | Comparison of subgrid-scale viscosity models and selective filtering strategy for large-eddy simulations[END_REF], that was however conducted for simpler configurations (namely, the Taylor-Green-vortex and an attached turbulent boundary layer). For such flows, it was shown that high-resolution schemes and sharp cutoff filters are mandatory to obtain a clear separation between resolved and unresolved scales. It is not clear however to which extent those findings can be directly extrapolated to a complex separated flow configuration such as the 2D-hill, and which criteria are the more suitable to characterize the numerical resolution. A series of increasingly fine grids is used to appraise several resolution criteria, including resolution in wall units, energy spectra and cross-correlations, and to identify sensitive regions where some criteria are difficult to fulfill. The ability to capture the underlying flow dynamics is then analyzed to explain discrepancies with the admitted reference values. Specifically, small compressibility effects characterizing the high-speed region at the hill crest and, most importantly, the early birth of the shear layer are shown to have a critical influence on the flow dynamics. Afterwards, the interactions between the numerical discretization and the subgrid model and their influence on the resolved flow features are investigated. For that purpose, several centered finite difference schemes and selective filters (both standard and optimized) are considered, which allows a separate study of the influence of dispersion and dissipation errors on LES quality. Special attention is paid to the choice of the filtering term, which is identified as the most influential numerical ingredient because of its interaction with SGS modeling. Filter-based SGS regularization is compared with explicit SGS models, based on classical eddy viscosity concept, to highlight the relative influence of small-scale modeling and numerical scheme resolution. Higher Reynolds number cases (Re = 19 000 and Re = 37 000) are also presented in the last part, with the double intent of: i) further supporting the important role of high-resolution schemes in capturing the flow dynamics and ii) demonstrating the ability of the proposed LES strategy to capture delicate Reynolds number effects despite the use of rather coarse grids. Even if the analyses are restricted to a particular family of schemes, the repeatability of the outcomes of [START_REF] Aubard | Comparison of subgrid-scale viscosity models and selective filtering strategy for large-eddy simulations[END_REF] for the present complex flow configuration provides interesting guidelines for the high-resolution schemes in general.

Problem Setup

Governing equations

The governing equations are the compressible Navier-Stokes equations written for a curvilinear domain by using a coordinate transform. The physical space (x, y) is mapped into a Cartesian regular computational space (ξ, η), and the third direction z is left unchanged. The set of equations for the unknown vector U = (ρ, ρu, ρv, ρw, ρE) T writes:

∂U ∂t + ∂F c ∂ξ + ∂G c ∂η + ∂H ∂z = f (1)
by defining the curvilinear fluxes as:

F c = 1 J F ∂y ∂η -G ∂x ∂η and G c = 1 J -F ∂y ∂ξ + G 1 J ∂x ∂ξ
where J is the Jacobian of the coordinate transform and F, G and H are the sum of the inviscid and visco-thermal fluxes (see [START_REF] Chicheportiche | Study of interpolation methods for high-accuracy computations on overlapping grids[END_REF][START_REF] Garrec | Multi-size-mesh multi-time-step algorithm for noise computation on curvilinear meshes[END_REF]). The vector f at the right-hand side of Eq. 1 is a forcing function. The specific total energy is defined as

E = p/[(γ -1)ρ] + (u 2 + v 2 + w 2 )/2
for an ideal gas. The viscous stress tensor is modeled for a Newtonian fluid as

σ ij = 2μS ij -(2/3)μS kk δ ij
, where S ij = (u i,j + u j,i )/2 is the strain rate tensor and μ is the dynamic viscosity, computed with Sutherland's law, and may be supplemented by the subgrid Reynolds stress tensor τ ij in case an explicit subgrid model is used (see Section 3.2).

Flow configuration

Sideviews of the streamwise-periodic-hill channel flow configuration are provided in Fig. 1.

The configuration consists in a plane channel with two-dimensional curved constrictions of height h. The exact shape functions of the hill can be found in the ERCOFTAC website [START_REF] Rapp | 2d periodic hill flow[END_REF].

The channel height is L y = 3.035h. A periodicity is imposed in the streamwise direction x with a sufficiently large distance L x = 9h between two consecutive hills in order to allow a partial recovery of the flow after its reattachment and an acceleration on the windward wall of the next hill. The third direction z is assumed homogeneous with L z = 4.5h, meaning that the span is very large and the flow can be considered as 2-D at mid-span. The xand zperiodicities are favorable insofar as the specification of a turbulent inlet condition is circumvented and the averaging along the span can be used to improve the statistical convergence.

Well-resolved incompressible LES have been published to serve as a reference for validation purposes [START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF][START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF]. Furthermore, proper experiments have been designed to reproduce this flow configuration using 10 hills in a water tunnel to mimic periodicity in the streamwise direction. The setup and the measurement results are documented in [START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF][START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF] and the data are made available as an ERCOFTAC database [START_REF] Rapp | 2d periodic hill flow[END_REF]. Convincing comparisons between different codes and between experiments and simulations have given confidence that this benchmark is well-posed. Its success in the CFD community is proven by the choice of this benchmark case in several European projects, e.g. DESIDER [START_REF]DESider -A European Effort on Hybrid RANS-LES Modelling; Results of the european-Union Funded Project[END_REF], ATAAC [START_REF] Schwamborn | Results and lessons learned from the eu-project Ataac[END_REF] and IDIHOM [START_REF] Kroll | IDIHOM: Industrialization of High-Order Methods -A Top-Down Approach; Results of a Collaborative Research Project Funded by the European Union[END_REF], to investigate the reliability of RANS/LES strategies or the efficiency of high-order codes in industrial-like configurations, respectively. The flow was also used as a test case in three ERCOFTAC/IAHR/COST workshops on high-order methods for compressible flows [START_REF]st International Workshop on High-Order CFD Methods[END_REF][START_REF]2nd International Workshop on High-Order CFD Methods, sponsored by DLR, AIAA and AFOSR[END_REF][START_REF]3rd International Workshop on High-Order CFD Methods, sponsored by NASA, AIAA[END_REF]. A Reynolds number is built from the hill height h, the kinematic viscosity of the fluid ν, and the bulk velocity at the top of the hill U b , defined by:

U b = 1 2.035h 3.035h h u (0, y) dy
where u is the mean streamwise velocity averaged over time and the spanwise direction z. Moderate values of Re = U b h/ν are selected to achieve a fully turbulent state while maintening the number of grid points required for a well-resolved LES computationally tractable. Many studies availabled in the literature have been conducted for Re = 10 595, a value originally used by Mellen et al. [START_REF] Mellen | Large eddy simulation of the flow over periodic hills[END_REF]. Validation experiments [START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF] were conducted in a water channel, so that an incompressible regime was considered in the reference numerical studies [START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF][START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF]. Since our flow solver is compressible, we imposed a low value for the Mach number, M=0.2 based on average conditions in the restricted section. Since periodic conditions are imposed in the streamwise direction, a forcing function f = (0, f u , 0, 0, f u u) T is added to the right-hand side of the Eq. 1 to enforce a constant mass-flow rate. Following [START_REF] Benocci | The role of the forcing-term in large-eddy simulation of equilibrium channel flow[END_REF][START_REF] Lenormand | Large eddy simulation of subsonic and supersonic channel flow at moderate Reynolds number[END_REF] the update of the term f u is defined as:

f n+1 u = f n u + L x V 2 Δt Q n+1 -Q 0 - 0.2 Δt Q n -Q 0 with Q n+1 = Q n - Δt L x Vf n u + D n + α Q n eff -Q 0
where L x is the streamwise length, V the total volume, Δt the timestep between iterations n and n + 1, Q the mass-flow rate and D the drag force exerted on the upper and lower walls. Q 0 is the target mass-flow rate, defined as For the present simulations we use a set of four grids of increasing density (with a refinement factor of about 2 in each space directions), with a number of points ranging from 64 × 33 × 32 to 512 × 256 × 256. The grids are adapted from those prescribed for the 3rd High-Order Workshop [START_REF]3rd International Workshop on High-Order CFD Methods, sponsored by NASA, AIAA[END_REF]. A view of the coarsest grid is presented in Fig. 2. Finer grids have the same topology. Parallel execution is achieved by cutting the computational meshes into block of 32 3 using MPI libraries. A number of computational cores included between 2 and 1024 is used, according to the grid density. For all cases, a time step corresponding to a CFL number of about 0.5 is used. After 23 flow-through times L x /U b , mean quantities are collected over a time period of 55 to 80 flow-through times, and are also averaged in the spanwise direction.

ρ ref U b × 2.035hL z .

LES Methodology

Numerical schemes

The numerical solver is the same as the one used in [START_REF] Aubard | Comparison of subgrid-scale viscosity models and selective filtering strategy for large-eddy simulations[END_REF], to which we refer for the formulation of the filtered equations and of the subgrid-scale terms. The governing (1) are integrated in time using an explicit low-storage six-step Runge-Kutta scheme optimized in the wavenumber space [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for noise computation[END_REF].

The spatial discretization is made up of two elements: a centered finite-difference (FD) scheme, which introduces a purely dispersive error, and a centered selective filter (SF), which is purely dissipative. A spatial derivative of a quantity u is approximated with centered finite differences as:

∂u ∂x j = 1 Δx N l=-N a l u ((j + l)Δx) (2)
The coefficients a j for standard schemes on a stencil of 2N + 1 points are obtained by cancelling terms up to order 2N in Taylor's expansions (hereafter referred to as o2 to o10 for stencils using 3 to 11 points). Some coefficients a j can be used to minimize the dispersion error in the wavenumber space. The Dispersion-Relation-Preserving (DRP) optimization used in the present study is described in [START_REF] Garrec | Multi-size-mesh multi-time-step algorithm for noise computation on curvilinear meshes[END_REF][START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for noise computation[END_REF]. Our baseline scheme will be an elevenpoint-stencil optimized finite-difference scheme [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for noise computation[END_REF], referred to as DRP11. The choice of an 11-point stencil has been proven to be a good trade-off in term of efficiency for wave propagation phenomena [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for noise computation[END_REF], and for large eddy simulations of developed turbulence [START_REF] Berland | A study of differentiation errors in large-eddy simulations based on the EDQNM theory[END_REF].

In the Appendix (see also [START_REF] Gloerfelt | Investigation of high-order methods in large-eddy simulation of separated flow in a channel with periodic constrictions[END_REF]), standard finite-difference schemes of various order for the convective terms are compared to the DRP11 scheme in terms of accuracy and computational cost to show the effect of numerical resolvability (intended as the minimal number of grid points required to achieve a given accuracy) on the quality of the LES simulations. Dispersion errors in the wavenumber space are depicted in Fig. 3a for various schemes of the family (2). Standard fourth-order finite differences are used in all cases for the viscous and heat fluxes.

A centered selective filter is used to eliminate grid-to-grid oscillations originating from unresolved solution modes and perturbations due to, e.g., boundary conditions or very stiff gradients. These must be eliminated to prevent the divergence of the numerical procedure. A filtered quantity u f j is computed as : ; fourth-order (o4) ; sixth-order (o6)

u f j = u j -χD f j with D f j = N l=-N d l u ((j + l)Δx) (3)
; eighth-order (o8) ; tenth-order (o10) and optimized 11-point (DRP11) with the coefficient χ chosen between 0 and 1. The filter has symmetric coefficients d l , which can be determined by cancelling terms up to order 2N in Taylor's expansions (hereafter referred to as o2 to o10 for stencils using 3 to 11 points). Some coefficients can also be used to minimize the dissipation error in the wavenumber space. The baseline filter considered in the following is an optimized filter with an eleven-point stencil [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for noise computation[END_REF] (referred to as the DRP11 filter). The dissipation errors introduced by the different filters tested in the present study are presented in Fig. 3b. Selective filtering is applied in each direction at the end of each time step with a constant value of the filtering coefficient χ . If not otherwise specified, a value of 0.2 is used in the following, which is the minimal value to prevent divergence. The influence of this parameter is discussed in Section 6.1.

At solid boundaries, the no-slip condition u = v = w = 0 is imposed, along with ∂p/∂n = 0 for the Eulerian part (where n is the direction normal to the solid surface), and an isothermal wall condition for the temperature. The finite-difference stencil for the convective terms is progressively reduced down to the second order when approaching solid walls, while viscous stress terms are evaluated from the interior points by using fourth-order backward differences. When the DRP11 filter is used, non-centered 11-point-stencil filters are used near the walls. The values of the coefficients of the non-centered filters can be found in Berland et al. [START_REF] Berland | High-order, low dispersive and low dissipative explicit schemes for multiple-scale and boundary problems[END_REF]. When standard filters are used, the order is reduced down to the second order as approaching the walls.

Subgrid modeling

In the following discussion, we restrict the attention to four classical SGS models already used in [START_REF] Aubard | Comparison of subgrid-scale viscosity models and selective filtering strategy for large-eddy simulations[END_REF]. A brief description of each model and some implementation issues are given hereafter.

Smagorinsky Model (SM)

In the Smagorinsky model, the deviatoric part of the subgridscale (SGS) tensor is modeled as

τ d ij = τ ij -1/3τ kk δ ij = -2ρν t S ij
, where S ij is the resolved strain rate tensor and ν t is the turbulent eddy viscosity, defined as:

ν t = C 2 S Δ 2 c (2S ij S ij ) 1/2 (4)
where Δ c = (ΔxΔyΔz) 1/3 and C S is the Smagorinsky constant. Moreover, the isotropic part of the SGS tensor τ ij is modeled, according to Yoshizawa [START_REF] Yoshizawa | Statistical theory for compressible turbulent shear flows, with the application to subgrid modeling[END_REF], by τ kk = 2C I ρΔ 2 c (2S ij S ij ) 1/2 , where C I = 0.066 [START_REF] Erlebacher | Toward the large-eddy simulation of compressible turbulent flows[END_REF]. For wall-bounded flows, the characteristic length scale Δ c is weighted near the walls using the classical formulation of the van Driest damping function [START_REF] Moin | Numerical investigation of turbulent channel flow[END_REF], Δ c = {1exp(y + /A + )}Δ c , where A + = 25 and y + = yu τ /ν is the wall-normal coordinate. The latter is chosen because of its large spread in the literature, even if more advanced formulations exist (e.g., [START_REF] Piomelli | Model consistency in large eddy simulation of turbulent channel flows[END_REF]), which better capture the asymptotic behavior of ν t toward the wall.

Dynamic Smagorinsky Model (DSM)

In the dynamic version, a test filter is introduced to determine a local value of the Smagorinsky constant, noted C d :

ν t = C 2 d Δ 2 c (2S ij S ij ) 1/2
The calculation of C d is based on Germano's identity [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF] and the least-square procedure of Lilly [START_REF] Lilly | A proposed modification of the Germano subgrid-scale closure method[END_REF]. The same method is applied to have a dynamic estimation of C I for the isotropic part. In the numerical implementation the test filter is a Laplacian. The positiveness of the total viscosity ν+ν t is enforced and the dynamic constants are filtered with a Laplacian filter to avoid too fast variations. The presence of solid walls, laminar flow regions or backscatter effects is taken into account thanks to the dynamic adaptation of ν t .

Multiscale Smagorinsky Model (MSM)

The multiscale model has been proposed by Hughes et al. [START_REF] Hughes | The multiscale formulation of large eddy simulation : decay of homogeneous isotropic turbulence[END_REF], where it is referred to as variational multiscale method. Since the effect of SGS on the resolved scales is limited to interactions with smaller scales, a separation ab initio into three types of scales (large, small and unresolved) is realized by using a test filter. The SGS eddy viscosity model is computed as a function of the smallest resolved scales. In the present study, the test filter is an 11-point filter with a cutoff at kΔx = π/3 [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for noise computation[END_REF]. Two versions of the multiscale method are considered, both based on the SM model. In the small-small version (MSM-ss), both the SGS tensor and the eddy viscosity are expressed as functions of the test-filtered rate of strain tensor Ŝij (obtained by applying the test filter to the resolved strain tensor S ij ):

τ d ij = -2ρ νt Ŝij and ν = C 2 S Δ 2 c (2 Ŝij Ŝij ) 1/2
In the large-small version (MSM-ls), the eddy viscosity is computed as in the standard model (4), while the SGS tensor is written as a function of the test-filtered rate of strain tensor, i.e.

τ d ij = -2ρν t Ŝij Note that the van-Driest weighted width Δ c is used in conjunction with SM model.

Regularization term (RT)

In this approach, no explicit SGS model is used and the dissipative effect of the subgrid motions is mimicked through a regularization procedure. Several authors propose the use of hyperviscosity [START_REF] Borue | Self-similar decay of three-dimensional homogeneous turbulence with hyperviscosity[END_REF][START_REF] Cook | A high-wavenumber viscosity for high-resolution numerical methods[END_REF] or spectral vanishing viscosity [START_REF] Karamanos | A spectral-vanishing viscosity method for large-eddy simulations[END_REF] concepts to provide sufficient dissipation at the smallest scales without recourse to an explicit model. The regularization can also originate directly from the numerical procedure, as in MILES (monotonically-integrated LES), where the dissipative properties of the numerical scheme constitutes an implicit model [START_REF] Boris | New insights into large eddy simulation[END_REF][START_REF] Fureby | Large eddy simulation of high-Reynolds number free and wall-bounded flows[END_REF]. Following this idea, a methodology combining low-dissipation schemes with the use a high-order filtering has been developed by Stolz et al. [START_REF] Stolz | Large-eddy simulation of high-Reynolds-number supersonic boundary layers using the approximate deconvolution model and rescaling and recycling technique[END_REF], Rizzetta et al. [START_REF] Rizzetta | A time-implicit high-order compact differencing and filtering scheme for large-eddy simulation[END_REF], or Bogey & Bailly [START_REF] Bogey | Large eddy simulations of round jets using explicit filtering with/without dynamic Smagorinsky model[END_REF]. In the present study, the selective filter (3), defined as part of the numerical discretization, acts as the regularization term. It also provides a smooth defiltering by removing the fluctuations at wavenumbers greater than the finite-difference scheme resolvability. As demonstrated by Mathew et al. [START_REF] Mathew | An explicit filtering method for large eddy simulation of compressible flows[END_REF], the effect of the Approximate Deconvolution Model (ADM) is globally similar to the convolution with an explicit filter. The effect of SGS motions is taken into account implicitly in the high-wavenumber range thanks to the smooth truncations of the filter and finite-difference schemes in the wavenumber space. When the selective filter is used alone, the strategy is referred to as RT (for Regularization Term). Such a procedure does not imply any additional cost for the explicit computation of SGS terms.

Discussion on LES Resolution at Re = 10 595

A detailed investigation is carried out for Re=10 595, which is the nominal configuration in most of the published works. In this part, some criteria are presented to quantify the LES resolution following the guidelines proposed by Davidson [START_REF] Davidson | Large eddy simulations: how to evaluate resolution[END_REF][START_REF] Davidson | How to estimate resolution of an LES of recirculating flow[END_REF]. First, a grid convergence study is conducted for the mean profiles. Second, the resolution in wall units, commonly used for wall-bounded turbulence [START_REF] Piomelli | Large-eddy simulations: theory and applications[END_REF], is discussed. Energy spectra in time are then shown to identify the LES cutoff and, finally, two-point correlations are used to assess the resolution in the separated region.

Grid convergence study for mean profiles

A grid sensitivity study is carried out for the baseline scheme (FD DRP11/SF DRP11) and the RT strategy by using the four grids described in Section 2.2. In Figs. 4 and 5, 10 vertical profiles at locations x/ h ∈ {0.05; 0.5; 1; 2; 3; 4; 5; 6; 7; 8} are presented for the mean velocities and turbulent intensities. Results obtained for the four grid levels are compared with the PIV data of Rapp and Manhart [START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF][START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF]. For the mean streamwise velocity, a very good match with the measurements is already observed for the coarsest grid (64 × 33 × 32). Small differences can be noticed for the grid 128 × 64 × 64 in the recovery region (5 ≤ x/ h ≤ 7). The discrepancies are more visible for the wall-normal component, which has a lower amplitude. The vertical velocity at x/ h=8 is particularly sensitive to grid refinement, even if the two finest grids are almost superimposed. The v /U b profile at the birth of the shear layer (x/ h=0.5) is also prone to errors. The profile obtained with finest grid is in very good agreement with PIV measurements, whereas it can be noted that the near-wall peak is significantly overpredicted for the grid 128 × 64 × 64.

Turbulent intensity profiles (Fig. 5), and especially the streamwise component u u , are more sensitive to grid resolution. For the finest grid of 512 × 256 × 256, profiles of v v and u v are in excellent agreement with the experimental data, but the levels of u u are slightly underpredicted, notably in the recovery of the bottom boundary layer for 4 ≤ x/ h ≤ 7. One striking aspect of Fig. 5 is that the coarse-and the two fine-grid results are close to the reference data whereas results with the 128×64×64 grid exhibit significant deviations, with a clear underestimation of the turbulent intensities. This suggests a non monotonic grid convergence, as also noted by other authors e.g. [START_REF] Xia | Constrained large-eddy simulation of separated flow in a channel with streamwise-periodic constrictions[END_REF]. Possible reasons for this non monotonic behavior are discussed in Section 5, and more specifically in Section 5.2.

Finally, for a more detailed assessment of the present solution, enlarged views of the more sensitive quantities (namely, the crosswise velocity v , and the stresses u u , v v ) are presented in Fig. 6 at two sensitive location (x/ h = 0.5 and x/ h=4). Our solutions on the finest grid are compared with the reference LES of Breuer et al. [START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF] and with the recent incompressible DNS of Krank et al. [START_REF] Krank | Direct numerical simulation of flow over periodic hills up to Re H = 10595[END_REF], based on a high-order discontinuous Galerkin code and a grid of about 2 × 10 8 degrees of freedom. A significant discrepancy with the [START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF][START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF] experimental data is also observed for the v v component at x/ h=0.5. At this position, the PIV is not able to capture the peak in the early shear layer. Our simulations are in good agreement with the DNS [START_REF] Krank | Direct numerical simulation of flow over periodic hills up to Re H = 10595[END_REF], thus supporting the conclusion that a high-resolution scheme and/or a very fine grid are required to capture the peak of Reynolds stresses in the shear layer. On the other hand, some dispersion is observed at x/ h = 4, namely for u u , both for the numerical solutions and the experimental data. Rapp and Manhart [START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF] have discussed the reliability of the measurements by comparing PIV and LDA (Laser Doppler Anemometry) techniques at this location. LDA was found to provide slightly smaller values of the streamwise fluctuations, in better compliance with our LES. However, given the rather small values taken by the observed quantities at this station, we conclude that an overall good agreement is achieved, namely among the numerical solutions.

Resolution in wall coordinates

The resolution in terms of wall coordinates Δx + i = u τ Δx i /ν is illustrated in Fig. 7. The streamwise evolution of the cell sizes along the lower wall is provided in Fig. 7a [START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF][START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF] reported in Fig. 7b for the four grids, ranges from O [START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF] wall units for the coarsest grid to O(1) for the finest one, with local maxima in the reattachment region and in the windward part of the hill (acceleration region).

The averaged values over x/ h and the maximal value are reported in Table 1. The recommendation of Piomelli and Chasnov [START_REF] Piomelli | Large-eddy simulations: theory and applications[END_REF] Δy + ≈ 1 is only satisfied in average for the finest grid. The criteria for a well-resolved LES of attached flows are Δx + ≈ 50-100 and Δz + ≈ 15-30 in the streamwise and spanwise directions, respectively. Since resolutions in x and z are similar for the present grids, the limiting values would be for Δz + . Table 1 indicates that the criterion for Δz + is only satisfied in average for the two finest grids. However, it will be shown in the following that wall unit criteria are overly pessimistic when a high-resolution scheme is used.

Velocity spectra

Figure 8 shows spectra for the three components of velocity fluctuations in the shear layer at (x, y)/ h=(2. [START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF]. The black dashed line has a slope of f -5/3 , and the vertical lines represent the estimated cut-off limit of the numerical scheme is consistent with their numerical strategy with a second-order solver and a similar number of cells 196 × 128 × 186 but clustered in the lower part of the domain. Spectra for other locations (not reported) have similar characteristics.

Davidson observed for channel flow [START_REF] Davidson | Large eddy simulations: how to evaluate resolution[END_REF] and for a recirculating flow in a plane asymmetric diffuser [START_REF] Davidson | How to estimate resolution of an LES of recirculating flow[END_REF] that energy spectra (in time or in spanwise direction) are not a good measure of LES resolution because the peak of SGS dissipation takes place at lower frequency/wavenumber than the grid cutoff (length scale corresponding to 10 cells or more). It is our opinion that the observed behavior is partly related to the fact that the grid acts as an 'implicit' filter, as commonly done in traditional LES. The cutoff due to the numerical methods is somewhat different, yielding subfilter scales [START_REF] Carati | On the modelling of the subgrid-scale and filtered-scale stress tensors in large-eddy simulations[END_REF][START_REF] Domaradzki | A subgrid-scale model based on the estimation of unresolved scales of turbulence[END_REF][START_REF] Gullbrand | The effect of numerical errors and turbulence models in large-eddy simulations of channel flow, with and without explicit filtering[END_REF]. In the present calculations, the explicit filter (used as part of the numerical discretization) removes subfilter scales, so that the cutoff in energy spectra corresponds better to the actual peak dissipation. This was illustrated for instance by Lund & Kaltenbach [START_REF] Lund | Experiments with explicit filtering for LES using a finite-difference method[END_REF]. Explicit filtering provides a better scale separation, so that energy spectra provide a more reliable estimate of the LES resolution capabilities.

Two-point correlations

By analyzing criteria for the assessment of the LES flow resolution, Davidson [START_REF] Davidson | Large eddy simulations: how to evaluate resolution[END_REF][START_REF] Davidson | How to estimate resolution of an LES of recirculating flow[END_REF] concluded that the most reliable one for a general flow (e.g. with recirculation regions) is provided by two-point correlations. The latter show by how many cells the largest scales are resolved. Since the most critical direction identified in the study of the resolution in wall units is the spanwise one, two-point correlations in this direction are presented in Fig. 9. They are computed for the three velocity components at different locations, using the time dimension to enhance their statistical convergence. Two representative plots for the streamwise velocity are shown. The point just after the separation in Fig. 9a is identified as the most critical location, as also observed by Fröhlich et al. [START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF]. A slower correlation decay is observed further in the shear layer, as illustrated by Fig. 9b. The decay of correlation coefficients from 1 to 0 is resolved by [START_REF] Meyers | Optimal model parameters for multi-objective large-eddy simulations[END_REF][START_REF] Davidson | Large eddy simulations: how to evaluate resolution[END_REF][START_REF] Xia | Constrained large-eddy simulation of separated flow in a channel with streamwise-periodic constrictions[END_REF][START_REF] Piomelli | Large-eddy simulations: theory and applications[END_REF] grid points for the first location and (6, 12, 26, 58) grid points for the second one for the grids with 64 × 33 × 32 to 512 × 256 × 256 points. Using a criterion of 6 to 10 points to describe the decay [START_REF] Davidson | Large eddy simulations: how to evaluate resolution[END_REF], it can be inferred that the two coarsest grids are not sufficient to resolve the rapid short-range correlation decay a b

Fig. just after the separation. The value of 4 points obtained for the 64 × 33 × 32 grid, which corresponds to the maximum scheme resolvability, indicates that the size of the largest scales are imposed by the numerical scheme. A similar study in the streamwise direction was conducted by Mokhtarpoor et al. [START_REF] Mokhtarpoor | Dynamic unified RANS-LES simulations of high Reynolds number separated flows[END_REF], who showed that the streamwise direction is not critical even with coarse grids.

Investigation of Sensitive Flow Features at Re = 10 595

The preceding section indicates that enhanced criteria like spectra and correlations allow to better evaluate LES resolution, notably when a separated flow region is present. Nonetheless, the level of validation is most of the time restricted to averaged quantities and a detailed analysis of unsteady flow features is required to better understand why and where the resolution criteria are undermined.

For that purpose, results obtained for various grid resolutions at Re=10 595 are extensively compared with the available data, to highlight the more sensitive flow regions and quantities of interest. We focus more specifically on the topology of the separation bubble at the hill crest and on the dynamics of the turbulent shear layer and its early development after the separation, whose resolution plays a key role for the quality of LES simulations, as shown hereafter.

Influence of compressibility effects on flow separation

The mean friction coefficient

C f = τ w /(0.5ρ ref U 2 b )
at the lower wall provides an overview of the separated flow region. The flow separates in the vicinity of the hill crest and reattaches near 4.2h, leading to a large recirculation bubble. A recovery region is then encountered followed by a secondary recirculation zone. The latter is located at the windward foot of the hill, close to x/ h=7, where the skin friction becomes locally negative. The wall shear stress then strongly increases due to the flow acceleration along the windward side of the next hill. Results corresponding to various grid resolutions are reported in Fig. 10. A clear convergence trend is observed: the maximum of the wall-shear stress on the windward face increases as the grid is refined, the two finest grid levels leading to closer solutions. The negative values in the recirculation for 2 ≤ x/ h ≤ 3 are more pronounced as the resolution is increased. Another very sensitive region corresponds to the separation above the hill crest. The inset in the figure provides a close-up view of the top of the hill. The three finest grids predict similar flow topology, with a mean separation location before the top, at x/ h ≈ -0.1, followed by small oscillations that are the imprint of small instantaneous separation bubbles. The fact that C f remains negative in the mean indicates that the flow is predominantly separated. Note that the coarser grid cannot capture such small separated regions and predicts a monotonically decreasing skin friction up to x/ h ≈ 0.21. Despite the sensitivity of the separation point to the grid resolution, the computed reattachment points, located at x/ h={4.31; 4.49; 4.32; 4.305} for the coarsest to finest grids, respectively, are in good agreement with the experimental value, highlighted with a vertical bar in the C f plots.

The friction coefficient is first compared in Fig. 11 to the LES of Breuer et al. [START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF] and Fröhlich et al. [START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF], both performed with incompressible flow solvers. The results of Breuer et al. [START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF] were obtained with the second-order finite-volume LESOCC solver using the DSM model on a grid of approximately 13.1 Mpts. The incompressible flow solvers used by Fröhlich et al. [START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF], LESOCC and STREAMLES, are also based on second-order central differencing and their computational grid (roughly 4.7 Mpts) is somewhat coarser than Breuer's. However, Fröhlich et al. used wall-functions near the upper wall, so that grid points could be clustered close to the bottom wall in their computations. Our finest grid results are in good agreement with those two references but exhibit a slightly higher peak on the windward face. Another striking difference is the behavior in the separation region shown in the inset. A small region of positive values is visible on both references, indicating a precursory separation bubble, and a main separation at x/ h ≈ 0.2. Similar oscillations are observed on our fine grid results but with lower levels, so that mean flow does not recover from the precursory separation before the hill crest. LES of Breuer et al. [START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF]; LES of Fröhlich et al. [START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF] In Fig. 12, we present further comparisons of our results with the LES of Ziefle et al. [START_REF] Ziefle | Large-eddy simulation of separated flow in a channel with streamwiseperiodic constrictions[END_REF], Xia et al. [START_REF] Xia | Constrained large-eddy simulation of separated flow in a channel with streamwise-periodic constrictions[END_REF] and Diosady & Murman [START_REF] Diosady | DNS Of flows over periodic hills using a discontinuous-Galerkin spectralelement method[END_REF], who all used a compressible flow solver at a low Mach number. Ziefle et al. [START_REF] Ziefle | Large-eddy simulation of separated flow in a channel with streamwiseperiodic constrictions[END_REF] used a compressible solver with M = 0.2, as in the present study, based on a fourth-order skew-symmetric central scheme, along with a standard fourth-order Runge-Kutta algorithm and the ADM SGS model. They chose a relatively coarse grid 128 × 69 × 72 (0.64 Mpts). As with our coarsest grid, a monotonic decrease of C f is seen in the inset, with a single separation slightly after x/ h ≈ 0.2. Nonetheless, the remaining of friction distribution on the lower wall is in fair agreement with our finest grid. A Mach number of 0.1 has been chosen in the more recent calculations of Xia et al. [START_REF] Xia | Constrained large-eddy simulation of separated flow in a channel with streamwise-periodic constrictions[END_REF] and Diosady & Murman [START_REF] Diosady | DNS Of flows over periodic hills using a discontinuous-Galerkin spectralelement method[END_REF]. The former used a finite-volume compressible solver of secondorder accuracy with different SGS modeling. We show here the results with their finest grid (176 × 69 × 72, i.e. 0.87 Mpts) and the CLES-SA model (Constrained LES using Spalart-Allmaras RANS model as constraint). The latter, used a space-time Discontinuous Galerkin (DG) formulation with a 4th-order time discretization and an 8th-order spatial discretization on a mesh of 128 × 64 × 64 cells. Note that the results in the separation region of Xia et al. [START_REF] Xia | Constrained large-eddy simulation of separated flow in a channel with streamwise-periodic constrictions[END_REF] and Diosady & Murman [START_REF] Diosady | DNS Of flows over periodic hills using a discontinuous-Galerkin spectralelement method[END_REF] are intermediate between our results for the finest grids and the incompressible LES of Fig. 11, in so far as the friction curves show a small plateau in 0 ≤ x/ h ≤ 0.15 reflecting flow separation (for Xia et al.) or reattachment (for Diosady & Murman) at the hill crest. There is hence reason to suspect that compressibility effects can affect the tiny recirculation zone at the hill top. A new simulation at M = 0.1 on the finest grid 512 × 256 × 256 has been realized to investigate the potential compressibility effects. The results are included in Figs. 11 and12 (orange dash-dotted line). Our distribution of C f for M = 0.1 is now almost superimposed with the results of Xia et al. [START_REF] Xia | Constrained large-eddy simulation of separated flow in a channel with streamwise-periodic constrictions[END_REF] and Diosady & Murman [START_REF] Diosady | DNS Of flows over periodic hills using a discontinuous-Galerkin spectralelement method[END_REF] at the same Mach number, although the C f curve remains slightly below the horizontal axis in our case. Small differences with the other numerical solutions at the same Mach can be due to the somewhat coarser grid resolution near the wall of the present simulations and to small residual averaging errors. We can conclude that the magnitude of C f is sensitive to the compressibility effects. Comparisons of the present LES of Ziefle et al. [START_REF] Ziefle | Large-eddy simulation of separated flow in a channel with streamwiseperiodic constrictions[END_REF]; LES of Xia et al. [START_REF] Xia | Constrained large-eddy simulation of separated flow in a channel with streamwise-periodic constrictions[END_REF]; LES of Diosady et al. [START_REF] Diosady | DNS Of flows over periodic hills using a discontinuous-Galerkin spectralelement method[END_REF] LES on the finest grid using M = 0.2 and M = 0.1 (not reported for brevity) show that the mean velocity and turbulent intensity profiles are almost unchanged. Generally speaking, compressibility has a relatively small effect on turbulent eddies in wall-bounded flows. Nonetheless, a known effect of compressibility is to cause earlier separation as compared with the incompressible flow [START_REF] Howarth | Concerning the effect of compressibility on laminar boundary layers and their separation[END_REF]. That is why the compressibility effects are essentially limited to the tiny recirculation bubble above the hill crest.

Shear layer dynamics and collective interaction phenomenon

In Section 4.1 we pointed out the non monotonic convergence behavior of the numerical solution, also encountered in other studies available in the literature. To shed some light on the physical mechanisms that may be responsible for the observed behavior, an analysis of the most influential flow structures is carried out. A thorough study of the vortex systems characterizing the 2D hill flow can be found in Fröhlich et al. [START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF].

By inspecting instantaneous snapshots of the fluctuating pressure, we could identify vortical structures by local minima but a great variability was noted for different time instants or grid resolutions. Nevertheless, the trace of small-size vortices is observed recurrently for x/ h < 1, i.e. at the beginning of the shear layer, which is also characterized by the highest levels of turbulent kinetic energy. Shortly downstream, large scales emerge rapidly, which is reminiscent of a collective interaction phenomenon, first described by Ho and Huang [START_REF] Ho | Subharmonics and vortex merging in mixing layers[END_REF] for mixing layers excited at low frequencies. In collective interaction, sudden formation of a large-scale vortex arises from the fusion of a number of smaller vortices shed at the natural frequency of the shear layer. This phenomenon has been observed experimentally for impinging jets [START_REF] Ho | Dynamics of an impinging jet. Part 1. The feedback phenomenon[END_REF] or in numerical simulations of rectangular cavity flows [START_REF] Gloerfelt | Large eddy simulations of high Reynolds number cavity flows[END_REF]. We may wonder if such mechanism is also present in the 2D hill flow and, in any case, if the early development of the shear layer plays a role in the non monotonic convergence.

For that purpose, the mean streamwise velocity field u is used to characterize the early growth of the mixing layer just after the hill restriction. Profiles at successive streamwise locations are extracted by interpolating the velocity field on vertical lines. Then the centerline of the mixing layer is determined by searching loci y 0.5 where u =0.5U b . The u -profiles are plotted as a function of the reduced coordinate (yy 0.5 )/ h in Fig. 13a. In

a b

Fig. [START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF] Shear-layer development. a Vertical extractions of mean streamwise velocity profiles from x/ h=0.05 to x/ h=1.05 every 0.1 (from blue to red) for grid 512×256×256. For each profile, the upward and downward black triangles respectively denote the upper and lower bounds, y up and y low , for the integration of the momentum thickness, Eq. order to examine the growth of the shear layer after separation from the hill crest, Fig. 13b shows the evolution of the momentum thickness defined as :

δ θ = y up y low u U b 1 - u U b dy (5)
where the integration bounds y up and y low are plotted in Fig. 13a. Due to the presence of the main recirculation below the shear layer, the lower bound is determined as the location below the centerline where d u /dy=0. The upper bound is defined as the location of the inflexion point above the centerline where d 2 u /dy 2 =0. The momentum thickness growth is linear by parts and three regions can be distinguished: the first is located between x/ h = 0 and 0.3 and is characterized by a growth rate dδ θ /dx 0.03 (cyan dash-dotted line), the second one (between x/ h = 0.3 and 0.8) has dδ θ /dx 0.06 (magenta dash-dotted line). The third region extends after x/ h = 0.8, where the boundaries of the shear layer become hard to define. Some discrepancies can be noted between the four grid resolutions. Coarser grids predict greater initial momentum thickness and earlier saturation. This effect is particularly evident for the 64 × 33 × 32 grid. Once again, results on the two finest grids are very close to each other. The initial growth rate is close to the value of 0.035 observed for turbulent mixing layer [START_REF] Oster | The forced mixing layer between parallel streams[END_REF]. Since the channel center is fulfilled by turbulent structures turning into the computational box, the shear layer experiences a highly turbulent environment and therefore spreads as a free turbulent mixing layer. After this initial phase, the growth rate is even higher. High values of dδ θ /dx 0.06 have been previously noted for a mixing layer [START_REF] Oster | The forced mixing layer between parallel streams[END_REF] or high-Reynolds-number cavity flows [START_REF] Gloerfelt | Numerical evidence of mode switching in the flow-induced oscillations by a cavity[END_REF]. It is believed that this is the indication of a collective interaction generating large-scale rolls in the mixing layer.

This first region should correspond to the formation of Kelvin-Helmholtz vortices at the natural frequency of the shear layer given by the linear theory of stability of parallel flows. Following Michalke [START_REF] Michalke | On the inviscid instability of the hyperbolic-tangent velocity profile[END_REF], the frequency is f 0.82U b /δ ω , where the vorticity thickness is δ ω = U b ∂ u ∂y max . The self-similarity of the velocity profiles in the early development of the shear layer is shown in Fig. 14a by dividing the reduced coordinate by the momentum thickness. Vertical velocity spectra along the centerline of the mixing layer are presented in Fig. 14b to check the presence of a tonal component at the Kelvin-Helmholtz frequency. In the region where the shedding of Kelvin-Helmholtz vortices is expected, an b a Fig. 14 a Self-similarity of the streamwise velocity profiles in the early mixing layer for x/ h=0.1 to x/ h=0.35 every 0.025 (from blue to red). b Velocity spectra E vv for grid 256 × 128 × 128 on the centerline of the mixing layer at locations x/ h ∈ {0; 0.03; 0.06; 0.10; 0.13; 0.16; 0.19; 0.23; 0.26; 0.29; 0.33} (blue to red) enrichment of the spectra can be noticed for the vertical component but no trace of tonal component is visible. The shape is rather characteristic of a fully turbulent flow, as already observed for a downstream location in Fig. 8. We conclude that it is not possible to determine a preferred size for the coherent structures observed on instantaneous snapshots and that the great variability in time induced by the channel turbulence makes hard to identify the collective interaction process. Indeed, the cavity between the hills can be categorized as closed [START_REF] Gloerfelt | Cavity Noise[END_REF], with flow reattachment on the floor. Consequently, no feedback loop is possible to fix the size of the large-scales. As shown by Fröhlich et al. [START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF], low-frequency modulations are linked to the return time of the periodic flow. So the characteristic length of the low-frequency content is L x .

The wavelengths corresponding to Michalke's frequency are depicted with white circles in Fig. 15, by using λ = U c /f and using the rough approximation U c = 0.5U b . The figure shows close-up views of u contours in the shear layer. Grid lines have been superimposed to appreciate the number of grid points in the shear layer. Even though the qualitative behavior of the shear layer is rather well captured on all grids considered in this study, a resolution of 256 × 128×128 or higher is required to partly resolve the early mixing layer downstream of the separation point, which in turn drives the behavior of the separated region. Given the extremely small number of grid points available to describe this sensitive region, using high-resolution methods appears to be essential for capturing at least in part the early shear layer dynamics. Even in this case, the initial layer thickness is overpredicted on the two coarser grids, since only larger structures can be captured. We believe that the insufficient resolution of the initial layer, along with the formation of highly chaotic large structures in the second part of the mixing layer, is responsible for the non monotonic grid convergence. The coarser grid used in the study is indeed too coarse to capture the flow behaviour around the hill crest and at the beginning of the shear layer. Additionally, the predicted incoming boundary layer is thicker, leading as a consequence to a thicker shear layer, as clearly seen from Fig. 13b. As a consequence, the characteristic size of the structures populating the shear layer is also larger. These large structures may be sufficiently persistent to lead to high velocity fluctuation levels, and to values of the turbulent intensities that are close to those obtained on the finest grids, although for a wrong physical mechanism. Indeed, although the 64 × 32 × 32 grid under-resolves the hill top region and the initial shear layer, it provides a reasonable resolution elsewhere, at least for the present high-resolvability scheme, leading to unexpectedly good predictions. The 128 × 64 × 64 grid is also too coarse to capture in detail the dynamics of the crest region. Nevertheless, it captures at least qualitatively the same physics as the two finest grids. This is shown both by the C f plot in Fig. 10, which exhibits a similar qualitative behaviour as the 256 × 128 × 128 and 512 × 256 × 256 grids, and by the evolution of the shear layer thickness of Fig. 13b. This figure shows in particular that the shear layer is much thinner in this case, compared to the previous one. As a consequence, the structures generated are smaller. Once again, these smaller structures are not resolved well enough, as shown by the correlation plot at x/ h ≈ 1, y/ h ≈ 1, where just only 6 points are observed in the initial correlation decay. These rather ill-resolved smaller structures grow at a somewhat slower rate than the one observed on the finer grids, probably because of the explicit filter damping close to grid cutoff, finally resulting in the lower turbulent intensities observed for this grid. When grid resolution is increased furthermore (on the two finest grids), not only the initial shear layer dynamics is better represented, but additionally the grids are fine enough to ensure a good resolution of the generated structures and prevent premature damping. In addition to the preceding discussion, we stress that the ill-resolved small structures in the mixing layer are coherent, and such energetic structures can not be mimicked by SGS models. As a consequence, introducing an SGS model does not alleviate the above-mentioned inadequacies on the coarsest grids. Using high-resolution schemes allows capturing at least in part this subtle behavior using 4 Million points without any particular adaptation of the grid to the shear layer, which represents a clear advantage over standard second-order flow solvers.

Sensitivity to the Numerical Schemes and SGS Models

In this section, we focus on the effect of the spatial scheme resolvability and of the chosen SGS treatment on LES quality. All of the calculations are carried out at Re = 10 595 on coarsest grid (64 × 33 × 32) to better highlight the differences among the various schemes.

As shown in the Appendix (see also [START_REF] Gloerfelt | Investigation of high-order methods in large-eddy simulation of separated flow in a channel with periodic constrictions[END_REF] for more details), the choice of the finite difference scheme for the convective terms has a relatively small influence on the numerical solution, provided that a high-order scheme (order 4 or higher) is selected. Theoretical studies based on Fourier analysis show that the number of grid points required to resolve a given solution wavelength with a prescribed error level (purely dispersive error for the present central schemes) decreases significantly when increasing the scheme order (for standard FD). Even higher gains are obtained for the optimized scheme DRP11. Furthermore, the computational cost per iteration and per point increases slowly when enlarging the discretization stencil, so that large stencil/high-resolution schemes provide a better cost/accuracy tradeoff than low order ones. This is why only the DRP11 scheme was retained.

On the other hand, the dissipative part of the numerical discretization, here ensured by a selective filter, plays a crucial role in the representation of SGS scales and interacts with SGS modeling, as discussed in the next sections.

Influence of numerical dissipation

In the following, the influence of the filtering method is investigated by selecting DRP11 as the spatial discretization scheme (Fig. 16). The results are found to be extremely sensitive to the filter order. With the 2nd-order standard centered filter, all turbulent structures are dissipated, leading to pseudo-laminar velocity profiles and negligibly small turbulent intensities. The 4th-order and 6th-order filters are also highly overdissipative. Contrary to the results obtained for finite-difference schemes, a major influence of the filters is observed up to the 10th-order. Furthermore, the filter amplitude, χ in Eq. 3, is a tuning parameter and, similar to artificial dissipation coefficients used sometimes in conjunction with central schemes [START_REF] Jameson | Numerical solutions of the Euler equations by finite volume methods using runge-Kutta time-stepping schemes[END_REF][START_REF] Kim | Optimized compact finite difference schemes with maximum resolution[END_REF][START_REF] Tam | Dispersion-relation-preserving finite difference schemes for computational acoustics[END_REF], may affect the results. It is expected that, for filters with a sufficiently high cutoff wavenumber (and specifically with a cutoff close to the grid one), and for sufficiently fine grids, the influence of this parameter remains rather small [START_REF] Aubard | Comparison of subgrid-scale viscosity models and selective filtering strategy for large-eddy simulations[END_REF][START_REF] Berland | Filter shape dependence and effective scale separation in large-eddy simulations based on relaxation filtering[END_REF][START_REF] Fauconnier | On the performance of relaxation filtering for large-eddy simulation[END_REF]. Nevertheless, a good practice is to choose the filter amplitude as small as possible to prevent any overdamping of the smallest resolved structures. The typical value of 0.2 was selected because it was the minimal one for preventing numerical divergence in all of the considered cases. The effect of the filter amplitude is assessed in Fig. 17 for the baseline filter DRP11 and two meshes (64 × 33 × 32 and 256 × 128 × 128). The results show that increasing the filter amplitude from 0.2 to 0.6 does not affect the mean flow quantities significantly but tends to damp the turbulent intensity profiles. This effect is more visible on the v v -profiles in Fig. 17a. The damping effect is suppressed by applying the filter at a lower frequency (i.e. every k time steps with k ∈ N, instead of after each time step). For instance, a calculation with χ = 0.6 and k = 3 gives equivalent results to the baseline calculation χ = 0.2 and k = 1, as illustrated in Fig. 17b. Finally, the sensitivity of the solution to χ becomes lower on more refined grids as shown in Fig. 17c [START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF][START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF] resolved part. Specifically, we observe that: i) using a finer grid allows to achieve numerically stable results even using a lower filter amplitude (results for χ = 0.1 are reported); 2) results for χ = 0.1 and χ = 0.2 are close to each other, at least to within uncertainties associated with, e.g., statistical averaging for the chosen sensitive quantity of interest; 3) results obtained by using higher filter amplitudes (the worst case χ = 0.6 is reported in the figure) become closer to the reference than on the coarser grid. Results for other stations and/or quantities of interest (not reported for brevity) also follow the preceding trends.

Note that the coupled influence of finite-difference and filtering schemes with χ =0.2 (not reported for brevity) leads to results that are similar to those obtained by changing the filter alone, thus indicating that the accuracy of the LES solution depends mostly on the sharpness of the numerical filter.

Influence of subgrid-scale models

Afterwards, the influence of SGS models is investigated by using the DRP11 scheme for both space approximation and numerical stabilization for Re = 10 595 on the 64 × 33 × 32 grid. We stress that, since we use a compressible solver, a numerical filtering step is always necessary to ensure stability, even when an explicit SGS model is applied. The filter is expected however to play a lesser role in this case, due to the regularizing effect of the explicit SGS term. Of course, due to the high nonlinearity of the problem, a clear separation of the effects of the SGS model and numerical filter can never be achieved.

Results for various SGS models are presented in Fig. 18 for selected flow property profiles at x/ h=3. At first sight, the mean streamwise velocity is almost not affected, indicating that SGS modeling choices are less influential than the numerical schemes. Differences are nonetheless observed for the vertical velocity and turbulent intensity profiles at this position. Similar results are obtained at other locations.

The Smagorinsky model with C S =0.18, either in standard or dynamic formulation, is overly dissipative. Choosing a lower constant (C S =0.1), as recommended by several authors (e.g. [START_REF] Meyers | Optimal model parameters for multi-objective large-eddy simulations[END_REF]), or dividing by a factor of 2 the value of the dynamical constant, yields results in better agreement with the reference. Similarly, using multiscale models significantly improves solution quality, although the results are not more accurate than those obtained with the RT strategy alone. The activity of eddy-viscosity SGS models is well illustrated experiments Rapp & Manhart [START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF]; DG LES of de la Llave Plata et al. [START_REF] De La Llave Plata | On the use of a high-order discontinuous Galerkin method for DNS and LES of wall-bounded turbulence[END_REF]. Mean streamwise u a and vertical v b velocities the dominant role played by numerical filtering. The activity of the SGS models is moderate already for the coarsest grid, and a lower activity is expected by refining the grid. On the other hand, the CPU over-cost of using explicit SGS models, relative to the RT model, is equal to 6.5, 17.4, 17 and 33.4 % for the SM, MSM-ls, MSM-ss, and DSM, respectively. As a consequence, the RT approach appears to offer an attractive tradeoff between accuracy and computational cost.

It could be argued that the ratio ν t /ν is not a reliable measure of the LES resolution since it compares LES with DNS and it does not provide an independent estimate of the LES resolution [START_REF] Davidson | Large eddy simulations: how to evaluate resolution[END_REF][START_REF] Davidson | How to estimate resolution of an LES of recirculating flow[END_REF]. Furthermore, the present results are obtained for a rather low Reynolds number flow. We will show in Section 7 that the above considerations still hold at Re=37 000 on the finest grid.

Investigations at Higher Reynolds Numbers

In this section we apply our high-resolution numerical strategy to investigate higher Reynolds number cases, namely Re=19000 and 37000. On the one hand, this allows to assess the capability of the present methodology to reproduce subtle Reynolds number effects against experimental results by Rapp and Mahnart [START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF]. On the other hand, the present simulations provide a numerical database for these flow cases, to be used in the future for the benchmarking of numerical methods or models. The calculations presented in the following are based on the DRP11 finite difference scheme and filter, the latter being also used as regularization for LES (RT strategy), unless otherwise specified. [START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF]; DG LES of de la Llave Plata et al. [START_REF] De La Llave Plata | On the use of a high-order discontinuous Galerkin method for DNS and LES of wall-bounded turbulence[END_REF]. Turbulent intensities u u a, v v b, u v c

Results at Re=19000

LES are carried out using the two finest grids of Section 2.2, which ensure an average height of the first cell close to the wall Δy + 3 on the coarser grid and about 1.5 on the finer one, which approximately corresponds to the resolution achieved on the grid 256 × 128 × 128 at Re=10595.

Figures 20 and21 show profiles at different streamwise stations for the mean velocity components and the turbulent intensities, respectively. The profiles are compared with the PIV measurements at Re=19 000 of Rapp and Manhart [START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF] and with the recent compressible Discontinuous Galerkin (DG) simulation of de la Llave Plata et al. [START_REF] De La Llave Plata | On the use of a high-order discontinuous Galerkin method for DNS and LES of wall-bounded turbulence[END_REF]. The latter performed LES with WALE model and the DG method with a polynomial approximation of degree 3 (DG-p3 of order 4). They used 4. [START_REF] Ziefle | Large-eddy simulation of separated flow in a channel with streamwiseperiodic constrictions[END_REF] [START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF]; MGLET LES of Manhart et al. [START_REF] Manhart | Assessment of eddy resolving techniques for the flow over periodically arranged hills up to Re=37,000[END_REF]. Mean streamwise u a and vertical v b velocities with the optimized DRP11 scheme and the DRP11 selective filter. A very good compliance is obtained for the mean velocity components. The u -profiles obtained on the finest grid (512 × 256 × 256) are closer to the DG ones [START_REF] De La Llave Plata | On the use of a high-order discontinuous Galerkin method for DNS and LES of wall-bounded turbulence[END_REF] than to those of the RT-LES with 256 × 128 × 128 points. Small discrepancies are noted for the vertical velocity profiles at x/ h = 2 and 6 with PIV measurements. At these locations, the 256 × 128 × 128 RT-LES is closer to the experimental points, and the 512 × 256 × 256 RT-LES is closer to the DG-LES [START_REF] De La Llave Plata | On the use of a high-order discontinuous Galerkin method for DNS and LES of wall-bounded turbulence[END_REF]. The LES with the 512 × 256 × 256 grid is however in better agreement with PIV data than other LES at x/ h = 8, which was identified as a location sensitive to grid refinement.

The Reynolds stresses, reported in Fig. 21, are in reasonably good agreement with the experimental reference. The different LES give similar results in the shear layer. At almost all locations the finest grid predicts levels closer to the PIV data than the 256 × 128 × 128 grid, except at x/ h=2, where u u /U 2 b and v v /U 2 b are slightly overestimated with the finest grid. At x/ h = 2 (inside the recirculation bubble), the experimental points for u u lie in between the two RT-LES and the DG results are closer to PIV. Note also that the peak near the bottom wall for the wall-normal fluctuations at x/ h=0.5 is not present in the measurements, probably due to a lack of resolution of the PIV in this region. Unfortunately, DG results are not shown for v v , which is the most difficult to compare with PIV results.

Results at Re=37000

Two LES calculations are carried out on the finest grid level of 512 × 256 × 256 points. The first one relies on the RT strategy, whereas an explicit SGS model, namely the multiscale model in small-small version (MSM-ss), has been added in the second one. experiments Rapp & Manhart [START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF]; MGLET LES of Manhart et al [START_REF] Manhart | Assessment of eddy resolving techniques for the flow over periodically arranged hills up to Re=37,000[END_REF]. Turbulent intensities u u a, v v b, u v c Figures 22 and 23 compare the results for the mean profiles and Reynolds stresses, respectively. The PIV data of Rapp and Manhart [START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF] and the MGLET LES of Manhart et al. [START_REF] Manhart | Assessment of eddy resolving techniques for the flow over periodically arranged hills up to Re=37,000[END_REF], both available in QNET database, are plotted to judge the quality of the results. Manhart et al. [START_REF] Manhart | Assessment of eddy resolving techniques for the flow over periodically arranged hills up to Re=37,000[END_REF] used the second-order-accurate Cartesian IBM (immersed boundary method) incompressible solver MGLET with a grid of 216 × 168 × 104 points (3.8 Mpts) with a maximum Δy + of 30 and the WALE SGS model. The present LES are finer approximately by a factor two in each direction. The maximum Δy + on the windward slope of the hill is 7.6 and the average value over the bottom wall is 2.6.

The main features of the mean streamwise velocity are captured by all simulations, as shown in Fig. 22a, the RT-LES with the 512 × 256 × 256 grid being the closest to experiments. A greater scatter is observed for mean vertical velocity profiles in Fig. 22b, bearing in mind that the magnitude of v is low. The results for RT-LES are closer to measurements and the profiles for MSM-LES comply with those of MGLET LES. The comparisons of For the v v (Fig. 23b), all the simulations dramatically overpredict the levels found in the experiments. The RT-LES and MSM-LES yield similar results, whereas MGLET profiles correspond to sometimes higher and sometimes lower levels than our LES. No clear conclusion can be made. The differences with experiments may be partially explained by uncertainties in the measurements but also by choices made in the numerical setup. Finally, the prediction is rather good for the u v intensity in Fig. 23c. Some slight variations are here again difficult to interpret. The use of an explicit SGS model, MSM-ss, does not improve the results and leads to results generally closer to the coarser LES with MGLET.

Influence of the Reynolds number

The friction coefficient along the bottom wall is compared for Reynolds numbers 10595, 19000 and 37000 in Fig. 24. The RT-LES cases with the 512 × 256 × 256 grid have been retained for the following analysis of Reynolds-number effects. The distribution of wall shear stress has a similar shape for the three Reynolds numbers. The magnitude of the peak on the windward-slope of the hill is reduced as the Reynolds number is increased, which is in agreement with the observations of Breuer et al. [START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF]. The close-up view around the hill crest in the inset indicates that the precursory bubble is more intense as Re increases.

The flow separates at x/ h=0. [START_REF]DESider -A European Effort on Hybrid RANS-LES Modelling; Results of the european-Union Funded Project[END_REF] To get insight into the Reynolds number effects on the mean quantities, some comparisons are presented in Fig. 25 at the location x/ h = 4, where LDA measurements are also provided in the study of Rapp & Manhart [START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF]. Overall, the differences for the streamwise component are rather small and agree well with experimental results. At x/ h=4, the flow is still separated for Re=10595 but has already reattached for Re=19000 and 37000. The component in the vertical direction is a more sensitive quantity. The trends follow fairly the measurements of Rapp & Manhart [START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF] at the different locations, even if a slight underestimation of v /U b is found at x/ h=4 for all cases but the hierarchy is still correct. For the Reynolds stresses, a very good agreement both in terms of level and hierarchy is found for u u and u v at x/ h=4. LDA measurements were also carried out at this sensitive position, and are superimposed for Re=19000 (filled triangles). The PIV, LDA and LES are in good agreement to each other for this Reynolds number. The most striking feature for the vertical velocity fluctuations is the lower magnitude observed for Re=37000 in PIV measurements, whereas similar levels are obtained at 10595 and 19000. If the trend is reasonably well reproduced by the simulations, the decrease of v v for the highest value of the Reynolds number is less pronounced.

A hint about the resolution capability at the higher Reynolds numbers is provided by the spanwise two-point correlations, calculated just after separation (at x/ h ≈ 1 and y/ h ≈ 1). The number of points in the correlation decay from 1 to 0 is estimated at 17 and 35 points respectively for grids 256 × 128 × 128 and 512 × 256 × 256 at Re = 19000. A value of 41 points within the decay of correlation is noted for Re = 37000, which is similar to the one obtained for Re = 19000 with the same grid. We can conclude that largest scales are correctly resolved on the finest grid for the higher Reynolds numbers.

High-resolution large eddy simulations have been conducted for the separated flow in a channel with streamwise-periodic constrictions at various Reynolds number. This configuration is well documented both experimentally and numerically and represents a good benchmark for investigating LES quality. The use of high-order algorithms minimizing dissipation and dispersion errors is seen to be essential for obtaining high-quality solutions while using a reduced number of grid points. Specifically, an extensive parametric study of the effect of numerical approximations has been conducted for coarse-grid LES with 65000 grid points, which have been found to show a rather fair agreement with the references for this benchmark case, provided that a suitable approximation and filtering method is chosen. In this respect, it is shown that the dissipative part of the numerical discretization, here represented by an explicit selective filter, plays a crucial role in the quality of the results. In particular, the filter selectivity in the wavenumber space is of the outmost importance to increase the resolved part of the dynamics. It allows a clear separation between resolved and subfilter scales by removing ill-resolved structures. A properly shaped filter also represents a good surrogate for SGS modeling. Indeed, an investigation of the influence of various explicit SGS models shows that the solution is weakly affected by the SGS model in use, provided that the modeled dissipation is low enough. This result is still valid for the highest Reynolds number case on a fine grid, showing that implicit modeling strategy is cost-effective and avoids the computational overhead introduced by the explicit SGS model. Some criteria have been evaluated to judge independently the quality of LES resolution. In the present calculations, the explicit filter (used as part of the numerical discretization) removes subfilter scales, so that the cutoff in energy spectra is a good criterion to estimate the LES cut-off. Two-point correlations constitute also a reliable resolution indicator, notably in the presence of recirculation regions.

Afterwards, details of the flow dynamics are investigated to get some insight into the flow features responsible for sensitive regions where criteria are difficult to fulfill. A first important outcome of the present analysis is the role of compressibility effects on the tiny recirculation bubble that appears at the top of the hill. For Re = 10 595, the precursory bubble is merged with the main recirculation at M = 0.2, whereas an incipient reattachment is observed at the hill crest using M = 0.1, and a clear precursory recirculation is obtained with incompressible solvers. The mean profiles remains almost unaffected by compressibility. The intensity of the precursory bubble is also seen to increase as the Reynolds number is higher. Furthermore, the early birth of the separated shear layer is found to be highly sensitive to local grid resolution (namely in the streamwise direction). The mechanism leading to large vortices in the shear layer indeed bears similarities with the collective interaction phenomenon, which has been encountered for mixing layer excited by a low-frequency component. In the 2D hill flow, this complex process is however highly disrupted by the turbulent fluctuations turning into the periodic box. A high spreading rate of the mixing layer has been observed but no preferred frequency is noticeable in the spectral analysis. This complex dynamics at the birth of the shear layer is hardly resolved with coarse grids and can be a hint to explain non monotonic grid convergence observed at Re = 10 595 for the coarsest grids.

Finally, new databases have been established for the higher Reynolds number cases, Re = 19 000 and 37 000, for which reliable PIV measurements are also available for validation. Our LES computations are in nice accordance with the experimental data. Furthermore, Reynolds-number effects are fairly well taken into account by the new calculations. The new high-resolution results can constitute useful databases for these higher Reynolds numbers, for which very few numerical solutions have been made available up to now.

The fine grid LES data for the different cases are made available for download on the following website: https://www.researchgate.net/publication/315413324 Benchmark database 2D periodic hill flow Acknowledgments This work was granted access to the HPC resources of IDRIS and CCRT under the allocation 1736 made by GENCI (Grand Equipement National de Calcul Intensif).

Compliance with Ethical Standards

Conflict of interests The authors declare that they have no conflict of interest.

Appendix: Influence of the Spatial Discretization Resolvability on CPU Costs

In this Appendix we report selected numerical results illustrating the role of the spatial discretization scheme for the convective terms for achieving a given LES quality/computational cost compromise. Specifically, results are presented for the 2D-hill flow at Re = 10 595 using various FD schemes defined by Eq. 2. As in Section 6, we focus on the coarsest grid 64 × 32 × 32.

Two families of central finite-difference schemes are tested, namely the standard versions (second-to tenth-order scheme, denoted FDo2 to FDo10), and the optimized DRP11 scheme. Numerical stabilization and RT modeling are achieved by applying the baseline filter, i.e. the DRP11 filter, in all cases. Profiles of u , v and v v at x/ h = 3 are reported in Fig. 26. Standard finite-difference schemes of various orders provide closeby results, except for the 2nd-order scheme, which introduces significant phase errors on the considered grid. This suggests that, beyond 2nd-order of accuracy, the quality of the solution becomes weakly dependent of the consistent part of the spatial scheme. However, for the 11-point schemes, and more specifically the optimized DRP11, the results tend to be more accurate.

The preceding results are in line with theoretical studies based on Fourier analysis, showing that the accuracy limit, i.e. the number of grid points required to represent a given wavelength with a prescribed error level, decreases significantly when increasing the order of accuracy (for standard schemes) or by optimizing the coefficients (for DRP schemes). 2 From top to bottom: accuracy limits (points per wavelength λ/Δx ) of the standard and optimized finite differences (FD) to achieve a dispersion error of 5 × 10 -4 ; accuracy limits multiplied by the number of points p = 2N + 1 of the stencil; CPU costs (Re = 10 595; grid 64 × 33 × 32, RT model) relative to the baseline scheme FD DRP11 /SF DRP11 ; cost criterion corresponding to the CPU cost, multiplied by the minimal number of points per wavelength to achieve the prescribed dispersion error The accuracy limits required by the present family of schemes is reported in Table 2. The error criterion corresponds to the upper value considered in [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for noise computation[END_REF]. Of course, the price to pay is an increased computational complexity due to the larger and larger stencil in use, leading to a higher computational complexity. This criterion, also reported in Table 2, tends to saturate and even slightly increases for orders higher than 6. Note that the normalized resolvability (i.e. the product of the number of points per wavelength by the number of points in the stencil) remains always better for the DRP11 scheme. However, the computational complexity of the FD scheme is an over-simplistic cost estimate, not well suited to characterize the overall CPU cost associated with a given scheme. Indeed, the discretization of the convective terms is only one of the many ingredients involved in the flow solver. This is why, in Table 2, we report the effective computational costs of 2D-hill simulations based on various FD schemes and the DRP11 filter (normalized with respect to the CPU cost of the baseline DRP11 scheme), along with a cost criterion corresponding to the CPU cost per iteration and per point multiplied by the accuracy limit. The table shows that: i) the brute CPU cost is weakly dependent on the chosen FD scheme; ii) the cost criterion is much lower for the high-resolution schemes, due to the considerably reduction in terms of grid points. The preceding results were obtained by applied the DRP11 selective filter. In order to discuss the computational cost associate for the filter, in Table 3 we report the accuracy limits (points per wavelength required to achieve a dissipation error threshold of 5×10 -4 /χ ) of the various selective filters, along with the effective CPU times per iteration and per grid point relative to the baseline scheme. In this case, all the simulations are based on the DRP11 FD scheme for the convective term. A cost criterion corresponding to CPU multiplied by the filter accuracy limit is also reported. The results were obtained for the typical filter amplitude χ = 0.2. The table shows that the brute CPU cost is even less sensitive to the chosen filter (applied only one time per time step) than to the FD scheme, Table 3 From top to bottom: accuracy limits (points per wavelength λ/Δx ) of the standard and optimized selective filters (SF) (χ = 0.2) to achieve a dissipation error of 5 × 10 -4 /χ ; CPU costs (Re = 10 595; grid 64×33×32, RT model) relative to the baseline scheme FD DRP11 /SF DRP11 ; cost criterion corresponding to the CPU cost, multiplied by the minimal number of points per wavelength to achieve the prescribed dissipation error while the dissipation accuracy limits vary even more sharply when increasing the filter resolution than the dispersion error limits. As a consequence, the normalized CPU times decreases greatly when increasing the filter resolution, achieving a minimum in the case of the DRP11 filter.

FD

Selective

Fig. 1

 1 Fig. 1 Periodic hill test case. The computational domain of length L x = 9h is replicated 3.5 times. The mean streamwise velocity and streamlines are shown in the first part and then the instantaneous streamwise velocity for Re=10 595 is depicted. In the last part (22.5 ≤ x/ h ≤ 31.5), high-resolved PIV (Particle Image Velocimetry) results from Kähler et al. [23] at Re=8 000 are shown

  The reference density is p ref /(rT ref ) with p ref =1 atm and T ref =298.15 K the wall temperature. A relaxation term based on the effective integrated mass-flow rate Q neff at the hill crest is added with a small coefficient α = 10 -3 to avoid too fast variations of the forcing term.

Fig. 2

 2 Fig. 2 View of the topology of the grid with 64 × 33 × 32 points

Fig. 3

 3 Fig. 3 Fourier analyses of central finite-difference and filter schemes. a Effective wavenumber k * x = 2 N l=1 a l sin (lk x) of finite-difference schemes as function of the reduced wavenumber k x. The magenta dashed line represents the exact relationship. b damping function D(k x) = d 0 + N l=1 2d l cos(lk x) of filters as function of k x. The different schemes are standard second-order (o2); fourth-order (o4) ; sixth-order (o6); eighth-order (o8) ; tenth-order (o10) and optimized 11-point (DRP11)

Fig. 4

 4 Fig. 4 Mean velocity profiles at Re=10595, u /U b a and v /U b b: 64×33×32; 128 × 64 × 64; 256 × 128 × 128; 512 × 256 × 256; PIV measurements of Rapp & Manhart[START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF][START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF]

Fig. 5

 5 Fig. 5 Turbulent intensity profiles atRe=10595, u u /U 2 b a, v v /U 2 b b and u v /U2 b c: 64 × 33 × 32; 128 × 64 × 64; 256 × 128 × 128; 512 × 256 × 256; PIV measurements of Rapp & Manhart[START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF][START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF]

Fig. 6 Fig. 7

 67 Fig. 6 Profiles of v a,d; u u b,e and v v c,f at x/ h = 0.5 a-c and x/ h = 4 d-f for Re=10595: present 512 × 256 × 256 ; LES of Breuer et al. [10] ; DNS of Krank et al. [55] ; PIV and LDA ( ) measurements of Rapp & Manhart [24]

Fig. 10

 10 Fig. 10 Friction coefficient along the bottom wall at Re=10 595: 64 × 33 × 32; 128 × 64 × 64; 256 × 128 × 128; 512 × 256 × 256

Fig. 11

 11 Fig. 11 Friction coefficient along the bottom wall at Re=10 595: 512 × 256 × 256 (M = 0.2); 512 × 256 × 256 (M = 0.1);LES of Breuer et al.[START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF]; LES of Fröhlich et al.[START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF]

Fig. 12

 12 Fig. 12 Friction coefficient along the bottom wall at Re=10 595: 512 × 256 × 256 (M = 0.2); 512 × 256 × 256 (M = 0.1);LES of Ziefle et al.[START_REF] Ziefle | Large-eddy simulation of separated flow in a channel with streamwiseperiodic constrictions[END_REF]; LES of Xia et al.[START_REF] Xia | Constrained large-eddy simulation of separated flow in a channel with streamwise-periodic constrictions[END_REF]; LES of Diosady et al.[START_REF] Diosady | DNS Of flows over periodic hills using a discontinuous-Galerkin spectralelement method[END_REF]

Fig. 15

 15 Fig. 15 Close-up views of the mixing layer growth for the different grid levels: a 64×33×32; b 128×64×64; c 256 × 128 × 128; d 512 × 256 × 256. The colormap represents the mean streamwise velocity with the grid depicted; the white bullets represent the centroids of the mixing layer. and the white circle the theoretical wavelengths of Kelvin-Helmholtz instabilities

Fig. 16

 16 Fig.[START_REF] Temmerman | Investigation of wall-function approximations and subgrid-scale models in large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF] Influence of selective-filtering schemes supplemented by the DRP11 finite differences at Re = 10 595: Standard SFo2 , SFo4 , SFo6 , SFo8 , SFo10 and SF DRP11. PIV measurements of Rapp & Manhart[START_REF] Breuer | Flow over periodic hills -numerical and experimental study in a wide range of Reynolds numbers[END_REF][START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF]. Profiles of u a; v b and v v c at x/ h = 3

Fig. 18 Fig. 19 Fig. 20

 181920 Fig. 18 Influence of subgrid-scale models at Re = 10 595: SM C S =0.18 , SM C S =0.1 , DSM C d , DSM C d /2 , MSM-ls , MSM-ss and RT . PIV measurements of Rapp & Manhart [10, 24]. Profiles of u a; v b and v v c at x/ h = 3

Fig. 21

 21 Fig. 21 Periodic hill flow at Re = 19000 : RT 256 × 128 × 128; RT 512 × 256 × 256; experiments Rapp & Manhart [24]; DG LES of de la Llave Plata et al. [75]. Turbulent intensities u u a, v v b, u v c

Fig. 23

 23 Fig. 23 Periodic hill flow at Re = 37000: RT 512 × 256 × 256; MSM-ss 512 × 256 × 256;experiments Rapp & Manhart[START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF]; MGLET LES of Manhart et al[START_REF] Manhart | Assessment of eddy resolving techniques for the flow over periodically arranged hills up to Re=37,000[END_REF]. Turbulent intensities u u a, v v b, u v c

Fig. 24

 24 Fig.[START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF] Influence of Reynolds number on friction coefficient (512×256×256 grid): Re=10595; Re=19000; Re=37000

Fig. 25

 25 Fig. 25 Influence of Reynolds number for grid 512 × 256 × 256 at x/ h = 4: Re = 10595 ; Re = 19000 ; Re = 37000 ; experiments of Rapp & Manhart [24] at Re = 10595 ; Re = 19000 ; Re = 37000 . Mean velocity profiles u a, v b and turbulent intensities u u c, v v d, u v e. Filled symbols for u u /U 2b profiles are LDA measurements[START_REF] Rapp | Flow over periodic hills: an experimental study[END_REF]

Fig. 26

 26 Fig. 26 Influence of finite-difference schemes supplemented by the DRP11 filter at Re = 10 595: Standard FDo2 , FDo4 , FDo6 , FDo8 , FDo10 and FD DRP11 . Profiles of u a; v b and v v c at x/ h = 3

Table 1

 1 Grid resolution in wall units for LES at Re = 10595

	a Δx + i for grid 512 × 256 × 256;

  [START_REF] Vreman | Comparison of numerical schemes in large-eddy simulation of the temporal mixing layer[END_REF][START_REF] Ghosal | An analysis of numerical errors in large-Eddy Simulations of turbulence[END_REF]. The spectra are computed through the Welch method applied to signal segments with the same length (approximately 70 time units), leading to a sampling frequency of 0.014U b /h. The duration of the signals in non-dimensional time units and the number of averaging segments are given in the figure caption for each grid. Each segment is windowed by a Hann function with no overlap and the results are presented as onesided spectra with amplitude corrected by the window factor (8/3). The approximate cutoff frequency of the numerical scheme, u local /(4Δx), is also reported. We can see that the effective frequency fall-off is imposed by the cell size, and that all simulations are indeed LES. For this Reynolds number, the inertial subrange with a slope of -5/3 is well-defined over more than one decade of frequencies. The present spectra are in good agreement with those of Fröhlich et al.

[START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF]

at (x, y)/ h=

(2.23,1.13)

in the low and medium frequencies, and the grid cut-off of

[START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF]

roughly corresponds to that of our 256 × 128 × 128 grid. This a b d Fig. 8 Velocity spectra E uu a; E vv b and E ww c at (x, y)/ h=(2.2,1) for Re=10595. 64 × 33 × 32 (T =244 divided in 4 segments); 128 × 64 × 64 (T =762 divided in 9 segments); 256 × 128 × 128 (T =202 divided in 3 segments); 512 × 256 × 256 (T =133 divided in 2 segments); Fröhlich et al.

  Markers + are added on the lines for 64 × 33 × 32 and 128 × 64 × 64 grids to indicate the resolution

9 Spanwise two-point correlation coefficients at two locations (shown in the inset) for the different grid resolutions: 64 × 33 × 32; 128 × 64 × 64; 256 × 128 × 128; 512 × 256 × 256.

  [START_REF] Tessicini | Approximate near-wall treatments based on zonal and hybrid RANS-LES methods for LES at high Reynolds numbers[END_REF] Influence of filtering coefficient χ for periodic hill at Re = 10 595. Profiles of v v at x/ h = 3

	a grid 64 × 33 × 32	χ=0.2;	χ=0.4;	χ =0.6; b grid 64 × 33 × 32
	χ =0.2;	χ=0.4 every two iterations;	χ =0.6 every three iterations; c
	grid 256 × 128 × 128	χ=0.1;	χ=0.2;	χ =0.6.	PIV measurements of
	Rapp & Manhart

  millions of degrees of freedom, corresponding to 64 × 16 × 64 elements. This is equivalent to the resolution of our grid 256 × 128 × 128

	a
	b
	Fig. 22 Periodic hill flow at Re = 37000:	RT 512 × 256 × 256;	MSM-ss 512 ×
	256 × 256;	experiments Rapp & Manhart

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.