domain: spi phys

text:
Francesca Basile
email: francesca.basile@onera.fr

Jean-Baptiste Chapelier
email: jean-baptiste.chapelier@onera.fr

Romain Laraufie
email: romain.laraufie@airbus.com

Pascal Frey
email: pascal.frey@sorbonne-universite.fr

hp-adaptive hybrid RANS/LES simulations for unstructured meshes with the discontinuous Galerkin method

In this paper, an hp-adaptation strategy designed for discontinuous Galerkin methods is applied to hybrid RANS/LES simulations. The 3D hp-adaptive strategy is suited for tetrahedral and hybrid prismatic/tetrahedral meshes, and relies on a metric-based remeshing approach. The metric field and the polynomial map of the adapted meshes are built from an a posteriori error estimator which couples the measure of the energy associated with the highest-order modes and the inter-element jumps, combined with a smoothness sensor which guides the choice between hand p-adaptation. The turbulence modelling relies on a Zonal Detached Eddy Simulation approach (ZDES mode 1). The resolution requirements in terms of minimum wall spacing Δ𝑦 + are first assessed from steady RANS flat plate simulations at Reynolds 5 • 10 6 and the value of the DES constant is assessed for DG from 3D Taylor-Green vortex flow simulations at Reynolds 5000. The developed hp-adaptation algorithm is then applied to hybrid RANS/LES simulations of the PPRIME nozzle configuration at diameter-based Reynolds 10 6 , starting the adaptive process from previously RANS-adapted meshes.

I. Introduction

In the design of engine nozzles of the future, scale-resolving simulations are necessary in order to improve the understanding of turbulence and the noise generation mechanisms. Despite their extensive and well-assessed use for industrial purposes, Reynolds-averaged Navier-Stokes (RANS) simulations fail in predicting accurately non-equilibrium turbulent flows in which the turbulent large scales strongly affect the mean flow quantities [START_REF] Chaouat | The state of the art of hybrid RANS/LES modeling for the simulation of turbulent flows[END_REF], as for example with strongly detached and transitional flows [START_REF] Sagaut | Multiscale & multiresolution approaches in turbulence[END_REF]. In fact in RANS, one solves the fluid dynamics equations only for the averaged quantities, while the effects of all the scales of turbulence are modelled. Large Eddy Simulation (LES) is designed to capture the unsteady features in transitional flows, gas turbine combustors, nozzles. Nonetheless, LES yields prohibitive costs to capture the wall-bounded turbulent dynamics at high Reynolds numbers. Hybrid RANS/LES models aim at combining the most attractive properties of the two methods: the attached boundary layer is predicted by RANS models, while the unsteady 3D large eddies are captured by LES, decreasing the computational cost of separated turbulent flows with respect to LES. This allows for a strong reduction of the degrees of freedom (dofs) which wall-resolved LES would need to capture the smaller structures developed in the boundary layer.

The original Detached Eddy Simulation (DES) technique [START_REF] Spalart | Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[END_REF] is a hybrid RANS/LES non zonal approach which, without requiring the user to manually select zones to be treated with RANS or LES, uses the RANS method in the boundary layer and switches continuously to LES far from the walls. If the mesh is an "ambiguous mesh" not compliant with the DES requirements, the switch between the two models can be located in the boundary layer, resulting in an early switch to LES. As a consequence, the velocity fluctuations from the LES field do not balance the loss of the modelled Reynolds stresses. Spalart et al. [START_REF] Spalart | A new version of detached-eddy simulation, resistant to ambiguous grid densities[END_REF] presented an improvement of the original DES, called Delayed DES (DDES), which detects the boundary layer with a sensor including information on gradients and viscosity, in order to solve it with the RANS equations.

An alternative approach is the Zonal Detached Eddy Simulation (ZDES), which is based on a zonal formulation of the hybrid length scale that allows for combining the zonal approach with DES/DDES [START_REF] Deck | Recent improvements in the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF]. Within ZDES, the user (or the model itself for ZDES mode 2) selects both RANS and DES domains and the hybrid length scale corresponding to the problem of interest (flows where the separation is fixed by the geometry, or induced by a pressure gradient on a curved surface, or strongly influenced by the dynamics of the incoming boundary layer). In particular, the ZDES formulation used in this work is the so-called mode 1, proposed by Deck [START_REF] Deck | Zonal-detached-eddy simulation of the flow around a high-lift configuration[END_REF] in the framework of flows with separation fixed by the geometry, for which in RANS regions the model is forced to behave as a RANS model, while in DES regions, the model can switch from the RANS mode to the LES mode. This approach has already been applied to a wide range of cases [START_REF] Deck | Zonal-detached-eddy simulation of the flow around a high-lift configuration[END_REF][START_REF] Brunet | Zonal-detached eddy simulation of transonic buffet on a civil aircraft type configuration[END_REF][START_REF] Brunet | Zonal-detached eddy simulation of a civil aircraft engine jet configuration[END_REF], including jet simulations [START_REF] Gand | Zonal Detached Eddy Simulation (ZDES) using turbulent inflow and high order schemes: Application to jet flows[END_REF][START_REF] Chauvet | Zonal Detached Eddy Simulation of a controlled propulsive jet[END_REF].

The simulation of turbulent unsteady phenomena requires as well accurate enough numerical schemes, presenting low dissipation and dispersion properties, for which discontinuous Galerkin (DG) methods [START_REF] Karniadakis | Spectral/hp element methods for computational fluid dynamics[END_REF][START_REF] Dolejší | Discontinuous galerkin method[END_REF] are good candidates. These methods are based on the variational formulation of the governing equations and combine features of Finite Volume (FV) and Finite Element (FE) methods. In recent years, these methods have become very popular for the solution of nonlinear convection dominated flow problems [START_REF] Cockburn | Discontinuous Galerkin methods for convection-dominated problems[END_REF]. Many aspects make DG methods particularly attractive, such as their high-order of accuracy achieved on arbitrary unstructured meshes, accurate description of curved boundaries and suitability to parallel computing thanks to a compact stencil. Another advantage of DG methods is the possibility to exploit a posteriori jump and spectral error indicators for mesh adaptation. These error indicators are convenient thanks to their efficiency, locality, simplicity and low computational cost [START_REF] Mavriplis | A posteriori error estimators for adaptive spectral element techniques[END_REF][START_REF] Persson | Sub-cell shock capturing for discontinuous Galerkin methods[END_REF][START_REF] Gassner | Explicit Discontinuous Galerkin Schemes with Adaptation in Space and Time[END_REF][START_REF] Remacle | An adaptive discontinuous Galerkin technique with an orthogonal basis applied to compressible flow problems[END_REF].

In particular the growing interest around DG lays in their hp-adaptivity, which means that one can not only locally adapt the size h of the mesh, but also the degree of the polynomials, p, within the element. In practice hp-adaptive methods lead to the concentration of the degrees of freedom (dofs) in regions of interest of the flow, in order to optimize both the computational cost and the accuracy of the simulations [START_REF] Dolejší | hp-DGFEM for nonlinear convection-diffusion problems[END_REF][START_REF] Leicht | Error estimation and hp-adaptive mesh refinement for discontinuous Galerkin methods[END_REF][START_REF] Wang | Adjoint-based h-p adaptive discontinuous Galerkin methods for the 2D compressible Euler equations[END_REF].

As regards the h-adaptive strategy for our algorithm, a metric-based remeshing approach [START_REF] Loseille | Continuous mesh framework part I: well-posed continuous interpolation error[END_REF][START_REF] Frey | Anisotropic mesh adaptation for CFD computations[END_REF] has been chosen for this work, showing promising results in the literature thanks to its flexibility to prescribe a precise size (and anisotropic features as aspect ratio adaptatiand orientation in an anisotropic context) to each element of the mesh. The employed remeshing tool is MMG [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF], which supports both 2D and 3D remeshing (MMG2D and MMG3D) and relies on primitive operators like edge splitting, edge collapsing, edge swapping, and node relocation. The MMG library has been already used extensively in the context of FV and FE methods [START_REF] Daviller | A mesh adaptation strategy to predict pressure losses in LES of swirled flows[END_REF][START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF][START_REF] Bourasseau | Contribution to a mesh refinement method based on the adjoint vector for the computation of aerodynamic outputs[END_REF][START_REF] Sadaka | Parallel finite-element codes for the simulation of two-dimensional and three-dimensional solid-liquid phase-change systems with natural convection[END_REF].

Regarding hp-adaptivity, the local error estimator indicates the elements that should be refined, but does not inform whether to refine the element by h or p. A method for making that choice is called an hp-decision strategy. Generally this choice is made according to an estimate of the solution smoothness in an element: if it is sufficiently smooth the adaptive algorithm opts for p-enrichment, while non-smooth zones are h-refined. Some different strategies have been proposed concerning the hp-decision [START_REF] Mavriplis | A posteriori error estimators for adaptive spectral element techniques[END_REF][START_REF] Leicht | Error estimation and hp-adaptive mesh refinement for discontinuous Galerkin methods[END_REF][START_REF] Burgess | An hp-adaptive discontinuous Galerkin solver for aerodynamic flows on mixed-element meshes[END_REF][START_REF] Wang | Discretization Error Estimation for Discontinuous Galerkin Methods using Error Transport Equations[END_REF][START_REF] Chalmers | A parallel hp-adaptive high order discontinuous Galerkin method for the incompressible Navier-Stokes equations[END_REF]. In particular Mavriplis [START_REF] Mavriplis | A posteriori error estimators for adaptive spectral element techniques[END_REF] determined if the solution is locally smooth or non-smooth by computing the decay rate of the Legendre expansion coefficients of the solution under the assumption that for non-smooth solutions, the discontinuities in the solutions deteriorate this decay rate. In this work we employ a smoothness indicator similar to the one proposed by Mavriplis and used in our previous work on steady flows [START_REF] Basile | Unstructured h-and hp-adaptive strategies for discontinuous Galerkin methods based on a posteriori error estimation for compressible flows[END_REF]. We employ here a slightly different formulation of the power decay of the modes, which is found more appropriate in the context of scale-resolving turbulent flows. The threshold value between smooth and non-smooth behaviour is chosen as the theoretical value of the slope of the turbulent spectrum -5/3 [START_REF] Chapelier | A spectral-element dynamic model for the Large-Eddy simulation of turbulent flows[END_REF].

In our work we apply an hp-adaptation algorithm, to a turbulent flow configuration relevant in an industrial context, namely a subsonic isothermal jet issued from a round nozzle at 𝑀 𝑗 = 0.9 and 𝑅𝑒 𝐷 = 10 6 . A hybrid RANS-LES framework based on the ZDES mode 1 is employed to model the turbulence in the jet region. An isotropic mesh refinement strategy is chosen, given the fundamental isotropic character of free-shear, small-scale turbulence [START_REF] Pope | Turbulent flows[END_REF] which characterizes the analyzed flow configuration. As a matter of fact, it can be expected that the ideal fully adapted mesh would be isotropic for this kind of physical phenomenon. In the literature we can find examples of isotropic remeshing strategies, mostly for unsteady complex applications [START_REF] Daviller | A mesh adaptation strategy to predict pressure losses in LES of swirled flows[END_REF][START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF].

The anisotropy of the wall-attached flow inside the nozzle is treated with a fixed pseudo-structured layer of prisms, easily adaptable with a pure 𝑝-adaptive approach in this area. In fact the use of prismatic or hexahedral elements for the boundary layer region is beneficial for computing the gradients with high-accuracy, in particular for RANS and hybrid RANS/LES simulations, where the near wall resolution constraint 𝑦 + ≃ 1 is fundamental to obtain a relevant solution. To this end it is common that unstructured meshes involve a structured (extruded regular quadrilaterals) or pseudo-structured (extruded unstructured quadrilaterals or triangles) boundary layer around the body, which blends with the outer tetrahedral field regions [START_REF] Marcum | Unstructured mesh generation using advancing layers and metric-based transition for viscous flowfields[END_REF][START_REF] Antoniadis | Assessment of high-order finite volume methods on unstructured meshes for RANS solutions of aeronautical configurations[END_REF][START_REF] Tomac | Steps towards automated robust RANS meshing[END_REF]. In particular, in the field of nozzle/jet configurations an unstructured mesh approach has been extensively used in the literature, both for RANS simulations [START_REF] Dippold | Generating a Grid for Unstructured RANS Simulations of Jet Flows[END_REF][START_REF]Acoustic Reference Nozzle with Mach 0.97, Unheated Jet Flow[END_REF] and LES of resolved jets

A. The Reynolds-averaged Navier-Stokes equations with the Spalart-Allmaras turbulence model

The motion of a compressible fluid in a three-dimensional domain Ω ∈ R 3 is described by the compressible Navier-Stokes (NS) equations with source terms, which take the form

𝜕 𝑡 u + ∇ • (F 𝑐 (u) -F 𝑣 (u, ∇u)) = s(u, ∇u), ∀x ∈ Ω, 𝑡 > 0 u(x) = u 0 , ∀x ∈ Ω (1)
with appropriate boundary conditions prescribed on 𝜕Ω. The vector u represents the conservative variables, with 𝜌 being the density, v being the velocity vector and 𝐸 =

𝑝 (𝛾 -1) 𝜌 + v • v 2 the specific total energy. 𝑝 is the static pressure, 𝛾 = 𝐶 𝑝 𝐶 𝑣
is the ratio of specific heats. F 𝑐 (u) and F 𝑣 (u, ∇u) are the convective and diffusive fluxes defined, respectively, as:

u =        𝜌 𝜌v 𝑇 𝜌𝐸        , F 𝑐 (u) =        𝜌v 𝑇 𝜌vv 𝑇 + pI 𝜌𝐸v 𝑇 + pv 𝑇        and F 𝑣 (u, ∇u) =        0 𝝉 v 𝑇 𝝉 -q 𝑇        (2)
with:

𝝉 = 2𝜇S 𝐷 = 𝜇 ∇v + (∇v) 𝑇 - 2 3 (∇ • v)I , q = -𝜆∇𝑇 (3)
p = 𝜌𝑅𝑇 ( 4
)
where 𝜇 is the dynamic viscosity, S 𝐷 is the deviatoric component of the strain-rate tensor S = 1 2 (∇v + (∇v) 𝑇 ), 𝑇 is the temperature, 𝑅 is the specific gas constant, 𝜆 = 𝜇 𝐶 𝑝 𝑃𝑟 is the thermal conductivity, with 𝑃𝑟 the Prandtl number. In this work, the value 𝑃𝑟 = 0.72 is used.

In the RANS formulation the mean flow equations are coupled with the one-equation turbulence model of Spalart-Allmaras (SA) [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] with the modification described in [START_REF] Allmaras | Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model[END_REF].

The turbulent conservative variable 𝜌 ν is added to the vector of conservative equations and the conservative and diffusive fluxes F 𝑐 (𝜌 ν), F 𝑣 (𝜌 ν, ∇𝜌 ν) are defined as:

F 𝑐 (𝜌 ν) = 𝜌 νv 𝑇 , F 𝑣 (𝜌 ν) = 1 𝜎 (𝜇 + 𝑓 𝑛1 𝜌 ν)∇ ν𝑇 (5)
In the diffusive fluxes, the turbulent stress tensor 𝜏 𝜏 𝜏 t and the turbulent heat fluxes q t are added respectively to 𝜏 𝜏 𝜏 and q defined in equations (3):

𝜏 𝜏 𝜏 t = 2𝜇 𝑡 S 𝐷 , q t = - 𝜇 𝑡 𝑃𝑟 𝑡 𝐶 𝑝 ∇𝑇 (6)
where 𝑃𝑟 𝑡 = 0.9 is the turbulent Prandtl number, 𝜇 𝑡 is the turbulent dynamic viscosity:

𝜇 𝑡 = 𝜌 ν 𝑓 𝑣1 ( 𝜒) for ν ≥ 0 0 for ν < 0 , 𝑓 𝑣1 ( 𝜒) = 𝜒 3 𝜒 3 + 𝑐 3 𝑣1 , 𝜒 = 𝜌 ν 𝜇 (7)
The source terms act on the conservation equation only for the turbulent variable 𝜌 ν and read

S(𝜌 ν, ∇(𝜌 ν)) = -𝜌(𝑃 -𝐷) - 𝑐 𝑏2 𝜎 𝜌∇ ν + 1 𝜎 (𝜈 + 𝑓 𝑛1 ν) ∇𝜌 • ∇ ν (8)
where the production and destruction terms 𝑃 and 𝐷 are defined by:

𝑃 = 𝑐 𝑏1 (1 -𝑓 𝑡2 ) ω ν for ν ≥ 0 𝑐 𝑏1 (1 -𝑐 𝑡3 )𝜔 ν for ν < 0 𝐷 =            𝑐 𝑤1 𝑓 𝑤 - 𝑐 𝑏1 𝜅 2 𝑓 𝑡2 ν 𝑑 2 for ν ≥ 0 -𝑐 𝑤1 ν 𝑑 2 for ν < 0 (9)
and

𝑓 𝑛1 =        1 for ν ≥ 0 𝑐 𝑛1 + 𝜒 3 𝑐 𝑛1 -𝜒 3 for ν < 0 , 𝑓 𝑡2 = 𝑐 𝑡3 exp(-𝑐 𝑡4 𝜒 2 ), (10)
with

𝑔 = 𝑟 + 𝑐 𝑤2 (𝑟 6 -𝑟), 𝑟 = min 𝑟 𝑙𝑖𝑚 , ν 𝜔𝜅 2 𝑑 2 , 𝑓 𝑤 = 𝑔 1 + 𝑐 6 𝑤3 𝑔 6 + 𝑐 6 𝑤3 1 6 (11)
𝑑, the distance to the nearest wall and 𝜔 the vorticity magnitude. The modified vorticity magnitude ω is given by

ω =        𝜔 + 𝜔 for 𝜔 > -𝑐 𝑣2 𝜔 𝜔 + 𝜔(𝑐 2 𝑣2 𝜔 + 𝑐 𝑣3 𝜔) (𝑐 𝑣3 -2𝑐 𝑣2 )𝜔 -𝜔 for 𝜔 < -𝑐 𝑣2 𝜔 (12)
where Ω and 𝑓 𝑣2 are given by

𝜔 = ν 𝜅 2 𝑑 2 𝑓 𝑣2 , 𝑓 𝑣2 = 1 - 𝜒 1 -𝑓 𝑣1 (13)
The values of the constants can be found in [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF][START_REF] Allmaras | Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model[END_REF].

B. The Zonal Detached Eddy Simulation Mode 1

The Zonal Detached Eddy Simulation (ZDES) which is used in this work is a hybrid RANS/LES technique, formalized by Deck [START_REF] Deck | Zonal-detached-eddy simulation of the flow around a high-lift configuration[END_REF] and initially based on the SA-RANS model, as an efficient solution to achieve a faster decay of the eddy viscosity in the LES mode with respect to classical DES/DDES.

In the ZDES framework, the wall distance 𝑑 𝑤 of the SA turbulent variable transport equation, is replaced with the hybrid length-scale d𝑍𝐷𝐸𝑆 . In practice we apply this modification in the destruction term in equation [START_REF] Gand | Zonal Detached Eddy Simulation (ZDES) using turbulent inflow and high order schemes: Application to jet flows[END_REF]. The eddy viscosity in the LES region scales then with the length-scale and the vorticity magnitude 𝜌 ν ∼ ωΔ 𝑣𝑜𝑙 . The first version of ZDES chosen in this work, is generally employed for flow configurations where the position of the separation is known a priori from the geometry, namely for massively separated flows.

In particular a user-defined zonal decomposition of the computational domain in RANS and DES areas can force RANS areas to act in RANS mode [START_REF] Deck | Numerical simulation of transonic buffet over a supercritical airfoil[END_REF]. This prevents the so-called "grid-induced-separation" phenomenon [START_REF] Spalart | A new version of detached-eddy simulation, resistant to ambiguous grid densities[END_REF] which can appear when a mesh is too fine in the longitudinal direction, where the original DES approach would reduce the RANS viscosity.

In the current work, where the studied test case is a jet issued from a nozzle, the model will be forced to RANS mode in the domain internal to nozzle, as depicted in figure 1.

Fig. 1 Forcing of the interior of the nozzle to RANS mode.

The hybrid length entering the SA equations is defined as: d𝑍𝐷𝐸𝑆 = 𝑑 𝑤 in RANS zones min(𝑑 𝑤 , 𝐶 𝐷𝐸𝑆 Δ𝐷𝐸𝑆 ) in DES zones [START_REF] Mavriplis | A posteriori error estimators for adaptive spectral element techniques[END_REF] where 𝐶 𝐷𝐸𝑆 is a coefficient usually set to 0.65 which has been calibrated in the decaying homogeneous turbulence [START_REF] Shur | Detached-eddy simulation of an airfoil at high angle of attack[END_REF]. Two aspects differentiate this formulation from the classical DES [START_REF] Spalart | Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[END_REF]. The first is that the sub-grid length scale Δ𝐷𝐸𝑆 is defined as the cubic root of the volume of the cell: Δ𝐷𝐸𝑆 = Δ 𝑣𝑜𝑙 = 3 √ 𝑉 𝑜𝑙, and not as the maximum grid extension Δ𝐷𝐸𝑆 = max(Δ 𝑥 , Δ 𝑦 , Δ 𝑧 ).

The second is that the near-wall corrections 𝑓 𝑣1 , 𝑓 𝑣2 , 𝑓 𝑤 are the same as the RANS Spalart-Allmaras model in RANS mode, while in LES mode they are modified:

𝑓 𝑣1 = 1, 𝑓 𝑣2 = 0, 𝑓 𝑤 = 1 ( 15
)
which are their asymptotic values far from the wall, avoiding a drop of subgrid viscosity that could be caused by the damping functions of the RANS model in resolved LES zones presenting low eddy viscosity levels.

C. DG Discretization

The DG discretization used in this work is based on a modal approach that relies on the use of a hierarchical and orthogonal polynomial basis for the Galerkin projection. In this work unsteady solutions are obtained by employing an explicit three-stages Runge-Kutta scheme. The DG method implemented in the CODA solver is briefly outlined below.

We start by defining a shape-regular partition of the domain Ω, into a tessellation T 𝐾 of 𝑁 non-overlapping and non-empty simplicial elements 𝐾 of characteristic size ℎ. We also define the sets E 𝑖 and E 𝑏 of interior and boundary faces in

T 𝐾 , such that E ℎ = E 𝑖 ∪ E 𝑏 . Let V 𝑝 ℎ = {𝜙 ℎ ∈ 𝐿 2 (Ω) : 𝜙 | 𝐾 ∈ P 𝑝 (𝐾)
, ∀𝐾 ∈ T 𝐾 } be the functional space formed by piece-wise polynomials of degree at most 𝑝, and (𝜙

u ℎ (x, 𝑡) = 𝑁 𝑝 ∑ 𝑙=1 𝜙 𝑙 𝐾 (x)u 𝑙 𝐾 (𝑡), ∀x ∈ 𝐾, 𝐾 ∈ T 𝐾 , ∀𝑡 ≥ 0 (16)
The polynomial coefficients (u 𝑙 𝐾 ) 1≤𝑙 ≤ 𝑁 𝑝 represent the degrees of freedom of the discrete problem in element 𝐾. The shape functions are polynomials that can be chosen arbitrarily. A methodology developed by Bassi et al. [START_REF] Bassi | On the flexibility of agglomeration based physical space discontinuous Galerkin discretizations[END_REF] consists in defining a starting set of monomial basis functions in each (arbitrarily shaped) element and applying a modified Gram-Schmidt orthonormalization procedure. The resulting basis yields a diagonal mass matrix in each element of the discretization, simplifying the resolution of the variational formulation. We indicate with Φ 𝐾 = {𝜙 1 𝐾 , ..., 𝜙 𝑁 𝑝 𝐾 } the generic basis of the polynomial space P 𝑝 (𝐾). Φ 𝑝 𝐾 is a hierarchical basis if it is contained in higher-order basis, that is Φ 𝑝 𝐾 ⊂ Φ 𝑝+1 𝐾 . We can rewrite equation ( 16) for an element 𝐾 of degree 𝑝 𝐾 as:

u ℎ (x, 𝑡) = 𝑝 𝐾 ∑ 𝑞=0 ∑ 𝑙 ∈𝑑 𝑞 𝜙 𝑙 𝐾 (x)u 𝑙 𝐾 (𝑡), ∀x ∈ 𝐾, 𝐾 ∈ T 𝐾 , ∀𝑡 ≥ 0 ( 17
)
where

𝑑 0 = {1} and 𝑑 𝑞 = 𝑙 ∈ 2...𝑁 𝑝 (𝐾) | 𝜙 𝑙 𝐾 ∈ P 𝑞 𝐾 \ P 𝑞-1 𝐾
is the set of indices of the basis functions of degree 𝑞. The conservation law is discretized in physical space by using a discontinuous Galerkin method and the semi-discrete variational form of the system of equations (1) thus reads: find

u ℎ in V 𝑝 ℎ such that ∀𝜙 ℎ ∈ V 𝑝 ℎ we have ∫ T 𝐾 𝜙 ℎ 𝜕 𝑡 u ℎ 𝑑𝑉 + L 𝑐 (u ℎ , 𝜙 ℎ ) + L 𝑣 (u ℎ , 𝜙 ℎ ) = 0 (18)
In equation ( 18) L 𝑐 and L 𝑣 represent the weak form of the convective and viscous terms respectively. The following notations are introduced: for a given interface 𝑒 in E 𝑖 we define the average operator as {{u}} = (u + + u -)/2, the jump operator is defined as [[u]] = u + ⊗ nu -⊗ n where u + and u -are the traces of the variable u at the interface between elements 𝐾 + and 𝐾 -. The DG discretization of the convective terms then reads

L 𝑐 (u ℎ , 𝜙 ℎ ) = - ∫ T 𝐾 F 𝑐 (u ℎ ) • ∇ ℎ 𝜙 ℎ 𝑑𝑉 + ∫ E 𝑖 [[𝜙 ℎ ]]h 𝑐 (u + ℎ , u - ℎ , n)𝑑𝑆 + ∫ E 𝑏 𝜙 + ℎ F 𝑐 (u 𝑏 ) • n𝑑𝑆 (19)
where the boundary values u 𝑏 = u 𝑏 (u + ℎ , u 𝑒𝑥𝑡 , n), with u 𝑒𝑥𝑡 a reference external state computed such that the boundary conditions are satisfied on E 𝑏 . The numerical flux h 𝑐 is chosen such that it is consistent and conservative. In this work we use the Roe flux [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF] with an entropy fix similar to that of Harten [START_REF] Harten | Self adjusting grid methods for one-dimensional hyperbolic conservation laws[END_REF] for all simulations.

The discretization of the viscous terms in this work is performed using the BR1 [START_REF] Bassi | A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations[END_REF] approach of Bassi and Rebay. This approach relies on the definition of the conservative variables gradients as auxiliary variables 𝝈 = ∇u which verify the following equations:

𝝈 -∇u = 0 ( 20
)
𝜕 𝑡 u + ∇ • F 𝑐 (u) + ∇ • F 𝑣 (u, 𝝈) = 0 ( 21
)
This leads to the introduction of the so called global lifting operator L ℎ such that:

𝝈 ℎ = ∇ ℎ u ℎ + L ℎ ( 22
)
and L ℎ satisfies the following condition

∫ T 𝐾 𝜙 ℎ L ℎ 𝑑𝑉 = - ∫ E 𝑖 {{𝜙}}[[u 𝑢 ]]𝑑𝑆 - ∫ E 𝑏 𝜙 + 2 (u + ℎ -u 𝑏 ) ⊗ n𝑑𝑆 (23)
The discrete variational form of the viscous term for the BR1 method therefore takes the form

L 𝑣 (u ℎ , 𝜙 ℎ ) = - ∫ T 𝐾 F 𝑣 (u ℎ , ∇ ℎ u ℎ + L ℎ ) • ∇ ℎ 𝜙 ℎ 𝑑𝑉 - ∫ E 𝑖 [[𝜙 ℎ ]]{{F 𝑣 (u ℎ , ∇ ℎ u ℎ + L ℎ )}} • n𝑑𝑆 - ∫ E 𝑏 𝜙 + ℎ F 𝑣 (u 𝑏 , ∇u 𝑏 + L ℎ ) • n𝑑𝑆 (24)

III. hp-adaptive algorithm

The hp-adaptive strategy used in this work, thoroughly described in [START_REF] Basile | Unstructured h-and hp-adaptive strategies for discontinuous Galerkin methods based on a posteriori error estimation for compressible flows[END_REF] for steady flows, is briefly reminded in this section and adapted to unsteady flows. Accurate DG-based a posteriori indicators extracted from the DG flow solver CODA are coupled with the remeshing library MMG. This error estimator identifies the regions lacking in accuracy, improving their resolution by either decreasing the size of the element or increasing the polynomial degree which approximates the solution. A smoothness indicator guides the hp-decision, leading to p-enrichment for smooth regions and h-refinement for non-smooth regions.

The hp-adaptive algorithm, initially based on fully simplicial meshes (triangles in 2D and tetrahedra in 3D) and then extended to hybrid-meshes, permits to take into account elements which cannot be remeshed. A variable polynomial degree setting adapts the resolution of the fixed elements by increasing their polynomial degree, overcoming the constraint that elements with a fixed geometry impose in a pure h-adaptive context. In these regions not subject to h-adaptation, the hp-choice is relaxed, and an element requiring increased resolution does not need anymore to be smooth to be p-adapted: a value of the error estimator greater than the target error estimator is sufficient to mark it for p-refinement.

The adaptation strategy can be applied similarly for static ℎ/ℎ𝑝-adaptation of unsteady flows. In this case the adaptive algorithm is applied once a pseudo-steady state of the flow is reached. Instead of the instantaneous values, the time-average of the error estimator accumulated for a given period is considered.

A. Error estimator

We use here an accurate and simple error estimator suited for various polynomial degrees. In previous work [START_REF] Basile | Unstructured h-and hp-adaptive strategies for discontinuous Galerkin methods based on a posteriori error estimation for compressible flows[END_REF], we followed the same idea of Colombo et al. [START_REF] Colombo | A p-adaptive implicit discontinuous Galerkin method for the under-resolved simulation of compressible turbulent flows[END_REF] and Bassi et al. [START_REF] Bassi | A p-adaptive matrix-free discontinuous Galerkin method for the implicit LES of incompressible transitional flows[END_REF] by combining two error estimators. One is based on the energy of the highest-order modes, the Small Scale Energy Density (SSED) 𝜖 SSED [START_REF] Naddei | A comparison of refinement indicators for p-adaptive simulations of steady and unsteady flows using discontinuous Galerkin methods[END_REF], while the other is based on the jumps across element interfaces, 𝜖 JUMP [START_REF] Bernard | High-order h-adaptive discontinuous Galerkin methods for ocean modelling[END_REF]. The reasons of this choice lie in the fact that an error indicator based on the highest order modes of the solution is only reliable for a high p, while a jump-based error estimator is accurate for every polynomial degree. The formulation of the error estimator presented in [START_REF] Basile | Unstructured h-and hp-adaptive strategies for discontinuous Galerkin methods based on a posteriori error estimation for compressible flows[END_REF], and used in the present work, is here recalled.

The SSED error estimator is based on a measure of the error E 𝐾 of the numerical solution u with respect to the exact solution u 𝑒𝑥 in the element by computing the norm of the difference between the numerical solution u ℎ, 𝑝 and the projection of the numerical solution on the reduced-order space V 𝑝-1 ℎ , namely u ℎ, 𝑝-1 , and dividing it by the volume of the element

𝜖 2 SSED,𝐾 = ∫ 𝐾 ||(𝜌v) ℎ, 𝑝 -(𝜌v) ℎ, 𝑝-1 || 2 𝑑𝑉 |𝐾 | = ||(𝜌v) ℎ, 𝑝 -(𝜌v) ℎ, 𝑝-1 || 2 𝐿 2 (𝐾) |𝐾 | (25)
where (𝜌v) ℎ is the momentum vector and |𝐾 | is the volume of the element.

The second error estimator is based on the assumption that, at a given interface 𝑒 of an element 𝐾 (an edge in 2D, a face in 3D), the error on the interface 𝑒 can be defined as half the jump of the variable traces across the interfaces.

Therefore we build an averaged error estimator over each interface 𝑒, for the norm of the momentum vector

𝜖 2 JUMP,𝑒 = ∫ 𝜕𝐾 𝑒 ||(𝜌v) + ℎ -(𝜌v) - ℎ || 2 𝑑𝑆 4|𝜕𝐾 𝑒 | = ||(𝜌v) + ℎ -(𝜌v) - ℎ || 2 𝐿 2 (𝜕𝐾 𝑒 ) 4|𝜕𝐾 𝑒 | (26)
and we normalize it by the area of each interface |𝜕𝐾 𝑒 | for consistency with the dimensions of the SSED indicator in equation [START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF].

For each element 𝐾 we compute the error estimator by the following simple rule:

𝜖 2 JUMP,𝐾 = 1 𝑁 𝑒 𝑁 𝑒 ∑ 𝑒=1 𝜖 2 JUMP,𝑒 (27)
with 𝑁 𝑒 the number of faces of the element in 3D (or edges in 2D

B. Smoothness indicator

The difficulty in hp-adaptive methods lies in the choice whether to adapt an element with h-refinement or penrichment. In this work we revise the approach adopted in previous work [START_REF] Basile | Unstructured h-and hp-adaptive strategies for discontinuous Galerkin methods based on a posteriori error estimation for compressible flows[END_REF], initially proposed for 1D by Mavriplis in [START_REF] Mavriplis | A posteriori error estimators for adaptive spectral element techniques[END_REF], exploiting the properties of turbulent flows. The assumption is that the decay rate of the spectrum of DG modal coefficients is related to the convergence rate of the solution. This information is exploited to evaluate the smoothness of the solution in the elements in the mesh. It is assumed that for a 1D Legendre expansion of coefficients 𝑞 = 0, ..., 𝑝 𝐾 the power decay of the modes 𝑞 can be expressed as:

𝐸 (𝑞) ≃ 𝐶𝑞 -𝜎 (29)
where 𝐸 (𝑞) is the energy associated with the mode 𝑞, 𝐶 and 𝜎 are constants determined by a least-squares best fit of log(𝑎 (𝑞) ) 2 vs. log(𝑞). The decay coefficient 𝜎 is then used as smoothness indicator. In particular, the modal energy is proportional to the square of the momentum 𝐸 (𝑞) ∝ 𝑎 (𝑞) 2 , where 𝑎 (𝑞) are the polynomial coefficients associated with the norm of the momentum for the 𝑞-th mode. High decay rates imply that the solution is smooth, while the solution deviates from analytical behaviour in presence of low decay rates.

In the context of a spectral dynamic modelling procedure for Large-Eddy Simulation, Chapelier et al. [START_REF] Chapelier | A spectral-element dynamic model for the Large-Eddy simulation of turbulent flows[END_REF] observed that high values of 𝜎 could be found in laminar or well-resolved regions, while low values are likely to appear in elements presenting an intense but poorly resolved small-scale activity. They used this evaluation of the energy decay as an estimation of the quality of the resolution in each element, allowing for adapting the intensity of the sub-grid dissipation locally. They evaluated a priori the threshold value of under-resolution 𝜎 𝑡 ℎ𝑟 from Taylor-Green Vortex (TGV) direct numerical simulation data at 𝑅𝑒 = 5000 and they obtained values very close to the theoretical slope -5/3 found in Fourier space in the inertial range. Therefore they were able to link the smoothness indicator to the universal -5/3 power law which uniquely describes the shape of the energy spectrum in the inertial subrange. We can exploit the same observations and use the value 𝜎 𝑡 ℎ𝑟 = 5/3 as the threshold between smooth and non-smooth behaviour in our work.

Elements with 𝜎 𝐾 > 5/3 are assumed to be smooth, and subject to p-enrichment; if this condition is not fulfilled h-refinement is used.

Since for 2D and 3D computations, several modal coefficients can contribute to the energy at index 𝑞, we need to gather the coefficients of the modes to retrieve one single value per index 𝑞.

The approach we follow in our work consists in computing the coefficient 𝑎 (𝑞) as the L 2 -norm of 𝑢 (𝑞) of the coefficients of the polynomial basis as:

𝑎 (𝑞) 𝐾 = √ ∑ 𝑙 ∈𝑑 𝑞 𝑢 (𝑙) 2 𝐾 ∀𝑞 ∈ (1, 𝑝 𝐾 ) (30)
The coefficient associated to 𝑞 = 0 represents the cell-averaged solution, and it can severely bias the decay rate. Therefore we chose not to employ the coefficient 𝑎 (0) in the log-log regression log(𝑎 (𝑞) ) 2 vs. log(𝑞). This means that we can compute the smoothness indicator only for high-order elements 𝑝 ≥ 2. Therefore in this work we decided to always fictitiously mark 𝑝 = 1 elements as smooth elements, always requiring 𝑝-adaptation. Moreover we use the smoothness indicator also as indicator of the convergence rate of the solution, which will be explained in detail in the following sections.

C. Strategy for unsteady flows

The strategy is to initially adapt the mesh with respect to a relatively affordable steady RANS simulation, which already captures some important features of the unsteady simulation. This mesh is then used to start the more expensive hybrid RANS/LES adaptation. This approach, where the accurate and expensive hybrid RANS/LES adaptation is performed starting from a mesh adapted to the cheaper RANS simulation, is preferred over adapting the mesh starting with hybrid RANS/LES simulations. In fact, an initial very coarse mesh would prevent the turbulent structures of the flow from developing, and could yield numerical instabilities and very poor quality results. Moreover it would dramatically increase the computational time for the whole adaptation process.

Starting from a RANS solution is also natural in the present hybrid RANS/LES context, as the wall-attached flow (and eventually other critical zones using a zonal DES approach) is solved in RANS mode, and a ZDES statistically-steady flow can be easily established from a RANS solution, after passing a transient period 𝑇 trans .

Once a statistically-steady solution has been reached, we are ready to accumulate the error estimator for a relatively short time period 𝑇 𝜖 𝜎 , with respect to the total simulation time needed for data analysis and statistics collection 𝑇 stats . This allows the user to perform the successive adaptation steps without accumulating flow statistics for a long time. In fact the error estimator is very sensitive to poorly-resolved zones, and does not need a very long time to become representative in such zones needing an improved resolution. Both the error estimator/smoothness indicator and the flow data are sampled at a constant sampling time step Δ𝑡 sampl .

After each adaptation and projection of the previous solution on the newly adapted mesh, the error and the smoothness estimator and flow data statistics are collected after a transient 𝑇 trans . This transient is necessary in order to let the solution adapt to the new mesh.

The transient period 𝑇 trans , the collection period of the error estimator and the smoothness indicator 𝑇 𝜖 𝜎 , and the sampling time interval Δ𝑡 sampl are generally case-dependent, as well as the total simulation duration where flow statistics are collected 𝑇 stats . Thus they must be chosen accordingly to features depending on the studied flow.

D. Degree and metric prescription

The method consists in performing successive simulations, with a mesh and polynomial degree refinement at each step. The methodology used in the present work is suited for hybrid meshes: during the adaptation process, boundary layer structured or pseudo-structured elements are kept unchanged, and the remeshing adaptive algorithm acts only for tetrahedra. No specific development is made to the MMG library, which preserves prisms (quadrilaterals in 2D) by default. For hybrid meshes we split the mesh into two parts

T 𝐾 ,𝑖 = T free,𝑖 ∪ T fixed,𝑖 (31)
where T free,𝑖 is the tetrahedral zone subject to remeshing, and T fixed,𝑖 is the prismatic layer, which cannot be remeshed. We define as well three non overlapping sets of mesh elements:

T 𝐾 ,𝑖 = T ℎ,𝑖 ∪ T 𝑝,𝑖 ∪ T =,𝑖 (32)
where T ℎ,𝑖 is the subset of the elements marked for h-adaptation, T 𝑝,𝑖 = T 𝑝,free,𝑖 ∪ T 𝑝,fixed,𝑖 is the subset of the elements marked for p-adaptation, which can include tetrahedral elements or fixed prismatic elements, and T =,𝑖 is the subset of the elements not marked for any adaptation. An element cannot be selected for both h-and p-adaptation.

The size prescription is based on the assumption that the local error 𝜖 𝐾 converges asymptotically to zero at a certain convergence rate 𝑚. The presence of a physical or geometrical discontinuity, yields the loss of the smooth theoretical convergence order 𝑚 = 𝑝 + 1. Here we exploit the fact that for non-smooth regions the convergence order of the solution is fixed by the discontinuity 𝜖 𝐾 ≃ 𝑂 (1).

We summarize here the steps to perform the adaptive procedure starting from a RANS-adapted mesh T 0 , for which the variable degree distribution 𝑝 0 and the characteristic lengths field ℎ 0 are available. We recall that the index 𝑖 identifies the step of the adaptive procedure.

1) A statistically-steady state solution is established on T 𝐾 ,𝑖 .

2) The a posteriori error estimator 𝜖 𝑖 𝐾 is computed on T 𝐾 ,𝑖 and collected for a period 𝑇 𝜖 𝜎 .

3) The smoothness indicator 𝜎 𝑖 𝐾 is computed on T 𝐾 ,𝑖 and collected for a period 𝑇 𝜖 𝜎 . 4) The target error 𝜖 * 𝑖+1 is chosen in order to fulfill the given refinement criterion with an iterative procedure. 5) Elements presenting fixed geometry are marked for p-enrichment:

T 𝑝,fixed,𝑖 = {𝐾 ∈ T fixed,𝑖 | 𝜖 𝐾 ,𝑖 > 𝜖 * 𝐾 ,𝑖+1 and 𝑝 𝐾 ,𝑖 + 1 ≤ 𝑝 max } (33)
6) Elements free to be h-or p-adapted, are marked for 𝑝-enrichment:

T 𝑝,free,𝑖 = {𝐾 ∈ T free,𝑖 | 𝜖 𝐾 ,𝑖 > 𝜖 * 𝐾 ,𝑖+1
9) The size of elements marked for ℎ-adaptation is adapted. The input size-field that MMG requires is node-based, however our error estimator, smoothness indicator and polynomial degree are element-based. Therefore a volume-weighted average is applied to these quantities. In particular, they are only computed from the elements surrounding the node 𝑛 which have not been marked for 𝑝-enrichment at the current step.

• The polynomial degree 𝑝 𝐾 ,𝑖 of elements 𝐾 ∈ T ℎ,𝑖 ∪ T =,𝑖 is averaged to the nodes 𝑝 𝑛,𝑖 .

• The error estimator 𝜖 𝐾 ,𝑖 of elements 𝐾 ∈ T ℎ,𝑖 ∪ T =,𝑖 is averaged to the nodes 𝜖 𝑛,𝑖 .

• The smoothness indicator 𝜎 𝐾 ,𝑖 of elements 𝐾 ∈ T ℎ,𝑖 ∪ T =,𝑖 is averaged to the nodes 𝜎 𝑛,𝑖 .

• The size ℎ * 𝑛 is imposed at the node 𝑛:

ℎ * 𝑛,𝑖+1 = ℎ 𝑛,𝑖 𝜖 * 𝑖+1 𝜖 𝑛,𝑖 1 𝑚 𝑛,𝑖 , ∀𝑛 ∈ T ℎ,free,𝑖 (37)
• A convergence rate 𝑚 𝑛 = 𝑝 𝑛 + 1 is imposed if the solution is smooth (𝜎 𝑛 > 5/3), and 𝑚 𝑛 = 1 if the solution is non-smooth (𝜎 𝑛 ≤ 5/3). • A user-defined factor, the maximum refinement factor 𝑟 ℎ is necessary to avoid uncontrolled size changes between two consecutive adaptation steps. In one adaptation step the edge can decrease its size by at most 𝑟 ℎ times, i.e. ℎ * 𝑛,𝑖+1 ≥ ℎ 𝑛,𝑖 /𝑟 ℎ . 10) The new metric is given to MMG, which generates the refined mesh T 𝐾 ,𝑖+1 . 11) The new polynomial degree map 𝑝 * 𝐾 ,𝑖+1 is projected from T 𝐾 ,𝑖 to T 𝐾 ,𝑖+1 thanks to a nearest neighbor interpolation. 12) The solution of T 𝐾 ,𝑖 is projected on the newly adapted mesh T 𝐾 ,𝑖+1 with an 𝐿 2 projection. We perform these steps until we reach a certain threshold, which in our case is the maximum number of dofs we want to compute.

Adopting this strategy, the regions presenting a higher error estimator value than the target, with smooth solution, 𝑝 = 1 approximation, or fixed elements, are 𝑝-enriched, while the regions with non-smooth solution (or already at their maximum polynomial degree 𝑝 max , where h-adaptation is enforced) are refined of a factor depending on the ratio between the target error and the value of the error estimator.

E. Choice of the target error

As target error 𝜖 * 𝑛,𝑖+1 , common to all nodes, we use the error which corresponds to a fixed increase of the number of degrees of freedom from one step to the following one. For this purpose, we introduce a parameter 𝑓 𝑟 .

Let N 𝐾 ,𝑖+1 be the number of dofs in the i-th mesh, and N * 𝐾 ,𝑖+1 be the number of dofs in the new (i+1)-th mesh. In the (i+1)-th adaptation step, the desired number of dofs for the new mesh is defined as N * 𝑖+1 = 𝑓 𝑟 • N 𝑖 dofs. In order to reach this number, we assume that the number of dofs 𝑁 * 𝐾 ,𝑖+1 in the new mesh can be computed thanks to the prescribed reduction of element sizes, and keeping into account a variable polynomial degree for the element 𝐾.

N * 𝑖+1 = ∑ 𝐾 N 𝑝 ( 𝑝 * 𝐾 ,𝑖+1 ) ℎ 𝐾 ,𝑖 ℎ * 𝐾 ,𝑖+1 𝑑 (38)
where 𝑁 𝑝 ( 𝑝 * 𝐾 ,𝑖+1 ) is the number of degrees of freedom associated to the desired polynomial degree in the element 𝐾, 𝑑 is the number of dimensions of the mesh, ℎ 𝐾 ,𝑖 is the actual characteristic sizes of the element 𝐾 at the adaptation step 𝑖 and ℎ * 𝐾 ,𝑖+1 is the desired characteristic size of the element 𝐾 at the step 𝑖 + 1. These element-based characteristic sizes are computed as an average of the lengths associated to the nodes 𝑛 of the element 𝐾.

The target error 𝜖 * 𝑖+1 is then defined with an iterative procedure. We use the bisection method from the Python package SciPy to solve the equation in 𝜖 * 𝑖+1 :

N * 𝑖+1 (𝜖 * 𝑖+1 ) -𝑓 𝑟 • N 𝑖 = 0 ( 39
)
which will ideally provide us an adapted mesh with N * 𝑖+1 dofs. In practice the adapted mesh will have a different number of dofs from the theoretical desired value, however we can expect this number to be close.

IV. Validation of ZDES on basic test cases

A. Influence of the 𝑦 + value on the flat plate case at Reynolds 5 • 10 6 Our goal is to determine the acceptable near-wall resolution in terms of 𝑦 + for high-order DG methods and provide a comparison with the best practices for classical second-order schemes. In this section we investigate the influence of the minimum 𝑦 + value of the mesh for three discretization orders 𝑝 = 1, 2, 3 on the flat plate test case. For this purpose RANS simulations are sufficient to assess the problem also for ZDES mode 1 simulations, since the boundary layer is fully solved in RANS mode.

The turbulent flat plate case is chosen as the "2D Zero Pressure Gradient Flat Plate Verification Case" for turbulence modelling by Nasa [START_REF]2D Zero Pressure Gradient Flat Plate Verification Case[END_REF]. The case is run at 𝑀 ∞ = 0.2, at a Reynolds number of 𝑅𝑒 𝐿 = 5 • 10 6 based on length 𝐿 = 1. The plot in figure 2 shows the layout of the flat plate grids used for this study, with the boundary conditions employed. The stagnation pressure 𝑝 𝑡 /𝑝 ∞ = 1.02828 and the stagnation temperature 𝑇 𝑡 /𝑇 ∞ = 1.008 are imposed at the inflow, non reflecting farfield Riemann conditions are imposed at the top boundary, and the static pressure 𝑝/𝑝 ∞ = 1 is imposed at the outflow. The flat plate at the bottom boundary is modelled as a solid adiabatic wall, while symmetric boundary conditions are used in the bottom wall part preceding the flat plate. The flat plate extends from 𝑥 = 0 to 𝑥 = 2𝐿 at 𝑦 = 0, the farfield is imposed at 𝑦 = 𝐿, while symmetry conditions extend from 𝑥 = -0.24𝐿 to 𝑥 = 0. The relevant minimum height of the first cell for a DG method is considered as Δ𝑦 𝐷𝐺 = Δ𝑑/( 𝑝 + 1), where Δ𝑑 is the height of the first element. For a FV method, the minimum wall normal spacing is computed as the distance of the cell center from the wall Δ𝑦 𝐹𝑉 = Δ𝑑/2. Since p1 simulations are second order accurate, the computation of the relevant height of the first cell is the same as classical FV.

Simulations with different wall normal spacings Δ𝑑 of the first element on the wall are performed. In the wall-normal direction the mesh is stretched with a linear progression, and the stretching factor 𝛼 is the parameter governing the wall normal spacing of the first element Δ𝑑 = 𝐿/(

𝑁 𝑦 -1 𝑖=0 𝛼 𝑖 ).
The study has been assessed on two mesh configurations. The first, called hereafter G1, consists in structured 𝑁 𝑥 × 𝑁 𝑦 = 224 × 97 grids, represented in figure 2. The streamwise spacing is kept at a constant value Δ𝑥 = 0.01 (equivalent to an average Δ𝑥 + ≃ 1800 over the flat platenotice that no DG scaling has been applied to this value, as it will be explained for 𝑦 + in the following) along the domain, as well as the number of elements in the mesh. G1 simulations are plotted in straight lines and circles in figures 3 and 4. For a DG method, several dofs are present in each element in each direction depending on the employed discretization order, thus G1 simulations provide a different number of dofs for each discretization (table 1).

The different global resolution between DG discretizations is likely to affect the study, as the results depend not only on the DG-relevant height of the first cell Δ𝑦 + 𝐷𝐺 , but also on the DG-relevant streamwise spacing Δ𝑥 + 𝐷𝐺 = Δ𝑥/( 𝑝 + 1). Therefore, we decided to conduct the study also defining the same number of dofs for the three discretizations, denominated here as G2. The p1 224x97 mesh is kept the same, with around 65K dofs, while a 158x69 mesh is *Note that p3 simulations in figure 3, both G1 and G2, present a numerical undershoot of the error on the 𝐶 𝐷 𝑓 𝑟𝑖𝑐 , for 𝑦 + 𝐷𝐺 exactly equal to 1. This behaviour is spurious and not representative of the real behaviour of the error on the 𝐶 𝐷 𝑓 𝑟𝑖𝑐 . Since results obtained by p2 and p3 simulations appear to converge to a precise error level for each discretization for 𝑦 < 1, we decided to employ in the last column of table 1, for both p2 and p3, the values of the error on the 𝐶 𝐷 𝑓 𝑟𝑖𝑐 at which each simulation converges, namely the error values obtained at 𝑦 + 𝐷𝐺 < 0. employed for p2 simulations and a 123x53 for p3 simulations for the G2 configuration. Also in this case the Δ𝑥 spacing is constant for each set of p1, p2, p3 meshes and the different wall normal spacing of the first cell are achieved by changing the stretching factor of the linear progression. Therefore all the G2 simulations are performed with around 65K dofs and are represented in dotted lines and crosses in figures 3 and4.

Different simulations are performed for p1, p2, p3 discretizations varying the height of the first cell, and the errors on the value of the friction drag coefficient over the flat plate, and on the skin friction coefficient at two different locations of the flat plate, 𝑥 = 0.1 and 𝑥 = 0.97, are monitored in figure 3 and4 for different values of the minimum 𝑦 + 𝐷𝐺 . In particular we use the minimum 𝑦 + 𝐷𝐺 averaged over the flat plate. Since for this test case the pressure drag coefficient is zero, the friction is the only contribution to the drag coefficient.

Typically, for wall-resolved RANS simulations, in order to ensure a proper representation of the velocity profile with accurate gradients, the chosen value for FV computations is Δ𝑦 + 𝐹𝑉 = 1 for the first layer. In our work we want to assess the possibility to employ higher values of Δ𝑦 + 𝐷𝐺 for higher-order methods. The first observation regarding figures 3 and 4 concerns the dependence of the 𝐶 𝑓 𝑟𝑖𝑐 , and consequently the 𝐶 𝐷 𝑓 𝑟𝑖𝑐 , on the number of dofs in the streamwise direction x. While the mesh G1 coincides with G2 for p1 computations, we can compare G1 and G2 curves for p2 and p3 computations and analyze the differences. For all the three plots, we notice that straight (G1) and dotted (G2) curves are superimposed for both p2 and p3 discretizations, for higher values of the Δ𝑦 + 𝐷𝐺 than 1. Below this threshold, the discrepancies between G1 and G2 results are due to the higher streamwise resolution of G1 simulations. This means that the convergence of the 𝐶 𝑓 𝑟𝑖𝑐 (and 𝐶 𝐷 𝑓 𝑟𝑖𝑐 ) value is guided only by the height of the first cell, except when dealing with very low values of the error, where a higher streamwise resolution leads to better results. Globally, looking at values of the error for Δ𝑦 + 𝐷𝐺 ≥ 1 no distictions are made between G1 and G2 meshes. This range of values is sufficient in this study, since we are not interested in very small cells with values of Δ𝑦 + 𝐷𝐺 lower than the classical value 1. In fact employing very refined meshes at the wall would impose severe constraints on the time step. The dependence of 𝐶 𝑓 𝑟𝑖𝑐 only on the height of the first cell is thus demostrated and the following observations will be made independently from the use of G1 or G2 meshes.

Looking at the blue curves in figures 3 and 4, in order to accurately capture the gradients in the near-wall region and achieve 1% errors on the 𝐶 𝑓 𝑟𝑖𝑐 (and 𝐶 𝐷 𝑓 𝑟𝑖𝑐 ), the requirement Δ𝑦 + 𝐷𝐺 ≃ 1 seems an adequate choice for the second order method. It follows that the second order DG discretization p1 imposes basically the same 𝑦 + constraint as a FV discretization.

While for Δ𝑦 + 𝐷𝐺 > 2.5 no substantial difference is seen among the three DG discretizations, the 𝐶 𝑓 𝑟𝑖𝑐 at two locations 𝑥 = 0.1 and 0.97, and the integrated 𝐶 𝐷 𝑓 𝑟𝑖𝑐 values, obtained with the higher order discretizations p2 and p3, appear to be less sensitive to the value Δ𝑦 + 𝐷𝐺 for lower values of Δ𝑦 + 𝐷𝐺 . The p2 and p3 discretizations achieve the same accuracy of the integrated 𝐶 𝐷 𝑓 𝑟𝑖𝑐 obtained by the second-order discretization at Δ𝑦 + 𝐷𝐺 = 1, with higher wall spacing respectively of Δ𝑦 + 𝐷𝐺 = 1.8 and Δ𝑦 + 𝐷𝐺 = 2.1 (figure 3). If we look at the error on the skin friction coefficient 𝐶 𝑓 𝑟𝑖𝑐 at 𝑥 = 0.1 in the left panel of figure 4, a value of Δ𝑦 + 𝐷𝐺 = 2.3 is sufficient for p2 and p3 discretizations to achieve the same accuracy compared to p1 simulations with Δ𝑦 + 𝐷𝐺 = 1, while Δ𝑦 + 𝐷𝐺 = 2.1 is sufficient for the 𝐶 𝑓 𝑟𝑖𝑐 at 𝑥 = 0.97. This leads in turn to a gain of almost one order of magnitude of accuracy for p2 and p3 simulations at Δ𝑦 + 𝐷𝐺 = 1 with respect to p1 at the same Δ𝑦 + 𝐷𝐺 = 1 . The last two columns of table 1 show the Δ𝑦 + 𝐷𝐺 values necessary to obtain a 1% error on the 𝐶 𝐷 𝑓 𝑟𝑖𝑐 and the error on the 𝐶 𝐷 𝑓 𝑟𝑖𝑐 at Δ𝑦 + 𝐷𝐺 = 1 for the three discretizations. The curve representing errors in figures 3 and 4 stagnate below Δ𝑦 + 𝐷𝐺 ≃ 1 due to other sources of errors (mostly the singular point of the flat plate) which do not depend only on the height of the first cell, but also on the streamwise resolution. The spurious undershoots that can be seen in the three figures (for G1p3 in figure 3, G1p2 in the left panel of figure 4 and G2p2 in the right panel of figure 4) are due to the coefficient oscillating around the reference value. This can lead to a fictitious yet localized sudden decrease of the error for very low error levels.

In conclusion, when handling high order discretizations, it appears from our analysis that the classical constraint of Δ𝑦 + 𝐷𝐺 = 1 can be relaxed, and a higher value of Δ𝑦 + 𝐷𝐺 ≃ 2.5 seems adequate in order to retrieve results similar to second order discretizations with Δ𝑦 + 𝐷𝐺 ≃ 1. Equivalent DG wall normal spacings for high-order methods can be then chosen more than twice than those employed for second order methods, while maintaining a good stability of time-explicit gradient-accurate simulations.

B. Isotropic turbulence decay on the Taylor-Green vortex at Reynolds 5000 on unstructured grids

In this section we assess the value of the constant of the DES model 𝐶 𝐷𝐸𝑆 in the context of DG simulations for unstructured tetrahedral grids. Shur et al. in [START_REF] Shur | Detached-eddy simulation of an airfoil at high angle of attack[END_REF] found the optimal value of the 𝐶 𝐷𝐸𝑆 which verifies a Kolmogorov inertial range up to the grid cutoff for simulations of an isotropic turbulent case and using a centered, fourth-order accurate differencing scheme. In the same spirit, we choose here to compute the 3D Taylor-Green vortex case to assess the capacity of DES to predict the physics of turbulence for different discretization orders of a modal DG method. This case describes freely decaying isotropic turbulence in a periodic box with length of each side 𝐿 𝑥 = 2𝜋𝐿 with 𝐿 = 1. The initial conditions with Mach 𝑀 0 = 0.1 are given by:

𝑝 = 𝜌 0 𝑉 2 0 1 𝛾𝑀 2 0 + 1 16
• (cos (2𝑥/𝐿) + cos (2𝑦/𝐿)) • (cos (2𝑧/𝐿) + 2)

𝑢 = 𝑉 0 • sin (𝑥/𝐿) • cos (𝑦/𝐿) • cos (𝑧/𝐿) (40
𝜌 ν = 𝜌(𝐶 𝑆 • Δ) 2 • √ 𝑆 𝑖 𝑗 𝑆 𝑖 𝑗 (41)
where 𝐶 𝑆 = 0.17 is the Smagorinsky constant, Δ = 3 √ 𝑉𝑜𝑙 𝑝+1 is the size of the mesh divided by the number of 1D degrees of freedom for a DG method, and 𝑆 𝑖 𝑗 is the strain rate tensor obtained analytically from equation [START_REF] Lorteau | Turbulent jet simulation using high-order DG methods for aeroacoustic analysis[END_REF].

Resolution NDOFS ≃ 72 3

We decided to perform simulations at relatively low fixed number of degrees of freedom 𝑁 𝐷𝑜𝐹𝑠 ≃ 72 3 (number is not exact due to the tetrahedral unstructured meshes for which an exact control of the number of elements is not straightforward) for different discretization orders, at 𝑝 = 1, 𝑝 = 2, 𝑝 = 3, 𝑝 = 4, 𝑝 = 5.

Different values of the constant 𝐶 𝐷𝐸𝑆 are tested for the 𝑝 = 3 discretization. Given that p1, p2, p3, p4, p5 discretizations count respectively 4, 10, 20, 35, 56 degrees of freedom per element, meshes including respectively around 72 3 /4, 72 3 /10, 72 3 /20, 72 3 /35, 72 3 /56 elements are employed for each discretization (figure 5). Note that the number of dofs per element is different from section IV.A because here the problem is three-dimensional. The quantities we monitor in order to assess the quality of the representation of the turbulence are the enstrophy 𝜉 and the kinetic energy 𝐸 𝑘 .

𝜉 = 1 |Ω| ∫ Ω 𝝎 • 𝝎 2 𝑑Ω ( 42
)
𝐸 𝑘 = 1 |Ω| ∫ Ω v • v 2 𝑑Ω ( 43
)
The reference DNS simulations have been performed on a mesh employing 480 3 dofs [START_REF] Chapelier | A spectral-element dynamic model for the Large-Eddy simulation of turbulent flows[END_REF] with a eighth order Spectral Difference method. For a fair comparison between the DES and the reference DNS computations, especially for the levels of enstrophy which are strongly affected by the resolution, the DNS is filtered using the transfer function associated with the DG discretization introduced by Chapelier et al. [START_REF] Chapelier | Development of a multiscale LES model in the context of a modal discontinuous Galerkin method[END_REF]. The DNS velocity and vorticity fields are filtered and the corresponding spatial averages are computed with a sampling Δ𝑡 𝑐 = 1 between 𝑡 𝑐 = 1 and 𝑡 𝑐 = 20. We compute the energy spectra of the solution at time 𝑡 = 10, when all scales of the flow are developed.

When analysing LES results obtained with DG methods it is assumed that the resolution limit is defined by the cut-off wavenumber 𝑘 1/3 = 3 √ NDOFS/3. This is justified by the fact that the DG discretization represent accurately the spectrum for wavenumbers up to 𝑘 1/3 = 3 √ NDOFS/3, followed by a steep decrease of the energy for higher wavenumbers [START_REF] Naddei | Spectral and modal energy transfer analyses of LES using the discontinuous Galerkin method and their application to the Variational Multiscale approach[END_REF]. This is only due to the properties of the DG approximation. The presence of extra numerical dissipation might lead to additional discrepancies between DNS and LES also at lower wavenumbers. The theoretical slope of the energy spectrum 𝑘 -5/3 is clearly observed in the inertial range for 𝑝 = 3, 4, 5 simulations, with a good fit with the DNS reference computation. The damping of energy occurs near a wavenumber of 14 for the three discretizations 𝑝 = 3, 4, 5, before the limit value 𝑘 1/3 = 72/3 = 24. The more the discretization order is increased, the later the cut-off occurs and the decay of the spectrum is slower. However the three curves are very close to each other, meaning that the maximum wave length resolution limit is the almost the same for 𝑝 ≥ 3, and no further significant improvements can be expected increasing the discretization order.

The p2 simulation presents an earlier damping than p3 but still maintaining the theoretical slope, while p1 simulation provides results far from the theoretical behaviour, damping a significant part of the energy associated with the large scales. This is due to the intrinsic numerical dissipation that low order DG schemes provide. The highest the order of the DG discretization, the lower the numerical dissipation. This means that the numerical dissipation is too strong for a low-order scheme to correctly represent the turbulent behaviour, unless the resolution is dramatically increased.

It is thus preferable to employ high-order 𝑝 ≥ 2 discretizations when an accurate representation of turbulence is sought with a fixed number of degrees of freedom.

For hp-adaptation to be attractive from the computational time point of view, an ad hoc load balancing technique is mandatory to fully exploit the potential of variable order degree simulations. Since this goes beyond the scope of this paper, and no sophisticated load balancing techniques are employed here, we choose to employ a maximum polynomial degree discretization 𝑝 = 3 in the adaptive algorithm for the PPRIME nozzle in section V. As shown in figure 6, the high-order 𝑝 = 3 discretization provides a good representation of turbulent phenomena, on par with higher-order discretizations.

All the following analysis are then conducted on a 𝑝 = 3 discretization, representative of our final test case.

The sensitivity of the model to the constant 𝐶 𝐷𝐸𝑆 is studied and the corresponding results displayed in figure 7. The results of DES simulations are presented in terms of the evolution of the enstrophy and the kinetic energy along the time of the simulation, with respect to filtered DNS data, and the 3D energy spectrum with respect to the DNS spectrum. Analyzing the plots of the enstrophy over the course of the simulation in the left panel of figure 7, the decrease of the 𝐶 𝐷𝐸𝑆 constant from the original value 𝐶 𝐷𝐸𝑆 = 0.65 to 𝐶 𝐷𝐸𝑆 = 0.01, 0.1, 0.2, 0.4, 0.5 leads to a slightly better representation of the peak of enstrophy, due to a weaker dissipation introduced by the DES model. At the same time, the kinetic energy evolution in the middle panel of figure 7 shows a low sensitivity to variations of the DES constant, meaning that the model predicts accurately the decay of the kinetic energy employing any constant.

The value of the constant 𝐶 𝐷𝐸𝑆 seems to affect marginally the energy spectra in the right panel of figure 7. The simulations performed with different 𝐶 𝐷𝐸𝑆 values provide indeed a similar decay of the turbulent energy spectrum.

The method seems robust enough according variations of 𝐶 𝐷𝐸𝑆 , and the standard value 𝐶 𝐷𝐸𝑆 = 0.65 does not appear to introduce an excessive amount of dissipation.

Mesh sensitivity study

In this subsection we assess the dependence of DG-DES simulations on the mesh resolution for a fixed polynomial degree. We perform simulations on different resolutions, and compare them with filtered DNS data at the respective cut-off wave number 𝑘 1/2 . The 𝑝 = 3 discretization is considered, and simulations counting NDOFS = 48 3 , 96 3 , 144 3 are performed in addition to NDOFS = 72 3 , on tetrahedral meshes. Three different 𝐶 𝐷𝐸𝑆 = 0.1, 0.4, 0.65 are compared.

We made comparisons for the enstrophy development for the three resolutions. Kinetic energy is not shown here because only very small differences can be found between the different computations. Note the different scale of the enstrophy 𝜉 levels for the three resolutions.

As already observed for the resolution 72 3 , increasing 𝐶 𝐷𝐸𝑆 leads to a stronger dissipation and lower levels of the enstrophy. A low value of the contant 𝐶 𝐷𝐸𝑆 = 0.1 does not provide a sufficient dissipation after the peak of enstrophy has been reached. In particular higher levels of the enstrophy than the respective filtered DNS are observed for 𝑡 𝑐 > 12 for the two coarse resolutions 48 3 (left panel in figure 8) and 72 3 (left panel in figure 7). On the other hand, when increasing the resolution to 96 3 and 144 3 , the model always dissipates enough the energy due to the unresolved scales, without needing high values of the 𝐶 𝐷𝐸𝑆 .

The energy spectrum is also analyzed to provide a better insight in the mechanism of the turbulent decay. As already pointed out for the resolution 72 3 , the energy spectra at 𝑡 𝑐 = 10 obtained with 3 different values of the 𝐶 𝐷𝐸𝑆 are almost indistinguishable, for each resolution. This confirms the overall low criticity of this constant for the generation and development of turbulent scales, as demonstrated by Shur et al. [START_REF] Shur | Detached-eddy simulation of an airfoil at high angle of attack[END_REF] for decaying homogeneous isotropic turbulence.

According to Sagaut et al. in [START_REF] Sagaut | Multiscale and multiresolution approaches in turbulence-LES, DES and Hybrid RANS/LES Methods: Applications and Guidelines[END_REF] the DES approach can be considered analogous to the well-known Smagorinsky model for free-shear flows, and a value of 𝐶 𝐷𝐸𝑆 = 0.65 is equivalent to the Smagorinsky constant 𝐶 𝑆 ≃ 0.2. Due to the dissipation of DG methods, a lower 𝐶 𝑆 = 0.1 is often employed [START_REF] De La Llave Plata | On the performance of a high-order multiscale DG approach to LES at increasing Reynolds number[END_REF] for LES/DG simulations. We could think that decreasing the constant to 𝐶 𝐷𝐸𝑆 = 0.3 or 0.4 would be an appropriate choice to allow the subgrid model to provide a sufficient amount of dissipation, without over-dissipating the turbulent structures. This value would correspond to a 𝐶 𝑆 ≃ 0.1 according to Sagaut's relations. However, for the purposes of this work we prefer to use the classical constant value 𝐶 𝐷𝐸𝑆 = 0.65, since very few differences are brought by changing the 𝐶 𝐷𝐸𝑆 value, and no differences are found out in the decay of the turbulent spectrum. A much wider analysis would be necessary to assess the effective benefits that changing this constant would bring, and eventually propose another value.

V. Application to a turbulent configuration: the isothermal subsonic turbulent jet flow issued

from the PPRIME nozzle at 𝑅𝑒 𝐷 = 10 6 .

In this section, we apply the hp-methodology initially introduced in [START_REF] Basile | Unstructured h-and hp-adaptive strategies for discontinuous Galerkin methods based on a posteriori error estimation for compressible flows[END_REF] and extended to unsteady flows in section III, to ZDES computations of a nozzle configuration. The nozzle test case analysed in this work corresponds to the configuration described in [START_REF] Brès | Importance of the nozzle-exit boundary-layer state in subsonic turbulent jets[END_REF], for which experiments have been performed at the Bruit & Vent jet-noise facility of the Institut PPRIME, Poitiers, France, with reference data available. Numerical results are available as well in the literature: this configuration has been simulated with LES by Brès et al. [START_REF] Brès | Importance of the nozzle-exit boundary-layer state in subsonic turbulent jets[END_REF] to analyse features as near-wall adaptive mesh refinement, synthetic turbulence and wall modelling inside the nozzle in LES, with hybrid RANS/LES by Gand et al. in [START_REF] Gand | On the generation of turbulent inflow for hybrid RANS/LES jet flow simulations[END_REF] to assess the generation of a turbulent inflow and with RANS by Neifeld et al. [START_REF] Neifeld | Jet noise prediction with Eddy relaxation source model[END_REF] for jet noise prediction purposes with eddy relaxation source model.

A. Computational set-up

The operating conditions are defined in terms of the total pressure ratio 𝑝 𝑡 /𝑝 ∞ = 1.7 and total temperature ratio 𝑇 𝑡 /𝑇 ∞ = 1.15. We refer to jet properties with the subscript 𝑗. The jet is isothermal (𝑇 𝑗 /𝑇 ∞ = 1.0), the jet Mach number is 𝑀 𝑗 = 𝑈 𝑗 /𝑐 𝑗 = 0.9, and the Reynolds number is 𝑅𝑒 𝐷 = 𝜌 𝑗 𝑈 𝑗 𝐷/𝜇 𝑗 ≃ 1 • 10 6 , where 𝑈 𝑗 is the mean jet exit longitudinal (x-direction) velocity, 𝑐 is the speed of sound, 𝜌 is the density and 𝜇 is the dynamic viscosity. The characteristic time scale is 𝑡 𝑐 = 𝐷/𝑈 𝑗 . A slow coflow at Mach number 𝑀 ∞ = 0.01 is imposed outside the nozzle (the flow is at rest in the experiment) to prevent spurious recirculation and facilitate flow entrainement.

All solid surfaces of the nozzle are treated as no-slip adiabatic walls. Farfield boundary conditions are set to stagnation pressure/stagnation temperature inflow at the inlet, static pressure outflow at the outlet of the domain, and non-reflecting free-stream to the lateral surfaces of the external cylindrical domain. The axisymmetric computational domain extends from approximately -10D (the longitudinal length of the nozzle in the x-direction) to 50D in the longitudinal direction x, and from -30D to 30D in the radial direction. The nozzle geometry is fully considered and its exit is centred at (𝑥, 𝑦, 𝑧) = (0, 0, 0).

A prismatic boundary layer is extruded in the direction normal to the surface of the nozzle, and the rest of the domain is filled with tetrahedra. The prismatic boundary layer is kept fixed during the adaptation process, since pure metric-based mesh adaptation can handle only simplicial meshes (triangles in 2D and tetrahedra in 3D). These elements can only be p-adapted.

In addition to the prismatic boundary layer, it has been chosen to preserve also the tetrahedral elements internal to the nozzle. Tetrahedra free to change size inside the nozzle but constrained by the fixed prismatic surfaces can severely affect the quality of the mesh. This can be easily done thanks to the ability of MMG to preserve tetrahedra specified by the user. These elements will be subject to p-adaptation only.

The hybrid RANS/LES adaptive simulations are initialized from RANS adaptive simulations: this means that some steps of steady RANS adaptations are needed before performing the first unsteady turbulent simulation. The steady RANS adapted mesh will be used as starting point for unsteady adaptations, for which very coarse initial meshes would prevent the turbulent structures of the flow from developing.

The initial mesh hpG1 employed in this work is an hp-mesh obtained from a steady hp-adaptation using RANS equations to solve the nozzle/jet flow [START_REF] Basile | Unstructured h-and hp-adaptive strategies for discontinuous Galerkin methods based on a posteriori error estimation for compressible flows[END_REF]. Three previous steps of RANS adaptation have been performed starting from a very coarse mesh counting 1.5 million dofs. The hpG1 mesh, counting 11.6 million dofs and attaining a sufficient level of convergence with RANS equations, is used as starting point for ZDES mesh adaptation. Two adaptation steps are performed, yielding a discretization hpG2 with 18.9 million dofs, and a discretization hpG3 with 31.7 million dofs.

Regarding unsteady turbulent simulations, a sponge layer is created for elements presenting their cell centers more than 35 diameters far from the nozzle exit in the streamwise direction 𝑥 𝑐 > 35𝐷, and 7 diameters from the axis of the jet 𝑟 𝑐 > 7𝐷 in the radial direction. The mesh elements in these regions have a very large size which can compromise the stability of the simulations when turbulent structures or high amplitudes sound waves pass through them, if a high order of accuracy is defined in those cells. An elegant way fitting to the present hp-methodology consists in forcing these elements to a low order p1 discretization. This creates a sponge layer which further damps the aforementioned problematic phenomena as well as reducing the overall computational burden of the simulations. This forcing is kept active through the adaptive procedure (figure 10). Flow statistics are collected on a post-processing cloud of points. At each Δ𝑡 sampl the solution is reconstructed on the post-processing mesh and the temporal statistical sampling is performed using previous data. The profiles shown in the next sections are averaged over 400 post-processing points in the azimuthal direction, for each radial coordinate. This can be done due to the azimuthal symmetry of the geometry, for which these locations are statistically equivalent. This azimuthal average eases statistical convergence and provides improved results. Only the jet centerline velocity profiles cannot be azimuthally averaged, and this makes this quantity harder to converge, needing longer simulation times.

B. Simulation/adaptation parameters

The time intervals required by the transient and the collection of error estimator, introduced in section III.C and expressed in characteristic time 𝑡 𝑐 , are chosen as follows:

• Δ𝑡 sampl = 0.0067𝑡 𝑐 , which means 150 samplings for each

𝑡 𝑐 • 𝑇 𝜖 𝜎 = 30𝑡 𝑐 • 𝑇 stats = 150𝑡 𝑐
• 𝑇 trans = 50𝑡 𝑐 In order to assess the quality of the intermediate adaptation steps, the flow statistics are accumulated as well for the ℎ𝑝-meshes hpG1 and hpG2 for 𝑇 stats = 150𝑡 𝑐 , whereas in order to obtain only the final discretization hpG3 this time could be reduced to 𝑇 𝜖 𝜎 = 30𝑡 𝑐 which is sufficient to obtain statistical convergence for the error estimator.

The parameters used by MMG and the adaptation module are:

• ℎ 𝑔𝑟 𝑎𝑑 = 1.5, showing a good compromise between excessive refinement induced by a low value of ℎ 𝑔𝑟 𝑎𝑑 , and the loss of isotropy and mesh quality induced by a high value. • ℎ 𝑚𝑖𝑛,tetra is set to ℎ 𝑚𝑖𝑛 = 0.016𝐷 for tetrahedra, which is the length of the smallest edge of the boundary layer prisms which constitutes the interface with tetrahedra. The height of the first element of the prismatic layer is ℎ 𝑚𝑖𝑛,prisms = 0.0003𝐷, and is fixed along the geometry and over the course of the adaptive procedure. • ℎ 𝑚𝑎𝑥 is set as the size of the domain where far-field boundary conditions are imposed.

• 𝑟 ℎ = 4 such that between two steps of the adaptation process, the edge of the equilateral triangle can at most divide by four its size. • 𝑓 𝑟 = 2 such that the desired number of dofs in the newly adapted mesh is twice the number of dofs present in the previous mesh. We remind that the adaptive procedure, together with the remeshing step performed by MMG, is not always able to respect the exact desired number of dofs, while however providing a good guess. • 𝑝 min = 1 is the minimum polynomial degree of the spatial discretization allowed for hp-adaptation.

• 𝑝 max = 3 is the maximum polynomial degree that the spatial discretization is allowed to reach for hp-adaptation.

C. Results

Computational meshes

Figure 11 shows the three ℎ𝑝-meshes employed in this study: hpG1, the initial hp-mesh of the ZDES adaptive procedure, adapted from a previous RANS adaptation; hpG2, the first adapted hp-mesh from a ZDES simulation; hpG3, the second adapted hp-mesh from a ZDES simulation.

The first hp-mesh is almost entirely discretized with p2 and p3 elements in the jet plume, except for the exit of the nozzle, which is discretized with p1 elements. The external part of the shear layer present p2 elements since the smoothness estimator in previous RANS simulations had detected non-smooth features in this region, and had prevented the algorithm from increasing the polynomial degree in this zone.

We see a concentration of the elements around the potential core. The zone relevant for unsteady ZDES simulations is the mixing layer zone where the vortices develop from the lip of the nozzle, and it is discretized with the maximum polynomial degree 𝑝 = 3. As seen from TGV tests in section IV.B, this property is desirable because high-order DG discretizations describe the turbulence more accurately than low order discretizations with an equivalent number of degrees of freedom.

hp-adapted meshes hpG2 and hpG3 reduce the number of p2 elements in the jet plume, increasing their polynomial order to the maximum p3, especially on the axis of the jet. Globally, the polynomial degree distribution remains the same for the three meshes, since the RANS adaptation had already individuated the important features of the flow and provided a good initial discretization.

The most important difference between the three hp-meshes resides in the size of the elements in the jet plume. The mixing layer, presenting strong vortical structures needing very fine resolutions to be well captured, is progressively h-refined, since the maximum polynomial degree has already been attained. Higher degree computations would need ad-hoc load balancing techniques, which are not the subject of the present paper, but will be addressed in future work.

Qualitative instantaneous field analysis

Before quantitatively assessing the numerical results, the visualizations of the instantaneous flow field for the three ℎ𝑝-meshes hpG1, hpG2, hpG3, are presented in figure 12 and 13 and compared to perform a first qualitative assessment.

Figure 12 describes the instantaneous fields of temperature (colour) and pressure fluctuations (grey scale), for the meshes hpG1 with 11.6 millions dofs (top), hpG2 with 18.9 millions dofs (middle) and hpG3 with 31.7 millions dofs (bottom). These pictures show how the most resolved computation captures smaller turbulent structures compared to the first and the second resolutions, thanks to the hp-adaptation in the zones of interest of the jet. Moreover we notice how the use of a finer grid resolution in the jet flow leads to a slower jet development (and a longer potential core). This trend will be further discussed in the next section.

In figure 13, vorticity snapshots are shown for three different sections 𝑥/𝐷 = 1, 2, 3 downstream the nozzle (plots from the left to the right), for hpG1, hpG2, hpG3 computations (top, middle, bottom rows). hpG1 resolves larger structures than hpG2 and hpG3 resolutions for 𝑥/𝐷 = 1, while solutions on hpG2 and hpG3 are almost indistinguishable. This means that the adaptive algorithm has converged to the same solution in hpG2/hpG3 simulations for the zone close to the exit nozzle. More consistent differences can be found for locations more downstream the exit nozzle. While for hpG1 at turbulent structures are barely captured at 𝑥/𝐷 = 2, 3, small structures are well predicted by hpG2 and, to a larger extent, by hpG3.

Quantitative data analysis

The boundary layer velocity profile is extracted in figure 14 at 𝑥/𝐷 = 0.04. For the three simulations, the velocity profiles for hpG1, hpG2, hpG3, respectively in blue, red, black straight lines, match experimental data in squared symbols and the RANS simulation on hpG1 in dotted lines. Note that in the proximity of the nozzle exit, no substantial differences can be found between hpG1, hpG2, hpG3 hybrid RANS/LES results, because the flow features here mainly depend on zones solved with RANS equations, which have already grid-converged.

Even though a turbulent velocity profile is obtained at the nozzle exit due to the RANS modelling inside the nozzle, no turbulent injection has been used in this work, thus some discrepancies between computational and experimental results (where turbulence is triggered inside the nozzle) are expected.

The mean streamwise velocity profiles on the jet axis 𝑟/𝐷 = 0 and on the lipline 𝑟/𝐷 = 0.5, are compared to experimental and numerical results obtained by Brès et al. [START_REF] Brès | Importance of the nozzle-exit boundary-layer state in subsonic turbulent jets[END_REF] in green lines, Gand et al. [START_REF] Gand | On the generation of turbulent inflow for hybrid RANS/LES jet flow simulations[END_REF] in magenta lines, Shur et al. [START_REF] Shur | LES-based evaluation of a microjet noise reduction concept in static and flight conditions[END_REF] in orange lines, and Lorteau et al. [START_REF] Lorteau | Turbulent jet simulation using high-order DG methods for aeroacoustic analysis[END_REF] in cyan lines, whose parameters and numerical set-up can be found in table 2. Dotted lines represent the coarse grid for each reference, while dashed lines represent fine grids. Shur and Lorteau show the peak of rms velocity instead of the lipline velocity. However the peak of the rms velocity is found very close to 𝑟/𝐷 = 0.5, and the two quantities can be directly compared.

Brès et al. used a second order Finite Volume approach on unstructured isotropic hexahedral-dominant grids, obtained with a a priori mesh adaptation procedure, where it is not an error estimator that guides the process, but an a priori knowledge of the physical phenomenon. Starting from a structured mesh with a low number of elements, several embedded zones of refinement with specific target length scale are defined by the user and enforced iteratively by the adaptation tool. Moreover further mesh refinement is added inside the nozzle in the near-wall zone. Their LES approach is based on the Vreman subgrid model to account for the physical effects of underresolved turbulence and on an equilibrium wall model inside the nozzle, which models the near-wall turbulence and allows for selecting element sizes larger than 𝑦 + = 1 at the wall. They also used synthetic-turbulence boundary conditions to model the boundary-layer trip present in the experiment on the internal nozzle surface.

Gand et al. adopted hexahedral structured meshes with a second-order FV scheme, modelling the turbulence with ZDES mode 2 and ZDES mode 3 approaches. Here we consider only their ZDES mode 2 results, since both ZDES mode 1 employed in this work and ZDES mode 2 aim at modelling the boundary layers inside the nozzle with a RANS approach without a tripping of turbulent structures and thus account for a similar modelling approach.

Shur et al. employed high-order FV schemes on structured meshes. They used a two-stage simulation procedure in which the inflow velocity profile was imposed from a RANS solution previously obtained on a coupled nozzle/jet plume RANS simulation, excluding the geometry from the ILES computation. This allowed them to reproduce the effect of the boundary layer considerably reducing the high cost of a coupled nozzle/plume LES. Because of this, a direct comparison between the number of dofs employed in their simulations and in our adaptive DG hp-meshes cannot be done in a straightforward manner, however it is still interesting to compare our results with those obtained by Shur as in both cases the nozzle exit profile is obtained from RANS modelling approaches.

Lorteau et al. adopted a fourth-order DG approach on fully tetrahedral meshes. No tripping procedure was employed, and a LES-Smagorinsky model was used to model turbulence. Their DG simulation results (G1) are shown in figures For each of the four numerical references, results obtained on two grid resolutions (except for Lorteau for which we have only one DG resolution) are extracted and compared to present hp-adaptive results. The first observation that can be made is that refining (in h and p) the jet plume, the velocity decay on the jet axis becomes slower, leading to longer potential cores (as already observed in figure 12). Surprisingly, we can observe how the adaptive simulation hpG2 seems to provide closer results to the experiment, while the more refined hpG3 mesh overestimates the length of the potential core. Even though this results could seem unexpected, this behaviour has been already pointed out by Shur et al. in [START_REF] Shur | LES-based evaluation of a microjet noise reduction concept in static and flight conditions[END_REF] (in orange lines), and can be found in results from Brès et al. [START_REF] Brès | Importance of the nozzle-exit boundary-layer state in subsonic turbulent jets[END_REF] (in green lines) and Gand et al. [START_REF] Gand | On the generation of turbulent inflow for hybrid RANS/LES jet flow simulations[END_REF] (in magenta lines) as well. The trend is clear, even though for Shur and Gand this is more marked than for Brès. A similar overestimation of the length of the potential core is also observed by Lorteau, and Uzun et al. [START_REF] Uzun | High-fidelity numerical simulation of a chevron nozzle jet flow[END_REF] for very fine grids of a chevron nozzle jet flow.

The reasons for this overestimation of the jet potential core are still not clear. Shur hypothesized two reasons for this behaviour: one is that it might be caused by the use of ILES, while the other is that it might be caused by the feedback loop between jet turbulence and nozzle not accurately reproduced by the system of boundary conditions. The second observation would be valid also for our ZDES mode 1 model, because the small-scale variations of the turbulent flow cannot go back upstream and interact with the interior of the nozzle, which is modelled in RANS. Lorteau et al. attributed turbulence within the nozzle, as it happens in the experiment. Even though a formal answer has not been found, we think that a combination of several effects could be the cause of this behaviour.

In particular the lack of turbulent injection: the turbulent transition here occurs in the shear layer, outside the nozzle. This consists in the processes of vortex rolling-up and pairing, after which three-dimensional turbulence appears. This may affect both the flow and the acoustic solution in different ways [START_REF] Lorteau | Turbulent jet simulation using high-order DG methods for aeroacoustic analysis[END_REF][START_REF] Brès | Importance of the nozzle-exit boundary-layer state in subsonic turbulent jets[END_REF][START_REF] Gand | Investigation of turbulence development in incompressible jets with zonal detached eddy simulation (ZDES) and synthetic turbulent inflow[END_REF]. We remind here that the experiment present a fully turbulent boundary layer, after triggering the turbulence inside the nozzle with a strip. A similar tripping technique could be considered as future work in our simulations in order to elucidate its effect on the jet plume flow dynamics and check if the simulations converge towards the experimental results when the mesh is refined or adapted.

It is also possible that the resolution is still insufficient. In this case, a fourth more refined simulation (not performed because of limited computational resources) should be performed in order to assess the flow properties in terms of mean and fluctuating profiles and verify the trends observed in the previous three hp-simulations.

Compared to reference simulations, the axial mean velocity is in very close accordance with Gand and Shur more refined simulations, while for 𝑥/𝐷 > 10 the decay of the velocity is faster for hpG3. This is due to the a lower resolution in this zone with respect to these two studies, suggesting that a further adaptation step providing additional refinement for 𝑥/𝐷 > 10 would eventually slow the decay of the centerline velocity.

An observation of the mean velocity profile on the lipline in the right panel of figure 15, reveals the impact that grid refinement has in the development of the shear layer is evident, especially for 𝑥/𝐷 > 2. Here grid refinement yields to an increased velocity along the lipline, obtaining closer results to the experiment. In figure 16 we show rms streamwise velocity profiles in the centerline (left) and in the lipline (right). Looking at the centerline velocity profiles, we notice that the profile is not smooth, but presents discontinuities. These discontinuities are evident for the coarse simulation hpG1, while they tend to decrease in intensity thanks to the adaptive grid refinement. This aspect can originate from the discontinuous nature of the DG scheme, where the solution is discontinuous across elements, but also possibly from the statistical sampling which is low at the centerline due to the impossibility of performing a spatial azimuthal averaging at this location. Also, especially in the proximity of interfaces of very coarse elements, the polynomial approximating the solution can have very sharp shapes. This behaviour is improved thanks to mesh adaptation, and peaks in the solution gradually disappear for hpG2 and hpG3. The eventual discontinuities in the solution are smoothed out for lipline plots in the right panels of figures 15 and 16, and velocity profiles in figures 17 and 18, by averaging the solution in the azimuthal direction, which is not possible for the centerline profiles.

The small peak in the lipline velocity rms profile around 𝑥/𝐷 = 0.5 is very similar for hpG2 and hpG3, meaning that in the zone close to the nozzle exit the two meshes display the same flow features and have reached a grid-converged solution. Poor resolution just downstream of the jet exit seems to strengthen the mixing-layer transition, resulting in higher axial (in the centerline) and peak values (in the lipline) of turbulence intensities for the hpG1 simulation compared to the more refined simulations hpG2 and hpG3, as can be seen both in the right and the left panel of figure 16 for 𝑥/𝐷 > 0.5. This generates an underestimation of the turbulence intensities with respect to experiments in the early phases of the jet development for 𝑥/𝐷 < 10. This behaviour can be found again in reference simulations. Dotted curves, representing coarser meshes, show lower rms levels than dashed curves, representing the finer mesh for each author. This is particularly evident for Shur and Gand. The decrease in axial/peak rms velocities for 𝑥/𝐷 < 10 can be directly associated with the longer potential core length for fine simulations [START_REF] Bogey | Grid sensitivity of flow field and noise of high-Reynolds-number jets computed by large-eddy simulation[END_REF][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]. Moreover, a lower resolution also leads to lower rms axial/peak velocities downstream of 𝑥/𝐷 = 10.

For 𝑥/𝐷 > 10 peak rms levels in the lipline match more closely the experiment when the mesh is adaptively refined, while axial rms levels for 𝑥/𝐷 > 10 suggest that statistics should be collected for a longer sampling period to provide smooth results. In fact, from 𝑥/𝐷 > 10 the flow varies at lower frequencies than the zone close to the exit nozzle. Moreover rms quantities need more time to converge than mean ones, and axial quantities cannot be azimuthally averaged. However, despite a higher intensity than GandG1/G2 and LorteauG1, the local peaks of axial rms for 𝑥/𝐷 > 10 are reproduced by the hpG3 simulation. Figures 17 and18 show the evolution of the mean and rms streamwise velocity at different axial downstream the nozzle exit. Looking at the figure 17, in particular 𝑥/𝐷 = 1, 5 locations, the shape of the mean velocity profile is better captured with mesh adaptation for the most significant zone 𝑟/𝐷 < 0.5, and a lower shear layer thickness of the jet is observed for hpG2 and hpG3 simulations. Moving downstream to 𝑥/𝐷 = 10, 15 profiles, similar observations to those made for the centerline profile in figure 15 can be made. hpG1 clearly underestimates the maximum velocity, due to a faster decay of the centerline velocity. hpG2 seems to provide results in close agreement with experimental results, while hpG3 simulation overestimates the velocity for 𝑟/𝐷 < 0.5 in both the locations. However, hpG3 profiles at 𝑥/𝐷 = 10 closely matches to GandG2 simulation, and has an intermediate value between GandG1 and GandG2 for 𝑥/𝐷 = 15. This confirms once again the trend of overestimation of the axial velocity when refining the mesh, already observed by other authors who performed numerical simulations on two different grid resolutions. From the profiles of rms velocity, we see that hpG2 and hpG3 provide lower levels of turbulent intensities with respect to hpG1 results, for 𝑥 = 1, 5, 10. Since the simulation does not present explicit turbulent structures at the nozzle exit, lower levels of rms velocity are expected in the vicinity of the nozzle. The higher levels of the G1 simulation, apparently closer to the experimental results than hpG2 and hpG3, are attributed to a lack of resolution, which fictitiously provides a better solution (as observed in the right panel of figure 16). The same behaviour is observed for 𝑥/𝐷 = 5 and 𝑥/𝐷 = 10, yet the hpG1 simulation shows a substantial overestimation of rms levels when getting closer to the axis. All these observations fall in line with Brès and Gand results, respectively displayed with green and magenta lines, comparing the coarse (dotted lines) and the fine (dashed lines) simulations. The rms velocity profile at 𝑥/𝐷 = 15 matches more closely the experiment as the mesh is refined, as can be observed also for Brès and Gand simulations. This emphasizes the ability of the hp-adaptation strategy to detect and refine the turbulent structures in the far wake and improve in turn the estimation of the turbulent fluctuations.

Overall good results, in line with numerical references, have been obtained with the present unstructured hp-adaptive strategy with a reasonable amount of degrees of freedom. The most refined hp-adapted mesh employed here counts a lower number of dofs than both the unstructured DG-p3 simulations by Lorteau [START_REF] Lorteau | Turbulent jet simulation using high-order DG methods for aeroacoustic analysis[END_REF] and the structured FV simulations by Gand, while achieving very close results.

Observations and comparisons with Brès and Shur in terms of number of dofs employed and achieved accuracy are quite tricky, because the former employed turbulent injection, changing the physics of the simulation, while the latter did not mesh the interior of the nozzle.

To conclude the assessment of the present simulations, the near field pressure fluctuations are analyzed. We directly record the pressure time history at two probes P1(𝑥/𝐷 = 2.62, 𝑟/𝐷 = 1.07) and P2(𝑥/𝐷 = 7.65, 𝑟/𝐷 = 1.78) downstream the nozzle exit for 36 equally spaced locations in the azimuthal direction. The Power Spectral Density is then computed from the pressure fluctuations signal and the resulting spectra are azimuthally averaged.

In figure 19 we compare hpG1, hpG2, hpG3 computations to the experiment and the coarser G1 simulation from Brès for which results are immediately available.

The cut-off Strouhal number, defined as 𝑆𝑡 𝑐𝑜 = 𝐷/(𝑐 𝑝𝑤 𝑣𝑜𝑙 1/3 𝑀 𝑎 ), where 𝑐 𝑝𝑤 is the number of cells per wavelength necessary to accurately propagate an acoustic wave, and 𝑀 𝑎 = 𝑈 𝑗 /𝑐 ∞ is the acoustical Mach number, depends on the resolution of the mesh in each sampled location. This means that more refined grids can resolve the PSD spectra up to increasing 𝑆𝑡 𝑐𝑜 numbers. This expected behaviour can actually be observed in both panels in figure 19.

Moving to the medium frequency range of the spectra, we notice a global overestimation of near field noise levels for the discretization hpG1 with respect to the experiment. While the discretization hpG2 provides spectra in better agreement with experimental results in the low-frequency range, an overestimation of the noise level is still present in the Strouhal range 1 < 𝑆𝑡 < 3. This overestimation is more marked for the probe located closer to the nozzle exit.

Using the finer hpG3 discretization, the shear-layer development is characterized by lower pressure spectra levels, in accordance to the lower turbulence intensities already pointed out in figure 16. This results in a very good prediction of the near-field noise levels in the medium frequency range, and a slight underestimation in the low-frequency range of the spectra. However, it should be noticed that some mismatch between adaptive numerical results and numerical/experimental results can be also caused by the relatively short duration of the adaptive simulations with respect to reference simulations. The lack of turbulence injection method may be the cause of the underestimation of the levels of the pressure spectra.

VI. Conclusions

In the present paper, an hp-adaptation strategy suitable to discontinuous Galerkin methods on hybrid prismatic/tetrahedral meshes has been proposed for solving turbulent fluid flow problems. An error estimator based on both the measure of the energy contained in the highest order polynomial modes and the jumps at the element interfaces, used in combination with a smoothness estimator, has been assessed in the framework of metric-based hp-adaptation for turbulent flows. A sensitivity study on the minimum wall spacing Δ𝑦 + for RANS has been carried out for a flat plate at Reynolds 5 • 10 6 , showing that increasing the order of accuracy of the DG approximation allows for relaxing the near-wall resolution. The value of the DES constant has been assessed for DG from a test case measuring the decay of isotropic homogeneous turbulence on grids composed of tetrahedra. It has been shown that the modulation of the DES constant does not impact significantly the enstrophy levels, however the choice of the polynomial degree has a significant impact and at least a fourth order approximation is required to obtain a good representation of turbulent structures in fully-developed turbulent regions. hp-adapted DG computations of the PPRIME nozzle at 𝑅𝑒 𝐷 = 10 6 have been carried out for a ZDES configuration, and the impact of increased resolution was analysed for three adapted meshes hpG1, hpG2, hpG3 presenting respectively 11.6, 18.9, 31.7 million dofs. Quantitative results compared fairly well to numerical references, obtaining overall close results to classical FV schemes on structured meshes and non-adaptive DG methods, with a reduced number of degrees of freedom. The improvements of the solution that the hp-adapted meshes provide with respect to initial RANS-adapted meshes are assessed on flow features as well as on the near pressure acoustic field.

Future work will include the acoustic analysis of the farfield propagated jet noise.

Fig. 2 RANS

 2 Fig. 2 RANS Flat plate at 𝑅𝑒 = 5 • 10 6 . The computational mesh G1 with stretching factor 𝛼 = 1.05 (left) and Mach contour on the whole domain of the flat plate (right)

Fig. 3 RANS

 3 Fig. 3 RANS Flat plate at 𝑅𝑒 = 5 • 10 6 . Error on the value of the integrated friction drag coefficient 𝐶 𝐷 𝑓 𝑟𝑖𝑐 along the flat plate.

Fig. 4 RANS

 4 Fig. 4 RANS Flat plate at 𝑅𝑒 = 5 • 10 6 . Error on the value of the friction coefficient 𝐶 𝑓 𝑟𝑖𝑐 at x=0.1 (left) and at x=0.97 (right) of the flat plate.

)

  𝑣 = -𝑉 0 • cos (𝑥/𝐿) • sin (𝑦/𝐿) • cos (𝑧/𝐿) 𝑤 = 0 representing large vortical structures, which break up into smaller structures and generate an energy cascade. The flow is initialized to be isothermal (𝑝/𝜌 = 𝑝 0 /𝜌 0 = 𝑅𝑇 0 ). The flow is computed at Reynolds number 𝑅𝑒 = 𝜌 0 𝑉 0 𝐿/𝜇 0 = 5000, and the unsteady simulation is performed for 𝑇 = 15𝑡 𝑐 , where 𝑡 𝑐 = 𝐿/𝑉 0 is the characteristic convective time.The pseudo-eddy viscosity field ν is initialized with the turbulent viscosity obtained from the Smagorinsky model

Fig. 5 ZDES

 5 Fig. 5 ZDES Taylor-Green vortex at 𝑅𝑒 = 5000. Computational mesh (left) for a 𝑝 = 3 simulation with resolution 72 3 , and respective iso-surface of 𝑄 = 5(𝑉 0 /𝐿) 2 at 𝑡 𝑐 = 10, coloured by the streamwise velocity component (right).

Fig. 6 ZDES

 6 Fig. 6 ZDES Taylor-Green vortex at 𝑅𝑒 5000. The 3D energy spectrum within inertial sub-range at time 𝑡 𝑐 = 10 compared to the Kolmogorov scaling 𝐸 (𝑘) ∝ 𝑘 -5/3 . p1, p2, p3, p4, p5 simulations with 𝐶 𝐷𝐸𝑆 = 0.65.

65 Fig. 7 ZDES

 657 Fig. 7 ZDES Taylor-Green vortex at 𝑅𝑒 = 5000. Volume-averaged enstrophy (left) and kinetic energy (middle) over the course of the simulation for a DGp3 discretization. The 3D energy spectrum within inertial sub-range at time 𝑡 𝑐 = 10 compared to the Kolmogorov scaling 𝐸 (𝑘) ∝ 𝑘 -5/3 (left). Different values of 𝐶 𝐷𝐸𝑆 .

65 Fig. 8 ZDES 65 Fig. 9 ZDES

 658659 Fig. 8 ZDES Taylor-Green vortex at 𝑅𝑒 = 5000. Volume-averaged enstrophy. Different values of 𝐶 𝐷𝐸𝑆 for three resolutions NDOFS = 48 3 , 96 3 , 144 3 .

Fig. 10

 10 Fig. 10 PPRIME. Slice of the computational domain for the hpG1 discretization.

Fig. 11

 11 Fig. 11 PPRIME nozzle. Zoom on the nozzle exit region and early jet plume for the three hp-adapted meshes, hpG1 at the left, hpG2 in the middle, hpG3 at the right.

Fig. 12

 12 Fig. 12 PPRIME nozzle. Iso-contours of temperature 𝑇/𝑇 ∞ in colored-scale and ( 𝑝 -𝑝 ∞ )/𝑝 ∞ in greyscale. Zoom of the hpG1 mesh (top), the hpG2 mesh (center) and the hpG3 mesh (bottom).

Fig. 13

 13 Fig. 13 PPRIME nozzle. Iso-contours of vorticity at 𝑥/𝐷 = 1, 2, 3 (from left to right) for hpG1(top), hpG2(middle), hpG3(bottom).

  15 and 16 corresponding to dotted lines. The shear layer is here laminar close to the nozzle exit, due to the lack of a turbulent injection technique to trigger the transition to turbulence within the nozzle.

Fig. 15

 15 Fig. 15 PPRIME nozzle. Mean streamwise velocity profiles on the jet axis 𝑟/𝐷 = 0 (left) and on the lipline 𝑟/𝐷 = 0.5 (right). Comparison between hp-adaptive and reference results.

Fig. 16

 16 Fig. 16 PPRIME nozzle. Root-mean square of the streamwise velocity profiles on the jet axis 𝑟/𝐷 = 0 (left) and on the lipline 𝑟/𝐷 = 0.5 (right). Comparison between hp-adaptive and reference results.

Fig. 17 PPRIME

 17 Fig. 17 PPRIME nozzle. Mean streamwise velocity profiles on different locations in the jet plume 𝑥/𝐷 = 1, 5, 10, 15 (from left to right). Comparison between hp-adaptive and experimental results.

Fig. 18 PPRIME

 18 Fig. 18 PPRIME nozzle. Root-mean square of the streamwise velocity profiles on different locations in the jet plume 𝑥/𝐷 = 1, 5, 10, 15 (from left to right).

Fig. 19

 19 Fig. 19 PPRIME nozzle. Power Spectral Density of the pressure signal recorded at locations P1(𝑥/𝐷 = 2.62, 𝑟/𝐷 = 1.07) at the left, and P2(𝑥/𝐷 = 7.65, 𝑟/𝐷 = 1.78) at the right.

  𝜖 𝐾 = 𝜖 SSED,𝐾 ,𝑛𝑜𝑟 𝑚 + 𝜖 JUMP,𝐾 ,𝑛𝑜𝑟 𝑚 (28) = 𝜖 SSED,𝐾min(𝜖 SSED ) T 𝐾 max(𝜖 SSED ) T 𝐾min(𝜖 SSED ) T 𝐾 + 𝜖 JUMP,𝐾min(𝜖 JUMP ) T 𝐾 max(𝜖 JUMP ) T 𝐾min(𝜖 JUMP ) T 𝐾

	normalization) before the coupling:
	).
	Both indicators are normalized by their respective maximum and minimum values over the whole domain T 𝐾 (min-max

  T ℎ,free,𝑖 = {𝐾 ∈ T free,𝑖 | 𝜖 𝐾 ,𝑖 > 𝜖 * 𝐾 ,𝑖+1 and [𝜎 𝐾 ,𝑖 ≤ 5/3 or (𝜎 𝐾 ,𝑖 > 5/3 and 𝑝 𝐾 ,𝑖 + 1 > 𝑝 max )] and 𝑝 𝐾 ,𝑖 ≠ 1} (35) 8) The polynomial degree of elements marked for 𝑝-adaptation is adapted as:

	𝑝 * 𝐾 ,𝑖+1 = max( 𝑝 𝐾 ,𝑖 + 1, 𝑝 𝑚𝑎𝑥 ),	∀𝐾 ∈ T 𝑝,𝑖

and [(𝜎 𝐾 ,𝑖 > 5/3 and 𝑝 𝐾 ,𝑖 + 1 ≤ 𝑝 max ) or 𝑝 𝐾 ,𝑖 = 1]} (34) 7) Elements free to be h-or p-adapted are marked for ℎ-adaptation:

Table 1 RANS

 1 Flat plate at 𝑅𝑒 = 5 • 10 6 .

	3.

Resolution of the employed grids G1 and G2. The 𝑦 + 𝐷𝐺 value needed to achieve 1% error on

  𝐶 𝐷 𝑓 𝑟𝑖𝑐 , and the achieved error on 𝐶 𝐷 𝑓 𝑟𝑖𝑐 at 𝑦 + 𝐷𝐺 ≃ 1 * for

the three discretizations (in the last two columns).

Fig. 14 PPRIME nozzle. Boundary layer velocity profile at

  this overestimation to the lack of an appropriate turbulent injection technique to trigger the transition to 𝑥/𝐷 = 0.04.

		0.60
		0.55					experiment hpG1-ZDES hpG2-ZDES hpG3-ZDES hpG1-RANS
		0.50
	r/D	0.45
		0.40
		0.35
		0.30	0.0	0.2	0.4	0.6	0.8	1.0
					u/U j

Comparison between hp-adaptive and reference results.

		#dofs (M)	spatial discretization	turbulence modelling	turbulent injection	simulated (𝐷/𝑈 𝑗 ) duration	mesh type
	Brès G1 [62]	16	FV-2o	LES-Vreman	yes	2000	unstruct hexa
	Brès G2 [62]	69	FV-2o	LES-Vreman	yes	1150	unstruct hexa
	Gand G1 [63]	48	FV-2o	ZDES mode 2	no	300	struct hexa
	Gand G2 [63]	154	FV-2o	ZDES mode 2	no	300	struct hexa
	Shur G1 [65]	8.4	FV-ho	RANS profile+ILES no	250	struct hexa
	Shur G2 [65]	23	FV-ho	RANS profile+ILES no	250	struct hexa
	Lorteau G1 [40] 78	DG-4o	LES-Smagorinsky	no	180	unstruct tetra
	Lorteau G2 [40] 165	FV-2o	LES-Smagorinsky	no	250	unstruct tetra
	hpG1 [present]	11.6	hp-DG	ZDES mode 1	no	150	unstruct tetra+prisms
	hpG2 [present]	18.9	hp-DG	ZDES mode 1	no	150	unstruct tetra+prisms
	hpG3 [present]	31.7	hp-DG	ZDES mode 1	no	150	unstruct tetra+prisms

Table 2 PPRIME nozzle. Comparison between parameters in simulations in the literature and present hp-adapted simulations.

 2

Acknowledgments

The authors would like to thank Dr Fabien Gand and Dr Mathieu Lorteau from ONERA for providing access to numerical data and fruitful discussions on jet simulations. Dr Fabio Naddei from ONERA is acknowledged as well for his valuable help in DG-adaptive strategies. This work is performed using HPC resources from GENCI-CINES (Grant 2021-A0102A11470) and internal ONERA and Airbus HPC resources.