Design and Mass Optimization of Numerical Models for Composite Wind Turbine Blades

Abstract

In this paper, a constrained optimization by linear approximation (COBYLA) algorithm is used to optimize the design of a 5 MW wind turbine blade. In the process of blade material modeling, the actual manufacturing conditions are considered, and the load of blades under 50 m/s wind conditions is analyzed based on the blade element momentum (BEM) method. Mass optimization was achieved by removing material from the shear webs. In addition, constraints such as tip displacement, stress, and frequency during blade design were considered. The results show that the mass is reduced by about 1.7% after removing material from blade webs, while the structural response of the blade remains unchanged. This case provides a practical reference for commercial wind turbine blades.

1. Introduction

Clean energy such as wind, solar, and hydrogen have been widely developed as the global energy crisis has intensified. As the main pneumatic components of wind turbine, the structural performance of wind turbine blades will directly affect the operating stability. In the structural design of a wind turbine, the geometrical shape of blades is determined by aerodynamic factors and generally cannot be modified [1]. To ensure sufficient structural stiffness and reliability while reducing the weight and cost of blades has become a research hotspot [2,3,4,5]. As the blades of modern large-scale wind turbines are made from a variety of composite materials, including glass fiber, foam, resin, and others, the optimization of their structure will be an extremely complicated problem.

Excitingly, the structural design of composite shells and wind turbine blades has attracted more and more researchers’ interest. Ma et al. [6] proposed a combinatorial optimization design method for wind turbine blades by PSO algorithm, which not only improves the power of wind turbines, but also reduces the manufacturing cost of blades. Based on simulation optimization and artificial neural network optimization methods [7], Albanesi et al. [8] determined the optimal layup sequence, layup quantity, and laydown decline scheme for lamination, while taking into account the displacement constraints after large elastic deformation, ultimately reducing blade mass and computational costs.

The topology optimization method [9] has also been applied to the structural design of wind turbine blades in related research. Albanesi et al. [10] constrained elastic displacement, stress, and natural frequency in the design and manufacture of 28.5 m composite blades, used genetic algorithms to code material positions to determine the optimal laminate stacking sequence in the blade skin, and used topological optimization to achieve lower blade mass. Taking rotor power coefficient and blade structure flexibility as aerodynamic performance and structure optimization objectives, Wang et al. [11] characterized shape optimization by the position of the NURBS curve, and realized topology optimization by the material densities assigned to blade cross-sections finite elements modeling. Sohouli et al. [12,13] proposed a discrete material optimization (DMO) method on the basis of a multi-phase topology optimization idea and made it possible to solve discrete shell optimization problems. Although the method is still remote from the actual manufacturing of laminates, it provides a different solution for the optimal design of composite shells. Song et al. [14] adopted the variable density topology optimization method to perform the internal structure design of the blade. In addition, they also used the CFD method to calculate the output torque and power of the wind turbine. The above studies show that the topological optimization method has a certain potential in solving the problems of material layout and structural optimization. However, topological optimization requires the interpolation of material properties (including Young’s modulus and density), which brings challenges to computational resources; because wind turbine blades are stacked with carbon fiber, glass fiber, foam, and other materials, the maximum number of layers in the design process may reach hundreds of layers [15,16].

In addition, the genetic algorithm is also widely used in wind turbine blade and composite structure optimization design. Barnes and Morozov [17] studied the geometrical structure of the blade thoroughly by modifying the width of the beam cap, the trailing edge stiffeners material, and the number and position of the shear web, coded different positions and thicknesses of the material, and used a genetic algorithm for parameter optimization to effectively achieve mass reduction. Wang et al. [18] considered five constraints (including stress, tip deformation, vibration, and continuity of laminate layups) in wind turbine design optimization based on finite element analysis and the genetic algorithm model; a kind of structural optimization approach for composite wind turbine blades was developed and successfully applied to 30 kW wind turbine. Nicholas [19] coded composite lamination angle, used a neural network instead of finite element analysis, and combined an artificial neural network with a genetic algorithm to achieve structural adjustment. The research of the above scholars effectively realizes the mass and layup optimization of wind turbine blades; since genetic algorithm is a random optimization algorithm, the local optimal situation may occur, and the prediction effect of artificial neural network can also be affected by the training samples.

The constrained optimization by linear approximation (COBYLA) algorithm does not require the evaluation of the objective function gradient to converge, and can reduce the errors caused by gradient evaluation. Therefore, a lot of researchers began to use the COBYLA algorithm for optimization design [20,21]. Selimefendigil and Öztop [22] used the COBYLA optimization algorithm to find the optimum location and size of the porous object in order to optimize the process of step convective heat transfer. Altieri et al. [23] proposed a reliability optimization method for viscous dampers in building frames subjected to a stochastic earthquake input with uncertain intensity, duration, and frequency characteristics. Mam et al. [24] adopted the COBYLA algorithm to solve the shape optimization problem of high-rise structures under wind and gravity loads. It can be seen that the COBYLA algorithm has great potential in shape optimization.

Benefiting from the advantages of the COBYLA algorithm in shape optimization, the present work targets commercially available wind turbine blades with complex layers. The shape optimization method is used to remove the material from the blade shear webs, and the mass optimization is quickly achieved without changing the structural performance of the wind turbine. In this case, the material radius, which was removed, is used as a variable, the tip displacement and maximum stress are used as optimization constraints, and the blade mass is designed as the optimization objective to conduct this research.

2. The Structure and Load of Blade

2.1. Geometry Description

The blade geometry model adopts the 5 MW model [25,26] of the National Renewable Energy Laboratory (NREL), which is composed of six types of airfoil, each of them with different twist angles along the spanwise position to adjust the aerodynamic performance. The overall length of the blade is 61.5 m, which is made of glass fiber, carbon fiber, and foam after stacking. The shell structure is divided into a leading edge, leading edge panel, spar cap, trailing edge panel, and other different areas to arrange the laminate layout reasonably. Meanwhile, to ensure that the blade has strong shear resistance, we designed a double web inside the blade. The structure of the blade is shown in Figure 1a, and the orientation of the laminate fiber is shown in Figure 1b. As can be seen from Figure 1a, the web divides the shell into three regions in spanwise direction, among which the spar cap of the blade is located between two shear webs. Considering that the laminates are laid inward along the mold during the manufacture of wind turbine blades [27], the blade materials are prescribed to be stacked inward along the shell structure in order to be closer to production.

2.2. Layup Design and Aerodynamic Loads

In order to obtain a more reasonable layup design scheme, this paper refers to a layup scheme of a 5 MW wind turbine which is close to commercial, so that it is closer to the engineering application. In the composite material stacking, the blade materials including gelcoat, triaxial fibers (SNL Triax), uniaxial carbon fiber (Carbon), uniaxial glass fiber (E-LT-5500), foam, and so on [28]; specific properties are shown in

Table 1. Material properties of the composite materials.

Material	Thickness (mm)	EX (MPa)	Ey (MPa)	Gxy (MPa)	Prxy	Dens (kg/m3)	UTSL (MPa)	UCSL (MPa)
Gelcoat	0.05	3440	\	1380	0.3	1235	\	\
E-LT-5500 (Uniaxial)	0.47	41,800	14,000	2630	0.28	1920	972	702
SNL Triax	0.94	27,700	13,650	7200	0.39	1850	700	\
Saertex (Uniaxial)	1	13,600	13,300	11,800	0.49	1780	144	213
Foam	1	256	256	22	0.3	200	\	\
Carbon (Uniaxial)	0.47	114,500	8390	5990	0.27	12,200	1546	1047

.

Wind turbine blades use gelcoat on the outermost layer, and the spar cap is reinforced with carbon fiber, which is lighter and stronger than glass fiber. Although it may cost more, it is now also used in some large wind turbine blades. It has been shown that adding core to the composite layering can improve the buckling performance of the blade [29], so the leading edge and trailing edge panels are filled with foam to form a sandwich structure. For the triaxial fabric design, according to reference [30], the triaxial fabric is equivalent to $[+45/-45/0]$ stacking with a thickness ratio of 35%, 35%, and 30%, respectively, while uniaxial fabric can be equivalent to $[0/90]$ stacking with thickness ratio of 95% and 5%. Since the thickness of the triaxial fabric is 0.94 mm, the corresponding thickness of the three angles is designed to be 0.32, 0.32, and 0.3 mm, respectively. In this paper, the uniaxial fabric is composed of 0° fibers. In addition, blades are usually designed to decrease in thickness from root to tip [31], thus the foam thickness in this design is 150 mm from the root to the DU40 airfoil, and gradually decreases to 10% from the DU40 airfoil to the tip, resulting in a blade mass of 17,211 kg. The specific layering scheme can be referred to [15,28]. The blade layering structure designed in this paper is shown in Figure 2.

2.3. Load Scenarios over the Blade

As the state response under extreme wind conditions is commonly considered in the design process of a wind turbine, the IEC6.2 design standard is selected in this paper. The wind model is a Kalman turbulence model, turbulence class is B, and the reference wind speed is 50 m/s. The turbine is full-span pitch and the blades are feathered, so the pitch angle is 90 degrees in OpenFast. At the same time that the turbine is shut down, the blades are subject to both aerodynamic and gravitational loads.

The blade element momentum (BEM) theory [32] is a widely used method for the aerodynamic load calculation of blades. The inlet velocity of the blade is composed of the inlet wind speed and the rotation velocity, the local flow angle follows from the expression:

$$\varphi =\arctan \frac{U_{\infty }(1-a)}{\mathsf{\Omega }r(1+a^{\prime })}$$

(1)

where $\varphi$ is the local flow angle, $U_{\infty }$ is the free-stream wind speed, and $\mathsf{\Omega }$ represents the rotor rotational speed, $r$ is the distance from the cross-section to the center of turbine hub. Further, $a$ and $a^{\prime }$ are the axial and angular induction factors, respectively, and their relationship to $\varphi$ can be obtained from the following formula.

$$f(\varphi )=\frac{\sin \varphi }{1-a}-\frac{\cos \varphi }{{\lambda }_r(1+a^{\prime })}=0$$

(2)

where ${\lambda }_r$ is the local tip-speed ratio, which can be calculated by $\mathsf{\Omega }r/U_{\infty }$. When the local inflow angle $\varphi$ is determined, tangential force $F_t$ and normal force $F_n$ can be calculated. Furthermore, the angle of attack also needs to be calculated, given by $\alpha =\varphi -\beta$, where $\beta$ is the twist angle, with its distribution along the blade length given in [17].

Reynolds number $R_e$ is necessary to characterize fluid flow. To obtain the Reynolds number, the local inflow velocity $V$, which, in accordance with Figure 3, needs to be calculated first.

$$V=\sqrt{V_x^2+V_y^2}=\sqrt{U_{\infty }^2{(1-a)}^2+{(\mathsf{\Omega }r)}^2{(1+a^{\prime })}^2}$$

(3)

Subsequently, the resulting value of $V$ is used together with the chord length $l$ for computing the local Reynolds number as $\mathrm{Re}=\rho Vl/\mu$, where the density $\rho =1.225\mathrm{kg}/{\mathrm{m}}^3$ and dynamic viscosity $\mu =1.793\times {10}^{-5}\mathrm{Ns}/{\mathrm{m}}^2$.With the specific angle of attack $\alpha$ and $R_e$ being calculated, the lift coefficient $C_L$ and the drag coefficient $C_D$, which are two coefficients related to the shape and surface properties of the object, need to be calculated. Then the lift force $L$ and drag force $D$ can be computed by $L=C_L\rho W^{2}l/2$ and $D=C_D\rho W^{2}l/2$. Finally, the normal force $Fn$ and tangential force $Ft$ in the plane of rotation can be obtained.

$$\{\begin{array}{c}Fn=\frac{C_L\rho W^{2}l}{2}\cos \varphi +\frac{C_D\rho W^{2}l}{2}\sin \varphi \\ Ft=\frac{C_L\rho W^{2}l}{2}\sin \varphi -\frac{C_D\rho W^{2}l}{2}\cos \varphi \end{array}$$

(4)

Since the aerodynamic center of the airfoil and the position of the pitch axis are generally different, the normal and tangential forces need to be transferred from center $O_1$ to pitch axis $O_2$. This transformation generates the following loads:

$$\{\begin{array}{c}Fx=Ft,Fy=Fn\\ M_Z=-d(Ft\sin \beta +Fn\cos \beta )\end{array}$$

(5)

where $d$ is the distance between $O_1$ and $O_2$, $M_Z$ is a torsional moment about pitch axis $O_2$. According to reference [10], the AeroDyn module provided by OpenFast was used to calculate aerodynamic load on blades. The mass of the element is computed with the term ${\left[{\displaystyle {\sum }_{j=1}^n{\rho }_j{h}_j}\right]}_i{A}_i$, where $A_i$ is the area of the element, ${\rho }_j$ and $h_j$ are layer and density and thickness, respectively.

In the AeroDyn module, analysis nodes are placed at 9 locations along the dimensionless span of the blade, which can be shown as [l/R] 0.022, 0.066, 0.11, 0.233, 0.366, 0.5, 0.633, 0.766, and 0.933. Tangential force and normal force of per unit length in the plane of the blade airfoil are applied to the blade. For a detailed coordinate diagram, refer to [33]. Figure 4 shows the extreme load conditions on the blade global coordinate system.

3. Design Optimization

3.1. The Optimization Model

Manufacturing costs in general are closely related to the weight of the blade, and lower blade mass means less cost. Therefore, the mass of the blade is taken as the optimization objective in this paper. It is considered feasible to remove material from the blades while ensuring adequate blade stiffness, as the laminates at the web are stacked more simpler and only were used to borne shear force. In this paper, the shape of the material removed from the shear web was set as a circle, and the mass was optimized by adjusting the circle’s radius. Simple shapes can reduce the use of computer resources to achieve optimization in a short time.

In this case, the circles were designed to remove material at five locations along the wingspan. The radii of the five circles were used as design variables, which can be expressed as:

$${\mathrm{R}=[\mathrm{R}}_1{ \mathrm{R}}_2 \dots {\mathrm{R}}_n{ ]}^T,n=5$$

(6)

Here, we set the maximum diameter of five circles to be less than 60% of the webs’ height and the corresponding maximum values for R₁, R₂, R₃, R₄, R₅ are 0.6, 0.8, 0.6, 0.4, 0.4 m, respectively.

To keep blade strength and stiffness unchanged, we constrained the maximum stress during optimization to be smaller than the maximum stress of the blade without shape optimization. Meanwhile, to avoid tower collision caused by over-tip displacement during optimal design of the wind turbine, the tip displacement will be constrained. According to the reference [19], the constrained maximum tip displacement in this case is lower than 7% of the radius of the wind turbine, that is, 4.3 m. Compared with 13.5% in reference [10], this case has a higher standard.

To sum up, the mathematical optimization model can be established:

$$\mathrm{Objective} \mathrm{function}: F=\min M(R)$$

$$\mathrm{subject} \mathrm{to}: Mises\le Mise{s}_{initial}$$

$$dis{p}_{\max }\le 0.07R_{Blade}$$

$$0\le 2R\le 0.6H_{SW}$$

where $M$ is the weight of the blade,$R$ is the radius of the material removed from two shear webs, $Mises$ represents the maximum Mises stress on wind turbine blades. $dis{p}_{\max }$ and $R_{Blade}$ correspond to the maximum tip displacement and wind turbine radius, respectively. $H_{SW}$ is the height of the shear web.

3.2. Optimization Algorithm

First, the finite element simulation was carried out for the wind turbine with intact shear webs to determine the positions that bear low loads. The maximum Mises stress of the blade was about 218 MPa, and the maximum elastic deformation was 1.09 m. In the initial simulation, due to the large load on the transition zone of the blade circular airfoil and DU40 airfoil, no material was removed from the web here.

The method used in this study is a gradient-free algorithm named COBYLA, which is optimized by adjusting the trust region radius. The detailed introduction can be found in reference [34,35]. In the initial optimization, we designed the values of five variable parameters as 0.3, 0.4, 0.25, 0.15, and 0.15 m, respectively. The composite layup was realized in COMSOL, and the gravity and aerodynamic load were applied. During iteration, the blade’s maximum stress and displacement were extracted as a constraint. When design conditions were met, the iteration was stopped. Figure 5 shows the specific optimization process.

4. Optimization Results

After iterative optimization, the optimal radii of material removed from the web were finally obtained, with no change in the maximum stress and tip displacement of the blade. The values of blade mass, maximum stress, and maximum tip displacement were recorded during the optimization process. Figure 6 shows the changes in tip displacement and radii of removed material. It can be seen from the figure that the tip displacement tends to be stable after 33 iterations in the optimization process, where the tip displacement stays at about 1.1 m, and the radii of the five resected holes (R1, R2, R3, R4, R5) were adjusted from 0.3, 0.4, 0.25, 0.15, 0.15 m to 0.6, 0.8, 0.6, 0.4, 0.4 m, respectively. Since the blade mass is inversely proportional to the radius of the removed hole, the adjusted R has reached the maximum value of the constraints, so it can be determined that the blade has reached optimal design.

During iteration, the hole diameter of the removed material contiguously changes, which will directly lead to a change of blade mass. Meanwhile, to ensure constant blade strength after optimization, the maximum stress is monitored during optimization. It should be noted that since the radii of removed material were specified in the initial optimization design, the initial blade mass recorded by monitoring was smaller than that (17211 kg) without material removal. In this case, the initial radius was about 15% of the web height, which is 0.3, 0.4, 0.25, 0.15, and 0.15 m as mentioned above.

From Figure 7a, it can be seen that the ratio between the maximum stress under optimization and that without optimization is extremely close to 1.0, which indicates that the whole stiffness of a blade does not be weakened during the process. As shown in Figure 7b, the blade mass gradually decreased, and after 25 steps of iteration, the mass tended to be stable, about 16,919 kg.

To compare the overall stress distribution of the blade in initial and optimization design, the stress nephogram of the blade is shown in Figure 8. It can be seen that the maximum stress of the DU40 airfoil is distributed at the spar cap and trailing edge. Therefore, structural strength of the transition zone from circular airfoil to aerodynamic airfoil should be considered in the design process of a wind turbine. After optimization, the maximum stress of the webs increased from 43.7 MPa to 56.9 MPa. Since the main load of the shell structure is borne by the spar cap and the main shear load is borne in the web, a slight change of stress in the web is considered acceptable under the condition that the whole stress of blade is constant.

Since the wind turbine has a long and thin structure, it can easily cause vibration. In the design process of a wind turbine, the natural frequency is also considered. To verify the availability of the optimized blade, a frequency analysis is carried out in this paper, as shown in

Table 2. Comparison between the reference and optimal blade designs.

Mode	Frequentry Reference [28]	Initial Design	Optimal Design	Mode	Frequentry Reference [28]	Initial Design	Optimal Design
1	0.87 Hz	0.91 Hz	0.88 Hz	4	3.91 Hz	3.46 Hz	3.42 Hz
2	1.06 Hz	1.09 Hz	1.06 Hz	5	5.57 Hz	5.35 Hz	5.26 Hz
3	2.68 Hz	2.59 Hz	2.55 Hz	6	6.45 Hz	6.73 Hz	6.49 Hz

. It can be seen that the first six-order natural frequencies of the blade do not change considerably in initial and optimal design, neither of which exceeds 10% of the reference value.

Figure 9 shows the vibration state of the blade after optimization. It can be seen from Figure 9 that the first six orders of the novel turbine blade include: 1st flapwise bending, 1st edgewise bending, 2nd flapwise bending, 2nd edgewise bending, 3rd flapwise bending, and 1st torsion, which are similar to reference [28].

5. Conclusions

An optimal design model was established according to the design parameters of a commercial 5 MW wind turbine blade under IEC6.2 design load. The aerodynamic load of wind turbine blades was calculated by AeroDyn, and COBYLA algorithm was applied to the structural parameter design of wind turbine blades. The blade mass was reduced by about 1.7%.

In the optimization process, the blade stress and tip displacement were constrained to minimize the blade mass. In addition, to further approximate the actual manufacturing, the material is stacked inwards along the manufacturing mold.

In this case, different laminates were stacked on the leading edge, leading edge panel, spar cap, trailing edge, and trailing edge panel of a wind turbine blade. The overall mass, natural frequency, and tip displacement of the blade all meet the design requirements and the cost is lower. At the same time, carbon fiber with great strength and low weight is used in the spar cap, which is in line with the current design trend of wind turbine blades. This case provides an engineering reference for wind turbine design. Future work will consider multi-objective optimization of blade mass and stiffness for commercial wind turbines.