Impact of Underlying RANS Turbulence Models in Zonal Detached Eddy Simulation: Application to a Compressor Rotor

: The present study focuses on the impact of the underlying RANS turbulence model in the Zonal Detached Eddy Simulation (ZDES) method when used for secondary ﬂow prediction. This is carried out in light of three issues commonly investigated for hybrid RANS/LES methods (detection and protection of attached boundary layer, emergence, and growth of resolved turbulent ﬂuctuations and accurate prediction of separation front due to progressive adverse pressure gradient). The studied conﬁguration is the ﬁrst rotor of a high pressure compressor. Three different turbulence modelings (Spalart and Allmaras model (SA), Menter model with (SST) and without (BSL) shear stress correction) are assessed as ZDES underlying turbulence model and also as turbulence model of unsteady RANS simulations. Whatever the underlying turbulence model, the ZDES behaves well with respect to the ﬁrst two issues as the boundary layers appear effectively shielded and the RANS-to-LES switch is close downstream of trailing edges and separation fronts leading to a quick LES treatment of wakes and shear layers. Both tip leakage and corner ﬂows are strongly inﬂuenced by the Navier–Stokes resolution approach (unsteady RANS vs. ZDES) but the underlying turbulence modelling (SA vs. SST vs. BSL) impacts mainly the junction ﬂow near the hub for both approaches. ZDES underlying turbulence model choice appear essential since it leads to quite different corner ﬂow separation topologies and so to inversion of the downstream stagnation pressure radial gradient.


INTRODUCTION
Within a turbomachinery, secondary flows play a significant role on performances and surge margin. The tip leakage vortex flow is one of the main secondary flows and can be responsible for rotating instabilities, rotating stall or even surge. The corner separation is another main secondary flow. Thus, the accurate prediction of the kind of flows is essential for a turbomachinery manufacturer. As the flow is unsteady and vortical, and turbulence is anisotropic, Reynolds-Average Navier-Stokes (RANS) turbulence model such as Spalart-Allmaras or Menter SST ones, failed to predict accurately the tip leakage flow. On the contrary, Large Eddy Simulation (LES) is a good candidate. Nonetheless, as the Reynolds number is high and 90% of kinetic energy of turbulent fluctuations must be resolved, the cost is still too expensive for design purpose. Thus LES, especially Wall Resolved LES, are used for moderate Reynods number. Hybrid RANS/LES (HRLES) approach is a trade-off between accuracy and computational costs as the attached boundary layer is treated by RANS model while the separated boundary layer and wake are simulated with LES. This Emergence and growth of resolved turbulent fluctuations across the interface from RANS to LES regions, c Accurate prediction of the separation point due to progressive adverse pressure gradient, and not a geometric singularity. The present study aims at evaluating and confronting several RANS turbulence models and ZDES method [Deck, 2012] with the same underlying turbulence models: Spalart-Allmaras (named SA model hereafter, Spalart and Allmaras [1994]), Menter without Shear Stress Transport correction (named BSL model hereafter for BaSeLine version, Menter [1994]) and with SST correction (named SST model hereafter, Menter [1994]), in light of these issues. After explaining in detail the experimental facility, the numerical setup and the ZDES approach, the study focuses on flow near walls such as on tip leakage and corner flows. All simulation are compared with experimental data.

EXPERIMENTAL FACILITY
The simulated rotor is the first one of the 3.5 stage axial research compressor CREATE located at the Laboratory of Fluid Mechanics and Acoustics (LMFA) in Lyon, France and designed by SAFRAN AIRCRAFT ENGINES [Ottavy et al., 2012]. The Reynolds number based on the rotor chord at midspan is close to 800,000. The tip Mach number is equal to 0.92 and the rotational speed is 11,543 rpm. The casing diameter is 0.52 m. The meridian view shown in figure 1 highlights measurement planes between rotor and stator. These measurements are carried out with both pneumatic four-hole probes and unsteady stagnation pressure probes allowing the calculation of overall performances (pressure ratio, efficiency), stall line position included, throug azimuthal and radial average, and radial profiles via azimuthal average. The unsteady stagnation pressure probes developed by VKI have a very high bandwidth (250 kHz) and their diameter is 2.5 mm (2.65% of the spatial periodicity of the whole compressor at casing).

NUMERICAL TEST BENCH
In the present study, only the first rotor of the compressor (R1) is investigated at the nominal rotational speed and the nominal operating point. The inlet boundary is set to section 25A of the experimental test rig as measurements are available for different operating points at this section. Besides, the possible comparisons with experimental data at section 25A make the numerical setup more reliable for the evaluation of the method. The outlet boundary is defined at two axial chords downstream the R1 in order to have the experimental section 26A within the computational domain and to have a quite axisymmetric outflow as boundary condition is based on this hypothesis (radial equilibrium). The main simplifications concern the two axial gaps at the hub which are not included in the computational domain. Although these gaps are responsible for secondary flows, as this paper focuses on tip leakage flow and on behavior of turbulence models in the vicinity of the blade-hub junction, these simplifications can explain discrepancies to experimental data but do not modify the comparison between RANS and ZDES approaches. Computational grid is obtained with a multi block approach using an O4H topology. The rotor tip clearance is meshed with an additional O-H block. Two passages are simulated in order to fit the spatial periodicity of the IGV. The whole grid comprises 88 million points in order to model 1/32 of the compressor. The normalized wall cell dimension normal to the wall fulfils ∆y + of the order of 1 in every zone. In the vicinity of blade walls, ∆x + and ∆z + are respectively of the order of 200-300 and 100. Downstream of the blades, the mesh is progressively coarsened to 1700 for ∆x + and 150 for ∆z + to avoid numerical reflections. The normalized wall cell dimensions ∆x + and ∆z + are detailed in Riéra et al. [2016]. The IGV effects are reproduced with specific inlet boundary conditions based on 2D cartography resulting from a previous IGV-R1 computation. It consists in a rotating cartography of stagnation pressure, stagnation enthalpy, flow direction and turbulent variables, based on a Fourier decomposition with 60 harmonics of the two-dimensional map of the flow. Similarly to the phase-lag approach, for each radius, the decomposition in Fourier series is performed along azimuthal axis in order to build the inflow condition at each time step. For the outlet boundary condition, the static pressure is specified at the hub and then a simplified radial equilibrium law. A classic rotation periodicity condition is set at the azimuthal boundaries and a no-slip adiabatic wall condition is applied at all wall surfaces. More details on boundary conditions are given in Riéra et al. [2016]. The numerical simulations have been performed using the elsA software [Cambier et al., 2013], developed at ONERA and co-owned by AIRBUS, SAFRAN and ONERA. This software is based on a cell centered finite volume technique and structured multiblock meshes. The spatial discretization scheme for the inviscid fluxes [Mary et al., 2000] is based on the second order accurate Advection Upstream Splitting Method for low Mach numbers (AUSM+P) associated with third order flux reconstruction with MUSCL technique. The viscous fluxes are computed with a second-order centered scheme. For efficiency, an implicit time integration is employed to deal with the very small grid size encountered near the wall. The time discretization scheme is the second order accurate Gear scheme. At each time step, an approximate Newton method based on the LU factorization solves the non-linear problem. The time step is set to 1.6 10 −7 s which leads to a Courant-Friedrich-Levy number lower than 1 except for the boundary layers. It corresponds to 1000 time steps per IGV passing period. In the present study, all computations are performed with a unique numerical setup: mesh, scheme and boundary conditions.

ZONAL DETACHED EDDY SIMULATION
The hybrid RANS/LES method called Zonal Detached Eddy Simulation developed at ONERA [Deck, 2012] is a global approach inherited from Detached Eddy Simulation [Spalart, 1997] and originally formulated with the turbulence model of Spalart and Allmaras [1994]. In the computational domain, in order to reduce the computational cost, the boundary layer is solved using RANS equations while the LES is active far from the wall (the RANS model becomes a subgrid scale model (SGS)) and in separated boundary layer. The ZDES approach is used to investigate different flow problems relative to boundary layer separation which can be fixed by the geometry and massive (mode 1), due by progressive adverse pressure gradient and massive (mode 2), due by adverse pressure gradient and influenced by the dynamics of the incoming boundary layer (mode 3, not used in this study) as shown in figure 2. This hybrid RANS/LES approach is a modification of the underlying RANS turbulence model. Whatever the ZDES mode, the turbulent length scale appearing in turbulent transport equations is replaced by a hybrid length scale L DES . This length scale depends on the subgrid scale ∆ and on local flow properties (e.g. velocity gradient, vorticity). The idea is to increase the destruction of modeled turbulence (e.g. eddy viscosity) in order to obtain the levels produced by a SGS model. The substitution of this turbulent length scale is performed for some source terms of the transport equation of ξ turbulent variable (e.g. production P ξ and/or destruction D ξ ). It should be noted that ξ can be the pseudo eddy viscosityν (Spalart-Allmaras model), the turbulent kinetic energy k, the turbulent dissipation ε, the turbulent specific dissipation ω or another turbulent variable with transport equation. The modes of ZDES differ mainly by the choice of hybrid length scale L DES as shown in table 1 which summarizes the different length scales (L RANS , L LES , L DES ), different subgrid scales (∆ I DES and ∆ II DES ) and the boundary layer protection function (generic name f p ). In this table, the mode number corresponds to the one defined in figure 2 and all terms refered to the original transport equations of RANS model. It should be noticed that f w andS terms of SA model are relative to wall dumping functions. For all modes, the C DES calibration is based on homogeneous isotropic turbulence with free decay [Shur et al., 1999]. The original formulation of ZDES is based on Spalart-Allmaras turbulence model (named ZDES-SA hereafter). More details are given in Deck [2012]. In a previous work, Uribe et al. [2017] investigate the formulation of ZDES based on the k − ω SST turbulence model of Menter [1994], especially the substitution terms and the subgrid scale model. The formulation retained for the present study is the one of Strelets [2001] and is summarized in table 1. With k − ω turbulence models, the RANS length scale is no more the wall distance but the ratio √ k β ω . Due to the formulation of SST model, some terms are decomposed in two parts and the correct value is obtained by ponderation with F 1 function. Both BaSeLine and SST formulations are assessed in the present study (named ZDES-BSL and ZDES-SST hereafter, respectively). Figure 3 shows the radial profile of absolute deviation angle α and stagnation pressure p ia in the experimental plane 26A located downstream the first rotor. Five numerical simulation predictions are confronted to experimental data obtained with pneumatic probes and unsteady pressure probes [Ottavy et al., 2012]. Two unsteady RANS simulations are performed for comparison: RANS-SA and RANS-SST while three ZDES simulations are performed with SA, BSL and SST underlying turbulence models. As all simulations are unsteady, these radial profiles are obtained after a time-average over twenty IGV periods and then an azimuthal average. The comparison of numerical results to experimental data obtained using pneunmatic probes (4-hole probes) show a good agreement between them. Levels and gradients are well captured by CFD. The most signif- Table 1: Differences between SA and SST turbulence model in the ZDES framework (mode numbers are relative to those defined in figure 2). d w is wall distance and all terms refer to original RANS formulation (e.g. D k or Dν) SA model SST model

ROTOR PERFORMANCES
icant discrepancies are observed near the casing and the hub. Close to the hub, this is due to the absence of leakage flow caused by axial gap between fixed and rotating parts of the hub. This gap is responsible for a flow injection characterized by a strong azimuthal component which interacts with the main flow and modifies locally the deviation angle. Near the casing, the discrepancies are due to the different prediction of tip leakage flow, especially the tip leakage vortex and its behavior across the weak shock, and its role in the double leakage phenomenon, as shown later. All ZDES predict an underdeviation of a few degrees. The agreement of absolute stagnation pressure profile with experimental data is not as good as the absolute deviation angle. The phenomena of tip leakage flow, double leakage flow, thickening and separation of blade boundary layer, investigated in next sections, influence more on the stagnation pressure than the deviation angle. Nevertheless, discrepancies remain quite small (less than 1.5%). The thickening and then separation predicted by ZDES-BSL and ZDES-SST are responsible for stagnation pressure losses above 30% h/H relative height. Nonetheless, despite the loss overestimation, these simulations are able to well capture the radial gradient of stagnation pressure, especially the negative gradient between 30% h/H and 50% relative height. Near the hub, the more spread separated boundary layer predicted by RANS-SST and ZDES-SA lead to an overestimation of losses. It should be noticed that, due to the radial redistribution of massflow and load over the span, a small compression near the casing is linked to a high compression near the hub and vice-versa. Thus radial profile pass each other, global positive gradient for RANS-SA, RANS-SST and ZDES-SA and global negative gradient for ZDES-SST and ZDES-BSL. The analysis of radiale profiles shows that an overall good agreement is obtained between numerical results and experimental data and validates the simulation. In next sections,

TIP LEAKAGE FLOW
The first problematic mentioned in introduction a concerns the boundary layer protection which must be treated with RANS equation system in their whole thickness while the flow far from walls and separated flows must be solved with LES behavior. In order to visualize the kind of used equations (RANS or LES) and the grey zones with intermediate treatment, the λ sensor is defined as follows: This sensor is equal to zero in RANS regions, one in LES regions and ranges from zero to one in grey zones. The contours of λ sensor near the casing are plotted in figure 4 for the three ZDES. As expected, the flow far from walls is treated by LES part of ZDES approach. For the simulation ZDES-SA, the RANS region is very small from the first axial plane, as this is limited to the casing boundary layer and the tip region. The tip leakage vortex is treated by LES from its appearance, once the flow exits the tip gap. The coherent aspect of the tip leakage vortex in the first plane is due to the necessary time for eddy viscosity destruction. Thus the protection function f d works well for ZDES-SA. For both ZDES-SST and ZDES-BSL, the tip leakage vortex belongs also to the LES region (see plane 1). Nonetheless, RANS regions are more spread, beyond the tip flow and the casing boundary layer, especially in the vicinity of pressure side. This is due to (i) the protection function 1 − F F which is more conservative as the boundary layer thickness is overestimated and (ii) the behavior of L RANS depending on advected turbulent kinetic energy k and specific turbulent dissipation ω while for ZDES-SA this length scale depends only on wall distance. Whatever the discrepancies between ZDES approaches, they behave well as boundary layer are well protected, the flow within the tip (from the third plane) and tip leakage vortex belong to LES region. Nevertheless, the ZDES-BSL, the flow in the tip region is treated in LES only from the fourth plane because of the higher level of eddy viscosity due to the absence of SST correction and its impact on RANS length scale. The comparaison between the two blades highlights qualitatively the influence of external perturbations (IGV wakes here) on the protection function and on RANS and LES length scales. The approach ZDES-SA is clearly less sensitive than ZDES-SST and ZDES-BSL to these perturbations.  Figure 5 depicts the isosurface of Q criterion colored by normalized helicity and entropy variation field in 31% x/c plane, for the five unsteady simulations, at the instant t = nT (multiple of IGV period). Due to vortex dissipation, especially across the weak shock at 31% x/c, all unsteady RANS simulations are unable to capture any interaction between the tip leakage vortex and the tip flow of the adjacent blade. There is no double leakage flow phenomenon. On the contrary, all ZDES capture the double leakage flow even the amplitude depends noticeably on the underlying RANS turbulence model. With SA and SST models, ZDES predict a vortex breakdown in the vicinity of the adjacent blade, leading to numerous secondary vortices interacting with tip flow. Some of them are convected upstream to interact with the tip leakage vortex of the adjacent blade. For BSL model, as the eddy viscosity is higher, the dissipation is increased leading to a vortex breakdown located more upstream and secondary vortices are partially dissipated. Fewer secondary vortices interact with the tip leakage vortex. This analysis shows that the underlying turbulence model has a small impact on tip leakage flow by comparison to the switch from fully unsteady RANS to ZDES approaches.

HUB CORNER FLOW
The contours of the λ sensor close to the rotor blade junction with the hub is depicted in figure 6, at the instant t = T /2 and for the three ZDES simulations. Isolines of entropy variation colored by ZDES mode are also plotted to highlight the tip leakage vortex. Whatever the formulation, the attached boundary layer over the hub or the rotor blade is well protected for all ZDES (see planes 1 to 4). As shown by entropy variation isolines, the boundary layer separation at the hub corner provokes the switch to LES behavior of ZDES as wanted (planes 5 to 7). Nonetheless, for the BSL underlying turbulence model, the RANS region is more spread and the separation is treated with RANS equations (planes 5 and 6). This can be due to the absence of SST correction and higher eddy viscosity and to a too thin separation (it is recall that mode 2 is well suited for massive separation i.e. with a height significantly higher than the boundary layer thickness). When the separated region is sufficiently thick, the RANS behavior is replaced by the LES one (plane 7). Over the suction side, as the protection function is more conservative, RANS region are more spread with SST and BSL turbulence models. The comparison between the two blades shows that the ZDES-SA is quite insensitive to external perturbations while both SST and BSL are significantly impacted by these perturbations.
As all walls are adiabatic, the entropy variation is a good sensor of losses within a compressor. The entropy variation field in axial plane located at 91% x/c is depicted in figure 7 at the instant t = 3T /4. This field allows the visualization of the tip flow, the hub corner flow and IGV vortices observable near the left rotor blade. Near the tip, the dissipation of tip leakage vortex by RANS approach and the breakdown by the ZDES one are clearly observed. The comparison between the two rotor blades highlight the influence of IGV tip vortex on the tip leakage vortex, especially its convection and breakdown. As shown previously, the separation is much more radially spread for SST and BSL models. The boundary layer is thicker with k − ω turbulence model. Near the hub, the ZDES-SA separation is more massive than the RANS one and is more sensitive to the IGV vortex. The interaction between this vortex and the junction flow is responsible for a more massive separation. This sensitivity is also observed with SST model. Moreover the ZDES-SST  Figure 6: Contours of λ sensor in different axial planes near the rotor blade junction with the hub at the instant t = T /2, isolines of entropy variation (to highlight the tip leakage vortex) colored by ZDES mode and wall streamlines over rotor blade separation is larger than the RANS one. The entropy field shows the confinement of separation for ZDES-SA by comparison to ZDES-SST and ZDES-BSL. Near the hub, the separation predicted by ZDES-SST and ZDES-BSL is less azimuthally and more radially spread than the ZDES-SA one. The comparison between ZDES-SST and ZDES-BSL shows that the higher level of eddy viscosity in the attached boundary layer with BSL model is responsible for a RANS to LES switch more downstream which delays the chaotic degeneration and the turbulent fluctuations appearance in the separated flow.
The previous analysis is quite qualitative as only one instant is investigated and this is not representative of time-averaged flow field. The time-averaged field of entropy variation in several axial planes, 21.7%, 31%, 44.8%, 67.8%, 91% and 109% x/c, with c the tip axial chord is shown in figure 8. It should be noticed that, due to the time-averaged, the influence of external perturbations is not observable. As in instantaneous data, the tip leakage vortex predicted by RANS approach, for both SA and SST turbulence models, is clearly dissipated from the plane located at 31% x/c. Near the tip, losses are very similar between these two models (levels and high entropy area). Near the hub, losses are significantly higher for the SST model due to a larger separated region coming from the blade-hub junction. Moreover, the boundary layer is thicker for this model. As simulations are based on RANS approach, similar conclusions to instantaneous analysis were expected. For ZDES computations, despite the time average, the tip leakage vortex is clearly visible, especially its trajectory. This trajectory depends more on the approach (RANS vs. ZDES) than the turbulence model (SA vs. SST). With ZDES, the vortex is more oriented towards the adjacent blade which promotes the double leakage flow phenomenon. The absence of this phenomenon in RANS simulations allows the visualization of the tip leakage vortex footprint in the downstream plane (foreground plane) while, in ZDES, the vortex can not be distinguished from the rotor wake. Near the hub, the analysis performed for only one instant is still valid. With SA model, the ZDES amplifies the hub corner separation with respect to RANS approach. On the contrary, with SST model, the most massive separation is obtained with the RANS approach. Thus, the underlying turbulence model has a strong influence on the junction flow and on corner separation. The comparison of all simulations highlights a hierarchy between them in terms of massive hub corner separation: RANS-SA < ZDES-BSL < ZDES-SST < ZDES-SA < RANS-SST. This is consistent with the radial profile of absolute stagnation pressure ( figure 3(b)).

CONCLUSION
The present study focuses on the flow around a rotor blade of HP compressor and aims at evaluating the ZDES approach based on different underlying turbulence models (SA, SST, BSL) regarding three major problematics: boundary layer protection, emergence of resolved turbulent content and accurate prediction of separation due to progressive adverse pressure gradient. All The main discrepancies come from the absence of axial gaps. The two ZDES based on k − ω model improves the agreement in terms of radial gradient of stagnation pressure, especially the negative gradient far from endwalls. Near the tip, the impact of the approach is clearly shown on the tip leakage vortex, whatever the underlying turbulence model which has a small effect. The tip leakage vortex is treated as a LES region allowing emergence of resolved turbulent content. The boundary layers over tip, pressure and suction sides are well protected. Near the hub, the underlying turbulence model has a strong influence on the junction flow development, especially the corner separation. This is also impacted by the eddy viscosity level as shown by the comparison between ZDES-BSL and ZDES-SST. The approach has also a significant influence on the hub corner separation as a hierachy can be defined in terms of extent of separated region. The comparison between RANS and ZDES approaches show that only the ZDES is able to accurately predict the tip leakage vortex in the whole passage while RANS approach dissipates this one. This error is noticeable in the time-average field as RANS and ZDES predict different flows. Only the ZDES approach is well suited for the tip leakage flow prediction. The impact of underlying turbulence model is smaller near the tip. Near the junction, the underlying turbulence model has a strong influence which modifies the corner separation prediction. In the present study, with ZDES, the SA model predicts a larger separation than both SST and BSL models. This must be further investigated as the second mode of ZDES is dedicated to massive separation which is not experimentally observed. Subsequent to this investigation, one of the research field dedicated to ZDES approach is the ability to predict corner flows, especially separations. Several items will be investigated: (i) underlying turbulence model, with Differential Reynolds Stress Model for example, (ii) influence of resolved turbulent dynamics in incoming boundary layer, with the mode 3 of ZDES and (iii) the impact of leakage flow due to axial gap which influences the flow near the hub, especially the flow angles and stagnation pressure.