Numerical Investigation of the Wake Vortex-Related Flow Mechanisms in Transonic Turbines

As the core equipment of the power generation system, a gas turbine is an indispensable energy-converting device in the national industry. The flow inside a high-pressure turbine (HPT) is highly unsteady, which has a great influence on the aerothermal performance and structural strength. To better clarify the flow mechanism and guide the advanced design, the basic flow characteristics of transonic turbines are investigated in the paper by a modified scale-adaptive simulation (SAS) model based on the shear stress transport (SST) turbulence model. The numerical results reveal the formation and development of the secondary flow structures such as wake vortex, pressure wave, shock wave, and the interactions among them. The length and frequency characteristics of wake are in good agreement with the large eddy simulation (LES) and the experimental data. Based on the detailed flow information, the local loss analysis is performed using the entropy generation rate. In summary, the wake vortex-related flow is the main origin of unsteadiness and entropy loss in high-pressure turbine cascade.


Introduction
Due to highly complex geometries and viscous effect, various secondary flows which cover a range of length scales and frequencies are formed in a gas turbine. The understanding of its flow mechanism is critical to the advanced design. The high cost of experimental equipment and the limitations of measurement technology make numerical research increasingly important.
Simulation technology is an important means to support the research and development of high-performance HPT. Over the past 30 years, significant progress has been made in computer technology and computational fluid dynamics (CFD) technology, especially in the full-scale three-dimensional CFD capability for turbomachinery. As a result, the design cycle of a gas turbine engine has been greatly reduced. In order to further improve the aerothermal performance and reduce the fuel consumption, the total pressure ratio and the inlet temperature of HPT are increasing, but its size becomes more compact. Therefore, in precise design of modern turbines, more attention should be paid to the flow details, rather than the average parameters obtained by Reynolds averaged N-S simulation. However, the existing theory for turbine aerodynamic design is mainly based on statistical modeling, ignoring the inherent unsteady flow characteristics. High-fidelity numerical simulation is an effective means to study the complex flow phenomena inside turbine components, and accurate simulation provides a solid foundation for accurate loss prediction. In recent years, the global hybrid models such as scaleadaptive simulation and delayed detached-eddy simulation have received considerable attentions, especially in engineering practices. Because they combine the advantages of RANS and LES methods, the pure LES and DNS are costly to be applied to high Reynolds number flows due to unacceptable CPU requirements [1].
The systematic research on vortex shedding started from the explanation of von Kármán, and the phenomenon is since associated to his name. The alternate vortex shedding is considered as responsible for a remarkable number of collapses of civil structures and damage to industrial equipment. Lienhard et al. analyzed a large number of experimental results of circular cylinders and summarized the law of vortex shedding at different Re [2]. An important feature of the flow past a circular cylinder is that the flow pattern depends on the Reynolds number. The knowledge of the mechanism of cylinder wake improved the understanding of turbine wake. Gehrer et al. [3] carried out experimental research by means of laser Doppler velocimetry (LDV) on the wake vortex of a linear turbine cascade and simulated it with different unsteady Reynolds averaged N-S simulation (URANS) models. As a result, the URANS method had insufficient resolution for unsteady wake vortices. Experiments of a linear turbine cascade by Sieverding et al. [4] revealed that the distances between the adjacent wake vortices are not equal, which is quite different from the flow past a circular cylinder. This is because the flow conditions are the same on both sides of a cylinder. By using the DDES method based on the S-A one-equation model [5], Lin et al. carried out a detailed numerical study on an HPT guide vane [6]. The periodic wake vortex shedding is an important origin of unsteadiness, but the unsteady effects related to the wake vortex have not been fully understood, because of the high Reynolds number, transonic Mach numbers, and high frequency, which impose great challenges for both experimental facilities and CFD solvers [6]. Their numerical results confirmed the findings of Sieverding et al., that is, the distance between an upstream pressure-side shedding vortex and a suction-side shedding vortex is always smaller than the distance between the latter and the downstream pressure-side shedding vortex. Egorov and Menter [7,8] applied the SST-SAS model to a simplified blade geometry of MTU aeroengines. The strong mixing process of gas and cold jets was correctly predicted, and the film cooling effectiveness of the SAS was more consistent with the experiments, compared with URANS. Morata et al. [9] compared the performance of the LES approach in HPT using structured and unstructured grids, respectively. The LES results can be used as an assessment criterion for the global hybrid model.
A modified technique based on the Menter SST-SAS model was proposed in the previous work [10] and improved the spatial resolution of the original model in turbomachinery. As a continuation of the above work, in this paper, an in-house CFD solver integrated with the improved SAS model is employed to delineate the unsteady phenomena and local loss related to vortex shedding, and the underlying flow mechanisms of transonic HPT vanes are analyzed in different working conditions.

SST-SAS Model Description.
In the course of the DESIDER project, Egorov and Menter [11] proposed the SST-SAS model and implemented in the commercial ANSYS CFX solver. The governing equations of the SAS model differ from those of the SST RANS model by an additional source term Q SAS in the transport equation for the turbulence eddy frequency ω. The governing equations can be expressed as The key to the SAS method is the introduction of the von Kármán length scale, L vK , into the scale-determining equation of the RANS turbulence model. The second length scale L vK depends on the local flow topology, which has different responses in steady flow regions and unsteady flow regions. The von Kármán scale is defined as the ratio of the first velocity derivative divided by the second velocity derivative, and the definition given by Menter is as follows: where U i represents the velocity vector, and the average strain rate S replaces the first velocity derivative U ′ . But when the flow unsteadiness is "insufficient," the scale-resolving capability of the original SST-SAS model may not be fully activated [12]. In order to solve the response problem to unsteadiness, a redefined von Kármán scale is proposed by theoretical analysis. According to the Helmholtz velocity decomposition theorem [13], in a three-dimensional vector field, U ′ can be resolved into a sum of an irrotational vector and a solenoidal vector. The average strain rate S represents the tensile and shear deformation, and the average vorticity Ω represents the rotation. Thus, the mathematical formulas of U ′ should contain both effects, and the Decaying Homogenous Isotropic Turbulence (DHIT) case [14] is used to quantitatively determine its two constants, C 1 and C 2 .

2
International Journal of Aerospace Engineering The complete equation for L vK in the present SAS can be written as where V CV is the control volume size.   [15]. A large number of experimental and numerical studies have been carried out on the transonic turbine guide vane, and its design parameters are given in Table 1. In the present study, the computational domain and mesh details at the leading and the trailing edges are shown in Figure 1. The total mesh cells are 7 million, which are divided into 76 blocks for parallel computation. The mesh has 60 evenly spaced points in the spanwise direction with a span length equal to 14.8% of the chord length, and the first layer near wall satisfies the simulation requirement of y + less than 1.0.
2.2.2. Case 2: CNRC Cascade. The second research object is from the Gas Turbine Laboratory of the National Research Council of Canada (CNRC) [16], and the parameters of the CNRC cascade are given in Table 1. Since there is no corresponding experimental schlieren in Case 1, the experimental schlieren images are used for quantitative verification. The CNRC cascade consisted of a set of 6 blades, whose profile represented the midspan section of a high-pressure turbine vane originally designed for a highly-loaded, single stage turbine with a stage pressure ratio of 3.8. The experimental tests were carried out in a large-scale, continuous flow, transonic linear cascade wind tunnel, and the goal of the CNRC project was to study the important flow phenomena that affect the performance of turbines, such as transition, wake and calmed region, shock/boundary layer interaction, and vortex shedding.
The computational domain is shown in Figure 2. Because of the periodicity in the span, in order to save computing resources, a span length equal to 4.9% of its chord length is selected for simulations. The total mesh cell size is 1.45 million, where the mesh has 22 evenly spaced points in the spanwise direction. In the first layer near wall, the maximum y + is less than one.

Computational Setup.
For Case 1, based on the SST turbulence model, the RANS, the URANS, Menter SST-SAS, Menter SST-DDES, and the present SAS were used for simulations with the same structured mesh. Except for the RANS which adopted the WENO third-order scheme, because the difference between the third-order and the fifth-order scheme is very small for RANS, all other simulations used the WENO fifth-order scheme to reconstruct the inviscid flux at interfaces. The sixth-order central difference scheme was used to calculate the viscous flux. The total temperature, total pressure, and velocity direction were given at the inlet, with the mean static pressure defined at the outlet based on the flow conditions for the VKI experiments as in Table 2.
The blade wall was set as adiabatic wall with a no-slip condition. Periodic boundary condition was used in the circumferential direction and the spanwise direction. The dual-time step method was used for unsteady calculations, and the time step size Δt is 0.25 μs. For Case 2, the numerical method and the boundary condition were given in the same way as in the first case. The isentropic exit Mach number under different operating conditions was 0.8, 1.0, and 1.16, respectively, which was obtained by adjusting the back pressure.  Figure 3 shows the comparisons of isentropic Mach number distributions  Table 2 obtained by the present SAS, the experimental data, and the time-averaged result of LES. The numerical results are basically consistent with the experimental data. In addition, the fluctuation of the velocity at certain downstream locations of the trailing edge is monitored for MUR129, represented by power spectral density (PSD) obtained using a Fourier transform. The PSD curves in Figure 4 follow the -5/3 law, indicating that the inertial region of the turbulence is correctly resolved by the present SAS. Moreover, the peak frequency is 37.5 kHz, which correspond to the vortex shedding frequency, consistent with the shedding frequency of the reference [6]. The grid quantity, the turbulence model, and the numerical scheme adopted in the current case are feasible.

An
Overview of Flow Field Structure. The density gradient contours can clearly delineate the flow fields, and its definition is as follows: Figure 5 shows the dimensionless density gradient for the flow condition MUR235, reference to Table 2, obtained with the RANS, the URANS, Menter SAS, Menter DDES, the present SAS, and the time averaging of the present SAS in sequence, followed by the LES result cited from the reference. The RANS model cannot predict the existence of the pressure wave and the vortex shedding at all. The trailing edge pressure wave calculated by the URANS is weak, and the shape of wake vortex is quite different from the actual flow structure. There is a lack of resolution in the wake vortex for Menter SAS, indicating that the ability of scale resolving is not fully excited. As shown in Figures 5(d), 5(e), and 5(g), the result of the present SAS is as capable as the large eddy simulation, and the DDES model also performs well: (1) there is a weak shock wave ① at approximately 84% axial chord length of the suction surface; (2) the wake vortices are predicted in pairs; (3) with the shedding of the wake, the pressure waves are induced, and the pressure wave ④ formed on the pressure side reflects after hitting the   International Journal of Aerospace Engineering suction surface of the adjacent vane; (4) there are interactions among the shock wave, the pressure waves, and the wake. In particular, the suction-side pressure wave ③ dissipate the wake flow, and the core of wake vortices basically disappears, which qualitatively shows the broken process of wake vortex during transport. Therefore, the following sections use the results of the modified SAS model for flow mechanism analysis.
The contours in Figure 6 indicate the similar flow details for the flow condition MUR129. The difference from MUR235 is that there is no shock wave in this flow condition. As the exit Mach number decreases, the dissipation of pressure waves on the wake becomes weaker, so the core of wake vortices can maintain its stability and be transported further downstream. show the scale characteristics during vortex transportation. The vortices shedding from the pressure surface and the suction surface are denoted as VPS (pressure-side shedding vortex) and VSS (suction-side shedding vortex), respectively, and are numbered in chronological order. The vortex core is defined as the position with the lowest pressure. In each vortex pair, the intensity of VPS is stronger than that of VSS, and the streamwise distance between VPS 1 and VSS 1 , d 1 , is smaller than the streamwise distance between VSS 1 and VPS 2 , d 2 , which are quite different from the vortices of a circular cylinder. So when studying the effect of upstream wake vortices on the current blade row, the feasibility of using circular cylinders to generate upstream wakes is worth reexamination.
To explain the asymmetry phenomena, it is necessary to start from the generation process of wake vortex. The coherent structures displayed by Q-criterion at different moments are shown in Figure 8, and the pressure coefficient is used for   International Journal of Aerospace Engineering coloring. First, define the period of a vortex pair as T wake , and define the moment when the suction-side near-wall vortex starts to appear as t = 0 T wake . Then, the vortex continuously absorbs the low-momentum fluid near the trailing edge and grows over time as shown in Figure 8(a). At the same time, due to the blocking effect of VSS 1 on the main flow, a highpressure zone appears near the trailing edge. The wake is a velocity deficit, especially just downstream of the trailing edge, and the downstream VSS 1 enters the range of the main flow on the suction side. But the main flow maintains a high speed due to inertia, thus forming a high-pressure region similar to the stagnation flow. Therefore, due to the pressure difference and the blade geometry, the near-wall vortex is divided into two parts, denoted as M and N, respectively. Soon at 0.53 T wake in Figure 8(b), the vortex N is divided into two parts by the same pressure difference, denoted as N 1 and N 2 , respectively. Then, the downstream VPS 1 gradually enters the main flow on the pressure side in Figure 8(c). As expected, these two pressure differences have the same squeezing effect on N 1 and N 2 , and at about 0.68 T wake , they are divided separately, namely, N 11 , N 12 , N 21 , and N 22 . Figure 8(d) indicates that the N 22 becomes a part of VPS 1 and moves downstream following the previous vortex pair; in addition, the N 12 is related to the formation of the suction-side pressure wave. Eventually, the N 11 and N 21 vortices merge to become the main part of the new VSS, and the vortex M is gradually dissipated. It is worth noting that VPS does not undergo such a complicated splitting process, so its strength is stronger than VSS. For the same reason, the asymmetric flow mechanism also causes d 1 to be smaller than d 2 .
The pressure coefficient is calculated by  As the vortices transported to the downstream, the distance d 1 gradually decreases. On the one hand, this is because the reverse pressure gradient received by VSS is larger than VPS as shown by the black dotted line in Figure 9, so that the motion of VSS is hindered and decelerated. On the other hand, the double-vortex structure formed in the vortex pair makes the two vortices closer to each other as shown in Figure 10. Ignoring the effect of the main flow, the spatial distance between the two vortices in the same vortex pair is small, so there is an interaction between them. There is an induced rotating flow field around each concentrated vortex, which means that VPS will produce an induced velocity for VSS, namely, U ab in Figure 10, and vice versa. But the intensity of VPS is stronger, that is, the induced velocity U ab is also greater than U ba . Therefore, VSS moves upstream relative to VPS, and the distance between the two vortices in the main flow direction is further reduced. Figure 11 shows the formation and development of suctionside pressure waves which are generated by the intense vortex shedding. More specifically, it should be that the vortex N 12 induces the generation of the suction-side pressure wave. As mentioned above, the vortex N 1 is divided into two parts at 0.68 T wake , namely, N 11 and N 12 , and the threedimensional view at this instant is displayed in Figure 11(a). As the vortex N 12 forms, it gradually moves away from the trailing edge and enters the main flow in a direction almost perpendicular to the main flow. The blockage effect on the main flow will cause pressure fluctuations to propagate in all directions, since the freestream is a subsonic flow, but the waves travel at the local speed of sound, which means that the Mach number relative to the motion of pressure waves is 1.0. It is worth noting that the main flow has absolute velocity and direction, and its Mach number near the trailing edge is around 0.9, so the pressure waves observed from Figures 11(b) and 11(c) have different moving speeds in different directions. These pressure waves move upstream over the suction surface at about 11% of the free-stream velocity, affecting the aft suction surface. In addition, since the suction-side pressure wave is closely related to the vortex N 12 , its frequency is exactly equal to the frequency of the VSS. In the same way, it can be concluded that the formation of pressure-side pressure waves is closely related to VPS.

Unsteady
Interactions Related to Wake Vortex. Combined with the instantaneous three-dimensional perspective of Figure 12, the evolution of wake vortex is explored. The shedding vortex just leaving the trailing edge exhibits periodicity in the spanwise direction, and the shape of the vortex core is similar to that of a cylinder with two-dimensional features. The tube around the vortex core exhibits threedimensional features but maintains periodicity in the span. When the wake vortex moves to the area affected by the shock wave, the spatial size of the vortex pair increases, and the presence of the shock causes the vortex tube to appear three-dimensional rolled up along the span, marked with the white dotted line. As a comparison in Figure 12(b), there is almost no radial mixing process in the same region under the flow condition MUR129. Then, the vortex is transported to the influence area of suction-side pressure waves. As the pressure wave passes, accompanied by continuous violent oscillations, the large-scale tube is strongly dissipated and broken into multiscale structures, regardless of MUR129 or MUR235. It can therefore be concluded that the interactions of the shock wave, pressure wave, and wake vortex accelerate the instability and fragmentation of the vortex.
3.1.6. Time-Averaged Flow Analysis. In order to provide a design reference for the HPT vane in the existing RANS design system, it is necessary to quantify the impact of unsteady flow on the aerodynamic performance. First, the time-averaged flow field of the present SAS under the flow condition MUR129 is compared with the RANS result as shown in Figure 13, and the centerline direction of the    Figures 5(a) and 5(f). It can be concluded that the main difference between the unsteady time-averaged flow and the RANS flow is in the wake region. On the one hand, the width of the unsteady time-averaged wake is wider, indicating that the mixing with the main flow is stronger, and on the other hand, the unsteady time-averaged wake has obvious asymmetry about the centerline. The trailing edge and the free shear layers on both sides of the trailing edge form a triangular area, which is called the base region, and the pressure in the base region is closely related to the wake vortices; increasing the base pressure coefficient can effectively reduce the trailing edge loss [17]. Therefore, in the design process of turbine blades, correct estimation of the base pressure coefficient is essential for predicting profile loss. Figure 14 shows the base pressure coefficient distribution obtained by the RANS and the time averaging, where the horizontal axis represents the arc length from the blade surface to the center point of the trailing edge (dimensionless by the trailing edge diameter). There are three main differences.      [16,20]. Firstly, under different flow conditions, the RANS results are always higher than the time-averaged results, which can underestimate the importance of wake loss. Secondly, the RANS have a base pressure plateau, but the timeaveraged results show significant nonuniformity. Thirdly, the RANS results are symmetrical about the trailing edge centerline, whereas the time-averaged results are not. The profile loss can be quantified based on the stagnation pressure loss coefficient. For the flow condition MUR235, the loss coefficients of the RANS and the present SAS are 5.86% and 6.32%, respectively, and the relative error is 7.85%. For another flow condition MUR129, the RANS and the unsteady time-averaged loss coefficients are 4.50% and 4.89%, respectively, and their relative error is 8.67%. Combined with the analysis in Section 3.1.7, it can be concluded that the error mainly comes from the shedding and development of wake vortex, because the RANS system overestimates the base pressure coefficient and ignores the nonuniformity and asymmetry of the base pressure. The stagnation pressure loss coefficient is calculated by 3.1.7. Local Loss Analysis. There are many definitions of loss coefficient in regular use for blade rows, such as 14 International Journal of Aerospace Engineering stagnation pressure loss coefficient, enthalpy loss coefficient, and entropy loss coefficient [17]. But those coefficients are used for global loss analysis, which represent the total losses of the whole region. The classical breakdown of loss into profile loss, end-wall loss, and leakage loss continues to be widely used in advanced design, so it is necessary to reflect the information about local loss generation in high-fidelity simulation. Essentially, losses are the results from viscous irreversibility and thermal irreversibility, and the entropy generation rate can quantify irreversibility and reveal the specific location of local loss [18,19]. The definition of entropy generation rate is Figures 15(a) and 15(b) are the local entropy generation rates corresponding to viscous irreversibility and thermal irre-versibility of the flow condition MUR235, respectively. Similarly, Figures 16(a) and 16(b) show the viscous entropy generation rate and the thermal entropy generation rate of MUR129, respectively. For linear cascades, the only loss is profile loss, and its components come from wake vortex, boundary layer, pressure wave, shock wave, etc. In particular, the periodic wake vortex is the most important origin of profile loss, and the pressure waves induced by wake vortex strongly dissipate the downstream vortices. In addition, these graphs visually compare viscous loss and thermal loss, the former is significantly larger than the latter. In order to quantitatively analyze these two losses, the Bejan number (Be) is introduced here, which represents the percentage of thermal loss in total local losses [19]. The distribution of Be in Figures 17 and 18illustrates that in areas with high local losses, such as wake vortex and boundary layer, the Be is between 0.1% and 1%, while in other low local-loss areas, the Be is about 90%. It can be summarized that the viscous effect plays a leading role in the loss of the VKI LS89 cascade.

Energy Separation in Wakes.
In addition to the entropy loss, vortex shedding may cause adverse effects which include high-frequency sound propagation, vibration, and local high heat transfer and energy separation into hot and cold regions [20]. The stagnation temperature on the wake centerline is found to be 30 K lower than that of the incoming fluid, while the stagnation temperature at the edges of the wake can be up to 30 K higher than incoming fluid, as shown in Figure 19(a), and Figure 19(b) is the stagnation temperature of the A-A profile. It was demonstrated that the temperature separation is a thermoacoustic effect where the wake centerline emerges colder than the surrounding fluid, which leads to the redistribution of the stagnation temperature and the generation of local hot spots. The nonuniform total temperature distribution is a source of entropy production, and more importantly, it should be carefully considered in the turbine cooling design.  [16,20] are selected to quantitatively analyze the relative error of simulation results and experimental data in different flow conditions. The density gradient contours are displayed in Figures 20-22.
For the supersonic condition in Figure 20, the exit Mach number is equal to 1.16, the downstream disturbance cannot propagate upstream, and there is no pressure wave in the turbine nozzle passage. Instead, there are two strong oblique shock waves generated at the trailing edge, and the passage shock wave impinges on the suction surface of the adjacent blade. The incident point of passage shock is well predicted, and the relative error of its length and angle characteristics is within 2%. After experiencing the shock/boundary-layer interaction, a reflected shock wave is formed. Under the supersonic flow condition, the origin of the vortex shedding moved from the vane trailing edge to the confluence of the two trailing edge shear layers, and the predicted vortex shedding frequency is about 25.9 kHz, which is consistent with the experimental frequency in the literature [16].
For Ma ex = 1:0 in reference [16] (or Ma ex = 0:97 in reference [20]), the experimental flow structure is shown in Figure 21(b). There is a normal shock wave at the throat, and supersonic flow decelerates to subsonic flow through this shock wave. The pressure waves and their reflected waves are located downstream of the normal shock wave. The present SAS model reproduces the flow details well, such as the distribution of wake vortex and wave system. The distances between different vortices are marked in the figures, that is, d 1 , d 2 , and d 3 in Figure 21(a) and l 1 , l 2 , and l 3 in Figure 21(b), and the maximum error between the distances is about 2%. For Ma ex = 0:8 in Figure 22, the flow structure similar to the flow condition MUR129 of LS89 can be observed. The wake vortices on the SS and PS appear in pairs and shed off alternately. With the shedding of the wake vortices, the pressure waves on both sides are induced and they propagate upstream, and the pressure-side pressure waves reflect on the adjacent vane.
The power spectral density curves in Figure 23 reveal that the turbulent inertial regions are correctly resolved. All the above results illustrate that the numerical methods can reproduce the flow field of real experimental flow conditions and it is reliable for analyzing the flow mechanisms inside transonic HPT.

Local Loss Analysis.
In this section, the local entropy generation under supersonic condition is analyzed. Figures 24(a) and 24(b) are the entropy generation rates due to viscous effect and thermal effect, respectively, and Figure 24(c) shows the corresponding local Bejan number distribution. It can be found that although the Ma ex is greater than 1, the conclusions are consistent with that obtained from the subsonic condition. First, the entropy generation rates are on the same order of magnitude as in Case 1, and then, the Be is between 0.1% and 1% in wakes and boundary layer, while Be is close to 90% in other areas.

Conclusions
A modified SAS model based on the Menter SST-SAS model was applied to the numerical simulation and flow mechanism analysis of the transonic turbine. The conclusions are as follows: (1) The wake vortex-related flow is the main origin of unsteadiness and losses in transonic HPT cascade, and local entropy generation analysis shows that the viscous loss is dominant (2) Comparing the RANS results with the time-averaged results of the present SAS, it can be found that the RANS overestimates the base pressure coefficient and cannot predict the nonuniformity and asymmetry of the base pressure, which leads to errors in the prediction of profile loss (3) The formation of suction-side shedding vortex has gone through a complex process of division and reunion, whereas the pressure-side wake vortex has not, due to the asymmetric geometries and pressure conditions. The vortices VPS and N 12 induced the generation of pressure-side and suction-side pressure waves, respectively (4) There is energy separation in the wake vortex, and the stagnation temperature difference between the hot region and the cold region is as high as 15% (5) The high-fidelity simulation based on the present SAS model can reveal the length scale and time scale characteristics of the high-pressure turbine. The CPU hours of the present SAS is about 2.5% of the structured LES in the reference Future works include continued flow mechanisms in other turbomachinery applications, in particular for multistage turbomachinery design and aerothermal performance analysis. 16 International Journal of Aerospace Engineering