Towards robust unstructured turbomachinery large eddy simulation

Abstract Industrial legacy codes usually have had long pedigrees within companies, and are deeply embedded into design processes. As the affordability and availability of computing power has increased, these codes have found themselves pushed into service as large eddy simulation solvers. The approximate Riemann solver of Roe, which is frequently used as the core method in such legacy codes, is shown to need much user care when adopted as the discretisation scheme for large eddy simulation. A kinetic energy preserving (KEP) scheme—which retains the same advantageous stencil and communications halo as the original Roe scheme—is instead implemented and tested. The adaptations of code required to switch between the two schemes were found to be extremely straightforward. As the KEP scheme intrinsically bounds the growth of the kinetic energy, it is significantly more stable than the classical non-dissipative schemes. This means that the expensive smoothing terms of the Roe scheme are not always necessary. Instead, an explicit subgrid scale turbulence model can be sensibly applied. As such, a range of mixed linear–non-linear turbulence models are tested. The performance of the KEP scheme is then tested against that of the Roe for canonical flows and engine-realistic turbine blade cutback trailing edge cases. The new KEP scheme is found to perform better than the original in all cases. A range of mesh topologies: hexahedral; prismatic; and tetrahedral; are also tested with both schemes, and the KEP scheme is again found to perform significantly better on all mesh types for these flows.


Introduction
In many industrial codes, the emphasis is placed on the ability to conform to complex geometries. This necessitates the use of unstructured formulations, and frequently results in second order algorithms being employed. Many of these industrial codes have long pedigrees stretching back many years, and were originally designed to perform Reynolds-Averaged Navier-Stokes (RANS) calculations. The extreme advances in computational power and corresponding decrease in hardware costs have brought the significantly more computationally intensive-but higher fidelity-Large Eddy Simulation (LES) and hybrid RANS/LES methods in from the realm of academic and intellectual investigation into the purview of industrial design calculations.
From the point of view of both the end-user and the support team, it is attractive to reuse as much of the RANS code as possible in constructing an LES solver, as much of the structural framework already exists, requiring little or no adaptation to file formats, boundary conditions, parallelisation libraries, and substantial code optimisation has already been conducted. To take advantage of this, an industrial legacy RANS code is here applied to LES simulations, its performance assessed, and any adaptations or improvements which are found to be necessary are made.
The code used in this exercise is HYDRA, a staple of the Rolls-Royce design process for many years. Variants of the scheme used in HYDRA have, in the past, been used to successfully perform a range of large eddy simulation calculations. Typically, these have either involved relatively high Mach number-very compressibleflows, or codes which have been heavily, and often on a problem-by-problem basis, modified. Ciardi et al. [1] developed a low dissipation version of the Roe scheme based on suppressing the appearance of dispersive ''wiggles'' in the local flow solution. A set of user-supplied constants were used to control the trade-off between minimising excessive dissipation and ensuring solution stability.
Other attempts to reduce any inherent dissipation have been made over the years. Page and McGuirk [19] and Li et al. [11] chose to calculate and use a Ducros et al. switch [3] on the smoothing terms to try and keep the diffusion under control-the subgrid effects themselves were then separately modelled with more explicitly added terms. O'Mahoney et al. [18] were also compelled to reduce the strength of the stabilising terms to avoid the http suppression of turbulence when modelling ingestion by rim seals. The dangers of excessive smoothing have been extensively discussed by Tristanto et al. [23].
Examples of some successful simulations in more specific applications include the transonic Mach number jets flows of Eastwood et al. [5], and the chevroned nozzles of Xia and Tucker [25]. In both of these cases, the smoothing constant was geometrically sculpted a priori, to act as a numerical turbulence model in areas of interest, and to provide stability in the far field. This has proved a successful approach, as long as the general nature of the flow is understood before the sculpting is carried out.
Despite these successes and advances, the results from the use of the Roe scheme as a discretisation for an LES solver for wider problems-or problems involving a range of flow conditions-have been disappointing at times, particularly when flows have been dominated by regions with low Mach number. In this paper, it is hoped that the problems HYDRA encounters as a legacy RANS solver applied to LES calculations can be explained, quantified, and mitigated as painlessly as possible.
The Euler equations are given by: This equation can be discretised in many ways, but for these purposes, discretisations are restricted to two-point edge based schemes, for ease of development, efficiency of parallelisation, and simplicity of application to unstructured problems. Here, the established Roe scheme is compared to a kinetic energy preserving (KEP) ''Jameson'' formulation. To solve the full Navier-Stokes equations, Eq. (1) must be augmented with a viscous flux term, F V . When these equations have been filtered onto a grid, residual stress tensors emerge, necessitating a model for the subgrid scale effects of turbulence. This paper first discusses the two different numerical approaches, after which the various test cases are introduced, before the numerical results are presented, before finally discussing what these findings suggest for industrial large eddy simulation.

Numerical methods
The solver used in this work, in its original RANS formulation, is built around the approximate Riemann solver of Roe [21]. In smooth regions of the flow, away from sharp gradients such as shocks, the inviscid flux calculation takes on the form of a central difference of the end points of each node, smoothed by some function of the Laplacian of the conserved variable vector.
Near to sharp gradients, the inviscid flux calculation decomposes to first order, and takes the form of a central difference smoothed by some function of the conserved variable vector itself.
The weighting of these two functions is given by a local pressure switch, intended to ensure the full second order accuracy of the former in areas of smooth flow, and to allow this accuracy to degrade for stability in the region of strong gradients. The matrix A ij is the inviscid flux Jacobian, A ij ¼ @F I =@Q . The variables e 2 and e 3 represent user supplied constants which control the smoothing levels. These variables, which are constants in the classical Roe scheme, were modified by Eastwood [4] and Ciardi et al. [1] to vary in space, and to vary in both time and space, respectively. An explicit 5-stage Runge-Kutta (RK) method is used to integrate the Roe scheme in time.
The unsmoothed fluxes themselves are given by the average of the flux each side of the control volume face.
This scheme brings a number of advantages to the RANS solution process. It offers considerable stability, whilst retaining second order accuracy in smooth regions. On industrial grids, a Courant number of 2.00 is found to be sufficient for stability (but not, for LES, for accuracy in time), even with the smoothing and viscous flux calculations updated only every other RK sub-step. This leads to a considerable acceleration in steady convergence, whilst maintaining second order accuracy away from shocks.

Conservation of kinetic energy
It is understood that as a bounding quantity of incompressible flow, the global growth rate of kinetic energy will dictate the stability of the scheme-with a global kinetic energy growth rate of zero, the scheme will be stable. For compressible flows, conservation of kinetic energy does not formally guarantee the stability of a solution, but evidence suggests that it does impart some resilience to schemes. Global kinetic energy is given by: A transport equation can then be written for the kinetic energy, k. @k @t ¼ u @ðquÞ @t À u 2 2 Following Jameson's one dimensional analysis through, a criterion is found for preserving kinetic energy [7]. Assuming boundary conditions are set appropriately, the integral in Eq. (5) can be written in discretised form as: This is satisfied if: As Jameson suggests, for a second order scheme, this can be satisfied if, as specified, we set: To implement this change within a Roe solver, it is simply necessary to change the order of the averaging of the inviscid flux terms. Instead of multiplying the primitive variables at each end of the edge to give the fluxes and then averaging them, the primitive variables are averaged at each end first, and then multiplied. This gives: The stability increase that the enforcing of the conservation of kinetic energy supplies means no artificial dissipation terms are required, so these can be removed:

Operator skew self-symmetry
By considering the one dimensional Burgers' equation, it is possible to explore the relevance of skew-symmetric operators to kinetic energy preserving schemes [20]. The full, viscous, Burgers' equation is: Global kinetic energy is conserved under the action of convection, but it is dissipated by viscous action in both the viscous Burgers' and the Navier-Stokes equations. Therefore, here the conservation criterion is applied to solely the inviscid parts. For clarity, this leads to considering the inviscid Burgers' equation, which can be written in both convection (Eq. (15)) and divergence (Eq. (14)) form: @u @t þ @ @x Combining these two forms can give the skew-symmetric form of the Burgers' equation: where B C u denotes a continuous ''Burgers' operator'', acting on u. Discretising, and making use of matrix-vector notation, where u represents the solution vector, and D the transport discretisation operator, gives: The global kinetic energy is then given by: Differentiating, the rate of change of global kinetic energy is: So, for non-trivial cases, the discrete global kinetic energy is conserved if and only if D ¼ ÀD Ã , i.e., if the transport discretisation operator is skew self-symmetric. This symmetry is a property of the continuous equations, so should be reflected in the discretisation scheme. Denoting the central difference and the KEP scheme flux transport operators by D R and D K respectively, the form of the two can be compared: It can be seen that the transport discretisation operator for the central difference scheme, D R , does not form a skew-symmetric operator. In contrast, the KEP scheme transport discretisation operator, D K , is skew self-symmetric. This failure of global conservation on behalf of the CD scheme requires the introduction of the smoothing terms which make up the full Roe scheme. This deviation of the Roe scheme from the physical continuous equations means that Roe transport terms a-physically influence the quantity of global kinetic energy. The KEP scheme, on the other hand, is capable of correctly capturing the skew-symmetric property of the continuous equations.

Subgrid scale modelling
The Roe scheme, when it is applied to LES calculations, is often run without an explicit turbulence model. The artificial diffusion terms are, in these cases, considered to provide sufficient dissipation to act as a numerical subgrid scale (SGS) model. The dissipation-free KEP scheme, however, requires the use of explicit SGS terms to account for the effects of the flow below the grid scale. A number of these SGS models have therefore been tested. By making use of the formulations of Liu et al. [12], a set of mixed linear-non-linear turbulence models have been introduced. These mixed models are somewhat more advanced than the simple static Smagorinksy, being able to capture some of the effects of the reverse energy cascade, but without being as expensive as the dual-filtering operation required for dynamic turbulence models, such as the Germano et al. [6]. Again, introducing these models when developing from a Roe solver is extremely straightforward-all the gradient terms required by the turbulence models are already available in the relevant subroutines, having been previously evaluated for the smoothing. The mixed models themselves are formed by the division of the turbulence terms into a purely dissipative linear eddy-viscosity term-here captured by a Smagorinsky model-and a non-linear term. The subgrid stress tensor, s R ij , can thus be written: The models which were used in this study are listed by their constituent terms in Table 1.

Isotropic decaying turbulence
The classic LES test case is that of isotropic decaying turbulence. Here, the experiments of Comte-Bellot and Corrsin [2] were used to provide both the initial conditions and the data for final comparison. Two fully orthogonal hexahedral meshes were prepared with triply-periodic spatial boundary conditions, one of 32 nodes along each side, the second with 64 nodes along each.
Secondly, a set of meshes were generated with different topologies-each being constructed from different volume elements, hexahedrals, triangular prisms, and tetrahedrals. Each of these meshes was designed to have approximately the same edge count (rough computational cost) as the 64 3 node simulation.
The simulation was initialised at time t 0 ¼ 0, which corresponded to the experimental data at non-dimensional measurement location tUo M ¼ 42, from the 5.08 cm grid experiment. The calculation could then be advanced through time to correspond with the experimental data at measurement location tUo M ¼ 171, at which time the resulting three-dimensional energy spectrum, EðkÞ, was extracted for comparison. The cost to carry out a simulation on the 32 3 case was around 5 core minutes on a single core of an AMD Phenom II 2.8 GHz processor. For the 64 3 , this cost was around 40 core minutes on the same architecture.

Imperial College case
The Imperial College case was taken from a series of experiments carried out at Imperial College in the late 1960s and early 1970s. The full geometry description is given in Kacker and Whitelaw [10], although a large series of experiments were carried out [8][9][10]. This represents a more challenging and industrially relevant case than the isotropic decaying turbulence. It comprises a pair of co-flowing jets, one of which is bound by a wall. The two jets are separated by a narrow lip. As the shear layer grows, it eventually impinges on the wall, resulting in a wall-jet type of flow.
This case was simulated at a high velocity ratio of 2.32. Although the mean profiles at inlet were supplied to the flow, the inflow was set as laminar-a deviation from the experimentally measured turbulence values of around 0.5%. This was not considered to be an important difference in this case, as some turbulent inflows had previously been tested on this wall-jet with a structured high order LES solver [17], including a Lund recycled inflow [13], and box turbulence imposed onto a mean profile. These imposed turbulent boundary conditions were found to have little influence on the prediction of either the mean flow or the turbulent statistics.
A triangular prismatic mesh was generated to represent this geometry, consisting of around 20 million nodes. The prisms were aligned to run in the spanwise direction. The mesh was generated to conform to the near wall distances in Table 2. Close to the lip, there were 35 nodes across the shear layer.
The large size of the domain (one metre long), and the relatively high element count necessitated very long run times to both mature and to average the flow. This took around 80,000 core hours on AMD Interlagos Opteron 2.3 GHz processors. The averaging time was set as the time required for the mainstream flow to travel two hundred slot heights downstream.

Karlsruhe case
The Karlsruhe case was taken from a set of work carried out in Karlsruhe University on engine-realistic geometries for cutback trailing edge film cooling [15,14,16]. This case represents an industrially relevant, and geometrically and aerodynamically complex, problem. The case consists of a basic geometry, with a coolant injection angle of 10°. The geometry specifications are shown below in Fig. 1, and given fully in [16].
These experiments varied both the turbulator layout within the coolant cavity, and the blowing ratio of the coolant and mainstream flow. The domain was again meshed using prismatic elements aligned in the spanwise direction, with a node count of 19 million. The mesh was generated to conform to the near wall distances in Table 2. Near the lip, there were 55 nodes across the lip thickness. The boundary conditions used in this solution were simple planar laminar conditions, as it was assumed that the turbulators would overwhelm inflow turbulence conditions.
To run a single simulation of the Karlsruhe geometry took around 20,000 core hours on AMD Interlagos Opteron 2.3 GHz processors. The averaging time once initialisation transients had been convected through was set as the time taken for the mainstream to travel one hundred slot heights downstream.

Isotropic decaying turbulence
To assess the difficulties faced by the Roe scheme when applied to LES calculations, it was initially tested upon the isotropic

Alpha
À2q decaying turbulence case. The performance was first tested by the simulation of the Comte-Bellot and Corrsin [2] experiments on triply periodic regular hexahedral grids. These were carried out without explicit SGS turbulence models, to determine the ability of the Roe solver to act as a numerical LES algorithm. The results of these simulations are shown in Fig. 2. It is apparent that the scheme tends to suppress excessive amounts of resolvable kinetic energy at higher wavenumbers. This fundamentally influences the nature of the turbulent cascade, with energy ''piling up'' at low to moderate wavenumbers. This remains true even with the quite fine 64 3 node grid, although as more of the energy is resolved, the wavelength of the first over-damped wavenumber is shortened. To more fully test the scheme, its ability to cope with more skewed grids was considered. As one of the reasons for making use of unstructured codes like HYDRA is their ability to cope with complex geometries, there is little point if their performance on non-perfectly orthogonal hexahedral meshes is poor. To this end, the simulations were repeated on the hexahedral, triangular prismatic, and tetrahedral meshes. The ability of the code to capture the evolution of the turbulent energy spectrum as a function of grid topology could then be determined. These results can be seen in Fig. 3. The hexahedral meshes, although some way from the experimentally determined values, capture the spectra most accurately. As the orthogonality of the topology deteriorated, the performance of the solver fell. The prismatic mesh is capable of performing only slightly more poorly than the hexahedral. However, the tetrahedral grid-despite its isotropy and near-regu larity-performed remarkably poorly, suppressing the resolution of turbulence at even moderate frequencies. This is worrying, as tetrahedral meshes are by far the easiest to rapidly fit around complex geometries. It is theorised that this dissipation factor is a function of the orthogonality of the mesh and its median control volume dual. The hexahedral mesh and its dual are perfectly orthogonal, while the prism mesh and its are very nearly so. However, it is geometrically impossible for tetrahedrals to behave in the same way: they are very in-orthogonal with their duals, hence the extremely high levels of dissipation experienced with the tetrahedral mesh. This may go some way to explaining the received wisdom within the CFD community that tetrahedral meshes are generally inappropriate for simulations, despite the ease of their generation.
We turn now to the KEP scheme. Although it is now stable without smoothing terms being included, it no longer has artificial dissipation to act as a numerical turbulence model. Initially, the Smagorinsky dissipative terms [22] are explicitly included within the viscous flux calculations. Again, the performance of the new scheme on the 32 3 and 64 3 orthogonal hexahedral meshes is compared to the Comte-Bellot and Corrsin experimental data. The results of the new scheme are more encouraging than those of the original solver, with both grids performing well-highlighting the ability of the new scheme and turbulence model to capture the behaviour of the turbulent cascade well, even on coarse meshes. These new HDT spectra are shown in Fig. 4. The moderate wavenumber excess and high wavenumber deficit have been mostly eliminated.   To assess the relative performance of the new scheme against the original on the more industrially relevant meshes, the different topology grids were again used. Once again, a considerable advantage was found when using the new method. The relevant spectra for the three topologies are shown in Fig. 5. Much the same trends are shown with the new scheme as with the original-tetrahedral meshes are found to be significantly more dissipative than their hexahedral counterparts, with the prism elements somewhere in between. It is, however, worth noting that the worst performing (tetrahedral) grid-when using the KEP scheme-was in fact less dissipative than the Roe scheme on the hexahedral grid. This comparison is shown in Fig. 6. It is possible that making use of this scheme could allow the more widespread use of tetrahedral grids, which are more easily automatically mapped to highly complicated geometries. This discussion, however, is left to a later date.
These basic studies make it apparent that the KEP scheme, with an appropriate turbulence model, is capable of significantly outperforming the original Roe scheme.
More advanced SGS models were then applied to the isotropic decaying turbulence test case. Their relative performance on the 64 3 node grid is shown in Fig. 7. As expected, the inclusion of the non-linear terms does not significantly influence the performance of the basic Smagorinsky model in this case. Indeed, these models are-by design-intended to have little influence here, as the Smagorinsky model is specifically calibrated for this basic test case. These models are intended to introduce more realism as the nature of the physical turbulence moves away from the simple HDT case, improving the prediction of subgrid turbulent behaviour in more complex situations, such as free shear flow, wall bound flows, and the like. To perform the more complex industrial-type calculations, the LANS-a mixed model was selected.

Imperial College case
The ability of the new code to predict the turbulent behaviour of somewhat more challenging flows was then tested by the simulation of some flows from the Imperial College experiments. These flows are characterised by a shear layer developing between two co-flowing jets. One jet is bounded by a wall, onto which the shear layer eventually grows to impinge. The resulting flow could be considered as something of a cross between a boundary layer and a free shear flow. A snapshot of the developing turbulent shear layer between the co-flowing jets is shown in Fig. 8.
Even simple wall-jets, such as these, have, in the intervening years since the experiments were conducted, proved remarkably intractable to RANS simulation. Liu showed a very large degree of scatter in results when varying the RANS turbulence model. Both the original scheme and the new KEP scheme are tested and    compared to the experimental data, to assess the performance gain on a more challenging LES flow than the simple isotropic decaying turbulence. Comparing the predictions of the time averaged x-direction velocity components at 10 slot-heights downstream of the ejection plane, the mean flow fields are not substantially different. It appears from this as though the Roe scheme has more thoroughly mixed the jets together. This was the effect of the smoothing-at this point in the flow, the smearing effect of the artificial viscosity terms of the Roe scheme was significantly greater than the turbulent mixing of the shear layer predicted by the KEP. That this is the case is borne out by an examination of the turbulent statistics at this point, Fig. 9. Significantly less u 0 u 0 turbulence is found to have been generated by the Roe scheme than by the KEP scheme. The statistics for the v 0 v 0 turbulence terms are even more striking-practically none is captured by the Roe scheme compared to that by the KEP. The discrepancies between the experimentally collected turbulent statistics and the KEP predictions are believed to be related to the laminar inflow conditions, which delay the development of a fully turbulent shear layer. In all cases, however, the KEP scheme dramatically outperforms the Roe. A qualitative examination of the flow predictions shows almost no turbulent mixing between the two layers when using a Roe discretisation, with the coolant wall-jet remaining very much intact without the development of large scale turbulent motions to break it up.

Karlsruhe University case
To assess the ability of the codes to perform fully industrially relevant simulations, calculations on the engine-realistic film cooling experiments of Martini and Schulz et al. were carried out. RANS calculations on these geometries are known to be potentially both unrealistic and heavily model-dependent. The main driving mechanism of these trailing edge film cooling flows is the shedding of vortices from the lip of the cutback-RANS methods are unable to comprehend this fundamentally unsteady and large-scale coherent phenomenon. This family of flows, then, is an area of active research, being difficult to simulate computationally, and having industrial relevance to engine design.
Initially, the original Roe discretisation scheme was used to perform an LES calculation on the 19 million node prism grid. The unsteady interactions between the coolant flow through the turbulators and the mainstream flow were visibly weak, and appeared qualitatively non-physical, even on such a relatively fine mesh- Fig. 10. In some of the more benign cases, it was sometimes difficult to get any turbulence to form in the shear layer at all. By lowering the smoothing constants, it was possible to allow some turbulence to form within the shear layer, but significantly lowering the value of e 2 caused increasing problems with stabilitywhich could prove very expensive, wasting thousands of core hours of computational effort.
It is likely that continued refinement of the mesh would eventually produce appropriate solutions, as the truncation error terms approach zero. Unfortunately, in an industrial context, a 19 million node LES calculation is already an expensive undertaking. If, by a small change to the discretisation scheme, acceptable results can be derived without resorting to quasi-DNS meshes, this will represent a significant step forward with the use of LES as a potential design tool.
In the Roe-scheme based simulations which failed to generate any shear layer turbulence at all, the situation is problematicalthough these simulations were LES-cost, they produced RANS-accuracy; a poor return for what represents a significant computational investment. The LES equations collapse onto the RANS equations in the absence of any time-dependency, with a mixing length model with an non-physical length scale coming from the artificial diffusion terms-effectively a mixing length model based on grid spacing. This leads to extreme grid-dependence of solution, and the likelihood of obtaining any useful results is slim. In some cases, these LES calculations offered worse performance than that of some RANS turbulence models, at hugely more expense- Fig. 11.
In contrast, the KEP scheme's performance was much more positive. Qualitatively, turbulent fluctuations and the expected large scale coherent von Karman vortex street can be clearly seen- Fig. 12. Much more turbulent mixing is observed between the Fig. 8. Wall-jet shear layer development. Fig. 9. Imperial College case turbulence profiles, x=y c ¼ 10.
two streams when using the KEP scheme than in the original formulation. Comparing the turbulent statistics of the flow, it is clear that the Roe scheme is, as in the isotropic decaying turbulence test case, suppressing the development of turbulence considerably.
The Roe scheme did not capture any of the experimental behaviour of the film cooling distribution, as the flow behaviour is not physically represented by the solution. However, the new KEP scheme, on the same grid, at the same running conditions, performed significantly better-providing some predictions which better agree with the experimental results. A comparison of the experimental and the two contrasting schemes can be seen in Fig. 13.
A full set of simulations were then carried out on the Karlsruhe problems, to assess the relative ability of the KEP code to predict effects on the adiabatic wall film cooling effectiveness of both blowing ratio and planform geometry. The results for the planform variation are shown in Fig. 14, and those for the blowing ratio variation in Fig. 15. The new KEP scheme performs remarkably wellfar superior to anything RANS has been shown to achieve. Noticeable, however, is the common offset-mixing was systematically predicted to be very slightly greater than is seen in the experimental data. However, the trends were captured remarkably accurately.

Discussion
It is clear that for these two co-flowing jet cases, as well as for the more general isotropic decaying turbulence test case, the KEP scheme seems to perform much better than the Roe when being run in ''LES-mode''. However, it should be pointed out that these cases were chosen specifically for study because the Roe scheme was found to capture them particularly poorly.
The promise that this scheme has shown is high-potentially even opening up some more rapid mesh generation techniques that have previously been considered inappropriate-particularly the use of tetrahedral elements to allow more efficient use of grid nodes, and more complex geometries to be quickly discretised, although this needs substantially more thorough studies.
The ease of implementation of the new scheme is also considered to be attractive. In terms of coding, at least, most legacy URANS codes would require very little adaptation to make use of     this discretisation scheme rather than the common Roe solver. In terms of both effectiveness and efficiency when applied to more realistic cases, the new scheme appears to offer considerable advantages. It has already been noted that it is likely that the original scheme could predict the flow behaviour for these cutback trailing edges with a significantly finer grid. Without carrying out extensive and expensive studies, it is difficult to know how much grid refinement would be required for this. Based on the results for the isotropic decaying turbulence, however, it is likely that the minimum number of nodes to produce a valid Roe scheme solution would be very considerably more than the minimum number required for a valid KEP solution. It appears that this saving could be of an order of magnitude, or more.
Within the current coding framework, there is a small performance benefit to running the new scheme. The very expensive calculations for the pressure switch need not be calculated, and the smoothing matrix A ¼ LKL À1 is not used. This also means that the Laplacian terms do not have to be determined and carried around. On the other hand, the KEP scheme requires calculation of the residual stress tensor, s R ij . Although the new code has not been through the extensive streamlining of the original, there is an approximately 20% speed benefit to the new scheme (architecture and problem dependent).
One of the real benefits of the KEP scheme is that the behaviours of the dissipation terms are now known: the solution quality is no longer reliant on a user-specified smoothing constant that can also influence stability. Instead, explicit subgrid scale turbulence models can be used, without having other dissipative terms acting in concert and potentially compromising the ability of the code to accurately capture the influence of the subgrid scale fluctuations.
We have focussed here on flows with relatively low Mach number. As a note of caution, one of the strengths of the Roe scheme is its ability to handle high Mach flows, and, with an appropriate pressure switch, to be able to capture shock waves without excessive dispersion error. It is unlikely that the kinetic energy preserving scheme-which has been shown here to outperform Roe at low Mach-will be able sustain its improved stability at high Mach numbers unaided. As a result, a blend of the two schemes, based on local Mach number, is suggested, to favour the KEP discretisation at low Mach, and the Roe in trans-and super-sonic regions.
Schemes such as the one tested here have been criticised for use in solving the full viscous Navier-Stokes equations [24]. The high gradients which are induced near viscous wall boundaries can generate substantial numerical dispersion in these unsmoothed equations. In the problems tested here, we have not found this to adversely affect the stability, or the accuracy of the solutions. It is believed that the latter is due to the tested flows being fundamentally dominated by large scale geometrically generated unsteadiness coming from upstream wakes and Kelvin-Helmholtz rollers, hence the near wall performance of the scheme is not paramount. It can be seen from Fig. 9 that in this context, the scheme is adequate for the uses shown here.
Furthermore, from an industrial perspective, the types of high Reynolds number flows where this may be a problem for accuracy and stability are highly unlikely to be directly simulated using wall resolved LES. Instead, it is likely a hybrid RANS-LES scheme would be employed, in which the inviscid flux calculation could also be changed in the near wall regions to, say, Roe, which is more suited to RANS calculations.

Conclusions
Jameson's KEP scheme has been shown to offer a number of significant advantages over the common Roe approximate Riemann solver scheme for use with large eddy simulations, particularly at low Mach number.
These advantages have been demonstrated for both test cases and industrially relevant cases. The reduced dissipation of the KEP scheme does not excessively suppress the development of turbulence in the same way as the Roe scheme on equal grids. This allows the use of much coarser grids than would be required to achieve the same level of accuracy with the Roe solver. Equally, the elimination of inherent dissipation in the solution means that an advanced turbulence model can be applied, likely to be better able to model the effects of subgrid scale fluctuations than that of the artificial dissipation terms originally designed for stability.
Importantly, the ease of switching between the two in an industrial unstructured code makes the decision to employ the KEP scheme inexpensive. The returns in performance-both in terms of accuracy and speed-for the very small outlay in effort can be exceptional.
It is hoped that this ease of development will go some way towards ameliorating the discomfort of the continuing and increasing adoption of LES as a design tool, which is to some degree based around the perceived complexity of developing, maintaining, and supporting a new LES code. This has been shown to be unnecessary, as existing legacy codes can perform well as LES simulators with only minimal modification to their existing structures.
The issue of cost still hampers the acceptance of LES as an industrial tool; at the moment it is still important to pick and choose where to employ it industrially for maximal returns in flow understanding and performance. This is likely to remain the case for the immediate future.