OPTIMIZATION OF THE LS89 AXIAL TURBINE PROFILE USING A CAD AND ADJOINT BASED APPROACH

The LS89 high pressure axial turbine vane was originally designed and optimized for a downstream isentropic Mach number of 0.9. This proﬁle has been widely used for CFD validation in the open literature but very few attempts have been made to improve the already optimized design. This paper presents the results of a single point CAD-based adjoint CFD optimization of the LS89. The gradients used for the optimization are obtained via a novel approach in which both the CAD kernel and grid generation are differentiated using Algorithmic Differentiation techniques. A reduction of almost 12% entropy generation is achieved, which is equivalent to a 16% total pressure loss reduction. The entropy generation is reduced whilst keeping the exit ﬂow angle as a ﬂow constraint, which is enforced via the penalty formulation. The resulting unconstrained optimization problem is solved by the L-BFGS-B algorithm. The trailing edge thickness and the axial chord length are kept ﬁxed as manufacturing constraints. The ﬂow is governed by the Reynolds-averaged Navier-Stokes equations and the one-equation transport Spalart-Allmaras turbulence model. The optimal proﬁle is compared and benchmarked against the baseline case.

dJ d α -performance sensitivity vector dJ d X -adjoint sensitivity vector d X d α -grid sensitivity vector ϕ P S -pressure side trailing edge wedge angle ϕ SS -suction side trailing edge wedge angle σ -solidity θ -momentum thickness ζ 2 -downstream loss coefficient ζ p -profile losses

INTRODUCTION
Designing and optimizing a turbine vane is a complex iterative design process that can take significant time and effort. The use of relatively low-cost numerical shape optimization methods has become more popular and have been widely used in turbomachinery applications. In particular, the adjoint method enables the efficient computation of the sensitivities required by gradient-based optimizers, at a cost independent of the number of design variables (Peter and Dwight, 2010).
Typically, the grid point coordinates have been used as design variables (Jameson, 2004). This approach offers a very rich design space but has some drawbacks that need to be considered. First, the connection to the CAD geometry is lost. This means that an additional step is required to transform the optimal shape defined by grid points back to smooth CAD shape. Secondly, although there are several successful methods that can be used to convert a cloud of grid points to CAD (Becker et al., 2011;Leal et al., 2010), the fitting error can impair the optimality of the shape. Furthermore, this additional step can take significant time and it is not guaranteed that the final approximated CAD shape will meet the design requirements and constraints (Samareh, 2001;Braibant and Fleury, 1984). Also, the use of a smoother is inevitable to avoid irregular shapes. To avoid these problems, in this work it is proposed to keep the CAD geometry in the optimization loop.
Changing the design variables from the grid coordinates to CAD-based parameters, such as axial chord length, stagger angle, trailing edge wedge angle, etc., requires an additional sensitivity computation step. The grid sensitivities, which are defined as the partial derivative of the grid coordinates with respect to the design variables, will need to be computed as well. This could be done by using finite differences but the accuracy of the gradients is dependent upon the chosen step-size for each design variable. If the step-size is too small the gradients are inaccurate due to limited arithmetic precision, whereas too big a step-size would yield truncation errors. This work aims to circumvent these issues by making use of Algorithmic Differentiation (AD) (Griewank and Walther, 2008) for the CAD kernel and the grid generation (Sanchez Torreguitart et al., 2016). AD allows to compute the derivatives of an output of a program with respect to the inputs by applying the chain rule in an automatic fashion throughout the evaluation of the code. This means that the grid sensitivities will be accurate to machine accuracy.
In this work, it is proposed to use a CAD-based adjoint-based optimization method for the LS89 (Arts et al., 1990;Van den Braembussche et al., 1990) high pressure axial turbine nozzle guide vane profile. The LS89 was originally designed and optimized at the Von Karman Institute for Fluid Dynamics for a subsonic isentropic outlet Mach number of 0.9 by the inverse method (Van den Braembussche et al., 1990), which is an iterative design method based on both potential and Euler type solvers that uses the difference between the calculated velocity distribution and the required one to modify the profile geometry. Montanelli et al. (2015) performed both single-point and multi-point optimizations on the LS89 to reduce the total pressure losses by constraining the outlet mass flow, whilst keeping the leading and trailing edge geometries and the profile thicknesses fixed and using the L-BFGS-B algorithm to explore different geometries with a wing parametrization model. Also, they solved the Reynolds-Averaged Navier-Stokes (RANS) equations and the one-equation transport Spalart-Allmaras turbulence model and computed the gradients assuming that the eddy viscosity and thermal conductivity are constant. Montanelli et al. (2015) achieved a total pressure loss reduction below 1% for the nominal case (outlet M ise = 0.927). The question arises whether the LS89 is indeed optimal, and cannot be optimized any further. This paper addresses this and presents an adjoint optimization of the turbine profile, by using a novel CAD-based approach in which: 1. The optimal shape remains defined within the CAD tool. The optimization problem herein is expressed by CAD parameters that are directly used in defining the CAD geometry by means of Bézier and B-spline curves.
2. The in-house CAD and grid generation tools are automatically differentiated in forward mode to obtain the exact derivatives of the grid coordinates with respect to the CAD-based design parameters. This allows to accurately predict the sensitivities and circumvent the errors introduced by finite differences.
3. The trailing edge thickness and axial chord length are kept fixed as manufacturing constraints and the exit flow angle is considered as an aerodynamic constraint.
A computer aided design and optimization tool for turbomachinery applications (CADO) (Verstraete, 2010) is used throughout this work.

CAD-BASED PARAMETRISATION
The construction of the turbine profile shown in Fig. 1(a) and 1(b), is based on the parametrization described by Pierret (1999).The profile is defined by a set of CAD-based or engineering based design parameters that are relevant to the aerodynamic performance (e.g. solidity, stagger angle, etc.) and the manufacturing requirements (e.g. axial chord length, trailing edge radius). First, a camber line is constructed. The points P LE ,P mid , P T E define the control points of the 2nd order Bézier curve describing the camber line. The camber line is used to define the position of the control points of the suction side (SS) and pressure side (PS) B-spline curves relative to it. The profile is constructed by two B-spline curves for the SS and PS and a circular arc at the TE to close the profile. The normal distance of the first control point relative to the camber line is calculated in such a way that G2 geometric continuity (i.e. equal curvature) is maintained between the SS and PS B-splines at the leading edge.

OPTIMIZATION
The purpose of the optimization algorithm ( Figure 2) is to reduce the entropy generation J 1 (Equation 1) whilst maintaining the exit flow angle J 2 (Equation 2) above or equal to the baseline value of J 2,ref = 74.83 deg, by modifying the design vector α. The flow state U and the design vector α are coupled via the primal flow governing equations. The constrained optimization problem is handled with the penalty method. The penalty term becomes active only when the exit flow angle is below the target value. The J 2,f lag parameter is either set to 1.0 or 0.0 in order to activate or deactivate the penalty term respectively. The weighed cost function of the single point J SP optimization problem is given by Equation 3. The L-BFGS-B algorithm (Zhu et al., 1997) available in the python SciPy package (Jones et al., 2001) is used to find the new design vector.

Grid generation
A multi-block structured grid is generated for every optimization iteration. A meshindependence study was carried out in order to select an appropriate mesh for the optimization. One O-grid block is placed around the profile and there are six additional H-grid blocks, four of them distributed around the profile and the remaining two being used as the inlet and outlet blocks. Sanchez Torreguitart et al. (2016) describes in more detail how the grid was generated.

Flow and Adjoint solvers
The flow solver employs a cell-centered finite volume discretization on multiblock structured grids. The three dimensional compressible RANS equations are solved with an implicit time integration scheme accelerated by local time-stepping and multigrid. The fluid is considered as a calorically perfect gas and the eddy-viscosity hypothesis is used to account for the effect of turbulence. Convective fluxes are computed using Roe's approximate Riemann solver (Roe P.L., 1981) with a MUSCL-type reconstruction (Van Leer, 1979) of primitive variables to attain second order accuracy. Oscillations near shocks are suppressed by a van-Albada type limiter (Venkatakrishnan V., 1993). The numerical dissipation of the scheme is controlled by the entropy correction by Harten and Hyman (Harten A., Hyman J.M., 1983) for both linear and non-linear eigenvalues. Viscous fluxes are calculated with a central discretization scheme, while the negative Spalart-Allmaras turbulence model (Allmaras S.R.; Johnson F.T.; Spalart P.R., 2012) is used for the turbulence closure problem assuming fully turbulent flow from the inlet (Re inlet ≈ 2 · 10 5 ). Boundary conditions are imposed weakly by utilizing the dummy cell concept Blazek J. (2001). The hand derived discrete adjoint solver uses constant eddy viscosity assumption which is a valid approach for engineering design applications. The flow solver and its discrete adjoint counterpart use the stabilization scheme described by Xu et al. (2015). Sanchez Torreguitart et al. (2016) showed the first results of differentiating the CAD kernel and grid generation tools using the AD tool ADOL-C (Walther and Griewank, 2014). In the work herein, the optimization algorithm computes the grid sensitivities d X/d α for every optimization step, using the forward vector mode as opposed to the forward scalar or reverse modes as it allows to compute d X/d α in one single evaluation of the primal at a relatively low cost. The performance sensitivities dJ/d α for each cost function J 1 or J 3 are computed by do- ing a scalar product of the adjoint calculated sensitivities dJ/d X (ie. the derivative of the cost function with respect to the grid coordinates X) with the grid sensitivities d X/d α. A suitable step-size for each design variable was selected to compute the gradients with finite differences and compare them against the gradients obtained by the discrete adjoint. Figure 3(a) and 3(b) show a good agreement with finite differences for both the entropy generation and exit flow angle. The largest sensitivities shown in Figure 3(a) and 3(b) correspond to the design variables α j = 1 (c ax ) and 4 (R T E ), which do not change during the optimization because they are considered as manufacturing constraints. Figure 4(a) shows significant shape differences between the baseline and the optimal shape, the later being a more aft-loaded profile than the baseline design. The isentropic Mach number distributions are shown in Figure 4(b), which shows a fairly good agreement between the M ise predicted by CFD and the experimental data (MUR45 test condition, with downstream M ise = 0.875, see (Arts et al., 1990)) for the baseline geometry.

RESULTS
The L-BFGS-B algorithm converges within 20 optimization iterations. More cycles were performed but no further decrease in the objective function was obtained. Despite the fact that the mass flow for the optimal geometry has increased slightly by 1.2%, the entropy generation is reduced by 11.6%, whilst maintaining the exit flow angle within ±0.1% of the baseline value. As the aft-loaded profile has significant rear suction side curvature, the flow is rapidly accelerated towards a higher peak Mach number than the baseline design. This is followed by a deceleration from X/c ax = 0.8 to 0.88 and a small acceleration from X/c ax = 0.88 to 0.91. At X/c ax = 0.91 the flow is rapidly decelerated again and leaves the trailing edge at a significantly lower exit isentropic Mach number. The strong deceleration towards the end of the suction side suggests that the boundary layer will be thicker and one would expect higher profile losses. However, the pressure losses figures shown in Table 3 show the opposite. The total pressure loss reduction between the inlet and the fully mixed-out flow at the outlet plane, expressed via P 01 − P 02 or by the ζ 2 coefficient (see Equation 4), is of the order of 16.3%. The downstream total pressure loss profile, at the plane X/c ax = 1.433, is shown in Figure 5. In the vertical axis the total pressure loss P 01 − P 0X between the inlet and a downstream plane at x/c ax = 1.433 is expressed as a a percentage of the downstream dynamic head q X at the X/c ax = 1.433 plane. This figure shows that the total pressure loss generated in the wake is reduced significantly.
In order to understand why the optimal shape is more efficient than the baseline, the following subsections will look into the loading, boundary layer parameters, base pressure and profile losses for each profile.
Zweifel loading coefficient The total loading, which can be expressed as the integral of the p P S − p SS over the axial chord length, has increased by 1% for the optimal shape as a result of the 1% increase in the pitch and having the same exit flow angle. The Zweifel loading coefficient Z w (Zweifel, 1945) is a measure of how efficiently the profile is loaded. An efficient loaded profile would be one in which the pressure is stagnated over the pressure side and the velocity on the suction side is constant and equal to the downstream velocity. The Zweifel's design rule says that loss is minimized for 0.8 < Z w < 1. Equation 5 was used to calculate the Zweifel loading coefficient, which increased from 0.636 for the baseline to 0.646 for the optimal.

Boundary layer parameters
The state of the boundary layer close to the trailing edge is deemed to have an impact on the profile losses (see Equation 6). Table 1 compares the boundary layer parameters between the baseline and optimal shapes for two locations on the suction side: at the peak M ise location and very close to the TE (approximately at x = 0.03m and x = 0.037m respectively, see Figure  4(b)). The later location is defined herein as the point of the suction side slightly upstream of the TE circle. The boundary layer velocity profiles are shown in Figures 6(a) and 6(b) for both locations, showing that there is no boundary layer separation. At the peak M ise location, the boundary layer thickness δ for the optimal design is smaller than for the baseline, as expected due to the higher acceleration. When the flow decelerates from the peak M ise location to the point just upstream of the trailing edge circle, the boundary layer thickness almost doubles for the baseline and triples for the optimal profile. The resulting δ of the optimal shape is 18% thicker than the baseline. The shape factor H is higher for the optimal than for the baseline design in both locations, which means that the optimal profile has stronger adverse pressure gradient. The shape factor for both designs are between the typical laminar Blasius boundary layer value (H = 2.59) and the turbulent values (H = 1.3 − 1.5). From a boundary layer perspective, the optimal design would have higher profile losses due to the higher thickness displacement and momentum thickness. However, there is at least one more factor to be considered when determining the profile losses: the base pressure.

Base pressure
The base pressure, which is the static pressure at the mid point of the trailing edge is known to play an important role in the profile losses. Higher base pressures contribute to the reduction of the profile losses. The trailing edge pressure distribution shown in Figure 7(a) shows that  Table 1: Boundary layer parameters the optimal shape has got 10% higher base pressure. It is generally thought that higher values of boundary layer thickness result in higher values of base pressure (Xu and Denton, 1988). Also, higher base pressures are expected of an aft-loaded profile with significant rear suction side curvature and higher values of unguided or uncovered turning (Sieverding et al., 1979). The angle between tangents to blade suction side in the throat and at the TE increased by 6.76 degrees when going from the baseline to the aft-loaded optimal profile. According to the base pressure correlation given by Sieverding et al. (1979), the optimal profile is expected to have 3% higher base pressure than the baseline. The percentage increase predicted by the correlation is mainly due to the higher suction side rear curvature because the trailing edge wedge angle difference between the baseline and optimal profiles is negligible. However, Sieverding et al. (1979) showed that, for 80% of the tested blades, there was a discrepancy between the predicted base pressures (using the correlation) and the measured ones of the order of +/-5% of the downstream pressure. Although there are no measurements available for the optimal profile yet, by taking this error into account the minimum and maximum values for the base pressure of the optimal profile can be estimated to be 5% lower and 16.7% higher than the baseline respectively. Despite the 10% higher base pressure predicted by RANS falls inside this range and is in line with what would be expected of an aft-loaded profile with significantly higher rear suction side curvature (Sieverding et al., 1979), only unsteady analysis and experimental tests will be able to provide enough evidence to support these results.   Sieverding et al. (2003). LES results are from (Vagnoli et al., 2015) Figure 7: Trailing edge pressure distributions Predicting the absolute values of the base pressure correctly can be very difficult due to the highly complex nature of the flow in the trailing edge. Vagnoli et al. (2015) showed that steady state simulations usually predict the wrong shape or distribution of pressure around the trailing edge, whereas LES or URANS capture much better the shape and the pressure values when compared to experimental data. The pressure distribution shown in Figure 7(a) is a nearlyuniform pressure around the trailing edge, which is suspected to be different to the real pressure distribution profile. However, it is believed that the base pressure delta difference between two steady state simulations, one for the optimal and the other for the baseline profiles, is a fairly good estimate of the real difference. The validity of this assumption has been investigated by performing two RANS simulations for the turbine profile that was tested by Sieverding et al. (2003) at M ise = 0.79 and 0.97 and comparing the base pressure delta between these two operating conditions against the experimental data and the LES results performed by Vagnoli et al. (2015) (Figure 7(b)). The base pressure delta, which is taken at S/D T E = 0 from Figure  7(b), is 0.18 for both the experiments and RANS and 0.20 for LES. This confirms that, despite the limitations of CFD in predicting the base pressure absolute values, RANS is able to capture very well the relative differences in base pressure between two operating conditions. Denton (1993) showed that the profile losses for incompressible flow can be related to the base pressure coefficient and the boundary layer parameters, as shown in Equation 6. The relative effect of the base pressure coefficient is increased in compressible flow (Corriveau and Sjolander, 2003). For blades operating in the transonic range from 0.8 to 1.2 the major source of loss is the trailing edge loss at these speeds and the base pressure plays a dominant role in determining the loss (Xu and Denton, 1988). The profile loss split is shown in Table 2. The higher δ * and θ values of the optimal shape contribute to increasing the profile loss. However, the reduction in profile losses due to the 302.2% increase in the trailing edge loss term dominates and explains the overall 30% reduction in profile losses. However, if the base pressure increase was equal or below 6.54%, the optimal shape would have similar or higher profile losses than the baseline, respectively. The question arises whether using RANS for the optimization of an axial turbine operating at transonic speeds is a valid approach, since the base pressure plays such a dominant effect on the profile losses for such speeds and the flow field around the TE of a transonic turbine blade is too complex to be captured by RANS. Future LES simulations and or experimental test are going to be necessary to support or reject the aerodynamic improvements reported herein for the optimal shape.   Table 3: Comparison between the baseline and optimal CONCLUSIONS This paper presents an single point optimization of the LS89 axial turbine cascade vane profile for the design downstream isentropic Mach number of 0.9. The results show that the total pressure losses and entropy generation can be reduced by 16% and nearly 12% respectively by going from a front to a rear loaded profile and by keeping the exit flow angle fixed. Despite the increase in mass flow, loading, boundary layer thickness displacement and momentum thickness, the base pressure coefficient plays a dominant role herein and reduces the profile losses by 30%. The substantial aerodynamic improvements reported herein have not been observed previously, probably due to the fact that the current CAD-based parametrization allowed the exploration of a richer design space. The relative delta in base pressure predicted by the steady simulation between the optimal and baseline profiles is deemed to be representative of the real value, as RANS tends to capture well the relative differences between two designs. However, future unsteady studies or experimental investigations would be beneficial in order to confirm these findings.