1. Introduction
In the design phase of a wind turbine blade, manufacturers often deal with the situation that their new blade design withstands extreme operational loads but fails to withstand loads due to survival wind, when in parked or idling mode. Parked or idling rotors, experiencing extreme wind speeds, are likely to encounter high angles of attack within the post-stall region [
1]. Whether these high angles of attack will appear depends on the inclination of the topography, the tilt angle of the nacelle, as well as the yaw misalignment of the inflow, static or dynamic (due to the wind turbulence). Such high angles of attack may give rise to stall-induced, edgewise vibrations on the blades, which, on many occasions, drive the design loads. These load-driving design load cases (DLCs) are designated in the IEC 61400-1 [
2] standard by the code name DLC-6.x. They include parked or/and idling cases at extreme wind speeds, which are either combined (DLC-6.2) or not (DLC-6.1), with some failure of the grid or malfunctioning of the control system. A typical fault for wind turbines is the loss of the grid connection, combined with a failure in the Uninterrupted Power Supply (UPS) system. In this circumstance, the rotor fails to track the wind and, therefore, high yaw misalignment angles and high angles of attack in deep stall may occur. A similar condition can be encountered during the installation of the turbine before its electrification.
Aeroelastic simulations of idling rotors indicate that maximum blade loads (combined bending moments) appear for yaw angles in the ranges (−40°, −15°) or (15°, 40°) [
1]. Yaw misalignment angles of −/+15° are likely to occur even in normal idling operation (DLC-6.1), while higher yaw angles can only be encountered as a result of some fault in the wind-tracking system (DLC-6.2). Aerodynamic loads in such conditions (i.e., in deep stall) can be accurately estimated using engineering tools that rely on tabulated airfoil data (i.e., the so-called polars), only if a valid dynamic stall model is employed. Nevertheless, the uncertainty in the prediction of the aerodynamic loads in dynamic stall is a well-known fact within the wind energy sector [
3,
4,
5] and it is the main cause of the consequent uncertainty in the prediction of the stall-induced edgewise vibrations. Different state-of-the-art dynamic stall models (e.g., Beddoes–Leishman [
6] and ONERA [
7]) can give significantly different load results in the onset of dynamic stall due to the quite different aerodynamic damping predicted by the stall models. Therefore, extreme loads of an idling rotor might substantially deviate, especially when the overall aerodynamic damping of certain modes of the system is very low or even negative. This explains why a blade that has been designed using a specific dynamic stall model may appear not to be able to withstand the loads of DLCs-6.x when a different dynamic stall model is considered. However, it is not easy to identify which model is the most appropriate for idling rotor analyses or which model provides the most conservative predictions. Dynamic stall measurements at very high angles of attack that could be used as a means to calibrate engineering state-of-the-art dynamic stall models are scarce.
One possible remedy for mitigating the above-discussed uncertainty in load prediction is to tailor the blade design in such a way that the damping of the poorly damped edgewise modes is enhanced. The results of the different dynamic stall models tend to converge when the damping of the rotor modes is well above the onset of the instability. An effective way to enhance the damping of the poorly damped edgewise modes is by means of Flap Edge Coupling (FEC). Previous studies [
8] have shown that, under certain conditions, the coupling of the edgewise and flapwise motion, in the low-damped edgewise modes, increases their aerodynamic damping. The most effective way to enhance FEC is to increase the blade’s structural twist angle, which is achieved by enhancing the cross-bending stiffness along the blade span.
When manufacturers come across the abovementioned problem, i.e., that DLC-6.x loads exceed the strength limits of their design, they usually apply ad hoc, local modifications, e.g., locally reinforcing the inner-blade structure in order to increase strength margins. Such modifications, although solving the problem, increase the weight of the blade, which is by no means optimum. Even if they decide to increase the blade’s structural twist in order to enhance the damping of the edgewise rotor modes, this is also usually done non-optimally, through local interventions in the inner structure of the blade, e.g., by increasing the thickness of the composite material towards the leading edge on the suction side and towards the trailing edge on the pressure side. However, an unbounded increase in FEC could entail implications in the operational loads. Blades with enhanced FEC, under the action of extreme flapwise loads, usually tend to undergo high leading edgewise deflections, which make the blade behave as a virtually forward swept blade. As a result of the forward sweep, the blade tends to twist towards stall and, therefore, to further increase the flapwise loading in normal operation.
A consistent and effective way to tailor the wind turbine blades is to determine the parameters in the design problem within the course of an optimization process. Within such an optimization loop, the best combination of design variables is sought that mitigates stall-induced vibrations without deteriorating normal operation loading, using, as objective functions, some design cost parameter. This is the approach adopted in the present work. Two different FEC methods are employed for the idling DTU-10MW Reference Wind Turbine (RWT) [
9] in order to mitigate stall-induced edgewise vibrations. The first FEC method considered is the one typically employed by wind turbine manufacturers. The thicknesses of the blade section walls are varied asymmetrically between the suction and the pressure side. Higher thicknesses of the uniaxial (UNIAX) and triaxial (TRIAX) composite materials is considered towards the leading edge (LE) of the section on its suction side (segment “leading” of L/P side in
Figure 1, left), while the thickness of the walls on the pressure side is reinforced towards the trailing edge (TE) (segment “trailing” of H/P side in
Figure 1, left). The second method is first addressed by the authors in [
10] and it consists of shifting the two spar caps in opposite directions. The spar cap on the suction side is displaced towards the LE, while the spar cap on the pressure side is displaced towards the TE (as illustrated in
Figure 1, right).
The design parameters in both methods are optimized within a MDAO loop, in which the objective function is the overall
LCoE of the wind turbine. The assessment of the potential design solutions considers both the critical survival wind load cases but also the most critical operational load cases, with the aim to ensure that design modifications that improve the loading of the idling rotor do not entail a deterioration in the operational loading characteristics. The last step of the optimization analyses is to combine both FEC methods and, at the same time, apply a moderate material Bend Twist Coupling (BTC) [
11,
12] through the introduction of an offset angle to the composite plies of the uniaxial (UNIAX) material over the spar caps of the blade (as illustrated in
Figure 2), aiming at compensating any increase in the operational loads due to FEC. The constraint in the optimization process is the requirement that the maximum obtained Tsai–Hill criterion values in the optimum design [
13] do not exceed those of the reference blade considering only operational ultimate loads. This is a deliberate choice, which entails a more conservative design than the one that satisfies the requirement that Tsai–Hill values do not exceed the unit failure threshold. However, it is a choice that ensures that the optimum blade design is directly comparable to the reference blade.
The MDAO framework used in the present study consists of the following in-house or freely available modules:
- -
The in-house hydro-servo-aero-elastic tool hGAST [
14,
15], which simulates the full wind turbine system at all operational conditions of interest and provides ultimate resultant internal loads (reaction forces and moments over the various components in the wind turbine) and power yield.
- -
A cross-sectional analysis tool [
10] based on thin lamination theory [
16], which provides stress distributions and values of the Tsai–Hill failure criterion over the various cross-sections in the structure. The cross-sectional tool takes, as input, the resultant ultimate internal loads generated by hGAST. The same module calculates beam cross-sectional properties in the different components, considering, as input, the material structural properties, thickness and distribution of stresses over every cross-section.
- -
A cost model that determines the cost of the full wind turbine system based on existing data available in the literature [
17], as well as on empirical formulae, which specify certain cost parameters as a function of the component dimensions.
- -
An optimization framework based on ready-made functions from the publicly available scipy library in Python [
18].
The contribution of the present paper lies in the design of wind turbine blades within an MDAO loop that are optimized for LCoE and account for the design-driving loads due to storm conditions. The aero-elastically tailored rotor blades combine both FEC and BTC passive control techniques, the first for enhancing stability in idling operation at survival wind speeds and the second for reducing loads in normal operation. Furthermore, two approaches to introduce FEC (i.e., material and geometrically based) are considered, separately and together. To effectively determine the optimum design parameters in all the adopted control techniques, both normal operation and idling cases are considered within the optimization process.
The paper is structured as follows.
Section 2 briefly describes the above-listed numerical tools and further details the optimization methodology followed.
Section 3 presents and discusses the results of the various optimization analyses performed and, finally,
Section 4 provides some overall conclusions comparing and assessing the different passive load control methods employed in the work.
2. Numerical Tools and Methodology
In the present section, the numerical tools constituting the multidisciplinary aeroelastic optimization framework used in the present work are briefly outlined. Relevant references, in which there are detailed description and validation of the models used by the optimization framework, are provided in the subsequent sections.
2.1. Multibody Servo-Aero-Elastic Analysis Tool hGAST
Nonlinear time domain aero-elastic simulations of the full wind turbine are performed using NTUA’s in-house servo-aero-elastic hGAST [
11,
14,
15]. The solver has been extensively used in numerous (national and EU-funded) research projects over the last 30 years, as well as by the industry. Its predictions for offshore applications have been verified through code-to-code comparison in the OC4 annex [
19,
20]. hGAST is formulated in the context of multibody method. The full wind turbine configuration is configured as an assembly of independent components, i.e., the blades, the drive train and the tower. The above components are modelled as linear Timoshenko beam elements, which are inter-connected through appropriate kinematic and dynamic constraint conditions. The multibody formulation is also extended to the component level. Highly flexible components, such as the wind turbine blades, are further divided into a number of rigidly interconnected “sub-bodies”. Large deflections and rotations gradually build up and nonlinear dynamics are introduced by imposing to each sub-body the deflections and rotations of preceding sub-bodies as rigid-body nonlinear motions. This approach allows for capturing the geometrical nonlinear effects due to large deflections and rotations using linear beam theory at the element level, but considering nonlinear effects at the sub-body level.
Although several aerodynamic options are available in hGAST [
15] (i.e., Blade Element Momentum Theory-BEMT, free wake vortex and finite volume CFD), in the present work, the BEMT approach is adopted, which is the only feasible in terms of computational cost for optimization analyses. The BEMT model used in hGAST is an elaborate version of the standard implementation in which dynamics inflow, yaw misalignment and sectional unsteady aerodynamics and dynamic stall effect are taken into account. Two stall models are available in hGAST, the ONERA [
7] and the Beddoes–Leishman [
6] model. A separate module that provides turbulent inflow boxes is also available, formulated on the basis of Veers’ model [
21]. For the calculation of the turbulent inflow wind time series either the von Karman or the Kaimal spectra, specified in the IEC 61400-1 standard, can be used.
2.2. Cross-Sectional Analysis Tool
An in-house cross-sectional analysis tool, based on thin lamination theory [
16], is employed for the calculation of the distribution of shear stresses over the cross-sections of the different components of the wind turbine (e.g., blades and tower) as well as for the calculation of Tsai–Hill failure criterion [
13]. The cross-sectional tool uses as input for the stress analyses, ultimate resultant loads with respect to a reference point on the cross-section, provided by the time domain aeroelastic analyses with hGAST. The same tool is also used for the calculation of the equivalent beam structural characteristics of the cross-sections, which become input to hGAST for the time domain aeroelastic analyses of the design load cases. It calculates equivalent beam inertia (e.g., mass and moments of inertia) and stiffness (e.g., bending and torsional) characteristics of hollow tubular sections with laminated composite skins and shear webs, while it is capable of treating the anisotropic behaviour of material plies (as shown in
Figure 1 and
Figure 2). Therefore, it can provide fully populated stiffness matrices, taking into account all material-driven coupling effects and it can be used for the modelling of the FEC and BTC blades dealt with in the present work. In a recent upgrade of the model [
10], the possibility to include up to three shear webs within a cross-section is added as an extra modelling feature (necessary for the modelling of the DTU-10 MW Reference blade, which comprises a third small TE shear web). In the same work, a more consistent computation of the shear rigidity properties is performed through the definition of appropriate Timoshenko shear factors, depending on the geometric characteristics in the section. The reader is directed to [
10] for detailed description and verification of the cross-sectional tool.
2.3. Wind Turbine Cost Model
A cost model of the full wind turbine was set up by combining models from the literature that estimate the cost of the materials in the various components along with the cost of the manufacturing processes. The detailed cost model by NREL [
22] is adopted for the blades, while the cost of the other parts of the rotor (i.e., hub, spinner and pitch mechanism) is based on empirical expressions provided in [
23]. The cost estimation of the generator, the gearbox and the tower is derived through up-scaling [
24] the data of a reference wind turbine of nominal power of 1.5 MW, which is described in detail in [
23]. For detailed description and verification of the cost model, the reader is directed to [
17].
2.4. Optimization Framework
The definition of the general optimization problem is:
Minimization of the objective function
with design variables
and fixed parameters
, subject to the following geometric constraints.
and to the following loading constraints
The scalar objective function is the LCoE of the entire turbine that requires as input the material distribution and mass of the turbine components as well as the annual energy production (AEP). The objective function depends on a number of structural design variables, , (e.g., thickness of skin and webs, position and angle of the shear webs, fibre angle, etc.) and aerodynamic design variables, , (e.g., chord, twist, relative thickness and sweep distributions), a set of fixed structural parameters, , (e.g., lay-up sequence on different regions on the section) and aerodynamic parameters, , (e.g., airfoil shapes and lift-drag polars selected from a database).
During the optimization loop, the following types of constraints are satisfied, as illustrated in
Figure 3:
- (i)
Geometric constraints that concern (a) the outer geometry of the blade and directly affect the aerodynamics, , (e.g., maximum chord of the blade or maximum/minimum relative thickness), (b) the inner-blade structure, , (e.g., maximum plausible displacement of the spar caps or the offset ply angle of the UD material over the spar caps) and, finally, (c) the overall turbine characteristics, , (e.g., maximum deflection of the blade tip or maximum shift of the blade natural frequencies).
- (ii)
Loading constraints that are (a) maximum stresses, , along the blade and (b) overall turbine loads, , as, for example, maximum rotor thrust or maximum blade root bending moment.
Communication interfaces are established between the different modules in the optimization loop. An inner-structure geometry parameterization routine is integrated into the cross-sectional analysis tool. This routine is responsible for defining the inner shape of the blade through the global set of design and fixed variables and . A global, spanwise parameterization of the geometry is defined on the basis of Bezier curves for inner-structure parameters, such as the thickness of the skin walls or the position/orientation of the shear webs. A similar parameterization routine is defined for the representation of the external blade shape on the basis of and using Bezier curves.
In the present work, the objective of the optimization framework is to minimize the LCoE of the reference wind turbine through the optimal use of different passive load control techniques (i.e., FEC and BTC) applied on the rotor blades. Apart from the rotor, the remaining wind turbine components remain unchanged, while only loading constraints concerning the maximum stresses along the blades are considered.
The
LCoE is given by:
where
ICC is the initial capital cost provided by the cost model,
BoP is a global approximation of the balance of plant, considered equal to 281 USD/kW of installed capacity,
OPEX denotes the operating expenses, considered equal to 5% of
ICC and
AEP is the annual energy production, calculated by the aeroelastic tool. The
LCoE in the present analyses is calculated for a life-time of the investment of
N = 20 years and a discount rate of
i = 6%.
The
AEP is calculated through a series of time domain aeroelastic simulations with hGAST (see
Section 2.1) over the whole range of the operating wind velocities (from the cut in to the cut out). Additional time domain simulations of selected load-driving DLCs of the IEC standard, including both normal operation and parked-idling operation cases, provide ultimate design loads for the candidate optimum solutions. The ultimate loads are then transformed to stresses and Tsai–Hill criterion values (equivalent stresses) by the thin lamination tool (see
Section 2.2) over the cross-sections in every component (e.g., blade and tower). The imposed constraint to the optimization process requires that the Tsai–Hill criterion values do not exceed the corresponding values of the reference blade. In this way, the thickness of the component walls and the design parameters of the passive load control methods tested by the optimizer are determined. Knowing the dimensions and the material distribution in the load carrying elements of the components the overall capital expenditure (CAPEX) of the candidate solution can be estimated using the cost model (see
Section 2.3).
The minimization problem is addressed through the free-gradient optimizer ‘COBYLA’, which is publicly available as Python module of the scipy library [
18] and is divided into two loops with the aim to moderate the computational cost of the optimization loop [
10,
12]. The outer loop specifies the parameters of FEC and BTC. In the present work, this includes the definition (i) of the thickness distribution along the blade span of the L/P “leading” and H/P “trailing” segments of the cross-sections in the material FEC method (see
Figure 1, left), (ii) of the displacement of the caps in the geometrical FEC method (see
Figure 1, right) and (iii) of the ply offset angle of the material BTC blades (see
Figure 2). Furthermore, the outer loop adjusts the geometric characteristics of the blade (blade planform, i.e., the twist angle distribution) and evaluates the cost function (
LCoE). The inner loop specifies an overall, uniform thickness coefficient for the walls of every cross-section that varies along the span. Its purpose is to estimate the thickness of the blade walls in such a way as to maintain the maximum values of Tsai–Hill criterion of the reference blade (constraint of the optimization process). This part of the optimization process is the most time consuming because it includes the computation of the ultimate resultant loads with hGAST, as well as the calculation of the beam properties and stress distributions of the modified blades based on the FEC and BTC with the cross-sectional tool.
2.5. Methodology
In this section, the different optimization scenarios addressed in the present work are detailed, while further specific information concerning the simulations performed for the assessment of candidate optimum solutions is provided.
In
Table 1, the five optimization cases (CASE A, B, C, D and E) examined in the paper are listed along with the passive control designs applied.
CASE A, concerns material FEC which is based on the thickening of the L/P “leading” and the H/P “trailing” element of the cross-section walls, respectively (see
Figure 1, left). In this optimization study, a different coefficient for the wall thickness increase is applied to each element that varies along the span of the blade. The radial distributions in the two coefficients (i.e., for the L/P and H/P sides) are parameterized using global Bezier interpolation functions with three control points (CP) per distribution. In the optimization process, three parameters per distribution (six in total) are considered as design variables: (i) the value of the coefficient of the first CP, which is fixed at the blade root, (ii) the value of the coefficient of the second CP and (iii) its radial position. The third CP is considered fixed at the blade tip with a fixed value of one (i.e., no thickness change is performed). In addition to the above two coefficients applied to specific segments of the cross-section, a global thickness coefficient is also considered for all the elements in the cross-section. It varies radially, in a discrete manner at the 21 radial sections analysed. In CASE A the planform of the blade is considered fixed and identical to that of the reference blade.
CASE B concerns geometric FEC which is based on the displacement of the spar caps of the two sides of the cross-section in opposite directions (see
Figure 1, right). In this case study, a different percentage shifting (relative to the section chord length) of every cap is considered that varies along the span of the blade. The radial distributions of the percentage displacements (for the L/P and H/P sides) are parameterized using again global Bezier interpolation functions with three CPs. As previously, the same three parameters in each distribution are considered as design variables and the remaining three parameters are kept fixed. Further, a global thickness coefficient is applied to all the elements in the cross-section at 21 radial stations while the planform of the blade is considered fixed.
CASE C concerns the combined application of both the above two methods. It is sought whether a combination of the two methods could be more flexible in tailoring the structural twist angle distribution. In order to reduce the total number of design variables (of the combined application of the two methods) from 12 to 8, the radial positions of the intermediate (mid-span) Bezier CPs of all curves are maintained fixed to the positions that simulations of CASES A and B converged to. The above choice implies that the overall shapes of the radial distribution curves would be quite insensitive to the exact radial placement of the intermediate CP.
CASE D concerns the combined application of both the above two methods (i.e., CASE C) in conjunction with the redesign of the blade twist. As mentioned in the introduction section, FEC entails an indirect, unfavourable BTC effect, induced by the forward sweep deflection. The above BTC effect causes a nose-up twisting of the blade sections that gives rise to increased flapwise loads. Re-twisting of the blade could be a means to mitigate the above effect. The twist distribution of the blade is parameterized using global Bezier interpolation functions with three CPs. The radial position and the twist of the first CP are fixed at r/R = 0.15 and θtwist = −13.1°. The radial position of the second and the third CP are fixed at r/R = 0.25 and 1, respectively, while the twist angle of these CPs is considered as a design variable.
CASE E concerns the combined application of both the above two methods (i.e., CASE D) in conjunction with material BTC (application of offset angle θ to the UD plies over the spar caps, see
Figure 2). As discussed in the introduction section, the aim of BTC in this case is to alleviate any negative effects of FEC on operational loads. BTC begins at r/R = 0.2 and consists of only one design variable (a constant ply offset angle which is applied from r/R = 0.2 up to the blade tip). The mechanism through which BTC alleviates loads is the nose-down twist rotation that takes place whenever the blade is subjected to excessive flapwise deflections, acting as a passive control of the flapwise loads of the blade [
11]. The induced torsional deformation due to BTC suggests that effective optimization of the loading and of the performance requires a redesign of the blade aerodynamic twist. The same parameterization of the blade twist distribution, as well as of the FEC parameters is employed as in CASE D.
As discussed in
Section 2.4, a subset of the most critical DLCs of IEC 61400-1 are simulated within the inner-optimization loop using hGAST, in order to determine the design loads of the components. In terms of normal operation conditions, simulations of DLC-1.3 (extreme turbulence wind conditions in normal operation) at a wind speed of 13 m/s are performed. Earlier studies showed [
10] that the above DLC is the design-driving DLC for the ultimate loads of the reference rotor. It is important to note that rationalization of the computational cost is absolutely critical at this point. To this end, the duration of the normal operation time domain simulation is limited to 50 s. The simulation window is centred on the time instant for which maximum loads of the reference rotor in a simulation with a total duration of 30 min are obtained for most of the cross-sections. It is noted that even for the cross-sections where maximum load occurs at a different time instant, still the load within the selected window is very close to its maximum. The above choice, dictated by the requirement for affordable computer cost, does not ensure that ultimate load is always tracked down. However, previous studies [
10] indicated that following the above approach, loads do not significantly deviate from the actual ultimate loads. Besides, optimized designs are eventually verified on the basis of detailed turbulent wind simulations over a range of wind and operation conditions, as demonstrated in
Section 3.4. As far as idling operation is concerned, the simulations of DLC-6.2 (turbine set in idling mode during the occurrence of 50-year return period wind speed and a grid loss) are at a yaw angle of 30° and at a wind speed of 50 m/s (50-year return period wind speed for the class IA, DTU 10 MW reference turbine). A larger time window of 300 s is considered for the idling turbine simulations, given that the uncertainty of the prediction of the ultimate load in idling mode is significantly higher compared to the one in normal operation mode. This is especially true in the event of flutter instability.
Within the inner-optimization loop, design of the inner structure wall thickness is performed under the constraint that the maximum values for the Tsai–Hill criterion in every cross-section of the reference blade (only due to normal operation extreme loads) are not exceeded. It is noted that the above constraint is rather conservative compared to the requirement that Tsai–Hill criterion values should remain below the failure limit of one. The reason is that, in this way we neglect any possibility that the reference blade may be overdesigned. This could be the explanation of why its Tsai–Hill values are lower than the failure threshold. It is also noted that in the present study, equivalent stress constraints are only applied to the blade loads. Optimization of the tower is out of the scope of the present study. However, in the final assessment of the optimum designs for the three case studies considered, it is showcased that eventually the designed interventions performed on the blades do not harm tower loads.