Tilting-Pad Bearings -- Contact Flexibility of the Pivot

1. Introduction

Most readers are likely familiar with the geometry and operation of tilting-pad journal bearings (TPJB). These bearings are mechanically complex, having multiple pads with various loading arrangements, e.g., load-on-pad (LOP) or load-between-pads (LBP). They are available in many commercial variations, e.g., spherical pivots (point contact), rocker pivots (line contact), and flexure pivots. It has become well known that the flexibility of the pivots can significantly affect the dynamic properties of tilting-pad bearings. However, based on experimental measurements, the current theory for pivot stiffness determination results in pivots that are too stiff. One reason for this deficiency is likely to be neglect of surface roughness at the contact.

When two bodies with rough surfaces are pressed together, the true contact area formed between the two bodies {\displaystyle A}is much smaller than the apparent or nominal contact area. This is attributed to n. atural and machined surfaces that typically exhibit roughness features, known as asperities, that occur across a broad range of length scales. These asperities create the principal paths for the flow of heat, electricity, and force. The parameters describing the asperity distribution are statistically based on asperity size (width and height) and spacial density.

The objective of this paper is to use one of the prominent rough-surface theories combined with an extension to compute stiffness for application to the pivot stiffness problem confronting researchers in tilting-pad bearing technology. A schematic of the typical contact region for a TPJB pad is shown in Figure 1.

Although most of the ideas presented here apply to all rough-surface contact situations, this paper is primarily concerned with the line or cylindrical contact.

2. Literature Survey

2.1. Research on Pivot Flexibility

The recognition of the relative independence of each pad of a TPJB has led to a method that would allow computations for complete bearing performance based upon the assembly of results for each pad. In 1964, Lund [1] published the first paper focused on TPJB dynamic coefficients using such an assembly method. This method was known as “Lund’s pad assembly method,” and now is commonly known as the KC method (K for stiffness, and C damping). In 1977, Shapiro and Colsher [2] provided further documentation, considering eigenvalue reduction for the real part, and included example analyses. In 1983, Allaire et. al. [3] presented a similar pad assembly method, and extended the Shapiro and Colsher complex eigenvalue reductions to include both the real and imaginary parts, an important improvement related to rotor stability computations. Both of these methods involved shaft translation and individual pad rotation degrees of freedom (DOF), as well as coupling between these DOF. The KC model is often referred to as the “full bearing coefficient model,” and mentioned here as the “pad assembly method.”

In the late 1970’s and early 1980’s, researchers began to recognize the significance of the pivot supporting each pad of a TPJB, and how the pivot flexibility affected bearing dynamic characteristics. Rouch [4] implemented the pivot flexibility parameter into the pad assembly method by including a radial DOF for the pad. Reduction of the full stiffness and damping matrices for nonsynchronous frequencies indicated the significance of pivot flexibility, for rotor stability determination in particular. This work of Rouch provides an excellent description of the KC or pad-assembly method.

Following the work of Rouch, researchers involved with the study of TPJB pivot effects began using the formulas for pressure and deflection between two elastic surfaces based on the results of H. Hertz from an 1895 German paper. These formulas have generally been used to model the pad contact stiffness, and can be found in Roark [5], and are categorized in a paper by Kirk and Reedy [6] specifically for engineers involved with TPJB design and specification.

For about the next two decades, most of the follow-on works involving pivot stiffness were centered on experimentation, and there is insufficient space to discuss them all. The experiments were not specifically related to pivot stiffness, but to bearing coefficients in general. One of the earliest experiments with published results was documented in the 1999 paper by Pettinato and De Choudhury [7]. They measured TPJB dynamic properties, and found that pivot flexibility reduced effective bearing stiffness and damping, and that a key-seat (cylindrical) pivot exhibited higher stiffness and damping than did a spherical pivot.

Dmochowski [8] performed experiments on a five-pad bearing and concluded that pivot flexibility and pad inertia can have significant effects on TPJB’s dynamic properties. Pivot flexibility was noted to provide either an increase or decrease in bearing stiffness coefficients, depending on operating conditions and bearing design.

Harris and Childs [9] directly measured pivot stiffness for a single pad of a four-pad ball-in-socket TPJB. They found that unlike the Hertzian result, the load versus deflection plot was linear over the range of loads considered.

In a 2011 paper, Dimond, Younan, and Allaire presented an excellent review of the tilting-pad journal bearing literature [10]. Theoretical and modeling issues were the primary concerns, with the note that some modeling options rely on the results of experimentation. This paper provided a section on review of TPJB dynamic models, with an overview on the approximate KCM (K stiffness, C damping, M mass) experimentally-identified model. This model determines constant (frequency-independent) coefficient values, and is based on curve fitting experimental results, not first principles, thus being subject to the potential vagaries of the experimental process. The KCM model is useful in many applications, but some important parameters, typically related to single pad DOF’s, cannot be explicitly included in the analysis, and from the perspective of the current paper, pivot flexibility is one of these parameters.

Wilkes and Childs presented a paper in 2012 [11] that encompassed many aspects of TPJB’s, including pivot flexibility. The work concludes that the inclusion of pivot flexibility is a requirement for both static and dynamic property calculations.

San Andres and Tao [12] investigated aspects of pivot flexibility via parametric analysis. It was shown that there are regions of pivot stiffness where either the pivot dominates effective dynamic properties, or the fluid film dominates.

In 2019, Dang, et. al., [13] showed by comparison to experiment that Hertz contact theory overestimates pivot stiffness, particularly for a rocker-backed bearing that experienced 30% lower stiffness than determined by the Hertz formula. The measurements also showed that both TPJB effective stiffness and damping varied significantly over the frequency range.

Shi, Jin and Yuan reported numerical results on the effects of pivot design on nonlinear rotordynamics in 2019 [14]. They used a Hertzian pivot model to show that an elastic pivot results in larger journal orbits, and that pivot deformations nonlinearly increase with increasing unbalance. They showed that a spherical pivot is more compliant than a cylindrical pivot.

In 2022, Wagner and Allaire [15] studied the frequency dependency related to TPJB pivots. They showed a significant effect of pivot stiffness on the variation of TPJB dynamic properties with whirl frequency. With reference to the most important range of frequencies, up to and below synchronous, increasing pivot flexibility monotonically reduced both effective stiffness and damping curves; with a tendency to “flatten” the curves over the frequency range. It was noted that this pivot behavior is likely a reason for measured stiffness, and damping, in particular, to be claimed as constants.

The most significant conclusion drawn from the literature concerning pivot effects is that, although it has been shown that there is universal concern for the significance of pivot flexibility, there has been no related theoretical examination beyond the theory of Hertz.

2.2. Rough Surface Research

The early work of Hertz was based on elastic material behavior, allowing determination of relative deformation and contact area between two ideally smooth bodies as a result of an applied load. However, it has long been known that all surfaces must be considered “rough” on a microscopic scale.

When two rough-surfaced bodies are pressed together, only the asperities or surface peaks will be in contact, and these peaks will carry very high loads. The high loads will also cause many of the asperities to yield, thus invalidating models based on purely elastic contact. Additionally, the actual area of contact is much smaller than the Hertzian computed, or nominal area.

Although not related to line contact, but to rough spherical contact, the work of Greenwood and Tripp [16] claimed that, at low loads, actual contact pressure in the presence of roughness is one-third that of the theoretical smooth Hertzian contact pressure. Additionally, Shi and Polycarpou [17] state that contact stiffness at light loads will be approximately one-third of the Hertzian result irrespective of the roughness level, and that this approximation is often adopted by researchers. These statements give support to the need to abandon Hertzian theory for accurate determination of surface contact stiffness.

One of the first, and best known, models of asperity contact was documented by Greenwood and Williamson [18]. They extended the elastic Hertzian solution to the contact of a population of asperities having a Gaussian height distribution. The rough surface was assumed to contact a rigid, flat plane. The asperities were taken to be independent, and the bulk material beneath the asperities did not deform.

In 1986, Kagami, et. al. [19], published an approximate method for dealing with rough surfaces having cylindrical curvature, including the effects of the bulk material deformation and asperity plasticity. This paper presented a procedure and results for rough-surface compliance, and concluded that the Hertz results did not conform to either the theory or experiment.

One of the first practical rough-surface models was developed by Jackson and Green (JG) [20], and was guided by the results of finite element analysis. This method, based on separation of two flat surfaces, provides a more realistic model than that of Greenwood and Williamson by permitting asperity deformation to be elastic, elastic-plastic, or fully plastic.

The rough surface statistical line-contact model used in the current work was developed in 2012 by Beheshti and Khonsari (BK) [21], who used the JG developments, along with the theories of several other researchers, to determine contact properties for curved rough surfaces having initial contact along a line. This extension to the JG work, in addition to the asperity deformation that included plasticity, considered elastic deformation of the bulk surfaces. Predictive expressions were developed for pressure distribution, contact width, and real area of contact, and the basic theory has been followed as the basis for the current work. However, BK presented no results for gross deformation along the load line; these being results required for stiffness determination.

3. Modeling the Pivot Line Contact

Much of the presentation in this section follows the works of [20] and [21], with some variations involving the technique of solution, and the extensions to provide compliance and stiffness.

3.1. Contact of Two Flat Rough Surfaces

The modeling used here follows a statistical approach that first requires establishment of a relationship for the contact of a single asperity, and then extends this to a population of asperities with heights following a Gaussian distribution. The principal assumptions are that the asperities all have the same spherical radius at the peaks, that the peak heights vary with a Gaussian distribution, and that they deform without interacting with the neighboring asperities. Additionally, it was shown in [16] that the contact of two rough surfaces can be modeled by an equivalent single rough surface contacting an ideally smooth plane, or vice versa. The effects of sliding or adhesion have been ignored.

The interest in the present paper involves the contact of curved surfaces, but the method will first compute and tabulate the pressures caused by bringing flat surfaces together under various loads using the JG theory [20]. The pressures that are then developed due to varied surface separations at different points caused by curvature will be sampled from the tabulation.

When dealing with contacting rough surfaces, two reference planes are defined. The first is the mean of the asperity heights, and the second is the mean of the surface heights. In Figure 2, z and d denote, respectively, the asperity height, and the separation of the asperity mean and ideally flat surfaces. The surface separation h is measured from the plane defined by the mean of the original surface heights. The interference ω is defined as ω = z – d.

There are basically six parameters required to define a single rough surface for the JG computations; three for the material, and three for the surface geometry. For the material, we need,

a)

elastic modulus E,

b)

Poisson’s ratio ν, and

c)

yield stress S_y,

and for geometry,

d)

standard deviation σ of the asperity peak height,

e)

radius value β assigned to the tip of each asperity, and

f)

the density of peaks η (number of peaks per unit area).

Using a surface roughness parameter of ξ = ηβσ, McCool [22] established the relationship between the standard deviations of the two reference surfaces, i.e.,

$$\frac{{\sigma }_s}{\sigma }=\sqrt{1-\frac{3.717\times {10}^{-4}}{{\xi }^2}}\mathrm{and}\frac{y_s}{\sigma }=\frac{0.04594}{\xi }$$

(1)

where σ_s is the standard deviation of the summit (asperity) heights.

The asperity height function is made dimensionless via σ, and $\overline{\mathrm{z}}$ = z/σ. The dimensionless Gaussian probability function $\mathsf{\phi }\left(\overline{\mathrm{z}}\right),$ required to statistically define the asperity height functions is given by,

$$\mathsf{\phi }\left(\overline{\mathrm{z}}\right)=\frac{1}{\sqrt{2\mathsf{\pi }}}\left(\frac{\mathsf{\sigma }}{{\mathsf{\sigma }}_{\mathrm{s}}}\right)\exp \left[-\frac{1}{2}{\left(\frac{\mathsf{\sigma }}{{\mathsf{\sigma }}_{\mathrm{s}}}\right)}^2{\overline{\mathrm{z}}}^2\right]$$

(2)

The effective elastic modulus ${\mathrm{E}}^{\mathrm{\text{'}}}$ for the two surface materials, also used by Hertz, is defined by,

$$\frac{1}{{\mathrm{E}}^{\mathrm{\text{'}}}}=\frac{1-{\mathsf{\nu }}_1^2}{{\mathrm{E}}_1}+\frac{1-{\mathsf{\nu }}_2^2}{{\mathrm{E}}_2}$$

(3)

where E_i and ν_i are elastic moduli and Poisson’s ratio for the two materials.

The JG model considers the variation of material hardness during asperity deformation with the inclusion of material yield strength, and defines a critical interference ${\overline{\mathsf{\omega }}}_{\mathrm{c}}$that marks the transition from elastic to elastic-plastic deformation, i.e.,

$${\overline{\mathsf{\omega }}}_{\mathrm{c}}={\left(\frac{{\overline{\mathrm{S}}}_{\mathrm{y}}\mathsf{\pi }\mathrm{C}}{2}\right)}^2\frac{\mathsf{\beta }}{\mathsf{\sigma }}$$

(4)

with ${\overline{\mathrm{S}}}_{\mathrm{y}}={\mathrm{S}}_{\mathrm{y}}/\mathrm{E}\mathrm{\text{'}}$, and C = 1.295exp(0.736ν). The effective yield strength and Poisson’s ratio are associated with the softer material.

The following relationship for pressure between rough flat surfaces has been provided with the JG model, namely,

$$\mathrm{p}\left(\mathrm{h}\right)=\mathsf{\eta }\mathsf{\beta }\mathsf{\sigma }{\mathrm{E}}^{\mathrm{\text{'}}}\mathsf{\Phi }\left(\mathrm{h}\right)$$

(5)

where,

$$\mathsf{\Phi }\left(\mathrm{h}\right)=\frac{4}{3}\left\{\sqrt{\frac{\mathsf{\sigma }}{\mathsf{\beta }}}{\int }_{\overline{\mathrm{d}}}^{\overline{\mathrm{d}}+1.9{\overline{\mathsf{\omega }}}_{\mathrm{c}}}{\overline{\mathsf{\omega }}}^{3/2}\mathsf{\phi }\left(\overline{\mathrm{z}}\right)\mathrm{d}\overline{\mathrm{z}}+\frac{\mathrm{C}\mathsf{\pi }{\overline{\mathrm{S}}}_{\mathrm{y}}}{2}{\int }_{\overline{\mathrm{d}}+1.9{\overline{\mathsf{\omega }}}_{\mathrm{c}}}^{\mathrm{\infty }}\left[{\mathsf{\Gamma }}_1+{\mathsf{\Gamma }}_2{\mathsf{\Gamma }}_3\right]\mathsf{\phi }\left(\overline{\mathrm{z}}\right)\mathrm{d}\overline{\mathrm{z}}\right\}$$

with,

$${\mathsf{\Gamma }}_1\left(\overline{\mathsf{\omega }}\right)=\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{1}{4}{\left(\frac{\overline{\mathsf{\omega }}}{{\overline{\mathsf{\omega }}}_{\mathrm{c}}}\right)}^{5/12}\frac{{\overline{\mathsf{\omega }}}^{3/2}}{\sqrt{\overline{\mathsf{\omega }}}}\right]$$

$${\mathsf{\Gamma }}_2\left(\overline{\mathsf{\omega }}\right)=1-\mathrm{e}\mathrm{x}\mathrm{p}\left[-0.82{\left(\sqrt{\overline{\mathsf{\omega }}}\sqrt{\frac{\mathsf{\sigma }}{\mathsf{\beta }}}{\left(\frac{\overline{\mathsf{\omega }}}{1.9{\overline{\mathsf{\omega }}}_{\mathrm{c}}}\right)}^{\mathrm{D}/2}\right)}^{-0.7}\right]$$

$${\mathsf{\Gamma }}_3\left(\overline{\mathsf{\omega }}\right)=\frac{11.36}{\mathrm{C}}\overline{\mathsf{\omega }}\left\{1-\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{1}{25}{\left(\frac{\overline{\mathsf{\omega }}}{{\overline{\mathsf{\omega }}}_{\mathrm{c}}}\right)}^{5/9}\right]\right\}$$

and $\mathrm{D}=0.14\mathrm{e}\mathrm{x}\mathrm{p}\left(23{\overline{\mathrm{S}}}_{\mathrm{y}}\right)$, $\overline{\mathsf{\omega }}=\overline{\mathrm{z}}-\overline{\mathrm{d}}$, and $\overline{\mathrm{d}}$= d/σ = (h-y_s)/σ

3.2. Application of the Theory to Line Contacts

When considering contact of two rough curved surfaces, accuracy requires statistical consideration of the asperity layer (via the JG model) as well as bulk deformation of the material supporting the asperities.

With the contact of two surfaces of different radii, the equivalent rough surface is characterized by an effective asperity surface curvature 1/R (effective radius R) which is a combination of the curvatures of the two rough surfaces, namely,

$$\frac{1}{\mathrm{R}}=\frac{1}{{\mathrm{R}}_1}\pm \frac{1}{{\mathrm{R}}_2}$$

(6)

where the “+” sign is used if both surfaces are convex, and the “–” sign is used for a convex surface in a groove (R₂ > R₁). As applied to the current TPJB study, R₂=R_h is the “housing” radius, and R₁=R_p is the “pad-back” radius. Equation (6) is also used in development of Hertzian theory, but for purposes of understanding, and as also used with the JG theory, a cylinder with effective radius R contacts a flat, rough plane (see Figure 3).

Including both the curvature and bulk deformation arising from the distribution of contact pressures gives an equation for mean separation,

$$\mathrm{h}\left(\mathrm{x}\right)={\mathrm{h}}_{00}+\frac{{\mathrm{x}}^2}{2\mathrm{R}}-\frac{2}{\mathsf{\pi }{\mathrm{E}}^{\mathrm{\text{'}}}}\mathrm{J}\left(\mathrm{x}\right)$$

(7)

with $\mathrm{J}\left(\mathrm{x}\right)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{p}\left(\mathrm{s}\right)\ln \left(\left|\mathrm{x}-\mathrm{s}\right|\right)\mathrm{d}\mathrm{s}$, h₀₀ is a constant to be determined, and x is the coordinate perpendicular to the contact line, with origin at the apex of the contact (see Figure 3). The integral J(x) is based on a line load applying pressure to an elastic half space [23].

The force balance must also be satisfied, viz,

$$\mathrm{W}=2\mathrm{L}{\int }_0^{\mathrm{\infty }}\mathrm{p}\left(\mathrm{x}\right)\mathrm{d}\mathrm{x}$$

(8)

where W is the applied load, L is the contact length, and the pressure distribution has been assumed symmetrical about the contact line. Via the iteration process, this force balance will help to identify the constant h₀₀.

In Figure 3, δ is the relative motion of approach along the loading axis of two points, one on each body; both remote from the contact zone. This value is termed the “compliance,” and will be used during the numerical computation of stiffness.

For convenience in comparison and presentation of results, it is desirable to nondimensionalize these equations using the related Hertzian parameters. $\overline{\mathrm{W}}$ is the dimensionless load, p_H is the peak Hertzian pressure, and x_H is the Hertzian contact half width.

$$\overline{\mathrm{W}}=\frac{\mathrm{W}}{\mathrm{R}\mathrm{L}{\mathrm{E}}^{\mathrm{\text{'}}}},{\mathrm{p}}_{\mathrm{H}}=\mathrm{E}\mathrm{\text{'}}\sqrt{\frac{\overline{\mathrm{W}}}{\mathsf{\pi }}},{\mathrm{x}}_{\mathrm{H}}=\sqrt{\frac{4\overline{\mathrm{W}}{\mathrm{R}}^2}{\mathsf{\pi }}}$$

With these definitions, the dimensionless pressure will always be less than 1, and the dimensionless contact half width will always be greater than 1. Other required dimensionless parameters are,

$$\overline{\mathrm{h}}=\mathrm{h}/\mathrm{R},\overline{\mathsf{\sigma }}=\mathsf{\sigma }/\mathrm{R},\overline{\mathsf{\beta }}=\mathsf{\beta }/\mathrm{R},\overline{\mathsf{\eta }}=\mathsf{\eta }{\mathrm{R}}^2$$

Employing these parameter modifications results in a non-dimensional set of equations that defines the system to be solved for the cylindrical contact problem, as follows,

$$\overline{\mathrm{p}}\left(\overline{\mathrm{h}}\right)=\mathrm{p}\left(\mathrm{h}\right)/{\mathrm{p}}_{\mathrm{H}}=\overline{\mathsf{\eta }}\overline{\mathsf{\beta }}\overline{\mathsf{\sigma }}\overline{\mathsf{\Phi }}\left(\overline{\mathrm{h}}\right)$$

(9)

$$\overline{\mathrm{h}}\left(\mathrm{X}\right)={\overline{\mathrm{h}}}_{00}+\frac{4\overline{\mathrm{W}}}{\mathsf{\pi }}\left(\frac{{\mathrm{X}}^2}{2}-\frac{1}{\mathsf{\pi }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\overline{\mathrm{p}}\left(\mathrm{s}\right)\ln \left(\left|\mathrm{X}-\mathrm{s}\right|\right)\mathrm{d}\mathrm{s}\right)$$

(10)

$$\frac{\mathsf{\pi }}{4}={\int }_0^{\mathrm{\infty }}\overline{\mathrm{p}}\left(\mathrm{X}\right)\mathrm{d}\mathrm{X}$$

(11)

The function $\overline{\mathsf{\Phi }}\left(\overline{\mathrm{h}}\right),$ is a non-dimensional form of $\mathsf{\Phi }\left(\mathrm{h}\right)$, that can be found in [20]. The X coordinate used in these equations is nondimensionalized using the Hertzian half width, i.e., X = x/x_H.

This system of equations (9 -11) can be discretized and solved for pressure distribution using an iterative numerical scheme. The dimensionless surface separation at point i along the X axis is given by,

$${\overline{\mathrm{h}}}_{\mathrm{i}}={\overline{\mathrm{h}}}_{00}+\frac{4\overline{\mathrm{W}}}{\mathsf{\pi }}\left(\frac{{\mathrm{X}}_{\mathrm{i}}^2}{2}-\frac{1}{\mathsf{\pi }}{\sum }_{\mathrm{j}=-\mathrm{N}}^{\mathrm{N}}{\mathrm{M}}_{\mathrm{i}\mathrm{j}}{\overline{\mathrm{p}}}_{\mathrm{j}}\mathsf{\Delta }\mathrm{X}\right)$$

(12)

where ΔX are equally spaced increments, and M_ij is determined from multilevel solution techniques [21,23] shown as,

$${\mathrm{M}}_{\mathrm{i}\mathrm{j}}=\left(\mathrm{i}-\mathrm{j}+\frac{1}{2}\right)\left[\ln \left(\left|\mathrm{i}-\mathrm{j}+\frac{1}{2}\right|∆\mathrm{X}\right)-1\right]-$$

$$\left(\mathrm{i}-\mathrm{j}-\frac{1}{2}\right)\left[\ln \left(\left|\mathrm{i}-\mathrm{j}-\frac{1}{2}\right|∆\mathrm{X}\right)-1\right]$$

Load equilibrium is assured by satisfying,

$$\frac{\mathsf{\pi }}{4}={\sum }_{\mathrm{k}=1}^{\mathrm{N}}{\mathrm{D}}_{\mathrm{k}}\overline{\mathrm{p}}\left({\mathrm{X}}_{\mathrm{k}}\right)$$

(13)

where D_k = (1,4,2,4, … ,2,4,1)$∆\mathrm{X}$/3 are the weights for integration by Simpson’s Rule.

The full system of equations tends to be ill-conditioned, and cannot generally be solved in a few iterations. The principal difficulty lies in the computation of the bulk deformation represented by J(x) in equation (7), and the sum in equation (12). If the cylinder is initially set to be rigid (sum term = 0), but with rough surface, the system solves readily. In the approach used here, this is the starting point, and a first solution can be found by giving estimates for h₀₀ and the pressure distribution. The next iteration starts with the results of the first iteration, but selects a proportion of the sum term in (12), and finds a new pressure solution. This procedure is followed until the full value of the sum term has been applied and convergence for J(x) has been attained. A flowchart for this procedure is provided in Figure 4. It is straight-forward to compute compliance from the results for pressure distribution (see equation (17)).

Pivot stiffness determination for a TPJB pad first requires knowledge of the compliance or relative deformation between the two bodies. The results from a solution of the subject system provide the statistical pressure distribution and the extent of the contact half-width. These results are essential for the compliance calculation, and a comparison indicating partial validation for the pressures and the contact half width is shown in Figure 5. Both axes are normalized relative to the Hertzian result.

The surface parameters used for this comparison are displayed in Table 1 (shown later) for the rough surface No. 4, and the two loads are 500N and 5000N. The material is steel at room temperature with a 460 MPa yield stress. The effective radius is 20 mm with a 45 mm length. The current results are shown relative to the approximate curve-fitted results from the BK theory of [21], and to the Hertzian results (used for normalization). This plot shows the relative increase in pressures as load is increased, and the relative reduction of contact width with increased load, both trending toward the Hertzian distribution as a high-load limit. The influence of smoothing the surface is similar to the load increasing behavior.

3.3. Contact Compliance and Stiffness

The computation for the change of center distances for a loaded cylinder in a groove, based upon Hertzian theory, is given by Roark [4] for two contacting bodies of the same material (E=E₁=E₂). Some algebraic manipulation simplifies the Hertz deflection equation in Roark to be the following,

$${\mathsf{\delta }}_{\mathrm{H}}=\frac{2\mathrm{W}\left(1-{\mathsf{\nu }}^2\right)}{\mathsf{\pi }\mathrm{L}\mathrm{E}}\left\{\frac{2}{3}+\mathrm{l}\mathrm{n}\left[1.731\frac{\mathrm{E}\mathrm{L}({\mathrm{R}}_2-{\mathrm{R}}_1)}{\mathrm{W}}\right]\right\}$$

(14)

A Hertzian-type formula for the related contact stiffness is found, firstly, by differentiation of equation (14) with respect to load W, and by then taking the inverse, giving the result,

$${\mathrm{K}}_{\mathrm{H}}=\frac{\mathsf{\pi }\mathrm{L}\mathrm{E}}{2\left(1-{\mathsf{\nu }}^2\right)}{\left\{-\frac{1}{3}+\ln \left[1.731\frac{\mathrm{E}\mathrm{L}({\mathrm{R}}_2-{\mathrm{R}}_1)}{\mathrm{W}}\right]\right\}}^{-1}$$

(15)

Computation of rough-surface stiffness based on the JG model must begin with displacement determination δ for a given load W (see Figure 3). An accurate numerical result for stiffness determined at load W_i can then be found using a central difference formula, namely,

$${\mathrm{K}}_{\mathrm{p}\mathrm{i}}=\frac{{\mathrm{W}}_{\mathrm{i}+1}-{\mathrm{W}}_{\mathrm{i}-1}}{{\mathsf{\delta }}_{\mathrm{i}+1}-{\mathsf{\delta }}_{\mathrm{i}-1}}$$

(16)

where W_i-1 and W_i+1, are load applications each separated from W_i by a small load step ΔW, with related displacement values δ.

The computed surface separation values h(x) of eq’n. (7) cannot be used directly to determine the gross movement of the two cylindrical bodies under load. However, at the edge of the contact area, at a position x=x_E, the asperity deformation becomes nearly zero, and the total relative deformation between bodies can be expressed in the following [19],

$$\mathsf{\delta }=\frac{{{\mathrm{x}}_{\mathrm{E}}}^2}{2\mathrm{R}}+\frac{2}{\mathsf{\pi }{\mathrm{E}}^{\mathrm{\text{'}}}}{\int }_{-{\mathrm{x}}_{\mathrm{E}}}^{{\mathrm{x}}_{\mathrm{E}}}\mathrm{p}\left(\mathrm{s}\right)\ln \left(\left|{\mathrm{x}}_{\mathrm{E}}-\mathrm{s}\right|\right)\mathrm{d}\mathrm{s}$$

(17)

This result includes only the initial geometry and bulk deformation influences, and is further discussed by Beheshti and Khonsari [24]. Unlike the integral of eq’n. (7), this integral can be numerically integrated directly because the computed pressures are known. Since the integral limits ideally extend to infinity, the choice for the position x_E is practically selected for a normalized pressure $\overline{\mathrm{p}}$ that has become very small, here chosen consistently as 0.001.

3.4. Surface Parameter Identification

As opposed to nanoscale measurements used by physicists, surface measurement techniques based on mechanical or optical methods are generally sufficient for the development of the set of parameters required for definition of most manufactured surfaces. However, though accepted, and often experimentally qualified for use in current rough surface analysis, the three parameters required to approximately define the topology of a rough surface (σ, β, η) are not generally available in the common literature.

Surface roughness most commonly refers to the variations in the height of the surface relative to a reference plane. The most often used and measured surface height parameter is R_a, defined as the average absolute deviation of a profile height from its mean line. The value of R_a is an official standard in most industrialized countries, and has the relationship σ ≈ $\sqrt{\mathsf{\pi }/2}$R_a for Gaussian surfaces. In many cases, σ and R_a are interchangeable.

McCool [25] provides a method for establishing the three required contact parameters (σ, β, η) by determination of the surface “spectral moments” that are derivable from measurements. More detailed discussion related to surface finish can be found in [26,27].

Most practical contact situations do not involve a single rough surface pressing against a smooth surface; usually the contact is with both surfaces being rough. Equivalent two-surface values for σ, β, and η can be found from the following equations [21],

$${\mathsf{\sigma }}_{\mathrm{e}\mathrm{q}}=\sqrt{{\mathsf{\sigma }}_1^2+{\mathsf{\sigma }}_2^2}$$

(18)

$$\frac{1}{{\mathsf{\beta }}_{\mathrm{e}\mathrm{q}}}=\sqrt{\frac{1}{{\mathsf{\beta }}_1^2}+\frac{1}{{\mathsf{\beta }}_2^2}}$$

(19)

$$\frac{1}{{\mathsf{\eta }}_{\mathrm{e}\mathrm{q}}}=\frac{1}{{\mathsf{\eta }}_1}{\left(\frac{{\mathsf{\beta }}_{\mathrm{e}\mathrm{q}}}{{\mathsf{\beta }}_1}\right)}^2+\frac{1}{{\mathsf{\eta }}_2}{\left(\frac{{\mathsf{\beta }}_{\mathrm{e}\mathrm{q}}}{{\mathsf{\beta }}_2}\right)}^2$$

(20)

These expressions are used to define the pivot contact conditions considered in the next section.

For use in this study, four “rough” surfaces with related parameters were selected from the literature and are shown in Table 1. The “smooth” to “rough” ranking is arbitrary; used only as a relative designation. The surface No. 1 may be roughly representative of a grinding operation, while surface No. 3 could represent a milled or broached finish [25].

4. Comparison to Experiment

4.1. Individual Pad Pivot Stiffness

All of the previous discussion has focused on the pivot flexibility of individual bearing pads. The dynamic performance of a complete TPJB requires consideration of an assembly of these pads, each having different conditions of loading. Each bearing pad must generally support a different load, and since pivot flexibility has a dependency on pad load, pivot stiffness will be different for each pad. The aim of this comparison is to show that the inclusion of pivot flexibility due to surface roughness improves the computational results relative to Hertzian contact theory, and in fact, gives results that compare well with experimental results for bearing dynamic properties.

Kulhanek and Childs [28], presented an experimental investigation to measure the dynamic coefficients of a five-pad, load-between-pads TPJB, with an interest in both static and frequency dependent properties. It is clear, and in fact noted in the study, that the measured results must include pivot flexibility. Many of the results from this study are tabulated in the Kulhanek thesis [29]; conveniently accessible for the current work. It was stated by Wilkes and Childs [11] that the method given by Kulhanek provided 95% confidence bounds for the measured bearing impedances.

The bearing under consideration has a diameter of 101.6 mm, an aspect ratio (L/D) of 0.6, central pivots, a machined (cold) pad clearance of 112 μm, a relatively small preload ratio of 0.27, and with the lubricant being ISO VG32 oil supplied at 44°C. Each pad has a mass of 0.44 kg and a mass moment of inertia about the pivot of 273.9 kg-mm². Detailed pivot characteristics were not included in either of the references. However, a photo of the bearing pad and pivot region is shown in the Kulhanek thesis [29], permitting estimation of some geometric details of the pivot.

The pivot used for this bearing is of cylindrical form; assumed to have the same length as the bearing pad. The “cylinder” sits in a groove that appears to be only slightly larger in radius. The pivot cylinder radius is scaled (very roughly) as 20 mm, and the groove radius is taken as 23 mm. The assumed material is the same steel for both contacting surfaces, with a Poisson’s ratio of 0.29, a yield stress of 360 MPa, and an elastic modulus of 190 GPa (slightly reduced from room temperature value for operating temperature).

Using this estimated pivot data in eq’n. (15), Hertzian pivot stiffness data as a function of pad load was computed with the results shown in Figure 6. Also shown, using eq’n. (16), are the numerical results of the BK extended formulation for the four surfaces of Table 1. This plot clearly displays the inadequacy of the Hertz formulation, especially for lower loads that show one-third the stiffness relative to the Hertz result.

The particular bearing chosen here was also considered in the work of Wagner and Allaire [15]. In that study, a rotational speed of 7000 rpm with a specific load (W/LD) of 1.723 MPa was used to estimate the magnitude of the pivot stiffness by comparison to measured data. The load applied to each of the loaded pads was 6500 N, and it was estimated that the stiffness of the film for each loaded pad was about 585 MN/m. It was shown that the best match to the measured data was obtained with a pivot stiffness of approximately 2× film stiffness, or 1170 MN/m. According to Figure 6, this is about one-half the Hertzian value, and is in the center of the stiffness range for the rough surfaces on the plot.

4.2. Complete Bearing Properties

Inclusion of pivot stiffness in TPJB modeling requires explicit consideration of the parameters governing the dynamics of each individual pad. The pad assembly or KC method would be appropriate for such analyses. The computer code used here is a KC type code based on a thermohydrodynamic (THD) method developed from theoretical formulations provided by Suganami and Szeri [30] and the text of Szeri [31]; further extended to include a pad radial degree-of-freedom for consideration of pivot flexibility [4]. Though undoubtedly of significance at higher speeds and loads, mechanical and thermal pad bending and growth effects have not been considered.

The mentioned KC-THD code was used to compute synchronous stiffness and damping coefficients for the subject bearing at both 50% and 60% pivot offsets,1 and with computations at a rotational speed of 7000 rpm over a range of representative bearing loads. Although the results for other speeds can be presented, the intent here is not to provide data for design, but to illustrate the impact of a significant effect that has not been previously considered. Furthermore, computations at higher speeds would not be as accurate since changes in clearances and effective preloads due to more pronounced thermal effects could not be considered.

No information was presented in the Kulhanek thesis [29] regarding the roughness state of the pivot contacting surfaces. After surveying the full-bearing dynamic coefficient results, using the four surfaces identified in Table 1, it was found that Surface No. 1 provided results that agreed most closely with measured data. This surface description, approximately conforming to a surface of medium grind, has been used for all pads of the subject bearings. Table 2 displays the individual pad loads and related pivot stiffness values (50% pivot offset) for two selected specific loads of the 7000 rpm computations. Pad numbering is CCW following rotation, and pads 1 and 2 straddle the load line. The bearing coordinate system is right handed, and uses positive X in the direction of the bearing load, as was generally favored by Lund [1].

The computed stiffness and damping coefficients for this bearing are shown in Figure 7, Figure 8, Figure 9 and Figure 10 for both 50 and 60 percent pivot offset positions. Only the direct coefficient results are shown; the cross-coupled terms are generally small and ignored. The solid lines on these plots apply to the results for the rough surface pivot contact, while the dashed curves refer to the results for Hertzian contact. The individual plotted points relate to the measured and tabulated impedance results from [29]. The pad numbering scheme is shown on the plots.

The stiffness results shown on the plots are the real values of the frequency-dependent complex impedance as evaluated at rotational frequency Ω. These stiffness coefficients are defined differently than those obtained from curve fits to the impedance functions as used by proponents of KCM modeling. The KCM fits are based on an approximate Taylor series expansion, and Kulhanek [29] has presented “stiffness” coefficients k_t as the constant part of the real impedance that has been decomposed into a constant term, and an Ω² term with a constant “mass” coefficient m_t. The measured points on the stiffness plots of Figure 7 and Figure 9 have been reconstructed as k = k_t - m_tΩ² in order to match the original measured impedance real part. The damping coefficients have been determined by dividing the imaginary part of the impedance by Ω.

It is acknowledged that the computed results for the complete TPJB, though reasonably comprehensive, do not consider all of the effects that may have an impact on the results (e.g., clearance change effects with temperature, convective and temporal inertia effects). However, with reference to the effects of the pivot, it is clear from the examples shown that the use of Hertzian pivot stiffness overestimates the complete bearing properties. For these examples, the overestimation is up to 20% for the effective stiffness coefficients, and 30% for the effective damping coefficients.

5. Discussion

This study has concentrated on cylindrical or rocker-back pivots. However, many bearings used in industry have spherical (also termed “circular”) or elliptical contact areas. For a given load, these contacts function with much higher pressure densities than do cylindrical pivots, thus causing full plastic deformation for many of the asperities. Additionally, a more significant portion of the total contact deformation may be associated with elastic bulk deformation. This situation is thus very similar to elastic Hertzian contact, and this has been noted in the literature [13], where a plot shows relatively small difference between measured circular pivot stiffness and the related Hertzian stiffness.

A question comes to mind concerning the paucity of research related to pivot stiffness. Firstly, this issue is of concern only to those interested in TPJB’s, and principally only for dynamic influences. Secondly, many researchers and designers may have considered the “problem solved,” especially after the work of Kirk and Gordon [6] that outlines Hertzian theory for application to pivot stiffness. Since many TPJB designs use the spherical pivot design, and that spherical contact stiffnesses are much closer to Hertzian stiffness results than are line stiffnesses, the idea may have been prominent that there was no need for further research. Finally, it has only been within the relatively recent years that experimentalists began to realize that measured TPJB dynamic properties poorly matched computations, and that the deficiency was likely to be pivot stiffness related to line contacts [11,12,13,14,15].

A paper of significance to TPJB designers; concerned with modeling circular and elliptical contacts, has been published by Beheshti and Khonsari [24]. This is a presentation of a modeling concept for such contacts, along with curve-fitted formulas based on the numerical simulations. These formulas can be readily used for the prediction of the maximum contact pressure, contact dimensions, contact compliance, real area of contact, and pressure distribution. Since the formulas give continuous results, it is a simple matter to numerically determine compliances for given loads, thus permitting estimation of contact stiffness. The same authors presented similar predictive formulas for cylindrical contacts [21], but did not include formulas for compliance; equations that could be quite useful to bearing designers, and a direction for future work.

In lieu of practical contact stiffness formulas for cylindrical contacts, both this study, and that of Dang et. al. [13], imply that a reasonable approximation for TPJB line-pivot stiffness; one better than Hertzian stiffness, would be to simply use one-half of the computed Hertzian value. Such an approximation, though not strongly justified, would be quite useful, particularly to the rotordynamic community.

The results of this study can also be applied to tilting-pad thrust bearings, particularly (but not exclusively) of the equalizing type, and especially if axial dynamics are of concern. Such bearings contain many moving parts, including leveling plates or linkages that are meant to distribute the load equally among pads. The pad/linkage system is composed of multiple contact points, each having flexibility affected by surface contact conditions. Just as determined for the TPJB, these contacts may significantly decrease the stiffness of the thrust bearing system for axial dynamics.

6. Conclusions

The focus of this study has been to investigate the pad pivot support for tilting-pad bearings in response to the recognition by multiple researchers that the pivot may be significantly more flexible than determined theoretically. Line or cylindrical pivots are the principal concern. Based on this effort, and the material drawn from related references, the following conclusions have been established.

Conclusions focused on TPJB performance are:

• The Hertz formula for stiffness, traditionally used for pivot stiffness determination, has been shown to greatly overestimate pivot stiffness, particularly for line contacts.
• Computed dynamic coefficient results for a complete TPJB, using rough surface pivot stiffness values, show significant improvement relative to results based on Hertzian pivot stiffnesses.
• Inclusion of pivot stiffness in TPJB calculations requires a method that considers the individual pad parameters, such as the “pad assembly method,” or the “KC method.”

Conclusions related to contact mechanics are:

• Either an increase in surface smoothness, or an increase in pad loading, will cause an increase in pivot contact stiffness.
• All machined surfaces over the nominal contact area are “rough”, on a microscopic scale, such that the actual support area is less than the nominal. This actual contact area is primarily composed of the asperity tips.
• The deformations causing the pivot flexibility are related to both the asperities and to the bulk deformations of the bodies supporting the asperities.
• The maximum statistical pressure for smooth surfaces is greater than that for rough surfaces, and it is closer to the Hertzian result. This trend also holds as the load is increased. Additionally, the contact area, relative to Hertzian contact, becomes smaller for smoother surfaces and higher loads.