A Three-Dimensional Semianalytical Model for Elastic-Plastic Sliding Contacts

A three-dimensional numerical model based on a semianalytical method in the framework of small strains and small displacements is presented for solving an elastic-plastic contact with surface traction. A Coulomb’s law is assumed for the friction, as commonly used for sliding contacts. The effects of the contact pressure distribution and residual strain on the geometry of the contacting surfaces are derived from Betti’s reciprocal theorem with initial strain. The main advantage of this approach over the classical ﬁnite element method (FEM) is the computing time, which is reduced by several orders of magnitude. The contact problem, which is one of the most time-consuming procedures in the elastic-plastic algorithm, is obtained using a method based on the variational prin-ciple and accelerated by means of the discrete convolution fast Fourier transform (FFT) and conjugate gradient methods. The FFT technique is also involved in the calculation of internal strains and stresses. A return-mapping algorithm with an elastic predictor/plastic corrector scheme and a von Mises criterion is used in the plasticity loop. The model is ﬁrst validated by comparison with results obtained by the FEM. The effect of the friction coefﬁcient on the contact pressure distribution, subsurface stress ﬁeld, and residual strains is also presented and discussed.


Introduction
The rolling contact fatigue phenomenon is involved in many mechanical components, such as rolling element bearings, gears, cams, and continuously variable transmissions. The fatigue life of the contacting surfaces is directly related to the stress state found at the surface and within the material in the vicinity of the contact. There exists a certain stress threshold below which the life of the contact is infinite, sometimes called the endurance limit. Lamagnère et al. ͓1͔ proposed an endurance limit concept based on the comparison between the maximum shear stress and the microyield stress. This implies an accurate knowledge of the stress field history and of the yield stress of the material. In some high demanding applications, rolling element bearings or gears are designed so that the applied load stays below the endurance limit. A purely elastic approach is then sufficient. However, it is almost impossible to guarantee that the yield stress will never be reached locally or accidentally all over the life of the component. The knowledge of the plastic strains and hardening state of the material is then of prime importance to evaluate the remaining life of the contacting surfaces. An attempt is made here to account for the contribution of the tangential loading on the stress and strain states found in an elastic-plastic sliding contact, in addition to purely normal effects.
An analytical relationship between the surface profile and contact pressure exists only for a limited number of ideal geometries ͑Hertz's theory͒. The Hertzian pressure distribution is strongly modified by the presence of a geometrical defect, such as roughness, ridge, furrow, or a dent such as those produced by a debris when entrapped within the contact conjunction. High local pressure peaks appear around the dent, producing high concentration of stress localized in the vicinity of the dent. Usually, the yield stress is exceeded and plastic flow occurs, and should be added to the initial plasticity introduced during indentation ͓2,3͔.
Recently, Jacq et al. ͓4͔ proposed a semianalytical elasticplastic method to solve 3D contact problems, fast enough to study a vertical loading/unloading or a moving load, for example, to simulate the rolling of a body on a surface defect. In this model, the contact pressure distribution is found to be modified from the purely elastic case mostly by plasticity induced change of the contacting surface conformity, which tends to limit the contact pressure to 2.8 or 3 times the yield stress, approximately, depending on the hardening law used and the contact geometry. Both subsurface strains and stresses are also found to be strongly modified, but this time due to a combination of three effects: ͑i͒ the occurrence of plastic strains, ͑ii͒ the material hardening, and ͑iii͒ the change of the contact pressure distribution. Based on the original algorithm the work of Jacq et al. has been improved in two ways: first, by also considering thermal effect in the elastic-plastic algorithm ͓5͔ to account for a surface heat source: second, by introducing the return-mapping algorithm with an elastic predictor and a plastic corrector scheme in the plasticity loop ͓6͔. Based on the same algorithm, Wang and Keer ͓7͔ studied the effect of various hardening laws on the elastic-plastic response of the indented material.
The presence of small size surface defects in a larger contact area requires a fine mesh of the contact area ͑up to 10 6 surface grid points͒. To reduce the computing time significantly, the contact module, originally based on a multilevel technique in Ref. ͓4͔, was replaced by a module based on a single-loop conjugate gradient method ͑GM͒ developed by Polonsky and Keer ͓8͔. A discussion on the efficiency of this method initially used to solve elastic rough contact problems is given by Allwood ͓9͔. This method was improved by the implementation of discrete convolution fast Fourier transform ͑DC-FFT͒ approaches as presented by Liu et al. ͓10͔ for the calculation of surface displacements and internal stress field.
The effect of a tangential loading of the surface, not considered in the original code, has been introduced in the elastic-plastic formulation and implemented in the computer code. Surface traction is included in the model as a shear stress proportional to the normal pressure by the use of a friction coefficient, suitable for sliding contacts. This paper presents first the theory and the numerical algorithm proposed to solve the elastic-plastic contact problem with surface traction. The distributions of contact pressure, residual stress and plastic strain are obtained and compared with the results of the finite element method ͑FEM͒. Finally, the effect of the traction coefficient is presented and discussed.

Theory and Numerical Procedure
The theory is based on the formulation of Betti's reciprocal theorem. The numerical scheme is based on the use of DC-FFT and CGM.
Hypotheses. The two bodies in contact are considered as halfspaces so the contact area is small with regard to the dimensions of the bodies in contact. The normal and tangential effects are treated separately and are uncoupled in the contact solver. It means that the tangential displacement of the surface points is not known; the boundary condition related to traction is expressed in terms of shear stress. Small strains and small displacements are assumed for the plasticity model, which allows us to limit the plastic analysis to the volume where yielding occurs while superposing residual strains to the elastic part. Higher residual strains as those resulting from a Rockwell indentation could be calculated using a finite element software and then introduced as an initial state, as far as the overrolling of the surface does not produce large additional strains ͓2,4͔.

Contact Model
Theory. According to Johnson ͓11͔, one can solve the problem of two nonconforming bodies in contacts B1 and B2 using two different methods: the direct or matrix inversion method and the variational method. Kalker ͓12͔ proposed a variational principle in which the true contact area and contact pressure are those that minimize the total complementary energy ͑V * ͒, subjected to the constraints that the contact pressure is everywhere positive, and there is no interpenetration, where ⌫ c is the surface on which p acts and U E * is the internal complementary energy of two stressed bodies, numerically equal for elastic materials with the elastic strain energy U E , which can be expressed in terms of the surface pressure and the normal displacements of both bodies by Finally, the problem is reduced to a set of equations ͑two equalities and two inequalities͒ that should be solved simultaneously.
h͑x 1 ,x 2 ͒ ജ 0 and p͑x 1 ,x 2 ͒ ജ 0 ͑5͒ If h͑x 1 ,x 2 ͒ Ͼ 0 then p͑x 1 ,x 2 ͒ = 0 ͑6͒ The first system's equation is the static equilibrium condition ͑Eq. ͑3͒͒. The two bodies are in equilibrium under the external force W and the resultant contact force due to the pressure p͑x 1 , x 2 ͒ over the contact surface ⌫ c . The final separation h͑x 1 , x 2 ͒ is expressed in Eq. ͑4͒ as a summation of h i ͑x 1 , x 2 ͒, the initial distance between the bodies, ␦, the rigid body displacement, and u 3 ͑B1+B2͒ ͑x 1 , x 2 ͒, the surface normal displacement of the two bodies B1 and B2. The elastic surface displacement u 3 ͑B1+B2͒ ͑x 1 , x 2 ͒ can be expressed by the use of the Boussinesq relation where x and are two points on the surface, and U͑x , ͒ is the displacement at point x due to a unit load at point ͑influence functions or Green's functions͒. Further, the reciprocal theorem will be used to include the surface residual displacement due to subsurface plasticity for the elastic-plastic contact problem. The set of Eqs. ͑5͒ and ͑6͒, also known as contact criteria, expresses the fact that the bodies cannot interpenetrate one another and pressure is nil outside the contact area.
Numerical Resolution. As ⌫ c is unknown, the procedure to solve the system of Eqs. ͑3͒-͑6͒ is iterative. Several iterative methods, such as Jacobi, Gauss-Seidel, and CGM, are available to solve this system. Among them, the CGM has important advantages.
͑1͒ A rigorous mathematical proof of the method convergence exists; ͑2͒ It offers a very high rate of convergence ͑superlinear͒; ͑3͒ It requires a very modest storage capacity that is extremely advantageous when very large systems of equations are involved.
The CGM is an iterative method that generates a sequence of approximations of the solution starting from an arbitrary initial approximation. The recurrence formula of the CG is where i =0,1,2..; r i ͑residue͒ and d i ͑direction͒ are vectors of N elements; P 0 is an arbitrary start vector ͑for example, a uniform pressure͒; and d 0 = This scheme, originally presented by Polonsky and Keer ͓8͔ and then based on the CGM and a multilevel multisummation method ͑MLMSM͒, was used to build the computer code. An important particularity of the iteration process is that the contact area is established in the course of the pressure iteration, so that there is no need for further iteration with respect to the contact area. Another distinctive feature of the iteration scheme used is that the force balance equation is enforced during each iteration for the contact pressure. The rigid deflection ␦ is not explicitly used in the iteration process but may easily be determined, as long as the pressure distribution is already known.
The most time-consuming works in the CGM are the multiplication operations between the influence coefficient matrix U by the pressure vector p and direction vector d. These multiplication operations require O͑N 2 ͒ operations, and if N is large, these need a large amount of computing time. To reduce the computing time, the FFT is used within each iteration of the CGM for the task of multiplying both pressure and direction vectors by the influence coefficient matrix. The number of operations is reduced to O͑N log N͒, where N is the number of grid points involved.

Plasticity Model.
Plasticity is an irreversible phenomenon that requires an incremental description. In a general incremental for-mulation of plasticity, a plastic strain increment depends on the preexisting plastic strain, the stress, the stress increment, and the hardening parameters, This general formulation is used for theoretical development. The subsurface residual stresses are calculated following the method proposed by Chiu ͓13,14͔, considering a cuboidal zone with uniform initial strains or eigenstrains and surrounded by an infinite elastic space ͓13͔ or a half-space ͓14͔. The calculation of the elastic stress field due to the contact pressure is more classical. The influence coefficients that give the stresses induced by a rectangular cell on the boundary surface of a half-space submitted to a uniform pressure are given in the Appendix.
Betti's Reciprocal Theorem. Consider two displacement fields, u and u * , continuous and two times derivable into an elastic body of volume ⍀ and of boundary ⌫. u corresponds to an elastic state with initial inelastic strains 0 and u * currently unknown. The elastic properties are the same for both states. Thus, one can write the following for stress and strain: It is to be noticed, that 0 is the summation of all the inelastic strains, i.e., thermal strains, plastic strains, etc. The reciprocal theorem ͓15͔ is generalized for the elastic case with initial inelastic strains in Refs. ͓16,17͔. For a detailed development of the formulation of the reciprocal theorem in the case of elastoplasticity, the reader can refer to Ref. ͓4͔, and in the case of thermoelastoplasticity, the reader can refer to Ref. ͓5͔. These formulations use the theoretical developments of Chiu ͓13,14͔ for calculating residual stresses and the theoretical developments of Liu and Wang ͓18͔ for calculating thermally induced stresses.
Betti's reciprocal theorem could be applied for an elastic state with initial strain ͑Brebbia ͓19͔; Mayeur et al. ͓20͔͒. Finally, according to Jacq et al. ͓4͔, the reciprocal theorem can be written with indicial notation as follows: Application to the Calculation of Surface Displacements. Equation ͑11͒ is now applied to both bodies in contact, where each of them is considered as a half-space ⍀ whose boundary ⌫ is loaded on a part ⌫ c and with a uniform coefficient of friction. Since the body forces are neglected ͑f i =0͒ and ij n j =−p i ͑values of p i corresponding to either the pressure for i = 3 or the shear for i =1 or 2͒, Eq. ͑11͒ becomes where ⍀ p is the initial volume occupied by the initial strains. Let us consider the state ͑u * , * , * , f i * ͒ as the application of a unit normal force at point A of the contact area, along x 3 , neglecting from now on the body forces ͑f i * =0͒. There is thus a uniform pressure p * = ͑A͒ =1/d⌫ at point A within an elementary surface d⌫ = dx 1 dx 2 , and p * ͑M A͒ = 0 elsewhere. Then, the first term on the right-hand side of Eq. ͑12͒ equals the normal displacement at point A for the nonstar state. Due to the coefficient of friction, a surface shear stress proportional to the normal pressure exists, t͑A͒ = p͑A͒ and t͑M A͒ = 0. Johnson ͓11͔ showed that if the two bodies have the same elastic constants, there is no influence of the surface shear stress over the sum of the normal displacements. Therefore, where M is point of the integration surface or volume and where the plastic strain p corresponds to the initial state 0 . The notation ͓M, p * ͑A͔͒ means the value at point M due to the application of the pressure p * applied at point A. For more details on how Eq. ͑13͒ is found, the reader may refer to Refs. ͓4,5͔.
The surface normal displacement of each body can then be expressed as a function of the contact pressure and of the plastic strain existing in the considered body. Hence, Application to Stress Calculation. Consider the reciprocal theorem applied again to a half-space ⍀ whose boundary ⌫ is loaded on a part ⌫ c . Because the plastic strains are related to stresses, the stress field must be evaluated from elastic-plastic contact conditions. Consider now the state ͑u ** , ** , ** , f i ** ͒ corresponding to a body force applied at point B of the half-space in the direction k and of magnitude 1, Hence, Stresses can be related to the displacements via Hooke's law. Therefore, the stress at every point of the half-space can be divided in three parts, In the first two terms, the pressure stress is linked to the normal contact pressure and surface shear stress ͑or tangential loading͒, while the residual stress is related to the plastic strains in the last term. The residual stress is the stress induced by the strain nuclei. It is also the stress produced by plastic strain remaining after unloading.
This relation shows that the stress field changes with plastic flow, primarily due to the appearance of plastic strains, but also because of the modification of contact pressure due to geometry changes.

Coupling Elastic and Plastic Parts.
Since plasticity is an irreversible phenomenon, the relation between plastic strain and contact pressure must also be incremental. Therefore, an incremental formulation of the elastic-plastic contact problem must be used.

͑19͒
Load balance: Surface separation: the elastic displacement due to contact pressure and tangential loading, the displacement due to plastic strain, where the displacement due to the plastic strain increment is ␦ p .
Plasticity model: ␦ p = f͑,␦, hardening parameters͒ Stress calculation: Contact conditions: The algorithm developed to solve the incremental elastic-plastic contact with friction problem is similar to the one proposed by Jacq et al. ͓4͔. The initial state may include residual strains. The elastic contact is first solved using the CGM, with any initial surface separation. The plasticity model is then used to calculate the plastic strain increment, considering also the effect of surface shear stress, enabling the calculation of the residual displacement increment. The residual surface displacement increment, which is a function of the plastic strain, is then calculated and compared to the one found during the previous step. This process is repeated until the residual displacement increments converge. Plastic strains, normal and tangential loads, contact pressure, residual surface displacement, and yield strength are then increased by their increment to define a new initial condition for the next loading step.
Plasticity Loop Improvement. The plasticity loop for the equivalent plastic deformation calculation was improved for convergence and accuracy needs. The proposed method is based on the work of Fotiu and Nemat-Nasser ͓21͔, who built a universal algorithm for the integration of the elastic-plastic constitutive equations. Isotropic and kinematic hardening, as well as thermal softening, may be used in the formulation. This method is unconditionally stable and accurate. The return mapping algorithm with an elastic predictor/plastic corrector scheme that was implemented in the code is briefly presented in Ref. ͓6͔. Compared to the previous algorithm, formerly based on the Prandtl-Reuss model ͓4͔, no more plasticity loop is needed, the computation of the plastic strains being reduced to a few number of iterations ͑typically one to four iterations for the return-mapping procedure using the Newton-Raphson scheme͒, resulting in a drastic reduction of the CPU time ͑at least by one order of magnitude͒ while improving the quality of the solution.

Some Results
The aim of this section is, first, to validate the model based on the semianalytical method ͑SAM͒ through comparison with FEM results in case of normal and tangential loadings and, second, to discuss the effect of surface traction on stress and strain fields. For convenience, all results will be presented for a circular point contact between a sphere of radius 10 mm and a half-space. The sphere will be first assumed rigid when comparing SAM and FEM numerical results. The sphere will behave elastically when investigating the effect of surface traction on stress and strain states. The elastic properties are E = 210 GPa and = 0.3 for Young's modulus and Poisson's ratio, respectively. The isotropic hardening behavior and the von Mises stress criterion are used in the plasticity model. The half-space properties are those of the M50 steel. The corresponding hardening law is described by Swift's law ͑Eq. ͑27͒͒, whose parameters are B = 1280 MPa, C = 4, and n = 0.095 ͓2͔. In Eq. ͑27͒, the equivalent plastic strain e p is expressed in the 10 −6 deformation.
Some results will be presented as normalized by the microyield stress e , defined as the stress corresponding to a proof strain of 20ϫ 10 −6 . The microyield stress threshold is useful in rolling contact fatigue applications, where it results in a satisfactory estimation of the endurance limit of high strength steels ͓1͔. Equation  ͑27͒ gives a microyield stress of 1731 MPa for the M50 steel, in contrast to the more conventional yield stress of 2636 MPa corresponding to a proof strain of 2 ϫ 10 −3 for the same steel.

Validation by Comparison With Finite Element Method
Results. To validate the results supplied by the code, a comparison with the commercial finite element code ABAQUS is made. For simplicity, the upper surface is modeled here as a rigid sphere. The 3D FE model includes 113,035 C3D8 bricks ͑linear elements with eight nodes͒ with a refined mesh in the contact zone ͑more than 300 elements in contact͒. The meshed volume is extended to 20 times the Hertz contact radius a in each direction ͑−20a to 20a along x and y, 0 to 20a along z͒, assumed to be sufficient for contact analysis by FEM ͓22͔. Boundary conditions are zero displacements for x = ±20a, y = ±20a, and z =20a. The normal load is imposed by a vertical displacement of the upper surface of the rigid sphere, up to a load of 800 N corresponding to a Hertz pressure of 4.35 GPa and a contact radius of 296 m for an elastic analysis. Figures 1 and 2 present a comparison between SAM and FEM results when the loaded half-space behaves elastically. As a whole, a good agreement is found for the pressure distribution and the stress profile ͑Figs. 1 and 2, respectively͒. However, some differences can be noticed. First of all, it is obvious that the pressure distribution ͑Fig. 1͒ corresponds to the Hertz solution with the SAM code, i.e., independent of the friction coefficient, due to the fact that normal and tangential effects are assumed uncoupled in the present model. Conversely, the FEM results exhibit an asymmetric pressure distribution when increasing surface traction, with a shift of the maximum contact pressure in the direction of the traction force. The von Mises stress profile along the geometrical axis of symmetry is shown in Fig. 2. The difference for a frictionless contact ͑Fig. 2͑a͒͒ between SAM and FEM numerical results is attributed to the boundary conditions in the FE model, which do not correspond exactly to an infinite half-space. Differences observed between the two methods of analysis when increasing the friction coefficient to 0.2 and 0.4 ͑Figs. 2͑b͒ and 2͑c͒͒ are explained by the two above mentioned reasons: an asymmetry of the pressure distribution and boundary conditions given at a finite distance from the contact for the FEM analysis. Finally, the results also confirm a classical result, which is that the maximum stress is found at the surface when the friction coefficient exceeds 0.3 ͓11͔.
The effects of plasticity are presented in Figs. 3-5, where the half-space behaves now as an elastic-plastic media, the load being still transmitted through a rigid sphere. Again, SAM and FEM results are globally in good agreement, with some slight differences that have the same origins as for the elastic simulations. The pressure distribution ͑Fig. 3͒ is affected by the hardening of the material, the maximum contact pressure being lowered. At the same time the contact area tends to increase, keeping the integral of the contact pressure constant. This time the pressure distribution is no more found to be symmetric with the SAM approach while increasing the friction coefficient ͑Figs. 3͑b͒ and 3͑c͒͒. However, the asymmetry is more pronounced with the FEM analysis, where tangential and normal effects are implicitly coupled. The von Mises stress found under loading and the equivalent plastic strain are given in Figs. 4 and 5, respectively, along the depth at the center of the contact. The main difference, which is observed in Fig. 5͑c͒ at the surface of the contact, is again attributed to a more pronounced shift of the maximum contact pressure with FEM compared to SAM, in addition to nonideal boundary conditions. This is a consequence of uncoupling the tangential and normal effects when solving the contact problem. Despite some differences, the fairly good agreement between FEM and SAM numerical results found in Figs. 1-5 validates the elastic-plastic contact solver being currently developed. The identification of the origins of theses differences will be used in the near future to improve the SAM code. The advantage of the SAM approach in terms of computing time should be recalled, with less than 1 min CPU time for all results presented here, compared with approximately 1 d of CPU time on the same personal computer for an equivalent FEM analysis.

Influence of Surface Traction on Stress and Strain (Semianalytical Method Results).
Numerical simulations are now performed to investigate the effect of normal and tangential loadings on the contact pressure distribution as well as on subsurface stress and strain fields found under load or after unloading ͑residual͒. From this point, the sphere will be modeled as an elastic body, still loaded against an elastic-plastic half-space. The friction coefficient will range from 0 to 0.5, and the normal load from 1000 N to 5000 N, the latter corresponding to a Hertz ͑elastic͒ contact pressure of 5.05 GPa and a normalized interference / c0 of 3.2. The critical interference c0 and the critical load W c0 are the normal deflection and the corresponding normal load at the onset of yielding, as introduced by Chang et al. ͓23͔. Most results are presented for the maximum normal load, i.e., 5000 N.
It should be noted that the critical interference and load, c0 and W c0 , respectively, are usually defined for a pure normal loading ͓23͔. Considering now the effect of a tangential loading superimposed to the normal load, one may easily plot the ratios c / c0 and W c / W c0 , where c and W c are the values found for a contact with friction. In Fig. 6, a very important effect of the friction coefficient on these critical values is found.
The maximum of the equivalent plastic strain versus the Hertz pressure, normalized by the microyield stress, is presented in Fig.  7 for different friction coefficients . The corresponding contact pressure distributions at the normal load of 5000 N are presented in Fig. 8. It should be noted that for = 0.4 and 0.5, the maximum of the plastic strain is found at the surface of the half-space, whereas it is found at the Hertzian depth for the frictionless contact. Another interesting result is the magnitude of the plastic strain, which reaches 5.34% for = 0.5 and a normal load of 5000 N, remaining acceptable with regard to the assumption of small strains. Finally, the marked effect of the friction coefficient should be noted, which increases drastically the maximum of the plastic strain from 0.34% for the frictionless contact at 5000 N up to 5.34% at = 0.5.
The von Mises stress profile at the center of the contact is shown in Fig. 9 for friction coefficient ranging from 0 to 0.5, first, under load with a normal load of 5000 N ͑Fig. 9͑a͒͒ and, second, after unloading ͑Fig. 9͑b͒͒. Three contributions of the plasticity are associated with the decrease of the maximum von Mises stress observed at the Hertzian depth in Fig. 9͑a͒: first, a change of the surface conformity due to the subsurface residual strain, second, the modification of the pressure distribution, and, third, the hardening of the material. Outside the plastic zone the stress profile follows the elastic solution ͑visible in Fig. 9͑a͒ for the frictionless contact͒. The location of the point where the maximum stress is found is presented in Fig. 10. Once again, the maximum stress is found at the surface of the half-space when the friction coefficient exceeds 0.32 ͑see Fig. 10͒. In addition, Fig. 10 indicates that this point moves away from the contact center along the traction direction ͑see curve x / a versus ͒ up to the critical value of 0.32, for which the maximum reaches the surface, but in the opposite direction. The discontinuity in the plots of x / a and z / a versus the friction coefficient is explained by the fact that two local maxima are competing, one in the Hertzian region located slightly ahead of the contact center ͑i.e., x Ͼ 0, see left part of the curve͒ and another closer to the surface and at the trailing edge of the contact. The maximum found near or at the surface is mostly related to the tangential load, whereas the maximum found in the Hertzian region is mostly related to the normal load. A further increase of the friction coefficient will finally move back that point to the contact center. More interesting is the residual stress profile found after unloading ͑Fig. 9͑b͒͒. One may observe two zones where residual stress is present, one at the Hertzian depth and another at the surface of the contact, including the frictionless contact. The same trend was previously observed by Jackson et al. ͓24͔ and Kadin et al. ͓25͔ for frictionless hemispherical contacts. The magnitude of the residual stress found at the surface increases with the friction coefficient to become higher than the one found at the Hertzian depth when the friction coefficient exceeds 0.3. In contrast, note a pronounced decrease of the maximum residual stress found at the Hertzian depth when the friction coefficient becomes higher than 0.4. It can be seen in Fig. 9͑b͒ that jagged lines are found for the lowest coefficients of friction. This has two reasons. First, the residual stress level is very low. Second, the mesh is probably not fine enough to capture the very high gradient found at this location, as shown with a map view in Fig. 11.
A profile of the equivalent plastic strain is given in Fig. 12 at the vertical of the contact center. One may observe that the plastic zone reaches the surface of the elastic-plastic body for a friction coefficient of approximately 0.2. Other simulations have shown that this critical friction coefficient is dependent on the normal load ͑or plasticity level͒, decreasing from 0.3 to 0 when the normal load increases from the value corresponding to first yielding.  This result is of prime importance for wear or running-in modeling when the material removal is based on a strain threshold criterion ͓6͔.
When residual stresses are present, it is sometimes important to know if it corresponds to tensile or compressive zones. Typically, for rolling contact fatigue applications, it is well known that compressive residual stress will close the crack faces, whereas tensile residual stress will favor the crack initiation and, later, its propagation. An interesting stress quantity for that identification is the hydrostatic pressure P hydr , defined here as where a positive value means a compressive zone and a negative one, a tensile zone. Figure 13 presents the residual hydrostatic pressure profile along the depth at the contact center for the frictionless contact and for various normal loads ranging from 1000 to 5000 N. An interesting result is the succession of compressive and tensile zones: compressive at the surface and at the Hertzian depth, tensile between the surface and the Hertzian depth, and tensile below the Hertzian depth. This is coherent with the observation of Kogut and Etsion ͓26͔ in terms of plastic strains for an axisymmetric contact. The amplitude of the variation of the hydrostatic pressure increases with the normal load. The tensile zone found between the surface and the Hertzian depth, sometimes called the "quiescent zone" for a rough contact ͓27͔, may explain why cracks initiated at the surface or at the Hertzian depth may propagate toward the Hertzian depth or toward the surface, respectively ͓28͔. In a similar manner, Fig. 14 presents the effect of surface traction on the same hydrostatic pressure profile after unloading for the highest normal load of 5000 N. Similar comments as those given for Fig. 9 could be made about the two local maxima found at the surface and at the Hertzian depth. Note also the effect of the  friction coefficient on the stress state in the quiescent zone, which becomes in compression above a critical friction coefficient found here between 0.2 and 0.3.
The examples presented here illustrate the performances of the SAM proposed by the authors. This method is an alternative to the use of the FEM, which remained until now almost the only tool for studying the elastic-plastic effects of a geometrical surface defect on the fatigue of the contacting materials ͑see, for example, Ref. ͓29͔͒.

Conclusion
A three-dimensional semianalytical elastic-plastic contact code, taking into account both normal and tangential loadings, has been developed. The contact solver, which is based on the CGM and FFT technique, allows to solve the transient elastic-plastic contact problem within a reasonable computing time even when a fine mesh is required, i.e., within a few minutes up to a few hours ͑for 10 6 grid points͒ on a PC, mostly depending on the number of cells inside the plastic zone. The proposed method is an alternative to the use of the FEM, for example, in the study of the effects of a geometrical surface defect on the fatigue of the contacting materials.
The model has been first validated by comparison with 3D FEM results obtained with the commercial software ABAQUS. The comparison has pointed out the side effect of solving the contact problem without coupling normal and tangential effects.
The stress and strain states after the vertical loading/unloading of an elastic-plastic half-space by an elastic sphere have been investigated for a steady-state problem with combined normal and tangential loadings and compared to the purely normal loading case. The results presented have shown a significant effect of the tangential loading not only on the magnitude and location of the maximum von Mises stress found under loading, but also on the residual stresses and strains that remain after unloading. An interesting point is the existence of a residual tensile zone between the Hertzian depth and the surface, which was identified by means of the hydrostatic stress.

Convention
a i,j ϭ ‫ץ‬a i / ‫ץ‬x j

Appendix: Stress Calculation in a Half-Space Loaded on Surface
The main purpose of this section is to give expressions for the subsurface stress field due to a uniformly distributed load over a rectangle area. The rectangle, with sides 2a ϫ 2b and centered at the origin, is subjected to uniform pressure p in the normal direction and uniform tractions t x and t y in the tangential directions.