Shape optimization of a bidirectional impulse turbine via surrogate models

ABSTRACT In this work, a bidirectional impulse turbine used in a wave energy device is simulated using a CFD technique and a shape optimization performed with multiple surrogates-assisted multi-objective evolutionary algorithm. The surrogates have been generated using the validated CFD technique. The objectives for the optimization are to maximize the shaft power and to minimize the pressure drop across the turbine. The hub and the tip thickness of the rotor blade profiles are modified and considered as the design variables. The Latin hypercube sampling technique generated design points are evaluated in the CFD solver, the objective responses are used to construct multiple surrogates, and a Pareto optimal front is generated through a multi-objective optimization approach. The turbine efficiency, which is the function of pressure drop and shaft power, is relatively increased by 10.4%. Abbreviations: BPS: best PRESS; GV: Guide vane; KRG: kriging; LE: leading edge; MOO: multi-objective optimization; OWC: oscillating water column; PRESS: predicted error sum of squares; PS: pressure side or pressure surface; RB: rotor blade; RBNN: radial basis neural network; Ref: reference; RMS: root mean square; RSA: response surface approximation; SS: suction side or suction surface; TE: trailing edge; WAS: weighted average surrogates


Introduction
Ocean wave generated by the action of wind over the water surface is a source of enormous renewable energy. Harnessing wave energy is difficult because of uncertainty in ocean waves, natural calamities, etc. The oscillating water column (OWC) (Figure 1), a widely used device to extract the power from the waves, has a cylindrical column in which wave rises and replicates the action of piston. This pressurizes the trapped air in the cylinder, and the air drives a bidirectional turbine. The airflow is reciprocating in nature. The turbine rotates in one direction irrespective of the airflow direction (Setoguchi, Santhakumar, Maeda, Takao, & Kaneko, 2001). Hence, it is called a bidirectional turbine or bidirectional flow turbine. There are several types of bidirectional turbines, and impulse turbine is one among them.
Usually a bidirectional turbine consists of a couple of guide-vane (GV) rows and a rotating blade row. One GV row deflects the flow toward the rotor blade (RB) inlet while the other directs the RB exit flow. An impulse turbine with GVs gives better performance than a turbine without GV (Maeda et al., 1999). Several works are done for optimization to modify parameters to improve the CONTACT Abdus Samad samad@iitm.ac.in performance. The parameters include blade thickness by camber line modification (Gomes, Henriques, Gato, & Falcão, 2012), number of RBs and GVs (Badhurshah & Samad, 2015), tip clearance (Thakker & Dhanasekaran, 2004), RB solidity (Cui & Liu, 2014), blade setting angle (Liu, Cui, Kim, & Shi, 2016), guide vane shape (Thakker, Dhanasekaran, & Ryan, 2005), hub-to-tip ratio (Thakker, Khaleeq, Takao, & Setoguchi, 2003), etc. for the impulse turbine. Gomes et al. (2012) optimized the RB thickness distribution by modifying the camber line and they generated the stacking optimized sections along the span. Badhurshah and Samad (2015) optimized the number of RBs and GVs of the turbine for single objective and found that the turbine with 38 rotor blades and 24 guide vanes as the optimal design. The RB with 0.25% tip clearance performs similar to the turbine without any tip clearance (Thakker & Dhanasekaran, 2004). Regarding the performance improvement of a turbomachinery blade, there are several other parameters, which were tried in other types of turbines, and the thickness is one of the important parameters. Takasaki, Takao, and Setoguchi (2014) discussed the effect of blade shape on the performance of a reaction turbine, and found that the blade with linearly varying thickness from the mean radius to the tip gives a higher efficiency, but it does not overcome the weakness of stalling. Samad, Kim, Goel, Haftka, and Shyy (2008) found that the adiabatic efficiency and pressure ratio of the compressor increases if stacking and thickness is modified in a compressor. Sarraf, Nouri, Ravelet, and Bakir (2011) found that the nominal point is shifted toward lower flow rate for the thick blade and the aerodynamic characteristics of an axial fan are steeper and the pressure fluctuations are lower. The flow at the exit is axial for a wide range of flow coefficients for thick blades and the efficiency is lower than thin blades. The thickness of NACA00xx increases the boundary layer thickness at low incidence angle (Sarraf, Djeridi, Prothin, & Billard, 2010). A thicker blade has more trailing edge (TE) drag and more mixing loss (Roelke & Haas, 1983).
Any engineering optimization problem has multiple variables and multiple objectives along with or without design constraints. The surrogate-based techniques help reducing evaluation of multiple designs in any CFD or experimental method. Multiple surrogate models are used because a single surrogate may give a higher approximation error (Viana, Simpson, Balabanov, & Toropov, 2014). Wu, Chau, and Li (2009) used the multiple data preprocessing techniques to improve neural network performance and they found that the neural network with moving average is optimal for flow predictions. The surrogates such as response surface approximation (RSA), weighted average surrogate (WAS), radial basis neural network (RBNN) and krigging (KRG) etc. are used approximate initial objective response data. There are several optimal solutions in multi-objective optimization (MOO). Many optimal solutions within the design space are required to find the Pareto-optimal front (PoF). A hybrid MOO consists of elitist non-dominated sorting of genetic algorithm (NSGA-II) (Goel et al., 2007).
Multiple surrogates modeling is implemented in many problems (Husain & Kim, 2010;Samad et al., 2008). The predicted error sum of squares (PRESS) through the kfold cross validation strategy finds the RMS error and filters out the inaccurate surrogates (Viana, Haftka, & Steffen, 2009). The error is calculated as: where, p is the number of points andẽ is the PRESS vector. For a given design space, a surrogate with the lowest PRESS RMS is defined as the best PRESS surrogate (BPS). The flow behavior through a turbine annulus at the tip is different from the hub. Varying the thickness from hub-to-tip varies the RB pitch from hub-to-tip, and it influences the flow pattern. The tip region is more liable to separation. The blade shape optimization is necessary to compromise between the reduction of frictional loss and the enhancement of shaft power (Takasaki et al., 2014).
This article focuses CFD-based optimization through surrogate-assisted multi-objective optimization approach to optimize the thickness of the RB profile from hub-totip of a bidirectional impulse turbine. The simulations are done with an in-house optimization code. Design parameters are modified, solved for flow, PoF are prepared and finally flow analysis is done to compare the original and optimal designs. Detailed approach along with the reason behind the performance improvement is presented in this article.

Numerical model and problem description
The reference geometry and profile parameters for bidirectional impulse turbine geometry are taken from an article of Maeda et al. (1999). The detailed specifications of the turbine and the blade profile are shown in Table 1   and Figure 2(a). The RB has circular profile on the pressure side (PS) and elliptic profile at the suction side (SS), and the maximum blade thickness located at 50% of the chord length. An unstructured grid with prism layers near the wall is constructed in the flow domain using ANSYS ICEM-CFD. Figure 3 shows the computational domain, mesh and design variables. The computational domain was extended from inlet of GV by 8.5 times of the chord length of GV to achieve the fully developed flow (Thakker & Dhanasekaran, 2004). In this work, Reynolds-averaged Navier-Stokes (RANS) equations are solved. The standard k-turbulence model with scalable wall function is used for the turbulence closure problem. The wall functions activate the local usage of the log law in regions where the y + is sufficiently small, in conjunction with the standard wall function approach in coarser y + regions. The inlet and the outlet boundary conditions are defined in Table 2. The performance of the impulse turbine under steady conditions is computed by the following non-dimensional parameters (Maeda et al., 1999): Input coefficient, Torque coefficient, For a higher pressure drop ( P), the performance of the turbine is poor because of the high separation on the SS, and the operating range of the turbine is reduced (Badhurshah & Samad, 2015). Hence, the minimization of was chosen as an objective function (f 1 ). The primary objective of any turbine is to produce more shaft power for a given flow coefficient or wave condition. Therefore, the shaft power was another objective function (f 2 ) to be maximized. As the flow is bidirectional, the airflow over the turbine is approximately sinusoidal based on regular wave theory. The wave height reaches its crest and changes its direction. Hence, the velocity of the air at the crest and the trough of wave is zero (Figure 2(b)), and the turbine does not get any power. In this case, the turbine dissipates energy. The minimization of pressure drop will increases the operating range which will increases the power extracting capability (Badhurshah & Samad, 2015). Hence, the objectives f 1 and f 2 are justified.
The blade-profile modification from the hub-to-tip changes the flow area, and the thickness of the blade at 50% chord length is changed by modifying the major radius of the elliptic profile (Figure 3(c)). Hub thickness (V 1 ) and tip thickness (V 2 ) are the design variables. Although there are many design parameters, the authors considered only thickness and checked if it can improve the performance of the turbine. The effect of each design variable on the objective functions is evaluated by simulating some preliminary designs. Using these results, a design space was created (Table 3) and the sample points were selected by the Latin hypercube sampling technique.

Optimization procedure
The CFD analysis requires long time to converge for each simulation, and evaluating many designs takes huge amount of time and computer resources. The systematic optimization methodology helps reducing the number of simulations and finds the optimal design among several alternatives. Hence, multiple-surrogate technique was used to create the initial population for MOO. RSA, KRG, RBNN and WAS were the surrogates constructed for this study. For four surrogates and two objective  functions (f 1 and f 2 ), eight surrogate models were produced.

Response surface approximation (RSA)
The response surface approximation methodology (Myers & Montgomery, 2002) is used to evaluate the responses obtained from the CFD model by fitting a simple polynomial to it. The response surface model can be expressed as: where, y(ξ ) is a unknown function depends on the controllable factors ξ 1 , ξ 2 , . . . .., ξ k . The symbol ε represents the other sources of variability that are inherent in the system. In this article, a second-order polynomial function is used to approximate the responses and it is given as: Where, β s are the regression coefficient and X s are the variable vectors.

Kriging model (KRG)
The KRG model estimates the unknown function as a combination of a constant global model g(x) and a smallscale systematic departure F(x). The F(x) creates localized deviation at some unknown point x using Gaussian correlation with variance (σ 2 ) (Simpson, Mistree, Korte, & Mauery, 1998).
The covariance matrix of Z(x) is given as, where, R(x i , x j ) is the spatial correlation function between any two sample points x i and x j and R is the correlation matrix.

Radial basis neural network (RBNN)
The RBNN uses a radial basis function as the activation function. It consists of a hidden input layer and an output layer. In this study, the Gaussian function (ψ) is used as an activation function and it is given as: where, σ is the Gaussian width, y is the design points, and y c is the known sample point. The approximation built by    radial basis function is given as: where, ω i is the weight coefficient and ω 0 is the bias term.

Weighted Average Surrogate (WAS)
WAS estimates a weighted function (y WAS ), which is a weighted average of multiple surrogates.
y WAS = ω KRG y KRG + ω RSA y RSA + ω RBNN y RBNN (13) where, ω KRG , ω RSA and ω RBNN , are the weights of the corresponding surrogates, and the sum of the weights is equal to one. The weights of the surrogates are estimated from the PRESS (Viana, 2011) and are calculated as: where, E i is the PRESS calculated error and N is the number of surrogates. Ref C

Best PRESS surrogate (BPS)
The BPS strategy finds the lowest PRESS surrogate. The BPS can be different for different problems, because quality of fit relies on the nature of data. For each surrogate, the PRESS for f 1 and f 2 were calculated.

MOO and clustering
MOO finds relationship among conflicting objectives and finds a set of optimal solutions. Every solution set has dominated and non-dominated solutions. It is expressed as: Where f i (x) is a vector of M objective functions andx is the vector of N design variables. Vector x i dominates over vector x j when x i is strictly better than x j for at least one objective and it is no worse than x j for all objectives. If this condition is not valid then x j is not dominated by x i . Between the two solutions, x i is the non-dominated solution if x i dominates x j . Global Pareto optimal set is obtained by using hybrid NSGA-II. Multiple surrogate models supplied the initial population and approximate Pareto optimal solutions are obtained (Deb, 2009). After that, k-means clustering algorithm is used to split the Pareto optimal solutions into five clusters. The detailed MOO procedure is shown in Figure 4.

Results and discussion
Initially, the turbine was designed in CAD software, the flow domain was extracted for a single blade passage and periodic conditions were used. The mesh was generated several times and checked if the results matches well with the experimental results of Maeda et al. (1999), and gives stable results. Based on the grid dependency test, the optimal number of elements was 1.28 million, as the efficiency of the turbine remained same between 1.2 million and 2.3 million for three different flow coefficients ( Figure 5).
The non-dimensional parameters such as input coefficient (C A ), torque coefficient (C T ) and efficiency (η) were validated with experimental results (Figure 6). The rotational speed of the RB was constant (600 rpm), and the inlet axial velocity was varied to get different flow coefficients. When ϕ < 1.2, there is a small degree of mismatch between the numerical and the experimental efficiencies because of the slight pressure difference occurs due to the surface finish issue. In the reference model for higher flow coefficients, the efficiency decreases and torque increases because inlet the axial velocity is high, which increases the torque. At the same time, the pressure drop and the energy loss are also high.
Nine design points for the two variables (V 1 and V 2 ) were evaluated by the RANS solver at ϕ = 1.25 and PoF was obtained (Figure 7). For each objective function, four surrogates were constructed and a total of eight surrogates generated the PoF via NSGA-II. The number of generations of NSGA-II was set to 250 (stopping criteria). The RSA model has the least PRESS for both the objectives (Table 4) and the RSA-RSA model was used as the BPS-BPS for further discussions. The KRG model has the highest PRESS and is ignored for further analysis in this work.
The results of the cluster points from each surrogate model are presented in Table 5. The cluster points obtained through k-means clustering algorithm for BPS-BPS model and the corresponding CFD results are shown in Table 6. The points A, B, C, D and E implied the cluster points or the designs obtained after clustering. The point A shows the point at extreme end point while the point E shows the opposite end in the PoF. To get a higher f 2 , the values of both the variables should be increased because increase in f 2 reduces the flow area and accelerates the flow as shown in Figure 8. At the same time, the thickness increases the frictional losses and pressure drop (f 1 ), and reduces the efficiency (f 0 ) and operating range.  The reference blade has a constant thickness from the hub-to-tip. The optimal design shows that a lower V 1 and a higher V 2 results in a better performance (Table 6). A higher V 2 reduces the area near the blade tip, and the flow remains attached to the blade surface as shown in Figure 8. The case A produces 6.85% higher shaft power and 4.23% lower pressure drop as compared to the Ref case. The case E increased the shaft power by 13.79% (Table 6). The efficiency (Eq 5) of the turbine is increased by 10.38% for C. Figure 9 shows the pressure distribution on the SS of the blade. The pressure on the tip and at the 50% chord length is higher for the Ref for ϕ = 2.75. The exit pressure of the Ref is higher than the optimized blade (case C) results in a high pressure region from the blade mid-chord to the TE, and the boundary layer separation occurs. In the optimized blades, static pressure is decreased at the 50% chord. This leads to the increase in flow and torque at the exit. The high-pressure regions surround more area of the PS for the case C implying the flow attachment ( Figure 10). Figure 11 shows the streamlines, and the separation lines are marked with dotted lines. For ϕ = 1.25, the flow is attached and the streamlines are uniform over the SS except the near the tip and the hub region. In higher ϕ,     (Figure 8). The separation region, which occupies about 40% of the span from the tip near the TE, is comparatively lower in the case of optimized blade C. In A and C, the separation is delayed, and the separated area is relatively less.
The flow near the hub is almost uniform in all the optimized cases. Uniform increase in thickness from hub-to-tip results in high frictional losses. Hence, the hub thickness is reduced and tip thickness is increased in A and C, the flow area near the tip is confined and the formation of passage vortex is reduced. Streamlines on the PS for ϕ = 2.75 is shown in Figure 11. The separation lines are visible on the PS of the case E as the thickness of the blade profile is uniformly increased. Hence, frictional losses on the blade E are also increased.
For ϕ = 2.75, streamlines at 90% of blade span is shown in Figure 8. In Ref, the passage vortex appears at 50% chord and enlarges itself as it moves toward the TE. In C, the passage vortex is reduced and shifted resulting in the higher torque in C.
A turbo line drawn from the hub to the tip on the SS midchord is shown in Figure 12. Figure 13 shows the static pressure distribution along the turbo line. At the mid-chord, the designs C and E show the minimum variation in pressure distribution along the span. There is a strong span-wise pressure-gradient in C and E, this drives the boundary layer fluid and secondary flows near the end wall toward the mid span direction near the SS. At the higher flow coefficients, a lower pressure region is spotted near the hub ( ∼ 10% span), where the velocity reaches the maximum and then decrease because of the spanwise adverse pressure gradient. The velocity distribution on the SS mid-chord is shown in the Figure 14.
The SS velocity at the mid-chord increases for C and E as compared to that of the reference case. This velocity in the end wall region is lower as compared to the midspan region because of the higher pressure and strong secondary flows near the end wall. Figure 15 shows the performance characteristics of the reference and the optimized blades for a wide range of flow rates. The pressure drop (f 1 ) is lower in C and A than in Ref, and it is the highest for E (Figure 15(a)). Figure 15(b) shows that the optimized cases produce higher shaft power (f 2 ). Figure 16 shows that the efficiency (f 0 ) of the turbine is increased for a wide range of flow coefficient for C.

Conclusions
In this article, a bidirectional impulse turbine is optimized for pressure drop minimization and shaft power maximization using CFD analysis and multiplesurrogate-based multi-objective optimization algorithm. The conclusions are: • Among the surrogates, RSA gave the least PRESS and was the BPS for both the objectives. It produced the least RMS error after evaluation by NSGA-II and CFD analysis. • A maximum pressure-drop of 3.6% and increase in shaft-power of 12.33% were obtained at the extreme end designs in the Pareto optimal front. The shaft-power enhancement was due to the low pressure at the blade exit and the flow attachment near the leading edge of the SS tip. This results in a steep velocity gradient on the pressure and SS. • The decrease in hub thickness and increase in tip thickness improves the turbine performance in terms of efficiency by compromising between the frictional loss and shaft power. • Increase in the tip thickness confines the flow area near the tip and reduces the formation of passage vortex on the suction side. • The efficiency, which is the function of both objectives, is relatively increased by 10.38% for case C. The final values of V 1 and V 2 are 15.4 and 17.6 mm respectively, where the reference value is 16 mm for both variables.

Symbols
b blade height (m) C A input coefficient C T torque coefficient E rms root mean square error f 0 combined objectives (η) f 1 objective function 1( P) f 2 objective function 2(P sh ) lr RB chord length (m) n span normalized P sh shaft power (W) Q flow rate(m 3 /s) r R mean radius (m) To output torque (N-m) U R circumferential velocity at mean radius (m/s) V a mean axial velocity (m/s) w i weight for a surrogate z number of RBs P pressure drop between settling chamber and atmosphere (Pa) η turbine efficiency ω angular velocity of RB (rad/s)

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was supported by Dr Abdus Samad Acknowledges Korean Federation of Science and Technology Societies (KOFST) for the 'Brain Pool' foreign scientist invitation award: [Grant Number none].