On Using the Transfer Matrix Formulation for Transient Analysis of Nonlinear Rotor Bearing Systems

For many years transfer matrices have been used to evaluate the steady-state vibration response of linear rotor bearing systems. More recently, they have been used to evaluate the steady-state periodic vibration response of nonlinear rotor bearing systems. For quasi-periodic and chaotic response, a transient solution is mandated and transient solution software can also be gainfully used to evaluate the stability of the above-mentioned periodic solutions. To date, transient solutions generally necessitate a different lumped parameter discretization of the rotor and involve solving simultaneously the differential equations for every degree of freedom. This article shows how transient analysis can be performed while maintaining the transfer matrix lumped parameter discretization. The technique is illustrated for a non symmetric unbalanced flexible rotor supported on hydrodynamic journal bearings or deep-groove ball bearings.


INTRODUCTION
Transfer matrix (TM) methods have been widely used to evaluate the vibration behavior of rotating machinery where forces due to nonlinear elements such as hydrodynamic bearings are approximated by stiffness and damping coefficients, thereby implicitly assuming that motions are in close proximity to the steady-state equilibrium position (Lund, 1974).Should motions deviate significantly from this position, e.g., in the presence of relatively large unbalance, but nevertheless achieve a steady state periodic response, it has been possible to use a modified TM ap-proach to determine this response, using harmonic balance principles (Liew et al., 2001;Hahn and Chen, 1994).In such cases, the stability of the equilibrium solutions needs to be addressed, and transient solution capability is a powerful means for doing this.Should the response not be periodic but quasi-periodic or chaotic, a full transient analysis is mandated.Thus, efficient transient analysis capability is necessary for fully investigating the vibrations of nonlinear rotor bearing systems.To date, such analysis generally necessitates a different discretization of the rotor as dictated by finite element formulation (Nelson, 1980), and requires solving simultaneously the differential equations for every degree of freedom, a formidable task which frequently involves system condensation (such as Guyan reduction) to reduce the computational effort.It would therefore be most desirable if transient solutions could be achieved while maintaining the TM model formulation, thereby avoiding the need for model reformulation and the problems associated with system condensation.This article develops such a transient solution technique and evaluates its versatility by analyzing flexible rotor bearing systems with hydrodynamic and rolling element bearings.
While other TM transient solution techniques exist (Gu, 1992;Rao et al., 1987), they are of a more complex nature than the technique developed here and utilize different rotor discretization and/or additional state variables.

TRANSIENT TRANSFER MATRIX TECHNIQUE
The transient TM method developed here matches the usual approach for transient solutions in that ordinary differential equations (ODEs) are formulated and integrated using some convenient technique (e.g., Runge-Kutta).Typically, the n equations of motion, namely; are reformulated to explicitly determine the accelerations, These second-order ODEs are then recast into 2n first-order ODEs, giving: The technique presented here follows the same approach in the formulation of functions f and g.However, rather than simultaneously solving the n equations of motion in Equation (1), transfer matrix relationships are used to solve the four equations of motion at each rotor station.
The TM approach typically discretizes the rotor into massless field and point mass elements (Rao, 1996).Figure 1 shows the forces acting on these elements at the i'th station.Note that such discretization necessarily assumes a lumped mass model.As in the traditional rotordynamic TM formulation, one can write the state vector at the i'th station as: except now the elements of S are time-dependent variables rather than vibration amplitudes.This state vector is transferred through massless spans using a field transfer matrix F and across lumped masses using a point transfer matrix P, namely: and Rao (1996) presents the details of the elements of F and P.

FIGURE 1
Free body diagrams for dynamic equilibrium in the (a) horizontal x-z plane and (b) vertical x-y plane.
To formulate the function f, system accelerations need to be calculated from system displacements and velocities.The first step is to use the displacements and velocities at time t to calculate shear force and bending moments at time t at the lumped masses.This is done by rearranging and partitioning the field transfer equations as follows: where and F is rearranged appropriately.In this formulation, s d are known (arguments of function f ) and s f are the required shear forces and bending moments.Manipulating gives: Note that the boundary conditions of the ends of the rotor are also needed.In the examples provided here, there are lumped masses at the unconstrained ends of the rotor where the shear forces and bending moments are zero.
Once the shear forces and bending moments are explicitly expressed as functions of the displacement state variables, the point transfer matrices are used to determine the system accelerations.Rather than using the full point transfer matrices, only relevant equations are used in order to avoid the problem of having time-dependent point matrices.
Thus the equations of motion for the i'th lumped mass become: 11) consists of second-order ODEs similar to Equation (2) and can be solved in the same manner.

SAMPLE CALCULATIONS
The capabilities of the transient TM technique were evaluated by comparing the transient responses obtained using the proposed TM-based technique with that obtained by existing in-house transient analysis software for the simple two bearing flexible rotor system shown in Figure 2 using both hydrodynamic journal bearings and rolling element bearings.Fourthorder Runge-Kutta integration was used for both the standard and the TM-based solutions.Figure 3 shows the line types used in Figures 5 through 8 and 10 through 13 which compare the TM-based solutions with the standard solutions.

Hydrodynamic Journal Bearings
Figure 4 shows the hydrodynamic journal bearing geometry and the relevant bearing parameters for the two identical bearings.The following assumptions were made.
• External excitation results solely from unbalance and gravity.
• The short bearing approximation (SBA) to Reynolds equation (with constant fluid properties) is valid.
• The pressures at the ends of the journal bearing are atmospheric.
• Only positive pressures contribute to the fluid film forces (i.e., π film or half Sommerfeld condition).
Note that this example is purely for verification of the solution technique.For real situations with larger L/D ratios the SBA may be inappropriate and a more accurate bearing force determination, such as is obtainable using a finite difference solution to Reynolds equation, may need to be used.Referring to Figure 4, are the forces exerted by the bearings on the rotor may then be  written as (Booker, 1965): where: and A typical periodic response is shown in Figures 5 and 6.The system is loaded by an unbalance of 1×10 −4 kgm applied to the central disc and run at 300 rad/s.As can be seen, the transient

FIGURE 6
Journal bearing #2 orbit 300 rad/s.TM technique produces an identical steady-state solution orbit to the standard transient solution.
Figures 7 and 8 show the transient response with zero initial conditions for the first 20 rotor cycles.The rotor speed here is 800 rad/s and no unbalance is applied to the system.Obviously the system is unstable.Again, the transient TM technique produces an identical transient to the standard solution.

Rolling Element Bearings
For this calculation the hydrodynamic journal bearings are replaced by the 7908C ball bearings for which the bearing geometry and relevant parameters are shown in Figure 9.
A fairly simple two-degree of freedom rolling element bearing model is used (Fukata et al., 1985).Referring to Figure 9, the overall contact deformation for the j'th rolling element, δ j , is given by: Summing the contact forces for each rolling element in the y and z directions gives: where FIGURE 9 Rolling bearing geometry.
Figures 10 and 11 show a predominantly periodic solution obtained at a rotor speed of 300 rad/s with an unbalance of 1×10 −4 kgm applied to the central disc while Figures 12 and  13 show the first 20 cycles of a more chaotic solution obtained at 600 rad/s for the same unbalance.Note that if allowed to run further than 20 cycles, the orbit at 600 rad/s produces an annulus shape and the comparison of solutions becomes difficult.Hence only 20 cycles are plotted from an all zero initial condition.
Figure 14 shows the spectra of the results comparing the more chaotic solution at 600 rad/s with the more periodic solution at 300 rad/s.100 cycles were used for the spectra with Hanning windowing.

DISCUSSION
As shown by the numerical examples, the transient TM method produces identical solutions to the standard transient method regardless of the forcing function nonlinearities.There was little difference in the computational time between the proposed and standard transient approaches, with the TM transient approach taking slightly longer.On the other hand, storage requirements are reduced as one is dealing with smaller matrices.The major advantage of the proposed transient TM approach

FIGURE 12
Ball bearing #1 orbit 600 rad/s. is that it complements existing TM techniques using the same lumped mass discretization.
The disadvantage of the proposed transient TM approach is its present restriction to lumped mass discretization and hence a large number of elements may be required for acceptable accuracy.Note, however, that this is not expected to result in numerical difficulties as problem size increases (common to standard

FIGURE 14
Displacement spectra at 300 rad/s and 600 rad/s.TM methods) because there is no continued matrix multiplication.Application to systems with flexible foundations also requires further development.

CONCLUSIONS
A transient TM technique has been successfully developed for nonlinear rotor bearing systems and successfully applied to systems with hydrodynamic journal and rolling element bearings.The major advantage of this technique is its ability to easily existing TM techniques with the same rotor modeling.While there is no significant difference in computation time compared with standard transient analysis software, the advantage of reduced storage requirements inherent to TM approaches has been maintained.Currently the technique is restricted to lumped mass discretization of the rotor and requires further development to include foundation effects.

ACKNOWLEDGMENTS
This work was supported by the Aeronautical and Maritime Research Laboratories DSTO (Australia).

FIGURE 4
FIGURE 4Hydrodynamic journal bearing geometry.
NOMENCLATUREC damping and gyroscopic matrix.[Ns m −1 ] E Y oung's modulus [Pa] f second-order system function F field transfer matrix g first-order system function G shear modulus [Pa] H Forcing function containing gravity, unbalance, and nonlinear bearing forces [N] K stiffness matrix [N m −1 ] M y,z bending moment about y and z axes, respectively [Nm] M mass matrix [kg] n number of system degrees of freedom P y,z applied load (gravity, unbalance, bearing force) in the y and z directions, respectively [N] P point transfer matrix s d displacement elements of S = {−z, θ, y, φ} T s f force elements of S = {M y , −V z , M z , V y } T S state vector = {−z, θ, M y , −V z , y, φ, M z , V y } T t time [s] V y,z shear force in y and z direction, respectively [N] bearing fluid [Pa s] ς angular coordinate along the bearing surface measured from the point of maximum film thickness [rad] ς 1 value of ς at which the pressure region becomes positive [rad] ψ angular position of the line of centers with respect to the y-axis [rad] .,' differentiation with respect to t and τ , respectively Rolling Element Bearing Notation c inner to outer race clearance [m] D b rolling element diameter [m] D p bearing pitch diameter [m] f z,y bearing force in x and y directions, respectively [N] K p load-deflection factor for point contact [Nm −1.5 ] n b number of rolling elements δ contact deformation or deflection [m] φ j angular location of the j'th rolling element [rad] φ 0 initial angular location of the rolling elements [rad] ω c cage angular velocity [rad s −1 ]