On the properties of energy stable ﬂux reconstruction schemes for implicit large eddy simulation

We begin by investigating the stability, order of accuracy, and dispersion and dissipation characteristics of the extended range of energy stable ﬂux reconstruction (E-ESFR) schemes in the context of implicit large eddy simulation (ILES). We proceed to demonstrate that subsets of the E-ESFR schemes are more stable than collocation nodal discontinuous Galerkin methods recovered with the ﬂux reconstruction approach (FRDG) for marginally- resolved ILES simulations of the Taylor–Green vortex. These schemes are shown to have reduced dissipation and dispersion errors relative to FRDG schemes of the same polynomial degree and, simultaneously, have increased Courant–Friedrichs–Lewy (CFL) limits. Finally, we simulate turbulent ﬂow over an SD7003 aerofoil using two of the most stable E- ESFR schemes identiﬁed by the aforementioned Taylor–Green vortex experiments. Results demonstrate that subsets of E-ESFR schemes appear more stable than the commonly used FRDG method, have increased CFL limits, and are suitable for ILES of complex turbulent ﬂows on unstructured grids. 2016 The an access wavenumbers, while less dissipation error at high wavenumbers. Also, they have reduced group velocities for the majority of wavenumbers. Importantly, they have signiﬁcantly less dispersion error for all but the highest resolved wavenumbers, when compared to their corresponding reference FRDG schemes. These results demonstrate that modifying the correction


Introduction
The flux reconstruction (FR) approach, first introduced by Huynh [1,2], is a high-order accurate numerical method that can be used with mixed-element unstructured grids. It was subsequently extended to two-dimensional simplex elements as the lifting collocation penalty (LCP) formulation [3], and to three-dimensions as the correction procedure via reconstruction (CPR) scheme [4]. High-order methods, such as the FR approach, are particularly appealing for a variety of advection-diffusion problems including the Euler and Navier-Stokes equations. They have been found to provide more accurate solutions with fewer degrees of freedom and reduced computational cost, when compared to conventional lowerorder schemes [5,6]. The FR approach has been shown to be particularly accurate for scale-resolving simulations of complex turbulent flows, including both direct numerical simulation (DNS) and large eddy simulation (LES) [6][7][8][9].
Interestingly, FR does not describe a single scheme, but is in fact a family of schemes. Huynh [1] described several schemes that are stable for linear advection, including one equivalent to a collocation nodal discontinuous Galerkin method, henceforth referred to as FRDG, another equivalent to an energy stable spectral difference (SD) method, and the so-called g 2 method. Subsequently, Vincent et al. [10] discovered a continuous range of single-parameter energy stable FR (ESFR) schemes, also referred to as the Vincent-Castonguay-Jameson-Huynh (VCJH) schemes, which are provably energy stable for 1D linear advection. We refer to this original range of schemes as the original ESFR (O-ESFR) schemes. Properties of these

Formulation
The flux reconstructions (FR) approach, first introduced by Huynh [1], can be used to spatially discretise a general conservation law of the form ∂u ∂t where u = u(x, t) is the solution with a given initial distribution u(x, 0) = u 0 , t is time, f = f (u) is the flux, and x is the spatial coordinate. We split the domain into a mesh where n is one out of a total of N e elements in the domain. For the sake of efficiency and simplicity, we perform all operations in a computational space where each element exists on the interval I ∈ [−1, 1]. For the current study we use the linear mapping where x n is the mesh node corresponding to the left face of n , n is the mapping function, and ξ is the location in reference space. This mapping can be inverted according to which is also linear. The solution within each element is represented by a degree k polynomial, which is allowed to be discontinuous at the interface between elements. This polynomial is supported by nodal basis functions generated at k + 1 solution points.
Therefore, the solution within each element in computational space can be approximated as u δ n = k l=0 u δ n,lφ n,l , (5) where û δ n =û δ n (ξ, t) is the polynomial representation of the solution within an element, u δ n,l = u δ n,l (t) is the value of the solution at solution point l, and φ n,l =φ n,l (ξ ) is its corresponding nodal basis function in reference space. For the onedimensional case these basis functions are the Lagrange polynomialŝ φ n,l (ξ ) = k m=0,m =l ξ − ξ n,m ξ n,l − ξ n,m . (6) The discrete form of the system of equations in computational space can then be written as [1,10] ∂û δ where f δ n is a polynomial representation of the flux in computational space.

Inviscid fluxes
The polynomial representation of a discontinuous flux function f D n =f D n (x, t) can be supported by the same polynomial basis as the solution according to [1] where the superscript D denotes that this flux, like the solution, is allowed to be discontinuous at the interface between elements and f δ n,l = f δ n,l (t) is the flux evaluated at the solution points.
We notice that the flux between elements must be continuous to maintain global conservation. However, since the solution is generally discontinuous between elements, so to is the flux. Huynh [1], proposed that we generate a globally C 0 continuous flux by applying a correction denoted by f δC n to the discontinuous flux in each element such that where f δ n is the globally continuous flux function referred to in Equation (7). Huynh [1] proposed we compute the flux correction aŝ f δC where f δ D n,L =f δ D n (−1, t), and f δ D n,R =f δ D n (1, t).
fluxes computed at the flux points between elements by an appropriate Riemann solver using the extrapolated values u − L , u + L , u − R , and u + R of the solution from the neighbouring elements at each edge. The functions g L = g L (ξ ) and g R = g R (ξ ) are correction functions that are degree k + 1 polynomials with the constraints Since g L and g R are of degree k + 1, so to is the continuous flux function f δ n . Its gradient is [1,10] which is of degree k and in the same polynomial space as û δ n .

Original range of energy stable schemes
Vincent et al. [10] identified a single parameter family of schemes, referred to here as the O-ESFR schemes, that are provably energy stable for linear advection. For these schemes the left correction function is and by symmetry the right correction function is where L k is the degree k Legendre polynomial. These correction functions will yield energy stable schemes provided where and c is a constant such that and Vincent et al. [10] also demonstrated that several well-known schemes can be recovered using this approach including a nodal discontinuous Galerkin scheme when c = 0, an energy stable spectral difference scheme when c = c S D , and Huynh's g 2 scheme when c = c HU .

Extended range of energy stable schemes
Vincent et al. [16] later introduced an extended range of energy stable-symmetric-conservative flux reconstruction correction functions (E-ESFR), of which the O-ESFR schemes are a subset. Like the O-ESFR schemes, all of the E-ESFR schemes are provably energy stable for linear advection. Schemes with k = 1 and k = 2 are single parameter families of schemes and, therefore, are an equivalent set to the O-ESFR schemes. However, for k ≥ 3 the E-ESFR schemes were found to have multiple parameters. For the current study we are interested in schemes with k = 1, 2, 3, and 4, while schemes with k ≥ 5, which have three or more parameters, are left for future work.
Here we will summarise the results of Vincent et al. [16]. First we define D as the nodal differentiation operator according to where φ j is the corresponding nodal basis function at solution point j. Similarly, the Vandermonde matrix V can be formulated as were L j (ξ i ) is a Legendre polynomial of degree j, normalised to unity at ξ = 1. Using the nodal differentiation operator and Vandermonde matrix, we can generate the modal differentiation operator Also, a modal mass matrix M can be generated according tõ We denote gradients of the correction functions given as vectors of modal coefficients for the left and right hand side, respectively, as g ξ L and g ξ R . These can be converted to nodal coefficients by and Vincent et al. [16] proved that FR correction functions are energy stable provided andg and Q is a real square matrix of dimension k + 1 that satisfies They are also symmetric if and conservative if Due to symmetry, Following this procedure yields an under-determined system of equations, which results in families of suitable solutions. All of the correction functions that are members of these families will be energy stable, symmetric, and conservative.

k = 1
Following the procedure outlined by Vincent et al. [16], we generate the k = 1 E-ESFR correction functions, which are found to have one parameter q 0 . The most general form of the correction function that satisfies the given constraints is where g ξ L [i] is the coefficient for mode i of the correction function derivative. These correction functions will yield an energy stable scheme provided Setting q 0 equal to q 0DG = 0 recovers the nodal discontinuous Galerkin scheme, q 0S D = 1/3 recovers the energy stable SD scheme, and q 0HU = 4/3 recovers the g 2 scheme of Huynh. The minimum value of q 0 for which the scheme is stable is

k = 2
Following the procedure outlined by Vincent et al. [16], we generate the k = 2 E-ESFR correction functions, which are found to have one parameter q 0 . The most general form of the correction function that satisfies the given constraints is These correction functions will yield an energy stable scheme provided Setting q 0 equal to q 0DG = 0 recovers the nodal discontinuous Galerkin scheme, q 0S D = 4/15 recovers the energy stable SD scheme, and q 0HU = 3/5 recovers the g 2 scheme of Huynh. The minimum value of q 0 for which the scheme is stable is q 0− = − 2 5 .

k = 3
For k = 3 the extended range of schemes has two parameters q 0 and q 1 [16]. The most general form of the correction function that satisfies the given constraints is These correction functions will yield a energy stable scheme provided By setting q 1 = 0 we recover the O-ESFR single parameter family of schemes. Subsequently, setting q 0 equal to q 0DG = 0 recovers the nodal discontinuous Galerkin scheme, q 0S D = 3/14 recovers the energy stable SD scheme, and q 0HU = 8/21 recovers the g 2 scheme of Huynh. The minimum value of q 0 for which the scheme is stable is q 0− = − 2 7 when q 1 = 0.

k = 4
For k = 4 the extended range of schemes has two parameters q 0 and q 1 [16]. The most general form of the correction function that satisfies the given constraints is These correction functions will yield a energy stable scheme provided 8 15 By setting q 1 = 0 we recover the O-ESFR single parameter family of schemes. Subsequently, setting q 0 equal to q 0DG = 0 recovers the nodal discontinuous Galerkin scheme, q 0S D = 8/45 recovers the energy stable SD scheme, and q 0HU = 5/18 recovers the g 2 scheme of Huynh. The minimum value of q 0 for which the scheme is stable is q 0− = − 2 9 when q 1 = 0.

Overview
To determine the influence of correction function choice on the stability, order of accuracy, and dispersion and dissipation behaviour of the E-ESFR schemes we perform Von Neumann analysis following Huynh [1] and Vincent et al [11]. For this study, we consider the linear advection equation with fully upwind fluxes used at the interfaces between elements.
For Von Neumann analysis we follow the notation of Vincent et al. [11]. For a general linear advection problem which admits plane wave solutions of the form provided the temporal frequency ω = ω(θ) satisfies the dispersion relation and dissipation relation where θ is the wave number and I = √ −1. Following Vincent et al. [11], we consider a mesh where all elements have constant width h = 1. The flux reconstruction scheme can then be cast in matrix-vector form as where û δ is the vector of solution point values, g ξ L is the gradient of the correction function evaluated at the solution points, is the boundary extrapolation operator, and D is the nodal gradient operator. We seek Bloch wave type solution to this of the form where θ δ is a prescribed baseline wavenumber within the range −π ≤ θ ≤ π and ω δ is the resulting temporal frequency of the scheme. The upwind interface flux can then be defined aŝ where r[l] =φ n,l (1).
By substituting Equation (64) and Equation (65) into Equation (62) we obtain where Equation (67) is a classical eigenvalue problem, where v is an eigenvector with one of k + 1 valid eigenvalues ω δ . The operator Q is a function of θ δ , and it follows that so too are v δ and ω δ . We use the approach outlined by Vincent et al. [11] to identify the true wavenumber from the k + 1 admissible solutions. The combined dispersion and dissipation error can then be computed from where ω is the exact temporal wavenumber for the chosen θ δ , grid spacing, and advection speed.

Stability
It has been demonstrated that the form of FR correction function can have a significant effect on the Courant-Friedrichs-Lewy (CFL) time-step limit associated with a given FR scheme [11]. While not the only factor, the CFL limit, which determines the maximum stable explicit time-step, can be an important determinant of overall simulation cost.
To determine CFL limits for the E-ESFR schemes, we cast the system of equations as a full space-time discretisation following Vincent et al. [11]. This can be written aŝ where û δ(m) is the approximate solution at t (m) and û δ(m+1) is the approximate solution at t (m+1) where t (m+1) = t (m) + τ and R is the fully discrete linear operator. For the standard four-stage fourth-order explicit Runge-Kutta scheme (referred to as RK44) this has the form For the scheme to remain stable, the spectral radius of R must remain less than unity for all wavenumbers θ δ in the range −π ≤ θ δ ≤ π . The maximum stable time-step τ C F L is defined as the maximum value of τ that satisfies this condition. Plots of τ C F L as a function of q 0 are shown in Fig. 1 for the O-ESFR k = 1 and k = 2 schemes and two-parameter E-ESFR k = 3 and k = 4 schemes with q 1 = 0, recovering their corresponding O-ESFR schemes. For all schemes, τ C F L tends to 0 as q 0 tends to q 0− . As q 0 increases τ C F L also increases to a maximum at q 0τ , with values provided in Table 1. Interestingly, the k = 1 scheme has a range of values that yield a maximum τ C F L = 2 for q 0 ∈ [3.62, 8.82]. As the polynomial degree is increased the maximum τ C F L decreases, which is consistent with previous studies [1,11].
Contours of τ C F L as a function of q 0 and q 1 are shown in Fig. 2 for the two-parameter k = 3 and k = 4 E-ESFR schemes. There are bands of schemes that yield large values of τ C F L , passing from the bottom left to top right of each plot. The O-ESFR schemes, defined by the line q 1 = 0 do not intersect this band through a global maximum. Instead, the E-ESFR schemes have maxima in regions of both positive q 0 and q 1 . τ C F L was found as q 0 → ∞ with corresponding values of q 1 given in Table 1. Interestingly, the maximum τ C F L for the k = 3 E-ESFR scheme is equal to the maximum τ C F L for the k = 2 O-ESFR scheme. Similarly, the maximum τ C F L for the k = 4 E-ESFR scheme is equal to the maximum τ C F L for the k = 3 O-ESFR scheme. This behaviour can be explained by examining the formulation of the correction functions in the limit q 0 → ∞. For k = 3 it can be shown that Table 1 Values of q 0τ for the O-ESFR schemes and q 0τ and q 1τ for the E-ESFR schemes that yield maximum τ C F L with RK44.
and for k = 4 We can readily compare Equations (72) to (75) for k = 3 with Equations (39) to (41) for k = 2. It is clear that these have an identical form, with only a change of variable from q 0 to q 1 . The same can be shown for Equations (76) to (80) for k = 4, by setting q 1 = 0 in Equations (43) to (46) for k = 3. Therefore, as q 0 → ∞, the k = 3 and k = 4 E-ESFR schemes have correction functions equivalent to the k = 2 and k = 3 O-ESFR schemes with q 1 = 0, respectively, with only a change of variable from q 0 to q 1 . This explains why they achieve the same τ C F L as the O-ESFR schemes that are one polynomial degree lower. However, these schemes fail to energise the highest degree polynomial mode and are unfavourable, as shown by Vincent et al. [10] for the c ∞ schemes described therein.
Nonetheless, the maximum achievable τ C F L of the E-ESFR is greater than those of the same degree O-ESFR schemes by approximately 86% for k = 3 and 55% for k = 4. They are also greater than their corresponding FRDG schemes by 491% and 384%, respectively. The correction functions for the E-ESFR schemes that achieve τ C F L are shown in Fig. 3 alongside the correction functions for the O-ESFR schemes that similarly achieve their corresponding τ C F L and, for reference, the FRDG scheme of the same polynomial degree. For both k = 3 and k = 4 the E-ESFR schemes that achieve τ C F L are shallower than the corresponding O-ESFR schemes, and their maximum gradients are subsequently reduced. This is consistent with the descriptions of Huynh [1], who noted that correction functions with shallower gradients should have less restrictive time-step stability limits. Finally, we note that whilst identification of FR schemes with an increased CFL limit is interesting, and potentially useful, an increased time-step limit is not the only determinant of simulation cost. The overall accuracy of a scheme is clearly also important. Previous studies indicate that FR schemes with larger CFL limits are theoretically less accurate. If finer grids are required to overcome such a reduction in accuracy for real-world problems this could negate any benefits associated with a larger CFL limit.

Order of accuracy
It has been shown previously that higher-order schemes can yield more accurate results on a per degree of freedom basis [6] and, in particular, when performing ILES using the FR scheme [9]. The order of accuracy A T of a scheme, based on dispersion and dissipation errors, can be determined from von Neumann analysis and is defined as [1] where θ δ R is a wavenumber chosen such that E T remains small and is within the well-resolved range. We choose θ δ R = π/16, π/8, π/4, and π/3 for the k = 1, 2, 3, and 4 schemes, respectively. Plots of the order of accuracy as a function of q 0 are  shown in Fig. 4 for single-parameter k = 1 and k = 2 schemes and two-parameter k = 3 and k = 4 schemes with q 1 = 0, recovering the O-ESFR schemes. All of the schemes achieve a level of super accuracy. Each achieves a maximum order of accuracy A T ≈ 2k + 1 at q 0 = q 0DG , which then plateaus to A T ≈ 2k near this point and tends to A T ≈ 2k − 1 as q 0 → ∞.
These results are consistent with the findings of Huynh [1] and Vincent et al. [11], who observed the same behaviour for these families of schemes.
Plots of order of accuracy as a function of q 0 and q 1 for the two-parameter k = 3 and k = 4 E-ESFR schemes are shown in This is apparent from the thin region about q 1 = 0 in the plots, which correspond to the space of the O-ESFR schemes. Again, peak order of accuracy is observed for the FRDG schemes at q 0 = 0 and q 1 = 0 and is A T = 2k + 1. For larger magnitudes of q 1 and nearby q 0 = 0, both schemes achieve an order of accuracy of A T = 2(k − 1). In the limit of q 0 → ∞, the schemes achieve A T = 2(k − 1) + 1 for q 1 = 0, a plateau of A T ≈ 2(k − 1), and A T ≈ 2(k − 1) − 1 as q 1 → ∞ concurrently. This is consistent with these schemes having the same correction functions as O-ESFR schemes one polynomial degree lower, as discussed previously.

Dispersion and dissipation
The ILES hypothesis relies on numerical dissipation concentrated at high wavenumbers to dissipate the smallest resolved turbulent scales [9]. The dispersive behaviour of a scheme is also of primary interest, since it will determine the observed wave-speeds via the group-velocity, and will also affect aliasing driven instabilities for non-linear flux functions. Dispersion and dissipation profiles are shown in Fig. 6 and Fig. 7, respectively, for k = 1 to k = 4 schemes with q 1 = 0, which recovers the O-ESFR schemes. By fixing q 1 , we can examine the influence of varying q 0 on the behaviour of each scheme. We plot profiles with several values of q 0 , specifically q 0 = q 0− /2, q 0− /4, q 0DG , q 0S D , and, q 0HU , allowing us to examine the behaviour of a wide range of schemes. From the dispersion profiles, it is clear that increasing q 0 decreases the dispersion relation at all wavenumbers for each scheme. Also, the group velocities of well-resolved waves is correspondingly decreased.
From the dissipation profiles, we observe that increasing q 0 increases the dissipation error up to θ δ ≈ kπ . After this point it decreases the amount of numerical dissipation. Interestingly, all of the dissipation profiles intersect near the point θ δ ≈ kπ .
Plots of dispersion and dissipation profiles are shown in Fig. 8 and Fig. 9, respectively, for k = 3 and Fig. 10 and Fig. 11 for the k = 4 two-parameter E-ESFR schemes, respectively. Each set of plots are shown for a fixed q 0 , specifically q 0 = q 0− /4, q 0DG , q 0S D , and q 0HU . Profiles are then shown for q 1 = q 1− /2, q 1− /4, 0, q 1+ /4, and q 1+ /2, where q 1− and q 1+ are the minimum and maximum stable values of q 1 for a given q 0 , respectively. Similar to the influence of q 0 , increasing q 1 tends to decrease the dispersion relation, but only up to θ δ ≈ (k − 0.5)π . After this point it tends to increase the dispersion relation up to the highest resolved wavenumbers. Therefore, unlike varying q 0 , varying q 1 has an opposite effect on the low and high wavenumber behaviour of the scheme. Increasing q 1 decreases the numerical dissipation at low wavenumbers and increases numerical dissipation at high wavenumbers, similar to decreasing q 0 . However, it also tends to shift the dissipation curve in general towards lower wavenumbers. This is apparent from the shifting intersection point in the dissipation profile near θ δ ≈ kπ . It is also apparent that large values of q 1 can yield non-physical dispersion profiles, where high-frequency waves are approximately stationary and nearly undamped. This behaviour is similar to schemes with large values of q 0 , as shown by Vincent et al. [11] for c → c ∞ , and is clearly unfavourable.

Taylor-Green vortex
To investigate the behaviour of the E-ESFR schemes for ILES we consider under-resolved simulations of the Taylor-Green vortex test case with k = 3 and k = 4 in the limit Re → ∞ by using the Euler equations, and at Re = 6400 with the Navier-Stokes equations. We consider coarse mesh resolutions that are representative of the type required for practical high-Reynolds number simulations of turbulent flows with complex geometries. For these types of studies, computational cost necessitates coarse meshes relative to the Kolmogorov microscale. We also do not use any anti-aliasing strategies in the current study, as these significantly increase the computational cost for industrial-scale simulations. This approach also allows us to investigate the ability of E-ESFR schemes to implicitly suppress aliasing driven instabilities, as was shown previously for one-dimensional problems by Vincent et al. [10]. The Taylor-Green vortex is particularly appealing because it initially contains only large-scale laminar structures, which then undergo transition to fully turbulent flow at the later stages  of the simulation [19][20][21][22]. This means that the initial flow field is well resolved, even with coarse meshes. However, as the simulation advances, it is expected to become progressively under-resolved as more energy is transferred to the smallest length scales. It can be expected that more stable numerical schemes will be able to simulate the Taylor-Green vortex test case for a greater amount of time, which corresponds to a wider range of turbulent length scales. Therefore, we consider total achievable simulation time in this test case as a proxy for the stability of a particular scheme for such under-resolved simulations, such as ILES.
where u, v, and w are the velocity components, p is the pressure, and ρ is the density. The terms T 0 and U 0 are constants specified such that the flow Mach number is Ma = 0.1 based on U 0 , effectively incompressible. The domain is a periodic cube with the dimensions −π L ≤ x, y, z ≤ +π L. These parameters yield a reference time scale of t c = L/U 0 . Both sets of simulations were run with approximately 32 3 degrees of freedom, using 8 3 and 6 3 elements for k = 3 and k = 4, respectively. It is well known that aliasing errors can cause non-linear instabilities for under-resolved ILES simulations. Therefore, schemes that can implicitly suppress aliasing driven instabilities are of particular interest, and will be identified in this study.
The study was performed by systematically varying q 0 and q 1 such that q 0 ∈ [−0.4, 1.0] and q 1 ∈ [−0.4, 0.6] in increments of 0.01. An independent simulation was run for each q 0 and q 1 pair, resulting in a total of 14, 241 simulations for each order, and each was run until the solution diverged. All simulations were run using a modified version of PyFR 1.1.0 [23] (modifications supplied as a patch in Electronic Supplementary Material). The maximum achieved simulation time was then recorded for each data point. All simulations were performed using an adaptive five-stage fourth-order Runge-Kutta (RK45) scheme with absolute and relative error tolerances of 10 −6 . Contours of the maximum achievable simulation time for the inviscid cases are shown in Fig. 12 for both schemes. It is clear from both of these plots that the stability of a particular scheme for such under-resolved simulations is highly-dependent on the choice of q 0 and q 1 . In general, schemes near the energy stability boundaries perform poorly, and typically advance the solution less than one t c for both k = 3 and k = 4 schemes. However, both plots show that particular subsets of schemes perform significantly better. For k = 3 there is a wide range of schemes, primarily with q 0 ∈ [−0.2, 0.6] and q 1 ∈ [−0.1, 0.2] that are significantly more stable. For k = 4 there is a much smaller subset of schemes with q 0 ∈ [−0.1, 0.8] and q 1 ∈ [−0.1, 0.1] that are found to be particularly stable.
Another observation is that the FRDG methods, although common and widely used, were not found to be the most stable for this test case.
All of the Taylor-Green vortex simulations were then repeated for the Re = 6400 case using the Navier-Stokes equations with k = 3 and k = 4. Contours of the maximum achievable simulation time for these cases are shown in Fig. 13 for both polynomial degrees. Although the Re = 6400 simulations are typically more stable than the inviscid simulations, due to the addition of physical viscosity, their overall behaviour is generally similar. Schemes near the boundaries of the energy-stable region tend to perform quite poorly, advancing the simulation often by less than one t c for both k = 3 and k = 4. However, there are particular subsets of schemes that perform significantly better and, importantly, these tend to be the same schemes that performed well in the inviscid cases.
Plots of the dissipation and dispersion profiles for the most stable k = 3 and k = 4 E-ESFR schemes from the inviscid Taylor-Green vortex simulation are shown in Fig. 14 and Fig. 15, respectively. These are shown alongside the FRDG schemes of the same polynomial degree for reference. These E-ESFR schemes behave similarly to their reference FRDG scheme, in the context of the range of possible schemes presented previously. However, they have increased numerical dissipation at low wavenumbers, while less dissipation error at high wavenumbers. Also, they have reduced group velocities for the majority of wavenumbers. Importantly, they have significantly less dispersion error for all but the highest resolved wavenumbers, when compared to their corresponding reference FRDG schemes. These results demonstrate that modifying the correction function could be a possible route to stabilising under-resolved ILES simulations. In addition, several of the E-ESFR schemes outperform the commonly used FRDG scheme, at least in terms of stability, for this test case and exhibit consistent dispersion and dissipation properties. Interestingly, these most stable schemes also had significantly larger τ C F L values when compared to the FRDG scheme, as shown previously in Fig. 2 for RK44. For example, the k = 3 scheme with q 0 = 0.3 and q 1 = 0.02 had a 62.6% greater τ C F L and the k = 4 scheme with q 0 = 0.05 and q 1 = 0.01 had a 12.8% greater τ C F L , when compared to corresponding FRDG schemes of the same polynomial degree.

Turbulent flow over an SD7003 aerofoil
To further assess the behaviour of E-ESFR schemes for ILES of the Navier-Stokes equations, we investigate transitional and turbulent flow over an SD7003 aerofoil [24] using the k = 3, q 0 = 0.30, and q 1 = 0.02 scheme and the k = 4, q 0 = 0.05, and q 1 = 0.01 scheme. These are the most stable schemes for p = 3 and p = 4, respectively, identified in the previous inviscid Taylor-Green vortex study. Each simulation is run using a modified version of PyFR 1.1.0 [23] (modifications supplied as a patch in Electronic Supplementary Material) at a Reynolds number Re = 60,000, Mach number Ma = 0.2, and angle of attack α = 8 • . This test case is commonly used to examine the suitability of numerical schemes for predicting separation, transition, and turbulent flow. It has been studied previously by, for example, Visbal and collaborators including Visbal et al. [25] and Garmann et al. [26], and Beck et al. [27] using a DG spectral element method (DGSEM). The characteristic features of the flow include laminar separation on the upper surface of the aerofoil, which then reattaches further downstream forming a laminar separation bubble. The flow transitions to turbulence partway along this separation bubble, creating a turbulent wake behind the aerofoil.
We use an unstructured hexahedral mesh with quadratically curved boundaries to match the aerofoil geometry, as shown in Fig. 16. The domain extends to 10c above and below the aerofoil, 20c downstream, and 0.2c in the span-wise direction, where c is the aerofoil chord length. A structured mesh is used in the boundary layer region, with a fully unstructured and refined wake region behind the aerofoil to capture the turbulent wake. The boundary layer resolution gives y + ≈ 0.6 and y + ≈ 0.4 at the first solution point off the surface for k = 3 and k = 4, respectively, where y + = u τ y/ν, u τ = C f /2U ∞ , U ∞ is the free-stream velocity magnitude, and C f ≈ 8.5 × 10 −3 is the maximum skin friction coefficient in the turbulent region  reported by Garmann et al. [26]. The mesh has a total of 138,024 hexahedral elements with 12 elements in the span-wise direction. An adaptive RK45 temporal scheme was used with relative and absolute error tolerances of 10 −6 [28][29][30], LDG [31] and Rusanov type [23] interface fluxes, and Gauss points for both solution and flux points in each element. The ratio of specific heats is γ = 1.4, the Prandtl number is Pr = 0.72, and constant viscosity is used due to the relatively low Mach number. The simulation is initially run to t = 20t c , where t c = c/U ∞ , to allow the flow to develop, separate, and transition.
Statistics are then extracted over an additional 20t c , including span-wise averaging where appropriate.
Instantaneous isosurfaces of q-criterion coloured by u/U ∞ , where u is stream-wise velocity magnitude, are shown in Fig. 17 for both simulations. There is a relatively short laminar separation bubble that forms on the upper surface of the aerofoil. The flow then undergoes transition near the 1/4 chord location, creating a turbulent wake extending downstream behind the aerofoil. This behaviour is consistent with observations from previous studies, such as the similar results presented by Garmann et al. [26] and Beck et al. [27]. Plots of u/U ∞ , where u is the time and span-averaged stream-wise   velocity magnitude, are shown in Fig. 18 for both cases. A relatively short separation bubble on the upper surface is clearly visible in both sets of results. This separation bubble is slightly longer in the k = 4 simulation, and is generally consistent with previous studies in terms of separation point (x sep ) and reattachment point (x rea ) as shown in Table 2. Table 2 also shows the time-averaged lift and drag coefficients for the current study. Again, these are in agreement with values reported in previous studies. A plot of the time and span-averaged surface pressure coefficients is provided in Fig. 19, alongside the results of Garmann et al. [26] using a finite-difference scheme and Beck et al. [27] using a DGSEM scheme at two polyno- Table 2 Results from the current SD7003 test cases using k = 3 with q 0 = 0.30, q 1 = 0.02 and k = 4 with q 0 = 0.05, q 1 = 0.01, and reference datasets for comparison.  mial degrees. The current k = 3 simulation under-predicts the pressure in the separation bubble relative to the reference datasets. However, the transition point is still in excellent agreement with the results of Garmann et al. [26], and the pressure coefficient in the fully turbulent aft portion of the aerofoil agrees well with both reference data sets. The higher-order k = 4 results show excellent agreement with both reference datasets throughout the pressure curve. In particular, when compared to the k = 7 DGSEM results of Beck et al. [27], the current k = 4 simulation shows excellent agreement even in the transition region of the separation bubble at x/c ≈ 0.3.

Conclusions
We have investigated the behaviour of the E-ESFR range of energy-stable symmetric conservative FR correction functions introduced by Vincent et al. [16]. Insights from von Neumann analysis identified values of q 0 and q 1 that yield the largest stability limit when using the explicit RK44 scheme. Additionally, it was shown that as q 0 → ∞ the two-parameter k = 3 and k = 4 schemes recovered the O-ESFR set of correction functions of degree k − 1. Schemes with the highest orders of accuracy based on dispersion and dissipation analysis were found on the line q 1 = 0, with a global peak corresponding to the FRDG scheme with q 0 = 0 and q 1 = 0 with A T = 2k + 1.
Further analysis demonstrated that both q 0 and q 1 have a significant effect on the dispersion and dissipation behaviour of the two-parameter schemes. While increasing q 0 decreases the dispersion relation for all wavenumbers, increasing q 1 tends to decrease the dispersion relation at low wavenumbers and, conversely, increases the dispersion relation for high wavenumbers. Increasing q 0 also decreases the group velocity of well-resolved waves, while increasing q 1 tends to decrease the group velocity for well-resolved waves, and increase it for marginally-resolved waves. Similar to decreasing q 0 , increasing q 1 tends to increase the amount of numerical dissipation at high wavenumbers and decreases dissipation at low wavenumbers. This shows that the dispersion and dissipation characteristics of E-ESFR schemes, which are important for the Fig. 19. Pressure coefficient C P as a function of x/c for current SD7003 simulations using k = 3 with q 0 = 0.30, q 1 = 0.02 (top) and k = 4 with q 0 = 0.05, q 1 = 0.01 (bottom), and reference results from Garmann et al. [26] and Beck et al. [27].
suppression of aliasing-driven instabilities when using ILES [9,10], can be controlled via judicious choice of the correction function.
Subsequently, simulations of the Taylor-Green vortex showed that the choice of correction function can significantly influence the stability of under-resolved simulations that are affected by aliasing driven instabilities. Bands of particularly stable schemes were identified, primarily with small positive q 0 and q 1 for both k = 3 and k = 4. Importantly, several schemes with non-zero q 0 and q 1 were found to be more stable than the commonly used FRDG scheme for this test case. Analysis of the dissipation and dispersion characteristics of these E-ESFR schemes showed that they have more dissipation than the FRDG scheme for well-resolved wavenumbers and less dissipation at high wavenumbers. Importantly, they have significantly less dispersion error for the majority of wavenumbers. These results demonstrate that the choice of correction function has an important influence on the stability of FR schemes for ILES. Finally, favourable k = 3 and k = 4 schemes identified from the Taylor Green vortex simulations with q 0 = 0.30 and q 1 = 0.02 and q 0 = 0.05 and q 1 = 0.01, respectively, were used to simulate turbulent flow over an SD7003 aerofoil, achieving good agreement with previous studies of this test case [26,27]. These results demonstrate that E-ESFR schemes can be more stable than FRDG in practice for ILES.
To put these results into context, high-order ILES simulations using the FR and DG methods are typically performed with de-aliasing to improve stability [32]. This requires non-linear flux functions to be represented in a higher polynomial degree space, and then projected down onto the same space as the solution. While de-aliasing has been shown to improve stability, it also incurs significant additional computational cost in terms of both memory and number of operations per degree of freedom [33]. In contrast, we have shown that by simply changing the FR correction function we can observe improved stability relative to FRDG, while simultaneously increasing the CFL limit. This appears to be a simple approach for stabilising ILES simulations using the FR approach, and could also be used in conjunction with de-aliasing if additional stabilisation is required. Future work will focus on the development of stability proofs for the E-ESFR discretisations of the Navier-Stokes equations, and on understanding the impact of correction function choice on the accuracy of ILES simulations across a range of canonical test cases.

Data statement
Data Statement: Data relating to the results in this manuscript can be downloaded as Electronic Supplementary Material under a CC-BY-NC-ND 4.0 license.