Rotor-Dynamic Analysis of a Small Steam Turbine Using Finite Element Method

Rotor-Dynamic Analysis of a Small Steam Turbine Using Finite Element Method

Mehmet ParlakHamdi Taplak 

Mechanical Engineering Department, Nevşehir Hacı Bektaş Veli University, Nevşehir 50300, Turkey

Mechanical Engineering Department, Erciyes University, Kayseri 38280, Turkey

Corresponding Author Email: 
mehmetparlak@nevsehir.edu.tr
Page: 
68-72
|
DOI: 
https://doi.org/10.18280/mmep.070108
Received: 
15 January 2019
|
Revised: 
21 December 2019
|
Accepted: 
28 December 2019
|
Available online: 
31 March 2020
| Citation

© 2020 IIETA. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).

OPEN ACCESS

Abstract: 

Turbine rotors have complex geometry that is difficult to determine their dynamical behaviour. Strong approaches like finite element methods is used to analyze the small steam turbine system to eliminate difficulty of complex geometry problems. Finite element method saves time by allowing reductions in equations beside it is possible to obtain faster solutions by using software that can solve these equations. In this paper, rotor-dynamic analysis of a small steam turbine rotor with certain geometrical and mechanical properties was performed and critical speeds, modal shapes and Campbell diagram of the system was obtained by using FEM based Dynrot software.

Keywords: 

critical speed, finite element method, mode shapes, steam turbine

1. Introduction

Rotating machines are significant sub-systems of engineering systems and extensively used in the industry. So, it is very important issue to understand dynamical behaviour of rotating systems. These machines, especially turbines, must be cautiously designed by taking into consideration their dynamic characteristics. Therefore there are many research that was carried out in the field of rotordynamic analysis.

Sinou [1] investigated the response of a rotor. Non-linear dynamic analysis of a rotor which is supported by roller bearings was performed. He studied on a system comprised of a disk with a single shaft, two flexible bearing supports and a roller bearing. He found that the reason of the exciter is imbalance. He used a numerical method called Harmonic Balance Method for this study.

Das et al. [2] aimed to develop an active vibration control scheme to control the transverse vibrations on the rotor shaft arising from imbalance and they performed an analysis on the vibration control and stability of a rotor-shaft system which has electromagnetic exciters.

Whalley and Abdul-Ameer [3], calculated the rotor resonance, critical speed and rotational frequency of a shaft whose diameter changes by the length. For this purpose, they used basic harmonic response method.

Brusa and Zolfini [4], GGG (Galileo Galilei on the Ground) evaluated and compared the experimental results and the results they obtained by using FEM based Dynrot program in terms of imbalance and stability.

Lei and Palazzolo [5] have analyzed a flexible rotor system supported by active magnetic bearings and synthesized the Campbell diagrams, and eigenvalues to optimize the rotor-dynamic characteristics and obtained the stability at the speed range. They also investigated the rotor critical speed, mode shapes, frequency responses and time responses.

Taplak and Parlak [6] studied on a small gas turbine rotor with certain geometrical and mechanical properties by using FEM based Dynrot program. They obtained the critical speeds, Campbell diagram and the response of the rotor system which has non-homogen compressor.

Jalali et al. [7] carried out full dynamic analysis of a high speed rotor using 3D finite element model, one-dimensional beam type model and experimental modal test. They obtained Campbell diagram, critical speeds, operational deflection shapes, and unbalance response of the rotor. Their finite element model results had good agreement with experimental results.

Chatelet et al. [8] studied on the modeling approaches for dynamic analyses of the rotating assemblies of turbomachines. They compared the results obtained from the 3D finite element modeling of two case studies with those obtained from a software named ROTORINSA which is based on beam type 1D finite elements.

Jeon et al. [9] created the rotordynamic model of a fuel turbopump and performed the dynamic analysis to predict natural frequencies, critical speeds and instability of the system.

Wang et al. [10] modelled a dual-rotor system by using 1D finite element model and 3D finite element model based on actual test rig. The modal frequencies and modal shapes of the dualrotor in static state were calculated and the results were verified by those of modal test.

Qin et al. [11] developed the finite element models to predict dynamic characteristics of rotating bolted disk-drum type rotors.

Wang et al. [12] investigated the comparison of vibration patterns and unbalance responses between rigid shaft with flexible supports like fan rotor system in aero-engine and the flexible rotor with rigid supports by using Euler beam element.

Ma et al. [13] developed a novel rotor-blade model, in which the rotating blades were simplified using Timoshenko beams.

Campbell diagrams and deformations caused by critical speeds is very important in investigation of the dynamic behaviours of rotors. For this reason, obtaining the Campbell diagrams of the rotating systems and determining the critical speeds are very useful for us to observe the system behaviours. There are several programs based on finite element method to obtain these data and diagrams. In this study, a program named Dynrot was used to make dynamic analysis and the evaluation of the results. For this purpose, a small steam turbine rotor certain geometrical and mechanical properties modeled and critical speeds and modal shapes of rotor were calculated by using Dynrot program.

2. Constitutive Relations and Theory (The Finite Element Method)

The finite element method is a general discretization method for the solution of partial derivative differential equations, and consequently, it finds its application in many other fields beyond rotordynamic analysis. The FEM is based on the subdivision of the structure into finite elements, i.e., into parts whose dimensions are not vanishingly small. Many different element formulations have been developed, depending on their shape and characteristics (beam elements, shell elements, plate elements, solid elements, and many others); however, the elements commonly used in elementary rotordynamics are just beam, mass, and spring elements.

Each element is essentially the model of a small deformable solid in which a limited number, usually quite small, of degrees of freedom is substituted for the infinity of degrees of freedom typical of continuous models. Inside each element, the displacement u(x, y, z) of the point of coordinates x, y, z is approximated by the linear combination of a number n of functions, the shape functions, which are assumed arbitrarily,

$\mathbf{u}(\mathrm{x}, \mathrm{y}, \mathrm{z}, \mathrm{t})=\mathbf{N}(\mathrm{x}, \mathrm{y}, \mathrm{z}) \mathbf{q}(\mathrm{t})$  (1)

where, u is the displacement, written as a vector of order 3 in the three-dimensional space (sometimes of higher order, if rotations are also considered), q is a vector in which the n generalized coordinates of the element are listed, and N is the matrix containing the shape functions. There are as many rows in N as in u and as many columns as the number n of degrees of freedom.

Usually the degrees of freedom of the elements are the displacements at given points, which are referred to as nodes. In this case, Eq. (1) is usually reduced to the simpler form

$\left\{\begin{array}{l}U_{x}(x, y, z, t) \\ U_{y}(x, y, z, t \\ U_{z}(x, y, z, t)\end{array}\right\}$$=\left[\begin{array}{ccc}\boldsymbol{N}(x, y, z) & 0 & 0 \\ 0 & \boldsymbol{N}(x, y, z) & 0 \\ 0 & 0 & \boldsymbol{N}(x, y, z)\end{array}\right]\left\{\begin{array}{l}\boldsymbol{q}_{x}(t) \\ \boldsymbol{q}_{y}(t) \\ \boldsymbol{q}_{z}(t)\end{array}\right\}$  (2)

where, the displacements in each direction are functions of the nodal displacements in the same direction only. Matrix N, in this case, has only one row and as many columns as the number of nodes of the element. Eq. (2) has been written for a three-dimensional element; a similar formulation can also be easily obtained for one- or two-dimensional elements.

The shape functions are, as already stated, arbitrary. The freedom in the choice of such functions is, however, limited, because they must satisfy several conditions. A first requirement is a simple mathematical formulation, which is needed to lead to developments that are not too complex. Usually a set of polynomials in the space coordinates is assumed. To get results that are closer to the exact solution of the differential equations, which constitute the continuous model discretized by the FEM, while reducing the size of the elements, the shape functions must

–be continuous and differentiable up to the required order, which depends on the type of element.

–be able to describe rigid-body motions of the element leading to vanishing elastic potential energy.

–lead to a constant strain field when the overall deformation of the element dictates so.

–lead to a deflected shape of each element that matches the shape of the neighboring elements. This means that when the nodes of two neighboring elements displace in a compatible way, all of the interface between the elements must displace in a compatible way.

Another condition, which is not always satisfied, is that the shape functions must be isotropic, i.e., must not show particular geometrical properties that depend on the orientation of the reference frame. Sometimes not all of these conditions are completely met; in particular, there are elements that fail to completely satisfy the matching of the deflected shapes of neighboring elements.

The nodes are usually located at the vertices or on the sides of the elements and are common to two or more of them, but points that are internal to an element can also be used. The equation of motion of each element can be written through Lagrange equation, but there are plenty of alternative formulations, all leading to the same final expressions. The strains can be expressed as functions of the derivatives of the displacements with respect to space coordinates. In general, it is possible to write a relationship of the type

$\mathbf{\epsilon}(\mathrm{x}, \mathrm{y}, \mathrm{z}, \mathrm{t})=\mathbf{B}(\mathrm{x}, \mathrm{y}, \mathrm{z}) \mathbf{q}(\mathrm{t})$  (3)

where, ϵ is a column matrix in which the various elements of the strain tensor are listed (it is commonly referred to as a strain vector, but it is such only in the sense that it is a column matrix) and B is a matrix containing appropriate derivatives of the shape functions with respect to the x, y, z coordinates. B has as many rows as the number of components of the strain vector and as many columns as the number of degrees of freedom of the element.

If the element is free from initial stresses and strains and the behavior of the material is linear, the stresses can be directly expressed from the strains

${\boldsymbol{\sigma}}(\mathrm{x}, \mathrm{y}, \mathrm{z}, \mathrm{t})=\mathbf{E} \boldsymbol{\epsilon}=\mathbf{E}(\mathrm{x}, \mathrm{y}, \mathrm{z}) \mathbf{B}(\mathrm{x}, \mathrm{y}, \mathrm{z}) \mathbf{q}(\mathrm{t})$  (4)

where, E is the stiffness matrix of the material. It is a symmetric square matrix whose elements can theoretically be functions of the space co-ordinates but are usually constant within the element. The potential energy of the element can be easily expressed as

$U=\frac{1}{2} \int_{v}^{0} \epsilon^{T} \boldsymbol{\sigma} d V=\frac{1}{2} \boldsymbol{q}^{T}\left(\int_{v}^{0} \boldsymbol{B}^{T} \boldsymbol{E} \boldsymbol{B} d V\right) \boldsymbol{q}$  (5)

The integral in Eq. (5) is the stiffness matrix of the element

$K=\left(\int_{v}^{0} \boldsymbol{B}^{T} \boldsymbol{E} \boldsymbol{B} d V\right)$  (6)

Because the shape functions do not depend on time, the generalized velocities can be expressed as

$\boldsymbol{u}(x, y, z, t)=\boldsymbol{N}(x, y, z) \dot{\boldsymbol{q}}(t)$  (7)

If all generalized coordinates are related to displacements, the kinetic energy and the mass matrix of the element can be expressed as

$\tau=\frac{1}{2} \dot{\boldsymbol{q}}^{T}\left(\int^{·}_{v} \rho \boldsymbol{N}^{T} \boldsymbol{N} d V\right) \dot{\boldsymbol{q}}$

$\boldsymbol{M}=\left(\int^{·}_{v} \rho \boldsymbol{N}^{T} \boldsymbol{N} d V\right)$  (8)

When some generalized displacements are physically rotations, Eq. (8) must be changed to introduce moments of inertia, but its basic structure remains the same. In the case of a non-natural (gyroscopic) system, the gyroscopic matrix can be obtained together with the mass matrix by taking also the spin speed into account when computing the kinetic energy.

As already stated, the FEM is often used just to compute the stiffness matrix to be used in the context of the lumped-parameters approach. In this case, the consistent mass matrix Eq. (8) is not computed and a diagonal matrix obtained by lumping the mass at the nodes is used. The advantage is that of dealing with a diagonal mass matrix, whose inversion is far simpler than that of the consistent mass matrix. The accuracy is, however, reduced or, better, a greater number of elements is needed to reach the same accuracy; as a consequence, the convenience between the two formulations must be assessed in each case. Generally speaking, the consistent approach leads to values of the natural frequencies that are higher than those computed using the elastic continuum model, whereas those obtained using the lumped-parameters approach are smaller.

If a force distribution p(x, y, z, t) acts on the body, the virtual work dL linked with the virtual displacement δu = Nδq and the nodal force vector f can be expressed in the form

$\delta L=\int_{v}^{0} \delta \boldsymbol{q}^{T} \boldsymbol{N}^{T} \boldsymbol{p}(x, y, z) d V$

$\boldsymbol{f}(t)=\int_{v}^{0} \boldsymbol{N}^{T} \boldsymbol{p}(x, y, z) d V$  (9)

In a similar way, it is possible to obtain the nodal force vectors co-responding to surface force distributions or to concentrated forces acting on any point of the element. The equation of motion of the nonrotating element is then the usual one for discrete undamped nonrotating systems

$M \ddot{q}+K q=f(t)$  (10)

where, vector f contains all forces acting on the element [14].

3. Results

3.1 Small steam turbine system and components

In this section, dimensions and mechanical properties of small steam turbine’s components are presented (Table 1). di and dd are inner and outer diameter of the elements respectivelly. ℓ denotes the axial length of the elements. Small steam turbine consists two sections, compressor and turbine whose mechanical properties and dimensions are same. As seen in Figure 1, finite elemet model of the system was modeled by using 1d beam element mesh in Dynrot software.

Figure 1. Finite element model of the small steam turbine

It can be showed that in Figure 1, the system consist of 7 beam elements that was created by 8 nodes. For instance, element 1 has two nodes, first nodes is beginning point and the second point is end point of the element 1. The system comprise two sections whose lumped mass is on second and seventh nodes and compressor and turbine can be considered as rigid bodies. Parameters for the mass elements are presented in Table 2. In addition bearings were modeled by using spring and damper elements on the fourth and fifth node.

Table 1. Mechanical and geometric properties of the elements

Element Number

1

2

3

4

5

6

7

di (mm)

0

0

0

0

0

0

0

dd (mm)

120

120

63

63

63

120

120

l (mm)

68.3

34.2

100

150

100

34.2

68.3

ρ (kg/m2)

0

0

7810

7810

7810

0

0

E (MN/m2)

2100

 

Table 2. Mass and inertia parameters of the compressors

Mass Number

1

2

Node Number

2

7

m (kg)

20.81

20.81

Jp (kg.m2)

0.285

0.285

Jt (kg.m2)

0.174

0.174

 

3.2 Critical speeds of the system and the system behavior at the critical speeds

Critical speed of a rotating system is the angular velocity that excites the natural frequency in the field of rotordynamics. While angular velocity of the system approaches the system’s critical speed (natural frequency), the system subject to excessive vibration and this cause to resonance. Thus, angular velocity, that cause resonance, is the critical speed of the rotating system.

The critical speed base on stiffness of the shaft and its support, mass of the shafts and attached parts, amount of damping in the system (damping of the shafts and support), magnitude of unbalance mass in the system, the shape of the rotor system.

The main equation (Eq. 11) is used to determine the critical speed for a linear isotropic rotor system, as used in the our model:

$\operatorname{det}\left(-\omega^{2}([M]-[G]-[K \omega])+[K]\right)=0$  (11)

where, [M], [G] and [K] is mass matrix, gyroscopic matrix and stiffness matrix, respectively. ω is also critical speed value of the system.

By solving equation 11, we obtain the critical speed of the small steam turbine respectively as 362,254993 rad/s, 818, 354116 rad/s, 13784,025470 rad/s and 36500,392864 rad/s. Figures 2-9 show the mode shape of the system at these critical speeds respectively.

Figure 2. The mode shape of the rotor at first critical speed

Figure 3. Three-dimensional view of the mode shape at first critical speed

Figure 4. The mode shape of the rotor at second critical speed

Figure 5. Three-dimensional view of the mode shape at second critical speed

Figure 6. The mode shape of the rotor at third critical speed

Figure 7. Three-dimensional view of the mode shape at third critical speed

Figure 8. The mode shape of the rotor at fourth critical speed

Figure 9. Three-dimensional view of the mode shape at fourth critical speed

As seen in the figures related with the mode shapes, the maximum value of the deflection is about 0.05 meters. In addition, because of the symmetrical configuration of the system, mode shapes occur like symmetric. System of Taplak and Parlak’s [6] previous study has a compressor and turbine whose mass and geometric parameters differ from each other, so mode shapes of that system are like non-symmetric.

Another way to identify the critical speed of the system is to use the Campbell diagram which presents the relationship between the natural frequencies and rotational speed. Campbell diagram illustrates the first and second natural frequencies of the turbine as a function of rotor speed [15]. Campbell diagram for the small steam turbine system is shown Figure 10.

Figure 10. Campbell diagram of the rotating system

The first and second critical speed can be found on Campbell diagram (Figure 10). First critical speed was obtained at first forward mode (intersection point of orange curve and black line) and the second critical speed was obtained at the forward mode (intersection point of blue curve and black line).

4. Conclusions

Because of many rotors have complex geometry, using finite element method saves time and provide quite accurate results for rotating machinery. Reasonable results are obtained from the rotor system modeled by FEM based Dynrot software. Results show that the deflection on the bearings area at the first critical speed is relatively small and the deflection on the bearings increase with increment of critical speed. Reason of this result is being low value damping and stiffness coefficient of the bearings. The problem can be solved by optimizing these values.

The studies in this paper helped in understanding the effectiveness of finite element models to predict natural frequencies, modal shapes of a turbine. The results show that finite element analysis for rotating system such as turbine is a powerful and effective tool for rotordynamic analysis.

  References

[1] Sinou, J.J. (2009). Non-linear dynamics and contacts of an unbalanced flexible rotor supported on ball bearings. Mechanism and Machine Theory, 44(9): 1713-1732. https://doi.org/10.1016/j.mechmachtheory.2009.02.004 

[2] Das, A.S., Nighil, M.C., Dutt J.K., Irretier, H. (2008). Vibration control and stability analysis of rotor-shaft system with electromagnetic exciters. Mechanism and Machine Theory, 43(10): 1295-1316. https://doi.org/10.1016/j.mechmachtheory.2007.10.007 

[3] Whalley, R., Abdul-Ameer, A. (2009). Contoured shaft and rotor dynamics. Mechanism and Machine Theory, 44(4): 772-783. https://doi.org/10.1016/j.mechmachtheory.2008.04.010 

[4] Brusa, E., Zolfini, G. (2002). Dynamics of multi-body rotors: Numerical and experimental FEM analysis of the scientific earth experiment Galileo Galilei ground. Meccanica, 37(3): 239-254. https://doi.org/10.1023/A:1020147020815 

[5] Lei, S., Palazzolo, A. (2008). Control of flexible rotor systems with active magnetic bearings. Journal of Sound and Vibration, 314(1–2): 19-38. https://doi.org/10.1016/j.jsv.2007.12.028 

[6] Taplak, H., Parlak, M. (2012). Evaluation of gas turbine rotor dynamic analysis using the finite element method. Measurement, 45(5): 1089-1097. https://doi.org/10.1016/j.measurement.2012.01.032 

[7] Jalali, M.H., Ghayour, M., Ziaei-Rad, S., Shahriari, B. (2014). Dynamic analysis of a high speed rotor-bearing system. Measurement, 53: 1-9. https://doi.org/10.1016/j.measurement.2014.03.010 

[8] Chatelet, E., D’Ambrosio, F., Jacquet-Richardet, G. (2005). Toward global modeling approaches for dynamic analyses of rotating assemblies of turbomachines. Journal of Sound and Vibration, 282: 163-178. https://doi.org/10.1016/j.jsv.2004.02.035 

[9] Jeon, S.M., Kwak, H.D., Yoon, S.H., Kim, J. (2013). Rotordynamic analysis of a high thrust liquid rocket engine fuel (Kerosene) turbopump. Aerosp Sci Technol, 26: 169-175. https://doi.org/10.1016/j.ast.2012.03.005

[10] Wang, N., Jiang, D., Xu, H. (2017). Dynamic characteristics analysis of a dual-rotor system with inter-shaft bearing. Proc IMech Part G: J Aerospace Engineering, 0(0): 1-12. https://doi.org/10.1177/0954410017748969

[11] Qin, Z., Han, Q., Chu, F. (2016). Bolt loosening at rotating joint interface and its influence on rotor dynamics. Engineering Failure Analysis, 59: 456-466. https://doi.org/10.1016/j.engfailanal.2015.11.002

[12] Wang, M., Han, Q., Wen, B., Zhang, H., Guan, T. (2017). Modal characteristics and unbalance responses of fan rotor system with flexible support structures in aero-engine. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 231(9): 1686-1705. https://doi.org/10.1177/0954410016658076

[13] Ma, H., Lu, Y., Wu, Z.Y., Tai, X., Li, H., Wen, B.C. (2015). A new dynamic model of rotor-blade systems. Journal of Sound and Vibration, 357: 168-194. https://doi.org/10.1016/j.jsv.2015.07.036

[14] Genta, G. (2005). Dynamics of Rotating Systems. Mechanical 1st Edition, Springer-Verlag, New York, USA.

[15] Manwell, J.F., McGowan, J.G., Rogers, A.L. (2002). Wind Energy Explained-Theory, Design and Application. John Wiley & Sons Ltd., West Sussex, England.