A characterization of unsteady effects for transonic turbine airfoil limit loading

Airfoil description

Two transonic turbine airfoils, shown in Figure 2 were used in this study: the first from the thesis of Sooriyakumaran at Carleton University in Ottawa, Canada, ⁸ and the second from the work of Sieverding at the European Research Project BRITE/EURAM CT96-0143 on “Turbulence Modeling of Unsteady Flows in Axial Turbines”. ⁹ The profile adapted from Sooriyakumaran is the main focus of this paper; however, the airfoil profile from Sieverding was also analyzed to ensure adequate validation of the transient nature of the flow dynamics around the trailing edge. Table 1 shows the relevant geometric parameters for each profile. As shown, the airfoil from Sieverding’s study has a much larger trailing edge thickness as well as higher unguided turning and stagger angle than the profile from Sooriyakumaran; however, the pitch-to-chord ratios for both profiles are similar. The airfoil testing aspect ratio for the results for the wind tunnel testing performed by Sooriyakumaran were over 1.5 to ensure end-wall effects could be ignored and the midspan flow assumed to be fully two-dimensional. ¹⁰ OPEN IN VIEWER

OPEN IN VIEWER

Table 1.

Airfoil geometric characteristics.

Airfoil	Sooriyakumaran	Sieverding
Aspect ratio, h/C t	1.5	0.7
Pitch-to-chord ratio, s/C t	0.66	0.70
Blockage ratio, t/o	0.05	0.22
Stagger angle, ζ	27°	50°
Inlet metal angle, β1	18°	0°
Unguided turning, θ UG	9°	14.5°
Zweifel number, Z	0.827	0.692

Figure 3 shows the computational domain used for this study. The inlet plane was positioned 1.5 axial chord lengths upstream of the airfoil leading edge to ensure uniform and fully developed flow. The extent of the outlet boundary, positioned seven axial chord lengths downstream of the trailing edge, was chosen specifically for the limit and supercritical loading conditions to ensure numerical reflections did not influence predictions near the evaluation area. ¹¹ The evaluation plane is taken to replicate the experiment used for validation and is positioned 1.5 axial chord lengths downstream of the airfoil leading edge. The top and bottom surfaces located a half pitch above and below the airfoil surface were periodic with a one-to-one mesh connectivity. Similarly, the positive and negative spanwise planes used for 3D predictions were periodic with a total thickness of 10% of the airfoil chord and one-to-one connectivity.OPEN IN VIEWER

Boundary and initial conditions

Boundary conditions used were set to replicate the experimental conditions of Sooriyakumaran. ⁸ The inlet boundary was prescribed a constant total pressure and total temperature, the outlet boundary an average static pressure, the top and bottom along with the positive and negative spanwise extents were set as periodic, and the airfoil surface a smooth, no-slip adiabatic wall. Incidence angle was set to 14°, and inlet turbulence intensity and length scale were set to 4% and 10% of airfoil pitch, respectively. Although uniform artificial inlet turbulence is unlikely in actual turbine engines; the authors choose this approach to ensure that an analogous comparison of downstream modeling approaches could be made, specifically around the trailing edge where shocks are most prevalent at these conditions. The development of true inflow turbulence is considered outside the scope of this paper and will be considered in future work. Outlet average static pressure was adjusted depending on the desired inlet-total to outlet-static pressure ratio. Reynolds number based on the velocity at the evaluation plane and the true chord was Re _c ∼ 1.0 × 10⁶.

It should be noted that although the simulations were setup to replicate the experimental conditions of Sooriyakumaran, not all exact physical conditions can be precisely reproduced. The facility in Ottawa, Canada is a blow-down configuration, and can sustain stable transonic flow for up to 60 seconds. Both Reynolds number and exit Mach number can be varied independently due to an ejector-diffuser assembly. ¹⁰ The cooled, dry air is discharged through a wire mesh that generates a known turbulence characteristic and is then directed through the test section. For a detailed description of the experimental facility, an interested reader is referred to either the work of Sooriyakumaran ⁸ or Taremi. ¹²

A few key differences between the experiment and CFD are listed next. First, the dimensions of the airfoils as well as the working fluid properties have been scaled from engine conditions to appropriate environments for the experimental facility. This is not of concern as the flow Reynolds number was used to match the conditions between simulation and experiment. Second the simulations assume a single midspan blade profile, which is certainly not the case in reality. However, a previous study was performed by Corriveau ¹⁰ that showed that an aspect ratio, h/c, of 1.5 is sufficient to allow for a region of essentially two-dimensional flow at the midspan. Lastly, the blades used in the experiment were manufactured using a direct metal laser sintering rapid prototyping process. This additive process was shown to result in a surface roughness of about 10 um, which was not accounted for in the simulations.

The inlet conditions were used to initialize the RANS simulations. However, it should be noted that these simulations employed grid sequencing, so the domain initialization did not significantly influence the convergence rate. Each transient modeling approach was initialized with the respective RANS solution. For example, the URANS, DDES and model free simulations for the critical loading condition were initialized using the RANS converged critical loading solution. A detailed description of the grid sequencing process has been documented by Owen. ¹¹

Modeling strategies

Although it is well known that the flow within a turbine blade passage is intrinsically unsteady due to vortex formation around the trailing edge, RANS modeling is the most extensively used simulation method in the gas turbine industry due to lower computational cost. ¹³ The unsteady RANS (URANS) modeling approach provides an improvement to the RANS method by stepping the solution in time to capture vortex formation and other transient phenomena with the smallest amount of increased computational effort. A further improvement in modeling capability is the Detached Eddy Simulation (DES) approach which combines URANS and Large Eddy Simulation (LES) strategies. ¹⁴

As mentioned, DES is a hybrid modeling approach that uses a RANS solution within the boundary layer and irrotational flow regions while the detached and freestream flow areas are solved using a modified LES. The LES methodology used is dependent on the specific DES model; here the SST k-omega formulation, obtained by modifying the dissipation term in the transport equation for turbulent kinetic energy, was employed to capture the unsteady separated regions. The delayed DES formulation (DDES) was used to enhance the ability of the model to distinguish between the LES and RANS regions. ¹⁴

Lastly, to avoid the empirical nature of the turbulence formulations as well as to ensure simulation fidelity, we tried to solve the flow field using an approach which does not rely on any turbulence modeling approaches. This methodology was similar to a direct numerical simulation approach; however, the simulation domain was kept as two-dimensional. Therefore, the authors have termed this as a pseudo-DNS or model free strategy. Using the results from the URANS simulation the Kolmogorov and Taylor length scales were estimated to be of the order 10⁻⁵ and 10⁻⁴ m respectively. Using these estimates it was determined that the average grid size was roughly half the Taylor length scale and roughly 10 times the Kolmogorov length scale. This was determined to be sufficient based on previous numerical investigations.^15,16

All simulations employed a coupled, implicit compressible algorithm to solve the equations for mass, momentum, and energy. Discretization was second order for the source and diffusion terms for all methods except DDES, which was a hybrid second-order upwind/bounded central-differencing scheme. All models used a second-order upwind scheme for calculation of the convection terms for the momentum and turbulence equations. The inviscid, convective fluxes were calculated using Liou’s AUSM+ flux-vector splitting scheme. ¹⁷ The Algebraic Multi-Grid (AMG) linear solver with a V-cycle, Gauss-Seidel relaxation scheme and bi-conjugate gradient stabilizer acceleration method were employed. ¹⁸

Convergence criterion and time step

All precursor RANS solutions were run until the flow had reached full convergence such that the normalized residuals of all transport quantities were steady and dropped to four orders of magnitude; in addition, the local total pressure loss coefficient did not show fluctuations over ±0.0001 for 1500 iterations which typically occurred around 3500 iterations.

To ensure minimal influence of numerical errors as well as to improve the precision and total run-time of each transient simulation, an adaptive time stepping algorithm was used for each transient simulation. The adaptive time-stepping algorithm adjusts the time-step automatically during the run to obtain a user specified temporal resolution throughout the domain.

The local time-step (Δτ) is computed considering both the Courant number and the Von Neumann stability condition, equation shown below:

Δ τ = min ( C F L * V ( x ) λ max , V N N * Δ x 2 ( x ) ν ( x ) )

(1)

Here CFL is the user-specified Courant number (dimensionless), V is the cell volume (m³), VNN is the Von Neumann number (VNN = 1), Δx is a characteristic cell length scale (m) and ν is kinematic viscosity (m²/s). λ_max is the maximum eigenvalue of the system as follows:

λ max = | u ⋅ a | + c | a |

(2)

Here an upper limit to the convective CFL number of 0.9 was selected as the criterion to determine the solution time step. Table 2 shows a summary of the mean time-step sizes for each simulation. The required time step to maintain sufficient temporal accuracy is dependent both on the loading condition and modeling strategy employed. The time step size was found to be in the range of 2.6 × 10⁻⁸ to 7.0 × 10⁻⁸ seconds.

OPEN IN VIEWER

Table 2.

Average time step.

Model strategy	Loading state	Avg. Δt (sec)
URANS	Critical	7.04 × 10−8
	Sublimit	6.45 × 10−8
	Limit	5.73 × 10−8
	Supercritical	4.92 × 10−8
Model-free (MF)	Critical	6.21 × 10−8
	Sublimit	6.14 × 10−8
	Limit	5.24 × 10−8
	Supercritical	4.87 × 10−8
DDES	Sublimit	2.85 × 10−8
	Limit	2.69 × 10−8

Grid description and convergence

As previously mentioned all transient simulations were initialized using a RANS solution. The grid employed for each RANS simulation was generated through an adaptive griding internal algorithm detailed next.

If | ∇ ϕ − ∇ ϕ min ∇ ϕ max − ∇ ϕ min | > 0.01 Cell Area or Volume → max ( G o 2 , G L P L )

(3)

Equation (3) shows the governing process of the adaptive mesh refinement scheme. Here ϕ denotes the scalar quantity of interest (absolute total pressure, total temperature, density, turbulent kinetic energy, turbulent eddy viscosity and the specific dissipation rate of turbulence), G _o the original cell area or cell volume prior to adaptation, and G _LPL is the area or volume of the last prism layer. This was used as the limiter of refinement to ensure a one-to-one ratio between the prism layers and the core grid.

As shown previously, ¹¹ the adapted grid algorithm proved robust in capturing the blade loading of four different transonic turbine airfoils at sublimit loading conditions, proving the ability to accurately capture the trailing edge shock system and the subsequent influence on the adjacent airfoil suction side boundary layer. Figure 4 shows an example of the grid density evolution throughout the AMR process.OPEN IN VIEWER

Each simulation began with a baseline mesh of roughly 70,000 elements. The core grid comprised unstructured, polyhedral elements with prismatic layers along the airfoil surface. The distribution of prisms normal to the airfoil remained constant with 40 layers, a hyperbolic growth rate, a near wall thickness achieving a max y+ less than 0.5, and a total thickness equal to 40% of the trailing edge thickness. It should be noted that the distribution along the tangential direction of the airfoil surface was changed with increasing refinement to maintain an aspect ratio of less than two. Nine AMR cycles were used for each RANS simulation, and additional AMR cycles only proved to increase computational resources unnecessarily as the increase in number of elements was less than 1%.

Although this grid methodology provided a grid independent solution for the RANS simulation approach there was no guarantee that it would suffice for the other models. Therefore, an additional grid independent study was performed for each loading condition and transient modeling strategy. Using the adapted mesh from the RANS solution as a baseline, the grid was isotropically refined by 50 and 25%, roughly doubling the total cell count in each case.

Figure 5 illustrates the differences in downstream predicted total pressure loss for each grid level. As can be seen, the change in total pressure loss distribution is minimal between grid refinements, producing a change in mixed-out total pressure loss coefficient of less than 2%. It should be noted that a grid verification procedure was applied to each loading condition for each transient strategy; the image below shows an example. As there was no significant improvement in using the refined grids the adapted grid of the RANS solution was used for each transient simulation.OPEN IN VIEWER

2D/3D RANS and URANS comparison

Prior to the main results section, the authors would like to highlight the differences between 2D and 3D RANS and URANS simulations. Given that the computational domain is modeled as periodic in the spanwise direction, the flow is expected to be similar for both 2D and 3D approaches. Figure 6 shows the surface isentropic Mach number distribution as well as the downstream pitch-wise total pressure loss distribution for both 2D and 3D simulations. Solid lines are the RANS predictions, dashed lines are the URANS predictions, and solid black squares are experimental data taken from the work of Sooriyakumaran.OPEN IN VIEWER

Looking at the blade loading distribution, it is clear there is no noticeable difference between 2D and 3D simulations. On the pressure side the flow accelerates consistently until the trailing edge expansion while the suction side shows a slight deceleration at the leading edge followed by nearly constant acceleration until the point of shock impingement. Following the shock/boundary layer interaction there is a period of nearly constant pressure until the trailing edge. This flow pattern is consistent between modeling domains.

Similarly, the downstream total pressure loss distribution shows only minor deviations between 2D and 3D approaches. The RANS simulations predict a significantly higher loss distribution in the middle of the wake with smaller width compared to URANS. As no major differences were seen between 2D and 3D results for both RANS and URANS modeling approaches, the remainder of the paper will only include the results from the 2D simulations.

Validation

As previously mentioned, two transonic turbine airfoils were used for numerical validation. Both airfoil geometries have been previously described and are shown in Figure 2. The total inlet to outlet static pressure ratio used by Sooriyakumaran was 2.25, a pressure ratio within the sublimit loading range, corresponding to an isentropic exit Mach number of 1.15 and for the profile used by Sieverding PR = 1.5, a subcritical loading condition corresponding to an isentropic exit Mach number of 0.79.

Figure 7 illustrates the predicted airfoil surface isentropic Mach number distribution for each profile compared to their respective experimental data. Extracted test information is shown as solid black squares while all simulation results are shown as various line distributions; RANS - solid lines, URANS - dashed lines, model free - dotted lines and DDES - dashed dotted lines.OPEN IN VIEWER

Overall distributions are predicted well for both profiles and loading conditions. Results for the profile taken from the work of Sooriyakumaran show nearly identical predictive capability for the full pressure surface as well as the majority of the suction surface. Differences arise along the suction surface at roughly 0.9 axial chord lengths (C _ax ) aft of the airfoil leading edge, at which point the model free approach predicts a smearing of the shock impact on the suction surface boundary layer. This is a product of an oscillatory effect caused by the unsteady vortex shedding at the airfoil trailing edge, and additional details will be provided. Similarly, the predicted results for the profile used by Sieverding agree quite will with the experimental test data for the majority of the airfoil surface. Discrepancies between modeling approaches arise at the point of shock impingement along the suction surface similar to what is observed for the profile used by Sooriyakumaran. This is believed to be primarily a product of the unsteady propagation of vortices near the trailing edge.

Next a look at the downstream total pressure loss distribution for the profile used by Sooriyakumaran as well as the base pressure distribution of the profile used by Sieverding are presented in Figure 8. Focusing first on the loss distribution for the profile used by Sooriyakumaran, it is clear that neither the URANS nor the model free approach predict the correct minima of total pressure loss. The total wake width is slightly increased compared to the experimental results and the free stream losses are significantly over predicted. For this loading condition, this suggests an inability to accurately predict the time averaged vortex composition. This is due to a combination of effects: such as base pressure region size, subsequent vortex core size and distribution as well as downstream shock interactions and is presented in more detail in the transient section. The over prediction of freestream losses is similar for all modeling approaches; however, in contrast both the RANS and DDES strategies show good agreement with the location and magnitude of the loss minima. The DDES approach shows the best overall agreement of the methods considered, predicting a slightly wider wake width than the RANS approach while maintaining the correct magnitude of loss at the wake centerline.OPEN IN VIEWER

As described by Sieverding the pressure distribution along the trailing edge is characterized by three pressure minima: two caused by the Prandtl-Meyer expansion around the suction and pressure sides prior to flow separation and a third at the trailing edge center caused by the enrolment of the separating shear layers into a vortex core. ¹⁹ In this case, none of the modeling approaches seem to be able to replicate the complex flow phenomena occurring at the trailing edge exactly. RANS modeling shows the largest discrepancy, predicting first a sharp rise in base pressure following the separation of the suction and pressure side shear layers followed up by a large pressure plateau towards the trailing edge centerline. The URANS approach is the only method to successfully predict the location of the three pressure minima, although the magnitude is significantly lower than the experimental data. The model free approach was unable to predict the pressure side expansion location and subsequent pressure minima, but the location of the centerline and suction side minima are correctly captured. It should be noted that the predicted base pressure distributions shown above are similar to previous numerical results using URANS and RANS methods. ¹³

Mean flow

Figures 9 and 10 show the time averaged, mixed-out total pressure loss coefficient as well as the mass flow averaged total flow turning at the evaluation plane for each of the loading and modeling strategies considered. All methods, with the exception of the URANS approach, predict a gradual rise in mixed-out total pressure loss coefficient from critical to limit loading. The sharp rise in loss prediction for the URANS model at the critical loading condition is explored in detail later; however, in summary this is due to the inability of the URANS method to accurately predict the transition from trailing edge dominated vortex formation to shock dominated vortex formation.OPEN IN VIEWER

OPEN IN VIEWER

Figure 10.

Time averaged, mass flow averaged total flow turning at each loading condition.

For all methods considered there is a sharp decline in loss coefficient from limit to super critical loading. This was unanticipated, as the total pressure loss was expected to increase with pressure ratio.^1,3,5 However, as the passage becomes fully choked, the oblique Trailing Edge Pressure Side Shock (TEPSS) caused by the rapid expansion around the trailing edge, completely clears the adjacent airfoil base pressure region resulting in reduced total pressure loss as well as flow turning. The interactions of the shock system with the blade boundary layer and wake structures becomes small and thus the predicted total pressure loss is reduced. Although a decrease in total pressure loss is desirable, there is a significant detriment in the ability of the airfoil to produce the desired work output, as this is directly related to the total flow turning. Furthermore, the large change in outlet flow angle can drastically decrease the efficiency of downstream blade rows through a change in incidence angle. A quick interpolation of results from Owen, ⁶ shows that for a change of incidence by less than 2°, such as the change from critical to limit loading, total pressure loss is expected to rise by roughly 1%. However, as the incidence angle is significantly changed ( > 9 ° ) , such as is the case for supercritical loading, the estimated loss increase is roughly 25%. Clearly this is not desirable, highlighting the importance for aero designers to understand the limit and supercritical loading conditions.

Shown in Figure 11, are the predicted blade loading distributions for each modeling strategy at each loading condition. The available experimental data of Sooriyakumaran has been included for completeness. As a reminder experimental test results are only available for the sublimit loading condition.OPEN IN VIEWER

Blade loading distributions along the airfoil pressure surface as well as the suction surface prior to the location of shock interaction ( ∼ 0.8 C a x for critical loading and ∼ 0.9 C a x for sublimit) are predicted to be identical across loading conditions and modeling approaches. This is a result of flow choking cause by the trailing edge pressure surface shock (TEPSS). The effects of unsteady shock propagation and trailing edge vortex shedding are found to have a minimal effect on the flow field and subsequent blade loading. Downstream of the pressure side shock, highlighted in Figure 12, there are clear differences in the ability of each modeling approach to predict the shock propagation and subsequent shock/boundary layer interaction.OPEN IN VIEWER

Unsteady propagation of the TEPSS caused by the interaction of the shock system with the trailing edge vortex shedding process causes the surface isentropic Mach number distribution to change. As the pressure surface shear layer is shed and rolled into a vortex core the shock system is weakened, resulting in an oscillatory shock effect. This system is then propagated downstream towards the adjacent airfoil causing a subsequent shock-boundary layer interaction which, in a time averaged sense, results in a smoothing of the shock impact.

This is most significant at lower pressure ratios where the predicted base pressure region is smaller, such as the critical loading condition. As illustrated through numerical Schlieren flow visualizations in a later section, the pressure and suction side shear layers do not achieve sufficient momentum to create a well-established base pressure region. This causes the trailing edge pressure side shock to form nearer the airfoil trailing edge resulting in significant influence from the vortex formation process. As the pressure ratio is increased, the base pressure region becomes established and causes vortex formation to occur downstream of the trailing edge shock system. All modeling approaches, with the exception of the model free strategy show minimal oscillation of the shock location as the base pressure region becomes established. Therefore, RANS/URANS and DDES all predict similar loading distributions for the sublimit loading conditions and above. As will be shown later, the model free approach predicts a much smaller base pressure region than the others, as is apparent in the base pressure distribution as well as flow visualizations. At supercritical loading all models predict near identical surface pressure distributions as the flow upstream of the TEPSS is fully choked and temporally stationary.

To further investigate the causes of the differences in vortex formation and subsequent downstream total pressure distributions, plots of the wall normal velocity on both the pressure and suction sides at 98% chord are presented in Figures 13 and 14. Due to the strong acceleration along the pressure surface, the boundary layer near the trailing edge of the pressure side is predicted to be laminar and near identical across loading conditions. The flow in this region is fully choked; therefore, any increase in pressure ratio only changes the flow downstream of the trailing edge pressure surface shock. The total boundary layer thickness, estimated as 99% of the freestream velocity, was predicted to be roughly 0.92 mm for the RANS, URANS and DDES simulations and 0.31 mm for the model free approach.OPEN IN VIEWER

OPEN IN VIEWER

Figure 14.

Suction side trailing edge boundary layer profiles.

Unlike the pressure side, the suction side of the airfoil surface shows a strong dependence on load condition and modeling strategy, specifically for the limit loading condition. A summary of the predicted boundary layer thickness for each condition and strategy is shown in Table 3. At the lowest pressure ratio, the flow is clearly attached and turbulent with the exception of the model free strategy which alludes to a transitionary state. Similar to the pressure side RANS, URANS and DDES approaches all show similar boundary layer thickness and distribution for the sublimit loading condition with a total thickness of roughly 1.5 mm.

OPEN IN VIEWER

Table 3.

Suction side boundary layer thickness summary (mm).

Model	Critical	Sublimit	Limit	Supercritical
RANS	2.35	1.49	2.80	1.46
URANS	2.35	1.49	2.16	1.46
MF	0.54	0.60	3.46	0.54
DDES		1.49	2.67

As the pressure ratio is increased to limit loading the boundary layer profile shows a period of separation. This is a product of the influence of the trailing edge pressure side shock interacting with the suction side shear layer that forms one side of the base pressure region. The shock/base pressure interaction causes acoustic waves to propagate upstream creating a brief period of instability in the boundary layer. Although each modeling strategy shows separation, the magnitude is significantly different between strategies with the model free approach showing the largest region of flow recirculation. The approximate size and location of the period of separation requires a deeper look at the full flow picture and is investigated further with flow visualizations later in this report.

Figure 15 shows the predicted base pressure distribution for each modeling strategy and loading condition. For critical loading, the flow at the trailing edge is characterized by attached vortices formed by the alternating shedding of the pressure and suction side shear layers. The trailing edge pressure surface shock gains sufficient expansion to create an established base pressure region, resulting in a period of nearly constant pressure towards the trailing edge centerline. Local pressure minima are seen at the locations of the shear layer separation points and are a result of the Prandtl-Meyer expansion around the trailing edge. All modeling approaches predict the general trend of the pressure distribution similarly, although the unsteady formulations show a significantly lower base pressure than the steady RANS simulation for the critical loading condition. As the pressure ratio is increased towards limit loading the magnitude of the base pressure is lowered and is consistent for all modeling approaches. A summary of the normalized time and surface averaged base pressure is listed in Table 4 and presented alongside the well-known base pressure correlation created by Sieverding ²⁰ in Figure 16.OPEN IN VIEWER

OPEN IN VIEWER

Table 4.

Area averaged base pressure summary.

Model strategy	Loading condition	Ps2/Po1	P b /Po1
RANS	Critical	0.498	0.466
	Sublimit	0.444	0.374
	Limit	0.408	0.300
	Supercritical	0.257	0.194
URANS	Critical	0.496	0.352
	Sublimit	0.444	0.370
	Limit	0.408	0.302
	Supercritical	0.257	0.194
MF	Critical	0.496	0.362
	Sublimit	0.444	0.323
	Limit	0.408	0.257
	Supercritical	0.259	0.176
DDES	Sublimit	0.444	0.390
	Limit	0.408	0.302

OPEN IN VIEWER

Figure 16.

Comparison with the base pressure correlation proposed by Sieverding.

Results are distinguished as follows: RANS simulations as red squares, URANS simulations as purple diamonds, model free simulations as green triangles and DDES as orange circles. The solid lines show the correlation by Sieverding for the geometric properties nearest the investigated airfoil and the dashed line signifies equal base and downstream static pressures to help with analysis.

As the total pressure ratio is increased the predicted base pressure decreases for all modeling approaches. During critical loading each unsteady model predicts a significantly lower base pressure than steady RANS. As the pressure ratio is increased to sublimit and limit loadings the base pressure is predicted similarly for each of the wall modeled simulations, suggesting the flow is statistically stationary within the base pressure region. The model free approach predicts a much lower base pressure than the others for each of the loading conditions. This is a product of the size of the base pressure region, which is directly related to the magnitude and state of the pressure and suction side boundary layers near the trailing edge. As shown previously, the predicted model free boundary layers are much smaller than the other approaches. This causes the shear layer separation at the trailing edge to occur more rapidly, creating a stronger pressure side expansion. As the expansion becomes stronger the subsequent base pressure becomes lower.

Figure 17 shows the predicated time averaged total pressure loss distributions for each pressure ratio and modeling strategy. As observed in the mixed-out mass flow averaged total pressure loss coefficient, downstream losses increase with an increase in pressure ratio until supercritical loading. This is due to a combination of increased wake and freestream losses. During critical loading freestream losses are predicted to be small compared to wake losses as a result of how the vortices are being shed from the trailing edge. When the base pressure region is attached to the trailing edge, vortices shed quickly leading to increased mixing and wake losses.OPEN IN VIEWER

Increase of the pressure ratio to sublimit loading causes the flow expansion around the trailing edge to become stronger, leading to a larger and more established base pressure region that pushes vortex formation downstream as well as decreases the frequency of the vortex shedding. Although there is a reduction in wake losses, shock and freestream losses are increased. This is a product of the strong oblique trailing edge shock system, with the pressure side shock interacting with the suction surface boundary layer and the suction side shock propagating through the downstream wake. For limit loading the loss is further increased due to the complex flow near the trailing edge and base pressure region. As the TEPSS passes through the base pressure region acoustic waves propagate upstream along the suction surface followed by flow separation and increased wake width and depth. Also included is the impact of the suction side shock on the downstream wake profile. Predictions for super critical loading show that as the TEPSS departs the base pressure region, the influence of vortex formation becomes small and freestream losses are reduced for the majority of the domain.

Looking at the difference in modeling approaches, it is clear that the RANS simulations predict the smallest wake width and largest wake depth of the modelling strategies considered for all loading conditions. This is due to a lack of vortex/shock interactions. With the exception of critical loading, the URANS and model free approaches show a much smaller wake depth; however, due to vortex formation they predict an increase in wake width compared to RANS. The DDES approach, shown only for the sublimit and limit loading conditions, shows the widest wake prediction as well as best correlation with available experimental data.

During critical loading the URANS method predicted an unrealistic downstream loss profile. It should be noted the authors preformed a rigorous grid refinement study as well as manual time step adjustment to repair the downstream loss prediction. It was concluded that the URANS method failed to correctly predict the transition from near trailing edge vortex formation to shock trailing edge vortex formation at the edge of the base pressure region. Further analysis results are presented and discussed later in this paper.

Transient phenomena

Shown in Figure 18 are the normalized power spectral densities of the mass flow averaged total pressure loss coefficient at the evaluation plane for the critical, sublimit and limit loading conditions. Supercritical loading is not presented as the loss coefficient for each modeling approach remained temporally constant. This is believed to be a result of the choking of the flow caused by the trailing edge pressure surface shock. As can be seen there is a peak energy that clearly defines the vortex shedding frequency for each modeling approach.OPEN IN VIEWER

A summary of the peak vortex frequency, calculated Strouhal number and vortex shedding period is given in Table 5. The Strouhal number was computed using equation (4) with f _vs being the peak vortex shedding frequency, d _te the trailing edge diameter and V_2,is the mass flow averaged isentropic flow velocity at the evaluation plane.

S t = f v s d t e V 2 , i s

(4)

OPEN IN VIEWER

Table 5.

Power spectrum analysis summary.

Model strategy	Loading	Peak freq. (kHz)	Strouhal #	Vortex period (ms)
URANS	Critical	60.3	0.241	0.017
	Sublimit	58.1	0.217	0.017
	Limit	60.7	0.217	0.016
MF	Critical	70.2	0.279	0.014
	Sublimit	84.7	0.316	0.012
	Limit	63.0	0.226	0.016
DDES	Sublimit	35.9	0.134	0.028
	Limit	46.2	0.165	0.022

For each modeling strategy, the average adaptive time step required to produce a convective CFL number of less than one was smaller for higher total pressure ratios. Although there is a wide range of predicted Strouhal numbers between modeling strategies, they are all within the bounds of similar experimental and computational investigations. ¹⁹ The model free simulation shows a much higher predicted vortex shedding frequency than the other unsteady methods for all loading conditions, with a peak difference occurring during sublimit loading. As discussed previously, this is a product of the state of the trailing edge boundary layer. The results of the URANS simulations predict lower Strouhal numbers than the model free approach; however, the DDES simulations predict Strouhal numbers to be much lower than either of the other methods. To help clarify the causes of the differences in vortex shedding frequency between modeling strategies flow visualizations of the normalized density gradient at various phase states are presented in the following section.