Aerodynamic and Structural Integrated Optimization Design of Horizontal-Axis Wind Turbine Blades

Abstract

A procedure based on MATLAB combined with ANSYS is presented and utilized for the aerodynamic and structural integrated optimization design of Horizontal-Axis Wind Turbine (HAWT) blades. Three modules are used for this purpose: an aerodynamic analysis module using the Blade Element Momentum (BEM) theory, a structural analysis module employing the Finite Element Method (FEM) and a multi-objective optimization module utilizing the non-dominated sorting genetic algorithm. The former two provide a sufficiently accurate solution of the aerodynamic and structural performances of the blade; the latter handles the design variables of the optimization problem, namely, the main geometrical shape and structural parameters of the blade, and promotes function optimization. The scope of the procedure is to achieve the best trade-off performances between the maximum Annual Energy Production (AEP) and the minimum blade mass under various design requirements. To prove the efficiency and reliability of the procedure, a commercial 1.5 megawatt (MW) HAWT blade is used as a case study. Compared with the original scheme, the optimization results show great improvements for the overall performance of the blade.

1. Introduction

Blades are regarded as the key components of the Horizontal-Axis Wind Turbine (HAWT) system and have been paid much attention by most of the leading wind turbine manufacturers to develop their own blade design. As world wind energy market continuously grows, a number of blade manufacturers have emerged recently, especially in China. However, their independent design and manufacturing capabilities are weak, and the wind blade market is still dominated by the leading wind turbine system manufacturers [1]. Thus, to enter the market successfully and be more competitive, improving the fundamental technology on design and production of multi-megawatt (MW) class blades is indispensable.

The design process of the blades can be divided into two stages: the aerodynamic design and the structural design [2]. From the perspective of aerodynamic design, aerodynamic loads, power performance, aerodynamic efficiency, and Annual Energy Production (AEP) are important, and from the perspective of structural design, composite material lay-up, mass, stiffness, buckling stability, and fatigue loads are concerned [3]. A successful blade design should take into account the interaction between the two stages and satisfy a wide range of objectives, so the design process is a complex multi-objective optimization task characterized by numerous trade-off decisions. Nevertheless, in order to simplify the process, the aerodynamic design and the structural design are separated by the conventional methods. As of now, most of the research is focused mainly on the optimization of either aerodynamic or structural performances, which are unable to get the overall optimal solutions [4,5,6,7,8,9]. Only a limited number of works are concerned with the optimization of both aerodynamic and structural performances, and the relative procedure for this purpose is scarce.

Grujicic [10] developed a multi-disciplinary design-optimization procedure for the design of a 5 MW HAWT blade with respect to the attainment of a minimal Cost of Energy (COE), and several potential solutions for remedying the performance deficiencies of the blade were obtained. Bottasso [11] described procedures for the multi-disciplinary design optimization of wind turbines. The optimization was performed by a multi-stage process that first alternates between an aerodynamic shape optimization to maximize the AEP and a structural blade optimization to minimize the blade weight, and then combined the two to yield the final optimum solution. However, the problems only take a single objective function into account at each time, thus they can not obtain the trade-off solutions among conflicted objectives.

Benini [12] described a two-objective optimization method to design stall regulated HAWT blades, the aim was to achieve the trade-off solutions between AEP per square meter of wind park and COE. Wang [13,14] applied a novel multi-objective optimization algorithm for the design of wind turbine blades by employing the minimum blade mass and the maximum power coefficient, the maximum AEP and minimum blade mass as the optimization objectives, respectively. However, the blades are treated as beam models to calculate the structural performances in the above research, and the material layups were not considered.

This paper describes a procedure for the aerodynamic and structural integrated multi-optimization design of HAWT blades to maximize the AEP and minimize the blade mass. The scope is to find the balance between the design process to obtain the optimum overall performance of the blades, by varying the main aerodynamic parameters (chord, twist, span-wise locations of airfoils, and the rotational speed) as well as structural parameters (material layup in the spar cap and the position of the shear webs) under various design requirements.

2. Modeling of the Blade

2.1. Finite Element Method (FEM) Model of the Blade

As a starting point for the optimization, an initial definition of the blade aerodynamic and structural configuration is required. In this paper, a commercial 1.5 MW HAWT blade with a length of 37 m and a mass of 6580.4 kg is used for a case study.

Figure 1 shows the geometry shape of the blade, which can be divided into three areas: root, transition region and aerodynamic region. The root area is normally circular in cross section in order to match up with the bolted flange, the aerodynamic region uses aerofoil section to capture the wind, and the transition region from the root section to the aerofoil section should be a smooth one for structural reasons.

Figure 1. The geometry shape of the blade.

Figure 1. The geometry shape of the blade.

The geometrical shape of the blade can be described by the airfoil series and the chord, twist and percent thickness distributions. The main geometrical features are summarized in Table 1. Seven control points (CP) with fixed radial locations are used for the chord and twist distributions, as shown in Figure 2. The radial locations of the CP3-7 are defined using half-cosine spacing, CP1 is at root, CP3 is at the maximum chord station and CP7 is at the tip. Then, the chord and twist are defined with an 8th order Bezier curve and a 5th Bezier curve, except the twist remains constant inboard of the maximum chord. As Figure 3 illustrates, the percent thickness distribution is represented by the span-wise location of airfoils listed in Table 1. Six more control points with fixed percent thickness are defined. Then, a cubic polynomial is fit through the control points.

Table 1. The main geometrical features of the blade.

Location (m)	Airfoil	Chord (m)	Twist (°)	Percent Thickness (%)
0	Circle	1.88	10.00	100
1.0	Circle	1.88	10.00	100
6.8	DU400EU	3.02	10.00	40
9.3	DU300EU	2.98	7.30	30
13.7	DU91_W2_250	2.51	4.35	25
29.8	NACA_64_618	1.68	−0.33	18
36.5	NACA_64_618	1.21	−1.13	18

Table 1. The main geometrical features of the blade.

Location (m)	Airfoil	Chord (m)	Twist (°)	Percent Thickness (%)
0	Circle	1.88	10.00	100
1.0	Circle	1.88	10.00	100
6.8	DU400EU	3.02	10.00	40
9.3	DU300EU	2.98	7.30	30
13.7	DU91_W2_250	2.51	4.35	25
29.8	NACA_64_618	1.68	−0.33	18
36.5	NACA_64_618	1.21	−1.13	18

Figure 2. (a) Chord and (b) twist distribution of the blade.

Figure 2. (a) Chord and (b) twist distribution of the blade.

Figure 3. Percent thickness distribution of the blade.

Figure 3. Percent thickness distribution of the blade.

Figure 4 shows a typical structural cross section of the blade, which is composed of four parts: leading edge, trailing edge, spar cap and shear webs. The spar cap mainly consists of glass fiber composite materials, while the other three parts consist of glass fiber composite materials with Balsa and PVC core materials. The material properties are listed in Table 2.

Figure 4. A typical structural cross section of the blade.

Figure 4. A typical structural cross section of the blade.

Table 2. The material properties of the blade.

Material	E1 (GPa)	E2 (GPa)	G12 (GPa)	v12	ρ (kg/m3)	Thickness (mm)
UD-1200-xxx	42.2	12.5	3.5	0.24	1910	0.873
2AX-1200-0660	11.5	11.5	11.8	0.61	1909	0.878
3AX-1200-6330	26.9	13.4	7.5	0.47	1910	0.906
3AX9060-1200-0336	31.8	11.4	6.5	0.49	1910	0.906
Balsa	3.5	0.8	0.16	0.39	151	25.4
PVC	0.05	0.05	0.02	0.09	60	10/15/20/25

Table 2. The material properties of the blade.

Material	E1 (GPa)	E2 (GPa)	G12 (GPa)	v12	ρ (kg/m3)	Thickness (mm)
UD-1200-xxx	42.2	12.5	3.5	0.24	1910	0.873
2AX-1200-0660	11.5	11.5	11.8	0.61	1909	0.878
3AX-1200-6330	26.9	13.4	7.5	0.47	1910	0.906
3AX9060-1200-0336	31.8	11.4	6.5	0.49	1910	0.906
Balsa	3.5	0.8	0.16	0.39	151	25.4
PVC	0.05	0.05	0.02	0.09	60	10/15/20/25

The FEM model of the blade is generated by using ANSYS Parametric Design Language (ADPL) in ANSYS software, which enables the creation of various FEM models by changing the main aerodynamic and structural parameters. The model is entirely created with SHELL99 element for modeling the spar cap and SHELL91 element for modeling the thick sandwich structures, namely the leading edge, trailing edge and shear webs. The two types of elements model the layup as orthotropic in given layers by real constants, which contain the material properties, orientation and thickness. The number of real constants should be reasonably determined, as a small number of real constants could result in big calculation errors, while a large number of them would be time- consuming. Based on a reasonable simplification of the layup, 290 real constants are defined through adjusting the model many times, as shown in Figure 5, each color represents a real constant. In order to prevent erroneous results, a regular quadrilateral mesh generation method is used to generate elements with low aspect ratios. The FEM model with a mass of 6555.2 kg is shown in Figure 6. The validation process of the FEM model had been carried out in our previous work [15] to guarantee the reliability of the numerical simulation.

Figure 5. Real constant distribution of the blade.

Figure 5. Real constant distribution of the blade.

Figure 6. Finite Element Method (FEM) model of the blade.

Figure 6. Finite Element Method (FEM) model of the blade.

In order to determine the flap-wise, edge-wise and torsional rigidities of the FEM model, two concentrated forces are applied at the tip in the corresponding directions, respectively. From the deformed results, the bending deflections and the node displacements are recorded, the bending angles and bending rate per unit length, the twist angle and corresponding rate of twist are calculated following the same approach used in [16]. Then, the flap-wise and edge-wise rigidities can be computed as the ration between the moment M acting on a given cross section and the rate of rotation dθ/dz; the torsional rigidity can be computed by the ration between the torque T transmitted across a given section and the corresponding rate of twist dφ per unit of length dz [8,17]. Figure 7 shows the comparisons between the rigidities of the FEM model and the original blade. As can be seen, the rigidities are in good agreement in the flap-wise, edge-wise and torsional directions.

Figure 7. (a) comparison of the flap-wise rigidity; (b) comparison of the edge-wise rigidity; (c) comparison of the torsional rigidity.

Figure 7. (a) comparison of the flap-wise rigidity; (b) comparison of the edge-wise rigidity; (c) comparison of the torsional rigidity.

2.2. Aerodynamic Loads of the Blade

Aerodynamic loads for design are computed using Blade Element Momentum (BEM) theory [18,19], which include the operational and the ultimate load cases. The operational load case considers the wind condition corresponding to the maximum root flap bending moment under operational state. The range of flow speeds from cut-in to cut-out is divided into different values every 0.5 m/s, and the flap bending moment under each wind speed is calculated using Equation (1):

$$M_{flap}={\displaystyle {\int }_0^{R}4\pi \rho u^{2}a}(1-a)F{r}^{2}dr$$

(1)

where $\rho$ is the air density, $u$ is the wind speed, $a$ is the axial induction factor, F is the Prandtl’s correction factor and R is the blade length.

The maximum value is derived from all the flap bending moments experienced over the range of flow speeds, then the corresponding wind speed $u$ can be obtained. Afterwards, the force normal to the rotor plane $d{p}_N$ and the force tangential to the rotor plane $d{p}_T$ are calculated as follows:

$$\{\begin{array}{l}d{p}_N=\frac{1}{2}\rho W^{2}c(r)(C_l\mathit{cos}\varphi +C_d\mathit{sin}\varphi )\mathit{dr}\\ d{p}_T=\frac{1}{2}\rho W^{2}c(r)(C_l\mathit{sin}\varphi -C_d\mathit{cos}\varphi )\mathit{dr}\end{array}$$

(2)

Where W is relative velocity composed of the axial velocity and the tangential velocity, $c(r)$ is the local chord, C_l and C_d are the lift and drag coefficients, $\varphi$ is the angle between the plane of rotation and the relative velocity.

The ultimate load case takes the gust speed with a return period of 50 years into account to make the wind turbine last for the design period. In this case, the blades are parked or idling and the ultimate load can be calculated approximately by empirical formula [20]:

$$p(r)=q_{2s}{C}_{f}c(r)$$

(3)

where C_f is a force coefficient, and $q_{2s}=\frac{1}{2}\rho V_{2s}^2$ is the dynamic pressure from an extreme wind speed time averaged over 2 s. For a homogeneous terrain, V_2s can be computed using:

$$V_{2s}=V_b{k}_t\left[\mathit{ln}\left(\frac{h_{\mathit{hub}}+2/3R}{z_0}\right)+3\right],$$

(4)

where $V_b=27m/s$ is a basis wind speed, $h_{\mathit{hub}}$ is the hub height, $z_0$ is the roughness length and $k_t$ is a terrain factor. If the surrounding landscape has no nearby obstacles and very low vegetation, the roughness length is approximately 0.01 m and the terrain factor is 0.17.

The operational and the ultimate load distributions of the FEM model are shown in Figure 8.

Figure 8. (a) Operational load distribution of the FEM model; (b) Ultimate load distribution of the FEM model.

Figure 8. (a) Operational load distribution of the FEM model; (b) Ultimate load distribution of the FEM model.

3. Formulation of the Optimization Problem

3.1. The Objective Functions

A successful blade design must satisfy a wide range of objectives, such as maximization of the AEP and power coefficient, resistance to extreme and fatigue loads, restriction of tip deflections, avoiding resonances, and minimization of weight and cost, some of these objectives are in conflict [2]. In order to make wind energy more competitive with other energy sources, manufacturers are concentrating on increasing the energy output capacity and bringing down the cost of blades and other components at the same time, so minimizing the COE is an attractive target pursued by many researchers. However, the cost of the blades involves many factors such as cost of the materials, production tooling, manufacturing labor and overland transportation, and some of them are hard to estimate. According to some studies [21,22], reducing materials to minimize the blade mass can not only cause cost reduction of the blade but also have a multiplier effect through out the system including the foundation. Thus, the maximum AEP and the minimum blade mass are taken as the objective functions.

The energy that can be captured by a wind turbine depends upon the power versus wind speed characteristic of the turbine and the wind-speed distribution at the turbine site. The wind-speed distribution is represented by the Weibull function. The AEP is defined as in Equations (5)–(7):

$$f(1)=\mathit{AEP}={\displaystyle \sum _{i=1}^N\frac{1}{2}}[P(u_i)+P(u_{i+1})]·f(u_i<u_0<u_{i+1})\times 8760$$

(5)

$$P(u_i)={\displaystyle {\int }_0^{R}4\pi r^3\rho u_i{\omega }^2{a}^{\prime }(1-a)\mathit{Fdr}}$$

(6)

$$f(u_i<u_0<u_{i+1})=\mathit{exp}\left(-{\left(\frac{u_i}{A}\right)}^k\right)-\mathit{exp}\left(-{\left(\frac{u_{i+1}}{A}\right)}^k\right),$$

(7)

where $P(u_i)$ is the power under wind speed $u_i$, $f(u_i<u_0<u_{i+1})$ is the probability that the wind speed lies between $u_i$ and $u_{i+1}$, $\omega$ is the rotational speed of the blade, $a^{\prime }$ is the tangential induction factor, $k$ is a form factor and $A$ is a scaling factor.

The blade mass depends upon the consumption of materials, which can be expressed as:

$$f(2)=\mathit{mass}={\displaystyle \sum _i{\rho }_i\times V_i}$$

(8)

where ${\rho }_i$ is the material density, V_i is the volume of the material.

3.2. Design Variables

The geometrical shape and the rotational speed of the blade contribute directly to the aerodynamic performances. Therefore, the values of the control points in Figure 2 and Figure 3 and the rotational speed are selected as aerodynamic design variables, while the blade length and the airfoil series are employed as the basic parameters. In order to match up with the bolted flange conveniently, the root diameter remains unchanged, which means the chord value of CP1 is fixed. Moreover, the chord of CP2 is equal to CP1, so the chord values of CP3-7 are defined as five aerodynamic variables ($x_1$ to $x_5$). As the twist remains constant inboard of the maximum chord, which means the twists of CP1-3 are the same, then the twist values of CP3-7 are defined as another five aerodynamic variables ($x_6$ to $x_{10}$). The span-wise locations of CP8 and CP13 are also fixed, one is at root and the other is at tip, thus the locations of CP9-12 (airfoils with relative thickness of 40%, 30%, 25% and 18%) are employed as four more aerodynamic variables ($x_{11}$ to $x_{14}$). In addition, the rated rotational speed of the rotor is selected as the last aerodynamic variable ($x_{15}$).

The spar cap of the blade is primarily designed with a relatively large thickness to carry the aerodynamic loads, so it makes a major contribution to the blade mass [23]. In addition, the blade mass can further decrease by positioning the shear webs appropriately. Therefore, the number and location of layers in the spar cap, the width of the spar cap, and the position of the shear webs are selected as the structural design variables.

Figure 9 shows the original material layup in the spar cap. The region from 4.4 m to 25.3 m (shown in green) will be optimized as it has a much greater number of layers than the other regions. Eight control points are used to simulate the layup in the selected region, and the number of layers are defined to change linearly between the control points, as shown in Figure 10. CP16-19 each has two parameters that are the number and the location of layers, while the other CPs each have one parameter that is the number of layers. In addition, CP17 and CP18 have the same number of layers. Thus, the number of layers of CP14-21 is defined as seven structural variables (x₁₆ to x₂₂), the location of layers of CP16-19 are defined as another four structural variables (x₂₃ to x₂₆). Two more structural variables (x₂₇ and x₂₈) are used to define the width of the spar cap and the position of the shear webs, respectively, as shown in Figure 4.

Figure 9. Original material layup of the spar cap.

Figure 9. Original material layup of the spar cap.

Figure 10. Parametric material layup of the spar cap.

Figure 10. Parametric material layup of the spar cap.

Twenty-eight variables in total are defined, which can be expressed in the following form:

$$X={[\begin{array}{cccc}{x}_1& x_2& \dots & x_n\end{array}]}^T,n=28.$$

(9)

3.3. Constraint Conditions

The blade design is a multi-criteria constrained optimization problem, and the aerodynamic and structural requirements should be well satisfied [24,25]. In this paper, the following constraint conditions are mainly taken into account: the strain, the tip deflection, the vibration and the buckling constraints. These constraints represent the strength, stiffness, dynamic behavior, and stability design requirements, respectively.

(1) Generally, the stress and the strain criterions should be considered to reflect the strength of the blade. However, the maximum stress generated in the considered 1.5 MW blade under different loads is far less than the allowable stress, while the maximum strain is much closer to the design value. Therefore, only the strain criterion is used to verify that no elements in the model exceeded the design strains of the material [22]. This is expressed as follows:

$${\epsilon }_{\mathit{max}}\le {\epsilon }_d/{\gamma }_{S1},$$

(10)

where ${\epsilon }_{\mathit{max}}$ is the maximum equivalent strain, ${\epsilon }_d$ is the permissible strain, and γ_s1 is the strain safety factor.

(2) In order to prevent the collision between the blade tip and the tower, the maximum tip deflection should be less than the set value. This can be expressed as follows:

$$d_{\mathit{max}}\le d_d/{\gamma }_{S2},$$

(11)

where d_max is the maximum tip deflection, d_d is the allowable tip deflection, and γ_s2 is the tip deflection safety factor.

(3) In order to avoid resonance, the natural frequency of the blade should be separated from the harmonic vibration associated with rotor rotation. This is expressed in the inequality form:

$$\left|F_{\mathit{blade-1}}-3F_{\mathit{rotor}}\right|\ge \Delta ,$$

(12)

where F_blade-1 is the first natural frequency of the blade, F_rotor is the frequency of the rotor rotation, and Δ is the associated allowable tolerance.

(4) Since the blade is a thin-walled structure, its surface panels near the root are particularly vulnerable to elastic instability, and the buckling problem must be addressed [23]. In order to avoid buckling failure, the buckling load should be greater than ultimate loads. This can be expressed as follows:

$${\lambda }_1\ge 1.0\times {\gamma }_{S3},$$

(13)

where λ₁ is a ratio of the buckling load to the maximum ultimate load, called lowest buckling eigenvalue, and γ_s3 is the buckling safety factor.

In addition, the lower and upper bounds of design variables should be set appropriately to control the change range, shown in Equation (14a). The twist, chord, and relative thickness distributions are all required to be monotonically decreasing, so the corresponding variables should be satisfied with the inequality form in Equation (14b,c). Considering the manufacturing maneuverability and the continuity of the material layup, the variables that represent the number of layers are required to be monotonically increasing to a maximum value and then monotonically decreasing, shown in Equation (14d,e).

$$\{\begin{array}{lll}{x}_i^L\le x_i\le x_i^U& i=1,2,\mathrm{...},28& (a)\\ x_j-x_{j+1}>0& j=1,2,3,4,6,7,8,9& (b)\\ x_k-x_{k+1}<0& j=11,12,13& (c)\\ x_g-x_{g+1}\le 0& g=16,17,18& (d)\\ x_h-x_{h+1}\ge 0& h=19,20,21& (e)\end{array},$$

(14)

where $x^L$ is the lower bound of the variables, and $x^U$ is the upper bound of the variables.

The lower and upper bounds of the variables and the constraint conditions are shown in Table 3.

Table 3. Lower and upper bounds of the variables and the constraint conditions.

Parameter	Lower Bound	Upper Bound	Unit
x1-x5	1.0	3.3	m
x6-x10	−2.0	12.0	°
x11-x14	6.0	30.0	m
x15	10	25	rpm
x16	28	38	-
x17	28	48	-
x18	33	58	-
x19	35	65	-
x20	33	55	-
x21	28	45	-
x22	28	40	-
x23	7.0	9.0	m
x24	10.0	13.0	m
x25	16.0	19.0	m
x26	20.5	22.0	m
x27	0.13	0.25	m
x28	0.50	0.70	m
εmax	-	0.005	-
dmax	-	5.5	m
λ1	1.2	-	-
Fblade-1	≤3Frotor − 0.3 or ≥3Frotor + 0.3		Hz

Table 3. Lower and upper bounds of the variables and the constraint conditions.

Parameter	Lower Bound	Upper Bound	Unit
x1-x5	1.0	3.3	m
x6-x10	−2.0	12.0	°
x11-x14	6.0	30.0	m
x15	10	25	rpm
x16	28	38	-
x17	28	48	-
x18	33	58	-
x19	35	65	-
x20	33	55	-
x21	28	45	-
x22	28	40	-
x23	7.0	9.0	m
x24	10.0	13.0	m
x25	16.0	19.0	m
x26	20.5	22.0	m
x27	0.13	0.25	m
x28	0.50	0.70	m
εmax	-	0.005	-
dmax	-	5.5	m
λ1	1.2	-	-
Fblade-1	≤3Frotor − 0.3 or ≥3Frotor + 0.3		Hz

4. The Integrated Optimization Design Procedure

The purpose of the present work is to improve the overall performance of HAWT blades by means of an optimization of its aerodynamic and structural integrated design, a procedure that interfaces both the MATLAB optimization tool and develops the finite element software ANSYS. Three modules are used in the procedure: an aerodynamic analysis module, a structural analysis module and a multi-objective optimization module. The former two provide a sufficiently accurate solution of the aerodynamic performances and structural behaviors of the blade; the latter handles the design variables of the optimization problem and promotes functions optimization. The non-dominated sorting genetic algorithm (NSGA) II [26,27,28] is adapted for the integrated optimization design. It is one of the most efficient and well-known multi-objective evolutionary algorithms and has been widely applied to solve complicated optimization problems. According to the method, the Pareto optimal front can be obtained considering the set of Non-Dominated Solutions.

Figure 11 shows the flowchart of the integrated optimization process. After inputting the original parameters, an initial population is generated randomly in MATLAB with the aerodynamic and structural variables written in a string of genes. The aerodynamic variables are used to define the geometry shape of the blade while the structural variables are used to define the internal structure. Then, the BEM theory is applied to evaluate the aerodynamic performance, such as power, power coefficient, aerodynamic loads, etc. Meanwhile, a macro file that can transfer the variables from MATLAB into ANSYS is created using APDL language. Through some specific commands, the procedure opens the software ANSYS, calls the macro file to generate a parametric FEM model of the blade and simulates the load cases in order to obtain the structural behaviors, namely, the strain level, the tip deflection, the natural frequencies, the bulking load factors, etc. After several constraint conditions have been checked, two appropriate fitness functions are evaluated. The next step is to classify the solutions according to a fast non-dominated sorting approach and assign the crowding distance. Finally, a new population is created and the process can restart until the optimization procedure converges.

Figure 11. Flowchart of the integrated optimization process.

Figure 11. Flowchart of the integrated optimization process.

5. Optimization Application and Results

The basic parameters of the rotor, wind condition and NSGA-II algorithm are listed in Table 4. The Weibull form factor and scaling factor are determined from local meteorological data of inland China with an annual average wind speed of 6 m/s, namely, k = 1.91 and A = 6.8 m/s. The probability distribution of the specified wind speed is shown in Figure 12.

Table 4. Parameters of rotor, wind condition and non-dominated sorting genetic algorithm (NSGA)-II algorithm.

Parameter	Value	Unit
Rotor diameter	77	m
Number of blades	3	-
Hub diameter	3	m
Hub height	75	m
Rated wind speed	12	m/s
Rated rotational speed	19	rpm
Rated power	1500	kW
Cut-in wind speed	4	m/s
Cut-out wind speed	25	m/s
Air density	1.225	kg/m3
Weibull form factor k	1.91	-
Weibull scaling factor A	6.8	m/s
Number of individuals	40	-
Number of iterations	30	-
Probability of crossover	0.8	-
Probability of mutation	0.05	-

Table 4. Parameters of rotor, wind condition and non-dominated sorting genetic algorithm (NSGA)-II algorithm.

Parameter	Value	Unit
Rotor diameter	77	m
Number of blades	3	-
Hub diameter	3	m
Hub height	75	m
Rated wind speed	12	m/s
Rated rotational speed	19	rpm
Rated power	1500	kW
Cut-in wind speed	4	m/s
Cut-out wind speed	25	m/s
Air density	1.225	kg/m3
Weibull form factor k	1.91	-
Weibull scaling factor A	6.8	m/s
Number of individuals	40	-
Number of iterations	30	-
Probability of crossover	0.8	-
Probability of mutation	0.05	-

Figure 12. Weibull probability distribution of wind speed.

Figure 12. Weibull probability distribution of wind speed.

Figure 13 shows the Pareto front obtained by taking the maximum AEP and the minimum blade mass as the optimization objectives. The points on the left hand of the front identify design solutions having low mass yet low energy production, whereas the points on the right hand of the front present design solutions having high energy production but also high mass, which indicates there are some conflicts between the two objectives. It cannot be said which point is the best, the choice of the solution in the practical design should be made according to the designer’s favor.

In order to explain the formation of the Pareto front better, four optimized schemes extracted at different positions on the Pareto front are analyzed, as marked in Figure 13. The values of the design variables of schemes A, B, C, D and the original scheme are listed in Table 5.

Figure 13. Pareto front of two objectives.

Figure 13. Pareto front of two objectives.

Table 5. Values of the design variables.

Parameter	Original Scheme	Scheme A	Scheme B	Scheme C	Scheme D	Unit
x1	3.08	2.70	2.60	2.47	2.44	m
x2	2.88	2.64	2.56	2.43	2.41	m
x3	2.30	2.21	2.15	2.10	2.08	m
x4	1.82	1.66	1.64	1.61	1.60	m
x5	1.21	1.16	1.16	1.19	1.17	m
x6	8.73	10.69	10.38	9.70	9.56	°
x7	6.64	6.72	6.61	6.59	6.60	°
x8	3.14	3.84	3.49	3.80	3.58	°
x9	0.43	0.85	0.82	0.88	0.97	°
x10	−1.13	−0.49	−0.81	−0.22	0.02	°
x11	6.75	7.73	6.96	7.10	7.03	m
x12	9.50	9.81	9.32	10.51	9.73	m
x13	14.10	14.18	14.32	14.34	14.21	m
x14	28.95	24.32	25.16	24.35	24.20	m
x15	19.0	15.5	15.6	15.2	14.9	rpm
x16	33	32	30	30	29	-
x17	43	38	36	35	34	-
x18	53	48	46	44	44	-
x19	62	53	51	50	48	-
x20	53	41	39	38	38	-
x21	43	33	31	30	30	-
x22	33	29	29	29	28	-
x23	7.8	7.6	8.0	7.6	7.8	m
x24	11.0	12.1	11.8	11.9	11.8	m
x25	18.0	17.9	17.9	17.7	17.7	m
x26	21.4	21.3	21.3	21.3	21.2	m
x27	0.188	0.195	0.180	0.193	0.191	m
x28	0.620	0.589	0.558	0.553	0.551	m

Table 5. Values of the design variables.

Parameter	Original Scheme	Scheme A	Scheme B	Scheme C	Scheme D	Unit
x1	3.08	2.70	2.60	2.47	2.44	m
x2	2.88	2.64	2.56	2.43	2.41	m
x3	2.30	2.21	2.15	2.10	2.08	m
x4	1.82	1.66	1.64	1.61	1.60	m
x5	1.21	1.16	1.16	1.19	1.17	m
x6	8.73	10.69	10.38	9.70	9.56	°
x7	6.64	6.72	6.61	6.59	6.60	°
x8	3.14	3.84	3.49	3.80	3.58	°
x9	0.43	0.85	0.82	0.88	0.97	°
x10	−1.13	−0.49	−0.81	−0.22	0.02	°
x11	6.75	7.73	6.96	7.10	7.03	m
x12	9.50	9.81	9.32	10.51	9.73	m
x13	14.10	14.18	14.32	14.34	14.21	m
x14	28.95	24.32	25.16	24.35	24.20	m
x15	19.0	15.5	15.6	15.2	14.9	rpm
x16	33	32	30	30	29	-
x17	43	38	36	35	34	-
x18	53	48	46	44	44	-
x19	62	53	51	50	48	-
x20	53	41	39	38	38	-
x21	43	33	31	30	30	-
x22	33	29	29	29	28	-
x23	7.8	7.6	8.0	7.6	7.8	m
x24	11.0	12.1	11.8	11.9	11.8	m
x25	18.0	17.9	17.9	17.7	17.7	m
x26	21.4	21.3	21.3	21.3	21.2	m
x27	0.188	0.195	0.180	0.193	0.191	m
x28	0.620	0.589	0.558	0.553	0.551	m

Figure 14 shows the chord distributions of the original scheme and optimized schemes. Due to the blade root diameter remaining the same, the chords of the optimization schemes change slightly in the root area. After the root area, the chords of the optimized schemes show an obvious decrease when compared to the original scheme, especially in the maximum chord region and near the blade tip. This indicates that the original scheme is possibly designed by the conventional methods (Wilson method or Glauert method), which leads to a larger chord. The reduction of the chord could reduce the power production to some extent, but it can result in a lighter blade and a lower thrust at the same time.

Figure 14. Comparison of chord distribution.

Figure 14. Comparison of chord distribution.

Figure 15 shows the twist distributions of the five schemes. The twists of the optimized schemes also change slightly in the root area for structural reasons. Then, the twists increase from the maximum chord region, but the distribution trends almost remain the same as the original scheme. Due to the decreasing of the rotational speed ($x_{15}$ listed in Table 5), the angle between the mean relative velocity and tangential direction will increase. With the incensement of the twist, the angle of attacks can almost remain the same, which guarantees the airfoils still have higher lift-to-drag ratios to partly make up for the power loss caused by the reduction of the chord.

Figure 15. Comparison of twist distribution.

Figure 15. Comparison of twist distribution.

The percent thickness distributions are shown in Figure 16. The locations of the airfoils with 40% and 30% thicknesses for the optimized schemes move toward the tip. From the structural point of view, this is good for increasing the section moment of inertia. The location of the airfoil with 25% thickness almost remains the same, while that of the airfoil with 18% thickness moves toward the root. As the airfoil with 18% thickness has a higher lift-to-drag ratio, the aerodynamic region composed by this airfoil is the main part of the blade to capture wind energy. Moving the location of the airfoil with 18% thickness towards the root can increase the length of this region, which is beneficial to capture more energy, thus improving the power efficiency from the aerodynamic point of view. The blade shapes become much smoother after optimization, which is more convenient for manufacturing. The decreasing of the rotational speed could reduce the rotor rotation times during a 20 year life span, and thus improve the fatigue life of the blade.

Figure 16. Comparison of percent thickness distribution.

Figure 16. Comparison of percent thickness distribution.

Figure 17 and Figure 18 show the comparison of power coefficients and powers between the original scheme and optimized schemes, respectively. When compared to the original scheme, the power coefficients of the optimized schemes increase significantly in the speed range of 4 to 8 m/s, and decrease slightly in the peed range of 10 to 13 m/s. The power coefficients of the optimized schemes all reach the maximum values of about 0.49 at 8 m/s wind speed (a 7.7% increase compared to the original scheme), and keep relatively high values in the speed range of 6 to 9 m/s. With the reduction of chord, the rated wind speeds of the 4 optimized schemes gradually change from 12 to 15 m/s, as shown in Figure 15. The comparisons in Figure 17 and Figure 18 reveal that the optimized schemes have better aerodynamic performances than the original scheme at low wind speeds, but the performances gradually become worse when the wind speed increases. As can be seen from Figure 12, for a turbine site with an annual average wind speed of 6 m/s, the probabilities of the wind speed lying between 4 m/s and 8 m/s are much higher than those lying between 10 m/s and 13 m/s, so the optimized schemes are more reasonable as they can utilize more wind energy resources.

Figure 17. Comparison of power coefficient under different wind speed.

Figure 17. Comparison of power coefficient under different wind speed.

Figure 18. Comparison of power under different wind speed.

Figure 18. Comparison of power under different wind speed.

Figure 19 shows the comparison of material layup in the spar cap between the original scheme and optimized schemes. From Figure 19 and Table 4, it can be seen that the number of layers changes slightly from 4.4 m to 8 m and 21.5 m to 25.3 m while it decreases obviously from 8 m to 21.5 m along span-wise of the blade, and the thickest region becomes smaller. This indicates that the two regions from 4.4 m to 8 m and 21.5 m to 25.3 m have less impact on the blade mass than the middle part. The reason is that the original number of layers in these two regions is less, and the relatively large lower bounds limit it to change significantly. As the region from 4.4 m to 8 m withstands higher loads, the number of layers in this region is a bit more than it in the region from 21.5 m to 25.3 m after optimization.

The position of the shear webs decrease after optimization, which means that the shear webs move toward the centerline of the spar cap. According to the sensitive analysis in our previous work [29], although the blade mass will increase from 6555.2 kg to 6560.2 kg when the position was reduced by 20%, the maximum strain can reduce from 0.00429 to 0.00412, i.e., improve the blade strength on the premise that the mass is almost unchanged. Thus, the number of layers in the spar cap could decrease further, which makes the blade much lighter. The width of the spar cap also decreases after optimization, and a smaller width of the spar cap can reduce the amount of materials, which is beneficial for reducing the blade mass.

Figure 19. Comparison of material layup of the spar cap.

Figure 19. Comparison of material layup of the spar cap.

The AEP, blade mass and structural performances of the schemes are presented in Table 6. Compared with the original scheme, the AEP of schemes A, B, C, D increase by 11.40%, 10.99%, 9.48% and 7.56%, respectively, while the blade mass decreases by 10.39%, 13.53%, 15.33% and 15.96%, respectively. The maximum strain increases, the tip deflection decreases firstly and then increases, and the lowest buckling load factor decreases, but they still satisfy the constraint conditions set in the procedure. The structural stiffness reduces with the decreasing of the blade. However, the reduced ratio of the stiffness is greater than that of the mass of the blade. Therefore, the first natural frequency decreases, but there is no occurrence of resonance. Compared with the best result having a blade mass of 6064.6 kg in [15], the reduction of chords and spar width can further decrease the blade mass obviously in the case of a small difference between the material layups.

Table 6. Comparison of the annual energy production (AEP), blade mass and structural performance.

Scheme	AEP (GWh/Year)	Blade Mass (kg)	Maximum Strain	Maximum Tip Deflection (m)	The Lowest Buckling Eigenvalue	The First Natural Frequency (Hz)
Original	3.175	6555.2	0.00429	4.60	2.024	1.027
A	3.537	5874.1	0.00481	4.13	1.491	0.969
B	3.524	5668.4	0.00485	4.36	1.358	0.938
C	3.476	5550.3	0.00491	4.75	1.252	0.907
D	3.415	5509.1	0.00498	4.92	1.216	0.884

Table 6. Comparison of the annual energy production (AEP), blade mass and structural performance.

Scheme	AEP (GWh/Year)	Blade Mass (kg)	Maximum Strain	Maximum Tip Deflection (m)	The Lowest Buckling Eigenvalue	The First Natural Frequency (Hz)
Original	3.175	6555.2	0.00429	4.60	2.024	1.027
A	3.537	5874.1	0.00481	4.13	1.491	0.969
B	3.524	5668.4	0.00485	4.36	1.358	0.938
C	3.476	5550.3	0.00491	4.75	1.252	0.907
D	3.415	5509.1	0.00498	4.92	1.216	0.884

6. Conclusions

This paper describes a multi-objective optimization method for the aerodynamic and structural integrated design of HAWT blades. The method is used to obtain the best trade-off solutions between the two conflict objectives of maximum AEP and minimum blade mass.

A procedure uses BEM theory and an FEM model, and the NSGA II algorithm is developed for this purpose. The BEM theory and FEM model are utilized to determine the aerodynamic and structural performances of HAWT blades, and the optimization fitness functions as well. The NSGA II algorithm is adopted to handle the design variables chosen for optimization and search for the group of optimal solutions following the basic principles of Genetic Programming and Pareto concepts.

The procedure has been applied successfully to a 1.5 MW commercial HAWT blade, and satisfactory schemes to increase the AEP and decrease the blade mass are achieved under a specific annual average wind speed. The results indicate that the maximum AEP requires blades having large chords and layer thicknesses, thus high masses, while the requirement of the maximum blade mass is just the opposite. The further aerodynamic and structural analysis of the optimization schemes show great improvements for the overall performances of the blade, which indicates the efficiency and reliability of the proposed procedure.

In future work, an attempt will be made to define the chord and the twist distributions more appropriately with Non-Uniform Rational B-Splines (NURBs) curves. The optimization of the material layup in the whole spar cap will also be carried out to further improve the overall performances of the blade.