Coupled bending-bending-torsion vibration of a rotating pre-twisted beam with aerofoil cross-section and flexible root by finite element method

The purpose of this paper is to extend a previously published beam model of a turbine blade including the centrifugal force field and root flexibility effects on a finite element model and to demonstrate the performance, accuracy and efficiency of the extended model for computing the natural frequencies. Therefore, only the modifications due to rotation and elastic root are presented in great detail. Considering the shear center effect on the transverse displacements, the geometric stiffness matrix due to the centrifugal force is developed from the geometric strain energy expression based on the large deflections and the increase of torsional stiffness because of the axial stress. In this work, the root flexibility of the blade is idealized by a continuum model unlike the discrete model approach of a combination of translational and rotational elastic springs, as used by other researchers. The cross-section properties of the fir-tree root of the blade considered as an example are expressed by assigning proper order polynomial functions similar to cross-sectional properties of a tapered blade. The correctness of the present extended finite element model is confirmed by the experimental and calculated results available in the literature. Comparisons of the present model results with those in the literature indicate excellent agreement.


Introduction
Pre-twisted beams with aerofoil cross-section are used in several types of engineering structures, such as turbine blades, helicopter blades, aircraft propellers, and wind turbine blades.If the centroid and shear center are not coincident as in the pre-twisted beam with an aerofoil cross-section, the flexural vibrations in two planes and torsional vibrations are inevitably coupled.Blade failure due to vibrations at or near a resonant condition has led to appreciable research in this field to avoid such undesired results.Many researchers have reported results on the vibrations of rotating pre-twisted beams of rectangular cross-section but few have discussed the same problem for an asymmetrical aerofoil cross-section.Furthermore, coupled bending-bending-torsion vibration of a blade with elastic support has received considerably less attention.Since there is no analytical solution for the vibration of rotating pre-twisted beam with an aerofoil cross-section, semi-analytical or numerical methods are utilized to solve this problem.
Houbolt and Brooks [11] developed the differential equations of motion for the lateral and torsional deformations of twisted rotating beams and used a Rayleigh-Ritz approach to compute the solution.Carnegie [3] formulated the total potential energy and kinetic energy of a rotating cantilever blade for small vibrations.Later, Carnegie [5] presented a derivation of the equation of motion (dynamics) of a pre-twisted cantilever blade mounted on the periphery of a rotating disc allowing for shear deflection, rotary inertia and torsion bending by use of the variational procedure.Also, Carnegie [6] reported a more detailed derivation for the kinetic energy of a rotating cantilever blade.Montoya [15] derived the equations of motion for a rotating pre-twisted cantilever beam with an aerofoil cross-section, including shear center and higher order effects, and solved these by Runge-Kutta numerical integration method.In his theoretical and experimental study, he obtained the natural frequencies of a blade with flexible root and showed satisfactory agreement between his experimental and calculated results.Fu [9] transformed the Carnegie's formulation for a rotating pre-twisted non-uniform Timoshenko beam in coupled bending-bending-torsion vibration into a set of recursion formulas as the basis of a lumped parameter approach.Karadag [12] developed a thick beam finite element with eighteen degrees of freedom including the shear center effect to determine the vibration characteristics of a rotating non-uniform aerofoil cross-sectioned blade.Abbas and Kamal [1] presented a finite element model having twenty-four degrees of freedom for the same problem, additionally, taking into account the root flexibility idealized as translational and rotational elastic springs.Sabuncu and Thomas [17] studied the vibration characteristics of pre-twisted aerofoil cross-section blade packets under rotating conditions using thin beam finite element model.As a new approximate method, Surace et al [18] presented an integral approach based on Green functions following Houbolt and Brooks' work.As a relevant contribution on this topic, Choi and Chou [8] proposed the modified differential quadrature method for vibration analysis of elastically supported turbomachinery blades, referring to the studies presented by Carnegie [4] and Fu [9] for the energy expressions.Using the mode expansion method, Lin et al [13] derived a closed form solution of the dynamic and static systems related to the differential equations for the coupled bending-bending vibration of a rotating non-uniform beam with tip mass, arbitrary pre-twist and an elastically restrained root.Recently, Yardimoglu and Inman [20,21] proposed a finite element model having fourteen degrees of freedom for coupled bending-bending-torsion vibration of a pre-twisted thick beam with varying aerofoil cross-section, following the Timoshenko beam finite element model for untwisted case [19].In the present paper, as an extension of the previous paper [20,21], the modifications due to the centrifugal force field and elastic root are presented in great detail.Considering the shear center effect, a geometric stiffness matrix due to the centrifugal force field is developed from the geometric strain energy [10] based on both large deflections [16] and the increase of the torsional stiffness [15] because of the axial stress.The elastic root of the blade is modelled by using the continuum model approach.Thus, the same finite element model formulation derived for a blade is employed to determine the stiffness and mass matrices of elastic root.An illustrative example is considered allowing a comparison to the experimental and theoretical results by Montoya [15] in order to demonstrate the performance, accuracy and efficiency of the present finite element model, and also to clarify the influence of the root flexibility effects on the natural frequencies.In this example, cross-section properties of the fir-tree root of the blade are expressed by assigning proper order polynomial functions similar to the cross-section properties of a tapered blade as reported in the previous paper.Comparisons of the present finite element results with the experimental and theoretical results in the literature (for simpler cases) indicate excellent agreement.

Rotational effects on finite element modelling
A blade with fir-tree root mounted on the periphery of a rotating disk as shown in Fig. 1 is considered.The view of the pre-twisted blade from tip toward root is also given in Fig. 2 to help visualize the configuration.The notation used through this paper is listed in Appendix A.
The longitudinal displacement of the beam is ignored [14] and the gyroscopic forces are neglected [10] to add the rotational effects on the finite element model derived in [20].The axial load acting on the cross-section at the distance R r , shown in Fig. 1, due to the centrifugal force field of the beam is obtained by using the following equation: Therefore, the axial load acting at any section Z due to the centrifugal force field of the beam segment between the section and the tip of the blade is written depending on the location of the section as follows: In Sections 4 and 6, the axial load P (Z) will be expressed in the element local co-ordinate system as P (z).However, the co-ordinate transformation based on Z = Z el + z is used in computer program.
On the other hand, the moment about the blade axis due to the axial stress arising from the centrifugal force is given in [15] by In this work, all parameters in Eq. ( 4) are functions of Z.

Root flexibility effect on finite element modelling
The root of a blade have a little flexibility due to the disk on which blade is mounted and the root attachment (such as fir-tree, T, pin-joint, etc) used to fix the blade onto the disc.These (disk and root attachment) can all be idealized as continuum models or discrete models to take the additional flexibility into account in finite element modelling.In this section, only the fir-tree type root attachment shown in Figs 1 and 3 is considered.
In the continuum model approach utilized in this step, the portion of the fir-tree type root attachment from R r to R b is treated in the usual manner used in the blade portion.Hence, cross-sectional properties of this portion can be formulated for finite element modelling by using proper order polynomial functions considering the dotted curve indicated in Fig. 3 for geometrical approximation.However, in the discrete model approach, the above mentioned portion is assumed as springs which offer flexibility in or about certain directions depending on the nodal freedoms chosen in finite element model.

R r R b
Fig. 3. Fir-tree root of the blade.

Strain energy
The elastic strain energy of a pre-twisted Timoshenko beam element derived in [20] is written as follows: where The symbol " " used in Eqs ( 5) and ( 6) represents differentiation with respect to z.
The geometric strain energy [10] of a pre-twisted beam element due to centrifugal force is written by considering Eqs (1) and (4) as follows: where in which r x (z), r y (z) and their differentiation r x (z) and r y (z) are given explicitly in [20].

Kinetic energy
Kinetic energy of an element is given in [20] as follows: where ) ) The overdot used in Eqs ( 12) to ( 16) is the usual compact notation for differentiation with respect to time.

Finite element formulation
The element employed in this section, derived in [20], has two nodes, with seven degrees of freedom at each node.The nodal variables are transverse displacements, cross section rotations and shear angles in two planes and torsional displacement.Hence, the nodal displacement vector is given by Therefore, the element displacement vector is expressed as follows Substituting the assumed polynomial functions for the displacements with coefficients expressed in terms of nodal displacements into Eqs (5), ( 7) and (11), the elastic strain energy and kinetic energy are written as follows: where and Here the matrices [k e ] and [m e ] are reported in reference [20].
The element geometric stiffness matrix is found by following the same procedure used in order to obtain the element elastic stiffness and mass matrices as follows: where  The polynomial vectors [P u ], [P ν ], and [P θz ] in Eq. ( 25) are given explicitly in [20] using the compact form [P dr ].All parameters in Eq. ( 25) are functions of z.

Dynamic equilibrium equation
To form the global matrices, element elastic stiffness, geometric stiffness and mass matrices given in the previous section are assembled in the usual way, then boundary conditions are applied as clamped at R r -free.In order to obtain the natural frequencies, the dynamic equilibrium equation is reduced to eigenvalue problem given below, The eigenvalue problem is then solved numerically.

Analysis and discussion
As an illustrative example, the rotating aerofoil cross-sectioned pre-twisted beam with rectangular cross-sectioned fir-tree root, shown in Figs 1 and 2, analyzed experimentally and theoretically (by means of Runge-Kutta numerical integration) by Montoya [15] is considered.Cross-section properties of the pre-twisted part of the blade given in tabular form at nine discrete cross-sections along the Z-axis in [15,20] are formulated by assigning the polynomial functions as given in [20] (and therefore not repeated here).However, the cross-sectional properties of the fir-tree root from R r to R b considering the dotted curve shown in Fig. 3 are given in Table 1 at five discrete cross-sections along the Z-axis.The cross-sectional data of the root are expressed as fourth-order polynomial functions determined by curve fitting.These functions are carefully checked by computing the coefficients of determination [7], which  are equal to unity, and also by inspecting the plots of functions and the data.The fixing radius of the blade denoted by R r is assumed to be at the narrowest section of the first scallop.For the sake of accuracy, the plot given for experimental and theoretical blade frequencies in terms of speed in [15] is scanned, and then this image is resized by utilizing an image editor to read the data as 0.25 Hz per pixel.Developing a computer program in MATLAB , confirmation of the present finite element model is achieved.The mass, elastic and geometric stiffness matrices in this program are computed by Gauss-Legendre numerical integration [2,22].
For sake of comparison the present results obtained by employing eight elements for pre-twisted part and one element for root part of the blade are computed.These results are then compared to the experimental and calculated results plotted by Montoya [15] in Tables 2 and 3, which compare the frequencies numerically.Also, the same numerical values are plotted and shown in Figs 4 and 5 to observe the trend of discrepancies visually.It is clear from Tables 2 and 4, and also from Figs 4 and 5 that the present results are in excellent agreement with the experimental and theoretical results reported in [15] for a simpler system.
Further validation of the present model taking the elastic root effect into consideration is accomplished by comparing the present results found for a non-rotating blade with rigid and elastic root to the results reported by Montoya [15] in Table 4. Again, excellent agreement can be seen from Table 4.

Conclusion
The present finite element is an excellent model for coupled bending-bending-torsion vibration analysis of a rotating varying aerofoil cross-sectional pre-twisted beam with flexible root.The comparisons of the present model results with previously published experimental and theoretical results show the superior performance of the present model.Also, the advantage of the present finite element formulation is the utilization a quite modest number of the degrees of freedom.Furthermore, the continuum model approach for elastic root of the pre-twisted blade in the finite element model is validated by considering the non-rotating blade frequencies.Linear velocities of centroid in xz and yz planes, respectively x, y Co-ordinate system through the shear center at root section X, Y, Z Global co-ordinate system z local co-ordinate distance measured along beam element

Fig. 5 .
Fig. 5. First harmonic frequencies of blade in terms of the rotation speed.

Table 4
Effect of root flexibility on frequencies (Hz) Polar moments of inertia about centroid and shear center, respectively I Gxx , I Gyy Area moments of inertia of the cross-section about Gxx and Gyy axes, respectively I Gxy Product moment of inertia of the cross-section about Gxx-Gyy axes I T Saint-Venant torsion constant J x , J y , J Coefficients of coupling about xx, yy and zz axes, respectively k Centrifugal force at the distance R r [P dr ] Polynomial vector for nodal variable dr (see Eq. (30) in [20]) {q} Global displacement vector {q n } Nodal displacement vector {q n1 }, {q n2 } First and second nodal displacement vectors {q e } Element displacement vector r x , r y Centroid co-ordinates with respect to shear center through xx, yy axes, respectively (see Fig. 2 in [20]) R b Minimum radius of pre-twisted portion (see Fig. 1) R r Minimum radius of elastic portion of attachment (see Fig. 1) Transverse displacement of the shear center in xz plane U e , U g Elastic and geometric strain energies ν Transverse displacement of the shear center in yz plane ν Gx , ν Gy