Hybrid RANS-LES of Shaped Hole Film Cooling on an Adiabatic Flat Plate at Low Reynolds Number

Hybrid RANS-LES methods offer a means of reducing computational cost and setup time to simulate transitional flows. Several methods are evaluated in ANSYS CFX, including Scale-Adaptive Simulation (SAS), Shielded Detached Eddy Simulation (SDES), Stress-Blended Eddy Simulation (SBES), and Zonal Large Eddy Simulation (ZLES), along with a no-model laminar simulation. Each is used to simulate an adiabatic flat plate film cooling experiment of a shaped hole at low Reynolds number. Adiabatic effectiveness is calculated for Blowing Ratio (BR) = 1.5 and Density Ratio (DR) = 1.5. The ZLES method and laminar simulation most accurately match experimental lateral-average adiabatic effectiveness along the streamwise direction from the trailing edge of the hole to 35 hole diameters downstream of the hole (X/D = 0 to X/D = 35), with RMS deviations of 5.1% and 4.2%, and maximum deviations of 8% and 11%, respectively. The accuracy of these models is attributed to the resolution of turbulent structures in not only the mixing region but in the upstream boundary layer as well, where the other methods utilize RANS and do not switch to LES.


Introduction
Modern gas turbine engine performance is driven through technologies enabling higher turbine inlet temperatures. One critical enabling technology for increasing turbine temperatures is the ability to cool the surfaces of rotating and statio-nary components in the turbine section by injecting a film of lower temperature air over the surface to shield it from the extreme gas temperatures exiting the combustor. With the advantages of film cooling in gas turbine engines, there comes the added complexity of the fluid dynamics and heat transfer involved.
Film cooling flows are complex and boundary-layer dominated, and injecting film cooling into high-pressure turbine sections requires extensive computational resources to accurately predict. Substantial research has been invested in better understanding and optimizing film cooling flows [1] [2]. Shaped film cooling holes are often used [3], which exhibit an improved spreading of the film over the surface while also keeping the film attached to the surface. Despite the breadth of work studying shaped hole film cooling, there is a lack of published research that demonstrates accurate CFD of the surface temperature and effectiveness. This is in no small part due to the highly turbulent and complex flow of shaped hole film cooling. Although limited research has been able to predict adiabatic effectiveness along the hole centerline, conventional Reynolds-averaged Navier-Stokes (RANS) turbulence models have been shown to be inadequate in capturing the lateral spreading of the jet in crossflow [4] [5] [6]. Recent work has also explored the use of a Lattice-Boltzmann solver to capture the turbulent flow physics for both adiabatic and conjugate heat transfer cases [7] [8]. With all things considered, some form of a scale-resolving simulation with a subgrid-scale model is desired in order to accurately predict the effectiveness of film cooling schemes.
Subgrid-scale turbulence models are needed when small-scale physics, smaller than what can be resolved by the computational mesh, are important to the overall solution [9]. Large eddy simulation (LES) employs subgrid-scale modeling and is commonly used in such situations. However, a proper LES is computationally more expensive than unsteady RANS (URANS), leading to LES being relatively impractical in engineering applications. This creates a need for hybrid RANS-LES methods, where part of the domain may be solved using RANS, but the areas of high turbulence are solved using LES. Several such hybrid methods that have been developed are available in the commercial software CFX and are investigated in this paper.
The hybrid models are evaluated in their ability to simulate the experimental film cooling flow of a flat plate with a shaped hole at low Reynolds number, with a focus on their ability to predict the adiabatic effectiveness. The relative advantages and weaknesses of each model in predicting such film cooling flow are then apparent, and a preferred hybrid model for predicting low Reynolds number film cooling is identified. The added understanding of the benefits and shortcomings of hybrid models is important in the design of industrial applications of film cooling, where the designer may prefer a unified model to reduce setup time when predicting the effectiveness of a particular cooling scheme. This work also bears importance to research and academia, as multiple models in development have seen limited evaluation in the application of film cooling. The rapid-

Reynolds-Averaged Navier-Stokes
The Navier-Stokes momentum equation [10] requires closure of the nonlinear convective term. Equation (1) of this term is commonly accomplished by the decomposition and averaging first suggested by Reynolds [11], leading to the Reynolds-Averaged Navier-Stokes momentum equation, shown in Equation (2).
Here, u is the velocity vector, x is the direction, t is time, p is pressure, ρ is density, µ is dynamic viscosity, and S is the strain rate tensor, and the overbar represents the mean-flow property.
The Boussinesq eddy viscosity approximation [12] can then be employed. For incompressible flows this is: where t µ is the turbulent eddy viscosity, k is the turbulent kinetic energy, and ij δ is the Kronecker delta.
In order to calculate t µ , the k-ω model developed by Wilcox [13] is used, where t t k µ ρ ν ω = = , and k is the turbulent kinetic energy and ω is the specific turbulence dissipation rate. The following equations are used to calculate k and ω: and where Ω is the mean rotation tensor.
The shear stress transport (SST) model, developed by Menter [14], improves the k-ω model by applying a limiter to the eddy viscosity term to account for the transport of turbulent shear stress: where * 0. This SST model is used in the RANS regions of all the hybrid models as well as in the URANS model.

Large Eddy Simulation
Large eddy simulation, first proposed by Smagorinsky [15] and later developed by Deardoff [16], offers another approach to the Navier-Stokes equations by filtering the time-dependent equations. Turbulent eddies larger than the filter are resolved and eddies smaller than the filter are modelled using a subgrid-scale model. The filter allows any field variable to be decomposed into filtered and sub-filtered components, and leads to the filtered momentum equation: where ij τ defines the subgrid-scale stress: The sub-grid scale model then models the deviatoric part of ij τ , The Smagorinsky model does this with: where t µ is the turbulent eddy viscosity, modeled in the Smagorinsky model by: where S is defined by The wall-adapting local eddy-viscosity (WALE) model [17] improves upon the Smagorinsky model with the ability to correctly compute asymptotic behavior near the wall and laminar to turbulent transition. In the WALE model, the turbulent eddy viscosity is modeled by: where g is velocity gradient tensor, Additionally, w C is the Smagorinsky constant which is set to 0.5. This LES WALE model is used in the LES region of the Zonal LES method.

Models Evaluated
Four different hybrid models are evaluated: scale adaptive simulation, shielded detached eddy simulation, stress-blended eddy simulation, and zonal eddy simulation. One variation of the zonal eddy simulation is examined, along with the URANS result and a no-model laminar result. The RANS model used in the hybrid methods as well as the URANS model is the k-ω shear stress transport model [14].
A brief description of each model is given. Additional details of each model can be found in the referenced sources.

Scale Adaptive Simulation (SAS)
The scale-adaptive simulation model is an improved URANS formulation which dynamically adjusts from the URANS modeled turbulence to resolved structures using the introduction of a von Kármán length scale in the turbulence scale equation. This results in an LES-like behavior in unsteady regions of the flow and allowing resolution of the turbulent spectrum, while maintaining the RANS in the rest of the domain. The model was first described by Menter and Egorov [18], with additional detail provided by Egorov, Menter, Lechner, and Cokljat [19].
The SAS model modifies the RANS-SST turbulent eddy frequency transport equation, shown in Equation (5), by the addition of a source term, SAS Q , to the RHS of the ω transport equation. This source term is computed using the von Kármán constant, κ : where 2 3.51 where c µ is a replacement in the k transport equation for k β , 1 S is the scalar invariant of the strain rate tensor, and U ′′ is the second velocity derivative.

Shielded Detached Eddy Simulation (SDES)
Detached eddy simulation (DES) was the first proposed hybrid approach combining LES and RANS. The concept behind DES is to maintain RANS modeling near the wall where eddies are attached, but to switch to LES where the flow separates and the eddies become detached. The DES model replaces the distance term d in the RANS model with a distance function d  : where DES C is an empirical constant with a value of 0.65.
To incorporate DES into the k-ω SST model, the dissipation term in the turbulent kinetic energy transport equation, shown in Equation (4), is modified to include the function DES F : where DES F is defined by the turbulence length scale, t L , the maximum grid length in any direction, max ∆ , and DES C : In order to prevent this limiter from becoming active in the attached portion of the boundary layer which can happen when max ∆ is less than the boundary layer thickness, causing a grid-induced separation, the delayed DES model (DDES) [20] adds a function inside of DES F : + Ω (22) in which S is the magnitude of the strain rate tensor, Ω is the magnitude of the vorticity tensor. Improving upon DES, the improved delayed DES (IDDES) model [ This SDES model, described by Gritskevich, Garbaruk, Schütze, and Menter [22], represents the current state of the art for detached eddy simulation and is the version evaluated in this study.

Stress-Blended Eddy Simulation (SBES)
An even more recent development is the stress-blended eddy simulation (SBES), which is recommended to supersede SDES. It employs the same shielding function as SDES, but also adds the ability to blend between the RANS and LES regions. The modeled stress tensor of the RANS and LES regions is blended: This model has demonstrated faster development of turbulent scales over SDES and has predicted a better match of data in turbulent mixing problems.
The SBES model was developed and described by Menter [23].

Zonal Large Eddy Simulation (ZLES)
The last hybrid method studied is zonal LES (ZLES). In this model an LES zone is specified by the user. An interface is created at the boundary of the zone, where the modelled unresolved turbulence of the RANS model is converted to resolved fluctuations for the LES zone. To do this, a source term is added to the momentum equations and a sink term is added to the modelled turbulent kinetic energy equation: , , where Δt is the timestep size and where: where the length scale is determined by 0.
This interfacing method for zonal LES is described by Menter, Garbaruk, and Smirinov [24]. This is used to combine the WALE LES model k-ω SST RANS model for the ZLES model.

Zero Smagorinsky Coefficient ZLES
To explore the impact of the Smagorinsky constant on the LES solution, the ZLES model is run again with this constant set to zero instead of using the WALE value. Setting the Smagorinsky constant to zero fixes the eddy viscosity to zero, effectively eliminating the subgrid-scale model. This would be analogous to DNS; however, the RANS zones are still specified, the LES advection scheme is still being used, and near-wall treatment is still applied. This model is referred to herein as Smag0 ZLES.

Laminar
In addition to the hybrid LES methods, an unsteady no-model case is explored.
No turbulence model or sub-grid scale model are used; only the laminar equations are employed. This can be considered an implicitly filtered LES as the grid acts as a low pass filter. This can also be referred to as a coarse or under-resolved DNS.

Model Comparison
All the models studied are executed on the same mesh, with the same boundary conditions, initial conditions, timestep size, and time duration. The differences in the results are wholly due to the differences in the models, allowing examination of the similarities and differences amongst them and their relative utility for film cooling CFD. A brief description of each model is given in Table 1.

Experimental Setup
The CFD model is based on the flat plate film cooling experiments conducted by Schroeder and Thole [25] [26] [27] [28]. A closed-loop wind tunnel is used to flow air that is maintained at 295 K through a test section as shown in Figure 1.
A portion of the air is bled off to be cooled and then is injected back into the mainstream through a single row of five film cooling holes. The cooling holes are on a flat plate made out of expanded polystyrene for its low thermal conductivity. The mass flow rate of the coolant is controlled to provide a specific density ratio and blowing ratio. A range of density ratios and blowing ratios are experimented, but this CFD study explores the experimental case with DR = 1.5 and BR = 1.5. In order to spatially measure the surface temperature on the cooled plate, infrared imaging is taken through a window in the test section. From this, the adiabatic film cooling effectiveness over the surface is calculated.
The cooling holes are a baseline fan-shape geometry that is publicly defined as the "7-7-7" shaped hole [27] [28], referring to the 7˚ layback angle and symmetric 7˚ lateral angles, though also referred to as the "30-7-7" shaped hole to include Open Journal of Fluid Dynamics the 30˚ injection angle and assume symmetry. The geometry of the hole is shown in Figure 2 and the details of the geometric parameters are given in Table 2. Open Journal of Fluid Dynamics  This hole shape provides challenges in resolving the inherit unsteadiness of the flows, making it well-suited to for this study.

Computational Domain
The mesh used in this study is a structured grid generated in the commercial CFD mesher ANSYS ICEM. The hole is meshed using an O-grid structure, and the test section and plenum are block-structured. The mesh was initially generated for half the domain then mirrored to guarantee symmetry and rule out any possible mesh-induced flow asymmetry. The first layer thickness of the cooled wall obtained a nominal y + of 3, and a nominal x + and z + of 12. A y + less than one was desirable but traded off for maintaining a low aspect ratio and reducing x + and z + to ensure minimal filter shape biasing of the flow within the boundary layer. The higher y + is mitigated by this being only an adiabatic mixing problem with scalable wall functions that are intrinsic to the k-ω SST model in CFX. The  where L Kol = (ν 3 /ε) 1/4 and t Kol = (ν/ε) 1/2 [29], and spatially-accurate values for ε and ν are from the time-averaged URANS solution output. By the metrics shown in Table 3, the mesh used is sufficient for wall-resolved LES. The only exception taken related to the mesh is of CFL, defined as u t x ∆ ∆ . The CFL requirement of being close to 1 is not obtained in the simulation, with a nominal CFL value of 3 in the test section which drops below 2 within the boundary layer and grows to 5 inside the hole. A CFL greater than 1 may dampen turbulence, but is accepted as CFX is a fully implicit code and, as shown later, the turbulent eddy frequencies

Mesh Validation
Having established that the mesh meets notional wall-resolved LES quality requirements, a grid sensitivity study is subsequently performed by applying varying refinement levels to the grid. The levels are incremented in quarter increments of the base length of the original mesh, from 1/4 to 5/4, and the timestep size is held constant. Although grid convergence for LES is the DNS solution, which is outside the scope and computational resources for this work, stability in the simulated flow field can be observed as the mesh is refined, suggesting grid insensitivity is achieved by the original mesh. In Figure 4  To observe the convergence of the adiabatic effectiveness on the cooled surface, Figure 5 shows the area-average effectiveness for the entire surface versus the grid refinement level. The number of elements for each grid is shown next to each point, along with the percent change in the area-average surface temperature from the preceding refinement level. As the mesh is refined from the 3/4 to the 1/1 grid level, the percent change is +6.13%. As it is refined from the 1/1 level to the 5/4 level, the area-average surface temperature changes by −1.72%. This verifies the level of grid insensitivity that adiabatic effectiveness has for the 1/1 level mesh used.

Numerical Methodology
Ideal gas properties for air are used with Sutherland's formula for dynamic viscosity and thermal conductivity is applied with a reference Prandtl number for of 0.707 at 293.15 K and 1 atm. The numerical equations for total energy including viscous work, momentum and mass, wall scale, and the applicable turbulence model were iterated in the solver. For all of hybrid models, the 2 nd order high-resolution advection scheme, developed by Barth and Jespersen [31], was used in the k-ω SST RANS regions, which switched in the LES region to a central differencing scheme, based on the normalized variable diagram approach of Leonard [32], and the convection boundedness criterion developed by Jasak, Weller, and Gosman [33]. The transient scheme is the high-resolution scheme which uses second order backward Euler scheme. The turbulence numerics are also the second order high resolution scheme.

Initialization
To initialize each of the hybrid simulations, a steady state RANS solution was first attained, which was then used to initialize a large timestep size (1e−3 s) URANS simulation which ran for one complete flow-through. That was then

Boundary Conditions
The inlet was placed 12D upstream of the hole and set to a velocity inlet at a temperature of 295 K. The experimental 2D velocity and turbulence intensity profiles measured 2.3D upstream of the hole were applied to the inlet without adjustment. The measured turbulence intensity was 0.5% across most of the test section and increases to 13% near the wall. The turbulence length scale was set to 22% of the inlet boundary layer thickness. The cooling plenum inlet is located 5D from the hole entrance. The coolant flow is set to a temperature of 196.7 K to accomplish a density ratio of 1.5 and the mass flow rate is set to achieve a blowing ratio of 1.5. The turbulence intensity at the plenum inlet is set to 5%, and turbulent viscosity ratio is set to 10. The outlet was located 40D downstream of the hole and is set to a static pressure opening at 1 atm, with backflow enabled.
The test section height was set to 5D, with a symmetry condition applied to the far wall, understanding the value of reducing the domain size versus the possible source for inaccuracy since the entire test section height is not modelled. The lateral sides of the domain were set to translational periodic interfaces. The remaining boundaries were set to no-slip adiabatic walls. A diagram of the computational domain is provided in Figure 6, and the computational setup is summarized in Table 4.

Results
The time-averaged adiabatic effectiveness contours on the cooled surface are presented in Figure 7 for each of the models along with the experimental data.
By inspection, a few differences between models can be noted. The    Additional contours of adiabatic effectiveness are presented at flow-normal slice planes located at distances from the hole trailing edge of X/D = 0, 5, 10, 15, 20, 25, 30, and 35 in Figure 8, showing how the coolant film is spreading laterally and growing as it progresses down the test section. The asymmetry of the SBES, SDES, and SAS models is evident in these contours as well.
To quantify the accuracy of each model in predicting film cooling effectiveness, the adiabatic effectiveness along the hole centerline and the lateral-average adiabatic effectiveness are plotted in Figure 9. The maximum and RMS of the percent deviation from the experimental data of the centerline and lateral-average effectiveness are shown in Figure 10. From these, it is shown that the ZLES and SDES models most accurately capture the centerline effectiveness, with a maximum absolute deviation of 8% occurring at X/D = 1 for the SDES, and 11% occurring at X/D = 23 for the ZLES. The RMS of the centerline deviation is also lowest for the SDES and ZLES, at 5% and 8%, respectively.
More important to the overall effectiveness of the cooling scheme, the lateral-average effectiveness on the surface is best captured by the ZLES and laminar Open Journal of Fluid Dynamics   and an RMS ≤ 2.5% for the centerline effectiveness, and a maximum ≤ 3.5% and an RMS ≤ 2.6% for the lateral-average effectiveness. The low deviation values suggest that the simulations were run for a sufficient duration for a well-converged time-averaged simulation.
In Figure 12, the Q-criterion is shown at the last timestep for each model at a The turbulence spectra for each model are presented in Figure 13. The spectra are of the sampled total turbulent kinetic energy, the sum of the resolved and unresolved, over the sampled run history. The resolved turbulent kinetic energy is calculated from velocity statistics by where the overbar denotes the time-average value. The unresolved turbulent kinetic energy is calculated from the subgrid-scale model by   fast Fourier transform of the time signal with a Hann filter. The spectra may be difficult to discern from one another, but it should be noted that the turbulence simulated by ZLES and Smag0 ZLES model spectra are well above those of the other models. It can also be seen that the −5/3 slope of the inertial range is only attained for a short range, roughly between 500 and 2000 Hz. The relatively short inertial subrange is consistent with low Reynolds number spectra [29]. At frequencies above the inertial subrange, dissipation is dominant and the energy quickly drops off.
In Figure 14, the ZLES spectrum is plotted along with the calibrated model spectra function proposed by Pope [29]. The wavenumber on the bottom axis is Integrating under the model spectra from wavenumber 6 representing the energy containing range to wavenumber 445 where the mesh resolution filters out higher frequencies and dividing that by the area under the entire curve yields that 99.77% of the κE area is contained by the wavenumbers below the filter cutoff. This suggests the timestep and grid size allow a sufficient level of turbulence to be resolved. Also shown in Figure 14 spectra from DNS simulations conducted by Kim, Moin, and Moser [36]. The DNS spectra are for an Re τ of 180 and 395, and at y + = 30. This experiment has an Re τ of 315, and the spectra are sampled at y + = 25, so the DNS spectra serve as a reference only. Albeit, the LES spectra falls between the lower and higher shear velocity Reynolds number spectra as expected. It also appears to achieve a lower maximum energy level, which is attributable to the lower y + location.
In order to quantify the energy contained by the spectra of each model, Figure   15 shows the integrated area under the spectra for each model, κE. Perhaps the most telling figure, this shows that the ZLES and Smag0 ZLES models calculate a Open Journal of Fluid Dynamics  significantly higher overall turbulence spectrum at the hole leading edge. Furthermore, even though the laminar model carries no subgrid turbulence, it still is capturing a significantly larger κE than the SBES, SDES, and SAS models.
The lower turbulence calculated by the SBES and SDES models is attributed to the boundary layer shielding functions forcing a RANS zone upstream of the hole near the wall, which does not allow turbulence to develop upstream of the hole. This behavior of attached eddies until geometry-induced separation occurs is consistent with detached eddy simulation theory. Unfortunately, in this case it means that the near-wall turbulence upstream is not being sufficiently developed by either SBES or SDES. Additionally, despite the SAS model not having such shielding functions and having a larger scale-resolving zone upstream of the hole, the turbulence is not developing nearly as much as in the ZLES model or even the SDES and SBES.

Conclusion
A computational study of hybrid RANS-LES models was performed to examine the ability of modern hybrid models in predicting film cooling effectiveness of a low Reynolds number flat plate film cooling experiment. It was found that the Open Journal of Fluid Dynamics shielding of the RANS boundary layer by the SBES and SDES models precluded the development of upstream turbulence that impacted the film cooling effectiveness. The SBES, SDES, and SAS models did not resolve a sufficient amount of turbulence and were rendered ineffective for low Reynolds film cooling simulation. Zonal LES and even the no-model laminar simulation exhibited higher upstream turbulence and were more accurate in capturing the film cooling effectiveness. Thus, despite the convenience of generic RANS-LES hybrid models, based on the results, they are not herein recommended for low Reynolds number film cooling. Further research is needed to examine their usefulness for higher Reynolds flow, closer to actual engine conditions, although the shielding functions will still act to shield the RANS boundary layer upstream of the cooling hole jet. Should a hybrid RANS-LES approach be desired, zonal LES is the recommended choice for low Reynolds film cooling applications by these conclusions.