A comparative study on polynomial dealiasing and split form discontinuous Galerkin schemes for under-resolved turbulence computations

This work focuses on the accuracy and stability of high-order nodal discontinuous Galerkin (DG) methods for under-resolved turbulence computations. In particular we consider the inviscid Taylor-Green vortex (TGV) flow to analyse the implicit large eddy simulation (iLES) capabilities of DG methods at very high Reynolds numbers. The governing equations are discretised in two ways in order to suppress aliasing errors introduced into the discrete variational forms due to the under-integration of non-linear terms. The first, more straightforward way relies on consistent/over-integration, where quadrature accuracy is improved by using a larger number of integration points, consistent with the degree of the non-linearities. The second strategy, originally applied in the high-order finite difference community, relies on a split (or skew-symmetric) form of the governing equations. Different split forms are available depending on how the variables in the non-linear terms are grouped. The desired split form is then built by averaging conservative and non-conservative forms of the governing equations, although conservativity of the DG scheme is fully preserved. A preliminary analysis based on Burgers' turbulence in one spatial dimension is conducted and shows the potential of split forms in keeping the energy of higher-order polynomial modes close to the expected levels. This indicates that the favourable dealiasing properties observed from split-form approaches in more classical schemes seem to hold for DG. The remainder of the study considers a comprehensive set of (under-resolved) computations of the inviscid TGV flow and compares the accuracy and robustness of consistent/over-integration and split form discretisations based on the local Lax-Friedrichs and Roe-type Riemann solvers...


Introduction
In computational fluid dynamics (CFD), high-fidelity simulations have been recognised as a crucial area of investigation for understanding the physics underlying high Reynolds number turbulent flows, thereby providing the key tools for off-cycle aerodynamic design by means of physics-based simulations [77].However, current numerical schemes adopted in industry might not be suitable to provide an effective solution for these problems due to a highly dissipative and dispersive behaviour that damps and distorts turbulent structures, thus altering the energy transfer across scales and ultimately corrupting flow physics.An attractive class of numerical methods that is now widely spread in academia and is perhaps taking its first steps into industry is constituted by the so-called spectral element methods (SEM).Among these, one has the well-known continuous and discontinuous Galerkin (CG/DG) methods [35,43], but also more recently developed formulations such as spectral volume [88], spectral difference [55] and flux reconstruction [40] schemes.The present work focuses on the DG method, but its results may also be of special interest to those working with flux reconstruction, as the latter is able to recover certain variants of DG [17,58,61].
Spectral element methods combine the advantageous properties of finite element/finite volume and spectral methods, namely geometric flexibility and reduced dispersion/dissipation errors (see e.g.[21,59,62] on the latter topic).Nevertheless, the application of spectral element methods to challenging problems such as turbulent flows over complex geometries is still somewhat hindered by numerical instability issues.These are twofold: (a) the overall dissipation of high-order spectral element approaches might be insufficient to stabilize high Reynolds number turbulence computations, which are normally under-resolved due to computational cost constraints and therefore have significant energy at the smaller scales captured [30,67]; (b) the most efficient spectral element methods that help alleviate the cost requirements commonly rely on under-integration of the non-linear terms, see e.g.[36], which induces aliasing-driven instabilities that can also lead to numerical divergence [45].If the discretisation order is not very high, dissipative effects might overcome polynomial aliasing errors and suppress instabilities, but accuracy is often compromised.At higher orders, however, aliasing and (lack of) dissipation issues require careful consideration and sometimes robustness may not be guaranteed even with dealiasing [50,60,67].
Issue (a) above can be addressed through the use of (dissipative) subgrid-scale models, in what would be effectively a large-eddy simulation (LES) methodology.Promising candidates within this category are based on the so-called variational multiscale (VMS) approach introduced by Hughes et al. [38,39].This approach confines modelling to the smaller scales and can be naturally accommodated in a spectral element setting if one takes into account the hierarchical nature of the polynomial expansions used to represent the numerical solution [69,89].More recently, DG-based VMS approaches have received renewed attention, see [1,12,16].The effects of polynomial dealiasing on this kind of approach have been considered in [4], which showed that under-integration can partially mask the physics embedded in the VMS model employed.An alternative strategy where dissipation is also confined to the small scales consists in the addition of spectral vanishing viscosity (SVV), see e.g.[42,46,47].SVV can be understood as a higher-order viscosity, designed to affect only the highest solution wavenumbers or polynomial modes [66,71].This allows for the exponential convergence property of spectral element methods to be maintained, which helps achieving high accuracy at well-resolved laminar and transitional flow regions.Successful applications of this approach, more commonly used with CG methods, can be found in [48,56,63,74].This SVV-based eddy-resolving approach can be considered an implicit LES (iLES) methodology, where, broadly speaking, numerical stabilization techniques are relied upon to dissipate small scales in lieu of a subgrid-scale turbulence model [34].
Implicit LES approaches based on discontinuous spectral element methods have also received significant attention in recent years [3,70,82,85].This is partially because the upwind dissipation peculiar to these schemes already provides good stability for advection-dominated flows.As a result, valid solutions can be obtained without additional stabilization or modelling at moderate Reynolds numbers.The accuracy of these methods with regards to iLES is in great part due to a favourable (SVV-like) dissipative behaviour in wavenumber space, which does not affect the large scales directly and is only significant at high wavenumbers/frequencies [67,68,83,84,87].Still, as mentioned previously, this stabilization can be insufficient for high Reynolds number flows, especially at higher discretisation orders, even with over-integration.Polynomial filtering techniques can also be used for dealiasing.With the proper choice of the filter parameters the overall robustness is improved by suppressing spurious oscillations [24,30].
In this work, we focus on issue (b) above.We consider consistent/over-integration as a reference dealiasing strategy, which in principle should only improve upon the accuracy and stability of DG-based iLES, and compare it against a novel class of split-form DG discretisations developed over recent years [28,29,31,32].The latter are expected to have positive built-in dealiasing characteristics and have already demonstrated remarkable potential in terms of robustness.More specifically, in [32], certain split forms proved capable of stabilizing the inviscid Taylor-Green vortex (TGV) problem [75] even at very high polynomial orders.This test case is extremely demanding in terms of stability due to the absence of molecular viscosity and can easily diverge numerically regardless of consistent/over-integration [67].However, the key goal beyond robustness is, of course, accuracy.The present study offers a detailed analysis of a comprehensive set of inviscid TGV test cases, and demonstrates that split form DG approaches can offer results of quality very similar to that achieved by consistent/over-integration.Although split forms are formulated by averaging standard conservative and advective forms of the governing equations, previous results derived in the context of finite difference methods [22] can be invoked to guarantee that certain nodal split form DG discretisations remain conservative.From a numerical perspective, therefore, this approach is believed to aggregate significant advantages for eddy-resolving computations of high Reynolds number compressible flows.However, the main goal of this study is not to advocate DG-based iLES approaches.Instead the goal is, for the first time, to investigate key aspects of solution quality and robustness the two dealiasing approaches for nodal DG methods might exhibit when viscous effects are negligibly small.This paper is organised as follows.Section 2 presents a brief overview of a standard nodal DG formulation along with the two dealiasing strategies under consideration.Next, Section 3 provides a preliminary analysis on the built-in dealiasing mechanisms of relevant split form DG methods in the context of Burgers' turbulence.In Section 4 we introduce the inviscid TGV test cases and discuss the stability of the two dealiasing techniques considered.Section 5 directly compares the solution quality of the two approaches for various polynomial orders and different Riemann solvers.Concluding remarks are given in Section 6.

Discontinuous Galerkin discretisations
This section introduces the different DG formulations considered in this study.We first briefly review the standard and over-integrated DGSEM approaches in Section 2.1.We then discuss the split form DG discretisations which are relevant to our work in Section 2.2.Lastly, in Section 2.3 we describe in more detail the numerical interface fluxes used with the different formulations being considered.

Standard and over-integrated DGSEM approaches
We consider a system of equations given by conservation laws in the form where q(x, y, z, t) is the vector of conserved variables and f (q), g(q) and h(q) are the flux vectors governing the transport of q in a physical domain Ω.For the spectral element approximation, we subdivide Ω into non-overlapping elements Ω e such that e Ω e is a mesh of the computational domain.For simplicity, we restrict the following presentation to straight-sided hexahedral elements, but remark that the extension to more general curvilinear meshes is available, see e.g.[49].
In the nodal DGSEM approach considered, the solution is approximated within each element Ω e through a polynomial expansion such as in which the summation index is assumed to span a polynomial space of degree N , and where time-dependent coefficients qn (t) are associated to the nodal basis functions ϕ n (x).The spatial dependence usually relies on a simple trilinear mapping relation x = x(ξ), cf.[32], that connects the physical space Ω e to the standard domain Ω 0 = [−1, 1] 3 .This allows for the basis functions to be defined with regards to Ω 0 directly, namely where, for example, index i = i(n) references the N -th order Lagrange polynomial i (ξ), which is based on an appropriate set of nodes {ξ l } N l=0 defined within [−1, 1].The nodal DG methods in this work uses a collocation approach, e.g.[49], where interpolation and quadrature nodes are selected to be the Legendre-Gauss-Lobatto (LGL) nodes.The choice of LGL quadrature is important, because it ensures that the discrete DG derivative matrix and the discrete mass matrix satisfies the summation-by-parts (SBP) property for any polynomial order [28].The SBP property is crucial to define a split form DGSEM that remains conservative [22,23].Finally, we note the cardinal property of Lagrange polynomials (2.4), whereby i (ξ l ) = δ il , with δ il = 1 for i = l and δ il = 0 otherwise.
We then require that (2.1), when evaluated with (2.2), should vanish locally in a projection sense, i.e.
for all the basis functions, ϕ m , used in (2.2).Having in mind the simulations to be discussed in this study, we now restrict our discretisation to the case of Cartesian meshes, where element-wise Jacobian factors and metric terms simplify considerably to J = 1 8 ∆x∆y∆z, x ξ = 1 2 ∆x, y η = 1 2 ∆y, z ζ = 1 2 ∆z, in which ∆x, ∆y, and ∆z represent the element side lengths along the three Cartesian directions.In this case, casting (2.5) as an integration over Ω 0 yields [49] Ω0 where F = [F, G, H] T is a vector containing the contravariant fluxes, which incorporate the element-wise constant metric terms, namely Finally, we integrate (2.6) by parts and replace the boundary terms stemming from the divergence theorem with numerical interface fluxes, F * .The resulting statement for the weak DG formulation reads in which the differential surface vector dS points toward the outside of Ω 0 .The DG schemes rely on well-known Riemann solvers [79] for F * to resolve discontinuities in the approximation at element surfaces.
There is additional flexibility for a split form DGSEM which can use special numerical fluxes to recover certain secondary properties like kinetic energy preservation [32], as will be discussed in Section 2.3.Each integral in (2.8) is evaluated numerically with a Gauss type quadrature, given the number Q of integration nodes, see e.g.[43].The standard DGSEM relies on collocation with Q = N + 1 integration points per direction, usually LGL nodes [35,49].While this choice only guarantees exact integration of flux functions depending linearly on the solution, the same choice is adopted, in many cases, for problems involving non-linear flux functions.This collocated approximation is advantageous in terms of computational cost, cf.[36], but can also be numerically unstable and less accurate due to polynomial aliasing errors introduced by the under-integration of non-linear flux functions, especially for under-resolved computations [30,45,60].The straightforward way to alleviate this problem is to simply use a larger number of quadrature points, i.e.Q > N + 1, in what would be an over-integrated DGSEM.
The over-integration mentioned above becomes a consistent integration when the quadrature order is chosen to be consistent with the nonlinearity of flux functions.For example, one typically requires Q ≈ 3 2 (N + 1) for quadratic and Q ≈ 2(N + 1) for cubic nonlinearities.These are respectively common choices for the consistent integration (of the advective terms) of the incompressible and compressible Navier-Stokes equations [43].For compressible flows with small density variations, using Q ≈ 3 2 (N + 1) might also reproduce the effects of consistent integration [3,4].However, we remark that Gauss type quadratures can only guarantee exact results for the integration of polynomial functions, and that compressible flow equations (written in the usual conservative variables) have flux vectors whose components are actually rational functions.As a result, strictly speaking, the associated quadratures are not exact even for Q ≈ 2(N + 1).As asymptotic convergence to exact results is nevertheless guaranteed [80], quadrature errors can be expected to decay numerically to negligible levels for Q sufficiently large, as exemplified in Section 4.
The standard DGSEM is perhaps more commonly known in the so-called strong form variant of (2.8), which is obtained with an additional integration by parts and reads which in general differs from (2.8) when quadratures are inexact.For standard DGSEM, however, it is possible to show that forms (2.8) and (2.9) are actually discretely equivalent [51], despite any quadrature errors.Also, due to the cardinal property of the Lagrange polynomials employed and the collocation of quadrature and solution nodes, the mass matrix stemming from the leftmost integral above has a diagonal structure and is given by [35,49] where {ω l } N l=0 is the set of quadrature weights associated to the set of LGL nodes {ξ l } N l=0 adopted.We note that when consistent integration is employed, mass matrices of nodal approaches are actually full.This is what causes dispersion-diffusion characteristics of standard and over-integrated DGSEM approaches to differ at higher wavenumbers [35].
In the case of standard DGSEM, the diagonal mass matrix allow for the semi-discrete evolution equation of each nodal coefficient to be written independently, in a form that resembles (2.1), namely in which the modified (tilde) fluxes above incorporate the interface contributions and are given by where D is the standard polynomial derivative matrix, see e.g.[49], defined as for m, n = 0, . . ., N. (2.15) In the following section, different split form approaches will be derived based on simple modifications of the rightmost terms in (2.12)-(2.14).

Split form DG discretisations
Here we discuss the so-called split form approaches to the discretisation of the advective terms governing general fluid flow.The structure of any split formulations is equation dependent because knowledge of the type of non-linearity (e.g.integer powers of the fluid velocity) is necessary to create an average of the conservative and advective forms of the advective terms [20,32,76].In this work we consider the compressible Euler equations, most commonly stated in conservation form as (2.1) with q = [ρ, ρu, ρv, ρw, ρe] T and where ρ e = p/(γ − 1) + ρ (u 2 + v 2 + w 2 )/2 while ρ, u, v, w, p, e and γ = 7/5 stand respectively for density, the three velocity components along the Cartesian directions, pressure, specific total energy and the usual ratio of specific heats.Split formulations are alternative discretisations of the non-linear transport terms typical of many partial differential equations (PDEs).Essentially, splitting techniques are built by averaging conservative and advective forms of a PDE.Relevant split forms usually improve upon the robustness of numerical methods by reducing aliasing errors.Split formulations have been often used in conjunction with high-order finite difference schemes [20,44,52,54,64] and spectral methods [5,91].There are many ways to rewrite the non-linear advective terms of a PDE.A good overview of different split form approaches can be found in [72].Depending on how one interprets the non-linearity of the Euler fluxes (quadratic, cubic or rational) there are several ways to rewrite the equations in an equivalent split form.To simplify the discussion, we introduce notation for specific instances of the quadratic splitting and of the cubic splitting explored by Kennedy and Gruber [44] (a Some splitting strategies are motivated by the type of flow physics that the numerics are expected to capture.For example, Ducros et al. [20] considered flows with small density variations and used the quadratic splitting (2.17) for each term in the Euler fluxes.Alternatively, Kennedy and Gruber [44] wanted to account for density variations and hence used both (2.17) and (2.18).These two splitting strategies, hereafter denoted respectively as DU and KG, are given, e.g. in the x-direction, by and or, by using the compact notation introduced in (2.17) and (2.18), (2.21) We note that while the DU splitting does not lead to discretisations that are formally kinetic energy preserving [72], the KG splitting does allow for that possibility [41].
With the concept of split forms, it is possible to build alternative DGSEM discretisations with special properties guaranteed at the discrete level, such as kinetic energy preservation or entropy consistency [28,29,31,32].However, because these are constructed by averaging conservative and non-conservative forms of the governing equations, it is important that the resulting DG scheme remains locally conservative.This can be guaranteed for the forms considered here by the following.We start by modifying the computation of the volume integral terms of the standard DGSEM in its strong variant, cf.(2.9) and (2.12)-(2.14).The collocated mass (2.10) and derivative (2.15) matrices are combined into the matrix Q := MD, which has the summation-by-parts (SBP) property [9], namely that is used to mimic integration-by-parts at a discrete level, by manipulating the derivative matrix as A remarkable result presented in Fisher et al. [23] and Fisher and Carpenter [22] showed that diagonal norm SBP matrices such as D can be reinterpreted as sub-cell finite volume type differencing operators.This result is what guarantees local conservation of diagonal norm SBP discretisations in the sense of Lax-Wendroff due to a telescoping flux differencing [22].For more details and proofs, the interested reader is referred to [32].
As further explored in [32], the flux differencing formulation mentioned above can be recast into a matrix type product form by the introduction of auxiliary sub-cell numerical flux functions F # , G # , and H # , which are required to be symmetric and consistent, e.g. (2.24) Finally, this matrix type product form based on sub-cell numerical fluxes can be linked (via transitivity) to the original element-wise differentiation operation, cf.(2.12)-(2.14),such that where the original differentiation is recovered for fully central sub-cell numerical fluxes, e.g.F # (q ijk , q njk ) = 1 2 [F(q ijk ) + F(q njk )], but also different split forms can be implemented effortlessly with other choices.For instance, if we use instead the product of two averages or the product of three averages for the sub-cell numerical flux we recover discrete versions of the quadratic (2.17) and cubic (2.18) split forms, see [32] for complete details.This creates a dictionary where one can easily translate between the continuous split forms and their discrete counterparts.For example, in the x-direction, the split forms in (2.21) can be recovered respectively through the sub-cell numerical fluxes F # DU (q ijk , q njk ) and F # KG (q ijk , q njk ) given by where {{•}} represents the arithmetic mean -in this case, between the relevant quantities from q ijk and q njk .The dot splitting notation (2.17) and(2.18)elucidates the usage of symmetric sub-cell fluxes: one simply replaces the dots with products of arithmetic means.Reference [32] discusses many other sub-cell flux functions, e.g.those related to other well-known split forms like that of Morinishi [64] or Pirozzoli [72].Lastly, we replace the rightmost terms in (2.12)-(2.14)with (2.25)-(2.27)and, following [32], we connect the choice of the sub-cell flux to that of the interface numerical flux.For example, if the KG splitting is adopted, the local Lax-Friedrichs (LLF) interface flux function will have, e.g. for the x-direction, the adapted form where λ max is an estimate of the fastest wave speed (in absolute value) at the interface between the left and right solution states q L and q R .More details about the interface flux functions employed in the different DGSEM approaches under consideration will be discussed in the next section.

Selection of the numerical interface flux function
We briefly outline the two types of interface flux functions that will be used with the over-integrated and split form DGSEM approaches considered.For over-integrated DGSEM, we use the classical forms of LLF and Roe-type schemes [79].For example, in the x-direction, the LLF flux formula reads where terms are denoted analogously to (2.29).As for the classical Roe solver, the dissipative (rightmost) term relies on a characteristic decomposition of the flux Jacobian evaluated at Roe's average state [73].
In contrast, for split form DGSEM, we couple the choices of sub-cell and interface numerical fluxes, as shown for the LLF case in (2.29).When the Roe flux is employed with the DU splitting, the symmetric part of Roe's flux formula is replaced with F # DU (q L , q R ) and the original dissipation term is maintained.In case the KG splitting is adopted, however, a slight modification is made to the dissipation term.We note that the specific Roe averaging remains identical to the classical scheme, but the dissipation term is modified to ensure that the KG scheme is kinetic energy stable, see [11,32,41].To do so, we alter the diagonal matrix of eigenvalues so that the first eigenvalue matches the last one [11], namely where a represents the speed of sound.
Finally, for the split form schemes we note a somewhat remarkable property.It is possible, in some cases, to obtain numerical solutions of the compressible Euler equations which nearly conserve the total kinetic energy [32].This can be achieved, for example, with a "fully central" interface flux such as which lacks the usual interface dissipation term.Simulations based on this central flux will be discussed in Section 5.2, where solutions obtained with very high polynomial orders will be compared to lower order solutions that conserve total kinetic energy.All simulation results in Secs. 4 and 5 are integrated in time with an explicit five stage, fourth order accurate low storage Runge-Kutta scheme [10], where a stable time step is computed according to an adjustable coefficient CFL∈ (0, 1], the local maximum wave speed and the relative grid size, e.g.[27].

Analysis of polynomial aliasing errors of split form approaches
We now undertake a preliminary assessment of how inexact quadratures can affect the energy of different solution coefficients and how the split form approach tends to keep this energy "under control," thus increasing numerical robustness.The objective is to explain, at least partially, why certain split form DG approximations have positive stabilisation and dealiasing properties.A similar analysis has been carried out for pseudo-spectral methods (in Fourier space) by Blaisdell et al. in [5], where it was demonstrated that split form discretisations suppress aliasing errors as if they had a built-in dealiasing mechanism.Here we simply illustrate through numerical experiments that this convenient feature seems to hold for nodal DG-based spectral element methods.More specifically, we consider the Burgers' equation in one spatial dimension, given by (2.1) with flux function f (q) = q 2 /2.Similar numerical experiments have been used in previous works [45,47] to illustrate the effect of different dealiasing techniques in controlling the solution's modal energy.
The inner-product formulation of the strong split form version of the Burgers' equation reads [28] where (•) denotes differentiation in space and •, • Q stands for the quadrature of the product between two functions using Q integration points (LGL nodes).Furthermore, test functions i are the usual Lagrange polynomials of degree N , cf. (2.4), and f * is the numerical flux required at the right and left elemental interfaces, denoted above as ⊕ and , respectively.The parameter α ∈ [0, 1] is used to introduce biased averages between conservative (α = 1) and non-conservative (α = 0) discretisations.
In order to focus on the terms requiring quadratures, we work instead with the weak split form version of the Burgers' equation, where boundary contributions can be more easily dismissed.In the context of DGSEM discretisations, this version is entirely equivalent to (3.1), see [28].If we disregard boundary contributions by setting f * = 0, this weak form becomes Note that in both formulations above, numerical differentiation of quadratic terms like (q 2 ) and (q i ) is performed in accordance with the collocated nodal DGSEM framework: arguments are evaluated via point-wise product at the Q integration points and treated as Lagrange polynomials of degree N = Q − 1 (instead of 2N ), whose derivative is then taken in a standard way.
Although our focus is on the characteristics of collocated nodal DGSEM, the assessment of solution coefficients in modal space is particularly insightful as it allows for interpretation by analogy with Fourier space, although it is recognized that this analogy is somewhat limited.This will require the appropriate transformations between nodal and modal spaces to be employed.In the following, we restrict ourselves to the reference domain [−1, 1] and consider a "frozen" numerical solution defined by the modal expansion where L j (ξ) is the orthonormal Legendre polynomial of degree j.By the Fourier/modal space analogy, the above expansion is representative of a turbulent velocity field, since q2 j ∝ (j + 1) −5/3 as typical of turbulence energy spectra.We note that such scaling is achievable even with the Burgers' equation by means of a specific forcing [68,92].The coefficients of a nodal expansion equivalent to that of (3.3) are given by q = Vq, where V ij = L j (ξ i ) is the Vandermonde matrix evaluated at the (N + 1) LGL nodes ξ i .
The quadratures in (3.2) are evaluated with the same set of Q = N + 1 points ξ i mentioned above, in accord with the collocated nodal DGSEM approach, as follows: where ∂ ξ (•) = ∂(•)/∂ξ, ω k are the quadrature weights and I LGL is the LGL interpolant of the product q δ i , i.e.I LGL (ξ) is a Lagrange polynomial of degree N defined as q δ (ξ k ) i (ξ k ) at the N + 1 quadrature nodes.Note that q δ (ξ) i (ξ) is actually a polynomial of degree 2N , and so the approximation associated with the interpolant I LGL is presumably related to the positive dealiasing properties observed in relevant split form DGSEM discretisations.
The split form obtained from (3.2) with α = 1/2 is of particular interest because it is analogous to the quadratic split (2.17) used to create stable split forms for the compressible Euler equations.This form is considered here for the Burgers' equation in an attempt to anticipate the behaviour of relevant split form DGSEM as applied to Euler-based turbulence, which will be addressed in the remaining sections of the present work through the inviscid Taylor-Green vortex problem.The conservative and non-conservative forms of (3.2) are also considered due to their obvious significance.These three forms are compared in Fig. 1, which shows the transformed (modal space) values of the right-hand side of (3.2), namely TRHS i ≈ ∂ t qi = V −1 ∂ t q, for different modal coefficients i = 0, . . ., N .
We stress that TRHS values do not take into account boundary contributions, as they are intended to represent the effect of the volumetric terms on the numerical solution.Also, for the one-dimensional test addressed, the interface flux contribution does not depend on Q or even α.The expansion order, N , varies for the different plots in Fig. 1.The exact values of TRHS i have been computed via consistent integration of the conservative term in (3.2) and are also shown in the plots for reference.We note that even higher values of N were considered and the main trends remained the same.Random perturbations of up to 50% in the coefficients of (3.3) have been added, and again the main tendencies observed in Fig. 1 remained.This demonstrates that the results obtained are not caused by a fortuitous, particular choice of parameters.Note also that the overall profile shown in Fig. 1 is physically consistent in the sense that energy is flowing from low-order modes (low frequencies) to the high-order ones (high frequencies), as happens in the energy cascade mechanism typical of turbulence.The most relevant observation one can draw from Fig. 1 is that conservative DGSEM tends to over-predict (intensify) the coefficients' variation rate TRHS i ≈ ∂ t qi whereas non-conservative DGSEM under-predicts it, this being especially significant at the highest-order modes.In terms of modal energy, one can expect these results to be even more significant for the relative variation rates (∂ t q2 i )/q 2 i , which, in the frozen solution analysis considered, should scale like In other words, higher-order modes suffer a larger relative deviation (over/under-prediction) in terms of modal energy variation rate since they are inherently less energetic.As a result, standard conservative and non-conservative DGSEM approaches without over-integration are expected to grossly miscalculate the dynamics of small turbulent scales, eventually inducing numerical instabilities.The over-prediction of modal energy at higher modes in conservative DGSEM approaches has been reported in previous studies [45] and correctly attributed to polynomial aliasing errors due to inexact integration of the non-linear terms.Here we have shown that the quadratic split form tends to predict energy variations more accurately due to the averaging of opposite conservative and non-conservative tendencies.In a sense, this can be regarded as a built-in dealiasing mechanism as it suppresses aliasing errors.This result is similar to that reported for pseudo-spectral methods [5], although the reasons are possibly different because aliasing errors here originate from inexact polynomial quadratures, whereas in [5] the aliasing errors stem from the handling of Fourier transforms with an inconsistent number of modes.
Another relevant point is that energy variation levels yielded here by split form DG approaches tend to be slightly below the correct levels, which can be advantageous for robustness and should help to keep energy "under control," suppressing excessive growth of small-scale energy and preventing numerical divergence in under-resolved turbulence simulations.It has also been found (not shown) that when the number of quadrature nodes Q is increased for the conservative form, energy variation levels approach the exact values monotonically from above, as over-integration gradually reduces the TRHS levels over-predicted by the conservative form and bring them to the correct values.This qualitatively explains why split forms can be more stable than conservative discretisations even when quadratures are consistent (with non-linearities of the problem) albeit still inexact.This is in line with a recent study [57] that proved over-integration unable to entirely eliminate aliasing errors in advection problems with non-constant speed functions, whereas appropriate split forms were able to.
We recall that when the compressible Euler equations are written using conserved variables, the flux functions to be integrated numerically involve rational functions which, in general, cannot be integrated exactly with a finite number Q of quadrature nodes.Integrations are however expected to converge to the exact values as Q increases, which could grant higher accuracy for over-integration approaches at high enough Q, if they are stable (and affordable).The accuracy and robustness of split and over-integrated conservative forms as applied to under-resolved turbulence are assessed and compared in the following sections.

TGV simulations and numerical stability
The Taylor-Green vortex flow was introduced in [78] as a model problem for the analysis of transition and turbulence decay.The test case was originally proposed for the incompressible Navier-Stokes equations in a cubic domain with (triply-)periodic boundary conditions.As in previous studies [19,75], we adopt a modified version of the initial conditions, which is suited for compressible flow solvers.The following expressions have been used as the initial state within Ω = [−π o , π o ] 3 , respectively, for the density, the three velocity components, and the static pressure: the total energy per unit volume being E = p/(γ − 1) + ρ(u 2 + v 2 + w 2 )/2.For the reference quantities, we have adopted the values o = ρ o = V o = 1 and c o = 10, leading to a Mach number of 0.1, which makes our cases nearly incompressible.A non-dimensional time t is adopted based on the scale o /V o = 1.A Reynolds number could be defined as Re = ρ o V o o /µ o , but only the inviscid problem is considered in this study, where the compressible Euler equations have been simulated directly, with γ = 1.4.Similar to Gassner et al. [32] we take the final time of the TGV flow evolution to be T = 14.The Euler-based simulations considered in the present study are expected to be representative of viscous TGV solutions at very high Reynolds numbers, provided that numerical dissipation mimics the dissipative character of Navier-Stokes turbulence in the limit of vanishing viscosity [25].The equivalent scenario in traditional LES would be to have zero molecular viscosity and rely solely on the regularization of a subgridscale model.Here, upwind dissipation is expected to play this latter role.In [67], the inviscid TGV was simulated with different Riemann solvers in the context of DGSEM with consistent integration.The study showed in particular that Roe-based discretisations were more robust than LLF-based ones, which crashed more easily -a surprising result in itself, which will be discussed below.However, at higher polynomial orders (N ≥ 6) both fluxes produced test cases that lacked stability.Subsequent DGSEM simulations [32] based on suitable split form discretisations demonstrated remarkable robustness and were capable of yielding stable solutions of the inviscid TGV with both Roe and LLF-type fluxes, even at very high-orders (e.g.N = 15).In the following we summarize the test cases addressed in [32,67], consider new (lower-order) test cases, and discuss the inherent instability issue of the inviscid TGV flow.
The base set of test cases considered are given in Table 1.Each column corresponds to the number of one-dimensional polynomial modes m = N + 1 used, N being the polynomial order.The number of elements varies by a factor of √ 2 between adjacent rows to reduce computational costs, as choosing a factor of 2 would requires a large number of degrees of freedom in the last line.The values of n el were chosen so that the degrees of freedom N dof = (n el m) 3 are kept (approximately) constant along a given row.Equispaced grids of cubic elements have been employed.Values in the Table's core represent the number of elements in each spatial direction, n el .Crossed out numbers denote simulations that crashed with consistent integration.The KG and DU split form discretisations were however able to stabilize all the tests cases, see [32] for a complete discussion on the robustness of other available split forms.Spectral element codes N ektar++ [8] and F LU XO (http://www.github.com/project-fluxo)have been used respectively for the computations based on consistent integration and split forms.As pointed out in [67], the crashes are probably not related to insufficient integration.We acknowledge that exact integration is never achieved with a finite number of quadrature points since the flux vector components are not polynomials, but rational functions, given that our consistent integration approach is based on the conservative form of the Euler equations.However, since quadrature errors still tend to zero as the number of integration points is increased [80], they are expected here to be negligible in the case of consistent integration, given the low Mach number adopted.In Table 2, we consider the time of crash, t c , versus the number of integrations points Q employed (per element and per dimension) for test case m = 8, n el = 14.Consistent integration is performed both for the volume and surface terms, which take Q 3 and Q 2 points, respectively.We note that this scaling is what makes consistent integration prohibitive at high polynomial orders.The asymptotic behaviour of t c for LLF (particularly as Q is increased from 16 to 32) indicates that consistent/over-integration has already contributed all it could towards stabilisation.Further, we note that increasing the value of Q does offer stabilisation for Roe and entries in Table 2 marked with "-" ran successfully.
For Legendre-Gauss-Lobatto quadratures, the number of nodes required for the consistent integration of linear, quadratic and cubic terms are respectively Q ≈ m, Q ≈ 3m/2 and Q ≈ 2m.We see from Table 2 that LLF's t c increases with Q until up to Q = 3m/2 and then remains practically unaffected.The same is true for the Roe flux, except that Roe-based cases are stabilized for Q ≥ 3m/2.This is consistent with the nearly incompressible nature of the TGV cases considered, since density variations are small and the terms being integrated are essentially quadratic.Despite of this, consistent integration of all the cases in Table 1 assumed a cubic non-linearity for the compressible Euler equations and used Q = 2m.The differences between LLF and Roe are likely due to strong over-upwinding effects induced by the former, which increase the likelihood of TGV instabilities, as explained in the following.
There is an open debate in the literature as to whether or not the inviscid TGV flow might develop singularities which lead to the actual collapse of the solution [7,15,33,37].However, as emphasized in [67], this possibility is only considered for the exact, energy conserving solution of the Euler equations.This is in contrast with the character of Navier-Stokes turbulence in the limit of vanishing viscosity, which remains dissipative and is not expected to develop singularities, see [25] Secs.5.2 and 9.3.We stress that any LES-like approach should follow the latter behaviour when viscosity is set to zero as subgrid dissipation (explicit or implicit) remains active.Nevertheless, an energy-conserving bias has been found to be partially induced by (some of) the discretisations considered due to the very sharp spectral dissipation expected at higher orders [68] and also from over-upwinding [65].It is known that a sharp dissipation in Fourier space can induce an energy-conserving behaviour in turbulent simulations [2].This is because, under certain circumstances, a cutoff-like spectral dissipation may prevent the formation and subsequent destruction of small scales near the cutoff wavenumber, acting effectively as a barrier to the inertial cascade and causing a pile-up of small-scale energy.The LLF solver is expected to have a particularly sharp dissipation due to its over-upwind bias for the momentum equations, owing to the disparity between acoustic and convective wave speeds, especially at low Mach numbers, see [65].As explained in [67], the main evidence that this behaviour is taking place (even if partially) is the so-called "energy bump" observed in the energy spectra of LLF-based DGSEM simulations, see e.g., [18,90].It has been demonstrated [2,14,26] that pre-dissipative bumps emerge as the solution begins to follow an energy-conserving dynamics when only a finite number of Fourier modes are retained (limit of increasingly sharp dissipation).More evidence of this artificial energy-conserving character will be given in the next section.
Besides the occurrence of the energy-conserving bias mentioned above, it's been suggested in [67] that simulations following this character (higher-order discretisations, especially LLF-based ones) are also inducing the emergence of the physical singularities long conjectured for the inviscid energy-conserving TGV problem.This second point is however somewhat speculative and may be very difficult to prove.Still, note that while suitable split form discretisations are able to prevent instabilities with both Roe and LLF solvers even at very high orders, fully central split forms which are nearly kinetic energy conserving are only stable at low polynomial orders (N ≤ 3).This further supports our claim that the TGV instabilities observed are likely of physical origin, as their likelihood of emergence increases in higher-order (better resolved) discretisations of energy-conserving character.
In summary, although the instabilities observed are not entirely understood at this point, it is clear that suitable split form discretisations have superior non-linear stability characteristics.Robustness, however, can be easily obtained at the cost of solution quality, and so the accuracy of split form DG approaches require further investigation.Note that such solutions are still prone to aliasing errors because they do not rely on standard consistent/over-integration.The evaluation of the accuracy and fidelity of suitable split form DG approaches by comparison to consistent/over-integration for the invisicid TGV flow is the main focus of the next section.To the authors' knowledge this is the first comparison between these two dealiasing approaches to be conducted in the context of under-resolved turbulence computations.observed as a general trend in our comparisons and is consistent with the estimates of Section 3. The close proximity between results based on split forms and consistent integration have been found to hold also at higher orders (cf.Fig. 4(a)).However, for low-order discretisations (m ≤ 3), results exhibited non-negligible differences, which is also consistent with Section 3. Low-order solutions will be discussed in Section 5.2, along with some results obtained at very high-orders.The results discussed so far indicate that high-order split form discretisations can provide solutions of quality very similar to that of consistent integration.An additional assessment of solution quality for both discretisations is presented in Fig. 3, namely, the QR diagrams of the solutions whose spectra have been considered in Fig. 2.These consist of joint PDFs of invariants Q and R of the velocity gradient tensor, which is evaluated at a large set of points within the flow domain.The resulting diagram provides specific insight about the kinematic topology of the flow [13].More details about QR diagrams and their construction is given in Appendix B. Typical QR diagrams exhibit a characteristic teardrop profile, as seen, e.g., in Fig. 3(b), in various turbulent scenarios, which is considered a qualitatively universal feature of turbulence [81].In particular, a prolonged extension of the profile over the bottom-right quadrant of the diagram, known as the "Viellefosse tail" [86], is expected.The canonical shape for a teardrop profile of isotropic turbulence can be found in e.g.[53].More generally, the overall profile shape observed in a QR diagram can reveal how closely numerical results are able to reproduce physical turbulent features.Therefore, it is a good indicator to assess the quality of the simulations run for the two dealiasing strategies under consideration.It is also useful to check how the lack of resolution or the polynomial order of the approximation can affect the QR diagram.
An interesting result is shown in Fig. 3, namely, that Roe-based QR profiles very closely resemble the canonical teardrop shape [53], whereas LLF-based ones do not conform so well to that shape and feature a much less pronounced tail.This further supports that LLF-based solutions are less physical than Roe-based ones and indicates that energy bumps contribute to the overall turbulent dynamics with kinematic states that are artificial.In fact, because LLF-based QR profiles are more symmetrically oriented with regards to the origin of the diagrams, bump-related kinematic states are probably somewhat random, while LLF-based QR profiles even resemble an artificially generated (Gaussian) turbulent state considered in [81].
Assuming that energy bumps are in fact emerging from an energy-conserving dynamics (as discussed in Section 4), theory predicts that bump-related scales should reach an equilibrium (loosely called a "thermalised" state) where energy equipartition is favoured [26], which is consistent with the more symmetrical distribution of kinematic states of LLF-based QR profiles.This hypothesis will receive further confirmation in Section 5.2, where very high-orders (with sharper dissipative behaviours) are considered.A final remark from Fig. 3 is that QR profiles of split forms are only slightly wider than those of consistent integration, but otherwise practically indistinguishable.This supports our claim that high-order split form discretisations provide solutions of quality very similar to that of consistent/over-integration for the inviscid TGV flow.

Solution quality at lower and higher polynomial orders
We start by considering Roe-based solutions at higher orders, namely m = 8 and m = 16 in Fig 4 .Note that these have the same number of DOFs and that the spectrum of the equivalent test case at m = 5 can be seen in Appendix A, cf.Fig. A.8(b).At m = 16, consistent integration is not sufficient to suppress TGV instabilities and hence only the result from the KG scheme is shown here.One can see that an energy bump can emerge even with Roe-type numerical flux as the polynomial order is increased.This is consistent with our claim that energy bumps are caused by a sharper dissipation in wavenumber space.Also, it seems from Fig. 4 that the higher order solutions only achievable with split forms follow the trends that would have been obtained with consistent/over-integration if these simulations were stable.
To confirm that the very high-order limit is representative of an energy conserving dynamics, we compare in Fig. 5 the QR profiles of case m = 16, n el = 7 with that of the kinetic energy preserving KG split form without interface dissipation, cf.(2.32), which is stable for m = 4.It is possible that such a stable solution is only available because of non-negliable dispersion errors that might be staggering the localisation of the thin shear layers traditionally regarded as connected to inviscid TGV crashes [7].As previously noted, the stability of simulations that use split forms without interface dissipation is delicate and crash for m > 4, see [32].In any case, as one can see from Fig. 5, the very high order Roe-based QR profile is symmetric and very similar to that of the solution computed with a "fully central" spilt form DG scheme without additional dissipation at interfaces.Clearly, both solutions are highly affected by over-energetic bump-related scales which favour equipartition of energy and a symmetrical distribution of kinematic states.The QR profile of the LLF-based solution for case m = 16, n el = 7 (not shown) is similar to those of Fig. 5.
Finally we consider some results obtained with low-order discretisations, namely, m ≤ 3.In Fig. 6, we show the m = 3 energy spectra obtained both for split forms and consistent integration.As one can see, non-negligible differences are present.We stress that LLF results seem bump-free, but actually the additional small-scale energy of LLF spectra (when compared to Roe spectra) might be related to a mild energy bump.This is supported by the evolution of the Roe and LLF energy spectra at later times (not shown), when a    small-scale spectral region seems to persist with notably high energy in the LLF-based spectrum during the third phase of the TGV flow (nearly homogeneous decay).Additional insight on the quality of low-order cases can be gained from QR diagrams, which are given in Fig. 7, now, for m = 2.These confirm that, at sufficiently low orders, split forms and consistent integration results can differ.More importantly, they show that the turbulent kinematics of low-order DGSEM solutions   are not as clean and accurate when compared to solution at (moderately) higher orders for the same DOFs.However, the overall shape of the low-order QR profiles is still reasonably correct, and these seem even more physical than those obtained by discretisations of very high order.

Conclusions
In this work we described two strategies to dealias the computation of under-resolved turbulent flows using the DGSEM framework.The first was the well-known technique of consistent/over-integration (sometimes referred to as polynomial dealiasing) where the approximation of variational forms in the DG framework are enriched with additional quadrature points to remove aliasing errors associated with non-linear terms.The other dealiasing strategy was to reformulate the PDE into an equivalent, at the continuous level, split form expression that employs an average of the conservative and non-conservative forms of the equation.Relevant DGSEM split forms can however be shown to remain conservative.It is important to note that both consistent/over-integration and the split form DG schemes require additional computational effort compared to the standard DGSEM.The exploration of optimising implementations of each dealiasing strategy is a subject of our ongoing research.
A preliminary investigation about the built-in dealiasing mechanism of DGSEM split forms has been conducted in one dimension through a frozen Burgers' turbulence scenario.This analysis showed that, without consistent/over-integration, standard conservative DG formulations tend to overestimate the energy content of high-order polynomial modes, whereas advective formulations tend to under-estimate them.Due to the averaging of these two opposite tendencies, the relevant quadratic split form DGSEM was shown to balance aliasing errors and to keep energy content slightly below the correct levels.This latter feature is believed to also help improve robustness in under-resolved computations by keeping the growth of energy levels "under control".
Subsequently, the robustness of DGSEM approaches based on consistent integration and split forms for the inviscid TGV problem was discussed in detail.This brought together results from two recent works [32,67] as well as new solutions obtained with lower polynomial orders and clarified how the observed TGV instabilities can be either induced or suppressed by the underlying numerics.Although TGV instabilities are not entirely  understood at this point, it was argued that discretisations that promote an energy-conserving bias for the TGV flow are prone to develop instabilities.In any case, relevant split forms clearly display superior robustness and only lack stability when central fluxes are used in place of a Riemann solver.Still, the novel part of the TGV analysis concerned the quality of the turbulent solutions obtained with stable split forms.This assessment was conducted by comparing these to solutions obtained with consistent/over-integration, which, as the standard dealiasing technique, is believed to yield the best accuracy DG-based model-free turbulence computations can offer.
In the accuracy assessment of consistent integration and split forms for the inviscid TGV flow, kinetic energy spectra and QR diagrams have been analysed.At moderately high orders (i.e.around sixth order), the performance of relevant split forms was shown to be very close to that obtained with consistent integration.In this case, energy spectra and QR diagrams both exhibited the expected physical trends.At lower orders (i.e.third order or below), however, the two DGSEM approaches exhibited non-negligible differences and inferior solution quality (for the same number of DOFs).Most surprisingly, at very high orders (i.e.eighth order or above), stable solutions showed unphysical features attributed to the energy-conserving bias induced by a sharper dissipative behaviour in wavenumber space.The latter has been found to be especially prominent for the local Lax-Friedrichs flux, owing to its over-upwind character at low Mach numbers.Although this energy-conserving bias had already been conjectured in a previous short note [65], new and stronger evidence for it was presented here.In summary, DG-based model-free turbulence computations at very high Reynolds numbers are more likely to yield stable and accurate solutions when complete Riemann solvers (such as Roe's flux) are employed at moderately high orders.For challenging simulations, the use of dealiasing strategies is of key importance.This work highlighted the role of high-order split form discretisations as a robust alternative compared to consistent/over-integration.More importantly, it showed that high-order solutions from both strategies were very similar in terms of solution quality.An interesting theme left for subsequent studies concerns the performance of split form DGSEM for cases involving curved elements, whose type of aliasing error, i.e. geometric aliasing, has a different source than that of under-resolved turbulence.While these two forms are expected to differ for flows with stronger compressibility effects [32], the results in Fig. A.8 encouraged us to consider only one of the two forms in the analyses in Section 5.The KG form has been chosen therefore as the representative split form in this study, since it is known to be more robust for properly compressible flows [32].

Fig. 1 .
Fig. 1.Comparison of integration errors for the 1D Burgers' equation with turbulent-like solution.The transformed (modal space) right-hand side of (3.2) is shown vs. the polynomial mode i for different polynomial orders N .Conservative, quadratic split form and non-conservative DGSEM forms are compared against reference values obtained with exact integration.

Fig. 5 .
Fig. 5. Comparing QR diagrams for (a) very high-order KG solution with Roe-type numerical flux and (b) Energy-preserving KG solution without any additional interface stabilisation in the numerical flux, see Eq. (2.32), which is stable for m ≤ 4.

Fig. 6 .
Fig.6.Comparison between energy spectra of consistent integration and split forms for case m = 3, n el = 37.
Evolution of the enstrophy.
Energy spectrum at t = 9.QR diagram for KG split form at t = 9.Pointwise difference between KG and DU.

Fig. A. 8 .
Fig. A.8.Comparison of split forms DU and KG for the test case m = 5, n el = 23 with Roe-type numerical flux.

Table 1
Summary of cases -crossed out cases crashed with consistent/over-integration, whereas all test cases ran to the final time with the KG and DU split forms.

Table 2
Time of crash (tc) vs. quadrature points Q for test case m = 8, n el = 14