Turbomachinery simulation challenges and the future

Dramatic changes to aircraft design are on the horizon to meet ambitious efficiency and emissions targets. This will bring about changes in the role and design of aero-engines, requiring a greater degree of high fidelity, coupled modelling. To bridge a wide range of spatial and temporal scales, modern turbomachinery design employs a wide range of modelling fidelities during different design stages and for different engine components. High fidelity methods such as Large Eddy Simulation (LES) have been successfully applied to numerous complex flows where traditional Reynolds-Averaged Navier-Stokes (RANS) modelling does not sufficiently represent flow physics to be consistently accurate. Using LES and, for high Reynolds number flows, hybrid LES-RANS, critical knowledge has been extracted and exploited both to inform designs and to improve lower order modelling. It seems certain that uptake of high fidelity methods will continue at a rapid pace. Here we explore some of the future challenges facing turbomachinery simulation using LES. These include use of higher order schemes, internal and external zonalisation and coupling, exploitation of hardware and pre and post-processing. Although these pose technological barriers, these will enable unexplored design space to be traversed with confidence, resulting in cleaner, quieter, and more efficient aircraft.


Introduction
The challenge to reduce pollutants such as CO 2 , NOx and noise in aviation continues to be a key technology driver. Aside from the airframe, civil turbofan aeroengines have been the focus of intense research, to reduce these emissions. The Advisory Council for Aeronautical Research in Europe (ACARE) set the target of a reduction in perceived noise level emission of 50% (to those of the year 2000) by the year 2020 [1]. The EU Flightpath 2050 report [2] states a 75% reduction in CO 2 , a 90% reduction in NOx and a 65% perceived noise level reduction relative to year 2000 levels as a goal. Improving recent aero-engine technologies has enabled significant reductions to be made, largely through the use of CFD and supporting experimental data.
A wide range of future aircraft designs are under consideration for the 2030-2035 timeframe, including: slender winged airframes with trusses, blended wing-flaps with morphing capability, blended wingbodies with podded or embedded boundary layer ingesting engines and subsonic and supersonic jets. Fig. 1 shows some concepts focused on reducing fuel burn and pollutants including noise. Most will include some form of gas turbine power system as direct propulsion or as a means to generate electricity for propulsion or auxiliary systems. Hall and Crichton [3] present a propulsion system consisting of an ultra-high bypass ratio turbofan driving two additional fans, embedded within an S-duct, for an all-lifting body airframe. As part of the Silent Aircraft project, fuel burn is predicted to be up to 45% less than a Boeing 777 [9]. Hybrid gas-turbine and electric propulsion systems coupled with airframes including large laminar flow regions are also under consideration by NASA and Boeing. These provide estimates of reductions in CO 2 (69-81%), NOx (21-42%) and noise (16-37 dB) in the 2030-2035 timeframe [10]. Many enabling technologies are presented by Ashcraft et al. [4].
Small improvements in efficiency from each component broadly add up for the whole engine. The number of passenger and cargo flights per year is forecast to increase over 20 years at 5% and 5.8% respectively (from 2011 as [5]), magnifying this effect over decades of service. However, the role that the entire engine, as well as the relative impact of the core, bypass or potentially open or closed rotor fans will change. For a modern bypass ratio the core accounts O 1/10 th the total thrust. It is clear that given the losses generated in such a complex flow path, the improvements made to the core flow will have a reduced relative contribution to the overall aircraft efficiency, yet must still efficiently drive ever larger fans. Large components such as the fan and installation effects however, now have a significant impact on overall engine and aircraft performance. New low speed fans may be geared rather than directly driven by the LPT. Great efforts are being made to improve low speed fan performance, requiring new understanding of flow, structural and acoustic interactions. The complex 3D shape of modern fans indicates multi-fidelity, multi-physics numerical simulation will be heavily relied upon. Therefore, coupling at the internal (core), engine (fan-intake) and external (engine-airframe) scales will later be considered. Similarly, critical core areas with important coupling such as heat transfer in the turbine will be discussed.
For context, we look back at how we arrived at today's status of turbomachinery CFD. Early use of velocity triangles in the 1950s proceeded to throughflow models such as Wu [11] and gradually led to 2D and 3D Euler and Navier-Stokes solvers. The main move from thousands to millions of cells was made in the mid 1980s to the mid 1990s using steady RANS, enabled by consistently improving CPU speeds. It could be argued, this is where reliance on advances in hardware to improve solutions was imprinted within the minds of a CFD generation. This has lead to poor exploitation of available computational resources today and the delay of algorithmic advances [12]. In the mid-late 1990s, unsteady CFD saw use, but turbulence modelling had become limiting, unable to model important flow physics such as secondary flows, transition and separation. In addition, the importance of coupling of external and internal flows became apparent. In Abhari and Epstein [13], hot gas was found to enter turbine blade cooling holes due to high external pressure fluctuations, underscoring the physical unsteadiness and coupled nature of such flows which was not captured by numerical models. Denton and Dawes [14] summarise turbomchinery CFD practice in the late 1990s by which time 3D multi-stage calculations were possible, albeit with the limitations of mixing planes (inter-row circumferential averaging) and approximations including interactions with secondary gas systems and leakage paths. Looking back these were not well modelled, the reality being that they are highly unsteady coupled systems that cannot be designed or modelled in isolation. Denton [15] discusses several gross limitations in more detail. Denton [16] studied loss in terms of entropy for a wide variety of mechanisms. As noted by Denton, a physical understanding may be more valuable than quantitative prediction based on empiricism. Later Denton  Pullan [17] study blade endwall secondary flow loss. Cascade measurements were shown to be of limited use. They proposed full stage calculations with well defined inlet conditions and investigations into the specific loss sources for a particular design, before design changes. Today much design is still based on 1D correlations and steady CFD being performed with little inflow data if any, in an attempt to verify improvements. A saturation point has been reached whereby uncertainties in modelling mask well founded design alterations. Current design systems have high uncertainties, so designs are often compromised by adding safety margins. These uncertainties range from problem definition (do we really know the true boundary conditions?) to modelling errors (discretisation and turbulence modelling). Current technology makes use of higher favourable and adverse pressure gradients and often high curvature blades, invalidating standard RANS calibrations. For example, shock boundary layer interaction modelling has two key deficiencies. (I), the shock foot location and strength and (II), the subsequent separation behaviour, all of which are unsteady. Clearly neither of these are in the comfort zone of typical steady RANS models, historically developed for attached boundary layer flows. The influence of coupling zones with RANS content also raises fundamental RANS modelling issues. For example, labyrinth seals generate low Mach number wake turbulence from seal teeth but a high rotational Reynolds number and hence a high Re boundary layer also develops [18]. Such Mach number disparities also raise numerical smoothing issues. In rim seals, the circumferential integral length scale is many times greater than that of the meridional plane. The implied lengthscales from RANS can also be wildly inaccurate. Jefferson-Loveday et al. [19] use a Hamilton-Jacobi differential equation to modify RANS length scales, improving drag coefficient prediction by 20-30%. An inflow-outflow condition is also developed to replicate pumping of fluid in and out of a rim seal [20]. Hence modest increases in geometrical complexity can have profound flow and modelling impacts. Inflow conditions in turbomachines are seldom well defined for modelling, even by gross quantities, let alone detailed fields. Even if they were known, as noted by Spalart [21], the k , k and S-A models are insensitive to freestream conditions, which is obviously incorrect from a physical standpoint. As an engine cycles through different operating conditions, the geometry changes, so not even the designed geometry is well defined. These may be seen as problems to modelling, or alternatively, an opportunity to realise large gains! User experience and the decisions made are often overlooked as one of the key controls over such uncertainties, which lead to the previous systems suggested to provide consistency for industrial use [22]. Sandberg and Michelassi [23] suggest extremely accurate CFD is required to realise order of efficiency 0.1-1.0% gains. From a flow resolution standpoint (ignoring problem definition fidelity and uncertainty), the ideal of Direct Numerical Simulation (DNS), is currently out of reach for most industrial flows due to high Reynolds numbers and impractical grid demands as DNS requires ORe 2.65 nodes [24]. Many detailed data sets have provided physical insights at low Re [25][26][27][28]. Most current CFD is still comprised of steady RANS modelling and acheiving improvements of the order of 0.1% seems out of reach. Modern turbomachinery design employs a wide range of modelling fidelities during different design stages and for different engine components. The required spectral gap to enable URANS instead of RANS often doesn't exist [29] and current trends -smaller, faster rotating components; higher frequency vortex shedding from thinner trailing edges; greater degree of transonic to supersonic flow with shock-BL interaction [23] point towards fewer flows able to utilize URANS. Indeed, turbulence is broadband, and the notional separation of resolved and modelled scales makes little sense. Currently in gas turbines 80% losses stem from BL flow [23]. The above suggests (U)RANS may completely fail for current applications which could result in redesigns costing millions to hundreds of millions in additional fuel and maintenance costs to manufacturers, airlines and adversely impacting the environment. Most eddy-resolving methods on the other hand, without largely resorting to DNS, will exceed current accuracy levels with ease where RANS fails, bringing about great benefits. Extreme resolution is not necessarily required.
High fidelity methods such as Large Eddy Simulation (LES) have been successfully applied to a wide range of complex flows where traditional (U)RANS modelling does not sufficiently represent flow physics to be consistently accurate [30]. Using LES and, for high Reynolds number flows, hybrid LES-RANS, critical knowledge has been extracted and exploited both to inform designs and to improve lower order modelling. Although LES is often cheaper than rig testing, it is still considerably more costly than RANS. Hence this ability to improve lower order models in a more general sense is key.
With adequate problem definition and mesh resolution, low Reynolds number and wake dominated flows can now be used to replace a substantial amount of rig testing. At higher Reynolds numbers, jet aeroacoustics is one area of application using hybrid LES-RANS [31]. When installed, ultra-high high bypass ratio jets interact strongly with the airframe [32,33], requiring these to be considered together. Noise from ultra-high bypass ratio jets installed near 3D wings under flight conditions has been accurately predicted [32]. Such coupled flows are difficult to accurately reproduce experimentally. Using LES, sound sources, length and time scales have been located and quantified in 3D [34]. Three-dimensional mean flow, turbulence and unsteady data recorded provide a wealth of potential knowledge. The most ideal flows can even be optimised using LES. For example using a genetic algorithm, turbine trailing edge cut-backs [35] (600 geometries), internal  [4], (c) NASA/ Boeing SUGAR concept using a Truss-Braced Wing (TBW) [5], (annotations added to [6][7][8]). cooling geometry [36] (10 geometries) and labyrinth seals [37] (80 geometries). The turbomachinery industry is now actively bringing LES into their design cycles, however cost will remain the long term constraint for eddy resolving methods. Areas of high impact are transitional, separated, and mixing flows and those with heat transfer. LES is expected to allow progress to be made for the more complex geometries found in industry. Best practices are hence critical and these are being established using robust and mature technology such as second order finite volume (FV) methods [31].
Transition modelling can also impact, for example, LPT performance or HPT heat transfer prediction. For the LPT Vadlamani [38] matches the unsteady transitional flow physics of an LPT with passing wakes. Pichler et al. [39] study LPT rotor stator axial gap effects performing a loss breakdown and identifying unexpected kinetic energy growth away from the blade. For the HPT, high heat flux unsteadiness can be observed even before transition due to Klabanoff streaks. If hybrid LES-RANS is to be used around such leading edges, either significant advances in RANS modelling to account for flow acceleration, curvature and pre-transitional behaviour, or novel LES/RANS zonalisation procedures require development.
Heat transfer is becoming more important in the adequate modelling of compressors and turbine temperatures have continually risen [40]. To prevent blades deteriorating, internal and external cooling flows are introduced. There has been a long standing demand for accurate heat transfer modelling under these chaotic flow conditions. Such demands can be largely met by LES, yet industrial use remains low, instead relying on costly experiments that may deliver data late in the design, where changes are expensive. For example, in 2005, Guo et al. [41] study film cooling holes injection angle, blowing ratio and resulting flow anisotropy, naturally resolved using LES. Sewall and Tafti [42] and Tyacke and Tucker [43] study internal cooling duct ribbed passages using LES. Such LES can be run for less than $500 compared to a rig test costing approximately $50,000. Trailing edge cutbacks are optimised using a perfectly parallel genetic algorithm by Watson and Tucker [35]. All of these use meshes of the order of 1-10 million cells, which in the shadow of Exascale computing (demonstrated 2018 [44], Argonne Labs, Aurora due 2021 [45]) should now be commonplace.
Tyacke and Tucker use flow categorisation to suggest wider application of eddy-resolving methods [22].
Although flow over isolated blades or mid span sections is performed routinely in academia and is making inroads into industry, multi-blade and multi-row effects are critical in revealing the true performance of real systems. Passing blade wakes and tip leakage vortices can heavily influence the main flow as shown by McNulty et al. [46]. Forward sweep was found useful to improve stall margin and reduce loss, highlighting the need for reliable 3D design and validation tools. Typical RANS performs poorly outside of a working range and Uncertainty Quantification could be used in future to provide confidence and quantitative evidence of risk in design. Saito et al. [47] study a two stage compressor with stator LE and TE leakage flow using DES with a 450 million cell mesh (104 blades). Yamada et al. [48] study stall in a seven stage axial compressor using 2 billion cells (650 million in the first two stages). Here stall cells in the first two stages grow to 1/ 3 rd annulus. In the 5 6 th stages, leading edge separation leads to rotating stall and large scale blockage in up and downstream rows. Clearly full annulus eddy-resolving calculations meet significant cost constraints and accurately predicting the sensitive flow physics in a compressor is essential to predict design and off-design performance.
We are now at a point where turbomachines and their interaction with future airframes need to radically change to meet new targets, putting designers in unknown territory. With ever shrinking design cycles and rapidly dropping computational cost, high fidelity simulation is now being absorbed into industry to supplant experiments, inform design tools and reveal unsteady flow physics in astonishing detail. In this paper we review some promising methods, note key challenges and identify paths forward.
The paper is structured as follows. Firstly, aspects of zonalisation and coupling are discussed along with some significant challenges facing the use of high fidelity methods in turbomachinery. The potential of high order schemes in terms of accuracy, cost, SGS modelling and its current uses is considered. Validation of high-fidelity methods and their use in lower order model development are examined. Barriers and suggestions on how to move forward are then discussed before summarising.  Each requires some variation in modelling methodology. Zonalisation and coupling of these zones is hence a critical issue in terms of both turbulence and geometry modelling fidelity [22]. Both these aspects are used to deal with multi-scale flow. Geometry representation is here considered more fluid in terms of fidelity, using immersed boundaries and body forces as well as traditional body fitted meshes. Fig. 3 shows such a hierarchy. Furthermore, design now requires multiphysics modelling. As noted by Dawes [49], geometry links different disciplines, perhaps requiring geometry definition in terms of level sets, allowing simple geometry modification during design. Combined physical modelling may include conjugate heat transfer for film, cut back trailing edge and internal cooling as [50][51][52][53], finite element analysis and aeroelasticity, acoustics [31,34,[54][55][56][57][58][59][60][61][62][63], particle transport (volcanic ash [64,65], sand), water droplets (rain, steam turbines), ice accretion and shedding [66] and reacting flows [67] among others.
Where it does not make sense to integrate software aimed at different physics, coupling software can be advantageous. CHIMPS, used for the Stanford whole engine simulations [68] where RANS is used for the compressor and turbine and coupled with and LES for the combustor is one example. The MpCCI library [69] also offers multiple grid functionality including interpolation with user defined methods. Other available software includes CORBA [70], designed for heterogeneous computing and interdisciplinary coupling and OpenPALM developed by Cerfacs and ONERA. Using OpenPALM, conjugate heat transfer for a combustor is performed by Duchaine et al. [71] and for a cooled turbine blade [72], combustor LES coupled with radiation modelling is also performed by Poitou et al. [73]. Aspects of zonalisation and coupling are now presented. Note, here, the term coupling does not refer to linking of different codes.

Hybrid RANS-LES turbulence modelling
LES can be selectively used for low Reynolds number flows such as those found in the LPT. The use of LES for high Reynolds number flows is generally still restricted to academia. For these, hybrid LES-RANS methods appear to be a clear path forward. The most well known and commonly used is Detached Eddy Simulation (DES) with its variants such as Delayed DES (DDES) [74], Improved DDES (IDDES) [75]. This has seen success for massively separated (wake type) flows. For more general cases Deck [76] suggests a zonalised approach to widen its use for flows with no clear separation point and later improvements [77]. However, as DES relies on RANS for the entire boundary layer, it inherits many of the deficiencies in the underlying RANS model which is typically a linear eddy viscosity model and can rely heavily on the grid filter. Aside from this, turbine blade film cooling is achieved with complex arrays of jets emanating from the wall. These inject locally large scale unsteadiness into an assumed (under DES) quasi-steady boundary layer flow. Recently, Sreekanth and Woodruff [78] introduced a RANS-LES zonalisation, whereby discontinuities in turbulence quantities can be avoided by substitutions in the governing equations. This avoids the well known log-layer mismatch but requires further testing.
Other well established zonal approaches such as Tucker and Davidson [79] require more input via the user or a control system, but allow explicit turbulence modelling and filter definition. Here, use of RANS in the inner layer only, brings a significant reduction in computational cost, a more diluted RANS contribution (< 1% of the boundary layer) and allows more fine tuned control over RANS and LES content if required. This may be required in transitional flow or zones prone to separation, however, this approach has been quite robust in its basic form [80,81]. The recent implementation of this zonalisation using a blending function [57] naturally enables the definition of RANS, LES and hybrid LES-RANS zones and synthetic turbulence generation between RANS and LES zones. For sensitive flow regions such as low surface curvature separation and transition in compressors and turbines, on-the-fly modifications could be applied using zonal information and flow metrics. In the HPT where the Reynolds number is high and hybrid LES-RANS is required, embedded LES zones may also be useful. For example near the leading edge where there is rapid eddy distortion, high flow curvature and leading edge transition, or downstream of film cooling holes where significant large scale mixing takes place near the heated surface. Sagaut et al. [82] use LES only in the transition, trailing edge and wake region of a T106 LPT blade. RANS is used elsewhere, with Non-Linear Disturbance Equations (NLDE) [83] used to extract fluctuations from RANS to provide to the LES zone. Nonreflecting boundary conditions are also applied at the LES interface to avoid spurious acoustic fields.
Significant performance penalties have been observed for compressor blade roughness [84], especially near the leading edge, the thickened boundary layer effecting blockage due to shock interaction. The fine scale topology of roughness elements can have significant effects on transition, drag and heat transfer [85,86]. Such topology and degree of roughness will vary over each component and over the life of the engine due to manufacturing variation, abrasion and accretion of debris from pollutants and dust, determined by flow paths and temperatures. Any near wall modelling should attempt to reflect this. As noted, local boundary layers are coupled with the passage flow, the effect being at the sub-component level.

Turbulence interfaces
Another aspect is the availability of data indicating inflow conditions. For example, even for carefully carried out academic experimental turbomachinery studies, even mean flow data is often lacking. The supply of critical integral scales of turbulence is severely lacking. This forces us to attempt to insert suitable mean flow and turbulent content at boundaries of domains and interfaces between modelling regions.
Sufficiently accurate inlet turbulence can be critical in some zones. Thomas et al. [92] study surface temperatures of high pressure turbine vanes downstream of a combustor. Unsteady combustor turbulence was required, constant 2D inlet conditions producing persistent secondary structures generating incorrect surface temperatures. Inflow turbulence should be of sufficient form and detail to generate the correct downstream flow. An example for an LPT in Fig. 4, shows types of turbulence content needed. This includes freestream turbulence based on a turbulence spectrum, wake and boundary layer turbulence, with the correct level and distribution of correlated fluctuations. These can be generated separately and combined [38], however this is not ideal for everyday use. A key inflow requirement, particularly for turbomachinery, where components are axially close, is that real turbulence develops rapidly and in a short distance. The inclusion of synthetic turbulence must ensure realistic conditions at, for example, the leading edge of a blade rather than the inlet face itself, or an arbitrary measurement location. Ideally methods also adhere to physical requirements such as avoiding continuity errors (being divergence free for incompressible flow) or acoustic contamination. The capability to reproduce anisotropic Reynolds stresses is also desirable, particularly for  Table 1 indicates advantages and disadvantages of common approaches. Previous reviews of inflow generation methods are provided by Sagaut et al. [82], Keating et al. [93] and Tabor and Baba-Ahmadi [94]. The most pragmatic and general approach appears to be the use of synthetic turbulence rather than use of precursor simulations, random perturbations or boundary layer flow where Lund type recycling (Lund et al. [91]) is often used. For recycling methods, spurious periodic signals from recycling ( f U L / r c r ) can arise. Other limitations relate to appropriate scaling for different flows, for example for compressible flow or those with strong pressure gradients. In practical terms, there are two broad considerations -(I) the turbulence generation method including its capabilities and limitations, and (II), how this is applied to the flow. The latter point also eludes to the application of these methods to RANS/LES modelling interfaces. Schlüter et al. [95] interface RANS-LES-RANS for a compressor-combustor-turbine simulation. On the fly time averaging and scaling of precursor LES fluctuations using RANS is used to interface upstream and downstream regions. Shur et al. [96] focus on matching boundary layer stresses by overlapping RANS to LES regions, enabling a fast switch to resolved LES from upstream RANS. Forcing terms have also been used to reduce the development length for resolved turbulence, however, this must be used with care, especially if considering acoustics, where forcing could contaminate the acoustic field. In this regard, the STG approach of Shur et al. [97] with a damping layer shows promise. A potentially infinite time series based on correlated Fourier modes is generated a priori. This avoids frequencies appearing in the  flow related to recycling of inflow turbulence. Forcing can be applied over a volume rather than at a plane or face. At the downstream end of the volume, statistics should match the defined statistics. This has been shown to greatly reduce acoustic contamination for compressible flow when passing from RANS to an LES region [97]. Sponge regions can be used to isolate traditional characteristic boundary conditions from flow unsteadiness. Where domains can be numerically extended, this could avoid use of non-reflective treatments which are often sensitive to tuning parameters. One aspect which appears not to have been tackled is synthesis of inflow turbulence generated by upstream passing wakes. These introduce significant mean and turbulence profile variation. It is possible predicted RANS lengthscales and mean velocities can be incorporated into VSTG methods to reflect this. For example, Shur et al. [97] provide lengthscale functions for channel, boundary layer and shear layer flows. An alternative, as noted in Section 2.4 is to model additional components at lower geometrical and/or turbulence modelling fidelity to generate suitable inflow or outflow conditions. For example, upstream bypass geometry for jet flow or downstream guide vanes for Fan-inlet coupling or Fan-OGV interaction. This also enables adequate coupling of different components and zones.

Blade rows
The flow through single or multiple blade rows contains both deterministic and chaotic flow leading to tonal and broadband noise. Rotor tip leakage and wake flow can impinge on downstream stators leading to upstream and downstream propagating acoustic waves. These, along with entropy and vorticity waves, can create a feed back loop within the flow. Prime numbers are often used for rotors and stators to avoid resonance and fatigue, requiring, ideally, full annulus calculations for real geometries. Changes to blade counts to enable single passage domain periodicity to be imposed changes throat areas, loading and associated flow frequencies. Aside from this, long wavelength disturbances such as compressor stall cells may span several passages making full annulus calculations the only solution. Using URANS and GPUs, Pullan et al. [98] revealed the route to compressor spike stall, this being caused by flow scales much larger than modelled turbulence.
For blade row calculations, cost reduction using ideally single passages remains challenging for multi-blade and multi-row calculations. These usually require sliding planes, where suitable interpolation is necessary to couple different zones whilst propagating numerous scattered waves. Reducing the cost by modelling a single passage is possible but not without issues. Standard periodic boundaries produce non physical synchronous vortex shedding. This makes some sense for URANS, but for eddy-resolving approaches, there is no physical synchronicity. Phase lag boundary conditions can corrupt turbulent spectra as the entire broadband flow cannot be adequately represented and stored [99]. Blade row, secondary flow, wake, boundary layer and shock interactions all contribute to the complex flow field. A detailed investigation is undertaken by Mouret [100]. Alternatively the computational time plane can be inclined, however this also has physical and numerical limitations, potentially requiring more than a single passage [101] and requires awkward post processing. Again, time inclining indicates a lagged time shift which is still non-physical. Hollow, 'soft vanes' which contain Helmholtz resonators may be used in future to control rotor-stator noise [4], which again may invalidate attempts to reduce modelled section size. Clearly the cost of resolving full annulus calculations and multi-stage calculations is limiting, indicating component and system level studies are necessary at different fidelity levels. With time, each will naturally progress up the geometry and turbulence modelling hierarchy.

Use of body forcing
The use of immersed boundary methods (IBM) and similar approaches such as IBM-smeared geometry (IBMSG), assuming an infinite number of blades [102], enable greater annulus sectors to be modelled at reduced cost. This allows the full annulus to be modelled at lower fidelity and cost, or could be used to enclose a body fitted mesh region, placing the effect of periodic boundaries far from the region of interest (ROI). Fig. 5, shows examples of increasing fidelity, which could be used to reduce flow development cost or place numerical boundaries far from the ROI. Here multiple arrows represent the equivalent infinite series of blades (IBMSG), dashed lines approximate geometry (IBM) and full lines wall resolved geometry.
Alternatively, positioning of key components can be studied such as fan-inlet or rotor-stator axial location without re-meshing. Cao et al. use an IBM to study fan-intake interaction under high incidence [103]. Either a suppressed level of post-separation distortion or an increase in the separation-free operating range was found. Fan inlet distortion can also be modelled using a localised version of IBMSG. Fig. 6(a) shows a body fitted mesh contrasted with a background mesh with localised IBMSG (Fig. 6(b)) to generate blade passages at reduced cost. A 1/3rd annulus distortion generator is also added using a classical IBM ( = u 0). An instantaneous flow for the rotor with inlet distortion is shown in Fig. 6(c). Fig. 7 shows total pressure distribution at the rotor inlet for two fan speeds. IBMSG is in good agreement with experiments, body fitted mesh and IBM results.
An aeroengine mounted to an airframe poses a huge range of scales from wall streaks and blade wakes, to the largest jet plume scales of the order of the jet diameter. Such multi-scale simulation is enabled with multi-fidelity turbulence and geometry modelling. Internal turbulence is generated using body forcing in the bypass duct of an aeroengine mounted on a pylon-wing-flap geometry (Fig. 2 frame D) to study the effect on jet flow [104]. It was found to promote the weaker, inner shear layer transition to turbulence. The effect of which is much weaker than the dominant outer shear layer. This is particularly so with increasing bypass ratio nozzles, which become closer and closer to single stream jets. As shown in Fig. 8, using body forcing, plasma actuators have been modelled, indicating the voltage required to eliminate intake separation is outside current technology levels [105]. Body forcing clearly has wide application in scoping studies, active flow control, introducing complex geometry into high fidelity simulations, representing geometry of disparate scales and imposing engine realistic boundary conditions. Although body forcing itself may be considered a low fidelity method, it is a key enabler for high-fidelity simulation of complex systems.

Heat transfer
Turbine inlet temperatures can reach 1500-1750 K, exceeding the melting point of the blade material, and have risen by approximately 7 K per year since the 1950s [40]. Since the 1960s, air from the compressor has been used to cool blades internally via turbulence inducing ribs and pins, the air of which exits to envelop the blade externally via film cooling holes. This has enabled vastly extended periods of operation, at the expense of compressor air and efficiency. It is desirable to lower the temperature at the blade root, where the stress is high and the tip, which is easily damaged by abrasion. Further damage can result from creep, fracture, fatigue and corrosion (enhanced by high temperature). These numerous coupled elements are indicated in Fig. 9. For creep, a rise of 10 K can halve the life of a turbine blade [40]. Careful control of cooling and temperature distribution throughout the engine cycle is therefore critical.
Ceramic matrices incorporating supporting fibres are likely to be used to help protect turbine blades [4], presenting new composite material properties, life and surface roughness challenges. This will require the accurate modelling of thermal and structural interfaces between materials and advances in turbulence modelling of a variety of surface finishes, which may change throughout blade life. Here, smaller more canonical cases can inform models used for design verification at larger scale.
For many years, it has been possible for LES to be used in industry to study internal cooling of turbine blades due to large scale separation and mixing [30,43]. For a two-pass duct, less than 10 million cells were required. Note that inflow turbulence is not needed here due to rapid development of large scales. This is shown in the left of Fig. 10 in addition to complex 3D corner flows. LES can be used to obtain accurate  heat transfer coefficients, removing the large uncertainty in measurements or RANS modelling. This is displayed to the right of Fig. 10, where the range of Nusselt number prediction using different RANS or SGS models, is reduced by a factor of four. Mayo et al. [107] study high rotation number effects on flow and heat transfer using simplified ducts of only 1 million cells. Real multi-pass geometries are highly complex, adapted to fit within blade structural constraints. These have been tackled successfully [108] using meshes up to 100 million cells, LES responding physically to the rapidly changing duct geometry and blockages such as pedestal banks in the trailing edge cut back. This removes the need for re-calibration of RANS models, which will not be general and may lead to poor heat transfer prediction. The same arguments can be made for trailing edge cut backs, an optimisation of which is reported by Watson and Tucker [35] using less than 1 million cells per case. LES should at least be used to check a final design before rig testing, which may not even be required.
Harnieh et al. [109] study film cooling of a Nozzle Guide Vane (NGV) in a fully coupled sense with internal passages producing coolant jets with the correct swirl and mixing with the external flow. Also, in contrast, a constant steady mass flow boundary condition is imposed at the coolant holes (with no resolved internal flow), adversely affecting heat transfer prediction. This is attributed to lack of swirl at the coolant holes exits and mixing at the internal-external flow interface. As noted in the introduction, hot gas can enter cooling holes, hence a representative plenum or inlet/outlet boundary condition with the correct mean flow characteristics may be required to capture this more effectively without resolving detailed internal geometry.
Conjugate heat transfer analysis is becoming important not only for turbines (for example, internal and film cooling) but for compressors [99], where adiabatic boundary conditions are no longer enough to accurately predict performance. The disparity in solid and fluid time scales ( = c c h k / T p 2 ) is often O (1000 10000) posing significant challenges in computational cost. One potential solution is that of solving the fast varying fluid and the slow varying solid in the frequency and time domains or disparate spatial scales using the block-spectral method of He [110]. The selective use of different methods in focused regions can hence bring about large overall benefits in accuracy, cost and flow detail, often unobtainable experimentally. The key is to make the combination of these methods more integrated within the solver, requiring a commitment to development, such that minimal user intervention is needed.

Particle transport
In a study to measure volcanic ash deposition of real engine geometries, Kim et al. [64] find turbine inlet and vane surface temperature to be key parameters. This is also influenced by blockage of cooling holes by dust that may or may not be deposited. This seems to invalidate any steady modelling, as the blade may deteriorate resulting in geometric changes and particles may rebound. Also important is that ash was also deposited on the temperature probes. Instrumentation effects on turbine blades have recently been studied using LES by Ubald [111,112]. This could be extended to study debris and fouling effects on measurement uncertainty. Cheng et al. [113] study turbine cooling hole  dust deposition, noting the difficulty of particle resuspension and predeposited particles that may significantly affect results. Grant and Tabakoff [114] use a Monte Carlo method to model particle rebound dynamics material removal from stationary and rotating blades. Particle size, relative location and restitution ratio are chosen from statistical models. Key erosion of the compressor blades includes the cutting back of the leading edge, thinning of the trailing edge and tip damage. This poses geometrical representation challenges and questions URANS validity at the trailing edge. To deal with geometry changes, Dawes [115] has suggested the use of implicit geometry via level sets to allow rapid re-meshing using unstructured meshes. For structured meshing, Ali et al. [116] use level sets to determine mesh blocking topologies. Alternatives may also include the use of immersed boundaries in wear regions which could be eroded. This would also require additional corrections such as wall distance calculation for turbulence modelling.

Computational benefits
Higher-order ( > O 2) methods offer the potential to achieve the same or better accuracy than second order methods for a reduced computational cost, which is highly attractive. For example, using fewer cells and local, compute intense operations as in FR methods. This more efficiently utilises available floating point operations (FLOPs) on modern computing architectures. For comparative purposes, in d dimensions, FV methods have one DOF per cell, and nodal high-order methods have O p ( ) d [117]. A side effect of higher computational density is that parallel partitions can be as small as one cell. This also better enables typical MPI latency hiding where computation is carried out whilst communication is performed. Fig. 11(a) displays the vast difference between the estimated number of DOF, required throughout a civil aeroengine for a 2 nd and  4 th -order finite difference (FD) scheme. The use of the 4 th -order scheme reduces the LES cost to below that of a 2 nd -order hybrid LES-RANS solution. Use of the proposed IBM however are not currently compatible with higher order methods, the external flow generating numerically problematic oscillations within IBM zones. Here additional smoothing may be sufficient where low geometry fidelity is required. Current finite volume solvers running on CPUs may only attain 5 10% of their peak performance. This is the theoretical maximum work throughput, for numerical simulation, generally measured in floating point operations per second (FLOPS). On idealised geometries, FR methods can reach 58% peak performance [118], potentially offering a large cost saving. It should be noted however, although this clearly showcased the potential scalability and computational efficiency of FR methods, very little useful data was recorded. The Reynolds number in the Low Pressure Turbine (LPT) is low and accurate prediction of unsteady phenomena are already tractable for industry using FD or FV methods. In practical terms, the more beneficial route for the LPT would be to reduce the computational cost (using fewer degrees of freedom at high order) to obtain the same or better accuracy than other methods.

Accuracy and cost
The accuracy and cost of a simulation is linked to the numerical algorithms employed. Here two codes are contrasted. These are a compact finite difference scheme (code A), and an optimised flux reconstruction scheme (code B). Particulars of each code are detailed in Table 2.
To contrast these, the well known Taylor-Green Vortex from the high order workshop is chosen. This reflects flow transition and decay found in many turbomachinery zones and highlights sensitivities in the spatial schemes chosen.
As far as is practically possible, cases are setup identically using 64 3 degrees of freedom and run for 20 time units. The 4-stage Runge-Kutta scheme is used for each code with a small time step of 10 3 s to isolate differences due to the spatial scheme. Fig. 12 plots the variation of dk dt / and enstrophy with time ( Fig. 12(a) and (b) respectively). The benefit from using higher order schemes is clear. Dissipation is greatly reduced, which may become more important as greater regions of the engine are modelled, requiring resolution of eddies over long distances. On a single CPU core, Code A using an explicit 4 th -order FD takes 0.52 s/iteration and 0.92 s/iteration using the 4 th -order compact scheme. Code B takes 0.97 s/iteration on a CPU core and 0.019 s/iteration on a GPU.

SGS modelling
As with RANS, common SGS models are based on local flow quantities and are of similar form. The Smagorinsky eddy-viscosity model, which is similar to the RANS mixing length model, is provided in Equation (1). This forms the basis for the most widely used SGS models. It contains a 2 term, akin to a dissipative truncation error. Moin et al. [119,120] show numerical errors can dominate SGS modelling at low order. A non-dissipative numerical scheme admits the sensible use of an SGS model. At high order, SGS modelling can dominate the numerical scheme, destroying any accuracy advantage. For industrial codes, high numerical smoothing can occur, an example being the Roe smoother at low Mach number [121]. In this case, an additional SGS model is inappropriate, and Implicit LES (ILES) is preferable [30].
With regards to SGS modelling itself, the standard Smagorinsky model has several negative traits including, being purely dissipative (dependent on C s ) and a non-vanishing eddy viscosity at walls. This can delay or even prevent transition to turbulence [122]. It is soley based on local variables, does not generate an anisotropic eddy viscosity to represent anisotropic turbulence and along with many other turbulence models has questionable validity in the inertial subrange at low Reynolds numbers [123]. Typically a range of length scales of at least = L / 1000 i is needed before any such range can be detected. Large regions of a gas turbine do not meet this criterion [124]. Even so, the Smagorinsky model forms the basis for many widely used models due to its simplicity. Examples include the Germano dynamic model to modify C s based on resolution, using a test filter (usually 2 ), requiring averaging procedures on µ SGS for stability [125]. Hughes [126] presents the Variational Multiscale Method (VMS) in which an additional filter is used to define the interaction between the smallest resolved scales and sub-grid scales. Vreman modifies the Smagorinsky model to better suit transitional flows [127]. The Wall-adapting local eddy-viscosity (WALE) model [128] and σ-model [129], fix gross failures of the Smagorinsky model. These are contrasted for transitional and endwall flow of a compressor blade by Scillitoe [122,130]. Many other functional and structural SGS models exist, some allowing non-local and non-equilibrium effects to be introduced via transport equations [131]. However, these often do not display good agreement with true subgrid stresses [132]. Non-linear SGS models allow modelling of energy backscatter and anisotropy, but are not widely used due to instability without additional dissipative elements that form mixed models [133].
It can be shown mathematically that some high order numerical schemes contain both dissipative terms and non-linear terms akin to the Clark model of explicit nonlinear SGS models [134]. Typically they also dissipate energy at high wave numbers, moving higher still with order. This appears a desirable trait, whereby sufficient dissipation occurs near the grid limit, allowing maximum resolution of grid scales as demonstrated by Wiart and Hillewaert [135]. This is particularly important for wave propagation in areas such as aeroacoustics where dispersion and aliasing can contaminate fields. An explicit SGS model, with questionable validity at low Reynolds numbers, is thus likely to create more problems than it solves. Most high order solvers appear tuned not to require an explicit SGS model, becoming Implicit LES (ILES) or Montonicity preserving ILES (MILES) [191].

Current uses
An insight into academic high order method penetration is gleaned from the International High Order Workshop [136]. Here, flows are classed into Verification and Validation, Advanced, and Computational Challenges as shown in Table 3. Flow physics challenges include transition, laminar separation, transition and reattachment, transonic flow and shock capturing, and inlet turbulence. These can all be sensitive to numerical properties, so careful scheme selection is required.
Based on FR, Dawes et al. [137,138] provides examples of more industrially relevant flows with relative computational costs to achieve similar accuracy to second order methods. The time step is local and adaptive through the duration of the simulation requiring close attention to parallel load-balance. Table 4 summarises the cases. The speed up of local time stepping over uniform time stepping is denoted ψ. The aeroacoustic models obtain a large speed up but the Octree mesh distribution is clearly non-ideal for the jet, but convenient for the time stepping scheme used. The jet case benefits from an additional 50% Table 2 Code details. (Note, FD = finite difference, E = explicit, C = compact, FR = flux reconstruction, DG = discontinuous Galerkin, SD = spectral difference, RK = Runge-Kutta). increase in far-field sound cut-off frequency moving from 2 nd to 3 rd order. Similarly C p distributions on the landing gear indicate a benefit at 3 rd order. By use of mesh, hardware and points-per-wavelength scalings, Dawes estimates the HOTnewt solver is approximately 2.9 × faster to an equivalent second order solver such as ONERA/ CEDRE. Given the efficient hardware, the cost is estimated to be 29 × cheaper. This is an example where building from the ground up using new technologies, progress can be made. It is far more difficult however to take mature codes and re-invent them. Marty et al. [139] use the elsA solver of ONERA with AUSM+(P) with a high-order MUSCL extension to reach 3 rd and 5 th order. This is applied to study the laminar separation bubble of the T106C LPT in conjunction with the WALE SGS model. Zonal-DES is also applied to a transonic compressor with a shock near the blade tip. Zonal-DES is also applied to a fan rotor in order to assess fan wake-OGV interaction and noise. Higher order MUSCL interpolation resolved more small scales, however, whatever the order of MUSCL interpolation, ZDES compared well with averaged and unsteady statistics, indicating larger scales were most important. This aligns with the fact that large scales are more efficient at generating noise [140].
Vadlamani et al. [141] use the COMP-SQUARE solver to study surface roughness effects on transitional boundary layers using a 6 th order compact FD scheme. Numerical instabilities are treated using a 10 th order non-dispersive filter. Using similar numerics [142], study non-equilibrium turbulence of a backward facing ramp. Non-equilibrium flow is found to be dependent on scales relating to boundary layer and wake turbulence. Suggested mechanisms are phase de-coherence and inter-scale spreading of disturbances.

Validation
The flow physics targeted by high-fidelity simulation requires more stringent validation than mean or bulk quantities. Tyacke and Tucker [22] indicate turbomachinery LES lags airframe validation levels. Turbomachinery mainly reaches mean flow with some Reynolds stresses and airframe validation predominantly reaches second order moments and spectral analysis. Integrated quantities were are also lacking for turbomachinery, likely due to access. High-fidelity CFD provides three-dimensional fields of high order statistics unlimited by access. Rather than replicating mean profiles from experiments, highfidelity simulation data now far exceeds that provided by experiments [32]. Measurements simply cannot match the quantity of detailed data obtained due to access limitations and disturbance of the flow. This is problematic as detailed multi-variable correlations are difficult if not impossible to obtain. There is also variation between experimental facilities. Mc Croskey [147] compares 40 NACA 0012 data sets from different facilities, concluding there is unacceptable variation and that "no single existing experiment is adequate". Difference between noise test facilities can vary by several decibels [31]. It seems useful to compare with multiple data sets if possible. Measured mean flow and Reynolds stresses can have 5% and 8% errors respectively [148]; 10% [149] for skin friction; 1.76% [153], 3% [150] for Reynolds stresses; 7% for pressure [151]. Heat transfer measurements are notoriously difficult to define and control, errors often in the range of 2-10% [152]; 12% [149]; 7-12% [150]; 8% [151]; 3.4-6.66% [153]. Even so, cases are often not sufficiently well defined for LES to be rigorously validated.  Table 2. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) Table 3 High order workshop test cases, with basic details in brackets. For example, inflow turbulence is rarely recorded. If provided it is incomplete, stating only turbulence intensity, or turbulence grill details. Experimental facilities can have various limitations, including unknown upstream flow within pipes; wake generating bars instead of blades; wind tunnel contractions and walls; additional flight stream shear layers in aeroacoustics; heat conductance into substrates (O10% [148,154]); and small scale geometrical irregularities from manufacture, tape and surface roughness. High-fidelity practitioners usually aim to model the real flow. However, it is impossible for CFD to identically match these due to computational resource, nor should the CFD be forced to inherit such deficiencies. Researchers more widely use high-fidelity simulation but industry holds significant but inaccessible measurement databases. Publicly available detailed inflow traverses are hence required to properly capture inflow turbulence dynamics. These could then be fed into synthetic turbulence generators or reduced order models to more precisely match experiment conditions. Duchaine et al. perform uncertainty quantification on a cooled turbine blade. Varying the coolant temperature in all holes, the LES data range with 95% confidence intervals overlaps measured data. It is hence important for both experiments and CFD to incorporate errors, confidence intervals or error probability distributions routinely. Rather than picking one superior dataset, we can be confident in validation where data overlaps.

Lower order model development
The limited design space covered due to computational cost of highfidelity methods at high Reynolds numbers limits what is achievable in design. High fidelity methods can be used to supplement experiments to gain deeper flow physics insight. Lower fidelity geometry modelling has been discussed in previous sections. Here, we consider lower fidelity turbulence modelling (RANS), which forms the basis for the majority of current CFD. Typically RANS provides some trend data, even if absolute accuracy is poor. A key application of detailed high-fidelity data is then to improve lower fidelity methods such as RANS and rapid design tools that may be lower than 2D.
Many RANS defects are known, but there are not many general fixes. Spalart discusses numerous turbulence modelling issues and fallacies [155] which should narrow potential fruitful avenues. Tucker [29] provides a broad review of current turbomachinery turbulence modelling defects and treatments. Issues include forwards and backwards transition; compressive/extensive strain, which can lead to suppression of leading edge separation and over-predicted heat transfer [156] and downstream accuracy degradation; curvature and rotational effects, eddy-viscosity models requiring sensitisation to curvature and additional body forces for rotation, both using a variety of terms and implementations indicating unreliable performance; freestream turbulence decay, which can incorrectly grow downstream in the absence of suitable destruction terms [132]; roughness, a simple model by Aupoix and Spalart [157] offsetting the wall distance, however surface roughness topology is known to be of significant importance, skin friction being 25% higher than sand grain roughness. This would appear difficult to effectively model in a meaningful way. Lack of the assumed alignment between the Reynolds stress and mean strain-rate tensors is found for compressor secondary flow [158] and labyrinth seal flows [18]. Monier [158] tests the quadratic constitutive relation (QCR) to improve modelled anisotropy, mainly seeing benefits upstream of the blade. Applied to a tip gap leakage flow, Monier [159] finds little improvement using the QCR. Case dependency is typical of RANS and its variety of proposed fixes. The use of URANS requires the existence of a spectral gap between unsteady large scales and modelled stochastic turbulence. Tucker provides many examples where such a gap is tenuous in turbomachinery [80]. Rapid eddy distortion can also occur around turbine blade and fan leading edges with inlet distortion or boundary layer ingestion.
Non-local non-linear effects are clearly important from a physical standpoint, however most RANS closures are based on local quantities, often with some convective element. In combination with the numerous issues noted above, manual analysis seems out of reach.
Further advances in RANS modelling may be achieved in the near future through the use of machine learning. The temptation is to use ML as a black box to try and force an LES solution out of RANS equations. We should instead be gaining and exploiting an understanding of the spatio-temporal changes in turbulence for different flows. Wu et al. [160] find the difference between the RANS and LES Reynolds stresses. This is then used to correct the RANS stress field. Forward propagation is then used (i.e. re-solve RANS with turbulent PDEs frozen) to get an improved velocity/pressure field. Here, training flows need to be close to test flows. How this is defined is an open research area. Wu et al. [161] later shows that any errors in the stress field are amplified at higher Reynolds numbers, hence modification to the eddy viscosity appears more general.
Weatheritt and Sandberg [162] use Gene Expression Programming (GEP) to build algebraic equations based on eddy resolving data whose terms can be studied. Weatheritt et al. [163] use similar techniques to suggest a non-linear expression for the Reynolds stress to improve prediction of the wake of an HPT blade where the Boussinesq approximation is found to be invalid. Great care is needed however, to prevent non-realizable or unstable solutions. However, regions of flow can be chosen and others ignored to train such models, in a sense calibrating a fix to a given range of flows. In many situations the complexity of the turbulence may render RANS too simple a palliate, even with AI for any substantial design parameter range.
A pragmatic approach is to identify trust regions based on flow physics to identify regions where RANS is known to be untrustworthy, for example breakdown of the Boussinesq hypothesis as Ling [164]. Later Ling [165] explores more general breakdown of assumptions using support vector machines, Adaboost decision trees and random forests. This gave a point by point breakdown of specific RANS assumptions with high or low uncertainty. Classifiers were then able to generalise to flows on which they were not trained. This may be more meaningful to a design engineer, trying to minimise computational cost. Unreliable RANS zones would then be more suited to LES modelling or lead to naturally flow adaptive model corrections. Further to this Ling [166] embeds Galilean invariance into a model to correct anisotropy trained on LES of channel flow, a duct, square cylinder, perpendicular and inclined jets in cross-flow and a converging-diverging channel. This improved predictions over baseline RANS on new cases including a wavy wall with separation and duct secondary flows. ML is still in its Table 4 Industrially oriented test cases ( [137,138] infancy and turbulence is highly non-linear and complex. Rather than blindly relying on ML, it seems prudent to take smaller steps towards improvements guided by knowledge.

Pre-processing
Parallel solvers and the use of High Performance Computing (HPC) is now common place. In 2009, Dawes [49] notes key changes required to CFD workflows, most notably the use of massively parallel processing of all aspects, including pre-processing, solution and post-processing. Still today, the solver is often the main parallel part of the process chain. As a rule of thumb, the largest cases usually translate to a run time of one to two months depending on HPC occupancy rates, larger cases use larger facilities and so the rule holds. Most human time is now spent in pre-processing and post-processing. Zonalisation and coupling at geometry level and turbulence modelling shifts significant time and effort into pre-processing. For example, one needs to perform mesh generation, and potentially IBM/IBMSG definition, synthetic turbulence inflow setup, RANS-LES zoning, numerical smoothing zoning and numerical flow instrumentation. Currently this requires significant manual effort which may outweigh the runtime of the simulation unless automation is increased significantly. Meshing for example is generally a human-led operation. The more complex the simulation, the longer it takes to generate, with specific attention paid to important flow regions for LES using a range of resolution criteria. For example, for boundary layers = + x 100, = + z 20, = N 2000 grid points per 3 , more generally, resolution of 90% TKE. For a wide range of flows this is not simple to automate, as it is not based solely on geometry and requires flow information. Fortunately, turbomachinery offers fairly distinct geometrical components and flows, where topological similarities can be leveraged and rules defined in advance. Many pre-processing operations also require nodal searches. For large meshes, which continue to grow, these can become prohibitively expensive at run time, where mesh partitioning and parallel reductions complicate searches. These operations should clearly be separated from the solver and can be faster sequentially than in parallel where MPI communication can saturate. K-dimensional-trees (KDTREEs) can help in many circumstances. It may be beneficial to completely separate the solver from input and initialisation via a greater degree of pre-processing as solvers are set for significant and rapid change.

Eddy resolving meshes
Chapman [167] estimated that when computing power reaches 10 14 FLOPS, LES would begin to compliment and perhaps replace some rigtesting. This was exceeded by a factor of ten in 2008 with IBM's Roadrunner. It should be noted that the use of overset meshes in regions of interest to avoid high cell counts away from boundaries is required to meet the grid point requirements of Chapman [167,168]. Most current meshes do not achieve this, for example structured meshes generally fan out from a boundary and do not coarsen in 3D, becoming larger than necessary. Choi and Moin [24] also recently reached similar requirements for wall resolved LES. 1 Choi and Moin estimate the number of grid points N required for a boundary layer scale approximately as Re L 2.65 , for DNS, Re L 2 , for wall resolved LES and Re L , for hybrid LES-RANS (wall modelled LES). Clearly, boundary layer adapted meshes will have the greatest impact at the highest Reynolds numbers where wall modelled LES is needed such as the engine fan or airframe wing, or at lower Reynolds numbers where wall resolved LES is necessary to capture transition such as the LPT.
The roadmap in the NASA 2030 report ([12] Fig. 1) study indicates, both wall modelled and wall resolved LES as low TRL for complex 3D flows at realistic Re when written in 2014. Large scale parallel mesh generation is indicated as high TRL, however the authors would argue this is not the case for LES of coupled turbomachinery-airframe flows. In this case, specifically flow adapted meshes are required to capture but not over-resolve key flow regions. In any case, an efficient computational mesh and modelling approach go hand in hand. Complexity is only set to increase, requiring more human input on each modelled element. To aid structured mesh generation Ali et al. [116], leverage level sets to gain useable block topologies for complex seals, an installed engine and an aircraft. Still, the mesh generation process needs much more parallelism and automation, possibly including adaptive methods that maintain mesh quality. Quality itself is poorly defined in the LES context and is related to underlying numerical schemes and modelling.
For current FV, large meshes mostly stem from connectivity (structured and unstructured) related to the boundary layer mesh. For hybrid structured-unstructured meshes, the surface mesh streamwise and cross-stream spacings are linked to final wall normal spacing and chosen aspect ratio at the structured-unstructured interface. For developing boundary layers, the estimates of Chapman [168] assumes hanging nodes, without which mesh density propagates outwards. This suggests use of non-conformal/hanging node/polyhedral meshes to match cell size with developing boundary layers and localised wake and transitional flow. For LES, the cell aspect ratio for the inner layer is normally assumed using wall units (for example, = + x 100, = + z 20 giving an aspect ratio of 5). Streamwise and cross-stream spacings could then be linked to the wall normal distribution, connectivity being another issue. Alternatively, the inner and outer layers can be considered separately. Fig. 13(a) shows a schematic indicating the disparity in mesh required through a boundary layer. The transition between these is non-trivial and example mesh topologies are indicated, though the transition would be more gradual. Two mesh requirements, particularly for boundary layers, are anisotropic elements and a gradual three-dimensional growth to isotropic cells. Where triangular or tetrahedral elements are required for connectivity, cell skewness limits aspect ratio and growth rate. Hence hanging nodes are attractive if they adhere to the above requirements to avoid sudden changes in cell size as in Octree meshes. Such a mesh is generated by Addad et al. [169] for a flat plate. For developing boundary layers, the boundary layer thickness and hence mesh requires a local scaling as indicated in Fig. 13(b) where a number of notional boundary layer mesh units are distributed based on the boundary layer thickness. For more general cases, mesh adaptation may be required at run time. Ideally the level of resolution in each direction should be the same, implying that for anisotropic flow, there is an optimal filter width in each direction. Toosi and Larsson [170] estimate sources of error based on two different grids. The optimum grid is assumed to have an optimal mesh and cell anisotropy distribution while achieving a uniform error. This provides a systematic procedure without resorting to adjoint methods, which are problematic for unsteady flow.
An example boundary layer mesh schematic and mesh filter scale to resolve 90% TKE with typical zonalised boundary layer mesh distribution for a compressor endwall are shown in Fig. 13. For boundary layers with correctly varying cell size, savings on the order of 10-100 times may be achievable [171]. In addition, boundary layer mesh may not be required in separated zones, offering further savings.
One benefit of FR methods is a greater resilience to mesh distortion. For an inviscid convecting vortex passing through a distorted mesh shown in Fig. 14(a), using a 4 th order FR scheme the mean error (point averaged l 2 norm) shown in Fig. 14(b) maintains order of accuracy, rather than degrading as with finite volume. Although promising, the use of higher order methods has major challenges before they become widely adopted. Not least of these is the requirement to use higher order meshes. At boundaries, a p th order scheme requires the boundary 1 We note Chapman's outer layer estimates are optimistic.
mesh to be C p continuous [172]. This may also require further improvements to geometry representation such as CAD. Unfortunately, even for established FV methods, complex geometry LES mesh generation can take weeks or months. This is then exacerbated with high order requirements. This bottleneck will not be easy to eliminate. Subsequent processing of lower order meshes may offer a non-ideal answer. One could mesh a compatible FV mesh and morph this to high order geometry, requiring only one extra pre-processing step. Such an approach which does not require CAD is presented by Ims et al. [173].

Time integration
Stability is a critical requirement for industrial use. A key aspect after spatial decomposition is the sequential bottle neck of time integration for thousands or millions of time steps. This limitation is also noted by Löhner [174]. Higher order methods generally have even more stringent CFL limits than second order. For explicit time stepping, the time step scales roughly as p 1/ after allowing for increases in cell  size [175]. Spiegel et al. [117] note p 2 CFL scaling for explicit time integration and p 5 becomes too expensive. These exacerbate the number of time steps required for large simulations or those with long or disparate time scales. For example, the blade passing frequency in a compressor is 1000 times greater than a surge [99], or for conjugate heat transfer, the thermal time constant of the solid is generally much greater than the fluid. For implicit time stepping, memory utilization may be limiting at high order as the bandwidth of matrices grows, scaling with p 6 [176]. Using implicit time stepping near walls, explicit elsewhere (as Shoeybi et al. [177]) and dynamic load balancing (at relatively few time instants), may in future help alleviate higher order memory requirements and achieve good load balance. Vermeire and Nadarajah [178] obtain speedups of over × 50 and a memory reduction of over 50%. Space-time cells or local time stepping may also be useful to allow more efficient time integration on a per cell basis. Lu et al. [179] achieve a speed up of up to × 100 over uniform time stepping. Other avenues are to perform parallel time integration such as the Parallel method [180] which involves a fine time step parallel correction to a sequential coarse time step solution. Alternatively the multigrid reduction in time (MGRIT) algorithm [181] involves a relatively non-intrusive wrapper to obtain a multigrid method in time.

Post-processing and data storage
A common issue for high fidelity simulation is that computational cost limits the simulation period. This is obvious for jet aeroacoustics, where low frequency resolution is limited by the number of periods data is collected. Similarly for high order correlations, low sampling times can limit confidence. Criteria for statistical convergence under a range of flow conditions are required. For example, initial transients need to be separated from periodic or pulsatile behaviour. Morton et al. [182] use deep dynamical modelling based on Koopman theory to extend time series with low relative error, based on short time samples, to assess potential flow control. This may allow the relatively short simulations to be extended somewhat, drastically reducing cost, or at least made more useful. For example, noise reduction concepts could be assessed based on dynamic models. However, further testing for more general 3D cases and caution is required.
Clearly we can and do store vast amounts of data when performing simulations. 2 Using neural networks, data encoding/decoding may reduce storage requirements by orders of magnitude [183]. The useable data is unfortunately largely wasted, as uncompressed, large datasets are computationally expensive and cumbersome to analyse. In the near term, with billions of DOF, visualisation and data co-processing may be required in real time owing to storage limitations. Furthermore, tensorial quantities such as 4 th order space-time correlations used in aeroacoustics can multiply dimensionality of data compounding storage issues and impeding comprehension. Tens of Terabytes are easily exceeded. Targeting Exascale, Ross et al. [184] from Argonne National Laboratory indicate the growing I/O gap between memory and compute performance and solutions ranging from hardware to low and high level software. Aids to HPC data analysis include Adaptable IO System (ADIOS) [185], a framework for I/O and data analysis, and Do It Yourself (DIY) [186], a library of functions for developing parallel analysis functionality. Some degree of on-the-fly modal representation of energetic scales may be useful as suggested by Dawes [137] to be reconstructed later. We tentatively suggest a notional Use Factor, UF = knowledge/cost. This is independent of the magnitude of data stored, emphasising the importance of effective analysis techniques. Due to sheer volume, any data that is stored may ultimately only be exploited using machine learning for such large data sets, to provide intelligible engineering insights and outputs.
For long term or large storage requirements, a glimmer of hope has appeared perhaps in using DNA as a storage medium. Depending on encoding technology, this enables × 2.2 10 3 PB/Kg of DNA [187] with a raw limit of 10 6 TB mm / 3 [188], transforming the ultimate storage limits and stable for centuries. The process to read and write has recently been automated by Microsoft and the University of Washington, taking around a day to write the word 'hello' [189], set to reduce significantly but likely requiring tiered storage and databases depending on data size and usage frequency. As data ages, access frequency usually drops as newer, better, more detailed data becomes available. Integrated databases were indicated as low TRL by the NASA 2030 report. Data from multi-fidelity simulation and experiments also needs merging. The source or age of data is relatively unimportant -it is the accuracy and uncertainty characteristics that should inform the engineer.

Hardware and algorithmic coupling
The rapid development of computing hardware opens up the possibility of using these tools not only on high performance computing (HPC) platforms, but on local desktops and distributed networks using combinations of CPUs, GPUs, Application-Specific Integrated Circuits (ASICs) and Field-Programmable Gate Arrays (FPGAs). The hierarchical nature of most computing systems involves compute units, multi-level memory systems and inter-connects between these at the chip and system level as shown in Fig. 11(b). This presents a highly case and algorithm dependent load balancing and code abstraction problem. There is then a link between the equations to be solved, the algorithms employed and the hardware that can efficiently be used.
Data structures are critical to code efficiency. A benefit of structured data (for example, i j k , , indexing) is that vectorisation is easily leveraged, multiplying the rate at which computation takes place on a CPU. However, for unstructured meshes with various elements, matrix operations are more efficient. To simplify code structure, NASA have used Fortran derived datatypes for different elements so that higher level operations are identical [117]. Matrix sparsity is evaluated and they are stored in a compressed format, Fortran 2003 type-bound procedure pointers then utilise compatible sparse matrix algorithms. Unfortunately it is not easy for compilers to automatically optimise such routines with derived types and at least code profiling, and most likely manual optimisations must be performed.
Abstraction for current hardware has been successfully demonstrated with projects such as PyFR [144] where low level modules target different hardware architectures and are linked with high level matrix operations. The OpenSBLI project [190] allows high level PDE specification, with parallel code generated for a range of architectures. Given the rapidly changing architectures of today, code modularity and maintainability is key. To release the potential of all computational power, low level alterations to assembly code can be required to surrounding libraries, which may not be sustainable. The coupling of multi-physics codes may offer ways to utilise different compute architectures more efficiently by tailoring codes and algorithms to different types of chip. Computational scientists therefore need to be involved throughout the entire code development cycle to help make bold choices early on and to further optimise in later stages.

Paths forward
In the near term, we should consider the benefits and disadvantages of each method and apply them where suited. Fig. 15 suggests where different methods may be applied indicated using FV, FD and FR. As the current workhorse, FV offers a great deal of flexibility and robustness. A key enabler is the increased use of hybrid LES-RANS to alleviate grid demands. Zonalisataion and treatment of transitional, accelerating and separated flows on low curvature surfaces, require attention to remedy deficiencies present in any RANS content. Methods to alleviate computational demands should be incorporated, such as highly adapted boundary layer meshing and improvements to time stepping procedures. Mixed geometrical fidelity offers a great range of flexibility in use, accuracy and computational cost.
Currently we may use high order with reduced cost for exploratory flows and flow physics and limit use of FV to industrial cases, where more modelling is used, leveraging the hierarchy in Fig. 3. Wave propagation for jets/fan acoustics seems an obvious use for higher order methods, potentially reducing mesh requirements for FV based LES. FR could also be used to explore numerical requirements in different zones, with these being implemented in similar more tractable schemes. The reduction in cost may also be useful for improving RANS modelling by studying simple geometries but where RANS is challenged by flow physics.
In the medium term, 4 th order appears to provide a significant increase in accuracy for minimal effort in mesh modification and coding. It is also less severe than > O 4 in regards to stability issues and memory requirements (particularly if implicit time stepping is employed). If FR is not used, there is a limited increase in stencil size and hence MPI data halo. When higher order methods mature, order can be increased, and this should be planned in current code iterations with increased modularity. For FR this should be fairly trivial.
Long term, the use of higher order methods largely depends on development of CFD as a whole. As noted, pre-processing is currently limiting the use of these methods along with practical limitations in the amount of RAM or storage available (i.e. cost). It is hoped time integration (or even parallelism) will see significant advances, which would counter or alleviate the number of (sequential) computed time steps required. Storage and post-processing of large simulation data is also becoming a bottle neck. Data encoding may drastically reduce storage requirements and allow greater details to be extracted with only a slight loss in fidelity. To extract knowledge from vast databases, data mining and machine learning, although young in CFD analysis, seem the only likely way to inform designs by reliable quantitative data.

Summary and outlook
The role and design of aero-engines is set to change rapidly with future aircraft. With this, highly coupled systems arise throughout at a range of scales. To model these requires great flexibility and mixed fidelity predictive simulation. Although current high fidelity methods have become mature and more widespread, FV methods in particular are more limited in accuracy but rich in application and user experience. Moving forward, higher order methods offer tantalising advantages. Great hurdles, particularly in the context of turbomachinery application, remain. However, great computational savings can still be made by optimising LES mesh and time stepping to obtain the required accuracy at reduced cost. The variety in modelled physics, and corresponding greater internal and external coupling, demands a degree of flexibility in the code base for use of different algorithms and hardware. This runs counter to current industrial practice in which a single code is developed over decades. However it seems certain that uptake of high fidelity methods will continue at a rapid pace. Over the next 10 years, these and future predictive methods will enable unexplored design space to be traversed with confidence.