# The Impact of Critical Flutter Velocity in Composite Wind Turbine Blade with Prebend Condition

Authors: H. S. Sunil Kumar; K. R. Jagadeesh; R. B. Anand; T. Rangaswamy; Srikanth Salyan; Jenoris Muthiya Solomon; Joshuva Arockia Dhanraj; Joshua Stephen Chellakumar Isaac Joshua Ramesh Lalvani
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

The present work focuses on the effect of flutter in prebend 100 m horizontal axis wind turbine blade (HAWT) within the stability limits. The study was carried out with an advanced beam model for idyllic structure in a DU-97-W-300 cross-sectional area. A Galerkin type of approach has been applied to derive the equations, and the analysis was performed using a standard FEA code which involves the PK method and double lattice method for calculating flutter solution and aerodynamic loads respectively. The results reveal the significance of inducing prebending to improve the stability of the blade structures, and hence, the flutter velocity has moved from 11 m/s to 23 m/s. Furthermore, the output highlights the effect of prebending on the structural stability and also the flutter limit was found to be lengthened.


---

## Body

## 1. Introduction

The tremendous increase in energy demand and the factors affecting environmental degradation due to production of energy has resulted in dependence on clean energy resources. Among various resources wind energy has outperformed its counterparts for reliable source for renewable energy resources and vast amount of research and implementation is on rise in the recent decades. The design optimization of wind turbine and blades is essential to extract more energy to meet increasing energy demand which results in study of increasing in blade size while exposing them to high load conditions; at the same time wind blades should be flexible and slender. Because of such slenderness issues there are multiple phenomena which tend to make the blades unstable. Figure1 and Figure 2 illustrate the wind turbine blade with prebend configuration and prebend blade terminology.Figure 1
Wind turbine blade with prebend configuration [1].Figure 2
Prebend blade terminology.The preliminary sciences of such instabilities like aeroelasticity, flutter instability, and relative factors have been clearly explained in [2]. The structural stability and its mathematical relation with the boundary and load conditions have been derived in [3] while, various types of dynamic instabilities and relations are detailed in [4]. The flutter instability has been found to be capable of becoming a prime design driver for future wind turbine blades and thus has emerged as an area of research, where researchers and scientists have come up with several methods to detect, control, and eliminate such instabilities. The flutter instability is considered more important as it is much difficult to detect when compared to resonance. One such method for enhancing the stability is what is called prebending of the blade structure. To tackle such issues without sacrificing the expectations of energy output, prebending of the blade structures is performed to reduce blade weight, improve overall energy generation capacity per annum, and increase tip deflection to a certain extent. The aim of this work is to investigate the wind turbine with prebend condition and its flutter velocity. The preliminary studies investigate the aeroelastic behaviour of blade structures revealing that the flutter limit of a blade depends upon natural frequencies of torsional as well as flapwise direction of a blade [5]. Further, studies on the influence of system parameters and its uncertainties for the wind turbine blades on onset of flutter suggest that structural randomness can make the blade more unstable by decreasing the critical speed for fluttering below the operational speed designated for the wind turbine [6]. After detecting the occurrence of flutter instability, researchers moved further to find optimum methods which show the closeness in approximation using reliability evaluation methods [7]. As further studies kept revealing more and more flutter control parameters, the analysis techniques gained more correctness and reliability [8]. The aeroelastic stability in a wind turbine is presented by another simplified method using Eigen value approach [9] and linear aerostatic turbine model [10] to find the critical frequency of three-bladed rotor arrangements. Similarly, [11, 12] study the fluttering of two-dimensional as well as three-dimensional nonlinear cyclic oscillations and its effect on aeroelastic instabilities. An approach to flutter analysis study by nonlinear damping for aerofoil using equivalent linearization [10] and Range Kutta method [13] for various tip speed ratio is reported in literature.The above analysis investigated the structural instability and flutter analysis of large wind turbine blades; there aroused a necessity to induce structural modification techniques and study their effect on the blade reliability and longer life under the desired operating conditions. Among such modifications, prebending is one where the blade is bent towards aerodynamic pressure to ensure interference margins with the tower. Although this technique has provoked more difficulties in transportation and storage, the studies show that this method has overall impact in the performance and reliability of the structure. Prebending being an effective technique has also increased the complexity of design and evaluation as the induced bend has an impact on length of the blade and on angle of attack, thereby affecting the aerodynamic load. Generally, a coupled aerodynamic-structural iterative analysis will predict the aerodynamic loading and structural response in a better way. Paper [14] has conducted a parametric study on a swept blade for a power output of 750 kW wind turbine along with modification in design conditions. Papers [15, 16] investigated the flutter speed for wind turbines with two different types of straight and swept blades. Paper [17] has established several methods for predicting the prebend shape of the wind turbine blades such as when the blades are straightened into their design configuration when they are subjected to both wind and inertial loads. Paper [1] studied the comparison of multiple prebend positions and performed an optimization analysis based on multiobjective particle swarm optimization algorithm (MOPSO) to study the optimized design of blade in order to perform well in reliability, cost, energy efficiency, and stability. The frequency along with deflection change in the geometry with respect to backward sweep is presented in [18]. References [19, 20] made investigations and concluded that backward swept shapes also help in load alleviation without greatly increasing blade root torsion and life-time equivalent fatigue moment. Mehrvarz [21] addressed a Timoshenko beam’s asymptotic stability with respect to the effect of piezoelectric actuators. Sahoo [22] compared numerical results computed using MATLAB with experimental data for higher order model of a composite layered shell structure. Meanwhile, Nahavandian [23] investigated the linear stability of thixotropic fluid obeying Moore model using a temporal stability analysis with infinitesimally small perturbations. Fakhraei [24] studied chaotic behaviour of a ground oscillating system by modeling a full nonlinear seven-degree freedom system with an addition of freedom for passengers in one degree. Lobitz [25] studied 3.5 MW turbine blades and concluded that the natural frequencies are subjected to change because of torsional and bending stiffness changes in blades. Meng [26] investigated the flutter frequency as well as critical speed for NREL 5 MW turbine blade for coupling of first torsional and second flapwise modes (instead of the third). Owens et al. [27] developed BLAST, an aerostatic stability design tool, which is an extension of the NASTRAN for flutter analysis for the study with better flutter trends, and studied the impact of composite blade to study the start of instability in blade in two different blades (Wind PACT of 1.5 MW and SNL 100–00 blade). The flutter onset was detected at 13.05 rpm and 40.6 rpm for the SNL 100–00 blades and Wind PACT, respectively. Jonkman et al. [28] developed a design tool called HAWCStab stability tool to study the properties of the blade by using eigenvalue analysis and linearization of nonlinear aeroelastic methods. Shakya et al. [29] observed that flap-torsional stiffness had elevated impact on critical flutter speed which is improved by 40% due to the unbalance in the complete segment of the blade with symmetric skin. Tamrakar et al. [30] found that the crack had a considerable effect on increasing the natural frequency of the rotor system in contrast to the uncracked rotor system. This is due to the buildup of the strain energy in the surrounding area of the crack. Yu Hua et al. [31] discussed the analysis of joint state and unknown input evaluation for linear systems with an unknown input and proposed two novel filters as advanced algorithmic control. Further Kai Wang et al. [32] proposed the Kalman filter algorithm for the estimation of state of charge (SOC) in supercapacitors and lithium batteries of electric vehicles. Chunli Liu et al. [33] suggested the new method of estimating the remaining useful life prediction of supercapacitors in energy storage devices. Further the occurrence of flutter instability in delaminated wind turbine blades has been studied [34]. Although these studies have analysed and revealed how prebending of the structure is altering the structural properties and stability, very few studies have been made to understand its influence on flutter limit. But, from these studies, it is clear that prebending seems to have a good impact on the blade’s flutter performance. The change of dependency of modes on load conditions has not been studied. This study highlights the change in dependency of a mode on aerodynamic load due to prebending induced on the structure.


## 2. Mathematical Modeling

The equations of motion for a force applied on a body by airflow are dependent on the body’s elastic deformation and it obeys both elasticity and aerodynamics theories. Figures3 and 4 show the flow direction in a blade arrangement and representation of the blade aerofoil DU-97-W-300, respectively. The fluid flows demonstrate the relationship between elastic and aerodynamic parameters. The configuration is supposed to be such that the blade's centre axis is perpendicular to the flow direction and U represents the available free stream air velocity.Figure 3
A view of blade arrangement and flow direction.Figure 4
2D diagram of aerofoil DU-97-W-300.A blade can be considered a cantilever beam structure because of its high aspect ratio structure. The current study focuses on the effects of flutter instability, as well as the deformations that result in both torsion and flap directions. The deformation is measured by a deflection ’X’ and rotation ’θ,’ where X indicates downward displacement and θ indicates a twist with the leading edge up.(1)∂2∂x2EI∂2X∂x2+m∂2X∂t2+myθ∂2θ∂t2+FL=0,−∂∂xGJ∂θ∂x+Iθ∂2θ∂t2+myθ∂2w∂t2+FM=0.Bending and torsional rigidities are represented by EI and GJ, respectively, whereas lift force and twisting moment are represented byFL and FM, which are as follows:(2)FL=ρU22cC1,FM=ρU22c2C2.The lift coefficient C1 and the moment coefficient C2 are calculated using quasi-steady subsonic aerodynamic theory.(3)C1=dC1dxθ+1Udxdt+1U34C−x0dθdt,(4)C2=−cπ8Udθdt−14C1.The lift force is one-fourth (1/4) of the chord value, as shown in equation (4).The moment and aerodynamic force of the prebend blade are given by quasi-steady double lattice approach in the aforementioned equations. By using (3) and (4), we get(5)∂2∂x2EI∂2X∂x2+m∂2X∂t2+myθ∂2θ∂t2+ρU22cdC1dθθ+1UdXdt+1U34c−x0dθdt=0,(6)−∂∂xGJ∂θ∂x+Iθ∂2θ∂t2+myθ∂2w∂t2+ρU22c2−cπ8Udθdt−14dC1dθθ+1UdXdt+1U34c−x0dθdt=0.The boundary conditions to be used in the equation for the displacements ‘θ’, ‘X’ are given below:(7)X=∂X∂x=0;forxvalues equal to0,(8)∂2X∂x2=∂3X∂x3=∂θ∂x=0;forxvalues equal toL.In the corresponding (5) and (6) the X, Y are considered as two independent equations for the values of U and yθ  = 0, where U and yθ represent the aerodynamic and inertial couplings and determined using the following equation: (9)Xx,t=Axeλtandθx,t=Bxeλt.The ODE equation is obtained by dividing the equation byeλt. By applying equation (9) in (5) and (6), we get a complex DE.(10)EIA'′'′+ρU22cdC1dθB+λρU2cdC1dθA+c34−yθcB+λ2mA+yθB=0,(11)−GJB−ρU22c2yθc−14dC1dθB−λρU2c2yθc−14dC1dθA+cyθc−1434−yθcdC1dθ−π8B+λ2myθA+IθB=0,where ( )’ =  d/dx and ( )’’ =  d2/dx2The boundary conditions are applied and remain the same as before, except the valuesX and θ, which are replaced by A and B. The partial derivatives of X are replaced by total derivatives. Since there is no way to solve (9) and (10) in closed form, an approximate answer is used in the study. The impact of airflow velocity U on wind blade is examined by eigenvalue which determines the structure’s stability condition. The parameter (λ) depends upon the air velocity U in a continuous form. For the nonzero values of U, a small value takes an exponent value and becomes a complex value λ = α + iω. This value is examined using non-self-nearby (adjoint) values. For values which are not equal to zero but are very small, both the U and (dC1/dθ) with a condition of less than< 2π, the losses of blade are more with most of the energy being lost to the surrounding air. This causes the oscillatory motion to be damped and instability exhibiting asymptotically for which ‘α’ takes negative value.As the value of U rises, “α” changes from negative to positive, and the “α” sign changes and becomes unstable at a point where α becomes positive. Also, when the critical velocity α = 0 is the air velocity, that corresponds to 0 and is represented by Uc. There critical velocity values for U have numerous different values; however the lowest value for critical velocity value is used to determine the limit or operational range. The value of U can be used to distinguish between two critical circumstances. The critical divergent state is given by the values for α = 0 and ω = 0 for the wing and critical flutter condition for blade is given for the conditions of α = 0 and ω ≠ 0.


## 3. Methodology

The Galerkin approach is used to solve the coupled partial differential equations of motion, which discretizes the equations spatially and provides an analysis of flexible structural instabilities owing to fluid flow.The torsional and flapwise displacements using Galerkin approach are given by(12)A=∑j=1nφ1x.q1t,B=∑j=n+1n′φ2x.q2t.In this equation,φ1 and φ2 represent the torsion and flap wise mode whereas n and n1 represent the mode numbers in both torsional and flapwise directions. Similarly, q1, q2 are synchronized with respect to the directions. The nondimensional lift and moment equations are substituted in the equation and multiplied by the values of Eigen values. The values are integrated with respect to the length of the blade. The obtained equation corresponds to the eigenvalue problem, solving which gives the stability analysis of the blade structure. Figure 5 explains the computation of unsteady aerodynamic coefficient using double lattice method and solving the eigenvalue problem to calculate the flutter velocity using PK method.Figure 5
Flow process of the analysis.The eigenvalue problem equation upon integration for the values of 0 toL, for the length of the blade,(13)K+U2Q+λ2Ma=0.The equation is solved for the increased values ofU to solve the value of Uc. The λ value will be negative for the small values of U and it is approximately equal to 0 for the first value of the critical air speed Uc. On approximating A and B, Uc values are obtained for the values of n and m equal to 1. For the values of λ = iω, Uc will be equal to U, which gives an equation in the form of determinant.(14)iω0−100iω0−1k11U2h12iωm11+Ul11iωm12+Ul120k22+U2h22iωm12+Ul21iωm22+Ul22=0.By substituting into the real part of (13), we can write the quadratic equation in U2 as(15)CU4+DU2+F=0.On substituting the values in (13),(16)U2=−D2C±12CD2−4CF.The iterations are carried out for a number ofU values, with the air speed with the least positive value taken as critical air speed for which the structural stability exists for the values below this critical air speed. The critical air speed is represented by UC and is compared to that with prebend and flat blade to study the impact of flutter limit analysis in a wind turbine blade.


### 3.1. Estimation of Onset for the Flutter Instability

The predication of flutter limits in wind turbine blades is carried out by time dependent aeroelastic analysis or eigenvalue analysis. The eigenvalue analysis involves a model of steady state equilibrium, which is a linear model for blade. This eigenvalue problem is used to study the instability mechanisms. Standard FEA blade model code is used for different blade arrangements and its properties are given in Table1. The 100-meter blade is used in this study with the boundary conditions applied as stated in (6) and (7), where the blade behaves similar to the cantilever beam. A 5 m tip deflection and an initial bend are also included in the structure. At an altitude of 121.92 meters, the airspeed is gradually increased from 0 m/s to 50 m/s, where the air density corresponds to replicating the blade's surrounding conditions. The blade DOF is depicted in 2D diagram (Figure 4). The mathematical model is solved by PK method for flutter analysis and validates equation by developing aerodynamic matrices for variables such as air speed and frequency. The blade model imported into the standard FEA code is assigned with the composite layup arrangement whose characteristics are listed in Table 1.Table 1
Ply stiffness properties.
MaterialExy (GPa)Eyy (GPa)Gxy (GPa)υxyDensity (kg/m3)Gel coat3.443.441.380.31235Resin3.53.51.40.31100Foam0.2560.2560.0220.3200TRIAX27.713.657.20.391850E-LT-550041.814.02.630.281920Saertax13.613.311.80.511780


### 3.2. Parametric Study

Now, the parameters like aerodynamic load are varied by varying the wind speed and from its output, the damping curves for the prebend structure are presented with the help of Campbell diagram. Two graphs are obtained for damping ratio due to dynamic as well as inertial forces acting due to the speed of the rotor and change in wind velocity. Also, the comparison shows the influence of aerodynamic load and inertial forces on the modes can be understood. Similarly, the change in stability of the blade due to prebending can also be observed. And also, the flutter speed comparisons with previous research and present one have been shown in Table2 and further mode frequency comparisons have been tabulated in Table 3.Table 2
Flutter speed comparison.
ReferencesBladeBlade length (m)Critical flutter speed (rpm)Lobitz [25]WindPACT 1.5 MW3342.3Hansen [9]NREL 5 MW6324.0Owens and Griffith [27]SNL 10010013.05Current studySNL 10010014.12Table 3
Mode frequency comparison.
FrequencyMeng [26]Jonkman et al. [28]Paryaya et al. [5]Hansen et al. [9]Present study1st flapwise0.720.70.640.70.62nd flapwise2.052.021.861.81.753rd flapwise4.37—4.343.63.491st torsional5.62—5.3986.04


## 3.1. Estimation of Onset for the Flutter Instability

The predication of flutter limits in wind turbine blades is carried out by time dependent aeroelastic analysis or eigenvalue analysis. The eigenvalue analysis involves a model of steady state equilibrium, which is a linear model for blade. This eigenvalue problem is used to study the instability mechanisms. Standard FEA blade model code is used for different blade arrangements and its properties are given in Table1. The 100-meter blade is used in this study with the boundary conditions applied as stated in (6) and (7), where the blade behaves similar to the cantilever beam. A 5 m tip deflection and an initial bend are also included in the structure. At an altitude of 121.92 meters, the airspeed is gradually increased from 0 m/s to 50 m/s, where the air density corresponds to replicating the blade's surrounding conditions. The blade DOF is depicted in 2D diagram (Figure 4). The mathematical model is solved by PK method for flutter analysis and validates equation by developing aerodynamic matrices for variables such as air speed and frequency. The blade model imported into the standard FEA code is assigned with the composite layup arrangement whose characteristics are listed in Table 1.Table 1
Ply stiffness properties.
MaterialExy (GPa)Eyy (GPa)Gxy (GPa)υxyDensity (kg/m3)Gel coat3.443.441.380.31235Resin3.53.51.40.31100Foam0.2560.2560.0220.3200TRIAX27.713.657.20.391850E-LT-550041.814.02.630.281920Saertax13.613.311.80.511780


## 3.2. Parametric Study

Now, the parameters like aerodynamic load are varied by varying the wind speed and from its output, the damping curves for the prebend structure are presented with the help of Campbell diagram. Two graphs are obtained for damping ratio due to dynamic as well as inertial forces acting due to the speed of the rotor and change in wind velocity. Also, the comparison shows the influence of aerodynamic load and inertial forces on the modes can be understood. Similarly, the change in stability of the blade due to prebending can also be observed. And also, the flutter speed comparisons with previous research and present one have been shown in Table2 and further mode frequency comparisons have been tabulated in Table 3.Table 2
Flutter speed comparison.
ReferencesBladeBlade length (m)Critical flutter speed (rpm)Lobitz [25]WindPACT 1.5 MW3342.3Hansen [9]NREL 5 MW6324.0Owens and Griffith [27]SNL 10010013.05Current studySNL 10010014.12Table 3
Mode frequency comparison.
FrequencyMeng [26]Jonkman et al. [28]Paryaya et al. [5]Hansen et al. [9]Present study1st flapwise0.720.70.640.70.62nd flapwise2.052.021.861.81.753rd flapwise4.37—4.343.63.491st torsional5.62—5.3986.04


## 4. Results and Discussion

The modes of the prebending are studied for natural frequencies and analysed with the same flat blade parameters. The frequency of all the modes in the study is in proximity to each other and thereby the prediction behaviour can be inferred due to the participation of all modes in flutter analysis of the blade structure. Table4 compares the modal analysis output obtained for flat and prebend configuration of the blade structure.Table 4
Natural frequency of prebend and flat modes.
ModePrebend blade mode frequency (Hz)Flat blade mode frequency (Hz)1st bending mode0.280.261st torsion mode0.520.452nd bending mode0.990.903rd bending mode2.211.944th bending mode2.242.032nd torsion mode4.073.681st chordwise bending mode4.914.663rd torsion mode5.615.244th torsion mode6.465.872nd chordwise bending mode9.238.38The results are plot showing the Eigen frequencies with respect to the increasing rotor velocity of blade; in Figure6, it is observed that the rotor velocity has a significant influence among 1st, 2nd, and 3rd bending mode and 1st torsional mode. Thus, this Campbell diagram can act as evidence to predict the modes that would contribute in the flutter analysis of the blade.Figure 6
Eigen frequencies for the blade with respect to rotor velocity.From Figure7, the damping curve for bending modes is nonlinear and indicates the stability of the bending modes when the wind velocity is increased. The 3rd and 4th bending modes show a higher sensitivity to the wind velocity changes than the 1st and 2nd bending mode.Figure 7
Damping curves of the blade bending modes with respect to the wind velocity.Figure8 shows the damping curves with respect to wind velocity. It is also noted that 1st and 2nd torsion mode are capable of easily initiating instability as their curves lay close to zero throughout the graph. But the 1st torsion mode seems to be more stable while the 2nd mode is not as stable as the first one. The graph also indicates the dominance of the 2nd torsion mode when it comes to contributing to the flutter mode as change in dampening ratio shows the sign change in the flutter velocity.Figure 8
Damping curves in blade torsion modes with respect to the wind velocity.The graph in Figure9 with curves plotted for damping with respect to rotor velocity clearly depicts the drastic change in slope of all bending modes, where 1st bending mode shows exception for a critical rotor speed in the range between 20 and 30 rpm. The 3rd and 4th bending mode show a mild and sensitive variation with respect to the rotor velocity. The more stable conditions are obtained for low rotor velocities and it can be concluded that, with changing the design to improve the flutter limit, more ideal rotor speed is the velocity of less than 20 rpm for aerodynamic properties with an onset of flutter stability.Figure 9
Damping curves of blade bending modes with respect to rotor velocity.Figure10 shows the damping curves variation along the rotor velocity. The graph shows that the 1st torsion mode is quite stable compared to the other torsion modes. Meanwhile, torsion modes 2 and 3 are found to be highly sensitive and nonlinear to rotor velocity.Figure 10
Damping curves of blade torsion modes with respect to rotor velocity.Figures11 and 12 reveal real part of eigenvalue for the bending and torsion modes and the plot is in acceptance with the inference gained from the damping ratio plots.Figure 11
Impact of wind velocity for real eigenvalue in bending modes.Figure 12
Impact of wind velocity in torsion modes.Figures13 and 14 of flutter analysis show recorded behaviour of flap and torsional modes with the onset of flutter in the prebend blade.Figure 13
Imaginary part of the eigenvalue with respect to wind velocity.Figure 14
Prebend blade with initial tip deflection.In Figures15 and 16, the initial tip deflection and simultaneous bending and torsion contribution of bending and torsion mode are clearly observed. It can be inferred that the deformation of blades is both flapwise and torsional revealing the excitation in both modes under the same frequency, and structure dampening capacity is exceeded due to the higher absorbed energy.Figure 15
Bending mode of the prebend blade.Figure 16
Torsion mode of the prebend blade.


## 5. Conclusion

The investigation was carried out with a focus to detect the occurrence of flutter instability in a horizontal axis large wind turbine blade with prebend (HAWT). The literature presents and evaluates the mathematic model for a prebend blade for wind turbines at natural frequencies. The study was carried out by advanced beam model for ideal structure in a DU-97-W-300 cross sectional area. A Galerkin type of approach has been used to derive the equations and the analysis was performed involving a standard FEA code which involves PK method and double lattice method for calculating flutter solution and aerodynamic loads, respectively. The study is applied to structural stability, estimation of flutter region, nonlinear influence of blade length, and rotor velocity. The main conclusions from this paper are summarised as follows:(1)
The results reveal the significance of inducing prebending to improve the stability of the blade structures and hence, the flutter velocity has moved from 11 m/s to 23 m/s(2)
The results also ensure that prebending in blades can be used significantly to modify its performances and structural stability without changing its aerodynamic properties(3)
New technique of obtaining damping curves with and without aerodynamic loads helps the designer to decide whether to modify the aerodynamic properties or structural properties for improved flutter limitsFor future work, the methods and theories used here can estimate the flutter region but cannot provide any details about interaction of fluid and the dynamic behaviour of the structure beyond the flutter point to study requires a more comprehensive model. The recommendation for further studies is that a wind turbine blade can be manufactured in sections (split blades) and it can be easy for transportation as well as for installation. Thus, integration of joints and structural nodes should be evaluated for performance parameters. It can be extended to determine the impact of prebending of the blade on AEP (annual energy production) which can also be studied as the prebending will cause the turbine to have a higher AEP than rotors having straight blade due to the increase in the angle of attack along the length of the blade. And also, the impact of glue delamination of blade structure on flutter limits can be analysed [34].



---

*Source: 2050821-2022-03-02.xml*