Dynamic Coefficients of Tilting Pad Bearing by Perturbing the Turbulence Model

Abstract

Tilting pad bearings are appropriate for the trend of high efficiency and reliability design of rotating machinery due to their high stability. The laminar and turbulent flow states exist in the lubricating oil film of high-speed and heavy-load tilting pad bearings simultaneously. By perturbing the multiple flow state lubrication model with a partial derivative method, together with the pad-pivot structural perturbations, the frequency-dependent stiffness and damping coefficients of tilting pad bearings, embracing the effect of dynamical variations of both turbulence and pressure-viscous, were numerically solved in this research. The importance of each perturbed variable was studied, and the results indicate that the perturbed film thickness included in turbulence coefficients perturbations is significant enough to be taken into account otherwise the equivalent stiffness coefficients will be obviously overestimated. Unlike the perturbed film thickness, the consideration of the perturbed viscosity is optional, because it makes the stiffness and damping coefficients larger at both laminar and turbulent flow states. For a simplified simulation and conservative prediction results, the perturbed viscosity can be neglected.

1. Introduction

Tilting pad bearings have been widely applied in many kinds of rotating machines, such as gas turbines [1], compressors [2], and so on, due to their ability to improve stability [3]. Tilting pad bearings are also suitable for oil-free operations [4,5,6], so they have the potential for new fields such as electrical submersible pumps [7,8].

The dynamic stiffness and damping coefficients of bearings prominently affect the dynamic behavior of a rotor system [9], including critical speed and stability [10,11]. The pad assembly method was firstly established by Lund [12] for calculating the synchronous dynamic coefficients of tilting pad bearings. Lund’s method requires that the excitation frequency of the pad and shaft is assumed synchronous with the shaft rotation frequency. However, the use of synchronous coefficients tends to overestimate stability, especially in the case of a high Sommerfeld number [13] so the sub-synchronous dynamic coefficients are quite important. Suh studied the effects of thermal boundary condition on the stiffness and damping coefficients for tilting pad bearings [14]. Experimental studies of tilting pad bearings also indicated the frequency dependence of dynamic coefficients [15]. Yang advanced the partial derivative method for the prediction of a complete set of dynamic coefficients of tilting pad gas bearings [16]. Though this method successfully solved the perturbation frequency effect of dynamic coefficients, it is not suitable for turbulent lubrication.

The turbulence effect on static and dynamic properties of bearings cannot be neglected in heavy loads or high rotational speeds [17]. White developed an approximate one-dimensional finite element method for Reynolds equation with turbulence correction factors based on Hirs bulk-flow theory. Then the dynamic characteristics under half speed whirl were investigated [18]. Dousti and Allaire presented a comprehensive thermohydrodynamic model for double-film floating disk thrust bearings based on Reynolds and energy equations considering turbulence effect [19]. Lv investigated heavy-duty bearings by the transient mixed-lubrication model considering both turbulence in the large film thickness area and contact of asperities near the minimum nominal film thickness [20]. Yan developed the formula for the transition region of the turbulent lubrication model [21].

However, the dynamical effect of turbulence on stiffness and damping coefficients of high-speed tilting pad bearings was not considered in the turbulent lubrication studies above. The dynamical variations of turbulence coefficients along with the shaft vibration at different frequencies were incorporated into the methods for dynamic coefficients of fix pad bearings lubricated by organic refrigerant [22] and supercritical CO₂ [23]. In this research, by introducing the pad-pivot structural perturbations to the partial derivative method in Ref. [23], it was furthered to calculate the frequency-dependent stiffness and damping coefficients of oil-lubricated tilting pad bearings with laminar and turbulent states. Both the perturbed viscosity and the perturbed turbulence coefficients were embraced in the dynamic model of tilting pad bearings and then the importance of each perturbed variable was studied.

2. Mathematical Method

2.1. The Transient Reynolds Equation

The Reynolds equation considering variable viscosity and coexistence of laminar and turbulent flow is written below.

$$\frac{1}{R^2}\frac{\partial }{\partial \theta }\left(\frac{h^3}{k_x\mu }\frac{\partial p}{\partial \theta }\right)+\frac{\partial }{\partial z}\left(\frac{h^3}{k_z\mu }\frac{\partial p}{\partial z}\right)=\frac{\omega }{2}\frac{\partial h}{\partial \theta }+\frac{\partial h}{\partial t}$$

(1)

The k_x and k_z are circumferential and axial turbulence coefficients, respectively.

$$k_x=12+a_{1}R{e}^{b_1}; k_z=12+a_{2}R{e}^{b_2}$$

(2)

where the local Reynolds number Re = ρωRh/μ. The values of constants a₁, b₁, a₂ and b₂ in the Equation (2) are determined by flow state in the lubricating film. For turbulence (Re > 2000), the a₁ = 0.0136, b₁ = 0.900 and a₂ = 0.0043, b₂ = 0.960 are given by the Ng-Pan turbulence model [24]. Otherwise, all the constants above are equal to zero and k_x = k_z = 12 corresponding to laminar flow condition.

Considering λ = 2z/L, $\overline{\mu }$ = μ/μ_a, $\overline{p}$ = ψ²p/μ_aω, ψ = C_p/R, d = 2R, $\overline{h}$ = h/C_p and $\overline{t}$ = ωt, the dimensionless equation was obtained.

$$\frac{\partial }{\partial \theta }\left(\frac{{\overline{h}}^3}{k_x\overline{\mu }}\frac{\partial \overline{p}}{\partial \theta }\right)+{\left(\frac{d}{L}\right)}^2\frac{\partial }{\partial \lambda }\left(\frac{{\overline{h}}^3}{k_z\overline{\mu }}\frac{\partial \overline{p}}{\partial \lambda }\right)=\frac{1}{2}\frac{\partial \overline{h}}{\partial \theta }+\frac{\partial \overline{h}}{\partial \overline{t}}$$

(3)

2.2. Static and Perturbed Film Thickness for Tilting Pad Bearings

The shaft whirling and the pads tilting around the static equilibrium position at small amplitude can be described by perturbation of eccentricity ratio, attitude angle (${\tilde{\epsilon }}_d$, ${\tilde{\phi }}_d$) and tilting angle ${\tilde{\delta }}_d$, respectively. Then the film thickness $\overline{h}$ of the lth pad of tilting pad bearing is composed of the static ${\overline{h}}_0$ between the shaft and the pads at equilibrium position and the perturbed film thickness ${\tilde{h}}_d$ [16].

$$\overline{h}={\overline{h}}_0+{\tilde{h}}_d{e}^{i\mathsf{\Omega }\overline{t}}$$

(4)

$${\overline{h}}_0=1-m\cos \left({\beta }_l-\theta \right)+{\epsilon }_0\cos \left(\theta -{\phi }_0\right)+{\overline{\delta }}_{0\_l}\sin \left({\beta }_l-\theta \right)$$

(5)

$${\tilde{h}}_d=\stackrel{shaft whirling}{\overbrace{{\tilde{\epsilon }}_d\cos \left(\theta -{\phi }_0\right)+{\epsilon }_0{\tilde{\phi }}_d\sin \left(\theta -{\phi }_0\right)}}+\stackrel{pads tilting}{\overbrace{{\tilde{\delta }}_{d\_l}\sin \left({\beta }_l-\theta \right)}}$$

(6)

where the β_l is the circumferential angle of the pivot position of the lth pad, as shown in Figure 1, and the subscript _l in the ${\overline{\delta }}_0$, ${\tilde{\delta }}_d$ represents the lth pad as well. The preload factor m = 1 − C_b/C_p and the $\overline{\delta }=\delta /\psi$. The dimensionless perturbation frequency Ω is the ratio of shaft perturbation frequency υ to the rotating frequency ω.

2.3. Static and Perturbed Reynolds Equation

With the shaft and pads perturbation, the pressure is expressed as static pressure ${\overline{p}}_0$ and complex perturbed pressure ${\tilde{p}}_d$.

$$\overline{p}={\overline{p}}_0+{\tilde{p}}_d{e}^{i\mathsf{\Omega }\overline{t}}$$

(7)

The partial derivative method was extended in [23] embracing dynamic variations of complete variables. For the lubricating oil film, the dynamic variation of both film thickness and viscosity originate from turbulence coefficients k_x, k_z as well as ${\overline{h}}^3/\overline{\mu }$. The linearized complex perturbation expansions of turbulence coefficients are combined responses to perturbed film thickness ${\tilde{h}}_d$ and perturbed viscosity ${\tilde{\mu }}_d$.

$$\begin{gathered} \frac{1}{k_x}=\frac{1}{k_{x0}}\left(1+b_1\left(\frac{{\tilde{\mu }}_d}{{\overline{\mu }}_0}-\frac{{\tilde{h}}_d}{{\overline{h}}_0}\right)e^{i\mathsf{\Omega }\overline{t}}\right) \\ \frac{1}{k_z}=\frac{1}{k_{z0}}\left(1+b_2\left(\frac{{\tilde{\mu }}_d}{{\overline{\mu }}_0}-\frac{{\tilde{h}}_d}{{\overline{h}}_0}\right)e^{i\mathsf{\Omega }\overline{t}}\right) \end{gathered}$$

(8)

where k_x0 and k_z0 are static turbulence coefficients at the equilibrium position. When the Re < 2000, the k_x0 = k_z0 = 12 for the laminar state and the corresponding b₁ and b₂ in Equation (8) are equal to zero.

$$\begin{gathered} k_{x0}=12+a_1{\left(R{e}_c\frac{{\overline{h}}_0}{{\overline{\mu }}_0}\right)}^{b_1} \\ {k}_{z0}=12+a_2{\left(R{e}_c\frac{{\overline{h}}_0}{{\overline{\mu }}_0}\right)}^{b_2} \end{gathered}$$

(9)

where the Re_c = ρωRC_p/μ_a is the average Reynolds number.

The relation between perturbed viscosity ${\tilde{\mu }}_d$ and perturbed pressure ${\tilde{p}}_d$ was constructed through a dynamic mapping [23] expressed by the dimensionless partial derivative of viscosity to pressure at ${\overline{p}}_0$.

$$\overline{\mu }={\overline{\mu }}_0+{\tilde{\mu }}_d{e}^{i\mathsf{\Omega }\overline{t}}={\overline{\mu }}_0+{\overline{\partial }}_{\overline{p}}\overline{\mu }\left(p_0,T\right){\tilde{p}}_d{e}^{i\mathsf{\Omega }\overline{t}}$$

(10)

Taking Equations (4), (7), (8) and (10) into Equation (3), the static and perturbed Reynolds equation is separated, which method of derivation and solution was verified in [23].

$$\frac{\partial }{\partial \theta }\left(\frac{{\overline{h}}_0^3}{k_{x0}{\overline{\mu }}_0}\frac{\partial {\overline{p}}_0}{\partial \theta }\right)+{\left(\frac{d}{L}\right)}^2\frac{\partial }{\partial \lambda }\left(\frac{{\overline{h}}_0^3}{k_{z0}{\overline{\mu }}_0}\frac{\partial {\overline{p}}_0}{\partial \lambda }\right)=\frac{1}{2}\frac{\partial {\overline{h}}_0}{\partial \theta }$$

(11)

$$\begin{array}{ll}\frac{\partial }{\partial \theta }\left(\frac{{\overline{h}}_0^3}{k_{x0}{\overline{\mu }}_0}\frac{\partial {\tilde{p}}_d}{\partial \theta }\right)& +{\left(\frac{d}{L}\right)}^2\frac{\partial }{\partial \lambda }\left(\frac{{\overline{h}}_0^3}{k_{z0}{\overline{\mu }}_0}\frac{\partial {\tilde{p}}_d}{\partial \lambda }\right)\hfill \\ & +\left(b_1-1\right)\frac{\partial }{\partial \theta }\left(\frac{{\overline{\partial }}_{\overline{p}}\overline{\mu }\left(p_0,T\right)}{{\overline{\mu }}_0}\frac{{\overline{h}}_0^3}{k_{x0}{\overline{\mu }}_0}\frac{\partial {\overline{p}}_0}{\partial \theta }{\tilde{p}}_d\right)\hfill \\ & +\left(b_2-1\right){\left(\frac{d}{L}\right)}^2\frac{\partial }{\partial \lambda }\left(\frac{{\overline{\partial }}_{\overline{p}}\overline{\mu }\left(p_0,T\right)}{{\overline{\mu }}_0}\frac{{\overline{h}}_0^3}{k_{z0}{\overline{\mu }}_0}\frac{\partial {\overline{p}}_0}{\partial \lambda }{\tilde{p}}_d\right)\hfill \\ & =\left(b_1-3\right)\frac{\partial }{\partial \theta }\left(\frac{{\overline{h}}_0^2}{k_{x0}{\overline{\mu }}_0}\frac{\partial {\overline{p}}_0}{\partial \theta }{\tilde{h}}_d\right)\hfill \\ & +\left(b_2-3\right){\left(\frac{d}{L}\right)}^2\frac{\partial }{\partial \lambda }\left(\frac{{\overline{h}}_0^2}{k_{z0}{\overline{\mu }}_0}\frac{\partial {\overline{p}}_0}{\partial \lambda }{\tilde{h}}_d\right)+\frac{1}{2}\frac{\partial {\tilde{h}}_d}{\partial \theta }+i\mathsf{\Omega }{\tilde{h}}_d\hfill \end{array}$$

(12)

The static Reynolds equation Equation (11) governs the static pressure ${\overline{p}}_0$. The perturbed equation Equation (12) is for solving the dynamic coefficients by means of a partial derivative method [16,23]. The boundary conditions of Equations (11) and (12) for oil-lubricated tilting pad bearings are shown in Equations (13)–(15). The oil is not able to bear the tensile stress. For both static and perturbed equations, the ambient pressure at two axial ends and the leading edge, and trailing edge of the pads is set at zero.

$${\overline{p}}_0\left(\theta ,\lambda =\pm 1\right)={\tilde{p}}_d\left(\theta ,\lambda =\pm 1\right)=0$$

(13)

$$\begin{array}{l}{\overline{p}}_0\left(\theta ={\theta }_L,\lambda \right)={\overline{p}}_0\left(\theta ={\theta }_T,\lambda \right)=0\\ {\tilde{p}}_d\left(\theta ={\theta }_L,\lambda \right)={\tilde{p}}_d\left(\theta ={\theta }_T,\lambda \right)=0\end{array}$$

(14)

When a trend of negative pressure occurs in the pads of bearing, the oil film breaks and the cavitation occurs. It is governed by the Reynolds boundary condition in this research shown as Equation (15). The static and perturbed pressure become zero at the fracture area of the oil film.

$$\begin{array}{ll}\begin{array}{l}{\overline{p}}_0\left(\theta ,\lambda \right)=0\\ {\tilde{p}}_d\left(\theta ,\lambda \right)=0\end{array}& if\end{array} {\overline{p}}_{0\_trend}\left(\theta ,\lambda \right)<0$$

(15)

2.4. The Dynamic Coefficients of Tilting Pad Bearings

For tilting pad bearings, without considering the radial displacement of pivot, the perturbation variables ${\tilde{\epsilon }}_d$, ${\tilde{\phi }}_d$ and ${\tilde{\delta }}_d$ are involved in perturbed equation Equation (12). Taking partial derivative of Equation (12) to ${\tilde{\epsilon }}_d$, ${\epsilon }_0{\tilde{\phi }}_d$ and ${\tilde{\delta }}_d$ yields three partial differential equations for $P_{\epsilon }$, $P_{\phi }$ and $P_{\delta }$ instead of ${\tilde{p}}_d$. Meanwhile, according to Equation (6) the corresponding ${\tilde{h}}_d$ is replaced by $H_{\epsilon }$ = $\cos \left(\theta -{\phi }_0\right)$, $H_{\phi }$ = $\sin \left(\theta -{\phi }_0\right)$ and $H_{\delta }$ = $\sin \left({\beta }_l-\theta \right)$, respectively. Then the numerical integration below is suitable for the coordinate system in Figure 1.

$$\begin{gathered} -\frac{1}{2}{{\displaystyle \int }}_{-1}^1{{\displaystyle \int }}_{{\theta }_L}^{{\theta }_T}{P}_{•\_l}\cos \theta d\theta d\lambda =k_{x•\_l}+i\mathsf{\Omega }{c}_{x•\_l} \\ -\frac{1}{2}{{\displaystyle \int }}_{-1}^1{{\displaystyle \int }}_{{\theta }_L}^{{\theta }_T}{P}_{•\_l}\sin \theta d\theta d\lambda =k_{y•\_l}+i\mathsf{\Omega }{c}_{y•\_l} \end{gathered}$$

(16)

where the subscript ∙ represents $\epsilon$, $\phi$, and $\delta$. The leading edge of the lth pad θ_L has a relation with pivot offset γ, pad’s angular extent θ_pad, and β_l as θ_L = β_l − γθ_pad. The relation between the lth pad’s trailing edge θ_T and its leading edge θ_L is θ_T = θ_L + θ_pad.

Equation (17) is employed to get dimensionless dynamic coefficients caused by shaft whirling in the coordinate system in Figure 1.

$$\begin{gathered} \left[\begin{array}{l}{k}_{xx\_l}\\ k_{xy\_l}\end{array}\right]=A_T\left[\begin{array}{l}{k}_{x\epsilon \_l}\\ k_{x\phi \_l}\end{array}\right]=\left[\begin{array}{ll}\cos {\phi }_0& -\sin {\phi }_0\\ \sin {\phi }_0& \cos {\phi }_0\end{array}\right]\left[\begin{array}{l}{k}_{x\epsilon \_l}\\ k_{x\phi \_l}\end{array}\right] \\ \left[\begin{array}{l}{k}_{yx\_l}\\ k_{yy\_l}\end{array}\right]=A_T\left[\begin{array}{l}{k}_{y\epsilon \_l}\\ k_{y\phi \_l}\end{array}\right];\left[\begin{array}{l}{c}_{xx\_l}\\ c_{xy\_l}\end{array}\right]=A_T\left[\begin{array}{l}{c}_{x\epsilon \_l}\\ c_{x\phi \_l}\end{array}\right];\left[\begin{array}{l}{c}_{yx\_l}\\ c_{yy\_l}\end{array}\right]=A_T\left[\begin{array}{l}{c}_{y\epsilon \_l}\\ c_{y\phi \_l}\end{array}\right] \end{gathered}$$

(17)

The dimensionless dynamic equation of the pad corresponding to moment to the pivot is

$$\overline{J}{\ddot{\overline{\delta }}}_{d\_l}={\overline{F}}_{y\_l}\cos {\beta }_l-{\overline{F}}_{x\_l}\sin {\beta }_l$$

(18)

The subsequent process for deriving the dimensionless equivalent dynamic coefficients of the tilting pad bearing $K_{ij}$ and $C_{ij}$ (I, j = x, y) is similar to [16]. The dimensionless definition form is taken as [25].

$$\begin{array}{l}{K}_{ij}=\frac{{\mu }_a\omega RL}{W{\psi }^2}{\displaystyle {\displaystyle \sum }_{l=1}^5}{k}_{ij\_l}-p\left(k_{i\delta \_l}{k}_{\delta j\_l}-{\mathsf{\Omega }}^2{c}_{i\delta \_l}{c}_{\delta j\_l}\right)-{\mathsf{\Omega }}^{2}q\left(k_{i\delta \_l}{c}_{\delta j\_l}+k_{\delta j\_l}{c}_{i\delta \_l}\right)\\ C_{ij}=\frac{{\mu }_a\omega RL}{W{\psi }^2}{\displaystyle {\displaystyle \sum }_{l=1}^5}{c}_{ij\_l}-p\left(k_{i\delta \_l}{c}_{\delta j\_l}+k_{\delta j\_l}{c}_{i\delta \_l}\right)-{\mathsf{\Omega }}^{2}q\left(c_{i\delta \_l}{c}_{\delta j\_l}-\frac{k_{i\delta \_l}{k}_{\delta j\_l}}{{\mathsf{\Omega }}^2}\right)\end{array}$$

(19)

where

$$p=\frac{k_{\delta \delta \_l}-\overline{J}{\mathsf{\Omega }}^2}{{\left(k_{\delta \delta \_l}-\overline{J}{\mathsf{\Omega }}^2\right)}^2+{\left(\mathsf{\Omega }{c}_{\delta \delta \_l}\right)}^2}; q=\frac{c_{\delta \delta \_l}}{{\left(k_{\delta \delta \_l}-\overline{J}{\mathsf{\Omega }}^2\right)}^2+{\left(\mathsf{\Omega }{c}_{\delta \delta \_l}\right)}^2}$$

(20)

The pad-shaft coupling dynamic coefficients of the lth pad in Equations (19) and (20) can be obtained according to the dynamic equation of pad moment to the pivot.

$$\begin{gathered} \begin{array}{l}{k}_{\delta x\_l}=k_{xx\_l}\sin {\beta }_l-k_{yx\_l}\cos {\beta }_l\\ c_{\delta x\_l}=c_{xx\_l}\sin {\beta }_l-c_{yx\_l}\cos {\beta }_l\end{array} \\ \begin{array}{l}{k}_{\delta y\_l}=k_{xy\_l}\sin {\beta }_l-k_{yy\_l}\cos {\beta }_l\\ c_{\delta y\_l}=c_{xy\_l}\sin {\beta }_l-c_{yy\_l}\cos {\beta }_l\end{array} \\ \begin{array}{l}{k}_{\delta \delta \_l}=k_{x\delta \_l}\sin {\beta }_l-k_{y\delta \_l}\cos {\beta }_l\\ c_{\delta \delta \_l}=c_{x\delta \_l}\sin {\beta }_l-c_{y\delta \_l}\cos {\beta }_l\end{array} \end{gathered}$$

(21)

The dimensional dynamic coefficients can be converted as follows (taking the x direction as an example).

$$KXX=\frac{W}{C_p}{K}_{xx}; CXX=\frac{W}{\omega C_p}{C}_{xx}$$

(22)

2.5. Numerical Method and Mesh Discretization

The static equation Equation (11) was discretized into a five-point difference scheme and solved by the iterative method. The static viscosity should be updated along with every iteration step of static pressure based on the Barus pressure–viscosity relationship [26] in this research. The local Reynolds number at each node need to be calculated for identifying the flow state. At the turbulence zone of lubricating film, the coefficients k_x0 and k_z0 should also be calculated during the iteration of pressure ${\overline{p}}_0$. If the static pressure ${\overline{p}}_0$ lower than zero occurs at a certain node during the iteration, its value would be set as zero. When the static pressure ${\overline{p}}_0$ reaches converge judge Equation (23), the Reynolds boundary condition Equation (15) for ${\overline{p}}_0$ is satisfied.

$$\frac{{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n\left|{\overline{p}}_{0_{\left(i,j\right)}}^{new}-{\overline{p}}_{0_{\left(i,j\right)}}^{old}\right|}}}{{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n\left|{\overline{p}}_{0_{\left(i,j\right)}}^{old}\right|}}}\le {10}^{-5}$$

(23)

The solution of the static equation is repeated with the iteration of the static tilting angle ${\delta }_0$ of each pad.

$$\overline{M}={{\displaystyle \int }}_{-1}^1{{\displaystyle \int }}_{{\theta }_L}^{{\theta }_T}\left({\overline{p}}_0-1\right)\sin \left({\beta }_l-\theta \right)d\theta d\lambda$$

(24)

When the dimensionless pivoted moment $\overline{M}$ shown in Equation (24) is lower than 10⁻⁵, the static tilting angle ${\delta }_0$ of the lth pad at its static equilibrium position is obtained.

All the variables solved from the static function Equation (11) are essential for the perturbed Reynolds equation Equation (12). Then the complex pressure distribution $P_{\epsilon }$, $P_{\phi }$ and $P_{\delta }$ were also solved numerically by finite difference method with over-relaxation iterative. At the node where the static pressure ${\overline{p}}_0$ equals zero, the values of the $P_{\epsilon }$, $P_{\phi }$ and $P_{\delta }$ should be set as zero for satisfying the Reynolds boundary condition Equation (15). The converge judges for the $P_{\epsilon }$, $P_{\phi }$ and $P_{\delta }$ are given by Equation (25).

$$\frac{{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n\left|P_{{•}_{\left(i,j\right)}}^{new}-P_{{•}_{\left(i,j\right)}}^{old}\right|}}}{{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^n\left|P_{{•}_{\left(i,j\right)}}^{old}\right|}}}\le {10}^{-5}$$

(25)

where the subscript ∙ represents $\epsilon$, $\phi$, and $\delta$.

The m and n in Equations (23) and (25) are grid numbers in the circumferential and axial direction, respectively. When the number of grids of each pad is larger than 32 × 26, the grid independence is satisfied. That is to say, the difference between results with grids of 32 × 26 and 64 × 52 is less than 2%. In the same way, the results with residuals of 10⁻⁵ and 10⁻⁶ have a small difference.

3. The Validation of Program

The program for the static characteristic of tilting pad bearing by authors was verified in Ref. [9]. For validating the dynamic characteristic program in this research, the equivalent dynamic coefficients of an oil-lubricated tilting pad bearing with the same parameters to No. 53 calculation data in the Ref. [25] were calculated. Regarding the viscosity as constant in both static and perturbed equations and the dimensionless perturbation frequency $\mathsf{\Omega }$ = 1, the results are compared with [25] and shown together in

Table 1. The validation of the dynamic characteristics program for the tilting pad bearings.

ε’	Constant μ and without kx, kz				Relative Deviation
	Kxx	Kyy	Cxx	Cyy	Kxx	Kyy	Cxx	Cyy
0.1	14.3	13.7	18.8	18.4	0.69%	1.44%	1.05%	1.08%
0.2	8.18	6.82	9.88	8.98	1.68%	0.44%	1.20%	1.43%
0.3	6.74	4.46	7.13	5.68	1.17%	2.19%	1.66%	1.73%
0.4	6.62	3.18	5.95	3.90	0.60%	2.15%	1.65%	1.76%
0.5	7.25	2.31	5.40	2.71	0.14%	1.70%	0.92%	1.81%
0.6	8.63	1.62	5.23	1.82	0.12%	2.41%	0.76%	2.15%
0.7	11.1	1.06	5.48	1.14	0.28%	1.85%	3.52%	2.56%
0.8	15.8	0.593	6.46	0.608	0.63%	2.79%	0.78%	3.18%
0.85	20.3	0.399	7.11	0.400	2.87%	2.44%	3.49%	2.91%

. The difference in results is less than 5%, which verified the correctness of the dynamic coefficients program with the partial derivative method. The ε′ represents the assembled eccentricity ratio of tilting pad bearing, which has relation with ε₀ as ε₀ = ε′(1 − m).

4. Results and Discussion

A number of studies have carried out thermohydrodynamic lubrication studies on tilting pad bearings [27,28,29]. However, the isothermal method, which considers the thermal effect by average viscosity (corresponding to mean temperature) of the oil film estimated from the overall heat balance [25], is still widely used for design purposes [30,31,32]. Such is introduced to advanced rotor dynamics software such as Ref. [33]. In this paper, we focus on the influence of dynamical variation of turbulence coefficients and viscosity on stiffness and damping coefficients. Thus, the lubricating film is considered isothermal. For the same reason, the impact of inertia of the pad is not considered. So the cross-coupling equivalent stiffness and damping coefficients of tilting pad bearing are close to zero.

The structure and operating parameters of the tilting pad bearing are shown in

Table 2. Parameters of tilting pad bearing.

Parameters	Value
Bearing radius (R)	70 mm
Bearing length (L)	56 mm
Pad arc angle (α)	60°
Preload factor (m)	0.5
Pivot offset (γ)	0.5
Density of lubrication oil (ρ)	871 kg/m3
Average viscosity of oil film (μa)	9.37 mPa·s

unless otherwise stated, and the assembled radius clearance C_b has values of 109.5 and 143.5 μm which were taken from a rotor system of a certain modeled gas turbine test rig.

4.1. Analysis of Coexistence of Laminar and Turbulent Flow

Figure 2 shows the changes in the static tilting angle of each pad (in degree) of tilting pad bearing with assembled eccentricity ratio ε′ at a rotating speed of 18,000 r/min. With the increase of ε′, the tilting angle of both the 5# pad and the 4# loaded pad increases, and that of the 2# pad and 3# loaded pad decreases. The results obtained by the model considering turbulence and the pure laminar model were compared. The corresponding distributions of the local Reynolds number of each pad are shown in Figure 3. The larger the ε′ is, the lower the local Reynolds number Re of oil film of 3# and 4# loaded pads are, and the higher the Re of other pads are. The reason is that a larger ε′ makes a closer distance between the shaft and loaded pads, which means a smaller oil film thickness of loaded pads.

For C_b = 109.5 μm, the flow state in all the pads is laminar at ε′ = 0.1. When the ε′ increases, turbulence appears in the upper reaches of the 1# top pad (see Figure 3a), which makes the tilting angle of the 1# pad larger (see Figure 2a) because the turbulence in the upper reaches increases the oil film force. However, the difference between the tilting angle of the 1# top pad obtained by the model considering turbulence and the pure laminar model increases first and then decreases at ε′ beyond 0.5 because the proportion of turbulence zone of the 1# pad becomes larger than pivot offset and the turbulence existed also in pivot downstream. When the ε′ increases to 0.7, flow in the whole oil film of the 1# pad is turbulent, so the influence of turbulence on its tilting angle becomes weak. The flow state in the 3# and 4# loaded pads is laminar at arbitrary ε′ for C_b = 109.5 μm. The tilting angle of both the 2# and the 5# pad obtained by the model considering turbulence becomes larger than that by the pure laminar model when ε′ is over 0.6 because turbulence appears in the upper reaches of the oil film of 2# and 5# pad.

The relation between the tilting angle of the five pads for C_b = 143.5 μm has some difference to smaller assembled radius clearance when the coexistence of laminar and turbulent flow is considered. The turbulence in the 2# pad makes its static tilting angle close to that of the 3# pad. The proportion of turbulence zone of the upper reaches of both the 2# and the 5# pad is about equal to pivot offset at arbitrary ε′ (see Figure 3b), so the turbulence effect increases their tilting angle obviously (see Figure 2b). The flow state in all the pads contains laminar and turbulence together at ε′ = 0.1, 0.2, and 0.3. When the ε′ increases to 0.4, the flow state in the 3# and 4# loaded pads become pure laminar flow. The influence of turbulence on the tilting angle of the 3# and 4# pads is also not very obvious at small ε′ because their lubricating film is thin. The turbulence exists in the 1#, 2#, and 5# pads at arbitrary ε′. The influence of turbulence on the tilting angle of the 1# pad is weak at ε′ larger than 0.4, the reason for which is the same as that of C_b = 109.5 μm at ε′ larger than 0.7.

4.2. The Influence of the Perturbed Viscosity on Dynamic Coefficients

Figure 4 illustrates the dimensionless equivalent dynamic coefficients with the change of ε′ for different assembled clearance C_b. The difference between results obtained by the model considering turbulence and the pure laminar model is compared, and the influence of perturbed viscosity was also investigated.

For either the model considering turbulence or the pure laminar model, each dimensionless stiffness or damping coefficient considering perturbed viscosity is larger than that without perturbed viscosity. Especially the stiffness coefficients and damping coefficients C_xx for C_b = 109.5 μm at ε′ larger than 0.7 (see Figure 4a,c,e), the influence of perturbed viscosity is notably obvious. A reason is that the smaller oil film thickness of the 3# and 4# loaded pads enhances the pressure-viscous effect.

For C_b = 109.5 μm, with the increase of ε′, the difference between results obtained by the model considering turbulence and the pure laminar model increases first and then declines after about ε′ = 0.4. Such is corresponding to the local Reynolds number (see Figure 3a) and the turbulence influence on the tilting angle of the 1# top pad (see Figure 2a).

For C_b = 143.5 μm, both the dimensionless equivalent stiffness coefficients and the damping coefficients considering turbulence are prominently larger than that obtained by the pure laminar model at ε′ = 0.1, 0.2, and 0.3, which reason is that the upper reaches of the oil film of 3# and 4# pads are in a turbulent state. When the ε′ is larger than 0.4, the difference between results obtained by the two models becomes small and decreases gently with ε′. A reason is that the flow state in the whole oil film of the 1# top pad is turbulent (see Figure 3b) and the influence of turbulence on the tilting angle of the 1# pad is weak (see Figure 2b). The turbulence has nearly no effect on dimensionless dynamic coefficients at large ε′ though the oil film of 2# and 5# pads contains both laminar and turbulent flow because the oil film thickness of the 3# and 4# loaded pads is small and the viscous effect is dominant.

4.3. The Influence of Perturbed Film Thickness Included in Turbulence Coefficients

For assembled radius clearance C_b = 143.5 μm and at assembled eccentricity ratio ε′ = 0.1, the variations of dimensionless equivalent dynamic coefficients of a tilting pad bearing with rotating speed under dimensionless perturbation frequency Ω = 0.5 and 1 are shown in Figure 5. Various kinds of solving situations were presented together for reflecting the contributions of perturbed film thickness included in turbulence coefficients perturbation. The results shown in red and black lines were obtained by the partial derivative method embracing turbulence perturbation. Though the turbulence was considered in the results shown in the blue line, the influence of dynamical variation of turbulence coefficients is not considered.

Comparison of red (or black) line to blue line indicates that neglecting the perturbed film thickness ${\tilde{h}}_d$ included in turbulence coefficients perturbation k_x, k_z leads to obviously excessive results of equivalent stiffness coefficients at rotating speed beyond a certain value. Under different dimensionless perturbation frequencies Ω, the regulation of the influence of turbulence perturbation on a certain stiffness or damping coefficient is similar.

The dimensionless equivalent stiffness coefficients under dimensionless perturbation frequency Ω = 0.5 are larger than that under Ω = 1. However, the dimensionless equivalent damping coefficients present the opposite relationship. Under different Ω, the variation regulation of each dynamic coefficient with rotating speed is consistent.

5. Conclusions

For studying the dynamical effect of turbulence on frequency-dependent stiffness and damping coefficients of high-speed oil-lubricated tilting pad bearings, the partial derivative method was used for perturbing the turbulent lubrication model. In this research, the perturbations of both viscosity and turbulence coefficients were embraced in the dynamic model of tilting pad bearings. The importance of each perturbed variable was studied, and the main conclusions from the results are as follows:

(1)

The influence of turbulence on dimensionless dynamic coefficients is obvious at a medium assembled eccentricity ratio (about 0.3 to 0.5) for small assembled radius clearance, and at a small assembled eccentricity ratio (about 0.1 to 0.3) for large assembled radius clearance. The perturbed viscosity takes effect mainly at the large assembled eccentricity ratio for small assembled radius clearance.

(2)

For large assembled radius clearance and at high rotating speed, neglecting the perturbed film thickness included in turbulence coefficients perturbations will lead to obviously excessive results of equivalent stiffness coefficients. Thus, the dynamic effect of turbulence needs to be considered.

(3)

The perturbed viscosity makes the stiffness and damping coefficients larger at both laminar and turbulent flow states, so the consideration of the perturbed viscosity is optional. For a simplified simulation and conservative prediction results, the perturbed viscosity can be neglected.