On the tribo-dynamics of lubricated point contacts: modeling and experiment

This study proposes a modeling methodology for the description of tribo-dynamic behavior of generic point contacts, considering surface defect/failure type of excitation. A twin-disk contact set-up is developed to experimentally determine the acceleration induced by lab-controlled groove surface feature present on one of the mating surfaces. The periodic passage of this artificial surface defect through lubricated contact zone leads to contact force fluctuation, and consequently normal direction dynamic relative displacement and relative velocity of the contact bodies. These dynamic responses impact the lubrication film thickness and flow, altering the tribological performance that governs the contact force. Implementing this interaction loop between dynamics and tribology, the computational simulations are shown to be able to yield fairly accurate predictions in comparison to the measured acceleration signals.


Introduction
Dynamic responses, such as noise and vibration, are common in applications of high-speed rolling mechanical elements (gears and bearings for instance), whose cyclic contact actions are deemed to be the main sources of excitations [1][2][3][4][5][6][7]. From an interdisciplinary point of view, dynamic vibration has been shown to have large influence on tribological performance occurring at the interface of contacting components, including fluctuating normal force and lubricant viscosity [8,9]. Although the physical interactions between tribology and dynamics are evident and substantial [9], a majority of works in the literature focused on only one way of impact. In the early 80s of last century, Wang and Cheng [10] numerically evaluated the thermal elastohydrodynamic lubrication (EHL) solutions of spur gear contacts subject to dynamic tooth force. Considering both normal and tangential vibrations, Yang et al. [11] investigated the responses in terms of lubrication film thickness, pressure, and temperature of a line contact. Both studies assumed perfectly smooth surfaces. Including surface roughness effects, Mohammadpour et al. [12] examined friction and efficiency of planetary gear systems through mixed thermal EHL analyses under dynamic condition. Joining a lumped parameter vibratory model to EHL governing equations, Li and Kahraman [8] demonstrated the unique transient gear EHL behavior under nonlinear dynamic load and measured surface roughness condition. All of these works showed strong dependence of tribological responses on dynamics. However, the reaction of dynamics to tribology was excluded. Along the opposite direction, tribology was introduced into the field of gear dynamics in a relatively shallow way, where the complex physics at the lubricated contact interface was simply represented by an empirical friction coefficient. Using such an approach, friction induced gear vibration was extensively studied, while leaving the reaction of friction to dynamics untouched [13][14][15][16]. Expanding to other aspects of tribology, Chen et al. [17] investigated the effect of lubricantinduced backlash reduction on gear dynamics, which was shown to be evident for high-precision gearing applications.
The first work that attempted to model the mutual relationship between gear dynamics and gear tribology was presented by Li and Kahraman [9]. Through an iterative numerical method, the equations of motion and the gear EHL governing equations were coupled together. The former yielded alternating surface velocities and contact force, which were used as inputs for the latter to determine friction forces and viscous shearing induced damping, both taking place along the tangential direction of contact surfaces. These tribological responses were then fed back into the dynamic analysis to update the tooth force and velocities till a convergence was reached. Employing this iterative practice, power loss [18], flash temperature [19], and contact fatigue [20] under tribodynamic condition were examined for cylindrical gear contacts. Considering hypoid gear applications, where the complex geometry imposes overwhelming challenges, such tribo-dynamic coupling was realized by implementing an interpolation based parametric surface approach [21].
Although gear tribology and gear dynamics have been claimed to be physically bridged [9,[18][19][20][21], they were bridged partially. These studies considered the tribological responses only in the direction that is tangent to the contact surfaces. The normal direction response [22] was overlooked. In a recent investigation, Li and Kolivand [23] discussed the interaction between dynamics and tribology in the normal direction of a line contact. It was shown surface irregularities, such as pit/defect and roughness, alter contact pressures and thus normal contact force. The resultant imbalance between the contact force and applied normal load initiates normal direction movement of contact bodies, leading to alternating displacements, velocities, and accelerations. Along the normal direction, dynamic relative displacement directly contributes to the thickness of lubrication film; and dynamic relative velocity is well linked to the squeeze term of Reynolds equation that governs fluid flow. Through this way, dynamics get back into action, and speaks loudly for its role in tribological responses. Utilizing the tribo-dynamic approach proposed, Li and Kolivand [23] quantified and demonstrated the relationship between fatigue pit dimension and acceleration response, providing a practical and effective way for the detection of fatigue pits of interested sizes. In this study, the assumption of a line contact, however, limits its applications. For instance, needle roller bearings or spur gears commonly have axial direction crown implemented to avoid any edge loading condition, where a line contact is no longer applicable. Experimental assessment of the model was also missing.
As discussed above, the physics-based tribo-dynamic description of lubricated rolling contact is yet not complete. Considering generic point contacts, a three degree-of-freedom (DOF) dynamic system is constructed to model the contact of an experimental set-up. In the experiment, a cylindrical disk pair is used to form the contact. A circular crown is applied along the axial direction on one of the disks to produce an elliptical point contact. On the other disk, no crown is applied, while a groove, which is aligned along the axial direction, is fabricated on the surface. This groove introduces cyclic surface topography variation within the contact zone, serving as lab-controlled excitation. Acceleration signals are captured through an accelerometer for different groove sizes under two loading levels, and compared to model predictions to demonstrate its effectiveness. In view of the limited liquid-solid interface slippage [24], no-slip boundary condition is adopted in the modeling process. This study addresses the tribo-dynamic behavior along the normal direction of contact, considering more practical point contact condition.

Experimental set-up and modeling methodology
To allow assessment of the proposed modeling methodology for the tribo-dynamic description of lubricated point contacts, this study employs and simulates the experimental set-up as shown in Fig. 1, where the contact is formed by two cylindrical disks of the same radius of 28.575 mm. Disk 1, which sits on the left, is crowned along the face width direction with the radius of 76.2 mm. Its counterpart, Disk 2, has no crown implemented. However, a groove is cut along its axial direction as shown in Fig. 2a to produce periodic surface topography excitation during operation. Figure 2b displays the groove cross section shape. For the investigation of groove size effect, three different sizes as defined in Table 1 are adopted in this work. Both disks of the contact pair are highly polished (after the groove is fabricated on Disk 2), such that their surface roughness, whose RMS amplitudes are 0.06 lm (outside the groove area for Disk 2), can be treated to be negligible in relative to the groove surface feature dimension and excluded in the computational process for faster convergence rate. Disks 1 and 2 are mounted onto their respective shafts, each of which is supported by four high-precision ball bearings accommodated within the bearing housing. As illustrated in Fig. 3, the shafts extend all the way through the housings and are belt driven by two independent motors. The motor and bearing housing associated with Disk-assembly 1 are firmly attached to the test platform, while those related to Disk-assembly 2 are allowed to move freely along the horizontal rails as displayed in Fig. 1. To provide normal force, F 0 , a hydraulic cylinder, which sits on the right end (not shown in Fig. 1), is utilized to push the bearing housing through a loading plate (that is bolted to the former) to the left, producing contact between the disk specimens. For the control of F 0 , a load cell (1210ACK-1K-B) is implemented between the hydraulic cylinder and the loading plate. A singleaxis ICP accelerometer (353B18) is fastened to the movable bearing housing as shown in Fig. 1 to measure the vibration signal along the horizontal direction, which is normal to the contact surfaces of the disks. A portable signal analyzer (Quattro DSPcentric signal processing hardware equipped with Sig-nalCalc ACE dynamic signal analysis package) is employed for the data acquisition. This study considers two loading levels of F 0 ¼ 900 and 1700 N, which result in Hertzian contact pressures of p h ¼ 1:37 and 1.69 GPa, respectively. The rolling velocity of the contact pair is set at 6 m/s that corresponds to the rotation frequency of 33.4 Hz. Sliding velocity is excluded to avoid large frictional heat and any scuffing related failure in the presence of significant groove surface feature [19]. It is noted the rotational speed  controllability under steady state condition of the machine is AE2 RPM. As such, the speed control related surface velocity error is AE0:01 m/s for the disk size employed, which is believed to be negligibly small. For lubrication, a typical turbine fluid [24], whose temperature is set at 40 C, is supplied into the mesh via a jet.
To describe the tribo-dynamic behavior of the above contact set-up, a three-DOF lumped parameter model is constructed in Fig. 4, where Disk-assemblies 1 and 2, and the motor-housing-assembly associated with Disk 2 are represented by masses of m 1 , m 2 , and m 3 , respectively. The shaft-bearing supports of the disk-assemblies are modeled as spring-damper elements, denoted by k 1 -c 1 and k 2 -c 2 . The stiffness, k 3 , that represents the compliance of the loading plate is set in between m 3 and the right wall of the test platform. The corresponding damping component, c 3 , is also put in place. The magnitudes of these mass and stiffness elements are listed in Table 2, where the adopted damping ratios are provided as well.

Governing equations
When the surface groove passes through the contact zone, it introduces fluctuation in lubrication film thickness, and consequently oscillation in contact pressure and contact force, F c . The imbalance between F c and the applied normal force, F 0 , results in the movements of the masses, whose displacements along the normal direction of contact are denoted as x d 1 , x d 2 , and x d 3 in Fig. 4. The equations of motion are written as where c m is related to rolling contact damping generated by mechanisms other than EHL film [25]. The associated damping ratio is estimated to be 0.002. The damping effect caused by fluid film squeezing is included inexplicitly below. In Eq. (1), the contact force, F c , is time-varying according to the transient pressure distribution, p, within the lubricated contact zone as where x points into the direction of rolling, and y is directed along the axial direction of the contact disks. The EHL computational domain is circumscribed by the start and end coordinates, x s and x e in the x direction, and y s and y e along the y direction. Additionally, t represents time. The contact pressure that dictates F c through Eq. (2) is well coupled with the lubrication film thickness, h, via the two-dimensional Reynolds equation of [24] o ox / x qh 3 12g where u r is the average of the surface tangential velocities, u 1 and u 2 as depicted in Fig. 4, and is commonly referred as the rolling velocity; and q and g are lubricant density and viscosity, which are dependent on pressure, respectively, as [24] q with q 0 and g 0 representing the density and viscosity under ambient pressure, and a 1 and a 2 denoting the pressure-viscosity coefficients. Together with the other lubricant property related constants involved in Eqs. (4) and (5), their values are referred to Ref. [24]. Lastly, in Eq. (3), / x and / y are flow coefficients, which can be utilized to include any lubricant non-Newtonian shear-thinning behavior. In this study, both flow coefficients take unity value under the pure rolling condition. By means of Eqs. (2) and (3), the time-varying contact force, which introduces perturbation into the equation of motion, is linked to the lubrication film thickness. This key parameter in the field of tribology is commonly composed of surface curvature gap, g 0 , elastic deformation, V, surface topography (such as surface roughness, defects, fatigue pits, or artificial/controlled surface features) heights of surfaces 1 and 2, denoted as S 1 and S 2 , respectively, and constant rigid body approach, h 0 [26][27][28]. In this study, one additional term is devised to describe the normal direction interaction between dynamics and tribology as where ðx d 1 À x d 2 Þ is the dynamic relative displacement of the contact bodies. With such an addition, the squeezing and compressibility of lubrication film, elastic deformation of contact, and supporting shaftbearing deformation are fully incorporated in the compliance of contact. Revisiting Eq. (3), the squeeze term is written as Utilizing Eq. (6), the time derivative of film thickness is decomposed into Therefore, in addition to ðx d 1 À x d 2 Þ in Eq. (6), the dynamic relative velocity, ð _ x d 1 À _ x d 2 Þ, in Eq. (8), contributes directly to the establishment of EHL lubrication film and consequently F c , engaging dynamics and tribology in the normal direction of contact by providing damping mechanism in an inexplicit way.
The Reynolds equation above governs the fluid flow within the EHL conjunction. When surface irregularities are significant or insufficient lubrication film is encountered, local asperity contacts take place. Assuming infinitesimal and constant film thickness within these local boundary lubrication spots, the following equation describes the contact behavior [26,28] oh ox ¼ 0 ð9Þ To couple the above dynamic [Eq. (1)] and tribological [Eqs. (3)(4)(5)(6)(7)(8)(9)] governing equations, the numerical procedure that has been shown to be effective in a line contact problem [23] is implemented and illustrated in Fig. 5. At the start, EHL simulation under quasi-static condition [28] is carried out to find pressure and film thickness information, which values serve as initial guesses of the EHL solutions of Eqs. (3)(4)(5)(6)(7)(8)(9). Setting the disk dynamic displacements and velocities at zero, i.e., x d i ¼ 0 and _ x d i ¼ 0 (i = 1, 2), the EHL simulation is performed to determine the contact pressure distribution. Via Eq. (2), the contact force is obtained and passed to the dynamic simulation to determine the target tribo-dynamic response at this time instant, i.e., the instantaneous bearing housing (m 3 ) acceleration. Preparing for the next time instant, relative displacement and velocity of the disks are updated and routed back to the EHL stage. This time-stepping computational process is executed for a large number of rotation cycles (100 disk rotations) to allow steady state dynamic signal to be reached (According to tryouts, steady state signal can be reached after as many as 50 rotation cycles.). On an i7-PC, the involved computational time for each simulation is on the order of 200 h. Only the last 2 cycles of signal are used to quantify the acceleration response.

Simulations and discussion
To demonstrate the capability of the proposed tribodynamic modeling methodology, the predicted dynamic responses of the contact pair outlined in Sect. 2 are analyzed and compared to the measurements under different loads (900 N and 1700 N) and surface groove excitations (Table 1). For these six simulations, the polished surface roughness heights are neglected in the film thickness equation, since they are negligibly small in comparison to the groove topography. In view of the pure rolling condition adopted, the non-Newtonian shear thinning effect does not apply, and / x ¼ 1 and / y ¼ 1 in Eq. (3).
According to Sect. 3, the dynamic relative displacement (x d 1 À x d 2 ) and dynamic relative velocity of the disk pair, respectively, impacts the lubrication film thickness and the film squeezing behavior, contributing to the tribological performance. In Fig. 6, the phase plots of these dynamic parameters are constructed for two consecutive disk rotation cycles. It is shown that both the displacement and velocity amplitudes increase with the groove size under the high loading condition of 1700 N (left column). When the load is reduced to 900 N (right column), similar trend is found, while such variations are relatively small, especially for the amplitude of ( _ . Comparing the right column with the left one, the reduction in normal load is seen to allow largely elevated dynamic displacement amplitudes under all groove size conditions. Such load effect on the velocity amplitude is found to be non-tangible. Figure 7 depicts the impact of these dynamic responses on the average film thickness within the nominal Hertzian contact zone under the high load scenario. The two high peaks in the plots are due to the passages of the deep groove valley though the contact zone during the two disk rotation cycles. When the Fig. 5 Numerical procedure employed for the solution of tribodynamic response groove feature is out of the Hertzian zone, the film thickness variation can be better observed in Fig. 8, where zoomed-in images are displayed. In general, more positive relative displacement, i.e., the disk surfaces get closer to each other, leads to smaller film thickness. Regarding the relative velocity, its influence on film thickness is relatively small and difficult to isolate owing to its dependence on the relative displacement, which varies with it at the same time. When the lubrication film thickness fluctuates, the contact force that is tightly related to the tribological performance (EHL pressure and asperity contact pressure if any) oscillates as illustrated in Fig. 9. The higher relative displacement is seen to result in higher F c , which is in line with the average film thickness behavior (Smaller film thickness corresponds to higher contact pressure, and therefore higher contact force). In each of the plots, the two large drops in F c correspond to the film thickness jumps in Fig. 7 as the groove valley travels through the contact zone. When the groove size gets larger, more significant perturbation is introduced in the contact force, which is responsible for the more evident dynamic responses in Fig. 6. Revisiting Fig. 8, the film thickness is seen to continue its dip after the groove exists the contact zone. This is due to the positive relative velocity that drives the further approach of the disk surfaces, accompanied by the increase of the contact force. The Yet, one noticeable difference in the variation of F c can be observed between the two loading cases. In Fig. 9, a high frequency oscillation is seen in the contact force as it increases with (x d 1 À x d 2 ). Such behavior only appears at the high end of (x d 1 À x d 2 ) in Fig. 11. The reason could be attributed to the film thickness difference between the two cases. Using the large groove excitation as an example in Fig. 12, the average film thickness under the high load condition  ranges from 0.42 to 0.46 lm (when the groove is out of the contact zone), which is below the film thickness range of 0.496-0.561 lm under the 900 N normal force. When the film thickness is thinner, the contact force is more sensitive to film thickness fluctuation, resulting in the high frequency oscillation of F c in Fig. 9. As (x d 1 À x d 2 ) increases, the approach of the contact surfaces leads to the film thickness decrease. Therefore, at the high end of (x d 1 À x d 2 ) in Fig. 11, where the film thickness becomes sufficiently small, such contact force oscillation becomes evident. In Fig. 9, the ranges of F c is recorded to be on the order of 200 N, 250 N, and 300 N for the small, medium, and  large grooves, respectively. Correspondingly, these respective ranges are on the order of 350 N, 400 N, and 450 N when F 0 is reduced from 1700 to 900 N in Fig. 11. The larger peak-to-peak values of F c under the low load condition correlates well with the muchelevated displacement amplitudes observed in Fig. 6, when F 0 is reduced. The experimental measurements as described in Sect. 2 are carried out. And the acceleration signals are compared with the predictions in Fig. 13 under the high load condition. Three resonance peaks are observed at 142 Hz, 1370 Hz, and 1500 Hz in the black prediction curve, which agree very well with those of the measured signals. Although the two high frequency peaks of the measurement performed with the small size groove are not very distinctly separated from each other, they are easily distinguishable when larger groove sizes are employed in Fig. 13b and c. In view of the acceleration magnitude, the black curve predictions and the measurements are in good agreement at the low frequency resonance of 142 Hz. In the vicinity of the two high frequency peaks, the differences between the predictions and the measurements are somewhat noticeable. The measured signals clearly have wider bands, yet are still reasonably close to the predictions in view of the noises associated with the rotations of ball bearings, belt drives and electric motors. Such dynamic signals can be transmitted to the bearing housing and recorded by the accelerometer, however, are not considered in the modeling process. In Fig. 14, where the smaller load of 900 N is applied, the measured acceleration amplitudes at the low frequency resonance are seen to be increased by 122%, 97%, and 92% for the small, medium, and large grooves, respectively, in comparison with the 1700 N load case. Such significant raises in vibration are in line with the large amplitude jumps   Fig. 6, and larger contact force peak-to-peak values recorded by comparing Fig. 11 with Fig. 9. As for the two high frequency resonances, the changes in acceleration amplitude are limited when the load is reduced. Overall, the model predictions are seen to agree reasonably well with the measurements.

Conclusion
In this study, a modeling methodology for the tribodynamic description of a point contact is proposed, focusing on surface defect/failure induced vibration. Modeling an experimental twin-disk contact set-up, the predictions are compared with the measured acceleration signal to show reasonably good agreements. To introduce surface defect feature that can be easily fabricated on disk surfaces, lab-controlled grooves of different sizes are employed. In the model, a three-DOF dynamic system is coupled with the EHL governing equations that dictate the tribological behavior taking place at the lubricated point contact interface. The contact force fluctuation due to the cyclic passage of the surface groove through the EHL zone is used to evaluate the dynamic relative displacement and velocity of the disk pair, which impact the lubrication film thickness and flow, and therefore the contact force. Because the surface defect induced contact force fluctuation excites the dynamic behavior, an interaction loop between tribology and dynamics is established. Through the computational and