Modal Analysis of a Thick-Disk Rotor with Interference Fit Using Finite Element Method

This paper is concerned with the modal analysis of a thick-disk rotor, which consists of an elastic shaft with a rigid thick disk assembled by interference fit, and the width of the thick disk is not negligible. Firstly, the friction moment on the contact surface of disk and shaft is deduced in terms of elastic theory, and a new enhanced coefficient of bending stiffness of assembly body is proposed and calculated for the first time. Secondly, the effect of the width of thick disk on diametrical moment of inertia, as well as the enhanced coefficient of bending stiffness of interference-fit part between disk and shaft, is included in the motion equations of thick-disk rotor, which are established based on finite elementmethod, and the natural frequencies of rotor are obtained by solving the motion equations. Then the modal analysis is performed to get the natural frequencies in ANSYS Workbench, in which the friction coefficient and interference fit are set to be the same as those of the finite element calculation method. At last the modal experiment is done to verify the accuracy of calculation and simulation.The results show that the calculation values using enhanced stiffness of assembly part are in good agreement with those of ANSYS Workbench and experiment, and the percent errors of the first natural frequency and the second natural frequency are down to about 0.32% and 0.83%, respectively.


Introduction
Rotor is the key component of rotating machinery, which is widely used in mechanical devices including gas turbines, aeroengines, industrial compressors, and motors etc., and the accurate prediction of natural frequencies of rotor is important for the safe operation of rotating machineries.There are many studies on dynamic analysis of rotor system for the solution of critical speeds and natural frequencies.In the studies by Nelson and Zorzi et al. [1,2], dynamic analysis of rotor system was investigated based on finite element method, and the critical speeds were calculated, which included the effect of various parameters such as rotary inertia, gyroscopic moment, and internal damping.Ku et al. [3] developed a Timoshenko beam finite element model and studied the combined effects of shear deformations and internal damping on forward and backward whirl speeds.Kalita et al. [4] and Chouksey et al. [5] analyzed the effect of fluid film bearing on whirl speeds of rotor-bearing system and found that a half whirl exists in the rotation caused by hydrodynamic bearings.In addition, Lee et al. [6] investigated two different approaches (Floquet theory and coordinate transformation) to develop the complex modal analysis for periodically time-varying linear rotor systems and found that the results of coordinate transformation are much efficient.Khulief et al. [7] investigated the modal characteristics of complex rotor-bearing systems using finite element method and came to the conclusion that the results including the rotational effects were more accurate.Forrai [8] dealt with the stability analysis of self-excited bending vibrations of linear symmetrical rotor-bearing systems with internal damping.
In the establishment of the above finite element model of single-disk rotor or complex rotor system, the width (or thickness) of rotor disk is usually neglected, and the bending stiffness of interference-fit part between disk and shaft equals the stiffness of shaft.In fact, the width of thick disk of rotor has an obvious effect on the diametrical moment of inertia of disk and bending stiffness of interference-fit part, which will change the calculation accuracy of natural frequencies of rotor system.
Liu et al. [9] analyzed the transient vibration of a single thick-disk rotor, whose radial moment of inertia of disk was bigger than the axial one.Pan et al. [10] studied the influence of quality and position of disk on critical speeds of rotor system, which was verified by theoretical calculation and experiment.Petrescu et al. [11] concluded the mass moments of inertia corresponding to different geometric shapes, objects, and profiles.Chen et al. [12] studied the influence of a hot-fit rotor on the local stiffness of the hollow shaft and obtained the accurate natural frequencies using contact elements.Huang et al. [13] optimized the contact stiffness factor of flexible rotor system of magnetic levitation motors and obtained the very close results of natural frequencies by ANSYS Workbench and modal experiment.Kim et al. [14] analyzed the forced vibration of a Timoshenko beam subjected to stationary and moving loads using the modal analysis method.Jalali et al. [15] performed the dynamic analysis of a high speed rotor-bearing system and got the natural frequencies and mode shapes of the rotor at rest under free-free boundary conditions.Giannopoulos et al. [16] analyzed the local flexibility of a shaft with a breathing crack, which is reciprocal of stiffness.In addition, the influences of normal contact stiffness on vibration modes of rod fastening rotors were studied and analyzed in [17,18].
In this study, the modal characteristics of a thick-disk rotor with interference fit is analyzed, whose radial moment of inertia of disk is smaller than the axial one.The friction moment of contact surfaces between disk and shaft is firstly deduced based on elastic theory, and the enhanced coefficient of bending stiffness of thick disk on shaft is proposed and calculated.Then the finite element motion equation of the thick-disk rotor is established, and the enhanced coefficient of stiffness and effect of thick-disk width on diametrical moment of inertia are included in the solution of motion equation for natural frequencies.In addition, the modal analysis is performed on ANSYS Workbench for the verification of calculation results, and the friction coefficient is set to be in accordance with that of finite element calculation.At last the modal experiment is carried out by the application of LMS Test.Lab, and the experiment results are in good agreement with those of theoretical calculation and ANSYS Workbench containing enhanced stiffness of assembly part.The research results not only provide a theoretical basis for the modal analysis of thick-disk rotor with interference fit, but also are suitable for the accurate calculation of natural frequencies of other rotors with interference fit.

Calculation of Friction Moment and Enhanced Coefficient of Bending Stiffness
. .Friction Moment.The mechanical model of disk-shaft assembly body is shown in Figure 1, in which the disk and shaft are assembled by interference fit.When the assembly body is subject to external bending moment, the bending of assembly body will produce friction moment on the contact surfaces because of the relative sliding of disk against shaft.In Figure 1,  is the bending angle of assembly body,  is the external bending moment,  is half of the width of assembly body,  is the internal diameter of hollow disk, which is equal to the diameter of shaft, R 0 is the radius of shaft (d=2R 0 ), and D and R 2 are the external diameter and radius of hollow disk (D=2R 2 ), respectively.
For the convenience of calculation and analysis, the dimensionless parameters are supposed to be as follows: where  is the magnitude of interference fit along diameter direction,  is a coordinate point of assembly body in x-axis direction, and  is the dimensionless parameter corresponding to .The bending stiffness (  ) of assembly body is expressed as where  is the elastic modulus of disk and shaft material,  c is the equivalent section moment of inertia of assembly body,  1 is the section moment of inertia of shaft, and (, , ) is the enhanced coefficient of bending stiffness, which depends on the dimensionless parameters , , .
According to Lame equation [19,20], the interference fit between disk and shaft can be expressed as where  is the compressive stress on the contact surfaces of disk and shaft, which is produced by interference fit.Substituting Eq. ( 1) into Eq.( 3), the compressive stress can be obtained as The pressure of upper and lower contact surfaces between shaft and disk are expressed as where () is the variation of interference fit in radius direction as the assembly body is bent under external moment and   ,   are the pressure of the upper and lower contact surface, respectively.Multiplication coefficient sin  represents   and   which are the contact pressure component of y direction and the contact pressure component of x direction, which can cancel each other out.
Figure 2 shows the friction moment on the contact surfaces.In Figure 2,  0  represents unit area of shaft surface, and the multiplication coefficient sin  represents the unit area perpendicular to the pressures   and   .When the assembly body of disk and shaft is bent under external bending moment, there will be the trend of relative sliding on the contact surfaces, and it is supposed that  is the friction coefficient of contact surfaces, so    0  and    0  are the produced friction forces on the unit area of upper or lower contact surfaces, and they can be expressed as From the integral of Eqs. ( 7) and ( 8), the friction forces can be expressed as follows: So the friction moment can be written Substituting Eqs. ( 9) and ( 10) into Eq.( 11) and solving the integral of Eq. ( 11) along  axis, the friction moment   can be obtained as follows: . .Enhanced Coefficient of Bending Stiffness.In Figure 1, the bending angle of shaft is expressed as The maximum external bending moment of shaft occurs when the stress of cross section reaches yield stress where   is the maximum external bending moment,  is the radius of shaft, and   is the yield stress of shaft material.
Then the external bending moment can be expressed as where the coefficient  satisfies the condition 0<  ⩽1.
To bend the assembly body of disk and shaft, extra external bending moment is needed to overcome the friction moment, so Eq. ( 13) can be expressed as According to Eqs. ( 12), (15), and ( 16), the enhanced coefficient of bending stiffness can be obtained as follows: where the enhanced coefficient depends on the parameters such as friction coefficient, yield stress, Young's modulus, , ,  etc.
Figure 3: Model of single disk rotor system.

Motion Equation of Rotor System
The mathematical model of rotor system is composed of single rigid thick disk, flexible shaft, and bearings, which is illustrated in Figure 3.
. .Rigid ick Disk.In Figure 3, the displacement vector of disk is given by And the motion equation of the rigid disk is expressed as where Ω is the spin angular velocity of disk,   is the vector of generalized forces, and   and   are mass matrix and gyroscopic matrix, respectively, which are defined as where  is the mass of disk and   and   are diametrical and polar moments of inertia of disk, respectively.When the width of thick disk is included in the calculation of moment of inertia as in Figure 4, the diametrical and polar moments of inertia of thick disk can be expressed as follows: [11] where   ,   are the external and inner radius of thick disk, respectively, and  is the width of thick disk. . .Sha Finite Element.The shaft finite element has eight degrees of freedom, and it is defined by two nodes, as shown in Figure 5.Its parameters are denoted as follows:  is the length of the shaft element,  is the area of the cross section of the shaft element,  SD and  SP are diametrical and polar moments of inertia of shaft element, and ,  are the Young's modulus and density of the shaft material, respectively.The nodal displacement vector of the th shaft element is given by When the rotor shaft is discretized into "N" finite elements, the total number of nodes becomes "N+1".And the motion equation of the shaft finite element without internal viscous damping can be expressed as [1,2,5]: where   ,   ,   are the mobile mass matrix, rotating mass matrix, and moment of inertia, respectively, and   is the stiffness of shaft element, which is written as The enhanced stiffness of the assembly part of shaft will be enhanced for the shrink fit of thick disk; according to Eqs.
(2) and ( 17), the enhanced stiffness can be got as follows: . .Bearings.The linear isotropic damped bearings can be modeled by the equation where   is the displacement vector at the bearing location,   is the bearing force vector, and subscript "" represents the bearing.And the damping and stiffness matrices of bearings are defined as where subscript "xx" represents x direction, "yy" represents y direction, and "xy" and "yx" represent the couple direction of x and y.
. .Motion Equation of Rotor System.Based on Eqs. ( 19), (25), and (28), the motion equation of the rotor-bearing system can be written as follows: For the solution of natural frequencies and mode shapes of rotor, the external force and damping can be ignored, and Eq.(30) becomes the free motion equation of thick-disk rotor as follows: where the stiffness matrix  has included the enhanced stiffness.
It is supposed that the solution of Eq. ( 31) is as follows:  =  exp(), and the characteristic equation can be got as follows: In Eq. (32), as not all amplitudes are equal to zero for free vibration rotor, the only condition for nonzero solution is that the coefficient of the determinant is equal to zero.
Then the characteristic values of Eq. ( 33) are the natural frequencies of rotor, and the nonzero vectors of Eq. ( 32) respond to the corresponding vibration modes.

Calculation of Natural Frequencies and Modal Analysis
. .Calculation of Natural Frequencies.The calculation model of a thick-disk rotor is made up of disk and shaft via an interference fit, and the structural parameters of the thickdisk rotor are shown in Table 1.
The finite element structure of the thick-disk rotor is shown in Figure 6.As shown in Figure 6(a), there are 9 nodes and the mass of disk acts on node 5; the stiffness matrix of thick-disk rotor is equal to the stiffness of shaft, so the model solution is not affected by the enhanced coefficient of stiffness, and the results will be the same with those of conventional thin-disk rotor.In Figure 6(b), the rotor is discretized into 11 nodes, the mass of disk acts on node 6, and the enhanced coefficient of stiffness is taken into account in the model establishment, so the stiffness of shaft between nodes 5 and 7 should multiply the enhanced coefficient for the solution of natural frequencies according to Eq. ( 27).
By solving Eq. ( 33), the natural frequencies of thick rotor are obtained using the parameters in Table 1, and the calculation results are shown in Table 2. Comparing the natural frequencies of 9-node rotor with those of 11 nodes, it can be seen that the natural frequencies are greatly influenced by the enhanced coefficient of stiffness, and the percent error of first frequency is about 5.22%, and the second is about 3.15%.

. . Modal Analysis of ick-Disk Rotor in ANSYS Workbench.
The 3D finite element solid model of thick-disk rotor is established in ANSYS Workbench, as shown in Figure 7.The finite element grid of rotor model is automatically generated, and the calculating model contains 2908 nodes and 526 units.By contrast, the natural frequencies obtained in ANSYS Workbench will be more accurate than those by finite element calculation.
The influence of interference fit of thick disk and shaft on stiffness is included in the ANSYS Workbench.So the contact mode is set to be of frictional type in ANSYS Workbench, and the magnitude of interference fit and friction coefficient are set to be in accordance with those of theoretical finite element calculation as in Table 1.In addition, the contact stiffness should be considered in ANSYS Workbench, and the conventional value is from 0.01 to 0.1 for bending analysis.When the contact stiffness is set to be 0.01, solving the model in ANSYS Workbench, the modal shape of the first         simplified to a point, which is applied at the middle point of the rotating shaft, and the mass and moment of inertia of this point are exactly equal to those of the thick disk.Compared to the solid model in Figure 8, the modal solution is not subject to the influence of the enhanced stiffness of assembly part.
Solving the rotor model in ANSYS Workbench, the modal shape of the first two natural frequencies of rotor is obtained as shown in Figure 10.It can be seen that the natural frequencies are 1 =189.03Hz, 2 =482.74Hz,respectively.The results are very close to those of finite element calculation without enhanced stiffness, which verifies that the calculation results of finite element method are consistent with the results . .Modal Experiment and Results Analysis.In order to accurately test the natural frequencies of the thick-disk rotor, the modal experiment with hammering method is done by using LMS Test.Lab.The modal experiment system is shown in Figure 11.The rotor is suspended by a set of bungees to simulate the free-free boundary conditions of the rotor, and the acceleration sensors are fixed on the different positions of the rotor surface (total 30 points).In order to ensure the validity of the amplitude signal and reduce the measurement error, each point is hammered three times, and the average value is adopted to obtain the excitation spectrum of the rotor.
The natural frequencies of modal experiment of the thickdisk rotor are shown in Figure 12, the horizontal coordinate represents the frequency range, the longitudinal coordinate represents the response amplitude, and o, f, v, and d represent the response points of the excitation.In Figure 12, o indicates the mechanical and electrical instability, f indicates the frequency stability, v indicates the pole vector stability, d indicates the damping and frequency stability, and s indicates that the three parameters are stable.And the first two natural frequencies of the rotor obtained by modal experiment are  1 =198.83Hzand  2 =503.58Hz.
The first two natural frequencies of thick-disk rotor with interference fit obtained by different methods are shown in Table 3. Comparing the natural frequencies of different methods, the influence of enhanced stiffness on natural frequencies of thick-disk rotor is very obvious, and the percent errors of natural frequencies with or without enhanced stiffness are about 5.32% and 3.16% in finite element calculation.The results without enhanced stiffness by finite element calculation with 9 nodes are very close to those of rotor with shaft and single mass point in ANSYS Workbench.The percent errors of the first two natural frequencies are down to 0.02% and 0.77%.The results with enhanced stiffness by finite element calculation with 11 nodes are close to those of 3D model in ANSYS Workbench with contact stiffness 0.01.When the contact stiffness is 0.1, the first two natural frequencies are  1 =203.31Hz, 2 =499.39Hz.The results of modal experiment are close to the results with enhanced stiffness.The first two natural frequencies approach the results of 3D solid model in ANSYS Workbench, the contact  stiffness of which is 0.1, 0.01, respectively.And the error percent is down to 0.32% and 0.84%.The variation of contact stiffness is caused by the deformation of rotor under different modes.

Conclusions
In this paper, the enhanced coefficient of stiffness for thickdisk rotor with interference fit is firstly put forward, which is deduced based on elastic theory, and the natural frequencies of thick-disk rotor with enhanced stiffness are calculated using finite element method.In addition, the modal analysis of thick disk is performed by ANSYS Workbench and experiment to verify the accuracy of calculation results, and the following conclusions are obtained.(i) For the thick-disk rotor with interference fit, the width of thick disk has obvious effect on the stiffness of shaft, and the enhanced coefficient of stiffness can be deduced according to elastic theory, which depends on the parameters such as friction coefficient, yield stress, Young's modulus, , , , etc.
(ii) The influence of enhanced stiffness on natural frequencies of thick-disk rotor is very obvious, and the percent errors of natural frequencies with or without enhanced stiffness are about 5.32% and 3.16% in finite element calculation.
(iii) The finite element calculation results without enhanced stiffness are very close to those of model composed of shaft and single mass point in ANSYS Workbench, and the results with enhanced stiffness are consistent with the 3D model in ANSYS Workbench and experiment.

2 MathematicalFigure 1 :
Figure 1: Mechanical model of disk shaft under external bending moment.(a) Axial plane view of model and (b) cross section view of model.

Figure 2 :
Figure 2: Friction moment on contact surfaces.

Figure 9 :
Figure 9: Rotor model with shaft and single mass point.

Figure
Figure B: Modal

Table 1 :
Structure parameters of thick-disk rotor.

Table 2 :
Natural frequencies of thick-disk rotor.

Table 3 :
Natural frequencies of rotor by different methods.