An Inlet Computation Zone Optimization for EHL Line Contacts

Most EHL numerical calculation methods considering both starved and flooded conditions, employ a fixed multiple of the Hertzian radius for the normalization of the computational domain. These methods are often used to investigate the influence of the lubricant supply on friction etc., but the solutions obtained might be numerically starved. The present numerical calculation method utilizes an optimized normalization of the computational domain to ensure that the solutions obtained are not numerically starved. With this normalization method, the computational domain can be appropriately meshed, regardless of the variability in the inlet length due to changes in the operating conditions. This method can, therefore, be used to obtain accurate EHL film thickness and pressure data over a wide range of operating conditions.


Introduction
Many moving machine assemblies involve rolling-sliding parts, such as cam mechanisms, rolling bearings, gears, etc. and these typically operate under elastohydrodynamic lubrication (EHL) conditions. There are numerous examples of investigations of lubricant film formation and friction in the interfaces between the rolling-sliding parts.
Most studies have been focused on high-load, low-speed conditions, but actual components also encounter lowload, high-speed conditions. Since severe wear and flaking may occur on the rolling-sliding surfaces due to insufficient oil film thickness, controlling the oil film thickness distribution is necessary for reliable performance over a wide range of operating conditions. Johnson [1] divided the lubrication regimes in EHL into four parts. These regimes are the piezo viscous rigid (PR), the isoviscous rigid (IR), the piezo viscous elastic (PE), and the isoviscous elastic (IE). The analysis method to calculate oil film thickness for IR was described by Martin [2], PR by Grubin [3], PE by Dowson and Higginson [4], and IE by Herrebrugh [5]. Moes [6,7] suggested a minimum film thickness formula that could be used over a wide range of operating conditions, which is an approximation formula including the transition regions PE, PR, IE, and IR. This approximation formula was curve-fitted with asymptotic lines for obtaining the minimum oil film thickness in PE, PR, IE, and IR regimes. For the transition regions between PE and PR, IE, and IR, the numerical calculation of the oil film thickness, does in many cases, experience convergence issues, which leads to inaccurate predictions.
Regarding numerical computation methods for EHL line contacts, Lubrecht [8], Venner [9], and Venner et al. [10] introduced a multigrid method based on a Gauss-Seidel solver for the Reynolds equation, applicable to the isothermal line contact problem. Venner [9] pointed out that "in some lowly loaded situations a larger inlet region was needed in order to avoid numerically starved lubrication". Furthermore, the numerical computation method for the transition regions between IE and IR was not described in detail. Alakhramsing et al. [11], Habchi et al. [12,13], and Hultqvist [14], investigated the effect of oil film thickness on friction based on a full-system finite element approach. They focused on high-load, low-speed conditions, under fully flooded lubrication. To this end, they applied the Hertzian radius as the normalization for the IE and PE regions, and the reduced radius of curvature for the IR and PR regions.
As Venner [8] concluded, "numerically starved lubrication" under lower load in the case of an inappropriate inlet distance, causes the divergence of the numerical results. This issue is a primary factor in the transition regions between the PE, IE regions, and the PR, IR regions, that causes instabilities in the numerical solution procedure.
Regarding the starved-and fully flooded conditions in EHL, Wedeven et al. [15] experimentally investigated the occurrence of starvation. Later, Hamrock et al. [16] suggested a simple expression of a computational inlet distance required to obtain a fully flooded simulation for EHL point contacts operating in the PE region. By means of this measure of the inlet distance, they also presented an approximative expression to estimate the oil film thickness for both starved and fully flooded EHL contacts operating within the PE region. For continued work devoted to EHL and starvation effects in the PE region, see [17][18][19][20][21]. Examples of work considering starved and fully flooded conditions in the IR region are the ones by Anandan [22,23] and Biboulet [24], where the reduced radius of curvature was employed as a normalization of the computational domain that they used to numerically calculate the oil film thickness.
In this paper, an expression of the inlet distance that ensures fully-flooded conditions is employed to normalize the computational domain. The numerical EHL line contact calculation model, for isothermal conditions, is implemented in a fully-coupled finite element-based approach. By means of the present normalization method, the corresponding computational domain can be appropriately meshed, regardless of the variability of the inlet length, to accommodate changes in the operating conditions. This method can, therefore, be used to obtain accurate EHL film thickness and pressure data over a wide range of operating conditions, including the transition regions between PE, PR, IE, and IR.

The EHL Line Contact Model
This section presents the numerical EHL line contact model that incorporates an optimized inlet computational domain, which provides solutions that are not numerically starved.

Inlet Computation Domain Optimization
Let us first revisit the experimental research results by Wedeven et al. [15], who investigated starvation in ball bearings, by means of optical analysis. Figure 1 shows a cross-section of the contact geometry with particular emphasis on the inlet zone of the contact. The length of the inlet zone is crucial for the numerical computation of oil film thickness and pressure, since a too short inlet computational domain may lead to numerically starved conditions. This is, for instance, a big issue when simulating the conditions that are typical for the iso-viscous rigid (IR) regime. The numerically starved condition induces a reduction of oil pressure and converged solutions cannot be obtained.
Wedeven et al. [15] presented an expression for the inlet zone length S f required to obtain fully flooded conditions, i.e.
where is the Hertzian radius, R � −1 = R −1 1 + R −1 2 the reduced radius of curvature, h 0 the central film thickness, and where E ′ is the reduced elastic modulus, and E i and v i are Young's modulus and the Poisson ratio of bodies i = 1 and 2 , respectively.
They also discovered that the onset of starvation occurs if the ratio of h b ∕h 0 of the oil film thickness at the inlet boundary h b = h x b to the central oil film thickness h 0 = h(0) falls below 9.
As previously mentioned, most EHL numerical calculation methods use a multiple of the Hertzian radius a to define the length of the computational domain. This may, however, lead to numerical starvation if the resulting inlet zone does not encompass the inlet distance S f , as defined in [15]. This issue becomes a complication when performing calculations over a large range of operating conditions, which is often required when simulating an actual application, as it leads to less accurate oil film thickness and pressure estimates. If we can correctly adjust the inlet distance S f for fully flooded conditions independent of numerical calculation engineers' technique, e.g. adjustment of required mesh size along with the inlet distance, the numerical calculation method can be used universally.

Governing Equations
In this section, the theoretical foundation of the fully-coupled finite element approach, in which the Reynolds equation and the equations governing linear elasticity are solved simultaneously, is presented. For line contacts, the Reynolds equation governing stationary flow is given by: where is the oil density, 0 the oil viscosity, h the oil film thickness, p the oil film pressure, U e = U 1 + U 2 ∕2 the mean entrainment speed and where is a penalty term, see also [14]. The computational domain is defined by where x b is the location where h = h b and x o is specified to ensure that it is larger than the point where the film ruptures and cavitation occur. The fluid-structure interaction in the EHL contact is manifested by the appearance of the total vertical displacement v Ω of the boundaries of the deforming solids' Ω , in the oil film thickness equation, i.e.
where h 00 is a measure of rigid body separation. This parameter is determined via force balance in the y-direction, i.e.
where w is the applied load per unit width. For a schematic illustration of the EHL line contact using present nomenclature, see Fig. 2. The system of equations governing the mechanics of the linear elastic body is where u and v are the displacements in the x -and y-direction, respectively, and c i , i = 1, 2, 3 , are the equivalent components of the stiffness matrix given by where E eq and eq are the equivalent Young's modulus and the equivalent Poisson's ratio of the contacting bodies, respectively. These are defined as in [25], i.e., As an approximate equation of the pressure-viscosity relationship, Roelands' equation is used. This equation can be expressed as, where p 0 = 1.9608 × 10 8 Pa and is the Pressure-viscosity coefficient. To describe the density's variation with pressure, Dowson and Higginson's relationship is adopted here. That is, where p r = 5.9 × 10 8 Pa and 0 is lubricant's density at zero pressure. This value for p r and the value 1.34 for the asymptote at infinite pressure are used to enable comparison with frequently available results, but it should be noticed that a much closer fit may be obtained by fitting to experimental data for a specific lubricant, see e.g. [26]. At high load conditions, the oil pressure distribution has fluctuations. This instability can be avoided with the Galerkin least-squares finite element method (GLS). For details, the reader is referred to [12].

Dimensionless Model
In this section, the dimensionless model of the elastohydrodynamic governing equations is introduced. In order to compare the results with the conventional model, the final calculation results in terms of minimum oil film thickness will be expressed in Moes's dimensionless parameters ( M, L ). However, a special non-dimensional model for numerical calculation is desired to avoid numerical starvation. A system of dimensionless equations based on b ′ should be expressed by To ensure fully-flooded conditions, the computational domain for the EHL line contact model employed in the present analysis is normalized with the inlet distance a + S f , instead of the Hertzian radius a . Moreover, since the minimum oil film thickness occurs at the center of the contact in the IR and PR regions, h 0 in Eq. (1) was replaced by Moes' approximate expression of the dimensionless minimum oil film thickness, h Moes min (M, L) , presented in [7], and thus: Based on this we define the formula for the inlet distance b ′ as As previously mentioned, the majority of the studies devoted to numerical calculations of EHL contacts are using a multiple of the Hertzian radius a to define the computational domain. The computational domain in the present approach is normalized by using an optimized definition of the inlet distance, i.e. b ′ , and the following normalization to transform the system of equations into the dimensionless form has been used: Moes min (M, L) M . where H 00 is the dimensionless rigid-body displacement and v Ω the dimensionless total displacement of the solid's upper surface in the vertical direction.
The dimensionless form of the force balance equation can be described as, and, in dimensionless form, the system of equations for the elastic displacements is given by, where c 1 , c 2 , and c 3 are defined in Sect. 2.1. Notice that for the dimensionless representation Eq. (22), E eq in Eq. (9) has been redefined and is now given by In the case of a line contact, a 2D computational domain can be used for the linear elasticity equations, by assuming plane strain condition in the z-direction. The non-dimensional size of the computational domain Ω D is taken as 60 × 60, based on the study by Habchi et al. [12], which revealed that this is sufficiently large to capture the stresses and displacements with adequate accuracy. A graphical illustration of Ω D and its boundary Ω D = Ω E ∪ Ω W ∪ Ω S ∪ Ω N is given in Fig. 3, where Ω C ⊂ Ω N .
The dimensionless computational domain Ω C for the EHL contact problem [where the Reynolds' Eq. (18) and the film thickness Eq. (20) are defined], which is part of the upper boundary of Ω D , is also depicted in Fig. 3. As illustrated in the insert to the right, the EHL contact spans the inlet region −4.5 < X < −1.5 , the central region −1.5 < X < 1.5, and the outlet region 1.5 < X < 2.5 . The same dimensionless domain was also used by Hultqvist [14], but in contrast to the present work where the inlet zone length b ′ is used to transform the z = ln 0 + 9.67, (26) (P) = P r + 1.34P P r + P . where n and t are the normal-and tangential stress, respectively.

Results and Discussion
This section describes the results produced using the numerical calculation model in Sect. 2. For verification of this numerical analysis, we confirmed the convergence of the numerical calculation model, oil pressure distribution, and oil film thickness distribution. We can describe Moes' dimensionless load parameter M and dimensionless

Convergence
In this section, we describe how the mesh convergence of the current model was investigated in terms of accuracy. Regarding the accuracy of the elastic deformation along the boundary Ω C , Habchi et al. [12] showed that could be computed with adequate accuracy even on a rather coarse mesh of the elastic body with a computational domain Ω D . This result was confirmed herein by comparing the results obtained on a normal and an extremely coarse mesh for Ω D , when the maximum mesh size of Ω C was specified as 0.00066, 0.001, and 0.00325. The normal and extremely coarse meshes when the maximum mesh size in Ω C was specified as 0.001 are depicted in Fig. 5. Table 1 shows the number of boundary elements N C and the number of domain elements N D , of the extremely coarse mesh, corresponding to the three maximum boundary element sizes 0.00066, 0.001, and 0.00325.
The corresponding three extremely coarse domain meshes are shown in Fig. 6. As the maximum element size only affects the mesh of the upper boundary along the EHL contact, Ω C , the change in the depth direction in the mesh of Ω D is marginal.
The resulting oil film pressure distributions from an investigation of the accuracy of the representation of the pressure spike for (L, M) = (10, 200), obtained with three domain meshes corresponding to the maximum element size of Ω C ; 32,864, 21,638, and 7258, are depicted in Fig. 7. The right side presents a magnification of the pressure spike region, and a convergent trend in the location of the pressure spike can be observed. The location and magnitude of the pressure spike are both almost the same with 21,638 Moreover, the mesh with 7258 domain elements seems not well enough to converge in the location of the pressure spike. Therefore, the mesh with 21,638 seems well enough to converge in the location of the pressure spike. This suggests that the mesh with 21,638 domain elements renders the best possible prediction, given the solver settings used. For this reason, this mesh is also employed when obtaining the reference solutions used in the comparison presented below. The dimensionless minimum oil film thickness h Moes min obtained using the meshes with N D = 7258 and N D = 21638 elements, as well as the relative difference between them (using the solution at the finer mesh as reference) are shown in Table 2. The average computing time for N D = 7258 and N D = 21638 elements on an intel(R) core(TM) i7-8700 CPU was 23 s and 1 m 13 s, respectively. The calculation conditions were chosen using representative values of Moes' dimensionless parameters M and L for each of the IR, PR, PE, and IE regions. The reason to focus on the minimum oil film thickness is that it is a key parameter in the expression of the inlet distance (15) which is investigated herein. As shown in Table 2, the dimensionless minimum oil film thickness decreases when the number of boundary elements N C (used to discretize the domain for the EHL line contact) decreases. It is, however, reassuring that the maximum relative error of the dimensionless minimum oil film thickness h Moes min is less than 0.28% (in the case L = 10 and M = 200 representing the PE region) with the coarser mesh having 7258 domain elements. Table 3, presents the dimensionless maximum oil pressure P max , and for the PE region the pressure spike P s , obtained using the meshes with N D = 7258 and N D = 21638 elements, as well as the relative difference between them (using the solution at the finer mesh as reference). As the table shows, the dimensionless maximum oil pressure P max and the pressure spike P s decrease along with the decrease of boundary elements number on the EHL contact line Ω C . It also shows that, the maximum relative difference of the dimensionless   , the relative difference of the dimensionless maximum oil pressure ΔP max ∕P max is less than 0.03%. These relative differences are very small, but it is worth noticing that the results presented herein were obtained from the mesh with 21,638 domain elements, to ensure an adequate representation of the oil pressure distribution.

Oil Pressure Distribution
In this section, the oil pressure distributions are shown for the computation domain to confirm the flooded conditions. As a verification, we confirmed both ends of the oil pressure distribution. If the oil pressure at the boundaries of the contact is zero, regardless of various inlet distances, the present numerical calculation method is shown to ensure the flooded condition.
The dimensionless oil pressure distribution on the computation domain is shown in Figs. 8,9,10,11. Each of these figures shows the 20 numerical computations for M  (enumerated from 1 to 20) starting from 200 and going down to 0.1, for L = 1, 2.5, 5, and 10 , respectively. It means that these results represent the oil pressure distributions spanning from the isoviscous elastic (IE) to the isoviscous rigid (IR) lubrication region. As shown in Figs. 8,9,10,11, the derivative of the oil pressure distribution is 0 at both the inlet and exit, verifying that the inlet region in the present numerical calculation method effectively prevents starved conditions. It can also be seen that the oil pressure distribution expands towards the inlet boundary as the value of M reduces.

Oil Film Thickness
In this section, the Moes' approximate expression h Moes min (M, L) for the minimum oil film thickness is compared with the minimum oil film thickness obtained using the present model, and the comparison is conducted over a large portion of the Johnson chart [1].
The calculation results of the dimensionless oil film thickness h Moes , for spanning the same range of the dimensionless load parameter M and material parameter L are depicted in Fig. 12. As shown here, the dimensionless minimum oil film thickness increases along with a decrease in M . Moes' approximate expression results [7] are in good agreement with the dimensionless minimum oil film thickness trend of the present model over a large portion of the Johnson chart [1]. A comparison of results of h Moes min (M, L) calculated by Moes' approximate expression vis a vis the present model is presented in Table 4, which shows that the relative difference between the two values is less than 2.7%

Conclusions
In this paper, we have presented a numerical EHL calculation method for oil film thickness that ensures fully flooded conditions by employing the present inlet distance normalization. The results obtained with the present model were compared with results from Moes' approximate minimum oil film thickness expression.
The main conclusions are: • With the definition of the inlet distance established herein, we confirmed that the oil pressure distributions obtained over the whole Johnson chart correspond to fully flooded conditions. • The maximum relative difference found in the mesh convergence study occurs in the PE region [represented here by (L, M) = (10, 200) ] and a mesh with 21,638 domain elements was employed to get a sufficiently accurate representation of the pressure spike.
-The maximum relative difference between dimensionless minimum oil film thickness obtained with the coarsest mesh (with 7258 domain elements) and the reference solution (obtained with the mesh with 21,638 domain elements) is less than 0.28%. -The maximum relative difference of the dimensionless maximum oil pressure obtained with the coarsest mesh (with 7258 domain elements) and the reference solution (obtained with the mesh with 21,638 domain elements) is less than 0.5%.
Regarding the oil film thickness, the relative difference between the prediction of the minimum oil film thickness using the present model and that by Moes' approximate expression is less than 3% over the whole Johnson chart.