Towards Large Eddy Simulation of gas turbine compressors

With increasing computing power, Large Eddy Simulation could be a useful simulation tool for gas turbine axial compressor design. This paper outlines a series of simulations performed on compressor geometries, ranging from a Controlled Diffusion Cascade stator blade to the periodic sector of a stage in a 3.5 stage axial compressor. The simulation results show that LES may offer advantages over traditional RANS methods when off-design conditions are considered – flow regimes where RANS models often fail to converge. The time-dependent nature of LES permits the resolution of transient flow structures, and can elucidate new mechanisms of vorticity generation on blade surfaces. It is shown that accurate LES is heavily reliant on both the near-wall mesh fidelity and the ability of the imposed inflow condition to recreate the conditions found in the reference experiment. For components embedded in a compressor this requires the generation of turbulence fluctuations at the inlet plane. A recycling method is developed that improves the quality of the flow in a single stage calculation of an axial compressor, and indicates that future developments in both the recycling technique and computing power will bring simulations of axial compressors within reach of industry in the


Introduction
In modern engineering practice an increasing emphasis is being placed on the role of Computational Fluid Dynamics (CFD) in the design process. In order for a CFD tool to be an effective part of the design environment it must provide highly accurate solutions to the flow problem in question within an acceptable timeframe -typically results are required within an overnight run. For gas turbine compressor components the CFD methods employed in industry have been centred around solution of Reynolds-Averaged Navier Stokes (RANS) equations. Whilst the capability of RANS to produce rapid flow solutions over a range of performance characteristics is highly desirable for industry, it is known to be deficient in key areas; first, RANS solutions often fail to converge at off-design conditions; the lack of confidence in the CFD solutions can restrict the design space. Second the ensembleaveraged RANS data does not give any information on the timedependent structure within the flow. Experimental evidence suggests that organised streamwise vorticity may play an important role in the evolution of the near-wall flow [1,2]. This is frequently referred to as Görtler vorticity and this array of streamwise vorticity is more easily observed on the concave pressure surface of a turbine. The absence of this structure from RANS solutions may lead to an inaccurate description of the boundary layer flow and therefore the loss mechanisms in the turbomachine.
In order to improve the understanding of the physical processes which lead to loss in gas turbines, a computational method which provides an accurate representation of the flow physics is required. Such a methodology demands both a high spatial resolution to capture the structure in the flow and a timeaccurate description of the physical processes. Direct Numerical Simulation (DNS), in which all of the scales of motion are resolved, produces the most accurate numerical representation of fluid flow phenomena, but this high level of accuracy is achieved at a high computational cost. The Reynolds numbers involved in practical gas turbine designs renders the use of DNS impractical in these flow configurations for many years to come, hence other time-dependent simulation techniques must be used. The two most viable alternatives are Large Eddy Simulation (LES) and Detached Eddy Simulation (DES). In the DES method, the boundary layer is modelled using a RANS method, and the larger scales of motion away from walls are solved explicitly using a LES-type approach. Whilst DES is viable for high Reynolds number external aerodynamics where the boundary layers are very thin compared to the geometry, its application to internal flows in turbomachinery is less beneficial. As was stated above the role of organised vorticity near to the wall appears to play an important role in the development of the boundary layer on blade surfaces, and DES will fail to resolve this structure as the boundary layer is assumed to be steady. In order to model these correlated fluctuations an interface treatment between the RANS and LES layers is required, which must be capable of generating correlated fluctuations that accurately model the organised vorticity. In Large Eddy Simulation, however, all scales of motion above a characteristic filter width are resolved, hence a well-refined grid near to the wall will capture the organised vorticity in the boundary layer. This paper, therefore, will exclusively consider numerical methods where the boundary layer is spatially wellresolved.
Early studies into flows relevant to turbomachinery using LES were restricted to simple geometries at modest Reynolds numbers. The laminar separation and transition to turbulence of a boundary layer have been simulated both on a flat plate [3] at a Reynolds number of 350,000 and on a geometry that includes a curved leading edge [4] at a Reynolds number of 3500 based on leading edge radius. In both cases the simulation results compared well with experiment and elucidated the transition mechanism in the flow. Wake-induced transition of a flat-plate boundary layer has been simulated numerically [5] at a Reynolds number of 150,000 with the transition mechanism captured in the simulation comparing very favourably with experiment. More recent studies into natural transition have shown that mode interaction plays an important role in the transition process [6], but the meshing requirements for these simulations mean that the simulation of natural transition will be restricted to simple geometries for the foreseeable future.
Published research which has focused on simulations of flows in turbomachinery geometries has generally been restricted to the mid-span of blades in linear cascades. Whilst computationally expensive, DNS has been performed for turbine cascades at transitional Reynolds numbers [7]. LES of turbine cascades have been performed at conditions that match reference DNS data [8], and it was found that the LES predicted a transition on the suction side of the blade due to impinging wakes that was delayed by some 10% when compared to the DNS. This was attributed to the  Cartesian components of position x þ ,y þ ,z þ distance in wall units coarse grid used in the simulation and highlights the resolution requirements for wall-bounded simulations. In all of the above simulations the incoming wakes were included in the simulations through the use of precursor calculations [9]. The effect of smallscale turbulence embedded within the wakes was studied by Wissink et al. [10], who concluded that whilst the large-scale motion of the wake triggers the Kelvin-Helmholtz (K-H) instability on the suction side of the blade, the small scale disturbances embedded into the wake flow seed the transition to turbulence in the shear layer flow.
Few published studies of compressor-type flows have been performed using LES. Early work by You et al. [11] established a computational framework within which Large Eddy Simulations of the loss mechanisms in rotor tip-clearance flows were investigated [12], including the effect of varying tip-gap size [13]. Based on the findings of this work, further efforts were made to improve turbulence modelling of the flow in turbomachinery [14]. The influence of freestream turbulence on transition on turbine blades has been investigated using LES [15], highlighting the need for well defined inlet conditions. Simulations of a linear compressor cascade at a moderate Reynolds number were studied by Zaki et al. [16]-the periodically passing wake caused the suction-side flow to separate and roll-up into discrete K-H vortices. On the pressure side the flow underwent a bypass transition due to the passing of the wake.
In this paper a series of Large Eddy Simulations of flows relevant to axial compressors are detailed. The simulations are performed to assess the applicability of the LES technique to flows that are of practical interest to axial gas turbine compressor design. The initial application of the LES method to compressor flows was performed through the simulation of an idealised linear cascade. This cascade, housed at the Naval Postgraduate School, Monterey, California (hereafter referred to as the Monterey Cascade) has been studied extensively and has a large database of results for a wide range of operating conditions. Simulations of a more realistic axial compressor rig were performed using the Cranfield low-speed cantilevered stator rig. These calculations were performed over a single stage of the machine, incorporating a rotor in the relative frame of reference, and a stator in the absolute frame. Experiments have recently been carried out on a research 3.5 stage compressor rig at the Whittle Laboratory, Cambridge University which represents a modern gas turbine compressor. The method is applied to this experiment to model a periodic sector containing three stators feeding into two rotors. This is to demonstrate the feasibility of conducting LES on a realistic modern compressor geometry. This paper is organised as follows: the numerical methods employed in the research code are outlined in Section 2. As a wide variety of flows were considered, each distinct geometry and the simulations performed on it are given their own self-contained section. The idealised linear compressor cascade simulations described in Section 3, the Cranfield cantilevered stator rig simulations are outlined in Section 4, and the Cambridge rig simulation is presented in Section 5. A discussion of the results in the context of the use of LES in an industrial environment is presented in Section 6. Concluding remarks are drawn in Section 7.

Numerical method
This work uses the Rolls-Royce CFD code Hydra [17] which is an unstructured, mixed element, compressible, density-based Navier-Stokes solver. The unstructured mixed element approach gives flexibility in grid generation to be able to handle complex geometry of industrial interest. The drawback is that, compared to a block-structured methodology, there is an overhead in computer memory and run-time, and unstructured spatial discretisation schemes are generally limited to second order spatial accuracy. The use of an unstructured industrial code that has been extended to handle Large Eddy Simulation means that the outcomes can be easily transferred to an industrial context.
The smoothing in the spatial discretisation has been reduced so as to avoid excessive dissipation of resolved eddies, and subgrid scale models incorporated. The important features are summarised below, and further details of the discretisation and testing on simpler LES flow problems can be found in Tristanto et al. [18]. The method has also been validated on multiple impinging jet problems [19] and the complex geometry capability demonstrated by the calculation of a Harrier aircraft in groundeffect [20].

Governing equations
The spatially filtered, Favre-averaged compressible Navier-Stokes equations can be expressed, for the conservative variables ðr,ru i ,EÞ, as @ @t and G(Q) contains viscous and conduction flux terms. The finite volume discretization provides an implicit filter for the large eddies and denotes unweighted filtered variables and~density weighted filtered variables. The finite volumes are created from the median-dual of the original unstructured mesh which may contain tetrahedra, hexahedra, pyramids and prisms.

Discretisation
The fluxes through the median dual control volume faces are accumulated to each node by looping over all the edges connecting the nodes. For an edge ij that connects nodes i and j, the flux is computed using the second-order accurate scheme of Moinier [21]: and the smoothing term is defined as [17] smoothing where L lp is the pseudo-Laplacian and For LES it is essential that the smoothing term should be kept as small as possible so as to avoid unphysical dissipation of the resolved eddies. Previous work by Ciardi et al. [22] on an earlier version of Hydra experimented with a 'wiggle detector' to automatically control the smoothing coefficient and tested this on the decay of isotropic turbulence and a fully developed channel flow. The self-adaptive scheme could reduce the smoothing coefficient to a value of 0.08 but would typically be higher when wiggles were detected. Similarly, Tristanto et al. [18] used the Ducros vorticity and divergence monitor [23] to control smoothing with results for a fully developed pipe flow and a Mach 0.9 free jet. In this case the Hydra code was used in comparison with a second unstructured method using Roe Flux Difference Splitting and linear-reconstruction with a least-squares approach gradient computation. The Ducros monitor includes a term which sets the minimum value of smoothing coefficient and for Tristanto et al.'s results this was set to 0.1, but again would be higher in many regions. In the present work, setting the smoothing coefficient to a fixed value of 0.2 was found to be more robust and reliable, and achieved the same low level of dissipation as the more complex Ducros monitor.
Temporal discretisation uses a standard third order accurate, three-stage Runge-Kutta algorithm.

Sub grid scale model
The standard Smagorinsky SGS model [24] defines the subgrid scale viscosity as where the Smagorinsky length scale is the instantaneous strain rate is and the filter width D is determined from the cube root of the median control volume surrounding a node and C s is a model constant.
A typical hexahedral mesh used for engineering calculations resolves the gradients in the boundary layer by decreasing spacing normal to the wall. This results in moderate to high aspect ratios close to the wall (hundreds or more). The consequence is that the volume based filter width stays relatively large and excessively high values of m sgs are found close to the wall. This could be avoided by using element aspect ratios that are significantly reduced to be closer to unity -but in practical terms the large increase in the number of elements would render the calculation impractical.
In the inner region of the Baldwin-Lomax or Cebeci-Smith mixing length RANS model, the turbulent viscosity takes a similar form to the subgrid scale viscosity of the Smagorinsky model: where the mixing length is in this case Since the Smagorinsky SGS model only accounts for the modelled part of the total stress, this should be smaller than the Reynolds stress predicted by a mixing length model. Hence, we have a way of restricting the excessive length scale in the Smagorinsky model: This can be interpreted as performing a similar operation to van Driest damping which reduces the length scale to zero near to the wall. Although the limiter employs a RANS-type model to compute an upper bound on the length scale near the wall, the subgrid-scale viscosity is calculated at every time step, and the simulation remains as an LES and is not LES coupled to a near wall RANS model. The Smagorinsky constant is set to C s ¼0.1 for all simulations and is typical of that used for wall bounded flows.

Parallel implementation
The complex geometry and the requirement to resolve small scales, necessarily leads to a large number of grid points (typically tens of millions or more). When coupled with the need to run for a large number of time steps (typically 500,000 or more), because of the disparity in time scales between the smallest and largest eddies, means that these calculations are unfeasible on current single processor machines and so must be run on large scale parallel facilities to achieve a reasonable turnaround time. This is a particularly important aspect if LES is to be feasible in a design environment.
Whilst structured multiblock CFD codes are relatively straightforward to implement in parallel using domain decomposition by block, the unstructured solver requires an efficient partitioning strategy and careful handling of the message passing to achieve good efficiency on large numbers of processors. The unstructured solver uses the OPLUS library [25] with message passing subsequently implemented in MPI. The partitioning is carried out in parallel using the ParMetis library. More information is provided by Hills [26] on how the parallel implementation has been tuned for large scale problems and near linear speed-up is demonstrated up to 1024 processors on an IBM Power5 system. Calculations presented here have been run on the HECToR Cray XT-4 system using up to 256 AMD Opteron processor cores, and the Loughborough University HPC facility using up 360 Intel Xeon processor cores.

Monterey cascade
Initial testing of the applicability of LES for compressor flows was carried out on the Naval Postgraduate School linear compressor cascade at Monterey. The rig models a Controlled Diffusion (CD) compressor cascade incorporating a stator blade of chord length c¼0.1273 m. The rig has been the subject of an extensive experimental testing programme over a wide range of flow inlet angles, b, and over a range of Reynolds numbers Re ¼ 400; 0002700; 000 based on the freestream velocity and chord length [27][28][29][30][31]. The current simulations were based on a series of experiments performed on the rig, the parameters of which are outlined in Table 1. The inflow angle b is defined as zero on the horizontal axis and is positive in a clockwise sense with the flow from left to right (see Fig. 1). Based on experimental data the design angle, where the minimum loss occurs, is b % 341. Cases 1-3, where the inflow angle is smaller than the design angle, are referred to as negative incidence simulations. Similarly, cases 6-9 are referred to as positive incidence simulations as the inlet flow angle is higher than the design angle. Earlier work published on this configuration can be found in McMullan and Page [40,44,46].
A two-dimensional plane through the computational domain is shown in Fig. 1. Hexahedral elements are generated using an industrial structured grid generator currently used for blade flows. The CFD solver views this as an unstructured set of hexahedral elements and does not take advantage of any implicit connectivity within the grid. The domain extends 0.25c upstream of the blade leading edge, and 0.75c downstream of the trailing edge. An outline of the mesh is included in the figure, with every other grid line shown for clarity. The blade surface is discretised with 370 nodes on both the suction and pressure surfaces, with grid clustering towards the leading and trailing edges. In nondimensional wall units, the axial spacing varies from x þ % 15250. The spacing normal to the wall is set so that the first cell height yields a minimum non-dimensional wall distance of y þ % 1:6. The spanwise spacing of the grid is set to produce z þ % 16. Two distinct domains were considered in this study, each having the same grid spacing in each direction, but differing in span. The first domain has a spanwise extent of L Z ¼ 0:2c, with 400 nodes placed uniformly across the span in order to satisfy the z þ criterion described above. This grid is denoted 'LES-W' and contains approximately 39 million nodes in the computational domain.
In order to assess the effect reducing the number of computational nodes by limiting the spanwise extent, a further computational domain is considered, denoted 'LES-N'; this grid has a spanwise domain extent of L Z ¼ 0:04c and has 80 nodes placed uniformly across the span in order to maintain a constant z þ between the domains. This leads to a total of 7.8 million nodes. Whilst the span of the domain is small and may lead to confinement issues, it is common practise in industry to use such streamtube representations of linear blades in RANS modelling, and it is therefore appealing to quantify the effect that the narrow span may have on the flow. Because of the large computational resource requirements of LES, it is difficult to carry out a formal mesh dependency analysis. The mesh spacing parameters are based on those used by other workers computing DNS and LES of turbomachinery [7,8,16] and preliminary calculations with coarser grids. The meshes used here have approximately half the number of nodes used in Direct Numerical Simulations of turbine cascades [7,16]. Whilst the authors believe that this is a well resolved simulation, some caution is needed when drawing conclusions from Large Eddy Simulations.
In all simulations a steady total pressure and temperature condition is applied at the inlet plane to provide the inlet boundary condition for the simulations. As the cascade did not feature transient wake passage, the background turbulence level remained constant at roughly 1.4%, and this is modelled in the present simulations through the use of pseudo-random white noise with a x / c   Gaussian distribution, which is superposed onto the inlet plane at each time step. As will be seen later, when applied to stage calculations this approach is insufficient and more realistic inflow turbulence required. However, for this relatively clean cascade inflow this small random perturbation is sufficient. A subsonic pressure outflow condition is specified at the outlet plane. As the Axial Velocity Density Ratio (AVDR) of the domain is set to 1.025 the use of periodic boundary conditions in the spanwise direction is not possible as the walls are not parallel; hence the spanwise boundaries of the domain are specified as inviscid walls. The AVDR is set to match that found in the experiment of Case 8, where b ¼ 401. The upper and lower boundaries of the domain are periodic.
The computational time step is 5:8 Â 10 À5 c=V ref , where V ref is the inlet velocity for case 8. Based on the freestream velocity and the smallest length scale present in the grid, the convective CFL number (that is ignoring the speed of sound) is 0.25. In order to avoid the high computational cost involved in propagating starting transient flow features through the domain the initial LES flow field in each case is obtained from a partially converged RANS solution. Calculations were typically run for 200,000 time steps before statistical samples were gathered every ten computational time steps over a further period of 300,000 time steps. For the LES-N grid, the simulations required approximately 24 wall hours on 512 AMD Opteron processor cores, whilst the LES-W simulations took 48 h when run on 1024 AMD Opteron processor cores.

Mean flow
Mean flow solutions of the blade are computed at mid-span and presented in a co-ordinate system that is local to the blade surface, with u t defined as the velocity component tangential to the blade. Instantaneous flow visualisations are presented in the global co-ordinate system with global velocity components. To be consistent with the experimental presentation, velocity data are normalised by the inlet velocity, V in , of the simulation.
The surface pressure coefficient, C p , is defined by where p in is the static pressure at the inlet boundary. Reynolds averaged solutions of Chen et al. [32], with a linear low-Reynolds number variant of the Launder and Sharma model [33] are used for comparison with the current simulation results where these are available.
Surface pressure distributions at mid-span for Cases 1, 6-8 are shown in Fig. 2. Good agreement is found between simulations on both LES grids and the experimental data. However, the pressure distribution on the suction surface between 0:0 o x=c o0:06 indicates that a separation bubble exists near the leading edge on the suction surface in both LES cases. The experimental data, however, suggests that the separation bubble is not as pronounced in the real flow.
Boundary reverse flow at this location. However, in the LES-N simulation the flow has reattached at this station. It is evident that LES-W produces the closest prediction to the experimental boundary layer data, with the LES-N data predicting a boundary layer that is somewhat thinner than the experimental data. The boundary layer RANS data of Chen et al. appears to have similar accuracy to the LES-N simulation.
The velocity profiles were processed to compute integral boundary layer properties, as shown in Fig. 4 for Case 7. The boundary layer displacement thickness on the suction surface is well-predicted for the LES-W domain, whilst the LES-N domain shows an under-prediction. As the boundary layer thickness at the trailing edge was reported as d ¼ 0:126c in the experiment, it is likely that the narrow domain with a spanwise extent of L z ¼ 0:04c will fail to capture the development of the boundary layer at this incidence. The LES underpredicts the displacement thickness towards the trailing edge, but is similar to trends reported by other workers [29,34]. High shape factors for the LES pressure surface (Fig. 4c) indicate that the simulations are retaining laminar flow, whereas in the experiment the flow underwent natural transition to turbulence at approximately 40% of chord. As the grid used here is not sufficiently well-refined to resolve the instability waves that precipitate the transition to turbulence this result is to be expected. Refinement of the grid to capture this feature would effectively turn this into a Direct Numerical Simulation and so be computationally impractical for a flow of this Reynolds number. The RANS data of Chen et al. highlights a limitation of RANS models-all of the scales of motion in the flow are assumed to be turbulent and hence the boundary layer on the pressure surface is turbulent over the entire blade surface. On the suction side the shape factor distributions computed from both LES presented here agree well with the experimental data.
The loss in the simulations is defined by where P in and P out are the averaged total pressures at the inlet and exit planes respectively. The computed loss coefficient for the simulations are shown in Fig. 5 and are compared with experimental data as well as the RANS simulations of Chen et al. A constant value for the AVDR was used in the current simulations, whereas the experimental data indicates a small variation of the AVDR with flow inlet angle. Near the design angle the LES-W simulations produce a reasonable estimate of the loss, whilst at high negative incidence the loss estimate is too large. At high positive incidence the computed loss from the LES-W grid is too low when compared to the experiment, a trend that is common with other numerical studies of the cascade [34,29]. Whilst the trend of the loss estimation from the LES-N simulations appear to be in good agreement with experiment, it should be emphasised that the flow in these simulations is compromised by confinement of the narrow span and the parasitic effects of the inviscid spanwise boundaries, which artificially thickens the boundary  layer and wake flow. The Large Eddy Simulations gives an improved prediction of loss as compared to RANS.

Flow visualisation
The mean flow from the LES-W and LES-N simulations show that the imposition of a narrow span has a significant impact on the accuracy of the flow prediction. To understand what causes these differences the instantaneous flow structure is analysed in detail. The Case 7 global axial velocity distribution close to the suction surface is shown in Fig. 6 for both grids. Both flow-fields have common features: the flow separates near the leading edge, characterised by the region of reverse flow between 0:0 o x=c t o0:08. It then reattaches to the blade surface and becoming turbulent. The reattachment creates axially-orientated streaks in the flow, indicating contra-rotating vortex pairs within the turbulent boundary layer on the convex surface. The streaks persist along the blade surface up to the trailing edge. The streaks are similar to those found in other numerical studies of separation bubbles on flat plates [4,35]. Whilst the two images are similar, only a small number of these streaks occur in the LES-N grid simulation (Fig. 6b), showing that the narrow domain may be restricting the development of the boundary layer. In comparison, LES-W grid visualisation in Fig. 6a shows several streaks across the span. Fig. 7 shows slices through the suction-side boundary layer near the trailing edge. For the LES-W simulation the boundary layer is dominated by the presence of 'mushroom-shaped' eruptions. The source of these eruptions correspond to the spanwise locations where the low-speed streaks are visible in Fig. 6a. Several scales of structure are also apparent, with small structures residing near the surface of the blade, with progressively larger structures dominating the much of the boundary layer. These flow patterns are suggestive of contra-rotating vortex pairs in the boundary layer and have been widely observed in both experimental and numerical investigations of flat-plate boundary layers. [35,36] For the narrower LES-N simulation, these flow structures are not visible. The larger structures in the boundary layer of the LES-W simulation have a spanwise scale of approximately 0.05cgreater than the spanwise extent of the LES-N simulations. It is therefore reasonable to conclude that the development of the boundary layer is adversely affected by the narrow span of the LES-N simulations when the boundary layer thickness approaches d ¼ 0:04c. To determine a minimum spanwise extent of a calculation, an estimate of the largest scales can be found from the maximum expected boundary layer thickness and the span should then be several times this scale in order to capture these structures.
Of particular interest in the study was a discovery made in the simulation of Case 1 on the LES-W grid. Fig. 8 shows the variation of turbulence intensity along the stagnation streamline upstream of the leading edge of the blade. The turbulence intensity increases dramatically just upstream of the blade, before approaching zero as the leading edge is approached. This type of behaviour has been observed in experiments [37] and in other simulations of the flow around an elliptical leading edge [38]. A possible mechanism for this dramatic ramping in turbulence intensity is interpreted through the use of a Q-criterion [39] isosurface in the region of the leading edge, shown in Fig. 9. There are toroidal structures upstream of the leading edge of the blade, which aperiodically shed from the stagnation region and are stretched over the leading edge of the blade, resulting in pairs of contra-rotating vortices. These streaks significantly disrupt the spanwise coherence of the Kelvin-Helmholtz vortices that emerge from the separation bubble, and substantially shortens the length of the separated region of flow. As the streaks form due to the  shedding of the toroidal structures from the stagnation region, their spanwise position is transient -a feature observed in other numerical simulations [38]. Therefore in a temporal average of the flow there would be no evidence of the existence of these structures, as they would be smeared out by the averaging process [40].

Cranfield BRR axial compressor
Owing to its geometric simplicity the Monterey cascade described above represents an idealised test case for Large Eddy Simulation. However, the simple geometry lacks many of the features that are present in real machines, such as hub-gaps in the stator, or upstream components such as a rotor. For LES to be used regularly in industrial practice, the simulation methodology must be capable of handling complex geometry features which may include cavities, shrouded stators and gaps between the blade tips and the machine end-walls. A further set of simulations were   performed on an axial flow compressor with representative real geometry features included in the computational domain. The rig in question is the Cranfield low-speed compressor rig [41], with the cantilevered stator build considered here. The Cranfield BRR compressor rig is a low-speed, 4-stage axial compressor with a 75:96 blade count per stage. Experiments were performed over a range of performance characteristics, with traverse and surface pressure data produced for what is denoted the peak efficiency of the machine.
For the LES calculations the component of interest in the rig is the 3rd stage stator blade. The Reynolds number of the stator flow, based on chord length, is Re % 180; 000. The computational domain of the stator is meshed with the same non-dimensional wall parameters and general node distribution as that of the Monterey Cascade simulations. 300 nodes are placed on each surface of the blade, and 600 nodes are distributed uniformly across the span.

Stage calculation
In a real compressor the flow at an arbitrary stage in the machine is very complex, being comprised of the wake from the upstream component, the remnants of the wakes that have propagated downstream from stage(s) further forward in the compressor, and background turbulence fluctuations embedded in what is nominally the freestream in the blade passage. A simple approach to model the flow at the component of interest is to incorporate more of the machine upstream in the computational domain. This strategy is tested here by introducing the upstream rotor (rotor 3) into the computational domain. A complication of introducing the upstream rotor into the BRR simulation is the 75:96 blade count of the stage. A full annulus meshed to the fidelity required for LES will result in a computational grid totalling some 3 billion nodes-a computational exercise that is beyond the capability of current resources. A lower integer blade count periodic sector is also not computationally feasible, and even an approximate 3:4 blade count would produce a mesh of 123 million nodes for a geometry that does not match the experiment. Therefore an approximation is made whereby the blade count is modified to 96 rotors and 96 stators per stage, permitting the simulation of one rotor and one stator. The rotor is geometrically scaled by a factor of 75/96 reducing the thickness and chord such that the solidity is equivalent to the original 75:96 blade count configuration.
The rotor domain is meshed to the same fidelity as the stator blade, resulting in a grid of approximately 36 million nodes. The time step of the simulation is Dt ¼ 4:85 Â 10 À8 s, and as the rotational speed of the machine is o ¼ 115:171 rad s À1 , 11,482 time steps are required to simulate one blade passage period. The initial flow field is obtained from a converged RANS solution (using a mixing plane interface) of the flow obtained at the peak efficiency of the machine. The flow is permitted to develop for 16 blade passage periods, in order to obtain a statistically stationary flow-field. Statistical samples are obtained over a further 80 blade passage periods, corresponding to ten flow-through times for the entire stage. A steady inlet condition, corresponding to the peak efficiency from experimental data, is imposed at rotor 3 inlet. No turbulent fluctuations are imposed onto the base inlet flow.
Surface pressure distributions near to the hub and casing are shown in Fig. 10. Near to the hub the prediction is reasonable, with the turbulent boundary layer propagating from the upstream rotor helping to keep the flow attached to the blade in this region. Near to the casing, however, the prediction of surface pressure has deteriorated significantly, with the distribution indicative of a flow that is close to stall. A plot of axial whirl angle at stator inlet and exit shows that the predicted flow angle near to the casing at stator 3 inlet is six degrees higher than that found in the experiment (Fig. 11). This increased angle of incidence of flow onto the stator blade results in the stall of the stator in the casing region. To further understand the cause of this discrepancy the flow is visualised through iso-surfaces of Q-criterion in Figs. 12 and 13. In Fig. 12 the suction surface of the rotor and pressure surface of the stator are shown, with the flow from right to left. The development of the flow over the rotor near to the hub results in a turbulent boundary layer which convects into the stator domain, and this turbulent boundary layer prevents the spurious separation of the flow on the stator blade in this region. The passage of the rotor wake into the stator frame of reference is also visible. Near to the casing, however, the laminar flow reaching the rotor leads to the formation of a tip-gap vortex with laminar upstream conditions. In Fig. 13   tip-vortex can be seen on the pressure surface of the rotor near the casing-the flow impinges on the blade and causes a large region of turbulent flow contaminating this region. This large, turbulent vortex propagates into the stator frame of reference and effectively stalls the blade near the casing, resulting in the poor mean predictions described above.
The imposition of a poor representation of the unsteady inlet condition at rotor 3 inlet has a parasitic influence on the flow at the casing endwall, resulting in very poor predictions of the flow in this area. A simple solution to this problem would be to simulate the whole machine from its inlet, where the ambient conditions will be largely time-independent, through to the component of interest. Whilst superficially attractive the meshing and simulation run-time demands of LES renders this option intractable, particularly for machines with non-integer blade counts. Therefore, it is more practical to develop a means of producing an unsteady inflow condition that produces turbulent flow at the inlet plane to the simulation.

Recycling boundary condition
The imposition of a time-dependent, turbulent inflow boundary condition is one of the most challenging aspects of Large Eddy Simulation. As has been shown above the imposition of boundary conditions that represent the mean conditions of the reference experiment are inadequate for LES. For flows where a high degree of unsteadiness is expected at the inlet plane it is necessary to impose a transient boundary condition which represents the instantaneous turbulent flow environment, and produces averaged statistical information that is a good approximation to the mean experimental data.
Several turbulent inflow condition generation techniques have been developed, ranging from digital-filtering approaches [42] to recycling/rescaling methods [43]. In digital filtering methods, target mean velocity profile and second order velocity statistics are imposed, along with a prescribed turbulence length scale. A convolution method then produces a series of time-dependent flow data with statistical information that matches the prescribed profiles. In a compressor flow, however, it is not clear if one prescribed length scale provides an accurate description of the turbulence in all regions of the flow, such as the wake from the upstream rotor, the endwall boundary layers, and the background fluctuations in the freestream. In addition, the unstructured nature of the CFD code renders the efficient computation of correlated fluctuations at the inlet plane a challenging task. In recycling methods turbulent fluctuations are extracted from a sampling plane within the computational domain, imposed onto a mean flow distribution upstream in the domain, and the fluctuations rescaled to fit a target statistical distribution. This method produces correlated turbulence fluctuations from within the computational domain at a low additional computational cost. This method is adopted in the current simulation.
The unstructured nature of the research code means that the arbitrary specification of a sampling plane in the computational domain is non-trivial, hence in the current simulation the sampling plane is placed at the exit boundary of the computational domain. A schematic of the recycling scheme is shown in Fig. 14.
As the fluctuations extracted from the exit boundary are in the absolute frame of reference, it is numerically trivial to superpose fluctuations onto a base mean flow in the same frame of reference. Therefore, a small stub domain which corresponds to the region downstream of the trailing edge of stator 2 is placed upstream of rotor 3.
In order to extract fluctuations from the sampling a plane, a low-pass filter is applied: The smoothing factor a is a small number chosen so that frequencies due to turbulence fluctuations are removed, whilst frequencies similar to the blade passing frequency are retained.  For the computational time steps used here and a required cut-off frequency off twice the blade passing frequency this determines a smoothing factor of approximately 2 Â 10 À4 for the BRR cases and 1:5 Â 10 À5 for the Cambridge case. The filter is applied at every time step and the low pass filtered running mean, Q m , is accumulated. The instantaneous fluctuations at the exit plane are then extracted from the instantaneous solution by These extracted fluctuations are then added onto a mean inlet flow distribution at every time step. Owing to the lack of available experimental inflow data, the base mean flow data is obtained from a RANS solution of the compressor at peak efficiency flow conditions. No experimental information is available concerning the r.m.s. statistics at the inlet plane to the simulation. However, the repeating stage nature of the BRR rig suggests that the fluctuations at one location in a given stage in the machine should be statistically equivalent to the same point in another stage, hence no rescaling of the extracted fluctuations is performed in this study.
The instantaneous velocity field at the inlet plane is calculated by the addition of the recycled fluctuation, u 0 onto the base mean flow variable, u for each individual for each velocity component: Given that the code is compressible, the instantaneous velocity fields are transformed into total pressure, total temperature and flow angles, from which the characteristic boundary conditions are computed. As the fluctuations that are extracted from the flow will be correlated both temporally and spatially, this technique generates realistic turbulence at the inlet plane of the simulation. Fig. 15 shows the axial velocity distribution at an arbitrary instant in time from the LES. When compared to the base inlet flow in Fig. 16, it is evident that the superposition of the fluctuations significantly perturbs the base inflow.
The stub domain is meshed to the same fidelity as that of stator 3, and adds a further 1.9 million nodes to the computational domain. All other simulation parameters remain unchanged with respect to the single stage calculation described in Section 4.1. It should be noted that the RANS mean inlet flow distribution for the recycling simulation is not an exact match to the mean experimental peak efficiency data. This, combined with the change in rotor blade count, means that the characteristic of the stage has changed slightly and the simulation is not an exact representation of the machine. In the description of the results that follows the experimental data is shown as a guide, with the main comparison drawn between the single stage and recycling LES calculations.
Surface pressure distributions near to the hub and case are shown in Fig. 17. At 5% height the surface prediction from the recycling simulation is very similar to that of the single stage calculation, indicating that the boundary layer flow at the hub endwall in stator 3 is turbulent in both simulations. Near to the casing however, the recycling simulation shows a significant improvement over the single stage calculation with a steady inlet condition. A suction peak is now present in the pressure distribution, which shows that the flow development in this region is now approaching what might be expected for the peak efficiency conditions that are being simulated. Instantaneous visualisations of the vortex structure in the recycling simulations are shown in Figs  show that the instantaneous flow development on the rotor pressure surface is widely different in both simulations; in the stage calculation with a steady inlet condition the tip vortex impinging on the blade surface near the casing produces a large region of turbulent flow, whereas the recycling simulation produces a tip vortex with a reduced development angle (grey line) which does not impinge so far down the rotor pressure surface. As the oncoming casing endwall boundary layer contains some level of turbulence fluctuations, the violent tip vortex transition is suppressed, resulting in an improved prediction of the rotor flow near the casing. This has the consequent effect of improving the stator flow-field near the casing, so improving the surface pressure distribution in this region. The axial whirl angle at stator 3 inlet and exit are shown in Fig. 20. The incidence angle at stator 3 inlet in the recycling simulation is now markedly lower than that of the stage calculation with steady inlet, and is in much better agreement with the experimental data. This data confirms that the recycling technique improves the flow prediction in stator 3 by improving the onset flow from the upstream rotor.

Cambridge axial compressor rig
The final strand of research involved the simulation of a full periodic sector of a stage in a new 3.5 stage axial compressor, housed at the Whittle Laboratory, Cambridge University. The rig is designed to run at low Mach numbers, with the Reynolds number of the flow, based on chord length, being Re¼350,000 for the stator blades. The rig was conceived with numerical simulation in mind, and the blade count of the machine has been set to permit low blade count periodic sectors. Each stage has 50 rotors and 75 stators, hence a periodic sector of a stage can be modelled through the use of two rotors and three stators.
The focus of this numerical study is the simulation of the research components in the machine, namely stator 2 and rotor 3. As the recycling boundary condition described above is also incorporated in the simulation of the Cambridge Rig, the computational domain contains the five blades which comprise a periodic sector and a small stub domain upstream of stator 2, which is in the relative frame of reference. As was mentioned above this permits the transfer of fluctuations from the recycling  sampling plane onto the inlet plane whilst keeping the flow at both planes in the same frame of reference.
The blades are meshed to a fidelity that is reasonable for Large Eddy Simulation, with y þ % 1, x þ % 15À50, and z þ % 32 on both the rotor and the stator. The resultant mesh for each stator totals 34.6 million nodes, while each rotor mesh totals 24.6 million nodes. The inclusion of the stub domain produces a mesh which has a total of 160 million nodes. A schematic of the geometry used in the simulation is shown in Fig. 21.
The time step for the simulation is set to Dt ¼ 2:45 Â 10 À8 s based on the same criterion used for the earlier calculations. In order for the rotors to pass completely through the periodic sector, it is necessary to compute 122,000 time steps. A notional fluid element requires 300,000 time steps to travel through the entire computational domain. The simulation is performed over 60 Intel Xeon hex-core processors, giving a total of 360 processor cores. A total of 190 wall-hours are required to simulate the periodic sector passage, and 475 h are required to simulate a flow-through of the domain. As up to 10 flow-through times must be simulated in order to produce converged flow statistics, it is apparent that this simulation represents a very challenging test case for Large Eddy Simulation.
Representative instantaneous iso-surfaces of vorticity magnitude obtained from the simulation are shown in Figs. 22 and 23. The value of the iso-surface is chosen such that the boundary layer flow can be seen, thus excluding the flow structure in the wake regions. The flow around the stator leading edge (Fig. 22) displays evidence of axially orientated streaks in what is ostensibly laminar flow on the convex stator suction surface. Further downstream the blade it can be seen that the streaks become more prominent, as the flow undergoes a transition to turbulence. The appearance of these streaks may be Klebanoff modes, which appear owing to the interaction of small-scale turbulence with the boundary layer flow. These small scale fluctuations originate in the wake upstream of the stator, which is supplied by the imposed recycling inlet condition.
The flow on the pressure surface of stator 2 and the suction surface of rotor 3 is shown in Fig. 23. The concave stator pressure surface displays evidence of Görtler vorticity, which arises due to the concavity of the pressure side. This vorticity manifests itself as axially orientated streaks which indicate the presence of contrarotating vortex pairs in the boundary layer flow. The rotor suction surface displays evidence of streaks in the laminar region of the flow near to the leading edge, in a similar manner to that described above for the stator suction surface.
The elucidation of these flow features in real axial turbomachine geometries using LES points to the potential of the simulation method to improve understanding of the flow in these machines. The extremely intensive computational demands of LES, however, mean that simulations of the type presented here will be limited to academic study for some time in the future. A systematic study of a machine over its entire performance characteristic will require the next generation of massively parallel computers that will arrive within three to five years.

Discussion
The simulations presented here cover many aspects relevant to the accurate simulation of turbomachinery flows, and it is pertinent to consider some of the outcomes in the broader context of time-dependent flow simulation and its applicability to problems of practical engineering interest.
The Monterey Controlled Diffusion Cascade, with its simple geometry and benign inlet conditions, provides an ideal test case for the capabilities of LES to predict flows in turbomachinery. The prediction of global parameters such as surface pressure are readily captured by LES if the grid is of sufficient refinement to prevent spurious laminar separation [44]. RANS models perform equally as well at on-design conditions at predicting global parameters as the models have been tuned for many years to produce excellent results. At off-design conditions, where separated flow and transition can be expected to occur, RANS performs less well and the current LES shows that there is some potential for LES to offer improvements in flow prediction at off-design conditions. In an industrial context, however, the use of LES cannot be considered with the same approach to meshing and domain restriction as can be afforded with RANS solution techniques. As LES is a time-dependent simulation of the flow, the computation should be considered as a numerical experiment, and sufficient care must be taken to ensure that the computational domain and boundary conditions are an accurate representation of the real flow. Given that the mesh  requirements for wall-bounded LES are particularly stringent it is attractive to reduce the computational cost of the simulation by reducing the spanwise domain length. Whilst superficially attractive, the results presented here show that a small spanwise domain length confines the flow and results in inaccurate flow development, compromising the loss estimation in the cascade. For linear cascades it is important that the boundary layer is not confined, hence the spanwise domain extent should be greater than the maximum boundary layer thickness expected in the flow.
The Monterey simulations have also shown that the capability of LES to predict boundary layer transition is dependent on the method by which the boundary layer undergoes transition. Natural transition, in which Tollmien-Schlichting (T-S) waves precipitate the transition, is extremely difficult to predict using LES as the nearwall grid spacing required to capture the development of the T-S waves results in a mesh that would be comparable to DNS. As the cases simulated here have Reynolds numbers of approximately 700,000 producing a grid that could capture these waves would result in a simulation run-time that is well beyond current computational resources. However, boundary layer transition that is precipitated by a laminar separation of the flow can be successfully captured by LES, and has been observed in all of the off-design Monterey simulations. LES is particularly suited to the simulation of free-shear layer type transition [45] where fundamental K-H vortices are inviscidly unstable, and similar conditions can be expected on the top of a laminar separation bubble. The shear layer transition promotes a turbulent reattachment of the boundary layer, and the LES presented here successfully captures these events, producing near-wall flow structure that is very similar to those observed in flat-plate boundary DNS [35].  Numerical simulation of real axial compressor components is currently an extremely challenging task. As described above, the mesh must be well-refined near to the wall, resulting in computational grids that are an order of magnitude larger than RANS meshes for the same component. For machines that have noninteger rotor/stator blade counts the simulation of a full annulus of a stage is presently well beyond current computing power, and a full compressor LES will not be considered feasible for some years to come. For compressors which do permit the simulation of a full periodic sector with a small number of blades, it is likely that these simulations will be performed on a stage that is embedded within the machine. As stated above LES should be considered as a numerical experiment, where the boundary conditions are an accurate representation of the flow at the location of interest. For LES of a stage embedded in a compressor it is necessary to generate a time dependent inlet condition that captures the passage of the wakes from the upstream components, and the background turbulence in the freestream. In order to produce an inflow condition of this type, sufficient experimental data on which to base the magnitude of the fluctuations is highly desirable, which includes the phase-locked mean velocity components and the second moment turbulence statistics in a two-dimensional distribution at the location of interest of the simulation. Experimentally, it is an extremely challenging task to obtain such an extensive set of flow statistics, as the volume of information that must be recorded to produce the relevant information is very large indeed. However, for accurate numerical simulation of compressor flows to be performed in the future this type of dataset is essential, and experimentalists must bear these issues in mind if the experiment is under consideration as the subject of a numerical investigation.
The Cranfield BRR calculations presented here are one of the first efforts to compute the flow within an axial compressor using a pure LES method. It has been shown that the inflow boundary condition is of critical importance, and that the time-dependent inlet condition outlined above is necessary to produce reasonable computational results. The recycling method here is a simple framework in which correlated turbulence fluctuations are added to the inlet boundary, and further development of this method will no doubt significantly improve the flow prediction.
The Cambridge rig, specifically designed to be amenable to numerical simulation through its integer blade count, has shown that the simulation of periodic sectors of a realistic axial compressor is within reach. The mesh resolution required for LES results in a computational cost that is currently at the upper limit of what can be achieved with current computing resources, but it can be expected that simulations of such a magnitude will become more commonplace in the near future. The simulation has shown that many interesting flow features are captured by the LES, and that the flow in axial compressors is much more complex than the authors had previously assumed. An improved understanding of the physical processes occurring in the flow could help to improve compressor design, resulting in more efficient gas turbines.

Concluding remarks
Large Eddy Simulations of a series of flow configurations relevant to axial compressors have been performed. The objective of the study was to assess the capability of LES to predict the features of axial compressors, and if the simulation technique could offer advantages over the RANS methodologies commonly used in the industrial design process. Simulations of an idealised controlled diffusion cascade have shown that LES can produce results which offer some advantages over RANS, and that the LES captures separation-induced transition in the flow. The capability of LES to elucidate new physical processes has been shown through the discovery of toroidal vortices in the stagnation region of the flow, which form contra-rotating vortex pairs around the leading edge. The computational domain requirements for LES have been demonstrated, as narrow span simulations effectively confine the flow and consequently produce poor simulation results.
Simulations of real axial compressors have highlighted the computational complexity of performing LES of real machines. In particular, an accurate description of the unsteady flow at the inlet plane of the simulation is essential in order to produce reasonable simulation results. It is recommended that experimental measurements are made with the demands of numerical simulation in mind, as an extensive set of turbulence statistics is required in order to produce accurate inflow conditions for LES. A simulation of a full periodic sector of a stage in an axial compressor has shown that such simulations are feasible with current computational resources.