Optimal Design and Aerodynamic Performance Prediction of a Horizontal Axis Small-Scale Wind Turbine

Yin, Rui; Xie, Jian-Bin; Yao, Ji

doi:https://doi.org/10.1155/2022/3947164

Mathematical Problems in Engineering

On this page

Abstract Introduction Conclusions Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2022 | Article ID 3947164 | https://doi.org/10.1155/2022/3947164

Optimal Design and Aerodynamic Performance Prediction of a Horizontal Axis Small-Scale Wind Turbine

Rui Yin,¹Jian-Bin Xie,^1,2,3and Ji Yao¹

Academic Editor: Rodolfo Bontempo

Received21 Jul 2021

Revised23 Feb 2022

Accepted13 Apr 2022

Published09 May 2022

Abstract

The improved aerodynamic design of a horizontal axis small-scale wind turbine blade is crucial to increasing the efficiency and annual energy production of the turbine. One of the vital stages in aerodynamic design is the selection of the airfoil. Using the existing airfoils for a blade design which results in higher turbine characteristics is tedious. Consequently, this paper provides an optimal design strategy for a horizontal axis small-scale wind turbine blade through the multiobjective optimization of the airfoil using the Nondominated Sorting Genetic Algorithm II (NSGA-II). The latter outperforms the other commonly used genetic algorithms (GAs), as well as the Computational Fluid Dynamics (CFD) investigation of the different airfoil types and the wind turbine rotors on the steady or unsteady state aerodynamic performance. An NACA4412 airfoil with higher aerodynamic efficiency is considered as a baseline for the optimization in order to increase the lift coefficient and lift to drag ratio while avoiding excessive variations in the maximum relative thickness and area. The optimized airfoil (NACA4412-OPT) is used as the cross-sectional profile in the design procedure for a novel 1.15 m diameter three-bladed wind turbine rotor at a wind speed of 11.5 m/s, tip speed ratio of 4.65, and pitch angle of 0.2° by the Wilson design method. The two-dimensional analysis demonstrates that the optimized airfoil outperforms the other airfoils yielding the highest lift coefficient and lift to drag ratio, as well as a larger pitch range. The three-dimensional analysis shows that the time-averaged power coefficient value (0.33) of the new wind turbine is almost 26% higher and more stable than that of the original wind turbine while avoiding a high increase of the axial thrust.

1. Introduction

The small-scale wind turbines (SWTs) are decent to overcome some sustainable development issues of the large-scale wind turbines (LWTs) associated with taking up space, visual interference, and noise pollution [1]. More precisely, they can provide cost-effective electricity in some areas where the average wind speed sometimes does not meet the requirements of an LWT and the demand for electricity is relatively low [2]. Consequently, SWTs are widely used in recent years, while most of the studies focus on the wind turbine blade design.

For instance, Sugathapala et al. [3] propose a design and optimization procedure of simplified wind turbine rotors for small-scale applications to design optimum rotors suitable for wind characteristics at specific locations while enhancing the local manufacturing capacities. The section profile along the blade is uniform as the NACA4412 airfoil. The performance evaluation using the Blade-Element Momentum (BEM) theory demonstrates that the performance reductions of the simplified rotor designs are not significant, and therefore they will be effective in enhancing the value addition through the local manufacture.

Abdelsalam et al. [4] study the aerodynamic performance of two different types of horizontal axis SWT rotors by experiments. The first one is a classical rotor with nonlinear chord and twist distributions, while the second one is a new linearized rotor design using the BEM theory. The section profile along the blade is uniform as the RISØ-A-24 airfoil. The new blades have a reduced size volume by 26%, approach the power coefficient achieved by nonlinear blades, and have higher starting ability and extended operating range at lower wind speed compared with the classical model.

Lee et al. [5] study the aerodynamic performance of two different types of horizontal axis SWT blades by experiments, as well as numerical simulations based on CFD. The first one is a typical type designed using the BEM theory, while the second one is a nontwisted type having a constant chord length. The section profile is uniform along the blade, and the SD8000 airfoil is considered. Both the numerical simulations and measurements demonstrate that the maximum power coefficient of the BEM blade is increased by more than 50% compared with the baseline blade.

Pourrajabian et al. [6] designed a horizontal axis SWT blade while considering the starting time, output power, and blade allowable stress using the BEM theory, simple beam theory, and a GA. The section profile is uniform along the blade, and the SG6043 airfoil is considered. They deduce that a hollow blade can expedite the starting, and it is also more powerful compared with the solid ones.

Rahgozar et al. [7] designed a 1 m timber horizontal axis SWT blade while considering the starting time and output power for four possible combinations of linear/nonlinear distributions of the chord length and twist angle, using the BEM theory and a GA. The section profile along the blade is uniform as the SG6043 airfoil. They deduce that the use of linear distribution can improve the starting performance at a lesser compromise of output power.

Gupta et al. [8] study the design and aerodynamic performance of two different types of horizontal axis SWT rotors using the BEM theory and CFD, respectively. The first one is designed using the SG6050 airfoil and SG6043 airfoil, while the second one is designed using the E555 airfoil and E216 airfoil. The BEM theory and CFD studies show that the power coefficient of the more suitable rotor is 0.445 and 0.35687, respectively.

Hasan et al. [9] present a study between the BEM theory analysis and CFD analysis of a horizontal axis SWT blade. The cross section of the blade consists of the S823 airfoil (root), S833 airfoil (middle), and S822 airfoil (tip). The BEM theory and CFD studies demonstrate that the best power coefficient of performance at 6° angle of attack is 0.47 and 0.43, respectively.

Muhsen et al. [10] designed a horizontal axis SWT blade and enhanced the design using the BEM theory and the QBlade software in order to obtain a power coefficient higher than 40% at a low wind speed of 5 m/s. The cross section of the blade consists of two symmetric sections belonging to the S1210 airfoil and S1223 airfoil. The final design can obtain almost 650 W with a power coefficient of 0.445 at a wind speed of 5.5 m/s, reaching power of 1.18 kW and a power coefficient of 0.40 at a wind speed of 7 m/s.

Pholdee et al. [11] present an optimization method of new MEXICO blades for a horizontal axis wind turbine in order to maximize the power coefficient, while the design variables are twist angles on the blade radius and rotating axis positions on a chord length of the airfoils. This approach uses a GA based on the Kriging surrogate. The CFD analyses show that the Kriging model with linear regression leads to better results than the Kriging model with second-order polynomial regression. The optimum blade obtained in this study shows a better performance than the original blade at a low wind speed of 10 m/s.

Shen et al. [12] optimize the geometry of nonstraight SWT blades not only in terms of the distribution of the chord and twist angle but also with a three-dimensional (3D) stacking line using a microgenetic algorithm and a lifting surface method with a free wake model, in order to maximize the annual energy production and improve the starting performance. They deduce that the wind turbine blades with a properly designed 3-dimensional stacking line can increase the annual energy production and have a better starting behavior than the two-dimensional optimized blade geometries.

Yang et al. [13] optimize the chord and twist angle of a 5 MW NREL wind turbine blade using the NSGA-II and the BEM theory in order to maximize the annual energy production and minimize the blade mass. The optimization design method can provide a better blade while increasing the annual energy production by 2.48% and reducing the blade mass by 5.52%.

Maheri [14] uses an automated wind turbine blade modeler within an optimization process in order to minimize the mass, the maximum von Mises stress, and the blade tip deflection for an adaptive NREL 5 MW blade using NSGA-II and generate files required by ANSYS for high fidelity analysis. The automated wind turbine blade modeler can treat any parameter required to define the size, topology, structure, and material of a blade as a design variable and automatically change them within an optimization process. And therefore it conducts an integrated design optimization.

The reviewed studies investigating blades are presented in Table 1.

The performance evaluation method of wind turbines can be used to perform model experiments in the wind tunnel or simulations using the BEM theory and CFD. The former can be used to generate relatively real data. However, it has the disadvantages of being time-consuming and expensive [15]. On the contrary, the BEM theory and CFD can generate diverse and reliable results with a low cost. Refan and Hangan [16] experimentally and theoretically evaluate the aerodynamic performance of an upwind, three-bladed, small horizontal axis wind turbine rotor of 2.2 m diameter. A comparison between the theoretical and experimental results shows that the overall prediction of the BEM theory is within an acceptable range of accuracy. However, the BEM theory prediction for SWTs is not as accurate as that for LWTs. Plaza et al. [17] perform an aerodynamic analysis of the MEXICO wind turbine rotor while comparing between the results of the BEM theory, CFD, and measurements. The results demonstrate that BEM calculations outperform the CFD estimates at low wind velocities. However, they fail at higher velocities because of separated flow conditions. The blade tip loss and 3D effects are partly responsible for the calculation inaccuracies, especially for the BEM theory. In [5], the numerical simulations by CFD for the aerodynamic performance of two different types of horizontal axis SWT blades are coherent with the experimental data. Moussa [18] evaluates the evolution of the power coefficient of a horizontal axis 3-bladed SWT with respect to the tip speed ratio by experiments and CFD. Empirical and numerical models of the evolution of the power coefficient with respect to the tip speed ratio are, respectively, confronted with the experimental results showing a good agreement with an estimated maximum error of almost 5%. Thus, in this paper, the CFD method is used to evaluate the aerodynamics performance.

It can be deduced from the literature that although the Multiobjective Evolutionary Algorithms (MOEAs), such as GA and NSGA-II, are used for multiobjective optimization of wind turbine rotor, NSGA-II is not used in any previous optimization studies related to SWTs. Moreover, few studies tackle the multiobjective optimization of the airfoils shape. Furthermore, there is a household horizontal axis SWT which has been on the market for years. It has poor efficiency and has not been updated in China. In fact, it has hindered the development of SWTs, and therefore it is necessary to optimize the design of the original wind turbine. Thus, this paper provides an optimal design strategy for a horizontal axis SWT blade by the multiobjective optimization of the airfoil using NSGA-II and CFD investigation of the different airfoil types and wind turbine rotors on the steady or unsteady state aerodynamic performance. The optimization aims at increasing the lift coefficient and lift to drag ratio while avoiding excessive variations in the maximum relative thickness and area. A novel 1.15 m diameter three-bladed wind turbine rotor at a wind speed of 11.5 m/s, tip speed ratio of 4.65, and pitch angle of 0.2° is then designed using the optimized airfoil as the cross-sectional profile by the Wilson design method. Afterwards, the main aerodynamic characteristics of the optimized airfoil and the new wind turbine are compared with the other airfoils and the original wind turbine, respectively.

2. Multiobjective Optimization of the Airfoil

2.1. Parametric Expression

The traditional methods of airfoil parameterization are divided into polynomial fitting method [19], free-form deformation technique [20], and linear superposition method [21]. The Hicks–Henne shape functions are a widely used linear superposition method having a better effect on the local fine-tuning [22]. Therefore, the airfoil shape is formed by attaching 10 Hicks–Henne shape functions to the initial airfoil with fixed leading edge and trailing edge locations. The shape of airfoils can be expressed as follows [21]:where and are coordinates on the top and bottom surfaces of the new airfoil, respectively; and are the coordinates on the top and bottom surfaces of the initial airfoil, respectively; n is the number of design variables; are the Hicks–Henne shape functions; are the coefficients of shape functions, which are used as design variables to modify the shape of the airfoil. The limits of the 10 coefficients are [23] δ₁ ∈ [–0.006, 0.006], δ₂ ∈ [–0.01, 0.006], δ₃ ∈ [–0.006, 0.008], δ₄ ∈ [–0.005, 0.005], δ₅ ∈ [–0.005, 0.01], δ₆ ∈ [–0.008, 0.006], δ₇ ∈ [–0.01, 0.01], δ₈ ∈ [–0.005, 0.01], δ₉ ∈ [–0.005, 0.005], and δ₁₀ ∈ [–0.005, 0.015].

The Hicks–Henne shape functions arewhere x_i can work in accordance with mathematical model of particle swarm optimization.

2.2. Numerical Calculation

The calculated domain boundary is a circle with a radius about 20 times the chord length (c) of the airfoil. The trailing edge of the airfoil coincides with the center. A full structured grid is used. The grids for the calculation are shown in Figure 1, and the local grid around the airfoil is the finest. The boundary condition of the domain boundary is velocity-inlet and pressure-outlet, and the boundary condition of the airfoil is set to wall. The scheme of the pressure-velocity coupling is the COUPLED method, and the pressure term is calculated in a standard format. The second-order upwind scheme is used to calculate the momentum, turbulent kinetic energy, and specific dissipation rate. The turbulence model is the k-ω SST model [24], in which the convergence criterion is k and ω below 1 × 10⁻⁶.

(a)

(b)

If y + on the airfoil surface is less than 1, the number of grids has little effect on the calculation results, which are shown in Table 2, where the Reynolds number is 6.4 × 10⁵ and the angle of attack α is 4°. The calculation results of Grid 2 and Grid 3 are very close, and Grid 2 has fewer iteration steps. There is a certain gap in the calculation results of Grid 1. After considering the calculation accuracy and time cost, Grid 2 was selected.

In order to confirm the validity of the numerical calculation, the NACA4412 airfoil was taken as an example to compare the numerical calculation results by CFD with the results calculated by the XFOIL software and the experimental data [25]. The lift coefficient, drag coefficient, and lift to drag ratio C_l/C_d vary with the angle of attack at a Reynolds number of 6.4 × 10⁵ in the steady state, as shown in Figure 2. In terms of the lift coefficient, the curve trend obtained by numerical calculation and the XFOIL software is consistent with the experimental data, but the former is closer to the experimental value; especially if the angle of attack is less than 5°, the maximum error is about 3.55% (angle of attack: 0°). With the increase of the angle of attack, the error may be caused by airfoil stall or the increase of the turbulent kinetic energy, and the maximum error increases to about 9.71% (angle of attack: 10°). For the drag coefficient and lift to drag ratio, the numerical calculation results by CFD are evidently better than those calculated by the XFOIL software; its maximum drag coefficient error is about 8.60% (angle of attack: 7°). In general, the results showed that the numerical calculation results by CFD are accurate and agree well with the experimental data. The following numerical calculations for the airfoil in this study were performed using the same method mentioned above.

(a)

(b)

(c)

2.3. Optimized Object Selection

The original turbine is scanned by 3D coordinates, then the overall wind turbine blade model is obtained by the reverse engineering software UG, and the standard shape of the original airfoil (S2046) is obtained by the software Profili. The physical image and numerical model of the original wind turbine rotor are shown in Figure 3. As shown in Table 3, the rated power is 200 W, and the rated wind speed is 11.5 m/s.

(a)

(b)

The shapes of the S2046 airfoil and some other standard airfoils specially designed for wind turbines are shown in Figure 4, where x and y are the coordinates parallel to and perpendicular to the chord length, respectively.

A smaller blade size in SWTs results in a lower Reynolds number based on the blade section Re_r ranging from 1 × 10⁴ to 5 × 10⁵ [26], and the main aerodynamic characteristics of these Reynolds numbers in the steady state are shown in Figure 5. These characteristics reveal that the Reynolds number does not have a significant effect on the trend of the lift and drag coefficient when it changes in the range of 1.3 × 10⁵ to 2.1 × 10⁵. In general, the increase of the Reynolds number will increase the lift coefficient and lift to drag ratio, and the optimal angle of attack (α_opt) for airfoils is about 6°. Therefore, the optimization was performed at a Reynolds number of 1.3 × 10⁵ and an angle of attack of 6°. In addition, as mentioned in a previous study [27], the NACA airfoil is suitable for SWTs due to its high lift to drag ratio. Furthermore, the maximum lift to drag ratio of the NACA4412 airfoil is about 1.04 times that of the FFA-W1-128 airfoil and 1.53 times that of the S809 airfoil. Therefore, the airfoils commonly used in LWTs that operate at high Reynolds numbers are not necessarily suitable for SWTs.

(a)

(b)

(c)

(d)

Although the lift coefficient of the NACA4412 airfoil is almost similar to that of the S2046 airfoil and the lift to drag ratio of the S2046 airfoil is generally bigger, when the angle of attack is greater than 9°, the performance of the S2046 airfoil becomes unstable. Furthermore, the maximum relative thickness (9.02%) and the airfoil area (0.0568c²) of the S2046 airfoil are too small to be optimized. In other words, the S2046 airfoil is already optimized in terms of the original wind turbine. Considering the stability of the airfoil operation and the strength of the material, the NACA4412 airfoil is selected as the optimal design object, while the lift coefficient and lift to drag ratio are 0.931244 and 42.512657 at a Reynolds number of 1.3 × 10⁵ and an angle of attack of 6°, respectively.

2.4. Optimization Process

The main parameters of airfoil performance include lift coefficient and lift to drag ratio[28]. In general, the higher the lift coefficient and lift to drag ratio, the better the airfoil performance. Thus, the optimization objective is to maximize the lift coefficient and lift to drag ratio. The subflow value is used as the stopping criterion in this optimization algorithm. The design requirements are that the lift coefficient and lift to drag ratio are increased by at least 12 and 14% to the NACA4412 airfoil, respectively. In order to ensure that the material meets the requirement of strength and cost, the constraint conditions are that the maximum relative thickness should not decrease, and the airfoil area is increased by up to 15% to the S2046 airfoil. A MATLAB code is developed for optimization. The flowchart of the optimization process is shown in Figure 6.

2.5. Optimization Algorithm

As a random optimization method, the GA has good generalization, robustness, and portability, which makes it widely used in engineering optimization design [29]. A GA is proposed according to the law of biological evolution, which seeks the optimal solution through selection, crossover, and mutation. The flow of the GA used in the optimization is shown in Figure 7 [30, 31]. The GA is particularly effective for multiobjective optimization problems [32]. Many studies indicate that elitism (elite-preserving strategy) is a key to convergence and computational efficiency [33, 34]. Deb et al. [35] compared the performance of three elitist evolutionary algorithms, namely, NSGA-II, strength Pareto evolutionary algorithm (SPEA), and Pareto archived evolution strategy (PAES). They found that in most problems the NSGA-II is able to find much better spread of solutions and better convergence near the true Pareto-optimal front.

For NSGA-II options, the population size is 12, the number of generations is 20, the crossover probability is 0.9, the crossover distribution index is 10, and the mutation distribution index is 20. A total of 24 Pareto points were obtained by optimizing the airfoil. The Pareto front of the optimized airfoil, shown in Figure 8, reveals that there was no absolute optimal design. However, the lift coefficient and lift to drag ratio should reach their maximum values simultaneously. Therefore, the objective function value is the sum of the lift coefficient values and lift to drag ratio values. Thus, the Pareto point corresponding to the maximum value of the objective function is chosen.

In the case of the same subflow value (241), the latest GAs results, shown in Table 4, demonstrate that only NSGA-II meets the design requirements and the constraint conditions, while other GAs may take more time. Thus, the NSGA-II is used in this study.

2.6. Optimization Results

The shape of the optimized airfoil, which is named NACA4412-OPT, is shown in Figure 9, and the detailed parameters of the airfoils are shown in Table 5.

The main aerodynamic characteristics of the optimized airfoil at a Reynolds number of 2.1 × 10⁵ in the steady state, shown in Figure 10, reveal that the overall performance of the NACA4412 airfoil has been improved. The maximum lift coefficient of the NACA4412-OPT airfoil is about 1.08 times that of the NACA4412 airfoil and 1.1 times that of the S2046 airfoil. The maximum lift to drag ratio of the NACA4412-OPT airfoil is about 1.14 times that of the NACA4412 airfoil and 1.03 times that of the S2046 airfoil. Moreover, the NACA4412-OPT airfoil has a larger stall angle of attack than the S2046 airfoil, and its high lift coefficient and high lift to drag ratio can be maintained over a wider range of angle of attack. Overall, the optimized airfoil shows better steady aerodynamic characteristics.

(a)

(b)

(c)

2.7. Comparison of Unsteady Aerodynamic Characteristics

In order to further investigate the performance of the optimized airfoil, the unsteady state calculation of the airfoil in pitch motion was performed. The variation rule of angle of attack is described aswhere the average angle of attack α₀ = 9°, the pitch amplitude Δα = 9°, and the angular frequency of oscillation γ = 1.84 rad/s. Also, the position of rotation center x/c = 0.25 [36]; the numerical strategy is the same as that of the aforementioned steady computation; the time step size is 0.01s, and the number of time steps is 2,000. After 6 cycles of calculations, the monitoring parameters remained stable. Thus, the last pitching period was selected for statistical analysis.

The main aerodynamic characteristics of the optimized airfoil at a Reynolds number of 2.1 × 10⁵ in the unsteady state, shown in Figure 11, indicate that the variation law of the unsteady lift and drag coefficient is clearly different from that of the steady state. Also, as shown in Figure 2, the lift and drag coefficients in the steady state correspond to the angle of attack. However, for the unsteady state, the lift and drag coefficient form a closed curve during one period of the airfoil pitch. The lift coefficient is increased, and stall delay occurs. The drag coefficient changes quickly, especially when the airfoil is pulled down. The lift and drag coefficient values of the airfoil are higher and more stable when the airfoil is upturned than when it is pulled down. The NACA4412-OPT airfoil has a larger pitch range than the S2046 airfoil. When the airfoil is upturned, the lift coefficient of the NACA4412-OPT airfoil is greater than that of the S2046 airfoil. When the airfoil is pulled down, the drag coefficient of the NACA4412-OPT airfoil is smaller and more stable than that of the S2046 airfoil.

(a)

(b)

The pressure distribution of the airfoils at an angle of attack of 6.3° in both pitching conditions is shown in Figure 12. Although the angle of attack is the same, the differential pressure when the airfoil is upturned is larger than when it is pulled down. This explains why the lift curves of the airfoil do not coincide during the pitching motion. In addition, the greater the differential pressure, the greater the lift coefficient. This may be due to the fact that less streamline separation occurs on the surface when the airfoil is upturned.

(a)

(b)

The pressure around the NACA4412-OPT airfoil is lower during the pitching motion, especially on the upper surface. This may indicate that the velocity around the NACA4412-OPT airfoil is faster. When the NACA4412-OPT airfoil is upturned, the pressure on the lower surface slightly decreases. However, the pressure on the upper surface highly decreases, the differential pressure is larger, and therefore the lift coefficient is larger. When it is pulled down, the situation is reversed.

The turbulent kinematic energy distribution of the airfoils at an angle of attack of 6.3° in both pitching conditions is shown in Figure 13. Although the angle of attack is the same, the turbulent kinematic energy when the airfoil is pulled down is larger than when it is upturned, and the turbulence only occurs on the upper surface. This explains why the drag curves of the airfoil do not coincide during the pitching motion, and the greater the turbulent kinematic energy, the greater the drag coefficient. This may be due to the fact that the streamline changes on the upper surface are larger when the airfoil is pulled down.

(a)

(b)

The NACA4412-OPT airfoil has slightly more turbulent kinetic energy when it is upturned. Although the drag coefficient is increased, the contribution to the lift coefficient is greater. When the airfoil is pulled down, the turbulent kinematic energy of the S2046 airfoil is significantly larger. Thus, the S2046 airfoil has a larger drag coefficient, and the lift and drag coefficient highly fluctuates.

In general, the optimized airfoil has better unsteady aerodynamic characteristics.

3. New Wind Turbine Rotor Design

3.1. Design Principle

The power coefficient is limited to 0.593 by the Betz law [37], and it is defined aswhere is the radius of the wind wheel.

As shown in Table 6, the expected power is 255 W, and the estimated power coefficient is 0.33 according to (4). The new turbine blade is designed with the NACA4412-OPT airfoil.

The tip speed ratio is expressed aswhere is the angular velocity.

The design of a hub is more about its load characteristics; due to its small size, it has little influence on the final blade power output [38]. Therefore, this study is not concerned with the design of the shape of the hub, and it is replaced with a hemisphere with the radius of the hub and a cylinder, which is consistent with the size of the original wind turbine, namely, with . The new turbine design focuses on the blade design. The BEM theory has been usually used in the blade design [39]. Based on this theory, the influence of the vortex behind the rotor is considered by the Glauert design method [40]. The Wilson design method improves the Glauert design method, takes into account the energy loss caused by the blade tip and blade root, and is widely used in many wind turbine blade design models [41]. Therefore, this study uses the Wilson design method. The Prandtl correction factor F is described as [42]where and are the Prandtl tip and root correction factors at radius of the impeller respectively; is the inlet angle:where and are the axial and toroidal correction factors, respectively; is pitch angle. is the speed ratio at radius of impeller :

According to the BEM theory by ignoring the resistance and combining with equations (4), (5), and (9), we can get

As the rotor diameter is small, the loss of blade root is not considered:

The relationship between , , , , , and in Table 6 satisfies the following equation:

The calculation steps of the Wilson design method are as follows:(1)A section of the blades is taken every 25 mm along the radial direction, with a total of 23 sections(2)For each section, using (10) as the objective function and equations (7), (11), and (13) as the constraint conditions, the axial correction factor a and toroidal correction factors b and F of each blade element are calculated(3)The value of c is solved according to (12)(4)The value of is solved according to (8)

The solutions of a, b, and F can be obtained using MATLAB. In order to calculate the chord length of each section more accurately, the corresponding lift coefficient is inquired according to the Reynolds number of each section, while the determination of the Reynolds number depends on c; thus, the iterative method is used. The equation for the calculation of the Reynolds number of each section is as follows:where is the dynamic viscosity of air.

In addition, the curve graph of the lift and drag coefficient of the airfoil at different Reynolds number needs to be established. The lift coefficient corresponding to the Reynolds number of each section is determined by piecewise linear interpolation. Since there are too many lines, only some of them are shown in Figure 14.

(a)

(b)

As mentioned earlier, the optimal angle of attack for the NACA4412-OPT airfoil is about 6°, which can also be seen in Figure 14.

3.2. Design Results

The final design of the turbine is obtained based on the calculated results in Table 7. The output energy of the wind turbine is mainly contributed by the half part of the blade near the tip, but the calculated chord length (c_i,c_i+1) and pitch angle (θ′) at the blade root are too large, which leads to trouble in the fabrication. In practice, the blade shape can be simplified by drawing a straight line at points with 70 and 90% of the spread length, while the performance costs are small at high wind speeds [7, 43]. Therefore, the chord length and pitch angle should be determined by approximating linearly based on a distance of 30% from the blade tip [44], where c and θ are the optimum chord length and pitch angle, respectively. The chord length and pitch angle before and after modification are shown in Figure 15.

(a)

(b)

The point coordinates required to visualize a 3D model were obtained, and the frameworks were plotted as shown in Figure 16(a) and modeled as shown in Figure 16(b). The whole rotor is shown in Figure 16(c).

(a)

(b)

(c)

4. Wind Turbine Performance Comparison

4.1. Calculation Zones and Boundary Conditions

The rotation domain near the rotor and stationary domain showing the entire calculation zone is illustrated in Figure 17(a). Typically, a wind turbine blade needs a far-field domain of 10–20 times the rotor diameter [45]. However, considering the calculation cost, area immediately adjacent to the wake of the wind turbine [46], 3D helical vortex structure [47], empirical value [48–50], and comparison of the calculation zones of different sizes, the rotation and stationary domains are finally modeled into a cylinder shape with diameters of 1.1 and 4D and heights of 0.22 and 7.22D, respectively. The distance from the inlet to the rotation domain is 2D, and the distance from the outlet to the rotation domain is 5D. The rated speed condition is used as an inlet condition of the stationary domain, and the incoming velocity is 11.5 m/s. The boundary condition of the outlet in the stationary domain is free outflow. The boundary condition of the surface of the rotor and hub are set to a nonslip wall, but the latter rotates with the rotation domain. The interfaces between the rotation domain and stationary domain are set to sliding mesh. The surface outside the stationary domain is set to far field. The second upwind scheme and the second central scheme are used to solve convection and diffusion terms, respectively; the SIMPLE method is applied for the uncoupling of velocity and pressure. The k-ω SST model is used as the turbulence model [24]. Simulations are considered to be converged when the residual values are below 10⁻⁴.

(a)

(b)

(c)

(d)

In order to make the calculation faster and more accurate, a fully structured mesh, which is denser near the blades, is used. The mesh for the rotation domain, stationary domain, and surface of the rotor is shown in Figures 17(b)–17(d), respectively.

4.2. Reliability Check of the Numerical Simulation

First, tests of the mesh independence with mesh number and time step size are conducted. If y + on the surface of the rotor is less than 30, the number of mesh has little effect on the calculation results at the steady state, as shown in Table 8. The calculation results of Mesh 2 and Mesh 3 are very close, and Mesh 2 takes less time. There is a certain gap in the calculation results of Mesh 1. Considering the calculation accuracy and time cost, Mesh 2 is selected.

The area-weighted average dynamic pressure on the blade surface versus the time for different time step sizes, including Δt = 0.0001, 0.0005, and 0.001 s, is shown in Figure 18. The number of time steps is 10,000, 2,000, and 1,000, respectively, 1 second in total, which is equivalent to 16 rotor rotation cycles. The last 8 cycles were selected for statistical analysis. The area-weighted average dynamic pressure for Δt = 0.0005 and 0.001 s is 310.31 and 310.54 Pa, respectively, and the maximum fluctuation of the area-weighted average dynamic pressure for Δt = 0.0005 s and 0.001 s is 0.521 Pa and 0.533 Pa. It can be concluded that there was no obvious difference in calculation results for Δt = 0.0005 and 0.001 s. However, the maximum fluctuation for Δt = 0.0005 s is smaller; thus considering a compromise between accuracy and computation cost, Δt = 0.0005 s is used. The same selection of time step settings and statistical cycles is used below.

Second, in order to demonstrate the reliability of the numerical calculation, the numerical results were compared with the nominal power. The steady-state calculation is taken as the initial value of the transient state. The torque coefficient in relation to time, shown in Figure 19, reveals that the time-averaged torque coefficient value is 0.0556, and the maximum fluctuation of the torque coefficient is 0.0001132, 0.2% of the average value, which indicates that the torque output is very stable at constant velocity. The resulting power is 200.126 W. The error from the nominal value is only 0.063%, which indicates that the numerical results are reliable. At the same time, the time-averaged torque coefficient value (0.0556) at the transient state is quite different from the torque coefficient (0.031877) at the steady state, which indicates that the calculation results of the transient state are more consistent with the actual results.

4.3. Comparison of Aerodynamic Characteristics

In order to confirm the superiority of airfoil optimization, it is necessary to compare the aerodynamic characteristics of the two wind turbines at the transient state.

The power coefficient values of different rotors in relation to time, shown in Figure 20, reveal that the time-averaged power coefficient values of the original turbine and new turbine are 0.2585 and 0.3254, respectively. The power output is about 252 W, which indicates that the new turbine basically achieves the design goal and increases power by about 26%. The maximum fluctuation values of the original turbine and new turbine are 0.0005264 and 0.0005268, which are 0.2 and 0.16% of the average value, respectively, that the power output of the two turbines is also very stable at constant velocity.

The axial thrust coefficients values of different rotors in relation to time, shown in Figure 21, reveal that the time-averaged power coefficient values of the original turbine and new turbine are 0.6715 and 0.6751, respectively. The axial thrust output values are about 56.498 and 56.801 N, respectively, which indicates that the new turbine does not significantly increase the axial thrust. The maximum fluctuation values of the original turbine and new turbine are 0.00128 and 0.000563, which are 19.06 and 8.34% of the average value, respectively, indicating that the axial thrust output of the new turbine is more stable than that of the original turbine at constant velocity.

The pressure distribution on the surface of the two rotors at t = 0.5 s, shown in Figure 22, reveals that the pressure on the pressure surface is basically positive, and the pressure on the suction surface is almost negative. The pressure values are higher near the leading edge, especially at the tip of the blades, which confirms that the output power of the wind turbine mainly depends on the tip of the blade [44]. The pressure contour variable range on the surface of the original turbine and new turbine is from −2,865.36 to 1,629.63 Pa and −3,267.02 to 1,696.75 Pa. The larger differential pressure on the blade surface of the new turbine is also the reason for the larger power coefficient and axial thrust coefficient.

(a)

(b)

The turbulent kinematic energy distribution on the surface of the two rotors at t = 0.5 s, shown in Figure 23, indicates that the turbulent kinematic energy on the pressure surface is basically smaller than that on the suction surface. The turbulent kinematic energy values are higher near the leading edge, especially at the tip of the blades. The turbulent kinematic energy contour variable ranges on the surface of the original turbine and the new turbine are 0.000181553 to 73.3571 m²/s² and 0.000206362 to 90.7744 m²/s², respectively. The turbulent kinematic energy values on the blade surface of the new turbine are larger but more well-distributed. This may be the reason for the more stable operation of the new turbine.

(a)

(b)

The pathlines colored by the turbulent kinematic energy around the two rotors, shown in Figure 24, show that the pathlines expand when passing through the turbine, and then the downstream pathlines slightly shrink and eventually become nearly parallel to the rotation axis. The pathlines at the hub have a discernible spiral shape. The pathlines of the new turbine rotor at the hub are smaller, while the pathlines of the original turbine rotor at the far downstream are more turbulent, which indicates that the flow field of the new turbine is more stable.

(a)

(b)

5. Conclusions

Currently, most of the studies on the multiobjective optimization of wind turbine rotor focus on the configuration of the blade and the efficient optimization algorithm (NSGA-II is mostly used for LWTs, e.g.), which may be one of the limitations of the SWTs design process. Therefore, this paper develops a horizontal axis SWT optimization technique through the multiobjective optimization of the airfoil in order to improve the efficiency of power generation for local renewable energy applications in China. The NACA4412 airfoil is considered as a baseline for the optimization using NSGA-II, which outperforms the other commonly used GAs. The optimized airfoil (NACA4412-OPT) has a larger maximum lift coefficient, larger maximum lift to drag ratio, larger stall angle of attack, and wider range of angle of attack. Its high lift coefficient and high lift to drag ratio can be maintained over the baseline airfoil as well as the original airfoil (S2046). It shows better steady and unsteady aerodynamic characteristics while avoiding excessive variations in the maximum relative thickness and area. It is used as cross-sectional profile in the design procedure for a new 1.15 m diameter three-bladed wind turbine rotor at a wind speed of 11.5 m/s, tip speed ratio of 4.65, and pitch angle of 0.2° by the Wilson design method. Finally, the unsteady simulations of the whole turbine are performed to analyze the flow field around the rotor and compare it with the original wind turbine rotor. The numerical simulation results based on CFD demonstrate that the time-averaged power coefficient value (0.33) of the new wind turbine is almost 26% higher and more stable than that of the original wind turbine while avoiding the significant increase of the axial thrust.

Data Availability

In order to prove the validity of the numerical calculation, NACA4412 was taken as an example to compare the numerical calculation results with the results calculated by XFOIL software and the experimental data [25].

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The work described in this paper was funded by the National Natural Science Foundation of China (no. 11862024) and the Scientific Research Fund Project of the Education Department of Yunnan Province (no. 2021J0823). The authors are grateful for their financial support.

References

N. Karanikas, S. Steele, K. Bruschi et al., “Occupational health hazards and risks in the wind industry,” Energy Reports, vol. 7, pp. 3750–3759, 2021.
View at: Publisher Site | Google Scholar
J. Leary, M. Czyrnek-Delêtre, A. Alsop et al., “Finding the niche: a review of market assessment methodologies for rural electrification with small scale wind power,” Renewable and Sustainable Energy Reviews, vol. 133, p. 110240, 2020.
View at: Publisher Site | Google Scholar
T. M. Sugathapala, S. Boteju, P. B. Withanage, and S. Wijewardane, “Aerodynamic modeling of simplified wind turbine rotors targeting small-scale applications in Sri Lanka,” Energy for Sustainable Development, vol. 59, pp. 71–82, 2020.
View at: Publisher Site | Google Scholar
A. M. Abdelsalam, W. A. El-Askary, M. A. Kotb, and I. M. Sakr, “Experimental study on small scale horizontal axis wind turbine of analytically-optimized blade with linearized chord twist angle profile,” Energy, vol. 216, p. 119304, 2021.
View at: Publisher Site | Google Scholar
M.-H. Lee, Y. C. Shiah, and C.-J. Bai, “Experiments and numerical simulations of the rotor-blade performance for a small-scale horizontal axis wind turbine,” Journal of Wind Engineering and Industrial Aerodynamics, vol. 149, pp. 17–29, 2016.
View at: Publisher Site | Google Scholar
A. Pourrajabian, P. A. Nazmi Afshar, M. Ahmadizadeh, and D. Wood, “Aero-structural design and optimization of a small wind turbine blade,” Renewable Energy, vol. 87, pp. 837–848, 2016.
View at: Publisher Site | Google Scholar
S. Rahgozar, A. Pourrajabian, S. A. A. Kazmi, and S. M. R. Kazmi, “Performance analysis of a small horizontal axis wind turbine under the use of linear/nonlinear distributions for the chord and twist angle,” Energy for Sustainable Development, vol. 58, pp. 42–49, 2020.
View at: Publisher Site | Google Scholar
R. K. Gupta, V. Warudkar, R. Purohit, and S. Singh Rajpurohit, “Modeling and aerodynamic analysis of small scale, mixed airfoil horizontal axis wind turbine blade,” Materials Today Proceedings, vol. 4, no. 4, pp. 5370–5384, 2017.
View at: Publisher Site | Google Scholar
M. Hasan, A. El-Shahat, and M. Rahman, “Performance investigation of three combined airfoils bladed small scale horizontal axis wind turbine by BEM and CFD analysis,” Journal of Power and Energy Engineering, vol. 05, no. 05, pp. 14–27, 2017.
View at: Publisher Site | Google Scholar
H. Muhsen, W. Al-Kouz, and W. Khan, “Small wind turbine blade design and optimization,” Symmetry, vol. 12, no. 1, p. 18, 2019.
View at: Publisher Site | Google Scholar
N. Pholdee, S. SujinBureerat, and W. Nuantong, “Kriging surrogate-based genetic algorithm optimization for blade design of a horizontal AxisWind turbine,” Computer Modeling in Engineering and Sciences, vol. 126, no. 1, pp. 261–273, 2021.
View at: Publisher Site | Google Scholar
X. Shen, H. Yang, J. Chen, X. Zhu, and Z. Du, “Aerodynamic shape optimization of non-straight small wind turbine blades,” Energy Conversion and Management, vol. 119, pp. 266–278, 2016.
View at: Publisher Site | Google Scholar
Y. Yang, C. Li, W. Zhang et al., “A multi-objective optimization for HAWT blades design by considering structural strength,” Journal of Mechanical Science and Technology, vol. 30, no. 8, pp. 3693–3703, 2016.
View at: Publisher Site | Google Scholar
A. Maheri, “Multiobjective optimisation and integrated design of wind turbine blades using WTBM-ANSYS for high fidelity structural analysis,” Renewable Energy, vol. 145, pp. 814–834, 2020.
View at: Publisher Site | Google Scholar
J. L. Moreira, A. L. Mesquita, L. F. Araujo, M. A. Galhardo, J. R. Vaz, and J. T. Pinho, “Experimental investigation of drivetrain resistance applied to small wind turbines,” Renewable Energy, vol. 153, pp. 324–333, 2020.
View at: Publisher Site | Google Scholar
M. Refan and H. Hangan, “Aerodynamic performance of a small horizontal axis wind turbine,” Journal of Solar Energy Engineering, vol. 134, no. 2, 2012.
View at: Publisher Site | Google Scholar
B. Plaza, R. Bardera, and S. Visiedo, “Comparison of BEM and CFD results for Mexico rotor aerodynamics,” Journal of Wind Engineering and Industrial Aerodynamics, vol. 145, pp. 115–122, 2015.
View at: Publisher Site | Google Scholar
M. O. Moussa, “Experimental and numerical performances analysis of a small three blades wind turbine,” Energy, vol. 203, p. 117807, 2020.
View at: Publisher Site | Google Scholar
K. D. Lee and S. Eyi, “Aerodynamic design via optimization,” Journal of Aircraft, vol. 29, no. 6, pp. 1012–1019, 1992.
View at: Publisher Site | Google Scholar
J. Samareh, “Aerodynamic shape optimization based on free-form deformation,” in Proceedings of the 10th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, p. 4630, Albany, New York, September 2004.
View at: Google Scholar
R. M. Hicks and P. A. Henne, “Wing design by numerical optimization,” Journal of Aircraft, vol. 15, no. 7, pp. 407–412, 1978.
View at: Publisher Site | Google Scholar
S. Zhou, H. Zhou, K. Yang, H. Dong, and Z. Gao, “Research on blade design method of multi-blade centrifugal fan for building efficient ventilation based on Hicks-Henne function,” Sustainable Energy Technologies and Assessments, vol. 43, p. 100971, 2021.
View at: Publisher Site | Google Scholar
Z. Xiaoping, D. Jifeng, L. Weiji, and Z. Yong, “Robust airfoil optimization with multi-objective estimation of distribution algorithm,” Chinese Journal of Aeronautics, vol. 21, no. 4, pp. 289–295, 2008.
View at: Publisher Site | Google Scholar
F. R. Menter, “Two-equation eddy-viscosity turbulence models for engineering applications,” AIAA Journal, vol. 32, no. 8, pp. 1598–1605, 1994.
View at: Publisher Site | Google Scholar
E. N. Jacobs and A. Sherman, “Airfoil section characteristics as affected by variations of the Reynolds number,” NACA Technical Report, vol. 586, pp. 227–267, 1937.
View at: Google Scholar
Y. Liu, P. Li, W. He, and K. Jiang, “Numerical study of the effect of surface grooves on the aerodynamic performance of a NACA 4415 airfoil for small wind turbines,” Journal of Wind Engineering and Industrial Aerodynamics, vol. 206, p. 104263, 2020.
View at: Publisher Site | Google Scholar
A. Musyafa, A. Adi, H. Justiono, and H. Cordova, “Implementation of fuzzy logic control (FLC) in horizontal Axis wind turbine prototype with airfoil profile NREL standard S83 at low rate wind speed,” Australian Journal of Basic and Applied Sciences, vol. 8, no. 13, pp. 61–68, 2014.
View at: Google Scholar
A. Björck, Coordinates and Caluclations for the FFA-W1-Xxx, FFA-W2-Xxx and FFA-W3-Xxx Series of Airfoils for Horizontal Axis Wind Turbines, Flygtekniska Försöksanstalten, the Aeronautical Research Institute of Sweden, Ulvsunda, Sweden, 1990.
A. Belegundu, E. Constans, R. Salagame, and D. Murthy, “Multi-objective optimization of laminated ceramic composites using genetic algorithms,” in Proceedings of the 5th Symposium on Multidisciplinary Analysis and Optimization, p. 4363, Panama City Beach, FL, September, 1994.
View at: Google Scholar
D. E. Goldberg, Genetic Algorithms in Search Optimization and Machine learning, Addison-Wesley Publishing Company Inc, Boston, 1989.
M. Mitchell, An Introduction to Genetic Algorithms, MIT press, Cambridge, MA, 1998.
K. Deb, “Multi-objective genetic algorithms: problem difficulties and construction of test problems,” Evolutionary Computation, vol. 7, no. 3, pp. 205–230, 1999.
View at: Publisher Site | Google Scholar
E. Zitzler, K. Deb, and L. Thiele, “Comparison of multiobjective evolutionary algorithms: empirical results,” Evolutionary Computation, vol. 8, no. 2, pp. 173–195, 2000.
View at: Publisher Site | Google Scholar
J. D. Knowles and D. W. Corne, “Approximating the nondominated front using the Pareto archived evolution strategy,” Evolutionary Computation, vol. 8, no. 2, pp. 149–172, 2000.
View at: Publisher Site | Google Scholar
K. Deb, S. Agrawal, A. Pratap, and T. Meyarivan, “A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: nsga-II, Parallel Problem Solving from Nature PPSN VI,” in Proceedings of the International Conference on Parallel Problem Solving from Nature, pp. 849–858, Springer, Paris, France, September, 2000.
View at: Publisher Site | Google Scholar
A. Björck and S. E. Thor, Dynamic Stall and 3D Effects, The Aeronautical Research Institute of Sweden, Bromma, Sweden, 1995.
A. Betz, Schraubenpropeller mit geringstem Energieverlust, Gottinger Nachrichten, Göttingen, Germany, pp. 193–213, 1919.
P. Freere, M. Sacher, J. Derricott, and B. Hanson, “A low cost wind turbine and blade performance,” Wind Engineering, vol. 34, no. 3, pp. 289–302, 2010.
View at: Publisher Site | Google Scholar
J. F. Manwell, J. G. McGowan, and A. L. Rogers, Wind Energy Explained: Theory, Design and Application, John Wiley & Sons, New Jersey, USA, 2010.
H. Glauert, “Airplane propellers,” Aerodynamic Theory, Springer, Berlin, Heidelberg, vol. 4, pp. 169–269, 1936.
View at: Google Scholar
R. E. Wilson and P. Lissaman, Applied Aerodynamics of Wind Power Machines, Oregon State University Corvallis, Corvallis, OR, USA, 1974.
L. Prandtl, Applied hydro-and aeromechanics: based on lectures of L. Prandtl, Dover Publications, Mineola, New York, 1957.
T. Burton, N. Jenkins, D. Sharpe, and E. Bossanyi, Wind Energy Handbook, John Wiley & Sons, New Jersey, USA, 2011.
C. h. Jo, J. y. Yim, K. h. Lee, and Y. h. Rho, “Performance of horizontal axis tidal current turbine by blade configuration,” Renewable Energy, vol. 42, pp. 195–206, 2012.
View at: Publisher Site | Google Scholar
A. P. Schaffarczyk, Introduction to Wind Turbine Aerodynamics, Springer Nature, Berlin, Heidelberg, 2020.
L. J. Vermeer, J. N. Sørensen, and A. Crespo, “Wind turbine wake aerodynamics,” Progress in Aerospace Sciences, vol. 39, no. 6-7, pp. 467–510, 2003.
View at: Publisher Site | Google Scholar
V. L. Okulov, I. V. Naumov, R. F. Mikkelsen, I. K. Kabardin, and J. N. Sørensen, “A regular Strouhal number for large-scale instability in the far wake of a rotor,” Journal of Fluid Mechanics, vol. 747, pp. 369–380, 2014.
View at: Publisher Site | Google Scholar
C. Yigit, “Effect of air-ducted blade design on horizontal axis wind turbine performance,” Energies, vol. 13, no. 14, p. 3618, 2020.
View at: Publisher Site | Google Scholar
W.-Q. Wang, R. Yin, and Y. Yan, “Design and prediction hydrodynamic performance of horizontal axis micro-hydrokinetic river turbine,” Renewable Energy, vol. 133, pp. 91–102, 2019.
View at: Publisher Site | Google Scholar
M. Carrión, R. Steijl, M. Woodgate, G. N. Barakos, X. Munduate, and S. Gomez-Iradi, “Aeroelastic analysis of wind turbines using a tightly coupled CFD-CSD method,” Journal of Fluids and Structures, vol. 50, pp. 392–415, 2014.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2022 Rui Yin et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1206

Downloads

767

Citations