Numerical Simulation of Slip-Stick Elastic Contact

Fretting defines a condition in which mechanical contacts are subjected to alternating tangential displacements, small compared to dimensions of contact area, due to oscillating loading conditions. Fretting wear and fretting fatigue are between the most important factors responsible for contact failure, especially when high loads are transmitted through non-conforming contacts, leading to highly localized stress concentrators in the vicinity of the contact region. Prediction of life span of machine elements working in such conditions requires assessment of stress and strain in the contacting bodies, which is the main subject of Contact Mechanics. Although fretting is intrinsically a multidisciplinary process, involving adhesion, oxidation, abrasion and pitting, modern approach suggests that contact stresses play a chief role.


Introduction
Fretting defines a condition in which mechanical contacts are subjected to alternating tangential displacements, small compared to dimensions of contact area, due to oscillating loading conditions.Fretting wear and fretting fatigue are between the most important factors responsible for contact failure, especially when high loads are transmitted through non-conforming contacts, leading to highly localized stress concentrators in the vicinity of the contact region.Prediction of life span of machine elements working in such conditions requires assessment of stress and strain in the contacting bodies, which is the main subject of Contact Mechanics.Although fretting is intrinsically a multidisciplinary process, involving adhesion, oxidation, abrasion and pitting, modern approach suggests that contact stresses play a chief role.
While analytic solutions in this research field lead to complex mathematical models, many without closed-form solution, numerical approach reveals itself as a useful engineering tool, capable of extending the few existing analytical results to technologically important contact scenarios.A numerical study may advance the understanding of fretting contact and provide assistance to the design of contacts with improved load-carrying capacity.
Elastic contact analysis considering interfacial friction and slip-stick behaviour originated in the works of Cattaneo (Cattaneo, 1938) and Mindlin (Mindlin, 1949).They proved independently that, even when the contacting bodies are globally sticking, a peripheral region of slip is to be assumed in order to remain in the frame of Linear Theory of Elasticity and to obey the Coulomb's law of friction.Based on these results, Johnson (Johnson, 1985) advanced the closed-form solution for the contact between similarly elastic materials undergoing a fretting loop.
In case of dissimilarly elastic materials, when the effects of normal and tangential tractions are coupled, an iterative solution has been achieved (Hills et al., 1993) only for the plane (i.e.cylindrical) contact.Many authors employ the so-called Goodman approximation (Goodman, 1962) when dealing with this type of contact, which neglects the influence of shear tractions on pressure, but retains that of pressure on tangential tractions.As proved in (Hills et al., 1993), this approximation is satisfactory in case of plane (cylindrical) contacts if Poisson's ratio is large enough, but the inaccuracy introduced in the simulation of the three dimensional contact between dissimilarly elastic materials cannot be a priori assessed.
In order to overcome this obstacle, recent works aimed to solve the problem numerically, using a method derived from the boundary element method, also referred to as semianalytical (SAM) in a review paper by Renauf et al. (Renauf et al., 2011).The strong point of this technique is that only a small region of the boundary of the contacting bodies, enclosing the contact area, is to be meshed, leading to a dramatic decrease in computational complexity compared to finite element method, in which discretization of the entire bulk is required.
Chen and Wang (Chen & Wang, 2008) advanced an algorithm for the non-conforming contact of dissimilarly elastic materials, and predicted the additional effect of an increasing tangential loading.Wang, Meng, Xiao, and Wang (Wang et al., 2011) investigated numerically the supplementary effect of a torsional moment, while Wang et al. (Wang et al., 2010) applied the algorithm advanced in (Chen & Wang, 2008) to contact of elastic layered half-spaces.However, the loading history was not accounted for in these studies, i.e. the full load was applied in one step.Gallego, Nélias, and Deyber (Gallego et al., 2010) applied numerical analysis in an incremental approach to study different fretting modes, and concluded that assumptions adopted in existing analytical models lead to arguably inaccurate results.It is asserted in (Gallego et al., 2010) that, due to irreversibility of friction, which is a dissipative process, loading history should be considered although a purely elastic contact analysis is intended.
An incremental iterative algorithm for the fully coupled elastic contact with slip and stick is advanced in this work.Existing algorithms for the uncoupled normal or tangential contact problems are adapted for modeling of transient contact, and combined in an iterative approach based on the mutual adjustment between contact tractions, resulting in a three level nested loop algorithm.The use of modern numerical methods allows for a fine discretization in both spatial and temporal domain, leading to well converged numerical solutions.

Formulation of continuous slip-stick elastic contact problem
In contact problem formulation, it is convenient to describe the initial geometries ( ) 1 2 ( , ) i hi x x of the two contacting bodies 1,2 i  in a Cartezian coordinate systems having its 1 x and 2 x axes contained in the common plane of contact (i.e. the plane passing through the first point of contact, chosen as to separate best the bounding surfaces).The direction of 3 x -axis will be referred to as the normal direction, while the other two are tangential.In the three dimensional case, forces and moments transmitted through the contact have components along all three axes, namely the normal force W , the tangential force , M M and the torsional moment 3 M .Superscripts denote the contacting body, and subscripts are used for the direction of the referred quantity.Under load, the bodies deform unless assumed rigid, leading to elastic displacements ( ) j i u , and move as rigid-bodies with translations ( ) j i  and rotations ( ) In Contact Mechanics, it is also assumed that contact area dimensions are small compared to extents of the contacting bodies, and therefore stresses in the contact region are independent of other boundary conditions.This assumption is well verified in case of non-conforming contacts, when stresses induced by the contact process are highly localized in the vicinity of the contact region.
Once a contact area is established, the imposed forces and moments lead to contact tractions, i.e. pressure ( ) j p in the normal direction and shear traction ( ) 1 2 ( , ) j q q q in the tangential direction.The latter appears only if interfacial friction is assumed, leading to three possible cases, in relation to the magnitude of the tangential load: full stick, partial slip (or slip-stick), or gross slip.The latter case is trivial, as shear tractions are related to pressure through Coulomb's law on all contact area.On the other hand, the works of Cattaneo (Cattaneo, 1938) and Mindlin (Mindlin, 1949) prove that the full-sticking contact cannot be solved in the Frame of Linear Theory of Elasticity, as it leads to infinite stresses at the boundary of the contact area.The study of the partial slip contact, which is found in fretting contact processes, concluding with assessment of contact tractions, is the main goal of this work.
The temporal dimension t is included in this model along the spatial dimensions 1 2 , x x to provide basis for reproduction of the loading history.Consequently, ( ) C t  denotes the contact area, and 1 2 ( , , ) h x x t the surface separation, both established at a specified point t on the loading curve.When 0 . The other relative (composite) quantities, lacking the superscript indicating the contacting body, are defined in a similar way: (1) (2) , and (1) Computation of the relative normal displacement 3 1 2 ( , , ) u x x t will be discussed in Section 3.2.
The complementarity conditions in Eq. ( 4) show that only compressive normal traction (i.e.pressure) is allowed on the contact area, meaning adhesion is not accounted for.While adhesion cannot be ruled out in case of rubber, the metallic materials are found to show little adhesion effects, as the actual contact area, established between the peaks of the inherent surface microtopography (i.e.roughness), is much smaller than the theoretical one.Therefore, study of adhesion effects is beyond the point of this study.
The framework leading to Eq. ( 3), discussed in detail in (Johnson, 1985), is depicted in Fig. 1.The tilting angles, resulting from application of flexing moments, are omitted for brevity.The dashed lines show the initial (i.e. at 0 t  , in undeformed state) profile of the contacting bodies, but in positions (relative to the initial point of contact O ) corresponding to rigidbody translations (1)  3 ( ) t  and (2) 3 ( ) t  , which have opposite signs due to the fact that the contacting bodies are compressed.In order to accommodate the interpenetration distance 3 ( ) t  (it should be remembered that the bodies are assumed impenetrable in the frame of Linear Theory of Elasticity), the bodies deform elastically, resulting in normal displacements (1) 3 ( ) u t and (2)  3 ( ) u t pointing inward the corresponding body, as both normal contact tractions are compressive.Superposition of these processes yield the profiles on the contacting bodies in deformed state (at a specified time t ), depicted using continuous lines in Fig. 1.

The contact model in the tangential direction
An analogous model can be established for the contact in the tangential direction, when a slip-stick regime is assumed.The model for this contact process, also developed in (Johnson, 1985;Chen & Wang, 2008), consists in the following equations and inequalities: 1.The static force equilibrium (at any time t ): , ) .
C t M t q x x t x q x x t x dx dx 2. The geometrical condition of deformation in the time frame 1 2 [ , ] t t , in which the tangential load is not allowed to change sign: 3. The contact complementarity conditions: Here, S  is the stick area, C S    the slip region,  the frictional coefficient and 1 2 ( , ) s s s the relative slip distances.The base for Eq. ( 7) is presented in Fig. 2, which depicts the slipstick contact process in the direction of 1 x  .Rigid-body motion and elastic deformation due to torsion are omitted for brevity.The time parameter is also omitted, meaning all quantities are bound to the considered time frame 1 2 [ , ] t t .Let us consider two points 1 P and 2 P on the contacting bodies (1) and (2), respectively, located at 1 t on the same axis normal to the common plane of contact.The depths of these points are large enough as to assume that the corresponding tangential deformations can be neglected.The axis intersects the bounding surfaces in the points 1 A and 2 A , respectively, where the bodies deform due shear tractions.
Firstly, the above mentioned axis is assumed to pass through the initial point of contact, therefore A also undergo tangential displacements.If for these points the composite parameters cancel each other, i.e. 1 2  and consequently a stick regime is established in O .
Secondly, the 1 2 P P axis intersects the contact area in the peripheral region, where, although . This position corresponds to a region in relative slip (also referred to as microslip).It should be noted that any point in the current contact area is either in stick, where the norm of the shear traction is smaller than the limiting friction, or in slip, where the shear traction norm equals the limiting friction.The existence of slip is intrinsically conditioned by an increase or decrease in the level of tangential load, and therefore a purely static model is not appropriate.

Formulation of discrete slip-stick elastic contact problem
The main obstacles in solving the continuous contact problem consist in the fact that neither the contact and stick area, nor the tractions distributions are known in advance.A trial-anderror approach is therefore necessary, but this infer integration of arbitrary functions (contact tractions distributions) over irregular domains (contact or stick area) in displacements computation.As this cannot generally be performed analytically, numerical integration, although computationally intensive, is preferred.The basic principles of contact problem discretization are discussed in this section, and the advantages of this approach are outlined.Computation of displacements fields using fast numerical methods derived from the theory of Digital Signal Processing (DSP) is also discussed.

Principles of problem discretization
Numerical resolution of elastic contact problem relies on considering continuous distributions as piecewise constant on the elementary cells of a mesh established in the common plane of contact, enclosing the contact area at any point on the loading path.In case of non-conforming contacts, the Hertz contact parameters (Hertz, 1895) provide a good guess value for the estimated domain.If during application of additional loading increments the current contact area reaches the boundaries of the initial mesh, the contact simulation has to be restarted with a larger domain.However, only a small surface domain P  of the contacting bodies needs to be considered, which constitutes an important advantage of SAM over other numerical techniques.In P  , contact geometry should be known, or can be extrapolated from existing data.The directions of the grid sides are aligned with those of the Cartesian coordinate system in the continuous problem formulation.The elementary cell area . The grid control points (centroids of rectangular elementary patches) are identified by a pair of indices ( , ) i j , with x is assumed to be constant over each patch, and equal to value computed in the control point.A simplified notation can thus be used, assuming that , where * * 1 2 , x x are the coordinates of the control point of cell ( , ) i j .
The limiting surfaces of the contacting bodies are sampled in two height arrays corresponding to grid control points.Such data can be obtained from an optical profilometer or can be generated numerically.The sum of the two heights at node ( , ) i j yields the composite initial surface height ( , ) hi i j .For the half-space approximation to remain valid, the slope of the initial separation should be small everywhere over P  .
To simulate the loading history, additional discretization is performed in the temporal domain, meaning the load is imposed in small steps 1, k N   .Consequently, ( , , ) p i j k denotes the nodal (elementary) pressure at the intersection of the line i with the column j of the rectangular grid, achieved after application of k loading increments.When only two indexes are employed, the referred quantity does not vary with the loading level, e.g.coordinates of grid nodes; one index is used for parameters varying with the loading level only, e.g.rigid-body translations; when no index is present, the denoted quantity is a constant in the numerical program, e.g. the grid steps.Digitization of Eqs. ( 1) -( 4) lead to the following numerical model, which will be referred from now on as the NC: ( , , ) 0 ( , , ) 0, ( , ) ( ); ( , , ) 0 ( , , ) 0, ( , ) ( ).
In these equations, , A A and S A are the discrete counterparts of P  , C  and S  , respectively, consisting in sets of elementary patches.Consequently, in the numerical formulation, the contact or the stick area can only vary in equal increments  .A fine mesh is thus required for accurate estimation of contact domains.

Numerical computation of displacement fields
In Contact Mechanics, it is common practice to assimilate the contacting bodies with elastic half-spaces.This assumption holds if the dimensions of contact area are small compared to significant dimensions of contacting bodies, and allows expressing displacements according to superposition principle applied to Green functions for the elastic half-space: where induced by a unit point force applied in origin along direction of j x  .The functions G for the elastic half-space, also referred to as fundamental solutions or Green functions, were derived by Boussinesq, (Boussinesq, 1969) and Cerruti (Cerruti, 1882).Integral in Eq. ( 17) cannot be performed analytically except for a few cases; however, its discrete counterpart can be computed numerically for any contact area and any distribution of tractions: where ( ) ( , )  is the influence coefficient, expressing the displacement in the direction of x   , induced in the cell ( , ) i j of the body ( ) n by a unit contact traction, uniformly distributed in the cell ( , ) k  , acting along direction of x   .The influence coefficients yield from integration of Green functions over rectangular elementary patches: , ( ) .
Eq. ( 18) is in fact a two-dimensional convolution product, and can be written in an equivalent manner using the symbol "  " to denote discrete cyclic convolution: The assembly of contributions of all contact tractions to displacements can be expressed as: ( When the second subscript is omitted, the referred displacement is the sum of all contributions.The positive sign in computation of 3 u is related to the fact that normal displacements are computed in coordinate systems linked to each body, having the 3 x axes pointing inward.It should also be noted that from pressure definition (2)   (1) , 1,2 , as the mutual actions of two bodies upon each other are equal, but directed to contrary parts.The influence coefficients ( )   n m K  can be computed using the following relations: , , ln ln ; , , 2 tan ln ; 2 ) ln ; , , ; , , , , , , ; ,  , , ; ,  , , ; , , ; where ( ) are the elastic constants (Young modulus and Poisson's ratio, respectively) of materials of the two contacting bodies.If (1)   (2) ( , , 1,2,3, and Eq. ( 21) takes a simplified form: The most efficient way to compute the two-dimensional convolution products in Eqs. ( 21 However, when transferring to frequency domain via fast Fourier transform (FFT), presumption of signal periodicity is automatically assumed, which leads to contamination of convolution result due to spurious neighbouring periods.Liu & al. (Liu & al., 2000) advanced a fast and robust algorithm to circumvent the periodicity error, which is also applied here.This algorithm is able to assess linear convolution by computing the discrete cyclic convolution of the two terms on an extended domain, with a special, "warp-around" arrangement of influence coefficients.The base for this approach is discussed in detail in (Press et al., 1992).

Algorithm description
Considering the similarity between NC and TC, the same type of algorithm could be used to solve either model considered independently.Indeed, in both cases, a linear system arising from Eqs. ( 11) and ( 15) respectively, is to be solved, having as unknowns the normal or the tangential tractions.The main difficulties consist in the fact that the systems are essentially non-linear, due to presence of terms containing rigid-body motions, and their size, related to the number of cells in C A and S A respectively, is not known in advance.Moreover, resolution of a linear system with N unknowns has an order of computation of points is required to reproduce the detailed features of rough contact behaviour.Consequently, only fast numerical methods, based on iterative approaches, are suitable for this type of problem.According to a review paper by Allwood (Allwood, 2005), the Conjugate Gradient Method (CGM) has the fastest convergence.Implementation of CGM to this type of problem is authorised by the fact that the influence coefficients matrix is symmetrical and positive definite, therefore convergence is assured.
The algorithm used in this study is based on the modified CGM advanced by Polonsky and Keer (Polonsky & Keer, 1999) for the study of frictionless rough contact problems under normal loading.Existence and uniqueness of the solution is discussed by Kalker and van Randen (Kalker & van Randen, 1972).This algorithm was later enhanced by Spinu and Diaconescu (Spinu & Diaconescu, 2008) to allow for a bending moment on the contact area, and it was recently applied by Spinu and Frunza (Spinu & Frunza, 2011a) to solve numerically the Cattaneo-Mindlin problem.

Numerical solution of uncoupled contact problems
Eq. ( 21) suggests that, in case of dissimilarly elastic materials (i.e. (1)  (2)    ), computation of displacements in either direction require the knowledge of contact tractions on every direction.From this yields the coupling between the NC and the TC.However, in case of similarly elastic materials, Eq. (32) proves that normal displacements are independent of shear tractions, and tangential displacements are independent of pressure.In the latter case, solution of NC is decoupled from that of TC; however, pressure distribution is still needed in the TC for assessment of shear tractions in the slip zone, according to Coulomb's law.In the framework developed herein, we refer to the solution of uncoupled NC or TC as being the solution of NC when shear tractions are known, but otherwise arbitrarily distributed, as well as the solution of TC when pressure is known, but otherwise arbitrarily distributed.
It should be noted that in some cases, matching names will be used for variables having the same role in the algorithm structure, although their content might be different in the NC from that in the TC.Let us assume that a new loading increment k is applied.In both NC and TC, all parameters are specific to the achieved loading level; therefore, the third formal parameter, related to the loading history, will be omitted for brevity.An additional symbol "  " is used to denote the increment of the referred quantity after application of the k-th loading increment, i.e. 1 1 1 ( , ) ( , , ) ( , , 1) The contact tractions corresponding to the new loading level ( ) p k or ( ) k q can be assessed using the same type of algorithm, consisting in the sequences described in the following paragraphs.
Here, , 1,2,3 i d i  is the initial descend direction in the CGM.The role of the other variables will be discussed later.
2. Adopt the initial guess for contact tractions, using the static force equilibrium.In a first approximation, it will be considered in the NC that all cells are in contact, i.e. .The contact tractions are sought in the following form: NC p i j a a x i j a x i j i j A q i j a x i j a TC i j A q i j a x i j a where the unknown coefficients ( ), 1,2,3 i a k i  are computed by plugging Eq. (34) into Eqs.( 9), ( 10) and ( 13) , ( 14), respectively, yielding: :

A i jA i j A i j A i j A i j A i j A i j A
x i j x i j x i j x i j x i j x i j 3. Compute the relative displacement field over the grid domain using Eq. ( 21).It should be remembered that in the NC, pressure was adopted in step 2 and shear tractions are assumed as known, but arbitrarily distributed, while in the TC, shear tractions were adopted in step 2, and pressure is assumed known but arbitrarily distributed.It should also be noted that in case of TC, the increment of displacement field, 1 2 ( , ), ( , ) , is needed.4. Estimate the rigid-body motions in order to linearize Eqs. ( 11) and ( 15).According to complementarity conditions in Eqs. ( 12) and ( 16), the following relations hold: , ) 0, ( , ) ; ( , ) ( , ) 0 : , ( , ) , ( , ) ( , ) 0 resulting in linear systems having a number of equations equal the number of cells in the contact and in the stick area, respectively.When convergence is reached, the equations in every of these two systems are not all independent; however, during iterations, they appear to be independent.To get the best possible estimates of rigidbody motions, the systems are treated as over-determined and the best fit is sought using least square minimization.The functional  is defined as the square sum of the residuals: , ) .
The rigid-body motions that minimize the goal function in Eq. ( 38) are found by setting the partial derivatives of  to zero, leading to the following solution: With this approximation, Eqs. ( 11) and ( 15) finally reduce to linear systems having the nodal contact tractions as unknowns.The systems sizes are given by the extents of the contact and of the stick area, respectively, both a priori unknown.A trial-and-error approach is used, in which the systems are solved on successive contact and stick areas, until static equilibrium equations and complementarity conditions are fully satisfied for every cell in the grid.
5. Start the conjugate gradient loop by computing the residuals ( , ), 1,2,3 r i j    : , ), ( , ) ; ( , ) ( , ) ( , ) : , ( , ) , ( , ) ( , ) ( , ) P S NC r i j hi i j u i j x i j x i j i j A r i j u i j x i j TC i j A r i j u i j x i j and the square sum of the residuals on the indicated domains: ,) .
7. Assess the length of the step to be made in the computed descend direction: where , 1,2,3 i c i  denotes the convolution of the influence coefficients matrix with the descend direction: (2) (1) 22  22 : ; : .
8. Memorize the current contact tractions in a new variable for subsequent error computation: , 9. Update the system solution by making a step of length  , computed using Eq. ( 43), along direction , 1,2,3 , derived in Eq. ( 42): C S NC p i j p i j d i j i j A q i j q i j d i j TC i j A q i j q i j d i j 10. Impose complementarity conditions to adjust the size of the system.In the NC, the contact or non-contact status of every cell in the contact domain is to be determined, while in the TC, the stick or slip regime of every cell in the contact area must be assessed.
In the NC, cells having negative pressure are removed from the current contact area, and the corresponding nodal pressures are set to zero.It should be remembered that the current model cannot incorporate adhesion; therefore, normal contact tractions can only be compressive, not tensile.This assumption leads to specific cells passing from the contact zone to the non-contact area, in order to obey the non-adhesion hypothesis.
On the other hand, as bodies are considered impenetrable in the frame of Linear Theory of Elasticity, specific cells having vanishing pressure but negative gap ( , ) h i j (which equals the residual), are reinstated in the current contact area.To this end, the corresponding nodal pressures are adjusted with a quantity proven (Polonsky & Keer, 1999) to be positive at all times, i.e.
In the TC, cells for which Coulomb's friction law is not verified are removed from the stick region, and the corresponding shear tractions are set to the value of limiting friction.Also, cells having micro-slip ( , ) i j s not opposite to shear tractions ( , ) i j q pass from the slip to the stick region, yielding an adjusted stick area:

S S S S i j A A i j i j p i j i j i j i j p i j i j
In either the NC or the TC, if any cell is removed from or re-enters the contact or the stick area, auxiliary variable  is set to zero, otherwise it is set to unity.This variable allows resetting the conjugate directions i d once new cells enter or leave the system.Indeed, these 12. Verify convergence to the imposed precision  : , and return to step 2 if the precision goals are not met.The algorithm can be summarized in the flow-chart depicted in Fig. 3.

Numerical solution of coupled contact problems
In the general case, the NC and the TC cannot be solved independently, as the two problems are coupled, i.e. solution of each one requires resolution of the other.While the NC can be decoupled from the TC by neglecting the normal displacements induced by tangential tractions, as in case of Goodman approximation, solution of NC is always required to solve the TC, as shear tractions are linked to pressure in the slip area through Coulomb's law.The solution of the partial slip-stick elastic contact problem, involving the coupled NC and TC, can be obtained in an iterative manner, using the algorithm consisting in the following steps.
. The latter constitutes a better approximation to coupled problem solution compared to shear tractions adopted in step 2. 6. Solve the uncoupled NC again, using as input the shear tractions computed in step 5.
The resulting pressure distribution ( , , ) p i j k is in its turn a better approximation to coupled problem solution than the one obtained in step 3. 7. Verify if pressure distributions resulted in two subsequent computations, ( , , ) p i j k and ( , , ) old p i j k , vary within an imposed precision goal.If convergence is not met, return to step 4.This algorithm can be used to simulate any loading history in the fully coupled three dimensional elastic contact with slip and stick.Essentially, three levels of iterations are employed.The inner level is responsible for solving either the normal (NC), or the tangential (TC) unconnected contact problem.The intermediate level mutually adjusts these solutions, based on the contribution of tractions in each direction to displacement fields.The outer level is employed to reproduce the loading history, and is related to friction being a dissipative process.Spinu and Frunza (Spinu & Frunza, 2011b) recently proved that results obtained by simulating the loading history are a closer match to existing analytical solutions than when the full load is applied in one step.

Results and discussions
The newly advanced algorithm is validated in this section by comparison with existing closed-form solutions for the contact of similarly elastic materials.The results are then extended by simulating the three-dimensional contact between dissimilarly elastic materials, proving the capacity of the program to solve a large variety of problems incorporating the slip-stick processes in elastic mechanical contacts.

Simulation of a fretting loop
The newly proposed algorithm is first validated by comparison with the closed-form solution advanced in (Johnson, 1985) for the spherical contact undergoing a fretting loop.A steel ball of radius 18 R mm  is pressed against an elastic half-space having the same elastic properties, with a normal force 1 W kN  . In this case, the NC is uncoupled from the TC, as shown by Eq. ( 32).A tangential force 1 T oscillating between two limiting values lim T and lim T  , where lim 0.9 T W   , is subsequently applied.During a fretting contact process, it is expected that friction vary on the contact area, as well as with accumulation of debris particles resulted from additional wear.However, for validation purposes, a frictional coefficient uniform over all contact area and constant during load application is assumed in this study.The numerical approach can equally handle mapped distributions of  , if such information is available.The loading history for fretting processes is depicted in Fig. 4, in which the load is applied incrementally in 700 N   equal steps.N  is chosen differently in the following simulations.When studying the contact between similarly elastic materials, accurate (i.e.concurring with the existing closed-form solution) results can be obtained with even a small number of increments ( 42 N   increments in simulations depicted in Figs. 5 and 6).Moreover, it is found that every different state on the loading path can be obtained by simulating only the states when the tangential load changes its sign of, e.g.every state M on the trajectory DF require computation of three states: the ones corresponding to points B and D, as well as the current one (corresponding to point M).This is not the case when simulating the contact between dissimilarly elastic materials.It is found that a large number of increments is required to obtain the detailed contact behaviour, due to non-overlapping stick regions from one loading step to the next.This feature was also observed by Gallego et al. (Gallego et al., 2010), who pointed out that a waved shear tractions profile is predicted numerically when the number of increments is small.Although some perturbations still appear, the large number of increments used in the current simulations 700 N   led to well converged numerical solutions., which are used as normalizers.Shear tractions profiles corresponding to equal loading levels but laying on different trajectories are depicted in Fig. 5, proving that the current state depends not only on the present load level, but also on its past evolution, thus showing a hysteretic behaviour or memory effect.Only profiles corresponding to states on AD are shown in Fig. 5. States on the DF trajectory can be obtained from those on AD by symmetry with respect to the 1 xaxis, and states on FH overlap those on BD corresponding to the same loading level.As depicted in Fig. 6, the force-displacement curve forms a hysteretic loop, also referred to as a fretting loop.The rigid-body tangential displacement is normalized by its value lim  corresponding to lim T .The analytical model predicts that the states corresponding to points B and F on the loading curve should overlap, which is also obtained through numerical simulation, within an imposed precision.
From a mechanical point of view, regions undergoing the most severe stresses are of interest, leading to prediction of yield inception or crack nucleation in the contacting bodies.An algorithm to assess stresses due to known, but otherwise arbitrarily distributed contact tractions is readily available (Liu & Wang, 2002).A typical distribution of von Mises equivalent stress VM H p  in the plane 2 0 x  , corresponding to point B on the loading path, is depicted in Fig. 7. Evolution of magnitude and of depth of von Mises maximum stress in a fretting loop is depicted in Figs.8-9, for various frictional coefficients. .However, larger friction processes lead to a competition between the in-depth maximum, which fluctuates around its position in the frictionless contact (corresponding to point A in Fig. 9), and a second extremum developing on the surface, at the trailing edge of the contact.The latter can seriously diminish the loadcarrying capacity of the contact, especially in case of surfaces with poorly controlled surface quality, due to superimposition of microtopography-induced stress perturbations.

Simulation of torsional contact
Program validation is subsequently performed in case of a spherical contact undergoing torsion applied simultaneously to a normal constant force.Based on the results of this author, (Mindlin, 1949), Johnson (Johnson, 1985) presents the closed-form solution for this contact scenario when a partial slip regime is established on the contact area (i.e. when the torsional moment 3 M is smaller than a limiting value 3 lim 3 1 6 ).The solution is later reviewed and enhanced for the case of viscoelastic materials by Dintwa et al. (Dintwa et al., 2005).
A spherical indenter of radius 18 R mm  is pressed with a normal force 1 W kN  against an elastic half-space, having the same elastic parameters, and subsequently an increasing torsional moment 3 3l i m M M  is applied, leading to shear tractions 1 q and 2 q depicted in Figs.10a and b, respectively.The norm of these tractions 1 1 ( , ) ( , ) ( , ) i j q i j q i j   q is compared in Fig. 11 with analytical results, using cylindrical coordinates.
Figure 12 depicts the dimensionless stick radius when the torsional moment 3 M is varied in the domain corresponding to partial slip.In all investigated cases, a good agreement between analytical results and numerical predictions is found, giving further confidence in the newly advanced computer program.

Extension of results
The numerical program advanced in this study is subsequently used to simulate the fretting contact between dissimilarly elastic materials.To our best knowledge, an analytical solution to this contact process has not been advanced yet.The ball in the previous example is considered rigid, and the loading history is simulated using 700 N   equal increments.Numerical simulations predict that the investigated fretting contact exhibits a unique behaviour in the first two loading cycles, after which it stabilizes to a fixed path.Figures 13 and 14 suggest that the shear tractions profiles corresponding to the BD and FH trajectories no longer match, as in case of similarly elastic materials.However, states D and H overlap, beside a few perturbations induced by the former boundaries of the stick area.Presumably, these perturbations are related to the discretization error, and the D and H states can be considered as concurring.Subsequent oscillating loading is expected to lead to states following the same fixed path, as in case of similarly elastic materials.An analogous behaviour was observed by Wang et al. (Wang et al., 2010) when studying numerically the partial slip contact of elastic layered half-spaces.

Conclusions
The work reported herein advances a numerical model for the fretting contact between dissimilarly elastic materials.A numerical approach is required to simulate this type of contact process, as analytical models can incorporate neither the loading history, which must be reproduced when friction is accounted for, nor the coupling of normal and tangential effects.
The implemented algorithm is based on three levels of iterations, fully incorporating the interconnectivity between normal and tangential tractions.The innermost level solves, in a conjugate gradient type loop, either the normal, or the tangential uncoupled contact problem, while the intermediate loop iterates between these solutions, until pressure convergence is reached.The outer level reproduces the loading history, and is based on the assumption that irreversibility of friction requires simulation of all previous states in an incremental load process.
The strong points of this algorithm consist in reduced computational complexity, compared to finite element analysis, as well as in its ability to handle arbitrarily shaped contact geometries or imposed frictional coefficient maps.Comparison with existing closed-form solutions for the spherical contact undergoing a fretting loop or torsion, when a uniform frictional coefficient is assumed, gives confidence in the newly advanced model.Evolution of stress state during a fretting loop in the spherical contact between similarly elastic materials is assessed.It is found that for specific combinations of loading levels and frictional coefficients, the most severely stressed region can be found on the bounding surfaces, at the trailing edge of the contact.
Numerical simulations suggest that the fretting contact between dissimilarly elastic materials exhibits a unique path in the first two loading cycles, which stabilizes with subsequent oscillating loading to a fixed trajectory, as in case of similarly elastic materials.
The newly advanced algorithm is expected to solve a large variety of contact problems involving interfacial friction, leading to a better understanding of complex multidisciplinary phenomena like fretting wear and fretting fatigue, in which transient contact tractions and induced subsurface stresses play an important role, as well as to a more accurate prediction of contact failure through yield inception or crack nucleation.Study of partial slip elasticplastic contact is anticipated for future contributions, by addition of the residual term, related to development of plastic strains.

Figure 1 .
Figure 1.Geometrical condition of deformation in the normal direction

Figure 2 .
Figure 2. Geometrical condition of deformation in the tangential direction ) or (32) is by using spectral methods.According to convolution theorem, the convolution of two signals, each having N samples, requires 2 ( ) O N operations in time/space domain, but only ( log ) O N N operations in the frequency domain, where it reduces to element-wise product.
like Gaussian elimination are used.In contact problems, grids with 3 10 N  nodes on the contact area are usually considered for nominal contact surfaces, in order to get an accurate estimation of contact area extents and to minimize the discretization error (i.e. the error introduced by digitization of problem parameters).When deterministic rough surfaces are studied, a minimum of 6 10 N  TC that all cells in the contact area are in stick, i.e. S C A A 

Figure 3 .
Figure 3. Flow-chart for the uncoupled NC or TC 1. Acquire the state after 1 k  loading increments and establish the new increment.2. Adopt vanishing guess value for the increment of shear tractions ( , ) i j q , i.e. ( , , ) ( , , 1), ( , )( 1)C i j k i j k i j A k     q q .3.Solve the uncoupled NC using the algorithm described in the previous section, and obtain pressure ( , , ),( , ) ( )C p i j k i j A k  .4.Memorize the obtained pressure for subsequent error computation: ( , , ) ( , , ) old p i j k p i j k  . 5. Solve the uncoupled TC having as input the previously computed pressure, and obtain the shear tractions ( , , ),( , ) ( )

Figure 4 .
Figure 4.The loading history in simulation of a fretting loop; all moments are assumed to vanish

Figure 5 .
Figure 5. Profiles of shear tractions in the plane 2 0x  corresponding to different points on the loading curve,

Figure 13 .Figure 14 .
Figure 13.Shear tractions profiles in the fretting contact between dissimilarly elastic materials, corresponding to different points on the loading path, 0.3   The Hertz frictionless theory for the similarly elastic contact scenario predicts a central pressure 1.996