Elastohydrodynamic lubrication of coated finite line contacts

Abstract

In this work, a finite element-based model is presented that simulates elastohydrodynamic lubrication in coated finite line contacts. Using this model, the film thickness and pressure distributions, between a straight roller with rounded edges on a plate, were analyzed. The model was successfully validated against representative results reported in literature. Parameter studies were conducted to study the influence of varying operating conditions, axial surface profile parameters and coating mechanical properties on the overall elastohydrodynamic lubrication behavior of the contact. It was found that in contrast with typical elastohydrodynamic lubrication behavior, the maximum pressure and minimum film thickness, which are located at the rear of the contact, are largely influenced by variations in load. Results also reveal that axial surface profile parameters and coating mechanical properties may act as amplifiers to the effect of load on pressure and film thickness distribution and can thus, if smartly chosen, significantly enhance lubrication performance.

Introduction

From past literature, dedicated to elastohydrodynamic lubrication (EHL), one may retrieve that a lot of work has been done on elliptical point and infinite line contact problems. In the latter, a uniform pressure distribution in axial direction is assumed. However, significantly less work has been done on finite line contact problems despite the high importance of the topic. Due to the finite lengths of components, high stresses are generated towards the extremities of the contact, often referred as edge loading. Typical examples include cam-roller followers pairs, rolling element bearings, gear teeth, etc. Axial surface profiling of components is therefore often utilized to minimize edge loading.

Depending on the type of surface profiling, the pressure and film thickness distribution may deviate significantly from that predicted using the infinitely long line contact assumption. The foundation of these observations was laid by experiments performed in 1967–1974 by Gohar and Cameron ¹ and Wymer and Cameron. ² In the latter, they measured the film thickness distribution for tapered rollers on a glass plate using optical interferometry. They were able to show that film shapes near the ends were very different from those at the central plane. Moreover, the absolute minimum film thickness always occurred near the extremities of the roller. The minimum film thickness and maximum pressure are crucial design parameters as they significantly affect wear and fatigue and thus ultimately the service life of the component.

Mostofi and Gohar ³ were one of the first to present a numerical solution to the finite line contact EHL problem. The type of rollers used was those with a straight length and rounded edges. However, the numerical results near the position were profiling starts were physically inconsistent. Using the same profiled rollers as in Mostofi and Gohar, ³ Park and Kim ⁴ obtained improved contour plots for the pressure and film thickness distribution using an improved numerical scheme as described in Park and Kim. ⁵ They also concluded that the maximum pressure and minimum film thickness always occur near the position were the profiling starts. The aforementioned numerical studies were rather limited to low or moderate loads. Extension to higher loads was made by Kushuwaha et al. ⁶

Shirzadegan et al. ⁷ recently presented a finite element-based model applicable to finite line contacts. The model developed by Shirzadegan et al. ⁷ is an extension of the pioneering work of Habchi et al. ⁸ and can easily cope with highly loaded situations by means of numerical stabilization. In Shirzadegan et al. ⁷ , different types of axial profiling were considered, i.e. rounded edges, logarithmic and crowing, and their influence on lubricant performance.

The aforementioned studies serve sufficient knowledge to perform more in-depth investigations in order to gain a fundamental understanding into the design limits of finite line contact problems. In most practical engineering applications operating conditions, such as entrainment velocity, radii of curvature and load vary with time. It is therefore important that the axial surface profile shape should suffice over the full range of operating conditions. Moreover, nowadays an increasing trend in the use of surface coatings in lubricated contacts is observed, and from past studies,^9–12 one may conclude that smart use of surface coatings can significantly enhance lubrication performance. However, according to the authors’ knowledge past studies, concerning lubricated-coated contacts, are rather limited to infinite line contacts and circular contacts.

Therefore, this paper presents a finite element method (FEM)-based finite line contact model that includes the possibility of having a coating on interacting solids. In this work, the lubricated conjunction of an axially surface profiled roller on a plate is analyzed. The numerical predictions are first quantitatively validated with benchmark results found in literature. Furthermore, the influence of operating conditions, roller profiling and coating mechanical properties on the tribological behavior of the contact is investigated.

Mathematical model

The model presented herein is similar to the circular coated contact model presented by Habchi, ¹² but then slightly modified in order to simulate finite line contacts. Furthermore, isothermal conditions are assumed to simplify the analysis. The model relies on a full system finite element resolution of the EHL governing equations, which are the Reynolds, linear elasticity and the load balance equations.

The equivalent computational domain Ω, for the EHL problem, is presented in Figure 1 where above the substrate a coating is present with a unit dimensionless thickness. Furthermore, both coating and substrate share a dimensionless axial length of two. The aforementioned geometrical dimensions, for both coating and substrate, are irrespective of the actual coating thickness $t_{\text{c}}$ and axial length L due to unique definitions of Y and Z (see Notation). As suggested in Habchi et al., ⁸ a dimensionless thickness of 60 for the substrate is sufficient to approximate a half-space for the calculation of elastic displacement field. The two-dimensional Reynolds equation should be solved for the pressure distribution on computational boundary Ω _f , representing the geometrical contact region. The geometrical dimensions for Ω _f are $-4.5\le X\le 1.5$ and $-1\le Y\le 0$ to satisfy zero pressure boundary conditions at the borders. ⁶ Note that in order to reduce computational effort, the advantage of symmetry around boundary Ω _s , which lies in the XZ-plane, has been taken.

Figure 1.

Equivalent geometry for EHL analysis of the finite line coated contact problem. Note that the dimensions are exaggerated for the sake of clarity.

The Reynolds equation, which is a convection-diffusion type equation, is strongly convection dominated in highly loaded situations and will thus exhibit superious oscillations in the solution when solved using a typical Galerkin formulation. ¹³

The inconsistent (non-residual based) artificial diffusion method is one of the oldest and simplest methods as it directly adds an artificial diffusion term to the Reynolds equation in regions where the solution is strongly convection dominated. In Habchi et al., ⁸ it was shown that the solution was not significantly affected with the use of isotropic artificial diffusion. Hence, in this work the inconsistent artificial diffusion stabilization technique is used. The slightly modified Reynolds equation can now be written as follows

$$\frac{\partial }{\partial X}((-\frac{\tilde{\rho }{H}^3}{\tilde{\eta }\lambda }+k_X)\frac{\partial P}{\partial X}+H\tilde{\rho })+\frac{\partial }{\partial Y}((-\frac{a^2}{(2L)2}\frac{\tilde{\rho }{H}^3}{\tilde{\eta }\lambda }+k_Y)\frac{\partial P}{\partial Y})=0$$

(1)

where the dimensionless speed parameter λ is defined as $\lambda =\frac{12u_{\text{m}}{\eta }_0{R}_x^2}{a^3{p}_{\text{h}}},k_{\text{AD},X}$ and $k_{\text{AD},Y}$ are the artificial diffusion coefficients. a and L are the Hertzian contact width and roller axial length, respectively. The variation of viscosity $\tilde{\eta }$ with pressure is simulated using Roeland’s viscosity-pressure relation, ¹⁴ while the density $\tilde{\rho }$ of lubricant is assumed to be dependent on pressure according to the Dowson–Higginson density–pressure relation. ¹⁵

The free boundary cavitation problem that arises at the exit of the lubricated contact is treated according to the penalty formulation of Wu. ¹⁶ This method adds an additional (penalty) term to the Reynolds equation that only acts in the negative pressure zones. The penalty term enforces the arising negative pressure in the solution towards zero. It is important to note that the equivalent diffusion tensor $\{k_X,k_Y\}=\{\frac{\tilde{\rho }{H}^3}{\tilde{\eta }\lambda },\frac{a^2}{(2L)2}\frac{\tilde{\rho }{H}^3}{\tilde{\eta }\lambda }\}$ of equation (1) is anisotropic due to different definitions of X and Y. The stabilizing terms will therefore have to be amended for the anisotropic nature accordingly. For this reason, $k_{\text{AD},X}$ and $k_{\text{AD},Y}$ are defined separately and constructed in the following manner ¹³

$$[k_{\text{AD},X}{k}_{\text{AD},Y}]={\alpha }_0{\alpha }_{\text{AD}}{h}_e|u|[{\rho }_{\text{AD},X}{\rho }_{\text{AD},Y}]$$

(2)

where $\{{\rho }_{\text{AD},X},{\rho }_{\text{AD},Y}\}$ are tuning parameters, h_e is a typical element size and $u=H\frac{\partial \tilde{\rho }}{\partial P}$ is the equivalent convection coefficient. The constant α₀ and upwind function ${\alpha }_{\text{AD}}$ are defined as follows ⁸

$${\alpha }_0=\frac{1}{2l}$$

(3a)

$${\alpha }_{\text{AD}}=\coth ({\mathrm{Pe}}_{\text{mean}})-1/{\mathrm{Pe}}_{\text{mean}}$$

(3b)

where l is the interpolation order. The mean diffusion coefficient $k_{\text{mean}}$ and scaled cell Peclet number ${\mathrm{Pe}}_{\text{mean}}$ (according to the formulation in Gale $\tilde{\text{a}}$ o et al. ¹⁷ ) are computed as follows

$$k_{\text{mean}}=(\frac{1}{k_X}+\frac{1}{k_Y}{)}^{-1}$$

(4a)

$${\mathrm{Pe}}_{\text{mean}}=\frac{{\mathrm{uh}}_e}{2k_{\text{mean}}l}$$

(4b)

As can be extracted from equation (2), the only difference between the definitions of $k_{\text{AD},X}$ and $k_{\text{AD},Y}$ lies in the choice of tuning parameters ${\rho }_{\text{AD},X}$ and ${\rho }_{\text{AD},Y}$ . The tuning parameters should be minimal just to suppress the generated oscillations and not so large to introduce excessive damping.

The EHL film thickness expression H for a straight cylindrical roller with axial dub-off profiling (see Figure 2) can simply be written as Park and Kim ⁴

$$H(X,Y)=H_0+\frac{X^2}{2}+\frac{L^2{R}_x}{4a^2}\frac{(Y+Y_{\text{d}}{)}^2}{2R_{\text{d}}}(Y<-Y_{\text{d}})+\frac{L^2{R}_x}{4a^2}\frac{(Y-Y_{\text{d}}{)}^2}{2R_{\text{d}}}(Y>Y_{\text{d}})-W_{\text{d}}(X,Y)$$

(5)

where H₀ is the rigid body displacement and $W_{\text{d}}$ is the contribution due to elastic deformation. The second term in equation (5) represents the static separation due to the geometry of the roller in undeformed state. Note that $Y_{\text{d}}=\frac{l_{\text{s}}}{L}$ is the dimensionless axial position where axial profiling starts, $R_{\text{d}}$ is the round corner radius and $l_{\text{s}}$ is the straight roller length. $(Y<-Y_{\text{d}})$ and $(Y>Y_{\text{d}})$ are Boolean functions, equal to one if true and zero if not true.

Figure 2.

Schematic of considered roller axial profile.

The conservation law states that the applied load should be balanced by the hydrodynamically generated force. Assuming that acceleration forces are negligible, the following load equation should hold for the contact

$${\int }_{{\Omega }_f}P(X,Y)\text{d}\Omega =\frac{1}{2}\pi$$

(6)

Note that equation (6) already takes into account symmetrical boundary conditions at plane ${\Omega }_{\text{s}}$ . Equation (6) is balanced by iteratively adjusting H₀ until the Reynolds equation, i.e. the pressure solution, converges.

For the elasticity problem three assumptions are made, namely:

both substrates are coated

the substrates of both interacting bodies share similar mechanical properties, i.e. $E_{\text{s},1}=E_{\text{s},2}=E_{\text{s}}$ and ${\nu }_{\text{s},1}={\nu }_{\text{s},2}={\nu }_{\text{s}}$

the coating materials on both substrates also share similar mechanical properties, i.e. i.e. $E_{\text{c},1}=E_{\text{c},2}=E_{\text{c}}$ and ${\nu }_{\text{c},1}={\nu }_{\text{c},2}={\nu }_{\text{c}}$

where subscripts “s” and “c” denote substrate and coating, respectively. Note that the aforementioned assumptions are purely made here to simplify the analysis. Extension to other problems, such as usage of different coating and/or substrate materials, can straightforwardly be taken in to account.

For the elastic deformation calculation, we again make use of an equivalent elastic model, see Habchi et al. ⁸ for more details. In the equivalent elastic model, one of the interacting bodies, thus substrate with coating, is rigid while the other body has equivalent material properties to compensate for the total elastic deformation. As both substrates and coatings share similar mechanical properties, the equivalent dimensionless material properties for substrate ( ${\tilde{E}}_{\text{s},\text{eq}},{\nu }_{\text{s},\text{eq}}$ ) and coating ( ${\tilde{E}}_{\text{c},\text{eq}},{\nu }_{\text{c},\text{eq}}$ ) simplify to

$$\{{\tilde{E}}_{\text{s},\text{eq}}=\frac{E_{\text{s}}}{2}\frac{a}{{\mathrm{Rp}}_{\text{h}}}{\nu }_{\text{s},\text{eq}}={\nu }_{\text{s}}\text{forthesubstrate}$$

(7a)

$$\{{\tilde{E}}_{\text{c},\text{eq}}=\frac{E_{\text{c}}}{2}\frac{a}{{\mathrm{Rp}}_{\text{h}}}{\nu }_{\text{c},\text{eq}}={\nu }_{\text{c}}\text{forthecoating}$$

(7b)

In the current case, since the two substrates and coatings are made of the same material it means that the use of the equivalent material properties defined in equation (7) leads to a total elastic deflection that is twice that of each solid body (substrate + coating) taken individually.

In fact, now two interdependent sets of the system of 3D elasticity equations need to be solved to calculate the elastic displacement field in both coating and substrate. The 3D elasticity equations are applied to dimensionless domain Ω to compute the total elastic deformation. The sets of equations described by eqns.8 and 9 describe the elastic deformation of substrate and coating, respectively. For the substrate we get

$$\zeta \frac{\partial }{\partial X}[({\tilde{\lambda }}_{\text{s}}+2{\tilde{\mu }}_{\text{s}})\frac{\partial U_{\text{d}}}{\partial X}+\psi {\tilde{\lambda }}_{\text{s}}\frac{\partial V_{\text{d}}}{\partial Y}+{\tilde{\lambda }}_{\text{s}}\frac{\partial W_{\text{d}}}{\partial Z}]+\psi \zeta \frac{\partial }{\partial Y}[{\tilde{\mu }}_{\text{s}}(\psi \frac{\partial U_{\text{d}}}{\partial Y}+\frac{\partial V_{\text{d}}}{\partial X})]+\zeta \frac{\partial }{\partial Z}[{\tilde{\mu }}_{\text{s}}(\frac{\partial U_{\text{d}}}{\partial Z}+\frac{\partial W_{\text{d}}}{\partial X})]=0,\zeta \frac{\partial }{\partial X}[{\tilde{\mu }}_{\text{s}}(\psi \frac{\partial U_{\text{d}}}{\partial Y}+\frac{\partial V_{\text{d}}}{\partial X})]+\zeta \psi \frac{\partial }{\partial Y}[{\tilde{\lambda }}_{\text{s}}\frac{\partial U_{\text{d}}}{\partial X}+\psi ({\tilde{\lambda }}_{\text{s}}+2{\tilde{\mu }}_{\text{s}})\frac{\partial V_{\text{d}}}{\partial Y}+{\tilde{\lambda }}_{\text{s}}\frac{\partial W_{\text{d}}}{\partial Z}]+\zeta \frac{\partial }{\partial Z}[{\tilde{\mu }}_{\text{s}}(\frac{\partial V_{\text{d}}}{\partial Z}+\psi \frac{\partial W_{\text{d}}}{\partial Y})]=0,\zeta \frac{\partial }{\partial X}[{\tilde{\mu }}_{\text{s}}(\frac{\partial U_{\text{d}}}{\partial Z}+\frac{\partial W_{\text{d}}}{\partial X})]+\zeta \psi \frac{\partial }{\partial Y}[{\tilde{\mu }}_{\text{s}}(\frac{\partial V_{\text{d}}}{\partial Z}+\psi \frac{\partial W_{\text{d}}}{\partial Y})]+\zeta \frac{\partial }{\partial Z}[{\tilde{\lambda }}_{\text{s}}\frac{\partial U_{\text{d}}}{\partial X}+\psi {\tilde{\lambda }}_{\text{s}}\frac{\partial V_{\text{d}}}{\partial Y}+({\tilde{\lambda }}_{\text{s}}+2{\tilde{\mu }}_{\text{s}})\frac{\partial W_{\text{d}}}{\partial Z}]=0$$

(8)

And for the coating

$$\frac{\partial }{\partial X}[({\tilde{\lambda }}_{\text{c}}+2{\tilde{\mu }}_{\text{c}})\frac{\partial U_{\text{d}}}{\partial X}+\psi {\tilde{\lambda }}_{\text{c}}\frac{\partial V_{\text{d}}}{\partial Y}+\zeta {\tilde{\lambda }}_{\text{c}}\frac{\partial W_{\text{d}}}{\partial Z}]+\psi \frac{\partial }{\partial Y}[{\tilde{\mu }}_{\text{c}}(\psi \frac{\partial U_{\text{d}}}{\partial Y}+\frac{\partial V_{\text{d}}}{\partial X})]+\zeta \frac{\partial }{\partial Z}[{\tilde{\mu }}_{\text{c}}(\zeta \frac{\partial U_{\text{d}}}{\partial Z}+\frac{\partial W_{\text{d}}}{\partial X})]=0,\frac{\partial }{\partial X}[{\tilde{\mu }}_{\text{c}}(\psi \frac{\partial U_{\text{d}}}{\partial Y}+\frac{\partial V_{\text{d}}}{\partial X})]+\psi \frac{\partial }{\partial Y}[{\tilde{\lambda }}_{\text{c}}\frac{\partial U_{\text{d}}}{\partial X}+\psi ({\tilde{\lambda }}_{\text{c}}+2{\tilde{\mu }}_{\text{c}})\frac{\partial V_{\text{d}}}{\partial Y}+{\tilde{\lambda }}_{\text{c}}\zeta \frac{\partial W_{\text{d}}}{\partial Z}]+\zeta \frac{\partial }{\partial Z}[{\tilde{\mu }}_{\text{c}}(\zeta \frac{\partial V_{\text{d}}}{\partial Z}+\psi \frac{\partial W_{\text{d}}}{\partial Y})]=0,\frac{\partial }{\partial X}[{\tilde{\mu }}_{\text{c}}(\zeta \frac{\partial U_{\text{d}}}{\partial Z}+\frac{\partial W_{\text{d}}}{\partial X})]+\psi \frac{\partial }{\partial Y}[{\tilde{\mu }}_{\text{c}}(\zeta \frac{\partial V_{\text{d}}}{\partial Z}+\psi \frac{\partial W_{\text{d}}}{\partial Y})]+\zeta \frac{\partial }{\partial Z}[{\tilde{\lambda }}_{\text{c}}\frac{\partial U_{\text{d}}}{\partial X}+\psi {\tilde{\lambda }}_{\text{c}}\frac{\partial V_{\text{d}}}{\partial Y}+\zeta ({\tilde{\lambda }}_{\text{c}}+2{\tilde{\mu }}_{\text{c}})\frac{\partial W_{\text{d}}}{\partial Z}]=0$$

(9)

where $\{U_{\text{d}},V_{\text{d}},W_{\text{d}}\}$ are the X, Y and Z-components of the solid’s elastic displacement field, $\psi =\frac{a}{2L}$ and $\zeta =\frac{a}{t_{\text{c}}}$ . Equations (8) and (9) are derived analogously to the elasticity equations presented in Habchi. ¹² The dimensionless (equivalent) Lamé’s coefficients for substrate and coating are calculated as follows

$$\{{\tilde{\mu }}_{\text{s}}=\frac{{\tilde{E}}_{\text{s},\text{eq}}}{2(1+{\nu }_{\text{s},\text{eq}})}{\tilde{\lambda }}_{\text{s}}=\frac{{\nu }_{\text{s},\text{eq}}{\tilde{E}}_{\text{s},\text{eq}}}{(1-2{\nu }_{\text{s},\text{eq}})(1+{\nu }_{\text{s},\text{eq}})}\text{forthesubstrate}$$

(10a)

$$\{{\tilde{\mu }}_{\text{c}}=\frac{{\tilde{E}}_{\text{c},\text{eq}}}{2(1+{\nu }_{\text{c},\text{eq}})}{\tilde{\lambda }}_{\text{c}}=\frac{{\nu }_{\text{c},\text{eq}}{\tilde{E}}_{\text{c},\text{eq}}}{(1-2{\nu }_{\text{c},\text{eq}})(1+{\nu }_{\text{c},\text{eq}})}\text{forthecoating}$$

(10b)

where the material properties ( ${\tilde{E}}_{\text{s},\text{eq}},{\nu }_{\text{s},\text{eq}}$ ) and ( ${\tilde{E}}_{\text{c},\text{eq}},{\nu }_{\text{c},\text{eq}}$ ) are evaluated by means of equations (7a) and (7b), respectively.

In order to obtain a unique solution for the EHL problem, proper boundary conditions (BCs) need to be imposed. These are summarized as follows:

For the Reynolds equation

$$\{P=0\text{on}\partial {\Omega }_f\nabla P·\text{n}=0\text{on}{\Omega }_s$$

(11)

For the elastic model

$$\{U_{\text{d}}=V_{\text{d}}=W_{\text{d}}=0\text{on}{\Omega }_D{\sigma }_n={\sigma }_{\mathrm{ZZ}}=\zeta [{\tilde{\lambda }}_{\text{c}}\frac{\partial U_{\text{d}}}{\partial X}+\psi {\tilde{\lambda }}_{\text{c}}\frac{\partial V_{\text{d}}}{\partial Y}+({\tilde{\lambda }}_{\text{c}}+2{\tilde{\mu }}_{\text{c}})\zeta \frac{\partial W_{\text{d}}}{\partial Z}]=-P\text{on}{\Omega }_f{V}_{\text{d}}=0\text{on}{\Omega }_s{\sigma }_n=0\text{elsewhere}$$

(12)

Additionally, a continuity BC on the common boundary of coating and substrate is imposed.

Results

Summarizing, the complete model consists of the modified Reynolds equation (1), the load balance equation (6), two sets of interdependent elasticity equations for substrate and coating equations (8) and (9), and their respective boundary conditions equations (11) and (12).

The dry dimensionless Hertzian pressure profile was taken as initial guess for the pressure distribution, while for the elastic displacement field the solution corresponding to the dry Hertzian contact was taken as initial guess. The load balance equation (6), which is associated with the unknown H₀, is added to the complete system of equations formed by equations (1), (8) and (9).

The developed model is solved using the FEM with a general purpose finite element analysis software. ¹⁸ The resulting system of non-linear equations is solved using a monolithic approach where all the dependent variables $(P,U_{\text{d}},V_{\text{d}},W_{\text{d}},H_0)$ are collected in one vector of unknowns and simultaneously solved using a modified Newton–Raphson iteration scheme. For specific numerical details, concerning the weak coupling resolution of the Reynolds and elasticity equations, the reader is referred to Habchi et al. ⁸ as only the main features of the model are recalled here.

Lagrange quintic elements were used for the hydrodynamic part, while quadratic elements were used for the elastic part. A custom-tailored mesh, similar to that detailed in Habchi, ¹² was deployed to reach an optimum combination between accuracy and calculation speed. The maximum element size in the contact zone in X-direction was chosen smaller than 0.06 while for Y-direction the maximum element size was chosen to be smaller than 0.03. This is because steeper gradients are expected in Y-direction. The element size was allowed to increase gradually as the distance from the contact boundary increased. The aforementioned corresponds to the usage of approximately 350,000 degrees of freedom for the uncoated case and approximately 400,000 degrees of freedom for the coated case. Converged solutions to relative errors ranging between ${10}^{-3}\text{ and }{10}^{-4}$ are typically reached within 12 iterations, corresponding to a computation time of approximately 2 minutes on an Intel(R) Xenon(R) CPU E5-2640 processor. Much less iterations (2) to (5) are required when calculations are continued from a previously obtained solution with a somewhat different set of input parameters.

The results herein are presented in twofold. The first part of the results corresponds to uncoated contacts, while the second part discusses results corresponding to coated contacts. Note that the particular case when the Young’s moduli of coating and substrate are identical, i.e. $E_{\text{c}}=E_{\text{s}}$ , corresponds to an uncoated contact.

Uncoated case

Influence of numerical stabilization on overall solution

The tuning factors $[{\rho }_{\text{AD},X},{\rho }_{\text{AD},Y}]$ are chosen different due to anisotropic nature of the diffusion tensor $\{k_X,k_Y\}$ of equation (1). In fact, in terms of the amount of artificial diffusion added for streamline stabilization, the effect is similar to that described in Habchi et al. ⁸ and therefore not detailed here, i.e. small deviations around the pressure spike are observed and the minimum film thickness is not significantly effected. It is thus safe to choose ${\rho }_{\text{AD}X}$ smaller than 0.5.

It is widely known, for infinite line and point contacts, that with increasing loads the Petrusevich pressure spike(s) shift more towards the exit of the lubricated contact and becomes smaller in magnitude at the same time. The maximum pressure then occurs at the center of the contact and is approximately the same as the maximum dry Hertzian contact pressure $p_{\text{h}}$ . However, for the finite line contact the maximum pressure $P_{\text{max}}$ and minimum film thickness $H_{\text{min}}$ occur near the edges of the contact. Therefore, the amount of artificial diffusion in crosswind direction needs to be chosen carefully as this also is the direction that causes the anisotropic nature of the diffusion tensor.

Figure 3 shows the pressure distribution for a straight roller with rounded edges with and without stabilization. The operating conditions and roller profiling geometrical parameters are given in Table 1. Note that the Young’s moduli of coating and substrate are identical, i.e. $E_{\text{c}}=E_{\text{s}}$ , which corresponds to an uncoated contact. As can be retrieved from Figure 3, the pressure distribution without crosswind stabilization is not smooth whereas with crosswind stabilization a smooth solution is obtained.

Figure 3.

Height expressions for the pressure distribution (a) with crosswind stabilization and (b) without crosswind stabilization.

Table 1.

Reference operating conditions and geometrical parameters for a straight roller with rounded edges.

Parameter	Value	Unit
F	3150	N
$E_{\text{s}}$	210	GPa
${\nu }_{\text{s}}$	0.3	–
$E_{\text{c}}$	210	GPa
${\nu }_{\text{c}}$	0.3	–
α	1.89E−8	Pa– 1
η 0	0.013	Paċs
$U_{\text{m}}$	1	m/s
Rx	0.008	m
L	0.01	m
$R_{\text{d}}$	0.127	m
$l_{\text{s}}$	0.0085	m
U	7.3891E−12	–
W	1.7904E−4	–
G	4150	–
$p_{\text{h}}$	1.17	GPa

From “numerical experiments” we ascertained that the influence of crosswind artificial diffusion on the absolute minimum film thickness is much more amplified when compared to the influence of streamline artificial diffusion. This is mainly due to the fact that axial pressure gradients $\frac{\partial P}{\partial Y}$ are affected in a much more amplified fashion with crosswind diffusion due to the anisotropic nature of the diffusion tensor. This is also the reason why crosswind diffusion has a greater influence on the absolute minimum film thickness $H_{\text{min}}$ as compared to the central plane minimum film thickness $H_{\text{min},\text{central}}$ , which seems to remain unaffected (see Figure 4(a)). The exact value for $H_{\text{min}}$ corresponds to a value of ${\rho }_{{\text{AD}}_Y}=0$ , but then the solution for pressure distribution is not smooth. This can be extracted from Figure 4(b) in which the pressure and film thickness distributions are plotted along the line X = 0. Note that the absolute maximum pressure $P_{\text{max}}$ and the maximum pressure at the central plane $P_{\text{max},\text{central}}$ are negligibly affected by the introduced amount of crosswind diffusion (see Figure 4(c)). Nevertheless, from “numerical experiments” we conclude that, as a rule of thumb, ${\rho }_{\text{AD},Y}$ should always be chosen smaller than 0.1.

Figure 4.

Influence of crosswind stabilization on (a) minima film thicknesses and (b) on axial pressure and film thickness distributions and (c) maxima pressures.

The present authors also attempted to implement more consistent (residual based) stabilizing methods, such as described in Codina ¹⁹ and Do Carmo and Galeão, ²⁰ to stabilize the solution in crosswind direction. These include shock-capturing techniques which aim to eliminate effects such as overshoots and undershoots close to discontinuities. Unfortunately, these techniques make the system of equations more non-linear and consequently induce convergence issues. This is of course an incentive for more detailed investigation into the implementation of more consistent techniques for the “anisotropic convection-dominated convection diffusion problem.”

Quantitative validation

Park and Kim ⁴ have presented benchmark results for an uncoated straight roller (with straight length $l_{\text{s}}$ ) with rounded edges (with dub-off radius $R_{\text{d}}$ ). They also compared their results qualitatively with experimental results obtained using optical interferometry. ²¹ The operating conditions and roller profile parameters are given in Table 2.

Table 2.

Operating conditions for reference case. Partly adapted from Park and Kim. ⁴

Parameter	Value	Unit
F	570.24	N
$E_{\text{s}}$	200	GPa
${\nu }_{\text{s}}$	0.3	–
$E_{\text{c}}$	200	GPa
${\nu }_{\text{c}}$	0.3	–
α	1.364E−8	Pa– 1
η 0	0.0528	Paċs
$U_{\text{m}}$	1	m/s
Rx	0.012	m
L	2 × 0.7 × Rx × 1.07	m
$R_{\text{d}}$	0.3 × Rx	m
$y_{\text{d}}$	0.7 × Rx	m
$l_{\text{s}}$	2 ×  $y_{\text{d}}$	m
U	2E−11	–
$W_{\text{Kim}}=\frac{F}{E\text{'}{R}_x^2}$	1.8E−5	–
G	3000	–
$p_{\text{h}}$	0.304	GPa

Note that the presented coated finite line contact model herein can be numerically validated with the results presented in Park and Kim ⁴ if the material properties of coating and substrate are set to be the same (see Table 2).

In Park and Kim ⁴ a slightly different definition for the load parameter ( $W_{\text{Kim}}=\frac{F}{E\text{'}{R}_x^2}$ ) is used than what is “usual” for infinite line contacts ( $W=\frac{F/L}{E\text{'}{R}_x}$ ). The axial length L is not given in Park and Kim, ⁴ but can somehow be estimated. In the present analysis, the axial length L was estimated on the basis of the contact footprint length, i.e. for an equivalent infinite line contact problem the load acting over the unit footprint length can be used as input data. ²² In Park and Kim, ⁴ the profiling starts at a position of $y_{\text{d}}=0.7$  × R_x from the central plane, meaning that for the infinite analysis $L=l_{\text{s}}=2$  × 0.7 × R_x would be used as input data. This statement does not always hold as for higher loads situations the lubricated contact footprint length becomes larger. This will be shown in the next subsection. For dry contact analysis, the footprint length is usually taken to be the same as the straight roller length $l_{\text{s}}$ . However, for a lubricated contact, the pressure distribution is slightly extended. ⁶ For this reason, the axial length here is assumed to be a factor of 1.07 times larger than the contact dry footprint length in Park and Kim, ⁴ i.e. $L=2$  × 0.7 × R_x × 1.07.

Comparisons are made according to different sections of the contact. These correspond to the streamline and axial sections through the maximum pressure and absolute minimum film thickness. The sections are defined as shown in Figure 5.

Figure 5.

Definition of sections through lubricated contact. Reproduced from reference Park and Kim. ⁴

Section 1-1 is the central line in streamline direction (plotted against Y = 0) where the minimum film thickness $H_{\text{min},\text{central}}$ occurs. Sections 2-2 and 3-3 correspond to contour sections where the absolute maximum pressure $P_{\text{max}}$ and the absolute minimum film thickness $H_{\text{min}}$ occur, respectively. Sections 4-4 and 5-5 are contour sections in transverse direction where the absolute maximum pressure and absolute minimum film thickness occur, respectively.

Figure 6 presents the zoomed-in contour plots for pressure and film thickness. Note that different dimensionless variables are used for the sake of comparison. It is clear that the maximum pressure and minimum film thickness occur near the region where profiling starts (thus near a dimensionless position of $y/R_x=-0.7$ ). In fact the secondary pressure peak occurs near the side constriction and rear exit of the lubricated contact area. Flow continuity demands that the pressure gradients should be coupled with local restrictions of minima film thickness. Hence, the secondary pressure peaks near the side constriction also inhibit lubricant flow in their vicinity. Consequently small islands (iso-film thickness contours) are formed at the rear of the contact. These are commonly referred as end closure films in literature. ⁶

Figure 6.

Zoomed-in contour plots of (a) the film thickness and (b) pressure distribution at the rear of the contact. Note that different dimensionless variables have been used for the sake of comparison.

It can readily be concluded from the figures that traditional EHL infinite line contact solutions do not reveal the tribological behavior at the extremities of the contact in terms of pressure and film thickness distributions. These findings are in line with previous studies.^4,6,23

In Figure 7, the results for pressure and film thickness distributions are presented according to the defined sections in Figure 5. The results are to be compared with those reported in Park and Kim. ⁴ Overall, good agreement is obtained for the results obtained using the current approach. The minimum film thickness and maximum pressure and their positions are accurately predicted. In Figure 7(e), there is some minor difference in obtained axial pressure distribution for Section 4-4, although the maximum pressure is accurately predicted. This is mainly due to the assumed axial length as earlier discussed. The assumed axial length might shift the location of maximum pressure a little, which in turn may lead to this discrepancy.

Figure 7.

Results for pressure and film thickness profiles, for the different sections, using the current approach (left column). Note that here $H=\frac{h}{R_x}$ and $P=\frac{p}{E\text{'}}$ . The figures are to be compared with those reproduced from Park and Kim ⁴ (right column).

Parameter study

It is of interest to see how the finite line contact responds to varying operating conditions such as load, speed and material properties. The dimensionless parameters of the aforementioned operating conditions are represented by the load parameter W, speed parameter U and material property parameter G. The operating conditions for the current cases are detailed in Table 1.

Figure 8 shows how the pressure and film thickness distribution (plotted along X = 0 and Y = 0) vary as a function of these three parameters, while keeping two constant at a time. Variation of W, U and G all result to some minor variation around the pressure spike along the central plane and/or central film thickness variation. These variations are much more explainable from traditional EHL solutions for line and/or elliptical contacts. The most remarkable observation is that the pressure and film thickness distribution, especially at the side extremities, are highly affected by varying loads. To be more specific, the variation in pressure distribution at the extremities seems to be highly amplified with increasing load in the sense that the secondary pressure peak smears out and the covered (lubricated) axial length becomes larger. These results are in line with previous theoretical ²¹ and experimental findings. ⁴

Figure 8.

Influence of varying operating conditions such as speed, elasticity and contact load on streamline and axial pressure shapes.

To make things more clear, Figure 9 plots the variation $H_{\text{min}}/H_{\text{min},\text{central}}$ as function of W, U and G from which it is clear that the absolute minimum film thickness is highly affected by variations in load. Especially from low to moderate loads, this phenomenon is much more visible. For variations in U and G, the ratio $H_{\text{min}}/H_{\text{min},\text{central}}$ seems to remain constant. Note that the behavior of $H_{\text{min},\text{central}}$ is much more explainable using traditional EHL solutions for infinite line contacts. ⁴ It is therefore much more interesting to study the behavior of ratio $H_{\text{min}}/H_{\text{min},\text{central}}$ from a designers perspective. In practice, one would like to maximize the value of $H_{\text{min}}/H_{\text{min},\text{central}}$ as $H_{\text{min},\text{central}}$ can fairly be estimated using the infinite line contact assumption. This would drastically decrease computational overhead. Figure 10 plots the contours of the film thickness for increasing loads. It can be seen that for all conditions $H_{\text{min}}$ occurs at the rear of the contact and that for increasing loads the lubricated area extends (also see Figure 8(e)). Also note that also, similar as for elliptical contacts, with increasing load the contact center gets more flattened and the end closure films (small islands of minimum film thickness) become smaller. Consequently, the ratio $H_{\text{min}}/H_{\text{min},\text{central}}$ also increases as can be seen from Figure 9(a). This however also depends on the axial profile itself, i.e. the straight roller length and dub-off radius as will be explained now.

Figure 9.

Variation of $H_{\text{min}},H_{\text{min},\text{central}}$ and ratio $H_{\text{min}}/H_{\text{min},\text{central}}$ with dimensionless speed, material and load parameters while keeping two fixed at a time. The minimum film thickness $H_{\text{min}}$ seems to behave different with increasing load.

Figure 10.

Contour plots of film thickness H for increasing values of dimensionless load parameter W showing how the position of minimum film thickness shifts as the contact center flattens.

So, apart from varying operating conditions, it is also interesting to take a look at the influence of geometrical parameters on the pressure and film thickness distributions. In fact, for the axial profile of the roller one may vary the straight roller length and dub-off radius to optimize the pressure distribution, i.e. to make it more uniform by reducing edge stress concentrations and consequently increase $H_{\text{min}}/H_{\text{min},\text{central}}$ .

Figure 11 presents the variations of pressure and film thickness profiles as a function of the dub-off radius $R_{\text{d}}$ . From Figure 11(a), it is clear that a higher relief radius smears out the secondary pressure peak, resulting to a larger contact area. Furthermore, the ratio $H_{\text{min}}/H_{\text{min},\text{central}}$ seems to increase with increasing $R_{\text{d}}$ .

Figure 11.

Influence of round corner radius on (a) pressure and film thickness distributions and (b) minima film thicknesses.

The aforementioned is amplified with increasing loads. One would then think that choosing a larger $R_{\text{d}}$ results in a more uniform the pressure profile and thus a better design. However, there seems to be an optimum range for minimum film thickness versus dub-off radius mapped against the range of loads. In fact, for a too large $R_{\text{d}}$ the ratio $H_{\text{min}}/H_{\text{min},\text{central}}$ starts to decrease after a certain applied load. This is mainly due to the fact that there is no space available for the pressure profile to extend as a zero boundary condition is imposed at the extremities. Consequently, the pressure gradient $\frac{\mathrm{dP}}{\text{d}Y}$ at the extremities increases and thus $H_{\text{min}}$ decreases (see Figure 11(b)).

A similar statement can be made about variations in the straight length of the roller (see Figure 12(a)). There is an optimum range for choosing an appropriate straight roller length as a too large value for $Y_{\text{d}}$ results to a decrease of $H_{\text{min}}/H_{\text{min},\text{central}}$ with increasing loads (see Figure 12(b)). If a too small value for $Y_{\text{d}}$ is chosen, the contact area is reduced and the maximum pressure consequently increases.

Figure 12.

Influence of roller straight length on (a) pressure and film thickness shape and (b) minima film thicknesses.

All these findings make it really hard to develop a robust correlation between absolute minimum film thickness and operating conditions, as axial profile design parameters play an equally important role. In fact, a good understanding of finite line contact behavior as a basis will lead to a better design of mechanical components in terms of film thickness and pressure distributions, and as a result, increase in service life.

Coated case

The influence of mechanical properties, in terms of coating stiffness and thickness, on the overall EHL behavior of a finite line contact will be discussed in this section. The reference operating conditions for the present results are identical to those presented in Table 1. In addition the coating thickness $t_{\text{c}}$ and Youngs modulus $E_{\text{c}}$ will subsequently be defined for the cases studied. Substrate mechanical properties are kept fixed and are given in Table 1.

In Figure 13(a) and (b), the influence of coating stiffness on the pressure and film thickness distribution is illustrated. Note that the case when $E_{\text{c}}$  = 210 GPa (gray line), corresponds to the uncoated contact case.

Figure 13.

Influence of coating elasticity $E_{\text{c}}$ and thickness $t_{\text{c}}$ on pressure and film thickness distribution. The graphs are plotted on lines Y = 0 (left column) and X = 0 (right column). The contact area increases for softer coatings. This phenomenon is amplified with coating thickness. Exact the opposite occurs for stiffer coatings with increasing coating thickness.

One can directly observe that with increasing stiffness of the coating the maximum pressure increases, while the contact width decreases. Furthermore, the pressure spike at the central plane and the secondary pressure peak at the rear of the contact increase in magnitude with increasing coating hardness. It seems obvious that with stiffer coatings higher pressures are expected, and to compensate for this (with fixed applied load), the contact area has to decrease. Taking a more detailed look at the behavior of minima film thicknesses, $H_{\text{min}}$ and $H_{\text{min},\text{central}}$ depicted in Figure 14(c), it is clear that over the whole range of coating stiffness the minimum film thickness increase slightly (less than 10%) with increasing coating stiffness. Obviously, stiffer coatings mean less deformation and thus slight increase in minimum film thickness.

Figure 14.

Influence of (a) coating thickness $t_{\text{c}}$ , (b) dimensionless load parameter $\overline{W}$ and (c) coating elasticity on minima film thicknesses.

Now taking a look at Figure 13(c) to Figure 13(f), we see that the aforementioned phenomena, i.e. higher pressure with increasing coating elasticity, are amplified with increasing coating thickness. To be more specific, we see that for elastic coatings ( $E_{\text{c}}$  = 70 GPa), the maximum pressure is further reduced with increasing coating thickness (see Figure 13(c) and (d)), while for stiff coatings ( $E_{\text{c}}$  = 410 GPa) the maximum pressure is further increased with increasing coating thickness (see Figure 13(e) and (f)). To elaborate a bit on this, if we have a really elastic coating, for example $E_{\text{c}}$  = 70 GPa, and we decrease the thickness of this coating, the influence of the stiff substrate ( $E_{\text{s}}$  = 210 GPa) increases. This means that there will be less deformation, i.e. the contact area decreases, thus the pressure will increase. Exactly the opposite occurs when considering a very stiff coating (e.g. $E_{\text{c}}$  = 410 GPa).

In line with previous findings, for coated contacts, we also see that for more elastic coatings the lubricated contact area is increased. This effect is further amplified with increasing coating thickness. The same also applies for increasing contact force, i.e. the lubricated contact area is also expanded with increasing load. Careful attention should be paid when dealing with elastic and thick coatings and high loads, in terms of minimum film thickness $H_{\text{min}}$ , as all three aforementioned factors lead to an increase in contact area. From Figure 14(b), it is clear that at relatively high loads the ratio $H_{\text{min}}/H_{\text{min},\text{central}}$ dramatically decreases. This is mainly due to the fact that at the rear of the contact, the pressure distribution does not have sufficient space to expand, and thus the secondary pressure peak again grows in magnitude. Consequently, pressure gradients at the extremities increase and thus will the minimum film thickness $H_{\text{min}}$ decrease. The influence of coating thickness, for elastic and stiff coatings, on $H_{\text{min}}$ and $H_{\text{min},\text{central}}$ , are depicted in Figure 14(a). At first sight, it can readily be concluded that the trends of minimum film thicknesses, for soft and stiff coatings, respectively, are opposite. In fact, for soft coatings, the minimum film thicknesses increase with increasing coating thickness, while for hard coatings, the minimum film thicknesses slightly decrease with increasing coating thickness. This behavior is not as expected as we saw that the minimum film thicknesses increased with increasing coating stiffness (refer to Figure 14(c) again). One would expect that for example, given coating which is harder than the substrate, the minimum film thickness would increase if the coating thickness is increased. Again, as can be extracted from Figure 14(a), this is not true due to the fact that when coating thickness is increased, the behavior of the minima film thicknesses are more likely governed by the pressure gradients ( $\frac{\mathrm{dP}}{\text{d}X}$ and $\frac{\mathrm{dP}}{\text{d}Y}$ ) at the exits as per flow continuity demand. Meaning that for hard coatings, the pressure gradients increase with increasing coating thickness, resulting to a decrease in film thickness. Similarly, for more elastic coatings the pressure gradients decrease with increasing coating thickness resulting to an increase in minimum film thicknesses.

Conclusions

In the present work, a FEM-based coated finite line contact model was developed. The lubricated conjunction between a roller with rounded edges on a plate was analyzed.

The developed model was quantitatively validated by means of representative results reported in literature. Good agreement between the results was obtained.

Parameter studies were carried out to investigate the influence of operating conditions, geometrical parameters (of axial surface profile) and coating mechanical properties on the overall EHL behavior of the contact. In line with previously reported findings, it is shown that the pressure and film thickness distributions for finite line contacts vary significantly different with applied load as compared to infinite line contact models. At increasing loads, the pressure distribution becomes more uniform in axial direction as long as there is space available for contact area expansion. When no space is left to compensate for higher loads, the secondary pressure peak at the extremities grows again and hence the absolute minimum film thickness decreases as per flow continuity demand.

Large round corner radii, large straight roller lengths, too elastic and thick coatings, all amplify the effect of increasing loads, i.e. the pressure profile expands in all cases to compensate for the applied load. When no space is left to compensate for higher loads, the secondary pressure peak at the extremities grows again and hence the absolute minimum film thickness decreases.

All these findings make it really hard to develop a robust correlation between absolute minimum film thickness, maximum pressure and operating conditions, as coating and axial profile geometrical parameters play an equally important role. The present results certainly contribute to a better understanding of lubricated and coated finite line contacts. This model can effectively be used for improved designs of finite line contact applications in terms of film thickness and pressure distributions, and thus ultimately, contributing to longer service life of the components.