Vibration analysis of paper machine ’ s asymmetric tube roll supported by spherical roller bearings

This paper presents a simulation method that is used to study subcritical vibrations of a tube roll in a paper machine. This study employs asymmetric 3D beam elements based on the Timoshenko beam theory. An asymmetric beammodel accounts for varying stiffness and mass distributions. Additionally, a detailed rolling element bearing model defines the excitations arising from the set of spherical roller bearings at both ends of the rotor. The results obtained from the simulation model are compared against the results from the measurements. The results indicate that the waviness of the bearing rolling surfaces contributes significantly to the subcritical vibrations while the asymmetric properties of the tube roll have only a fractional effect on the studied vibrations. 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
In manufacturing processes of rotors and bearings, unwanted non-idealities cannot be avoided. As a consequence a rotorbearing system may exhibit undesirable subcritical superharmonic resonances when the rotating speed of the rotor is a fraction of the natural frequency of the rotor [1]. Superharmonic excitations can be caused by geometrical inaccuracies such as roundness errors of the bearing surfaces, or a combination of the gravity and uneven or asymmetric stiffness or mass distributions. In high-end processes such as papermaking, the subcritical vibrations in the process may reduce the quality of the final product or disturb the manufacturing process itself.
The asymmetric rotor has an uneven mass and stiffness distribution resulting in uneven mass moments of inertia with respect to the principal axis. This means that the bending mode pairs of the unsupported shaft have unequal frequencies.
In practice, the cross-section of the asymmetric shaft often varies along the length of the shaft due to manufacturing https://doi.org/10.1016/j.ymssp.2017. 11.030 0888-3270/Ó 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). methods. As for a bearing, manufacturing inaccuracies of shaft geometry, adapter sleeve and inner ring of a bearing assembly contribute mainly to the waviness of the rolling path.
Kang et al. [2] have presented a generalized finite element formulation to model asymmetric rotors. They considered the effect of deviatory inertia and stiffness in an asymmetric shaft with an attached disk on the dynamic behavior of the system. Using the harmonic balance method, they showed that subcritical resonance occurred at speeds equal to the fraction of critical speeds. Ganesan [3] has demonstrated the effect of bearing and shaft asymmetries on the stability of the rotor when the excitation frequency is close to the system's natural frequencies. Two different values for the stiffness of the rotor and bearing in the two principal directions were used. Oncescu et al. [4] have modified the finite element procedure for symmetric shafts to account asymmetric shaft properties. The model is limited due to linearized equations, but it can be used to study weakly non-linear periodic systems.
The dynamics of asymmetric rotors using solid models has been presented by Rao and Sreenivas in [5]. They introduced a transient rotor dynamic analysis for an asymmetric shaft using solid elements that account for stress stiffening and spin softening effects. Lee and Lee [6], and Joh and Lee [7] have proposed a modal testing theory for asymmetric rotor systems. They introduced a unidirectional random excitation technique to estimate directional frequency response functions. Hsieh et al. [8] have studied the combined lateral and torsional vibration in asymmetric rotating systems. They introduced a modified transfer matrix from a symmetric system to an asymmetric one. They learned that synchronous lateral modes were split and new peaks occurred at 2X original frequencies.
Recently, Meng et al. [9] studied the stability of an asymmetric anisotropic system. They used 3D solid elements to define the rotor that enables the more accurate modeling of rotors that have a complex geometry. In comparison with experimental measurement, the reduced linear model is accurate and time efficient and it can be implemented in many realistic case studies. Özs ßahin et al. [10] have introduced an analytical model for asymmetric rotor-bearing systems. As expected, the proposed analytical model is computationally more efficient than the conventional FEM. Frequency response functions were calculated by using the receptance coupling method and modification techniques. The instability of an asymmetric shaft has been studied by Srinath et al. [11]. They proposed an analytical modeling method the stability analysis of a shaft where the equation of motion was formulated on the inertial reference frame.
The effect of the polynomial chaos (PC) expansion order on an asymmetric rotor system response is proposed by Sinou and Jacquelin [12]. The asymmetric properties mainly came from coupling attached to one end of the rotor. By implementing the PC expansion methodology, they studied the first four steady state harmonics, and the effect of PC expansion around the critical speed. They used the Timoshenko beam element with four degrees of freedoms at each node to model the rotor, and they calculated the mean and deviatoric stiffness matrices from the system stiffness matrix.
The objective of this paper is to study subcritical resonances by employing a computationally efficient model. The existing literature shows methods that are used in analyzing asymmetric rotors. The examples are simple and emphasize asymmetric properties of shafts and attached components such as disks. In contrast, this paper introduces an industrial system where the asymmetric properties are fractional, but still might lead to undesired vibrations at the operational speeds. To this end, this study utilizes a finite element method to model an asymmetric paper machine tube roll supported with two non-ideal spherical roller bearings with surface waviness. Waviness in roller bearings usually comes from manufacturing and installation errors, inner and outer race deformation during the operation, or abrasive wear. Several studies are available related to the model of the ball bearing with waviness, such as those by Jang and Jeong [13] and Xu and Li [14]. Recently, Zhang et al. [15] have proposed the stability analysis of a rotor bearing system due to multiple excitations includes bearing waviness. System excitations are considered utilizing the discrete estate transition matrix method (DSTM).
In this study, according to measured data, the wall thickness of the tube roll varies, which is the main cause for asymmetry in the rotor. A generalized equation of motion including the deviatoric inertia and stiffness of an asymmetric rotor is proposed. The waviness of the rolling path of the assembled bearing is measured from the inner rings. Mainly, the emphasis is on the second order waviness (i.e. ellipcity, ovality) that contributes to half of the critical rotational speeds. The third and fourth waviness components are accounted for, but they are expected to have only a minor impact on the vibrations in the studied rotational speed range. All of the measured values are in the range of commonly used manufacturing tolerances. A transient analysis is performed for the system, and the results are converted into the frequency domain, which allows comparison with the experimental responses. The results, displayed as waterfall spectra, indicate that the bearing waviness dominates the superharmonic frequencies, and the uneven mass and stiffness distribution of the rotor play minor roles. However, to study the subcritical responses more accurately, both bearing waviness and rotor asymmetry should be considered.

Mean and deviatoric matrices
The stiffness matrix of an asymmetric beam in a coordinate system that is rotating with the rotor can be formulated as follows ð1Þ a y ¼ 12EIy where E is the elastic modulus of the material, I y and I z are second moments of area with respect to y-and z-coordinates, L is the length of the element, k s is the Timoshenko shear correction factor, G is the shear modulus, and A is the cross-section area. The shear correction factor for a hollow circular cross-section is [17]: where m is Poisson's ratio of the material and R is the ratio between the inner and outer radius of rotor element. If the inner radius equals to zero, Eq. (2) returns the shear correction factor for a solid circular cross-section. In general, rotor-bearing systems can be analyzed either in a rotating or fixed stationary coordinate system. A common choice is to adopt the fixed coordinate system, since the description of bearing and support properties are more straightforward. However, in case of asymmetric rotor, the rotor matrices are not time invariant in the fixed coordinate system. The transformation between the rotating and fixed coordinate systems can be described as follows: where u and u are the displacement vectors in the fixed and rotating coordinate system, respectively. Transformation matrix T 1 can be written as where h is rotor rotation angle. For the asymmetric beam element matrices, the transformation from rotating to fixed coordinate system can be accomplished as follows where the transformation matrix T is constructed from the following block matrices T ¼ where 0 is a 2 Â 2 null-matrix. As can be seen in Eq. (5), rotor stiffness matrix is varying as a function of time as the rotor rotates. It is possible to divide the rotor's stiffness matrix to time invariant and time varying parts as follows where k m is the mean stiffness matrix describing the average stiffness. Correspondingly, k c and k s are the deviatoric stiffness matrices which describe the time varying part of the stiffness. As a conclusion, the asymmetric cross-sectional properties can be accounted for with the help of the mean and deviatoric stiffness and mass matrices, similarly as was shown, for example, in [2,18]. In this study, the mean and deviatory matrices are derived using similar approach as proposed in Kang et al. [2]. However, contrary to their formulation, in this study, the shear deformation is described using the inclusion technique [19], while in formulation proposed by Kang et al. [2] additional two degrees of freedom where added to the beam for the shear deformation description. The mean stiffness matrix can be constructed as follows where With the assumption I z > I y , the deviatoric stiffness matrices k c and k s can be expressed as follows The mean and deviatory mass matrices are defined using the same approach as defining the deviatory stiffness matrices using the assumption I z > I y . Expression of beam elements mass matrices can be found in Appendix A. The area moments of inertia are defined around principal axes instead of mean axes. The angle between the principal and mean axes needs to be accounted for in the description of the equation of motion. Furthermore, it is assumed in this study that the asymmetric shaft properties do not significantly contribute to the gyroscopic matrix, and thus, the mean values of shaft's cross-sectional properties are used in the gyroscopic matrices.

Principal angles and moments of inertias
In asymmetric rotor studies, the principal moments of area and the principal angle are a preferred choice to define the asymmetry. Eqs. (11)-(13) depicts the calculation routine for second moments of area and polar moment of inertia of an arbitrary body about the general coordinate in Cartesian and polar coordinates [20] where I z and I y are the second moment of area about the z and y axes, respectively, and I zy is the polar moment of inertia. The principal moments of inertia and the principal angle can be calculated as [21]: cos 2b À I zy sin 2b ð14Þ tan 2b ¼ ÀI zy The equations above are used in definition of the real constants of the used elements.

Waviness in spherical roller bearing
A spherical roller bearing model has been introduced by Ghalamchi et al. [22]. Principal dimensions which have been used in the simulation are presented in Fig. 2. The bearing forces are calculated as a function of elastic deformation and relative displacement between rollers and inner and outer races. Raceways are assumed to be rigid and only the rollers contact deformations are considered based on Hertzian contact deformation. A single roller force can be calculated as where K c is the roller contact stiffness coefficient and d ij is the elastic deformation of the roller i in row j. Finally, the total bearing forces can be calculated as a sum of all individual roller forces. Fig. 3 shows the waviness pattern on the inner and outer rings. The effect of the waviness profile of the bearing inner ring can be added to the elastic deformation of the roller i in the row j; d ij , as follows where d w ij is the elastic deformation of the roller i in the row j including bearing ring waviness, b ij is the angular position of the roller, h is the rotation angle of the inner ring, k harmonic waviness order, c k is the kth order waviness amplitude, / k is the phase angle of kth order waviness, respectively. The amplitudes and phase angles of the harmonic components can be obtained from measurements by analyzing the results with the FFT.

Governing equation of motion and transient analysis
The equation of motion for an asymmetric rotor is presented as where M m ; M s ; M c are the global mean and deviatoric mass matrices, C is damping matrix, G is the gyroscopic matrix (see Appendix A for mass and gyroscopic matrix extractions), K m ; K c ; K s are the global mean and deviatoric stiffness matrices which are shown in Eqs. (8)-(10), Q is the external force vector that includes the bearing and gravity forces, and Q c and Q s are the force components due to unbalance load.

Verification measurement methods and equipment
Steel tube roll of a paper machine was used for simulation verification. A schematic drawing of the roll with the principal dimensions is presented in Fig. 4.
The rotor was supported by roller element bearings during the measurements. The bearings were spherical roller bearings, SKF 23026CCK/W33. The bearing has 25 rollers in two rows. The bearing was installed on the rotor shaft with a conical adapter sleeve. The parameters of the bearings which was shown in Fig. 2 is presented in Table 1.
During the dynamic measurements, the test rotor was driven with an electric motor. The electric motor was coupled to the rotor through a constant velocity drive shaft to reduce the effect of misalignment.

Bearing measurement set-up
The bearing's inner ring rolling paths were measured to determine waviness components. The excitation frequency caused by the rolling path is the waviness number multiplied by the rotating frequency [23]. The effect of outer ring waviness components was excluded from the study. In this paper, emphasis was on the second order waviness (ovality) that contributes at the half critical rotational speed, but also the third and fourth waviness components were measured and accounted for in the simulation model.
In the measurement procedure, the bearing was first installed on the rotor shaft. Next, the bearing's outer ring and rolling elements were removed, while the inner ring remained fixed to the shaft. This procedure is vital in discovering the inner ring geometry after the installation of the bearing, including the shaft roundness, adapter sleeve thickness variations, inner ring thickness variations and their deformation during the installation process.
The inner ring rolling path geometry was measured with a hybrid-roundness measurement method [24,25], that utilizes four sensor signals to separate the roundness profile and the center point movement of the measured bearing. Four HBM W5TK tactile inductive displacement sensors were used to measure the bearing geometry. The schematic and the actual test set-up of the four-point hybrid method are presented in Fig. 5. Each rolling path was measured in the centre of the path (Fig. 6). A standard roundness measurement machines cannot be used for this measurement, since the inner ring has to be on the rotor shaft during the measurement.  After the bearing roller path measurement, the roundness profile was calculated with hybrid method algorithm [25]. The roundness profile was further analysed with FFT to obtain the amplitudes and phase angles of the waviness components.

Rotor shell thickness measurement set-up
The rotor shell thickness was measured with an ultrasonic measurement device [26] to determine the thickness variation of the tube roll resulting in the asymmetric properties of the rotor. The device consisted of the 5 MHz ultrasonic probe and receiver, Krautkramer H5KF, and the PC-interface card, Krautkramer USPC2100. The thickness of the rotor shell was determined by measuring the time of flight between an emission and the reception of a short ultrasonic pulse, when the velocity of the sound in the material is known [27]. Moreover, water as a sound transmitting medium was used between the probe and the roll surface. The measuring principle (Fig. 7) is based on reflection of the soundwave from all acoustic boundary surfaces.
The roll thickness was measured at 720 equally distributed measurement points with a cylindrical grid of 36 Â 20. In other words, 36 measurement sections in the direction of the roll axis and 20 measurement points circumferentially in each cross-section were used. The roll was measured in a lathe and the ultrasonic probe was fixed on the tool carriage. Rotational position of the roll and longitudinal position of the probe were measured with a rotary encoder. The measurement setup is presented in Fig. 8. The ultrasonic probe was calibrated with a steel step gauge. The uncertainty of the measurement was estimated to be less than 0.1 mm.

Dynamic run-out measurement set-up
The dynamic run-out of the roll was measured with two triangulation reflective laser displacement sensors (Matsushita NAIS LM 300) in the horizontal and vertical directions (see Fig. 4). The measurement system consisted of a PC-based data acquisition unit, two laser sensors and their amplifiers. A rotary encoder was used to trigger the measurement. The displacement was converted to a voltage in laser sensor amplifiers. The amplifier voltage was first low-pass filtered with a cut-off frequency of 500 Hz and then sampled with a sample rate of 2 kHz to avoid aliasing. The cut-off frequency was about five times higher than the highest frequency under interest. The number of samples at each rotor rotating velocity was 16,384 points, which leads to an 8.192  Ultrasonic sensor measurement procedure. The first echo is from the sensor surface, the second from the roll outer surface and the third from the roll inner surface. ToF is the measured soundwave time of the flight. [26]. equipment at a slow reference speed to find out the static run-out of the roll. The static run-out was removed from the actual measurement signals to determine the dynamic response of the roll. An example of unfiltered displacement data in time domain from a tube roll is presented in Fig. 9.
The following signal processing operations were executed in Matlab to produce the frequency spectrums of the measured data: 1. Digitized signals after data acquisition. Static run-out is removed from the signals.

Simulation and measurement results
The rotor thickness was measured from 720 points in total. The thickness varies along the length and the rotation angle of the tube roll as shown in Fig. 10. In the simulation model, the tube roll was divided into 12 beam elements in the axial direction. In each element, 3 Â 20 ¼ 60 points were taken into account to calculate the principal moments of inertia and the rotation angle of the principal axis from the mean axis.
In each angular section of a rotor element, the mean value of the section thickness was calculated from three consecutive points, as illustrated in Fig. 11. Equation of polynomial around the circle in polar coordinate is generated using curve fitting method. Table 2 presents the calculated principal moments of inertia and principal angle for the tube roll elements. The measured rolling path waviness components of the spherical roller bearings are presented in Table 3.
Free-free natural frequencies and modes of the studied tube roll are obtained and compared with measured free-free modes in a conference paper by Ghalamchi et al. [28]. For the completeness of the presented work, the results are reproduced in Table 4. The results show a good agreement and the difference is in the scale of 0.5%, 3%, and 4% for the three first bending mode pairs, respectively.

Dynamic analysis
The obtained peaks from the studies are gathered in Table 5 showing the rotational speeds, corresponding frequencies and amplitudes of 2Â response peak values. Figs. 12 and 13 show the measured and simulated vertical responses, respectively. Since the spectrum lines are obtained at discrete rotation speed steps with a 0.2 Hz increment, the maximum ampli-

x 20
Points for one element Cut 12 elements Fig. 11. Measurement points for one element [28].   The comparison between the measurements and the simulations shows a fair agreement in the horizontal direction. In vertical direction, the simulations clearly differ from the measurements in terms of both the frequency and the amplitude. The large differences in amplitudes are possibly due to the fact that the pedestal and the bearing housing are not measured accurately. Thus, the stiffness and specifically the damping of the support is not emulated in the model as accurately that would be required in order to get the same level of response amplitudes in both measurement and simulations. However, the results indicate that the simulation model captures the responses arising from the studied non-idealities that may cause problems in the operational speed range. The amplitude of the vibrations cannot be repeated, but the simulation model is clearly in resonance at the rotational speed that is half of the first bending mode frequency. As can be seen from the measured results, the real system has excitations that are not considered in the simulation models. The most dominating response type that is missing from the simulations is the 1Â response. Responses that are synchronous with the rotational speed arise mainly from two sources: mass unbalance and the eccentricity of the bearing rings. In this particular case, the tube roll is balanced from five balancing planes and the unbalance is almost negligible in the measurements. Observing the measured 1Â responses confirms this because the responses do not increase considerably as functions of rotational speed even though the forces from the unbalance would be proportional to the rotational speed square.
The complex real life systems has multiple sources for excitations. This paper shows a modeling method for analyzing asymmetric systems. In order to verify the simulation method, Figs. 16 and 17 show the waterfall plots of the tube roll responses in a simulation model where all the other disturbances except the uneven mass and stiffness distribution are neglected. The tube roll experiences excitations twice in a revolution due to a combination of uneven mass distribution and gravity. In addition, the tube roll has lower bending stiffness in the direction of gravity twice in a revolution due to asymmetry. In practice, the study will reveal the share of the subcritical vibration excited by the uneven mass and stiffness distribution. Previous studies in the literature [1,29] have shown that the symmetric model is able to capture subcritical  vibration arising from the bearing waviness. Therefore, it can be assumed that the majority of the responses are from bearing excitations.
In Figs. 16 and 17 the tube roll hits the resonances in both the vertical and the horizontal direction. The frequencies of the highest peaks correspond with the results obtained in the model where the bearing disturbances were included. As expected, the amplitudes of the subcritical vibrations are considerably lower, being 30 and 20 lm in the vertical and horizontal directions, respectively. Still, the responses of the asymmetric tube roll without external disturbances are clearly observable, confirming that using the asymmetric system matrices makes the dynamic model of the rotor more accurate and enables reaching the dynamic responses at half of the rotational speed frequency.

Conclusions
In this study, dynamic responses of an asymmetric rotor structure were analyzed using a computational procedure based on the finite element analysis that employs Timoshenko beam theory. The description of the asymmetric shaft accounts for the uneven mass and stiffness distributions, resulting in a realistic description of the rotor. The example case is a complex industrial component with a number of inaccuracies that are difficult to measure or predict. However, the thickness of the  roll and the bearing rolling path waviness were measured in the laboratory and emulated in the simulation model. The simulation results were presented as frequency spectra that were compared against the measured spectra. The comparison between the simulation and the measurements shows a good agreement, although the simulated responses in the horizontal direction were more accurate. Most likely the differences were due to the properties of the support that were not accurately captured. The simulation method was studied further by removing the bearing waviness from the simulation model when the only excitations in the system were the uneven mass and stiffness distributions. The asymmetric shaft was assumed to cause vibration at the 2Â rotational speed when the gravity excites the system twice per revolution. Nevertheless, the amplitude of the response arising from the asymmetric mass and stiffness was less significant. The responses with the bearing waviness components were tenfold and closer to the measured values. The future studied should include the rotational speeds covering the 3Â and 4Â subcritical speeds in order to evaluate the effects of higher harmonic excitations to the system dynamics.

Acknowledgement
The authors would like to thank Academy of Finland for funding the ViDROM project (Grant No. 277502).  where the unbalance vector for node i can be written as where m u is the unbalance mass, e is unbalance eccentricity, X rotor's rotation speed and a initial angular location of the unbalance mass, respectively.