## Abstract

In this paper we investigate the performance of eddy viscosity turbulence models for high-pressure nozzle guide vane (NGV) flows with engine-realistic turbulence. Using metrics such as integrated kinetic energy (KE) loss and mixing rate, we compare simulation results using turbulence models with high-fidelity experimental data from fully-cooled NGVs, operating at engine-matched conditions of Reynolds number and Mach. For the widely-used k–ω shear stress transport (SST) model, the aerodynamic and thermal wakes of the NGVs were undermixed. We compare simulation results with those using other common turbulence models including the standard k–ω, baseline k–ω, Spalart–Allmaras, and baseline explicit algebraic Reynolds stress model. Alternative turbulence models formulations are explored, with an emphasis on modifications to the implementation of the shear stress limiter. We also investigate the sensitivity of models to the specified inlet turbulence intensity. We demonstrate that predictions of integrated KE loss, as well as the decay trends of the aerodynamic and thermal wakes are a closer match to experimental data for the baseline algebraic Reynolds stress model than for other turbulence models in more common use for NGV predictions.

## 1 Introduction

Although recent advances in numerical algorithms and an ever-increasing availability of computing resources have facilitated one-off large eddy simulations of fully-cooled turbine stages (see, e.g., Ref. [1]), such simulations remain prohibitively expensive for use in design practice, meaning that engineers continue to rely on Reynolds-averaged Navier–Stokes (RANS) models to predict the performance of turbine components. RANS models rely on experimental data for calibrating closure coefficients. Calibration is often performed using measurements from flows with very low freestream turbulence, leading to unpredictable behavior in high-turbulence flows such as those encountered at the combustor–turbine interface and in the high-pressure (HP) nozzle guide vane (NGV).

For rich-burn combustor architectures, based on the limited data available from reacting gas turbine combustors (e.g., Refs. [2,3]) turbulence intensities between 10% and 30% have been reported at the HP NGV inlet plane. Interestingly, these studies reported only a mild sensitivity of the turbulence intensity (TI) at combustor outlet to the presence of combustion as opposed to geometry. Recent literature from non-reacting low-temperature experiments [4–6] for rich-burn architectures shows average TI levels in the range 10–25%. The relative ease of instrumenting low-temperature experiments has also facilitated detailed investigations of turbulence length scales at combustor outlet. Barringer et al. [4], Cha et al. [5], and Folk et al. [6] all reported length scales similar to the combustor dilution port diameter. For typical configurations, this corresponds to integral length scales around 0.20–0.40 times the NGV axial chord length.

Turning to lean-burn combustors, Schroll et al. [7] measured turbine inlet turbulence intensities between 8% and 20% under reacting conditions. A number of low temperature studies have successfully used non-reacting combustor simulators to generate lean-burn representative swirl and temperature profiles (see, e.g., Refs. [8–10]). In similarity to low temperature rich-burn simulators, non-reacting simulators allow for more detailed characterization of turbulence than their reacting counterparts (in which measurements are invariably limited). Peak turbulence intensities around 25% were reported by Bacci et al. [10]. For the same combustor simulator, Koupper et al. [11] measured integral length scales equal to approximately 10% of the swirler diameter. Further studies are necessary to determine the factors affecting length scales in lean-burn combustors, though intuitively we might expect the size of the turbulent structures to depend on both the swirler blade spacing and the overall size of the swirler.

For the *k*–*ω* shear stress transport (SST) model [12], the recommended *maximum* turbulent-to-molecular viscosity ratio at inlet is 10^{−2}. This value is several orders of magnitude smaller than the values typical for combustor-outlet flows (see, e.g., Cha et al. [5], who measured turbulent-to-molecular viscosity ratios around 10^{4}). This is one example of a high-turbulence flow which falls well outside the calibration range of most turbulence models. It is of practical interest for designers to understand how different turbulence models behave in these high turbulence environments.

To understand the behavior of various common turbulence models in a high-turbulence flow environment, we compare the standard *k*–*ω* (SKO) model [13], baseline *k*–*ω* (BSL) model [12], *k*–*ω* shear stress transport (SST) model [14], the baseline explicit algebraic Reynolds stress model (BSL-ARSM) [15,16], and the Spalart–Allmaras (SA) model [17] in terms of their ability to accurately predict the aerodynamic and thermal field downstream of fully-cooled NGVs with high inlet turbulence. The computational fluid dynamics (CFD) results are compared with experimental data from an annular cascade of film-cooled parts operating at engine-matched Reynolds and Mach number.

## 2 Shear Stress Transport Model

The *k*–*ω* shear stress transport (SST) turbulence model is the de facto industry standard for most turbomachinery flows. The popularity of the SST model [14] can be attributed—in part—to its zonal treatment of boundary layer and free stream regions. The SST model uses SKO [13] near the walls (log-layer and viscous sublayer), and the standard *k*–*ɛ* (SKE) [18] model in the freestream. The zonal approach is advantageous because SKO is reasonably accurate for boundary layer flows but exhibits unrealistic sensitivity to turbulence boundary conditions (see, e.g., Ref. [19]). The *k*–*ɛ* model is considered more robust in terms of its sensitivity to inlet conditions of turbulence but is less accurate for wall-bounded shear because it relies on wall functions to model the viscous sublayer. The zonal approach of the *k*–*ω* SST model exploits the strengths of two constituent models in the regions in which they are most applicable. We now briefly consider the fundamentals of the underlying SKO and SKE models.

### 2.1 Eddy Viscosity and Turbulence Scales.

*τ*

_{ij}, is given by

*S*

_{ij}, is given by

*μ*

_{t}is eddy viscosity, ∂

*U*

_{k}/∂

*x*

_{k}is the velocity divergence,

*ρ*is volumetric mass density,

*k*is turbulence kinetic energy (TKE; defined as $uiui\xaf/2$), and

*δ*

_{ij}is the Kronecker delta.

*U*and

*u*are the mean-flow and the fluctuating velocity components, respectively, using Einstein notation for indexing. By convention,

*τ*

_{ij}is equal to the

*negative*Reynolds stress tensor, $\rho uiuj\xaf$, where the overbar represents Favre averaging. Equation (1) closes the RANS equation system by relating the turbulent stresses to the mean-flow strain rates through

*turbulent eddy viscosity*,

*μ*

_{t}. The eddy viscosity represents momentum diffusivity through turbulence, in analogy to the dynamic viscosity which represents the momentum exchange though non-turbulent shear. In analogy to the kinematic viscosity, it is possible to define a

*kinematic eddy viscosity*,

*v*

_{t}, where

*v*

_{t}=

*μ*

_{t}/

*ρ*. Prandtl proposed calculating

*v*

_{t}using characteristic properties of the largest turbulent eddies. Based on dimensional arguments, kinematic eddy viscosity (units of [m

^{2}s

^{−1}]) could be regarded as proportional to the product of a characteristic turbulence velocity scale,

*v*(with units of [m s

^{−1}]), and a turbulence length scale,

*ℓ*(with units of [m]). In this abstraction, turbulence length scale might be viewed as the mean eddy diameter, and the turbulence velocity scale may be thought of as the time-average

*unsteady*velocity of the fluid circulating in the eddy. Linear scaling of diffusivity with eddy size and circulation velocity makes sense because larger eddies (i.e.,

*ℓ*) will exchange mean flow momentum between two more distant points, and a greater circulation velocity (i.e.,

*v*) increases the time rate of transfer of momentum (i.e., diffusion). In the SKO model kinematic eddy viscosity is given by

*k*is the TKE (units of [m

^{2}s

^{−2}]), and

*ω*is the specific turbulence dissipation rate (units of [s

^{−1}]). In continuation of the earlier dimensional arguments,

*k*can be thought of as the square of the turbulence velocity scale, and

*ω*can be viewed as the reciprocal of a turbulent timescale. As ever with dimensional arguments, kinematic eddy viscosity can be expressed using many combinations of variables which relate to the

*energy*and

*scale*of the turbulence. In the SKE model [18] kinematic eddy viscosity is expressed as

^{2}s

^{−3}]). Equation (4) is based on the observations of Taylor [20], who showed (using hotwire measurements of decaying grid turbulence) that

*ℓ*was proportional to $k1.5/\u03f5$.

*k*and

*ω*, using the identity

*k*and $\u03f5$ into a model for the transport of

*k*and

*ω*(see full derivation in Ref. [21]). We now consider the transport equations of the combined SKO–SKE model.

### 2.2 Shear Stress Transport Equations.

*U*

_{i}/∂

*x*

_{j}is the velocity Jacobian,

*D*/

*Dt*is the material derivative,

*ρ*is volumetric mass density,

*μ*is dynamic viscosity, and

*σ*

_{k}is a closure coefficient. On the right-hand side of Eq. (6), turbulence production term (II) represents the energy gained by turbulent eddies through interactions with the velocity gradients of the mean-flow, and dissipation term (III) represents the dissipation of TKE at the Kolmogorov microscale. Term (IV) models the combined effects of laminar viscous diffusion, pressure diffusion, and turbulent diffusion of TKE. We now consider the turbulence production term (II) in detail.

*τ*

_{ij}, we get

*S*, is a coordinate-independent scalar measure of the strain rate, given by

*C*is nominally equal to unity but is often treated as an empirical correction factor, which is designed to prevent excessive turbulence production through shockwaves (default value of 3; see cfx theory guide [22]). For divergence-free (i.e., incompressible) flow, it can be shown that Eq. (9) is equal to the eddy viscosity times the square of the strain invariant,

*S*. This simplification is convenient for discussion purposes, where the compressible definition (used in simulations) is less intuitive.

*ω*. We recall that a physical interpretation of the specific dissipation rate is that its reciprocal is proportional to the turbulence time scale. Specific dissipation rate,

*ω*, can be modeled using the following

*empirical*transport equation

*α*,

*β*, $\sigma \omega $, and

*σ*

_{ω2}are closure coefficients. In analogy to the TKE equation, the production (II), dissipation (III), and diffusion terms (IV), appear on the right-hand side. Term (V) represents cross-diffusion of turbulent kinetic energy and specific dissipation rate. Cross-diffusion term (V) appears because of the transformation of the transport equation for $\u03f5$ (see Ref. [18]) into a transport equation for

*ω*. The cross-diffusion term alleviates the well-known problem of the SKO model, where a small change in the freestream value of

*ω*can give rise to a non-physical change in the solution (see, e.g., Refs. [12,23]). Equation (11) is justified by dimensional considerations and empirical calibration (see discussion by Wilcox [24]), and physical interpretations of the individual terms do not necessarily exist. This contrasts with the equation for TKE transport (Eq. (6)), which can be derived directly from the exact RANS equations, even though in practice closure of most terms involves modeling assumptions.

### 2.3 Model/Closure Coefficient Blending.

*ϕ*represents a closure coefficient (e.g.,

*σ*

_{k}), which varies from

*ϕ*

_{1}(SKO regime coefficient

*σ*

_{k1}) to

*ϕ*

_{2}(SKE regime coefficient

*σ*

_{k2}) as a blending function,

*F*

_{1}, varies from 1 to 0.

*F*

_{1}is defined as

_{1}is given by

*y*is distance to the nearest wall, and where

*CD*

_{kω}is the cross-diffusion term from Eq. (11) defined by

^{−10}is imposed to ensure

*CD*

_{kω}remains positive. Blending function

*F*

_{1}increases smoothly from 0 (SKE-regime) to 1 (SKO-regime) as arg

_{1}varies from 0 to infinity (

*F*

_{1}> 0.99 for arg

_{1}> 1.28). Consider now the definition of arg

_{1}. The first sub-argument (I) of arg

_{1}represents the turbulence length scale $(\u2113=k/\beta *\omega )$, divided by the distance to the nearest wall,

*y*(see Ref. [12]). This term is meant to

*activate*the SKO-regime (i.e.,

*F*

_{1}= 1) in the log-layer, based on the empirical observation that

*ℓ*≈ 2.5

*y*in the log-layer. The second sub-argument (II) of arg

_{1}is tailored to

*activate*the SKO-regime (i.e.,

*F*

_{1}= 1) in the viscous sublayer, based on the empirical observation that 1/

*ωy*

^{2}is approximately constant in this region [12]. The third sub-argument (III) of arg

_{1}is designed to

*deactivate*the blending function (i.e.,

*F*

_{1}= 0) for large positive values of

*CD*

_{kω}, which are considered indicative of free shear (i.e., SKE-regime) flow. In addition to interpolating the closure coefficients,

*F*

_{1}is also used to deactivate the cross-diffusion term in Eq. (11) near the walls, where it is not valid (for the same reasons as the SKE model). This is achieved multiplying the cross-diffusion term (V) in Eq. (11) by (1 −

*F*

_{1}).

### 2.4 Shear Stress Limiter.

*realizability*limiter), from which the model derives its name. The shear stress limiter improves predictions of separation in flows with moderate adverse pressure gradients. The apparent cause for the poor predictions of separation in models

*without*a stress limiter is the overprediction of turbulent shear stress in the outer boundary layer [12]. The shear stress limiter—designed to prevent these overpredictions—is implemented by modifying the algebraic expression for eddy viscosity

*a*

_{1}is the Bradshaw constant (

*a*

_{1}= 0.31; see Ref. [25]),

*S*is the familiar strain invariant, and

*F*

_{2}is a second blending function defined as:

_{2}is given by

_{2}is equal to 2

*ℓ*/

*y*, and the second sub-argument of arg

_{2}is identical to that of arg

_{1}. The first sub-argument of

*F*

_{2}ensures that

*F*

_{2}remains active (i.e.,

*F*

_{2}= 1) slightly further into the freestream than

*F*

_{1}, which ensures that the stress limiter can operate in the outer boundary layer. We now consider the operation of the stress limiter (i.e., Eq. (16)) in detail.

Outside of the boundary layer (i.e., *F*_{2} = 0), the SKO model definition of eddy viscosity is used (*μ*_{t} = *ρk*/*ω*; see Eq. (3)). When *F*_{2} = 1 (i.e., within the boundary layer), the limiter is active when *S* > *a*_{1}*ω*. Taking the square of both sides, and multiplying by *μ*_{t} (the inequality holds because all variables are positive), it can be shown (using Eq. (3) and the approximations that $P\u2248\mu tS2$ and $a1\u2248\beta *$) that the inequality is equivalent to the condition that $P>E$. The condition $P>E$ is considered indicative of an adverse pressure gradient [12]. To understand how Eq. (16) limits the shear stress, we consider the boundary layer momentum equations, for which the primary shear stress can be approximated as *τ*_{xy} ≈ *μ*_{t}*S* (where *x* is the streamwise direction and *y* is the wall-normal direction). Substituting *τ*_{xy} into Eq. (16), we see that when *F*_{2} = 1, the shear stress obeys *τ*_{xy} ≤ *ρa*_{1}*k*, where *a*_{1} is taken from the experiments of Bradshaw et al. [25]. Accurate detection of boundary layers is essential for correct function of the limiter because the inequality is not valid in the freestream.

An interesting historical detail is that the *original* SST model [12] (often referred to as SST1994) used the magnitude of vorticity, *Ω*, instead of the strain rate invariant, *S*, for limiting shear stress in Eq. (16). This was found to be problematic, however, because *Ω* amplifies numerical noise, resulting in a speckled eddy viscosity distribution (see, e.g., Ref. [26]). An additional advantage of using an *S*-based limiter is that the inverse relation between eddy viscosity and the strain rate ensures *realizability* of the turbulent stresses. *Realizability* implies that the turbulent stresses are mathematically achievable (see Ref. [27]). This imposes the constraint that each principal stress of *τ*_{ij} is no greater than zero (considering that the autocorrelation $\u2212\rho uu\xaf$ cannot be positive) and no less than −2*ρk* (considering that $\u2212\rho uiui\xaf=\u22122\rho k$). Presumably, it is for these two reasons that the *S*-based limiter is more widely used in modern solvers [14]. We refer to the model using a strain-rate-based stress limiter as SST2003 (to distinguish from SST1994). As we will see in the current study, the stress limiter of the SST2003 model is too strict, and results in undermixing for engine-realistic turbulence.

### 2.5 Standard *k*–*ω* and Baseline Model.

We will not formally introduce the SKO and BSL transport equations, but note that from Eqs. (6)–(18), the SKO model [13] can be defined by setting *F*_{1} = 1 (constant) and by calculating eddy viscosity using *μ*_{t} = *ρk*/*ω* (i.e., by disabling the stress limiter and omitting the cross-diffusion term). Model coefficients are also modified (see Ref. [13]). Similarly, the baseline (BSL) model [12] can be defined by disabling the stress limiter, setting closure coefficient *σ*_{k1} = 2 (versus 1.176 for SST), but retaining the blended zonal treatment and cross-diffusion term. We perform simulations using these simplified models to investigate the effects on NGV flows of model blending and the shear stress limiter.

### 2.6 Production Limiters.

*stagnation point anomaly*is explained by the linear nature of the Boussinesq hypothesis, by which normal strains erroneously contribute to TKE production (see explanation of Moore and Moore [27]). For NGV flows, production limiters are also relevant for the flow around the suction peak, where high strain rates otherwise result in TKE over-production (recall $P\u2248\mu tS2$). We use a production limiter with $Clim=10$, unless explicitly stated otherwise.

## 3 High Turbulence Flat Plate Test

The flat plate test is perhaps the simplest example of boundary layer flow and has been used extensively for turbulence model calibration. Good likeness to experimental data in the flat plate test can be regarded as a prerequisite for accurate predictions in more complex flows. We briefly consider the literature on flat plate tests specific to the effects of freestream turbulence.

Folk et al. [6] recently investigated boundary layer development for low and high freestream TI at the same length scale (approximately 40 mm). The study builds on earlier work of Hancock and Bradshaw [29], as well as Thole and Bogard [30], who reported elevated skin friction at high TI, but did not report/control the momentum thickness Reynolds number, $Re\theta $. Folk et al. [6] used a turbulence grid and a rich-burn combustor simulator, to generate turbulence intensities of approximately 2.3% and 18.7%, reporting a 20% increase in skin friction coefficient and a 22% increase in dissipation coefficient at high TI (for the same $Re\theta $). These effects were shown to result from deep penetration of freestream turbulence into the boundary layer.

The flat plate data of Folk et al. [6] are plotted in Fig. 1 alongside equivalent simulations using the SST2003 turbulence model, which were performed based on the test setup description in Ref. [31]. We now review the accuracy of SST2003 for predicting the skin friction coefficient, *C*_{f}, and an auxiliary boundary layer shape factor (ratio of momentum-to-energy thickness), *H*_{23}. The auxiliary shape factor is relevant because it can be used to estimate the dissipation coefficient, *C*_{d}, using the approximation *C*_{d} ≈ *C*_{f}/4*H*_{23} (see Ref. [6]). The CFD-predicted velocity profiles exhibit *peaking* behavior close to the leading edge of the flat plate, as well as freestream acceleration due to blockage. In the current study, the reference velocity was defined as the peak velocity of the local profile, and displacement, momentum, and energy thickness were calculated by integrating up to the point with peak velocity.

First consider the skin friction coefficient, *C*_{f}, which is shown in Fig. 1(a). The experimental data of Folk et al. [6] indicates a 20% increase in skin friction (average across all points) going from low TI (2.3%) to high TI (18.7%). By comparison, SST2003 CFD predicts an increase of 35% (on average over the range from $1000<Re\theta <2000$). The overall agreement between experiments and CFD is excellent, with an average overprediction of 2.8% at low TI and an average overprediction of 16.1% at high TI. A better match at low TI is expected, because turbulence models are typically calibrated using low TI data.

Now consider the auxiliary shape factor (ratio of momentum-to-energy thickness), *H*_{23}, which is shown in Fig. 1(b). The experiments of Folk et al. [6] indicate a 4% reduction of *H*_{23} for the high TI case. A similar result is given by CFD, but the shape factor is lower than corresponding experimental data by, on average, 2.7%. Once again, the agreement between experiments and CFD is excellent at low TI, and acceptable at high TI.

Folk et al. show that the dissipation coefficient can be approximated as *C*_{d} ≈ *C*_{f}/4*H*_{23}. Experimentally-measured and CFD-predicted dissipation coefficient data (using the approximated formula) are shown in Fig. 1(c). At low TI, the CFD-predicted dissipation coefficient is—on average—accurate to within 3.1% of experimentally-measured values. Going from low to high TI, *C*_{d} increases by 22% in experiments, and by 46% in SST2003 CFD. At high TI, we therefore expect SST2003 to overpredict kinetic energy loss by 20% for a given $Re\theta $. These results are promising, because they show that SST2003 is both relatively *accurate* in absolute terms and predicts the correct *sensitivity* to freestream TI. Average dissipation coefficient data for all turbulence models considered in the current study are summarized in Appendix A.

## 4 High Turbulence Nozzle Guide Vane Flows

Amend et al. [32] and Burdett and Povey [33] used the SST2003 model to simulate the flow of fully-cooled NGVs with turbulence of high intensity (13%) and large length scale (approximately equal to the NGV axial chord) at inlet. Both studies investigated the spreading and decay rates of the aerodynamic wake, and reported significant undermixing, when comparing SST2003 results to high-fidelity experimental validation data. Simplistically, an underprediction of spreading and decay rates indicates an underprediction of and/or incorrect distribution of eddy viscosity. Consider the contours of eddy viscosity, *μ*_{t}, from the SST2003 simulations of Amend et al. [32], which are shown in Fig. 2 (cross-section at NGV midspan). The cascade inlet conditions were characterized using hot-wire anemometry, pressure probes, and thermocouples (for details see Ref. [32]). The integral length scale at cascade inlet was of the order of one axial chord (1.0 *C*_{x}), and the TI was approximately 13%. Most studies in which rich-burn combustor simulators have been used report integral length scales in the range 0.2 *C*_{x}–0.4 *C*_{x} (see, e.g., Refs. [4–6]). In comparison to these studies we have a relatively large length scale, and this makes it particularly interesting for exploring limitations of existing turbulence models at conditions of both high turbulence and high length scale (with correspondingly lower dissipation rate with axial distance).

The most striking feature in Fig. 2 is a sudden reduction of eddy viscosity along a sinusoidal front, located approximately 0.8 *C*_{x} upstream of the leading edge. This reduction is *not* the result of a transport phenomenon, but rather is *algebraically enforced* through the shear stress limiter of the SST model (see Eq. (16)). We recall that the shear stress limiter was designed to improve predictions of separation in flows with strong adverse pressure gradients (see Ref. [12] for details). This is achieved by limiting eddy viscosity in the outer boundary layer, on the condition that $P>E$. The region of *intended* operation of the limiter is demonstrated on the uncovered suction side (SS). For the current case, the shear stress limiter is also *erroneously* active from far upstream of the leading edge, through the NGV passage, and extending up to the geometric throat. Erroneous activation throughout the vane passage (i.e., well beyond the outer boundary layer) is related to the combined effects of relatively large length scales in the freestream (on the order of 1.0 *C*_{x} at domain inlet) and high cross-passage strain (limiter active for *SF*_{2} > *a*_{1}*ω*, with *S* high on account of the flow field, and *F*_{2} = 1 due to the large turbulence length scales in the freestream). In effect, the limiter artificially shields the near-wall region (on the vane surface) from incoming freestream turbulence, resulting in undermixing of the thermal mixing layer and aerodynamic boundary layer. The shielding phenomenon also explains why *k–ω* SST is sometimes reported to have a very low sensitivity to changes in the inlet turbulence. This may be desirable for external aerodynamic flows, but presents a problem in turbomachinery applications, where incoming turbulence is generally higher, and known to affect heat transfer and kinetic energy loss.

The first blending function is also sensitive to the large freestream length scales (see arg_{1} in Eq. (14)), such that SKO-regime is active for approximately half of the NGV passage. As such, the intended zonal treatment (SKO in boundary layer; SKE in freestream) is not achieved. The effect is secondary to the more pressing shielding effect of the stress limiter, because its influence on results appears to be relatively small in comparison, as can be verified by modifying the first blending function.

## 5 Modified Blending and Stress Limiter

For engine-realistic turbulence, neither the closure coefficient blending function, *F*_{1}, nor the shear stress limiter of the SST2003 model function as intended. The result is that almost the entire NGV passage is treated in the same way as the boundary layer. Both problems are related to the large turbulence length scales in the freestream. Remedial model modifications and alternatives are now considered.

### 5.1 Modified Closure Coefficient Blending.

_{1}(see Eq. (14)) can cause unreliable detection of boundary layers. This occurs because the large freestream length scales (representative of combustor outlet flow) render the log-layer indistinguishable from the freestream. The problem can be mitigated by using an alternate definition for

*F*

_{1}, where only the viscous sublayer is modeled using SKO (with freestream

*and*log-layer modeled using SKE). Menter [12] suggests the following modified form for arg

_{1}

Equation (20) is equivalent to Eq. (14), but with the length scale-dependent term removed. We implement this modification in cfx with a UserFortran junction box routine, which overwrites vertex values of *F*_{1} (files and instructions available on request). Results using this modification are referred to as SSTF1V (i.e., SST with *F*_{1} active in the viscous sublayer only). Results for this modification are presented in Appendices A and B.

### 5.2 Modified Shear Stress/Realizability Limiter.

Now consider the *erroneous* activation of the shear stress limiter in the freestream (see Fig. 2). Problematic behavior of the shear stress limiter has recently been reported by Wei [34], who found that stress limiter-induced undermixing could lead to premature corner separation in compressors. Wei proposed disabling the shear stress limiter entirely, reporting improved predictions using this method. In doing so, Wei effectively recovered the BSL model with slightly different closure coefficients (*σ*_{k1} = 1.176 instead of 2.0).

*original*SST model [12], which uses the magnitude of vorticity,

*Ω*, instead of the strain rate invariant,

*S*, i.e.

Compared to Eq. (16), Eq. (21) reduces the artificial suppression of eddy viscosity in the freestream, because *Ω* ≈ 0 in the freestream, whereas *Ω* ≈ *S* in boundary layers. In theory, the *Ω*-based shear stress limiter should therefore retain its *intended* operation in boundary layers, while eliminating the *erroneous* activation in the freestream. Where we present results using this modification we refer to them as SST1994. For consistency with the original SST model from 1994, the *Ω*-based limiter is paired with a production limiter clipping factor $Clim=20$ (versus $Clim=10$ in SST2003). Another minor difference is that the original (SST1994) model applies the production limiter to the TKE equation *only*, whereas cfx limits the production terms in *both* the TKE and specific dissipation rate equations (in accordance with SST2003). An SST1994-type limiter can be implemented in cfx by adding a source term in the *ω* equation, which makes up the difference between the *limited* and the *unlimited* source term. This appears to make little difference in practice, however, because production-to-dissipation ratio (as modeled) exceeds 20 only very locally (SS peak) and by a small margin.

*S*

_{ij}), and ∂

*U*

_{k}/∂

*x*

_{k}is the familiar velocity divergence. The inequality follows from the compressible Boussinesq approximation in principal coordinates (i.e., coordinates for which

*S*

_{ij}is purely diagonal). The condition is enforced using the same clipping approach as is used in Eq. (16). We refer to this modification as SSTPSR (SST with exact principal stress realizability).

*directly*increases computational cost. It may therefore be preferable to

*estimate*the largest principal strain, $\lambda max$. Durbin [28] derived an

*upper bound*for the largest eigenvalue of the strain rate tensor, $\lambda max$. For incompressible flow, realizability is guaranteed for

*S*which is a factor of $2$ smaller than in the current study; the inequality has been scaled accordingly for consistency with Eq. (10). The condition in Eq. (23) is equivalent to the familiar SST2003 shear stress limiter (see Eq. (16)) with $a1=1/3\u22480.577$, provided that

*F*

_{2}= 1 (true for large length scales). We refer to this modification as SSTDR (Durbin realizability). For boundary layers in particular, Bradshaw et al. [25] suggest that

*a*

_{1}= 0.31 (this latter value being adopted in the SST model). This does not necessarily contradict the

*inequality*of Durbin [28]. However, since the Bradshaw assumption does not hold in the freestream, we might expect that the SST2003 shear stress limiter is too strict in this region. Setting $a1=1/3$ may partially alleviate this issue, although this can still result in condition that is

*locally*too restrictive, considering that the value for $\lambda max$ is conservatively estimated as its

*upper bound*. Park and Park [36] derived a

*lower bound*for $\lambda max$, which is equivalent to $a1=4/3$. This method does not reduce eddy viscosity below the realizability limit, but also does not guarantee realizability. Such a model is conceptually useful, but arguably of limited practical value. This realizability modification was not tested.

In summary, we consider the following approaches for mitigating the artificial suppression of turbulence in SST2003:

Disabling the stress limiter completely (recovering the BSL and SKO models; tested with $Clim=10$).

*Ω*-based limiter (tested with $Clim=20$ for consistency with the original SST1994 model [12]); we refer to this as SST1994.Principal stress realizability (Eq. (22); tested

*without*a production limiter); we refer to this as SSTPSR.Durbin realizability [28] (equivalent to $a1=1/3$; tested

*without*a production limiter); we refer to this as SSTDR.

Methods (iii) and (iv) were tested *without* a production limiter because they are inherently realizable (SSTDR only where *F*_{2} = 1). Considering that turbulence production scales with the difference of the principal stresses (see Ref. [27]), we expect that this approach will result in a slight over-production of turbulence, because experimental data shows that the principal stresses obey *τ* ≤ −0.2*ρk* [27] for plane-strain problems, rather than the slightly less restrictive *mathematical* condition of realizability, *τ* ≤ 0. The point is that the resulting turbulent stresses are formally *realizable*, but not necessarily *physically* achievable in practice. This is clearly a compromise but avoids sensitivity to the somewhat arbitrary production limiter clipping factor.

## 6 Additional Turbulence Models

The baseline algebraic Reynolds stress model, and the Spalart–Allmaras turbulence model were also tested. Both turbulence models were used in their documented forms only; for brevity we therefore introduce them only at a high level.

### 6.1 Explicit Algebraic Reynolds Stress Model.

*τ*

_{ij}is the turbulent stress tensor,

*a*

_{ij}is the anisotropy tensor, and

*δ*

_{ij}is the Kronecker delta. The anisotropy tensor,

*a*

_{ij}, is modeled using a polynomial of the strain rate (

*S*

_{ij}) and vorticity (

*Ω*

_{ij}) tensors (see Wallin and Johansson explicit ARSM [15]), which is derived from the exact Reynolds stress model equations by neglecting the advection and diffusion terms. The non-linear stress–strain relation (i.e., Eq. (24)) is typically tailored to guarantee stress realizability, which eliminates the need for additional production and realizability limiters.

A modified version of the Wallin and Johansson model [15] is available in cfx, which combines a surrogate BSL model with a third-order tensor polynomial for *a*_{ij} [16]. We refer to this model as BSL-ARSM. Given that BSL-ARSM is inherently realizable, it is perhaps surprising that cfx defaults to a production limiter $(Clim=10)$ with this model. In the current study, the production limiter was disabled for BSL-ARSM.

### 6.2 Spalart–Allmaras Model.

The Spalart–Allmaras (SA) turbulence model [17] is an empirically-motivated transport equation for an eddy viscosity-like variable. SA is available as a partially-developed feature in ansys cfx. We include results from the SA model due to its historical use in turbomachinery applications. Results using the *default* wall-resolving model version are presented (as opposed to the version using wall functions, which performed poorly on the flat plate test case, and underpredicted loss for the NGV test case). Of note, SA does not model the decay of turbulent viscosity in the freestream, which is considered advantageous for external aerodynamics, but might be problematic for internal flows. For the current NGV test case this may be acceptable, however, considering that large length scales result in *naturally* low decay rates in the freestream.

## 7 Experimental Data and Setup

The experimental data of this study has previously been reported in Ref. [32], which considered the flow through an annular cascade of fully-cooled high-pressure NGVs at a Reynolds number of 7.5 × 10^{5} (based on NGV axial chord) and a Mach number of 0.92. We briefly review the test section layout and instrumentation. For a detailed description, see Refs. [32,37].

### 7.1 Test Section Layout.

A meridional view of the test section is shown in Fig. 3. The mainstream flow passes through a turbulence grid into a contracting inlet duct upstream of the NGVs. Both the vane airfoil and the endwalls are film-cooled. The airfoil coolant and hub-endwall coolant are supplied from a plenum at the hub. The casing endwall coolant is supplied from a separate plenum at the casing. Downstream of the NGVs, the partially mixed flow exhausts to atmospheric pressure through a parallel constant-area exit duct, which extends approximately three axial chords downstream of the NGV trailing edge before exhausting to a large plenum.

### 7.2 Instrumentation.

Detailed traverse data of total pressure, total temperature, TI, and integral length scale were taken at the mainstream inlet plane (plane 1; see Fig. 3). At plane 1, total pressure and total temperature were also monitored continuously using fixed Pitot total pressure rakes and thermocouple rakes. The mainstream mass flow rate was measured using a high-precision calibrated sonic nozzle (see high-level facility layout in Ref. [37]). Coolant stream mass flow rates were metered individually using sonic nozzles. Total pressure and total temperature in the hub and casing coolant plena were monitored. The NGV wake was characterized at three downstream measurement planes, which were located at 0.25 *C*_{x} (plane 2), 0.50 *C*_{x} (plane 3) and 0.75 *C*_{x} (plane 4) downstream of the midspan trailing edge. At each downstream plane, the flow was surveyed over two vane pitches circumferentially, and from hub to case radially. An aerodynamic five-hole probe was used to measure flow angle, total pressure and Mach number, enabling the calculation of mass-flux-weighted averages. A thermocouple probe was used to survey the downstream thermal field.

### 7.3 Cascade Inlet Conditions.

Experimentally-measured and CFD-predicted conditions of total pressure, total temperature, turbulence intensity, and length scale at plane 1 are shown as radial profiles in Fig. 4. The experimental profiles were obtained by circumferentially averaging area traverse data from a 34-deg sector. The probes and measurement repeatability are discussed in Ref. [38]. Due to a probe size limitation, a small region near the hub could not be traversed. Similarly, the turbulence data near the casing needed patching due to a local unsteady interaction with the probe access port. The patching method was informed by a *k–ɛ* CFD simulation of the upstream flow conditioning system and turbulence grid (for details see Ref. [32]).

First consider the experimentally-measured profiles in Fig. 4. Normalized total pressure (*p*_{0,norm}) was very uniform, with a variation of ±0.1% about the mean. Normalized total temperature (*T*_{0,norm}) was also relatively constant across most of the span, except for the near-wall regions, where thermal boundary layers developed due to convective heat transfer from the hot flow to the test section walls. The turbulence intensity (TI) profile was reasonably flat, with an average of 13% and slightly elevated intensity in the lower half of the span. The integral length scale (defined from integral of autocorrelation; see Appendix D) was on average 0.99 *C*_{x}, with a value of zero enforced at the walls of the patched profile to ensure suitability as a CFD inlet condition.

Now consider the CFD-predicted profiles at plane 1. The total pressure profile is similar in shape and magnitude; the temperature profile is flat because the walls were modeled as adiabatic; both turbulence parameters are reasonably well matched. Overall, this suggests that the patched experimentally-measured boundary conditions are appropriate for use as inlet conditions in CFD.

## 8 Numerical Methods

Simulations of the high-pressure NGVs were performed for eight turbulence models, using the fully-resolved vane geometry, operating at engine-matched conditions of Mach and Reynolds numbers. We now briefly consider the computational domain, grid, model settings, and boundary conditions.

### 8.1 Domain and Grid.

The computational domain and grid are shown in Fig. 5. The mainstream inlet is located at plane 1, and the outlet is located at the end of the exit duct (see test section layout in Fig. 3). The film cooling holes, internal cooling passages, and trailing-edge slot were fully-resolved, with coolant introduced from inlets inside the respective plena at the hub and casing.

The 312 million-cell computational grid is hex-dominant, with refinement near walls, in proximity to small features, and inside the exit duct. A 23-cell prism layer was used for solid boundaries. A grid sensitivity study was performed for the SST2003 model (see Ref. [32]). This was considered to be the most pessimistic, because the other models are more diffusive by nature. For the SST2003 simulation, average *y*^{+} values on the pressure side, suction side, endwalls, and vane internals were 0.7, 1.3, 0.6, and 0.8, respectively.

### 8.2 Model Settings and Discretization.

Simulations of the NGVs were performed in ansys cfx using the SST2003, SKO, BSL, SST1994, SSTPSR, SSTDR, BSL-ARSM, and SA turbulence models. The working fluid was dry air, which was modeled as an ideal gas with constant transport properties. The compressible energy equation (including viscous work) was used. The turbulent Prandtl number was decreased from its default value of 0.90 (calibrated for wall heat-transfer) to a value of 0.50 (calibrated for thermal mixing layers; see Ref. [39]). Advective fluxes were discretized using a blended first and second-order upwind-biased scheme, optimized to locally ensure boundedness while minimizing diffusiveness. Diffusion and pressure terms were evaluated using finite element shape functions. Model-specific settings for the production limiter and realizability modification/limiter are summarized in Table 1.

### 8.3 Boundary Conditions.

Inlet conditions of total temperature, total pressure, length scale and TI were specified as the experimentally-measured and patched profiles at plane 1 (see Fig. 4). In ansys cfx, the length scale is defined as $\u2113=k/\beta *\omega $ (see Ref. [22]). Good collapse of the experimentally-measured integral length scales at inlet using this definition was demonstrated in Fig. 4 (see also the discussion about length scale definitions in Appendix D). For the SA model, the experimental length scales and turbulence intensities were translated into eddy viscosity using Eq. (3) (with *k* and *ℓ* calculated using the same method as for two-equation models). Hub and casing coolant inlets were mass flow inlets with constant total temperature, a TI of 5% and a viscosity ratio of 10. The domain outlet was a radial equilibrium condition, with an identical average pressure for all models. All remaining boundaries are walls, which were modeled as adiabatic and hydraulically smooth.

## 9 Performance Metrics

In this section, metrics of aerodynamic efficiency, and wake spreading/decay rates are summarized.

*p*

_{02}and

*p*

_{2}are the total and static pressure downstream of the NGVs, respectively, which vary as a function of

*r*(radius), and

*θ*(circumferential coordinate).

*p*

_{01}is the upstream mass-flux-averaged total pressure, and exponent

*χ*is equal to (

*γ*− 1)/

*γ*, where

*γ*is the ratio of specific heats (assumed constant and equal to 1.40). The local KE loss coefficient,

*ζ*

_{1}, represents the local KE deficit of the real flow relative to an isentropic expansion from the upstream total pressure to the local downstream static pressure.

*T*

_{01}, and

*T*

_{0c}. $p02\xaf$ and $p2\xaf$ are, respectively, the mass-flux-weighted total pressure, and the area-averaged static pressure downstream of the NGVs. For experiments, these averages are evaluated over two vane pitches, following the same pair of wakes in each measurement plane. Finally,

*p*

_{0c}is the mass-flow-weighted total pressure of the coolant feeds. The plane-averaged KE loss coefficient,

*ζ*

_{2}, represents the fractional KE deficit of the bulk flow, relative to mass-weighted isentropic expansions of mainstream and coolant streams to the area-averaged downstream static pressure.

The plane-averaged KE loss coefficient can be decomposed into profile-region and endwall-region KE loss coefficients. In the current study, the profile-region KE loss coefficient, *ζ*_{p}, is calculated using Eq. (26), with the mass-flux-weighted total pressure, $p02\xaf$, evaluated over the central 80% of span (i.e., excluding endwall boundary layers). The endwall-region KE loss coefficient, *ζ*_{e}, can then be evaluated as the difference between the plane-averaged KE loss coefficient and the profile-region KE loss coefficient (i.e., *ζ*_{e} = *ζ*_{2} − *ζ*_{p}).

*ψ*, which is defined by:

*T*

_{02}(

*r*,

*θ*) is the downstream field of total temperature, with

*T*

_{01}and

*T*

_{0c}the mass-flux-weighted total temperatures of the hot mainstream at inlet, and coolant streams at inlet, respectively. Due to spatial non-uniformity of total temperature at inlet, it is possible for

*ψ*to be locally negative.

We also consider the aerodynamic and thermal wakes in terms of their spreading and decay rates, so as to compare to predicted mixing rates using various turbulence models for engine-realistic inlet turbulence. Amend et al. [32] studied spreading and decay rates in terms of the peak intensity and the full width at half maximum (FWHM) of *ζ*_{1} and *ψ*. FWHM can be misleading for measuring spreading, however, because it does not account for wake *shape*. As an example, a wake profile with high excess kurtosis (i.e., *tailedness*) can still have low FWHM, because the distribution at levels below half the maximum does not affect the FWHM value. The degree of spreading can also be quantified in terms of the standard deviation of the wake distribution (i.e., assuming the wake profile is a continuous probability density function for *θ*, normalized such that integrated/total probability is unity). In the current study, both the standard deviation of the wake distribution and its peak intensity were calculated based on the radially-averaged wake profiles of *ζ*_{1} (10–90% span) and *ψ* (20–80% span). For annular cascades, the wake generally deforms due to radial variations in the circumferential progression rate of the wake. In a direct radial averaging approach, this can result in smearing of the peaks. For this reason, the wakes were *restacked* (peaks were circumferentially aligned), prior to averaging radially (averaging performed in the transformed coordinate system).

## 10 Results and Analysis

Simulations using eight turbulence models (see Table 1) are compared to experimental data from fully-cooled NGVs at engine-matched conditions of Mach number, Reynolds number, and turbulence. The purpose is to compare the performance of turbulence models in terms of total KE loss and wake mixing characteristics. The aim is to help inform turbulence model selection for aerodynamic simulations in conditions of high TI and large length scale (e.g., HP NGVs).

### 10.1 Kinetic Energy Loss and Cold Gas Mass Fraction at Downstream Measurement Plane.

In this section we consider the aerodynamic and thermal wakes at the downstream measurement plane (plane 4 in Fig. 3) 0.75 *C*_{x} downstream of the midspan trailing edge. In Fig. 6 we compare experimental data with CFD predictions using eight different turbulence models. The downstream plane was chosen because it emphasizes turbulence model differences the most (greater development distance). The plane is close to the HP rotor leading edge plane. The local kinetic energy loss coefficient, *ζ*_{1}(*r*, *θ*), is shown in the left column of Fig. 6, and the cold gas mass fraction, *ψ*(*r*, *θ*), is shown in the right column of Fig. 6. Radially, the data spans from hub to casing, covering two vane pitches in the circumferential direction. The data are viewed from downstream, and in normalized cylindrical coordinates, where *r*_{norm} is the span-normalized radius (datumed to the hub) and *θ*_{norm} is the vane-pitch-normalized circumferential location. For each turbulence model, comparisons are drawn to the corresponding experimental data, and the SST2003 reference model.

We now consider each row of Fig. 6 in turn, examining the wake shapes, peak intensity, and the width of both the aerodynamic and thermal wakes.

*Experimental data; row 1.*First consider the contours of KE loss (aerodynamic wake), shown on the left-hand-side. The bowed wakes (regions of total pressure deficit) can be seen at approximately*θ*_{norm}= 0 and 1. Near*r*_{norm}= 0 and 1, loss increases steeply on account of the endwall boundary layers. Loss cores coalescing with the wake around 15% and 75% span correspond to secondary flows. The hub secondary flow has greater intensity on account of a lower static pressure and correspondingly higher Mach number. Now consider the contours of cold gas mass fraction (thermal wake), which are shown on the right-hand-side. As the thermal wake is dominated by (relatively spanwise uniform) coolant ejection from films on the vane surface and from the trailing edge slot, the thermal wake shape is similar to the aerodynamic wake shape. Thermal secondary flows appear as accumulations of cold fluid (high*ψ*) at*θ*_{norm}= 0 and 1, around 15% and 90% span. These coalesce with the near-wall thermal deficit (increase of*ψ*) of the platform cooling flows (see cooling configuration in Fig. 3).*SST2003 model; row 2.*Predictions using the SST2003 reference model give rise to aerodynamic and thermal wake*shapes*which match the experimental data reasonably well, although the CFD-predicted wakes are significantly undermixed compared to experiments (for detailed analysis and discussion see Ref. [32]). Likewise, features of the aerodynamic and thermal*endwall flow*are adequately captured but appear undermixed compared to experimental data. Undermixing is related to the shear stress limiter (see Eq. (16)), which is active throughout the vane passage (see Fig. 2) and reduces eddy viscosity significantly below the realizability limit. This inhibits turbulence production and stresses (via eddy viscosity), consequently reducing mixing of the aerodynamic boundary layer and thermal mixing layer with the freestream.*SKO and BSL models; rows 3 and 4.*These models are mathematically very similar to each other; the only distinction is that the*k–ω*cross-diffusion term is omitted in SKO, and that closure coefficients are not blended in the SKO model. For this reason, and because predictions are almost identical, we treat the results together. We recall that SKO and BSL differ from SST2003 in that they do not incorporate a shear stress/realizability limiter. The effect is twofold because—without a realizability limiter—the incoming eddy viscosity is fully advected (i.e., not suppressed algebraically), which has the knock-on effect of amplifying TKE production around the SS peak $(P\u2248\mu tS2)$. The additional TKE enhances diffusion of the airfoil boundary layers, which leads to aerodynamic wakes that are broader and more radially uniform than using the SST2003 model. Peak KE loss in the profile-region (central 80% of span) is not significantly reduced, however, because the greater turbulence also increases KE dissipation at the wall (see flat plate data). In terms of the aerodynamic wake shape, the casing secondary flow loss cores in SKO and BSL predictions migrate closer to the casing, which leads to reduced likeness with experiments. The near-endwall KE loss distributions are very similar to that of SST2003, and experiments. Now consider the predicted thermal wake for the SKO and BSL models. The thermal wake is significantly more mixed-out than for the SST2003 model, with lower peak*ψ*in the profile region of the wake (central 80% of span). Secondary flows marked by cold gas (roll-up of the platform cooling flow) are also more mixed-out than the SST2003 model, such that their shape is a better match with experimental data. In between the wakes, the predicted penetration depth of endwall coolant into the freestream is significantly greater than the SST2003 prediction and slightly greater than experiments (i.e., overmixed).*SST1994 model; row 5.*This model uses a vorticity-based shear stress limiter (see Eq. (21)) which is only active in boundary layers. In other words, the model is not realizable in the freestream, which leads to over-production of TKE and excessive spreading of both aerodynamic and thermal wakes, compared to experiments. The increase in spreading (when compared to SKO and BSL) can be explained by the relaxed production limiter clip factor of 20 (used for consistency with Ref. [12]) versus the default value, $Clim=10$ (used for SST2003, SKO, and BSL). In combination with a non-realizable stress–strain relation, $Clim=20$ allows for greater freestream production of TKE around the SS peak than any other model. The additional TKE increases eddy viscosity and KE loss. Predictions from the SST1994 model feature a speckled eddy viscosity distribution around the SS peak (not shown), which occurs due to numerical error in the strain rate tensor and a feedback loop with the realizability limiter (see also Ref. [26]).*SSTPSR and SSTDR models; rows 6 and 7.*These models are based on the same transport equations as the SST model but use modified definitions for the realizability/stress limiter. SSTPSR enforces*exact*stress realizability (see Eq. (22)) by calculating principal strains directly, and SSTDR*estimates*principal strains using incompressible arguments (see Eq. (23)). Both models were tested*without*a TKE production limiter $(Clim=999)$. This corresponds to the condition of maximum TKE production, given the requirement that principal Reynolds stresses cannot be negative. This contrasts with the SST2003 model, which is also realizable, but enforces a much stricter limit on TKE production and turbulent stresses. The SSTPSR and SSTDR models are discussed collectively because they are conceptually alike (exact stress realizability), and because their aerodynamic and thermal wake profiles at 0.75*C*_{x}are very similar. Compared to the SST2003 predictions, the predicted aerodynamic and thermal wakes were slightly broader, with enhanced radial mixing. Peak intensity of KE loss was very similar to the SST2003 model. The predicted thermal wake shapes and secondary flows marked by cold gas (platform cooling roll-up) were considerably more mixed-out than using the SST2003 model, but slightly less mixed than in experiments. Overall the results for SSTPSR and SSTDR were similar; users may prefer to use SSTDR because of the relative ease of implementation and reduced computational cost.*BSL-ARSM model; row 8.*In this model the Boussinesq approximation is replaced with a non-linear stress–strain relation (i.e., Eq. (24)) with inbuilt realizability. For this reason, BSL-ARSM can be used without a TKE production limiter (i.e., $Clim=999$). The predictions using the BSL-ARSM model show significantly wider and less intense wakes than the SST2003 model. Overall, the BSL-ARSM model appears to give the best match with experiments among the models tested. The thermal field matches experiments well both in terms of profile-region coolant spreading and the shape of secondary flows marked by cold gas, but is slightly overmixed near the endwalls.*Spalart–Allmaras (SA) model; row 9.*This model differs from the other models in that it is based on a single transport equation for an eddy viscosity-like variable (instead of separate transport equations for*k*and*ω*). The SA model predicts an aerodynamic wake that is similar to that from the SSTPSR and SSTDR models, both in terms of shape and peak intensity. The corresponding thermal wake is slightly less mixed-out than for the SSTPSR and SSTDR models, especially in terms of the endwall cooling flow. Overall, SA may be an attractive choice for aerodynamic calculations because its predictions are similar to the SSTDR/SSTPSR models but the model is easier to recalibrate than its two-equation counterparts. Notably, by nature of the construction of the Boussinesq approximation (see Eq. (1))*without*the TKE term, the SA model*always*produces non-realizable stresses. Nevertheless, diffusivity does not*explode*(like in SST1994, for example) because the positive feedback loop through the TKE production term is not present. Despite non-physical Reynolds stresses, turbulent momentum diffusion is adequately captured because the effect depends on*spatial gradients*of the stresses, not their absolute values. Non-realizable stresses*may*affect diffusion of the thermal field, however, because the turbulent viscous work term in the RANS energy equation depends on the*absolute*values of the turbulent stresses. This might also explain why the thermal field of the SA model appears slightly less mixed than for the SSTPSR and SSTDR models, despite having a very similar aerodynamic field.

### 10.2 Plane-Averaged, Profile-Region, and Endwall-Region Loss.

High freestream TI has been shown to increase the KE loss in both the flat plate environment [6,29,30] and NGV [6,41,42]. We therefore consider the impact of turbulence model on integrated KE loss. The trends of plane-averaged KE loss coefficient, *ζ*_{2}, profile-region KE loss coefficient (central 80% span), *ζ*_{p}, and endwall-region KE loss coefficient (top and bottom 10% span), *ζ*_{e}, are shown in Figs. 7(a)–7(c), respectively, as a function of normalized axial distance, *x*/*C*_{x}. CFD data were interrogated at axial intervals of approximately 0.4 mm. Experimental data are shown for planes 2, 3, and 4 (0.25, 0.50, and 0.75 *C*_{x}, respectively); the same pair of wakes was considered at each plane (i.e., the same region of flow). Experimental bias uncertainty for *ζ* was approximately ±0.0016 (absolute; 95% confidence interval). The bias uncertainty was calculated from manufacturer-quoted instrument accuracy for combined linearity, hysteresis, and repeatability error.

First, consider the trends of plane-averaged KE loss coefficient (*ζ*_{2}) in Fig. 7(a). The CFD-predicted trends are relatively similar to each other in terms of *shape* but are offset in absolute value of the coefficient. We consider the CFD predictions in ascending order of *ζ*_{2}. Plane-averaged KE loss was lowest for the SST2003 model, with an 11.9% underprediction relative to experiments (on average across planes 2–4). Losses were very similar between SSTDR (1.4% lower than experiment), BSL-ARSM (1.2% lower than experiments), SSTPSR (1.5% greater than experiment), and SA (3.3% greater than experiment). The predictions of these four models are within the estimated experimental bias uncertainty range at all three planes. The remaining three models (BSL; SKO; SST1994) significantly overpredict plane-averaged KE loss: by 21.1% for the BSL model, by 23.3% for the SKO model, and by 56.1% for the SST1994 model. This result is consistent with the expectation that turbulence models with less restrictive turbulence production terms (realizability limiter and production clip factor) will predict greater KE losses. Indeed, only the non-realizable models (BSL, SKO, and SST1994) *overpredicted* the experimentally-measured plane-averaged KE loss. Conversely, in the presence of large length scales, SST2003 enforces the most restrictive condition of realizability in the freestream and was the only model to *underpredict* plane-averaged KE loss. The agreement of the realizable models (SSTDR, SSTPSR, and BSL-ARSM) with the experimental data is—of course—a good outcome, but the performance would be expected to worsen if surface roughness effects were modeled. The effect was investigated by performing an SST2003 simulation, using a sandgrain roughness of 2.5 × 10^{−4} × *C*_{x} for all wetted NGV surfaces (based on representative profilometry data of a similar new vane; see Ref. [43]). The roughened walls gave rise to a 15.0% increase in plane-averaged KE loss, compared to the SST2003 simulation with smooth walls. Another effect which is known to affect KE loss is the state of the boundary layer (turbulent versus laminar). The turbulence models assume fully turbulent boundary layers by design (high dissipation coefficient) across the entire vane surface. In the real flow—although the leading edge showerhead cooling flows are expected to cause immediate transition—there is likely local relaminarization on the SS [44], which has an associated dissipation coefficient lower than that for an equivalent turbulent flow (see dissipation coefficient correlations [45]). If this relaminarization could be accurately modeled in CFD, we would expect this to reduce the predicted profile-region KE loss. That is, the effect of adding this in the modeling would act in a contrary direction to the effect of adding surface roughness.

Now consider the trends of profile-region KE loss coefficient, *ζ*_{p}. These are shown in Fig. 7(b). The general trend of both the experimental measurements and CFD-predicted data for profile-region KE loss coefficient is a gently asymptotic form which remains approximately constant for *x*/*C*_{x} > 0.50. Asymptotic behavior indicates that secondary kinetic energy has fully mixed out (see also discussion in Refs. [32,33]), meaning that the approximately constant plateau value can be viewed as the mixed-out profile-region KE loss. As with the plane-averaged KE loss coefficient data (Fig. 7(a)) the CFD trends for all eight turbulence models are similar in *shape* to each other but offset in absolute value of the coefficient.

Consider now the CFD predictions of profile-region KE loss in ascending order of profile-region KE loss coefficient: SST2003 predicts a trend of *ζ*_{p} that is (average of planes 2–4) 11.8% lower than experimentally-measured values; SSTPSR, BSL-ARSM, and SSTDR have average overpredictions of 1.5%, 1.7%, and 1.7%, respectively; SA has greater profile-region loss, with an average overprediction of 8.4%; for the non-realizable two-equation models, overpredictions of *ζ*_{p} were 29.4% (BSL), 31.4% (SKO), and 73.0% (SST1994). Once again, the models with *exact* realizability matched experimental data best, with significant overpredictions for the non-realizable models, and a slight underprediction using SST2003. In superposition with surface roughness effects, we expect profile-region loss to increase for all models. This could well worsen the good agreement between experimental data and the SSTDR, SSTPSR, and BSL ARSM models. Indeed, the profile-region loss coefficient increased by 15.1% in the SST2003 simulation including surface roughness, which corresponds to a 1.6% overprediction, compared to experiments. Similarly, we might expect a slight reduction of the CFD-predicted profile KE loss, if SS boundary layer relaminarization had been modeled (boundary layers are currently modeled as fully turbulent).

Finally, consider the trends of endwall-region KE loss coefficient, *ζ*_{e}. These are shown in Fig. 7(c) (note scaled *y*-axis). All CFD-predicted trends of *ζ*_{e} essentially overlay. Considering that *ζ*_{2} = *ζ*_{p} + *ζ*_{e}, this indicates that the variations of plane-averaged (i.e., total) KE loss are primarily driven by an increase of profile-region KE loss. This suggests that models can be tuned to match profile loss in isolation. A possible explanation for the low sensitivity of endwall-region KE loss to the turbulence model is that the CFD inlet boundary condition (experimentally-measured at plane 1; see Fig. 3) has much lower eddy viscosity near the endwalls, than in the freestream. Recalling that TKE production is proportional to eddy viscosity, this reduces the influence of the production limiter and/or realizability modification in this region. Endwall-region KE loss predictions from all models fall inside the experimental bias uncertainty range at planes 3 and 4, with a mis-match at plane 2 (0.25 *C*_{x}). We note that some aspects of the experimental measurement uncertainty are hard to quantify, particularly in relation to measuring very thin boundary layers (exacerbated at plane 2) with a probe of finite tip diameter.

For completeness, we note the presence of small variations in the BSL-ARSM trends of profile-region and plane-averaged KE loss. These occur due to non-monotonic decrease of mass-flux-averaged total pressure. This effect was investigated using a simplified test case. The behavior was insensitive to grid coarsening, time-step refinement, and floating-point precision, but appeared to be amplified/damped depending on the inlet turbulence. This suggests that the phenomenon is model-related. The BSL-ARSM model is a fairly recent development, and this may be the result of some model-instability at extreme boundary conditions of turbulence.

### 10.3 Aerodynamic and Thermal Wake Shapes.

Undermixing of NGV wakes is a frequently-observed phenomenon affecting CFD predictions using *steady* RANS models. Lower-than-experimental spreading-rate and peak-decay-rate occur because trailing edge *unsteady* vortex shedding and the associated mixing cannot be captured by *steady* RANS. The studies of Burdett and Povey [33], as well as Kopriva et al. [46] demonstrated that using *unsteady* RANS increases the CFD-predicted spreading/decay rates slightly, when compared to *steady* RANS, resulting in marginally better collapse with experimental data.

An additional mechanism which affects the wake width downstream of the NGVs is SS boundary layer growth. Using a linear turbine cascade, Folk et al. [6] investigated the effects of high inlet TI, reporting a 37% increase of kinetic energy loss, and a 28% increase of wake width (compared to low TI inlet conditions). Blade-surface boundary layer measurements using hotwire probes showed that the increase in kinetic energy loss and wake width for high TI conditions were primarily driven by enhanced mixing between the freestream and the SS boundary layer. The boundary layer diffusion effect is separate from (but possibly interacting with) the unsteady trailing edge effects but does not appear to be inherently unsteady in itself. It should be possible to predict the non-coupled component of this effect with *steady* RANS.

Before examining the spreading and decay rates of the aerodynamic and thermal wakes, we analyze the *shapes* of the wakes. Radially-averaged profiles of the *restacked* aerodynamic and thermal wakes at 0.75 *C*_{x} (plane 4; see Fig. 3) are shown in Fig. 8. The restacked aerodynamic and thermal wakes were averaged over span ranges of 10–90%, and 20–80% respectively. The averaging window for the thermal wake was slightly narrower than for the aerodynamic wake to reduce the influence of endwall coolant penetration, which was sensitive to the choice of turbulence model (see Fig. 6). Experimental data are marked using circles. The bias uncertainty in *ζ*_{1} was estimated to be ±0.0016 (absolute; 95% confidence interval). The bias uncertainty in *ψ* was estimated to be ±0.04 (absolute; 95% confidence interval; propagated from K-type thermocouple bias uncertainty [47]). CFD data are shown for the eight turbulence models in Table 1.

First consider the circumferential aerodynamic wake profiles (restacked and radially-averaged) of local KE loss coefficient, *ζ*_{1}, shown in Fig. 8(a). The experimentally-measured wake profile is approximately symmetric, with peak value *ζ*_{1} = 0.103. All CFD-predicted wake profiles have slight asymmetry with greater weight toward the SS, and with significantly greater peak loss than the experiment (between 20.0% and 49.0% greater). CFD predictions using the SST2003, SKO, BSL, SSTPSR, SSTDR, and SA models all have similar peak loss, which is between 40.7% and 49.0% higher than the experimentally-measured value. This suggests that the realizability modification differences between these models do not have a strong influence on peak loss. SST1994 and BSL-ARSM have closer peak loss to the experimental value, with respective over-predictions of 23.5% and 20.0%. The turbulence models differentiate themselves from each other mostly in terms of their KE loss distribution on the SS of the wake (see also contours of Fig. 6). This is consistent with Folk et al. [6], who found that freestream turbulence primarily affects the SS boundary layer, and the corresponding SS of the wake. The pressure side (PS) wake profile was less affected by the choice of turbulence model (except for the SST1994 model). Over-spreading for the SST1994 model (on both PS and SS) is related to the relaxed TKE production limiter clip factor $(Clim=20)$, and the lack of a freestream realizability limiter, which give rise to excessive spreading rates and KE losses.

Now consider the circumferential profiles of the restacked and radially-averaged thermal wake, which are shown in terms of the cold gas mass fraction, *ψ*, in Fig. 8(b). The experimentally-measured thermal wake is broadly symmetric about the peak value of cold gas mass fraction (0.211) and similar in shape to the corresponding aerodynamic wake. The CFD-predicted thermal wakes are also generally symmetric, but the wake widths are sensitive to the turbulence model. The thermal wake peak intensity was greatest for SST2003 (0.352 cold gas mass fraction; +66.6% compared to experiments). The remaining models, by descending order of peak cold gas mass fraction, were SSTDR (0.264; +24.9%), SA (0.256; +21.1%), SSTPSR (0.243; +15.1%), BSL (0.204; −3.4%), BSL-ARSM (0.200; −5.4%), SKO (0.197; −6.9%), and SST1994 (0.138; −34.6%). Perhaps surprisingly, the model ordering by peak cold gas mass fraction does not correlate one-to-one with the ordering by peak KE loss. The reason is that cold gas mass fraction is a conserved quantity, such that greater diffusivity (i.e., eddy viscosity) will always lead to additional spreading *and* decay. For the aerodynamic wake, greater model diffusivity increases wake width, but is also accompanied by a rise in the profile-region KE loss coefficient (due to additional turbulent shear stress near the walls). As such, peak KE loss coefficient values do not necessarily decrease when the wake width increases. This implies that the thermal field is easier to optimize, because the cost function is simpler.

### 10.4 Aerodynamic Wake Spreading and Decay.

We now consider the spreading and decay trends of the aerodynamic wake for experiments and eight turbulence models. Spreading and decay rates merit attention because the dynamic loading of the HP rotor blade depends on the intensity of the NGV wake. Aerodynamic wake width and intensity are studied using radially-averaged profiles (10–90% span) of the restacked wake. A wide span range was chosen to reduce bias from local gaps or hotspots in the wake, and to improve precision uncertainty of the experimental data. We consider the axial range 0.0 *C*_{x} and 1.0 *C*_{x}. We quantify wake width in terms of the pitch-normalized standard deviation of the wake distribution, *σ*/*s*, and we quantify wake intensity in terms of the peak value of (radially-averaged) local KE loss, $\zeta 1,max$. The pitch-normalized standard deviation *σ*/*s* has been multiplied by $12$ (the standard deviation of a uniform distribution over the interval 0 to 1 in pitch is $1/12$) such that the “width parameter” is unity for uniformly distributed KE loss.

First consider the trends of pitch-normalized and rescaled aerodynamic wake width (i.e., standard deviation of the wake distribution), *σ*/*s*, which are shown in Fig. 9(a). The experimental trend is for an approximately linear increase in wake width with axial distance. The CFD trends are broadly similar, but with periodic inflections, owing to varying wake whirl angle, which affects the wake width when evaluated on an axial plane (see discussion in Ref. [32]). We now consider the CFD models in ascending order of average wake width.

SST2003 predicted the narrowest wake with a width 55.0% lower than experiments (average across planes 2–4). SSTDR, SSTPSR, and SA closed the gap with experiments slightly, with respective underpredictions of 42.1%, 37.5%, and 32.2%. Wake width for BSL-ARSM was on average 30.7% lower than for experiments. SKO and BSL had respective underpredictions of 5.0% and 4.0%, and were the best match with experiments. As illustrated by the circumferential profiles of Fig. 8, the increase in wake width (compared to SST2003) is biased toward the SS, with relatively little change in the PS wake. The point is that the good agreement of both SKO and BSL model predictions with experimental data occurs due to excessive spreading of the SS boundary layer, which compensates for the underprediction of PS wake spreading. The SS spreading is precipitated by an overprediction of TKE production in the vane passage (both models are non-realizable). For both models, the excessive width of the SS-wake is accompanied by an approximately 30% overprediction of profile-region KE loss. SST1994 was the only model to overpredict aerodynamic wake width (+30.3%). Excessive spreading occurs because the model is non-realizable (i.e., susceptible to turbulence overproduction) and has a relaxed clip factor $(Clim=20)$. This was accompanied by a 73.0% overprediction of profile-region KE loss. The takeaway is that—for high TI NGV flows—turbulence model choice can significantly affect predictions of aerodynamic wake width through the degree of boundary layer spreading on the SS of the airfoil. SS boundary layer spreading correlates with increased profile-region KE loss, such that models that more accurately predict wake width generally overpredict loss.

Let us now consider the *growth rate* of the aerodynamic wake, which we define as the increase in width between 0.25 *C*_{x} and 0.75 *C*_{x}. Except for SST1994 and BSL-ARSM, CFD-predicted wake growth rates are fairly similar to each other and fall in the range between 38% and 61% of the experimentally-measured growth rate. Lower-than-experimental growth rates in CFD may partly be explained by the absence of unsteady vortex shedding effects, which were not modeled. It is therefore especially surprising that the SST1994 and BSL-ARSM models predicted growth rates that were in relatively close agreement with experiments. We consider the cause of this behavior now. For SST1994, the wake growth rate was 91.5% of the experimentally-measured values. Compared to the other *linear* eddy viscosity models, the spreading rate in SST1994 is amplified by the abnormally high TKE in the *freestream*, which originates from the region of TKE over-production around the SS-peak. The point is that while this may *imitate* unsteady effects, the greater spreading rate is an artifact of excessive production of TKE in the vane passage. We recall that this also led to a 73.0% overprediction of the profile-region KE loss coefficient. For BSL-ARSM, the wake growth rate was 121% of experimentally-measured values. Higher spreading rates in BSL-ARSM are potentially helpful because they improve the collapse with the experimental data *without* compromising the model accuracy in terms of its integral KE loss predictions. The mechanism for accelerated wake spreading in BSL-ARSM is not currently understood, but is likely related to how the constitutive relation for turbulent stress (i.e., Eq. (24)) responds to free shear.

Now consider the axial evolution of peak aerodynamic wake intensity, $\zeta 1,max$, which is shown in Fig. 9(b). The experimental trend is quasi-exponential decay. The CFD-predicted trends are similar in shape, but with higher peak wake intensities. Periodic inflections in the CFD trends appear due to the varying streamwise pressure gradient, which increases/decreases KE loss production rates (through mixing) as the wake passes through the NGV downstream potential field (see discussion in Ref. [32]). As already discussed in the context of Fig. 8, peak aerodynamic wake intensity was over-predicted by all turbulence models, with most models falling in the range +30.0% to +41.0% (average across planes 2–4). The SST1994 and BSL-ARSM models are the exception, with respective wake depth overpredictions of only 14.7% and 19.2% (average across planes 2–4). Slightly greater peak wake intensity *decay rates* in SST1994 and BSL-ARSM occur for the same reasons as the elevated *growth rates*.

In summary we conclude that the generic eddy viscosity models are unlikely to yield accurate predictions of both NGV integral KE loss *and* aerodynamic wake development *simultaneously* because the free shear response of the model drives both the production of TKE in the vane passage and the mixing of the wake. Therefore, any closure coefficient adjustments designed to increase free shear diffusivity in the wake (to *imitate* unsteady phenomena) are likely to also lead to an overproduction of TKE around the SS-peak, ultimately resulting in an overprediction of integral KE loss and increased spreading of the SS boundary layer. A solution might be separate zonal treatment of the wake region and the vane passage region, as has been explored by Pacciani et al. [48] and Marioni et al. [49]. In this type of approach, wake mixing can be tuned independently from KE dissipation and SS boundary layer spreading. The disadvantage is that such models lack generality, because they require *a priori* knowledge of the flow.

### 10.5 Thermal Wake Spreading and Decay.

We now consider the spreading and decay trends of the thermal wake for experimental measurements and predictions from eight turbulence models. Thermal wake width and intensity are presented in terms of cold gas mass fraction, *ψ*. Wake width and peak intensity are calculated from radially-averaged profiles of the restacked wake (20–80% span). Compared to the analysis of the *aerodynamic* wake, the radial averaging window was reduced to 20–80% span to mitigate the influence of endwall coolant spreading into the 10–90% span range for *some* turbulence models (see Fig. 6). We consider the axial range between 0.0 *C*_{x} and 1.0 *C*_{x}. Wake width was quantified using the pitch-normalized standard deviation of the thermal wake distribution (scaled by $12$), and wake depth is defined as the peak cold gas mass fraction, $\psi max$.

First consider the trends of pitch-normalized and rescaled thermal wake width, *σ*/*s*, which are shown in Fig. 10(a) as a function of the normalized axial distance, *x*/*C*_{x} (datumed to trailing edge). In similarity to the *aerodynamic* wake, the experimental *thermal* wake width increases approximately linearly with axial distance. CFD trends of the thermal wake width are approximately correlated with corresponding CFD predictions of aerodynamic wake width. This is expected, because the turbulent enthalpy flux is modeled using eddy diffusivity (turbulent Prandtl number of 0.50), such that the diffusivities of enthalpy and momentum are proportional. Accordingly, SST2003 predicted by far the narrowest thermal wake, with an average width 45.9% less than the experimental measurements (average across planes 2–4). Thermal wake width was also underpredicted by the SA (−32.3%), SSTDR (−25.6%), and SSTPSR (−17.5%) models. The gradual increase in average thermal wake width from model-to-model is related to relaxation of the stress limiting conditions and TKE production limiters, which gives rise to elevated freestream turbulence. The BSL, BSL-ARSM, SKO, and SST1994 models overpredicted the thermal wake width by +7.2%, +9.3%, +19.4%, and +46.1%, respectively. Notably, the proportion of models overpredicting *thermal* wake width (four out of eight models) was greater than the number overpredicting *aerodynamic* wake width (one out of eight models). The reason is that the thermal mixing layer (i.e., film coolant on the airfoil) is generally thicker than the attendant aerodynamic boundary layer, such that coolant is more exposed to the freestream turbulence and also more sensitive to its modeling. This may also explain why the CFD-predicted thermal wakes are broadly symmetric, whereas the corresponding aerodynamic wakes are weighted toward the SS.

Now consider the axial trends of peak cold gas mass fraction, $\psi max$, which are shown in Fig. 10(b). Because cold gas mass fraction is a conserved quantity, we expect the trends of peak thermal wake intensity to decay with the inverse of wake width, such that *width times intensity* is approximately constant. Indeed, both the experimentally-measured and CFD-predicted trends of $\psi max$ are related approximately inversely to the normalized axial distance. Overpredictions of $\psi max$ (compared to experiments; average across planes 2–4) occurred for the SST2003 (+62.5%), SSTDR (+21.9%), SA (+20.0%), and SSTPSR (+13.5%) models. The ordering is slightly different to that of the thermal wake *width*, due to the compounding effect of the wake *shape* (i.e., width times intensity is only *approximately* constant). Although it is a useful representation, this points to a limitation of quantifying wake width using a single scalar metric. Next, consider the BSL-ARSM, SKO, and BSL models, which predicted thermal wake intensity exceptionally well, with respective differences from experiments of only +0.5%, −2.2%, and −2.5%. Unfortunately, this good agreement is likely the result of compensating error, whereby the absence of unsteady vortex shedding effects (not modeled) is offset by excessive turbulence model diffusivity in the vane passage. This can be verified by inspecting the thickness of the endwall coolant layers for these models (see Fig. 6), which over-penetrate the freestream when compared to experimental measurements. If we take the—probably reasonable—assumption that unsteady effects are negligible for the real (i.e., experimental) endwall flow, this indicates that thermal field diffusivity for otherwise steady turbulent phenomena is overpredicted. In applications where the accurate prediction of endwall coolant spreading is a priority, the SSTPSR and SSTDR models may offer a better compromise (see Fig. 6).

## 11 Conclusion

The performance of eight eddy viscosity turbulence models was investigated for a high-pressure NGV flow with engine-realistic turbulence (large length scale, high intensity). The models considered were SST2003 [14], SKO [13] (*k–ω* without cross-diffusion), BSL [12], SST1994 [12], SSTPSR (SST with exact principal stress realizability; see Ref. [35]), SSTDR (SST with Durbin realizability condition [28]), BSL-ARSM [16], and Spalart–Allmaras [17]. The purpose was to evaluate model accuracy for predicting KE losses, boundary layer spreading, and wake mixing in applications with high levels of freestream turbulence. Each model was benchmarked against the experimental flat plate data of Folk et al. [6], and the fully-cooled NGV flow of Amend et al. [32].

For the flat plate test, CFD predictions were compared to experimentally-measured dissipation coefficient data for inlet conditions of low TI (2.3%) and high TI (18.7%). To emulate turbulence conditions representative of a rich-burn combustor, Folk et al. [6] performed their flat plate experiments with exceptionally large freestream length scales (roughly 40 times the momentum thickness). For low TI, all turbulence models except SKO predicted average dissipation coefficients within 10% of the experimentally-measured value. Good agreement between CFD and experiments at low TI may have been expected on the basis that turbulence models are typically calibrated at low TI conditions, but we note that freestream length scales were considerably larger than most turbulence model calibration data. At high TI, all models correctly predicted an increase of the dissipation coefficient (compared to low TI), but systematically overpredicted the absolute values, when compared to experiments (maximum error of +34.2%; see Appendix A). We conclude that turbulence models tend to overpredict skin friction losses for high TI wall-bounded flows.

For the fully-cooled NGV test case, CFD-predicted trends of KE loss, and aerodynamic/thermal wake mixing were compared to experimental measurements for inlet conditions of large length scale, high intensity turbulence (*ℓ* ≈ *C*_{x}; $TI\u224813%$). Simulations using the *de facto* industry standard SST2003 model revealed that SST blending functions effectively malfunction for large freestream length scales due to the explicit appearance of the length scale in the blending function arguments. This resulted in the vane passage being effectively treated as a boundary layer both in terms of closure coefficient blending and stress limiter activation (see Figs. 11 and 2). This highlights the potential pitfalls of automatic zonal treatments. Remedial turbulence model modifications and alternatives were explored.

The first SST blending function (responsible for closure coefficient blending) was modified to activate in the viscous sublayer only. This led to better-behaved zonal treatment (see Fig. 11) but did not affect KE loss or wake development significantly. Considering the low sensitivity of results, the modification was not carried forward into the subsequent adaptations of the shear stress limiter, which we consider now.

Erroneous activation of the SST2003 shear stress limiter throughout the vane passage (see Fig. 2) is conceptually similar to the intended operation of a *realizability limiter*. While stress realizability (i.e., the condition that turbulent stresses are mathematically achievable) can generally be regarded as a favorable model attribute, the stress limiting condition imposed by SST2003 is almost twice as restrictive as the classical realizability limit of Durbin [28], and therefore not necessarily valid in the freestream. Artificial reduction of turbulence in the freestream is especially problematic at high TI conditions, because flat plate data shows that freestream turbulence has a significant effect on predictions of KE dissipation [6]. Turbulence-driven diffusion of the aerodynamic boundary layer, and thermal mixing layer *on the vane surface* is also expected to play a significant role in the eventual width of the corresponding wakes. To investigate the role of freestream turbulence, results from turbulence models with various types of stress/realizability limiters (including SST2003) were compared against experimental measurements of KE loss and aerodynamic/thermal wake development. Other common eddy viscosity models, such as SA, BSL, SKO, and BSL-ARSM were also benchmarked. For completeness, we note that the experimental validation data includes the additional effect of unsteady trailing-edge vortex shedding, which also enhances aerodynamic and thermal field diffusivity, but which is not captured by *steady* RANS simulations. If modeled using *unsteady* RANS, we would expect for these effects to increase aerodynamic/thermal wake width slightly while reducing peak intensity proportionately, *without* affecting integral KE loss (see Ref. [33]).

Integrated loss was studied in terms of the plane-averaged KE loss coefficient, which was decomposed into contributions from the endwall (0–10% and 90–100% span) and profile regions (central 80% span). Endwall-region KE loss was relatively insensitive to the turbulence model, which led us to study integrated loss by proxy of the profile-region KE loss coefficient alone. Model-to-model variations of profile-region KE loss were as high as 196%. The CFD-predicted profile loss correlated inversely with the restrictiveness of the production and/or realizability limiter. The reason is that the models largely differentiate themselves from each other in terms of the TKE produced around the SS peak, which drives subsequent SS boundary layer losses. Accordingly, SST2003 predicted the lowest KE loss, with an average deficit of 11.8%, when compared to experiments. Profile-region KE loss predictions from the realizable SSTDR, SSTPSR, and BSL-ARSM models were within 2% of the experimentally-measured values. We should not place too much emphasis on the good absolute accuracy, considering that the effects of surface roughness (which increases loss) and potential boundary layer relaminarization (which would reduce loss) were not modeled. Non-realizable two-equation models (SKO, BSL, and SST1994) should be avoided at high TI conditions, because they significantly overpredict profile-region KE loss (29–73%) on account of TKE *overproduction* around the SS peak. Their results also exhibit dependence on the arbitrary production limiter threshold, which is used as a safeguard against uncontrolled TKE production.

Aerodynamic wake development was studied in terms of the radially-averaged profiles (10–90% span) of the local KE loss coefficient. All turbulence models overpredicted aerodynamic *wake depth* (i.e., peak KE loss) by between 15% and 41%, compared to experiments, presumably due to the absence of unsteady vortex shedding effects in the steady RANS simulations. Predictions of aerodynamic *wake width* varied between 45% and 130% of the experimentally-measured average (across planes 2–4; see Fig. 3). Model-to-model differences in aerodynamic *wake width* were largely dominated by the degree of SS boundary layer diffusion, which was driven by levels of TKE production around the SS-peak. When benchmarking CFD predictions against experimental data taken *downstream of the vanes*, this mechanism of boundary layer spreading is easily conflated with unsteady vortex shedding effects, which were not modeled. Steady RANS simulations using the SKO and BSL models, for example, predicted similar aerodynamic wake width as experiments (underpredictions of 5%, and 4%, respectively). In these cases, the SS spreading effect compensated for the absent unsteady effect, but came at the expense of an overprediction of KE loss (approximately 30%), and greater wake asymmetry than for experimental data (see Fig. 8). Of course, it would be desirable—in principle—if eddy viscosity models could accurately mimic unsteady effects (considering the lower computational cost). The current results emphasize, however, that SS spreading and wake growth rates downstream of the trailing edge are inherently linked, because they both depend on the free shear response of the model. It is therefore not possible to independently enhance wake spreading rates, because the necessary model modifications would also trigger additional TKE production around the SS-peak, resulting in excessive KE dissipation and diffusivity throughout the vane passage. Simultaneous predictive accuracy for KE loss, SS boundary layer diffusion, and wake spreading would likely require the separate zonal treatment of the wake and the vane passage (see, e.g., Refs. [48,49]).

Thermal field development was investigated using radially-averaged profiles (20–80% span) of the cold gas mass fraction. Because cold gas mass fraction is a conserved quantity, thermal wake width and peak intensity were approximately inversely related (both as a function of axial distance and from model-to-model). This is in contrast to the aerodynamic wake (i.e., KE loss), for which the models with widened wakes also experienced a rise of integrated loss, causing peak KE loss decay trends to remain largely unchanged. Because turbulent enthalpy fluxes were modeled using eddy diffusivity (turbulent Prandtl number of 0.50), we expect thermal field mixing to scale with aerodynamic diffusivity. Indeed, the thermal field of the SST2003 simulation was the most undermixed, with a deficit in wake width of 45.9% and an excess in peak intensity of 62.5%, compared to experiments (average across planes 2–4). Predictions of peak cold gas mass fraction from the SA, SSTDR, and SSTPSR models closed the gap with experiments by approximately 70% of the initial difference. SKO, BSL, and BSL-ARSM achieved the best absolute agreement with experiments, although over-penetration of endwall coolant into the freestream (when compared to experiments) suggests that this only occurs because excessive diffusivity in the vane passage compensates for the fact that unsteady vortex shedding effects were not modeled. Where endwall coolant mixing is a priority, the SSTPSR or SSTDR models might be a good choice. Finally, we note that the turbulent Prandtl number of 0.5, used in the current simulations, is not suitable for wall heat transfer problems. Here it may be helpful to blend the turbulent Prandtl number (0.5 in the freestream; 0.9 near the wall) using the *modified* first blending function (see Appendix B).

Turbulence model sensitivity to a change of inlet TI was investigated in terms of CFD-predicted profile-region KE loss coefficient and aerodynamic wake width (see Appendix C). As expected, all models predicted an increase of profile-region KE loss and aerodynamic wake width for inlet TI increasing from 1.3% to 13% (constant length scale; *ℓ* ≈ *C*_{x}). Surprisingly, the profile-region KE loss coefficient and wake width rose by only 2.7% and 2.1% (respectively) using SST2003: a result which might be regarded as non-physical (too low) when compared to representative experimental studies. The realizable models achieved a better match with the reported sensitivity in the literature. Notably, for all *linear* eddy viscosity models, the increase of wake width (for TI increasing from 1.3% to 13%) was approximately proportional to the corresponding increase of profile KE loss. This reiterates the strong correlation between SS spreading and profile KE loss for this type of turbulence model.

An unresolved ambiguity in modeling large length scale turbulence is the translation of experimentally-measured *integral* length scales into *turbulence length scales* (to estimate turbulence dissipation rate). Especially at high TI, the chosen definition has a significant impact on predictions of SS spreading and KE loss. Considering that the *turbulence length scale* is largely a modeling concept, the proportionality constant between turbulence length scale and integral length scale might be viewed as an optimization parameter. Surprisingly, the values of the proportionality constant found in the literature vary by an order of magnitude, depending on the CFD package. We recommend explicitly reporting the method used to estimate turbulence dissipation rate from experimental measurements.

## Acknowledgment

The financial support of Innovate UK, the Aerospace Technology Institute, the Engineering and Physical Sciences Research Council, and Rolls-Royce plc are gratefully acknowledged, as is the technical support of Rolls-Royce plc with preceding programs from which we inherit data used in this paper. We also thank Folk, Miller, and Coull for providing additional unpublished information relating to Ref. [6].

## Conflict of Interest

The authors attest that all data for this study are included in the paper.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## Nomenclature

### Romans

*a*_{1}=Bradshaw constant (–)

*a*_{ij}=Reynolds stress anisotropy (–)

- arg
_{1/2}= blending function 1/2 argument (–)

*C*=compressible production coefficient (–)

*C*_{d}=dissipation coefficient (–)

*C*_{f}=skin friction coefficient (–)

- $Clim$ =
production limiter clip factor (–)

- $E$ =
turbulence dissipation rate (kg m

^{−1}s^{−3})*F*_{1/2}=blending function 1/2 of SST model (–)

*H*_{23}=momentum-to-energy thickness ratio (–)

*k*=turbulence kinetic energy (m

^{2}s^{−2})*ℓ*=turbulence/integral length scale (m)

- $m\u02d9$ =
mass flow rate (kg s

^{−1})*p*=pressure (Pa)

- $P$ =
shear production of turbulence (kg m

^{−1}s^{−3})- $Re\theta $ =
momentum thickness Reynolds number (–)

*r*=radius (m)

*s*=angular vane pitch (deg)

*S*=strain rate invariant (s

^{−1})*S*_{ij}=strain rate tensor (s

^{−1})*T*=temperature (static or total) (K)

*u*=fluctuating velocity component (m s

^{−1})*U*=mean-flow velocity component (m s

^{−1})*v*=turbulence velocity scale (m s

^{−1})*x*=Cartesian coordinate (m)

*y*=wall distance (m)

### Greeks

*α*=closure coefficient for production of

*ω*(–)*β*=closure coefficient for dissipation of

*ω*(–)- $\beta *$ =
closure coefficient for dissipation of

*k*(–)*γ*=ratio of specific heats (–)

*δ*_{ij}=Kronecker delta (–)

- $\u03f5$ =
turbulence dissipation rate (m

^{2}s^{−3})*ζ*_{1}=local KE loss coefficient (–)

*ζ*_{2}=plane-averaged KE loss coefficient (–)

*ζ*_{p}=profile-region KE loss coefficient (–)

*ζ*_{e}=endwall-region KE loss coefficient (–)

*θ*=circumferential coordinate (deg)

*λ*=strain rate tensor eigenvalue (s

^{−1})*μ*=dynamic viscosity (Pa s)

*μ*_{t}=eddy (turbulent) viscosity (Pa s)

*v*_{t}=kinematic eddy viscosity (m

^{2}s^{−1})*ρ*=volumetric mass density (kg m

^{−3})*σ*=standard deviation of circumferential wake profile when treated as a probability density function (deg)

*σ*_{k}=eddy diffusivity of

*k*(–)- $\sigma \omega $ =
eddy diffusivity of

*ω*(–)*χ*=isentropic expression (–)

*ψ*=cold gas mass fraction (–)

*ω*=specific turbulence dissipation rate (s

^{−1})*Ω*=vorticity magnitude (s

^{−1})*Ω*_{ij}=vorticity tensor (s

^{−1})

### Subscripts

## Appendix A: Flat Plate Test Data

RANS simulations of a flat plate were completed for conditions of low and high freestream TI at inlet. The purpose was to evaluate the accuracy of common turbulence models in a high-turbulence environment by comparing them with experimental validation data. The experimental data and test section geometry are from Folk et al. [6,31]. The simulations were completed on a nominally two-dimensional CFD grid, which was extruded by a single cell in the spanwise direction of the flat plate. The inlet was located 120 mm upstream of the elliptical leading edge (6:1 aspect ratio) of the flat plate, and was specified as a constant total gauge pressure of 340 Pa, a turbulence length scale of 40 mm $(\u2113=k/\beta *\omega )$, and a total temperature of 293.15 K. The low and high turbulence cases correspond to inlet TI levels of 2.3% and 18.7%, respectively. The outlet condition was an absolute static pressure of 1.0 atm (ambient pressure). Model-specific production and realizability limiter settings are summarized in Table 1. Dissipation coefficient data for low and high TI inlet conditions are summarized in Table 2.

Low turbulence 10^{3} × C_{d} (versus exp.) | High turbulence 10^{3} × C_{d} (versus exp.) | |
---|---|---|

Experiment (exp.) Folk et al. [6] | 1.601 | 2.024 |

SST2003 | 1.651 (+3.1%) | 2.416 (+19.4%) |

SSTF1V | 1.726 (+7.8%) | 2.614 (+29.2%) |

SKO | 1.891 (+18.1%) | 2.609 (+28.9%) |

BSL | 1.735 (+8.4%) | 2.620 (+29.4%) |

SST1994 | 1.699 (+6.1%) | 2.716 (+34.2%) |

SSTPSR | 1.743 (+8.9%) | 2.552 (+26.1%) |

SSTDR | 1.736 (+8.4%) | 2.521 (+24.6%) |

BSL-ARSM | 1.715 (+7.1%) | 2.691 (+33.0%) |

SA | 1.734 (+8.3%) | 2.051 (+1.3%) |

SKE (C_{lim} = 10) | 1.697 (+6.0%) | 2.631 (+30.0%) |

Low turbulence 10^{3} × C_{d} (versus exp.) | High turbulence 10^{3} × C_{d} (versus exp.) | |
---|---|---|

Experiment (exp.) Folk et al. [6] | 1.601 | 2.024 |

SST2003 | 1.651 (+3.1%) | 2.416 (+19.4%) |

SSTF1V | 1.726 (+7.8%) | 2.614 (+29.2%) |

SKO | 1.891 (+18.1%) | 2.609 (+28.9%) |

BSL | 1.735 (+8.4%) | 2.620 (+29.4%) |

SST1994 | 1.699 (+6.1%) | 2.716 (+34.2%) |

SSTPSR | 1.743 (+8.9%) | 2.552 (+26.1%) |

SSTDR | 1.736 (+8.4%) | 2.521 (+24.6%) |

BSL-ARSM | 1.715 (+7.1%) | 2.691 (+33.0%) |

SA | 1.734 (+8.3%) | 2.051 (+1.3%) |

SKE (C_{lim} = 10) | 1.697 (+6.0%) | 2.631 (+30.0%) |

Note: CFD data were resampled at experimental values of $Re\theta $.

At low TI, all models except SKO predict the dissipation coefficient accurately to within 10% of the experimental data of Folk et al. [6]. The significant overprediction of *C*_{d} in SKO (+18.1%) presumably arises due to the well-known sensitivity of SKO to freestream values of *ω*, which are unusually low in the current simulations (see discussion in Ref. [19]).

At high TI, the dissipation coefficient was systematically overpredicted, with a maximum error of +34.2% (for SST1994). The model-to-model variation was also somewhat greater at high TI, which presumably occurs because the models differentiate themselves from each other in how they handle extreme levels of turbulent stress (realizability and production limiters). For the flat plate test case, SA was the most accurate model.

## Appendix B: Modified Blending Function

For the NGV test case, the first blending function (see Eq. (13)) does not function as intended, because the first sub argument of the closure coefficient blending function (see Eq. (14)) is sensitive to the large freestream turbulence length scales. As a result, the intended zonal treatment is not achieved. We consider a modified definition of the first blending function, *F*_{1v}, which is calculated using the function argument, arg_{1v}, defined in Eq. (20) and activates SKO-regime in the viscous sublayer only. We now compare results using the *default* blending function (SST2003) and *modified* blending function (SSTF1V) in terms of their regional treatment (SKO versus SKE), as well as their KE loss, and wake mixing character. Both simulations correspond to the high inlet turbulence condition ($TI\u224813%$; *ℓ* ≈ 1.0 *C*_{x}). Consider now the blade-to-blade midspan section of the NGVs, shown in Fig. 11. Inside the contoured regions, the closure coefficients are closer to SKO-regime than to SKE-regime (i.e., *F*_{1} or *F*_{1v} > 0.50), and outside the contoured regions the closure coefficients are closer to SKE-regime (i.e., *F*_{1} or *F*_{1v} ≤ 0.5). The edge of the contoured regions can therefore be viewed as the location where the turbulence model is switched. We now review the zonal treatment in SST2003 and SSTF1V.

For the SST2003 model, the region upstream of the NGVs, and the PS-half of the vane passage use SKO, whereas the SS-half of the vane passage and the wake use SKE. The gaps in the SKO-regime zone near the leading edge and on the PS are explained by the film cooling holes (internals obscured). For the modified blending function (SSTF1V model), only the near-wall regions use SKO-regime (see detail view of trailing edge, leading edge and PS). Overall, the behavior of the SSTF1V model appears more predictable and well-behaved in the presence of film cooling and large length scales. We briefly consider the effect of the differing coefficient blending on predictions of KE loss and wake mixing.

The plane-averaged KE loss coefficient increased by 1.8% for SSTF1V, when compared to SST2003 (average across planes 2–4). The increase was primarily driven by the profile-region KE loss coefficient, which increased by 2.1%, while the endwall-region KE loss coefficient only rose by 0.2%. Changes in wake width/depth were +1.8%/−0.7% for the restacked aerodynamic wake, and +1.2%/−1.0% for the restacked thermal wake. We conclude that KE loss and wake spreading/decay predictions are not significantly affected by the different zonal treatment.

## Appendix C: Inlet Turbulence Intensity

Turbulence model sensitivity to inlet TI was investigated by performing additional simulations of the NGVs for an inlet TI of 1.3% (reduced from $TI\u224813%$; length scale unchanged). For the sake of brevity, we consider the TI sensitivity of *profile-region* KE loss coefficient, and *aerodynamic* wake width only. Endwall losses were not significantly affected, and *thermal* wake width was approximately proportional to the *aerodynamic* wake width.

Turbulence model predictions of profile-region KE loss coefficient, *ζ*_{p}, and the pitch-normalized and rescaled aerodynamic wake width, *σ*/*s*, are shown in Fig. 12 for inlet conditions of low $(TI\u22481.3%)$ and high ($TI\u224813%$) turbulence. The values are shown for an axial plane located 0.75 *C*_{x} downstream of the midspan trailing edge (plane 4; see Fig. 3). The results at this station were representative of the overall trend. The vertical dashed lines correspond to the loss/width predicted by the SST2003 model at low TI.

First consider the profile-region KE loss coefficient, *ζ*_{p}, which is shown on the left side of Fig. 12. Regrettably, no experimental validation data was taken at the low TI inlet condition. Consulting experimental studies concerning the effect of inlet TI on the profile loss of turbine cascades, we expect an increase of profile loss at high turbulence conditions. Ames and Plesniak [41] measured an 8.4% rise in profile loss (for inlet conditions changing from a TI of 1.1% with *ℓ* ≈ 0.84 *C*_{x} to a TI of 12.0% with *ℓ* ≈ 0.43 *C*_{x}), and Folk et al. [6] reported a 37% increase of profile loss (for TI increasing from approximately 1.5–10%; constant length scale, *ℓ* ≈ 0.18 *C*_{x}). Surprisingly, the SST2003 model predicted only a 2.7% increase of the profile-region KE loss coefficient for the current NGV test case (for TI increasing from 1.3% to 13%; *ℓ* ≈ *C*_{x}). The outcome was unexpected, because flat plate tests using the SST2003 model (see Fig. 1) showed a 46% increase of the dissipation coefficient between the low and high TI condition, which corresponded to a slight overprediction, when compared to the experimental flat plate data of Folk et al. [6]. This reiterates the dominant role of stress limiting conditions in the predictions of NGV loss. For the realizable models, the increase of profile-region KE loss at high TI was 12.5% (SSTDR), 16.8% (BSL-ARSM), and 17.3% (SSTPSR), which is somewhat better-aligned with the existing literature. SA predicted a similar increase of 17.1%, while BSL, SKO, and SST1994 saw profile-region KE loss rise by 38.7%, 40.7%, and 96.9%, respectively. The conclusion is that all turbulence models correctly predict increasing KE loss with increasingly TI, but with increasing variation in absolute value between models as TI increases (9.6% of the model average at low TI; 73.9% of the model average at high TI).

Now consider the pitch-normalized and rescaled aerodynamic wake width, *σ*/*s*, which is shown on the right side of Fig. 12. At high TI, Ames and Plesniak [41] reported an 11.3% increase of wake width, while Folk et al. [6] found that wake width was increased by approximately 28%. In the framework of an eddy viscosity model, this can be explained by recalling that kinematic eddy viscosity (i.e., diffusivity of momentum) scales linearly with TI for a constant length scale. For the current test case, SST2003 had the lowest sensitivity to inlet TI, with wake width increasing by only 2.1%. For the other models, wake width increased by between 29% and 189%. The increases were approximately proportional with the corresponding change of the profile-region KE loss coefficient (for TI increasing from 1.3% to 13%). This is in line with our earlier observation that wake width and profile-region KE loss are inherently linked for linear eddy viscosity models, since both are driven by the TKE production around the SS peak. Indeed, the only *non-linear* eddy viscosity model (i.e., BSL-ARSM) is the exception to this rule, with a disproportionate increase of wake width, compared to the rise of KE loss. The underlying mechanism is not currently understood but is likely related to spatially-varying degrees of *coherence* between the turbulent stress and mean flow strain tensors, which allows for independent modulation of TKE production (see discussion by Moore and Moore [27]).

## Appendix D: Length Scale Definition

*integral*length scales into modeled

*turbulence*length scales (and vice versa). This is relevant because length scales are typically used to estimate the turbulence dissipation rate when setting RANS inlet boundary conditions. While the expression of Taylor [20]

*A*remains debated. Cha et al. [5], Folk [31], and Amend et al. [32] used

*A*= 1, when comparing RANS-predicted turbulence length scales to experimentally-measured integral length scales.

*A*= 1 is also used by the ansys cfx solver [22]. The outcomes from the studies using

*A*= 1 are inconclusive, however, with Amend et al. [32], Folk [31], and Cha et al. [5] reporting that CFD-predicted length scales were 15%, 108%, and approximately 200% greater, respectively, than experimentally-measured integral length scales. It is interesting to note that the study of Cha et al. [5] would have achieved a good collapse had the proportionality constant of Moore and Moore [27] been used (they estimated

*A*= 0.32 based on hotwire measurements of integral length scale, and

*k*, with $\u03f5$ estimated from turbulence spectra). The ansys fluent solver uses

*A*≈ 0.164, while Wilcox [24] suggests that a value

*A*= 0.09 is most appropriate. The point is that while the chosen value of the proportionality constant

*A*can have a strong influence on CFD predictions of KE loss and wake width, values in common use vary by an order of magnitude.

Additional uncertainty arises due to variation in the definition of the experimentally-measured integral length scale. In the current study, the integral length scale was calculated as the product of the local velocity magnitude and an average integral time scale. The average integral time scale was defined as the arithmetic mean of the components in *x*, *r*, and *θ*, each of which was calculated from the integral of the normalized autocorrelation function of the corresponding velocity component from zero time-shift to the first zero crossing. To build understanding in this field, and where back-to-back comparisons are desired, we emphasize the importance of reporting both the method for calculating experimentally-measured *integral* length scales, and the proportionality constant used to calculate $\u03f5$ from *k* and *ℓ*.