# Advanced Modeling of Lubricated Interfaces in General Curvilinear Grids

Suhaib Ardah  Department of Mechanical Engineering, Imperial College London,
London, SW7 2AZ, UK  Corresponding author. Email
<EMAIL_ADDRESS>Francisco J. Profito  Department of
Mechanical Engineering, Polytechnic School of the University of São Paulo, São
Paulo, Brazil  Tom Reddyhoff  Department of Mechanical Engineering, Imperial
College London, London, SW7 2AZ, UK  Daniele Dini  Department of Mechanical
Engineering, Imperial College London, London, SW7 2AZ, UK

###### Abstract

Tackling fluid-flow problems involving intricate surface geometries has been
the catalyst for a plethora of numerical investigations aimed at accommodating
curved complex boundaries. An example is the application of body-fitted
curvilinear coordinate transformation, where the one-to-one correspondence of
grid points from the physical to the computational domain is achieved. In
lubricated interfaces, such conversion is challenging due to the complex
governing equations in the mapped-grid, the numerical instabilities exhibited
by their non-linearities and the severity of operating conditions. The present
contribution proposes a Reynolds-based, finite volume fluid-structure
interaction (FSI) framework for solving thermal elastohydrodynamic lubrication
(TEHL) problems mapped onto non-orthogonal curvilinear grids in the
computational domain. We demonstrate how the strong conservation form of the
pertinent governing equations can be expressed in three-dimensional
curvilinear grids and discretised using finite volume method to ensure fluid-
flow conservation and enforce mass-conserving cavitation conditions. Numerical
and experimental benchmarks showcase the robustness and versatility of the
proposed framework to simulate a diverse range of lubrication problems, hence
achieving a predictive computational tool that would enable a shift towards
tribology-aware design.

###### keywords:

Thermal elastohydrodynamic lubrication, Body-fitted curvilinear coordinates,
Non-orthogonal grids, Finite volume method

##  1 Introduction

The utilisation of structured rectangular grids for numerically solving
Reynolds-type equations that govern the fluid flow behaviour in lubricated
contacts has been used extensively due to its straightforward discretisation
procedure, flexibility in imposing the boundary conditions, and lower
computational effort to solve the discrete system of equations. Nevertheless,
a poorly constructed grid may yield discretisation errors, particularly in
problems with complex geometries, and can therefore influence the accuracy of
the numerical solution. Moreover, mathematical models described in Cartesian
coordinate systems may suffer from weaknesses when modelling flow at solid
boundaries that are either curved or not oriented along the coordinate axes
[1, 2]. Although the Reynolds model is considered to be an appropriate
approximation of the full Navier-Stokes equations for solving film lubrication
problems, yet the normalisation technique which traditional Reynolds-based
approaches adopt aimed at imposing the hydrodynamic solution onto orthogonal
rectangular grids can often result in inaccurate solutions due to the
negligence of any curved features within the lubricated contact. This is more
profound near the inlet and outlet regions of the contact, which are usually
characterised by complex curvatures, as reported in [3], where it was shown
that CFD and Reynolds models yielded slightly different thermal solutions at
the contact boundaries. To overcome these limitations and therefore improve
the performance of fluid-structure interaction (FSI) solvers, commercial codes
such as those based on the Computational Fluid Dynamics (CFD) approach have
been extensively used, owing to their capabilities to capture complex geometry
details in the solution of the fluid flow problem. Another strategy to deal
with complicated geometries is mapping the physical flow domain onto a
computational structured domain using specific coordinate transformations or
body-fitted grids built via Computer-Aided Design (CAD) models [2].

The application and development of coordinate transformation methods,
particularly the technique of body-fitted curvilinear coordinate systems with
coordinate lines coincident with the boundaries of arbitrarily shaped bodies,
gained popularity in solving CFD problems because of their versatility in
simulating flow around arbitrary physical geometries. Such coordinate
transformation methods enable irregular physical domains to be mapped onto
rectangular computational domains where the governing equations can be solved
more efficiently, the discretisation errors minimised, and the boundary
conditions imposed more accurately. Ogawa and Ishiguro [4] solved the two-
dimensional Navier-Stokes equations in body-fitted coordinates using the
finite difference method to describe the incompressible flow of blood in the
human ventricle and the dynamic stall process on an oscillating airfoil.
Later, Yang et al. [5] followed by Raithby et al. [6] studied heat transfer
and fluid flow in complex geometries by applying non-orthogonal curvilinear
coordinates. More recently, Cannata et al. [7] simulated wave propagation
using an integral formulation of the contravariant Navier-Stokes equations in
an attempt to reproduce the complex geometries found in coastal regions.
Jarrah and Rizwan-Uddin [8] extended the applicability of the nodal
integration method scheme by deriving it in two-dimensional curvilinear
coordinates to solve the convection-diffusion equation discretised by
quadrilateral elements.

The implementation of the aforementioned coordinate transformation methods
requires the governing equations to be written in their corresponding
contravariant components form [9, 10, 11, 12, 13]. However, according to an
extensive review by Thompson [14] summarising the applications of grid
generation techniques in fluid dynamics problems, the governing equations in
the transformed domain can present additional source terms due to the
differentiation of metric coefficients remaining in the gradient terms [14].
Nevertheless, such additional terms can be handled using conservative
discretisation methods, such as the finite volume method (FVM), commonly used
for solving fluid dynamics problems due to their inherent conservative
characteristic and ability to ensure local and global flow conservation in the
discrete formulation. In this case, the transformed equations written in
contravariant components should be further expressed in their conversation
form to properly apply the FVM. This is particularly important when the FVM is
used to discretise the equations as such problem is non encountered if other
discretisation methods (_e.g._ finite differences and finite elements) are
used, as such techniques do not require the strong conservation form of the
equations, which is inherently their limitation. Readers may refer to [2, 13,
14] for a detailed description of coordinate transformation techniques and
their advantages in solving engineering applications.

To the authors’ knowledge, the derivation of the governing equations and
solution of the fully coupled thermo-elastohydrodynamic lubrication problem
using coordinate transformation methods has rarely been exploited. Han and
Paranjpe [15] proposed a finite volume model to study the thermohydrodynamic
performance of journal bearings using an energy equation formulated in a body-
fitted coordinate system, yet unrealistic thermal boundary conditions were
adopted for the bushing-oil interface, with maximum generated pressures not
exceeding 100 MPa due to the conformality of the geometry. Moreover, Roberts
et al. [16] implemented a finite-element-based approach for modelling
multiphase lubrication problems using curvilinear shell elements; however,
their analysis neglected any thermal and cavitation-related effects that may
incur changes to the behaviour of the confined fluid.

The present contribution provides a new integrated formulation to solve point
contact TEHL problems on a body-fitted normalised curvilinear grid. The
transformation of the governing equations from physical Cartesian coordinates
to non-orthogonal normalised curvilinear coordinates is achieved
systematically using the Jacobian and the metric tensor of the coordinate
transformation in order to convert the derivatives into the newly defined
normalised curvilinear system in the computational domain. Moreover, the
contravariant velocity components are derived to avoid potential numerical
errors associated with the alignment of the surfaces with staggered velocity
components, as suggested by [17], and to obtain the conservation forms of the
transformed governing equations for the appropriated application of the finite
volume method. The coupled governing equations, including the solid
deformation equation, the $(p-\theta)$ generalised Reynolds equation and the
fluid energy equation, are solved iteratively within an efficient and fully
conservative finite volume based fluid-structure interaction (FSI) framework.
The conservative aspect of the proposed framework is advantageous because it
allows the more straightforward incorporation of the mass-conserving
$(p-\theta)$ cavitation model in the discretisation process. Numerical
benchmarks and experimental thermal data measured using infrared spectroscopy
are used to validate the developed framework. Results show that the developed
framework successfully predicts lubrication performance in terms of contact
temperature, generated fluid pressures, traction coefficients and film
thickness profiles. Although the present work focuses on counterformal
contacts, the proposed method is applicable and can be easily adapted to solve
any lubricated contact problem in the presence of complex geometries.

The blueprint of the current study is as follows. The governing equations of
the fluid flow problem, including the $(p-\theta)$ generalised Reynolds
equation, the fluid energy equation, and the fluid velocity equations, are
derived in Section 2 in the Cartesian and normalised curvilinear coordinates
in their respective strong conservation forms. Section 3 illustrates how the
transformed differential equations are solved using the finite volume method,
while Section 3.3 discusses the architecture of the FSI framework employed to
solve the derived equations numerically. The validity of the developed
framework is assessed in Section 4 against numerical and experimental results
before outlining in Section 5 the main findings and potential applications of
the present contribution.

##  2 Mathematical Formulations

This section presents a unified derivation of the TEHL governing equations in
a normalised curvilinear coordinate system that allows the solution of the
full TEHL problem in a single structured computational domain while preserving
the conservative nature of the equations. A concise overview of the
transformation relationships between Cartesian and general curvilinear
coordinates is shown in Appendix A. Furthermore, the general strong
conservation form of the transport equation for an arbitrary scalar physical
quantity written in these coordinate systems is summarised in Appendix B.

The proposed TEHL framework is developed based on the geometric configuration
illustrated in Fig. 1A, which is typical of a circular elastohydrodynamic
(EHD) contact. The framework implicitly considers the dependency of the
lubricant rheological properties (_e.g._ dynamic viscosity, density, thermal
conductivity and specific heat capacity) on temperature, pressure and shear-
rate effects. However, and for the sake of brevity, the rheological models
adopted to simulate the isothermal and thermal cases discussed in this study
will only be reported where relevant. For the derivation of the TEHL governing
equations in normalised curvilinear coordinates, the following coordinate
transformation from the physical Cartesian space $x^{i}\left(x,y,z\right)$ to
the mapped normalised computational space $X^{i}\left(X,Y,Z\right)$ was
adopted,

$\displaystyle x$ $\displaystyle=\left(L_{x}\right)X,\quad\quad\quad$
$\displaystyle 0\leq x\leq L_{x}\quad$ $\displaystyle\text{and}\quad 0\leq
X\leq 1,$ (1) $\displaystyle y$
$\displaystyle=\left(L_{y}\right)Y=\left(\frac{L_{x}}{r_{xy}}\right)Y,\quad\quad\quad$
$\displaystyle 0\leq y\leq L_{y}\quad$ $\displaystyle\text{and}\quad 0\leq
Y\leq 1,$ $\displaystyle z$ $\displaystyle=hZ+z_{1}=\left(\varepsilon
L_{x}\right)\left(HZ+Z_{1}\right),\quad\quad\quad$ $\displaystyle z_{1}\leq
z\leq\left(h+z_{1}\right)\quad$ $\displaystyle\text{and}\quad 0\leq Z\leq 1,$

where $L_{x}$ and $L_{y}$ are the characteristic lengths of the lubricated
interface in the $x$\- and $y$-directions, respectively, $h(x,y)$ is the
lubricant film thickness, and $z_{1}$ is the height of the lower surface
($z_{1}=0$ is assumed in the current study). Readers may refer to Appendix G
for a summary of the parameters and variables employed throughout this work.

Figure 1: Characteristics of the geometry considered by the proposed
framework. (A) Depicts a deep groove ball bearing, where a point contact is
formed between the ball and the race, while (B) provides a magnified schematic
of the contact illustrating the confined lubricant flowing in between the
deformed surfaces and projected onto a curvilinear grid system ($\zeta$,
$\eta$, $\gamma$). (C) Presents the transformed normalised computational
domain (_X_ , _Y_ , _Z_) representing the lubricated contact, and (D)
elucidates the properties of the control volume (CV) with respect to the
central node of the CV "P", the surrounding nodal positions and the CV faces
bounding "P".

###  2.1 Generalised Reynolds Equation with Elrod-Adams $(p-\theta)$
Cavitation Model

The generalised Reynolds equation (GRE) is coined as being a revamped form of
the conventional Reynolds equation for considering the variation of the
lubricant properties across the film thickness [18]. Due to the impact of
fluid-film rupture on the lubricant flow within the contact, the cavitation
phenomenon is accounted for in this work by considering the Elrod–Adams
$(p-\theta)$ cavitation model [19, 20] which automatically enforces the
complementary Jakobsson, Floberg and Olsson (JFO) boundary conditions for
cavitation [21, 22, 23]. This entails modifying the GRE by introducing a flow
parameter $\theta$ to account for the biphasic mixture within the cavitation
zones [24, 25]. Accordingly, the $(p-\theta)$ generalised Reynolds equation
can be written in the strong conservation vector form in Cartesian coordinates
as

$\displaystyle\underbrace{\boldsymbol{\nabla}\cdot\left({\varepsilon}\boldsymbol{\nabla}p\right)}_{\text{P}oiseuille\;Term}=\underbrace{\boldsymbol{\nabla}\cdot\left[\theta\left({\rho^{*}_{e}}\boldsymbol{u_{m}}+{\rho^{*}_{1}}\boldsymbol{u_{1}}\right)\right]}_{\text{C}ouette\;Term}\;+\underbrace{\frac{\partial\left(\theta\rho_{e}\right)}{\partial
t}}_{\text{T}ransient\;Term},$ (2)

where the coefficients which incorporate the variations of the lubricant
properties across the film thickness are defined as

$\displaystyle\varepsilon=\frac{\eta_{e}}{\eta^{{}^{\prime}}_{e}}\rho^{{}^{\prime}}-\rho^{{}^{\prime\prime}},\quad\quad\rho^{*}_{e}=2\eta_{e}\rho^{{}^{\prime}},\quad\quad\rho^{*}_{1}=\rho_{e}-\rho^{*}_{e},\quad\quad\rho_{e}=\bigintssss_{z_{1}}^{z_{2}}\rho\;dz.$

In the above equations,
$\boldsymbol{\nabla}(\cdot)=\left[\dfrac{\partial}{\partial
x},\dfrac{\partial}{\partial y}\right]^{T}$ is the gradient operator in
Cartesian coordinates, while $\boldsymbol{u_{m}}=\left[u_{m},v_{m}\right]^{T}$
and $\boldsymbol{u_{1}}=\left[u_{1},v_{1}\right]^{T}$ are the entrainment and
lower surface velocity vectors in the $x$\- and $y$-directions, respectively.
The integral terms $\eta_{e}$, $\eta^{{}^{\prime}}_{e}$, $\rho^{{}^{\prime}}$,
$\rho^{{}^{\prime\prime}}$ and $\rho^{*}_{e}$ are defined in Appendix D.

The JFO boundary conditions for cavitation are preferable to other cavitation-
related formulations as they impose mass conservation throughout the
lubricated domain, including at the ruptured and reformation boundaries.
According to this formulation, the magnitude of the fractional film content
($\theta$) is dictated by the following complementarity condition [26]

$\displaystyle(p-p_{cav})(1-{\theta})=0\quad{\rightarrow}\quad\left\\{\begin{array}[]{lr}p>p_{cav}\quad{\rightarrow}\quad{\theta}=1&\text{pressured
zones,}\\\ p=p_{cav}\quad{\rightarrow}\quad
0\;{\leq}\;{\theta}<1&\text{cavitated zones.}\end{array}\right.$ (5)

Considering the coordinate transformation of Eq. 1 in combination with the
transformation relationship for spatial derivatives of Eq. 35, along with the
normalisation parameters of the dependent variables summarised in Appendix G,
the dimensionless $(p-\theta)$ generalised Reynolds equation can be expressed
in the strong conservation vector form in normalised coordinates as

$\displaystyle{\boldsymbol{\nabla_{\scriptscriptstyle
X}}\cdot\bigl{(}\boldsymbol{\Gamma_{D}}\boldsymbol{\nabla_{\scriptscriptstyle
X}}\bar{p}\bigr{)}}=\boldsymbol{\nabla_{\scriptscriptstyle
X}}\cdot\bigl{[}\theta\bigl{(}\boldsymbol{\Gamma_{CM}}\boldsymbol{U_{m}}+\boldsymbol{\Gamma_{C1}}\boldsymbol{U_{1}}\bigr{)}\bigr{]}+{\frac{\partial\left(\bar{\rho}_{e}\theta
H\right)}{\partial\tau}},$ (6)

where

$\displaystyle\boldsymbol{\Gamma_{D}}$
$\displaystyle=\bar{\varepsilon}\begin{bmatrix}1&0\\\
0&r^{2}_{xy}\end{bmatrix},$ $\displaystyle\boldsymbol{\Gamma_{CM}}$
$\displaystyle=\bar{\rho^{*}_{e}}\begin{bmatrix}H&0\\\
0&Hr_{xy}\end{bmatrix},$ $\displaystyle\boldsymbol{\Gamma_{C1}}$
$\displaystyle=\bar{\rho^{*}_{1}}\begin{bmatrix}H&0\\\
0&Hr_{xy}\end{bmatrix}.$

Similarly, in the above equations, $\boldsymbol{\nabla_{\scriptscriptstyle
X}}(\cdot)=\left[\dfrac{\partial}{\partial X},\dfrac{\partial}{\partial
Y}\right]^{T}$ is the gradient operator in normalised coordinates, and
$\boldsymbol{U_{m}}=\left[U_{m},V_{m}\right]^{T}$ and
$\boldsymbol{U_{1}}=\left[U_{1},V_{1}\right]^{T}$ are the dimensionless
entrainment and lower surface velocity vectors in the $X$\- and
$Y$-directions, respectively. All relevant dimensionless parameters are found
in Appendix D and Appendix G.

###  2.2 Film Thickness Equation

The geometry of the lubricant gap in point contact EHL can be expressed in
Cartesian and normalised coordinates as [27]

$\displaystyle h(x,y)$
$\displaystyle=h_{0}+\frac{x^{2}}{2R_{x}}+\frac{y^{2}}{2R_{y}}+{\delta}{(x,y)},$
(7a) $\displaystyle H(X,Y)$ $\displaystyle=H_{0}+{A_{\scriptscriptstyle
X}X^{2}+A_{\scriptscriptstyle Y}Y^{2}}+\bar{\delta}{(X,Y)},$ (7b)

where $h_{0}$ is the rigid separation between the contacting bodies, $R_{x}$
and $R_{y}$ are the equivalent radii curvature of the contacting solids at the
centre of contact in the $x$\- and $y$-directions, respectively, and $\delta$
is the pressure-induced surface normal displacement, which is calculated by
the following Boussinesq convolution integral for a 3-D elastic half-space

$\displaystyle{\delta}(x,y)$ $\displaystyle=\left(\frac{2}{\pi
E^{\prime}}\right)\bigintsss\\!\\!\\!\\!\\!\bigintsss_{\Omega}\frac{{p}(x^{\prime},y^{\prime})}{\sqrt{\left(x-x^{\prime}\right)^{2}+\left(y-y^{\prime}\right)^{2}}}\,dx^{\prime}dy^{\prime},$
(8a) $\displaystyle\bar{\delta}(X,Y)$
$\displaystyle=\left(\frac{2L_{x}p_{0}}{\pi
E^{\prime}h_{0}r_{xy}}\right)\bigintsss\\!\\!\\!\\!\\!\bigintsss_{\bar{\Omega}}\frac{\bar{p}(X^{\prime},Y^{\prime})}{\sqrt{\left(X-X^{\prime}\right)^{2}+{r_{xy}^{-2}}\left(Y-Y^{\prime}\right)^{2}}}\,dX^{\prime}dY^{\prime},$
(8b)

with
$E^{\prime}=2\left[(1-{\nu_{1}}^{2})/E_{1})+(1-{\nu_{2}}^{2})/E_{2}\right]^{-1}$
being the effective elastic modulus.

###  2.3 Fluid Energy Equation

The general energy equation for fluid film lubrication with compressible
viscous fluids and internal heat transfer by diffusion governed by Fourier’s
law in an isotropic medium can be written in terms of temperature in the
strong conservation vector form in Cartesian coordinates as [17, 28]

$\displaystyle\begin{split}\frac{\partial\left(\rho c_{p}T\right)}{\partial
t}+\boldsymbol{\nabla}\cdot\left(\rho
c_{p}\boldsymbol{\upsilon}T\right)&=\boldsymbol{\nabla}\cdot\left(k\boldsymbol{\nabla}T\right)+Q_{p}+Q_{cp}+Q_{\Phi}+\dot{q}_{v},\\\
\\\ Q_{p}=\underbrace{{{\beta}T\left(\frac{\partial p}{\partial
t}+u\frac{\partial p}{\partial x}+v\frac{\partial p}{\partial
y}\right)}}_{\begin{subarray}{c}\text{Compressive}\\\
\text{Heating/Cooling}\end{subarray}},\quad\quad
Q_{cp}&=\underbrace{{\rho}T\left(\frac{\partial c_{p}}{\partial
t}+\boldsymbol{\upsilon}\cdot\boldsymbol{\nabla}c_{p}\right)}_{\begin{subarray}{c}\text{Enthalpic}\\\
\text{Heating/Cooling}\end{subarray}},\quad\quad
Q_{\Phi}=\underbrace{{{\eta}\left[\left(\frac{\partial u}{\partial
z}\right)^{2}+\left(\frac{\partial v}{\partial
z}\right)^{2}\right]}}_{\begin{subarray}{c}\text{Shear}\\\
\text{Heating}\end{subarray}},\end{split}$ (9)

where $\rho$, $c_{p}$, $k$, $\beta$ and $\eta$ are the lubricant density,
thermal heat capacity, thermal conductivity, thermal compressibility and
dynamic viscosity, respectively. Moreover,
$\boldsymbol{\upsilon}=\left[u,v,w\right]^{T}$ is the fluid velocity vector,
and $\dot{q}_{v}$ is the rate of heat source or sink within the fluid volume
due to an arbitrary source term, such as heat released as a result of
asperity-asperity contact, which is neglected in the current work.

Using the coordinate transformation relationships summarised in Appendix A in
combination with the strong conservation form of a scalar transport equation
for general non-orthogonal curvilinear coordinate systems shown in Appendix B,
the dimensionless energy equation can be written in the strong conservation
vector form in normalised coordinates as

$\displaystyle\text{Pe}\left[\frac{\partial\left(H\bar{\rho}\bar{c}_{p}\bar{T}\right)}{\partial\tau}\right]+\text{Pe}\left[\boldsymbol{\nabla_{\scriptscriptstyle
X}}\cdot\left(\bar{\rho}\bar{c}_{p}\boldsymbol{\tilde{V}}\bar{T}\right)\right]=\boldsymbol{\nabla_{\scriptscriptstyle
X}}\cdot\left(\left[\boldsymbol{\bar{\tilde{\kappa}}}\right]\boldsymbol{\nabla_{\scriptscriptstyle
X}}\bar{T}\right)+\bar{\tilde{Q}}_{p}+\bar{\tilde{Q}}_{cp}+\bar{\tilde{Q}}_{\Phi}+\bar{\tilde{\dot{q}}}_{V},$
(10) $\displaystyle\begin{aligned}
\bar{\tilde{Q}}_{p}&=\text{Br}^{*}\bar{\beta}\bar{T}\left(H\frac{\partial\bar{p}}{\partial\tau}+\tilde{V}^{X}\frac{\partial\bar{p}}{\partial
X}+\tilde{V}^{Y}\frac{\partial\bar{p}}{\partial
Y}\right),\quad\quad\bar{\tilde{Q}}_{cp}=\text{Pe}\bar{\rho}\bar{T}\left(H\frac{\partial\bar{c}_{p}}{\partial\tau}+\boldsymbol{\tilde{V}}\cdot\boldsymbol{\nabla_{\scriptscriptstyle
X}}\bar{c}_{p}\right),\\\\[4.30554pt]
\bar{\tilde{Q}}_{\Phi}&=\frac{\text{Br}\bar{\eta}}{H^{3}}\left[\left(\frac{\partial\tilde{V^{X}}}{\partial
Z}\right)^{2}+\left(\frac{1}{r_{xy}}\frac{\tilde{V}^{Y}}{\partial
Z}\right)^{2}\right],\quad\quad\quad\quad\bar{\tilde{\dot{q}}}_{V}=H\bar{{\dot{q}}}_{V},\end{aligned}$
$\displaystyle\begin{aligned}
\left[\boldsymbol{\bar{\tilde{\kappa}}}\right]&\approx\bar{k}\begin{pmatrix}\varepsilon^{2}H&0&-\varepsilon^{2}\left(\dfrac{\partial
H}{\partial X}Z\right)\\\ 0&\left(\varepsilon
r_{xy}\right)^{2}H&-\left(\varepsilon r_{xy}\right)^{2}\left(\dfrac{\partial
H}{\partial Y}Z\right)\\\ -\varepsilon^{2}\left(\dfrac{\partial H}{\partial
X}Z\right)&-\left(\varepsilon r_{xy}\right)^{2}\left(\dfrac{\partial
H}{\partial
Y}Z\right)&\left(\dfrac{1}{H}\right)+\varepsilon^{2}\left[\left(\dfrac{\partial
H}{\partial X}\right)^{2}+\left(r_{xy}\dfrac{\partial H}{\partial
Y}\right)^{2}\right]\dfrac{Z^{2}}{H}\end{pmatrix}\\\
&\approx\bar{k}\begin{pmatrix}\varepsilon^{2}H&0&0\\\ 0&\left(\varepsilon
r_{xy}\right)^{2}H&0\\\
0&0&\left(\dfrac{1}{H}\right)+\varepsilon^{2}\left[\left(\dfrac{\partial
H}{\partial X}\right)^{2}+\left(r_{xy}\dfrac{\partial H}{\partial
Y}\right)^{2}\right]\dfrac{Z^{2}}{H}\end{pmatrix},\end{aligned}$

where
$\boldsymbol{\tilde{V}}=\left[\tilde{V}^{X},\tilde{V}^{Y},\tilde{V}^{Z}\right]^{T}$
is the dimensionless fluid velocity contravariant vector in the normalised
coordinate system, and $\left[\boldsymbol{\bar{\tilde{\kappa}}}\right]$ is the
dimensionless diffusion matrix. All other pertinent dimensionless parameters
are found in Appendix D and Appendix G. According to the lubrication theory,
the magnitude of the lubricant film thickness is much smaller than the other
dimensions of the contact. In this case, both the fluid pressure gradient
across the film thickness and the terms of the governing equations
proportional to the scaling factor $\varepsilon=\left(h_{0}/L_{x,y}\right)$
are assumed negligible, i.e., $\partial p/\partial z\approx 0$ and
$\boldsymbol{\mathcal{O}}(\varepsilon^{n})<<1$, hence why the non-diagonal
elements of the diffusion matrix have been omitted. However, those neglected
terms can play an important role at contact boundaries with pronounced
curvatures or in regions where the scale factor $\varepsilon$ is not tiny
(_i.e._ at the inlet and outlet regions). Moreover, the heat convection across
the film thickness and the heat diffusion along the contact (_i.e._ , the
first two elements of the main diagonal of
$\left[\boldsymbol{\bar{\tilde{\kappa}}}\right]$) are usually neglected in
lubrication problems but were kept for convenience to maintain the generality
of the proposed simulation framework.

###  2.4 Fluid Velocity Components

The Cartesian velocity components
$\boldsymbol{\upsilon}=\left[u,v,w\right]^{T}$ of the lubricant flow on the
contact interface obtained in the derivation of the generalised Reynolds
equation from the Navier-Stokes equation can be written as

$\displaystyle u(x,y,z)$ $\displaystyle=u_{1}+\frac{\partial p}{\partial
x}I_{p}+u_{s}I_{s},$ (11a) $\displaystyle v(x,y,z)$
$\displaystyle=v_{1}+\frac{\partial p}{\partial y}I_{p}+v_{s}I_{s},$ (11b)
$\displaystyle w(x,y,z)$
$\displaystyle=-\frac{1}{\rho}\frac{\partial}{\partial
x}\bigintssss_{z_{1}}^{z}\rho
u\;dz^{\prime}-\frac{1}{\rho}\frac{\partial}{\partial
y}\bigintssss_{z_{1}}^{z}\rho v\;dz^{\prime}$ (11c)

where $u_{1,2}$ and $v_{1,2}$ are the surface velocities in the $x$\- and
$y$-directions, respectively, while $u_{s}$ and $v_{s}$ are the sliding
velocities in the respective directions. The integral coefficients $I_{p}$ and
$I_{s}$ utilised in the above expressions are given in Appendix D. Moreover,
the Cartesian velocity component across the film thickness direction is
calculated by integrating the mass conservation equation shown in Appendix C.

Considering the coordinate transformation relationships summarised in Appendix
A and the normalisation parameters of the dependent variables summarised in
Appendix G, the dimensionless fluid velocity contravariant components
$\boldsymbol{\tilde{V}}=\left[\tilde{V}^{X},\tilde{V}^{Y},\tilde{V}^{Z}\right]^{T}$
in the normalised coordinate system can be expressed as

$\displaystyle\tilde{V}^{X}\left(X,Y,Z\right)$
$\displaystyle=H\underbrace{\left(U_{1}+\frac{\partial\bar{p}}{\partial
X}\bar{I}_{p}+U_{s}\bar{I}_{s}\right)}_{U}=HU,$ (12a)
$\displaystyle\tilde{V}^{Y}\left(X,Y,Z\right)$
$\displaystyle=r_{xy}H\underbrace{\left(V_{1}+r_{xy}\frac{\partial\bar{p}}{\partial
Y}\bar{I}_{p}+V_{s}\bar{I}_{s}\right)}_{V}=r_{xy}HV,$ (12b)
$\displaystyle\tilde{V}^{Z}\left(X,Y,Z\right)$
$\displaystyle=W-\left(\frac{\partial H}{\partial X}U+r_{xy}\frac{\partial
H}{\partial Y}V\right)Z,$ (12c)

where the dimensionless integral coefficients $\bar{I}_{p}$ and $\bar{I}_{s}$
are given in Appendix D.

##  3 Numerical Solution

This section demonstrates how the $(p-\theta)$ generalised Reynolds and fluid
energy equations are solved on the mapped normalised computational domain
using the finite volume method (FVM) [2, 17, 29, 30]. The FVM approach was
adopted to numerically solve the partial differential equations (PDEs)
governing the lubricant flow in TEHL contacts derived in Section 2 as it
automatically enforces the mass flow conservation [29], mainly when the strong
conservation form of the PDEs is applied [1], as carried out in this study.
Regardless of the chosen coordinate system, providing the strong conservation
form of the governing transport equations is crucial for finite volume
solutions because of the Gauss-divergence theorem that easily converts the
volume integral of divergence into a surface integral of fluxes. Furthermore,
the FVM has the relative advantage of being mathematically straightforward
and, in particular, is excellent for problems in which quantity conservation
is essential. This inherent conservative nature of the FVM is paramount in
predicting the cavitation regions in lubricated contacts via the application
of mass-conserving cavitation models. With regards to the Elrod-Adams
$(p-\theta)$ cavitation model adopted in this work, finite volume solutions
are favoured since they automatically enforce the complementary JFO boundary
conditions for cavitation via the upwind-based treatment of the convective
terms present in the $(p-\theta)$ generalised Reynolds equation [26, 29].
Moreover, applying FVM schemes on mapped-grids has proven effective in
considering the contribution of the convection and diffusion cross-derivative
terms resulting from non-orthogonal coordinate transformations [31, 32, 33].

###  3.1 Finite Volume Solution of the $(p-\theta)$ Generalised Reynolds
Equation

The finite volume discretisation of the dimensionless steady-state
$(p-\theta)$ generalised Reynolds equation begins by integrating Eq. 6 over
the control volume (CV) illustrated in Fig. 1D followed by applying the Gauss-
divergence theorem to convert the volume integrals of divergence into surface
integrals of fluxes. Accordingly:

$\displaystyle\oiiint_{\mathcal{V}}{\boldsymbol{\nabla_{\scriptscriptstyle
X}}\cdot\left(\boldsymbol{\Gamma_{D}}\boldsymbol{\nabla_{\scriptscriptstyle
X}}\bar{p}\right)}\,d\mathcal{V}$
$\displaystyle=\oiiint_{\mathcal{V}}\boldsymbol{\nabla_{\scriptscriptstyle
X}}\cdot\bigl{[}\theta\bigl{(}\boldsymbol{\Gamma_{CM}}\boldsymbol{U_{m}}+\boldsymbol{\Gamma_{C1}}\boldsymbol{U_{1}}\bigr{)}\bigr{]}\,d\mathcal{V},$
(13a)
$\displaystyle\oiint_{\partial\mathcal{S}}\left(\boldsymbol{\Gamma_{D}}\boldsymbol{\nabla_{\scriptscriptstyle
X}}\bar{p}\right)\cdot\vec{\boldsymbol{n}}\,d\mathcal{S}$
$\displaystyle=\oiint_{\partial\mathcal{S}}\bigl{[}\theta\bigl{(}\boldsymbol{\Gamma_{CM}}\boldsymbol{U_{m}}+\boldsymbol{\Gamma_{C1}}\boldsymbol{U_{1}}\bigr{)}\bigr{]}\cdot\vec{\boldsymbol{n}}\,d\mathcal{S},$
(13b)

where $\mathcal{V}$ is the CV volume of unit height (since the dimensionless
GRE is two-dimensional defined in the $XY$-plane) enclosing node $P$ as shown
in Fig. 1D, and $\vec{\boldsymbol{n}}$ is the unit vector normal to the CV
faces $\mathcal{S}$. The integrals over the CV faces can be approximated using
the midpoint rule such that the diffusive and convective fluxes are evaluated
at the centroid of each face [17, 26, 29]. This entails expressing the surface
integrals as a summation of the diffusive and convective fluxes travelling
through the CV faces, knowing that the surface unit vectors on opposite sides
of the CV have opposite signs. Therefore, the discrete form of Eq. 13b
formulated in terms of the central CV node ($P$) and the surrounding CV faces
($f$) becomes

$\displaystyle\sum_{f\sim
nb(P)}\left(\boldsymbol{\Gamma_{D}}\boldsymbol{\nabla_{\scriptscriptstyle
X}}\bar{p}\right)_{f}\cdot\boldsymbol{\mathcal{S}}_{f}=\sum_{f\sim
nb(P)}\bigl{[}\theta\bigl{(}\boldsymbol{\Gamma_{CM}}\boldsymbol{U_{m}}+\boldsymbol{\Gamma_{C1}}\boldsymbol{U_{1}}\bigr{)}\bigr{]}_{f}\cdot\boldsymbol{\mathcal{S}}_{f},$
(14)

where $\boldsymbol{\mathcal{S}}_{f}=({S_{f}}\vec{\boldsymbol{n}}_{f})$
represents the surface vectors normal to the CV faces of unit heights, which
include the orientation of the fluxes and the face lengths/widths they pass
through. The diffusive and convective fluxes at the CV faces were approximated
using the central (CDS) and upwind (UWD) differencing schemes, respectively,
based on the quantities of the neighbouring CV centroids.

The semi-system approach proposed initially by [34] is utilised in the present
work to improve the stability of the iterative numerical solution of the EHL
problem. This is carried out by expressing the film thickness in the Couette
term (R.H.S) of the GRE at a nodal position ($i,j$) as a function of the
neighbouring nodal pressures using the influence coefficients
$\left(D^{i,j}_{k,l}\right)$ [35, 36] given by the discretisation of the
Boussinesq integral and assuming a piecewise zero-order polynomial pressure
distribution, as shown in Eq. 15. Such an approach ensures the diagonal
dominance of the coefficient matrix in the solution of the GRE, ergo avoiding
an ill-conditioned matrix which stabilises the iterative solution of the
linear system of equations, particularly in thin film problems or under severe
loading conditions where fluid viscosity could rise to several orders of
magnitude.

$\displaystyle
H_{i,j}=H_{0}+\frac{X_{i}^{2}+Y_{j}^{2}}{2R^{\prime}}+\biggl{[}D^{i,j}_{k-1,l}\left(\bar{p}_{k-1,l}\right)+D^{i,j}_{k,l}\left(\bar{p}_{k,l}\right)+D^{i,j}_{k+1,l}\left(\bar{p}_{k+1,l}\right)\biggr{]}.$
(15)

Therefore, the linear system of equations derived from the finite volume
discretisation of the dimensionless $(p-\theta)$ generalised Reynolds equation
for semi-system EHL solution can be written as

$\displaystyle\boldsymbol{\alpha}\left(\bar{p}_{W}\right)+\boldsymbol{\beta}\left(\bar{p}_{P}\right)+\boldsymbol{\gamma}\left(\bar{p}_{E}\right)+\boldsymbol{\zeta}\left(\bar{p}_{N}\right)+\boldsymbol{\lambda}\left(\bar{p}_{S}\right)=\boldsymbol{\varphi},$
(16)

where $\boldsymbol{\alpha}$, $\boldsymbol{\beta}$, $\boldsymbol{\gamma}$,
$\boldsymbol{\zeta}$, $\boldsymbol{\lambda}$ and $\boldsymbol{\varphi}$ are
defined in Appendix E.

###  3.2 Finite Volume Solution of the Fluid Energy Equation

The finite volume discretisation of the dimensionless steady-state energy
equation for fluid film lubrication follows the same steps adopted for the
GRE, except that for the fluid energy equation, the discretisation is
conducted in a three-dimensional grid. Hence, by integrating Eq. 10 over the
control volume (CV) illustrated in Fig. 1D followed by applying the Gauss-
divergence theorem one has

$\displaystyle\text{Pe}\oiiint_{\mathcal{V}}\left[\boldsymbol{\nabla_{\scriptscriptstyle
X}}\cdot\left(\bar{\rho}\bar{c}_{p}\boldsymbol{\tilde{V}}\bar{T}\right)\right]\,d\mathcal{V}$
$\displaystyle=\oiiint_{\mathcal{V}}\boldsymbol{\nabla_{\scriptscriptstyle
X}}\cdot\left(\left[\boldsymbol{\bar{\tilde{\kappa}}}\right]^{T}\boldsymbol{\nabla_{\scriptscriptstyle
X}}\bar{T}\right)\,d\mathcal{V}+\oiiint_{\mathcal{V}}\bar{\tilde{Q}}_{T}\,d\mathcal{V},$
(17a)
$\displaystyle\text{Pe}\oiint_{\partial\mathcal{S}}\left(\bar{\rho}\bar{c}_{p}\boldsymbol{\tilde{V}}\bar{T}\right)\cdot\vec{\boldsymbol{n}}\,d\mathcal{S}$
$\displaystyle=\oiint_{\partial\mathcal{S}}\left(\left[\boldsymbol{\bar{\tilde{\kappa}}}\right]^{T}\boldsymbol{\nabla_{\scriptscriptstyle
X}}\bar{T}\right)\cdot\vec{\boldsymbol{n}}\,d\mathcal{S}+\oiiint_{\mathcal{V}}\bar{\tilde{Q}}_{T}\,d\mathcal{V}.$
(17b)

In addition to the midpoint rule for approximating the flux integrals at the
centroid of each CV face, the source term can be approximated as the product
between its magnitude at the CV centre and the CV volume. Hence [17]:

$\displaystyle\text{Pe}\sum_{f\sim
nb(P)}\left(\bar{\rho}\bar{c}_{p}\boldsymbol{\tilde{V}}\bar{T}\right)_{f}\cdot\boldsymbol{\mathcal{S}}_{f}=\sum_{f\sim
nb(P)}\left(\left[\boldsymbol{\bar{\tilde{\kappa}}}\right]^{T}\boldsymbol{\nabla}\bar{T}\right)_{f}\cdot\boldsymbol{\mathcal{S}}_{f}+\left(\bar{\tilde{Q}}_{T}\right)_{P}\cdot\left(\Delta\mathcal{V}\right)_{P},$
(18)

where $\Delta\mathcal{V}$ is the CV volume whereas
$\boldsymbol{\mathcal{S}}_{f}$ is the surface vectors normal to the CV faces
that include the orientation of the fluxes and the face areas they pass
through. The convective and diffusive fluxes were also approximated using the
upwind and central differencing schemes, respectively. Therefore, the system
of algebraic equations derived from the finite volume discretisation of the
dimensionless fluid energy equation can be written as

$\displaystyle\boldsymbol{\alpha}\left(\bar{T}_{W}\right)+\boldsymbol{\beta}\left(\bar{T}_{P}\right)+\boldsymbol{\gamma}\left(\bar{T}_{E}\right)+\boldsymbol{\zeta}\left(\bar{T}_{N}\right)+\boldsymbol{\lambda}\left(\bar{T}_{S}\right)+\boldsymbol{\phi}\left(\bar{T}_{T}\right)+\boldsymbol{\psi}\left(\bar{T}_{B}\right)=\boldsymbol{\varphi},$
(19)

where the above coefficients are defined in Appendix F.

###  3.3 Fluid-Structure Interaction Framework

The FSI framework developed by the authors in [37] for simulating two-
dimensional TEHL problems has been scaled up in this work to accommodate
three-dimensional TEHL solutions following the proposed mathematical and
numerical formulations described in the previous sections. The framework
begins by initialising the EHL solver to calculate the lubricant pressure and
film thickness while neglecting temperature changes. Afterwards, the results
of the EHL solver are fed into the TEHL solver to update the temperature
distributions in the lubricant film and solids. The overall convergence is
obtained when relative pressure and thermal solutions errors are minimised
below a certain convergence criterion. Although a brief description of the EHL
and TEHL solvers will be provided here, readers may refer to Section 4.1 from
[37] for a detailed description of the simulation framework.

The EHL solver calculates the lubricant pressure and film thickness through an
iterative process. Using an iterative Gauss-Seidel line relaxation method, the
lubricant pressure $\left(\bar{p}\right)$ and film fraction
$\left(\theta\right)$ are computed by solving the dimensionless $(p-\theta)$
generalised Reynolds equation (Eq. 16). Concurrently, lubricant properties,
such as viscosity $\left(\bar{\eta}\right)$ and density
$\left(\bar{\rho}\right)$, are updated due to pressure and shear-rate effects
using pre-defined rheological models. At the same time, the pressure-induced
surface deformation $\left(\bar{\delta}\right)$ and its influence on film
thickness $\left(H\right)$ are evaluated using Eq. 7, where the discrete
convolution and fast Fourier transform (DC-FFT) algorithm [35, 38] are adopted
to solve the Boussinesq integral efficiently. The EHL iterative process
terminates when the relative error between the computed pressures
$\left(\bar{p}\right)$ at successive iterations is minimised, and the
equilibrium of forces over the contact domain ${\Omega}$, given in the
dimensionless form in Eq. 20 for point contact problems, is achieved.

$\displaystyle\iint_{\Omega}\bar{p}(X,Y)\,dX\,dY=\frac{2\pi}{3}.$ (20)

Subsequently, the converged EHL results initialise the thermal loop,
simultaneously calculating the lubricant $\left(\bar{T}\right)$ and solids
temperatures $\left(\bar{T}_{s{1,2}}\right)$ iteratively and updating the
fluid rheological behaviour. The dimensionless fluid energy equation (Eq. 19)
is solved using a Gauss-Seidel line relaxation, whereas the solids
temperatures are calculated similarly but neglecting any source term. The
coupling between the thermal-fluid and thermal-solid solvers is established by
imposing heat flux continuity at the fluid-solid interfaces using the
conjugate heat transfer method [39]. The numerical convergence of the thermal
loop is achieved when the errors of the lubricant and solid temperatures, the
heat flux transmitted and the temperatures at the fluid-solid interfaces are
minimised below a specified tolerance between successive iterative steps.

Once the EHL and thermal solvers converge, a final convergence check on the
pressure and temperatures is performed to evaluate the global convergence of
the TEHL solution. The EHL solver is then recalled if global convergence is
not achieved, thus repeating the steps mentioned above until the global
convergence is reached. Point Gauss-Seidel Method with Aitken Acceleration
(PGMA) is integrated into the simulation framework to improve the convergence
of the iterative solvers. Aimed at coupling black-box solvers in fluid-
structure interaction problems [40], PGMA is based on the dynamic variation of
the relaxation parameters and has proved to be more effective in simulating
two-dimensional TEHL problems than Fixed Point Gauss-Seidel Method (PGMF),
where the latter employs constant relaxation parameters [37]. The PGMA method
updates the magnitude of the relaxation factors (${\omega_{\psi}^{n}}$) for an
arbitrary field variable $\psi$ using residues of the current
$\left(\vec{\boldsymbol{\mathfrak{R}}}_{\psi}^{n}\right)$ and previous
$\left({\vec{\boldsymbol{\mathfrak{R}}}_{\psi}^{n-1}}\right)$ iterations as
follows

$\displaystyle\begin{split}{\omega_{\psi}^{n}}&={\omega_{\psi}^{n-1}}\left[\frac{\left({\vec{{\boldsymbol{\mathfrak{R}}}}_{\psi}^{n-1}}\right)^{{T}}\cdot\left({\vec{{\boldsymbol{\mathfrak{R}}}}_{\psi}^{n-1}}-{\vec{{\boldsymbol{\mathfrak{R}}}}_{\psi}^{n}}\right)}{{\left({\vec{{\boldsymbol{\mathfrak{R}}}}_{\psi}^{n}}-{\vec{{\boldsymbol{\mathfrak{R}}}}_{\psi}^{n-1}}\right)}^{{T}}\cdot\left({\vec{{\boldsymbol{\mathfrak{R}}}}_{\psi}^{n}}-{\vec{{\boldsymbol{\mathfrak{R}}}}_{\psi}^{n-1}}\right)}\right],\\\\[4.30554pt]
{\omega_{\psi}^{n}}\;&{\xleftarrow{\;}\text{min}\left[\text{max}\left({\omega_{\psi}^{n}},\;{\omega_{min}}\right),\;{\omega_{max}}\right]},\end{split}$
(21a) where $\omega_{min}$ and $\omega_{max}$ are usually set as 0.001 and
1.0, respectively, and are imposed to limit the magnitude of
${\omega_{\psi}^{n}}$. The residue and the update of a field variable in a
given iteration $n$ are calculated as
$\displaystyle\begin{split}{\vec{\boldsymbol{\mathfrak{R}}}^{n}_{\psi}}&={\vec{\boldsymbol{\psi}}^{n}}-{\vec{\boldsymbol{\psi}}^{n-1}},\\\\[4.30554pt]
{\vec{\boldsymbol{\psi}}^{n}}&={\vec{\boldsymbol{\psi}}^{n-1}}+\omega_{\psi}^{n}\left({\vec{\boldsymbol{\mathfrak{R}}}^{n}_{\psi}}\right).\end{split}$
(21b)

Readers may refer to [37] for a detailed description of PGMA implementation
and a comprehensive analysis of the improved convergence rates that the PGMA
offers compared to the PGMF partitioned technique.

##  4 Results and Discussion

The accuracy of the developed TEHL simulation framework is assessed using
numerical results representative of thermal point-loaded contacts found in the
literature. Furthermore, temperature measurements obtained using infrared
spectroscopy and experimental traction data are analysed to corroborate the
proposed framework’s correctness further. The rheological models used to
characterise the lubricant behaviour due to changes in pressure, temperature
and shear-rate are specified according to each relevant study considered.

###  4.1 Numerical Benchmark

The first numerical benchmark is based on a series of investigations carried
out by Kaneta et al. [41] aimed at studying the influence of the slide-roll
ratio (SRR) on the lubrication performance of a $\mathrm{Si_{3}N_{4}}$-steel
contact. The oil viscosity change due to pressure and temperature effects is
described using the thermal Roelands model (Eq. 22a) and the Eyring model (Eq.
22b) is adopted to account for the lubricant non-Newtonian behaviour, while
the oil density is calculated using the Dowson and Higginson relationship (Eq.
22c). All the parameters and operating conditions required for the TEHL
simulation can be found in [41].

$\displaystyle\eta$
$\displaystyle=\eta_{0}\;\exp\Biggl{\\{}(\ln\eta_{0}+9.67)\times\bigg{[}-1+\left(1+5.1\times
10^{-9}p\right)^{\dfrac{\alpha}{5.1\times
10^{-9}\left(\ln\eta_{0}+9.67\right)}}\left(\frac{T-138}{T_{0}-138}\right)^{-\dfrac{\beta(T_{0}-138)}{\left(\ln\eta_{0}+9.67\right)}}\bigg{]}\Biggl{\\}},$
(22a) $\displaystyle\eta^{*}$
$\displaystyle=\frac{\eta\left(\tau_{e}/\tau_{0}\right)}{\sinh{\left(\tau_{e}/\tau_{0}\right)}},$
(22b) $\displaystyle\rho$ $\displaystyle=\rho_{0}\left(1+\frac{0.6\times
10^{-9}p}{1+1.7\times
10^{-9}p}\right)\left[1-0.00065\left(T-T_{0}\right)\right].$ (22c)

Figure 2: Comparison of the mid-plane ($Y=0$) dimensionless pressure and film
thickness profiles and lubricant film, steel (lower solid) and
$\mathrm{Si_{3}N_{4}}$ (upper solid) surface temperatures along the
entertainment $X$-direction calculated with the proposed TEHL framework
(lines) against their counterparts (dots) in Kaneta et al. [41]. The cases are
representative of $p_{h}=1\;\text{GPa}$, $R^{{}^{\prime}}=0.02\;\text{m}$,
$T_{0}=303\;\text{K}$, $u_{e}=1.32\;\text{m/s}$ and slide-roll-ratios of (A)
SRR = 0.1, (B) SRR = 0.5, (C) SRR = 1.0, (D) SRR = 1.5 and (E) SRR = 1.5.

It can be deduced from Fig. 2 that a good agreement is achieved for pressure,
film thickness and temperature profiles. As verified in the figure, the film
temperature increases with SRR due to the higher heat generated by shear near
the mid-plane of the film, thus explaining the higher lubricant temperature
compared to the surface temperatures regardless of the SRR. Because of the
temperature rise in the contact, fluid viscosity tends to decrease, resulting
in lower film thickness for the highest SRRs. Furthermore, the difference in
the surface temperatures can be explained by the different thermal
conductivities of the surfaces, wherein the higher temperatures are found for
the $\mathrm{Si_{3}N_{4}}$ surface ($k=23\;\text{W/mK}$) as it carries less
heat away compared to steel ($k=46\;\text{W/mK}$). This is the case for low
SRR, while the situation differs for higher SRRs where convection effects
dominate.

The second numerical benchmark employed to validate the developed TEHL
simulation framework is based on another study conducted by Kaneta et al. [42]
aimed at assessing the thermal behaviour of lubricated contacts when solids of
dissimilar thermal properties are considered. The point contact configurations
are that of steel-SiC, steel-$\mathrm{Al_{2}O_{3}}$, and
steel-$\mathrm{Si_{3}N_{4}}$ contacts where each solid pair is separated by a
continuously variable transmission (CVT) fluid and loaded with a maximum
Hertizan pressure $p_{h}=0.75\,\text{GPa}$. The input parameters required for
the TEHL simulations, such as the thermophysical properties of the solids and
lubricant, can be found in [42]. The viscosity- and density-pressure-
temperature dependencies of the lubricant are described via the following
Ohno’s relationships [43]

$\displaystyle\log_{10}\eta$
$\displaystyle=7-\frac{11.3\left(T-T_{VE}\right)\left(\frac{206}{T_{VE}}\right)}{35.9+\left(T-T_{VE}\right)\left(\frac{206}{T_{VE}}\right)},$
(23a) $\displaystyle\ln\rho$
$\displaystyle=\ln\rho_{0}+0.0012T\left(1-\frac{206}{T_{VE}}\right)+\frac{p}{12.8\times
10^{8}},$ (23b)

where $T_{VE}=206+209.4\ln\left(1+1.445\times 10^{-9}p\right)$, while the non-
Newtonian lubricant behaviour is described using the Eyring shear-thinning
model given in Eq. 22b.

Figure 3: Comparison of the mid-plane ($Y=0$) dimensionless surface
temperature profiles along the entertainment $X$-direction calculated with the
proposed TEHL framework (lines) against their counterparts (dots) in [42].
(A), (B) and (C) portray the surface temperatures when the SiC,
$\mathrm{Al_{2}O_{3}}$ and $\mathrm{Si_{3}N_{4}}$ surfaces are travelling
faster than the steel surface, respectively, while (D), (E) and (F) correspond
to the thermal results during which the steel surface is faster than the SiC,
$\mathrm{Al_{2}O_{3}}$ and $\mathrm{Si_{3}N_{4}}$ surfaces, respectively.
Dimensional thermal maps (in Kelvin) pertinent to the bounding surfaces and
mid-layer of the confined film are illustrated for each contact pair and speed
condition. The numerical investigations were carried out for
$p_{h}=0.75\,\text{GPa}$, $T_{0}=353\,\text{K}$, SRR = 1 and entrainment
speeds of $u_{e}|_{SiC-Steel}=4.4\;\text{m/s}$,
$u_{e}|_{\mathrm{Al_{2}O_{3}}-Steel}=4.1\;\text{m/s}$ and
$u_{e}|_{\mathrm{Si_{3}N_{4}}-Steel}=3.8\;\text{m/s}$.

The robustness of the developed TEHL simulation framework can be inferred from
the comparisons illustrated in Fig. 3. Inevitably, the predicted temperature
of the SiC surface is always higher than that of steel regardless of how fast
the solids move relative to each other due to the significant difference in
thermal conductivity of the two materials ($k_{SiC}=200\,\text{W/mK}$ while
$k_{steel}=21\,\text{W/mK}$). On the other hand, switching the surface
velocity directions has little influence on the magnitude of the maximum solid
temperatures in the case of steel-$\mathrm{Al_{2}O_{3}}$ contact, where the
thermal conductivity of $\mathrm{Al_{2}O_{3}}$
($k_{\mathrm{Al_{2}O_{3}}}=29\,\text{W/mK}$) is similar to that of steel;
however, higher temperatures are always found on the slower surface, which
carries less heat away from the contact, as depicted in Figs. 3B and 3E.
Although the thermal conductivity of steel and $\mathrm{Si_{3}N_{4}}$ are
almost identical ($k_{\mathrm{Si_{3}N_{4}}}=23\,\text{W/mK}$), yet it is
evident from Figs. 3C and 3F that a significant change in solid temperatures
can be observed depending on which surface is travelling at a higher speed. To
further corroborate the proposed TEHL framework, Fig. 4 portrays how the
dimensionless lubricant velocity profiles across the film thickness at
($X=Y=0$) calculated using Eqs.12a, 12b and 12c agree very well with those
evaluated in [42] under thermal conditions. Furthermore, to gain a more
quantitative understanding of those results, Fig. 5 illustrates the mid-plane
($Y=0$) temperature, shear-rate and lubricant velocity distributions for the
different pair of materials.

Figure 4: Dimensionless lubricant velocity profiles across the film thickness
at ($X=Y=0$) calculated with the proposed TEHL framework (lines) against their
counterparts (dots) in [41] for SRR = 1 and entrainment speeds of
$u_{e}|_{SiC-Steel}=4.4\;\text{m/s}$,
$u_{e}|_{\mathrm{Al_{2}O_{3}}-Steel}=4.1\;\text{m/s}$ and
$u_{e}|_{\mathrm{Si_{3}N_{4}}-Steel}=3.8\;\text{m/s}$. Figure 5: Mid-plane
($Y=0$) contour maps corresponding to the operating conditions in Figs. 3 and
4. (A) Lubricant temperature, (B) shear-rate and (C) lubricant velocity
distributions.

###  4.2 Experimental Benchmark

In addition to validating against numerical-based results, the use of
experimental data enhances the reliability of the proposed set of governing
equations and evaluates its accuracy in predicting performance of lubricated
contacts. Thermal and traction results measured in an experimental study
carried out by Reddyhoff et al. [44] are used in the current work. It is
important to note that the use of temperature data for corroboration of
numerical results is unconventional and that the majority of validation
exercises are often limited to friction and film thickness measurements. The
authors in [44] determined the average temperature of a ball and disk point
contact using infrared spectroscopy technique where emissivity radiated from
the oil, ball and disk is captured and isolated using band pass optical
filters and a chromium coating applied to the disk. The obtained radiation is
then converted to thermal data according to Planck’s Law and several
calibration tests at different temperature and film thickness [44]. The tests
were conducted on a steel ($k_{\mathrm{Steel}}=21\,W/m.K$), silicon nitride
($k_{\mathrm{Si_{3}N_{4}}}=30\,W/m.K$) and zirconia
($k_{\mathrm{Zirconia}}=3\,W/m.K$) ball loaded against a sapphire disk and
separated using the traction fluid Santotrac 50 and PAO 4. Operating
conditions such as applied load and entrainment speeds were varied during the
tests [44] depending on the material of the ball in order to ensure a constant
central film thickness of 100 nm required for calibration under the conditions
of 0.583 GPa average contact pressure and 40 oC inlet oil temperature for all
the tests. In addition to the solid temperatures, traction coefficients
measured independently using MTM2 rig are used in the current validation. To
accurately simulate the thermal cases, the Roelands and Eyring models (Eqs.
24a and 24b respectively) [45] are used to characterise the viscosity of
Santotrac 50 due to pressure, temperature and shear-thinning effects, while
the Yasutomi model coupled with the Carreau equation (Eqs. 24c and 24d,
respectively) [46] are used to model the PAO 4 viscous behaviour. The density
variation of both fluids are numerically describe using the Dowson and
Higginson model [27] while their thermal conductivity and heat capacity are
correlated using the models by Larsson and Anderson [47]. Readers may refer to
[45, 46, 47] for the input parameters of the corresponding relationships,
while the solids thermal and mechanical properties can be found in [44].

$\displaystyle\eta$
$\displaystyle=\eta_{0}\;\exp\bigg{(}(\ln\eta_{0}+9.67)\times\bigg{[}-1+\left(1+5.1\times
10^{-9}p\right)^{\dfrac{\alpha}{5.1\times
10^{-9}\left(\ln\eta_{0}+9.67\right)}}\left(\frac{T-138}{T_{0}-138}\right)^{-\dfrac{\beta(T_{0}-138)}{\left(\ln\eta_{0}+9.67\right)}}\bigg{]}\bigg{)},$
(24a) $\displaystyle\eta^{*}$
$\displaystyle=\eta\left(\tau_{e}/\tau_{0}\right)\bigg{/}\sinh{\left(\tau_{e}/\tau_{0}\right)},$
(24b) $\displaystyle\mu$
$\displaystyle=\mu_{g}\exp\left[\frac{-2.303C_{1}\left(T-\left[T_{g0}+A_{1}\ln\left(1+A_{2}p\right)\right]\right)\left(1+b_{1}p\right)^{b_{2}}}{C_{2}+\left(T-\left[T_{g0}+A_{1}\ln\left(1+A_{2}p\right)\right]\right)\left(1+b_{1}p\right)^{b_{2}}}\right],$
(24c) $\displaystyle\eta$
$\displaystyle=\mu\left[1+\left(\frac{\tau_{e}}{18\times
10^{6}}\right)^{2}\right]^{-6/13},$ (24d) $\displaystyle\rho$
$\displaystyle=\rho_{0}\left[1+\frac{0.6\times 10^{-9}p}{1+1.7\times
10^{-9}p}\right]\left[1-0.00065\left(T-T_{0}\right)\right],$ (24e)
$\displaystyle k$ $\displaystyle=k_{0}\left(1+\frac{c_{1}p}{1+c_{2}p}\right),$
(24f) $\displaystyle\rho c_{p}$ $\displaystyle=\left(\rho
c_{p}\right)_{0}\left[1+\beta(p)(T-T_{0})\left(1+\frac{k_{1}p}{1+k_{2}p}\right)\right].$
(24g)

Figure 6: Average surface temperatures and traction coefficients for a ball on
disk contact simulated using the finite volume framework (FVM) and validated
using experimental results measured via infrared microscopy and MTM2 test rig
for different ball materials and a sapphire disk; both solids bounding either
Santotrac 50 or PAO 4 oil. The solid-fluid pairs and their respective
entertainment speeds ($\mathrm{u_{e}}$) are: (a) steel-sapphire-Santotrac 50
($\mathrm{u_{e}=0.308\;m/s}$), (b) silicon nitride-sapphire-Santotrac 50
($\mathrm{u_{e}=0.293\;m/s}$), (c) zirconia-sapphire-Santotrac 50
($\mathrm{u_{e}=0.308\;m/s}$), (d) steel-sapphire-PAO 4
($\mathrm{u_{e}=0.929\;m/s}$), (e) silicon nitride-sapphire-PAO 4
($\mathrm{u_{e}=0.806\;m/s}$) and (f) zirconia-sapphire-PAO 4
($\mathrm{u_{e}=0.929\;m/s}$). All cases are carried out for SRR = [0, 0.5,
1.0, 1.5] and applied loads $\mathrm{W_{steel-sapphire}=20\;N}$,
$\mathrm{W_{zirconia-sapphire}=20\;N}$ and $\mathrm{W_{steel-
silicon\;nitride}=12\;N}$, regardless of the entraining oil.

It is evident from Fig. 6 that the current framework correctly predicts an
increase in solid temperature with SRR due to the increase in heat generation.
Vis-à-vis the contacts involving steel and silicon nitride with Santotrac 50
and PAO 4 (Figs. 6a, b, d and e), a considerable good match is obtained for
the solid temperatures between those obtained numerically using the framework
and their counterparts measured experimentally. However, the temperatures of
the steel ball in Fig. 6d and sapphire disk in Fig. 6e at SRR = 1.5 are found
to be numerically overestimated with a 10-11 % deviation from their
experimental results. This could be attributed to the film thickness and its
impact on the proportion of generated heat in the contact; unlike the tests
during which the central film thickness was kept constant at 100 nm, the
developed framework computes a 110 nm central thickness for both cases as it
is unfeasible to impose a constant thickness in the simulation due to the
influence of pressure, temperature and shear-rate effects on the film
geometry. On the other hand, it is obvious that significant discrepancies
between the measured and computed average solid temperatures are found for the
cases involving the zirconia ball, with either Santotrac 50 or PAO 4 (Figs, 6c
and f, respectively). This is not surprising as zirconia is characterised as a
highly effective thermal resistance material, thus possessing an ability to
hinder the flow of heat through its depth, which consequently results in high
temperatures at the outer layer of the solid where a significant proportion of
heat is accumulated. As a result of the near-infrared radiative properties of
zirconia, the utilised infrared camera in [44] was incapable of accurately
measuring the temperature distributions which led the authors to implement
Carslaw & Jaeger’s equation, as reported in [44].

Furthermore, traction coefficient data illustrated in Fig. 6 confirm the
validity of the developed numerical models in predicting the lubrication
performance. The only noticeable discrepancy between experiments and numerical
results is obtained under pure-rolling conditions in cases where Santotrac 50
is applied (Figs. 6a, b and c). This could be attributed to the limitations
which could occur when controlling the thickness at a prescribed level in
cases where surfaces travel at similar speeds, especially considering that a
good match in friction coefficients is obtained at the same speed conditions
but when PAO 4 is used (Figs. 6d, e and f). Moreover, the good match in
friction coefficients obtained with the zirconia ball shows that the
numerically predicted heat generation is accurately predicted, ergo
elucidating the aforementioned reasoning to identify what may have caused the
mismatch in measure and computed average temperatures.

Figure 7: Three-dimensional cross-sectional view depicting the mid-plane
distribution of the film temperature ($\mathrm{T}$), film dynamic viscosity
($\mathrm{\eta}$) and a two-dimensional illustration of the film content
($\mathrm{\theta}$) at slide-roll-ratio SRR = 1.5 for a silicon nitride ball
loaded against a sapphire disk (W = 20 N), where the two solids are separated
using (A) Santotrac 50 oil at entrainment speed $\mathrm{u_{e}=0.293\;m/s}$
and (B) PAO 4 oil at entrainment speed $\mathrm{u_{e}=0.806\;m/s}$. According
to the complementary JFO conditions enforced by the Elrod-Adams cavitation
algorithm, film rupture occurs when the fractional film content
($\mathrm{\theta}$) drops below unity, hence yielding in a biphasic mixture of
gas/vapour, as shown near the contact outlet.

To achieve a clear visualisation of the FVM thermal results presented in Figs.
6, and particularly those related to the silicon nitride and sapphire contact
case at SRR = 1.5, Fig. 7 illustrates three-dimensional cross-sectional maps
of the film temperature, dynamic viscosity as well as a two-dimensional
representation of the fluid film fraction dictated by the solution of the
$(p-\theta)$ generalised Reynolds equation, where Figs. 7A and 7B correspond
to Santotrac 50 and PAO4 oils, respectively. Regarding the temperatures of the
two films, it is evident that shear heating effects are more dominant in the
case of PAO4, hence yielding in slightly higher temperatures compared to
Santotrac 50. The direct influence of temperature on film viscosity ergo
yields much lower PAO4 viscosity values compared to Santotrac 50, the latter
being characterised by a more uniform viscosity distribution across the
contact length. However, it is important to note that different rheological
models were used to account for the behaviour of Santotrac 50 and PAO4 as
previously described, due to the limited availability of relevant parameters;
hence this can potentially play an important role in these comparisons. Film
fraction maps, which are rarely discussed in the field of TEHL simulations,
illustrate how film rupture commences near the outlet of the contact,
particularly where the divergent zone is found, where a higher liquid-
vapour/gases mixture can be noticed within the cavitation zone in the case of
PAO4 in Fig. 7B (12% liquid compared to 88% mixture of liquid and
gases/vapour).

##  5 Conclusion

We have presented a finite volume fluid-structure interaction (FSI) framework
for solving thermal-EHL problems mapped on non-orthogonal curvilinear grids
using the strong conservation form of the pertinent governing equations.
Numerical and experimental benchmarks demonstrated the accuracy of the
developed TEHL simulation framework in predicting the performance of
lubricated interfaces, where pressure, film thickness, friction forces and
fluid and solid temperature profiles have been determined. The finite volume
implementation facilitates the numerical discretisation of the governing
equations in the computational domain, including the convective and diffusive
cross-derivative terms of the fluid energy equation resulting from the non-
orthogonal grid mapping, while maintaining the virtue of conservation property
in the discrete level of finite volume methods, which is expensive to other
discretisation methods. Furthermore, the Elrod-Adams $(p-\theta)$ cavitation
model is combined with the generalised Reynolds equation to account for fluid
film cavitation and detect the formation of biphasic mixtures within the
contact. The effectiveness of the developed framework for coupling the TEHL
solvers was demonstrated by simulating lubricated counterformal contacts under
various operating conditions. The architecture of the framework and the
additional implementation of the semi-system approach and Point Gauss-Seidel
Method with Aitken Acceleration (PGMA) for the EHL solution revamped the
numerical convergence and significantly decreased the computational cost of
the simulations.

The current contribution may be regarded as an enhancement to traditional
Reynolds-based solvers due to the limitations encountered in considering any
curvature effects when modelling FSI problems in the presence of complex
geometries. Although commercial CFD codes may excel in tackling such problems
due to their powerful resources, yet they suffer from higher computational
overhead time and additional complexities in generating surface meshes. The
proposed method of mapping the strong conservation form of the TEHL governing
equations on a computational domain using Reynolds-based solvers overcomes
these limitations, which will be further tested by analysing the details of
fluid recirculation and temperature gradients at the contact inlet compared
with both CFD and modified Reynolds solvers using standard Cartesian grids.
This method holds promise and will be further investigated in applications
where curvature effects and cross-derivative terms are more profound than
counterformal contacts, such as journal bearing conformal systems, which the
authors will study in future contributions.

##  Appendix A Coordinate Transformation from Cartesian to Curvilinear Grid
System

The time-dependent coordinate transformation of the physical Cartesian space
$\chi^{3}_{t}\subset\mathbb{R}^{3}$ from the transformed generalised
curvilinear space $\varepsilon^{3}_{\tau}\subset\mathbb{R}^{3}$, where
$\mathbb{R}^{3}$ is the three-dimensional Euclidian space, can be expressed as
[12, 13]

$\displaystyle\boldsymbol{x\left(\xi,\tau\right)}:\varepsilon^{3}_{\tau}\rightarrow\chi^{3}_{t},\quad\quad\boldsymbol{\xi}=\left(\xi,\eta,\gamma,\tau\right)\in\varepsilon^{3}_{\tau},\quad\quad\boldsymbol{x}=\left(x,y,z,t\right)\in\chi^{3}_{t},$
$\displaystyle
x^{i}=x^{i}\left(\xi^{j},\tau\right)=x^{i}\left(\xi,\eta,\gamma,\tau\right),\quad\left(i,j=1,2,3\right),$
(25)

where $\boldsymbol{x\left(\xi,\tau\right)}$ is assumed to be smooth and
invertible with respect to $\xi^{j}$ and $\tau$. Therefore, the time-dependent
inverse transformation can be written as

$\displaystyle\boldsymbol{\xi\left(x,t\right)}:\chi^{3}_{t}\rightarrow\varepsilon^{3}_{\tau},$
$\displaystyle\xi^{i}=\xi^{i}\left(x^{j},t\right)=\xi^{i}\left(x,y,z,t\right),\quad\left(i,j=1,2,3\right).$
(26)

In the above equations, $x^{i}$ are the Cartesian coordinates of a point in
the physical space which at time $t$ has curvilinear coordinates $\xi^{j}$.
The time-dependent coordinate transformation (Eq. A) allows the solution of
transport equations in time-dependent moving grids in physical space.

The coordinate transformation $\boldsymbol{x\left(\xi,\tau\right)}$ is
characterised by the transformation tensor $A^{i}_{j}$, also known as the
covariant basis matrix, and the determinant known as Jacobian $J$ defined as

$\displaystyle A^{i}_{j}=\left(\frac{\partial
x^{i}}{\partial\xi^{j}}\right)\quad\Rightarrow\quad\left[\boldsymbol{A}\right]=\begin{pmatrix}x_{\xi}&x_{\eta}&x_{\gamma}\\\
y_{\xi}&y_{\eta}&y_{\gamma}\\\ z_{\xi}&z_{\eta}&z_{\gamma}\end{pmatrix},\quad
J=det\left(\frac{\partial
x^{i}}{\partial\xi^{j}}\right)=det\left[\boldsymbol{A}\right],$ (27)

where $x_{\xi}$ and $z_{\eta}$ are shorthand notations for the partial
derivatives $\dfrac{\partial x}{\partial\zeta}$ and $\dfrac{\partial
z}{\partial\eta}$, respectively. Similarly, the inverse coordinate
transformation $\boldsymbol{\xi\left(x,t\right)}$, known as the contravariant
basis matrix, is characterised by the inverse transformation tensor and
Jacobian as follows

$\displaystyle\bar{A}^{i}_{j}=\left(\frac{\partial\xi^{i}}{\partial
x^{j}}\right)\quad\Rightarrow\quad\left[\boldsymbol{\bar{A}}\right]=\begin{pmatrix}\xi_{x}&\xi_{y}&\xi_{z}\\\
\eta_{x}&\eta_{y}&\eta_{z}\\\
\gamma_{x}&\gamma_{y}&\gamma_{z}\end{pmatrix},\quad\frac{1}{J}=det\left(\frac{\partial\xi^{i}}{\partial
x^{j}}\right)=det\left[\boldsymbol{\bar{A}}\right].$ (28)

The transformation tensors are correlated by the identity
$A^{i}_{j}\bar{A}^{k}_{j}=\delta^{k}_{j}$, or in matrix array notation,
$\boldsymbol{\left[A\right]\left[\bar{A}\right]=\left[I\right]}$ $\Rightarrow$
$\boldsymbol{\left[\bar{A}\right]=\left[A\right]^{-1}}$. Therefore, the
relationship between the elements of $\boldsymbol{\left[A\right]}$ and
$\boldsymbol{\left[\bar{A}\right]}$ can be written as

$\displaystyle\frac{\partial\xi^{i}}{\partial x^{j}}$
$\displaystyle=\frac{1}{J}\left(\frac{\partial
x^{j+1}}{\partial\xi^{i+1}}\frac{\partial
x^{j+2}}{\partial\xi^{i+2}}-\frac{\partial
x^{j+1}}{\partial\xi^{i+2}}\frac{\partial x^{j+2}}{\partial\xi^{i+1}}\right),$
(29a) $\displaystyle\frac{\partial x^{i}}{\partial\xi^{j}}$
$\displaystyle=J\left(\frac{\partial\xi^{j+1}}{\partial
x^{i+1}}\frac{\partial\xi^{j+2}}{\partial
x^{i+2}}-\frac{\partial\xi^{j+1}}{\partial
x^{i+2}}\frac{\partial\xi^{j+2}}{\partial x^{i+1}}\right).$ (29b)

In the above equations, the superscript indices follow the cyclic permutation
order (1,2,3). With this convention, any index, say $l$, is identified with
$l\pm 3$; thus, for instance, $\dfrac{\partial
x^{4}}{\partial\xi^{5}}=\dfrac{\partial
x^{1}}{\partial\xi^{2}}=\dfrac{\partial x}{\partial\eta}$.

The covariant base vectors $\boldsymbol{g_{i}}$ which describe how a point in
space changes as the coordinates change can be expressed as

$\displaystyle\boldsymbol{g_{i}}=\left(\frac{\partial
x^{j}}{\partial\xi^{i}}\right)\boldsymbol{i_{j}}=A^{j}_{i}\boldsymbol{i_{j}}\Rightarrow\begin{bmatrix}\boldsymbol{g_{1}}\\\
\boldsymbol{g_{2}}\\\
\boldsymbol{g_{3}}\end{bmatrix}=\left[\boldsymbol{A}\right]^{T}\begin{bmatrix}\boldsymbol{i_{1}}\\\
\boldsymbol{i_{2}}\\\ \boldsymbol{i_{3}}\end{bmatrix},$ (30)

while the contravariant base vectors $\boldsymbol{g^{i}}$ that describe how
the coordinates change as the point in space changes can be written as

$\displaystyle\boldsymbol{g^{i}}=\left(\frac{\partial\xi^{i}}{\partial
x^{j}}\right)\boldsymbol{i_{j}}=A^{j}_{i}\boldsymbol{i_{j}}\Rightarrow\begin{bmatrix}\boldsymbol{g^{1}}\\\
\boldsymbol{g^{2}}\\\
\boldsymbol{g^{3}}\end{bmatrix}=\left[\boldsymbol{\bar{A}}\right]\begin{bmatrix}\boldsymbol{i_{1}}\\\
\boldsymbol{i_{2}}\\\ \boldsymbol{i_{3}}\end{bmatrix}.$ (31)

The covariant and contravariant metric coefficients, $g_{ij}$ and $g^{ij}$
respectively, can be calculated by the scalar products of the covariant and
contravariant bases in general curvilinear coordinates with non-orthonormal
bases [48] as follows

$\displaystyle g_{ij}$
$\displaystyle=\boldsymbol{g}_{i}\cdot\boldsymbol{g}_{j}=\frac{\partial
x^{k}}{\partial\xi^{i}}\frac{\partial
x^{k}}{\partial\xi^{j}}\Rightarrow\boldsymbol{\left[G\right]}=\boldsymbol{\left[A\right]^{T}}\boldsymbol{\left[A\right]}=\boldsymbol{\left[G\right]^{T}},$
(32a) $\displaystyle g^{ij}$
$\displaystyle=\boldsymbol{g}^{i}\cdot\boldsymbol{g}^{j}=\frac{\partial\xi^{i}}{\partial
x^{k}}\frac{\partial\xi^{j}}{\partial
x^{k}}\Rightarrow\boldsymbol{\left[\bar{G}\right]}=\boldsymbol{\left[\bar{A}\right]}\boldsymbol{\left[\bar{A}\right]^{T}}=\boldsymbol{\left[\bar{G}\right]^{T}},$
(32b)

in addition to the physical covariant base vectors which can be evaluated
using

$\displaystyle\boldsymbol{e_{i}}=\left(1/\sqrt{g_{ii}}\right)\boldsymbol{g_{i}}.$
(33)

The metric coefficients expressed in Eq. 32 can be employed to alternate
between the covariant and contravariant components of a vector. Suppose a
vector field $\boldsymbol{\upsilon}$ is a function of position vector
$\boldsymbol{r}$ at an arbitrary point P in space, $\boldsymbol{\upsilon}$ can
be expressed in terms of components in the: global Cartesian basis
$\\{\boldsymbol{i}_{i}\\}$, local contravariant basis
$\\{\boldsymbol{g}^{i}\\}$, local covariant basis $\\{\boldsymbol{g}_{i}\\}$,
or physical covariant basis $\\{\boldsymbol{e}_{i}\\}$. Using the Einstein
summation notation, this implies

$\displaystyle\begin{split}\boldsymbol{\upsilon}&=u_{i}\boldsymbol{i}_{i}=\upsilon_{i}\boldsymbol{g}^{i}=\upsilon^{i}\boldsymbol{g}_{i}=\hat{\upsilon}^{i}\boldsymbol{e}_{i},\\\
\boldsymbol{\upsilon}&=u\boldsymbol{i}+v\boldsymbol{j}+w\boldsymbol{k}\\\
&=v_{1}\boldsymbol{g}^{1}+v_{2}\boldsymbol{g}^{2}+v_{3}\boldsymbol{g}^{3}\\\
&=v^{1}\boldsymbol{g}_{1}+v^{2}\boldsymbol{g}_{2}+v^{3}\boldsymbol{g}_{3}\\\
&=\hat{v}^{1}\boldsymbol{e}_{1}+\hat{v}^{2}\boldsymbol{e}_{2}+\hat{v}^{3}\boldsymbol{e}_{3},\end{split}$
(34)

where $u_{i}$, $v_{i}$, $v^{i}$ and $\hat{v}^{i}$ are the Cartesian,
covariant, contravariant, and physical components of $\boldsymbol{\upsilon}$,
respectively. The relationships between these components can be written as

$\displaystyle v_{i}=\left(\frac{\partial
x^{j}}{\partial\xi^{i}}\right)u_{j}=A^{j}_{i}u_{j}=g_{ij}v^{j}\Rightarrow\begin{bmatrix}v_{1}\\\
v_{2}\\\ v_{3}\end{bmatrix}=\left[\boldsymbol{A}\right]^{T}\begin{bmatrix}u\\\
v\\\ w\end{bmatrix}=\left[\boldsymbol{G}\right]\begin{bmatrix}v^{1}\\\
v^{2}\\\ v^{3}\end{bmatrix},$

$\displaystyle v^{i}=\left(\frac{\partial\xi^{i}}{\partial
x^{j}}\right)u_{j}=\bar{A}^{i}_{j}u_{j}=g^{ij}v_{j}\Rightarrow\begin{bmatrix}v^{1}\\\
v^{2}\\\
v^{3}\end{bmatrix}=\left[\boldsymbol{\bar{A}}\right]\begin{bmatrix}u\\\ v\\\
w\end{bmatrix}=\left[\boldsymbol{\bar{G}}\right]\begin{bmatrix}v_{1}\\\
v_{2}\\\ v_{3}\end{bmatrix},$
$\displaystyle\hat{v}^{i}=\sqrt{g_{ii}}v^{i}\Rightarrow\begin{bmatrix}\hat{v}^{1}&\hat{v}^{2}&\hat{v}^{3}\end{bmatrix}^{T}=\sqrt{diag\left[\boldsymbol{G}\right]}\begin{bmatrix}v^{1}&v^{2}&v^{3}\end{bmatrix}^{T},$
$\displaystyle\begin{split}u_{i}&=\left(\frac{\partial\xi^{j}}{\partial
x^{i}}\right)v_{j}=\bar{A}^{j}_{i}v_{j}\Rightarrow\begin{bmatrix}u&v&w\end{bmatrix}^{T}=\left[\boldsymbol{\bar{A}}\right]^{T}\begin{bmatrix}v_{1}&v_{2}&v_{3}\end{bmatrix}^{T}\\\
&=\left(\frac{\partial
x^{i}}{\partial\xi^{j}}\right)v^{j}=A^{i}_{j}v^{j}\Rightarrow\begin{bmatrix}u&v&w\end{bmatrix}^{T}=\left[\boldsymbol{A}\right]^{T}\begin{bmatrix}v^{1}&v^{2}&v^{3}\end{bmatrix}^{T}\\\
&=\left(\frac{\partial
x^{i}}{\partial\xi^{j}}\right)\frac{\hat{v}^{j}}{\sqrt{g_{ii}}}=\left(\frac{A^{i}_{j}}{\sqrt{g_{jj}}}\right)\hat{v}^{j}.\end{split}$

To represent the conservation equations in the strong conservation form using
generalised coordinates, it is convenient to express the flow velocity field
with the Cartesian components $u_{i}\left(u,v,w\right)$ and use the
contravariant components $v^{i}\left(v^{1},v^{2},v^{3}\right)$ for the
advection velocity [4, 10, 11]. Therefore, it is only necessary to transform
the spatial derivatives with respect to Cartesian coordinates into generalised
coordinates. For a general scalar field $\varphi$, this transformation can be
written as

$\displaystyle\frac{\partial\varphi}{\partial
x^{i}}=\frac{\partial\xi^{j}}{\partial
x^{i}}\frac{\partial\varphi}{\partial\xi^{j}}=\bar{A}^{j}_{i}\frac{\partial\varphi}{\partial\xi^{j}}\Rightarrow\begin{bmatrix}\dfrac{\partial\varphi}{\partial
x}&\dfrac{\partial\varphi}{\partial y}&\dfrac{\partial\varphi}{\partial
z}\end{bmatrix}^{T}=\left[\boldsymbol{\bar{A}}\right]^{T}\begin{bmatrix}\dfrac{\partial\varphi}{\partial\xi}&\dfrac{\partial\varphi}{\partial\eta}&\dfrac{\partial\varphi}{\partial\gamma}\end{bmatrix}^{T}.$
(35)

##  Appendix B General Strong Scalar Conservation Equation

The general conservation equation for an arbitrary scalar physical quantity
$\phi$ can be written in coordinate-free vector form as [29]

$\displaystyle\frac{\partial\left(\rho\phi\right)}{\partial
t}+\boldsymbol{\nabla}\cdot\left(\rho\boldsymbol{\nu}\phi\right)=\boldsymbol{\nabla}\cdot\left(\Gamma\boldsymbol{\nabla}\phi\right)+S_{\phi},$
(36)

where $\rho$ is the fluid density, $\boldsymbol{\nu}$ is the velocity vector
field of the fluid, $\Gamma$ is the diffusion coefficient, and $S_{\phi}$
denotes the source/sink term. This equation is often referred to as a strong
conservation (or divergent) form of the governing transport equation of
$\phi$. In Cartesian coordinates $x_{i}(x,y,z)$, this equation is expressed in
tensor notation as [29]

$\displaystyle\frac{\partial\left(\rho\phi\right)}{\partial
t}+\frac{\partial}{\partial x^{j}}\left(\rho
u_{j}\phi\right)=\frac{\partial}{\partial
x^{j}}\left(\Gamma\frac{\partial\phi}{\partial x^{j}}\right)+S_{\phi},$ (37)

where $u_{i}(u,v,w)$ are the Cartesian components of the velocity vector,
i.e., $\boldsymbol{\nu}=u\boldsymbol{i}+v\boldsymbol{j}+w\boldsymbol{k}$,
$\\{\boldsymbol{i},\boldsymbol{j},\boldsymbol{k}\\}$ being the unit vectors in
the direction of the Cartesian axes (Cartesian basis).

Considering a generalised curvilinear coordinate system defined by the time-
independent coordinate transformation given in Eq. A, the general scalar
conservation equation may be written in the contravariant components form as
[4, 9, 10, 11, 12, 49]

$\displaystyle\frac{\partial\left(J\rho\phi\right)}{\partial
t}+\frac{\partial}{\partial\xi^{j}}\left(J\rho\nu^{j}\phi\right)=\frac{\partial}{\partial\xi^{j}}\left[Jg^{kj}\Gamma\left(\frac{\partial\phi}{\partial\xi^{k}}\right)\right]+JS_{\phi},\quad\quad\quad
g^{kj}=\frac{\partial\xi^{k}}{\partial x^{m}}\frac{\partial\xi^{j}}{\partial
x^{m}},$ (38)

where $\nu^{j}$ are the contravariant components of the flow velocity vector
in the local basis of the tangential vectors $\\{\boldsymbol{g_{j}}\\}$
(covariant basis). These components are proportional to the physical
components of the velocity vector tangential to the coordinate curve $\xi^{j}$
at any point in space.

##  Appendix C Mass Conservation Equation (Continuity Equation)

The mass conservation equation for compressible fluids is expressed in
coordinate-free vector form as [50]

$\displaystyle\frac{D\rho}{Dt}+\rho\boldsymbol{\nabla}\cdot\boldsymbol{v}=0\Leftrightarrow\frac{\partial\rho}{\partial
t}+\boldsymbol{\nabla}\cdot\left(\rho\boldsymbol{v}\right)=0,$ (39)

where $\rho$ is the fluid density and $\boldsymbol{v}$ the velocity vector
field of the fluid. This implies that the continuity equation in Cartesian
coordinates can be written in tensor notation as

$\displaystyle\frac{D\rho}{Dt}+\rho\frac{\partial u_{i}}{\partial
x_{i}}=0\Leftrightarrow\frac{\partial\rho}{\partial
t}+\frac{\partial\left(\rho u_{i}\right)}{\partial x_{i}}=0.$ (40)

Therefore, considering the coordinates transformation described in Appendix A
and the normalisation of the dependent variables summarized in Appenidx G, the
dimensionless continuity equation derived in the generalised transformed
domain can be expressed in tensor notation as

$\displaystyle\frac{\partial\bar{\rho}}{\partial\tau}+\nabla_{X}\cdot\left(\bar{\rho}\boldsymbol{\tilde{V}}\right)=0,$
(41)

where the components of
$\boldsymbol{\tilde{V}}=\left[\tilde{V}^{X}\;\tilde{V}^{Y}\;\tilde{V}^{Z}\right]$
are evaluated using Eqs. 12a, 12b and 12c.

##  Appendix D Integral Terms of the Generalised Reynolds Equation and
Velocity Components

The integral terms of the p-$\theta$ modified generalised Reynolds equation
and velocity components expressed using Cartesian coordinates are defined as

$\displaystyle\frac{1}{\eta_{e}}=\bigintssss_{z_{1}}^{z_{2}}\frac{dz}{\eta},$
(42a)
$\displaystyle\frac{1}{\eta_{e}^{{}^{\prime}}}=\bigintssss_{z_{1}}^{z_{2}}\frac{z\;dz}{\eta},$
(42b)
$\displaystyle\rho^{{}^{\prime}}=\bigintssss_{z_{1}}^{z_{2}}\rho\;\left(\bigintssss_{z_{1}}^{z}\frac{dz^{\prime}}{\eta}\right)\;dz,$
(42c)
$\displaystyle\rho^{{}^{\prime\prime}}=\bigintssss_{z_{1}}^{z_{2}}\rho\;\left(\bigintssss_{z_{1}}^{z}\frac{z^{\prime}dz^{\prime}}{\eta}\right)\;dz,$
(42d) $\displaystyle
I(x,y,z)=\bigintssss_{z_{1}}^{z}\frac{z^{\prime}\;dz^{\prime}}{\eta},$ (42e)
$\displaystyle J(x,y,z)=\bigintssss_{z_{1}}^{z}\frac{dz^{\prime}}{\eta}.$
(42f) $\displaystyle
I_{p}(x,y,z)=\left(I-J\frac{\eta_{e}}{\eta_{e}^{{}^{\prime}}}\right),$ (42g)
$\displaystyle I_{s}(x,y,z)=J\eta_{e}.$ (42h)

On the other hand, the coefficients and integral terms of p-$\theta$ mass-
conserving generalised Reynolds equation and velocity components expressed
using normalised transformed coordinates are given as

$\displaystyle\bar{I}(X,Y,Z)=\bigintssss_{0}^{Z}\frac{Z^{\prime}\;dZ^{\prime}}{\bar{\eta}},$
(43a)
$\displaystyle\bar{J}(X,Y,Z)=\bigintssss_{0}^{Z}\frac{dZ^{\prime}}{\bar{\eta}},$
(43b)
$\displaystyle\bar{I}_{1}(X,Y)=\bar{I}(X,Y,1)=\bigintssss_{0}^{1}\frac{Z\;dZ}{\bar{\eta}},$
(43c)
$\displaystyle\bar{J}_{1}(X,Y)=\bar{J}(X,Y,1)=\bigintssss_{0}^{1}\frac{dZ}{\bar{\eta}},$
(43d)
$\displaystyle\bar{I}_{p}(X,Y,Z)=H^{2}\left[\bar{I}(X,Y,Z)-\bar{J}(X,Y,Z)\frac{\bar{I}_{1}(X,Y)}{\bar{J}_{1}(X,Y)}\right],$
(43e)
$\displaystyle\bar{I}_{s}(X,Y,Z)=\frac{\bar{J}(X,Y,Z)}{\bar{J}_{1}(X,Y)},$
(43f)
$\displaystyle\bar{R}_{1}=\bigintssss_{0}^{1}\bar{\rho}\;\left(\bigintssss_{0}^{Z}\frac{dZ^{\prime}}{\bar{\eta}}\right)\;dZ,$
(43g)
$\displaystyle\bar{R}_{2}=\bigintssss_{0}^{1}\bar{\rho}\;\left(\bigintssss_{0}^{Z}\frac{Z^{\prime}\;dZ^{\prime}}{\bar{\eta}}\right)\;dZ,$
(43h) $\displaystyle\bar{R}_{3}=\bigintssss_{0}^{1}\bar{\rho}\;dZ,$ (43i)
$\displaystyle\bar{\varepsilon}(X,Y)=H^{3}\left[\frac{\bar{I}_{1}(X,Y)\bar{R}_{1}(X,Y)}{\bar{J}_{1}(X,Y)}-\bar{R}_{2}(X,Y)\right],$
(43j)
$\displaystyle\bar{\rho}^{*}_{e}(X,Y)=\frac{2\bar{R}_{1}(X,Y)}{\bar{J}_{1}(X,Y)},$
(43k)
$\displaystyle\bar{\rho}^{*}_{1}(X,Y)=\left[\bar{R}_{3}(X,Y)-\frac{2\bar{R}_{1}(X,Y)}{\bar{J}_{1}(X,Y)}\right],$
(43l) $\displaystyle\bar{\rho}_{e}(X,Y)=\bar{R}_{3}(X,Y).$ (43m)

##  Appendix E Coefficients of the Discrete $(p-\theta)$ Generalised Reynolds
Equation in the Normalised Transformed Domain

Assuming the variables computed at the CV faces to be equivalent to the
average of the pertinent variables at the surrounding nodes, _e.g._
$\left(\bar{\varepsilon}\right)_{w}=\biggl{[}\left(\bar{\varepsilon}\right)_{W}+\left(\bar{\varepsilon}\right)_{P}\biggr{]}\;\bigg{/}\;2=\biggl{[}\left(\bar{\varepsilon}\right)_{k,i-1}+\left(\bar{\varepsilon}\right)_{k,i}\biggr{]}\;\bigg{/}\;2$,
hence the coefficients in Eq. 16 are defined as

$\displaystyle\alpha_{D}$
$\displaystyle=\bigg{[}\left(\bar{\varepsilon}\right)_{w}\cdot\left(\delta
Y\right)_{w}\big{/}\left(\delta X\right)_{w}\bigg{]}$ (44a)
$\displaystyle\alpha_{C}$ $\displaystyle=-\left(\delta
Y\right)_{w}\cdot\bigg{[}\left(\bar{\rho}^{*}_{e}U_{m}+\bar{\rho}^{*}_{1}U_{1}\right)_{i}\cdot\theta_{i}\cdot
D^{i,j}_{k-1,l}-\left(\bar{\rho}^{*}_{e}U_{m}+\bar{\rho}^{*}_{1}U_{1}\right)_{i-1}\cdot\theta_{i-1}\cdot
D^{i-1,j}_{k-1,l}\bigg{]}$ (44b) $\displaystyle\alpha$
$\displaystyle=\alpha_{D}+\alpha_{C}$ (44c)
$\displaystyle\begin{split}\beta_{D}&=-\bigg{[}\left(\bar{\varepsilon}\right)_{w}\cdot\left(\delta
Y\right)_{w}\big{/}\left(\delta
X\right)_{w}\bigg{]}-\bigg{[}\left(\bar{\varepsilon}\right)_{e}\cdot\left(\delta
Y\right)_{e}\big{/}\left(\delta X\right)_{e}\bigg{]}\\\
&-r_{xy}^{2}\bigg{[}\left(\bar{\varepsilon}\right)_{n}\cdot\left(\delta
X\right)_{n}\big{/}\left(\delta
Y\right)_{n}\bigg{]}-r_{xy}^{2}\bigg{[}\left(\bar{\varepsilon}\right)_{s}\cdot\left(\delta
X\right)_{s}\big{/}\left(\delta Y\right)_{s}\bigg{]}\end{split}$ (45a)
$\displaystyle\beta_{C}$ $\displaystyle=-\left(\delta
Y\right)_{w}\cdot\bigg{[}\left(\bar{\rho}^{*}_{e}U_{m}+\bar{\rho}^{*}_{1}U_{1}\right)_{i}\cdot\theta_{i}\cdot
D^{i,j}_{k,l}-\left(\bar{\rho}^{*}_{e}U_{m}+\bar{\rho}^{*}_{1}U_{1}\right)_{i-1}\cdot\theta_{i-1}\cdot
D^{i-1,j}_{k,l}\bigg{]}$ (45b) $\displaystyle\beta$
$\displaystyle=\beta_{D}+\beta_{C}$ (45c) $\displaystyle\gamma_{D}$
$\displaystyle=\bigg{[}\left(\bar{\varepsilon}\right)_{e}\cdot\left(\delta
Y\right)_{e}\big{/}\left(\delta X\right)_{e}\bigg{]}$ (46a)
$\displaystyle\gamma_{C}$ $\displaystyle=-\left(\delta
Y\right)_{e}\cdot\bigg{[}\left(\bar{\rho}^{*}_{e}U_{m}+\bar{\rho}^{*}_{1}U_{1}\right)_{i}\cdot\theta_{i}\cdot
D^{i,j}_{k+1,l}-\left(\bar{\rho}^{*}_{e}U_{m}+\bar{\rho}^{*}_{1}U_{1}\right)_{i-1}\cdot\theta_{i-1}\cdot
D^{i-1,j}_{k+1,l}\bigg{]}$ (46b) $\displaystyle\gamma$
$\displaystyle=\gamma_{D}+\gamma_{C}$ (46c) $\displaystyle\zeta_{D}$
$\displaystyle=0$ (47a) $\displaystyle\zeta_{C}$ $\displaystyle=0$ (47b)
$\displaystyle\zeta$ $\displaystyle=\zeta_{D}+\zeta_{C}$ (47c)
$\displaystyle\lambda_{D}$ $\displaystyle=0$ (48a) $\displaystyle\lambda_{C}$
$\displaystyle=0$ (48b) $\displaystyle\lambda$
$\displaystyle=\lambda_{D}+\lambda_{C}$ (48c) $\displaystyle\varphi_{D}$
$\displaystyle=\left(\bar{p}\right)_{N}\cdot
r^{2}_{xy}\cdot\bigg{[}\left(\bar{\varepsilon}\right)_{n}\cdot\left(\delta
X\right)_{n}\big{/}\left(\delta
Y\right)_{n}\bigg{]}+\left(\bar{p}\right)_{S}\cdot
r^{2}_{xy}\cdot\bigg{[}\left(\bar{\varepsilon}\right)_{n}\cdot\left(\delta
X\right)_{s}\big{/}\left(\delta Y\right)_{s}\bigg{]}$ (49a)
$\displaystyle\begin{split}\varphi_{C}&=\left(\bar{\rho}^{*}_{e}U_{m}+\bar{\rho}^{*}_{1}U_{1}\right)_{i}\cdot\theta_{i}\cdot\Bigg{[}H_{i}-\bigg{(}\big{[}D^{i,j}_{k-1,l}\cdot\bar{p}_{k-1,l}\big{]}-\big{[}D^{i,j}_{k,l}\cdot\bar{p}_{k,l}\big{]}-\big{[}D^{i,j}_{k+1,l}\cdot\bar{p}_{k+1,l}\big{]}\bigg{)}\Bigg{]}\\\
&-\left(\bar{\rho}^{*}_{e}U_{m}+\bar{\rho}^{*}_{1}U_{1}\right)_{i-1}\cdot\theta_{i-1}\cdot\Bigg{[}H_{i-1}-\bigg{(}\big{[}D^{i-1,j}_{k-1,l}\cdot\bar{p}_{k-1,l}\big{]}-\big{[}D^{i-1,j}_{k,l}\cdot\bar{p}_{k,l}\big{]}-\big{[}D^{i-1,j}_{k+1,l}\cdot\bar{p}_{k+1,l}\big{]}\bigg{)}\Bigg{]}\end{split}$
(49b) $\displaystyle\varphi$ $\displaystyle=\varphi_{D}+\varphi_{C}$ (49c)

##  Appendix F Coefficients of the Discrete Energy Equation in the Normalised
Transformed Domain

Using the $\big{[}\big{[}\;\big{]}\big{]}$ operator of Pantankar where
$\big{[}\big{[}A,B\big{]}\big{]}$ denotes the greater of $A$ and $B$ [51], and
assuming the variables computed at the CV faces to be equivalent to the
average of the pertinent variables at the surrounding nodes, _e.g._
$\left(\bar{k}\right)_{w}=\biggl{[}\left(\bar{k}\right)_{W}+\left(\bar{k}\right)_{P}\biggr{]}\;\bigg{/}\;2=\biggl{[}\left(\bar{k}\right)_{k,i-1}+\left(\bar{k}\right)_{k,i}\biggr{]}\;\bigg{/}\;2$,
hence the coefficients in Eq. 19 are defined as

$\displaystyle\alpha_{D}$
$\displaystyle=-\left(\varepsilon\right)^{2}\cdot\bar{k}_{w}\cdot
H_{w}\cdot\mathcal{S}_{w}\big{/}\left(\delta X\right)_{w}$ (50a)
$\displaystyle\alpha_{C}$
$\displaystyle=-\;\big{[}\big{[}\left(\bar{\rho}\bar{c}_{p}\tilde{V}^{X}\right)_{w}\cdot\mathrm{Pe}\cdot\mathcal{S}_{w},\;0\big{]}\big{]}$
(50b) $\displaystyle\alpha$ $\displaystyle=\alpha_{D}+\alpha_{C}$ (50c)
$\displaystyle\gamma_{D}$
$\displaystyle=-\left(\varepsilon\right)^{2}\cdot\bar{k}_{e}\cdot
H_{e}\cdot\mathcal{S}_{e}\big{/}\left(\delta X\right)_{e}$ (51a)
$\displaystyle\gamma_{C}$
$\displaystyle=-\;\big{[}\big{[}-\left(\bar{\rho}\bar{c}_{p}\tilde{V}^{X}\right)_{e}\cdot\mathrm{Pe}\cdot\mathcal{S}_{e},\;0\big{]}\big{]}$
(51b) $\displaystyle\gamma$ $\displaystyle=\gamma_{D}+\gamma_{C}$ (51c)
$\displaystyle\zeta_{D}$
$\displaystyle=-\left(r_{xy}\varepsilon\right)^{2}\cdot\bar{k}_{n}\cdot
H_{n}\cdot\mathcal{S}_{n}\big{/}\left(\delta Y\right)_{n}$ (52a)
$\displaystyle\zeta_{C}$
$\displaystyle=-\;\big{[}\big{[}-\left(\bar{\rho}\bar{c}_{p}\tilde{V}^{Y}\right)_{n}\cdot\mathrm{Pe}\cdot\mathcal{S}_{n},\;0\big{]}\big{]}$
(52b) $\displaystyle\zeta$ $\displaystyle=\zeta_{D}+\zeta_{C}$ (52c)
$\displaystyle\lambda_{D}$
$\displaystyle=-\left(r_{xy}\varepsilon\right)^{2}\cdot\bar{k}_{s}\cdot
H_{s}\cdot\mathcal{S}_{s}\big{/}\left(\delta Y\right)_{s}$ (53a)
$\displaystyle\lambda_{C}$
$\displaystyle=-\;\big{[}\big{[}\left(\bar{\rho}\bar{c}_{p}\tilde{V}^{Y}\right)_{s}\cdot\mathrm{Pe}\cdot\mathcal{S}_{s},\;0\big{]}\big{]}$
(53b) $\displaystyle\lambda$ $\displaystyle=\lambda_{D}+\lambda_{C}$ (53c)
$\displaystyle\phi_{D}$
$\displaystyle=-\left(\bar{k}_{t}\cdot\mathcal{S}_{t}\right)\cdot\bigg{[}1+\left(\varepsilon^{2}\cdot\left[\left(\nabla_{X}\cdot
H\right)_{t}\right]^{2}\cdot Z_{t}^{2}\right)+\left(\left(\varepsilon
r_{xy}\right)^{2}\cdot\left[\left(\nabla_{Y}\cdot H\right)_{t}\right]^{2}\cdot
Z_{t}^{2}\right)\bigg{]}\bigg{/}\left(H\cdot\left(\delta Z\right)_{t}\right)$
(54a) $\displaystyle\phi_{C}$
$\displaystyle=-\;\big{[}\big{[}-\left(\bar{\rho}\bar{c}_{p}\tilde{V}^{Z}\right)_{t}\cdot\mathrm{Pe}\cdot\mathcal{S}_{t},\;0\big{]}\big{]}$
(54b) $\displaystyle\phi$ $\displaystyle=\phi_{D}+\phi_{C}$ (54c)
$\displaystyle\psi_{D}$
$\displaystyle=-\left(\bar{k}_{b}\cdot\mathcal{S}_{b}\right)\cdot\bigg{[}1+\left(\varepsilon^{2}\cdot\left[\left(\nabla_{X}\cdot
H\right)_{b}\right]^{2}\cdot Z_{b}^{2}\right)+\left(\left(\varepsilon
r_{xy}\right)^{2}\cdot\left[\left(\nabla_{Y}\cdot H\right)_{b}\right]^{2}\cdot
Z_{b}^{2}\right)\bigg{]}\bigg{/}\left(H\cdot\left(\delta Z\right)_{b}\right)$
(55a) $\displaystyle\psi_{C}$
$\displaystyle=-\;\big{[}\big{[}\left(\bar{\rho}\bar{c}_{p}\tilde{V}^{Z}\right)_{b}\cdot\mathrm{Pe}\cdot\mathcal{S}_{b},\;0\big{]}\big{]}$
(55b) $\displaystyle\psi$ $\displaystyle=\psi_{D}+\psi_{C}$ (55c)
$\displaystyle\beta_{D}$
$\displaystyle=-\alpha_{D}-\gamma_{D}-\zeta_{D}-\lambda_{D}-\phi_{D}-\psi_{D}$
(56a)
$\displaystyle\begin{split}\beta_{C}&=\big{[}\big{[}\left(\bar{\rho}\bar{c}_{p}\tilde{V}^{X}\right)_{e}\cdot\mathrm{Pe}\cdot\mathcal{S}_{e},\;0\big{]}\big{]}+\big{[}\big{[}\left(\bar{\rho}\bar{c}_{p}\tilde{V}^{Y}\right)_{n}\cdot\mathrm{Pe}\cdot\mathcal{S}_{n},\;0\big{]}\big{]}+\big{[}\big{[}\left(\bar{\rho}\bar{c}_{p}\tilde{V}^{Z}\right)_{t}\cdot\mathrm{Pe}\cdot\mathcal{S}_{t},\;0\big{]}\big{]}\\\
&+\big{[}\big{[}-\left(\bar{\rho}\bar{c}_{p}\tilde{V}^{X}\right)_{w}\cdot\mathrm{Pe}\cdot\mathcal{S}_{w},\;0\big{]}\big{]}+\big{[}\big{[}-\left(\bar{\rho}\bar{c}_{p}\tilde{V}^{Y}\right)_{s}\cdot\mathrm{Pe}\cdot\mathcal{S}_{s},\;0\big{]}\big{]}+\big{[}\big{[}-\left(\bar{\rho}\bar{c}_{p}\tilde{V}^{Z}\right)_{b}\cdot\mathrm{Pe}\cdot\mathcal{S}_{b},\;0\big{]}\big{]}\end{split}$
(56b) $\displaystyle\beta$ $\displaystyle=\beta_{D}+\beta_{C}$ (56c)
$\displaystyle\varphi$
$\displaystyle=\biggl{[}\left(\bar{\tilde{Q}}_{p}\right)_{P}+\left(\bar{\tilde{Q}}_{cp}\right)_{P}+\left(\bar{\tilde{Q}}_{\Phi}\right)_{P}+\left(\bar{\tilde{q}}_{V}\right)_{P}\biggr{]}\cdot\left(\Delta\mathcal{V}\right)_{P}$
(57a)

##  Appendix G Transformation and Dimensionless Parameters

Parameters used in the transformation of the governing equations from the
Cartesian coordinate system into the normalised transformed curvilinear grid
system, as well as other dimensional and dimensionless parameters relevant to
the TEHL framework.

Table 1: Pertinent symbols and their respective definitions and mathematical expressions. Symbol |  Parameter |  Value
---|---|---
$A^{i}_{j}$, $\bar{A^{i}_{j}}$ |  Coordinate transformation and inverse coordinate transformation matrix |  Eq. 27, 28
$A_{\scriptscriptstyle X}$ |  Normalised curvature in the $x$-direction |  $A_{\scriptscriptstyle X}=\left({L_{x}}\right)^{2}\left({2R_{x}h_{0}}\right)^{-1}$
$A_{\scriptscriptstyle Y}$ |  Normalised curvature in the $y$-direction |  $A_{Y}=\left({L_{y}}\right)^{2}\left({2R_{y}h_{0}}\right)^{-1}$
$a$ |  Hertz contact half-width of circular point contact [m] |  $a=\sqrt[3]{\left(3W_{L}R^{{}^{\prime}}\right)\big{/}\left(2E^{{}^{\prime}}\right)}$
$\mathrm{Br}$, $\mathrm{Br^{*}}$ |  Brinckman and modified Brinckman number |  $\mathrm{Br}=\dfrac{\eta_{0}u^{2}_{e}}{k_{0}T_{0}},\mathrm{Br^{*}}=\left(\beta_{0}T_{0}Br\right)$
$c_{p}$, $\bar{c}_{p}$ |  Dimensional $\mathrm{[J/kg.K]}$ and normalised specific heat capacity |  $\bar{c}_{p}={c_{p}}*\left({c_{p_{0}}}\right)^{-1}$
$D^{i,j}_{k,l}$ |  Elements of influence coefficient matrix |  -
$E_{1}$, $E_{2}$ |  Lower (1) and upper (2) solid elastic modulus [Pa] |  -
$E^{{}^{\prime}}$ |  Effective elastic modulus [Pa] |  $E^{\prime}=2\left[(1-{\nu_{1}}^{2})/E_{1})+(1-{\nu_{2}}^{2})/E_{2}\right]^{-1}$
$\bar{G}$ |  Dimensionless material parameter |  $\bar{G}=\alpha*E^{{}^{\prime}}$
$\boldsymbol{g_{i}}$, $\boldsymbol{g^{i}}$ |  Covariant and contravariant base vectors |  Eq. 30, 31
$\boldsymbol{g_{ij}}$, $\boldsymbol{g^{ij}}$ |  Covariant and contravariant metric coefficients |  Eq. 32
$h$, $H$ |  Dimensional [m] and normalised film thickness |  $H=h*\left(h_{0}\right)^{-1}=\left(z_{2}-z_{1}\right)*\left(h_{0}\right)^{-1}$
$h_{0}$, $H_{0}$ |  Dimensional [m] and normalised rigid separation between two bodies prior to deformation |  $H_{0}={h_{0}}/{R}=2.69\bar{U}^{0.67}\bar{G}^{0.53}\bar{W}^{-0.067}\left(1-0.61e^{-0.73}\right)$
${J}$ |  Jacobian |  Eq. 27
$k$, $\bar{k}$ |  Dimensional [W/m.K] and normalised thermal conductivity |  $\bar{k}={k}*\left({k_{0}}\right)^{-1}$
$L_{x}$, $L_{y}$ |  Characteristic length of the contact in the sliding and transverse direction, respectively |  ($L_{x}$ = $L_{y}=a$ for circular contact)
$N_{x}$, $N_{y}$ |  Mesh density in the $x$\- and $y$-directions, respectively |  -
$\vec{\boldsymbol{n}}$ |  Normal unit vector |  -
$nb$ |  Neighbouring nodes |  -
$\mathrm{Pe}$ |  Peclet number |  $\mathrm{Pe}=\left(\rho_{0}c_{p_{0}}u_{e}L_{x}\varepsilon^{2}\right)\big{/}k_{0}$
$\mathrm{Pr}$ |  Prandtl number |  $\mathrm{Pr}=\left({\eta_{0}c_{p_{0}}}\right)\big{/}k_{0}$
$p$, $\bar{p}$ |  Dimensional [Pa] and normalised pressure |  $\bar{p}=p*\biggl{[}\left(\varepsilon^{2}L_{x}\right)\big{/}\left(\eta_{0}u_{e}\right)\biggr{]}$
$p_{cav}$ |  Cavitation pressure [Pa] |  -
$\dot{q}_{V}$, $\bar{\dot{q}}_{V}$ |  Dimensional $\mathrm{[W/m^{3}]}$ and normalised rate of heat source/sink |  $\bar{\dot{q}}_{V}=\dot{q}_{V}\left(\varepsilon^{2}L_{x}^{2}\right)\left(k_{0}T_{0}\right)^{-1}$
$R^{\prime}$ |  Effective radius of curvature [m] |  ($R^{\prime}=R_{x}=R_{y}$ for circular contact)
$R_{x}$ |  Equivalent radii of curvature of the solids (1,2) at the centre of contact in the $x$-direction [m] |  $R_{x}={(R_{x1}R_{x2})}*{(R_{x1}+R_{x2})^{-1}}$
$R_{y}$ |  Equivalent radii of curvature of the solids (1,2) at the centre of contact in the $y$-direction [m] |  $R_{y}={(R_{y1}R_{y2})}*{(R_{y1}+R_{y2})^{-1}}$
$r_{xy}$ |  Length aspect ratio |  $r_{xy}=L_{x}*\left(L_{y}\right)^{-1}$
$SRR$ |  Slide-to-roll ratio |  $SRR=u_{s}\big{/}u_{e}$
$SRR_{x}$, $SRR_{y}$ |  Ratio of sliding speed in $x$\- and $y$-direction to entrainmemnt speed |  $SRR_{x}=u_{s}\big{/}u_{e}$, $SRR_{y}=v_{s}\big{/}u_{e}$
$T$, $\bar{T}$ |  Dimensional [K] and normalised temperature |  $\bar{T}=T*\left(T_{0}\right)^{-1}$
$T_{0}$ |  Ambient temperature [K] |  -
$t$, $\tau$ |  Dimensional [s] and normalised time |  $\tau=t*\left(u_{e}/L_{x}\right)$
$U_{1}$, $V_{1}$ |  Normalised lower solid velocity in the $x$\- and $y$-directions |  $U_{1}=u_{1}/u_{e}$, $V_{1}=v_{1}/u_{e}$
$\bar{U}$ |  Dimensionless speed parameter |  $\bar{U}=u_{m}\eta_{0}*\left(E^{{}^{\prime}}R^{{}^{\prime}}\right)^{-1}$
$u$, $v$, $w$ |  Flow velocity components in the physical Cartesian space, $x$-, $y$\- and $z$-directions, respectively [m/s] |  Eq. 11a, 11b, 11c
$u_{e}$, $v_{e}$ |  Entertainment velocity in sliding and transverse directions [m/s] |  $u_{e}=\left(u_{1}+u_{2}\right)/2$, $v_{e}=\left(v_{1}+v_{2}\right)/2$
$u_{m}(x,y)$ |  General mean velocity of surfaces [m/s], $u_{m}(x,y)=u_{e}$ if $u_{1}$ and $u_{2}$ are constant in space |  $u_{m}(x,y)=\bigl{[}u_{1}(x,y)+u_{2}(x,y)\bigr{]}\bigg{/}2$
$u_{s}$, $v_{s}$ |  Sliding velocity in sliding and transverse directions |  $u_{s}=\left(u_{2}-u_{1}\right)$, $v_{s}=\left(v_{2}-v_{1}\right)$
$u_{1}$, $u_{2}$ |  Lower ($1$) and upper ($2$) solid velocity in the $x$-direction [m/s] |  -
$\tilde{V}^{X}$, $\tilde{V}^{Y}$, $\tilde{V}^{Z}$ |  Contravariant velocity components in the normalised transformed domain, $X$-, $Y$\- and $Z$-directions, respectively |  Eq. 12a, 12b, 12c
$v_{1}$, $v_{2}$ |  Lower ($1$) and upper ($2$) solid velocity in the $y$-direction [m/s] |  -
$\bar{W}$ |  Dimensionless load parameter |  $\bar{W}=W_{L}*\left(E^{{}^{\prime}}{R^{{}^{\prime}}}^{2}\right)^{-1}$
$W_{L}$ |  Applied load $\mathrm{[N]}$ |  -
$x$, $X$ |  Dimensional [m] and normalised length in sliding direction |  $X=x*\left(L_{x}\right)^{-1}$
$y$, $Y$ |  Dimensional [m] and normalised length in transverse direction |  $Y=y*\left(L_{y}\right)^{-1}$
$Z_{1}$, $Z_{2}$ |  Lower (1) and upper (2) normalised solid height with respect to the reference origin |  $Z_{1}=z_{1}*\left(h_{0}\right)^{-1}$, $Z_{2}=z_{2}*\left(h_{0}\right)^{-1}$
$z$, $Z$ |  Dimensional [m] and normalised length in film thickness direction |  $Z=\left(z-z_{1}\right)*\left(h(x,y)\right)^{-1}$
$z_{1}$, $z_{2}$ |  Lower (1) and upper (2) solid height with respect to the reference origin |  -
$\alpha$ |  Pressure-viscosity coefficient $\mathrm{[m^{2}/N]}$ |  -
$\beta$, $\bar{\beta}$ |  Dimensional $\mathrm{[m^{2}/N]}$ and normalised compressibility |  $\bar{\beta}=\beta*\left(\beta_{0}\right)^{-1}$
$\dot{\gamma}$, $\bar{\dot{\gamma}}$ |  Dimensional $\mathrm{[s^{-1}]}$ and normalised shear-rate |  $\bar{\dot{\gamma}}=\dot{\gamma}*\biggl{[}\left(h_{0}\right)/\left(u_{e}\right)\biggr{]}$
$\delta X_{i}$, $\delta Y_{i}$, $\delta Z_{i}$ |  Normalised length, width and height of face $i$, respectively |  -
$\delta$, $\bar{\delta}$ |  Dimensional [m] and normalised surface elastic deflection |  $\bar{\delta}=\delta*\left(h_{0}\right)^{-1}$
$\varepsilon$ |  Normalised length scale |  $\varepsilon=h_{0}*\left(L_{x}\right)^{-1}$
$\zeta$, $\eta$, $\gamma$ |  Curvilinear coordinate system |  -
$\eta$, $\bar{\eta}$ |  Dimensional [Pa.s] and normalised dynamic viscosity |  $\bar{\eta}={\eta}*\left({\eta_{0}}\right)^{-1}$
$\eta_{0}$, $\rho_{0}$, $k_{0}$, $c_{p_{0}}$, $\beta_{0}$ |  Ambient dynamic viscosity [Pa.s], density $\mathrm{[kg/m^{3}]}$, conductivity [W/m.K], specific heat capacity $\mathrm{[J/kg.K]}$ and compressibility $\mathrm{[m^{2}/N]}$, respectively |  -
$\theta$ |  Fluid film fraction |  -
$\nu_{1}$, $\nu_{2}$ |  Lower (1) and upper (2) solid Poisson ratio |  -
$\rho$, $\bar{\rho}$ |  Dimensional $\mathrm{[kg/m^{3}]}$ and normalised density |  $\bar{\rho}={\rho}*\left({\rho_{0}}\right)^{-1}$
$\sigma_{0}$ |  Diffusivity $\mathrm{[W.m^{2}.J^{-1}]}$ |  $\sigma_{0}={k_{0}}*\left({\rho_{0}c_{p_{0}}}\right)^{-1}$
$\tau_{e}$ |  Equivalent shear stress $\mathrm{[Pa]}$ |  $\tau_{e}=\sqrt{\tau_{xz}^{2}+\tau_{yz}^{2}}$
$\tau_{0}$ |  Eyring shear stress $\mathrm{[Pa]}$ |  -
$\omega_{min}$, $\omega_{max}$ |  Minimum and maximum PGMA relaxation factor threshold |  -
$\omega_{\psi}^{n}$, $\omega_{\psi}^{n-1}$ |  Relaxation factor at current ($n$) and preceding ($n-1$) iterative step |  -
$\mathcal{S}_{f}$ |  Surface vector normal to the control volume (CV) face $f$ |  -
$\Delta\mathcal{V}$ |  Volume of the control volume (CV) |  $\Delta\mathcal{V}=\left(\delta X\right)\left(\delta Y\right)\left(\delta Z\right)$
$\boldsymbol{\vec{{\mathfrak{R}}}_{\psi}^{n}}$, $\boldsymbol{\vec{{\mathfrak{R}}}_{\psi}^{n-1}}$ |  Residue field at current ($n$) and preceding ($n-1$) iterative step |  -

## References

  * [1] J. Blazek, Computational Fluid Dynamics: Principles and Applications, Butterworth-Heinemann, Oxford, England, UK, 2015.  doi:10.1016/C2013-0-19038-1.
  * [2] H. Versteeg, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Pearson, London, England, UK, 2007.
  * [3] D. Lee, Exploring elastohydrodynamic lubrication using finite-volume computational modelling techniques, Ph.D. thesis, Imperial College, London, England, UK (2018).  doi:10.25560/81590.
  * [4] S. Ogawa, T. Ishiguro, A method for computing flow fields around moving bodies, J. Comput. Phys. 69 (1) (1987) 49–68.  doi:10.1016/0021-9991(87)90155-0.
  * [5] H. Q. Yang, K. T. Yang, J. R. Lloyd, Buoyant flow calculations with non-orthogonal curvilinear co-ordinates for vertical and horizontal parallelepiped enclosures, Int. J. Numer. Methods Eng. 25 (2) (1988) 331–345.  doi:10.1002/nme.1620250205.
  * [6] G. D. Raithby, P. F. Galpin, J. P. Van Doormaal, Prediction of heat and fluid flow in complex geometries using general orthogonal coordinates, Numerical Heat Transfer 9 (2) (1986) 125–142.  doi:10.1080/10407788608913469.
  * [7] G. Cannata, C. Petrelli, L. Barsi, F. Gallerano, Numerical integration of the contravariant integral form of the navier–stokes equations in time-dependent curvilinear coordinate systems for three-dimensional free surface flows, Continuum Mech. Thermodyn. 31 (2) (2019) 491–519.  doi:10.1007/s00161-018-0703-1.
  * [8] I. Jarrah, Rizwan-uddin, Nodal integral methods in general 2d curvilinear coordinates - applied to convection–diffusion equation in domains discretized using quadrilateral elements, Int. J. Heat Mass Transfer 187 (2022) 122559.  doi:10.1016/j.ijheatmasstransfer.2022.122559.
  * [9] T. Kajishima, K. Taira, Computational Fluid Dynamics, Springer International Publishing, Cham, Switzerland, 2017.
  * [10] H. Q. Yang, S. D. Habchi, A. J. Przekwas, General strong conservation formulation of Navier-Stokes equations in nonorthogonal curvilinear coordinates, AIAA Journal 32 (5) (1994) 936–941.  doi:10.2514/3.12077.
  * [11] M. Vinokur, Conservation equations of gasdynamics in curvilinear coordinate systems, J. Comput. Phys. 14 (2) (1974) 105–125.  doi:10.1016/0021-9991(74)90008-4.
  * [12] V. D. Liseikin, Grid Generation Methods (Scientific Computation), Springer, Berlin, Germany, 2017.
  * [13] M. Farrashkhalvat, Basic Structured Grid Generation: With an Introduction to Unstructured Grid Generation, Butterworth-Heinemann, Oxford, England, UK, 2003\.
  * [14] J. Thompson, A survey of grid generation techniques in computational fluid dynamics, in: 21st Aerospace Sciences Meeting, Vol. 1983-447, American Institute of Aeronautics and Astronautics, 1983, pp. 1–36.  doi:10.2514/6.1983-447.
  * [15] T. Han, R. S. Paranjpe, A finite volume analysis of the thermohydrodynamic performance of finite journal bearings, J. Tribol. 112 (3) (1990) 557–565.  doi:10.1115/1.2920293.
  * [16] S. A. Roberts, D. R. Noble, E. M. Benner, P. R. Schunk, Multiphase hydrodynamic lubrication flow using a three-dimensional shell finite element model, Comput. Fluids 87 (2013) 12–25.  doi:10.1016/j.compfluid.2012.08.009.
  * [17] F. Moukalled, L. Mangani, M. Darwish, The Finite Volume Method in Computational Fluid Dynamics, Springer International Publishing, Cham, Switzerland, 2016.
  * [18] D. Dowson, A generalized reynolds equation for fluid-film lubrication, International Journal of Mechanical Sciences 4 (2) (1962) 159–170.  doi:10.1016/S0020-7403(62)80038-1.
  * [19] H. Elrod, M. Adams, A computer program for cavitation and starvation problems, in: Leeds-Lyon Symposium on Tribology, Vol. 1, 1975, pp. 37–41.
  * [20] H. G. Elrod, A cavitation algorithm, J. Lubr. Technol. 103 (3) (1981) 350–354.  doi:10.1115/1.3251669.
  * [21] B. Jakobsson, The finite journal bearing considering vaporization, Trans. Chalmers Univ. of Tech, Sweden 190 (1965).
  * [22] L. Floberg, On hydrodynamic lubricationwith special reference to sub-cavity pressures and number of streamers in cavitation regions, Acta Polytechnica Scandinavica Mechanical Engineering Series 19 (1965).
  * [23] L. Floberg, On the Tensile Strength of Liquids, Transactions of Machine Elements Division / Lund Technical University, Lund, Sweden, Lund Technical University, 1973.
  * [24] R. F. Ausas, M. Jai, G. C. Buscaglia, A mass-conserving algorithm for dynamical lubrication problems with cavitation, J. Tribol. 131 (3) (2009).  doi:10.1115/1.3142903.
  * [25] M. Miraskari, F. Hemmati, A. Jalali, M. Y. Alqaradawi, M. S. Gadala, A robust modification to the universal cavitation algorithm in journal bearings, J. Tribol. 139 (3) (2017).  doi:10.1115/1.4034244.
  * [26] F. J. Profito, M. Giacopini, D. C. Zachariadis, D. Dini, A general finite volume method for the solution of the reynolds lubrication equation with a mass-conserving cavitation model, Tribology Letters 60 (18) (2015) 1–21.  doi:10.1007/s11249-015-0588-0.
  * [27] D. Zhu, J. Wang, Interfacial mechanics: theories and methods for contact and Lubrication, CRC Press, Boca Raton, FL, USA, 2019.
  * [28] W. Habchi, Finite Element Modeling of Elastohydrodynamic Lubrication Problems, Wiley, Hoboken, NJ, USA, 2018.
  * [29] J. H. Ferziger, M. Perić, R. L. Street, Computational Methods for Fluid Dynamics, Springer International Publishing, Cham, Switzerland, 2020.
  * [30] T. J. Chung, Computational fluid dynamics, Cambridge University Press, Cambridge, England, UK, 2002.  doi:10.1017/CBO9780511606205.
  * [31] I. Demirdzic, A finite volume method for computation of flow in complex geometries, Ph.D. thesis, Imperial College London, UK (1982).
  * [32] M. Peric, A finite volume method for the prediction of three-dimensional fluid flow in complex ducts, Ph.D. thesis, Imperial College London, UK (1985).
  * [33] D. Drikakis, W. Rider, High-Resolution Methods for Incompressible and Low-Speed Flows, Springer, Berlin, Germany, 2005.
  * [34] X. Ai, Numerical analyses of elastohydrodynamically lubricated line and point contacts with rough surfaces by using semi-system and multigrid methods (volumes 1 and 2), Ph.D. thesis, Northwestern University, US (1993).
  * [35] Q. J. Wang, L. Sun, X. Zhang, S. Liu, D. Zhu, FFT-Based methods for computational contact mechanics, Frontiers in Mechanical Engineering 6 (2020).  doi:10.3389/fmech.2020.00061.
  * [36] L. A. E. Hough, The stress produced in a semi-infinite solid by pressure on part of the boundary, Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 228 (659-669) (1929) 377–420.  doi:10.1098/rsta.1929.0009.
  * [37] S. Ardah, F. J. Profito, D. Dini, An integrated finite volume framework for thermal elasto-hydrodynamic lubrication, Tribology International 177 (2023) 107935\.  doi:10.1016/j.triboint.2022.107935.
  * [38] S. Liu, Q. Wang, G. Liu, A versatile method of discrete convolution and FFT (DC-FFT) for contact analyses, Wear 243 (1) (2000) 101–111.  doi:10.1016/S0043-1648(00)00427-0.
  * [39] T. Verstraete, S. Scholl, Stability analysis of partitioned methods for predicting conjugate heat transfer, International Journal of Heat and Mass Transfer 101 (2016) 852–869.  doi:10.1016/j.ijheatmasstransfer.2016.05.041.
  * [40] J. Degroote, Partitioned simulation of fluid-structure interaction, Archives of Computational Methods in Engineering 20 (3) (2013) 185–238.  doi:10.1007/s11831-013-9085-5.
  * [41] M. Kaneta, P. Yang, I. Krupka, M. Hartl, Fundamentals of thermal elastohydrodynamic lubrication in si3n4 and steel circular contacts, Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology 229 (8) (2015) 929–939.  doi:10.1177/1350650114565679.
  * [42] M. Kaneta, K. Matsuda, H. Nishikawa, Effects of thermal properties of contact materials and slide-roll ratio in elastohydrodynamic lubrication, Journal of Tribology 144 (6) (2022).  doi:10.1115/1.4053095.
  * [43] N. Ohno, High-pressure behavior of toroidal cvt fluid for automobile, Tribol. Int. 40 (2) (2007) 233–238.  doi:10.1016/j.triboint.2005.09.015.
  * [44] J. Lu, T. Reddyhoff, D. Dini, A study of thermal effects in ehl rheology and friction using infrared microscopy, Tribol. Int. 146 (2020) 106179.  doi:10.1016/j.triboint.2020.106179.
  * [45] M. Hartinger, T. Reddyhoff, Cfd modeling compared to temperature and friction measurements of an ehl line contact, Tribol. Int. 126 (2018) 144–152.  doi:10.1016/j.triboint.2018.05.012.
  * [46] S. Bair, S. Flores-Torres, The viscosity of polyalphaolefins mixtures at high pressure and stress, J. Tribol. 141 (2) (2019).  doi:10.1115/1.4041124.
  * [47] R. Larsson, O. Andersson, Lubricant thermal conductivity and heat capacity under high pressure, Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology 214 (4) (2000) 337–342.  doi:10.1243/1350650001543223.
  * [48] H. Nguyen-Schäfer, J.-P. Schmidt, Tensor Analysis and Elementary Differential Geometry for Physicists and Engineers, Springer, Berlin, Germany, 2017.
  * [49] H. Daiguji, B. R. Shin, Some numerical schemes using curvilinear coordinate grids for incompressible and compressible navier-stokes equations, Sadhana 18 (3) (1993) 431–476.  doi:10.1007/BF02744365.
  * [50] F. White, Viscous Fluid Flow, 3rd Edition, McGraw-Hill, New York, US, 2007.
  * [51] S. Patankar, Numerical Heat Transfer and Fluid Flow, CRC Press, Boca Raton, FL, USA, 1980.