One-Stage Power Turbine Preliminary Design and Analysis

ABSTRACT: Design and development of gas turbine components are a complex multidisciplinary process. At the beginning of the power class definition and engine configuration it is necessary to conduct a market study. The results obtained are used in gas turbine thermodynamic cycle calculations and analysis in order to define the gas turbine design point. Several possible design points are evaluated during this procedure. After this step, the gas turbine components are designed, including: compressor, combustion chamber and turbine. For industrial gas turbine purposes, it is common to use a free turbine after the gas generator, also commonly named power turbine. In this work, a power turbine was initially designed by meanline techniques, considering internal loss mechanisms, to obtain the main dimensions. The geometries of the components were generated in a 3-D environment to make possible the mesh generation, process to discretize the physical domain into a computational domain and use a 3-D Computational Fluid Dynamics tool. The results from the meanline approach and from the 3-D turbulent flow numerical simulations were compared to verify the turbine operational conditions and its predictions at design and offdesign conditions. The gas turbine under study is a project, derived from a low thrust turbojet previously developed by Instituto de Aeronáutica e Espaço. The power turbine project uses the same turbojet gas generator, already designed and currently under tests.


INTRODUCTION
The Instituto de Aeronáutica e Espaço (IAE) is one of the Brazilian government institutes responsible for conducting research and development related to aerospace engineering. IAE has been developing in conjunction with TGM industry and Instituto Tecnológico de Aeronáutica (ITA) a gas turbine project able to be used in two types of applications: a turbojet -in thrust class of 5 kN, and a turboshaft -in a shaft power class of 1.2 MW for power generation application.
The aim is to develop a gas turbine in which the economic feasibility can be justified by its different usages, since its two classes employ the same gas generator. The first version of the designed turbojet is installed in a test facility at IAE and its test protocol is under progress. The turboshaft class was studied after the design of the turbojet engine version. The power turbine preliminary design (Silva, 2012) was finished, but improvements in the aerothermodynamics, heat transfer, stress and strain calculations of the components are being performed.
This work presents the design of the power turbine, its aerothermodynamics, which was performed using a reducedorder numerical design tool and Computational Fluid Dynamics (CFD) simulations to verify the results from the preliminary design made by Silva (2012), providing the data needed for structural and heat transfer analysis, to be done afterwards. An in-house code written in FORTRAN language (Silva, 2012) was used to supply input data for commercial turbomachinery design software AxSTREAM® and to determine a preliminary design of a free power turbine. With the use of AxSTREAM® software the turbine Nozzle Guide Vane (NGV) and rotor blades were generated by an airfoil stacking procedure and the 3-D geometry was obtained. This 3-D geometry was carefully treated to start the mesh generation process using the software ANSYS TurboGrid®. The CFD solver used in this work was the ANSYS CFX® v.13.

POWER TURBINE DESIGN
The preliminary design process was performed using two different meanline numerical tools. The first is an in-house tool developed by (Silva, 2012) and the second one is commercial software called AxStream®. Meanline means that all fluid properties including the turbomachinery dimensions are calculated at mean blade height. This technique neglects the viscous effects that are accounted using loss modelling. There are several loss models, each one adequate for each case; therefore, there is no general loss model. The average values of all fluid properties from hub-to-tip are very close to the values obtained at mean radius. Thereby, the use of meanline technique is very interesting to obtain fast results and solutions for a first turbomachinery design process.
The main results from the meanline tool developed by (Silva, 2012) are presented in Table 1, in which the loss model developed by Kacker and Okapuu (1982) was used. The power turbine preliminary design was performed using AxSTREAM® software and it was conducted with the Craig and Cox (1970) loss model. Both loss models can be used at design and off-design conditions. For the power turbine studies, both loss models are adequate. For turbomachinery design several parameters must be varied and imposed by the designer during the preliminary sizing. Some design parameters and values were determined using an in-house code developed by Silva (2012) and the results used as an inlet data at the commercial software.
Along the years, several loss models were developed, for example: Ainley and Mathieson (1951), Denton (1978), Dunham and Smith (1968), Dunham and Came (1970), Dunham and Panton (1973), Rody et al. (1980) and Smith (1965). There are other models, but only the most commonly used were mentioned. In the present work, the loss models used are not described due to the high number of information, mainly about its calibrations, validation ranges, loss sources, constraints and so on. The reader that would like to understand the models in details it is highly recommended to search for the publication sources above mentioned.
The two-dimensional view of the power turbine stage, in the streamwise direction, calculated by AxSTREAM® making the use of the inlet data, calculated by the in-house code developed by Silva (2012), as presented in Table 1, is shown in Fig. 1a. The blue component is the stator row (or NGV) and the red component is the rotor row.

THE 3-D TURBINE FLOW CALCULATIONS
Since the aim of this study is to discuss the preliminary design results obtained previously Silva (2012), in this paper there are no discussions about numerical schemes, turbulence models, spatial and time marching methods, discretization issues and other topics around the CFD numerical implementations and its programming techniques. Certainly, all software and hardware used in numerical simulations demand high CPU time, currently called high performance computing (or HPC).
In this paper, the methods and models used are available in the open turbomachinery literature and additional details can be found in Craig and Cox (1970), Kacker and Okapuu (1982), ), Bringhenti et al. (2001, Tomita (2009),  and Murari et al. (2011).
The turbine geometry calculated based on meanline technique was appropriately treated to generate its 3-D geometrical view making use of Computer-Aided Design (CAD) software. It is recognized that the CAD treatments, mesh generation processes and the flow calculations using CFD technique are very time consuming processes.
The use of CFD techniques is very important during the turbomachinery design processes, because some flow features cannot be visualized by the designer from meanline calculations, for example: boundary layer separation, high pressure gradients regions, shock wave position and expansion waves. Thus, in order to calculate the 3-D flow characteristics, the commercial CFD software ANSYS CFX® was used, with the calculations based on the Reynolds-Averaged Navier-Stokes equations (RANS). The turbulence effects were accounted using eddy viscosity turbulence model.

MESH GENERATION PROCESS
The mesh was generated using ANSYS TurboGrid® software. As this software is part of the environment, it also allows a direct link with other ANSYS® package software, such as CFX® for CFD calculation.
Different mesh configurations can be found in the TurboGrid® software. It is important to consider the control volume quality as its orthogonality, skew angles, aspect ratios, smoothing and others as presented in the references ,  and Menter (1983).
Four meshes with different refinements were generated to evaluate the mesh dependency for the particular turbine geometry studied. These meshes compromise hexahedral elements with the H/J/C/L configurations to better accommodate the turbine geometry. At this point, it is important to mention that the y + (position from blade surface at perpendicular direction) requirements must be considered, depending on the turbulence model that will be used to calculate the flow eddy viscosity. The Reynolds number was defined relative to the chord, being 4.03 x 10 5 at stator row and 2.68 x 10 5 at rotor row.
The y + value consists in determination of the distance between the first centroid (or vertex) of the control volume at wall from the wall surface, ensuring that this distance is consistent with the turbulence model. The y + is defined by the equation: where: ρ and μ: density and dynamic viscosity, respectively; Δy: y variation; τ w : shear stress at wall.
Generally, for y + ≤ 2, a very refined mesh must be generated. For y + > 2, the CFD software determines the velocity field using wall functions. It is recommended by CFD solver that the best practice to use the wall functions is for y + values around 30 ≤ y + ≤ 300. Different ranges can be found in the literature, but it is highly dependent on the wall functions implemented in each computational program. Table 2 shows four mesh configurations generated to perform the mesh dependency study. The meshes were generated and the fluid mechanics equations were calculated until the point that there are no fluid properties variation in the solution from CFD calculations. These fluid and flow characteristics variations were considered in specific streamwise locations (section 1, 2 and 3 as shown in Fig. 1b) and from blade hub-to-tip positions. For all cases, the numerical settings as spatial and time integration used in simulations are the same and the two-equation turbulence   . . model developed by Menter et al. (2004), called Shear-Stress Transport (SST), was used to calculate the flow eddy viscosity. The details are described in the next section. The Mach number, total pressure and static pressure were compared for all mesh configurations, and the values were calculated using mass averaging at each surface (section). Moreover, the isentropic efficiency, mass flow, expansion ratio and y + values were monitored and analyzed for each mesh as shown in Table 3. The reference pressure value is 101,325 Pa.   In Fig. 2, at section 1, it is possible to observe the impact of the velocity variations in different mesh configurations and its influence on Mach number and consequently on total pressure. The meshes R and R+ presented closely results. Note that, for 0 -20% of blade span and for 80 -100% of blade span, the results from different meshes present some interesting differences. This is due to the inadequate mesh resolution close to the wall, where boundary-layer is very important and pronounced.
From Fig. 3, at section 2, the same behavior and explanation can be performed. The meshes R and R+ presented closely results including the mesh F. At this point, it is possible to mention that the mesh G is not recommended for the simulations like the ones used in this work.
In Fig. 4, at section 3, the mesh F presented small differences in the solution when compared with the solutions from meshes R and R+. Again, the meshes R and R+ presented closely results including regions close to the turbine endwall and casing.

Relative blade span
From these results it is possible to define that the mesh R is more adequate for the numerical simulation of the turbine designed in this study. The mesh R+ is also an alternative, but it is more expensive (CPU time) than mesh R due to the high number of control volumes, becoming the mesh R the best compromise. Figure 5 shows the surface mesh characteristics at stator blade row. Regions close to the wall surfaces need more refinement, the smoothing at wall regions in the computational domain should be carefully built to ensure a good resolution of the boundary-layer and acceptable flow characteristics resolution. method for the mass terms, similar to the methodology proposed by Majumdar (1998).
For the discretization of convective terms from the Navier-Stokes equations, first and second order upwind numerical schemes are available. The limiters proposed by Barth and Jesperson (1989) are used to control the numerical discretization order, mainly for discontinuously regions, and to avoid numerical instabilities in shock regions when upwind based methods are used. The diffusion terms are discretized making the use of the shape functions. This methodology is usual in finite elements method (FEM) and has shown to be very robust for this purpose. The discretization model applied for the convective terms from Navier-Stokes equations was the second-order upwind scheme.
To accelerate the numerical procedure, the anisotropic algebraic multigrid technique developed by Raw (1996) was used to avoid error propagation due to the irregular elements. Conservation requirements are imposed by the method during the coarsening process. CFX® has only implicit time-marching scheme for time integration. An automatic time-step calculation method was used to ensure good numerical stability during the simulation. Fast convergence with few numbers of iterations can be achieved using this implicit scheme. All simulations were done for steady-state regime.
To account the turbulence effects, the flow eddy viscosity was calculated based on the two-equation Shear Stress Transport (SST) model Menter (1983). This model is adequate for the present analysis according to the recommendations and analysis from Menter et al. (2003Menter et al. ( , 2004, Denton (2010), and CFD Online (2011). The model uses the k-ω formulation at the internal region of the boundary layer and the k-ε model away from the wall. This control is determined using blending functions developed by Menter (1983).
The boundary-conditions were set at several surfaces as described:

NUMERICAL ISSUES AND BOUNDARY CONDITIONS
The commercial software CFX® v.13 was used as a CFD solver to perform the flow calculations based on the continuity, Navier-Stokes and energy equations. The compressible version of the solver was set-up. CFX® uses the coupling between the equations and the modified Rhie and Chow (1982)  Inter-rows: mixing-plane approach. The solution from CFD solver was monitored by convergence criteria based on the residue decayment, number of iterations and mass-flow variations (if no changes in the mass-flow was detected, the solution can be considered achieved; this is common in the turbomachinery CFD applications). Figure 6a shows the residual decayment of the continuity and momentum equations monitored during numerical iterations.
The mass-flow at turbine inlet was monitored during numerical iterations. Figure 6b shows the mass-flow behavior during the numerical simulation. With 150 iterations the solution was considered achieved. An important numerical issue to mention is the very effective initial condition, also called initialization process used in the CFX® solver, in which, before the first iteration, the solver estimates a good distribution of pressure, temperature and velocity fields along the turbine streamwise and blade span directions, to guarantee a good initial condition.

DESiGN-POiNT ANAlYSiS
At design-point condition, the turbine flowfield was calculated using CFD technique and each blade row (NGV and rotor) was analysed separately. The flow characteristics in important and/or critical regions were studied, as near to hub and tip zones, as well as the average values for some parameters at stage inlet and outlet sections. The static pressure distribution along each blade profile is plotted against its chord to evaluate the behavior of the flow along the blade pressure and suction sides. The blade profile is a multi-circular arc (MCA).
The NGV Mach number and entropy contours for 20, 50 and 80% spanwise are shown in Figs. 7 to 9. It is important to mention that the scales considered in each figure are not held constant; this was done for better visualization. From Fig 7, it is possible to observe high flow acceleration from NGV leading-to-trailing edge, but at 68% of the blade chord a sudden increase in static pressure is observed, causing flow deceleration. The result is an increase in the entropy. Figure 8 show the same analysis, but for 50% of the NGV blade span. Note that, at suction side, there are quite differences in static pressure values when compared with the 20% of the NGV spanwise. This difference is because the flow velocities are higher in regions close to the NGV hub. The degree of flow deceleration at 50% of the NGV spanwise is lower when compared with 20% of the NGV spanwise. Hence, the internal losses are lower and, consequently, the entropy generation is decreased too. Figure 9 show the results for 80% of NGV spanwise. Comparing the results of static pressure for different NGV blade span, it is possible to observe that the flow acceleration region, close to the blade leading edge, at blade suction side, changes from around 12 to 18% of the blade chord. The blade angles and thickness chord ratio distributions are responsible for these variations. Figure 9b shows small Mach number values when compared with Fig. 8b. This is a common flow behavior along the NGV spanwise. The entropy generation are indications that some improvements in the blade design can be necessary to improve the flow characteristics through the NGV channel.
For the blade pressure side, the results show a satisfactory pressure variation. At suction side, close to the hub, there is a supersonic bubble, covering the full flow passage. Except for this bubble influence, the results are in agreement with the prediction . Figure 10a provides a visualization of the velocity vectors related with the flow acceleration at suction side at 20% of the NGV spanwise.   The total pressure and Mach number distributions, in a midplane between blades of the same row, are shown in Figs. 10b and 10c. In these figures, it is possible to observe the radial variation of the flow velocity and the high velocity region at NGV hub close to the trailing edge. The gas expansion is observed from the total pressure decreasing along the NGV blade chord.

Number os iterations
The result also indicates that higher velocities at trailing edge, near to the hub, increase the internal losses due to the association with skin friction.
The total pressure distribution shows the influence of boundary-layer in the turbine hub and casing walls. The total pressure, Mach number and entropy distributions are shown at stator inlet and outlet sections (Figs. 11 to 13).
The influence of boundary-layer should be carefully analyzed during the turbomachinery design process. Depending on the boundary-layer thickness, the influence on the velocity field is expressive, because an increase in velocity changes all flow properties distributions and generates more internal losses, degrading the machine efficiency. The variation of the velocity along the streamwise direction is important due to the influence on the gas expansion process related with the density variation.
Entropy generation sources are generated at walls, due to the flow viscous effects and high velocities that are found in the blade-to-blade passage area; the velocity increases because the blade passage area is a convergent nozzle and high flow deflections are needed to improve the energy transfer. Moreover, different outlet blade angles and blade thickness distributions have strong influence on the energy transfer process. The losses are associated with the flow velocities and their growth is proportional to the velocity magnitude squared. As aforementioned, at the NGV domain, the flow has higher velocities close to the hub regions and its losses can be observed in Fig. 13b, in which high entropy generation can be observed.
The same analysis was performed for rotor blade row as presented in the Figs. 14 to 16. The static pressure distribution along the blade chord at suction and pressure sides, Mach number and entropy contours are presented. The same spanwise percentages were used: 20, 50 and 80%. Figure 14a shows a good variation of static pressure at blade surfaces for 20% of the rotor spanwise.
It is possible to observe that the Mach number is lower than unity and the turbine is unchoked at the design-point. This is a preliminary design strategy to maintain a certain design flexibility during the preliminary analysis. Sometimes, to obtain more power or pressure ratio, it is necessary to handle some possibilities. Hence, it is recommended that the design parameters not start with the limit values. Unlike NGV, at rotor hub, the velocity magnitude is lower than rotor tip values due to the tip speed. Figure 14b shows a good Mach number behavior of the flow along the rotor passage. Figure 14c shows good characteristics around the entropy evaluation. There is no separation or serious problems close to the rotor trailing edge at 20% of rotor spanwise. Figure 15 present the results for 50% of rotor spanwise. The flow velocity across the rotor passage for this position is greater than 20% of rotor spanwise. Mach number reaches the unity. At this position the flow is choked and a high energy transfer from fluid to rotor is achieved.
Though the value of Mach number for 50% of rotor spanwise reaches the unity, the behavior of the flow for this blade position is acceptable. The entropy increases only at rotor trailing edge. This is normal condition due to the flow mixture in this region.
The characteristic of the static pressure at 80% of the rotor spanwise is also acceptable, as shown in Fig. 16a. High Mach number can be observed at 80% of the rotor spanwise; this is due to the high tip speed values at this position, and supersonic flow can be reached as presented in Fig. 16b. In the present work, the tip Mach number reached the value of 1.2. Some geometrical modifications at rotor blade in the suction side, around at 80% of rotor spanwise, can be necessary to improve the flow characteristics in this region. The wake formation can be minimized changing the blade outlet angle and its thickness.
The entropy contours verified in Fig. 16c is due to the mixing region. In general, the results from CFD simulations at 80% of the rotor spanwise are acceptable.
For the suction side, as at stator, the results show satisfactory behavior for the flow property distributions. At rotor, there is a larger difference in flow behavior between the three radial positions, as shown in Figs. 14 to 16. Near the hub, the flow is almost subsonic and the parameter variations are softly and small. This is a consequence of the use of low degree of reaction at this region.
Considering the analysis at hub, the flow properties increase along the surface, becoming very smooth until approximately 60% of the chord, where it starts to suffer influence from the suction side, as also occurs at NGV row. However, for the rotor, this variation increases until the trailing edge, what is a consequence of the smaller blade-to-blade passage due to the higher number of blades. The highest Mach number value, approximately 1.0, is obtained close to the tip. All these results are in agreement with the prediction described by Silva (2012).
At blade suction side, mainly close to the tip, flow properties present great variation until the throat region (approximately 60% of the chord). At this point, the flow properties start to oscillate, with a sequence of deceleration, acceleration and deceleration again, as can be seen in Figs. 16a and 16b. The maximum Mach number found in both accelerations does not exceed 1.2. These results are also in agreement with the predictions made by Silva (2012).
In Fig. 17a it is possible to see the influence of the tip clearance on the rotor flowfield. The radial variation of the pressure and the velocity in this location are high. The clearance must be evaluated with a very careful process due to the high machine degradation related with the leakage flow. This leakage flow creates a secondary flow that interacts with the main flow. This mixing process generates high entropy and energy dissipation. Figure 17a shows the streamlines in this region. It is very common the structures formation, as scrapping vortex, that generally causes the drop in turbine efficiency. Figure 17b shows the entropy contours at rotor row outlet with interesting vertical formations at rotor tip clearance region. In the same figure it is possible to observe the secondary flow    growing and acting in the main flow, close to the hub, and creating at rotor leading edge horseshoe vortex. The effects of tip clearance in axial flow turbines that operate with high pressure ratios are severe for efficiency and gas expansion process. Several studies have been conducted to improve turbine performance decreasing the effects of tip clearance (Silva and Tomita, 2013;Ameri et al., 1998;Tallman, 2004;Lee and Choi, 2010;Tallman and Lakshminarayana, 2001;Azad et al., 2000;Bunker, 2014;Thulin et al., 1982;Sang and Byoung, 2008;Krishnababu et al., 2009). The results shown in the previous figures present a satisfactory flowfield within the axial flow turbine domain. The power turbine stage performance can be evaluated through the parameters calculated from the CFD techniques, as shown in Table 4. For efficiency, power and mass flow rate in parenthesis, the AxSTREAM® predicted different values and are also shown the error.
Generally, 1-D methods use loss modeling developed for reduced-order design tools, as meanline. These models should be calibrated for each design purpose. Some three-dimensional effects are not possible to account using 1-D tools, but, in general, the results agree for both techniques.

OFF-DESiGN POiNT ANAlYSiS
The off-design performance analysis was conducted only at the nominal speed of 27,920 rpm for comparison with the predicted results, varying the expansion ratio during the simulations. The results are presented in Fig 18 together with the preliminary study and its respective values. Figure 18a shows the variation of mass-flow with expansion ratio, in which the results from meanline and 3-D CFD are quite different. The CFD results show that the turbine choke condition is close to 6.85 kg/s of mass flow. This value is lower than that the choke condition determined in the meanline tool. This is due to the boundary-layer growth along the turbine streamwise direction that is, generally, underestimated by reduced-order tools. With high boundary-layer thickness, the velocity in the main flow increases, causing choking conditions. Meanline tool generally uses very simple models to account the blockage along the streamwise directions and in some tools this value is estimated by the designer. The expansion ratio for both tools at design-point presented slightly values with a difference around 0.11. Based on the CFD results, the turbine is choked at designpoint conditions. This is normal and is in agreement with turbine design methodology and means that high energy transfer occurs in the gas expansion process.
It is clear that a calibration of some coefficients used in the loss modeling should be adjusted to improve the turbine prediction for performance calculations. This is a very common action into the turbomachines design, in which, for several reasons, these  (c) (b) (a) calibration issues are considered as proprietary information for each industry. Figure 18b shows that the peak efficiency at design-point operation condition is slightly different for both tools with a difference of 0.6, in which the efficiency from CFD solution is lower. This is because the loss mechanisms determined by meanline are underpredicted. Generally, 3-D CFD considering viscous flows presents more accurate results due to the better representation of boundary-layers and fluid friction. The efficiency determined by CFD tool is calculated based on mass-averaging process. Figure 18c shows a good agreement between both numerical tools until the expansion rate of 1.98. The power determined by CFD tool is lower at design-point (around 160 kW lower than meanline tool). This can be explained due to the internal losses underpredicted by meanline tool. In general, all obtained results were well represented using the meanline tools, that are very fast to perform the calculations, comparing CPU time to calculate the 3-D turbulent flow by CFD and simple correlations.

CONCLUSION
In this study an in-house code written in FORTRAN language Silva (2012) was used to supply input data for commercial turbomachinery design software AxSTREAM® and to determine a preliminary design of a free power turbine. With the use of AxSTREAM® software the turbine NGV and rotor blades were generated by an airfoil stacking procedure and the 3-D geometry was obtained. This 3-D geometry was carefully treated to start the mesh generation process using the software ANSYS TurboGrid®. Four meshes with different number of control volumes were generated to study the mesh dependence from CFD solver. The CFD solver used in this work was the ANSYS CFX® v.13. Based on the mesh dependence study the mesh R was chosen to perform all other simulations of the turbine design-point and off-design point operational conditions.
The 3-D turbulent flow within the turbine computational domain was performed using RANS method applied on the fluid mechanics equations (continuity, momentum and energy). The numerical issues and boundary conditions were described in details. The simulation results from meanline and CFD tools were compared at design and off-design conditions. For both, the results are physically correct and are in agreement with the behavior of turbine flowfield.
Improvements in blade profile can be performed to increase the turbine efficiency decreasing the internal losses as shown from CFD results. Future work can be developed to calibrate the loss model to be applied in this free power turbine design. The CFD tool can be used aiming at improvements in the loss modelling calibration process.
The numerical tools and the methodology used in the present work can be used for axial flow turbine design purposes.