E ﬀ ects of Tilting Pad Journal Bearing Design Parameters on the Pad-Pivot Friction and Nonlinear Rotordynamic Bifurcations

: This study numerically analyzes and investigates the e ﬀ ects of the bearing design parameters of a tilting pad journal bearing (TPJB) on the pad-pivot friction-induced nonlinear rotordynamic phenomena and bifurcations. The bearing parameters were set to the pad preload, pivot o ﬀ set, spherical pivot radius, and bearing length to diameter (L / D) ratio. The Stribeck curve model (SCM) model was applied at the contact surface between the pad and the pivot, which varied to the boundary-mixed-ﬂuid friction state depending on the friction condition. The rotor-bearing model was set up with a symmetrical ﬁve-pad TPJB system supporting a Je ﬀ cott type rigid rotor. The ﬂuid repelling force generated in the oil ﬁlm between each pad and the shaft was calculated using a ﬁnite element method. The simulation recurrently conducted the transient numerical integration to obtain the Poincar é maps and phase states of the journal and pad with various bearing design variables, then the nonlinear properties of each condition were analyzed by expressing the bifurcation diagrams. the original ﬁndings of this study are: (1) The pad preload and pivot o ﬀ set signiﬁcantly inﬂuenced the emergence of Hopf bifurcations and the associated limit cycles. In contrast, (2) the pivot radius and L / D ratio contributed relatively less to the friction-induced instability. Resultantly, (3) all the e ﬀ ects diminished when the rotor operated under the larger mass eccentricity of the disc.


Introduction
A tilting pad journal bearing (TPJB) is known as a structure in which the lower surface of the pad executes a tilt motion on the upper surface of the pivot, with respect to the rotation of the shaft. This tilt motion suppresses the cross-coupled stiffness of the bearing and improves stability [1]. As a result, TPJBs have been widely selected in many modern turbomachinery.
The typical structure of a pivot can be categorized into the rocker back (cylindrical) type and the ball-socket (spherical) type. Various researchers have analyzed the differences in the performance of the two structures. Wygant et al. [2,3] conducted a series of experiments to measure the eccentricity ratios, dynamic stiffness, and damping force coefficients of the two types of pivots in the TPJB system. They found variations in the attitude angles and cross-coupled stiffness between the two types of TPJB and concluded that the higher pad-pivot friction in the ball and socket pivot was the primary contributor to the differences. Pettinato and De Choudhury [4,5] measured the differences in the pad temperature profiles and the power loss in a five-shoe TPJB of the two types of pivots. Mehdi et al. [6] developed a numerical model to investigate the static and dynamic characteristics of TPJBs by considering the (2) In this study, the friction coefficient, μf, was determined using the Stribeck curve model; therefore, the surface roughness condition and sliding speed of the pad-pivot were considered. Figure 1. Schematics of the pad-pivot friction mechanism in a spherical pivot-type tilting pad journal bearing .

Stribeck Curve Model (SCM)
The SCM distinguishes the friction states as the boundary, mixed and hydrodynamic regimes, and thus suggests a more realistic friction coefficient based on the relative motion between the surfaces. There are various factors to consider in the friction coefficient, such as the surface roughness, static load, and lubricant state. Nonetheless, it was assumed that the friction coefficient was determined only by the sliding speed of the pad under the specific surface and lubricant conditions, as described in Tables 1 and 2. This is a common scheme to interpret the SCM. Lu et al. examined the relation between the mathematical derivation and experimental results of this approach in [14].
The resulting friction coefficient, as shown in Figure 2, was used in every time step in the numerical integration process for nonlinear response analysis. This is according to the bearing design parameters discussed in the next section.  In this study, the friction coefficient, µ f , was determined using the Stribeck curve model; therefore, the surface roughness condition and sliding speed of the pad-pivot were considered.

Stribeck Curve Model (SCM)
The SCM distinguishes the friction states as the boundary, mixed and hydrodynamic regimes, and thus suggests a more realistic friction coefficient based on the relative motion between the surfaces. There are various factors to consider in the friction coefficient, such as the surface roughness, static load, and lubricant state. Nonetheless, it was assumed that the friction coefficient was determined only by the sliding speed of the pad under the specific surface and lubricant conditions, as described in Tables 1 and 2. This is a common scheme to interpret the SCM. Lu et al. examined the relation between the mathematical derivation and experimental results of this approach in [14].
The resulting friction coefficient, as shown in Figure 2, was used in every time step in the numerical integration process for nonlinear response analysis. This is according to the bearing design parameters discussed in the next section.     Figure 3 represents the schematics of the five-pad TPJB model employed in this study. The hydrodynamic pressure of the lubricant film on each pad, p, is able to be estimated using the Reynolds equation for iso-viscous, incompressible, and laminar flow conditions as follows: Figure 3. Five-pad tilting pad journal bearing (TPJB) (load on pad) diagram and the respective coordinates. 3 3

Rotor-Bearing Model
3.1. Five-Pad Tilting Pad Journal Bearing (TPJB) Figure 3 represents the schematics of the five-pad TPJB model employed in this study. The hydrodynamic pressure of the lubricant film on each pad, p, is able to be estimated using the Reynolds equation for iso-viscous, incompressible, and laminar flow conditions as follows: where z and θ are the axial and angular positions of the bearing, respectively; µ v is the dynamic viscosity of the lubricant oil; ω is the angular velocity of the journal; and h is the film thickness on the pad at the {z, θ} position: where e is the vector that designates the current position in the reference frame {x, y, z}; θ p is the pivot position in the θ axis; δ is the tilt angle of the pad; and C b and C p are the clearances of the bearing and pad, respectively. Superscript (j) denotes the jth pad, and subscripts x, y, z, pitch and roll represent the relevant axes and motions. The solution of the Reynolds equation was calculated by applying a finite element model, as shown in Figure 4. It consisted of a three-node simplex, triangular-type mesh created on the overall fluid film layer on the pad. The fluid repelling force between the pad and shaft was obtained by integrating the pressure throughout the mesh for each pad as follows: where F (j) is the hydrodynamic force on the jth pad, L is the bearing length, θ B (j) is the beginning angle of the jth pad, and θ E (j) is the end angle of the jth pad. The boundary conditions was incorporated with an assumption of ambient pressure at the bearing ends (p amb = p| z = ± L/2 ) and supply pressure at the leading and trailing edges (p sup = p| θ = θB (j) , θE (j) ). In addition, Reynolds cavitation boundary conditions were applied in the calculation. The integrated pad pressure functions as the normal force, which induces the friction moment in Equation (2).
Hydrodynamic (|dδ/dt| > 70 rad/s) : μf > 0.03  Figure 3 represents the schematics of the five-pad TPJB model employed in this study. The hydrodynamic pressure of the lubricant film on each pad, p, is able to be estimated using the Reynolds equation for iso-viscous, incompressible, and laminar flow conditions as follows:   where e is the vector that designates the current position in the reference frame {x, y, z}; θp is the pivot position in the θ axis; δ is the tilt angle of the pad; and Cb and Cp are the clearances of the bearing and pad, respectively. Superscript (j) denotes the jth pad, and subscripts x, y, z, pitch and roll represent the relevant axes and motions. The solution of the Reynolds equation was calculated by applying a finite element model, as shown in Figure 4. It consisted of a three-node simplex, triangular-type mesh created on the overall fluid film layer on the pad. The fluid repelling force between the pad and shaft was obtained by integrating the pressure throughout the mesh for each pad as follows: (5) where F (j) is the hydrodynamic force on the jth pad, L is the bearing length, θB (j) is the beginning angle of the jth pad, and θE (j) is the end angle of the jth pad. The boundary conditions was incorporated with an assumption of ambient pressure at the bearing ends (pamb = p|z = ± L/2) and supply pressure at the leading and trailing edges (psup = p|θ = θB (j) , θE (j) ). In addition, Reynolds cavitation boundary conditions were applied in the calculation. The integrated pad pressure functions as the normal force, which induces the friction moment in Equation (2).

Five-Pad Tilting Pad Journal Bearing (TPJB)
The overall hydrodynamic force from entire pads to the journal can be calculated as: The moment on the pad about its pivot owing to the pressure distribution can be expressed as: where r is the vector from the pivot contact point on pad j to the location of the differential force on pad j.  The overall hydrodynamic force from entire pads to the journal can be calculated as: The moment on the pad about its pivot owing to the pressure distribution can be expressed as: where → r is the vector from the pivot contact point on pad j to the location of the differential force on pad j.

Rigid Jeffcott Rotor-TPJB Model
As can be seen in Figure 5, a "Jeffcott"-type symmetric rigid rotor supported on the five-pad TPJB was employed as a mechanical model to investigate the rotordynamic bifurcation analysis owing to the pad-pivot friction. The governing equations of motion for the rotor and TPJB pads in the system can be written as: m J ..
where m J and I p are the rotor mass and pad moment of inertia, respectively; and W s and W d are the static and dynamic loads on the rotor, respectively.
.. δ roll are second derivatives of pad pitch and roll motion respectively. Generally, the static load is the weight of the rotor or a side load, and the dynamic force can be unbalance force due to the mass eccentricity on the disc. The separate moments of inertia equations were written for each pad (j). The numerical values of the disc, bearing, pad parameters, and pressure boundary condition have been provided in Table 3. In particular, the amount of unbalance on disc, bearing diameter to length ratio, pad preload, pad offset, and pivot radius are the major design parameters in the TPJB performance [15,16]. It is noteworthy that the present rotor model is assumed to have a rigid shaft and symmetrical geometry. This implies that the two journals have the same attitude and no roll motion in the y-z plane. Hence, the second equation in Equation (7) and fourth equation in Equation (8)

Rigid Jeffcott Rotor-TPJB Model
As can be seen in Figure 5, a "Jeffcott"-type symmetric rigid rotor supported on the five-pad TPJB was employed as a mechanical model to investigate the rotordynamic bifurcation analysis owing to the pad-pivot friction. The governing equations of motion for the rotor and TPJB pads in the system can be written as: where mJ and Ip are the rotor mass and pad moment of inertia, respectively; and Ws and Wd are the static and dynamic loads on the rotor, respectively. ̈ and y are second derivatives of rotor translational motion in x and y direction respectively. ̈ℎ ̈ are second derivatives of pad pitch and roll motion respectively. Generally, the static load is the weight of the rotor or a side load, and the dynamic force can be unbalance force due to the mass eccentricity on the disc. The separate moments of inertia equations were written for each pad (j). The numerical values of the disc, bearing, pad parameters, and pressure boundary condition have been provided in Table 3. In particular, the amount of unbalance on disc, bearing diameter to length ratio, pad preload, pad offset, and pivot radius are the major design parameters in the TPJB performance [15,16]. It is noteworthy that the present rotor model is assumed to have a rigid shaft and symmetrical geometry. This implies that the two journals have the same attitude and no roll motion in the y-z plane. Hence, the second equation in Equation (7) and fourth equation in Equation (8)

Numerical Results
The lubricant film pressure, force, and moment of inertia were numerically calculated using the finite element TPJB model, and the MATLAB ® (R2020a, MathWorks, Natick, MA, USA, 2020) routine "ode15s" was utilized for numerical integration (NI) of Equation (8) with a relative tolerance of 10 −9 to secure convergence and accuracy. The initial conditions for the transient NI were basically set as a free falling from the center (i.e., x(0) = y(0) = 0, dx(0)/dt = dy(0)/dt = 0) and δ pitch = 0, dδ pitch /dt = 0 for all pads. For consecutive run-up operation, the final states of the current operation speed were defined as the initial value of the next operation speed.

Anlaysis
The use of the transient NI can yield a Poincaré section, which is a hypersurface in the state space that is transverse to the flow of the governing equation of a system [17]. As shown in Figure 6, the consecutive collections of Poincaré dots clearly indicate the stability of the response in terms of frequency components and bifurcation occurrences if it combines with some control parameters; in this case, the bearing parameters. To obtain the Poincaré dots on the sections, NI was performed for 400 revolution periods for each selected TPJB design parameter, and the steady states were assumed to ensure for the last 100 revolutions. The projected bifurcation diagrams in Figures 7-10 plot the non-dimensional vertical journal motion (i.e., y/C b versus the selected design parameter), with the revolution speed ranging from 0 rpm to 20 krpm. The simulations were recurrently conducted with various disc imbalance eccentricities from 0.0 to 0.2 C b . The windows in the figures are placed to introduce reference plane to compare instability onsets.

Pivot Radius
The pivot radius might be intuitively designed for the pad-pivot friction moment because it is directly related to the amount of the friction moment, as shown in Equation (2). The selected pivot radius ranged from 0.01 m to 0.02 m. As can be seen in Figure 7, a larger pivot radius indicates a higher dynamic eccentricity of the friction-induced instability. On the contrary, the large size of the pivot delayed the occurrence of the Hopf bifurcation marginally; and became more pronounced in the form of a subcritical type of Hopf bifurcation. As the disc had unbalance eccentricity, the response appears to have a Neimark-Sacker (N-S) bifurcation in the high speed range. This results in the quasi-periodic motion suddenly changing to a 1/3 sub-synchronous response, and vice-versa.

Numerical Results
The lubricant film pressure, force, and moment of inertia were numerically calculated using the finite element TPJB model, and the MATLAB ® (R2020a, MathWorks, Natick, MA, USA, 2020) routine "ode15s" was utilized for numerical integration (NI) of Equation (8) with a relative tolerance of 10 −9 to secure convergence and accuracy. The initial conditions for the transient NI were basically set as a free falling from the center (i.e., x(0) = y(0) = 0, dx(0)/dt = dy(0)/dt = 0) and δpitch = 0, dδpitch/dt = 0 for all pads. For consecutive run-up operation, the final states of the current operation speed were defined as the initial value of the next operation speed.

Anlaysis
The use of the transient NI can yield a Poincaré section, which is a hypersurface in the state space that is transverse to the flow of the governing equation of a system [17]. As shown in Figure 6, the consecutive collections of Poincaré dots clearly indicate the stability of the response in terms of frequency components and bifurcation occurrences if it combines with some control parameters; in this case, the bearing parameters. To obtain the Poincaré dots on the sections, NI was performed for 400 revolution periods for each selected TPJB design parameter, and the steady states were assumed to ensure for the last 100 revolutions. The projected bifurcation diagrams in Figures 7-10 plot the non-dimensional vertical journal motion (i.e., y/Cb versus the selected design parameter), with the revolution speed ranging from 0 rpm to 20 krpm. The simulations were recurrently conducted with various disc imbalance eccentricities from 0.0 to 0.2 Cb. The windows in the figures are placed to introduce reference plane to compare instability onsets.

Pivot Radius
The pivot radius might be intuitively designed for the pad-pivot friction moment because it is directly related to the amount of the friction moment, as shown in Equation (2). The selected pivot radius ranged from 0.01 m to 0.02 m. As can be seen in Figure 7, a larger pivot radius indicates a higher dynamic eccentricity of the friction-induced instability. On the contrary, the large size of the pivot delayed the occurrence of the Hopf bifurcation marginally; and became more pronounced in the form of a subcritical type of Hopf bifurcation. As the disc had unbalance eccentricity, the response

Pad Preload
The tilt pad preload is the most tunable design parameter for the TPJB. The preload can be defined as This represents the ratio between the radii of curvature of the bearing and pad. The amount of preload is known to influence the damping aspect of the bearing [16]. The values of the preloads were set from 0.5 to 0.6, which is within the general application range.
As shown in Figure 8, the proper application of a pad preload can induce a significant change in terms of bifurcation emergence, such that the Hopf bifurcation occurrence near 14 krpm is notably set aside, while maintaining the amplitude of responses. In addition, it can be seen that this tendency remains even in a disc unbalance state. Similar to the adjustment of the pivot radius, it is possible to confirm the appearance of the N-S bifurcation, where the quasi-periodic is altered mutually in a sub-synchronous response, in the high-speed operation state.

Pivot Offset
Another practical parameter available to the TPJB geometry is the pivot offset, which is known to affect the load capacity of the bearing. The centrally pivoted pad, χp/χ = 0.5, has a pivot on the center back. The typical offset range in this study was set as 0.5 to 0.6 (50% to 60% offset). As seen in Figure 9, the larger pivot offset tends to display a lower response amplitude, and delays the occurrence of Hopf bifurcation near 14 krpm. In cases with lower unbalance, the offset significantly pushes back the Hopf bifurcation and stabilize the response states, whereas, the Neimark-Sacker bifurcation emerges during the high speed operation near 17-19 krpm. If the rotor system has a higher imbalance on the disc, the period doubling bifurcation and some peculiar occurrences, such as 1/2 sub-synchronous response, were observed before the Hopf bifurcation event; this trend is more prominent when the TPJB has more pivot offset.

Pivot Offset
Another practical parameter available to the TPJB geometry is the pivot offset, which is known to affect the load capacity of the bearing. The centrally pivoted pad, χp/χ = 0.5, has a pivot on the center back. The typical offset range in this study was set as 0.5 to 0.6 (50% to 60% offset). As seen in Figure 9, the larger pivot offset tends to display a lower response amplitude, and delays the occurrence of Hopf bifurcation near 14 krpm. In cases with lower unbalance, the offset significantly pushes back the Hopf bifurcation and stabilize the response states, whereas, the Neimark-Sacker bifurcation emerges during the high speed operation near 17-19 krpm. If the rotor system has a higher imbalance on the disc, the period doubling bifurcation and some peculiar occurrences, such as 1/2 sub-synchronous response, were observed before the Hopf bifurcation event; this trend is more prominent when the TPJB has more pivot offset.  It is known that a higher bearing L/D ratio enhances the effective damping by further expanding the fluid layer. Here, the L/D ratio was selected from the range of 0.5 to 1.0, which is the practical range for turbomachineries. As shown in Figure 10, unlike the other parameters in this study, the L/D ratio does not show a consistent and notable effect on the friction-induced instability except for the fully balanced disc condition. This indicates that the amplitude maintains a similar magnitude, and the onsets of Hopf bifurcations are slightly moved to the higher operating speed. On analyzing the overall response bifurcation diagram, it is tedious to determine a clear tendency in the unbalanced cases; therefore, the L/D ratio should be selected carefully in consideration of the bearing load and operating speed. (a)

Pivot Offset
Another practical parameter available to the TPJB geometry is the pivot offset, which is known to affect the load capacity of the bearing. The centrally pivoted pad, χ p /χ = 0.5, has a pivot on the center back. The typical offset range in this study was set as 0.5 to 0.6 (50% to 60% offset). As seen in Figure 9, the larger pivot offset tends to display a lower response amplitude, and delays the occurrence of Hopf bifurcation near 14 krpm. In cases with lower unbalance, the offset significantly pushes back the Hopf bifurcation and stabilize the response states, whereas, the Neimark-Sacker bifurcation emerges during the high speed operation near 17-19 krpm. If the rotor system has a higher imbalance on the disc, the period doubling bifurcation and some peculiar occurrences, such as 1/2 sub-synchronous response, were observed before the Hopf bifurcation event; this trend is more prominent when the TPJB has more pivot offset.

Bearing Length to Diameter (L/D) Ratio
It is known that a higher bearing L/D ratio enhances the effective damping by further expanding the fluid layer. Here, the L/D ratio was selected from the range of 0.5 to 1.0, which is the practical range for turbomachineries. As shown in Figure 10, unlike the other parameters in this study, the L/D ratio does not show a consistent and notable effect on the friction-induced instability except for the fully balanced disc condition. This indicates that the amplitude maintains a similar magnitude, and the onsets of Hopf bifurcations are slightly moved to the higher operating speed. On analyzing the overall response bifurcation diagram, it is tedious to determine a clear tendency in the unbalanced cases; therefore, the L/D ratio should be selected carefully in consideration of the bearing load and operating speed. ratio does not show a consistent and notable effect on the friction-induced instability except for the fully balanced disc condition. This indicates that the amplitude maintains a similar magnitude, and the onsets of Hopf bifurcations are slightly moved to the higher operating speed. On analyzing the overall response bifurcation diagram, it is tedious to determine a clear tendency in the unbalanced cases; therefore, the L/D ratio should be selected carefully in consideration of the bearing load and operating speed.   Figures 11-14 provide representative examples of the response of the rotor-TPJB that might demonstrate the characteristics of the design parameter to the effects of pivot friction. The results are shown in the orbit, frequency, phase portrait, and Poincare attractor formats, which correspond to the selected conditions in Figures 7-10. Figure 11a shows that the response is 1/3 sub-synchronous motion based on the attractor and frequency spectrum. It can be seen that when the radius of the pivot was increased (R pvt = 0.01 0.02), under the corresponding operating conditions, the response changed to quasi-periodic type with a higher amplitude, as can be seen in Figure 11b. Figure 12 effectively shows the effect of pad preload. An appropriately increased preload (m p = 1/2 3/4) can alter the limit cycle, as shown in Figure 12a, to the equilibrium state, as in Figure 12b; it contributes to the improvement in rotordynamic stability. Figure 13 shows that the quasi-periodic motion, as seen in Figure 13a, changes to a quenched 1×-synchronous response with a lower amplitude, as shown in Figure 13b, if the pivot offset is applied properly, e.g., χ p /χ = 0.5 0.6, even in a large unbalanced eccentricity on the disc. Lastly, Figure 14 shows an example where the extended fluid film layer (L/D = 0.5 1.0) stabilizes the friction-induced quasi-periodic motions not substantially but partially to 1/2×-subsynchronous response.

Orbits and Pad Motions
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 18 seen in Figure 13a, changes to a quenched 1×-synchronous response with a lower amplitude, as shown in Figure 13b, if the pivot offset is applied properly, e.g., χp/χ = 0.5 0.6, even in a large unbalanced eccentricity on the disc. Lastly, Figure 14 shows an example where the extended fluid film layer (L/D = 0.5 1.0) stabilizes the friction-induced quasi-periodic motions not substantially but partially to 1/2×-subsynchronous response.

Conclusions
In this study, TPJB design parameters such as the pivot radius, pad preload, pivot offset, and L/D ratio were selected and their effects on the pad-pivot friction-induced rotordynamic instability and the associated bifurcations were studied. Numerical techniques such as Poincaré sections were employed for rotor spin speed ranges (bifurcation diagram) accumulated with regard to the TPJB design parameters. The primary observations from this numerical study are summarized as follows: Pivot radius:  An increase in the pivot radius induces higher vibration amplitude;  The Hopf bifurcation event was marginally delayed;  The higher disc mass eccentricity condition undermined the effect of the pivot radius. Pad preload:  An increase in the pad preload significantly delayed the outbreaks of Hopf bifurcation points;  The amplitude of the response remained relatively constant;  In the larger disc unbalance condition, the preload stabilized the instability.

Conclusions
In this study, TPJB design parameters such as the pivot radius, pad preload, pivot offset, and L/D ratio were selected and their effects on the pad-pivot friction-induced rotordynamic instability and the associated bifurcations were studied. Numerical techniques such as Poincaré sections were employed for rotor spin speed ranges (bifurcation diagram) accumulated with regard to the TPJB design parameters. The primary observations from this numerical study are summarized as follows: Pivot radius: • An increase in the pivot radius induces higher vibration amplitude; • The Hopf bifurcation event was marginally delayed; • The higher disc mass eccentricity condition undermined the effect of the pivot radius.
Pad preload: • An increase in the pad preload significantly delayed the outbreaks of Hopf bifurcation points; • The amplitude of the response remained relatively constant; • In the larger disc unbalance condition, the preload stabilized the instability.
Pivot offset: • An increase of the pivot offset delayed the outbreaks of Hopf bifurcation points; • The amplitude of the response decreased; • The larger disc mass unbalance undermines the effect of pivot offset.
L/D ratio: • A higher L/D ratio tended to stabilize the response; however, it did not display any conspicuous effect. Nevertheless, the fully balanced condition was clearly observed; • An increase in the disc mass eccentricity undermined the effect of the L/D ratio; • Nonetheless, a higher L/D ratio led to an enhanced damping effect, which stabilized the quasi-periodic to the 1/2 sub-synchronous responses.
The pad preload and pivot offset play significant roles in decreasing the rotordynamic instability induced by the pad-pivot friction effects. Accurate selection of pad design parameters is essential, especially for high-speed, high-unbalance operation environments in rotor-TPJB systems. These conclusions serve a general purpose. However, they do not essentially hold true for all turbomachineries.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
The majority of symbols and notations used throughout the paper are defined below for quick reference. Others are clarified with their appearance in case of need.