Abstract
This work presents a novel framework for the aerodynamic design and optimization of blades for small horizontal axis wind turbines (WT). The framework is based on a state-of-the-art blade element momentum model, which is complemented with the XFOIL 6.96 software in order to provide an estimate of the sectional blade aerodynamics. The framework considers an innovative nested-hybrid solution procedure based on two metaheuristics, the virtual gene genetic algorithm and the simulated annealing algorithm, to provide a near-optimal solution to the problem. The objective of the study is to maximize the aerodynamic efficiency of small WT (SWT) rotors for a wide range of operational conditions. The design variables are (1) the airfoil shape at the different blade span positions and the radial variation of the geometrical variables of (2) chord length, (3) twist angle, and (4) thickness along the blade span. A wind tunnel validation study of optimized rotors based on the NACA 4-digit airfoil series is presented. Based on the experimental data, improvements in terms of the aerodynamic efficiency, the cut-in wind speed, and the amount of material used during the manufacturing process were achieved. Recommendations for the aerodynamic design of SWT rotors are provided based on field experience.
1. Introduction
There has been an increased necessity for the development of small grid-independent energy applications [1] (e.g., domestic, street lighting, and rural electrification) in which wind energy plays an important role. However, the integration of small wind energy technologies is challenging, mostly because of the low average wind speed experienced at the locations where small wind turbines (SWT) are commonly installed [2]. In addition, SWT show reduced aerodynamic efficiencies as compared to large-scale wind turbines (WT) [3] given that they generally operate at low Reynolds numbers and under higher turbulent conditions. SWT designers and manufacturers have shown considerable interest in improving the overall efficiency of their prototypes in order to increase their competitiveness against other small-scale renewable energy technologies. Hence, in order to address this need, improved SWT aerodynamic, acoustic, structural, mechanical, and electrical technologies are required.
Several contributions of variable complexity have been reported in the literature with regard to the design and optimization of WT blades [3–13]. Even though the problem has been extensively studied, no general theories [10] or general conclusions for the optimal design of WT blades have been proposed. In addition, only few works have made major contributions to the technoeconomical modeling and optimization of WT blades while considering all-encompassing schemes [5–8, 10]. The classical blade design methodology adopted in such works consists of three steps: (1) the definition of the blade performance metric (e.g., the aerodynamic efficiency, the cost of energy, etc.) and the objective function to be extremized, (2) the definition of the decision variables and their respective bounds, and (3) the choice of the optimization strategy. Most works determined the blade’s aerodynamic efficiency through the use of a blade element momentum (BEM) model, considered complex and highly constrained decision variables, and considered a population-based metaheuristic (e.g., the simple genetic algorithm (GA)) or a calculus-based procedure (e.g., sequential quadratic programming) to provide a solution to their optimization problem. The common pitfall of such works is that the adopted solution procedures (e.g., calculus-based optimization) are likely to reach local maxima due to the complexity of their objective functions. In this regard, the use of state-of-the-art mathematical tools of soft computing may provide a better solution to this rather complex problem. Moreover, the majority of the works have evaluated only a few scenarios and have failed in providing a comprehensive description of their optimized results and, therefore, do not provide appropriate guidelines and/or recommendations to WT blade designers.
This research contributes to the objective of improving SWT aerodynamic technologies by developing a pragmatic framework for the aerodynamic design and optimization of SWT blades and by conducting a validation study, which in turn allowed identifying key features during the SWT blade design. The proposed framework can be considered as a constituent part of a more general, holistic, and multidisciplinary framework known as the wind turbine blade optimization (WTBO) problem. The WTBO problem consists of the design and manufacture of WT blades that exhibit high aerodynamic efficiency, low noise generation, and proper structural design and manufacture to withstand the dynamic loads generated by the expected wind field at the minimum cost and while considering a control/regulation system [11] (e.g., active stall or furling systems [14], among others). The WTBO problem has been applied in the design of blades for medium and large-scale pitch-regulated WT; however, it has been scarcely applied in the design of blades for small-scale stall-regulated or furling-regulated WT [3, 4, 10], hence denoting an important subfield of study that requires more attention.
The starting point in the design and optimization of SWT rotors are the existing guidelines and international standards which give a clear overview of the important factors and safety regulations that affect the decision making of the WT designer. For SWT, requirements may be altered if the SWT safety is not compromised as stated by the International Electrotechnical Commission (IEC) standards [15], which is often the case given the highly constrained scenarios where SWT operate. In this regard, SWT must be as flexible as possible in order meet the energy demand; therefore, any other relevant secondary objectives, such as noise generation, may be relaxed in order to maximize the cost-effectiveness of SWT. This work is focused on the study of the aerodynamic behavior of optimized SWT rotors. Secondary objectives related to the minimization of acoustic, structural, mechanical, and electrical issues were not explicitly considered in the proposed framework with the aim of studying the key geometrical factors influencing the aerodynamic efficiency of SWT rotors.
The remaining of this work is structured as follows. Section 2 presents the theoretical foundations and the development of the SWT blade optimization framework. Section 3 describes the development of an optimized SWT prototype which is used to perform a validation study of a set of optimized rotors. Section 4 discusses the results. Finally, Section 5 presents the overall conclusions and outlines future research.
2. Theoretical Development
2.1. Two-Dimensional Steady-State Airfoil Aerodynamics
A large range of studies have focused on determining the experimental aerodynamics of two-dimensional airfoils [16–20], which are the base for the design of WT rotors. However, the characterizations are commonly performed for a limited set of angles-of-attack (α) and medium-to-large Reynolds numbers (Re). In addition, it is unusual to find aerodynamic characterizations of complete series of airfoils. To overcome these issues, recent approaches in the literature have focused on the development of automatized procedures for the creation of new airfoils and its aerodynamics are studied theoretically through the implementation of computational fluid dynamics (CFD) approaches [19, 21–23]. Despite the advances in airfoil shape development and airfoil aerodynamics prediction, only a few studies have addressed which types of airfoils are appropriate for the construction of efficient SWT rotors [4, 18, 24, 25]. The result thereof is the fact that SWT aerodynamic technologies have been scarcely explored; thus, a formal study to identify the geometric blade features that improve the SWT aerodynamic efficiency is of practical and theoretical relevance. In order to meet this objective, the present work considers the complete National Advisory Committee for Aeronautics (NACA) 4-digit airfoil series as a test case.
There is an infinite set of two-dimensional airfoils that might be used in the development of WT blades. The NACA 4-digit series of airfoils [16] are a small subset of analytically developed airfoils that have been used in different applications, including wind energy [26]. Each of the four digits describes different geometrical features in the set. The first digit describes the maximum camber as a percentage of the airfoil chord length. The airfoil camber is defined as a measure of the asymmetry between the top and the bottom surfaces of the airfoil. Therefore, airfoils without camber are defined as symmetric airfoils. The second digit describes the distance of the maximum camber from the airfoil leading edge in percentage tens of the chord length. The last two digits describe the maximum thickness of the airfoil as a percentage of the airfoil chord length. For example, the NACA 7512 is an airfoil having a maximum camber of 7% located at 50% from the leading edge and having a maximum thickness of 12% with respect to the airfoil chord length.
Only a small set of the NACA 4-digit airfoils have been experimentally characterized [16, 17, 20]. Therefore, their aerodynamics is approximated through the implementation of CFD approaches. In the present work, the two-dimensional, incompressible, steady-state aerodynamic forces, in terms of the aerodynamic coefficients of lift and drag , are obtained by means of the XFOIL 6.96 software. XFOIL [27] is a Fortran-based software, created by Drela in 1989, for the analysis of the subsonic aerodynamics of isolated airfoils. XFOIL computations are based on a combined panel method and an integral boundary layer formulation, which in turn is complemented with an laminar-to-turbulent transition method. Given an initial inflow condition, the flow velocity distribution around the airfoil is computed from the panel method while accounting for viscous forces and the induced vorticity from the airfoil surface. The resultant boundary layer and wake are interacted with a surface transpiration model. The resultant flow field is incorporated into the fluid mechanics viscous equations, yielding a nonlinear elliptic system of equations which is solved by a Newton-Raphson algorithm, resulting in both a complete pressure and velocity distributions in the airfoil vicinity. The lift coefficient is calculated by direct surface pressure integration, as viscous contributions to the lift force can be neglected, and the pressure coefficient is calculated using the Karman-Tsien compressibility correction. The drag coefficient is determined from the wake momentum thickness far downstream and calculated with use of the Squire-Young formulation. The methods, corrections, and the boundary layer formulation used in XFOIL are extensively described in [27, 28]. The XFOIL 6.96 software was incorporated into the BEM model (described in next section) and the hybrid solution procedure, both of which were developed in a MATLAB 2012a environment.
The airfoil aerodynamics can be computed only for pre-stall conditions (e.g., at low angles-of-attack) since the complex flow behavior of the stall state cannot be accurately modeled by XFOIL. Stall is a local boundary layer flow detachment from the blade surface experienced when the local effective angle-of-attack is large and the flow stream lines cannot follow the blade surface due to the large induced adverse pressure gradients. The stall condition produces dramatic reductions in the local lift force while increasing the local drag force, thus, reducing the aerodynamic efficiency of the WT rotor. Nevertheless, induced stall conditions are often useful for control and regulation purposes which, however, are not considered in this work. A complete database of aerodynamic coefficients was constructed from XFOIL 6.96 for the range of airfoil thicknesses [%] (i.e., the last two digits of the NACA 4-digit series), the range of angles-of-attack , and the range of Reynolds numbers for each combination of the maximum camber value and the maximum camber position (i.e., the first two digits of the NACA 4-digit series). The database quality was inspected and a few corrections (e.g., related to numerical divergences) were performed at intermediary angles-of-attack and Reynolds numbers by implementing bicubic spline interpolations. The resultant aerodynamic coefficients were extrapolated up to by means of the Viterna-Corrigan semiempirical method [29, 30] in order to model the stall state. The database was stored and specific aerodynamic values were made available to the BEM model via three-dimensional bicubic spline interpolations during the rotor aerodynamic performance evaluation.
2.2. Blade Element Momentum (BEM) Model
The WT aerodynamic efficiency is quantified through the use of a state-of-the-art blade element momentum (BEM) model [29–31], which has shown good accuracy when predicting the aerodynamic efficiency of large-scale WT blades. The main differences between the BEM models proposed in the literature are related to its numerical solution procedures and the semiempirical correction factors used to account for the effects that were not explicitly considered during the derivation of the model. The steady-state BEM model relates the blade geometry with the rotor aerodynamic performance depending on the magnitude of the free-stream wind speed approaching the WT and the WT rotational speed ( [RPM]).
The BEM model is composed of two different parts. The first part is focused on the description of the local aerodynamics of the different blade elements (e.g., blade sectional airfoils), which are subject to a force balance analysis that is dependent on the local inflow condition. The second part is focused on the description of the global processes of mass and momentum conservation in which the integral versions of the governing equations of the fluid mechanics are used to quantify and balance the changes in linear and angular momentum of the air flow with the forces and torques on the WT rotor. The resultant set of coupled equations are reinterpreted using the concepts of axial and tangential induction factors, which describe the effective local axial and tangential wind speed magnitudes in a normalized fashion. The Prandtl factor , as shown in (2), is included in the BEM nonlinear model to describe additional losses caused by the generation of vortices at the tip of the blades.
The BEM model variables are the axial and tangential induction factors, which are readily solved using a Newton-Raphson algorithm for a given inflow condition and certain blade geometric characteristics while considering a predefined maximum number of iterations and a finite numerical tolerance. The blade geometric characteristics are included in the BEM model both explicitly, through the calculation of the chord radial solidity , and implicitly, during the calculation of the effective lift and drag aerodynamic coefficients at each of the blade span positions. The Reynolds number effects are incorporated into the BEM model during the calculation of the aerodynamic coefficients. The effective Reynolds number at each blade span position is quantified using (3).
BEM model:
Prandtl factor:
Effective Reynolds number:
Power coefficient:
Thrust coefficient:
Total torque at rotor and torque coefficient:
From the solved induction factors, normalized variables describing the WT rotor aerodynamic performance can be derived, such as the power coefficient , the thrust coefficient , and the torque coefficient , which in turn can be transformed into the dimensioned variables of mechanical power , axial force , and torque force , respectively, as shown in (4)–(6). The BEM model assumptions, limitations, and complete derivation can be consulted in [29, 30] and its nomenclature is presented in Table 1.
2.3. Problem Definition and Optimization Strategy
The blade variables to be optimized are the blade sectional airfoil shapes and the geometrical distributions of chord length ( [m]), twist angle ( [°]), and thickness ( [%]) along the WT blade span. However, there is no general theory (i.e., there are no analytical closed-form expressions or concise guidelines that are not based on crude assumptions) describing the optimal radial variation of these variables. Therefore, robust and flexible analytical functions that depend on a reduced set of continuous parameters were developed to describe such geometric distributions, as shown in (8)–(10). Each of the different analytical functions describing the geometrical distributions contains four variables. The variables define the initial and final value (i.e., the value at and ) of the geometrical distributions, respectively. The variables define the concavity and vertical displacement of the geometric distributions, respectively. The blade radial discretization was performed in a nonuniform way, as described in (7a) and (7b), in order to reduce the amount of discrete blade segments, which proved to reduce the demand of computational resources during the optimization process without incurring numerical inaccuracies during BEM calculations.
After performing a sensitivity analysis, it was concluded that 17 discrete radial positions were enough to obtain accurate numerical values of the BEM model predictions. This value differs from the values considered in previous studies [7, 13] and special attention is required since considering less than 17 discrete positions might lead to numerical inaccuracies that in turn might influence the optimization results. The parametric expressions used for the representation of the normalized radial coordinate, the local chord, twist angle, and thickness are the following:
The explicit variables to be optimized are the continuous values of . Each of the 12 variables to be optimized must be linearly bounded to avoid unphysical solutions. The WT rotor aerodynamic efficiency, in terms of the BEM-derived power coefficient, is the performance metric to be maximized in the present work. However, in order to account for different operational states, the mono-objective performance metric to be maximized is the integral of the power coefficient over a range of tip speed ratios, as defined in (11). The selection of this performance metric has the ultimate effect of improving the annual energy production of the SWT even without explicitly considering a site-specific wind resource.
Additional metrics are monitored during the optimization process, including the integral of the thrust coefficient over the same range of tip speed ratios, as defined in (12), the cut-in wind speed, and the total volume of the blade. The cut-in wind speed is approximated from the simple balance between (1) the product of the free-stream dynamic pressure and the blade solidity, which results in an aerodynamic force over the rotor that is dependent on the blade geometry, and (2) the torque needed to overcome the static friction and/or the cogging torque of the electric generator , as shown in (13). Consider
The complete WTBO model, as presented, is a mixed-integer, linearly constrained, nonconvex, and nonlinear optimization model. The number of elementary operations needed to compute the aerodynamic efficiency of a WT rotor increases in an approximately linear way as a function of the number of blade discrete segments . However, the number of possible combinations of blade sectional airfoils increases in an exponential way as a function of the number of discrete blade segments and the number of available airfoils . Therefore, the WTBO model can be categorized as an NPO-complete problem [32], one of the hardest problems to solve in computational optimization.
Calculus-based optimization approaches such as linear/quadratic programming, gradient methods, and Newton-Raphson, among others cannot be used to solve the problem to optimality. In addition, it is considered inadequate to implement relaxation procedures (i.e., procedures that transform the original problem into another equivalent problem of reduced complexity) over the BEM model in order to use tree search-based procedures (e.g., branch and bound, branch and cut, etc.) since the already simplified physics will be greatly influenced, hence, significantly affecting the original problem and its optimal solution. To attempt solving optimization problems of such complexity, heuristic and metaheuristic approaches [33] are commonly employed. However, heuristic and metaheuristic algorithms cannot guarantee optimality. Therefore, multiple executions of the optimization procedure must be performed to obtain a near-optimal solution to the problem instance being solved, as shown in Section 4. A wide range of metaheuristics can be used to perform the optimization task. However, evolutionary-based algorithms have proven to be efficient while optimizing WT blades [5, 10, 13].
This work extends the optimization effectiveness of previous approaches by developing a novel hybrid optimization procedure in which the geometrical shape of the blade is optimized in a nested procedure, while including more decision variables (e.g., number of airfoils and number of discrete segments). The proposed hybrid procedure consists of two parts: the first part focuses on optimizing the geometrical distributions of chord length, twist angle, and thickness, given a combination of sectional airfoils, by employing a local search metaheuristic based on the simulated annealing (SA) algorithm. The combination of sectional airfoils is proposed by a leading population-based metaheuristic based on the virtual gene genetic algorithm (vgGA). Hence, the SA algorithm is nested into the vgGA during the iterative optimization process. In the next subsections, a detailed description of the metaheuristic algorithms and the hybrid procedure is presented.
2.4. Local Search Metaheuristic Based on the Simulated Annealing Algorithm
The present work implements the SA algorithm to optimize the blade geometric distributions, shown in (8)–(10), in order to maximize the objective function, shown in (11). The SA algorithm [33–35] is based on the natural process of annealing of solids, which is a thermal process for obtaining low energy states in solids by exposing them to a heat bath. The annealing process consists of the following two steps: (1) increase the temperature of the heat bath until the solid melts and (2) slowly decrease the temperature of the heat bath until the particles arrange themselves in the low energy ground state of the solid. In the melt phase, the constituent particles arrange themselves randomly, and in the cooling stage the particles arrange themselves in an organized lattice. The analogy between the annealing process and the optimization process is that potential solutions of the studied problem are comparable to states of a physical system undergoing a process of annealing and the quality of a potential solution is equivalent to the energy of a state.
While traditional local search (and variable local search) procedures will generally accept transitions between feasible neighbor solutions only if its quality is better than the current solution, the SA algorithm has the ability of allowing transitions between feasible solutions whose qualities are worse than the current solution in the hope of avoiding getting stuck at a local extreme point. The transition to a worse solution is modeled as a probabilistic process, which relies on the Metropolis criterion [34] and depends on the actual state and the new state quality. The probability of rejection of states with lower quality increases as the search continues (i.e., the SA algorithm behaves as a greedy algorithm), therefore, simulating the arrangement of particles in an organized lattice during the cooling process.
The SA’s ease of implementation, convergence properties, and ability to escape local extreme points have made it a popular technique during the past two decades. Therefore, the SA algorithm is considered as a suitable solution procedure for solving the WTBO model presented in this work. The interested reader is referred to [33–35] for more detailed descriptions of the SA algorithm. The SA’s state is represented by a 4 × 3 matrix containing the continuous values of the geometrical variables . The initial state is computed randomly within the geometrical variables bounds. During the optimization process, neighbor solutions are generated according to where represents the actual state and is the continuous uniform probability density function that computes random neighbor matrices within the interval . After computing a neighboring solution, a feasibility check is performed to verify that the new state does not contain values outside the defined bounds. In the case of unfeasibility, the elements outside the bounds are set to the nearest boundary value. The SA implementation allows for the variation of the temperature within a Markov chain every time a transition is accepted, thereby allowing the algorithm to escape local maxima by relaxing the Metropolis criterion [33]. A schematic diagram of the SA-based procedure is shown in Figure 1. The SA’s parameters were tuned by conducting a statistical design of experiments, as described in Section 2.6.
2.5. Population Metaheuristic Based on the Virtual Gene Genetic Algorithm (vgGA) and the Hybrid vgGA-SA Solution Procedure
The present work implements the virtual gene genetic algorithm (vgGA) [36] to optimize the selection and combination of different airfoil shapes along the WT blade span. The GA refers to a class of adaptive search procedures based on the principles derived from natural evolution and genetics. In the GA, genes represent the design variables and potential solutions (i.e., individuals) are represented by chromosomes containing genes. Like other population-based methods, as the GA proceeds through generations (i.e., iterations), it holds a set of possible solutions referred to as the population. Based on the fitness value (i.e., the quality of a solution), which is determined through the evaluation of the population’s individuals, two individuals (i.e., the parents) are selected and are combined through genetic operators (crossover and mutation) resulting in two new individuals (e.g., the offspring). It is expected that the population fitness will improve through the generations (iterations) until a near-optimal solution is obtained.
Many different variants of GAs have been considered in the literature [33] differing in the way genetic operators act over the population and the way individuals are represented. Research has demonstrated that the performance of a GA relies upon the proper choice of the selection and crossover mechanisms and the proper design of the individual representation. In the present work, a binary tournament selection, a uniformly distributed one-point crossover, and a uniformly distributed mutation methods were considered. The individual’s representation is based on a chromosome containing binary digits (genotype). The combination of binary digits in turn defines the airfoil shape to be used at each of the blade span positions (phenotype). One of the issues with considering this kind of representation is related to the estimation of the minimum number of bits needed to completely describe the range of values each variable can assume. In order to quantify this value, the chromosome is divided into segments, corresponding to the number of blade discrete points. Each segment will contain MAB bits, quantified using (15), which are sufficient to represent the number of possible airfoils that each blade section can assume. Consider where the ceil function rounds the argument towards infinity and is the number of available airfoils. Each of the segments in the chromosome (genotype) is reinterpreted as a real number (phenotype), which is rounded to obtain a feasible integer value representing an airfoil identifier. For example, a hypothetical first chromosome segment [000⋯000] represents that the first available airfoil (i.e., airfoil #1) will be located at the first discrete segment of the blade, which is the closest to the blade root. Another hypothetical first chromosome segment [111⋯111] represents that the last available airfoil will be located at the first discrete segment of the blade, which is the closest to the blade root. A complete example showing the genotype and the phenotype representations of a hypothetical individual, considering 8 available airfoils and 4 discrete segments, is shown next:
Note that in the previous example the number of discrete segments was conveniently limited to . As discussed before, the minimum number of segments needed to avoid numerical inaccuracies during BEM calculations is approximately . In the remaining of the work, discrete segments were considered.
Once an entire population is constructed, its evaluation proceeds by executing five times the SA-based procedure proposed in the previous section, in which the near-optimal radial distributions of the geometrical variables of chord length, twist angle, and thickness are obtained for each individual. The quality of the best found solution, obtained through the five SA-based executions, is considered as the fitness of the individual. A schematic diagram of the vgGA-SA hybrid solution procedure is shown in Figure 2. The SA’s and the vgGA’s parameters were tuned by conducting two independent statistical designs of experiments, as presented in the next subsection.
2.6. The vgGA-SA’s Parameters Tuning through Statistical Designs of Experiments
Optimization procedures based on heuristic and metaheuristic algorithms rely on their parameters to obtain high-quality/near-optimal solutions to the problem being investigated. It has been observed that the use of poorly tuned parameters may limit the ability of optimization procedures to find high-quality solutions. The parameter tuning of complex algorithms is a well-known problem; however, there are no standard procedures to perform the tuning task.
In order to properly tune the vgGA-SA’s parameters, two independent full factorial designs of experiments (DOEs) were conducted. The outcomes of the DOEs allowed (1) understanding which factors (parameters) are significant during the optimization process, (2) finding the best combination of values for the parameters of both the SA and the vgGA algorithms, and (3) quantifying the effects of using poor-tuned parameters. Both DOEs were processed with the aid of the statistical software Minitab v16.
SA’s Parameters Tuning. For the SA’s parameters tuning process, the rotor parameters described in Section 3 and the NACA 75XX airfoil were considered. Table 2 summarizes the SA’s factors to be tuned and the considered levels. The DOE considered 10 replicates for each combination of parameters’ levels, resulting in 80,000 SA-based simulations.
The DOE evaluation involves performing an analysis of variance (ANOVA), the aim of which is to determine what factors are significant during the optimization process. The ANOVA model adequacy verification did not exhibit issues with regard to (1) the residual’s normality behavior, (2) the experiments independency, and (3) the constant variance assumptions. The validations of the first and third assumptions were performed by conducting an Anderson-Darling normal probability test and a Bartlett’s and Levene’s constant variance tests, respectively.
With a significance level of 5%, the ANOVA results suggested that there is relative significance between all the factors, which means that specific combinations of them will allow the SA-based procedure to reach high quality solutions. To enhance the understanding of how these factors influence the optimization results, Figure 3 presents the main effects of the factors on the response (i.e., the objective metric, quantified using (11)). From the main effects charts shown in Figure 3, it can be observed that some factors are more important than others on average (although all factors are statistically significant, as suggested by the ANOVA).
From Figure 3(a) it can be observed that by varying the system’s temperature within Markov chains, a slightly better solution quality can be obtained on average. The temperature variation along Markov chains allows the system’s temperature to increase when high quality states are explored. The result of increasing the system’s temperature is that the SA greedy behavior is relaxed in the cooling stage, and thus the algorithm is granted with the ability to continue searching for better states and avoids the possibility of being trapped at local maxima.
Figure 3(b) indicates that the initial temperature has a significant impact on the performance of the SA algorithm. The initial system’s temperature must be low in order to obtain a better solution quality on average. The SA algorithm tends to be more selective (i.e., behaves as a greedy algorithm) if the system’s temperature is low enough because the Metropolis criterion will reject transitions to low quality states more often.
From Figures 3(c), 3(d), and 3(e) it can be noted that by increasing the maximum number of accepted trials, the maximum number of trials, and the number of Markov chains, a better solution quality can be reached on average. This behavior is to be expected since the algorithm is allowed to continue its search for the best solution. Nevertheless, no significant improvements in solution quality are observed when allowing more than 50 trials.
Figure 3(f) shows that the neighborhood size, as a percentage of the search space size, has great influence on the performance of the SA algorithm. By allowing large state transitions, an improved solution quality can be obtained on average. This implicitly means that the solution space is not smooth, which is a common characteristic of nonconvex search spaces. To overcome the issue of optimizing variables on nonconvex spaces, the SA algorithm performs abrupt changes in the state in order to properly explore the solutions space at the cost of increasing the variance of the state quality during the optimization process.
From Figure 3(g), it can be observed that the cooling rate has little effect on the mean response. By lowering the magnitude of the cooling rate a slightly better solution quality can be obtained on average. As the system’s temperature decreases, the SA algorithm becomes more selective regarding the visited states; thus, low quality states are discarded by the Metropolis criterion.
After a careful inspection of the DOE results and after performing a detailed study of the best solution among the 80,000 simulations, the combination of parameters that maximizes the performance of the SA-based procedure was identified and is summarized next: neighborhood size , system’s initial temperature , system’s cooling rate , maximum number of accepted trials , maximum number of trials , and number of Markov chains . The matrices and refer to the upper and lower bounds of the 12 variables to be optimized and its difference defines the size of the search space.
vgGA’s Parameters Tuning. The vgGA’s parameters were tuned following the previously described procedure. Once again, the rotor parameters described in Section 3 were considered. In addition, the nested SA-based procedure considered the previously tuned parameters. Table 3 resumes the vgGA’s factors to be tuned and the considered levels. The DOE considered 10 replicates for each combination of parameters’ levels, resulting in 810 vgGA-SA-based simulations.
The ANOVA model adequacy verification did not exhibit issues with regard to (1) the residual’s normality behavior, (2) the experiments independency, and (3) the constant variance assumptions. With a significance level of 5%, the ANOVA results suggested that the significant factors are the population size, the number of generations, and the mutation probability. It is worth nothing that the crossover probability was not significant.
Figure 4 presents the main effects of the factors on the response. From Figures 4(a) and 4(b), it can be noted that by increasing the number of generations and the population size a much better solution quality can be reached on average. This behavior is to be expected since the algorithm is allowed to continue its search for the best solution. Although both figures suggest that additional improvements in the average solution quality may be obtained when increasing the value of both parameters’ levels (e.g., beyond the proposed levels), it is considered appropriate to maintain such levels since the trade-off between the computational resources spent and the obtained solution quality worsens as the levels’ values are increased.
Figure 4(c) indicates that the mutation probability has a significant impact on the performance of the vgGA algorithm. On the one hand, small values of the mutation probability inhibit the ability of the vgGA to explore the search space. On the other hand, large values of the mutation probability lead to large variance values on the population fitness, thus the vgGA algorithm cannot identify and combine the characteristics of the best individuals.
Although the crossover probability is not a significant factor (as suggested by the ANOVA), it is observed from Figure 4(d) that by increasing its value a slightly better solution quality can be obtained on average.
After a careful inspection of the DOE results and after performing a detailed study of the best solution among the 810 simulations, the combination of parameters that maximizes the performance of the vgGA-based procedure was identified and is summarized next: population size , number of generations , crossover probability , and mutation probability .
Both the SA’s and the vgGA’s best combinations of parameters values are highly dependent on the problem size. Therefore, a specific parameter tuning process must be performed for when considering different scenarios.
3. Experimental Development
A validation study of a set of optimized rotors based on the NACA 4-digit airfoils, as described in next subsections, was performed at the subsonic closed-return wind tunnel of the Center for Research and Innovation in Aerospace Engineering (CIIIA) [37]. The dimensions of the wind tunnel’s rectangular test section are 1 × 1 × 1.5 m3. A maximum wind speed of 60 m/s can be reached at the test section with the aid of a 300 kW frequency-controlled fan. All wind speed measurements were made using a hot-wire anemometer located at the center position in the upwind part of the test section.
In order to avoid flow-wall blockage effects, the rotor radius was set to 0.175 m during the rotor optimization process, thus limiting the swept rotor area to 10% of the wind tunnel cross sectional area. The radius of the blade root (i.e., the first blade position that contributes to the energy conversion process, ) was set to 0.02 m. The integration limits of the objective function, (11), were set to and with the purpose of (1) forcing the vgGA-SA hybrid procedure to find efficient rotors that have lower cut-in wind speeds and (2) reducing the magnitude of the rotational speed needed to operate at the maximum efficient state, therefore, improving the WT operational safety and reducing the magnitude of the rotational speeds needed during the aerodynamic characterization procedure. The following lower and upper bound geometric matrices were used during the optimization process:
The air properties (e.g., density and kinematic viscosity) were calculated as a function of the local pressure, temperature, and relative humidity as suggested by Picard et al. [38].
3.1. Small Wind Turbine Prototype
To perform the validation study of the optimized rotors, a SWT prototype was constructed. The constituent parts of the SWT prototype are (1) the optimized NACA 4-digit-based rotors, (2) a 14-pole neodymium permanent magnet, brushless, radial flux, synchronous, three phase EMP6354-KV200 electric generator, (3) an electrical system composed of a six-phase Guerte KBPC3510 rectifier, two voltage, and two current sensors (Allegro systems ACS758 and ACS715, resp.) located in both the AC and DC parts of the system, and (4) a set of electric resistances based on Kanthal A1 wire. The voltage and current measurements were monitored using a National Instruments data acquisition system (NI DAQ-USB-6211) and a LabVIEW interface, which reported and logged the real-time operation of the SWT prototype.
The optimized rotors were manufactured using a rapid prototyping Dimension Elite 3D printer. A high density ABSplus P430 thermoplastic material, printed in layers of 0.178 mm in thickness, was used during the rotor manufacturing process. After the manufacture, a rotor balance postprocess was performed, which consisted in (1) performing a slight sanding of the blade surface to reduce roughness and (2) adding thin layers of an epoxy material in an homogeneous way along the blade surface in order to balance the weight of the rotor and to improve its structural strength.
The EMP6354-KV200 electric generator was characterized following the theoretical and experimental procedures proposed in [29, 39, 40], where a set of characteristic curves of mechanical torque, rectified voltage, and rectified current as a function of both the rotational speed (η [rpm]) and the effective electric load were obtained by means of a customized electromechanical test bench. The electromechanical characterization allowed determining the aerodynamic efficiency of the SWT at any operational condition by comparing the available free-stream power and the mechanical power delivered to the electric generator. The SWT rotational speed ( [RPM]) was controlled by changing the effective electric load and calculated by measuring the AC voltage frequency , as shown in (18). In addition, the generator’s no-load opposing torque value was experimentally determined to be 0.0437766 Nm. Consider
The aerodynamic characterization was performed only in steady states, which were reached after a few seconds of operation. All the relevant variables were logged at a frequency of 14 Hz over a two-minute period for each operational state.
4. Results and Discussion
Two different types of evaluation procedures were conducted in order to demonstrate the capabilities of the proposed framework. The first procedure focused on evaluating the ability of the SA-based solution procedure for finding near-optimal radial geometrical distributions for blades composed of one type of airfoil. The second procedure focused on studying and comparing the best result obtained by the vgGA-SA hybrid solution procedure with the best results obtained by the first evaluation procedure. Both procedures considered the evaluation of the complete NACA 4-digit airfoil series (i.e., varying the maximum camber from 0 to 9% and the position of the maximum camber from 0 to 90% of the chord length, resulting in 100 different base airfoils).
First Evaluation. For the first evaluation procedure, the SA-based algorithm was executed thirty times, considering that the resultant blade was composed of only one type of NACA 4-digit airfoil. The average performance of the SA-based algorithm, considering the specific case of optimizing the geometrical distributions of a rotor based on the NACA 75XX airfoil, is shown in Figure 5. The SA-based algorithm was able to reach a high-quality solution in the first 80 iterations, while additional iterations were required to obtain a high quality near-optimal solution. This is consistent with the findings reported in Section 2.6. Similar results were obtained for the remaining 99 NACA airfoils. The average time required to perform an evaluation of the objective function was 0.1229 seconds, considering that the computational experiments were performed on a Sager NP8230-P151SM1 workstation with an Intel Core i7-4800MQ processor operating at 4 GHz, 16 GB of memory DDR3 @1600 MHz on an Intel HM87 chipset motherboard.
The best result among the thirty SA-based executions, for each of the 100 NACA 4-digit airfoils, was stored and is reported in Figure 6. From Figure 6 it can be noted that the aerodynamic efficiency of the optimized rotors is significantly affected by the selected base airfoil shape. Rotors based on commonly studied NACA 4-digit airfoils (e.g., NACA 44XX) are not the best performing rotors. It can be observed that by increasing the maximum airfoil camber value better aerodynamic efficiencies can be achieved. However, the position of the airfoil maximum camber also has great influence in the rotor efficiency. On the one hand, poor rotor efficiencies were obtained when the position of the maximum camber was located close to the leading edge or close to the trailing edge of the base airfoil. The worst efficiencies were obtained from the NACA 0XXX and the NACA X0XX airfoils, which are symmetrical airfoils. In addition, low efficiencies were obtained from the NACA [7-9][1-2]XX airfoils with the notable example of the NACA 92XX airfoils (note that for visualization purposes the NACA X1XX-based rotors performance metrics were not plotted). This is due to the high asymmetric features they exhibit close to the leading edge which promotes an early flow detachment from the blade surface (stall condition) that prevails over a large range of operational conditions. On the other hand, the best efficiencies were obtained when the position of the maximum camber was in the range between 40% and 60% of the chord length (e.g., NACA [7-9][4-6]XX airfoils).
The thrust force increases with increasing maximum camber value and when the maximum camber is located around 40% to 60% of the chord length. This is consistent with the fact that the rotor efficiency is high when such geometrical features are considered. In addition, it was observed that for most symmetric airfoils the thrust coefficient was slightly lower in contrast with asymmetric airfoils. This is consistent with the fact that lift and drag forces are generally larger for asymmetric airfoils, which results in larger axial forces on the rotor.
The cut-in wind speed decreases with increasing maximum camber. However, since the simple cut-in wind speed model is fully dependent on the geometric distributions (e.g., (8), (9), and (13)), it is not possible to identify explicitly if the position of the maximum camber influences the cut-in wind speed value. However, since the torque coefficient is directly related to the power coefficient, as shown in (6), the position of the maximum camber positively affects the cut-in wind speed when its values are in the range between 40% and 60% of the chord length, as described later in the experimental validation. The obtained results suggest that improvements in both the power coefficient and the SWT rotor cut-in wind speed can be achieved by properly selecting the base airfoil shapes of the SWT blade. This finding contradicts previous work [10], in which it is stated that it is necessary to trade some power-producing capability in order to obtain improvements in starting performance. It is conspicuous that for asymmetric airfoils (e.g., NACA [7-9]5XX) the SA algorithm identifies as best solution WT blades with high solidity, as noted during the calculation of the cut-in wind speed. This behavior is explained since high solidity rotors typically operate at larger effective Reynolds numbers, thus improving the aerodynamic efficiency of the rotor.
Finally, the blade volume decreases with increasing maximum camber. This result can be understood since blades composed of asymmetric airfoils (e.g., NACA 75XX) have lower cross sectional areas. However, the total volume was highly influenced by the resultant airfoil thickness distribution, which resulted in low values for blades having asymmetric airfoils (e.g., the NACA 75XX optimized rotor has a thickness distribution in the range [%]). This observed geometrical trend can limit the structural integrity of the WT rotor. Therefore, detailed structural studies must be performed for larger rotors experiencing larger axial forces.
Second Evaluation. For the second evaluation procedure, the hybrid vgGA-SA solution procedure was executed thirty times. The best result among the thirty executions is presented in Figure 7, in which a hypothetical fixed rotational speed of 2,500 RPM was assumed in order to construct the versus , the versus , and the versus performance curves. The hybrid procedure selected and combined different NACA 4-digit airfoils along the WT blade span. For this test case, a total of 10017 possible combinations of airfoils along the blade span can be constructed.
Interestingly, the vgGA-SA procedure considered the NACA 6[4-6]XX, the NACA 7[4-6]XX, the NACA 85XX, and the NACA 95XX airfoils in the best found solution, which were also the best performers in the first evaluation procedure. In addition, the vgGA-SA procedure determined nearly constant distributions of the chord length and thickness distributions along the blade span. The obtained high chord length value is attributed to the aerodynamic benefits of operating at larger Reynolds numbers while the low thickness values help reducing the side effects of operating in near stall conditions.
Figure 8 compares the geometry and the aerodynamic performance curves for different optimized NACA 4-Digit rotors. Based on the results shown in Figures 6 and 8, three optimized rotors based on the NACA hybrid, the NACA 44XX, and the NACA 75XX airfoils were manufactured, as described in Section 3, and tested at the CIIIA wind tunnel. Figures 9, 10, and 11 summarize the theoretical and experimental evaluation of the optimized rotors for different inflow conditions and electrical loads .
From Figures 9–11 several conclusions can be highlighted. The BEM model predicts the rotor aerodynamics reasonably well. However, a slight systematic overprediction of the power coefficient was found at most of the operational conditions. The overprediction was high at low Reynolds numbers and decreased monotonically with increasing Reynolds number. The observed overprediction can be partially attributed to an inadequate estimation of the aerodynamic coefficients by XFOIL 6.96 at low Reynolds numbers [22]. Depending on the operational condition (i.e., the tip speed ratio and electrical load) the efficiency of the WT can fall up to 10% due to Reynolds number effects. The maximum observed power coefficient was 45%, 46.7%, and 47% for the NACA 44XX, the NACA 75XX, and the NACA hybrid rotors, respectively. The theoretical cut-in wind speeds were 6.988 m/s, 6.018 m/s, and 5.848 m/s while the experimental cut-in wind speeds were 8.60 m/s, 6.26 m/s, and 6.03 m/s for the NACA 44XX, the NACA 75XX, and the NACA hybrid rotors, respectively. Additional improvements in the SWT starting performance can be obtained in practice by reducing the generator’s shaft static friction, the generator’s cogging torque, or by increasing the WT rotor diameter.
The rotor aerodynamics was predicted reasonably well in stall conditions. A metastable zone, related to aerodynamic instabilities, was observed during the transition from stall to nonstall operational states. Depending on the previous operational state of the WT, different operational paths can be observed upon changing the inflow speed condition. The size of the metastable zones increases when increasing the electric load as observed in Figures 7–9. Control systems relying on load-induced stall must account for the adverse effects of aerodynamic instabilities in order to (1) improve the SWT aerodynamic efficiency in low wind regimes and (2) properly regulate the WT under unwanted conditions.
Finally, the BEM model predictions are poor at large tip speed ratios. This is because the validity of the BEM model is limited to operational steady states where global flow detachment is not experienced. Global flow detachment occurs at high axial induction factors (e.g., when the WT rotational speed is large, so the rotor swept area is perceived as a solid disk by the free-stream wind) giving rise to vortex ring states. Empirical relationships are commonly used to model the conditions where the BEM theory breaks down [29, 30]. However, these generic relationships are not always valid and improvements to these expressions should be sought. Global flow detachment is also caused by high rotor solidity values. Therefore, SWT rotors must be designed to have moderately high rotor solidity values, since solidity increases the SWT efficiency by increasing the effective Reynolds number as previously discussed and must be designed to operate at lower rotational speeds. Lower rotational speeds can be achieved by considering hybrid rotors based on asymmetric airfoils (e.g., NACA [7-9]5XX) that perform well at low and medium tip speed ratios.
5. Conclusions and Future Research
This paper presented a novel and comprehensive framework for the aerodynamic design of small hybrid wind turbine rotors and proposed a hybrid solution procedure based on two metaheuristics (the genetic algorithm and the simulated annealing algorithm) to solve this complex problem.
The rotor performance was quantified through the implementation of a state-of-the-art blade element momentum model, which was coupled with the XFOIL 6.96 software and the Viterna-Corrigan extrapolation procedure for the determination of the local blade aerodynamics. It was observed that the steady-state blade element momentum model provides accurate predictions of the SWT rotor aerodynamic efficiency as long as stall conditions (observed at low tip speed ratios) or global flow detachment conditions (observed at large tip speed ratios) are not experienced and the supplied airfoil aerodynamics information (lift and drag coefficient as a function of both the Reynolds number and the angle-of-attack) is accurate.
The hybrid solution procedure optimized both the radial geometric distributions of the blade, by using a simulated annealing algorithm, and the selection of the airfoil shape at each of the blade discrete stations, by using the virtual gene genetic algorithm. The objective was to maximize the aerodynamic performance of the rotor for a wide range of tip speed ratios. Two different evaluation procedures were conducted to determine the ability of the proposed hybrid optimization strategy to find high-quality solutions, while considering the entire NACA 4-digit airfoil series as a test case. Several wind tunnel tests were performed to validate the achieved theoretical results. The validation procedure showed good results and demonstrated the suitability of the proposed framework for the development of SWT rotors. From the obtained results, several key geometric characteristics that influence the aerodynamic efficiency of SWT rotors were identified and studied. The findings are believed to be useful for SWT rotor development and increase the feasibility of using SWT in urban areas or at sites with a low wind resource. In addition, they allow WT designers to improve their prototypes with the aim of increasing their profitability prospections.
In addition to the specific objective of proposing a novel framework for the design of SWT rotors, this work studied the aerodynamic performance of the NACA 4-digit airfoil series both theoretically and experimentally with the aim of identifying the best performers for the construction of optimized SWT rotors. The theoretical and experimental results demonstrated that the NACA [7-9][4-6]XX airfoils perform very well for small wind energy applications. During the validation study, significant Reynolds number effects were observed. At low Reynolds numbers the power coefficient can be 5% up to 10% lower depending on the operation regime of the WT. Additionally, it was observed that certain asymmetric features in the airfoil shape geometry (e.g., large values of the maximum camber and the location of the maximum camber at a distance between 40% and 60% from the leading edge of the airfoil) increases the aerodynamic efficiency of the WT rotor. Three major improvements were achieved while using specific NACA 4-digit airfoils: (1) a reduction of the cut-in wind speed (up to 2.57 m/s in contrast with other NACA 4-digit airfoil-based rotors), (2) increased aerodynamic efficiency (the maximum achieved efficiency was 47% at tip speed ratios close to four), and (3) a reduction on the amount of material used in the manufacturing process, in contrast with other NACA 4-digit airfoil-based rotors.
Regarding future lines of research, the BEM model’s ability to predict the rotor aerodynamics at low and large tip speed ratios provides room for improvement and a fertile area of inquiry. In the stall regime, advanced computational fluid dynamics (CFD) tools that solve the Reynolds Averaged Navier-Stokes (RANS) equations complemented with turbulence models may provide better predictions of the airfoil aerodynamics, albeit at the cost of increased computational requirements. Improvements to the empirical relationships used to model the global flow detachment states are required in order to improve the predictions of the rotor aerodynamic efficiency at large tip speed ratios.
The optimized rotors did not incorporate technologies such as spoilers, winglets, stall barriers, vortex generators, or flaps. In case of incorporating these technologies, it might be possible to improve even more the SWT aerodynamic efficiency. Research on the performance improvements by the inclusion of these technologies is currently under way. Moreover, other families of airfoils are being tested as well.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
Support from the Research Chair for Wind Energy at Instituto Tecnológico y de Estudios Superiores de Monterrey (ITESM), under contract CAT158 is gratefully acknowledged. Special thanks are due to Dr. Patricia Zambrano and to all staff of the Center for Research and Innovation in Aerospace Engineering (CIIIA) of the Universidad Autónoma de Nuevo León for providing access to their wind tunnel facility and supporting the experimental section of this work.