1. Introduction
A widely used strategy to increase efficiency and specific power output of gas turbines consists of increasing the Turbine Entry Temperature (TET). Therefore, metal temperatures approach the melting point of the alloys used in high-pressure nozzles, which are greatly loaded as shown in [
1]. Advanced cooling systems such as film cooling are necessary to protect nozzles from high temperatures. Such technology consists of the injection of coolant flow spilled from the high-pressure compressor stage to create a protecting fluid film on the blade surface. The interaction of the coolant with the main flow increases complexity in the flow’s structures. The insurgence of counter-rotating vortices (called “kidney-shaped vortices”) with a velocity component perpendicular to the surface is described in [
2]. These vortices are generated by the redistribution of the vorticity content in the boundary layer and play a key role in the performance of the cooling devices. Kidney vortices cause a lift-off effect moving the coolant away from the wall and reducing the capability of film cooling in covering and protecting the blade surface. These three-dimensional structures are numerically investigated over flat plates in [
3]. They noticed that the kidney vortices are greatly influenced by the boundary layer development inside the coolant channel. An additional whirling structure, the so-called “tornado effect”, is described in [
4], where the presence of a core-wise, cross-flow transport in hairpin vortices (created by coolant injection) in the laminar boundary layer is shown. This core-wise phenomenon continuously moves the fluid from the wall to the free stream. Similar structures are described in [
5] for a turbulent boundary layer and in [
6] in presence of interaction between a shock and the boundary layer.
The complexity of the fluid structures is further increased when the film cooling device is considered in a realistic turbomachinery configuration. The effect of the cooling position on the losses of a turbine stage is analyzed in [
7]. They highlighted that the position and the inclination of a cooling jet strongly affects both the losses and the capability of the boundary layer to resist to an adverse pressure gradient. Moreover, in high-pressure stages oblique shocks shed from the Trailing Edge (TE) impinge onto the suction side of the adjacent vane. In case of presence of a row of film cooling holes their effectiveness is affected by the three-dimensional interaction between the cooling flow, the boundary layer development, and the shock itself. Several authors demonstrate a non-negligible variation in the local value of adiabatic effectiveness
when the hole is located a few diameters before the shock impingement zone [
8,
9,
10].
More recently Uncertainty Quantification (UQ) methodologies have been coupled with CFD simulations, providing the opportunity to quantify the stochastic variations associated with the inevitable manufacturing deviations and the fluctuation of boundary conditions in an engine-like environment [
11]. The sensitivity of film cooling performance to relevant non-dimensional parameters is discussed in [
12] without a statistical analysis. A UQ analysis is used in [
13] to demonstrate that small variations in the coolant duct geometry affect the coolant flow redistribution in a non-negligible way. These kinds of deviations are generally known as aleatory uncertainties. The different capability to propagate aleatory uncertainties using different CFD approaches such as Reynolds-Averaged Navier-Stokes (RANS), unsteady RANS or Large Eddy Simulation (LES), provides identification on the accuracy of different numerical approach as shown by [
14]. These deviations are known in open literature as epistemic uncertainty and are more related to the limitations of the numerical model. The quantification of the epistemic uncertainties is still an open problem.
In this work the effect of an impinging shock is studied numerically using a test case designed by the Institut für Thermische Ströemungsmaschinen in Karlsruhe [
9,
15]. A shaped plate positioned at the mid-height of a converging nozzle creates a diverging region where the Mach number reaches
, while at the end of the plate the flow slows down through a shock system. Two geometrical parameters are perturbed to introduce the aleatory uncertainty associated with manufacturing deviations, namely the diameter of both the plate trailing edge (D
) and the cooling hole (D
). The deterministic maps of adiabatic effectiveness
close to the hole exit section are obtained using two different turbulence models (the SST model [
16] and the Reynolds Stress Model [
17]) and compared with the available experimental data. The deterministic values and the stochastic distribution of the performance parameters (blowing ratio
, density ratio
, momentum ratio
and discharge coefficient
) are also calculated using two Probabilistic Collocation Methods (PCM) with Hermite polynomials and Padè–Legendre approximants. Present work underlines the limitations of turbulence modelling in the analysis of
maps and shows that the impact of manufacturing uncertainty is not relevant for the calculation of the performance parameters except for
, to which is associated a standard deviation up to
.
2. Test-Case Description
Numerical simulations are widely used in cooperation with experimental tests, but only recently CFD simulations have been coupled with UQ to study a wide range of turbomachinery problems [
13,
18,
19]. Concerning the present study, a complete description of the experimental facility is reported in [
9] and a detailed numerical campaign using the k-
turbulence model [
20] has been reported in [
10]. The proposed simulations have been carried out considering deterministic boundary conditions
Figure 1 shows a sketch of the control volume that includes the converging nozzle, the central (shaped) plate, the cylindrical cooling hole, and the plenum. The main flow reaches sonic speed at the throat while the shape of the lower part of the contoured plate allows the flow to accelerate further to supersonic velocities. At the trailing edge, an oblique shock wave is generated, which impinges on the lower wall in a region immediately downstream of the cooling-hole exit position. The flow physics (including the effectiveness distribution obtained in [
21]) is depicted in
Figure 2, which is obtained for geometry and boundary conditions corresponding to the experiments. The coolant interacts with the main flow generating a weak shock, which merges with the oblique shock shed by the plate trailing edge that impinges on the lower wall. The
distribution is modified by the presence of the adverse pressure gradient that is responsible for a local growth of the boundary layer.
In the experimental apparatus, five injection holes are located upstream of the shock impingement position and have a pitch-to-diameter ratio of
[
9]. The control volume selected for the present numerical campaign includes only a portion of the entire domain, being the channel symmetric along the hole mid-plane and between two adjacent holes. Therefore, the computational domain is limited in the z-direction by two symmetry planes that are set normal to the x-z plane (see
Figure 1 for reference). This choice neglects the formation of coolant oscillations on planes parallel to the x-z plane and is coherent with the steady nature of the computations presented in the present work. The theoretical free stream Mach number at the shock location is around
and the blowing ratio value of
has been selected. The non-dimensional coolant total pressure is
with respect to main-flow inlet total pressure and the non-dimensional coolant total temperature is
with respect to main-flow inlet total temperature. For each investigated case, the numerical value has been calculated using Equation (
1) in accordance with the experimental definition:
In Equation (
1)
is the coolant density,
is the main-flow density,
is the coolant velocity and
is the main-flow velocity. The value of the adiabatic effectiveness is essentially a non-dimensional representation of the adiabatic wall temperature
and is calculated using Equation (
2):
where
is the coolant stagnation temperature (evaluated approximately at the hole exit section) and
is the recovery temperature of the main flow. In case of perfect coverage of the wall by the coolant flow
is equal to
and
is
; in case of absence of coolant
is equivalent to
and
is
Calculations with an adiabatic condition for the cooled plate are used to evaluate the adiabatic wall temperature distribution. The evaluation of the main-flow recovery temperature has been performed using Equation (
3):
The value of the recovery temperature is a function of the recovery factor
r that in case of turbulent flows can be estimated as
[
22].
3. Numerical Methodology
The computational domain assumes flow symmetry at the hole centerline and between two adjacent holes. That assumption is justified both by the fact that the boundary conditions are uniform and that unsteady effects in the lateral direction are expected to be negligible. As shown in
Figure 1, the control volume includes the main-flow inlet, a coolant supply plenum, and the cylindrical cooling channel, the end-walls that define the nozzle, the central shaped plate, and the outlet. The presence of the plenum in film cooling simulations is crucial to correctly evaluate the aerodynamic losses and the development of secondary flows [
23,
24,
25,
26] and then it is included in the present computational control volume. Steady RANS equations are solved using a pressure-based approach. Pressure-velocity equations are coupled, and the calculation is fully second-order accurate. The Shear Stress Transport (SST) model [
16] and Reynolds Stress Model (RSM) [
17] already implemented in the commercial ANSYS
® FLUENT
® code are used as turbulence closures to try to overcome the limitations of the standard k-
model [
20]. The main difference between the two selected turbulence closures consists of the calculation of the stress-tensor. RSM can evaluate the non-isotropic features of turbulence while SST model considers the turbulence as an isotropic quantity. The comparison of the two different approaches is aimed at assessing the role of the non-isotropic turbulence features due to the shock/coolant interaction in the evaluation of the distribution of adiabatic effectiveness. Inlet turbulence level is set at
while a turbulence length scale of 1 mm has been considered. Walls are set as adiabatic and viscous heating is considered to make the heat transfer evaluation accurate. Calculation is converged when all the residuals are constantly below
and the mass-flow error is below
.
The uncertainty quantification analysis has been performed by considering two geometrical deviations as aleatory parameter: (
i) the diameter of the cooling channel; (
) the diameter of trailing edge of the central plate. Both deviations are considered to obtain a uniform distribution. The deviations related to the geometrical parameters have been assumed as
of the nominal value for the coolant diameter and
for the trailing-edge diameter [
11]. Selected deviations are slightly higher than the realistic ones to obtain a wider response surface as result. The test matrix of the numerical campaign is summarized in
Table 1 using non-dimensional
D values with respect to the nominal configuration.
The meshes used for the nine calculations (necessary for the UQ analysis) are generated using the commercial hybrid grid generator Centaur
TM by Centaursoft. The mesh is denser in the region where the shock-boundary layer interaction is expected (see
Figure 3a) and the overall quality of the mesh is much higher than in [
10]. The number of tetrahedral elements in the coolant hole is increased as shown in
Figure 3b, with high resolution of the channel flow. The generated grids are particularly accurate in the near-wall region, where a value of
always below
is obtained. The final mesh used for the nominal case contains approximately
elements, which trebles the number of elements used for the analysis showed in [
10] for the same case. The main reason for the large increase of elements occurring in the present work is the refinement of the portion of control volume where shock/coolant interaction is expected. To include a smooth transition of the elements’ dimensions all over the domain, also the other regions have been refined. All the used meshes include a similar number of elements and are generated using the same rules for elements sizing. Since the simulations in [
10] have been demonstrated to be grid independent, the ones of the present work are also considered to be reliable for the present analysis.
The aims of the present numerical campaign can be summarized in the following list:
to correctly capture the characteristic flow features of the present test case;
to compare the deterministic distributions of obtained with different turbulence models;
to propose a reliable strategy for the calculation of non-dimensional performance parameters;
to compare the performance of two different polynomials bases in the evaluation of aleatory uncertainty;
to quantify the uncertainty associated with manufacturing deviations for the calculation of non-dimensional performance parameters (except for );
to see how uncertainty propagates depending on the turbulence model.
5. Data Processing
To evaluate the performance of film cooling for the current configuration, the parameters defined in the previous section must be calculated using the deterministic results of the numerical simulations. Considering the control volume defined in
Figure 1, it is possible to assign a specific number to each section:
section 1 is the main-flow inlet;
section 2 is the cooling-hole inlet section;
section 3 is the cooling-hole outlet section;
section 4 is the main-flow outlet section.
Among the necessary values for the calculation of the parameters, the ones in the following list are known:
fluid properties:
gas constant ;
the ratio between the specific heats ;
control volume geometrical characteristics:
main-flow inlet area ;
coolant inlet area ;
hole diameter D;
coolant hole section area ;
fluid dynamics boundary conditions:
stagnation pressure and stagnation temperature at main-flow inlet section;
stagnation pressure and stagnation temperature at coolant inlet section;
values obtained from each converged simulation:
coolant mass-flow ;
main-flow pressure calculated on the plate before the hole exit section ;
coolant pressure calculated on the plate before the hole exit section .
The isentropic relations and the perfect gas law are used to complete the definition of the non-dimensional parameters [
27]:
where
,
and
are the stagnation values,
T,
p and
are the static values and
is the Mach number. More information can be found in the nomenclature section. The Equation (
7) is exact for a steady, adiabatic process while Equations (
8) and (
9) are valid for a steady, isentropic process only. Under these hypotheses, the non-dimensional parameters can be calculated as follows.
5.1. Blowing Ratio
The blowing ratio has been defined in Equation (
1). The main-flow mass flux
is determined through Equation (
11):
where the
value is evaluated using Equation (
12):
where
is the stagnation temperature in section 1. The isentropic Mach number
is calculated using Equation (
13):
The coolant mass flux
is determined through Equation (
14):
5.2. Density Ratio
The density ratio has been defined in Equation (
4). The main-flow density value
is calculated using Equation (
15):
where the
is the one calculated using Equation (
13). The value of the stagnation density
is calculated using Equation (
16):
The coolant density value
is calculated using Equation (
17):
where the value of the stagnation density
is calculated using Equation (
18):
while the
is the one calculated using Equation (
19):
5.3. Momentum Ratio
The momentum ratio is defined according to Equation (
5), which can be rewritten as reported in Equation (
20):
The main-flow momentum
is calculated as reported in Equation (
21):
where the
value is calculated using Equation (
13). The coolant momentum
is calculated as reported in Equation (
22):
where the
value is calculated using Equation (
19). Since both the main-flow and the coolant are modelled using air as perfect gas, they have the same value for
. Therefore, considering Equations (
21) and (
22) the momentum ratio can be rewritten as in Equation (
23):
5.4. Discharge Coefficient
All the necessary values for the calculation of the discharge coefficient as defined in Equation (
6) has been defined in the previous sections. Then, it is sufficient to substitute
with
to obtain Equation (
24):
7. Uncertainty Quantification Methods for CFD
The effect of manufacturing deviations on film cooling performance in presence of shocks has been studied by means of a non-intrusive statistical analysis. The proposed approach is based on the assumption that the perturbation of either boundary conditions or geometric features propagate linearly in the solution. Deviations in the prediction of a System Response Quantity (SRQ) depends on: (i) a Probability Distribution Function (PDF) that describes the deviation of the aleatory parameter; () the adopted turbulence closure which defines the epistemic uncertainty; () the physics involved in the flow structures. The statistical evaluation of the selected SRQ is obtained by performing a set of deterministic simulations with different boundary conditions, which are representative of the PDF of the aleatory uncertainty.
The efficiency and the reliability of the UQ approach depends on the adopted sampling technique. In the context of uncertainty quantification used along numerical simulations, the sampling process consists of collecting a set of results of deterministic analyses obtained for different boundary conditions. That sampling could be obtained using Monte Carlo approach, which implies many simulations [
34]. As a consequence of the high computational demand, Monte Carlo procedure is not reliable when simulations are obtained using Computational Fluid Dynamics (CFD). A more suitable approach is based on the polynomial description of the stochastic variation. That approach allows for the determination of the PDF of the selected aleatory parameters by performing simulations only for a relatively small amount of boundary conditions. The investigated working points are obtained using the PCM procedure based on the representation of uncertainties on a linear space where the orthonormal basis is formed by a set of polynomials. According to this assumption, the propagation of an uncertainty parameter
associated with a PDF
within a physical system can be described in a polynomial form such as in Equation (
25):
In Equation (
25)
and
are deterministic functions and
is a set of multidimensional orthogonal polynomials in the random variables
. Such methodology is usually known as Polynomial Chaos (PC) [
35,
36]. An extensive review of the methods is detailed in [
37].
7.1. Hermite Polynomials
When aleatory uncertainty is associated with a normal distribution and the flow field is not affected by singularities, the basis
on which the space of solutions can be reconstructed is formed by Hermite polynomials. The generic expansion of the stochastic response in Equation (
25) can be expressed as reported in Equation (
26):
In a 1-D space, the form of the Hermite polynomials is reported in Equation (
27) up to the second order, where
is the
order polynomial. A detailed description of a stochastic analysis based on the use of Hermite’s polynomials for the study of aleatory uncertainty in turbomachinery-relevant problems is reported in [
14]. In that work the inlet flow velocity in an internal cooling system is assumed as aleatory uncertainty. The inlet velocity is associated with a normal PDF and the boundary conditions are selected according a PCM based on Hermite’s polynomials. Results demonstrated a correlation between aleatory and epistemic uncertainty in presence of non-linear problems dominated by turbulence. That kind of approach is the first one selected for the present activity.
7.2. Padè–Legendre Approximants
A different choice of polynomial basis is needed when the flow structure is characterized by discontinuities such as shocks, which is the case described in [
6]. Shock-dominated flow structures can be modelled by Padè approximation as shown in [
38]. The Padè approximation is a generalization of PC where a discontinuous response surface is described by a ratio of PC expansions. As this expansion can have poles, it can describe discontinuities. Numerator and denominator of the rational function are determined through a finite sum of orthogonal polynomials bases whose coefficients can be calculated from the function’s values for a predefined set of points. Examples of this include Padè–Jacobi, Padè–Chebischev and Padè–Legendre (PL) approximants. The latter represent the second approach used in the present work following the prescriptions found in [
38,
39].
A complete basis of Legendre polynomials is defined by Equation (
28):
Polynomials are defined uniquely as
and the expansion can be continued to any desired order. Every function
u can be defined as a linear combination of Legendre’s basis (Equation (
29)):
The coefficients
in the
Legendre coefficient of
u are defined as in Equation (
30):
The function is known in specific discrete points and the Gaussian quadrature provides the discrete formulation where each scalar product in Equation (
30) is defined as in Equation (
31):
The nodes
in Equation (
31) are the quadrature points and
are the quadrature weights of the Gauss–Lobatto quadrature rule as described in [
40]. The nodes are given as the roots of
. Consequently, the nodes are chosen as the roots of the polynomial
. The Gauss–Lobatto quadrature rule is chosen because it requires less function evaluations than standard Gauss quadrature. Given two integers
M and
L, a Padè–Legendre approximation of a function
u is the ratio of two approximating polynomials
and
based on the Legendre basis. The overall order of the reconstruction is
. Once the polynomial representation of the underlying function is obtained, it is possible to obtain mean
and variance
directly from the functions in Equations (
32) and (
33):
In Equation (
32)
is the random parameter associated with the PDF
,
represents the domain of the random parameters and
x is the physical space.
8. Discussion on the Statistical Analysis
The aim of this section is to evaluate how manufacturing deviations affect non-dimensional parameters of a cooling device in supersonic conditions depending on turbulence modelling and UQ methodology. Deterministic values of
,
,
and
are reported in
Table 2 and
Table 3. Response surfaces for the parameters are visible from
Figure 6,
Figure 7,
Figure 8,
Figure 9,
Figure 10,
Figure 11,
Figure 12 and
Figure 13. Finally, their mean values and standard deviations are reported in
Table 4.
Data included in
Table 2 and
Table 3 show the dependence of the parameters on both the cooling-hole diameter and the plate’s trailing-edge dimension. The
value is almost constant (between
and
) for the investigated configurations, which is an expected result. In fact, considering the definition of
in Equation (
24) and increasing the coolant channel diameter D
both the actual and the isentropic mass-flows increase by similar amounts, with a limited impact by the variation of thermodynamic losses. Concerning the other parameters,
shows a variation up to
with respect to the nominal value (which is the one obtained for the nominal configuration, Geom. 5) and
and
are characterized by variations in the order of
, irrespective of the selected turbulence model.
The outcome obtained using the
values is showed from
Figure 6a to
Figure 7b. As expected, the RSs are almost planar for all the investigated cases and there are limited differences between the SST and the RSM data. It can be observed that the selected approach based on the Hermite polynomials is limited to a second-order accuracy when applied to a set of 9 deterministic values, then the generation of a RS that could match all the deterministic values is not possible (see as example
Figure 6). On the contrary, Padè–Legendre approximants can generate discontinuous RSs at sixth order with the same set of values, which gives the outcomes showed in
Figure 7. The latter RSs accurately follows the deterministic values thanks to the superior ability to change the surface curvature. The same outcome can be observed by comparing
Figure 9 and
Figure 11 with
Figure 8 and
Figure 10 respectively. The ability of Padè–Legendre approximants in the generation of a RS that almost matches all the deterministic values obtained for
and
is guaranteed by the possibility to generate surfaces with quick curvature variations.
Concerning the
RSs showed from
Figure 12a to
Figure 13b, no remarkable differences are individuated by changing either the turbulence or the UQ methodology. The latter outcome is ascribable to the smooth distribution of the deterministic points. Looking at the RSs it is possible to note that
varies by a non-negligible factor (up to
as previously underlined) mainly depending on the D
value.
Generally speaking, the impact of the variation of on the selected parameters seems to be negligible for all the investigated cases. That outcome could be explained considering that the calculations have been performed with steady flow assumption, which neglects the possibility that a variation in the value would affect the shedding frequency and the time-averaged shock intensity and impact location. However, more analyses are necessary on a symmetric control volume with periodic conditions to verify that guess.
From a technological point of view, it is of major interest to determine the mean value and the standard deviation (STD) associated with the statistical analysis. For that reason, the values obtained for the nominal case and using the two UQ methodologies are reported in
Table 4. Values obtained using the UQ methodologies are quite close (within an accuracy of
) to the deterministic values. What is primarily interesting are the associated STD values that represent the numerical uncertainty associated with manufacturing deviations. As can be observed, STD values associated with
,
and
are below
. Concerning
, STD values are between
and
, thus demonstrating a non-negligible impact of manufacturing deviations on that fundamental parameters. It is worth underlining that the values obtained using Hermite polynomials and Padè–Legendre approximants are quite close with each other. It can be concluded that despite the higher ability in the reconstruction of RSs, to use Padè–Legendre approximants does not provide a relevant improvement in the evaluation of mean values and STDs (at least in a case where the RSs do not present steep curvature variations).
9. Conclusions
The main topic of the present activity is the evaluation of the impact of manufacturing deviations on the performance parameters (, , and ) of a film cooling configuration in presence of oblique shocks. A total amount of 18 steady simulations are used to take into consideration the geometrical uncertainty associated with the cooling-hole diameter and the dimension of the central plate trailing edge. Two different turbulence closures (SST and RSM) are used to underline the impact of the anisotropy on the calculation of the 2-D maps. Two different UQ methodologies based on Hermite polynomials and Padè–Legendre approximants are used to determine the most suitable approach.
Concerning the deterministic results obtained for the nominal configuration, it is observed that the coolant is mainly confined close to the centerline and its mixing with the main flow is underestimated if compared with the available experimental map. The use of the RSM approach for turbulence modelling does not provide significant improvements in the calculation of
. Those outcomes are in line with computations by other authors [
29] and are associated with the under-prediction of the vertical mixing in the wake of a jet for high
values. Also, the
and
parameters are compared with the experimental values and are found to be within the experimental accuracy and in line with literature data.
Concerning the UQ analysis, it is observed that , and are not almost unaffected by the selected manufacturing deviations while greatly depends on D. No relevant impact is found using either SST model or RSM in combination with Hermite polynomials and Padè–Legendre approximants. The latter can generate Reynolds stresses that match the deterministic points with greater accuracy with respect to Hermite polynomials, but it does not mean that relevant changes are found in the calculation of mean values and STD. In fact, all the data obtained using the UQ methodologies are close (within an accuracy of ) to the corresponding deterministic ones. STD values for are between and , thus confirming a non-negligible impact of manufacturing deviations.