Comparison of beam theories for characterisation of a NREL wind turbine blade flap-wise vibration

Abstract

The harsh weather conditions offshore wind turbine blades are working under have necessitated the need for a more effective approach for detecting their failures. Dynamic characteristics e.g., vibration, carry useful information related to blades structural health. Critical analysis of the blade’s free vibration represents one of the fundamental steps for assessing structural dynamics. When a rotating blade deflects, either in the plane of rotation or perpendicular to it, the centrifugal force exerts inertia force along the blade span, which increases the blade’s natural frequencies compared with the stationary ones. However, the influence of different blade parameters on the flap-wise vibrations is not very well understood. In this paper, the effects of such parameters on dynamic characteristics of National Renewable Energy Laboratory (NREL) 5-MW wind turbine blades are investigated using different beam theories. The examined models have been used to determine the natural frequencies and mode shapes of the National Renewable Energy Laboratory (NREL) 5-MW wind turbine. Results demonstrate that increasing angular velocity significantly impacts the natural frequencies and mode shapes. The rotary inertia is found to influence the free vibration responses of the studied blades. Moreover, increasing hub radius, pre-cone and pitch angles are found to have less influence on the natural frequencies. Mode shape comparisons were also carried out using MSF (modal scale factor) and MAC (modal assurance criteria). Compared to the other investigated methods, Bernoulli’s based algorithms are found to produce less accurate results.

Introduction

The blade is considered one of the most critical components in a wind turbine due to the unpredictable external forces such as the atmospheric environment and varying operational conditions. The challenge is further complicated as the size of the blade becomes larger and as the wind turbines are constructed in more difficult terrains.^1,2 Blades are a major contributor to the downtime of wind turbines and account for over 41% of total failures. ³ Therefore, vibration analysis has been and continues to be a subject of immense importance because it leads to understanding the dynamic characteristics of wind turbine blades. Natural frequencies and mode shapes play an essential role in the design and control of these structures. The blade is usually modelled as a rotating beam. The out-of-plane bending vibration of the Euler-Bernoulli beam is examined in, ⁴ which uses the variational iteration technique. Natural frequencies and mode shapes of a rotating beam are studied for various angular velocities and different taper ratios. The flexural free vibration analysis of a double tapered Euler-Bernoulli is investigated in. ⁵ The differential transform method technique is used to solve the governing differential equation. The centrifugal force is incorporated in order to find its effect on the dynamic characteristics. A composite Timoshenko beam (TB) with single delamination travelled by a constant amplitude moving force is modelled in. ⁶ The governing differential equation of motion is derived, and the effects of shear deformation and rotary inertia have been accounted for. Eigen-values technique and Ritz method were used to obtain the dynamic response under the action of moving force. The mechanical behaviour of the delamination layers is modelled by a piecewise-linear spring foundation. ⁶ The flap-wise bending vibration of a rotating tapered Rayleigh beam is analysed in. ⁷ The integral equation method was developed to determine the natural frequencies. The influences of hub radius, rotational speed, taper ratio and rotary inertia on the natural frequencies are elucidated. The horizontal axis wind turbine blades were developed with varying cross-section that gradually twists from circular at the root to different aerofoil profile along the blade span. These types of blades are usually installed on the hub with a pre-cone angle and a pitch angle during the operation process. The contributions of pre-cone angle and the pre-twist angle variation on the kinetic and strain energy, the influences of centrifugal force and rotary inertia are addressed in. ⁸ Hamilton’s principle in combination with differential transform method is used, ⁹ the free vibration analysis of a Timoshenko beam featuring in the direction of the transverse bending-lateral bending-torsion coupling is analysed based on energy expressions. Thus, the kinetic and the potential energy are used to incorporate the rotary inertia, shear deformation, rotational speed, hub radius and bending-bending-torsion coupling deformation. The work in ¹⁰ reports the development of a vibration monitoring system to estimate the deflection at the tip of the blade, especially in cases when the blades begin to vibrate excessively. Analytical methods are used in ¹¹ to investigate the forced vibration of the symmetric cracked Timoshenko beam under harmonic excitations. The response displacement and the rotating angles of a coupled bending and torsional vibrations of the cantilever Timoshenko beam are studied. Analytical methods are used in ¹¹ to investigate the forced vibration of the symmetric cracked Timoshenko beam under harmonic excitations. The response displacement and the rotating angles of a coupled bending and torsional vibrations of a cantilever Timoshenko beam are studied. The TB is divided into several segments due to the existence of cracks. ¹¹ The general form of Green’s functions of each beam segment with the unknown boundary conditions is derived. The study described the matching conditions of the cracked cross-section and the relationship between the boundaries of the beam segments. The damaged cross-section is modelled as a mixture line-spring, and the compliances are determined through the Paris equation. The study found that the depth, location and number of cracks have significant effects on the natural frequencies of the TB. Also, the location of cracks can be easily positioned by the discontinuity of the beam response. If the wind turbine blade has a long span, it is appropriate to simplify the blade as a non-uniform beam. In each blade section, the shear centre may not be coincident with the beam centroid, which leads to a coupled bending and torsional vibration. There have been several studies examining the influences of coupled vibrations on the dynamic response. The differential transform method is applied in ¹¹ to solve the coupled partial differential equations. ¹² A non-linear beam model based on the geometrically exact beam theory is used in ¹³ to build the wind turbine blade. Taking the NREL 5MW wind turbine blade as an example, the free vibration analysis is used to investigate the dynamic analysis under periodic unsteady inflows.

As stated, due to the externally applied force from a precise atmospheric environment and varying operational conditions, the free vibration analysis of the wind turbine blade is the first issue that should be addressed accurately to find the modal parameters. Numerous studies with different approaches have been carried out to investigate the flap-wise vibration of rotating beams. Some studies have adopted Euler-Bernoulli rotating-beam theory^4,5 or even non-rotating Bernoulli beam due to its simplicity. Timoshenko beam theory which added the rotary inertia and shear deformation was also considered for vibration analysis.^6,7 Rayleigh beam theory^14,15 is among the main methods used as it investigates free vibration by considering the rotary inertia, which is found to have an influence on the dynamic characteristics and forced response of wind turbine blades. The work in ⁸ uses Rayleigh beam theory and incorporates the influence of the pitch and pre-cone angles on the dynamic characteristics of the horizontal axis wind turbine blade in the flap wise direction. The study in ¹³ is adopted the geometrically exact beam theory to find the natural frequencies and corresponding mode shapes in the flap-wise and edge-wise directions to find the deflection response in both directions. Despite the significant contribution being made to analyse the flap-wise vibration, there are some shortcomings. The aforementioned studies^4,5,6,7 did not achieve accurate free vibration results that would predict the dynamic response of the blade. However, unless a comprehensive dynamic model takes into account parameters such as centrifugal forces, rotary inertia, gravity, angular velocity, hub radius, pre-cone and pitch angles, the supplying input data to model-based control would be inaccurate. A significant advance was made in ⁸ by scrutinising the effects of harmonic excitation due to the gravitational force, pre-cone and pitch angles, in addition to the other parameters. However, a shortcoming in the study was not applying the comparison of methods on a blade with the same geometry. Furthermore, despite mentioning the significant effects of pitch and pre-cone angles with increasing the rotating speed, the results did not clarify these effects or the mentioned factors' effect on the high natural frequencies. The work in ¹³ is found the natural frequencies of a single blade at zero rotation speed. However, the natural frequencies without considering the influences of centrifugal force and different values of rotation speed would not be accurate and sufficient as a comprehensive model.

In this study, the comparison between different beam theories is established for the model of modern wind turbine blades. The finite element code is written using the specification of 5MW NREL blade for structural modelling and simulation studies. Modal properties, including the natural frequencies, translational mode shapes, angular mode shapes, modal scale factor and modal shape correlation coefficient, are compared. The effect of changing the blade parameters, including hub radius, pitch and pre-cone angles, on the blade performance is investigated. The effect of increasing the angular velocity on the modal parameters is also studied.

Derivation of the governing equations

In this study, the horizontal axis wind turbine blade under the flap-wise vibration is considered. The effect of rotary inertia, centrifugal force, and the effect of pre-cone and pitch angles are included. The blade deformation can be expressed by the superposition of static displacement-time independent- and dynamic displacement due to the influence of angular velocity and axial component of gravity force. The blade is modelled by a cantilever beam with variable cross-sections. One end of the blade is fixed to a rigid hub with radius R, and the other is free. The blade is considered to rotate at an angular velocity $\mathrm{\Omega }$ . The blade undergoes flap-wise bending vibration, and the governing equation can be obtained by using the variational principle.

The blade being considered in this study is shown in Figure 2. The Z-axis is along the hub’s central axis; the X-axis is along the spanwise of the blade, and Y-axis follows the right-hand rule. The total kinetic energy of the horizontal axis wind turbine blade due to the flap-wise vibration is

$${\rm K}=K_1+{{\rm K}}_2+{{\rm K}}_3$$

(1)

The first term is due to the flap-wise bending, expressed as

$$K_1=\frac{1}{2} {\displaystyle \underset{0}{\overset{L}{\int }}\rho A\left(x\right){\left(\frac{\partial w(x,t)}{\partial t}\right)}^{2}dx}$$

(2)

The second expression of the kinetic energy results from the rotary inertia of the blade, as follow

$$K_2=\frac{1}{2}{\displaystyle \underset{0}{\overset{L}{\int }}[\rho I\left(x\right){\left(\frac{{\partial }^{2}w(x,t)}{\partial x\partial t}\right)}^2+\rho I^{\ast }\left(x\right){\mathrm{\Omega }}^2{\left(\frac{\partial w(x,t)}{\partial x}\right)}^2]dx}$$

(3)

where

$I^{\ast }=I\left(x\right){\text{cos}}^2\left(\varphi \right)$ and $\varphi$ is the pre-cone angle.

The last term of the kinetic energy is resulted from adding the effects of pre-cone and pitch angles, and comes from the fact that the axis of rotation is not parallel with the flap-wise direction of the blade. ⁸

$${{\rm K}}_3=\frac{1}{2}{\displaystyle \underset{0}{\overset{L}{\int }}\left[-\frac{1}{8}\rho A\left(x\right){\mathrm{\Omega }}^2\left(2\cos \left(2\varphi \right)+\cos \left(2\left(\varphi -{\theta }_p\right)\right)+2\cos \left(2{\theta }_p\right)+\cos \left(2\left(\varphi +{\theta }_p\right)\right)-6\right)w{\left(x,t\right)}^2\right]dx}$$

(4)

In the above equations: $\rho ,$ represents the blade density, $A\left(x\right)$ its cross-sectional area in the flap-wise direction at distance $x$ relative to the blade span, $I\left(x\right)$ blade moment of inertia, $L$ is the blade length, $\mathrm{\Omega }$ angular velocity of the blade, $\varphi$ pre-cone angle, ${\theta }_p$ pitch angle, $w(x,t)$ is the flap-wise bending displacement.

The first term of the kinetic energy equation (1) is due to the flexural vibration of the flap-wise direction and is identical to the Bernoulli, Timoshenko and Rayleigh beam theories. The second and third terms are due to the rotary inertia of the blade. These terms are ignored in the Euler-Bernoulli theory. The second term is found in the Timoshenko beam theory. The third term comes from the fact that the axis of rotation is not parallel with the flap-wise span by the effect of pre-cone and twist angle of the blade. ⁸

The potential energy of the HAWT blade is

$$U=U_1+U_2$$

(5)

The strain energy of the blade due to flap-wise bending is

$$U_1=\frac{1}{2} \underset{0}{\overset{L}{{\displaystyle \int }}}EI\left(x\right){\left(\frac{{\partial }^{2}w}{\partial x^2}\right)}^2\text{d}x$$

(6)

The potential energy due to the centrifugal force is

$$U_2=\frac{1}{2} \underset{0}{\overset{L}{{\displaystyle \int }}}T\left(x\right){\left(\frac{\partial w}{\partial x}\right)}^2\text{d}x$$

(7)

with

$$T\left(x\right)=\frac{1}{2} \underset{x}{\overset{L}{{\displaystyle \int }}}\rho A\left(x\right)\left[\left({\mathrm{\Omega }}^2\left(R+x \text{cos}\left(\varphi \right)\right)\right)-g\text{cos}\left(\text{Ωt}\right)\right]\text{d}x$$

(8)

where $T\left(x\right)$ is the axial force due to the centrifugal tension force at a distance $x$ from the centre of rotation. The influence of the centrifugal force on the static displacement of the turbine blade almost concentrates on the axial and flap-wise direction. $E$ is the modulus of elasticity, $g$ is the gravitational acceleration, $R$ is the hub radius.

The virtual work of the applied distributed force, $f\left(x,t\right)$ , is given by

$$\delta \mathrm{w}\left(t\right)={\displaystyle \underset{0}{\overset{L}{\int }}f\left(x,t\right)\delta w\text{d}x}$$

(9)

where, $f\left(x,t\right)$ is the external distributed load along the blade span.

Noting that the operation of variation is commutative with respect to both integration and differentiation. So, the variation of kinetic and potential energy can be written as

$$\text{δ}{\rm K}=\delta {\displaystyle \underset{t1}{\overset{t2}{\int }}\left(K_1+{{\rm K}}_2+{{\rm K}}_3\right)\text{d}t}$$

(10)

$$\text{δ}U=\delta {\displaystyle \underset{t1}{\overset{t2}{\int }}\left(U_1+U_2\right)\text{d}t}$$

(11)

Hamiton’s principle can be stated as

$$\delta {\displaystyle \underset{t1}{\overset{t2}{\int }}\left(K-U+\mathrm{w}\right)\text{d}t}=0$$

(12)

$$\delta w=0 at t=t_{1 }and t=t_2.$$

Thus,

The governing equation of motion is obtainable as

$$\rho A\left(x\right)\frac{{\partial }^{2}w}{\partial t^2}+\frac{{\partial }^2}{\partial x^2}\left( EI\left(x\right)\frac{{\partial }^{2}w}{\partial x^2}\right)-\frac{\partial }{\partial x}\left(T\left(x\right)\frac{\partial w}{\partial x}\right)-\frac{\partial }{\partial x}\left(\rho I\left(x\right)\frac{{\partial }^{3}w}{\partial x\partial t^2}\right)+{\mathrm{\Omega }}^2\frac{\partial }{\partial x}\left(\rho I^{\ast }\left(x\right)\frac{\partial w}{\partial x}\right)-k\left(x\right)w=f\left(x,t\right)$$

(13)

With the boundary conditions at x=0

$$w=0, \frac{\partial w}{\partial x}=0{|}_{x=0}=0$$

and at x=L

$$EI\frac{{\partial }^{2}w}{\partial x^2}{|}_{x=L}=0,$$

$$\frac{\partial }{\partial x}\left(EI\left(x\right)\frac{{\partial }^{2}w}{\partial x^2}\right)-\rho I\left(x\right)\frac{{\partial }^{3}w}{\partial x\partial t^2}+{\mathrm{\Omega }}^2 \rho I^{\ast }\left(x\right)\frac{\partial w}{\partial x}{|}_{x=L}=0$$

(14)

Free vibration analysis

There are many methods that may be used to find the free vibration analysis of rotating beams. Most of those methods are applied for beams with certain cross-sections. Therefore, they cannot be applied to the blade with various and complicated cross-sections. Approximate methods such as Galerkin residual approach can be adopted. The governing equation (13) together with approximate solution $w\left(x,t\right)$ is employed to formulate the weak form as

$$\begin{array}{l}{\displaystyle \underset{0}{\overset{L}{\int }}\rho A\left(x\right)}\frac{{\partial }^{2}w\left(x,t\right)}{\partial t^2}{\phi }_{i}dx+{\displaystyle \underset{0}{\overset{L}{\int }}\frac{{\partial }^2}{\partial x^2}}\left( EI\left(x\right)\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right){\phi }_{i}dx-{\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial }{\partial x}}\left(T\left(x\right)\frac{\partial w\left(x,t\right)}{\partial x}\right){\phi }_{i}dx\\ -{\displaystyle \underset{0}{\overset{L}{\int }}\frac{\partial }{\partial x}}\left(\rho I\left(x\right)\frac{{\partial }^{3}w\left(x,t\right)}{\partial x\partial t^2}\right){\phi }_{i}dx -{\displaystyle \underset{0}{\overset{L}{\int }}k\left(x\right)w\left(x,t\right){\phi }_{i}dx}+{\displaystyle \underset{0}{\overset{L}{\int }}{\mathrm{\Omega }}^2}\frac{\partial }{\partial x}\left(\rho I^{\ast }\left(x\right)\frac{\partial w\left(x,t\right)}{\partial x}\right){\phi }_{i}dx -{\displaystyle \underset{0}{\overset{L}{\int }}f\left(x,t\right){\phi }_{i}dx}=0\end{array}$$

(15)

where, ${\phi }_i$ is the trial function, and i =1,2,…N

The flap-wise deflection and the slop due to flexural vibration of the wind turbine blade are approximated by a weighted linear combination of trial functions ¹⁶

$$w\left(x,t\right)={\displaystyle \sum _{j=1}^N{\phi }_j\left(x\right)u_j\left(t\right)}$$

(16)

where, $u_j$ is the joint displacement, and j=1,2,…N

The second derivatives of the trial solution are expressed as

$$\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}={\displaystyle \sum }_{j=1}^N\frac{{\partial }^2{\phi }_j\left(x\right)}{\partial x^2}{u}_j\left(t\right)$$

(17)

$$\frac{{\partial }^{2}w\left(x,t\right)}{\partial t^2}={\displaystyle \sum }_{j=1}^N{\phi }_j\left(x\right)\frac{{\partial }^2{u}_j\left(t\right)}{\partial t^2}$$

(18)

substituting equations (14), (15) and (16) into equation (13) and integrating with respect to $x$ from 0 to $L$ , the general form of the matrix equation of an N degree of freedom system subjected to excitation force may be obtained with the following expression:

$$\left[M\right]\left\{\ddot{U}\right\}+\left[C\right]\left\{\dot{U}\right\}+\left[K\right]\left\{U\right\}=\left\{F\right\}$$

(19)

where $\left[M\right] and$ $\left[K\right]$ are positive-definite mass and stiffness matrices respectively, $\left[C\right]$ is the damping matrix, $\left\{U\right\}, \left\{\dot{U}\right\} \text{and} \left\{\ddot{U}\right\}$ are columns matrices of linear displacement, velocity and acceleration, respectively, $\left\{F\right\}$ a column matrix of the excitation forces or moments. The natural frequencies in the flap wise direction and the corresponding mode shapes would be calculated by assuming:

$$u_j=X_j{e}^{\text{i}\omega t} , {\dot{u}}_j=\text{i}{X}_j{e}^{\text{i}\omega t} , {\ddot{u}}_j=-{\omega }^2{X}_j{e}^{\text{i}\omega t}$$

(20)

where $X_j$ is a vector of order $n$ representing the modal function, $\text{i}=\sqrt{-1}$ , $t$ is the time variable, and $\omega$ is a constant identified to represent the frequency of flap wise vibration of the vector $X_j$ . Substituting equations (20) into equation (19), a set of ordinary differential equations of modal functions and frequencies are obtained. For free vibration, the blade equation (19) can be written as (Equation. (21)):

$$M\ddot{U}+C\dot{U}+KU=0$$

(21)

From equations (20), we get

$$\left(K+\text{i}\omega C-{\omega }^{2}M\right)X=0$$

(22)

If the damping matrix $C$ is neglected, the undamped modal equation can be expressed as (Equation (23)):

$$-{\omega }^{2}X+M^{-1 }KX=0$$

$$M^{-1}KX=\lambda X$$

(23)

where the eigenvalues ${\lambda }_j={\omega }_j^2$

The previous equations yield the $n$ eignsolutions:

$$\left({\lambda }_1,X_{1i}\right), \left({\lambda }_2,X_{2i}\right), \left({\lambda }_3,X_{3i}\right), \dots \dots \dots \dots \left({\lambda }_n,X_{ni}\right)$$

The governing equation of the flap-wise deformation equation (15) is transformed to the general form of the matrix equations for free vibration of N degree of freedom system as indicated in equation (23) which it is represented the fundamental equation to form the pertinent equations to determine the natural frequencies and mode shapes configuration of the wind turbine blades. The elements of mass and stiffness matrices are

$$M_{ij}={\displaystyle \underset{0}{\overset{L}{\int }}\rho A\left(x\right){\phi }_j\left(x\right){\phi }_i\left(x\right)dx}+{\displaystyle \underset{0}{\overset{L}{\int }}\rho I\left(x\right)}\frac{\partial {\phi }_j}{\partial x}\frac{\partial {\phi }_i}{\partial x}dx$$

(24)

$$K_{ij}={\displaystyle \underset{0}{\overset{L}{\int }}EI\left(x\right)}\frac{{\partial }^2{\phi }_i\left(x\right)}{\partial x^2}\frac{{\partial }^2{\phi }_j\left(x\right)}{\partial x^2}dx+{\displaystyle \underset{0}{\overset{L}{\int }}T\left(x\right)}\frac{\partial {\phi }_i}{\partial x}\left(\frac{\partial {\phi }_j}{\partial x}\right)dx-{\displaystyle \underset{0}{\overset{L}{\int }}k\left(x\right){\phi }_j{\phi }_{i}dx}-{\displaystyle \underset{0}{\overset{L}{\int }}{\mathrm{\Omega }}^2\rho I^{\ast }\left(x\right)}\frac{\partial {\phi }_i}{\partial x}\left(\frac{\partial {\phi }_j}{\partial x}\right)dx$$

(25)

where

$$k\left(x\right)=-\frac{1}{8}\rho A\left(x\right)\left(2\cos \left(2\varphi \right)+\cos \left(2\left(\varphi -{\theta }_p\right)\right)+2\cos \left(2{\theta }_p\right)+\cos \left(2\left(\varphi +{\theta }_p\right)\right)-6\right){\mathrm{\Omega }}^2$$

$I^{\ast }, T\left(x\right)$ as indicated in equations (3) and (8), respectively.

The global mass and stiffness matrices can be constructed once the elementary value of mass and stiffness from equations (24) and (25) were obtained.

Figure 1 shows the computational procedure that followed in order to programming the free vibration analysis.

Figure 1.

Free vibration computational procedure.

Applicability of Euler-Bernoulli, Rayleigh and Timoshenko beam theories in fap-wise deformation

If the first two terms of equation (13) are considered, and the rest parts are neglected, the Euler-Bernoulli beam’s governing equation for the forced flap-wise vibration of non-uniform and non-rotating beam is obtained as

$$\rho A\left(x\right)\frac{{\partial }^{2}w\left(x,t\right)}{\partial t^2}+\frac{{\partial }^2}{\partial x^2}\left( EI\left(x\right)\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right)=f\left(x,t\right)$$

(26)

If the third term of equation (13) is added, the Euler-Bernoulli rotating beam’s governing equation for the forced flap-wise vibration of non-uniform beam is obtained as

$$\rho A\left(x\right)\frac{{\partial }^{2}w(x,t)}{\partial t^2}+\frac{{\partial }^2}{\partial x^2}\left( EI\left(x\right)\frac{{\partial }^{2}w(x,t)}{\partial x^2}\right)-\frac{\partial }{\partial x}\left(T\left(x\right)\frac{\partial w(x,t)}{\partial x}\right)=f\left(x,t\right)$$

(28)

where, $f\left(x,t\right)$ is the external distributed load along the beam span.

The governing equation for a rotating Rayleigh beam can be expressed directly from equation (13) by neglecting the terms of the effect of pre-cone and pitch angles to be as

$$\rho A\left(x\right)\frac{{\partial }^{2}w}{\partial t^2}+\frac{{\partial }^2}{\partial x^2}\left( EI\left(x\right)\frac{{\partial }^{2}w}{\partial x^2}\right)-\frac{\partial }{\partial x}\left(T\left(x\right)\frac{\partial w}{\partial x}\right)-\frac{\partial }{\partial x}\left(\rho I\left(x\right)\frac{{\partial }^{3}w}{\partial x\partial t^2}\right)+{\mathrm{\Omega }}^2\frac{\partial }{\partial x}\left(\rho I\left(x\right)\frac{\partial w}{\partial x}\right)=f\left(x,t\right)$$

(29)

If the shear deformation and rotary inertia are considered, the procedure, presented by Timoshenko beam theory. The governing equation for the forced flap-wise vibration of non-uniform beam is obtained as ¹⁶

$$\frac{{\partial }^2}{\partial x^2}\left( EI\left(x\right)\frac{{\partial }^{2}w(x,t)}{\partial x^2}\right)+\rho A\left(x\right)\frac{{\partial }^{2}w(x,t)}{\partial t^2}-\rho I\left(x\right)\left(1+\frac{E}{\mathit{\varkappa }G}\right)\frac{{\partial }^{4}w}{\partial x^2\partial t^2}+\frac{{\rho }^{2}I}{\mathit{\varkappa }G}\frac{{\partial }^{4}w}{\partial t^4}+\frac{EI}{\mathit{\varkappa }AG}\frac{{\partial }^{2}f}{\partial x^2}-\frac{\rho I}{\mathit{\varkappa }AG}\frac{{\partial }^{2}f}{\partial t^2}=f\left(x,t\right)$$

(30)

For free vibration equation (29) reduces to

$$\frac{{\partial }^2}{\partial x^2}\left( EI\left(x\right)\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}\right)+\rho A\left(x\right)\frac{{\partial }^{2}w\left(x,t\right)}{\partial t^2}-\rho I\left(x\right)\left(1+\frac{E}{\mathit{\varkappa }G}\right)\frac{{\partial }^{4}w}{\partial x^2\partial t^2}+\frac{{\rho }^{2}I}{\mathit{\varkappa }G}\frac{{\partial }^{4}w}{\partial t^4}=0$$

(31)

where $G$ denotes the modulus of rigidity of the material and $\mathit{\varkappa }$ is a constant, known as Timoshenko’s shear coefficient, which depends on the shape of the cross-section.

Shear correction value of wind turbine blade

Timoshenko beam theory is based on the concept of a characteristics response to shear force. The theory adopts the correction shear factor that compensates for the assumption of a constant shear strain over the cross-section. The shear factor is a reciprocal value that primarily depends on the shape of structure cross-section, on the other hand, Timoshenko suggested constant value of shear coefficient depending on the shape of cross-section; for the symmetric beam with a rectangular section, it is 0.83, whilst with a circular cross-section it is 0.9. ¹⁷ This correction value is defined as the average shear strain within a section to the shear strain at its centroid. ¹⁸ The study ¹⁹ adopted many methods to evaluate the shear correction values. All the values are estimated according to the aspect ratio of the cross-section; the value of the circular cross-section is 0.9 and this value decreases with decreasing the ratio of aspect ratio to 0.286. ¹⁹

The primary design drivers for large wind turbine blades are blade fatigue, as wind blades are required to be certified for 20 years of life span. Moreover, the natural modes of vibration, which are determined by its physical behaviour are inherent to the dynamic system, need to be predicted with accurate values. The forces on the blade element are assumed to calculate by means of the dimensional of each aerofoil characteristic along the spanwise direction of the blade. ²⁰ Based on the technical report, ²¹ the blade is divided into 17 different elements considering the coincident aerofoil in each element, ²² as shown in Figure 2. According to the Composite Material Fatigue Database, ²³ the elastic material property data and mass density are listed in Table 1. Properties for triaxial material, which are denoted as SNL Triax, were determined by averaging the test-derived data for the uni-axial and double bias material. ²³ Splines were drawn for mass and stiffness along the blade span according to data, ²¹ as shown in Figure 3 and Figure 4 are clarified the complicated distribution of blade mass and stiffness along the blade span. In this study, the correction shear factor starting at the value of 0.9 as the first section is circle cross-section and decreasing according to the chord and thickness of different blade aerofoils. ²¹

Figure 2.

Division of the blade into 17 elements. ²¹

Table 1.

Material Property Data Selected from Doe/MSU Database. ²³

Laminate Definition	Volume fraction %	$E_L$ (GPa)	$E_T$ (GPa)	${\upsilon }_{LT}$	$G_{LT}$ (GPa)	p (Kg/ $m^3$ )
E-LT-5500/EP-3	54	41.8	14.0	0.28	2.63	1920
Saertex/ep-3	44	13.6	13.3	0.51	11.8	1780
SNL Triax		27.7	13.65	0.39	7.2	1850

$E_L$ and $E_T$ : longitudinal and transverse modulus of elasticity, $\upsilon :$ Poisson’s ratio, G: shear modulus, $\rho$ : mass density.

Figure 3.

the mass per unit length (ρA) distribution along the spanwise of NREL 5MW blade. ²¹

Figure 4.

the flap-wise flexural rigidity (EI) distribution along the spanwise of NREL 5MW. ²¹

Comparison and correlation between mode shapes

It is possible to identify the structure condition from the variation of modal-based indices such as the frequency and mode shapes. The aim is to obtain the baseline of structure modal parameters when it is perfectly healthy. However, when any changes of these parameters occur, it may be possible to derive the causes which bring about the change.

The mode shapes are spatially distributed characteristics that include information about the structure. The first index for quantifying the comparison between the pairs of mode shapes is Modal Scale Factor(MSF). This technique represents the slope of the individual points of mode shapes which these points should be close to the straight line through which they will have a slope of $\pm 1$ . The MSF index for the mode shape and its counterpart is as ²⁴

$$MSF\left(X,A\right)=\frac{{\displaystyle {\sum }_{j=1}^N{\left(\text{ΨX}\right)}_j{\left(\mathrm{\Psi }A\right)}_j}}{{\displaystyle {\sum }_{j=1}^N{\left(\mathrm{\Psi }A\right)}_j{\left(\mathrm{\Psi }A\right)}_j}}$$

(32)

Where $\text{ΨX},\text{ΨA}$ are the mode shapes vectors that used for comparing, N is the number of degrees of freedom for both $X,A$ , j is the node number, and depending on which mode is taken as the reference one:

$$MSF\left(A,\text{X}\right)=\frac{{\displaystyle {\sum }_{j=1}^N{\left(\mathrm{\Psi }A\right)}_j{\left(\text{ΨX}\right)}_j}}{{\displaystyle {\sum }_{j=1}^N{\left(\text{ΨX}\right)}_j{\left(\text{ΨX}\right)}_j}}$$

(33)

If $\left\{\text{ΨA}\right\}\equiv \left\{\text{ΨA}\right\}$ then $MSF\left(X,A\right)=MSF\left(A,X\right)=1,$ and the two modes will be perfectly correlated.

In the second case, they differ by a simple scalar multiplier $\left\{\text{ΨA}\right\}=\gamma \left\{\text{ΨA}\right\}$ and

$MSF\left(X,A\right)=\gamma$ whilst $MSF\left(A,X\right)=\frac{1}{\gamma }$

The second index is the Modal Shape Correlation Coefficient (MSCC) and this provides the correlation of two pairs of the mode shapes.

$$MSCC\left(A,X\right)=\frac{{\left|{\displaystyle {\sum }_{i=1}^N{\left(\text{ΨX}\right)}_j{\left(\mathrm{\Psi }A\right)}_j}\right|}^2}{\left({\displaystyle {\sum }_{i=1}^N{\left(\text{ΨX}\right)}_j{\left(\text{ΨX}\right)}_j}\right)\left({\displaystyle {\sum }_{i=1}^N{\left(\mathrm{\Psi }A\right)}_j{\left(\mathrm{\Psi }A\right)}_j}\right)}$$

(34)

The $MSCC$ values of two mode shapes are always between 0 and 1. If the two modes are perfectly correlated, then $MSCC=1,$ and for different modes, a value close to 0 should be obtained due to the orthogonality of the mode shapes.

Results and discussion

The flowchart of the algorithm that used in this study is shown in Figure 1. The National Renewable Energy Laboratory (NREL) 5MW reference wind turbine is chosen to validate the previous equations and analysis. It is an upwind three-blade horizontal-axis multi-megawatt wind turbine. Based on the technical report ²¹ the blade is divided into 17 different elements considering the coincident aerofoil in each element, as shown in Figure 2. The rotor radius and the hub diameter of the blade are about 63m and 3m, respectively. The pre-cone and pitch angles are ${2.5}^{°}$ and $0^{°}$ , respectively. And the rotor speed is 12.1 RPM ²¹ . According to the Composite Material Fatigue Database, ²³ the elastic material property data and mass density are listed in Table 1. Properties for triaxial material, which are denoted as SNL Triax, were determined by averaging the test-derived data for the uni-axial and double bias material. ²³ Splines were drawn for mass and stiffness distributions according to data, ²¹ as shown in Figures 3 and 4.

The first 10 natural frequencies results of a single 5MW NREL blade with adopting different theories at zero rotation speed and without aerodynamic forces are compared together and with data from^8,13,25 are listed in Table 2. The results showed a good agreement of five models with the literature results, which indicated the appropriateness of the presented models using the finite element code. ²⁶

Table 2.

Natural frequencies (Hz) for the first 10 modes at angular speed V =0 RPM Rayleigh beam, rotating Bernoulli beam, non-rotating Bernoulli beam and Rayleigh beam with adding the pitch and pre-cone angles effects and results from Refs. 8,13,25 of NREL-5MW HAWT blade.

Method Mode NO.	Natural freqs With blade parameters Present	Rayleigh (Hz) Present	BModes(Hz) NREL	FAST(Hz) NREL	[13] (Hz)	[25] (Hz)	[8] (Hz)	Bernoulli Present	Rot-Bernoulli Present	Timoshenko Present
1	0.680	0.680	0.69	0.68	0.673	0.68	0.727	0.681	0.681	0.678
2	1.985	1.985	2.00	1.94	1.926	1.98		3.060	3.059	1.960
3	4.543	4.543	4.69	4.43	4.427	4.66		3.915	3.914	4.427
4	8.132	8.132						8.168	8.167	7.808
5	12.674	12.674						12.763	12.762	11.980
6	18.031	18.031						18.210	18.209	16.774
7	24.214	24.214						24.529	24.528	22.185
8	31.323	31.323						31.856	31.855	28.199
9	39.565	39.565						40.160	40.159	36.876
10	48.194	48.194						49.974	49.973	40.022

The data in Table 2 shows the first 10 natural frequencies of the NREL 5MW at zero rotational speed. The comparison is made between the natural frequencies by applying the different beam theories: Rayleigh, rotating Bernoulli, non-rotating Bernoulli beams, and with the blade parameters such as pitch, pre-cone angles and gravitational force component, for the same geometry and material properties of NREL-5MW HAWT blade. The comparison is made against the BModes and FAST codes, and from,^8,13,25 it was found that the results obtained by the present study agree well with the numerical results calculated from other investigations. The fundamental frequencies by adopted different theories in this study have approximately the same value, which indicates that ignoring the influence of centrifugal force at the stationary value, lead to the same results in the fundamental mode. In the same Table 2, there is a significant change in the second and third frequencies for both rotating and non-rotating Bernoulli theory. According to equations (26) and (28) the rotary inertia effect is not included, and, as a result, it has a significant impact on the second and third modes. However, this impact becomes less in the higher modes.

To investigate the effect of hub radius and rotor angular velocity on the natural frequencies and corresponding mode shapes of the blade equations (24) and (25) are utilised. The results in Table 3 show that the natural frequencies in the first 10 modes are increased slightly and barely noticeable with increasing the hub radius from 0 to 10m at a constant angular velocity = 12.1RPM. According to equation (8), it could be interpreted as increasing of hub radius leads to change in the centrifugal force.

Table 3.

The effect of increasing hub radius (0–6 m) on the first 10 natural frequencies (Hz) at angular speed Ω =12.1 RPM of NREL-5MW HAWT blade.

Radius (m) Mode No	0	2	3	4	5	6
1	0.762957	0.766002	0.767519	0.769034	0.770545	0.772053
2	2.071668	2.075022	2.076697	2.07837	2.080042	2.081713
3	4.626922	4.630292	4.631975	4.633659	4.635341	4.637023
4	8.214347	8.217784	8.219502	8.221219	8.222936	8.224653
5	12.75037	12.7537	12.75537	12.75703	12.75869	12.76036
6	18.10221	18.10547	18.1071	18.10873	18.11036	18.112
7	24.27965	24.28283	24.28443	24.28602	24.28762	24.28921
8	31.38531	31.38849	31.39008	31.39167	31.39326	31.39485
9	39.62018	39.62329	39.62484	39.62639	39.62795	39.6295
10	48.25511	48.2586	48.26034	48.26209	48.26383	48.26558

The second parameter, whose variation results in a considerable effect, is the angular velocity. The results in Table 4 show the influence of angular velocity from 0 to 25 RPM on the first 10 natural frequencies of the Rayleigh model. The angular velocity directly affects the centrifugal force and rotary inertia and, in general, on the stiffness matrix, as indicated in equations (13) and (25). In the case of the 5MW wind turbine, where the blade rotating at 12.1 RPM, Table 4 shows that the fundamental mode frequency increases by 10% compared with the fundamental frequency for a non-rotating blade at rotating speed 0 RPM in the same Table 4. Typically, centrifugal stiffening results in an increase of the first natural frequency in the out-of-plane direction. This effect becomes low in the higher modes, as shown in the same table. It indicates that the influence of centrifugal stiffening becomes progressively less in the higher modes.

Table 4.

The first 10 frequencies (Hz) VS the changing in rotor speed Ω (0–25 RPM) of NREL-5MW HAWT blade.

Angular speed Mode No	0	2.5	5	10	12	12.1	15	20	25
1	0.680874	0.684626	0.695748	0.738362	0.762118	0.765242	0.803784	0.886296	0.980913
2	1.98517	1.988989	2.0004	2.045396	2.071327	2.074184	2.118209	2.215993	2.335441
3	4.543338	4.547043	4.558138	4.602226	4.627877	4.629449	4.674693	4.774116	4.898674
4	8.132802	8.136483	8.147513	8.191461	8.217115	8.216925	8.264097	8.36454	8.491616
5	12.67489	12.6784	12.68893	12.73093	12.7555	12.75287	12.80056	12.89723	13.02018
6	18.03164	18.03505	18.04528	18.08616	18.11008	18.10465	18.15402	18.24849	18.36905
7	24.21415	24.21748	24.22746	24.26733	24.29068	24.28204	24.3336	24.42601	24.5442
8	31.32392	31.32723	31.33714	31.37675	31.39996	31.38769	31.44263	31.53459	31.65236
9	39.56532	39.56853	39.57814	39.61654	39.63904	39.62251	39.68041	39.76953	39.88364
10	48.19454	48.19817	48.20906	48.25262	48.27818	48.25772	48.32523	48.42686	48.55748

The effects of the blade parameters such as the pitch and pre-cone angles are shown in Table 5 and Table 6. However, although these angles have contributions on the blade stiffness matrices, their effects on the natural frequencies are barely noticeable. Moreover, the minimum effects of pre-cone and pitch angles have approximately the same influence, which decreases in the higher modes.

Table 5.

Effects of pitch angle (−30 to 30) on the first ten natural frequencies (Hz) at angular speed Ω =12.1 RPM of NREL-5MW HAWT blade.

Pitch angle Mode No	−30	−20	−10	0	10	20	30
1	0.758405	0.761911	0.764193	0.764984	0.764193	0.761911	0.758405
2	2.071056	2.072284	2.073088	2.073367	2.073088	2.072284	2.071056
3	4.626638	4.627145	4.627478	4.627594	4.627478	4.627145	4.626638
4	8.213143	8.213396	8.213564	8.213622	8.213564	8.213396	8.213143
5	12.74748	12.74762	12.74771	12.74774	12.74771	12.74762	12.74748
6	18.09724	18.09731	18.09735	18.09737	18.09735	18.09731	18.09724
7	24.27222	24.27224	24.27226	24.27227	24.27226	24.27224	24.27222
8	31.37506	31.37506	31.37506	31.37506	31.37506	31.37506	31.37506
9	39.60662	39.60659	39.60657	39.60657	39.60657	39.60659	39.60662
10	48.23834	48.23832	48.23831	48.23831	48.23831	48.23832	48.23834

Table 6.

Effects of pre-cone angle (−5 to 5) on the first 10 natural frequencies (Hz) at angular speed Ω =0 RPM of NREL-5MW HAWT blade.

Pre-cone angleMode NO	−5	−3	−1	0	1	3	5
1	0.764733	0.764861	0.764926	0.764934	0.764926	0.764861	0.764733
2	2.073275	2.073323	2.073347	2.07335	2.073347	2.073323	2.073275
3	4.627554	4.627575	4.627586	4.627587	4.627586	4.627575	4.627554
4	8.2136	8.213612	8.213618	8.213619	8.213618	8.213612	8.2136
5	12.74773	12.74773	12.74774	12.74774	12.74774	12.74773	12.74773
6	18.09736	18.09737	18.09737	18.09737	18.09737	18.09737	18.09736
7	24.27226	24.27227	24.27227	24.27227	24.27227	24.27227	24.27226
8	31.37506	31.37506	31.37506	31.37506	31.37506	31.37506	31.37506
9	39.60656	39.60657	39.60657	39.60657	39.60657	39.60657	39.60656
10	48.23831	48.23831	48.23831	48.23831	48.23831	48.23831	48.23831

Figure 5 shows the first six mode shapes for Rayleigh, rotating Bernoulli, non-rotating Bernoulli and Timoshenko beams, and blade parameters such as pitch, pre-cone angles, and gravitational force component at angular velocity = 0 RPM. The mode shapes from different theories are superimposed on the same figure to show the shape comparison. The figures show that the fundamental mode has the same path line in all adopted theories. The discrepancy was clearly noticed in the second and third modes. However, this agrees with the discrepancy in natural frequencies in the same modes and confirms that the rotary inertia significantly impacts the second and third modes. As in Table 2 and Figure 5, it can be seen from mode number four that the effect of rotary inertia is significantly decreasing, and its influence on the dynamic characteristics become less important in all modes higher than mode four. In Figure 5 it is possible to see that not only the degree of corresponding between the pairs of mode shapes is achieved, but also the nature of any discrepancies which do exist. It can be pointed out that the blade parameters such as pitch and precone angles have not affected the mode shapes. Furthermore, the rotating and non-rotating- Bernoulli methods have also the same mode shapes, which indicate that the existence of rotary inertia with their counterparts theories has an impact on the mode shapes. This elucidates that the influence of adding centrifugal force has no effect at a stationary value of angular velocity, whilst the rotary inertia has an impact on the first three of natural frequencies and mode shapes.

Figure 5.

Mode shapes of the first six modes of NREL-5MW HAWT blade at rotational speed=0 RPM.

Table 7 shows the first 10 natural frequencies by the Rayleigh method and with adding the effect of pitch and pre-cone angles as functions with increasing rotor speed. It can be noticed that the influences of pitch and pre-cone angles on the first 10 natural frequencies are barely noticeable. On the complex structure such as the blade in this study, the natural frequencies are closely spaced, so it is more difficult to identify the discrepancies or the correlation between the modes and their counterparts. Figure 6 shows the effect of increasing the angular speed on the first four mode shapes in the case of translational and angular mode shapes. There is a significant impact of rotating speed on the first and second mode shapes. The deviations of the points from the expected line are systematic along the blade span in modes three and four in translational mode shapes. At the same time, the increase in angular speed has more effect on modes three and four in the case of angular mode shapes. The additional information of mode shapes is represented in Figure 7. The lowest four angular mode shapes are shown in the figure for the Rayleigh blade with adding pitch and pre-cone effect at angular speed= 0 RPM. As in figure 5, the absence of the rotary inertia effect has a significant influence in modes two and three. While this effect shrinkages at mode four.

Table 7.

The first 10 natural frequencies (Hz) variation as a function of rotor speed Ω = (0–25 rpm) of NREL-5MW HAWT blade.

The first 10 natural frequencies with increasing the rotor speed of HAWT blade with effects of pre-cone& twist angles
rotating Speed Mode NO.	0	2.5	5	10	12	15	20	25
1	0.680874	0.684734	0.69617	0.739946	0.76432	0.807028	0.891473	0.988143
2	1.98517	1.9891	2.000843	2.047126	2.073786	2.121962	2.222351	2.344832
3	4.543338	4.547153	4.558576	4.603959	4.630356	4.678524	4.780768	4.90877
4	8.132802	8.136594	8.147957	8.193226	8.219648	8.268028	8.371428	8.502185
5	12.67489	12.6785	12.68936	12.73264	12.75795	12.80438	12.90396	13.03056
6	18.03164	18.03515	18.0457	18.08783	18.11249	18.15777	18.25511	18.37931
7	24.21415	24.21758	24.22787	24.26897	24.29304	24.33727	24.4325	24.55428
8	31.32392	31.32733	31.33755	31.37838	31.4023	31.44629	31.54107	31.66245
9	39.56532	39.56863	39.57854	39.61814	39.64134	39.68398	39.77585	39.89347
10	48.19454	48.19828	48.2095	48.25441	48.28076	48.32926	48.43402	48.56868
The first 10 natural frequencies with increasing the rotor speed of HAWT blade by Rayleigh-Ritz theory
rotating Speed Mode NO.	0	2.5	5	10	12	15	20	25
1	0.680874	0.684626	0.695748	0.738362	0.762118	0.803784	0.886296	0.980913
2	1.98517	1.988989	2.0004	2.045396	2.071327	2.118209	2.215993	2.335441
3	4.543338	4.547043	4.558138	4.602226	4.627877	4.674693	4.774116	4.898674
4	8.132802	8.136483	8.147513	8.191461	8.217115	8.264097	8.36454	8.491616
5	12.67489	12.6784	12.68893	12.73093	12.7555	12.80056	12.89723	13.02018
6	18.03164	18.03505	18.04528	18.08616	18.11008	18.15402	18.24849	18.36905
7	24.21415	24.21748	24.22746	24.26733	24.29068	24.3336	24.42601	24.5442
8	31.32392	31.32723	31.33714	31.37675	31.39996	31.44263	31.53459	31.65236
9	39.56532	39.56853	39.57814	39.61654	39.63904	39.68041	39.76953	39.88364
10	48.19454	48.19817	48.20906	48.25262	48.27818	48.32523	48.42686	48.55748

Figure 6.

The first, second, third and fourth mode shapes of NREL-5MW HAWT blade at different angular velocities.

Figure 7.

The first, second, third and fourth angular mode shapes of NREL-5MW HAWT blade at rotational speed=0 RPM

There is a range of blade pitch angles required for power control. Table 8, Table 9 and Table 10 show the effect of different pitch angles on the first 10 frequencies. Table 8 shows the influence of increasing the pitch angle on the natural frequencies at angular velocity = 0 RPM. The table shows fixed frequencies values for each mode that are exactly the same as the values of Table 2. At angular velocity = 0 RPM, there is no effect of centrifugal force on the natural frequencies as expected.

Table 8.

Effects of pitch angle (-30°to 30°) on the first ten natural frequencies (Hz) at angular speed Ω =0 RPM of NREL-5MW HAWT blade.

Pitch angle Mode No	−30	−20	0	10	20	30
1	0.680255	0.680255	0.680255	0.680255	0.680255	0.680255
2	1.983365	1.983365	1.983365	1.983365	1.983365	1.983365
3	4.539208	4.539208	4.539208	4.539208	4.539208	4.539208
4	8.125409	8.125409	8.125409	8.125409	8.125409	8.125409
5	12.66336	12.66336	12.66336	12.66336	12.66336	12.66336
6	18.01524	18.01524	18.01524	18.01524	18.01524	18.01524
7	24.19214	24.19214	24.19214	24.19214	24.19214	24.19214
8	31.29545	31.29545	31.29545	31.29545	31.29545	31.29545
9	39.52935	39.52935	39.52935	39.52935	39.52935	39.52935
10	48.15072	48.15072	48.15072	48.15072	48.15072	48.15072

Table 9.

Effects of pitch angle (−30°to 30°) on the first ten natural frequencies (Hz) at angular speed Ω =15 RPM of NREL-5MW HAWT blade.

Pitch angle Mode No	−30	−20	0	10	20	30
1	0.796578	0.801763	0.806295	0.805129	0.801763	0.796578
2	2.116362	2.118316	2.120033	2.11959	2.118316	2.116362
3	4.672611	4.673495	4.674271	4.674071	4.673495	4.672611
4	8.259576	8.260074	8.260512	8.260399	8.260074	8.259576
5	12.79214	12.79246	12.79274	12.79266	12.79246	12.79214
6	18.14084	18.14106	18.14126	18.14121	18.14106	18.14084
7	24.31483	24.315	24.31515	24.31511	24.315	24.31483
8	31.41746	31.41759	31.4177	31.41767	31.41759	31.41746
9	39.64772	39.64782	39.6479	39.64788	39.64782	39.64772
10	48.28517	48.28525	48.28532	48.2853	48.28525	48.28517

Table 10.

Effects of pitch angle (-30°to 30°) on the first ten natural frequencies (Hz) at anngular speed Ω =25 RPM of NREL-5MW HAWT blade.

Pitch angle Mode No	−30	−20	0	10	20	30
1	0.965086	0.976939	0.987245	0.984599	0.976939	0.965086
2	2.333463	2.338383	2.342701	2.341588	2.338383	2.333463
3	4.899912	4.902251	4.904307	4.903777	4.902251	4.899912
4	8.49193	8.493274	8.494456	8.494151	8.493274	8.49193
5	13.01707	13.01795	13.01871	13.01852	13.01795	13.01707
6	18.36145	18.36207	18.36261	18.36247	18.36207	18.36145
7	24.5311	24.53156	24.53196	24.53185	24.53156	24.5311
8	31.633	31.63335	31.63366	31.63358	31.63335	31.633
9	39.85668	39.85696	39.8572	39.85714	39.85696	39.85668
10	48.5241	48.52433	48.52452	48.52447	48.52433	48.5241

Table 9 shows the influence of changing the pitch angle at angular velocity = 15 RPM. The frequency values are increasing by increasing the rotational speed. Furthermore, these frequency values are slightly changing with changing the pitch angle, and this change becomes less and less in the higher modes. According to Table 10, the natural frequencies become higher by around 30% at angular speed = 25 RPM comparing to Table 8. As in Table 9, the changes in frequencies with changing the pitch angle are barely noticeable in the lower and higher modes at rotational speed =25 RPM.

Figure 8 shows the comparisons between the first 10 frequencies values for all the adopted theories against the frequencies of the Rayleigh theory for each of the modes included in the comparison. From the figure, it can be seen the degree of correlation between the counterparts of the results. Furthermore, it can be noticed that the Rayleigh model without the effect of control angles is closely correlated with the Rayleigh model with the effect of blade control angles. Furthermore, there is a high level of correlation in the Timoshenko model compared with the Rayleigh model.

Figure 8.

first 10 natural frequencies of the 5 MW wind turbine blade by different theories.

Whilst Figure 9 shows the natural frequencies at different rotation speeds versus counterparts frequencies at rotational speed = 0 rpm. As illustrated, the frequencies are coincident with the stationary values along the line 45°, which indicates no apparent discrepancies between them. However, some positive identification of each mode with its counterpart can be achieved clearly by the mode shape correlation method.

Figure 9.

frequencies of various rotational speed vs. the stationary values

The results of MSF are demonstrated in Figure 10 for the first three mode shapes in the case of rotational velocity = 012.1 and 15 RPM. Although the three modes in different pairs of modes are correlated and close to one for counterparts values, the figure shows the scalar multiplier γ becomes less when comparing the modes at 12.1 RPM rotating speed with 15 RPM. It reveals that this parameter gives small sensitivity of MSF index as it represents the slope of the best straight line through the pairs of comparing mode shapes points.

Figure 10.

The deviation of first three mode shapes by applying Modal Scale Factor index, blade at rotating speed=12.1 RPM, blade at rotating speed =0 RPM comparing with rotational speed =15 RPM.

Figure 11 shows the Mode Shape Correlation Coefficient index results of the pairs of the first three mode shapes for quantifying the comparison between them at different rotational speeds. The MSCC results are not significantly deviated which all values are close to unity. The MSCC values of mode shapes at 0 RPM counterpart of mode shapes at 15 RPM exhibit more deviation comparing with values of mode shapes at 12.1 rpm counterpart of mode shapes at 15 RPM. However, MSCC results are very close to the unity, and any small sensitivity may bring a change in the structural characteristics.

Figure 11.

The deviation of first three mode shapes by applying Mode Shape Correlation Coefficient index, blade at rotating speed=12.1 RPM, blade at rotating speed =0 RPM comparing with rotational speed =15 RPM.

An alternative approach to comparing the mode shapes is by making an x-y plot to represent each element in the mode shape vector, such as is shown in Figure 12. The first six mode shapes of the Rayleigh model with and without the effect of pitch and pre-cone angles are represented in this figure. Each individual type of point on this figure is related to a specific degree of freedom. It could be noted that all the points in all the first six modes lie close to a straight line passing through the origin, which indicates the two models are perfectly correlated.

Figure 12.

The comparison of 5 MW blade mode shapes, Rayleigh model with control angles effect vs Rayleigh model without control angles effect.

The data is represented in Figure 13 to compare the first six mode shapes of the Rayleigh and Timoshenko models in the case of the translational mode shapes Figure 13 a and angular mode shapes Figure 13 b. It is clear that there are discrepancies from the straight line in angular mode shapes. The deviations of the points from the expected line are systematic in some way in the second, three and fourth modes which indicates the existence of a correlation between them. There is a deviation in mode one, which infers that the effect of shear deformation may affect the angular mode shape due to the existence of shear angle. The discrepancy is increased in mode six in the same figure, which indicates the effect of shear in the higher modes.

Figure 13

. a. The comparison of 5 MW blade of translational mode shapes, Rayleigh model vs Timoshenko model b. The comparison of 5 MW blade of angular mode shapes, Rayleigh model vs Timoshenko model.

Figure 14 and Figure 15 represent the MAC for the translational and angular mode shapes data of the Rayleigh model comparing with adding the blade parameters such as pre-cone and pitch angle effect. The figures show that the data from the two models are perfectly correlated in the case of translational mode shapes Figure 14, and angular mode shapes Figure 15. From Figures 14 and 15, it could be noticed that only the first six modes are included. The MAC matrix is symmetric, and there are non-zero off-diagonal numbers that indicate that some of the mode shapes appear to exhibit a degree of correlation with each other, which may not be expected since the orthogonality between non-diagonal modes is predictable. However, the orthogonality property can only be translated to a perfect condition when all the degrees of freedom are included in the comparison.

Figure 14.

MAC index of translational mode shapes of 5 MW blade, correlation of Rayleigh and Blade with adding the control angles models.

Figure 15.

MAC index of angular mode shapes of 5 MW blade, correlation of Rayleigh and Blade with adding the pitch and pre-cone angles effect.

Figure 16 and Figure 17 represent the MAC matrices of the Rayleigh model with rotating speed = 0 RPM and 25 RPM in the case of translational and angular mode shapes. The correlation was slightly different in both figures. There are discrepancies in the matrix values around the diagonal comparing with the MAC matrix of the Rayleigh model in the case of translational and angular mode, as shown in Figures 16(b) and 17(b). It could be interpreted as the modes appear to exhibit a different degree of correlation between each other.

Figure 16.

a. MAC index of translational mode shapes of 5 MW blade, correlation of Rayleigh models at 0 RPM and 25 RPM.b, correlation of Rayleigh model with and without the blade control angles effect at 0 RPM.

Figure 17.

a. MAC index of angular mode shapes of 5 MW blade, correlation of Rayleigh models at 0 RPM and 25 RPM. b, correlation of Rayleigh model with and without the blade control angles effect at 0 RPM.

Figure 18 represents the MAC index of the rotating Bernoulli model comparing with the Rayleigh model at rotating speed = 0 RPM. In Figure 18 a, the matrix represents the translational mode shapes whilst 18 b represent the angular mode shapes. Both of them exhibit poor correlation in modes two and three. The MAC index of Timoshenko theory comparing with Rayleigh theory is showed in Figure 19 a and b in case of translational and angular mode shapes. The figures show a perfect correlation between the two models.

Figure 18.

a. MAC index of translational mode shapes of 5 MW blade, correlation of Rayleigh model and rotating Bernoulli model at 0 RPM.b, correlation of angular mode shapes of Rayleigh model and rotating Bernoulli model at 0 RPM.

Figure 19.

MAC index of 5 MW blade. a, correlation of translational mode shapes of Timoshenko and Rayleigh models; b, correlation of angular mode shapes of Timoshenko and Rayleigh models.

Conclusion

The flap-wise vibration of the NREL-5MW HAWT blade was investigated. All the blade parameters are included in this study, such as the centrifugal force, gravitational axial component, rotary inertia, pitch angle, pre-cone angle and hub radius. Subsequently, the Rayleigh, Timoshenko, rotating and non-rotating Bernoulli theories are used with the same model, same blade geometry, same material properties and the same choice of degrees of freedom as an effective way to make a clear comparison between them. The aim is to develop an accurate model that has clearly explained that for precise condition monitoring, a comprehensive dynamic model ensure the accuracy in determining the natural frequencies and corresponding mode shapes. The comparison procedures between the chosen models are based on the modal data. According to the obtained results, the following points are concluded:

1. Comparison of natural frequencies is made by tabulation of all chosen models results. The similarity of natural frequencies has been confirmed. The results clearly indicate that the free vibration analysis in the flap-wise direction is approximately the same in the case of the Rayleigh method and the model, which includes the blade control angles such as the pitch and pre-cone angles. Also, the results showed that increasing the angular velocity has a significant effect on the natural frequencies, whilst the influence of increasing the hub radius has less effect. The pitch and pre-cone angles have minimal influence on the natural frequencies. An important conclusion is that Bernoulli theory does not produce accurate results.

2. The second way to compare the theories is by plotting the natural frequencies for each mode included in the comparison. The frequencies values plot of the Rayleigh method with blade control angles effect versus the Rayleigh without the control angles effect are exactly coincident along the lone ${45}^{°}$ , which indicates that there are no discrepancies between them. On the other hand, there is some discrepancy in the frequencies plot of Timoshenko versus Rayleigh.

3. The comparison between the mode shapes is performed by plotting the deformed shape of each adopted theory along the blade span and overlaying the plots on each other for the first six mode shapes to show the points of similarities and differences. The mode shapes figures show a similarity in the fundamental mode where it has the same path line in all adopted theories. The discrepancy was clearly noticed in the second and third modes of Bernoulli and rotating-Bernoulli. However, this agrees with the discrepancy in natural frequencies in the same modes and confirms that the rotary inertia significantly impacts second and third modes.

4. Another way that is adopted to compare the different mode shapes is by making x-y plot. This comparison was convenient to verify the Rayleigh model, where all individual points on the plot lie close to the straight line passing through the origin. Also, there is the inconsistency of Timoshenko versus Rayleigh in the first six mode shapes plot in the case of angular mode shapes, which indicates the influence of shear deformation.

5. In the goal to reinforce the conclusion and to validate the most accurate model, many different types of comparison were achieved and not just to rely on one. The MSF indicate sensitivity as an index, and the deviation of the corresponding mode shapes can be adopted to compare the effect of different parameters. The comparison between the mode shapes using MAC provides a measure of similarity of the pair-wise eigenvectors/mode shapes. Using The MAC index gives the same conclusion, and the shortcoming of Bernoulli theory in the first and second modes is clearly noticed. The MAC index indicates the correlation in the Timoshenko mode shapes values compared with Rayleigh. Useful means of quantifying the degree of correlation between different theories. The Timoshenko theory has the influence of angular mode shapes, which indicates the effect of shear deformation to predict the vibrational behaviour of wind turbine blades and their dynamic characteristics.