Study on coupled mode flutter parameters of large wind turbine blades

As the size of wind turbine blades increases, the flexibility of the blades increases. In actual operation, airflow flow can cause aerodynamic elastic instability of the blade structure. Long blades may experience coupled mode flutter due to the bending torsion coupling effect, leading to blade failure. Based on Euler Bernoulli beam theory combined with Theodorsen non directional aerodynamic loads, a blade flutter characteristic equation is established through finite element method. Taking NREL 5 MW wind turbine blades as an example, analyze the influence of parameter changes in different regions of the blades on flutter characteristics. Research has found that paramter changes in the tip region of blade have the greatest impact on flutter characteristics. The vibration frequency shows an overall upward trend with the increase of waving stiffness and torsional stiffness. The flutter velocity of the three regions tends to stabilize as the bending stiffness decreases. The blade flutter speed increases with the increase of torsional stiffness. The radius of gyration is inversely proportional to the flutter frequency and flutter velocity. The impact of centroid offset on blade structure flutter frequency is minimal, but the centroid offset in the tip region has a greater impact on flutter velocity. Increasing the torsional frequency can prevent coupled mode flutter and provide a theoretical basis for blade flutter prevention design.


Structural model
Assuming that the blade main beam is an isotropic beam, based on the Euler beam motion differential equation derived by Hodges 16 , the structural model is shown in Fig. 1.The blade rotates at a constant angular velocity, and the blade root is fixed at the origin O. Establish two sets of coordinate systems: the overall blade rotation coordinate system XYZ, with the x-axis along the length direction and coincident with the elastic axis of the blade when it is not deformed, the y-axis in the rotation plane, and the z-axis perpendicular to the rotation plane.The local coordinate system ξ ηζ is fixed on the blade, ξ is tangent to the elastic axis of deflection, and η and ζ are the main coordinate axes after the deformation of the blade cross-section.O is a point on the elastic axis before deformation, O ′ is a point on the elastic axis after deformation, and u, v, and w are the displacements in the axial, waving, and oscillation directions of the blade.β represents the geometric pitch angle of the blade.
The blade structure model includes two directions of motion: wave bending and twisting.The blade motion equation is established based on the Hamiltonian principle, as shown in the equation: δU represents the variation of kinetic energy, δT represent the variation of strain energy, and δW represents virtual work.If only wave bend and torsion is considered, the differential equation of blade control can be expressed as: (1) where L and M is aerodynamic forces and variable pitch; EI represents the waving bending stiffness, GJ represents the torsional stiffness, and w and φ represent the displacement in the waving and twisting directions, respectively; e is the offset of the centroid; is the speed of the wind turbine; ρ is the leaf density; A is the cross-sectional area of the blade; m is the mass of the blade per unit length; K m is the polar radius of gyration around the elastic axis, and K m 1 and K m 2 are the axis mass radius of gyration around the main neutral axis and perpendicular to the elastic axis, respectively.

Theodorsen aerodynamics
As shown in Fig. 2, the three-dimensional blade rotates at an angular velocity , assuming that the blade rotates in stationary air, the incoming flow velocity is controlled by the blade rotation.When blade flutter occurs,the airflow effectively adheres to the surface of the blade.Assuming that the initial twist angle of the blade is zero, without considering the nonlinear term and the warping effect of the blade, it is assumed that the center of mass position of the rotating beam does not coincide with the elastic axis.AC is the aerodynamic center, E is the elastic axis, and CG is the center of gravity axis.The coupling form of blade bending and torsion degrees of freedom can be described by the two-dimensional airfoil's two degree of freedom coupled vibration model in the figure.
Theodorsen theory 17 uses the potential flow assumption and Kutta condition to determine the aerodynamic forces of oscillating airfoils.The Theodorsen function C(k) is a complex function based on the reduce frequency k, which can be transformed into a quasi-static aerodynamic model by transforming the function into an identity matrix.Lift L and variable pitch M can be expressed as: where C Lα represents the lift coefficient; ρ ∞ represents the air density; reduce frequency k, where k = 2πf f b u and f is the flutter frequency; b is the half chord length; a is the distance from the elastic axis to the midpoint of the chord; U is the incoming flow velocity.

Calculation and model validation of NREL 5MW blade flutter speed
Substitute Eqs. ( 5) and ( 6) into Eqs.( 2) and (3), use finite element method and treat the blade as a variable crosssection linear Euler Bernoulli beam, with three degrees of freedom: w, δω/δx, φ at each node.The variational form of Eqs. ( 2) and (3) requires the interpolation function of the unit to be continuous, while the first derivative of φ and the second derivative of w are non-zero.The structural region is discretized into a 1-dimensional 2-node element, assuming that the φ(x) and w(x) of each node's degrees of freedom are linear functions and cubic polynomials, respectively.The displacement of unit nodes satisfies: P e i φ e i .Using the Galerkin method 18 , the second-order ordinary differential equation of the blade coupled beam element is obtained by weighted integration on the element: Among them, M e s and M e s are the structural mass matrix and aerodynamic mass matrix; C e a represents the aerodynamic damping matrix; K e s and K e a represents the structural stiffness matrix and aerodynamic stiffness matrix, respectively.Equation (8) can solve the complex eigenvalue problem through iterative methods to obtain the flutter frequency and damping of the blade.Calculate the damping ratio by dividing the negative real part of the complex eigenvalue by the modulus.When the damping ratio is zero, it indicates the occurrence of flutter.Define the frequency and velocity of flutter as flutter frequency and flutter velocity.The flutter frequency occurs at the moment of coupling between the first torsional and third bending modes, at which point the aerodynamic damping ratio corresponding to the flutter velocity is zero.
Using the airfoil parameters and section stiffness characteristics of NREL 5 MW blades in reference 19 , a MATLAB finite element flutter analysis program was developed to calculate blade flutter frequency and velocity.To verify the correctness of the program design, the natural frequency of the blade was calculated and compared with existing literature, as shown in Table 1.From the table, it can be seen that the program design is correct.
To verify the correctness and reliability of the aeroelastic model of wind turbine blades, the blade flutter frequency and velocity were calculated and compared with existing literature models, as shown in Table 2. From the table, it can be seen that the calculated results are close to those in reference 7 .

The influence of blade structural parameters on flutter
The aerodynamic elastic instability phenomenon of wind turbine blades belongs to modal coupling.When the first torsional mode of the blade is coupled with a certain waving mode, the blade occurs classical flutter.When the angle of attack and wind speed are the same, the flutter frequency of the blade will change with the variation of structural parameters.Large wind turbine blades have different airfoil cross-sections along the spanwise direction, and each airfoil has different composite material layers on the beam cap, web, leading edge, and trailing edge.Divide the wind turbine blades into three regions according to different airfoils, as shown in Fig. 3.
Assuming that the geometric shape of the blade and the aerodynamic forces acting on the blade remains unchanged, and considering the special characteristics of the airfoil section of the blade, a comparative study is conducted on the effects of structural parameter changes in different regions on the coupled modal flutter characteristics from the aspects of blade bending stiffness, torsional stiffness, center of mass displacement, and radius of gyration.www.nature.com/scientificreports/

The influence of waving stiffness
The classic manifestation of blade flutter is the coupling of bending mode and torsion mode, so the blade waving stiffness has a significant impact on the frequency of the bending dominated mode.Define the non dimensional waving stiffness as: Among them, EI p is the cross-sectional waving stiffness of the blade after changes, and EI o is the cross- sectional waving stiffness of the original blade.Simplify the blade as a cantilever structure, which can be seen as composed of upper and lower beam caps and front and rear web plates, with no bending torsion coupling effect on the front and rear web plates 22 .According to the theory of cantilever beams, Eq. ( 10) holds.
According to Eq. ( 10), as the bending stiffness of the blade decreases and the mass and length of the blade increase, the flutter frequency of the blade will decrease.The variation trends of flutter frequency and flutter speed are shown in Figs. 4 and 5. From Fig. 4, it can be seen that when EI * < 0.4 , the stiffness change in the blade tip region has a significant impact on the flutter frequency, reducing it by 30%.The blade flutter frequency is about 71% of the reference frequency.When 0.4 < EI * < 0.6 , the stiffness change in the blade root area has a significant impact on the flutter frequency, reducing by about 5%.When 0.6 < EI * < 1 , the change rule of the tip and middle regions of the blade is similar, with a maximum reduction rate of about 6% at the tip.When EI * > 1 , the trend of changes in the three regions tends to be consistent, with the tip and middle regions of the blade having a significant impact on the flutter frequency.As shown in Fig. 5, the frequency of the wave mode decreases with the decrease of bending stiffness, and the coupling time with the torsional mode is delayed.When EI * < 1 , the flutter speed increases relative to the reference speed.When EI * < 0.4 , the influence of the tip region is significant, and the trend of changes in the middle and root regions of the leaves is consistent.When 0.4 < EI * < 0.8 , the changes in the tip and middle regions of the blade have a relatively small impact on the flutter speed, while the changes in the root region have a significant impact on the flutter speed.When 0.8 < EI * < 1 , the blade flutter speed decreases to the reference speed as the stiffness increases, and the stiffness continues to increase.The flutter speed slowly decreases and tends to stabilize.

The influence of torsional stiffness
The dimensionless torsional stiffness is defined as: From equation ( 11), it can be concluded that: The torsional frequency of the blade decreases with the decrease of torsional stiffness, and is no longer coupled with the torsional and third-order wave modes, but rather coupled with the torsional and second-order or first-order wave modes to produce flutter.The variation trends of flutter frequency and flutter speed are shown in Figs. 6 and 7.
(11)  www.nature.com/scientificreports/From Fig. 6, it can be seen that the blade flutter frequency decreases with the decrease of torsional stiffness, and the influence of the tip and mid blade regions on the flutter frequency is similar.From Fig. 7, it can be seen that the flutter speed in the tip and middle regions of the blade decreases rapidly with the decrease of torsional stiffness, while the decrease in the root region is slow.The torsional stiffness of the blade tip region has a significant impact on flutter velocity.

The influence of centroid shift
The distance e * between the center of gravity and the elastic axis in the blade section is a dynamic instability parameter that affects the overall structure of the blade.The centroid position of the airfoil section of the blade is defined as the dimensionless centroid offset in Eq. ( 13): Among them, e p and e o are the displacement of the center of mass of the changed blade and the reference blade, respectively.Flutter only occurs when the center of mass moves towards the trailing edge of the blade.The variation trends of flutter frequency and flutter speed are shown in Figs. 8 and 9.
From Fig. 8, it can be seen that the flutter frequency slightly increases with the increase of Centroid offset, and the influence of the blade tip egion is slightly greater.The maximum flutter frequency change rate is 2%.From the flutter speed curve in Fig. 9, it can be seen that the influence of the blade tip regin is bigger.As e * decreases, the center of mass moves backward towards the trailing edge of the blade, resulting in a decrease in flutter speed and a greater likelihood of flutter occurring.

Conclusion
The blade is simplified as an Euler Bernoulli beam, and combined with the unsteady aerodynamic theory of Theodorson, the coupled flutter dynamic equation of the blade considering the Characteristics of coupled motion between waving and twisting.The blade flutter characteristic equation is solved using the complex eigenvalue method.The frequency and velocity of NREL 5 MW blade flutter was solved, and the correctness and reliability of the wind turbine blade aeroelastic model were verified by comparing it with existing literature models.
A study was conducted on the coupled modal flutter parameters of NREL 5 MW blades, and it was verified that the coupling of the first torsional mode and a certain wave bending mode of the blades can cause classical flutter phenomena.The research results indicate that: By analyzing the structural parameter changes in the blade tip, middle, and root regions, it was revealed that the parameter changes in the blade tip region have the greatest impact on flutter frequency and velocity.The flutter frequency shows an overall upward trend with the increase of waving stiffness and torsional stiffness.The flutter velocity of the three regions tends to stabilize as the bending stiffness decreases.The blade flutter speed increases with the increase of torsional stiffness.The radius of gyration is inversely proportional to the flutter frequency and flutter velocity.The impact of centroid offset on blade structure flutter frequency is minimal, but the centroid offset in the tip region has a greater impact on flutter velocity.By analyzing the parameter changes in different regions, guide the aerodynamic elastic design of blades and provide theoretical basis for preventing blade flutter.In the problem of blade flutter in megawatt wind turbines, wave frequency frequency and torsional frequency has become the main considerations for future design.In the design of blade anti flutter, the influence of torsional stiffness should be emphasized.How to maintain a high torsional frequency and avoid coupling with a certain bending mode has become a key factor in improving the flutter limit.

Figure 2 .
Figure 2. Aerodynamics on blade elements under bending torsion coupling mode.

Figure 3 .
Figure 3. Three Areas of the Blade.

Figure 8 .
Figure 8. Flutter frequency variation with center of mass.

Figure 9 .
Figure 9. Flutter speed changes with center of mass.

Table 2 .
Comparison of flutter speed and frequency with existing research.