Determination of Natural Frequencies and Mode Shapes of a Wind Turbine Rotor Blades using Timoshenko Beam Elements

In the simulation of a wind turbine, the lowest eigenmodes of the rotor blades are usually used to describe their elastic deformation in the frame of a multibody system. In this paper, a finite element beam model for the rotor blades based on the transfer matrix method is proposed. Both static and kinetic field matrices for the 3D Timoshenko beam element are derived by numerical integration of the differential equations of motion using RUNGE KUTTA 4 th order procedure. The 10 beam reference axis in the general case is at an arbitrary location in the cross section. The inertia term in the motion differential equation is expressed using appropriate shape functions for the Timoshenko beam. The kinetic field matrix is built by numerical integration applied on the approximated inertia term. The beam element stiffness and mass matrices are calculated by simple matrix operations from both field matrices. The system stiffness and mass matrices of the rotor blade model are assembled in the usual finite element manner in a global coordinate system with the accounting for structural twist 15 angle and possibly pre-bending. The program developed for the above calculations and the final solution of the eigenvalue problem is accomplished using MuPAD, a symbolic math toolbox of MATLAB ® . The calculated natural frequencies using generic rotor blade data are compared with the results proposed from FAST and ADAMS software.


Introduction
Vibration of an elastic system refers to a limited reciprocating motion of a particle or an object of the system.Wind turbines operate in an unsteady environment which results in a vibrating response (see Manwell et al., 2009).They consist of long slender structures (rotor blades and tower), which result in resonant frequencies that should be taken into account during the initial design and construction.When the excitation frequency of the vibrating system is near any natural frequency, an undesirable resonant state occurs with large amplitudes, which may lead to damage or even the collapse of the wind turbine or its components.The vibration response, especially of the rotor blades, depends on the stiffness which is a function of the materials used, the design and the size (see Jureczko et al., 2005).Therefore, the estimation of nat-ural frequencies in the early design phase plays an important role in avoiding resonance.
The eigenmodes associated with the lowest natural frequencies are employed as shape functions to describe the elastic deformation of the rotor blade beam model in the frame of the usual simulation of the wind turbine as a multibody system.Generally, the first two bending eigenmodes in each direction (flapwise and edgewise) and optionally the two additional torsional eigenmodes are used.The determination of the lowest eigenmodes with sufficient numerical accuracy is the first step with respect to the modal superposition applied to the deformational motion of the rotor blades.
Due to the geometrical complexity of the blade cross section profiles and the extended use of composite materials, the exact calculation of natural frequencies in the design process takes a considerable amount of time and financial expense Published by Copernicus Publications on behalf of the European Academy of Wind Energy e.V. E. Stanoev and S. Kusuma Chandrashekhara: Natural frequencies and mode shapes of wind turbine rotor blades with regards to the 3-D modelling of the rotor blade using CAD software; hence a simplified finite element beam model is necessary.A twisted rotor blade is simplified into a cantilever beam with a non-uniform cross section.The structural twist angle is implemented by changing the orientation of the principal axis of the blade cross section along the length of the blade.
In the finite element formulation of beams two linear beam theories are established, the Euler-Bernoulli beam model and the Timoshenko beam model.Although the Euler-Bernoulli beam theory is widely used, the Timoshenko beam theory is considered to be better as it incorporates the effects of transverse shear and the rotational inertia on the dynamic response of the beam (see Kaya, 2006).In the classical approach of the finite element formulation for the free vibration analysis of beams, the stiffness and mass matrices are derived using interpolation functions derived from second-and fourth-order Hermite polynomials.The stiffness matrix is derived from the following equation (see Wu, 2013): where K (e) is the element stiffness matrix, B is the strain matrix and D m is the elasticity matrix for the beam.The element mass matrix of the beam (see Wu, 2013) is derived using the following equation: where M (e) is the element mass matrix, ρ is the mass density, v is the volume and a is the matrix of interpolation functions.
Using the above-mentioned standard relations and appropriate shape functions for the Euler-Bernoulli beam and the Timoshenko beam, the stiffness matrix and consistent mass matrix for the finite beam element can be derived.However, an alternative to this procedure, based on the transfer matrix method for the beam theory (see Graf andVassilev, 2006:69-88 andStanoev 2007) is developed in the present article.The element stiffness matrix can be derived on the basis of numerical integration of the first-order ordinary differential equation system for the differential beam element.The associated mass matrix can be developed by numerical integration of the inertia term in the differential equation of motion.The numerical integration results in static and kinetic field matrices, from which the element stiffness and mass matrices can be easily assembled.
In the present article, the above-mentioned procedure is used to determine the element stiffness and element mass matrix for the Timoshenko beam.The interpolation functions used for deriving the mass matrix are based on Hermite polynomials according to Bazoune and Khulief (2003).The system stiffness and mass matrices for the rotor blade are assembled in a global coordinate basis in the usual finite element manner.The numerical solution of the associated eigenvalue problem for the system without damping is computed using computer algebra software (in the frame of MATLAB ® ).
In Sects. 2 and 3, the proposed method is described in detail; in Sect. 4 the method is applied on a rotor blade structure with 288 degrees of freedom (DOFs).The results for the natural frequencies and the corresponding eigenmodes are compared with the results calculated using the FAST and ADAMS software.
2 Theoretical framework for the 3-D Timoshenko beam

Kinematic relationships
The general assumptions in the linear beam theory are as follows: a.The beam reference (longitudinal) axis is at any arbitrary location of the cross section.
b.The only stresses that occur are normal stresses σ x and shear stresses τ xy , τ xz .
c. Cross section planes, which are initially normal to the longitudinal axis, will remain plane after deformation.
The geometrical representation of the deformation state of a beam cross section is shown in the Fig. 1.The deformations u p , v p and w p at a cross-sectional point P are determined by three displacement functions u(x), v(x) and w(x) and three cross-sectional rotation angles ϕ x (x), ϕ y (x) and ϕ z (x) -all of them are a function of the beam axis coordinate x.The differential equation system is derived in accordance with Stanoev (2013): As seen in Fig. 1, the displacement vector u p can be expressed at any cross section point P as The three components of the strains occurring in the beam can be expressed as the derivatives of the displacement functions u p , v p and w p .The axial strain ε xx and the two shear strain components γ xz and γ xy are given by the following: where γ z and γ y are the constant shear strains that are not neglected in Timoshenko beam theory.

Principle of virtual work for internal forces
The virtual work components for internal forces corresponding to normal stresses and shear stresses are given by the following: where N is the axial force, M z and M y are bending internal moments, Q y and Q z are the corresponding shear forces and M TP is the St. Venant torsional moment.
With the introduction of the constitutive relations of Hooke's material law for the stresses corresponding to the internal forces in Eq. ( 6) and expressing the strains using Eqs.(4a)-( 4c) and (5a)-(5b), the virtual work relationship is reformulated as Here, A sz = α sz •A and A sy = α sy •A are the shear areas in the z and y directions respectively, α sz , α sy are the corresponding shear correction coefficients, A is the cross-sectional area, A y is the static moment with respect to the z axis, A z is the static moment with respect to the y axis, A yy is the moment of inertia with respect to the z direction, A zz is the moment of inertia with respect to the y direction, A yz is the deviation moment of inertia and I T is the torsional moment of inertia.
After the coefficient comparison in Eqs. ( 6) and ( 7) the internal forces corresponding to the normal stresses can be expressed by introducing the cross sectional stiffness matrix EA: The shear stress components in Eqs. ( 6) and ( 7) lead to the following relations: www.wind-energ-sci.net/4/57/2019/Wind Energ.Sci., 4, 57-69, 2019 The relation described in Eq. ( 10a) and (10b) implies that the chosen reference axis coincides with the shear centre due to neglected shear-torsion coupling terms in Eq. ( 7).
For the special case where the cross section coordinate system coincides with the principal axes, the deformation relationship in Eq. ( 8) reduces to

Differential equation system
The virtual work relation in Eq. ( 7) is rewritten for the case where the beam coordinate system coincides with the principal axis of the cross section -see Eq. ( 11).
After the partial integration of Eq. ( 12) the cross section deformation relations and differential equilibrium conditions for the Timoshenko beam element are compiled in a firstorder differential equation system (see Eqs. 13,14).For the Timoshenko beam with an arbitrary beam reference axis at any point on the cross section (see Eq. 8), the system of differential equations can be expressed in the following form: The differential equation system for the Timoshenko beam with the beam reference axis on principal axes can be represented in the following matrix form: The entries S ij in Eq. ( 13) are determined by the inversion of the cross-sectional stiffness matrix in Eq. ( 8):

Alternative finite element formulation
In the classical finite element formulation, the beam stiffness matrix and the consistent beam mass matrix are derived by developing an approach for the displacement functions through shape (interpolation) functions, which consist of first-and third-order Hermite polynomials.In this section, an alternative finite element procedure is presented, based on the numerical Runge-Kutta fourth-order integration of the differential motion equations.The integration of the static part (the coefficient matrix in Eqs. 13 and 14 respectively) leads to the static field matrix, while the integration of the inertia terms in the equation of motion Eq. ( 16) results in a kinetic field matrix.In the last step of the integration, using both field matrices, the respective element stiffness and the element mass matrices can be calculated using simple matrix operations.

The differential equations of motion
The differential equations of motion for the differential beam element (see Fig. 2) can be written in a matrix form according to Stanoev (2007) and Müller (2012) as follows: In the above-mentioned equation, matrix A is the coefficient matrix (see Eqs. 13, 14) and vector z = z 1 z 2 is the state vector, where T is the vector of the displacement functions, and (16a) T is the vector of internal (section) force functions.(16b) The vector b = b 1 b 2 contains the known excitation forces; however, for an eigenvalue problem b = 0.The coefficient matrix A and the state vector z in addition to excitation force b constitute the static part of the motion equation.The kinetic part of the motion equation can be expressed using a matrix of the interpolation functions and nodal acceleration vector as follows: where z1 = ü v ẅ φx φy φz T is the vector of ac- 6×12) is the matrix of interpolation functions (see Sect. 3.3), V (a) , V (b) R (6×1) is the vector with nodal accelerations and m R (6×6) is the inertia matrix of the differential beam element (see Sect. 3.2).

The inertia matrix term
The inertia matrix in Eq. ( 18) implies that the distributed mass µ (x) (kg m −1 ) is applied eccentrically at any location (y, z) in the cross section.The inertia matrix is expressed according to Stanoev (2013): where p , y and z in (kg m) are the mass moments of inertia for the cross section.

Shape functions for the Timoshenko beam element
The acceleration term z1 in Eq. ( 17) is expressed using the product of Hermite interpolating polynomials and the nodal acceleration vectors V (a) and V (b).Shape functions for axial and torsional deformations u (ξ ) and ϕ x (ξ ), respectively, are derived using first-order polynomial as follows: To express the coefficients a j in terms of the nodal displacements u (ξ = 0) = u a , u (ξ = 1) = u b or in terms of the  torsion ϕ x (ξ = 0) = ϕ x a , ϕ x (ξ = 1) = ϕ x b the following relations are applied to Eq. ( 20): Substituting Eq. ( 21) into Eq.( 20) results in the shape function for axial deformation and in the following function for torsional deformation ϕ x : The relationships in Eq. (10a), Eq. ( 11) and the corresponding part of Eq. ( 14) are a starting point to derive approximation functions for bending deformation in the xz plane: Using the above-mentioned relation the expression for w is given by The translational deformation function w (ξ ) is approximated by a cubic polynomial function: Using the constant shear strain relation in Eqs. ( 5b) and ( 26) the polynomial expression for constant shear strain can be deduced as follows: By including Eqs. ( 27) and ( 26) in Eq. ( 25) the polynomial expression for ϕ y (ξ ) results in To determine the coefficients c j the following boundary conditions are applied: The inversion of Eq. ( 29) yields The interpolation functions for w (x, y, z) and ϕ y (x, y, z), Eqs. ( 26) and ( 28), can be expressed by employing Eq. ( 30): where the products of both matrices N w and G w are introduced as functions H w j (j = 1, . .., 4).
A similar method is used for determining the approximation functions v (ξ ) and ϕ z (ξ ) for bending deformation in the xy plane.

Numerical integration
The special form of the numerical Runge-Kutta fourth-order integration method applied here is described in detail in Müller and Wolf (1975) and Schenk (2012).The integration operator is applied to the equations of motion in Eq. ( 16), i.e.
In order to gain sufficient numerical precision the beam axis needs to be divided into at least 20 integration intervals.The integration operator transfers the known state vector at beginning of the integration interval to the end of the interval.
The integration procedure starts with the state vector at the first node a, i.e. at location (x = 0): The integration operator is subsequently applied to the evaluated coefficient matrix A at each interval by excluding the initial state vector z (a).The result are static field matrices F (x, a), multiplicatively linked to z (a) see Eq. ( 40).Therefore, each F (x, a) matrix "transfers" the state vector at location (x = 0) to the end x of the integration interval considered (transfer matrix method).In the frame of the integration procedure the actual field matrix F (x, a) serves column wise as an initial vector for the next interval, and the components of the state vector z (a) remain excluded.The beam "load" vector b 1 b 2 , Eq. ( 38), evaluated in the actual interval, yields the  column in the F (x, a) matrix after integration -Eq.( 40).
The numerical integration of the inertia term b m in Eq. ( 38) is carried out column wise analogously to the "load" vector, by excluding the nodal accelerations V R (e, t) -the results are kinetic field matrices H (x, a) at the end of each integration interval (at location x) (see Eq. 40):  + This type of numerical integration allows (slightly) varying values of the coefficients of the A matrix, the b m inertia term and the b vector along the beam axis -i.e.all stiffness, mass and external load values of the beam element may vary.After the last integration step at the second node b, at location (x = L), static L (e) and kinetic H (e) field matrices are obtained: where According to Eqs. ( 13) and ( 14) the state variable z 1 represents six-component displacement vectors for V (a) and V (b), respectively, and the state variable z 2 represents sixcomponent internal forces vectors S (a) and S (b), at locations (x = 0) and (x = L) respectively.
one can derive the element stiffness matrix K (e), the element mass matrix M (e) and the element forces and moments F 0 using simple matrix operations as shown in Eq. ( 43).Then the beam element relationships for the internal forces can be formulated as

Single masses at eccentric positions
The numerical integration according to Runge-Kutta, described in Sec.3.4, offers the possibility of including single load or mass quantities within a beam element.Single eccentric masses can be taken into account at the integration interval boundaries.In the local coordinate system of the beam element, at a general position vector x AE = y E z E T , the eccentric mass and the vector representation of dynamic equilibrium, Eq. ( 44a)-(44b), is as shown in Fig. 4. The beam reference axis is at point A, and vector V E represents the acceleration vector at the point of application (x, x AE ).With the help of the dynamic equilibrium conditions Eq. ( 44a)-(44b), additional inertia forces and moments due to eccentric mass can be determined (see Li, 2016): where N L and N R are internal force vectors on the left and right in differential proximity to the point at location x, M L and M R are internal moment vectors on the left and right in differential proximity to the point at location x (see Fig. 4), m E is the eccentric single mass, Ei are the mass moments of inertia of the single mass and φi are the angular accelerations at location x.The additional inertia matrix M E (m E , Ei , x, x AE ), derived from Eq. ( 44a)-(44b), is analogous to the inertia matrix due to distributed mass in Eq. ( 18).During the numerical integration within the beam element (see Sect. 3.4) an additional eccentric inertia term has to be added to the kinetic field matrix at the end x (the point of application of eccentric mass) of the corresponding integration intervalsee Eqs. ( 40) and (45).
Single masses do not usually appear in a rotor blade model, but the same finite element may be used for the modelling of wind turbine towers.In this case single masses within a finite beam element could represent bolted ring flange connections or the mass of any equipment, such as lifts etc.

The eigenvalue problem
The system matrices for a rotor blade beam model are assembled in the usual finite element manner employing the developed element matrices K (e) and M (e) from Eq. ( 43).
In the case of free damped oscillation, the linear homoge- where M R (n×n) is the system mass matrix, K R (n×n) is the system stiffness matrix, D R (n×n) is the system damping matrix and q(t) R (n×1) is the nodal displacement vector.
The system matrices are symmetric and positively definite for finite element structures.For a free undamped system, the equation of motion is reduced to By introducing the following solution approach which is given by q (t) = q e iω 0 t , q (t) = q (iω 0 ) 2 e iω 0 t (48) into the equation of motion (Eq.47) the following eigenvalue problem is obtained: where I is a unity matrix.The condition for non-trivial solution for Eq. ( 49) is given by The n-grade characteristic polynomial p ω 2 0k has n real solutions ω 0k , (k = 1, . .., n) (eigenfrequencies) and n associated eigenvectors qk , calculated from Eq. ( 49).For real-life tasks the solution is usually acquired by use of eigensolver software.

Numerical example
The programming code for the procedure and the graphic plots described/shown below was written in MuPAD, a symbolic math toolbox in MATLAB ® (see Kusuma Chandrashekhara, 2018).The code was verified using realistic data for a wind turbine rotor blade.The blade structural data belong to a 5 MW reference wind turbine designed for offshore development (Jonkman et al., 2009).The blade is 63 m in length and is divided into 48 beam elements.The blade structural data consist of distributed mass (m L ), blade extensional stiffness (EA), flapwise stiffness (EA zz ), edgewise stiffness (EA zz ), torsional stiffness (GI T ), flapwise mass moment of inertia y and edgewise mass moment of inertia ( z ).Due to the lack of shear stiffness data in Jonkman et al. (2009), the values of (GA sz ) and (GA sy ) (used for the calculations in Table 1) -the respective edgewise and flapwise shear stiffness -are estimated to be 10 % and 20 % of www.wind-energ-sci.net/4/57/2019/Wind Energ.Sci., 4, 57-69, 2019 extensional stiffness (EA), respectively.The values of the above-mentioned stiffness and mass moment of inertias are specified at spanwise locations along the blade pitch axis and about the principal axes of each cross section as oriented by a twist angle (γ ) defined in the input data.The twist angle is incorporated using the rotational transformation of each local element stiffness and mass matrices into the global coordinate system.The results of first three (flapwise and edgewise) eigenfrequencies calculated using the Timoshenko beam model (see Kusuma Chandrashekhara, 2018), are compared with the proposed results from FAST and ADAMS in Jonkman et al. (2009).The results are as shown in Table 1.
The mode shapes and the corresponding eigenfrequencies for the first flapwise and edgewise eigenmodes as well for two torsional eigenmodes are as shown in  In Table 2 the calculated natural frequencies for three different variants of the shear correction coefficients are shown, approximated as GA sz EA or GA sy EA ratios.The comparison with the frequencies calculated using the Bernoulli beam model outlines the tendency toward a stiffer structure due to the presupposed infinite shear stiffness in this case.Natural frequencies f Bern are 0.5 %-1.0 % higher on average than f Timin the (30 %, 60 %) case.The natural frequencies remain unchanged for both beam models for the purely torsional modes only.The reason for this is that the equations for torsion and bending are uncoupled (for the case of the principal axes, see Eq. 14) and remain the same in both models.

Conclusion and outlook
The proposed Timoshenko beam element in the 3-D description has been developed on the basis of the transfer matrix method.Both static and kinetic field matrices for the beam element are derived by applying a Runge-Kutta fourth-order numerical integration procedure on the differential equations of motion in a special way.Appropriate shape functions for the Timoshenko beam have been used to approximate the inertia forces in the motion differential equation.The beam element stiffness and mass matrices are assembled by matrix operations from the derived element field matrices.The usual finite element equations of motion for the rotor blade model are cast in the general case accounting for the structural twist angle and possible pre-bending.Therefore, in the general case the rotor blade beam model represents a polygonal approximated space curve.
For the sake of verification, the natural frequencies and associated eigenmodes are calculated using real-life rotor blade data with the incorporation of realistic twist angle data.The first two edgewise and flapwise eigenfrequencies obtained are compared with the proposed results from FAST and ADAMS software given in Jonkman et al. (2009).It can be observed that the deviation of the results of the Timoshenko beam model from FAST is comparatively lower and is in good agreement with FAST; thus, it can be stated that the approach presented regarding alternative finite element formulation works well.
One key input parameter for the Timoshenko beam model is the shear stiffness.As it was not the main goal of the present work to determine an appropriate shear correction coefficient for realistic rotor blade data, the numerical example was performed with a very rough approximation for GA sz and GA sy .This was done in order to simply demonstrate the performance and the differences to the Bernoulli beam model.However, if detailed data for the complex multilayer design of a rotor blade were available, a more realistic estimation of the shear stiffness could be expected.A workable method for the determination of the shear correction coefficient of a real-life rotor blade represents an important topic for further research.

Figure 2 .
Figure 2. The finite beam element: internal forces and nodal DOFs.

Figure 3 .
Figure 3. Definition of the dimensionless coordinate (ξ ) of a beam element.

Figure 4 .
Figure 4. Eccentrically applied mass m E at the x AE point of the beam in the 3-D case.

Table 2 .
Comparison between Timoshenko and Bernoulli beams with three variants for shear stiffness values.