{'fraction_non_alphanumeric': 0.03947562958947316, 'fraction_numerical': 0.012838204129909813, 'mean_word_length': 4.590714974514396, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the ability of discontinuous Galerkin methods to simulate under-resolved turbulent flows

We investigate the ability of discontinuous Galerkin (DG) methods to simulate under-resolved turbulent flows in large-eddy simulation. The role of the Riemann solver and the subgrid-scale model in the prediction of a variety of flow regimes, including transition to turbulence, wall-free turbulence and wall-bounded turbulence, are examined. Numerical and theoretical results show the Riemann solver in the DG scheme plays the role of an implicit subgrid-scale model and introduces numerical dissipation in under-resolved turbulent regions of the flow. This implicit model behaves like a dynamic model and vanishes for flows that do not contain subgrid scales, such as laminar flows, which is a critical feature to accurately predict transition to turbulence. In addition, for the moderate-Reynolds-number turbulence problems considered, the implicit model provides a more accurate representation of the actual subgrid scales in the flow than state-of-the-art explicit eddy viscosity models, including dynamic Smagorinsky, WALE and Vreman. The results in this paper indicate new best practices for subgrid-scale modeling are needed with high-order DG methods.


Introduction
Over the past few years, discontinuous Galerkin (DG) methods for large-eddy simulation (LES) have emerged as a promising approach to solve complex turbulent flows. DG methods allow for high-order discretizations on complex geometries and unstructured meshes. This is critical for the simulation of transitional and turbulent flows in industrial applications, which require that small-scale small-magnitude features are accurately propagated over complex three-dimensional geometries. In addition, DG methods are well-suited to emerging computing architectures, including graphics processing units (GPUs) and other many-core architectures, due to their high flop-to-communication ratio [1,2]. The use of DG methods for LES is being further encouraged by successful numerical predictions [8,20,21,22,23,26,33,59,69,76,81].
Despite the significant research investment on discontinuous Galerkin methods for large-eddy simulation [9,22,25,51,58,82], the precise roles of the subgrid-scale (SGS) model (if any) and of the Riemann solver used by the DG scheme in the ability to predict transitional and turbulent flows remains unclear. In this paper, we investigate the ability of discontinuous Galerkin methods to simulate a variety of flow regimes, including transition to turbulence, under-resolved 1 wall-free turbulence, and under-resolved wall-bounded turbulence. The Taylor-Green vortex problem [75] and the turbulent channel flow [46] at various Reynolds numbers are considered to that end.
The remainder of the paper is structured as follows. In Section 2, we summarize the test problems and studies performed. Numerical results for the Taylor-Green vortex and the turbulent channel flow are presented in Sections 3 and 4, respectively. We conclude with some remarks in Section 5. The details of the DG discretization are presented in Appendix A.

Summary of test cases and studies performed
We consider the nearly incompressible Taylor-Green vortex [75] and turbulent channel flow [46] problems at various Reynolds numbers. The focus in the Taylor-Green vortex is to investigate the effect of the Riemann solver and the SGS model on the dissipation of kinetic energy. The focus in the turbulent channel flow is to investigate the effect on the turbulent transport. In both test problems, the fluid is assumed to be Newtonian, calorically perfect, in thermodynamic equilibrium, with Fourier's law of heat conduction, and with the Stokes' hypothesis, as discussed in Appendix A. The dynamic viscosity is constant, the Prandtl number = 0.71, and the ratio of specific heats = / = 1.4. The complete description the problems is presented in Sections 3.1 and 4.1.
We focus on third-order DG methods. For this accuracy order, the large scales in the flow are affected less by the numerical dissipation of the DG scheme than by the viscous dissipation due to the eddy viscosity in the explicit SGS models considered 2 . A discussion on how our results are expected to extend to higher accuracy orders is presented in Section 5. Due to their lower computational cost for moderately high accuracy orders, we use hybridized DG methods [21], a class of discontinuous Galerkin methods that includes the HDG [60,67], the EDG [68] and the IEDG [19] methods. The details of the hybridized DG methods for the compressible Euler and Navier-Stokes equations are presented in Appendix A. The third-order, three-stage -stable diagonally implicit Runge-Kutta DIRK (3,3) scheme [4] is used for the temporal discretization, and the time-step size is chosen sufficiently small so that the spatial discretization error dominates the time discretization error [22].
For each test problem, we perform two studies: One for the Riemann solver and another for the SGS model. For the Riemann solver studies, we consider the following stabilization matrices in the hybridized DG scheme where = ( · )/ is the Jacobian matrix of the inviscid flux normal to the element face, denotes the maximum-magnitude eigenvalue of , | · | is the generalized absolute value operator, is the identity matrix, and̂︀ ℎ is the approximation to the trace of the solution on the element faces, as described in Appendix A. We note that the stabilization matrix implicitly defines the Riemann solver in 1 As is customary, we use the term under-resolved to refer to simulations in which the exact solution contains scales that are smaller than the grid resolution (the so-called subgrid scales). 2 The numerical dissipation in an ℓ-th order DG method vanishes at the rate ( 2ℓ ) in the small wavenumber limit → 0 [3]. The viscous dissipation of a second-order operator, such as the viscous operator of the Navier-Stokes equations, vanishes at the rate ( 2 ). This is, the decay rate of a signal is proportional to the 2ℓ-th power of its wavenumber and to the square of its wavenumber, respectively. Regarding dissipation at high wavenumbers, the numerical dissipation in high-order DG methods is more localized near the grid Nyquist wavenumber than the viscous dissipation [25]. hybridized DG methods. Additional details on these stabilization matrices and on the relationship between the stabilization matrix and the resulting Riemann solver are presented in [22]. For the purpose of this paper, we note that the two first stabilization matrices lead to Roe-type Riemann solvers, and the third one yields a Lax-Friedrichs-type solver. No explicit SGS models are used for the Riemann solver studies.
For the SGS model studies, we focus on state-of-the-art eddy viscosity models. These models are based on the Boussinesq eddy viscosity assumption and enter the governing equations through an augmented viscous operator. In particular, the static Smagorinsky [72], dynamic Smagorinsky [48], WALE [62] and Vreman [79] models are considered, in addition to implicit LES without an explicit model. The explicit models are further equipped with the Yoshizawa model [80] for the isotropic part of the SGS stress tensor, the Knight model [41] for the turbulent diffusion, and an SGS eddy Prandtl number approach with = 0.7 for the subgrid-scale heat transfer. All other extra terms arising in the Favre-filtered Navier-Stokes equations [32] are not modeled due to their negligible magnitude compared to the previous terms [52]. The length scale involved in the SGS models is set to ∆ = 1/3 / , where is the volume of the element and the polynomial order of the DG approximation. For the dynamic Smagorinsky model, projection onto the space of polynomials of degree ′ = ⌊ /2⌋ is used for the coarse-graining, where ⌊ /2⌋ is the greatest integer less than or equal to /2, i.e. ′ = 1 in our case. The stabilization matrix (1a) is used for the SGS studies.
Compressibility effects will be neglected for the two following purposes: (i) Incompressible DNS results will be used as reference data, and (ii) spatial filtering and Favre filtering will be assumed to be equivalent in the presentation of the numerical results. We finally note that, for some of the Reynolds numbers considered, the mesh resolution is intentionally insufficient to match the DNS data, thus allowing for a more meaningful analysis of the role of the implicit and explicit models in under-resolved simulations.

Case description
The Taylor-Green vortex (TGV) [75] is a canonical problem in fluid mechanics commonly used to study vortex dynamics, turbulence transition, and turbulence decay. It contains several flow regimes in a single construct and is therefore an excellent test case for our purpose. In particular, the TGV problem describes the evolution of the flow in a cubic domain Ω = [− , ) 3 with triple periodic boundaries, starting from the smooth initial condition where , and = ( , , ) denote density, pressure and the velocity vector, respectively, and 0 , 0 , 0 > 0 are some reference density, pressure and velocity magnitude. The large-scale eddy in the initial condition leads to smaller and smaller structures through vortex stretching. Near ≈ 7 / 0 , the vortical structures undergo structural changes, and around ≈ 8 − 9 / 0 they break down and the flow transitions to turbulence 3 , where the exact times depend on the Reynolds number = 0 0 / . The Reynolds numbers = 1600 case, the smallest turbulent structures and the maximum dissipation rate of kinetic energy occur at ≈ 9 / 0 . After this time, the turbulent motion dissipates all the kinetic energy and the flow eventually comes to rest through a decay phase similar to that in decaying homogeneous isotropic turbulence, yet not isotropic here. At this Reynolds number, the subgrid scales are moderate compared to the resolved scales. In the inviscid TGV, however, the smallest turbulent scales become arbitrarily small, and the range of subgrid scales arbitrarily large, as time evolves.
The reference Mach number is set to 0 = 0 / 0 = 0.1, where 0 denotes the speed of sound at temperature 0 . This completes the non-dimensional description of the problem.

Details of the numerical discretization
The computational domain is partitioned into a uniform 64 × 64 × 64 Cartesian grid and the third-order Embedded DG (EDG) scheme [68] is used for the spatial discretization. The solution is computed from 0 = 0 to = 15 / 0 . Three different phases exist in the simulation. Before ≈ 4 / 0 , the flow is laminar and with no subgrid scales. This is followed by an under-resolved laminar phase (i.e. with subgrid scales) that lasts until ≈ 7 − 9 / 0 . From then on, the flow is turbulent and under-resolved.

Numerical results
3.3.1. Riemann solver study Figure 1 shows the time evolution of kinetic energy dissipation rate, in the viscous (left) and inviscid (right) Taylor-Green vortex for the three Riemann solvers considered. Figure 2 shows the time evolution of the quantity where = ∇ × denotes vorticity. The second term in the right-hand side of (4) corresponds to the dissipation of kinetic energy due to physical mechanisms. Π is therefore the contribution of the numerical scheme to the dissipation of kinetic energy, and is referred to as the numerical dissipation of kinetic energy. We note that Π should approximately account for the transfer of kinetic energy from resolved scales to subgrid scales 4 . Also note that Eq. (4) is derived from the incompressible kinetic energy equation since numerical results suggest the incompressible kinetic energy equation is more appropriate than the compressible kinetic energy equation to assess numerical dissipation in nearly incompressible flows [22]. Tables 1, 2, 3 and 4 collect the average absolute-value jump across elements on the periodic plane = − for each conservation variable = 1, ..., 5, defined as at = 1600 and = 3 / 0 , = 1600 and = 8 / 0 , = ∞ and = 3 / 0 , and = ∞ and = 8 / 0 , respectively. Here, ℎ, = + ℎ, | − − ℎ, | and ⟨ ℎ, ⟩ = ( + ℎ, | + − ℎ, | )/2 denote the face jump and face average operators. The term | is used to emphasize that the DG solution varies inside each element and is to be evaluated on the element face. From these figures and tables, several remarks follow.  (i) The DG scheme adds numerical dissipation when the exact solution contains subgrid scales, that is, after ≈ 4 / 0 . This suggests discontinuous Galerkin methods have a built-in (implicit) subgrid-scale model and introduce additional dissipation in the presence of under-resolved turbulence.
(ii) The numerical dissipation is negligible before subgrid scales appear in the flow. Unlike the explicit SGS models in the next section, the implicit model succeeds to detect there are no subgrid scales and does not add numerical dissipation under those conditions. Observations (i) and (ii) imply the numerical dissipation of the DG scheme is positive for under-resolved turbulent flows and vanishes for laminar flows, that is, the implicit model behaves like a dynamical model. This dynamic behavior is justified as follows: On the one hand, if the exact solution does not contain subgrid scales, it is well represented in the DG approximation space and the scheme is in the asymptotic convergence regime. This implies the inter-element jumps are small and in particular || ℎ || = (ℎ +1 ) [60,67,68], where ℎ denotes the element size and is the polynomial order of the DG approximation. Indeed, we note the very small magnitude of the inter-element jumps in Tables 1 and 3. Since the amount of numerical dissipation per unit element face area in DG methods is of order (|| ℎ || 2 ) [7, 23,24], it is therefore  negligible when there are no subgrid scales. On the other hand, when the exact solution contains subgrid scales (i.e. when the simulation becomes under-resolved), the inter-element jumps grow and stabilize the scheme by adding numerical dissipation. Hence, the Riemann solver plays the role of a dynamic SGS model and accounts for the effect of the subgrid scales in a similar way as explicit models do. We note that, while subgrid scales are most commonly encountered in the simulation of turbulent flows, they may also exist in laminar flows, such as in the Taylor-Green vortex between ≈ 4 / 0 and ≈ 7 − 9 / 0 .
(iii) Despite under-resolution, no significant differences between Riemann solvers are observed in the viscous case. Even in the inviscid limit the role of the Riemann solver is still moderate. We note that the stabilization matrix = yields much smaller jumps in the momentum fields than the other two stabilization matrices. This is due to the over-upwinding for the momentum equations provided by Lax-Friedrichstype Riemann solvers at low Mach numbers [58], and indicates that DG methods have a (nonlinear) autocorrection mechanism that adapts the magnitude of the inter-element jumps to partially compensate for overshoots in the Riemann solver. This auto-correction mechanism in turn justifies the minor role of the Riemann solver.
(iv) All the Riemann solvers slightly underestimate the peak dissipation with respect to the direct numerical simulation (DNS) data [77]. That is, the numerical dissipation is smaller than the true SGS dissipation (only) when the smallest turbulent structures appear. Under those conditions, the DG scheme could benefit from the addition of an explicit SGS model. However, as will be discussed in Section 3.3.2, this is not the case in practice since the implicit model is partially inhibited by the use of an explicit model. Figure 3 shows the one-dimensional kinetic energy spectra at = 8 / 0 and = 9 / 0 for the viscous and inviscid Taylor-Green vortex, respectively. In all cases, we observe an inertial range in which the spectrum follows a power law with exponent close to the theoretical value of −5/3 [42]. The inertial range is followed  by a dissipative range, until the grid Nyquist wavenumber = 96.0/ is achieved and no smaller scales exist in the discretization. In the viscous case, the inertial range extends up to ≈ 20/ , whereas DNS results [77] indicate it extends until ≈ 40/ . This numerically-induced premature end of the inertial range is predicted by eigenanalysis [55], Fourier analysis [5] and non-modal analysis [25] theory. Also, note that no significant differences in the spectra are observed between Riemann solvers in the viscous case (neither at the times shown nor at other times); which shows that at this Reynolds number the viscous dissipation has a much larger impact on the dynamics of all the scales than the Riemann solver. We also note the presence of an energy pileup at large wavenumbers in the inviscid case for the Riemann solver = . These pre-dissipative bumps are predicted by eigenanalysis [55] and non-modal analysis [25] theory, and are consistent with results in the literature for inviscid low Mach numbers flows when using Riemann solvers that are based on the maximum-magnitude eigenvalue of [56,58]. From these results and the insights from eigenanalysis and non-modal analysis, the Riemann solver has some impact on the dynamics of the smallest resolved scales, particularly at low Mach numbers and high Reynolds numbers (more precisely, at low Mach numbers and high cell Péclet numbers); whereas the large-scale dynamics are affected to a much lesser extent by the Riemann solver. As the Mach number increases and the cell Péclet number decreases, the choice of Riemann solver is expected to have a smaller impact on the behavior of the small scales.  Tables 1, 2, 3 and 4. Also, Figure 6 shows snapshots of the vorticity norm on = − for the viscous case at = 8 / 0 . We recall that the stabilization matrix (1a) is used for the SGS study. These results can be summarized as follows:

Subgrid-scale model study
(i) The eddy viscosity in all the explicit SGS models fails to vanish when there are no subgrid scales in the flow, and produces unphysical dissipation during this phase and also during the under-resolved laminar phase. Among the explicit models, dynamic Smagorinsky introduces the least amount of dissipation in these two phases. We note that a larger dissipation of kinetic energy corresponds to effectively solving a lower Reynolds number flow; which is consistent with the maximum dissipation rate occuring at an earlier time that is characteristic of the Taylor-Green vortex at lower Reynolds numbers [13].
We emphasize that the built-in stabilization due to inter-element jumps (i.e. the implicit subgrid-scale model) in the DG scheme provides a more accurate mechanism to detect the absence of subgrid scales, in which case very little dissipation is added, than the explicit models. This is a critical advantage of the implicit model to simulate transitional flows.
(ii) When subgrid scales appear in the flow, the amount of dissipation introduced by the implicit model is closer to the true SGS value than that introduced by the explicit models. Only in the inviscid case in fully turbulent regime, the dynamic Smagorinsky model performs similarly to the implicit model. We note that, whenever the simulation is under-resolved, the inter-element jumps with an explicit SGS model are much smaller than with no model (see Tables 1−4); which is due to the additional stabilization provided by the eddy viscosity. The use of an explicit model therefore partially inhibits the implicit model, and in fact the dissipation of kinetic energy in the turbulent regime is smaller with the explicit models.
(iii) While the inter-element jumps in the numerical solution are smaller and the vorticity norm field is smoother with an explicit SGS model (due to the eddy viscosity dissipation), ILES provides more accurate results. Therefore, lack of smoothness in the DG solution is not an indicator for low solution quality -and, as discussed previously, the inter-element jumps are actually responsible for the built-in model in the scheme.
(iv) Regarding the relative performance of the explicit models, no major differences are observed between static Smagorinsky and Vreman. The WALE model led to nonlinear instability and the simulation breakdown at ≈ 4.59 / 0 and 2.75 / 0 in the viscous and inviscid cases, respectively. This lack of robustness is due to the high nonlinearity of the WALE model and may limit its applicability with high-order DG methods. The dynamic Smagorinsky model provides the most accurate representation of the subgrid scales among the explicit models.    Figure 7 shows the one-dimensional kinetic energy spectra for ILES, the static Smagorinsky, dynamic Smagorinsky and Vreman models at = 8 / 0 and = 9 / 0 for the viscous and inviscid cases, respectively. The WALE model crashed before these times. We recall that the grid Nyquist wavenumber is = 96.0/ . The kinetic energy spectrum in ILES agrees with the theoretical −5/3 slope of decay of the inertial range for a larger range of wavenumbers than the explicit models. Small differences are observed between the static Smagorinsky, WALE and Vreman models, both at the times shown as well as at all other times. The dynamic Smagorinsky spectrum is in between that of the implicit LES and those of the static models. All the explicit models, especially the static ones, dissipate kinetic energy at larger scales than the implicit model, and in particular at scales that are much larger than the grid Nyquist wavenumber. This is consistent with eigenanalysis [3,55], Fourier analysis [5] and non-modal analysis [25]. The fact that explicit models dissipate energy at larger scales than the implicit model may have important consequences in practice, as discussed in Sections 4 and 5.

Summary
The Taylor-Green vortex results indicate that: • Discontinuous Galerkin methods have a built-in (implicit) subgrid-scale model and numerically dissipate kinetic energy in under-resolved turbulence simulations. The inter-element jumps and the Riemann solver are responsible for the implicit model.
• The amount of dissipation implicitly introduced by the DG scheme is closer to the actual dissipation in the subgrid scales than that with the explicit models considered.
• The implicit dissipation is more localized near the grid Nyquist wavenumber (i.e. more localized in the smallest resolved scales) than that introduced by the explicit models.
• The implicit model does not add dissipation when there are no subgrid scales, whereas the explicit models do. That is, the implicit model behaves like a dynamic model.
Except for the second remark above, these numerical observations have been justified by theoretical studies of DG methods. The second observation could be explained by relating the DG stabilization due to the inter-element jumps with the subgrid-scale closure terms arising from variational multiscale (VMS) [37] or Mori-Zwanzig (MZ) [63,64] approaches. For one-dimensional linear convection, the MZ-VMS procedure with the assumptions of finite memory and linear quadrature, referred to as -MZ-VMS, actually leads to a subgrid-scale closure term that is equivalent to the standard upwind flux [65], i.e. the implicit model with standard upwinding is the same as that given by -MZ-VMS. This analogy is less straightforward in the case of the Navier-Stokes equations and has not been investigated in this paper.

Case description
We consider the turbulent channel flow [46]  As is customary in channel flows, the three velocity components are denoted by ( , , ), the time-averaged velocity by a capital letter, the fluctuations by a prime, and the ensemble average by ⟨ · ⟩. Thus, = ⟨ ⟩ and = + ′ , and similarly for the wall-normal and spanwise velocity components.

Details of the numerical discretization
The channel flow is simulated in a doubly-periodic domain Ω = [0, 4 ) × [0, 2 ) × [0, 2 ). Periodicity is imposed along the (streamwise) and (spanwise) directions. The computational domain is partitioned into a 48 × 32 × 40 Cartesian grid and the third-order Interior Embedded DG (IEDG) scheme [19] is used for the spatial discretization. The element size is constant along the streamwise and spanwise directions. The high-order nodes in the wall-normal direction are uniformly distributed in a mapped coordinate that is related to through = sin( /2) sin( /2) The distance between high-order nodes (in wall units) is summarized in Table 5. As is customary, the nearwall velocity, time and length scales for non-dimensionalization are , / 0 2 and / 0 , respectively, and the superscript + is used to indicate that a quantity is expressed in wall units. A run-up time + 1 = 2000 is used for the flow to achieve its stationary distribution on the chaotic attractor. The flow statistics are then collected over a time window + 2 = 1000 to ensure statistical convergence of the mean velocity and the Reynolds stresses.
We note that the pressure gradient is such that the mass flow is the one that led to the desired in DNS [46]. As a consequence, the computed from the wall stress in LES (denoted by hereinafter) may y + not exactly agree with the target (denoted by ) if the resolution is not sufficiently fine to match the wall stress in DNS. This will be the case at = 544 and, for some of the explicit models, also at 182.

Numerical results
Prior to presenting the numerical results, we note that homogeneity in the and directions holds pointwise in the exact solution, but only elementwise in the DG solution. We shall omit this nuance and compute ensemble averages (i.e. time averages) as time-, streamwise-and spanwise-averages. This accelerates convergence of the statistics of the flow, and allows comparison with DNS results for which the ensembleaveraged quantities are only a function of the wall-normal coordinate. Figure 8 shows the mean velocity profile at = 182 and 544 for the Riemann solvers considered. The Reynolds stresses at = 182 are shown in Figure 9, and the values of at both Reynolds numbers are collected in Table 6. There are no significant differences between Riemann solvers in terms of wall friction, mean velocity and Reynolds stresses. Further analysis of these results is presented in the comparison between implicit and explicit models below.  Figures 10 and 11 show the mean velocity profile at = 182 and 544, respectively, with ILES, the static Smagorinsky, dynamic Smagorinsky, WALE and Vreman models. Two different types of nondimensionalization are used on the left and right images. On the left images, the inner boundary layer scales, namely the near-wall velocity and length scales, are used for non-dimensionalization. These scales are computed from the wall friction in the simulations, as opposed to the target wall friction. This is the proper choice to investigate if the viscous, buffer and log sublayers of the boundary layer are accurately resolved, as we will discuss below. On the right images, the outer boundary layer scales, namely and , are used for non-dimensionalization so that all simulations are non-dimensionalized with respect to the same reference values. This way, the scaling factor between dimensional and non-dimensional data is the same in all cases, and the dimensional velocity profiles can be directly compared. Both choices of nondimensionalization complement each other and help understand the performance of the models and the reasons for the mismatch (if any) between LES and DNS. Figure 12 shows the Reynolds stresses at = 182, Figure 13 the dynamic Smagorinsky constant at both Reynolds numbers, and Table 6   these results, we introduce the two following quantities

Subgrid-scale model study
From these definitions, it follows that / = − 0 ⟨ ′ ′ ⟩ and ( + * ) / = (1− / ). Note that the mean shear stress varies linearly across the channel due to the uniform pressure gradient and -momentum conservation, and therefore (1 − / ) = ⟨ ⟩. Note also that − 0 ⟨ ′ ′ ⟩ ≈ −⟨ ′′ ′′ ⟩, where ′′ and ′′ are the Favre fluctuating velocities, and thus this term corresponds to the net resolved -momentum turbulent transport, per unit area, across -planes. From these considerations, can be interpreted as an eddy viscosity due to the resolved turbulent motion, and will be referred to as the resolved eddy viscosity; whereas * can be regarded as a total eddy viscosity or an effective eddy viscosity in the simulation, and accounts for the turbulent transport due to the explicit model, the implicit model and the resolved turbulence. Note that, for the exact solution, it holds that * = . Also, a modeled eddy viscosity, accounting for the implicit and explicit models, is not straightforward to define in the context of DG methods and in fact it cannot be computed as * − .
The resolved eddy viscosity in the viscous and buffer layers for the SGS models considered is shown in the top of Figure 14. The total eddy viscosity is shown in the bottom of the figure. Except for static Smagorinsky, * vanishes as expected in the viscous layer, and grows in the buffer layer as the turbulent transport increases and dominates the molecular transport. We emphasize that these eddy viscosities inform of the turbulent transport, and not of the dissipation of kinetic energy which was the focus in the Taylor-Green vortex. While eddy viscosity SGS models account for the subgrid scale contribution to these two terms through an (the same) explicit eddy viscosity, these terms are of a fundamentally different nature, as apparent from the filtered momentum and kinetic energy equations [32].
All the ingredients are now in place to discuss the numerical results. We focus first on the viscous sublayer, then on the buffer layer, and finally on the log layer.
Viscous sublayer ( + 10): This layer is characterized by ⟨ ⟩ ≈ and the fact that molecular transport dominates turbulent transport; which yields the well-known result + ≈ + . Since this velocity law is independent of the Reynolds number, all simulations should ideally match the DNS data on the left images of Figures 10−11, despite the slightly different Reynolds number in LES. From these two figures, the viscous sublayer is accurately resolved at both Reynolds numbers by implicit LES, the dynamic Smagorinsky, the WALE and the Vreman models. Provided that there is enough resolution to capture the mean flow and that no eddy viscosity is added by the SGS model, this is expected due to the lack of subgrid scales in this laminar sublayer. The eddy viscosity, however, does not vanish near the wall with static Smagorinsky (see * in Figure 14) and this leads to inaccurate wall friction and mean velocity in this layer. While no significant differences were observed between static Smagorinsky and the other static models, namely WALE and Vreman, for the wall-free flow conditions in the Taylor-Green vortex, this is not the case for this wall-bounded flow and, as expected, points to the inability of the static Smagorinsky model to vanish near walls as one the main limitations of the model.
Buffer layer (10 + 40): The mean velocity + = + ( + ) in the buffer layer is also approximately independent of the Reynolds number [10,18,36,46,50,54], and therefore LES should ideally match DNS on the left images of Figures 10−11. At = 182, the buffer layer is only accurately resolved with ILES. At = 544, the buffer layer is not accurately resolved in any of the LES simulations. From the top images in Figure 14, the use of an explicit model reduces the resolved turbulent transport. As discussed before and justified by linear analysis techniques [5,55,25], the explicit models damp the large scales more than the implicit model does. Since the large turbulent scales are responsible for most of the turbulent transport, this justifies the smaller . From the bottom images in Figure 14, the total eddy viscosity * is y + also smaller with an explicit model, and so is the total turbulent transport.
Log layer ( + 40): This layer is characterized by ⟨ ⟩ ≈ and the fact that turbulent transport dominates molecular transport. This, combined with the mixing-length hypothesis ℓ = , where ℓ is the mixing length, is the von Kármán constant and is the distance to the wall, leads to + ≈ −1 log + + , where , > 0 are positive constants that, for channel flows, are approximately independent of the Reynolds number. DNS [10,18,36,46,50,54] and experimental [53,71] data have confirmed the existence of this log layer, and in particular ≈ 0.384 and ≈ 4.27 [46]. The left images in Figures 10 and 11 show that all simulations succeed to predict this logarithmic dependence with the correct mixing length ℓ = , ≈ 0.38. The matching velocity between the buffer and log layers, however, does not agree with DNS, and thus ̸ = 4.27, if the buffer layer is not accurately resolved. This leads to an incorrect outer velocity seen by the viscous sublayer; which is responsible for the misprediction of the wall friction and thus in Table 6.
To conclude, we note that the Reynolds stresses have a stronger dependence on the Reynolds number [46]. Since the simulations with different SGS models predict slightly different wall frictions, and thus solve slightly different Reynolds numbers, it is challenging to infer additional conclusions, beyond those already y +

Conclusions
We investigated the ability of discontinuous Galerkin methods to predict under-resolved turbulent flows. The Taylor-Green vortex and the turbulent channel flow at various Reynolds numbers were considered to this end. Numerical results showed that DG methods without an explicit subgrid-scale model (implicitly) introduce numerical dissipation in under-resolved turbulence simulations. This implicit subgrid-scale model is due to the inter-element jumps and the Riemann solver, and behaves like a dynamic model in the sense that it vanishes for laminar flows that do not contain subgrid scales; which is a critical feature to accurately simulate transitional flows. In addition, for the moderate-Reynolds-number turbulence problems considered, the implicit model provided a more accurate representation of the actual subgrid scales than state-of-the-art explicit eddy viscosity models. Theoretical results for DG methods were used to justify these numerical observations. Some premises that are widely accepted in the LES community do no longer hold in the context of high-order DG methods. First, the built-in subgrid-scale model in the DG scheme is partially inhibited when using an explicit model, and the total amount of dissipation does not necessarily increase with an explicit model. Second, since eddy viscosity models dissipate kinetic energy at larger scales than the implicit model, they reduce more significantly the energy content of scales that are larger than the grid Nyquist wavenumber; which may have important consequences in practice. In particular, explicit eddy viscosity models may not allow taking advantage of the low numerical dissipation at large scales of high-order DG methods [25,55]; which is critical for transition prediction and moderate-Reynolds-number turbulence. From these two observations, a change in the current best practices for subgrid-scale modeling may be required with DG methods. This is not completely surprising considering that the state-of-the-art subgrid-scale models, based on the Boussinesq eddy viscosity assumption and an augmented viscous operator, have been developed and successfully applied with discretization schemes whose built-in stabilization (if any) and dissipation characteristics are different from those in discontinuous Galerkin methods.
To conclude, we briefly discuss how our results are expected to extend to higher accuracy orders. For accuracy orders up to about seven, similar results to those in this paper are expected. For accuracy orders beyond about eight, the numerical dissipation in the DG scheme becomes very small also at high wavenumbers [25], and some form of regularization/model may be required to enhance stability and accuracy. It will be critical that the explicit regularization/model used with high-order DG schemes localizes dissipation at the desired wavenumbers and is consistent with the expected SGS dissipation spectrum [14,43,45]; which could be achieved with Spectral Vanishing Vicosity [74,39,40,57], Variational Multiscale [16,37,38] and Mori-Zwanzig [63,64,65]  where is the thermodynamic pressure, the viscous stress tensor, the heat flux, and the Kronecker delta. For a calorically perfect gas in thermodynamic equilibrium, = ( − 1) (︀ − || || 2 /2 )︀ , where = / > 1 is the ratio of specific heats and in particular ≈ 1.4 for air. and are the specific heats at constant pressure and volume, respectively. For a Newtonian fluid with the Fourier's law of heat conduction, the viscous stress tensor and heat flux are given by where denotes temperature, the dynamic (shear) viscosity, the bulk viscosity, = / the thermal conductivity, and the Prandtl number. In particular, ≈ 0.71 for air, and additionally = 0 under the Stokes' hypothesis. The unsteady, compressible Euler equations are obtained by dropping the viscous flux in Eq. (A.1b).

Finite element mesh
We denote by ℎ a collection of stationary, non-singular, conforming, -th degree curved elements that partition 5 Ω, and set ℎ := { : ∈ ℎ } to be the collection of the boundaries of the elements in ℎ . For an element of the collection ℎ , = ∩ Ω is a boundary face if its − 1 Lebesgue measure is nonzero. For two elements + and − of ℎ , = + ∩ − is the interior face between + and − if its − 1 Lebesgue measure is nonzero. We denote by ℰ ℎ and ℰ ℎ the set of interior and boundary faces, respectively, and we define ℰ ℎ := ℰ ℎ ∪ ℰ ℎ as the union of interior and boundary faces. Note that, by definition, ℎ and ℰ ℎ are different. More precisely, an interior face is counted twice in ℎ but only once in ℰ ℎ , whereas a boundary face is counted once both in ℎ and ℰ ℎ .

Finite element spaces
Let ( ) denote the space of polynomials of degree at most on a domain ⊂ R , let 2 ( ) be the space of Lebesgue square-integrable functions on , and 0 ( ) the space of continuous functions on . Also, let denote the -th degree parametric mapping from the reference element to an element ∈ ℎ in the physical domain, and be the -th degree parametric mapping from the reference face to a face ∈ ℰ ℎ in the physical domain. We then introduce the following discontinuous finite element spaces in ℎ and the following finite element space on the mesh skeleton ℰ ℎ where ℰ E ℎ is a subset of ℰ ℎ . Different choices of ℰ E ℎ lead to different discretization methods within the hybridized DG family that have different properties in terms of accuracy, stability, and number of globally coupled unknowns [21,61]. In particular, the HDG, EDG and IEDG methods are obtained by setting ℰ E ℎ = ∅, ℰ E ℎ = ℰ ℎ and ℰ E ℎ = ℰ ℎ , respectively. We finally define several inner products associated with these finite element spaces. In particular, given , ∈ ℎ , , ∈ ℎ and , ∈ ℳ ℎ , we write where : denotes the Frobenius inner product of two matrices.

Hybridized DG discretization
The hybridized DG discretization of the (unfiltered) unsteady compressible Navier-Stokes equations reads as follows: Find (︀ ℎ ( ), ℎ ( ),̂︀ ℎ ( ) )︀ ∈ ℎ × ℎ × ℳ ℎ such that for all ( , , ) ∈ ℎ × ℎ × ℳ ℎ and all ∈ [0, ), as well as for all ∈ ℎ . Here, denotes the unit normal vector pointing outwards from the elements,̂︀ ℎ is the boundary condition term (whose precise definition depends on the type of boundary condition), and̂︀ ℎ and︀ ℎ are the inviscid and viscous numerical fluxes defined aŝ︀ ℎ (̂︀ ℎ , ℎ ) = (̂︀ ℎ ) · + (̂︀ ℎ , ℎ ; ) · ( ℎ −̂︀ ℎ ), (A.8a)︀ ℎ (̂︀ ℎ , ℎ , ℎ ) = (̂︀ ℎ , ℎ ) · . (A.8b) We note that this form of the numerical flux does not involve an explicit Riemann solver on the element faces. Instead, it is the so-called stabilization matrix ∈ R × that implicitly defines the Riemann solver in hybridized DG methods. The interested reader is referred to [22] for a discussion on the relationship between the stabilization matrix and the resulting Riemann solver. The hybridized DG discretization of the unsteady compressible Euler equations is obtained by dropping Eq. (A.7a) and the viscous terms in Equations (A.7b)−(A.7c). Finally, in the context of explicit LES, we augment the viscous stress tensor and heat flux in (A.3) with the modeling terms where is the kinematic eddy viscosity, the isotropic part of the subgrid-scale stress tensor, and the SGS eddy Prandtl number, and are computed using the subgrid-scale models discussed in Section 2. Additional details on the hybridized DG methods for compressible flows are presented in [21].