Research article 19 Jan 2022
Research article  19 Jan 2022
Landbased wind turbines with flexible railtransportable blades – Part 2: 3D finite element design optimization of the rotor blades
 ^{1}Wind Energy Technologies Department, Sandia National Laboratories, Albuquerque, NM 87185, USA
 ^{2}National Wind Technology Center, National Renewable Energy Laboratory, Golden, CO 80401, USA
 ^{1}Wind Energy Technologies Department, Sandia National Laboratories, Albuquerque, NM 87185, USA
 ^{2}National Wind Technology Center, National Renewable Energy Laboratory, Golden, CO 80401, USA
Correspondence: Ernesto Camarena (ecamare@sandia.gov)
Hide author detailsCorrespondence: Ernesto Camarena (ecamare@sandia.gov)
Increasing growth in landbased 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 fieldassembled joints result in added mass and loads, as well as increased reliability concerns in operation. An alternative to this methodology is to design slender flexible blades that can be shipped on rail lines by flexing during transport. However, the increased flexibility is challenging to accommodate with a typical glassfiber, upwind design. In a twopart paper series, several design options are evaluated to enable slender flexible blades: downwind machines, optimized carbon fiber, and active aerodynamic controls. Part 1 presents the systemlevel optimization of the rotor variants as compared to conventional and segmented baselines, with a lowfidelity representation of the blades. The present work, Part 2, supplements the systemlevel optimization in Part 1 with highfidelity blade structural optimization to ensure that the 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 laminatelevel blade optimization and an interface to the International Energy Agency Wind Task 37 blade ontology. Transporting long, flexible blades via controlled flapwise bending is found to be a viable approach for blades of up to 100 m. The results confirm that blade mass can 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, deflection, and resonant frequencies is presented. An analysis of driving load cases revealed that the downwind designs are dominated by loads from sudden, abrupt events like gusts rather than fatigue. Finally, an analysis of carbon fiber spar caps for downwind machines finds that, compared to typical carbon fibers, the use of a new heavytow carbon fiber in the spar caps is found to yield between 9 % and 13 % cost savings.
 Article
(5641 KB)  Companion paper
 BibTeX
 EndNote
This paper has been authored by Sandia National Laboratories under contract no. DENA0003525 with the US Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a nonexclusive, paidup, irrevocable, worldwide license to publish or reproduce the published form of this paper, 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/doepublicaccessplan, last access: 14 January 2022).
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 lowspecificpower turbines has increased wind plant capacity factors, which is a measure of how 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 lowwindspeed sites. Current blades are approaching 80 m in length for landbased installations and over 100 m for offshore turbines. Bolinger et al. (2020) have shown that the ongoing reductions in specific power are likely to continue. While not as much of an issue for offshore wind turbines, landbased machines are currently constrained by transportation logistics. Current landbased transportation constraints limit monolithic blade lengths to about 80 m in length and 4.75 m in width and height. Thus, continued growth in rotor sizes for landbased turbines will likely require design changes.
In 2018, the US Department of Energy (DOE) funded the Big Adaptive Rotor (BAR) project to study the design drivers of future highcapacityfactor landbased 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 tools and methodologies. After an initial study and downselection 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 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.

BARUAG (upwind – air transport – glassfiber spar caps) – a baseline upwind design composed primarily of fiberglass composites with 4 m of blade prebend in a manner similar to current industrystandard designs;

BARDRG (downwind – rail transportable – glassfiber spar caps) – a straight, slender, and downwind design composed primarily of fiberglass composites, intended to be transportable by train on existing railways;

BARDRC (downwind – rail transportable – carbon fiber spar caps) – a straight, slender, and downwind design using industrystandard carbon fiber composite for the spar caps, intended to be transportable by train on existing railways;

BARUSC (upwind – segmented – carbon fiber spar caps) – a segmented upwind design, composed of two sections using industrystandard carbon fiber composite for the spar caps and attached rigidly together by a mechanical joint;

BARURC (upwind – rail transportable – carbon fiber spar caps) – an upwind design using industrystandard carbon fiber composite for the spar caps, intended to be transportable by train on existing railways and installed with 8^{∘} of nacelle tilt and 4^{∘} of rotor precone and no blade prebend.
All blades had a total length, L, of 100 m. Further detail of the BAR concepts and initial design optimization is found in Part 1 (Bortolotti et al., 2021). Since the use of traditional aerospace carbon fibers has contributed limited benefit to wind turbines due to high cost, a new lowcost carbon fiber (Ennis et al., 2019), referred to here as heavytow (HT) carbon fiber, is also evaluated. Thus, the designs with carbon fiber spar caps have the following additional cases:

BARDRCHT – a variation of the BARDRC design, using heavytow carbon fiber composite for spar caps instead of the industry standard,

BARUSCHT – a variation of the BARUSC design, using heavytow carbon fiber composite for spar caps instead of the industry standard,

BARURCHT – a variation of the BARURC design, using heavytow 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 structural fidelities. The National Renewable Energy Laboratory (NREL) led the lowerfidelity 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 twopart 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 wind fields and the interaction of numerous turbine components. These types of models are used at the expense of accurately quantifying the 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 applies to most panels in a given blade. Inadequate 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 wellresolved. Note that linear crosssectional analysis tools (Feil et al., 2020) were utilized in the lowerfidelity optimization. Nonlinear crosssectional 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.
In an effort to ensure the BAR designs were neither over nor underdesigned, the systemlevel optimization was followed by a highfidelity structural optimization. The increase in fidelity required keeping the aerodynamic shape unaltered and relying on the lowerfidelity model for the loading information. The present contribution details the highfidelity 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, landbased rotors are revealed. Solid finite elements (FEs) with a layerwise 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 degreeoffreedom requirement would still be too high for practical design. Shell elements (ANSYS SHELL181) were exclusively utilized in this work because they are 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, with added structural analyses and capability to accept input from the International Energy Agency Wind Task 37 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^{®} gradientbased optimization methods and various calls to ANSYS^{®} to perform meshing and structural analyses. The code and documentation can be obtained from https://github.com/sandialabs/NuMAD (last access: 14 January 2022).
The remainder of this paper summarizes the methods used to derive the material properties needed in the structural analysis and the application of aeroelastic loading to a highfidelity 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. Sections 3 and 4 detail the updates to NuMAD. Section 3 shows how the loading information was transferred from the lowerfidelity aeroelastic simulations to the higherfidelity case. Section 4 describes how the highfidelity structural optimization was performed. Finally, results are presented in Sect. 5 and conclusions are made in Sect. 6.
A common difficulty for the composite 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 were found to have been fully characterized experimentally with elastic, rupture, fatigue, and mass properties (Samborsky et al., 2022). The 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 Tables 1 and 2. Note that γ_{σ} is the factor of safety for material rupture.
^{a} Test result from a smaller fiber volume fraction. (Industry baseline pultrusion tests at 62 % fiber volume fraction; Miller et al., 2019.) ^{b} Uniaxial glass composite value.
The uniaxial glass laminate was assumed to comprise Vectorply ELT 5500 with EPIKOTE MGS RIMR 135–EPIKURE RIMH 1366 epoxy resin from Samborsky et al. (2022). The biaxial glass composite was assumed to comprise [±45^{∘}]_{6} PPGDevold DB810E05 – a fabric infused with EPIKOTE MGS RIMR 135–EPIKURE RIMH 1366 epoxy resin. The triaxial glass composite was assumed here to comprise [(±45^{∘}) (0^{∘})_{2}]_{s} Saertex U14EU92000940T1300 and Saertex VU900790083001270 – fabrics infused with Vantico TDT 177155.
2.1 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 give summaries of how homogenization theory was used for the rest of the composites.
2.1.1 Experimentally determined properties
All nine of the elastic properties of the uniaxial glass laminate were obtained directly from experimental values (Samborsky et al., 2022). Since the tensile and compressive effective elastic moduli, E^{*}, were not equal, the arithmetic mean of the respective tensile and compressive value was utilized. For the shear moduli, G^{*}, the arithmetic mean of ${G}_{\mathrm{12}}^{*}$ and ${G}_{\mathrm{21}}^{*}$ was reported as ${G}_{\mathrm{12}}^{*}$. Likewise, the arithmetic mean of ${G}_{\mathrm{13}}^{*}$ and ${G}_{\mathrm{31}}^{*}$ was reported as ${G}_{\mathrm{13}}^{*}$ and the arithmetic mean of ${G}_{\mathrm{23}}^{*}$ and ${G}_{\mathrm{32}}^{*}$ was reported as ${G}_{\mathrm{23}}^{*}$. The approach taken for the shear moduli is not applicable for the Poisson ratios, since in general, ${\mathit{\nu}}_{ij}^{*}\ne {\mathit{\nu}}_{ji}^{*}$, where ${\mathit{\nu}}_{ij}^{*}$ denotes Poisson ratios. ${\mathit{\nu}}_{\mathrm{12}}^{*}$, ${\mathit{\nu}}_{\mathrm{13}}^{*}$, and ${\mathit{\nu}}_{\mathrm{23}}^{*}$ were used directly since ${\mathit{\nu}}_{\mathrm{21}}^{*}$, ${\mathit{\nu}}_{\mathrm{31}}^{*}$, and ${\mathit{\nu}}_{\mathrm{32}}^{*}$ were obtained by loading the laminate in unconventional ways.
2.1.2 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 homogenized as a stack of homogeneous layers. For such a case, Yu (2012) has rigorously proved that the inplane strains and transverse stresses are uniform throughout the thickness of the laminate. Thus, the effective properties of the stack of layers can be obtained from a hybrid rule of mixtures in an exact manner. Thus, the approach of Yu (2012) 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 [$\mathrm{45}/\mathrm{45}$^{∘}]. The thickness fraction for each layer, defined as the layer thickness divided by the total thickness, was 0.5. Note that actual thicknesses are not required for the analysis. As for the triaxial glass unit cell, it comprised three layers of the uniaxial glass laminate with a stacking sequence of [$\mathrm{45}/\mathrm{0}/\mathrm{45}$^{∘}], where the ±45^{∘} layers had a thickness fraction of 0.25 and the 0^{∘} layer had a 0.5 thickness fraction. All thickness fractions were derived from examining the SNL/MSU/DOE Composite Materials for Wind Database (SNL/MSU/DOE, 2022).
2.1.3 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 fiberlevel details for the homogenization. The microstructure considered here was a hexangular packing of circular fibers, as shown in the repeating unit cell (RUC) of Fig. 1. The local stress and strain fields in this microstructure are not amenable to the simplifications of Sect. 2.1.2. Thus, the effective elastic properties of the carbon fiber composites 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 wellknown 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 industrystandard carbon fiber and 234.69 GPa for the heavytow 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 were used.
2.2 Rupture properties
Table 1 also shows the strength values for each laminate. Note that X is the longitudinal tensile strength, Y is the transverse (inplane) tensile strength, and their primed counterparts are the compressive strengths. S is the inplane shear strength, and R and 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 inplane strengths. The experimental data for the remaining laminates came from the SNL/MSU/DOE Composite Materials for Wind database (SNL/MSU/DOE, 2022). The outofplane shear strength properties, R and T, for the biaxial and triaxial glass laminates, are set equal to the outofplane properties of the 0^{∘} laminate because they 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 inplane shear strength, S, for biaxial and triaxial glass laminates. These values were both determined 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 DNVGLST0376 (2015), whereas for the carbon fiber laminates, they are from Ennis et al. (2019).
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 are representative of the most critical loads from the beam model's time history analysis. The designs identified by the lowfidelity optimization were run through the aeroservoelastic solver OpenFAST (NREL, 2021a). Load equivalency 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. Note that we assume that airfoil shape deviations due to material deformations do not affect the loading on each blade design. The ultimate and fatigue limit states were evaluated from the International Electrotechnical Commission (IEC) standard (IEC614001, 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.
It is important to note that these crosssectional 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}_{i}^{\left(j\right)}$ (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.
A small angle between ${v}_{\mathrm{3}}^{\left(j\right)}$ (j=1, 2, 3, …, k) and w_{3} exists but is unaccounted for in NuMAD. Thus, it is assumed that ${v}_{\mathrm{3}}^{\left(j\right)}={w}_{\mathrm{3}}={x}_{\mathrm{3}}$ (j=1, 2, 3, …, k). The remaining angle between ${v}_{\mathit{\alpha}}^{\left(j\right)}$ 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. 3a and are related by
where C_{ij} denotes 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, and fatigue failures. Both cases utilized the OpenFAST resultant axial forces, ${F}_{\mathrm{3}}^{w}$; both bendingmoment components, ${M}_{\mathit{\alpha}}^{w}$; and the torsional moment, ${M}_{\mathrm{3}}^{w}$. 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 ${P}_{i}=\left[{F}_{\mathrm{3}}^{w}\right(t\left)\phantom{\rule{0.25em}{0ex}}{M}_{\mathrm{1}}^{w}\right(t\left)\phantom{\rule{0.25em}{0ex}}{M}_{\mathrm{2}}^{w}\right(t\left)\phantom{\rule{0.25em}{0ex}}{M}_{\mathrm{3}}^{w}\right(t){]}^{T}$. Note that P_{i} varies with spanwise 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 tipdeflection case, the time of the maximum tip deflection from the beam model, t^{*}, was determined. Then P_{i}, at a given cross section, was defined by
The components were then transformed to the x_{i} coordinate system.
IEC (IEC614001, 2005) allows for lower factors of safety if numerous analysis directions are considered. Eight analysis directions were considered here. Thus, the loads used to evaluate blade failure were obtained by letting θ, as defined in Fig. 3b, vary from 0 to 360^{∘} in load angle increments, Δθ, of 45^{∘}; yielding eight FE load cases. The results for the load directions $\mathrm{0}{}^{\circ}\le \mathit{\theta}<\mathrm{180}{}^{\circ}$ were obtained by
where
Load directions from $\mathrm{180}{}^{\circ}\le \mathit{\theta}<\mathrm{360}{}^{\circ}$ were, however, obtained by
Note in Eq. (6) that the minimum of ${M}_{\mathrm{1}}^{y}$ 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 bendingmoment resultants along the span could possibly come from different times. Thus, it is an artificial distribution suitable for design use. Also, unlike the deflection case, the two bendingmoment components (i.e., flapwise and edgewise bending) at a span location were projected onto eight directions, as defined by y_{1} in Fig. 3b. 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 for evaluating failure, the forces to be applied to the blade were obtained by writing mechanical equilibrium expressions for each spanwise position of interest. Figure 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 at $\stackrel{\mathrm{\u203e}}{{z}_{i}}$, where $\stackrel{\mathrm{\u203e}}{{z}_{i}}=\frac{\mathrm{1}}{\mathrm{2}}({z}_{i}+{z}_{i+\mathrm{1}})$. The forces to be applied are found by solving the following linear system of equations for F_{i}:
where k is the number of cross sections with resultant moment data. Because OpenFAST only outputs these moment data at a maximum of 10 blade cross sections, the moments were with a shapepreserving piecewise cubic interpolation. In this way, force distributions parallel to w_{1} and w_{2} are obtained. Finally, these force distributions were transformed into nodal forces on the blade outer shell by following Berg et al. (2011) but with an extension to apply axial forces and twisting moments.
Once the preprocessing steps were completed for each design, the main structural optimization process was performed in NuMAD. The following sections describe the definition of the optimization objective function, design variables, and algorithms and methods used.
The overall objective of the highfidelity structural optimization in NuMAD was to minimize the mass of the blade structure subject to the following five performance constraints: (1) resistance against material rupture, (2) resistance against fatigue damage over the operating life span of the turbine, (3) resistance against structural buckling, (4) maximum tip deflection that does not exceed the limit for tower clearance, and (5) natural vibration frequencies in the flap and edge directions that do not coincide with the main rotor harmonics. Minimization of mass is ideal for nearly any mechanical device if only for the sake of reducing consumption and the cost of materials for the component itself. For the present application, there is additional importance from a systemlevel perspective.
4.1 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 were first computed. Then the load data that are associated with the maximum tip deflection from OpenFAST were 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 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 criterion that neglects any effects due to transverse normal stress (ANSYS, 2017) (here, transverse normal stress refers to normal stress perpendicular to the laminate plane). 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 of analysis directions and that the default value of −1 in ANSYS was used for the three shear coupling coefficients in the Tsai–Wu criterion. The buckling load factor associated with the lowest mode, β_{θ} (i=1, 2, 3, …, n_{θ}), was then stored from a subsequent eigenvalue buckling analysis of the entire blade. This process is repeated for every load direction. It is notable that geometric and material nonlinearity is unaccounted for. These constraints are summarized mathematically as follows:
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_{flap} 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 listed in Table 3. Note that an explicit rail transport constraint is not implemented in the present analysis. Since this constraint is accounted for in Part 1 (Bortolotti et al., 2021), equating the tip deflections between high and lowfidelity analyses was considered to be a suitable proxy constraint for the present analysis. Thus, this constraint is satisfied here through iteration with the lowfidelity optimization. Also, it should be mentioned that this study does not account for stress concentrations due to ply drops and other abrupt changes in material thickness. In practice, more material and, therefore, mass would be required than what is concluded in this study.
4.1.1 Fatigue damage postprocessor
NuMAD 3.0 now includes the capability to evaluate fatigue damage for every material layer at a discrete number of spanwise stations. First, stress data in the local layer coordinate system for every element and every layer are read into MATLAB^{®}, as well as the element location along x_{3}. This is done for both the y_{1}(θ) and ${y}_{\mathrm{1}}(\mathit{\theta}+\mathrm{90}{}^{\circ})$ directions to consider the additive fatigue damage brought by two orthogonal load directions. The script then loops through each blade cross section (a total of n_{sec} cross sections), comprising the root and each OpenFAST gage location. Because a time history analysis in the shell model was avoided, at each of these locations Markov matrices for both orthogonal directions are constructed from rainflow cycle counting on the ${M}_{\mathrm{1}}^{y}({M}_{\mathrm{1}}^{v},{M}_{\mathrm{2}}^{v},\mathit{\theta})$ and ${M}_{\mathrm{1}}^{y}({M}_{\mathrm{1}}^{v},{M}_{\mathrm{2}}^{v},\mathit{\theta}+\mathrm{90}{}^{\circ})$ time histories (see Fig. 4). Thus, the NREL Crunch (Buhl, 2003) program was set up to analyze each wind speed file with “calculated channels” for ${M}_{\mathrm{1}}^{y}$ at each spanwise location and direction. Because the time history data from OpenFAST comprise 10 min of simulated time, the cycle counts are extrapolated to account for the design life of 30 years by considering a Rayleigh distribution of the annual wind speed. Fatigue damage is then calculated for all layers in elements within 0.75 m of each cross section.
In a given element and a given layer, fatigue damage is calculated by converting the Markov matrix from being defined with moment amplitudes and means to the amplitudes and means of the predominant normal stress component with pseudoconstitutive relations:
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 momentbased Markov matrices are converted to stressbased Markov matrices suitable for damage calculations. The wellknown 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, κ denotes the cycle counts, and ${\mathit{\sigma}}_{\mathrm{a}}^{i}$ and ${\mathit{\sigma}}_{\mathrm{m}}^{j}$ are the stress amplitudes and mean stresses for the layer, respectively. N_{f} is the number of cycles to failure and was computed with a shifted Goodman approach as
from the DNVGL standard (DNVGLST0376, 2015). γ is the fatigue factor of safety and was taken to be equal to γ_{σ} because more than 2 load directions were being analyzed but fewer than 12.
Everything in the preceding paragraph was also done for the ${y}_{\mathrm{1}}(\mathit{\theta}+\mathrm{90}{}^{\circ})$ 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}_{\mathrm{1}}(\mathit{\theta}+\mathrm{90}{}^{\circ})$ directions. A fatigue damage fraction is assigned to the ith cross section by taking the maximum damage of every element within 0.75 m of the cross section and the maximum of each layer in those elements as in
The process repeats for the remaining cross sections. Then, the current direction is incremented by 45^{∘} and the process repeats for ${n}_{\mathit{\theta}}/\mathrm{2}$ times.
4.2 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 gradientbased 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 gradientbased 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}, over n_{θ} transverse loading directions, as explained in Sect. 3. For fatigue damage, the results are summed over the number of blade cross sections analyzed for fatigue damage, n_{sec}; the number of blade components, n_{comp}; and over onehalf of the number of loading directions, as the loading directions are analyzed in pairs for fatigue analysis. D_{θ,ij} represents the maximum damage fraction in cross section i and blade component j, under load direction θ. The constants c_{1} and c_{2} are casedependent, and ideally chosen so that the constraint aggregation terms are negligible for a design that fully satisfies all constraints but become dominant when any constraint is violated. Values of c_{1}=60 000 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 the 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).
4.3 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, leadingedge panels and trailingedge panels (each on both the pressure side and the suction side of the airfoil cross section), leadingedge reinforcement, trailingedge reinforcement, and two shear webs connecting the suction side and pressure side spar caps. There is also a layer of skin on the outer and inner surfaces along the entire spanwise length of the blade. These main components are illustrated in a cross section view in Fig. 6.
The data structure that NuMAD uses to define a blade includes geometric specifications such as width and thickness for each of these components along the span of the blade. The design variables defined for structural optimization in the present study control the spanwise thicknesses of all the components, the width of the spar caps (defined as a single nominal value for the entire length), and the spanwise extent of the shear webs. Note that each component often comprised different materials. Table 4 shows the stack of materials in each component. The thickness of each material was allowed to vary individually. Note that the only difference in the stack definition between the carbon blades listed in Table 4 and the corresponding heavytow designs was the carbon material in the spar caps. Every laminate orientation was fixed such that ${E}_{\mathrm{1}}^{*}$ “pointed toward” the blade tip and remained tangential to the shell reference surface.
The thickness of each component along the span of each blade was defined at 30 equally spaced points from the root to the tip. Six key points along the length were defined as design variables for each component, with the intermediate points interpolated as a smooth spline to define the spanwise distribution.
4.4 Algorithms and other methodology
As described in Sect. 1, the initial blade designs were received from collaborators at NREL, following the lowfidelity optimization (Bortolotti et al., 2021). Initially, the designs generally did not satisfy the structural performance constraints under the loads determined by the OpenFAST aeroelastic analysis. The highfidelity 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 builtin function, fmincon, a general nonlinear gradientbased 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 WindPlant Integrated System Design and Engineering Model (WISDEM^{®}) tool (NREL, 2021b) did not meet all constraints, they were carefully developed based on experience and previous triedandtrue designs. In the interest of staying within reasonably manufacturable parameters and preserving a smooth iteration cycle between highfidelity and lowfidelity optimization, it was desirable to find an optimized structural design as close to the initial state as possible while still achieving the objectives. Gradientbased 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 optimization was completed, some final postprocessing steps were taken, as described in the next section.
The geometry was discretized with SHELL181 elements that were nominally 0.4 m long. The reference surface of the elements coincided with the outer mold line of the blades. Note that when certain shell elements, like SHELL181, are used in this way, they have been observed to overpredict the torsional angle of twist by about 30 % (Branner et al., 2007; Fedorov, 2012). Using SHELL181 elements was deemed appropriate since the angle of twist was not a quantity of interest, nor did it affect the loads since they came from the lowfidelity optimization. Furthermore, the torsional response of the blades was expected to be minimal in comparison to the flexural response.
4.5 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 the 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 (BARDRG, BARDRC, and BARURC), and modifications were made as necessary. The global optimization cycle was repeated, beginning with systemlevel design and optimization until converging on a final design.
5.1 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. Identical .yaml files were used for the stiffness and mass verification.
Because BeamDyn (NREL, 2021c) is the highestfidelity beam model from the NREL tools and because loads from BeamDyn inform both lowfidelity and highfidelity 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 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 crosssectional translations, rotations, and warping and is 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 crosssectional rotations and warping effects from the shell model. Figure 7 shows the flapwise and edgewise deflection due 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 BARUAG. Figure 7a shows a good comparison for the flapwise deflection. The comparison of the edgewise deflections in Fig. 7b is not as good as that of the flapwise deflections but satisfactory nonetheless. Note that the arithmetic mean displacement of a cross section, ${\stackrel{\mathrm{\u203e}}{u}}_{i}$, acts parallel to x_{i} (see Fig. 2).
During the course of the project, the blade mass between WISDEM and NuMAD for BARUAG was also verified. Table 5 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.
5.2 Optimized blade designs
The results of the blade design optimization provide insight into understanding the advantages, disadvantages, and potential challenges in different approaches for lowspecificpower wind turbine design. The final design states following highfidelity optimization for all the blade models, with regard to the mass and state of performance constraints, are summarized in Table 6. 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 BARDRG, BARDRC, BARDRCHT, BARURC, and BARURCHT. These designs then became driven mainly by material rupture and buckling constraints. In contrast, the upwind designs – BARUAG, BARUSC, and BARUSCHT – took on a higher final mass and were mainly driven by the constraint of tip deflection. Although the flexible downwind designs were found to have lower blade mass, Part 1 of this study found that they also produced less power due to decreased cutout speed and reduced rotorswept areas due to loading.
However, the flexible models were seen to generally have a more limited range of feasible design space. For BARDRG, 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.
BARDRC has an advantage over BARDRG in that the higher stiffnesstomass ratio of the carbon fiber design brings the flap and edge frequencies up well into the safe range. There are still the material rupture and buckling constraints; however, they remain in opposition to the railtransportability 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 BARDRG, at a maximum damage factor of 0.597. Though on the surface this could be interpreted as a positive result, there may be reason to investigate the analysis approach and further verify it in future work. It is widely known in the 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 applications like this could be the most productive course of action.
Table 7 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 7 highlights the source of the maximum bending moment at the root in the flap and edgewise directions, as well as the maximum tip deflection. It is evident that the dominant loads generally come from certain case simulations, in particular 1.4, with some from 1.3 from the IEC standard design load cases. This trend is not universal but is consistent among the examined designs.
Regarding the individual designs, a clear observation is that, for the downwind designs (BARDRG, BARDRC, and BARDRCHT), loads are distinctly dominated by design load case 1.4, whereas the upwind design loads are more balanced between simulation cases. Although DLC 1.4 is prominent for all designs, when it dominates in upwind designs, it does so by a notably smaller margin compared to downwind designs. The implication is that downwind designs will generally be more specifically driven by sudden, abrupt events like gusts, whereas upwind designs would be somewhat more dictated by normal operation. Furthermore, the dominance of DLC 1.4 in downwind designs may suggest that if some measures were taken to alleviate loads in sudden gusts, such as active aerodevices or bend–twist coupling, it could perhaps allow the blade mass to come down yet further or increase the robustness of the design.
Table 7 also indicates any components of the blade that are limiting with regard to critical driving constraints. For blades driven by material rupture, the shear web skin and trailingedge reinforcement commonly showed peak failure indices, particularly near the root of the blade where the shear webs begin and near relatively sudden changes in the cross section. Figures 8–10 show the Tsai–Wu failure index distribution for the critical loading for BARDRG, BARDRC, and BARURC, respectively, as these were the blades most governed by material rupture. These same designs also showed buckling as a driving constraint, mainly in leadingedge panels and spar caps, about 20 m from the root of the blade. Figures 11–13 show the buckling modes nearest to the critical load factor for BARDRG, BARDRC, and BARURC, respectively.
In a comparison of the optimized design for blades featuring carbon fiber composite spar caps using industrystandard unidirectional carbon fiber composite vs. heavytow carbon fiber composite, only minor differences were seen in the final designs between the two materials. A basic trend was seen that for designs with the heavytow composite spar caps, which have a slightly lower ultimate strength, the spar cap thickness generally increased marginally compared to the equivalent designs using the industrystandard carbon fiber composite. Figures 14–16 show the final optimized spar cap thickness profile for BARDRC–BARDRCHT, BARUSC–BARUSCHT, and BARURC–BARURCHT, respectively.
A detailed cost analysis was performed for each blade design by collaborators at NREL (Bortolotti et al., 2021). The results for blades featuring carbon fiber composite spar caps are summarized in Table 8. Due to the significantly reduced cost of heavytow carbon fiber composite compared to the industrystandard composite, combined with the slightness of change in optimized designs, the analysis indicates roughly a 10 % cost saving in blade manufacture, which is made possible by employing the heavytow 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.
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 systemlevel 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 designs. Numerous updates to NuMAD were performed, which included an objectoriented approach with 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 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. The use of the heavytow carbon in the spar caps was found to yield a 9 %–13 % cost saving compared to an industrystandard carbon fiber.
Since it is important to ascertain how each blade will interact with the rest of the turbine, a more detailed aeroelastic analysis will be performed in a future work. The focus there will be on system dynamics and aeroelastic instabilities in the frequency domain. Furthermore, unlike the present study, which had one tower for all designs, the tower will be redesigned for each blade design to explore tradeoffs and opportunities for cost savings. The impact of the location of the center of gravity of the rotor–nacelle assembly will be investigated, especially for downwind rotors where aerodynamic thrust and rotor gravity loads align and generate larger moments along the tower compared to equivalent upwind configurations. Better tuning of the controller will also be achieved through a controls codesign.
The NuMAD code as well the .yaml files and other inputs is publicly available at https://doi.org/10.5281/zenodo.5851606 (Camarena et al., 2022).
JP and NJ performed conceptualization, funding acquisition, project administration, and supervision. EC and EA implemented the new NuMAD methodology and prepared the original draft. PB and RF assisted with the BeamDyn investigation, and PB assisted with the conceptualization and cost analysis. All authors performed editing and review.
The contact author has declared that neither they nor their coauthors have any competing interests.
The views expressed in the article do not necessarily represent the views of the DOE or the US Government.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the US Department of Energy's National Nuclear Security Administration. The authors would like to acknowledge Chris Kelley, Brandon Ennis, and Ryan Clarke for their meaningful technical contributions to the project.
This manuscript has been authored by Sandia National Laboratories under contract no. DENA0003525 with the US Department of Energy (DOE) and by the National Renewable Energy Laboratory, operated by the Alliance for Sustainable Energy, LLC, for the US DOE under contract no. DEAC3608GO28308. Funding was provided by the US DOE Office of Energy Efficiency and Renewable Energy Wind Energy Technologies Office.
This paper was edited by Sandrine Aubrun and reviewed by two anonymous referees.
ANSYS: Mechanical, Release 18.1, ANSYS Documentation, ANSYS, Inc., available at: https://www.ansys.com/ (last access: 14 January 2022), 2017. a
Berg, J., Paquette, J., and Resor, B.: Mapping of 1D beam loads to the 3D wind blade for buckling analysis, Structures, Structural Dynamics, and Materials and Colocated Conferences, American Institute of Aeronautics and Astronautics, https://doi.org/10.2514/6.20111880, 2011. a
Berg, J. C. and Resor, B. R.: Numerical Manufacturing And Design Tool (NuMAD v2.0) for Wind Turbine Blades: User's Guide, Sandia National Laboratories, SAND20127028, https://doi.org/10.2172/1051715, 2012. a
Bertolini, P., Sarhadi, A., Stolpe, M., Eder, M. A., and Power, L. W.: Comparison of stress distributions between numerical crosssection analysis and 3D analysis of tapered beams, in: ICCM22 2019, Engineers Australia, Melbourne, VIC, 539–550, 2019. a
Bir, G. S.: Computerized method for preliminary structural design of composite wind turbine blades, J. Solar Energ. Eng., 123, 372–381, https://doi.org/10.1115/1.1413217, 2001. a
Bir, G. S.: User's guide to PreComp (PreProcessor for Computing Composite Blade Properties), National Renewable Energy Laboratory, Technical Report, https://doi.org/10.2172/876556, 2006. a
Bolinger, M., Lantz, E., Wiser, R., Hoen, B., Rand, J., and Hammond, R.: Opportunities for and challenges to further reductions in the “specific power” rating of wind turbines installed in the United States, Wind Eng., 45, 0309524X19901012, https://doi.org/10.1177/0309524X19901012, 2020. a, b
Bortolotti, P., Gaertner, E., and Dykes, K.: System modeling frameworks for wind turbines and plants, in: Wind Energy Systems Engineering Workshop, NREL, available at: https://www.nrel.gov/docs/fy20osti/75658.pdf (last access: 14 January 2022), 2019. a
Bortolotti, P., Johnson, N., Abbas, N. J., Anderson, E., Camarena, E., and Paquette, J.: Landbased wind turbines with flexible railtransportable blades – Part 1: Conceptual design and aeroservoelastic performance, Wind Energ. Sci., 6, 1277–1290, https://doi.org/10.5194/wes612772021, 2021. a, b, c, d, e, f
Bottasso, C. L., Campagnolo, F., Croce, A., Dilli, S., Gualdoni, F., and Nielsen, M. B.: Structural optimization of wind turbine rotor blades by multilevel sectional/multibody/3DFEM analysis, Multibod. Syst. Dynam., 32, 87–116, https://doi.org/10.1007/s1104401393943, 2014. a
Branner, K., Berring, P., Berggreen, C., and Knudsen, H. W.: Torsional performance of wind turbine blades – Part II: Numerical validation, in: International conference on composite materials (ICCM16), July 2007, Kyoto, Japan, 8–13, 2007. a
Buhl, M. L.: Crunch user's guide, National Renewable Energy Laboratory, NREL/EL50030122, available at: https://www.nrel.gov/docs/gen/fy02/30122.pdf (last access: 14 January 2022), 2003. a
Camarena, E., Anderson, E., Ruehl, K., Clarke, R. J., Paquette, J., and Ennis, B. L.: NuMAD v3.0 (3.0), Zenodo [code], https://doi.org/10.5281/zenodo.5851606, 2022. a
Carron, W. S. and Bortolotti, P.: Innovative rail transport of a supersized landbased wind turbine blade, J. Phys.: Conf. Ser., 1618, 042041, https://doi.org/10.1088/17426596/1618/4/042041, 2020. a
Couturier, P. J. and Krenk, S.: Nonlinear collapse of general thinwalled crosssections under pure bending, AIAA J., 54, 3956–3964, https://doi.org/10.2514/1.J054838, 2016. a
Davila, C. G. and Camanho, P. P.: Failure criteria for FRP laminates in plane stress, National Aeronautics and Space Administration, TM2003212663, available at: https://ntrs.nasa.gov/citations/20040001363 (last access: 14 January 2022), 2003. a
DNVGLST0376: Rotor blades for wind turbines, https://www.dnv.com/Publications/anewdnvstandardforrotorblades88803 (last access: 14 January 2022), 2015. a, b
Ennis, B. L., Kelley, C. L., Naughton, B. T., Norris, R. E., Das, S., Lee, D. L., and Miller, D. A.: Optimized carbon fiber composites in wind turbine blade design, Tech. Rep. SAND201914173, Sandia National Laboratories, https://doi.org/10.2172/1592956, 2019. a, b, c
Fedorov, V.: Bendtwist coupling effects in wind turbine blades, DTU Wind Energy, Denmark, available at: https://backend.orbit.dtu.dk/ws/portalfiles/portal/54637711/PhD_Thesis_Vladimir_Fedorov.pdf (last access: 14 January 2022), 2012. a
Feil, R., Pflumm, T., Bortolotti, P., and Morandini, M.: A crosssectional aeroelastic analysis and structural optimization tool for slender composite structures, Composite Struct., 253, 112755, https://doi.org/10.1016/j.compstruct.2020.112755, 2020. a
Harursampath, D. and Hodges, D. H.: Asymptotic analysis of the nonlinear behavior of long anisotropic tubes, Int. J. NonLin. Mech., 34, 1003–1018, https://doi.org/10.1016/S00207462(98)000705, 1999. a
Herráez, M., Bergan, A. C., Lopes, C. S., and González, C.: Computational micromechanics model for the analysis of fiber kinking in unidirectional fiberreinforced polymers, Mech. Mater., 142, 103299, https://doi.org/10.1016/j.mechmat.2019.103299, 2020. a
IEC614001: Wind turbines – Part 1: Design requirements, International Electrotechnical Commission, Geneva, Switzerland, 2005. a, b, c
Johnson, N., Bortolotti, P., Dykes, K. L., Barter, G. E., Moriarty, P. J., Carron, W. S., Wendt, F. F., Veers, P. S., Paquette, J., Kelly, C., and Ennis, B.: Investigation of Innovative Rotor Concepts for the Big Adaptive Rotor Project, https://doi.org/10.2172/1563139, 2019. a
Kreisselmeier, G. and Steinhauser, R.: Systematic Control Design by Optimizing a Vector Performance Indicator, Pergamon, 113–117, https://doi.org/10.1016/B9780080244884.50022X, 1980. a
Miller, D. A., Samborsky, D. D., and Ennis, B. L.: Mechanical testing summary: Optimized carbon fiber composites in wind turbine blade design, Tech. Rep. SAND201910780, Sandia National Laboratories, https://doi.org/10.2172/1562792, 2019. a, b
Ning, A. and Petch, D.: Integrated design of downwind landbased wind turbines using analytic gradients, Wind Energy, 19, 2137–2152, https://doi.org/10.1002/we.1972, 2016. a
NREL: OpenFAST, v2.5.0, available at: https://github.com/OpenFAST/openfast (last access: 14 January 2022), 2021a. a
NREL: WISDEM, v3.2.0, available at: https://github.com/WISDEM/WISDEM (last access: 14 January 2022), 2021b. a
NREL: BeamDyn, available at: https://openfast.readthedocs.io/en/dev/source/user/beamdyn/index.html (last access: 14 January 2022), 2021c. a
Samborsky, D. D., Mandell, J. F., and Agastra, P.: 3D static elastic constants and strength properties of a glass/epoxy unidirectional laminate, available at: https://www.montana.edu/composites/documents/3D Static Property Report.pdf, last access: 14 January 2022. a, b, c
SNL/MSU/DOE: Composite material fatigue database, Version 29.0, available at: https://www.montana.edu/composites/, last access: 14 January 2022. a, b
Yu, W.: An exact solution for micromechanical analysis of periodically layered composites, Mech. Res. Commun., 46, 71–75, https://doi.org/10.1016/j.mechrescom.2012.08.007, 2012. a, b
Yu, W.: Multiscale structural mechanics: Topdown modelling of composites using the structural genome, John Wiley and Sons, Limited, ISBN 9781119092674, 2019. a