3D finite element analysis of a two-surface wear model in fretting tests

This article aims at developing a computationally efficient framework to simulate the erosion of two contact surfaces in three-dimensional (3D), depending on the body resistance. The framework involves finite element (FE) resolution of a fretting problem, wear computation via a non-local criterion including a wear distribution parameter (WDP), as well as updating of the geometry and automatic remeshing. Its originality is based on the capability to capture the damage on each surface and obtain local and global results for a quantitative and qualitative analysis. Numerical simulations are carried out for two 3D contact specimens with different values of WDP. The results highlight the importance of correctly modelling wear: One-surface wear model is sufficient from a global point of view (wear volume), or whenever the wear resistance for a body is much higher than that of another one, whereas a 3D two-surface wear model is essential to capturing local effects (contact pressure, wear footprint, etc.) related to the difference in wear resistance of the bodies.


Introduction
A large number of engineering systems such as rivets, turbine blade fixings, shafts, and steam generators are subjected to small amplitude motions that generate a fretting problem. This tribological phenomenon can lead to the degradation of material surfaces, such as wear on the two contact surfaces, and then reduces the service life of the systems.
It is therefore an interesting goal to aim at developing a robust numerical environment in order to predict surface damage.
First of all, the wear model must be based on the experimental observations, which highlight, on the one hand, a deterioration mechanism on one or both contact surfaces, depending on the wear resistance of the antagonistic materials, and on the other hand, the need for a three-dimensional (3D) wear model to capture the local effects.
The second point concerns the choice of the wear law. Classically, many authors use a derived version of Archard's law to predict wear effect, which means that they refer to a criterion of dissipated energy as Ref. [1], even if more than one hundred models can be listed [2].
The straight application of this criterion only allows to model the deterioration of a contact surface or a balanced distributed wear loss between the contact surfaces. Therefore, this law must be appropriately adapted in the context of our study.
A review of the literature reveals that although most of studies are limited to modelling wear on one-surface [3][4][5][6][7][8][9], a few papers deal with wear on two contact surfaces [10][11][12][13][14][15]. For this purpose, McColl et al. [10] have adapted the wear criterion by simply introducing an arbitrary wear coefficient for each surface, without worrying about the differences of resistance between the bodies. This last point seems to be present in Ref. [11] since the authors define a new parameter called the wear-mode index that indicates the relative magnitude of the wear coefficients of the contact pair. This index depends on the wear coefficients, the geometry, and the sliding distance. However, since the wear coefficients are considered independent of each other, the study of the different wear resistance is not exhaustive. To the best of authors' knowledge, the only work that takes into account the body resistance is presented in Ref. [13]. The authors propose a predictive wear model that characterise the volume loss of the entire contact pair by means of a wear distribution parameter (WDP,  ). But actually, this model uses a local description of the Archard's wear law and is limited to a global analysis of a two-dimensional (2D) problem.
The goal of this work is to develop a robust numerical environment to simulate the 3D wear of fatigue-fretting, when two contact bodies do not have the same wear resistance. The originality of the present work consists in introducing computational strategies for damage post-processing in order to obtain both local and global quantities and to analyse the results qualitatively and quantitatively. This numerical procedure makes it possible to estimate the wear volume for each surface in contact while constructing a precise vision of the footprint created by surface erosion. A 3D numerical investigation is performed for several values of wear resistance with two types of geometries. The material is a titanium-based alloy Ti-6Al-4V.
The paper is organized as follows. Section 2 consists of the strategy of simulation, the description of the new numerical wear model, and finally the mesh update procedure. Section 3 describes the finite element (FE) model, namely geometry and meshes, material parameters, and boundary conditions. The numerical results are assessed in the rest of the paper. The effect of WDP in return is examined. The hysteresis loops, the wear profiles, the wear volume, the contact pressure, and the comparison between wears are brought up and discussed in respect to the literature results.

Wear numerical model 2.1 Strategy of the wear simulation technique
Section 2.1 is devoted to the general description of the procedure for simulating wear evolution. The incremental process of the wear modelling requires a three-step resolution. First, the calculation is performed with the FE solver of Abaqus © Ref. [16]. The contact algorithm uses an impactor/target technique, and the friction behaviour is modelled by a Coulomb's law. Then, the contact pressure and the sliding distance can be evaluated for the purpose of wear calculation on the slave surface. Secondly, the implementation of the 3D wear model is developed in the numerical environment of the z-set suite [17], which controls mesh generation, FE calculation, and post-processing. The post-processing allows to determine the wear depth to be applied at each node of the contact surfaces, using the non-local version of the Archard's law [18]. The concept of a WDP is then introduced to estimate the percentage rate of the erosion for each contact surface. The updating of the geometry is made by shifting the contact nodes with a vertical translation of the wear depth. Finally, an updated mesh is generated by means of the code GMSH [19]. Thus, after p incremental process, the pth updating of the contact surfaces leads to a new geometry configuation called "profile p".
The strategy of the simulation process is illustrated in Fig. 1. Only the implementation of the step 2 is set out in Section 2.2.

Wear model implementation
The wear model is based on an average friction energy approach. This approach is implemented in Ref. [20] for a one-surface wear model. For completeness and conciseness, this model is briefly recalled.
At each node i of the slave surface, the wear depth The exponent s refers to the slave, and ij d represents the distance between the nodes i and j. The calculation is made for predefined values (smoothing parameters) of  and R. The identification of these parameters requires to carry out a preliminary parametric study, whose outcome is a set ( , Then, after N numerical fretting cycles (N is the profile number), the cumulative wear at node i on each contact surface becomes For each wear profile p , the contact nodes are moved depending on the wear depth  to the contact surface ( Fig. 1), and their new position is updated by Eq. (7): Finally, the volume loss ( , ) , m,s p V    of each specimen after p profiles is calculated by a postprocessing in z-set and allows to get the change in volume due to the adaptative meshing. Figure 3 illustrates the computation process. Each FE wear profile consists of three steps: (1) the contact solver with the increment sliding loops; (2) the evaluation of wear depth and of the wear volume; and (3) the geometry update and the remeshing.

FE mesh and boundary conditions
Two fretting wear tests are considered: (1) C/P with a cylinder radius of 10 mm; (2) P/P with a flat punch surface of 3 mm and a curvature radius of 1.35 mm. The lateral widths of the master and slave are 2 and 3 mm, respectively. Only one half of the geometries are considered due to the symmetry in respect of the mid plane z 0 (Fig. 4). As mentioned in Ref. [22], the choice of the mesh size at contact is crucial to produce satisfactory stress results, since they are of great importance in the estimation of the lifetime. Therefore, a free mesh is made of four-node linear elements and refined towards the contact area with the element size of around 50 μm. This choice of the element size in the contact zone is motivated by Ref. [18].
The loading history, modelled to simulate the fretting cycle, is illustrated in Fig. 4(c), and is carried out as follows: (1) A normal force P at point A and a multi-point constraint condition to have the same vertical displacement on the top of the master are imposed, while the side and bottom edges of the plate are fixed; (2) a periodic horizontal cyclic displacement   is imposed on the plate when the normal force is kept constant.
The boundary conditions correspond to a mechanical response under gross slip conditions. The value of the normal force P is 300 N/mm, and the displacement amplitude   is 25 μm for the C/P configuration; for the P/P configuration, P  2,470 N/mm, and    150 μm. For each profile, the updated mesh takes into account the wear of the master/slave surfaces.

Material parameters
The specimens are made of Ti-6Al-4V titanium alloy. The behaviour is linear elastic, with a Young's modulus E of 119 GPa and a Poisson's ratio  of 0.3. For the sake of simplicity, the friction and wear coefficients will be taken as constants.  = 0.8 is used to represent the average value of a dry Ti 6 Al 4 V-Ti 6 Al 4 V contact, as measured in the experimental tests [23,24]. An average value of the wear parameter k w = 1.5 × 10 −5 mm 3 /J is representative of the fretting contact pair [25].

Results and discussion
The model is able to provide qualitative and quantitative information for arbitrary geometries. Specifically, the wear volume and the wear footprint are analysed on master and slave contact surfaces for www.Springer.com/journal/40544 | Friction several wear distributions. Particular attention is paid to the influence of the wear distribution factor on the contact conditions, namely numerical loops, contact pressure, wear profile, and wear volume. Since the numerical process introduces several parameters to reduce the computational cost, the convergence analysis is a prerequisite to be able to guarantee the quality of the results. A study is therefore conducted for each simulation, and a set of parameters is determined in the same way as Ref. [20]. The values provide stable wear profile evolution without requiring long simulations.
Note that in order to validate this implementation of the wear distribution model, the results with 0   are compared to those of Ref. [18].

Hysteresis loops
Simulated tangential force-displacement loops are illustrated for two contact geometries (C/P and P/P), three FE wear profiles (  1, 25, 50 p ), and various values of WDP (     0, 0 25, 0 5, 0 75, 1  ) (Fig. 5). The hysteresis loops for 0   (wear only on the plate) are consistent with Ref. [20]. For a given numerical profile p, the curves coincide whatever the value of WDP. A large hysteresis with a parallelogram shape is systematically observed and confirms the gross sliding conditions (Fig. 5(a)). The conclusions are the same for the two specimen shapes. Nevertheless, the horizontal part of the loop corresponding to the sliding period does not tilt during the fretting cycles because just a few cycles are simulated. This is consistent with the fact that the present study is mainly devoted to the influence of the distribution parameter on damage response. During the fretting test, the area of hysteresis, i.e., the dissipated energy, increases. These simulations allow to reproduce the experimental observations for similar loading conditions qualitatively, e.g., for C/P case [3,6,10,25,26], and for P/P case [27,28], which are very encouraging for the wear simulation technique.
In conclusion, the wear distribution factor has no influence on the tangential force-displacement loops, demonstrating that the energy dissipated is the same regardless of the wear process.

Evolution of the geometry of the contact surfaces
Section 4.2 is devoted to the comparison of wear footprints according to the wear distribution factor. First, 3D plots give an overview of the shape of the eroded area, and allow us to demonstrate the smoothness of the profiles (Fig. 6). The plots are shown for two values of  (0 and 0.5), and for both C/P and P/P geometries. For the sake of clarity, a translation along z-axis has been applied to the profile of the indenter (cylinder or punch), so that the traces are not superimposed. Of course, there is no wear on the indenter for 0   , and in this particular case, the erosion reaches a maximum in the middle of the print for C/P case, and at both contact ends for  www.Springer.com/journal/40544 | Friction P/P. This trend is preserved for 0 5    , but now the erosion is present on the plane and on the indenters. The same results are represented on the maps for several values of the wear distribution factor  at the 50th profile for the cylinder-plate case (Fig. 7) and at the 10th profile for the punch/plate case (Fig. 8). In both cases, the plate is in the left column, meanwhile the right column shows the results on cylinder or punch. The contact width does not depend on wear resistance of the contact bodies and the  parameter, but the damage on each body strongly depends on  .
For 1   (the first row), the only eroded body is the cylinder (Fig. 7) or the punch (Fig. 8) The footprint has similar shapes for the two bodies. More precisely, wear is well balanced for 0 5    , in the sense that the footprint has the same shape and depth on both bodies. As far as the cylinder-plate configuration is concerned, the maximum wear depth takes place  systematically in the middle of the print in the xy-plane whether on the plate or the cylinder: A "U-shape" defect is obtained in the direction perpendicular to the sliding direction z (Fig. 7). On the plate, wear is lower in plane ( 1 z  ) than in the symmetry plane ( 0 z  ). The maximum wear for the cylinder is located at the contact edge, i.e., in plane 1 z  . In the punch-plate configuration, a typical "W-shape" profile on the xy-plane is combined with a "U-shape" profile in the yz-plane. The maximum of the wear depth is always noticed at the lateral borders in the xy-plane for the two contact bodies. Wear is the maximum in plane 1 z  (Fig. 8)    on the cylinder or the punch. For every map, the shape of the damaged area remains smooth. These results can also be drawn along profiles for the C/P case (Fig. 9) and for the P/P case (Fig. 10). In  Figs. 9(a) and 10(a), the wear depth is plotted for all the points in the xy-plane as a function of their abscissa x. Both the "U-shape" and "W-shape" profiles are highlighted. This type of profile is classically observed in Ref. [3]. The reference curves obtained for the 2D plane stress calculations are also shown. The agreement with the case 0   is very good for the C/P case. For punch-plate, the wear depth is about 15% larger for the 2D plane stress calculation than that in 3D. This is consistent with the fact that there is a true 3D effect for the P/P case (damage heterogeneity in the footprint, with a lower wear depth in the middle). The 3D calculation may then differ from 2D, which is not the case for C/P conditions, where the gradient in z-direction is low. This matter is demonstrated by Figs. 9(b) and 10(b), which show wear depths for points in the symmetry plane 0 x  as a function of z. The left side 0 z  corresponds to the symmetry plane of the system, and 1 z  to the lateral free surface.
Even if in some cases (for small  ), wear depth is lower on the cylinder free surface, and profiles are almost constant in the direction z for the C/P case. In constrast, for the P/P case, wear depth is significantly higher on the lateral free surface, involving a strong gradient in z-direction. For a given value of  , the apparent scatter is linked to the fact that the points have a different location in the footprint (z values between 0 and 1). This scatter is more emphasized for the P/P case. This is consistent with the previous maps, which exhibit a large heterogeneity of the wear depth between the center of the footprint (   0, 0 x z ) and the corners of the punch ( 1 z  , x at the contact limit).
In conclusion, the analysis of the previous results shows that the footprints and wear profiles in different planes are strongly dependent on the WDP parameter. It is therefore essential to use a 3D two-surface wear model to best predict the erosion of the contact pair and capture the local phenomena.

Damage kinetics and evolution of the contact pressures
After having discussed the shape of the footprints for a given profile in Section 4.2, the present paragraph illustrates three different moments in the test, namely the first profile, the final profile in Section 4.2, and an intermediate profile (25 for C/P case, 5 for P/P). Figures 11 and 12 deal with the evolution of the wear depth in the xy-plane (the left column) and Fig. 11 Influence of the  (= 0, 0.5, and 1) during the fretting test on wear depth. The values above (below) the horizontal axis concern the cylinder (the plate) (the left column), and the contact pressure distribution in the xy-plane for cylinder-plate configuration are shown in the right column. The following points can be deduced from the cylinder-plate configuration. During the fretting test, the wear depth increases, and the contact area extends in the same way irrespectively of the value of  (the left column in Fig. 11). The two-surface wear model perfectly redistributes the total volume of wear between the contact surfaces. For 0 5    , the wear profile of the plate and the cylinder are almost identical. The point scatter enables to visualize the wear depths in the perpendicular direction to the sliding. Larger wear depths correspond to the 0 z  plane for the plate and to the 1 z  plane for the cylinder (Fig. 9).
As far as contact pressure is concerned, the profile reduces and flattens during the numerical fretting test. For the cylinder-plate case (the right column in Fig. 11), the initial pressure profile follows a Hertzian type distribution that evolves toward a distribution Fig. 12 Influence of the  (= 0, 0.5, and 1) during the fretting test on wear depth. The values above (below) the horizontal axis concern the punch (the plate) (the left column), and the contact pressure distribution for punch-plate configuration are shown in the right column.
www.Springer.com/journal/40544 | Friction predicted in Ref. [29]. The successive curves depart more slowly from this profile when the value of WDP is higher. Obtaining a flattened shape requires a higher number of cycles.
For the punch-plate configuration, the main trends in terms of wear and pressure are very similar for the various values of  . The set of points allows to visualize the wear depth (the pressure) in the direction perpendicular to the sliding. We note more dispersed results for this type of specimen. The evolution of the wear depth is in good agreement with the corresponding experiments. The maximal wear depth is located at the lateral borders, as observed for pressure (the left column in Fig. 12). More specifically, this maximum wear is reached on the 1 z  edge for the two contact solids. This position remains the same during the fretting test. The contact area cannot be extended due to the geometry of the impactor. The initial distribution of the contact pressure matches the expressions proposed by Ciavarella et al. [30]. The evolution of the pressure is weak during the cycles and is close for all  .
It is well known that the evaluation of contact pressure is not an easy task, and the field may be highly perturbated by the local geometrical defects. The results in Figs. 11 and 12 are raw values, and it may be interesting to discriminate between the effect of the perturbation and the real gradient of the field. An estimation of the difference is provided in the Appendix. For the sake of readability, only one profile is shown for each geometric case (namely 50 for C/P and 10 for P/P). In brief, the  does not influence the contact area but can play a significant role in the evolution of pressure according to contact geometries.

One-or two-surface wear model
The objective in Section 4.4 is to compare one-and two-surface wear models in terms of resulting wear depth, wear volume, contact width, and stress, in order to assess the benefits of the two-surface wear model. For the two-surface wear approach, an equivalent wear profile tot h  is determined by adding for each x-value of the computed wear depth  , or m s h    on master and on slave surfaces. It is compared with the h   at each point in the case of one-surface wear.
For a better clarity, we only present the equivalent wear profiles on the xy-plane for 0 z  and on the yz-plane for 0 x  (Fig. 13). The conclusions are very close for the examined geometries. The sum of the master and slave profiles closely correlates regardless of the WDP. Thus, the relations on the xy-and yz-planes are where L is the lateral width of the master, and a is the semi-contact width. These results, a direct consequence of using the Coulomb friction model and Archard's wear law, highlight the robustness of the numerical environment. Moreover, the contact widths are comparable (Fig. 14). Relative differences do not exceed 4.2% for the maximum wear depth and 1% for the contact width. As shown in Fig. 14(a) for C/P case and Fig. 14(b) for P/P case, the volume losses of master m V and slave s V linearly depend on the WDP. In fact, this point is in agreement with Refs. [13,31]. It only seems valid in the case of very small stroke amplitude, which is consistent with the assumptions of our study. For the case of 0 5    , the wear volume of the master and the slave are identical. If the wear volumes are summed up all together, the result is remarkably stable for all  values.
where tot V represents the total wear volume. The same observation can be made with a contact size, which is reported on the same plot (scale on the right).
Lastly, the analysis focuses on the stress state at the contact interface. For the two-surface approach, this stress state remains similar to that of the one-surface case. These observations are valid for the two studied geometries. To illustrate these comments, the comparison of some results is presented for the P/P case (Fig. 15). The normal zz  and the von Mises stresses mises  on the xy-plane are only displayed for 0 z  and 1 z  for better clarity (the left column). The evolution of these stresses during a fretting cycle is also plotted for a node at the contact limit, x a  (the right column). These obtained curves by the both wear modellings are overlaid in the direction perpendicular to the sliding direction z and all along the time.
These results confirm the validity of reducing the simulations to one-surface wear, provided that only . The left column shows the results for the C/P configuration after 50 profiles, and the right column shows the results for the P/P configuration after 10 profiles. the total depth or total volume of wear is targeted. Consequently, one-surface wear model is valid for a global approach of the phenomenon.

Discussion
The purpose of Section 4.5 is to compare our results to those in the literature. This comparison seems complex for two reasons. Firstly, to the best of authors' knowledge, there are just a few papers about two-surface wear [12][13][14][15][31][32][33][34]. Secondly, the selected geometries, loading conditions, wear models, or resolution techniques usually differ from one article to the next.
Our C/P configuration is found in two other articles [12,13]. Nevertheless, these studies concern 2D element finite simulations, and the conditions are not quite the same. Then, comparisons can only be made on a qualitative basis, as shown in Table 1 in the second one. The fretting test is achieved under an amplitude   lower than the a in the first research, while this amplitude is equal to or greater than the a in the second one. Figure 16 illustrates the equivalent wear profiles from these 2D studies and from our 3D case in the xy-plane, 0 z  . The results in the literature show that the relative difference of the maximum equivalent wear is small (about 1%), and the contact widths are similar. According to the authors, the simulation of one-surface wear can then be considered as a successful strategy to reduce the computation time. It should be stressed that this conclusion is only valid when the amplitude is smaller or equal to the a. As for our results, we get the same trends as those resulting from the literature, highlighting a good agreement between one-and two-surface wear. Thus, Fig. 15 Example of the stress state during a fretting cycle in the xy-plane (the left column) and for a border contact node (the right column) for punch-plate configuration. if we stick to a quantitative analysis (wear volume and wear depth), we can simply settle for a one-surface wear model.

Conclusions
The purpose of this study is the implementation of a consistent and robust numerical environment to predict 3D wear evolution on each contact surface, taking into account the difference of resistance between the bodies.
Based on average friction energy formulation, this new model gives the ability to predict erosion on two contact surfaces through a WDP. The framework focuses on describing both local and global data, so that a qualitative and quantitative analysis can be performed. Numerical simulations are carried out for 3D Ti-6Al-4V-Ti-6Al-4V contact specimens in sliding conditions.
The further analysis reveals some salient points, whose main contents are listed below: 1) Numerical tangential force-displacement loops are similar whatever the WDP, and their areas increase during the fretting test.
2) The evolution of contact width 2a is unchanged during the fretting cycles whatever the value of WDP.
3) The wear volume of each contact bodies depends on the WDP, but the sum of these volumes is kept constant. Finally, the one-surface wear model is sufficient as long as the global quantities are regarded, or the difference in wear resistance between the bodies is great. It should be noted that the wear volume is well balanced for 0 5 4) The footprints are specific to the difference between the resistance of the bodies, and their shapes differ according to the situation. Thus, the two-surface wear model enables to precisely describe the wear footprint observed in the experimental tests. Accordingly, this wear model is essential when a detailed description of the mechanisms is desired. 5) On the xy-and yz-planes, the wear profile shape is specific to the WDP. However, the equivalent wear profile, i.e., their sum of the master and slave profiles, is identical.
6) WDP can play a fundamental role in the evolution of the contact pressure according to the geometric  specimen. An initial profile of the Hertzian type is systematically obtained, but reduced more or less quickly for C/P configuration, whereas no change is observed in P/P configuration. These various remarks underline the fundamental role of WDP on the qualitative analysis of the data. Depiste the gain in computing time, a one-surface wear model is no longer sufficient to best capture the local mechanisms, and does not take the body resistance into consideration. Regarding the local quantities, the 3D two-surface wear becomes essential. Moreover, the result comparison with the experimental studies in the literature highlights the robustness of our numerical environment. For C/P configuration, we observe an evolution of the wear depth and an extension of the contact area during the fretting test (the left column in Fig. 17). The profile of contact pressure reduces and flattens more or less quickly depending on the value of  . The  only plays a function on the contact pressure distribution (the right column in Fig. 17).

Appendix Complementary results
In contrast, for P/P configuration, the maximum of the wear peak is always located at the same position during the fretting test (at the contact border), and there is no extension of the contact area (Fig. 18). At last, the contact pressure changes slightly. Figure 19 completes the analysis of the contact pressure and highlights the impact of geometric perturbations.

Declaration of competing interest
The authors have no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which concern the cylinder (the plate) (the left column), and the contact pressure distribution in the xy-plane for cylinder-plate configuration are shown in the right column.

Fig. 18
Effect of the  (=0.25 and 0.75) during the fretting test on wear profiles. The values above (below) the horizontal axis concern the punch (the plate) (the left column), and the contact pressure distribution for punch-plate configuration are shown in the right column. Fig. 19 Illustration of 3D effect on contact pressure for profile 50 (left column, C/P case) and profile 10 (right column, P/P case). The pressure values are fitted for five classes depending on the z-coordinate.
www.Springer.com/journal/40544 | Friction permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.
The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.