Land-based wind turbines with flexible rail-transportable blades – Part 2: 3D finite element 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 field-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 flexible blades that can be shipped on rail lines by flexing during transport. However, the increased flexibility is challenging to accommodate with a typical glass-fiber, upwind design. In a two-part 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 system-level optimization of the rotor variants as compared to conventional and segmented baselines, with a low-fidelity representation of the blades. The present work, Part 2, supplements the system-level optimization in Part 1 with high-fidelity 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 laminate-level 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 heavy-tow carbon fiber in the spar caps is found to yield between 9 % and 13 % cost savings.

1 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⁻²), is diminishing. Using low-specific-power 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 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) have 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 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 US 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 tools and methodologies. After an initial study and down-selection 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.

• 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;

• 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 using industry-standard carbon fiber composite for the spar caps and attached rigidly together by a mechanical joint;

• BAR-URC (upwind – rail transportable – carbon fiber spar caps) – an upwind design using industry-standard 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 low-cost carbon fiber (Ennis et al., 2019), referred to here as heavy-tow (HT) carbon fiber, is also evaluated. Thus, the designs with carbon fiber spar caps have the following additional cases:

• BAR-DRCHT – a variation of the BAR-DRC design, using heavy-tow carbon fiber composite for spar caps instead of the industry standard,

• 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 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 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 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.

In an effort to ensure the BAR designs were neither over nor 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 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 degree-of-freedom 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 object-oriented 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 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 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. Sections 3 and 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 Sect. 6.

2 Material properties

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.

Table 1 Mechanical properties of composite laminates by fiber type. All strength properties are design values (not the characteristic values) and are obtained with γ_σ=1.74.

^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.

Table 2 Mechanical properties of the isotropic constituents. All strength properties are design values (not the characteristic values) and are obtained with γ_σ=1.92. (A larger γ_σ is required by the International Electrotechnical Commission – IEC – IEC-61400-1, 2005, for the foam.)

The uniaxial glass laminate was assumed to comprise Vectorply E-LT 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^∘]₆ PPG-Devold DB810-E05 – a fabric infused with EPIKOTE MGS RIMR 135–EPIKURE RIMH 1366 epoxy resin. The triaxial glass composite was assumed here to comprise [(±45^∘) (0^∘)₂]_s Saertex U14EU920-00940-T1300 and Saertex VU-90079-00830-01270 – fabrics infused with Vantico TDT 177-155.

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 in-plane 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 fiber-level 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 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 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 were used.

Figure 1 Representative volume element utilized to determine the effective elastic properties of the carbon fiber composites.

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 (in-plane) tensile strength, and their primed counterparts are the compressive strengths. S is the in-plane shear strength, and R and T are the shear strengths associated with σ₂₃ and σ₁₃, 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, 2022). 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 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 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).

3 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 are 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 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 (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.

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_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₁, points toward the leading edge instead of the flap direction.

Figure 2 Comparison of the FAST results' coordinate system with that of the NuMAD and ANSYS coordinate systems.

Figure 3 Coordinate systems used to transfer time histories from OpenFAST to the ANSYS model.

A small angle between $v_{\mathrm{3}}^{\left(j\right)}$ (j=1, 2, 3, …, k) and w₃ 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

$$\begin{array}{}\text{(1)}& w_i=C_{iq}{v}_q^{\left(j\right)},\end{array}$$

where C_ij denotes the direction cosines that can be stored in a matrix as shown below:

$$\begin{array}{}\text{(2)}& C=\left[\begin{array}{ccc}\cos \left(\mathit{\mu }\right)& \sin \left(\mathit{\mu }\right)& \mathrm{0}\\ -\sin \left(\mathit{\mu }\right)& \cos \left(\mathit{\mu }\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right].\end{array}$$

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 bending-moment 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)M_{\mathrm{1}}^w\right(t\left)M_{\mathrm{2}}^w\right(t\left)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 tip-deflection 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

$$\begin{array}{}\text{(3)}& P_i={\left\{F_{\mathrm{3}}^w\left(t^{*}\right)M_{\mathrm{1}}^w\left(t^{*}\right)M_{\mathrm{2}}^w\left(t^{*}\right)M_{\mathrm{3}}^w\left(t^{*}\right)\right\}}^T.\end{array}$$

The components were then transformed to the x_i coordinate system.

IEC (IEC-61400-1, 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

$$\begin{array}{}\text{(4)}& P_i=\left\{\begin{array}{l}\max \left(F_{\mathrm{3}}^w\left(t\right)\right)\\ \max \left(M_{\mathrm{1}}^y\left(t\right)\right)\cos \left(\mathit{\theta }\right)\\ \max \left(M_{\mathrm{1}}^y\left(t\right)\right)\sin \left(\mathit{\theta }\right)\\ \max \left(M_{\mathrm{3}}^w\left(t\right)\right)\end{array}\right\}\mathrm{0}{}^{\circ }\le \mathit{\theta }<\mathrm{180}{}^{\circ },\end{array}$$

where

$$\begin{array}{}\text{(5)}& M_{\mathrm{1}}^y=M_{\mathrm{1}}^w\left(t,w_{\mathrm{3}}\right)\cos \left(\mathit{\theta }\right)+M_{\mathrm{2}}^w\left(t,w_{\mathrm{3}}\right)\sin \left(\mathit{\theta }\right).\end{array}$$

Load directions from $\mathrm{180}{}^{\circ }\le \mathit{\theta }<\mathrm{360}{}^{\circ }$ were, however, obtained by

$$\begin{array}{}\text{(6)}& P_i=\left\{\begin{array}{l}\max \left(F_{\mathrm{3}}^w\left(t\right)\right)\\ \min \left(M_{\mathrm{1}}^y\left(t\right)\right)\cos \left(\mathit{\theta }\right)\\ \min \left(M_{\mathrm{1}}^y\left(t\right)\right)\sin \left(\mathit{\theta }\right)\\ \max \left(M_{\mathrm{3}}^w\left(t\right)\right)\end{array}\right\}\mathrm{0}{}^{\circ }\le \mathit{\theta }<\mathrm{180}{}^{\circ }.\end{array}$$

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 bending-moment 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 bending-moment components (i.e., flapwise and edgewise bending) at a span location were projected onto eight directions, as defined by y₁ 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{‾}}{z_i}$, where $\stackrel{\mathrm{‾}}{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:

$$\begin{array}{}\text{(7)}& M_i=\sum _{j=i}^k{F}_j\left({\stackrel{\mathrm{‾}}{z}}_j-z_i\right),\end{array}$$

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 shape-preserving piecewise cubic interpolation. In this way, force distributions parallel to w₁ and w₂ 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.

Figure 4 Free body diagram used to determine transverse loads from a given distribution of moments along the blade span.

4 Main optimization loop

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 high-fidelity 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 system-level perspective.

Figure 5 FE operations' workflow in a single iteration of the optimization.

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:

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }\mathit{Material}\mathit{rupture}.{\mathit{\sigma }}_{\mathit{\theta }}<{\displaystyle \frac{\mathrm{1}}{{\mathit{\gamma }}_{\mathit{\sigma }}}}\Rightarrow {\mathit{\sigma }}_{\mathit{\theta }}{\mathit{\gamma }}_{\mathit{\sigma }}<\mathrm{1}.\text{(9)}& {\displaystyle }& {\displaystyle }\mathit{Fatigue}\mathit{damage}.D<{\displaystyle \frac{\mathrm{1}}{{\mathit{\gamma }}_D}}\Rightarrow D{\mathit{\gamma }}_D<\mathrm{1}.\text{(10)}& {\displaystyle }& {\displaystyle }\mathit{Buckling}.{\mathit{\beta }}_{\mathit{\theta }}>{\mathit{\gamma }}_{\mathit{\beta }}\Rightarrow {\displaystyle \frac{{\mathit{\gamma }}_{\mathit{\beta }}}{{\mathit{\beta }}_{\mathit{\theta }}}}<\mathrm{1}.\text{(11)}& {\displaystyle }& {\displaystyle }\mathit{Maximum}\mathit{tip}\mathit{deflection}.{\mathit{\gamma }}_u{u}_{\mathrm{tip}}<u_{\max }\Rightarrow {\displaystyle \frac{{\mathit{\gamma }}_u{u}_{\mathrm{tip}}}{u_{\max }}}<\mathrm{1}.\text{(12)}& {\displaystyle }& {\displaystyle }\mathit{Minimum}\mathit{flap}\mathit{frequency}.f_{\mathrm{flap}}>{\mathit{\gamma }}_f\left(\mathrm{3}{f}_{\mathrm{rot}}\right)\Rightarrow {\displaystyle \frac{{\mathit{\gamma }}_f\left(\mathrm{3}{f}_{\mathrm{rot}}\right)}{f_{\mathrm{flap}}}}<\mathrm{1}.\end{array}$$

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }\mathit{Minimum}\mathit{edge}\mathit{frequency}.f_{\mathrm{edge}}>{\mathit{\gamma }}_f\left(\mathrm{4}{f}_{\mathrm{rot}}\right)\Rightarrow {\displaystyle \frac{{\mathit{\gamma }}_f\left(\mathrm{4}{f}_{\mathrm{rot}}\right)}{f_{\mathrm{edge}}}}<\mathrm{1}.\text{(14)}& {\displaystyle }& {\displaystyle }\mathit{Flap}\mathit{and}\mathit{edge}\mathit{frequency}\mathit{separation}.f_{\mathrm{edge}}>{\mathit{\gamma }}_f{f}_{\mathrm{flap}}\Rightarrow {\displaystyle \frac{{\mathit{\gamma }}_f{f}_{\mathrm{flap}}}{f_{\mathrm{edge}}}}<\mathrm{1}.\end{array}$$

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 low-fidelity analyses was considered to be a suitable proxy constraint for the present analysis. Thus, this constraint is satisfied here through iteration with the low-fidelity 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.

Table 3 Utilized factors of safety.

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₃. This is done for both the y₁(θ) 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 pseudo-constitutive relations:

$$\begin{array}{}\text{(15)}& {\mathit{\sigma }}_{\mathrm{a}}=\mathit{\alpha }{M}_{\mathrm{a}},{\mathit{\sigma }}_{\mathrm{m}}=\mathit{\alpha }{M}_{\mathrm{m}},\end{array}$$

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

$$\begin{array}{}\text{(16)}& \mathit{\alpha }={\mathit{\sigma }}_{\mathrm{11}}/M_{\mathrm{r}}\left(\mathit{\theta },x_{\mathrm{3}}\right),\end{array}$$

where σ₁₁ is the normal stress in the fiber direction of a given layer and M_r(θ,x₃) 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₁(θ) direction as

$$\begin{array}{}\text{(17)}& d^{\mathit{\theta }}=\sum _{i=\mathrm{1}}^{n_{\mathrm{a}}}\sum _{j=\mathrm{1}}^{n_{\mathrm{m}}}{\displaystyle \frac{\mathit{\kappa }\left({\mathit{\sigma }}_{\mathrm{a}}^i,{\mathit{\sigma }}_{\mathrm{m}}^j\right)}{N_{\mathrm{f}}}},\end{array}$$

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

$$\begin{array}{}\text{(18)}& N_{\mathrm{f}}={\left({\displaystyle \frac{X+\left|X^{\prime }\right|-|\mathrm{2}{\mathit{\sigma }}_{\mathrm{m}}^j{\mathit{\gamma }}_{\mathit{\sigma }}-X+|X^{\prime }\left|\right|}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{a}}^i\mathit{\gamma }}}\right)}^m\end{array}$$

from the DNVGL standard (DNVGL-ST-0376, 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

$$\begin{array}{}\text{(19)}& d=d^{\mathit{\theta }}+d^{\mathit{\theta }+\mathrm{90}{}^{\circ }},\end{array}$$

where d is the fatigue damage of the layer due to cycles along the orthogonal y₁(θ) 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

$$\begin{array}{}\text{(20)}& D_{\mathit{\theta },i}=\max \left(d\right)\left(\mathit{\theta }=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots },n_{\mathit{\theta }}/\mathrm{2}\right).\end{array}$$

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 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 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:

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& L=m_b+c_{\mathrm{1}}\\ & \left(\begin{array}{c}\frac{\mathrm{1}}{n_{\mathit{\theta }}}{\sum }_{\mathit{\theta }=\mathrm{1}}^{n_{\mathit{\theta }}}\left\{\frac{\mathrm{1}}{n_{\mathrm{el}}}{\sum }_{i=\mathrm{1}}^{n_{\mathrm{el}}}{\left({\mathit{\sigma }}_{\mathit{\theta },i}{\mathit{\gamma }}_{\mathit{\sigma }}\right)}^{c_{\mathrm{2}}}+\frac{\mathrm{1}}{n_{\mathrm{modes}}}{\sum }_{i=\mathrm{1}}^{n_{\mathrm{modes}}}{\left(\frac{{\mathit{\gamma }}_{\mathit{\beta }}}{{\mathit{\beta }}_{\mathit{\theta },i}}\right)}^{c_{\mathrm{2}}}\right\}\\ +\frac{\mathrm{2}}{n_{\mathit{\theta }}}\frac{\mathrm{1}}{n_{\mathrm{comp}}}\frac{\mathrm{1}}{n_{\sec }}{\sum }_{\mathit{\theta }=\mathrm{1}}^{\frac{n_{\mathit{\theta }}}{\mathrm{2}}}{\sum }_{i=\mathrm{1}}^{n_{\sec }}{\sum }_{j=\mathrm{1}}^{n_{\mathrm{comp}}}{\left(D_{\mathit{\theta },ij}{\mathit{\gamma }}_D\right)}^{c_{\mathrm{2}}}\\ +{\left(\frac{u_{\mathrm{tip}}}{u_{\max }}\right)}^{c_{\mathrm{2}}}+{\left(\frac{{\mathit{\gamma }}_f\left(\mathrm{3}{f}_{\mathrm{rot}}\right)}{f_{\mathrm{flap}}}\right)}^{c_{\mathrm{2}}}+{\left(\frac{{\mathit{\gamma }}_f\left(\mathrm{4}{f}_{\mathrm{rot}}\right)}{f_{\mathrm{edge}}}\right)}^{c_{\mathrm{2}}}+{\left(\frac{{\mathit{\gamma }}_f{f}_{\mathrm{flap}}}{f_{\mathrm{edge}}}\right)}^{c_{\mathrm{2}}}\end{array}\right).\end{array}\end{array}$$

Table 4 Materials used in each component. LE denotes leading edge; TE denotes trailing edge.

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 one-half 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₁ and c₂ are case-dependent, 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₁=60 000 and c₂=2 were found to be suitable, and c₁ 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₂ 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, 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 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.

Figure 6 Cross-sectional component breakdown of a wind turbine blade.

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 heavy-tow 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 low-fidelity 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 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 and Engineering Model (WISDEM^®) tool (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 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 low-fidelity optimization. Furthermore, the torsional response of the blades was expected to be minimal in comparison to the flexural response.

Figure 7 Flapwise and edgewise static deflection comparison between NuMAD and BeamDyn during a design iteration for BAR-UAG.

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 (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.

5 Results

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 highest-fidelity beam model from the NREL tools and because loads from BeamDyn 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 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 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 cross-sectional 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₂ (acts parallel to x₂). It was applied along the length of the reference axis during a design iteration for BAR-UAG. 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{‾}}{u}}_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 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.

Table 5 Blade mass comparison between WISDEM and NuMAD during a design iteration for BAR-UAG.

Table 6 Final mass and state of design constraints for optimized blade designs. Constraint index (unitless) values expressed as defined in Eqs. (8)–(13), with limiting constraints denoted with bold font.

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 low-specific-power wind turbine design. The final design states following high-fidelity 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 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. 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 cut-out speed and reduced rotor-swept areas due to loading.

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 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 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 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 Driving design load cases and limiting components for optimized blade designs. Included is the margin by which the indicated design load case exceeds the others.

Figure 8 Tsai–Wu failure index for the optimized BAR-DRG design under critical loading, shown from the (a) pressure side, (b) suction side, and (c) cut-away view, showing the leading shear web. Vulnerable areas circled in red.

Figure 9 Tsai–Wu failure index for the optimized BAR-DRC design under critical loading, shown from the (a) pressure side, (b) suction side, and (c) cut-away view, showing the leading shear web. Vulnerable areas circled in red.

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.

Figure 10 Tsai–Wu failure index for the optimized BAR-URC design under critical loading, shown from the (a) pressure side, (b) suction side, and (c) cut-away view, showing the leading shear web. Vulnerable areas circled in red.

Figure 11 Displacement plot of the most critical buckling mode for the optimized BAR-DRG design viewed from the (a) suction side and (b) trailing edge.

Figure 12 Displacement plot of most critical buckling mode for optimized BAR-DRC design viewed from the (a) suction side and (b) trailing edge.

Figure 13 Displacement plot of most critical buckling mode for optimized BAR-URC design viewed from the (a) suction side and (b) trailing edge.

Regarding the individual designs, a clear observation is that, for the downwind designs (BAR-DRG, BAR-DRC, and BAR-DRCHT), 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 aero-devices or bend–twist coupling, it could perhaps allow the blade mass to come down yet further or increase the robustness of the design.

Figure 14 Comparison of optimized spar cap thickness distribution for the BAR-DRC design using industry-standard carbon fiber unidirectional composite and heavy-tow fiber composite.

Figure 15 Comparison of optimized spar cap thickness distribution for the BAR-USC design using industry-standard carbon fiber unidirectional composite and heavy-tow fiber composite.

Figure 16 Comparison of optimized spar cap thickness distribution for BAR-URC design using industry-standard carbon unidirectional composite and heavy-tow fiber composite.

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 trailing-edge 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 BAR-DRG, BAR-DRC, and BAR-URC, 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–13 show the buckling modes nearest to the critical load factor for BAR-DRG, BAR-DRC, and BAR-URC, respectively.

In a 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 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 using the industry-standard carbon fiber composite. Figures 14–16 show the final optimized spar cap thickness profile for BAR-DRC–BAR-DRCHT, BAR-USC–BAR-USCHT, and BAR-URC–BAR-URCHT, 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 heavy-tow 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 saving 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.

6 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 designs. Numerous updates to NuMAD were performed, which included an object-oriented 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 heavy-tow carbon in the spar caps was found to yield a 9 %–13 % cost saving compared to an industry-standard carbon fiber.

Table 8 Cost comparison of final blade designs using industry-standard carbon fiber composite spar caps vs. heavy-tow carbon fiber composite spar caps.

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 co-design.