The Effect of Working Parameters upon Elastohydrodynamic Film Thickness Under Periodic Load Variation

There are a number of widely used machine components, such as rolling element bearings, gears and cams, which operate in the lubrication regime known as Elastohydrodynamics (EHD), where lubricant film thickness is governed by hydrodynamic action of convergent geometry, elastic deformation between non-conformal contacting surfaces, and the increase of lubricant viscosity with pressure. Variable loading conditions occur not only in all the machine components mentioned above, but also in natural joints such as hip or knee joints of humans or many vertebrates. Experimental studies of the behaviour of EHD films under variable loading are scarce and to authors’ knowledge systematic studies of the evolution of lubricant film thickness in EHD contacts subjected to forced harmonic variation of load are even less common. The aim of the present study is to explore the effect of load amplitude on the EHD film behaviour. This is done in alternating cycles with the load varying about a fixed, preset value at various amplitudes. Experimental results are compared with a simple theoretical analysis based on the speed of change of contact’s dimensions, a semi-analytical solution which includes both speed variation and squeeze effect, and finally with a full numerical solution.


Introduction
It is well known that the formation of EHD films is determined by the mechanisms which take place in the contact inlet region according to the classical lubrication theory. In the meantime, EHD films most often work under transient condition of varying load, geometry or speed. These make the behaviour of the EHD film different from steady-state conditions. Classical EHD theory indicates that central film thickness h c increases with a product of entrainment speed U and dynamic viscosity , e.g. h c ∝ (U ) 0.67 with the assumption of fully flooded conditions in the contact inlet zone [1]. Details about normal contact of elastic solids and rolling contact of elastic bodies are addressed in [2]. It has to be agreed that the formation and behaviour of EHD films are fairly well understood in steady-state conditions from both theoretical and experimental point of view [3][4][5][6][7]. The influence of transient conditions of speed or load upon the behaviour of the EHD contacts has been investigated in the past, both experimentally and theoretically [8][9][10]; however, experimental studies dedicated to the effect vibrations upon the formation of the EHD film are few. Glovnea and Spikes [11] have shown that rapid change of speed induces fluctuations on the EHD film thickness while Wijnant et al. [12] showed that the same effect is caused by a rapid variation of load. The amplitude of these fluctuations decreases with time, after a few cycles, very much alike the dynamic response of a mass/spring/damper system. Kilali et al. [13] reported a new test rig capable of loading an EHD contact both as an impact load and as a continuously variable load. Ciulli and Bassani [14] also used optical interferometry to study the film thickness in a contact subjected to random vibrations. Cann and Lubrecht [15] studied the effect of lowfrequency loading/unloading upon the track replenishment of a rolling EHD contact. They noticed the beneficial effect of these conditions upon the lubricant film thickness. Nagata et al. [16] carried out film thickness measurements in an EHD contact subjected to lateral vibrations, that is, along a direction contained in the plane of the film. They tested greases and their base oils and found that lateral vibrations help to replenish the EHD contact.
No systematic experimental study on the effect of load amplitude on EHD films has been published before. Therefore, in the present study, the authors have conducted experiments on the measurement of lubricant film thickness in a EHD contact subjected to forced harmonic vibration; the effect of load amplitude on EHD films is subsequently extracted from the experimental results and compared to theoretical analyses.

Experimental Setup and Method
The experimental method for measuring film thickness in the present study is the optical interferometry technique. This technique was firstly used by Cameron and Gohar [4] and then extensively applied by many other researchers. There are currently a number of techniques for measuring the lubricant film thickness in elastohydrodynamic contacts based on optical interferometry developed in various laboratories around the world: Ultra-thin film interferometry (UTFI) [17], Spacer Layer Imaging Method (SLIM) [18], Relative Optical Interference Intensity (ROII) [19], and differential colorimetry [20].
The principle of optical interferometry, as it is used in this research, to evaluate the EHD film thickness is seen in Fig. 1. The EHD contact under study is formed between a glass disc and a super-polished steel ball. The contacting surface of the disc is coated with a semi-reflective chromium layer (5-10 nm thick) plus a spacer layer (approximately 150 nm thickness). White light source is used to generate film thickness maps in the form of coloured interferometry images. The optical interferometry technique shown in Fig. 1 is incorporated in a test rig which is capable to produce harmonic vibrations normal to the EHD contact. A schematic of the test rig is seen in Fig. 2. The contacting disc is attached to a shaft driven by an electrical motor. The contacting ball is attached to a shaft which is supported by suitable bearings placed in a carriage. The ball is free to rotate by the friction force generated within the EHD contact.
The ball carriage is supported by a bellow which allows the oscillatory motion to be transmitted to the ball, but avoiding leakage of the lubricant from the working chamber. The sinusoidal force acting normal to the EHD contact was precisely measured by a load cell placed on a plunger which drives the ball carriage. The electrical signal of the load cell was recorded on a memory oscilloscope which was triggered at the same time with the CCD camera such that the precise correspondence between the load value and the interference frames was obtained. An additional ball and carriage are placed under the disc on the side opposite to the test ball to minimise the deflection of the disc under loading.
The detailed procedures of images' recording and film thickness calibration are provided in [21]. The lubricant used for testing in this paper is a synthetic hydrocarbon PAO 40 (Poly-Alpha-Olefin) oil, with kinematic viscosity of 339.8 mm 2 /s (339.8 cSt) at 40 °C and 35.8 mm 2 /s (35.8 cSt) at 100 °C. The density at 40 °C is 840 kg/m 3 . Throughout all experiments, the temperature of the lubricant chamber was controlled and kept at 40 ± 1 °C by suitable fitted thermocouples, while the entrainment speed was set at 0.1 m/s. The frequency of the harmonic load was 60 Hz.
The load was varied on alternating fashion where it changed harmonically about a preset mean value of 30 ± 1 N, at various amplitudes, as shown in Fig. 3. The amplitude was changed between tests at 13 N, 20 N and 29 N. For the maximum amplitude of the alternating cycles, load varies from 5.4 N to 62 N which leads to a Hertzian pressure varying from about 0.3 GPa to about 0.7 GPa. Figure 4 shows interferometry images of the contact, at various amplitudes, during the alternating cycles of load variation.

Experimental Results-Alternating Cycles
Central film thickness measured during these cycles together with the theoretical film thickness calculated for the same time instants using the corresponding values of load are seen in Figs. 5, 6, and 7. The error bars indicate the spread of the experimental results extracted from approximately 60 cycles of harmonic loading, while the markers indicate the averaged value. The steady-state theoretical film thickness was calculated with the wellknown Hamrock-Dowson relationship [1].
In this equation the following common notations were used: U-entrainment speed, that is average speed of the surfaces in the direction of rolling x, η 0 , η-viscosity and pressure-viscosity coefficient of the lubricant at ambient pressure, respectively, W-load, E*-reduced elastic modulus of the solids E * = 2 1 − 2 It can be seen that the measured film thickness is larger than the theoretical, steady-state value for the load increasing phase and smaller in the load decreasing phase. This result is somehow expected, from experiments carried out with shock loading, and can be attributed to the squeeze effect in the fluid film [e.g. 10,12]. It would also be expected that the larger the amplitude of the load variation, the larger the deviation from the steady-state condition.

Fig. 2 Schematic of experimental rig
During the load increasing phase, the central film thickness at four instants of time is shown in Fig. 8. It can be seen that the transient film thickness depends on the amplitude of the load, increasing with this parameter.

Theoretical Analysis Based on Effective Entrainment Speed
The rapid variation of the load, accompanied by a similar change of the contact dimensions, has two effects: squeezing of the fluid, trapping it inside the contact due to the rapid increase of the load, and altering the entrainment speed due to rapid change of the contact dimensions. This approach is supported by the numerical analysis of the dynamic behaviour of elastohydrodynamic contacts due to Hooke [22] who states that the increase of the contact dimensions, due to load variation, is equivalent to an increase of the effective entrainment speed. There are full numerical solutions which are able to generate film thickness taking into account these effects; however, those require considerable time, firstly to build up the numerical code and secondly to run it. In [23] the authors present a simple analytical solution for a sinusoidal varying load based on the effective entrainment speed due to change of contact dimensions. Only the main steps are shown here. If the terms not including the speed and the load are grouped in a constant factor termed A, then the film thickness given by Eq. (1) becomes as follows: Consider the load varying in a sinusoidal fashion, W(t) = W m + W 0 sin t where W m is the mean load, W 0 is the amplitude of the cycle and ω is the circular frequency of the loading cycle. As the load increases, the contact radius increases rapidly which changes the entrainment speed as the contact edge moves in opposite direction to the entrainment. Denoting the contact radius with a, the effective entrainment speed is given by the following equation: where u 0 is the absolute velocity of the surfaces of the ball and disc and da/dt is the relative velocity that is given by the speed of the change of the contact radius. Introducing this entrainment velocity into Eq. (2), the transient film thickness can be written as follows: The radius of the Hertzian, circular contact, which can be assumed equal to the EHD contact radius, is given by relationship (5) [2]: The contact radius is obviously a function of time as the load is a function of time as seen above. Carrying out the time derivation and introducing in (4), after simple re-arrangements, the transient film thickness is given by Eq. (6).
In this equation R is the reduced radius of curvature of the contact (in this case the radius of the ball). The term Au is the constant film thickness corresponding to a load equal to the amplitude of the cycle; thus, the term on the left-hand side of relationship (6) is the non-dimensional, transient film thickness. It was preferred to normalise the transient film thickness in this way because it is easier to compare to the experimental values. Figure 9 shows the ratio of the theoretical transient and steady-state film thickness for a frequency of 60 Hz, 0.1 m/s, for load variation similar to that shown in Fig. 4, that is, the mean load of 30 N and amplitudes of 29 N, 20 N and 13 N.
As it can be seen in this figure, the larger the load amplitude, the larger the departure of the film thickness from the steady-state value. This is qualitatively similar to the behaviour shown in Figs. 5, 6, and 7.
A direct comparison between theoretical steady-state (with variable load) and transient film thickness is illustrated in Fig. 10, for the three cycles of the load variation. The frequency is 60 Hz, the load amplitude 29 N, and the mean load is 30 N, while the entrainment speed is 0.1 m/s in these simulations. The load is normalised by the amplitude. The simplified analysis is not intended to give a full, quantitative prediction of the film thickness; however, a qualitative comparison can be made. During the load increasing phase, the transient film thickness exceeds the steady-state values due to the large effective entrainment speed; during the load decreasing phase, the contact dimensions decrease rapidly making the effective entrainment speed smaller than the set value, thus resulting in a smaller transient film thickness. Fig. 8 Central film thickness variation, mean load of 30 N and alternating load cycles

Theoretical Analysis Based on Squeeze and Entrainment
The changing entrainment speed at the entrance of the contact due to the Hertzian width variation when the load fluctuates is in reality only one of the two important effects determining the film thickness variation in the EHD contact, as discussed in [24,25]. The other effect is the squeeze film influence introduced by the time-dependent Reynolds equation. For the sake of simplicity a line contact is assumed, its behaviour might only be representative of what happens in the centre of the point contact when it comes to film fluctuation propagations, as discussed in [26]. So the governing Reynolds equation is, where u e0 represents the entrainment speed in the contact under stationary conditions. Now, assuming that there is a fluctuation of the film thickness given by the function, where (t) is the normal approach fluctuation of the two surfaces and v(x,t) is the elastic displacement associated with the resulting pressure. This pressure should always equate the fluctuating load, Equations (7-9) represent the full EHD problem to be solved. This problem can be solved full numerically or in a semi-analytical way. For the sake of simplicity, let us first start with a semi-analytical solution as described in [24]. Since the complete semi-analytical solution has been described elsewhere, here, only a summary is presented.
The following dimensionless groups are introduced, where q = 1 1 − exp (− p) which represent the reduced pressures introduced from an Ertel-Grubin analysis. Thus Eq. (7) can be written as, Equation (10) is then integrated using the boundary condition that for X = 1 then Q = 1/K 1 with K 1 = αp h , equivalent to the Ertel-Grubin condition p = ∞ . After which a second boundary condition is introduced, for X = L then Q = 0, with L ≈ ∞ . The result is, Discretising the time derivative for one time step as, then Eq. (11) can finally be written in an implicit scheme as Equation (12) can be solved analytically for the inlet film thickness H * by following the work of Wolveridge et al. [27].
Finally, as discussed in [25], in reality the constant K is not constant and varies in time because a varies in time and that affects the entrainment speed also; thus, u e0 is replaced by u e (t). As described in [25], this solution is limited to small amplitude variations in load and contact size. Notice that Hooke and Morales-Espejel [28] developed an improved squeeze film solution that covers much wider application range.  Figure 11 shows results for the changes in central film thickness as function of sinusoidal fluctuations in the applied contact force. For completeness, also included are the results obtained by solving Eqs. (7)-(9) fully numerically employing multi-grid and multi-level integration techniques as described in Appendix A. Both the film thickness and the force variations are normalised by the corresponding values at stationary conditions. As previously observed in the measurements, the effect of including the time-dependent squeeze film term is reflected in these results noting that when increasing the load amplitude, the minima of the film fluctuations deviates less from the stationary values when compared to the maximum values.

Figures 12 shows a comparison of the experimental results
with the results obtained by three sets of theoretical results that is, the simple analytical solution, the semi-analytical solution and the full numerical multi-grid multi-level integration method. The full numerical approach and the semianalytical model do not need explicitly the value of the effective entrainment speed (as defined in Sect. 5.1) and its effect is already considered in the formulation, they only need the steady-state effective speed. But it is clear from good agreement in the comparison between the experiments and the two models that if the time-dependent variation of the contact dimension is ignored the model will be incorrect. The analytical solution based only on the effective entrainment velocity follows the general trend, but gives large deviations from the experiments especially during the load decreasing phase. The deviation is about 9 percent at lowest load amplitude and up to 29 percent at the largest load amplitude employed in these experiments. During the load increasing phase, it shows relatively good agreement. The fully numerical results follow closely the experimental results, within five percent over the whole cycle, for all three amplitudes.
Finally, the semi-analytical method shows very good agreement with experiment and also with the full numerical values. It underestimates both the measured and fully numerical film thickness at the largest load amplitude, at the beginning of the loading cycle. However, this is to be expected since the semi-analytical method applies to relatively small entrainment velocity fluctuations as explained in [25] corresponding to either relatively low-load amplitudes and/or low-load frequencies. For the largest load amplitude, with the nomenclature of reference [25], A u ≈ 0.31 is still Fig. 12 Comparison between experimental and theoretical central film thickness ▸ lower than the applicability limit of 2/3; thus, the semi-analytical method can still be applied for this case despite its slightly lower accuracy.
This comparison shows that the squeeze effect plays the dominant role on the EHD film behaviour during the decreasing phase of the load cycle, while during the load increasing phase, squeeze and effective entrainment play equal roles.

Conclusions
This paper fills the gap in further understanding the effect of load amplitude on the behaviour of EHD film. Film thickness was successfully measured by optical interferometry technique based a ball-on-disc arrangement experimental rig, which is capable of incorporating harmonic load variation exerted normal to the EHD contact plane.
Regarding the experimental results, it was found that the measured film thickness is larger than the theoretical, steadystate values for the load increasing phase and smaller during the load decreasing phase of the loading cycle. It was also found that the larger the amplitude of the load variation, the larger the deviation from the steady-state condition.
The experimental findings were successfully compared with theoretical predictions from three different approaches: a simple transient model based on the variation of effective entrainment speed due to contact dimension changes, a semianalytical model including squeeze and effective entrainment and finally full numerical analysis. It was observed that both the theoretical transient and measured film thicknesses deviate from the steady-state values, but the squeeze film effect makes the real film thickness vary less than the theoretically predicted thickness by the transient analysis.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A. Numerical Solution
Equations (7)-(9) describe a typical line EHD contact solved fully numerically using a well-established iterative procedure [29]. A detailed description of the method can be found in the reference; therefore, it is not repeated here. Only the most relevant details are included. Namely, the equations were made dimensionless using the following groups: The derivatives in Reynolds equation were approximated using second-order finite differences for the pressure term and second-order backward finite differences for the film thickness terms. The numerical solution process made use of multi-grid relaxation methods for improved convergence and of multi-level integration techniques for fast evaluation of the elastic displacements. In order to facilitate comparison with the analytical method, the lubricant was assumed as Newtonian incompressible and its variations of viscosity with pressure were assumed to follow the well-known Barus law.
In all cases presented, the numerical solution domain was set to −2.5 ≤ X ≤ 1.5 with 17 and 513 discretization points, respectively, at the coarsest and finest levels. The time step was conveniently set equal to the spatial discretization length at the finest grid. The number of time steps was in turn set to have three complete oscillations of the applied contact load. The residual norm of Reynolds equation was kept in the order of 1e−4 or smaller at all time steps. The possible effects of mesh density and of artificial numerical starvation due to fluctuations in contact length were investigated by solving the time-dependent case with the largest load amplitude (29 N) using a numerical domain −4.5 ≤ X ≤ 1.5 with 1025 discrete points at the finest level. The accuracy of the film thickness fluctuations predicted in both settings was then within 1%.