Part I — Theoretical Model for a Synchronous Thermal Instability Operating in Overhung Rotors

A theoretical model has been developed for a synchronous thermal instability that is caused by differential viscous shearing in bearings of overhung rotors. This model employs an unbalance threshold criterion for instability instead of utilizing a traditional frequency-domain stability analysis. The current model will be used to investigate several case studies for both plain and tilting pad journal bearing rotors in the second part of this article.


INTRODUCTION
Within the last decade much interest has been shown in the area of thermal bending in overhung rotors.Recent theoretical (Keogh and Morton, 1994) and experimental (deJongh and Morton, 1994;Faulkner et al., 1997) studies have focused on this problem.So far, most of the theoretical work has been done using complex analysis, and almost no comparison between experimental and theoretical research has been made.This article develops a simple model of rotor thermal bending without using complex analysis.
Rotor thermal bending is a complicated phenomenon which can lead to unstable spiral vibrations (Schmied, 1987).This phenomenon is primarily due to a temperature difference developing across the journal.The temperature difference can either be caused by the journal rubbing against stationary components or by viscous shearing within the lubricant.The former mechanism was first noticed by Newkirk in the 1920s and later analyzed by Dimarogonas (1973).This mechanism is called the Newkirk Effect.The latter mechanism has recently been studied (Keogh and Morton, 1994) and is often referred to as the Morton Effect.This article concentrates on the Morton Effect.
The Morton Effect occurs when the journal is executing a synchronous orbit around the bearing center.This centered orbit causes one portion of the journal surface to always be at the minimum film thickness, while a diametrically opposite section of the journal surface is always at the maximum film thickness.Lower film thickness areas are generally associated with higher viscous shear stresses which produce higher temperatures.As a result, a hot spot will develop on the journal surface exposed to the minimum film thickness region and a cold spot will be formed on the surface at maximum film thickness.This leads to a temperature gradient developing across the journal.Such a gradient creates a thermal unbalance which is in the direction of the temperature gradient.Thermal bending will then occur if the resultant of the thermal unbalance and the overhung mechanical unbalance is adequate.Under these conditions the bent shaft will decrease the bearing clearance and elevate the thermal gradient.The increased temperature gradient will then initiate more thermal bending.These actions describe a positive feedback mechanism which will drive the system unstable.
Most of the prior models of the Morton Effect have involved complex stability analyses (Keogh and Morton, 1994;Larsson, 1999aLarsson, , 1999b)), and the typical outline of these models is shown Figure 1.
In these models, an initial deformed rotor configuration is assumed and then used as an input in a dynamic model of the rotor system.This dynamic model produces a synchronous orbit which causes a circumferential temperature distribution to develop on portions of the shaft within the bearings.A thermal bend, which depends on the temperature distribution and the rotor thermal characteristics, is subsequently initiated.This thermal bend affects the rotor dynamics and gives rise to the The initial mechanical and thermal bends, which are used in the models by other researchers, produce unbalances that affect the response of the rotor system and determine the size of the synchronous orbit.Hence, the bends in the previous models can be replaced with equivalent unbalances.
In real machines satisfactory performance is established by using specified vibration, or unbalance, threshold levels that are given in standards published by various organizations.One method of specifying such levels is to restrict the maximum allowable centrifugal force that acts on the machine.This restriction then limits the allowable unbalance for each speed value.Such a criterion determines whether the rotor will be "stable" or "unstable" during operation, i.e., if the unbalance level exceeds the threshold level given by the allowable unbalance then the rotor will be "unstable" and will have to be re-balanced.This notion of practical stability will be utilized in the current model.
The current model for the Morton Effect is shown in Figure 2 and first requires an estimate for the initial mechanical unbalance in the rotor.This unbalance is then put into the dynamic rotor model to obtain the synchronous orbit.A fluid film bearing model is then used to establish the temperature distribution which leads to the thermal unbalance.Finally, an unbalance threshold criterion is used to determine the system stability.

INITIAL MECHANICAL UNBALANCE
Each rotor system has its own weight distribution, stiffness, and initial shaft deformation.As a result, the amount and location of the mechanical unbalance is very difficult to predict in an arbitrary rotor system.Therefore, a nominal initial mechanical unbalance (U m ) will be assumed.This mechanical unbalance will be defined as the unbalance created from a centrifugal force equal to 10% of the total static rotor weight (W).U m will occur when the rotor is running at its maximum continuous operating speed (ω MCOS ), and can be mathematically defined as: This unbalance will act at an angle of zero degrees with respect to a coordinate system on the rotor, and will be located at the center-of-gravity of the overhung mass.

FILM THICKNESS EXPRESSIONS
In order to establish the static equilibrium position for the synchronous orbit, the film thickness expressions must first be calculated.Cameron (1966) gives the film thickness (h) for a plain journal bearing as: The bearing clearance is C b = R b − R j (R b = bearing radius, R j = journal radius) and the eccentricity ratio ε = e/C b .A film thickness expression was derived for the tilting pad journal bearings: This expression agrees well with the one provided by Monmousseau et al. (1997).In Equation (3), R p = pad radius of curvature, θ c = angle to line joining journal center to bearing center, δ = tilt angle, θ p = pivot angular position, and t p = pad thickness.

STATIC EQUILIBRIUM POSITION
The synchronous orbit develops about a static equilibrium position that was obtained for both the plain journal bearings and the tilting pad ones.A narrow bearing assumption was made for the plain journal bearings.Such an assumption enabled the use of an analytical solution (Cameron, 1966) for the static eccentricity ratio and static attitude angle of these components.However, the complexity of the tilting pad bearings made it necessary to numerically compute the required static equilibrium position.This finite element computation involved the balancing of vertical and horizontal forces acting on the journal and equilibrium of the moments acting on the pads.The resulting system of equations was then solved by using the Newton-Raphson method.

SYNCHRONOUS ORBIT AND HOT/COLD SPOT LOCATIONS
During its synchronous orbit, the journal traverses an elliptical path about its static equilibrium position, O j0 (see Figure 3).Once this configuration is achieved, a general point, P, on the surface of the journal will also travel along an elliptical path.With reference to Figure 10, the path taken by the general point

FIGURE 3
Schematic diagram of the synchronous orbit.
From the geometry in Figure 3, the location of the hot and cold spots can be determined.The hot spot (H) will be defined as the point closest to the bearing sleeve at time t = 0.In this configuration, it is also assumed that γ = ωt + λ so that O j0 , O j , and P = H all lie on the same straight line.Hence, at time t = 0: The time-dependent locations of the hot spot (H (x H , y H )) and the diametrically opposite cold spot (C (x C , y C )) can now be obtained by using the definition of γ 0 : (a) x H = e 0 cos(θ j0 ) + A x cos(ωt + φ x ) + R j cos(ωt + γ 0 ) [8] (b) y H = e 0 sin(θ j0 ) + A y sin(ωt It should be noted that the parameters A x , A y , φ x , and φ y can be obtained directly from the forced response program in VT-FAST (Virginia Tech-Front-end Automated Simulation of Turbomachinery) which is a rotordynamics software package that is currently used in both the classroom and industry.

TEMPERATURE DISTRIBUTION
The next step in the model is to obtain the circumferential temperature distribution at the journal surface within the bearing.In order to obtain this temperature distribution, an energy equation must first be formulated and solved.The physical domain, over which this energy equation is to be applied, can be represented as shown in Figure 4. Figure 4 represents a simple model for energy flow in the bearing.It assumes negligible axial heat flow and steady-state conditions within the control box.It should be noted that the control box encompasses a portion of the journal surface.Hence, the temperature within the box will be the same as the circumferential journal surface temperature.

FIGURE 4
Energy flow in bearing.
Another key assumption in this model is that the final temperature of the journal is obtained by averaging the distributions from each of the dynamic positions.As a result, the average journal temperature distribution can be obtained by considering several points in the synchronous orbit.
Conservation of energy can generally be stated as: Energy Generation Rate = Energy Outflux Rate + Energy Accumulation Rate The steady-state assumption means that the accumulation rate is zero.Therefore, applying conservation of energy to Figure 4 yields: Substituting the appropriate terms into the above equation gives: The rate of viscous energy dissipation, Ėvisc , depends on the speed of the journal surface and the viscous shear stress, τ , which acts over an area dxdz (z = axial dimension).This energy source term heats up the advecting lubricant which has a density of ρ l and a specific heat capacity of c l .A fraction, f, of the remaining heat is transferred to the journal while the rest is assumed to be lost to the bearing housing and surroundings.It is also assumed that the total heat loss is equal to H dxdz(T − T amb ), where H = heat transfer coefficient and T amb = average ambient journal temperature = average ambient bearing temperature.After simplification the above energy equation reduces to: If the lubricant is Newtonian, then the shear stress will be given by τ = µdu/dy.Furthermore, a linear velocity profile can be assumed, i.e., the shear strain rate du/dy = ωR j /h.This Petroff type of simplification has been found to give reasonable accuracy (Cameron, 1966).It should also be noted that since h R j (journal radius), x ∼ R j ξ .The net result of these substitutions yields the following equation: Equation ( 13) can now be solved for the temperature distribution in either a plain journal bearing or in a tilting pad journal bearing.

FIGURE 5
Thermally-induced bend in a rotor with an overhung mass.

THERMAL UNBALANCE
A thermally induced bend in the overhung rotor can be represented as shown in Figure 5.The static bending moment (M) in the journal can be expressed as: where E = Young's Modulus, I = area moment of inertia, and ψ = bend angle as shown in Figure 5.
The maximum stress (σ max ) in the shaft can be written as: If this maximum stress is primarily caused by thermal effects, then: where α = thermal conductivity of the journal and T is the overall mean temperature difference between the hot and cold spots.It is assumed that T is independent of axial position within the bearing.Furthermore, the small deflections allow the thermal and mechanical problems to be uncoupled.Equation ( 14) can be substituted into Equation ( 15) and Equations ( 15) and ( 16) can be equated to give: For a given journal speed and orbit, the right-hand-side of Equation ( 17) is essentially constant.Therefore, Equation ( 17) can be integrated between the start of the bearing-where it is assumed that the bend angle is zero-and some arbitrary axial distance, z (see Figure 5), to give: where y = deflection of the journal.This equation can be evaluated at the end of the bearing to give the bend angle at this location: Also, Equation ( 18) can be integrated over the bearing to give the deflection (y b ) at the end of the bearing: The additional deflection (y d1 ) at the overhung mass location can then be calculated using the small angle assumption: By substituting Equation ( 19) into Equation ( 21) and by adding up Equations ( 20) and ( 21), the total thermal deflection at the overhang can be obtained: Equation ( 22) gives the thermal deflection of the overhung mass.This deflection can be converted into a thermal unbalance, U t : where m d = overhung mass.

INSTABILITY CRITERION
The resultant unbalance (U) from U m and U t can be represented as shown in Figure 6.
With reference to Figure 6, the mechanical and thermal unbalances can be added vectorially to produce the resultant unbalance (U) which can be represented as follows: The angle (ωt − θ CH ) is actually calculated by averaging similar angles at several different points in the synchronous orbit.
A threshold unbalance (U thr ) for this synchronous thermal instability has to be defined.This unbalance was assumed to be where W = rotor weight and ω = variable angular journal speed.The value of 15% was chosen because it gave the best results with the plain journal bearing cases.The rotor will be unstable whenever U exceeds U thr .As a result, the threshold speed for instability (ω thr ) occurs when U = U thr .This instability criterion can be obtained graphically

FIGURE 7
Flowchart for plain and tilting pad journal bearing programs.from the intersection of the U vs. ω and the U thr vs. ω curves, and it determines the onset of the Morton Effect in the rotor system.

PROGRAM FLOWCHART
The flowchart shown in Figure 7 is applicable to both the plain and tilting pad journal bearing cases.At the beginning of the program, the rotor speed is updated and the effective viscosity is calculated.This information is then used to establish a static equilibrium position.The mechanical unbalance (U m ) is next used in VT-FAST to get the dynamic orbit which is combined with the static position to obtain the total synchronous orbit.The temperature difference between the hot and cold spot is then calculated at equal time intervals in the orbit.By averaging these temperature differences, the overall mean temperature difference ( T) between the hot and cold spot is acquired.Next, the unbalances (U t , U, and U thr ) are computed and the process is repeated for each rotor speed.U and U thr can then be compared over the speed range to determine whether the Morton Effect will occur in the rotor.

CONCLUSIONS
A theoretical model was developed for a synchronous thermal instability that acts primarily in overhung rotors.This model utilizes a threshold unbalance criterion for stability instead of employing the traditional frequency-domain feedback analysis.The computer program, which is based on this model, will be used to investigate different case studies in the second part of this paper.

NOMENCLATURE
Dimensions: M = mass, L = length, t = time, T = temperature A x , A y amplitude parameters for elliptic orbit (L) c l lubricant specific heat capacity (L 2 t −2

FIGURE 1
FIGURE 1 Schematic diagram of the Morton Effect model used by other researchers.

FIGURE 2
FIGURE 2Schematic diagram of current Morton Effect model.