Land-based wind turbines with flexible rail-transportable blades – Part 2: 3D finite element design optimization of the rotor blades
- 1Wind Energy Technologies Department, Sandia National Laboratories, Albuquerque, NM 87185, USA
- 2National Wind Technology Center, National Renewable Energy Laboratory, Golden, CO 80401, USA
Correspondence: Ernesto Camarena (email@example.com)
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.
This paper has been authored by Sandia National Laboratories under contract no. DE-NA0003525 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 non-exclusive, paid-up, 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/doe-public-access-plan, 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 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.
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.
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∘]6 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∘)2]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 and was reported as . Likewise, the arithmetic mean of and was reported as and the arithmetic mean of and was reported as . The approach taken for the shear moduli is not applicable for the Poisson ratios, since in general, , where denotes Poisson ratios. , , and were used directly since , , and 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 [∘]. 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 [∘], 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.
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 σ23 and σ13, respectively. Note that all of these are characteristic layer strengths.
All of the strengths for the uniaxial glass laminate were directly obtained from physical test data. For the remaining glass laminates, test data were only available for the in-plane strengths. The experimental data for the remaining laminates came from the SNL/MSU/DOE Composite Materials for Wind database (SNL/MSU/DOE, 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).
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, (j=1, 2, 3, …, k), with those of the NuMAD loads system, wi, and the ANSYS coordinate system, xi, which are invariant. Note that the xi coordinate system is like the wi system, but the first axis, x1, points toward the leading edge instead of the flap direction.
A small angle between (j=1, 2, 3, …, k) and w3 exists but is unaccounted for in NuMAD. Thus, it is assumed that (j=1, 2, 3, …, k). The remaining angle between 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 Cij denotes the direction cosines that can be stored in a matrix as shown below:
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, ; both bending-moment components, ; and the torsional moment, . The superscript indicates the reference frame. These forces and moments will be referred to collectively with Pi (i=1, 2, …, 4), a generalized load vector, where . Note that Pi 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 Pi, at a given cross section, was defined by
The components were then transformed to the xi 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 were obtained by
Load directions from were, however, obtained by
Note in Eq. (6) that the minimum of 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 y1 in Fig. 3b. The loads in these eight directions were also transformed to the xi 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, Mi acting at a distance, zi, from the blade root and the forces to be solved and Fi acting at , where . The forces to be applied are found by solving the following linear system of equations for Fi:
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 w1 and w2 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 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.
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, utip is the highest tip deflection of the blade in the direction normal to the rotor's sweep plane and umax is the highest allowable blade deflection to ensure tower clearance. fflap and fedge are the first natural vibration frequencies in the flapwise and edgewise directions, respectively, and frot 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.
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 x3. This is done for both the y1(θ) and 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 nsec 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 and 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 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:
where σa and σm are predominant normal stress amplitudes and means, respectively, and Ma and Mm 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 Mr(θ,x3) 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 y1(θ) direction as
where na and nm are the number of bins for the stress amplitude and mean stress, respectively, κ denotes the cycle counts, and and are the stress amplitudes and mean stresses for the layer, respectively. Nf is the number of cycles to failure and was computed with a shifted Goodman approach as
Everything in the preceding paragraph was also done for the 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 y1(θ) and 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 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:
For material rupture and buckling, the key results are summed or, in other words, aggregated over the number of elements, nel, and the number of buckling modes analyzed, nmodes, 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, nsec; the number of blade components, ncomp; 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 c1 and c2 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 c1=60 000 and c2=2 were found to be suitable, and c1 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. c2 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.
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 “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.
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.
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 q2 (acts parallel to x2). 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, , acts parallel to xi (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.
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 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 (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.
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.
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.
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.
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 co-authors 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. DE-NA0003525 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. DE-AC36-08GO28308. 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.
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 Co-located Conferences, American Institute of Aeronautics and Astronautics, https://doi.org/10.2514/6.2011-1880, 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, SAND2012-7028, 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 cross-section analysis and 3D analysis of tapered beams, in: ICCM22 2019, Engineers Australia, Melbourne, VIC, 539–550, 2019. 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.: Land-based wind turbines with flexible rail-transportable blades – Part 1: Conceptual design and aeroservoelastic performance, Wind Energ. Sci., 6, 1277–1290, https://doi.org/10.5194/wes-6-1277-2021, 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/3D-FEM analysis, Multibod. Syst. Dynam., 32, 87–116, https://doi.org/10.1007/s11044-013-9394-3, 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 (ICCM-16), July 2007, Kyoto, Japan, 8–13, 2007. a
Davila, C. G. and Camanho, P. P.: Failure criteria for FRP laminates in plane stress, National Aeronautics and Space Administration, TM-2003-212663, available at: https://ntrs.nasa.gov/citations/20040001363 (last access: 14 January 2022), 2003. a
DNVGL-ST-0376: Rotor blades for wind turbines, https://www.dnv.com/Publications/a-new-dnv-standard-for-rotor-blades-88803 (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. SAND2019-14173, Sandia National Laboratories, https://doi.org/10.2172/1592956, 2019. a, b, c
Fedorov, V.: Bend-twist 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 cross-sectional 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 non-linear behavior of long anisotropic tubes, Int. J. Non-Lin. Mech., 34, 1003–1018, https://doi.org/10.1016/S0020-7462(98)00070-5, 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 fiber-reinforced polymers, Mech. Mater., 142, 103299, https://doi.org/10.1016/j.mechmat.2019.103299, 2020. a
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/B978-0-08-024488-4.50022-X, 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. SAND2019-10780, Sandia National Laboratories, https://doi.org/10.2172/1562792, 2019. a, b
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.: 3-D 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
Yu, W.: Multiscale structural mechanics: Top-down modelling of composites using the structural genome, John Wiley and Sons, Limited, ISBN 9781119092674, 2019. a