Artificial Compressibility Approaches in Flux Reconstruction for Incompressible Viscous Flow Simulations

Several competing artificial compressibility methods for the incompressible flow equations are examined using the high-order flux reconstruction method. The established artificial compressibility method (ACM) of \citet{Chorin1967} is compared to the alternative entropically damped (EDAC) method of \citet{Clausen2013}, as well as an ACM formulation with hyperbolised diffusion. While the former requires the solution to be converged to a divergence free state at each physical time step through pseudo iterations, the latter can be applied explicitly. We examine the sensitivity of both methods to the parameterisation for a series of test cases over a range of Reynolds numbers. As the compressibility is reduced, EDAC is found to give linear improvements in divergence whereas ACM yields diminishing returns. For the Taylor--Green vortex, EDAC is found to perform well; however on the more challenging circular cylinder at $Re=3900$, EDAC gives rise to early transition of the free shear-layer and over-production of the turbulence kinetic energy. This is attributed to the spatial pressure fluctuations of the method. Similar behaviour is observed for an aerofoil at $Re=60,000$ with an attached transitional boundary layer. It is concluded that hyperbolic diffusion of ACM can be beneficial but at the cost of case setup time, and EDAC can be an efficient method for incompressible flow. However, care must be taken as pressure fluctuations can have a significant impact on physics and the remedy causes the governing equation to become overly stiff.


Introduction
There is a growing desire across many industries to gain detailed insight into transient flow phenomena. In many cases, however, this is difficult to achieve with experimentation due to high costs and limited access to observe key physics. Over the recent decades, there has been growing interest in high-order methods in engineering applications as, compared to low order approaches, the resolution and computational efficiencies achievable can make many intractable problems tractable [47]. A further factor in this success has a shifting desire from modelled simulations, of which Reynolds averaged Navier-Stokes (RANS) may be considered representative, to the scale-resolving methods such as large eddy simulation (LES). Although RANS has seen great success in several industrially relevant problems, such as transonic wing shape optimisation [22] and stall prediction in transonic fans [49]. The approximations of the model limit its accuracy for various problems involving transition, strong streamline curvature, and relaminarization. Durbin [9] discusses recent advances, remaining and new challenges to RANS turbulence modelling as application areas widen. Scale resolving simulations, on the other hand, attempt to directly resolve the majority of the physical length scales on the grid. Such approaches can elucidate complex flow physics, especially those rooted in the small scale motions. In explicit LES methods, the effect of unresolved length scales are modelled using a sub-grid scale (SGS) model, often coupled to a filter. An alternative paradigm that has gained wider adoption, is implicit LES (ILES) where no SGS model is used and dissipation of the small scales is handled implicitly by the numerical dissipation of the discretisation. The spatial scheme we will use to perform ILES is the high-order flux reconstruction method of Huynh [19], in the PyFR implementation [52]. * Corresponding author wtrojak@ic.ac.uk (W. Trojak); nrv@smail.iitm.ac.in (.N.R. Vadlamani); James.Tyacke@brunel.ac.uk (.J. Tyacke) ORCID(s): 0000-0002-4407-8956 (W. Trojak) A challenge occurs in the limit as Mach number, , tends to zero. The compressible flow equations become increasingly stiff, ultimately resulting in an elliptic equation for the pressure field at = 0. If the compressible Navier-Stokes equations are used in this low Mach limit, the stiffness can lead to prohibitively expensive calculations. In the case of the method of Roe [34] it also leads to excessive numerical dissipation. These problems can be alleviated by the use of a low Mach preconditioners, such as that of Turkel [45]. Alternatively, in the incompressible regime a separate solver is used to solve the pressure Poisson equation. However, it is non-trivial to produce a scalable solver for the Poisson equation, although new methods such as that of Fortunato and Townsend [13] are beginning to confront aspects of this. Furthermore, within finite element methods ensuring solution compatibility requires detailed analysis that is dependent on order and element type, among other factors [11].
A single solver is preferable to make use of established tools and optimisations, and to simplify work flows. Artificial compressibility approaches allows the use of established compressible tool to calculated state solutions of incompressible flows. The first such method was the artificial compressibility method (ACM) of Chorin [5], which can be interpreted as assuming constant entropy together with an artificial compressibility to relax the pressure and velocity field onto a divergence free solution. ACM was extended to the calculation of unsteady flows by Rogers et al. [35] by using it as a method to solve the implicit equations for each physical time step. Later Jameson [20] interpreted the relaxation as a pseudo-time dimension, using explicit time stepping to perform the relaxation. More recently, Nishikawa [29] proposed a general technique where diffusion terms are hyperbolised. This has the advantage of stability scaling with ℎ −1 -rather than ℎ −2 -for some mesh spacing ℎ. This approach fits naturally with ACM and has previously been investigated by Ahn [2] and Trojak et al. [43], who introduced novel techniques to optimise the computational implementation of ACM with hyperbolised diffusion (ACM-HD).
A major issue with ACM is the requirement to converge the pressure and velocity field for each time step, which can be costly, although some convergence acceleration methods have proved successful [24,26]. An alternative is the entropically damped artificial compressibility (EDAC) method introduced by Clausen [6]. Here entropy is not fixed, but density fluctuations are minimised. This results in a similar system of equations as ACM, with the primary difference being a pressure diffusion term. This approximation leads to time dependent equations which produce an almost divergence free flows thus enabling explicit time stepping to be used. Further studies have shown the method is effective on both model problems and turbulent channel flows [1,8,42].
These three methods, ACM, ACM-HD, and EDAC, have benefits and drawbacks. We examine there relative costs and effectiveness for relevant unsteady turbulent LES. To the authors knowledge the EDAC system has not previously been used with FR and we wish to understand its effectiveness and the effect of varying the compressibility parameter within a high-order approach. Furthermore, as the EDAC system will be applied as a conservative equation, we wish to further the understanding of the Riemann problem as it forms an important part of the FR method [28]. We also wish to better understand the effectiveness of the ACM approaches and how the runtime of the various methods compare on GPU hardware.
To this end, in Section 2 we introduce the high-order FR approach used in this work. In Section 3, we detail the systems of governing equation for incompressible flow and the artificial compressibility approaches studied in this work. We also explore aspects of the eigen-structure of these equations. Some additional details on the Riemann problem are included in Sections A and B. Then, in Section 4, the main numerical results on unsteady turbulent test cases are presented. Finally, conclusions are drawn in Section 5.

Preliminaries
In this work we consider artificial compressibility approaches solved via the high-order method flux reconstruction (FR) [19], as implemented in the PyFR tool [52]. The original FR method of Huynh [19] has been adapted to handle problems including advection-diffusion equations on element typologies such as simplicies, hypercubes, prisms, and affine pyramids. For completeness the FR method is summarised here, where for brevity we restrict the statement of the method to one dimension. [48,51] and references therein are recommended for applications to alternative topologies.
Characteristic of finite element methods, FR uses a partition of the domain into conformal sub-domains such that = ∪ =1 and ∩ = ∅ if ≠ . Within each sub-domain, a number of nodes are positioned such that a Lagrange finite element can be formed for a conservation equation of the form where in this description of FR we will assume is periodic. It is typical to use a reference domain̂ with the projection ∶ ↦̂ as this makes all the operators the same for a given element topology which has some clear computational benefits. For a line, quadrilateral, or hexahedral element a common choice of reference domain is: For the sub-domain and points ∈ , a Lagrange finite element can be formed giving the solution and flux polynomials as: where is the number of nodal value within each element, and is the j-th Lagrange polynomial defined as for the set of reference points { 1 , … , } ∈̂ . In Eq. (2),• symbolises this is in the reference domain and indicates this corresponds to a piece-wise discontinuous approximation. The flux reconstruction algorithm provides a method to calculate the gradient of the flux function, , corresponding to a  0 approximation of in . This gradient approximation is given as where • and • are the projection to the left and right interfaces. Then is the common interface value formed by using , and −1, and similarly for . Some degree of upwinding should be applied when calculating the common inviscid interface flux in order to stabilise the method. This may be provided from one of a number of approximate Riemann solvers [41]. These common interface values are applied to the element via polynomial correction functions ℎ ( ) and ℎ ( ). To enforce the common value, these functions have the conditions that ℎ (−1) = ℎ (1) = 1 and ℎ (1) = ℎ (−1) = 0.
Once the approximate flux divergence is calculated, one of a number of ODE integration techniques can be used. In this work, Runge-Kutta (RK) time integration will be used for explicit time stepping. When solving the dual time systems that arise in the ACM system, it is logical to use an implicit time scheme. In this case, the BDF2 method is used, coupled to an explicit RK smoother in pseudo-time. Throughout this work, we will use the adaptive low storage RK procedure of Kennedy et al. [21], particularly the RK3(2)4[2R+] method. This is a third order, four stage method where the embedded scheme is second order. The embedded system is used to predict the error [17] and a PI controller can be used to set either a global or local time step [26,51] for the physical time step or pseudo-time step, respectively.

Incompressible flow
The incompressible Navier-Stokes equations for a single phase and constant density can be written as where = [ , , ] is the velocity vector, is the pressure, and is the kinematic viscosity. Taking the divergence of this equation and enforcing a solenoidal velocity field gives the pressure Poisson equation This equation gives a closed form for pressure up to a constant offset. The difficulty arises from the global nature of solutions to the Poisson equation and the complexities of solving this in parallel with the same efficiency as hyperbolic and parabolic equations.

Artificial Compressibility
The artificial compressibility method (ACM) introduced by Chorin [5] can be derived by starting from the compressible Navier-Stokes equations and imposing constant entropy to close an incompressible model [6]. A pseudo-time derivative can be introduced, which interprets a physical time step as converging the field variables in pseudo-time to a steady state. This ACM system of equations can be written in conservative form as: where is the pseudo-time variable, and is the physical time. Here, is a measure of the compressibility of the flow in pseudo-time, defining an artificial Mach number of the flow as = 1∕ 2 , which leads to velocity divergence being balanced by pressure gradients in the pseudo-time dimension. This equation shows the key behaviour of incompressible flow, as → 0 the velocity field converges on to a divergence free solution ∇ ⋅ → 0. Importantly though there is no explicit elliptic equation to solve, and a solver developed for conservation equations may be readily applied to ACM.
To understand the prorogation of information in the system further, consider the flux vector in is defined as = [ , 2 + , , ] , then the inviscid flux Jacobian can be found to be: where = [ , , , ] . This is used to highlight the difference between the conserved variables in pseudo-time and real-time. The eigenvalue of this Jacobian can be found to be: for 2 = 2 + . The repeated eigenvalue indicates that the inviscid system is linearly degenerate, and through calculation of the eigenvectors, it can be seen that this results in a contact discontinuity in the Riemann problem with discontinuities in the tangential velocity components, and . As the contact is not stationary, exact or structure approximating Riemann solvers can be developed without complication via approaches such as that of Elsworth and Toro [12]. However, in this work only Rusanov approximate Riemann solvers were applied to give a common interface flux via: where the Davis type maximum wavespeed estimate is used as: Importantly this shows that the maximum absolute eigenvalue scales with √ , ie the stiffness will scale with one over the fictional Mach number.

Artificial Compressibility -Hyperbolised Diffusion
The hyperbolic diffusion method of Nishikawa [29] aims to remove parabolic terms in governing equations through additional auxiliary equations which, once converged, yield the gradients of the conserved variables. As an example consider the linear advection-diffusion equation: the diffusion term can be hyperbolised by adding an auxiliary equation to give the system: Here pseudo-time is again used to converge the system, where is a preconditioning parameter to account for the differing stiffness of the equations. This system can be understood by considering → 0, as this happens → . As the hyperbolic diffusion method and the artificial compressibility method can both be formulated to use pseudotransient continuation, ACM is a good candidate for hyperbolic diffusion. The new ACM-HD governing conservation equations can then be expressed in three dimensions as: This yields the eigenvalues: where: Consequently, a Davis estimate of the maximum absolute eigenvalue can be taken as: Consideration must be given as to how the preconditioning parameter, , is set. If = ( ) is used then the stiffness will be approximately independent of the viscosity, a similar argument has been proposed by Nishikawa and Liu [30]. The repeated eigenvalues also highlight another feature of this governing equation, namely that there will be a linear degeneracy in the Riemann problem, resulting in a stationary contact discontinuity. This introduces some difficulties in using an exact Riemann solver or an HLLC type approximate solver to produce the common interface flux, as a correctly resolved problem would have a contact discontinuity at the interface. The alternative is to again use a Rusanov type Riemann approximate solver with Davis wavespeed estimate for [ , , , ] . Central differencing was used for the remaining terms owing to the diffusive nature of the auxiliary equations. More details of the Riemann problem relating to ACM-HD are given in Section A.
An advantage of this system should be apparent, we solely have first order derivatives in the system and consequently the numerical stability is only dependent on the advection scheme, which scales with ℎ −1 . This is at the cost of additional equations, but this can be offset by only requiring the simpler advection FR algorithm, as well as an increased rate of convergence reported for hyperbolic diffusion [30].

Entropically Damped Artificial Compressibility
Clausen [6] introduced the entropically damped artificial compressibility (EDAC) method where closure was achieved by minimising density variations, rater than constant entropy as in ACM. This leads to the following evolution equation for the pressure where is a Reynolds number. The key difference is that parabolic regularisation introduced to the mass equation removes the need for pseudo-transient continuation, ie explicit time stepping can be used. The name EDAC then stems from the relationship between pressure diffusion and entropy. If entropy is defined via the functional = log( ), then pressure diffusion terms will provide damping to the entropy field. To see this consider the companion entropy equation which may be calculated explicitly, following the general form of Dafermos [7], for the entropy-flux pair ( , ) = (log( ), ): Not only does this show how pressure diffusion effects the entropy, but it also shows that areas of low pressures will move the solution away from the the physical condition of constant entropy if accompanied by a second derivative of pressure. For example, the canonical case of laminar flow around a cylinder. In order to apply the EDAC method in the PyFR framework, it must be cast as a conservation law. If it is assumed that ∇ ⋅ = 0, then in three-dimensions we obtain for where 1∕ 2 = and 1∕ = . This formulation shows a second aspect of EDAC, that pressure fluctuations are resolved spatially whereas in ACM they could be resolved in pseudo-time.
To understand how the propagation of information differs in this system, we will again consider the inviscid flux Jacobian in , then the eigenvalues can be found to be: for 2 = 2 ∕4 + + . A Davis estimate may again be constructed as: Having defined these speeds, we will mainly use the Rusanov approximate Riemann solver throughout this work for EDAC. However, in Section B we explore the Riemann problem in more detail and define an HLLC approach which will also be investigated numerically. A further insight given from Section B is that the EDAC system cannot support a Riemann problem which leads to a solution containing two rarefaction waves. Although the concept of a shock in an artificial compressibility system is strange, this property can be thought of as consistent with the entropy dissipation of the method. Comparison of these eigenvalues to ACM shows a new pressure dependency which is interpreted as a consequence of resolving velocity divergence as pressure waves in space. Furthermore, assuming that max scales with the true maximum absolute eigenvalue, it can be seen that the EDAC system will often be significantly stiffer and, given that the permitted the use of explicit time stepping, this could lead to smaller stable time steps.

TGV at = 1600
The Taylor-Green [38] vortex is a well studied test case which exhibits three regimes: an inviscid laminar regime, turbulent transition through vortex sheet roll-up, and full homogeneous isotropic turbulent decay. The last regime is generally observed for ≳ 1000 [4]. The initial condition is taken as: = sin cos cos , where a typical Mach number is = 0.08. When performing ACM-HD calculations, these terms are differentiated to give the initial gradients. The domain is a fully periodic cube Ω = [0, 2 ] 3 , which was partitioned into a mesh of regular hexahedral elements with approximately 128 3 solution points depending on order. The solution and flux point locations in the reference domain were set using the tensor product of the Gauss-Legendre quadrature. For this case, unless otherwise stated, the common inviscid flux is calculated using a Rusanov type approach as detailed in Section 3. For the ACM and EDAC methods, the common viscous interface flux in calculated using an LDG approach with a small penalty, = 0.1. Two separate time integration schemes were used. For ACM and ACM-HD, BDF2 was coupled to an adaptive RK34 explicit smoother, details of which can be found in Section 2. This same RK method was used for the explicit time stepping of the EDAC method; however, with global adaptation of the explicit physical step rather than locally for the pseudo stepping. To further accelerate the convergence of the ACM and ACM-HD calculations, Pmultigrid was used [25]. The cycles used are shown in Fig. 1 as they were found to give rapid convergence and follow the suggested asymmetry of Trojak and Witherden [44]. The exact configuration files can be found in the attached electronic supplementary material.
Several functionals were used to access the numerical performance of the AC methods. The first is enstrophy, defined as: where the vorticity is defined as = ∇ × . Reference DNS data for this functional was provided by Van Rees in private communication, containing results for a longer time period than that detailed in van Rees et al. [33]. Two further functionals were examined to assess the quality of the incompressible solution: These are the volume averaged velocity divergence and absolute velocity divergence. In order to evaluate the EDAC method, a sweep was performed for ∈ {3, … , 100} at orders 2, 3, and 4. Studying  Fig. 2, it can be observed that does not have a significant effect on the enstrophy production for the EDAC method at this resolution, and the largest differences are seen at the peak dissipation, shown here. Comparison is made in      3a between EDAC and the other ACM methods. This shows that after peak enstrophy production, there is a noticeable deviation between the results. A further comparison can be made from Fig. 3 between EDAC and the compressible Navier-Stokes equations (NSE) at = 0.08. This makes it clear that EDAC and NSE are comparable, caused by both methods dissipating entropy due to compressibility and both systems having a similar Mach number. In contrast, the ACM and ACM-HD systems do not cause entropy dissipation and the results are evidently different from the entropically damped systems. There is also a notable disparity between the ACM and ACM-HD systems at = 2. The difference is due to Eq. (26) being evaluated with the converged gradient terms in ACM-HD, which are one order higher than the reconstructed gradients used in standard ACM. This difference is more notable at lower as the polynomial space is more restricted.        ACM-HD led to lower velocity divergence than ACM, which is attributable to the improved convergence properties of hyperbolic diffusion methods. Furthermore, lower values of could lead to lower values of velocity divergence, but as increased, this effect diminished. This is indicative of the increased stiffness in the mass equation leading to it becoming dominant in the convergence of the system. An interesting phenomenon is observed in Fig. 5a for EDAC. Here, after peak dissipation, the ratio of absolute divergence to enstrophy becomes approximately linear. Outside of the enstrophy production regime, the source of this can be understood from the vorticity form of the momentum equation. A term that scales with both vorticity and velocity divergence, that would ordinarily cancel, is present in the vorticity equation. The conclusion, which may have been anticipated, is that vortex dominated flows may induce stronger divergence clearing from the method, which may in turn be problematic for EDAC as this manifests as larger pressure fluctuations. At lower order, due to lower enstrophy production and higher numerical dissipation, this relationship is not as clear.
In several places we have considered the maximum eigenvalue due to its importance in the efficacy of the explicit Runge-Kutta time integration. A lower absolute maximum eigenvalue will increase the maximum stable time step, but clearly from Section 3.4 increasing will increase the maximum absolute eigenvalue. Therefore, a trade-off between stiffness, runtime and divergence, exists.   For the EDAC method, a globally adaptive time stepping procedure was used. To understand the trade-off that occurs, the statistics of the time steps used were collected and the average time step size is shown in Fig. 5b. This shows the same linear relationship at high √ values as was observed for the divergence. In the case of Δ , the origin of this relationship can be clearly seen form the eigenvalues of the inviscid flux for the EDAC system of equation. In Figs. 4a and 5b results are also presented for HLL and HLLC as well as Rusanov. The HLL and HLLC schemes are derived in Section B. These data show that HLLC can give a sizeable reduction in the divergence average, and both HLL and HLLC lead to an increase in the maximum usable Δ . The cause of this is a better model of the maximum stable eigenvalue, as the Davis wavespeed estimate more often over predicts this wavespeed [40], which will lead to a stricter CFL condition. Based on these results there is a clear benefit to using the HLLC Riemann solver at the interfaces, and given the availability of FLOPs on GPU, the additional computation can be largely hidden by memory latency caused by bandwidth limitations.

Circular Cylinder at = 3900
The circular cylinder has previously been studied in great detail, a comprehensive review was presented by Williamson [50]. When = 3900, the flow is in a sub-critical regime and is considered to be both computationally challenging and physically interesting. This is due to the presence of multiple phenomena, namely: a laminar boundary-layer, separation, a free shear-layer, turbulent transition of the free shear-layer, and a turbulent wake. The difficulties in simulating this case are exemplified by the spread of reported drag coefficients in the literature [32], with a value in the range of ∈ [1, 1.4] being not uncommon. Subsequently reported by Lehmkuhl et al. [23], a cause for this is that the time averaged wake has two distinct modes, a low energy L-mode and a high energy H-mode, with long non-dimensional times separating transition event between the modes, typically of the order 10 3 . Owing to this long time between transitions, a simulation will typically only capture a single mode.
In the work of Vermeire et al. [46] the PyFR solver was compared to other tools. It was shown that with an appropriate mesh simulations of the compressible Navier-Stokes equations at = 0.2 gave good agreement with the DNS results of Lehmkuhl et al. [23]. For this reason we will use the same mesh in this investigation which covered the domain Ω = [−9 , 25 ] × [−9 , −9 ] × [0, ] for diameter . More recently, Dzanic et al. [10] performed a DNS of this case. Comparison is made later to their data, as well as to additional TKE budget data collected by the present authors rerunning their configuration.  For this test, seven configurations were tested: EDAC with ∈ {4, 20, 100}, ACM with ∈ {3, 4.5} and ACM-HD = 3 with ∈ {10 , 100 }. The case was set with a free stream velocity and pressure of = 0.2 √ and 0 = 1. Constant velocity inlet and constant pressure outlet boundary conditions were used for all cases. Although, Riemann invariant boundary conditions are available, they were deemed unnecessary in this case. An initialisation period was set from 0 tô = ∕ = 47 at = 2 followed by = 4 until̂ = 100. Afterwards, the time average statistics were collected until̂ = 200. To accelerate the convergence of the ACM and ACM-HD dual time stepping, the same P-multigrid method was used as described in Section 4.1.
The instantaneous Q-criterion for the EDAC method at̂ = ∕ = 201 is shown in Fig. 6. Some differences in the wake structure are visible and notably at lower values the free shear layer appears to transition earlier. Focusing on EDAC initially, the time and span averaged velocity profiles -Figs. 7 and 8 -show that there is a pronounced difference between the different values, with the highest value of giving the best agreement with the the DNS H-mode. It is interesting that as is increased monotonic convergence is not observed. From the Q-criterion plots, it seems as though transition is delayed for = 20 compared to = 4; however, small structures in the wake look as though they have coalesced into larger scale structures. Further comparison can be made with additional Q-criterion plots in Section C.
Studying the progression of the streamwise Reynolds stress through the downstream slices, Figs. 7c and 8c, it is observed that EDAC initially overestimates the Reynolds stresses, leading to the downstream stresses at ∕ = 2.02 being under-predicted. Clearly this will have an impact on the turbulent kinetic energy budget, a detailed study of which is presented by Tian and Xiao [39], who showed that at ∕ = 1.06, a balance exists between the convection, production, pressure transport, turbulent transport, and dissipation. Closer still to the cylinder, Tian and Xiao [39] showed the production term is dominant, and we hypothesise that over-production would lead to a downstream energy deficit, and that this is happening for EDAC. The turbulent production is defined as where we have used Einstein notation. Plotting the production term at ∕ = 1.06, Fig. 9a, over-production is clearly visible in comparison to DNS data. This is followed by a drop and subsequent under-prediction of production downstream, Fig. 9c. The question arises, what is driving this over production? From the average pressure profiles of Fig. 10 it can be seen that the pressure in the wake region is significantly lower than the DNS of Dzanic et al. [10]. This is a result of the artificial compressibility in the pressure field. From the interpretation that = 1∕ 2 , it can be understood that this pressure variation will increase as decreases. The pressure gradient this leads to is the driving force of the increased production. Downstream the pressure gradient reduces and there is an associated drop in production. This highlights the difference between the EDAC method and ACM. In ACM, the velocity divergence is balanced by pressure fluctuations resolved in pseudotime, whereas the EDAC method uses the alternative mechanism of spatial pressure variations. The spatial pressure fluctuation are able to have a significant impact on the observed physics, whereas in pseudo-time they are diffused and convected out of the solution by the explicit smoother. A further effect of the pressure reduction is an increase in the entropy dissipation, as seen from Eq. (26), which helps to explain the energy deficit.
Studying the ACM and ACM-HD results, it is interesting to see that from Fig. 7, the ACM simulation with = 3 is in the L shedding mode, whereas all the other results are in the H shedding mode, which acts to highlight the sensitivity  of this case and which mode is initially captured. The ACM-HD results for = 100 outperform those of = 10 , given that a higher reduces the stiffness while in the asymptotic limit a higher is preferred. The ACM results with = 4.5 are comparable to ACM-HD with = 3 and = 100 . However, from the case setup, this ACM configuration was on the limit of what is stable. This can also be understood when considering the results of Fig. 4c.
As was discussed in Section 4.1, the increased stiffness of the EDAC method at higher values will lead to lower maximum stable time steps which will negatively impact the runtime. The data in Table 1 gives a runtime comparison for one flow over diameter of the cylinder, from̂ = 201 tô = 202. This shows that increasing by a factor of 5 gives approximately doubles the runtime of EDAC, which is approximately in line with the increase in the value of predicted from section Section 3.4. Furthermore, the runtimes for EDAC were observed to be lower than those of ACM; however, to achieve comparable results to ACM, for EDAC would need to increased and is unlikely to be faster in that case. For the cylinder case with ACM it was found that Δ = 2 × 10 −3 and Δ ∕Δ = 20 were stable, whereas due to the increased stability of ACM-HD Δ = 4 × 10 −3 and Δ ∕Δ = 20 could be used. This is the cause of reduced runtime of ACM-HD compared to ACM observed in Table 1, and is one of the attractive features of hyperbolic diffusion. Although more equations are solved, the simpler advection algorithm coupled to the stability improvement lead to a speedup even at lower Reynolds numbers [43].   Table 1 Runtime for circular cylinder at = 3900 and = 4 for one flow over diameter and 12 partitions run using NVIDIA V100 GPUs.

SD7003 at
= 6 × 10 4 The SD-7003 aerofoil [36] is an asymmetric aerofoil that has been explored in several studies, for example the studies by Beck et al. [3], Galbraith and Visbal [14], and Garmann et al. [15]. The particular configuration that will be tested here is that of = 6 × 10 4 and an angle of attack of = 8 • . In this configuration, a laminar separation bubble is formed on the suction surface followed by turbulent transition of the boundary layer. From numerical studies [3, 14,   15], it has been shown that the lift and drag coefficient are sensitive to the separation location, which in turn will be dependent on the surface pressure distribution. This configuration was previously used by Loppi et al. [26] to assess the effectiveness of locally adaptive pseudo-time stepping and P-multigrid in the convergence acceleration of ACM within the PyFR framework.
In order to allow direct comparison to there results, the same mesh was used in this investigation. As it was found in Loppi et al. [26]  with the divergence and global force metrics output periodically as well as data on the time step size used. As this case is at a higher , aliasing is of greater concern. To ameliorate this, over-integration anti-aliasing of the flux at the solution points is used [37]. All methods on this case were run at = 4 and the over-integration order was set to = 11, similar to the setup of Loppi et al. [26]. The effect of different terms for the EDAC scheme can be qualitatively understood from the Q-criterion isosurfaces shown in Fig. 11. From these plots, it looks as though the boundary layer is significantly thicker in the = 4 case. There also seems to be some variation in the transition point. Then moving on to study the coefficient of lift and drag, presented in Fig. 12, it is observed that = 4 significantly over-predicts drag and under-predicts lift, consistent with the boundary layer thickening observed in the instantaneous flow fields. At the high values of , good agreement is achieved in comparison to the numerical results of Vermeire et al. [46] and Loppi et al. [26]. Subjectively, the lower = 20 results seems to be more similar to the compressible results of Vermeire et al. [46], whereas = 100 compares more favourably with the ACM results of Loppi et al. [26].
To gain greater insight into the cause of the poor results at = 4 we studied the time and span averaged stream-wise velocity and pressure at the trailing edge, shown in Fig. 14, averaged over̂ ∈ [0, 60]. We were unable to find any detailed DNS data for this particular configuration and location; however, it is clear that there are significant differences  in the physics observed. As was observed in Section 4.2, the spatial pressure variations of EDAC are providing the mechanism driving the inaccuracy. Here, rather than early transition and over-production, Fig. 14a shows that there is recirculating region on the suction surface, indicative of a larger adverse pressure gradient. This recirculating flow is then responsible for a thicker boundary layer, which in turn will lead to higher . As the thicker boundary layer will reduce the flow turning, this will also lower the . As was stated earlier, can have a significant effect on the stiffness of EDAC and Δ and the runtime. The average Δ for the different EDAC configurations are presented in Table 2. Also included in this table is the ratio predicted from using the approximate absolute eigenvalue based on the inflow and outflow. This approximation under-predicts the impact of ; however, it does give a good first approximation if trying to predict simulation cost and is consistent with Fig. 5b.
Similarly to the cylinder test case, this case was run from̂ = 45 tô = 46, and the runtime recorded, similar to the test performed by Loppi et al. [26]; however, here NVIDIA V100 GPUs are used. The results are presented in Table 3 and there is a clear benefit to EDAC over ACM in this case, with ACM-HD being significantly slower than ACM. The reason for this is, although a larger time step can be taken, the anti-aliasing used greatly increases the cost. This is due to the large number of global reads and writes required for the flux calculation. In this case at = 4, = 11 anti-aliasing is used, given the flux for ACM-HD on a hexahedral element requires (1 + + 2 )( + 1) values per element, flux anti-aliasing with a high degree will require substantially more bandwidth than ACM.

Conclusions
In this work we explored the properties of EDAC, the widely used ACM technique, and the hyperbolic-diffusion approach ACM-HD. Novel insight into the Riemann problem for the conservative EDAC system and ACM-HD system was provided and displayed the effect the various parameters of the schemes would have on stiffness. This analysis identified some unique challenges when formulating either maximum absolute eigenvalues estimations or Riemann invariant boundary conditions in EDAC and ACM-HD.
The main advantage of the EDAC system is the ability to use explicit time stepping when simulating the unsteady incompressible Navier-Stokes equations. It was shown in simulations of the Taylor-Green vortex that as the EDAC compressibility parameter, tended to infinity, the divergence scaled linearly with √ . Furthermore, it was found that HLLC common interface calculations could significantly reduce divergence and increase the maximum stable time step compared to Rusanov. In contrast, for ACM and ACM-HD as tended to infinity this linear relation was not observed, instead improvements in divergence diminished. It was shown that in more challenging LES cases EDAC can lead to erroneous physics if incorrectly configured. For the circular cylinder at = 3900, the free shear layer transition could be triggered early by the spatial pressure variation that balances divergence in EDAC. Errors in transitions were also observed with EDAC if poorly configured on the SD7003 aerofoil at = 6 × 10 4 . Again if was too small, the spatial pressure variations could prevent boundary layer reattachment. ACM and ACM-HD simulations on the other hand did not face this problem, with ACM-HD having a comparable or faster run time than EDAC on the cylinder case. The excessive stiffness of the EDAC system required to correctly simulate transition was found to be responsible for long runtime.
It is concluded that EDAC can be an effective alternative to ACM for unsteady problems, although care must be taken in transitional cases. Furthermore, the improved stability and convergence properties of ACM-HD can make more favourable compared to ACM.

A. ACM-HD Riemann Problem
To understand the structure and limitations of the Riemann problem for the ACM-HD system we focus on the 2D system. Calculating the inviscid flux Jacobian in the direction as: This has the eigenvalues: for 2 = 2 ∕4 + ∕ and 2 = 2 + + ∕ . The associated eigenvectors are where = √ 2 + . This shows that the linear degeneracy indicated by the repeated eigenvalues 1 = 2 = 3 , results in a stationary contact discontinuity. The discontinuous properties are , , as well as and . These final two are part of the same wave and so are related through a Riemann invariant. A stationary contact discontinuity makes formulating an exact Riemann solver more complex, although possible. If the purpose of the exact Riemann solver is to form an upper bound on the maximum absolute eigenvalue [16], then this is not an issue; however, it is not possible to produce a common interface value. A similar issue would be confronted by the HLLC method and a linearisation of the central contact would be necessary. The structure of the 2D ACM-HD Riemann problem is shown in Fig. 15, given that:

B. EDAC Riemann Problem
To understand the structure of the wave fan for EDAC and devise Riemann solvers, we will first state the Riemann problem: Using a transformation to the reference problem̂ = [1, 0, 0] , we can form general solutions to the Riemann problem by only studying the first column vector for the invsicid flux, = inv ⋅̂ . The Jacobian of the conservative formulation of EDAC can be straightforwardly found as for 2 = 2 ∕4 + + . The associated right eigenvectors are: This combination of eigenvalues and eigenvectors tells use that the Riemann problem for this case does have a linear degeneracy, in that it can support a contact discontinuity in and , the structure can be seen more clearly in Fig. 16. With this established, the strategy to form an exact Riemann solver will follow a similar procedure to Elsworth and Toro [12]. We will include the main steps of the procedure as there are some interesting differences compared to the standard ACM formulation.

B.1. Exact Solution
The first stage of the exact Riemann solver is to find the non-linear equation that governs * , this has the general form where ( * ) = * − , and ( * ) = − * . Solving this will give * and * . The -waves 3 and 4 can be either rarefactions or shock waves and will depend on the wave type. We will start by finding for the case when the waves are rarefactions.    Figure 18: Centreline average quantities, with experimental data from Parnaudeau et al. [31] and DNS data from Lehmkuhl et al. [23].