Large eddy simulation of a low-pressure turbine cascade with turbulent end wall boundary layers

We present results of implicit large eddy simulation (LES) and different Reynolds-averaged Navier-Stokes (RANS) models of the MTU 161 low pressure turbine at an exit Reynolds number of 90,000 and exit Mach number of 0.6. The LES results are based on a high order discontinuous Galerkin method and the RANS is computed using a classical finite-volume approach. The paper discusses the steps taken to create realistic inflow boundary conditions in terms of end wall boundary layer thickness and free stream turbulence intensity. This is achieved by tailoring the input distribution of total pressure and temperature, Reynolds stresses and turbulent length scale to a Fourier series based synthetic turbulence generator. With this procedure, excellent agreement with the experiment can be achieved in terms of blade loading at midspan and wake total pressure losses at midspan and over the channel height. Based on the validated setup, we focus on the discussion of secondary flow structures emerging due to the interaction of the incoming boundary layer and the turbine blade and compare the LES to two commonly used RANS models. Since we are able to create consistent setups for both LES and RANS, all discrepancies can be directly attributed to physical modelling problems. We show that both a linear eddy viscosity model and a differential Reynolds stress model coupled with a state-of-the-art correlation-based transition model fail, in this case, to predict the separation induced transition process around midspan. Moreover, their prediction of secondary flow losses leaves room for improvement as shown by a detailed discussion turbulence kinetic energy and anisotropy fields.


Introduction
Many large eddy simulation (LES) studies of turbomachinery flows focus on the low-pressure turbine (LPT) due to its low Reynolds number regime of 10 5 .To further reduce the simulation costs, the assumption of a statistically 2D flow at the midspan of a turbine blade is often made by applying periodic boundary conditions in the spanwise direction.However, a significant amount of the aerodynamic losses is generated in the secondary flow regions influenced by the interaction of end wall and blade boundary layers [1].Hence, the next logical step to evaluate the performance of the cascade is to conduct 3D simulations, including the endwall boundary layers, e.g.[2].
The MTU T161, considered in this work, is representative of high lift lowpressure turbine airfoils used in modern jet engines [3].The blades with a chord length of C = 0.069935 m and an average aspect ratio of 2.65 are staggered at an angle of 61.72 • .The cascade is arranged with pitch to chord ratio of t/C = 0.956.It features diverging end walls at an angle of 12 • , such that the flow cannot be studied using a simple spanwise periodic setup.Its geometry and boundary conditions are in the public domain and it has been the subject of both experimental [4] and numerical [5,6,7] investigations.The numerical investigations have focused on operating points with a Mach number of 0.6 and Reynolds numbers of 90,000 and 200,000 based on isentropic exit conditions.Müller-Schindewolffs et al. [5] performed a direct numerical simulation of a section of the profile in which the effect of the diverging end walls was modelled using inviscid walls (termed quasi 3D, Q3D).Recently, results obtained with a second order FV code were presented, which included the end wall boundary layers but the analysis was focussed on the flow physics in the mid-section [8].Various computations of the full 3D configuration were conducted using high-order codes during the EU project TILDA [6,7].However, due to the specification of laminar end wall boundary layers and no freestream turbulence at the inflow, no satisfactory results could be obtained.With a full 3D LES including appropriate turbulent end wall boundary layers and freestream turbulence obtained with a high order Discontinuous Galerkin (DG) method, we aim to provide a highquality reference dataset for this configuration.

Numerical method
We use DLR's flow solver for turbomachinery applications, TRACE, to perform an implicit LES with a kinetic-energy-preserving DG scheme for the spatial discretisation of the implicitly filtered Navier-Stokes equations [9].The scheme is based on the collocated nodal Discontinuous Galerkin Spectral Element Method (DGSEM) on Legendre-Gauss-Lobatto nodes with a polynomial order of 5.The anti-aliasing is performed by the split-formulation of Kennedy and Gruber [10], cf.[11].Due to the non-uniqueness of the solution at the element interfaces, Roe's approximate Riemann solver is applied for the advective part and the viscous terms are discretised by the Bassi-Rebay 1 scheme [12].To advance in time, a third-order explicit Runge-Kutta scheme of [13] has been used.Resolved turbulent scales are injected at the inflow boundary using a synthetic turbulence generation (STG) method based on randomized Fourier modes [14].The fluctuating velocity is introduced using a Riemann boundary condition [15] and, in contrast to the inner faces, no numerical flux function is used at the inflow faces.

Inflow boundary conditions
The operating conditions are specified by centerline inlet flow angle of α 1 = 41 • , total pressure of p t,1 = 11636 Pa and total temperature of T t,1 = 303.25 K as well as outlet Mach number Ma 2,th = 0.6 and Reynolds number Re 2,th = 90, 000 based on isentropic relations computed with the centerline outflow pressure.Furthermore it is required to match the momentum thickness of the incoming end wall boundary layers as well as the decay of freestream turbulence.The boundary conditions are then set by a spanwise varying profiles of p t,1 (z) and T t,1 (z) with a constant α 1 as well as spanwise varying Reynolds stress tensor u ′ i u ′ j (z) and turbulent length scale L T (z).All these profiles will be obtained by scaling a boundary layer profile in wall units at Re θ = 670 (https://www.mech.kth.se/˜pschlatt/DATA/, [16]) and combining it with freestream turbulence values.We perform a precursor channel flow simulation to determine the distance from the inlet required for the boundary layer to develop and meet the specifications.
Starting from the above centerline inflow conditions, we use the isentropic relations T t (z) to obtain a spanwise variation.Here, we assume a constant static temperature T 1 and pressure p 1 throughout the boundary layer.These two quantities can be computed from (1) by introducing a centerline Mach number Ma 1,center , which essentially controls the value of the T t and p t at the end wall.It has to be chosen such that the total pressure is always greater than the static pressure resulting from the simulation.Hence, an iteration might be necessary.It has to be noted, though, that small adaptations of Ma 1,center mainly influence the total pressure profile close to the wall and do not have a significant impact on the quantities of interest in the simulation.For this concrete case, we chose Ma 1,center = 0.362.The discrete DNS velocity profile u + i,DNS and the corresponding distances to the wall z + i,DNS are now transformed to obtain a Mach number profile Ma(z) simply using the definitions of u + and z + via and the DNS freestream velocity u + ∞,DNS .The speed of sound a 1 and viscosity ν 1 can be computed from T 1 and p 1 using ideal gas and Sutherland laws, re- spectively.This allows to describe the spanwise variation of T t and p t using (1).The non-dimensional turbulent stress tensor is scaled analogously via u Finally, the integral turbulent length scale in the boundary layer is estimated from the turbulent kinetic energy and dissipation rate as L + T = (k + ) 3/2 /ε + and dimensionalised as the z-coordinate in (2).
Optimally, the freestream Reynolds stresses and turbulent length scale should be chosen such that the measured turbulent decay is matched.However, in this case, this would result in a length scale in the order of cascade pitch.This could, in principle, be accommodated for by simulating more than one blade at the respective expense of computational effort.Instead, we chose to decrease the turbulent length scale and scale the non-zero components of the Reynolds stress tensor at the inflow of the domain, to reproduce the turbulence intensity in the blade leading edge plane according to the experiment.This leads to a stronger decay of turbulent structures and has to be considered when assessing the quality of the results.
The boundary layer profile and freestream values of the turbulence quantities are combined where they intersect at the edge of the boundary layer.Because of the large freestream turbulent length scale compared to the smaller length scale in the boundary layer, we use L T = max(L T,BL , L T,freestream ) for distances to the wall greater than δ 99 .Finally, the Reynolds stress tensor is rotated from the streamline- aligned to a Cartesian coordinate system.Fig. 1 (left) shows both the development of momentum thickness Reynolds number Re θ and freestream turbulence intensity Tu in comparison to fits to the measured data.As discussed, the turbulent decay is  steeper in our setup but the turbulence intensity at the blade leading edge (dashed vertical line) is well captured.

Validation and analysis of flow structures
The mesh consisting of 876,960 hexahedral elements (geometry polynomial order q = 4, solution polynomial order p = 5) with 189.4M degrees of freedom (DoF) for the final computation was designed to meet widely accepted criteria for wall resolution required for LES.This is demonstrated in Fig. 1 (right) which shows the streamwise, wall normal and spanwise solution point distances to be below 30, 1 and 30, respectively.The resolution in the free stream was ensured by the ratio of solution point distance and estimated Kolmogorov scale below 6 along a mid passage streamline.To be able to assess mesh dependence, results from a preliminary study with a significantly coarser mesh in spanwise direction with 312,424 elements (q = 4, p = 5) and 67.9M DoF are also included.After clearing the initial transient, the statistics on the fine mesh were sampled for 100 convective time units based on chord length and outflow velocity.Fig. 2 shows a comparison of the midspan blade surface pressure distribution and wake total pressure losses with the experiment.Both simulations agree very well within the 99% confidence intervals [17] for these first order statistics.While LES and experiment agree very well on the pressure side, a difference in surface pressure can be found on the suction side between x/C ax = 0.1 and 0.6.Similar offsets have been found in previous studies of this configuration, e.g.[5].The laminar separation bubble and subsequent turbulent reattachment indicated by the pressure plateau between 0.7 and 0.9 and recovery of the base pressure, on the other hand, is captured very well.The total pressure loss coefficient ζ is computed from time-averaged primitive variables in the wake plane at x/C ax = 1.4 and the upstream reference plane at x/C ax = −0.099.The pitch coordinate y is given with respect to the point y LE on the blade at its minimum axial coordinate x/C ax = 0.The simulation on the coarse mesh shows larger confidence intervals in the highly turbulent region of the wake due to its shorter averaging time of only 31 convective time units.Both simulations and the experimental data agree very well within the statistical confidence intervals.
With the simulation setup validated against the experiment, we can now discuss the system of secondary flows developing at the intersection of blade and end wall.It is shown in Fig. 3 (left and middle) in terms of time-averaged surface streaklines visualised with line integral convolution and vortices visualised using an isosurface of Q = 10 7 s −2 coloured with streamwise vorticity ω sw to indicate the sense of rota- tion.Limiting surface lines have been added by hand to ease the interpretation.The horseshoe vortex develops due to a roll-up of the incoming boundary layer.While its suction side leg (blue) dissipates well before the suction peak, its pressure side leg follows the passage cross flow towards the suction side of the next blade lifting off the end wall and becoming the passage vortex (PV).Behind the blade, a second large structure can be identified as the trailing shed vortex (TSV), which rotates in the opposite direction.Both vortices persist well downstream of the blade and become visible as the two regions of strong loss in Fig. 3 (right).The loss distribution in the wake plane also agrees very well with the measured data.Another small but very distinct vortex developing along the pressure side of the blade is the leading edge corner vortex.The pressure side features a short separation bubble close to the leading edge.Due to the diverging end walls, the backflow within the bubble is driven towards the endwalls where it rolls up, lifts off and mixes with the newly developing boundary layer in the passage between the passage vortex and the pressure side of the blade.

Conclusions
We have presented a new well-resolved dataset of the MTU T161 at Ma 2,th = 0.6, Re 2,th = 90, 000 and α 1 = 41 • with focus on appropriate reproduction of inflow turbulent boundary layers and freestream turbulence.Average blade loading and losses distribution due to secondary flow features agree well with the experiment underlining the validity of the presented approach.After the brief discussion of the secondary flow structures in this paper, more in depth analysis of the unsteady data set is ongoing.Future analyses will introduce a numerical experiment featuring unsteady wakes at engine-relevant Strouhal numbers.

Fig. 1
Fig.1Momentum thickness Reynolds number and turbulence intensity over axial distance from inlet plane (left) and average solution point distances in streamwise, wall normal and spanwise direction on the blade at midspan (right).

Fig. 3
Fig. 3 Mean secondary flow structures viewed towards end wall (left) and towards suction side (middle).Pressure loss coefficient at x/c ax = 1.4 (right) -LES contour colours vs. experiment contour lines.