An advanced efficient model for adhesive wear in elastic—plastic spherical contact

A finite element (FE) model combining submodel technique is presented for the adhesive wear in elastic—plastic spherical contact. It consists of a global model, showing the potential location of fracture under combined normal and tangential loading, and a refined mesh submodel covering only the region near the potential fracture. This allows to describe the morphology of wear particle more accurately than that in a previously developed model by the authors. A range of normal loading is studied to show its effect on the shape and volume of wear particles. Two main regimes of mild and severe wear (along with a relatively narrow transition region between them) are found, which show almost linear and power-law dependency of wear rate on normal loading, respectively. Such behavior agrees with published experimental observations. However, the transition region is theoretically predicted here for the first time.


Introduction
Among the various wear types, the adhesive wear is argued as the least avoidable one [1,2], which results in material transfer between contacting asperities of the mating surfaces due to strong adhesive bonds there [3]. According to the well-known Archard equation [4], the amount of material lost from a surface during sliding contact (wear volume) is proportional to the normal load applied to the surface. However, many following experimental studies found that the linear relation is maintained only in a certain range of normal loads and it fails for lower or higher load limits [5][6][7][8]. Above the higher load limit, a transition from mild to severe wear occurs. This experimentally observed transition is commonly reported in the relevant literature as an abrupt one see, for example, Fig. 1 in Ref. [9]. However, somewhat different experimental results showing a narrow transition region are reported in Ref. [10]. Unfortunately, due to the formidable com-plexity of adhesive wear mechanism [2,11], accurate modeling including the above-mentioned transition still remains one of the challenges in tribology.
Suh [12] proposed a qualitative explanation for adhesive wear indicating that wear occurs by "delamination" of sheets caused by subsurface deformation, crack nucleation, and crack propagation. In Refs. [12,13] such flake-like wear sheets were experimentally observed. Then Suh and co-workers proposed quantitative models for subsurface crack nucleation and propagation [14,15]. However, in Ref. [14] an a priori imaginary inclusion was assumed and in Ref. [15] the linear elastic fracture mechanic (LEFM) was used, which ignored possible plastic deformation at the crack tip [16].
Inspired by the pioneering work of Rabinowicz [17], Aghababaei et al. [18] used atomistic simulations and found a characteristic length scale d * that controls the adhesive wear mechanism at asperity level. Junction sizes below d * produce plastic flattening for asperity, hardness of the sphere l P length of wear particle L c critical load at yield inception in full stick contact condition P normal load P * dimensionless normal load, P * = P/L c R radius of sphere s sliding distance t P thickness of wear particle u P tangential displacement at formation of wear particle u x tangential displacement V 0 volume of original hemisphere, V 0 = 2/3πR 3 V p volumes of wear particle Y yield strength


Poisson's ratio ω Interference associated with low wear regime and sizes above d * produce brittle fracture induced wear, associated with mild wear regime. With further simulations [19], they investigated the microscopic origin of wear transition at the asperity level and found that interaction between subsurface stress fields of neighboring contact spots promotes a transition from mild to severe wear. In the mild wear regime, asperities are detached in the form of tiny particles, the size of which is comparable to the junction size and hence, it follows the linear relation of Archard equation. In the severe wear regime, large wear debris occurs due to deep crack propagation below surface contact. A single particle size in this regime corresponds to the apparent contact area of multiple asperities and hence, the linear relation of Archard equation fails. Following the work of Rabinowicz [17] and Aghababaei et al. [18], Popov et al. [20] introduced an "asperity-free" wear criterion to numerically solve the wear evolution of rough surfaces. They identified two types of wear: One leading to the formation of a modified topography, which is too smooth to wear further and one showing continuously going on wear. A power-law dependency of wear volume on normal load was found with exponent larger than one.
Recently, the authors developed a finite element (FE) model for adhesive wear in elastic-plastic spherical contact [21]. This model was based on realistic ductile failure criteria in a three-dimensional (3D) contact problem and failed elements are deleted. Flake-like wear particles were observed and the wear coefficient of Archard equation was obtained. However, due to extremely long computing time, only two cases of normal load were investigated in Ref. [21], incapable to reveal the actual dependency of adhesive wear on normal load. In addition, due to the failed elements deletion, unphysical material loss (at the slip interface) occurred, leading to overestimated interference of the sphere. The "submodel" technique [22], used in Ref. [23] for fretting contact, may be a better choice to avoid the above deficiencies in Ref. [21].
The main goal of the present study is to investigate theoretically the dependency of adhesive wear in spherical contact on normal load using an advanced, more efficient, 3D FE model than the previously developed one in Ref. [21]. A range of normal loading is studied to show its effect on wear rate. The results are compared with published experimental results available in the literature and the transition region from mild to severe wear regime in adhesive wear is theoretically predicted for the first time. Figure 1 presents a schematic of the contact problem between a homogeneous deformable sphere and a rigid flat under combined normal and tangential loading. Initially a load-control normal loading P is applied on the rigid flat and then a displacement-control tangential loadingu x is applied in a stepwise manner. Increasing u x is associated with increasing interference ω and a slip interface, as indicated by the dotted line in Fig. 1(a), associated with the formation of wear particle is obtained at interference ωω p . At this instant the | https://mc03.manuscriptcentral.com/friction volume of wear particle above the slip interface is V p . The cross section (plane y = 0) of the resultant wear particle with its length l P and thickness t P are presented in Fig. 1(b). In the final step of tangential displacement loading u x = s, the wear particle, attached to the rigid flat, will be fully detached from the sphere bulk (schematically described by dashed lines). In Ref. [21], the adhesive wear of the spherical contact presented in Fig. 1 was studied by the authors with the FE method. To model the fracture, ductile material failure criteria were applied using two steps: The Johnson-Cook (JC) criterion [24] for damage initiation in the first step, and fracture energy criterion [25] for damage evolution in the second step. The layer of failed elements (indicating a slip interface) according to these criteria was deleted from the mesh. Thus, the active elements remaining between the flat and slip interface specify the wear particle.

Theoretical model
In the mesh design of the sphere in Ref. [21], a Zone I near the contact region (see global model in Fig. 2) was created with refined cubic elements of size a 1 to better describe the fracture and wear particle. The size  of Zone I was determined by trials, which is inefficient in terms of computer time. This deficiency can be resolved by implementing the submodel technique in the FE method [22]. The submodel technique is used to study a locality of a global model with a refined mesh based on interpolation of the solution of the relatively coarse mesh global model. Hence, it is very efficient in obtaining an accurate detailed solution in a local region of the global model [22]. With this technique, a submodel inside Zone I of Ref. [21] can be created just after the results of the global model (the whole hemisphere) are obtained. The size of the submodel can be properly determined, without trials, based on the coarse distribution of plastic strain in the global model. Thus, the whole solution process could be more efficient.
The element deletion used in the FE model in Ref. [21] also causes overestimation of the interference ω due to the nonphysical material removal at the slip interface, which actually should have zero thickness. The submodel technique helps to resolve this issue of element deletion, which affects the predicted volume of wear particle. To accomplish this, the element deletion in the global model is disabled. Then, using the obtained results from the global model, the submodel with refined mesh and resumed element deletion can better predict the fracture evolution and the final wear particle volume.
Elastic perfectly plastic material behavior is assumed for the global model as in Ref. [21] and its loading process is identical to that described in Ref. [21], but the ductile material failure criteria are not used since element deletion is not needed. Full stick contact condition is assumed for the interface of sphere and flat, which means that relative displacement of points engaged in contact is prevented.

The submodel
The location of the submodel, inside Zone I of the global model in Fig. 2(b), and its mesh design are presented in Fig. 3. To predict the fracture evolution, which requires much smaller element size than a 1 of Zone I, a plate-like submodel (see Fig. 3(a)) with radius, r s , and thickness, t s (see Fig. 3(b)), is created near the sphere tip inside Zone I of the global model. The boundary conditions for the submodel at every loading step are known from the prior obtained deformations distribution of the global model. The values of r s and t s are selected to be just enough to cover the potential fracture location, based on the already obtained location of maximum plastic strain in Zone I of the global model (see Ref. [21]). The cubic element size a s (a s < a 1 ) in the submodel is selected to be small enough to obtain accurate volume of the wear particle as explained at the end of this Section.
The material of the submodel is assumed elastic perfectly plastic with ductile material failure criteria same as those in Ref. [21]: JC criterion [24] for damage initiation and fracture energy criterion [25] for damage evolution. A failed element according to the failure criteria is completely deleted from the mesh. As indicated in Ref. [21], the deficiency of element deletion can be minimized by using sufficiently small element size, which requires longer computing time. The submodel, has fewer such small elements and hence, provides a more efficient way for higher accuracy at a shorter computing time.
In the submodel, the evolved slip interface (the dotted line in Fig. 1(a)) due to deleted failed elements inside the material is assumed frictionless based on Ref. [26] (see also [21]). The "loading" of the submodel is obtained by the deformations of its boundaries, which are transferred from the global model. With the progressing deformations of the submodel, the slip interface and the wear particle are formed when increasing interference due to the increasing tangential loading reaches ω = ω p (see Fig. 1(a)). The sliding distance s for a fully detached wear particle from the sphere is assumed, as was done in Ref. [21] for simplicity, as the sum of u x at the instant of the wear particle formation and the wear particle length l P . The volume of the wear particle, V p (see Fig. 1(a)), is obtained by summing up the volumes of all active elements within it and its dimensionless value * The accuracy of the global model was already validated in Ref. [21] by comparing the numerical results with the analytical solution of Hertz [27] for an elastic contact. Convergence of the numerical solution for the submodel is assured by sequentially reducing its elements size a s until the change of V p becomes less than 5%.

Results and discussion
A range of dimensionless normal loads P * from 1 to 200 was investigated. Using dimensionless parameters for presenting results renders the model a sort of universality. Here P * is defined as P * = P/L c where L c is the critical normal load at yield inception of a sphere under full stick contact condition as given in Ref. [28]: Aluminum 2024 T351 is selected for the sphere since the failure criteria for this particular material have already been studied by experiments [29][30][31] and are well documented in Refs. [32][33][34]. As in Ref. The results for P * = 100 are discussed first as a typical | https://mc03.manuscriptcentral.com/friction case in order to demonstrate the efficiency of the present model compared to that in Ref. [21]. The values of r 1 , t 1 , a 1 , r s , t s , and a s are normalized by R and indicated with the superscript *. For the global model, the values of * 1 r = 0.2 and * 1 t = 0.05 are used for Zone I and provide large enough region to cover the potential fracture. A value of * 1 a = 0.004 (compared to the much smaller * 1 a = 0.001 in Ref. [21]) is used, which provides high enough mesh density for accurate deformations (the interference ω 0 is only 2.1% different compared to ω 0 with *  V at a s /R approaches zero when the element size decreases. The corresponding cross sections in green color (see online version) of the wear particle at its middle plane (plane y = 0, see Fig. 2) are also shown, along with the associated computing times to solve the submodel.

Effect of submodel element size
As can be seen from Fig. 4, smaller * s a takes longer computer time and increases * P V with a decreasing rate. Smaller * s a means smaller material loss at the slip interface due to deleted elements. It is expected that for infinitely small * s a , the material loss will tend to zero and hence, * P V will converge to its expected accurate value. This corresponds to the bright circle at the intersection of the best fit curve with the vertical axis showing the value of * P V = 3.74e-5. Compared with this value, * P V for the other five cases of * s a between 0.002 and 0.0005 presents differences of 13.7%, 9%, 4.7%, 1.4%, and 0.51%, respectively. Therefore, it can be summarized that in the submodel the element size * s a = 0.001 could be the preferred choice for P * = 100 providing accurate enough results for volume and shape of the wear particle with relatively short computing time.
For * s a = 0.001 the dimensionless wear volume in Fig. 4 is * P V = 3.57e-5, which is 13% less than the predicted value of * P V = 4.1e-5 in Ref. [21] using the same element size for * 1 a there. This difference is due to the element deletion used in Ref. [21], which overestimated the sphere final interference and hence, the wear particle size. Also, the present computing time is only 18 h compared to the 90 h required for the previous model in Ref. [21] using the same computer with the same element size in Zone I. This comparison clearly demonstrates the better accuracy and efficiency of the present advanced model.

Effect of normal load P *
The effect of normal loading on the wear volume is studied for a range of P * between 1 and 200. The input parameters used for the global model and the www.Springer.com/journal/40544 | Friction submodel for all the studied cases were selected, as was described for P * = 100 above, and are listed in Table 1. The output parameters for the dimensionless diameter of the contact area prior to tangential loading, * 0 d = d 0 /R, dimensionless length, thickness and volume of the wear particle, * P l = P l /d 0 , * P t = P t /d 0 , and * P V , respectively, are listed in Table 2. The obtained numerical results for the dimensionless sliding distance at the wear particle formation, * P u = P u /d 0 , the ratio of thickness to length of the wear particle, * P t / * P l , are also listed in Table 2.
It can be seen from Table 2 that both the length * P l and thickness * P t increase as P * increases, representing increasing size of the wear particle. However, * P t is always much smaller than * P l ( * P t / * P l < 0.07), Table 1 Input parameters used in the present model.  representing slender flake-like wear particles for all the 7 normal load cases. The ratio of * P t / * P l increases as P * increases, representing thicker wear particle for larger P * . It can also be seen that the wear volume * P V increases with increasing P * as expected.
Instead of the wear volume, the wear rate w representing the wear volume per unit sliding distance, is used in many prior studies (see Ref. [8]) to investigate adhesive wear. This parameter is given by: where s = P u + P l is the sliding distance when the wear particle is just fully detached from the sphere [21]. Figure 5 presents the numerical results of wear rate w for the various values of P * in Table 1. The solid dots correspond to the numerical results with their best fit curve presented by the solid line. As shown in Fig. 5, when P * increases, w increases exponentially with a power-law dependency on normal load, similar to the finding in Ref. [18], having the form: where C and n are constants depending on material properties and contact conditions. The goodness of fit R 2 in Fig. 5 is 0.998. However, for P * ≤ 20, the average error is 41%. Hence, Eq. (3) is not valid over the entire range of w, which covers three orders of magnitude (see the log-log scale in Fig. 6). To resolve this, the entire range of P * is divided into three regimes: for P * ≤ 20, 20 < P * < 30 and P * ≥ 30, respectively, see 2.7e-6( *) , 20 * 30 3.1e-5( *) , * 30 Equation (4) is a much better fit showing goodness of fit R 2 larger than 0.999 and average error less than 2% for each of the three regimes of P * . In the range of 20 < P * < 30, a narrow transition region is observed connecting the two main wear regimes.
It can be seen that for P * ≤ 20, with n = 1.08 the wear rate is almost linearly dependent on the normal load, as in the Archard equation. For P * > 30 with n = 1.56 there is an obvious deviation from the Archard equation.   Similar power-law dependency of wear volume on mean contact pressure for contacting rough surfaces is theoretically predicted in Ref. [20] with power n in the range between 1.36 to 1.77.
In the narrow transition region of 20 < P * < 30 from mild to severe wear, an increase of w occurs similar to the experimental finding in Ref. [10] for oxidative wear (see Fig. 7 in Ref. [10]). It is important to notice here that the power n = 2.22 in this transition region is even higher than that in the severe wear regime. In the present model, the increase of w is associated with similar increase of wear particle slenderness, * P t / * P l . As shown in Table 2, for P * < 20 the values of * P t / * P l are less than 0.01 indicating extremely thin wear particles, but for P * > 30 the values of * P t / * P l are 5 times larger.
Some experimental results are available in Refs. [8,35,36] showing wear rate vs. load for Al alloys. Similar two regimes of mild and severe wear are observed with power n in the range of 0.56-0.6 for the mild wear regime, and 1.36-2.1 for the severe wear regime. The present theoretical values of n fairly agree with these experimental results. The differences may be due to the different contact types and material hardness values H being: block on ring with H = 450 MPa in Ref. [8], shaft on bearing with H = 380 MPa in Ref. [35], and indenting flat-end pin on disk with H = 95 MPa in Ref. [36], compared to the present flattening spherical contact with H = 910 MPa.

Conclusions
1) An advanced finite element (FE) model was presented for the adhesive wear of spherical contact under combined normal and tangential loading. It consists of a global model, showing the potential location of fracture and a submodel covering only the region near the potential fracture with refined mesh. Damage initiation and evolution in the submodel are considered and failed elements are deleted from the mesh, allowing to more accurately simulate the wear particle formation.
2) The present model eliminated some deficiencies found in a previously developed model [21]. With smaller number, and properly selected size of smaller elements, in the submodel, the results for volume and shape of the wear particle are improved and the computing time is substantially reduced compared to the older model [21].
3) For different values of dimensionless normal load P * , flake-like wear particles were observed and their volume increased with increasing P * . Two main wear regimes along with a transition region between them were found showing almost linear and powerlaw dependency of wear rate, w, on P * , respectively, as presented in Eq. (4), similar to such experimental observations. 4) It should be noted here that only a single case of sphere radius, material properties, and contact geometry was considered in the present study. Investigating the effects of sphere radius, material properties, as well as sliding distance, etc., on the constants C and n of Eq. (3) requires a major effort, which is beyond the scope of the present study. However, a wide prospect now exists for future work based on the present advanced model.