Effect of Mesh Quality on Flux Reconstruction in Multi-Dimensions

Theoretical methods are developed to understand the effect of non-uniform grids on Flux Reconstruction (FR) in multi-dimensions. The analysis reveals that the same effect of expanding and contracting grids is seen in two dimensions as in one dimension. Namely, that expansions cause instability and contractions cause excess dissipation. Subsequent numerical experiments on the Taylor-Green Vortex with jittered elements show the effect of localised regions of expansion and contraction, with an initial increase in the kinetic energy observed on non-uniform meshes. Some comparison is made between second-order FR and second-order finite volume (FV). FR is found to be more resilient to mesh deformation, however, FV is found to be more resolved when operated at second order on the same mesh. In both cases, it is recommended that a kinetic energy preserving/conservation formulation should be used as this can greatly increase resilience to mesh deformation.


I. Introduction
Since the inception of Spectral Volume methods by Wang, 1 the trajectory of high order methods has trended towards the Flux Reconstruction method of Huynh2 and Vincent et al. 3 This approach idiomatically draws on the work of Finite Elements, see Brenner and Ridgway-Scott, 4 enabling the high performance of Flux Reconstruction on heterogeneous computing -as can be seen in the highly efficient use of vast computing resources by Vincent et al. 5 However, the move towards high order was not born out of a need for more efficient use of modern HPC environments.For example, Brandvik and Pullan 6 showed that high throughput could be obtained using second order Finite Volume (FV) methods.Instead, the main motivating factor has been the increased uptake by industry of turbulence resolving methods, such as Large Eddy Simulation (LES), as this allows for far better exploration of flow physics and moves towards the aim of computational wind tunnels.The main feature of LES is the modelling of the very smallest scales of motion, which removes the need for DNS levels of resolution.However, Chow and Moin 7 and Ghosal 8 showed that, for LES, the need to keep the truncation error small to enable the sensible use of sub-grid scale models meant that the grid requirements were demanding.A move to higher order would mean that the scaling of the truncation error was to a higher power of grid spacing -thus lowering the grid requirements and decoupling the scaling of aliasing error and truncation error.Hence, for wall resolved LES, calculations are often impractical unless the improved mesh resolution requirements of high order methods are considered.
The analytical understanding of Flux Reconstruction has been explored to large extent in the work of Vincent et al., 9 Jameson et al., 10 and Castonguay et al., 11 where the stability of linear advection, advectiondiffusion, and non-linear problems was presented.The key findings were the energy stability of FR on linear problems, and the condition for energy stability on non-linear problems.In addition, by investigating the dispersion and dissipation characteristics of FR, the existence of superconvergence after temporal integration and the corresponding CFL limits were found.This work was limited to one dimension -although still applicable, the investigation of the exact behaviour of FR in higher dimensions has been limited, such as that of Williams and Jameson 12 and Sheshadri and Jameson. 13This work focused primarily on the proof of the Sobolev type energy stability in 2D in a similar manner to that of Hesthaven and Warburton, 14 alongside some numerical studies performed for validation.
The advantage of FR -one that leads to high performance on heterogeneous and massively parallel architectures -is its unstructured and sub-domain nature.Unstructured grids also allow far more complex geometries to be considered, but the resulting meshes experience deformation, expansion, and contraction of the elements.We wish to characterise the performance of FR under these conditions, and so far the effect of linear mesh deformation on FR has only been considered in one dimension by Trojak et al. 15 Therefore, we make use of the seminal work of Lele,16 in which the dispersion and dissipation of finite difference methods were considered in both one and two dimensions.We wish to repeat this process for FR, but extend it to also consider deformed grids.
In this paper, we present an extension to the one dimensional analytical work of Vincent et al. 9 and Trojak et al. 15 This extension will be shown for a two dimensional case on quadrilaterals with rectilinear mesh stretching, but could also be performed on higher dimensional hypercubes.From the basis of this more general von Neumann analysis, the behaviour of FR on linearly mapped meshes can be explored.The aim of this work is to understand the effect of moving to higher dimensionality on key metrics governing scheme performance, such as CFL limit, dispersion, and dissipation.Finally the Taylor-Green Vortex will be used to understand how deformed meshes affect full Navier-Stokes calculations, with reference calculations performed by an industrial second order finite volume method.

II. Flux Reconstruction
Flux Reconstruction 2, 17 (FR) applied to the linear advection equation will form the basis of the initial investigation to be carried out, and for the reader's convenience an overview of the scheme is presented here.However, for a more detailed understanding the reader should consult Castonguay 17 or Huynh. 2 This 1D scheme can be readily converted to two and three dimensions for quadrilaterals and hexahedrals, respectively.First, let us consider the one dimensional advection equation: The FR method is related to the Discontinuous Galerkin (DG) method 18 and makes use of the same subdivision of the domain into discontinuous sub-domains: Within the standardised sub-domain, Ω ∈ R d , computational spatial variables are defined.When d = 1, Ω = [−1, 1], using ξ to denote the value taken.This computational space is discretised with (p + 1) d solution points, and 2d(p + 1) d−1 flux points, placed at the edges of the sub-domain.The solution and flux point locations are determined using a tensor grid of a 1D quadrature.Fig. 1a shows a 1D example of this.To transform from Ω n → Ω, a Jacobian J n is defined such that: With this domain set up, we now proceed with defining the steps to construct a continuous solution from the discontinuous segments.The first stage is to define a local solution polynomial in Ω using Lagrange interpolation.
Repeating the interpolation for the discontinuous flux in Ω: Now, using the Jacobian and the solution polynomials, the primitive and flux values can be calculated in the physical domain Ω n : Defining ûδ l = ûδ (−1) and ûδ r = ûδ (1), this can be repeated for the flux values.Once the primitive values at the interface, I, have been interpolated, a common interface flux, f c I , can be calculated in the physical domain.For a general case this is done using a Riemann solver on the primitives at the interface, such as: Roe; 19 flux vector splitting; 20 or HLL. 21In order to get a spatially continuous solution over Ω the common interface flux has to be incorporated into solution.For FR this is done by using a correction function to propagate the corrected flux into Ω n .Figure 1a shows the g 2 correction function proposed by Huynh. 2 Generally the left and right correction functions are defined as h l (ξ) and h r (ξ).
The procedure to correct the discontinuous flux in the sub-domain, Ω, is: Finally, the solution is advanced in time following Eq.( 9) -this can be done via a sensible choice of temporal integration.
The definition of correction functions is an evolving field, however in the present study we are concerned with the effect of grid deformation on FR.Therefore, we will restrict ourselves to the simplest one parameter family of correction functions: original stable FR (OSFR).

III. Two-Dimensional Von Neumann Analysis
The procedure for investigating the dispersion and dissipation properties of finite element methods has been laid out in some detail by Huynh, 2 Hesthaven and Warburton, 14 and Vincent et al. 9 and is broadly classified as a von Neumann analysis.The procedure was, however, only performed in 1D, with critical insight into the analytical performance of FR when applied to more realistic problems overlooked.Extension of the analysis to higher dimension domains was performed by Lele 16 for various finite difference schemes.This did, however, avoid the increased complexity of finite element von Neumann analysis.To begin our extension we introduce the 2D linear advection equation: Flux reconstruction then uses the superposition of the discontinuous and corrected flux divergence, meaning Eq. ( 11) can be rewritten as: Taking the following definition of the Jacobian, the computational-physical domain transformation can be defined: where we use ∇ to mean [ ∂ ∂ξ , ∂ ∂η ] T in 2D.From the work of Huynh, 2 Castonguay, 17 and Sheshadri et al., 22 Eq.( 8) can be written in two dimensions in the following way: where we use l, r, b, and t subscripts to mean left, right, bottom, and top respectively.We have also ensured that grid transformations are purely rectilinear, i.e.G 2 = G 3 = 0. We may now use Eq. ( 16) and convert it into a matrix form: To apply the correction function, we need to calculate the the interface values around the element.For the case of generalised central/upwinding with upwinding ratio α, the common interface fluxes may be written as: where α = 1 gives rise to upwinding and α = 0.5 produces central difference.Hence the divergence correction can be written as: where h l is the gradient of the left correction function at the solution points and l l are the values of the polynomial basis at the left interface and so on for r, t, and b.Therefore, by grouping terms by their cell indexing and transforming each term into the physical domain: where: Finally, we are here interested in the frequency response of the system, and, importantly to engineers and technicians, how the cell's orientation relative to an oncoming wave affects performance.Therefore, we impose a trial solution of the form: and by substitution into Eq.(11), the advection velocity, a, can be found, which is shown schematically in Fig. 2a.
The plane wave can then be projected into the computational domain and discretised as: where, for brevity, δ i = x i −x i−1 and δ j = y j −y j−1 are defined.Inserting Eq.( 29) into Eq.(24), an Eigenvalue problem can be obtained as: where (c(k)) and (c(k)) are equal to the dispersion and dissipation, respectively, for FR.Alternatively, this can be cast in the form of an update equation.If initially Eq.( 24) is combined with Eq.( 29), then a new matrix, Q i,j , can be defined: The update equation takes the linear FR operator Q and imposes some temporal discretisation, so we may write: where the superscript denotes the time level, and our update matrix is R. Shown here is also an example definition for R for a 3-step 3 rd -order Runge-Kutta time integration scheme.Finally, in keeping with von Neumann's theorems 23,24 and Banach's fixed point theorem, 25 the spectral radius of R has to be less than or equal to 1 for stability.ρ(R) 1 ∀ k ∈ R. Extending this analysis further, the linear FR operator matrix, Q i,j , can be diagonalised as: where W is a matrix of eigenvectors and Λ is a diagonal matrix of normalised eigenvalues.From this the weights of modes used to reconstruct the solution can be found as: where β β β is an array of mode weights used to project u i,j into the functional space of FR.The posedness of the projection can then be measured for wavenumbers using the condiontion number of the matrix of modes defined as: where λ(W i,j ) is an eigenvalue of W i,j , with the matrix approaching singular as κ → ∞.
The results of this section can then be extended to n-dimensions and, in particular, the three dimensional case will be investigated to show the continuation of trends with higher dimensionality.The analysis can broadly be repeated and is excluded for brevity, but importantly the prescribed solution is taken as: where the angles are as shown in Fig. 2b, and hence the 3D convective velocities for linear advection are:

IV. Analytical Findings
The analytical methods presented in Section III allow us to investigate many properties of FR, however from Eq. (32-34) it can be seen that the functional space of Q is 8 dimensional, leading to the functional space of ρ(R) being 9 dimensional (τ, γ x , γ y , ∆ x , ∆ y , k, θ, ι, p).Therefore we need to restrict our investigation to some key results relating to grid deformation.Firstly, understanding the dispersion and dissipation ( (ω) & (ω)) in 2D for both uniform and stretched grids will be important.Secondly, we wish to understand how higher dimensionality and grid deformations affect the temporal stability of FR through evaluation of the CFL limits. 26Here the dispersion and dissipation relations will be useful in explaining the trends seen, and will aid in linking this work to that of Trojak et al. 15 The definition of the CFL limit in higher dimensions will be taken as: where d is the dimensionality and a i and ∆ i are the characteristic velocity and grid spacing in the i th dimension, respectively.Finally, we wish to understand if correction functions can be used to alleviate any effects of deformation by understanding how the scheme properties vary with correction function.Within this study the link to the quality of posing of the linear system will be explored alongside how this relates to the other properties.

IV.A. Effect of Grid Expansion on Dispersion and Dissipation
We begin by considering the dispersion and dissipation on a uniform grid in two dimensions.We are concerned here with the primary mode -as FR has multiple modes, this is the one that physically represents the wave, although, as was found by Asthana, 27 this may not be how the energy distributes itself.The physical mode is identified as the dispersion relation that goes through zero and had a dissipation that is also zero and remains constant for small wavenumbers.The results of this are shown in Figs. 3 and 4. It is clearly seen that for all orders FR becomes more dispersive and dissipative at θ = 45 • .Furthermore, there doesn't seem to be any widening in the range of angles over which FR becomes more dispersive and dissipative as the order in reduced.By comparison with the the results of Lele, 16 where a similar test is performed for standard and compact difference schemes, FR shows a comparatively smaller change in performance as angle is varied.It is thought that this due to the method of polynomial fitting used by FR, namely that this implementation of FR used a tensor grid of monomials i.e., the number of solution points is (p + 1) d and hence the monomials in the interpolation go from (ξ 0 η 0 , ξ 1 η 0 . . .ξ p η p ).By contrast, finite differences do not include the cross product terms, which will become increasingly dominant as the angle is increased.
Moving on, we then consider the impact of non-uniform grids on the character of the dispersion and dissipation.There are two cases that have been identified from previous work as being of interest.Firstly, when the grid is expanding, does this give the same positive dissipation in higher dimensions as seen in the 1D case?Secondly, if positive dissipation is seen in the 2D case, will an orthogonal contraction help to stabilise the grid?
The first of these questions is explored in Figs.5a and 5b.It can be observed in Fig. 5b that expanding grids do cause positive dissipation in higher dimensions.It was previously explained in one dimension that as a wave moves through elements of different size the group velocity (dω/dk) will change.For an expanding grid this leads to low wavenumber energy collecting in elements.The same mechanism looks to be responsible in higher dimensions.The next question of whether a contraction orthogonal to the expansion will help to stabilise the situation is the considered in Figs.5c and 5d.From Fig. 5d, it can be seen that the answer is yes, however as the incidence angle of the wave approaches θ = 0 the stabilistion brought about by the contraction decays.The θ = 0 case is identical to the case of γ y = 1.

IV.B. Effect of Grid Aspect Ratio on CFL Limit
Beyond the resolution of the scheme is the question of setting up a case and running it on some machine.For this, knowledge of the temporal stability is key, and we will begin by looking at the effect of the relative size of an element in x and y on the CFL limits of FR.In this case the grid is not expanding or contracting, merely ratio of ∆ x to ∆ y is varied.
As can be seen from Fig. 6 there is a clear impact on the CFL limit of FR when elements are rectangular, with no change in the CFL limit when the waves are aligned with the grid.What is evident is that the angle of incidence where the CFL limit is smallest, for a given size ratio, is when a wave is incident at an angle of tan −1 (∆ y /∆ x ).This angle corresponds to the maximum length across the element and hence, when the wave is decomposed into x and y components, the wavenumbers are lowest -i.e., inf θ∈R (sup {∆ y cos (θ), ∆ x sin (θ)}) when θ = tan −1 (∆ y /∆ x ).As the wavenumbers will be at their lowest necessary to form the wave, from the 1D dissipation of FR, the dissipation will also be at its lowest.Therefore, there is less dissipation in the spatial scheme available to counteract the negative dissipation of the temporal scheme, hence reducing temporal stability at θ = tan −1 (∆ y /∆ x ).

IV.C. Effect of Grid Expansion on CFL Limit
Now we introduce to the grid an expansion or contraction in x and y, with varying incident angles.The results of these distortions are shown in Fig. 7. What can be seen broadly from Fig. 7 is that the minimum CFL limit is at θ = 45 • , with temporal performance peaking as expected at θ = 0 • , 90 • .This result agrees with that of Fig. 6.However, this also shows that the angle of minimum CFL is only dependent on the local element shape, and in the case investigated here the central element is always square.Furthermore, when the CFL limit in the quasi-1D case (θ = 0 • , 90 • ) is compared to the results of Trojak et al., 15 the CFL limit is found to be lower than the 1D case.This may be due in part to the increased modes of the system and their coupling leading to a less stable system.This decrease is corroborated by numerical tests, and the analytical reasoning will follow shortly.
A second point to note is that for non grid-aligned waves the expansion or contraction in both components affects stability, with a contraction orthogonal to an expansion again helping to stabilise the scheme.This point is subtle, because if the decomposition of the wave into x and y were linearly independent then the lowest would dominate.However, this result demonstrates that there is coupling between the x and y components -which could be used advantageously, as was discussed earlier.
Throughout the analytical tests in which waves were injected at incidence on a square central element, the Nyquist wavenumber was found to be k nq,θ = k nq,0 / cos(θ) for 0 • θ 45 • .This result can be understood in one of two ways.Firstly that a wave at an angle can draw on more points in the normal direction to form a fit of higher wavenumbers.Or that the wave can be thought of as being a coupled decomposition of the wave into the x and y directions, and although it has been said that these are not independent, this does mean that higher wavenumbers can be supported.To understand further why expanding meshes are less stable and contracting meshes more stable it can be illuminating to consider the 1D linear advection equation, which for upwinded FR can be written as: using Euler's method to temporally integrate this equation, we find: where the superscript n denotes the time step.If u n−m j−m−1 (∀ m ∈ N) is then recursively substituted, the final form is then: where we assume the solution is on a geometrically expanding grid in order to substitute for the Jacobian -hence being only valid for linearly transformed elements.If we consider that C 0 and C −1 are linear operators, then rather than prescribing a solution, the dynamics of linear operators can be used. 28So, if the mesh extends infinitely downwind, then it is sufficient to say that Eq.( 43) is stable when it is a hypercylic orbit.Hence, the stability criterion is that −2γτ C −1 = T is a matrixable linear hypercylic operator.The definition of which is that sup T n 1, ∀ n ∈ N, which in turn implies that ρ(T) = 1.What this aims to show is that the stability criterion is dependent on the product of τ and γ, as well as on C 0 .Therefore, as γ increases the maximum stable τ decreases, for constant C −1 .This also explains that although the scheme may be formally unstable, correct setting of τ for a given γ can lead to a T that is still hypercyclic and give a bounded solution.However because of the (−1) m this does mean that if U n−1 ⊂ R is the set of solutions then at some time step n − 1 the solution u n j U n−1 i.e, for a sinusoidal solution the computed valued may exceed the prescribed magnitude.

IV.D. Effect of Correction Functions with Grid Expansion
As was mentioned in Section I, a series of correction functions with peak temporal stability and spatial accuracy was proposed by Vincent et al., 9 defined by a correction parameter, ι + .These correction functions exhibited the superconvergence expected of Nodal DG (Cockburn et al. 29 ), however ι + = ι DG as they account for the variation in the dispersion caused by the discrete temporal integration.We wish to investigate whether the advantage of this family of correction functions is maintained in 2D or if an analogue can be found.Therefore, the correction function parameter is varied for different angles and grid expansion rate, the results of which are presented in Fig. 8. Using Fig. 8c as an example, the peak CFL at ι + can be seen clearly in the case of θ = 0 • , with its peak value reduced in comparison to the 1D case.(This was discussed in Section IV.C.However, the clear peak at θ = 0 • , 90 • does not significantly persist as the wave angle increases to θ = 45 • , with the peak becoming substantially flattened.Therefore, the balancing effect that the modification of correction function has on the dispersion of the scheme seems to have limited scope.There is also no other correction function that seems able to achieve the same effect at the intermediate range of angles.Furthermore, the persistence of the peak CFL limit is not seen as the expansion rate and order is varied, with ι + in fact suffering the most appreciable decay in performance as angle is varied, when compared to the other correction functions.
To investigate further the impact of varying correction functions on higher dimensional problems, we will consider the projection of the solution into the functional space of FR for linear advection.The method for understanding this is by studying how well posed the linear operators are.The process of projection was outlined in Section III.Presented here are the several results showing how the posedness varies with angle, correction function, and order.
Several insights into the different multidimensional behaviour of FR can be gained by studying Fig. 9.

By comparison of Fig 9a & 9b
it can be seen that Huynh's g 2 correction function causes the projection to be more ill-posed on average compared to Nodal DG methods.Hence, the superconvergent DG recovering scheme of Vincent et al. 9 has decreased temporal stability.This is because the ill-posedness indicates the sensitivity of the reconstruction to change, a more sensitive reconstruction means that error can result in energy being transferred to more dissipative modes.Furthermore, in both cases, across all wavenumbers, the condition number can be seen to increase with incident angle.When compared with Fig. 7 and Fig. 8c  it can be seen that for a given correction function the CFL limit reaches its minimum at θ = 45 • , therefore a stark increase in condition number can also cause a decrease in temporal stability, potentially due to too much transport between modes.A result presented by Trojak et al. 15 was that for FR the best points per wavelength (PPW) performance was seen for p = 4 with the PPW increasing for orders higher than this.A result that is exhibited by Fig .9c is that the the condition number for p = 5 schemes is higher than for lower orders and may have passed a point where increased order is out weighted by inaccuracy in ill-conditioning.Which may explain the optimal result seen, and is touching on the fundamental problem characterised by Runge's phenomenathat high-order may introduce accuracy through high-order but may also introduce inaccuracy through poor conditioning.A result that cannot be clearly seen in either Fig 9a or 9b however, was exhibited by FR.For θ = 0 the condition number was significantly higher than that found for 1D FR owing to the natural poor conditioning of a 2D system acting as a quasi-1D one.This, therefore, explains the lower CFL number for this case, as was shown in Fig. 7. To show that the results found in this section extend to higher dimensions, the von Neumann analysis was repeated for 3D 'hexi-linear' grids.The primary result of interest is the illposedness of the functional projection and this can be seen in Fig. 10.The message is that, as expected, the ill-posedness of the reconstruction increases with order, and, while small increases in the condition number can give increased temporal stability, larger increases in condition number tend to act to destabilise the coupled spatial-temporal scheme.

V. Non-linear Navier-Stokes Equations with Randomised Grids
In the preceding sections, the focus has been on analysis of the linear advection equation.However, this has limited scope when applied to the equation sets that are considered practically.Therefore, we wish to understand what the effect of grid deformation is when we move to a full non-linear implementations of the Navier-Stokes equations, and how the behaviour observed for linear advection carries across.
It is very common within the CFD community to use the canonical Taylor-Green Vortex (TGV) 30 test   case to assess the numerics of a solver applied to the Navier-Stokes equations with turbulence -and to that end there is a plethora of DNS data available for comparison. 31,32 owever, this case is quite contrived and ultimately will favour spectral or structured methods due to the Cartesian and periodic domain, whilst also being unrepresentative of engineering flows that are often wall bounded and/or have complex geometries.Hence, we propose linearly deforming the elements of the mesh by jittering the corner nodes to be be more representative of real mesh conditions.Importantly, these deformation will introduce cross multiplication into the Jacobian, as well as local regions of expansion and contraction.
The initial conditions of the TGV being used here are those of DeBonis, 32 where the character of the flow is controlled by the non-dimensional parameters defined as: where we will use the standard set of free-variables for the velocity, density, pressure, and gas characteristics: Here, due to the solver implementation, we use a specific gas constant of unity and hence, to achieve the required Reynolds and Prandtl numbers, the dynamic viscosity and thermal conductivity can be set appropriately.
As has been stated, we will take the uniform periodic mesh on the domain Ω ∈ [−π, π] 3 , and jitter the corner nodes of the elements that are interior to the domain.The amount of jitter is calculated using a time seeded random number shifted to be centred about zero and scaled by a global factor between zero and unity.The scaling factor is such that zero gives a uniform mesh and unity could lead to edges of zero length.After jittering, the solution points are then linearly positioned within the element using the thin plate spline radial basis function together with the mapping from the uniform to jittered corner nodes.This gives mapping of solution points to solution points within the jittered elements.Finally, a quality metric is needed to describe, in a single number, the relative quality of the meshes produced.We opted for a volume ratio shape factor, slightly redefined as: where S h is the surface area of the hexahedral element and V h is the volume of the hexahedral elements.The quality metric, q h , is then defined as the ratio of the volume of the element to the volume of a sphere with the same surface area, with q h = π/6 for a perfect cube.To put this parameter into context, some example meshes are shown in Fig. 11.
(a) q h = π/6 ≈ 0.7236 (b) q h = 0.7201 (c) q h = 0.7016 The statistics that will be studied here are the decay of the kinetic energy and the enstrophy decay rate, which are defined respectively as: where ω ω ω = ∇ × [u, v, w] T is vorticity and |Ω| is the domain volume.q h = 0:0 q h = 0:7157 ref  q h = 0:0 q h = 0:7000 ref  Fig. 13 shows the first of these results.It may be observed that as the grid quality decreases there appears a period of time near the start of the simulation where the turbulent kinetic energy increases.As time progresses, dissipation is again seen and the point of peak dissipation arrives early, moving from t ≈ 8.5 to t ≈ 7.5.
The explanation of this is believed to be that initially the regions in the mesh that are locally expanding cause an increase in the energy due to the positive dissipation.It was discussed in Section IV.A, that for the linear advection equation, dissipation is positive at low wavenumbers for expanding grids and negative for contracting grids at the same wavenumber.It is thought that as the simulation progresses, the energy cascade of large scales to small scales then means that more of the solution lies in the more dissipative higher wavenumber region for both expanding and contracting grids.Therefore, the net dissipation at a later time is higher than the uniform case and hence the peak dissipation is earlier.This is consistent with a lower R e and hence higher global dissipation.We can conclude that the stability of this case is brought about by the physics of the Navier-Stokes and cascade of energy from low to high wavenumbers, which side steps the problem of positive dissipation of low wavenumbers on expanding meshes.This result is interesting as it is in slight contradiction to the result of Trojak et al., 15 where FR was found to be more robust to grid for Euler's equations.However, in that investigation the Isentropic Convecting Vortex was considered where there was a large convective velocity.This may expose slightly different properties, although this is not yet fully understood.
To provide some reference as to how FR performs relative to an established method we will use an edge-based Finite Volume (FV) method for comparison.The FV method is a standard central second order method with L2 Roe smoothing 33 for stabilisation, which has been validated previously. 34The particular FR scheme used in this comparison is p = 1, giving second order, the same as the FV scheme.However, this puts FR at a significant disadvantage as its numeric characteristics at low order are particularly poor.For example, consider the dispersion and dissipation relations in Fig. 14, which, by comparison to the result of Lele,16 show that FR has noticeably lower resolving abilities when compared against a second order FD scheme.

0
:=4 :=2 3:=4 :  With this in mind, we present the results of tests on various jittered grids with a total of 170 3 degrees of freedom in Fig. 15.For the uniform case the enstrophy clearly shows that FR is underresolved compared to FV, which is also shown by a slightly increased rate of dissipation earlier -indicating that the implicit filter is too narrow.If we now consider the effect of jittering, several things may be concluded.Firstly, we were unable to run FR with j f = 0.5 as the simulation quickly became unstable.Secondly, for −dE k /dt it seems that the peak value is less sensitive with FR than with FV, with central FV seeing some very large amplitude oscillation in −dE k /dt.This is likely to be rooted in the central differencing at the interfaces, as if we change to a kinetic energy preserving formulation 35 , 36 as is displayed in Fig. 16, these oscillations are removed and the sensitivity to jitter is reduced.The enstrophy (Fig. 16b) seems to indicate that a large amount of what seemed to be resolved energy may have in fact been dispersion induced fluctuations.However, in both cases FV was able to run with grids up to j f = 0.9 and q h = 0.6382 -not shown -it appears that in these cases the added stability of the smoothing has greatly helped FV, especially in the central difference case where running without smoothing caused the case to fail even at low levels of jitter.
Before moving on it must be noted that in FR we once again see a dip in dE k /dt, which is also present to a lesser extent in both versions of FV tested.This is again linked to the instability caused by locally expanding grids, but in the case of FV the dip is smaller and is aided by by the use of smoothing.The conclusion for this comparison is that FV is somewhat resilient to degradation to the mesh quality, with the resilience coming from smoothing for central differenced FV.This allowed highly warped meshes to be run, but at the expense of accuracy with excess dissipation affecting the solution.KEP was found to be far more resilient and could even run without smoothing.FR, when run at low order, was unsuited to this problem,  but did see less degradation compared to central FV and there is the potential for p = 1 FR to equally benefit on poorer quality meshes from smoothing via a different Riemann solver, such as Roe, or from the use of KEP. 37

VI. Conclusions
Through this work we have presented a theoretical extension of the von Neumann analysis of FR to higher dimensions.This allowed improved understanding of the character of the dispersion and dissipation relations of FR as the incident angle of a wave was varied.Differences were noted between the behaviours of FR and finite differencing methods, primarily that FR saw lower variation in character with angle.The effect of higher dimensionality on the CFL limit was also investigated, with higher dimensionality causing a reduction in the CFL limit.
Investigations were then performed on deformed meshes, and the same theoretical behaviour was seen in two dimensions as in one -specifically that expanding meshes cause instability while contracting meshes cause excess dissipation.However, when coupled together, the effects can act to cancel each other out.Numerical experiments were subsequently performed using the Taylor-Green Vortex case, but with the element corner node positions jittered.Tests showed 5 th order to be more resilient to mesh poor quality meshes than 3 rd order, but in both cases the effect of localised regions of expansion are thought to be responsible for an initial increase in the kinetic energy of the solution.The appearance of smaller turbulent scales within the TGV solution as time progressed then counteracted this effect, as high wavenumbers on locally contracting regions experience excess dissipation.Lastly, comparison was made between FR and a second order FV method.It was found that FR was more resilient to mesh deformation that FV methods, however FR is far from optimal when running at second order.In both cases it recommended that kinetic energy preserving methods should be used, as it is likely that they will greatly increase resilience to mesh quality.

3
Flux and solution point layout for p = 3 in Ω, with corresponding left and right Huynh correction functions.Discontinuous primitive polynomials and the interpolated values at the flux points for adjacent sub-domains.

Figure 1 :
Figure 1: Point layout in Ω for p = 3 and cell interface topology.

( a )
Schematic showing inclined plane wave passing through a cell with a geometrically transformed rectilinear stencil.Schematic showing inclined plane wave passing through a cell with a geometrically transformed rectilinear stencil.

Figure 3 :
Figure3: Dissipation of the primary mode for 2D upwinded FR, with Huynh g 2 corrections, at various orders.The radial distance is the normalised wavenumber (including the effect of angle), and the azimuthal distance is the angle of incidence to the element.

Figure 4 :
Figure4: Dissipation of the primary mode for 2D upwinded FR, with Huynh g 2 corrections, at various orders.The radial distance is the normalised wavenumber (including the effect of angle), and the azimuthal distance is the angle of incidence to the element.

9 Figure 5 :
Figure 5: Dissipation and dissipation of the primary mode for 2D upwinded FR, with Huynh g 2 corrections, for different grid expansion factors.The radial distance is the normalised wavenumber (including the effect of angle), and the azimuthal distance is the angle of incidence to the element.The solid black line on the dissipation plots is the contour of zero dissipation.

Figure 7 :
Figure 7: CFL limit for 2D linear advection with FR (p = 4) using Huynh correction function, showing variation with θ and γ x for some set values of γ y .Time integration is RK44.

Figure 8 :
Figure 8: CFL limit for 2D linear advection, at several orders.Varying correction function parameter, ι, and angle θ.Time integration is RK44.The VCJH correction function is shown as a dashed black line.

( a )
Variation of mode condition number with angle and wavenumber for upwinded 2D FR, p = 2, with Huynh g2 correction functions.(b) Variation of mode condition number with angle and wavenumber for upwinded 2D FR, p = 2, with Nodal DG correction functions.Variation of mode condition number against order, p, for upwinded 2D FR, θ = 45 • , with Huynh g2 correction functions.Variation of mode condition number against order, p, for upwinded 1D FR, with Huynh g2 correction functions.

Figure 9 :
Figure 9: Condition number of various linear FR configurations

Figure 11 :
Figure 11: Example slices through a 3D hexahedral mesh to illustrate the mesh quality metric.
(a) Selected turbulent kinetic energy dissipation.(b) Variation of turbulent kinetic energy dissipation with jitter.Dashed contour at zero dissipation.

Figure 12 :
Figure 12: Effect of jitter on turbulent kinetic energy dissipation of the TGV (R e = 1600) for FR, p = 2, with Huynh g 2 correction functions on a 40 3 element mesh.Explicit time step size is ∆t = 1 × 10 −3 .

( a )
Selected turbulent kinetic energy dissipation.(b) Variation of turbulent kinetic energy dissipation with jitter.Dashed contour at zero dissipation.

Figure 13 :
Figure 13: Effect of jitter on turbulent kinetic energy dissipation of the TGV (R e = 1600) for FR, p = 4, with Huynh g 2 correction functions on a 24 3 element mesh.Explicit time step is ∆t = 1 × 10 −3 .

Figure 15 :
Figure15: Comparison of FR, p = 1 with DG correction functions with a second order central FV scheme with L2 Roe smoothing both with 170 3 degrees of freedom and ∆t ≈ 5 × 10 −4 .A reference DNS solution is provided from.31

Figure 16 :
Figure16: Comparison of FR, p = 1 with DG correction functions with a second order KEP FV scheme with L2 Roe smoothing both with 170 3 degrees of freedom and ∆t ≈ 5 × 10 −4 .A reference DNS solution is provide from.31