# Full and Hybrid Multiscale Lubrication Modeling

Arthur Francisco, Noël Brunetière

► To cite this version:

Arthur Francisco, Noël Brunetière. Full and Hybrid Multiscale Lubrication Modeling. Lubricants, 2022, 10 (12), pp.329. 10.3390/lubricants10120329. hal-04223334

HAL Id: hal-04223334

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

Submitted on 29 Sep 2023

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

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

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

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

Article

# Full and Hybrid Multiscale Lubrication Modeling

Arthur Francisco <fnref n="1" />,*,† and Noël Brunetière <fnref n="2" />,†

**Abstract:** The numerical solution for the lubrication of parallel rough surfaces cannot be obtained using the well-known flow factors of Patir and Cheng. Nor can it be determined using homogenization techniques. Is there an alternative, besides a purely long-term deterministic way of solving the problem? The present paper aims at proposing a multiscale approach in order to reduce the computing time, specific to deterministic resolutions, while maintaining good accuracy. The configuration is a parallel rough surface slider, with imposed hydrodynamic operating conditions. The domain consists of independent macro-elements, on which the Reynolds equation is solved. Then, the macro-element boundaries are adjusted to ensure global mass conservation. In its hybrid version, the algorithm replaces some well-chosen macro-elements by simple linear finite elements. The results clearly show the potential of our method. Because the lubrication of each macro-element can be processed independently, the multicore architecture of the processor is exploited. Even if the performance depends on the ratio roughness/height, the computing time is half than for the classical deterministic method, with a few percent errors. The work concludes with some recommendations on the configurations for which the multiscale method is best suited, such as surfaces with short correlation lengths.

**Keywords:** lubrication; multiscale modeling; rough surfaces

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

check for
updates

**Citation:** Francisco, A.; Brunetière, N. Full and Hybrid Multiscale Lubrication Modeling. *Lubricants* **2022**, *10*, 329. https://doi.org/10.3390/lubricants10120329

Received: 27 October 2022

Accepted: 19 November 2022

Published: 23 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

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

## 1. Introduction

In a lubricated contact, when the fluid enters a converging zone, the pressure increases and reaches a maximum near the minimum of the film thickness. Assuming that a diverging area follows the converging entrance, the fluid is stretched and the pressure decreases. The Reynolds equation as originally expressed does not account for the cavitation area in a realistic way and abundant literature has followed. Among the methods that guarantee the Reynolds conditions, the penalty method and the well-known Elrod's algorithm are quite easy to implement. The penalty method promoted by Wu [1] treats both the film area and the cavitation area by introducing a penalty term into the Reynolds equation. When the pressure falls below the cavitation pressure $P_c$, its value is forced to $P_c$. Whilst providing an accurate pressure boundary, the method fails in predicting the film reformation boundary. The seminal Elrod and Adams algorithm [2,3] splits the Reynolds equation by the means of a switch function to ensure the equation validity in both the gaseous and liquid areas. Mass conserving is guaranteed in both areas; however, the location of the boundaries is mesh dependent.

Let us describe the cavitation phenomenon a little more. If the pressure falls below the ambient pressure, three situations can be considered [4]: (a) the dissolved gas is released and forms bubbles, (b) the bubbles already present expand, and (c) the fluid, which contains no gas, evaporates. The cases are referred to as a (a) gaseous cavitation, (b) pseudocavitation, and (c) vaporous cavitation. In a more recent study, Bai et al. [5] propose a review on the cavitation in a thin liquid layer, in which a section is dedicated to the hydrodynamic cavitation. However, the authors do not distinguish the different cavitation phenomena as reported by Braun and Hannon [4]. The simplified taxonomy (a)-(b)-(c) is that the

<footnote n="1">¹ Institut Pprime, CNRS, Université de Poitiers, ISAE-ENSMA, 16021 Angoulême, France</footnote>

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

<footnote>* Correspondence: arthur.francisco@univ-poitiers.fr; Tel.: +33-545-673-249</footnote>

<footnote>† These authors contributed equally to this work.</footnote>
latter present does not include adsorption effects, as studied by Belova et al. [6], on the heterogeneous cavitation. Although rough surfaces are likely to adsorb much more gas than smooth surfaces because of their fractal nature, the related phenomenon significantly complicates the lubrication modeling of rough surfaces and therefore is not addressed in the present work.

Regarding the gaseous cavitation in oil, Li et al. [7] and Song and Gu [8] consider the lubricant as a mixture of pure liquid and air, as the result of the dissolution and the release of air in an instantaneous equilibrium saturation state. The model satisfies the classical Jakobson–Floberg–Olsson (JFO) conditions, but it does not account for the gaseous cavitation rate, as performed by Hao and Gu [9] and more recently by Ding et al. [10].

Ransegnola et al. [11] model both the vaporous and gaseous cavitation in an oil-lubricated bearing, predicting not only the cavitation area but also the distribution of the gas, vapor, and liquid states. However, with respect to water as a lubricant, Magaletti et al. [12] provide a graph that shows the cavitation pressure of ultra-pure water as a function of temperature. The remarkable data are the cavitation pressure of the water at ambient conditions, about –120 MPa, because of its high tensile stress. This very low value makes one think that behind the common cavitation in water, there is essentially a gaseous cavitation or even a pseudocavitation.

In the present work, the fluid is assumed to be a biphasic mixture: air bubbles are present in an incompressible liquid, which leads to the pseudocavitation in the diverging zones. Considering the mixture as a homogeneous media allows for the use of the Reynolds equation in the whole contact zone: the fluid rheology is modified in the depressurized zones and is constant in the full-film zones. Three advantages are brought with this model presented by Brunetière [13] and referred to as the “Lubricant General Model” (LGM). (1) It is a handy model because once the fluid density and the fluid viscosity are defined, the Reynolds equation is solved in the same way, whatever the zone. (2) The cavitation area transition is smooth, making it particularly suitable for multiscale meshes. (3) The varying density model is a mass conserving one; as such, it accurately accounts for the film rupture and reformation.

Grützmacher et al. [14] propose a review of the multiscale approaches about texturing in tribology. The classification that is proposed is well suited to distinguish the different strategies related to the lubrication of rough contacts. In particular, two classes are of importance here: the analytical multiscale methods and the numerical multiscale modeling.

Lubricated contacts may require multiscale approaches for various reasons. Heavily loaded contacts can lead to high gradient pressures, such as in the outlet spike in point contacts. The problem can be addressed with fine meshes, but without speedup convergence algorithms, such as multigrid methods, the computing time becomes prohibitive. The multiscale approach is then a means to accelerate the convergence of the Reynolds equation resolution. The mesh can also exhibit different element scales based on the pressure gradient: the steeper the pressure field, the finer the mesh. In both situations, the Reynolds equation is not modified. The above common techniques belong to the numerical multiscale modeling class.

The analytical multiscale methods involve treatments applied to the Reynolds equation. When roughness influences the fluid flow because of the small film thicknesses, instead of a deterministic resolution, the equation can be enriched with flow factor modifiers [15,16]. The flow factors are determined with a few of the statistical properties of the rough surface, which limit the accuracy of the results. However, once the flow factors are determined, the surfaces being considered smooth, coarse meshes can be used, decreasing the computing effort. This stochastic approach has since been improved to cope with the micro-cavitation to handle a broader variety of surface roughnesses [17], but it remains global: the local effects of the roughness on the pressure cannot be captured. In order to get rid of the stochastic approximation of the roughness, the Reynolds equation can be viewed as a set of equations solving different wavelength pressure problems [18–20]. Indeed, an asymptotic expansion of the pressure is written with respect to a scalar related to the roughness
wavelength, leading to a modified Reynolds equation. The process, based on rigorous mathematical developments, is called homogenization and is compatible with any periodic roughness. Homogenization techniques have received more attention in recent years than flow factor modifiers with the works of, among others, de Boer et al. [21,22], de Boer and Almqvist [23], and Han et al. [24]. In Rom et al. [25], the reader can find the homogenization advantages that explain its greater development. However, flow factor modifiers seem to be more widely used to date, such as in [26–30].

Computing flow factors and homogenization suffer from a common drawback: for flat parallel rough surfaces, no pressure build-up exists, although it is experimentally observed [31–34]. A compromise between the deterministic and the stochastic methods has been proposed by Brunetière and Wang [35]. The Reynolds equation is filtered: above an $f_c$ roughness frequency, averaging is used; and below $f_c$, a deterministic solution is computed on a coarse grid. For this analytical multiscale method, the computational effort is much less than for the deterministic case, but the micro-cavitation is not taken into account. Pei et al. [36,37] use a Guyan reduction to condensate finite element cells. In doing so, the general linear system is smaller than the one obtained with the usual finite element method. A further reduction is obtained while defining the master/slave nodes on the cell boundary. Even if the bandwidth is larger, the computing time is up to five times smaller than for the conventional finite element method (FEM). The major drawback of this numerical multiscale modeling is, however, that the cavitation is not taken into account.

The present paper, belonging to the numerical multiscale approaches, is a step beyond the work of Brunetière and Francisco [38]: the domain is divided into macro-cells inside which a deterministic FEM is used to solve the Reynolds equation; then, the macro-cells are linked together using a mass-conserving principle to cover the whole domain. In addition, when the average film thickness over a macro-element is sufficient to ignore the roughness, the macro-element is not submeshed but rather replaced by a linear finite element: it is the hybrid version of the algorithm.

## 2. Materials and Methods

A flat lubricated rough contact is modeled with the functioning parameters presented in Table 1.

<table><caption>**Table 1.** Contact characteristics.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Minimum of the film thickness</td><td>0.05 to 10.0 µm</td></tr><tr><td>Fluid type</td><td>mixture (liquid and gas)</td></tr><tr><td>Top-scale mesh, N<sub>e</sub><sub>x</sub> × N<sub>e</sub><sub>y</sub> elements</td><td>128 × 128</td></tr><tr><td>Bottom-scale meshes, n<sub>e</sub><sub>x</sub> × n<sub>e</sub><sub>y</sub> elements</td><td>8 × 8</td></tr><tr><td>Sliding speed along x</td><td>5.0 m s<sup>-1</sup></td></tr><tr><td>Fluid viscosity</td><td>10<sup>-3</sup> Pa s</td></tr><tr><td>Fluid density</td><td>10<sup>3</sup> kg m<sup>-3</sup></td></tr><tr><td>Ambient pressure</td><td>10<sup>5</sup> Pa</td></tr></tbody></table>

The rough surface used in the simulation is a numerically generated 1025 × 1025 points bi-Gaussian surface, which is typical of surfaces used in tribological applications, Figure 1. It is half the one used by Brunetière and Francisco [38]. The surface characteristics are presented in Table 2.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Bi-Gaussian rough surface used for the numerical simulations.</figcaption></figure>

<table><caption>**Table 2.** Rough surface characteristics.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Surface type</td><td>numerically generated</td></tr><tr><td>Numerical size n × n points (nodes)</td><td>1025 × 1025</td></tr><tr><td>Physical size, L<sub>x</sub> × L<sub>y</sub></td><td>1 cm × 1 cm</td></tr><tr><td>RMS roughness, Sq</td><td>0.1 µm</td></tr><tr><td>Roughness skewness, SSk</td><td>-1.4</td></tr><tr><td>Roughness kurtosis, SKu</td><td>6.3</td></tr><tr><td>Autocorrelation length, Sal(s = 0.2)</td><td>124 µm</td></tr><tr><td>Ratio of domain length L<sub>x</sub> to correlation length Sal, Rcl</td><td>80.4</td></tr></tbody></table>

## 2.1. *Model Equations*

Under the usual lubrication assumptions, the Reynolds equation writes:

$$ \nabla \left( \frac{\rho h^3}{\mu} \nabla p \right) = \nabla (6\mu\rho h) \qquad (1) $$

where the left-hand side is the diffusive part of the equation, and the right-hand side is the convective part if $\rho$ is pressure-dependent.

As fully detailed in the supplementary material (SM), Section S1, the Bubnov-Galerkin method leads to spurious oscillations if the convective part becomes greater than the diffusive part. As a solution, whereas the Poiseuille term is weighted by the shape functions $N_i$, the Couette term is weighted by upwind functions $\tilde{N}_i$ defined as:

$$ \tilde{N}_i = N_i + \frac{\alpha}{2} l_e \frac{\partial N_i}{\partial x} \qquad (2) $$

where

* $l_e$ is the streamline length of an element;

* $\alpha = \coth(P_e) - \frac{1}{P_e}$, $P_e$ is the Peclet number of an element;

* $P_e = 3 \frac{\mu u (\partial \rho / \partial p)}{\rho h^2} l_e$.
It is the Petrov-Galerkin scheme for which different functions than the shape functions are used for the weak formulation. Here, the shape and the weighting functions are chosen linear.

Hence, after summation over the entire domain $\Omega$, one obtains the weak form of the Equation (1) for each internal node $i$:

$$ -\iint_{\Omega_i} \frac{\rho h^3}{\mu} \nabla N_i \cdot \nabla p \, d\Omega + \iint_{\Omega_i} \overline{\rho h} \mathbf{u} \cdot \nabla N_i \, d\Omega = 0 \quad (3) $$

where

* $\Omega_i$ is the subdomain where $N_i$ is not null;

* $\overline{\rho h} = (\rho h)_k \bar{N}_k$;

* $\bar{N}_k = N_k - \frac{\alpha}{2} l_e \frac{\partial N_k}{\partial x}$.

As explained in the SM, Subsection S2.1 and S2.2.

During the resolution process, the left part of Equation (3) is actually not null but equal to a residual $R_i$. The Newton-Raphson method leads to the determination of the increment $\delta p_j$ at the node $j$, such that $R_i(p + \delta p_j) = 0$, hence:

$$ R_i + \frac{\partial R_i}{\partial p_j} \delta p_j = 0 \quad (4) $$

Thus, for all $\Omega$ internal nodes, the pressure increments are iteratively determined after the resolution of the system $[K]\{\delta p\} = -\{R\}$, detailed in SM, Subsection S2.3.

The implementation of the lubricant general model is straightforward because it is a mixture of an incompressible liquid—subscript “l”—and a gas—subscript “g”. Noting $r_g$, the specific gas constant, $\lambda$ the gas mass fraction, and $T$ the ambient temperature:

* For an incompressible liquid:

$$ \begin{aligned} \rho &= \rho_l \\ \frac{\partial \rho}{\partial p} &= 0 \end{aligned} $$

* For a perfect gas:

$$ \begin{aligned} \rho &= \frac{p}{r_g T} \\ \frac{\partial \rho}{\partial p} &= \frac{1}{r_g T} \end{aligned} $$

* For a mixture of both:

$$ \begin{aligned} \rho &= \frac{1}{\frac{1-\lambda}{\rho_l} + \frac{\lambda r_g T}{p}} \\ \frac{\partial \rho}{\partial p} &= \frac{\lambda \rho_l^2 r_g T}{(p(1-\lambda) + \lambda r_g T \rho_l)^2} \end{aligned} $$

## 2.2. Tested Cavitation Algorithms

As mentioned earlier, in the present work, a varying density model is used for the fluid—a mixture of water and air bubbles—ensuring mass conserving throughout the contact. Even if this pseudocavitation model proves to be accurate [13], it is nevertheless compared to the penalty method and to the Elrod algorithm: the former, for its ease of use, and the latter, for its recognized accuracy.
### 2.2.1. Elrod's Algorithm, the Reference

The full-film area and the cavitation area ($p = P_c$) are explicitly separated: in the former, the Reynolds equation applies, whereas in the latter, a pure Couette flow is modeled with the lubricant fraction $θ$: $∇(6ρθuh) = 0$. Both areas are handled in a single modified Reynolds equation thanks to a switch function that de/activates the Poiseuille contribution.

For each node $i$, if $p_i > P_c$, then $θ_i = 1$; otherwise, if $p_i = P_c$, then $θ_i < 1$. It is undoubtedly the most used algorithm because it guarantees the JFO conditions, ensuring mass conserving everywhere inside $Ω$, and no parameter needs to be tuned. Elrod's algorithm is therefore used here as the reference to compare the penalty method and the general model.

### 2.2.2. The Penalty Method

A means to force sub-cavitation pressure to remain at the cavitation pressure $P_c$ consists of adding a penalty term in the equation to solve. The penalty correction must act as a source term whenever the pressures fall below $P_c$ and it must vanish whenever the pressures are above $P_c$. Thus, modifying the Equation (1), one obtains:

$$ \nabla \left( \frac{\rho h^3}{\mu} \nabla p \right) = \nabla (6\mu\rho h) + \kappa (p - P_c)^{-} \quad (5) $$

where $(p - P_c)^-$ is the $(p - P_c)$ negative part and $κ$ an arbitrary high coefficient that makes negligible the other terms. If $p > P_c$, then $(p - P_c)^- = 0$; otherwise, $(p - P_c)^- \simeq 0$.

For each finite element $e$, a local linear system is written in the form $[K_e]\{\delta p\} = \{B_e\}$ and $κ$ is set to $τ \times \max\{B_e\}$ and $τ = 10^k$. For a fair comparison of the penalty method to the Elrod algorithm, the parameter $τ$ must be well chosen.

To do so, the values $k \in [2, 8]$ are tested on a rough parabolic slider, Figure 2, and the results are compared to Elrod's results. The parabolic slider highlights the pressure build-up ending with a large cavitated zone; the multi-lobed slider highlights the film reformations; and the flat slider is the challenging configuration for which only the deterministic model gives satisfactory results, so far.

<figure><img src="image_6.png" /><figcaption>Figure 2. Lubricated contact domain $Ω$, $(x, y) \in [0, 10] \times [-5, 5]$, and slider profiles at $y = 0$.</figcaption></figure>
Figure 3 shows the differences between the penalty method pressures and JFO pressures for the parabolic slider (blue curve, scale on the right). Several orders of magnitude of τ parameter are represented, yet no difference between the curves is observed in the full-film zone. However, as expected, tiny differences are located in the cavitated area, zoomed in on Figure 4.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Influence of $τ = 10^k$ on the pressure for the parabolic slider at $y = 0$. The film thickness is the blue curve whose scale is shown on the right.</figcaption></figure>

If closer attention is paid to the cavitated area, Figure 4, the penalty results become indistinguishable for $k > 4$.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Influence of $τ = 10^k$ on the pressure for the parabolic slider at $y = 0$. Zoom around $x = 8$ mm.</figcaption></figure>

As suggested by Table 3, the computing time increases with $k$, and values above 5 bring no additional benefit. Therefore, $k = 5$ is chosen to test the method further.

It is to be noted that no additional variable is required, unlike Elrod's algorithm, and the implementation is quite easy: it surely explains why many authors choose the penalty method for Hertzian lubricated contact.
<table><caption>**Table 3.** Influence of $\tau$ on the results. $n_{it}$ stands for the number of Newton-Raphson iterations until convergence, and cpu is the computing time.</caption><thead><tr><th></th><th>JFO</th><th>k = 2</th><th>k = 3</th><th>k = 4</th><th>k = 5</th><th>k = 6</th><th>k = 7</th><th>k = 8</th></tr></thead><tbody><tr><td>n<sub>it</sub></td><td>33</td><td>29</td><td>34</td><td>39</td><td>43</td><td>47</td><td>52</td><td>63</td></tr><tr><td>cpu (s)</td><td>317.7</td><td>266.7</td><td>309.5</td><td>361.9</td><td>402.2</td><td>434.7</td><td>606.3</td><td>614.1</td></tr><tr><td>Load (N)</td><td>480.7</td><td>481.7</td><td>486.2</td><td>487.3</td><td>487.5</td><td>487.5</td><td>487.6</td><td>487.6</td></tr></tbody></table>

When it comes to film reformation, to what extent is the penalty method not advantageous because of its non-mass-conserving character? Figure 5 answers the question by highlighting the large differences between both methods on a multi-lobed rough surface. The penalty method does not suit our needs, namely a rough flat slider with scattered film reformations.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Comparison of pressure profiles obtained in JFO conditions (JFO) and with the penalty method (PEN), with a multi-lobed slider at $y = 0$.</figcaption></figure>

### 2.2.3. The Lubricant General Model

The LGM model—also referred to as MIX, as in 'MIXture', in the present paper—lays on the assumption that the lubricant contains a constant mass fraction $\lambda$ of gas. It is a rather simple model inspired by physical considerations. When the lubricant is stretched, the gas expansion makes the lubricant volume increase. Following Braun and Hannon classification [4], the LGM belongs to the pseudocavitation models.

One way of choosing a suitable value for the parameter $\lambda$ is to fit numerical results on experimental data, as proposed by Brunetière [13] with $m = 7.5$ in the relationship $\lambda = m \times 10^{-5}$, $10^{-5}$ being the common order of magnitude.

Another way of determining $m$ is to consider Henry's law that states that the amount of dissolved gas in a liquid is proportional to its partial pressure above the liquid. The proportionality factor is called Henry's law constant [39].

Then, $C_a = H_a^{cp} p$, where $C_a$ (mol m$^{-3}$) is the concentration of a species in the liquid phase, $p$ (Pa) is the partial pressure of that species in the gas phase under equilibrium conditions, and $H_a^{cp}$ is Henry's constant.

Considering the amount of dissolved oxygen and nitrogen in water, we can write:

$$ \begin{cases} C_{O_2} = H_{O_2}^{cp} P_{atm} & \text{with} \quad H_{O_2}^{cp} = 1.3 \times 10^{-5} \text{ mol m}^{-3} \text{ Pa}^{-1} \\ C_{N_2} = H_{N_2}^{cp} P_{atm} & \text{with} \quad H_{N_2}^{cp} = 6.5 \times 10^{-6} \text{ mol m}^{-3} \text{ Pa}^{-1} \end{cases} $$
Because partial pressures add, the mass in kilogram of dissolved air per liter of water is $\lambda'$:

$$ \lambda' = (0.21 \times C_{O_2} \times 32 + 0.79 \times C_{N_2} \times 28) \times 10^{-6} \\ = 2.3 \times 10^{-5} \text{ kg L}^{-1} $$

where for oxygen—molar mass 32 g mol⁻¹, atmospheric pressure fraction 21%—the mass dissolved per liter of water is 0.21 × C₀₂ × 32 × 10⁻⁶. Finally, the mass fraction $\lambda$ is $\lambda'\rho$, which yields the same value.

$H_a^{cp}$ values are taken from the huge database provided by Sander [40]. The result is of the same order as the previous one (7.5 × 10⁻⁵). It can however be objected that dissolved gas implies gaseous cavitation, which is not the assumption for the LGM. However, according to Grando et al. [41], absorption usually occurs at a much slower rate than release, and the liquid may not be able to absorb the gas in the flow time available during the positive pressure region. Even if the latter statement relates to a journal bearing, we suppose it to be applicable here. The scenario proposed is hence that once the dissolved gas is released in a cavitation area, the lubricant remains a two-phase liquid.

To figure out the importance of $\lambda$, a parametric study is carried out to quantify its effect on the numerical results.

Figures 6 and 7 show the pressure profile along $y = 0$ for different values of $m$, with $\lambda = m \times 10^{-5}$ and JFO’s model. Because the global behavior is the same, the detailed view highlights the subatmospheric pressures in a cavitated zone. High values of $m$ mean more dissolved gas in the lubricant, so the volume increases sooner and the pressure decreases at a slower rate. As a consequence, with increasing values of $\lambda$, the pressure has higher values below the atmospheric pressure in the cavitated zones. It is worth noting that the cavitation pressure chosen for the JFO model is the atmospheric pressure, which is higher than observed values. However, compared to the full-film pressures, the impact remains limited. The same applies to the resulting load, for which the cavitated areas account for a little.

<figure><img src="image_10.png" /><figcaption>Figure 6. General model (MIX) pressure with different values of $m$ in $\lambda = m \times 10^{-5}$ for the multi-lobed slider at $y = 0$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Zoom on a cavitated zone of the multi-lobed slider at $y = 0$—the general model (MIX) pressure is plotted with different values of $m$ from $\lambda = m \times 10^{-5}$ and compared to Elrod's algorithm.</figcaption></figure>

$\lambda$ is chosen so that it provides a small computing time while keeping the normal load $N$ close to JFO's; as Table 4 suggests, $\lambda = 5.0 \times 10^{-5}$, which is rather close to $7.5 \times 10^{-5}$.

<table><caption>**Table 4.** Influence of $\lambda$ on the results. $n_{it}$ stands for the number of Newton-Raphson iterations until convergence, and $cpu$ is the computing time.</caption><thead><tr><th></th><th>JFO</th><th>m = 1</th><th>m = 3</th><th>m = 5</th><th>m = 7</th><th>m = 9</th></tr></thead><tbody><tr><td>n<sub>it</sub></td><td>22</td><td>19</td><td>17</td><td>16</td><td>18</td><td>22</td></tr><tr><td>cpu (s)</td><td>231.9</td><td>181.4</td><td>161.6</td><td>152.2</td><td>171.1</td><td>226.8</td></tr><tr><td>Load (N)</td><td>127.9</td><td>125.7</td><td>127.1</td><td>127.8</td><td>128.0</td><td>127.8</td></tr></tbody></table>

To conclude this chapter, the LGM pseudocavitation model, once tuned, is chosen here, as it offers accuracy and efficiency close to those of a classical JFO model. Moreover, the cavitation boundaries are mesh-independent, which gives it a definite advantage.

## 2.3. The Multiscale Approach

The multiscale approach introduced by Brunetière and Francisco [38] is detailed hereafter. Let us consider a domain $\Omega$ as illustrated in Figure 8. $\Omega$, the hatched zone, is discretized according to $n_e^t \times n_e^t$ elements mesh. The superscript $^t$ refers here to “Top Scale” (TS). Each TS element is called a macro-element because it is also discretized, thus providing a “Bottom Scale” (BS) $n_e^b \times n_e^b$ elements mesh. At the BS level, the whole $\Omega$ mesh is $(n_e^t \cdot n_e^b) \times (n_e^t \cdot n_e^b)$ elements: it is the fine mesh on which deterministic simulations are carried out. In the present case, it is a square regular mesh, but the method can apply to any configurations.

To help the understanding, a variable value is located on the grid with two indices (row, column). Once the pressures $P_{ij}$ are initialized at the TS level, a TS element, “D”, for instance, in Figure 8a, is considered. The parameters are the four nodal pressures $\{P_{32}, P_{42}, P_{43}, P_{33}\}$, interpolated on the macro-element boundary $\Gamma_D$. The variables $p_{ij}$, the inner pressures, are the unknowns determined with the Reynolds equation resolution, Figure 8d.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Overview of the multiscale approach. (a) The top-scale (TS) mesh is a 4 × 4 linear quadrangle mesh. Each TS element (macro-element) is subdivided to form a bottom-scale (BS) mesh, e.g., macro-elements A–D. The TS pressure is interpolated on BS elements (b,c), and then, the Reynolds equation is solved to determine the BS pressure (d). (e) The mass flow is calculated on the boundary $\Gamma_1 \cup \Gamma_2 \cup \Gamma_3 \cup \Gamma_4$, and it is distributed on the corner nodes. Once the macro-elements A, B, C, D are treated, the mass flow balance can be performed on the red node (f).</figcaption></figure>

Each macro-element can be treated independently, which allows for multithreaded computations. When the whole set of macro-element computations is ended, the macro-element boundary mass flows are calculated, Figure 8e, and nodal mass flows are distributed on the macro-element node corners. Thereafter, for each TS node $IJ$, the residual nodal mass flow $Q_{IJ}$ is calculated. At the TS level, the pressures are considered correct if $Q_{IJ} \approx 0$.

One of the key ideas of the multiscale iterative scheme is to perturb the TS nodal pressures to assess the effects on the TS nodal mass flows, Figure 9a,b. Hence, properly modifying the nodal pressures makes the TS nodal mass flow increase: it is a classical Newton-Raphson procedure where the derivatives of the mass flow balance with respect to the pressure must be determined. Here, numerical derivatives are used with $\delta P = P/100$.

When the system $[\partial Q/\partial P]\{\delta P\} = \{Q\}$ is determined and solved, the TS pressures are updated, and the process starts a new iteration until a stop criterion is met. Because all macro-elements are treated in the same way, the multiscale approach is called full multiscale (FMS).
<figure><img src="image_13.png" /><figcaption>**Figure 9.** How the pressure affects the mass flow at the macro-element corner. (a) The mass flow is determined on the macro-element boundary and distributed on the corners. (b) A small variation is introduced in the corner pressure and the mass flows are updated. When all macro-elements are treated, the TS matrix [$\partial Q/\partial P$] can be assembled.</figcaption></figure>

For better multiscale time computing performances and accuracy, Brunetière and Francisco [38] recommend two macro-elements per unit *Rcl* as a rule of thumb. As 1024 is a power of 2, the nearest values are $Ne_x = Ne_y = 128$ and $Ne_x = Ne_y = 256$. As the results do not differ much for both configurations, the chosen grid is $Ne_x = Ne_y = 128$.

The results of the deterministic method—classical FEM resolution of the Reynolds equation—are used as the reference. Two accuracy criteria are chosen: the normal load *N* as a global criterion and the maximum of the film pressure $P_{max}$ as a local criterion. Two computing times are recorded: the total computing time—sum over all threads—$T_t$ and the apparent computing time $T_a$. The computing performance is assessed dividing the computing times by the computing time of the deterministic simulation. The same goes for the computing accuracy which is obtained dividing the normal *N* load and the maximum of the film pressure $P_{max}$ by the deterministic related values.

# 3. Results and Discussion

## 3.1. Deterministic Case

The solver MUMPS v5.4.1—Multifrontal Massively Parallel Solver [42,43]—is used to iteratively solve the linear systems {$R$} + [$\partial R/\partial P$]{$\delta P$} = 0. As expected, the computing time required for the deterministic model decreases when the minimum of the film thickness $h_{min}$ increases; see Figure 10. It is worth noting that the computing time decreases almost linearly (in a log-log scale) until there is no more carrying capacity, *N* = 10.0 N—the atmospheric pressure loading. From a numerical point of view, it also appears that the ratio of the total computing time to the apparent computing time is a constant around 3. It is as if three threads were simultaneously used during the whole case resolution. Actually, there are eight available threads, but the numerical treatments do not fully use them, as shown in Figure 11.
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Minimum of the film thickness $h_{min}$ influence on the computing time of the deterministic model. The left scale is used for the load carrying capacity *N*, and the right scale is used for the computing time (total $T_t$ and apparent $T_a$).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11.** CPUs usage of MUMPS when solving a deterministic case. The lines represent the different threads.</figcaption></figure>

## 3.2. FMS, Full Multiscale

In the FMS configuration, the whole set of the macro-elements is subdivided. It is reminded that there are 128 × 128 macro-elements; each macro-element is subdivided into 8 × 8 elements, for an overall number of 1025 × 1025 nodes.

Figure 12 clearly shows that, first, the errors on *N* and $P_{max}$ do not exceed 6%. The computing time gain is roughly between 20 and 70%, which is less promising than the earlier results, obtained in incompressible conditions [38]. This matter of fact is not surprising: modeling a gas-liquid mixture leads to many more iterations at both levels. The interesting point is that for severe conditions—$h_{min}$ = 0.05 µm, $Sq$ = 0.1 µm—the computing time is a third of the deterministic one, with nonetheless a rather good accuracy.

Concerning the curves, the evolution is not monotonous. Even if the convergence criterion is set to an order of magnitude lower, the results remain unchanged. Therefore, the variations are not explained by the convergence criteria on both levels. In addition, there are no more variations in the maximum of the film pressure $P_{max}$ than on the normal load *N*; it suggests that the variations are not due to the local character of a result. By now, we do not have a satisfactory explanation for the load and the pressure shapes. As for the computing time, things are different. Depending on the pressure set at the corners of a macro-element, more or less iterations are needed: if the corner pressure values differ largely from a top-scale iteration to the next, a subsequent number of iterations will occur at the bottom scale.
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Minimum of the film thickness $h_{min}$ influence on the performance of the FMS model. The results are the differences between the FMS and the deterministic models. The left scale is used for the results of relative differences between the models, and the right scale is used for the apparent computing time $T_a$ relative differences.</figcaption></figure>

Compared to the deterministic solution, the FMS approach qualitatively exhibits the same pressure field. As shown in Figure 13, and despite a 6% $P_{max}$ difference, the pressure fields are very close to each other.

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Pressure fields: deterministic model (DET) and full multiscale (FMS). $h_{min} = 0.05$ µm.</figcaption></figure>

As for the residual maps, Figure 14, it is interesting to remark that the trends are qualitatively the same: the higher residuals zones roughly match. Note the grid on the FMS picture: the residuals are set to null on the lines, except the intersection nodes. Indeed, the intersection nodes are the top-scale nodes, but elsewhere on the grid, the nodes are boundary conditions for the inner macro-element nodes. It is important to keep in mind that the convergence criterion, whilst being of the same kind—as often, based on $\frac{\delta P}{P}$ values—cannot be thoroughly compared: in the deterministic model, the whole pressure field is updated at once, whereas in the FMS model, it is updated by parts and, more specifically, for each macro-element. That explains why the residuals are qualitatively compared.
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Residual field: deterministic model (DET) and full multiscale (FMS). $h_{min}$ = 0.05 µm.</figcaption></figure>

### 3.3. *Beyond FMS: Hybrid Multiscale (HMS)*

When dealing with the HMS—Hybrid Multiscale—the macro-elements are filtered according to the related film thickness. The basic idea is to only refine the macro-elements for which the roughness may influence the lubrication. Simply put, a film thickness threshold $h_T$ is set: below $h_T$, the macro-elements are refined as in the FMS, whereas each macro-element above $h_T$ is treated as a single finite element. The film thickness threshold is deduced from the percentage of the bottom-scale elements set by the user.

Figure 15 is a part of the parametric study that aims to guide the user regarding the percentage of the BS elements to set as a function of the minimum of the film thickness $h_{min}$. The whole study can be found in SM, Section S3.

Unfortunately, there is no clear strategy rising from the study. To help the understanding of the plots in Figure 15, in the light green areas, the error ($N$ and $P_{max}$) is below 10%, whereas in the dark green areas, the error is below 5%. In both cases, the computing time is at least half the deterministic computing time. When the four configurations are compared—$h_{min}$ = 0.05, 0.15, 1.00, and 5.00 µm—it is hard to give any advice on the choice of the percentage of the BS elements to set. Surprisingly, in the harsher case, no BS element is needed to reach very satisfactory results, whereas 100% is necessary for the triple of $h_{min}$. So far, different tracks were followed in vain to solve this paradox, and the next explores deals with the BS criterion based on the average film thickness. The parametric study contains other values of $h_{min}$; see SM Section 3 for the complete study.

Anyway, if approximated results are sought—up to a 10% error on $P_{max}$ or $N$—then the HMS is the right tool because of its fastness (a tenth of the deterministic computing time). As an example, the case $h_{min}$ = 0.05 µm is handled about ten times faster with half BS elements. The counterpart is obviously the precision; however, it remains below 10%. Figure 16 shows that the pressure field is in good agreement with the deterministic case. Moreover, the convergence residuals are uniformly spread, like the BS elements in Figure 17, which suggests that the pressure field is a fairly good approximation everywhere.
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Influence of the percentage of bottom scale (BS) on N, Pmax, and Tr. (a) hmin = 0.05 µm, (b) hmin = 0.15 µm, (c) hmin = 1.00 µm, (d) hmin = 5.00 µm. The values are normalized with respect to the deterministic case. The green areas correspond to errors below 5% and a computing time at least half the deterministic one. The light green areas correspond to errors below 10% and a computing time at least half the deterministic one.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 16.** Pressure (a) and residual (b) fields: HMS model, 50% BS-elements, hmin = 0.05 µm.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 17.** BS-elements distribution, 50% BS elements, $h_{min}$ = 0.05 µm. The black squares are the BS elements.</figcaption></figure>

The FMS, and even more the HMS, provides for the user with fast approximations of the pressure field. However, a natural question arises: what about simply reducing the mesh size in the deterministic case?

### 3.4. Mesh Size Reduction

Four meshes are used to assess the mesh influence on $P_{max}$, $h_{min}$, and the relative apparent computing time differences: the 513 × 513 nodes, 641 × 641 nodes, 769 × 769 nodes, and 897 × 897 nodes. The four models are deterministic. The recurring question about coarsening is, what kind of interpolation would be best? Unfortunately, there is no good answer because not all the information from a fine grid can be transferred to coarser grids. Specifically, the roughness height extrema are likely to disappear during the coarsening step. Therefore, between the three interpolation methods that were tested—namely the approximation, linear, and cubic spline—the linear method is chosen because it is simply the most common. The reader is referred to SM, Figure S8, for being convinced that the linear interpolation remains a good compromise.

As shown in Figure 18, when using the coarsest mesh 513 × 513, the computing time falls to values of the same order as those for the HMS. The spikes located at $h_{min}$ = 0.15 µm; 3 µm are explained by slightly higher computing times that are emphasized by the relative differences calculus. The error on $P_{max}$ and $h_{min}$ is even better than for the HMS, Figures 19 and 20.

<figure><img src="image_22.png" /><figcaption>**Figure 18.** Computing time reduction as a function of $h_{min}$ for different meshes.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 19.** Load error as a function of $h_{min}$ for different meshes.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 20.** Pressure error as a function of $h_{min}$ for different meshes.</figcaption></figure>

If Figure 13—DET—and Figure 21 are compared, it appears that indeed the results are qualitatively very close. The explanation is as follows. With six elements per autocorrelation length—a roughly 513 × 513 nodes mesh, because $Rcl \approx 80$—the mean shape of the surface is well caught. Yet, the main contributors to the lubrication are long wavelengths, which explains that finer meshes bring few additional details.
<figure><img src="image_25.png" /><figcaption>**Figure 21.** Pressure field: deterministic model, 513 × 513, $h_{min}$ = 0.05 µm.</figcaption></figure>

**4. Conclusions**

In the incompressible case, the FMS model provides for a very fast tool to solve the Reynolds equation. However, strong departures from linearity affect the FMS/HMS model efficiency for compressible fluids. Indeed, with the cavitation phenomenon, the computing process iterates much more: iterations are not only needed at the top-scale level but also at the bottom-scale level. In addition, some relaxation coefficients must be tuned on both levels, and particularly when the ratio $\sigma/h_{min}$ exceeds 1.

Another a priori negative point is that reducing the mesh density allows for the same computing time reduction along with the same result accuracy. And yet, the slider is flat, meaning that the pressure build-up is only generated by the roughness ($Sq$ = 0.1 µm). Despite the preceding points, it is worth noting that the computing times—when the HMS is used with a few BS elements—remain attractive, compared to the deterministic model, whether it is reduced or not.

The benefits of the FMS/HMS approach are much more expected in the following cases:

1. Shorter wavelengths—$Rcl$ > 100 and $L_{x,y}$ < 10mm;

2. Rough texturized surfaces—square dimples being modeled with simple FEM macro-elements, and rough contacting parts discretized with BS elements;

3. As many threads as there are macro-elements—which leads to a GPU implementation of the numerical code.

As for the numerical aspects, additional computing time can be saved.

1. The TS element boundaries are updated once the whole TS element batch is processed. However, some TS elements converge slowly—mainly because of a local narrow slider gap. Therefore, a new criterion has to be set up to locally guarantee the best compromise accuracy/iteration number.

2. The TS element boundary update must be monitored because important changes in the TS pressure field affect the four other connected macro-elements—in particular, oscillations are undesirable.

3. The heights of the element boundaries are the result of the domain division, and this can lead to rough relief for some of them, with convergence problems. In a future work, a sensitivity analysis will study the effect of numerically smoothed boundaries on the slider lubrication. We are confident that the results will not change much and that the Reynolds equation will be solved faster on TS elements.
To conclude the present work, the FMS/HMS approach is accurate and fast due to a highly parallelizable structure, and promising keys to improvement make this technique a good candidate for the lubrication of rough parallel surfaces.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/lubricants10120329/s1, Section S1: The SU-PG method; Subsection S2.1: The weak formulation using the Bubnov-Galerkin method; Subsection S2.2: Couette term upwind weighting; Subsection S2.3: The Newton-Raphson scheme; Section S3: HMS parametric study for the flat rough slider; Figure S8: Influence of the interpolation type on the coarse mesh results related to the fine mesh.

**Author Contributions:** Conceptualization, N.B. and A.F.; methodology, N.B. and A.F.; software, N.B. and A.F.; validation, N.B. and A.F.; writing—original draft preparation, A.F.; writing—review and editing, N.B. and A.F.. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work pertains to the French government program Investissements d'Avenir (LABEX INTERACTIFS, reference ANR-11-LABX-0017-01, and EUR INTREE, reference ANR-18-EURE-0010).

**Data Availability Statement:** Not applicable

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

**Abbreviations**

The following abbreviations are used in this manuscript:

<table><tr><td>BS</td><td>Bottom Scale</td></tr><tr><td>FEM</td><td>Finite Element Method</td></tr><tr><td>FMS</td><td>Full Multiscale</td></tr><tr><td>HMS</td><td>Hybrid Multiscale</td></tr><tr><td>LGM</td><td>Lubricant General Model</td></tr><tr><td>MIX</td><td>MIXture, same as GLM</td></tr><tr><td>TS</td><td>Top Scale; TS element = a macro-element.</td></tr></table>

**Nomenclature**

Reynolds equation

e $e^{th}$ finite element

h film thickness

p pressure

l_e streamline length of an element e

N slider normal load

$N_i$ shape and weighting function at node i

$P_c$ cavitation pressure

$R_i$ Reynolds equation residual at node i

u slider speed vector

$[B_e]$ $e^{th}$ element right-hand side

$[K_e]$ $e^{th}$ element elementary matrix

$\delta p_j$ pressure increment at node j

$\mu$ fluid dynamic viscosity

$\Omega$ lubricated contact domain

$\rho$ fluid density

Penalty method

$\kappa$ arbitrary chosen coefficient for the penalty method

$\tau$ penalty coefficient, $\tau = 10^k$

Elrod algorithm

$\theta$ fluid volume fraction
General model

$p$ partial pressure

$r_g$ the specific gas constant; for dry air $r_g = 287.0 \text{ J kg}^{-1} \text{ K}^{-1}$

$C_a$ liquid concentration of a species a

$H_a^{cp}$ Henry's constant of a species a

$P_{atm}$ atmospheric pressure, $P_{atm} = 1.01 \times 10^5 \text{ Pa}$

$\lambda$ mass fraction

Multiscale

$h_T$ film thickness threshold

$n = ne_x = ne_y$, total number of BS nodes

$n_e^b$ number of bottom-scale (BS) elements of a TS element, along x or y

$n_e^t$ number of top-scale (TS) elements, or macro-elements, along x or y

BS Bottom scale

P TS pressure

Q TS nodal mass flow

TS Top scale

*{b}* subdivided macro-element FEM right-hand side

*[k]* subdivided macro-element FEM system matrix

References

1. Wu, S.R. A Penalty Formulation and Numerical Approximation of the Reynolds-Hertz Problem of Elastohydrodynamic Lubrication. *Int. J. Eng. Sci.* **1986**, *24*, 1001–1013. [CrossRef]

2. Elrod, H.G.; Adams, M.L. A Computer Program for Cavitation and Starvation Problems. In Proceedings of the 1st Leeds—Lyon Symposium on Tribology, Leeds, UK, September 1974; Mechanical Engineering Publications Ltd.: Leeds, UK, 1974.

3. Elrod, H.G. A Cavitation Algorithm. *J. Lubr. Technol.* **1981**, *103*, 350–354. [CrossRef]

4. Braun, M.J.; Hannon, W.M. Cavitation Formation and Modelling for Fluid Film Bearings: A Review. *Proc. Inst. Mech. Eng. Part J. Eng. Tribol.* **2010**, *224*, 839–863. [CrossRef]

5. Bai, L.; Yan, J.; Zeng, Z.; Ma, Y. Cavitation in Thin Liquid Layer: A Review. *Ultrason. Sonochemistry* **2020**, *66*, 105092. [CrossRef]

6. Belova, V.; Krasowska, M.; Wang, D.; Ralston, J.; Shchukin, D.G.; Möhwald, H. Influence of Adsorbed Gas at Liquid/Solid Interfaces on Heterogeneous Cavitation. *Chem. Sci.* **2012**, *4*, 248–256. [CrossRef]

7. Li, X.S.; Song, Y.; Hao, Z.R.; Gu, C.W. Cavitation Mechanism of Oil-Film Bearing and Development of a New Gaseous Cavitation Model Based on Air Solubility. *J. Tribol.* **2012**, *134*, 031701. [CrossRef]

8. Song, Y.; Gu, C.W. Development and Validation of a Three-Dimensional Computational Fluid Dynamics Analysis for Journal Bearings Considering Cavitation and Conjugate Heat Transfer. *J. Eng. Gas Turbines Power* **2015**, *137*, 122502. [CrossRef]

9. Hao, Z.R.; Gu, C.W. Numerical Modeling for Gaseous Cavitation of Oil Film and Non-Equilibrium Dissolution Effects in Thrust Bearings. *Tribol. Int.* **2014**, *78*, 14–26. [CrossRef]

10. Ding, A.; Ren, X.; Li, X.; Gu, C. A New Gaseous Cavitation Model in a Tilting-Pad Journal Bearing. *Sci. Prog.* **2021**, *104*, 00368504211029431. [CrossRef]

11. Ransegnola, T.; Sadeghi, F.; Vacca, A. An Efficient Cavitation Model for Compressible Fluid Film Bearings. *Tribol. Trans.* **2021**, *64*, 434–453. [CrossRef]

12. Magaletti, F.; Gallo, M.; Casciola, C.M. Water Cavitation from Ambient to High Temperatures. *Sci. Rep.* **2021**, *11*, 20801. [CrossRef]

13. Brunetière, N. A General Model for Liquid and Gas Lubrication, Including Cavitation. *J. Tribol.* **2018**, *140*, 021702. [CrossRef]

14. Grützmacher, P.G.; Profito, F.J.; Rosenkranz, A. Multi-Scale Surface Texturing in Tribology—Current Knowledge and Future Perspectives. *Lubricants* **2019**, *7*, 95. [CrossRef]

15. Patir, N.; Cheng, H.S. An Average Flow Model for Determining Effects of Three-Dimensional Roughness on Partial Hydrodynamic Lubrication. *J. Lubr. Technol.* **1978**, *100*, 12–17. [CrossRef]

16. Patir, N.; Cheng, H.S. Application of Average Flow Model to Lubrication Between Rough Sliding Surfaces. *J. Lubr. Technol.* **1979**, *101*, 220–229. [CrossRef]

17. Harp, S.R.; Salant, R.F. An Average Flow Model of Rough Surface Lubrication with Inter-Asperity Cavitation. *J. Tribol.* **2000**, *123*, 134–143. [CrossRef]

18. Elrod, H.G. A General Theory for Laminar Lubrication with Reynolds Roughness. *J. Lubr. Technol.* **1979**, *101*, 8–14. [CrossRef]

19. Bayada, G.; Faure, J.B. A Double Scale Analysis Approach of the Reynolds Roughness Comments and Application to the Journal Bearing. *J. Tribol.* **1989**, *111*, 323–330. [CrossRef]

20. Bayada, G.; Martin, S.; Vázquez, C. An Average Flow Model of the Reynolds Roughness Including a Mass-Flow Preserving Cavitation Model. *J. Tribol.* **2005**, *127*, 793–802. [CrossRef]
21. de Boer, G.N.; Gao, L.; Hewson, R.W.; Thompson, H.M.; Raske, N.; Toropov, V.V. A Multiscale Method for Optimising Surface Topography in Elastohydrodynamic Lubrication (EHL) Using Metamodels. *Struct. Multidiscip. Optim.* **2016**, *54*, 483–497.
[CrossRef]

22. de Boer, G.N.; Gao, L.; Hewson, R.W.; Thompson, H.M. Heterogeneous Multiscale Methods for Modelling Surface Topography in Elastohydrodynamic Lubrication Line Contacts. *Tribol. Int.* **2017**, *113*, 262–278. [CrossRef]

23. de Boer, G.; Almqvist, A. On the Two-Scale Modelling of Elastohydrodynamic Lubrication in Tilted-Pad Bearings. *Lubricants* **2018**, *6*, 78. [CrossRef]

24. Han, Y.; Meng, Q.; de Boer, G. Two-Scale Homogenization of Hydrodynamic Lubrication in a Mechanical Seal with Isotropic Roughness Based on the Elrod Cavitation Algorithm. *Proc. Inst. Mech. Eng. Part J J. Eng. Tribol.* **2022**, *236*, 359–385. [CrossRef]

25. Rom, M.; König, F.; Müller, S.; Jacobs, G. Why Homogenization Should Be the Averaging Method of Choice in Hydrodynamic Lubrication. *Appl. Eng. Sci.* **2021**, *7*, 100055. [CrossRef]

26. Sen, O.T.; Akalin, O. A Multiscale Cylinder Bore Honing Pattern Lubrication Model for Improved Engine Friction. *SAE Int. J. Fuels Lubr.* **2019**, *12*, 155–166. [CrossRef]

27. König, F.; Rosenkranz, A.; Grützmacher, P.G.; Mücklich, F.; Jacobs, G. Effect of Single- and Multi-Scale Surface Patterns on the Frictional Performance of Journal Bearings—A Numerical Study. *Tribol. Int.* **2020**, *143*, 106041. [CrossRef]

28. Shisode, M.; Hazrati, J.; Mishra, T.; de Rooij, M.; van den Boogaard, T. Mixed Lubrication Friction Model Including Surface Texture Effects for Sheet Metal Forming. *J. Mater. Process. Technol.* **2021**, *291*, 117035. [CrossRef]

29. Annadi, R.R.; Ismail, S. Tribological Performance of Textured Parallel Sliding Contact under Mixed Lubrication Condition by Considering Mass Conservation Condition and Couple-Stress Parameter. *Proc. Inst. Mech. Eng. Part J J. Eng. Tribol.* **2021**, *235*, 410–422. [CrossRef]

30. Annadi, R.R.; Syed, I. Tribological Study of Multi-Scale Textured Parallel Sliding Contacts by Considering a Mixed Lubrication Regime and Mass Conservation Condition. *Proc. Inst. Mech. Eng. Part J J. Eng. Tribol.* **2022**, *236*, 807–822. [CrossRef]

31. Lebeck, A.O. Parallel Sliding Load Support in the Mixed Friction Regime. Part 1—The Experimental Data. *J. Tribol.* **1987**, *109*, 189–195. [CrossRef]

32. Lebeck, A.O. Parallel Sliding Load Support in the Mixed Friction Regime. Part 2—Evaluation of the Mechanisms. *J. Tribol.* **1987**, *109*, 196–205. [CrossRef]

33. Nyemeck, A.P.; Brunetière, N.; Tournerie, B. A Multiscale Approach to the Mixed Lubrication Regime: Application to Mechanical Seals. *Tribol. Lett.* **2012**, *47*, 417–429. [CrossRef]

34. Brunetière, N.; Francisco, A. Lubrication Mechanisms Between Parallel Rough Surfaces. *Tribol. Lett.* **2019**, *67*, 116. [CrossRef]

35. Brunetiere, N.; Wang, Q.J. Large-Scale Simulation of Fluid Flows for Lubrication of Rough Surfaces. *J. Tribol.* **2014**, *136*, 011701. [CrossRef]

36. Pei, S.; Ma, S.; Xu, H.; Wang, F.; Zhang, Y. A Multiscale Method of Modeling Surface Texture in Hydrodynamic Regime. *Tribol. Int.* **2011**, *44*, 1810–1818. [CrossRef]

37. Pei, S.; Xu, H.; Shi, F. A Deterministic Multiscale Computation Method for Rough Surface Lubrication. *Tribol. Int.* **2016**, *94*, 502–508. [CrossRef]

38. Brunetière, N.; Francisco, A. Multiscale Modeling Applied to the Hydrodynamic Lubrication of Rough Surfaces for Computation Time Reduction. *Lubricants* **2018**, *6*, 83. [CrossRef]

39. Wikipedia. Henry's Law. 2020. Available online: https://en.wikipedia.org/wiki/Henry's_law (accessed on 26 October 2022).

40. Sander, R. Compilation of Henry's Law Constants (Version 4.0) for Water as Solvent. *Atmos. Chem. Phys.* **2015**, *15*, 4399–4981. [CrossRef]

41. Grando, F.P.; Priest, M.; Prata, A.T. A Two-Phase Flow Approach to Cavitation Modelling in Journal Bearings. *Tribol. Lett.* **2006**, *21*, 233. [CrossRef]

42. Amestoy, P.R.; Duff, I.S.; L'Excellent, J.Y.; Koster, J. A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling. *SIAM J. Matrix Anal. Appl.* **2001**, *23*, 15–41. [CrossRef]

43. Amestoy, P.R.; Buttari, A.; L'Excellent, J.Y.; Mary, T. Performance and Scalability of the Block Low-Rank Multifrontal Factorization on Multicore Architectures. *ACM Trans. Math. Softw.* **2019**, *45*, 1–26. [CrossRef]