Fatigue Damage of an Asperity in Frictionless Normal Contact with a Rigid Flat

Abstract

Surface fatigue wear widely exists, and it occurs as long as a sufficient number of loading–unloading cycles are applied. Slowing down surface fatigue wear requires understanding the evolution of fatigue damage in the surface. Real surfaces are composed of many asperities; therefore, it is important to study the fatigue damage of a single asperity. A finite element model of an asperity subjected to cyclic elastic–plastic normal loading was developed under frictionless contact condition. The asperity can be either completely or partially unloaded in a loading cycle. For the sake of completeness, both cases were investigated in the present study. The multiaxial Fatemi-Socie fatigue criterion was adopted to evaluate the fatigue damage of the asperity in elastic shakedown state, which was achieved after several loading cycles. For the case of complete unloading, severe fatigue damage was confined in a subsurface ridge starting from the edge of the maximum loaded contact area. The shape and volume of the wear particles were predicted based on a fundamentally valid assumption. For the case of partial unloading, the fatigue damage was much milder. Finally, potential research directions to expand the current study are suggested.

1. Introduction

Surface fatigue wear caused by contact loading is harmful for the reliability and durability of rough surfaces, which are very common in many applications, including false teeth [1], artificial hip joints [2], magnetic storage devices [3], and mechanical components such as gears [4], as seen in Figure 1a. There are myriad asperities on a rough surface [5]; therefore, it is crucial to study the fatigue of one single asperity subjected to contact loading. In practice, a single asperity can be subjected to very complex contact loading such as cyclic combined normal and tangential loading with varying amplitude and phase difference [6]. However, to gain basic physical insights, one may start with investigation of an asperity in normal contact with a rigid flat, as shown in Figure 1b.

Jackson et al. [7] performed a finite element study in a sphere and found that the residual stresses at contact edge increase during unloading and secondary plasticity might occur. Etsion et al. [8] studied the unloading contact behavior of a sphere and obtained expressions describing the relationships between various contact parameters such as contact load, contact area, and interference in one complete loading–unloading cycle. Zhao et al. [9] and Zhao et al. [10] considered the effect of strain hardening on the loading–unloading contact behavior under slip (frictionless) and stick contact conditions, respectively. Jana et al. [11] analyzed the loading–unloading behavior of spheres of elastically and plastically graded material and showed the effects of material gradation parameters on the stress field evolution and the energy dissipation due to plasticity.

The aforementioned studies [7,8,9,10,11] only investigated the mechanical behavior in one loading–unloading cycle. To evaluate the fatigue damage in the asperity, multiple cycles should be applied. Kadin et al. [12] studied the multiple loading–unloading process of an elastic–plastic sphere under frictionless condition. They found that after completion of the first loading cycle, no additional plastic deformation was accumulated in the subsequent cycles reaching elastic shakedown. Zait et al. [13] reported that under full stick contact condition, the accumulation of plastic deformation stopped after a few cycles. Wang et al. [14] employed the Drucker–Prager plasticity model for spheres made of geomaterials, and observed that with increasing cycles, spheres of material with relatively low yield strength experience ratchetting at first, then plastic shakedown, and, eventually, elastic shakedown.

Although useful information about an asperity subjected to cyclic normal loading has been provided [12,13,14], little attention has been paid to asperity fatigue that leads to microcracks and eventually material wear. Xu and Komvopoulos [15,16] explored the fatigue microcrack growth mode in an elastic asperity with initial crack subjected to cyclic adhesive normal contact and sliding contact. However, the position and orientation of the initial crack in asperities, which was unknown a priori, were assumed arbitrarily in their studies. Zhao et al. [17] introduced a damage variable that characterizes the material degradation with increasing loading cycles. The crack initiation site was identified where the damage variable first reached a predefined critical value, and the initiation sites were intensively located slightly below the contact area. Unfortunately, the damage variable based approach can only predict the microcrack initiation position but not orientation, which is important to the crack growth mode, as evidenced elsewhere [15,16]. Meanwhile, existing work [15,16,17] on the contact fatigue in an asperity only considers the plane strain condition, which is only relevant to surfaces with uncommon 2D roughness. Fatigue analysis of a 3D asperity is still missing.

In this study, based on an elastic–plastic finite element model, the multiaxial Fatemi-Socie fatigue criterion was employed to investigate the fatigue microcrack position and potential fatigue wear particle in an asperity subjected to cyclic normal loading. The asperity can either be completely or partially unloaded in a loading cycle. For the sake of completeness, both cases are discussed.

2. Elastic–Plastic Spherical Contact

Figure 1b shows a deformable and frictionless spherical asperity of radius R in contact with a rigid flat on which a normal load P is applied. The normal load P contributes to a circular contact area A with radius a and interference ω. With the normal load P increasing, the asperity deforms elastically until reaching the critical load P_c. Once P rises above P_c, the plastic deformation will occur. The critical contact load P_c, interference ω_c and contact radius a_c under frictionless condition can be calculated as follows [18]:

$$P_{\mathrm{c}}=\frac{{\pi }^3{C}_v^3{Y}^3{R}^2}{6E^2}$$

(1a)

$${\omega }_{\mathrm{c}}=\frac{{\pi }^2{C}_v^2{Y}^{2}R}{4E^2}$$

(1b)

$$a_{\mathrm{c}}=\frac{\pi C_{v}YR}{2E}$$

(1c)

where v, E and Y are the Poisson’s ratio, Young’s modulus, and yield strength of the material, respectively, C_v = 1.30075 + 0.87825v + 0.54373v², and the critical contact area A_c = πa_c². Their corresponding contact parameters P, ω and A can be normalized using these critical values. Therefore, universal dimensionless relationships for P^* vs. ω^* and A^* vs. ω^* are obtained [19]. These relationships for unloading depend on the maximum dimensionless load in a cycle, P^*_max. For subsequent loading cycles with unchanged P^*_max, these relationships are the same as those for the first cycle.

3. Finite Element Model

A finite element (FE) model was established to simulate the asperity contact using the commercial software ANSYS 19.2 developed by Ansys, Inc. based in Canonsburg, PA, USA. The 3D asperity was simplified as a 2D axisymmetric model [20,21,22], which was divided into four separate zones. The mesh density decreases as the distance from the sphere summit increases, because high mesh density is necessary near the sphere summit where a high stress/strain gradient exists. The strip in zone I beneath the sphere surface is to guarantee the full coverage of the loaded contact area. Key dimensions of the FE model are illustrated in Figure 2a and the element sizes in zones I–IV are 2 × 10⁻⁴R, 3.2 × 10⁻³R, 0.0125R and 0.05R, respectively. The FE model consists of 18,657 six-node axisymmetric triangular solid elements (PLANE 183) comprising 37,319 nodes, while the rigid flat and sphere surface were modeled by a target element (TARGE 169) and several contact elements (CONTA 172), respectively.

The material was assumed to be isotropic and elastic–perfectly plastic. The elastic and plastic deformation behaviors were governed by Hooke’s and Prandtl–Reuss laws, respectively. These constitutive relationships were settled using the built-in functionality of ANSYS. The transition from elastic to plastic deformation was determined by the von Mises criterion. The material properties of the asperity were Young’s modulus E = 200 GPa, Poisson’s ratio v = 0.32, and yield strength Y = 210 MPa that are typical for low carbon steels. The dimensionless results are independent of R as long as the theory of continuum mechanics holds [23]; therefore, the sphere radius R was set as 10 mm without losing generality. The contact interface was assumed to be frictionless and the encastré boundary condition was applied at the sphere bottom.

Figure 2b shows the cyclic triangular-wave normal loading history applied on the rigid flat. The normal load in the 2nd and subsequent cycles changed from bP^*_max to P^*_max, where 0 ≤ b < 1. b = 0 means that the normal load decreased to 0 at the end of each loading–unloading cycle, while other values of b correspond to partial unloading cases. The increment for each load step is 0.5 P_c to ensure convergence. A parameter of 1 ≤ P^*_max ≤ 700 was selected for investigation to cover a wide enough loading range. A thorough model validation of the FE model is given in Appendix A.

4. Fatemi-Socie Multiaxial Fatigue Criterion

Material points in the asperity undergo multiaxial stress state when the normal load is applied. Moreover, the principal directions constantly change during the loading process. Under this circumstance, the critical plane-based fatigue criteria are most reliable in predicting the position and orientation of fatigue microcracks [24]. The Fatemi-Socie multiaxial fatigue criterion [25] correlates well with both low and high cycle fatigue of materials for gears [26], bearings [27], etc., and is thus adopted in this paper. According to this criterion, the fatigue microcrack at one material point initiates along the critical plane that undergoes the maximum shear strain range over one loading cycle. If there are several critical planes for one material point, the critical plane is defined as the plane where the fatigue damage is the largest, which is quantified by the Fatemi-Socie (FS) damage parameter (it hereafter will be denoted by FSDP for brevity) [25]:

$$\mathrm{FSDP}=\frac{\Delta {\gamma }_{\mathrm{cp}}}{2}\left(1+\frac{{\sigma }_{\mathrm{cp},\max }}{Y}\right)$$

(2)

where Δγ_cp and σ_cp,max are the shear strain range and maximum normal stress on the critical plane over one cycle, respectively, and Y is the material yield strength. The initiation life (the number of cycles required for the fatigue microcrack initiation) is a monotone decreasing function of FSDP. As can be seen from Equation (2), a positive σ_cp,max increases FSDP. This is because a normal tensile stress can reduce the friction between the crack surfaces and accelerate crack generation.

5. Framework of Determining the Critical Plane and FSDP

To illustrate the framework of determining the critical plane orientation and FSDP for every material point in the asperity, the numerical procedure in Figure 3 is given and the details are explained as follows.

1.

Extract the stabilized stress–strain response in shakedown state at each material point over one cycle from the FE model (see Section 3 for details);

2.

Seek the critical plane at each material point.

To seek the critical plane of a material point requires examining the shear strain range of all planes passing through this point. The plane to be examined is denoted as the x₂ plane in Figure 4a. The orientation of the x₂ plane can be described by two Euler angles, α and β. The coordinate x₀y₀z₀ is the original from which the stress–strain response is obtained from the FE model and the origin locates at the material point’s position. α and β denote the rotation angle of coordinate x₀y₀z₀ around z₀ axis and the rotation angle of the new coordinate x₁y₁z₀ around y₁ axis, respectively. Both α and β change from 0° to 180° with an interval of 1° so that 181² = 32,761 orientations of the x₂ plane are obtained, which are enough to represent arbitrary planes passing through a material point. The strain components in coordinate x₂y₁z₁ (denoted by superscript ’) can be obtained using the following coordinate transformation.

$${\epsilon }^{\prime }=M\epsilon M^{\mathrm{T}},{\sigma }^{\prime }=M\sigma M^{\mathrm{T}}$$

(3a)

$$\epsilon =\left[\begin{array}{ccc}{\epsilon }_{xx}& {\gamma }_{xy}/2& 0\\ {\gamma }_{yx}/2& {\epsilon }_{yy}& 0\\ 0& 0& {\epsilon }_{zz}\end{array}\right],\mathsf{\sigma }=\left[\begin{array}{ccc}{\sigma }_{xx}& {\tau }_{xy}& 0\\ {\tau }_{yx}& {\sigma }_{yy}& 0\\ 0& 0& {\sigma }_{zz}\end{array}\right]$$

(3b)

$$M=\left[\begin{array}{ccc}\cos \beta & 0& -\sin \beta \\ 0& 1& 0\\ \sin \beta & 0& \cos \beta \end{array}\right]\left[\begin{array}{ccc}\cos \alpha & \sin \alpha & 0\\ -\sin \alpha & \cos \alpha & 0\\ 0& 0& 1\end{array}\right]$$

(3c)

where ε and σ are the strain and stress matrices in coordinate x₀y₀z₀, respectively. M is the transformation matrix, and superscript ^T denotes the transpose operation.

Performing coordinate transformation for the strain matrix of a material point at every load step in a cycle, the trajectory of the tip of shear strain vector (γ_xy’, γ_xz’) on x₂ plane can be obtained as shown in Figure 4b. As can be seen, this trajectory is usually not a simple straight line due to the non-proportionality of contact loading. To determine the shear strain range, the Minimum-Circumscribed Circle proposal [28,29,30,31] can be used. By searching among the 32,761 orientations of the x₂ plane, the critical plane having the maximum Δγ’, i.e., Δγ_cp, can be found.

3.

Calculate the FSDP at each material point.

The normal stress history on the critical plane in a cycle can be obtained using Equation (3) once the critical plane is settled. Taking the maximum normal stress on the critical plane, σ_cp,max, and the shear strain range on the critical plane, Δγ_cp, into Equation (2), FSDP is obtained.

6. Results and Discussion

6.1. Shakedown Analysis

As mentioned in Section 4, it is a premise for using the Fatemi-Socie criterion that the material points enter the shakedown state where the stress–strain responses are stabilized in a cycle. The relationships between various contact parameters in the unloading process are only dependent on the maximum dimensionless load in a cycle, P^*_max, as mentioned in Section 2; therefore, the shakedown analysis is only executed for complete unloading cases (b = 0) without losing generality. The stress–strain responses up to eight loading cycles were extracted to check whether shakedown was achieved.

Figure 5a shows the increment of equivalent plastic strain Δε^p_eq in the 1st, 3rd, 5th, and 8th cycles on the surface for P^*_max = 700. Δε^p_eq mainly takes place in the 1st cycle, and the increment gradually decreases to zero after a few cycles so that the material points on the sphere surface reach shakedown state. Different decay rates can be observed for Δε^p_eq. Points at the edge of maximum loaded contact area A_max (x/a_max equal or close to 1) spend most cycles reaching the shakedown state among all surface points. The stress–strain response of the surface point at x/a_max = 1 is shown in Figure 5b for the first eight cycles. The response gradually becomes linear and reversible in the last few cycles, indicating that all material points enter the shakedown state after eight cycles. Figure 5c,d show the same content as in Figure 5a,b for P^*_max = 10. In this case, elastic shakedown occurs in the second cycle, much faster than the case of P^*_max = 700.

A general rule was found that for P^*_max ranging from 1 to 700, the number of loading cycles required to reach elastic shakedown increased with rising P^*_max. Therefore, calculation for the stress–strain response in the first eight cycles can ensure that all material points entered the elastic shakedown state, and fatigue analysis was conducted based on the results of the 8th loading cycle thereafter.

6.2. Complete Unloading Cases (b = 0)

Complete loading–unloading widely exists in gears, bearings, and other mechanical components. In addition, it is a special type of cyclic normal loading. Therefore, complete unloading cases are discussed firstly, while the results of partial unloading cases (b ≠ 0) will be illustrated later in Section 6.3.

Using the numerical procedure illustrated in Figure 3, Δγ_cp and FSDP can be obtained for all material points in the asperity. Figure 6a–c show the contours of Δγ_cp/(Δγ_cp)^max at the sphere tip for three representative P^*_max ranging from 1 to 700, where d = R − y denotes the depth from the sphere summit in the axial direction. The area of high Δγ_cp surrounded by red lines in Figure 6a–c expands in the subsurface as P^*_max increases. Figure 6d–f show the contours of FSDP/(FSDP)^max at the sphere tip under the same loading cases as in Figure 6a–c. A high FSDP area can be observed in the subsurface within 0.5b away from the sphere summit. The maximum FSDP among all the material points, (FSDP)^max, occurs near the edge of A_max, where x/a_max = 1.

To better illustrate the 2D contours in Figure 6, 3D contours for P^*_max = 10 and P^*_max = 700 are shown as examples in Figure 7. For both cases, a ridge-like region in the Δγ_cp contour (see Figure 7a,b) can be observed that starts from the edge of A_max, extends to the subsurface, and intersects the axis of symmetry at a certain depth from the surface. On the other hand, the FSDP contours for P^*_max = 10 and P^*_max = 700 look different, as shown in Figure 7c,d. Although high FSDP occurs in the ridge-like region in the subsurface of the asperity, the number and depth of the ridges varies for the two cases. When P^*_max = 10, there is only one single ridge in the FSDP contour. When P^*_max = 700, two ridges exist. The FSDP values for points on these ridges are significantly higher; therefore, fatigue microcracks are much more prone to initiate here.

Moreover, the FSDP contours for P^*_max = 10 and P^*_max = 700 are representative of two typical FSDP distribution patterns for P^*_max < 30 and P^*_max ≥ 30, respectively. Figure 8 shows FSDP at the ridge ends, which are denoted by points 1, 2 and 3 (see Figure 7c,d), as functions of P^*_max. Point 3 only exists for P^*_max ≥ 30. Moreover, FSDP at point 3 is significantly lower than that at point 2, as are other points on the lower ridge compared to those on the upper one (see Figure 7d). Hence, for P^*_max ≥ 30, the upper ridge is much riskier in terms of fatigue microcrack initiation. Thus, the following discussion mainly focuses on the single ridge for P^*_max < 30 and the upper one for P^*_max ≥ 30.

To illustrate the joint effect of the shear strain range and the maximum normal stress on FSDP, the ridge in the Δγ_cp contour and the (upper) ridge in the FSDP contour are presented in the same picture for 1 ≤ P^*_max ≤ 700 (see Figure 9). The FSDP ridge is located closer to the surface for the whole load range. With increasing P^*_max, the gap between the two ridges becomes wider, indicating that the effect of σ_cp,max is more significant for higher P^*_max. This makes sense because a higher P^*_max is associated with stronger compression on the surface. This compression “pushes” the subsurface material away from the sphere tip and results in tensile normal stresses on the critical planes of subsurface points close to the surface. Consequently, the FSDP ridge is located closer to the surface.

An interesting phenomenon is that both (FSDP)^max and (Δγ_cp)^max occur at the edge of A_max, where x/a_max = 1 for P^*_max > 1. To explain this, the residual radius of curvature of the asperity after complete unloading was investigated. As shown in Figure 10a, the portions of the surface within A_max (x/a_max ≤ 1) and outside (x/a_max > 1) that are represented by solid and dash-dotted curves, respectively, can be fitted by two arcs with different radii of curvature. Figure 10b shows that the increase in ρ₁ becomes more salient with increasing P^*_max. This is because of the residual plastic deformation beneath the loaded contact area. Meanwhile, ρ₂/R is equal to 1 for x/a_max > 1. In other words, this portion of surface after unloading coincides with the original sphere surface. The difference between ρ₁ and ρ₂ leads to a sharp change of radius of curvature at the edge of A_max, where x/a_max = 1. Therefore, the deformation is concentrated in the region close to the edge of A_max where (FSDP)^max occurs.

After undergoing a certain number of loading–unloading cycles, the fatigue microcrack will first initiate at the point with (FSDP)^max along its critical plane. As the cycle number further increases, more microcracks will subsequently initiate at points in the descending order of their FSDP. These microcracks will propagate and coalesce to long cracks, as observed in extensive experiments [32,33]. Unfortunately, there does not seem to be a unified theory that can accurately describe the trajectory of both crack initiation and propagation process [27]. It is worth noting that the points on the (upper) ridge in the FSDP contour have much higher FSDP values than those on both sides. Therefore, it is reasonable to assume that the crack is formed along this ridge. This long crack results in a fatigue wear particle detached from the bulk of the asperity. Figure 11a shows the (upper) ridge in the FSDP contour for 1 ≤ P^*_max ≤ 700. With increasing P^*_max, the ridge moves deeper into the subsurface, as would be expected. Figure 11b shows the dimensionless volume of the wear particle normalized by the total volume of the asperity under different P^*_max. The detached volume accounts for no more than 0.1% of the total volume of asperity, indicating that the fatigue wear in an asperity subjected to cyclic normal loading is very mild, even if the normal load is much higher than its critical value P_c.

6.3. Partial Unloading Cases (b ≠ 0)

In most engineering problems, asperities may undergo complex loading histories. Although complete unloading cases where b = 0 have been investigated in detail above, partial unloading cases are more general and need to be investigated.

Figure 12a–d show the 2D contours of FSDP when b = 0.1, 0.25, 0.5 and 0.75, respectively for a typical case of P^*_max = 700. As in the complete unloading cases, (FSDP)^max occurs at the edge of A_max. However, the high FSDP area is restricted beneath the edge of A_max, while a low FSDP area occurs in the subsurface within A_max. Higher values of b accompany a larger area with low FSDP. This is because the material points within A_max undergo lower Δγ_cp than those in the case of b = 0. According to Equation (2), this brings about lower FSDP.

Figure 13 shows the (FSDP)^max as functions of P^*_max under different values of b. The positive correlation between (FSDP)^max and P^*_max is a similarity between partial and complete unloading cases. Meanwhile, the value of b has a dominant effect on (FSDP)^max under various values of P^*_max, ranging from 1 to 700. Higher values of b not only restrict the high FSDP area, but also decrease (FSDP)^max.

Further work can be done to supplement and improve the present work. Firstly, because deviations from reality are inevitable due to friction and time-varying temperature, experiments are strongly recommended as the solid validation of the theoretical analysis results. Other complex loading histories can be applied in addition to triangular-wave normal loadings. The constitution law of the material should be carefully selected to better simulate the practical surface wear. Moreover, the fatigue crack propagation can be investigated to predict the volume of fatigue wear particle accurately. Finally, the relationship between the fatigue in a single asperity with the overall fatigue damage of a rough surface remains to be studied.

7. Conclusions

In this study, the fatigue damage of a frictionless asperity subjected to cyclic normal loading was evaluated using the multiaxial Fatemi-Socie fatigue criterion. The shakedown analysis based on an elastic–plastic FE model showed that all material points in the asperity entered the elastic shakedown state after a few loading cycles. The cyclic loading applied to the asperity can be divided into two categories, i.e., the asperity is completely unloaded and partially unloaded in a loading cycle. For the case of complete unloading, ridge-like regions were observed from both the Δγ_cp and FSDP contours at the sphere tip. The maximum fatigue damage, (FSDP)^max occurred at the edge of the maximum contact area A_max. This resulted from the sharp change of radius of curvature at the edge of A_max after unloading. Two typical patterns for FSDP contours were found for P^*_max < 30 and P^*_max ≥ 30. Fatigue microcracks were prone to occurring along the single ridge in the FSDP contour for P^*_max < 30 and the upper one for P^*_max ≥ 30. The (upper) ridge in the FSDP contour was assumed as the long crack that eventually resulted in a wear particle detaching from the asperity. The wear particle only occupied a tiny portion of the whole asperity, indicating that contact fatigue wear is very mild. For the case of partial unloading, the region of high FSDP was only restricted below the edge of A_max. Moreover, the ridge-like regions of high FSDP were not observed in this case, and higher values of b led to lower (FSDP)^max, as would be expected. Finally, directions of future work have been suggested.