1D finite element for modelling of turbine blade vibration in the field of centrifugal forces

The paper deals with the modelling of turbine blade vibrations by means of a novel 1D finite element that has only 16 degrees of freedom. Assuming linear elastic behaviour of the blade material and considering small displacements and strains, the derived blade finite element takes into account the effects of tension, torsion and bending in accordance with the Bernoulli’s hypothesis. Additionally, the finite element interlinks bending and torsion, and respects membrane forces acting on the blade. The derivation of matrices and vectors describing the blade finite element is provided in detail by using the Lagrange’s equations while the effect of membrane forces is included via the virtual work principle. For modelling purposes, the mathematical model of a turbine blade requires only the knowledge of cross-section contour points at several selected sections along the turbine blade axis. On the basis of these points, cross-section characteristics including the warping function are approximated along the blade axis by means of cubic splines. The advantage of this approach lies in the fact that all the blade cross-section parameters are identified before running numerical simulations. The warping function introduced in this paper and derived by variational principle describes cross-section warping caused just by torsion of a prismatic rod. For the verification of the proposed 1D finite element, an analysis of modal properties of the turbine blade M6 L-1 manufactured by Doosan Škoda Power is performed. This is achieved by comparing the lowest natural frequencies and corresponding mode shapes computed by the 1D and 3D models for a standing blade. The results revealed good agreement between both models despite the significant difference in their degrees of freedom. The applicability of the 1D finite element is further demonstrated by analyzing the dependence of natural frequencies on rotor speed. c © 2019 University of West Bohemia. All rights reserved.


Introduction
Many papers and books focused on problems concerning turbine blade dynamics and vibrations were published in the past. At the beginning of research, entire blades were modelled as a 1D continuum and described by equations of motion including equilibrium conditions of an infinitesimal element. The solution of derived partial differential equations describing motion of the blade was sought by approximate methods, because the use of computers was very limited. One of the earlier works published in the field of vibrations was written by Filipov [4]. Modelling of a blade packet by means of simple beams coupled with discrete massless springs was shown in paper [1], where each blade was modelled as a massless beam with mass concentrated at its ends. Here, the number of degrees of freedom (DOF) corresponded to a number of blades. This very simplified model takes into account only the first mode shapes of individual blades.
Progress in the finite element (FE) method led to application of the aforementioned modelling approach to turbine blade dynamics and related fields. Selected theoretical works dealt with interaction of a bladed wheel with steam aerodynamic forces [7] or with mutual interaction between individual blades caused by frictional forces [9]. The last mentioned work dealt with a very interesting approach that uses cyclic symmetry for the reduction of the DOF. This approach allows to respect even other nonlinear forces between individual blades, while the solution is performed in frequency domain. The recent papers [8] and [6] presented non-linear dynamic analysis of a friction member between blades. However, the used models were simple, e.g., the model mentioned in [8] is only a 2-DOF system.
The majority of current works considers models that use the FE analysis by means of 3D elements. Not only due to high computational costs, the complex models usually include only several blades, bladed packets or an entire bladed disc. To simulate the behaviour of a whole turbo-machine with minimal costs, another approach is necessary. For this purpose, the present paper introduces a novel 1D rotating FE model that despite its relatively low number of DOF is able to well approximate the behaviour of a 3D blade. The description of the finite element, which will be derived in the following sections, includes all beam properties and the pre-processor for warping function calculation and takes into account the variability of blade cross-section along the blade axis. Additionally, the following kinematical considerations are applied: • The cross-section projection into the plane perpendicular to the blade axis remains unchanged. • The kinetic and potential energies of the blade finite element respect warping of the cross-section. • The material of the blade finite element satisfies the Hooke's law.
• The material of the blade finite element is assumed to be homogeneous and isotropic.
• Bending deflections of the blade are considered to be small and corresponding rotations can be expressed as their derivatives. • The cross-section characteristics are approximated by cubic polynomials along the blade axis. • The blade model respects membrane forces.

Mathematical description of the 1D finite element
In this section, the mathematical background of the proposed 1D blade finite element is given. This is achieved firstly by introducing the kinematics of displacements, expressing the kinetic and deformation potential energies of the finite element and outlining the calculation of the crosssection warping function. Consequently, the matrices and vectors describing the finite element are derived by applying the Lagrange's equations and by assuming the effect of membrane force. Lastly, the equations of motion of one finite element as well as those of the entire blade are assemblied.

Kinematics of displacements
For a better understanding of the following sections, let us introduce a coordinate system ξ, η, ζ, where ξ is an arbitrarily chosen axis in the direction of the longest blade dimension (the blade axis mentioned above). The axis η is parallel with the vector of angular speed ω and the third axis ζ is completed in such a way to get a right-handed coordinate system shown in Fig. 1. Note that this system ξ, η, ζ defines the rotating coordinate system of the rotor. The coordinate system ξ, η, ζ rotates at the constant angular speed ω about the axis η which is identical to the axis of the turbine rotor. The intersection of ξ axis with the blade cross-section is marked as point P . Let us consider the beginning of the finite element in the distance r from the axis of rotation. Then, the motion of an arbitrary point of the cross-section at the coordinate ξ can be decomposed into the primary sliding motion of the reference point S * from the initial state point S (displacements of the cross-section shear centre u S , v S , w S ) and the secondary spherical motion of the rigid cross-section caused by bending (ψ, ϑ) and torsion (ϕ) and the tertiary motion in the x-direction caused by warping. Consequently, the radius-vector of the shear centre in coordinate system ξ, η, ζ can be expressed in the form where u 0 corresponds to the displacement of all cross-section points in the ξ direction caused by the tension and η S , ζ S and η C , ζ C denote the coordinates of the shear centre and the crosssection centroid, respectively. The symbols w 0 and v 0 denote the displacements in the ζ and η directions caused by the bending in planes ξζ and ξη, respectively. Due to the assumption of small displacements, the angles associated with the bending can be written as Taking the shear centre as a reference point, we can express the radius vector r of an arbitrary point of the blade cross-section in the coordinate system ξ, η, ζ in the following form wherer is the radius vector of the arbitrary point in the coordinate systems x, y, z and T is the transformation matrix from the x, y, z coordinate system to the ξ, η, ζ coordinate system T = ⎡ ⎣ cos ψ cos ϑ, sin ϕ sin ϑ − cos ϕ cos ϑ sin ψ, sin ϑ cos ϕ + cos ϑ sin ϕ sin ψ sin ψ, cos ψ cos ϕ, − cos ψ sin ϕ − cos ψ sin ϑ, cos ϑ sin ϕ + sin ϑ cos ϕ sin ψ, cos ϑ cos ϕ − sin ϑ sin ϕ sin ψ ⎤ ⎦ . (4) The radius vectorr has the formr which means that the overall warping corresponds to a product of the warping function g and the relative twist ϕ . For practical reasons, we express the velocity vector in the rotating coordinate system ξ, η, ζ asṙ where the angular speed vector ω = [0 ω 0] T corresponds to the rotation of the coordinate system ξ, η, ζ in relation to the stationary coordinate system.

Kinetic and deformation potential energies of one finite element
The kinetic energy in the form where ρ is the density, can be derived by substituting (3), (4), (5) into (6). Due to the complexity of the resulting relation, it is not presented in this paper. For the derivation of the deformation potential energy, the Bernoulli's beam hypothesis is taken into consideration. According to this hypothesis, the deformation tensor consists only of the following non-zero components where u is the displacement of the arbitrary point of the cross-section in the ξ direction and g denotes partial derivative of the warping function with respect to the coordinate ξ. By assuming that the warping function varies slowly along the coordinate ξ, we can neglect the last term of (8) 1 . The Hooke's law for the Bernoulli's beam theory takes the form where σ = [σ ξ τ ξζ τ ξη ] T is the stress vector, ε = [ε ξ γ ξζ γ ξη ] T is the deformation vector and E = diag(E G G) is a diagonal matrix, where E is the Young's modulus and G is the shear modulus. Then, the relation for the deformation potential energy of one finite element can be written in the following form or after substitution where is the so called torsion constant and The bar over subscripts in Dηζ denotes a quantity in the coordinate system with the origin in the cross-section centroid.

Calculation of the warping function
In this paper, the cross-section warping caused only by blade torsion is taken into account. The calculation of the warping function is performed supposing clear and free torsion of a prismatic rod. Let us assume that the cross-section rotates about the yet unknown point called the shear centre. Taking into consideration clear torsion, the first three terms at the right hand side of (8) 1 are equal to zero because these terms correspond to the deformations caused by tension and bending. The fourth term is zero because the torsion angle ϕ is linearly dependent on the coordinate ξ, and thus, ϕ = 0. The last term in (8) 1 is also equal to zero because the warping function g is constant in the case of clear torsion of the prismatic rod. Then, for the twisted prismatic rod, the potential energy of an infinitesimal element of length dξ can be written in the form 1 dξ where γ ξη and γ ξζ express shear strains given by relations (8) 2 and (8) 3 . The substitution of (8) into (14) and the addition of external effect yield where M ξ is the external torsional moment. Relation (15) contains the total potential energy of deformation and external forces E P . To obtain a unique solution of (15), it is necessary to add the following three conditions: ensuring zero normal force and two bending moments related to axes η and ζ. From the second and third conditions in (16), it is possible to get D ϕη = 0, D ϕζ = 0. The solution of the warping function will correspond to a minimum of the following functional where λ * 1 , λ * 2 , λ * 3 are the Lagrange's multipliers. In (17), there are seven independent quantities ϕ , g, η S , ζ S , λ * 1 , λ * 2 , λ * 3 because the shear centre coordinates are still unknown. Then, a functional variation takes the form Because the individual variations are independent of each other, the terms contained in these variations must be equal to zero, resulting in a system of seven equations Equation (19) is a torque equilibrium condition independent of the other ones. The remaining six equations (20)-(25) are used for the calculation of unknown quantities g, The FE discretization of (20)-(25) leads to a system of (n + 5) linear algebraic equations, where n corresponds to the number of discretization points of the cross-section area as well as to the number of function values of the warping function. Then, the vector of unknowns can take the following form where g ∈ R n,1 is the function value vector of the warping function at the discretization points.

Approximations of displacements and of cross-section characteristics
The displacements in the ξ, η, ζ directions and the corresponding torsion and bending angles ϕ, ψ, ϑ can be approximated as where Prime derivative for tensile deformations was used due to the variability of the blade profile and because it is appropriate to use higher order polynomials instead of linear ones. Because the use of quadratic polynomials without derivatives would require an extra mid-site node for each finite element, the present study applies cubic polynomials and considers derivatives at the ends of the elements. The same degree of approximation (cubic polynomials) are used for the approximation of cross-section characteristics along the coordinate ξ, leading to the following relations: where

Matrices and vectors
In (30), the vector of generalized displacements has the following form The symbolM e represents the FE mass matrix, ωG e is the FE matrix of gyroscopic effects,K e is the FE stiffness matrix,K e D is the FE circulation matrix andf De is the FE vector of centrifugal forces. For better clarity, the specific forms of these matrices and the vector of centrifugal forces can be found in the Appendix of this paper. Note that compared to the cross-section characteristics listed in (28), which depend on the coordinate ξ, the material parameters ρ, E and G appearing in the aforementioned matrices and the vector are assumed to be constant within one finite element.

Membrane generalized force vector
In this section, the effect of the membrane force is included into the mathematical model of the blade finite element. For this purpose, it is necessary to express the action of the force S(ξ) in the transverse directions η and ζ, see Fig. 2.
In accordance with Fig. 2, the force resultant in the direction η has the form By assuming the angle ψ to be small and taking (2) 2 into consideration, (31) can be rewritten as Analogously, the force resultant in the direction ζ can be derived, i.e., by considering the angle ϑ to be small and respecting (2) 1 , To be able to apply relations (32) and (33), the determination of the axial force S(ξ) has to be performed. As apparent from Fig. 3, the origin of the coordinate ξ for the blade finite element lies in the distance r from the axis of rotation. The axial force S o acting on the finite element from the rest of the blade can be expressed in the form S o = Δmeω 2 , where Δm is the mass of the remaining part of the blade in the direction ξ and e is its centroid eccentricity. Then, using the approximation of A(s) from (28), the axial force can be written as The integration of (34) results in and its corresponding derivative with respect to ξ By substituting (35) and (36) into (32) and (33), the membrane forces acting on the infinitesimal element in the directions η and ζ can be obtained In (37) and (38), the terms in curly brackets represent loads (force per unit length) caused by membrane forces. Their virtual work over the whole finite element can be expressed by relations The substitution of the approximations (27) into (39) and (40) and their integration over ξ yield Note that for the sake of clarity, special integrating matrices were introduced in (41) and (42) by taking the following substitutions into consideration Finally, the entire vector of membrane forces acting on a blade finite element can be written as whereK e M = ω 2Me M ∈ R 16,16 and the individual submatrices have the form and

Total equations of motion
With reference to Section 2.5, the total equations of motion for a blade finite element can be expressed asM eq e (t) + ωG eq e (t) wheref Ze (t) is the vector of remaining external forces. With respect to (A5) and (45), the matrixK e D and the vectorf Me in (48) can be modified leading to the equations of motion in the following form M eq e (t) + ωG eq e (t) From the practical point of view, the sequence of general displacements inq e (and forces) can be rearranged in such a way that the first 8 vector components correspond to the first node of finite element and the remaining 8 vector components to the second one. This rearrangement can be performed with the help of a permutation matrix J, i.e., all matrices and vectors in (49) are transformed as M e = J TMe J and q e = J Tq e and f Ze (t) = J Tf Ze (t) etc. Then, the total equations of motion for the blade finite element have the final form To mathematically describe the motion of the entire blade, it is necessary to assembly the total matrices and vectors in a way that is standard in the FE analysis. Then, the equations of motion of the blade can be obtained in the following form, similar to (50), These equations describe a transient vibration of the blade.

Verification of the proposed blade finite element on the basis of modal properties
At the beginning, it should be pointed out that the calculation of the warping function outlined in Section 2.3 was performed supposing clear and free torsion of a prismatic rod. However, the blade can be generally of variable cross-section and be exposed to constrained torsion. Nevertheless, for the derived blade finite element, we consider simplifying assumptions that the warping function is identical for free and constrained torsion as well as for prismatic and non-prismatic rods. It means that this function depends only on the cross-section shape. To be able to use the blade finite elements for modelling purposes, it is necessary to identify blade parameters a priori to running any numerical simulation. Firstly, for each known section (always perpendicular to the blade axis), the cross-section characteristics, the coordinates of the center of gravity and the shear centre and the warping function have to calculated. Consequently, these parameters are approximated by means of cubic polynomials along the blade axis and polynomial coefficients a ij in (28) are determined for each blade finite element. Finally, it is possible to assembly the corresponding matrices and vectors introduced in Section 2 and listed in Appendix. The verification of the proposed 1D blade finite element is carried out on the model of the turbine blade M6 L-1 manufactured by Doosan Š koda Power. In this case, the geometry of this blade is given by individual cross-section contour points obtained at 10 different sections along the turbine blade axis. Considering the zero angular speed (ω = 0), the modal behaviour of the turbine blade is analyzed by means of 1D and 3D models, particularly of interest is the behaviour corresponding to the lowest eigenfrequencies. For the 1D model based on the proposed finite element, all required input parameters are calculated as described in the previous paragraph. In the case of the 3D model, the geometry of the turbine blade is reconstructed using cubic splines and the resulting volume is discretized with tetrahedral finite elements. The modal analysis, which includes the calculation of natural frequencies Ω i , i = 1, 2, . . . , n (n ∼ blade DOF), is based on the following equations where v is the eigenvector. The boundary conditions prescribed in both 1D and 3D models correspond to a clamped low end and to a free upper end of the turbine blade. Although these boundary conditions do not faithfully simulate the actual blade attachment, they are sufficient to verify suitability of the proposed 1D finite element. In Fig. 4, the dependence of natural frequencies on the number of DOF is shown as obtained by the 3D model. As apparent, a quite high number of used DOF (here about 400 thousand DOF) is necessary to achieve reasonable results. Note that the figure also contains the natural frequencies calculated by the 1D model, which is, however, described only by 154 DOF and 19 finite elements, i.e., no dependence on the number of DOF is considered here. A quantitative   comparison between the results computed by the 3D and the 1D models is provided in Table 1.
The absolute difference is observed to be around 10 percent, which may seem to be high, but considering the difference in DOF (400 000 vs. 154), it is acceptable. Additionally, it should be noted that a numerical test with lower number of finite elements in the 1D model was also performed (10 finite elements with 81 DOF). For the first few natural frequencies, the differences between the 154 and 81 DOF were found to be less than one percent, making both models comparable. For illustration, Figs. 5 and 6 show the first two mode shapes as computed by the 3D and 1D models. The blade before deformation is depicted on the left of the figures, the central position corresponds to the mode shapes of the 3D model and the right position shows the mode shapes obtained by applying the presented methodology. In both cases, these figures indicate a good agreement between the mode shapes from both computational blade models.
In the remainder of this section, the behaviour of the blade in a rotating coordinate system (ω = 0) is studied. Because of the complex implementation of the 3D model (commercial software usually do not offer the possibility to include effects of the rotating system), only the blade model described by the 1D finite elements is assumed. Then, the corresponding mathematical model has the form given by (51) with natural frequencies dependent on the angular speed ω. Compared to the previous paragraphs, the calculation of the frequencies cannot be performed according to (52), but the system of n second order differential equations (51) has to be converted to a system of 2n first order differential equations by using the following trivial identity Next, let us consider only the eigenvalue problem which involves the solution of the extended system of equations where The natural frequencies Ω i of the rotating blade correspond to the imaginary parts of eigenvalues λ i , i = 1, 2, . . . , 2n, i.e., Ω i ≈ Im {λ i }, obtained by the solution of (53). Fig. 7 shows the dependence of the three lowest natural frequencies on the rotor speed which, in this case, is expressed by revolutions per minute (ω = 2πn ω ) in the interval 0-3 000 rpm. Note that the values corresponding to the stationary state (ω = 0) are listed in Table 1. From the curves in Fig. 7, it is very well apparent that the first two natural frequencies, which correspond to bending mode shapes, increase with increasing rotor speed. By contrast, the third natural frequency, representing torsion mode shape, seems to be minimally influenced by the change in rotor speed within the studied rpm interval. It should be pointed out that the knowledge of the aforementioned dependencies in turbo-machines is closely related to design and various operating states of their turbines.

Conclusions
The aim of the present paper was to introduce and verify a novel 1D blade finite element with 16 DOF, whose application in vibration simulations can bring a significant reduction of DOF of the total mathematical model and, at the same time, provide results that are comparable to those obtained from 3D FE models with a much higher number of DOF. The proposed finite element can be implemented into any commercial software, where its connection to other types of elements can be performed by means of transformation matrices between coupled coordinates.
Lastly, it should be pointed out that the blade finite element can serve as a starting point for vibration simulations that consider a turbo-machine as a whole, i.e., by including all the relevant parts such as shaft, generators, blade rings, bandages, bearings and sealing. In this case, one should realize that the total mathematical model will be periodically time-dependent because of stationary parts (bearings and sealing), whose stiffness and damping will vary in the rotating coordinate system during one revolution and be linearly dependent on the angular speed of the turbine. This is also the reason why the Coriolis forces are taken into consideration in the matrix of gyroscopic effects ωG. By contrast, the rotating parts of the turbo-machine (turbine shaft, generator non-symmetrical shaft, blade rings and bandages) will be time-invariable in the rotating coordinate system fixed to the turbo-machine. For response determination and stability assessment, it is possible to use the approach and methodology described by the authors of this paper in [2,3,10].