M Catchirayer

J.-F Boussuge

P Sagaut

) M Montagnac

X Garnaud

Extended integral wall-model for large-eddy simulations of compressible wall-bounded turbulent flows

all-modeling is required to make large-eddy simulations of high-Reynolds number wall-bounded turbulent flows feasible in terms of computational cost. Here, an extension of the integral wall-model for large-eddy simulations (iWMLESs) for incompressible flows developed by Yang et al. ["Integral wall model for large eddy simulations of wall-bounded turbulent flows," Phys. Fluids 27(2), 025112 (2015)] to compressible and isothermal flows is proposed and assessed. The iWMLES approach is analogous to the von KármÍ

. INTRODUCTION

Computational fluid dynamics (CFD) is widely used as a research or engineering tool in fluid mechanics. Reynolds-Averaged Navier-Stokes (RANS) modeling, which is based on statistical turbulence models, is the most popular paradigm. RANS simulations allow mean quantities to be computed at a reasonable accuracy and computational cost even for complex engineering flows. [START_REF] Gourdain | RANS and LES for the heat transfer prediction in turbine guide vane[END_REF] But RANS modeling shows its limits in predicting highly unsteady phenomena like separation. Large-Eddy Simulation (LES) is a way to overcome these limitations and the next step toward achieving more accurate simulations. [START_REF] Sagaut | Large-eddy Simulation for Incompressible Flows[END_REF][START_REF] Garnier | Large-Eddy Simulation for Compressible Flows[END_REF][START_REF] Sagaut | Large eddy simulation for aerodynamics: Status and perspectives[END_REF] In RANS simulations, the whole turbulent spectrum is modeled, while LES captures the large energetic scales of the flow, and only the small scales are modeled. However, nearwall turbulence contains very small scales (near-wall streaks [START_REF] Smits | High-Reynolds number wall turbulence[END_REF] ) and simulating them would require prohibitively fine meshes. A first estimation of LES computational cost, with or without wall-modeling, was given by Chapman. [START_REF] Chapman | Computational aerodynamics development and outlook[END_REF] Using a more accurate correlation for high Reynolds number flows, Choi a) b) c) Electronic mail: pierre.sagaut@univ-amu.fr and Moin [START_REF] Choi | Grid-point requirements for large eddy simulation: Chapman's estimates revisited[END_REF] recently updated these estimations: the number of mesh points to solve a boundary layer flow is approximately proportional to Re 13/7 in the inner layer and Re 1 in the outer layer with Re being the Reynolds number. In typical aeronautical applications, one has Re ∼ 10 7 ; thus, the inner layer requires more than 99% of the boundary layer grid points as illustrated by Piomelli and Balaras. [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF] Hybrid RANS/LES methods are a way to reduce this computational cost. Several reviews of these methods exist. [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF][START_REF] Sagaut | Multiscale and Multiresolution Approaches in Turbulence[END_REF][START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF][START_REF] Fröhlich | Hybrid LES/RANS methods for the simulation of turbulent flows[END_REF] The two most popular are the Detached-Eddy Simulation (DES [START_REF] Spalart | Detached-eddy simulation[END_REF] ) and the Wall-Modeled LES (WMLES [START_REF] Sagaut | Large-eddy Simulation for Incompressible Flows[END_REF][START_REF] Larsson | Large eddy simulation with modeled wall-stress: Recent progress and future directions[END_REF][START_REF] Bose | Wall-modeled large-eddy simulation for complex turbulent flows[END_REF] ). Both methods have pros and cons but the same goal: to tackle LES limitations to high Reynolds number flows by using RANS modeling near walls. The number of required grid points, N, was estimated [START_REF] Deck | Zonal detached eddy simulation (ZDES) of a spatially developing flat plate turbulent boundary layer over the Reynolds number range 3150 ≤ Re θ ≤ 14000[END_REF][START_REF] Deck | Improvements in zonal detached eddy simulation for wall modeled large-eddy simulation[END_REF] to be N ∼ Re 2.17 for wall-resolved and wall-modeled LES of the flat plate boundary layer.

On the one hand, the DES takes advantage of the similarity of the LES and RANS equations: in the near-wall region, the RANS equations are solved, and when the mesh is fine enough, the turbulent viscosity switches to a subgrid-scale model and LES behavior is obtained. The DES was historically developed with the Spalart-Allmaras turbulence model, [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] but other models are now available. [START_REF] Strelets | Detached eddy simulation of massively separated flows[END_REF] However, the transition region between RANS and LES modeling can be a cause an-Pohlhausen integral method for laminar flows: the velocity profile is parameterized, and unknown coefficients are determined by matching boundary conditions obeying the integral boundary layer momentum equation. It allows non-equilibrium effects such as pressure gradient and convection to be included at a computing cost similar to analytical wall-models. To take into account density variations and temperature gradients, the temperature profile is also parameterized and the integral compressible boundary layer energy equation is considered. Parameterized profiles are based on the usual logarithmic wall functions with corrective terms to extend their range of validity. Instead of solving a set of differential equations as wall-models based on the thin boundary layer equation approach, a simple linear system is solved. The proposed wall-model is implemented in a finite-volume cell-centered structured grid solver and assessed on adiabatic and isothermal plane channel flows at several friction Reynolds and Mach numbers. For low Mach number cases, mean profiles, wall fluxes, and turbulent fluctuations are in agreement with those of Direct Numerical Simulation (DNS). For supersonic flows, the results are in good agreement with the DNS data, especially the mean velocity quantities and the wall friction, while standard analytical wall-models show their limits.

Electronic mail: catchirayer@cerfacs.fr Electronic mail: boussuge@cerfacs.fr D. Papadogiannis, of errors. Indeed non-physical terms are generated near the transition region. [START_REF] Hamba | Analysis of filtered Navier-Stokes equation for hybrid RANS/LES simulation[END_REF] To cure this issue and depending on the case, stochastic forcing [START_REF] Piomelli | The inner-outer layer interface in large-eddy simulations with wall-layer models[END_REF][START_REF] Keating | A dynamic stochastic forcing method as a wall-layer model for large-eddy simulation[END_REF] methods can be calibrated and used in the transition region to cancel these terms and accelerate the transition between RANS and LES.

On the other hand, WMLES solves LES equations in all the computational domain, but as opposed to a resolved LES, the mesh is coarsened in the vicinity of walls. Then, a wallmodel is used to compute the wall fluxes to take into account the effects of the inner layer part of boundary layers. Wallmodels use LES data near walls as the input. So the near-wall treatment in WMLES is very localized which eases implementation, efficient parallel computation, and application to complex geometries. However, due to a coarse mesh near walls, the numerical methods and subgrid-scale model must be carefully chosen [START_REF] Kawai | Wall-modeling in large eddy simulation: Length scales, grid resolution, and accuracy[END_REF][START_REF] Bocquet | A compressible wall model for largeeddy simulation with application to prediction of aerothermal quantities[END_REF] to avoid the development of non-physical turbulent structures. Indeed, by definition, LES [START_REF] Sagaut | Large-eddy Simulation for Incompressible Flows[END_REF][START_REF] Garnier | Large-Eddy Simulation for Compressible Flows[END_REF] solves the largest eddies and models the smallest ones, thanks to the subgrid-scale model. The frequency cutoff is generally given by the cell size of the computational grid and should be located within the inertial range of the turbulent energy spectrum. Due to this restriction, near-wall cells in resolved LES must be fine as the turbulent structures inside boundary layers are very small. [START_REF] Smits | High-Reynolds number wall turbulence[END_REF][START_REF] Jiménez | Cascades in wall-bounded turbulence[END_REF] Eddies not captured by the grid are then taken into account by the subgrid-scale model, often assuming they are isotropic. In WMLES, the local LES filter is not expected to fall in the inertial range because of the coarse near-wall cells; thus, there is no reason to accurately predict the subgrid-scale viscosity. Corrections of the near-wall subgrid-scale viscosity have been proposed by several authors. [START_REF] Bocquet | A compressible wall model for largeeddy simulation with application to prediction of aerothermal quantities[END_REF][START_REF] Cabot | Approximate wall boundary conditions in the largeeddy simulation of high Reynolds number flow[END_REF][START_REF] Wu | A constraint for the subgrid-scale stresses in the logarithmic region of high Reynolds number turbulent boundary layers: A solution to the log-layer mismatch problem[END_REF] Bocquet, Sagaut, and Jouhaud [START_REF] Bocquet | A compressible wall model for largeeddy simulation with application to prediction of aerothermal quantities[END_REF] showed in a WMLES of a plane channel flow that without correction the subgrid-scale term is highly overestimated in the first off-wall cell.

Another consequence of the use of a coarse grid near the wall in WMLES is the potential under-resolution of the turbulent streaky structures 5 located in the logarithmic region. Indeed in WMLES, inner layer effects are modeled by the wallmodel and only the outer layer is directly solved. Therefore, streaks in the viscous sublayer are not captured. The overlap (logarithmic) region is mainly composed of large-scale structures: Large-Scale Motions (LSM [START_REF] Kovasznay | Large-scale motion in the intermittent region of a turbulent boundary layer[END_REF] ) which have streamwise sizes of l x ∼ 2-3δ and the Very-Large-Scale Motions (VLSM [START_REF] Kim | Very large-scale motion in the outer layer[END_REF][START_REF] Hutchins | Evidence of very long meandering features in the logarithmic region of turbulent boundary layers[END_REF] ) whose sizes scale with l x ∼ 10δ with δ being the boundary layer thickness. Since, in the outer region, a major contribution of the turbulent kinetic energy and Reynolds stress is due to these large structures, they must be well resolved. The number of grid points per wavelength required to resolve these structures depends on the numerical methods employed. [START_REF] Kawai | Wall-modeling in large eddy simulation: Length scales, grid resolution, and accuracy[END_REF] In the wall-normal direction, the wall-model input data should be placed in the overlap region, i.e., at y 0.05-0.2δ. If the turbulent motions are under-resolved, non-physical turbulent structures are present near walls like the super-streaks encountered in DES or other hybrid RANS/LES methods. [START_REF] Larsson | Large eddy simulation with modeled wall-stress: Recent progress and future directions[END_REF][START_REF] Piomelli | The inner-outer layer interface in large-eddy simulations with wall-layer models[END_REF] However, if the grid is too fine, under-resolved near-wall streaks start to be captured and are used as the input to steady wall-model equations. [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF] Therefore, the physical and numerical modeling of the flow must be carefully chosen to accurately compute the skin friction. Otherwise the well-known log-layer mismatch appears. [START_REF] Larsson | Large eddy simulation with modeled wall-stress: Recent progress and future directions[END_REF] In the following, only WMLES is considered. Two types of wall-models exist: analytical and numerical. On the one hand, analytical wall-models are the simplest ones. Both the velocity and the temperature at a given point (usually the first off-wall point) are assumed to follow a mean solution. The treatment's cost is negligible compared to a LES time step, but the range of validity of an analytical wall-model is restricted to simple flows. Note that analytical wall-models using the standard logarithmic law of the wall 10 are also called equilibrium wall-models since this law is derived assuming a constant shear stress in the inner layer part of the turbulent boundary layer. Moreover, the logarithmic law of the wall is in theory only valid for a steady incompressible flow along a fully turbulent boundary layer at zero pressure gradient.

On the other hand, numerical wall-models solve a set of differential equations usually based on the thin boundary layer equation approximation under the RANS formalism on an implicit or separated [START_REF] Park | Numerical aspects and implementation of a two-layer zonal wall model for les of compressible turbulent flows on unstructured meshes[END_REF] grid between the wall and the wallmodel interface. Therefore, the effects of the turbulent structures located in the inner layer are taken into account globally. [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF] Numerical wall-models have a larger domain of validity than analytical wall-models since they can take into account more physical terms: advection, unsteadiness, pressure gradient, etc. The cheapest and simplest numerical wall-models (called equilibrium wall-models) only consider the wall-normal diffusion. Unsteadiness is assumed to be negligible, and pressure gradients and convection terms balance exactly. Nevertheless, this hypothesis can lead to inaccurate estimation of the wall fluxes in non-equilibrium flows. To improve the wall flux prediction, non-equilibrium wall-models are thus required. To go beyond the equilibrium assumption, non-equilibrium wall-models take also into account the pressure gradient [START_REF] Bocquet | A compressible wall model for largeeddy simulation with application to prediction of aerothermal quantities[END_REF] and possibly the time derivative. [START_REF] Hoffmann | Approximate wall boundary conditions for large eddy simulations[END_REF] Bidimensional convective terms can also be added. [START_REF] Balaras | Two-layer approximate boundary conditions for large-eddy simulations[END_REF] Most complete numerical wall-models consider the full three-dimensional unsteady RANS thin boundary layer equations. [START_REF] Cabot | Approximate wall boundary conditions in the largeeddy simulation of high Reynolds number flow[END_REF][START_REF] Wang | Dynamic wall modeling for large-eddy simulation of complex turbulent flows[END_REF][START_REF] Kawai | Dynamic non-equilibrium wall-modeling for large eddy simulation at high Reynolds numbers[END_REF][START_REF] Park | An improved dynamic non-equilibrium wall-model for large eddy simulation[END_REF] Indeed, it is important to keep all terms of the thin boundary layer equations since the pressure gradient and the convective terms are coupled. [START_REF] Hickel | A parametrized nonequilibrium wall-model for large-eddy simulations[END_REF] However, the computational cost of numerical wall-models is also higher due to the numerical resolution of complex two-or three-dimensional equations. Besides, to add non-equilibrium terms such as wallparallel diffusion or convective terms, each wall-adjacent cell should exchange data with its neighbors, reducing the parallel efficiency of the CFD solver especially for unstructured solvers. [START_REF] Park | Numerical aspects and implementation of a two-layer zonal wall model for les of compressible turbulent flows on unstructured meshes[END_REF] To take into account compressibility effects and/or temperature gradients, the current state of the art in WMLES is to use an equilibrium numerical wall-model. Indeed, such a kind of wall-model has been successfully used on

• incompressible and compressible plane channel flows with or without heat transfer; 23 • supersonic boundary layer flows; [START_REF] Kawai | Wall-modeling in large eddy simulation: Length scales, grid resolution, and accuracy[END_REF] • oblique shock/boundary-layer interaction; [START_REF] Bermejo-Moreno | Confinement effects in shock wave/turbulent boundary layer interactions through wall-modelled large-eddy simulations[END_REF] • supersonic Couette flows with heat transfer. [START_REF] Yang | Aerodynamic heating in wall-modeled large-eddy simulation of high-speed flows[END_REF] However, in flows with many walls (e.g., in turbomachinery flows), the computational overhead of such wall-models can be important. This is especially the case at very high Reynolds numbers since the computational cost of numerical wall-models is often non-linear with respect to the Reynolds number due to the increased number of degrees of freedoms required to their resolution.

In the frame of this work, the aim is to develop a wallmodel capable of taking into account the physical phenomena encountered in turbomachinery systems. Usually, in turbomachinery flows, Mach numbers are up to 1.5, [START_REF] Weber | 3-d transonic flow in a compressor cascade with shock-induced corner stall[END_REF][START_REF] Hah | Large eddy simulation of transonic flow field in nasa rotor 37[END_REF][START_REF] Ottavy | Experimental and computational methods for flow investigation in high-speed multistage compressor[END_REF][START_REF] Léonard | Steady/unsteady Reynolds-averaged Navier-Stokes and large eddy simulations of a turbine blade at high subsonic outlet mach number[END_REF][START_REF] Gourdain | Prediction of the unsteady turbulent flow in an axial compressor stage. Part 1: Comparison of unsteady RANS and LES with experiments[END_REF] Reynolds numbers [START_REF] Tyacke | Large eddy simulation for turbines: Methodologies, cost and future outlooks[END_REF][START_REF] Tucker | Eddy resolving simulations in aerospaceinvited paper (numerical fluid 2014)[END_REF] between 10 5 and 10 7 , and temperature gradients in the order of ten to hundred Kelvin. [START_REF] Arts | Aero-thermal performance of a two-dimensional highly loaded transonic turbine nozzle guide vane: A test case for inviscid and viscous flow computations[END_REF] As pointed out by Tyacke and Tucker, [START_REF] Tyacke | Future use of large eddy simulation in aero-engines[END_REF] if the flow in low-pressure turbines can be simulated with LES due to their relatively low-Reynolds number, the flows in the other components require a wallmodeling in order to reduce computational requirements at high Reynolds numbers. The interest of a wall-modeling for compressor flows is also highlighted by Gourdain et al. [START_REF] Gourdain | Large eddy simulation of flows in industrial compressors: A path from 2015 to 2035[END_REF] Recently, a hybrid wall-model, namely, integral Wall-Model for Large-Eddy Simulation (iWMLES), has been developed for adiabatic incompressible flows. [START_REF] Yang | Integral wall model for large eddy simulations of wall-bounded turbulent flows[END_REF] In iWMLES, the longitudinal velocity profile is parameterized and imposed to satisfy some boundary conditions and among them the vertically integrated momentum thin boundary layer equation. This method can include more physical phenomena than an analytical wall-model, like pressure gradients, roughness effects, or curvature effects while keeping the computational cost very low compared to a standard numerical wall-model as only a weak solution is sought. iWMLES has been applied to attached boundary layer flows with a rough surface [START_REF] Yang | Integral wall model for large eddy simulations of wall-bounded turbulent flows[END_REF] or separated flow along a smooth surface. [START_REF] Cadieux | Wallmodeled large eddy simulation of laminar and turbulent separation bubble flows[END_REF] In this study, an extension of iWMLES for compressible and isothermal flows on smooth walls is proposed. In Sec. II, the flow solver and the wall-model implementation are described. Then, in Sec. III, the iWMLES approach is explained along with its numerical resolution. Finally, in Sec. IV, the integral wall-model developed in this study is tested on quasi-incompressible and supersonic bi-periodic plane channel flows, whose Mach numbers and temperature gradients are representative of typical turbomachinery flows. Wall fluxes, velocity, temperature mean profiles, and fluctuations are compared to the Direct Numerical Simulation (DNS) data.

II. DESCRIPTION OF THE FLOW SOLVER

Even though wall-models can be seen as black boxes which compute wall fluxes using LES data such as the velocity, pressure, or temperature at some distance from the wall, they are strongly coupled with the flow solver. Indeed numerical choices must be made on how (i) the LES data are used as inputs of the wall-model; (ii) the outputs of the wall-model are used by the LES flow solver; and (iii) to take into account the coarse near-wall cells.

A. Physical and numerical modeling of the flow

The CFD solver used for this study (elsA-ONERA [START_REF] Cambier | The onera elsA CFD software: Input from research and feedback from industry[END_REF] ) is based on a cell-centered finite-volume approach to solve LES equations. The fluid is assumed to be a perfect gas, and the viscosity µ varies with the temperature T according to Sutherland's law,

µ(T ) = µ ref T T ref 3/2 T ref + S T + S , (1)
with Sutherland's constant S = 110.4 K. µ ref and T ref are, respectively, the reference viscosity and the reference temperature.

Compressible LES equations are written under the Favre averaging and on implicit mesh filtering formalism. In the following, RANS and LES variables would be denoted, respectively, with a bar (f ) and a tilde ( f ).

The subgrid-scale model is estimated by the Wall-Adapting Local Eddy-viscosity (WALE [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF] ) subgrid-scale model which is based on the velocity gradient,

ν sgs = (C w ∆) 2 (S d ij S d ij ) 3/2 ( S ij S ij ) 5/2 + (S d ij S d ij ) 5/4 , (2)
with the constant being C = 0.5, ∆ being the cell volume, u i being the velocity components, S ij = (∂ u i /∂x j + ∂ u j /∂x i )/2 being the components of the strain rate tensor, and

S d ij = 1 2 ∂ u i ∂x l ∂ u l ∂x j + ∂ u j ∂x l ∂ u l ∂x i + 1 3 ∂ u m ∂x l ∂ u l ∂x m δ ij , (3)
where δ ij is the Kronecker symbol.

In the finite volume framework, thanks to the Gauss divergence theorem, the divergence of the convective and diffusive terms of the Navier-Stokes equations is reduced to a sum of fluxes. Convective fluxes are estimated using a second-order centered scheme without artificial dissipation in the skewsymmetric form, which is observed to reduce aliasing errors. [START_REF] Kravchenko | On the effect of numerical errors in large eddy simulations of turbulent flows[END_REF] To determine the diffusive fluxes, the laminar viscosity, the subgrid-scale model, the velocity gradient, and the temperature gradient at each interface between cells are required. The laminar viscosity and subgrid-scale model at the interface are evaluated by averaging the two adjacent cells' values. Velocity and temperature gradients are computed using the Gauss divergence theorem on a shifted control volume (centered on the interface). This discretization leads to a three-point stencil second-order scheme. At walls, two fictitious cells are introduced in order to keep the same schemes.

When dealing with moderate or high Reynolds number flows, the convective fluxes are expected to be dominant over the diffusive fluxes, except in the walls' vicinities. However, in WMLES, the velocity and temperature gradients' discretization is erroneous at wall interfaces since the near-wall coarse mesh used cannot capture high gradients. Therefore, velocity and temperature gradients at each wall interface are evaluated using a wall-model as explained in Subsection II B.

Time advancement is performed explicitly using a lowstorage Runge-Kutta scheme with four stages. For a differential equation ∂u/∂t = f (u, t), the steps are given by [START_REF] Sagaut | Large-eddy Simulation for Incompressible Flows[END_REF] , u n+1 = u n + R [START_REF] Garnier | Large-Eddy Simulation for Compressible Flows[END_REF] , (4) with ∆t being the time step and R (i) being the intermediate residuals defined by

                         u (0) = u n , u (1) = u n + ∆t 4 R (0) , u (2) = u n + ∆t 3 R (1) , u (3) = u n + ∆t 2 R
                   R (0) = f (u (0) , t n ),
R [START_REF] Gourdain | RANS and LES for the heat transfer prediction in turbine guide vane[END_REF] = f (u [START_REF] Gourdain | RANS and LES for the heat transfer prediction in turbine guide vane[END_REF] , t n + ∆t 4 ), R [START_REF] Sagaut | Large-eddy Simulation for Incompressible Flows[END_REF] = f (u [START_REF] Sagaut | Large-eddy Simulation for Incompressible Flows[END_REF] , t n + ∆t 3 ), R [START_REF] Garnier | Large-Eddy Simulation for Compressible Flows[END_REF] = f (u [START_REF] Garnier | Large-Eddy Simulation for Compressible Flows[END_REF] , t n + ∆t 2 ).

(

) 5
This scheme is fourth-order if f is linear, second-order otherwise. In all the following simulations, a constant time step is chosen such that the Courant-Friedrichs-Lewy (CFL) number is around 0.7. This numerical framework has already been proposed and validated, [START_REF] Bocquet | A compressible wall model for largeeddy simulation with application to prediction of aerothermal quantities[END_REF] except here the WALE subgrid-scale model is employed instead of the selective Smagorinsky model to better predict supersonic flows.

B. Relation between the flow solver and the wall-model

Due to the coarse mesh used in WMLES, high velocity and temperature gradients at walls cannot be properly captured. That is why at each time step, during the viscous flux computation, the wall fluxes are computed by the wall-model. Given information from the LES field at the matching point (the first off-wall cell in the frame of this work) like the velocity vector, pressure, density, or temperature, the wall friction vector τ w and the wall heat flux φ w are returned by the wallmodel. Note that τ w and φ w are not exactly LES variables as they are the outputs of RANS-based wall-model equations.

The wall friction vector τ w is assumed to be aligned with the velocity vector at the first off-wall point projected in the plane parallel to the wall, u 1 . Therefore, only two-dimensional equations can be considered in wall-models, reducing their complexity and computational cost. The projection step is briefly reminded here. At each call of the wall-model, the velocity vector at the first off-wall point u 1 is projected in the plane parallel to the wall as shown in Fig. 1 the wall-normal vector, y 1 being the distance between the wall and the first off-wall cell, and t = u 1 / u 1 being the local velocity direction. The wall-model equations are solved in the local boundary layer frame (t, n), and the wall friction modulus τ w and the wall heat flux are returned to the LES solver. As τ w is assumed to be aligned with u 1 , the wall friction vector coordinates in the flow solver frame can be computed and used as a boundary condition by the flow solver. It can be noticed that from the wall-model point of view, the boundary layer is always attached.

C. Analytical wall-model

A fair comparison between wall-models can only be made in a given flow solver since there is a strong coupling between the two. Therefore, a standard analytical wall-model is implemented in order to assess the wider range of validity of numerical wall-models with respect to analytical ones. It is based on the turbulent boundary layer equations under zero pressure gradient for an incompressible flow. Turbulence is assumed to be statistically bidimensional and the mean flow steady. Under these hypotheses, the wall-normal momentum equation implies that the wall-normal variation of pressure is of the order of the boundary layer thickness. So this variation can be neglected with respect to the streamwise pressure gradient, and the pressure is assumed to be constant along the wall-normal direction. Moreover, by an order-of-magnitude analysis, it can be shown that in the limit of infinite Reynolds number, the shear stress τ is constant in the inner layer part of the boundary layer and therefore equal to the wall shear stress τ w . This result is the key point to justify the use of wall-laws because a single velocity scale u τ = τ w /ρ w in the whole inner layer can be defined and a velocity profile u + = u/u τ , independent of the outer layer, can be sought. With this velocity scale, a linear profile u + = y + in the linear sublayer and a logarithmic profile u + = (1/κ) ln y + + B in the inertial sublayer are found. Similarly, for an isothermal wall, the heat flux φ is constant and equal to the wall heat flux φ w in the inner layer region and a temperature scale T τ = φ w /(ρ w c p u τ ) can be defined with c p being the specific heat capacity at constant pressure and ρ w being the wall density. A self-similar temperature profile T + = (T -T w )/T τ following a linear law in the linear sublayer and a logarithmic law in the inertial sublayer is found, where T is the wall temperature.

Based on these results, Reichardt [START_REF] Reichardt | Complete representation of the turbulent velocity distribution in smooth pipes[END_REF] and Kader [START_REF] Kader | Temperature and concentration profiles in fully turbulent boundary layers[END_REF] proposed smooth velocity and temperature wall-laws in the whole inner layer. They provide an explicit relation between the mean velocity and temperature and the wall fluxes τ w and φ w . From these two laws, an analytical wall-model [START_REF] Bocquet | A compressible wall model for largeeddy simulation with application to prediction of aerothermal quantities[END_REF][START_REF] Bras | Development of compressible large-eddy simulations combining high-order schemes and wall modeling[END_REF] can be built,

u + 1 = u + 1 = u 1 u τ = 1 κ ln(1 + κy + 1 ) + (B - 1 κ ln κ) × (1 -exp( -y + 1 11 ) - y + 1 11 exp( -y + 1 3 )), (6)
T + 1 = T + 1 = T 1 -T w T τ = Pry + 1 exp(-Γ) + ( 1 κ t ln(1 + y + 1 ) + β) exp(-1/Γ), (7)
with Pr being the Prandtl number and κ = 0.41, (8a)

B = 5.25, (8b)
κ t = 1 2.12 , (8c)
Γ = 10 -2 (Pry + 1 ) 4 1 + 5Pr 3 y + 1 , (8d)
β = (3.85Pr 1/3 -1.3) 2 + 1 κ t ln(Pr). ( 8e
)
As explained in Subsection II B, the bidimensional unsteady velocity u 1 and temperature T 1 at the first offwall cell are imposed to follow a mean wall-law at each time step. Given LES data at the first off-wall point, the following procedure is applied in order to compute the wall fluxes:

1. The pressure is assumed to be constant in the wall-normal direction, allowing the wall density to be determined, thanks to the perfect gas relation ρ w = p 1 /(rT w ) with p 1 being the pressure at the first point and r being the perfect gas constant. In the case of an adiabatic boundary condition, the wall temperature is supposed to be equal to the friction temperature

T τ = T e (1+Pr 1/3 (γ -1)/2 M e 2
). As the values at the outer edge of the boundary layer like temperature T e and Mach number M e = u e / γr T e are not accessible in complex geometries, their values at the first off-wall cell are taken. 2. Before being used as the input by the wall-model, the velocity vector is projected along the plane parallel to the wall as explained in Subsection II B. 3. The friction velocity u τ = τ w /ρ w is computed, thanks to Eq. ( 6) and the Newton-Raphson algorithm. 4. For an isothermal wall, the friction temperature T τ = -φ w /(ρ w c p u τ ) is computed, thanks to Eq. ( 7) and the Newton-Raphson algorithm. 5. The wall friction vector components in the flow solver frame are computed by switching from the boundary layer frame to the flow solver frame.

In the following, this analytical wall-model will be denoted as aWMLES.

III. COMPRESSIBLE INTEGRAL WMLES APPROACH

In this section, an integral wall-model for the LES (iWM-LES) of compressible flows along adiabatic or isothermal smooth plates is presented. First, the key idea of the wall-model is described. As it is based on the integral thin boundary layer equations, these equations are formulated before explaining how to solve them. Finally, the overall wall-model resolution is summarized.

A. Description of the iWMLES approach

For a laminar boundary layer, Kármán [START_REF] Kármán | Über laminare und turbulente reibung[END_REF] and Pohlhausen [START_REF] Pohlhausen | Zur näherungsweisen integration der differentialgleichung der iaminaren grenzschicht[END_REF] proposed a method to estimate the local velocity profile without solving differential equations. Instead of directly solving the boundary layer equations, the longitudinal velocity profile is approximated by a fourth-order polynomial. Then, the only unknowns to fully determine the velocity profile are the polynomial coefficients. They are computed to impose the local longitudinal velocity to verify physical boundary conditions and among them the vertically integrated boundary layer momentum equations. Thus, a velocity profile taking into account the pressure gradient is built at a very low computational cost. In the case of the laminar boundary layer with zero pressure gradient, the friction coefficient computed presents 2% of error compared to the exact Blasius solution.

Here, a formulation of iWMLES for compressible and isothermal flows on smooth walls is proposed and assessed. This wall-model is analog to the von Kármán-Pohlhausen integral method except the boundary layer is assumed to be fully turbulent. To handle these kinds of flows, the velocity and the temperature profiles are parameterized. Then, the unknown coefficients are computed so that velocity and temperature profiles satisfy some boundary conditions and among them the vertically integrated bidimensional turbulent boundary layer momentum and energy equations. Therefore the computed profiles represent a weak solution of the integral boundary layer equations and take into account more physics than analytical wall-models. Instead of solving directly the thin boundary layer equations and using expensive numerical algorithms, a weak solution is sought and a simple scalar equation system has to be solved. So iWMLES can be seen as a hybrid wallmodel: either an analytical wall-model with dynamic coefficients or a numerical wall-model where only a weak solution is searched.

B. Integral compressible thin boundary layer equations

In this subsection, the vertically integrated boundary layer momentum and energy equations required in iWMLES are formulated. In the following, the longitudinal and normal velocities are, respectively, u and , and the longitudinal and wall-normal directions are, respectively, the x and y axes. The other variables follow conventional nomenclature. As opposed to the original formulation of iWMLES for incompressible flows, [START_REF] Yang | Integral wall model for large eddy simulations of wall-bounded turbulent flows[END_REF] Pr is the molecular thermal conductivity. µ t and λ t are the turbulent viscosity and the thermal conductivity, respectively, and is the fluctuating part of the wall-normal velocity. Equation (9a) is related to the conservation of mass, and Eqs. (9b) and (9c) are the longitudinal and wall-normal momentum conservation equations. Equation (9d) is the energy equation, and the perfect gas law (9e) closes the equation system. In the above equations, the diffusion term along the longitudinal direction x is neglected.

In the following, for the sake of simplicity, any vertically integrated quantities f between y = 0 and y = y 1 will be denoted as

L f = y 1 0 f dy.
According to Eq. (9c), the wall-normal pressure gradient is of the order of magnitude of the boundary layer thickness and therefore is small compared to the streamwise pressure gradient. This involves an important result in the boundary layer theory: the pressure can be assumed to be constant along the wall-normal direction, which means that it is an input of the wall-model extracted from the LES field. Taking advantage of this approximation and the perfect gas law (9e), Eqs. (9a), (9b), and (9d) are integrated along the wallnormal direction between the wall (assumed without the loss of generality to be located at y = 0) and the first off-wall point at y = y 1 ,

ρ 1 v 1 = - y 1 0 ∂ ρ ∂t + ∂ ρu ∂x dy, ( 10a
)
∂L ρu ∂t -u 1 ∂L ρ ∂t + M x = τ 1 -τ w , ( 10b
)
c p r -1 ∂ p 1 ∂t y 1 -c p T 1 ∂L ρ ∂t + M T x = L τ ∂u ∂y -φ 1 + φ w , (10c)
with u 1 , T 1 , and p 1 being the LES data at y 1 . The convective terms are expressed as follows:

M x = ∂ p 1 ∂x y 1 + ∂L ρu 2 ∂x -u 1 ∂L ρu ∂x , (11)
M T x = c p r p 1 ∂L u ∂x -c p T 1 ∂L ρu ∂x - ∂ p 1 ∂x L u . ( 12
)
Scalar equations (10b) and (10c) are, respectively, the vertically integrated momentum and energy equations. Velocity and temperature profiles in iWMLES will be constrained to be solutions of these equations. Note that these equations require only the values of turbulence models at the first offwall cell and of the vertically integrated shear stress. Thus, iWMLES is less sensitive to turbulence models than numerical wall-models. Indeed, in numerical wall-models, turbulence models are involved over the entire wall-normal direction and need to be carefully chosen. [START_REF] Bocquet | A compressible wall model for largeeddy simulation with application to prediction of aerothermal quantities[END_REF][START_REF] Cabot | Approximate wall boundary conditions in the largeeddy simulation of high Reynolds number flow[END_REF][START_REF] Wang | Dynamic wall modeling for large-eddy simulation of complex turbulent flows[END_REF][START_REF] Kawai | Dynamic non-equilibrium wall-modeling for large eddy simulation at high Reynolds numbers[END_REF] Therefore, standard mixing-length models are used,

µ t (y) = ρ(κyVD 26 (y)) 2 | ∂u ∂y |(y), (13)
λ t (y) = ρc p (κy) 2 VD 26 (y)VD 35 (y)| ∂u ∂y |(y), (14)
with VD being the Van-Driest damping function defined by

VD a (y) = 1 -exp -y + 3 a , ∀a ∈ R * , ( 15
)
y + 3 (y) = y ρ τ w /ρ µ = yu τ ρ w ρ µ . ( 16
)
An y + 3 scaling (also denoted as y * in the literature) is used to improve compressibility effect prediction in Eqs. ( 13) and (14). Using this scaling, a better prediction of the wall fluxes on supersonic boundary layers is obtained [START_REF] Aupoix | Compressible turbulent boundary layer modelling[END_REF] for either adiabatic or isothermal walls. Better behavior is also observed for velocity fluctuations and the Reynolds stress [START_REF] Coleman | A numerical study of turbulent supersonic isothermal-wall channel flow[END_REF][START_REF] Huang | Compressible turbulent channel flows: DNS results and modeling[END_REF] with the y + 3 scaling compared to y + .

C. Parameterized profiles

Instead of directly solving Eqs. (10b) and (10c), a weak solution is sought by parameterizing velocity and temperature profiles. Analytical results show that for a steady incompressible fully turbulent boundary layer flow with zero pressure gradient, the longitudinal velocity follows a logarithmic law. Even if the existence and the universality of this law are still debated, [START_REF] Nagib | Variations of von Kármán coefficient in canonical flows[END_REF][START_REF] George | Is there a universal log law for turbulent wall-bounded flows?[END_REF] it is observed to hold even for low Mach number flows with or without pressure gradients [START_REF] Houra | Effects of adverse pressure gradient on heat transfer mechanism in thermal boundary layer[END_REF][START_REF] Johnstone | The resilience of the logarithmic law to pressure gradients: Evidence from DNS[END_REF] or curvature effects. [START_REF] Gillis | Turbulent boundary-layer flow and structure on a convex wall and its redevelopment on a flat wall[END_REF] Thus, the velocity and temperature profiles are parameterized based on the logarithmic laws of the wall defined in Eqs. ( 6) and (7) with corrective terms,

u(y; u τ ; δ ν ; A) = u τ 1 κ ln(1 + κy + ) + (B - 1 κ ln κ) × (1 -exp( -y + 11 ) - y + 11 exp( -y + 3
))

+ Ay + (1 -exp( -y + δ ν )) , (17)
             T (y; u τ ; T τ ; δ c ; A T ) = T w + T τ Pry + exp(-Γ) + ( 1 κ t ln(1 + y + ) + β) exp(-1/Γ) + A T y + (1 -exp( -y + δ c
)) for isothermal wall,

T (y; u) = T 1 1 + Pr 1/3 γ-1 2 M 2 1 (1 -u 2 u 2 1
) for adiabatic wall,

with M 1 = u 1 / γr T 1 being the local Mach number based on the velocity projected in the plane parallel to the wall ( u 1 ) and A, A T , δ c , δ ν , u τ , and T τ being six scalar unknowns.

The correction terms A and A T aim to simulate the deviation from the logarithmic law and build a local dynamic wall-law which satisfies the vertically integrated momentum and energy equations (10b) and (10c). A linear correction term is added to the temperature profile, and a square root correction term is added to the velocity profile as previous studies have shown that a square root region appears when the boundary layer is not at equilibrium. [START_REF]Analysis of a turbulent boundary layer subjected to a strong adverse pressure gradient[END_REF][START_REF] Skote | Direct numerical simulation of a separated turbulent boundary layer[END_REF] Since the viscous sublayer seems to be only slightly affected by non-equilibrium effects, correction terms are applied only in the inertial sublayer, thanks to the Van-Driest damping function 1exp(-y + /δ ν ).

For an adiabatic boundary condition, the wall heat flux is zero and only the wall friction vector needs to be determined by the wall-model. Moreover, the temperature profile does not follow the logarithmic law of the wall. Thus, in this case, the temperature is determined from Walz's law. [START_REF] Walz | translation in Boundary Layers of Flow and Temperature[END_REF] In theory, this relation needs the values at the outer edge of the boundary layer, but these quantities are not easily accessible in complex geometries or in unstructured solvers. Local values at first offwall cell are then used, and this approximation will be validated a posteriori. Besides, as the wall heat flux does not need to be computed, the vertically integrated energy equation (10c) is not solved as it will be explained in Subsection III D.

As the y + scaling does not take into account the wake region, selected profiles and turbulence models ( 13) and ( 14) are only valid in the inner layer of the boundary layer. From a physical viewpoint, improved stability and accuracy are obtained when the wall-model interface is located in the overlap region between the inner and outer layers. [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF][START_REF] Larsson | Large eddy simulation with modeled wall-stress: Recent progress and future directions[END_REF] Thus, the wall-model interface is aimed to be at about 5%-15% of the boundary layer thickness and y + 50. This is also the case of the analytical wall-model defined by Eqs. ( 6) and (7).

D. Evaluation of the profile parameters: Selected boundary conditions

To completely define the velocity and temperature profiles given by Eqs. ( 17) and ( 18), three scalar unknowns must be computed in the case of an adiabatic wall: u τ , δ ν , and A. If the wall is considered to be isothermal, there are three additional scalar unknowns: T τ , δ c , and A T . They are determined according to boundary conditions and vertically integrated boundary layer equations:

1. the velocity profile defined in Eq. ( 17) must match with the LES velocity field at y 1 in order to be continuous,

u 1 (y 1 ; u τ ; δ ν ; A) = u 1 ; (19)
2. the diffusive layer thickness δ ν is supposed to be constant,

δ ν = 11, (20)
and this value corresponds to the intersection point between the linear law u + = y + and the logarithmic law of the wall in the case of an incompressible flow; 3. the velocity profile defined in Eq. ( 17) must be the solution of the vertically integrated momentum equation (10b); 4. in the case of an isothermal boundary condition, three additional unknowns A T , δ c , and T τ must be computed: (a) the temperature profile defined in Eq. ( 18) must match with the LES temperature field at y 1 in order to be continuous,

T 1 (y 1 ; u τ ; T τ ; δ c ; A T ) = T 1 , (21)
(b) as with the diffusive layer thickness δ ν , the convective layer thickness δ c is supposed to be constant and

δ c = 11, (22)
(c) the velocity and temperature profiles defined in Eqs. ( 17) and ( 18) must be the solution of the vertically integrated energy equation (10c).

All these equations form a system of three or six (depending on whether the wall is supposed to be adiabatic or isothermal) scalar equations,

                               u τ 1 κ ln(1 + κy + 1 ) + (B - 1 κ ln κ)(1 -exp( -y + 1 11 ) - y + 1 11 exp( -y + 1 3 )) + Ay + 1 (1 -exp( -y + 1 δ ν )) = u 1 , δ ν = 11, ∂L ρu ∂t -u 1 ∂L ρ ∂t + M x = τ 1 -τ w , (23a)
and if the boundary condition is isothermal,

                             T w + T τ Pry + 1 exp(-Γ) + ( 1 κ t ln(1 + y + 1 ) + β) exp(-1/Γ) + A T y + 1 (1 -exp( -y + 1 δ c )) = T 1 , δ c = 11, c p r -1 ∂ p 1 ∂t y 1 -c p T 1 ∂L ρ ∂t + M T x = L τ ∂u ∂y -φ 1 + φ w , (23b)
which is much easier to solve than numerical wall-model differential equations.

For an adiabatic boundary condition, the list of unknowns is reduced to three as only the wall friction vector has to be computed by the wall-model. Indeed in this case, the friction temperature T τ is zero as well as the wall heat flux. As explained in Subsection III C, Walz's law [START_REF] Walz | translation in Boundary Layers of Flow and Temperature[END_REF] is used to link the temperature to the velocity profile and the correction term A T is not required. Therefore, in this case, only the velocity profile needs to be determined with the friction velocity u τ , the diffusive layer thickness δ ν , and the correction term A.

E. Numerical resolution of the integral equations

The method to solve the system of Eqs. (23) written in Subsection III D is now detailed in the case of an isothermal wall. For an adiabatic wall, the resolution is simpler as the temperature profile is determined from the velocity profile using Eq. (18). In this section, variables at the current time step will be denoted with a superscript n and n -1 for those at the previous time step.

At each wall interface, using LES data located at the first off-wall cell y 1 , thanks to Eqs. ( 19)-( 22), the scalar unknowns A n , δ ν , A n T , and δ c are expressed as functions of u n τ and T n τ ,

                                       A n (u n τ ; T n τ ) = u n 1 u n τ -1 κ ln(1 + κy + 1 n ) + (B -1 κ ln κ)(1 -exp( -y + 1 n 11 ) - y + 1 n 11 exp( -y + 1 n 3 )) y + 1 n (1 -exp( -y + 1 n δ ν ))
,

δ ν = 11, A n T (u n τ ; T n τ ) = T n 1 -T w T n τ -Pry + 1 n exp(-Γ) + ( 1 κ t ln(1 + y + 1 n ) + β) exp(-1/Γ) y + 1 n (1 -exp( -y + 1 n δ c ))
,

δ c = 11. (24)
Thus, system ( 23) is reduced to a non-linear system of two equations [(10b) and (10c)] with two unknowns: u n τ and T n τ . The temporal derivatives in these equations are discretized with an explicit scheme and the spatial derivatives with a second-order centered scheme. An upwind scheme is used if the wall interface is adjacent to another boundary condition. Then Eqs. (10b) and (10c) become

                 L n ρu (u n τ ; T n τ ) -u n-1 1 L n ρ (u n τ ; T n τ ) = L n-1 ρu (1 -u n-1 1 ) + ∆t(-M n-1 x + τ 1 n-1 -τ w n-1 ), (25a)
c p T n-1 1 L n ρ (u n τ ; T n τ ) = c p r -1 y 1 ( p n 1 -p n-1 1 ) + c p T n-1 1 L n-1 ρ , -∆t(-M n-1 T x + L n-1 τ ∂u ∂y -φ 1 n-1 + φ w n-1 ), (25b)
where the terms depending on the unknowns u n τ and T n τ are specified and with ∆t being the time step of the LES solver. Given data located at the first off-wall point y 1 at the current and previous time steps along with the wall fluxes at previous iteration u n-1 τ and T n-1 τ , u n τ and T n τ are computed by solving Eqs. (25a) and (25b), thanks to bidimensional Newton-Raphson's algorithm whose corresponding Jacobian matrix is given in Appendix A. u n-1 τ and T n-1 τ are chosen as initial guesses of the algorithm. At the first iteration of the flow solver, i.e., at n = 0, wall fluxes are estimated by a first-order finite difference scheme like in wall-resolved LES.

Integral terms L f in Eqs. (25a) and (25b) cannot be determined analytically unless the density is constant. [START_REF] Yang | Integral wall model for large eddy simulations of wall-bounded turbulent flows[END_REF] Besides, since the pressure is assumed to be constant along the wallnormal direction, the density profile can be determined from the temperature profile using the perfect gas law. Thus, a Gauss-Legendre quadrature is used to compute the L f terms: definite integrals are approximated as a weighted sum of function values at N points. For a function f : [0; y 1 ] → R, the quadrature implies

y 1 0 f (y) dy y 1 2 N i=0 w i f y 1 2 (1 + ξ i ) , (26)
with (w i ) 0≤i ≤N and (ξ i ) 0≤i ≤N , respectively, being the weights and the Gauss nodes. In the case of the Gauss-Legendre quadrature, i and ξ i are chosen in such a way that the quadrature is exact for Legendre polynomials of degree less than or equal to 2N -1. Figure 2 shows the observed computational overhead of iWMLES and an equilibrium numerical wall-model [START_REF] Bocquet | A compressible wall model for largeeddy simulation with application to prediction of aerothermal quantities[END_REF] (denoted as TBLE in the following) with respect to an analytical model aWMLES in the cases of an isothermal or adiabatic wall. The iWMLES computational cost is much lower compared to the TBLE one, although the former takes into account FIG. 2. Ratio of the observed computational cost T wall-model of the iWMLES and a numerical wall-model for LES (TBLE) based on the equilibrium thin boundary layer equations with respect to the cost T aWMLES of an analytical wall-model per cell interface for adiabatic and isothermal walls. For iWMLES, corresponding linear regression is in dashed lines. FIG. 3. Scheme of the iWMLES resolution algorithm. convective terms. Besides, and as opposed to TBLE, the iWM-LES algorithm scales linearly with the number of Gauss nodes N used. The iWMLES algorithm is nearly twice as fast when the wall is assumed to be adiabatic instead of isothermal. Indeed, in this case, the integral energy equation is not solved and a velocity-temperature relation is used. Such improvement is also observed for TBLE but in a reduced proportion as the temperature profile still needs to be computed. In the following, N is taken equal to 10.

The iWMLES resolution algorithm is illustrated in Fig. 3 and can be summarized as follows:

1. At the first iteration, wall fluxes are computed by a nearest-cell approximation as in a standard resolved LES. 2. Then, at each wall interface, the LES data at matching point y 1 at the previous and current iterations n along with the wall fluxes at previous iteration n -1 are used as the input to the iWMLES model in order to compute the wall fluxes at the current iteration. (a) As with the analytical wall-model described in Subsection II C, wall density ρ n w is computed from the perfect gas relation. (b) The velocity at the current iteration n is projected along the plane parallel to the wall as explained in Subsection II B. (c) At each wall interface and given data at the previous solver iteration, the scalar parameters A n-1 , A n-1 T , δ ν , and δ c are computed from u n-1 τ and T n-1 τ , thanks to Eq. ( 24) in order to define the velocity and temperature profiles given by Eqs. ( 17) and ( 18). (d) These profiles are used to compute the right-hand side terms of Eqs. (25a) and (25b). Integral terms L n-1 f are approximated with a Gauss-Legendre quadrature. Note that the density profile is determined from the temperature profile using the perfect gas relation Eq. (9e). Convective terms M n-1

x and M n-1 T x are computed using a second-order centered scheme. Knowing the form of the velocity and temperature profiles, their wall-normal derivatives can be analytically derived. They are given in Appendix B. Then, turbulent viscosity µ n-1 t and thermal conductivity λ n-1 t are given by Eqs. ( 13) and ( 14).

IV. APPLICATION TO BI-PERIODIC PLANE CHANNEL FLOWS

The integral wall-model for compressible flows described in Sec. III is assessed on plane channel flows, and results are compared to the available DNS data. After presenting the selected DNS configurations, the numerical setup used to launch WMLES computation is described. Finally, the results on quasi-incompressible and supersonic cases are presented.

A. Test cases

In order to evaluate the capacity of the integral wallmodel to accurately predict wall fluxes, several DNS of biperiodic plane channel flows are considered as summarized in Table I.

For incompressible DNS, several friction Reynolds numbers Re τ are available (1020, [START_REF] Abe | Surface heat-flux fluctuations in a turbulent channel flow up to Re τ = 1020 with Pr = 0.025 and 0.71[END_REF] 2003, [START_REF] Hoyas | Scaling of the velocity fluctuations in turbulent channels up to Re τ = 2003[END_REF] 4179, [START_REF] Lozano-Durán | Effect of the computational domain on direct simulations of turbulent channels up to Re τ = 4200[END_REF] and 5186 [START_REF] Lee | Direct numerical simulation of turbulent channel flow up to Re τ 5200[END_REF] ). For the Re τ = 1020 DNS, a temperature gradient is imposed by a constant and uniform heat flux at walls. To reproduce this DNS in WMLES, a constant wall temperature T is applied and the flow is forced such that T b /T = 1.1 with T b = ∫ 2h 0 ρuT dy / ∫ 2h 0 ρu dy being the bulk temperature and h being the half-height channel. In the other DNS, the wall

with Re f = 2hρ f u b /µ f being the Reynolds number based on the film density ρ f = 2ρ ρ b /(ρ + ρ b ), the film viscosity µ f computed from Sutherland's law, and the film temperature For all plane channel flow simulations, the flow is forced to prevent the decrease in kinetic energy due to dissipation. The compressible DNS are reproduced here by applying a constant wall temperature T and forcing the bulk temperature T b to get the same ratio T b /T as in the corresponding DNS.

T f = (T b + T )/2. ρ b = ∫
In all cases, it has to be noted that the friction Reynolds number is an output of the computation as it depends on the friction velocity u τ . It is the bulk Reynolds number Re b = ρ b u b h/µ b which defines the flow, where µ b is the bulk viscosity obtained from Sutherland's law and the bulk temperature T b . The wall heat flux coefficient B q is defined by B q = φ w ρ w c p u τ T w and the friction Mach number M τ by M τ = u τ / √ γrT w .

B. Simulation setups

The cell-centered finite volume flow solver elsA-ONERA described in Subsection II A is used to carry out all the simulations. In all adiabatic WMLES cases, the Prandtl number Pr is set to 0.72. In the other cases, Pr is chosen according to the DNS values: at Re τ = 1020, Pr is set at 0.71, and in all supersonic cases, Pr = 0.7. For all cases, the turbulent Prandtl number Pr t is assumed constant and equal to 0.9. In fact, using Germano-like dynamic calibration, one can observe that this parameter is not constant and is spatially varying. In practice, it appears to have a minor influence only. [START_REF] Li | Revisiting the subgrid-scale Prandtl number for large-eddy simulation[END_REF][START_REF] Moin | A dynamic subgrid-scale model for compressible turbulence and scalar transport[END_REF] The computational domain's size is defined by (x, y, z) ∈ Ω = [0; 2πδ] × [0; 2δ] × [0; πδ] with δ being the boundary layer thickness which is equal to the channel half-height h and x, y, and z, respectively, being the coordinates in the longitudinal, wall-normal, and spanwise directions. For high-Reynolds number flows like those considered in this study, the computational domain size is big enough to capture turbulent eddies and has been used in previous studies. [START_REF] Bocquet | A compressible wall model for largeeddy simulation with application to prediction of aerothermal quantities[END_REF][START_REF] Bras | Development of compressible large-eddy simulations combining high-order schemes and wall modeling[END_REF][START_REF] De Wiart | Implicit les of free and wall-bounded turbulent flows based on the discontinuous galerkin/symmetric interior penalty method[END_REF] Two uniform grids G1 and G2 are considered: 25 × 21 × 21 and 49 × 41 × 41. On the former mesh, the first off-wall cell is located at y 1 = 0.05δ, while on the latter, y 1 = 0.025δ. Since cell sizes are approximately given by ∆x = 5.2y 1 and ∆z = πy 1 , the cell aspect ratio follows the outer-layer turbulent structure's size. All cases' parameters are summarized in Tables II-IV.

In order to enforce the correct mass flow rate and bulk temperature, energy is injected in the domain through two volumic source terms S x and Q x which are, respectively, added to the streamwise momentum and energy equations, S x aims to force the flow to reach a given bulk Reynolds ber by imposing the correct mass flow rate (ρu) b target , while Q x forces the bulk temperature. These two terms are computed at each iteration of the flow solver. Constant parts can be added to S x and Q x in order to accelerate flow convergence. [START_REF] Cabrit | Modélisation des flux pariétaux sur les tuyères des moteurs à propergol solide[END_REF] The source terms are also taken into account in iWMLES equations (25a) and (25b). Moreover, numerical experiments showed that identical results are obtained when M x and M T x are taken into account if the wall-model input data are filtered in time. [START_REF] Yang | Integral wall model for large eddy simulations of wall-bounded turbulent flows[END_REF][START_REF] Meneveau | A Lagrangian dynamic subgridscale model of turbulence[END_REF] Indeed in plane channel flows, convective terms are statistically zero, but due to the unsteady nature of LES, non-filtered M x and M T x terms are always oscillating even after reaching flow convergence. In order to reduce the computational cost of the simulations, convective terms are neglected in the following.

S x (t) = (ρu) b target - 1 Ω Ω ρu dΩ, (28)
Q x (t) = T b target - ∫ Ω ρuT dΩ ∫ Ω ρu dΩ + uS x . (29)
The initial condition is built from a steady state based on a power-law velocity profile,

∀(x, y, z) ∈ Ω,                        u(x, y, z, t = 0) = 8 7 u b target 1 -|1 - y h | v(x, y, z, t = 0) = 0 w(x, y, z, t = 0) = 0 T (x, y, z, t = 0) = T w , (30)
which is perturbed with an additional random noise applied to the velocity field. The noise maximum amplitude is equal to 10% in the longitudinal direction and 5% in the other directions of the bulk velocity. It can be noted that random noise is enough to force turbulent transition of the flow because coarse grids are considered in this study. For finer grids as encountered in resolved LES, vortex rings are added to the initial flow to mimic turbulent structures. [START_REF] Bras | Development of compressible large-eddy simulations combining high-order schemes and wall modeling[END_REF][START_REF] Kremer | Étude numérique d'écoulements de paroi compressibles: Méthodes d'intégration temporelle semi-implicites et application au canal plan turbulent[END_REF] However, WMLES grids are not fine enough to discretize these kinds of structures.

For all simulations, flow convergence toward a statistically steady state is assumed to be reached once the kinetic energy's variation in the whole domain does not exceed 0.01% of its mean value. Then, statistics are obtained by time averaging the FIG. 4. Mean velocity profiles in wall units for adiabatic quasiincompressible cases. Profiles are shifted by multiples of 10. From bottom to top: G1 M0.2 R2003 adia, G2 M0.2 R2003 adia, G1 M0.2 R4179 adia, G2 M0.2 R4179 adia, G1 M0.2 R5186 adia, G2 M0.2 R5186 adia.

flow over 200 × 2πδ/u b periods through the channel. Space averaging along the longitudinal and spanwise directions is also done before comparing the mean profiles with DNS data.

C. Results for quasi-incompressible flows

Before evaluating the capacity of iWMLES to take into account complex physical phenomena such as compressibility effects, it is required to check whether iWMLES is able to recover standard results where aWMLES is known to be accurate. As explained in Subsection IV A, four incompressible DNS are considered. Besides, a small temperature gradient is also applied on all adiabatic cases in order to evaluate iWMLES on quasi-incompressible isothermal flows. Therefore, seven quasi-incompressible cases in terms of Reynolds numbers and wall boundary conditions are computed with aWMLES and iWMLES. All numerical and physical case parameters are summarized in Tables II andIII.

Mean profiles in wall units for adiabatic and isothermal cases are, respectively, shown in Figs. 4 and5. In all cases, iWMLES recovers the logarithmic law of the wall and is in agreement with the DNS profiles, which means that the corrective terms added in the parameterized velocity and temperature profiles cancel out for different Reynolds numbers and grids. Similar profiles are obtained with aWMLES. No log-layer mismatch is observed unlike hybrid RANS/LES methods, [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF][START_REF] Larsson | Large eddy simulation with modeled wall-stress: Recent progress and future directions[END_REF] which justifies the numerical setup employed in this study.

Velocity and temperature fluctuations at Re τ = 1020 and 5186 are shown in Figs. 678. In all cases, aWMLES and iWMLES predict almost identical turbulent fluctuations whose overall levels are either in good agreement with DNS when it is available or coherent otherwise. The near-wall peak on streamwise velocity is not captured due to the coarse grid used and seems to start to be captured when the grid is refined. The only major difference between iWMLES and aWMLES occurs in the temperature fluctuation profiles at the first offwall cell. As only DNS data at Re τ = 1020 are available, it cannot be said whether iWMLES's temperature fluctuation estimation at the first cell is more accurate than that of aWMLES. Relative errors on the bulk friction coefficient and the Nusselt number are shown in Table V. As observed with the mean profiles, accurate wall friction fluxes are predicted by both iWMLES and aWMLES. In adiabatic cases, a maximum error of 3.1% is obtained. In isothermal cases, errors on the friction coefficient increase with the Reynolds number with maximum errors of 9.1% for aWMLES and 8.5% for iWM-LES. For Re τ ≥ 2003, the reference wall friction is taken from incompressible adiabatic DNS. In WMLES, as a compressible code is considered and due to the imposed temperature gradient, a density gradient exists. Therefore, the assumption that the DNS friction coefficient is still valid may not hold when the Reynolds number increases. The largest errors on the Nusselt number occur at the lowest Reynolds number. This result may be due to the position of the first off-wall cell approaching the buffer layer. At higher Reynolds number cases, iWMLES and aWMLES are both in agreement with Kays's correlation.

To sum up, iWMLES is able to recover the standard logarithmic laws for both the velocity and temperature profiles. Velocity and temperature profiles (mean and fluctuations) are in good agreement with adiabatic and isothermal DNS data. It is interesting to note that the WMLES using either aWM-LES or iWMLES predicts nearly identical fluctuations even if they are based on two different approaches. From a numerical point of view, differences could be observed using high-order schemes. From a physical point of view, velocity and temperature fluctuations are expected to be much more dependent on the features of the numerical methods (convective flux scheme, subgrid-scale model, etc.) used than on the wallmodel. Indeed, this result is supported by recent studies of Cossu and Hwang. [START_REF] Cossu | Self-sustaining processes at all scales in wallbounded turbulent shear flows[END_REF] They show that outer layer large-scale structures of wall-bounded turbulent flows can self-sustain without near-wall smaller-scale, by extracting energy from the mean flow. In WMLES, near-wall grids are too coarse to capture near-wall streaks but only large turbulent structures. Therefore, turbulence fluctuations would depend mostly on the resolution of these large-scale structures. Unlike near-wall small-scale structures, they are not expected to depend heavily on wall fluxes. This could explain the relative independence of the velocity fluctuations with respect to the wall-model used.

D. Results for isothermal supersonic flows

More discriminant test cases, namely, supersonic channel flows, are now considered. Three different cases in terms of Reynolds and Mach numbers are computed as shown in Table IV. In all cases, a constant wall temperature is imposed in order to cool the flow and evacuate the injected energy through forcing terms as pointed out by Trettel and Larsson. [START_REF] Trettel | Mean velocity scaling for compressible wall turbulence with heat transfer[END_REF] Usual logarithmic laws of the wall are not expected to hold on such flows. Thanks to its correction terms and computation of a weak solution of the vertically integrated boundary layer equations, iWMLES is expected to be more accurate than aWMLES.

Mean velocity profiles in wall units are shown in Fig. 9. As expected, aWMLES imposes the first off-wall cell to follow Reichardt's law which is incorrect here. Due to compressibility effects, velocity profiles are shifted upwards. This effect has also been observed on boundary layer flows. [START_REF] Lagha | A numerical study of compressible turbulent boundary layers[END_REF] Therefore, the friction velocity is overestimated and the whole velocity profiles in wall units are below the DNS profiles. On the contrary, iWMLES predicts correctly the wall friction. In all cases, the first off-wall point is in agreement with the DNS data as well as the whole profile. Only the second and third off-wall cells are shifted under the DNS velocity profiles for the case G2 M1.7 663 (Re τ = 663 on the finer grid G2). This can be due to the position of the wall-model interface which is in the buffer layer.

Likewise, mean temperature profiles in wall units for the aWMLES are not in agreement with the DNS profiles as shown in Fig. 9. The wall heat flux is systematically underestimated as the first off-wall cell is above the DNS. Temperature profiles are known to be more affected than velocity profiles. [START_REF] Houra | Effects of adverse pressure gradient on heat transfer mechanism in thermal boundary layer[END_REF][START_REF] Huang | Law of the wall for turbulent flows in pressure gradients[END_REF] The iWMLES model predicts more accurate wall heat flux than aWMLES, and more coherent temperature profiles are obtained. As observed for the velocity in case G2 M1.7 663, the temperature profile does not follow the usual trend observed in other cases: the wall heat flux is slightly overestimated in iWMLES and the temperature profile is below the DNS reference.

Observations on mean profiles are confirmed by the measured errors on the friction Mach number M τ and the heat flux coefficient B q in Table VI. Indeed larger errors are obtained with aWMLES compared to iWMLES. The maximum error on the friction Mach number is of 6.7% with the former and only 2.7% with the latter. For the heat flux coefficient, larger errors are computed as it is underestimated by at most 21.3% with aWMLES and 15.2% with iWMLES. Here, the aWMLES model shows its limits as the standard logarithmic laws are no longer valid on these kinds of flows. In all cases, iWMLES predicts more accurately wall fluxes than aWMLES.

Velocity and temperature fluctuations for Re τ = 972 are shown, respectively, in Figs. 10 and11. Other cases are not displayed for conciseness, as equivalent results and conclusions are obtained. Despite predicting different mean profiles, velocity fluctuations are very similar between iWMLES and aWMLES even if the former takes into account more physics than the latter in order to compute wall fluxes. Indeed, wallmodels implemented here influence only the mean profiles

Case

Wall-model ∆(M τ ) in % ∆(B q ) in % G1 M1. data, except the first rows of cells above the wall where the fluctuations are overestimated. Compared to quasi-incompressible flow DNS, the near-wall peak is more spread out in supersonic flows. This can explain the overestimation of the fluctuations near walls. Indeed when the grid is refined, the near-wall peak is better captured and better agreement with the DNS data is obtained.

Finally, in contrast to the velocity fluctuations, the temperature fluctuations are not similar between aWMLES and iWMLES. Indeed, temperature fluctuations with iWMLES are in correct agreement with the DNS profiles on the coarser grid, while aWMLES strongly overestimates it. Moreover, when the grid is refined, better agreement is obtained as the first off-wall cell is located at the near-wall peak position and the iWMLES approach predicts temperature fluctuations close to the DNS reference.

To sum up, iWMLES has a wider domain of validity and is able to compute more accurate wall fluxes than aWMLES, especially the wall friction. Better mean profiles are obtained with iWMLES. Furthermore, as shown in Subsection III E, iWMLES cost is only slightly more expensive compared to aWMLES.

V. CONCLUDING REMARKS

In this study, an attempt is made to develop a valid wallmodel for a typical turbomachinery flow range, i.e., flows with Mach numbers up to about 1.5, Reynolds numbers between 10 6 and 10 7 , and temperature gradients in the order of ten to hundred Kelvin. For this purpose, a compressible extension of the integral wall-model approach originally introduced for incompressible flows has been proposed, namely, the iWMLES. To this end, both the mean velocity and the mean temperature profiles have been parameterized using logarithmic laws of the wall with additional correction terms. These profiles are coupled to integral relations for compressible thin boundary layers in order to determine unknown parameters. Therefore, iWM-LES takes into account more physics than standard equilibrium wall-models. Besides, the computational cost is by a factor of about 100 compared to wall-models based on the numerical resolution of thin boundary layer equations since only a local simple scalar system is solved. An analytical wall-model, namely, the aWMLES, based on the standard logarithmic laws of the wall, is also implemented and used as a basis of comparison.

iWMLES is validated on both isothermal and adiabatic wall cases and is observed to accurately predict both subsonic and supersonic flows. Indeed, in subsonic flows, mean velocity and temperature profiles are in agreement with the DNS reference data, as well as the turbulent fluctuation profiles. iWMLES is therefore able to recover the standard logarithmic laws of the wall. In supersonic flows, as expected, aWMLES fails to predict the mean flow. On the contrary, iWMLES, by estimating more accurately the wall fluxes, is able to obtain results in agreement with the DNS data. Interestingly, both wall-models considered in this study predict nearly identical velocity fluctuations. This result indicates that the turbulent fluctuations are weakly dependent on the wall-model.

To go further, and before applying iWMLES to turbomachinery applications, it could be interesting to apply iWMLES in flows with (i) strong non-equilibrium effects, (ii) strong compressibility effects, and (iii) complex wall boundaries. These effects could be studied by simulating an oblique shock/boundary-layer interaction, hypersonic flows, or wavy walls. [START_REF] Nam | A ghost-cell immersed boundary method for large-eddy simulations of compressible turbulent flows[END_REF][START_REF] Sun | Direct numerical simulation of a fully developed compressible wall turbulence over a wavy wall[END_REF] For the latter case, taking into account curvature effects could be considered. However, this case has, for the moment, not been simulated at a Reynolds number high enough for WMLES. Furthermore, DNS of spatially evolving supersonic boundary layers with isothermal walls [START_REF] Shahab | Statistical characteristics of an isothermal, supersonic developing boundary layer flow from DNS data[END_REF][START_REF] Chu | Effect of wall temperature on hypersonic turbulent boundary layer[END_REF][START_REF] Hadjadj | Effect of wall temperature in supersonic turbulent boundary layers: A numerical study[END_REF][START_REF] Shadloo | Statistical behavior of supersonic turbulent boundary layers with heat transfer at M ∞ = 2[END_REF][START_REF] Shadloo | Laminar-turbulent transition in supersonic boundary layers with surface heat transfer: A numerical study[END_REF] has been recently performed. Friction Reynolds numbers reached are increasing, and these cases could later be used as reference cases in order to evaluate wall-models for heat transfer prediction. Nonetheless, iWMLES, by taking into account more physics, is able to extend the domain of validity observed with a standard analytical wall-model with similar computational cost.

θ

  and N ∼ Re1.

17

 17

  θ

  FIG. 1. Scheme of the projection of the first off-wall point velocity onto the local boundary layer frame.

  (e) The wall fluxes at the current time step u n τ and T n τ are determined in such a way that the left-hand sides of Eqs. (25a) and (25b) equal their right-hand sides using the Newton-Raphson algorithm. u n-1 τ and T n-1 τ are used as the initial guess. For each evaluation, the scalar parameters A n , A n T , δ c , and δ ν are computed and the corresponding velocity and temperature profiles are used to approximate integral terms L n f , thanks to the Gauss-Legendre quadrature. (f) Once the Newton-Raphson algorithm has converged, the wall friction modulus τ n w and wall heat flux φ w n are fed back to the LES solver. The wall friction vector components in the flow solver frame are finally determined by switching from the boundary layer to the flow solver frame as explained in Subsection II B.

2h 0 ρ

 0 dy/(2h) is the bulk density. For the compressible DNS under consideration, Re τ is equal to 1015, 663, or 972. At Re τ = 1015, the Mach number M b = u b / √ γrT w based on the bulk velocity u b = ∫ 2h 0 ρu dy/(2h ρ b ) is equal to 1.5. In the other two cases, M b = 1.7. For all cases, only non-adiabatic boundary conditions exist.

FIG. 5 .

 5 FIG.5. Mean velocity (a) and temperature (b) profiles in wall units for isothermal quasi-incompressible cases. Profiles are shifted by multiples of 10. From bottom to top: G1 M0.2 R1020, G2 M0.2 R1020, G1 M0.2 R2003, G2 M0.2 R2003, G1 M0.2 R4179, G2 M0.2 R4179, G1 M0.2 R5186, G2 M0.2 R5186.

FIG. 7 .FIG. 8 .

 78 FIG. 7. Velocity fluctuations for isothermal cases. For each case, from bottom to top: + , + , and u + . (a) Case G1 M0.2 R1020 isot. (b) Case G2 M0.2 R1020 isot. (c) Case G1 M0.2 R5186 isot. (d) Case G2 M0.2 R5186 isot.

FIG. 9 .

 9 FIG. 9. Mean velocity (a) and temperature (b) profiles in wall units for isothermal supersonic cases. Profiles are shifted by multiples of 10. From bottom to top: G1 M1.5 R1015, G2 M1.5 R1015, G1 M1.7 R663, G2 M1.7 R663, G1 M1.7 R972, G2 M1.7 R972.

FIG. 10 .

 10 FIG. 10. Velocity fluctuations for isothermal supersonic cases. For each case, from bottom to top: + , + , and u + . (a) Case G1 M1.7 R972 isot. (b) Case G2 M1.7 R972 isot.

  only bidimensional equations are considered here due to the projection step explained in Subsection II B. The two-dimensional turbulent thin boundary layer equations in Cartesian coordinates for a perfect gas are

	with the shear stress τ = (µ + µ t ) ∂u ∂y and the heat flux
	φ = -(λ + λ t ) ∂T ∂y , and λ =	c p µ
				∂ ρ ∂t	+	∂ ρu ∂x	+	∂ ρv ∂y	= 0,	(9a)
			∂ ρu ∂t	+	∂ ρuu ∂x	+	∂ ρvu ∂y	= -	∂p ∂x	+	∂τ ∂y	,	(9b)
							∂p ∂y	+		∂ ρv 2 ∂y	= 0,	(9c)
		∂ ρc p T ∂t	+	∂ ρuc p T ∂x	+	∂ ρvc p T ∂y	=	∂p ∂t	+ u	∂p ∂x	+ τ	∂u ∂y	-	∂φ ∂y	,
															(9d)
												p = ρrT ,	(9e)

TABLE I .

 I Plane channel flow DNS considered.

	DNS		(N x , N y , N z )	Boundary conditions	M b	Re τ	Re b	M τ	T b /T	B q / Nu
	Abe, Kawamura, and Matsuo 70	(2048, 448, 1536)	Uniform heat-flux	. . .	1020	20 721	. . .	. . .	. . . /83.0
	Hoyas and Jiménez 71		(6144, 633, 4608)	Adiabatic	. . .	2003	43 590	. . .	. . .	. . .
	Lozano-Durán and Jiménez 72	(-,1081,-)	Adiabatic	. . .	4179	98 302	. . .	. . .	. . .
	Lee and Moser 73		(10240, 1536, 7680)	Adiabatic	. . .	5186	124 862	. . .	. . .	. . .
	Modesti and Pirozzoli 75	(2048, 512, 1024)	Isothermal	1.5	1015	17 000	0.065	1.35	0.038/123.5
	Trettel and Larsson 76		(800, 246, 400)	Isothermal	1.7	663	10 000	0.077	1.45	0.053/68.9
	Trettel and Larsson 76		(896, 384, 480)	Isothermal	1.7	972	15 500	0.073	1.45	0.050/112.7
	is considered adiabatic. For the WMLES, a small temperature
	gradient T b /T = 1.1 will also be applied in order to assess the
	accuracy of the wall heat flux computation. The adiabatic DNS
	friction coefficients are still considered as reference. The ref-
	erence Nusselt number Nu = 2hφ /(λ(T -T b )) is computed
	from Kays's correlation, 74
	Nu =	0.023Re 0.8 f Pr 0.88 + 2.03(Pr 2/3 -0.78)Re -0.1

f ,

TABLE II .

 II Adiabatic quasi-incompressible plane channel flow cases in WMLES.

	Case name	(N x , N y , N z )	M b	Re τ	Re b	∆ + x	y + 1 (at cell)	∆ + z
	G1 M0.2 R2003 adia	(25, 21, 21)	0.2	2003	43 590	524	100	315
	G2 M0.2 R2003 adia	(49, 41, 41)	0.2	2003	43 590	262	50	157
	G1 M0.2 R4179 adia	(25, 21, 21)	0.2	4179	98 302	1094	209	656
	G2 M0.2 R4179 adia	(49, 41, 41)	0.2	4179	98 302	547	104	328
	G1 M0.2 R5186 adia	(25, 21, 21)	0.2	5186	124 862	1358	259	815
	G2 M0.2 R5186 adia	(49, 41, 41)	0.2	5186	124 862	679	130	407

TABLE III .

 III Isothermal quasi-incompressible plane channel flow cases in WMLES.

	Case name	(N x , N y , N z )	M b	Re τ	Re b	T /T b	∆ + x	y + 1 (at cell)	∆ + z
	G1 M0.2 R1020	(25, 21, 21)	0.2	1020	20 721	1.1	267	51	160
	G2 M0.2 R1020	(49, 41, 41)	0.2	1020	20 721	1.1	134	26	80
	G1 M0.2 R2003	(25, 21, 21)	0.2	2003	43 590	1.1	524	100	315
	G2 M0.2 R2003	(49, 41, 41)	0.2	2003	43 590	1.1	262	50	157
	G1 M0.2 R4179	(25, 21, 21)	0.2	4179	98 302	1.1	1094	209	656
	G2 M0.2 R4179	(49, 41, 41)	0.2	4179	98 302	1.1	547	104	328
	G1 M0.2 R5186	(25, 21, 21)	0.2	5186	124 862	1.1	1358	259	815
	G2 M0.2 R5186	(49, 41, 41)	0.2	5186	124 862	1.1	679	130	407

TABLE IV .

 IV Isothermal supersonic plane channel flow cases in WMLES.

	Case name	(N x , N y , N z )	M b	Re τ	Re b	T b /T	∆ + x	y + 1 (at cell)	∆ + z
	G1 M1.5 R1015	(25, 21, 21)	1.5	1015	17 000	1.35	266	51	159
	G2 M1.5 R1015	(49, 41, 41)	1.5	1015	17 000	1.35	133	25	80
	G1 M1.7 R663	(25, 21, 21)	1.7	663	10 000	1.45	174	33	104
	G2 M1.7 R663	(49, 41, 41)	1.7	663	10 000	1.45	87	17	52
	G1 M1.7 R972	(25, 21, 21)	1.7	972	15 500	1.45	254	49	153
	G2 M1.7 R972	(49, 41, 41)	1.7	972	15 500	1.45	127	24	76

TABLE V .

 V Relative errors in % on the bulk friction coefficient C f b = τw

					1 2 ρ b u 2 b
	and the Nusselt number Nu = channel flow cases in WMLES.	2hφw λ(Tw -T b ) for quasi-incompressible plane
	Case	Wall-model	∆(C f b ) in %	∆(Nu) in %
	G1 M0.2 R2003 adia	aWMLES	0.3	. . .
	G2 M0.2 R2003 adia	aWMLES	0.7	. . .
	G1 M0.2 R4179 adia	aWMLES	1.3	. . .
	G2 M0.2 R4179 adia	aWMLES	2.5	. . .
	G1 M0.2 R5186 adia	aWMLES	2.3	. . .
	G2 M0.2 R5186 adia	aWMLES	3.1	. . .
	G1 M0.2 R2003 adia	iWMLES	0.3	. . .
	G2 M0.2 R2003 adia	iWMLES	0.1	. . .
	G1 M0.2 R4179 adia	iWMLES	2.0	. . .
	G2 M0.2 R4179 adia	iWMLES	1.4	. . .
	G1 M0.2 R5186 adia	iWMLES	2.5	. . .
	G2 M0.2 R5186 adia	iWMLES	2.3	. . .
	G1 M0.2 R1020	aWMLES	3.8	13.7
	G2 M0.2 R1020	aWMLES	5.1	14.6
	G1 M0.2 R2003	aWMLES	5.3	1.9 a
	G2 M0.2 R2003	aWMLES	5.3	1.3 a
	G1 M0.2 R4179	aWMLES	7.1	0.04 a
	G2 M0.2 R4179	aWMLES	8.0	1.0 a
	M0.2 R5186	aWMLES	8.2	0.7 a
	G2 M0.2 R5186	aWMLES	9.1	1.8 a
	G1 M0.2 R1020	iWMLES	1.2	7.4
	G2 M0.2 R1020	iWMLES	7.1	14.0
	G1 M0.2 R2003	iWMLES	4.7	4.5 a
	G2 M0.2 R2003	iWMLES	6.5	2.2 a
	G1 M0.2 R4179	iWMLES	7.7	1.8 a
	G2 M0.2 R4179	iWMLES	6.8	2.2 a
	G1 M0.2 R5186	iWMLES	7.9	1.5 a
	G2 M0.2 R5186	iWMLES	8.5	0.7 a

a Reference Nusselt number is computed from Kays's correlation.

ACKNOWLEDGMENTS

The authors acknowledge A. Trettel and J. Larsson for kindly providing the DNS data of their supersonic channel flows. This work was supported by the Association Nationale de la Recherche et de la Technologie (ANRT; Grant No. 2015/0193).

APPENDIX A: JACOBIAN MATRIX OF THE iWMLES EQUATIONS

As explained in Subsection III E, the Newton-Raphson algorithm is used to solve the system of Eqs. (25). In this appendix, the corresponding Jacobian matrix is detailed in the case of an isothermal boundary condition. If the wall is considered to be adiabatic, the Jacobian matrix is reduced to a one-dimensional function.

Let f and g be two functions defined by

-∆t(-M n-1

x

Variables at time step n -1 are known and can be considered as constant scalars of the equation system as well as the pressure p n which is an input from the LES to the wall-model equations. Friction velocity u n τ and temperature T n τ at time step n are the unknowns. The system of Eqs. ( 25) can be simply written as

and the corresponding Jacobian matrix J is by definition

J components are given by

Thus, J is completely defined by the derivatives of integral terms L n ρ and L n ρu with respect to u n τ and T n τ . As the fluid is supposed to be a perfect gas, density and temperature profiles are linked with pressure ρ n (y) = p n 1 /(rT n (y)) and

dy, (A6c)

Temperature T n and velocity profile u n between the wall assumed to be located at y = 0 and the first off-wall point y 1 are given by Eqs. ( 17) and ( 18), so )) , (A7a)

and

with κ, κ t , β, and Γ defined in Eq. ( 8). So ∂Γ ∂u n τ = 10 -2 (Pry +n ) [START_REF] Sagaut | Large eddy simulation for aerodynamics: Status and perspectives[END_REF] 4 + 15Pr 3 y +n u n τ (1 + 5Pr 3 y +n ) 2

. (A9)

The derivatives of A n and A n T are obtained from Eq. ( 24),

where Γ 1 = Γ(y 1 ). Note that ∂A n

∂T n τ = 0. Therefore, using Eqs. (A7), (A8), and (A10), the integrated terms in Eqs. (A6) and so the Jacobian matrix components (A5) can be expressed analytically. From an initial guess (u n τ 0 , T n τ 0 ), the Newton-Raphson algorithm consists in computing a solution to Eq. (A3) by iteration

As explained in Subsection III E, the wall friction velocity and the wall friction temperature at previous time step are taken as initial guess. At each iteration of the algorithm, a Gauss-Legendre quadrature is used to compute the Jacobian matrix components. The algorithm is stopped once f (u n τ ; T n τ ) and g(u n τ ; T n τ ) are both below a given threshold. In this study, the threshold is set at 10 -13 .

In the case of an adiabatic boundary condition, T τ is zero. Therefore, ∂f ∂u n τ is the only component of J and

where ∂ū n ∂u n τ is given by Eq. (A7a).

APPENDIX B: WALL-NORMAL DERIVATIVES OF THE PARAMETERIZED PROFILES

In iWMLES equations (25), flux values at the first offwall point τ