Theoretical Evaluation on Lubrication Performance of Thrust-Type Foil Bearings in Liquid Nitrogen

1. Introduction

The recovery and reutilization of liquid-propellant rockets can significantly reduce the launch cost and play a critical role in promoting large-scale high-frequency space activities [1,2]. Whether the high-speed turbopumps, the heart part of liquid rockets, can truly achieve maintenance-free reutilization determines the actual launch cost [3]. Taking the space shuttle as an example, NASA initially aimed to control the cost under 30 million dollars per launch. However, due to the long-duration complicated maintenance of the liquid rocket turbopumps, the actual cost of space shuttle per launch is as high as 500 million dollars [4]. Till now, all the launched liquid-propellant rocket turbopumps have utilized rolling element bearings to support their high-speed shaft systems, which transmit external loads by point or line contacting form [5,6]. Due to the limitation of wear resistance, it is challenging to realize a long service life for the rolling bearings under high rotational speed. The reported maximum DN value (the product of inner diameter and rotational speed) of rolling bearings in rocket turbopumps is around 3×10⁶ mm·r·min^-1 [7]. To overcome the contradiction between service life and rotational speed in rolling bearings, major space powers have explored the application feasibility of fluid-film bearings to improve high-speed performance. The pressurized film inside fluid-film bearings, usually generated by hydrodynamic or hydrostatic effect, separate the rotor journal from the bearing surface, which can significantly decrease the wear degree [8,9,10,11,12]. As early as 1983, Hanum et al. [13] from NASA (National Aeronautics and Space Administration) designed a combination bearing system with two angular contact rolling bearings and a hydrostatic bearing. They tested the combination bearing in the MK48 turbopump rotor system and found that the wear inside the rolling bearings was slight even after 2337 start-ups and running at a maximum rotational speed of 70000 r·min^-1. The introduction of hydrostatic fluid-film bearing could significantly extend the wear life of rolling bearings. Since then, researchers have begun exploring the feasibility of completely abandoning rolling bearings and using fluid-film bearings only. Ohta et al. [14] designed and tested a hydrostatic bearing with ten cavities for LE-5 liquid hydrogen (LH₂) turbopumps. The inner diameter of the bearing is 52.1 mm and the average radial clearance is 50 μm. During the start-up stage, the rotor journal could separate from the bearing surface through the hydrostatic pressure of LH₂. When the rotational speed reached 50000 r·min^-1 and the supply of high-pressure liquid hydrogen was stopped, the rotor can also realize a stable running state through the generated hydrodynamic pressure. Around 2010, French national space agency and Swedish national space board implemented the TPX full-size LH₂-turbopump demonstration program and hydrostatic fluid-film bearings were employed for frequently tests [15]. The maximum tested rotational speed is 39500 r·min^-1.

However, the severe friction and wear during the start-up and stopping stages of hydrodynamic fluid-film bearings cannot be ignored, and it is almost impossible to add a high-pressure fluid source in rocket turbopumps for hydrostatic fluid-film bearings [16]. Therefore, new bearing schemes have been proposed and investigated successively. Xu et al. [17] proposed the scheme of superconducting compound bearings, which combine superconducting magnetic field and hydrodynamic fluid field, for reusable rocket turbopumps. The superconducting compound bearings can simultaneously ensure friction-free operation during the start-up and stopping stages through superconducting magnetic field, and provide high stability during the steady working stage through the hydrodynamic fluid field [18]. Specially, due to the outstanding reliability of gas foil bearings in turbomachinery, NASA evaluated the performance of foil bearings in cryogenic fluids and the feasibility for rocket application [19]. The mean time between failures of gas foil bearings in aircraft turbo-compressors can exceed 60000 hours [20]. The demonstration experiments of foil bearings in LH₂ and LO₂ (liquid oxygen) rocket turbopumps [21] were successfully conducted in 1992 and 1993, respectively. The maximum test speeds reached 91000 r·min^-1 and 25000 r·min^-1, respectively. After over 100 times of start-up/stopping cycles, the wear of the foil bearings was acceptable for continuing experiments.

Generally, the existed theoretical investigations about the foil bearings, including journal and thrust bearings, mainly concentrate on the working environment of gas lubrication. Since thrust foil bearings are required to adapt the variation of axial loads, it is necessary to reveal the interaction mechanism between journal and the bearing surface [22]. Shi et al. [23] established an analytical model of aerostatic thrust bearing with 3-DOF considering the interaction among perturbation Reynolds equations corresponding to the stiffness and damping. Xu et al. [24] proposed an exhaustive model of gas foil thrust bearings combining isoparametric elements with contact mechanics, which enables the analysis of arbitrarily shaped foils. Chen et al. [25] established a dynamic model of a foil thrust bearing with an angular swing pad and the influence of angular disturbance on the loading capacity was analyzed. Feng et al. [26] employed link-spring structural model to calculate the equivalent vertical stiffness of bump-type gas foil thrust bearing. For engineering applications, it is important to ensure the load capacity of thrust foil bearings. Somaya et al. [27] evaluated the performance of an aerodynamic foil thrust bearings at a high speed of 350,000 r·min^-1 and the load capacity coefficient was improved to 5.36×10^-6 N/(mm³·kr·min^-1). Besides the increase in the operational speed, structural modifications are more frequently adopted to improve the load capacity, including the hybrid bearing combining air foil bearing and hydrostatic air bearing [28], multi-layer foil structure [29,30,31], and Ryleigh step air foil structure [32]. These investigations provide a well theoretical foundation for the analysis of thrust foil bearings in various lubrication conditions.

In this study, with the application prospects in liquid rocket turbopumps, theoretical evaluations on lubrication performance of thrust-type foil bearings in liquid nitrogen were conducted. Liquid nitrogen was adopted to simulate the cryogenic fluid in rocket turbopumps for safety and the further experimental verifications in the Lab. The static and dynamic characteristics of the bearings were analyzed and the application feasibility in rocket turbopumps were discussed.

2. Theoretical Analysis Procedures

2.1. Fluid Lubrication Equations

Consider liquid nitrogen as an ideal incompressible medium. Using the method illustrated in Figure 1 to establish a cylindrical coordinate system, and based on the continuity equation and the Navier-Stokes equations, the dimensionless form of the 2D incompressible Reynolds equation applicable to foil thrust bearings can be derived as follows [33]:

$$\{\begin{array}{l}\frac{1}{\overline{r}}\frac{\partial }{\partial \overline{r}}\overline{r}{H}^3\frac{\partial P}{\partial \overline{r}}+\frac{1}{{\overline{r}}^2}\frac{\partial }{\partial \theta }{H}^3\frac{\partial P}{\partial \theta }=\mathrm{\Lambda }\frac{\partial H}{\partial \theta }+2\mathrm{\Lambda }\gamma \frac{\partial H}{\partial \overline{t}}\\ H=\frac{h}{h_2},\overline{r}=\frac{r}{r_2},P=\frac{p}{p_s},\overline{t}=w_{e}t,\mathrm{\Lambda }=\frac{6{\mu }_0\omega }{P_s}{\left(\frac{r_2}{h_2}\right)}^2,\gamma =\frac{{\omega }_e}{\omega }\end{array}$$

(1)

Where, $h_2$ represents the minimum initial fluid film thickness, $r_2$ represents the outer diameter of the thrust pad, $H$ represents the dimensionless fluid film thickness, $\theta$ represents the circumferential coordinate, $\overline{r}$ represents the dimensionless axial coordinate, $p_s$ represents the supply pressure of liquid nitrogen, ${\mu }_0$ represents the dynamic viscosity, ${\omega }_e$ represents the external excitation frequency, and $\gamma$ represents the ratio of external excitation frequency to shaft rotational frequency.

During actual installation, rotor eccentricity is an unavoidable fault in rotating machinery. This eccentricity changes the distance between the surface of the thrust foil bearing and the rotor disk, as shown in Figure 2. The thrust pads are numbered counterclockwise, and due to the inclination of the thrust disk relative to the bearing, the initial fluid film thickness distribution on different thrust pads becomes uneven. This uneven distribution affects the performance of the foil thrust bearing. Therefore, it is necessary to consider the impact of rotor misalignment on the performance of foil journal bearings.

The fluid film thickness of a foil thrust bearing is composed of the initial gap $h_2+g\left( r,\theta \right)$ between the rotor and the thrust disk surfaces, the rotor misalignment gap$h_{\phi }\left( r,\theta \right)$, and the deformation $\delta \left( r,\theta \right)$of the bearing surface. Using the cylindrical coordinate system shown in Figure 1, the dimensionless expression for fluid film thickness, considering foil deformation and rotor misalignment, can be derived as follows [34]:

$$H=1+\overline{g}\left(\overline{r},\theta \right)+\overline{\delta }\left(\overline{r},\theta \right)+{\overline{h}}_{\phi }\left(\overline{r},\theta \right)$$

(2)

The initial wedge gap is:

$$\overline{g}\left( \overline{r},\theta \right)=\{\begin{array}{c}\left(\frac{h_1}{h_2}-1\right)\left(1-\frac{\theta }{b\beta }\right)0<\theta <b\beta \\ 0b\beta <\theta <\beta \end{array}$$

(3)

The gap caused by rotor misalignment is:

$$\{\begin{array}{l}{\overline{h}}_{\phi }\left( \overline{r},\theta \right)=r_2\left({\overline{\phi }}_y\overline{r}\cos \left(\theta \right)-{\overline{\phi }}_x\overline{r}\sin \left(\theta \right)\right)\\ {\overline{\phi }}_x=\frac{{\phi }_x}{h_2} ;{\overline{\phi }}_y=\frac{{\phi }_y}{h_2}\end{array}$$

(4)

2.2. Foil Deformation Equation

To consider the effect of foil deformation on bearing performance, we simplify each convex foil into a link-spring structure [35] consisting of two rigid links and a horizontal spring, as shown in Figure 3. This model fully considers the interaction forces between adjacent wavy foils and the frictional forces on the contact surfaces, providing a comprehensive simulation of the actual working state of the wavy foils. Since the top foil has a small thickness, a thin plate finite element model with three degrees of freedom is established, considering both bending and shear stiffness. Each convex foil is equivalent to a spring, and the total stiffness matrix$\left[k_f\right]$ is obtained by adding the equivalent vertical stiffness $\left[k_v\right]$of each convex foil calculated by the link-spring model and the stiffness matrix $\left[k_t\right]$ of the top foil calculated by the finite element method at the corresponding nodes.

The deformation of the top foil can be calculated using the direct stiffness method as$\left[k_f\right]\left[\mathrm{U}\right]=\left[\mathrm{F}\right]$, where $\left[\mathrm{U}\right]$ and $\left[\mathrm{F}\right]$ are the total displacement vector and the fluid film pressure vector of the foil structure, respectively. Here, $\left[\mathrm{U}\right]={\left[{\delta }_1,{\theta }_{x1},{\theta }_{y1},\dots ,{\delta }_i,{\theta }_{xi},{\theta }_{yi},\dots \right]}^T$ represents the total displacement vector, $\left[\mathrm{F}\right]={\left[P_1,0,0,\dots ,P_i,0,0,\dots \right]}^T$, ${\delta }_i$ is the deformation of the foil, and ${\theta }_{xi},{\theta }_{yi}$ is the rotation angle about the X and Y axes.

2.3. Dynamic Characteristics Calculation

Given a displacement and velocity disturbance$\left(\Delta \overline{z},\Delta \overline{\dot{z}},\Delta {\overline{\phi }}_x,\Delta {\overline{\dot{\phi }}}_x,\Delta {\overline{\phi }}_y,\Delta {\overline{\dot{\phi }}}_y\right)$ in the axial direction of the thrust disk, the oil film pressure P, oil film thickness H, and foil deformation $\overline{\delta }$ are expanded in a first-order Taylor series about the equilibrium position$\left(\Delta \overline{z},\Delta \overline{\dot{z}},\Delta {\overline{\phi }}_x,\Delta {\overline{\dot{\phi }}}_x,\Delta {\overline{\phi }}_y,\Delta {\overline{\dot{\phi }}}_y\right)$ as follows:

$$\{\begin{array}{l}P=P_0+P_z\Delta \overline{z}+P_{\dot{z}}\Delta \overline{\dot{z}}+P_{{\phi }_x}\Delta {\overline{\phi }}_x+P_{{\dot{\phi }}_x}\Delta {\overline{\dot{\phi }}}_x+P_{{\phi }_y}\Delta {\overline{\phi }}_y+P_{{\dot{\phi }}_y}\Delta {\overline{\dot{\phi }}}_y\\ H=H_0+H_z\Delta \overline{z}+H_{\dot{z}}\Delta \overline{\dot{z}}+H_{{\phi }_x}\Delta {\overline{\phi }}_x+H_{{\dot{\phi }}_x}\Delta {\overline{\dot{\phi }}}_x+H_{{\phi }_y}\Delta {\overline{\phi }}_y+H_{{\dot{\phi }}_y}\Delta {\overline{\dot{\phi }}}_y\\ \overline{\delta }={\overline{\delta }}_0+{\overline{\delta }}_z\Delta \overline{z}+{\overline{\delta }}_{\dot{z}}\Delta \overline{\dot{z}}+{\overline{\delta }}_{{\phi }_x}\Delta {\overline{\phi }}_x+{\overline{\delta }}_{{\dot{\phi }}_x}\Delta {\overline{\dot{\phi }}}_x+{\overline{\delta }}_{{\phi }_y}\Delta {\overline{\phi }}_y+{\overline{\delta }}_{{\dot{\phi }}_y}\Delta {\overline{\dot{\phi }}}_y\end{array}$$

(5)

When the rotor undergoes small movements near the equilibrium position, the disturbed oil film pressure distribution can be obtained as follows:

$$\{\begin{array}{l}H=H_0+\Delta \overline{z}+r_2\left(\Delta {\overline{\phi }}_y\Delta \overline{r}\cos \left(\theta \right)-\Delta {\overline{\phi }}_x\Delta \overline{r}\sin \left(\theta \right)\right)+\Delta \overline{\delta }\\ H_0=1+\overline{g}+r_2\left({\overline{\phi }}_y\Delta \overline{r}\cos \left(\theta \right)-{\overline{\phi }}_x\Delta \overline{r}\sin \left(\theta \right)\right)+{\overline{\delta }}_0\\ \Delta \overline{\delta }={\overline{\delta }}_z\Delta \overline{z}+{\overline{\delta }}_{{\phi }_x}\Delta {\overline{\phi }}_x+{\overline{\delta }}_{{\phi }_y}\Delta {\overline{\phi }}_y\end{array}$$

(6)

The disturbance of foil deformation${\overline{\delta }}_z,{\overline{\delta }}_{\dot{z}},{\overline{\delta }}_{{\phi }_x},{\overline{\delta }}_{{\dot{\phi }}_x},{\overline{\delta }}_{{\phi }_y},{\overline{\delta }}_{{\dot{\phi }}_y}$ can be obtained using the direct stiffness method as follows:

$${\overline{\delta }}_{\dot{\zeta }}=P_{\dot{\zeta }}/k_f\left(\zeta =z,{\phi }_x,{\phi }_y\right)$$

(7)

Where, $k_f$ represents the total stiffness matrix of foil deformation, the equation for the disturbed oil film distribution can be obtained by combining the above equations:

$$\{\begin{array}{l}{H}_z=1+{\overline{\delta }}_z\\ H_{{\phi }_x}=-r_2\overline{r}\sin \left(\theta \right)+{\overline{\delta }}_{{\phi }_x}\\ H_{{\phi }_y}=r_2\overline{r}\cos \left(\theta \right)+{\overline{\delta }}_{{\phi }_y}\\ {\left[\begin{array}{ccc}{H}_{\dot{z}}& H_{{\dot{\phi }}_x}& H_{{\dot{\phi }}_y}\end{array}\right]}^T={\left[\begin{array}{ccc}{\overline{\delta }}_{\dot{z}}& {\overline{\delta }}_{{\dot{\phi }}_x}& {\overline{\delta }}_{{\dot{\phi }}_y}\end{array}\right]}^T\end{array}$$

(8)

Substitute into the transient Reynolds equation, neglect higher-order terms, and combine and rearrange the equations for displacement and velocity disturbances to obtain a steady-state equation for the static equilibrium position and 6 disturbance equations for the disturbance quantities$\left(\Delta \overline{z},\Delta \overline{\dot{z}},\Delta {\overline{\phi }}_x,\Delta {\overline{\dot{\phi }}}_x,\Delta {\overline{\phi }}_y,\Delta {\overline{\dot{\phi }}}_y\right)$.

The steady-state equation is as follows:

$$\frac{1}{\overline{r}}\frac{\partial }{\partial \overline{r}}\overline{r}{H_0}^3\frac{\partial P_0}{\partial \overline{r}}+\frac{1}{{\overline{r}}^2}\frac{\partial }{\partial \theta }{H_0}^3\frac{\partial P_0}{\partial \theta }=\mathrm{\Lambda }\frac{\partial H_0}{\partial \theta }$$

(9)

The equations for the displacement disturbance$\left(\Delta \overline{z},\Delta {\overline{\phi }}_x,\Delta {\overline{\phi }}_y\right)$ are:

$$\{\begin{array}{l}\frac{1}{\overline{r}}\frac{\partial }{\partial \overline{r}}\overline{r}{H_0}^3\frac{\partial P_z}{\partial \overline{r}}+3\overline{r}{H_0}^2\left(1+\frac{P_z}{k_f}\right)\frac{\partial P_0}{\partial \overline{r}}+\frac{1}{{\overline{r}}^2}\frac{\partial }{\partial \theta }{H_0}^3\frac{\partial P_z}{\partial \theta }+3\overline{r}{H_0}^2\left(1+\frac{P_z}{k_f}\right)\frac{\partial P_0}{\partial \theta }\\ =\mathrm{\Lambda }\frac{\partial \left(1+\frac{P_z}{k_f}\right)}{\partial \theta }-2\mathrm{\Lambda }\gamma \frac{P_{\dot{z}}}{k_f}\\ \frac{1}{\overline{r}}\frac{\partial }{\partial \overline{r}}\overline{r}{H_0}^3\frac{\partial P_{{\phi }_x}}{\partial \overline{r}}+3\overline{r}{H_0}^2\left(-r_2\overline{r}\sin \left(\theta \right)+\frac{P_{{\phi }_x}}{k_f}\right)\frac{\partial P_0}{\partial \overline{r}}+\frac{1}{{\overline{r}}^2}\frac{\partial }{\partial \theta }{H_0}^3\frac{\partial P_{{\phi }_x}}{\partial \theta }+3\overline{r}{H_0}^2\left(-r_2\overline{r}\sin \left(\theta \right)+\frac{P_{{\phi }_x}}{k_f}\right)\frac{\partial P_0}{\partial \theta }\\ =\mathrm{\Lambda }\frac{\partial \left(-r_2\overline{r}\sin \left(\theta \right)+\frac{P_{{\phi }_x}}{k_f}\right)}{\partial \theta }-2\mathrm{\Lambda }\gamma \frac{P_{{\dot{\phi }}_x}}{k_f}\\ \frac{1}{\overline{r}}\frac{\partial }{\partial \overline{r}}\overline{r}{H_0}^3\frac{\partial P_{{\phi }_y}}{\partial \overline{r}}+3\overline{r}{H_0}^2\left(r_2\overline{r}\cos \left(\theta \right)+\frac{P_{{\phi }_y}}{k_f}\right)\frac{\partial P_0}{\partial \overline{r}}+\frac{1}{{\overline{r}}^2}\frac{\partial }{\partial \theta }{H_0}^3\frac{\partial P_{{\phi }_y}}{\partial \theta }+3\overline{r}{H_0}^2\left(r_2\overline{r}\cos \left(\theta \right)+\frac{P_{{\phi }_y}}{k_f}\right)\frac{\partial P_0}{\partial \theta }\\ =\mathrm{\Lambda }\frac{\partial \left(r_2\overline{r}\cos \left(\theta \right)+\frac{P_{{\phi }_y}}{k_f}\right)}{\partial \theta }-2\mathrm{\Lambda }\gamma \frac{P_{{\dot{\phi }}_y}}{k_f}\end{array}$$

(10)

The equations for the velocity disturbance$\left(\Delta \overline{\dot{z}},\Delta {\overline{\dot{\phi }}}_x,\Delta {\overline{\dot{\phi }}}_y\right)$ are:

$$\{\begin{array}{l}\frac{1}{\overline{r}}\frac{\partial }{\partial \overline{r}}\overline{r}{H_0}^3\frac{\partial P_{\dot{z}}}{\partial \overline{r}}+3\overline{r}{H_0}^2\frac{P_{\dot{z}}}{k_f}\frac{\partial P_0}{\partial \overline{r}}+\frac{1}{{\overline{r}}^2}\frac{\partial }{\partial \theta }{H_0}^3\frac{\partial P_{\dot{z}}}{\partial \theta }+3\overline{r}{H_0}^2\frac{P_{\dot{z}}}{k_f}\frac{\partial P_0}{\partial \theta }=\mathrm{\Lambda }\frac{\partial \frac{P_{\dot{z}}}{k_f}}{\partial \theta }+2\mathrm{\Lambda }\gamma \left(1+\frac{P_z}{k_f}\right)\\ \frac{1}{\overline{r}}\frac{\partial }{\partial \overline{r}}\overline{r}{H_0}^3\frac{\partial P_{{\dot{\phi }}_x}}{\partial \overline{r}}+3\overline{r}{H_0}^2\frac{P_{{\dot{\phi }}_x}}{k_f}\frac{\partial P_0}{\partial \overline{r}}+\frac{1}{{\overline{r}}^2}\frac{\partial }{\partial \theta }{H_0}^3\frac{\partial P_{{\dot{\phi }}_x}}{\partial \theta }+3\overline{r}{H_0}^2\frac{P_{{\dot{\phi }}_x}}{k_f}\frac{\partial P_0}{\partial \theta }=\mathrm{\Lambda }\frac{\partial \frac{P_{{\dot{\phi }}_x}}{k_f}}{\partial \theta }+2\mathrm{\Lambda }\gamma \left(-r_2\overline{r}\sin \left(\theta \right)+\frac{P_{{\phi }_x}}{k_f}\right)\\ \frac{1}{\overline{r}}\frac{\partial }{\partial \overline{r}}\overline{r}{H_0}^3\frac{\partial P_{{\dot{\phi }}_y}}{\partial \overline{r}}+3\overline{r}{H_0}^2\frac{P_{{\dot{\phi }}_y}}{k_f}\frac{\partial P_0}{\partial \overline{r}}+\frac{1}{{\overline{r}}^2}\frac{\partial }{\partial \theta }{H_0}^3\frac{\partial P_{{\dot{\phi }}_y}}{\partial \theta }+3\overline{r}{H_0}^2\frac{P_{{\dot{\phi }}_y}}{k_f}\frac{\partial P_0}{\partial \theta }=\mathrm{\Lambda }\frac{\partial \frac{P_{{\dot{\phi }}_y}}{k_f}}{\partial \theta }+2\mathrm{\Lambda }\gamma \left(r_2\overline{r}\cos \left(\theta \right)+\frac{P_{{\phi }_y}}{k_f}\right)\end{array}$$

(11)

The four boundaries of the thrust pad of the foil thrust bearing are in contact with the external environment, and the boundary conditions can be derived as follows:

$$\begin{array}{c}\begin{array}{c}{P}_0=P_s\\ P_z=P_{{\phi }_x}=P_{{\phi }_y}=0\end{array}\\ P_{\dot{z}}=P_{{\dot{\phi }}_x}=P_{{\dot{\phi }}_y}=0\end{array}\} \mathrm{at}\overline{r}=\frac{r_1}{r_2},1\mathrm{or}\theta =\beta i\left(i=1,2,\dots ,6\right)$$

(12)

By coupling the foil deformation equation with the disturbance equation set, and using finite difference method along with Newton-Raphson iteration, the oil film pressure P₀, the oil film thickness H₀ at the static equilibrium position, and the disturbance quantities$\left(P_z,P_{\dot{z}},P_{{\phi }_x},P_{{\dot{\phi }}_x},P_{{\phi }_y},P_{{\dot{\phi }}_y}\right)$.

By integrating over the solution domain, the axial load, frictional torque, axial dynamic stiffness, and damping of the foil thrust bearing can be obtained. The integration formulas are as follows [36]:

$$F_{\mathrm{load}}=p_s{r_2}^2\underset{0}{\overset{\beta }{{\displaystyle \int }}}\underset{r_1/r_2}{\overset{1}{{\displaystyle \int }}}(P_0-P_s)\overline{r}d\overline{r}d\theta$$

(13)

$$T_c=p_s{r_2}^2{h}_2\underset{0}{\overset{\beta }{{\displaystyle \int }}}\underset{r_1/r_2}{\overset{1}{{\displaystyle \int }}}(\frac{\overline{r}H}{2}\frac{\partial P_0}{\partial \theta }+\frac{{\overline{r}}^3}{6H}\mathrm{\Lambda })\overline{r}d\overline{r}d\theta$$

(14)

$$\left[\begin{array}{ccc}{K}_{zz}& K_{z{\phi }_x}& K_{z{\phi }_y}\\ K_{{\phi }_{x}z}& K_{{\phi }_x{\phi }_x}& K_{{\phi }_x{\phi }_y}\\ K_{{\phi }_{y}z}& K_{{\phi }_y{\phi }_x}& K_{{\phi }_y{\phi }_y}\end{array}\right]=\frac{p_s{r_2}^2}{h_2}\underset{0}{\overset{\beta }{{\displaystyle \int }}}\underset{r_1/r_2}{\overset{1}{{\displaystyle \int }}}\left[\begin{array}{c}-1\\ r_2\overline{r}sin\theta \\ -r_2\overline{r}cos\theta \end{array}\right]\left[\begin{array}{ccc}{P}_z& P_{{\phi }_x}& P_{{\phi }_y}\end{array}\right]\overline{r}d\overline{r}d\theta$$

(15)

$$\left[\begin{array}{ccc}{C}_{zz}& C_{z{\phi }_x}& C_{z{\phi }_y}\\ C_{{\phi }_{x}z}& C_{{\phi }_x{\phi }_x}& C_{{\phi }_x{\phi }_y}\\ C_{{\phi }_{y}z}& C_{{\phi }_y{\phi }_x}& C_{{\phi }_y{\phi }_y}\end{array}\right]=\frac{p_s{r_2}^2}{h_2\omega \gamma }\underset{0}{\overset{\beta }{{\displaystyle \int }}}\underset{r_1/r_2}{\overset{1}{{\displaystyle \int }}}\left[\begin{array}{c}-1\\ r_2\overline{r}sin\theta \\ -r_2\overline{r}cos\theta \end{array}\right]\left[\begin{array}{ccc}{P}_z& P_{{\dot{\phi }}_x}& P_{{\dot{\phi }}_y}\end{array}\right]\overline{r}d\overline{r}d\theta$$

(16)

3.1. Structure Parameters of Foil Thrust Bearing

For the evaluation of the lubrication performance of foil thrust bearings in liquid nitrogen, this study analyzes the changes in static and dynamic characteristics under different operating and structural parameters. The structural parameters of the foil thrust bearing studied in this paper are shown in Table 1. It should be noted that smaller top foil and wave foil thicknesses were chosen to better study the effect of foil deformation on bearing performance, allowing for better comparison with traditional ball bearings. Unless otherwise specified, the predicted results in this section are for the performance evaluation of a thrust pad under the condition that the thrust disk is not tilted.

3.3. Dynamic Characteristics Analysis

The stiffness and damping coefficients of the bearing are directly affected by the tilt angle of the thrust collar. When the thrust collar is untilted, the dynamic stiffness and damping coefficients of the bearing exhibit a skew-symmetric distribution. The cross-coupling stiffness and damping coefficients ($Kz{\phi }_x,K{\phi }_{x}z,Kz{\phi }_y,K{\phi }_{y}z,Cz{\phi }_x,C{\phi }_{x}z,Cz{\phi }_y,C{\phi }_{y}z$) are both zero, the rotational cross-coupling stiffness-damping coefficients are equal ($K{\phi }_x{\phi }_y=-K{\phi }_y{\phi }_x,C{\phi }_x{\phi }_y=-C{\phi }_y{\phi }_x,$ ), and the direct rotational stiffness-damping coefficients are also equal ($K{\phi }_x{\phi }_x=K{\phi }_y{\phi }_y,C{\phi }_x{\phi }_x=C{\phi }_y{\phi }_y$). However, compared to the direct translational stiffness-damping coefficients ($K_{zz}{C}_{zz}$), the amplitude of the rotational stiffness-damping coefficients is smaller. The dynamic performance of the bearing is mainly determined by the direct stiffness-damping coefficients ($K_{zz}{C}_{zz}$) because when the tilt angle of the thrust foil bearing is zero, the film thickness distribution is symmetric with respect to the X and Y axes, resulting in very small coupling effects between translational and tilt motions. Therefore, unless otherwise specified, the dynamic analysis of thrust foil bearings is conducted assuming no tilt of the thrust collar, resulting in the calculation of the direct translational stiffness-damping coefficients of a thrust pad.

Figure 11 investigates the effect of different excitation frequency ratios on the direct translational stiffness and damping coefficients of the bearing. From the figure, it can be observed that as the excitation frequency increases, the direct translational stiffness coefficient increases, while the direct translational damping coefficient decreases. This is because the increase in excitation frequency leads to a greater squeeze effect on the lubricating film, thereby increasing the pressure distribution.

Figure 12 illustrates the effect of different friction coefficients on the direct translational stiffness and damping coefficients. It can be seen that an increase in the friction coefficient leads to an increase in the equivalent stiffness of the waves, resulting in an increase in both the direct translational stiffness and damping coefficients, approaching a constant value. Increasing the speed also increases the rate of growth of the stiffness and damping coefficients. This indicates that at high speeds, the dynamic stiffness-damping coefficients are more sensitive to the friction coefficient. This is because at higher speeds, there is a higher load-bearing capacity, forcing the foil structure to undergo larger deformations. In contrast, at lower speeds, the lower load-bearing capacity means that an increase in the friction coefficient does not cause significant deformation of the foil structure compared to the original film thickness, thus not significantly affecting the film pressure distribution.

Figure 13 shows the effect of the minimum film thickness on the direct translational stiffness and damping coefficients. It can be seen that at a certain speed, increasing the minimum film thickness reduces the static load of the bearing, leading to a continuous decrease in both the direct translational stiffness and damping coefficients, with the rate of decrease slowing down. Additionally, increasing the speed causes the direct translational stiffness coefficient to gradually increase and the direct translational damping coefficient to gradually decrease, with the growth rate of both coefficients decreasing with speed.

Figure 14 displays the effect of wedge height on the direct translational stiffness and damping coefficients. It can be observed that at the same speed, increasing the wedge height results in a rapid decrease in the direct translational stiffness and a gradual decrease in the direct translational damping coefficient, approaching a certain value. When increasing the speed, the direct translational stiffness coefficient increases, and the direct translational damping coefficient decreases, but the growth rate of both coefficients decreases with speed.

In order to determine the effect of the tilt angle of the thrust collar on the bearing stiffness, we selected ${\phi }_x$ as the variable and studied its effect on the overall dynamic stiffness-damping coefficients of the thrust foil bearing. The results in Figure 15 indicate that as the tilt angle of the thrust collar increases, the direct translational stiffness increases and the direct translational damping coefficient decreases. However, the magnitude of these changes is small, reflecting the adaptive response characteristics of the foil structure, allowing it to deform to improve the composition of the film thickness in response to load changes. Additionally, due to the small thickness of the foils used in this study, they can undergo significant deformation. Therefore, compared to rigid bearings, the thrust foil bearings used in this study show a smaller change in stiffness with tilt angle, indicating a greater stability of the foil bearings.

Additionally, the combined effect of static and dynamic misalignment about the X-axis leads to non-recoverable and non-dissipative forces in the film, manifested as the rapid increase of the cross-coupling stiffness and damping coefficients ($Kz{\phi }_x,K{\phi }_{x}z,Cz{\phi }_x,C{\phi }_{x}z$) from zero. However, they are not symmetrically equal ($Kz{\phi }_x=K{\phi }_{x}z,Cz{\phi }_x=C{\phi }_{x}z$), and their difference increases with the tilt angle. Nevertheless, compared to the direct translational coefficients ($Kzz,Czz$), their values remain small, thus maintaining the stability of the bearing. Furthermore, from the figure, it can be seen that since there is no static deviation around the y-axis (${\phi }_y$ = 0), the cross-coupling stiffness and damping coefficients ($Kz{\phi }_y,K{\phi }_{y}z,Cz{\phi }_y,C{\phi }_{y}z$) are not as sensitive to the static deviation factor (${\phi }_x$). Therefore, the increase in the cross-coupling stiffness and damping coefficients ($Kz{\phi }_y,K{\phi }_{y}z,Cz{\phi }_y,C{\phi }_{y}z$) is due to the increase in bearing pressure caused by the tilt of the thrust collar.

The trend of the rotational direct stiffness and damping coefficients ($K_{{\phi }_x}{\phi }_x,K_{{\phi }_y}{\phi }_y,C_{{\phi }_x}{\phi }_x,C_{{\phi }_y}{\phi }_y$) is similar to that of the direct translational stiffness and damping coefficients ($Kzz,Czz$). As the tilt angle (${\phi }_x$) increases, the rotational direct stiffness coefficient ($K_{{\phi }_x}{\phi }_x,K_{{\phi }_y}{\phi }_y$) increases, while the rotational direct damping coefficient ($C_{{\phi }_x}{\phi }_x,C_{{\phi }_y}{\phi }_y$) decreases. The rotational cross-coupling stiffness and damping coefficients ($K_{{\phi }_x}{\phi }_y,K_{{\phi }_y}{\phi }_x,C_{{\phi }_x}{\phi }_y,C_{{\phi }_y}{\phi }_x$) are not as sensitive to the static deviation factor (${\phi }_x$), and their values are also very small and can be ignored. Therefore, even with a slight tilt of the thrust collar, the rotational stiffness-damping coefficients and cross-coupling stiffness-damping coefficients remain very small, and the direct translational stiffness-damping coefficients play a major role.

Figure 16 shows the effect of the tilt angle of the thrust collar on the direct translational dynamic stiffness and damping coefficients of the thrust foil bearing for different foil thicknesses. It can be seen that the direct translational stiffness coefficient increases with the tilt angle and tends to stabilize with increasing foil thickness. The direct translational damping coefficient, however, is influenced by both the tilt angle and the foil thickness. An increase in foil thickness changes the trend of the direct translational damping coefficient from decreasing with tilt angle to increasing. This is because an increase in foil thickness enhances the stiffness of the foils, gradually reducing the deformation of the foils to nearly zero, increasing the static load of the bearing, and thus increasing the direct translational stiffness and damping coefficients of the bearing.

4. Conclusions

This study simulates the ultra-low temperature environment of a liquid rocket turbopump using liquid nitrogen as the medium, and investigates the performance characteristics of foil thrust bearings in liquid nitrogen lubricant. By utilizing a connecting rod spring model and a finite element model of the foil, the deformation equation of the foil is established. The perturbation method is used to compute the time-varying Reynolds lubrication equation, resulting in one steady pressure control equation and six pressure perturbation control equations. The finite difference method and Newton-Raphson iteration method are employed for solving, analyzing the effects of parameters such as film thickness, rotational speed, and tilt angle on the static and dynamic performance of the thrust bearing. Based on the simulation results, the following conclusions are drawn:

In the analysis of static characteristics, an increase in the friction coefficient will increase the static load and friction torque of the bearing, tending towards a stable value. At high speeds, the influence of the friction coefficient is more significant. As the minimum initial film thickness decreases, both the static load and friction torque of the bearing increase rapidly. With an increase in wedge height, the static load of the bearing first increases and then decreases, while the friction torque decreases slowly. This indicates that there is an optimal wedge height to achieve a greater load capacity while also minimizing friction torque and reducing power loss. When the thrust collar is tilted, the overall static load and friction torque of the foil thrust bearing increase with the increasing tilt angle of the thrust collar and the thickness of the foil. When the foil thickness is small, the provided equivalent stiffness is small, allowing the foil structure to undergo a greater deformation, which improves the composition of the film thickness. As a result, the overall static load and friction torque of the bearing show a smaller change with the increasing tilt angle of the thrust collar, demonstrating the adaptive load response characteristics of the foil thrust bearing and proving its good stability.

In the analysis of dynamic characteristics, as the excitation frequency increases, the direct translational dynamic stiffness of the foil thrust bearing increases while the damping coefficient decreases. Both the direct translational dynamic stiffness and damping coefficient increase with an increase in the friction coefficient and decrease with an increase in the minimum initial film thickness and wedge height. When the thrust collar is tilted, the direct translational stiffness coefficient increases with the tilt angle, and when the foil thickness is small, the direct translational damping coefficient decreases with the tilt angle. However, as the foil thickness increases, the trend of the direct translational damping coefficient changing from a decrease to an increase with the tilt angle reverses, and the increase in foil thickness gradually increases the magnitude of change in the direct translational stiffness and damping coefficients.

In summary, compared to ball bearings, foil thrust bearings exhibit excellent damping characteristics. They can adapt to changes in the minimum film thickness between the rotor and bearing surfaces through foil deformation, greatly improving film thickness and reducing the risk of bearing-rotor contact and abrasion. The adaptive response characteristics of foil thrust bearings, along with their completely non-contact lubrication, allow turbopumps to maintain stable operation even under harsh conditions (such as rotor misalignment, impact loads, etc.). Additionally, foil thrust bearings have high operating speeds, compact structures, require no maintenance, and offer high load-carrying capacity, all of which meet the requirements for repeated recovery. Therefore, based on the findings of this study, the excellent performance characteristics of foil thrust bearings make them a viable alternative to ball bearings for meeting the requirements of repeated recovery in liquid rocket turbopumps.