Optimization Study on the Blade Profiles of A Horizontal Axis Tidal Turbine Based on BEM-CFD Model

In order to increase the performance of horizontal tidal turbines, a multi-objective optimization model was proposed in this study. Firstly, the prediction model for horizontal tidal turbines was built, which coupled the blade element momentum (BEM) theory and the CFD calculation. Secondly, a multi-objective optimization method coupled the response surface method (RSM) with the multi-objective genetic algorithm NSGA-II was applied to obtain the optimal blade profiles. The pitch angle and the chord length distribution were chosen as the design variables, while the mean power coefficient and the variance of power coefficient were chosen as the objective functions. With the mean power coefficient improved by 4.1% and the variance of power coefficient decreased by 46.7%, results showed that both objective functions could be improved.


Introduction
In recent decades, high demand and consumption of energy have caught a serious problem around the world (Butala and Novak, 1999). Under the pressure of economy development and environment protection, to develop renewable energy is a practical way for sustainable development (Wang et al., 2011). Oceans contain an enormous amount of renewable energy in different forms such as thermal, wave, and tidal current energy etc (Uihlein and Magagna, 2016). With such characteristics as vast, reliable, regular and predictable, tidal current energy is one of the most advanced forms which is expected to contribute significantly to future energy supply (Uihlein and Magagna, 2016;Ponta and Jacovkis, 2008).
In coastal area with reasonable tides, a tidal turbine is one of the commercial devices to utilize tidal current energy with low environmental impact. Two types of turbines are usually adopted: horizontal axis and vertical axis (Derakhshan et al., 2017). As for wind turbines, horizontal axis wind turbines (HAWT) are more efficient than vertical axis wind turbines (VAWT) (Ferdoues et al., 2017). Since the design of tidal turbines was based on existing wind turbines, horizontal axis tidal turbines (HATT) have wider applications than vertical axis tidal turbines (VATT) for the high efficiency. However, certain differences between HAWT and HATT need to be investigated. The most obvious difference between HAWT and HATT is that water is approximately 800 times denser than air, and moving water holds considerably more energy than that of air. Consequently, with shorter and thicker blades, etc., tidal turbines should be modified to sustain heavy loadings and energy (Ouro et al., 2017).
During recent years, HAWTs have attracted much attention, and much development has taken place including model testing of prototype's performance and installation (Goundar and Ahmed, 2013). Tests were carried out by Bahaj et al. (2007) in cavitation tunnel and in towing tank of an 800 mm diameter HATT model to study the turbine's power and thrust coefficients in ranges of tip speed ratio and pitch angle. Results showed the operation of a single turbine under various conditions and simultaneously provided detailed experimental data for the validation of numerical models later. Huang et al. (Huang and Kanemoto, 2015;Huang et al., 2016) designed a bi-directional counter-rotational type HATT using fully symmetrical hydrofoils. Experimental tests and computational fluid dynamic (CFD) simulations were both implemented. Results indicated that the counter-rotating rotors performed much more excellent in power coefficient than a single rotor did. The fully symmetrical hydrofoils, blade pitch angle distribution and axial distance between the front and rear blades were optimized to enhance the performance of the counter-rotating turbine (Wei et al., 2015). A high performance and cavitation free hydrofoil named HF-Sx was designed by Goundar et al. (Goundar et al., 2012;Goundar and Ahmed, 2013) and a 10 m diameter, 3-bladed HATT was designed using these new hydrofoils to remove from cavitation under expected operating conditions. Besides, a blade element momentum (BEM) model was developed by Batten et al. (2008) to predict the performance of tidal turbines. By comparing the predicted results with experimental data obtained from tests of an 800 mm diameter rotor mentioned in Bahaj et al. (2007), the theoretical predictions were demonstrated to be in good agreement with the experiments (Batten et al., 2008). Lee et al. (2015) investigated the effect of distance between the dual rotors on the performance of a counter-rotating current turbine by using CFD simulations and experimental methods. It was found that the dual rotor produced more power than the single one and the blade gap distance can be used as a parameter for counter rotating turbine design. However, due to the limited computational resources, only one blade was performed numerical simulations with periodic boundary conditions. A coupled BEM-CFD model for predicting performance of tidal turbines was developed by Malki et al. (2013). This method could take the influence of upstream hydrodynamics and free surface on rotor performance into calculation. Results from this model fitted well against the experimental data and the classical BEM theory, which made it a useful method to consider more complex conditions such as tidal turbine array. Belloni et al. (2017) compared the performances of bare, ducted, and open-center turbines for axial flow. A model coupled the Reynolds-averaged Navier-Stokes (RANS) and the BEM theory was applied. Results showed that moderate increases in power density can be increased through the ducted and open-center devices, however the power generated by these two kinds of turbines was lower than the bare one.
The hydrodynamic performance of an HATT mainly depended on the blade geometry and the operating conditions. The blade geometry parameters included twist angle, chord length and thickness distributions while the operating conditions mainly referred to the tip speed ratio (TSR) which was the ratio of blade tip speed to inflow velocity. Many parameters were used to describe the performance of turbines, such as power coefficient, startup time, cost of energy, thrust coefficient etc. The purpose of optimization was to find the optimal solution for all objectives which was multiobjective optimization. Owing to the fact that there are so many objectives to optimize, designers usually just focused on two or three objectives. Multi-objective optimization was widely used in wind turbine optimization and gradually used in tidal turbines in recent years. Benini and Toffolo (2002) described a multi-objective optimization of a stallregulated HAWT. They found an optimal solution to achieve the best tradeoff performance between annual energy production per square meter of wind park and energy consumption. Genetic Algorithms (GA) was chosen to find the final optimal solution. Hwang et al. (2013) and Abdulrahman and Wood (2017) adopted GA to optimize wind turbines in different objectives. Guo et al. (2013) applied the second version of non-dominated sorting genetic algorithm (NSGA-II) to optimize the pitch angle distribution for a higher power coefficient and a lower thrust coefficient. Besides, the response surface technology was used to establish the relationships between the variables and the objective functions to save the computational cost. Four different multi-objective optimization algorithms, NSGA-II, multi objective particle swarm optimization algorithm (MOPSO), multi objective cuckoo search algorithm (MOCS) and multi objective flower pollination algorithm (MOFPA), and their performances were evaluated in combination with BEM theory by Tahani et al. (2015). Results showed that the Pareto fronts obtained from the MOFPA and NSGA-II had better quality than the others.
Few studies focused on automatic multi-objective optimization combined with the prediction method of BEM-CFD for tidal turbines. In this study, an automatic multi-objective optimization process was carried out trading off with the variance to obtain the maximum average power in a range of TSRs. The tidal turbine designed by Bahaj et al. (2007) was chosen as the prototype and the BEM-CFD model was built on MATLAB and ANSYS CFX to predict the performance of the turbine. Bezier curve technology and response surface technology were introduced to the optimization process to reduce the design variables. The whole optimization process was established on ISIGHT, an excellent multi-disciplinary optimization design platform.

Governing equations
The CFD simulation is to get solutions of the Navier-Stokes equations which represent the conservation of mass and momentum. For steady state fluid, they can be expressed as: where is the density of fluid, is the velocity vector, u i is the i-th component of the velocity vector, and are the laminar and turbulent dynamic viscosities respectively, and S i includes additional sources such as the sources generated by the rotating turbines.
There are many turbulence models having been developed (Malki et al., 2013;Edmunds et al., 2017). Being simple, stable and having low computational cost, tur-k − ε bulence model is most widely used for simulating the mean flow (Launder and Spalding, 1972). However, with a single length scale, it is insufficient to calculate the turbulence viscosity. As a result, the diffusion occurring at time scales is not taken into consideration. Focusing on this problem, the RNG model has a better performance by using statistical analysis towards different scales of motion. However, the turbulence is considered as isotropic in these models, which is an unreasonable approximation in rotational flows. Higher order turbulence models can produce better prediction like the Reynolds stress model and large eddy simulation although they are computational costly. k − ε ε As the main object of this study was a turbine rotor rather than the fluid structure in high turbulence region, model was chosen as a tradeoff between computational cost and accuracy. Two equations were to be solved in this model, k represents the turbulent kinetic energy and represents the dissipation of this energy: Turbulent viscosity is calculated: In Eqs.
(3)-(5), , , , and are constants, and G represents the turbulence generation rate. The effective viscosity used in the diffusion equation in the momentum equation, the k equation and equation is the sum of the laminar and turbulent viscosities.

BEM-CFD method
The BEM-CFD method has been used in tidal turbines for performance prediction and the results show good agreement with experimental data (Shives and Crawford, 2017;Guo et al., 2015;Belloni et al., 2017). It is a method coupled the BEM theory with CFD model through source items which indicates the feature of blades.
BEM theory is a hybrid method to analyze the aerodynamic performance of rotors. It combines momentum theory and blade element theory. In this method, the rotor is simplified by actuator disk theory. The influence of a rotor is time-averaged so that the force exerted on the rotor is the same on a given radius of the plane. With knowing lift and drag coefficients (C L and C D ) database of each blade element, the hydrodynamic performance at the given radius can be determined. Fig. 1a shows the discretization of a rotor, usually discretized into 10-20 elements. For each blade element, the geometric parameters include the chord length c, the thickness t and radial width dr, as shown in Fig. 1b. Fig. 1c shows the force analysis on a blade element.
The lift F L , and the drag F D , forces are defined as follows: α where, C L and C D are the lift and drag coefficients as the function of the angle of attack , and the relative velocity W can be expressed as: where a and are the axial and tangential induction factor respectively, V 1 is the inflow velocity and is the angular velocity of the rotor. The axial and tangential forces on the blade can be determined by projecting the lift and drag forces onto the axial and tangential directions as shown in Fig. 1c: ϕ where is the flow inclination angle defined by: Substituting Eqs. (6) and (7) into Eqs. (9) and (10) yields: The two components were the source terms S i in Eq.
(2) converted into forces per volume, as shown in Eq. (14). The process of BEM-CFD method was shown in Fig. 2. The CFD model was settled with appropriate initial and boundary conditions. Source terms obtained from BEM theory were introduced to CFD model which were dependent on the blade characteristics and the TSR. New source terms in the next iteration can be obtained by substituting the new velocity and pressure field into the BEM model. Thus, the rotor's geometry can be simplified into a disc with exact thickness. As a result of that, high quality meshes can be created.

Definitions
The tip speed ratio (TSR) is the ratio of the blade tip speed to the free stream velocity: The power coefficient C P is an integral of Eq. (10) from the hub to tip, which is defined as: where A represents the swept area. In this study, C P was chosen as the objective to assess the performance of a rotor. A high average value and a low variance of power coefficient over a range of TSRs were the aims of this optimization. It stood for a good performance in a range of operating conditions.

Computational domain and boundary conditions
In this study, the turbine designed by Bahaj et al. (2007) was chosen as the object. Experiments in a towing tank and a cavitation tunnel had been carried out to investigate the operating performance. This rotor is a three-bladed rotor with the diameter (D) of 0.8 m. The simulation domain is a rectangle similar to the experimental tank. By extending 2D upstream and 8D downstream, the width and depth of the tank cross section is 3D and 1.75D respectively, where 1.5D from the bottom is the interface of water and air. Fig. 3a illustrates the simulation domain. The thickness of rotor e was set to 0.08 m according to Guo et al. (2015). The nacelle and the support structure were also modeled with the diameter of 0.1 m and 0.08 m, respectively. The hub was symmetrical at the front and back cones.
The commercial code ANSYS ICEM was used to generate a structured hybrid mesh that consisted of hexahedrons and quadrilateral. There were 5.7×10 5 nodes, 5.5×10 5 hexahedral meshes and 4.3×10 4 tetrahedral meshes as shown in Fig. 3b.
The uniform velocity at inlet is 1.73 m/s with 5% turbulence intensity and the reference pressure at outlet was set to zero which means that the flow was discharged into the atmosphere. The wall at the top was specified as an opening condition and became free surface with the relative pressure of 0 Pa. The rest of the three walls, the hub, the nacelle and the support structure, were modeled as walls without slip. The model was solved via ANSYS CFX 17.0, a commercial CFD finite volume solver. The simulation was solved with a high resolution advection scheme until it reached the max iterations of 50 or the RMS residual target of 1.0×10 -5 .

Variables and objectives
According to the report by Batten et al. (2006), the hydrodynamic design parameters of a rotor as a whole basically consisted of the rotor diameter, rotational speed and the blade geometry. The blade design included the hydrofoils, the chord length distribution and the pitch angle distribution or 'twist' of the blade. Since the optimization design based on an existing turbine with certain hydrofoils, the The model turbine blade was discretized into 17 elements, meaning that the pitch angle distribution and chord length distribution were determined by 17 points respectively. The geometries of 9 elements among them are shown in Table 1. According to Bezier curve, curves were fitted through control points. It was suitable for the geometric description (Kharal and Saleem, 2012). In this study, the cubic Bezier curve with 4 control points was employed to get the desirable distributions as shown in Fig. 4. The angles and the chord lengths c 1 -c 4 determined the pitch angle and the chord length c in the radial direction, which were the eight variables adopted in the optimization process.
At a given tidal farm area, it was of primary interest to obtain the maximum power output under a stable operating condition of a tidal turbine. The power output was determined by the kinetic energy available in the moving water, the coefficient of energy conversion by the rotor and other losses in the whole system. The coefficient of energy conversion had been defined in Sec.2.3 and the maximum of which was 16/27 according to Betz limitation. In this study, the power coefficients at TSR=4, 6 and 8 were taken into consideration while TSR = 6 was the optimum operating condition. Two objectives are shown as follows: where C P, 4 , C P, 6 and C P, 8 represent the power coefficient at TSR = 4, 6, 8, respectively; C PV is the variance of these three values. The aim of this paper is to obtain the optimal pitch angle and chord length distribution to achieve the maximum of C PM trading off with C PV .

Multi objective genetic algorithm
Genetic algorithm is a method for solving optimization problems which is inspired by the process of natural selection. It is one of the most widely used and researched optimization methods. It differs from typical optimization method in that it generates a population of solutions at each iteration rather than a single point. At each step, the genetic algorithm selects individuals at random from the current population named parents and then generates the next generation by crossover and mutation operators (Tahani et al., 2015;Kharal and Saleem, 2012;Selig and Coverstone-Carroll, 1996).
As for a successive computation, the populations will lead to optimal solutions which are called Pareto solutions. For a typical process of genetic algorithm, a population of potential solutions for optimization problem is created and initialized randomly at first. Each individual is assessed according to its fitness or measure of goodness. The higher an individual's fitness is, the larger the probability that the individual is to live and become the best solution in this generation which are called the elitism. They are going to propagate to the next generation. The elitism ensures the survival of the most suitable individual from one generation to the next without altering its string representation. Then the crossover and mutation are employed for the purpose of exploitation and exploration respectively. The crossover operator combines two parents selected randomly and creates offspring. The mutation operator apply random changes to parents in the values of genes at one or more positions and then form the next generation. Mutation also provides a mechanism to avoid getting trapped in a local extreme value (Tahani et al., 2015;Kharal and Saleem, 2012;Selig and Coverstone-Carroll, 1996).
Among the multi objective genetic algorithm, NSGA-II is a popular method to optimize in both wind and tidal turbines as mentioned in Section 1. Fig. 5 shows the structure of NSGA-II.
There are two concepts used in NSGA-II to sort solutions, non-dominated sorting and crowding distance sorting. Firstly, the children population (Q t ) and the parent population (P t ) are combined together to form a new population (R t ). The size of R t is 2 N, summation of P t and Q t . Secondly, the population R t is sorted according to non-domination. The best solutions are classified into the best nondominated set F 1 . These are the most important solutions in R t . If the size of F 1 is smaller than that of N, all solutions in set F 1 are chosen for the new population P t+1 . Then set F 2 are chosen to form P t+1 , followed by F 3 until the population P t+1 is filled with enough solutions at set F t . If the population of P t+1 is larger than N, set F t is going to be sorted by crowded -comparison operator to remove the excess solutions. At last the population of P t+1 of size N is used to create an offspring population Q t+1 (Deb et al., 2002).

Optimization methodology
The multi objective optimization coupled the Bezier curve and the response surface method (RSM) with the multi objective genetic algorithm (NSGA-II). With these techniques introduced to the platform ISIGHT, the optimization model can be automatic and low in computational cost.
The flow chart of the optimization model is shown in Fig. 6. The objectives and design variables has been defined in Section 3.1. Then sample points were extracted in the range of variables according to Latin hypercube test design. For each sample point, new pitch angle and chord length distribution would be fitted by Bezier curve. Besides, the objective functions with new blade profiles were calculated through BEM-CFD. After that, the fitting relationships between the eight variables and the two objective functions were established by applying the response surface technique. Four types of error analysis were applied to examine the accuracy of the approximation. On the basis of the relationships between the objective functions and variables presented by quadratic polynomial regression equations, the multi-objective optimization of the blade profiles was executed using the algorithm NSGA-II, and a final optimal solution would be chosen from the Pareto frontier solution obtained from the optimization.

Results and discussion
α Blade element theory gave the forces on each blade element as functions of blade geometry and lift and drag coefficients which were related to the angle of attack . The two-dimensional C L and C D used in the BEM-CFD model in this study were calculated by CFD simulation. Since the rotor operated at different TSRs during the numerical program, the values of C L and C D also had to be calculated with different Reynolds numbers. Therefore, each foil profile was simulated at Reynolds numbers of 1×10 5 , 2×10 5 , 3×10 5 , and 4×10 5 to provide an extensive data base for interpolation at a later stage. Fig. 7 illustrates sample C L and C D values for 5 different NACA 63-8×× foils generated at a Reynolds number of 3×10 5 .
Based on the database of C L and C D , the power coefficient from BEM-CFD was compared with the test results carried out in a towing tank as shown in Fig. 8. It showed that the simulation could provide reasonable predictions. Specifically, the simulation agreed well with the experi-  (Deb et al., 2002).  ZHANG Da-hai et al. China Ocean Eng., 2019, Vol. 33, No. 4, P. 436-445 441 mental data near the optimum operating conditions of TSR = 6-8. However, the BEM-CFD model seemed to over pre-dict the power coefficient in the high TSR region and to under predict in the low TSR region. It might be attributed to the accuracy of tip loss factor used in BEM theory, neglect of mechanical losses and wake model in the CFD simulations etc.
Since the eight variables were determined in Section 3.1, the ranges of these variables are shown in Table 2 where the angle and chord length were adjusted between the upper and the lower bounds. 90 experimental samples were generated according to the range with the help of Latin hypercube experimental design method as shown in Fig. 9, with the pitch angle and chord length distribution of these samples being fitted by Bezier curve. The objective functions of the 90 experimental samples calculated through BEM-CFD method are shown in Fig. 10. Then the response surface method was employed to establish the fitting relationships between the eight optimization variables and objective functions in the form of quadratic polynomial regression equations. Table 3 shows four different types of error analysis of this approximation. It indicates that all four types of error analysis of the approximation achieve the acceptance level, which means that the fitting relationships of C PM and C PV are reliable.
The aim of this optimization was to obtain the optimal pitch angle and chord length distribution which achieved a tradeoff between the C PM and C PV , where C PM represented the mean performance near the optimum operating condition of TSR = 6 and C PV represented the stability of operation. Since the rotor could not remain operating at the optimum condition, a high C PM ensured that the rotor operated efficiently near the optimum condition while a low C PV kept it in low vibration when the operating condition changed.
Then the multi-objective genetic algorithm of NSGA-IIwas carried out to obtain the optimal solution. In this study, the population size and the number of generations were both set to 100, and the remaining parameters were default values. Fig. 11 plots the optimization results where each solution presented a possible blade profile and the final optimal solution was determined by the tradeoff between C PM and  C PV . Three optimal solutions were chosen from the Pareto frontier solution and were named as optimal blade A, B and C respectively. Performance comparison between C PM and C PV is shown in Table 4, showing that optimal blade A has a better hydrodynamic performance in C PM while optimal blade C had less operating fluctuation near the optimum condition. The pitch angle and chord length distribution of these three optimal blades are shown in Fig. 12. It can be seen that the pitch angle at the hub of all these three optim-al blades increased about 3 degrees compared with the original blade. It also shows that optimal blade A has a larger twist angle range than the original and other two optimal blades. Otherwise a smaller slope of pitch angle distribution ensured the rotor operated in a relatively stable condition shown by optimal blade C. As for the chord length distribution, three optimal blades almost kept the same and were longer along the whole radial direction compared with the original blade. The power coefficients of these four blades were verified by BEM-CFD method at the range of TSR=3-10. It can be seen from Fig. 13 that the results agree well with the optimization. Optimal blade A with larger twist angle range has a better hydrodynamic performance considering the mean power coefficients. Specially, in the high TSR range, the power coefficients were larger than the others. And optimal blade C has a stable operating performance especially considering the TSR of 5, 6 and 7. Optimal blade B performed well both in mean power coefficients and operating stability.
In summary, the hydrodynamic performance and the operating stability of the rotor were enhanced at the same time by altering the pitch angle and chord length distribution.
As shown in Figs. 14-17, the velocity fields of original and optimal blades are different. It can be seen from the figures that the hub part and supporting structure of the turbine have a great influence on the downstream and near the surface flow field of the turbine, which cannot be considered by the BEM method.

Conclusions
In this study, the blade profiles including the pitch angle and chord length distribution of a horizontal axis tidal turbine was optimized automatically, and the conclusions have been drawn as follows.    (1) A BEM-CFD model was developed to predict the performance of the horizontal tidal turbine. This model took the effect of nacelle, upright and free-stream into account and it improved computational efficiency compared with 3D CFD simulation. It could provide reasonable predictions in the operating regions especially near the optimum operating condition.
(2) An automatic multi-objective optimization method has been presented which coupled the Bezier curve parameterization, response surface method and multi-objective genetic algorithm based on the optimization platform of ISIGHT. Comparisons between the optimal and original blade profiles showed that a larger twist angle range can improve the mean energy coefficient. However, it will decrease the operating stability.
(3) In consideration of the tradeoff between the mean power coefficients and the operating stability, optimal blade B with a reasonable higher twist angle range and chord length performs better, whose C PM and C PV achieve to 0.4050 and 0.0033, respectively.
(4) However, this study did not take cavitation and mechanical losses into account which was an important problem for tidal turbine. Besides, the accuracy of BEM-CFD method also should be improved.