Efficient Sub-Modeling for Adhesive Wear in Elastic–Plastic Spherical Contacts

Abstract

This paper presents a novel approach for simulating adhesive wear in elastic–plastic spherical contacts using an improved finite element sub-model. Initially, a global model with a coarse mesh identifies the potential wear region under combined normal loading and tangential displacement. Subsequently, a refined mesh sub-model simulates the crack initiation and propagation until the formation of a wear particle. This refined sub-model efficiently handles a wide range of spherical radii and normal loads. An expression is derived relating the dimensionless wear volume and wear rate to the dimensionless normal load, revealing the limited effect of the sphere radius on the wear rate. The effect of the mechanical properties on the wear particle morphology is also analyzed.

1. Introduction

Wear is the key factor causing material loss and mechanical failure between contacting surfaces, which typically occurs in various tribo-components. There are mainly five types of wear in tribology: abrasive, erosive, corrosive, fatigue and adhesive [1]. Among them, adhesive wear is the least avoidable between contacting surfaces that are subjected to strong adhesive bonds, in which the material removal occurs at the contacting asperities of the mating surfaces [2]. Although very common, an accurate prediction of adhesive wear is still a challenging question in tribology [3]. Hence, developing an effective model to reveal the mechanism of adhesive wear is of great significance.

Archard [4] replaced Holm’s concept [5] for the removal of atoms with the removal of spherical wear particles of the same radius as their contact area to describe the adhesive wear process. This alternative approach was supported by experimental evidence [6], showing that the material removal is not at an atomic level.

Archard [4] proposed a model for the adhesive wear, based on a linear relation to calculate the wear volume, which was expressed as $W=KPL/H$ where W is the wear volume, P is the normal load, L is the sliding distance, H is the hardness and K is the wear coefficient.

Due to the simplicity and practicality of Archard’s model in engineering applications, it became a well-accepted wear model in the following decades [7]. However, it has inherent limitations, one of which is the lack of a physical understanding of the mechanism of the wear particle formation. Therefore, Archard’s model can only give a reliable prediction for cases where the wear coefficient $K$ was already found experimentally, implying that several wear experiments imitating the correct tribological system are required to obtain wear predictions.

Indeed, a lot of work has been done to perform the integrated analysis of the fracture behavior to enable adhesive wear modeling [7]. The pioneering study performed by Hills and Ashelby [8] introduced the classical fracture mechanics to explore the wear behavior between contacting surfaces. In Ref. [8], the crack was initially generated in a pre-selected location. Subsequently, Suh and coworkers [9,10] developed the delamination theory of wear, in which they assumed the crack nucleation and subsurface crack propagation. Here also, the crack location was predetermined. Although the wear particle formation process can be investigated by a predetermined initial crack location, it limits the crack to propagate from this predetermined location, which may cause an unreliable prediction. Hence, for proper wear modeling, a model that can simulate the crack initiation is required.

The finite element method (FEM) is an efficient method to determine the crack initiation based on the damage mechanism concept. By introducing the Johnson–Cook strain-based ductile fracture criterion [11], Wu and Shi [12] were able to simulate a 2D crack initiation and propagation using a finite element model for a cylindrical contact. According to [12], the crack initiation was at the contact surface, which is in contradiction to the assumption made in the previously mentioned studies [8,9,10].

To simulate the adhesive wear of 3D spherical contacts an FE model, Zhang and Etsion [13] used the JC criterion, similar to [12], for the crack initiation. However, a different approach was used for the crack propagation in [13]. Specifically, a fracture energy criterion [14] was used, where the fracture energy, $G_{\mathrm{f}}$, described the energy needed to open a unit area of the crack [14]. By incorporating this physical approach into the adhesive wear simulation, the crack propagation was modeled as a crack opening process rather than the same way as its generation.

In Ref. [13], the static friction and the adhesive wear of a spherical contact under combined normal and tangential loading were presented for both elastic and elastic–plastic regimes for the normal loading of a sphere with a 10 mm radius. In Zhang’s study [13], an extremely refined mesh was required to obtain an accurate wear prediction, which was very time consuming, limiting the study to one material and two normal loads. However, some insights regarding the formation of the adhesive wear particle were obtained in both [12,13]. The research showed that the crack initiation started at the trailing edge of the contact surface and then, another crack was initiated at the leading edge. As the tangential displacement increased, the cracks propagated and eventually connected with each other, resulting in the formation of a wear particle. It can be inferred that investigating the adhesive wear using predetermined cracks may not provide accurate results.

In order to improve the computing efficiency of the FE wear model [13], Zhang and Etsion [15] adopted a sub-model technique which enabled them to efficiently calculate the wear volume and wear rate for different normal preloads. While this sub-model was appropriate to characterize the wear particle morphology, the model was not able to present the realistic mechanical response during the wear process, which is required to better analyze the mechanism of the wear particle formation, tangential stiffness and friction behavior.

In the present study, an improved sub-model technique is suggested which innates the advantages of both the models previously proposed by Zhang and Etsion [13,15]. Thus, the present improved sub-model is able to efficiently predict the mechanical response in a spherical contact and the wear particle morphology for different radii, normal preloads and material properties.

2. Theoretical Background

Figure 1, taken from [13], presents a schematic representation of the contact problem of a rigid flat and a homogeneous deformable sphere under combined normal and tangential loading. The contact condition between the rigid flat and the outer surface of the sphere is a full stick. The loading process consists of two stages. At the beginning of the simulation, a normal load $P$ is applied to the rigid flat, resulting in a vertical displacement of the rigid flat, also termed as the interference and denoted by ${\mathsf{\omega }}_0$. Due to the deformation of the sphere tip, a circular contact area $\mathrm{A}$ of diameter ${\mathrm{d}}_0$ is obtained. The characteristics of the contact area, friction coefficient and corresponding critical load, $L_{\mathrm{c}}$, were studied previously by Brizmer et al. [16,17]. The normal load is usually normalized by the critical normal load at the yield inception under the full stick contact condition, $P^{*}=P/L_{\mathrm{c}}$, given in Ref. [16].

$$L_{\mathrm{c}}=\stackrel{-}{L_{\mathrm{c}}}\frac{{\mathsf{\pi }}^3\mathrm{Y}}{6}{C}_{\mathsf{\nu }}^3{\left(R\left(1-{\mathsf{\nu }}^2\right)\frac{Y}{E}\right)}^2$$

(1)

where $\stackrel{-}{L_{\mathrm{c}}}=\left(8.88\mathsf{\nu }-10.13\left({\mathsf{\nu }}^2+0.089\right)\right)$ and $C_{\mathrm{v}}=1.234+1.256\mathsf{\nu }$, and $R,E,Y$ and $\mathsf{\nu }$ are the radius of the sphere, Young’s modulus, the yield strength and Poisson’s ratio, respectively. With this normalization, the contact is elastic–plastic for $P^{*}>1$. For $P^{*}<1$, the contact is elastic.

Following the normal preload, a tangential displacement ${\mathrm{u}}_{\mathrm{x}}$ is applied to the rigid flat until a formation of the wear particle is achieved.

Similar to Ref. [13], the adhesive wear propagation is studied using the FEM. To model the fracture process, two fracture criteria are used. First, the Johnson–Cook (JC) criterion [11] is used for the damage initiation of the first crack. According to the JC criterion, the crack occurs when the strain at a material point reaches the equivalent plastic strain ${\mathsf{\epsilon }}^{\mathrm{f}}$ defined in [11] as the following equation.

$${\mathsf{\epsilon }}^{\mathrm{f}}=\left(D_1+D_2{\mathrm{e}}^{D_3{\mathsf{\sigma }}^{\mathrm{*}}}\right)\left(1+D_4\mathrm{l}\mathrm{n}{\dot{\mathsf{\epsilon }}}^{\mathrm{*}}\right)\left(1+D_5{T}^{\mathrm{*}}\right)$$

(2)

As shown in Equation (2), the plastic strain consists of three parts, which describes the effect of the stress, the strain rate and the temperature, respectively. For a quasi-static and isothermal analysis, the effects of the strain rate and temperature are absent $(D_4=D_5=0)$. Hence, Equation (2) can be reduced to the following.

$${\mathsf{\epsilon }}^{\mathrm{f}}=D_1+D_2{\mathrm{e}}^{D_3{\mathsf{\sigma }}^{\mathrm{*}}}$$

(3)

The equivalent plastic strain ${\mathsf{\epsilon }}^{\mathrm{f}}$ in Equation (3) depends on the stress triaxiality ${\mathsf{\sigma }}^{\mathrm{*}}$, which is given in the following form.

$${\mathsf{\sigma }}^{\mathrm{*}}={\mathsf{\sigma }}_{\mathrm{m}}/\overline{\mathsf{\sigma }}$$

(4)

where ${\mathsf{\sigma }}_{\mathrm{m}}$ is the hydrostatic stress and $\overline{\mathsf{\sigma }}$ is the von Mises equivalent stress.

Equation (3) and first bracket of Equation (2) follow the form presented by Hancock and Mackenzie [18]. The coefficient $D_1$ is the crack nucleation strain for the materials that allow for a considerable plastic flow prior to the crack nucleation. $D_1$ will be zero otherwise. $D_2$ is a material constant and $D_3$ expresses how the equivalent plastic strain is affected by the stress triaxiality ${\mathsf{\sigma }}^{*}$. These three parameters can be obtained from the tensile tests [19,20].

The fracture energy ${\mathrm{G}}_{\mathrm{f}}$ [14] for the damage evolution is used to enable the crack propagation. The opening of a crack requires the expenditure of energy to overcome the associated stress. Therefore, the crack opening process leads to the absorption of energy. The amount of energy needed per unit of the crack area in opening the crack from zero to a given width $\mathrm{c}$ is calculated using the following equation [14].

$${\mathrm{G}}_{\mathrm{f}}={\int }_0^{\mathrm{c}}\mathsf{\sigma }\mathrm{d}\mathsf{\omega }$$

(5)

With this approach, the softening response after the damage initiation, which is characterized by the reduction in the associated stress needed to open a crack as the crack length increases, is expressed by a stress-displacement response rather than by a stress–strain response. In the FEM, this stress-displacement concept requires a definition of a characteristic width $\mathrm{c}$ and the characteristic width has size effect on the fracture energy. The fracture energy given by Hu and Wittmann [21] from the tensile tests is a function of the dimensionless crack size, normalized by the specimen width, which shows that the measured fracture energy is significantly influenced by the dimensionless crack size. In order to eliminate the crack size effect, size-independent fracture energy should be used. In a further study the relation for the size-independent fracture energy was given as the following equation [22].

$$G_{\mathrm{F}}=\frac{K_{\mathrm{I}\mathrm{C}}^2}{E}$$

(6)

where $K_{\mathrm{I}\mathrm{C}}$ is the size-independent specific fracture energy. This equation enables the extraction of fracture energy without experiments and is free of the crack size effect.

In this study, the ABAQUS/explicit 2020 platform was used to predict the adhesive wear behavior. In the given wear modeling problem, this fracture energy is associated with the length of the mesh size in ABAQUS [23] since, in order to open a crack, the crack needs to go through an entire element at the crack tip. The sensitivity of the mesh was evaluated in [15] to determine the impact of the mesh size on the wear volume. It was shown in [15] that below a certain mesh size, a further reduction had a minor effect on the results but significantly increased the computation time.

Using the sub-model technique presented in [15], the boundary conditions of all the surfaces can be obtained from the deformation distribution of the global model. This means that all the surfaces are pre-deformed in the same way as they are in the global model. However, the rigid flat is not included in this sub-model, meaning that the normal load and tangential displacement cannot be calculated. Hence, it is impossible to capture the frictional behavior at the contacting surface. Moreover, although the less efficient global model presented in Ref. [13] was able to capture the friction behavior, as the predefined normal load decreases, a much smaller mesh size is needed, which requires a longer computation time. Therefore, a more efficient sub-model that can overcome the above limitations is suggested hereafter.

3. FE Model

The solution process consisted of two steps. First, a simulation using the global model without element deletion was performed to derive the size and deformation distributions for the sub-model, which is similar to the process described in [15]. In the second step, an improved sub-model was used with the obtained deformation distribution and element deletion to simulate the wear and friction behavior.

Figure 2 presents the global model with different sections consisting of different mesh densities and the location and assembly of the improved sub-model. The mesh design in the global model was the same as the one described in Ref. [13]. As shown in Figure 2a, three different sections are shown. Section I in red was where the potential wear particle was expected and had the finest mesh. Section II and III highlighted in yellow and blue, respectively, had a gradually increasing mesh size with an increasing distance from section I.

Using the sub-model technique [15], only a part of section I was used for the second solution step. This sub-model had a refined mesh compared to the mesh in the global model for better accuracy. As shown in Figure 2b, a plate-like sub-model highlighted in yellow with definitions of its thickness $t_s$ and radius $r_s$ in contact with a rigid flat is shown. The contact surface between the rigid flat and the improved sub-model is highlighted in blue.

The thickness $t_s$ and the radius $r_s$ of the improved sub-model were chosen based on the location of the maximum stress in section I of the global model. Due to the fact that the investigated material model is assumed to be perfectly elastic–plastic, these two parameters indicate the size of the obtained yield region. The thickness $t_s$ and the radius $r_s$ are defined as the vertical and tangential distances from the tip of the sphere to where the yield region is included.

These two parameters were slightly larger than those in [15], as their values were selected based on the location of the maximum strain found in the global model. Furthermore, using the geometry of the sub-model similar to [15] will result in the crack propagating beyond the sub-model boundary. This happens, since the element deletion was not included in the deformation distribution applied on all the surfaces of the sub-model in [15]. Thus, the interference in [15] was reduced and the thickness of the sub-model necessary to prevent the crack from penetrating the sub-model boundary was smaller.

Investigation on the effect of these two parameters, $t_s$ and $r_s$, was performed by separately conducting simulations with twice the thickness and radius. Although the differences in the wear volume and static friction coefficient were less than 7%, the computational time almost doubled. This indicated that the sub-model’s size determined by the maximum stress location was sufficient to effectively capture the wear process.

Figure 2c shows the boundary conditions applied to the improved sub-model. The contact surface was left free of constraints. A symmetric constraint was imposed on the symmetrical surface, highlighted in red, while the deformation distribution from the global model was applied to the side and bottom of the surface, highlighted in green. The deformation distribution from the global model serves two purposes: primarily, to produce contact between the plate-like sub-model and the rigid flat and to eliminate the effects caused by the geometry change, as the stress field depends on the geometry and the improved sub-model is not in a spherical shape.

It should be mentioned that an alternative solution was attempted, where instead of using the deformation distribution from the global model, fixed boundary conditions were used. However, these boundary conditions resulted in a significantly different crack initiation and propagation compared to Ref. [13]. As a result of the fixed boundary condition on the bottom, the sub-model was unable to deform in the normal loading direction, which led to a 50% reduction in the interference. Although the tangential behavior was primarily dictated by the behavior of the contact surface, it was also heavily influenced by the fixed boundary condition on the side, as this condition significantly altered the stress field. Moreover, the comparison presented in the results section demonstrates that the friction behavior obtained using the improved sub-model, which incorporated the deformation distribution on both the bottom and side, closely matched the results obtained in [13] using the global model. Therefore, the deformation distribution boundary conditions are necessary for accuracy and efficiency.

Finally, a full stick contact condition was defined between the rigid flat and the sphere. After the crack initiates, one more interaction must be defined, as there is a newly formed surface due to the failed elements deletion from the mesh. Eventually, a wear particle may form between the rigid flat and the slip interface, which is the interface between the newly formed wear particle and the bulk of the sphere. The interaction property applied on the slip interface is assumed to be frictionless to simplify and simulate. Thus, a shear crack cannot transmit shear or normal tensile stresses [24].

The elements with an aspect ratio close to one were used to reduce the mesh sensitively to the crack propagating direction; representing the uniform mesh as was also shown in [15]. In the present study, the mesh independence check was performed by increasing the mesh density and guaranteeing the relative error within a small, predetermined tolerance (10%) for all the normal loads and sphere radii. For instance, using R = 10 [mm] and $P^{*}$ = 100, changing the mesh size from 0.01 mm to 0.005 mm resulted in a change of the results below 9%. However, the calculating time of the latter was approximately twice that of the former. Therefore, the mesh size of 0.01 mm was used for R = 10 [mm] and $P^{*}$ = 100.

4. Results and Discussion

The adhesive wear and friction behaviors of the spherical contact under combined loading were studied by carrying out wear simulations using the improved sub-model for the sphere radius ranging from 5 [mm] to 50 [mm], with dimensionless normal loads $P^{*}$ ranging from 15 to 150. It was found in the current model that the tangential displacement $u_x=9{\mathsf{\omega }}_0$ was sufficient to achieve a wear particle for all the different simulations. It was also observed in [13] with the global model. Hence, to ensure the full formation of the wear particle, a value of $u_x=12{\mathsf{\omega }}_0$ was used for all the simulations. Similarly, in Refs. [13,15], the sphere material was aluminum 2024 T351, as this particular material was thoroughly studied both numerically and experimentally [25,26,27,28]. The material properties of aluminum 2024 T351 are: $E=74 \mathrm{G}\mathrm{p}\mathrm{a},v=0.33,Y=325 \mathrm{M}\mathrm{p}\mathrm{a}$ and a density of $\rho =2780 \mathrm{K}\mathrm{g}/{\mathrm{m}}^3$. The JC criterion coefficients of the damage initiation for aluminum 2024 T351 are (Ref. [28]): $D_1=0.13,D_2=0.13 \mathrm{a}\mathrm{n}\mathrm{d} D_3=-1.5$. For the damage evolution, the fracture energy was $G_f=20 \mathrm{K}\mathrm{J}/{\mathrm{m}}^2$ [28].

In the following, the verification of the improved sub-model is presented, followed by the investigation of the tangential mechanical response and the wear particle morphology for the different normal loads and sphere radii. The effect of changing the different mechanical properties ($D_1,D_2,D_3{,G}_f,E/Y$) on the obtained results is also discussed.

4.1. Verification of the Improved Sub-Model

The results from Ref. [13] were used to verify the adequacy of the present model. Figure 3 shows the five instants for the fracture evolution obtained by the improved sub-model, corresponding to the five instants obtained in Ref. [13] for the $R=10 \left[\mathrm{m}\mathrm{m}\right]$ and $P^{*}=100$ case with the definition of the wear particle parameters. The black color represents the free edge formed during the fracture evolution process; the deeper blue color represents the crack formed below the contact surface. Figure 4 shows the friction behavior and the interference for the same case, as shown in Figure 3. Both the dimensionless tangential force $Q/P$ and the dimensionless interference $\mathsf{\omega }/{\mathsf{\omega }}_0$ are plotted with respect to the dimensionless tangential displacement $u_x/{\mathsf{\omega }}_0$. The five instants of the fracture evolution are marked in vertical dash line.

For instant $I_A$, the crack initially appeared on the trailing edge of the contact surface. The term ‘trailing edge’ refers to the edge located in the direction opposite to the tangential displacement $u_x$, as indicated by the arrow shown in Figure 3a pointing in the $u_x$ direction. With an increasing tangential displacement, the crack spread along the edge of the growing contact interface, and at instant $I_B$, the tangential force reached its maximum, resulting in a sliding inception. At instant $I_C$, the first crack at the leading edge under the contact surface appeared. With the fracture extending below the contact surface, the fracture from the trailing edge linked with the one from the leading edge at instant $I_D$, as shown in Figure 3d. A region of un-failed elements that were not deleted according to the JC and fracture energy criteria was surrounded by the linked cracks. Eventually, at instant $I_E$ the un-failed region vanished leading to the wear particle creation. These five stages for the fracture evolution using the improved sub-model matched well with the results given in [13]. The differences in the corresponding tangential displacement for the five instants, as shown in Figure 4, were 4%, 6.8%, 5.4%, 2.3% and 1.5%, respectively.

Along with the crack initiation and propagation, the tangential force behavior of the improved sub-model corresponded well to the results from [13]. The tangential stiffness reduced with the increasing tangential displacement $u_x$ until instant $I_B$, at which the tangential stiffness completely vanished. The static friction coefficient represented by $Q_{max}/P$ was 0.36 for the improved sub-model and 0.35 in [13]. The relative error of the static friction coefficient between these two models was less than 5%.

Regarding the comparison between the dimensionless interference shown in Figure 4, the difference between [13] and the current result increased with the dimensionless displacement, which was due to the element deletion, as the element size used in the improved sub-model was larger than in [13]. However, the largest difference was less than 7%.

The wear particles formed in [13] and the one formed using the improved sub-model were similar in shape, having the shape shown in Figure 3f. The maximum thickness and length of the wear particle in [13] were ${\left(t_p\right)}_{max}=0.012R$ and $l_p=0.16R$, with a particle volume of $V_p=4.3e^{-2} \left[\mathrm{m}{\mathrm{m}}^3\right]$. The corresponding value using the improved sub-model were ${\left(t_p\right)}_{max}=0.0133R,l_P=0.155R$ and $V_p=4.36e^{-2} \left[\mathrm{m}{\mathrm{m}}^3\right]$. The differences between these three parameters were less than 10%, especially for the wear particle volume having a 2% difference. An additional parameter was also denoted in Figure 3f, which was the wear particle half width $W_p$, being $W_p=0.063R$ for this case.

To conclude, the results predicted by the improved sub-model were in a good agreement with the predictions made in [13], validating the improved sub-model for predicting the tangential behavior wear particle geometry. These results verified that the new sub-model can provide a reliable prediction for the wear morphology and the mechanical response, thereby validating the effectiveness of the present model.

The computation time was reduced from 90 h to 30 h, compared to the model used in [13]. A further reduction in the computation time can be achieved by reasonably increasing the loading rate of the simulation, both in the normal and tangential directions. Moreover, as shown in Figure 3, it can be observed that some non-contact and non-wear regions existed, which had no effect on the wear predictions when the tangential loading was terminated, indicating that the current sub-model simulating area can be further reduced for the purpose of a parametric study. Therefore, using a smaller sub-model simulating region compared to the present one, which fully covers the potential wear and contact regions and reasonably increases the loading rate, can further reduce the simulating time. For R = 10 [mm] and P* = 100, using the above accelerating techniques reduced the computing time from 30 h to 15 h, and the relative errors of the wear particle volume and the static friction coefficient compared to the current improved sub-model were less than 3%.

4.2. The Effect of the Normal Load $P^{*}$

Table 1. Geometrical and mesh parameters used in the global and improved sub-model.

${\mathit{P}}^{\mathit{*}}$	${\mathit{r}}_1^{\mathit{*}}({\mathit{r}}_1/\mathit{R})$	${\mathit{t}}_1^{\mathit{*}}({\mathit{t}}_1/\mathit{R})$	${\mathit{a}}_1^{\mathit{*}}({\mathit{a}}_1/\mathit{R})$	${\mathit{r}}_{\mathit{s}}^{\mathit{*}}({\mathit{r}}_{\mathit{s}}/\mathit{R})$	${\mathit{t}}_{\mathit{s}}^{\mathit{*}}({\mathit{t}}_{\mathit{s}}/\mathit{R})$	${\mathit{a}}_{\mathit{s}}^{\mathit{*}}({\mathit{a}}_{\mathit{s}}/\mathit{R})$
15	0.1	0.015	0.0015	0.07	0.008	5 × 10−4
20 to 30	0.1	0.018	0.002	0.08	0.012	6 × 10−4
50	0.13	0.022	0.002	0.1	0.02	7 × 10−4
75	0.2	0.05	0.003	0.14	0.026	0.001
100	0.2	0.05	0.004	0.15	0.038	0.001
150	0.25	0.07	0.004	0.19	0.05	0.002

summarizes the parameters used in the sub-model to study the effect of the normal load on the wear volume and the tangential load. The elements $r_1$ and $t_1$ are the radius and thickness of section I in the global model (Figure 2a), respectively. The elements $a_1$ and $a_s$ are the mesh sizes of section I in the global model and in the improved sub-model, respectively. Superscript $*$ is used to indicate the dimensionless parameter, normalized by the sphere radius $R$.

For accuracy, the mesh size $a_1$ should be decreased as the applied normal load decreases. On the other hand, in order to reduce the number of elements and the computation time, the radius and thickness of section I in the global model were reduced as well since decreasing the normal load resulted in a smaller potential region for the wear particle.

The parameters used in the sub-model were adjusted according to the stress field in the global model. As the simulations were conducted for the different radii, ranging from 5 to 50 mm, the parameters were normalized using the sphere radius R.

Figure 5 presents the results for the dimensionless tangential load $Q/P$ and the dimensionless interference $\omega /{\omega }_0$ vs. the dimensionless tangential displacement $u_x/{\mathsf{\omega }}_0$ for different dimensionless normal loads P* and different radii R. Figure 5 demonstrates that the transient dimensionless tangential forces and interferences exhibited a similar trend for the given material properties under different sphere radii and normal loads.

The dimensionless tangential force curves exhibited a certain non-monotonic behavior, particularly following the formation of the wear particles. These non-monotonic fluctuations were caused by the inherent numerical error. More specifically, the deletion and recontacting of the elements on the newly generated surface during the tangential loading led to the observed non-monotonic fluctuations. While it was impossible to entirely eliminate these fluctuations, they were reduced by implementing a more refined mesh. Similar fluctuations were also reported in Ref. [13]. Figure 4 provides a comparison that demonstrates that the fluctuation levels were within an acceptable range when compared to the previously published results.

As shown in Figure 5, the curves for the larger sphere radius were more consolidated. This behavior was caused by the dimensionless sub-model parameters ($r_s^{*}$, $t_s^{*}$). Although these two parameters allowed for the wear particle formation, the sub-model geometry might not have been large enough to exclude the influence of the applied boundary conditions. A few simulations of R = 5 [mm] were carried out using larger geometry, and the friction behavior was close to the friction behavior of R = 50 [mm]. It should be noted that the wear volume and wear particle parameters had a slight difference, below 10%. Since, as will be shown, the dimensionless results were independent of the radius, it was more beneficial to look at the larger radii as long as the adhesion, which was associated with nano and micron scale, was not considered.

Furthermore, it can be observed that the dimensionless interference during tangential loading increased with the increase in the dimensionless normal load. This can be explained by the junction growth under the combined loading [29].

Figure 6 illustrates the maximum dimensionless tangential forces, which indicate the static friction coefficient corresponding to the different sphere radii plotted against the dimensionless normal load. In general, the static friction coefficient decreased with the increasing dimensionless normal load, which showed a consistent result with [17]. A comparison between the results presented in [17] is also shown in percentages by the dashed lines. The differences for the small normal load, $15<P^{*}<25$, were lower than 10%; however, for the larger normal loads, $50<P^{*}<150$, the differences became larger, but were no more than 35% with respect to [17].

As shown in Figure 7, the maximum tangential force and the corresponding normal load are shown for the different sphere radii, where the different colors represent the different sphere radii.

A numerical fitted curve, in the black solid line, with $Q_{max}={\mu }_s\ast P$ resulted in ${\mu }_s$ = 0.37 for all the simulations conducted in the current study, regardless of the sphere radius. Interestingly, in Rabinowicz [30], for 50 separate experiments with different material combinations, a lower limit for the friction coefficient, approx. 0.37, was reported.

4.3. Wear Characteristics

The effects of the normal load intensity and the sphere radius on the wear particle morphology and wear rate were also analyzed. The geometry of a typical wear particle is presented in Figure 3f and the obtained results for the wear particle length $l_p$, thickness $t_p$ and half width $W_p$ are shown in Figure 8, in both dimensional and dimensionless manners. As shown in Figure 8, the dimensional and dimensionless parameters are indicated by the solid and the dashed lines, respectively. The same normalization as in [15] was adopted in the current study. All three wear particle parameters were normalized using the diameter of the contact area after the normal loading, $d_0$, which is a function of the radius and the interference in the following formula [16].

$$d_0=2\sqrt{{\omega }_{0}R}$$

(7)

It can be seen in Figure 8 that all the dimensional parameters of the wear particle increased as $P^{*}$ increased, indicating that increasing the normal load led to an increase in the wear particle size.

A consolidation of the dimensionless results can be observed. With the increases in the dimensionless normal load, the dimensionless wear particle length remained constant at a value of one. The behavior of the dimensionless wear particle length in this study was in contrast to that reported in [15], where $l_p^{*}$ was observed to increase with the dimensionless normal load caused by the junction growth. This difference can be attributed to the deletion of elements on the trailing edge of the contact surface in the current sub-model. Following the initiation of the crack, the element adjacent to the crack, as shown in Figure 3a, continued to be deleted from the wear particle. This deletion of the adjacent element ultimately resulted in a reduction in the wear particle length. Although some deletion occurred due to the small thickness at the trailing edge, it didn’t affect the prediction of the wear volume.

As shown in Figure 8b, the dimensionless wear particle thickness increased with the dimensionless normal load, and this behavior corelates well with [15]. The numerical curve fitting shown in the black solid line provided a function $t_p^{*}$ with respect to the $P^{*}$ with a goodness of $R^2=0.97$.

$$t_p^{*}=7.1e^{-4}{P}^{*}+0.021$$

(8)

It can be observed that the predicted thickness was slightly higher compared to the results in Ref. [15]. The differences in $t_p$ were caused by the differences in the obtained interferences in both sub-models. Including the rigid flat in the sub-model resulted in a slightly increased interference. As the interference increased, so did the size of the plastic region [31], forcing the crack to propagate more in the vertical direction and resulting in a thicker wear particle.

As shown in Figure 8c, the dimensionless wear particle half width was almost constant with a value of $W_p^{*}=0.5$, meaning half of the contact diameter $d_0$. This was due to the growth of the junctions primarily affecting the contact area in the tangential loading direction, which was perpendicular to the width direction of the wear particle. As a result, the wear particle width remained unchanged.

The wear volume $V_p$ is defined as the volume of the wear particle. The wear rate $w$ represents the wear volume per unit of tangential displacement, and it could be used for surface wear evaluation.

The wear rate is given by the following equation.

$$w=V_p/s$$

(9)

where $s$ is the sliding distance when the wear particle is fully detached from the sphere and is calculated as the following.

$$s=u_p+l_p$$

(10)

where $u_p$ is the tangential displacement of the rigid flat when the wear particle is formed. It is noteworthy that the simulation ran until the tangential displacement of $u_p$ and not s was achieved. With this method, the computational time can be reduced while ensuring a more realistic result. Obtaining a complete detachment of the wear particles by defining the sliding distance for the simulation required an extensive computational time and led to an unphysical wear volume. This was because the element deletion driven by the JC criterion and fracture energy continued after the wear particle was fully formed (at tangential displacement $u_p$). This led to an additional element removal from the wear particle during the detachment. Consequently, the wear volume for the complete detachment was smaller than the wear particle volume $V_p$ when it was just formed. The previous studies on the adhesive wear also utilized the same concept [13,15,32].

Figure 9 presents the results of the wear volume $V_p$ and the wear rate $w$ for the various dimensionless normal loads $P^{*}$ and different sphere radii R. The wear volume was normalized by the hemisphere volume, $V_p^{*}=V_p/V_0$, where $V_0=2\pi R^3/3$. The wear rate was normalized term by term, where the nominator was normalized by the hemisphere volume and the denominator was normalized by ${\omega }_0$.

$$w^{*}=V_p^{*}/\left({(u}_p+l_p)/{\omega }_0\right)$$

(11)

This normalization led to a general relation between the dimensionless wear rate $w^{*}$ and the dimensionless normal load $P^{*}$, allowing for the prediction of the wear rate and the wear volume independent of the radius for aluminum 2024 T351. The only parameter required to obtain the dimensionless wear rate was ${\mathsf{\omega }}_0$, which can be calculated by following equation [33].

$$P^{*}={\left({\omega }^{*}\right)}^{3/2}\left(1-\exp \left(\frac{1}{1-{\left({\omega }^{*}\right)}^{\beta }}\right)\right)$$

(12)

where $\beta$ is the linear function of Poisson’s ratio, $\beta =0.174+0.08\nu$, and ${\omega }^{*}$ is normalized by the critical interference ${\delta }_c$ for the stick contact condition, given by [33].

$${\delta }_c/{\mathsf{\omega }}_c=6.82\nu -7.83\left({\nu }^2+0.0586\right)$$

(13)

where ${\omega }_c$ is the critical interference at the yield inception for a slip contact condition [33].

$${{\omega }_c=\left(C_v\frac{\pi \left(1-{\nu }^2\right)}{2}\left(\frac{Y}{E}\right)\right)}^{2}R$$

(14)

Figure 9a shows that the wear volume increased as $P^{*}$ and R increased. The dimensionless wear volume curves were consolidated. In order to demonstrate the relation between $V_p^{*}$ and $P^{*}$, a numerical curve fitting with a goodness of $R^2=0.997$ was given as following equation.

$$V_p^{*}=3.819{\times 10}^{-9}{\left(P^{*}\right)}^2$$

(15)

By using Equation (15) and the normalization $V_p^{*}=V_p/V_0$, the dimensional wear volume can be calculated. The obtained result corelated well with the result in [15], with a difference less than 9%.

Figure 9b shows a similar behavior for the wear volume. The relation between $w^{*}$ and $P^{*}$ were given by a numerical fitted curve with a goodness of $R^2$ = 0.989.

$$w^{*}=6.39{\times 10}^{-12}{\left(P^{*}\right)}^{2.5}$$

(16)

By substituting Equations (12)–(14) into Equation (16), the dimensional wear rate can be calculated. The presented model deals with the creation of a wear particle from a well-defined spherical asperity. The other shapes of the contacting asperities or post-abrasive asperities were not examined in the present study.

Since there is an exchange of asperities in in the contact between rough surfaces in relative motion during the tangential displacement, and assuming that all these asperities are spherical, the wear rate for the entire rough surface was derived from the uniformed dimensionless result of one asperity obtained in this study.

The power-law dependency on the dimensionless normal load for the wear rate reported in [15] is given by the following equation.

$$w=3.26{\times 10}^{-5}{\left(P^{*}\right)}^{1.54}$$

(17)

The dimensionless relation obtained in the current study had the power of 2.5 in $P^{*}$ term. The reason for this difference was the normalization used for the dimensionless wear rate. The normalization used for the dimensionless wear rate was $w^{*}=(V_P/V_0)/(s/{\omega }_0)$. The denominator was normalized by ${\mathsf{\omega }}_0$, which was affected by the dimensionless normal load, as shown in Equation (12). The numerator was normalized by $V_0$, which was a constant. Hence, in Equation (16), a higher power on $P^{*}$ term was obtained.

To verify the expression obtained above, a simulation using the sphere radius, $R=15 \left[\mathrm{m}\mathrm{m}\right]$, was carried out for a dimensionless normal $P^{*}=100$. The dimensionless wear rate was $6.754\times$10⁻⁷, and the relative difference between the simulation result and the predicted results from Equation (16) was approx. 2.5%.

In Ref. [15], a transition region was observed for the wear rates. Three regimes of the wear rate with different exponents of $P^{*}$ were identified, which were $P^{*}\le 20$, $20<P^{*}\le 30$ and $P^{*}\ge 30$, representing the mild wear, transition region and severe wear, respectively. However, with the improved sub-model, this transition region was not observed. The power-law dependency on the dimensionless normal load for the wear rate $w$ in the current study is given by following equation.

$$\left\{\begin{array}{c}w=4.6{\times 10}^{-5}{\left(P^{*}\right)}^{1.5},\hspace{1em}{P}^{*}<30\\ w=2.74{\times 10}^{-5}{\left(P^{*}\right)}^{1.6},\hspace{1em}{P}^{*}>30\end{array}\right.$$

(18)

The change in the power of $P^{*}$ was not as significant between the two regions compared to what was reported in [15]. This difference was due to the different boundary condition at the contact surface of the sub-model. The deformation distribution used on the contact surface in [15] did not include the crack initiation and propagation, resulting in a smaller interference for the same normal load in the transition region. For instance, for R = 10 [mm] with $P^{*}=25$, the dimensionless interference when the wear particle was formed ${\omega }_p^{*}={\omega }_p/{\omega }_0$ was 1.2 in [15], while in the current study ${\omega }_p^{*}=2.1$.

4.4. Parametric Study for the Material Properties

The results presented above were limited to a specific material, aluminum 2024 T351. To obtain a more general understanding of how the material properties affected the wear particle formation and mechanical behavior, a parametric study was utilized. The initiation and propagation of a crack was influenced by the three coefficients $D_1,D_2 \mathrm{a}\mathrm{n}\mathrm{d} D_3$ in the JC failure model and by the fracture energy $G_f$.

The range of the studied coefficients $D_1$ and $D_2$ was from 0.05 to 0.2, and $D_3$ ranging from −1 to −2. These ranges covered the changes of approx. 50% of the original values for the aluminum analyzed in the current study. However, different materials can have values outside this range, for example FeCoNiCr high entropy alloy [34] had values of $D_1=0.004$, $D_2=0.665$, and $D_3=-1.5$ and free-cutting steel 50SiB8 [35] had values of $D_1=0.0733$, $D_2=0.7204$, and $D_3=-1.5643$.

The effect of $D_1$ on the friction behavior and wear volume $V_p$ was very small, and the maximum difference of the static friction coefficient and wear volume was less than 5% and 6%, respectively, for the range of the tested $D_1$ values. $D_2$ had a relatively large influence on the wear volume. The maximum change in the wear volume was 30% for the smallest value of $D_2$.

For $D_1$ and $D_2$, the most affected result was the wear particle thickness $t_p$. As these two parameters decreased, $t_p$ increased. As shown in Equation (3), $D_1$ and $D_2$ influenced the equivalent plastic strain in the same way. As $D_1$ and $D_2$ decreased, ${\mathsf{\epsilon }}_{\mathrm{f}}$ decreased. For the spherical contact under normal loading, the stress triaxiality ${\sigma }^{*}$ had a negative value as the element was being compressed, and it increased with the distance from the contact surface. This implies that the equivalent plastic strain for the crack generation decreased with an increasing distance from the contact surface. By reducing ${\mathsf{\epsilon }}_{\mathrm{f}}$, the crack propagation in the vertical direction became higher. Hence, the wear particle thickness increased as $D_1$ and $D_2$ decreased.

For the influence of $D_3$, the most affected values were $u_p$ and $V_p$. A maximum increases in $u_p$ about 60% was observed, compared to the original values of $u_p$ and $V_p$ for the smallest value of $D_3$. As $D_3$ decreased, the equivalent plastic strain increased for the negative value of ${\sigma }^{*}$, meaning that a larger strain is need for the element deletion, which explains the increase in $u_p$. The wear volume change was caused by an additional element deletion during the larger tangential displacement., which may have resulted in fewer physical results.

Figure 10 displays the friction behavior and dimensionless interference for various values of the fracture energy. While the JC parameters were more available in the literature, the same was not observed for the fracture energy. Hence, for the current parametric study on $G_f$, a range spanning two orders of magnitude around the fracture energy of aluminum was used. It is important to note that the fracture energy values shown in Figure 10 were normalized by the fracture energy of aluminum 2024 T351.

It can be observed from Figure 10 that the static friction coefficient increased with the fracture energy. As shown in Figure 3, the maximum value of the dimensionless tangential force, or the static friction coefficient, was achieved after the crack initiation. This suggests that the fracture energy had an impact on the static friction coefficient. Since the force required to open a crack increased with the fracture energy, the static friction increased as well. However, it should be noted that the influence of the fracture energy was relatively small, with a maximum difference of only 6.2% observed between the two orders of magnitude of the fracture energy.

In contrast, the dimensionless interference decreased with the fracture energy. The change in the interference consisted of two contributors, namely the junction growth and element deletion. With an increase in the fracture energy, the number of elements deleted decreased, leading to a reduction in the dimensionless interference.

While changing $D_1$, $D_2$, $D_3$ and $G_f$ caused some changes to the wear particle morphology, the wear particle formation process for the different parameters was similar to those presented in Figure 3.

However, for a different material, Ti-6Al-4V, a difference in the crack propagation behavior was observed, as presented in Figure 11. The material properties were [36]: $E=110 \mathrm{G}\mathrm{p}\mathrm{a}$, $\nu =0.33$, $Y=880 \mathrm{M}\mathrm{p}\mathrm{a}$ and a density of $\rho =4430 \mathrm{K}\mathrm{g}/{\mathrm{m}}^3$. The JC criterion [36] and the fracture energy [37] were: $D_1=-0.09$, $D_2=0.25$, $D_3=-0.5$ and $G_f=47 \mathrm{K}\mathrm{J}/{\mathrm{m}}^2$.

As shown in Figure 11, at instant $I_A$, a crack was generated on the trailing edge of the contact surface. At instant $I_B$, the tangential force reached its maximum. These two stages are similar to Figure 3. However, at instant $I_C$, instead of a crack generated below the leading edge of the contact surface, the crack propagated along the edge of the contact surface and reached the leading edge. Interestingly, at $I_D$, a second crack was generated within the contact surface surrounded by the first crack. As the tangential displacement increased, the second crack linked with the crack at the trailing edge of the contact surface and led to the formation of a wear particle at instant $I_E$.

Since it was not observed that changing the JC coefficients or fracture energy led to such changes in the wear particle formation process, it was assumed that this new crack propagation may have been caused by the change in $E/Y$. For aluminum 2024 T351 $E/Y=227$ and for Ti-6Al-4V, this value was $E/Y=125$.

5. Conclusions

An efficient FE model was presented for adhesive spherical contact under combined normal loading and tangential displacement. The model consisted of a global model that helped determine the potential fracture location and a sub-model with a refined mesh for the precise prediction of the wear particle formation. The improved sub-model included the rigid flat and two criteria were employed to simulate the crack initiation and propagation. Thus, the present sub-model inherited the efficiency of the sub-model introduced in the literature and the capability to predict the tangential mechanical response.

The current study primarily focused on the effects of the sphere radius and a range of relatively high normal loads on the wear particle evolution. However, the developed method can be extended to the wear simulations under various parameter combinations, including the materials, operating conditions, etc. The main contributions of the present study were the unified dimensionless relation between the dimensionless wear rate and the dimensionless normal load for the sphere radius varying from 5 to 50 [mm], as well as the insights gained into how the JC parameters and material properties (E/Y) affect the friction and wear behaviors. The obtained relations for the wear volume and wear rate may be generalized for the wider range of mechanical properties with the presented efficient sub-model.

Furthermore, the present method offers the possibility of extending the single spherical wear model to the adhesive wear prediction for rough surfaces and for fretting wear modeling. Nevertheless, the modeling of the adhesive wear can be improved by introducing a physical approach for the crack propagation, which is not based on element deletion. Such an approach would provide a more realistic and physically based simulation of the crack propagation and wear particle formation.