An EHL Extension of the Unsteady FBNS Algorithm

Many engineering applications rely on lubricated gaps where the hydrodynamic pressure distribution is influenced by cavitation phenomena and elastic deformations. To obtain details about the conditions within the lubricated gap, solvers are required that can model cavitation and elastic deformation effects efficiently when a large amount of discretization cells is employed. The presented unsteady EHL-FBNS solver can compute the solution of such large problems that require the consideration of both mass-conserving cavitation and elastic deformation. The execution time of the presented algorithm scales almost with Nlog(N)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\log (N)$$\end{document} where N is the number of computational grid points. A detailed description of the algorithm and the discretized equations is presented. The MATLAB©\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\copyright }$$\end{document} code is provided in the supplements along with a maintained version on GitHub to encourage its usage and further development. The output of the solver is compared to and validated with analytical, simulated, and experimental results from the literature to provide a detailed comparison of different discretization schemes of the Couette term in presence of gap height discontinuities. As a final result, the most favorable scheme is identified for the unsteady study of surface textures in ball-on-disc tribometers under EHL conditions.


Introduction
In the case of lubrication flows in narrow gaps, the Reynolds equation [1] is a handy tool to determine the hydrodynamic pressure distribution in a simpler way than by using the full Navier-Stokes Equations (NSE) [2,Ch. 7]. Since cavitation commonly occurs in lubrication flows, various models have been developed to describe this phenomenon [3]. Especially when it occurs within surface textures, massconserving properties of the cavitation model are required to properly describe the flow's transition from the cavitation region to the full-film region, because this full-film reformulation interface has a great effect on the extension of the cavitated area and the subsequent downstream rise in pressure within the full-film region [4]. The required mass-conserving properties can be taken into account with the Jakobsson-Floberg-Olsson (JFO) [5,6] cavitation model [3]. Starting from the cavitation algorithm of Elrod [7], Giacopini et al. [8] developed a one-dimensional finite element method (FEM) solver that couples the Reynolds equation with the mass-conserving JFO cavitation model through a complementarity formulation. This work was extended by Bertocchi et al. [9] to consider two-dimensional problems with compressible, piezoviscous, and shear-thinning fluid behavior. The arising complementarity problem was reformulated to be expressed by an unconstrained equation system by Woloszynski et al. [10], resulting in the Fischer-Burmeister-Newton-Schur (FBNS) algorithm. As demonstrated by Woloszynski et al., the FBNS algorithm is of remarkable computational efficiency also for high spatial resolutions.
In many cases, the hydrodynamic pressure can deform the lubricated surfaces notably leading to the regime of Elasto-Hydrodynamic Lubrication (EHL) [11]. Various solvers have been developed to tackle EHL problems, some of the most prominent ones are the finite difference method (FDM) Multigrid solver of Venner and Lubrecht [12] and the FEM solver of Habchi [13]. Some algorithms are also capable of simulating surface contact along with the Reynolds equation 80 Page 2 of 25 [14][15][16][17][18]. Since the full-film reformulation interface is often not of relevance in EHL problems, many EHL solvers do not employ mass-conserving cavitation models. However, in some cases -such as starved lubrication-mass-conserving cavitation is crucial and has been considered in several works [19,20]. Among them, the coupling of pressure, mass-conserving cavitation, elastic deformation, a roughness asperity contact model and the FBNS algorithm was achieved by Ferretti [21,22]. In contrast to Ferretti's work, the FBNS algorithm is coupled with the elastic deformation of an elastic half-space within this paper, thus presenting the new EHL-FBNS algorithm. Due to the half-space assumption, the elastic deformation is a linear convolution of a kernel function with the hydrodynamic pressure field. This allows exploiting the fast Fourier transformation (FFT) to speed up the computation of the elastic deformation [23]. Furthermore, a proportional integral derivative (PID) controller is employed to meet the load balance equation through adjustment of the rigid body displacement as already introduced by Wang et al. [24]. Eventually, the EHL-FBNS algorithm is capable of efficiently computing the solution of large problems that require the consideration of both mass-conserving cavitation and elastic deformation at the same time.
In the beginning of this paper, the basic equations and the general procedure of the EHL-FBNS algorithm are summarized. The algorithm is implemented in MATLAB © with the finite volume method (FVM) and a generic order spatial discretization scheme for the Couette term of the Reynolds equation. The discretized equations are supplied in Appendix. Then, the performance of the steady EHL-FBNS implementation is compared to the original FBNS algorithm of Woloszynski et al. [10]. Afterward, one-and two-dimensional literature reference cases [9,25] of a convergent slider with rectangular pocket are used to compare the EHL-FBNS output to analytical and simulated reference results, assess the influence of the discretization order of the Couette term in the presence of gap height discontinuities through grid convergence studies and give an example case where both mass-conserving cavitation and elastic properties of the solver are required. Moreover, unsteady EHL-FBNS simulations are performed for the set-up of a single texture that passes through the EHL contact of a ball-on-disc tribometer. The results are compared to experimental and simulated data of Mourier et al. [26]. This allows to demonstrate the stability of the EHL-FBNS algorithm under EHL operating conditions with discontinuous surface textures and to provide recommendations about the most suitable discretization scheme. Eventually, the EHL-FBNS algorithm is validated for surface texture investigations with ball-on-disc tribometers under unsteady EHL operating conditions. The MATLAB © code, set-up, and visualization scripts are provided in the supplements. The MATLAB © scripts are thoroughly commented to encourage their usage and further development. A maintained and publicly available version of the code can also be found on GitHub: https:// github. com/ ErikH ansen Git/ EHL.

Governing Equations
The lubrication flow in a narrow gap (schematically depicted in Fig. 1) is governed by the Reynolds equation considering mass-conserving cavitation with the JFO model [5,6] at any set of spatial coordinates x 1 and x 2 and time t [9]: In this equation, h denotes the gap height, l is the density of the liquid phase and l describes the dynamic viscosity of is composed of U up as the velocity of the upper surface and U low as the velocity of the lower surface in x 1 -direction. The relative pressure p = p hd − p cav and the cavity fraction = 1 − l are the solution variables, where is the mixture density of the flow. The hydrodynamic pressure p hd is prevented from falling below the cavitation pressure p cav by adding the following complementary constraints [9]: Depending on whether the liquid phase is modeled as isoor piezoviscous through the Barus or Roelands model, its dynamic viscosity reads [12,Ch. 1.3.3], [27]: [26] is the pressure viscosity index, B and R denote the pressure viscosity coefficient and p 0,R is a constant in the Roelands equation. The dynamic viscosity of the liquid phase at ambient pressure is 0 . Moreover, depending on whether the liquid phase is assumed to be of constant density or to be compressible according to the Dowson-Higginson model, the liquid phase density is given by [12,Ch. 1.3.4], [9,27]: where 0 is the density of the liquid phase at ambient pressure and C 1 and C 2 are constants.
The gap height h can be constructed as a superposition of the rigid body displacement of the two surfaces h d , the variation of the gap height due to the rigid geometry of the surfaces h g and the elastic deformation of the gap height h el due to the hydrodynamic pressure [2,Ch. 19.2]: Depending on whether the upper and lower surfaces are assumed to be rigid or elastic half-spaces, the elastic deformation of the gap height can be expressed as [

EHL-FBNS Algorithm
The set of equations described above can be solved numerically with the EHL-FBNS algorithm presented in the following. This new algorithm is based on the FBNS algorithm developed by Woloszynski et al. [10] and extends it by taking elastic surface deformation and the load balance equation into account. First of all, the dimensionless Reynolds equation considering mass-conserving cavitation is defined as:  (9) and (12) are provided in Appendix A.1. The set of discrete dimensionless Reynolds equations ⃗ G can be expressed in matrix-vector notation through the pressure coefficients contributed by the Poiseuille term A Po , the discretized dimensionless relative pressures ⃗ p * , the cavity fraction coefficients contributed by the Couette and unsteady term B, the discretized cavity fractions ⃗ , and the remaining constants from the Couette and unsteady term ⃗ c [10]: The set of discrete dimensionless Fischer-Burmeister Equations (12) are denoted by ⃗ F( ⃗ p * , ⃗ ) . The non-dimensional properties * and * at each discrete point are computed according to the respective Eqs. (3, 4 and 10).
The gap height h at each discrete point is computed according to Eqs. (5 and 6). It is prevented from becoming lower than 1 nm by using truncation at this instant. (10) Afterward, the non-dimensional gap height h * is determined through Equation (10). If the surfaces are chosen to be elastic, Equation (6) is discretized by assuming a constant pressure over the rectangular discretization cell [30,Ch. 3.3], [31,Ch. 3.1], the discretized equation is also provided in Appendix A.2. Since the resulting equation is a linear convolution of a kernel function with the hydrodynamic pressure field, it is computed in Fourier space by means of FFT to speed up the computation. Attention is paid to double the size of the kernel in each direction and to zero pad the hydrodynamic pressure field such that a linear instead of a circular convolution is obtained. After the convolution, the deformation and pressure fields are resized to their original size [23,32,33]. After computing ⃗ G and ⃗ F , the Newton-Raphson method is used to determine the updates of non-dimensional relative pressure ⃗ p * and cavity fraction ⃗ [10]: The most important extension of the EHL-FBNS algorithm compared to the original FBNS algorithm is the approximation of the pressure Jacobian J G,p * of the dimensionless Reynolds equation when elastic deformation is taken into account. The idea is to consider the dependence h * (p * ) by inserting it into ⃗ c , thus creating the matrix A h . Due to the kernel function, this would result in A h being a full matrix which is prone to lose its diagonal dominance and therefore being unfeasible to invert and likely to cause unstable behavior in the iteration process. This is rectified by approximating A h only by some of its diagonals as already done in the literature for other EHL algorithms: for example, Venner and Lubrecht [12, Chs. C.1, C.3.2] who combine it with distributive relaxation and multigrid methods or Wang et al. [15] who employ the semi-system method. In case of the EHL-FBNS algorithm, A h is reduced to the 5 diagonals that correspond to the South, West, Center, East, and North cells. Eventually, the Jacobians of ⃗ G read J G,p * = A Po + A h and J G, = B . The boundary conditions of ⃗ p * are considered in A Po and ⃗ c and the boundary conditions of ⃗ in ⃗ F and J F, . If Neumann boundary conditions are used for ⃗ , the Jacobian J F, would contain several diagonals. In this case, it is approximated only by its main diagonal. It is worthwhile to note that this approximation of the Jacobians eventually only affects the updates of non-dimensional relative pressure ⃗ p * and cavity fraction ⃗ , but never the computation of ⃗ G and ⃗ F .
The discrete formulations of the Jacobians J G,p * = A Po + A h and J G, = B are provided in Appendix A.1 and A.3. The center entries of the Jacobians of the dimensionless Fischer-Burmeister equation J F,p * ,C and J F, ,C for each discrete point are [10]: Here, p * C,aux and * C,aux are the auxiliary dimensionless pressure and cavity fraction which are adjusted such that J F,p * and J F, do not become singular [10]. To prevent them from having center entries close to zero within the range (− , ) , they are adjusted as: where machine epsilon is given by ≈ 2.2204 ⋅ 10 −16 . As already done in the original FBNS algorithm, the corresponding columns of the Jacobian J and rows of the updates ⃗ are swapped if J F,p * ,C < J F, ,C to obtain a reordered system [10]: Due to the swapping, A F is better conditioned than J F,p * which is exploited when Equation system (19) is solved [10]: After obtaining ⃗ a and ⃗ b , the earlier performed row swapping is reversed to get the updates of non-dimensional pressure ⃗ n p * and cavity fraction ⃗ n [10]. The new values ⃗ p * ,n and ⃗ n at iteration n are obtained by means of relaxation: where p * , ⃗ p * ,n−1 , and ⃗ n−1 are the relaxation factors and previous solutions of non-dimensional relative pressure and cavity fraction. Depending on the simulated case, relaxation coefficients between 0.05 and 1 resulted in good tradeoffs between convergence speed and stability. Preventing ⃗ p * ,n from having values below 0 and ⃗ n from having values below 0 or above 1 through truncation furthermore enhances favorable convergence properties.
If a constant load force is prescribed, the dimensionless rigid body displacement h * d = h d ∕h ref is adjusted through a PID controller to meet the load balance Equation (8) as already done by Wang et al. [24]. This is done by first determining the resulting normal load force F n N through the discretized load balance equation: where N x 1 , Δx 1 , N x 2 , and Δx 2 are the amount and spacing of the discretization cells in x 1 -and x 2 -direction and p n hd,C is the hydrodynamic pressure at the center of each discrete cell. The residual of the load balance equation is defined as: where F N,ref is a reference normal force that is usually just set equal to the imposed normal force F N,imp . Note that r n F N can be either positive or negative, depending on whether F n N is larger or smaller than F N,imp . This is required for the PID controller to work properly. Finally, r n F N is fed into the PID controller to determine h * ,n+1 d of the next iteration step [24]: Note that K I is only multiplied with the sum up until r n−1 F N , since r n F N is already considered by K P . For all of the later considered simulations with imposed normal load force, K P = 0.001 , K I = 0.02 , and K D = 0.001 worked well. At last, the following residuals are computed: Note that the residuals r n max, G and r n max, F are directly affected by the relaxation factors and ⃗ G n and ⃗ F n are computed through the solutions ⃗ p * ,n−1 and ⃗ n−1 . The EHL-FBNS algorithm is repeated as long as r n EHL-FBNS and in case of an imposed normal load force abs(r n F N ) are above the tolerance tol = 10 −6 .
(28) r n EHL-FBNS = max r n max, p * , r n max, , r n max, G , r n max,G , r n max, F , r n max,F .

Swap columns of J and rows of δ
Compute δ a and δ b Reswap rows of δ a and δ b to obtain δ n p * and δ n θ Compute p * ,n and θ n , truncate p * ,n below 0 and θ n outside of 0 and 1 Yes Yes Page 7 of 25 80 The most important steps of the algorithm structure are also visualized in Fig. 2. The initial guess is always a zero cavity fraction field and a pressure field at ambient pressure. If an unsteady simulation is performed, the solution at t = 0 is obtained through the steady problem caused by the geometry at t = 0 . Furthermore, the PID controller also takes the residuals of the load balance equation of the previous time steps into account if t > 0.

Results and Discussion
In order to assess the performance of the presented EHL-FBNS algorithm, it is firstly employed in a numerical literature test case of a textured parallel slider. Then, the results of the EHL-FBNS algorithm are compared to the analytical solution of a rigid one-dimensional convergent slider with rectangular pocket to show the effect of the discretization order of the Couette term on the accuracy of the simulation result when gap height discontinuities are present. Subsequently, the slider is extended to a two-dimensional geometry and an elastic model is employed to give an example case where both mass-conserving cavitation and elasticity show relevant effects. Afterward, another experimental-numerical literature case is simulated with the EHL-FBNS algorithm to validate the code for textured ball-on-disc investigations, evaluate the code's stability in unsteady EHL conditions, and compare different spatial discretization schemes. For all considered cases, the second-order midpoint rule is used for the evaluation of integrals arising from the FVM and the Poiseuille term is always discretized with a second-order central scheme. Consequently, the resulting order of the dimensionless Reynolds equation discretized with the FVM in the steady case is first order with the UI and second order with the QUICK scheme. In the unsteady case, only first order is achievable for both UI and QUICK since the firstorder Euler implicit scheme is employed. All of the EHL-FBNS simulations are performed with MATLAB © R2020a.

Parallel Slider with a Varying Amount of Trapezoidal Pockets
The parallel slider with a varying amount of trapezoidal pockets used by Woloszynski et al. [10] serves as first test case. This set-up is chosen because it can cause a generic amount of distinctive cavitation regions. This is to demonstrate the good performance and stability properties of the EHL-FBNS algorithm since the simulation of such cases showed to be difficult or inefficient with other codes from the literature. Furthermore, the solid properties are chosen such that noticeable effects due to elastic deformations on the pressure profile occur. Thereby, it is shown that the consideration of elasticity does not alter the performance scaling. While being of little physical interest, this numerical set-up allows a comparison of the performance scaling to the original FBNS algorithm of Woloszynski et al.
The variation of the gap height due to the rigid geometry of the surfaces h g is constructed by assembling several unit geometries. Each unit geometry is composed of 20 ⋅ 20 cells with h g = h min + h p . This square is surrounded by a one cell thick layer with h g = h min + h p ∕2 , which is in turn surrounded by a layer of four cells with h g = h min , resulting in a square of 30 ⋅ 30 cells. Depending on the desired size of the computational domain, a certain amount K x 1 ⋅ K x 2 of these unit cells is attached to each other and finally surrounded by a layer of cells for the Dirichlet boundary conditions of hydrodynamic pressure at p amb and zero cavity fraction. At last, the coordinates of each cell center are set such they are in the range of . The resulting array of cells in the exemplary case of K x 1 = K x 2 = 2 is visualized in Fig. 3. The rigid body displacement h d is set to 0 since it is already considered in h g .
The values of the parameters used in the EHL-FBNS simulations are summarized in Table 1. Piezoviscosity and compressibility of the liquid phase are not considered. The  Table 2 Summary of the code execution times of the FBNS algorithm given by Woloszynski et al. [10] and the ones of the EHL-FBNS algorithm with UI or QUICK discretization and rigid or elastic geometry to simulate the parallel slider with a various amount of trapezoidal pockets The FBNS simulations of Woloszynski et al. were performed on a workstation with 32 GB RAM and an Intel Xeon 3.3 GHz processor while the computations of the EHL-FBNS algorithm were conducted on a workstation with a 64 GB RAM and an AMD Ryzen 9 3900X 12-Core 3.8 GHz processor is always enforced such that the resulting mesh is always quadratic. By increasing the amount of unit geometries, the total amount of discretization cells N is increased. Each resulting geometry is simulated with UI and QUICK for the rigid and elastic case. In the elastic case, the solid bodies' Young's modulus E and Poisson ratio are set such that the elastic deformation shows notable effects even though the overall pressures are low. At the same time, this set-up produces many distinctive cavitation regions. The rigid simulations use relaxation factors of = 1 . Since the elastic simulations are generally more unstable, the relaxation factors have to be reduced to 0.5.
To quantify the algorithm's performance independently of the hardware's computational power, the non-dimensional code execution time is defined as: where t ex is the code execution time for a certain total amount of discretization cells N. The non-dimensional code execution time t * ex of the EHL-FBNS algorithm in the rigid UI case is compared with the one of the original FBNS algorithm of Woloszynski et al. [10] in Fig. 4. Their code was also implemented in MATLAB © , considered rigid geometries and was discretized with the FVM, where the Poiseuille term was discretized with second-order central differences and a first-order upwind scheme was employed for the cavity fraction. While Woloszynski et al. [10] use different tolerances for different residuals ( 10 −3 for r n max,G and 10 −6 for r n max, G ), the EHL-FBNS algorithm employs an even stricter convergence criterion ( 10 −6 for r n EHL-FBNS given by Equation (28)). Furthermore, three reference curves for the scaling of the non-dimensional code execution time are displayed: linear t * ex,lin = M 1 N , logarithmic t * ex,log = M 2 N log(N) and quad- It can be seen that the original FBNS algorithm has a performance scaling close to the linear reference while the EHL-FBNS algorithm performs a little bit slower than the N log(N) reference for large N but is always much faster than the quadratic reference. The difference between the FBNS and EHL-FBNS performances might be due to the fact that the EHL-FBNS algorithm constructs the matrices A Po and B and vector ⃗ c at each iteration step, while the FBNS algorithm might exploit that they are constant in a rigid and isoviscous simulation with an incompressible liquid phase. However, the exact details of the implementation of the original FBNS algorithm are unknown and the difference in time scaling cannot be pinned down rigorously.
Next, it is investigated how combinations of UI or QUICK discretization and rigid or elastic geometry affect the performance. The code execution times t ex are provided in Table 2. The corresponding non-dimensional code execution times t * ex are displayed in Fig. 5. Because all curves have almost the same inclinations in the double logarithmic diagram, it can be deduced that the performance scaling of the EHL-FBNS algorithm stays similar in all cases. To prove that the operating conditions were chosen such that discretization scheme and elastic model have a noticeable impact on the results while the performance scaling remains unchanged, exemplary results are considered in the following for the geometry with K x 1 = 8 . The pressure profiles of the UI simulations are visualized in Fig. 6 for the rigid (a)   and elastic (b) case. The contour of the regions where the hydrodynamic pressure p hd reaches the cavitation pressure p cav is visualized by orange lines. The elastic deformation drastically reduces the resulting pressure profile in comparison to the rigid simulation. The results of the same cases but with the QUICK scheme are shown in Fig. 7 for the rigid (a) and elastic (b) case for comparison. Differences between the UI and the QUICK scheme are noticeable.

Convergent Slider with Rectangular Pocket
The next test case is a convergent slider with a single rectangular pocket that introduces discontinuities in the gap height. For this set-up, a mass-conserving cavitation model is essential to predict the full-film reformulation properly. The aim of the simulations is to show the effect of the spacial discretization order on the pressure distribution when gap height discontinuities are present. For this steady case, the UI scheme eventually results in first and the QUICK scheme in second-order accuracy. The investigation is firstly done for a rigid one-dimensional geometry because it has the analytical solution of Fowell et al. [25] for comparison. Next, the two-dimensional set-up of Bertocchi et al. [9] is used on the one hand to demonstrate that the algorithm of Bertocchi et al. and the EHL-FBNS algorithm give consistent results in the rigid case and on the other hand to show that the additional consideration of elastic deformations even at traditional hydrodynamic operating conditions is of great relevance. A sketch of the one-dimensional configuration is depicted in Fig. 8a while its extension to the two-dimensional geometry is described in Fig. 8b. The analytical solution of a rigid one-dimensional converging slider with rectangular pocket and incompressible isoviscous liquid phase was derived by Fowell et al. [25] and was also used for code verification by Giacopini et al. [8]. The parameters used in the current study are summarized in Table 3. The one-dimensional geometry was replicated    [9]. Simulations were performed for the small ( L x 2 = 10 mm , w = 7 mm ) and large ( L x 2 = 30 ⋅ 10 mm , w = 30 ⋅ 7 mm ) geometry by a pseudo one-dimensional grid with three discretization points in x 2 -direction and Neumann boundary conditions for pressure and cavity fraction at the south and north boundary.
To investigate the grid convergence properties of the spacial schemes for this kind of set-up, the resulting pressure profiles of the UI and QUICK schemes are shown in Figs. 9 and 10 for different grid resolutions alongside the analytical solution. While both schemes converge toward the same analytical solution, the first-order UI scheme converges faster at much lower resolutions than the second-order QUICK scheme. It is therefore concluded that for rigid geometries, lower order discretization schemes are preferable when gap height discontinuities are present.
Next, the two-dimensional set-up of a converging slider with rectangular pocket is considered. It is the same set-up that was used by Bertocchi et al. [9] to show the agreement of their code with the formulation of Ausas et al. [4] and Giacopini et al. [8]. This set-up is simulated to demonstrate the good agreement of the EHL-FBNS algorithm with the aforementioned algorithms in the rigid case and to point out that the resulting pressure distribution is greatly affected when common elastic bodies instead of rigid ones are used.
The parameters employed in the EHL-FBNS simulations are summarized in Table 4. In the rigid case, Bertocchi et al.  Figure 11 shows the pressure distribution of the rigid UI EHL-FBNS simulations along the center line for both geometry widths next to the results of Bertocchi et al. [9]. The curve of Bertocchi et al. was read in with the software Engauge Digitizer and can therefore be subject to minor deviation from the original data. Still, the results of the EHL-FBNS algorithm and Bertocchi et al. are similar. Furthermore, Bertocchi et al. found their results being in close agreement to Ausas et al. [4] and Giacopini et al. [8]. This indicates that the EHL-FBNS algorithm is consistent with all the mentioned works.
In the following, only the small geometry ( L x 2 = 10 mm , w = 7 mm ) is considered to show that even at traditional hydrodynamic operating conditions with low pressures between 1 and 10 MPa , the employment of common elastic parameters can induce small variations in the gap height which in turn severely alter the resulting hydrodynamic pressure and cavity fraction profiles. The resulting pressure profiles of the rigid and elastic simulations are shown in Fig. 12 for the UI scheme. The boundary of the cavitation region is indicated by an orange line and the center line is marked in green. The resulting pressure profiles are severely weakened due to the introduced elastic deformation and the cavitation region has a larger downstream extension.
The gap height h, hydrodynamic pressure p hd and cavity fraction along the center line are compared in Fig. 13 for the rigid, elastic, UI, and Quick simulations. It becomes Fig. 12 Distribution of hydrodynamic pressure p hd in the inclined slider with pocket for UI discretization and rigid a and elastic b model visible that the differences between rigid and elastically deformed gap height are small compared to the corresponding differences in the hydrodynamic pressure. The maximum hydrodynamic pressure is not only reduced to about 20% , but also the general shape and peak change drastically. In the elastic case, the cavitation region extends over the whole length of the pocket, thus reducing the area of positive pressure gradient in front of the downstream end of the pocket and eventually causing only a small pressure increase compared to the rigid case. Comparing the UI and QUICK schemes shows a tendency of the higher order scheme to cause oscillations at the downstream end of the cavitation region and minor discrepancies in the pressures profiles at the end of the pocket.
From these results, it is concluded that even at low pressures between 1 and 10 MPa , small elastic deformations can significantly alter the hydrodynamic pressure and cavity fraction profiles. Consequently, results of traditional hydrodynamic simulations with mass-conserving cavitation but no elastic deformation model must be handled with care when making statements about the possible increase of the load carrying capacity due to surface textures. Furthermore, the performed simulations show that the EHL-FBNS algorithm is a useful tool to build up on the considered reference investigations since it can take the additionally required elastic deformation effects into account. Moreover, the implementation of the solver allows to switch conveniently between first-and second-order discretization schemes of the Couette term. This allows to quickly identify the preferable scheme through grid convergence studies.

Ball-on-Disc Tribometer with a Single Texture
The second set-up is the ball-on-disc tribometer with a single texture used in the simulations and experiments of Mourier et al. [26]. These were performed to investigate the unsteady effect of isolated dimples passing through an EHL contact at pure rolling and rolling-sliding conditions. Mourier et al. state that they used shallower dimples in the simulations than in the experiments because the deep dimples compromised convergence. In the following, both geometries are simulated with the EHL-FBNS algorithm to show that the presented algorithm can provide converged results in either case. Firstly, grid convergence studies are performed to identify the preferable discretization scheme. Then, the gap height measurements of Mourier et al. are used to validate the EHL-FBNS algorithm for ball-on-disc tribometers. Lastly, differences in the simulated results of Mourier et al. and the EHL-FBNS algorithm are discussed and additional deductions about the discretization schemes are drawn.
The rolling-sliding condition of the EHL contact is characterized through the slide-to-roll ratio SSR as provided by Mourier et al. [26]: The time dependent variation of the gap height due to the rigid geometry of the surfaces is expressed as [26]: where d is the dimple depth, r is the dimple radius and D is the distance of any position to the moving dimple center [26]: The dimple center has the coordinates x 1,d and x 2,d with [26]: The initial position of the dimple center is described by x 1,0 and x 2,0 at time t = 0 s . The dimple moves with velocity u tex = U low . The parameters used in the EHL-FBNS simulations are summarized in Table 5. All simulations are performed with the same mean velocity u m , but different SSR and therefore different U up , U low and number of time steps N t . In the case of SSR = 0 , N t = 257 time steps and a dimple radius of r = 15.5 μm are used. To replicate the experiment, a dimple depth of d = 7 μm is employed, while d = 0.175 μm is used to be consistent with the numerical set-up of Mourier (33) et al. [26]. For SSR = −0.5 , the dimple radius r = 21.5 μm is considered. Since the dimple moves more slowly in this case,    the FDM multigrid solver of Venner and Lubrecht [12]. As stated by Mourier et al., the accuracy of this solver is of second order in space and time. Due to the Euler implicit discretization of the unsteady term in the EHL-FBNS solver, the achieved accuracy is generally limited to first order for both the UI and the QUICK scheme.

Param. Value
In order to perform a grid convergence study, simulations with the UI and QUICK scheme at SSR = 0 , respectively, for the deep and shallow dimple geometries, are employed to assess both schemes with and without gap height discontinuities. For this grid convergence study, the less strict definition r n EHL-FBNS = mean abs ⃗ n p * is used instead of Equation (28). This is necessary in order to perform the QUICK simulations at low resolutions, otherwise a stall in the convergence of the residual r n max,G can occur at some time steps for the used value of the underrelaxation coefficient. Nonetheless, for the highest grid resolution, the stricter residual of Equation (28) does not lead to a stall and the maximum deviation in the resulting gap heights along the center line when the dimple is at position x 1,d = 0 is less than 0.1% when compared for both residual definitions. Note also that for any other simulation in this publication, the stricter definition given by Equation (28) (and thus following the suggestion by Woloszynski et al. [10]) is employed. The effect of the residual definition on the amount of required iterations for the highest resolution is discussed at the end of this section.
For the deep dimple at position x 1,d = 0 , the resulting gap height profiles along the center line obtained for different resolutions are displayed in Fig. 14 for the UI and in Fig. 15 for the QUICK scheme. The comparison of both figures shows that the QUICK scheme has a worse convergence behavior around the dimple than the UI scheme and even converges toward a different solution at the downstream end of the dimple. It is therefore concluded that in this case, the UI scheme with the first-order discretization of the Couette term is preferable over the QUICK scheme.
For the shallow dimple at position x 1,d = 0 , the resulting gap height profiles along the center line obtained for different resolutions are displayed in Fig. 16 for the UI and in Fig. 17 for the QUICK scheme. For this geometry with no discontinuity, both schemes converge toward the same solution, but the second-order discretization of the Couette term with the QUICK scheme delivers the better convergence properties and is therefore preferable.
The following simulations are performed at the highest resolution level and with the residual definition of Equation (28) to further support the previously identified arguments about the preferable discretization scheme of the Couette term. The results of the EHL-FBNS simulations are provided as videos in the supplements. Exemplary results of gap height h, hydrodynamic pressure p hd and cavity fraction at t = 0 s , SSR = 0 and d = 7 μm are shown for UI and QUICK simulations in Fig. 18. The contour line of p hd = p cav is marked in orange while the center line is displayed in green. Apart from more pronounced spikes in the hydrodynamic pressure at the downstream end of the EHL-contact zone with the more accurate QUICK scheme, both simulations produce almost the same results at first glance.
In the following, the EHL-FBNS results are compared to the experimental and simulated counterparts of Mourier et al. [26]. The data of Mourier et al. was read in with the software Engauge Digitizer (https:// marku mmitc hell. github. io/ engau ge-digit izer/) and can therefore be subject to minor deviation from the original data.   Most of the domain is still unaffected by the dimple and basically corresponds to the steady solution without a dimple. Small differences in the gap height h between UI and QUICK can be observed while the more accurate QUICK scheme closely fits the experimental results in the center of the domain. The systematic difference in the gap height between UI and QUICK is mostly due to a different value of the rigid body displacement h d which is set such that the load balance equation is eventually met. For the dimple positions 2-4, the QUICK scheme produces large deviations in the gap height compared to the experimental results at the Fig. 18 Distribution of gap height h a,b, hydrodynamic pressure p hd c,d and cavity fraction e,f in the ball-on-disc tribometer with single texture at t = 0 s for UI (left) and QUICK (right) discretization discontinuity at the downstream rim of the dimple. There, the UI scheme manages to fit the experimentally measured gap height better. These results are consistent with the finding of LeVeque [34,Ch. 11] that for discontinuous problems, first-order methods give smoother results while second-order methods cause oscillations. Both schemes deviate from the experimental results at the downstream rim of the dimple when it is leaving the EHL contact at position 5. At all five Fig. 19 At SSR = 0 (a) and SSR = −0.5 (b): distribution of gap height h (top) and hydrodynamic pressure p hd (bottom) along the center line of the ball-on-disc tribometer with single texture for UI and QUICK discretization against the experimental results of Mourier et al. [26] positions, the hydrodynamic pressure distributions of UI and QUICK are mostly in close agreement, but the UI scheme produces higher pressure spikes at the downstream rim of the dimple which in turn cause larger elastic deformations.
This explains why the UI scheme shows better agreement in the gap height since this more diffusive scheme eventually tends to smooth the discontinuity at the rim of the dimple by adjusting the hydrodynamic pressure accordingly.  [26] Consequently, the lower order UI scheme is recommendable close to discontinuities in the gap height while the QUICK scheme is more advantageous at smoother parts of the geometry.
The UI and QUICK results along with experimental values of Mourier et al. in case of the deep dimple with SSR = −0.5 are depicted in Fig. 19b. In this case, the downstream area of the dimple also gets deformed because the dimple moves at a lower speed than the mean velocity u m . This means that some of the fluid that is initially within the dimple leaves the texture behind which causes a deformation since more volume is occupied outside of the dimple. This behavior is principally also replicated by the EHL-FBNS results but in a more pronounced way than in the experiments. The higher order QUICK scheme matches the experimental results closer than the UI scheme in the vicinity of this effect as depicted at dimple position 3. The UI scheme causes higher hydrodynamic pressures downstream of the dimple which this time cause a larger overestimation of the occurring deformation than done by the QUICK scheme. However, the stronger this effects becomes, the larger the deviation between experiment and simulation becomes as shown at dimple position 4. Still, experimental and EHL-FBNS results generally show a good agreement in the gap height distribution. The EHL-FBNS algorithm is thereby validated for simulations of discontinuous textures in ballon-disc tribometers under EHL operating conditions.
Next, the EHL-FBNS results are compared to the simulated results of Mourier et al. [26].  UI results slightly differ at some points. The reason for the closer agreement of the EHL-FBNS results with Mourier et al. is expected to be due the dimple moving at a lower speed, thus introducing slower unsteady effects. Since the discrete time steps are of the same length as in the case of SSR = 0 , the time resolution is relatively higher for SSR = −0.5 , thus enabling all schemes to deliver almost the same results.
In order to better understand the performance of the EHL-FBNS algorithm in unsteady EHL conditions, the required amount of iterations N n at each position of the dimple center x 1,d is displayed in Fig. 21. The most iterations are needed to compute the steady solution at the initial position of the dimple and once the dimple reaches the EHL-contact zone. Since it produces more irregular results, the deep dimple generally requires more iterations than the shallow one. UI and QUICK simulations need a similar amount of iterations for each respective pair of simulation. Each of the eight simulations required between 2 and 6 h code execution time on a desktop computer.
The amount of required iterations can be significantly reduced if the residual definition of Equation (28) is changed to r n EHL-FBNS = mean abs ⃗ n p * . While the difference between the resulting gap height profiles is negligible as demonstrated in the beginning of this section, the less strict residual definition reduces the required amount of iterations noticeably, as displayed in Fig. 22. Furthermore, this reduces the code execution time per simulation to a range between 0.5 and 2.5 h. Summarizing, the EHL-FBNS algorithm manages to deliver converged results even in unsteady EHL operating conditions with deep dimples with strong discontinuities at their rim. Moreover, the first-order UI scheme gives closer agreement to the experimental results of Mourier et al. [26] in the vicinity of deep dimples than the higher order QUICK scheme. Therefore, lower order spatial discretization schemes are recommended close to strong gap height discontinuities while higher order schemes are recommended at smoother parts of the geometry due to their higher accuracy in these areas. Moreover, the extreme pressure spikes at the rim of the deep dimples raise the question whether the elastic half-space assumption is still valid in this area. Nonetheless, the good agreement of the gap height distribution with the experimental data of Mourier et al. validates the EHL-FBNS algorithm for simulations of textures in ball-ondisc tribometers under unsteady EHL operating conditions.

Conclusion
The EHL-FBNS algorithm is presented in this paper. It allows the simulation of the unsteady hydrodynamic pressure build up under consideration of mass-conserving cavitation with the JFO model within lubrication gaps. Furthermore, its versatile implementation enables the simulation of various combinations of isoviscous or piezoviscous flows, incompressible or compressible liquid phases, rigid or elastic surfaces, first-or second-order spatial discretizations of the Couette term, imposed rigid body displacements or normal forces and Dirichlet or Neumann boundary conditions. The algorithm is explained in detail within this paper and the implemented MATLAB © code with the corresponding setup and visualization scripts is provided in the supplements. That way, the results can be reproduced by downloading the code and repeating the simulations. Moreover, the usage and further development of the thoroughly commented code is encouraged through its public accessibility and maintenance on GitHub: https:// github. com/ ErikH ansen Git/ EHL.
Within this paper, the EHL-FBNS algorithm results were compared to analytical, simulated and experimental literature data of Woloszynski et al. [10], Fowell et al. [25], Bertocchi et al. [9] and Mourier et al. [26]. The key findings are: • The performance of the EHL-FBNS code almost scales with N log(N) in simulations with N discretization cells and a constant rigid body displacement h d . • The EHL-FBNS code can deliver converged results even when extreme gap height discontinuities are present. • Higher order spatial discretizations of the Couette term can cause large errors in the gap height distributions when gap height discontinuities are present in the EHL contact. Therefore, lower order spatial discretization schemes are recommended close to strong gap height discontinuities while higher order schemes are recommended at smoother parts of the geometry due to their higher accuracy in these regions. • In order to evaluate whether surface textures in a geometry introduce discontinuities in the gap height that in turn worsen the result quality of higher order discretizations of the Couette term, performing a grid convergence study with different discretization orders is recommended. • The EHL-FBNS algorithm is validated for the investigation of deep dimples with discontinuous rims in ball-ondisc tribometers under EHL operating conditions. • Even at traditional purely hydrodynamic operating conditions with low pressures between 1 and 10 MPa , the employment of elastic models is recommended because the resulting pressure profiles are strongly influenced by slight elastic deformations.
• The EHL-FBNS algorithm is a convenient tool to perform grid convergence studies with different discretization orders of the Couette term while taking both elastic deformation and mass-conserving cavitation effects into account.
Since it was shown that both higher and lower order discretization schemes are beneficial in different parts of a geometry, a hybrid approach is suggested for future implementation (for example by means of flux limiters [29,Ch. 4.4.6]) to achieve maximum accuracy for profiles that show both smooth and discontinuous features.

Discretized Dimensionless Reynolds Equation Considering Mass-Conserving Cavitation and Discretized Dimensionless Fischer-Burmeister Equation
Applying the FVM with a two-dimensional grid as schematically depicted in Fig. 23 to Equation (9) eventually results in the discrete dimensionless Reynolds equation considering mass-conserving cavitation at each cell center G C : where the two-dimensional cell is of area A * , has the boundary A * and its outward pointing normal vector ⃗ n is of length 1. Since a rectangular mesh aligned with a Cartesian coordinate system is used, the components n 1 and n 2 can take the values of ±1 or 0. The line increment of the cell boundary is denoted by L * . The scalar product is implied by ⋅ . The integrals are discretized with the second-order midpoint rule. Values or derivatives at the boundaries of the Poiseuille term are discretized with a second-order central interpolation or differential scheme. Its coefficients A Po read: The dimensionless cell spacings are defined as Δx * 1 = Δx 1 ∕x 1,ref and Δx * 2 = Δx 2 ∕x 2,ref . At the cell boundaries, * Po is determined as:  Special attention has to be paid close to the boundaries if no WestWest neighbor cell exists. In this case, the first-order upwind interpolation scheme is used for the approximation of w = W . Then, the coefficients of B Co read: At the cell boundaries, * Co is determined as: except when there is no WestWest neighbor cell, such that * Co,e stays the way it is but * Co,w becomes: The unsteady term is discretized with the first-order Euler implicit scheme. The coefficient of B Ti reads:    where * ,prev Ti,C and prev C correspond to the values of the previous time step. The discrete non-dimensional Fischer-Burmeister Equation (12) at each cell center reads:

Discretized Elastic Deformation
In its non-dimensional form, the gap height equation (5) Fig. 24 Schematic sketch of the alignment of the mirrored dimensionless kernel function K * with the dimensionless pressure field p * hd to obtain the dimensionless elastic deformation h * el,C via convolution This convolution can be interpreted as the alignment of the non-dimensional hydrodynamic pressure field p * hd below the mirrored non-dimensional kernel function K * as shown in Fig. 24. This is explained in detail in the following. The non-dimensional kernel function is mirrored in x * 1 -and x * 2 -direction. Afterward, the center entry K * (0, 0) is aligned with p * hd,C = p * hd (x * 1,C , x * 2,C ) . That way, aligned pairs of K * and p * hd are created. The product of each respective pair is computed, for example K * (0, 0) ⋅ p * hd,C or K * (−Δx * 1 , 0) ⋅ p * hd,E . The sum over all of the products is equal to h * el,C . If instead of h * el,C = h * el (x * 1,C , x * 2,C ) the dimensionless elastic deformation at the West cell h * el,W = h * el (x * 1,W , x * 2,W ) is desired, the center entry K * (0, 0) of the mirrored non-dimensional kernel function has to be aligned with p * hd,W = p * hd (x * 1,W , x * 2,W ). Assuming the hydrodynamic pressure being constant over the rectangular discretization cell of area A = Δx 1 Δx 2 , the kernel K can be expressed as [30,Ch. 3.3], [31,Ch. 3.1]: where the certain terms are consolidated as:

Discretized Pressure Jacobian of the Dimensionless Reynolds Equation
In the rigid case, the following procedure is obsolete and the pressure Jacobian of the dimensionless Reynolds Equation is simply J G,p * = A Po . In order to consider the relationship between h * and p * in J G,p * in the elastic case, the coefficients of the Couette and unsteady term are reformulated: with Using the previously mentioned generic order upwind interpolation scheme for the Couette term, this results in: Except if no WestWest neighbor cell exists and * Co,w has to be replaced by: Afterward, these expressions and Equations (59) and (60) are inserted into c C of Equation (34) (thus equations (56) and (57)) and only the coefficients of p * S , p * W , p * C , p * E and p * N are considered to obtain J G,p * = A Po + A h = A Po + A Co + A Ti . Eventually, a generic scheme can be found to find the coefficients of an arbitrary diagonal D of matrix A Co :