Tanguy Messager
email: uy.messager@univ-brest.fr

-Mariusz Pyrz
email: mariusz.pyrz@eudil.fr

Pierre Chauchot
email: chauchot@lfremer.fr

Optimized laminations for submarine composite hulls

Keywords: cylinders, buckling, composite materials, laminations, optimization, genetic algorithm MOTS-CLES: cylindres, flambage, materiaux composites, stratifications, optimisation, algorithme genetique

This work deals with the development of a numerical tool for the design of filament wound submersible composite cylinders. A genetic algorithm coupled with analytical shell models allows stacking sequences to be optimized in order to increase the buckling pressures. The examples presented concern IFREMER's current studies and characteristic lamination patterns have been obtained. FEM calculations confirm the gains corresponding to the solutions obtained. The robustness of the optimized laminations is attentively studied. The peiformances of the optimization procedure are also discussed. RESUME: Ce travail concerne /'elaboration d'un outil numerique d'aide a Ia conception d 'enceintes so us-marines cylindriques composites realisees par enroulement filamentaire. Un algorithme genetique couple a des modeles analytiques des coques permet d'en optimiser Ia stratification a fin d 'accroitre Ia limite de flambage. Les applications presentees sont liees aux travaux actuels de l'IFREMER. Les calculs d'optimisation ont conduit a /'obtention d'un type de stratification caracteristique. Les calculs MEF confirment /'interet des solutions optimisees. Une attention particuliere est accordee a /'etude de Ia robustesse de ces solutions. De plus, les peifonnances de Ia procedure d' optimisation sont analysees.

Introduction

The use of composite materials in the manufacturing of submersible structures allows low weight to displacement ratios and long endurance for devices of limited energy carrying capability [YIN 92] [GRA 95]. Various potential applications, such as autonomous underwater vehicles and submarine exploration hulls, are currently developed by IFREMER with participation in European BRITE and MAST programs [DA V 96]. Figure 1 presents an application with a towed acoustic sonar. The cylindrical hulls used for such structures are composed of multilayered thin or medium-thick composite cylinders, obtained by the filament winding process [DVO 96], and closed by rigid metallic end-closures as shown in Figure 2. Experimental tests of cylindrical highly pressurized hulls (60 MPa service pressure) have been carried out by IFREMER. These tests, validated by numerical studies using FEM models, have confirmed the risk of failure due to buckling caused by the external hydrostatic pressure [DAY 96].

All cylinders investigated so far used orthotropic [ +55° N] (angles measured with respect to the cylinder axis) stacking sequences. This choice is due to the ratio between the circumferential and the axial load components resulting from the pressure [GAY 97]. Such a lamination pattern, classical for thin tubes submitted to internal hydrostatic pressure, does not take into account the instability risk.

Figure 1. Example of application

Figure 2. Cylindrical hull

This paper deals with an optimal design tool for submersible composite cylinders, allowing the search for lamination patterns that increase the buckling pressure Pa. The approach developed consists in coupling analytical models of cylindrical composite shell buckling with an optimization procedure based on a genetic algorithm directly manipulating the integer design parameters. The numerical applications presented, including specific manufacturer constraints, are connected with IFREMER's current studies. FEM models are applied to verify the stability limit increases corresponding to optimized stacking sequences. The performances of the optimization procedure are studied. Moreover, the robustness of the optimized solution patterns is also enhanced by analyzing the influences of the material characteristic values and the geometrical imperfections.

Cylinder modeling

The geometry of the cylindrical shell is characterized by the length L, its meanradius R and the wall-thickness h. As shown in figure 3, the u, v and w displacements are expressed in the (x,y,z) coordinate system related to the cylindrical mid-surface. Each k 1 h composite layer characterized by the filament winding angle noticed fA (measured in degrees with respect of the cylinder axis) is cross-ply, i.e. made up of equal amounts of fibers in the +Bk and -Bk directions. According to previous modeling studies [GUG 95] [PAP 98], the cylinders are assumed to be simply supported at their ends. The geometrical imperfections are the main source of buckling load reductions for both isotropic and composite thin-walled cylinders under hydrostatic pressure [TAB 97]. These geometric defects lead to considerable discrepancies (sometimes exceeding 50%) between experimental results and perfect cyJinder model predictions [BAR 83]. For several years, many authors have studied and developed analytical models of imperfect laminated cylinders subjected to buckling. Those formulations have been based on two main features [TEN 71]:

-The geometrical imperfection is modelled by an axisymetric, waving, modal deflection of the mid-surface of the cylinder. This assumption is an extension of previous analysis, concerning the buckling of isotropic imperfect cylinders.

-The geometrical defect effects are assumed to be of the second order. Only the nonlinear terms of deformation take into account this modal imperfection. This assumption leads to complex nonlinear formulation, then requiring specific and expensive solving procedures. A new formulation was developed for our study: by analogy to their real shapes [DVO 96], each ply is assumed to have a thickness defect modeled by an axisymetric modal waviness defined as:

w =a sin (mnxj L)
where a is the equivalent maximum imperfection amplitude related to each ply. As for the classical imperfection modeling [TEN 71], [BAR 83], the number of the longitudinal half waves m is related to m. The global geometrical imperfection of the cylinder is then assumed to be generated by all of these ply defects. The effects of the ply imperfections are introduced in the problem by correcting the stiffness coefficients of the laminate. The final linear form of the eigen-value problem of buckling is obtained using Galerkin's method. The formulation of these geometrical imperfection models, tests and comparisons with FEM results are detailed in [MES 99].

FE models

The finite element modeling has been performed in order to verify the critical load increases for optimized stacking sequences. A linear analysis of stability has been carried out using the Samcef FEM code for three types of composite laminated shell elements: hybrid (neglecting the shear transverse effects), Mindlin and volume elements. The rigid end-closures were modelled by rigid body elements leading to simply supported boundary conditions [GUG 95] [PAP 98]. Moreover, the hybrid shell element model allows the buckling analysis of geometrically imperfect cylinders. As for the most part of the standard FEM codes, the global geometrical imperfection is modeled by introducing a modal deviation in the coordinate values of the nodes. Independently, calculus has been performed using the Adina code and FE model, similar to the Samcef Mindlin one. Moreover, some comparisons have been realized applying a laminated shell element (neglecting the transverse shear effects) of the Solvia FEM code [PAP 98].

Optimization procedure

GAs principle

Among several computational techniques developed to deal with discrete engineering optimization problems, the promising applications of Genetic Algorithms (GAs) have been explored intensively in recent years. GAs are a class of computational models inspired by evolution, which may be used directly to solve non differentiable and discontinuous problems. They are particularly effective in the case of many locally optimal points. Numerous recent applications of GA to various engineering problems presented in the literature have shown the usefulness of this class of search procedures. The standard GA manipulates a fixed number of individuals, each corresponding to a potential solution of the optimization problem. A measure of performance, called fitness and corresponding to the objective function, is associated with each of them. The starting population is randomly created. The process of evolution is simulated using a set of biologically inspired operators (selection, crossover and mutation). Accordingly to evolutionary theories, the most suited elements, having the higher fitness values, are likely to survive and to produce next generations of "children" by recombining the features of "parents". Then, each iterative step creates a new generation of possible solutions. Finally, the GA sends back the best individual. This simulated evolution process is detailed for example in [GOL 94], [MIC 96] and [SEB 96].

Variable encoding and genetic operators

For a large majority of GAs applications in structural optimization, the design parameters are binary encoded and manipulated as bit strings called chromosomes [GOL 94]. This technique seems to be quite natural in many problems, for example in the case of topological optimization, where a value of binary gene can indicate the presence or absence of material in a given part of a structure [PYR 98]. Some authors [SEB 96] mention however a risk of processing of binary encoded variables due to epistasis phenomenon [MIC 96]. For the present study an "integer" version of GA, manipulating directly integer parameters, has been created. The whole arithmetical crossover operator is applied. This technique is defined as follows, where x 1 and i 2 are "parents", and .X' 1 and i' 2 are the corresponding "children": {

x

' 1 = I nt [ r x 1 + (I -r) x 2 ] with 2 2 1 /nt = mterger part r randomly generated x' = Int [r x + (I -r) x J j .
The selection procedure of the genetic algorith It)tStr $\1~5mqtM:tep}6Jluctive chances to the more fit individuals. The created GA applies the tournament procedure [SEB 96]: the bettern one of two individuals randomly selected is then chosen. This technique, inducing a constant selective pressure [SEB 96], can enhance efficiently GA search [Y AN 97]. The elitist selection technique is used and the best chromosome is automatically reproduced in the next generation. Moreover, the random uniform mutation [MIC 96] is applied.

Numerical applications

Cylinders considered

The applications presented concern IFREMER' s current studies carried out in connection with the manufacturer CNIM (Constructions Navales Industrielles de Ia Mediterannee). The structures considered are small scale cylinders having a 400 mm length and a 175 rnm internal diameter. The material is a carbon fiber reinforced epoxy resin. According to previous experimental and numerical studies [

Solutions including manufacturing constraints

The optimization problem is formulated as follows:

to maximize the objective function: F = Per

According to the manufacturer constraints, the filament winding angle values are subjected to the five following values (expressed in degrees):

8; E { 30 ,' 45 ,' 60 ,' 75 ,' 90}

For each thickness case investigated and for each analytical model used, 20 runs of GA program were performed. The following values of GA operator probabilities [GOL 94] have been applied in the optimization code: 80% for crossover and 0.8% for mutation. The numbers of individuals and numbers of generations are detailed in § 4.3 . The best stacking sequences (leading to greater Per values) obtained using the Sanders and the CT3 models coupled to the GA are presented in figure 5 and6, respectively. In these figures, the lamination plies are numbered from the inner to the outer surface of the cylinder. For each optimized stacking sequences, a table presents the critical pressures values calculated using different analytical and Samcef models. The shapes of buckling modes are similar to those presented in Figure 4. The increase of critical pressure, deducted from optimized stacking sequences (with respect to the [55N] initial solutions), is evaluated using the gain G and the nondimensionalized gain factor y defined as follows:

where Per ( [55 N]) corresponds to the limit of stability presented in table 3. The buckling pressure increase calculated using the FE models appears to be substantial.

As shown in Figures 5 and6, it is worth noting that the optimization procedure leads to typical, similar [90NI I a I 30N 2 I f3 I 90N 3 ] laminations, where a and f3 indicate plies of transition zones. However, the optimized solutions obtained using the Sanders-type model provide lower buckling pressure values than those obtained using the CT3 model. Thus, in spite of small differences on the buckling pressure values between the Sanders-type and the CT3 models (see Table 2 and Figures 5 and6), the optimization results appear to be sensitive to the shear transverse effects, especial1y for the thicker cylinders.

Additional buckling calculus, concerning the optimized stacking sequences presented in Figure 6 have also been performed using the Adina and Salvia FEM codes. The results, presented in Tables 4 and5, confirm significant buckling pressure increases. Moreover, we can notice that these critical loads are greater for all cylinders than the stability limits of In the case of cylinder number 6, an additional manufacturing requirement has been included in the optimization problem: five successive plies whose angle values are less than 75° are compulsorily followed by a 90° ply in order to compact the resin and so therefore to regularize the fiber ratio. Figure 7 presents the corresponding result, obtained using the CT3 model. In spite of this additional constraint, the [90N 1 I a I 30N 2 I f3 I 90N 3 ] pattern is preserved. Additional calculus have been carried out in order to improve the influence of the manufacturing constraints. The optimization of cylinders was then performed for ply angles changing from 0° to 90° by step of 1 o. The additional fiber regularization manufacturing constraint previously considered for the 32-ply cylinder, was not included. Figure 8 presents the best stacking sequences obtained (after 15 GA runs) using the CT3 model. The indicated y values are calculated using the Samcef Mindlin shell model. These optimized laminations are very similar to the ones obtained previously, always leading to [90Nl I a I ON2 I {3 I 90N 3 ] laminations. The y values are slightly greater than those obtained including manufacturing limitations. As noticed in § 4.2, the results corresponding to the Sanders-type model (not presented) also lead to lower buckling pressure values, especially for the thicker cylinders. Figure 9 presents the similarity between a [90Nt I ON 2 I 90N 3 ] lamination pattern and an 1-section beam. Thus, considering this analogy, we can conclude that the optimized laminations act as circumferential stiffeners.

Figure 9. Analogy with an /-section beam

The search for the worst stacking sequences, providing the lowest buckling pressures, was performed as well in order to verify how distant they are from the optimized solutions. For ply angles changing from 0° to 90° by step of I 0 , the worst solutions obtained always give [ON] laminations. Table I 0 presents the corresponding buckling pressure and the y values (then negative) calculated using the Samcef Mindlin shell model. These results complete the previous interpretation and the analogy presented in Figure 9: the worst laminations [ON] provide the lowest bending stiffness in the (y, z) plane and consequently lead to the lowest buckling pressure values.

cylinder n°

Performances of GA

This paragraph concerns the study of performances of the GA procedure for calculations performed in § 4.2.I and the reproducibility of the corresponding solutions presented in Figures 5 and6. The manufacturer angle constraints (cf. § 4.2.I ), limit the space design of the optimization problem to 5N different lamination solutions for a N-plies cylinders. Table 7 presents these cardinality values (i.e. the number of different possible solutions) for each cylinder thickness case. The number of individuals and generations used for the GA code and the corresponding CPU time (on a Sun Ultrasparc workstation) are also indicated. It is worth noting that the ratio between the number of buckling analyses realized during a GA run and the cardinality of the discrete optimization problem is about 0.1 % for the cylinder number 1 and less than 3.5 • 1 o-16 % for the cylinder number 6. 8 presents the average fitness (and the corresponding standard deviation) values of the best individuals obtained for these 20 optimization calculations performed using the CT3 model. The average values appear to be very close to the ones presented in Figure 6. The corresponding st~ndard deviation values indicated very small discrepancies. Moreover, all of the best laminations obtained for each of the 20 runs are always [90N 1 I a I 30N2 I {3 I 90N3] lamination patterns, having few numbers of angle value differences with respect to the solutions presented in Figure 6. The results deduced from the calculation using the Sanders-type model are similar. Thus, the GA procedure always provides similar solutions by exploring a very small number of possible configurations.

Sensitivity to material parameters

The experimental difficulties in characterizing properly the material of the structures under consideration [ DA V 96] induce doubts about the exactness of the orthotropic characteristic values presented in § 4.1. In order to evaluate numerically the influence of these uncertainties, 50 runs of GA, including the angle constraints (cf. § 4.2.1) have been carried out for each cylinder case. The CT3 shell model is used. Each of the nine orthotropic parameters is randomly and independently generated for every calculation within the following range of values (modulus in GPa):

102.9 The interval values are related to the study presented in [CHA 99] and correspond to variations of the fiber ratio and the longitudinal fiber Young's modulus. Figure 10 presents the average and the corresponding standard deviation values of the lamination angles for the 50 best solutions obtained for each cylinder case. As shown, for example, for the optimized lamination of the cylinder 1, only the plies number 3 and 9 vary (equal to 45, 60 or 75°). Thus, the [90NI I a I 30N 2 1 f3 I 90N 3 ] typical laminations are always preserved. Such a lamination pattern appears to be stable and robust in spite of uncertainties over material characteristics.

Conclusion

An optimal design tool for submarine composite, laminated cylinders exposed to instability risk has been developed using a genetic algorithm coupled with analytical buckling models. It allows the optimization of the stacking sequences in order to increase the buckling pressure values. Numerical applications concerning IFREMER's studies have been performed including the manufacturing constraints.

The optimized cylinder laminations are based on a typical pattern [90NI I a I 30N 2 I f3 I 90N 3 ] similar for all analyzed examples. The FE modeling has confirmed the significant increases (exceeding several tens of percent) of buckling loads of optimized stacking sequences with respect to classical [55N] laminations. Additional studies have confirmed that such laminations of [90NI I a I 30N2 I f3 I 90N 3 ] type are robust and preserved even if mechanical properties are uncertain and if structures are geometrically imperfect. Moreover, the numerical applications have shown the usefulness of the GA optimization procedure. Large discrete design spaces were processed and solved in acceptable CPU times.

An experimental study concerning the carbon reinforced epoxy resin and having the optimized laminations presented in this paper, is currently being performed in IFREMER laboratories. Some buckling tests on glass/epoxy cylinders with [55N] and optimized stacking sequences are also being performed on a special testing device described in [GIN 99]. The first results indicate a substantial increase of buckling pressure for the optimized laminations.

FigureFigure 4 .

 4 Figure 3. Geometry and coordinate system

cylinder

  6~ . . . . . . . . . . . . . . . . . . . . .

  Figure 8. Unconstrained optimized solutions

9Figure 10 .plies

 10 Figure 10. Mechanical characteristic influence

Figure 11 .

 11 Figure 11. Optimized solutions including geometrical imperfections

Table 1 .

 1 DA V 96] [PAP 98], the following orthotropic mechanical properties are considered (modulus in GPa): £ 1 =156~ £2=9.65; £ 3 =6.57~ G 12 =5.473; G 13 =2.8; G 23 =3.925; v 12 =0.27~ v 13 =0.34; v 23 =0.492. Each cross-ply is about 0.624 mm thick [PAP 98]. Six cylinders, thin-walled and half-thick, are investigated. Table1presents the corresponding numbers of plies and thicknesses.Table2shows values of critical pressure, calculated using both the analytical and FE models, for the initial cylinders having [55N] lamination pattern. As mentioned in § 2.1, buckling modes obtained from FE and analytical models always match one of the shapes presented in Figure4. For the thicker cylinders, the shear transverse effects appear to have a great influence on the buckling pressure calculated using FE models. The results corresponding to the Sanders-type model are always similar to the CT3 ones. However, these values are lower than the FE results for the thinner cylinders and higher for the thickest. Thickness of the cylinders

	cylinder number	1	2	3	4	5	6
	number of plies N	10	13	16	19	26	32
	thickness h (mm)	6.2	8.1	10.0	11.9	16.2	20.0

Table 2 .

 2 Buckling pressures of the [55Nl cylinders 4.2. Optimization results

Table 4 .

 4 [90M I eM] laminations studied in [PAP 98]. Results obtained using Adina

	cylinder Per (MPa) G (MPa) f(%)	cvlinder Per (MPa) G (MPa) y(%)
	1	30.0	8.7	40.8	1	29.1	I 0.1	53.2
	2	54.7	13.3	32.1	2	56.4	17.8	46.1
	3	89.6	21.4	31.4	3	95.6	29.8	45.3
	4	133.4	31.4	30.8	4	145.0	45.8	46.2

Table 5 .

 5 Results obtained using Salvia

Table 6 .

 6 Buckling pressure and y loss values for the [ON] cylinders

		1	2	3	4	5	6
	Per (MPa)	4.7	8.1	13 .0	19.3	39.6	63.4
	r (%)	-74.5	-77.0	-77.4	-76.9	-72.7	-69.4

Table 7 .

 7 Cardinality and population parameters of the GA

	cylinder n°	1	2	3	4	5	6
	Cardinality	9.8•10 6 1.2•1 0 9 1.5•10 11 1.9•10	13 1.5•10 18 2.3•10	22
	number of individuals	50	60	70	80	90	100
	number of generations	200	300	400	500	600	800
	CPU time (minutes)	9	18	31	46	64	98
	As mentioned in § 4.2.1, 20 successive runs of GA have been performed for each
	cylinder case. Table

Table 8 .

 8 The [90N 1 I a I 30N 2 I {3 I 90N:d laminations appear to benefit from a high reproducibility. Statistical results of the GA runs

	cylinder n°	1	2	3	4	5	6
	Fitness average	31.46	53.21	80.95	118.70 246.41 404.52
	standard deviation	0.14	0.18	0.66	0.95	1.02	1.07
	4.4. Robustness of the optimized solutions
	4.4.1.

Acknowledgments

The authors wish to thank the "Conseil Regional de Bretagne" and the MAS3-CT97 -0091 project of the European Community MAST program for their support.


Marta De La Llave Plata
email: marta.de_la_llave_plata@onera.fr

Vincent Couaillier

Marie-Claire Le Pape

On the use of a high-order discontinuous Galerkin method for DNS and LES of wall-bounded turbulence

Keywords: Discontinuous Galerkin method, large-eddy simulation, direct numerical simulation, wall-bounded turbulence

We assess the performance of a high-order discontinuous Galerkin (DG) approach to simulate wall-bounded turbulent flows for a range of Reynolds numbers. First, a plane channel flow configuration with turbulent heat transfer at Re τ = 640 is considered. Second, a detached flow configuration is investigated. This configuration represents the flow over periodically arranged hills. Four bulk Reynolds numbers Re b are considered: 2 800, 10 595, 19 000, and 37 000. For Re b = 2 800 direct numerical simulation (DNS) is used. Large-eddy simulations (LES) based on the WALE subgrid-scale model are carried out for the higher Re b . The simulation results are compared to reference data from CFD and experiment. hp-convergence analyses are performed which demonstrate the superior performance of increasing the polynomial order as compared to refining the mesh. It appears from this work that the use of a subgrid modelling approach together with local hp-adaptation could greatly improve the accuracy of the solution without penalising the computational cost of the simulation.

Introduction

The development of mathematical and numerical tools for the solution of the Navier-Stokes (N-S) equations in the turbulent flow regime is a topic of acknowledged practical importance to scientists and engineers. The pioneering work of Orszag and Patterson in 1972 [START_REF] Orszag | Numerical simulation of threedimensional homogeneous isotropic turbulence[END_REF], proved that it is possible to perform computer simulations of fully developed turbulent flows, called direct numerical simulations (DNS). DNS provides a deterministic time-dependent numerical solution of the fluid flow equations, in which all the scales of motion are accurately resolved. Unfortunately, DNS is computationally very expensive, and is restricted to low-to-moderate Reynolds numbers. This fact has motivated the scientific community to explore alternative approaches capable of reducing the number of degrees of freedom (DOFs) necessary to resolve the turbulence.

The large-eddy simulation (LES) approach provides a way of reducing the size of the numerical problem at hand by solving explicitly only for the largescale motions while modelling the effect of the small (or subgrid) scales on the resolved field. The separation of scales between large and small can be achieved either by applying a spatial filter to the N-S equations or by means of a projection operation [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows: An Introduction[END_REF][START_REF] Garnier | Large Eddy Simulation for Compressible Flows[END_REF].

The accurate and efficient simulation of turbulent flows imposes two main requirements on the discretization scheme. First, it must be able to provide very low levels of dissipation and dispersion. The aim is to resolve accurately over long integration periods the finest scales of turbulence that can be represented on the computational grid. Second, it must be highly parallelizable so that we can achieve high parallel performance on many-core architectures.

The high-order discontinuous Galerkin (DG) method [START_REF] Lesaint | On a finite element method for solving the neutron transport equation[END_REF][START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF][START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for scalar conservation laws II: general framework[END_REF][START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF]] meets these two requirements, i.e. high-order accuracy and excellent parallel capabilities [START_REF] Biswas | Parallel adaptive finite element methods for conservation laws[END_REF]. This is the reason why the scientific community is showing a growing interest in this type of approach [START_REF] Kroll | ADIGMA -A European Initiative on the Development of Adaptive Higher-Order Variational Methods for Aerospace Applications[END_REF][START_REF] Kroll | IDIHOM: Industrialization of High-Order Methods -A Top-Down Approach[END_REF].

The DG method is based on the variational projection of the N-S equations onto a hierarchy of polynomial basis functions. They combine some features of finite volume (FV) methods and finite element (FE) methods. It is characterised by being locally conservative and by providing high-order spectral accuracy not only in the interior of the computational domain but also at the physical boundaries, and on irregular meshes. Another attractive feature of the DG method is the possibility to enforce boundary conditions weakly, which has been shown by a number of authors to improve the quality of the solution [START_REF] Collis | Discontinuous Galerkin methods for turbulence simulation[END_REF][START_REF] Abbas | Weak versus strong no-slip boundary conditions for the NavierStokes equations[END_REF][START_REF] Nordström | Weak and strong wall boundary procedures and convergence to steady-state of the NavierStokes equations[END_REF]. Modal DG approximations also provide a natural framework for the implementation of adaption algorithms based on local p-refinement, by locally increasing (or decreasing) the polynomial order within the elements. Furthermore, the compacity of the DG scheme greatly simplifies the implementation of h-refinement techniques involving hanging nodes [START_REF] Houston | hp-adaptive discontinuous Galerkin finite/ element methods for first-order hyperbolic problems[END_REF].

The overall flexibility provided by the DG approach therefore makes this type of discretization a very appealing tool for the simulation of inhomogeneous high-Reynolds-number flows.

A small number of (well-resolved) DNS computations using DG have been reported in the literature. We can cite the work of Wei and Pollard [START_REF] Wei | Direct numerical simulation of compressible turbulent channel flows using the discontinuous Galerkin method[END_REF][START_REF] Wei | Interactions among pressure, density, vorticity and their gradients in compressible turbulent channel flows[END_REF] who applied the DG method to study the compressible turbulent channel flow for a range Mach numbers. In [START_REF] Chapelier | Inviscid and viscous simulations of the Taylor-green vortex flow using a modal discontinuous Galerkin approach[END_REF][START_REF] Chapelier | Evaluation of a high-order DG method for the DNS of turbulent flows[END_REF] Chapelier et al. have used the DG method to perform DNS of freely decaying and wall-bounded turbulent flows, highlighting the benefits of the DG approach over traditional numerical schemes in terms of number of DOFs. Diosady et al. [START_REF] Diosady | DNS of Flows over Periodic Hills using a Discontinuous Galerkin Spectral-Element Method[END_REF][START_REF] Garai | DNS of Flow in a Low-Pressure Turbine Cascade with Elevated Inflow Turbulence Using a Discontinuous-Galerkin Spectral-Element Method[END_REF] have also performed DNS simulations of a number of turbulent flow configurations. In [START_REF] Diosady | DNS of Flows over Periodic Hills using a Discontinuous Galerkin Spectral-Element Method[END_REF], the authors demonstrate the superiority of high-order over low-order methods at equal number of DOFs. In [START_REF] Bassi | On the development of an implicit high-order Discontinuous Galerkin method for DNS and implicit LES of turbulent flows[END_REF], Bassi et al. use an implicit high-order DG method to perform a DNS of the incompressible flow past a sphere at Re D = 1 000.

A greater number of under-resolved DNS, often also called no-model or im-plicit LES (ILES), computations can be found in the recent literature (see e.g. [START_REF] Bassi | On the development of an implicit high-order Discontinuous Galerkin method for DNS and implicit LES of turbulent flows[END_REF][START_REF] Gassner | On the accuracy of high-order discretizations for underresolved turbulence simulations[END_REF][START_REF] Uranga | Implicit Large-Eddy Simulation of transition to turbulence at low Reynolds numbers using a discontinuous Galerkin method[END_REF][START_REF] Carton De Wiart | Implicit LES of free and wall-bounded turbulent flows based on the discontinuous Galerkin/symmetric interior penalty method[END_REF][START_REF] Beck | High-order discontinuous Galerkin spectral element methods for transitional and turbulent flow simulations[END_REF][START_REF] Beck | High-order discontinuous Galerkin spectral element methods for transitional and turbulent flow simulations[END_REF][START_REF] Bolemann | High-order discontinuous Galerkin schemes for Large-Eddy simulations of moderate Reynolds number flows[END_REF][START_REF] Flad | Simulation of underresolved turbulent flows by adaptive filtering using the high order discontinuous Galerkin spectral element method[END_REF]). In the absence of an explicit subgrid-scale model, the philosophy behind ILES is to exploit the dissipation properties of the DG scheme to achieve a stable simulation. In order to reduce the numerical errors due to aliasing when very coarse meshes are used, some authors propose to use additional stabilisation techniques, such as polynomial de-aliasing (over-integration) or polynomial filtering (see e.g. [START_REF] Kirby | De-aliasing on non-uniform grids: algorithms and applications[END_REF][START_REF] Mengaldo | De-aliasing techniques for high-order spectral element methods on regular and irregular grids[END_REF][START_REF] Gassner | On the accuracy of high-order discretizations for underresolved turbulence simulations[END_REF]).

As regards to the application of the DG method to LES based on subgridscale modelling we can cite a number of works. Van Der Bos et al. [START_REF] Van Der Bos | Computational error-analysis of a discontinuous Galerkin discretization applied to Large-Eddy Simulation of homogeneous turbulence[END_REF] have performed LES of homogeneous isotropic turbulence (HIT) at low Reynolds number using the DG approach. They performed a parametric study in terms of the approximation order, the type of convective flux function and the value of the Smagorinsky constant. They found the best results when using the Smagorinsky model and a high order of accuracy. Marek et al. [START_REF] Marek | Large-eddy simulation of incompressible free round jet with discontinuous Galerkin method[END_REF] have performed an LES of an incompressible free round jet at Re D = 20 000 using an eddy-viscosity approach. Abbà et al. [START_REF] Abbà | Dynamic models for Large-Eddy Simulation of compressible flows with a high-order DG method[END_REF] have developed an anisotropic dynamic model in the context of a high-order DG method for compressible LES. To this end, the authors exploit the hierarchical nature of the FE basis to define the grid and test filters via projection operators. In [START_REF] Chapelier | Development of a multiscale LES model in the context of a modal discontinuous Galerkin method[END_REF], Chapelier et al. have proposed a variational multiscale approach based on a modal DG method and applied it to the Taylor-Green vortex configuration at Re = 3 000. Finally, Beck et al. [START_REF] Beck | On the Influence of Polynomial De-aliasing on Subgrid Scale Models[END_REF] have studied the influence of over-integration on the performance of subgridscale models. The authors demonstrate the benefits of using de-aliasing when a subgrid-scale closure is used at high polynomial degrees.

In this work we will focus on the use of a high-order modal DG approach to perform DNS and LES of two wall-bounded flow configurations for a range of Reynolds numbers. Firstly, under-resolved DNS simulations of a turbulent plane channel flow configuration with heat transfer at Re τ = 640 are performed. The adopted strategy consists in modifying the mesh size and the polynomial degree so that the number of DOFs remains constant. The results from the DG simulations are compared against the reference DNS data available in [START_REF]DNS database of Kawamura[END_REF].

A detached flow configuration, the periodic flow over a 2D hill, is also investigated. This configuration represents the flow over periodically arranged hills. This test case is of interest for the evaluation of high-order methods in the context of DNS and LES, due to the periodic boundary conditions and the 2D character of the geometry which significantly reduces the computational cost. A number of DG simulations of this configuration at Re b = 2800 and 10 595 have been reported in the recent literature (see e.g. [START_REF] Diosady | DNS of Flows over Periodic Hills using a Discontinuous Galerkin Spectral-Element Method[END_REF][START_REF] Beck | High-order discontinuous Galerkin spectral element methods for transitional and turbulent flow simulations[END_REF][START_REF] Abbà | Dynamic models for Large-Eddy Simulation of compressible flows with a high-order DG method[END_REF]). Here, we extend the study by also considering the higher Reynolds numbers Re b = 19 000 and 37 000.

Four different bulk Reynolds numbers, Re b based on the bulk velocity at the hill crest u b and the hill height h, are therefore considered: 2 800, 10 595, 19 000, and 37 000. The lowest Reynolds number case is computed using DNS and the three other cases are computed using LES based on the WALE subgrid-scale model [START_REF] Nicoud | Subgrid-Scale Stress Modelling Based on the Square of the Velocity Gradient Tensor[END_REF]. An additional no-model LES of the Re b = 10 595 configuration is performed and the results are compared to the corresponding model-based LES.

The DG solutions are compared to the CFD data provided by Breuer et al. [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Re[END_REF] (Re b = 2800 and Re b = 10, 595) using the second-order FV code LESOCC on a curvilinear grid and by Manhart et al. [START_REF] Manhart | Assessment of eddy resolving techniques for the flow over periodically arranged hills up to Re[END_REF] (Re b = 37, 000) using the secondorder FV code MGLET on a Cartesian mesh. All LES results are also compared to the experimental data provided by Rapp [40].

This work is aimed at demonstrating the capability of the DG method to produce accurate results for physically complex problems with a reduced number of DOFs. We also highlight the beneficial effect of employing an LES subgrid modelling approach (as opposed to a no-model approach) on coarse meshes when we move towards higher Reynolds numbers.

Governing equations

The equations governing the evolution of a Newtonian fluid are the N-S equations

∂U ∂t + ∇ • (f c -f v ) = S (1)
where, U = (ρ, ρV, ρE) T is the vector of conservative variables. The vectors f c and f v are the convective and viscous fluxes, respectively,

f c = ρV, ρV ⊗ V + p Ī, (ρE + p) V T f v = (0, τ, τ • V -q) T (2)
and S is a source term. In (2) ρ represents the density, E the total energy per unit mass, V is the velocity vector, and ⊗ denotes the tensor product. The total energy per unit volume is given by

ρE = p γ -1 + 1 2 ρ V 2 (3)
where the adiabatic index γ = cp cv is the ratio of specific heat capacities at constant pressure and constant volume conditions, respectively. The static pressure p is related to the static temperature T through the equation of state p = ρRT . Here the gas constant is defined as R = c p -c v . The stress tensor τ is related to the velocity gradients by the Stokes relations

τ = µ 1 2 ∇V + ∇V T - 2 3 ∇ • V Ī (4)
and the heat-flux vector is given by

q = -κ∇T (5)
where κ = cpµ P r is the thermal conductivity. The Prandtl number P r is assumed to be constant and equal to 0.71. The Sutherland law provides an expression for the dynamic viscosity µ as a function of the temperature.

LES modelling approach

In this section, the LES approach adopted in this work is briefly outlined. For a detailed description of the compressible LES formalism used in this research the reader is referred to [START_REF] Lesieur | Large-eddy simulations of turbulence[END_REF][START_REF] Chapelier | Développement et évaluation de la méthode de Galerkin discontinue pour la simulation des grandes échelles des écoulements turbulents[END_REF][START_REF] Chapelier | Development of a multiscale LES model in the context of a modal discontinuous Galerkin method[END_REF].

In what follows we use Favre (or density-weighted) filtering to define a resolved quantity φ in terms of the resolved (projected) conservative variables,

φ = ρφ ρ (6)
In order to model the effect of the unresolved scales on the resolved field, a subgrid-scale stress term τ sgs ij is added to the momentum equations. We make use of an eddy-viscosity assumption to express this term as a function of the known resolved field, namely,

τ sgs ij = ρν t 2 Sij - 2 3 Skk δ ij ; Sij = 1 2 ∂ ũi ∂x j + ∂ ũj ∂x i ( 7
)
where ρ is the resolved density, ν t is the eddy viscosity, and Sij denotes the resolved rate-of-strain tensor. In [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convection-dominated problems[END_REF], ũi denotes the i-th component of the velocity vector V. Based on the same hypothesis considered by Ducros et al. [START_REF] Ducros | LES of a spatially growing boundary layer over an adiabatic flat plate at low Mach number[END_REF] and Lesieur et al. [START_REF] Lesieur | Large-eddy simulations of turbulence[END_REF] a subgrid heat-flux vector term is also added to the energy equation,

Q i = ρc p ν t P r t ∂ T ∂x i (8)
in which P r t is the turbulent Prandtl number which is set to 0.6 . The LES simulations presented in this paper have been performed using the Wall-Adapting Local Eddy-viscosity (WALE) model proposed by Nicoud and Ducros [37]. This model is based on a tensor invariant and is able to represent the proper scaling at the wall, ν t = O(y 3 ). In the WALE approach, the eddy viscosity ν t is expressed as

ν t = (C W ∆) 2 Sd ij Sd ij 3 2 Sij Sij 5 2 + Sd ij Sd ij 5 4 (9)
where Sd ij is defined as

Sd ij = 1 2 g2 ij + g2 ji - 1 3 δ ij g2 kk ( 10
)
and it represents the symmetric part of the tensor g2 ij = gik gkj , where gij = ∂ ũi /∂x j . The filter width ∆ is defined in terms of the DG discretization parameters h and p as ∆ = h/(p + 1), and the model constant takes a value of 0.55 ≤ C W ≤ 0.6 as recommended in [START_REF] Nicoud | Subgrid-Scale Stress Modelling Based on the Square of the Velocity Gradient Tensor[END_REF].

The discontinuous Galerkin method

The simulation results presented here have been performed using the highorder DG solver Aghora [START_REF] Renac | Aghora: A High-Order DG Solver for Turbulent Flow Simulations, IDI-HOM: Industrialisation of High-Order Methods -A Top Down Approach[END_REF][START_REF] Chapelier | Evaluation of a high-order DG method for the DNS of turbulent flows[END_REF] developed at Onera. This solver is designed to solve the full set of compressible N-S equations on unstructured grids.

The DG discretization used in this work is based on a modal approach that relies on the use of a hierarchy of orthogonal polynomial functions as basis for the Galerkin projection. The modal basis used is that proposed in [START_REF] Bassi | On the flexibility of agglomeration based physical space DG discretizations[END_REF] by Bassi et al. The discrete orthogonal polynomial spaces are directly computed in physical space. The methodology consists in defining a starting set of monomial basis functions in each (arbitrarily shaped) element and applying a modified Gram-Schmidt orthogonalization procedure. This ensures the orthogonality of the basis, and thereby the diagonality of the mass matrix. Numerical integration is efficiently performed by means of a Gaussian quadrature. The time integration is performed by using a strong stability preserving third-order Runge-Kutta scheme [START_REF] Gottlieb | Strong stability-preserving high-order time discretization methods[END_REF]. The DG method implemented in the Aghora solver is briefly outlined below.

We start by defining a shape-regular partition of the domain Ω, into N nonoverlapping and non-empty cells κ of characteristic size h.We also define the sets E i and E b of interior and boundary faces in Ω h , such that

E h = E i ∪ E b . Let V p h = {φ ∈ L 2 (Ω h ) : φ| κ ∈ P p (κ)
, ∀κ ∈ Ω h } be the functional space of piecewise polynomials of degree at most p, and (φ 1 κ , . . . , φ Np κ ) ∈ P p (κ) a hierarchical and orthonormal modal basis of V p h , of dimension N p , confined to κ. The solution in each element is thus expressed as

u h (x, t) = Np l=1 φ l κ (x)U l κ (t), ∀x ∈ κ, κ ∈ Ω h , ∀t ≥ 0, (11)
in which the polynomial coefficients (U l κ ) 1≤l≤Np represent the DOFs of the discrete problem in element κ. The semi-discrete variational form of system of equations (1) thus reads: find

u h in V p h such that ∀ φ h ∈ V p h we have ˆΩh φ h ∂ t u h dV + L c (u h , φ h ) + L v (u h , φ h ) = ˆΩh φ h S h (12)
In Eqn. ( 12) L c and L v represent the variational projection of the convective and viscous terms, respectively, onto the functional space V p h . Similarly, the right-hand-side of Eqn. ( 12) is the variational projection of the source term S onto V p h . We now introduce the following notation: for a given interface e in E i we define the average operator {u} = (u + + u -)/2 and the jump operator [[u]] = u + n + -u -n -, where u + and u -are the traces of the variable u at the interface between elements κ + and κ -(see Fig. 1). The DG discretization of the convective terms then reads

L c (u h , φ h ) = - ˆΩh f c (u h ) • ∇ h φ h dV + ˆEi [[φ h ]]h c (u + h , u - h , n)dS + ˆEb φ + h f c u b (u + h , n) • ndS (13)
The numerical flux h c is chosen so that it is consistent and compact. For the simulations presented in this paper the local Lax-Friedrichs (LLF) flux has been employed.

h c (u + h , u - h , n) = {f c (u h )} + α [[u h ]] 2 (14) (15)
α = max{ρ s ∂(f c (u) • n) ∂u : u = u ± h } (16)
in which ∂ ∂u (f c (u) • n) is the Jacobian matrix of the convective fluxes in the n direction and ρ s is the spectral radius of the Jacobian matrix. The discretization of the viscous terms is performed using the symmetric interior penalty (SIP) method [START_REF] Arnold | Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems[END_REF], described in [START_REF] Hartmann | An optimal order interior penalty discontinuous Galerkin discretization of the compressible Navier-Stokes equations[END_REF],

L v (u h , φ h ) = ˆΩh f v (u h , ∇ h u h ) • ∇ h φ h dV - ˆEi [[φ h ]] f v u h , ∇ h u h • ndS - ˆEb φ + h f v u b (u + h , n), ∇u b (u + h , n) • ndS - ˆEi [[u h ]] G T u h ∇ h φ h • ndS - ˆEb (u + h -u b (u + h , n)) G T u b ∇ h φ + h • ndS + ˆEi η IP [[u h ]][[φ h ]]dS + ˆEb η IP u + h -u b (u + h , n) φ + h dS (17)
where G = ∂f v /∂ ∇ h u h is the so-called homogeneity tensor. The penalty parameter η IP has to be chosen sufficiently large to ensure the coercivity of the bilinear form and thus the numerical stability of the simulation [START_REF] Arnold | Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems[END_REF].

Turbulent heat transfer in a plane channel at Re τ = 640

Following the work of Abe et al. [START_REF] Abe | Surface heat-flux fluctuations in a turbulent channel flow up to Re τ = 1020 with P r = 0.025 and 0.7,1[END_REF] we consider a fully developed turbulent flow in a plane channel with a passive scalar field representing the temperature in an incompressible flow. The size of the computational domain is 12.8h × 2h × 6.4h, where h is half the height of the channel. Periodicity and homogeneity of the flow is assumed in the streamwise and spanwise directions. To study the heat transfer in the channel, the authors consider an additional equation governing the transport of the temperature

∂ρθ ∂t + ∇ • (ρθu) -∇ • κ c p ∇θ = f θ (18)
The temperature θ thus behaves as a passive scalar and f θ is a source term that ensures a constant heat-flux at the wall Q w . Among the different simulations reported in [START_REF] Abe | Surface heat-flux fluctuations in a turbulent channel flow up to Re τ = 1020 with P r = 0.025 and 0.7,1[END_REF], we focus on the DNS at Reynolds number Re τ = 640 and Prandtl number Pr = 0.71, performed using an incompressible solver on a grid composed of 1024 × 256 × 1024 nodes.

We also define the bulk Reynolds number based on the bulk velocity u b as

Re b = u b 2h ν = 24 428 u b = 1 2h ˆ2h 0 u 1 (y)dy
A forcing term is introduced in the x-momentum equation, -(∇ x p) o , and in the energy equation, -u(∇ x p) o . The pressure gradient (∇ x p) o is considered spatially constant and ensures a target mass flow rate ṁ0 corresponding to the specified Re b . Thus, the value of the Re τ is an outcome of the simulation that will allow us to assess the accuracy of the employed numerical approach. Based on the algorithm proposed by Deschamps [START_REF] Deschamps | Simulation numérique de la turbulence inhomogène incompressible dans un écoulement de canal plan, Note technique ONERA[END_REF] in the incompressible case, and extended to the compressible case by Lenormand et al. [START_REF] Lenormand | Contribution à la simulation des grandes échelles d'écoulements turbulents compressibles en milieux confinés[END_REF][START_REF] Lenormand | Large-eddy simulation of compressible channel flow at moderate Reynolds number[END_REF] the pressure gradient is computed at each time step as

(∇ x p) n+1 o = (∇ x p) n o + 1 L y L z ∆t ṁ0 -2 ṁn + ṁn-1 ( 19
)
where ṁn is the mass flow rate computed at time step n, and L y and L z represent the dimensions of the computational domain in the normal and spanwise directions, respectively.

To maintain consistency with the work presented in [START_REF] Abe | Surface heat-flux fluctuations in a turbulent channel flow up to Re τ = 1020 with P r = 0.025 and 0.7,1[END_REF], Eqn. ( 18) is solved at each time-step by using the DG approach presented in Sec. [START_REF] Lesaint | On a finite element method for solving the neutron transport equation[END_REF]. By doing so we obtain a scalar field devoid of compressibility effects and therefore a truthful comparison with the DNS results reported in [START_REF] Abe | Surface heat-flux fluctuations in a turbulent channel flow up to Re τ = 1020 with P r = 0.025 and 0.7,1[END_REF] can be made. The source term f θ is defined in terms of the imposed wall heat-flux1 Q w as

f θ = -ρu 1 (∇ x θ) o ; (∇ x θ) o = - Q w c p (ρu) b h (20)

Computational details

Under-resolved DNS computations of the configuration described above have been performed using the Aghora DG solver. Three different levels of hprefinement have been considered in this study, for which the number of DOFs is the same and equal to approximately 4.2 million (see Table 1). This resolution is four times coarser in each space direction than that used in the reference DNS. To ensure the stability of the simulation at this coarse resolution, overintegration is used as de-aliasing technique by considering an extra-quadrature point in each spatial direction. For a pth-order simulation we thus have (p + 1) 3 DOFs in each element, and (p + 2) 3 quadrature nodes.

The resolution at the wall is evaluated in wall units, i.e. ∆y + = Re τ (∆y/h), where ∆y denotes the height of the first cell at the wall, and Re τ = u τ h/ν is the friction Reynolds number. The friction velocity is given by u τ = τ w /ρ w , being τ w the wall-shear stress and ρ w the density at the wall. Note that in the DG approach, an estimate of ∆y + can be obtained by dividing the cell height by the number of DOFs within the element. For a polynomial approximation of degree p we therefore define y + = ∆y + /(p+ 1)/2. For the simulations presented in Sec. 5.2 this leads to the following values : y + < 2, ∆x + < 34, and ∆z + < 17. These values are well within the recommendations given by Piomelli and Chasnov [START_REF] Piomelli | Large-Eddy Simulations: Theory and Applications[END_REF] for wall-resolving LES: y + < 2, ∆x + = 50 -150, and ∆z + ≈ 15 -40.

The normalised values of the time step used in the different computations are compiled in Table 2. The CPU time per convective-time unit t c = 12.8h/u b required by each simulation is also provided in this table. Note that these times have been normalised by that of the cheapest simulation DG -p2.

Simulation results

Figures 2 and3 show a comparison between the profiles of mean velocity and velocity correlations obtained from the DG computations and the reference DNS data provided in [START_REF]DNS database of Kawamura[END_REF]. It is important to bear in mind that the value of u τ used to normalised the DG profiles is determined from the value of Re τ yielded by the simulation (see Table 2). This allows us to assess the quality of our results in an objective manner.

We clearly see that the best accuracy is obtained when using a fifth-order approximation on a coarse mesh. These results highlight the benefit of increasing the polynomial order with respect to reducing the mesh size, as regards the accuracy of the solution. This point has also been demonstrated in previous work [START_REF] Chapelier | Inviscid and viscous simulations of the Taylor-green vortex flow using a modal discontinuous Galerkin approach[END_REF][START_REF] Chapelier | Développement et évaluation de la méthode de Galerkin discontinue pour la simulation des grandes échelles des écoulements turbulents[END_REF][START_REF] Chapelier | Evaluation of a high-order DG method for the DNS of turbulent flows[END_REF].

Important discrepancies are found for the DG-p2 computation with respect to the two other computations. These are probably due to the fact that the corresponding value of Re τ is in this case significantly lower than the reference value of 640, leading to an overestimation of the streamwise velocity. It is worth noting that the slight oscillations observed on the profiles of the root-meansquare (rms) of the velocity fluctuations are due to the discontinuous character of the DG solution. Indeed, we have found that the wavelength associated with these oscillations corresponds to the distance between elements.

As regards the solution to the passive scalar equation, only the DG -p4 discretization has been considered. The reason is the very slow statistical convergence of this quantity, specially when using a compressible solver, which leads to a high computational cost. Figure 4 shows the profiles of the mean scalar and the correlations of the scalar and velocity fluctuations. These data have been normalised by u τ and the friction temperature θ τ = Q w /ρ w c p u τ . We can observe the very good agreement found with the reference data. This is despite the great reduction in the number of DOFs with respect to the reference DNS.

Turbulent flow over a 2D periodic hill

In this section, we consider a fully developed turbulent flow over a 2D periodic hill. This test case was initially proposed in an ERCOFTAC/IAHR workshop in 1995. A detailed description of the geometrical configuration of the channel as well as of the flow conditions can be found on the ERCOFTAC QNET-CFD website [START_REF]The ERCOFTAC QNET-CFD Wiki Home Page: qnet-ercoftac[END_REF].

This configuration represents the flow over periodically arranged hills. The flow presents separation from the curved surface, reattachment, and post reattachment recovery.

The dimensions of the channel are as follows : L x = 9h in the streamwise direction, L y = 3.035h in the vertical direction, and L z = 4.5h in the spanwise direction, where h is the hill height at the crest. Periodicity of the flow is assumed in the streamwise and spanwise directions, and isothermal no-slip boundary conditions are imposed on the lower and upper walls. The flow is driven by a mean pressure gradient in the streamwise direction. This body force ∇ x p is implemented as a forcing term in the x-momentum and energy equations, as explained in Sec. [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF].

The Mach number is set to M 0 = 0.1 so that the flow regime is quasiincompressible. The bulk Reynolds number Re b is defined based on the bulk velocity u b and on the hill height, h, namely,

Re b = u b h ν u b = 1 2.035 ˆ3.035h h u 1 (y)dy (21)
Four different values of Re b have been considered: 2 800, 10 595, 19 000, and 37 000. As already mentioned in Sec. 1, the case at Re b = 2 800 is considered as a DNS. For the higher values of Re b we employ an LES approach based on the WALE model [START_REF] Nicoud | Subgrid-Scale Stress Modelling Based on the Square of the Velocity Gradient Tensor[END_REF]. An additional no-model LES computation of the Re b = 10 595 configuration is also performed.

The following sections provide a summary of the computational parameters and the grids used in these simulations, as well as a summary of the main results obtained.

Computational details

Two levels of mesh refinement have been used in this study, shown in Fig. 5. A coarse fourth-order mesh composed of 8 192 elements and a twice finer version of this mesh composed of 65 536 elements. Both meshes were generated using the 3D finite element mesh generator Gmsh [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF].The details of the high-order grids used for the DG simulations as well as those used in the reference FV computations in [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Re[END_REF] (LESOCC code) and [START_REF] Manhart | Assessment of eddy resolving techniques for the flow over periodically arranged hills up to Re[END_REF] (MGLET code) are provided in Table 3. The polynomial degree of the simulations is set to p =3, 4 up to 5 for the DNS on the coarse grid, and p =3 and 4 for all other simulations.

As for the channel flow computations presented in the Sec. 5, the resolution at the wall is evaluated in terms of the normalised values ∆y + = ∆y u τ /ν. In this case, the friction velocity, u τ , is measured at the point of maximum wall-shear stress, which is located on the windward side of the hill as we will see in Sec. 6.2.1. Figure 6 shows the distribution of wall distances y + and grid spacings (∆x + , ∆z + ) along the channel. The values of y + at the point of maximum wall-shear stress are compiled in Table 4. We observe that for the fifth-order DNS computation on the fine mesh a value of y + < 1 is found which is sufficient for a wall-resolved simulation. On the coarser meshes this value was close to 2 (not included in Table 4). For the higher-Reynolds-number cases the values of y + are above the recommended value for a wall-resolving LES [START_REF] Piomelli | Large-Eddy Simulations: Theory and Applications[END_REF], although in line with the values reported in [START_REF] Manhart | Assessment of eddy resolving techniques for the flow over periodically arranged hills up to Re[END_REF] for the MGLET simulations at Re b = 10 595 and Re b = 37 000. As pointed out by the author, the distance over which this maximum wall stress is acting on the flow is very short (≈ 0.5h), and at all other positions the wall resolution can be considered as sufficient, as can be seen in Fig. 6.

The normalised time step used in the different computations is of the order of 10 -4 . It is worth noting that in this case the severe restriction on the time step arises from the convective time scale of the acoustic waves due to the low value of the Mach number (M 0 = 0.1), and not from the viscous time scale associated with the penalty term in the SIP method.

The solution is advanced in time until a statistically steady state is reached. From this point, the flow statistics are gathered over a sufficient number of convective times t c = L x /u b . The time-averaged data is also averaged in the spanwise direction to achieve more rapidly statistical convergence.

Simulation results

The results from the different DG simulations are analysed in detailed in the following sections. To this end, the profiles of mean streamwise and vertical velocity, as well as the turbulent stresses are compared against the available reference data at different stations along the streamwise direction. The locations 4). The streamlines of the time-averaged velocity field for Re b =2 800 up to 37 000 are depicted in Fig. 7. As expected, the reattachment point moves upstream with increasing Reynolds number, reducing the size of the recirculation bubble. The position of the reattachment point at the lowest Reynolds number is x R /h ≈ 5.4 which is in agreement with the value reported in [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Re[END_REF]. For Re b =10 595 the reattachment point is located at x R /h ≈ 3.9, which is slightly different from the experimental value (x R /h ≈ 4.21) [START_REF] Rapp | Flow over periodic hills -an experimental study[END_REF] and from the reference LESOCC simulation (x R /h ≈ 4.69). Finally, for the highest Reynolds number Re b =37 000 the reattachment point is located at x R /h ≈ 3.2 not far from the value of x R /h ≈ 3.68 reported by Chaouat and Schiestel [START_REF] Chaouat | Hybrid RANS-LES Simulations of the Turbulent Flow over Periodic Hills at High Reynolds Number[END_REF] in their hybrid RANS/LES simulation on a grid composed of 1 M points.

The profiles of time-averaged wall-shear stress at the bottom wall for Re b =2 800 are depicted in Fig. 8(a). We observe that all discretizations, except for the DG -p3 simulation on the coarse mesh, are able to reproduce the peak of maximum wall-shear stress. The strong oscillations near the separation point exhibited by the DG -p3 and DG -p4 coarse-mesh simulations are clearly due to a lack of resolution in this area. Indeed, this oscillating behaviour disappears when the polynomial degree is increased. This is shown by the results obtained from the DG-p5 coarse-mesh simulation, which closely match the data obtained from the best resolved simulations for a reduced number of degrees of freedom (1.77 M in the DG -p5 with respect to 4.19 M and 8.19 M in the DG -p3 and DG -p4 fine-mesh simulations, respectively). This outcome emphasises the greater advantage of using local hp-adaptation for flow configurations featuring highly localised phenomena.

In Fig. 8(b) we compare the profiles of wall-shear stress for all four Reynolds numbers. The trends are in agreement with the results reported in [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Re[END_REF]. Note, that this cost is normalised by that of the cheapest DNS computation, i.e. DG -p3 on the coarse mesh.

Figures 10 and11 show the profiles of averaged velocity and turbulent stresses obtained from the DG-DNS at the different levels of resolution considered (see Table 5). They are compared to the DNS data from the LESOCC code. It can be observed that all DG simulations are able to reproduce the main features of the flow (separation, recirculation and reattachment). Naturally, the best match with the reference DNS is obtained for the DG -p4 simulation on the fine grid (8.19 M DOFs), as we reach numerical convergence in terms of number of degrees of freedom.

However, if we look carefully at the solution from the DG -p5 coarse-mesh simulation we observe that the match with the DG -p4 fine-mesh simulation and the reference data is remarkably good. This is taking into account the great reduction in the number of DOFs with respect to the two other simulations (8.19 M and 13.1 M, respectively). It is also worth noting that the computational cost of the DG -p5 is significantly lower than that of the DG -p4 simulation and only slightly higher than that of the DG -p3 simulation on the fine mesh (see Table 5).

As already discussed in the previous section, the DG -p3 and DG -p4 coarse-mesh simulations exhibit an oscillatory behaviour in the vicinity of the separation point (see Fig. 10).These oscillations disappear when the polynomial order is increased to p = 5. We can also observe the presence of spurious overshoots in the profiles of u u at practically all stations. The locations of these peaks actually coincide with the grid element interfaces and are a consequence of the discontinuous character of the DG solutions at this coarse resolution. These jumps are smoothed out when the polynomial degree is increased, as can be observed from the u u profiles in the DG -p5 simulation.

We can conclude from this analysis that an optimal performance of the DG method could be obtained by appropriate local hp-adaptation. This would allow us to greatly reduce the cost of the simulation while keeping the same level of accuracy as that provided by uniform h and/or p-refinement.

LES results at Re b = 10 595

Figures 12 and13 show the profiles of mean velocity and turbulent stresses obtained from the DG -p3 fine-mesh simulations using the no-model and the WALE approaches. They are compared to the LES data produced by the LESOCC code as well as to the experimental data from Rapp [START_REF] Ch | Experimentelle Studie der turbulenten Strömung über periodische Hügel[END_REF].

It can be seen from these plots that at this relatively low Reynolds number the results yielded by the no-model and the WALE approaches are not fundamentally different. The WALE approach, however, seems to provide slightly better results in the region close behind the separation point (x = 0.5h). This is specially the case for the mean vertical velocity profile as well as for the turbulent stresses, for which the no-model simulation presents a more marked overshoot near the wall. The profiles obtained using the WALE approach appear to be closer to the reference data. These comments also hold for the area located in the middle of the separation zone (corresponding to x = 2h), although the differences become gradually less significant as we exit the recirculation zone (x = 4h). Finally, in the area well behind the recirculation bubble (x = 6h) both simulations provide similar results. However, the WALE simulation appears to represent with higher fidelity the mean vertical velocity profile which is closer to the experiment, though significantly far from the simulation.

It is also worth noting the smoothing effect that the addition of the subgrid model has on u u . Indeed, at all stations, we observe that the unphysical peaks generated by the DG method at the element interfaces are considerably damped when the WALE approach is used.

Figures 9(a Figures 14 and15 show the profiles of mean velocity and turbulent stresses obtained from the DG -p3 fine-mesh simulation using the WALE approach. They are compared to the experimental data from Rapp [START_REF] Ch | Experimentelle Studie der turbulenten Strömung über periodische Hügel[END_REF].

In view of the results from the previous section, only the WALE approach has been considered at this higher Reynolds number. In particular, due to the fact that the raw DG scheme did not appear to provide a sufficient level of dissipation to achieve a stable simulation. This was also the case for the Re b = 37 000 simulation.

An overall good agreement is found between the results from this simulation and the experimental data. The main differences are found in the vicinity of the separation point. At this location the profiles of u v from the DG simulation present a strong overshoot near the wall. This result might suggest a need to improve the resolution in this area. In the recirculation region (x = 2h) however, the numerical solution is in relatively good agreement with the experimental data. Right after the recirculation bubble (x = 4h and x = 6h) the match between the mean vertical velocity profiles deteriorates, while at the final stations (x = 8h) the agreement between the DG profiles and the experiment improves considerably.

It is difficult discern, however, whether the discrepancies observed at some stations in the channel are due to numerical errors, or to measurement uncertainties during the experiment. Further DG simulations at this higher Reynolds number, ideally using a locally hp-adapted grid, could help understand the nature of these differences. 6.6. LES at Re b = 37 000 Figures 16 and17 show the profiles of mean velocity and turbulent stresses obtained from the DG -p3 fine-mesh simulation using the WALE approach. They are compared to the LES data produced by the MGLET code [START_REF] Manhart | Assessment of eddy resolving techniques for the flow over periodically arranged hills up to Re[END_REF] as well as to the experimental data from Rapp [START_REF] Ch | Experimentelle Studie der turbulenten Strömung über periodische Hügel[END_REF].

An overall good agreement is found between the mean streamwise velocity profiles and the reference data. Some discrepancies between the experimental data and the DG solution are found however for the mean vertical velocity profiles, which appear in general to be closer to the reference CFD data. The match between the turbulent stresses u u and u v from the DG simulation and the experiment is relatively good. Note that the overshoot that appears near the bottom wall at x = 0.5h in both profiles is also observed in the MGLET profiles as well as in the numerical results reported by Chaouat and Schiestel in [START_REF] Chaouat | Hybrid RANS-LES Simulations of the Turbulent Flow over Periodic Hills at High Reynolds Number[END_REF].

Conclusions

The accuracy of a modal DG approach to simulate wall-bounded turbulent flows for a range of Reynolds numbers has been assessed in this work.

Firstly, under-resolved DNS computations of a plane channel flow configuration with turbulent heat transfer at Re τ = 640 have been carried out. The DG simulations involve four times less DOFs in each space direction than the reference DNS simulation of Abe et al. [START_REF] Abe | Surface heat-flux fluctuations in a turbulent channel flow up to Re τ = 1020 with P r = 0.025 and 0.7,1[END_REF]. In the absence of a subgrid model, a simple de-aliasing of the solution via overintegration is sufficient to achieve a stable simulation at this relatively low Reynolds number. The hp-convergence analyses clearly demonstrate the greater advantage of using high-order approximations on very coarse meshes, as opposed to mesh refinement in combination with low-order schemes. We actually find the best match with the reference DNS data for the fifth-order simulation on the coarsest mesh.

In a second stage, the ability of the DG method to capture the Reynoldsnumber dependent features in a detached flow configuration has been investigated for Reynolds numbers ranging from Re b = 2 800 up to 37 000. As for the channel flow configuration, hp-convergence studies have been carried out for the DNS at Re b = 2 800. Once more, the outcome from this analysis demonstrates the superior performance of increasing the polynomial order as compared to refining the mesh, specially in the area close to separation. A very good agreement with the reference simulation of Breuer et al. [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Re[END_REF] is actually obtained by considering a sixth-order discretization (p = 5) on a very coarse mesh, involving only 1.77 M DOFs with respect to the 13.1 M DOFs used in the reference FV computation. This result emphasises the great advantage of using local hpadaptation for flow configurations featuring highly localised phenomena, as is the case here.

For Re b = 10 595, 19 000, and 37 000 fourth-order LES (p = 3) based on the WALE approach have been carried out on a grid composed of 65 536 elements, which leads to 4.19M DOFs. All simulations are able to capture the main features of the flow. In particular, a good agreement with the reference data is found for Re b = 10 595 and Re b = 19 000. The discrepancies found for Re b = 37 000 between the DG simulation and the experiment are in line with those reported in [START_REF] Manhart | Assessment of eddy resolving techniques for the flow over periodically arranged hills up to Re[END_REF] using the MGLET FV code on a mesh composed of 4.1M points and in [START_REF] Chaouat | Hybrid RANS-LES Simulations of the Turbulent Flow over Periodic Hills at High Reynolds Number[END_REF] using a hybrid RANS/LES approach on a 1 M point grid. An additional no-model LES simulation has been performed at Re b = 10 595. The solutions from the no-model and the WALE approaches do not appear to be fundamentally different. The WALE approach, however, seems to provide slightly more accurate results, specially in the vicinity of the separation point and in the recirculation area. The smoothing effect that the addition of the subgrid model has on the u u profiles has also been highlighted. Indeed, at all stations, we observe that the unphysical peaks generated by the interface jumps of the DG solution are considerably damped by the use of a subgrid model. Note that for Re b = 19 000 and 37 000 the no-model approach did not yield stable simulations. This is a manifestation of the very low dissipation level provided by the DG scheme, and the need to use subgrid-scale modelling when we move towards higher Reynolds numbers. We see from these results that the combined used of local hp-adaptation with a subgrid-scale modelling approach could fully exploit the benefits of DG methods for industrial LES. Furthermore, in order to reduce the computational cost associated with the use of high polynomial degrees it would be advantageous to use a hybrid MPI/openMP approach, instead of the pure MPI approach used in this work. This is the subject of current research in the framework of the Aghora project.

x/h = 0.05

x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (a) Mean streamwise velocity U/u b profiles

x/h = 0.05

x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (b) Mean vertical velocity V /u b profiles x/h = 0.05

x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (a) Reynolds stress u u /u 2 b profiles

x/h = 0.05 x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (b) Shear stress u v /u 2 b profiles x/h = 0.05

x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (a) Mean streamwise velocity U/u b profiles

x/h = 0.05

x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (b) Mean vertical velocity V /u b profiles x/h = 0.05

x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (a) Reynolds stress u u /u 2 b profiles

x/h = 0.05 x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (b) Shear stress u v /u 2 b profiles x/h = 0.05

x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (a) Mean streamwise velocity U/u b profiles

x/h = 0.05

x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (b) Mean vertical velocity V /u b profiles

29

x/h = 0.05

x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (a) Reynolds stress u u /u 2 b profiles

x/h = 0.05 x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (b) Shear stress u v /u 2 b profiles x/h = 0.05

x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (a) Mean streamwise velocity U/u b profiles

x/h = 0.05

x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (b) Mean vertical velocity V /u b profiles x/h = 0.05

x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (a) Reynolds stress u u /u 2 b profiles

x/h = 0.05 x/h = 0.5

x/h = 2

x/h = 4

x/h = 6

x/h = 8 (b) Shear stress u v /u 2 b profiles

Figure 1 :

 1 Figure 1: Trace of a variable on a face.

Figure 2 :

 2 Figure 2: Profiles of mean velocity U + (left) and rms of the velocity fluctuations u + rms , v + rms , w + rms (right) for DG computations of the channel flow. Comparison against reference DNS data.

Figure 3 :

 3 Figure 3: Profiles of cross correlations -uv + for DG computations of the channel flow. Comparison against reference DNS data.

Figure 4 :

 4 Figure 4: Profiles of mean scalar θ + (left) and rms θ + rms and correlations of the scalar and velocity fluctuations uθ + and vθ + (right) for DG -p4 computation of the channel flow. Comparison against reference DNS data.

(a) Coarse mesh 32 × 32 Figure 5 :

 32325 Figure 5: Fourth-order meshes used in the DG simulations.

Figure 6 :

 6 Figure 6: Streamwise distribution of y + (right) and ∆x + , ∆z + (left) in the DG simulations. Re b LES approach Simulation # DOFs y + max 2800 -DG-p4 fine 1.02 M 0.9

( a )

 a Re b = 2800 (b) Re b = 10 595 (b) Re b = 19 000 (c) Re b = 37 000

Figure 7 :

 7 Figure 7: Streamlines of the time-averaged streamwise velocity field U/u b from DG simulations (see Table4).

( a )

 a Re b = 2800 (b) Re b = 2800 up to Re b = 37 000

Figure 8 :

 8 Figure 8: Profiles of wall-shear stress at bottom wall for the DG simulations.

6. 3 .

 3 DNS results at Re b = 2 800 Table 5 summarises the five DNS simulations performed at Re b = 2 800. The CPU time required to advance one convective time is also shown in this table.

( a )Figure 9 :

 a9 Figure 9: Instantaneous field of vorticity magnitude at Re b = 2 800 and 10 595.

  ) and (b) show the field of instantaneous vorticity magnitude displayed by the DG solutions at Re b = 2 800 and 10 595, respectively. We can see on these plots the finer structures present in the flow at this higher Re b . 6.5. LES results at Re b = 19 000

Figure 10 :

 10 Figure 10: DNS at Re b = 2 800. Comparison of DG simulations with reference data. • : LESOCC O(2); --: DG -p3 coarse; --: DG -p4 coarse; --: DG -p5 coarse; -: DG -p3 fine; -: DG -p4 fine.

Figure 11 :

 11 Figure 11: DNS at Re b = 2 800. Comparison of DG simulations with reference data. • : LESOCC O(2); --: DG -p3 coarse; --: DG -p4 coarse; --: DG -p5 coarse; -: DG -p3 fine; -: DG -p4 fine.

Figure 12 :

 12 Figure 12: LES at Re b = 10 595. Comparison of DG simulations with reference data. • : Experiment; • : LESOCC O(2) (13.1 M dofs); -: No-model DG -p3 fine (4.19 M dofs); -: WALE DG -p3 fine (4.19 M dofs).

Figure 13 :

 13 Figure 13: LES at Re b = 10 595. Comparison of DG with reference data. • : Experiment; • : LESOCC O(2) (13.1 M dofs); -: No-model DG -p3 fine (4.19 M dofs); -: WALE DG -p3 fine (4.19 M dofs).

Figure 14 :

 14 Figure 14: LES at Re = 19 000. Comparison of DG simulations with experimental data. • : Experiment; -: WALE DG -p3 fine (4.19 Mdofs).

Figure 15 :

 15 Figure 15: LES at Re = 19 000. Comparison of DG simulations with experimental data. • : Experiment; -: WALE DG -p3 simulation on fine mesh (4.19 Mdofs).

Figure 16 :

 16 Figure 16: LES at Re = 37 000. Comparison of DG simulations with reference MGLET simulation and experimental data. • : Experiment; • : MGLET O(2) (4.1 Mdofs); -: WALE DG -p3 fine (4.19 Mdofs).

Figure 17 :

 17 Figure 17: LES at Re = 37 000. Comparison of DG simulations with reference MGLET simulation and experimental data. • : Experiment; • : MGLET O(2) (4.1 Mdofs); -: WALE DG -p3 fine (4.19 Mdofs).

Table 1 :

 1 Computational grids.

	Simulation	Order (p + 1)	# Elements	# DOFs
	Abe et al. [49]	4	1024 × 256 × 1024 268.4 M
	DG-p2	3	84 × 22 × 84	4.19 M
	DG-p3	4	64 × 16 × 64	4.19 M
	DG-p4	5	52 × 12 × 52	4.06 M
	Simulation ∆t (u b /h)/10 -4 CP U time/tc Re τ
	DG-p2	5.0	1.00	619
	DG-p3	8.0	1.16	676
	DG-p4	6.7	1.36	647

Table 2 :

 2 Computational details for the plane channel flow DG simulations.

Table 3 :

 3 Grid resolutions used in the reference and DG computations.

		Grid resolution	# Elements Grid order
	LESOCC	-	12.4 M	1
	MGLET	216×168×104	3.77 M	1
	Coarse mesh	32×16×16	8 192	4
	Fine mesh	64×32×32	65 536	4

Table 4 :

 4 Wall resolution y + max at point of maximum wall-shear stress for the DG simulations.

  -4 #Cores CPU time/tc

	LESOCC FV-p1	13.1 M	20	-	-
	DG-p3 Coarse	0.52 M	2.5	264	1.0
	DG-p4 Coarse	1.02 M	2.5	264	2.8
	DG-p5 Coarse	1.77 M	1.7	264	13.1
	DG-p3 Fine	4.19 M	1.8	2208	10.8
	DG-p4 Fine	8.19 M	1.4	2208	39.9

Table 5 :

 5 Computational details of DG-DNS computations at Re b = 2 800.

	Rapp [40].

6.2.1. Qualitative comparison of DG solutions at Re b =2 800 up to 37 000

The heat-flux associated with the passive scalar Qw should not be mistaken for the heatflux q that appears in the compressible energy equation.

Acknowledgments

This research was partly funded by the European H2020 TILDA project, with additional support from Safran Group. The authors would like to thank Dr. E. Martin for his help concerning parallel computing issues. We also thank Dr. J.-B. Chapelier for his help with the plane channel configuration during his postdoctoral work at Onera. This work was performed using HPC resources from GENCI (Grant 2015 -c20152a7441).