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

Authors: Rui Yin; Jian-Bin Xie; Ji Yao
Journal: Mathematical Problems in Engineering
 (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.


---

## Body

## 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 Table1.Table 1
Summary of the reviewed studies investigating blades.
ReferenceBlade typeModel numbers of the airfoilScaleOptimization goalsOptimization methodsOptimization variablesAerodynamics prediction methods[3]UniformNACA4412Small———BEM[4]UniformRISØ-A-24Small———ExperimentBEM[5]UniformSD8000Small———ExperimentBEMCFD (Fluent)[6]UniformSG6043SmallStarting timeGAChordBEMOutput powerTwist angleBlade allowable stressShell thickness[7]UniformSG6043SmallStarting timeGAChordBEMOutput powerTwist angle[8]MixedSG6050, SG6043, E555, and E216Small———BEMCFD (Fluent)[9]MixedS823, S833, and S822Small———BEMCFD (CFX)[10]MixedS1210 and S1223Small———BEMCFD (QBlade)[11]ModifiedDU91-W2-250 and RISØ-A1-21SmallPower coefficientGAChordCFD (Fluent)NACA 64-418Twist angle[12]ModifiedS809SmallAnnual energy productionGAChordCFDStarting performanceTwist angle(A lifting surface method)3D stacking line[13]ModifiedDU series airfoilsLargeAnnual energy productionNSGA-IIChordBEMNACA 64-618Blade massTwist angle[14]ModifiedDU series airfoilsLargeBlade massNSGA-IIBlade sizeCFD (ANSYS)NACA 64 A17Maximum von Mises stressTopologyBlade tip deflectionStructureMaterialThe 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]:(1)yu=yu0+∑i=1n/2δfixi,yl=yl0+∑i=n/2+1nδifix,where yu and yl are coordinates on the top and bottom surfaces of the new airfoil, respectively; yu0 and yl0 are the coordinates on the top and bottom surfaces of the initial airfoil, respectively; n is the number of design variables; fix are the Hicks–Henne shape functions; δi 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] δ1 ∈ [–0.006, 0.006], δ2 ∈ [–0.01, 0.006], δ3 ∈ [–0.006, 0.008], δ4 ∈ [–0.005, 0.005], δ5 ∈ [–0.005, 0.01], δ6 ∈ [–0.008, 0.006], δ7 ∈ [–0.01, 0.01], δ8 ∈ [–0.005, 0.01], δ9 ∈ [–0.005, 0.005], and δ10 ∈ [–0.005, 0.015].The Hicks–Henne shape functions are(2)fix=x0.251−xe−20x,i=1,6,sin3πxei,i=2,⋯,5,7,⋯,10,ei=ln0.5lnxi,xi=0.2,0.4,0.6,0.8,0.2,0.4,0.6,0.8,i=2,⋯,5,7,⋯,10,where xi 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−6.Figure 1
Computational zone and grids: (a) whole grids and (b) finest local grid around hydrofoil.
(a)(b)Ify + 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 × 105 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.Table 2
The calculation results with different number of grids (Re = 6.4 × 105).
ItemCellsConvergence stepsLift coefficientClDrag coefficientCdGrid 1112703000.782400.0152Grid 2273703500.776590.01459Grid 3411705000.774630.01445In 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 Cl/Cd vary with the angle of attack at a Reynolds number of 6.4 × 105 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.Figure 2
Lift coefficient, drag coefficient, and lift to drag ratio of the NACA4412 airfoil (Re = 6.4 × 105): (a) lift coefficient, (b) drag coefficient, and (c) lift to drag ratio.
(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 Figure3. As shown in Table 3, the rated power is 200 W, and the rated wind speed is 11.5 m/s.Figure 3
The original wind turbine rotor: (a) physical image and (b) numerical model.
(a)(b)Table 3
The original turbine blade design parameters.
Design parametersValuesB: blade number (EA)3D: rotor diameter (m)1.15P: rated power (W)200V1: rated wind speed (m/s)11.5N: angular speed (rpm)880ρ: air density (kg/m3)1.225Cp: estimated power coefﬁcient0.259The shapes of the S2046 airfoil and some other standard airfoils specially designed for wind turbines are shown in Figure4, where x and y are the coordinates parallel to and perpendicular to the chord length, respectively.Figure 4
The shape of standard airfoils.A smaller blade size in SWTs results in a lower Reynolds number based on the blade section Rer ranging from 1 × 104 to 5 × 105 [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 × 105 to 2.1 × 105. 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 × 105 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.Figure 5
The main aerodynamic characteristics of standard airfoils in a steady state: (a) lift coefficient (Re = 1.3 × 105); (b) lift to drag ratio (Re = 1.3 × 105); (c) lift coefficient (Re = 2.1 × 105); (d) lift to drag ratio (Re = 2.1 × 105).
(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.0568c2) 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 × 105 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.Figure 6
Optimization process.


### 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.Figure 7
The flow of the GA.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 Figure8, 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.Figure 8
The Pareto front of the optimized airfoil.In the case of the same subflow value (241), the latest GAs results, shown in Table4, 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.Table 4
The detailed parameters of airfoils.
Airfoil/algorithmCl (rate of increase)Cl/Cd (rate of increase)Maximum relative thickness (%)Area (rate of increase)NAGA-II1.0506 (12.82%)48.765 (14.71%)10.160.0649c2 (14.26%)NCGA1.0381 (11.47%)48.422 (13.9%)10.310.0672c2 (18.31%)AMGA1.0327 (10.89%)48.148 (13.26%)10.950.0695c2 (22.36%)Multi-Island GA0.97855 (5.08%)46.271 (8.84%)11.110.0723c2 (27.29%)


### 2.6. Optimization Results

The shape of the optimized airfoil, which is named NACA4412-OPT, is shown in Figure9, and the detailed parameters of the airfoils are shown in Table 5.Figure 9
The shape of optimized airfoil.Table 5
The detailed parameters of airfoils.
AirfoilsMaximum relative camber (position)Maximum relative thickness (position)Leading edge radiusAreaS20469.02% (28.6%)2.58% (54.5%)0.8007%0.0568c2NACA4412-OPT10.16% (26.1%)4.24% (37.5%)1.37%0.0649c2NACA441212% (30%)4% (40%)%0.081c2The main aerodynamic characteristics of the optimized airfoil at a Reynolds number of 2.1 × 105 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.Figure 10
The main aerodynamic characteristics of the optimized airfoil in a steady state (Re = 2.1 × 105): (a) lift coefficient, (b) drag coefficient, and (c) lift to drag ratio.
(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 as(3)α=α0+Δαsinλt,where the average angle of attack α0 = 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 × 105 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.Figure 11
The main aerodynamic characteristics of the optimized airfoil in an unsteady state (Re = 2.1 × 105): (a) lift coefficient; (b) drag coefficient.
(a)(b)The pressure distribution of the airfoils at an angle of attack of 6.3° in both pitching conditions is shown in Figure12. 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.Figure 12
Pressure distribution of airfoils (Pa): (a) upturned; (b) pulled down.
(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 Figure13. 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.Figure 13
Turbulent kinematic energy distribution of airfoils (m2/s2): (a) upturned; (b) 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.


## 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]:(1)yu=yu0+∑i=1n/2δfixi,yl=yl0+∑i=n/2+1nδifix,where yu and yl are coordinates on the top and bottom surfaces of the new airfoil, respectively; yu0 and yl0 are the coordinates on the top and bottom surfaces of the initial airfoil, respectively; n is the number of design variables; fix are the Hicks–Henne shape functions; δi 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] δ1 ∈ [–0.006, 0.006], δ2 ∈ [–0.01, 0.006], δ3 ∈ [–0.006, 0.008], δ4 ∈ [–0.005, 0.005], δ5 ∈ [–0.005, 0.01], δ6 ∈ [–0.008, 0.006], δ7 ∈ [–0.01, 0.01], δ8 ∈ [–0.005, 0.01], δ9 ∈ [–0.005, 0.005], and δ10 ∈ [–0.005, 0.015].The Hicks–Henne shape functions are(2)fix=x0.251−xe−20x,i=1,6,sin3πxei,i=2,⋯,5,7,⋯,10,ei=ln0.5lnxi,xi=0.2,0.4,0.6,0.8,0.2,0.4,0.6,0.8,i=2,⋯,5,7,⋯,10,where xi 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−6.Figure 1
Computational zone and grids: (a) whole grids and (b) finest local grid around hydrofoil.
(a)(b)Ify + 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 × 105 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.Table 2
The calculation results with different number of grids (Re = 6.4 × 105).
ItemCellsConvergence stepsLift coefficientClDrag coefficientCdGrid 1112703000.782400.0152Grid 2273703500.776590.01459Grid 3411705000.774630.01445In 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 Cl/Cd vary with the angle of attack at a Reynolds number of 6.4 × 105 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.Figure 2
Lift coefficient, drag coefficient, and lift to drag ratio of the NACA4412 airfoil (Re = 6.4 × 105): (a) lift coefficient, (b) drag coefficient, and (c) lift to drag ratio.
(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 Figure3. As shown in Table 3, the rated power is 200 W, and the rated wind speed is 11.5 m/s.Figure 3
The original wind turbine rotor: (a) physical image and (b) numerical model.
(a)(b)Table 3
The original turbine blade design parameters.
Design parametersValuesB: blade number (EA)3D: rotor diameter (m)1.15P: rated power (W)200V1: rated wind speed (m/s)11.5N: angular speed (rpm)880ρ: air density (kg/m3)1.225Cp: estimated power coefﬁcient0.259The shapes of the S2046 airfoil and some other standard airfoils specially designed for wind turbines are shown in Figure4, where x and y are the coordinates parallel to and perpendicular to the chord length, respectively.Figure 4
The shape of standard airfoils.A smaller blade size in SWTs results in a lower Reynolds number based on the blade section Rer ranging from 1 × 104 to 5 × 105 [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 × 105 to 2.1 × 105. 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 × 105 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.Figure 5
The main aerodynamic characteristics of standard airfoils in a steady state: (a) lift coefficient (Re = 1.3 × 105); (b) lift to drag ratio (Re = 1.3 × 105); (c) lift coefficient (Re = 2.1 × 105); (d) lift to drag ratio (Re = 2.1 × 105).
(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.0568c2) 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 × 105 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.Figure 6
Optimization process.


## 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.Figure 7
The flow of the GA.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 Figure8, 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.Figure 8
The Pareto front of the optimized airfoil.In the case of the same subflow value (241), the latest GAs results, shown in Table4, 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.Table 4
The detailed parameters of airfoils.
Airfoil/algorithmCl (rate of increase)Cl/Cd (rate of increase)Maximum relative thickness (%)Area (rate of increase)NAGA-II1.0506 (12.82%)48.765 (14.71%)10.160.0649c2 (14.26%)NCGA1.0381 (11.47%)48.422 (13.9%)10.310.0672c2 (18.31%)AMGA1.0327 (10.89%)48.148 (13.26%)10.950.0695c2 (22.36%)Multi-Island GA0.97855 (5.08%)46.271 (8.84%)11.110.0723c2 (27.29%)


## 2.6. Optimization Results

The shape of the optimized airfoil, which is named NACA4412-OPT, is shown in Figure9, and the detailed parameters of the airfoils are shown in Table 5.Figure 9
The shape of optimized airfoil.Table 5
The detailed parameters of airfoils.
AirfoilsMaximum relative camber (position)Maximum relative thickness (position)Leading edge radiusAreaS20469.02% (28.6%)2.58% (54.5%)0.8007%0.0568c2NACA4412-OPT10.16% (26.1%)4.24% (37.5%)1.37%0.0649c2NACA441212% (30%)4% (40%)%0.081c2The main aerodynamic characteristics of the optimized airfoil at a Reynolds number of 2.1 × 105 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.Figure 10
The main aerodynamic characteristics of the optimized airfoil in a steady state (Re = 2.1 × 105): (a) lift coefficient, (b) drag coefficient, and (c) lift to drag ratio.
(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 as(3)α=α0+Δαsinλt,where the average angle of attack α0 = 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 × 105 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.Figure 11
The main aerodynamic characteristics of the optimized airfoil in an unsteady state (Re = 2.1 × 105): (a) lift coefficient; (b) drag coefficient.
(a)(b)The pressure distribution of the airfoils at an angle of attack of 6.3° in both pitching conditions is shown in Figure12. 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.Figure 12
Pressure distribution of airfoils (Pa): (a) upturned; (b) pulled down.
(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 Figure13. 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.Figure 13
Turbulent kinematic energy distribution of airfoils (m2/s2): (a) upturned; (b) 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 coefficientCP is limited to 0.593 by the Betz law [37], and it is defined as(4)CP=P1/2πρR2V13,where R is the radius of the wind wheel.As shown in Table6, 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.Table 6
The new turbine blade design parameters.
Design parametersValuesB: blade number (EA)3D: rotor diameter (m)1.15P: expected power (W)255V1: rated wind speed (m/s)11.5Cp: estimated power coefﬁcient0.33λ0: tip speed ratio4.65η: estimated efﬁciency coefﬁcient0.8The tip speed ratio is expressed as(5)λ0=ωRV1,where ω 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 Rhub and a cylinder, which is consistent with the size of the original wind turbine, namely, with Rhub=65mm. 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](6)Ftip=2πarccoseBr−R/2RsinφFhub=2πarccoseBRhub−r/2RsinφF=FtipFhub,where Ftip and Fhub are the Prandtl tip and root correction factors at radius r of the impeller respectively; φ is the inlet angle:(7)φ=arctan1−aλ1+b,(8)θ=φ−αopt,where a and b are the axial and toroidal correction factors, respectively; θ is pitch angle. λ is the speed ratio at radius of impeller r:(9)λ=ωrV1.According to the BEM theory by ignoring the resistance and combining with equations (4), (5), and (9), we can get(10)dCPdλ=8λ02b1−aFλ3,(11)bb+1λ2=a1−aF,(12)c=8πraF1−aFsin2φBCl1−a2cosφ.As the rotor diameter is small, the loss of blade root is not considered:(13)F=Ftip=2πarccoseBr−R/2Rsinφ.The relationship betweenD, P, Cp, η, ρ, and V1 in Table 6 satisfies the following equation:(14)D=PCPηρ/2V13π/4.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 ofc is solved according to (12)(4)
The value ofθ is solved according to (8)The solutions ofa, 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:(15)Re=ρV11−acμsinφ,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 Figure14.Figure 14
The main aerodynamic characteristics of the NACA4412-OPT airfoil: (a) lift coefficient; (b) lift to drag ratio.
(a)(b)As mentioned earlier, the optimal angle of attack for the NACA4412-OPT airfoil is about 6°, which can also be seen in Figure14.


### 3.2. Design Results

The final design of the turbine is obtained based on the calculated results in Table7. 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 (ci,ci+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.Table 7
The new turbine blade design results.
r (mm)λabFφ (deg)θ′ (deg)Re/ × 104ci (mm)ci+1 (mm)Δ (mm)c (mm)θ (deg)250.220.301.870.9049.9143.915.6272.9472.940.000695.4324.02500.430.320.770.9142.7336.739.09109.56109.560.000192.4822.52750.650.330.430.9236.6830.6811.33122.61122.611 × 10−689.5321.021000.870.330.270.9331.7225.7212.79122.71122.716 × 10−786.5819.521251.090.330.190.9427.6921.6913.74116.82116.822 × 10−583.6318.021501.300.340.140.9524.4318.4314.38108.74108.74−1 × 10−580.6816.521751.520.340.100.9621.7715.7714.82100.26100.264 × 10−577.7315.022001.740.340.080.9619.5813.5815.1392.1892.183 × 10−674.7813.522251.960.340.060.9617.7511.7515.3684.8084.809 × 10−671.8311.652502.170.340.050.9616.2010.2015.5378.1978.194 × 10−668.8810.112752.390.340.040.9614.898.8915.6672.3072.302 × 10−665.938.833002.610.340.040.9613.757.7515.7667.0767.076 × 10−562.997.743252.830.340.030.9612.766.7615.8462.4262.42−3 × 10−560.046.813503.040.350.030.9611.895.8915.9058.2758.27−8 × 10−757.095.993753.260.350.020.9511.125.1215.9454.5454.54−9 × 10−754.145.254003.480.350.020.9410.424.4215.9751.1751.17−1 × 10−651.194.554253.700.350.020.939.783.7815.9848.1048.10−1 × 10−648.243.884503.910.360.020.919.183.1815.9745.2745.27−9 × 10−745.293.224754.130.370.020.888.602.6015.9142.5842.58−8 × 10−742.342.575004.350.380.020.838.022.0215.7839.9239.92−8 × 10−539.391.935254.570.410.020.757.371.3715.4737.0137.01−2 × 10−736.441.315504.780.450.020.606.500.5014.5832.8232.822 × 10−633.490.725755—————————30.540.20Figure 15
Chord length and pitch angle before and after modification: (a) chord length; (b) pitch angle.
(a)(b)The point coordinates required to visualize a 3D model were obtained, and the frameworks were plotted as shown in Figure16(a) and modeled as shown in Figure 16(b). The whole rotor is shown in Figure 16(c).Figure 16
The new wind turbine rotor: (a) cross-sectional profiles of all sections; (b) one blade; (c) whole turbine rotor.
(a)(b)(c)


## 3.1. Design Principle

The power coefficientCP is limited to 0.593 by the Betz law [37], and it is defined as(4)CP=P1/2πρR2V13,where R is the radius of the wind wheel.As shown in Table6, 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.Table 6
The new turbine blade design parameters.
Design parametersValuesB: blade number (EA)3D: rotor diameter (m)1.15P: expected power (W)255V1: rated wind speed (m/s)11.5Cp: estimated power coefﬁcient0.33λ0: tip speed ratio4.65η: estimated efﬁciency coefﬁcient0.8The tip speed ratio is expressed as(5)λ0=ωRV1,where ω 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 Rhub and a cylinder, which is consistent with the size of the original wind turbine, namely, with Rhub=65mm. 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](6)Ftip=2πarccoseBr−R/2RsinφFhub=2πarccoseBRhub−r/2RsinφF=FtipFhub,where Ftip and Fhub are the Prandtl tip and root correction factors at radius r of the impeller respectively; φ is the inlet angle:(7)φ=arctan1−aλ1+b,(8)θ=φ−αopt,where a and b are the axial and toroidal correction factors, respectively; θ is pitch angle. λ is the speed ratio at radius of impeller r:(9)λ=ωrV1.According to the BEM theory by ignoring the resistance and combining with equations (4), (5), and (9), we can get(10)dCPdλ=8λ02b1−aFλ3,(11)bb+1λ2=a1−aF,(12)c=8πraF1−aFsin2φBCl1−a2cosφ.As the rotor diameter is small, the loss of blade root is not considered:(13)F=Ftip=2πarccoseBr−R/2Rsinφ.The relationship betweenD, P, Cp, η, ρ, and V1 in Table 6 satisfies the following equation:(14)D=PCPηρ/2V13π/4.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 ofc is solved according to (12)(4)
The value ofθ is solved according to (8)The solutions ofa, 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:(15)Re=ρV11−acμsinφ,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 Figure14.Figure 14
The main aerodynamic characteristics of the NACA4412-OPT airfoil: (a) lift coefficient; (b) lift to drag ratio.
(a)(b)As mentioned earlier, the optimal angle of attack for the NACA4412-OPT airfoil is about 6°, which can also be seen in Figure14.


## 3.2. Design Results

The final design of the turbine is obtained based on the calculated results in Table7. 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 (ci,ci+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.Table 7
The new turbine blade design results.
r (mm)λabFφ (deg)θ′ (deg)Re/ × 104ci (mm)ci+1 (mm)Δ (mm)c (mm)θ (deg)250.220.301.870.9049.9143.915.6272.9472.940.000695.4324.02500.430.320.770.9142.7336.739.09109.56109.560.000192.4822.52750.650.330.430.9236.6830.6811.33122.61122.611 × 10−689.5321.021000.870.330.270.9331.7225.7212.79122.71122.716 × 10−786.5819.521251.090.330.190.9427.6921.6913.74116.82116.822 × 10−583.6318.021501.300.340.140.9524.4318.4314.38108.74108.74−1 × 10−580.6816.521751.520.340.100.9621.7715.7714.82100.26100.264 × 10−577.7315.022001.740.340.080.9619.5813.5815.1392.1892.183 × 10−674.7813.522251.960.340.060.9617.7511.7515.3684.8084.809 × 10−671.8311.652502.170.340.050.9616.2010.2015.5378.1978.194 × 10−668.8810.112752.390.340.040.9614.898.8915.6672.3072.302 × 10−665.938.833002.610.340.040.9613.757.7515.7667.0767.076 × 10−562.997.743252.830.340.030.9612.766.7615.8462.4262.42−3 × 10−560.046.813503.040.350.030.9611.895.8915.9058.2758.27−8 × 10−757.095.993753.260.350.020.9511.125.1215.9454.5454.54−9 × 10−754.145.254003.480.350.020.9410.424.4215.9751.1751.17−1 × 10−651.194.554253.700.350.020.939.783.7815.9848.1048.10−1 × 10−648.243.884503.910.360.020.919.183.1815.9745.2745.27−9 × 10−745.293.224754.130.370.020.888.602.6015.9142.5842.58−8 × 10−742.342.575004.350.380.020.838.022.0215.7839.9239.92−8 × 10−539.391.935254.570.410.020.757.371.3715.4737.0137.01−2 × 10−736.441.315504.780.450.020.606.500.5014.5832.8232.822 × 10−633.490.725755—————————30.540.20Figure 15
Chord length and pitch angle before and after modification: (a) chord length; (b) pitch angle.
(a)(b)The point coordinates required to visualize a 3D model were obtained, and the frameworks were plotted as shown in Figure16(a) and modeled as shown in Figure 16(b). The whole rotor is shown in Figure 16(c).Figure 16
The new wind turbine rotor: (a) cross-sectional profiles of all sections; (b) one blade; (c) whole turbine rotor.
(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 Figure17(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−4.Figure 17
Calculation zones and mesh: (a) calculation zones; (b) mesh for rotation domain; (c) mesh for stationary domain; (d) mesh for surface of rotor.
(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 Figures17(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. Ify + 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.Table 8
The calculation results with different number of mesh.
ItemElementsCalculation time (h)Torque coefficientCMAxial thrust coefficientCTMesh 13564665270.0298960.47433Mesh 26427530510.0318770.48897Mesh 37592941760.0316600.48717The 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.Figure 18
Time histories of dynamic pressure on blade surface for different time step sizes.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 Figure19, 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.Figure 19
Time histories of the torque coefficient.


### 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 Figure20, 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.Figure 20
Time histories of the power coefficient of different rotors.The axial thrust coefficients values of different rotors in relation to time, shown in Figure21, 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.Figure 21
Time histories of the axial thrust coefficient of different rotors.The pressure distribution on the surface of the two rotors att = 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.Figure 22
Pressure distribution on the surface of the rotors: (a) pressure surface; (b) suction surface.
(a)(b)The turbulent kinematic energy distribution on the surface of the two rotors att = 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 m2/s2 and 0.000206362 to 90.7744 m2/s2, 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.Figure 23
Turbulent kinematic energy distribution on the surface of rotors: (a) pressure surface; (b) suction surface.
(a)(b)The pathlines colored by the turbulent kinematic energy around the two rotors, shown in Figure24, 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.Figure 24
Pathlines around the rotors. (a) The original turbine rotor. (b) The new turbine rotor.
(a)(b)


## 4.1. Calculation Zones and Boundary Conditions

The rotation domain near the rotor and stationary domain showing the entire calculation zone is illustrated in Figure17(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−4.Figure 17
Calculation zones and mesh: (a) calculation zones; (b) mesh for rotation domain; (c) mesh for stationary domain; (d) mesh for surface of rotor.
(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 Figures17(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. Ify + 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.Table 8
The calculation results with different number of mesh.
ItemElementsCalculation time (h)Torque coefficientCMAxial thrust coefficientCTMesh 13564665270.0298960.47433Mesh 26427530510.0318770.48897Mesh 37592941760.0316600.48717The 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.Figure 18
Time histories of dynamic pressure on blade surface for different time step sizes.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 Figure19, 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.Figure 19
Time histories of the torque coefficient.


## 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 Figure20, 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.Figure 20
Time histories of the power coefficient of different rotors.The axial thrust coefficients values of different rotors in relation to time, shown in Figure21, 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.Figure 21
Time histories of the axial thrust coefficient of different rotors.The pressure distribution on the surface of the two rotors att = 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.Figure 22
Pressure distribution on the surface of the rotors: (a) pressure surface; (b) suction surface.
(a)(b)The turbulent kinematic energy distribution on the surface of the two rotors att = 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 m2/s2 and 0.000206362 to 90.7744 m2/s2, 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.Figure 23
Turbulent kinematic energy distribution on the surface of rotors: (a) pressure surface; (b) suction surface.
(a)(b)The pathlines colored by the turbulent kinematic energy around the two rotors, shown in Figure24, 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.Figure 24
Pathlines around the rotors. (a) The original turbine rotor. (b) The new turbine rotor.
(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.



---

*Source: 3947164-2022-05-09.xml*