Accuracy, Stability, and Performance Comparison between the Spectral Difference and Flux Reconstruction Schemes
Christopher Cox, Will Trojak, Tarik Dzanic, Freddie Witherden, Antony Jameson
AAccuracy, Stability, and Performance Comparison between the SpectralDi ff erence and Flux Reconstruction Schemes C. Cox a, ∗ , W. Trojak b , T. Dzanic b , F. D. Witherden b , A. Jameson a,b a Department of Aerospace Engineering, Texas A & M University, College Station, TX 77843, USA b Department of Ocean Engineering, Texas A & M University, College Station, TX 77843, USA
Abstract
We report the development of a discontinuous spectral element flow solver that includes theimplementation of both spectral di ff erence and flux reconstruction formulations. With this high orderframework, we have constructed a foundation upon which to provide a fair and accurate assessmentof these two schemes in terms of accuracy, stability, and performance with special attention to thetrue spectral di ff erence scheme and the modified spectral di ff erence scheme recovered via the fluxreconstruction formulation. Building on previous analysis of the spectral di ff erence and flux reconstructionschemes, we provide a novel nonlinear stability analysis of the spectral di ff erence scheme. Throughvarious numerical experiments, we demonstrate the additional stability a ff orded by the true, baselinespectral di ff erence scheme without explicit filtering or de-aliasing due to its inherent feature of staggeredflux points. This arrangement leads to favorable suppression of aliasing errors and improves stabilityneeded for under-resolved simulations of turbulent flows. Keywords: discontinuous spectral element, spectral di ff erence, flux reconstruction, implicit large eddysimulation
1. Introduction
Computational fluid dynamics presents practitioners with many challenges, chief among which isresolving the often wide range of length scales while keeping computational cost su ffi ciently low. This ∗ Corresponding author: C. Cox
Email addresses: [email protected] (C. Cox), [email protected] (W. Trojak), [email protected] (T. Dzanic), [email protected] (F. D. Witherden), [email protected] (A. Jameson)
Preprint submitted to Computers and Fluids a r X i v : . [ phy s i c s . c o m p - ph ] S e p is crucial if such simulations are to meaningfully impact engineering design cycles. Reynolds-AveragedNavier-Stokes (RANS) methods, the prevailing mode of choice in the industry, have exhibited significantshortcomings in simulating complex turbulent flows, and as such, there is considerable interest in thedevelopment of high-fidelity scale-resolving simulations. Although far superior in terms of accuracy,these scale-resolving simulations can be orders of magnitude more computationally expensive than theirRANS counterparts which makes them intractable for many practical engineering purposes. To addressthis challenge, various families of methods have emerged over several decades, one of which is thespectral element method (SEM), a set of high-order techniques that has been successfully used formany applications. These methods developed out of discontinuous techniques, such as that of Reed andHill [1], which forwent some solution continuity in favor of localizing the calculation to sub-domains.This sub-domain structure—with reduced inter-element communication—can increase the computationale ffi ciency through structured compute regions that are well suited to modern massively parallel computerarchitectures such as graphic processing units (GPU).Discontinuous SEM o ff ers geometric flexibility and reduced dissipation / dispersion errors for high-fidelity computations; however, application of these schemes to turbulent flow problems can beproblematic due to numerical instability issues. As the cost of resolving the finest physical length scalesgrows prohibitively large with increasing Reynolds number, scale-resolving simulations are typicallyrestricted to resolving only the statistically significant length scales. For a su ffi ciently high-order scheme,this lack of resolution can cause aliasing errors to occur and produce unstable simulations [2]. Theseerrors originate from the high-order of the flux function and / or the geometry and limit the space inwhich the approximate solution can reside [3]. To ameliorate these errors and achieve stability, varioustechniques have been introduced, such as spectral vanishing viscosity methods (SVV) [4, 5, 6, 7], modalfiltering [8, 9, 10, 11], and split skew-symmetric methods [12, 13, 14, 15]. However, these techniquesdo come with a notable computational cost and, in some cases, tunable parameters, and it has becomecommonplace to perform simulations without explicit filtering or de-aliasing applied to the solution. Onesuch approach in the context of solving turbulent flows with discontinuous SEM is implicit large eddysimulation (ILES) [16, 17, 18, 19, 20], from which high-fidelity solutions can be obtained without anyadded modeling or filtering traditionally used to account for sub-grid length scales by utilizing the inherentnumerical dissipation of the scheme. However, this dissipation may be insu ffi cient when using high-orderdiscretizations for high Reynolds number turbulent flows, and it is not yet evident which method is bestsuited for robustly achieving stable and accurate simulations for these flows. There is speculation thatcertain methods may have more favorable de-aliasing properties which can result in improvements instability, although it has not been thoroughly explored.In this paper, we investigate two nodal discontinuous spectral element methods with several similarities.The first method is the flux reconstruction (FR) method of Huynh [21] and Vincent et al. [22]. Thismethod uses a local polynomial approximation of the solution to form an approximation to the fluxsuch that continuity is enforced through inter-element communication and correction functions. Thismethod has been adapted for several element topologies [23] and has been applied to various equationsets including the Euler equations [24, 25], Navier–Stokes equations [24, 26], and their incompressiblecounterparts [27, 28]. Several implementations of FR are available that have demonstrated the possibilityto achieve high computational e ffi ciency and scalability on large problems [29, 30]. The second methodis the spectral di ff erence (SD) method originally put forth by Kopriva et al. [31, 32], where a staggeredarrangement of points is used within each element, with one set of points for the solution and anotherfor the flux and its gradient. The formal stability of this method for linear problems was explored byJameson [33], who found a Lobatto-type distribution for the flux points to be important. Furthermore,Huynh [21] found that the accuracy of the scheme is independent of the solution point locations for linearproblems. Similar to FR, this method has been successfully applied to non-linear equations [34, 35, 36] aswell as in the simulation of complex physics [37, 38, 39, 40].The SD method is of interest as the approximation of the flux function, which is projected into thesolution space through di ff erentiation, is one degree higher than the solution. It is conjectured that thisincreased order of the flux equips SD with a favorable amount of de-aliasing in comparison to FR. Inthe body of SD and FR literature, there has been little comparative study between these related methodsand the e ff ect that di ff erent techniques for the flux function approximation will have on the stability andaccuracy of the methods. We investigate the di ff erences and similarities for these schemes when used inILES, and show the e ff ects of the higher degree of the flux approximation on the stability of the method.To this end, this work is structured with the formulation of SD and FR schemes on hyper-cube elements inSection 2. Non-linear analysis of the SD method is presented in Section 3, where the instability mechanicsare considered as well as scaling arguments for the error. Section 4 sets forth the formulation used forthe Navier–Stokes equations and Section 5 details results from numerical experiments for a series of testcases. Finally, conclusions are drawn in Section 6.
2. Discontinuous Spectral Element Formulations on Hexahedral Elements
For the sake of completeness, we briefly describe in the following sections the SD and FR schemes ontensor product hexahedral elements such that a self-contained comparison of the di ff erent formulationscan be made. We will begin by prescribing the shared definitions for partitioning the domain, reference domain,and how transformation from the reference domain and physical domain are constructed. The arbitraryconnected solution domain Ω ⊂ R is partitioned into N e non-overlapping, conforming, hexahedralelements, each denoted by Ω e , such that Ω = N e (cid:91) e = Ω e , N e (cid:92) e = Ω e = ∅ . (1)Each three-dimensional physical element Ω e is mapped to a reference element Ω r = { ξ, η, β | − (cid:54) ξ, η, β (cid:54) } through a mapping of the form x ( ξ, η, β ) = K (cid:88) k = x k φ k ( ξ, η, β ) , (2)where K is the number of nodes per element Ω e , x k = ( x k , y k , z k ) are nodal Cartesian coordinates, and φ k ( ξ, η, β ) are the nodal shape functions. After transformation into the computational domain, thegoverning equations in Eq. (42) can be re-written in the form ∂ ˆ U ∂ t + ∂ ˆ f ∂ξ + ∂ ˆ g ∂η + ∂ ˆ h ∂β = U = | J | U , ˆ f ˆ g ˆ h = | J | J − fgh . (4)For stationary grids, the Jacobian is defined as J = ∂ ( x , y , z ) /∂ ( ξ, η, β ). This information is needed at boththe solution and flux points within each reference element in accordance with the spectral di ff erence andflux reconstruction methodologies described in Sections 2.2 and 2.3. ff erence Following the original work of Kopriva and Kolias [31, 32] and Lui et al. [41], we briefly describe herethe three-dimensional spectral di ff erence formulation for which the distribution of solution points in areference cube can be interpreted from the distribution of points in the reference square shown in Fig. 1a.In this two-dimensional representation, the number of solution points (blue circles) along each direction isfour—these points, representing a polynomial of order (cid:112) =
3, are located at Gauss–Legendre quadraturepoints. The number of flux points (black squares) along each direction is one higher than the number ofsolution points—these points are also located at Gauss–Legendre quadrature points in the interior plus thetwo end points at -1 and 1. Using the (cid:112) + (cid:112) + (cid:112) and (cid:112) + ξ direction can be built using (cid:108) i ( ξ ) = (cid:112) + (cid:89) s = s (cid:44) i (cid:18) ξ − ξ s ξ i − ξ s (cid:19) ∀ i ∈ { , . . . , (cid:112) + } , (5a) (cid:104) i + ( ξ ) = (cid:112) + (cid:89) s = s (cid:44) i (cid:18) ξ − ξ s + ξ i + − ξ s + (cid:19) ∀ i ∈ { , . . . , (cid:112) + } , (5b)with analogous definitions made for the η and β directions. Here it can be observed that (cid:108) i ( ξ s ) = δ is , andthe complete polynomial approximation can be obtained within Ω r through tensor products of the three (cid:112) degree one-dimensional Lagrange polynomials by U δ r ( ξ, η, β ) = (cid:112) + (cid:88) k = (cid:112) + (cid:88) j = (cid:112) + (cid:88) i = ˆ U δ r | i , j , k | J r | i , j , k | (cid:108) i ( ξ ) (cid:108) j ( η ) (cid:108) k ( β ) (6)where ˆ U δ r | i , j , k = ˆ U δ r ( ξ i , η j , β k ) are the nodal coe ffi cients of the solution in Ω r that represent the value ofthe approximate solution polynomial ˆ U δ r evaluated at the set of solution points. The values of the fluxvectors can be obtained in a similar manner, but instead by using the three (cid:112) + (cid:104) i + , (cid:104) j + and (cid:104) k + byˆ f δ r ( ξ, η, β ) = (cid:112) + (cid:88) k = (cid:112) + (cid:88) j = (cid:112) + (cid:88) i = ˆ f δ r | i + , j , k (cid:104) i + ( ξ ) (cid:108) j ( η ) (cid:108) k ( β ) , (7a)ˆ g δ r ( ξ, η, β ) = (cid:112) + (cid:88) k = (cid:112) + (cid:88) j = (cid:112) + (cid:88) i = ˆ g δ r | i , j + , k (cid:108) i ( ξ ) (cid:104) j + ( η ) (cid:108) k ( β ) , (7b)ˆ h δ r ( ξ, η, β ) = (cid:112) + (cid:88) k = (cid:112) + (cid:88) j = (cid:112) + (cid:88) i = ˆ h δ r | i , j , k + (cid:108) i ( ξ ) (cid:108) j ( η ) (cid:104) k + ( β ) . (7c)The nodal coe ffi cients, ˆ f δ r | i + / , j , k = ˆ f δ r ( ξ i + / , η j , β k ), of the approximate discontinuous fluxes, ˆ f δ r , arecomputed from the solution at the flux points ˆ U r | i + / , j , k obtained by Eq. (6). Similar expressions canbe defined for ˆ g δ and ˆ h δ . The gradients at the solution points are computed using the solution at theflux points with the derivative of the Lagrange polynomial approach (see Sun et al. [35]). The gradientscan then be interpolated from the solution points to the flux points using a similar Lagrange interpolationapproach given in Eq. (6) to obtain the terms ˆ ∇ ˆ U δ r | i + / , j , k , ˆ ∇ ˆ U δ r | i , j + / , k , and ˆ ∇ ˆ U δ r | i , j , k + / . These gradientsare needed only for evaluation of the viscous fluxes.The common inviscid flux ˆ f δ Ie ( U δ L , U δ R ) at an interface between elements in the reference space canbe computed using any suitable approximate or exact Riemann solver, where the subscripts L and R denote left and right states of an interface . Similar expressions can be defined for ˆ g δ Ie and ˆ h δ Ie . In theSD implementation, the common viscous fluxes such as ˆ f δ Iv ( U δ L , ∇ U δ L , U δ R , ∇ U δ R ) are computed using anapproach analogous to inviscid Riemann solvers. In this work, we use the simple averaging approachfrom Bassi and Rebay (BR1) [42]. Note that the fluxes in Eqs. (7a)-(7c) are continuous within eachelement, but discontinuous across element interfaces. Globally continuous fluxes can be achieved in SDby replacing the interpolated values of the fluxes at element interfaces (denoted by a 1 / (cid:112) + / (a) (b)Figure 1: Distribution of solution points (SP • ) and flux points (FP (cid:4) ) with (cid:112) = ff erence method and (b)flux reconstruction method inside a unit reference element Ω r . with the common fluxes such that derivatives of the continuous fluxes can then be written as ∂ ˆ f δ Cr ∂ξ = (cid:112) + (cid:88) k = (cid:112) + (cid:88) j = (cid:34) ˆ f δ Ir | , j , k d (cid:104) ( ξ )d ξ + ˆ f δ Ir | (cid:112) + , j , k d (cid:104) (cid:112) + ( ξ )d ξ + (cid:112) (cid:88) i = ˆ f δ r | i + , j , k d (cid:104) i + ( ξ )d ξ (cid:35) (cid:108) j ( η ) (cid:108) k ( β ) , (8a) ∂ ˆ g δ Cr ∂η = (cid:112) + (cid:88) k = (cid:112) + (cid:88) i = (cid:34) ˆ g δ Ir | i , , k d (cid:104) ( η )d η + ˆ g δ Ir | i , (cid:112) + , k d (cid:104) (cid:112) + ( η )d η + (cid:112) (cid:88) j = ˆ g δ r | i , j + , k d (cid:104) j + ( η )d η (cid:35) (cid:108) i ( ξ ) (cid:108) k ( β ) , (8b) ∂ ˆ h δ Cr ∂β = (cid:112) + (cid:88) j = (cid:112) + (cid:88) i = (cid:34) ˆ h δ Ir | i , j , d (cid:104) ( β )d β + ˆ h δ Ir | i , j , (cid:112) + d (cid:104) (cid:112) + ( β )d β + (cid:112) (cid:88) k = ˆ h δ r | i , j , k + d (cid:104) k + ( β )d β (cid:35) (cid:108) i ( ξ ) (cid:108) j ( η ) . (8c) Following the original work by Huynh [21, 43], we briefly describe here the three-dimensionalflux reconstruction formulation for which the distribution of solution points in a reference cube canbe interpreted from the distribution of points in the reference square shown in Fig. 1b. In this 2Drepresentation, the number of solution points (blue circles) along each direction is four—these points,representing a polynomial of order (cid:112) =
3, are located at Gauss–Legendre quadrature points. The fluxpoints (black squares) along each direction are located at the two end points at -1 and 1. Using the solutionat the (cid:112) + (cid:112) degree Lagrange interpolating polynomial along each ξ , η , and β directioncan be constructed using Eq. (5a). Tensor products may once again be applied on the one dimensionalLagrange polynomial to obtain a complete approximation of the solution and the fluxes by U δ r ( ξ, η, β ) = (cid:112) + (cid:88) k = (cid:112) + (cid:88) j = (cid:112) + (cid:88) i = ˆ U δ r | i , j , k | J r | i , j , k | (cid:108) i ( ξ ) (cid:108) j ( η ) (cid:108) k ( β ) , (9a)ˆ f δ r ( ξ, η, β ) = (cid:112) + (cid:88) k = (cid:112) + (cid:88) j = (cid:112) + (cid:88) i = ˆ f δ r | i , j , k (cid:108) i ( ξ ) (cid:108) j ( η ) (cid:108) k ( β ) , (9b)ˆ g δ r ( ξ, η, β ) = (cid:112) + (cid:88) k = (cid:112) + (cid:88) j = (cid:112) + (cid:88) i = ˆ g δ r | i , j , k (cid:108) i ( ξ ) (cid:108) j ( η ) (cid:108) k ( β ) , (9c)ˆ h δ r ( ξ, η, β ) = (cid:112) + (cid:88) k = (cid:112) + (cid:88) j = (cid:112) + (cid:88) i = ˆ h δ r | i , j , k (cid:108) i ( ξ ) (cid:108) j ( η ) (cid:108) k ( β ) . (9d)In FR, as in SD, the nodal coe ffi cients, ˆ f δ r | i , j , k = ˆ f δ r ( ξ i , η j , β k ), of the approximate discontinuous fluxes ˆ f δ r are computed from ˆ U δ r | i , j , k and ˆ ∇ ˆ U δ r | i , j , k , where the latter term is only required for the viscous fluxes. Similarexpressions can be defined for ˆ g δ and ˆ h δ . In accordance with the methodology of the flux reconstructionscheme, the continuous flux functions defined along ξ , η and β directions can be written asˆ f δ Cr = ˆ f δ r + (cid:2) ˆ f δ Ir − − ˆ f δ r ( − , η, β ) (cid:3) g LB ( ξ ) + (cid:2) ˆ f δ Ir + − ˆ f δ r (1 , η, β ) (cid:3) g RB ( ξ ) , (10a)ˆ g δ Cr = ˆ g δ r + (cid:2) ˆ g δ Ir − − ˆ g δ r ( ξ, − , β ) (cid:3) g LB ( η ) + (cid:2) ˆ g δ Ir + − ˆ g δ r ( ξ, , β ) (cid:3) g RB ( η ) , (10b)ˆ h δ Cr = ˆ h δ r + (cid:2) ˆ h δ Ir − − ˆ h δ r ( ξ, η, − (cid:3) g LB ( β ) + (cid:2) ˆ h δ Ir + − ˆ h δ r ( ξ, η, (cid:3) g RB ( β ) (10c)where g LB and g RB represent left boundary (LB) and right boundary (RB) correction functions inthe reference element, respectively. A stable correction function as defined by Huynh [21] andVincent et al. [22] can be generalised for the left boundary as g LB ( ξ ) = α (cid:82) R , (cid:112) + ( ξ ) + (1 − α ) (cid:82) R , (cid:112) ( ξ ) (11)where (cid:82) R , ( · ) ( ξ ) represents the right Radau polynomial [44]. The expression for a correction to theright boundary is obtained simply by reflection of g LB ( ξ ) such that g RB ( ξ ) = g LB ( − ξ ) on the interval Ω r = { ξ | − (cid:54) ξ (cid:54) } . Choosing α = α = ( (cid:112) + / (2 (cid:112) +
1) recovers a modified SD method—in the current work, it is this scheme to which we directly compare true SD. Another type of schemecan be obtained by setting α = (cid:112) / (2 (cid:112) + g scheme identified byHuynh [21] that collocates solution points with the Lobatto points. These three schemes are referred toherein as FR DG , FR SD , and FR , respectively. Lastly, Romero et al. [45] provided a simplified formulationof the FR scheme that substitutes a Lagrange interpolation operation for the correction functions. Theyo ff ered a proof of equivalence of their scheme to FR DG , provided that solution points are placed at thecorresponding Gauss–Legendre points. This method is referred to as direct FR (DFR) [45, 46].From Eqs. (10a)-(10c), we can obtain the derivatives of the continuous flux functions ∂ ˆ f δ Cr ∂ξ = (cid:112) + (cid:88) k = (cid:112) + (cid:88) j = (cid:112) + (cid:88) i = ˆ f δ r | i , j , k d (cid:108) i ( ξ )d ξ (cid:108) j ( η ) (cid:108) k ( β ) + (cid:2) ˆ f δ Ir − − ˆ f δ r ( − , η, β ) (cid:3) d g LB ( ξ )d ξ + (cid:2) ˆ f δ Ir + − ˆ f δ r (1 , η, β ) (cid:3) d g RB ( ξ )d ξ , (12a) ∂ ˆ g δ Cr ∂η = (cid:112) + (cid:88) k = (cid:112) + (cid:88) j = (cid:112) + (cid:88) i = ˆ g δ r | i , j , k (cid:108) i ( ξ ) d (cid:108) j ( η )d η (cid:108) k ( β ) + (cid:2) ˆ g δ Ir − − ˆ g δ r ( ξ, − , β ) (cid:3) d g LB ( η )d η + (cid:2) ˆ g δ Ir + − ˆ g δ r ( ξ, , β ) (cid:3) d g RB ( η )d η , (12b) ∂ ˆ h δ Cr ∂β = (cid:112) + (cid:88) k = (cid:112) + (cid:88) j = (cid:112) + (cid:88) i = ˆ h δ r | i , j , k (cid:108) i ( ξ ) (cid:108) j ( η ) d (cid:108) k ( β )d β + (cid:2) ˆ h δ Ir − − ˆ h δ r ( ξ, η, − (cid:3) d g LB ( β )d β + (cid:2) ˆ h δ Ir + − ˆ h δ r ( ξ, η, (cid:3) d g RB ( β )d β . (12c)In the FR implementation, the common viscous fluxes are computed using a BR2-type, second procedureof Bassi and Rebay [47] written using the flux reconstruction methodology [43] to achieve compactnessof the stencil in multiple dimensions.Once the divergence of the continuous flux is obtained by Eqs. (8a)-(8c) for the spectral di ff erencescheme or Eqs. (12a)-(12c) for the flux reconstruction scheme, an appropriate time stepping technique canbe applied to march the solution forward in time. The implementation of both schemes is done within asingle coding framework such that fair and proper comparisons of the two methodologies can be made interms of stability, accuracy, and performance.0
3. Nonlinear Stability of Spectral Di ff erence In the work of Jameson et al. [48], the non-linear stability of the flux reconstruction method wasinvestigated, and it was found that the solution decay could be decomposed into a stable component anda non-linear component which can cause instabilities. As this analysis was useful in understanding themechanism by which non-linearities a ff ect stability and how de-aliasing methods can mitigate this, wewill perform a similar analysis for the spectral di ff erence method in order to highlight the di ff erences thatarise between these two techniques.Consider a scalar conservation law in one dimension ∂ u ∂ t + ∂ f ∂ x = , (13)where the lower-case terms denote a scalar quantity. This may be cast in the reference domain as ∂ ˆ u δ ∂ t + ∂ ˆ f δ ∂ξ = . (14)As was introduced in the previous section, the approximate flux ˆ f δ in the FR and SD methodologies isreplaced by a corrected flux ˆ f δ C that enforces C continuity in the flux between elements. A similarexpression for the corrected flux used in the flux reconstruction method can be written for the spectraldi ff erence method as ˆ f δ C = ˆ f δ + (cid:16) ˆ f δ IL − ˆ f δ L (cid:17) (cid:104) + (cid:16) ˆ f δ IR − ˆ f δ R (cid:17) (cid:104) (cid:112) + . For brevity, we temporarily drop the subscript r and refer to left and right interfaces of a given elementwith the subscripts L and R . Therefore, ∂ ˆ u δ ∂ t = − ∂ ˆ f δ ∂ξ − (cid:16) ˆ f δ IL − ˆ f δ L (cid:17) d (cid:104) d ξ − (cid:16) ˆ f δ IR − ˆ f δ R (cid:17) d (cid:104) (cid:112) + d ξ . (15)To analyze stability in a norm that induces a Sobolev space, namely (cid:107) ˆ v (cid:107) W , (cid:112) ,ι = (cid:90) − ˆ v + ι (cid:0) ∂ (cid:112) ˆ v (cid:1) d ξ, (16)1we investigate the behavior of dd t (cid:107) ˆ v (cid:107) W , (cid:112) ,ι = dd t (cid:90) − ˆ v + ι (cid:0) ∂ (cid:112) ˆ v (cid:1) d ξ. (17)By taking Eq. (15) and, following the work of Jameson et al. [48], multiplying it by ˆ u δ and integrating,we obtain12 dd t (cid:90) − ( ˆ u δ ) d ξ = − (cid:90) − ˆ u δ ∂ ˆ f δ ∂ξ d ξ − (cid:16) ˆ f δ IL − ˆ f δ L (cid:17) (cid:90) − ˆ u δ d (cid:104) d ξ d ξ − (cid:16) ˆ f δ IR − ˆ f δ R (cid:17) (cid:90) − ˆ u δ d (cid:104) (cid:112) + d ξ d ξ. (18)Upon using the product rule, this may be rewritten as12 dd t (cid:90) − ( ˆ u δ ) d ξ = (cid:90) − ˆ f δ ∂ ˆ u δ ∂ξ d ξ + (cid:16) ˆ f δ IL − ˆ f δ L (cid:17) (cid:90) − (cid:104) d ˆ u δ d ξ d ξ + (cid:16) ˆ f δ IR − ˆ f δ R (cid:17) (cid:90) − (cid:104) (cid:112) + d ˆ u δ d ξ d ξ + (cid:16) ˆ f δ IL ˆ u δ L − ˆ f δ IR ˆ u δ R (cid:17) . (19)Furthermore, taking Eq. (15) and di ff erentiating it (cid:112) times gives ∂∂ t (cid:18) ∂ (cid:112) ˆ u δ ∂ξ (cid:112) (cid:19) = − ∂ (cid:112) + ˆ f δ ∂ξ (cid:112) + − (cid:16) ˆ f δ IL − ˆ f δ L (cid:17) d (cid:112) + (cid:104) d ξ (cid:112) + − (cid:16) ˆ f δ IR − ˆ f δ R (cid:17) d (cid:112) + (cid:104) (cid:112) + d ξ (cid:112) + . (20)A key di ff erence between this derivation and that for the flux reconstruction method is that ˆ f δ is apolynomial of degree (cid:112) +
1, and therefore the first term on the right-hand side is not zero, but a constant.If Eq. (20) is then multiplied by the (cid:112) th derivative of ˆ u δ and integrated, the following is obtained12 dd t (cid:90) − (cid:18) ∂ (cid:112) ˆ u δ ∂ξ (cid:112) (cid:19) d ξ = − (cid:90) − ∂ (cid:112) + ˆ f δ ∂ξ (cid:112) + ∂ (cid:112) ˆ u δ ∂ξ (cid:112) d ξ − (cid:16) ˆ f δ IL − ˆ f δ L (cid:17) (cid:90) − ∂ (cid:112) ˆ u δ ∂ξ (cid:112) d (cid:112) + (cid:104) d ξ (cid:112) + d ξ − (cid:16) ˆ f δ IR − ˆ f δ R (cid:17) (cid:90) − ∂ (cid:112) ˆ u δ ∂ξ (cid:112) d (cid:112) + (cid:104) (cid:112) + d ξ (cid:112) + d ξ, (21)which, in turn, may be written as12 dd t (cid:90) − (cid:18) ∂ (cid:112) ˆ u δ ∂ξ (cid:112) (cid:19) d ξ = − ∂ (cid:112) + ˆ f δ ∂ξ (cid:112) + ∂ (cid:112) ˆ u δ ∂ξ (cid:112) − (cid:16) ˆ f δ IL − ˆ f δ L (cid:17) ∂ (cid:112) ˆ u δ ∂ξ (cid:112) d (cid:112) + (cid:104) d ξ (cid:112) + − (cid:16) ˆ f δ IR − ˆ f δ R (cid:17) ∂ (cid:112) ˆ u δ ∂ξ (cid:112) d (cid:112) + (cid:104) (cid:112) + d ξ (cid:112) + . (22)2By combining Eqs. (20) and (22) and taking the norm as given by Eq. (16), we then obtain12 dd t (cid:107) ˆ u δ (cid:107) W , (cid:112) ,ι = (cid:90) − ˆ f δ ∂ ˆ u δ ∂ξ d ξ + (cid:16) ˆ f δ IL − ˆ f δ L (cid:17) (cid:90) − (cid:104) d ˆ u δ d ξ d ξ + (cid:16) ˆ f δ IR − ˆ f δ R (cid:17) (cid:90) − (cid:104) (cid:112) + d ˆ u δ d ξ d ξ − ι ∂ (cid:112) + ˆ f δ ∂ξ (cid:112) + ∂ (cid:112) ˆ u δ ∂ξ (cid:112) − ι (cid:16) ˆ f δ IL − ˆ f δ L (cid:17) ∂ (cid:112) ˆ u δ ∂ξ (cid:112) d (cid:112) + (cid:104) d ξ (cid:112) + − ι (cid:16) ˆ f δ IR − ˆ f δ R (cid:17) ∂ (cid:112) ˆ u δ ∂ξ (cid:112) d (cid:112) + (cid:104) (cid:112) + d ξ (cid:112) + + (cid:16) ˆ f δ IL ˆ u δ L − ˆ f δ IR ˆ u δ R (cid:17) . (23)To proceed further, we refer to the work of Huynh [21] who showed that for the linear case, SD couldbe recovered from FR for a given correction function. This correction function, which we denote as g ,recovered SD when the SD flux points were formed from the (cid:112) degree Gauss–Legendre quadrature pointswith points added at − (cid:104) = g L and (cid:104) (cid:112) + = g R . (24)Therefore, we may write Eq. (23) as12 dd t (cid:107) ˆ u δ (cid:107) W , (cid:112) ,ι = (cid:90) − ˆ f δ ∂ ˆ u δ ∂ξ d ξ + (cid:16) ˆ f δ IL − ˆ f δ L (cid:17) (cid:90) − g L d ˆ u δ d ξ d ξ + (cid:16) ˆ f δ IR − ˆ f δ R (cid:17) (cid:90) − g R d ˆ u δ d ξ d ξ − ι ∂ (cid:112) + ˆ f δ ∂ξ (cid:112) + ∂ (cid:112) ˆ u δ ∂ξ (cid:112) − ι (cid:16) ˆ f δ IL − ˆ f δ L (cid:17) ∂ (cid:112) ˆ u δ ∂ξ (cid:112) d (cid:112) + g L d ξ (cid:112) + − ι (cid:16) ˆ f δ IR − ˆ f δ R (cid:17) ∂ (cid:112) ˆ u δ ∂ξ (cid:112) d (cid:112) + g R d ξ (cid:112) + + ( ˆ f δ IL ˆ u δ L − ˆ f δ IR ˆ u δ R ) , (25)and from Vincent et al. [22], we use (cid:90) − g L ∂ ˆ u δ ∂ξ d ξ − ι ∂ (cid:112) ˆ u δ ∂ξ (cid:112) d (cid:112) + g L d ξ (cid:112) + = , (26a) (cid:90) − g R ∂ ˆ u δ ∂ξ d ξ − ι ∂ (cid:112) ˆ u δ ∂ξ (cid:112) d (cid:112) + g R d ξ (cid:112) + = , (26b)which reduces Eq. (25) to12 dd t (cid:107) ˆ u δ (cid:107) W , (cid:112) ,ι = (cid:90) − ˆ f δ ∂ ˆ u δ ∂ξ d ξ + (cid:16) ˆ f δ IL ˆ u δ L − ˆ f δ IR ˆ u δ R (cid:17) − ι ∂ (cid:112) + ˆ f δ ∂ξ (cid:112) + ∂ (cid:112) ˆ u δ ∂ξ (cid:112) . (27)3If the broken norm is then constructed from this for N elements on a periodic domain, we obtain12 dd t (cid:107) u δ (cid:107) W , (cid:112) ,ι = Θ + N − (cid:88) i = (cid:15) i − ι N − (cid:88) i = ∂ (cid:112) + f δ ∂ x (cid:112) + ∂ (cid:112) u δ ∂ x (cid:112) , (28)where (cid:15) i = (cid:90) Ω i ( f δ − f ) ∂ u δ ∂ x d x . (29)Here, the term Θ is the interface contribution to the stability for which a full derivation can be found in[48], and the reader is referred to that work for a more complete derivation. If the common interface valuesare set such that they form an E-flux [49, 50], then Θ (cid:54) ff ect for SD.To illustrate more clearly the e ff ect that the di ff erence between the schemes has on the approximationof the flux gradient, we will now examine the error scaling. Using theorems and corollaries presented byBernardi and Maday [51], we further analyze the behavior of the error in the flux evaluated in the L norm (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∂ f ∂ x − ∂ f δ C ∂ x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L . Throughout, we adopt a similar notation to Bernardi and Maday [51], where we define the Sobolev space H k ( X ) = (cid:110) v ∈ L ( X ) | ∀ m ∈ N , m (cid:54) k , ∂ m v ∈ L ( X ) (cid:111) , where X is an open, bounded, Lipschitz-continuous set of R , and the norm induced on the space H k is (cid:107) u (cid:107) H k = (cid:118)(cid:117)(cid:116)(cid:90) X k (cid:88) m = (cid:0) ∂ m u (cid:1) d x . (30)To analyze the behavior of the flux error, we establish two necessary theorems. Theorem 3.1. (See Bernardi and Maday [51], Thm. 13.2.) For some function u ∈ H k with k > / and theLagrange interpolation operator I (cid:112) g ∈ P (cid:112) such that I (cid:112) g ( ζ i ) = g ( ζ i ) for some points set of points { ζ i } i (cid:54) (cid:112) + ,the following estimate holds (cid:107) u − I (cid:112) u (cid:107) L (cid:54) C ( k )( (cid:112) + − k (cid:107) u (cid:107) H k , (31)4 for some constant C that is only dependent on k. Theorem 3.2. (See Bernardi and Maday [51], Thm. 13.4.) For some function u ∈ H k with real numbersk and r such that k (cid:62) and r < k and the Lagrange interpolation operator I (cid:112) g defined in Thm. 3.1, thefollowing estimates hold (cid:107) u − I (cid:112) u (cid:107) H r (cid:54) C ( (cid:112) + r / − k (cid:107) u (cid:107) H k if r (cid:54) , C ( (cid:112) + r − / − k (cid:107) u (cid:107) H k if r (cid:62) . (32)With these theorems established, we look to determine the bound on (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∂ f ∂ x − ∂ f δ C ∂ x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L = (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∂ f ∂ x − (cid:20) ∂ f δ ∂ x + (cid:16) f δ IL − f δ L (cid:17) d g L d x + (cid:16) f δ IR − f δ R (cid:17) d g R d x (cid:21)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L . (33)From the triangle inequality, this may be rewritten as (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∂ f ∂ x − ∂ f δ C ∂ x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L (cid:54) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∂ f ∂ x − ∂ f δ ∂ x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L + (cid:12)(cid:12)(cid:12) f δ IL − f δ L (cid:12)(cid:12)(cid:12)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) d g L d x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L + (cid:12)(cid:12)(cid:12) f δ IR − f δ R (cid:12)(cid:12)(cid:12)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) d g R d x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L . (34)We impose that the interface values take the form f δ Ir | L = κ L f δ r − | R + (1 − κ ) f δ r | L , for κ L ∈ [0 , , (35)where r denotes the element index, and impose similar behavior at the opposite interface. We may thenwrite (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∂ f ∂ x − ∂ f δ C ∂ x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L (cid:54) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∂ f ∂ x − ∂ f δ ∂ x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L + κ L (cid:12)(cid:12)(cid:12) f δ r − | R − f δ r | L (cid:12)(cid:12)(cid:12)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) d g L d x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L + κ R (cid:12)(cid:12)(cid:12) f δ r | R − f δ r + | L (cid:12)(cid:12)(cid:12)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) d g R d x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L . (36)As the true interface term is the same for both sides of the interface, we may write κ L (cid:16) f δ r − | R − f δ r | L (cid:17) = κ L (cid:16) f δ r − | R − f r − | R + f r | L − f δ r | L (cid:17) , (37)which can be generalized for the other interface. Under the assumption that f is a high-order function of u such that if u ∈ H k then f ∈ H mk for m (cid:62)
1, the interface correction will scale with the interpolation errorof f κ L (cid:12)(cid:12)(cid:12) f δ r − | R − f δ r | L (cid:12)(cid:12)(cid:12) (cid:54) κ L C ( (cid:112) + − mk (cid:107) f (cid:107) H mk . (38)5It is then straightforward to prove the following bound for a Lagrange polynomial (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) d h i + d x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L (cid:54) C ( (cid:112) + , and as the correction function for SD is a Lagrange polynomial, we may use this to give κ L (cid:12)(cid:12)(cid:12) f δ r − | R − f δ r | L (cid:12)(cid:12)(cid:12)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) d g L d x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L (cid:54) κ L C ( (cid:112) + − mk (cid:107) f (cid:107) H mk . (39)Considering the first term on the right-hand-side of Eq. (36), we can modify Thm. 3.2 to yield (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∂ f ∂ x − ∂ f δ ∂ x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L (cid:54) (cid:13)(cid:13)(cid:13) f − f δ (cid:13)(cid:13)(cid:13) H (cid:54) C ( (cid:112) + / − mk (cid:107) f (cid:107) H mk . Combining these results and taking into account that κ ∈ [0 , (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∂ f ∂ x − ∂ f δ C ∂ x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L (cid:54) C ( (cid:112) + / − mk (cid:107) f (cid:107) H mk . (40)Repeating these steps for FR, we find the similar and expected relation that (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∂ f ∂ x − ∂ f δ C ∂ x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L (cid:54) C ( (cid:112) + / − mk (cid:107) f (cid:107) H mk . (41)As a result, the error of SD can be lower due to the di ff erent scaling, the di ff erence being most evidentwhen the ratio ( (cid:112) + / ( (cid:112) +
1) is largest and k is large compared to (cid:112) (i.e. in under-resolved cases). Weremark that this result is separate from arguments concerning the study of the scheme’s asymptotic rate ofconvergence with respect to grid spacing. In that case, it is known that DG-type FR schemes can obtainsuper-convergence one degree higher than SD and other FR variants [52, 53].
4. Governing Equations
Consider the full three-dimensional compressible Navier–Stokes equations written in strongconservation form for a Cartesian coordinate system ( x , y , z ) ∂ U ∂ t + ∂ f ∂ x + ∂ g ∂ y + ∂ h ∂ z = . (42)6The vector of state variables, U ( x , y , z , t ), is defined for [ x , y , z ] ∈ Ω ⊂ R and t ∈ R + , with U = [ ρ ρ u ρ v ρ w ρ E ] T and the flux vectors f , g and h contain both inviscid terms, denoted by ( · ) e , and viscousterms, denoted by ( · ) v , where f = f e − f v , g = g e − g v , h = h e − h v . (43)The inviscid flux vectors can be written as f e = ρ u ρ u + p ρ uv ρ uw ( ρ E + p ) u , g e = ρ v ρ vu ρ v + p ρ vw ( ρ E + p ) v , h e = ρ w ρ wu ρ wv ρ w + p ( ρ E + p ) w (44)and the viscous flux vectors can be written as f v = τ xx τ xy τ xz κ ∂ T ∂ x + u τ xx + v τ xy + w τ xz , g v = τ yx τ yy τ yz κ ∂ T ∂ y + u τ yx + v τ yy + w τ yz , h v = τ zx τ zy τ zz κ ∂ T ∂ z + u τ zx + v τ zy + w τ zz . (45)The total energy is E = p / [ ρ ( γ − + ( u + v + w ) / κ = ( µ c p ) / Pr . UnderStokes’ hypothesis, the bulk viscosity is assigned a value of zero, leading to the second coe ffi cient ofviscosity taking the value λ = − / µ ; therefore, we can write τ xx = µ ∂ u ∂ x + λ ∇ · u , τ yy = µ ∂ v ∂ y + λ ∇ · u , τ zz = µ ∂ w ∂ z + λ ∇ · u , (46a) τ xy = τ yx = µ (cid:32) ∂ v ∂ x + ∂ u ∂ y (cid:33) , τ xz = τ zx = µ (cid:32) ∂ w ∂ x + ∂ u ∂ z (cid:33) , τ yz = τ zy = µ (cid:32) ∂ w ∂ y + ∂ v ∂ z (cid:33) . (46b)In the formulation above, u , v , and w are the components of velocity in the x , y , and z directions,respectively, and ρ represents the density, p the pressure, µ the dynamic viscosity, ν the kinematic viscosity,7 Pr the Prandtl number, γ the specific heat ratio, and c p the specific heat at constant pressure. Unless statedotherwise, the Prandtl number and specific heat ratio are set constant at Pr = .
72 and γ = .
5. Numerical Experiments
The results from a series of numerical experiments performed comparing SD and FR SD will now bepresented. We will begin with a 1D linear test case that can be modified such at aliasing is introduced. Given alinear advection equation with variable propagation speed, an equivalent scalar conservation form can bederived ∂ u ∂ t + (2 − sin x ) ∂ u ∂ x = ⇔ ∂ u ∂ t + ∂ (2 − cos x ) u ∂ x = u sin x . (47)In the latter form, the equation introduces aliasing errors in numerical calculations, and thus is a suitablecandidate for identifying de-aliasing properties of numerical schemes without the presence of non-linearities. Furthermore, when this equation is applied to a periodic domain Ω = [0 , π ], the solutionis shown to analytically have a time period of T = π/ √
3, allowing for exact calculations of the error [54].The initial condition for this test was chosen to be a reconstruction of the energy spectra E ( k , t = = Ck k exp (cid:18) − k k (cid:19) , where C = √ π , and k = , (48)which is similar to the condition used by Alhawwary et al. [55] and San [56]. A 1D scalar field was thenreconstructed from the spectra as u ( x , t = = k max (cid:88) k = (cid:112) E ( k ,
0) cos ( kx + Ψ ( k )) , (49)where k max = Ψ ( k ) ∈ (0 , π ] is a random phase angle forwavenumber k . With this initial condition, multiple modes are excited while E ( k ) → k → ∞ , whichmakes di ff erences in aliasing evident.8 t − − − k E ( k ) t , (cid:112) = t , (cid:112) = SD , (cid:112) = SD , (cid:112) = , (cid:112) = , (cid:112) = (a) 120 DoF. − − − − − k E ( k ) t , (cid:112) = t , (cid:112) = SD , (cid:112) = SD , (cid:112) = , (cid:112) = , (cid:112) = (b) 600 DoF.Figure 2: Energy spectrum comparison of FR SD and SD with centrally-di ff erenced interfaces after one time period averagedover 1 × initial conditions. A comparison of the numerical results is shown in Fig. 2 for various polynomial orders and gridresolutions after one time period. For this case, only the average spectra results of the experiments usingcentrally-di ff erenced interfaces are shown as negligible di ff erences between SD and FR SD were observedwhen upwinding was used. This e ff ect can be attributed to the numerical dissipation of upwinded schemesat high frequencies which can be su ffi cient to dampen aliasing errors in this case. When using central-di ff erencing on the coarse grid, instabilities were evident in the spectra of the calculations using the FRmethod, whereas the SD method was stable. As the grid was refined, the FR method became stable buthad notably more energy at high wavenumbers than the SD method. As it can be shown analytically forthis equation that aliasing will be introduced at the highest wavenumbers and propagated to the lowerwavenumbers, it is evident that pure SD is more stable due to less aliasing error. To demonstrate and compare super-convergence of the flux reconstruction and spectral di ff erenceschemes [57, 58] for the Euler equations within the current implementation, we solve the isentropic Eulervortex [59] in a free-stream flow for which there exists an exact analytical solution. Super-convergencefor these types of schemes is said to be achieved once the observed order of accuracy is greater than (cid:112) + r c and strength (cid:15) , positioned in the domain at ( x o , y o ), and here we9will consider a vortex advecting purely in the y -direction. The analytical solution at ( x , y , t ) for this testcase is given by ρ ( x , y , t ) = ρ ∞ (cid:32) − ( γ − (cid:15) M ∞ π exp (2 f ) (cid:33) γ − , (50a) u ( x , y , t ) = U ∞ (cid:32) (cid:15) ( y − y o − U ∞ t )2 π r c exp ( f ) (cid:33) , (50b) v ( x , y , t ) = U ∞ (cid:32) − (cid:15) ( x − x o )2 π r c exp ( f ) (cid:33) , (50c) p ( x , y , t ) = p ∞ (cid:32) ρρ ∞ (cid:33) γ . (50d)where f = (1 − ( x − x o ) − ( y − y o − U ∞ t ) ) / r c . To match the conditions of Vincent et al. [52] andWitherden et al. [29], we set the free-stream conditions to ρ ∞ = U ∞ =
1, and p ∞ = ( ρ ∞ U ∞ ) / ( γ M ∞ ),where the free-stream Mach number is M ∞ = .
4. We prescribe the size and strength of the vortex to be r c = . (cid:15) = .
5, respectively, and initially position the vortex at the center of the domain located at( x o , y o ) = (20 , Ω = { x , y ∈ R | (cid:54) x , y (cid:54) } is partitioned using four di ff erent meshes of120 , 140 , 160 , and 180 elements. The upper and lower boundaries are treated as periodic while theleft and right boundaries are prescribed free-stream conditions. These conditions result in modeling aninfinite array of coupled vortices; however, the impact of the vortex on the free-stream at the boundariesis negligible since the vortex size r c is small compared to the length L =
40 of the domain and thevortex strength exponentially decays from its origin [52]. Therefore, we are e ff ectively modeling a vortexpropagating through an infinite domain. We consider a polynomial order (cid:112) =
3, which gives 480 , 560 ,640 , and 720 DoF for the various meshes. We use Davis’ form of the Rusanov approximate Riemannsolver [60] to compute inviscid numerical fluxes at the interfaces between elements, and we use the low-storage, five-stage, fourth-order accurate Runge–Kutta scheme of Carpenter and Kennedy [61] with a timestep of ∆ t = . × − to explicitly march the solution through time. This time step is small enoughsuch that all truncation errors are dominated by the spatial discretization.0 (a) (b)Figure 3: Isentropic Euler Vortex: (a) (cid:112) = ρ after 45 advective flow cycles ( t = Ω = { x , y ∈ R | (cid:54) x , y (cid:54) } which is partitioned into 120 elements, (b) enlarged image of the vortex centered at the origin( x , y ) = (20 ,
20) at t = To assess the order of accuracy, we compute the L -norm of the density error (cid:107) e (cid:107) inside an integrationwindow Ω I = { x , y ∈ R | − (cid:54) x − x o (cid:54) , − (cid:54) y − y o (cid:54) } at each moment in time the vortex advectsthrough the entire computational domain and returns to the origin, which occurs when t = t (cid:63) L / u ∞ for t (cid:63) ∈ { , , . . . , } . The L -norm of the density error is defined as (cid:107) e (cid:107) = (cid:115)(cid:90) Ω I ( ρ n ( x , y ) − ρ e ( x , y )) d x (51)where ρ n ( x , y ) is the numerical density and ρ e ( x , y ) is the exact analytical solution given in Eq. (50a) at t =
0. To approximate the integrals in Eq. (51), we apply a more than su ffi cient high-strength quadraturerule. To compare our results against those obtained by [52] and [29], we plot the observed convergenceof the FR SD and SD schemes in Fig. 4, where the order of accuracy at any given time is determined bycomputing the slope of the line given by a least-squares fit of log( (cid:107) e (cid:107) ) as a function of log( h ). For thefour di ff erent meshes, we use grid spacings h ∈ { / , / , / , / } . For comparison, we also plot resultsfrom other FR schemes built into the current solver in Fig. 4 including FR DG , FR and DFR. We observean approximate 2 (cid:112) + DG at t = (cid:112) under FR SD . We also confirmthat the super accuracy of the DFR scheme is equivalent to that of FR DG since solution points are placedat corresponding Gauss–Legendre points.1 Figure 4: Isentropic Euler Vortex: super accuracy with (cid:112) = DG ,FR SD , FR , DFR and SD. We can recast the nodal form of the solution polynomial into its modal form by using a set of modal basisfunctions—orthogonal Legendre polynomials (cid:76) i ( ξ ), (cid:76) j ( η )—and their corresponding modal coe ffi cients c i , j [3]. Following the work of Spiegel et al. [62], we plot | c i , j | within each element (see Fig. 5), normalizingby the mean mode and zeroing all modes less than 1 × − . In these images, the values of | c i , j | in thelower left corner of each element correspond to the magnitude of the mean mode c , (cid:76) ( ξ ) (cid:76) ( η ). Thevalues in the upper right corner of each element correspond to the magnitude of the highest Legendremode c (cid:112) , (cid:112) (cid:76) (cid:112) ( ξ ) (cid:76) (cid:112) ( η ). From left to right and bottom to top, these modal coe ffi cients correspond to themagnitude of the Legendre modes of increasing order with respect to ξ and η , respectively, up to (cid:112) .Under FR SD , we demonstrate in Fig. 5 that the higher frequency modes in regions away from the vortexare more energized in comparison to SD. The larger magnitudes of the higher modes in FR SD can beattributed to aliasing errors. By comparison, the SD scheme is successful at suppressing this energy at thehigher modes, with the dominant modes away from the vortex being the lowest order mean mode, which isconsistent with analytic solution. In turn this produces a lower error in the solution, as demonstrated by thetime history plot of the L -norm of density shown in Fig. 6a. As a result, this causes rate of convergencehistory to initially increase sharply to a level above 2 (cid:112) between t = t =
480 s, then level o ff for theremaining portion of the simulation. This rapid approach to an order greater than 2 (cid:112) indicates favorableaccuracy properties of the SD scheme, thereby reducing contamination of the solution from aliasing errors.This result is consistent with the analytical findings presented in Eqs. (40) and (41). Ultimately, this o ff ersimproved stability when performing implicit large eddy simulations of turbulent flow problems such as2those studied in Sec. 5.4 and Sec 5.5. (a) (b)Figure 5: Isentropic Euler Vortex: modal coe ffi cients for (cid:112) = { x , y | < x < , < y < } on a 120 × t =
40 s). (a) FR SD , (b) SD.(a) (b)Figure 6: Isentropic Euler Vortex: L -norm of density error || e || as a function of time for (cid:112) = black ) and FR SD ( red ),(b) L -norm of density error as a function of grid spacing h at t = In this section, we simulate the steady, two-dimensional, inviscid, subsonic flow over a cylinder asgoverned by the compressible Euler equations. This test case is constructed to assess numerically-3generated entropy and was used in Mengaldo et al. [63] to test the e ff ectiveness of global de-aliasingfor the FR DG scheme at di ff erent polynomial orders. Ideally, zero entropy should be generated for aninviscid, subsonic simulation, however aliasing in the numerical method introduces a mechanism allowingthe build-up of entropy. To reduce numerical entropy generation due to the mesh representation of thecylinder wall, the curvature of the cylinder is represented with 176 quartic elements with 54 elements inthe radial direction. The mesh, shown in Fig. 7a, extends 10 d into the farfield and contains a total of176 × = M ∞ = . (cid:112) = (cid:112) = (cid:112) = (a) (b) (c)Figure 7: Inviscid, subsonic flow over 2D cylinder: mach number (a) mesh; (b) FR SD , (cid:112) =
6; (c) SD, (cid:112) = Mach number contours from the (cid:112) = SD and SD can be seen in Figs. 7b and 7c,respectively, appearing qualitatively identical. Results of numerically-generated entropy (J kg − K − ) for (cid:112) ∈ { , , } are shown in Fig. 8 and tabulated in Tab. 1. For (cid:112) = (cid:112) =
6, similar results for both FR SD and SD were observed, with entropy generation ranging between ± . × − throughout the entiredomain, with the di ff erence in results between the two schemes being negligible at these polynomialorders. However, for the (cid:112) = ∆ s min = − . × − and ∆ s max = . × − for FR SD and ∆ s min = − . × − and ∆ s max = . × − for SD. These results demonstrate reduced numerical entropy generation under SDby a factor of approximately three, indicating more favorable results for this particular under-resolved caseat (cid:112) = (cid:112) + / ( (cid:112) +
1) is greatest for the SD scheme.4 (a) (b) (c)(d) (e) (f)Figure 8: Inviscid, subsonic flow over 2D cylinder: entropy (a) FR SD , (cid:112) =
2; (b) FR SD , (cid:112) =
4; (c) FR SD , (cid:112) =
6; (d) SD, (cid:112) = (cid:112) =
4; (f) SD, (cid:112) = (cid:112) FR SD SD (cid:112) + (cid:112) + ∆ s min ∆ s max ∆ s min ∆ s max − . × − . × − − . × − . × − / = . − . × − . × − − . × − . × − / = . − . × − . × − − . × − . × − / = . − K − ) under FR SD and SD. = SD schemes for under-resolvedsimulations of turbulent flow.The initial conditions of velocity and pressure for the TGV are given by ρ ( x , y , z , = pRT o (52a) u ( x , y , z , = u o sin (cid:18) xL (cid:19) cos (cid:18) yL (cid:19) cos (cid:18) zL (cid:19) (52b) v ( x , y , z , = − u o cos (cid:18) xL (cid:19) sin (cid:18) yL (cid:19) cos (cid:18) zL (cid:19) (52c) w ( x , y , z , = p ( x , y , z , = ρ o u o (cid:32) γ M o + (cid:34) cos (cid:32) xL (cid:33) + cos (cid:32) yL (cid:33)(cid:35) (cid:34) cos (cid:32) zL (cid:33) + (cid:35)(cid:33) (52e)where the reference velocity, density, and Mach number are u o = ρ o =
1, and M o = .
1, respectively. Thequantity L defines a length scale for the problem; Reynolds number is defined as Re = ( ρ o u o L ) /µ , and isset at 1 600. The fluid is modeled as a perfect gas with a specific heat ratio of γ = . Pr = .
71. From the ideal gas law RT o = p o /ρ o , and if we initialize the flow field with the assumptionof isothermal flow, then p /ρ = p o /ρ o . This relationship allows the initial density field to be set accordingto Eq. (52a). The flow is computed inside a square domain Ω = { x , y , z | (cid:54) x , y , z (cid:54) π L } with periodicboundaries. A characteristic convective time scale can be defined as t c = L / u o . The non-dimensionalintegrated kinetic energy is K = ρ o u o V (cid:90) Ω ρ u · u d x (53)where V is the total volume of the domain and d x = d x d y d z . For this test case we choose L = V = π . The principal method of testing turbulent flow simulation methodologies6using the TGV test case is to compute and track the dissipation rate of the kinetic energy through time.The dissipation rate based upon the kinetic energy is (cid:15) ( K ) = − d K d t (cid:63) (54)where t (cid:63) = tu o / L . The non-dimensional integrated enstrophy is ζ = t c ρ o V (cid:90) Ω ρ ω · ω d x . (55)For strictly incompressible flow, the non-dimensional theoretical vorticity-based dissipation rate isproportional to ζ by (cid:15) ( ζ ) = µρ o u o t c ζ. (56)In a compressible fluid, the non-dimensional theoretical dissipation rate is based upon the summation ofthe following three terms (cid:15) ( S d ) = µ t c ρ o u o V (cid:90) Ω S d : S d d x , (57a) (cid:15) ( p ) = − t c ρ o u o V (cid:90) Ω p ∇ · u d x , (57b) (cid:15) ( µ b ) = µ b t c ρ o u o V (cid:90) Ω ( ∇ · u ) d x (57c)where (cid:15) ( S d ) and (cid:15) ( p ) are the dissipation terms based upon the deviatoric strain-rate tensor S d and pressuredilatation, respectively. Under Stokes’ hypothesis, the bulk viscosity µ b is assigned a value of zero, whichleads to the second coe ffi cient of viscosity taking the value λ = − / µ ; therefore, the dissipation due to thebulk viscosity is neglected. Furthermore, for low Mach number flows with negligible compressibilitye ff ects, the theoretical dissipation rate reasonably approximates the integrated enstrophy and can beestimated by (cid:15) ( S d ). In these simulations, we compute the theoretical dissipation rate as (cid:15) ( S d ) + (cid:15) ( p ). Allintegrals are approximated with a su ffi ciently high-strength quadrature rule. The measured dissipation rate7 (cid:15) ( K ) is computed during post-processing using second-order finite di ff erences to approximate the temporalderivative of the kinetic energy. A reference solution has been provided by van Rees et al. [72], which hasto be scaled by a factor of 1 / V to match the presentation of the current results. These authors performeda direct numerical simulation (DNS) at Re = . (a) (b) (c)Figure 9: TGV: SD result of Q -criterion ( QL / u o = .
5) colored by velocity magnitude at (a) t (cid:63) =
5, (b) t (cid:63) =
11 and (c) t (cid:63) = grid using (cid:112) = DoF).
First, we perform well-resolved simulations of the TGV using a 64 grid and (cid:112) =
3, giving a total of256 DoF, to show the ability of both SD and FR SD to accurately capture the flow physics of the TGVand its transition to and subsequent decaying of turbulence. All simulations for this test case are run usingDavis’ form of the Rusanov approximate Riemann solver such that a close comparison can be made to theresults from Vermeire et al. [73] who used FR DG with similar initial conditions. Figure 9 demonstrates theroll-up of the vortex sheets at t (cid:63) =
5, the transition to turbulence leading to the production of small-scalevortical structures at t (cid:63) =
11, and the subsequent decaying of these structures depicted at t (cid:63) =
20. Resultsof (cid:15) ( K ) and (cid:15) ( S d ) + (cid:15) ( p ) in Fig. 10a and Fig. 10b indicate little discrepancy between the measured andtheoretical dissipation rates, with the peak dissipation rate occurring near t (cid:63) =
9. The actual di ff erencebetween (cid:15) ( K ) and (cid:15) ( S d ) + (cid:15) ( p ) is plotted in Fig. 10c and can be attributed to numerical dissipation anddispersion, non-conservation in evaluating the derivative of the conservative variables since the schemeis only guaranteed to be C continuous [19], and numerical errors aliased from the higher modes to the8lower ones. We can observe that the maximum di ff erence under SD is approximately 60% of that exhibitedunder FR SD . The pressure dilatation-based dissipation rate—which measures compressibility e ff ects onthe dissipation of turbulent energy—among the two schemes is essentially identical and shown in Fig. 10d.Maximum values of (cid:15) ( p ) are approximately 2 × − .Following the procedure laid out in Brachet et al. [67], we compute the spherically-averaged energyspectra E ( κ ) at the peak dissipation rate ( t (cid:63) = SD exhibit an accumulation of energy near the cuto ff wavenumber κ =
128 due to the dissipation inherent to the Riemann solver [20]. Sharp dissipation is known to promotethis pile-up of energy prior to the dissipation range and induce a more pronounced bottleneck e ff ect [74].This build-up of energy at the smallest captured scales is related to contamination of the true physics bynumerical errors such as dispersion.9 (a) (b)(c) (d)Figure 10: TGV: (a) measured dissipation rate based on kinetic energy (cid:15) ( K ), (b) theoretical dissipation based on strain-rate (cid:15) ( S d ) and pressure dilatation (cid:15) ( p ), (c) di ff erence between (a) and (b) (cid:15) ( K ) − (cid:15) ( S d ) − (cid:15) ( p ), (d) pressure dilatation (cid:15) ( p ). Resultsare from a 64 grid using (cid:112) = DoF). DNS results have been provided by van Rees et al. [72]. Figure 11: TGV: energy spectra at t (cid:63) = grid using (cid:112) = DoF). The cuto ff wavenumber ( − ) and the -5 / −− ) are plotted in gray. DNS results have been provided by van Rees et al. [72]. We perform under-resolved simulations of the TGV using an 8 grid while increasing (cid:112) to see the e ff ectof higher polynomial orders on stability for true spectral di ff erence and the modified spectral di ff erencerecovered via the flux reconstruction formulation. We start the simulations at (cid:112) = (cid:112) =
8. Therefore,we are considering seven di ff erent levels of resolution: 24 , 32 , 40 , 48 , 56 , 64 and 72 DoF. To reducethe amount of numerical dissipation, we run all simulations for this test case using Roe’s scheme [75] forthe approximate Riemann solver. Results of (cid:15) ( K ) and (cid:15) ( S d ) + (cid:15) ( p ) are plotted in Fig. 12. In Fig. 12a, weobserve a large amount of numerical dissipation in the results computed using (cid:112) =
3, whereby the rateof kinetic energy loss is overestimated at earlier times in the simulation, where the flow is restricted to asmaller range of scales. The simulation from FR SD is quickly rendered unstable at (cid:112) =
4, largely due toaliasing errors produced at the higher wavenumbers when substantial roll-up of the vortex sheets occursnear t (cid:63) = (cid:112) . The simulationsfrom the SD scheme, on the other hand, demonstrate that as (cid:112) is increased further, the solution is stable andthe di ff erence between the measured dissipation rate due to kinetic energy and the theoretical dissipationrate becomes smaller, and the result from (cid:15) ( S d ) + (cid:15) ( p ) approaches the DNS result up to (cid:112) =
7. However,the SD solution does become unstable at (cid:112) = t (cid:63) =
5. Overall, these results indicate suppressed1aliasing errors in and enhanced stability of the SD scheme on coarse grids with higher polynomial orderswhen performing under-resolved turbulence simulations without any filtering, subgrid-scale modeling, orde-aliasing.2 (a) (b)(c) (d)(e) (f)Figure 12: TGV: measured dissipation rate based on kinetic energy (cid:15) ( K ) ( black ) and theoretical dissipation rate based on strain-rate and pressure dilatation (cid:15) ( S d ) + (cid:15) ( p ) ( red ) on a 8 grid; (a) (cid:112) =
3, (b) (cid:112) =
4, (c) (cid:112) =
5, (d) (cid:112) =
6, (e) (cid:112) =
7, (f) (cid:112) =
8. DNSresults have been provided by van Rees et al. [72]. =
60 000 , α = ◦ We perform implicit large eddy simulations of the transitional flow of a Selig–Donovan (SD) 7003airfoil [76, 77] at Re =
60 000, Mach number M = . α = ◦ . This test case iscommonly used to assess a numerical scheme’s ability to predict separation and transition in a turbulentflow [78, 79, 16, 80, 17], and we compare results from the flux reconstruction and spectral di ff erenceschemes without any filtering, subgrid-scale modeling, or de-aliasing. Laminar flow separation andreattachment occurs on the upper surface of the airfoil, forming a laminar separation bubble (LSB) nearthe leading edge. Lift and drag on an airfoil can be significantly a ff ected by an LSB, which can causestability and control issues. The flow experiences transition near reattachment in the unsteady solution,which causes a region of turbulence over a large portion of the airfoil’s upper surface and a turbulent wakedownstream of the airfoil. (a) (b)Figure 13: SD7003 at Re =
60 000, α = ◦ : (a) near wall region of mesh A, provided by Vermeire et al. [73] (b) near wall regionof mesh B. Figure 14: SD7003 at Re =
60 000, α = ◦ : isosurface of Q -criterion ( Qc / u ∞ = (cid:112) = To perform these simulations, we use two meshes of di ff erent resolution as shown in Fig. 13, the first(mesh A) of which was provided by Vermeire et al. [73]. We use these two di ff erent meshes to study theability of each scheme to simulate under-resolved transitional and turbulent flow at varying levels of (cid:112) .Mesh A contains a total of 137 916 hexahedral elements with 12 elements in the spanwise direction. Thedomain extends 10 c upstream and 20 c downstream of the airfoil and extends in the spanwise directionby 0 . c , where c is the chord length. This spanwise length is deemed su ffi cient for capturing spanwisestructures [16]. We use this mesh to verify our implementation and directly compare results to those from awell-established FR implementation in PyFR [73]. For this mesh, we set (cid:112) = . × DoF. The second mesh constructed (mesh B) is acoarser mesh that contains a total of 33 264 elements with 8 elements in the spanwise direction, whichprovides roughly the same number of degrees of freedom (1 . × ) using (cid:112) =
7. The upper surface ofthe airfoil in mesh A and mesh B is represented with 173 and 110 elements along the chord, respectively.This gives a total of 5 . × DoF on the upper surface in mesh A and 5 . × DoF in mesh B.To better capture the solid boundary curvature, the airfoil surface is represented by quartic elements. Ano-slip adiabatic boundary condition is used for the airfoil surface, Riemann invariant boundary conditionsare applied to the far field, and periodic conditions are applied in the spanwise direction. We use the low-storage, four-stage, third-order embedded pair time integration scheme (RK[4,3(2)]-2N) with adaptivetime-stepping to integrate in time. We march forward in time for 30 t c , where t c = c / u ∞ is one convective5time period. At 20 t c the flow is considered fully developed, and we collect time and spanwise-averagestatistics between 20 t c and 30 t c .Figure 14 displays an isosurface of the Q -criterion ( Qc / u ∞ = (cid:112) =
7. Time and spanwise-averaged plots of the pressure and skin friction coe ffi cientsare shown in Fig. 15. We report maximum skin friction values in the turbulent region above the airfoil usingSD of 8 . × − (mesh A, (cid:112) =
4) and 8 . × − (mesh B, (cid:112) =
7) and FR DG of 7 . × − (mesh A, (cid:112) = y + values of 8.95, 12.54 and 8.40, respectively. However, the corresponding y + values ofthe first solution point nearest the airfoil surface, y + | sp , are 0.42, 0.25 and 0.39. Table 2 demonstrates thataveraged values of the lift coe ffi cient C L and drag coe ffi cient C D as well as time and spanwise-averagedvalues of flow separation x s / c and reattachment x r / c locations of the laminar separation bubble are inagreement with various discontinuous spectral element results of implicit large eddy simulation found inthe literature. The ILES results from Garmann et al. [80], who used a 6th order finite di ff erence scheme,are also provided in the table. We report here that under SD, the simulation is stable on both the coarsemesh ( (cid:112) =
7) and fine mesh ( (cid:112) = DG , the simulation is rendered unstable only on the coarsemesh, and under FR SD , the simulation is unstable on both meshes. These findings demonstrate the extrastability a ff orded by the staggered arrangement of flux points inherent to the SD scheme for achieving astable under-resolved implicit large eddy simulation of transitional flow using a higher polynomial orderon a coarse grid. Furthermore, in light of the FR results for this test case, we recommend the use ofFR DG instead of FR SD when filtering or de-aliasing is not applied for these under-resolved simulations ofturbulent flows.6 (a) (b)Figure 15: SD7003 at Re =
60 000, α = ◦ : (a) pressure coe ffi cient C p , (b) upper surface skin friction coe ffi cient C f .Results corresponding to (cid:112) = (cid:112) = (cid:112) C L C D x s / c x r / cCurrent SD A 137 916 4 0.938 0.049 0.032 0.317FR SD A 137 916 4 (cid:55) (cid:55) (cid:55) (cid:55) FR DG A 137 916 4 0.942 0.051 0.031 0.330SD B 33 264 7 0.940 0.048 0.028 0.301FR SD B 33 264 7 (cid:55) (cid:55) (cid:55) (cid:55) FR DG B 33 264 7 (cid:55) (cid:55) (cid:55) (cid:55)
Vermeire et al. [73] FR DG A 137 916 4 0.941 0.049 0.045 0.315Romero [81] DFR - 202 500 4 0.950 0.045 0.035 -Beck et al. [17] DGSEM - 66 500 3 0.923 0.045 0.027 0.310Beck et al. [17] DGSEM - 8 900 7 0.932 0.050 0.030 0.336Garmann et al. [80] FD (6th order) - 12 549 120 - 0.969 0.039 0.023 0.259Table 2: SD7003 at Re =
60 000, α = ◦ : averaged lift coe ffi cient C L , drag coe ffi cient C D , separation location x s / c , andreattachment location x r / c . Unstable simulations are indicated by the symbol (cid:55) . Results from various authors are providedfor reference. Performance of the spectral di ff erence and the flux reconstruction schemes was measured using thesimulations on mesh A in terms of wall-clock time taken to compute the divergence of the flux ∇ · F = ∂ x f + ∂ y g + ∂ z h , normalized by the total degrees of freedom, number of equations to solve, and number7of stages k in the time stepping scheme, such that t (cid:63) wall = t wall / DoF / N eq / k . All simulations have beendone using double precision. The results shown in Tab. 3 demonstrate that, with the current high-orderframework of the solver, the performance of the spectral di ff erence and flux reconstruction schemes isapproximately identical on mesh A using (cid:112) = ff er a complimentary andmore supportive view on the e ffi ciency of the spectral di ff erence scheme. Scheme t (cid:63) wall (1 × − s)SD 0.5920FR 0.5924Table 3: Wall-clock time to compute ∇ · F in mesh A using 48 Intel Xeon E5-2680 v4processors, normalized by total degrees of freedom, number of equations, and numberof RK stages. All calculations are done using double precision.
6. Conclusions
We reported the development of various discontinuous spectral element methods within a single high-order coding framework such that a fair and impartial comparison among several numerical schemesmay be performed—most notably the true spectral di ff erence and flux reconstruction methods. Withthis construct, we were able to assess the accuracy, stability, and performance of these two schemes.Furthermore, we provided a novel nonlinear stability analysis of the spectral di ff erence scheme anddemonstrated that the error bound for this scheme can be smaller than the flux reconstruction scheme dueto the staggered nature of the flux points. We performed a number of numerical experiments to supportthis analysis, such as heterogeneous linear advection, isentropic Euler vortex, inviscid, subsonic flow overa cylinder, Taylor–Green vortex at Re = Re =
60 000. Theseresults highlighted the advantages of using the baseline SD scheme on coarse grids with higher polynomialorders and demonstrated the potential for extra stability a ff orded by the SD scheme in achieving stableunder-resolved implicit large eddy simulations of turbulent flow. Based on both numerical analysis andexperiments, we can conclude that the pure spectral di ff erence method can be more robust for nonlinearproblems than its flux reconstruction analog, incurring less of a need for de-aliasing.8 Acknowledgments
We would like to thank the support received under the Texas A&M Chancellor’s Research Initiative forpartially funding this work. We would also like thank Guido Lodato for helpful discussions and for sharingresults from his spectral di ff erence flow solver. Lastly, we thank the Texas Advanced Computing Centerand Texas A&M University’s High Performance Research Computing facility for providing the resourcesto perform these simulations. References [1] W. Reed, T. Hill, Triangular mesh methods for the neutron transport equation, Tech. Rep. LA-UR-73-479, Los Alamos Scientific Laboratory (1973).[2] J.-B. Chapelier, M. D. L. L. Plata, F. Renac, Inviscid and viscous simulations of the Taylor-Greenvortex flow using a modal discontinuous Galerkin approach, in: 42nd AIAA Fluid DynamicsConference and Exhibit, American Institute of Aeronautics and Astronautics, 2012. doi:10.2514/6.2012-3073 .[3] G. Karniadakis, S. Sherwin, Spectral / hp Element Methods for Computational Fluid Dynamics, 2ndEdition, Oxford University Press, 2005.[4] E. Tadmor, Shock capturing by the spectral viscosity method, Computer Methods in AppliedMechanics and Engineering 80 (1-3) (1990) 197–208. doi:10.1016/0045-7825(90)90023-f .[5] E. Tadmor, Super viscosity and spectral approximation of nonlinear conservation laws, ClarendonPress Oxford University Press, Oxford New York, 1993, Ch. 5, pp. 69–82.[6] G.-S. Karamanos, G. Karniadakis, A spectral vanishing viscosity method for large-eddy simulations,Journal of Computational Physics 163 (1) (2000) 22–50. doi:10.1006/jcph.2000.6552 .[7] Y. Maday, S. M. O. Kaber, E. Tadmor, Legendre pseudospectral viscosity method for nonlinearconservation laws, SIAM Journal on Numerical Analysis 30 (2) (1993) 321–342. doi:10.1137/0730016 .[8] G. J. Gassner, A. D. Beck, On the accuracy of high-order discretizations for underresolved turbulencesimulations, Theoretical and Computational Fluid Dynamics 27 (3-4) (2012) 221–237. doi:10.1007/s00162-011-0253-7 .9[9] J. Hesthaven, T. Warburton, Nodal Discontinuous Galerkin Methods - Algorithms, Analysis, andApplications, Springer-Berlin, 2008.[10] D. Gottlieb, C.-W. Shu, On the Gibbs phenomenon and its resolution, SIAM Review 39 (4) (1997)644–668. doi:10.1137/s0036144596301390 .[11] D. Gottlieb, J. Hesthaven, Spectral methods for hyperbolic problems, Journal of Computational andApplied Mathematics 128 (1-2) (2001) 83–131. doi:10.1016/s0377-0427(00)00510-0 .[12] A. Kravchenko, P. Moin, On the e ff ect of numerical errors in large eddy simulations of turbulentflows, Journal of Computational Physics 131 (2) (1997) 310–322. doi:10.1006/jcph.1996.5597 .[13] G. Blaisdell, E. Spyropoulos, J. Qin, The e ff ect of the formulation of nonlinear terms on aliasingerrors in spectral methods, Applied Numerical Mathematics 21 (3) (1996) 207–219. doi:10.1016/0168-9274(96)00005-0 .[14] G. J. Gassner, A. R. Winters, D. A. Kopriva, Split form nodal discontinuous Galerkin schemes withsummation-by-parts property for the compressible euler equations, Journal of Computational Physics327 (2016) 39–66. doi:10.1016/j.jcp.2016.09.013 .[15] A. Winters, R. Moura, G. Mengaldo, G. Gassner, S. Walch, J. Peiro, S. Sherwin, A comparativestudy on polynomial dealiasing and split form discontinuous Galerkin schemes for under-resolvedturbulence computations, Journal of Computational Physics 372 (2018) 1–21.[16] A. Uranga, P.-O. Persson, M. Drela, J. Peraire, Implicit large eddy simulation of transition toturbulence at low Reynolds numbers using a discontinuous Galerkin method, International Journalfor Numerical Method in Fluids (87) (2011) 232–261.[17] A. Beck, T. Bolemann, D. Flad, H. Frank, G. Gassner, F. Hindenlang, C.-D. Munz, High-orderdiscontinuous Galerkin spectral element methods for transitional and turbulent flow simulations,International Journal for Numerical Methods in Fluids 76 (2014) 522–548.[18] C. C. de Wiart, K. Hillewaert, DNS and ILES of transitional flows around a SD7003 using a highorder discontinuous Galerkin method, in: Seventh International Conference on Computational FluidDynamics, Big Island, HI, 2012.[19] B. Vermeire, S. Nadarajah, P. Tucker, Implicit large eddy simulation using the high-order correctionprocedure via reconstruction scheme, International Journal for Numerical Methods in Fluids 82(2016) 231–260.[20] R. Moura, G. Mengaldo, J. Peiró, S. Sherwin, On the eddy-resolving capability of high-orderdiscontinuous Galerkin approaches to implicit LES / under-resolved DNS of Euler turbulence,0 Journal of Computational Physics 330 (2017) 615–623.[21] H. Huynh, A flux reconstruction approach to high-order schemes including discontinuous Galerkinmethods, in: 18th AIAA Computational Fluid Dynamics Conference, Miami, FL, 2007.[22] P. Vincent, P. Castonguay, A. Jameson, A new class of high-order energy stable flux reconstructionschemes, Journal of Scientific Computing 47 (1) (2011) 50–72.[23] F. Witherden, B. Vermeire, P. Vincent, Heterogeneous computing on mixed unstructured grids withPyFR, Computers and Fluids 120 (2015) 173–186.[24] D. Williams, P. Castonguay, P. Vincent, A. Jameson, An extension of energy stable fluxreconstruction to unsteady, non-linear, viscous problems on mixed grids, in: 20th AIAAComputational Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics,2011. doi:10.2514/6.2011-3405 .[25] P. Castonguay, P. Vincent, A. Jameson, Application of high-order energy stable flux reconstructionschemes to the Euler equations, in: 49th AIAA Aerospace Sciences Meeting including the NewHorizons Forum and Aerospace Exposition, American Institute of Aeronautics and Astronautics,2011. doi:10.2514/6.2011-686 .[26] P. Castonguay, D. Williams, P. Vincent, A. Jameson, Energy stable flux reconstruction schemesfor advection–di ff usion problems, Computer Methods in Applied Mechanics and Engineering 267(2013) 400–417. doi:10.1016/j.cma.2013.08.012 .[27] C. Cox, C. Liang, M. Plesniak, Spectral di ff erence solution of incompressible flow over an inlinetube bundle with oscillating cylinder, in: ASME Pressure Vessels and Piping Conference, Toronto,Ontario, Canada, 2012, pp. 9–20.[28] N. Loppi, F. Witherden, A. Jameson, P. Vincent, A high-order cross-platform incompressibleNavier–Stokes solver via artificial compressibility with application to a turbulent jet, ComputerPhysics Communications 233 (2018) 193–205. doi:10.1016/j.cpc.2018.06.016 .[29] F. Witherden, A. Farrington, P. Vincent, PyFR: An open source framework for solving advection-di ff usion type problems on streaming architectures using the flux reconstruction approach, ComputerPhysics Communications 185 (2014) 3028–3040.[30] P. Vincent, F. Witherden, B. Vermeire, J. S. Park, A. Iyer, Towards green aviation with Pythonat petascale, in: SC16: International Conference for High Performance Computing, Networking,Storage and Analysis, IEEE, 2016. doi:10.1109/sc.2016.1 .[31] D. A. Kopriva, A conservative staggered-grid Chebyshev multidomain method for compressible1 flows. II. A semi-structured method, Journal of Computational Physics 128 (2) (1996) 475–488. doi:10.1006/jcph.1996.0225 .[32] D. A. Kopriva, J. H. Kolias, A conservative staggered-grid Chebyshev multidomain method forcompressible flows, Journal of Computational Physics 125 (1) (1996) 244–261. doi:10.1006/jcph.1996.0091 .[33] A. Jameson, A proof of the stability of the spectral di ff erence method for all orders of accuracy,Journal of Scientific Computing 45 (2010) 348–358.[34] Z. Wang, Y. Liu, G. May, A. Jameson, Spectral di ff erence method for unstructured grids II: Extensionto the Euler equations, Journal of Scientific Computing 32 (2007) 45–71.[35] Y. Sun, Z. Wang, Y. Liu, High-order multidomain spectral di ff erence method for the Navier-Stokesequations on unstructured hexahedral grids, Communications in Computational Physics 2 (2) (2007)310–333.[36] M. Yu, Z. Wang, H. Hu, A high-order spectral di ff erence method for unstructured dynamic grids,Computers and Fluids 48 (2011) 84–97.[37] A. Chan, P. Dewey, A. Jameson, C. Liang, A. Smits, Vortex suppression and drag reduction in thewake of counter-rotating cylinders, Journal of Fluid Mechanics 679 (2011) 343–382.[38] G. Lodato, L. Vervisch, P. Clavin, Direct numerical simulation of shock wavy-wall interaction:analysis of cellular shock structures and flow patterns, Journal of Fluid Mechanics 789 (2016) 221–258. doi:10.1017/jfm.2015.731 .[39] G. Lodato, Characteristic modal shock detection for discontinuous finite element methods,Computers and Fluids 179 (2019) 309–333.[40] J. Wang, M. Miesch, C. Liang, Convection in oblate solar-type stars, The Astrophysical Journal830 (1) (2016).[41] Y. Liu, M. Vinokur, Z. Wang, Spectral di ff erence method for unstructured grids I: Basic formulation,Journal of Computational Phyics 216 (2006) 780–801.[42] F. Bassi, S. Rebay, High-order accurate discontinuous finite element solution of the 2D Eulerequations, Journal of Computational Physics 138 (1997) 251–285.[43] H. Huynh, A reconstruction approach to high-order schemes including discontinuous Galerkin fordi ff usion, in: 47th AIAA Aerospace Sciences Meeting, Orlando, FL, 2009.[44] R. Radau, Étude sur les formules d’approximation qui servent à calculer la valeur numérique d’uneintégrale définie., Journal de Mathématiques Pures et Appliquées 6 (1880) 283–336.2[45] J. Romero, K. Asthana, A. Jameson, A simplified formulation of the flux reconstruction method,Journal Scientific Computing 67 (2016) 351–374.[46] H. Huynh, Discontinuous Galerkin via interpolation: the direct flux reconstruction method, Journalof Scientific Computing 82 (75) (2020).[47] F. Bassi, S. Rebay, A high-order discontinuous Galerkin method for compressible turbulent flows,in: G. K. B. Cockburn, C.-W. Shu (Eds.), Discontinuous Galerkin methods: Theory, Computation,and Application. Lecture notes in Computational Science and Engineering, Vol. 11, Springer, 2000,pp. 77–88.[48] A. Jameson, P. Castonguay, P. Vincent, On the non-linear stability of flux reconstruction schemes,Journal of Scientific Computing 50 (2) (2012) 434–445.[49] P. D. Lax, Hyperbolic systems of conservation laws II, Communications on Pure and AppliedMathematics 10 (4) (1957) 537–566. doi:10.1002/cpa.3160100406 .[50] S. Osher, Riemann solvers, the entropy condition, and di ff erence, SIAM Journal on NumericalAnalysis 21 (2) (1984) 217–235. doi:10.1137/0721016 .[51] C. Bernardi, Y. Maday, Spectral methods, in: Handbook of Numerical Analysis, Elsevier, 1997, pp.209–485. doi:10.1016/s1570-8659(97)80003-8 .[52] P. Vincent, P. Castonguay, A. Jameson, Insights from von Neumann analysis of high-order fluxreconstruction schemes, Journal of Computational Physics 230 (22) (2011) 8134–8154.[53] K. Asthana, J. Watkins, A. Jameson, On consistency and rate of convergence of flux reconstructionfor time-dependent problems, Journal of Computational Physics 334 (2017) 367–391. doi:10.1016/j.jcp.2017.01.008 .[54] W. Trojak, Generalised Sobolev Stable Flux Reconstruction (2018) 1–16 arXiv:1804.04714 .[55] M. Alhawwary, Z. Wang, Fourier analysis and evaluation of DG, FD and compact di ff erencemethods for conservation laws, Journal of Computational Physics 373 (2018) 835–862. doi:10.1016/j.jcp.2018.07.018 .[56] O. San, Analysis of low-pass filters for approximate deconvolution closure modelling in one-dimensional decaying Burgers turbulence, International Journal of Computational Fluid Dynamics30 (1) (2016) 20–37. doi:10.1080/10618562.2016.1155705 .[57] B. Cockburn, M. Luskin, C.-W. Shu, E. Süli, Post-processing of Galerkin methods for hyperbolicproblems, in: Lecture Notes in Computational Science and Engineering, Springer Berlin Heidelberg,2000, pp. 291–300. doi:10.1007/978-3-642-59721-3_24 .3[58] S. Adjerid, K. D. Devine, J. E. Flaherty, L. Krivodonova, A posteriori error estimation fordiscontinuous Galerkin solutions of hyperbolic problems, Computer Methods in Applied Mechanicsand Engineering 191 (11-12) (2002) 1097–1112. doi:10.1016/s0045-7825(01)00318-8 .[59] C.-W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes forhyperbolic conservation laws, in: A. Quarteroni (Ed.), Advanced Numerical Approximation ofNonlinear Hyperbolic Conservation Laws, Vol. 1697 of Lecture Notes in Mathematics, Springer,1998, pp. 325–432.[60] S. Davis, Simplified second-order Godunov-type methods, SIAM Journal on Scientific and StatisticalComputing 9 (3) (1988) 445–473.[61] M. Carpenter, C. Kennedy, Fourth-order 2N-storage Runge-Kutta schemes, Tech. Rep. 109112,NASA (1994).[62] S. Spiegel, H. Huynh, J. DeBonis, De-aliasing through over-integration applied to the fluxreconstruction and discontinuous Galerkin methods, in: 22nd AIAA Aviation Forum, Dallas, TX,2015.[63] G. Mengaldo, D. D. Grazia, D. Moxey, P. Vincent, S. Sherwin, Dealiasing techniques for high-order spectral element methods on regular and irregular grids, Journal of Computational Physics 299(2015) 56–81.[64] M. Carpenter, C. Kennedy, Third-order 2N-storage Runge-Kutta schemes with error control, Tech.Rep. 109111, NASA (1994).[65] C. Kennedy, M. Carpenter, R. Lewis, Low-storage, explicit Runge-Kutta schemes for thecompressible Navier-Stokes equations, Applied Numerical Mathematics 35 (3) (2000) 177–219.[66] G. Taylor, A. Green, Mechanism of the production of small eddies from large ones, Proceedingsof the Royal Society of London. Series A - Mathematical and Physical Sciences 158 (895) (1937)499–521.[67] M. Brachet, D. Meiron, S. Orszag, B. Nickel, R. Morf, U. Frisch, Small-scale structure of the Taylor-Green vortex, Journal of Fluid Mechanics 130 (1983) 411–452.[68] Z. Wang, K. Fidkowski, R. Abgrall, F. Bassi, D. Caraeni, A. Cary, H. Deconinck, R. Hartmann,K. Hillewaert, H. Huynh, N. Kroll, G. May, P.-O. Persson, B. van Leer, M. Visbal, High-order CFDmethods: current status and perspective, International Journal for Numerical Methods in Fluids 72(2013) 811–845.[69] G. Gassner, A. Beck, On the accuracy of high-order discretizations for underresolved turbulence4 simulations, Theoretical Computational Fluid Dynamics 27 221–237.[70] C. C. de Wiart, K. Hillewaert, M. Duponcheel, G. Winckelmans, Assessment of a discontinuousGalerkin method for the simulation of vortical flows at high Reynolds number, International Journalfor Numerical Methods in Fluids 74 (2014) 469–493.[71] J. Bull, A. Jameson, Simulation of the Taylor-Green vortex using high-order flux reconstructionschemes, AIAA Journal 53 (9) (2015).[72] W. van Rees, A. Leonard, D. Pullin, P. Koumoutsakos, A comparison of vortex and pseudo-spectral methods for the simulation of periodic vortical flows at high Reynolds numbers, Journalof Computational Physics 230 (2011) 2794–2805.[73] B. Vermeire, F. Witherden, P. Vincent, On the utility of GPU accelerated high-order methods forunsteady flow simulations: a comparison with industry-standard tools, Journal of ComputationalPhysics 334 (2017) 497–521.[74] G. Falkovich, Bottleneck phenomena in developed turbulence, Physics of Fluids 6 (4) (1994).[75] E. Toro, Riemann solvers and numerical methods for fluid dynamics, 3rd Edition, Springer, 1999.[76] M. Selig, J. Donovan, D. Fraser, Airfoils at Low Speeds, Vol. 1, H.A. Stokely, Virginia Beach, VA,USA, 1989.[77] M. Selig, J. Guglielmo, A. Broeren, P. Giguére, Summary of Low-Speed Airfoil Data, Vol. 1,SoarTech, Virginia Beach, VA, USA, 1995.[78] M. Visbal, R. Gordnier, M. Galbraith, High-fidelity simulations of moving and flexible airfoils at lowReynolds numbers, Experiments in Fluids 46 (5) (2009) 903–922.[79] M. Galbraith, M. Visbal, Implicit large eddy simulation of low-Reynolds-number transitional flowpast the SD7003 airfoil, in: 40th Fluid Dynamics Conference and Exhibit, Chicago, IL, 2010.[80] D. Garmann, M. Visbal, P. Orkwis, Comparative study of implicit and subgrid-scale model large-eddy simulation techniques for low-Reynolds number airfoil applications, International Journal forNumerical Methods in Fluids 71 (2013) 1546–1565.[81] J. Romero, On the development of the direct flux reconstruction scheme for high-order fluid flowsimulations, Ph.D. thesis, Stanford University, Stanford, CA (2017).[82] C. Liang, C. Cox, M. Plesniak, A comparison of computational e ffi ciencies of spectral di ffff