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

Silun Zhang, Mohamed-Amine Hassini, Mihai Arghir

► To cite this version:

Silun Zhang, Mohamed-Amine Hassini, Mihai Arghir. Accuracy and Grid Convergence of the Numerical Solution of the Energy Equation in Fluid Film Lubrication: Application to the 1D Slider. Lubricants, 2018, 6 (4), pp.95. 10.3390/lubricants6040095. [https://doi.org/10.3390/lubricants6040095](https://doi.org/10.3390/lubricants6040095). [hal-02291733](https://HAL.inistre.fr/hal-02291733)

HAL Id: hal-02291733

https://hal.science/hal-02291733v1

Submitted on 4 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

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

Silun Zhang 1,2, Mohamed-Amine Hassini 1, and Mihai Arghir 2,*

Received: 13 September 2018; Accepted: 22 October 2018; Published: 26 October 2018

<figure><img src="image_3.png" /></figure>

check for
updates

**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.

**Keywords:** energy equation; numerical solution; lobatto point quadrature; legendre polynomials

## 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

<footnote n="1">1 ERMES-IMSIA, EDF Lab Paris-Saclay, 91120 Palaiseau, France; silun-s.zhang@edf.fr (S.Z.); mohamed-amine.hassini@edf.fr (M.-A.H.)</footnote>

<footnote n="2">2 Institut Pprime, UPR CNRS 3346, Université de Poitiers, ISAE ENSMA, 86962 Chasseneuil Futuroscope, France</footnote>

<footnote>* Correspondence: mihai.arghir@univ-poitiers.fr; Tel.: +33-(0)549496540</footnote>
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.
<figure><img src="image_4.png" /><figcaption>Figure 1. The 1D inclined slider.</figcaption></figure>

<table><caption>Table 1. The data used for the 1D slider.</caption><thead><tr><th>Physical Characteristics</th><th>Values</th><th>Physical Characteristics</th><th>Values</th></tr></thead><tbody><tr><td>Density ρ<sub>0</sub></td><td>800 [kg/m<sup>3</sup>]</td><td>Inlet thickness h<sub>1</sub></td><td>1.8288e-4 [m]</td></tr><tr><td>Specific Heat Capacity C<sub>p</sub></td><td>2000 [J/kgK]</td><td>Outlet thickness h<sub>2</sub></td><td>0.9144e-4 [m]</td></tr><tr><td>Thermal Conductivity λ</td><td>0.14 [W/mK]</td><td>Slider length L</td><td>0.18288 [m]</td></tr><tr><td>Dynamic Viscosity η<sub>0</sub></td><td>0.081 [Pa·s]</td><td>Ambient temperature T<sub>a</sub></td><td>20 [°C]</td></tr><tr><td>Lower Wall Operating Speed U</td><td>31.946 [m/s]</td><td>Ambient pressure P<sub>a</sub></td><td>1 [bar]</td></tr></tbody></table>

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(uT)}{\partial x} + \frac{\partial(vT)}{\partial y} \right] = \lambda \frac{\partial^2 T}{\partial y^2} + \eta_0 \left( \frac{\partial u}{\partial y} \right)^2 \quad (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 $\bar{y}$ is introduced to take into account that the film thickness $h$ is not constant.

$$ y = \bar{y}h(x) \quad (2) $$

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

$$ \rho_0 C_p \left[ \frac{\partial(uT)}{\partial x} - \frac{\bar{y} dh}{h dx} \frac{\partial(uT)}{\partial \bar{y}} + \frac{1}{h} \frac{\partial(vT)}{\partial \bar{y}} \right] = \lambda \frac{1}{h^2} \frac{\partial T}{\partial \bar{y}} + \eta_0 \frac{1}{h^2} \left( \frac{\partial u}{\partial \bar{y}} \right)^2 \quad (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.
<figure><img src="image_5.png" /><figcaption>Figure 2. A typical 2D computational cell and the notation used for a 2D grid.</figcaption></figure>

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:

$$ \rho_0 C_p \int_{x_w}^{x_e} \int_{\bar{y}_s}^{\bar{y}_n} \left[ \frac{\partial(uT)}{\partial x} - \frac{\bar{y}}{h} \frac{dh}{dx} \frac{\partial(uT)}{\partial \bar{y}} + \frac{1}{h} \frac{\partial(vT)}{\partial \bar{y}} \right] dx d\bar{y} \\ = \lambda \frac{1}{h_p^2} \int_{x_w}^{x_e} \int_{\bar{y}_s}^{\bar{y}_n} \frac{\partial^2 T}{\partial \bar{y}^2} dx d\bar{y} + \eta_0 \frac{1}{h_p^2} \left( \frac{\partial u}{\partial \bar{y}} \right)_p^2 \Delta x \Delta \bar{y} \quad (4) $$

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

$$ \int_{x_w}^{x_e} \int_{\bar{y}_s}^{\bar{y}_n} \frac{\partial(uT)}{\partial x} dx d\bar{y} = [(uT)_e - (uT)_w] \Delta \bar{y} \quad (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:

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

This yields the discretization form of Equation (3):

$$ \begin{aligned} & (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 (T_N - T_S) \\ & = \frac{\lambda}{\rho_0 C_p h_p^2} (T_N - 2T_P + T_S) \frac{\Delta x}{\Delta \bar{y}} + \frac{\eta_0}{\rho_0 C_p h_p^2} \left( \left. \frac{\partial u}{\partial y} \right)_P \right)^2 \Delta x \Delta \bar{y} \end{aligned} \quad (7) $$

where:

$$ \begin{aligned} a_e &= \frac{u_e[\operatorname{sgn}(u_e) - 1]}{2} \Delta \bar{y}; & a_w &= \frac{u_w[\operatorname{sgn}(u_w) - 1]}{2} \Delta \bar{y}; \\ a_n &= \frac{v_n[\operatorname{sgn}(v_n) - 1]}{2} \frac{\Delta x}{h_p}; & a_s &= \frac{v_s[\operatorname{sgn}(v_s) - 1]}{2} \frac{\Delta x}{h_p}; \\ a_P &= \frac{u_P \bar{y}_P}{2} \left. \frac{dh}{dx} \right|_P \frac{\Delta x}{h_p} \end{aligned} \quad (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(x) = \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) \quad (9)$$

$$u(x, \bar{y}) = h^2 \frac{\bar{y}(\bar{y} - 1)}{2\eta_0} \frac{dp}{dx} + (1 - \bar{y})U \quad (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 \bar{y}} - \frac{\bar{y}}{h} \frac{dh}{dx} \frac{\partial u}{\partial \bar{y}} = 0 \quad (11)$$

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

$$v(x, \bar{y}) = -h \int_{0}^{\bar{y}} \frac{\partial u}{\partial x} d\bar{y} + \frac{dh}{dx} \int_{0}^{\bar{y}} \frac{\bar{y}}{\bar{y}} \frac{\partial u}{\partial \bar{y}} d\bar{y} \quad (12)$$

It is to be underlined that in the following approach only the boundary condition at $\bar{y} = 0$ is used, while the second boundary condition at $\bar{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^{\circ}$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 \bar{T}(x, \bar{y})}{\partial \bar{y}} = \frac{1}{T_a} \frac{\partial T(x, \bar{y})}{\partial \bar{y}} \quad (13)$$

where $T_a = 20^{\circ}$C is a reference temperature.

The results are depicted in Figures 3 and 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 \bar{T}_{Ny}(x_i, \bar{y})}{\partial \bar{y}} - \frac{\partial \bar{T}_{ref}(x_i, \bar{y})}{\partial \bar{y}} \right)^2}}{\sqrt{\frac{1}{n} \sum_{i=1}^{n} \left( \frac{\partial \bar{T}_{ref}(x_i, \bar{y})}{\partial \bar{y}} \right)^2}} \quad \text{for } \bar{y} = 0 \text{ or } \bar{y} = 1 \quad (14)$$

where $\bar{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 in Appendix B.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Dimensionless temperature gradient at the lower wall for the natural discretization method (NDM) solution and imposed wall temperatures ($h_1/h_2 = 2$, $Nx = 80$ cells).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Dimensionless temperature gradient at the upper wall for the NDM solution and imposed wall temperatures ($h_1/h_2 = 2$, $Nx = 80$ cells).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** (a) Relative difference and (b) computational time of the NDM solution for y direction grid refinements ($h_1/h_2 = 2$, $Nx = 80$ cells).</figcaption></figure>
### 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 [-1, 1], the following coordinate transformation is used:

$$y = \frac{(\zeta + 1)h}{2} \qquad (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(uT)}{\partial x} - \frac{(\zeta + 1)}{h} \frac{dh}{dx} \frac{\partial(uT)}{\partial \zeta} + \frac{2}{h} \frac{\partial(vT)}{\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 \qquad (16)$$

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

$$T(x, \zeta) = \sum_{j=0}^{N} \hat{T}_{j}(x) P_{j}(\zeta) \qquad (17)$$

where N is the highest order of the Legendre polynomials, $P_j(\zeta)$ is the Legendre polynomial of degree j, and $\hat{T}_j(x)$ 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 $\hat{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 $dP_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 [1, \dots, N-1]$. This leads to $N-1$ partial differential equations with the unknown $\hat{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 $\hat{T}_j$ is obtained.

Figures 6 and 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.
<figure><img src="image_9.png" /><figcaption>Figure 6. Dimensionless temperature gradient at the lower wall for the Lobatto Points Collocation Method (LPCM) solution with imposed wall temperatures ($h_1/h_2 = 2$, $Nx = 80$ cells).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7. Dimensionless temperature gradient at the upper wall for the LPCM solution with imposed wall temperatures ($h_1/h_2 = 2$, $Nx = 80$ cells).</figcaption></figure>

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 \bar{T}_{N_{Lobatto}}(x_i, \zeta)}{\partial \zeta} - \frac{\partial \bar{T}_{refNDM}(x_i, \bar{y})}{\partial \bar{y}} \right)^2}}{\sqrt{\frac{1}{n} \sum_{i=1}^{n} \left( \frac{d \bar{T}_{refNDM}(x_i, \bar{y})}{d \bar{y}} \right)^2}} \quad (18) $$

where $\zeta = 1$ or $\zeta = -1$ for LPCM, and $\bar{y} = 1$ or $\bar{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 Figures A1 and A2 in Appendix A, respectively. The obtained polynomial coefficients of temperature of degree j, $\hat{T}_j(\bar{X})$, are given in Table A2 in Appendix B. The present analysis indicates that the LPCM results are more accurate for the same grid density.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** (a) Relative difference and (b) computational time of the LPCM solution for Y direction grid refinements ($h_1/h_2 = 2$, Nx = 80 cells)</figcaption></figure>

# 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_1/h_2 = 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_1/h_2 = 4$ and $h_1/h_2 = 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_1/h_2 = 4$ are depicted in Figures 9–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.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Dimensionless temperature gradient at the lower wall for the LPCM solution with imposed wall temperatures ($h_1/h_2 = 4$, $Nx = 80$ cells).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Dimensionless temperature gradient at the upper wall for the LPCM solution with imposed wall temperatures ($h_1/h_2 = 4$, $Nx = 80$ cells).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** (a) Relative difference and (b) computational time of the LPCM solution for Y direction grid refinements ($h_1/h_2 = 4$, $Nx = 80$ cells).</figcaption></figure>
The results obtained for $h_1/h_2 = 8$ are depicted in Figures 12–14. The same remarks as for $h_1/h_2 = 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 Figures 11b and 14b that the computational effort of the NDM increases one order of magnitude with increasing the ratio $h_1/h_2$ 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 Figures A3–A6 in Appendix A and its numerical values of the wall temperature gradients are given in Table A1 in Appendix B.

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Dimensionless temperature gradient at the lower wall for the LPCM solution with imposed wall temperatures ($h_1/h_2 = 8$, $Nx = 80$ cells).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Dimensionless temperature gradient at the upper wall for the LPCM solution with imposed wall temperatures ($h_1/h_2 = 8$, $Nx = 80$ cells).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 14.** (a) Relative difference and (b) computational time of the LPCM solution for Y direction grid refinements ($h_1/h_2 = 8$, $Nx = 80$ cells).</figcaption></figure>

## 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, $(\partial T/\partial y)_{LowerW} = 0$ while the temperature of the upper wall is $T_{UpperW} = 30$ °C. The film thickness ratio is $h_1/h_2 = 4$ and the inlet temperature $T_{inlet} = 20$ °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). Figures 15 and 16 show the consistent resolution of the upper wall temperature gradient and of the lower wall temperature.

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Dimensionless temperature gradient at the upper wall.</figcaption></figure>

A different test case consists of imposing different wall temperatures, $T_{UpperW} = 30$ °C, $T_{LowerW} = 20$ °C, and $T_{inlet} = 20$ °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 Figures 17 and 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.
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Temperature at the lower wall.</figcaption></figure>

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

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Dimensionless temperature gradient at the upper wall.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 18.** Dimensionless temperature gradient at the lower wall.</figcaption></figure>
## 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 $η(T) = 0.13885e^{-0.045(T-T_{ambient})}$ was used to replace the constant viscosity $η_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 $N_x = 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].

<figure><img src="image_22.png" /><figcaption>Figure 19. (a) Pressure variation in the 1D slider ($h_1/h_2 = 2$) coupled thermo-hydrodynamic (THD) solution, and (b) outlet temperature variation across the film thickness in the 1D slider ($h_1/h_2 = 2$) coupled THD solution.</figcaption></figure>

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 $ε_K$ in terms of wall temperature gradients between two successive grids is defined as:

$$ε_K = \frac{\sqrt{\frac{1}{n} \sum_{i=1}^{n} \left( \frac{\partial T_K(y)}{\partial y} - \frac{\partial T_{K+1}(y)}{\partial y} \right)^2}}{\sqrt{\frac{1}{n} \sum_{i=1}^{n} \left( \frac{\partial T_{K+1}(y)}{\partial y} \right)^2}} \quad (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.

<figure><img src="image_23.png" /><figcaption>**Figure 20.** (a) Relative difference $ε_K$ of the NDM for different grid refinements, and (b) computational time of the NDM for different numbers of cells, over the film thickness for the 1D THD slider.</figcaption></figure>

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.

<figure><img src="image_24.png" /><figcaption>**Figure 21.** (a) Relative difference $ε_N$ between the LPCM and the reference NDM $Ny = 120$, and (b) computational time of the LPCM compared to the reference NDM $Ny = 120$ results, for the 1D THD slider.</figcaption></figure>

## 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 µm, while the radial bearing clearance was 143 µ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: ρ = 850 kg/m³, C<sub>p</sub> = 2000 J/(kgK), and λ = 0.13 W/(mK). The viscosity of the oil was 0.0416 Pa·s at 40 °C and 0.0191 Pa·s at 60 °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 × 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.

<figure><img src="image_25.png" /><figcaption>**Figure 22.** (a) Comparison of measured pressures and current numerical results in the mid plane of the loaded lobe (3500 rpm, 6 kN load), and (b) comparisons of measured temperatures and current numerical results in the mid plane of the loaded lobe.</figcaption></figure>

**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.

**Author Contributions:** Conceptualization, M.A.; Methodology, S.Z. and M.A.; Software, S.Z.; Validation, S.Z. and M.-A.H.; Formal Analysis, M.A.; Data Curation, S.Z.; Writing—Original Draft Preparation, S.Z. and M.A.; Writing—Review & Editing, M.-A.H. and M.A.; Visualization, S.Z.; Supervision, M.A. and M.-A.H.; Project Administration, M.A.; Funding Acquisition, M.-A.H.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors thank Électricité de France (EDF) R&D for their financial support.

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

## Appendix A. Temperature Fields

$$ \text{Case 1: } h_1/h_2 = 2, T_{UpperW} = 20 \,^{\circ}\text{C}, T_{LowerW} = 20 \,^{\circ}\text{C}, T_{inlet} = 20 \,^{\circ}\text{C}, T_{outlet} = 20 \,^{\circ}\text{C}. $$

<figure><img src="image_26.png" /><figcaption>Figure A1. LPCM, N = 12.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure A2. NDM, Ny = 80.</figcaption></figure>
Case 2: $h_1/h_2 = 4$, $T_{UpperW} = 20$ °C, $T_{LowerW} = 20$ °C, $T_{inlet} = 20$ °C, $T_{outlet} = 20$ °C.

<figure><img src="image_28.png" /><figcaption>Figure A3. LPCM, N = 16.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Figure A4. NDM, Ny = 160.</figcaption></figure>
Case 3: $h_1/h_2 = 8$, $T_{UpperW} = 20$ °C, $T_{LowerW} = 20$ °C, $T_{inlet} = 20$ °C, $T_{outlet} = 20$ °C.

<figure><img src="image_30.png" /><figcaption>Figure A5. LPCM, N = 16.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Figure A6. NDM, Ny = 160.</figcaption></figure>
Case 4: $h_1/h_2 = 4$, $T_{UpperW} = 30$ °C, $(\partial T/\partial y)_{LowerW} = 0$, $T_{inlet} = 20$ °C, $T_{outlet} = 20$ °C.

<figure><img src="image_32.png" /><figcaption>Figure A7. LPCM, N = 16.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Figure A8. NDM, Ny = 160.</figcaption></figure>
Case 5: $h_1/h_2 = 4$, $T_{UpperW} = 30$ °C, $T_{LowerW} = 20$ °C, $T_{inlet} = 20$ °C, $T_{outlet} = 20$ °C.

<figure><img src="image_34.png" /><figcaption>Figure A9. LPCM, N = 16.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Figure A10. NDM, Ny = 160.</figcaption></figure>
## Appendix B. Reference Values

<table><caption>Table A1. Result of the all cases obtained by the best NDM mesh.</caption><thead><tr><th></th><th>X&#772; [-]</th><th>0.00625</th><th>0.11900</th><th>0.24400</th><th>0.36900</th><th>0.49400</th><th>0.61900</th><th>0.74400</th><th>0.86900</th><th>0.99400</th></tr></thead><tbody><tr><td rowspan="2">Case 1</td><td>d&#772;T/d&#772;y&#772;_LowerW</td><td>1.55718</td><td>5.21441</td><td>6.77197</td><td>7.44975</td><td>7.53430</td><td>7.14248</td><td>6.36616</td><td>5.34304</td><td>4.35192</td></tr><tr><td>d&#772;T/d&#772;y&#772;_UpperW</td><td>-0.13495</td><td>-0.96439</td><td>-1.92252</td><td>-3.09524</td><td>-4.61156</td><td>-6.64528</td><td>-9.46179</td><td>-13.49330</td><td>-19.48457</td></tr><tr><td rowspan="2">Case 2</td><td>d&#772;T/d&#772;y&#772;_LowerW</td><td>1.43830</td><td>5.35156</td><td>7.30481</td><td>8.42416</td><td>8.87026</td><td>8.57660</td><td>7.37678</td><td>5.21935</td><td>3.57358</td></tr><tr><td>d&#772;T/d&#772;y&#772;_UpperW</td><td>-2.46419</td><td>-2.34816</td><td>-2.24448</td><td>-2.22769</td><td>-2.42486</td><td>-3.13238</td><td>-5.03152</td><td>-10.04817</td><td>-24.05486</td></tr><tr><td rowspan="2">Case 3</td><td>d&#772;T/d&#772;y&#772;_LowerW</td><td>0.830646</td><td>3.894458</td><td>5.609075</td><td>6.852488</td><td>7.746559</td><td>8.195229</td><td>7.845554</td><td>5.840858</td><td>3.052226</td></tr><tr><td>d&#772;T/d&#772;y&#772;_UpperW</td><td>-5.17655</td><td>-4.86225</td><td>-4.47603</td><td>-4.02043</td><td>-3.49895</td><td>-2.97106</td><td>-2.75667</td><td>-4.76198</td><td>-25.4916</td></tr><tr><td rowspan="2">Case 4</td><td>T&#772;_LowerW [°C]</td><td>20.276</td><td>22.293</td><td>24.839</td><td>27.649</td><td>30.669</td><td>33.738</td><td>36.470</td><td>38.132</td><td>38.967</td></tr><tr><td>d&#772;T/d&#772;y&#772;_UpperW</td><td>-0.41141</td><td>-0.33314</td><td>-0.28279</td><td>-0.32754</td><td>-0.59665</td><td>-1.37867</td><td>-3.38145</td><td>-8.58824</td><td>-23.9895</td></tr><tr><td rowspan="2">Case 5</td><td>d&#772;T/d&#772;y&#772;_LowerW</td><td>0.938225</td><td>5.39037</td><td>7.32993</td><td>8.39752</td><td>8.86711</td><td>8.55326</td><td>7.32938</td><td>5.15078</td><td>3.62624</td></tr><tr><td>d&#772;T/d&#772;y&#772;_UpperW</td><td>-0.42745</td><td>-0.34021</td><td>-0.28590</td><td>-0.32154</td><td>-0.59344</td><td>-1.36122</td><td>-3.32850</td><td>-8.46418</td><td>-23.39503</td></tr></tbody></table>

<table><caption>Table A2. Coefficients of Legendre polynomials $\hat{T}_j(\bar{X})$ describing temperature variations for Case 1, $N = 12$.</caption><thead><tr><th>X&#772; [-]</th><th>0.00625</th><th>0.11900</th><th>0.24400</th><th>0.36900</th><th>0.49400</th><th>0.61900</th><th>0.74400</th><th>0.86900</th><th>0.99400</th></tr></thead><tbody><tr><td>1</td><td>20.39</td><td>23.98</td><td>27.32</td><td>30.51</td><td>33.61</td><td>36.64</td><td>39.70</td><td>42.97</td><td>46.43</td></tr><tr><td>2</td><td>-0.09</td><td>-1.09</td><td>-1.70</td><td>-1.73</td><td>-1.04</td><td>0.54</td><td>3.24</td><td>7.36</td><td>12.66</td></tr><tr><td>3</td><td>-0.16</td><td>-2.41</td><td>-5.04</td><td>-7.85</td><td>-10.74</td><td>-13.60</td><td>-16.28</td><td>-18.67</td><td>-20.42</td></tr><tr><td>4</td><td>-0.07</td><td>-0.07</td><td>0.28</td><td>0.43</td><td>0.09</td><td>-1.00</td><td>-3.10</td><td>-6.52</td><td>-11.08</td></tr><tr><td>5</td><td>-0.10</td><td>-0.99</td><td>-1.76</td><td>-2.35</td><td>-2.77</td><td>-3.14</td><td>-3.62</td><td>-4.55</td><td>-6.18</td></tr><tr><td>6</td><td>0.03</td><td>0.78</td><td>1.18</td><td>1.20</td><td>0.93</td><td>0.47</td><td>-0.14</td><td>-0.90</td><td>-1.72</td></tr><tr><td>7</td><td>-0.04</td><td>-0.37</td><td>-0.43</td><td>-0.30</td><td>-0.11</td><td>0.06</td><td>0.17</td><td>0.20</td><td>0.12</td></tr><tr><td>8</td><td>0.06</td><td>0.32</td><td>0.24</td><td>0.12</td><td>0.03</td><td>0.00</td><td>0.00</td><td>0.05</td><td>0.12</td></tr><tr><td>9</td><td>-0.04</td><td>-0.19</td><td>-0.10</td><td>-0.02</td><td>0.02</td><td>0.03</td><td>0.03</td><td>0.04</td><td>0.05</td></tr><tr><td>10</td><td>0.04</td><td>0.07</td><td>0.01</td><td>-0.01</td><td>-0.01</td><td>-0.01</td><td>0.00</td><td>0.00</td><td>0.01</td></tr><tr><td>11</td><td>-0.04</td><td>-0.02</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>12</td><td>0.03</td><td>0.00</td><td>-0.01</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>13</td><td>-0.02</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr></tbody></table>
References

1. Frêne, J. *Lubrification Hydrodynamique: Paliers et Butées*; Eyrolle: Paris, France, 1990.

2. Elrod, H.G.; Brewe, D.E. *Thermohydrodynamic Analysis for Laminar Lubricating Films*; NASA Technical Report No. 88845; Elsevier: Amsterdam, The Netherlands, 1986.

3. Elrod, H.G. Efficient Numerical Method for Computation of the Thermohydrodynamics of Laminar Lubricating Films. *ASME J. Tribol. 1991*, *113*, 506–511. [CrossRef]

4. Moraru, L. Numerical Predictions and Measurements in the Lubrication of Aeronautical Engine and Transmission Components. Ph.D. Thesis, The University of Toledo, Toledo, OH, USA, 2005.

5. Feng, K.; Kaneko, S. Thermohydrodynamic Study of Multiwound Foil Bearing Using Lobatto Point Quadrature. *ASME. J. Tribol. 2009*. [CrossRef]

6. Mahner, M.; Lehn, A.; Schweizer, B. Thermogas- and thermohydrodynamic simulation of thrust and slider bearings: Convergence and efficiency of different reduction approaches. *Tribol. Int.* **2015**, *93*, 539–554. [CrossRef]

7. Ferziger, J.H.; Peric, M. *Computational Methods for Fluid Dynamics*; Springer: Berlin/Heidelberg, Germany, 2001.

8. Elrod, H.G. A cavitation algorithm. *ASME J. Tribol. 1981*, *103*, 350–354. [CrossRef]

9. Dominique, B.; Aurélian, F. *Paliers Hydrodynamiques 1 and 2, Équations, Modèles Numériques Isothermes et Lubrification Mixte*; Lavoisier: Paris, France, 2011.

10. Giraudeau, C.; Bouyer, J.; Fillon, M.; Hélène, M.; Beaurain, J. Experimental Study of the Influence of Scratches on the Performance of a Two-Lobe Journal Bearing. *Tribol. Trans. 2017*, *60*, 942–955. [CrossRef]

11. Woloszynski, T.; Podsiadlo, P.; Stachowiak, G.W. Efficient Solution to the Cavitation Problem in Hydrodynamic Lubrication. *Tribol. Lett.* **2015**. [CrossRef]

<figure><img src="image_36.png" /></figure>

© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).