Multiscale Modeling Applied to the Hydrodynamic Lubrication of Rough Surfaces for Computation Time Reduction

This paper presents a multiscale finite element method applied to the simulation of a lubricating film flowing between rough surfaces. The objective of this approach is to study flows between large rough surfaces needing very fine meshes while maintaining a reasonable computation time. For this purpose, the domain is split into a number of subdomains (bottom-scale meshes) connected by a coarse mesh (top-scale). The pressure distribution at the top-scale is used as boundary conditions for the bottom-scale problems. This pressure is adjusted to ensure global mass flow balance between the contiguous subdomains. This multiscale method allows for a significant reduction of the number of operations as well as a satisfactory accuracy of the results if the top-scale mesh is properly fitted to the roughness lateral scale. Furthermore the present method is well-suited to parallel computation, leading to much more significant computation time reduction.


Introduction
The lubrication regime between surfaces in relative motion is controlled by the duty parameter G = µVL W , where µ is the lubricant viscosity, V the sliding speed, L the contact length and W the applied load [1].When G is reduced, the friction decreases until asperity contact occurs and mixed lubrication begins.Reducing the viscosity of lubricants is a common solution used to decrease friction losses in engines [2] or other applications [3].The decrease in µ leads to a reduction of G and of the thickness of the fluid film between the surfaces.As a consequence, the mixed lubrication tends to become the standard lubrication regime.It is thus of importance to have tools able to simulate this lubrication regime in a reasonable computation time and with sufficient accuracy.
The first step is the analysis of the flow between rough surfaces.The fluid flow is generally described using the Reynolds equation, which was initially solved assuming smooth surfaces [4].The surface roughness must thus be taken into account in the model.One of the problems is that rough surfaces exhibit a multiscale nature, which makes their analysis and simulation difficult [5].To get rid of this complexity, the first models developed in the 1960s assumed simplified geometric configurations.For example, Tzeng and Saibel [6] assumed one-directional striated surfaces.Using an averaging method, they were able to include correction factors in the Reynolds equation.Their approach was improved by Tonder [7].To consider more complex surfaces with two-dimensional patterns, Patir and Cheng introduced the flow factor theory in the 1970s [8,9].They used very small samples of rough surface on which they performed direct numerical simulations also called deterministic simulations.They obtained correction factors for the fluid flow that they incorporated in the averaged Reynolds equation.Their method had much success and was improved by several authors to derive the flow factors theoretically [10,11] or to include micro-cavitation [12].Another method was developed meanwhile.Indeed, assuming that the roughness has a periodic pattern, it is possible to use the homogenization theory [13].This method was extended to take account of cavitation [14], as well as the topography of real surfaces [15,16].This approach was recently further improved with heterogenous homogenization [17,18].However, all these averaging methods suffer a limitation.They are unable to simulate the load experimentally observed between nominally parallel flat surfaces [1].
With the development of computer capacities, it has been possible to perform deterministic simulations.Deterministic means that the conventional Reynolds equation is used while considering the real surface roughness.A very fine mesh is then used to capture the details of the topography of the rough surfaces.It was first used in point contacts [19] and later applied to line contacts [20].More recently, the deterministic solution was used for conformal parallel surfaces [21].It has been made possible to simulate the roughness-induced load generation between parallel surfaces.However, because of the required computation time, the studied surface was of limited extent, and periodical boundary conditions were used, which is not realistic.
There has been a great research effort to reduce the computation time for such problems and consider the multiscale nature of rough surfaces.One of the first efficient methods was based on the multi-grid techniques [22].The principle is to use several mesh densities to remove low frequency residuals in order to accelerate the convergence.This method has been used by many authors [19,20] because of its efficiency.This method is, however, limited to structured regular meshes.A multiscale approach based on a solution combining two mesh levels was proposed by Neymeck et al. [23].The macro-mesh pressure solution serves as boundary conditions for the micro-mesh.This pressure is adjusted to ensure global mass conservation between micro-subdomains.The method was however limited to an unidimensional macro-mesh and cannot be extended to large rough surfaces.More recently, Brunetiére and Wang [24] proposed a mixed solution.The method consists of combining an averaging method for the high frequencies of the roughness and a deterministic solution for the low frequencies.This approach however relies on a filtering process, which can be difficult to perform.
The present paper proposes to extend the initial multiscale principle proposed by Nyemeck et al. [23] to a two-dimensional problem discretized with the finite element method.First, the method is presented.Then, some results showing its performance in terms of reducing the computation time are presented for the case of hydrodynamic lubrication between rigid rough surfaces.Finally, the accuracy of the method is discussed.

Configuration of the Problem
The configuration of the studied problem is presented in Figure 1.A rigid smooth surface is moving along the x-axis at a constant speed V.A stationary and rigid rough surface of size L x by L y is placed at the average distance h 0 from the bottom surface.A viscous fluid separates the two surfaces.The fluid flow between the surfaces is in steady state and governed by the following Reynolds equation: where h is the local film thickness, µ is the viscosity of the fluid, ρ is the density of the fluid, and p is the pressure, which is the unknown of the problem.The film thickness is defined as h = h 0 + h t , where h t is the roughness of the surface.The usual lubrication assumptions (h L x and L y , continuous media, laminar flow, fluid inertia neglected, constant properties through the film thickness, no slip on the walls) hold to obtain this equation.Some additional hypotheses are used: 1.
The fluid is incompressible; 2.
There is no film rupture or cavitation; 3.
The surfaces are rigid.
These assumptions, even if not realistic, are used in this initial work because the aim of the paper is to focus on the multiscale method.Considering a multiphysics problem would make the discussion on the model performance more difficult.

Finite Element Formulation
The Reynolds equation is multiplied by a weighting function N i and is integrated over the whole domain.It constitutes a residual: An integration by parts is then carried out to reduce the order of the derivatives with respect to the pressure.A weak formulation is then obtained: Flow rate along the boundary Γ The second term is the mass flow rate along the boundary weighted by the function N i .The domain Ω is split into n e isoparametric elements with n nodes.The Bubnov-Galerkin method is used, meaning that the weighting functions N i are chosen equal to the shape function of the elements [25].Using the n shape functions of the n nodes, n equations are obtained.Since the pressure is known along the boundary of the domain, it is not necessary to use equations for the nodes located on the boundary of the domain.Thus, the shape functions vanish along the boundary as well as the boundary integral in Equation (3).
Even if the problem is linear, it is useful to consider the more general nonlinear case for situations such as compressible fluids or cavitation.The Newton-Raphson method leads to the following set of equations: where p j is the pressure at node j and δp j a pressure variation.This system of equations is solved using two different direct solvers.The multi-frontal UMFPack (http://faculty.cse.tamu.edu/davis/suitesparse.html) solver in its single-threaded version is used as well as the parallel solver MUMPS (http://mumps.enseeiht.fr/).The two packages are solvers for systems of linear equations with sparse matrices.A Gaussian elimination is performed by transforming the matrix in the product of a lower triangular matrix L by an upper triangular matrix U using multi-frontal methods.The solution is obtained in three steps.First, the system is analyzed and a symbolic factorization is performed.For a given shape of the matrix-non zero locations-the analysis is done once.During the second step, the matrix is numerically factorized, leading to L and U matrices.Finally the system is solved.These two solvers are called from a unique interface, MSOLV (https://tribo-pprime.github.io/MSOLV/index.html).

Multiscale Approach
As shown in Figure 2, the regular method consists of solving the problem equations on a fine mesh containing n e elements and covering the whole domain.In the multiscale approach, the domain is first discretized with a coarse mesh.It is the macro scale or top-scale.The number of elements of this mesh is denoted by n t e .Each of the n t e elements is then finely meshed.This is the micro scale or bottom-scale.In the present paper, all the top elements have the same size, and all the subdomains have the same number of elements n b e .The following relationship between the number of elements is finally obtained: Macro scale mesh Full domain

Multiscale method
Micro scale mesh × number of macrocells The calculation principle and the solution procedure of the multiscale method are described in Figure 3.An initial guess of the pressure distribution is first assumed at the top-scale (it can be a uniform pressure distribution).This pressure distribution is then used as boundary conditions for the micro scale problems.The pressure solutions at the bottom-scale are obtained using the approach presented in Section 2.2.Note that there are n t e independent subproblems to be solved simultaneously.Knowing the pressure distributions at the bottom-scale, it is then possible to calculate the mass flow rate along the boundary of the subdomains: In this equation, N t i is the top-scale shape function.There is thus one flow rate value per top-scale node.These flow rates are exported to the top-scale where they are summed to obtain the residual flow rate for each top-scale node: By making small variations in the pressure at the boundaries of the bottom-scale problems, derivatives of the flow rate with respect to the pressure can be calculated.Since only the right-hand sides of the system of equations are modified, the solution is obtained very quickly.Finally, the pressure at the top-scale is adjusted to ensure mass flow balance between contiguous elements.This means that the flow rate residuals R t t must vanish: Since the problem is linear, the solution at the top-scale is obtained after one iteration at the bottom-scale (regardless of the initial pressure estimate).If the pressure at the bottom-scale is required, for example to calculate the load, a second simulation is performed at the bottom-scale with the new pressure distribution from the top-scale.Note that when cavitation is included, several iterations are necessary because of the nonlinear behavior of the problem.The scale transition is ensured by assuming a linear pressure distribution along the bottom-scale domains and a global mass conservation between the elements at the top-scale.Local mass conservation is not ensured between contiguous subdomains.These are the main simplifications involved in the present approach.They allow decreasing the computation time but they induce errors in the pressure solution compared to the global deterministic solution.
Finally, it should be noted that the algorithm was coded with modern Fortran language.

Results
The configuration used for the problem is presented in Figure 1 and the parameters are given in Table 1.The rough surface used in the simulation is presented in Figure 4.It is a numerically-simulated 2048 × 2048 bi-Gaussian surface, which is typical of surfaces used in tribological applications [26].This isotropic surface has a smooth plateau with deep valleys, as indicated by the values of SSk and SKu in Table 1.Simulated surfaces are used in this work because it is possible to accurately control their statistical properties.However, a measured surface could also be used.The film thickness is chosen to avoid any asperity contact, something which is not addressed in the present paper.The surface is divided in k = 2-512 subdomains per side for the multiscale simulation.Because of the square shape of the problem, k = n t e .The fluid is assumed to be incompressible and cavitation is neglected in this initial approach.The pressure is equal to zero along the boundary.The resulting pressure is proportional to a reference hydrodynamic pressure: where L = L x = L y is the width of the domain.Note that a square shape was used to simplify the analysis of the results but, more generally speaking, the finite element method used in the model allows for high flexibility in the shape of the domain.The pressure distributions obtained with the deterministic solution and with the multiscale solution with k = 8 are presented in Figure 5.The top-scale pressure gives the main trends of the pressure field obtained with the full solution, while the details are obtained at the bottom level.The bottom-scale solution appears to be very close to the deterministic solution even if there are some differences along the edges of the subdomain where linear pressure boundary conditions are imposed.The objective of the multiscale approach is to reduce the computation time.To asses the efficiency of the method, two types of calculation are performed.First, the mono-processor solver UMFPack was used to solve both the bottom-and top-scale problems.Second, the multiprocessor solver MUMPS was used for the top-scale whereas the mono-processor solver was used in concurrent threads for the bottom-scale problems.Indeed, since there were n t e independent bottom-scale problems to solve, this can be done in parallel using the Open-MP library.For the deterministic full solution (k = 1), only MUMPS was able to solve the problem because its 64 bits integer implementation allows for large memory allocation.The computation time obtained in this case was used as a reference.
Figure 6 presents the evolution of the computation time regarding the number of subdomains per side, for both the cases of sequential and parallel computations.The curves exhibit the same behavior: the computation time decreased when the number of subdomains was increased, then it reached a plateau.Based on the parallel computation curve, it can be shown that the present method could divide the computation time by about a factor of 13 when compared to the deterministic solution.The use of parallel computation could reduce the computing time by an additional factor.In the present case, where eight threads were used, this factor was about four.The multiscale approach was thus very effective at reducing the calculation time.Because of the linear pressure distribution assumed along the edges of the subdomains, some error on the pressure distribution was to be expected.The overall errors were assessed by using the calculated load, which was the integral of the pressure over the domain.Figure 7 presents the value of the load versus the number of subdomains per side, k.The load was normalized by the value obtained with the exact deterministic solution.At low values of k, there was a significant error, which could be higher than 50%.This error tended to decrease when k increased.A small error was obtained for the highest value of k.The results shown in Figure 7 are not surprising because when k increases the top-scale mesh is able to capture almost all the details of the exact solution.An illustration of this assertion is given in Figure 8.When k was small, the pressure distribution was almost uniform on the whole top-scale domain.When k was increased, the details of the pressure distribution appeared more clearly at the top-scale, leading to better boundary conditions for the bottom-scale problems, and hence a better accuracy.

Computation Time
The computation time to solve a set of equations is directly linked to the number of operations to factorize the system.According to Duff [27], the number of operations or flops n f for a sparse system of N unknowns can be approximated by Using an advanced method that includes an analysis step makes the number of flops decrease: [28] In the case of a large domain, such as in the present problem, N is approximately equal to the number of elements of the domain n e , and the number of flops is n f ∝ n 3/2 e .When the multiscale approach is used, there is one domain with n t e elements and n t e domains of n b e = n e n t e elements.Thus the number of flops is This equation is plotted in Figure 6 with n e = 2048 2 and k = n t e and after normalization by n f (k = 1).This equation is able to describe the significant decrease of the computation time observed when k is increased from 1 to 8. Based on this theory, the number of flops and thus the computation time should be reduced by about 2 orders of magnitude when the number of subdomains is set at its optimal value, n t e 230 2 .In the present case, the reduction is by one order of magnitude and the optimal value is n t e = 128 2 .The difference between the theory and the real computation time and more particularly the appearance of the plateau can be explained as follows.When k increases, the size of the bottom-scale matrices decreases.Thus, their sparsity (i.e., the ratio of the number of terms in the matrix to the number of non-zero terms) becomes lower and the efficiency of the solver algorithm is lowered with a number of flops tending to N 2 , which is higher than the expected N 3/2 .

Accuracy of the Method
It has been shown that the accuracy of the method depends on the number of subdomains used in the multiscale approach.The accuracy will be improved if the top-scale mesh is able to capture the main features of the pressure distribution.These features are expected to vary with the lateral size of the surface asperities.To analyze this assumption, several surfaces with different correlation lengths (expressed by Sal according to the ISO standard 25178) were generated and used in the simulations.The effect of Sal on the lateral size of the surface asperities is clearly visible in Figure 9.In addition, a second surface having the same statistical properties (Sal = 0.007L) as the reference surface was also generated.
For each of the six surfaces (the reference surface and five others), the deterministic solution was computed for a reference.Then, multiscale simulations were performed with different k values between two and 512, and the relative error of the computed load was calculated.In Figure 10, this relative error is presented as a function of the ratio of the correlation length Sal to the size L/k of the subdomain.Even if for a given surface, the evolution is not monotonic, a general trend can be observed.The error tends to decrease when the size ratio increases.Indeed, the top-scale mesh is able to capture more accurately the details of the pressure distribution.Moreover, a decrease in the slope is observed when the size ratio is approximately equal to one.This means that when the subdomains have a size less than the correlation length of the surface, the error decreases more rapidly.To reach an error of about 0.01, the size of the subdomain must be less than one-half the correlation length of the surface roughness.Thus, this gives a quantitative indicator for the discretization of the problem.This value can be compared to the optimal number of subdomains given by Equation ( 12) to find the best compromise in terms of accuracy and computation time.In the reference case, Sal = 0.007L, this rule suggests the use of k > 285.This corresponds to the plateau zone where the computation time is minimal (see Figure 6).

Conclusions
A multiscale finite element solution of the Reynolds equation was presented with the aim of studying the fluid flow between rough surfaces demanding very fine meshes.The domain is split into a number of subdomains (bottom-scale mesh) connected by a coarse mesh (top-scale).The pressure distribution at the top-scale is used as boundary conditions for the bottom-scale problems.This pressure is adjusted to ensure the global mass flow balance between contiguous subdomains.This multiscale method allows for a significant reduction in the number of operations.In the example presented, the computation time was reduced by a factor of 13 from that needed by the full deterministic solution.Because each bottom-scale problem can be solved independently, parallel computation can be easily performed to further reduce the calculation time.The simplifications used to connect the subdomains induce some error, which can be maintained below 1% if the size of the subdomains is less than one-half the correlation length of the surface roughness.
There are several ways to further develop this method.For example, it can be applied to unstructured meshes with triangular elements and triangular subdomains.The size of the subdomains can be automatically adjusted depending on the surface parameters.The next step is now to include the occurrence of cavitation, which is mandatory when the lubrication regime is mixed.Another point to consider is transient problems appearing when the two surfaces are rough, as well as the squeeze effect.

Figure 1 .
Figure 1.Configuration of the problem.

Figure 4 .
Figure 4. Rough surface used in the simulation.

Figure 5 .
Figure 5.Comparison of the deterministic and multiscale pressure distributions.

Figure 6 .
Figure 6.Computation time as a function of the number of subdomains per side k.

Figure 7 .
Figure 7. Calculated load as a function of the number of subdomains per side k.

Figure 8 .
Figure 8.Comparison of the deterministic and top-scale pressure distributions with different values of k.

Figure 9 .
Figure 9.The different surfaces used for the accuracy analysis.

Figure 10 .
Figure 10.Relative error in the load as a function of the ratio of the surface correlation length to the sub-domain size.

Table 1 .
Geometrical configuration of the problem.