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

Michael Boehler, Akshay Sudesh, Mark Turner

Department of Aerospace Engineering & Engineering Mechanics, University of Cincinnati, Cincinnati, USA.

DOI: 10.4236/ojfd.2020.104019                

Abstract

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.

Keywords

Film Cooling, Large Eddy Simulation, Detached-Eddy Simulation, Stress-Blended Eddy Simulation, Laidback Fan-Shaped Cooling Hole, Film Cooling Effectiveness

1. 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 stationary 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 rapidly-growing development of turbulence models offers many options, and this serves to catalogue the most recent of those models. Novelty is in the comparison of each model for film cooling and in the overall accuracy attained for the particular shaped hole film cooling case. This will illuminate existing literature by offering the results of each model for comparison and by describing the methodology used to obtain the results.

2. Theoretical Methodology

2.1. Reynolds-Averaged Navier-Stokes

The Navier-Stokes momentum equation [10] requires closure of the nonlinear convective term. Equation (1) shows the incompressible Navier-Stokes momentum equation in conservation form, where this nonlinear term is $\frac{\partial }{\partial x_j}\left(u_j{u}_i\right)$. Closure 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).

$\rho \frac{\partial u_i}{\partial t}+\rho \frac{\partial }{\partial x_j}\left(u_j{u}_i\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left(2\mu S_{ij}\right)$ (1)

$\rho \frac{\partial {\bar{u}}_i}{\partial t}+\rho \frac{\partial {\bar{u}}_i{\bar{u}}_j}{\partial x_j}=-\frac{\partial \bar{p}}{\partial x_i}+\frac{\partial }{\partial x_j}\left(2\mu {\bar{S}}_{ij}-\rho \overline{{u^{\prime }}_i{u^{\prime }}_j}\right)$ (2)

Here, u is the velocity vector, x is the direction, t is time, p is pressure, $\rho$ is density, $\mu$ 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:

$-\rho \overline{{u^{\prime }}_i{u^{\prime }}_j}={\mu }_t\left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}\right)-\frac{2}{3}\rho k{\delta }_{ij}$ (3)

where ${\mu }_t$ is the turbulent eddy viscosity, k is the turbulent kinetic energy, and ${\delta }_{ij}$ is the Kronecker delta.

In order to calculate ${\mu }_t$, the k-ω model developed by Wilcox [13] is used, where ${\mu }_t/\rho ={\nu }_t=k/\omega$, and k is the turbulent kinetic energy and ω is the specific turbulence dissipation rate. The following equations are used to calculate k and ω:

$\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho u_{j}k\right)}{\partial x_j}=P_k-{\beta }_k\rho \omega k+\frac{\partial }{\partial x_j}\left[\left(\mu +{\sigma }_k\frac{\rho k}{\omega }\right)\frac{\partial k}{\partial x_j}\right]$ (4)

$\begin{array}{l}\frac{\partial \left(\rho \omega \right)}{\partial t}+\frac{\partial \left(\rho u_j\omega \right)}{\partial x_j}\\ =\frac{\alpha \omega }{k}{P}_k-{\beta }_{\omega }\rho {\omega }^2+\frac{\partial }{\partial x_j}\left[\left(\mu +{\sigma }_{\omega }\frac{\rho k}{\omega }\right)\frac{\partial \omega }{\partial x_j}\right]+\frac{p{\sigma }_d}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}\end{array}$ (5)

where $P_k={\tau }_{ij}\frac{\partial u_i}{\partial x_j}$, and ${\tau }_{ij}={\mu }_t\left(2S_{ij}-\frac{2}{3}\frac{\partial u_k}{\partial x_k}{\delta }_{ij}\right)-\frac{2}{3}\rho k{\delta }_{ij}$.

Additionally, ${\beta }_k=9/100$, ${\sigma }_k=3/5$, $\alpha =13/25$, ${\sigma }_{\omega }=1/2$, and ${\sigma }_d=1/8$ if $\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}>0$, else ${\sigma }_d=0$. Also, ${\beta }_{\omega }=0.0708\frac{1+85{\chi }_{\omega }}{1-100{\chi }_{\omega }}$, where ${\chi }_{\omega }=\frac{{\Omega }_{ij}{\Omega }_{jk}{S}_{ki}}{{\left({\beta }_k\omega \right)}^3}$, and where $\Omega$ 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:

${\mu }_t=\frac{\rho k}{\omega }\frac{1}{\max \left[\frac{1}{{\alpha }^{*}},\frac{S{F}_2}{0.31\omega }\right]}$ (6)

where ${\alpha }^{*}=\frac{0.024+R{e}_t/6}{1+R{e}_t/6}$, $F_2=\tanh \left({\varphi }_2^2\right)$, and ${\varphi }_2=\max \left[2\frac{\sqrt{k}}{0.09\omega d},\frac{500\mu }{\rho d^2\omega }\right]$. Additionally, $R{e}_t=\frac{\rho k}{\mu \omega }$, and d is the distance to the nearest surface.

This SST model is used in the RANS regions of all the hybrid models as well as in the URANS model.

2.2. 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:

$\frac{\partial }{\partial t}\left(\rho {\bar{u}}_i\right)+\frac{\partial }{\partial x_j}\left(\rho {\bar{u}}_i{\bar{u}}_j\right)=-\frac{\partial \bar{p}}{\partial x_i}+\frac{\partial }{\partial x_j}\left(\mu \frac{\partial {\bar{u}}_i}{\partial x_j}\right)-\frac{\partial {\tau }_{ij}}{\partial x_j}$ (7)

where ${\tau }_{ij}$ defines the subgrid-scale stress:

${\tau }_{ij}\equiv \rho \overline{u_i{u}_j}-\rho {\bar{u}}_i{\bar{u}}_j$ (8)

The sub-grid scale model then models the deviatoric part of ${\tau }_{ij}$, ${\tau }_{ij}-\frac{1}{3}{\tau }_{kk}{\delta }_{ij}$. The Smagorinsky model does this with:

${\tau }_{ij}-\frac{1}{3}{\tau }_{kk}{\delta }_{ij}=2{\mu }_t{S}_{ij}$ (9)

where ${\mu }_t$ is the turbulent eddy viscosity, modeled in the Smagorinsky model by:

${\mu }_t=\rho {\left(C_S\Delta \right)}^2\left|\bar{S}\right|$ (10)

where $\left|\bar{S}\right|$ is defined by $\left|\bar{S}\right|\equiv \sqrt{2{\bar{S}}_{ij}{\bar{S}}_{ij}}$, $C_S$ is the Smagorinsky constant which typically ranges from 0.1 to 0.2, and $\Delta$ is the local grid scale, computed as the cube root of the element volume.

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:

${\mu }_t=\rho {\left(C_w\Delta \right)}^2\frac{{\left(S_{ij}^d{S}_{ij}^d\right)}^{3/2}}{{\left({\bar{S}}_{ij}{\bar{S}}_{ij}\right)}^{5/2}+{\left(S_{ij}^d{S}_{ij}^d\right)}^{5/4}}$ (11)

where $S_{ij}^d$ is defined by:

$S_{ij}^d=\frac{1}{2}\left({\bar{g}}_{ij}^2+{\bar{g}}_{ji}^2\right)-\frac{1}{3}{\delta }_{ij}{\bar{g}}_{kk}^2$ (12)

where g is velocity gradient tensor, ${\bar{g}}_{ij}=\frac{\partial {\bar{u}}_i}{\partial x_j}$. Additionally, $C_w$ 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.

2.3. 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.

2.3.1. 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, $Q_{SAS}$, to the RHS of the ω transport equation. This source term is computed using the von Kármán constant, $\kappa$ :

$Q_{SAS}=\max \left[\rho {\zeta }_2\kappa S_1^2{\left(\frac{L}{L_{vK}}\right)}^2-C\frac{2\rho k}{{\sigma }_{\varphi }}\max \left(\frac{1}{{\omega }^2}\frac{\partial \omega }{\partial x_j}\frac{\partial \omega }{\partial x_j},\frac{1}{k^2}\frac{\partial k}{\partial x_j}\frac{\partial k}{\partial x_j}\right),0\right]$ (13)

where ${\zeta }_2=3.51$, $C=2$, ${\sigma }_{\varphi }=2/3$, L is the modeled turbulence length scale and $L_{vK}$ is the von Kármán length-scale:

$L=\sqrt{k}/\left(c_{\mu }^{1/4}\omega \right)$ (14)

$L_{vK}=\frac{\kappa s}{\left|U^{″}\right|}$ (15)

where $c_{\mu }$ is a replacement in the k transport equation for ${\beta }_k$, $S_1$ is the scalar invariant of the strain rate tensor, and $U^{″}$ is the second velocity derivative.

2.3.2. 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 $\tilde{d}$ :

$\tilde{d}=\min \left(d,C_{DES}\Delta \right)$ (16)

where $C_{DES}$ 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 $F_{DES}$ :

${\epsilon }_k=\rho {\beta }_{k}k\omega F_{DES}$ (17)

where $F_{DES}$ is defined by the turbulence length scale, $L_t$, the maximum grid length in any direction, ${\Delta }_{\max }$, and $C_{DES}$ :

$F_{DES}=\max \left(\frac{L_t}{C_{DES}{\Delta }_{\max }},1\right)$ (18)

$L_t=\frac{\sqrt{k}}{{\beta }_k\omega }$ (19)

In order to prevent this limiter from becoming active in the attached portion of the boundary layer which can happen when ${\Delta }_{\max }$ is less than the boundary layer thickness, causing a grid-induced separation, the delayed DES model (DDES) [20] adds a function inside of $F_{DES}$ :

$F_{DES}=\max \left[\frac{L_t}{C_{DES}{\Delta }_{\max }}\left(1-f_{DDES}\right),1\right]$ (20)

$f_{DDES}=1-\tanh \left[{\left(C_{d1}{r}_d\right)}^{C_{d2}}\right]$ (21)

where $C_{d1}=20$, $C_{d2}=3$, and $r_d$ is computed by:

$r_d=\frac{{\mu }_t+\mu }{\rho {\kappa }^2{d}^2\sqrt{0.5\left(S^2+{\Omega }^2\right)}}$ (22)

in which S is the magnitude of the strain rate tensor, $\Omega$ is the magnitude of the vorticity tensor.

Improving upon DES, the improved delayed DES (IDDES) model [21] and shielded DES (SDES) model improve upon DES through modifications to the $f_{DDES}$ function that help shield the RANS wall boundary layer from the LES region and offering a faster transition to LES by applying a formula to the grid length scale:

${\Delta }_{SDES}=\max \left(\sqrt[3]{\text{Volume}},0.2{\Delta }_{\max }\right)$ (23)

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.

2.3.3. 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:

${\tau }_{ij}{}_{{}_{SBES}}=f_{SDES}{\tau }_{ij}{}_{{}_{RANS}}+\left(1-f_{SDES}\right){\tau }_{ij}{}_{{}_{LES}}$ (24)

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].

2.3.4. 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:

$F_{mom,i}=\frac{\rho u_{f,i}}{\text{Δ}t}$ (25)

$F_k=-0.5\frac{\rho u_{f,i}^2}{\text{Δ}t}$ (26)

where $\text{Δ}t$ is the timestep size and $u_{f,i}$ is the velocity fluctuation to be transferred. The velocity fluctuation is determined by a random flow generator that uses the RANS turbulent length and time scales for input:

$u_{f,i}=\sqrt{\frac{2}{3}k}\sqrt{\frac{2}{N}}\underset{n=1}{\overset{N}{{\displaystyle \sum }}}\left[p_i^n\cos \left(ar{g}^n\right)+q_i^n\sin \left(ar{g}^n\right)\right]$ (27)

where:

$p_i^n={\epsilon }_{ijk}{\eta }_j^n{d}_k^n$, $q_i^n={\epsilon }_{ijk}{\xi }_j^n{d}_k^n$ (28)

$ar{g}^n=2\pi \left(\frac{d_i^n{x}_i}{L_t}+\frac{{\omega }^{n}t}{{\tau }_t}\right)$ (29)

where the length scale is determined by $L_t=0.5\sqrt{k}/c_{\mu }\omega$, and the time scale is ${\tau }_t=L_t/\sqrt{k}$. Additionally, random distributions $N\left(\varphi ,\psi \right)$, where $\varphi$ is the mean and $\psi$ is the standard deviation, are the following:

${\eta }_i^n=N\left(0,1\right)$, ${\xi }_i^n=N\left(0,1\right)$, $d_i^n=N\left(0,0.5\right)$, ${\omega }_i^n=N\left(1,1\right)$ (30)

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.

2.3.5. 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.

2.3.6. 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.

2.4. 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.

3. 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

Table 1. Turbulence models.

Figure 1. Closed loop cooled flat plate experimental facility [25] [26] [27] [28].

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.

Table 2. Shaped hole geometry parameters.

Figure 2. Fan-shaped hole geometry.

This hole shape provides challenges in resolving the inherit unsteadiness of the flows, making it well-suited to for this study.

4. Computational Setup

4.1. 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 final mesh size was 18 million elements. The side view is shown with a closeup near the hole in Figure 3(a), and a view of the mesh on the cooled surface is shown in Figure 3(b).

Table 3 highlights the boundary layer grid requirements for LES in CFX, and where this mesh stands in terms of those requirements. Nx, Ny, and Nz are the number of nodes within the boundary layer in the x, y, and z directions, respectively. The Kolmorgorov scales required for DNS are also shown for reference, where L_Kol = (ν³/ε)^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 $\frac{u\Delta t}{\Delta 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

Table 3. LES grid and timestep requirements.

Figure 3. Multi-block structured symmetric mesh used in all simulations.

are resolved well beyond the inertial subrange for the relatively low Reynolds number flow.

4.2. 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, the velocity profile at a location 10D upstream of the leading edge of the hole and at the lateral center is shown for each of the grid levels. The y⁺ locations of the grid layers for each mesh level are shown by the markers on each line. The 1/4 level grid has 25 layers in the y direction, whereas the 5/4 level grid has 122 and the 1/1 level grid used has 97. The measured profile from the experiment is shown for reference with the sample locations shown with markers. This experimental profile is applied at the inlet of the computational domain. As the mesh is refined from the 1/4 level to the 5/4 level, the shear velocity, u_τ, decreases, which in turn increases u⁺ and slightly decreases y⁺. The calculated u_τ value for each mesh is shown in the legend. Also in the legend is the RMS of the percent change in the velocity u along the y direction from the preceding grid level. The curve for u⁺ = y⁺ is displayed, which all meshes follow closely up to y⁺ = 5. As can be seen, the velocity profile exhibits minimal change from the 3/4 to the 1/1 and 5/4 level grid, with

Figure 4. Grid convergence of velocity profile 10D upstream of hole.

an RMS of 1.76% from the 3/4 to the 1/1 level and an RMS of 0.68% from the 1/1 to 5/4 level grid. This demonstrates that the velocity upstream of the hole is grid insensitive for the 1/1 level grid used.

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.

4.3. 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.

4.4. 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

Figure 5. Convergence of area-average effectiveness.

used to initialize a to a smaller timestep size (1e−4 s) URANS simulation which ran for three more flow-throughs. The length of the domain from inlet to outlet is 0.424 m, and the nominal flow velocity is 10 m/s, so 424 timesteps were required for one flow-through at this timestep size. The URANS results presented are of the smaller timestep simulation, with data collection starting after the first smaller timestep flow-through to the end of the fourth flow-through, with the time-average results of this range exhibiting good convergence (<5% change between flow-throughs in centerline and lateral adiabatic effectiveness along the cooled plate). The final state of the URANS solution was then used to initialize each of the hybrid simulations. The hybrid simulations used a timestep size of 1e−4 s as well and were run for one complete flow-through before data sampling, and then run for five more flow-throughs, collecting the time-averaged data presented herein. The time-averaged results were well converged as shown later, and the timestep size was sufficient in resolving the inertial subrange frequencies.

4.5. 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.

5. 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 URANS surface contours of adiabatic effectiveness are symmetric with two sharp streaks of higher effectiveness, whereas all of the scale-resolving simulation models exhibit a single streak spread across the lateral direction. In comparing the Zonal LES to

Figure 6. Computational domain and boundary conditions.

Figure 7. Time-averaged adiabatic effectiveness contours on flat plate surface.

Table 4. Computational setup summary.

the other hybrid models, the SBES, SDES, and SAS models exhibit an asymmetry where the streak favors one side or the other that is not apparent to the same extent in the ZLES model or in the test data. Albeit, the ZLES model and test data are not perfectly symmetric, and symmetry is not necessary for this complex flow field. Furthermore, the experiments of Aghasi and Gutmark [34] demonstrated varying levels of asymmetry for differently manufactured test plate coupons. 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

Figure 8. Spanwise adiabatic effectiveness contours from X/D = 0 to X/D = 35.

Figure 9. Adiabatic effectiveness of (a) centerline and (b) lateral average.

Figure 10. Effectiveness deviation from experimental data of (a) centerline and (b) lateral average.

simulation. The largest absolute deviations from the test data for these are immediately downstream of the trailing edge of the hole, with a value of 8% for the ZLES and a value of 12% for the laminar simulation. The RMS of the lateral-average deviation is also lowest for the ZLES and laminar simulation, with values of 5% and 4%, respectively.

In order to verify that a sufficient number of flow-throughs have been run for transient averaging, the convergence of the centerline and lateral-average effectiveness after each complete flow-through is evaluated. For reference, plots of the centerline and lateral-average effectiveness convergence after each flow through of the ZLES model are shown in Figure 11. The running time-averaged effectiveness is plotted at the end of each flow-through. From the second-to-last to the last flow-through, all models have achieved a maximum deviation ≤ 4.1% 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 normalized level of 0.01. The instantaneous Q-criterion shows how the coherent turbulent structures develop in areas of high turbulence. The structures are most developed for the ZLES, Smag0 ZLES, and laminar solutions. The structures of the SDES and SBES solutions appear very similar, with no significant increase in coherent turbulence evident in the SBES model. The SAS model has substantially fewer developed structures of any of the models, as expected.

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 $1/2\left[{\left(u-\bar{u}\right)}^2+{\left(v-\bar{v}\right)}^2+{\left(w-\bar{w}\right)}^2\right]$, where the overbar denotes the time-average value. The unresolved turbulent kinetic energy is calculated from the subgrid-scale model by $1/0.3{\left(C_S\Delta S\right)}^2$. The spectra presented are sampled at a single point above the hole leading edge, on the centerline, and at a distance off the wall corresponding to the location where Re_θ becomes 670, to ensure relatively isotropic flow [35], which corresponds to a y⁺ of approximately 25 for this experiment. The time-accurate turbulent kinetic energy of the sampled time range is converted into the frequency domain using a

Figure 11. ZLES time-average effectiveness convergence for (a) centerline and (b) lateral average.

Figure 12. Q-Criterion isosurfaces, level 0.01, at the final timestep, colored by streamwise velocity.

Figure 13. Energy spectra of all models above hole leading edge, y⁺ = 25.

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 obtained by dividing frequency by the time-average velocity at that point. From the comparison of the ZLES spectrum and the model spectrum, it can be determined that the LES simulates the spectrum well at this location up to the grid cutoff frequency, which is approximately 4000 Hz or a wavenumber of approximately 450 m⁻¹. Beyond this, the higher wavenumber turbulence is modeled. 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

Figure 14. Model energy spectra of ZLES at y⁺ = 25 compared to available DNS at y⁺ = 30 [36].

Figure 15. Total spectral energy for each model, m²/s³.

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.

6. 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 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.

Acknowledgements

This work could not have been accomplished without the support of Dr. Karen Thole and Turbine Heat Transfer and Aerodynamics Group at Penn State, as well as the resources made available by the Ohio Supercomputer Center, and funding provided by the Ohio Federal Research Network.

Nomenclature

BR: blowing ratio, ρ_cu_c/ρ_∞u_∞

D: hole diameter

DR: density ratio

E(κ): energy spectrum function

k: turbulent kinetic energy

Re_θ: momentum thickness Reynolds number, θu_∞/ν

Re_τ: shear velocity Reynolds number, δu_τ/ν

T_c: coolant temperature

u_c: coolant jet velocity

u_τ: shear velocity

x⁺: streamwise dimensionless wall distance, u_τx/ν

y⁺: wall-normal dimensionless wall distance, u_τy/ν

z⁺: spanwise dimensionless wall distance, u_τz/ν

ε: turbulence dissipation rate

η: adiabatic effectiveness, $\frac{T_{\infty }-T}{T_{\infty }-T_c}$

κ: wavenumber

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Goldstein, R.J. (1971) Film Cooling. Advances in Heat Transfer, 7, 321-379. https://doi.org/10.1016/S0065-2717(08)70020-0
[2]	Bogard, D.G. and Thole, K.A. (2006) Gas Turbine Film Cooling. Journal of Propulsion and Power, 22, 249-270. https://doi.org/10.2514/1.18034
[3]	Bunker, R.S. (2005) A Review of Shaped Hole Turbine Film Cooling Technology. Journal of Heat Transfer, 127, 441-453. https://doi.org/10.1115/1.1860562
[4]	Walters, D.K. and Leylek, J.H. (1997) A Systematic Computational Methodology Applied to a Three-Dimensional Film-Cooling Flowfield. Journal of Turbomachinery, 119, 777-785. https://doi.org/10.1115/1.2841188
[5]	Hoda, A. and Acharya, S. (2000) Predictions of a Film Coolant Jet in Crossflow with Different Turbulence Models. Journal of Turbomachinery, 122, 558-569. https://doi.org/10.1115/1.1302322
[6]	Azzi, A. and Lakehal, D. (2002) Perspectives in Modeling Film Cooling of Turbine Blades by Transcending Conventional Two-Equation Turbulence Models. Journal of Turbomachinery, 124, 472-484. https://doi.org/10.1115/1.1485294
[7]	Wang, X., et al. (2017) Numerical Study on the Near-Wall Characteristics of Compound Angled Film Cooling Based on Hybrid Thermal Lattice Boltzmann Method. Applied Thermal Engineering, 129, 1670-1681. https://doi.org/10.1016/j.applthermaleng.2017.10.125
[8]	Simon, M., Gautier, S., Vanoli, E. and Auzillon, P. (2019) Aerothermal Simulations Comparisons of a Shaped-Hole Film Cooling Flow. Proceedings ASME Turbo Expo, Phoenix, 17-21 June 2019, V05AT12A006. https://doi.org/10.1115/GT2019-90458
[9]	Meneveau, C. (2010) Turbulence: Subgrid-Scale Modeling. Scholarpedia, 5, 9489. https://doi.org/10.4249/scholarpedia.9489
[10]	Stokes, G.G. (1842) On the Steady Motion of Incompressible Fluids. Transactions of the Cambridge Philosophical Society, 7, 439-453.
[11]	Reynolds, O. (1895) On the Dynamical Theory of Incompressible Viscous Fluids and the Determination of the Criterion. Philosophical Transactions of the Royal Society of London A, 186, 123-164. https://doi.org/10.1098/rsta.1895.0004
[12]	Boussinesq, J. (1877) Essai sur la Théorie des Eaux Courantes. Mémoires Présentés par Divers Savants à l’Académie des Sciences, 1-680.
[13]	Wilcox, D.C. (1988) Reassessment of the Scale-Determining Equation for Advanced Turbulence Models. AIAA Journal, 26, 1299-1310. https://doi.org/10.2514/3.10041
[14]	Menter, F.R. (1994) Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications. AIAA Journal, 32, 1598-1605. https://doi.org/10.2514/3.12149
[15]	Smagorinsky, J. (1963) General Circulation Experiments with the Primitive Equations. Monthly Weather Review, 91, 99-164. https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2
[16]	Deardorff, J. (1970) A Numerical Study of Three-Dimensional Turbulent Channel Flow at Large Reynolds Numbers. Journal of Fluid Mechanics, 41, 453-480. https://doi.org/10.1017/S0022112070000691
[17]	Nicoud, F. and Ducros, F. (1999) Subgrid-Scale Stress Modelling Based on the Square of the Velocity Gradient. Turbulence and Combustion, 62, 183-200. https://doi.org/10.1023/A:1009995426001
[18]	Menter, F.R. and Egorov, Y. (2010) The Scale-Adaptive Simulation Method for Unsteady Turbulent Flow Predictions. Part 1: Theory and Model Description. Flow, Turbulence and Combustion, 85, 113-138. https://doi.org/10.1007/s10494-010-9264-5
[19]	Egorov, Y., Menter, F.R., Lechner, R. and Cokljat, D. (2010) The Scale-Adaptive Simulation Method for Unsteady Turbulent Flow Predictions. Part 2: Application to Complex Flows. Flow, Turbulence and Combustion, 85, 113-138. https://doi.org/10.1007/s10494-010-9265-4
[20]	Spalart, P.R., Deck, S., Shur, M.L., Squires, K.D., Strelets, M.K. and Travin, A. (2006) A New Version of Detached-Eddy Simulation, Resistant to Ambiguous Grid Densities. Theoretical and Computational Fluid Dynamics, 20, 181-195. https://doi.org/10.1007/s00162-006-0015-0
[21]	Shur, M.L., Spalart, P.R., Strelets, M.K. and Travin, A.K. (2008) A Hybrid RANS-LES Approach with Delayed-DES and Wall-Modelled LES Capabilities. International Journal of Heat and Fluid Flow, 29, 1638-1649. https://doi.org/10.1016/j.ijheatfluidflow.2008.07.001
[22]	Gritskevich, M.S., Garbaruk, A.V., Schütze, J. and Menter, F.R. (2012) Development of DDES and IDDES Formulations for the k-ω Shear Stress Transport Model. Flow, Turbulence and Combustion, 88, 431-449. https://doi.org/10.1007/s10494-011-9378-4
[23]	Menter, F.R. (2018) Stress-Blended Eddy Simulation (SBES)—A New Paradigm in Hybrid RANS-LES Modeling. Notes on Numerical Fluid Mechanics and Multidisciplinary Design, 137, 27-37. https://doi.org/10.1007/978-3-319-70031-1_3
[24]	Menter, F.R., Garbaruk, A. and Smirinov, P. (2010) Scale Adaptive Simulation with Artificial Forcing. Notes on Numerical Fluid Mechanics and Multidisciplinary Design, 111, 235-246. https://doi.org/10.1007/978-3-642-14168-3_20
[25]	Schroeder, R.P. and Thole, K.A. (2016) Thermal Field Measurements for a Shaped Hole at Low and High Freestream Turbulence Intensity. Journal of Turbomachinery, 139, Article ID: 021012. https://doi.org/10.1115/GT2016-56967
[26]	Schroeder, R.P. and Thole, K.A. (2017) Effect of In-Hole Roughness on Film Cooling from a Shaped Hole. Journal of Turbomachinery, 139, Article ID: 031004. https://doi.org/10.1115/1.4034847
[27]	Schroeder, R.P. and Thole, K.A. (2017) Effect of High Freestream Turbulence on Flowfields of Shaped Film Cooling Holes. Journal of Turbomachinery, 138, Article ID: 091001. https://doi.org/10.1115/1.4032736
[28]	Schroeder, R.P. and Thole, K.A. (2014) Adiabatic Effectiveness Measurements for a Baseline Shaped Film Cooling Hole. Proceedings ASME Turbo Expo: Turbine Technical Conference and Exposition, Düsseldorf, 16-20 June 2014, V05BT13A036. https://doi.org/10.1115/GT2014-25992
[29]	Pope, S.B. (2000) Turbulent Flows. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511840531
[30]	ANSYS, Inc. (2015) ANSYS CFX Solver Theory Guide. Release 16.2.
[31]	Barth, T. and Jespersen, D. (1989) The Design and Application of Upwind Schemes on Unstructured Meshes. Proceedings 27th Aerospace Sciences Meeting, Reno, 9-12 January 1989, 366. https://doi.org/10.2514/6.1989-366
[32]	Leonard, B.P. (1991) The ULTIMATE Conservative Difference Scheme Applied to Unsteady One-Dimensional Advection. Computer Methods in Applied Mechanics and Engineering, 88, 17-74. https://doi.org/10.1016/0045-7825(91)90232-U
[33]	Jasak, H., Weller, H.G. and Gosman, A.D. (1999) High Resolution NVD Differencing Scheme for Arbitrarily Unstructured Meshes. International Journal for Numerical Methods in Fluids, 31, 413-419. https://doi.org/10.1002/(SICI)1097-0363(19990930)31:2<431::AID-FLD884>3.0.CO;2-T
[34]	Aghasi, P. and Gutmark, E. (2017) Dependence of Film Cooling Effectiveness on 3D Printed Cooling Holes. Journal of Heat Transfer, 139, Article ID: 102003. https://doi.org/10.1115/1.4036509
[35]	Stratton, Z.T. and Shih, T.I.-P. (2019) Identifying Weaknesses in Eddy-Viscosity Models for Predicting Film Cooling via Large-Eddy Simulations. Journal of Propulsion and Power, 35, Article ID: 011006. https://doi.org/10.2514/1.B37314
[36]	Kim, J., Moin, P. and Moser, R. (1987) Turbulence Statistics in Fully Developed Channel Flow at Low Reynolds Number. Journal of Fluid Mechanics, 177, 133-166. https://doi.org/10.1017/S0022112087000892