Multi-material and thickness optimization of a wind turbine blade root section

State-of-the-art wind turbines blades are large and complex structures, consequently making design challenging. The complexity is largely a result of us-ing laminated composites which are high-performing and versatile materials, but also inhomogeneous hence exhibiting anisotropic behavior as well as challenging failure mechanisms. Structural optimization has been shown to be an invaluable tool for solving such challenging design problems, however a framework for treating all relevant parameters, i.e., material choice, ﬁber orientation, stacking sequence, and layer thickness, on the current scale is not available. Thus, the motivation of the present work is to present an over-all gradient-based approach, where the key structural criteria are considered for multiple design load cases. The optimization framework is based on a variation of the Discrete Material and Thickness Optimization approach, where the thickness is directly parametrized, allowing for appropriately treating the sandwich parts of the blade. Application of the developed framework leads to an optimized design consisting of complex variable-thickness laminates, a good overall distribution of the structural criteria in the model, and a 23% reduction in mass compared to the initial design.


Introduction
Recent years have seen a significant increase in investment in renewable energy solutions (IEA, 2023).Wind energy has become particularly desirable due achieved reductions in Levelized Cost of Energy (LCoE), versatility, as turbines can be established almost anywhere, and general green profile of the technology (GWEC, 2023).Nevertheless, it is desirable to further reduced costs to enhance the investment appeal and provide more energy to satisfy the increasing demand of the populace.A key metric for maximizing energy harvesting is rotor span size, particularly for offshore turbines that have higher installation costs, and thus, increasing the size cost-efficiently is the main objective in design.As a result, the length of state-of-theart commercial offshore blades currently exceed 100 m, and this is expected to continue to increase in the future.Design of such blades is challenging as a consequence of the scale, complex shape, and use of highperforming laminated composite materials.The particular laminated composites utilized consist of fiberreinforced polymers and transversely-stiff low-density core materials such as balsa wood or structural foams.Because of the inhomogeneous nature of laminated composites, they exhibit anisotropic behavior and have several associated failure mechanisms which can make failure prediction particularly tedious, consequently slowing the design process.
Accelerating design can be achieved through development and application of structural optimization methods.Several attempts at structural optimization of wind turbine blades have been demonstrated in literature.One of the pioneering works by Jureczko et al. (2005) demonstrated full blade optimization under multiple structural criteria through the applica-tion of a zero-order optimization method.Numerous following works have also proposed a zero-order approach to solving the problem with more criteria and analysis-model fidelity as computational power has increased, see Chen et al. (2013); Monte et al. (2013); Barnes and Morozov (2016); Monte et al. (2017); Albanesi et al. (2018a,b).In these works, structures are parametrized with material choice, layer, or laminate thickness.
Although correct application of zero-order methods can lead to finding the global optimum of a design problem, the computational cost associated severely restricts the treatable size of the problem to, typically, less than a dozen design variables.An alternative that allows for solving optimization problems several orders of magnitude larger is a gradient-based approach.Early works on gradient-based optimization of wind turbine blades have typically utilized a multi-step optimization approach, initially applying topology optimization to design the internal configuration of the considered blades, followed by sizing optimization of the outer shell (Buckney et al., 2013;Forcier and Joncas, 2012;Albanesi et al., 2020).Other works focus exclusively on the sizing optimization part to allow for more fidelity in the analysis model.Such sizing optimization has been demonstrated on shell models in Bottasso et al. (2014); Hermansen et al. (2023) and solid-shell models in Sjølund and Lund (2018), all including numerous structural criteria directly in the gradient-based optimization.Noteworthy is also Hu et al. (2016), where authors demonstrate a reliabilitybased optimization approach applied to a blade to optimize for random fatigue loads.
Design of wind turbine blades is inherently a multidisciplinary undertaking with mutual dependence between aerodynamic and structural criteria.It is therefore relevant to attempt including both disciplines in optimization to minimize post-optimization modifications and several works have shown the potential of carrying this out.Such optimization has been performed using either a coupled (Bottasso et al., 2012;Bortolotti et al., 2016;Scott et al., 2020) or de-coupled (Bottasso et al., 2016;Scott et al., 2022;Mangano et al., 2022) workflow.It should be noted, that including aerodynamic considerations in the optimization problem increases the computational cost associated with analysis and computation of sensitivities.A reduction in fidelity is therefore necessary to be able to solve the problem compared to separate structural and aeroelastic optimizations.The structural part of the multidisciplinary optimizations is typically carried out as sizing optimization with layer thickness as design variables.
The state of the art review indicates that a general framework for treating material selection, fiber orientation, stacking sequence, and thickness simultaneously is not available for gradient-based wind turbine blade optimization.This work hence aims to extend the state of the art by development and application of such gradient-based multi-material and thickness optimization to wind turbine blades.Such problems are computationally challenging to solve due to the inherently many design variables, that represent the material candidates and thickness of every layer in the model.To manage the computational expense, optimization is carried out on the outer shell of the root section of a state-of-the-art offshore commercial blade, spanning the first 35m of the particular blade.Focusing on the root section is motivated by the majority of material being placed here, i.e., 60% of the total blade mass.Moreover, the region is subject to the largest bending loads and its shape is primarily governed by the root connection and not aerodynamic considerations.As such, the costly aerodynamic analysis and sensitivity analysis as well as coupling considerations to the structural analysis are neglected as part of this work.
Gradient-based multi-material and thickness optimization is enabled through the Discrete Material and Thickness Optimization (DMTO) approach (Stegmann and Lund, 2005;Lund and Stegmann, 2005;Sørensen and Lund, 2013;Sørensen et al., 2014), where the inherently discrete problem is made continuous and intermediate values are penalized.As opposed to the original DMTO formulation, a direct parametrization of the thickness (Sjølund et al., 2018) is used in this work to appropriately treat the sandwich laminates in the design.The relevant failure criteria, i.e., buckling, static and fatigue failure, are all considered, each for four design load cases.
The paper will continue with a presentation of the wind turbine blade root section model in Section 2 which includes the finite element strategy and how loads are applied.In Section 3 the Discrete Material and Thickness Optimization approach appropriate for parametrizing both monolithic and sandwich regions of the blade is introduced.The structural analyses for assessing the critical criteria are introduced in the following Section 4 along with the associated safety factors that are dictated by design standards.In Section 5 the optimization methods and strategy are presented with emphasis on how the problem can be solved efficiently.Results are then shown in Section 6 and are followed by a discussion in Section 7 of the overall approach adopted.The paper is then finalized with conclusions in Section 8.

Modeling the Blade
The optimization is to be carried out on the outer shell of a root section, i.e., the first 35 m of a state-of-theart industrial offshore wind turbine blade with a total length exceeding 100 m.Despite only considering the root section, the structural analysis and subsequent sensitivity analysis is still computationally challenging.Focusing on only the outer shell of the root section reduces the computational expense, and as this part of the blade contains the majority of material in the blade, it is a region of considerable potential for optimization.Moreover, the increase in size creates a stronger coupling between aerodynamic and structural effects and thus applying structural optimization on a full blade would necessitate more postprocessing for the aerodynamic design.The shape of the blade at the root is however governed mostly by the connection and the dependence on aerodynamics is therefore less pronounced.
The work is based on the model presented in Hermansen et al. (2023).However, as this work considers both material selection and thickness optimization, the blade is parametrized with a generic layup to demonstrate the potential of the approach for application during early design phases.Hence, the main work in parametrizing the structure is division into a number of element groups (patches), where the generic layup is associated.As outlined in the previous work, the partitioning of the blade is carried out in two steps.The first step is cross-sectional partitioning which is done based on a number of identified characteristic regions.These regions are shown in Fig. 1.Fig. 1: A cross section in the considered blade root section with characteristic regions indicated.Note that although the patches are split transversely for the trailing edge core region, these two regions share the same material design variable.
Material design variables are assigned to these regions, thus spanning the entire length of the blade in or-der to comply with the manufacturing, where plies are draped from a roll, see e.g.Krogh et al. (2023).The subsequent partitioning of the length of the blade creates a number of smaller patches, where the thickness design variables are associated, hence allowing the optimization to yield a variable-thickness design, see Fig. 2.Many materials are used to manufacture the real blade including different Uni-Directional (UD) Glass and Carbon Fiber Reinforced Polymers (GFRP and CFRP), oriented in various directions, as well as balsa wood and structural foams.The number of materials are simplified in the optimization to GFRP in UD, 0 • /90 • Cross-Ply (CP) and ±45 • biax for fiber candidates while PVC foam is used for the core material candidate.Note that balsa wood is typically used for the core material in the outer shell of blades with a lower-density PVC foam for the shear webs.Due to the lack of fatigue properties for both, Rohacell WF110, which is a higher-density closed-cell PMI foam, is used as a compromise for both.The material properties are summarized in tables 1 and 2 for static and fatigue respectively.OPTIDAT Nijssen et al. (2006) database.Thus, the applied properties are only representative of these particular materials and can vary depending on the supplier.Unfortunately, no data has been found for the fatigue properties of cross-ply laminates for all directions and R-values.The fatigue properties are therefore estimated from a scaling between the static strength values and available fatigue data.The foam fatigue  properties are based on that found in Zenkert andBurman (2009, 2011).It follows that the candidate materials are assigned differently for monolithic regions (spar cap and edges) and for the core regions.Monolithic areas are thus assigned only fiber candidate materials and the core regions are assigned a sandwich-type laminate, where the center layers have only core material candidates, and the face layers have fiber material candidates.For generality, the same initial layup has 20 layers in every patch.The fiber regions use a uniform layer thick-ness of 4mm, while the core regions use 2.5mm for the face layers and 6.25mm for the core layers.Note that actual fiber mats are substantially thinner, however modeling these with their "as-manufactured" thickness would significantly increase computational cost.Using thicker layers does imply reduced accuracy in the failure and fatigue analyses, where evaluation is done at the top and bottom of every layer, however this is deemed acceptable during preliminary design.
A shell-based finite element model is used in the analysis.Shell models are standard in the industry since a shell approach is particularly convenient for modeling blades.Using a layup-offset formulation, only the outer shape has to be represented to conduct analysis with layups offset inwards.Some discrepancies have been highlighted in literature regarding shell modeling, namely incorrect representation of the crosssectional area of closed and highly curved structures (Laird et al., 2005;Branner et al., 2007), and the influence of the artificial drilling degree of freedom on the torsional stiffness estimation (Tavares et al., 2022).These effects can be minimized through careful consideration when establishing the model.
The fidelity of the model used is selected based on a trade-off between accuracy and numerical efficiency.To verify the model, its strain distribution has been compared to that of a high-fidelity model that is used for certification of the blade considered.The model is meshed with 9-noded isoparametric shell finite elements.The mesh consists of 26,569 of such elements and 105,239 corresponding nodes, see Fig. 3.

Load Application
Analysis of wind turbine blades is a multi-disciplinary undertaking.Besides the structural considerations, aerodynamic analyses and syntheses are used to find an efficient shape to maximize the power output for given wind conditions.Such analysis hence yields the loads that are applied to the blade for the structural analysis.In this work, a load envelope with four load cases resulting from such analysis is used.The load cases  The criteria for assessing the structural integrity have to be evaluated for each of these load cases.This motivates the choice of using only four load cases in the optimization, which is the minimum prescribed by the design standards (DNV-GL, 2015) to minimize the computational cost.The trade-off is an increase in safety factors applied to each criterion.These are determined in Subsection 4.4 after introduction of the structural criteria.
The loads are given as a distribution of moment loads along the length of the blade.A spline is then fit to these distributions and from these spline fits, a distributed load is derivable, see Hermansen et al. (2023).The resulting loads are applied at the top and bottom of the shear webs that support the spar cap.However, since only the root section is being considered, the outer part has to be removed and its applied loads transferred to the remaining structure.In this work, the loads are applied directly as section forces on the cross section for simplicity.Consequently, local effects and deformation occur in the vicinity of the load introduction.To remove these effects, a reinforcement is added around the load introduction site and the patch adjacent to it is not parametrized for the optimization.To move a significant distance from the load site, the patch made larger than the other, see Fig. 2. The particular size of the de-parametrized patch is determined based on observations of the strain distribution.

Discrete Material and Direct Thickness Optimization
Optimal selection of materials in laminates is inherently a discrete programming problem.The discrete behavior of a given material candidate design variable x (p,l,c) in such a problem is as described in Eq. (1).
Here l is the given layer number and p is the patch number.For laminated composites, consisting of several material layers, the thickness design variables behave similarly, as they parametrize the presence of constant-thickness layers in each laminate composing the structure.The thickness design variables ρ (d,l) is hence described in Eq. (2).
Assignment of these variables to a structure is exemplified in Fig. 5.
Fig. 5: An example of a DMTO parametrization applied to optimize a cantilever plate.
The computational expense is excessive for solving discrete optimization problems of any modern industrial structure encountered in e.g., the wind turbine and aerospace sectors.Instead, if the optimization problem can be formulated with continuous design variables and functions, efficient gradient-based methods are available to solve such problems with available computational resources.The efficiency of gradient-based methods have been utilized to solve impressively large-scale problems see e.g.Aage et al. (2017); Baandrup et al. (2020), and are therefore also imperative in solving the present problem.This work utilizes the Discrete Material Optimization (DMO) approach (Stegmann and Lund, 2005;Lund and Stegmann, 2005) and its extension, the Discrete Material and Thickness Optimization (DMTO) approach (Sørensen et al., 2014) for enabling gradient-based methods.In DMO/-DMTO, the discrete design variables are replaced with continuous equivalents and discreteness is achieved by penalizing artificial intermediate design variable values.This penalization is achieved by utilizing an appropriate interpolation function that makes intermediate values expensive.In this work, the multi-phase version of the Rational Approximation of Material Parameters (RAMP) (Stolpe and Svanberg, 2001a;Hvejsel and Lund, 2011) is applied for interpolating both material and thickness design variables, see Eqs. ( 3 Here, q is the penalization factor, and penalization is achieved by using q > 0 in the weighting functions. As opposed to the popular alternative, Solid Isotropic Material with Penalization (SIMP) (Bendsøe, 1989;Zhou and Rozvany, 1991), RAMP allows the design variables to go to zero, as its gradient does not become singular.
The material candidate variables are incorporated into the structural problem through the constitutive matrix of each candidate material.In order to avoid solving an excessive amount of combinatorial problems, a weighted average C(e,l) of the candidate material constitutive matrices C (e,l,c) is utilized, see Eq. ( 5). (e,l,c)  (5) Here, e is the element number and N c denotes the number of candidates.Different strategies have been proposed for parametrizing the layer thickness.As this work concerns a wind turbine blade, which is majorly a sandwich structure, the alternative parametrization introduced in Sjølund et al. (2018Sjølund et al. ( , 2019)), henceforth named as the Discrete Material and Direct Thickness Optimization (DMDTO) approach, is utilized, see also L öffelmann (2021).In DMDTO the layer thickness t (e,l) is directly scaled with a density variable instead of the constitutive matrix and this parametrization strategy allows for defining different ply groups within each laminate in the model.Such ply groups can be defined for the face and core layers respectively allowing for variable thickness in each region that has a distinct ply group.Varying the thickness between the face and core ply groups is necessary to design sandwich structures, which utilize increased thickness of the core to increase the bending stiffness and hence reduce the amount of costly face materials.Such is not possible with the original DMTO, as a density variable is associated to the constitutive matrices and changing a layer thickness therefore introduces artificial weighting into the problem.The difference between the DMTO and DMDTO approaches is illustrated in Fig. 6.The thickness of each layer is thus parametrized with a continuous density variable ρ (d,l) through an interpolation function v(ρ), see Eq. (6).
Note here that element e is part of domain d.The parametrization is implemented with a dummy layer to offset the layup to the outer shell, see details in  Sjølund et al. (2019).In this context, the thickness of the outer layer is kept constant during optimization to ensure plies are dropped internally.
A number of additional constraints are necessary to make a DMTO approach yield realizable designs.The proposed multi-phase interpolation approach does not guarantee that only a single candidate material is chosen, and thus a set of linear constraints is included for each layer, requiring that the sum of candidate design variables does not exceed unity, see Eq. ( 7).
In Sørensen and Lund (2013), constraints on the change in thickness within the stack have been formulated to ensure that intermediate voids do not occur and these are given in Eq. ( 8).
It has been recognized that this straightforward formulation is insufficient as it forms through-the-thickness density bands, i.e., the optimizer places the same intermediate density value in all layers which is unaffected by penalization.The remedy introduced is a set of special move limits on the density variables, that ensures presence of material in underlying layers, see Eq. ( 9).
Here, T is a threshold value which enforces more density in underlying layers with decreasing T, and, hence, material is always removed from the top when applying the move limits, see Sørensen and Lund (2013) for further details.Using DMDTO, the thickness constraints are definable for each distinct sequence of the laminate, i.e. ply groups, allowing for using a different layer thickness in each ply group, see details in Sjølund et al. (2019).In the present work, such constraints are defined for the top, core and bottom ply groups respectively, see Fig. 7.
Fig. 7: The core-region laminate illustrated with application of the thickness constraints to each defined ply group.
Ply-drop constraints are also included in the optimization problem in the form presented in Sørensen and Lund (2013).Ply drops are known to cause stress concentrations as a consequence of an abrupt thickness change and build-up of resin pockets in the vicinity of the drop.The ply-drop constraints act to prevent this, and implementation can conveniently be made with offset in the patch parametrization, see Eq. (10).

−S ≤
Here, S indicates the maximum allowable number of plies to drop and N l is the number of layers.
Finally, the DMTO problem can be formulated, which is stated in the nested analysis and design formulation in Eq. ( 11).The structural constraints g(x, ρ) are defined in the subsequent section, where the associated structural analyses are presented.minimize l) , for top ply group ρ (d,l+1) ≤ ρ (d,l) , for core ply group l) , for bottom ply group

Structural Analysis
In this work, buckling, static failure and fatigue damage are all considered as part of the optimization.These structural criteria are all quantified through linear elastic finite element analysis.As such, the point of departure is solving the structural state equations, given in Eq. ( 12).

KU = F (12)
Here, K designates the global stiffness matrix, U is the displacement vector, and F is the work-equivalent nodal force vector.
From the resulting displacements, the strain vector in structural coordinates ε xy is determined on an element basis as shown in Eq. ( 13).
Here, B (e,l,m) is the strain-displacement matrix, for element e, layer l, evaluation position in the layer m, and u (e) is the element displacement vector.From the strain vector, the element stress vector in structural coordinates σ xy is calculated following Eq.( 14).Note that the constitutive matrix C(e,l) used here is the interpolated matrix calculated in Eq. ( 5).The strain and stress vectors are calculated at the top and bottom of each layer in every element, and hence m = 1, 2.

Buckling Analysis
Linear buckling analysis involves solving the eigenvalue problem shown in Eq. ( 15).
Here, λ BLF is the j th eigenvalue or buckling load factor, Φ (j) is the associated buckling mode, and N j is the number of buckling modes extracted when solving the eigenvalue problem.The stress-stiffness matrix K σ is computed on an element-level as shown in Eq. ( 16).
Mat G (e,l) dV (16) Here, G (e,l) is a matrix containing shape functions and these are paired with the corresponding stress components in the matrix S (e,l) Mat .The integral is evaluated numerically using Gauss quadrature.
The linear buckling problem is solved using a subspace iteration method with a static shifting factor of 1.5.For more details, see Bathe and Wilson (1972).Additional details on buckling in a DMO context can be found in Lund (2009).

Failure Analysis
For laminated composites, static failure is assessed for each layer in the laminate using an appropriate criterion.This work makes use of the maximum strain and stress criteria (G ürdal et al., 1999) for failure assessment, see also Lund (2018) for details on DMO-based failure optimization.These criteria are quite simple with no stress-interaction effects.Accurately assessing stress interaction is typically done using the Puck criterion (Puck and Sch ürmann, 2002;Deuschle and Puck, 2013), but including such in optimization is challenging due to the associated computational cost.Such criteria should hence be applied in later design phases.Moreover, the validity of simpler multi-axial failure criteria, such as Tsai-Wu and Tsai-Hill for general purpose failure evaluation has been questioned in recent research, see Talreja (2023).Max strain and stress criteria are therefore deemed suitable for preliminary design in general but is also accepted by design guidelines for wind turbine blades (DNV-GL, 2015).
For evaluation, the strain and stress vectors have to be rotated to the principal material directions in order to compare them to the material properties.The vectors are rotated by application of appropriate transformation matrices T, see Eqs. ( 17 denote strains and stresses in material coordinates (1, 2) respectively.Note, additional steps are necessary to derive strain transformation, see Jones (1999).The failure index FI (e,l,m) is then computed as the ratio of strain/stress in material coordinates to the strength of the material US, see Eq. ( 19).
FI (e,l,m) = max S (e,l,m) ⊘ US Here S (e,l,m) is represents either strain and stress, US is a vector containing the material ultimate strengths, and ⊘ denotes element-wise division (Hadamard division).

Fatigue Analysis
To ensure longevity of the blade, a high-cycle fatigue assessment is also included in the optimization.Offset is taken in an arbitrary random spectrum, which is quantified by using rainflow counting (Matsuishi and Endo, 1968).Rainflow counting yields a set of scaling factors c, that are used to scale the stress state according to the magnitude of the load, see Eq. (20).
Here, i indicates a discrete point in the quantified spectrum, and σ indicates the reference static stress state.Note that it is assumed that the loads act proportionally during the load history and as a result, the direction of the principal stress will remain constant.Proportional loading has profound implications on the computational efficiency, when solving the optimization problem, see  To assess fatigue damage, the alternating stress is compared to the material S-N curve.However, S-N curves derived from uni-axial tests and hence a single damaging equivalent stress σ (e,l,m,z,i) eqv is necessitated, which is to be calculated using an appropriate mean-stress correction method.In this work, the modified Goodman, and Constant-Life Diagram (CLD) interpolation methods are adopted for this purpose.The modified Goodman approach is used for evaluation in the core material candidates.Its expression is given in Eq. ( 23).
Here, z indicates the stress component.It should be emphasized, that the Goodman correction may yield non-conservative estimates of the influence of the mean stress for fiber-reinforced materials (Vassilopoulos and Keller, 2011).This motivates the choice of using CLDinterpolation to improve the accuracy for these materials.Inclusion in optimization has been demonstrated in Hermansen and Lund (2023) and reference is made to said work for derivation of the S-N curves for differing R-ratios, i.e., the ratio of valley to peak for a given cycle.
After determining the equivalent stress through either Goodman or properties derived through CLD interpolation, the number of cycles to failure N (e,l,m,z,i) can be computed using an appropriate S-N curve representation.Here the Basquin expression is adopted, which depicts a linear curve in a log-log plot, see Eq. (24).
Here, σ (z) f indicates the fatigue strength and b (z) is the slope on the curve for each stress component z.
Finally, damage D (e,l,m,z) is evaluated using the wellknown Palmgren-Miner expression to accumulate damage from every load combination in the spectrum N loads , see Eq. ( 25).
Here, n (i) is the number of loads with a particular combination of amplitude and mean, D lim is an adjustable critical fraction limit to introduce safety, and c L is a scaling factor that can be used to scale the damage in case of having a block spectrum where the same sequence of load cycles are repeated.The particular load spectrum used is generated using the Fortran intrinsic random number() function with N loads = 10 4 and c L = 10 3 .The limit D lim is implemented indirectly through the safety factors introduced in the following section and, hence, D lim = 1 is used.

Design Safety Factors
Wind turbine blade design standards given in DNV-GL (2015) dictates a safety factor that must be used when carrying out analysis.The magnitude of said safety factor depends on a number of partial factors assessing the criticality of various aspects related to each criterion.Calculation of the total safety factor γ m is given as the product of these partial factors, see Eq. ( 26).
Here, γ m0 is a base safety factor used for all analyses, γ mc is the factor accounting for criticality of failure, γ m1 is for long-term degradation of the materials used, γ m2 accounts for temperature effects, γ m3 is a factor for material and manufacturing effects, γ m4 accounts for the accuracy of analysis approach, and γ m5 is determined based on the number of load cases considered.The resulting safety factors for the present problem are given in Table 3. 5 Optimization Approach

Treating Singularities in Stress-based Optimization
The design space in stress-based optimization contains singular optima that are solutions present in a degenerate subspace and hence unreachable by gradientbased optimizers (Sved and Ginos, 1968;Kirsch, 1990;Rozvany and Birker, 1994).A couple of solutions to this problem has been proposed in literature.Here in the present work, the qp-method by Bruggi ( 2008) is adopted, since it has an additional benefit of further penalizing intermediate density.The approach entails using q < 0 in the RAMP interpolation functions of Eqs. ( 3) and (4).A negative penalty factor increases the stress artificially for design variables with intermediate density values, making them uneconomical for the optimizer to choose.The qp-method is applied the stress, failure index, and fatigue damage functions with a unique weighting function for each measure to allow for different degrees of penalization, see Eqs. ( 27

FI
(e,l,m,z,c) = w FI (x (p,l,c) )FI (e,l,m,z) (28) In this work, penalization is applied to the failure index and damage only.If adding penalization to the stress, it is more likely that the stress may exceed the ultimate stress of one of the candidate materials, making it impossible to perform mean-stress corrections and, hence, causing the optimization to end prematurely.An illustration of weighting factors with different penalization exponents is given in Hermansen and Lund (2023).The optimization functions are formulated from the weighted average of the candidate interpolated measures.For failure indices, the weighted average FI is calculated as shown in Eq. (30).
FI (e,l,m,z) = The interpolation is similarly applied to the damage, however, the measure has to be stabilized to combat its inherent nonlinearity.

Damage Scaling
Due to incorporating an S-N curve in the evaluation of fatigue damage, see Eq. ( 24), the measure becomes substantially nonlinear.Thus, using a constraint formulated with the damage as given in Eq. ( 25) makes the problem highly challenging to solve.Consequently, a stabilization of the damage measure is necessary to achieve acceptable results from the optimization.Numerous techniques for the purpose have been proposed in literature.This work will make use of an inverse P-mean norm scaling of the damage (Olesen et al., 2021) that combines the stabilization achieved, when scaling using a constant exponent (Lee et al., 2015;Zhang et al., 2019) and the accuracy of the true damage between zero and one, i.e., when the constraint is feasible.These features are achieved through a reformulation of the P-mean function and the introduction of a weighting factor c D , see Eq. ( 31).
Here, D s is the scaled damage, s is a scaling factor, and P D controls the accuracy of the measure.The factor c D governs the degree of stabilization and may also be changed using a continuation approach, where it is gradually decreased, hence approaching the original damage formulation, as convergence is approached.Such a continuation scheme is also adopted here and the details on it are presented later.After the damages have been scaled, a weighted average Ds is computed similarly to the constitutive matrix of Eq. ( 5) and failure index of Eq. ( 30), see Eq. (32). D(e,l,m,z) Note that different strategies can be adopted when scaling and interpolating the damage, however the approach presented here has performed the best in various numerical tests, see Hermansen and Lund (2023) for further discussions.

Local Aggregation
The failure indices and damages are locally defined measures, and the optimization problem therefore contains a large amount of both design variables and constraints making it both difficult and computationally expensive to solve.In this work, the number of constraints is reduced by P-norm aggregations, which compute a smooth, conservative, approximation of the maximum value g PN of each constraint, see Eq. (33).Here, g (y) is a constraint value, N y is the number of function values and P governs the accuracy of the approximation.A local aggregation approach (París et al., 2010;Le et al., 2010) is adopted where the Pnorm function is computed for each patch p in the problem such that less values are included in each P-norm in order to increase its accuracy in estimating the maximum value.Nevertheless, to achieve reasonable accuracy in the approximation, a relatively large value of P is necessary.As a result, the optimization problem becomes more nonlinear which is undesirable, as it potentially makes it difficult to solve.To make the aggregate approximation more accurate without using large exponents, the adaptive constraint scaling approach (Le et al., 2010) is used to re-normalize the constraints with information from previous iterations.The implementation is described in Oest and Lund (2017).It should also be noted that alternative techniques have been proposed in literature to increase the accuracy in aggregation itself, such as improved aggregation by Kennedy and Hicken (2015) or by using a maximum rectifier function in Norato et al. (2022).However, the techniques applied in this work are widely used, and are known to generally work well in the context of optimization of laminated composites (Lund, 2018;Hermansen and Lund, 2023).

Continuation Strategy
DMTO problems are difficult to solve in general because of the increase in non-convexity caused by interpolation of material candidates with penalization.Further exacerbating the non-convexity is the non-linear criteria used in the present optimization and it is therefore highly unlikely to find a satisfactory solution without a continuation approach.Continuation (Bendsøe and Sigmund, 2003) involves gradually increasing penalization during optimization which allows finding a solution to a less non-convex problem which can subsequently act as a stronger starting design for the next problem with increased penalization.However, modifying the penalization is a heuristic approach (Stolpe and Svanberg, 2001b), and the continuation scheme is therefore selected based on numerical experimentation.
In this work, continuation is applied to both the inverse P-mean function scaling factor and for design variable penalization.The strategies used for each factor are shown in Table 4.The penalty factors are modified after 30 iterations and c D continuation is enabled after 200 iterations where the factor is decreased every 15 iterations.Continuation is applied only after all constraints are feasible, as otherwise the optimizer may spend all iterations in a step on feasibility restoration instead of minimizing the objective function.It is therefore not guaranteed that the desired effect of approaching a strong minimum during each step is achieved.If the problem converges prior to reaching the set number of iterations between increments, the factors are modified immediately, and optimization is continued until all continuation steps are applied.

Design Sensitivity Analysis
The methods introduced in previous sections facilitates efficient computation of the gradients of every structural criterion through adjoint design sensitivity analysis.Common for the applied structural criteria is that differentiation with respect to a given design variable x (j) yields Eq. ( 34) as a result of applying the chain rule.
Here, f denotes a structural criterion function used in the optimization.The partial derivative of displacements with respect to the design variables dU dx (j) is determined by differentiating the finite element equations of Eq. ( 12).The resulting expression is shown in Eq. ( 35).
It is assumed that the load is independent of the design variables and hence ∂F ∂x (j) = 0. Isolating the desired derivative term dU dx (j) is achieved through inversion of the stiffness matrix, see Eq. ( 36).
This expression can at this point be solved, however it requires solving a system of equations for every design variable in the problem.In order to reduce the associated computational expense, an adjoint vector λ is defined as shown in Eq. ( 38).
The stiffness matrix is re-inverted to the left-hand side, and the terms are transposed.As the stiffness matrix is symmetric, it is unaltered by transposition, yielding Eq. ( 39).
The adjoint system of equations is solved once per constraint, which is convenient for the present problem formulation, where aggregation brings down the number of constraints.Note that the factored stiffness matrix from analysis is also directly reusable in this equation.The adjoint vector is then inserted into Eq.( 37), yielding the general adjoint sensitivity expression, see Eq. ( 40).
In this work, a semi-analytical approach is used where the adjoint vector is computed analytically, the stiffness matrix gradient ∂K ∂x (j) is found using central finite difference approximations and the partial derivative ∂ f ∂x (j) is found using forward finite difference approximations.For the finite difference calculations, a step size of ∆x = 10 −5 is used for all criteria.
The specific terms for each criterion can then be found by application of the adjoint method.Design sensitivity analysis of the buckling eigenvalues is available in literature and can be found in Sørensen et al. (2014).
Following the outlined adjoint approach, the final derivative expression is also achieved for this criterion, see Eq. ( 42).
Note that the partial derivative ∂F I ∂x (j) for candidate material variable is available analytically and is thus computed as such.
The derivative of fatigue damage is expanded in Eq. ( 43).
This is similarly reformulated using the adjoint approach to Eq. ( 44).
In this case, the partial derivative of the stress in structural coordinates with respect to design variables ∂σ xy ∂x (j)   is computed using forward finite difference for the thickness variables and analytically for the candidate material variables.For the remaining analytical partial derivatives in the failure and fatigue gradients, see Hermansen and Lund (2023).

Optimization Strategy and Auxiliary Methods
An efficient approach to formulate the optimization problem has now been proposed.Nevertheless, the DMTO problem is still difficult to solve, as it involves many design variables, structural, and non-structural constraints, see Eq. ( 11).The non-structural constraints are however all linear which motivates using a framework that is able to handle such well.This is offered in Sequential Linear Programming (SLP) which is thus adopted in this work.The optimization framework is built around the Sparse Nonlinear Optimizer (SNOPT) (Gill et al., 2005) that is used to solve the optimization problem.Auxiliary methods are implemented to successfully adapt SLP to the highly non-linear problem considered.The most basic addition is an adaptive move limit scheme.Such is essential to ensure that the optimizer does not oscillate between the two same points around a minimum due to the linear approximation of the optimization functions preventing the optimizer from converging.It also serves to improve convergence in general, as oscillations, that inevitably occur during the optimization, are treated promptly.
In the SLP framework, the structural functions are linearly approximated, which can lead to a poor representation during some iteration in the optimization.This may lead the optimization problem into infeasibility if not controlled properly, which will cause the optimizer to halt.In order to impose unconditional feasibility to the problem, a merit function approach (Svanberg, 2007) is adopted, see Sørensen et al. (2014) for details and settings.
Still, because of the inaccurate representation, it can be difficult to return to feasibility, because the linear approximations of the structural functions appear so to the optimizer.To control the infeasibility and also improve convergence, an SLP-filter (Chin and Fletcher, 2003;Sørensen and Lund, 2015;Hermansen and Lund, 2023) is introduced to reject poor solutions that may be found during the optimization.A poor solution is defined as having a larger objective function, in the context of minimization, or the infinity norm of the constraints being more infeasible than a previous solution.Moreover, solutions where any constraint ex-ceeds a set infeasibility limit of u lim = 0.01 are rejected and the optimization problem is solved again with reduced move limits.A new problem with tighter move limits is solved until a new solution that satisfies the criteria of the SLP filter is found or the convergence criteria are met.
To assess the quality of the achieved solution, the measure of non-discreteness (Sigmund, 2007) is utilized.This measure determines the amount of (non) discreteness in the design.Two measures are used, tracking the discreteness of the density variables M dnd and candidate variables M cnd individually (Sørensen et al., 2014).These expressions are given in Eqs. ( 45) and ( 46).
× 100% ( 45) Here, N d is the number of domains and V (e,l) is the volume of a layer.
Finally, a summary of settings for the various methods introduced in the paper is given in Table 5.

Results
This section will present the results attained through the approach described throughout the paper.As dictated by the safety factors determined in Subsection Here, m indicates the mass merit function, N p is the number of patches, see Fig. 2, and Manufacturing Constraints abbreviates the candidate sum, through-thethickness change, as well as ply-drop constraints of Eqs. ( 7)-(10).Note a constraint limit of 0.99 is used for the failure index FI PN and fatigue damage D PN constraints, such that the criteria never become violated, despite the allowance of u lim = 0.01 in the SLPfilter.Remark also, that these aggregated constraints are multiplied the adaptive constraint scaling factor described in Subsection 5.3.Convergence of the problems is defined as a relative change in design variables between iterations of 0.1%.In total, the presented problem contains 4848 design variables, 1744 structural constraints, and 6834 manufacturing constraints.

Problem Behavior
The problem converges in 204 iterations and the associated convergence history is shown in Fig. 8.It is observed that all criteria appear to be quite stable throughout optimization despite the applied continuation scheme.Particularly remarkable is the highly non-linear fatigue damage being controlled well by the applied inverse P-mean norm scaling strategy of Eq. ( 31).The mass has been reduced 23% from the starting point, see Fig. 8a, which is attributed to thickness removal, as the ply groups with fiber and core layers are predetermined as a consequence of the parametrization, see Fig. 7. Measures of non-discreteness of M cnd = 1.49% and M dnd = 4.06% are achieved for the resulting design, see Eqs. ( 45) and ( 46).Notable from all criteria convergence plots, Figs.8b-8d, is that the flap-wise load cases are critical, with associated constraints being active, or very close to, at convergence.The non-criticality of the edge-wise load cases may be a result of the leading and trailing edges are assigned the same initial layup as the spar cap.Due to the layup of these regions being overdimensioned in the initial design, the optimizer could have settled on a poor local minimum prematurely.An improvement could be achievable by modifying the initial design to better reflect the loading applied to these regions.
It is evident from the convergence plots, that particularly the edge-wise fatigue damage constraints are significantly reduced, when the inverse P-mean continuation scheme, see Table 4, is applied, see Fig. 8d.Unfortunately, the optimizer seems to have settled well on the optimum once the inverse P-mean continuation starts and thus the final design is not impacted much by the change.If improvements are desirable, the continuation scheme could be modified accordingly through numerical experimentation.

Optimized Laminate Design and Structural Response
The optimized laminates for the fiber regions are illustrated in Fig. 9, see also Fig. 1 for the denominations used.As observed, the spar caps have distinct material choice, with the majority being UD and the top layer being biax, see Figs. 9b and 9c.The top layer is subject to the largest bending moments and therefore needs the biax layer as UD does not exhibit sufficient transverse and shear strength.However, the failure and fatigue constraints are not quite active in the spar caps.The reason is likely that the thickness of the top layer is fixed during optimization indicating the outer biax layer could be reduced further.The failure constraint is however active in the following UD layer on the downwind side, see the failure distribution in Fig. 11c, due to significant shear and a biax layer could therefore be desirable here.A probable reason for not achieving an active constraint here is the optimizer settling on material choice faster than thickness distribution as a result of limitations imposed by the thickness constraints of Eq. ( 8).
The leading edge plies have been selected as mostly biax, see Figs. 9d and 9e.This choice seems to be motivated by a combination of the load case not being severe, an initial design that is far from the optimum, and the thickness constraints that slow the removal of material.It is observable that the leading     9a.Zero thickness indicates the outer mold, as the layup is offset from here in the shell element formulation.The limit on the y-axis is the thickness of the original layup, i.e. 80 mm.edge is only close to being critical locally in static failure, see Fig. 11c, which further supports, that more material could be removed in this region.Furthermore, the laminates are observed to contain some nondiscrete plies, which can also be a result of the loading being less severe in this area, thus requiring additional penalization to reach a discrete choice of material (or reducing the laminate thickness in the initial design as previously discussed).
In the trailing edge laminates, see Figs. 9f and 9g, the choice of material is distinct.The laminates however contain cross plies, particularly in the downwind side, which is attributed to a local minimum that is reached for the same reasons as the leading edge.However, the presence of cross plies is not as pronounced as the presence of biax in the leading edge.The level of static failure is only moderate in the trailing edge, however buckling is critical near the root, see Fig. 11b.
A general observation is that local thickness valleys, i.e., significantly lower thickness in a patch compared to those adjacent to it, are not particularly prominent in the optimized fiber regions.This is likely also a result of using the thickness constraints, which governs how much the thickness of a stack can change at a time.Thus, using the thickness constraints may help reduce the amount of post-modification of the design, if the design has to conform to e.g., ply-drop rules.On the other hand, the constraints may cause the design to end up in a poor local minimum, as observed with the edge regions, and if design exploration is the objective of the optimization, it may be undesirable to include these.
The optimized core laminates are illustrated in Fig. 10.As observed, the material choice is close to being distinct for every laminate, with the exception of a mostly cross-ply layer in the trailing edge core regions, see Figs. 10d and 10f.The laminates are observed to contain more local thickness valleys compared to the fiber-dominated regions, which is a consequence of the thicker core layers.The main design driver for the trailing edge core regions, shown in Figs.10d-10g, is buckling, which is pronounced in both the downwind side for load case 1, illustrated in Fig. 11a, and the upwind side for load case 3. Static failure is also observed to be critical near the root in these regions as seen from Fig. 11c.For the leading edge on the upwind side, static failure and fatigue damage seem to be the design driver as evident from Figs. 11c and 11d.There is however a portion of the blade that is not subject to significant static failure, and thus more material could likely be removed in the region.On the leading edge core downwind side there are a few local areas with static failure close to being critical, however it is evident that more material can be removed from this area.The non-criticality of these regions may also be attributed to the ply-drop constraints, which dictate that only one ply drop is allowed between adjacent plies, and thus the thickness of the spar cap may be the dimensioning factor.

Discussion
As evident from the paper, several methods and simplifications are necessary to be able to solve the presented problem.Although the problem treated is challenging, it behaved quite stable, as evident from the convergence plots, as a result of the applied approach.Nevertheless, a number of alternative approaches are available which could improve the results.
The analyses used in this work has been kept as simple as possible which can be acceptable for preliminary design.However, it could be interesting to include more detail through more sophisticated analyses to see if such leads to further improvements in the achieved design -either as additional reductions in objective or less tedious postprocessing.For buckling evaluation, simple linear analysis is applied in this work, however typical blade design is carried out using non-linear analysis.In fact, application of nonlinear buckling analysis allows for reducing the safety factors dictated by design guidelines, see Subsection 4.4.The inclusion should however be carefully considered with respect to the increase in computational cost associated with using a non-linear solver.
In terms of static failure, multi-axial effects are not accounted for by the simple max strain and stress criteria.Such could easily be included through e.g., the well-known Tsai-Wu or Hashin criteria.Additionally, the design guidelines also prescribe inter-fiber failure evaluation using Puck's criterion.The inclusion of Puck in the optimization is challenging, as it is essentially a critical-plane approach, which involves an exhaustive search for the plane where the maximum strain/stress are acting.Similar could be done in the fatigue analysis to further improve the fatigue damage estimation, see e.g.Dong et al. (2016).The computational cost of such a search is substantial and its inclusion is therefore not straightforward for the problem considered.Optimization with critical-plane methods has been demonstrated in Sartorti et al. (2023) of both academic and industrial metal structures.The authors apply the critical-plane approach to hotspots identified prior to the optimization, e.g. the rounding in the well-known L-beam example, and the fatigue analysis therefore typically involves 10 3 elements.As it can be difficult to know a priori, where the hotspots are on a  wind turbine blade, it is however likely not feasible to adopt such an approach.Note also that, the fatigue model could also be extended to take into account additional fatigue-specific effects, which is discussed in detail in Hermansen and Lund (2023);Hermansen et al. (2023).
As opposed to previous works (Sjølund and Lund, 2018;Hermansen et al., 2023), this work has used only four design load cases in the analysis.This is acceptable using up-scaled safety factors according to the design guidelines, however it may be desirable to include twelve load cases instead as it may change the design.Noticeable in Hermansen et al. (2023), where the same blade root section is considered, is that a mixed flap-and edge-wise load case is a critical load case and the associated constraints therefore quickly become active during optimization.This particular load case is not included in the present, where only extremal flap-and edge-wise load cases are included, which is acceptable according to design standards, but it is questionable if the design is actually conservative by increasing the safety factors with the prescribed amount for the given case.On the other hand, using four load cases may be sufficient for the preliminary design, which could be followed by e.g., a thickness sizing optimization, with additional load cases and higher-fidelity analyses.
A counterweight to the increase in computational cost of including more load cases, using non-linear criteria, etc., can be to use linear shell elements and a fully analytical design sensitivity analysis framework.The currently most costly part of the framework is computation of the finite difference approximations to the partial derivatives presented in Subsection 5.5.Deriving the analytical gradients could therefore allow for increasing fidelity elsewhere e.g., of the model, analyses, or problem.This is left for future work.
The inverse P-mean damage scaling used for stafatigue damage is observed to handle the non-linearity well, but the damage is not well-distributed in the optimized design unlike what has achieved in previous works e.g, Olesen et al. (2021); Hermansen and Lund (2023);Hermansen et al. (2023).The resulting damage distribution depends on finding the correct continuation strategy to the inverse P-mean scaling factors which can be difficult to determine.Through numerical testing, a better scheme could be found, but it is undesirable for this problem because of the associated computational cost of solving it.
The ply-drop constraints used in the work are simplistic and are formulated as allowable number of layers in adjacent patches, see Eq. ( 10).It is desirable to formulate more sophisticated ply-drops constraints that are independent of the patch parametrization.In Sjølund and Lund (2018); Hermansen and Lund (2023) these have been taken into account in an average sense over the design patches, but this implies postprocessing is necessary to realize the design.In the alternative multi-level optimization framework, the final step is a realization step, where numerous manufacturing constraints are taken into account, see e.g.review by Albazzan et al. (2019).However, this is typically done using gradient-free approaches, which limits the size of the problems treated.

Conclusion
This paper has presented a framework for carrying out efficient multi-material and topology optimization of wind turbine blades.Creation of the blade model root section for the purpose of optimization is shown along with application of four design load cases, which criticality are assessed in subsequent structural analysis.Parametrization of the structure is carried out using the Discrete Material and Direct Thickness Optimization (DMDTO), which is particularly suitable for sandwich structures that constitutes a majority of the blade structure.The structure is subject to linear buckling, static failure, and fatigue damage constraints for the four design load cases, and the associated analysis procedures are hence defined.The two latter structural criteria are particularly difficult to handle in optimization and a number of methods are therefore introduced to improve convergence of the problem, specifically the handling of singular optima, damage scaling, and constraint locality.The design sensitivity analysis is then presented along with a semi-analytical approach to compute the resulting gradients.A Sequential Linear Programming (SLP) framework is subsequently described, which uses a number of auxiliary techniques to further enhance the likelihood of finding a strong optimum.These are a continuation approach, adaptive move limits, a merit function, and an SLP-filter.The effectiveness of the approach is demonstrated, yielding a reduction in mass of 23% with respect to the initial design.Variable-thickness laminates are achieved in the optimized design, and these are almost fully discrete.Particularly good distributions of buckling and static failure are achieved, with fatigue damage being more localized due to its highly nonlinear behavior.The work emphasizes the potential of applying DMDTO in the early design phases, which could offer the designer an excellent point of departure for subsequent detail design.

Fig. 2 :
Fig. 2: The division of the blade into thickness design variable patches.

Fig. 3 :
Fig. 3: The finite element mesh and parametrization illustrated from different perspectives.

Fig. 4 :
Fig. 4: The four load cases considered in the optimization.

Fig. 6 :
Fig. 6: The difference between the DMTO and DMDTO parametrizations.The left is the topology optimization-based DMTO, where the constitutive matrices are parametrized and the left is DMDTO.
Zhang et al. (2019);Sartorti et al. (2023); Hermansen and Lund (2023) for details.The scaling factors are subsequently divided into amplitude and mean components, yielding a corre- (21) and (22) respectively.σ (e,l,m,i) a g (a) Convergence of the mass function.(b) Convergence of the buckling functions.(c) Convergence of the failure index functions.(d) Convergence of the fatigue damage functions.

Fig. 8 :
Fig.8: Convergence of the criteria functions.Modification of the penalty factors is observable, when the constraints suddenly increase in value, most notable for the static failure and fatigue damage constraints.
(a) Material designations for the fiber-dominated laminates.(b) Initial thickness distribution of the spar cap on the downwind side.(c) Optimized thickness distribution of the spar cap on the upwind side.(d) Initial thickness distribution of the leading edge on the downwind side.(e) Optimized thickness distribution of the leading edge on the upwind side.(f) Initial thickness distribution of the trailing edge on the downwind side.(g) Optimized thickness distribution of the trailing edge on the upwind side.

Fig. 9 :
Fig.9: Laminate design for regions containing majorly UD layers.Material designations are given in Figure9a.Zero thickness indicates the outer mold, as the layup is offset from here in the shell element formulation.The limit on the y-axis is the thickness of the original layup, i.e. 80 mm.
(a) The first buckling mode of load case one.(b) The second buckling mode of load case two.(c) Envelope plot of the maximum failure index for all load cases in the optimized blade.(d) Envelope plot of the maximum damage component for all load cases in the optimized blade.

Fig. 11 :
Fig. 11: Optimized structural response of the blade root section.Note the de-parametrized patch next to the load introduction is not included in the visualization.

Table 1 :
Elastic and static strength properties used.

Table 2 :
Fatigue strength properties used.

Table 3 :
Partial factors and resulting safety factors applied for each criterion (DNV-GL, 2015).

Table 4 :
Values applied in each weighting function, and how they are changed during optimization.

Table 5 :
An overview of the various parameters used to solve the optimization problem.
4.4, the right-hand side of each respective constraint is given, thus the optimization problem can be established, see Eq. (47).