Modelling of mixed lubrication using the heterogenous multiscale methods

Abstract

Surface topography is inherently multiscale. Macroscale topography is covered in microscale surface roughness, which has a substantial, quantifiable impact on overall performance. However, this difference in scales raises modelling challenges. Contemporary mixed lubrication research is dominated by two methods: the accurate but computationally expensive deterministic approach and the more efficient flow factors method that cannot predict local parameters (pressure, contact, etc.) accurately. The principal objective of this work is to develop a framework that will both predict local scale variables accurately and solve in a feasible amount of time.

Through the use of the Heterogeneous Multiscale Methods (HMM) the macroscale and microscale surface features can be modelled using separate domains, whilst maintaining the coupling between the scales. This homogenisation approach allows for both accurate predictions of local scale phenomena and computational efficiency. The HMM have been extensively validated for hydrodynamic/elastohydrodynamic lubrication. Within this work, the HMM are extended to model mixed lubrication and are validated against existing models and experimental data. To build further confidence, plots of simulated data are presented for various conditions on the Stribeck curve.

Both idealised and measured microscale topographies are applied and their impact quantified. The interaction between asperities in contact and lubricant flow has been fully resolved. This framework represents the first time that this microscale behaviour has informed the solution of a computationally efficient macroscale operating within the mixed lubrication regime. Average solution time is 8 h on a local PC with 16GB of RAM.

Introduction

Despite mixed lubrication being one of the three primary lubrication regimes, contemporary mixed lubrication models struggle to provide accurate information on local-scale variables (pressure, contact, etc.) in a computationally feasible amount of time. ¹

Within the hydrodynamic lubrication regime, contacting surfaces are fully separated by the hydrodynamic forces. Conversely, contacting surfaces operating within the mixed lubrication regime are partially separated by the hydrodynamic lubrication forces, but individual asperities will make contact. These asperities will sustain a non-negligible portion of the loading. Regardless of whether a surface operates within the hydrodynamic or mixed lubrication regime, surface topography has a measurable impact on bulk performance. ² The surface topography will be microscale (microns in length), whereas the overall domain will be macroscale (centimetres in length). This separation of scale is arguably the primary modelling difficulty for lubricated contacts.

The average flow approach is one of the most prevalent methods within contemporary lubrication modelling research. This approach applies factors to existing lubrication equations to alter the behaviour and allow non-smooth problems to be modelled. Patir and Cheng ³ published the first average flow model, wherein flow factors were applied to the Reynolds equation, thereby extending it to model the elastohydrodynamic behaviour. Zhu and Cheng ⁴ extended the Patir and Cheng equation to model mixed lubrication by incorporating the Greenwood and Tripp ⁵ asperity contact model. Thereby proving that the parameters that define the microscale roughness have a substantial impact on overall performance.

The work of Zhu and Cheng represented a step forward in the literature but possessed several key flaws: the stochastic asperity distribution is Gaussian, the model is inaccurate for high loading, and the surfaces are assumed to be possess identical roughness. Subsequent authors have taken the average flow approach for mixed lubrication modelling further. For example, Morales-Espejel et al. ⁶ introduced non-Gaussian roughness modelling and Sahlin et al. ⁷ used a homogenisation approach that generates flow factors more representative of real surfaces.

Average flow models are common approaches^3,4,6–8 for mixed lubrication modelling. These models generally solve quickly and predict macroscale variables (friction, loading, etc.) accurately. However, they cannot predict microscale variables with any degree of accuracy (pressure, contact, etc.). This is because the average flow method does not model the microscale and cannot predict how the asperities interact with the flow and deform in contact with the opposing surfaces. Instead, a statistical depiction of roughness is used to inform the macroscale parameters, but no individual microscale features are modelled. Flow factors are just applied to alter equations (typically Reynolds equation) so that they can predict macroscale performance in the mixed lubrication regime. This is a barrier to engineering design, as a database of measurements is required from which to derive flow factors. Typically, these flow factors are derived from experimental data or from microscale simulations.^3,8

Deterministic methods are another major method in contemporary research. A deterministic approach will apply a mesh of sufficient resolution to define every asperity across the macroscale. This meshing strategy results in accurate predictions at both macroscale and microscale levels, but the required number of mesh elements can be prohibitively large. Deterministic methods have been applied to mixed lubrication modelling successfully, but in many scenarios, this approach would lead to an infeasible problem due to the computational expense. Computational Fluid Dynamics (CFD) is a method with similar limitations to the deterministic approach. The mesh resolution required means that CFD is generally applied to single scale applications. ⁹ By utilising the Navier-Stokes equations as a basis, CFD has even lower computational efficiency than deterministic methods using the simpler Reynolds equation.

Jiang et al. ¹⁰ presented one of the first deterministic tribological models. Their predictions of local scale phenomena, like pressure and film thickness, were revolutionary at the time. However, despite serious efforts to increase efficiency, the authors were only able to model a single asperity. A more contemporary example is the work of Liang et al., ¹¹ wherein a water-lubricated thrust pad bearing is modelled deterministically.

Wang et al. ¹² developed a deterministic unified numerical model that accounts for thermal affects, lubricant starvation, and models cavitation across the scales. However, this work is rendered computationally viable by applying an average Reynolds equation over a mesh too coarse to fully resolve the fluid mechanics around asperities. Thus, the local scale phenomenon cannot be said to be fully explored.

Deterministic methods can be deployed at the microscale to calculate flow factors for use in the macroscale. For example, Gulzar et al. use a modified form of the Patir and Cheng model, wherein the flow factors are generated from deterministic microscale models. ¹³ This hybrid approach offsets most of the key disadvantages of the average flow models by actually modelling microscale surface topography. For these reasons, using deterministic microscale models to inform an average flow model via flow factors is an increasingly common technique.^13–15

Homogenisation methods have been proposed to provide both the speed of average flow models, and the accuracy of deterministic approaches. Rom et al. demonstrate that a homogenisation method can achieve an increased level of accuracy compared to Patir and Cheng model's, with the same computational cost. ¹⁶ Patir and Cheng's model is one of the most commonly used models in contemporary tribological research, and yet it demonstrably struggles to accurately model flow dynamics over asymmetric surfaces. ¹⁷ This due to the simplistic boundary conditions and generation of flow factors. Homogenisation can be used to generate accurate flow factors, regardless of the surface topography; as was demonstrated by Persson and Scaraggi. ¹⁸ Using homogenisation to generate flow factors and thereby account for microscale surface topography is an increasingly common technique.^19,20

By employing the Heterogenous Multiscale Methods (HMM), the impact of microscale surface topography can be determined in a feasible amount of time. ²¹ First defined by Weinan and Engquist in 2003, ²² the HMM facilitate the coupling of two disparate domains in a computationally efficient way. In a tribological context, this means a macroscale model that covers the entire solution domain coupled with microscale models that are representations of the microscale surface topography. In this paper, Microscale models are created that will model asperity contact and fluid mechanics around these asperities. Gao and Hewson were the first to apply the HMM to tribological modelling. ²³ Since that first tribological application, the HMM has been steadily expanded in scope by various authors: modelling Elastohydrodynamic Lubrication (EHL) ²⁴ ; adding a metamodel for increased efficiency ²⁵ ; and modelling representative surface topography. ²¹ Throughout all these advances, the core advantage to the HMM has remained its computational efficiency.

Other multiscale methods are applicable to tribology. Brunetière and Francisco develop an excellent multiscale lubrication model that is efficient and accurate, ²⁶ that is subsequently developed to add cavitation. ²⁷ However, their method is highly specialised for modelling parallel plates and likely could not be applied to non-parallel contacts. Multi-grid methods are a multiscale methods that alternate between applied meshes/grids to speed up convergence. These methods are still used in contemporary research, but can be applied to speed up the convergence of direct solvers, and are therefore not relevant here.^28,29

In this paper, a new framework is presented that will allow for accurate predictions of parameters across both microscale and macroscale domain sizes. This framework is computationally inexpensive through the implementation of asperity contact within the HMM. The HMM is far from the first multiscale homogenisation technique, but the efficiency of the methods allows for a much greater separation of scales.^30–32 Lubricant hydrodynamics are modelled with pressure-dependent density and viscosity across both scales, using the Dowson and Higginson density variation equation ³³ and the McEwen viscosity variation equation ³⁴ respectively. Cavitation is modelled via the approach of Söderfjäll et al. ³⁵ The asperity contact is modelled using the Greenwood and Tripp equation ³⁶ and fluid flow around the asperity contacts will cause localised deformation. This deformation is modelled as elastic and the range of input loads is calibrated to keep the stresses below the yield point, this was undertaken to ensure that neither plastic deformation or wear occur in the cases considered as these interaction fall outside the scope of this article. As high load mixed lubrication will lead to plastic deformation, this model cannot be deployed for the entire mixed lubrication regime without the addition of a wear model. Through the combination of these extensively validated approaches, the resulting framework is stable and efficient across the mixed and hydrodynamic lubrication regimes. All nomenclature is presented again in Appendix A.

The framework is successfully validated against models and experimental data from contemporary literature. Further confidence is built through the presentation of results across a Stribeck curve.

Materials and methods

The heterogenous multiscale methods

The Heterogenous Multiscale Methods (HMM) facilitate the separation of scales into disparate modelling domains. Hence the impact of surface roughness on lubricated contact can be determined at an efficient resolution (micrometres) and then homogenised to inform the convergence of the overall solution domain (meters).

The microscale will therefore model the surface features of the pad but not account for the overall morphology; whereas the macroscale will model the overall morphology but not account for the surface features.

Many equations and variables are defined in this section; to aid clarity, the nomenclature is summarised in Appendix A.

Implementation

The flowchart in Figure 1 illustrates the overall data flow and execution order for the modelling framework.

Figure 1.

Flowchart depicting the HMM framework and mathematical procedures deployed in this work.

Multiple variables feature in both the macroscale and microscale. For clarity, macroscale variables are capitalised, and microscale variables are not.

As discussed in the section on macrocscale fluid mechanics, microscale outputs are used to generate gradients that inform the convergence of the textured macroscale variables. Additional microscale variables are passed through to the macroscale via homogenisation, as defined in the homgogenisation section.

Macroscale

The macroscale model covers the entire solution domain. This means that it will always be equal in size to the overall geometry of whatever is being modelled.

Macroscale geometry

The presented framework is applicable to any contacting surface undergoing mixed/hydrodynamic lubrication. Within this work, the framework is applied to a tilted pad bearing made of ASTM A36 ³⁷ steel. This is because previous works that modelled tilted pad bearings validated the HMM extensively.^21,38 Figure 2 depicts the macroscale geometry and the parameters that define the macroscale are presented in Table 1. This pad is assumed to be operating under steady state conditions, and the lubricant is assumed to be isothermal.

Figure 2.

Diagram of the macroscale domain.

Table 1.

Parameters that define the macroscale.

Parameter	Value (Units)
Rotational speed $(\mathit{\omega })$	$1(rad/s)$
Young's modulus $(\mathit{E})$	$200(GPa)$
Yield stress $({\mathit{\sigma }}_{\mathit{y}})$	$250(MPa)$
Poisson ratio $(\mathit{v})$	$0.4$
Stiffness $(\mathit{k})$	$20(MPa/\mu m)$
Asperity height $({\mathit{h}}_{\mathit{A}})$	$3.2(\mu m)$
Inner radius $({\mathit{R}}_0)$	$0.05(m)$
Outer radius $({\mathit{R}}_1)$	$0.1(m)$
Pad angle $({\mathit{\theta }}_{\mathit{d}})$	$50(deg)$
Pad thickness $({\mathit{t}}_{\mathit{b}})$	$10(mm)$

Macroscale solid mechanics

The macroscale model is subject to localised deformation using a linear elastic model based on the infinitesimal strain theory. The governing equations are:

$$\mathrm{\nabla }.\mathit{\sigma }=0$$

(1)

$$\mathit{\sigma }=2Gϵ+\lambda tr(ϵ)\mathit{I}$$

(2)

$$ϵ={\displaystyle \frac{1}{2}[\mathrm{\nabla }\mathit{u}+{(\mathrm{\nabla }\mathit{u})}^T]}$$

(3)

$$\delta =(-\mathit{u}.\mathit{n})$$

(4)

where $\mathrm{\nabla }=\left(\frac{\mathrm{\partial }}{\mathrm{\partial }X},\frac{\mathrm{\partial }}{\mathrm{\partial }Y},\frac{\mathrm{\partial }}{\mathrm{\partial }Z}\right)$ is the gradient operator for the $X,Y,Z$ Cartesian coordinate directions; $\mathit{\sigma }$ is the stress tensor; G is the shear modulus defined by $G=\frac{E}{2(1+v)}$ ; $ϵ$ is the strain tensor; $\lambda$ is Lamé's first parameter, defined by $\lambda =\frac{Ev}{(1+v)(1-2v)}$ ; $\mathit{I}$ is the identity tensor; $\delta$ is the deformation where $z=h$ ; $\mathit{u}=(u,v,w)$ is the vector of deformation and $\mathit{n}$ is the normal vector.

Standard deformation constraints are applied such that the lubricated pad face $(abcd)$ can deform in the normal plane, the bottom pad face $(efgh)$ is fixed in place and cannot deform, and the remaining faces cannot deform normally.

Macroscale fluid mechanics

The key parameters that define the fluid behaviour are presented in Table 2. To facilitate cross regime modelling while only varying load, a hydrostatic jacking force is applied in the form of a hydrostatic pressure $P_{static}$ . This allows loading to vary by several orders of magnitude and a full Stribeck curve to be modelled. It is common for heavy duty bearings to have a pressurised lubricant. This prevents damage or excessive wear during start up by keeping the contacting surfaces separated. ³⁹

Table 2.

Parameters that define the fluid mechanics behaviour of the model.

Parameters	Value (Units)
Unpressurised lubricant density $({\mathit{\rho }}_{\mathbf{0}})$	$850(kg/m^3)$
Unpressurised lubricant viscosity $({\mathit{\eta }}_{\mathbf{0}})$	$0.5(Pa\cdot s)$
Pressure viscosity coefficient $({ϵ}_{\mathbf{0}})$	$16.5(GP{a}^{-1})$
McEwen exponent $(\mathit{q})$	8.38
Cavitation factor $(\mathit{\alpha })$	$0.01$
Max cavitation pressure $(\mathit{\beta })$	$2(MPa)$
Hydrostatic pressure $({\mathit{P}}_{\mathit{s}\mathit{t}\mathit{a}\mathit{t}\mathit{i}\mathit{c}})$	$30(bar)$

Note that $\beta$ is not the pressure at which cavitation begins, but $-\beta$ is the pressure at which cavitation occurs to the maximum possible extent. This is explained further in subsequent sections.

Rotational speed components U and V (in the X and Y coordinate direction respectively) are defined through the introduction of a cylindrical coordinate system. This coordinate system is vertically aligned to the cartesian coordinates such that $r=\sqrt{X^2+Y^2}$ and $\varphi ={\tan }^{-1}\left(\frac{Y}{X}\right)$ . The rotational speed components are then defined as $U=-\omega r\sin (\varphi )$ and $V=\omega r\cos (\varphi )$ .

To facilitate a robust convergence of the model across a wide range of input parameters, a linearised derivative approach is employed. Previous work used a Reynolds basis function with weighting coefficients, however this could not be extended for mixed lubrication. ²¹ This new derivative basis function requires more microscale data but will work across the lubrication regimes.

With the HMM framework, multiple variables are passed from the microscale to the macroscale. These variables, used to build the metamodel and couple the scales, are collected into a vector $\mathit{\zeta }$ with components:

$$\mathit{\zeta }=\left(P,{\displaystyle \frac{\mathrm{\partial }P}{\mathrm{\partial }X},{\displaystyle \frac{\mathrm{\partial }P}{\mathrm{\partial }Y},H,U,V}}\right)$$

(5)

where P is the macroscale hydrodynamic pressure; $\frac{\mathrm{\partial }P}{\mathrm{\partial }X}$ and $\frac{\mathrm{\partial }P}{\mathrm{\partial }Y}$ are the pressure gradients in the x and y direction respectively; H is the macroscale film thickness; and U and V are the respective velocities in the X and Y direction due to the rotation. $P,\frac{\mathrm{\partial }P}{\mathrm{\partial }X},\frac{\mathrm{\partial }P}{\mathrm{\partial }Y},H$ are necessary to incorporate microscale effects into the macroscale. U and V are not unknowns, but they are needed in both scale domains and reduce the metamodeling demands of the microscale. This is because at any given point the combination of U and V are unique. This is also true for $\omega$ and r. In this work, U and V were arbitrarily chosen over $\omega$ and r.

The introduction of ${\zeta }_j$ is primarily to simplify the nomenclature but also serves to demonstrate that the HMM coupling is a mathematical procedure not linked to these specific variables. Previous works have used different ${\zeta }_j$ variables. ²¹ The subscript j is to denote the number of $\mathit{\zeta }$ variables.

A mass conservation approach is used to resolve the mass fluxes across the lubricated region and calculate hydrodynamic pressures. The mass flux $\mathit{Q}$ and its divergence are calculated as

$$\mathit{Q}=f({\zeta }_j),j=1,\dots ,6$$

(6)

$$\mathrm{\nabla }\cdot \mathit{Q}=\sum _{j=1}^6\left({\displaystyle \frac{\mathrm{\partial }{Q}_X}{\mathrm{\partial }{\zeta }_j}\cdot {\displaystyle \frac{\mathrm{\partial }{\zeta }_j}{\mathrm{\partial }X}+{\displaystyle \frac{\mathrm{\partial }{Q}_Y}{\mathrm{\partial }{\zeta }_j}\cdot {\displaystyle \frac{\mathrm{\partial }{\zeta }_j}{\mathrm{\partial }Y}}}}}\right)$$

(7)

where $Q_X$ is the mass flux in the X coordinate direction, and $Q_Y$ is the mass flux in the Y coordinate direction.

The derivatives within Equation (7) are calculated by applying a small perturbation to the ${\zeta }_j$ variables and then using the central difference approximation to determine the gradient. This requires a microscale model with a positive perturbation and a microscale with a negative perturbation. The use of six ${\zeta }_j$ variables require twelve perturbed microscale models and one unperturbed model. Hence thirteen microscale models are needed for every node within the macroscale.

Macroscale film thickness

Macroscale film thickness, H, is calculated using

$$H=H_0+H_g+\delta$$

(8)

where $H_0$ is the contact separation, $H_g$ is a result of the pad geometry, and $\delta$ is the normal mechanical deformation of the pad as defined in Equation (4). $H_g$ is calculated using

$$H_g={\displaystyle \frac{H_{0}r}{R_2}{\displaystyle \frac{\sin (\varphi )-\sin \left({\displaystyle \frac{{\theta }_p}{2}}\right)}{\sin \left(-{\displaystyle \frac{{\theta }_p}{2}}\right)-\sin \left({\displaystyle \frac{{\theta }_p}{2}}\right)}}}$$

(9)

where ${\theta }_p$ is the arc length of the pad in radians and $\varphi ={\tan }^{-1}\left(\frac{Y}{X}\right)$ .

Macroscale density and viscosity variation

Variation in viscosity due to pressure is modelled using the McEwen equation:

$${\eta }_{vary}={\eta }_0{\left(1+{\displaystyle \frac{{ϵ}_0}{\gamma }P}\right)}^{\gamma }$$

(10)

where ${\eta }_{vary}$ is the increased viscosity due to pressure, ${\eta }_0$ is the initial viscosity before pressurisation, ${ϵ}_0$ is the pressure viscosity coefficient, p is the local macroscale hydrodynamic pressure, and $\gamma$ is the McEwen exponent.

Variation in density is modelled using the Dowson and Higginson density equation:

$${\rho }_{vary}={\rho }_0{\displaystyle \frac{(0.59\times {10}^9+1.34P)}{(0.59\times {10}^9+P)}}$$

(11)

where ${\rho }_{vary}$ is the increased density due to pressure and ${\rho }_0$ is the initial density before pressurisation.

Macroscale cavitation

Cavitation is modelled across both the macroscale and microscale model using the mass conserving Söderfjäll et al. ³⁵ cavitation model. There are subtle differences between the implementations, hence the macroscale and microscale equations will be presented separately. This is because the macroscale geometry is purely converging and will not cause cavitation. However, the microscale cavitation effects must be factored into the macroscale solution. Hence, the macroscale cavitation equations are:

$$\rho ={\rho }_{vary}{\displaystyle \frac{F*+\alpha }{1+\alpha }}$$

(12)

$$\eta ={\eta }_{vary}{\displaystyle \frac{F*+\alpha }{1+\alpha }}$$

(13)

where $\rho$ is the density including cavitation effects and the Dowson-Higginson density variation, $\eta$ is viscosity including the cavitation effects and McEwen's viscosity variation, $F*$ is the homogenised cavitation factor (introduced in the homgenisation section), and $\alpha$ is an arbitrary constant. A parametric sweep is performed for $\alpha$ to ensure that convergence is achieved rapidly when $\alpha$ is equal to 0.01.

Macroscale coefficient of friction

It is necessary to calculate a coefficient of friction to demonstrate the transition between the hydrodynamic and mixed lubrication regimes. Separate calculations of friction from hydrodynamic shear and asperity contact forces must be made. The friction forces can then be combined into the coefficient of friction:

$$\mu ={\displaystyle \frac{(T_{hydro}+T_{asp})}{F_0{\displaystyle \frac{R_0+R_1}{2}}}}$$

(14)

where $T_{hydro}$ is the total friction torque due to the hydrodynamic shear, $T_{asp}$ is the total friction torque due to the asperity contact, and $F_0$ is the loading applied to the pad. Within this work, the loading is varied across the framework executions to generate data over a range of Hersey numbers.

Hydrodynamic friction torques are found in the X and Y component directions through integrating the respective shear stresses across the macroscale domain:

$$T_{hydro.X}=\int \int R{\tau }_X*dXdY$$

(15)

$$T_{hydro.Y}=\int \int R{\tau }_Y*dXdY$$

(16)

where $T_{hydro.X}$ and $T_{hydro.Y}$ are the friction torques due to shear in the microscale X and Y coordinate respectively, R is the radius, and ${\tau }_X*$ and ${\tau }_Y*$ are the homogenised microscale shear stresses in the microscale X and Y coordinate respectively. These homogenised variables will be presented in the homgenisation section.

The $T_{hydro.X}$ and $T_{hydro.Y}$ values can then be combined as

$$T_{hydro}=\sqrt{T_{hydro.X}^2+T_{hydro.Y}^2}.$$

(17)

Friction torque due to the asperity contact can be calculated simply using:

$$T_{asp}={\mu }_{material}\int \int R{P}_{asp}*dXdY$$

(18)

where $P_{asp}*$ is the homogenised asperity contact pressure, and ${\mu }_{material}$ is the coefficient of friction inherent to dry contact between the contacting surface asperities.

Metamodel

Within this work, the moving least squares design of experiments approach of de Boer et al. ³⁸ is employed. This approach serves to increase the efficiency of the framework by reducing the number of microscale models required. If the inputs to a microscale model would be sufficiently close to data that already exists within the design of experiments, the output values can be interpolated from the design of experiments data, thereby ensuring a more rapid convergence. The metamodeling mathematics employed within this work are identical to those of Montgomery et al. ²¹ and will not be reintroduced. The full validation of this metamodel can be accessed via the cited paper.

Without this metamodel, every node within the macroscale, will require 13 microscale models for every iteration. Using the metamodel, the outputs can be interpolated where the inputs for a microscale model in a subsequent iteration are suitably close to previous inputs. Therefore, the overall required number of microscale models is reduced. This metamodel represents one of the principal advantages for the HMM compared to other homogenisation methods. The metamodel does not alter the solution overall solution procedure, it merely allows the interpolation of $\mathit{\zeta }$ parameters to find microscale outputs if the input $\mathit{\zeta }$ parameters are sufficiently close to existing data. Without the metamodel, the exact same solution would be achieved, but with an increased computational cost.

Microscale

The relationship between the macroscale and microscale is bounded by the constraint that microscale length l must be significantly smaller than macroscale length L. An order of magnitude difference is the minimum limit. ³⁸

Asperity contact generally leads to the wear of the contacting surface. However, as this model is non-transient, all deformation will be modelled as elastic. This means pressure must be below the yield stress of the material across the entire domain. Figure 3 presents a diagram of the microscale domain and key microscale parameters are presented in Table 3.

Figure 3.

Diagram of the microscale domain.

Table 3.

Key microscale parameters.

Parameter	Value (Units)
Microscale length in x $({\mathit{l}}_{\mathit{x}})$	$1(mm)$
Microscale length in y $({\mathit{l}}_{\mathit{y}})$	$1(mm)$
Asperities per microscale model	$1$
Young's modulus $(\mathit{E})$	$200(GPa)$
Poisson ratio $(\mathit{v})$	$0.4$
Stiffness $(\mathit{k})$	$20(MPa/\mu m)$
Asperity height $({\mathit{h}}_{\mathit{A}})$	$3.2(\mu m)$
Yield stress $({\mathit{\sigma }}_{\mathit{y}})$	$250(MPa)$
Material Friction coefficient $({\mathit{\mu }}_{\mathit{m}\mathit{a}\mathit{t}\mathit{e}\mathit{r}\mathit{i}\mathit{a}\mathit{l}})$	$0.8$

The material properties of the microscale are identical to those of the macroscale, but they are restated in Table 3 for clarity.

Microscale topography

Within this work, two different microscale topographies are applied, both sized for a single asperity. These are an idealised sinusoidal topography and a measured topography. The idealised sinusoidal topography is defined by:

$$h_t={\displaystyle \frac{h_A}{2}\cos \left({\displaystyle \frac{2\pi x}{l_x}}\right)\cos \left({\displaystyle \frac{2\pi y}{l_y}}\right)}$$

(19)

where $h_t$ is the applied idealised sinusoidal topography, $h_A$ is the asperity height, x and y are the coordinates along their respective axis, and $l_x$ and $l_y$ are the microscale lengths in their respective coordinate directions.

Modelling one asperity or many asperities makes no difference to the homogenised outputs, assuming that all the asperities are the same size. The authors found that increasing the microscale length five-fold (leading to twenty-five asperities) had a negligible impact on outputs. The normalised root mean square error between the model with one asperity per microscale and the model that featured an enlarged microscale with twenty-five asperities per microscale was $0.00008$ .

The process for measuring topography is defined subsequently. The raw data is then processed to remove tilt, ensure that any gaps in the capture are filled, and enforce the continuity constraints. Gaps in the capture are filled by assigning a value equal to the average of the adjacent data. The processed topography is then assigned to $h_t$ and implemented into the microscale film thickness.

The overall microscale film thickness equation is

$$h=H+2h_t+{\delta }_{micro}$$

(20)

where h is the microscale film thickness, H is the macroscale film thickness that will define the initial separation, and ${\delta }_{micro}$ is the microscale deformation. Both contacting surfaces are assumed to be rough, with identical and in-phase topography, hence Equation (20) multiplies the applied topography $h_t$ by two. Mirroring the microscale surface topography is not necessary for convergence, but it reduces complexity and makes minimal differences to the output. ${\delta }_{micro}$ is defined as

$${\delta }_{micro}=kp$$

(21)

where k is the material stiffness and p is the microscale hydrodynamic pressure. k is equal to $E/t_b$ . However, if asperity contact is occurring, the value of ${\delta }_{micro}$ will be equal to the overclosure at that point. Overclosure is calculated in Equation (22). The local deformation of the microscale due to asperity contact will necessarily be the minimum value required to fully separate the surfaces.

Asperity contact load

When macroscale loading is sufficient, the film will not fully separate the surfaces and asperity contact will occur. The overclosure of the surfaces is defined as

$$h_{over}=\{\begin{array}{cc}{h}_{min}-h& h\le h_{min}\\ 0& h>h_{min}\end{array}$$

(22)

where $h_{over}$ is the overclosure of the surface, h is microscale film thickness, and $h_{min}$ is the minimum film thickness before asperity contact is said to occur, as defined in the following section. If the overclosure is not zero, there can be no film. At this point h is effectively the separation of the contacting surfaces.

The asperity contact forces are calculated using

$$F_{asp}={\displaystyle \frac{4}{3}E{r}_{contact}^{\frac{1}{2}}{h}_{over}^{\frac{3}{2}}}$$

(23)

Where $F_{asp}$ is the asperity contact load, E is the material's Young's modulus, and $r_{contact}$ is the radius of the asperity contact. The radius $r_{contact}$ is calculated by assuming that the asperity contact region is circular and calculating the radius $r_{contact}$ from the overclosure area. Equation (23) is a reformulation of the Greenwood and Tripp contact equation for Hertzian contact, as presented by Taylor et al. ⁴⁰ Asperity contact pressure $p_{asp}$ can then be simply calculated by dividing the load by asperity contact area $A_{asp}$ such that $p_{asp}=\frac{F_{asp}}{A_{asp}}$ .

Asperity contact flow obstruction

The contact between the asperities will prevent lubricant flow within the contacting region. Modelling this phenomenon requires the implementation of zero flux boundaries within the microscale, such that no fluid can flow into a contacting region. Any region with asperity contact cannot have fluid flow.

This represents a complicated problem to resolve due to the extreme pressures that will occur in ultra-thin films. The piezoviscosity of the lubricant will lead to pressure increasing dramatically as the surfaces approach contact, causing extensive local deformation around the contact. The actual shape of the asperity-on-asperity contact can vary substantially; hence a robust solution mechanism is required.

Within this work, modelling of asperity contact begins by fully separating the surfaces and iteratively moving them towards contact. The hydrodynamic pressures of each iteration inform the deformation of the next. Once the hydrodynamic pressures have been calculated at the minimum film thickness before contact, the localised deformation can be applied to the asperity contact model and solved. This allows for a boundary to be applied that is reflective of the localised deformation.

The assigned value of minimum film thickness ( $h_m$ ) is set such that $h_m=0.5\mu m$ . This is because, below this value, continuum mechanics begins to be invalid, as indicated by the Knudsen number ⁴¹ $K_n=\frac{\lambda }{L_{char}}$ , where $\lambda$ is the mean free path length and $L_{char}$ is the characteristic length. The minimum film thickness is also needed for numerical stability. Negative film thicknesses are filtered out through the implementation of asperity contact boundaries. The minimum film thickness is needed to provide a buffer so that if film thickness reduces during an intermediary hydrodynamic model, it cannot become negative. In this case, $L_{char}$ effectively represents the minimum value of film thickness that can be achieved when using a threshold value of $K_n=0.001$ . ⁴²

Microscale fluid dynamics

The various microscale models use the Reynolds equation to predict lubricant mass flow rates. The microscale equations are:

$$q_x=-{\displaystyle \frac{\rho h^3}{12\eta }{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }x}+{\displaystyle \frac{\rho Uh}{2}}}}$$

(24)

$$q_y=-{\displaystyle \frac{\rho h^3}{12\eta }{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }y}+{\displaystyle \frac{\rho Vh}{2}}}}$$

(25)

$${\displaystyle \frac{\mathrm{\partial }{q}_x}{\mathrm{\partial }x}+{\displaystyle \frac{\mathrm{\partial }{q}_y}{\mathrm{\partial }y}=0}}$$

(26)

where $q_x$ and $q_y$ are the mass flow rates in the x and y direction respectively, $\rho$ is lubricant density, $\eta$ is lubricant viscosity, h is microscale film thickness, and p is the microscale hydrodynamic pressure. Mass conservation is enforced through Equation (26) to calculate hydrodynamic pressures.

Microscale deformation

Microscale deformation is more complicated than the macroscale deformation. For the sake of clarity, a summary will be provided here. Local deformation within the microscale is calculated differently based on the microscale region. In the fluid domain (as can be seen in Figure 4), local deformation is simply a function of the hydrodynamic pressures applied. Within the asperity contact domain, deformation is equal to the overlap between the asperities. This overclosure $h_{over}$ is used to calculate the asperity contact pressure in Equation (23). It is important to note that microscale deformation does not directly couple to macroscale deformation. Instead, the respective homogenised pressures are used to inform macroscale deformation.

Figure 4.

Explanation of the homogenisation boundaries.

Microscale constraints

Four constraints are required for the solution of Equations (24)–(26):

A near-periodicity constraint is applied around the perimeter of the microscale domain:

$$p(l_x,y)=p(0,y)+{\mathrm{\Delta }}_{x}p$$

(27)

$$p(x,l_y)=p(x,0)+{\mathrm{\Delta }}_{y}p$$

(28)

where ${\mathrm{\Delta }}_{x}p$ and ${\mathrm{\Delta }}_{y}p$ are the pressure differences in the corresponding coordinate directions from the macroscale, defined as ${\mathrm{\Delta }}_{x}p=\frac{\mathrm{\partial }P}{\mathrm{\partial }X}{l}_x$ and ${\mathrm{\Delta }}_{y}p=\frac{\mathrm{\partial }P}{\mathrm{\partial }Y}{l}_y$ . This constraint is necessary for the coupling of the scales and ensures that as the microscale length tends towards zero, the microscale tends towards a single point. The applied gradient means that the constraint cannot be described as purely periodic. This condition effectively means that each opposed boundary will have an identical pressure distribution, plus an offset.

• 2. The coupling of the scales necessitates that the average pressure in the microscale is representative of the pressure at the corresponding macroscale node. A constraining pressure $p_0$ is calculated using

  $$p_0=P+{\displaystyle \frac{1}{2}sgn(-{\mathrm{\Delta }}_{x}p)|{\mathrm{\Delta }}_{x}p|+{\displaystyle \frac{1}{2}sgn(-{\mathrm{\Delta }}_{y}p)|{\mathrm{\Delta }}_{y}p|}}$$

  (29)

where $sgn$ is the signum function.

• 3. The application of the second constraint can lead to negative values of pressure on the boundaries. These negative values are non-physical. To prevent this, a minimum required pressure $p_m$ is defined and applied via:

  $$p_m=min(p_0,p_0+{\mathrm{\Delta }}_{x}p,p_0+{\mathrm{\Delta }}_{y}p,p_0+{\mathrm{\Delta }}_{x}p+{\mathrm{\Delta }}_{y}p)$$

  (30)
  $$p(0,0)=\{\begin{array}{c}{p}_o\\ p_0-p_m\end{array}\begin{array}{c}{p}_m\ge 0\\ p_m<0\end{array}$$

  (31)

• 4. The implementation of asperity contact will also lead to microscale flow obstructions, formulated as:

  $$h=0|q_x{n}_x+q_y{n}_y=0.$$

  (32)

where $n_x$ and $n_y$ are normal to the boundary in the x and y component direction.

Microscale density and viscosity variation

Microscale variation in density and viscosity is identical to the macroscale implementation, except that the microscale pressure p is used instead. For clarity, the McEwen viscosity equation and the Dowson-Higginson equation are restated here for the microscale:

$${\eta }_{vary}={\eta }_0{\left(1+{\displaystyle \frac{{ϵ}_0}{\gamma }p}\right)}^{\gamma }$$

(33)

$${\rho }_{vary}={\rho }_0{\displaystyle \frac{(0.59\times {10}^9+1.34p)}{(0.59\times {10}^9+p)}.}$$

(34)

Microscale cavitation

Microscale cavitation modelling is substantially different from the macroscale implementation. The cavitation factor f is calculated using Equation (35) and then used in Equations (36)–(37):

$$f=\{\begin{array}{cc}1& p>0\\ 1-2{\left({\displaystyle \frac{p}{\beta }}\right)}^3-3{\left({\displaystyle \frac{p}{\beta }}\right)}^2& -\beta \le p\le 0\\ 0& p<-\beta \end{array}$$

(35)

$$\rho ={\rho }_{vary}{\displaystyle \frac{f+\alpha }{1+\alpha }}$$

(36)

$$\eta ={\eta }_{vary}{\displaystyle \frac{f+\alpha }{1+\alpha }}$$

(37)

The microscale cavitation factor f is homogenised and used in the macroscale as $F*$ .

Within each model, a parametric sweep for $\alpha$ is performed within COMSOL Multiphysics ⁴³ to ensure that convergence is achieved rapidly when $\alpha$ is equal to 0.01, effectively relaxing $\alpha$ to the target value for this cavitation model. ³⁵

Shear forces

Shear forces are calculated across every microscale model using

$${\tau }_x={\displaystyle \frac{U\eta }{h}\pm {\displaystyle \frac{h^2}{2}{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }x}}}}$$

(38)

$${\tau }_y={\displaystyle \frac{V\eta }{h}\pm {\displaystyle \frac{h^2}{2}{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }y}}}}$$

(39)

where ${\tau }_x$ and ${\tau }_y$ are the shear stresses in the x and y coordinate directions respectively and h is the microscale film thickness. Both equations are arbitrarily formulated to find the shear stress at the opposed contacting surface, where the distance to the surface equals h. When measuring the shear stress at the opposing surface, both $\pm$ should be taken as $+$ . These values for shear are then homogenised for later use.

Homogenisation

Through the homogenisation process, microscale variables are transformed and used as inputs at the macroscale. Through this approach, the impact of microscale surface topography will inform the convergence of the macroscale.

Homogenisation domains

Due to introduction of boundaries within the flow, the homogenisation domain is different for $P*$ and $P_{asp}*$ . Figure 4 presents these differing homogenisation domains.

$\mathrm{\Omega }$ covers the entire microscale domain, ${\mathrm{\Omega }}_s$ is the domain within which asperity contact is occurring, and ${\mathrm{\Omega }}_f$ is the lubricant flow domain. These domains are defined as specified as :

$$\{\begin{array}{c}{\mathrm{\Omega }}_s(x,y)\\ {\mathrm{\Omega }}_f(x,y)\end{array}\begin{array}{c}h(x,y)\le h_{min}\\ h(x,y)>h_{min}\end{array}$$

(40)

$$\mathrm{\Omega }(x,y)={\mathrm{\Omega }}_s(x,y)\cup {\mathrm{\Omega }}_f(x,y)$$

(41)

$${\mathrm{\Gamma }}_x={\mathrm{\Omega }}_f\left(\pm {\displaystyle \frac{l_x}{2},y}\right)$$

(42)

$${\mathrm{\Gamma }}_y={\mathrm{\Omega }}_f\left(x,\pm {\displaystyle \frac{l_y}{2}}\right)$$

(43)

The $\mathrm{\Gamma }$ fixed fluid flux boundaries specified in Equations (42) and (43) simply enforce that the flow rate into the microscale, is equal to the flow rate out of the microscale for each coordinate direction. Hence choosing a side is arbitrary and the value of ${\mathrm{\Gamma }}_x$ and ${\mathrm{\Gamma }}_y$ does not change.

Coupling the domains

Coupling between the scales is achieved via the following equations:

$$P*={\displaystyle \frac{\int \int pd{\mathrm{\Omega }}_f}{\int \int d{\mathrm{\Omega }}_f}}$$

(44)

$$P_{asp}*={\displaystyle \frac{\int \int p_{asp}d{\mathrm{\Omega }}_s}{\int \int d{\mathrm{\Omega }}_s}}$$

(45)

$$F*={\displaystyle \frac{\int \int fd{\mathrm{\Omega }}_f}{\int \int d{\mathrm{\Omega }}_f}}$$

(46)

$$Q_X={\displaystyle \frac{\int q_{x}d{\mathrm{\Gamma }}_x}{\int d{\mathrm{\Gamma }}_x}}$$

(47)

$$Q_Y={\displaystyle \frac{\int q_{y}d{\mathrm{\Gamma }}_y}{\int d{\mathrm{\Gamma }}_y}}$$

(48)

$${\tau }_X*={\displaystyle \frac{\int \int {\tau }_{x}d{\mathrm{\Omega }}_f}{\int \int d{\mathrm{\Omega }}_f}}$$

(49)

$${\tau }_Y*={\displaystyle \frac{\int \int {\tau }_{y}d{\mathrm{\Omega }}_f}{\int \int d{\mathrm{\Omega }}_f}}$$

(50)

where all variables are as defined in Table 4. These homogenised variables are used to inform macroscale convergence with the microscale outputs, as presented in Equation (7). Through this method, macroscale convergence is informed by microscale surface topography. ${\tau }_X*$ and ${\tau }_Y*$ are only used to calculate a coefficient of friction, as demonstrated in Equations (15) and (16). The structure of these equations ensures that as the microscale model tends towards a point, the overall solution would recover the smooth Reynold's equation solution, thereby mathematically justifying the approach.

Table 4.

The homogenised variables passed between scales.

Microscale variable	Macroscale homogenised equivalent
Hydrodynamic pressure $(\mathit{p})$	Homogenised hydrodynamic pressure $(P*)$
Asperity contact pressure $({\mathit{p}}_{\mathit{a}\mathit{s}\mathit{p}})$	Homogenised asperity contact pressure $(P_{asp}*)$
Cavitation factor $(\mathit{f})$	Homogenised cavitation factor $(F*)$
Microscale flow rate in x $({\mathit{q}}_{\mathit{x}})$	Macroscale flow rate in X $(Q_X)$
Microscale flow rate in y $({\mathit{q}}_{\mathit{y}})$	Macroscale flow rate in Y $(Q_Y)$
Microscale shear stress in x $({\mathit{\tau }}_{\mathit{x}})$	Macroscale shear stress in X $({\tau }_X*)$
Microscale shear stress in y $({\mathit{\tau }}_{\mathit{y}})$	Macroscale shear stress in Y $({\tau }_Y*)$

Measured surface topography

A comparison of the outputs for measured and idealised asperity topography is presented in Figure 5. The measured topography was captured from a tilted pad bearing using the Bruker NPFLEX optical profilometer. ⁴⁴ The captured topography data was then postprocessed to eliminate noise and ensure that any gaps in the capture were addressed. This asperity data can then be applied to the microscale and contrasted with the output data when the idealised topography is applied.

Computation

The models for this work were solved using the finite element method in COMSOL Multiphysics version 5.5 with MATLAB Livelink. ⁴³ All models were executed on a personal computer with 16GB of RAM and an 11^th Gen Intel Core i5-11400 processor. The PARDISO solver was used within COMSOL Multiphysics. ⁴⁵ Convergence was achieved at a solver error tolerance of ${10}^{-3}$ .

The macroscale mesh employed in this work contains 1606 elements. This mesh is generated using the triangular mesh function within COMSOL. Increasing the number of elements to 6923 has a negligible (less than 1%) impact on outputs, but results in a significant increase in computational cost. Conversely reducing the mesh density will make convergence less reliable and generally increase the number of iterations required to achieve convergence. Hence any potential gain in efficiency is offset by the increased number of iterations.

The meshing strategy employed for the microscale model also features a triangular mesh applied within COMSOL. It is important to note that the macroscale and microscale mesh are wholly separate. However, the mesh is also configured so that if additional boundaries are introduced due to asperity contact, the mesh resolution will increase around that boundary. This will result in most microscale models differing in element count.

For an exemplar microscale model with significant asperity contact, applying the meshing strategy results in a mesh of 245,420 elements. Applying a mesh of 979,692 elements changes the asperity contact area by 0.63%. Conversely, reducing the mesh resolution to 61,604 results in an asperity contact area change of 1.18%. Thus, the meshing strategy employed is appropriate. A mesh this fine is necessary to fully resolve flow boundaries in every microscale model.

Results

This section contains a comparison of idealised versus measured topography; exemplar models operating across the Stribeck curve, and a summary of the validation undertaken.

Idealised versus measured

The applied microscale topography is presented as Figure 5. Both topographies are required to be symmetrical around the midline of the x and y axes to ensure that the continuity constraint is enforced. Both topographies have identical peak-to-peak heights and average heights.

Figure 5.

Plots of the applied idealised (a) and measured (b) topography.

Microscale performance

Figure 6 presents the microscale hydrodynamic and asperity contact pressures with both idealised and measured topography, thus allowing for a direct comparison. Both microscale models are from the corner denoted as c in Figure 2 on the macroscale and have the same inputs. Plots of cavitation factor $f$ are provided as Figure 6(e) and (f) to demonstrate the effect topography has on micro-cavitation.

Figure 6.

A comparison of microscale outputs for idealised and measured topography. (a) and (b) are hydrodynamic pressures for idealised and measured topography respectively. (c) and (d) are asperity contact pressures for idealised and measured topography respectively. (e) and (f) are the cavitation factor f for the idealised and measured topography respectively.

Stribeck curve

19 different models were run with varying input loads for both idealised and measured microscale topographies. Figure 7(a) plots the coefficient of friction against Hersey number for these models. The Reynolds basis model of Montgomery et al. ²¹ is used to predict coefficients of friction, before transition to mixed lubrication, for comparison. In Figure 7(a), mixed lubrication is defined as occurring once the asperity contact pressure supports more than 1% of the load. Across the Stribeck curve of Figure 7(a), the average percentage difference in coefficient of friction between idealised and measured topography is $4.59\text{\%}$ . This number rises to $7.9\text{\%}$ across the transition region where asperity contact begins to dominate hydrodynamics. Figure 7(b) plots the $\lambda$ ratio against the Hersey number. The $\lambda$ ratio is commonly used to estimate the lubrication regime. It is calculated as film thickness over the asperity height. Whilst the values of $\lambda$ that denote lubrication regime transition will vary depending on the surface and loading scenario, it can safely be said that for boundary lubrication, $\lambda \le 1$ , and for hydrodynamic lubrication, $\lambda >5$ . ⁴²

Figure 7.

(a) is the logarithmic plot of the coefficient of friction against Hersey number.

Lubrication regime outputs

To build confidence in the modelling techniques, models from three different regions of the Stribeck curve are presented, as shown in Table 5. These models have idealised sinusoidal topography and were chosen to ensure that the presented models are evenly spread across the Stribeck curve.

Table 5.

Inputs to the exemplar models.

Lubrication regime	Hersey number	$\mathit{\lambda }$
Hydrodynamic	$3.48\times {10}^{-6}$	$19.31$
Mixed (Transition)	$1.76\times {10}^{-6}$	$3.58$
Mixed	$3.70\times {10}^{-7}$	$3.44$

Exemplar microscale models are presented alongside the macroscale models in subsequent sections. These exemplar microscale models are all in the same geometric location on the pad. This point is on the corner denoted as c in Figure 2. For all the models in this paper, the local stresses are below the yield stresses. This means that all deformation is elastic, the surfaces are not plastically deformed, and wear does not occur. ⁴⁶ This was confirmed by creating a 3-dimensional microscale model and determining that the sub surface stresses peak at $39.1MPa$ for the mixed lubrication model presented in Table 5. This is substantially below the yield stress of $250MPa$ . The method was the same as that of de Boer et al. ⁴⁷ but with the addition of contact mechanics. As this model features the highest macroscale loading and the highest microscale pressure, it is clear that all stresses in the models are below the yield limit and therefore neither plastic deformation or wear will take place.

Hydrodynamic

Figure 8 presents the performance of the model when operating under the hydrodynamic lubrication regime. As there is no asperity contact under these loading conditions, the plots of asperity contact pressure are not included.

Figure 8.

Key outputs when operating within the hydrodynamic lubrication regime, (a) and (b) are the hydrodynamic pressure across the macroscale and an exemplar microscale model and (c) and (d) are the film thicknesses across the macroscale and microscale. (b) and (d) depict the same microscale model, located at point c in Figure 3.

Mixed (transition)

Figure 9 presents the model performance when the contacting surfaces are operating within the mixed lubrication regime, but asperity contact is still relatively minor (but non negligible). The blank spaces at the centre of Figure 9(b) and (d) are because asperity contact physically blocks the flow, hence there can be no hydrodynamic pressure or film thickness.

Figure 9.

Key outputs when operating within the transition between the hydrodynamic and mixed lubrication regimes, (a) and (b) are the hydrodynamic pressure across the macroscale and an exemplar microscale model, (c) and (d) are the film thicknesses across the macroscale and a microscale model, and (e) and (f) depict the asperity contact pressure across the macroscale and a microscale model. (b), (d), and (f) all depict the same microscale model, located at point c in Figure 3.

Mixed lubrication

Figure 10 presents the model performance when the contacting surfaces are operating deep within the mixed lubrication regime.

Figure 10.

Key outputs when operating within the mixed lubrication regime, (a) and (b) are the hydrodynamic pressure across the macroscale and an exemplar microscale model, (c) and (d) are the film thicknesses across the macroscale and a microscale model, and (e) and (f) depict the asperity contact pressure across the macroscale and a microscale model. (b), (d), and (f) all depict the same microscale model, located at point c in Figure 3.

Computational cost

Table 6 lists the execution times of the presented HMM framework for each of the exemplar models.

Table 6.

Execution time of the models.

Lubrication regime	Hersey number	Execution time (hours)
Hydrodynamic	$3.48\times {10}^{-6}$	$5.46$
Mixed (Transition)	$1.76\times {10}^{-6}$	$5.90$
Mixed	$3.70\times {10}^{-7}$	$6.36$

Comparison to contemporary literature

In previous works, the HMM framework was validated against deterministic models. ²¹ However, the addition of asperity contact within this framework has rendered any attempt at a deterministic comparison infeasible. When attempted, the deterministic model required over $100GB$ of rapid access memory to just assemble the matrices, even after reducing the frequency of the microscale roughness. The required computational resources were not available to solve even this simplified deterministic model. Instead, a comparison will be made with multiple contemporary models and experimental data.

Montgomery et al.

The idealised sinusoidal topography results from Montgomery et al. ²¹ are replicated using this mathematical framework. This required that some features of the new framework be disabled (piezoviscosity and compressibility). Between both models, there is a mean percentage error of 0.1% and 4.8% for pressure and film thickness respectively. The root mean squared error for pressure and film thickness is $0.0292$ and $0.0599$ respectively. These results prove that the data of Montgomery et al. ²¹ has been reproduced.

Liang et al.

To further validate the framework, a model is constructed for direct comparison to the published works of Liang et al. ⁴⁸ This is a model of a water-lubricated thrust bearing. No changes are made to the framework, beyond adjusting the geometry of the pads, the material properties, and the fluid properties. Figure 11 presents a comparison between this framework and the work of Liang et al.

Figure 11.

A comparison of the measured and calculated film thickness of Liang et al. and the outputs predicted by the framework presented in this paper. (a) is at 150 kPa and 200 rpm, (b) is at 150 kPa and 600 rpm.

Table 7 compares the accuracy of this framework in predicting the experimental measurements of Liang et al., to the accuracy of the original predictions by Liang et al. ⁴⁸ This comparison is performed using a calculation for root mean squared error on the data presented in Figure 11.

Table 7.

A quantitative comparison between the predictions of the framework presented in this work, and the original model of Liang et al.

Film thickness predictions by	Root mean squared error at 200 RPM	Root mean squared error at 600 RPM
This framework	$0.0383$	$0.0230$
Liang et al. model	$0.112$	$0.0786$

The root mean squared errors are lower for the predictions of this framework. Thereby quantitively demonstrating an increased modelling accuracy.

Discussion

The capacity for the presented framework to model the hydrodynamic and mixed lubrication regime has been clearly demonstrated and validated. This represents the first time the HMM have been used to model mixed lubrication. The framework is both accurate at the local scale and computationally efficient.

Earlier lubricated contact models with the HMM used a Reynolds equation basis function for macroscale fluid dynamics. ³⁸ By moving to a linearised derivative-based approach, in this paper, lubricated contacts can be modelled no matter how far away they are from the Reynolds equation. This innovation is essential for implementing asperity contact. Moving to the linearised derivative basis does have one noticeable drawback: the required microscale data is increased by a factor of thirteen, thereby increasing the computational cost. This increased cost is partially offset as the framework converges significantly faster, typically only requiring two to three iterations. However, the execution time is still an average of two times longer than the Reynolds basis HMM formulation. ²¹ This increase in computational cost is offset by the ability to model mixed lubrication and lubricated contact scenarios that differ significantly from the Reynolds equation.

The computational cost of this framework is dramatically less than a deterministic approach. Exactly by how much is hard to define, because many authors do not provide data on execution time and computation hardware.^{11,14,49–51} The best point of comparison may well be that the authors of this paper attempted to solve a deterministic model to provide a comparison. However, the model could not be solved as it required over $100GB$ of RAM just to assemble the matrices.

Figure 7 presents the Stribeck curve for this framework and for the framework of Montgomery et al. ²¹ The advantages of the new framework are clear as the Reynolds basis framework cannot solve if the solution moves too far away from the Reynolds equation. This means that the Reynolds basis framework fails before mixed lubrication occurs. While each model was provided the same inputs for loading, neither of the Reynolds basis models could solve at a Hersey number below $2.5\times {10}^{-6}$ . As this threshold is approached, the Reynolds basis models increasingly diverge from the predictions of the new, derivative based framework. This further proves the advantages in stability and accuracy of prediction for the new framework. Figure 8 demonstrates the lubrication regime transition and provides additional evidence for the regime transition indicated in the Stribeck curve of Figure 7. Although the transition does not occur exactly at the value of $\lambda =5$ suggested by the literature as an approximation, the transition is nonetheless evident. The coefficient of friction increase can therefore be explained by the onset of mixed lubrication.

The HMM framework itself is already extensively validated for modelling lubricated contact.^21,23,38 By comparing the developed framework against the work of Montgomery et al., it is proven that the accuracy of the HMM is not compromised by forsaking the Reynolds basis for a derivative basis. The mixed lubrication implementation is validated against the experimental results of Liang et al. ⁴⁸ Confidence is further built in the model through the presentation of a wide array of outputs within various stages of the Stribeck curve.

Loading input conditions are selected to ensure that the yield stress of the pad material was not exceeded. This ensured that all deformation was elastic, and wear did not occur. The implementation of plastic deformation and wear is an excellent candidate for future work. This would allow the framework to be extended into the boundary lubrication regime. Therefore requiring the modelling of permanent deformation and volume removal of surface asperities at the micro-scale to be fully coupled into the macro-scale using homogenisation and the development of the corresponding computational algorithm.

Within this work, the HMM framework is successfully deployed to model mixed lubrication. Perhaps most significantly, the resulting framework is efficient and stable regardless of input conditions. This framework represents the first time that pressure variant viscosity and density have been implemented within the HMM. Critically, the key advantages associated with the HMM are maintained, namely computational efficiency and the ability to couple disparate domains. For these reasons, this framework is therefore beyond the capabilities of contemporary methods to model the hydrodynamic and mixed lubrication regimes. Future work will target a reduction in the amount of required microscale data. The presented framework requires an average of 7 times more microscale models compared to Montgomery et al. ²¹ If a similar data requirement for previous hydrodynamic models using the HMM framework can be matched for mixed lubrication, computational cost can be reduced by nearly an order of magnitude.

Conclusions

In this work, a combined hydrodynamic and mixed lubrication solution framework is presented. Mixed lubrication can be modelled as long as the stresses are within the elastic limit, because plastic deformation and wear are not considered. This framework leverages the efficiency of the HMM to model the effect of microscale surface roughness on overall performance, in a computationally viable amount of time. On a local PC of average specifications, average execution time averages to 8 h across the Stribeck curve. This represents a huge step forward, as the performance of macroscale lubricated contact, within the mixed lubrication regime, can be predicted using realistic microscale modelling, at a reasonable computational expense.

The need for modelling microscale roughness is clearly demonstrated through comparison of the macroscale performance with idealised or measured microscale topography. These topographies are identical in peak to trough distance, average asperity height, and asperity frequency, but still result in substantial differences in bulk property predictions. Across the Stribeck curve, coefficients of friction varied by an average of 5% between the idealised and measured topographies, rising to 9.57% as asperity contact begins to dominate over hydrodynamics. This highlights that accurate macroscale predictions require microscale considerations.

The presented framework is validated against contemporary literature and can model cavitation, piezoviscosity, fluid compressibility, asperity contact, flow obstruction, and local deformation.