Numerical analysis of single pad of thrust bearing with the energy equation solved by the characteristic-based split method

Abstract

The solution of the energy equation of thermo-elasto-hydrodynamic analysis for bearings by the finite element method usually leads to convergence difficulties due to the presence of convection terms inherited from the Navier–Stokes equations. In this work, the numerical analysis is performed with finite element method universally by adopting the characteristic-based split method to solve the energy equation. Five case studies of fixed pad thrust bearings have been set up with different geometries, loads, and lubricants. The two-dimensional film pressure is obtained by solving the Reynolds equation with pre-defined axial load on the pad. The energy equation of the lubricant film and the heat transfer equation of the bearing pad are handled by characteristic-based split method and conventional finite element method in three-dimensional space, respectively. Hot oil carry-over effect and variable lubricant viscosity are considered in the simulations. The results of the temperature distributions in the lubricant film and the bearing pad are presented. The possible usability of characteristic-based split method for future thermo-elasto-hydrodynamic analysis is discussed.

Introduction

Thrust bearings are the key factor for the pumps to achieve stable performance and long-term service. Unlike journal bearings, they often operate under conditions with relatively high surface velocity (larger working radius) and substantial axial load. As a consequence, thrust bearings usually operate in the conditions where lubricant temperature rise and thermal-elastic deformation of pad surface is non-ignorable. It is important to get these lubrication properties in the designing stage. This is usually done by the thermo-elasto-hydrodynamic (TEHD) analysis.

TEHD analysis is performed by researchers on thrust bearings with different structures and operating conditions. To the knowledge of the authors, Castelli et al. ¹ did the first TEHD analysis in late 1960s. They took pad deflection, temperature, and pressure into consideration. Fillon and colleagues^2,3 made enormous contribution to the research of the thermal effects of bearings in the last three decades. Monmousseau et al. ⁴ focused on the transient simulation of a tilting-pad journal bearing back to the 1990s. Fillon and colleagues^5,6 continuously made progress in the related fields. The fundamental theories of the TEHD analysis were thoroughly discussed in these works forming the base of this work.

Conventionally, in thrust and journal bearings, the lubricant film between the bearing pad and the rotating runner is modeled by (1) the Reynolds equation, which describes the pressure distribution of the two-dimensional film, and (2) the energy equation, which governs the temperature of the lubricant. The Reynolds equation can be easily discretized and solved by certain numerical methods, such as the finite element method (FEM) ⁷ and the finite difference method (FDM). ⁸ FEM is also used to calculate the deformation of the bearing pad due to both thermal stress and pressure on the pad surface. However, FEM can be applied when the governing equations of the physical system are self-adjoint. ⁹ Unfortunately, it is not the case when dealing with general energy equations either in the analysis of hydrodynamic lubrication or in other fluid dynamics problems. Therefore, when solving energy equations, FDM or similar methods are dominant in the numerical study of lubrication.^10–12 A novel method with the name of characteristic-based split (CBS) method has been developed in the past decades. This makes it possible to use FEM to solve fluid dynamics equations, including the energy equation.

Zienkiewicz and colleagues^13–15 introduced the CBS method in the 1990s. Zienkiewicz et al. ⁹ continued to explore the CBS method in order to solve fluid dynamics problems in different flow regimes including incompressible and compressible flow, turbulent flow, flow through porous media, and so on. According to their theory, if the computational mesh updates its position in an incremental Lagrangian manner along the dominant velocity or “characteristic,” the first-order convective terms in the Navier–Stokes (N-S) equations would disappear. Then, the remaining fluid dynamics problem can be solved by standard Galerkin spatial approximation, which is widely used in FEM approaches. CBS method allows the computational fluid dynamics (CFD) analysis to take advantages of various types of finite elements.

It is notable that some researchers also made efforts to use FEM universally in TEHD analysis. Kim et al. ⁷ utilized FEM techniques to analyze heat transfer effects in the fluid film of a journal bearing. However, the structure of their energy transport equation was not the same with the general energy equation inherited from the N-S equations. Brunetière and Tournerie ¹⁶ performed FEM with upwind schemes on the two-dimensional steady-state energy equation, which was an equation with convective terms as the same as the N-S equations. With the encouragement of the capability of FEM on handling complex geometry, Pahud ¹⁷ did TEHD analysis using streamline-upwind/Petrov–Galerkin (SU/PG) method developed by Brooks and Hughes. ¹⁸ Kucinschi et al. ¹⁹ employed FEM with upwind techniques in a transient TEHD analysis. Recently, Varela and Santos ²⁰ also solved the energy equation by SU/PG formulations. It is reported that the SU/PG method achieved better stabilizing ability than the standard Galerkin procedure in steady-state FEM analysis. This SU/PG method is stemmed from the PG weighting method. So, the SU/PG method achieves the stabilizing effect by modifying the standard weighting/shape functions. Only the steady-state solution can be obtained from SU/PG method.

In this work, the CBS method is applied to solve the energy equation in hydrodynamic lubricated thrust bearing. For calculation of the heat transfer in the bearing pad, the conventional FEM is adopted. The pressure field is given by solving the two-dimensional Reynolds equation, also with FEM. Thus, only FEM, instead of the combination with other methods, is implemented.

Energy equation and CBS method

Energy equation

Assume a coordinate system, in Figure 1, with its origin on the surface of the runner and its z-axis pointing to the bearing pad.

Figure 1.

Sketch of the bearing pad and the coordinate system.

The energy equation in hydrodynamic lubrication theory is in the form of equation (1) in Cartesian coordinates

$$\begin{array}{c}\multicolumn{1}{c}{\rho c(\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}+w\frac{\partial T}{\partial z})=}\\ \multicolumn{1}{c}{k_{\mathrm{f}}\frac{{\partial }^{2}T}{\partial z^2}+\eta ({\left(\frac{\partial u}{\partial z}\right)}^2+{\left(\frac{\partial v}{\partial z}\right)}^2)}\end{array}$$

(1)

where c, k _f, and η are the specific heat, the thermal conductivity, and dynamic viscosity of the lubricant, respectively. Equation (1) describes the temperature distribution of the lubricant film between the runner and bearing pad in three-dimensional space. However, it is common practice to model the pressure as two-dimensional lubricant film, which is governed by the well-known Reynolds equation. There should be some kind of transformation between the two physical fields with different dimensions of numerical treatment. This is done by modeling the velocity field. Equations (2) and (3) show the two velocity components in three-dimensional Cartesian coordinates

$$\begin{array}{c}\multicolumn{1}{c}{u=U_{\mathrm{r}}+(\underset{0}{\overset{z}{\int }}\frac{z\prime \mathrm{d}z\prime }{\eta }-\frac{{\eta }_e}{\eta {\prime }_e}h\underset{0}{\overset{z}{\int }}\frac{\mathrm{d}z\prime }{\eta })\frac{\partial p}{\partial x}-\frac{{\eta }_e{U}_{\mathrm{r}}}{h}\underset{0}{\overset{z}{\int }}\frac{\mathrm{d}z\prime }{\eta }}\\ \multicolumn{1}{c}{+\underset{0}{\overset{z}{\int }}\frac{1}{\eta }\underset{0}{\overset{z\prime }{\int }}{I}_{\mathrm{x}}\mathrm{d}z″\mathrm{d}z\prime -{\eta }_{\mathrm{e}}{\Gamma }_{\mathrm{x}}h\underset{0}{\overset{z}{\int }}\frac{\mathrm{d}z\prime }{\eta }}\end{array}$$

(2)

$$\begin{array}{c}\multicolumn{1}{c}{v=V_{\mathrm{r}}+(\underset{0}{\overset{z}{\int }}\frac{z\prime \mathrm{d}z\prime }{\eta }-\frac{{\eta }_e}{\eta {\prime }_e}h\underset{0}{\overset{z}{\int }}\frac{\mathrm{d}z\prime }{\eta })\frac{\partial p}{\partial y}-\frac{{\eta }_e{U}_{\mathrm{r}}}{h}\underset{0}{\overset{z}{\int }}\frac{\mathrm{d}z\prime }{\eta }}\\ \multicolumn{1}{c}{+\underset{0}{\overset{z}{\int }}\frac{1}{\eta }\underset{0}{\overset{z\prime }{\int }}{I}_{\mathrm{y}}\mathrm{d}z″\mathrm{d}z\prime -{\eta }_{\mathrm{e}}{\Gamma }_{\mathrm{y}}h\underset{0}{\overset{z}{\int }}\frac{\mathrm{d}z\prime }{\eta }}\end{array}$$

(3)

where U _r(x,y,t) and V _r(x,y,t) are the surface velocities of the rotating runner at specified location, p is the film pressure and h(x,y,t) is the film thickness. These formulae are based on those originally proposed by Dowson and Higginson. ²¹ The explanation of other symbols can be found in Appendix 2. With the help of continuity equation, equation (4), the velocity component w along z-direction can be derived

$$\frac{\partial \rho }{\partial t}+\frac{\partial (\rho u)}{\partial x}+\frac{\partial (\rho v)}{\partial y}+\frac{\partial (\rho w)}{\partial z}=0$$

(4)

Substitute equations (2) and (3) and the derived w velocity into equation (1) and follow the procedures developed by Yang and Rodkiewicz, ²² the energy equation can be transformed into the dimensionless form, equation (5). A bar over the symbol indicates a non-dimensional quantity. The relationships between every variable and the reference quantities are listed in Appendix 3. The brief instruction on the derivation of the non-dimensional form of the energy equation is also listed in Appendix 3

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{Pr}·{\mathrm{Re}}^{*}·(\overline{\rho }\frac{\partial \overline{T}}{\partial \overline{t}}+\overline{\rho }\overline{u}\frac{\partial \overline{T}}{\partial \overline{x}}+\overline{\rho }\overline{v}\frac{\partial \overline{T}}{\partial \overline{y}}+\overline{\rho }(-\frac{\overline{q}}{\overline{\rho }\overline{h}})\frac{\partial \overline{T}}{\partial \overline{z}})}\\ \multicolumn{1}{c}{-\frac{1}{{\overline{h}}^2}\frac{{\partial }^2\overline{T}}{\partial {\overline{z}}^2}}\\ \multicolumn{1}{c}{=\mathrm{Pr}·\mathrm{Ec}·\frac{\overline{\eta }}{{\overline{h}}^2}({\left(\frac{\partial \overline{u}}{\partial \overline{z}}\right)}^2+{\left(\frac{\partial \overline{v}}{\partial \overline{z}}\right)}^2)}\end{array}$$

(5)

where $\overline{h}$ is the dimensionless film thickness and $\overline{q}$ is defined as

$$\begin{array}{c}\multicolumn{1}{c}{\overline{q}=\frac{\partial }{\partial \overline{t}}(\overline{h}\underset{0}{\overset{z}{\int }}\overline{\rho }\mathrm{d}\overline{z}\prime )}\\ \multicolumn{1}{c}{+\frac{\partial }{\partial \overline{x}}(\overline{h}\underset{0}{\overset{\overline{z}}{\int }}\overline{\rho }\overline{u}\mathrm{d}\overline{z}\prime )+\frac{\partial }{\partial \overline{y}}(\overline{h}\underset{0}{\overset{\overline{z}}{\int }}\overline{\rho }\overline{v}\mathrm{d}\overline{z}\prime )}\end{array}$$

(6)

In equation (5), Pr, Re*, and Ec are the Prandtl number, the modified Reynolds number, and the Eckert number, respectively.

The CBS method

To illustrate the implementation of CBS method, let us first consider a generalized scalar variable ϕ({x},t), such as the temperature. The transport equation of ϕ is given as equation (7), in the indicial form

$$\frac{\partial \varphi }{\partial t}+\frac{\partial (U_i\varphi )}{\partial x_i}-\frac{\partial }{\partial x_i}\left(k\frac{\partial \varphi }{\partial x_i}\right)+Q=0$$

(7)

where x_i is the component of position tensor {x}, U_i is the component of velocity tensor {U}, and Q is the general source term. Since ϕ ({x},t) is a time-dependent variable, the difference of ϕ at the current time step n and the next time step n+1 can be expressed in a time-stepping form, equation (8)

$$\Delta \varphi ={\varphi }^{n+1}-{\varphi }^n$$

(8)

The time-varying unsteady solution of scalar ϕ can be found if the initial condition and the step length between the two time steps are well defined. A fundamental objective of CBS method is to derive the expression of Δϕ in terms that can, exclusively, be found at the current time step n. Zienkiewicz and colleagues^9,23 explained this method in detail in their works. This is achieved by CBS method in three major procedures. Let x′_i (t) be the coordinate component along the characteristic (velocity) component U _a,i , then ϕ = ϕ({x′},t). Equation (7) can be written as (the diffusion terms are averaged along the characteristic)

$$\frac{\partial \varphi }{\partial t}-\frac{\partial }{\partial x{\prime }_i}\left(k\frac{\partial \varphi }{\partial x{\prime }_i}\right)+Q(\{\mathbf{x}\mathbf{\text{'}}\})=0$$

(9)

During the two successive time steps, the moving coordinate {x′} travels a distance of {δ} = {U _a}Δt, where {U _a} is the averaged value of the characteristic vector and Δt is the time length between the two time steps, respectively. As illustrated in Figure 2, the time variance of ϕ could be written as

$$\begin{array}{c}\multicolumn{1}{c}{\frac{1}{\Delta t}({\varphi }^{n+1}-{{\varphi }^n|}_{(\{\mathbf{x}\}-\{\mathbf{\delta }\})})\approx \gamma {[\frac{\partial }{\partial x_i}\left(k\frac{\partial \varphi }{\partial x_i}\right)-Q]}^{n+1}}\\ \multicolumn{1}{c}{+(1-\gamma ){{[\frac{\partial }{\partial x_i}\left(k\frac{\partial \varphi }{\partial x_i}\right)-Q]}^n|}_{(\{\mathbf{x}\}-\{\mathbf{\delta }\})}}\end{array}$$

(10)

Figure 2.

Time variance of ϕ. ⁹

where γ is a constant. The value ϕⁿ |_({x}−{δ}) can be approximated by the Taylor expansion

$$\begin{array}{c}\multicolumn{1}{c}{{{\varphi }^n|}_{(\{\mathbf{x}\}-\{\mathbf{\delta }\})}\approx {\varphi }^n-{\delta }_i\frac{\partial {\varphi }^n}{\partial x_i}}\\ \multicolumn{1}{c}{+\frac{1}{2}\left(\left({\delta }_i\frac{\partial }{\partial x_i}\right)\left({\delta }_j\frac{\partial }{\partial x_j}\right){\varphi }^n\right)+O(\Delta t^3)}\end{array}$$

(11)

And if we assume γ = 0.5, we have

$$\begin{array}{c}\multicolumn{1}{c}{\frac{1}{2}{\frac{\partial }{\partial x_i}\left(k\frac{\partial \varphi }{\partial x_i}\right)|}_{(\{\mathbf{x}\}-\{\mathbf{\delta }\})}}\\ \multicolumn{1}{c}{\approx \frac{1}{2}\frac{\partial }{\partial x_i}{\left(k\frac{\partial \varphi }{\partial x_i}\right)}^n-\frac{1}{2}\frac{\partial }{\partial x_j}\left(k\frac{\partial }{\partial x_j}\left({\delta }_i\frac{\partial {\varphi }^n}{\partial x_i}\right)\right)+O(\Delta t^2)}\end{array}$$

(12)

$$\frac{1}{2}{Q|}_{(\{\mathbf{x}\}-\{\mathbf{\delta }\})}\approx \frac{Q^n}{2}-\frac{{\delta }_i}{2}\frac{\partial Q^n}{\partial x_i}$$

(13)

Note that distance component δ_i can be considered to be independent on coordinate, and then use equation (14), to approximate this distance. Then, the second term on the right-hand side of equation (12) could be written as equation (15)

$${\delta }_i=U_{a,i}\Delta t$$

(14)

$$-\frac{1}{2}\frac{\partial }{\partial x_j}\left(k\frac{\partial }{\partial x_j}\left({\delta }_i\frac{\partial {\varphi }^n}{\partial x_i}\right)\right)=-\frac{\Delta {\mathrm{tU}}_{a,i}}{2}\frac{\partial }{\partial x_j}\left(k\frac{\partial }{\partial x_j}\left(\frac{\partial {\varphi }^n}{\partial x_i}\right)\right)$$

(15)

In the same sense, the third term on the right-hand side of equation (11) will be

$$\begin{array}{c}\multicolumn{1}{c}{\frac{1}{2}\left(\left({\delta }_i\frac{\partial }{\partial x_i}\right)\left({\delta }_j\frac{\partial }{\partial x_j}\right){\varphi }^n\right)}\\ \multicolumn{1}{c}{=\frac{1}{2}\left({\delta }_i{\delta }_j\frac{\partial }{\partial x_i}\frac{\partial }{\partial x_j}{\varphi }^n\right)}\\ \multicolumn{1}{c}{=\frac{\Delta t^2}{2}\left(U_{a,i}{U}_{a,j}\frac{\partial }{\partial x_i}\frac{\partial }{\partial x_j}{\varphi }^n\right)}\end{array}$$

(16)

If the general diffusion coefficient, k, could be assumed to be constant, equation (15) could be transformed into equation (17)

$$-\frac{\Delta {\mathrm{tU}}_{a,i}}{2}\frac{\partial }{\partial x_j}\left(k\frac{\partial }{\partial x_j}\left(\frac{\partial {\varphi }^n}{\partial x_i}\right)\right)=-\frac{\Delta {\mathrm{tU}}_{a,i}}{2}\frac{\partial }{\partial x_i}\left(\frac{\partial }{\partial x_j}\left(k\frac{\partial {\varphi }^n}{\partial x_j}\right)\right)$$

(17)

Insert equations (11)–(13) into equation (10) and considering equations (14), (16), (17), and γ = 0.5, we would have

$$\begin{array}{c}\multicolumn{1}{c}{\frac{1}{\Delta t}({\varphi }^{n+1}-{\varphi }^n)\approx }\\ \multicolumn{1}{c}{-\{U_{a,i}\frac{\partial {\varphi }^n}{\partial x_i}-\frac{1}{2}(\frac{\partial }{\partial x_i}{\left(k\frac{\partial \varphi }{\partial x_i}\right)}^{n+1}+\frac{\partial }{\partial x_i}{\left(k\frac{\partial \varphi }{\partial x_i}\right)}^n)}\\ \multicolumn{1}{c}{+\frac{1}{2}(Q^{n+1}+Q^n)\}}\\ \multicolumn{1}{c}{+\frac{\Delta t}{2}\{U_{a,i}{U}_{a,j}\frac{\partial }{\partial x_i}\frac{\partial }{\partial x_j}{\varphi }^n-U_{a,i}\frac{\partial }{\partial x_i}\left(\frac{\partial }{\partial x_j}{\left(k\frac{\partial \varphi }{\partial x_j}\right)}^n\right)+U_{a,i}\frac{\partial Q^n}{\partial x_i}\}}\end{array}$$

(18)

Let us focus on the following term

$$\frac{\Delta t}{2}{U}_{a,i}{U}_{a,j}\frac{\partial }{\partial x_i}\frac{\partial }{\partial x_j}{\varphi }^n$$

(19)

The averaged characteristic component U _a,i could be written as

$$U_{a,i}=\frac{1}{2}(U_i^{n+1}+{U_i^n|}_{(\{\mathbf{x}\}-\{\mathbf{\delta }\})})$$

(20)

Again, employ the Taylor expansion ⁹

$${U_i^n|}_{(\{\mathbf{x}\}-\{\mathbf{\delta }\})}\approx U_i^n-\Delta {\mathrm{tU}}_j^n\frac{\partial U_i^n}{\partial x_j}+O(\Delta t^2)$$

(21)

Insert equations (20) and (21) into equation (19) and neglecting the higher order terms with respect to Δt, we can write

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\Delta t}{2}{U}_{a,i}{U}_{a,j}\frac{\partial }{\partial x_i}\frac{\partial }{\partial x_j}{\varphi }^n}\\ \multicolumn{1}{c}{\approx \frac{\Delta t}{2}\left(\frac{1}{2}(U_i^{n+1}+U_i^n)\right)\left(\frac{1}{2}(U_j^{n+1}+U_j^n)\right)\frac{\partial }{\partial x_i}\frac{\partial }{\partial x_j}{\varphi }^n}\end{array}$$

(22)

Furthermore, the following relationship can be assumed

$$\begin{array}{c}\multicolumn{1}{c}{\frac{1}{2}(\frac{\partial }{\partial x_i}{\left(k\frac{\partial \varphi }{\partial x_i}\right)}^{n+1}+\frac{\partial }{\partial x_i}{\left(k\frac{\partial \varphi }{\partial x_i}\right)}^n)}\\ \multicolumn{1}{c}{\approx \frac{\partial }{\partial x_i}{\left(k\frac{\partial \varphi }{\partial x_i}\right)}^n+O(\Delta t)}\end{array}$$

(23)

$$\frac{1}{2}(Q^{n+1}+Q^n)\approx Q^n+O(\Delta t)$$

(24)

$$\frac{1}{2}(U_i^{n+1}+U_i^n)\approx U_i^n+O(\Delta t)$$

(25)

From equations (20) to (25), and refer equation (14), the formulation for Δϕ is expressed based on equation (18) as

$$\begin{array}{c}\multicolumn{1}{c}{\Delta \varphi ={\varphi }^{n+1}-{\varphi }^n}\\ \multicolumn{1}{c}{\approx -\Delta t\{U_i^n\frac{\partial {\varphi }^n}{\partial x_i}-\frac{\Delta t}{2}{U}_j^n\frac{\partial U_i^n}{\partial x_j}\frac{\partial {\varphi }^n}{\partial x_i}-\frac{\partial }{\partial x_i}{\left(k\frac{\partial \varphi }{\partial x_i}\right)}^n+Q^n\}}\\ \multicolumn{1}{c}{+\frac{\Delta t^2}{2}\{U_i^n{U}_j^n\frac{\partial }{\partial x_i}\frac{\partial }{\partial x_j}{\varphi }^n-U_i^n\frac{\partial }{\partial x_i}\left(\frac{\partial }{\partial x_j}{\left(k\frac{\partial \varphi }{\partial x_j}\right)}^n\right)+U_i^n\frac{\partial Q^n}{\partial x_i}\}}\end{array}$$

(26)

In fact, because the above equation is expressed in tensor notation, the following simple relationship holds

$$U_j^n\frac{\partial U_i^n}{\partial x_j}\frac{\partial {\varphi }^n}{\partial x_i}=U_i^n\frac{\partial U_j^n}{\partial x_i}\frac{\partial {\varphi }^n}{\partial x_j}$$

(27)

Then, we can combine the following two terms

$$\begin{array}{c}\multicolumn{1}{c}{U_i^n\frac{\partial U_j^n}{\partial x_i}\frac{\partial {\varphi }^n}{\partial x_j}+U_i^n{U}_j^n\frac{\partial }{\partial x_i}\frac{\partial }{\partial x_j}{\varphi }^n}\\ \multicolumn{1}{c}{=U_i^n(\frac{\partial U_j^n}{\partial x_i}\frac{\partial {\varphi }^n}{\partial x_j}+U_j^n\frac{\partial }{\partial x_i}\frac{\partial }{\partial x_j}{\varphi }^n)}\\ \multicolumn{1}{c}{=U_i^n\frac{\partial }{\partial x_i}\left(U_j^n\frac{\partial }{\partial x_j}{\varphi }^n\right)}\end{array}$$

(28)

Finally, considering equation (28), equation (26) becomes

$$\begin{array}{c}\multicolumn{1}{c}{\Delta \varphi ={\varphi }^{n+1}-{\varphi }^n}\\ \multicolumn{1}{c}{\approx -\Delta t\{U_i^n\frac{\partial {\varphi }^n}{\partial x_i}-\frac{\partial }{\partial x_i}{\left(k\frac{\partial \varphi }{\partial x_i}\right)}^n+Q^n\}}\\ \multicolumn{1}{c}{+\frac{\Delta t^2}{2}\{U_i^n\frac{\partial }{\partial x_i}\left(U_j^n\frac{\partial {\varphi }^n}{\partial x_j}\right)-U_i^n\frac{\partial }{\partial x_i}\left(\frac{\partial }{\partial x_j}{\left(k\frac{\partial \varphi }{\partial x_j}\right)}^n\right)+U_i^n\frac{\partial Q^n}{\partial x_i}\}}\end{array}$$

(29)

Generally, equation (29) is the same with that proposed by Zienkiewicz et al. ⁹ At this point, ϕ^{n +1} can be approximated by solution at time step n. As in any other FEM implementations, the standard Galerkin spatial approximation, equation (30), can now be used to substitute the scalar unknown in equation (29)

$$\varphi =[\mathbf{N}]\{\stackrel{~}{}\varphi \}$$

(30)

where [N] is the shape function row vector and $\{\stackrel{~}{}\varphi \}$ is the nodal value vector. After integration by parts and taking the standard shape functions [N]^T as the weighting functions, the time-stepping form of the integrated Δϕ is

$$\begin{array}{c}\multicolumn{1}{c}{[\mathbf{M}]({\{\stackrel{~}{}\varphi \}}^{n+1}-{\{\stackrel{~}{}\varphi \}}^n)=}\\ \multicolumn{1}{c}{-\Delta t[([\mathbf{C}]{\{\stackrel{~}{}\varphi \}}^n+[\mathbf{K}]{\{\stackrel{~}{}\varphi \}}^n+{\{\mathbf{f}\}}^n)}\\ \multicolumn{1}{c}{-\Delta t([{\mathbf{K}}_{\mathrm{u}}]{\{\stackrel{~}{}\varphi \}}^n+{\{{\mathbf{f}}_{\mathrm{s}}\}}^n)]}\end{array}$$

(31)

where [M], [C], [K], and [K _u] are the matrices while {f} and {f _s} are the column vectors. Equations (32)–(37) are the expressions of these terms

$$[\mathbf{M}]=\underset{\Omega }{\int }{[\mathbf{N}]}^{\mathrm{T}}[\mathbf{N}]\mathrm{d}\Omega$$

(32)

$$[\mathbf{C}]=\underset{\Omega }{\int }{[\mathbf{N}]}^{\mathrm{T}}\frac{\partial }{\partial x_i}(U_i[\mathbf{N}])\mathrm{d}\Omega$$

(33)

$$[\mathbf{K}]=\underset{\Omega }{\int }\frac{\partial {[\mathbf{N}]}^{\mathrm{T}}}{\partial x_i}k\frac{\partial [\mathbf{N}]}{\partial x_i}\mathrm{d}\Omega$$

(34)

$$\left[{\mathbf{K}}_{\mathrm{u}}\right]=-\frac{1}{2}\underset{\Omega }{\int }\frac{\partial }{\partial x_i}\left(U_i{[\mathbf{N}]}^{\mathrm{T}}\right)\frac{\partial }{\partial x_j}(U_j[\mathbf{N}])\mathrm{d}\Omega$$

(35)

$$\{\mathbf{f}\}=\underset{\Omega }{\int }{[\mathbf{N}]}^{\mathrm{T}}Q\mathrm{d}\Omega +\mathrm{b}.\mathrm{t}.$$

(36)

$$[{\mathbf{f}}_{\mathrm{s}}]=-\frac{1}{2}\underset{\Omega }{\int }\frac{\partial }{\partial x_i}\left(U_i{[\mathbf{N}]}^{\mathrm{T}}\right)Q\mathrm{d}\Omega +\mathrm{b}.\mathrm{t}.$$

(37)

The term “b.t.” in equations (36) and (37) are the results of integration by parts and they stand for boundary conditions (Neumann boundaries). After choosing the time step length, Δt, and the initial value of $\{\stackrel{~}{}\varphi \}$ , the solution at every time step can be found from equation (8) by solving equation (31).

The time step length Δt has strong influence on the stability of CBS method. The largest Δt for each element is evaluated with respect to the local fluid velocity, the thermal conductivity, and the average element length (h _avg), along the element’s characteristic. Zienkiewicz et al. ⁹ discussed several methods to determine Δt. For this work, Δt is obtained by the following procedures based on the method proposed by Zienkiewicz et al. ⁹ Let Pe be the Peclet number, and C _c be the critical Courant number of an element

$$\mathrm{Pe}=\frac{{\mathrm{Uh}}_{\mathrm{avg}}}{2k}$$

(38)

$$C_{\mathrm{c}}={(\frac{1}{{\mathrm{Pe}}^2}+1)}^{0.5}-\frac{1}{\mathrm{Pe}}$$

(39)

The Δt of the element is chosen to satisfy

$$\frac{U\Delta t}{h_{\mathrm{avg}}}<C_{\mathrm{c}}$$

(40)

The global Δt is the smallest one among all the elements.

Remark on the energy equation

Let us rearrange the terms in equation (5) to transform it into the similar structure of equation (7). Then, we get equation (41)

$$\frac{\partial \overline{T}}{\partial \overline{t}}+(\overline{u}\frac{\partial \overline{T}}{\partial \overline{x}}+\overline{v}\frac{\partial \overline{T}}{\partial \overline{y}}+{\overline{w}}_{\mathrm{E}}\frac{\partial \overline{T}}{\partial \overline{z}})-{\overline{k}}_{\mathrm{E}}\frac{{\partial }^2\overline{T}}{\partial {\overline{z}}^2}+{\overline{Q}}_{\mathrm{E}}=0$$

(41)

where

$${\overline{w}}_{\mathrm{E}}=-\frac{\overline{q}}{\overline{\rho }\overline{h}}$$

(42)

$${\overline{k}}_{\mathrm{E}}=\frac{1}{\mathrm{Pr}·{\mathrm{Re}}^{*}·\overline{\rho }{\overline{h}}^2}$$

(43)

$${\overline{Q}}_{\mathrm{E}}=-\frac{\mathrm{Ec}·\overline{\eta }}{{\mathrm{Re}}^{*}·\overline{\rho }{\overline{h}}^2}({\left(\frac{\partial \overline{u}}{\partial \overline{z}}\right)}^2+{\left(\frac{\partial \overline{v}}{\partial \overline{z}}\right)}^2)$$

(44)

Compared to equation (7), only the z-direction diffusion term remains. Thus, equation (34) can be written as

$$[\mathbf{K}]=\underset{\Omega }{\int }\frac{\partial {[\mathbf{N}]}^{\mathrm{T}}}{\partial z}{\overline{k}}_{\mathrm{E}}\frac{\partial [\mathbf{N}]}{\partial z}\mathrm{d}\Omega$$

(45)

The temperature field now can be solved by a time-stepping manner after the selection of initial value and boundary conditions. The actual selected type of finite element and the technique of obtaining the explicit form of equation (31) by mass-lumping will be discussed in later section.

Reynolds equation and heat transfer equation

Reynolds equations

To illustrate the usability of CBS method in TEHD analysis of thrust bearing, we select the Reynolds equation utilized by Yang and Rodkiewicz. ²² Equations (46) and (47) are the Reynolds equation and its non-dimensional equivalent, respectively. The non-dimensional terms in these equations are listed in Appendices 2 and 3

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\partial }{\partial x}\left({\left(\frac{\rho }{\eta }\right)}_{\mathrm{e}}{h}^3\frac{\partial p}{\partial x}\right)+\frac{\partial }{\partial y}\left({\left(\frac{\rho }{\eta }\right)}_{\mathrm{e}}{h}^3\frac{\partial p}{\partial y}\right)=}\\ \multicolumn{1}{c}{6U\frac{\partial }{\partial x}({\rho }^{*}h)+6V\frac{\partial }{\partial y}({\rho }^{*}h)+12\frac{\partial }{\partial t}({\rho }_{\mathrm{e}}h)+12\rho \frac{\partial h}{\partial t}}\\ \multicolumn{1}{c}{+12\frac{\partial }{\partial x}({\Omega }_{\mathrm{x}}{h}^3)+12\frac{\partial }{\partial y}({\Omega }_y{h}^3)}\end{array}$$

(46)

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\partial }{\partial \overline{x}}\left({\overline{\left(\frac{\rho }{\eta }\right)}}_{\mathrm{e}}{\overline{h}}^3\frac{\partial \overline{p}}{\partial \overline{x}}\right)+\frac{\partial }{\partial \overline{y}}\left({\overline{\left(\frac{\rho }{\eta }\right)}}_{\mathrm{e}}{\overline{h}}^3\frac{\partial \overline{p}}{\partial \overline{y}}\right)=}\\ \multicolumn{1}{c}{6\frac{\partial }{\partial \overline{x}}({\overline{\rho }}^{*}\overline{h})+6\frac{\partial }{\partial \overline{y}}({\overline{\rho }}^{*}\overline{h})+12\frac{\partial }{\partial \overline{t}}({\overline{\rho }}_{\mathrm{e}}\overline{h})+12\overline{\rho }\frac{\partial \overline{h}}{\partial \overline{t}}}\\ \multicolumn{1}{c}{+12{\mathrm{Re}}^{*}\frac{\partial }{\partial \overline{x}}({\overline{\Omega }}_{\mathrm{x}}{\overline{h}}^3)+12{\mathrm{Re}}^{*}\frac{\partial }{\partial \overline{y}}({\overline{\Omega }}_y{\overline{h}}^3)}\end{array}$$

(47)

where h(x,y,t) is the film thickness, which could be expressed as equation (48) when taking the thermal-elastic deformation W into consideration

$$h(x,y,t)=h_0+r(x,y)·\sin \alpha ·\sin (\frac{{\theta }_0}{2}-\theta )+W(x,y,t)$$

(48)

In the above expression, r(x,y) is the radius at specific location and h ₀ is the film thickness at the pad center. α is one of the geometric parameters in Table 1 and Figure 1. The term W in equation (48) contains both the elastic deformation under the pressure of the lubricant film and the thermal deformation due to temperature rise in the bearing pad. It is noted here that the term W will not be considered in the demonstration computation performed later in this work.

Table 1.

Thrust bearing geometry and properties.

Description	Value
	Cases 1–4	Case 2
Geometry
Outer radius, r 2 (m)	2.6	0.45
Inner radius, r 1 (m)	1.625	0.16
Pad angle, θ 0 (°)	18.5	22.5
Pad thickness, b (m)	0.23	0.1
Pad inclination, α (rad)	0.0002	0.0005
Operating conditions
Rotational speed (r/min)	92	1780
Ambient lubricant temperature (K)	313.15	333.15
Axial load (N)	See Table 2
Lubricant properties
Reference temperature, T 0 (K)	313.15	333.15
Reference density, ρ 0 (kg m−3)	857	1000
Reference viscosity, η 0 (Pa s)	0.0047	0.00047
Thermal conductivity, k f (W m−1 K−1)	0.13	0.6617
Specific heat, c (J kg−1 K−1)	1783	4152
Pad material properties
Thermal conductivity, k p (W m−1 K−1)	50	7.96
Thermal transfer coefficient, h p (W m−2 K−1)	750	2000

Heat transfer equation

The temperature distribution inside the bearing pad is governed by the heat transfer equation. If we assume no source term is in this equation, it reduces into the Poisson’s equation in three-dimensional Cartesian coordinates

$$k_{\mathrm{p}}(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z_{\mathrm{p}}^2})={\rho }_{\mathrm{p}}{c}_{\mathrm{p}}\frac{\partial T}{\partial t}$$

(49)

where k _p, ρ _p, and c _p are the thermal conductivity, density, and specific heat of the bearing pad, respectively. z _p is the z-coordinate of the bearing pad with the origin on the film–pad surface. The non-dimensional form of equation (49) is equation (50). The non-dimensional reference values and procedures are listed in Appendix 3

$${\mathrm{Fo}}_{\mathrm{p}}(\frac{{\partial }^2\overline{T}}{\partial {\overline{x}}^2}+\frac{{\partial }^2\overline{T}}{\partial {\overline{y}}^2}+\frac{{\partial }^2\overline{T}}{\partial {\overline{z}}_{\mathrm{p}}^2})=\frac{\partial \overline{T}}{\partial \overline{t}}$$

(50)

where

$${\mathrm{Fo}}_{\mathrm{p}}=\frac{k_{\mathrm{p}}}{{\rho }_{\mathrm{p}}{c}_{\mathrm{p}}{B}_{\mathrm{m}}{U}_{\mathrm{m}}}$$

(51)

where B _m is the reference length and U _m is the nominal average velocity (Appendix 3).

Boundary conditions

Classical Reynolds boundary conditions are used when solving the Reynolds equation. It is assumed that the bearing is fully flooded and the pressure is equal to 0 at the boundaries.

As illustrated in Figure 3, the heat flux continuity is applied at the film–pad interface. On the other surfaces of the pad, except the leading edge, a universal heat transfer coefficient h _p is used. The heat exchange at the film/runner interface is set to be 0. This boundary condition will result in overestimation of the lubricant temperature. These boundary conditions could be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{{k\frac{\partial T}{\partial z}|}_{z=h}={k_{\mathrm{p}}\frac{\partial T}{\partial z_{\mathrm{p}}}|}_{z_{\mathrm{p}}=0}}\\ \multicolumn{1}{c}{-k_{\mathrm{p}}{\left(\frac{\partial T}{\partial {\mathbf{n}}_{\mathrm{vec}}}\right)}_{\mathrm{s}}=h_{\mathrm{p}}(T_{\mathrm{s}}-T_0)}\end{array}$$

(52)

Figure 3.

Surfaces and edges of the computational domain.

where T ₀ is the ambient lubricant temperature. T _s is the surface temperature at faces of the bearing pad other than the film–pad interface. n _vec is the surface normal vector. The non-dimensional form of equation (52) is

$${\frac{\overline{k}}{\overline{h}}\frac{\partial \overline{T}}{\partial \overline{z}}|}_{\overline{z}=1}={\frac{\partial \overline{T}}{\partial {\overline{z}}_{\mathrm{p}}}|}_{{\overline{z}}_{\mathrm{p}}=0}$$

(53)

$$-{\left(\frac{\partial \overline{T}}{\partial {\mathbf{n}}_{\mathrm{vec}}}\right)}_{\mathrm{s}}=\mathrm{Nu}({\overline{T}}_{\mathrm{s}}-1)$$

(54)

where Nu is the Nusselt number defined as

$$\mathrm{Nu}=\frac{h_{\mathrm{p}}{B}_{\mathrm{m}}}{k_{\mathrm{p}}}$$

(55)

The mixing effect between the hot lubricant from the exit of the preceding bearing pad and the relatively cold lubricant supply will rise the actual temperature at the leading edge of the film. This mixing effect is considered by introducing equation (56) ²⁴ in the single-pad simulation

$$T_{\mathrm{in}}=T_{sup}+\Delta T_{\mathrm{avg}}\left(\frac{k_{\mathrm{mix}}}{2-2k_{\mathrm{mix}}}\right)$$

(56)

where k _mix is the hot oil carry-over coefficient. In equation (56), T _in, T _sup, and ΔT _avg are the temperature at the leading edge, the temperature of the supplied lubricant, and the circumferential temperature rise across the bearing pad, respectively. The values of T _in and ΔT _avg are evaluated by averaging the temperature solutions on the leading edge and tailing edge of the lubricant film. The mixing effect takes place at very film–pad temperature iteration process illustrated in Figure 4 (next section).

Figure 4.

Solution flow chart. n ₁, n ₂, n ₃, and n ₄ are the indicators of the iteration step for different physical fields.

Some considerations should be taken before proceeding to the simulation stage. When dealing with thrust bearing with sector-shaped pads, Reynolds equation (46) is used to be defined in polar coordinates to grasp the geometry more accurately. However, in this work, we chose the equations in Cartesian coordinates to get a set of universally applicable formulations for finite element analysis. The isoparametric finite elements ²⁵ are applied to deal with the curvature of the geometry of the bearing pad.

Although the Reynolds equation is two-dimensional, the energy equation and the temperature distribution in the bearing pad are solved in three-dimensional space. The pressure field could be considered as the so-called middle-layer ²² of the lubricant film while the temperature distribution is across the apace between the bearing pad and the rotating runner.

Single-pad simulation

Two types of finite elements are utilized for different governing equations. The eight-node tri-linear isoparametric cubic element ²⁶ is used for discretization of the film temperature and the temperature in the bearing pad. For Reynolds equation, the four-node bi-linear isoparametric rectangular element is adopted. The shape function matrices, which are associated with those two finite element types, can be found in the work of Zienkiewicz et al. ²⁶ The integration listed from equations (32) to (37) are calculated by Gaussian integration. ²⁵

Five test cases are set up to evaluate the ability of the CBS method to solve the energy equation. As listed in Table 1, the bearing geometries (except the pad inclination) and properties of Cases 1–4 are the same with the large Itaipu thrust bearing studied by Wodtke et al. ²⁷ The lubricant is oil with its dynamic viscosity at T ₀ = 313.15 K (it is also the reference and oil supply temperature in the simulation) equals 0.047 Pa s. The dynamic viscosity of the oil at higher temperature is expressed by equation (57) ^28,29 with β = −0.03297. Case 5 takes the geometry from the water-lubricated thrust bearing in a pump. The viscosity versus temperature property of water is described by the expression of equation (58). ³⁰ In equation (58), η _293.15 is the dynamic viscosity of water at 293.15 K and the exponent A(T) is expressed by equation (59). The heat transfer coefficient, h _p, for Case 5 is assumed to be 2000 W m⁻² K⁻¹. Case 5 is set up to test the capability of the CBS method in simulating lubrication condition with relatively lower viscosity and thinner lubricant film. The potential effects of turbulence due to the low viscosity of water and the large radius of pad geometry are simply ignored. It is emphasized that this work is made to demonstrate the usability of CBS method in solving the energy equation in thrust bearings

$$\eta (T)={\eta }_0\exp (\beta (T-T_0))$$

(57)

$$\eta (T)={\eta }_{293.15}\times {10}^{A(T)}$$

(58)

$$\begin{array}{c}\hfill A(T)=\frac{293.15-T}{T-177.15}\times (1.2378-1.303\times {10}^{-3}(293.15-T)\hfill \\ \multicolumn{1}{c}{+3.06\times {10}^{-6}{(293.15-T)}^2+2.55\times {10}^{-8}{(293.15-T)}^3)}\end{array}$$

(59)

Two kinds of meshes are generated for the analysis. One is two-dimensional mesh for the Reynolds equation. The second is three-dimensional mesh for the energy equation. Therefore, the lubricant film in the analysis of temperature is modeled as layered structure, with several nodes and elements in the z-direction. Three-dimensional mesh is also used for the heat transfer computation. For Cases 1–4, the node number along every coordinate axis is set according to the referenced work. ²⁷ An additional finer mesh is deployed to test the dependence of the simulation on the mesh. For Case 5, the similar mesh structure is adapted. The mesh information and the axial load settings could be found in Table 2. The axial loads of Cases 1 and 2 have the same value described in the work of Wodtke et al. ²⁷ The hot oil carry-over coefficients, k _mix, for oil-lubricated thrust bearings (Cases 1–4) are assumed to be 0.589 based on the work done by Ettles. ²⁴ For Case 5, k _mix is assumed to be 0.6.

Table 2.

Test case settings.

Case	Lubricant	Number of elements			Axial load (N)	Hot oil carry-over coefficient, k mix
		Radius direction	Circumferential direction	Axial direction
1	Oil	30	48	20	1.731 × 106	0.589
2	Oil	45	72	30	1.731 × 106	0.589
3	Oil	30	48	20	2.5965 × 106	0.589
4	Oil	30	48	20	0.8655 × 106	0.589
5	Water	50	50	20	1.5 × 104	0.6

The general form of CBS method, equation (31), does not exhibit good efficiency since the solution of the full-scale linear equations should be found at every time step. And the numerical simulation maybe unstable due to some improper initial values. It is suggested that we can transform equation (31) into its explicit form, equation (60), by the mass-lumping method. Then, no actual decomposition of matrix is needed in the later iteration process

$$\begin{array}{c}\multicolumn{1}{c}{{\{\stackrel{~}{}\varphi \}}^{n+1}-\{\stackrel{~}{}\varphi {\}}^n=}\\ \multicolumn{1}{c}{-\Delta t{[{\mathbf{M}}_{\mathrm{l}}]}^{-1}[([\mathbf{C}]{\{\stackrel{~}{}\varphi \}}^n+[\mathbf{K}]{\{\stackrel{~}{}\varphi \}}^n+{\{\mathbf{f}\}}^n)}\\ \multicolumn{1}{c}{-\Delta t([{\mathbf{K}}_{\mathrm{u}}]{\{\stackrel{~}{}\varphi \}}^n+{\{{\mathbf{f}}_{\mathrm{s}}\}}^n)]}\end{array}$$

(60)

where [M _l] is the lumped mass matrix which is a diagonal matrix defined as

$$M_{\mathrm{l},\mathrm{ab}}={\delta }_{\mathrm{ab}}\underset{\Omega }{\int }{N}_{\mathrm{a}}\mathrm{d}\Omega (\mathrm{no}\mathrm{sum})$$

(61)

in which M _l,ab is the element at Row a and Column b in matrix [M _l], δ _ab is the Kronecker delta, and N _a is the shape function of node “a” in an element. Zienkiewicz et al. ⁹ discussed the validity of this method in steady-state analysis.

The Reynolds equation and the heat transfer equation are solved for their steady-state solutions by conventional FEM.

Several procedures are performed to solve the equations, as illustrated in Figure 4. There are two sub-iterations in the numerical solution procedures. One is the pressure–load iteration, and the other is the film–pad temperature iteration. The pressure distribution is first obtained by the initial guess of film thickness (h ₀ = h ₀₀ equals 1 × 10⁻⁴ m for Cases 1–4 and equals 6 × 10⁻⁵ for Case 5). Then, the integral of pressure over the area of the pad surface is compared with the axial load to get a correction of film thickness. An iterative process with film thickness correction continues until the axial load is balanced by the pressure integration over the pad surface. The above pressure–load iteration continues until the norm of the pressure difference of two successive iterations is smaller than the limit, e _p. Then, the velocity field is calculated by equations (2) and (3), in their dimensionless forms. In the film–pad temperature iteration, the temperature distribution is obtained through solving equation (41) by the CBS method with mass-lumping, equation (60). Then, the heat transfer equation of the pad is solved with the heat flux on the film–pad interface as the boundary condition. After the heat transfer solution, the temperature of the film–pad interface will be updated and the hot oil carry-over effect will take place. Then, the CBS method will be launched again to get the film temperature and the heat flux on the film–pad interface. The convergence of film–pad temperature iteration is achieved when the normalized difference of the successive temperature on the film–pad interface, $e_{\mathrm{int}}^{n3}$ , is smaller than the pre-defined limit, e _int. This is expressed by equation (62)

$$\frac{|{\{{\mathbf{T}}_{\mathrm{int}}\}}^{n_3+1}-{\{{\mathbf{T}}_{\mathrm{int}}\}}^{n_3}|}{\left|{\{{\mathbf{T}}_{\mathrm{int}}\}}^{n_3+1}\right|}=e_{\mathrm{int}}^{n_3}<e_{\mathrm{int}}$$

(62)

where {T _int} is the temperature of the film–pad interface and superscript n ₃ indicates the previous film–pad temperature iteration step. For the CBS method, the convergence is monitored through the normalized value expressed by equation (63)

$$\frac{|{\{\Delta \mathbf{T}\}}^{n_2+1}-{\{\Delta \mathbf{T}\}}^{n_2}|}{\left|{\{\Delta \mathbf{T}\}}^{n_2+1}\right|}=e_{\mathrm{cbs}}^{n_2}<e_{\mathrm{cbs}}$$

(63)

where $e_{\mathrm{cbs}}^{n_2}$ is the monitored convergence control value of the current time step and e _cbs is the limit. An under relaxation coefficient of 0.01 is used for the implementation of the heat flux boundary condition, equation (54), when solving the heat transfer equation, equation (50). After a complete film–pad temperature iteration, the normalized differences between the current and the previous solutions of film temperature ( $e_{\mathrm{T}}^{n_4}$ ) and the viscosity of the lubricant ( $e_{\eta }^{n_4}$ ) are botch checked for convergence by equations (64) and (65)

$$\frac{|{\{\mathbf{T}\}}^{n_4+1}-{\{\mathbf{T}\}}^{n_4}|}{\left|{\{\mathbf{T}\}}^{n_4+1}\right|}=e_{\mathrm{T}}^{n_4}<e_{\mathrm{T}}$$

(64)

$$\frac{|{\{\mathbf{\eta }\}}^{n_4+1}-{\{\mathbf{\eta }\}}^{n_4}|}{\left|{\{\mathbf{\eta }\}}^{n_4+1}\right|}=e_{\eta }^{n_4}<e_{\eta }$$

(65)

If the convergence limitations are not met, the viscosity of the lubricant will be modified by equation (57) (equation (58)) according to the local film temperature of the current solution and the whole solution procedures will be performed all over again from solving the Reynolds equation. All the convergence control parameters are listed in Table 3.

Table 3.

Convergence control parameters.

e p	e int	e cbs	e T	e η
10−4	5 × 10−3	10−5	10−4	10−2

Results and discussion

Figures 5 and 6 are the temperature distributions on the film–pad interface of Case 1 obtained from (1) the lubricant film solved by the CBS method and (2) the heat transfer in the bearing pad. These two figures show that the temperature in the lubricant film rises from the leading edge to the tailing edge. The maximum temperature appears at the location where the film thickness is relatively small and the velocity gradient along the z-direction is intense. Figure 7 shows the isothermals at the film–pad interface along the centerline in the direction of the bearing angle. The isothermal of Case 2, which uses the refined mesh resolution, is shown in Figure 8. Little deviation could be found between those two figures, indicating that the original coarser mesh can be used to obtain a mesh-independent solution of the film temperature. The coarser mesh was also used in the analysis done by Wodtke et al., ²⁷ in which the simulation was performed with FDM. This means that the CBS method could obtain steady-state solution of the energy equation with the same mesh density. It is worth to be noted that the vertical axes and the graph in both Figures 7 and 8 are separated by the dashed lines, which indicate the film–pad interface. The lower portions of these axes are non-dimensionalized by referencing the initial film thickness, h ₀₀. The upper portions take B _m as the reference value (see Appendix 3 for the expression of B _m). In order to show the temperature in a more clear way, in all these figures the unit of the temperature is chosen as °C.

Figure 5.

Film temperature profile on film–pad interface: Case 1. The unit of temperature is °C.

Figure 6.

Pad temperature profile on film–pad interface: Case 1. The unit of temperature is °C.

Figure 7.

Isothermal between the film and the pad: Case 1. The unit of temperature is °C.

Figure 8.

Isothermal between the film and the pad: Case 2. The unit of temperature is °C.

Figure 9 shows the temperature curves along the center of the bearing pad in the circumferential direction on the film–pad interface of Cases 1, 3, and 4. The temperature trends in all these figures are roughly the same. As expected, the more the axial load, the higher the peak temperature. Significant deviation exists when the result of Case 1 is compared with the analysis conducted by Wodtke et al. ²⁷ This may be attributed to the fact that instead of complete TEHD simulation of the tilting-pad thrust bearing, a thermohydrodynamic (THD) analysis of a fixed pad thrust bearing with the same geometry and lubricant property is performed. In this way, the CBS method could be implemented to solve the energy equation without calculating the tilting angle of the bearing pad.

Figure 9.

Temperature along the centerline in circumferential direction: Cases 1, 3, and 4. The unit of temperature is °C.

The film temperature solution of Case 5 is illustrated in Figure 10 in the form of isothermal diagram. The lubricant of Case 5 is liquid water, which has much lower viscosity in the temperature range of this study. Water-lubricated thrust bearing usually operates under thinner lubricant film. In this work, the oil film thickness at the center of the bearing pad of Case 1 is about two times the thickness of the water film of Case 5. Case 5 could be a demonstrating simulation of THD analysis with low viscosity lubricant. The result in Figure 10 indicates that the temperature rise across the bearing pad along the centerline in the circumferential direction is about 2.4°C.

Figure 10.

Isothermal between the film and the pad: Case 5. The unit of temperature is °C.

It is observed that the leading edge temperature rises substantially due to the hot oil carry-over effect in Cases 1–4. However, the temperature rise at the leading edge of the water-lubricated thrust bearing under the current simulation settings of this study is negligible.

The typical convergence histories of CBS method, equation (60), during the iteration process of Cases 1 and 5 are illustrated in Figure 11. The average time steps needed to achieve the convergence limitation are less than 3000. As shown in Figure 11, when the simulation approximately reaches the steady state, the value expressed by equation (63) is small, making the solution of {T} slightly oscillate about its target value. From Figure 11, it is also observed that the iteration processes tend to converge more rapidly when the steady-state solution is about to be reached.

Figure 11.

Convergence history of $e_{\mathrm{cbs}}^{n_2}$ .

Although the CBS method and SU/PG method both can be considered as being within the scope of FEM, significant difference can be found between them. SU/PG method mainly deals with the oscillatory solution of convection–diffusion problem by using the PG type of weighting in the FEM formulations. The self-adjointness is achieved by choosing the optimal parameter of the shape functions according to the element’s Peclet number. It is essentially a kind of guesswork. One can easily find out that the mass matrices resulting from SU/PG method are non-symmetric, ⁹ which will cause difficulty in subsequent solution process. The CBS method, on the other side, obtains its self-adjointness in fundamentally different manner. Consequently, no additional modification or justification is needed on the shape functions.

Conclusion

The energy equation of the analysis of hydrodynamic lubricated thrust bearing is solved by the finite element CBS method. The steady state of the three-dimensional temperature distribution across the lubricant film is obtained in a time-stepping manner.

The results from the case studies show that

1. The steady-state energy equation, with the convection terms present, can be solved by finite element CBS method.

2. Steady-state solution can be obtained with similar mesh density compared to FDM.

3. The CBS method, when used to solve the energy equation, can cope with lubrication with relatively lower lubricant viscosity and thinner lubricant film.

4. Typical convergence history of the CBS method shows that relatively faster converge behavior could be expected when the steady-state solution is about to be reached.

The CBS method is originally developed to deal with N-S equations. This work has extended its usability in lubrication analysis.