A General Finite Volume Method for the Solution of the Reynolds Lubrication Equation with a Mass-Conserving Cavitation Model

This contribution presents the development of a general discretization scheme for the solution of Reynolds equation with a mass-conserving cavitation model and its application for the numerical simulation of lubricated contacts to be discretized using irregular grids. Such scheme is based on a hybrid-type formulation, here named as element-based finite volume method that combines the flexibility of the FEM to deal with unstructured grids, while preserving the local and global fluid-flow conservation aspect of the FVM throughout the discretized domain. The accuracy and robustness of the method are successfully tested using several benchmark cases proposed in the recent literature. Simulations of fully or partially textured sliding bearings are finally employed to show the advantages of being able to adopt irregular meshes both in terms of flexibility for the discretization of complex surface features and computational speed.


Introduction
The exact close-form solution of the Reynolds lubrication equation exists only for very particular problems involving simple contact geometries, boundaries conditions and isothermal flows. Examples of such peculiar applications include short and long journal bearing components, as well as special formulations conceived for sliding and pure-squeeze bearings. However, when the evaluation of more complicated, realistic lubrication systems is desired, approximate solutions based on numerical techniques have to be considered. The numerical solution of Reynolds equation is traditionally performed by adopting the Finite Difference Method (FDM). The predominance of such discretization scheme is due to its relative simplicity in terms of formulation, along with the geometrically uncomplicated domains usually found in conventional lubrication systems, which may well be discretized using structured meshes. In this method, the partial derivatives of the equation are approximated by finite difference formulas obtained from the truncation of the high-order terms of Taylor expansions.
In contrast, in applications involving more complex geometries, when unstructured and/or irregular grids are necessary to accurately discretize the solution domain, the FDM has limited applicability. In the presence of irregular geometries and meshes, the Finite Element Method (FEM) is the most widespread method employed for the numerical solution of the Reynolds equation [1][2][3][4][5][6][7], especially due to its great flexibility to deal with distorted elements (interpolation and shape functions). However, when mass conservation is contemplated in the fluid-film cavitation modelling by the imposition of the so-called JFO (Jakobsson, Floberg and Olsson) conditions [8][9][10], the solution of the modified diffusionconvection Reynolds equation (see e.g. the − Elrod-Adams model [11][12]) is not straightforwardly accomplished with the FEM formulation. Essentially, the main difficulties arise in the discretization of the convective term of the modified equation, as well as in the enforcement of the flow conservation on the cavitation boundaries throughout the lubricated contact. Recently, a new FEM formulation based on the mathematical derivation of a linear complementary problem (LCP) has been proposed for solving the fluid-film lubrication with mass-conserving cavitation model [13][14]. Moreover, [15] have developed an efficient algorithm for fluid pressure calculation defined according to a reformulation of the complementary constraints imposed by the JFO conditions. Alternatively, the finite volume method (FVM), widely utilized in Computational Fluid Dynamic (CFD) simulations, has been proven to be very effective for solving lubrication problem, especially in thermohydrodynamic analysis [16][17], lubrication with discontinuous domains [18,19] and textured surfaces [20], as well as in piston-ring applications [21]. The main advantage of the FVM is its conservative characteristic, which in turn enables to impose local and global flow conservation in the discrete formulation. Such intrinsic conservative feature automatically enforces the complementary JFO conditions and allows a straightforward upwind-based treatment for the convective term of the modified Reynolds equation, hence facilitating the entire discretization procedure. The only disadvantage of conventional FVM schemes is the absence of a standardized strategy for the discretization of unstructured meshes, which is a prerogative of the FEM.
In this context, the main objective of the present contribution is to propose a new discretization scheme for solving Reynolds equation with the − Elrod-Adams cavitation model on irregular meshes. Furthermore, a generalization of the cavitation algorithm proposed by Ausas [20] was proposed for application to unstructured meshes. The new scheme is based on a hybrid-type formulation, hereafter named as the Element-Based Finite Volume Method (EbFVM), which combines the geometric flexibility of the FEM to deal with unstructured grids while preserving the local and global fluid-flow conservation aspect of the FVM throughout the discretized domain [22][23]. Essentially, the EbFVM is a FVM scheme that employs the concept of elements and their interpolation and shape functions adopted in the FEM. However, differently from the classical FEM methodologies, local and global conservation of the transport properties are exactly enforced in the EbFVM, since the construction of the discretized equations follows the FVM principles; this eliminates the potential issues arising from the use of the FEM mathematical foundation to derive the weak form of the formulation. In other words, in the EbFVM the discretized equations represent the physical balances over control volumes that are built as dual entities around all vertices of the primal grid, with contributions of neighboring elements. The motivation for the use of the EbFVM in lubrication applications is based on its flexibility, generality and suitability for computational implementations.
The reminder of the article is organized as follows. The mathematical formulation of the fluid-film lubrication problem (Reynolds equation) with the chosen mass-conserving cavitation model is presented in Section 2. The fundamentals of the EbFVM and their applications to the modified Reynolds equation is described in details in Section 3, where the procedure employed for the iterative solution of the discrete system of linear equations is also explained. Finally, in Section 4, a number of test cases based on classical examples taken from the recent literature are presented for validation purposes; the capabilities of the proposed formulation are then illustrated by simulating a textured sliding bearing to show the versatility and robustness of the proposed method to solve lubrication problems characterized by thin fluids films in the presence of cavitation, particularly when non-structured meshes are involved.

Mathematical Formulations
In this section, the mathematical formulations that describes the fluid-film lubrication problem is presented in details. For completeness, the most general form of the Equation of the Mechanics of Viscous Thin Films is considered [24][25], as well as its modified version obtained after the consideration of the mass-conserving Elrod-Adams cavitation ( − ) model.

Lubrication Theory (Reynolds equation)
The general Equation of the Mechanics of Viscous Thin Films (hereafter simply referred to as Reynolds equation) for isothermal, transient lubrication systems, in which the lubricant rheological properties (i.e. dynamic viscosity and density) are assumed constant across the film thickness, can be expressed in Cartesian coordinates as follows: where is the hydrodynamic pressure, ( 1 , 2 ) the geometric descriptions of the contact surfaces and and the lubricant dynamic viscosity and density, respectively. The velocities of the two mating surfaces in the and coordinate directions are designated as ( 1 , 2 ) and ( 1 , 2 ) respectively; the subscripts 1 and 2 denote bottom and top surfaces according to the reference definitions in Fig. 1.
It should be noticed at this point that in Eq. 1 the normal-squeeze velocities can only be produced by (local) temporal changes in the lubricant film thickness, i.e. ( 2 − 1 ) ≡ ( 2 − 1 ) . This assumption is generally valid for sliding and thrust bearings, but not necessarily for rotating-type systems as journal or rolling bearings. For the latter components, some modifications have to be taken into account in this equation, essentially by handling both the translation and normal squeeze terms in order to deal with the particular kinematic features of the systems. Such changes, however, by no means compromise the correctness and the applicability of the numerical scheme proposed in this work.

Mass-Conserving Cavitation Model
The cavitation (or fluid-film rupture) phenomenon plays an important role in the performance of lubrication systems and so must be included in the mathematical modelling of the problem. In order to address such cavitation effects, the mass-conserving formulation proposed by JFO (Jakobsson, Floberg and Olsson) [8][9][10] is here adopted. This cavitation formulation imposes proper complementary boundary conditions that ensure the mass-conservation principle throughout the lubricated domain. The Elrod-Adams [11][12] − model for cavitation is utilized to automatically satisfy the JFO's conditions. In this fashion, the so-called lubricant film fraction parameter is introduced into the right-hand side of Eq. 1, providing the following modified − version of the Mechanics of Viscous Thin Films: The film fraction field can be interpreted as an auxiliary variable associated to the proportion of liquid (lubricant) everywhere within the domain. The magnitude of varies in the interval [0,1], with = 1 in the pressured regions where ≥ , i.e. where the lubricant film is fully developed. However, in the cavitated zones, i.e. where 0 ≤ < 1, the lubricant film is broken (film rupture) and the medium is filled with a biphasic mixture of liquid and gases/vapors; the fluid pressure within these cavitation zones is assumed constant and equal to limit cavitation pressure .

Conservative Form of the Equation of the Mechanics of Viscous Thin Films
For the development of the EbFVM discretization scheme, it is convenient to re-write Eq. 2 in its conservative vector form, as follows: where and are the diffusivity and convective matrices. Moreover, the scalars , , and designate the translation squeeze, normal squeeze and temporal expansion source terms. Such quantities are defined as follows: Notice that the conservative form of Eq. 3 is general, in the sense that particular geometric and kinematic characteristics of different lubrication systems, as well as possible anisotropic effects on the roughness scale in mixed-lubrication applications, may be included in the formulation by easily modifying the definitions of the matrix and scalar terms of the equation.

Lubricant Rheology
The lubricant rheological properties, i.e. density and dynamic viscosity, are strongly affected by the temperature, pressure and shear-rate conditions. Since thermal effects have been neglected in the present contribution, only the isothermal density-pressure (Dowson-Higginson equation), viscosity-pressure (e.g. Barus and/or Rhoelands equation) and viscosity-shear-thinning (Eyring equation) effects will be considered for the corrections of the lubricant properties [26][27].

Fundamentals of the Element-Based Finite Volume Method (EbFVM)
Similarly to any numerical method used to approximate the solution of partial differential equations, the initial step of the EbFVM consists in the geometric discretization of the continuum domain in smaller, simple-shaped sub-regions denominated elements ( ), which are inter-connected by nodes ( ) that coincide with the vertices of the elements. The collection of these geometric entities is here denominated geometric grid and for the discretization of the Reynolds equation, it is composed by either triangular or quadrangular elements. In Fig. 2 a schematic illustration of the domain discretization and the geometric entities involved is shown. In contrast to the classical finite volume method commonly used for structured grids, where elements and control volumes are coincident ("cell-center" formulation) and the transport equations are integrated directly using the geometric grid, in the EbFVM scheme the elements are used as an auxiliary entity from which the control volumes are constructed ("cell-vertex" formulation). In this case, the control volumes constitute a secondary or computational grid, where the local conservation of the physical fluxes balances are effectively enforced. As depicted in Fig. 2, a control volume ( ) is associated to a given node ( ) of the geometric grid and is composed by portions of the neighboring elements. Those portions are designated as sub-control volumes ( ) and result from a subdivision of the elements. For a barycentric subdivision (or "median rule"), the division lines are obtained by joining the centroid of the element to the midpoints of the element edges. Those lines are denoted faces ( ) and are the basic geometric entity from which the fluxes are locally evaluated in the discretized equations. Particularly, by assuming the "midpoint rule" to approximate the fluxes over the faces, the midpoint of each face (where the fluxes are effectively calculated) and its respective normal vector are defined, respectively, as integration point ( ) and normal vector (⃗⃗ ). Following this strategy, each control volume associated to a given node may be seen as a sum of the nearest sub-control volumes of the elements that share , and the fluxes at one specific integration point can be calculated using data from the element in which the integration point is placed [28][29].
As already stated, due to the "cell-vertex" nature of the EbFVM, the unknowns of the problem are approximated at the nodes of the geometric grid. Furthermore, in order to deal with the distorted elements of irregular grids, as well as local variations in the flow and transport properties, pre-defined families of interpolation and shape functions can be employed. These functions are based on the definition of parametric elements that are described with respect to a local (transformed) coordinate system. Thus, all the required calculations previously undertaken in the physical domain can be performed more easily in the transformed domain. The reader should refer to Appendix 1 for details about the coordinate transformation and the associated shape functions of triangular and quadrangular elements.
Once the integration points are identified within each element, the fluxes of the transport properties through the faces of the sub-control volumes can be calculated according to the geometric and general parameters associated with the element nodes. This procedure allows elements, and so the geometric grid, to be treated as the major entity over which all the calculations are performed independently. Subsequently, the conservation equation for every control volume can be determined from the fluxes contributions coming from its respective communal surrounding sub-control volumes. Thus, after an element-by-element assembly of the global system of equations (analogous to the assembly procedure in the FEM), the total physical balance of the conservation equation are automatically satisfied by the inherent connectivity of the elements.
In the next sections, the temporal and spatial discretization of each term of the modified Reynolds equation (Eq. 3a) will be presented in details. Further explanations concerning the EbFVM scheme for a generic conservation equation can be found in [22][23]28].

Temporal Discretization
The temporal approximation of Eq. 3a can be established by integrating all terms of the equation in a time interval Δ , as follows: Assuming an implicit first-order time discretization scheme (backward Euler method), the integrands of Eq. 5 can be approximated to their values at the end of the time interval. The advantage of such formulation is its unconditional stability, so that the numerical accuracy of the solution is mainly controlled by the size of the time step Δ . Thus, by performing the temporal discretization using = Δ , where is the n-th time instant, Eq. 5 can be rewritten as:

Spatial Discretization
For the spatial discretization of Eq. 6 at each time instant = Δ , it is convenient to consider its integral formulation (or weak formulation), which may be obtained by adopting the Weighted Residual Method (WRM). The reason for using such integral approach is to ensure, mathematically, the minimization of the approximation errors throughout the solution domain [30]. In this case, Eq. 6 is expressed as follows: where corresponds to the weighting function that weights the discretization errors over the entire integration domain ( * ), and the superscript '~' denotes the approximated solutions of the transport properties and unknown variables. Particularly for finite volume methods, is assumed unitary within every control volume and null outside (sub-domain method) [31]. Thus, Eq. 7 can be restated for each control volume as: where is the sub-domain enclosed by each control volume centered at every node .
By applying the divergence (or Gauss-Ostrogradsky's) theorem to the first two integrals of Eq. 8, one obtains: where is the boundary that delimits each control volume and ⃗ ⃗ the unit vector normal to that boundary (see Fig. 2).
Finally, recalling that the fluxes over the faces are approximated by the "midpoint rule" and by assuming a second-order approximation scheme for the area integrals (in which the integral is evaluated as the product between the integrant function at the center and the area [32]), Eq. 9 can be expressed in its discretized form as follows: where ⃗⃗ * = ⃗⃗ Δ is the normal face vector that retains, simultaneously, the normal orientation (⃗⃗ ) and the length (Δ ) of each delimiting face associated to a given control volume; Δ is the area of the s-th that compound the control volume and Δ is the full area of the respective (see Appendix 1 for details about the computation of ⃗ ⃗ * for triangular and quadrangular elements). Moreover, in Eq. 10 ℰ denotes the set of elements that contribute to the formation of the around node and ℱ the set of faces inside the corresponding element that encloses .
The above discrete equation is valid for every control volume constructed surrounding each node of the geometry grid. The inner summations in the first two terms of the equation correspond, respectively, to the diffusive and convective fluxes through each sub-control volume, while the outer summations represent the entire flux balances over the associated . The last two components of the equation are the source terms evaluated explicitly at node .
As already pointed out, the assembly of the global system of equations is performed in an element-byelement fashion, similar to the procedure employed in the FEM. Furthermore, all local computations are carried out at element level using the transformed domain approach (parametric elements); this allows the elements to be treated independently, no matter how distorted they are with respect to the global coordinate system. Such procedure renders the definition of element matrices and vectors constituting each term of Eq. 10 straightforward, analogously to the stiffness and mass matrices and load vectors used in FEM formulations. In the following, those definitions will be exposed in details by adopting the index notation for compactness.

Diffusivity Matrix (Poiseuille Term)
The total diffusive flux (first term of Eq. 10) across the sub-control volumes of a given element can be represented as: where the components of the vector ⃗ ⃗ correspond to the diffusive fluxes through each sub-control volume, (⃗ ⃗ ) is the vector of the nodal pressures and is the diffusivity matrix of the element. See Appendix 2 for details about the derivation of .

Convective Matrix (Couette Term)
Similar to the diffusion term, the total convective flux (second term of Eq. 10) flowing through the subcontrol volumes of a given element can be expressed as: where the entries of the vector ⃗ are the convective fluxes through each sub-control volume, (⃗ ⃗ ) the vector of the nodal lubricant film fraction and is the convective matrix of the element. See Appendix 3 for details about the derivation of .

Source Term Vectors
The source terms (third and fourth terms of Eq. 10) associated to the sub-control volumes of a given element can be written as follows: where the members of the vectors ( ⃗ ) and ( ⃗ ) are the source terms of each element node, while  1), respectively. In the above equations, the operator " ⋅ " denotes scalar product.

Assembly of the Global Linear System of Equations
Subsequently to the calculation of the diffusivity and convective matrices, as well as the source term vectors for all elements of the geometric grid, the discrete equation associated with every control volume are assembled in a global system of linear equations. The rows of the resulting linear system correspond to the entire flux balance over each control volume constructed around every node of the grid. This can be mathematically represented as: where denotes the neighboring nodes around and the coefficients and are the entries of the linear system of equations computed during the assembly operation.

Solution of the Linear System of Equations
For each time instant , the global system of equations defined in the previous section (Eq. 14), embraces two unknown variables, namely the hydrodynamic pressure (̃), and the lubricant film fraction (̃) fields. Thus, a specific numerical procedure has to be considered for the simultaneous solution of ̃ and ̃. In the present contribution, such simultaneous calculations are accomplished by adopting a generalization of the cavitation algorithm proposed by [20], which was originally conceived for applications with structured grids. In this sense, the generalization here devised consists in an extension of the original algorithm for lubrication problems to be discretized by irregular meshes.
The main idea of the algorithm is based on the iterative solution of linear system of equations according to the well-known Gauss-Seidel method with Successive Over-Relaxation (SOR). In this case, by assuming that , and are known a priori, the nodal values of ̃ and ̃ can be computed iteratively, as follows: where is the r-th iteration of the SOR method, and the respective relaxation factors, and ̃, and Θ , the intermediate nodal values of ̃ and ̃.
The intermediate quantities ̃, and Θ , are approximated from Eq. 14 by isolating the respective nodal values and admitting the others fixed to the values calculated in the previous iteration, which mathematically corresponds to: Furthermore, the complementary conditions related to the − cavitation model displayed in Eq. 3b must be satisfied during the iterative procedure. In the discrete domain this reads: The discrete complementary conditions of Eq. 17 enforce the cavitation regions and boundaries to be positioned always at the nodes of the geometric grid. Moreover, as the discretization method is inherently conservative in the discrete domain (finite volume method), the local conservation of the lubricant flow is ensured locally for all grid nodes (control volumes), including those laying on the cavitation boundaries; this automatically satisfies the mass-conservation conditions imposed by the JFO cavitation model. The extended SOR algorithm proposed in this contribution is shown as a pseudo-code in Appendix 5.

Results and Discussion
In this section, the results of several simulations cases performed using the EbFVM discretization scheme and the iterative algorithm described above are presented. Initially, examples from the literature have been considered for validation purposes, followed by simulations of a textured sliding bearing with different dimples' distribution aimed to illustrate the flexibility of the method for dealing with irregular meshes.

Validation Cases
The example cases considered in this section have been chosen to test the correctness and accuracy of the proposed formulation and algorithmic implementation by comparison with other solution methodologies published in the recent literature, including both alternative mass-conserving Reynolds' solvers and full CFD simulations. Each considered example has been selected to demonstrate the various features considered in the newly developed computational framework. It should be emphasized here that the main focus of the work presented by the authors is the validation of the newly developed discretization scheme rather than the direct comparison of performance between different algorithms proposed in the literature.

o Single and double parabolic slider with density-pressure correction
This first example has been chosen to evaluate the effectiveness of the proposed EbFVM formulation in comparison with different methodologies for solving lubrication problems involving compressible fluids. To do so, two simulation cases characterized by single and double parabolic sliders similar to those tested in [33] have been analyzed. However, as proposed in [14], the lubricant film geometry has been modified with respect to that used in [33] in order to highlight the effects of the fluid compressibility. The main geometric and operational parameters of each parabolic slider are summarized in Table 1. The results for the single parabolic slider are examined against those calculated using the formulations proposed in [33], the linear complementary problem (LCP) solution in [14] and the incompressible case. Figure 3a confirms the agreement of the current EbFVM formulation with those obtained with the algorithms proposed by Sahlin et al. [33] and Bertocchi et al. [14], both in terms of hydrodynamic pressures (maximum difference in peak pressure of 0.76%) and location of the cavitation boundary separating the pressured and cavitated zones.
Additionally, the double parabolic slider case also proposed by [33] is used to assess the capability of the EbFVM for predicting the fluid-film reformation boundary. The results are compared both with respect to a CFD analysis based on the full solution of the Navier-Stokes equations developed using the OpenFOAM code [34], as well as with the same methods employed by Sahlin et al. [33] and Bertocchi et al. [14]. The consideration of CFD solutions as a comparative basis is aimed to (1) ensure the validity of Reynold equation for flow problems which satisfy the fundamental lubrication hypotheses, and further (2) to demonstrate that in common lubrication applications the intricate flow behaviour of the liquid-vapour/gases mixture within the cavitation zones is in general meaningless for the overall predictions of the system performance. Figure 3b illustrates an excellent agreement between the results calculated with the current EbFVM scheme compared to those obtained from the CFD simulations and using the methodologies by Sahlin et al. [33] and Bertocchi et al. [14]; the maximum deviation for the hydrodynamic peak pressure is ≤0.5% and the rupture and reformation cavitation boundaries are practically identical.  o Long journal bearing with viscosity-pressure correction In this second example, the capability of the proposed EbFVM formulation for solving problems with piezoviscous lubricants is tested. Particularly, the simulation case of an infinitely long journal bearing solved in [14] based on the LCP approach for an incompressible piezoviscous lubricant flow is considered. The main problem parameters are listed in Table 2. Two different operational conditions represented by the fixed eccentricity ratios of = 0.93 and = 0.95 are investigated. The Barus equation is adopted for the piezo-viscosity correction. Figure 4 compares the hydrodynamic pressure profiles along the bearing circumferential directions computed with the current EbFVM formulation against the results obtaining using the LCP method [14] for the two eccentricity ratios considered. Again, the results are in good agreement in both cases, with a maximum peak pressure difference of 0.9%. Excellent agreement is observed for the cavitation boundaries as well.  o Pure squeeze motion between circular plates In this example case, the validation of the proposed EbFVM scheme is assessed for a twodimensional domain under time dependent conditions. The classical lubrication problem of two circular plates with finite radius and under pure (sinusoidal) squeeze motion is considered. Moreover, any viscosity and density corrections are ignored. The results are compared with those obtained from the analytical solution developed by [3] and the linear complementary method of [14]. The main system parameters are displayed in Table 3. For the numerical simulation, an irregular grid with 2913 nodes has been used, and 3 oscillating periods with 576 times steps per cycle are adopted for the time solution. Figure 5 depicts the time variation of the extension of the cavitation zone calculated from the above-described methodologies chosen for comparison; notice that the internal "annulus-like" surface plots correspond to the simulated film fraction fields at 3 different times, which illustrate the progressive variation of the cavitation zone (in blue). Once again, the agreement of the results is excellent.  o Sliding pocket bearing with viscosity-pressure, viscosity-shear-thinning and density-pressure correction This second two-dimensional example deals with a sliding pocket bearing lubricated with a compressible, piezoviscous and shear-thinning fluid. Two geometric configurations characterizing a finite and an infinitely long bearing are investigated. All the geometric and operational parameters are listed in Table 4. Figure 6 compares the pressure profiles along the bearing mid-section calculated with the current EbFVM scheme and with the LCP formulation by Bertocchi et al. [14]. As can be observed, the results computed with the new finite volume scheme match very well those obtained with LCP technique for both bearing geometries, including the analytical solution for the infinitely long bearing case.

Textured Sliding Bearing Simulation
In order to illustrate the robustness and flexibility of the EbFVM to tackle irregular and dense meshes, a set of simulations for a parallel and flat textured sliding bearing have been undertaken for the geometric and operational parameters listed in Table 5. Full and partial texturing are considered by varying the number of dimples' cells added onto the stationary bearing surface from 10x2 to 10x10 cells; the dimples are spherical in shape and cover approximately 20% of each texture cell (see details of the dimples' dimensions in Table 5). Moreover, the lubricant density and viscosity are assumed constant as well.  The lubricated domains are discretized using both regular and irregular grids; in particular, the latter grid type are used for allowing mesh refinements in the neighborhoods of the texture features, where the pressure gradients are more pronounced. Such meshing control may well be used to reduce the number of nodes, and so the computational cost of simulations, needed to discretize the same textured pad when compared to structured grids. Figure 7 illustrates both the regular and irregular meshes of 3 representative texture variants with 10x10, 10x6 and 10x2 dimples cells, as well as the associated hydrodynamic pressure fields calculated; notice the considerable decrease in the number of nodes for the unstructured meshes.
The number of outer iterations of the extended SOR algorithm (see Appendix 5), the computational time for convergence and the hydrodynamic load capacity are employed for comparing the results calculated on the regular and irregular meshes. The optimum relaxation factors of the SOR method for pressure ( ) have been found empirically by numerical experimentation. The optimum values for each texture pattern and mesh topology are displayed in Table 6; the relaxation value for the lubricant film fraction (cavitation) is assumed fixed as = 1. All the simulations were carried out in MATLAB in a computer with 8 GB RAM and Intel Core i7-3630 CPU 2.40 GHz; the algorithm reported in Appendix 5 is implemented as a MATLAB MEX-File that calls a compiled C function for performance improvement. As can be verified in Table 6, the variation on the mesh topology hardly changed the overall load carrying capacity for all texture cases examined, but led to a significant speed-up in the simulation time, yielded essentially by the decrease in the number of degrees of freedom (DOFs) for the irregular grids. Such occurrence demonstrates the enhancements that can be achieved by using the current finite volume method; this is a result of the newly acquired capability of deploying non-structured meshes to simulate lubricated contact, which in turn can contribute either to reduce the number of DOFs in problems with complex geometries or to improve, by mesh refinements, the solution accuracy near regions characterized by special surface features in the contact domain. Although the absolute computational efficiency of the method is not the main concern at this point, since the idea of this numerical experiment is to illustrate how irregular meshes can be used to minimize the simulation time of lubrication problem with textured surfaces, the current formulation yields simulations times comparable to those reported for similar textured bearings in a recent publication [15].
We believe that performing a thorough and rigorous quantitative comparison between different solvers and methods in term of computational cost is a very important task. However, given the complexity of such task for extracting and quantifying the individual contribution to efficiency provided by cavitation algorithms, discretization schemes, hardware and code optimization methodologies followed by different researchers, this is deemed outside the scope of the present work and is the subject of the authors' current and future research efforts. At this stage it may be worth mentioning that the most efficient solvers can potentially be obtained combining reliable discretization schemes, such as the one proposed here for unstructured grids, and efficient solution methods (see e.g. [15]).
Finally, Figure 8 depicts a detailed view of the hydrodynamic pressure distributions calculated in on the regular and irregular meshes for the same texture variant. As expected, no noticeable difference is observed between the pressure fields. It is important to remark that the results displayed in Table 6, particularly the number of iterations, do not follow a monotonic pattern as the number of DOFs is reduced; this is somewhat counter-intuitive but is solely due to the fact that each dimple variant used to generate the results reproduced in Figure 7 corresponds to a different geometric domain. In fact, the monotonic pattern is expected be observed if the contact geometry (number of dimples rows) is kept constant and only the mesh topology is allowed to change. This has been here evaluated through a set of simulations for the 10x6 texture variant and regular meshes. In this case, each dimple cell was discretized locally by different grid densities, i.e. 20x20, 30x30, 40x40, 50x50, 60x60. The overall results are summarized in the Table 7. As can be seen, as the total number of mesh nodes increases, the amount of interaction required for convergence also increases accordingly.

Conclusion
This paper introduces Element-Based Finite Volume Method as a new paradigm in the simulation of lubrication problems characterized by thin fluid films in the presence of cavitation. The EbFVM algorithm developed to solve the conservative form of the equation of the mechanics of viscous thin films has been tested and successfully validated by comparison with existing alternative algorithms and full CFD simulations. The test cases presented have been used to demonstrate that the numerical approach proposed is robust and versatile and can incorporate compressibility, piezoviscosity and shear-thinning of the lubricant, while being applicable to any transient hydrodynamic lubrication problem in the presence of cavitation. By successfully combining the flexibility of the FEM to deal with irregular geometries and the inherent conservative nature of the FVM, one of the main features of the EbFVM is that it allows the use of irregular meshes to discretize geometrically complex domain. This is shown to confer enormous flexibility to the method and enables a significant reduction in the computational requirements for solving lubrication problems in the presence of textured surfaces. Thus, the proposed method offers a new promising route for the improvement and design optimization of tribological surfaces.

Appendix 1 -Elements for EbFVM Discretization
In the following, all the geometric and tensor properties of triangular and quadrangular elements are presented for the completeness of the EbFVM derivation. The properties are shown with respect to the standard transformed domain (parametric elements). (A1.2.9)