Abstract
In vertically oriented machines with journal bearing, there are no predefined static radial loads, such as dead weight for horizontal rotor. Most of the commercial software is designed to calculate rotordynamic and bearing properties based on machines with a horizontally oriented rotor; that is, the bearing properties are calculated at a static eccentricity. For tilting-pad bearings, there are no existing analytical expressions for bearing parameters and the bearing parameters are dependent on eccentricity and load angle. The objective of this paper is to present a simplified method to perform numerical simulations on vertical rotors including bearing parameters. Instead of recalculating the bearing parameters in each time step polynomials are used to represent the bearing parameters for present eccentricities and load angles. Numerical results are compared with results from tests performed in a test rig. The test rig consists of two guide bearings and a midspan rotor. The guide bearings are 4-pad tilting-pad bearings. Shaft displacement and strains in the bearing bracket are measured to determine the test rig’s properties. The comparison between measurements and simulated results shows small deviations in absolute displacement and load levels, which can be expected due to difficulties in calculating exact bearing parameters.
1. Introduction
Difficulties arise when simulating vertical machines with tiling-pad guide bearings. The cause of this is that, in systems with vertical rotors, there is no predetermined operating point. In vertical machines such as hydropower units and pumps, it is the mass unbalance, the generator’s properties, and the flow conditions in the turbine that determine the position of the rotor [1]. Commercial software uses algorithms designed for calculations of rotordynamic and bearing properties for machines with a horizontally oriented rotor. This means that the shaft has a set operating point in the bearing based on static radial loads caused by the dead weight of the rotor. Bearing properties from software for horizontal machines are calculated at a static eccentricity; these types of software cannot perform simulations of a vertical rotor’s behaviour. Assuming isotropic bearing parameters, that is, setting and , it is only possible to determine half of the natural frequencies and that should thus be avoided [2].
Historically both numerical simulations and experimental studies in test rigs have been performed to evaluate properties of tilting pad bearings [3, 4]. Glienicke et al. [5] performed tests in a horizontal test rig consisting of a rotor and bearings. The design of Glienickes test rig was similar to the test rig presented in this paper, except for the possibility to only operate the test rig in horizontal orientation. It consisted of two 4-pad tilting-pad bearings with a 50 mm inner diameter and midspan rotors of different masses; one of the rotors had a mass of approximately 25 kg.
Childs, Vance, San Andre’s, and Murphy et al. have performed an impressive amount of test rig testing at Texas A&M University. Several studies are based on bearing parameters when the shaft is at a stationary eccentricity. Examples of performed tests at Texas A&M are evaluation of bearing models [6], influence of excitation frequencies on the bearing parameters [7, 8], and difference in bearing parameters depending on whether the load is applied on pad (LOP) or between pad (LBP) [9]. Studies have also been performed when the shaft does not have a stationary eccentricity. San Andrés and De Santiago [10] studied bearing properties for plain radial bearings at high dynamic radial loading, and Wygant et al. [11] determined bearing parameters as a function of preload for tilting pad bearings at dynamic loading. Other universities and institutes have also analysed bearing behaviour in test rigs. De Castro et al. [12] analyzed instabilities in rotor bearing systems by comparing numerical models with a vertical power plant and a horizontal test rig with a midspan rotor and plain journal bearings. For hydrodynamic tilting pad bearings, the stiffness and damping properties are nonlinear [13]. To determine the dynamic properties of tilting pad bearings at large eccentricities, numerical calculations are required [14]; that is, analytical expressions do not exist. To perform simulations of unbalance response in vertical machines it requires that bearing parameters are updated in each time step with respect to eccentricity and load angle, which is time consuming. Cardinali et al. [1] have performed nonlinear simulations of a vertical hydropower unit with tilting-pad bearings and White et al. [15] performed simulations of a vertical pump with tilting-pad bearings.
The objective of this paper is to present a simplified method to perform numerical simulations on vertical rotors including bearing parameters, depending on the load angle, and to compare the numerical results with results from a test rig.
The structure of this paper is first to describe the design of the test rig, which is numerically simulated and used to verify the simulations. Following this the theoretical model regarding bearing representation and numerical models is presented. Subsequent sections present results and conclusions from simulations and measurements.
2. Test Rig Description
The test rig is approximately 1.1 m in height with a width of 0.42 × 0.42 m. It is designed to enable operation in both horizontal and vertical orientations. The test rig consists of two identical guide bearings and a midspan rotor; see Figure 1. The rotor is manufactured with SS-1550 steel and the bearings are 4-pad tilting-pad bearings manufactured by Kingsbury. The properties of the bearings are presented in Table 1. There is a stinger and jaw type coupling between the motor and the rotor. The jaw type coupling is connected between the motor and stinger to even out potential torque fluctuations from the motor and disengage the axial connection between the rotor and motor. The stinger that connects the jaw type coupling to the rotor is axially supported via a bearing which also supports the stinger radially. The stinger is slender with a diameter of 4 mm in order to minimise transference of radial loads from stinger misalignments but still strong enough to axially support the rotor and transfer torque from motor to rotor. The geometry of the midspan rotor and the position of the bearings are presented in later section where the numerical model is presented. The operating speed of the test rig is from 0 to 3000 rpm and the rotor rotates subcritically with respect to lateral bending modes. A detailed description of the test rig is presented in a pervious paper by Nässelqvist et.al. [16].
Radial load and shaft displacement are measured at each bearing using strain gauges (Kyowa KFG-5-350-D16-11L3M2S) and inductive displacement sensors (Contrinex DW-AD-509-M8); see Figure 2. The load-strain relationship was verified using the dead weight of the rotor when the test rig was in horizontal position. In the test rig the excitation of shaft motion is self-induced; that is, the shaft motion is caused by the dead weight of the rotor and applied mass unbalance. Acquisition of data is performed with three HBM Quantum X MX 840A units. The data is sampled at 1200 Hz and forth-order antialiasing filters at 480 Hz are used at each channel.
3. Theoretical Model
The numerical model of the test rig consists of a midspan rotor supported in two tilting-pad bearings and its bearing brackets. The equation of motion for the assembled rotor model, in terms of finite elements, can be written in matrix form at constant rotational speed as where is the mass matrix, is the damping matrix, is the stiffness matrix, is the gyroscopic matrix, is the force vector, is the rotational speed, and is the displacement vector in - and -directions.
The stinger’s influence on the rotor’s dynamic behaviour is insignificant and neglected in the numerical model. The only interconnections that influence the rotor’s behaviour are the radial bearings and bearing brackets. A schematic model of the representation of the test rig used in the analysis is shown in Figure 6.
3.1. Determination of Bearing Parameters
To be able to perform numerical simulations on the test rig, bearing properties for current eccentricities at the operational speed need to be known. The equation of motion for the bearings is where is the bearing’s stiffness matrix, is the bearing’s damping matrix, is the fluid inertia matrix, and are the forces that occur based on the displacements as well as their displacement velocities and accelerations. Figure 3 presents a schematic representation of stiffness, damping, and displacements in a journal bearing. The forces in tilting pad bearings caused by the fluid inertia properties are, in general, significantly lower than the forces caused by the bearing’s damping and stiffness properties and will be disregarded in the remainder of this paper [3].
For tilting-pad bearings, the stiffness and damping are dependent on the load angle; that is, Childs et al. showed a large difference in bearing parameter depending on wether the load is on pad (LOP) or between pad (LBP). A commercial software [17], using Navier-Stokes equation and considering thermal effects on the fluid, is used to calculate the bearing parameters presented in this paper. The bearing properties, for LOP and LBP, as a function of eccentricity, are calculated for the bearings in the test rig. The load is applied in -direction and the shaft eccentricity, , in a radial bearing is the relationship between the radial shaft displacement and the radial bearing clearance. Figure 4 presents calculated and as a function of eccentricity and depending on wether load is on or between pad. In the numerical simulations, 4-order polynomials are used to represent the calculated bearing properties as a function of eccentricity. The value for all polynomials was above 0.98 for eccentricities between 7 and 60%.
(a)
(b)
To perform simulations of the test rig in vertical position, the bearing properties need to be known for all load angles and all occurring eccentricities. To simplify the representation of the bearing parameters harmonic relation is assumed between LOP-LBP and load angle . Calculated stiffness and damping as a function of eccentricity and the assumed harmonic relation between LOP and LBP are used to declare the bearing properties according to the following: where is the load angle. Figure 5 presents as a function of load angle at 40% eccentricity; the stiffness is calculated using (3) and compared with stiffness calculated in the commercial software.
The bearing parameters must be transformed, according to (5), in each time step of the simulations due to the load angle dependency of the bearing parameters. Consider the following: where
3.2. Analysis Method
To simplify the analysis, the system describing the numerical model of the test rig is reduced to a 5-node finite element model [18]. A schematic description of the reduction of the model is presented in Figure 6. The disc of the midspan rotor is represented in the numerical model as a lumped mass, , in node 2. Nodes 1, 2, and 3 connect the beam elements of the rotor with consistent mass matrix. The bearing elements are connected to the rotor in nodes 1 and 3 and to the bearing bracket in nodes 4 and 5. The spring element of the bearing bracket is connected to the ground and the mass of the bearing brackets are modelled as lumped masses in nodes 4 and 5. The diameter of the rotor disc is 168 mm and the radius where the mass unbalance is applied is 70 mm. A commercial finite element software was used to determine mass and stiffness properties of the bearing bracket. The mass of the bearing bracket is 5 kg and it has an isotropic stiffness matrix, , and has no damping; that is, MN/m. The harmonic radial force, , representing the unbalance force is applied in node 2 in the model.
By using a state vector , assembled equation of motion can be written in state space as where
Assuming a homogenous solution in the exponential form and inserting this in (9), an eigenvalue problem can be obtained as where are the eigenvalues and are the corresponding eigenvectors. The obtained eigenvalues and the corresponding eigenvectors can be complex valued and appear in complex conjugate pairs. From the analysis of the eigenvalues, it is possible to study the stability and the eigenfrequencies of the system. The complex eigenvalues can be written as where and
In (12), is the undamped natural frequency of the th mode and is the modal damping ratio associated with the th mode. If the damping ratio, , is less than 1, the damped natural frequencies, , can be obtained from the imaginary part of the eigenvalue according to the following: From the analysis of eigenvalues it is also possible to investigate the stability of the system.
Using the initial conditions together with the sum of the homogenous and particular solution, the total solution of the assembled equation of motion, (1), can be obtained.
The stiffness matrix of the model includes bearing properties as a function of current eccentricity and load angle, (3) and (4). The stiffness and damping matrix of the bearing, and , must be updated in each time step with respect to load angle and shaft eccentricity, which results in a recalculation of in each time step. The angular discretization of the numerical simulations is one degree. The numerical model is implemented in Fortran and (7) is solved using variable step size Runge-Kutta integration the nodal displacements are determined. The displacement in nodes 1 and 3 represents the shaft’s absolute displacement, . Displacements in nodes 4 and 5, , represent compression of bearing bracket and the difference between absolute displacement and housing compression in the shaft’s displacement in the bearing; that is, . Radial load acting on the bearing bracket is calculated by multiplying bearing bracket compression by the bracket stiffness.
4. Results
Natural frequencies and stability for the test rig were calculated from 500 to 3000 rpm. The bearing load and bearing properties are dependent on rotational speed and load angle. Bearing properties used in calculations of natural frequencies and stability are calculated for load on pad at a mass unbalance of 1.7 kgm. Figure 7 presents calculated lateral bending modes and stability for the model presented in Figure 6. The results show that the test rig rotates subcritically in operating speeds of 500–3000 rpm.
(a)
(b)
Measurements and numerical simulations were performed at the operational speed of 2350 rpm. The residual unbalance of the “balanced” rotor is approximately 0.2 × 10−3 kg·m. To visualise the influence of differences in bearing properties due to load angle, three different mass unbalances were used to cause the total unbalances of approximately 1.4 × 10−3, 3.3 × 10−3, and 5.0 × 10−3 kg·m.
4.1. Calculated Bearing Parameters
Bearing parameters were determined for shaft eccentricities of 5–65% at 2350 rpm, for both LOP and LBP. The cross-coupled stiffness and damping are more than a decade lower than the direct stiffness and damping and is therefore neglected in the calculation of bearing load and shaft displacement. Figure 8 presents and as a function of load angle and eccentricity using (3) and (4).
(a)
(b)
4.2. Measured and Simulated Shaft Displacement and Radial Load
For all measurements the same test procedure was used. Oil and bearing temperatures were 20°C when the test was started and a predefined speed ramp file was used to reach 2350 rpm. Shaft displacements and bearing loads were recorded when 2350 rpm was reached.
The bearing loads and displacements presented in this section are measured at the upper bearing and bracket (nodes and ). Using the unbalance weight of 1.4 × 10−3 kg·m, the eccentricity is relatively small and at small eccentricities the difference in bearing parameters is also small, depending on LOP or LBP; see Figure 4. This results in a circular shape of the displacement and load orbit. Both measured and simulated results show the same characteristics (Figures 9 and 10).
(a)
(b)
(a)
(b)
At increased unbalance mass, the shaft eccentricity increases. Increased bearing load also causes increased difference in bearing parameters depending on the load angle. The difference in bearing parameters depending on load angle results in a square-shaped displacement and load orbit; see Figures 11, 12, 13, and 14 where simulation and measurements are performed using the unbalance of 3.3 × 10−3 and 5.0 × 10−3 kg·m. The measured and simulated results show similar characteristics of the orbits in these load cases as well.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
5. Discussion
The rotor in the numerical model consists of only two elements, which should cause a stiffer representation of rotor in numerical simulations than in reality. The current configuration of the test rig is configured to evaluate bearing properties, and a stiff rotor with natural frequencies above operational speed is used. The stiff rotor causes small deflection of the rotor, and the stiffer representation of the rotor in the numerical model will have an insignificant influence on calculated loads and displacements.
Stiffness and damping in tilting-pad bearings are dependent on oil temperature, bearing geometry, radial load, rotational speed, and so forth. determining an exact value of all of these parameters is not possible, which causes a deviation between calculated bearing parameters and real bearing parameters. Still, according to the results presented in this paper, the estimated bearing parameters are good enough to perform numerical simulations of vertical machines and achieve acceptable levels and correct characteristics of loads and displacements.
The bearing parameters calculated in the commercial software are calculated as a function of static eccentricity. The eccentricity occurring in the test rig due to the unbalance whirls around in the bearing. According to prior tests in the test rig, at constant rotational speed, the relation between the amplitudes of applied radial load and shaft eccentricity is independent of whether the radial load is static or dynamic. From these observations, it is assumed that functions presenting bearing parameters calculated at static eccentricity can be used as bearing parameters at a whirling eccentricity by transforming them to a rotating coordinate system (5).
6. Conclusion
The results presented in this paper show that it is possible to numerically simulate the unbalance response of the vertical test rig using the method presented in this paper. Comparison between measured and simulated results show small deviations between measured and simulated displacements and load levels, which can be expected due to difficulties in calculating exact bearing parameters. More importantly, the shape of the orbits corresponds well to both load and displacement at all unbalance levels. The method developed to simulate unbalance response in the test rig brings possibilities to perform unbalance response simulations of vertical machines with tilting pad bearings.
The computational times for the simulations are short because the bearing parameters are presented in polynomial. Calculation of 10,000 time steps (28 revolutions) takes less than 30 seconds on a standard laptop.
Conflict of Interests
This is to state that there is no potential and/or relevant competing financial or other interests between the authors that might be affected by publication of the results contained in this paper. There is therefore no potential conflict of interests to declare.
Acknowledgments
The research presented was carried out as a part of the activities at “Swedish Hydropower Centre (SVC).” SVC was established by the Swedish Energy Agency, Elforsk, and Svenska Kraftnät together with Luleå University of Technology, the Royal Institute of Technology, Chalmers University of Technology, and Uppsala University. See SVC’s website at http://www.svc.nu/. The authors would also like to express their appreciation to Vattenfall Hydropower for its financial support and the staff at Vattenfall R&D who helped them in this project.