Temporal Stabilisation of Flux Reconstruction on Linear Problems

Temporal Stabilisation of Flux Reconstruction on Linear Problems Will Trojak∗, Rob Watson† and Paul G. Tucker‡ arXiv:1805.06313v1 [math.NA] 16 May 2018 Department of Engineering, University of Cambridge, Cambridge, UK, CB2 1PZ Filtering is often used in Large Eddy Simulation with a global filter width, instead here a filter width in the reference domain of high order Flux Reconstruction is considered. It is shown via Von Neumann analysis how filtering effects the dispersion and dissipation of the scheme when spatially and temporally discretised. With it being shown that filtering stabilises the scheme temporally by upto 25% for forth order FR. The impact of filtering on error production is calculated, highlighting the reduction in convective velocity caused and showing numerically the impact on order of accuracy. Finally, the turbulent Taylor-Green case is used to understand the effect of reference domain filtering on the transition to turbulence, and a filter Reynolds number is defined that is shown to be useful in understanding the effect of filtering on simulations. Nomenclature Roman βn weight of nth mode β array of mode weights a convective velocity c(k) wavespeed at wavenumber k γ grid geometric expansion factor C+1 downwind cell FR matrix Γ diagonal eigenvalue matrix C0 centre cell FR matrix δj mesh spacing C−1 upwind cell FR matrix ι VCJH correction function parameter D first derivative matrix Λ normalised diagonal eigenvalue matrix ξ local computation spatial variable en nth element solution point error σ Gaussian filter width f flux variable ς Gaussian bump solution width hl & hr left and right correction functions τ time step Jn nth cell Jacobian solution domain knq solution point Nyquist wavenumber, (p + 1)/δj Ω Ωn nth solution sub-domain k̂ knq normalised wavenumber, [0, π] Ω̂ reference sub-domain Ma Mach number Subscript p solution polynomial order •l variable at left of cell Pr Prandtl number •r variable at right of cell Q spatial scheme matrix Superscript Re Reynolds number ˆ• transformed value R update matrix • locally fitted polynomial of value S filter matrix •c common value at interface u primitive variable •δ discontinuous value W eigenvector matrix Greek I. Introduction Flux Reconstruction, a high order method originally proposed by Huynh,16 has proved to be largely robust1 and is gaining interest as a commercial CFD tool — particularly as FR seems to overcome some of ∗ PhD Candidate, AIAA Student Member, wt247@cam.ac.uk Fellow, AIAA Member ‡ Professor, AIAA Assoc. Fellow † Research 1 of 16 American Institute of Aeronautics and Astronautics the computational limitations confronted by lower order methods when performing Large Eddy Simulation.2 LES itself is of industrial interest as it offers the possiblity of performing without the turbulence modelling uncertainty associated with Reynolds-Averaged approaches, giving improved aerodynamic and aerothermal predictions as a result. However, a hurdle that must be confronted by all high order methods is the limitations imposed by the stability of temporal integration, and the effect that this can have on the temporal integration. This can mean that explicit time steps, in the case of FR, are between half and one fifth of those expected from a low order method such as finite volume. These lower CFL limits introduced by FR have the potential to become limiting when running the larger cases that FR makes possible and LES/DNS demand. It has long been known that filtering can aid the temporal stability of schemes, for example Fischer and Mullen3 showed the stabilising effect of sharp spectral filtering on spectral element methods. The authors, in this paper, aim to present the generalised fully discrete von Neumann analysis of Flux Reconstruction and present how, on uniform grids, the harmonic modes of the system impact the stability of the function projection. The use of filtering as a method of instability amelioration is then introduced and the impact on the fully discretised scheme is investigated, along with the effect that reference domain filtering has on the implicit filter of the scheme. A key tool used in the assessment of filtering on both semi- and fully-discrete schemes is an analytical error study, with an extension made here so that the fully discretised version may be considered. Finally, numerical tests are performed to study the effect of filtering when implemented, with the aim of both validating the analytical results and showing how these methods impact the calculation of fully turbulent flows. II. Flux Reconstruction As discussed in Section I, the numerical method under investigation is Flux Reconstruction. We will now give a brief overview of the mechanics of the scheme by presenting a 1D first order example. However, further information on the derivation, stability, and implementation can be found in.4–12 Let the domain of the solution be defined as Ω. We then perform domain decomposition such that we have sub-domains, Ωi , that fulfil: N N \ [ Ωi = ∅ (1) Ωi and Ω= i=1 i=1 We then define a reference domain, Ω̂ ∈ [−1, 1], and associate a Jacobian transformation from one domain to the other such that Ji : Ωi 7→ Ω̂i . Within the sub-domains and reference sub-domain we position solution points such that xj,p = {x0 , . . . , xp : xi ∈ Ωj } and ξ p = {ξ0 , . . . , ξp : ξi ∈ Ω̂}. Collocated points on the edge of the sub-domain are also positioned and in 1D this gives ξ f = {−1, 1}. The positioning of points is of great importance and was investigated for tetrahedra by Witherden.13 With the domain of the solution established we now introduce both the equation to be solve and the method to solve it. ∂u ∂f + =0 ∂t ∂x (2) If we denote a variable in the reference sub-domain as ûi , then ûi = Ji ui , and for each sub-domain a discontinuous polynomial may be defined as: ûδi = p X ûδi (ξj )lj (ξ) (3) j=0 where p is the order of the polynomial fit and lj (ξ) is the j th Lagrange polynomial.  p  Y ξ − ξk lj (ξ) = (1 − δjk ) ξj − ξk (4) k=0 The flux polynomial may be similarly defined and subsequently interpolated to the boundary of the element. Having placed boundary nodes on adjacent elements such that they are collocated, a common interface flux may be calculated. The aim of this is to correct the discontinuous polynomial such that a C 0 continuous solution is formed. The common interface flux can be calculated by a number of methods, however when 2 of 16 American Institute of Aeronautics and Astronautics solving for hyperbolic equations it is important that the upwind direction is considered in the calculation, and a Riemann problem should be approximately solved. More details are available from Toro.14 With the common interface flux at the left and right interfaces defined(fˆlI and fˆrI respectively) the flux correction can be applied. This is performed using a correction flux defined such that: fˆC (ξ) = (fˆlI − fˆlδ )hl (ξ) + (fˆrI − fˆrδ )hr (ξ) (5) where hl and hr are the left and right correction functions. It is obvious that they have boundary conditions: hl (−1) = hr (1) = 1 and hl (1) = hr (1) = 0. Further information on the definition of the correction function may be found in.4–6 Upon differentiating Eq.(5) the transformed conservation equation may be written as: ∂ fˆj ∂ ûj =− ∂t ∂ξ p X I δ dhl (ξ) I δ dhr (ξ) δ dlj (ξ) − (fˆi,l − fˆi,l ) − (fˆi,r − fˆi,r ) =− fˆi,j dξ dξ dξ j=0 (6) (7) With the equation semi-discretised, the solution may now be advanced in time via a suitably selected temporal integration method. III. Analytical Methods We begin the investigation into the temporal stabilisation of FR via a theoretical numerical analysis. Two analytical techniques will be presented so that the character of FR under the effect of filtering may be assessed. These are von Neumann analysis and an accompanying error analysis. We start by introducing the von Neumann analysis which will be performed in a manner similar to that introduced by Lele,15 Huynh,16 and Vincent et al.17 The error study aims to extend the work of Hesthaven et al.18 and Asthana et al.,19 with an extension made to show the fully-discretised temporal behaviour of harmonic solutions. III.A. Von Neumann Analysis The von Neumann analysis present will be focused on the linear advection equation. ∂u ∂u +a =0 ∂t ∂x (8) A matrix form of this equation consistent with FR operators was presented by Huynh16 and Vincent et al.17 and we will follow this, but with the added generalisations for arbitrary meshes as in Trojak et al.,1 which leads to the following form: ∂ ūj −1 −1 = −Jj+1 C+1 ûj+1 − Jj−1 C0 ûj − Jj−1 C−1 ûj−1 ∂t C+1 = (1 − α)hr T ll T T C0 = D − αhl ll − (1 − α)hr lr T C−1 = αhl lr (9) (10) (11) (12) where D is the differentiation matrix at the solution points, hl and hr are the gradients of the left and right correction functions evaluated at the solution points, and ll and lr are interpolation vectors to the left and right interfaces. The spectral response of the system is then obtained by applying a trial Bloch wave solution and vectorising the result to give: u(x, t) = vei(kx−ωt) uj = vj e ik( ξ+1 2 δj+1 +xj −c(k)t) (13) (14) where δj = xj − xj−1 and for wavenumber k. Substituting this into Eq.(9), FR may be cast into the semi-discrete form: ∂ ūj = Qūj (15) ∂t 3 of 16 American Institute of Aeronautics and Astronautics −1 −1 C−1 e−ikδj . For the case of generalised interfaces it follows that Q = −Jj+1 C+1 eikδj+1 − Jj−1 C0 − Jj−1 This may then be transformed into the fully discrete form by application of a suitable temporal integration method, thereby giving what is often called the update equation: ūj (t + τ ) = R(Q)ūj (t) ūj (nτ ) = Rn (Q)ūj (0) R33 = I + (τ Q)2 (τ Q)3 (τ Q)1 + + 1! 2! 3! (16) (17) (18) where R is the update matrix and is a function of Q that is dependent on the temporal integration used. In a method similar to that used on the semi-discrete form, the Bloch wave solution can be applied and discretised, leading to the result: ikτ e|−ik(c−1)τ {z } v = e Rv e λ −ik(c−1)nτ v = eiknτ Rn v (19) (20) The result of Eq.(19) was first reported by Asthana and Jameson,20 but the accompanying analysis was limited, and was followed up by Vermeire and Vincent21 where the application of FR to implicit LES was considered. Equation (19) can be seen to be an eigenvalue problem of the update matrix, with v being the eigenmodes of the fully discretised scheme. The eigenvalues of the system are then temporally shifted and scaled solutions, with the fully discretised wave speed recovered as: i log (λ) +1 kτ ω ′ = c(k; τ )k c(k; τ ) = (21) (22) This can provide useful insight into FR, as previously the semi-discreted form led to a disconnect between the numerically realised and analytical dispersion/disspation in all but highly temporally resolved waves. In contrast, Eq. (21) gives a more complete picture. A further result is that scheme stability can then be understood to be when ℑ(c(k; τ )) < 0 ∀ k ∈ R, which is equivalent to |λ(k; τ )| 6 1 ∀ k ∈ R. This will give the same results as von Neumann theorem for stability,22 but will give added information about the mechanism driving the scheme unstable. III.B. Error Analysis Von Neumann analysis can be insightful when trying to understand how a method will react when applied to problem. However, this only gives a snapshot of the behaviour, and it can also be of use to understand how a solution will develop as the dispersion/dissipation transfer function is repeatedly applied and integration error has iterations in which to grow. We also consider, therfore, how the error developes with time, which was initially investigated for Nodal DG via FR by Asthana et al.19 In the semi-discrete sense this may be can be done by utilising exponential integration of Eq.(15), as Q forms a linear spatial differentiation operator, see Pope.23 The first step, following the method of,19 is to diagonalise the operator matrix: −1 −1 Q = W Q ΓQ W Q = ikWQ ΛQ WQ (23) where W is a matrix of eigenvectors, Γ is a diagonal matrix of eigenvalues and Λ is a diagonal matrix of complex convective velocities. From here the temporal integration can be performed: −1 uj (0) uδj (t) = exp (ctQ)u(0) = WQ exp (ikctΛQ )WQ (24) where the definition of the exponential has been used to simplify the form by considering W−1 W = I. By applying the initial condition as a sum of the eigenvectors, again following Asthana:19 uj (0) = exp (ikxj )WQβ = exp (ikxj ) p X n=0
wQ,n βn = exp ikJj (ξξ + 1) exp (ikxj ) 4 of 16 American Institute of Aeronautics and Astronautics (25) where β is an array of weights associated with the reconstruction of the initial solution from the eigenvectors wQ . Substituting this into Eq.(24) leads to: β = exp (ikxj ) uδj (t) = WQ exp (ikctΛQ ) exp (ikxj )β p X exp (ikctλn )βn wQ,n (26) n=0 Lastly, to understand the rate convergence, the error is calculated via analytically evolving the initial condition such that: p X βn wQ,n (27) uj (t) = exp (ik(xj − ct)) n=0 and hence the semi-discretised error is: ej (t) = uδj (t) − uj (t) = exp (ik(xj − ct)) p  X n=0 
exp ikct(λn + 1) − 1 βn wQ,n (28) From here the it is necessary to calculate the eigenvectors and values and then decompose the solution to calculate βi , from Eq.(25). If instead the diagonalised form of the discrete FR operator, Eq.(23), is substituted into the fully discretised version of the update equation, Eq.(17), then: !n r X (ikτ )m δ m −1 uj (nτ ) = WΛ W uj (0) (29) m! m=0 !n r X (ikτ )m δ m −1 β exp (ikxj ) WΛ W Wβ (30) uj (nτ ) = m! m=0 And therefore, the fully discrete error can be derived as: ej,r (nτ ) = uδj (nτ ) − uj (nτ ) = r X (ikτ )m WΛm W−1 m! m=0 !n ! β exp (ikxj ) − exp (−ikcnτ )I Wβ (31) where r is the number of RK temporal integration sub-steps and I is the identity matrix. What this shows is how the truncation of the temporal integration introduces error, and, interestingly that lim (ej,r )r→∞ = ej , by remembering the recursive definition of exp (z). An attempt can be made to form a temporal integration method that reduces the error via the introduction of factors at each RK step, however it was found that this method reduces the CFL limit by an order of magnitude due to the weights varying for each wavenumber, and does not seem to be feasible in the general case. IV. Filtering Filtering can be used in LES to suppress numerical errors by excluding erroneous high wavenumbers, such that sub-grid scale convergence is separated from numerical convergence, see Ghousal.24 This allows explicit LES to converge upon the spatially filtered NSE. This does, however, require the domain-level setting of filter width. This is not directly congruous with FR, where operations are defined by matrices that are spatially invariant due to the use of a reference domain. Therefore, it is proposed that reference domain filtering is used such that the computational benefit of spatially invariant operators may still be realised, whilst still using filtering to modify the underlying scheme’s behaviour. Filtering was recently investigated by Asthana et al.,25 where it was used with a globally defined filter width to suppress aliasing errors when non-linear equations are considered. It was proven that filtering of the solution in FR is equivalent to added dissipation and super-dissipation. It is therefore thought that local reference domain filtering may allow control over the linear advection dissipation. We begin by defining the filter to be studied and the possible means of application. There are a plethora of possible filters, with a list of important ones given by Pope,26 for which the inverse Fourier transform can be used to calculate the spatial filter kernel, s(r) = F −1 {S(k)}. We will primarily focus on the Gaussian 5 of 16 American Institute of Aeronautics and Astronautics filter kernel, as this is a canonical example which can also be presented in a C p continuous manner which is believed may make the result more stable. The filter matrix can be defined as: r   6|ξξ i − ξ j |2 6 (32) exp − Si,j = πσ 2 σ2 where the filter width is σ and the transformed solution point coordinates are ξ . To ensure that the filter does not act as a net source or sink the integrated effect of each solution point is used to renormalise the effect. The filter can be applied in three ways, with varying affect. The first approach is to apply the filter to the entire spatial scheme, as in: ∂ ūj = SQūj (33) ∂t The filter may also be just applied to the differentiation operator, defining C0 as: C0 = SD − αhl T ll − (1 − α)hr T lr (34) which would allow the correction function to introduce higher wavenumber content but suppress it from the differentiation operator. Finally, the filter can be exclusively applied to the correction function, defining C0 , C−1 and C+1 as: C+1 = (1 − α)Shr T ll T (35) T C0 = D − αShl ll − (1 − α)Shr lr T C−1 = αShl lr (36) (37) which would have the effect of convolving the correction function by the filter therefore preventing the correction function from passing high wavenumber information from adjacent elements. The advantage of the filter approaches proposed in Eq.(35-37) is that the scheme can be potentially filtered in specific areas that such that the ingress of spurious or unwanted modes maybe suppressed. V. V.A. Results Von Neumann Analysis The main aim of this paper is to investigate the temporal stability of the scheme and hence the dissipation of the scheme is of primary importance. However, as was observed by Vermiere and Vincent,21 the group velocity of waves near the Nyquist limit becomes very high when near to the CFL limit. Hence, a further result of interest would be the control of this. Initially evaluating the performance of the scheme without filtering, using RK33 temporal integration, it can be seen by comparing Fig.1 and Fig.2, with τ either side of the CFL limit for uniform grids, that, as expected, the temporal instability originates from higher wavenumbers. However, instability is spread to a wider band of wavenumbers via the other eigenmodes. However, for general solutions a multiplicity of modes will form a solution. Therefore, the instability shown Fig.2 will likely result in growth of errors and ultimately lead to failure. But, what can be seen is that the demanded harmonic mode only becomes anti-dissipative at high wavenumbers after a region of large dissipation. This makes the suppression of unstable high wavenumber modes via filtering viable, as solution propagation into this high wavenumber region is already bottle-necked by the high dissipation for k̂ ≈ 5π/8. Hence, limiting the effect of filtering on solution accuracy at high wavenumbers. Now we apply a Gaussian filter to FR as in Eq.(32 & 33). Analytically, from Fig. 3, with σ = 0.6, it can be seen that the scheme has been stabilised for a CFL number that was previously unstable. The impact of the filter is clearly visible, with high wavenumbers in primary harmonic showing added dissipation. Yet, the region of wavenumbers that experience high dissipation has now had that comparatively reduced. The Nyquist diagrams in Figs.2 & 3 are useful in visualising what has happened in this case. These show that although peak dissipation has reduced, the overall dissipation has increased, as might be expected. The trade off that has been made is evident from Fig. 1a & 3a, where the phase and group velocities for a central band of wavenumbers (π/4 < k̂ < π/2) has been reduced. This will impact the convective velocity of solution which will become more evident during later error studies. This numerical analysis is concluded by considering the von Neumann conditions for ensuring that the operator performs a contraction transformation when temporally discretised, i.e. ρ(R) 6 1. The aim of this 6 of 16 American Institute of Aeronautics and Astronautics π Primary 0 -2.5 0 ℑ(ω ′ ) ℜ(ω ′ ) π/2 −π/2 -7.5 spectral Primary −π 0 π/4 -5 π/2 3π/4 -10 π 0 π/4 k̂ (a) Dispersion π/2 k̂ 3π/4 π (b) Dissipation Figure 1: Fully discretised analysis of FR using Huynh correction functions, p = 4 with RK33 time integration, τ = 0.166, close to the temporal stability limit. π Primary 0 -2.5 0 ℑ(ω ′ ) ℜ(ω ′ ) π/2 −π/2 -7.5 spectral Primary −π 0 π/4 -5 π/2 3π/4 -10 π 0 π/4 k̂ (a) Dispersion π/2 k̂ 3π/4 π (b) Dissipation Figure 2: Fully discretised analysis of FR using Huynh correction functions, p = 4 with RK33 time integration, τ = 0.17. Showing instability due to exceeding the CFL limit. 7 of 16 American Institute of Aeronautics and Astronautics π Primary 0 -2.5 0 ℑ(ω ′ ) ℜ(ω ′ ) π/2 −π/2 -7.5 spectral Primary −π 0 π/4 -5 π/2 3π/4 -10 π 0 π/4 π/2 k̂ k̂ (a) Dispersion 3π/4 π (b) Dissipation Figure 3: Fully discretised analysis of FR using Huynh correction functions, p = 4 with RK33 time integration, τ = 0.17. Gaussian spatial filtering is applied as in Eq.(33) with σ = 0.6. (a) Eq.(33) (b) Eq.(34) (c) Eq.(35-37) Figure 4: FR, p = 4, RK33, showing variation of CFL limit with filter width, σ, and correction function parameter, ι for various methods of filter application. The dotted line is the ι+ correction of Vincent et al.17 that give peak CFL limit, with the dashed contour at the CFL limit associated with ι+ (CFL = 0.2118). 8 of 16 American Institute of Aeronautics and Astronautics analysis is to understand the extent to which filtering can improve temporal stability and to understand the maximum stable filter width for the various methods of filter application. For this we will consider only the single parameter energy stable flux reconstruction scheme.4 It was previously found by Vincent et al.17 that there is an optimal correction function that recovers super-convergence, with ι = ι+ , and gives peak the CFL limit of VJCH correction functions when spatial temporal discretisation is considered. Through Fig. 4 the temporal stability of the various methods of filter application are explored. It is clear that filtering the full scheme provides the most temporally stable results, with peak CFL also occurring at ι+ , giving a maximum 25% boost in CFL limit. What is made clear from Fig. 4b is that the application of filtering to solely the differentiation operator has only a limited impact on temporal stability. This is consistent with the order of the differentiation operator being lower than that of the correction and that instability is instigated in the range of wavenumbers that are resolved by the higher order. V.B. Analytic Error Study The effect of a Gaussian filter in the spectral domain is very well understood and it can be seen from Section V.A that filtering adds dissipation to the high wavenumbers. When considering the temporal stability of the scheme this is very beneficial, however, it is important to consider the effect of filtering on the accuracy of the scheme. Following on from the method introduced in Section III.B, the semi-discrete error can be calculated, the results of which are shown in Fig. 5. In addition to this analysis, the dissipative half life of (a) σ = 0 (b) σ = 0.6 Figure 5: Analytical error against time and normalised wavenumber. This is for FR with p = 4 and interface upwinding (α = 1) on a uniform mesh. The dashed line is the decay time period for the primary harmonic and the other time periods are shown in dotted lines. the modes is calculated in accordance with Asthana et al.:19 τ1/2,n = − 1 ckJj−1 ℑ(λn ) (38) where τ1/2,n is the half life of the nthL mode (the contours of which are shown in Fig. 5), and λn is the nth p diagonal entry of Λ. Such that Λ = n=0 λn . From initially studying the unfiltered scheme error, Fig.5a, the resolvable wavenuimbers can be broken down into three regions, firstly a region where the half-life of the primary mode is long. In this region waves are well resolved even in the long time integration. Secondly, at slightly higher wavenumbers the primary mode decay time decreases with out a large increase in the half-life of the other modes, in this region the initially well resolved solution will decay uniformly. Lastly at the highest wavenumbers there is a region where the primary mode half-life is low but the half-0lives of the the other modes is large. This leads to the the decay of the primary mode, but slower moving secondary modes do not decay as readily. Hence, stripping error in the error plot. Studying the filtered case, Fig.5b, much of the same behaviour can be seen, however there is an intermediate band in the range of wavenumbers from π/4 < k̂ < π/2 where the error seems to become very high. This, by comparison with Fig.3a, is the region in which there is dispersion undershoot. Hence in this region the phase and group velocities are lower and the error will increase as the phase difference increases. As the 9 of 16 American Institute of Aeronautics and Astronautics primary mode half-life is broadly unaffected in this region it is expected that at some time later the wave will have locked on again with this analytic solution and the error will decrease. For the semi-discrete case the remaining behaviour at higher wavenumbers looks much the same as the unfiltered case. (a) σ = 0, τ = 0.1 (b) σ = 0.6, τ = 0.1 (c) σ = 0, τ = 0.17 (d) σ = 0.6, τ = 0.17 Figure 6: Analytical error against time and normalised wavenumber, with various filters. This is for FR with p = 4 and interface upwinding (α = 1) on a uniform mesh, when temporally and spatially discretised with RK33 temporal integration. Two diffenrent time steps are used, which are stable and unstable for the Huynh correction functions used here. The dashed line is the decay time period for the primary harmonic and the other time periods are shown in dotted lines. Using Eq.(31) the fully discretised error can be explored, along with how discrete temporal integration will effect the error production mechanisms. Comparison of Fig.6a and Fig.6c shows the mechanism of solution divergence when the CFL limit is exceeded. There are five distinct zones of wavenumbers that become unstable and widen as the time progresses. These regions of instability align with the unstable dissipation that is displayed in Fig.2b. If the effect of filtering is then considered, as in Fig.6b&6d, similar behaviour of the error is observed as in the case of exponential filtering. However, at high wavenumbers the region where the dissipation in Fig.3b approaches zero is highlighted as a region where the primary mode half-life is slightly increased. Clearly this is benificial, as it delays the onset of spurious modes. What may also be noted at these higher wavenumbers is that the higher frequency spurious modes are less pronounced when discrete temporal integration is considered, which originates from the discrete temporal integration also behaving some what like a filter. V.C. Numerical Experiments Validation of the behaviour analytically exhibited is performed via a wave convection test case, in which a wave is input into 1D domain. Interface upwinding is utilised which permits the use of a blow out boundary condition at outlet. During this test three operating conditions are investigated, all of which are at a CFL number that is unstable for the baseline scheme. Two of the waves are for unfiltered schemes and the 10 of 16 American Institute of Aeronautics and Astronautics third has a Gaussian filtered applied with σ = 0.6. The results are displayed in Fig.7a and it is clear that application of the filter has stabilised this case, which is in agreement with the earlier findings. To highlight the effect of filtering on convective velocity, we present a test in which a Gaussian bump is convected via the linear advection equation through a periodic domain, see Fig.7c. If the bump is prescribed to be sufficiently narrow this will excite a significant proportion of the wavenumber space, highlighting any instability. The initial condition was therefore chosen to be:
u = exp − (x − 0.5)2 /ς 2 (39) where ς = 0.2. Comparison of the filtered and unfiltered bump cases clearly shows that the convective velocity has been reduced in the filtered cases, with the diffusive effect of the filter clearly displayed as a reduction in the peak value. This reduction in the convective velocity is consistant with the reduction in the group and phase velocities shown in Fig.3a. A filter width was also tested that exceeded that the maximum value found to be stable in Section V.A. This was motivated by the need to understand the mechanism driving the scheme unstable, and it is clear from Fig.7b that at x ≈ 0.1 the conserved variable has dipped below zero. Hence, at this level of filtering, the filter is redistributing mass faster than incoming mass can replace it. However, this does seem to be gradient dependant originating from the application of the filter to the differentiation operator, D and exemplified by the negative derivative on the upwind downwards slope seeming to lead to this region of negative conserved variable. The analytic solution of this problem permits the calculation of the error which is shown in Fig. 7c & 7d. The effect of filtering on the error of the scheme has not been significant, causing an upwards shift in the error together with an a small decrease in the order accuracy. This decrease in the order accuracy is to be expected as it has previously been shown by Asthana et al.25 that filtering is equivalent to additionally solving various orders of diffusive derivatives. Although it could be proposed that in order to increase the time step size the grid could be coarsened, in doing so the Nyquist wavenumber will also be reduced and hence the advantage of this method of filtering, due to the already high dissipation at k̂ ≈ 5π/8, little spectral information is lost and that which is lost, at the very highest wavenumbers, cannot have significant impact on solution due to the impassable dissipation gap seen in Fig. 1 & 2. What may also be observed are the diverging errors of the unstable filter width, which is consistent with the findings of Fig.4a. V.D. Taylor-Green Vortex Finally, we apply the methodology to an investigation more representative of engineering applications, the Turbulent Taylor-Green Vortex. The primitive variables are initialised as:       y z x cos cos (40) u = U0 sin L L L       x y z v = −U0 cos sin cos (41) L L L w=0 (42) p = p0 + ρ= p RT0 ρ0 U02 16  cos  2x L  + cos  2y L  cos  2z L  +2  (43) (44) The specific case investigated is: Re = ρ0 U 0 L , µ U0 = 1, Pr = 0.71 = ρ0 = 1, µγR , κ(γ − 1) p0 = 100, 27 Ma = 0.08 = √ R = 1, U0 γRT0 γ = 1.4 28 (45) (46) This case follows the case set up of Brachet et al. and DeBonis. Their DNS data is used as comparison, with the main concern of this investigation being the domain averaged kinetic energy and kinetic energy dissipation. Z 1 ρui ui dΩ (47) Ek = 2ρ0 Ω Ω 11 of 16 American Institute of Aeronautics and Astronautics 1 CFL = 0.25 CFL = 0.2, σ = 0.6 CFL = 0.2 u(x; t = 10) 0.75 u(x = 1; t) 1 0.5 1.7π 0 1.8π 0.5 0.25 exact CFL = 0.2, σ = 0.7 CFL = 0.2, σ = 0.6 CFL = 0.1 0 -0.5 0 0.25 0.5 -1 0 π/2 3π/2 π 2π t (a) Output in time at x = 1 with k = 1. σ=0 10 10−4 10−5 0.6 −6 10−7 10−8 0.8 −1 10−2 0.7 10−3 10−4 σ 0.7 10−3 |e(nt = 103 )|2 −1 σ |e(nt = 103 )|2 100 0.8 10−2 10 1 (b) Spatial slice of Gaussian bump, with periodic boundary conditions, after 10 cycles on a unit domain. σ=0 100 10 0.75 x 10−5 0.6 10−6 10−7 101 102 10−8 0.5 101 102 0.5 nx nx (c) Error calculated of Gaussian bump convection after 103 time iterations. p = 4 with Huynh correction, RK33 and constant CFL = 0.167, limit for Huynh correction. (d) Error calculated of Gaussian bump convection after 103 time iterations. p = 4 with Huynh correction, RK44 and constant CFL = 0.189, limit for Huynh correction. Figure 7: Numerical tests of FR using Huynh correction functions, p = 4. 12 of 16 American Institute of Aeronautics and Astronautics To isolate the effect of filtering for various various Re , the Taylor-Green test case is used with cell Re fixed at Re,cell = 40. (a) Re = 1600, N = 403 (b) Re = 800, N = 203 (c) Re = 400, N = 103 Figure 8: Taylor-Green Vortex turbulent kinetic energy dissipation. For p = 4, using Rusanov flux an the inviscid common interfaces and BR1 at the viscous common interfaces and DG correction functions. The method of temporal integration is RK44 with dt = 1 × 10−3 . Between each case the cell Re was held constant at Re,cell = 40. The reference data is that of Brachet et al.27 What should be clear from Fig.8a is that at this resolution the transition to turbulence is well captured , with only minor deviation from peak kinetic energy dissipation. This is a result that is shared by Bull and Jameson.29 However, as we move to lower Re , Fig.8c, the deviations grow, which may be explained by the increased sensitivity at lower Re to changes in the apparent viscosity. With a filter added the effect on lower Re cases also becomes more pronounced, this is again due to the increased sensitivity of lower Re flows to changes in apparent viscosity. However, there is another effect that impacts the results here, that the the filter spectral width is constant in the reference domain and so depends on the physical size of the element. The result of this is that as the grid becomes coarser the filter becomes more dissipative, i.e. the filter Re depends inversely on element length. From the filter viscosity of derived by Asthana et al.25 and the filter definition of Eq.(32) the filter Re can be defined as: Re,filter = 24ρuτ σ2 h (48) where τ is the time step size and h is the grid spacing. To investigate further the effect of grid spacing and filter width on the simulation of the full Navier-Stokes equations, a series of different grid spacings and filter width are considered for constant global Re = 1600. What can be seen from Fig.9 is that there is a marked impact of the filtering as Re,filter reduces below one, with exact transition dependent on the resolution of the grid relative to the flow. Furthermore, the 13 of 16 American Institute of Aeronautics and Astronautics (a) Re = 1600 on grids 403 , 323 , 243 , and 163 with σ = 0.3. (b) Re = 1600 with 323 elements with σ ∈ {0, 0.3, 0.35, 0.4}. Figure 9: Turbulent Taylor-Green vortex kinetic energy dissipation for FR, p = 4, using nodal DG corrections, with Rusanov inviscid and BR1 viscous common interface flux calculations. The temporal integration method is RK44 with τ = 1×10−3 . The legend format is [Re,filter , Re,cell ]. The reference data is that of Brachet et al.27 14 of 16 American Institute of Aeronautics and Astronautics definition of filter viscosity used in the definition of Re,filter takes into account only the lowest order diffusion derivative and hence there will be some higher order dependence of filter Reynolds number on grid spacing that will also affect the exact point of transition. The physical interpretation of this result is that as the filter Re decreases, the rate at which solutions are convected through an element reduces compared to the time scale of the diffusion caused by the filter. Therefore, a reduction of Re,filter below one implies that the time scale of diffusion has become more significant than that of inertia, and so the filter diffusion becomes far more apparent. VI. Conclusions It has been shown that filtering in the reference domain can be used to alter the behaviour of the underlying FR scheme by the addition of dissipation at higher wavenumbers. This dissipation can be used to increase the temporal stability, which, for the case of p = 4 with RK33 temporal integration, can increase the CFL limit by ≈ 25%. For temporal stabilisation the most successful method of filtering FR is complete filtering of the spatial scheme, as it provides a more coherent treatment of higher wavenumbers. Hence, this form of filtering in the reference domain can be thought of as modification of the numerics for implicit LES, rather than filtering for explicit LES. Numerical tests were used to validate analytical results, with an investigation of order accuracy showing that filtering only mildly reduces the order of accuracy of the method. Finally, a filter Reynolds number is defined and, through investigation with turbulent flows, filtering is shown to have a significant impact when the filter Reynolds number reduces below a critical value, which in this case that was found to be one. Acknowledgments The support of the Engineering and Physical Sciences Research Council of the United Kingdom is gratefully acknowledged. References 1 Trojak, W., Watson, R., and Tucker, P. G., “High Order Flux Reconstruction on Stretched and Warped Meshes,” AIAA SciTech, 55th AIAA Aerospace Sciences Meeting, Grapevine Texas, 2017, pp. 1–12. 2 Vincent, P. E., Witherden, F., Vermeire, B., Park, J. S., and Iyer, A., “Towards Green Aviation with Python at Petascale,” International Conference for High Performance Computing, Networking, Storage and Analysis, SC , No. November, 2017, pp. 1–11. 3 Fischer, P. and Mullen, J., “Filter-Based Stabilization of Spectral Element Methods,” Comptes Rendus de l’Académie des Sciences - Series I - Mathematics, Vol. 332, No. 3, 2001, pp. 265–270. 4 Vincent, P. E., Castonguay, P., and Jameson, A., “A New Class of High-Order Energy Stable Flux Reconstruction Schemes,” Journal of Scientific Computing, Vol. 47, No. 1, 2010, pp. 50–72. 5 Vincent, P. E., Farrington, A. M., Witherden, F. D., and Jameson, A., “An Extended Range of Stable-SymmetricConservative Flux Reconstruction Correction Functions,” Computer Methods in Applied Mechanics and Engineering, Vol. 296, 2015, pp. 248–272. 6 Trojak, W., “Generalised Sobolev Stable Flux Reconstruction,” Vol. Arxiv:1804, No. 04714, 2018, pp. 1–19. 7 Witherden, F. D., Vincent, P. E., and Jameson, A., High-Order Flux Reconstruction Schemes, Vol. 17, Elsevier B.V., 1st ed., 2016. 8 Jameson, A., Vincent, P. E., and Castonguay, P., “On The Non-Linear Stability of Flux Reconstruction Schemes,” Journal of Scientific Computing, Vol. 50, No. 2, 2012, pp. 434–445. 9 Witherden, F. D., Farrington, A. M., and Vincent, P. E., “PyFR: An open source framework for solving advection-diffusion type problems on streaming architectures using the flux reconstruction approach,” Computer Physics Communications, Vol. 185, No. 11, 2014, pp. 3028–3040. 10 Castonguay, P., Williams, D. M., Vincent, P. E., and Jameson, A., “Energy Stable Flux Reconstruction Schemes for Advection-Diffusion Problems,” Computer Methods in Applied Mechanics and Engineering, Vol. 267, No. 1, 2013, pp. 400–417. 11 Williams, D. M. and Jameson, A., “Energy stable flux reconstruction schemes for advection-diffusion problems on tetrahedra,” Journal of Scientific Computing, Vol. 59, No. 3, 2014, pp. 721–759. 12 Ou, K., Vincent, P. E., and Jameson, A., “High-Order Methods for Diffusion Equation with Energy Stable Flux Reconstruction Scheme,” AIAA, Vol. 49, No. January, 2011, pp. 1–19. 13 Witherden, F. D., Park, J. S., and Vincent, P. E., “An Analysis of Solution Point Coordinates for Flux Reconstruction Schemes on Tetrahedral Elements,” Journal of Scientific Computing, Vol. 69, No. 2, 2016, pp. 1–16. 14 Toro, E., Riemann Solvers and Numerical Methods for Fluid Dynamicsb - A Practical Introduction, Springer-Verlag Berlin Heidelberg, Dordrecht Berlin Heidelberg London New York, 3rd ed., 2009. 15 of 16 American Institute of Aeronautics and Astronautics 15 Lele, S. K., “Compact Finite Difference Schemes With Spectral-Like Resolution,” Journal of Computational Physics, Vol. 103, No. 1, 1992, pp. 16–42. 16 Huynh, H. T., “A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods,” 18th AIAA Computational Fluid Dynamics Conference, Vol. 2007-4079, 2007, pp. 1–42. 17 Vincent, P. E., Castonguay, P., and Jameson, A., “Insights From von Neumann Analysis of High-Order Flux Reconstruction Schemes,” Journal of Computational Physics, Vol. 230, No. 22, 2011, pp. 8134–8154. 18 Hesthaven, J. S. and Warburton, T., Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, Vol. 54 of Texts in Applied Mathematics, Springer New York, New York, NY, 1st ed., 2008. 19 Asthana, K., Watkins, J., and Jameson, A., “On Consistency and Rate of Convergence of Flux Reconstruction For Time-Dependent Problems,” Journal of Computational Physics, Vol. 334, 2017, pp. 367–391. 20 Asthana, K. and Jameson, A., “High-Order Flux Reconstruction Schemes with Minimal Dispersion and Dissipation,” Journal of Scientific Computing, Vol. 62, No. 3, 2014, pp. 913–944. 21 Vermeire, B. and Vincent, P., “On the Behaviour of Fully-Discrete Flux Reconstruction Schemes,” Journal of Computational Physics, Vol. 315, No. 1, dec 2017, pp. 1053–1079. 22 Isaacson, E. and Keller, H. B., Analysis of Numerical Methods, John Wilery & Sons Ltd., New York, 2nd ed., 1994. 23 Pope, D. A., “An Exponential Method of Numerical Integration of Ordinary Differential Equations,” Communications of the ACM , Vol. 6, No. 8, 1963, pp. 491–493. 24 Ghosal, S., “An Analysis of Numerical Errors in Large-Eddy Simulations of Turbulence,” Journal of Computational Physics, Vol. 125, No. 1, 1996, pp. 187–206. 25 Asthana, K., López-Morales, M. R., and Jameson, A., “Non-Linear Stabilization of High-Order Flux Reconstruction Schemes via Fourier-Spectral Filtering,” Journal of Computational Physics, Vol. 303, 2015, pp. 269–294. 26 Pope, S. B., “Large-Eddy Simulation,” Turbulent Flows, chap. 13, Cambridge University Press, Cambridge, 1st ed., 2010, pp. 558–639. 27 Brachert, M. E., Orszag, S. A., Nickel, B. G., Morf, R. H., and Frisch, U., “Small-Scale Structure of the Taylor-Green Vortex,” Journal of Fluid Mechanics, Vol. 130, 1983, pp. 411–452. 28 DeBonis, J., “Solutions of the Taylor-Green Vortex Problem Using High-Resolution Explicit Finite Difference Methods,” 51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, , No. February 2013, 2013. 29 Bull, J. R. and Jameson, A., “Simulation of the Compressible Taylor Green Vortex using High-Order Flux Reconstruction Schemes,” 7th AIAA Theoretical Fluid Meshanics Conference, No. June, Atlanta, 2014, pp. 1–17. 16 of 16 American Institute of Aeronautics and Astronautics