Land-based wind turbines with ﬂexible rail transportable blades – Part II: 3D FEM design optimization of the rotor blades

. Increasing growth in land-based wind turbine blades to enable higher machine capacities and capacity factors is creating challenges in design, manufacturing, logistics, and operation. Enabling further blade growth will require technology innovation. An emerging solution to overcome logistics constraints is to segment the blades spanwise and chordwise, which is effective, but the additional ﬁeld-assembled joints result in added mass and loads, as well as increased reliability concerns in operation. An alternative to this methodology is to design slender ﬂexible blades that can be shipped on rail lines by 5 ﬂexing during transport. However, the increased ﬂexibility is challenging to accommodate with a typical glass-ﬁber, upwind design. In a two-part paper series, several design options are evaluated to enable slender ﬂexible blades: downwind machines, optimized carbon ﬁber, and active aerodynamic controls. Part 1 presents the system-level optimization of the rotor variants as compared to conventional and segmented baselines, with a low-ﬁdelity representation of the blades. The present work, Part 2, supplements the system-level optimization in Part 1 with high-ﬁdelity blade structural optimization to ensure that the 10 designs are at feasible optima with respect to material strength and fatigue limits, as well as global stability and structural dynamics constraints. To accommodate the requirements of the design process, a new version of the Numerical Manufacturing And Design (NuMAD) code has been developed and released. The code now supports laminate-level blade optimization and an interface to the International Energy Agency Wind Task 37 blade ontology. Transporting long, ﬂexible blades via controlled ﬂapwise bending is found to be a viable approach for blades up to 100 m. The results conﬁrm that blade mass can 15 be substantially reduced by going either to a downwind design or to a highly coned and tilted upwind design. A discussion of active and inactive constraints consisting of material rupture, fatigue damage, buckling, deﬂection

manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

Introduction
Wind turbine rotors have been growing at a faster rate than generator capacity, enabling higher energy capture in low and moderate wind speeds (Bolinger et al., 2020). Thus, the ratio of the nameplate capacity to the swept area, or specific power (W m −2 ) is diminishing. Using low specific power turbines has increased wind plant capacity factors which is a measure of how 30 often maximum power is produced. This trend of growing capacity factors (i.e., lower specific power) allows modern turbines to both produce power more consistently and operate in low wind speed sites. Current blades are approaching 80 m in length for land-based installations and over 100 m for offshore turbines. Bolinger et al. (2020) has shown that the ongoing reductions in specific power are likely to continue. While not as much of an issue for offshore wind turbines, land-based machines are currently constrained by transportation logistics. Current land-based transportation constraints limit monolithic blade lengths 35 to about 80 m in length, and 4.75 m in width and height. Thus, continued growth in rotor sizes for land-based turbines will likely require design changes.
In 2018, the U.S. Department of Energy (DOE) funded the Big Adaptive Rotor (BAR) project to study the design drivers of future high capacity factor land-based turbines, and investigate potential technology solutions to the challenges that these designs create. The project used a 5 MW, 206 m rotor as a reference turbine platform to examine the limits of current design 40 tools and methodologies. After an initial study and down-select of turbine innovations (Johnson et al., 2019), the concept of slender flexible blades that can be transported by train via controlled flapwise bending (Carron and Bortolotti, 2020) was chosen as the primary concept to be evaluated. Three enabling technologies were selected as having high potential to enable this concept: downwind rotors, carbon fiber, and distributed aerodynamic controls. Downwind rotors, while having noise and cyclic loading concerns, experience their highest deflections away from the tower, lessening that constraint. Carbon fiber, while 45 expensive in its current form, can be used in blades to enable thinner airfoils to allow higher flexibility while maintaining the required structural integrity. Reference rotor models were designed and optimized for each of these technologies. Each design is listed below as well as the nomenclature adopted.
-BAR-UAG (Upwind -Air transport -Glass fiber spar caps): A baseline upwind design composed primarily of fiberglass composites with 4 m of blade prebend in a manner similar to current industry standard designs.

50
-BAR-DRG (Downwind -Rail transportable -Glass fiber spar caps): A straight, slender, and downwind design composed primarily of fiberglass composites, intended to be transportable by train on existing railways.
-BAR-DRC (Downwind -Rail transportable -Carbon fiber spar caps): A straight, slender, and downwind design using industry-standard carbon fiber composite for the spar caps, intended to be transportable by train on existing railways.
-BAR-USC (Upwind -Segmented -Carbon fiber spar caps): A segmented upwind design, composed of two sections 55 using industry-standard carbon fiber composite for the spar caps and attached rigidly together by a mechanical joint.
-BAR-USCHT: A variation of the BAR-USC design, using heavy-tow carbon fiber composite for spar caps instead of the industry standard.
-BAR-URCHT: A variation of the BAR-URC design, using heavy-tow carbon fiber composite for spar caps instead of the industry standard.
Optimal designs were obtained by performing numerous iterations between two levels of optimization: each with differing 70 structural fidelities. The National Renewable Energy Laboratory (NREL) led the lower-fidelity system optimization, wherein the levelized cost of energy was minimized while capturing the overall behavior of the system at the turbine level. This was performed first, whereby a preliminary blade design was created. Further details are presented in Part 1 of this two-part series (Bortolotti et al., 2021).
Beam models, in part, enabled the design of the aerodynamic shape and sizing of the spar cap while accounting for stochastic 75 wind fields and the interaction of numerous turbine components. These types of models are used at the expense of accurately quantifying ultimate and fatigue failure of the blades, as well as buckling instabilities. Some researchers have accounted for buckling instabilities with beams by using simple analytical buckling formulas (buckling is a form of instability) (Ning and Petch, 2016;Bir, 2001). These formulas, however, are often derived for flat, homogeneous, and isotropic panels with unrealistic boundary conditions (i.e., fixed, free, pinned, etc); neither of these conditions apply to most panels in a given blade. Inadequate 80 stress recovery for beam models has also been shown due to blade taper (Bertolini et al., 2019) effects and localized effects near the root are not able to be well-resolved. Note that linear cross-sectional analysis tools (Feil et al., 2020) were utilized in the lower-fidelity optimization. Nonlinear cross-sectional approaches exist (Harursampath and Hodges, 1999;Couturier and Krenk, 2016) and can account for flattening of the cross section under flexure, known as the Brazier effect. These, however, are limited to standard structural cross sections and will not detect localized buckling of a panel without simplifications.

85
In an effort to ensure the BAR designs were neither over or underdesigned, the system-level optimization was followed by a high-fidelity structural optimization. The increase in fidelity required keeping the aerodynamic shape unaltered and relying on the lower-fidelity model for the loading information. The present contribution details the high-fidelity methods and associated failure modes employed by Sandia National Laboratories (SNL). It differs from a related work (Bottasso et al., 2014) in that the structural design characteristics of long, slender wind turbine blades for large, land-based rotors are revealed. Solid finite elements (FEs) with a layer-wise discretization can provide the best approximation to the 3D elasticity boundary value problem. This is because no simplifications or assumptions on the kinematics, local fields, and constitutive response are required. Designing a blade with a model that exclusively uses solid elements would require far too many degrees of freedom due to the vast separation of scales between layer thicknesses and the overall blade size. Even if layers were grouped with a homogenization theory, the degree of freedom requirement would still be too high for practical design. Shell elements (ANSYS 95 SHELL181) were exclusively utilized in this work because it is a good compromise between fidelity and computational cost.
This design approach required numerous improvements to the MATLAB ® code named Numerical Manufacturing And Design (NuMAD) 2.0 (Berg and Resor, 2012). NuMAD 2.0 had a strong dependency on a graphical user interface (GUI), which interferes with optimization. This culminated in the current release of NuMAD; NuMAD 3.0. The new release incorporates optimization, added structural analyses, and capability to accept input from the International Energy Agency Wind Task 37 100 blade ontology (referred to here as a yaml file (Bortolotti et al., 2019)). It includes the GUI but with the addition of an objectoriented approach. This was done in order to automate the new optimization process, which utilizes MATLAB ® gradient-based optimization methods and various calls to ANSYS ® to perform meshing and structural analyses. The code can be obtained from https://github.com/sandialabs. The remainder of this paper summarizes the methods used to derive the material properties needed in the structural analysis 105 and the application of aeroelastic loading to a high-fidelity structural model, and provides details regarding the structural optimization of the wind turbine blade. Specifically, Sect. 2 details how a complete set of material properties required for the optimization and analysis were obtained. Section 3 -4 detail the updates to NuMAD. Section 3 shows how the loading information was transferred from the lower-fidelity aeroelastic simulations to the higher-fidelity case. Section 4 describes how the high-fidelity structural optimization was performed. Finally, results are presented in Sect. 5 and conclusions are made in 110 Sect. 6.

Material Properties
A common difficulty for the composites simulation community is that publicly available experimental data sets often provide a limited subset of inputs needed for accurate modeling. Experimental data for only a unidirectional glass composite was found to have been fully characterized experimentally with elastic, rupture, fatigue, and mass properties (Samborsky et al.). The 115 rest of the composites required micromechanical analyses and approximate rupture calculations for laminates. This section elaborates on the assumptions and calculations adopted to generate the material properties. The required properties for the composite laminates are shown in Table 1 and 2. Note that γ σ is the factor of safety for material rupture.

Elastic Properties
The only material found in a publicly available data set with fully characterized elastic properties was the uniaxial glass laminate. Other materials were largely limited to a few properties. The following three sections show how experimental values of the uniaxial glass laminate were distilled and summaries of how homogenization theory was used for the rest of the composites.

Experimentally Determined Properties
All nine of the elastic properties of the uniaxial glass laminate were obtained directly from experimental values (Samborsky et al.). Since the tensile and compressive effective elastic moduli, E * , were not equal, the arithmetic mean of the respective 130 tensile and compressive value was utilized. For the shear moduli, G * , the arithmetic mean of G * 12 and G * 21 was reported as G * 12 . Likewise, the arithmetic mean of G * 13 and G * 31 was reported as G * 13 and the arithmetic mean of G * 23 and G * 32 was reported as G * 23 . The approach taken for the shear moduli is not applicable for the Poisson ratios, since in general, ν * ij = ν * ji , where ν * ij are Poisson ratios. ν * 12 , ν * 13 , and ν * 23 were used directly since ν * 21 , ν * 31 , and ν * 32 were obtained by loading the laminate in unconventional ways.

Analytically Determined Properties
Homogenization, a formal subset of micromechanics, is routinely performed to enable efficient analysis and design. The biaxial and triaxial glass laminates both comprise uniaxial fabric placed at a number of angles that are then stitched together. Thus, they are both able to be respectively homogenized as a stack of homogeneous layers. For such a case, Yu (2012) has rigorously proved that the in-plane strains and transverse stresses are uniform throughout the thickness of the laminate. Thus, the effective 140 properties of the stack of layers can be obtained from a hybrid rule of mixtures in an exact manner. Thus, Yu's approach was adopted here for both of the laminates.
Note from Table 1 that one of the Poisson ratios of the biaxial glass is greater than 0.5. This is permissible since the commonly known maximum limit of 0.5 is exclusively for isotropic materials and this homogenized material is anisotropic. The repeating unit cell considered for the biaxial glass composite was two layers of the uniaxial glass laminate with a stacking sequence of

Numerically Determined Properties
The layers of the carbon fiber laminates are all in the same direction and have not been characterized as fully as the uniaxial glass laminate. Thus, it is necessary to consider fiber-level details for the homogenization. The microstructure considered here was a hexangular packing of circular fibers, as shown in the RUC of Fig. 1. The local stress and strain fields in this were obtained with a numerical model, known in micromechanics as representative volume element (RVE) analysis. The periodic boundary conditions were used because they are known to satisfy the well-known Hill-Mandel macrohomogeneity condition and only one unit cell is required in the RVE for converged properties. Further details of RVE analysis can be found in Yu (2019). The elastic modulus of the fiber along the fiber axis was calibrated from the corresponding tensile test data in Miller et al. (2019) and the Voigt rule of mixtures. We obtained 230.26 GPa for the industry standard carbon fiber and 234.69 160 GPa for the heavy-tow fiber. The rest of the fiber properties were obtained from AS4 material properties in Herráez et al. (2020). As for the isotropic epoxy, the elastic modulus of 3.2 GPa and Poisson ratio of 0.347 was used.

Rupture Properties
Table 1 also shows the strength values for each laminate. Note that X is the longitudinal tensile strength, Y is transverse (inplane) tensile strength, and their primed counterparts are the compressive strengths. S is the in-plane shear strength, and R and 165 T are the shear strengths associated with σ 23 and σ 13 , respectively. Note that all of these are characteristic layer strengths.
All of the strengths for the uniaxial glass laminate were directly obtained from physical test data. For the remaining glass laminates, test data were only available for the in-plane strengths. The experimental data for the remaining laminates came from the SNL/MSU/DOE Composite Materials for Wind database (SNL/MSU/DOE). The out-of-plane shear strength properties, R and T , for the biaxial and triaxial glass laminates, are set equal to the out-of-plane properties of the 0 • laminate because they 170 are assumed to be insensitive to fiber architecture and layup.
The SNL/MSU/DOE Composite Materials for Wind Database does not contain test data of the in-plane shear strength, S, for biaxial and triaxial glass laminates. These values were both determined with using a failure analysis technique that incorporated progressive damage analysis with classical lamination theory and LaRC03 (Davila and Camanho, 2003). Progressive damage analysis was implemented as ply stiffness degradation after the ply was determined to have failed. The ultimate strength was determined by the last ply failure method.
X and X for the uniaxial carbon laminates were obtained from the 95 % value (the same value used in Ennis et al. (2019)).
The rest of the strength values for both carbon fiber laminates were assumed to be the same as the uniaxial glass laminate. The fatigue exponents, m, for the glass are standard values from DNVGL-ST-0376 (2015), whereas for the carbon fiber laminates, they are from Ennis et al. (2019).

Design Loads
Computational resources today cannot accommodate for fluid-structure interaction models in conjunction with a dynamic shell model of a blade for each design load case (DLC). Therefore, it is necessary to construct a reduced set of static load cases that is representative of the most critical loads from the beam model's time history analysis. The designs identified by the low-fidelity optimization were run through the aeroservoelastic solver OpenFAST (NREL, 2021a). Load equivalency 185 between the OpenFAST beam models and the shell models constructed by NuMAD was established by equating the resultant forces and moments due to the stresses at various spanwise cross sections. Thus, the loading conditions were derived from the results of OpenFAST from NREL. The ultimate and fatigue limit states were evaluated from the International Electrotechnical Commission (IEC) standard (IEC-61400-1, 2005). Cases 1.1, 1.3, 1.4, 1.5, 5.1, 6.1, and 6.3 were considered for the ultimate limit states and 1.1 and 1.2 were considered for the fatigue limit states.

190
It is important to note that these cross-sectional resultants from OpenFAST are in coordinate systems that are in the deformed configuration and are aligned with the local principal axes (structural) of the cross section. Thus, those coordinate systems change with respect to span and time. Figure 2 contrasts the OpenFAST results coordinate basis vectors, v (j) i (j = 1, 2, 3, ..., k), with those of the NuMAD loads system, w i , and the ANSYS coordinate system, x i , which are invariant. Note that the x i coordinate system is like the w i system but the first axis, x 1 points toward the leading edge instead of the flap direction.

195
A small angle between v (j) 3 (j = 1, 2, 3, ..., k) and w 3 exists but is unaccounted for in NuMAD. Thus, it is assumed that v (j) 3 = w 3 = x 3 (j = 1, 2, 3, ..., k). The remaining angle between the v (j) α and w α is referred to as µ, which varies with span. Here and through the rest of this document, Greek indices assume 1 and 2 (except where explicitly indicated). These coordinates are illustrated in Fig. 3(a) and are related by where C ij are the direction cosines that can be stored in a matrix as shown below In this work, the principal axis angle was obtained by PreComp (Bir, 2006).  Overall, two load categories were deemed necessary to properly analyze the most critical loads; those needed to evaluate the maximum tip deflection and those needed for evaluating blade failure. Here, blade failure consists of ultimate failure, buckling, 205 and fatigue failures. Both cases utilized the OpenFAST resultant axial forces, F w 3 , both bending moment components, M w α , and the torsional moment, M w 3 . The superscript indicates the reference frame. These forces and moments will be referred to collectively with P i (i = 1, 2, .., 4), a generalized load vector, where varies with span-wise location and that resultant shear forces that are transverse to the blade were assumed to be negligible for establishing load equivalency for both load cases. For the tip-deflection case, the time at which the maximum tip deflection 210 from the beam model, t * , was determined. Then the P i , at a given cross section was defined by The components were then transformed to the x i coordinate system.

220
Load directions from 180 • ≤ θ < 360 • were, however, obtained by Note in Eq. (6) that the minimum of M y 1 is found instead of the maximum but θ still ranges from 0 • to 180 • . Unlike the resultants used for the deflection analysis, which all occurred at a specific time in the OpenFAST simulations, each of the axial force, torsion, and bending moment resultants along the span could possibly come from different times. Thus, it is an artificial 225 distribution suitable for design use. Also, unlike the deflection case, the two bending-moment components (i.e., flapwise and edgewise bending) at a span location were projected onto eight directions, as defined by y 1 in Fig. 3(b). The loads in these eight directions were also transformed to the x i system.
For any given distribution of bending moments on the blade, whether for the deflection analysis or evaluating failure, the forces to be applied to the blade were obtained by writing mechanical equilibrium expressions for each spanwise position of 230 interest. Fig. 4 shows the known spanwise bending moments, M i acting at a distance, z i , from the blade root and the forces to be solved, and F i acting atz i , wherez i = 1 2 (z i + z i+1 ). The forces to be applied are found by solving the following linear where k is the number of cross sections with resultant moment data. Because OpenFAST only outputs this moment data at a 235 maximum of 10 blade cross sections, the moments were with a shape-preserving piecewise cubic interpolation. In this way, force distributions parallel to w 1 and w 2 are obtained. Finally, these force distributions were transformed to nodal forces on the blade outer shell by following Berg et al. (2011) but with an extension to apply axial forces and twisting moments.

Main Optimization Loop
Once the preprocessing steps were completed for each design, the main structural optimization process was performed in 240 NuMAD. The following sections describe the definition of the optimization objective function, design variables, and algorithms/methods used.
The overall objective of the high-fidelity structural optimization in NuMAD was to minimize mass of the blade structure subject to the following four performance constraints: 1) resistance against material rupture, 2) resistance against fatigue damage over operating life span of the turbine, 3) resistance against structural buckling, 4) maximum tip deflection not to

Constraints
For each iteration in the optimization, NuMAD calls ANSYS to obtain the constraint values and stress data, as well as to compute blade mass for the objective function. Figure 5 illustrates the FE workflow for a single iteration of the optimization.
The total blade mass and frequencies are first computed. Then the load data that is associated to the maximum tip deflection from OpenFAST was applied to the blade for a single linear-static analysis to obtain the tip deflection. Then blade failure (comprising material rupture and buckling) was evaluated by considering each load direction, θ. For a given load angle, a 255 material rupture index was evaluated at every location and layer of the blade from a linear-static analysis. This was computed by utilizing a Tsai-Wu failure criteria that neglects any effects due to transverse normal stress (ANSYS, 2017). The maximum rupture index for the whole blade is stored, σ θ (i = 1, 2, 3, ...n θ ), as well as stress data for each layer of every element for the fatigue postprocessor described in Sect. 4.1.1. Note that n θ is the number analysis directions and that the default value of -1 in ANSYS was used for the three shear coupling coefficients in the Tsai-Wu criteria. The buckling load factor associated to 260 the lowest mode, β θ (i = 1, 2, 3, ...n θ ), was then stored from a subsequent eigenvalue buckling analysis of the entire blade.
This process repeats for every load direction. It is notable that geometric and material nonlinearity is unaccounted for. These constraints are summarized mathematically, as follows: Material Rupture: Fatigue Damage: Buckling: Maximum Tip Deflection: Minimum Flap Frequency: Minimum Edge Frequency: In the above, u tip is the highest tip deflection of the blade in the direction normal to the rotor's sweep plane, and u max is the highest allowable blade deflection to ensure tower clearance. f f lap and f edge are the first natural vibration frequencies in the flapwise and edgewise directions, respectively, and f rot is the frequency of the rotor at the rated speed. γ D , γ β , γ u , and γ f are the factors of safety applied to fatigue, buckling, tip deflection, and natural frequencies, respectively (γ σ defined in Sect. 2). Note that γ D is set to 1.0 because the fatigue factor of safety is implemented in Eq. (18). The rest of the design values are 275 listed in Table 3.

Fatigue Damage Postprocessor
NuMAD 3.0 now includes the capability to evaluate fatigue damage for every material layer at a discreet number of spanwise stations. First, stress data in the local layer coordinate system for every element and every layer is read into MATLAB ® , as  well as element location along x 3 . This is done for both the y 1 (θ) and y 1 (θ + 90 • ) directions to consider the additive fatigue  Fig. 4). Thus, the NREL Crunch (Buhl, 2003) program was setup to analyze each wind speed file with "calculated channels" for M y 1 at each spanwise location and direction. Because where σ a and σ m are predominant normal stress amplitudes and means, respectively, and M a and M m are the amplitude moments and mean moments, respectively. The scalar parameter α is obtained by where σ 11 is the normal stress in the fiber direction of a given layer, and M r (θ, x 3 ) is the magnitude of the resultant moment applied to the cross section. Thus, the moment-based Markov matrices are converted to stress-based Markov matrices suitable for damage calculations. The well-known Palmgren-Miner (1945) linear damage rule is then utilized to find the total damage in a layer due to cycles along the y 1 (θ) direction as where n a , and n m are the number of bins for the stress amplitude and mean stress, respectively, κ are the cycle counts, and σ i a and σ j m are the stress amplitudes and mean stresses for the layer, respectively. N f is the number of cycles to failure and were computed with a shifted Goodman approach as from the DNVGL standard (DNVGL- ST-0376, 2015). γ is the fatigue factor of safety and was taken to be equal to γ σ because 305 more than two load directions were being analyzed but less than 12.
Everything in the preceding paragraph was also done for the y 1 (θ + 90 • ) direction in parallel. The fatigue damage of the layer then becomes where d is the fatigue damage of the layer due to cycles along the orthogonal y 1 (θ) and y 1 (θ + 90 • ) directions. A fatigue

Objective Function
In defining the objective, or fitness function, it was important to consider that these constraints are all nonlinear and not expressible in closed form with respect to design parameters. Also, because gradient-based optimization was employed in the design process, it was advantageous to define an objective function that is as smooth and continuous in the design space as possible. A constraint aggregation approach has often been used to apply stress constraints to structures in gradient-based 320 optimization while preserving smoothness in the objective (Kreisselmeier and Steinhauser, 1980). With these things in mind, the objective function was defined with the performance constraints enforced as penalty terms, as follows: For material rupture and buckling, the key results are summed, or, in other words, aggregated over the number of elements, n el , and the number of buckling modes analyzed, n modes , respectively, over n θ transverse loading directions, as explained = 60000 and c 2 = 2 were found to be suitable, and c 1 was scaled to match the approximate mass of the blade so that the penalty terms have an effective balance with mass when the constraints become active. c 2 is an exponent in the penalty term, which has a much stronger impact on rate of growth and smoothness. Finally, at any given design state in the optimization process, all the terms in the objective function are evaluated by running FE analyses on the blade under the loads derived from OpenFAST using ANSYS (see Sect. 3 for more details).

Design Variables
The structural optimization in NuMAD was performed by modifying geometric aspects of each blade's design with regard to material thicknesses. Each blade comprises the following main components: spar caps, leading edge panels, and trailing edge panels (each on both the pressure side and the suction side of the airfoil cross section), leading edge reinforcement, trailing edge reinforcement, and two shear webs connecting the suction side and pressure side spar caps. There is also a layer of skin

Algorithms and Other Methodology
As described in Sect. 1, the initial blade designs were received from collaborators at NREL, following the low-fidelity optimization (Bortolotti et al., 2021). Initially, the designs generally did not satisfy the structural performance constraints under the 355 loads determined by the OpenFAST aeroelastic analysis. The high-fidelity structural optimization updated the design to a state in which the performance constraints were satisfied while simultaneously minimizing the mass of the blade. The optimization was performed using the MATLAB built-in function, fmincon, a general nonlinear gradient-based optimization solver. The interior point variation was used, obtaining the gradient/sensitivities of the objective with central differencing.
While the initial designs from the NREL Wind-Plant Integrated System Design & Engineering Model (WISDEM ® ) tool 360 (NREL, 2021b) did not meet all constraints, they were carefully developed based on experience and previous tried-and-true designs. In the interest of staying within reasonably manufacturable parameters and preserving a smooth iteration cycle between high-fidelity and low-fidelity optimization, it was desirable to find an optimized structural design as close to the initial state as possible while still achieving the objectives. Gradient-based optimization adjusts the design variables that most strongly affect the objective, moving to a more favorable state with minimal overall change in design. Once the main automated phase of 365 optimization was completed, some final postprocessing steps were taken, as described in the next section.

Optimization Postprocessing
After the automated structural optimization was completed for each blade design, the results were inspected visually for quality and verification using the ANSYS results viewing window. If needed, minor manual adjustments were made to the design where there could potentially be problems with fabrication or other issues. The model output was examined for factors such as 370 percentage of failed elements and exceeded buckling load factors to evaluate whether any key results were driven by poor mesh quality or fictitious physical effects. The updated design was then sent back to collaborators at NREL for further iteration. If applicable, the new design was evaluated for transportation by rail (BAR-DRG, BAR-DRC, and BAR-URC), and modifications were made as necessary. The global optimization cycle was repeated, beginning with system-level design and optimization until converging on a final design.

Verification
It was first necessary to verify that the shell model coincided with the beam model at NREL. The blade stiffness and total blade mass between each model were verified during one of the design iterations. An identical yaml file was used for the stiffness and mass verification, respectively.

380
Because BeamDyn (NREL, 2021c) is the highest-fidelity beam model from the NREL tools and because loads from Beam-Dyn inform both low-fidelity and high-fidelity optimizers, the deflections between BeamDyn and ANSYS were compared in lieu of direct stiffness comparisons. Since both models make use of a reference line to build the geometry, it was convenient to track its deflection. The displacement of the reference line is easily obtained from the NREL beam analysis because it is part of the formulation. The displacement of the reference line is not easily obtained from a shell (or solid model) since a 385 reference line is not an intrinsic part of the governing differential equations. The resulting displacement field from a shell (or solid) model is a combination of cross-sectional translations, rotations, and warping; an unknown combination at that. Merely using the arithmetic average of displacements from the nodes at a cross section is an approximate approach. Nonetheless, this approach was taken here and thus an approximate comparison to BeamDyn was made. No special effort was made to remove cross-sectional rotations and warping effects from the shell model. Figure 7 shows the flapwise and edgewise deflection due 390 to a uniform flapwise line load; labeled as q 2 (acts parallel to x 2 ). It was applied along the length of the reference axis during a design iteration for BAR-UAG. Figure 7(a) shows a good comparison for the flapwise deflection. The comparison of the edgewise deflections in Fig. 7(b) is not as good as the flapwise deflections but satisfactory nonetheless. Note that the arithmetic mean displacement of a cross section,ū i , acts parallel to x i (see Fig. 2).   During the course of the project, the blade mass between WISDEM and NuMAD for BAR-UAG was also verified. Table   395 4 shows the total mass of each material as well as the total mass of the blade. The codes are within 5 % of each other. The discrepancies were traced primarily to modeling inaccuracies of complex geometries and layup information.

Optimized Blade Designs
The results of the blade design optimization provide insight into understanding the advantages, disadvantages, and potential challenges in different approaches for low-specific-power wind turbine design. The final design states following high-fidelity 400 optimization for all the blade models, with regard to mass and state of performance constraints, are summarized in Table 5. In general, it is shown that the total mass of a blade can be reduced by adopting a flexible design and alleviating the constraint on tip deflection, as for designs BAR-DRG, BAR-DRC, BAR-DRCHT, BAR-URC, and BAR-URCHT. These designs then became driven mainly by material rupture and buckling constraints. In contrast, the upwind designs-BAR-UAG, BAR-USC, and BAR-USCHT-took on a higher final mass and were mainly driven by the constraint of tip deflection.

405
However, the flexible models were seen to generally have a more limited range of feasible design space. For BAR-DRG, for example, both the maximum section stiffness allowable for rail transportability (Bortolotti et al., 2021) and the minimum flap frequency constraint are near their prescribed limits. Because these two constraints are inherently competing, there exists only a small window of feasibility in the spar cap design that satisfies both under the assumptions used.
BAR-DRC has an advantage over BAR-DRG in that the higher stiffness-to-mass ratio of the carbon fiber design brings the 410 flap and edge frequencies up well into the safe range. There are still the material rupture and buckling constraints, however, they remain in opposition of the rail-transportability constraint. Still, if a flexible design is to be pursued, the evidence would support the adoption of a carbon fiber composite spar cap.
None of the models had fatigue damage as a driving factor in the optimization. The nearest model being affected is BAR-DRG, at a maximum damage factor of 0.597. Though on the surface this could be interpreted as a positive result, there may be 415 reason to investigate the analysis approach and further verify it in future work. It is widely known in industry, for example, that when a blade failure arises from fatigue damage, it nearly always nucleates from vulnerable areas of load concentration, such as bolts and joints in panels. These features drive the source and propagation of fatigue damage, yet they are not realistically captured in the shell models currently employed. Even if such features were modeled more realistically, fatigue life inherently has a high degree of uncertainty, and it may turn out that setting/verifying the appropriate safety factors to use for modeling 420 applications like this could be the most productive course of action. Table 6 shows the design load cases, as listed in Sect. 3, which dominated the maximum loads experienced by the blade in the aeroelastic simulations and used in the structural analysis. Although the load distribution applied to the FE model is derived from maximum forces and moments taken throughout the length of the blade and throughout all simulations as previously explained, Table 6 highlights the source of the maximum bending moment at the root in the flap and edgewise directions, as 425 well as the maximum tip deflection. It is evident that the dominant loads generally come from certain case simulations, in  respectively, as these were the blades most governed by material rupture. These same designs also showed buckling as a driving constraint, mainly in leading edge panels and spar caps, about 20 m from the root of the blade. Figures 11, 12, and 13 show the buckling modes nearest to the critical load factor for BAR-DRG, BAR-DRC, and BAR-URC, respectively.  In comparison of the optimized design for blades featuring carbon fiber composite spar caps using industry-standard unidirectional carbon fiber composite vs. heavy-tow carbon fiber composite; only minor differences were seen in the final designs 445 between the two materials. A basic trend was seen that for designs with the heavy-tow composite spar caps, which have a slightly lower ultimate strength, the spar cap thickness generally increased marginally compared to the equivalent designs us-   A detailed cost analysis was performed for each blade design by collaborators at NREL (Bortolotti et al., 2021). The results 450 for blades featuring carbon fiber composite spar caps are summarized in Table 7. Due to the significantly reduced cost of heavytow carbon fiber composite compared to the industry-standard composite, combined with the slightness of change in optimized designs, the analysis indicates roughly a 10 % cost savings in blade manufacture, which is made possible by employing the heavy-tow composite instead of the industry standard. It would seem to be worth investigating this avenue further to verify if the findings are correct, and what additional factors could affect the final cost that may not be accounted for in this study.

Conclusions
Long and flexible blades were designed in a multifidelity optimization collaboration between NREL and SNL. NREL focused on minimizing the levelized cost of energy of each design with a system-level optimization. Since beam FE models of the blades were utilized at the system level, higher structural fidelity was utilized by SNL to minimize mass while simultaneously producing realistic design. Numerous updates to NuMAD were performed, which included an object-oriented approach with 460 new capabilities comprising yaml file compatibility, optimization, and new structural analyses.
Blades that are transportable via controlled flapwise bending were found to be a viable approach to alleviate the logistical issues of oversized blades. For all designs with deflection and fatigue constraints, the deflection constraint was active and fatigue damage was inactive. Other constraints such as buckling and rupture were also found to be active; further confirming the need for structural fidelity greater than a beam FE model. The results confirm that blade mass can, in general, be reduced substantially 465 by going to a downwind design or by increasing nacelle tilt and rotor precone angles above conventional values. The downwind designs were found to be dominated by sudden, abrupt events like gusts, rather than normal operating conditions. If downwind designs are pursued, load-alleviation strategies, such as bend-twist coupling, and active aerodynamics are recommended for future research. The use of the heavy-tow carbon in the spar caps was found to yield a 9-13 % cost savings compared to an industry-standard carbon fiber.