An assessment of quantitative predictions of deterministic mixed lubrication solvers

The ability to simulate mixed lubrication problems has greatly improved, especially in concentrated lubricated contacts. A mixed lubrication simulation method was developed by utilizing the semi-system approach which has been proven to be highly useful for improving stability and robustness of mixed lubrication simulations. Then different variants of the model were developed by varying the discretization schemes used to treat the Couette flow terms in the Reynolds equation, varying the evaluation of density derivatives and varying the contribution of terms in the coefficient matrix. The resulting pressure distribution, film thickness distribution, lambda ratio, contact ratio, and the computation time were compared and found to be strongly influenced by the choice of solution scheme. This indicates that the output from mixed lubrication solvers can be readily used for qualitative and parametric studies, but care should be taken when making quantitative predictions.


Introduction
Contacting surfaces in mechanical components are important for transmitting power and motion. To reduce the power loss and improve the efficiency and life of mechanical components, the contacting surfaces are usually lubricated [1]. In engineering practice, rough surfaces are universal. and they affect the thickness of lubricant film within the contact. When the film thickness is not large enough to separate the two surfaces, solid to solid contact starts to occur and ultimately, the effect of roughness becomes dominant.
Under such conditions, the lubrication state is referred to as the mixed lubrication regime where fluid lubrication and solid to solid contact occur simultaneously. Mixed lubrication is inevitable in engineering applications. Therefore, modelling and predicting the performance of the mixed lubricated contacts is one of the key challenges and popular topics in lubrication science, and great efforts have been made in this field.
Generally, there are two types of mixed lubrication models: stochastic (roughness defined by stochastic parameters) and deterministic (roughness defined at each grid node).
One milestone in the stochastic mixed lubrication model development was made by Patir and Cheng [2]. Later, several improvements to the stochastic model were introduced [3][4][5].
Stochastic models use mean values, which results in loss of localized information, such as detailed pressure and film thickness variations which are important for studies on lubrication breakdown and failures. Therefore, stochastic models will not be discussed further in this study.
To solve the Reynolds equation, it is discretized to form a set of linear algebraic equations. The coefficient matrix of the discretized Reynolds equation is generally constructed from the Poiseuille flow terms alone. Deterministically solving a mixed lubrication model involves dealing with very thin lubricant films. Such thin films weaken the validity of Poiseuille flow and make the coefficient matrix lose its diagonal dominance, which results in loss of efficiency and accuracy of the linear algebraic solvers. A solution to this problem was proposed by Ai et al. [9] in the form of semi-system approach. The basic concept of the semi-system method is that both Poiseuille flow terms and Couette flow terms are used to construct the coefficient matrix, ensuring diagonal dominance even under severe contact conditions. Based on the semi-system method, Hu and Zhu [13] proposed a mixed lubrication model for point contact elastohydrodynamic lubrication (EHL) problems which could handle the complete lubrication behavior from boundary to full film lubrication regimes. According to this model, when the film thickness was small enough, localized contact occurred. Parametric [16] and experimental [17] studies have proved the capability of this model.
The Reynolds equation is solved by first discretizing (several different discretization schemes available) and then solving the resulting linear algebra problem. Until now, it is presumed that all the different solution methods and discretization schemes produce identical results under similar working conditions. In the past, Liu et al. [14] studied the effects of differential scheme and mesh density on the point contact EHL film thickness predictions and recommended that for ultrathin films (below 10~20 nm), dense meshes should be used. They also suggested that the first-order backward schemes are better when dealing with these ultrathin films. Further, Zhu et al. [15] pointed out that the accurate calculation of roughness derivatives is critical for ensuring numerical accuracy. But no study addresses the issue of repeatability of predictions from a mixed lubrication (EHL) solver due to the different implementations of the mathematical model. Therefore, the current study aims at developing this understanding by utilizing a series of key implementations. Various implementation cases are selected by considering, the information included into the coefficient matrix (main diagonal or tridiagonal), the different ways of dealing with the density derivative (numerical differentiation or differentiating its empirical equation) and finally linking these to the separate and combined implementations of the Couette terms. The results from these implementations are presented as pressure distribution, film thickness distribution, lambda ratio, and the contact ratio. A comparison among the individual cases was made foror the waviness surfaces as well as the numerically generated rough surfaces. To ensure the validity of the simulations, the results for the waviness surface are first compared against the work of Venner et al. [18]. The findings from the current study have direct implications on the use of semi-system method in lubrication simulations. This paper is organized by first presenting the equations representing the mixed lubrication regime. Both dimensional and non-dimensional equation sets are given. Then the different discretization cases are outlined, and the derivation of the density term is discussed. Finally, six representative cases are chosen to analyze the consistency of the solver and the results are presented for all these cases. The predictions from the current study are compared against Venner's work [18] to establish the validity of the current model.
Finally, results are presented for the rough surface simulations and discussion is made on the use of the findings from the lubrication solvers for quantitative predictions.

Details of the mixed lubrication model
A typical set of basic equations that formulate a rough surface point contact EHL problem are given in this section [13]. The definitions of all the symbols used below are given in the Nomenclature.
The steady state Reynolds equation is used to focus our attention on the discretization schemes and not the transient effects: The x and y represent the coordinate axes, p is the pressure distribution, h is the film thickness distribution,  is the density,  is the viscosity, and u1 and u2 are the velocity of ball and disc surfaces.
In Eq. (2), the first term on the R.H.S. gives the minimum undeformed gap, h0, the second and third terms represent the paraboloid Hertzian macro-geometry where Rx, Ry are radii of curvature in x and y direction, the fourth and fifth terms give the roughness (microgeometry) of the contacting bodies ( ) while the last term is the deformation, v.
Elastic deformation equation [19]: where Ee is the equivalent Young's modulus.

Load balance equation:
( ) =, w p x y dxdy   (4) where w is the applied load.
These four equations (Eqs. (1-4)) together form the complete system of equations describing the point contact EHL problems. Due to the fact that since extremely highpressure value is observed in the contact zone, the viscosity and density of the lubricant are also function of pressures. The Barus viscosity pressure equation [20] is used to represent the changes in viscosity and the Dowson-Higginson density pressure equation [21] represents changes in density in this study and are given below.
Barus viscosity pressure equation: where  is the pressure-viscosity coefficient, and  is the ambient viscosity.
Density pressure equation: where p is the pressure in GPa, and  is the ambient density.
The detailed information about the non-dimensionalization for these equations can be found in Appendix A. The resulting non-dimensional Reynolds equation is, It is important to mention that the Poiseuille flow terms are discretized by the second order central difference method and the discretization is given in Eq.
For the Couette flow terms, the first order backward difference method is used, as suggested by Liu et al. [14]. But it should be pointed out that the Couette flow terms, Two further cases result due to the density term, X  . Due to the fact that the density variation as a function of pressure is given by the formulation, Eq. (6) which is an empirical equation, the derivative of density can be estimated by either taking the derivative of this empirical equation (6) or using numerical differentiation. The latter method is given in Eq. (11).
Instead of calculating the derivatives of density by direct finite difference, Ai et al. [9] used the chain derivation rule, given in Eq. (12).  The non-dimensional elastic deformation V is determined by the pressure distribution.
The discrete equation used to compute the elastic deformation is 11 , , , , 11 In this equation, D is called the influence matrix and is calculated based on Eq. (3), M and N are the number of discretization points in the x and y directions, respectively. The detailed procedure for obtaining the influence coefficients, D can be found in work [22].
Eq. (14) can be expressed as a function of unknown pressures, Pij, Pi+1 j, and Pi-1 j, and the deformation is also expressed in the form of multiple diagonal terms which can then be added to Eq. (13) to implement the semi-system method. Whether the influence coefficient has been written in terms of the pressure, Pij alone or in terms of the three pressures, Pij, Pi+1 j, and Pi-1 j, two different formulations result: one with main diagonal terms only and the other with main as well as secondary diagonals, respectively. Huang et al. [23] only extracted the main diagonal terms, but Ai et al. [9] also extracted the secondary diagonal terms.
Based on the above discussion, different possibilities for formulating the system are possible by combining the single and multiple diagonal formulations, the separate and combined approaches and the different density derivatives. A total of six choices result for implementing the semi-system method and are summarized in Table 1. The numbers '0' and '1' in the table represent the different combinations. '1' means that the method described in the header of the particular column was used and '0'means that the method mentioned in the header of that particular column was not used. All the detailed expressions and formulae for these different implementations are provided in Appendix B.

Contact configuration and numerical simulation details
Two sets of simulation cases were considered. The first set of simulations were performed for the waviness surface and the second set of simulation cases used computer generated rough surfaces.
The input parameters of the EHL problem and waviness surface were extracted from the work of Venner et al. [18]. The equation for generating waviness surface is, Where Rw is the waviness surface data which is used as roughness data in Eq. (2), Ax is the non-dimensional waviness amplitude, Xs is the location of waviness start in x direction,  is the non-dimensional wavelength. The waviness amplitude is 0.08 m and wavelength is 59 m [18]. In this paper, the starting location of waviness is at the start location of the calculation domain. The calculation domain is Xs ≤ X ≤ Xe and Ys ≤Y ≤Ye, where Xs = -2.5, Xe = 1.5, Ys = -2, Ye = 2. The discrete grid density is 1024×1024, unless otherwise stated. This grid density has been proved to be dense enough in literature [14,15].
The parameters of the EHL problem are listed in Table 2. The values of dimensionless speed parameter U selected in this paper varies from 4×10 -16 to 4×10 -11 where the results for U = 4×10 -11 correspond to the simulations results of Venner et al. [18] and are used to validate the implementation methods and the computer codes used in this paper.
Besides the waviness surface, numerically generated Gaussian isotropic rough surfaces were used in this work. The rough surfaces were generated based upon the method from Wu [24], and the generated rough surface data, as used in this study, is also provided as supplementary data with this paper. The attached data is a random number matrix (1024×1024) following Gaussian distribution with a unit root-mean-square variance (Rq).
Rough surfaces with any other Rq values can be obtained by multiplying the data with that specific Rq value. In the current study, the Rq used is the same as the amplitude of the waviness surface. for both pressure distribution and load balance calculation were set as 1×10 -6 . Moreover, to improve the numerical calculation efficiency, as suggested by Liu et al. [14], when the dimensionless film thickness values were less than 1×10 -6 , solid contact was assumed to occur. It should be noted that the criterion used here for establishing a solid contact still lacks firm theoretical basis. Zhu [1] discussed this issue and tried to connect this threshold value to the molecular length of lubricant. To ensure that the selection of this criteria has no effect on the outcomes presented in the current study, the value of minimum film thickness for establishing solid contact was kept the same for the different implementations of numerical models presented in this study. Another important aspect about the asperity contact is that the contact behavior is resolution dependent [25]. Thus, the input roughness data and the size of solution domain are kept fixed for comparing corresponding simulations.
With ultrathin lubricant films, it is hard to plot fine film thickness distribution in the Hertzian contact zone by common plotting methods as the different regions with different film thickness values cannot be identified due to the very small overall film thickness values involved. Venner et al. [18] used pseudo interference graphs to tackle this problem.
The dimensionless film thickness data was represented by the corresponding intensity values based on Eq. (18).
where  is non-dimensional wavelength for which the default value used in this study is 0.05.
For the simulation results with rough surface, the lambda ratio () and contact ratio (Ac) are used to identify different contact configurations. In this work, the lambda ratio is not defined as the minimum film thickness value divided by the composite surface roughness as the minimum film thickness loses generality when rough surface is used.
Instead, the average film thickness value in the Hertzian contact zone (hac) is used. The corresponding equation is The contact ratio is defined as the number of contact grid points divided by the total number of grid points in the Hertzian contact zone as shown in Eq. (20).
where Nc is the number of contact points in the Hertzian contact zone, and NHertz is the total number of grid points of the Hertzian contact zone.

Waviness surface (validation and discussion of results)
In Fig. 1, the pressure and film thickness profiles are presented at Y = 0, from the present work for all six implementation cases. All these schemes produce similar results with slight differences in the pressure peak height towards the exit of the contact region. In were used to ensure that the differences were not related to the grid. The pressure and film thickness results are also extracted from reference [18] and plotted in Fig. 2 for comparison.
It can be readily seen that the qualitative and quantitative match is very good but, once again, slight differences exist in the prediction of pressure peak height towards the exit of the contacting region.
To further investigate the effect of the numerical implementations, the contours of corresponding pressure and film thickness distributions are shown in Fig. 3 and Fig. 4, respectively. It is evident from these contour plots that the converged results from all the six implementations are the same. It is important to note that in Fig. 4, the method proposed by Venner (given in Eq. (16) [18]) is used to plot film thickness contours.
Next, the speed parameter in the simulations was reduced from U = 4×10 -11 to lower values of U = 4×10 -12 , 4×10 -13 , and 4×10 -14 successively. Fig. 5 shows the corresponding pressure and film thickness profiles at Y = 0 for these different speeds. It shows that as the speed decreases, the discrepancies between the different implementing methods increases.
In order to show this effect more comprehensively, pressure and film thickness distributions in the Hertzian contact zone with U = 4×10 -14 are plotted in Fig. 6 and Fig. 7 respectively. These plots correspond to the profiles in Fig. 5. No big differences in the pressure distribution results with different implementing methods can be observed.
However, the implementation method used has dramatic influence on the film thickness distribution.
In the following, the results from these different implementations are analyzed to access the key reasons for the discrepancies among the output. To facilitate discussion, the Therefore, for the waviness surfaces, the output of the solver is not reliable when the effect of Couette terms is only included in the main diagonal terms.
An interesting phenomenon can be observed by looking at the output (c) and (d). These two methods only differ in the definition of the density derivative (differentiating the empirical equation or using finite differences) but the output from both the cases is different. Also, it seems that the output from (c) is not correct. Therefore, it can be inferred that the solution method can be stabilized if the density derivative is defined through the chain rule (see Eq. (12)). This stability can be linked to the additional pressure derivative that appears in the derivation of density.
The next task is to compare the time involved in each of these implementation methods. Table 3  With the relative computation time required and the accuracy for each of these cases in mind, it is believed that, for the waviness surfaces, the (0,0,0) implementation that considers separate form, effect of Couette flow included in the secondary diagonals as well and the density derivative implemented through finite differencing, is the most suitable for computations and is recommended.

Numerically generated isotropic Gaussian surface
Next, a critical assessment is made on the effect of these implementations on the output of mixed lubrication solvers when numerically generated rough surfaces are used. The 3dimensional visualization and roughness profiles along X = 0 and Y = 0 of the generated surface is shown in Fig. 8. In the first set of simulations, the speed parameter is fixed as U = 4×10 -13 . This speed value is selected to ensure that it is a full film case for all the implementations. The simulations with lower speed where asperity contact is impossible to avoid, will be presented later in this section.
The pressure and film thickness profiles at Y = 0 are presented in Fig. 9 for all the different implementation methods. Contrary to the waviness surface case, there are obvious differences among the different profiles corresponding to the each of the methods used.
The contours of pressure and film thickness within the Hertzian contact zone are further shown in Fig. 10 and Fig. 11. All these distributions were normalized to the interval between 0 and 1 to generate these graphs with the same color scale. The corresponding lambda ratio and relative computation time are also compared and are listed in Table 4.
The actual computation time of using (0,0,0) method is also provided.
The film thickness distribution plots in Fig. 11 are much clearer than the pressure distribution plots in Fig. 10. As the differences can be seen more clearly in Fig. 11, only the differences in the film thickness predictions are discussed. Once again, to facilitate discussion, the different implementations in Fig. 11 have been labelled as (a) to (f). First, it is important to note that the output of the solvers for the rough surface simulations is quite different from the output by the respective solvers for the case of waviness surfaces. To further provide evidence to support our argument, the lambda ratio values are listed in Table 4. These values also show that the output is solver or implementation dependent.
Corresponding to the film thickness distribution plots in Fig. 11 The trend is similar to that observed in the case of waviness surfaces at low speeds. This indicates that for the generated rough surfaces without asperity contact, only including the effect of Couette flow in the primary diagonal could reduce the computation efficiency.
In order to test the different implementing methods with asperity contact, an ultralow speed parameter value of U = 4×10 -16 is imposed, and the simulations are repeated. Fig. 12 shows the pressure and film thickness profiles at Y = 0 for all the different implementations.
The contours of contact area are plotted in Fig. 13 instead of contours for pressure and film thickness. This is to ensure that the effect of asperity contact is actually compared. Table 5 lists the lambda ratio, contact ratio, and relative computation time values with the actual computation time of using the (0,0,0) implementation. From Fig. 12, the discrepancies in pressure and film thickness profiles is less compared to the case of no asperity contact for the different implementations. The different implementations have been labelled again to facilitate discussion. Based on the contact area graphs (Fig. 13), the pair (a)  results, but the effect is less pronounced compared to the case of no asperity contact.
According to Table 5, the lambda ratio and contact ratio values could be divided into three groups as well, which is similar to the lambda ratio values in Table 4. Interestingly, as the effect of different formulations of the density derivative terms is less obvious compared to when no asperity contact occurs. This can be clearly seen when comparing case (a), (b), (c), and (d) as these result in nearly the same lambda ratio and contact ratio values. Moreover, it could be inferred that these outputs have convergent trend as the severity of contact increases. Then, the question comes to how to choose between the combined form and the separated form of the Couette term. Although some researchers preferred the separated form [9,10,13,15], there is still no solid evident that which one is the most accurate and should be recommended.
Next, the computation time in Table 5 is analyzed. Unlike the two cases presented above (waviness surface and the rough surface without asperity contact), the (0,0,1) implementation is the most efficient (shortest computation time). Interestingly, the computation time for the (0,0,0) implementation (shortest for the waviness and no asperity contact case) is even longer than the (0,1,0) and the (0,1,1) cases. Another important thing to notice is that the differences in relative time among all the different implementations is not as great as the case of waviness surface or the case of rough surface with no asperity contact.
In summary, the implementation method for the semi-system method has dramatic influence on the simulation output. A key point to stress is that the effect varies from case to case and is highly dependent upon the problem being simulated. For the waviness surface, like the rough surface, as the lubricant film gets thinner, the discrepancies become greater. However, when using numerically generated isotropic Gaussian rough surfaces, unlike the waviness surfaces, modifying the diagonal or secondary diagonal terms does not affect the output. In this case, other parameters like, the combined or separated form of the Couette terms and the chain derivation rule for the derivatives of density are the key factors affecting the output. In case of severe asperity contact, when the speed reduces, only the combined or separated derivatives of the Couette term matters and the effect of the different formulations for the density derivative is negligible.
A comparison of computational times further complicates the process of selection of the adequate implementation method as, once again, it is dependent upon the problem at hand. For the simulation cases presented above, as long as the waviness surfaces or rough surfaces (without asperity contact) are considered, the (0,0,0) implementation is computationally the most efficient but when rough surfaces with asperity contact are considered, the (0,0,1) implementation gives the least computation time.
Finally, combining the effect of different implementation methods on the accuracy and efficiency of simulation output, the (0,0,0) implementation method is recommended for waviness surfaces. However, it is hard to say which implementation method should be recommended when generated rough surfaces are used.

Conclusions
In this paper, mixed lubrication solvers were developed by considering the combined Therefore, when using EHL and mixed lubrication solvers care should be taken when using their output for quantitative predictions. And more attention should be paid to how to get enough accurate output through deterministic mixed lubrication solvers.    Table 1 Six implementation methods of the semi-system method Table 2 Parameters of the EHL problem Table 3 Relative simulation time (Waviness surface) Table 4 Lambda ratio and relative simulation time (Rough surface, U = 4×10 -13 ) The Hertzian contact formulas are used as follows: