Multi-Objective Structural Optimization of a Composite Wind Turbine Blade Considering Natural Frequencies of Vibration and Global Stability

: Aspects concerning resonance and global stability of a wind turbine blade must be carefully considered in its optimal design. In this paper, a composite wind turbine blade with an external geometry based on the NREL 5 MW model was subjected to multi-objective structural optimization considering these aspects. Four multi-objective structural optimization problems are formulated considering the blade mass, the maximum blade tip displacement, the natural frequencies of vibration, and the critical load factor as objective functions. The design variables are the number of plies, material, and ﬁber orientation. The design constraints are the materials’ margin of safety, the blade’s allowable tip displacement, and the minimum load factor. The blade model is submitted to the loads determined by the actuator lines theory and discretized in a ﬁnite element parameterized model using the Femap software according to geometric design variables. Among many multi-objective evolutionary algorithms available in the literature concerning evolutionary computation, the NSGA-II is the adopted evolutionary algorithm to solve the multi-objective optimization problems. Pareto fronts are obtained and performance indicators are used to evaluate the distribution of the non-dominated solutions. Multi-criteria decision-making is used to extract the solutions from the Pareto fronts according to the decision-maker’s preferences. The values of the objective functions, design variables, and constraints are presented for each extracted solution. The proposed study is expected to contribute to the multi-objective optimization and the structural design of wind turbine blades.


Introduction
Designing wind turbine blades or propellers is not a trivial task and involves multidisciplinary aspects of aerodynamic, structural, materials science, and optimized designs.Recently, several methods have been proposed for multidisciplinary optimization, including aerodynamic and structural performance.Generally, the formulations considering a single objective are the most usual and studied in the literature.
Optimizing the design of a wind turbine blade is extremely important in a more general context that can significantly impact projects such as wind farms, a strategic trend in many countries worldwide [1].A detailed review of optimization problem formulations, numerical models, strategies, and techniques concerning horizontal-axis wind turbines (HAWTs) is presented by Cherouri et al. in [2].
Usually, the objective functions considered in formulations of a single-or multiobjective optimization of a wind turbine are the minimization of the cost of energy (COE), maximization of the annual energy production (AEP), minimization of the blade mass, minimization of the rotor thrust or torque, and maximization of power, for instance.The design constraints can be geometric, such as ground clearance, blade tip deflection, and strains; aerodynamic, such as airfoil characteristics, maximum chord, noise levels; and physical, such as linear inequality, rated power, thrust, shaft torque, axial induction factor, stresses, the natural frequency of vibration, buckling, blade fatigue and damage, and static failure.The design variables can be grouped into two categories: the geometrical and physical parameters of the wind turbine blade.Among them, they can be listed such as rotor shape, airfoil characteristics and blade regulation, shell thickness and rotate rate, number of ribs and arrangement of the stiffening ribs, chord, twist, spar (thickness, location, and length), hub height, rotor speed, rotor diameter, and rated power.
Most of the structural optimization problems of wind turbine blades consider the total mass as an objective and constraints involving stresses, strains, and displacements (at the blade's tip) as constraints.And, most formulations consider only a single-objective: the blade's weight or mass [3,4].Also, commonly, the multi-objective structural optimization problems (MOSOPs) for wind turbine blades consider only two objectives.Generally, they seek to minimize the blade mass and the tip displacement.
Todoroki and Kawakami [5] adopted a multi-objective genetic algorithm (MOGA) to find the optimum design of a blade made of composite materials.The objective functions were the cost and mass of the blade.The design variables were skin thickness, number of carbon-fiber-reinforced polymer (CFRP), and glass-fiber-reinforced polymer (GFRP) layers in five regions.Ultimate strength, fatigue strength, and deflection were the considered constraints.To overcome the expensive computational cost of determining the constraints demanded by the analysis of the finite model of the blade, a Kriging model response surface approximation was used during the evolutionary process.
Wang et al. [6] proposed an improved NSGA to solve a MOSOP of a 5 MW wind turbine blade, in which the maximum power coefficient and the minimum blade mass were the optimization objectives.The design variables were the chord, the twist, and the thickness of five key sections.In addition to the geometry constraints, the power output was limited to the rated power output by pitch regulation when the wind speed is higher than the rated speed until the cut-out speed was achieved.Hu et al. [7] proposed a multiobjective approach considering the cost and mass of a HAWT composite blade as objective functions and layer thickness, material type, and orientation angle as design variables.The design constraints were derived from the ultimate strength analysis, fatigue failure, and critical deflection.A GA was adopted to solve the proposed MOsOP.A MOGA was used by Dal Monte et al. [8], in which the objective functions were the blade mass and the maximum deformation.The design variables were the choice of the employed materials and their placement in the layout of the blade skin.The constraints concerned the limits for the maximum mass and the maximum deformation of the blade.
The optimal structural design of the HAWT blade (1.5 MW), made of GRFP and CRFP materials, using the NSGA-II, was presented by Zhu et al. [9], in which the mass and the cost were the objective functions.The design variables were the number and location of layers in the spar cap, the spar cap's width, and the shear webs' position.In addition, the strain, blade/tower clearance, and vibration were set as design constraints.He and Agarwal [10] used the NSGA-II and jMetal [11] GA-based algorithms to optimize the shape of a well-known wind turbine airfoil S809.The idea was to improve its lift and drag characteristics by setting two objectives to increase its lift and lift-to-drag ratio subjected to the thickness, camber, and other geometric properties as constraints.A review concerning deterministic and metaheuristics algorithms, objective functions, design variables, design constraints, tools, and models for the performance evaluation of a HAWT was presented by Chehouri et al. in [2].
Gao et al. [12] adopted NSGA-II and fuzzy evaluation to solve the aerodynamic and structural coupling optimization of a 1.5 MW wind turbine blade.The AEP and blade mass were the objective functions.In addition, chord lengths and twist angles were set as design variables, whereas strength, stiffness, and stability of the blade as constraints.Using a weighted approach, the objective functions were combined into a single-objective optimization problem.Sessarego et al. [13] developed a concurrent-hybrid NSGA-II (hybrid NSGA-II) to the simultaneous optimization of the AEP, flapwise root-bending moment, and mass of the NREL 5 MW wind-turbine blade.Strain and displacements were set as constraints.Fagan et al. [14] used a NSGA-II and a finite element modeling to find the optimal structural design for a GFRP composite blade.The objective functions were the tip deflection and the blade mass.The number of bi-axial or unidirectional plies in the spar caps, the outer aerodynamic skin, the inner skin reinforcement, and the shear web were the design variables.
A multi-objective aerodynamic and structural optimization of a wind turbine blade using a novel adaptive game method was proposed by Meng et al. in [15].An integrated structural and aerodynamic optimization of a wind turbine blade was performed in a multi-objective optimization concerning aerodynamic and structural objective functions.The constraints were those of the blade's maximum equivalent strain, maximum displacement, and first-order natural frequency of vibration.An improved NSGA-II was proposed by Özkan and Genç [16].The objective functions were the blade mass and cost of the blade, both to be minimized.The composite material type and spar cap layer number were the design variables.The constraints were the deformation, strain, stress, natural frequency, and failure criteria.
Zhue et al. [17] presented an approach for the multi-objective aerodynamic and structural optimization of HAWT blades using the NSGA-II.The blade element momentum (BEM) theory was used to evaluate the aerodynamic performances, and the FEM method was adopted to evaluate the structural behaviors.The objective was to find a balance between the aerodynamic design and structural design to minimize the COE of HAWT blades.The AEP and the minimum wind turbine mass were considere as two conflicting objectives in the optimization problem discussed by Zhu et al. [18].Aerodynamic and structural parameters of the blade and tower were adopted as design variables.Design constraints included stress, strain, deflection, vibration and buckling limits.The NSGA-II was used to find the non-dominated solutions between the objectives.Additional references in the context of the multi-objective optimization of wind turbine blades can be found in [19][20][21][22][23][24][25][26].
This paper proposes the formulation of MOSOPs considering other objective functions that go beyond the traditional ones found in the literature, resulting in several approaches.In this sense, essential aspects, such as the dynamic behavior of the blade in the context of natural frequencies of vibration, avoiding resonance, and its global stability regarding the maximum loads that can cause some buckling, are considered in the multi-objective optimization problems proposed in this paper.The weight minimization can generate flexible structures with very low natural frequencies of vibration (less than 1 Hz).In this way, these structures are susceptible to the dynamic action of the wind.It can be seen, therefore, that minimizing weight and maximizing natural frequencies to avoid resonance are conflicting objectives but important for the design.In addition, maximizing the other two (second and third) can also prevent the possibility of undesirable vibration modes superposition.
Concerning aspects related to buckling of the structure, it can be necessary, or even mandatory, to investigate the its global stability more precisely [27][28][29].In Richardson et al. [27] the authors wrote in this context: "This is an important constraint in practice, which would make the method more relevant for practical application.".
It is important to remark that minimizing weight, maximizing the natural frequencies of vibration, and maximizing the first critical load factor can offer the decision-maker (DM) a range of possibilities to choose optimized blade designs for these objective functions, avoiding uncountable single-objective optimizations with various combinations of these.Thus, the main contributions of this paper are summarized as: (i) the consideration of natural frequencies of vibration and critical load factors concerning global stability as objective functions, generally neglected; (ii) the consideration of natural frequencies of vibration and critical load factors concerning global stability as constraints when these are not objective functions, generally neglected; (iii) the proposition of structural optimization problems with two, three and four objectives; and, (iv) the extraction of desired nondominated solutions according to the preferences of the DM.The formulations of the MOSOPs proposed in this paper have not yet been covered by studies found in the literature.In this sense, it is understood that this is the main objective of the discussion carried out here.
It is essential to emphasize that the formulations of the optimization problems with multiple objectives present a series of advantages compared to the formulations with only one objective.Pareto fronts (PFs) provide a set of non-dominated solutions obtained by avoiding formulating and solving a series of single-objective optimization problems independently, leading to a high computational cost.It is important to note that, although the main idea is to consider conflicting objective functions, one can also include "redundant" objective functions in formulating multi-objective problems.As a result, PFs provide the decision-maker with a set of solutions from which he can make his choices according to his criteria and preferences.
In the proposed MOSOPs, the design variables are the layups on different regions of the blade.The constraints are the maximum strain, the blade tip displacement, the natural frequencies of vibrations, and the load factor concerning the global stability of the blade when these are not considered objective functions.The blade is discretized in a finite element parameterized model using the Femap software, available at https://www.plm.automation.siemens.com/global/en/products/simcenter/femap.html(accessed on 22 September 2022), according to geometric design variables.The structural analysis of a wind turbine blade requires many loading conditions.Thus, the analysis of blades requires expensive and time-consuming processing, as the models are built using finite element meshes that are properly discretized.Therefore, the blade model is subjected to the loads determined by the actuator line model in a critical load condition [30] to reduce the excessive time demanded by these analyses.Among various multi-objective evolutionary algorithms MOEAs from the evolutionary computation literature such as MOEA/D [31], MOPSO [32] MOGWO [33], the non-dominated sorting genetic algorithm II (NSGA-II) [34] is the evolutionary algorithm, largely used in the literature, adopted to solve MOSOPs discussed in this paper.
The remainder of the paper is organized as follows: The formulations of the proposed multi-objective optimization problems, the optimization framework, the design variables, and the design constraints are provided in Section 2. The geometric model of the blade, detailing the finite element model and the load conditions, is presented in Section 3. The evolutionary algorithm NSGA-II adopted in this paper to solve the proposed multi-objective optimization problems, the performance indicators, and a multi-criteria decision-making (MCDM) used to extract solutions from the PF are provided in Section 4. The experiments, results, and discussions are presented in Section 5. Finally, the conclusions and extensions of this work are reported in Section 6.

The Multi-Objective Structural Optimization Problem
The basis that was used to propose and formulate the MOSOPs of the wind turbine blade, analyzed in this paper, was inspired by aspects related to its dynamic behavior and structural stability, such as the natural frequencies of vibration and load factors concerning the global stability of the blade.

MOSOPs Formulations
Four MOSOPs were defined and analyzed, with different objective functions and constraints.The blade's weight (W(x)), the tip displacement (δ max (x)), the load factor related to the critical buckling load (λ(x)), and the natural frequencies of vibration (ω i , i = 1, 2, 3), are considered as objectives (varying according to each MOSOP).Their formulations are described following.

3.
The MOSOP3 is given by: The MOSOP4 is given by: where ε(x) is the maximum failure strain.The blade tip displacement, load factor, and strain are limited by δ, ε, and λ, respectively.The limits of the design variables is defined by the lower x and upper x bounds.

The Optimization Framework
The blade's finite element model was developed in the Simcenter Femap software and NX Nastran solver [35].The optimization problem was implemented in Python [36].The algorithm comprises two parts, the first one being the optimization itself and the second making the link between the optimization code and the finite element model, updating the blade's structure properties, running the simulations, and collecting the results.The first part uses the pymoo library [37].The second part of the code is implemented using the Femap's application programming interface (API) [38].The optimization flowchart combining NSGA-II algorithm with finite element method is provided in Figure 1.

Design Variables
The optimization problem proposed has the design variables related to the laminate plies (number of plies, material, and fiber orientation) in each of the 30 blade regions detailed in Section 3. As material options for the plies, carbon fiber, fiberglass, and PVC foam (core material commercially named Divinycell) were considered, and the mechanical properties used are presented in Table 1.The thickness of the plies of each material was defined based on the mechanical testing studies used in [39,40].Carbon fiber plies have 2 mm thickness, fiberglass plies have 1 mm thickness, and Divinycell plies have 5 mm thickness.The angles considered for the fiberglass and carbon fiber plies were 0 • , 45 • , −45 • , and 90 • , and for the Divinycell ply, only 0 • was taken into account since it is an isotropic material.It is important to remark that elements that influence the mechanical properties of GFRP composites include fiber type, matrix, relative amounts of constituents, and fiber orientation [41][42][43].In Table 1 X ε t and X ε c are, respectively, the allowable tensile and compressive strains on the longitudinal direction; Y ε t and Y ε c are, respectively, the allowable tensile and compressive strain on the transverse direction; S ε is the allowable strain for in-plane shear; and ε 1 , ε 2 , and γ 12 are the longitudinal, transverse, and shear elastic strains, respectively.The values of X ε t , X ε c , Y ε t , Y ε c , and S ε of all materials are considered in the design.

Property
Carbon Fiber UD [39] Fiberglass UD [39] Divinycell H200 [40] ρ Some standard industrial design rules were adopted [44,45] to guarantee the stack-up order manufacturing, such as the following: The R3 rule is satisfied simply by studying half the laminate and mirroring it.The stackup order would be fixed to guarantee that all rules are respected and, at the same time, maintain a reasonable number of design variables, and the order defined is shown in Figure 2. First, to describe the design variable, consider only one laminate instead of the entire problem.Then, considering the stack-up sequence defined above, a total of nine types of plies are available for the design.Material and orientation of the plies are constants, but how many plies each ply type will have can be selected: N 1 , N 2 , ..., N 9 .However, to guarantee that rule R2 it's respected, the number of plies at +45 • and −45 • must be the same, reducing the number of design variables to seven: N 1 , N 2 , ..., N 7 .
Second, considering the complete optimization problem, two main laminates are defined, one laminate for the airfoil skin of the blade and the other for the spars.Therefore, those two laminates can be defined with seven design variables each, that can be written as N i A and N i S , for 1 ≤ i ≤ 7, referring respectively to the airfoil skin laminate and the spar laminate.
With those two main laminates defined, the blade's 30 regions, described in Section 3, are populated with laminates that are proportional to the main laminates.The number one, the i-th ply type on the j-th region laminate, assuming it's in the blade's spar, can be defined as , where β Rj is the proportion factor and 1 ≤ j ≤ 30.Finally, we can define the design variables array x, and their lower and upper bounds as

Design Constraints
The constraints are normalized, such as where m j is the number of load factors of the blade and γ m is the partial safety factor for materials, defined in [46] as 1.2 for stability analysis (Equation ( 9)) and 1.3 for rupture analysis of non ductile materials, such as composite (Equation ( 8)).
After solving the equilibrium equation obtained from a discrete system of finite elements, the nodal displacements {δ} are achieved, given by: where [K] is the stiffness matrix and {p} are the load components [47].
The solution of the eigenvectors and eigenvalues of the Equation ( 11) provides the natural frequencies of vibrations of the blade.
where [M] is the mass matrix, and φ m f is the m f -th eigenvector corresponding to the m f -th eigenvalue [47].
The load factors λ are obtained by the evaluation of the eigenvalues of the matrix where [K G ] is the geometric matrix of the structure, and λ m λ are the equivalent eigenvalues with respect to the m λ load factors of the structure.The lowest value λ cr of λ m λ gives the buckling load factor or critical load factor of the structure.The eigenvector ν m λ represents the corresponding instability mode of the blade [48].
The laminate failure theory adopted is the maximum strain criteria proposed by Waddoups [49].This theory limits the strain on the material's longitudinal, transverse, and shear directions and allows the analyst to identify the failure mode.The maximum strain failure index I f f is calculated as:

Geometric Model of the Blade
The blade structure is defined by a constant outer skin that gives the blade its aerodynamic airfoil shape and two internal spars positioned at 15% and 50% of the chord [50].The blade contains 17 aerodynamic profiles along its length, including the Cylinder1, Cylinder2, Cylinder3, DUO40_A17, DUO35_A17, DUO30_A17, DUO25_A17, DUO21_A17, and NACA64_A17.The blade's total length is 61.5 m, and the whole mass is 17,740 kg.The completed distributed blade aerodynamic properties can be found in [51].The blade is divided into thirty regions (Figures 3 and 4), the outer skin is divided into three parts along its chord, and each of the two spars completes the five different layups in each airfoil section, and span-wise, the blade is divided into six sections.Thus, in each region, a different composite layup can be defined.

Finite Element Model
The blade's structure was modeled using linear shell elements.A mesh convergence study was performed to determine the adequate element size.The element size was reduced while the value of the first natural vibration frequency (ω 1 ) was monitored, as shown in Figure 5.The final element size set was 300 mm, resulting in a mesh with 7610 elements and 35,033 nodes.From those elements, 7509, or 98.67%, are QUAD4 elements (three components of translation and two in-plane components of rotation at corners); a 2 × 2 set of integration points is used for membrane and bending strains in order to accommodate coupling, and 101, or 1.33%, are TRIA3 elements (combination of a constant strain triangular (CST) membrane with a composite bending triangle).Both finite elements are completely detailed in [52].Three load components and the torsional moment (M Z ) were applied to the structure at nine different locations along its length.In addition, an interpolation element (RBE3) was created to apply the loads on the structure for each load position.Its central node was located exactly on the aerodynamic center of the airfoil, as illustrated by Figure 6.
The structure was constrained at its root, where the blade is connected to the wind turbine hub.Every node on this interface region, highlighted in red in Figure 7 had its six degrees of freedom constrained.
After the finite element model it's defined the matrices [K], [M], and [K G ] can be assembled after the transformation to a global axis with the element rotation matrix R i for the i-th element, where [K i ], [M i ], and K G i are the elemental stiffness, mass, and geometric matrices, respectively.The components of strains in x, y, and xy directions concerning any point throughout the laminate thickness of each finite element are written as where ε are the strain components at an arbitrary point on the laminate thickness, ε o are the element's membrane strain components, κ are the surface curvature components of the element, z is the distance between the analyzed point of the laminate to its mid plane, δ e are the element displacements, and B e f is the bending-membrane strain matrix of the element.The xz and yz strain components can be obtained similarly.

Load Conditions
The loads used in the analysis are calculated according to the regulatory standard IEC 61400-1 [46] and using the aerodynamic module AeroDyn developed by NREL and embodied in FASTV8 packages developed by NREL [53].
AeroDyn uses the actuator line technique to calculate the aerodynamic loads on wind turbine blades.The program approximates the 3D flow around the blade by computing the 2D flow at cross-sections along the blade.The lift forces, pitching moments, and drag forces at each cross-section are then computed and lumped at a node in the 2D cross-section.
The 2D forces and moment at each node are distributed throughout each blade as distributed loads per unit length.Finally, the 2D distributed loads are integrated along the blade's length to obtain the total 3D aerodynamic loads on the blade.This method allows for a computationally efficient calculation of the aerodynamic loads on wind turbine blades, which is essential for the design and optimization of wind turbines [30].
The software only provides output for nine stations along the span of the blade.In each finite element analysis (FEA), the three components of internal forces and the torsional moment (M Z ) calculated were applied at each airfoil's aerodynamic center, which was assumed to be at 25% of the chord.
From all load conditions defined in [46], the most critical (using the tip deflection as a criterion) was the Design Load Case 1.3 (DLC 1.3), that is, at energy production situation on extreme turbulence model (ETM) and with the wind velocity at the cut-out wind speed (V out = 25 m/s).In addition, a yaw misalignment of 8 o was also considered as a deviation from optimal operational conditions.Therefore, that was the one used in the optimization process.Figure 8 shows the concentrated loads calculated using FASTv8 on the critical condition.The international standard [46] takes into account uncertainties and variability of loads by partial safety factors, to ensure the safety of the design, defined as: where F d is the design value of the internal loads, γ f is the partial safety factors for loads, and F k is the characteristic value for loads.The condition analyzed (DLC 1.3) is classified by the standard as a normal condition.Thus, the partial safety factor for loads (γ f ) is 1.35.

NSGA-II and Performance Indicators
NSGA-II is a multi-objective evolutionary algorithm that is widely used to solve MOOPs.It was proposed by Deb and Deb et al. in their paper [34], and it features three important components: (i) Non-dominated sorting approach: This is a fast approach used to identify the non-dominated solutions in a population.It assigns a rank to each solution based on how many other solutions dominate it.Solutions with a lower rank are better than solutions with a higher rank.(ii) Crowding distance estimation procedure: This measures the diversity of the solutions in a population.It calculates the distance between neighboring solutions in the objective space and uses this to determine which solutions should be preserved in the next generation.(iii) Crowding comparison operator: This is a simple operator used to compare solutions based on their rank and crowding distance.It selects the best solutions for the next generation based on these criteria.
By combining these three components, NSGA-II is able to efficiently search for a set of non-dominated solutions that represent a diverse range of trade-offs between the objectives.It has been shown to be effective in solving a wide range of multi-objective optimization problems and is widely used in research and industry.

Performance Indicators
Several performance indicators have been proposed to evaluate the quality of nondominated solutions on the PFs.Although the true PF of optimal solutions is unknown in most real engineering application problems, metrics to evaluate and compare the quality of the approximations obtained from different algorithms are critical to assessing the results.
Cardinality, convergence, distribution, and spread are some of the different characteristics measured by those indicators [54].In this paper, the hypervolume (HV) and spacing (S) metrics are adopted for all MOSOPs studied since both are quite popular in the literature and do not require the previous knowledge of the true PF.

Hypervolume
The HV indicator proposed by [55] can be described as the volume dominated by the PF approximation in the objective space.To calculate the HV, it is necessary to set a reference point (R HV ) in the objective functions space.The reference points selected for each MOSOP in this study are presented in Table 2.The HV helps to evaluate both the convergence and the distribution (or diversity) of the non-dominated solution throughout the PF.

Spacing
As the name indicates, this indicator measures how the non-dominated solutions on the PF spread.Proposed by [56], the spacing metric is calculated as: where SP is the spacing indicator value, S is the set of solutions that compose the PF approximation obtained, d i is the distance between a point s i S and the closest point of the PF approximation, d is the mean value of d i calculated for all solutions in S, and NOb is the number of objectives in the optimization problem.

Multi-Criteria Decision Making (MCDM)
It is of great interest to structural designers to extract desired solutions from the PF.This extraction must provide the values of the objective functions and the final design variables.On the other hand, this task is not trivial, and just a visual inspection of the PF maybe not be sufficient.
Several MCDM are provided in the literature to extract solutions from PFs.Among them, a technique for order of preference by similarity to ideal solution TOPSIS is a robust technique proposed by Hwang and Yoon [57] largely used in the literature, such as in [58][59][60][61][62], to find the best compromise solutions from the PF.Using pseudo-weigth vector method is also one alternative for a MCDM and one of them was adopted in this paper.
The pseudo-weight vector approach proposed in [63] gives a simple and practical way to extract the desired solution from a PF and perform a visual analysis and comparison between trade-off solutions: where f min i and f max i denote the minimum and maximum values of the i-th objective function among the solutions of the PF.This method generates, for each solution of the PF, a pseudo-weight vector w that represents the relative importance factor for each objective function.According to the DM weight preference (w p ) the solution with the closest pseudoweight vector is selected.
Instead of a strict value of w p , a good practice for the actual DM is to select ranges of values and further explore the different solutions extracted by the MCDM.However, for the sake of this study, a single preferred weight vector was selected for each MOSOP and presented in Table 2.

Experiments, Results and Discussion
After all independent runs, the solutions were combined, identifying the non-dominated solutions leading to the PF.Finally, performance indicators were obtained, and the MCDM was used to extract the preferred solution.
The parallel coordinates [64] were used to represent the normalized values of the objective functions in the range [0, 1], for the MOSOP3 and MOSOP4, both presenting four objective functions.The horizontal axis represents each objective function and the vertical axis represents its respective normalized value.
Figure 9 shows the non-dominated solutions of the PF of the MOSOP1 in which the extracted solutions is marked in red. Figure 10 details the schematic representation of the laminate selected for the extracted solution.In the similar figures of all MOSOPs, the regions were grouped length wise to the blade, and the regions labels refer to Figures 3 and 4. Figure 11 shows the FEA fiber failure index (I f f ) and buckling eigenvector obtained for the extracted solution of MOSOP1.Figure 12 depicts several views of the PF obtained for MOSOP2, Figure 13 shows a 3D visualization of the PF, Figure 14 provides the schematic representation of the laminate of the extracted solution, and Figure 15 details the FEA fiber failure index (I f f ), and the buckling eigenvector results of the extracted solution of MOSOP2.The parallel coordinates of the PF obtained for MOSOP3 are plotted in Figure 16.The pair-plots of the PF obtained for MOSOP3 are depicted in Figure 17. Figure 18 details the schematic representation of the laminate selected of the extracted solution, Figure 16 shows the PF represented in parallel coordinates [64], and Figure 19 provides the FEA fiber failure index (I f f ), and buckling eigenvector obtained for the extracted solution of MOSOP3.Figure 20 shows the PF represented in parallel coordinates, the pair-plots of the PF obtained for MOSOP4 are depicted in Figure 21, Figure 22 details the schematic representation of the laminate selected for the extracted solution, Figure 23 provides the FEA fiber failure index (I f f ) and buckling eigenvector results of the extracted solution, and Figure 24 illustrates the first three vibration modes of the extracted solution.The failure index (I f f ) on the linear static analysis and the displacement on the buckling results (the eigenvector), depicted in Figures 11,15,19 and 23, has four contours, from top to bottom: the first is the I f f in the entire model; the second is the same as the first but with the skin hidden so that the spars can be seen; the third is the displacement in the buckling analysis for the entire model; and the fourth is the same as the third but with skin elements hidden so that the spars can be seen.The color scale in those images refers to the I f f values.
The performance indicators calculated values, the reference point used in the HV computation (R HV ), and the preferred weight vector (w p ) used in the MCDM are displayed in Table 2. Analyzing the SP metric, one can notice that MOSOP3 had the most spread PF approximation, and MOSOP1 and MOSOP4 had considerably less distribution across the objective space.It's worth noting that both MOSOP2 and MOSOP44 are the only cases with δ max (x) as constraints, which could be causing the smaller SP.The HV, however, can not be compared between MOSOPS since the order of magnitude of each objective influences the indicator's final value.Observing the extracted solutions of each PF concerning all MOSOPS, the blade mass (W(x)) is quite similar, approximately 20 t.
When analyzing the PFs obtained, it can be noticed that across different MOSOPS, the lower bound (or maximum bound if the objective function was to be maximized) of values for the same objective function were all very similar, if not the same.Taking the W(x) objective as an example, every PF includes at least one solutions close to the W(x) = 20 t.That can indicate that the solutions are converging to the actual PF and the algorithm's robustness.
Through the pair-plots, it's possible to note correlations between the multiple objectives, providing insights for new design exploration and new optimization techniques to be implemented and serving as a verification of the model consistency to the expectations set by the knowledge of the physical behavior of the problem.Figures 12 and 17 demonstrate a clear non linear correlation between the increase in weight and the decrease in the maximum displacement.In Figure 17, a relation between the increases in W(x) and λ 1 (x) can be observed.In Figure 21 it is possible to see not so clear linear relations between all the natural frequencies and the blade's weight as well as between each other.
When analyzing the laminates selected (Figures 10, 14, 18 and 22), it is important to remark that all plies were plotted with the same thickness to facilitate the visualization.Actually, different materials have different thicknesses (carbon fiber plies are 2 mm, fiberglass are 1 mm, and core plies are 5 mm in thickness).Thus, a few observations can be made, such as: 1.
In the skin regions, the laminates use PVC foam cores, while in the spar regions, no solution has a core.That is because the moment of inertia of the blade in the direction of the main bending loads is heavily increased as the skin thickness gets higher, but is not influenced by the spars' thickness.

2.
In all laminates, the overall thickness of the carbon fiber plies is at least double the thickness of the fiberglass.Considering that the fiberglass is heavier and provides less stiffness, that was expected.The x L values set didn't allow the number of fiberglass plies to go down to zero.Thus, the fiberglass to carbon fiber ratio could be even lower.However, the blade's cost is an important industrial design parameter that was not considered in the optimization problem definition.Fiberglass woven is cheaper than carbon fiber and is the main material used in wind turbine blades.

3.
All solutions observed an unexpected increase in the thickness of the spar regions close to the tip of the blade.The change in the aerodynamic airfoils can explain that behavior through the length of the blade that causes the skin of the blade to be flatter and therefore more prone to local buckling, requiring the spars' reinforcement.However, it's also important to add that due to the blade's tapering, the increase in thickness in those areas causes a minimal increase in weight.

4.
Through the analysis of results, it can be concluded that adding ply staggering strategies to the laminate parameterization to avoid abrupt changes in stiffness could also provide better results with a better load path through the structure.
When analyzing the FEA results, an important observation that can be reported is the strain concentration points.It is clear that the change in stiffness in the transition between laminate regions generates high values of strain and stresses, as could be expected with significant laminate drop-offs.It was also possible to visualize that, except for MOSOP4, all other solutions had extremely localized instability events that also occurred in the laminate regions' transitions.
The design attributes ω 1 (x), ω 2 (x), ω 3 (x), δ max (x), λ 1 (x), and I f f max (x), of all extracted solutions through MCDM are presented in Table 3.The radar plot provided in Figure 25 was used to compare the solutions of the selected final blade's geometry of each MOSOP.It allowed us to quickly compare the different design aspects and perform visual tradeoff considerations.

Conclusions and Extensions
Structural optimization of wind turbine blades is proposed in this paper, considering four different sets of conflicting objective functions, such as the minimization of weight, maximization of the first three natural frequencies of vibration, minimization of the blade tip displacement, and maximization of the load factor referent to the first buckling mode of the structure.
The MOSOPs discussed in this paper have included several objective functions in the proposed formulations that are traditionally neglected.Those objectives were combined in two, three, and four MOSOPs.The constraints were the fiber failure index, considering the maximum principal strain as failure criteria, besides one or more of the objectives when they were not considered objective functions.
The objective and constraint functions were evaluated through FEA in the academic version of Femap.The optimization algorithm used was the evolutionary multi-objective NSGA-II [65].The HV and S performance indicators were computed and made available for comparison between MOSOPs.The MCDM [63] pseudo-weight method was used to extract the solutions from the PF optimization.For every MOSOP defined, the PF approximation was presented, including the solutions selected through MCDM, their laminate schematic representation, the fiber failure index, and buckling contour plots.
The results obtained with this study were considered competitive compared to currently available wind turbine blades, specifically the 5 MW NREL blade used as a reference.
In this sense, aspects related to dynamic behavior and global stability, important issues in the structural design of a wind turbine blade, were highlighted.Furthermore, the various proposed formulations provided sets of non-dominated solutions, which did not require the solution of several sequential single-objective structural optimization problems.In addition, these non-dominated solutions were available for the DM to make choices according to given preferences.Also, the formulations of these MOSOPs proposed in this paper have not yet been covered by studies in the literature.In this sense, this was the main objective of the discussion carried out in this paper.
The observations and conclusions carried out from the present work allow its continuity.Firstly, high-fidelity structural analysis comes with performance challenges.For example, each independent run of the optimization had a high computational time (16-18 h).To overcome this challenge, a future research extension will consider surrogate models, proper orthogonal decomposition (POD), and radial basis functions (RBF) [66] approximation while comparing its results to high-fidelity solutions.Large displacements, nonlinearity, fatigue analysis, and dynamic load cases are some aspects not considered for this study, but in the future, they will be incorporated into the model.Important aspects concerning the vulnerability of wind turbine blades to various damages and failures will be considered by evaluating their performance, especially after optimization [67,68].Finally, the inclusion of other objective functions in MOSOPs formulations will be considered in future studies , the use of MOEAs specially designed to solve many-objective optimization problems such as NSGA-III [69], and other MCDMs such as TOPSIS.

Figure 1 .
Figure 1.The optimization flowchart combining NSGA-II algorithm with finite element method.
• (R1) Disorientation rule: The maximum angle offset between two consecutive plies is 45 • .• (R2) Laminate balancing rule: The laminate must be balanced with respect to the 0 • direction, which implies the same number of 45 • and −45 • plies.• (R3) Symmetrical laminate rule: to avoid the in-plane traction and bending coupling.

Figure 2 .
Figure 2. Fixed stack-up sequence on the optimization problem.

Figure 3 .
Figure 3. Blade's outer skin regions divided in 18 regions containing a different composite layup.

Figure 4 .
Figure 4. Blade's spars regions is divided in six, from 19-30, symmetric regions containing a different composite layup.

Figure 5 .
Figure 5. Mesh convergence study to determine the element size.

Figure 6 .
Figure 6.Load application to the structure using RBE3 element.

Figure 7 .
Figure 7. Fixed constraint applied to the root of the structure.

Figure 8 .
Figure 8. Concentrated loads calculated using FASTv8 on the critical condition.

Figure 10 .
Figure 10.Schematic representation of the laminate selected through MCDM from the MOSOP1 PF.

Figure 11 .
Figure 11.FEA fiber failure index and buckling eigenvector results for MOSOP1 solution extracted.

Figure 12 .
Figure 12.Scatter plot of the PF obtained for MOSOP2.

Figure 13 .
Figure 13.3D Visualization of the PF obtained for MOSOP2.

Figure 14 .
Figure 14.Schematic representation of the laminate selected through MCDM from the MOSOP2 PF.

Figure 15 .
Figure 15.FEA fiber failure index and buckling eigenvector results for MOSOP2 solution extracted.

Figure 16 .
Figure 16.Parallel coordinates plot of the PF obtained for MOSOP3.

Figure 17 .
Figure 17.Scatter plot of the PF obtained for MOSOP3.

Figure 18 .
Figure 18.Schematic representation of the laminate selected through MCDM from the MOSOP3 PF.

Figure 19 .
Figure 19.FEA fiber failure index and buckling eigenvector results for MOSOP3 solution extracted.

Figure 20 .
Figure 20.Parallel Coordinates plot of the PF obtained for MOSOP4.

Figure 21 .
Figure 21.Scatter plot of the PF obtained for MOSOP4.

Figure 22 .
Figure 22.Schematic representation of the laminate selected through MCDM from the MOSOP4 PF.

Figure 23 .
Figure 23.FEA fiber failure index and buckling eigenvector results for MOSOP4 solution extracted.

Figure 24 .
Figure 24.FEA vibration modes eigenvectors of the extracted solution of MOSOP4.

Figure 25 .
Figure 25.Radar plot comparison of normalized objective functions.

Table 2 .
Parameters used in the optimization results post-processing and performance indicators.

Table 3 .
Design attributes of each extracted solution.