# A Review of Tilting Pad Bearing Theory

Authors: Timothy Dimond; Amir Younan; Paul Allaire
Journal: International Journal of Rotating Machinery
 (2011)

---

## Abstract

A theoretical basis for static and dynamic operation of tilting pad journal bearings (TPJBs)
has evolved over the last 50 years. Originally demonstrated by Lund using the pad assembly
method and a classic Reynolds equation solution, the current state of the art includes full
thermoelastohydrodynamic solutions of the generalized Reynolds equation that include fluid
convective inertia effects, pad motions; and thermal and mechanical deformations of the
pads and shaft. The development of TPJB theory is reviewed, emphasizing dynamic modeling. The paper begins with the early analyses of fixed geometry bearings and continues to
modern analyses that include pad motion and stiffness and damping effects. The development
of thermohydrodynamic, thermoelastohydrodynamic, and bulk-flow analyses is reviewed. The
theories of TPJB dynamics, including synchronous and nonsynchronous models, are reviewed.
A discussion of temporal inertia effects in tilting pad bearing is considered. Future
trends are discussed, and a path for experimental verification is proposed.


---

## Body

## 1. Introduction

Rotating machinery such as pumps, compressors, fans, turbines, and generators are ubiquitous in industrial settings. Fluid film journal bearings are suited to applications experiencing higher speeds and loads, often in excess of 90 m/s surface velocity and with typical bearing specific loads of 700 kPa–3.5 MPa. In these bearings, hydrodynamic action is used to support the rotor on a thin lubricating film. Typical film thicknesses between the rotor and the bearing surface are on the order of 100μm. This compares to shaft diameters on the order of 25 mm to 1 m, leading to typical diametral clearance ratios for fluid film bearings of 1-2 μm/mm. With the dimensions and rotational speeds for high-speed applications, lubricating flows are at the low end of the turbulent regime for oil lubricants. There is not a single accepted Reynolds number corresponding to the onset of turbulence in bearing lubrication. The prior literature indicates a range, with turbulence onset taken to occur over a range of Re from 1000 to 1500 [1–5]. Some applications use water or other low-viscosity process fluid as the bearing lubricant, which results in a highly turbulent lubricating flow.The hydrodynamic action in fluid film bearings is fundamentally a fluid-structure interaction effect. When these effects are linearized and perturbed in the two orthogonal radial directions relative to the shaft, they result in equivalent lateral stiffness and damping coefficients which can then be used in lateral vibration analysis of rotors.Tilting pad journal bearings are a source of both static support and dynamic stiffness and damping. Tilting pad journal bearings have a number of pads, typically four or five. A four pad bearing is shown in Figure1. Each pad in the bearing is free to rotate about a pivot and cannot support a moment. As a result, the destabilizing forces are greatly reduced or eliminated, and the bearings are no longer a potential source of rotordynamic instability. This feature has made tilting pad journal bearings the standard fluid-film bearing for most high-speed applications.Figure 1
Tilting pad journal bearing.High-speed rotordynamic applications often have rotors that pass through one or two bending critical speeds as the machines are accelerated to the operating speed. The damping from the fluid film bearings is required to safely pass through these bending critical speeds as the rotating element is accelerated. The damping also helps suppress potentially destabilizing forces from sources such as radial seals, balance pistons, impeller eye seals, internal friction fits, and unbalanced electromagnetic forces [6].Characterization of the dynamic response of tilting pad bearings is vital to successful design of high-speed rotating machinery. Theoretical models exist for prediction of the dynamic response. These are particularly important at the design stage of modern high-speed rotating machinery. These models have evolved from analytical solutions of the lubricating film of fixed geometry bearings to full finite element and finite difference numerical solutions that include analysis of the lubricating flow; the energy balance between the lubricant, the bearing, and the rotor; and mechanical and thermal deformations of the shaft and bearing pads.Additionally, there is still some controversy within the rotordynamic community over the proper dynamic modeling method for tilting pad journal bearings. Some researchers question whether consideration of excitation frequencies other than rotor operating speed is necessary. There is also continued discussion on the number of degrees of freedom to retain within the bearing, either implicitly or explicitly, when nonsynchronous dynamic models are considered. The issue of the presence of fluid temporal inertia effects also arises as part of these discussions. This paper will review and discuss those issues.This paper is organized into nine sections. In Section2, the history of the development of tilting pad bearing dynamic theory is briefly discussed, including a stiffness-damping bearing model. This discussion considers the history of lubrication theory with an emphasis on development of bearing dynamic models. Section 3 reviews initial developments of bearing dynamic models, including key developments in fixed geometry bearings and synchronously reduced tilting pad bearing dynamic models. Section 4 reviews later developments, including nonsynchronous tilting pad bearing dynamic models, thermohydrodynamic (THD) lubrication analysis, and thermoelastohydrodynamic (TEHD) lubrication analysis.Many of the investigators of bearing dynamics were concerned with the onset of rotordynamic instability, so stability assessments in the literature are common. Less common are predictions of critical speeds and unbalance response. The previous work considered in this paper is not comprehensive, but is representative of developments of bearing models. The works cited in this paper and their references do give a comprehensive treatment of bearing modeling developments.Section5 summarizes the current state of the art for the thermoelastohydrodynamic (TEHD) lubrication theory for bearings. The effects of the hydrodynamics due to fluid flow, the energy equation for heat transfer within the bearing, mechanical deflections, and thermal growth are summarized. Section 6 summarizes the averaged flow approaches, including the mixing length theory originally proposed by Constantinescu, and the bulk flow approach proposed by Hirs. In these approaches, the properties of the lubricant film are averaged. The TEHD analysis described in Section 5 is based on a differential equation approach to characterization of the lubricating film and the bearing components. As a result, the flow characteristics are modeled locally, including across the lubricant film. The averaged flow approaches described in Section 6, which can be extended to THD analyses, rely on averaging of flow properties across the lubricating film.Section7 summarizes the current state of the art in tilting pad bearing theory for bearing dynamics. Tilting pad bearings are the only practical bearings used in high-speed flexible rotor industrial machines. The three tilting pad bearing dynamic models that appear recently in the literature are also considered. These dynamic models, (1) the full stiffness-damping model (full KC); (2) the reduced order stiffness-damping-mass (KCM) model; (3) the synchronously reduced model, are described in Sections 5.4–7.4. The frequency-dependent stiffness damping (frequency-dependent KC), an implicit version of the full KC model, is also considered.Section8 discusses two approaches to fluid temporal inertia effects in fluid film bearings and the relative importance and limitations to their analysis.Discussion and conclusions are provided in Section9. This summarizes the evolution of bearing modeling and discusses future trends and opportunities for experimental confirmation.


## 2. Tilting Pad Bearing Modeling Development

### 2.1. Early Lubrication Theory and Nondimensionalization

The fundamental lubrication equation was originally formulated by Reynolds in 1886 [7]. He developed the theory that explained the experimental results of Tower and Petroff. Reynolds assumed that the flow could be treated as isoviscous and laminar. These assumptions were well justified due to typical bearing operating conditions in the late 1800s. He also assumed that since the lubricating film was thin in the radial direction compared to the circumferential and axial directions that there was no pressure gradient across the film radially. The dominant flow characteristic was then due to shearing effects. By simplification of the fluid Navier-Stokes and continuity equations consistent with the assumptions, a single equation describing lubricating flows was derived:(1)ddη(h3dpdη)+ddz(h3dpdz)=6μ[(Uo+U1)dhdη+2V1].Equation (1) is the classic Reynolds equation, expressed in terms of a single pad local coordinate system. In (1), η is the local coordinate direction along the slider, h is the film thickness, p is the developed pressure in the lubricating film, μ is the lubricant dynamic viscosity, z is the axial direction, Uo and U1 represent the relative velocity between slider and pad in the η direction, and V1 represents the squeeze velocity between the slider and the pad. The geometry and velocity vectors for a plain slider is shown in Figure 2. Reynolds noted that since the lubricating film thickness is small compared to the radius of curvature of typical bearings, (1) can be written in a Cartesian coordinate system without significant loss of accuracy. Reynolds was able to provide some series solutions to (1) to find the pressure field within the bearing, but the first closed-form solution to Reynolds equation was found by Sommerfeld [8]. Sommerfeld also developed a group of bearing parameters for nondimensionalization of the results. The Sommerfeld number is(2)So=μNLDW(Dcd)2.Figure 2
Plain slider.The Sommerfeld number recognizes the effect of the net force applied to the shaft at the bearing,W, the projected area of the bearing LD, the clearance ratio cd/D, the rotational speed in rev/s N, and the dynamic viscosity of the lubricating fluid μ as parameters that affect bearing operation in the laminar regime. Many of the theoretical results in the literature use (2) for nondimensionalization.One of the assumptions in arriving at the Sommerfeld number is that the lubricating flow is laminar. While appropriate for lower surface speed bearings where the laminar flow assumption can be justified, the assumption is increasingly violated due to the high rotational speeds typical for many modern bearings. An additional dimensionless group for evaluating bearings operating in the turbulent regime is the Reynolds number using the bearing diametral clearancecd as the characteristic length, and ωD as the characteristic velocity:(3)Re=ρωDcdμ,
where ρ is the density of the lubricant, ω is the shaft rotational speed, and D is the diameter of the shaft. However, the Reynolds number in this form overemphasizes the effect of the shaft diameter, so an alternative reduced Reynolds number is also typically defined:(4)Re*=ρωDcdμ(cdD)=ρωcd2μ.
Both Re and Re* consider the effect of the fluid density but do not include the load on the bearing. However, Re and Re* are the currently accepted methods of nondimensionalizing the lubricating flow turbulence characteristics in modern bearing treatments. Recent papers often report the bearing specific load W/(LD)in conjunction with Re or Re* when describing results instead of using the Sommerfeld number.The analyses of Reynolds and Sommerfeld focused solely on the solution to the hydrodynamic flow field problem, described by (1). Later research included simultaneous solutions of the energy equation, resulting in a thermohydrodynamic (THD) analysis. Another refinement to the problem was to include mechanical deformation effects due to applied pressures and thermal growth, resulting in the thermoelastohydrodynamic (TEHD) theory. The historical development of these theories is discussed briefly in Sections 3-4.


### 2.2. Geometric Considerations

Other bearing design geometric properties will enter the discussion. The dynamic properties of bearings are affected by the bearing design geometry. The geometric properties are as follows.(i)
Bearing preload,m=1-cb/cp. The difference between the bearing clearance cb and the pad clearance cp form a converging hydrodynamic wedge purely through geometry, even for a centered rotor.(ii)
Pivot location: pivot location relative to the leading edge of the pad expressed as a percentage of pad arc length(iii)
Load orientation: napplied load relative to the bearing pads. Load on pad and load between pad configurations are typical.


## 2.1. Early Lubrication Theory and Nondimensionalization

The fundamental lubrication equation was originally formulated by Reynolds in 1886 [7]. He developed the theory that explained the experimental results of Tower and Petroff. Reynolds assumed that the flow could be treated as isoviscous and laminar. These assumptions were well justified due to typical bearing operating conditions in the late 1800s. He also assumed that since the lubricating film was thin in the radial direction compared to the circumferential and axial directions that there was no pressure gradient across the film radially. The dominant flow characteristic was then due to shearing effects. By simplification of the fluid Navier-Stokes and continuity equations consistent with the assumptions, a single equation describing lubricating flows was derived:(1)ddη(h3dpdη)+ddz(h3dpdz)=6μ[(Uo+U1)dhdη+2V1].Equation (1) is the classic Reynolds equation, expressed in terms of a single pad local coordinate system. In (1), η is the local coordinate direction along the slider, h is the film thickness, p is the developed pressure in the lubricating film, μ is the lubricant dynamic viscosity, z is the axial direction, Uo and U1 represent the relative velocity between slider and pad in the η direction, and V1 represents the squeeze velocity between the slider and the pad. The geometry and velocity vectors for a plain slider is shown in Figure 2. Reynolds noted that since the lubricating film thickness is small compared to the radius of curvature of typical bearings, (1) can be written in a Cartesian coordinate system without significant loss of accuracy. Reynolds was able to provide some series solutions to (1) to find the pressure field within the bearing, but the first closed-form solution to Reynolds equation was found by Sommerfeld [8]. Sommerfeld also developed a group of bearing parameters for nondimensionalization of the results. The Sommerfeld number is(2)So=μNLDW(Dcd)2.Figure 2
Plain slider.The Sommerfeld number recognizes the effect of the net force applied to the shaft at the bearing,W, the projected area of the bearing LD, the clearance ratio cd/D, the rotational speed in rev/s N, and the dynamic viscosity of the lubricating fluid μ as parameters that affect bearing operation in the laminar regime. Many of the theoretical results in the literature use (2) for nondimensionalization.One of the assumptions in arriving at the Sommerfeld number is that the lubricating flow is laminar. While appropriate for lower surface speed bearings where the laminar flow assumption can be justified, the assumption is increasingly violated due to the high rotational speeds typical for many modern bearings. An additional dimensionless group for evaluating bearings operating in the turbulent regime is the Reynolds number using the bearing diametral clearancecd as the characteristic length, and ωD as the characteristic velocity:(3)Re=ρωDcdμ,
where ρ is the density of the lubricant, ω is the shaft rotational speed, and D is the diameter of the shaft. However, the Reynolds number in this form overemphasizes the effect of the shaft diameter, so an alternative reduced Reynolds number is also typically defined:(4)Re*=ρωDcdμ(cdD)=ρωcd2μ.
Both Re and Re* consider the effect of the fluid density but do not include the load on the bearing. However, Re and Re* are the currently accepted methods of nondimensionalizing the lubricating flow turbulence characteristics in modern bearing treatments. Recent papers often report the bearing specific load W/(LD)in conjunction with Re or Re* when describing results instead of using the Sommerfeld number.The analyses of Reynolds and Sommerfeld focused solely on the solution to the hydrodynamic flow field problem, described by (1). Later research included simultaneous solutions of the energy equation, resulting in a thermohydrodynamic (THD) analysis. Another refinement to the problem was to include mechanical deformation effects due to applied pressures and thermal growth, resulting in the thermoelastohydrodynamic (TEHD) theory. The historical development of these theories is discussed briefly in Sections 3-4.


## 2.2. Geometric Considerations

Other bearing design geometric properties will enter the discussion. The dynamic properties of bearings are affected by the bearing design geometry. The geometric properties are as follows.(i)
Bearing preload,m=1-cb/cp. The difference between the bearing clearance cb and the pad clearance cp form a converging hydrodynamic wedge purely through geometry, even for a centered rotor.(ii)
Pivot location: pivot location relative to the leading edge of the pad expressed as a percentage of pad arc length(iii)
Load orientation: napplied load relative to the bearing pads. Load on pad and load between pad configurations are typical.


## 3. Fixed Geometry and Synchronously Reduced Bearing Dynamics

The solutions by Reynolds and Sommerfeld were for the pressure field and the net forces of the lubricating oil film. Most analysts of the era considered the rotor to be simply supported at the bearings. As the understanding of the linearized bearing response improved, researchers recognized the equivalent stiffness and damping effects provided by the lubricating film.The first attempts to quantify the dynamic response of the lubricating film itself were made by Stodola [9] and by Stodola's student Hummel [10]. Stodola and Hummel were able to obtain a solution for the oil film stiffness and correctly obtained analytical linearized direct and cross-coupled stiffness terms based on the Sommerfeld closed form solution to (1). However, they did not recognize the damping effect of the oil film, and their predictions indicated that an unstable rotor would have vibration levels that increased without bound. This is a limitation of linear analysis that does not consider the nonlinear behavior of the oil film under large excursions or practical considerations such as contact between the rotor and the stator. Hummel acknowledged in his thesis that the rotor vibration remained finite but did not provide a specific mechanism.Another early analysis that recognized the effect of bearing flexibility on critical speeds was reported by Linn and Prohl [11]. While not referring to oil film flexibility explicitly, a general bearing flexibility was assumed and the resulting lowering of the critical speeds compared to a rigid support assumption was demonstrated.Fixed geometry radial bearings were standard in the first half of the 20th century, and tilting pad bearings only saw significant adoption begin during the 1960s. However, the tilting pad thrust bearing was invented independently by Kingsbury and Michell. Michell also invented the tilting pad journal bearing and installations of the tilting pad journal bearing appear as early as 1916 [12]. An installation of a combined tilting pad journal and thrust bearing on the H.M.S. Mackay of the British Royal Navy, which was launched in 1918, is shown in Figure 3. Fixed geometry bearings remained the standard in the first half of the 20th century due to reduced cost for fixed geometry installations, the higher parasitic losses associated with tilting pad bearings compared to fixed geometry bearings, lower load capacity of tilting pad bearings [13], and lower operating speeds that could tolerate the destabilizing effects of fixed geometry bearings. Boyd and Raimondi in particular stated [14]:Figure 3
Michell Combined Tilting Pad Journal and Thrust Bearing [12].“[T]he plain journal bearing compares favorably with the pivoted-pad bearing and by many criteria is somewhat superior to the latter.”Because of these factors, the perceived drawbacks to tilting pad bearings outweighed the benefits.The advantages of tilting pad bearings in removing the bearings as a source of self-excited vibrations was originally recognized by Hagg in 1946 [15]. Hagg presented experimental results for several tilting pad bearings, including 3-pad, 4-pad, 5-pad, and 6-pad bearings. However, the fluids model that Hagg employed to explain the stabilizing features of these bearings was fundamentally incorrect. A linear flow profile was assumed that did not account for the Reynolds equation, (1).While the development reported by Hagg was significant, many analysts continued to work with plain journal bearings. The benefit of improved stability margin was still not significant enough to designers of the era to overcome the perceived drawbacks discussed previously. Additionally, noncircular bearing bores were discovered to enhance the stability margin and were a lower cost option to a tilting pad arrangement.Sternlicht [16] presented a finite difference solution to the Reynolds equation based on an isoviscous lubricant. The finite difference solution was used to calculate the developed pressure field, which was then integrated to calculate forces. The force solution was then perturbed to determine eight stiffness and damping coefficients based on rotor motion at the journal. These eight coefficients are widely accepted as a proper model for fixed geometry journal bearings where temporal inertia is not important.Solutions to the perturbed Reynolds equation also began to appear in textbooks, including the ones by Smith [17], Pinkus and Sternlicht [18], and Tondl [19]. Smith's treatment of the subject was brief, but did include the eight stiffness and damping coefficients. Pinkus and Sternlicht investigated the stability of rotors supported in plain journals, but the solutions were performed in polar coordinates. Modern rotordynamics analyses are performed in cartesien coordinates for simplicity and for direct comparison to vibration measurements. Tondl's text was an investigation of various sources of rotordynamic instability, of which fixed geometry bearings were a significant contributor. Tondl's treatment accounted for direct stiffness, direct damping, and cross-coupled stiffness terms. His investigation included a treatment of the perturbed Reynolds equation for both linear and nonlinear rotor vibrations. Tondl also considered the benefits of noncircular fixed geometry bearing stator profiles, including lobed bearings.Even with the improvements to fixed geometry designs to enhance stability, there is still a limit where the destabilizing forces are high enough to overcome the damping and drive the rotor unstable. Typically, the limit is reached when the operating frequency is greater than twice the first bending natural frequency [15]. These limits began to be reached on a more consistent basis in the 1960s. The advantage of the stabilizing effect of the tilting pad bearing was then seen to overcome the drawbacks of decreased load capacity and higher parasitic losses.Adoption of tilting pad bearings was also made easier by analytical solutions that allowed designers to understand the dynamic properties. One of the major advances in understanding the dynamics of tilting pad bearings came from Lund's landmark paper in 1964 [20]. Based on analyses of fixed pads, which are essentially partial arc bearings, stiffness, and damping coefficients were calculated. The equations of motion for the pads were then considered based on the calculated fixed pad stiffness and damping values, which were then summed vectorially to calculate the full bearing coefficients. This pad assembly method was not a simultaneous solution of the lubrication problem for all the pads. The pad assembly method is less computationally expensive and more approximate. However, this approach was suitable for the computers available in the 1960s. The dynamic coefficients were then reduced synchronously, or using the shaft rotational frequency as the excitation frequency of interest, to obtain the eight stiffness and damping coefficients related to rotor motion. Results were presented for four-pad, five-pad, and six-pad tilting pad bearings. Lund recast rotations of the pads as equivalent translations. Later solutions treated pad rotations as rotational motion.The landmark work by Lund led to a significant research effort to extend the analyses of tilting pad bearings. Thermohydrodynamic and TEHD solutions and turbulence corrections for high rotational speeds also begin to appear for tilting pad bearings. The use of synchronously reduced coefficients to described the linearized dynamics became the norm, following the results reported by Lund [20].Orcutt [21] followed the same basic approach as Lund [20] by developing a partial arc bearing solution. He accounted for turbulence effects in the lubricating film using the analysis of Ng and Pan [22]. Similar to Lund, Orcutt solved the lubrication problem for each pad individually and then performed a synchronously reduced assembly method similar to Lund. While not a comprehensive formulation, modifications to (1) were included to account for turbulence effects in the lubricating film without resorting to a full Navier-Stokes solution. Orcutt considered different lubricants, different numbers of pads, and different pad preloads. Orcutt's analysis indicated that symmetry in the tilting pad bearing leads to isotropic bearing dynamic properties. This is not strictly correct, since a simultaneous TEHD solution and operating experience indicates differential heating of the tilting pads, but symmetric tilting pad bearings such as four-pad bearings in a load-between-pad configuration are nearly isotropic. Orcutt's results were plotted against the Sommerfeld number. He also suggested that the results showed that operation above the first bending critical was possible, which was generally avoided by designers of the era to avoid stability problems. Including pad preload was also claimed to improve dynamic characteristics. It was shown later [23, 24] that low preload leads to more stable systems.Nicholas et al. [25] developed stiffness and damping coefficients for the five pad tilting pad journal bearing. Several bearing configurations were considered, including load on pad and load between pad, different pivot offsets ranging from 0.5 to 0.55 and different bearing preloads from 0 to 0.5. The reported effective stiffness and damping coefficients were synchronously reduced. A finite element formulation of (1) was used, so it was an isoviscous, laminar analysis. The finite element method was used to produce single pad solutions with pad assembly similar to Lund [20]. The reduced stiffness and damping coefficients were plotted versus Sommerfeld number. The pad inertia effects were neglected in the analysis.Nicholas and Kirk [26] examined several fixed and tilting pad bearings, including four-pad and five-pad tilting pad bearings, for application to axial compressors. The synchronously reduced stiffness and damping coefficients were used for unbalance and stability analyses. This paper explored the effect of manufacturing tolerances on the performance of the several bearing types. For a four-pad tilting pad bearing the synchronously reduced stiffness was shown to vary by 50 percent and synchronously reduced damping by half an order of magnitude due to typical manufacturing tolerances, though that was an extreme case. Variations of 10 percent in the dynamic coefficients were more typical. The effect of bearing preload on the compressor stability margin was investigated, and lower preload was determined to enhance stability because the bearing stiffness was reduced, allowing more motion for the damping to be effective. The four-pad bearing configuration for axial compressors was explored further in [27]. Nicholas and Kirk again used the synchronously reduced bearing dynamic coefficients for both unbalance response and stability margin.Jones and Martin [28] performed another geometric study of tilting pad bearing characteristics, considering different preloads; bearing L/D ratios; 3, 5, and 7 tilting pads; and load orientation. The analysis was used to calculate minimum oil film thickness, average pad temperatures, bearing parasitic power losses, and the synchronously reduced stiffness and damping coefficients. The analysis was described as quasi-THD, since average pad temperature was used to calculate the average oil viscosity for each pad. The modeling included turbulence effects and was performed using finite difference methods.Ettles [13] developed a TEHD analysis of tilting pad bearings. The analysis included a generalized Reynolds equation solution using the turbulence model of Constantinescu [29] and the local calculated Reynolds number to obtain an effective viscosity. The energy equation was simplified into a 1D solution, and the relative error compared to a 2D solution was calculated to have a maximum value of 3.52 percent for an L/D ratio of 9.9. Elastic deformation of the pads due to applied loads and thermal expansion were also considered. Ettles' solution was a simultaneous, iterative solution for all the bearing pads. The nondimensional dynamic coefficients were reduced synchronously and were plotted as a function of Sommerfeld number. The results were compared to dynamic experiments by Malcher [30] for four pad bearings. The theoretical results were within 10 percent of the reported measured values. The reduction of effective film stiffness and damping due to pad flexibility was noted. A stability analysis of the pad motion as a check for pad flutter was also included.Hashimoto et al. [31] also developed a TEHD analysis suitable to large scale tilting pad bearings with two pads. The large generator application could be supported on two lower bearing pads in a load between pad configuration, so the top pads were eliminated. The generalized Reynolds equation with a turbulence model relating effective viscosity to local Reynolds number was solved simultaneously with a 1D energy equation and a deformation model for the bearing pads. Pad preloads of 0, 0.1, and 0.2 were considered, and the pads were centrally pivoted for all cases. The results assuming a laminar fluid were compared to the result obtained when the turbulence model was included. The bearing stiffness for the turbulent case was up to 20 percent lower for Sommerfeld number up to about 0.4. For higher Sommerfeld numbers, the stiffness for the turbulent case was higher than the laminar case by as much as 100 percent. A similar trend was observed in the damping coefficients, but the crossover point was at a Sommerfeld number of 0.1-0.2. The dynamic coefficient reduction method was not explicitly stated but the reported coefficients are consistent with synchronous reduction.Knight and Barrett [32] presented a THD analysis of four pad tilting pad bearings with central pivots in a load on pad configuration. No turbulence model was considered in calculating the results. The solution was based on a finite element solution to the classic Reynolds equation combined with a finite difference solution to the 2D energy equation. A simultaneous solution for the pads was developed. The effective viscosity for each bearing pad was based on the cross-film average temperature for each bearing pad. This average viscosity was then treated as constant in the Reynolds solution. Heat transfer through the pads was treated as radial conduction and shaft surface temperatures were based on the overall average film temperature. The full bearing coefficients were calculated, which were reduced synchronously for presentation of the results. When compared to an isothermal calculation, a difference of 10–35 percent in dynamic properties was reported. This demonstrated the effect of temperature and viscosity on the stiffness and damping coefficients.Brugier and Pascal [33] also investigated tilting pad bearings for large turbogenerator sets. The bearings considered had three tilting pads. The TEHD model accounted for thermal effects in the lubricant as well as mechanical and thermal deformations of the bearing pads. The model included a generalized Reynolds equation and energy equation similar to that described in detail in Section 5. The solution for the dynamic coefficients was based on numerical differentiation and simulated shaft perturbations within the code. The description of coefficient extraction method was unclear, but the plots in the paper are consistent with synchronously reduced coefficients. The solution was obtained using finite difference techniques with overrelaxation. There was no indication that a turbulence model was used. The paper also showed a drop in effective stiffness and damping due to pad deformation, which acts like a spring in series with the oil film.Ettles [34] presented another THD analysis of tilting pad bearings. The synchronously reduced coefficients were calculated and compared to results reported by Brockwell and Dmochowski [35]. Ettles considered the transition region for turbulence in the lubricating flow to be in the range of 1100–1400. Turbulence was accounted for by lowering the inlet oil temperature into the pad, which was justified by measurements showing the improved heat transfer in the bearing due to the onset of turbulence. As a result, the treatment of turbulence was purely empirical with no formal turbulence model. The model also accounted for thermal and mechanical deformations of the bearing pads and showed the resulting drop in effective stiffness and damping due to deformation. A finite difference solution was employed. Fillon et al. [36] also demonstrated the need to consider bearing element deformation to obtain accurate predictions of bearing behavior.Brockwell et al. [37] developed a THD solution similar to that of Ettles [34]. The two-dimensional THD model included pad thermal expansion and elastic deformations, along with pivot flexibility. Beam theory was used to model pad flexibility. Nondimensional forms of the hydrodynamic, energy, and heat transfer equations were presented. The viscosity terms were averaged across the film thickness, resulting in a bulk-flow approach to the lubrication problem. Dynamic coefficients were calculated for synchronous excitations based on a perturbed Reynolds solution. Static and dynamic results were compared to a test of a five-pad tilting pad bearing in a load-between-pad configuration. The running speeds tested were 900, 1800, 2700, and 3600 rpm. Qualitative trends were matched by the predictions, with differences of 10–15 percent between theory and experiment for equivalent bearing stiffness and damping. The agreement was improved by including the effect of shaft flexibility in the overall identification and analysis procedures.Hopf and Schüeler [38] performed investigations of the transition from laminar to turbulent flow in large turbine bearings. The THD analysis was not described in detail, but was able to predict measured temperatures in the bearing tilting pads. The sudden drop in bearing local temperature due to the onset of turbulence or Taylor vortices and the resulting enhanced heat transfer was documented.Hyun et al. [39] also developed a THD solution to the tilting pad bearing performance. The generalized Reynolds equation was solved, using Reichardt's wall formula and the turbulence model presented by Ng and Pan [22] to model turbulence. A three-dimensional model was used for the energy equation, with correction factors for cavitation. Predictions were compared to experimental measurements of film pressure, film temperature, pad temperature, and load capacity, with agreement within 5 percent for a four-pad tilting pad bearing in a load-between-pad configuration.Nicholas and Wygant [40] presented results and design considerations for highly loaded bearings, with specific loads of up to 3.45 MPa (500 psi). The paper focused on pivot design, specifically steel pivots with bronze pads and steel pivots with steel pads. The synchronously reduced stiffness and damping coefficients were used to facilitate the discussion on the effect of pivot stiffness on the dynamic coefficients. The effective stiffness and damping were lowered since the pivot acts as an additional spring in series with the oil film. Hertzian contact theory was used to calculate the effective pivot stiffness.Several researchers have investigated transient effects in tilting pad bearings. These studies were influenced in part by a review of bearing failures presented by Conway-Jones and Leopard [41]. In their review of failure modes, Conway-Jones and Leopard determined that thermal transients as a function of oil inlet temperature and rotational speed led to bearing seizure failures. A theoretical study of the phenomena was performed by Monmousseau and Fillon [42]. Their TEHD analysis was based on the generalized Reynolds equation described in Section 5, and transient forms of the energy and heat transfer equations. Thermal growth was modeled using a plane stress assumption and free boundary conditions on expansion. Finite difference techniques were used in the solution. The pivot flexibility model developed by Kirk and Reedy [43] was also used in the model. The theoretical study predicted seizing after 53 seconds of operation for a rotor was accelerated from 0 to 10,000 rpm over 5 s, with oil lubricant inlet temperature of 20°C. The time to seizure was extended or eliminated by reducing the rate of shaft acceleration and increasing the oil inlet temperature in the analysis.Transient effects were also considered by Monmousseau et al. [44] in the analysis of tilting pad bearings. The TEHD analysis again included the generalized Reynolds equation, and transient forms of the energy and heat transfer equations, with thermal growth included in the overall model. Finite difference techniques were used in the solution. A step change in bearing specific load was modeled. Predictions were compared to pad babbitt temperature measurements, with 10–15 percent difference using the TEHD model. The thermal transient period was on the order of 60 s, while the mechanical transient period was on the order of one shaft rotation at a running speed of 4,000 rpm. The authors concluded that thermal and elastic effects should be modeled to accurately capture bearing transient behavior.


## 4. Nonsynchronous Bearing Model Development

The work by Lund in 1964 [20] reported synchronously reduced bearing coefficients for the tilting pad journal bearing. While appropriate for unbalance response analysis since the unbalance forcing frequency is driven by shaft rotation, the synchronous coefficients are not appropriate in general. The proper reduction method is based on overall system excitation frequency which is in general nonsynchronous with shaft rotation. This distinction was made clear during the presentation of Nicholas et al. [23], which reported results for the stability analysis of an 11-stage compressor. Nicholas et al. used the synchronously reduced bearing coefficients to estimate the compressor stability margin. During the discussion of the paper, Lund told the presenter that the use of synchronously reduced coefficients for stability margin was mathematically incorrect and that reduction at the rotor natural frequency was correct [24]. This comment spurred much research into the development of nonsynchronous bearing dynamic models.Shapiro and Colsher [45] examined the effect of bearing preload on the dynamic response to tilting pad bearings. They addressed the fact that synchronous reduction was typical at the time but indicated that it would lead to erroneous stability predictions. This allowed for the possibility of nonsynchronous dynamic reduction of bearing coefficients. However, the analysis only considered the onset of instability, where the real part of the rotor-bearing system eigenvalues is zero. Only the imaginary part of the eigenvalue was considered in the reduction, which is not the most general linear solution. The most general solution includes the real part of the eigenvalue. Shapiro and Colsher presented the full stiffness and damping matrices for a five-pad tilting pad bearing, resulting in 28 stiffness and 28 damping coefficients.Allaire et al. [46] presented a pad assembly method for tilting pad bearings that explicitly included the motion of the pads and the resulting stiffness and damping coefficients. The solution was based on perturbing the shaft or pad as appropriate for several applied loads and static eccentricities to develop a table of data. The overall equilibrium point and resulting dynamic coefficients were then found through linear interpolation. This method is valid for isoviscous, laminar analyses but needs modification to account for thermal and turbulence effects. The analysis did not account for hot oil carryover, where oil exiting one pad affects the oil entering the next pad. The full coefficient matrix for tilting pad bearings was developed which was independent of the pad inertia and excitation frequency. Results were presented for a five-pad bearing with load on pad and zero-pad preload. Most of the plots were of the full stiffness and damping coefficients as a function of Sommerfeld number, but one set of synchronously reduced plots was presented.Parsell et al. [47] further explored the frequency effects in tilting pad bearings. The authors postulated that synchronously reduced bearing coefficients may be an acceptable engineering approximation for rotordynamic stability, though it was mathematically incorrect. The synchronous coefficients have since been shown to give a nonconservative estimate of stability, for example, [48]. The main purpose of the paper was to plot reduced bearing stiffness and damping coefficients as a function of whirl frequency ratio, which is the ratio of excitation frequency to running speed. Sommerfeld numbers of 0.1, 1, and 10 were considered for five-pad load between pad configurations. Preloads of 0 and 0.3 were also considered. The frequency dependence was reduced for low Sommerfeld number and for high preload. For high Sommerfeld number and zero preload, the reduced coefficients were highly frequency dependent, and effective stiffness and damping approached zero for whirl frequency ratios from 0.3 to 0.5. This effect was postulated later [49] to be due to light bearing load compared to running speed. In that case, the shaft runs nearly centered in the bearing, so negligible pressure forces act on the shaft.Rouch [50] presented a method for modeling pivot flexibility as a spring in series with the effective oil film dynamic coefficients. The pad assembly method was used to determine the oil film dynamic characteristics. A full matrix including the effect of pad rotations and pad translations was presented. Dynamic reduction was performed nonsynchronously. Plots of effective reduced stiffness and damping as a function of excitation frequency were presented. The pivot stiffness reduced the effective stiffness and damping by up to 50 percent due to a spring in series with the oil film. The effect was most pronounced when the effective pivot stiffness was the same order of magnitude as the oil film stiffness. A stability analysis of a flexible rotor was also performed using the nonsynchronous and synchronous bearing coefficients. It was shown that the nonsynchronous coefficients gave a more conservative estimate of rotor stability compared to synchronously reduced coefficients.Lund and Pedersen [51] presented an approximate method for including pad deformations and pivot flexibility in the overall dynamic response of tilting pad bearings. An isoviscous solution to (1) was obtained using a finite difference method. The pad deformations were approximated as the deformation of a beam under a distributed pressure load, and the pivot flexibility was obtained using Hertzian contact theory. The model incorporated the effect of pad deformation by an increased effective clearance between the pad and the shaft. The effective clearance was also treated as dynamic and harmonic, using the vibrational motion of the pad deflections to calculate the effective clearance. The pivot flexibility was treated as a spring in series, and the overall bearing coefficients were dynamically reduced. The authors advocated nonsynchronous reduction of bearing coefficients in general, but synchronously reduced bearing coefficients were presented in the results. The reduction in stiffness and damping due to flexibility effects by up to 50 percent was demonstrated.Branagan [52] presented a TEHD finite element solution for fixed geometry and tilting pad journal bearings. Polynomial profiles for the thermal and viscosity solutions were assumed in the axial and cross-film directions, leaving a 1D solution in the circumferential direction. A simultaneous, iterative solution procedure for the Reynolds equation, the energy equation, and the deformation model was performed to obtain accurate estimates of the changes in boundary conditions due to hot oil carryover. The full bearing coefficients were calculated, and a damped eigenvalue analysis was used to reduce the full bearing coefficients to the eight coefficients related to the shaft degrees of freedom. The dynamic reduction was performed to improve the run times for subsequent rotordynamic analyses. Reduction of the computational expense for rotordynamic models in this fashion is not important with modern computers, where a full eigenvalue analysis of a rotor beam model with full bearing coefficients has a run time of less than 10 s [53]. It was demonstrated that inclusion of pad and pivot flexibility effects could reduce the calculated stiffness and damping coefficients by up to 50 percent.Barrett et al. [54] provided an extension to reduction of full tilting pad bearing coefficients to the eight frequency-dependent stiffness and coefficients. In their analysis, Barrett et al. included the real part of the eigenvalue to allow for general damped analyses. The analysis considered bearings in a load between pad and treated pad inertia as negligible. Single pad solutions similar to Lund [20] were developed, and the pad assembly method was applied. Pad rotations were treated as equivalent translations due to small angle perturbations. The dynamic reduction was performed in the pad local coordinate system for each pad, with transformation to global coordinates as the final step. An excitation frequency dependence in reduced stiffness and damping coefficients was demonstrated. The bearings considered had effective stiffness that decreased by about 20 percent from 0.5X to 1X, where X is the excitation frequency corresponding to shaft rotational speed. A difference in effective bearing coefficients depending on the distance from the stability margin was also demonstrated. This demonstrated that consideration of both the real and imaginary parts of the eigenvalues give a more general treatment in stability analysis, since the only case that has a purely imaginary eigenvalue is the neutrally stable solution. The neutrally stable solution is appropriate for the onset of instability, but is not appropriate for stable or unstable systems.Earles et al. [55] developed a finite element solution for a single tilting pad including the lubricating film and pad deformation effects. The lubricating flow was treated as laminar, isoviscous, and incompressible. Thermal effects were not considered. The pad was modeled using plane strain isoparametric finite elements, and the pressure solution for the lubricant was solved simultaneously to determine the flow field and the final pad dimensions. The pad degrees of freedom were reduced to a single coordinate using Guyan techniques; the single coordinate represented the final pad radius of curvature. The stiffness and damping coefficients were found through numerical perturbation of the shaft and pad positions and deformations. The result for the single pad was then transformed to local coordinates, which implied usage of the results for pad assembly solutions. The mechanical deformations were shown to lower effective stiffness and damping coefficients and were within 5 percent of the coefficients reported by Lund and Pedersen. A damped eigenvalue solution similar to Barrett et al. [54] was employed. The single pad model was extended to a full tilting pad bearing model with lubricating film and pivot flexibility effects [56]. The full stiffness and damping matrices with flexibility effects were presented. Additional terms to account for pivot flexibility and pad deformation effects were incorporated into the global bearing stiffness and damping matrices. The full coefficients were used to perform a damped eigenvalue analysis.White and Chan [57] presented a finite element THD bearing analysis with turbulence correction. Turbulence correction factors were based on bulk flow theory proposed by Hirs [58]. The full stiffness and damping coefficients were found from the perturbed Reynold's equation and the reduction method of Parsell et al. [47] was used for reduction to the eight stiffness and damping coefficients. The Parsell et al. method only uses the imaginary part of the eigenvalue to perform the dynamic reduction. The method is correct for forced response analyses and analyses to determine the onset of instability, but is not correct for a general free response analysis with damped eigenvalues. White and Chan compared the effective stiffness and damping coefficients for synchronous and half-whirl reduction frequencies and showed a reduction in effective damping for half-frequency whirl of up to 20 percent. They also showed that the effective damping was reduced for off-center pivots. Nondimensional dynamic coefficients were plotted as a function of Sommerfeld number.Brockett and Barrett [59] presented a tilting pad bearing dynamic reduction method suitable for transfer matrix analyses. The reduction resulted in a second-order transfer function with a fourth-order residual frequency dependent stiffness. The reduction admitted damped eigenvalue solutions. Results of the bearing representation in a transfer matrix model were compared to a finite element solution for a flexible rotor with the full bearing coefficients modeled explicitly. Agreement within 1 percent of the system eigenvalues was obtained, but the transfer matrix method missed highly damped modes.Kim et al. [60] presented a nonsynchronous reduction of tilting pad bearing coefficients, including pad deformation effects. The pad deformation model reduction was achieved through modal representation of pad movement, with modal truncation. Further dynamic reduction of the oil film and mechanical flexibility effects was performed nonsynchronously to obtain the eight frequency-dependent stiffness and damping coefficients. The analysis accounted for thermal effects on lubricant viscosity as well as pivot flexibility, pad rotations, and pad deformations. A turbulence model was not included. Variable viscosity due to temperature changes was accounted for. The theoretical results were compared to results reported by Brockwell et al. [37], and the analysis predicted the measured drop in stiffness and damping due to mechanical deformations. The theory agreed with the Brockwell et al. data within 10 percent. Theoretical results were also compared to experimental results reported by Fillon et al. [61]. The effects of the TEHD model on the coefficients was shown in separate plots with synchronously reduced coefficients. A stability estimate for an eight-stage gas reinjection compressor was performed and compared to results presented by Wilson and Barrett [62]. It was shown that use of the frequency-dependent stiffness and damping coefficients in the rotor bearing model resulted in a lower stability margin compared to a synchronously reduced bearing model, which agreed with the Wilson and Barrett results.


## 5. Thermoelastohydrodynamic Tilting Pad Bearing Lubrication Theory

Tilting pad bearing lubrication theory has evolved, from fixed geometry isoviscous analytical solutions, to advanced finite element solutions including hydrodynamic, energy, and deformation effects. The modifications to (1) and additional equations to model the energy balance and turbulence modeling are summarized in Sections 5.1–5.3.


### 5.1. Generalized Reynolds Equation

Modern tilting pad bearing lubrication theory is based on thermoelastohydrodynamic models that include equations describing the hydrodynamic flows, heat transfer and shear heating, and mechanical deformations [63–66]. Whiles these solutions evaluate temperature effects in the lubricant film and pad deformations, they do not involve the advanced elastohydrodynamic solutions found for ball or roller bearings. The following discussion does not consider mechanical deformations, which are geometry specific. A very brief discussion of the key equations follows. A comprehensive derivation of the equations is provided in [63]. While more recent work has further refined the TEHD analysis, particularly with the inclusion of a 3D energy equation [67], these refinements have only been applied to fixed geometry bearing analyses.Reynolds' equation, (1), is the fundamental equation for lubricating flows assuming a laminar, isoviscous lubricant. The generalized Reynolds equation results from the fluid continuity and momentum equations, with the assumption that the pressure profile is constant across the lubricating film. The generalized form includes convective inertia effects through an eddy-viscosity model and allows for cross-film variations in viscosity. The formulation of the generalized Reynolds equation in terms of pad local coordinates is [63–66](5)∂∂η{h3Γ(η,z,Re*)∂p∂η}+∂∂z{h3Γ(η,z,Re*)∂p∂z}=-UG(η,z,Re*)∂h∂η.
In (5)), Γ and G represent generalized local effective viscosity functions with turbulence effects included, and U represents the motion of the journal relative to the bearing pad. Equation (5)) has been modified from [63] to include the effect of reduced Reynolds number Re*=(ρωh2)/μ, where ρ is the lubricant density and ω is the rotational speed. This is to account for turbulence effects, especially those due to the low viscosity of the lubricating fluid for some process fluid lubricated bearings. For example, the viscosity of water is two orders of magnitude less than the viscosity of oil [68].The effective cross-film viscositiesΓ and G, as a function of the journal radial and axial positions, are given by [63]
(6)Γ(η,z,Re*)=∫01[χ2(η,1,z,Re*)χ1(η,1,z,Re*)χ2(η,ξ,z,Re*)-χ2(η,1,z,Re*)χ1(η,1,z,Re*)χ1(η,ξ,z,Re*)]dξ,G(η,z,Re*)=1χ1(η,1,z,Re*)∫01χ1(η,ξ,z,Re*)dξ,χ1(η,ξ,z)=∫0ξ1μe(η,ξ′,z,Re*)dξ′,χ2(η,ξ,z,Re*)=∫0ξξ′μe(η,ξ′,z,Re*)dξ′,
where χi represent intermediate viscosity functions, μe represents the effective turbulent viscosity, ξ represents the local pad squeeze direction, and ξ′ represents the dummy variable of integration.The flow profile in the bearing is treated as a combination of Couette and Poiseuille flow, which is expressed as:(7)u(η,ξ,z)=∂p∂ηχ2(η,ξ,z,Re*)+{Uχ1(η,h,z,Re*)-∂p∂η[χ2(η,h,z,Re*)χ1(η,h,z,Re*)]}×χ1(η,ξ,z,Re*),(8)w(η,ξ,z,Re*)=∂p∂zχ2(η,ξ,z,Re*)-∂p∂z[χ2(η,h,z,Re*)χ1(η,h,z,Re*)]χ1(η,ξ,z,Re*),
where u represents local velocity in the sliding direction, w represents local velocity in the axial direction, and U represents the surface velocity of the shaft relative to the pad.The solutions to (5)–(8) result in the pressure field developed in the lubricating film as a function of radial and axial position. The pressure across the film is assumed constant since it is small compared to the radial and axial length scales. Integrating the pressure field over the area of the bearing surfaces yields the net forces in the bearing.


### 5.2. Energy Equation

The viscosity of many lubricants is a strong function of temperature. The developed pressures in hydrodynamic bearings are not large enough to significantly affect the viscosity. To account for temperature effects, the 2D energy equation, including shear heating terms, is considered in the model presented in [63–66].(9)ρCp(u∂T∂η+v∂T∂ξ)=∂∂η(k∂T∂η)+∂∂ξ(ke∂T∂η)+μe[(∂u∂ξ)2+(∂w∂ξ)2],
where v is the fluid velocity in the squeeze direction, Cp is the specific heat of the lubricant, T is the lubricant temperature, k is the heat conductivity of the lubricant, and ke is the effective heat conductivity of the lubricant corrected for turbulence. The relative importance of the shear heating term, μe[(∂u/∂ξ)2+(∂w/∂ξ)2], in (9) is dependent on the lubricant considered. For oil-lubricated bearings operated at high speed, shear heating effects are significant and viscosity variation due to temperature must be considered in the analysis. For many process fluid lubricants such as water, the shear heating effects can often be neglected due to the low lubricant viscosity, allowing for isoviscous analyses.


### 5.3. TEHD Solution Turbulence Modeling

The turbulence model implemented by He [63] was based on models developed by Elrod and Ng [69], with a modification by Suganami and Szeri [1] to account for transition flow. The turbulence is modeled using an eddy viscosity law. The effective viscosity is found by.(10)μe(η,ξ,z)=μ(1+ϵmν),
where ν is the kinematic viscosity. The eddy viscosity ϵm is calculated from(11)ϵmν=κ[ξ+-δl+tanh(ξ+δl+)].
In (11), based on turbulent boundary layer theory, κ=0.4 and δl+=10.7. These constants were found empirically and are reported in Elrod and Ng [69]. The nondimensional distance from the wall ξ+ is defined as:(12)ξ+=ξν|τ|ρ.
And the local shear stress in the lubricant in (12) is found from:(13)τ=μe(∂u∂ξ)2+(∂u∂ξ)2.
The inclusion of the effective viscosity μe is implicit in (10)–(13), so the shear stress distribution and the effective local viscosity in the turbulent regime is typically found iteratively, such as in the finite element code developed by He [63]. By including the eddy-viscosity model in the TEHD code, convective inertia effects are approximated in the solution. However, temporal inertia effects are not currently considered in the analysis.Equations (10)–(13) are valid for fully developed turbulent flows [69]. To account for low levels of turbulence and transitional flows, Suganami and Szeri proposed an additional factor γ, which modifies (10)) as [1].(14)μe(η,ξ,z)=μ(1+γϵmν).
The factor γ is dependent on the maximum Reynolds' number in the lubricating flow and is defined as.(15)γ={0,Remax⁡<ReL,1-[ReH-Remax⁡ReH-ReL]1/8,ReL≤Remax⁡≤ReH1,ReH<Remax⁡,,
where Remax=(ρuh)/μ is the maximum local Reynolds' number in the lubricant and u is the local fluid velocity, ReL is the critical Reynolds' number for the onset of transition flow from laminar flow and ReH is the critical Reynolds' number for the onset of turbulence. There is not a consensus in the literature on the critical Reynolds' numbers ReL and ReH, corresponding to the onset of transition flow and turbulent flow, respectively [63]. Proposed values of Reynolds' number for the onset of transition flow range from about 500–1000 and range from about 800–1500 for the onset of turbulent flow [1–5].


### 5.4. Perturbed Reynolds Equation

Once the pressure profile solution is found, the generalized Reynolds equation is perturbed. The first-order perturbation results in the equivalent stiffness and damping coefficients,kij and cij, respectively, which have the general form [63–66].(16)kij=-∂f∂uij,cij=-∂f∂u̇ij.
The specific stiffness and damping coefficients are defined by He [63].When (5)–(16) are considered, it is apparent that the turbulence model chosen significantly affects the predicted dynamic coefficients for the tilting pad bearing.


## 5.1. Generalized Reynolds Equation

Modern tilting pad bearing lubrication theory is based on thermoelastohydrodynamic models that include equations describing the hydrodynamic flows, heat transfer and shear heating, and mechanical deformations [63–66]. Whiles these solutions evaluate temperature effects in the lubricant film and pad deformations, they do not involve the advanced elastohydrodynamic solutions found for ball or roller bearings. The following discussion does not consider mechanical deformations, which are geometry specific. A very brief discussion of the key equations follows. A comprehensive derivation of the equations is provided in [63]. While more recent work has further refined the TEHD analysis, particularly with the inclusion of a 3D energy equation [67], these refinements have only been applied to fixed geometry bearing analyses.Reynolds' equation, (1), is the fundamental equation for lubricating flows assuming a laminar, isoviscous lubricant. The generalized Reynolds equation results from the fluid continuity and momentum equations, with the assumption that the pressure profile is constant across the lubricating film. The generalized form includes convective inertia effects through an eddy-viscosity model and allows for cross-film variations in viscosity. The formulation of the generalized Reynolds equation in terms of pad local coordinates is [63–66](5)∂∂η{h3Γ(η,z,Re*)∂p∂η}+∂∂z{h3Γ(η,z,Re*)∂p∂z}=-UG(η,z,Re*)∂h∂η.
In (5)), Γ and G represent generalized local effective viscosity functions with turbulence effects included, and U represents the motion of the journal relative to the bearing pad. Equation (5)) has been modified from [63] to include the effect of reduced Reynolds number Re*=(ρωh2)/μ, where ρ is the lubricant density and ω is the rotational speed. This is to account for turbulence effects, especially those due to the low viscosity of the lubricating fluid for some process fluid lubricated bearings. For example, the viscosity of water is two orders of magnitude less than the viscosity of oil [68].The effective cross-film viscositiesΓ and G, as a function of the journal radial and axial positions, are given by [63]
(6)Γ(η,z,Re*)=∫01[χ2(η,1,z,Re*)χ1(η,1,z,Re*)χ2(η,ξ,z,Re*)-χ2(η,1,z,Re*)χ1(η,1,z,Re*)χ1(η,ξ,z,Re*)]dξ,G(η,z,Re*)=1χ1(η,1,z,Re*)∫01χ1(η,ξ,z,Re*)dξ,χ1(η,ξ,z)=∫0ξ1μe(η,ξ′,z,Re*)dξ′,χ2(η,ξ,z,Re*)=∫0ξξ′μe(η,ξ′,z,Re*)dξ′,
where χi represent intermediate viscosity functions, μe represents the effective turbulent viscosity, ξ represents the local pad squeeze direction, and ξ′ represents the dummy variable of integration.The flow profile in the bearing is treated as a combination of Couette and Poiseuille flow, which is expressed as:(7)u(η,ξ,z)=∂p∂ηχ2(η,ξ,z,Re*)+{Uχ1(η,h,z,Re*)-∂p∂η[χ2(η,h,z,Re*)χ1(η,h,z,Re*)]}×χ1(η,ξ,z,Re*),(8)w(η,ξ,z,Re*)=∂p∂zχ2(η,ξ,z,Re*)-∂p∂z[χ2(η,h,z,Re*)χ1(η,h,z,Re*)]χ1(η,ξ,z,Re*),
where u represents local velocity in the sliding direction, w represents local velocity in the axial direction, and U represents the surface velocity of the shaft relative to the pad.The solutions to (5)–(8) result in the pressure field developed in the lubricating film as a function of radial and axial position. The pressure across the film is assumed constant since it is small compared to the radial and axial length scales. Integrating the pressure field over the area of the bearing surfaces yields the net forces in the bearing.


## 5.2. Energy Equation

The viscosity of many lubricants is a strong function of temperature. The developed pressures in hydrodynamic bearings are not large enough to significantly affect the viscosity. To account for temperature effects, the 2D energy equation, including shear heating terms, is considered in the model presented in [63–66].(9)ρCp(u∂T∂η+v∂T∂ξ)=∂∂η(k∂T∂η)+∂∂ξ(ke∂T∂η)+μe[(∂u∂ξ)2+(∂w∂ξ)2],
where v is the fluid velocity in the squeeze direction, Cp is the specific heat of the lubricant, T is the lubricant temperature, k is the heat conductivity of the lubricant, and ke is the effective heat conductivity of the lubricant corrected for turbulence. The relative importance of the shear heating term, μe[(∂u/∂ξ)2+(∂w/∂ξ)2], in (9) is dependent on the lubricant considered. For oil-lubricated bearings operated at high speed, shear heating effects are significant and viscosity variation due to temperature must be considered in the analysis. For many process fluid lubricants such as water, the shear heating effects can often be neglected due to the low lubricant viscosity, allowing for isoviscous analyses.


## 5.3. TEHD Solution Turbulence Modeling

The turbulence model implemented by He [63] was based on models developed by Elrod and Ng [69], with a modification by Suganami and Szeri [1] to account for transition flow. The turbulence is modeled using an eddy viscosity law. The effective viscosity is found by.(10)μe(η,ξ,z)=μ(1+ϵmν),
where ν is the kinematic viscosity. The eddy viscosity ϵm is calculated from(11)ϵmν=κ[ξ+-δl+tanh(ξ+δl+)].
In (11), based on turbulent boundary layer theory, κ=0.4 and δl+=10.7. These constants were found empirically and are reported in Elrod and Ng [69]. The nondimensional distance from the wall ξ+ is defined as:(12)ξ+=ξν|τ|ρ.
And the local shear stress in the lubricant in (12) is found from:(13)τ=μe(∂u∂ξ)2+(∂u∂ξ)2.
The inclusion of the effective viscosity μe is implicit in (10)–(13), so the shear stress distribution and the effective local viscosity in the turbulent regime is typically found iteratively, such as in the finite element code developed by He [63]. By including the eddy-viscosity model in the TEHD code, convective inertia effects are approximated in the solution. However, temporal inertia effects are not currently considered in the analysis.Equations (10)–(13) are valid for fully developed turbulent flows [69]. To account for low levels of turbulence and transitional flows, Suganami and Szeri proposed an additional factor γ, which modifies (10)) as [1].(14)μe(η,ξ,z)=μ(1+γϵmν).
The factor γ is dependent on the maximum Reynolds' number in the lubricating flow and is defined as.(15)γ={0,Remax⁡<ReL,1-[ReH-Remax⁡ReH-ReL]1/8,ReL≤Remax⁡≤ReH1,ReH<Remax⁡,,
where Remax=(ρuh)/μ is the maximum local Reynolds' number in the lubricant and u is the local fluid velocity, ReL is the critical Reynolds' number for the onset of transition flow from laminar flow and ReH is the critical Reynolds' number for the onset of turbulence. There is not a consensus in the literature on the critical Reynolds' numbers ReL and ReH, corresponding to the onset of transition flow and turbulent flow, respectively [63]. Proposed values of Reynolds' number for the onset of transition flow range from about 500–1000 and range from about 800–1500 for the onset of turbulent flow [1–5].


## 5.4. Perturbed Reynolds Equation

Once the pressure profile solution is found, the generalized Reynolds equation is perturbed. The first-order perturbation results in the equivalent stiffness and damping coefficients,kij and cij, respectively, which have the general form [63–66].(16)kij=-∂f∂uij,cij=-∂f∂u̇ij.
The specific stiffness and damping coefficients are defined by He [63].When (5)–(16) are considered, it is apparent that the turbulence model chosen significantly affects the predicted dynamic coefficients for the tilting pad bearing.


## 6. Averaged Flow Approaches

Another approach to the lubrication problem is the averaged flow method, where the properties of the lubricant across the film are averaged. First proposed by Constantinescu [70] in terms of turbulent mixing length theory, the model was refined in several follow-on papers [29, 71–73].Other authors have used predominantly empirical bulk flow approaches, including Hirs [58] and San Andrés [74–76]. In all of the averaged-flow formulations, the temporal inertia term ρ∂u/∂t is retained from the Navier-Stokes equations and is incorporated into the analysis.The averaged flow approaches rely on averaging of the fluid properties across the film, including average velocity and viscosity. For example, the formulations typically consider an average fluid velocity of the form:(17)1h∫0hudy=uavg.Section6.1 considers the mixing length theory approach. Section 6.2 considers the empirical approaches.


### 6.1. Mixing Length Theory

Constantinescu published a series of papers [29, 70–73] detailing the application of mixing length theory to the turbulent lubrication problem. His focus was on high Reynolds number flows, with Re≥1000. The Reynolds number is defined as Re=(ρωDh)/μ, where h is the local film thickness. The turbulence was modeled by treating all variables associated with the flow as a mean value, defined by (17) plus a fluctuation. Turbulent stress terms resulted from the simplification of the Navier-Stokes equations. To model these stresses, Prandtl's mixing length theory was employed, which produced results that agreed with experimental data on Poiseuille flow. Analytical solutions for average film properties were possible with this assumption. This can be extended for other average film properties, such as viscosity, in THD analyses similar to those proposed by San Andrés [74–76].The series of papers culminated in a journal bearing lubrication theory published by Constantinescu and Galetuse in 1982 [73]. In this paper, a modification to the Reynolds equation was proposed that included the contribution of turbulent stresses. In nondimensional form, the proposed model was(18)∂∂η̂(ĥ3Gη∂p̂∂η̂)+(RL)2∂∂ẑ(ĥ3Gz∂p̂∂ẑ)=12∂ĥ∂η̂-RecrR∇̂(ĥ2Î*).
Implicit in (18) is the shear stress at the boundaries. The term multiplied by Re* in (18) is the contribution of fluid inertia in terms of turbulent stresses. In (18), Gη and Gz are parameters dependent on the average Reynolds number of the flow Re, R is the radius of the bearing, L is the axial length of the bearing, and cr is the radial bearing clearance. The gradient operator ∇̂ is defined as(19)∇̂=i∂∂η̂+kRL∂∂ẑ
and the averaged flow profiles I have the following definitions:(20)Î*=iGηÎη+kGzÎz,Îη=1V2cr(∂Iηη∂η+∂Iηz∂z),Îz=1V2cr(∂Iηz∂η+∂Izz∂z),Iηη=∫0hu2dξ=αUm2h+βV2h-γUmVh,Iηz=∫0huwdξ=α′UmWmh-γ′WmVh,Izz=∫0hw2dξ=α′′Wm2h.
In (20), Um represents the mean fluid velocity in the slider direction, Wm represents the mean fluid velocity in the axial direction, u, v, and w represent local fluid velocities, and the coefficients α, β, γ, and δ and their primes are dependent on the Reynolds number of the flow. The assumed flow profile used in the definitions of the I terms is a parabolic profile which is interpreted as a Poiseuille mean velocity.For laminar flows,Gη=Gz=1/12, α=α′=α′′=6/5, β=δ=2/15, and γ=1/5. For turbulent flows (Re>5000) [73], the coefficients are functions of the mean Reynolds number of the flow:
(21)1Gη=12+0.0136Re0.9,1Gz=12+0.0043Re0.96,α=α′=α′′=1,β=0.885Re0.367,γ=γ′=0,δ=1.95Re0.43.
The Constantinescu approach allows for rapid solutions to the lubrication problem, but there are some drawbacks. The system of equations relies on a minimum of six empirical coefficients to characterize turbulence, which requires extensive experimental data to validate. The approach also relies on an assumption that the average of the product of flow profiles is the product of the averages, that is,(22)∫0hu2dy=hU2,∫0huwdy=hUW,∫0hw2dy=hW2
which cannot be justified as noted by Szeri [77]. Constantinescu acknowledged that (22) was at best only approximately correct [71].


### 6.2. Empirical Approaches

Hirs [58] proposed a model that was predominantly based on experimentally measured bulk-flow properties relative to a surface or wall and the corresponding shear stresses at the boundaries based on a set of flow conditions. The method does not consider the shape of internal flow profiles or fluctuations within the lubricating film.By solely considering the average flow properties and the boundary conditions, Hirs developed a set of pressure equations for sliding surfaces as(23)-h2μU∂p∂η(μρUh)1+m0=12n0{Uη(Uη2+Uz2)(1+m0)/2+(Uη-1)[(Uη-1)2+Uz2](1+m0)/2},-h2μU∂p∂z(μρUh)1+m0=12n0{Uz(Uη2+Uz2)(1+m0)/2+Uz[(Uη-1)2+Uz2](1+m0)/2},
where Uη and Uz represent dimensionless mean flow velocities, and the constants n0 and m0 are found empirically from representative flows. In a subsequent paper, Hirs [78] delineated the various flow regimes requiring experimental data to determine these constants and summarized results for experiments that were already available in the literature. The minimum Reynolds number for any of the empirical coefficients was 1000.In terms of friction factors, the Hirs approach can generally produce more accurate results for the lubrication transition flow regime because experimental data is fitted to the model, assuming such data is available for that bearing configuration and flow condition. The transition region is where many modern oil-lubricated bearings operate. The key drawback is that the method totally relies on empirical data. The types of experiments that would be required to obtain a complete set of empirical coefficients were alluded to by Hirs [78]. However, it is not clear how extensive the experimental support would have to be, especially when factors such as shaft eccentricity ratio, pivot offset, thermal effects and bearing preload are considered. Thermal effects make the estimation of an average Reynolds number difficult. Pivot offset and nonzero bearing preload will alter the shape of the converging wedge. Hirs indicated that high-eccentricity bearings gave less accurate results [58]. The combination of moderate operating eccentricity and high bearing preload may mimic a zero-preload, high-eccentricity bearing in terms of shape of the converging wedge.


### 6.3. Comparison of Approaches

Taylor and Dowson [79] directly compared the mixing length approach of Constantinescu, the eddy-viscosity model of Ng, Pan, and Elrod, and the empirical bulk-flow approach of Hirs. When comparing the predicted factors Gη, all three methods gave close agreement for Reynolds number greater than 2000. The Constantinescu model overpredicted Gη by up to 50 percent when compared to Elrod and Ng and Hirs. All three models deviated from each other in for 1000≤Re≤2000. Based on these results, Taylor and Dowson concluded that the eddy-viscosity model proposed by Elrod and Ng [69] was more accurate than the model proposed by Constantinescu [72].The transition region from laminar flow to turbulent flow presents challenges to the both the eddy-viscosity model and the bulk flow approaches. Suganami and Szeri [1] were able to address this challenge in part with an additional scaling factor in the effective viscosity obtained from the eddy-viscosity model, represented by (11). The Suganami and Szeri scaling factor was able to model temperature rise in a bearing that was run in the transition flow region more accurately than either the laminar or turbulent models. This scaling factor, although empirical, reduces the discrepancy between the Elrod and Ng model and the Hirs model in the transition region.Bouard et al. [80] compared three turbulence models using a finite difference solution to the generalized Reynolds equation, the energy equation, and the heat transfer equation. The three models compared were the Ng and Pan model, [22], the Elrod and Ng model [69], and the Constantinescu model [29]. The comparison was performed within a common finite difference framework, which is distinct from the bulk-flow approach used by Constantinescu. The three models were compared to experiments reported by Taniguchi et al. [4]. All three models overpredicted temperatures in the laminar flow regime, which was attributed to poor characterization of the experimental boundary conditions. All three models matched the experiment within 2 percent at operating speeds above 3,600 rpm. All three models gave similar predictions of power loss. The Constantinescu turbulence model gave predictions closest to measurement for film thickness and babbitt temperature. The authors concluded that all three models gave similar predictions and that Constantinescu should be used because the computational run times were shorter. The TEHD analysis performed by He gave better predictions in the laminar regime [63], using the Elrod and Ng model with the modification proposed by Suganami and Szeri [1] in the transition region.


## 6.1. Mixing Length Theory

Constantinescu published a series of papers [29, 70–73] detailing the application of mixing length theory to the turbulent lubrication problem. His focus was on high Reynolds number flows, with Re≥1000. The Reynolds number is defined as Re=(ρωDh)/μ, where h is the local film thickness. The turbulence was modeled by treating all variables associated with the flow as a mean value, defined by (17) plus a fluctuation. Turbulent stress terms resulted from the simplification of the Navier-Stokes equations. To model these stresses, Prandtl's mixing length theory was employed, which produced results that agreed with experimental data on Poiseuille flow. Analytical solutions for average film properties were possible with this assumption. This can be extended for other average film properties, such as viscosity, in THD analyses similar to those proposed by San Andrés [74–76].The series of papers culminated in a journal bearing lubrication theory published by Constantinescu and Galetuse in 1982 [73]. In this paper, a modification to the Reynolds equation was proposed that included the contribution of turbulent stresses. In nondimensional form, the proposed model was(18)∂∂η̂(ĥ3Gη∂p̂∂η̂)+(RL)2∂∂ẑ(ĥ3Gz∂p̂∂ẑ)=12∂ĥ∂η̂-RecrR∇̂(ĥ2Î*).
Implicit in (18) is the shear stress at the boundaries. The term multiplied by Re* in (18) is the contribution of fluid inertia in terms of turbulent stresses. In (18), Gη and Gz are parameters dependent on the average Reynolds number of the flow Re, R is the radius of the bearing, L is the axial length of the bearing, and cr is the radial bearing clearance. The gradient operator ∇̂ is defined as(19)∇̂=i∂∂η̂+kRL∂∂ẑ
and the averaged flow profiles I have the following definitions:(20)Î*=iGηÎη+kGzÎz,Îη=1V2cr(∂Iηη∂η+∂Iηz∂z),Îz=1V2cr(∂Iηz∂η+∂Izz∂z),Iηη=∫0hu2dξ=αUm2h+βV2h-γUmVh,Iηz=∫0huwdξ=α′UmWmh-γ′WmVh,Izz=∫0hw2dξ=α′′Wm2h.
In (20), Um represents the mean fluid velocity in the slider direction, Wm represents the mean fluid velocity in the axial direction, u, v, and w represent local fluid velocities, and the coefficients α, β, γ, and δ and their primes are dependent on the Reynolds number of the flow. The assumed flow profile used in the definitions of the I terms is a parabolic profile which is interpreted as a Poiseuille mean velocity.For laminar flows,Gη=Gz=1/12, α=α′=α′′=6/5, β=δ=2/15, and γ=1/5. For turbulent flows (Re>5000) [73], the coefficients are functions of the mean Reynolds number of the flow:
(21)1Gη=12+0.0136Re0.9,1Gz=12+0.0043Re0.96,α=α′=α′′=1,β=0.885Re0.367,γ=γ′=0,δ=1.95Re0.43.
The Constantinescu approach allows for rapid solutions to the lubrication problem, but there are some drawbacks. The system of equations relies on a minimum of six empirical coefficients to characterize turbulence, which requires extensive experimental data to validate. The approach also relies on an assumption that the average of the product of flow profiles is the product of the averages, that is,(22)∫0hu2dy=hU2,∫0huwdy=hUW,∫0hw2dy=hW2
which cannot be justified as noted by Szeri [77]. Constantinescu acknowledged that (22) was at best only approximately correct [71].


## 6.2. Empirical Approaches

Hirs [58] proposed a model that was predominantly based on experimentally measured bulk-flow properties relative to a surface or wall and the corresponding shear stresses at the boundaries based on a set of flow conditions. The method does not consider the shape of internal flow profiles or fluctuations within the lubricating film.By solely considering the average flow properties and the boundary conditions, Hirs developed a set of pressure equations for sliding surfaces as(23)-h2μU∂p∂η(μρUh)1+m0=12n0{Uη(Uη2+Uz2)(1+m0)/2+(Uη-1)[(Uη-1)2+Uz2](1+m0)/2},-h2μU∂p∂z(μρUh)1+m0=12n0{Uz(Uη2+Uz2)(1+m0)/2+Uz[(Uη-1)2+Uz2](1+m0)/2},
where Uη and Uz represent dimensionless mean flow velocities, and the constants n0 and m0 are found empirically from representative flows. In a subsequent paper, Hirs [78] delineated the various flow regimes requiring experimental data to determine these constants and summarized results for experiments that were already available in the literature. The minimum Reynolds number for any of the empirical coefficients was 1000.In terms of friction factors, the Hirs approach can generally produce more accurate results for the lubrication transition flow regime because experimental data is fitted to the model, assuming such data is available for that bearing configuration and flow condition. The transition region is where many modern oil-lubricated bearings operate. The key drawback is that the method totally relies on empirical data. The types of experiments that would be required to obtain a complete set of empirical coefficients were alluded to by Hirs [78]. However, it is not clear how extensive the experimental support would have to be, especially when factors such as shaft eccentricity ratio, pivot offset, thermal effects and bearing preload are considered. Thermal effects make the estimation of an average Reynolds number difficult. Pivot offset and nonzero bearing preload will alter the shape of the converging wedge. Hirs indicated that high-eccentricity bearings gave less accurate results [58]. The combination of moderate operating eccentricity and high bearing preload may mimic a zero-preload, high-eccentricity bearing in terms of shape of the converging wedge.


## 6.3. Comparison of Approaches

Taylor and Dowson [79] directly compared the mixing length approach of Constantinescu, the eddy-viscosity model of Ng, Pan, and Elrod, and the empirical bulk-flow approach of Hirs. When comparing the predicted factors Gη, all three methods gave close agreement for Reynolds number greater than 2000. The Constantinescu model overpredicted Gη by up to 50 percent when compared to Elrod and Ng and Hirs. All three models deviated from each other in for 1000≤Re≤2000. Based on these results, Taylor and Dowson concluded that the eddy-viscosity model proposed by Elrod and Ng [69] was more accurate than the model proposed by Constantinescu [72].The transition region from laminar flow to turbulent flow presents challenges to the both the eddy-viscosity model and the bulk flow approaches. Suganami and Szeri [1] were able to address this challenge in part with an additional scaling factor in the effective viscosity obtained from the eddy-viscosity model, represented by (11). The Suganami and Szeri scaling factor was able to model temperature rise in a bearing that was run in the transition flow region more accurately than either the laminar or turbulent models. This scaling factor, although empirical, reduces the discrepancy between the Elrod and Ng model and the Hirs model in the transition region.Bouard et al. [80] compared three turbulence models using a finite difference solution to the generalized Reynolds equation, the energy equation, and the heat transfer equation. The three models compared were the Ng and Pan model, [22], the Elrod and Ng model [69], and the Constantinescu model [29]. The comparison was performed within a common finite difference framework, which is distinct from the bulk-flow approach used by Constantinescu. The three models were compared to experiments reported by Taniguchi et al. [4]. All three models overpredicted temperatures in the laminar flow regime, which was attributed to poor characterization of the experimental boundary conditions. All three models matched the experiment within 2 percent at operating speeds above 3,600 rpm. All three models gave similar predictions of power loss. The Constantinescu turbulence model gave predictions closest to measurement for film thickness and babbitt temperature. The authors concluded that all three models gave similar predictions and that Constantinescu should be used because the computational run times were shorter. The TEHD analysis performed by He gave better predictions in the laminar regime [63], using the Elrod and Ng model with the modification proposed by Suganami and Szeri [1] in the transition region.


## 7. Review of Tilting Pad Bearing Dynamic Models

The tilting pad dynamics developed from various TEHD models are based on explicit modeling of the motion of the pads. The modeling procedure is summarized in the following section. The development of the bearing model reduced to the shaft degrees of freedom and an experimentally identified two-degree-of-freedom bearing model are also summarized.


### 7.1. Single Pad Bearing Dynamics

The lubricating film is typically represented with stiffness and damping coefficients in linear analyses. To illustrate this concept, two free-body diagrams are provided. The first, Figure4, shows rigid shaft interactions with a single tilting pad through the oil film. The springs and dampers in Figure 4 are shown schematically for clarity of the figure. The actual reactions are fluid structure interaction forces between the shaft, the lubricating film, and the pad [77].Figure 4
Free body diagram, dhaft translational degrees of freedom, and rigid pivots.The second free-body diagram, Figure5, shows the linearized fluid-structure interactions between the pad and the shaft, and between the pad and ground. The free body diagrams are shown separately because the linearized stiffness and damping coefficients are in general non-selfadjoint. In Figures 4 and 5, a single pad is shown for clarity of the figures. A typical bearing would have four or five pads.Figure 5
Free body diagram, pad rotational degrees of freedom, and rigid pivots.When a force balance is considered on the free body diagrams, the resulting equations of motion can be expressed in matrix form as [46, 81](24)[Ms000Ms000Jp]{η̈ξ̈θ̈}+[cηηcηξcηθcξηcξξcξθcθηcθξcθθ]{η̇ξ̇θ̇}+[kηηkηξkηθkξηkξξkξθkθηkθξkθθ]{ηξθ}={fηfξ0},
where kij represents lubricant equivalent stiffness, cij represents lubricant equivalent damping, Ms represents the mass of the shaft, Jp represents the mass moment of inertia of the pad about the pivot, η and ξ represent the orthogonal directions in the pad local coordinate system, x and y represent the global coordinate system, and θ represents rotations of the bearing pad about the pivot. For brevity, (24) can be rewritten in matrix form as(25)Mp′v̈p+Cp′v̇p+Kp′vp=fp′.
Primed matrices refer to the pad local (η,ξ) coordinate systems, and unprimed matrices refer to the fixed (x,y) coordinate system. The vector vp represents the rotor and pad motion in pad local coordinates. The tilting pad bearing dynamics are typically transformed to the global shaft coordinate (x,y) system for the purposes of rotordynamic analyses. The coordinate transformation from local to global coordinates for the local pad degrees of freedom is given by [46]:(26){xyθ}=[-sinψ-cos⁡ψ0cos⁡ψ-sinψ0001]{ηξθ},
where ψ is the angle between the x-axis and the pivot location. In matrix form, (26) can be rewritten as up=QTvp, or alternatively Qup=vp, where Q represents the single-pad coordinate transformation matrix and up represents the vector of single pad dynamics expressed in global coordinates, or up=[xyθ]T. The total transformation from local to global coordinates on a per-pad basis is then given by(27)fp=QTfp′=QTMp′Qüp+QTCp′Qu̇p+QTKpQup=Mpüp+Cpu̇p+Kpup.


### 7.2. Assembled Tilting Pad Equation of Motion

Once the individual pad equations of motion are transformed to global coordinates, the overall equations of motion can be assembled [46]. The fundamental equation of motion with no fluid temporal inertia effects, rigid pads, and rigid pivots expressed in terms of shaft degrees of freedom and pad rotations for an Np pad bearing is [46, 82, 83](28)[Ms000⋯00Ms00⋯000J10⋯0000J2⋱⋮⋮⋮⋮⋱⋱0000⋯0JNp]{ẍÿθ̈1θ̈2⋮θ̈Np}+[cxxcxycxθ1cxθ2⋯cxθNpcyxcyycyθ1cyθ2⋯cyθNpcθ1xcθ1ycθ1θ10⋯0cθ2xcθ2y0cθ2θ2⋱⋮⋮⋮⋮⋱⋱0cθNpxcθNpy0⋯0cθNpθNp]{ẋẏθ̇1θ̇2⋮θ̇Np}+[kxxkxykxθ1kxθ2⋯kxθNpkyxkyykyθ1kyθ2⋯kyθNpkθ1xkθ1ykθ1θ10⋯0kθ2xkθ2y0kθ2θ2⋱⋮⋮⋮⋮⋱⋱0kθNpxkθNpy0⋯0kθNpθNp]{xyθ1θ2⋮θNp}={fxfy00⋮0}.
Equation (28) is described explicitly for a five pad tilting pad bearing in Appendix A. The stiffness and damping terms in (28) are also described in terms of the local pad contributions and the coordinate transformations in Appendix A. Equation (28) can also be written in matrix notation as(29)Mü+Cu̇+Ku=f.


### 7.3. Tilting Pad Bearing Model Reduction

Equation (28) has not been traditionally used to describe tilting pad journal bearing behavior in rotordynamic analyses. It is typically reduced dynamically to the shaft degrees of freedom associated with the bearing, and the pad degrees of freedom are not explicitly considered. Historically, there are a few reasons for this choice. Older rotordynamics analyses were based on the transfer matrix method, originally described separately by Myklestad [84] and Prohl [85]. The method, based on beam theory, can be modified to include a discrete stiffness from a bearing as long as it is related to the appropriate beam degree of freedom. The transfer matrix method is not capable of admitting a full-coefficient representation of a tilting pad bearing unless it is transformed into an equivalent transfer function. This exact transfer function was developed by Brockett and Barrett [59]. However, it was published in 1993, when computer power was reaching a point where desktop finite element analyses of rotors were feasible. It was also published after the use of less sophisticated reduced-order bearing models had become the industry standard. Modern rotordynamic analysis packages using finite element formulations such as [53] are capable of using the full coefficient representation.There were also some fundamental misunderstandings of the tilting pad journal bearing results originally presented by Lund in 1964 [20]. The design curves presented by Lund were reduced synchronously, or using the shaft running speed as the reduction frequency. While not intended to be used in general [51], the synchronous coefficients as an excitation-frequency-independent representation of tilting pad bearing dynamics became the industry standard. Manufacturers developed design tools based on the synchronously reduced coefficients for rotating machinery. This has been encoded in industry standards such as API 617 for centrifugal compressors [86]. However, future editions of the API standards will reflect nonsynchronous tilting pad bearing coefficients [87].For the purposes of this discussion, (28) will be reduced to the shaft degrees of freedom. The shaft degree of freedom approach was also heavily influenced by the performance of fixed geometry fluid film bearings, which do not have the additional degrees of freedom associated with the pads. Equation (28) is partitioned into shaft and pad degrees of freedom. It can be rewritten in block matrix format as(30)[Ms00Jp]{üsθ̈}+[CuuCuθCθuCθθ]{u̇sθ̇}+[KuuKuθKθuKθθ]{usθ}={fs0},
where Ms is the diagonal matrix of shaft masses, Jp is the diagonal matrix of pad inertias, Cuu, Cuθ, Cθu, Cθθ represent the damping submatrices, Kuu, Kuθ, Kθu, Kθθ represent the stiffness submatrices, us represents the shaft translation degrees of freedom, θ represents the pad rotation degrees of freedom, and fs represents externally applied forces to the shaft. By expanding (30), the resulting tilting pad bearing model equations take the form:(31)Msüs+Cuuu̇s+Cuθθ̇+Kuuus+Kuθθ=fs,Jpθ̈+Cθuu̇s+Cθθθ̇+Kθuus+Kθθθ=0.Dynamic reduction is performed in the frequency domain by assuming a solution of the formus=Usest,θ=Θest,fs=Fsest. The damped excitation frequency s=p+jq is in general nonsynchronous except in the case of unbalance response. Use of the Laplace transform allows for a general damped solution to the reduction problem [47, 52, 54, 63]. By substituting the assumed solution into (31), the following equations in the frequency domain are found(32)(sCuu+Kuu)Us+(sCuθ+Kuθ)Θ=Fs-s2MsUs,(33)(sCθu+Kθu)Us+(s2Jp+sCθθ+Kθθ)Θ=0.
Next equation (33) is solved in terms of the pad rotations Θ as(34)Θ=-(s2Jp+sCθθ+Kθθ)-1(sCθu+Kθu)Us.
Then by back substitution of (34) into (32), the bearing coefficients are expressed in terms of the shaft degrees of freedom Us as(35)[sCuu+Kuu-(sCuθ+Kuθ)×(s2Jp+sCθθ+Kθθ)-1(sCθu+Kθu)]Us=Fs-s2MsUs.
Using the individual stiffness and damping coefficients from the full tilting pad representation, (35) can be written out in detail as(36)[scxx+kxx-Axxscxy+kxy-Axyscyx+kyx-Ayxscyy+kyy-Ayy]{XY}={Fx-s2MsXFy-s2MsY}.
The four A coefficients are defined for rigid pivots as(37)Axx=∑i=1Np(scxθi+kxθi)(scθix+kθix)s2Jpi+scθiθi+kθiθi,(38)Axy=∑i=1Np(scxθi+kxθi)(scθiy+kθiy)s2Jpi+scθiθi+kθiθi,(39)Ayx=∑i=1Np(scyθi+kyθi)(scθix+kθix)s2Jpi+scθiθi+kθiθi,(40)Ayy=∑i=1Np(scyθi+kyθi)(scθiy+kθiy)s2Jpi+scθiθi+kθiθi.
For a forced response analysis, where s=jΩ, the reduced direct horizontal stiffness as a function of excitation frequency is then(41)k¯xx(Ω)=kxx-Re[∑i=1Np(kxθi+jΩcxθi)(kθix+jΩcθix)kθiθi-Ω2Jpi+jΩcθiθi],
where Re is the real part, and the reduced direct horizontal damping as a function of excitation frequency is(42)c¯xx(Ω)=1ΩIm⁡[jΩcxx-∑i=1Np(kxθi+jΩcxθi)(kθix+jΩcθix)kθiθi-Ω2Jpi+jΩcθiθi],where Im denotes the imaginary part. The reduced stiffness and damping terms k¯xy, k¯yx, k¯yy, c¯xy, c¯yx, c¯yy are found similarly to (41), (42). Returning to the time domain via an inverse Fourier transform, the reduced model of bearing dynamics then becomes:(43)[Ms00Ms]{ẍÿ}+[c¯xx(Ω)c¯xy(Ω)c¯yx(Ω)c¯yy(Ω)]{ẋẏ}+[k¯xx(Ω)k¯xy(Ω)k¯yx(Ω)k¯yy(Ω)]{xy}={fxfy}.
For a free response analysis where fx=fy=0, the effective stiffness and damping can be found numerically in the form:(44)k¯xx(s)=Re(kxx+scxx-∑i=1Np(scxθi+kxθi)(scθix+kθix)s2Jpi+scθiθi+kθiθi),c¯xx(s)=Im⁡[1s(∑i=1Np(scxθi+kxθi)(scθix+kθix)s2Jpi+scθiθi+kθiθikxx+scxx-∑i=1Np(scxθi+kxθi)(scθix+kθix)s2Jpi+scθiθi+kθiθi)].If the perturbation frequencys=jΩ is taken to be the rotational speed of the machine, then (41), (42) represent the direct horizontal synchronously reduced stiffness and damping coefficients. The synchronously reduced k¯xy, k¯yx, k¯yy, c¯xy, c¯yx, c¯yy are found similarly. If the eigenvalue s=p+jq is treated as a general perturbation frequency not equal to the machine rotational speed, then the coefficients are the nonsynchronously reduced coefficients at the excitation frequency of interest. The resulting coefficient matrices are non-self-adjoint. The destabilizing tangential forces due to fluid structure interactions are represented as cross-coupled stiffnesses k¯xy and k¯yx. However, these cross-coupled stiffness terms are generally 3 orders of magnitude less than the direct stiffness terms in tilting pad bearings and typically neglected. In the frequency domain, the effective tangential forces due to damping, proportional to Ωc¯xx,Ωc¯yy, are also much greater than the cross-coupled stiffness terms, which indicates that the destabilizing forces are small and do not adversely affect bearing dynamic performance.The reduced-order model with pad dynamics considered implicitly, (43), is applicable to rigid pivot bearings. Representation of pivot flexibility in tilting pad bearings requires consideration of additional degrees of freedom. Treatments of the pivot flexibility case have been addressed by several authors, including [43, 50].


### 7.4. Reduced Order Nonsynchronous Bearing Models

An alternative experimental approach to characterizing TPJB behavior is based on an experimentally identified model in the frequency domain. The experimentally derived model is based on measurement of force inputs and bearing housing outputs. The shaft is held rigidly in rolling element bearings, and the bearing is allowed to move radially. The bearing housing is perturbed and displacements of the bearing relative to the shaft are measured. The method, originally applied to fixed pad hydrostatic bearings [88, 89], has been applied recently to flexible pivot bearings [90, 91] and four-pad and five-pad tilting pad journal bearings [92–95]. A brief description of the experimental identification method follows. A detailed description is available in the above references.The system identification method employed in [92–95] assumes that the bearing dynamic properties can be modeled as a two degree-of-freedom system based on tilting pad bearing housing (x,y) translations. The identification procedure described in [88–95] is accomplished in the frequency domain by applying a sinusoidal excitation of the form x=XejΩt,y=YejΩt,fx=FxejΩt,fy=FyejΩt, where the excitation frequency Ω is in general nonsynchronous. Excitation is accomplished by simultaneously applying several sinusoidal forces to the bearing housing with mass M, resulting in a pseudorandom perturbation [88, 89]. The resulting complex impedance is then determined as a function of excitation frequency. In the frequency domain, the net bearing response is expressed in terms of complex impedances Zij in the form:(45)[ZxxZxyZyxZyy]{XY}={Fx+Ω2MXFy+Ω2MY}.The real and imaginary parts of the complex impedance functionszij are then plotted as a function of frequency. The method is capable of discerning frequency dependence in both real and imaginary parts of the complex impedances, within the limits imposed by measurement uncertainty and repeatability. The subscripts i,j represent the appropriate rotor degree of freedom. Power spectral density functions are used to reduce the effects of noise on the measurements. For the TPJB bearings reported in [92–94], the trends in the data resulted in an adequate (r2≥0.95) model described by(46)Re(Zij)=k̃ij-Ω2m̃ij;Im⁡(Zij)=Ωc̃ij.
The complex impedance of the system then takes the form:(47)Zij=k̃ij-Ω2m̃ij+jΩc̃ij.
Then, by substituting (47) into (45), and performing an inverse Fourier transform to return (45) to the time domain, the resulting model for TPJB behavior is given by(48)[M+m̃xxm̃xym̃yxM+m̃yy]{ẍÿ}+[c̃xxc̃xyc̃yxc̃yy]{ẋẏ}+[k̃xxk̃xyk̃yxk̃yy]{xy}={fxfy},
where the m̃ij represent the identified lubricant mass coefficients, the c̃ij represent the identified damping coefficients, and the k̃ij represent the identified stiffness coefficients. Equation (48) will be referred to as the KCM model. The second-order representation has been proposed as a nonsynchronous representation of TPJB behavior with twelve frequency-independent dynamic coefficients and two degrees of freedom [92–94], which is in contrast with the frequency-dependent KC model presented in Section 7.3. The effect of pad dynamics on the overall tilting pad bearing dynamics is not considered explicitly in this formulation.These results were reviewed by Childs [96]. In this paper, the results from several bearing tests were reviewed, and the author stated that there was no apparent frequency dependency other than that captured by the KCM model, with the exception of the bearing originally reported in Childs and Harris [94], where an apparent frequency dependence in the damping coefficient was observed. However, this paper did not address the negative identified lubricant inertia coefficients reported for several tests in [92–94].When the frequency response data for flexible pivot bearings reported in [90, 91] and tilting pad bearings [92–94] is interpreted using (48) by the respective authors, several common themes emerge. The reported data are compared to two models: a tilting pad bearing model based solely on the classic Reynolds equation, (1), and a thermohydrodynamic bulk-flow analysis developed by San Andrés and presented in [74–76]. The model developed by San Andrés is based on theories developed by Hirs [58] and Constantinescu [29, 70] that include temporal and convective inertia terms from the Navier-Stokes equations averaged across the lubricating film. The quadratic behavior observed in the real part of the impedance occurs with both the Reynolds equation and the bulk flow models. The data generally have better agreement with the bulk-flow model than with the Reynolds equation model, for example, [90], so the improved agreement is attributed to both temporal and convective fluid inertia effects.


## 7.1. Single Pad Bearing Dynamics

The lubricating film is typically represented with stiffness and damping coefficients in linear analyses. To illustrate this concept, two free-body diagrams are provided. The first, Figure4, shows rigid shaft interactions with a single tilting pad through the oil film. The springs and dampers in Figure 4 are shown schematically for clarity of the figure. The actual reactions are fluid structure interaction forces between the shaft, the lubricating film, and the pad [77].Figure 4
Free body diagram, dhaft translational degrees of freedom, and rigid pivots.The second free-body diagram, Figure5, shows the linearized fluid-structure interactions between the pad and the shaft, and between the pad and ground. The free body diagrams are shown separately because the linearized stiffness and damping coefficients are in general non-selfadjoint. In Figures 4 and 5, a single pad is shown for clarity of the figures. A typical bearing would have four or five pads.Figure 5
Free body diagram, pad rotational degrees of freedom, and rigid pivots.When a force balance is considered on the free body diagrams, the resulting equations of motion can be expressed in matrix form as [46, 81](24)[Ms000Ms000Jp]{η̈ξ̈θ̈}+[cηηcηξcηθcξηcξξcξθcθηcθξcθθ]{η̇ξ̇θ̇}+[kηηkηξkηθkξηkξξkξθkθηkθξkθθ]{ηξθ}={fηfξ0},
where kij represents lubricant equivalent stiffness, cij represents lubricant equivalent damping, Ms represents the mass of the shaft, Jp represents the mass moment of inertia of the pad about the pivot, η and ξ represent the orthogonal directions in the pad local coordinate system, x and y represent the global coordinate system, and θ represents rotations of the bearing pad about the pivot. For brevity, (24) can be rewritten in matrix form as(25)Mp′v̈p+Cp′v̇p+Kp′vp=fp′.
Primed matrices refer to the pad local (η,ξ) coordinate systems, and unprimed matrices refer to the fixed (x,y) coordinate system. The vector vp represents the rotor and pad motion in pad local coordinates. The tilting pad bearing dynamics are typically transformed to the global shaft coordinate (x,y) system for the purposes of rotordynamic analyses. The coordinate transformation from local to global coordinates for the local pad degrees of freedom is given by [46]:(26){xyθ}=[-sinψ-cos⁡ψ0cos⁡ψ-sinψ0001]{ηξθ},
where ψ is the angle between the x-axis and the pivot location. In matrix form, (26) can be rewritten as up=QTvp, or alternatively Qup=vp, where Q represents the single-pad coordinate transformation matrix and up represents the vector of single pad dynamics expressed in global coordinates, or up=[xyθ]T. The total transformation from local to global coordinates on a per-pad basis is then given by(27)fp=QTfp′=QTMp′Qüp+QTCp′Qu̇p+QTKpQup=Mpüp+Cpu̇p+Kpup.


## 7.2. Assembled Tilting Pad Equation of Motion

Once the individual pad equations of motion are transformed to global coordinates, the overall equations of motion can be assembled [46]. The fundamental equation of motion with no fluid temporal inertia effects, rigid pads, and rigid pivots expressed in terms of shaft degrees of freedom and pad rotations for an Np pad bearing is [46, 82, 83](28)[Ms000⋯00Ms00⋯000J10⋯0000J2⋱⋮⋮⋮⋮⋱⋱0000⋯0JNp]{ẍÿθ̈1θ̈2⋮θ̈Np}+[cxxcxycxθ1cxθ2⋯cxθNpcyxcyycyθ1cyθ2⋯cyθNpcθ1xcθ1ycθ1θ10⋯0cθ2xcθ2y0cθ2θ2⋱⋮⋮⋮⋮⋱⋱0cθNpxcθNpy0⋯0cθNpθNp]{ẋẏθ̇1θ̇2⋮θ̇Np}+[kxxkxykxθ1kxθ2⋯kxθNpkyxkyykyθ1kyθ2⋯kyθNpkθ1xkθ1ykθ1θ10⋯0kθ2xkθ2y0kθ2θ2⋱⋮⋮⋮⋮⋱⋱0kθNpxkθNpy0⋯0kθNpθNp]{xyθ1θ2⋮θNp}={fxfy00⋮0}.
Equation (28) is described explicitly for a five pad tilting pad bearing in Appendix A. The stiffness and damping terms in (28) are also described in terms of the local pad contributions and the coordinate transformations in Appendix A. Equation (28) can also be written in matrix notation as(29)Mü+Cu̇+Ku=f.


## 7.3. Tilting Pad Bearing Model Reduction

Equation (28) has not been traditionally used to describe tilting pad journal bearing behavior in rotordynamic analyses. It is typically reduced dynamically to the shaft degrees of freedom associated with the bearing, and the pad degrees of freedom are not explicitly considered. Historically, there are a few reasons for this choice. Older rotordynamics analyses were based on the transfer matrix method, originally described separately by Myklestad [84] and Prohl [85]. The method, based on beam theory, can be modified to include a discrete stiffness from a bearing as long as it is related to the appropriate beam degree of freedom. The transfer matrix method is not capable of admitting a full-coefficient representation of a tilting pad bearing unless it is transformed into an equivalent transfer function. This exact transfer function was developed by Brockett and Barrett [59]. However, it was published in 1993, when computer power was reaching a point where desktop finite element analyses of rotors were feasible. It was also published after the use of less sophisticated reduced-order bearing models had become the industry standard. Modern rotordynamic analysis packages using finite element formulations such as [53] are capable of using the full coefficient representation.There were also some fundamental misunderstandings of the tilting pad journal bearing results originally presented by Lund in 1964 [20]. The design curves presented by Lund were reduced synchronously, or using the shaft running speed as the reduction frequency. While not intended to be used in general [51], the synchronous coefficients as an excitation-frequency-independent representation of tilting pad bearing dynamics became the industry standard. Manufacturers developed design tools based on the synchronously reduced coefficients for rotating machinery. This has been encoded in industry standards such as API 617 for centrifugal compressors [86]. However, future editions of the API standards will reflect nonsynchronous tilting pad bearing coefficients [87].For the purposes of this discussion, (28) will be reduced to the shaft degrees of freedom. The shaft degree of freedom approach was also heavily influenced by the performance of fixed geometry fluid film bearings, which do not have the additional degrees of freedom associated with the pads. Equation (28) is partitioned into shaft and pad degrees of freedom. It can be rewritten in block matrix format as(30)[Ms00Jp]{üsθ̈}+[CuuCuθCθuCθθ]{u̇sθ̇}+[KuuKuθKθuKθθ]{usθ}={fs0},
where Ms is the diagonal matrix of shaft masses, Jp is the diagonal matrix of pad inertias, Cuu, Cuθ, Cθu, Cθθ represent the damping submatrices, Kuu, Kuθ, Kθu, Kθθ represent the stiffness submatrices, us represents the shaft translation degrees of freedom, θ represents the pad rotation degrees of freedom, and fs represents externally applied forces to the shaft. By expanding (30), the resulting tilting pad bearing model equations take the form:(31)Msüs+Cuuu̇s+Cuθθ̇+Kuuus+Kuθθ=fs,Jpθ̈+Cθuu̇s+Cθθθ̇+Kθuus+Kθθθ=0.Dynamic reduction is performed in the frequency domain by assuming a solution of the formus=Usest,θ=Θest,fs=Fsest. The damped excitation frequency s=p+jq is in general nonsynchronous except in the case of unbalance response. Use of the Laplace transform allows for a general damped solution to the reduction problem [47, 52, 54, 63]. By substituting the assumed solution into (31), the following equations in the frequency domain are found(32)(sCuu+Kuu)Us+(sCuθ+Kuθ)Θ=Fs-s2MsUs,(33)(sCθu+Kθu)Us+(s2Jp+sCθθ+Kθθ)Θ=0.
Next equation (33) is solved in terms of the pad rotations Θ as(34)Θ=-(s2Jp+sCθθ+Kθθ)-1(sCθu+Kθu)Us.
Then by back substitution of (34) into (32), the bearing coefficients are expressed in terms of the shaft degrees of freedom Us as(35)[sCuu+Kuu-(sCuθ+Kuθ)×(s2Jp+sCθθ+Kθθ)-1(sCθu+Kθu)]Us=Fs-s2MsUs.
Using the individual stiffness and damping coefficients from the full tilting pad representation, (35) can be written out in detail as(36)[scxx+kxx-Axxscxy+kxy-Axyscyx+kyx-Ayxscyy+kyy-Ayy]{XY}={Fx-s2MsXFy-s2MsY}.
The four A coefficients are defined for rigid pivots as(37)Axx=∑i=1Np(scxθi+kxθi)(scθix+kθix)s2Jpi+scθiθi+kθiθi,(38)Axy=∑i=1Np(scxθi+kxθi)(scθiy+kθiy)s2Jpi+scθiθi+kθiθi,(39)Ayx=∑i=1Np(scyθi+kyθi)(scθix+kθix)s2Jpi+scθiθi+kθiθi,(40)Ayy=∑i=1Np(scyθi+kyθi)(scθiy+kθiy)s2Jpi+scθiθi+kθiθi.
For a forced response analysis, where s=jΩ, the reduced direct horizontal stiffness as a function of excitation frequency is then(41)k¯xx(Ω)=kxx-Re[∑i=1Np(kxθi+jΩcxθi)(kθix+jΩcθix)kθiθi-Ω2Jpi+jΩcθiθi],
where Re is the real part, and the reduced direct horizontal damping as a function of excitation frequency is(42)c¯xx(Ω)=1ΩIm⁡[jΩcxx-∑i=1Np(kxθi+jΩcxθi)(kθix+jΩcθix)kθiθi-Ω2Jpi+jΩcθiθi],where Im denotes the imaginary part. The reduced stiffness and damping terms k¯xy, k¯yx, k¯yy, c¯xy, c¯yx, c¯yy are found similarly to (41), (42). Returning to the time domain via an inverse Fourier transform, the reduced model of bearing dynamics then becomes:(43)[Ms00Ms]{ẍÿ}+[c¯xx(Ω)c¯xy(Ω)c¯yx(Ω)c¯yy(Ω)]{ẋẏ}+[k¯xx(Ω)k¯xy(Ω)k¯yx(Ω)k¯yy(Ω)]{xy}={fxfy}.
For a free response analysis where fx=fy=0, the effective stiffness and damping can be found numerically in the form:(44)k¯xx(s)=Re(kxx+scxx-∑i=1Np(scxθi+kxθi)(scθix+kθix)s2Jpi+scθiθi+kθiθi),c¯xx(s)=Im⁡[1s(∑i=1Np(scxθi+kxθi)(scθix+kθix)s2Jpi+scθiθi+kθiθikxx+scxx-∑i=1Np(scxθi+kxθi)(scθix+kθix)s2Jpi+scθiθi+kθiθi)].If the perturbation frequencys=jΩ is taken to be the rotational speed of the machine, then (41), (42) represent the direct horizontal synchronously reduced stiffness and damping coefficients. The synchronously reduced k¯xy, k¯yx, k¯yy, c¯xy, c¯yx, c¯yy are found similarly. If the eigenvalue s=p+jq is treated as a general perturbation frequency not equal to the machine rotational speed, then the coefficients are the nonsynchronously reduced coefficients at the excitation frequency of interest. The resulting coefficient matrices are non-self-adjoint. The destabilizing tangential forces due to fluid structure interactions are represented as cross-coupled stiffnesses k¯xy and k¯yx. However, these cross-coupled stiffness terms are generally 3 orders of magnitude less than the direct stiffness terms in tilting pad bearings and typically neglected. In the frequency domain, the effective tangential forces due to damping, proportional to Ωc¯xx,Ωc¯yy, are also much greater than the cross-coupled stiffness terms, which indicates that the destabilizing forces are small and do not adversely affect bearing dynamic performance.The reduced-order model with pad dynamics considered implicitly, (43), is applicable to rigid pivot bearings. Representation of pivot flexibility in tilting pad bearings requires consideration of additional degrees of freedom. Treatments of the pivot flexibility case have been addressed by several authors, including [43, 50].


## 7.4. Reduced Order Nonsynchronous Bearing Models

An alternative experimental approach to characterizing TPJB behavior is based on an experimentally identified model in the frequency domain. The experimentally derived model is based on measurement of force inputs and bearing housing outputs. The shaft is held rigidly in rolling element bearings, and the bearing is allowed to move radially. The bearing housing is perturbed and displacements of the bearing relative to the shaft are measured. The method, originally applied to fixed pad hydrostatic bearings [88, 89], has been applied recently to flexible pivot bearings [90, 91] and four-pad and five-pad tilting pad journal bearings [92–95]. A brief description of the experimental identification method follows. A detailed description is available in the above references.The system identification method employed in [92–95] assumes that the bearing dynamic properties can be modeled as a two degree-of-freedom system based on tilting pad bearing housing (x,y) translations. The identification procedure described in [88–95] is accomplished in the frequency domain by applying a sinusoidal excitation of the form x=XejΩt,y=YejΩt,fx=FxejΩt,fy=FyejΩt, where the excitation frequency Ω is in general nonsynchronous. Excitation is accomplished by simultaneously applying several sinusoidal forces to the bearing housing with mass M, resulting in a pseudorandom perturbation [88, 89]. The resulting complex impedance is then determined as a function of excitation frequency. In the frequency domain, the net bearing response is expressed in terms of complex impedances Zij in the form:(45)[ZxxZxyZyxZyy]{XY}={Fx+Ω2MXFy+Ω2MY}.The real and imaginary parts of the complex impedance functionszij are then plotted as a function of frequency. The method is capable of discerning frequency dependence in both real and imaginary parts of the complex impedances, within the limits imposed by measurement uncertainty and repeatability. The subscripts i,j represent the appropriate rotor degree of freedom. Power spectral density functions are used to reduce the effects of noise on the measurements. For the TPJB bearings reported in [92–94], the trends in the data resulted in an adequate (r2≥0.95) model described by(46)Re(Zij)=k̃ij-Ω2m̃ij;Im⁡(Zij)=Ωc̃ij.
The complex impedance of the system then takes the form:(47)Zij=k̃ij-Ω2m̃ij+jΩc̃ij.
Then, by substituting (47) into (45), and performing an inverse Fourier transform to return (45) to the time domain, the resulting model for TPJB behavior is given by(48)[M+m̃xxm̃xym̃yxM+m̃yy]{ẍÿ}+[c̃xxc̃xyc̃yxc̃yy]{ẋẏ}+[k̃xxk̃xyk̃yxk̃yy]{xy}={fxfy},
where the m̃ij represent the identified lubricant mass coefficients, the c̃ij represent the identified damping coefficients, and the k̃ij represent the identified stiffness coefficients. Equation (48) will be referred to as the KCM model. The second-order representation has been proposed as a nonsynchronous representation of TPJB behavior with twelve frequency-independent dynamic coefficients and two degrees of freedom [92–94], which is in contrast with the frequency-dependent KC model presented in Section 7.3. The effect of pad dynamics on the overall tilting pad bearing dynamics is not considered explicitly in this formulation.These results were reviewed by Childs [96]. In this paper, the results from several bearing tests were reviewed, and the author stated that there was no apparent frequency dependency other than that captured by the KCM model, with the exception of the bearing originally reported in Childs and Harris [94], where an apparent frequency dependence in the damping coefficient was observed. However, this paper did not address the negative identified lubricant inertia coefficients reported for several tests in [92–94].When the frequency response data for flexible pivot bearings reported in [90, 91] and tilting pad bearings [92–94] is interpreted using (48) by the respective authors, several common themes emerge. The reported data are compared to two models: a tilting pad bearing model based solely on the classic Reynolds equation, (1), and a thermohydrodynamic bulk-flow analysis developed by San Andrés and presented in [74–76]. The model developed by San Andrés is based on theories developed by Hirs [58] and Constantinescu [29, 70] that include temporal and convective inertia terms from the Navier-Stokes equations averaged across the lubricating film. The quadratic behavior observed in the real part of the impedance occurs with both the Reynolds equation and the bulk flow models. The data generally have better agreement with the bulk-flow model than with the Reynolds equation model, for example, [90], so the improved agreement is attributed to both temporal and convective fluid inertia effects.


## 8. Temporal Inertia Effects

The inclusion of temporal inertia effects in the bulk flow model is justified in [90–94] using the work of Reinhardt and Lund [97]. Reinhardt and Lund retained the temporal and convective inertia terms in their nondimensional formulation of the lubrication problem and obtained a solution that indicated that temporal inertia effects were important in lubricating flows with Re≥100. This is in contrast to results presented by Szeri et al. [98] for squeeze film dampers and Szeri [77] for fluid film bearings. A detailed discussion of the difference in the two models and the underlying physics is presented in [83] and is repeated in Section 8.1.The theory proposed by San Andrés is in contrast to the TEHD analysis based on the generalized Reynolds equation proposed by He [63] and summarized in Section 5. The generalized Reynolds equation presented by He accounts for convective inertia effects, but not temporal inertia effects, through the use of an eddy viscosity model. The eddy viscosity model used in the TEHD analysis by He [63] has an effect on both the fluid effective stiffness and damping, but the TEHD theory does not predict inertia coefficients.A discussion of the effect on stability analysis is also presented in [90–94]. All of the papers discuss the use of synchronous versus nonsynchronous coefficients in stability analyses and state that a frequency-dependent KC model requires an iterative solution to calculate the system eigenvalues. This is correct for older rotordynamic analyses, especially transfer matrix analyses. However, modern rotordynamic codes such as the one documented by [53] are finite element based and can easily accept the additional degrees of freedom required to represent pad motion for implementation of the full KC TPJB model. Since the pad degrees of freedom are explicit within this framework, an iterative eigenvalue solution is not required. It has been recently shown that the KCM model is not guaranteed to produce a conservative estimate of flexible rotor stability [99].


### 8.1. Comparison of Reinhardt and Lund to Szeri

There have been two distinct approaches to temporal inertia effects in hydrodynamic lubrication documented in the literature. The two approaches were compared originally in [83], and the discussion is repeated here for completeness.Both Reinhardt and Lund [97] and Szeri et al. [98] considered the effects of fluid inertia on journal bearing lubricating flows by investigation of the convective and temporal inertia terms in the Navier-Stokes equations. Szeri et al. [98] considered the analysis for squeeze film dampers, but expanded the analysis to journal bearings in [77].Both analyses agree on the nondimensional form of the perturbed N-S equations with inertia terms. The general approach to calculate rotordynamic coefficients is to perform a Taylor series expansion about the reduced Reynolds number, resulting in anx-direction force:(49)f̂x=f̂x0(0)+Re*f̂x0(1)+(k̂xx(0)+Re*k̂xx(1))Δx̂+(k̂xy(0)+Re*k̂xy(1))Δŷ+(ĉxx(0)+Re*ĉxx(1))Δx̂̇+(ĉxy(0)+Re*ĉxy(1))Δŷ̇+Re*(m̂xxΔx̂̈+m̂xyΔŷ̈)+O[(Re*)2].
The hat symbol in (49) indicates nondimensional quantities. The y-direction force is similar.Both analyses also essentially agree on the nondimensionalization of force, damping coefficients, and stiffness coefficients, with minor differences in expression of rotational speed:(50)f̂i=fi(μωLR2/πcr2),k̂ij=kij[(μLω/π)(R/cr)3],ĉij=cij[(μL/π)(R/cr)3].
There is a key difference in the two approaches to the inertia terms at this point. Reinhardt and Lund [97] considered the fluid density when nondimensionalizing the inertia terms, resulting in(51)m̂ij=mij[(ρπR2L/π2)(R/cr)].
When using (51) to dimensionalize results, Reinhardt and Lund derived a dimensional expansion with no influence of reduced Reynolds number on the inertia terms, resulting in [97].(52)fx=fx0(0)+Re*Fx0(1)+(kxx(0)+Re*kxx(1))Δx+(kxy(0)+Re*kxy(1))Δy+(cxx(0)+Re*cxx(1))Δẋ+(cxy(0)+Re*cxy(1))Δẏ+mxxΔẍ+mxyΔÿ.
However, Szeri et al. [98] and Szeri [77] considered the fluid viscosity when nondimensionalizing the inertia terms, resulting in:(53)m̂ij=mij[(μL/πω)(R/cr)3].
This choice of nondimensionalization results in a dimensional expansion that indicates that reduced Reynolds number is a coefficient on the dimensional inertia terms, or(54)fx=fx0(0)+Re*fx0(1)+(kxx(0)+Re*kxx(1))Δx+(kxy(0)+Re*kxy(1))Δy+(cxx(0)+Re*cxx(1))Δẋ+(cxy(0)+Re*cxy(1))Δẏ+Re*(mxxΔẍ+mxyΔÿ).When considering a laminar lubricating type flow, the flow results are dominated by fluid shear effects. As a result, the fluid viscosity is more fundamental than fluid density in nondimensionalizing and scaling the results. This would imply that for low reduced Reynolds number, the fluid inertia (added mass) effects are not significant. Using data derived in [97] for a plain journal bearing with a typical diameter of 127 mm and Re*=0.077, (54) results in an added mass coefficient of 1 kg, versus an added mass coefficient of 12 kg using (52) [83].The Taylor series expansion resulting in (50) is also based on Re*. For the Taylor series expansion to converge, Re*<1 [77]. Modern bearings are reaching surface speeds where Re*≥1, so the Reinhardt and Lund and Szeri analyses of the relative effect of temporal inertia on fluid film bearings are no longer valid as the series does not converge. These analyses are also not valid for low-viscosity process fluid lubricants such as water, where Re*≫1 for typical industrial applications. An alternative method will have to be found for linearized nondimensional estimates of turbulent, inertial flow.Temporal inertia terms will be important for bearings with low viscosity lubricants or highRe* oil bearings so temporal inertia terms will be present. However, those inertia terms are not very important in small amplitude linear modeling of rotor dynamic calculations unless there is some sort of significant radial acceleration in the shaft. There will have to be major research done as to whether an extension of Reynolds equation following Elrod and Ng, averaged film relations following Hirs and Constantinescu, or new forms of lubrication modeling will produce the most accurate predictions of bearing performance.


## 8.1. Comparison of Reinhardt and Lund to Szeri

There have been two distinct approaches to temporal inertia effects in hydrodynamic lubrication documented in the literature. The two approaches were compared originally in [83], and the discussion is repeated here for completeness.Both Reinhardt and Lund [97] and Szeri et al. [98] considered the effects of fluid inertia on journal bearing lubricating flows by investigation of the convective and temporal inertia terms in the Navier-Stokes equations. Szeri et al. [98] considered the analysis for squeeze film dampers, but expanded the analysis to journal bearings in [77].Both analyses agree on the nondimensional form of the perturbed N-S equations with inertia terms. The general approach to calculate rotordynamic coefficients is to perform a Taylor series expansion about the reduced Reynolds number, resulting in anx-direction force:(49)f̂x=f̂x0(0)+Re*f̂x0(1)+(k̂xx(0)+Re*k̂xx(1))Δx̂+(k̂xy(0)+Re*k̂xy(1))Δŷ+(ĉxx(0)+Re*ĉxx(1))Δx̂̇+(ĉxy(0)+Re*ĉxy(1))Δŷ̇+Re*(m̂xxΔx̂̈+m̂xyΔŷ̈)+O[(Re*)2].
The hat symbol in (49) indicates nondimensional quantities. The y-direction force is similar.Both analyses also essentially agree on the nondimensionalization of force, damping coefficients, and stiffness coefficients, with minor differences in expression of rotational speed:(50)f̂i=fi(μωLR2/πcr2),k̂ij=kij[(μLω/π)(R/cr)3],ĉij=cij[(μL/π)(R/cr)3].
There is a key difference in the two approaches to the inertia terms at this point. Reinhardt and Lund [97] considered the fluid density when nondimensionalizing the inertia terms, resulting in(51)m̂ij=mij[(ρπR2L/π2)(R/cr)].
When using (51) to dimensionalize results, Reinhardt and Lund derived a dimensional expansion with no influence of reduced Reynolds number on the inertia terms, resulting in [97].(52)fx=fx0(0)+Re*Fx0(1)+(kxx(0)+Re*kxx(1))Δx+(kxy(0)+Re*kxy(1))Δy+(cxx(0)+Re*cxx(1))Δẋ+(cxy(0)+Re*cxy(1))Δẏ+mxxΔẍ+mxyΔÿ.
However, Szeri et al. [98] and Szeri [77] considered the fluid viscosity when nondimensionalizing the inertia terms, resulting in:(53)m̂ij=mij[(μL/πω)(R/cr)3].
This choice of nondimensionalization results in a dimensional expansion that indicates that reduced Reynolds number is a coefficient on the dimensional inertia terms, or(54)fx=fx0(0)+Re*fx0(1)+(kxx(0)+Re*kxx(1))Δx+(kxy(0)+Re*kxy(1))Δy+(cxx(0)+Re*cxx(1))Δẋ+(cxy(0)+Re*cxy(1))Δẏ+Re*(mxxΔẍ+mxyΔÿ).When considering a laminar lubricating type flow, the flow results are dominated by fluid shear effects. As a result, the fluid viscosity is more fundamental than fluid density in nondimensionalizing and scaling the results. This would imply that for low reduced Reynolds number, the fluid inertia (added mass) effects are not significant. Using data derived in [97] for a plain journal bearing with a typical diameter of 127 mm and Re*=0.077, (54) results in an added mass coefficient of 1 kg, versus an added mass coefficient of 12 kg using (52) [83].The Taylor series expansion resulting in (50) is also based on Re*. For the Taylor series expansion to converge, Re*<1 [77]. Modern bearings are reaching surface speeds where Re*≥1, so the Reinhardt and Lund and Szeri analyses of the relative effect of temporal inertia on fluid film bearings are no longer valid as the series does not converge. These analyses are also not valid for low-viscosity process fluid lubricants such as water, where Re*≫1 for typical industrial applications. An alternative method will have to be found for linearized nondimensional estimates of turbulent, inertial flow.Temporal inertia terms will be important for bearings with low viscosity lubricants or highRe* oil bearings so temporal inertia terms will be present. However, those inertia terms are not very important in small amplitude linear modeling of rotor dynamic calculations unless there is some sort of significant radial acceleration in the shaft. There will have to be major research done as to whether an extension of Reynolds equation following Elrod and Ng, averaged film relations following Hirs and Constantinescu, or new forms of lubrication modeling will produce the most accurate predictions of bearing performance.


## 9. Discussion and Conclusions

Since the original development of the lubrication equation by Reynolds [7], there has been an increasing level of sophistication in the calculation of bearing properties, due to the inclusion of thermal heating effects, mechanical and thermal deformations, and turbulence corrections. The initial solutions only considered the fluid flow inside the bearing. The thermal effects were added through solutions to the energy equation and mechanical deformations were included with deflection analyses. One of the key differences between the TEHD model presented in Section 5 and the bulk flow models presented in Section 6 is the treatment of the lubricating film. The TEHD analysis proceeds from a differential approach to the flow field. The mathematics involved are more easily justified than the averaging approaches employed in mixing length theory model of Constantinescu, which rely on an approximation of the products of averages of the flow field. The Hirs model gives good agreement with friction data because it is fit to that data—hence it relies entirely on extensive empirical data to be implemented. As a result, the approach requires experimental data for each new application.Generally, consideration of more complex lubrication models followed the experience of industrial users. As the classical Reynolds solution diverged from user experience, the addition of energy and deformation effects into the analysis became necessary to allow for reliable designs. The factors requiring these modeling improvements, including increasing speeds and bearing specific loads, are demands by industrial users that continue to influence the need to improve tilting pad bearing models.There has been a similar evolution in the understanding of the bearing dynamics, especially with tilting pad journal bearings. Initially treated as simple supports, inclusion of stiffness effects and later damping effects improved the understanding of the bearing contribution to the overall rotordynamic system. These improvements came as user experience did not match with simpler bearing dynamic models.Tilting pad bearings were adopted to address self-excited vibrations from the fluid structure interactions within fixed pad bearings. The initial understanding was that the synchronous response was a sufficient representation of the bearing dynamics regardless of excitation frequency based on a misinterpretation of the work by Lund. More recent investigations, especially into rotordynamic stability, indicate that the dynamic response is excitation-frequency-dependent.The nonsynchronous modeling presented in Sections7.1–7.3 is not comprehensive since pivot flexibility and foundation flexibility effects are not considered. The discussion does cover the basic ideas in the current understanding of bearing dynamic theory.The KCM experimentally identified model is a fundamentally different model compared to the full bearing coefficient model. The full bearing coefficients are obtained from first principles. The KCM model is based solely on system identification experiments and arises from observation of the system. The observations are consistent with a 12-coefficient second-order nonsynchronous dynamic representation with frequency-independent stiffness, damping, and mass coefficients. This is a “black box” identification technique this is suitable for obtaining a tentative system model and is a technique that is also popular in the controls community for developing an approximate model of a plant to be controlled.The issue of the relative importance of temporal inertia effects is still an open area of discussion. There are conflicting treatments in the literature of the relative importance of the temporal inertia term in the Navier-Stokes equations, and the assumptions made in developing these treatments are being invalidated by current and projected operating speeds and loads in industrial bearings. Development of a new approach to the generalized Reynolds equation or another simplified form of the Navier-Stokes equations is an opportunity for future research.The proper dynamic model for tilting pad journal bearings is another area of research and discussion in the literature. This paper summarizes the two approaches, and new methods have been developed to directly compare the two approaches [100]. A new area of research is to experimentally identify the pad transfer functions instead of relying solely on implicit treatments based on measurements of rotor or bearing housing motion. Early work in this area has been recently reported [101]. It is anticipated that the inclusion of these measurements, in conjunction with multiple-input-multiple-output system identification techniques, will be useful in determining the correct TPJB dynamic model.



---

*Source: 908469-2011-08-22.xml*