# Dynamics of Flexible Rotor Systems with an Interim Mass Unbalanced Disk Using a Spectral Element Model

Authors: Sangkyu Choi; Usik Lee
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A frequency domain spectral element model is developed for a rotor system that consists of two spinning shafts and an interim disk or blade system. In this study, the shafts are represented by spinning Timoshenko beam models, and the interim disk system is represented by a uniform thick rigid disk with an unbalanced mass. In our derivation of the governing equations of motion of the disk system, the disk is considered to be wobbling about the geometric center of the disk at which the spinning shafts are attached. The high accuracy of the proposed spectral element model is evaluated by comparison with the natural frequencies obtained using the conventional finite element method (FEM). The spectral element model is then used to investigate the effects of the unbalanced mass on the natural frequencies and dynamic responses of an example rotor system.


---

## Body

## 1. Introduction

Motors, engines, turbines, and machine tools are typical examples of rotor systems. Rotor systems usually consist of multiple spinning shafts and interim disks (or multiple sets of blades) that are connected to each other to form one-dimensional (1-D) dynamic systems. Such systems are normally supported by multiple bearings. Mass imbalance or the misalignment of interim disks can induce significant vibration that can result in the structural failure of a rotor system. Thus, it is very important to accurately predict the dynamic characteristics of a rotor system in the early design phase. There have been many studies on the modeling and analysis of the dynamics of various rotor systems, and an extensive review on the subject has been reported by Nelson [1].The effects of transverse shear deformation and rotary inertia become more significant when the diameter-to-length ratio of the shaft is large and also when vibration occurs at a high frequency. Thus, Timoshenko beam theory has been widely applied to shafts that spin at very high speeds [2–5], whereas interim disks have been represented by thin or thick rigid disks [6]. Historically, many analytical methods have been applied to the dynamics of rotor systems. These include Rayleigh’s method, Dunkerley’s formula, Holzer’s method, and the transfer matrix method [7]. The transfer matrix method is a very efficient solution method for 1D systems (such as rotor systems) because the sizes of the transfer matrices to be handled in the method are not too large, in general. However, as the transfer matrix method provides dynamic responses only at the end point of a 1D system, postprocessing is normally required to compute dynamic responses at an arbitrary position in the system. The finite element method (FEM) is well recognized as one of the most powerful computational methods to deal with various complex structural engineering problems. Thus, many researchers, including Nelson [3] and Ruhl and Booker [8], have applied the FEM to stability and dynamic analyses of various rotor systems. However, the FEM can require a huge number of degrees of freedom (DOFs) to improve its accuracy, especially when a rotor system is flexible and operating at a very high spinning speed. This can result in a significant increase in computation burden. To address this problem, reduced order modeling techniques have been introduced (e.g., [9]). Although reduced order modeling techniques are effective in reducing the sizes of finite element models, they can degrade the solution accuracy, especially at high frequencies. In the literature, Xiang et al. [10–12] have presented the wavelet-based finite element methods for the analysis of rotor systems. As the alternative to the conventional FEM, the spectral element method (SEM) can be adopted for dynamic modeling and analyses of flexible rotor systems because it meets two requirements: high accuracy up to the frequency range of interest and high computational efficiency by using only a minimum number of DOFs [13, 14]. The SEM might have apparent advantages over other solution methods such as the transfer matrix method and the FEM, especially when it is applied to a 1D dynamic system such as a rotor system.Lee and Jang [5] and Lee [14] were the first to apply the SEM to the dynamic modeling and analysis of rotor systems represented by a spinning flexible Timoshenko shaft, with or without a mass-balanced interim disk. To the authors’ knowledge, the SEM has not been applied to rotor systems that have interim disks with an unbalanced mass.The goals of this study are (1) to develop a spectral element model for rotor systems that consist of multiple flexible shafts and interim disks with an unbalanced mass; (2) to evaluate the proposed spectral element model by comparison with the solutions obtained using the FEM; and (3) to use the proposed spectral element model to investigate the effects of an unbalanced mass on the natural frequencies and dynamic responses of an example rotor system.


## 2. Spectral Element Model for the Disk

Consider a rotor system that consists of two uniform shafts and an interim disk or a set of blades. The rotor system is bearing-supported at both ends of the rotor system, as shown in Figure1. The spectral element model of the rotor system can be derived by assembling the spectral elements for two spinning shafts and one interim disk. The spectral element model for the bearing-supported spinning flexible shafts is available from [14] and is summarized in Appendix A. The spectral element model for the interim disks with an unbalanced mass will be formulated from the equations of motion in this section.Figure 1
Simply supported spinning shaft with an interim thick disk.


### 2.1. Equations of Motion

Figure2 shows a disk system that consists of a uniform disk (thickness 2h and mass mD) and an unbalanced mass (mass mU). The disk is rigid, and the shafts clamped at the geometric centers of the left-hand side and right-hand side surfaces of the disk are always normal to the corresponding surfaces, while the disk is spinning. In Figure 2, (x,y,z) is the global reference coordinate system (denoted by “R”) and (1, 2, 3) is the moving coordinate system fixed on the disk (denoted by “D”). The origins of both coordinate systems are placed at the geometric center “O” of the disk. The unit direction vectors nx, ny, and nz are defined in the x-, y-, and z-directions, respectively, and n1, n2, and n3 are in the 1-, 2-, and 3-directions, respectively. Other unit direction vectors are also defined in Figure 2. The displacements in the y- and z-directions at the geometric center O of the disk are denoted by vo(x,t) and wo(x,t), respectively. The rotations about the x-, y-, and z-axes are represented by θo(x,t), ϕo(x,t), and ψo(x,t), respectively. dU1 and dU2 denote the location of an unbalanced mass mD. The mass center of the whole disk system, including the unbalanced mass, is denoted by “C” and its location is denoted by dC1 and dC2, where dC1=(mU/mD)dC1 and dC2=(mU/mD)dC2. The forces and moments acting on the surfaces of the disk are shown in Figure 3.Figure 2
Deformed configuration of a thick uniform disk with an unbalanced mass (mU).Figure 3
External forces and moments acting on the surfaces of a disk.Assuming that the disk is spinning at a constant speedΩ (radian/s) about its central axis, the kinetic energy of the disk system can be obtained as
(1)T=12mv·v+12ω·I·ω,
where
(2)m=mD+mU,
and v is the velocity of the mass center of the whole disk system, ω is the angular velocity of the disk, and I is the mass moment of inertia matrix of the whole disk system with respect to its mass center. These are given by
(3)v=[u˙o+dC2ϕ˙osin(θo+Ωt)aaa-dC2ψ˙ocos⁡(θo+Ωt)+dC2(θ˙o+Ω)aaa×{ψosin(θo+Ωt)+ϕocos⁡(θo+Ωt)}ϕ˙o]nx+[v˙o+dC1ψ˙o-dC2(θ˙o+Ω)×sin(θo+Ωt)]ny+[w˙o-dC1ϕ˙o+dC2(θ˙o+Ω)cos⁡(θo+Ωt)]nz,ω=(Ω+θ˙o-ϕoψ˙o)n1+[ϕ˙ocos⁡(θo+Ωt)+ψ˙osin(θo+Ωt)]n2-[ϕ˙osin(θo+Ωt)-ψ˙ocos⁡(θo+Ωt)]n3,I=[I11-I120-I12I22000I33],
where
(4)I11=JD+mDdC22+mU(dU2-dC2)2,I12=mDdC1dC2+mU(dU1-dC1)(dU2-dC2),I22=ID+mDdC12+mU(dU1-dC1)2,I33=ID+mD(dC12+dC22)+mU(dU2-dC2)2+mU(dU1-dC1)2,
where ID is the mass moment of inertia and JD is the polar moment of inertia of the uniform disk.The virtual work done by the forces and moments acting on the surfaces of the disk, which are displayed in Figure3, is given by
(5)δW=QxrDδur+QxlDδul+MxrDδθr+MxlDδθl+QyrDδvr+QylDδvl+MzrDδψr+MzlDδψl+QzrDδwr+QzlDδwl+MyrDδϕr+MylDδϕl,
where the subscripts r and l denote the quantities on the right-hand side and left-hand side surfaces of the disk, respectively, and the superscripts D denote the quantities related to the disk. The displacements and rotations at the geometric centers of two opposite surfaces of the rigid disk satisfy the following relationships:
(6)vl=vo-hψl=vo-hψo,vr=vo+hψr=vo+hψo,wl=wo+hϕl=wo+hϕo,wr=wo-hϕr=wo-hϕo,uo=ul=ur,θo=θl=θr,ϕo=ϕl=ϕr,ψo=ψl=ψr.
By using (6), the virtual work described by (5) can be rewritten as
(7)δW=QxDδuo+MxDδθo+QyDδvo+(hQyrD-hQylD+MzD)δψo+QzDδwo-(hQzrD-hQzlD-MyD)δϕo,
where
(8)QjD=QjlD+QjrD,MjD=MjlD+MjrD,(j=x,y,z).By applying the kinetic energy given by (1) and the generalized forces reduced from (7) to Lagrange’s equations [15], followed by a lengthy mathematical manipulation including some approximations, the equations of motion of the disk system can be derived as follows:
(9)mu¨o=QxD,mv¨o+mUdU1ψ¨o=fymU+QyD,mw¨o-mUdU1ϕ¨o=fzmU+QzD,Jθ¨o=MxD,Iψ¨o-JDΩϕ˙o+mUdU1v¨o=dU1fymU+MzD+h(QyrD-QylD),Iϕ¨o+JDΩψ˙o-mUdU1w¨o=-dU1fzmU+MyD-h(QzrD-QzlD).
By using (6), (9) can be rewritten in terms of the displacements and rotations defined at the center of the left-hand surface of the disk as follows:
(10)mu¨l=QxD,m(v¨l+hψ¨l)+mUdU1ψ¨l=fymU+QyD,m(w¨l-hϕ¨l)-mUdU1ϕ¨l=fzmU+QzD,Jθ¨l=MxD,(I+mUdU1h)ψ¨l-JDΩϕ˙l+(mh+mUdU1)(v¨l+hψ¨l)=(h+dU1)fymU+MzD+2hQyrD,(I+mUdU1h)ϕ¨l+JDΩψ˙l-(mh+mUdU1)(w¨l-hϕ¨l)=-(h+dU1)fzmU+MyD-2hQzrD,
where
(11)fymU=mUdU2Ω2cos⁡Ωt,fzmU=mUdU2Ω2sinΩt,I=ID+mUdU12,J=JD+mUdU22.IfdU1=dU2=0, (9) or (10) can be reduced to the equations of motion for the uniform thick disk without an unbalanced mass [14].


### 2.2. Formulation of the Spectral Element Model

By following the general procedure of the spectral element formulation [14], all displacement fields and forcing terms are represented in their spectral forms by using the discrete Fourier transform theory as follows:
(12){ui(t)θi(t)vi(t)ψi(t)wi(t)ϕi(t)}=1N∑n=0N-1{UinDΘinDVinDΨinDWinDΦinD}eiωnt,aaaaaaaaaaaaaaaaaaaaaaaaaaaiaaaaa(i=o,r,l),{QxjD(t)QyjD(t)QzjD(t)MxjD(t)MyjD(t)MzjD(t)}=1N∑n=0N-1{QxjnDQyjnDQzjnDMxjnDMyjnDMzjnD}eiωnt,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa(j=r,l),{fymU(t)fzmU(t)}=1N∑n=0N-1{fynmUfznmU}eiωnt,
where i=-1 is the imaginary unit. The quantities with subscript n are the spectral (Fourier) components of the corresponding time histories, and N/2 is the number of spectral components up to the Nyquist frequency considered in the fast Fourier transform (FFT)-based spectral analysis. The subscripts n used for the spectral components will be omitted hereafter for the sake of brevity.Substituting (12) into (10) gives the equations of motion in the frequency domain as
(13)-ω2mUlD=QxD,-ω2JΘlD=MxD,-ω2m(VlD+hΨlD)-ω2mUdU1ΨlD=fymU+QyD,-ω2(I+mUdU1h)ΨlD-iωJDΩΦlD-ω2(mh+mUdU1)(VlD+hΨlD)=(h+dU1)fymU+MzD+2hQyrD,-ω2m(WlD-hΦlD)+ω2mUdU1ΦlD=fzmU+QzD,-ω2(I+mUdU1h)ΦlD+iωJDΩΨlD+ω2(mh+mUdU1)(WlD-hΦlD)=-(h+dU1)fzmU+MyD-2hQzrD.
Similarly, the following relations can be obtained from (6) by using the spectral forms given by (12) as
(14)VlD-VrD+2hΨlD=0,WlD-WrD-2hΦlD=0,UlD-UrD=0,ΘlD-ΘrD=0,ΦlD-ΦrD=0,ΨlD-ΨrD=0.Equations (13) and (14) can be combined to obtain the matrix form of the equation as follows:
(15)[S11D0S21DS22D]{dlDdrD}={flD0},
where
(16)dlD={UlD,ΘlD,VlD,ΨlD,WlD,ΦlD}T,drD={UrD,ΘrD,VrD,ΨrD,WrD,ΦrD}T,flD={QxD,MxD,QyD+fymU,(h+dU1)fymU+MzD+2hQyrD,aQzD+fzmU,-(h+dU1)fzmU+MyD-2hQzrD2hQyrD,}T,(17)S11D=-[ω2M1000ω2M2iωG0-iωGω2M3],S21D=[I2000H1000H2],S22D=-I6,
where the matrices M1, M2, M3, G, H1, and H2 are defined in Appendix B. In (16), dlD and drD are the spectral components of the nodal degrees of freedom (or, “spectral nodal DOFs”) defined on the left-hand side and right-hand side of the disk, respectively. I2 and I6 are 2×2 and 6×6 identity matrices, respectively, and0 represents null matrices or vectors. In this paper, the T superscript is used to denote the transpose of a vector or matrix quantity. Equation (15) is the spectral element equation that represents the proposed spectral element model for uniform thick disks with an unbalanced mass.


## 2.1. Equations of Motion

Figure2 shows a disk system that consists of a uniform disk (thickness 2h and mass mD) and an unbalanced mass (mass mU). The disk is rigid, and the shafts clamped at the geometric centers of the left-hand side and right-hand side surfaces of the disk are always normal to the corresponding surfaces, while the disk is spinning. In Figure 2, (x,y,z) is the global reference coordinate system (denoted by “R”) and (1, 2, 3) is the moving coordinate system fixed on the disk (denoted by “D”). The origins of both coordinate systems are placed at the geometric center “O” of the disk. The unit direction vectors nx, ny, and nz are defined in the x-, y-, and z-directions, respectively, and n1, n2, and n3 are in the 1-, 2-, and 3-directions, respectively. Other unit direction vectors are also defined in Figure 2. The displacements in the y- and z-directions at the geometric center O of the disk are denoted by vo(x,t) and wo(x,t), respectively. The rotations about the x-, y-, and z-axes are represented by θo(x,t), ϕo(x,t), and ψo(x,t), respectively. dU1 and dU2 denote the location of an unbalanced mass mD. The mass center of the whole disk system, including the unbalanced mass, is denoted by “C” and its location is denoted by dC1 and dC2, where dC1=(mU/mD)dC1 and dC2=(mU/mD)dC2. The forces and moments acting on the surfaces of the disk are shown in Figure 3.Figure 2
Deformed configuration of a thick uniform disk with an unbalanced mass (mU).Figure 3
External forces and moments acting on the surfaces of a disk.Assuming that the disk is spinning at a constant speedΩ (radian/s) about its central axis, the kinetic energy of the disk system can be obtained as
(1)T=12mv·v+12ω·I·ω,
where
(2)m=mD+mU,
and v is the velocity of the mass center of the whole disk system, ω is the angular velocity of the disk, and I is the mass moment of inertia matrix of the whole disk system with respect to its mass center. These are given by
(3)v=[u˙o+dC2ϕ˙osin(θo+Ωt)aaa-dC2ψ˙ocos⁡(θo+Ωt)+dC2(θ˙o+Ω)aaa×{ψosin(θo+Ωt)+ϕocos⁡(θo+Ωt)}ϕ˙o]nx+[v˙o+dC1ψ˙o-dC2(θ˙o+Ω)×sin(θo+Ωt)]ny+[w˙o-dC1ϕ˙o+dC2(θ˙o+Ω)cos⁡(θo+Ωt)]nz,ω=(Ω+θ˙o-ϕoψ˙o)n1+[ϕ˙ocos⁡(θo+Ωt)+ψ˙osin(θo+Ωt)]n2-[ϕ˙osin(θo+Ωt)-ψ˙ocos⁡(θo+Ωt)]n3,I=[I11-I120-I12I22000I33],
where
(4)I11=JD+mDdC22+mU(dU2-dC2)2,I12=mDdC1dC2+mU(dU1-dC1)(dU2-dC2),I22=ID+mDdC12+mU(dU1-dC1)2,I33=ID+mD(dC12+dC22)+mU(dU2-dC2)2+mU(dU1-dC1)2,
where ID is the mass moment of inertia and JD is the polar moment of inertia of the uniform disk.The virtual work done by the forces and moments acting on the surfaces of the disk, which are displayed in Figure3, is given by
(5)δW=QxrDδur+QxlDδul+MxrDδθr+MxlDδθl+QyrDδvr+QylDδvl+MzrDδψr+MzlDδψl+QzrDδwr+QzlDδwl+MyrDδϕr+MylDδϕl,
where the subscripts r and l denote the quantities on the right-hand side and left-hand side surfaces of the disk, respectively, and the superscripts D denote the quantities related to the disk. The displacements and rotations at the geometric centers of two opposite surfaces of the rigid disk satisfy the following relationships:
(6)vl=vo-hψl=vo-hψo,vr=vo+hψr=vo+hψo,wl=wo+hϕl=wo+hϕo,wr=wo-hϕr=wo-hϕo,uo=ul=ur,θo=θl=θr,ϕo=ϕl=ϕr,ψo=ψl=ψr.
By using (6), the virtual work described by (5) can be rewritten as
(7)δW=QxDδuo+MxDδθo+QyDδvo+(hQyrD-hQylD+MzD)δψo+QzDδwo-(hQzrD-hQzlD-MyD)δϕo,
where
(8)QjD=QjlD+QjrD,MjD=MjlD+MjrD,(j=x,y,z).By applying the kinetic energy given by (1) and the generalized forces reduced from (7) to Lagrange’s equations [15], followed by a lengthy mathematical manipulation including some approximations, the equations of motion of the disk system can be derived as follows:
(9)mu¨o=QxD,mv¨o+mUdU1ψ¨o=fymU+QyD,mw¨o-mUdU1ϕ¨o=fzmU+QzD,Jθ¨o=MxD,Iψ¨o-JDΩϕ˙o+mUdU1v¨o=dU1fymU+MzD+h(QyrD-QylD),Iϕ¨o+JDΩψ˙o-mUdU1w¨o=-dU1fzmU+MyD-h(QzrD-QzlD).
By using (6), (9) can be rewritten in terms of the displacements and rotations defined at the center of the left-hand surface of the disk as follows:
(10)mu¨l=QxD,m(v¨l+hψ¨l)+mUdU1ψ¨l=fymU+QyD,m(w¨l-hϕ¨l)-mUdU1ϕ¨l=fzmU+QzD,Jθ¨l=MxD,(I+mUdU1h)ψ¨l-JDΩϕ˙l+(mh+mUdU1)(v¨l+hψ¨l)=(h+dU1)fymU+MzD+2hQyrD,(I+mUdU1h)ϕ¨l+JDΩψ˙l-(mh+mUdU1)(w¨l-hϕ¨l)=-(h+dU1)fzmU+MyD-2hQzrD,
where
(11)fymU=mUdU2Ω2cos⁡Ωt,fzmU=mUdU2Ω2sinΩt,I=ID+mUdU12,J=JD+mUdU22.IfdU1=dU2=0, (9) or (10) can be reduced to the equations of motion for the uniform thick disk without an unbalanced mass [14].


## 2.2. Formulation of the Spectral Element Model

By following the general procedure of the spectral element formulation [14], all displacement fields and forcing terms are represented in their spectral forms by using the discrete Fourier transform theory as follows:
(12){ui(t)θi(t)vi(t)ψi(t)wi(t)ϕi(t)}=1N∑n=0N-1{UinDΘinDVinDΨinDWinDΦinD}eiωnt,aaaaaaaaaaaaaaaaaaaaaaaaaaaiaaaaa(i=o,r,l),{QxjD(t)QyjD(t)QzjD(t)MxjD(t)MyjD(t)MzjD(t)}=1N∑n=0N-1{QxjnDQyjnDQzjnDMxjnDMyjnDMzjnD}eiωnt,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa(j=r,l),{fymU(t)fzmU(t)}=1N∑n=0N-1{fynmUfznmU}eiωnt,
where i=-1 is the imaginary unit. The quantities with subscript n are the spectral (Fourier) components of the corresponding time histories, and N/2 is the number of spectral components up to the Nyquist frequency considered in the fast Fourier transform (FFT)-based spectral analysis. The subscripts n used for the spectral components will be omitted hereafter for the sake of brevity.Substituting (12) into (10) gives the equations of motion in the frequency domain as
(13)-ω2mUlD=QxD,-ω2JΘlD=MxD,-ω2m(VlD+hΨlD)-ω2mUdU1ΨlD=fymU+QyD,-ω2(I+mUdU1h)ΨlD-iωJDΩΦlD-ω2(mh+mUdU1)(VlD+hΨlD)=(h+dU1)fymU+MzD+2hQyrD,-ω2m(WlD-hΦlD)+ω2mUdU1ΦlD=fzmU+QzD,-ω2(I+mUdU1h)ΦlD+iωJDΩΨlD+ω2(mh+mUdU1)(WlD-hΦlD)=-(h+dU1)fzmU+MyD-2hQzrD.
Similarly, the following relations can be obtained from (6) by using the spectral forms given by (12) as
(14)VlD-VrD+2hΨlD=0,WlD-WrD-2hΦlD=0,UlD-UrD=0,ΘlD-ΘrD=0,ΦlD-ΦrD=0,ΨlD-ΨrD=0.Equations (13) and (14) can be combined to obtain the matrix form of the equation as follows:
(15)[S11D0S21DS22D]{dlDdrD}={flD0},
where
(16)dlD={UlD,ΘlD,VlD,ΨlD,WlD,ΦlD}T,drD={UrD,ΘrD,VrD,ΨrD,WrD,ΦrD}T,flD={QxD,MxD,QyD+fymU,(h+dU1)fymU+MzD+2hQyrD,aQzD+fzmU,-(h+dU1)fzmU+MyD-2hQzrD2hQyrD,}T,(17)S11D=-[ω2M1000ω2M2iωG0-iωGω2M3],S21D=[I2000H1000H2],S22D=-I6,
where the matrices M1, M2, M3, G, H1, and H2 are defined in Appendix B. In (16), dlD and drD are the spectral components of the nodal degrees of freedom (or, “spectral nodal DOFs”) defined on the left-hand side and right-hand side of the disk, respectively. I2 and I6 are 2×2 and 6×6 identity matrices, respectively, and0 represents null matrices or vectors. In this paper, the T superscript is used to denote the transpose of a vector or matrix quantity. Equation (15) is the spectral element equation that represents the proposed spectral element model for uniform thick disks with an unbalanced mass.


## 3. Spectral Element Model for the Rotor System

Figure4 shows a bearing-supported rotor system that consists of two uniform spinning shafts and an interim disk with an unbalanced mass. The spectral element model for the rotor system can be derived by assembling spectral elements for two spinning uniform shafts and a disk system. The spectral element model for the disk system is given by (15). The spectral element model for the spinning uniform shafts is provided in Appendix A. The assembly can be derived by using the same procedure that is commonly used in the standard FEM.Figure 4
Spectral element model for a rotor system with an interim disk.The spectral element models for the left shaft (shaft 1), the disk (D), and the right shaft (shaft 2) can be represented by
(18)[S111S121S211S221]{dl1dr1}={fl1fr1}(shaft1),[S11D0S21DS22D]{dlDdrD}={flD0}(diskD),[S112S122S212S222]{dl2dr2}={fl2fr2}(shaft2),
where the subscripts l and r denote the quantities defined at the left-end node and right-end node of a shaft or disk element, respectively.For the assembly of three spectral element equations given by (18), the global spectral nodal DOFs are defined as
(19)dl=dl1={Ul1,Θl1,Vl1,Ψl1,Wl1,Φl1}T,d0=d0D={U0D,Θ0D,V0D,Ψ0D,W0D,Φ0D}T,dr=dr2={Ur2,Θr2,Vr2,Ψr2,Wr2,Φr2}T.
The connectivity conditions at the interfaces between the disk and the left-hand side and right-hand side shafts satisfy the following relations:
(20)dr1=dlD,dl2=drD.
From (6), we can obtain the following relationship in the frequency domain:
(21)dlD=Tl0d0D,drD=Tr0d0D,
where the transformation matrices Tl0 and Tr0 are defined in Appendix C.By using (19) and (20), three spectral element equations can be assembled in the following form:
(22)[S111S121Tl00Tl0TS211Tl0TS221Tl0+Tl0TS11DTl0+Tr0TS112Tr0Tr0TS1220S212Tr0S222]×{dld0dr}={flf0fr},
where
(23)fl=fl1,f0=Tl0Tfr1+fD+Tr0Tfl2,fr=fr2,
with
(24)fD={QxD,MxD,QyD+fymU,dU1fymU+MzD+2hQyrD,aaaiQzD+fzmU,-dU1fzmU+MyD-2hQzrD}T.Dynamic responses in the time domain can be obtained first by solving (22) for the spectral nodal DOFs and then by applying the inverse FFT algorithm to the spectral nodal DOFs to obtain the time histories of dynamic responses.


## 4. Numerical Results and Discussion

### 4.1. Validation of the Proposed Spectral Element Model

To validate the proposed spectral element model, a simply supported rotor system, as shown in Figure1, was considered as a numerical example. The geometric and material properties of the shafts are given by l1=l2=0.25 m, ρA=9.86 kg/m, ρI=9.86×10-4 kgm, ρJ=1.97×10-3 kgm, EI=2.51×104 Nm2, κGA=8.06×107 N, EA=2.51×108 N, and GJ=1.93×104 Nm2. The geometric and material properties of the interim uniform disk are given by h=0.02 m, rD=0.05 m, mD=4.93 kg, ID=3.08×10-3 kgm2, and JD=6.16×10-3 kgm2. The accuracy of the proposed spectral element model was evaluated and the results are shown in Tables 1 and 2 by comparing the natural frequencies obtained by the proposed spectral element model (denoted as SEM) and the conventional finite element model (denoted as FEM). The finite element model used in this study is described in Appendix D.Table 1
Natural frequencies (Hz) of the rotor system (mU = 0 kg, Ω = 100 Hz).
Modes
SEM(3)
FEM(n)
n
=
3
n
=
11
n
=
21
n
=
41
n
=
51
n
=
61
1st
Backward
158.806
158.901
158.807
158.806
158.806
158.806
158.806
Forward
160.245
160.336
160.246
160.245
160.245
160.245
160.245
2nd
Backward
858.810
916.713
859.182
858.891
858.827
858.819
858.810
Forward
891.508
956.441
891.933
891.601
891.528
891.519
891.508
3rd
Backward
1762.83
2390.46
1765.90
1763.51
1762.98
1762.91
1762.83
Forward
1781.60
2416.47
1784.75
1782.30
1781.75
1781.68
1781.60
4th
Backward
2913.20
3581.68
2924.14
2915.61
2913.71
2913.48
2913.20
Forward
3002.85
3676.30
3014.91
3005.51
3003.42
3003.17
3002.85
Note: (n): the total number of finite elements used in the analysis.Table 2
Natural frequencies (Hz) of the rotor system (mU = 0.493 kg, Ω = 100 Hz).
Modes
SEM(3)
FEM(n)
n
=
3
n
=
11
n
=
21
n
=
41
n
=
51
n
=
61
1st
Backward
153.879
153.888
153.880
153.879
153.879
153.879
153.879
Forward
155.229
155.238
155.230
155.230
155.229
155.229
155.229
2nd
Backward
858.096
862.037
858.467
858.177
858.113
858.105
858.096
Forward
890.763
895.183
891.186
890.856
890.783
890.774
890.763
3rd
Backward
1750.56
1779.63
1753.59
1751.22
1750.70
1750.64
1750.56
Forward
1769.11
1798.58
1772.22
1769.80
1769.26
1769.19
1769.11
4th
Backward
2906.57
3040.74
2917.39
2908.95
2907.08
2906.85
2906.57
Forward
2995.77
3143.81
3007.69
2998.40
2996.33
2996.08
2995.77
Note: (n): the total number of finite elements used in the analysis.Table1 shows the natural frequencies of the rotor system when the interim disk is mass-balanced (i.e., mU=0 kg). Table 2 shows the natural frequencies when the interim disk is mass-unbalanced (i.e., mU= 0.493 kg). For the results in Tables 1 and 2, we assumed that the spinning speed of the shaft is Ω= 100 Hz. The SEM results were obtained by using only three finite elements (two finite elements for the shafts and one finite element for the disk). On the other hand, for the FEM results, the number of finite elements was increased gradually until the first six digits of the lowest four natural frequencies had sufficiently converged. In Tables 1 and 2, “backward” and “forward” denote the backward whirling modes and forward whirling modes, respectively. Tables 1 and 2 show that the FEM results converge to the SEM results as the number of finite elements used is increased. Since FEM solutions should converge to exact solutions as the number of finite elements (or DOFs employed in the analysis) is increased to infinity, we conclude that the proposed SEM has the capability to provide very accurate solutions by using only a small number of finite elements in the analysis.


### 4.2. Effects of Unbalanced Mass on Natural Frequencies and Critical Speeds

Figure6 shows the three lowest natural frequencies when the disk system is uniform without an unbalanced mass (i.e., mU=0) and those when the disk system has an unbalanced mass (mU=0.1mD) at two locations specified by (dU1=0, dU2=0.5rD) and (dU1=0.5h, dU2=0.5rD). From Figure 6, the following issues can be investigated: (1) the unbalanced mass tends to lower the natural frequencies. This observation can be also made by comparing the natural frequencies provided in Tables 1 and 2. However, the effect of the unbalanced mass on the second natural frequency is negligible because a nodal point of the second vibration mode happens to coincide with the location of the disk. (2) The effects of an unbalanced mass on the natural frequencies seem to be dominant when it is located apart from the geometric center of the disk in the radial direction rather than in the thickness (or axial) direction.The lowest two critical speeds are also shown in Table3. A critical speed is defined as a spinning speed of the disk that is identical to a natural frequency of the rotor system. Table 3 shows that the existence of an unbalanced mass tends to decrease critical speeds for both backward and forward whirling modes.Table 3
The lowest two critical speeds (Hz) of the rotor system.
Modes
WhenmU = 0 kg
WhenmU = 0.493 kg
1st
Backward
158.481
153.528
Forward
160.717
155.614
2nd
Backward
756.573
755.977
Forward
1052.65
1051.12


### 4.3. Effects of an Unbalanced Mass on Dynamic Responses

Figure7 shows the dynamic responses at the geometric center of the disk (node 0 in Figure 5) when an unbalanced mass (mU=0.1mD) is located apart from the geometric center of the disk only in the radial direction (i.e., dU1=0 and dU2=0.5rD) for two spinning speeds of the shaft: Ω= 100 Hz and Ω= 153.528 Hz (the lowest backward whirling critical speed). Similarly, Figure 8 shows the dynamic responses at the middle of the right-hand side shaft (node 3 in Figure 5). From Figures 7 and 8, the following issues can be investigated: (1) when the unbalanced mass is located apart from the geometric center of the disk only in the radial direction, the rotor system vibrates in forms that are symmetric with respect to the middle of the disk due to the unbalanced mass-induced centrifugal forces. Accordingly, the rotations or wobbling angles ϕo and ψo at the center of the disk are found to be zero, whereas those at the middle of the right-hand side shaft are not zero (for both) when Ω=100 Hz and Ω= 153.528 Hz. (2) Both the transverse vibrations and wobbling angles increase significantly when the shaft is spinning at the critical speed of Ω= 153.528 Hz.Figure 5
Spectral element model for the identification of an unbalanced mass (mU) within an interim disk.Comparison of the lowest three natural frequencies of the rotor systems with and without an unbalanced mass (mU = 0.493 kg) at two different locations.
(a)
1st natural frequency
(b)
2nd natural frequency
(c)
3rd natural frequencyDynamic responses at the geometric center of the disk (node 0) when an unbalanced mass (mU = 0.493 kg) is located at du1=0 and du2=0.5rD.
(a)
ForΩ = 100 Hz
(b)
ForΩ = 153.528 Hz (backward whirling critical speed)Dynamic responses at the middle of the right-hand side shaft (node 3) when an unbalanced mass (mU = 0.493 kg) is located at du1=0 and du2=0.5rD.
(a)
ForΩ = 100 Hz
(b)
ForΩ = 153.528 Hz (backward whirling critical speed)Figures9 and 10 show the dynamic responses at the geometric center of the disk (node 0) and at the middle of the right-hand side shaft (node 3), respectively, when an unbalanced mass (mU=0.1mD) is located apart from the geometric center of the disk both in the radial and thickness (axial) directions (i.e., dU1=0.5h and dU2=0.5rD) for two spinning speeds: Ω=100 and 153.528 Hz. Unbalanced mass-induced moments are generated when the unbalanced mass is moved in the thickness direction by a distance dU1=0.5h. Thus, comparing Figures 7 and 9 shows that the unbalanced mass-induced moments result in nonzero wobbling angles ϕo and ψo at the center of the disk when Ω=100 and 153.528 Hz. Although the wobbling angles ϕo and ψo due to the unbalanced mass-induced moments are found to be insignificant for the present example problem, they can be important with respect to structural failure of bearings in the disk systems.Dynamic responses at the geometric center of the disk (node 0) when an unbalanced mass (mU = 0.493 kg) is located at du1=0.5h and du2=0.5rD.
(a)
ForΩ = 100 Hz
(b)
ForΩ = 153.528 Hz (backward whirling critical speed)Dynamic responses at the middle of the right-hand side shaft (node 3) when an unbalanced mass (mU = 0.493 kg) is located at du1=0.5h and du2=0.5rD.
(a)
ForΩ = 100 Hz
(b)
ForΩ = 153.528 Hz (backward whirling critical speed)Figure11 shows the dynamic response when the shaft is spinning at two critical speeds: the lowest backward whirling critical speed, Ω= 153.528 Hz, and the lowest forward whirling critical speed, Ω= 155.614 Hz. At the lowest backward whirling critical speed of Ω= 153.528 Hz, the shaft is whirling in the direction opposite to that of the spinning shaft (i.e., clockwise), as shown in Figure 11(a). However, Figure 11(b) shows that the shaft is whirling in the same direction as that of the spinning shaft (i.e., counterclockwise) at the lowest forward whirling critical speed of Ω= 155.614 Hz.Dynamic responses at the geometric center of the disk (node 0) when an unbalanced mass (mU = 0.493 kg) is located at du1=0.5h and du2=0.5rD.
(a)
ForΩ = 153.528 Hz  (backward whirling critical speed)
(b)
ForΩ = 155.614 Hz (forward whirling critical speed)


## 4.1. Validation of the Proposed Spectral Element Model

To validate the proposed spectral element model, a simply supported rotor system, as shown in Figure1, was considered as a numerical example. The geometric and material properties of the shafts are given by l1=l2=0.25 m, ρA=9.86 kg/m, ρI=9.86×10-4 kgm, ρJ=1.97×10-3 kgm, EI=2.51×104 Nm2, κGA=8.06×107 N, EA=2.51×108 N, and GJ=1.93×104 Nm2. The geometric and material properties of the interim uniform disk are given by h=0.02 m, rD=0.05 m, mD=4.93 kg, ID=3.08×10-3 kgm2, and JD=6.16×10-3 kgm2. The accuracy of the proposed spectral element model was evaluated and the results are shown in Tables 1 and 2 by comparing the natural frequencies obtained by the proposed spectral element model (denoted as SEM) and the conventional finite element model (denoted as FEM). The finite element model used in this study is described in Appendix D.Table 1
Natural frequencies (Hz) of the rotor system (mU = 0 kg, Ω = 100 Hz).
Modes
SEM(3)
FEM(n)
n
=
3
n
=
11
n
=
21
n
=
41
n
=
51
n
=
61
1st
Backward
158.806
158.901
158.807
158.806
158.806
158.806
158.806
Forward
160.245
160.336
160.246
160.245
160.245
160.245
160.245
2nd
Backward
858.810
916.713
859.182
858.891
858.827
858.819
858.810
Forward
891.508
956.441
891.933
891.601
891.528
891.519
891.508
3rd
Backward
1762.83
2390.46
1765.90
1763.51
1762.98
1762.91
1762.83
Forward
1781.60
2416.47
1784.75
1782.30
1781.75
1781.68
1781.60
4th
Backward
2913.20
3581.68
2924.14
2915.61
2913.71
2913.48
2913.20
Forward
3002.85
3676.30
3014.91
3005.51
3003.42
3003.17
3002.85
Note: (n): the total number of finite elements used in the analysis.Table 2
Natural frequencies (Hz) of the rotor system (mU = 0.493 kg, Ω = 100 Hz).
Modes
SEM(3)
FEM(n)
n
=
3
n
=
11
n
=
21
n
=
41
n
=
51
n
=
61
1st
Backward
153.879
153.888
153.880
153.879
153.879
153.879
153.879
Forward
155.229
155.238
155.230
155.230
155.229
155.229
155.229
2nd
Backward
858.096
862.037
858.467
858.177
858.113
858.105
858.096
Forward
890.763
895.183
891.186
890.856
890.783
890.774
890.763
3rd
Backward
1750.56
1779.63
1753.59
1751.22
1750.70
1750.64
1750.56
Forward
1769.11
1798.58
1772.22
1769.80
1769.26
1769.19
1769.11
4th
Backward
2906.57
3040.74
2917.39
2908.95
2907.08
2906.85
2906.57
Forward
2995.77
3143.81
3007.69
2998.40
2996.33
2996.08
2995.77
Note: (n): the total number of finite elements used in the analysis.Table1 shows the natural frequencies of the rotor system when the interim disk is mass-balanced (i.e., mU=0 kg). Table 2 shows the natural frequencies when the interim disk is mass-unbalanced (i.e., mU= 0.493 kg). For the results in Tables 1 and 2, we assumed that the spinning speed of the shaft is Ω= 100 Hz. The SEM results were obtained by using only three finite elements (two finite elements for the shafts and one finite element for the disk). On the other hand, for the FEM results, the number of finite elements was increased gradually until the first six digits of the lowest four natural frequencies had sufficiently converged. In Tables 1 and 2, “backward” and “forward” denote the backward whirling modes and forward whirling modes, respectively. Tables 1 and 2 show that the FEM results converge to the SEM results as the number of finite elements used is increased. Since FEM solutions should converge to exact solutions as the number of finite elements (or DOFs employed in the analysis) is increased to infinity, we conclude that the proposed SEM has the capability to provide very accurate solutions by using only a small number of finite elements in the analysis.


## 4.2. Effects of Unbalanced Mass on Natural Frequencies and Critical Speeds

Figure6 shows the three lowest natural frequencies when the disk system is uniform without an unbalanced mass (i.e., mU=0) and those when the disk system has an unbalanced mass (mU=0.1mD) at two locations specified by (dU1=0, dU2=0.5rD) and (dU1=0.5h, dU2=0.5rD). From Figure 6, the following issues can be investigated: (1) the unbalanced mass tends to lower the natural frequencies. This observation can be also made by comparing the natural frequencies provided in Tables 1 and 2. However, the effect of the unbalanced mass on the second natural frequency is negligible because a nodal point of the second vibration mode happens to coincide with the location of the disk. (2) The effects of an unbalanced mass on the natural frequencies seem to be dominant when it is located apart from the geometric center of the disk in the radial direction rather than in the thickness (or axial) direction.The lowest two critical speeds are also shown in Table3. A critical speed is defined as a spinning speed of the disk that is identical to a natural frequency of the rotor system. Table 3 shows that the existence of an unbalanced mass tends to decrease critical speeds for both backward and forward whirling modes.Table 3
The lowest two critical speeds (Hz) of the rotor system.
Modes
WhenmU = 0 kg
WhenmU = 0.493 kg
1st
Backward
158.481
153.528
Forward
160.717
155.614
2nd
Backward
756.573
755.977
Forward
1052.65
1051.12


## 4.3. Effects of an Unbalanced Mass on Dynamic Responses

Figure7 shows the dynamic responses at the geometric center of the disk (node 0 in Figure 5) when an unbalanced mass (mU=0.1mD) is located apart from the geometric center of the disk only in the radial direction (i.e., dU1=0 and dU2=0.5rD) for two spinning speeds of the shaft: Ω= 100 Hz and Ω= 153.528 Hz (the lowest backward whirling critical speed). Similarly, Figure 8 shows the dynamic responses at the middle of the right-hand side shaft (node 3 in Figure 5). From Figures 7 and 8, the following issues can be investigated: (1) when the unbalanced mass is located apart from the geometric center of the disk only in the radial direction, the rotor system vibrates in forms that are symmetric with respect to the middle of the disk due to the unbalanced mass-induced centrifugal forces. Accordingly, the rotations or wobbling angles ϕo and ψo at the center of the disk are found to be zero, whereas those at the middle of the right-hand side shaft are not zero (for both) when Ω=100 Hz and Ω= 153.528 Hz. (2) Both the transverse vibrations and wobbling angles increase significantly when the shaft is spinning at the critical speed of Ω= 153.528 Hz.Figure 5
Spectral element model for the identification of an unbalanced mass (mU) within an interim disk.Comparison of the lowest three natural frequencies of the rotor systems with and without an unbalanced mass (mU = 0.493 kg) at two different locations.
(a)
1st natural frequency
(b)
2nd natural frequency
(c)
3rd natural frequencyDynamic responses at the geometric center of the disk (node 0) when an unbalanced mass (mU = 0.493 kg) is located at du1=0 and du2=0.5rD.
(a)
ForΩ = 100 Hz
(b)
ForΩ = 153.528 Hz (backward whirling critical speed)Dynamic responses at the middle of the right-hand side shaft (node 3) when an unbalanced mass (mU = 0.493 kg) is located at du1=0 and du2=0.5rD.
(a)
ForΩ = 100 Hz
(b)
ForΩ = 153.528 Hz (backward whirling critical speed)Figures9 and 10 show the dynamic responses at the geometric center of the disk (node 0) and at the middle of the right-hand side shaft (node 3), respectively, when an unbalanced mass (mU=0.1mD) is located apart from the geometric center of the disk both in the radial and thickness (axial) directions (i.e., dU1=0.5h and dU2=0.5rD) for two spinning speeds: Ω=100 and 153.528 Hz. Unbalanced mass-induced moments are generated when the unbalanced mass is moved in the thickness direction by a distance dU1=0.5h. Thus, comparing Figures 7 and 9 shows that the unbalanced mass-induced moments result in nonzero wobbling angles ϕo and ψo at the center of the disk when Ω=100 and 153.528 Hz. Although the wobbling angles ϕo and ψo due to the unbalanced mass-induced moments are found to be insignificant for the present example problem, they can be important with respect to structural failure of bearings in the disk systems.Dynamic responses at the geometric center of the disk (node 0) when an unbalanced mass (mU = 0.493 kg) is located at du1=0.5h and du2=0.5rD.
(a)
ForΩ = 100 Hz
(b)
ForΩ = 153.528 Hz (backward whirling critical speed)Dynamic responses at the middle of the right-hand side shaft (node 3) when an unbalanced mass (mU = 0.493 kg) is located at du1=0.5h and du2=0.5rD.
(a)
ForΩ = 100 Hz
(b)
ForΩ = 153.528 Hz (backward whirling critical speed)Figure11 shows the dynamic response when the shaft is spinning at two critical speeds: the lowest backward whirling critical speed, Ω= 153.528 Hz, and the lowest forward whirling critical speed, Ω= 155.614 Hz. At the lowest backward whirling critical speed of Ω= 153.528 Hz, the shaft is whirling in the direction opposite to that of the spinning shaft (i.e., clockwise), as shown in Figure 11(a). However, Figure 11(b) shows that the shaft is whirling in the same direction as that of the spinning shaft (i.e., counterclockwise) at the lowest forward whirling critical speed of Ω= 155.614 Hz.Dynamic responses at the geometric center of the disk (node 0) when an unbalanced mass (mU = 0.493 kg) is located at du1=0.5h and du2=0.5rD.
(a)
ForΩ = 153.528 Hz  (backward whirling critical speed)
(b)
ForΩ = 155.614 Hz (forward whirling critical speed)


## 5. Conclusion

The results obtained in this study can be summarized as follows.(1)
A frequency domain spectral element model was developed for a bearing-supported flexible rotor system with a mass-unbalanced disk.(2)
The equations of motion of a disk system were derived by taking into account the wobbling motion of the disk at its geometric center.(3)
Through numerical studies, the following issues were investigated. (a) The proposed spectral element model provides very accurate dynamic characteristics (natural frequencies and critical speeds) by using only a minimum number of finite elements, especially for a high spinning speed. (b) The natural frequencies and critical speeds for the backward whirling modes are, in general, lower than those for the forward whirling modes. (c) The existence of an unbalanced mass tends to decrease critical speeds for both backward and forward whirling modes. (d) The transverse vibrations and wobbling angles were found to increase significantly when the shaft is spinning at a critical speed. (e) Our numerical simulation shows that the unbalanced mass-induced moments can result in wobbling motions (ϕo and ψo) of the disk.(4)
Proposed spectral element model is applicable to the identification of the unbalanced mass parametersmU, dU1, and dU2.



---

*Source: 325434-2014-04-28.xml*