Enhanced Delayed Detached Eddy Simulation for Cavities and Labyrinth seals

A range of popular hybrid Reynolds averaged Navier-Stokes - large eddy simulation (RANS-LES) methods are tested for a cavity and two labyrinth seal geometries using an in-house high-order computational ﬂuid dynamics (CFD) code and a commercial CFD code. The models include the standard Spalart-Allmaras (SA) and Menter SST versions of delayed detached eddy simulation (DDES) and the Menter scale adaptive simulation (SAS) model. A recently formulated, enhanced, variant of SA-DDES presented in the literature and a new variant using the Menter SST model are also investigated. The latter modify the original deﬁnition of the subgrid length-scale used in standard DDES based on local vorticity and strain. For all geometries, the meshes are considered to be hybrid RANS-LES adequate. Very low levels of resolved turbulence and quasi-two-dimensional ﬂow ﬁelds are observed for the standard DDES and SAS models even for the test cases here that contain obstacles, sharp edges and swirling ﬂow. Similar ﬁndings are observed for both the commercial and in-house high-order CFD codes. For the cavity simulations, when using standard DDES and SAS, there is a signiﬁcant under prediction of turbulent statistics compared with experimental measurements. The enhanced versions of DDES, on the other hand, show a signiﬁcant improvement and resolve turbulent content over a wide range of scales. Improved agreement with experimental measurements is also observed for proﬁles of the vertical velocity component. For the ﬁrst labyrinth seal geometry swirl velocities are more accurately captured by the enhanced DDES versions. For the second labyrinth seal geometry, the mass ﬂow coefﬁcient prediction using the enhanced models is signiﬁcantly improved (up to 30 %).

The use of eddy resolving methods in turbomachinery has the potential to dramatically increase the predictive accuracy in numerical modelling and it is becoming more widespread for the prediction of gas turbine flows in all areas of the engine. The secondary air system is one area where expected benefits are high given the highly complex flow physics driven by the highly variable geometry and components with rotating walls. These flows have large regions of recirculation, complex strains, are highly anisotropic and RANS models can perform poorly. RANS models can also demonstrate extreme sensitivity to model choice [1] [2]. It is common for ad-hoc corrections to be applied to cases of specific interest (see of example [2] [3]) to offer remedies for their inadequacy. Fully resolving flow features down to the smallest (Kolomogorov) scales (direct numerical simulation) is currently not feasible due to the high resolutions required and computational cost.
Alleviating this to some extent, large eddy simulation (LES) proposes to resolve the large energy containing scales within the flow and model the smallest eddies which display a much higher degree of isotropy. Further cost saving may be obtained by using hybrid RANS-LES approaches which can significantly alleviate grid demand. This is because attached flow regions (i.e. boundary layers), where significant grid resolution is required to capture streak-like flow structures, is modelled with RANS. For boundary layers RANS is reasonably well calibrated. LES is used for free-shear regions where it is expected that high gains in accuracy will be obtained. Such methods can save orders of magnitude [4] in grid requirements in many areas of turbomachines. One of the most popular hybrid RANS-LES methods is detached eddy simulation [5]. Since its inception in 1997, detached eddy simulation (DES) has received various modifications to help address deficiencies reported in the academic literature. One of these deficiencies, as encountered in the studies of Menter and Kuntz [6], is modeled stress depletion. This occurs in the grey area between the switchover from RANS to LES where the velocity fluctuation from the LES content is not sufficient to compensate the loss of modeled turbulent stresses. In severe cases of modeled stress depletion, an underestimation of skin friction can result, leading to early separation. Also it has recently been suggested by Deck [7], that the possible delay in formation of instabilities in mixing layers may be due to the advection of upstream eddy viscosity. This may be particularly prevalent in cavity and seal flows where large recirculation regions are observed. DES has been widely used in many areas of the gas turbine, see for example the compressor [8], combustor [9], turbine [10], cavities [1] turbine blade cooling [11] and jets [12]. The method provides seamless integration between the RANS and LES regions and is primarily grid controlled. DES is a non-zonal method and hence is advantageous because locations of separations do not need to be known in advance or flagged by the user, which for complex geometries can become somewhat of a difficult and arduous task. An additional benefit is the ease at which resolved regions can revert to RANS. Implicit (I)LES, for example, implies zero eddy viscosity. This can create problems for some RANS models (see for example [13]) where it became necessary to freeze an eddy viscosity field from a precursor RANS simulation.
Various formulations of DES exist (see for example [14]), but the original and most popular uses a Spalart-Allmaras GTP-19-1652, Jefferson-Loveday (SA) base RANS model. With this a new length scale d is defined asd = min(d,C DES ∆) where d is the nearest wall distance and ∆ is a filter width related to local cell dimensions. Boundary layers are covered by the original RANS formulation, and detached regions by LES. The filter width (∆) has traditionally been taken as the maximum local grid spacing of each coordinate direction for two reasons. The first is so that the model should return the RANS method in attached boundary layers and therefore represents the safest estimate to ensure that requirement [6]. This is to minimise the chance of modelled stress depletion which occurs when a mesh is not fine enough to support LES content after the switchover from RANS to LES (a common problem to most if not all hybrid RANS-LES methods). Secondly this filter definition is physically justified given that for isotropic eddies near the cut-off, the smallest dimensions of a strongly anisotropic grid cell do not contribute to the grid resolution capability [15]. Hence on strongly anisotropic grids the volume based filter definition (∆ vol = (∆x∆y∆z) 1/3 ) is not physically justified. The former filter definition has, however, contributed to the delay in development of shear layers and/or two-dimensional flow (see [7]).
To address the issue of slow shear layer development a recent variant of DDES is trialled in this work which has been formulated by Shur et al. [16]. It uses a modified definition of SGS lengthscale that is flow dependent. This follows development by several other authors on SGS length scales including [17] who propose a definition based on dissipation rate of isotropic turbulence on anisotropic grids and [18] who used a flow-dependent definition. The latter was based on the orientation of the local vorticity vector. In the early stages of a plane mixing layer, where the flow is highly two-dimensional, the vorticity vector is in the cross-stream direction. A plausible physical length scale was derived such that it became a function of the stream-wise and stream-normal mesh dimensions and in isotropic mesh regions reduced to the standard volume filter definition. The above has been re-formulated into a DDES context by Shur et al. [16] and is again based on the local orientation of the vorticity vector. However, [16] also use an additional term called the vortex tilt measure (VTM).
This is designed to detect areas of flow with largely two-dimensional flow structures, and maintain a value close to zero. In regions of fully developed turbulence it is designed to be close to one. The VTM term is introduced as a multiplication factor with the enhanced length scale described above. In this paper we also trial a new variant using the Menter SST turbulence model. The key to these methods is that in fully developed turbulence and on isotropic meshes the standard DDES model is retained where it is physically justified. In addition, a range of popular hybrid RANS-LES models are investigated for comparison and are described in the following sections.

TURBULENCE MODELS
For the standard Spalart-Allmaras (SA) based DDES model [19] the wall distance in the SA RANS model is replaced with the following DDES length scale: where ∆ = max(∆x, ∆y, ∆z) for standard DDES, C DES,SA = 0.65 and d is the nearest wall distance. The shielding function of the DDES formulation is given by: where C d1 = 8 and: where S is the magnitude of the strain rate tensor and Ω is the magnitude of the vorticity vector.

MSST-DDES
For Menter SST based DDES [20] the standard destruction term of the k equation is modified as follows: where the DDES length scale is given by: Here the function f d is the same as for SA based DDES with the exception that C d1 = 20. Again ∆ = max(∆x, ∆y, ∆z) for standard DDES and C DES,MEN = 0.61. The RANS length scale l MEN is given by: Here β * = 0.09 and: where F 1 is a blending function in the standard Menter SST model [21] and C DES,MEN1 = 0.78 and C DES,MEN2 = 0.61.

Enhanced DDES formulation
This method was recently presented by Shur et al. [16] using the SA turbulence model. Here, a new variant is also investigated using the Menter SST model. For both approaches, the subgrid length-scale definition ∆ is modified as: The quantity∆ Ω is a mesh based length scale sensitized to strongly anisotropic cells by taking into account the direction of the local vorticity vector. It reads:∆ where l n = n Ω × r n and n Ω is the unit vector of the direction of vorticity. Also r n (n = 1...8) are the locations of the cell vertices for a cell with the center vector r.
The V T M term provides a measure for the identification of early shear layer regions. It is based on the local strain tensor S and the vorticity vector Ω Ω Ω. It reads: where ν * = 0.2ν and tr() denotes trace. To ensure rapid change between regions the following piecewise-linear function is used:

SST-SAS
The transport equations for the SST-SAS [22] model differ from those of the SST RANS model by the addition of a SAS source term Q SAS in the transport equation for ω. It reads: where µ ω with c µ = 0.09. The von Kármán length scale L νK is given by: To control the damping of high wave numbers a limiter is introduced through analyzing the equilibrium eddy viscosity of the SST-SAS model. It reads: where ∆ vol = (∆x · ∆y · ∆z) 1/3 , C s = 0.1 and β and α are functions of the SST model.

SOLVER AND NUMERICS
Simulations are performed using the in-house multi-block solver 'DOLPHIN' [3] and the ANSYS FLUENT [23] CFD code (Version 18.2). DOLPHIN solves the flow on structured grids which can be curvilinear and parallelizes over blocks using MPI. For spatial discretization, a finite difference approach is used motivated by the relative ease of extension to highorder accuracy. A sixth order compact scheme is used for the inviscid terms and standard second-order central differences for the viscous terms. For stability, a tenth order implicit filtering scheme is employed for the conservative variables following Visbal [24]. For problems with extremely fine mesh resolution, such as in the cases here, the stability of the explicit timemarching scheme used in [3] was found to be too restrictive and very inefficient. Therefore, the approximately factored diagonal implicit scheme of Pulliam and Chaussee [25] is used. This has the advantage of solving five scalar systems of equations rather than a block system thus considerably saving computational work. To further reduce computational work viscous flux Jacobians are approximated to the viscous Jacobian spectral radius and included as a diagonal term on the implicit side. The physical time derivative is discretized using a second order backwards difference and the equations are GTP-19-1652, Jefferson-Loveday subiterated in pseudotime at each physical timestep until convergence is achieved. Full details of the implicit numerical scheme implemented in DOLPHIN may be found in [25].
For FLUENT based simulations, the cavity and first labyrinth seal geometries use the pressure based solver with the semi-implicit method for pressure linked equations (SIMPLE) algorithm. The second-order bounded central difference scheme is used for momentum terms. For time advancement the bounded second order implicit formulation is used. For the second labyrinth seal geometry the flow in the seal gap reaches a Mach number of (Ma = 0.5) and the density based solver was found to offer vastly superior convergence properties for this case. The third order MUSCL and implicit second order spatial and time discretization schemes were used. Full details of the numerical schemes may be found in [26].

CAVITY CASE SETUP
A cavity flow is investigated based on the experiments of Grace et al. [27] where detailed LASER Doppler velocimetry (LDV) and hot-wire measurements were made at a Reynolds number of 50, 000 based on the length of the cavity. The cavity was of length 0.1 m and depth 0.025 m giving a length-to-depth ratio of L/D = 4. The span-wise extent was the full width of the wind tunnel (0.9114 m). This resulted in a cavity that satisfied the conditions for two-dimensional mean flow within the central portion of the cavity specified by [28]. Flow measurements over a range of spanwise locations were also conducted to confirm this [27]. A schematic of the geometry considered is shown in Figure 1. In the current simulation the spanwise extent of the domain is three cavity widths (3L). The inlet velocity profile is specified using experimental measurements at an x location of (x/L = −0.1) upstream of the cavity leading edge. The inlet profile is laminar and the measured inlet turbulence intensity was less than 0.2%. Hence, transported turbulence quantities are obtained for a turbulent viscosity ratio (ν t /ν) of 1. In DOLPHIN, at the outlet, all flow variables are extrapolated using second order differences apart from the static pressure which is specified. At the inlet, the velocity, temperature and turbulent transport variables are specified and the pressure is extrapolated from interior nodes. In FLUENT the standard "velocity inlet" and "pressure outlet" boundary conditions are used. Towards the upper wall the mesh is rapidly expanded and a slip wall boundary condition is specified. is in the θ direction. Again the first off-wall grid node is located at y + < 1 and in terms of these values the mesh is close to wall-resolved requirements [30] and hybrid RANS-LES adequate [31].

LABYRINTH SEAL 2 CASE SETUP
Labyrinth seal 2 is based on the experiments of [32] and is an example of a straight-through labyrinth seal. The section of the seal is shown in Figure 3  to the grid resolution and the generation and growth of wake-like structures within the cavity shear layer as seen in Figure  4. This is due to the modified definition of filter length scale in these models. If two-dimensional flow regions are detected the model automatically adjusts the DDES length scale definition from the standard maximum local grid spacing to a more plausible definition of grid spacing based on the x and y directions where the vorticity vector is aligned with z. The VTM term is also close to zero in two-dimensional flow regions and close to one in fully developed 3D flow with random orientations of strain and vorticity. Hence if two-dimensional flow regions are detected, as observed for the standard DDES models, the length scale is appropriately adjusted.
Comparisons are made with experimentally measured mean flow streamlines [27] in Figure 5. Upstream of (x/L = 0.38) there is a stagnated flow zone where flow velocities are lower than the measurement uncertainty and have hence been omitted from the comparison [27]. The flow is characterized by a large recirculation region centered towards the rear of the cavity. using DOLPHIN and FLUENT. A key feature of this geometry is the extreme difference in local length scales between the tooth tips, where the passage is contracted, and the large cavity regions that reside in-between. For both codes, again, a stark difference is observed between the enhanced DDES formulations and the rest of the standard hybrid RANS-LES models, even for this highly featured geometry with significant swirl component. The former reveal highly complex turbulent flow fields containing intricate eddy flow structures with high levels of mixing. The seal teeth lead to the formation of instabilities which rapidly grow and generate highly 3-D turbulence. As the flow is accelerated through the seal gap high levels of turbulence are generated as vorticity is shed from the tooth tip. This flow impinges on the downstream step, the presence of which causes the flow to penetrate deep into the cavity. Examination of vorticity magnitude iso-surfaces for the standard DDES and SAS models revealed essentially quasi-two-dimensional fields. The rotating lower wall generates significant shear and production of eddy viscosity. It is thought that large recirculation regions formed within the cavities contribute to a reduction in resolved turbulent content through advection of eddy viscosity around the inside of cavities. Figure 9 shows profiles of swirl velocities at locations A and B and axial velocity at location A. The enhanced DDES formulations give the most accurate representations of swirl velocities across the seal gap at both locations whereas the other models over-predict the cavity swirl velocity by around 20% at location A. It should also be noted that both the enhanced DDES-SA and new enhanced DDES-MEN models perform comparably with one another at both locations. All models perform very similarly for the axial velocity predictions and general agreement is good. It is worth pointing out that in the experiments [29] bimodal velocity distributions were detected, and it is possible that significant errors exist in the experimental averages. In this case, as for the cavity case, very little unsteadiness was observed for the standard DDES methods and SAS. The enhanced DDES approaches on the other hand produce highly unsteady resolved turbulent content. Figure 10 shows instantaneous contours of vorticity magnitude for labyrinth seal 2 and the range of models tested for both codes. For compactness, only two of the cavities are shown in each case. The standard DDES and SAS models remained essentially steady and two-dimensional with very little unsteady resolved content. On the other hand, as can be seen in Figure 10, the enhanced DDES models produce significant resolved turbulent content over a range of length scales within the cavity zones. A wall-jet is generated from the tip clearance exit and drives a large recirculation zones within each cavity. The amount of shear generated by the jets is reduced due to the presence of the recirculating zones and hence energy is lost in these shear layers. Following this, the shear layers spread downstream to be accelerated into the next clearance gap.

LABYRINTH SEAL 2 RESULTS
A noticeable feature of this flow is the extreme acceleration of the fluid as it enters through the seal gap and it continues to accelerate until it reaches the next cavity. Following this it continues to develop rich turbulent structures. Table 1 gives percentage errors in flow coefficient compared with the experimental measurements of φ [32]. The flow coefficient is defined as: whereṁ is the mass flow rate through the seal, P in is the inlet pressure, A i is the clearance area and R is the universal gas constant.
The percentage errors are calculated using Equation 16. As can be seen, the standard DDES and SAS models significantly over-predict the mass flow coefficient by around 30 %. There is around a 5-10 % variation between codes for the standard DDES models and less than 1 % difference for the SAS model. The enhanced DDES-SA and DDES-MEN approaches show significant (up to 35 %) improvement compared to the standard DDES models and SAS models with 0.1 % and 5 % errors respectively. Interestingly, it is noted that the improvements observed with the enhanced models are due to the decreased flow coefficient indicating that there is insufficient blockage of the flow for the standard models. It should be noted that the enhanced simulations showed that the flow in the clearance region was essentially laminar (ν t /ν ≈ 2) with a very rapid transition to turbulent flow in the cavity regions. High levels of eddy viscosity in the wall-jet region for the standard DDES and SAS models damped out most resolved turbulent content. Although minimal, slightly higher levels of unsteadiness were observed in the standard DDES models using DOLPHIN with its high-order numerical schemes. This is thought to be the reason for the larger variation in flow coefficients observed between the standard models. There is some difference between standard DDES models between codes, which is attributed to different numerical schemes and level of artificial viscosity (lower in the case of the high-order scheme). It is interesting to note the higher flow coefficient for the higher order code, which could be due to a reduction in artificial viscosity in an essentially steady flow that is operating as neither RANS nor LES. The SAS model flow coefficients, between codes, are more comparable. This may be due to the SAS model operating in fully RANS mode, with higher levels of eddy viscosity and hence being less sensitive to code numerics.
It should be noted that for all of the standard approaches, the over-prediction of flow coefficient would indicate insufficient blockage provided by the modelled component. The flow here generates highly unstable shear layers which are ingested into the next seal clearance. This behaviour (including transition) is historically challenging for RANS models to cope with, and the enhanced models that do resolve turbulent content are able to respond in a much more physical way.

CONCLUSIONS
A range of popular and industrially available hybrid RANS-LES models have been tested for a cavity and two labyrinth seal geometries using a commercial CFD solver and an in-house high-order code. The models include Menter and Spalart For the cavity case, the standard DDES and SAS models produced essentially two-dimensional, steady-state results resulting in very low, if any resolved turbulence for both codes. The behavour is caused be excessive levels of eddy viscosity in the cavity shear layer and high turbulence production in the cavity recirculation region. The lack of resolved turbulence is reflected in the significant underprediction of streamwise turbulent velocity fluctuations and cross-stresses. In contrast, the enhanced DDES-SA and new enhanced DDES-MEN formulations lead to a rapid transition to fully 3D turbulence with a range of scales. Streamwise velocity profiles are reasonably predicted by all models but significant improvements are observed for cavity-normal velocity profiles using the enhanced DDES models. Comparison with mean flow streamlines indicate improved prediction of key flow structures compared with the standard DDES and SAS models. Menter based DDES models can yield lower levels of eddy viscosity when compared to equivilent steady RANS simulations which means that they are not operating as RANS or LES. It is recommended that hybrid RANS-LES models should be extensively tested for specific flow configurations and that special care is exercised by CFD practitioners when using many of the popular hybrid RANS-LES models that are currently available in commercial CFD packages.

ACKNOWLEDGEMENTS
The author acknowledges the computational resources provided by HPC Midlands Plus.