Normal form analysis of stator rub in rotating machinery

This work considers analysis of sustained impacting cycles of rotating shafts with potentially many disks. The insight that this is an internal resonance phenomena makes this an ideal system to be studied with the method of normal forms. However, the presence of arbitrary non smooth nonlinearities due to impact and rub mean that the method must be extended by incorporating an Alternating Frequency/Time (AFT) step to capture nonlinear forces. The process results in an elegant formulation that can model a very wide variety of rotor systems and is demonstrated by comparing against simulation of a contacting overhung rotor.


Introduction
Contact between rotating machinery and surrounding stators is an issue that can effect a wide variety of engineered systems, from pumps to turbines and drilling rigs. It has also been shown to lead to some highly complex dynamical behaviour [1][2][3] including internal resonance [4]. Recent work by the authors identified that sustained bouncing cycles can be seen as an internal resonance that is identifiable in the rotating coordinate system [5], unlike the phenomena in [4], which are explained in the stationary frame.
The work in [5] offers a useful insight for analysis of such asynchronous bouncing motions, because it identifies that bouncing cycles can primarily be seen as an interaction between just two modes of the system; this suggests that analysis can be performed by reducing the system to just its resonant modes. The method of normal forms [6] embodies this approach, because it works by reducing systems to just the resonant components; all nonresonant components are held within a near-identity transformation. The systems under consideration have the following form in a coordinate system that rotates with the shaft: However after transformation into complex modal rotating system coordinates, and then the normal forms transformation, the system typically becomes represented by just two complex modal amplitudes, which are solved with two transformed harmonic balance equations. Hence this is a very powerful and elegant method to reduce nonlinear problems in rotating shafts that can potentially handle large problems with complex nonlinearities. e-mail: a.d.shaw@swansea.ac.uk

Solution overview
Firstly, system (1) is transformed into first order form: where y = {q,q} T , A contains the underlying linear conservative terms, and therefore all nonconservative and nonlinear terms are contained within n y (y). The system is then transformed to complex modal form using the eigensolutions of A:ṗ = Λp + n p (p) where Λ is a diagonal matrix of eigenvalues and y = (Φp) and Φ contains eigenvalues. Note that this system has the same number of degrees of freedom as system (1) because only one of each conjugate pair of eigensolutions is used. Furthermore, the retained eigenvalues have form ± ω i where it is important that the sign reflects the direction of the physical whirling i.e. positive if it is the same as the shaft rotation.
The method of normal forms seeks a near-identity transformation p = u + h(u) such that the variable u results in a systemu that is much simpler equation than (3) to solve, and can be solved exactly. Note that form of the actual form of the transformation is decided in the frequency domain as a trial solution is implemented. If we subtract (4) from (3) we obtain:ṗ −u =ḣ = Λh + n p (p) − n u (u) A problem occurs in trying to eliminate p from this equation; this has to be done by approximation. Many texts apply the approximation n p (p) ≈ n p (u) to accuracy O(h).
In this method the more accurate approximation The general solution of (4) is assumed to be of form: i.e. a summation of vectors multiplied by exponential time functions, so that each transformed modal variable is complex Fourier series with unknown fundamental frequency ω r and signed harmonics. In the matrix form Ut, U will have dimensions N × n f where n f indicates the length of the Fourier transform, and t will be a vector of all the terms e i ω r t with length n f . The aim is a system where u is simple, and therefore it is desirable that U is as sparse as possible. Similar representations are made for all other variables: Differentiation can also be achieved by noting thaṫ where Ψ is an n f × n f diagonal matrix where the diagonal entries are all of form i ω r . Thus equation (7) may be written: In order to evaluate the trial nonlinear frequency components N p an Alternating Frequency/Time (AFT) is used. This consists of firstly approximating the time series of p at all time steps using p = (U + H −1 ) t (obtained from using the forms of (9) in (6)). Then, all nonlinear and nonconservative forces n p (p) are evaluated based on this time series. Finally a Fast Fourier Transform is used to return the frequency components of the forces N p = F(n p (p)). Note that because our state variables give velocity as well as displacement, there is no need to know ω r to evaluate this stage. Note that U is to be very sparse, whereas H −1 is a matrix of constants from a previous solution or initial guess, and is in general full except where U is nonzero.
The corresponding elements k, of the matrices in (11) can be compared separately because they each relate to a different harmonic of a different modal variable. Hence equation (11) (with t eliminated) can be considered term by term: In general, because we want to simplify equation (4), equation (12) is solved by choosing: However, if Ψ , ≈ Λ k,k , this will cause H k, to be large, violating the assumption of a near identity transformation.
These terms are known as resonant terms and must be solved by choosing: Typically interesting solutions (i.e. not simply synchronous whirling in or out of contact) occur when exactly two modes become internally resonant [5]. Therefore, expressing the harmonic components of the resonant equation of motion (4) will give a harmonic balance problem of the form: Hence our system of N degrees of freedom results in this greatly reduced form. Equation (15) has 5 unknowns; the real and imaginary parts of U i, j and U k, and also the unknown fundamental frequency ω r which is within Ψ j, j and Ψ , (see the definition of Ψ following equation (10)). However, we may impose that one of the transformed modal variables is purely real; this constrains the phase of this modal variable and therefore locks the solution in time. Therefore, with this imposition equation (15) is solvable. The method outlined is very general with regard to the form of nonlinearity in n q (q,q) due to the use of the AFT step. Furthermore, it is also very general with regard to the size of the system matrices, so should scale to rotor systems with many degrees of freedom with appropriate modal truncations.

System description
As shown in Fig. 1 (a), a disc of mass m, polar moment of inertia I p , and diametral moment inertia I d , is mounted on an inertialess rigid shaft of length . The shaft is pinned at point O, and rotations around this point are resisted by a linear isotropic rotational viscously damped spring, with rate k and damping coefficient c. The disc rotates at a constant angular speed Ω about its centre point C; however imperfections in its geometry cause its centre of mass M to be at a distance ε from C, resulting in out of balance forces. This rotor has two degrees of freedom which are the displacements of the centre point u and v as shown in Fig. 1 (b), which are directly coupled to rotations ψ and θ assuming small angles.
A snubber ring with clearance δ exists at a distance a along the shaft. This generates additional restoring forces when the displacement at this point is greater than δ, represented by the nonlinear contact force N x given by where k s is the stator contact stiffness, δ is the clearance of the stator and r = x 2 c + y 2 c . The system is nondimensionalised using the undamped natural frequency of the underlying nonrotating is the effective inertia, and the displacement required to make contact with the stator ∆ = δ a is used to scale displacement: This leads to equations of motion of the form of equation (1) with matrices as follows: where β = k s a k expresses the ratio of stator stiffness to the underlying linear stiffness.

Results
The above analysis was performed for a system with β = 13.2,f = 0.065,Î p = 0.14. The initial value for h −1 is obtained from the linear synchronous whirl solution. The solution was then iterated twice, each time using the current values of U to refine the estimate for h −1 by using equations (13) and (14). At each speed, four initial guesses where chosen -with the 2nd resonant component at initial phase differences of 0, π/2, π and 3π/2, and the amplitude sufficient to give some contact. A 1024 point FFT was used for the AFT stage. Figure 2 presents the underlying linear whirl velocities -although note that these are represented in the rotating coordinate system and hence show a downward trend, due to the subtraction of the drive speed. It can be seen that the linear whirl speeds form a 2:1 ratio at Ω ≈ 3.3; hence this is a good region in which to search for 2:1 internally resonant solutions. Figure 3 shows the results of the analysis, in terms of the amplitudes of the resonant components. At slightly above the drive speed for the onset of internal resonance as found from Fig. 2, the nontrivial solutions begin to be found.
These solution come in two branches, an upper and a lower for each resonant component. Figure 4 allows us to compare the two solutions atΩ = 3.8. It may be seen from Fig. 4 (a) that the upper branch solution is stable -the numerical simulation, initiated with the analytical solution, follows it perfectly (30 orbit cycles are simulated). Figure  4 shows that the lower branch solution is similar in magnitude and character to the upper branch, although the 'loop' in the orbit is oriented differently. However, while the numerical simulation follows the analytical solution initially, but then diverges away from it. This suggests that this solution branch consists of unstable solutions.