Accuracy and Grid Convergence of the Numerical Solution of the Energy Equation in Fluid Film Lubrication: Application to the 1D Slider

Abstract

The present work is focused on the numerical solution of the complete energy equation used in fluid film lubrication. The work was motivated by the fact that the complete energy equation has no analytical solution that can be used for validations. Its accuracy and computation time are related to the employed numerical method and to the grid resolution. The natural discretization method (NDM) applied on different grids is systematically compared with the spectral method (the Lobatto Point Colocation Method or LPCM) with different polynomial degrees. A one dimensional inclined slider is used for the numerical tests, and the energy equation is artificially decoupled from the Reynolds equation. This approach enables us to focus all the attention on the numerical solution of the energy equation. The results show that the LPCM is one or two orders of magnitude more efficient than the NDM in terms of computation time. The energy equation is then coupled with the Reynolds equation in a thermo-hydrodynamic analysis of the same 1D slider; the numerical results confirm again the efficiency of the LPCM. A thermo-hydrodynamic analysis of a two-lobe journal bearing is then presented as a practical application.

1. Introduction

The flow of thin lubricant films in journal and thrust bearings is most often described by the Reynolds equation of lubrication coupled with the energy transport equation. The numerical solution of these two coupled equations is a problem that has been solved for many decades [1]. However, it requires a computational effort that can render transient analyses very time consuming. Therefore, efficient solvers and adequate coupling strategies are of major importance to perform complex analyses in a reasonable amount of time.

If film rupture and reformation (traditionally designated as “cavitation”) are absent, and if the flow regime is laminar and isothermal, then the Reynolds equation is an elliptic linear differential equation. A direct solver can be used after its discretization in the thin film plane.

The energy equation contains convective transport terms, conductive transport terms across the thin film, and dissipative source terms arising from its coupling with the Reynolds equation. Its character is therefore not entirely elliptic, and its solver is different from the one used for the Reynolds equation. Moreover, the energy equation must also be discretized across the thin fluid film, and the number of discretization points in this direction must be large enough to capture temperature gradients near the walls. Therefore, the energy equation requires a substantially higher computational effort compared to the Reynolds equation. Developing an accurate and efficient solver for the energy equation is then an important step toward solving non-isothermal lubrication problems. The task is not simple because there is no analytical solution of the complete energy equation to be used for validations. For example, if the analytical solution of the laminar and isothermal thin film flow in a one dimensional (1D) slider can be used to validate the solver of the Reynolds equation, no similar solution exists for the energy equation. For validation, a numerical solver for the energy equation will have to be tested with different wall boundary conditions, and the numerical results must be checked for various grid densities. Such results are absent from the literature.

As mentioned above, the energy equation must also be discretized over the film thickness. The normal practice is to increase the number of discretization cells across the thin film until a grid independent numerical solution is reached, and this renders its solution time consuming. In the natural discretization methods (NDMs), such as the finite difference or finite volumes methods, the variation of the temperature between two cells is linear. This leads to a significant number of discretization cells across the thin film. An efficient approach for solving the energy equation was introduced by Elrod and uses the Legendre polynomials to describe the temperature variation across the thin film [2]. The coefficients of the Legendre polynomials expansion were first calculated by using a Galerkin approach. Later, Elrod used the collocation method at the Lobatto points; that is, the roots of the derivative of the highest order of the used Legendre polynomial. This method, known as the “Lobatto Points Collocation Method (LPCM)”, proved to be more efficient [3].

Despite its efficiency, the method is barely used, likely because of its complexity. Thus, Moraru [4] extended the LPCM to compressible films by describing the variation of the density across the thin film with Legendre polynomials. Feng and Kaneko [5] used the LPCM to solve the energy equation in aerodynamic foil bearings. Lehne et al. [6] presented a comprehensive review of the numerical solution strategies of the coupled Reynolds and energy equations, including the LPCM. Except for the cited references, the literature is scarce in references dealing with the LPCM.

This work presents a systematic comparison between the natural discretization method (NDM) of the energy equation and its LPCM approximation. In order to decouple the Reynolds and energy equation, the viscosity is supposed to be constant and the flow laminar. A one dimensional (1D) inclined slider is used for the numerical tests. The Reynolds equation then has an analytical solution, and the analysis is entirely focused on the energy equation. The results for the 1D slider compare both the number of points and the computational time required by the NDM and by the LPCM to obtain grid independent solutions. The results show how the NDM is excessively time consuming when high accuracy is sought, and the net economy of computational time brought by the LPCM approach.

A thermo-hydrodynamic analysis of the 1D slider with the coupled Reynolds and energy equation is then performed. The analysis highlights the same conclusions, namely the large superiority of the LPCM compared to NDM in terms of computational effort. Lastly, a two-lobe journal bearing with axial supply grooves is analyzed to demonstrate the accuracy that can be expected with simple, non-triggered boundary conditions for the energy equation.

2. The Numerical Solution of the Energy Equation Based on the Natural Discretization Method

The 1D inclined slider shown in Figure 1 is used in References [2,3] for the resolution of the Reynolds equation coupled with the energy equation. The current work is focused on the resolution of the energy equation decoupled from the Reynolds equation, for which an analytical solution is used. This gives an accurate numerical solution of the energy equation without the influence of the coupling with pressure and velocities affected by numerical uncertainties. The numerical data used is detailed in

Table 1. The data used for the 1D slider.

Physical Characteristics	Values	Physical Characteristics	Values
Density ${\rho }_0$	$800 \left[\mathrm{kg}/{\mathrm{m}}^3\right]$	Inlet thickness $h_1$	$1.8288e^{-4} \left[\mathrm{m}\right]$
Specific Heat Capacity $C_p$	$2000 \left[\mathrm{J}/\mathrm{kgK}\right]$	Outlet thickness $h_2$	$0.9144e^{-4} \left[\mathrm{m}\right]$
Thermal Conductivity $\lambda$	$0.14 \left[ \mathrm{W}/\mathrm{mK}\right]$	Slider length $L$	$0.18288 \left[\mathrm{m}\right]$
Dynamic Viscosity ${\eta }_0$	$0.081 \left[\mathrm{Pa}·\mathrm{s}\right]$	Ambient temperature $Ta$	$20 \left[°\mathrm{C}\right]$
Lower Wall Operating Speed $U$	$31.946 \left[\mathrm{m}/\mathrm{s}\right]$	Ambient pressure $Pa$	$1 \left[\mathrm{bar}\right]$

.

The conservative form of the energy equation for the lubrication thin film flow of an incompressible fluid writes [1]:

$${\rho }_0{C}_p\left[\frac{\partial \left(uT\right)}{\partial x}+\frac{\partial \left(vT\right)}{\partial y}\right]= \lambda \frac{{\partial }^{2}T}{\partial y^2}+{\eta }_0{\left(\frac{\partial u}{\partial y}\right)}^2$$

(1)

It contains the convection transport terms on its left hand side, and the temperature diffusion and dissipation terms on its right hand side. The terms on the right hand side are simplified according to the lubrication thin film assumption; that is, only derivatives across the thin film are taken into account.

The dimensionless coordinate $\overline{y}$ is introduced to take into account that the film thickness $h$ is not constant.

$$y=\overline{y}h\left(x\right)$$

(2)

Following this coordinate transformation, the energy Equation (1) for the 1D slider becomes:

$${\rho }_0{C}_p\left[\frac{\partial \left(uT\right)}{\partial x}-\frac{\overline{y}}{h}\frac{dh}{dx}\frac{\partial \left(uT\right)}{\partial \overline{y}}+\frac{1}{h}\frac{\partial \left(vT\right)}{\partial \overline{y}}\right]=\lambda \frac{1}{h^2}\frac{\partial ²T}{\partial \overline{y}²} + {\eta }_0\frac{1}{h^2}{\left(\frac{\partial u}{\partial \overline{y}}\right)}^2$$

(3)

The slider is discretized with 2D computational cells as depicted in Figure 2. Following the variable change described by Equation (2), the computational domain is rectangular and orthogonal. The computational cells have four plane faces, denoted by lower-case letters corresponding to their direction (e, w, n, and s) with respect to the central node P.

In the natural discretization numerical approach, the energy equation is solved using the finite volume method [7]. Equation (3) is therefore integrated over the control volumes corresponding to the 2D computational cells:

$$\begin{array}{cc}\hfill {\rho }_0{C}_p{\int }_{x_w}^{x_e}{\int }_{{\overline{y}}_s}^{{\overline{y}}_n}& \left[\frac{\partial \left(uT\right)}{\partial x}-\frac{\overline{y}}{h}\frac{dh}{dx}\frac{\partial \left(uT\right)}{\partial \overline{y}}+\frac{1}{h}\frac{\partial \left(vT\right)}{\partial \overline{y}}\right]dxd\overline{y}\hfill \\ & =\lambda \frac{1}{h_P^2}{\int }_{x_w}^{x_e}{\int }_{{\overline{y}}_s}^{{\overline{y}}_n}\frac{{\partial }^{2}T}{\partial {\overline{y}}^2}dxd\overline{y} + {\eta }_0\frac{1}{h_P^2}{\left(\frac{\partial u}{\partial \overline{y}}\right)}_P{}^2\Delta x\Delta \overline{y}\end{array}$$

(4)

For example, the convection transport terms in the x direction can be expressed as:

$${\int }_{x_w}^{x_e}{\int }_{{\overline{y}}_s}^{{\overline{y}}_n}\frac{\partial \left(uT\right)}{\partial x}dxd\overline{y}=\left[{\left(uT\right)}_e-{\left(uT\right)}_w\right]\Delta \overline{y}$$

(5)

An upwind discretization technique is used for the convective transport terms to avoid numerical instabilities [7]. For example, at the east face of the control volume, the temperature $T_e$ is up-winded based on the fluid flow direction. Mathematically, this can be expressed as:

$${\left(uT\right)}_e=u_e\left[-T_P\frac{1-sgn\left(u_e\right)}{2}+T_E\frac{1-sgn\left(u_e\right)}{2}+T_P\right]$$

(6)

This yields the discretization form of Equation (3):

$$\begin{array}{cc}\hfill (a_e+a_w+a_n+& a_s)T_P-a_e{T}_E-a_w{T}_W-a_n{T}_N-a_s{T}_S-a_P\left(T_N-T_S\right)\hfill \\ & \hfill =\frac{\lambda }{{\rho }_0{C}_p{h}_P^2}\left(T_N-2T_P+T_S\right)\frac{\Delta x}{\Delta \overline{y}}+\frac{{\eta }_0}{{\rho }_0{C}_p{h}_P^2} {\left({\frac{\partial u}{\partial y}|}_P\right)}^2 \Delta x\Delta \overline{y}\end{array}$$

(7)

where:

$$\begin{gathered} a_e =\frac{u_e\left[sgn\left(u_e\right)-1\right]}{2}\Delta \overline{y}; a_w=\frac{u_w\left[sgn\left(u_w\right)-1\right]}{2}\Delta \overline{y}; \\ {a}_n=\frac{v_n\left[sgn\left(v_n\right)-1\right]}{2}\frac{\Delta x}{h_P}; a_s=\frac{v_s\left[sgn\left(v_s\right)-1\right]}{2}\frac{\Delta x}{h_P}; \\ {a}_P=\frac{u_P{\overline{y}}_P}{2}{\frac{dh}{dx}|}_P\frac{\Delta x}{h_P} \end{gathered}$$

(8)

As stated in the introduction, the goal of the present work is to investigate the accuracy of the numerical solution of the energy equation given by the NDM and LPCM methods. Therefore, in order to avoid the uncertainties introduced by the numerical solution of the Reynolds equation, its analytic solution for the 1D slider is used.

Several other simplifying assumptions are needed for decoupling the Reynolds and the energy equations. Thus, the viscosity is considered constant and the flow regime laminar. Following the analytical solution, the pressure and the velocity in the 1D slider are [1]:

$$p\left(x\right)=\frac{6{\eta }_{0}UL}{h_1-h_2}\left(h^{-1}-\frac{h_1{h}_2}{h_1+h_2}{h}^{-2}-\frac{1}{h_1+h_2}\right)$$

(9)

$$u\left(x,\overline{y}\right)=h^2\frac{\overline{y}\left(\overline{y}-1\right)}{2{\eta }_0}\frac{dp}{dx}+\left(1-\overline{y}\right)U$$

(10)

The velocity $v$ component across the thin film is deduced by integrating the continuity Equation (11) over the film thickness:

$$\frac{\partial u}{\partial x}+\frac{1}{h}\frac{\partial v}{\partial \overline{y}}-\frac{\overline{y}}{h}\frac{dh}{dx}\frac{\partial u}{\partial \overline{y}}=0$$

(11)

With the boundary conditions $v\left(x,\overline{y}=0\right)=0$ and $v\left(x,\overline{y}=1\right)=0$, this yields:

$$v\left(x,\overline{y}\right)=-h{\int }_0^{\overline{y}}\frac{\partial u}{\partial x}d\overline{y}+\frac{dh}{dx}{\int }_0^{\overline{y}}\overline{y}\frac{\partial u}{\partial \overline{y}}d\overline{y}$$

(12)

It is to be underlined that in the following approach only the boundary condition at $\overline{y}=0$ is used, while the second boundary condition at $\overline{y}=1$ is used to check the accuracy of the numerical integration.

The discretized system given by Equation (7) is solved using any numerical procedure for a linear system with a positive definite matrix. Grids with 80 equidistant cells in the x direction and different numbers of cells across the film thickness were used. The number of equidistant grid cells across the film thickness are Ny = 10, 20, 40, 50, 80, 100, and 160. In a first test, the ambient temperature of 20 °C was imposed on the lower and upper walls and at the inlet. The dimensionless wall temperature gradient was monitored (This is equivalent to monitoring the wall heat fluxes). This dimensionless wall temperature gradient is defined by:

$$\frac{\partial \overline{T}\left(x,\overline{y}\right)}{\partial \overline{y}}=\frac{1}{T_a}\frac{\partial T\left(x,\overline{y}\right)}{\partial \overline{y}}$$

(13)

where $T_a=20 °\mathrm{C}$ is a reference temperature.

The results are depicted in Figure 3 and Figure 4, and show that the curves become superposed starting with Ny = 40 cells. An estimation of the accuracy is obtained by using the relative difference of the wall temperature gradients which is defined in Equation (14):

$$\epsilon =\frac{\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left(\frac{\partial {\overline{T}}_{Ny}\left(x_i,\overline{y}\right)}{\partial \overline{y}}-\frac{\partial {\overline{T}}_{ref}\left(x_i,\overline{y}\right)}{\partial \overline{y}}\right)}^2}}{\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left(\frac{\partial {\overline{T}}_{ref}\left(x_i,\overline{y}\right)}{\partial \overline{y}}\right)}^2}} for \overline{y}=0 or \overline{y}=1$$

(14)

where ${\overline{T}}_{ref}$ is the solution obtained with the finest grid Ny = 160.

The variation of these errors with the number of grid cells across the film thickness is depicted in Figure 5a, and the computational time is depicted in Figure 5b. A rapid increase in the computational time with the grid cells must be underlined. For example, the calculation case with Ny = 80 grid cells required 1.844 s, while the computational time for the case with Ny = 160 cells was one order of magnitude higher. However, for the purpose of accuracy, the mesh with Ny = 160 was chosen as the reference solution in the following and its numerical values of the wall temperature gradients are given in

Table A1. Result of the all cases obtained by the best NDM mesh.

	$\overline{\mathit{X}} [-]$	0.00625	0.11900	0.24400	0.36900	0.49400	0.61900	0.74400	0.86900	0.99400
Case 1	$\mathit{d}\overline{\mathit{T}}/\mathit{d}\overline{\mathit{y}}\_\mathit{L}\mathit{o}\mathit{w}\mathit{e}\mathit{r}\mathit{W}$	1.55718	5.21441	6.77197	7.44975	7.53430	7.14248	6.36616	5.34304	4.35192
	$\mathit{d}\overline{\mathit{T}}/\mathit{d}\overline{\mathit{y}}\_\mathit{U}\mathit{p}\mathit{p}\mathit{e}\mathit{r}\mathit{W}$	−0.13495	−0.96439	−1.92252	−3.09524	−4.61156	−6.64528	−9.46179	−13.49330	−19.48457
Case 2	$\mathit{d}\overline{\mathit{T}}/\mathit{d}\overline{\mathit{y}}\_\mathit{L}\mathit{o}\mathit{w}\mathit{e}\mathit{r}\mathit{W}$	1.43830	5.35156	7.30481	8.42416	8.87026	8.57660	7.37678	5.21935	3.57358
	$\mathit{d}\overline{\mathit{T}}/\mathit{d}\overline{\mathit{y}}\_\mathit{U}\mathit{p}\mathit{p}\mathit{e}\mathit{r}\mathit{W}$	−2.46419	−2.34816	−2.24448	−2.22769	−2.42486	−3.13238	−5.03152	−10.04817	−24.05486
Case 3	$\mathit{d}\overline{\mathit{T}}/\mathit{d}\overline{\mathit{y}}\_\mathit{L}\mathit{o}\mathit{w}\mathit{e}\mathit{r}\mathit{W}$	0.830646	3.894458	5.609075	6.852488	7.746559	8.195229	7.845554	5.840858	3.052226
	$\mathit{d}\overline{\mathit{T}}/\mathit{d}\overline{\mathit{y}}\_\mathit{U}\mathit{p}\mathit{p}\mathit{e}\mathit{r}\mathit{W}$	−5.17655	−4.86225	−4.47603	−4.02043	−3.49895	−2.97106	−2.75667	−4.76198	−25.4916
Case 4	$\mathit{T}\_\mathit{L}\mathit{o}\mathit{w}\mathit{e}\mathit{r}\mathit{W}$[°C]	20.276	22.293	24.839	27.649	30.669	33.738	36.470	38.132	38.967
	$\mathit{d}\overline{\mathit{T}}/\mathit{d}\overline{\mathit{y}}\_\mathit{U}\mathit{p}\mathit{p}\mathit{e}\mathit{r}\mathit{W}$	−0.41141	−0.33314	−0.28279	−0.32754	−0.59665	−1.37867	−3.38145	−8.58824	−23.9895
Case 5	$\mathit{d}\overline{\mathit{T}}/\mathit{d}\overline{\mathit{y}}\_\mathit{L}\mathit{o}\mathit{w}\mathit{e}\mathit{r}\mathit{W}$	0.938225	5.39037	7.32993	8.39752	8.86711	8.55326	7.32938	5.15078	3.62624
	$\mathit{d}\overline{\mathit{T}}/\mathit{d}\overline{\mathit{y}}\_\mathit{U}\mathit{p}\mathit{p}\mathit{e}\mathit{r}\mathit{W}$	−0.42745	−0.34021	−0.28590	−0.32154	−0.59344	−1.36122	−3.32850	−8.46418	−23.39503

in Appendix B.

3. The Numerical Solution of the Energy Equation Based on the Lobatto Points Collocation Method

The Lobatto Point Collocation Method (LPCM) is based on the approximation of the temperature with Legendre polynomials across the film thickness. Because Legendre polynomials are defined on the interval $\left[-1, 1\right]$, the following coordinate transformation is used:

$$y=\frac{\left(\zeta +1\right)h}{2}$$

(15)

where $\zeta$ is the new dimensionless coordinate across the film thickness.

Following the coordinate transformation, the energy Equation (1) becomes:

$$\rho C_p\left[\frac{\partial \left(uT\right)}{\partial x}-\frac{\left(\zeta +1\right)}{h}\frac{dh}{dx}\frac{\partial \left(uT\right)}{\partial \zeta }+\frac{2}{h}\frac{\partial \left(vT\right)}{\partial \zeta }\right]=\lambda \frac{4}{h^2}\frac{\partial ²T}{\partial \zeta ²} + {\eta }_0\frac{4}{h^2}{\left(\frac{\partial u}{\partial \zeta }\right)}^2$$

(16)

For an incompressible lubricant with constant viscosity, only the variable temperature $T$ is approximated across the film thickness:

$$T\left(x,\zeta \right)={\displaystyle \sum }_{j=0}^N{\widehat{T}}_j\left(x\right)P_j\left(\zeta \right)$$

(17)

where $N$ is the highest order of the Legendre polynomials, $P_j\left(\zeta \right)$ is the Legendre polynomial of degree j, and ${\widehat{T}}_j\left(x\right)$ is the polynomial coefficients of temperature of degree j.

This description holds for every point in the 2D computational domain. Compared to the NDM, which computes directly the temperature by solving the energy equation discretized over the film thickness, the LPCM calculates the polynomial coefficients of temperature ${\widehat{T}}_j$. Based on the work of Mahner et al. [6], the coefficients of the Legendre polynomials expansion can be obtained using different methods, but it is accepted that the most reliable is the collocation method using the Lobatto points. The Lobatto points are the roots of the derivative of the Legendre polynomial of the highest degree (i.e., the roots of $d{P}_N/d\zeta$).

For a given position in the x direction, the temperature is replaced by its approximation from Equation (17) across the fluid film, and the energy Equation (16) is enforced to hold true for each Lobatto point, ${\zeta }_j,j\in \left[1,\dots ,N-1\right]$. This leads to N − 1 partial differential equations with the unknown ${\widehat{T}}_j$. The boundary conditions are applied at the two walls, $\zeta =-1$ and $\zeta =1$, which leads to the other two equations. In total, a system of N + 1 equations for the N + 1 unknown ${\widehat{T}}_j$ is obtained.

Figure 6 and Figure 7 depict the dimensionless wall temperature gradient obtained with an increasing degree of the Legendre polynomials. The difference between the NDM solution and the LPCM solution at the lower wall (Figure 6) in the inlet section is reduced by increasing the degree of the Legendre polynomial.

The reference results used for comparison are given by the NDM with Ny = 160 grid cells. The relative difference between the reference results and the wall temperature gradient obtained with the LPCM is defined as:

$$\epsilon =\frac{\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left(\frac{2\partial {\overline{T}}_{N_{Lobbato}}\left(x_i,\zeta \right)}{\partial \zeta }-\frac{\partial {\overline{T}}_{re{f}_{NDM}}\left(x_i,\overline{y}\right)}{d\overline{y}}\right)}^2}}{\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left(\frac{d{\overline{T}}_{re{f}_{NDM}}\left(x_i,\overline{y}\right)}{d\overline{y}}\right)}^2}}$$

(18)

where $\zeta =1$ or $\zeta =-1$ for LPCM, and $\overline{y}=1$ or $\overline{y}=0$ for NDM.

The relative difference is depicted in Figure 8a. This relative difference decreases rapidly with the increase of the degree of the Legendre polynomial. However, the lower and upper wall results show different trends. A higher degree of the Legendre polynomial is required to reach a grid independent solution at the lower wall. This is due to the more predominant thermal effect at the driving wall; that is, the lower wall.

The computational effort is depicted in Figure 8b. Compared with the reference solution, the computational time of the LPCM solution is divided by ten, since only a limited expansion of the Legendre polynomials is necessary to obtain grid independent results.

For comparison, the temperature fields obtained with the LPCM (N = 12) and the NDM (Ny = 160) are depicted in Figure A1 and Figure A2 in Appendix A, respectively. The obtained polynomial coefficients of temperature of degree j, ${\widehat{T}}_j\left(\overline{X}\right)$, are given in

Table A2. Coefficients of Legendre polynomials ${\widehat{T}}_j\left(\overline{X}\right)$ describing temperature variations for Case 1, N = 12.

	$\overline{\mathit{X}}[-]$	0.00625	0.11900	0.24400	0.36900	0.49400	0.61900	0.74400	0.86900	0.99400
jth Degree
1		20.39	23.98	27.32	30.51	33.61	36.64	39.70	42.97	46.43
2		−0.09	−1.09	−1.70	−1.73	−1.04	0.54	3.24	7.36	12.66
3		−0.16	−2.41	−5.04	−7.85	−10.74	−13.60	−16.28	−18.67	−20.42
4		−0.07	−0.07	0.28	0.43	0.09	−1.00	−3.10	−6.52	−11.08
5		−0.10	−0.99	−1.76	−2.35	−2.77	−3.14	−3.62	−4.55	−6.18
6		0.03	0.78	1.18	1.20	0.93	0.47	−0.14	−0.90	−1.72
7		−0.04	−0.37	−0.43	−0.30	−0.11	0.06	0.17	0.20	0.12
8		0.06	0.32	0.24	0.12	0.03	0.00	0.00	0.05	0.12
9		−0.04	−0.19	−0.10	−0.02	0.02	0.03	0.03	0.04	0.05
10		0.04	0.07	0.01	−0.01	−0.01	−0.01	0.00	0.00	0.01
11		−0.04	−0.02	0.01	0.01	0.01	0.00	0.00	0.00	0.00
12		0.03	0.00	−0.01	0.00	0.00	0.00	0.00	0.00	0.00
13		−0.02	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00

in Appendix B. The present analysis indicates that the LPCM results are more accurate for the same grid density.

4. Further Comparison of Numerical Results Obtained by NDM and LPCM

The previous results were obtained for a linear converging 1D slider with an inlet/outlet film thickness ratio h₁/h₂ = 2, imposed wall temperatures (Dirichlet boundary condition), and completely decoupled from the Reynolds equation. Other simple cases of the converging 1D slider were also investigated and bring interesting conclusions.

4.1. Different Geometrical Configurations of the 1D Slider

Two other inlet/outlet film thickness ratios, h₁/h₂ = 4 and h₁/h₂ = 8, were investigated while keeping the same wall temperature boundary conditions and the decoupled Reynolds equation. The increased inlet/outlet film thickness ratio leads to a slower convergence of the NDM results with the y grid refinements. The best NDM solution obtained with Ny = 160 was considered as the reference results, and the number of computational cells in the x direction was kept the same (Nx = 80 cells).

The results obtained for h₁/h₂ = 4 are depicted in Figure 9, Figure 10 and Figure 11. Again, the lower and upper wall results show different trends. The upper wall temperature gradient reaches grid convergence starting with the Legendre polynomials of 11 degrees, while the resolution of the lower wall gradient needs N = 16 Legendre polynomials. However, Figure 11b shows that even when using a high number of Legendre polynomials, the LPCM requires one order of magnitude less computational time than the NDM for the same accuracy.

The results obtained for h₁/h₂ = 8 are depicted in Figure 12, Figure 13 and Figure 14. The same remarks as for h₁/h₂ = 4 can be drawn, except for the fact that this time, a higher order of the Legendre polynomials was needed for the grid convergence solution, and this is related to the increased inlet/outlet film thickness ratio. Figure 14b shows that the computational time of the NDM with Ny = 160 increases for this calculation case. It can be seen from Figure 11b and Figure 14b that the computational effort of the NDM increases one order of magnitude with increasing the ratio h₁/h₂ from 4 to 8. This was not the case for the LPCM, which required the same computational time for these two cases, one or two orders of magnitude lower than the NDM. Therefore, LPCM remains largely superior to NDM in terms of the computational time.

The temperature fields obtained with the LPCM (N = 16) and the NDM (Ny = 160) are depicted in Figure A3, Figure A4, Figure A5 and Figure A6 in Appendix A and its numerical values of the wall temperature gradients are given in Table A1 in Appendix B.

4.2. Different Thermal Boundary Conditions Applied to 1D Slider

The previous cases dealt with imposed wall temperatures, while the wall temperature gradient was a calculation result. In the following test, the lower wall of the 1D slider was adiabatic; that is, ${\left(\partial T/\partial \mathrm{y}\right)}_{\mathrm{LowerW}}=0$ while the temperature of the upper wall is ${\mathrm{T}}_{\mathrm{UpperW}}=30 °\mathrm{C}$. The film thickness ratio is h₁/h₂ = 4 and the inlet temperature ${\mathrm{T}}_{\mathrm{inlet}}=20 °\mathrm{C}$. For the LPCM, 10 and 14 Lobatto points were used and the results were compared to the NDM (Nx = 160 cells and Ny = 160 cells, all equidistant). Figure 15 and Figure 16 show the consistent resolution of the upper wall temperature gradient and of the lower wall temperature.

A different test case consists of imposing different wall temperatures, ${\mathrm{T}}_{\mathrm{UpperW}}=30 °\mathrm{C}$, ${\mathrm{T}}_{\mathrm{LowerW}}=20 °\mathrm{C}$, and ${\mathrm{T}}_{\mathrm{inlet}}=20 °\mathrm{C}$. The LPCM was performed with 10 and 14 Lobatto points, and the previous resolution was used for the NDM. The results are depicted in Figure 17 and Figure 18, and show that discrepancies exist between the NDM solution and the LPCM solution at the lower wall (Figure 18) in the inlet section. Again, this difference is reduced with the increase of the degree of the Legendre polynomials. It is necessary to highlight that the LPCM method approximates the temperature gradient at the walls based on the polynomial coefficients of temperature over the film thickness, while the NDM method calculates the temperature gradient with finite differences over a limited number of boundary cells. Thus, the two methods have different approximation orders. In some extreme configurations, it is necessary to increase the degree of the Legendre polynomial in order to reach the grid independent solution.

The temperature fields of these two last cases are depicted in Figure A7, Figure A8, Figure A9 and Figure A10 of Appendix A and its numerical values of the wall temperatures and temperature gradients are given in Table A1 in Appendix B.

5. Results for the Energy Equation Coupled with the Reynolds Equation

The above presented results showed the grid refinements required for obtaining accurate solutions of the energy equation when decoupled from the Reynolds equation. In reality, the viscosity of a liquid lubricant depends strongly on the temperature, and therefore decoupling is very difficult. For this reason, the data found in the literature always considered the complete thermo-hydrodynamic analysis.

This thermo-hydrodynamic (THD) coupled analysis was subsequently performed for the 1D slider case. The results of the LPCM were compared with the data from Ref. [3]. Their grid convergence and the computational time were analyzed by comparisons with NDM results. This step completes the validation of the decoupled energy equation described in the previous sections.

A temperature dependent viscosity following an exponentially decaying law $\eta \left(T\right)=0.13885e^{-0.045\left(T-T_{ambient}\right)}$ was used to replace the constant viscosity ${\eta }_0$. The rest of the geometrical and physical parameters were the same as in the previous case. The inlet and wall temperatures were equal to the ambient reference temperature. As for the variations of the temperature across the thin film, the terms in the generalized Reynolds equation with variable viscosity were also discretized by Legendre polynomials, and the coefficients were obtained by collocation at the Lobatto points.

As for Ref. [3], the computational domain of the 1D slider was discretized using Nx = 30 equidistant cells in the main flow direction and 10 Lobatto points across the film thickness. The Reynolds and the energy equations were numerically solved following a segregated approach.

Figure 19a compares the pressure variation in the 1D slider obtained with the LPCM and the one given in Ref. [3]. Figure 19b presents the variation of the outlet temperature difference across the film thickness. Both the pressure and the temperature variations obtained with the LPCM are in agreement with Ref. [3].

This 1D non-isothermal slider was also used to compare the efficiency of the LPCM with the NDM. Several calculations were performed with the NDM method in order to check the grid convergence and for obtaining results that could serve as a reference. These tests used seven grid refined in the y direction (Ny = 10, 20, 40, 60, 80, 100, and 120 equally spaced control volumes), while a constant number of 30 control volumes was used in the $x$ direction. The relative difference ${\epsilon }_K$ in terms of wall temperature gradients between two successive grids is defined as:

$${\epsilon }_K=\frac{\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left(\frac{\partial T_K\left(y\right)}{\partial y}-\frac{\partial T_{K+1}\left(y\right)}{\partial y}\right)}^2}}{\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left(\frac{\partial T_{K+1}\left(y\right)}{\partial y}\right)}^2}}$$

(19)

where the subscript K indicates the grid refinement level in the y direction.

Figure 20a shows that a minimum number of 40 cells in the y direction is necessary to reach a satisfactory grid-independent solution. The computational time is depicted in Figure 20b. The solution obtained by the NDM with 120 equally spaced cells over the film thickness is considered as the reference for the 1D thermo-hydrodynamic slider test case.

The LPCM results obtained with different numbers of Lobatto points were compared with the reference NDM solution. In Figure 21a, the relative difference drops rapidly and remains below 1%, starting with 13 Lobatto points. Figure 21b shows that the computational time for the LPCM does not exceed 2 s, while the reference method takes about 54 s.

6. Example of a Two-Lobe Journal Bearing

The coupled Reynolds and energy equations represent a solver for thermo-hydrodynamic problems in lubricated journal and thrust bearings. However, pure validations are more difficult in the context of lubricated bearings because of the additional effects, such as the film rupture and reformation (cavitation), which must be modeled and dealt with Ref. [8,9].

Not the least, several user-defined heat transfer parameters must be specified, and this is mainly done on a trial and error basis. The uncertainties in this kind of problem may therefore be quite important, and comparisons with experimental data are not exactly validations of the numerical approaches.

For the completeness of the present work, the numerical analysis of a real journal bearing is presented in the following. The recent experimental results published by Giraudeau et al. [10] in 2016 for a two-lobe journal bearing with axial supply grooves were used for comparing numerical predictions with experimental results. The length of the tested bearing was 68.4 mm and its diameter was 100 mm. The radial assembly clearance was 68 $\mathsf{\mu }\mathrm{m}$, while the radial bearing clearance was 143 $\mathsf{\mu }\mathrm{m}$. The bearing was lubricated by an ISO VG 46 oil supplied at a constant pressure of 0.17 MPa and a temperature of 43 °C. The following oil characteristics were used for the calculations: $\rho =850 \mathrm{kg}/{\mathrm{m}}^3$, $C_p=2000 \mathrm{J}/\left(\mathrm{kgK}\right)$, and $\lambda =0.13 \mathrm{W}/\left(\mathrm{mK}\right)$. The viscosity of the oil was $0.0416 \mathrm{Pa}·\mathrm{s}$ at $40 °\mathrm{C}$ and $0.0191 \mathrm{Pa}·\mathrm{s}$ at $60 °\mathrm{C}$. The variation is described by an exponentially decaying law. The angular speed was 3500 rpm and the static load was 6 kN. Only the results for the lower loaded lobe are shown here. Cavitation was dealt with by using the algorithm introduced in Ref. [11]. This algorithm uses a free boundary formulation of the incompressible Reynolds equation. The problem was then solved with an efficient solver based on the Fischer–Burmeister form, Newton algorithm, and Shur’s complement.

The shaft was considered to have a constant temperature estimated from experiments, while adiabatic wall conditions were imposed on the bushing. The computational domain was discretized using 32 $\times$ 16 cells in circumferential and axial directions, while 11 Lobatto points were used to describe the temperature variation across the fluid film.

Figure 22a,b depicts the pressure and the temperature variations in the circumferential direction of the bearing mid plane. The predicted pressures show good agreement with the measurements, while the predicted temperatures show a reasonable agreement. The quality of the prediction could be improved by including the thermal deformation of the bushing and refining the boundary conditions for the energy equation.

7. Conclusions

The present work focused on the numerical solution of the energy equation in very simple test cases as part of thin fluid films lubrication. The work was motivated by the fact that the complete energy equation has no analytical solution that can be used for comparisons in validations. Therefore, an energy equation decoupled from the Reynolds equation in the case of a 1D slider was imagined, for which the analytical solution of the velocity field was used. The reference values used in all cases presented are detailed in Appendix B. The second concern was an efficient solver for the energy equation. Two numerical methods were compared: The NDM based on a finite volume discretization, and the LPCM based on a Legendre polynomial approximation of the temperature across the film thickness. The LPCM proved to be one or two orders of magnitude more efficient than the NDM in terms of computation time. This is not negligible if one takes into account that when coupled, the Reynolds and the energy equations are numerically solved in a segregated and iterative manner. The thermo-hydrodynamic calculation of a 1D slider (i.e., for coupled Reynolds and energy equations) confirms this conclusion.

Comparisons of the LPCM thermo-hydrodynamic analyses with the experimental results obtained in a two lobe journal bearing were also presented. However, they are not considered as real validations of the numerical model, since very simple and intuitive boundary conditions of the energy equation were used. Alternatively, the results presented for the 1D slider may be used for validating the first development steps of any solver based on the energy equation for thin film flows.