Unbalance Responses of Rotor / Stator Systems with Nonlinear Bearings by the Time Finite Element Method

Since the early 1970s, major works in rotordynamics were oriented toward the calculation of critical speeds and unbalance responses. The current trend is to take into account many kinds of non-linearities in order to obtain more realistic predictions. The use of algorithms based on nonlinear methods is therefore needed. This article first describes the time finite element method. The method is then applied to nonlinear rotor/stator systems where bearings present a radial clearance.

is therefore not appropriate since one has to wait until after the transient phase has died which can be very time-consuming.
Thus, frequential methods like the incremental harmonic balance (Kim et al., 1991) or trigonometric collocation methods (Nataraj, 1989) are of great use thanks to their time efficiency.However, these methods become less attractive as the frequency content of the solution includes a high number of harmonics.On the other hand, methods in the time domain, such as the shooting method (Sundararajan and Noah, 1997), can also be used to find out periodic solutions however this one becomes time-consuming when dealing with large systems.
In this article, the time finite element method is described.This time-based method enables us to get steady-state solutions and to assess their stability.To prove its reliability in rotordynamics problems, two examples are addressed.Both consist of rotor/stator systems with radial clearance in the bearings.

Description of the Method
The aim of the time finite element method is to find out the periodic solutions of forced systems.It is based on Hamilton's Law of Varying Action (Bayley, 1975;Baruch and Riff, 1982): The principle of this method is to interpolate the displacement of all spatial degrees of freedom between given instants t i and t i+1 by polynomials (Wang, 1995(Wang, , 1997)).In this article, the Lagrange polynomials of order k are used, however all kinds of polynomials may be used (Park, 1996).The displacement of the 156 D. DEMAILLY ET AL.
jth spatial dof upon [t i , t i+1 ] is then expressed in function of its values at k + 1 time nodes, namely x j (t l i ), l = 0, k.For clarity purposes, instants t 0 i = t i and t k i = t i+1 are called limit nodal times, and instants t l i for l = 1, k − 1 are called internal nodal times.Thus: This expression for displacement [3] can also be used for the virtual displacement by replacing x j by δx j .It can then be derived to obtain ẋ j and δ ẋ j .Including all these approximations in Equation ( 2) with t 0 = t i and t f = t i+1 leads to where N = {N 1 • • • N k+1 }, and ⊗ stands for the right Kronecker product.Here, the vector X i [6] is composed of the values of the displacement of all spatial dof at each nodal time t l i , and the vector B i [7] contains momenta evaluated at limit nodal times only: Since the relation [5] has to be verified for any virtual displacement, the following equation can be obtained: The elementary matrices and vectors of the above expression are: The last task is then to assemble all elementary Equations (8), as done with spatial finite elements (Zienkiewicz and Taylor, 1994), keeping in mind that for a periodic solution x(t 1 ) = x(t n+1 ) and p(t 1 ) = p(t n+1 ).Thus, momentum terms in B i cancel each other, and the final Equation ( 11) can be solved for the nodal displacement vector (Equation ( 12)).

Study of the Stability
The stability of the solution is analyzed by means of the Floquet theory (Nayfeh and Balachandran, 1995).To construct the monodromy matrix, Equation ( 8) is first written in its incremental form by perturbing x(t l i ), p(t i ), and p(t i+1 ), and, then reordered to put limit and internal terms together (Equation ( 12)). with and Since the perturbed momenta are equal to zero for the internal nodal times, i.e., δ B I i = 0, a condensation procedure similar to the Guyan reduction (Cook et al., 1989) is used to obtain: where Equation ( 18) of element 1 is then combined with Equation (18) of element 2 in order to eliminate the common terms δx(t 2 ) and δp(t 2 ).This leads to a new equation involving δx and δp at instants t 1 and t 3 only, which can also be combined with Equation (18) of element 3, and so on.By going on this way until the last element, one finds the final expression (Equation ( 21)) involving the perturbed state at times t 1 = 0 and t n+1 = T .
Equation ( 21) is finally modified through a simple matrix manipulation to get the expression (22) in which the monodromy matrix can be recognized.The eigenvalues of this matrix allow for assessment of the solution's stability and the types of bifurcation encountered (Nayfeh and Balachandran, 1995).

APPLICATIONS
To demonstrate the validity of the time finite element method in nonlinear rotordynamics problems, two examples are addressed.Both include a ball bearing whose modelization is briefly described below (El-Saeidy, 2000).For all the calculations, 8 time finite elements with Lagrange polynomials of order 4 are used.The algorithm includes an arc-length continuation technique with step control (Blair et al., 1997).

Bearing Model
This model uses the Hertz theory to evaluate contact forces between balls and races (see Figure 1).Each ball is located by its angular position: where The relative radial distance d r between the inner and outer races of the bearing in the θ k direction can be expressed in terms of their horizontal and vertical displacements: Thus, based on the Hertz theory, the radial contact force generated by the kth ball on the races is The bearing reaction is obtained by summing all contributions of individual contact forces projected on global axes.
Example 1 This example consists of a 4 dof system (Demailly et al., 2001).The rotor and stator are modeled by a mass and stiffness (see Figure 2).Viscous damping is assumed for both of them.The forces acting on the system are the gravity and an unbalance.The parameters used for this model are given in Table 1.The bearing, composed by eight balls, is modeled as described above.
The calculation of the response of this system to an unbalance has been done for a rotational speed of 0-6000 rpm (100 Hz).Figures 3a and b    Figure 5 shows the stator orbit at 780 rpm (13 Hz), obtained by the time finite element method (continuous line) and the Runge-Kutta integration scheme (dashed line) for the case K H = 10 7 N/m, δ = 60 µm.Both results are in close agreement.It should be noted that the solution achieved by the temporal integrator is still a transient one, even after an elapsed time of 120 periods.This is due to the small amount of damping combined with the multiple losses of contact per period occurring at this frequency.

Example 2
This second example deals with the finite element model of a flexible rotor as represented in Figure 6.This rotor has an overhung disc and is supported by two rolling bearings.Its length is 1.70 m and its diameter is 40 mm.Bearing housings are taken into account through a spring/mass system.Translations in vertical and horizontal directions are considered.Rotational

FIGURE 6
Finite element model of the flexible rotor.

FIGURE 7
Amplitudes of vibration.

FIGURE 9a
Amplitudes of vibration of the bearing housings, 600-6000 rpm.
degrees of freedom at each node of the rotor are also included to take into account the gyroscopic effects.This leads to a total of 32 degrees of freedom.In Figure 7, the amplitudes of vibration at the disc location, due to an unbalance, are shown on the speed range of 600-4200 rpm (10-70 Hz).The main peak of vibration, approximately 1800 rpm (30 Hz), corresponds to the first flexural mode of the vibration of the rotor.It clearly exhibits a hardening spring effect, causing jump phenomena while passing through this critical speed: When slowly accelerating, the vibration suddenly reduces, whereas when decelerating, it abruptly increases.

FIGURE 9b
Amplitudes of vibration of the bearing housings, zoom, 2400-3900 rpm.
Theses jumps, occurring at 1900 and 1780 rpm (31.8 and 29.6 Hz), respectively, are well predicted by the stability analysis.As depicted in Figure 8, one of the Floquet multipliers leaves the unit circle by +1 indicating that cyclic fold bifurcations happen.Secondary Hopf bifurcations have also been detected.
The vibrations at the housings location (see Figure 9) are less important than those at the disc location.In Figure 9, the solid and dashed lines represent the left and right bearing housing, respectively.One can observe a second peak of vibration at 3120 rpm, nearly twice the critical speed, followed by two other minor peaks.As can be seen in Figure 9b, jump phenomena exist for this second peak and appear at 3090 and 3150 rpm
(51.5 and 52.5 Hz).Concerning the right bearing housing, the curve describes a loop meaning that the vibration suddenly reduces while slowly accelerating/decelerating through the rotational speed 3150/3090 rpm.This demonstrates the reliability of the finite element method coupled with an arc-length continuation scheme to calculate without any problem sharp peaks of response.
In Figure 9, the two continuous (dashed) lines, representing the vertical and horizontal displacements of the left (right) bearing housing, are almost superimposed.Nevertheless, it turns out that the orbits of the housings are not really circular as observed in Figures 10 to 13 which present the evolution of their

FIGURE 13
Housings orbits, 4570-4790 rpm.orbit between four frequency ranges.For a better visualization, the scale has been adapted for the orbit of the second housing: a magnification factor of five in Figures 10 and 11 and ten in Figure 12.Eight time finite elements with Lagrange polynomials of order 4 seem to be sufficient to trace out the complex orbits described by the housings, except around 2640 rpm (44 Hz) (cf.subplot 1-1 to 1-3 of Figure 12).

CONCLUSIONS
The present article has briefly presented the time finite element method and demonstrated its applicability in rotordynamics problems involving nonlinear ball bearings with a small amount of clearance.It has been pointed out that this method is perfectly able to track complex solutions as found at secondary peaks of vibration and that the study of the stability of solutions is clearly assessed by building the monodromy matrix with few matrix manipulations.This method does not yield the exact solution, since this one is approached upon each time finite element by polynomials.However, the approximated solution is good enough from an engineering point of view.If an extreme precision is required, either the number of elements or the degree of polynomials can be adjusted in consequence, at the price of an increase of the computational time.
This method is not well suited when the solution contain harmonics in the high frequency range, since an important number of elements should be used to match the solution.However, frequential methods suffer from the same drawback.It should be mentioned that the time finite element method is not restricted to the calculation of forced vibrations such as unbalance responses, but that it can also be applied to oil whirl problems encountered with rotors supported on hydrodynamic journal bearings.

FIGURE 1
FIGURE 1Location of the kth ball by θ k .
FIGURE 3aAmplitudes of vibration versus rotational frequency, clearance value of 20 µm, K H = 10 7 N/m.

FIGURE 5
FIGURE 5Stator orbit at 13 Hz.1200 rpm (20 Hz) exhibits a hardening effect.Again, jump phenomena occur, as predicted by the study of the stability.Curves shown in Figure4are actually those obtained for a Hertz contact stiffness K H of 10 8 N/m, but the same observations are also true for K H = 10 7 N/m.Figure5shows the stator orbit at 780 rpm (13 Hz), obtained by the time finite element method (continuous line) and the Runge-Kutta integration scheme (dashed line) for the case K H = 10 7 N/m, δ = 60 µm.Both results are in close agreement.It should be noted that the solution achieved by the temporal integrator is still a transient one, even after an elapsed time of 120 periods.This is due to the small amount of damping combined with the multiple losses of contact per period occurring at this frequency.