A Simple Semi-Analytic Contact Mechanical Model for Tangential and Torsional Fretting Wear of Axisymmetric Contacts

: Fretting wear of axisymmetric contacts is considered within the framework of the Hertz– Mindlin approximation and the Archard law for the linear wear. If the characteristic time scale for the wear is much larger than the duration of a single fretting oscillation, the proﬁle change due to wear during one fretting cycle can be neglected for the contact problem as a zero-order approximation. This allows to give an exact contact solution during each fretting cycle, depending on the current worn proﬁle, and thus for the explicit statement of an ordinary integro-differential equation system for the time-evolution of the fretting proﬁle, which can be easily solved numerically. The proposed method gives the same results as a known, contact mechanically more rigorous simulation procedure that also operates within the framework of the Hertz–Mindlin approximation, but works signiﬁcantly faster than the latter one. Tangential and torsional fretting wear are considered in detail. A comparison of the numerical prediction for the evolution of the worn proﬁle in partial slip torsional fretting of a rubber ball on abrasive paper shows good agreement with experimental results from the literature.


Introduction
Many tribological contacts in technical or biological systems are subject to oscillations because they are part of a periodically operating machine or a periodically changing environment. The displacement amplitudes are often too small to cause global sliding, i.e., the contacts consist of changing areas of local sticking and slipping/sliding. This phenomenon, called "fretting", causes wear as well as fatigue and, although highly localized, is of enormous importance for the service life of various tribological systems, e.g., electrical conductors [1], turbine blades [2], or artificial joints [3].
Depending on the maximum extent of the slip area during the fretting oscillation, often different fretting regimes are distinguished [4], namely, the partial slip regime (for small slip areas), the gross sliding regime, and a mixed regime in between. While in different regimes, either wear or fatigue may be the dominating damage mechanism, both always are interacting or even competing phenomena [5]. Therefore, the theoretical solution procedure for a given fretting problem should consist of at least three components: a contact mechanical solution (depending on the material model, the frictional regime, and so on) as the basis for the determination of local stresses and displacements [6,7], a wear routine [8][9][10], as well as a crack initiation and propagation routine [11][12][13]. Doing this in a coupled and rigorous way (e.g., based on FE simulations [14]) is often very costly from a computational perspective. Here, semi-analytic (contact mechanical) approaches-which already have various applications in tribology [15]-can be a good "trade-off" between effort and prediction quality [16].
Therefore, in the present manuscript, a simple and efficient semi-analytic contact mechanical procedure for the fast simulation of fretting wear is proposed and discussed, which could also be extended for the analysis of fretting fatigue. The model operates within the framework of the Hertz-Mindlin approximation to general axisymmetric tangential contact problems, which will be stated in the manuscript, and which allows for an exact solution, if the contact problem only slowly changes with time.
The manuscript is organized as follows. In Section 2, first, the physical assumptions of the model are stated and the necessary classic solutions of axisymmetric contact mechanics are re-iterated. These solutions are formulated in terms of Abel-type integral transforms, the efficient numerical implementation of which is also briefly summarized. In Section 3, the semi-analytical procedure for the fast simulation of tangential and torsional fretting wear is presented and its predictions are compared to a contact mechanically more rigorous model and to experimental results from the literature. A discussion and conclusive remarks finish the manuscript.

Fundamental Assumptions
Firstly, we assume that all deformations are elastic. It follows from the fundamental solution by Boussinesq [17] and Cerruti [18] (see, e.g., the book of Johnson [19]) for the static relative surface displacements between two elastic half-spaces with the shear moduli G 1 and G 2 and the Poisson ratios ν 1 and ν 2 , which act on each other with a point force at the origin of a common cartesian coordinate system {x,y,z}-z being the direction of the common normal of the half-spaces-that the normal and tangential contact problems are elastically decoupled, if the half-spaces are elastically similar, i.e., if For half-space contact mechanics to be valid, the contacting bodies must obey the restrictions of the half-space approximation, i.e., their macroscopic dimensions have to be much bigger than the largest characteristic contact size, and the surface gradients in the vicinity of the contact must be small. Moreover, surface roughness is neglected, and all timedependent processes are analyzed in the quasi-static limit, i.e., elastic wave propagation is not considered. Moreover, friction in the contact shall obey a local Amontons-Coulomb law with a constant coefficient of friction µ, and for the uniaxial tangential contact problem with friction, the lateral displacements u y resulting from a surface shear stress τ xz are neglected. The entirety of the latter assumptions is what we will refer to as Hertz-Mindlin approximation [20,21] to the axisymmetric tangential contact problem with friction.
For the determination of the wear dynamics, a local Archard [22] wear law is used. Note that the Archard law, in connection with the friction law assumed by the contact mechanical component of the model, is equivalent to an energy-based wear model, which is commonly used in the literature on fretting wear [23]. However, from the point of view of the semi-analytical contact model given in Section 3, any local wear law according to which the profile change is an explicit function of the local contact pressure and the local slip displacement (or other known quantities) could be applied.
The severity of all these assumptions will be discussed later in the manuscript.

Elementary Axisymmetric Contact Solutions
Let the gap between the contacting bodies at first contact be an axisymmetric, monotonous, continuous function f (r), with r being the polar radius in the contact plane. Then, the normal contact solution (with prescribed contact radius a) for the indentation depth d, normal force F, contact pressure distribution p(r), and relative displacements u z (r) outside the contact area, first found by Schubert [24] and Galin [25], is given by [26]  with the effective Young's modulus Usually, the prescribed quantities for the contact problem will rather be either the indentation depth (for displacement-controlled systems) or the normal force (for force-controlled systems), instead of the contact radius; however, the latter can easily be retained from the respective prescribed quantity by inverting the relevant equation in the above solution. Now, consider the simple Cattaneo-Mindlin [21,27] loading history for the contacting bodies: they are first pressed against each other by a normal force F and subsequently loaded in the tangential direction by an increasing force F x . As the tangential load increases, local slip propagates from the contact edge into the contact area. The radius of the remaining inner stick area at the end of the loading procedure shall be c. Ciavarella [28] and Jäger [29] were the first to recognize that the general axisymmetric Cattaneo-Mindlin problem (within the simplifications detailed in the previous subsection) is solved by a basis (indicated by the index "B") tangential stress distribution i.e., the problem can be reduced to the frictionless normal contact problem. The slip displacement distribution in the slip area after the complete loading process is given by [26] with the effective shear modulus and the same notations as in Equation (2). The argument c in Equations (4) and (5) refers to the normal contact solution if the contact radius were c. The basis relations between c; the macroscopic relative tangential displacement, u x0 ; and the tangential force, F x , are [26] The solution to the torsional contact problem is a little more elaborate numerically, because there is no equivalent to the Ciavarella-Jäger principle for tangential contact, i.e., the torsional contact problem with friction cannot be reduced to the frictionless normal contact. Nevertheless, conceptionally, the solution to the general axisymmetric torsional contact problem with friction poses no difficulties and has been found by Jäger [30]. The assumptions for that solution are comparable to the Hertz-Mindlin approximation, namely, the contact partners are assumed to be elastically similar, and the frictional interaction shall obey a local Amontons-Coulomb law.
For the elementary loading history (constant normal force, subsequently applied increasing torque; sometimes referred to as the "Lubkin [31] problem"), the solution goes as follows [32]: first, a function φ must be determined according to which gives the relation between the amplitude of the torsional oscillation, ϕ A , and the stick radius: Symmetry 2021, 13, 1582

of 11
After that, the slip displacement distribution in the slip area is given by The torque as a function of the stick radius can be calculated easily, but will not be needed in later stages of the present manuscript.

Numerical Implementation of the Abel-Like Integral Transforms
The efficient numerical implementation of the integral transforms appearing in the previous subsection is straight forward, but shall be described briefly in the following.
Let us start with the determination of the equivalent profile g(x) from the axisymmetric profile f (r): Following an idea by Benad [33], we integrate by parts to get rid of the singularity in the integrand: where the prime denotes the first derivative with respect to the given function argument.
The last integral is evaluated on equidistant arrays for the spatial coordinates. Applying the trapezoidal rule, the integral is reduced to a simple matrix-vector multiplication. All derivatives are calculated by second-order finite differences. Integration by parts for the pressure distribution p(r) results in where cosh −1 denotes the inverse hyperbolic cosine. For the additional transformations in the torsional case, integration by parts gives and For brevity, the discretization procedure is omitted, as it works completely similar to as in the first case.

Results
In this section, the proposed semi-analytical procedure for the simulation of fretting wear is shown. First, the model for tangential fretting wear is detailed and its predictions are compared to a contact mechanically more rigorous numerical model in Sections 3.1 and 3.2. In Section 3.3., the main idea of the model is applied to a different fretting mode, torsional fretting, and the results of a numerical simulation for the evolution of the worn profile in a partial slip torsional fretting contact of a rubber ball on abrasive paper are compared to experiments from the literature in Section 3.4.

Model for Tangential Fretting Wear
The Cattaneo-Mindlin loading only corresponds to the first loading of the fretting contact. As is known, the solution to the tangential contact problem depends on the loading history. So, how will the contact solution look like during the fretting cycles? When the direction of the tangential motion is reversed, there is a moment of complete stick, because, according to the local Amontons law, the whole slip area is at the boundary state |τ| = µp. Moreover, wear happens very slowly, i.e., on a much larger time scale than the fretting cycle duration. If there were no wear at all during two cycles, after reversing motion, a basis stress distribution can be linearly superposed [34]: to obtain the stress distribution with the current stick radius c*. At the end of the half-cycle, the stress distribution is the same as at its beginning, with the sign reversed. Therefore, on the one hand, the relations between the amplitudes of the tangential force or displacement oscillation, F A and u A , and the minimum stick radius c of the fretting cycle are also given by Equation (7). On the other hand, the slip displacement during the half-cycle is During the second half-cycle, the same happens again. Hence, assuming a local Archard wear law, a zero-order approximation for the profile change during one cycle is given by where the time t is measured in cycle durations, k wear is the dimensionless wear coefficient in the Archard law and σ 0 is the hardness of the worn material. "Zero-order approximation" in this context refers to the fact that the influence of the change of profile due to the wear in one cycle on the contact problem and the resulting change in the slip displacement is neglected.
Equations (17) and (18) are based on the supposition that there is a remaining stick area when reversing the tangential motion. At the gross slip transition, one obtains If the sliding distances remain small compared with the contact radius (so that the system retains its approximate axial symmetry), the model can also be used in the gross slip regime. In the first fretting half-cycle, gross slip will start at a macroscopic relative tangential displacement Hence, the local slip displacement distribution in the contact area after a full cycle is given by Moreover, it is apparent from Equation (7) that, for both fully displacement-controlled (i.e., with a fixed indentation depth and amplitude of the tangential displacement oscillation) and fully force-controlled (i.e., with a fixed normal force and amplitude of the tangential force oscillation) tangential fretting, the minimum stick radius during one fretting cycle is constant over time, if there are no profile changes for r < c. That is selfconsistent within the basic assumptions of the model (no wear without relative motion between contacting surface points), and c thus has the meaning of a permanent stick radius.
Hence, the proposed semi-analytical model works as follows: in every fretting cycle, the slip displacement distribution depending on the current profile and the current pressure distribution is calculated; based on the Archard law, the profile is updated and the new pressure distribution is calculated based on the normal contact solution. Finally, it should be noted that the second and third term in the slip displacement formula in Equation (5) only depend on the function g within the permanent stick area, which does not change over time.
Very common characteristics of fretting contacts are "fretting loops", i.e., hysteresis loops between the macroscopic relative tangential displacement and tangential force during the fretting process. The shape of the fretting loop and its time evolution are often used to characterize the fretting regime [35,36], because it is a simple and effective way to visualize the frictional energy dissipation in the slip area of the contact. The fretting loops can be calculated very easily in the model described above. From the basic contact solution in Equation (7) and the superposition in Equation (16), it is clear that the relative tangential displacement and the tangential force during a specific fretting cycle (without gross slip) with the current equivalent profile g(x,t) are given by where the different signs correspond to the first and second fretting half-cycles. For brevity, the time argument was omitted. In the gross slip regime, c equals zero and the tangential displacement during the partial slip stage of the oscillation must be modified according to The force relation can be taken from Equation (22) and the relations in the sliding stage are trivial. Note that, because of the Ciavarella-Jäger principle, the fretting loop can be determined simply based on the normal contact solution for the current fretting cycle and the amplitude of the tangential oscillation.

Comparison with a Contact Mechanically More Rigorous Model
Dimaki et al. [37] proposed a fast numerical algorithm for fretting wear of axisymmetric contacts, where the slip displacement during one fretting cycle is calculated based on a separate (in the sense that it is almost independent of the wear simulation) numerical model within the framework of the method of dimensionality reduction (MDR) [38]; the fretting profile in their algorithm is also updated only once at the beginning of each fretting cycle and the MDR solution to the tangential contact problem operates within the Hertz-Mindlin approximation. Therefore, as the main (and practically only) difference between this numerical model and the semi-analytical approach proposed in the present manuscript is the calculation of the slip displacement, it makes sense to compare the predictions of both models regarding the time-evolution of the worn profile. In Figure 1, the results of both algorithms for the worn profiles are shown for an initially parabolic profile f (r) = r 2 /(2R), with some radius of curvature R, in normalized variables (note that, in normalized variables, the results only depend on the parameters given in the figure). Both algorithms were implemented in MATLAB ® R2021a. The radial coordinate is normalized to the initial contact radius, a 0 = (dR) 1/2 , and the worn profile for the fixed indentation depth d. The radius of the permanent stick area was set to c = 0.5a 0 . The characteristic number of fretting cycles for the wear problem, N 0 , results from the wear law [37]: with the initial normal force F 0 at the beginning of the fretting, and the tangential displacement oscillation amplitude u A . proach and the contact mechanically more rigorous model. Hence, it can be concluded that the zero-order approximation in Equation (18) gives a very good result (within the Hertz-Mindlin framework, of course) if the profiles are worn slowly compared with the fretting cycle duration. Nevertheless, if the implementation of the Abel-like transforms is done the same way, the semi-analytical model works faster by a factor of the order of one, because the number of numerical operations for one fretting cycle is significantly reduced.

Model for Torsional Fretting Wear
The main idea of the semi-analytical approach described in Section 3.1-that the slip displacement during one fretting cycle can be calculated analytically as a zero-order approximation-is also applicable to other fretting modes, for example, torsional fretting.
If, once again, we neglect the profile change on the time scale of a single fretting cycle, the torsional stress distribution between two reversal points of the fretting oscillation can be obtained by linearly superposing a stress distribution for the elementary loading in a similar way as in Equation (16) [40]. During one fretting cycle, this happens twice, so the total slip displacement during one fretting cycle is four times the expression given in Equation (10). The slip displacement once again determines the profile change according to the linear Archard wear law. Note that, for the torsional problem, c = 0 only exists as a limit, as the center point of the contact always remains sticking. There are no visible differences between the prediction of the semi-analytical approach and the contact mechanically more rigorous model. Hence, it can be concluded that the zero-order approximation in Equation (18) gives a very good result (within the Hertz-Mindlin framework, of course) if the profiles are worn slowly compared with the fretting cycle duration. Nevertheless, if the implementation of the Abel-like transforms is done the same way, the semi-analytical model works faster by a factor of the order of one, because the number of numerical operations for one fretting cycle is significantly reduced.

Model for Torsional Fretting Wear
The main idea of the semi-analytical approach described in Section 3.1-that the slip displacement during one fretting cycle can be calculated analytically as a zero-order approximation-is also applicable to other fretting modes, for example, torsional fretting.
If, once again, we neglect the profile change on the time scale of a single fretting cycle, the torsional stress distribution between two reversal points of the fretting oscillation can be obtained by linearly superposing a stress distribution for the elementary loading in a similar way as in Equation (16) [40]. During one fretting cycle, this happens twice, so the total slip displacement during one fretting cycle is four times the expression given in Equation (10). The slip displacement once again determines the profile change according to the linear Archard wear law. Note that, for the torsional problem, c = 0 only exists as a limit, as the center point of the contact always remains sticking.
Hence, the semi-analytical algorithm for torsional fretting wear works the same as in the case of tangential fretting, except that now, four Abel-like transforms must be performed in every cycle. Interestingly, as the torsional contact problem cannot be reduced to the normal one, the minimal stick radius during one cycle is not constant while the profile is worn, but rather slightly increases. Nevertheless, the notation c shall still refer to the radius of the permanent stick region, i.e., the minimal stick radius during the first fretting cycle. The characteristic number of fretting cycles is given by with the amplitude of the torsional angular oscillation, ϕ A.

Comparison with Experimental Results
In their analysis of the limiting no-wear profile for multi-mode fretting, Dmitriev et al. [41] conducted experiments on torsional fretting of a rubber half-sphere on a flat covered by abrasive paper. The radius of the sphere is (approximately) 19 mm (the paper wrongly says 14.5 mm, but from the respective figure, the correct value can be estimated), the indentation depth is 2 mm, and the radius of the permanent stick area is 3.25 mm [41]. The theoretical prediction will only depend on the number of normalized fretting cycles, N/N 0 . In the experiments, the profiles were measured in consecutive steps of 2.160 cycles each. Based on the experiments, the value of N 0 was estimated as 14.400, which corresponds to a dimensional wear rate of approximately 7.5 mm 3 /m (which is quite high). In Figure 2, the experimental results and the theoretical prediction are shown. They agree well with each other; in the experiments, the wear seems to be slower than predicted in the beginning (which makes sense, as the wear rate in the simulation is high); later, the wear rate in the experiments seems to increase and the convergence to the limiting profile is even faster than in the numerical calculation. Moreover, for the elastic ball in the experiment, there is a visible displacement of the lower point of the wearing area (corresponding to the edge of the permanent stick zone), while the position of the upper point (corresponding to the edge of the initial contact area) is virtually unchanged. On the other hand, the inverse relationship is observed in the calculation by the semi-analytical model. fretting cycle. The characteristic number of fretting cycles is given by with the amplitude of the torsional angular oscillation, φA.

Comparison with Experimental Results
In their analysis of the limiting no-wear profile for multi-mode fretting, Dmitriev et al. [41] conducted experiments on torsional fretting of a rubber half-sphere on a flat covered by abrasive paper. The radius of the sphere is (approximately) 19 mm (the paper wrongly says 14.5 mm, but from the respective figure, the correct value can be estimated), the indentation depth is 2 mm, and the radius of the permanent stick area is 3.25 mm [41]. The theoretical prediction will only depend on the number of normalized fretting cycles, N/N0. In the experiments, the profiles were measured in consecutive steps of 2.160 cycles each. Based on the experiments, the value of N0 was estimated as 14.400, which corresponds to a dimensional wear rate of approximately 7.5 mm 3 /m (which is quite high). In Figure 2, the experimental results and the theoretical prediction are shown. They agree well with each other; in the experiments, the wear seems to be slower than predicted in the beginning (which makes sense, as the wear rate in the simulation is high); later, the wear rate in the experiments seems to increase and the convergence to the limiting profile is even faster than in the numerical calculation. Moreover, for the elastic ball in the experiment, there is a visible displacement of the lower point of the wearing area (corresponding to the edge of the permanent stick zone), while the position of the upper point (corresponding to the edge of the initial contact area) is virtually unchanged. On the other hand, the inverse relationship is observed in the calculation by the semi-analytical model.

Discussion
As stressed, the presented semi-analytical model of fretting wear is based on several simplifying assumptions that restrict its scope of applicability, at least in a quantitative sense. The most severe simplification is probably the local application of the linear Archard wear law. Moreover, neglecting surface roughness and non-elastic deformations can be unjustified simplifications in several circumstances.
Fouvry et al. [23] have pointed out that an elementary global-local, e.g., energy-based, description of the fretting wear process is applicable if the tribo-couple is dominated by non-adhesive wear mechanisms. This is confirmed by the above comparison of the present model proposal with the experimental results on torsional fretting wear of rubber on abrasive paper and in line with the findings of Varenberg et al. [42] that the role of oxide wear debris depends on the dominant fretting wear mechanism. However, in general, and especially in systems dominated by adhesive wear, an appropriate description of the third body in contact (i.e., the debris behavior) is vital for a comprehensive understand-9 of 11 ing of fretting wear. A deep overview on predictive wear models was given by Meng and Ludema [43]. For an up-to-date review on recent advances in the understanding of wear mechanisms and modelling, the reader is referred to the extensive work of Meng et al. [44]. However, quantitatively predictive wear equations are still extremely scarce and usually the main bottleneck for the predictive power of (even very elaborate) fretting wear simulations-which are often also based on a simple Archard-or energy-based wear law, despite their highly sophisticated material and contact models.
Another problematic characteristic of the model is that the macro-slip slope of the fretting hysteresis loop is zero, because it strictly is a Coulomb-based model, i.e., the friction coefficient is assumed to be constant (at least over one fretting cycle); however, in experimental fretting loops, the macro-slip slope is usually slightly positive (see, e.g., [45]).
Hu et al. [46] analyzed the effect of plastic deformations (due to the forming stress singularity at the boundary of the stick zone) on the fretting behavior and found that the main effect of plasticity is to allow the wear to penetrate the permanent stick region, and thus to never cease.
Moreover, the other assumptions of what in the present manuscript is referred to as the Hertz-Mindlin approximation can be a source of quantitative error. Generally, the tangential contact problem is not axisymmetric, and neither is, for example, the shape of the stick area. However, the error of the approximation is relatively small (at least compared with the uncertainty in the prediction of wear coefficients), as demonstrated by Munisamy et al. [47] by a self-consistent numerical solution of the Cattaneo-Mindlin problem. A deep and general analysis of tangential contact problems with elastic coupling can be found in the book of Barber [48].
On the other hand, several extensions of the proposed method are possible and planned for future work. Firstly, the Abel-like transforms can be written as explicit convolutions and, therefore, accelerated by the fast Fourier transform (FFT) [49]. Moreover, once the normal and tangential contact problems in the Hertz-Mindlin approximation are solved, the complete stress state inside the materials can be obtained via simple 1D-integrals [50], which allows for the fatigue analysis. Hence, the suggested semi-analytical model can be a part of a more global model of fretting in axisymmetric contacts. Finally, based on the tangential contact solution for power-law graded elastic materials by Heß and Popov [51], the proposed semi-analytical procedure can be generalized for power-law graded materials. This might be of high practical relevance, as this material class possibly offers a solution to the wear-fatigue dilemma in fretting [52].

Conclusions
A semi-analytic contact mechanical model for the fast simulation of fretting wear of axisymmetric contacts has been proposed. The main idea of the model is to neglect the change of the indenter profile during one cycle of the fretting oscillation, which-within the Hertz-Mindlin approximation-allows to analytically calculate the slip-length during one cycle in closed form, and hence gives a zero-order approximation for the time-evolution of the worn profile. This procedure gives the same results as contact mechanically rigorous simulations based on the Hertz-Mindlin approximation, but works significantly faster than those. A comparison of the numerical prediction for the evolution of the worn profile in partial slip torsional fretting of a rubber ball on abrasive paper shows good agreement with the experimental results from the literature.