Multi-Objective Optimization of an Axial Flow Turbine Design Using Surrogate Modeling and Genetic Algorithm

Turbines are used as power generation devices in the jet propulsion and space technology due to its high power-to-weight ratio and reliability. The widespread use of turbines prompted researchers to investigate different methods for improving their performance such as reducing energy wastage and green gas emissions. The performance of a turbine depends on thermal as well as geometrical factors. Studying the performance change by varying the geometrical factors is difficult due to the presence of components consisting of complex shapes, such as stationary and rotating blade profiles, and small axial and radial clearances.

Some of the pertinent literatures which have addressed the influence of geometrical parameters on efficiency are briefly discussed here. References [1,2] numerically investigated the effect of axial gap and found that an optimum axial gap improves turbine efficiency due to better aerodynamic performance. References [3–5] have investigated the effect of tip clearance between rotating blade and stationary casing on the efficiency of turbines. Experimental studies in Ref. [6] have shown that effective generation of vortex near suction surface reduces the tip clearance leakage. References [7,8] numerically investigated the combined effect of geometrical parameters and identified hub, shroud profile, and tip clearance as the main influential parameters affecting turbine efficiency.

The major drawback of turbomachinery optimization employing computational fluid dynamics (CFD) solver is the time required for solving the large number of computationally intensive simulations corresponding to a large number of input geometry and thermal variable combinations that are required by the algorithms employed for finding the global optimum [9]. In some of the investigations, this difficulty is overcome by (i) solving simplified models involving assumptions such as irrotationality and incompressibility [10] and (ii) developing regression (surrogate) models from a manageable number of numerical simulations selected using sampling techniques [11,12]. References [13,14] used the Kriging-based surrogate model for multi-objective optimization of turbomachinery and nozzle flows. These studies have demonstrated that Kriging with the use of spatial correlation functions for generating the model is able to provide higher accuracy for multi-objective optimization problems.

Literature review shows that optimization problems are often multidisciplinary in nature and there is a need for considering more than one objective function simultaneously. The objective functions in multi-objective problems may not simultaneously attain maximum for the same factor vector and there can be many solutions called tradeoff solutions with better performances and selecting one from those designs involves a tradeoff. References [15,16] have reported different methods for finding the optimal solutions for multi-objective optimization problems. A method for multi-objective optimization considers Pareto-optimal front consisting of many optimal solutions known as Pareto-optimal solutions.

This paper presents a surrogate model-based multi-objective optimization methodology using Kriging and genetic algorithm for optimizing the power and efficiency of a gas turbine. Here, the objective functions are to maximize the power output and the efficiency. The optimization methodology consists of (i) validation of the CFD model using an experimental data, (ii) sensitivity analysis to identify significant parameters affecting power and efficiency, (iii) the selection of design space-filling samples using Latin hypercube experimental design, (iv) genetic algorithm-based optimization of parameters of a radial basis function which acts as surrogate (Kriging) model, (v) the generation of Pareto-optimal front for power and efficiency using a nondominated sorting genetic algorithm (NSGA-II), and (vi) CFD analysis of the optimized turbine configuration. The details about the methodology, methods, the optimum design parameter, and the flow behavior are presented and discussed in this paper.

A schematic diagram of the two-stage impulse axial flow turbine consisting of nozzle, rotors, and stator is shown in Fig. 1. There are 47 full admission subsonic nozzles placed at a pitch circle diameter (D_m) of 220 mm. The first and second stage rotors having 73 and 69 blades convert the kinetic energy of the high-velocity fluid leaving the nozzles to useful work. The stator consisting of 72 blades enables guiding the flow between the first rotor and the second rotor. The values of various parameters of the basic turbine design are shown in Table 1. The geometric variables varied for observing the effect are the axial gap between stationary and rotating blades (γ), the tip clearance between rotating blade and stationary casing (τ), and the stagger angles of the blade profiles; β_n, β_fr, β_s, β_sr as shown in Fig. 1.

The CFD analysis for the basic design configuration is carried out with geometrical parameters in Table 1. The blades are straight and solid modeled by a constant blade profile in the radial direction over a distance equal to the height of the blade. Mesh for single flow passage in the nozzle, first rotor, stator, and second rotor are generated using ANSYS-Turbogrid consisting of 250,000 elements for the nozzle and first rotor flow passage. The number of elements in the stator and second rotor flow passage are nearly equal and is 210,000. The CFD analysis is carried out using ANSYS-CFX by employing the shear stress transport (SST) k–ω turbulence model and stage interface between rotating and stationary domains. A summary of the results of the CFD analysis is shown in Table 2. The comparison between CFD and experimental results shows a good match with an error of 6.21% for power and 0.55% for efficiency.

The CFD evaluations to determine the number of computational cells required for grid converged solution are obtained by adopting the stepwise strategy described in Ref. [18]. Three grid systems consisting of 920,000, 2,830,000, and 9,950,000 elements have been developed for this. Efficiency of the turbine (η) is selected for observing the effect of grid size. Parameters for grid independence as described in Ref. [18] have been computed and tabulated in Table 3. The values of approximate relative error, extrapolated relative error, and the grid convergence index obtained from the computations are 0.14%, 0.26%, and 0.35%, respectively. These low values of percentage error show the insensitivity of the fine grid on the solution. All the results reported for the design configuration in this paper are those obtained from a grid system consisting of 920,000 cells.

Uncertainty in the experimental results of power output and efficiency are evaluated by using the least count of instruments used for observing the measured quantities and by applying the error propagation rules of different mathematical operations in the equations for power output and efficiency. The computed values of errors in the experimental values of different quantities are tabulated in Table 4. It can be seen that the percentage error in the values of power output and efficiency is 0.70% and 1.14%, respectively.

In this section, studies are conducted to identify the influence of geometric variables on the objective functions that maximize power output (P) and efficiency (η) of the turbine. The range of values taken for axial gap and tip clearance is γ ∈ [1.2–14 mm] and τ ∈ [0.25–2.25 mm]. Stagger angle is changed by the rotation of blade profile with respect to a pivot point at the leading edge as shown in Fig. 2. The rotation in clockwise and anticlockwise directions is specified as “+x” and “−x,” respectively. The profiles obtained after clockwise and anticlockwise directions are denoted as C_x and AC_x, where “x” represents the angle of rotation. A fortran-based program with inputs as profile co-ordinate for design configuration and angle of rotation is used for generating the rotated blade profile co-ordinates. The angle of rotation for nozzle, first rotor, stator, and second rotor profiles is denoted by variables Δβ_n, Δβ_fr, Δβ_s, and Δβ_sr, respectively. Ranges for the variables are Δβ_n ∈ [−6 deg to 1 deg], Δβ_fr ∈ [−3 deg to 6 deg], Δβ_s ∈ [−4 deg to 2 deg], and Δβ_sr ∈ [−6 deg to 4 deg]. A total of 49 number of CFD simulations were carried out with values of variables tabulated in Table 5.

Power and efficiency obtained from CFD simulations have been shown in Fig. 3. The lowest value of the power and efficiency of turbine to meet the functional requirement is 130 kW and 50%, respectively. Figure 3 shows that the range of geometric variable satisfying the lowest requirements is γ ∈ [2.3–14 mm], τ ∈ [0.25–2 mm]. Δβ_n ∈ [−2 deg to 0 deg], Δβ_fr ∈ [−3 deg to 6 deg], Δβ_s ∈ [−4 deg to 2 deg], and Δβ_sr ∈ [−6 deg to 4 deg]. The sensitivity of variables is compared based on percentage variance (%VAR) defined as the ratio of the normalized variance of each variable to the total normalized variance of all variables. The %VAR for Δβ_n, τ, and γ as per Table 6 are 82.601, 37.904, and 4.195, respectively. We consider that these values are significant and conclude that these variables are sensitive to the functional requirements.

Prediction models for system performance are developed using data points generated either theoretically or experimentally. The accuracy of the models depends on the space-filling characteristics of the input points. The present study uses a random Latin hypercube design for selecting “n” numbers of experimental points having a maximum spread in the design space. The number of variables defining a point in the design space (k) is three, viz., Δβ_n, γ, and τ. The initial data points with n = 100 are generated with a random permutation and evolve through an iterating process to maximize the space-filling characteristics of the design expressed using a scalar point function. The variables are normalized to convert them in the range [0, 1] such that design points can be generated from a cube [0, 1]³. For the Latin hypercube design, developed with n = 100 and k = 3 which is shown in Fig. 4, we use the python libraries in Ref. [19]. The technique and algorithm for developing the libraries can be found in Ref. [20].

This is also called the likelihood of the parameter for the given set of data, where

I is n × 1 unit vector and µ is the mean of observed responses, and Cov(y, y) = σ²Ψ where “σ” is the standard deviation.

Genetic algorithm is used to determine values of θ_j and p_j for maximizing the log of likelihood function which is expressed as log(L). Maximizing the log-likelihood function is equivalent to minimizing the negative of the log-likelihood function. In general, Kriging parameters θ_j and p_j are estimated as a minimization problem by minimizing the negative of the log-likelihood function. PyKrige (Version 1.5.0) is employed in this work for generating Kriging models for power and efficiency which uses a genetic algorithm with θ_j ∈ [0–10] and p_j ∈ [1.5–2]. Dataset to train Kriging models for power and efficiency consists of 90 points obtained from the sampling plan. Remaining ten points are kept unseen for model development and used for testing the model. Table 7 shows optimum value of θ_j and p_j for predictor variables Δβ_n, γ, and τ of the Kriging models for power and efficiency. Convergence of fitness function −log(L) with generations of genetic algorithm for Kriging models of efficiency and power is shown in Figs. 5(a) and 5(b), respectively.

Comparison of results of ten data (CFD) and Kriging model for the output validation of power and efficiency is shown in Figs. 6 and 7, respectively. The figures show a 45 deg line and the deviation of Kriging from CFD results can be observed as the deviation from slope = 1 line. R-square value for the model of power and efficiency is 0.995672 and 0.966793, respectively. It is concluded that models have good prediction capability.

An NSGA-II for the multi-objective maximization of two variables, viz., power output and efficiency, is used in this study. Since both the power output and the efficiency may not simultaneously attain the maximum for the same factor vector, a tradeoff becomes essential for the optimization. NSGA provides a set of solutions in which a better tradeoff is obtained between the competing objectives of maximum power output and efficiency. NSGA operates on the principle of elitism and diversity preservation using nondominated sorting and crowding distance strategies [22]. In nondominated sorting, solutions are classified as a nondominated front by considering the goodness of solutions determined using the dominance of fitness values of objective functions. F_x represents nondominated fronts where “x” is the rank of the front in descending order of dominance with x = 1, 2, 3…. For a set of solutions y_n, solution y₁ dominates over solution y₂, if y₁ is not worse than y₂ in all objectives or y₁ is strictly better than y₂ in at least one objective. Nondominated sorting generates a set of candidate solutions that are not dominated by any member of the solution set. In NSGA-II, a population of size “N” is generated in every iteration using the concept of nondominated sorting from a population of size “2N” obtained by combining parent and offspring populations of sizes N each. To start with, members in the F₁ front (rank 1) are selected for generating a new generation. Then, members in the front having the next lower ranks F₂, F₃,…, F_L are selected for forming the new generation. In case the nondominated front up to F_L contains more than N candidates, the selection is carried out using a crowding distance parameter for ensuring diversity of points in the objective space.

In the present problem with two objectives, the crowding distance of an ith point, c_i is defined using the values of power and efficiency at the neighboring points.

c_i = power at i + 1th point − power at i − 1th point + efficiency at i + point − power at i − 1th point − efficiency at i − 1th point

A very high value is assigned as the crowding distance for terminal points in the front. This ensures that the terminal points obtain maximum priority for selection to the new population. NSGA improves the goodness of solutions as iterations progress and converges to a set of optimal points with size “N” known as a Pareto-optimal set. For a biobjective problem, a set of nondominated points is manually generated from Pareto-optimal set using a nondominated sorting. These points are the set of optimal points from which a suitable configuration is selected based on the experience of the decision-maker.

In the present problem, optimization is carried out by controlling the nozzle rotation angle, axial gap, and tip clearance in the constrained ranges −2 deg ≤ Δβ_n ≤ 0 deg, 2.3 ≤ γ ≤ 14, and 0.25 ≤ τ ≤ 2, so as to maximize power output and efficiency of the turbine. Initial population with N = 100 (Fig. 8) for NSGA is generated using random sampling of the three design variables (Δβ_n, γ, τ) within the ranges. The power and efficiency obtained as output from the Kriging models serve as the fitness functions for NSGA. It can be seen that the predictor points of the Pareto-optimal set shown in Fig. 9 have been shifted to a higher power (above 156 kW) with the efficiency ranging from 49.85% to 61.19%. Point 1 in Fig. 9 corresponds to the basic design of the turbine. Points 2 and 3 correspond to the design points identified from the sensitivity analysis as having maximum power and efficiency, respectively. In order to select an optimum configuration, a set of nondominated points is identified from the Pareto-optimal set. The domain of values of the predictor variables (γ, τ, Δβ_n) corresponding to the objective functions (η, P) of the nondominated points is shown in Table 8. These points are joined to obtain a Pareto front as shown by a black line in Fig. 9. The obtained domain of predictor variables provides a certain level of flexibility in selecting an optimum configuration for the decision-maker. In the present work, we have selected a higher tip clearance (τ) equal to 1.25 mm with an efficiency of 55.28%. An increase in tip clearance by 66.7% compared to the original design reduces the chances of rubbing between rotor and casing. Other design variables are axial gap (γ) and nozzle rotation angle (Δβ_n) with values of 6.03 mm and −0.19 deg, respectively. The CFD analysis of the selected design has been found to predict a power output and efficiency of 164.99 kW and 53.01%, whereas the corresponding values from the Kriging model are 164.05 kW and 55.2%. It can be seen that the percentage differences are only 0.58 and 4.1, respectively.

As mentioned in earlier sections, the three design variables, viz., axial gap (γ), tip clearance (τ), and nozzle rotation angle (Δβ_n) have been varied in this work to determine an optimum configuration that gives better performance in terms of power output and efficiency. The values of design variables in the basic design configuration and the optimized configuration are shown in Table 9. The inflow boundary conditions used for the analysis for both the configurations are stagnation pressure equal to 5.6 bar and mass flowrate equal to 5.316 kg/s. The boundary conditions on the walls are no-slip for the velocity components and zero normal pressure gradient. The results along with a discussion on the comparison of flow in the two design configurations are given in the following sections. The reasons for the better performance in optimized design are also discussed.

The contour plot of velocity in the basic design configuration is shown in Fig. 10. A flow separation is observed at the suction side of the first rotor at about 20% of the axial chord. Ahead of the separation point, a sudden increase in velocity equal to 393 m/s can be seen with a corresponding Mach number equal to 1.4. This is not a desirable condition in a subsonic turbine. The supersonic flow instantaneously changes to subsonic flow through the formation of a shock across the flow path as indicated in Fig. 10. A sudden increase in temperature from 200 K to 235 K at the suction side of the first rotor is shown as Detail-A in Fig. 11 for temperature distribution at mid-span of nozzle-first rotor flow regime confirms the formation of a suction side shock. The shock initiates from the suction side and extends to the flow in the axial gap region between the nozzle and the first rotor. The total pressure loss across the shock wave is found to be 0.8 bar and this local pressure drop is 14.3% of total pressure at the nozzle inlet. Based on the values of local pressure in this region, it can be concluded that flow separation and shock contribute to major losses in the flow path of basic design configuration.

The variation of the magnitude of relative velocity in the nozzle-first rotor axial gap region at the mid-span for one pitch of blade flow passage is shown in Fig. 12. It can be seen that the velocity increases from 123 m/s downstream of the blades called the wake region to 282 m/s at the middle of the passage called the jet region. Jet and wake from the nozzle are expected to mix within the axial gap region having a length of 3.5 mm to produce a uniform velocity at the first rotor inlet. Even though there is entropy generation as well as loss of availability associated with mixing, a certain optimum mixing is expected to develop a favorable rotor incidence angle that reduces the losses of energy in the rotor passages. The extent of mixing in the axial gap is assessed by comparing the magnitude of relative flow velocity at different axial locations within the axial gap as shown in Fig. 12. Locations 1 and 2 are at 2% and 50% of the axial gap from the nozzle trailing edge. The maximum velocity deviations at locations 1 and 2 are 159 m/s and 131 m/s, respectively. This indicates that velocity fluctuations at nozzle exit have not reduced much due to an axial gap of 3.5 mm length. The range of first rotor incidence calculated based on relative flow velocity at location 2 is between −0.31 deg and 1.98 deg which can lead to higher incidence losses.

Figure 13 shows the stream tracer plot in the tip clearance region of the first rotor. The computed value of tip clearance leakage flow is 6.86 g/s. By judging based on the density of streamlines, maximum leakage flow can be observed near the leading and trailing edge portions. The leakage generates a tip vortex which gradually develops from the mid-chord of the blade. The vortex strengthens as it moves toward the suction side by absorbing energy from leak flow having a high velocity of 250 m/s. Mixing of strong tip vortex disturbs main flow near the suction surface as seen in Fig. 13. This hindrance to the main flow by the tip leakage flow can be considered as a major source of the energy losses.

Figure 14 shows the contour plot of the magnitude of relative velocity. It can be seen that the sudden increase in relative velocity in the suction side of the first rotor observed for basic design is not found in the optimized design. The peak velocity of 393 m/s in the basic design at about 20% of the axial chord has been reduced to 374 m/s in the optimized design with a Mach number of 1.3. The supersonic flow changes to subsonic flow through a shock. From the distribution of temperature shown as Detail-A in Fig. 15 for the nozzle-first rotor flow passage, it can be seen that the size of the flow regime affected by shock wave is reduced for the turbine with an optimized design configuration. The shock initiated from the first rotor suction side could not penetrate to the nozzle-first rotor axial gap region of the optimized configuration. The total pressure loss across the shock wave is found to be 0.3 bar, which is only 5.36% of total pressure at the nozzle inlet compared to 14.30% in the basic design.

The variation of the magnitude of relative velocity in the nozzle-first rotor axial gap region at the mid-span for one pitch of blade flow passage is shown in Fig. 16. The jet and wake velocity downstream of the nozzle blade is 294 m/s and 106 m/s, respectively. The extent of mixing for the higher axial gap of γ = 6.031 mm is shown in Fig. 16. The maximum velocity deviations at location 2 (50% axial gap) are found to be 19 m/s, which is much less than the corresponding values in the basic design (131 m/s). This has resulted in the variation of first rotor incidence at location 2 to reduce considerably. The variation in the optimized design is between 0.259 deg and 0.518 deg compared to −0.31 deg and 1.98 deg in the basic design. This has contributed to a reduction of losses in the rotor.

Figure 17 shows the stream tracer plot in the tip clearance region of the first rotor. Increase in tip clearance to 1.2485 mm has been found to shift the location of the tip clearance vortex from the suction side of the first rotor to the top of the rotor in the tip clearance region as shown in Fig. 17(a). Increase in tip clearance aids the formation of vortex ahead of the mid-chord of the first rotor. Stream tracer plots for tip clearance flow at three sections of the first rotor blade; near the leading edge, 50% of axial chord, and near the trailing edge are given in Figs. 17(b)–17(d), respectively. These figures clearly show the presence of a vortex extending from leading edge to trailing edge of the blade. This vortex acts as a flow barrier between pressure and suction surface thereby reducing the tip clearance leakage to 5.1 g/s from 6.86 g/s in the basic design. The lower tip clearance leakage has been found to reduce the flow disturbances in the main flow path.

Based on the discussion given in previous sections, the increase in tip clearance to 1.249 mm for the optimum design configuration generated an optimum sized tip vortex which reduced the tip clearance leakage by 25.66% compared to the basic design configuration. A discussion on the effect of further increase of tip clearance on the tip vortex formation and tip clearance leakage is given in this section. A higher tip clearance of 2.25 mm is taken for this study. A stream tracer plot in Fig. 18(a) for first rotor tip clearance shows that the vortex remains on top of the rotor, extending from leading edge to trailing edge, which is similar to the one observed in the optimized design configuration. The stream tracer plots at two sections of the first rotor, one near the leading edge and the other at 50% of the axial chord are shown in Figs. 18(b) and 18(c). These plots show that the size of the vortex is lower than tip clearance leading to leak paths through the gap above the vortex. The additional leak paths increased the tip clearance leakage flow to 13.1 g/s compared to 5.1 g/s in the optimized design configuration. The higher leak rate as well as increased flow disturbance in the main flow path at the suction side of the rotor resulted in a reduction in the efficiency and power output to 40.72% and 155.60 kW, respectively, for the design configuration with higher tip clearance gap of 2.25 mm.

A set of design parameters required for achieving higher efficiency and power output in axial flow gas turbines has been obtained in this paper by developing the Kriging (surrogate) model and carrying out multi-objective design optimization using NSGA-II. An optimum combination of axial gap and tip clearance is found to produce a flow structure with a reduced amount of energy losses for making higher efficiency and power output. The major consequences in flow are adequate mixing of jet and wake flow, reduction in the effect of shock, and formation of an optimum sized vortex in the tip clearance flow path. The geometric variables that can influence the efficiency and power of the turbine are identified using sensitivity analysis in which the axial gap (γ), tip clearance (τ), and rotation angle of nozzle profile (Δβ_n) are varied in the ranges γ ∈ [2.3–14 mm], τ ∈ [0.25–2 mm], and Δβ_n ∈ [−2 deg to 0 deg]. Accurate data from ANSYS-CFD simulations employing the SST k–ω turbulence model are generated for a set of space-filling sampling designs selected using a random Latin hypercube design of experiment. The CFD data are used as input for the development of surrogate models for the efficiency and power by using Kriging interpolation on a radial basis function. A set of Pareto-optimal solutions for the turbine are generated using NSGA-II in which power and efficiency were set for maximization. The CFD analysis is carried out for a configuration selected from the Pareto-optimal solutions to understand the fluid dynamic reasons for the improvements in performance. The increase in the axial gap of the optimized design configuration has been found to improve the mixing of jet and wake flow in the nozzle-first rotor axial gap region resulting in a proper incidence angle favoring a reduction in rotor incidence and passage flow losses. A reduction in shock strength is observed in the first rotor flow passage of optimized configuration which is found to reduce the total pressure loss across shock by 62.5%. Increase in the tip clearance is found to shift the location of the tip clearance vortex from the suction side of the first rotor (for basic design configuration) to the top of the rotor in the tip clearance region (for optimized configuration) resulting in a reduction of tip clearance leakage by 25.66%. The percentage increase in the efficiency and power of the optimized design compared to the basic design are 5.1% and 1.02%, respectively.

There are no conflicts of interest. This article does not include research in which human participants were involved. Informed consent not applicable. This article does not include any research in which animal participants were involved.

The authors attest that all data for this study are included in the paper.

k =

number of predictor variables

m =

mass flowrate, kg/s

n =

number of samples

x =

sample data

y =

response data

Η =

efficiency, %

I =

n × 1 unit vector

L =

likelihood function

N =

speed of turbine, rpm

P =

power output, kW

c_i =

crowding distance

p₀₁ =

total pressure at inlet, bar

p₁ =

pressure at inlet, bar

p₂ =

pressure at outlet, bar

D_m =

pitch circle diameter, mm

T₀₁ =

total temperature at inlet, K

β_fr =

stagger angle of first rotor blade, deg

β_n =

stagger angle of nozzle blade, deg

β_s =

stagger angle of stator blade, deg

β_sr =

stagger angle of second rotor blade, deg

γ =

axial gap, mm

Δβ_fr =

rotation angle of first rotor profile, deg

Δβ_n =

rotation angle of nozzle profile, deg

Δβ_s =

rotation angle of stator profile, deg

Δβ_sr =

rotation angle of second rotor profile, deg

σ =

standard deviation

τ =

tip clearance, mm

Ψ =

correlation matrix

ψⁱ =

radial basis function