Investigation of Unsteady Flow Interactions in a Transonic High Pressure Turbine Using Nonlinear Harmonic Method

The performance of a transonic high pressure turbine is mainly influenced by the unsteady interactions associated with the passing blades. In this paper, the unsteady flow interactions in a transonic turbine have been numerically investigated using the nonlinear harmonic (NLH) method in comparison with the steady and unsteady Reynolds-averaged Navier–Stokes (RANS). The comparison shows that the NLH method using three harmonics could capture the main unsteady flow interactions efficiently with about seven times smaller computational cost than the unsteady RANS, resulting in a more accurate time-averaged flow than for steady RANS. However, the continuity of the flow variables across the rotor-stator interface has shown some discrepancies compared with the unsteady RANS, which can be further satisfied by increasing the numbers of harmonics. The unsteady interactions are analyzed in detail; the results show that the wake and trailing edge shock from the upstream stator are the major sources of unsteadiness in the downstream rotor passage. The stator trailing edge shock impinges on the suction side of the passing rotor blades and generates pressure waves. These pressure waves are periodically reflected back to trigger the stator wake shedding. These waves are strong enough to travel through the rotor passage, and eventually affect the flow at the rotor’s trailing edge. The stator wakes are chopped by the downstream rotor, and travel through the rotor passage. This significantly enhances the unsteadiness of the flow near the rotor trailing edge. Lastly, the deterministic stresses and enthalpy distributions extracted from the NLH method have revealed that the effects of the unsteadiness are relatively weaker in the axial direction. Furthermore, the deterministic correlations analysis has shown that, some empirical deterministic correlations models based on the decay concept of compressors are not suitable for turbines.


Introduction
The complex three-dimensional flows in high pressure turbines are mainly attributed to the inherent unsteadiness due to the relative motion between the stator and rotor blades. The understanding of these unsteady flows is essential to further reduce engine weight. Weight reducation can be achieved by using fewer turbine stages or fewer blades in a row. All these measures would lead to increased blade loading. A transonic high pressure turbine is a crucial component of a modern aero-engine. With fewer stages, transonic turbines operate at higher pressure ratios, The previous studies mainly focus on the validation of the NLH method on various types of compressors while only a few test cases are reported for the turbines [20,23]. Some preliminary work has been done on the Aachen 1.5 stage axial turbine to study clocking effects using the NLH method [20]. However, detailed investigations for the transonic high pressure turbines have been seldom conducted with the NLH method. The NLH method employs deterministic correlations to reproduce the unsteadiness due to the rotor-stator interactions. Modeling deterministic correlations is the key to close the average-passage equation system which could provide a rigorous mathematical framework for accounting for the unsteady blade row interaction in steady state calculation by introducing deterministic correlations. Liu et al. [24,25] studied the deterministic stress and enthalpy correlations in compressors, which could be valuable for modeling deterministic correlations. However, few studies focused on deterministic correlations in turbines. The deterministic analysis is important to evaluate the relative importance of various deterministic components particularly related to the unsteady blade row interactions. The deterministic components of particular interest can be used in the model development of the deterministic based codes for the transonic turbines.
The current study focuses on the unsteady flow interactions in a transonic high-pressure turbine stage using the unsteady time-accurate and the NLH methods. A widely used two-equation shear-stress transport (SST) model has been used for the detailed analysis in all the cases. The results from the NLH are compared with the steady and unsteady RANS. Then the unsteady interactions are analyzed in detail based on the unsteady RANS and the NLH. Lastly, the deterministic stresses and enthalpy distributions extracted from the NLH method are investigated in both the rotor and stator passages. The results can help designers to understand the physics of unsteady losses and their origins, and model deterministic correlations.

Nonlinear Harmonic Method
The NLH method is built around the idea that the instantaneous flow in a blade row is decomposed into a time-averaged part and a sum of periodic unsteady perturbations that are Fourier decomposed into N harmonics [18]: where U( → r ) is the time-averaged part and U ( → r , t) is the unsteady part due to unsteady disturbances. In the current study, the unsteady disturbances are generated by the passing blades. The periodic perturbation term (U ( → r , t)) is transformed into the frequency domain by the Fourier transform. The solution accuracy is ensured by the number of considered harmonics and the number of perturbations used in the decomposition [26].
The time-mean flow equation is obtained from time-averaging the unsteady equations. The conservative form of the time-mean flow equations are formulated from the finite-volume method as: Here, i represents the cell number of volume Ω i . F c and F v are discretized convective and viscous fluxes, respectively. S T , is the source term containing the Coriolis and centrifugal terms.
The complex amplitudes U k are made pseudo time dependent by casting the first order linearized equations into the frequency domain. The conservative form of these amplitudes is formulated in the finite-volume method as: Like Reynolds averaging, the NLH method includes additional stresses, with the form of time-averaged products of fluctuations. These terms are called deterministic stresses. The deterministic stresses represent the unsteadiness in the flow field and these stresses are directly computed from the harmonic values to achieve the closure of the system. The NLH method works in the frequency domain. That allows the calculation of only one blade passage per row. If the physical phenomenon includes only a few harmonics, the NLH method can be significantly cheaper than the conventional unsteady methods [27].

Test Case
The meridional section of the transonic turbine stage is shown in Figure 1a. The convergent-divergent vane accelerates the flow to supersonic speed. The hub line is cylindrical while the shroud line is conical in the rotor section. The minimum area of the flow path is near the trailing edge of the stator. The number of stator vanes and rotor blades are 24 and 36, respectively. The rotation speed is 10,600 revolutions per minute (RPM) and the pressure ratio is 3.5. The Reynolds number at stage inlet is about 5.25 × 10 6 based on the stator chord. The geometric turning angle for the stator vane is 70 • and for the rotor blade 107 • . The absolute Mach number at stator vane inlet is about 0.14 and at rotor inlet is about 1.12. The rotor tip clearance is 1 mm. More details could be found in [14]. As shown in Figure 1b, two planes are selected downstream of the stator and are named DS1 and DS2 for analysis of results. Similarly two planes are selected downstream of the rotor, named DR1 and DR2.
The complex amplitudes are made pseudo time dependent by casting the first order linearized equations into the frequency domain. The conservative form of these amplitudes is formulated in the finite-volume method as: Like Reynolds averaging, the NLH method includes additional stresses, with the form of timeaveraged products of fluctuations. These terms are called deterministic stresses. The deterministic stresses represent the unsteadiness in the flow field and these stresses are directly computed from the harmonic values to achieve the closure of the system. The NLH method works in the frequency domain. That allows the calculation of only one blade passage per row. If the physical phenomenon includes only a few harmonics, the NLH method can be significantly cheaper than the conventional unsteady methods [27].

Test Case
The meridional section of the transonic turbine stage is shown in Figure 1a. The convergentdivergent vane accelerates the flow to supersonic speed. The hub line is cylindrical while the shroud line is conical in the rotor section. The minimum area of the flow path is near the trailing edge of the stator. The number of stator vanes and rotor blades are 24 and 36, respectively. The rotation speed is 10,600 revolutions per minute (RPM) and the pressure ratio is 3.5. The Reynolds number at stage inlet is about 5.25 × 10 6 based on the stator chord. The geometric turning angle for the stator vane is 70° and for the rotor blade 107°. The absolute Mach number at stator vane inlet is about 0.14 and at rotor inlet is about 1.12. The rotor tip clearance is 1mm. More details could be found in [14]. As shown in Figure 1b, two planes are selected downstream of the stator and are named DS1 and DS2 for analysis of results. Similarly two planes are selected downstream of the rotor, named DR1 and DR2.

Numerical Method
The NUMECA Fine/Turbo CFD flow solver (version 10.2, NUMECA International, Brussels, Belgium) is used in the current study. It is a three-dimensional multi-block flow solver using the steady, unsteady and harmonic methods. This solver has wide applications in turbomachinery simulation and design. The flow solver is employed with the Navier-Stokes code using finite volume method. The second order central scheme with a second and fourth order artificial dissipation is used for spatial discretization. For temporal discretization, a fourth order Runge-Kutta scheme is employed. The local time-stepping approach and implicit residual smoothing are used to boost the convergence time.
For the NLH method, the number of frequencies per perturbation is set to 3, and the maximum number of perturbations per blade row is set to 1. Since the harmonic method uses a limited number

Numerical Method
The NUMECA Fine/Turbo CFD flow solver (version 10.2, NUMECA International, Brussels, Belgium) is used in the current study. It is a three-dimensional multi-block flow solver using the steady, unsteady and harmonic methods. This solver has wide applications in turbomachinery simulation and design. The flow solver is employed with the Navier-Stokes code using finite volume method. The second order central scheme with a second and fourth order artificial dissipation is used for spatial discretization. For temporal discretization, a fourth order Runge-Kutta scheme is employed. The local time-stepping approach and implicit residual smoothing are used to boost the convergence time.
For the NLH method, the number of frequencies per perturbation is set to 3, and the maximum number of perturbations per blade row is set to 1. Since the harmonic method uses a limited number of harmonics, thus continuity cannot be rigorously satisfied, but it is observed that the continuity increases with more harmonics considered in the calculation. For the unsteady time-accurate simulation, the computational domain consists of two stator vane passages and three rotor blade passages. This avoids the introduction of an error due to the geometry scaling. The dual time stepping approach proposed by Jameson is used [28]. The numbers of physical time steps are set to 120 for each cycle after careful monitoring of the convergence of some important parameters (mass flow, efficiency and pressure ratio). In order to ensure convergence, the difference between the inlet and outlet mass flow rate was ensured to be <0.1% in all simulations. The time step for the unsteady method is about 3.93 × 10 −6 s and the same time step has been used to construct the computed harmonics in the time domain for the NLH method. This makes it possible to compare the instantaneous data between the two methods at the same time steps. The time step size independence study was initially performed, based on that study, the above mentioned time step size was selected which provides a 0.25 • angular resolution between each time step. This angular resolution is fine enough to capture the unsteady flow effects.
The selection of an appropriate turbulence model is critical to accurately predict the flow in turbomachinery. There is no single universal turbulence model available for all engineering applications, and the existing turbulence models do not produce results with acceptable accuracy in turbomachinery applications [29]. The SST model was the most appropriate for the current problem since it effectively handles the separations and adverse pressure gradients. The SST model modifies the turbulence viscosity function to improve the prediction of separated flows and to avoid overestimating Reynolds stresses with the k-ω and k-ε models in adverse pressure gradients. This modification is done based on the Townsend-Lighthill design for the specification of turbulent eddy viscosity [19]. Additionally, the multigrid method was used to accelerate the convergence. The computational time required by the NLH method is about five times longer than that of the steady-state method and around seven times shorter than that of the unsteady method in the present investigation.

Numerical Grid
In the current investigations, the computational grid is generated using the IGG-AutoGrid5, which is a fully automated grid generation tool. Since the grid resolution in radial and circumferential directions has a significant impact on the result, the fully structured multi-block strategy is employed to ensure the optimal grid.
For each blade passage, an O4H-type grid is generated and a butterfly topology, which consists of two blocks with an inner H block and an external O block, is used for the grid in the rotor tip gap. The number of grid points in the circumferential direction must follow the guideline for the NLH method. At least 30 cells must be placed in the circumferential direction for the first harmonic and no less than 15 cells for the other harmonics [19]. The circumferential cell number is also determined by the number of blades on the each side of the rotor-stator interface. The recommended grid resolution in the circumferential direction is given by following equations [19]: where, N is the number of harmonics, N s is the number of stator blades and N r is the number of rotor blades. A grid independence study has been conducted with grid sizes of about 1.8, 2.9 and 3.8 million for a single flow channel and finally a 2.9 million grid size is selected that also satisfies the above mentioned grid requirements for the NLH method. The same grid distribution has been subsequently used for the steady and time-domain unsteady simulation to ensure the grid consistency. The domain scaling method for the R/S interface proposed by Rai [19] has been implemented for unsteady calculations. For this method, the pitch distance must be identical across the R/S interface. The main advantage of having the same pitch distance across the R/S interface is that it avoids the need to consider any time periodicity in the boundary condition To resolve the boundary layer on the blade and over the endwall, the y+ of the first cell away from the wall is less than 3. For most grid cells, the skewness angle is larger than 36 degrees; the aspect ratio and expansion ratio are less than 1000 and 1.8 respectively. The grid distributions near the leading edge of both the stator and the rotor are shown in Figure 2. distance must be identical across the R/S interface. The main advantage of having the same pitch distance across the R/S interface is that it avoids the need to consider any time periodicity in the boundary condition treatment. There are 24 stator vanes and 36 rotor blades in the stage. The computational domain consists of two stator vanes and three rotor blades for unsteady time-accurate simulation.
To resolve the boundary layer on the blade and over the endwall, the y+ of the first cell away from the wall is less than 3. For most grid cells, the skewness angle is larger than 36 degrees; the aspect ratio and expansion ratio are less than 1000 and 1.8 respectively. The grid distributions near the leading edge of both the stator and the rotor are shown in Figure 2.

Boundary Conditions
The absolute total pressure (3.5 bar) and total temperature (403 K) on the inlet boundary are specified. An average static pressure is imposed at the outlet boundary. The other dependent variables on the outlet are obtained from the interior field through extrapolation. The R/S interface in the current case is very close to the blades therefore, non-physical reflections of pressure waves, especially due to the shock, can occur at the interface. To handle this issue, a two dimensional nonreflecting boundary condition based on a characteristic analysis of the linearized Euler equations was applied at the R/S interface for the steady and NLH calculations [6]. The result of the steady calculation is used as the initial solution for both the unsteady and the NLH simulations.

Validation of Numerical Method
The Aachen 1.5 stage axial turbine test case provided by the RWTH Aachen University was employed for the validation of the numerical schemes used in the current study. The Aachen 1.5 stage axial turbine was tested at the Institute of Jet Propulsion and Turbomachinery. For validation, only the first stage was considered that is a subsonic axial turbine with untwisted blades. The first stage has 36 stator and 41 rotor blades. The turbine was operated at 3500 RPM. The turbine inlet temperature was kept at 308 K during experiments. The stator inlet Reynolds number was 6.8 × 10 5 and the inlet flow velocity was 45 m/s. The rotor tip clearance was 0.4 mm for this particular case. The detailed experimental procedure and data can be found in references [30][31][32].
The SST turbulence model has been employed for the validation case. The steady simulation was initially performed with the mixing plane approach. The steady results were subsequently used for the initialization of the unsteady RANS and NLH simulations. The original blade count ratio is 36:41; however to reduce the computational cost of the unsteady simulation, the blade count ratio of 6:7 has been used. The blade scaling factor of 0.976 was used for the rotor to offer the same blockage. The accuracy of unsteady simulation can be further increased by considering the full stage simulation but at the expense of much higher computational cost. The NLH method works in the frequency domain, therefore a single blade passage is required for the simulation. Three harmonics have been considered and the maximum number of perturbations has been set to 1 for the NLH simulation.
The rotor surface pressure distribution is presented in Figure 3. The time-averaged results from the unsteady RANS and NLH method are presented in comparison with the experimental data. The unsteady RANS and NLH results are in good agreement with the experimental data. The difference in predicted results on the rotor pressure side negligible, however the difference in results on the

Boundary Conditions
The absolute total pressure (3.5 bar) and total temperature (403 K) on the inlet boundary are specified. An average static pressure is imposed at the outlet boundary. The other dependent variables on the outlet are obtained from the interior field through extrapolation. The R/S interface in the current case is very close to the blades therefore, non-physical reflections of pressure waves, especially due to the shock, can occur at the interface. To handle this issue, a two dimensional non-reflecting boundary condition based on a characteristic analysis of the linearized Euler equations was applied at the R/S interface for the steady and NLH calculations [6]. The result of the steady calculation is used as the initial solution for both the unsteady and the NLH simulations.

Validation of Numerical Method
The Aachen 1.5 stage axial turbine test case provided by the RWTH Aachen University was employed for the validation of the numerical schemes used in the current study. The Aachen 1.5 stage axial turbine was tested at the Institute of Jet Propulsion and Turbomachinery. For validation, only the first stage was considered that is a subsonic axial turbine with untwisted blades. The first stage has 36 stator and 41 rotor blades. The turbine was operated at 3500 RPM. The turbine inlet temperature was kept at 308 K during experiments. The stator inlet Reynolds number was 6.8 × 10 5 and the inlet flow velocity was 45 m/s. The rotor tip clearance was 0.4 mm for this particular case. The detailed experimental procedure and data can be found in references [30][31][32].
The SST turbulence model has been employed for the validation case. The steady simulation was initially performed with the mixing plane approach. The steady results were subsequently used for the initialization of the unsteady RANS and NLH simulations. The original blade count ratio is 36:41; however to reduce the computational cost of the unsteady simulation, the blade count ratio of 6:7 has been used. The blade scaling factor of 0.976 was used for the rotor to offer the same blockage. The accuracy of unsteady simulation can be further increased by considering the full stage simulation but at the expense of much higher computational cost. The NLH method works in the frequency domain, therefore a single blade passage is required for the simulation. Three harmonics have been considered and the maximum number of perturbations has been set to 1 for the NLH simulation.
The rotor surface pressure distribution is presented in Figure 3. The time-averaged results from the unsteady RANS and NLH method are presented in comparison with the experimental data. The unsteady RANS and NLH results are in good agreement with the experimental data. The difference in predicted results on the rotor pressure side negligible, however the difference in results on the rotor pressure side is within 2.5%. The above test case shows that the NLH method has predicted the promising results in comparison with unsteady RANS and the experimental data. The NLH method has an added advantage over the unsteady RANS of considerably reduced computational cost. rotor pressure side is within 2.5%. The above test case shows that the NLH method has predicted the promising results in comparison with unsteady RANS and the experimental data. The NLH method has an added advantage over the unsteady RANS of considerably reduced computational cost.  Figure 4 presents the performance of the turbine stage predicted by the steady, unsteady and NLH simulations at various pressure ratios at about 96% of the designed RPM. No considerable difference can be seen in terms of mass flow rate predicted by the three methods. However, the stage efficiency predicted by the NLH method is closer to the unsteady method than the steady method. The same turbulence model and computational domain are used in all methods. The discrepancy is mainly attributed to the unsteady R/S interactions, that are effectively captured by the unsteady and NLH method. The unsteady RANS being the highly accurate method with a much higher computational cost, however the NLH method has predicted the more promising results in comparison with unsteady RANS at a much reduced cost. The isentropic Mach number distribution at the midspan is shown in Figure 5. Generally, the results between the unsteady time-averaged and the NLH method are in good agreement. On the stator blade, no considerable difference between the unsteady time-averaged result and the NLH method can be seen in Mach number distribution. The steady calculation predicts the results on the front-portion of the blade surface correctly, but it deviates at the rear part of the stator blade on the suction side, where the shock wave and unsteady stator-rotor interaction occur. For the rotor blade, some differences on the suction surface in the shock region can be observed in the steady calculation. Again, the discrepancy seen here is mainly attributed to the unsteady R/S interactions. Therefore, a detailed study of the unsteady R/S interactions is necessary to predict the performance accurately.  Figure 4 presents the performance of the turbine stage predicted by the steady, unsteady and NLH simulations at various pressure ratios at about 96% of the designed RPM. No considerable difference can be seen in terms of mass flow rate predicted by the three methods. However, the stage efficiency predicted by the NLH method is closer to the unsteady method than the steady method. The same turbulence model and computational domain are used in all methods. The discrepancy is mainly attributed to the unsteady R/S interactions, that are effectively captured by the unsteady and NLH method. The unsteady RANS being the highly accurate method with a much higher computational cost, however the NLH method has predicted the more promising results in comparison with unsteady RANS at a much reduced cost. rotor pressure side is within 2.5%. The above test case shows that the NLH method has predicted the promising results in comparison with unsteady RANS and the experimental data. The NLH method has an added advantage over the unsteady RANS of considerably reduced computational cost.  Figure 4 presents the performance of the turbine stage predicted by the steady, unsteady and NLH simulations at various pressure ratios at about 96% of the designed RPM. No considerable difference can be seen in terms of mass flow rate predicted by the three methods. However, the stage efficiency predicted by the NLH method is closer to the unsteady method than the steady method. The same turbulence model and computational domain are used in all methods. The discrepancy is mainly attributed to the unsteady R/S interactions, that are effectively captured by the unsteady and NLH method. The unsteady RANS being the highly accurate method with a much higher computational cost, however the NLH method has predicted the more promising results in comparison with unsteady RANS at a much reduced cost. The isentropic Mach number distribution at the midspan is shown in Figure 5. Generally, the results between the unsteady time-averaged and the NLH method are in good agreement. On the stator blade, no considerable difference between the unsteady time-averaged result and the NLH method can be seen in Mach number distribution. The steady calculation predicts the results on the front-portion of the blade surface correctly, but it deviates at the rear part of the stator blade on the suction side, where the shock wave and unsteady stator-rotor interaction occur. For the rotor blade, some differences on the suction surface in the shock region can be observed in the steady calculation. Again, the discrepancy seen here is mainly attributed to the unsteady R/S interactions. Therefore, a detailed study of the unsteady R/S interactions is necessary to predict the performance accurately. The isentropic Mach number distribution at the midspan is shown in Figure 5. Generally, the results between the unsteady time-averaged and the NLH method are in good agreement. On the stator blade, no considerable difference between the unsteady time-averaged result and the NLH method can be seen in Mach number distribution. The steady calculation predicts the results on the front-portion of the blade surface correctly, but it deviates at the rear part of the stator blade on the suction side, where the shock wave and unsteady stator-rotor interaction occur. For the rotor blade, some differences on the suction surface in the shock region can be observed in the steady calculation. Again, the discrepancy seen here is mainly attributed to the unsteady R/S interactions. Therefore, a detailed study of the unsteady R/S interactions is necessary to predict the performance accurately.

Stator Flow Field
The flow field in the stator passage will be discussed in detail in this section. The entropy and the Mach number distribution downstream of the stator at DS2 plane is shown in Figure 6. The stator wake can be clearly identified in the region with high entropy as indicated by the dashed oval. The stator suction side shock moving downstream is obvious in the Mach number distribution near the trailing edge as indicated by the dashed line.  Figure 7 presents the time-averaged velocity, yaw angle and entropy downstream of the stator at two locations (DS1 & DS2) over two stator pitches. Some differences in the stator wake region for the velocity magnitude, yaw angle and entropy are observed between the two methods. There is a significant decrease in velocity at DS1 near stator TE, where the stator wake has a strong influence on the flow field. The yaw angle (tan ) , is being calculated from the meridional (V ) and tangential velocity (V ), and sudden increase in yaw angle is observed across the shock. Entropy distribution shows that the NLH method predicts relatively higher entropy in the shock region for both DS1 and DS2 planes.

Stator Flow Field
The flow field in the stator passage will be discussed in detail in this section. The entropy and the Mach number distribution downstream of the stator at DS2 plane is shown in Figure 6. The stator wake can be clearly identified in the region with high entropy as indicated by the dashed oval. The stator suction side shock moving downstream is obvious in the Mach number distribution near the trailing edge as indicated by the dashed line.

Stator Flow Field
The flow field in the stator passage will be discussed in detail in this section. The entropy and the Mach number distribution downstream of the stator at DS2 plane is shown in Figure 6. The stator wake can be clearly identified in the region with high entropy as indicated by the dashed oval. The stator suction side shock moving downstream is obvious in the Mach number distribution near the trailing edge as indicated by the dashed line.  Figure 7 presents the time-averaged velocity, yaw angle and entropy downstream of the stator at two locations (DS1 & DS2) over two stator pitches. Some differences in the stator wake region for the velocity magnitude, yaw angle and entropy are observed between the two methods. There is a significant decrease in velocity at DS1 near stator TE, where the stator wake has a strong influence on the flow field. The yaw angle (tan ), is being calculated from the meridional (V ) and tangential velocity (V ), and sudden increase in yaw angle is observed across the shock. Entropy distribution shows that the NLH method predicts relatively higher entropy in the shock region for both DS1 and DS2 planes.    Figure 7 presents the time-averaged velocity, yaw angle and entropy downstream of the stator at two locations (DS1 & DS2) over two stator pitches. Some differences in the stator wake region for the velocity magnitude, yaw angle and entropy are observed between the two methods. There is a significant decrease in velocity at DS1 near stator TE, where the stator wake has a strong influence on the flow field. The yaw angle (tan −1 V t V m ), is being calculated from the meridional (V m ) and tangential velocity (V t ), and sudden increase in yaw angle is observed across the shock. Entropy distribution shows that the NLH method predicts relatively higher entropy in the shock region for both DS1 and DS2 planes.

Stator Flow Field
The flow field in the stator passage will be discussed in detail in this section. The entropy and the Mach number distribution downstream of the stator at DS2 plane is shown in Figure 6. The stator wake can be clearly identified in the region with high entropy as indicated by the dashed oval. The stator suction side shock moving downstream is obvious in the Mach number distribution near the trailing edge as indicated by the dashed line.  Figure 7 presents the time-averaged velocity, yaw angle and entropy downstream of the stator at two locations (DS1 & DS2) over two stator pitches. Some differences in the stator wake region for the velocity magnitude, yaw angle and entropy are observed between the two methods. There is a significant decrease in velocity at DS1 near stator TE, where the stator wake has a strong influence on the flow field. The yaw angle (tan ) , is being calculated from the meridional (V ) and tangential velocity (V ), and sudden increase in yaw angle is observed across the shock. Entropy distribution shows that the NLH method predicts relatively higher entropy in the shock region for both DS1 and DS2 planes.      Figure 9 shows the entropy downstream of the rotor at the DR2 plane. The high loss region is located near the shroud, where the flow is highly influenced by the rotor tip clearance. The rotor wake can be identified by the higher entropy region, as indicated by the dashed oval in the figure.

Rotor Flow Field
The upper passage vortex region, combined with the tip leakage vortex accommodates more loss (indicated by 1). The lower passage vortex (indicated by 2) is much smaller than the upper passage vortex. However, it still widens the wake area near the hub. The rotor suction side shock and reflected pressure side shock can be identified in the Mach number distribution. The rotor flow field is strongly influenced by the stator wakes. In the time-resolved velocity flow field, strong variations of the velocity magnitude can be found at the rotor exit. The evolution of these variations is strongly influenced by the highly turbulent stator exit flow. Figure 10 presents the time averaged velocity, yaw angle and entropy at two locations downstream of the rotor (DR1 & DR2) over three blade pitches. The suction side shock and reflected pressure side shock can be seen in the velocity distribution in Figure 10a. No obvious difference can be seen in the results predicted by the unsteady and NLH methods downstream of rotor.   Figure 9 shows the entropy downstream of the rotor at the DR2 plane. The high loss region is located near the shroud, where the flow is highly influenced by the rotor tip clearance. The rotor wake can be identified by the higher entropy region, as indicated by the dashed oval in the figure. The upper passage vortex region, combined with the tip leakage vortex accommodates more loss (indicated by 1). The lower passage vortex (indicated by 2) is much smaller than the upper passage vortex. However, it still widens the wake area near the hub. The rotor suction side shock and reflected pressure side shock can be identified in the Mach number distribution.   Figure 9 shows the entropy downstream of the rotor at the DR2 plane. The high loss region is located near the shroud, where the flow is highly influenced by the rotor tip clearance. The rotor wake can be identified by the higher entropy region, as indicated by the dashed oval in the figure.

Rotor Flow Field
The upper passage vortex region, combined with the tip leakage vortex accommodates more loss (indicated by 1). The lower passage vortex (indicated by 2) is much smaller than the upper passage vortex. However, it still widens the wake area near the hub. The rotor suction side shock and reflected pressure side shock can be identified in the Mach number distribution. The rotor flow field is strongly influenced by the stator wakes. In the time-resolved velocity flow field, strong variations of the velocity magnitude can be found at the rotor exit. The evolution of these variations is strongly influenced by the highly turbulent stator exit flow. Figure 10 presents the time averaged velocity, yaw angle and entropy at two locations downstream of the rotor (DR1 & DR2) over three blade pitches. The suction side shock and reflected pressure side shock can be seen in the velocity distribution in Figure 10a. No obvious difference can be seen in the results predicted by the unsteady and NLH methods downstream of rotor.  The rotor flow field is strongly influenced by the stator wakes. In the time-resolved velocity flow field, strong variations of the velocity magnitude can be found at the rotor exit. The evolution of these variations is strongly influenced by the highly turbulent stator exit flow. Figure 10 presents the time averaged velocity, yaw angle and entropy at two locations downstream of the rotor (DR1 & DR2) over three blade pitches. The suction side shock and reflected pressure side shock can be seen in the velocity distribution in Figure 10a. No obvious difference can be seen in the results predicted by the unsteady and NLH methods downstream of rotor.   Figure 9 shows the entropy downstream of the rotor at the DR2 plane. The high loss region is located near the shroud, where the flow is highly influenced by the rotor tip clearance. The rotor wake can be identified by the higher entropy region, as indicated by the dashed oval in the figure.

Rotor Flow Field
The upper passage vortex region, combined with the tip leakage vortex accommodates more loss (indicated by 1). The lower passage vortex (indicated by 2) is much smaller than the upper passage vortex. However, it still widens the wake area near the hub. The rotor suction side shock and reflected pressure side shock can be identified in the Mach number distribution. The rotor flow field is strongly influenced by the stator wakes. In the time-resolved velocity flow field, strong variations of the velocity magnitude can be found at the rotor exit. The evolution of these variations is strongly influenced by the highly turbulent stator exit flow. Figure 10 presents the time averaged velocity, yaw angle and entropy at two locations downstream of the rotor (DR1 & DR2) over three blade pitches. The suction side shock and reflected pressure side shock can be seen in the velocity distribution in Figure 10a. No obvious difference can be seen in the results predicted by the unsteady and NLH methods downstream of rotor.    Figure 11 shows the mean flow in the rotor passage. The high entropy region in the passage (indicated by 1) appears to be due to the travelling of chopped stator wakes through the passage. The rotor pressure side trailing edge shock, impinging on the suction side of the neighboring blade can be seen. The impinging shock interacts strongly with the rotor suction side boundary layer flow. The shock is then reflected and mixes with the rotor wake flow (dotted arrow) and decelerates the flow to subsonic conditions. There is also radial widening of the flow channel in the rotor section which also contributes to slow down the flow velocity at the exit of the rotor. Another weak shock on the trailing edge of rotor suction side is also evident in the Mach number distribution (solid arrow). This weak shock ultimately mixes with the rotor wake as well and causes a further flow deceleration in the rotor wake flow. The maximum Mach number appeared on 2/3 of the rotor chord on the suction side at 50% span. The presence of shocks generates entropy gradients in the flow passage. The shock-wake interactions at the rotor trailing edge causes significance flow unsteadiness downstream of the rotor. The higher unsteadiness will lead to the increased losses due to the wake and vortex shedding effects at the rotor trailing edge.
Energies 2018, 11, x 10 of 19 Figure 11 shows the mean flow in the rotor passage. The high entropy region in the passage (indicated by 1) appears to be due to the travelling of chopped stator wakes through the passage. The rotor pressure side trailing edge shock, impinging on the suction side of the neighboring blade can be seen. The impinging shock interacts strongly with the rotor suction side boundary layer flow. The shock is then reflected and mixes with the rotor wake flow (dotted arrow) and decelerates the flow to subsonic conditions. There is also radial widening of the flow channel in the rotor section which also contributes to slow down the flow velocity at the exit of the rotor. Another weak shock on the trailing edge of rotor suction side is also evident in the Mach number distribution (solid arrow). This weak shock ultimately mixes with the rotor wake as well and causes a further flow deceleration in the rotor wake flow. The maximum Mach number appeared on 2/3 of the rotor chord on the suction side at 50% span. The presence of shocks generates entropy gradients in the flow passage. The shockwake interactions at the rotor trailing edge causes significance flow unsteadiness downstream of the rotor. The higher unsteadiness will lead to the increased losses due to the wake and vortex shedding effects at the rotor trailing edge.  Figure 12 shows the time-space plots of entropy at the midspan, for the two locations DS2 & DR2, downstream of stator and rotor, respectively. The relative position of the stator vanes and the rotor blades is also shown in the figure. The moving blade effects appear under 45° structure format downstream of the rotor and the rotor wakes are indicated by the dashed line, whereas for the stator they appear horizontally. Within each blade passing period, the regions of high (indicated by dotted circle A) and low (indicated by dotted circle B) entropy appear. The high and low entropy levels are the evidence of the chopped stator wakes travelling through the rotor passages independently. The discrete spots of high and low entropy in the rotor wake (indicated by the dashed line) are the evidence of the wake shedding phenomena. The wake shedding is seemed to be phase locked to the blade passing period which will be further explained by the unsteady pressure distribution later in this section. Seven discrete structures of high entropy appear in each blade passing period (indicated by the dotted oval). If the wake shedding is not phase-locked to the passing blade then, there should be a continuous stripe of high entropy in the wake region. The wake shedding is also observed at the stator trailing edge where high entropy regions appear in discrete mode.   Figure 12 shows the time-space plots of entropy at the midspan, for the two locations DS2 & DR2, downstream of stator and rotor, respectively. The relative position of the stator vanes and the rotor blades is also shown in the figure. The moving blade effects appear under 45 • structure format downstream of the rotor and the rotor wakes are indicated by the dashed line, whereas for the stator they appear horizontally. Within each blade passing period, the regions of high (indicated by dotted circle A) and low (indicated by dotted circle B) entropy appear. The high and low entropy levels are the evidence of the chopped stator wakes travelling through the rotor passages independently. The discrete spots of high and low entropy in the rotor wake (indicated by the dashed line) are the evidence of the wake shedding phenomena. The wake shedding is seemed to be phase locked to the blade passing period which will be further explained by the unsteady pressure distribution later in this section. Seven discrete structures of high entropy appear in each blade passing period (indicated by the dotted oval). If the wake shedding is not phase-locked to the passing blade then, there should be a continuous stripe of high entropy in the wake region. The wake shedding is also observed at the stator trailing edge where high entropy regions appear in discrete mode.

Stator-Rotor Interaction
Energies 2018, 11, x 10 of 19 Figure 11 shows the mean flow in the rotor passage. The high entropy region in the passage (indicated by 1) appears to be due to the travelling of chopped stator wakes through the passage. The rotor pressure side trailing edge shock, impinging on the suction side of the neighboring blade can be seen. The impinging shock interacts strongly with the rotor suction side boundary layer flow. The shock is then reflected and mixes with the rotor wake flow (dotted arrow) and decelerates the flow to subsonic conditions. There is also radial widening of the flow channel in the rotor section which also contributes to slow down the flow velocity at the exit of the rotor. Another weak shock on the trailing edge of rotor suction side is also evident in the Mach number distribution (solid arrow). This weak shock ultimately mixes with the rotor wake as well and causes a further flow deceleration in the rotor wake flow. The maximum Mach number appeared on 2/3 of the rotor chord on the suction side at 50% span. The presence of shocks generates entropy gradients in the flow passage. The shockwake interactions at the rotor trailing edge causes significance flow unsteadiness downstream of the rotor. The higher unsteadiness will lead to the increased losses due to the wake and vortex shedding effects at the rotor trailing edge.  Figure 12 shows the time-space plots of entropy at the midspan, for the two locations DS2 & DR2, downstream of stator and rotor, respectively. The relative position of the stator vanes and the rotor blades is also shown in the figure. The moving blade effects appear under 45° structure format downstream of the rotor and the rotor wakes are indicated by the dashed line, whereas for the stator they appear horizontally. Within each blade passing period, the regions of high (indicated by dotted circle A) and low (indicated by dotted circle B) entropy appear. The high and low entropy levels are the evidence of the chopped stator wakes travelling through the rotor passages independently. The discrete spots of high and low entropy in the rotor wake (indicated by the dashed line) are the evidence of the wake shedding phenomena. The wake shedding is seemed to be phase locked to the blade passing period which will be further explained by the unsteady pressure distribution later in this section. Seven discrete structures of high entropy appear in each blade passing period (indicated by the dotted oval). If the wake shedding is not phase-locked to the passing blade then, there should be a continuous stripe of high entropy in the wake region. The wake shedding is also observed at the stator trailing edge where high entropy regions appear in discrete mode.   Figure 13 shows the each individual harmonic pressure amplitude on the stator and rotor surface at various spanwise positions (LE: leading edge; SS: suction surface; TE: trailing edge; PS: pressure surface; LE: leading edge). The pressure amplitude is dominated by the first and second order of harmonics, because the low frequency components are the dominant components of the periodic flow field. This also demonstrates that three orders of harmonic are well capable of predicting the main flow parameters. The pressure amplitudes on the suction surface are higher than on the pressure surface, since the rotor rotates relative to the stator from stator suction surface to pressure surface. This is different from that in compressors [33].

Stator-Rotor Interaction
Energies 2018, 11, x 11 of 19 Figure 13 shows the each individual harmonic pressure amplitude on the stator and rotor surface at various spanwise positions (LE: leading edge; SS: suction surface; TE: trailing edge; PS: pressure surface; LE: leading edge). The pressure amplitude is dominated by the first and second order of harmonics, because the low frequency components are the dominant components of the periodic flow field. This also demonstrates that three orders of harmonic are well capable of predicting the main flow parameters. The pressure amplitudes on the suction surface are higher than on the pressure surface, since the rotor rotates relative to the stator from stator suction surface to pressure surface. This is different from that in compressors [33]. The transport of the stator wakes through the rotor passages can be well explained by the entropy contour in Figure 14. The stator wakes are chopped by the rotor leading edge and travel independently through the rotor passage and influence the rotor trailing edge flow. The flow field exhibits strong periodic unsteady flow phenomena. There is a phase locking between the wake shedding and the blade passing period which is evident from the instantaneous entropy distribution. The chopped stator wakes accumulate towards the suction side of the rotor and interact strongly with the rotor trailing edge wake. With three orders of harmonic perturbation, the NLH method is fully capable of capturing the physical phenomena of the stator wake transport through the rotor passage. The transport of the stator wakes through the rotor passages can be well explained by the entropy contour in Figure 14. The stator wakes are chopped by the rotor leading edge and travel independently through the rotor passage and influence the rotor trailing edge flow. The flow field exhibits strong periodic unsteady flow phenomena. There is a phase locking between the wake shedding and the blade passing period which is evident from the instantaneous entropy distribution. The chopped stator wakes accumulate towards the suction side of the rotor and interact strongly with the rotor trailing edge wake. With three orders of harmonic perturbation, the NLH method is fully capable of capturing the physical phenomena of the stator wake transport through the rotor passage. The wake passing the rotor-stator interface is slightly wider with the NLH method as compared with the unsteady time-accurate method, but it still shows acceptable accuracy. In the NLH solution, a cut-off of flow variables occurs when the information passes across the rotor-stator interface. The periodic unsteady perturbations are Fourier decomposed at the interface therefore; the stator wakes and shocks are smeared. The continuity of the wakes and shocks across the interface would improve with higher numbers of harmonics. The smearing of the data at the rotor-stator interface is usual, as the numerical mixing takes place on these interfaces. However, for the unsteady time-domain computations, there is no such interface and transportation of flow quantities through the interface is smooth.
The reason for the stator wake phase locking can be explained by the unsteady pressure distribution in Figure 15. The stator trailing edge shock wave hits the passing by rotor blade at t/τ = 2/8, creating a strong circular pressure wave, which is periodically reflected back to the stator suction side (SS) at t/τ = 4/8, thus triggering the phenomena of wake shedding. A part of this pressure wave travels downstream through the rotor passage (t/τ = 6/8) and triggers the wake shedding at the rotor trailing edge as well. The pressure distribution downstream of the rotor-stator interface is smoother with the unsteady method than the NLH method. However, the NLH method has captured the pressure wave trend in good agreement with the unsteady method using only limited numbers of harmonics. The inaccuracies in the rotor-stator interface vicinity could be further reduced with the increased numbers of harmonics. The wake passing the rotor-stator interface is slightly wider with the NLH method as compared with the unsteady time-accurate method, but it still shows acceptable accuracy. In the NLH solution, a cutoff of flow variables occurs when the information passes across the rotor-stator interface. The periodic unsteady perturbations are Fourier decomposed at the interface therefore; the stator wakes and shocks are smeared. The continuity of the wakes and shocks across the interface would improve with higher numbers of harmonics. The smearing of the data at the rotor-stator interface is usual, as the numerical mixing takes place on these interfaces. However, for the unsteady time-domain computations, there is no such interface and transportation of flow quantities through the interface is smooth. The reason for the stator wake phase locking can be explained by the unsteady pressure distribution in Figure 15. The stator trailing edge shock wave hits the passing by rotor blade at t/τ = 2/8, creating a strong circular pressure wave, which is periodically reflected back to the stator suction side (SS) at t/τ = 4/8, thus triggering the phenomena of wake shedding. A part of this pressure wave travels downstream through the rotor passage (t/τ = 6/8) and triggers the wake shedding at the rotor trailing edge as well. The pressure distribution downstream of the rotor-stator interface is smoother with the unsteady method than the NLH method. However, the NLH method has captured the pressure wave trend in good agreement with the unsteady method using only limited numbers of harmonics. The inaccuracies in the rotor-stator interface vicinity could be further reduced with the increased numbers of harmonics.   Figure 16 represents the stator and rotor surface pressure distribution at 50% span. Timeaveraged data has been presented in comparison with the reconstructed instantaneous data at various time steps to monitor the flow unsteadiness. The time averaged and instantaneous pressure distribution on the stator pressure side is the same except at the trailing edge. However, on the suction side, fluctuations near the trailing edge can be seen. This is where the strong unsteady statorrotor interaction occurs. Furthermore, the harmonic pressure amplitude in Figure 13 also shows that the pressure fluctuations on the stator suction side are higher. This is mainly because the rotor rotates relative to the stator from stator suction surface to pressure surface. The boundary layer near the trailing edge of the suction side is thicker may also cause the potential flow interaction to spread wider in that region.
The fluctuations of the surface pressure on the rotor are large compared with the stator. The whole chord of the rotor is influenced by the unsteadiness caused by the stator wake. The pressure fluctuations on the rotor suction side are much higher compared with the pressure side, which is in line with the findings of the harmonic pressure amplitude in Figure 13. The presence of the shock wave and the accumulation of the chopped stator wake on rotor suction side cause more periodic fluctuations and unsteadiness compared with the pressure side. Overall, the stator-rotor interaction can affect about one third of the stator axial chord near trailing edge; while the whole chord of the rotor is influenced. This indicates that the stator-rotor interactions in this turbine stage are very significant.  Figure 16 represents the stator and rotor surface pressure distribution at 50% span. Time-averaged data has been presented in comparison with the reconstructed instantaneous data at various time steps to monitor the flow unsteadiness. The time averaged and instantaneous pressure distribution on the stator pressure side is the same except at the trailing edge. However, on the suction side, fluctuations near the trailing edge can be seen. This is where the strong unsteady stator-rotor interaction occurs. Furthermore, the harmonic pressure amplitude in Figure 13 also shows that the pressure fluctuations on the stator suction side are higher. This is mainly because the rotor rotates relative to the stator from stator suction surface to pressure surface. The boundary layer near the trailing edge of the suction side is thicker may also cause the potential flow interaction to spread wider in that region.
The fluctuations of the surface pressure on the rotor are large compared with the stator. The whole chord of the rotor is influenced by the unsteadiness caused by the stator wake. The pressure fluctuations on the rotor suction side are much higher compared with the pressure side, which is in line with the findings of the harmonic pressure amplitude in Figure 13. The presence of the shock wave and the accumulation of the chopped stator wake on rotor suction side cause more periodic fluctuations and unsteadiness compared with the pressure side. Overall, the stator-rotor interaction can affect about one third of the stator axial chord near trailing edge; while the whole chord of the rotor is influenced. This indicates that the stator-rotor interactions in this turbine stage are very significant. proposed the first model using overlapping grids for all blade rows without interface planes. In the past 30 years, many efforts have been made and several deterministic correlations models have been proposed and showed that the models apparently improve the numerical results in comparison with steady mixing plane method. However, none of these models has been widely used in the turbomachinery routine design because these models are either too complex, too temperamental or still too time consuming for routine use. Hence, more research work for modeling the deterministic correlations should be conducted.
In the NLH method the deterministic stresses appear through the time averaging of the periodic unsteadiness due to the relative rotor-stator movement. The level of the deterministic stresses in the flow field is representative of the unsteady effects like wakes and potential effects. The deterministic stress as well as the deterministic enthalpy both represent the unsteadiness. The deterministic stresses are directly computed from harmonic values. The detailed information about the flow unsteadiness appears in nine deterministic terms, of which six are related to velocity-velocity deterministic stresses (due to symmetry) and three refer to enthalpy-velocity deterministic fluxes as shown in the Equation (7): Figure 17 shows the three orders of harmonic deterministic stress DCS_ZZ and deterministic enthalpy DCH_Z at 50% span (only one component is presented for comparison). The levels at the stator inlet are very low, however at the exit of the stator, the levels appeared to be higher due to the unsteady potential field from the downstream rotating blades. The deterministic stress and enthalpy level are also higher in the rotor passage due to the unsteadiness from the upstream stator wake. Compared to the pressure side, the level of unsteadiness on the suction side of the rotor blade is higher. The reason is that the chopped stator wakes tends to accumulate towards the suction side of the rotor blades.  proposed the first model using overlapping grids for all blade rows without interface planes. In the past 30 years, many efforts have been made and several deterministic correlations models have been proposed and showed that the models apparently improve the numerical results in comparison with steady mixing plane method. However, none of these models has been widely used in the turbomachinery routine design because these models are either too complex, too temperamental or still too time consuming for routine use. Hence, more research work for modeling the deterministic correlations should be conducted.
In the NLH method the deterministic stresses appear through the time averaging of the periodic unsteadiness due to the relative rotor-stator movement. The level of the deterministic stresses in the flow field is representative of the unsteady effects like wakes and potential effects. The deterministic stress as well as the deterministic enthalpy both represent the unsteadiness. The deterministic stresses are directly computed from harmonic values. The detailed information about the flow unsteadiness appears in nine deterministic terms, of which six are related to velocity-velocity deterministic stresses (due to symmetry) and three refer to enthalpy-velocity deterministic fluxes as shown in the Equation (7): Figure 17 shows the three orders of harmonic deterministic stress DCS_ZZ and deterministic enthalpy DCH_Z at 50% span (only one component is presented for comparison). The levels at the stator inlet are very low, however at the exit of the stator, the levels appeared to be higher due to the unsteady potential field from the downstream rotating blades. The deterministic stress and enthalpy level are also higher in the rotor passage due to the unsteadiness from the upstream stator wake. Compared to the pressure side, the level of unsteadiness on the suction side of the rotor blade is higher. The reason is that the chopped stator wakes tends to accumulate towards the suction side of the rotor blades.   Figure 19 presents a one-dimensional variation of the deterministic stresses and deterministic enthalpy fluxes along the mean flow direction. The variations are rather complex and differ from one another. The amplitude of the deterministic terms related to the X velocity (DCS_XX, DCS_XY, DCS_XZ and DCH_X) appears much smaller than the deterministic terms in the Y and Z direction, for both stator and rotor passages. Similar behavior is observed in the region between stator and rotor near the R/S interface. This leads to the conclusion that the flow unsteadiness is stronger in Y and Z direction as compared with the X direction. Therefore, the deterministic components in the Y and Z direction are particularly related to the unsteady stator-rotor interactions. The maximum values appear in the region of the stator TE, the R/S interface and the rotor LE, which is the region mainly influenced by the unsteady stator-rotor interactions. There is a gap (discontinuity in deterministic terms) at the R/S interface because they are calculated at the interface in the stator and rotor passage respectively. Overall a decreasing trend in the deterministic terms appears from the rotor LE to the outlet, however, the large oscillation of the terms may be attributed to the chopped stator wakes, which travel independently through the rotor passage, creating fluctuations for the deterministic terms. Again, this demonstrates that the deterministic stresses and deterministic enthalpy terms in the NLH method are capable to predict the unsteady flow.
From the above analysis it is inferred that, the distributions of deterministic stress and deterministic enthalpy are quite different from that in compressors [25]. Usually in compressors, the deterministic stress and enthalpy have large magnitudes at the rotor-stator interface and decrease rapidly in the flow direction downstream of the interface; however, there is no obvious similar behavior in the turbine. This is because in compressors, there exists strong inviscid stretching due to the divergent passage, and the viscous decay for wakes. Both of these two mechanisms make the deterministic stresses and deterministic enthalpy decay rapidly. This indicates that, some empirical   Figure 18 demonstrates the results at 95% span. The results show higher levels of unsteadiness on the suction side of the rotor where the tip leakage vortex is the dominant source of unsteadiness. The higher levels appear roughly from the leading edge of the rotor and they are in-line with the tip clearance flow. The introduction of the deterministic terms in the NLH method produces the unsteady results with precision and good level of confidence. The other deterministic stresses and enthalpy components also show the same behavior. Only some of them are presented here due to the paper length limitations.  Figure 19 presents a one-dimensional variation of the deterministic stresses and deterministic enthalpy fluxes along the mean flow direction. The variations are rather complex and differ from one another. The amplitude of the deterministic terms related to the X velocity (DCS_XX, DCS_XY, DCS_XZ and DCH_X) appears much smaller than the deterministic terms in the Y and Z direction, for both stator and rotor passages. Similar behavior is observed in the region between stator and rotor near the R/S interface. This leads to the conclusion that the flow unsteadiness is stronger in Y and Z direction as compared with the X direction. Therefore, the deterministic components in the Y and Z direction are particularly related to the unsteady stator-rotor interactions. The maximum values appear in the region of the stator TE, the R/S interface and the rotor LE, which is the region mainly influenced by the unsteady stator-rotor interactions. There is a gap (discontinuity in deterministic terms) at the R/S interface because they are calculated at the interface in the stator and rotor passage respectively. Overall a decreasing trend in the deterministic terms appears from the rotor LE to the outlet, however, the large oscillation of the terms may be attributed to the chopped stator wakes, which travel independently through the rotor passage, creating fluctuations for the deterministic terms. Again, this demonstrates that the deterministic stresses and deterministic enthalpy terms in the NLH method are capable to predict the unsteady flow.
From the above analysis it is inferred that, the distributions of deterministic stress and deterministic enthalpy are quite different from that in compressors [25]. Usually in compressors, the deterministic stress and enthalpy have large magnitudes at the rotor-stator interface and decrease rapidly in the flow direction downstream of the interface; however, there is no obvious similar behavior in the turbine. This is because in compressors, there exists strong inviscid stretching due to the divergent passage, and the viscous decay for wakes. Both of these two mechanisms make the deterministic stresses and deterministic enthalpy decay rapidly. This indicates that, some empirical  Figure 19 presents a one-dimensional variation of the deterministic stresses and deterministic enthalpy fluxes along the mean flow direction. The variations are rather complex and differ from one another. The amplitude of the deterministic terms related to the X velocity (DCS_XX, DCS_XY, DCS_XZ and DCH_X) appears much smaller than the deterministic terms in the Y and Z direction, for both stator and rotor passages. Similar behavior is observed in the region between stator and rotor near the R/S interface. This leads to the conclusion that the flow unsteadiness is stronger in Y and Z direction as compared with the X direction. Therefore, the deterministic components in the Y and Z direction are particularly related to the unsteady stator-rotor interactions. The maximum values appear in the region of the stator TE, the R/S interface and the rotor LE, which is the region mainly influenced by the unsteady stator-rotor interactions. There is a gap (discontinuity in deterministic terms) at the R/S interface because they are calculated at the interface in the stator and rotor passage respectively. Overall a decreasing trend in the deterministic terms appears from the rotor LE to the outlet, however, the large oscillation of the terms may be attributed to the chopped stator wakes, which travel independently through the rotor passage, creating fluctuations for the deterministic terms. Again, this demonstrates that the deterministic stresses and deterministic enthalpy terms in the NLH method are capable to predict the unsteady flow.
From the above analysis it is inferred that, the distributions of deterministic stress and deterministic enthalpy are quite different from that in compressors [25]. Usually in compressors, the deterministic stress and enthalpy have large magnitudes at the rotor-stator interface and decrease rapidly in the flow direction downstream of the interface; however, there is no obvious similar behavior in the turbine. This is because in compressors, there exists strong inviscid stretching due to the divergent passage, and the viscous decay for wakes. Both of these two mechanisms make the deterministic stresses and deterministic enthalpy decay rapidly. This indicates that, some empirical deterministic correlations models based on the decay concept of compressors are not suitable for turbines.

Conclusions
The NLH method has proved to be very efficient compared to the time-domain unsteady method. The frequency casting allows the solving of the harmonic transport equations in a similar way to that of the steady flow equations with a single blade passage, thus the time required to solve the additional harmonic variables is reduced compared with the unsteady method. In the current study, the computational time required by the NLH method is seven times smaller than that of the unsteady RANS. The NLH method has proved to be a promising way to study the stator-rotor interaction and can be used as a fast method for unsteady flow predictions at the design stage.
(1). The stator wakes are chopped as they enter the rotor passage by the rotor leading edge. The chopped wake segments travel through the rotor passage independently. They also interfere with the rotor trailing edge shocks, creating a source of unsteadiness for the trailing edge flow. The R/S interface treatment in the NLH method has shown good capability for reproducing the unsteady signals across the interface. Some minor differences can be seen across the R/S interface, which can be further satisfied by increasing the numbers of harmonics. (2). Static pressure distribution has been used to explain the triggering mechanism and the phaselocking of the wake shedding phenomena. The pressure wave generated by the reflection of the stator trailing edge shock on the suction side of the passing rotor blade is responsible for the phase locking. These pressure waves travel upstream and enforce the wake shedding as they strike the rear section of the stator suction side, where it interacts with the boundary layer. These pressure waves are strong enough to travel downstream within the rotor passage and influence the rotor trailing edge flow. The impingement of the stator shock and the accumulation of chopped stator wakes on the rotor suction side have induced higher periodic fluctuations on the suction side compared with the pressure side.

Conclusions
The NLH method has proved to be very efficient compared to the time-domain unsteady method. The frequency casting allows the solving of the harmonic transport equations in a similar way to that of the steady flow equations with a single blade passage, thus the time required to solve the additional harmonic variables is reduced compared with the unsteady method. In the current study, the computational time required by the NLH method is seven times smaller than that of the unsteady RANS. The NLH method has proved to be a promising way to study the stator-rotor interaction and can be used as a fast method for unsteady flow predictions at the design stage.
(1) The stator wakes are chopped as they enter the rotor passage by the rotor leading edge.
The chopped wake segments travel through the rotor passage independently. They also interfere with the rotor trailing edge shocks, creating a source of unsteadiness for the trailing edge flow. The R/S interface treatment in the NLH method has shown good capability for reproducing the unsteady signals across the interface. Some minor differences can be seen across the R/S interface, which can be further satisfied by increasing the numbers of harmonics. (2) Static pressure distribution has been used to explain the triggering mechanism and the phase-locking of the wake shedding phenomena. The pressure wave generated by the reflection of the stator trailing edge shock on the suction side of the passing rotor blade is responsible for the phase locking. These pressure waves travel upstream and enforce the wake shedding as they strike the rear section of the stator suction side, where it interacts with the boundary layer. These pressure waves are strong enough to travel downstream within the rotor passage and influence the rotor trailing edge flow. The impingement of the stator shock and the accumulation of chopped stator wakes on the rotor suction side have induced higher periodic fluctuations on the suction side compared with the pressure side. (3) The distribution of the deterministic correlations shows that the flow unsteadiness is particularly stronger in the Y and Z directions. Furthermore, the distributions of the deterministic correlations indicate that, some empirical deterministic correlations models based on the decay concept of compressors are not suitable for turbines.