Sébastien Deck
email: sebastien.deck@onera.fr

Nicolas Renard

Towards an enhanced protection of attached boundary layers in hybrid RANS/LES methods

Keywords: Hybrid RANS/LES methods, Detached Eddy Simulation, ZDES, boundary layer shielding, Grid Induced Separation, Modelled Stress Depletion, Grey area mitigation

A robust strategy for the RANS shielding of attached boundary layers in a hybrid RANS/LES context is presented, addressing two major issues. Firstly, the attached boundary layers must be detected to ensure their RANS treatment, but the original shielding function used by automatic methods such as DDES (2006) or ZDES mode 2 (2012) fails for fine meshes and/or with adverse pressure gradients, motivating the present study. Secondly, even with a stronger shielding, there should not be an excessive delay in the formation of instabilities and resolved LES content in free shear layers. Both objectives cannot be simultaneously reached by only tuning a constant in the original shielding function. The proposed new ZDES mode 2 consequently involves three main ingredients, namely a second shielding function detecting the outer part of the wake layer including with strong adverse pressure gradients, an inhibition function which switches off the second shielding when flow separation is detected, and a significant enhancement of the destruction of eddy viscosity in grey areas. The calibration of the method relies on a set of boundary layers at various Reynolds numbers and pressure gradients conditions and is confirmed by a priori tests in three-dimensional flows around curved geometries. The resulting case-independent model, the new ZDES mode 2, is assessed on four test cases, namely a flat-plate boundary layer, a mixing layer, a backward facing step and transonic buffet over a supercritical airfoil. Overall, it is shown that the new ZDES mode 2 is relevant with respect to four objectives: protection of attached boundary layers for any grid cell size (including infinite mesh refinement) and pressure gradient, RANS shielding at least as broad as the original shielding in any situation, minimum delay in the formation of resolved LES content, and full compatibility of the resulting subgrid scale model in the LES branch with the other modes of ZDES (modes 1 and 3) ensuring a continuous treatment of resolved turbulence across zones treated with different ZDES modes. The new robust ZDES mode 2 consequently is a case-independent answer to the demand for a general automatic and robust RANS/LES treatment of attached and massively separated flows.

Introduction

Computational Fluid Dynamics (CFD) is playing an expanding role in aerospace design, and enhancing its capabilities especially for unsteady flow simulations is considered as a major objective in the field [START_REF] Slotnick | CFD Vision 2030 Study: a path to revolutionary computational aerosciences[END_REF][START_REF] Deck | High-fidelity simulations of unsteady civil aircraft aerodynamics: stakes and perspectives. application of zonal detached eddy simulation[END_REF]. Among turbulence-resolving methods, the hybrid RANS/LES paradigm is now widely acknowledged as one of the main strategies to drastically reduce computational cost in a wide range of complex industrial applications, compared with Direct Numerical Simulation (DNS) or standard Wall Resolved Large Eddy Simulation (WRLES) techniques. Indeed, the turbulent fluctuations may be resolved in LES only where needed, whereas the rest of the flow is modelled in RANS at a significantly lower cost. Most of the time, resolving turbulence is not needed in attached boundary layers. In this case, free shear layers and massively separated flows may be treated in LES whereas attached flows are described in RANS, sparing the otherwise very high cost of resolving wall-bounded turbulence in attached boundary layers. Even in the case when the user wishes to resolve turbulence in the outer layer of an attached boundary layer, hybrid RANS/LES methods still provide a significant saving compared with WRLES since a near-wall RANS treatment results in a Wall-Modeled LES behaviour [START_REF] Deck | Zonal Detached Eddy Simulation (ZDES) of a spatially developing flat[END_REF], but this less frequent case is out of the scope of the present paper.

Among the hybrid RANS/LES methods, the approach that has probably drawn most attention in the recent time frame is the Detached Eddy Simulation (DES97) which was proposed by Spalart et al. [START_REF] Spalart | Comments on the feasibility of LES for wings and on a hybrid RANS/LES approach[END_REF]. The idea is to simulate the attached boundary layer in RANS mode whereas the separated flow should be simulated in LES mode. Ideally, the hybrid model works in RANS mode inside the attached boundary layer and in LES mode outside. Starting from the one-equation Spalart-Allmaras turbulence model [START_REF] Spalart | A one equation turbulence model for aerodynamic flows[END_REF], the wall distance d w used as a turbulent length scale for near-wall destruction of the eddy viscosity in its transport equation is replaced in LES regions by the estimate C DES ∆ of the cell size, so that the model provides typical subgrid-scale viscosity levels in a Smagorinsky-like behaviour. The switch between the RANS and LES modes is simply determined by comparing between d w and C DES ∆ as the characteristic length scale d replaces d w in the transport equation for the eddy viscosity:

d = min(d w , C DES ∆) (1)
In the natural or "safe" use of DES illustrated in figure 1, the mesh in attached boundary layers is not too fine, so that C DES ∆ is greater than the local boundary layer thickness δ. Consequently, d is equal to d w in the whole attached boundary layer, leading to the desired RANS treatment, whereas farther from the wall, in free shear layers, d = C DES ∆ provides a LES behaviour.

The "DES-type" methods where the attached boundary layers are modelled in RANS mode can be considered as "weak" RANS-LES coupling methods since there is no mechanism to transfer the modelled turbulence energy into resolved turbulence energy. Consequently, these methods introduce a "grey-area" in which the solution is neither pure RANS nor pure LES since the switch of the model from RANS to LES does not imply an instantaneous change in the resolution level.

Practically, the eddy viscosity remains continuous across the RANS/LES interface but the rapid decrease of the level of RANS eddy viscosity enables the development of strong instabilities populating massively separated flows (see [START_REF] Sagaut | Multiscale and multiresolution approaches in turbulence-LES, DES and hybrid RANS/LES methods: applications and guidelines[END_REF][START_REF] Fröhlich | Hybrid RANS/LES methods for the simulation of turbulent flows[END_REF] for further discussion). The smaller the grey-area, the better the model since a rapid RANS to LES switch is desirable.

In practice, two weaknesses in the use of hybrid methods for applied aerodynamics are classically identified, as illustrated by fig. 1.

• Issue # (I): the first weakness occurs when a "grey area", i.e. an area where the model switches from RANS to LES as d [START_REF] Ashton | Recalibrating Detached-Eddy Simulation to eliminate modelled-stress depletion[END_REF] switches from d w to C DES ∆, is located inside the attached boundary layer. This happens typically when the mesh is very fine (C DES ∆ < δ). Properly resolving turbulence inside the boundary layer would however require an even finer mesh suitable for WMLES and the upstream injection of resolved turbulence into the boundary layer. If this is not the case, the velocity fluctuations are expected to be insufficiently developed to compensate for the loss of modelled turbulent stresses (also referred to as "model-stress depletion" or MSD in the literature). This can lead to non-physical outcomes, like an underestimation of the skin friction. At worst, the separation line moves too far forward leading to a Grid Induced Separation (GIS) [START_REF] Menter | A Scale-Adaptive Simulation Model for Turbulent Flow Predictions[END_REF]. In order to get rid of this latter drawback, Spalart et al. [START_REF] Spalart | A new version of Detached-Eddy Simulation, resistant to ambiguous grid densities[END_REF] proposed a modification of the model

so that the model is able to refuse the "LES" mode if the shielding function indicates the presence of a boundary layer:

f d = 1 -tanh C 1 ν t + ν √ u i,j u i,j κ 2 d 2 w C2 C 1 = 8 and C 2 = 3 (3)
where u i,j = ∂u i /∂x j is the gradient of the instantaneous velocity. √ u i,j u i,j may also be written more explicitly as i,j ∂ui ∂xj

2 1 2
.

Thanks to this modification, the model is supposed to stay in RANS mode in attached boundary layers and switch into LES mode only outside (where f d = 1).

An increasing number of studies ( [START_REF] Probst | A comparison of detached-eddy simulation and reynolds-stress modelling applied to the flow over a backward-facing step and an airfoil at stall[END_REF][START_REF] Ashton | Flow Dynamics Past a 30P30N Three-Element Airfoil Using Improved Delayed Detached Eddy Simulation[END_REF][START_REF] Jain | Assessment of shielding parameters in conventional DDES method under the presence of alternative turbulence length scales[END_REF] among others) have noted that the shielding provided by the f d function becomes insufficient when the grid spacing in the streamwise direction, say ∆x, becomes smaller than 0.3δ where δ is the local boundary layer thickness. An illustration is provided in fig. 2 where the classical DDES method of 2006 is shown to fail to protect the whole boundary layer profile on a very fine mesh (∆x/δ close to 0.1 in the finest region where the failure occurs, wheras the higher values of ∆x/δ in the upstream part of the domain enable a correct shielding), leading to a severe underprediction of mean skin friction. This is something that is on the increase with the use of even finer grids permitted by the current supercomputer facilities. Note that this situation occurs systematically in the frame of Cartesian/octree grid refinements classically encountered with Latice Boltzmann Model (LBM) solvers where very isotropic grids are generated. This can result in non-intuitive degradation of the results with mesh refinement. The issue is most problematic in adverse pressure gradient flows with the unphysically early flow separation that might follow. It should be reminded that in addition to grid refinement, a similar issue may affect the shielding of attached boundary layers because of the interaction with one of the transition functions of the Spalart-Allmaras model, as analysed by Vatsa et al. [START_REF] Vatsa | Grid sensitivity of SA-based Delayed-Detached-Eddy-Simulation model for blunt-body flows[END_REF]. Several temporary solutions have been proposed to keep the whole attached boundary layer in RANS mode. Firstly, it should be pointed out that the issue is less severe in a zonal framework such as the Zonal Detached Eddy Simulation technique [START_REF] Deck | Recent improvements of the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF], where the user has the possibility to define by hand a RANS zone where no undesirable LES switch can occur, using mode 0 of ZDES whereas the LES region is set to mode 1. In this framework, the authors strongly recommend to avoid interactions between RANS-modelled laminar-turbulent transition and the hybrid RANS/LES shielding. More precisely, the safest manner to prevent the troubles described by Vatsa et al. [START_REF] Vatsa | Grid sensitivity of SA-based Delayed-Detached-Eddy-Simulation model for blunt-body flows[END_REF] can be to devise RANS zones (in ZDES, the so-called "mode 0") which contain the laminar-turbulent transition locations. However, regarding the mesh refinement issue in hybrid RANS/LES zones, even ZDES is confronted with the insufficient shielding by the f d function because this function is used by the automatic mode 2 of ZDES, which aims to select the RANS and LES regions automatically, without any user input.

The most obvious temporary solution is to modify the constants involved in the shielding function definition (3), especially C 1 , in order to increase its protective capability. This has been investigated by several groups in the hybrid RANS/LES community. Probst et al. [START_REF] Probst | A comparison of detached-eddy simulation and reynolds-stress modelling applied to the flow over a backward-facing step and an airfoil at stall[END_REF] recommended an increased value C 1 = 16. In his recalibration study, Ashton [START_REF] Ashton | Recalibrating Detached-Eddy Simulation to eliminate modelled-stress depletion[END_REF] concluded that C 1 = 12 is a good compromise to restore the shielding functionality of the f d function without falling into an overly conservative shielding which would result in an excessive delay in the formation of instabilities in mixing layers in grey areas (see issue (II) below). Similarly, Jain et al. [START_REF] Jain | Assessment of shielding parameters in conventional DDES method under the presence of alternative turbulence length scales[END_REF] retain C 1 = 14 and emphasize that in their test case, over-shielding is observed for C 1 = 16. As for Siggeirsson et al. [START_REF] Siggeirsson | Sensitivity study of the SA-DDES shielding function[END_REF], they opt out for C 1 = 16, following the suggestion of Probst et al. [START_REF] Probst | A comparison of detached-eddy simulation and reynolds-stress modelling applied to the flow over a backward-facing step and an airfoil at stall[END_REF]. Note that though not reported in the original paper [START_REF] Spalart | A new version of Detached-Eddy Simulation, resistant to ambiguous grid densities[END_REF], the values of C 1 and C 2 have been assessed

by the authors of DDES on tests in the flat plate boundary layer. Spalart et al. indicate that "a value (of C 1 ) larger than 8 would delay LES in even larger regions, which would be safer in the sense of MSD, but is undesirable overall". In addition, the proposal of increasing the value of C 1 is heavily case-dependent and does not guarantee the robustness of the shielding for severe grid refinements as classically encountered in complex mesh geometries. As already mentioned, a too large value C 1 may delay too much the switch into LES mode which concerns the second important issue of any "DES-type" calculation (see issue (II) below). This clearly advocates for a new rather than recalibrated shielding function.

Probst et al. [START_REF] Probst | A comparison of detached-eddy simulation and reynolds-stress modelling applied to the flow over a backward-facing step and an airfoil at stall[END_REF] state that "a more general approach directly including the effects of pressure gradients without affecting equilibrium boundary layers is certainly desirable". After confirming the underprediction of the boundary layer thickness growth by f d (even recalibrated) near separation, they introduce the Algebraic Delayed Detached Eddy Simulation (ADDES) strategy [START_REF] Probst | Detached-Eddy Simulation of Aerodynamic Flows Using a Reynolds-Stress Background Model and Algebraic RANS/LES Sensors[END_REF][START_REF] Knopp | An algebraic sensor for the rans-les switch in delayed detached-eddy simulation[END_REF], which resorts to several sensors in order to identify the attached boundary layers.

One sensor relies on the compressible Bernoulli equation and is limited to external flows. The two other sensors resort to the d w Ω function, more precisely the extrema of it and of its first order derivative, which are somewhat sensitive quantities.

The existence of several sensors not always in agreement with each other is a consequence of the difficulty of devising a universal shielding valid for both internal and external flows including non-canonical situations.

One proposal of a new shielding function aiming at universality is embedded in the Shielded Detached Eddy Simulation (SDES) and Stress-Blended Eddy Simulation (SBES) methods available in the ANSYS CFD software [START_REF] Menter | Stress-Blended Eddy Simulation (SBES) -A New Paradigm in Hybrid RANS-LES Modeling[END_REF]. Unfortunately, this function is proprietary and its expression has not been published or communicated to date, meaning that it is not available to the hybrid RANS/LES community. Other works devoted to the shielding of boundary layers include the proposal of He et al. [START_REF] He | A dynamic delayed detached-eddy simulation model for turbulent flows[END_REF] and Dong et al. [START_REF] Dong | An enhanced version of delayed detached-eddy simulation based on the v 2 -f model[END_REF] who suggest that the shielded boundary layers should not be fully treated in RANS when the grid resolution is fine enough and resolved turbulence is already present, which somewhat deviates from the intended use of DDES where the whole attached boundary layer should be treated in RANS. Another proposal, by

Islam & Thornber [START_REF] Islam | A high-order hybrid turbulence model with implicit large-eddy simulation[END_REF], informs a RANS/LES blending function by means of a transport equation for a dedicated variable whose behaviour mimics the creation of boundary layers near walls. The physics behind this variable and transport equation seems however limited, while solving an additional transport equation implies a significant additional cost. One last study that one may mention is the work by Weihing et al. [START_REF] Weihing | Development of alternative shielding functions for Detached Eddy Simulations[END_REF], where two shielding strategies have been considered: one relying on a Bernoulli equation (B-DES) and the other based on the vorticity (Vorticity Induced Algebraic DES:

VIA-DES). However, the study is limited to external aerodynamics, which is especially obvious for the Bernoulli-based strategy. Using a Bernoulli equation may also be problematic in unsteady flows. The dependency on the reference frame could also be a problem, for example when the geometry features moving parts. Concerning the proposed alternative based on vorticity, Weihing et al. report calibration difficulties. Regarding the detection of LES content, they suggest that the σ operator of Nicoud et al. [START_REF] Nicoud | Using singular values to build a subgrid-scale model for large eddy simulations[END_REF] might indicate the presence of resolved 3D turbulence, whose detection is not straightforward.

• Issue # (II): the second difficulty for "DES-type" methods is the possible delay in the formation of instabilities in mixing layers due to the advection of the upstream RANS eddy viscosity, as illustrated by region (II) in figure 1. This problem deserves attention particularly due to its relevance to the field of acoustics [START_REF] Bogey | Large-eddy simulation of the flow and acoustic fields of a Reynolds number 10 5 subsonic jet with tripped exit boundary layers[END_REF]. Pioneering remedies relying on the definition of subgrid length scales ∆ ω based on the local vorticity field have been proposed in [START_REF] Chauvet | Zonal-Detached-Eddy Simulation of a Controlled Propulsive Jet[END_REF][START_REF] Deck | Recent improvements of the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF] and have recently inspired some other groups [START_REF] Shur | An Enhanced Version of DES with Rapid Transition from RANS to LES in Separated Flows[END_REF][START_REF] Mockett | Two non-zonal approaches to accelerate RANS to LES transition of free shear layers in DES[END_REF][START_REF] Fuchs | Recent results with grey-area improved DDES for a wide range of flows[END_REF]. It should also be emphasized that a workpackage of the European Union-funded

Go4Hybrid project [START_REF] Peng | Go4Hybrid: Grey Area Mitigation for Hybrid RANS-LES Methods -Free Shear Layer[END_REF] has been devoted to this issue of grey area mitigation. A challenging task consists in designing a method with an enhanced protection of attached boundary layers while preserving a rapid switch into LES mode in the case of free-shear flows. Conversely, enhancing the rapidity of switch from RANS to LES can have detrimental effects on the shielding of attached boundary layers if this shielding is not strong enough, as illustrated by the initial difficulties encountered by Shur et al. [START_REF] Shur | An Enhanced Version of DES with Rapid Transition from RANS to LES in Separated Flows[END_REF] who had to restrain the activation of their instability unlocking function F KH in order to avoid its corrupting the shielding of the attached boundary layer which would otherwise lead to a significant underestimation of skin friction. The presentation of the enhanced ZDES mode 2 will start with section §2 including a summary of the ZDES formulation currently in use (ZDES mode 2 (2012)) followed by a proposed new formulation for ZDES mode 2 (EP) defined in order to meet the above objectives. Among several features, this formulation involves two new sensor functions. Both functions are generally applicable to any hybrid RANS/LES method and will be proposed, calibrated (with case-independent constants) in section §3 and detailed for their use in ZDES in Appendix A. The method ZDES mode 2 (EP) is finally exercised and compared with the reference DDES (2006) and ZDES mode 2 (2012) models in section §4 on four demanding test cases, namely a compressible turbulent flat plate boundary layer, a spatially developing mixing layer, a backward facing step and the transonic buffet over a supercritical airfoil.

Method formulation

Formulation of the Zonal Detached Eddy Simulation approach currently in use

The Zonal Detached Eddy Simulation (ZDES) technique, developed by Deck from 2002 onwards, was first published in [START_REF] Deck | Zonal-Detached Eddy Simulation of the flow around a high-lift configuration[END_REF] and the complete formulation that proposes an efficient solution to prevent delay in the formation of instabilities has been more recently published in Ref. [START_REF] Deck | Recent improvements of the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF]. The latter formulation is the formulation of ZDES currently in use and its mode 2 will be designated ZDES mode 2 (2012) in the following. This approach takes full advantage of its zonal nature, not only to allow the user to specify RANS and LES regions, but also to make possible the use of various formulations (the so-called modes) within the same calculation. Besides, ZDES also provides an "automatic" operating option (referred to as mode 2 in the following) for which the switch between RANS and LES regions is dynamically set by the model itself.

Thus, ZDES offers an attractive flexibility in the treatment of turbulent flows in technical applications [START_REF] Deck | High-fidelity simulations of unsteady civil aircraft aerodynamics: stakes and perspectives. application of zonal detached eddy simulation[END_REF].

The starting point of ZDES is the one-equation Spalart-Allmaras turbulence model [START_REF] Spalart | A one equation turbulence model for aerodynamic flows[END_REF] in compressible formulation [START_REF] Deck | Development and application of Spalart Allmaras one equation turbulence model to three-dimensional supersonic complex configurations[END_REF] for the transport of the pseudo eddy-viscosity ν defined such that:

µ t = ρνf v1 (4)
in which the wall distance d w , which is the RANS turbulent length scale for near-wall destruction of ν, is replaced by the characteristic length d:

∂ t (ρν) + ∂ xi (ρνu i ) = c b1 Sρν + 1 σ (∂ xi ((µ + ρν) ∂ xi ν) + c b2 (∂ xi ν) (∂ xi (ρν))) -ρc w1 f w νd 2 (5)
The hybrid length scale d and the functions f v1 , S and f w in ZDES are defined zonally in order to provide either a RANS behaviour (original definitions, especially in the so-called mode 0) or a LES behaviour (the turbulence model is used as a subgrid scale model in the LES branches of modes 1 to 3), depending on the mode of ZDES to which the considered zone of the flow is set.

Three specific hybrid length scale formulations, also called modes, are then optimized to be employed on three typical flow field topologies (as sketched in Table 1). Mode 1 concerns flows where the separation is triggered by a relatively abrupt variation in the geometry; mode 2 is retained when the location of separation is induced by a pressure gradient on a gently-curved surface and mode 3 for flows where the separation is strongly influenced by the dynamics of the incoming boundary layer. This latter mode is often referred to as Wall-Modelled Large Eddy Simulation (WMLES) mode (see Ref. [START_REF] Deck | Zonal Detached Eddy Simulation (ZDES) of a spatially developing flat[END_REF][START_REF] Laraufie | A dynamic forcing method for unsteady turbulent inflow conditions[END_REF][START_REF] Renard | Improvements in the formulation of Zonal Detached Eddy Simulation for Wall-Modeled Large-Eddy Simulation[END_REF]) and also has a Wall-Resolved Large Eddy Simulation (WRLES) capability [START_REF] Deck | Zonal Detached Eddy Simulation of a spatially developing flat plate turbulent boundary layer[END_REF][START_REF] Deck | Large scale contribution to mean wall shear stress in high Reynolds number flat plate boundary layers up to Re θ = 13650[END_REF]. The ZDES method aims to treat all classes of flow problems indicated in Table 1 in a single model, with the hybrid length scale defined in the same table.

The current paper is mainly concerned by mode 2 of ZDES which can be considered as the "automatic" mode of ZDES conversely to modes 1 and 3 where the decision load of the user is increased. When the shielding is not active, the model reduces to a behaviour equivalent to mode 1, which is the natural mode for LES in the framework of ZDES. It should be noted that ZDES mode 3 also reduces to a behaviour equivalent to mode 1 in the outer region of attached boundary layers, confirming the compatibility between all three modes of ZDES (see objective 4 in the introduction). Mode 1 of ZDES is very similar to the original DES97 method and provides only a very limited near-wall RANS protection based on the cell size:

dI = min(d w , C DES ∆I ) (6)
ZDES mode 1 differs from DES97 on several aspects. The cell size ∆I is either based on the volume ∆I = ∆ vol = (∆x∆y∆z) 1/3 or based on the vorticity ∆I = ∆ ω = √ S ω , as introduced by Chauvet et al. [START_REF] Chauvet | Zonal-Detached-Eddy Simulation of a Controlled Propulsive Jet[END_REF] and generalised by Deck [START_REF] Deck | Recent improvements of the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF], where S ω is the average cross section of the cell normal to the vorticity vector ω (see [START_REF] Deck | Recent improvements of the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF] for more details). Using the maximum cell size ∆ max = max(∆x, ∆y, ∆z) is not recommended. Moreover, when a LES region is detected, i.e. when d w > C DES ∆I , the near-wall functions are deactivated:

if d w > C DES ∆I then f v1,mode 1 = 1, f v2,mode 1 = 0 and f w,mode 1 = 1 ( 7
)
Zonal Detached Eddy Simulation (ZDES)

dZDES =              dw (

Examples of Applications

Base flow [START_REF] Weiss | Control of the antisymmetric mode (m=1) for high reynolds axisymmetric separating/reattaching flows[END_REF][START_REF] Simon | Numerical simulation of the compressible mixing layer past an axisymmetric trailing edge[END_REF][START_REF] Simon | Zonal Detached Eddy Simulation of projectiles in the Subsonic and Transsonic Regimes[END_REF][START_REF] Mochel | Zonal Immersed Boundary Conditions: Application to a high Reynolds number afterbody flow[END_REF], Airfoil separation [START_REF] Deck | Numerical investigation of the flow dynamics past a threeelement aerofoil[END_REF][START_REF] Deck | High-fidelity simulations of unsteady civil aircraft aerodynamics: stakes and perspectives. application of zonal detached eddy simulation[END_REF][START_REF] Richez | Zonal Detached Eddy Simulation of Separated Flow Around a Finite-Span Wing[END_REF][START_REF] Zhang | Zonal Detached-Eddy Simulation of Turbulent Unsteady Flow over Iced Airfoils[END_REF], Turbulent Boundary Layers jets [START_REF] Chauvet | Zonal-Detached-Eddy Simulation of a Controlled Propulsive Jet[END_REF][START_REF] Deck | High-fidelity simulations of unsteady civil aircraft aerodynamics: stakes and perspectives. application of zonal detached eddy simulation[END_REF] jets [START_REF] Gand | Investigation of turbulence development in incompressible jets with Zonal Detached Eddy Simulation (ZDES) and synthetic turbulent inflow[END_REF][START_REF] Verrière | Zonal Detached Eddy Simulation of a dua-stream jet[END_REF] (WRLES [START_REF] Deck | Large scale contribution to mean wall shear stress in high Reynolds number flat plate boundary layers up to Re θ = 13650[END_REF][START_REF] Renard | On the scale-dependent turbulent convection velocity in a spatially developing flat plate turbulent boundary layer at Reynolds number Re θ = 13000[END_REF][START_REF] Deck | Zonal Detached Eddy Simulation of a spatially developing flat plate turbulent boundary layer[END_REF],WMLES [START_REF] Deck | Zonal Detached Eddy Simulation (ZDES) of a spatially developing flat[END_REF][START_REF] Laraufie | A dynamic forcing method for unsteady turbulent inflow conditions[END_REF][START_REF] Renard | Improvements in the formulation of Zonal Detached Eddy Simulation for Wall-Modeled Large-Eddy Simulation[END_REF][START_REF] Deck | A rapid and low noise switch from RANS to WMLES on curvilinear grids with compressible flow solvers[END_REF] User Load high low ("automatic mode") high The original ZDES mode 2 (2012) [START_REF] Deck | Recent improvements of the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF] clearly borrows ideas from DDES [START_REF] Spalart | A new version of Detached-Eddy Simulation, resistant to ambiguous grid densities[END_REF] since dII reads as:

dII = d w -f P max 0, d w -C DES ∆II (8)
where f P is the shielding function. Enhancements of the shielding function will be proposed in §3. In ZDES mode 2 (2012) it is simply given by the following expression:

f P = f d (r d ) (9)
where

f d (r) = 1 -tanh (C 1 r) C2 C 1 = 8 and C 2 = 3 (10)
and

r d = ν + ν √ u i,j u i,j κ 2 d 2 w (11)
In practice, the formulas of ZDES mode 2 (2012) differ from those of DDES in the definition of :

• the subgrid length scale ∆II :

∆II = (1 -test ∆ )∆ max + test ∆ (∆ vol or ∆ ω ) (12)
where ∆ max , ∆ vol and ∆ ω have already been defined when discussing mode 1, and

test ∆ = 0 if f d (r d ) ≤ f d0 , 1 otherwise (13)
with f d0 = 0.8. This definition of ∆II may be considered as zonal (though transparent for the user) because of the threshold f d0 = 0.8, which had been defined thanks to flat plate boundary layer calculations. The proposal of a new subgrid length scale ∆II was not a minor adjustment of the detached eddy simulation (DES) formulation, because the modified length scales depend not only on the grid (∆x, ∆y, ∆z), but also on the velocity gradients (u i,j ) and pseudo eddy-viscosity field (ν), since ∆II ≡ ∆II (∆x, ∆y, ∆z, d w , u i,j , ν, ν). Using equation 12 in LES areas for mode 2 solved the problem of delay in the formation of instabilities, which can have dramatic effects on the pressure (and, thus, acoustic) field [START_REF] Deck | Recent improvements of the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF].

• the treatment of the near-wall functions in LES mode. The RANS near-wall functions are used in the shielded areas whereas these functions switch to a LES-adapted behaviour identical to mode 1 in the regions that are not shielded:

f v1 = (1 -test ∆ )f v1 (ν/ν) + test ∆ f v1,mode 1 (14)
f v2 = (1 -test ∆ )f v2 (ν/ν) + test ∆ f v2,mode 1 (15)
f w = (1 -test ∆ )f w (ν, d w , S) + test ∆ f w,mode 1 (16)
S = |ω| + ν κ 2 d 2 w f v2 ( 17
)
where |ω| is the vorticity magnitude.

The formulation of ZDES mode 2 (2012) reminded above satisfies objective 4 of the introduction (compatibility with the other modes of ZDES in the LES regions) and objective 3 (quick formation of instabilities in free shear layers). However, the shielding of attached boundary layers provided by eq. ( 9) fails in the case of severe streamwise grid refinement, say ∆x/δ ≤ 0.3. This can be explained by the fact that equation [START_REF] Davis | Further experiments on supersonic turbulent flow development in a square duct[END_REF] uses the pseudo eddy-viscosity field ν to protect ν itself.

In other words, once ν is corrupted, the corrupted ν enters eq. ( 9) to protect the even more corrupted ν and so on. In the long-time range, the solution may suffer from severe Modelled Stress Depletion. The problem is even worse in presence of adverse pressure gradients.

The next subsection proposes a reformulation of ZDES mode 2 suitable for the introduction of a better shielding complying with objective 1 (full shielding even for infinitely fine meshes with pressure gradients), in a spirit meeting objective 2 (the original shielding function is kept) and 4 (the formulation reduces to ZDES mode 1 in LES regions), with a solution to preserve objective 3 (minimum delay in the formation of instabilities in free shear layers in spite of the stronger shielding). The new sensor functions introduced will be defined in section §3.

Salient principles of the enhanced formulation of ZDES mode 2

In order to reach all four objectives mentioned in the introduction, three key elements are added to ZDES mode 2 to obtain ZDES mode 2 (EP) (Enhanced Protection):

1. Second shielding function: a new sensor G ν is introduced in order to protect the regions where the original sensor r d [START_REF] Davis | Further experiments on supersonic turbulent flow development in a square duct[END_REF] fails to detect attached boundary layers. The main objective of G ν is consequently to identify the wake layer of attached boundary layers, as will be shown in the next section where the sensor definition will be provided. The second shielding function f d (G ν ) is directly obtained from the sensor G ν by applying the f d function [START_REF] Dacles-Mariani | Numerical/Experimental Study of a Wingtip Vortex in the Near Field[END_REF] to this argument. A protection function satisfying both objective 1 (protection in all cases thanks to the combination of the sensors r d and G ν ) and objective 2 (the protection is never less than the one provided by the original sensor r d ) may be obtained by multiplying both shielding functions to obtain a protection function

f P = f d (r d )f d (G ν ).
This expression would imply that the boundary layer is protected as soon as one of the sensors detects it.

2. Inhibition function: the increased shielding provided by f d (G ν ) might result in significant delays in the formation of instabilities in free shear layers in some cases. Indeed, the early stages of some free shear layers could be mistakenly detected as the wake layer of an attached boundary layer. In order to attain objective 3, the second shielding function should consequently be inhibited in these cases. To that purpose, an inhibition function f R (G Ω ) based on a second new sensor G Ω is introduced to release the instabilities in the free shear layers that would otherwise be stabilised by the shielding. The sensor G Ω detects detached flows and will be detailed in the next section together with the function f R applied to this sensor. The inhibition function f R (G Ω ) should be 0 when detached flows are detected and 1 in attached flows. In order to achieve objective 2, the inhibition function f R (G Ω ) does not deactivate the original shielding function f d (r d ) as it switches only the second shielding function f d (G ν ) off. The expression of the resulting protection function is consequently the following:

f P = f d (r d ) 1 -1 -f d (G ν ) f R (G Ω ) (18)
The test (test ∆ ) for the RANS/LES switch of the subgrid length scale and of the near-wall functions is redefined as follows:

test ∆ = 0 if f d (r d ) ≤ f d0 or if (1 -(1 -f d (G ν ))f R (G Ω )) ≤ f d0 , 1 otherwise (19)
3. Enhanced destruction in grey areas: assuming that attached boundary layers are perfectly protected, i.e. f P = 0 in their whole profiles, then the other regions of the flow, which are defined by f P > 0, should be treated in LES, i.e. f P = 1. If there are grey areas where 0 < f P < 1, the pseudo eddy-viscosity should be destructed in these areas in order to reach as quickly as possible LES levels. Since the protection function f P depends on the pseudo eddy-viscosity field (see the second shielding function f d (G ν ) introduced in the next section), the destruction of ν will indeed make f P closer to 1. In grey areas, the destruction term of ν is significantly enhanced by setting f w = 100.

This corresponds to a destruction 100 times as high as in the typical case of f w = 1 (let us be reminded that this is the value of f w in a RANS logarithmic layer and in the true LES regions), which is similar to the effect of a reduction of the characteristic length d by a factor 10 (see eq. ( 5)). The purpose is to induce a very strong destruction of ν as long as f P is not close to 1, so that the areas where f P is neither 0 nor 1 should almost disappear, which is the reason why the exact value (f w = 100) does not matter as long as it is a high value. The application of this strong enhancement of destruction must be isolated from the pure RANS and LES subgrid scale model behaviours. This is achieved by bounding its application between two thresholds. The threshold that guarantees that the LES subgrid scale model is not corrupted is simply the previously introduced test ∆ : if test ∆ = 1, then the model is considered as close to pure LES and the destruction must no be enhanced. A second threshold is introduced to ensure that the pure RANS model is not corrupted, defined as another test (test w ). In pure RANS, test w = 0, whereas when the destruction is enhanced the condition test w = 1 is satisfied. Assuming that f d0 is greater than 0.5 (typically true since the recommended value is f d0 = 0.8), test w is defined similarly to test ∆ :

test w = 0 if f d (r d ) ≤ 1 -f d0 or if (1 -(1 -f d (G ν ))f R (G Ω )) ≤ 1 -f d0 , 1 otherwise (20)
The destruction term is enhanced only when test w = 1 and test ∆ = 0. The resulting f w function is finally given by:

f w = (1 -test w )f w (ν, d w , S) + test w ((1 -test ∆ ) • 100 + test ∆ f w,mode 1 ) (21)
where

S = |ω| + ν κ 2 d 2 w f v2 (22)
This enhanced destruction in grey areas contributes to objective 3 (reduce the delay in the formation of instabilities in some free shear layers) without altering the protection of attached boundary layers thanks to the threshold provided by test w .

Concerning objective 4 (compatibility with the other modes of ZDES), the continuity of treatment of resolved turbulence in pure LES regions is easily achieved since the proposed ZDES mode 2 (EP) reduces to the behaviour of ZDES mode 1 when the protection is not active (i.e. f P = 1).

The second shielding function f d (G ν ) and the inhibition function f R (G Ω ) will be defined and calibrated in the next section, after which the detailed formulation and implementation of ZDES mode 2 (EP) will finally be provided in Appendix A.

Proposals for the shielding and inhibition functions

As introduced in the previous section, the enhancement of ZDES mode 2 requires two new functions: a second shielding function defined so that attached boundary layers are fully shielded when used in conjunction with the original function, and an inhibition function aiming to release instabilities in free shear layers in some cases where the second shielding function may be detrimental to the formation of instabilities. Such functions should be motivated by the flow physics rather than model-dependent. For this reason, general functions are introduced in section §3.1, which will be calibrated in §3.2 and adapted to the particular case of the Spalart-Allmaras-based ZDES technique in §Appendix A.

Motivation and expression of the shielding and inhibition functions

The present discussion is independent of the underlying RANS model used in the hybrid RANS/LES context. The only assumption made is the eddy viscosity assumption: it is assumed here that an eddy viscosity ν t is available and represents the expected RANS levels of turbulence in the flow. If the attached boundary layers are properly shielded, ν t will indeed be available using the model in RANS mode and will be able to enter the shielding functions as an argument.

In summary, the arguments available for the shielding functions are the instantaneous fields of velocity, density, pressure, (RANS) eddy viscosity, and the geometry of the problem, which includes the distance to the nearest wall d w .

It should be noted here that only fully turbulent boundary layers are considered. Indeed, as mentioned in the introduction, it is strongly recommended that the laminar-turbulent transition regions be treated in RANS mode (e.g. mode 0 in the case of ZDES) rather than by the hybrid RANS/LES branch of the model, so that there is no interaction between the description of the laminar-turbulent transition and the shielding of fully turbulent attached boundary layers. This is easy to achieve in zonal methods such as ZDES.

The original r d sensor, which will be kept (see objective 2 in the introduction), is defined in the general (modelindependent) case as follows:

r d = ν t + ν √ u i,j u i,j κ 2 d 2 w (23)
to which the function f d is applied:

f d (r) = 1 -tanh (C 1 r) C2 C 1 = 8 and C 2 = 3 (24)
The sensor r d represents the squared ratio between the turbulent length scale provided by the RANS model ν t / √ u i,j u i,j

and the turbulent length scale in the logarithmic layer of attached boundary layers κd w , completed with a contribution of the molecular viscosity so that the viscous region below the logarithmic layer is also detected. In the logarithmic layer,

r d is close to 1.
The f d function is very close to 0 for r d values greater than approximately 0.2. This implies that part of the wake layer is also shielded. However, the r d sensor is too low (more precisely, below 0.2) in the outermost part of the wake layer, which is consequently not protected. As mentioned in the introduction, one simple fix may be to modify the constants C 1 and C 2 so that the shielding is active for values of r d even lower than 0.2. However, this would excessively delay the formation of instabilities in some free shear layers. The literature review has also shown that there is no consensus about the recommended values of C 1 and C 2 for an enhanced shielding as they seem to be case-dependent.

Furthermore, it should be reminded that the motivation of the r d sensor is based on the properties of the logarithmic layer, making it perfectly suited for the log layer shielding, but not really adapted to the protection of the wake layer, especially because the latter layer is strongly dependent on the pressure gradient. For these reasons and also because of objective 2 (additive construction of the new shielding), the original shielding function f d (r d ) is kept with its constants C 1 and C 2 unchanged. This means that the boundary layer profile is shielded from the wall to the bottom of the wake layer by the original shielding function, whereas the second shielding function has to detect the rest of the wake layer.

Shielding function for the outer part of the wake layer: options not retained

The detection of the outermost part of attached boundary layers is a very difficult task, all the more as there are some similarities between the wake layer and free shear layers while the latter ideally should not be shielded. Several options have been investigated with no success. Some of them have been already considered by other groups and are mentioned in the literature review in the introduction together with some of their shortcomings.

The first class of sensors considered was based on the energy losses induced by boundary layers. This is for instance the case of a sensor based on the Bernoulli equation. However, such sensors suffer from a lack of generality since the situation is much more complex for internal and unsteady flows than for the simple case of steady external aerodynamics, especially regarding the definition of an outer state and the possible energy conservation outside of the boundary layers.

Reference frame issues may also arise in the case of bodies with relative motions.

A second class of sensors was sought based on kinematic properties, i.e. using mainly the velocity field. The most common sensors are based on the vorticity, but their use in general situations seems extremely difficult if one needs to precisely define where the edge of the boundary layer is. Indeed, the vorticity near the edge is typically small, it has to be compared with a reference value not easily chosen, and it may be mixed with other phenomena as soon as cases more complex than a canonical isolated boundary layer are considered.

Another quantity that is very interesting is the wall-normal velocity, whose wall-normal gradient ∂u n /∂n is affected by the displacement effect caused by the thickening of the boundary layer. This quantity is already commonly used for other purposes in laminar-turbulent transition modelling. Nevertheless, severe difficulties arise for the present objectives. The first difficulty is the Reynolds number dependence of ∂u n /∂n, which makes the design of a Reynolds-number invariant sensor a complex task. Another issue is the inviscid effect of pressure gradients associated to a non-uniform outer potential flow which can also include wall-normal gradients of the wall-normal velocity that may be confused with the viscous and turbulent displacement effect occurring inside the boundary layer. One other major problem is that the wall-normal gradient of wall-normal velocity caused by the displacement effect in boundary layers is not large compared with the principal stretching caused by the shear. This should not be a major problem in the boundary layer, considering that the principal stretching direction of the rate of strain tensor is tilted by 45 o with respect to the shear direction in the case of a simple shear, but a dramatic confusion can occur in the near wake of a body if the principal stretching direction coincides with the nearest wall direction (this may be for instance the case at one fraction of boundary layer thickness downstream the trailing edge of a flat plate). One last difficulty is the somewhat high sensitivity of the wall-normal velocity to numerical imperfections. For these reasons, a kinematic sensor was not adopted.

Proposed shielding function for the outer part of the wake layer

As no sensor based on energetic or kinematic properties could be found, it is proposed here that the second shielding function rely also on the eddy viscosity field, just as the original function based on r d . Since the objective is to work in pair with the original function and to cover the regions not shielded originally, the fact that both functions use the same kind of arguments is seen as a positive feature, especially given the successful use of f d (r d ) in many cases. The possible drawback is that the second shielding function will share with the original function the negative property that in absence of eddy viscosity (ν t = 0 everywhere), there is no shielding of the attached boundary layer. In other words, if the eddy viscosity has been totally destroyed in the boundary layer, then it may be impossible to recover it. This implies that similarly to DDES (2006) or ZDES mode 2 (2012), the enhanced model still has to be initiated from a RANS field with the corresponding ν t levels, which in practice is not an issue. Besides, the fact that the shielding of ν t depends on ν t itself has been shown to be a problem in the original shielding function since the original shielding decreases if ν t starts to decay, resulting in a self-sustained decay process. Nevertheless, this will not be problematic for the proposed second shielding function because it depends on the gradient of ν t rather than on ν t itself, as will be explained in the following.

The rationale of the new sensor is very simple: in the outer part of the wake layer, the eddy viscosity ν t decreases from high near-wall values to very small values outside the boundary layer. This should also be true in adverse pressure gradient boundary layers, where the new sensor is especially needed given the important role played by the wake layer in APG boundary layers, since the original sensor fails in theses regions. The sensor is consequently sensitive to negative wallnormal gradients of the eddy viscosity, ∂ν t /∂n < 0 and will be noted G νt . In order to obtain a Reynolds-number-invariant behaviour, it has to be non-dimensionalised. This is achieved using the same quantities as for the non-dimensionalisation of ν t in r d [START_REF] Delville | La décomposition orthogonale aux valeurs propres et l'analyse de l'organisation tridimensionnelle des écoulements turbulents cisaillés libres[END_REF], leading to the following expression:

G νt = C 3 max(0, -∂ν t /∂n) √ u i,j u i,j κd w ( 25
)
to which the same function f d ( 24) is applied to obtain the second shielding function f d (G νt ). Other choices could be considered for the non-dimensionalisation of ∂ν t /∂n, for instance involving ν t . However, the proposed choice has the positive feature that the sensor decays with the wall distance (far away from walls, no boundary layers are expected) and it focuses on the interesting feature of the eddy viscosity field, i.e. its gradient rather than its level.

It should be emphasized that although it is a necessary condition, non-dimensionalisation is not sufficient to obtain the Reynolds-number invariance of the sensor. Using the local skin friction velocity u τ and boundary layer thickness δ, one may however write:

G νt = C 3 max(0, -∂(νt/(uτ δ)) ∂(n/δ) √ u i,j u i,j /(u τ /δ) κd w /δ (26)
Assuming that ν t /(u τ δ) and √ u i,j u i,j /(u τ /δ) are Reynolds-invariant functions of d w /δ in the outer layer at high Reynolds number, which is a common assumption, then it may be concluded that the proposed sensor is Reynolds-number invariant in the outer layer at high Reynolds numbers. This invariant behaviour of the sensor G νt will be verified in the next section where the constant C 3 will also be calibrated by considering boundary layer profiles at various Reynolds numbers. It should be noted that building on the non-dimensionalisation exposed here, it is already possible to guess the order of magnitude of C 3 , as detailed in Appendix B. The calibration of C 3 proposed in the next section will include non-zero pressure gradients and a value of C 3 will be retained for a general use as a case-independent constant.

One last feature of the proposed sensor G νt (25) is worth emphasizing. The fact that it detects the regions close to the wall were the eddy viscosity decreases with the wall distance makes it capable of detecting the edge of boundary layers where the eddy viscosity would have decayed in the outer layer for some reason. Indeed, this would bring a gradient of eddy viscosity closer to the wall than is expected from the flow physics, probably leading to a higher value of G νt which reinforces the shielding of the boundary layer in that region. This is the reason why the new sensor, although it also depends on the ν t field which it is supposed to shield, can be more efficient than the original sensor for which a decay of the ν t field is associated with a decrease of the shielding. Here, a decay of the ν t field in the outer layer induces an increase of the shielding just below. The practical consequence is that using this sensor in the new ZDES mode 2 (EP), the authors have observed that the eddy viscosity in some regions of unsteady flows could be sometimes reconstructed after having been previously destructed (see especially §4.5). This is contrasting with ZDES mode 2 (2012) where it is very rare that eddy viscosity destructed in boundary layers can be produced again in unsteady flows.

The proposed second shielding function based on G νt (25) might be active in some free shear layers if they are close enough to a wall. This may be typically the case behind a backward facing step. This motivates the introduction of an inhibition function able to detect detached flow regions in order to deactivate the second shielding function there.

Inhibition function

A first attempt to detect detached flow regions relied on an estimation of the fraction of resolved turbulence compared with the modelled part. The idea was that if detached flow regions have the potential for the resolution of turbulence,

i.e. if the flow instabilities are strong enough and the mesh fine enough, then some resolved turbulence should appear, which could be enough for the sensor to detect that the region should be treated in LES rather than in RANS. Although this strategy did work to some extent, it was not efficient enough to properly unlock resolved turbulence in all free shear layers. This is not really surprising since the release of LES content would rely on the presence of LES content, leaving room for situations where the field stays in RANS.

Since the objective is to detect detached flow regions, the most obvious choice is to detect vorticity concentration away from the walls. The new sensor is consequently based on the wall-normal gradient of the vorticity magnitude, ∂(|ω|)/∂n.

In most attached boundary layers, this gradient is negative since the vorticity magnitude is high near the wall and decays away from the wall. In the case of flow separation on the contrary, the bottom of the free shear layer is characterised by a positive gradient, ∂(|ω|)/∂n > 0, because the core of the free shear layer carries high values of vorticity. A simple idea is consequently to activate the inhibition function when ∂(|ω|)/∂n > 0. However, this also happens in some attached boundary layers close to separation, where a small positive gradient is encountered. In boundary layers close to separation but still attached, the inhibition function should not be activated. As a consequence, a threshold on ∂(|ω|)/∂n must be applied, large enough for all attached boundary layers to be shielded, and small enough for the inhibition function to activate in separated flows. Since the first objective of the method is to protect attached boundary layers, the threshold should be chosen so that the protection of attached boundary layers is guaranteed. In order to minimise the negative consequences on objective 3 (delay in the formation of instabilities), the threshold will however be no larger than needed.

The consequence is that ∂(|ω|)/∂n will be non-dimensionalised and the resulting sensor calibrated by considering attached boundary layers rather than free shear layers in typical flow separations.

Several quantities could be considered for non-dimensionalising ∂(|ω|)/∂n and obtaining the sensor G ω . The choice made is to use ν t and √ u i,j u i,j , because of the good properties of the resulting sensor:

G ω = ∂(|ω|) ∂n ν t √ u i,j u i,j
sensor when either the rate of strain or the vorticity vanishes. The sensor does not decay with the wall distance but this is not an issue since the shielding sensors do decay with d w , meaning that far away from walls the inhibition sensor will have little effect whatever its value.

In addition to the detection of separated mean flow regions, the function G ω is also sensitive to the presence of resolved turbulence in the instantaneous field. Indeed, resolved eddies close to the wall may contain regions where ∂(|ω|)/∂n > 0.

The inhibition of the second shielding function in such cases is a positive feature since the main objective of this function is to protect the wake layer, which may be resolved rather than modelled if resolved turbulence is present in this region (for instance during the merging between a boundary layer and a free shear layer located above it). It may be concluded that the aim of the sensor G ω is to detect both the detached flow regions and some resolved turbulence regions that would be shielded otherwise.

The inhibition function f R is applied to the sensor G ω in order to be 0 when the second shielding function should be deactivated, i.e. when G ω is above the chosen threshold. Conversely, f R should be 1 when the second shielding is kept, that is when G ω is below the threshold. The expression for f R connecting smoothly these two behaviours is the following:

f R (G ω ) =            1 if G ω ≤ C 4 1 1+exp -6α 1-α 2 if C 4 < G ω < 4 3 C 4 , with α = 7 6 C4-Gω 1 6 C4 0 if G ω ≥ 4 3 C 4 ( 28
)
where C 4 is a constant related to the threshold value of G ω . It has to be calibrated so that it is as small as possible while being large enough in order to protect boundary layers even close to separation. A value for C 4 will be proposed in the next subsection for a general use as a case-independent constant. The function f R (G ω ) is illustrated in fig. 3.

Calibration using a boundary layer code

In order to calibrate the constants C 3 and C 4 involved in the second shielding function [START_REF] Fröhlich | Hybrid RANS/LES methods for the simulation of turbulent flows[END_REF] and in the inhibition function [START_REF] Gand | Investigation of turbulence development in incompressible jets with Zonal Detached Eddy Simulation (ZDES) and synthetic turbulent inflow[END_REF] respectively, a set of boundary layer computations covering a wide range of Reynolds numbers and pressure gradients is considered. This set has been generated using the ONERA boundary layer code CLICET ( [START_REF] Aupoix | Couches Limites Bidimensionnelles Compressibles. Descriptif et Mode d'emploi du Code CLICET -Version[END_REF]). This code solves the compressible boundary layer equations with a self adaptive mesh refinement algorithm ensuring proper grid convergence. It has been used in studies such as [START_REF] Aupoix | Extensions of the Spalart-Allmaras turbulence model to account for wall roughness[END_REF][START_REF] Aupoix | Extension of Lysak's Approach to Evaluate the Wall Pressure Spectrum for Boundary Layer Flows[END_REF][START_REF] Renard | A theoretical decomposition of mean skin friction generation into physical phenomena across the boundary layer[END_REF].

The set of boundary layers considered is based on the following law for the outer velocity: Firstly, values slightly above zero are observed in the buffer layer. This is explained by the behaviour of the r d sensor [START_REF] Delville | La décomposition orthogonale aux valeurs propres et l'analyse de l'organisation tridimensionnelle des écoulements turbulents cisaillés libres[END_REF] which is designed to be close to 1 in the logarithmic layer thanks to its dependence on ν t and to increase with decreasing wall distances in the viscous sub-layer so that it reaches high values very near the wall thanks to its dependence on ν. Indeed in the buffer layer, when getting closer to the wall, the quick decay of ν t is not immediately compensated by the increasing importance of viscous effects, leading to a local minimum of r d which explains why f d (r d ) is slightly above zero there. However, this is not an issue since the maximum value reached by f d (r d ) in the buffer layer is not large at all and no impact on the mean velocity profile in the buffer layer in methods such as DDES has ever been observed or reported, to the authors' best knowledge. It should also be noted that the maximum value of f d (r d ) in the buffer layer displayed in fig. 6 is well below the recommended threshold 1 -f d0 = 0.2 above which the enhanced destruction in grey areas is activated (see equation [START_REF] Deck | Large scale contribution to mean wall shear stress in high Reynolds number flat plate boundary layers up to Re θ = 13650[END_REF]).

u e = u 0 x -x 0 -x 0 m (29
The second region where the original shielding function f d (r d ) is above zero is the outer part of the wake layer, as illustrated by fig. 6. This second issue is not benign, contrary to the first issue. Indeed, the shielding function f Based on the observations of fig. 6, the objective of a better shielding of attached boundary layers is to protect the outer part of the wake layer in order to have a shielding function lower than approximately 0.05 across the full boundary layer profile. The simplest idea is to increase the C 1 constant involved in the original shielding function [START_REF] Delville | La décomposition orthogonale aux valeurs propres et l'analyse de l'organisation tridimensionnelle des écoulements turbulents cisaillés libres[END_REF][START_REF] Dong | An enhanced version of delayed detached-eddy simulation based on the v 2 -f model[END_REF], and this is worth investigating for a better understanding of the problem, although the literature review in the introduction indicates that this cannot lead to a method satisfying all objectives set to the present study. This simple attempt is illustrated in fig. 7 where the original shielding function modified with C 1 = 21 is displayed. This indicates that such a simple modification can be successful in favorable and zero pressure gradients but it is not sufficient to shield the outer part of the wake layer in the adverse pressure gradient cases, where high values (close to 1) of the shielding function are still observed. A value of C 1 significantly higher than 21 would consequently be needed. This modification cannot be satisfying since the considered value (C 1 = 21) is already much larger than the values experimented in the literature (typically a maximum C 1 = 16 since high values have been reported to lead to an excessive delay in the formation of instabilities in free shear layers, as reminded in the introduction). This demonstrates that introducing a new shielding function is necessary to reach all the objectives set for the new method.

In order to enhance the shielding of the attached boundary layers, the second shielding function introduced in §3.1 is combined with the original shielding function to obtain a protection

f P = f d (r d )f d (G νt ), as motivated in §2.2, whereas
G νt is defined in eq. ( 25). The constant involved is set to C 3 = 25 in figure 8. This value is larger than the value C 3 = 9 estimated for ZPG boundary layers in §3.1.2 because the APG cases are expected to require a stronger shielding than the ZPG boundary layers. This is confirmed by fig. 8 where the outer part of the wake layer is shielded with f P very close to 0 except in the strongest adverse pressure gradient case where f P is slightly above 0, justifying the need for a The second constant that needs calibration is the constant C 4 involved in the inhibition function f R (G ω ) [START_REF] Gand | Investigation of turbulence development in incompressible jets with Zonal Detached Eddy Simulation (ZDES) and synthetic turbulent inflow[END_REF]. As introduced in §3.1.3, the inhibition function should activate in some free shear layers were its role will be to deactivate the second shielding function. The value of the constant C 4 must be calibrated so that the inhibition function does not activate in attached boundary layers (where the second shielding function must obviously not be deactivated). In presence of adverse pressure gradients, the velocity profile of an attached turbulent boundary layer sometimes features inflection points and regions where G ω > 0, typically one such region very close to the wall and one further away from the wall.

The G ω > 0 region closest to the wall is not problematic here since the original shielding function f d (r d ) does not fail in the vicinity of the wall, so that deactivating the second shielding function if the inhibition function activates will not have any consequence. A G ω > 0 region further away from the wall can on the other hand be problematic, and this is why a constant C 4 strictly greater than zero is needed. More precisely, C 4 should be the smallest value possible that ensures that G ω is smaller than C 4 in attached boundary layers. Plotting G ω in fig. 9 suggests that C 4 = 0.03 is a satisfying value since all profiles of G ω remain well below this threshold including for the strongest adverse pressure gradient case considered, which is close to separation. One may also notice in fig. 9 that the inhibition function should not have any tendency to activate near the edge of the boundary layer since G ω takes very large negative values in this region. This is a very favourable property since one of the aims of the second shielding function is to prevent corruption of the eddy viscosity near the edge of the boundary layer by the low SGS values outside the boundary layer, so that the second shielding function must not be deactivated near the edge of the boundary layer. Overall, as illustrated by fig. 9 the choice C 4 = 0.03 introduces a reasonable safety margin to ensure the proper protection of attached boundary layers. In case of flow separation (of physical origin), the activation of the inhibition function will then be enabled thanks to this not too high value of C 3 in order to deactivate the second shielding function and unlock instabilities leading to resolved turbulence as desired for separated flows. The choice made here is consequently to switch as much as possible into LES for any detached flow whereas attached boundary layers are fully protected. This means that even mild flow separations will be treated as much as possible in LES. Another option could have been to treat mild separations in RANS and reserve LES for massive separations only, which would have been closer to the choices apparently made by Menter for his SDES/SBES methods [START_REF] Menter | Stress-Blended Eddy Simulation (SBES) -A New Paradigm in Hybrid RANS-LES Modeling[END_REF]. However, the limitations of RANS models to predict accurately mild flow separations and the relatively successful behaviour of methods such as ZDES mode 2 to predict mild separations with resolved turbulence in former studies encourage to choose a LES switch occurring in all detached flows. Concerning the very accurate prediction of mild flow separations, it is the authors' opinion that the safest option may be to resort to a WMLES strategy in the concerned regions, i.e. use mode 3 of ZDES in the mildly separated zones and mode 2 only in regions where massive separations occur.

Coming back to figure 9, with the recommended value C 4 = 0.03 for the inhibition function f R (G ω ), the latter will never activate (i.e. it will stay equal to 1) in the whole dataset considered. For this reason, the full protection function 18) applied to the present dataset is equivalent to the product of shielding functions The proposed shielding and inhibition functions may be used in any hybrid RANS/LES context as soon as the eddy viscosity field is available (i.e. whatever underlying RANS model is used, provided it predicts an eddy viscosity level, and whatever RANS/LES strategy is adopted). The formulation of the enhanced ZDES is detailed in Appendix A with the constant values and complete sensor expressions defining the new method ZDES mode 2 (EP). The generality of the formulation in three-dimensional flows will be examined in the next subsection considering a few flow configurations for a priori tests of the sensors. Simulations using the enhanced ZDES mode 2 (EP) will then be presented in section 4 for evaluation of the method on four test cases to determine to what extent the new shielding is much stronger than the original one, and at what cost in terms of delay in the formation of instabilities in free shear layers this has been obtained.

f P = f d (r d )(1 -(1 -f d (G νt ))f R (G ω )) (
f d (r d )f d (G νt )

Confirmation of the calibration based on a priori tests involving three-dimensional flows

The constants involved in the new shielding and inhibition functions have been calibrated in section 3.2 based on a wide set of two-dimensional flat plate boundary layers. The generality of the formulation in three-dimensional flows around curved walls will be assessed next by means of a priori tests of the functions evaluated on RANS fields. Three configurations are considered, namely a corner flow in a supersonic square duct, the cross-flow over a diamong wing and a turbulent separation at the trailing edge of a NACA 4412 airfoil.

Corner flow in a supersonic square duct

The first a priori test is based on the 3D supersonic square duct 3DSSD validation case of the NASA turbulence modeling resource database1 which relies on the experiment by Davis & Gessner [11]. This square duct flow includes four turbulent boundary layers spatially growing until they merge into a fully developed duct flow. The inflow Mach number is 3.9 and the flow develops over 50D, where D is the height and width of the cross-section. The a priori test is based on a RANS simulation with the Spalart-Allmaras turbulence model [START_REF] Spalart | A one equation turbulence model for aerodynamic flows[END_REF] and the QCR2000 closure relation [START_REF] Spalart | Strategies for turbulence modelling and simulations[END_REF], so that corner vortices are captured. The main purpose of this test is indeed to check that the vorticity associated to the secondary circulation near the corners does not activate the inhibition function. The whole corner area should be shielded to ensure a RANS treatment. .14). One may conclude that the new shielding is satisfying and that the inhibition function never induces a reduction of the shielding in this a priori test case, suggesting that the secondary circulation does not trigger an undesired switch to LES, leading to the whole corner flow being treated in RANS as expected from ZDES mode 2 in this case.

Cross-flow over a diamond wing

The next a priori test case is the flow over a diamond wing at Mach 0.

NACA 4412 airfoil trailing edge separation

The ability of the shielding and inhibition functions to detect a shallow flow separation is tested on a NACA 4412 airfoil trailing edge separation inspired by the 2DN44 validation case of the NASA turbulence modeling resource database2 but with different conditions. We consider here a NACA 4412 airfoil of chord c = 0.128 m at a Mach number of 0.177 simulated with the Spalart-Allmaras turbulence model [START_REF] Spalart | A one equation turbulence model for aerodynamic flows[END_REF] for three different conditions introduced in order to include a

Reynolds number effect and an angle of attack effect: the effect of the activated inhibition function in the detached flow region. This illustrates the capacity of the method to ensure RANS shielding until separation and LES switch in separated zones, as expected from ZDES mode 2.

The success on these a priori tests confirms the calibration of the constants of the model which is evaluated in the next section by running actual simulations involving the enhanced ZDES mode 2.

4. Assessment of enhanced ZDES mode 2

Overview

The proposed enhancement of ZDES mode 2 (especially functions A. with the standard value C 1 = 8 and with a higher value C 1 = 21 as identified in §3.2, and finally the enhanced method ZDES mode 2 (EP) proposed here. This latter method is assessed on the fourth test case consisting of transonic buffet over the supercritical OAT15A airfoil. This is a very challenging test case since it features a shock-wave / boundary layer interaction including a smooth-wall flow separation moving back and forth over time, which means that the model must be able to switch from LES back to RANS when the shock is moving into the downstream direction. All ZDES calculations evaluate the vorticity-based cell size ∆ ω rather than volume-based. The flow solver used (FLU3M, developed by ONERA) is a compressible finite-volume multi-block structured solver with second-order accuracy numerical schemes and has been used and validated on a wide range of both academic and industrial studies. More details on the flow solver may be found in [START_REF] Deck | Recent improvements of the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF].

4.2. Test case #1: spatially developing turbulent boundary layer

Test case

This first test case is a spatially developing zero-pressure-gradient turbulent boundary layer over a smooth flat plate of length L = 1 m. The free stream velocity is U 0 = 683 m.s -1 , the static pressure is set to P 0 = 19820 P a, the temperature equals 288K leading to a Reynolds number per meter Re = 9.15 × 10 6 m -1 and a freestream Mach number M 0 = 2. The choice of a supersonic boundary layer is a good test for robustness of the method. It should however be noted that results similar to those presented below have been obtained in a low subsonic boundary layer (not presented here for the sake of conciseness).

The streamwise evolution of the salient Reynolds numbers together with the major parameters of the grid resolution are summarized in figure 14. The first cell is at y + = 1 in the framework of a cell-centred code (i.e. the first vertex is at

y + = 2)
. The grid is designed to have a classical WMLES resolution for 0.35 ≤ x/L ≤ 0.8 and represents a challenging case for the RANS shielding of the boundary layer in a classical hybrid RANS/LES framework. The grid distribution in the streamwise direction is not uniform in order to maintain ∆x/δ ≈ 0.1 (or 100 + ≤ ∆x + ≤ 200 + ) in the WMLES region.

Conversely, the grid in the spanwise direction is constant ∆z + = 50. For x/L > 0.8, mesh cells are stretched again. The total number of points in the streamwise, spanwise and wall-normal directions are respectively N x = 713, N y = 127 and N z = 31 leading to a total number N xyz = 2.81 × 10 6 points. The physical time-step is set to ∆t CF D = 2.4 × 10 -7 s.

Enough time iterations are performed for the simulation to reach a stationary state which is analyzed next.

The sensitivity of the original DDES (2006) (with standard constants, especially C DES = 0.65) to the grid resolution has been illustrated in the introduction, where fig. 2 presents results obtained in the present test case. This suggests that further grid refinement must be included in the test, and the choice made is to simulate an infinitely fine mesh in order to check that the shielding can withstand any mesh refinement. This is achieved by setting the constant C DES = 0, which emulates a vanishing grid cell size. In the following, the methods will be tested with C DES = 0, which is a very demanding test whose successful outcome gives confidence that the new model is mesh-proof. In the compressible context of the test case, the results are rescaled as introduced by Hopkins & Inouye [START_REF] Hopkins | An evaluation of theories for predicting turbulent skin friction and heat transfer on flat plates at supersonic and hypersonic mach numbers[END_REF], especially the mean skin friction coefficient

C f = τ w /(1/2ρ ∞ U 2 ∞
) is reduced to the incompressible case C f i using the Van Driest II approach:

C f i = F c C f with F c = T w T ∞ -1 arcsin 1 - T ∞ T w 1/2 -2 (30)
which can be compared at a given Re δ2 = θU ∞ ρ ∞ / µ w with θ =

Results and discussion

The evolution of the mean skin friction coefficient displayed in fig. 15 2012) with an increased constant C 1 = 21 provides a prediction of skin friction identical to the RANS level. This is consistent with the information provided by fig. 7 which indicates that in a case with zero-pressure gradient like here, the choice of C 1 = 21 provides an almost perfect shielding of the entire boundary layer profile. It is however reminded that fig. 7 suggests that the shielding might fail in presence of a strong adverse pressure gradient. This is not the case of the proposed new ZDES mode 2 (EP) for which fig. 8 indicates that the shielding should be able to resist strong adverse pressure gradients. In the case of the very demanding test displayed in fig. 15, ZDES mode 2 (EP) does predict the same skin friction as the RANS simulation. This confirms that the new shielding is very robust since the grid size seen by the model is zero (because of C DES = 0).

In order to confirm that the new shielding is fully satisfying, turbulent profiles are displayed in figure 16, which displays the velocity and eddy-viscosity profiles at station X/L = 0.3 (respectively X/L = 0.7). At both stations, it can be seen that the original DDES (2006) and ZDES mode 2 (2012, original C 1 = 8) feature an almost vanishing eddy viscosity profile and a velocity profile akin to a laminar boundary layer. On the contrary, the original ZDES mode 2 with a modified constant C 1 = 21 and most interestingly the new ZDES mode 2 (EP) predict eddy viscosity and mean velocity profiles virtually indiscernible from the RANS profiles. This result gives a strong confidence in the ability of the new ZDES mode 2 (EP) to properly shield attached boundary layers since even an infinite grid refinement simulated by setting C DES = 0 does not alter the turbulent profiles. Together with the calibration performed in §3.2 on a dataset of boundary layers and the a priori 3D tests of section 3.3, these results provide strong elements suggesting that ZDES mode 2 (EP) reaches objective #1 of the introduction (shielding protecting the whole attached boundary layer profile including on infinitely fine meshes with pressure gradients). The following two test cases will for this reason be devoted to the other objective that has to be demonstrated, i.e. objective #3 (minimum delay in the formation of instabilities in free shear layers), on meshes not particularly fine since in the case of this objective, a coarse mesh is a more demanding test than a very fine mesh as needed for the evaluation of objective #1.

Test case #2: spatially developing mixing layer 4.3.1. Test case

The spatially developing mixing layer permits to study the behaviour of the hybrid subgrid-scale operating mode and especially to assess if the new proposal does not delay too much the switch into LES mode. The free shear flow studied here is generated at the confluence of two boundary layers developing from each side of a splitter plate which is modelled as a flat plate with zero thickness (see figure 17). The physical conditions correspond to the ones studied by Delville the high and low speed sides respectively). These numerical parameters (location of transition, grid, size of the domain illustrated in figure 17) are those proposed by Deck [START_REF] Deck | Recent improvements of the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF] and were later retained in the EU-funded Go4Hybrid project [START_REF] Peng | Go4Hybrid: Grey Area Mitigation for Hybrid RANS-LES Methods -Free Shear Layer[END_REF] because they offered a challenging test case for the grey area mitigation and enabled comparison with experimental data at an affordable computational cost.

The basic two-dimensional grid contains N xy ≈ 197 000 points (including 170 000 points for the mixing layer i.e. for

x > 0) and has been duplicated in the spanwise direction leading to a total number of points N xyz = 19.10 6 with N z = 96.

Near 30 points are clustered in the initial vorticity thickness. The time step is set to ∆t CF D = 10 -6 s.

Results and discussion

A first qualitative insight into the resolved turbulence in the mixing layer is provided by fig. 18 where a Q criterion is plotted. As expected, the original DDES (2006) simulation suffers from a severe delay in the formation of instabilities, with large unphysical almost-2D rollers developing downstream and no properly 3D resolved turbulence visible. This is caused by the excess of eddy viscosity in the wake of the splitter plate as illustrated in the inserts. The behaviour obtained is completely unphysical since the solution displays the behavior of a laminar shear layer: growth of a two-dimensional Kelvin Helmholtz instability followed by vortex pairing. Improvements over the baseline DDES exist, using for instance a shear layer adapted subgrid length-scale [START_REF] Shur | An Enhanced Version of DES with Rapid Transition from RANS to LES in Separated Flows[END_REF], but it should be reminded that such enhancements of the RANS/LES A more quantitative look at the results is given in fig. 19 where the vorticity thickness is classically defined as δ ω = (U high -U low )/ max y (∂U/∂y). Because of its extremely important delay in the formation of instabilities, the DDES (2006) simulation strongly underpredicts the streamwise growth of δ ω . All presented hybrid RANS/LES calculations underpredict the initial growth of δ ω near the trailing edge. As shown in [START_REF] Deck | Recent improvements of the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF], this could be improved if resolved turbulence were present in the incident boundary layers, using ZDES mode 3, but this is out of the scope of the present study where the capacity of the hybrid RANS/LES method to quickly recover a resolved LES content from a RANS upstream region is considered, which is vital for an automatic use of the methods with no manual injection of turbulence or RANS/LES manual switching. After the initial relaxation, the streamwise evolution of δ ω is similar to the experimental and RANS behaviours for all ZDES mode 2 tests considered. In fig. 19, there is no significant difference for δ ω between the original ZDES mode 2 (2012, C 1 = 8), ZDES mode 2 (2012, C 1 = 21) and the proposed ZDES mode 2 (EP). This quantitatively confirms the visual impression obtained in fig. 18 that the proposed new ZDES mode 2 (EP) does not induce a significant deterioration of the delay in the formation of fully resolved 3D turbulence in the present free shear layer.

A more detailed analysis of the turbulent profiles is provided in fig. 20 In order to fully confirm the good properties of the new ZDES mode 2 (EP), a case where a free shear layer is located close to walls is considered in the next section devoted to a backward facing step flow. This makes the detection of the free shear layer more difficult since there are similarities between the wake layer in attached flows and a separated shear layer (see [START_REF] Coles | The law of the wake in the turbulent boundary layer[END_REF]). For this reason, the inhibition function detecting flow separation will be crucial in preventing part of the separated region to be mistakenly detected as the wake region of an attached boundary layer. The following test case is consequently seen as mandatory for the development and validation of a new hybrid RANS/LES method. Testing the method in a backward facing step flow will show its capacity to resolve turbulence in a free shear layer in spite of its proximity to walls. This especially relies on the proper detection of flow separation using the inhibition function. An efficient destruction of the eddy viscosity in the resulting grey area (switch from RANS to LES) is also needed for success in this test case. It should be noted that during the development of the proposed new ZDES mode 2

(EP), it was found rather easy to succeed in the first two test cases (shielding of the boundary layer and resolution of the mixing layer) whereas obtaining properly resolved turbulence in the BFS flow with no significant delay in the apparition of LES content while keeping the very strong shielding demonstrated in test case #1 proved very challenging.

The chosen test case is a backward facing step of height h = 35 mm corresponding to the A3C combustion chamber which has been experimentally studied by Moreau et al. [START_REF] Moreau | Experimental and Numerical Study of a Turbulent Recirculation Zone with Combustion[END_REF] and described in detail by Sainte-Rose et al. [START_REF] Sainte-Rose | A DES method applied to a Backward Facing Step reactive flow[END_REF]. With a freestream velocity U 0 = 50 m/s and a static temperature of 520 K, the Reynolds number is Re h = 40 000. The ratio between the boundary layer thickness at the separation point and the step height is δ/h = 0.37. The grid counts N x = 803, N y = 136 and N z = 36 points in the streamwise, vertical and spanwise directions respectively, with ∆z/h = 0.057. The time step selected is ∆t CFD = 10 -6 s, i.e. ∆t CFD = 1.42 • 10 -3 h/U 0 . It should be noted that the mesh is not especially fine, which is justified by the objective of the present test case which is not designed to check the shielding capability of the method but rather to evaluate the possibility to recover quickly a LES content on a somewhat coarse mesh resolution typical of what is affordable in an industrial context. The transonic buffet is an aerodynamic phenomenon that results in a large-scale self-sustained motion of the shock over the surface of the airfoil. The numerical simulation of transonic buffet is a highly challenging and difficult case for hybrid RANS/LES methods because it presents a thin-layer separation that moves significantly in time at a time-scale that is much slower than the one of the wall-bounded turbulence. This is the most complex case, and a good outcome would be a positive indicator for the robustness of the new version of ZDES mode 2 intended to treat industrial applications.

Results and discussion

The supercritical airfoil OAT15A is computed here in the same flow conditions as in the experiment by Jacquin et al [START_REF] Jacquin | Experimental study of shock oscillation over a transonic supercritical profile[END_REF]. The profile has a chord equal to c = 230mm and a relative thickness of 12.3%. Transition is fixed near the leading edge at x/c = 7% on both sides of the airfoil. The Reynolds number, based on the freestream velocity and the chord, is equal to 3 × 10 6 , whereas the freestream Mach number is set to 0.73 and the angle of attack is equal to 3.5 deg.

The basic two-dimensional grid in (x, y) plane is based on a C-H topology and contains 403 000 points. The threedimensional grid is built by duplicating the two-dimensional grid in the spanwise direction; a constant grid extension is considered ∆z/c = 0.00139 using N z = 145, which leads to L z /c = 0.2 and a total number of 58.4 × 10 6 points.

The corresponding grid resolution in the spanwise (respectively streamwise) direction using classical wall-unit scaling is while the whole simulation lasted over 500 000 CPU hours.

Results and discussion

The salient features of the instantaneous flow field are highlighted by figure 27, which shows the isosurface of the Q criterion coloured by the streamwise velocity component for two extreme locations of the shock. The location of the shock is evidenced by plotting an iso-surface of pressure. One can first note that no resolved turbulent content can be seen in the flow upstream of the shock as this boundary layer is treated in URANS mode within ZDES mode 2. Downstream of separation, ZDES mode 2 switches rapidly in its scale resolving mode as no delay in the LES content formation can be depicted. To get another insight on how the present ZDES mode 2 operates, figure 28 displays the instantaneous eddy viscosity field for the two aforementioned locations of the shock. It is worth noting that no spurious drop of eddy The salient features of pressure field are then displayed in figure 30. The pressure coefficient C p = P -P0 1 2 γM 2 0 P0 distribution is compared with experimental data in figure 30(a). The motion of the shock is evidenced by the spreading aspect of the shock pressure distribution. This aspect is well reproduced by ZDES mode 2 together with the pressure plateau in the supersonic area and the pressure levels at the blunt trailing edge, the latter being a crucial point to capture the circulation around the airfoil. Nevertheless it is observed that the separation is located somewhat too upstream. Such a behavior has also been observed in WRLES simulations (see [START_REF] Garnier | Large-eddy simulation of transonic buffet over a supercritical airfoil[END_REF]) and can partly be attributed to the short duration of the calculation. The shock motion is not perfectly periodic in the experiment where signals were collected over 2300 periods. In the present calculation only about 8 periods were acquired due to the high CPU cost. The existence of buffet cells [START_REF] Plante | Similitude between 3D cellular patterns in transonic buffet and subsonic stall[END_REF][START_REF] Iovnovich | Numerical study of shock buffet on three-dimensional wings[END_REF] which cannot be captured by the calculation as due to the limited spanwise extent of the grid (L z /c = 0.2 here) is an other source of discrepancy with the experimental set-up affecting the numerical prediction of the pressure field. [7] has been used. Indeed, this method is well adapted to study short data that are known to consist of sinusoids in white noise (e.g. [START_REF] Trapier | Experimental study of supersonic inlet buzz[END_REF][START_REF] Deck | Numerical simulation of transonic buffet over a supercritical airfoil[END_REF][START_REF] Larchevêque | Large-eddy simulation of a compressible flow in a threedimensional open cavity at high reynolds number[END_REF]). The first sensor is located at the upper side close to the trailing edge at x/c = 0.9. The periodic nature of the large-scale self-sustained motion of the shock is evidenced by the strong harmonics peaks in the spectra.

On the experimental side, the main peak is at 69 Hz (or f c/U 0 = 0.065) and is very well reproduced by the present calculations together with the broadband spectral distribution at higher frequencies. The second sensor is located at the lower-side of the airfoil where ZDES mode 2 operates in URANS mode owing to its design so that only the high amplitude and low frequency fluctuations can be captured by the calculation.

A deeper investigation of the flow is beyond the scope of this paper (see for instance [START_REF] Deck | Numerical simulation of transonic buffet over a supercritical airfoil[END_REF][START_REF] Jacquin | Experimental study of shock oscillation over a transonic supercritical profile[END_REF][START_REF] Fukushima | Wall-modeled large-eddy simulation of transonic airfoil buffet at high reynolds number[END_REF]) but this calculation shows the robustness of the new ZDES mode 2 to preserve the (U)RANS behaviour upstream of separation (even on a very fine grid) in a highly challenging configuration namely when the location of separation moves in time.

Conclusion

A strategy to better protect attached boundary layers which should be treated in RANS mode in a hybrid RANS/LES context has been introduced. It addresses two issues typically encountered by hybrid RANS/LES methods such as the original DDES (2006) [START_REF] Spalart | A new version of Detached-Eddy Simulation, resistant to ambiguous grid densities[END_REF] or the original ZDES mode 2 (2012) [START_REF] Deck | Recent improvements of the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF] determining automatically whether a given region of the flow is treated in RANS or in LES. The first issue is the RANS shielding of attached boundary layers which should hold whatever mesh refinement and pressure gradient are encountered. The original shielding function shared by DDES (2006) and ZDES mode 2 (2012) has been reported in recent years to fail in some cases (fine mesh, adverse pressure gradient), motivating the present study. The second issue is the difficulty to quickly obtain a resolved LES content after a RANS-to-LES switch of the method, typically in free shear layers. Improving the RANS shielding of attached boundary layers without delaying too much this RANS-to-LES transition is extremely challenging and cannot be achieved by only tuning a constant in the original shielding function.

To address these issues, four objectives have been defined for the development of a new version of ZDES mode 2.

Objective #1 stipulates that the attached boundary layers should be fully shielding including on infinitely fine meshes with arbitrary pressure gradients. According to objective #2, the new method should never protect attached boundary layers less than what the original ZDES mode 2 (2012) did. This means that the new method will always be safer (in terms of RANS shielding) than the original method. Despite the strongly increased shielding, according to objective #3 there should be a minimum delay in the formation of instabilities in free shear layers, so that the new method should be almost as efficient as the original method to resolve LES content near RANS-to-LES transition regions.

Figure 1 :

 1 Figure 1: Main issues in "DES-type" methods. Issue (I): insufficient protection of the attached boundary layer region when C DES ∆ < δ, Issue (II): longitudinally induced grey-area leading to a possible delay in the formation of "LES-content".

Figure 2 :

 2 Figure 2: Streamwise evolution of the skin friction coefficient together with the streamwise mesh size ∆x/δ in a flat plate boundary layer test case (further described in §4.2).

3 . 4 .

 34 Minimum delay in the formation of instabilities in free shear layers: the new ZDES mode 2 (EP) should recover a proper LES-content almost as rapidly as the original ZDES mode 2 (2012) despite the strongly increased shielding. Full compatibility of ZDES mode 2 (EP) with ZDES modes 1 and 3: the resulting subgrid scale model in LES regions should be the same for all three modes. This ensures a continuity of treatment of resolved turbulence across zones treated with different modes of ZDES.

  Shallow separation or attached flow Location of separation Known in advance (geometry) Unknown a priori (APG) Unknown Treatment of the attached TBL U RAN S U RAN S W M LES or W RLES

Figure 3 :

 3 Figure 3: Function f R (Gω) (28) with C 4 = 0.03.

) where u 0

 0 fig.4, illustrating the diversity of the dataset and the proximity of the most adverse pressure gradient case with respect to flow separation.

Figure 4 :

 4 Figure 4: Inner-scaled mean velocity profiles of the dataset of turbulent boundary layers considered for calibration of the model constants.

Figure 5 :

 5 Figure 5: Outer-scaled wall distance as a function of the outer-scaled mean velocity in the turbulent boundary layers considered for calibration of the model constants.

  d (r d) reaches high values (close to 1) in this region, implying an actual failure of the shielding. This is especially true near the edge of the boundary layer, which explains the sensitivity of methods such as DDES to the cell size near the edge of the boundary layer and the tendency of the eddy viscosity field to be corrupted by LES SGS values coming from just outside the boundary layer. Besides, considering fig.6, it clearly appears that the shielding of attached boundary layers by the original function f d (r d ) is incomplete in the outer part of the wake layer irrespective of the cell size since the shielding function does not depend on the cell size. What depends on the cell size is the gravity of the consequences of the deficient shielding: the finer the mesh in the region not properly shielded and in the region just above, the stronger the decay of the eddy viscosity (because of local destruction and diffusion of low eddy viscosity values from the region above). As far as the pressure gradient effect is concerned, fig.6also reveals that the failure of the shielding is worse for stronger adverse pressure gradients. This explains why failure of shielding in methods such as DDES has been more frequently observed in presence of adverse pressure gradients than in zero pressure gradient boundary layers, which is very detrimental because of the risk of grid induced separation in such flows. It should however be emphasized that in the zero pressure gradient case the last ten percent of the velocity profile are not properly shielded by f d (r d ) which is already a large region, and even with favorable pressure gradients the shielding fails near the edge of the boundary layer.

Figure 6 :

 6 Figure 6: Original shielding function f d (r d ) (eq. (23-24)) applied to the turbulent boundary layers dataset and displayed as a function of the outer-scaled mean velocity.

Figure 7 :

 7 Figure 7: Original shielding function f d (r d ) (eq. (23-24)) modified with C 1 = 21 applied to the turbulent boundary layers dataset and displayed as a function of the outer-scaled mean velocity.

Figure 8 :

 8 Figure 8: New shielding function f P = f d (r d )f d (Gν t ) (25) with C 3 = 25 applied to the turbulent boundary layers dataset and displayed as a function of the outer-scaled mean velocity.

  already plotted in figure8. The final shielding proposed here with the recommended constants C 3 = 25

Figure 9 :

 9 Figure 9: Inhibition function sensor Gω applied to the turbulent boundary layers dataset and displayed as a function of the outer-scaled mean velocity (detail on the right hand side).

Figure 10 :

 10 Figure 10: A priori test of the shielding and inhibition functions on a corner flow in a supersonic square duct. Colormap: difference f Pf d (r d )f P 2 , pink isoline: f d (r d ) = 0.2, black isoline: f P = 0.2, violet isoline: f R = 0.8 (functions as detailed in equations A.4, A.7, A.10), white streamlines based on the in-plane velocity components.

15 and

 15 Re c = 4.11 • 10 6 based on the root chord c = 1.2 m at a 12 deg angle of attack as described in [18]. The 53 deg sweep angle of the leading edge of the diamond wing induces a significant crossflow and separation of a vortex sheet, as illustrated in fig. 11, making it a demanding 3D test case for the shielding of attached boundary layers and LES switch in detached flow regions. The a priori tests are performed on the basis of a RANS simulation using the Spalart-Allmaras turbulence model[START_REF] Spalart | A one equation turbulence model for aerodynamic flows[END_REF] with the rotation correction by Dacles-Mariani[START_REF] Dacles-Mariani | Numerical/Experimental Study of a Wingtip Vortex in the Near Field[END_REF] (SA-R variant with C rot = 4 in NASA terminology). More details may be found in[START_REF] Deck | Zonal Detached Eddy Simulation (ZDES) of the flow around the AVT-183 diamond wing configuration[END_REF].

Figure 11 :

 11 Figure 11: Skin friction lines and eddy viscosity fields of the SA-R RANS calculation of the diamond wing flow. Figure adapted from [18].

Figure 12 :

 12 Figure 12: A priori test of the shielding and inhibition functions in the flow over a diamond wing. Colormap: difference f P -f d (r d )f P 2 , pink isoline: f d (r d ) = 0.2, black isoline: f P = 0.2 (functions as detailed in equations A.4, A.9, A.12).

• a )

 a Re c = 2.28 • 10 6 , AOA = 12 deg, fully turbulent • b) Re c = 2.28 • 10 6 , AOA = 13.87 deg, transition forced at x/c = 0.025 on the upper surface and x/c = 0.103 on the lower surface • c) Re c = 2.28 • 10 5 , AOA = 12 deg, transition forced at x/c = 0.025 on the upper surface and x/c = 0.103 on the lower surface As illustrated by fig. 13, the three flow conditions lead to different separation bubble sizes, but the shielding and inhibition functions have the same qualitative behaviour in all cases, confirming their generality. In all cases, the new protection function f P (black isolines) (A.4,A.9, A.12) shields the whole attached boundary layer profiles, even in the thick boundary layer submitted to an adverse pressure gradient shortly before separation on the upper surface where the original shielding function f d (r d ) (pink isolines) fails to protect more than half of the profile. No undesired activation of the inhibition function (A.12) is seen upstream the separation point, whereas the f P -f d (r d )f P 2 field (colormap) indicates

(a) Rec = 2 . 28 •Figure 13 :

 22813 Figure 13: A priori test of the shielding and inhibition functions in a NACA 4412 airfoil trailing edge separation. Colormap: difference f P -f d (r d )f P 2 , pink isoline: f d (r d ) = 0.2, black isoline: f P = 0.2 (functions as detailed in equations A.4, A.9, A.12).

Figure 14 :,

 14 Figure 14: Turbulent Boundary Layer test case. Grid and streamwise evolution of relevant local Reynolds numbers Re θ = ρ∞U∞θ µ∞ , Re δ2 =

∞ 0 ρ u ρ∞U∞ 1

 1 u U∞ dy. The mean velocity profiles are Van Driest transformed according to u * = u 0 ρ / ρ w du and plotted as a function of y + = y u τ ρ w / µ w .

Figure 15 :

 15 Figure 15: Streamwise evolution of the skin friction coefficient (30) in the TBL test case.

  confirms the very strong failure of both the original DDES (2006) and original ZDES mode 2 (2012, original C 1 = 8) to shield the attached boundary layer under the infinite grid refinement simulated by setting C DES = 0. The skin friction is very low and the boundary layer almost relaminarizes. Conversely, the original ZDES mode 2 (

Figure 16 :Figure 17 :

 1617 Figure 16: Velocity and eddy viscosity profiles at stations x/L = 0.3 and x/L = 0.7.

[ 23 ]

 23 , namely U high = 41.54 m/s and U low = 22.4 m/s on both sides of the flat plate. At X = -0.01m, the boundary layers are turbulent with momentum thickness θ high = 10 -3 m and shape factor H high = 1.35 on the high speed side and θ low = 0.73 • 10 -3 m and H low = 1.37 on the low-speed side. The Reynolds numbers based on the momentum thickness are respectively equal to Re θ high = 2900 and Re θ low = 1200. To capture the correct boundary-layer profiles at the trailing edge, the flat plate has a length of 0.82 m on the high-speed side and 0.46 m on the low speed side and fine tuning has been done by varying the laminar-turbulent transition location (X transition = -0.708 m and X transition = -0.388 m on

  transition can also deteriorate the attached boundary layer shielding. Concerning the ZDES technique, the original mode 2 (2012, C 1 = 8) features a much faster growth of instabilities leading quickly to resolved 3D turbulence. As illustrated by the Q criterion in fig.18, the behaviour of the original ZDES mode 2 (2012) modified with C 1 = 21 is qualitatively similar to the original ZDES mode 2 (2012, C 1 = 8). The eddy viscosity in the wake of the splitter plate is not destructed as fast in the C 1 = 21 simulation as in the C 1 = 8 baseline simulation, but the impact in the resolved turbulence seen in fig.18is marginal. This can be explained by the relatively large distance between the free shear layer and the closest walls in this test case, which is very favourable to the sensor switching to LES mode. This observation motivated the introduction of the third test case described in the next section and meant as a more demanding test case. Concerning the proposed new ZDES mode 2 (EP), the Q criterion plotted in fig.18also suggests that the resolved turbulence is similar to the original ZDES mode 2 (2012). The new method is capable of resolving 3D turbulence after a short relaxation distance behind the trailing edge with no significant difference with the original ZDES mode 2 (2012) visible here. The eddy viscosity field features a slower destruction, which is not unexpected since the shielding is much stronger with ZDES mode 2 (EP) as illustrated in the flat plate boundary layer test case before. Comparing the eddy viscosity levels very near the trailing edge in ZDES mode 2 (2012, C 1 = 21) and in ZDES mode 2 (EP) reveals that the region of low values of the eddy viscosity (favourable to the development of LES content) is expanding faster in the vicinity of the trailing edge with ZDES mode 2 (EP) than with ZDES mode 2 (2012, C 1 = 21). This appears clearly if the levels appearing in blue are compared behind the trailing edge for ZDES mode 2 (EP), which is close to ZDES mode 2 (2012) and growing faster than ZDES mode 2 (2012, C 1 = 21) in the downstream direction. The explanation for this observation is that the modification in ZDES mode 2 (2012, C 1 = 21) increases the sensitivity of the sensor detecting RANS regions anywhere in the boundary layer. On the

Figure 18 :

 18 Figure 18: Isosurface of the Q criterion in the mixing layer, coloured by streamwise velocity. Insert: mean eddy viscosity field near the trailing edge.

Figure 19 :

 19 Figure 19: Streamwise evolution of the vorticity thickness in the mixing layer

  (mean velocity) and 21 (streamwise turbulent intensity), where ∆U = U high -U low is the shear velocity and η = (y -y 1/2 )/δ ω is a self-similarity parameter with y 1/2 the location where (U -U low )/(U high -U low ) = 1/2. The comparison of the mean velocity profiles confirms that the hybrid RANS/LES feature some delay in the initial stages of the development of the shear layer but do recover a proper profile further downstream (except for the original DDES (2006)). ZDES mode 2 (2012, C 1 = 21) seems to deviate from the

Figure 22 :

 22 Figure 22: Power Spectral Density of the streamwise velocity fluctuations in the mixing layer

4. 4 .

 4 Test case #3: Backward Facing Step 4.4.1. Test case

Figure 23 :

 23 Figure 23: A priori test of the shielding and inhibition functions on the backward facing step (test case #3). Colormap: difference f P -f d (r d )f P 2 , pink isoline: f d (r d ) = 0.2, black isoline: f P = 0.2 (functions as detailed in equations A.4, A.7, A.10).

Figure 24 :Figure 25 :

 2425 Figure 24: Instantaneous eddy viscosity field in the backward facing step flow.

Figure 26 : 4 . 5 .

 2645 Figure 26: Mean velocity and streamwise turbulent intensity profiles in the backward facing step flow.

  ∆z + = 100 (respectively ∆x + = 150) which correspond classically to WMLES grid resolution making the present grid very challenging to assess the robustness of the present ZDES mode 2 formulation. The time step is fixed at ∆t CF D = 2.10 -7 s, i.e. ∆ t = ∆t CF D (U 0 /c) = 2.1 • 10 -4 in terms of chord-based convective time scale, with six Newton inner iterations. The unsteady signals are acquired over almost 8 periods of the shock motion, which represents about 350 000 CPU hours,

Figure 27 :

 27 Figure 27: Isosurface of the Q criterion Q = 40 U 0 c

  The power spectral densities (PSD) of pressure fluctuations for two characteristic sensors are compared with the experiment in figure30(b). Due to the short duration of the ZDES simulation, an auto-regressive (AR) model based on Burg's method

5 Figure 30 :

 530 Figure 30: Salient features of the pressure field compared with experimental data [34].

  Finally, the new ZDES mode 2 should be, according to objective #4, fully compatible with the other modes (1 and 3) of ZDES, i.e. the subgrid scale model in LES regions is the same for all three modes, ensuring a continuity of LES treatment of resolved turbulence across zones set to different ZDES modes.The proposal to reach these four objectives is a new version of ZDES mode 2 with enhanced protection capabilities called ZDES mode 2 (EP). It involves three key ingredients. The first ingredient is a second shielding function intended for the detection of the regions of attached boundary layers not protected by the original shielding function. It relies on the wall-normal gradient of eddy viscosity, which is negative in the outer part of the wake layer where the second shielding function is needed since the original function fails to detect this region. The second ingredient is an inhibition function which is designed to switch off the second shielding in presence of flow separation. This avoids the RANS shielding of some detached flow regions and it relies on the wall-normal gradient of the vorticity magnitude as this quantity takes significant positive values in presence of flow separation. Finally, the third ingredient is a strong enhancement of the destruction term of eddy viscosity in so-called grey areas where the shielding and inhibition functions indicate that the model is in an intermediate state between RANS and LES. The rationale is to almost eliminate such grey areas by strongly destructing the eddy viscosity as soon as the RANS protection is becoming inactive in order to recover LES levels of subgrid scale viscosity as quickly as possible.The proposed new ZDES mode 2 (EP) has been calibrated in a set of flat plate boundary layers experiencing a wide variety of pressure gradients and Reynolds numbers. The model constants result from a compromise that should be suited for any flow configuration since special care has been taken to obtain non-dimensional Reynolds-number invariant sensors.The generality of the formulation and constants has been confirmed by a priori tests relying on RANS simulations of three-dimensional flow configurations. It is consequently not recommended to tune the constants in future use of the model since they must be case-independent, making the model general.The new ZDES mode 2 (EP) has then been tested by simulating four cases. The first test case shows that the considered flat-plate boundary layer under zero pressure gradient is perfectly shielded by ZDES mode 2 (EP) even under infinite mesh refinement as artificially simulated by setting C DES = 0. The second test case is a spatially developing mixing layer where the new ZDES mode 2 (EP) demonstrates a minimum delay in the formation of instabilities compared with the original ZDES mode 2 (2012). It is more efficient than the original ZDES mode 2 (2012) modified with C 1 = 21 for a better shielding, although this latter modification is likely to be insufficient in attached boundary layers with strong adverse pressure gradients. Indeed, ZDES mode 2 (EP) introduces a second sensor in order to shield attached boundary layers in a more general way, which goes far beyond the simple tuning of a constant of the original model, the latter option being strongly discouraged. The third test case demonstrates the capacity of the new ZDES mode 2 (EP) to resolve turbulence in a backward facing step with a performance very similar to the original ZDES mode 2 (2012) although this is a very challenging test case since the much stronger shielding for attached boundary layers must not be active in the free shear

  The development of a more robust shielding function is of outstanding importance for the safe use of hybrid RANS/LES methods on industrial configurations but no such function with a satisfying degree of universality is available in the open literature, as this short review has shown. In this paper, we propose to enhance the shielding properties of the original f d Performance of the shielding in ZDES mode 2 (EP): the new shielding should protect the whole attached boundary layer profile including on infinitely fine meshes with arbitrary pressure gradients.

	function in the frame of hybrid RANS/LES calculations in a case-independent manner, leading to a new version of ZDES
	mode 2 with Enhanced Protection of attached boundary layers, ZDES mode 2 (EP). Four major objectives are set to the
	study:
	1.

2. Additive construction of ZDES mode 2 (EP): the shielding should be at least as efficient as the original shielding used in ZDES mode 2 (2012) in all possible flow cases. This implies that the original shielding function will be kept and completed by a second function instead of being replaced by a different function. This choice is motivated by the wide range of successful applications of methods such as DDES or ZDES mode 2 based on the original f d function since 2006.

Table 1 :

 1 Classification of typical flow problems and associated modes of ZDES. I: separation fixed by the geometry, II: separation induced by a pressure gradient on a curved surface, III: separation strongly influenced by the dynamics of the incoming boundary layer. whereas f v1 , f v2 and f w are computed as in the original RANS model if d w ≤ C DES ∆I . Deactivating the near-wall functions in LES mode prevents an excessive destruction of subgrid scale viscosity on very fine meshes. In this strategy, the low-Reynolds number correction function Ψ(ν t /ν) [62] is consequently not needed.

  4, A.7, A.10) will be assessed in four test cases.The first test case is a spatially developing turbulent boundary layer, where the aim is to assess the capacity of the new model to properly shield the attached boundary layer under infinite grid refinement (this is deliberately a very demanding test). The mitigation of the delay in the development of resolved turbulence is evaluated in the next two test cases, namely a spatially developing mixing layer and a backward facing step flow. The latter case is very demanding because of the proximity between the free shear layer and the bottom wall which makes crucial the detection of flow separation.Four different methods are compared on these test cases: the original DDES (2006), the original ZDES mode 2 (2012)

[START_REF] Fukushima | Wall-modeled large-eddy simulation of transonic airfoil buffet at high reynolds number[END_REF] Indeed, the activation of the obtained sensor is more likely in presence of high levels of ν t , which may typically be the case in free shear layers still treated in RANS mode at the initialisation of the simulation. The sensor will also tend to be activated in regions of low velocity gradients, which may be encountered in separated flows. The whole velocity gradient rather than only its symmetric or antisymmetric parts is used in order to reduce the risk of undesired activation of the

Available at https://turbmodels.larc.nasa.gov/3dsqduct val.html [retrieved 31 July

2019].

Available at https://turbmodels.larc.nasa.gov/naca4412sep val.html [retrieved

July 2019].

Acknowledgments

The authors wish to thank all the people involved in the past and present evolution of the FLU3M code.

experiment slightly more than ZDES mode 2 (EP) in the early stages of the development (X/L = 0.2). This is confirmed by the peak in the streamwise turbulent intensity profiles at X/L = 0.2, where ZDES mode 2 (2012, C 1 = 21) is below ZDES mode 2 (EP), which itself is slightly below the original ZDES mode 2 (2012, C 1 = 8). This illustrates that a slight delay in the formation of instabilities is introduced by ZDES mode 2 (EP) but it remains small and it is less important than the delay introduced by setting C 1 = 21 in ZDES mode 2 (2012). After the full development of resolved turbulence in the streamwise direction, the turbulent intensity profiles of ZDES mode 2 (EP), ZDES mode 2 (2012, C 1 = 21) and the original ZDES mode 2 (2012, C 1 = 8) are virtually indiscernible at X/L = 0.8. In contrast, the turbulent intensity profiles obtained with the original DDES (2006) are severely underestimated at X/L = 0.2 and overestimated at X/L = 0.8 (because of the excessive coherence and length scale of resolved turbulence at the latter station).

In order to check that the spectral content of the resolved turbulence is physically correct, spectra of the streamwise velocity fluctuations are presented in fig. 22. These confirm that the resolved turbulence in the original DDES (2006) has no physical background. On the contrary, the spectra obtained with all versions of the ZDES mode 2 are close to each other and in fair agreement with the experimental data, especially at X/L = 0.8, except of course in the highest frequency range which is beyond the LES cutoff. A significant inertial spectral region is even resolved by ZDES at X/L = 0.8. This layer in spite of the proximity of the walls. The results gathered for these three cases confirm the relevance of ZDES mode 2 (EP) with respect to the four objectives set for the new method. Finally, the ability of the model to deal with a flow separation whose location and amplitude is changing over time has been demonstrated by simulating the transonic buffet over the supercritical OAT15A airfoil. The results on this very demanding test case are extremely encouraging and confirm that the model can reconstruct a RANS treatment of a boundary layer which was mostly treated in LES during the preceding flow separation.

In view of the presented results, the use of the new ZDES mode 2 (EP) with the constants given in this paper is recommended for any case where an automatic hybrid RANS/LES treatment of attached and massively separated flows is needed. The authors strongly discourage the tuning of constants, including in the original models, since this lacks generality and leads to case-dependent treatments. In contrast, the new robust ZDES mode 2 (EP) is a case-independent answer to the demand for a general automatic and robust RANS/LES treatment of attached and massively separated flows.

One final note is that in the zonal context of the ZDES technique, it is still possible to manually set a zone of attached flow to the RANS mode (ZDES mode 0), which is the safest shielding possible. This shows the interest of keeping a zonal framework within which ZDES mode 2 (EP) itself is an automatic mode. This is crucial regarding mild flow separations.

Such situations will be treated mostly in LES with ZDES mode 2 (EP). For a more accurate resolution of these flows, defining a ZDES mode 3 (WMLES) zone in the concerned region may sometimes be desirable. If alternatively one wishes a RANS treatment of the mild flow separation, this may be sought by deactivating the inhibition function (at the cost of losing resolved LES content especially in flow situations close to walls) or ideally obtained by manually defining a ZDES mode 0 (RANS) zone around the mild flow separation location if the geometry is not too complex. This also holds for the interaction of a resolved shear layer with the wake layer of an attached boundary layer, which is a very ambiguous situation where the user might desire a RANS or LES behaviour depending on the focus of the study. Indeed, the distinction between attached flows, mild and massive flow separations is probably too simplistic to represent the diversity of flow situations treated by hybrid RANS/LES approaches, meaning that an automatic approach benefits a lot from being informed by the user, which is what a zonal approach does.

Appendix A. Practical formulation and implementation of ZDES mode 2 (EP)

For the application of the proposed second shielding function and inhibition function to ZDES mode 2, it is most convenient to use the pseudo eddy-viscosity ν as the estimate of the eddy viscosity in the sensors ( 25) and [START_REF] Fukushima | Wall-modeled large-eddy simulation of transonic airfoil buffet at high reynolds number[END_REF]. Indeed, ν is the transported variable directly available from the Spalart-Allmaras turbulence model, and it provides a more robust shielding than the eddy viscosity ν t which may decay very quickly (because of the f v1 function) in the regions of low eddy viscosity. It should be noted that whenever ν/ν is significantly greater than 1, there is practically no difference between ν and ν t and the sensors will be used as intended in their study detailed in §3.1-3.2. The resulting formulation of ZDES mode 2 (EP) based on the Spalart-Allmaras turbulence model [START_REF] Spalart | A one equation turbulence model for aerodynamic flows[END_REF] is detailed in the following.

The eddy viscosity field is related to the transported pseudo eddy-viscosity by the classical relation:

where the transport equation reads as usual:

The model length scale is defined by:

where d w is the wall distance and the protection function reads as:

The sensor of the first (original) shielding function is:

to which the function f d is applied:

The second shielding function may be defined as:

where the sensor reads:

Alternatively, a stronger second shielding function may be obtained as follows:

with β = 2.5 in order to procure a total shielding (in the sense f P = 0 across the whole boundary layer profile) in strong adverse pressure gradient cases close to separation such as illustrated in fig. 8. Such a reinforcement of f P 2 has not been retained for the simulations presented in section 4 in order to isolate the effect of the new sensor G ν . It may be needed in flows gradually approaching separation. Its impact on the formation of LES content in the mixing layer and backward facing step cases (test cases #2 and #3) is smaller than the delay caused by setting C 1 = 21 in the original ZDES mode 2 (2012).

The inhibition function is formulated as follows:

and C 4 = 0.03

with the following sensor:

where |ω| is the vorticity magnitude.

If for numerical reasons the vorticity field is noisy, limiting the activation of the inhibition function may be needed in order to avoid undesired switching to LES mode caused by numerical fluctuations of vorticity. Such a limited inhibition function may be obtained by replacing f R with the following expression:

By definition of the vorticity magnitude, if ζ > √ 2, then the limited inhibition function is always f R,lim = 1 (this may be used to seek a very robust RANS shielding of attached boundary layers including in presence of resolved turbulence from a shear layer located above the boundary layer, or to seek to avoid LES switching in a mild flow separation, but this would be at the cost of losing LES resolved content in some flow situations close to walls even though an acceptable LES resolution of free shear layers far away from walls should still be obtained). The recommended value is ζ = 0.8, which ensures that numerical fluctuations of vorticity are ignored as soon as they are small compared with the norm of the whole velocity gradient, while the original inhibition function, which is physically correct, is not altered in separated flows where the sheared regions are characterised by values of |ω| and √ u i,j u i,j of the same order. Limiting the inhibition function by setting ζ = 0.8 has no significant impact on the mixing layer and backward facing step cases (test cases #2 and #3) and it may enhance the numerical robustness of simulations where the convergence of the vorticity field is not perfect. It should be noted that the limited f R,lim has not been used for the simulations presented in section 4 because it was not needed in these cases. The limitation of f R is purely motivated by numerical contingencies and should not alter the physical model which is satisfying in its non-limited form.

Two boolean tests are introduced for the switch to LES and for the enhancement of destruction in the grey areas respectively:

The relation between the pseudo eddy-viscosity field and the eddy viscosity is completed by the definition of the function involved:

whereas the cell size estimate is the following:

As for the production term, it involves the following expression:

with

Finally, the destruction term function is:

where the functions from mode 1 are the following:

Concerning the evaluation of the wall normal derivative ∂/∂n involved in the new sensors G ν and G Ω , it is obtained as the scalar product between the gradient of the considered field and the local wall normal n. This latter vector is classically obtained from the wall distance field d w :

Similarly to the SA turbulence model, this implies that the hybrid model includes a dependency to the nearest wall even in free shear flows far away from walls, which is not an issue. Indeed, both shielding sensors r d (A.5) and G ν (A.8) decay with the wall distance. This especially implies that the second shielding function will tend not to switch to RANS far away from walls. Consequently, the value of the inhibition function f R , which can vary because the sensor G Ω (A.11) does not decay with d w , tends to have no effect in such regions (see eq. (A.4) with f P 2 = 1). As for the numerical robustness of the method, the discontinuity of the local wall normal n field along the separation surfaces between the nearest walls has not been found to cause any numerical difficulty during the simulations. If however more regularity were wished in some cases on physical grounds, smoothing the wall normal field would still be an option since the n field may be seen as an input to the present model.

It is recommended to initialize the pseudo eddy-viscosity field from a RANS field since the instantaneous values of this field will be involved in the RANS shielding of attached boundary layers. This is not an issue since ZDES computations √ u i,j u i,j /(u τ /δ) κd w /δ (B.1) Indeed, based on typical behaviours in Zero Pressure Gradient canonical boundary layers in the outer layer, the maximum eddy viscosity is reached near d w = 0.5δ and is of the order of u τ κ 0.5δ/2 (assuming a parabolic profile in the outer layer with an initial slope matching the logarithmic law as an approximate guess), so that a fair approximation of the order of magnitude of -∂ν t /∂n in the outer part of the wake layer is (u τ κ 0.5δ/2)/(δ -0.5δ) = u τ κ/2. As for the velocity gradient, the difference between U + ∞ and u + (d w = 0.5δ) is of the order of 3 (see typical wake laws), corresponding to a distance δ -0.5δ = 0.5δ. The typical value of G νt in the outer part of the wake layer (where the average d w is 0.75δ) is consequently of the order of C 3 • (u τ κ/2)/((3u τ /(0.5δ))κ(0.75δ)) = 1/9 C 3 . This suggests that choosing C 3 of the order of 9 could be sufficient for ZPG boundary layers. This is however a very approximative evaluation and the pressure gradient cases must also be considered. This is why the calibration of C 3 proposed in section 3.2 includes non-zero pressure gradients and a larger value of C 3 has been retained for a general use as a case-independent constant.