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

Daniel Nélias, Eduard Antaluca, Vincent Boucly, Spiridon Cretu

► To cite this version:

Daniel Nélias, Eduard Antaluca, Vincent Boucly, Spiridon Cretu. A Three-Dimensional Semianalytical Model for Elastic-Plastic Sliding Contacts. Journal of Tribology, 2007, 129 (4), pp.761-771. 10.1115/1.2768076 . hal-01962873

HAL Id: hal-01962873

https://hal.utc.fr/hal-01962873

Submitted on 28 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 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 finite 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 principle 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 first validated by comparison with results obtained by the FEM. The effect of the friction coefficient on the contact pressure distribution, subsurface stress field, and residual strains is also presented and discussed.

**Keywords:** elastic-plastic contact, SAM method, contact mechanics, stress analysis, dry contact

**Daniel Nélias**

LaMCoS, INSA Lyon,
CNRS UMR 5259,
Villeurbanne F69621, France
e-mail: daniel.nelias@insa-lyon.fr

**Eduard Antaluca**

LaMCoS, INSA Lyon,
CNRS UMR 5259,
Villeurbanne F69621, France;
Technical University “Gh. Asachi,”
Iasi 700050, Romania

**Vincent Bouly**

LaMCoS, INSA Lyon,
CNRS UMR 5259,
Villeurbanne F69621, France

**Spiridon Cretu**

Technical University “Gh. Asachi,”
Iasi 700050, Romania

## 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 elastic-plastic 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 half-spaces 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,

$$V^* = U_E^* + \int_{\Gamma_c} p(h - \delta)ds \quad (1)$$

where $\Gamma_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

$$U_E^* = U_E = \frac{1}{2} \int_{\Gamma_c} p(u_3^{B1} + u_3^{B2})d\Gamma \quad (2)$$

Finally, the problem is reduced to a set of equations (two equalities and two inequalities) that should be solved simultaneously.

$$W = \int_{\Gamma_c} p(x_1, x_2)d\Gamma \quad (3)$$

$$h(x_1, x_2) = h_i(x_1, x_2) + \delta + u_3^{(B1+B2)}(x_1, x_2) \quad (4)$$

$$h(x_1, x_2) \ge 0 \quad \text{and} \quad p(x_1, x_2) \ge 0 \quad (5)$$

$$\text{If } h(x_1, x_2) > 0 \text{ then } p(x_1, x_2) = 0 \quad (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 $\Gamma_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, $\delta$, 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

$$u(x) = \int_{\Gamma_c} U(x, \xi)p(\xi)d\Gamma \quad (7)$$

where $x$ and $\xi$ are two points on the surface, and $U(x, \xi)$ is the displacement at point $x$ due to a unit load at point $\xi$ (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 $\Gamma_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

$$\mathbf{p}_{i+1} = \mathbf{p}_i - \frac{\mathbf{r}_i^T \mathbf{r}_i}{\mathbf{d}_i^T U \mathbf{d}_i} \mathbf{d}_i \quad (8a)$$

$$\mathbf{r}_{i+1} = \mathbf{r}_i - \frac{\mathbf{r}_i^T \mathbf{r}_i}{\mathbf{d}_i^T U \mathbf{d}_i} U \mathbf{d}_i \quad (8b)$$

$$\mathbf{d}_{i+1} = -\mathbf{r}_{i+1} + \frac{\mathbf{r}_{i+1}^T \mathbf{r}_{i+1}}{\mathbf{r}_i^T \mathbf{r}_i} \mathbf{d}_i \quad (8c)$$

where $i=0, 1, 2, ...$; $\mathbf{r}_i$ (residue) and $\mathbf{d}_i$ (direction) are vectors of $N$ elements; $\mathbf{P}_0$ is an arbitrary start vector (for example, a uniform pressure); and $\mathbf{d}_0 = \mathbf{r}_0 = B - U\mathbf{P}_0$, with $B = h - h_i - \delta$ (see Eq. (4)).

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 $\delta$ 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 $\mathbf{p}$ and direction vector $\mathbf{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,

$$ \delta\varepsilon^p = f(\varepsilon^p, \sigma, \delta\sigma, \text{hardening parameters}) \quad (9) $$

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 $\Omega$ and of boundary $\Gamma$. $u$ corresponds to an elastic state with initial inelastic strains $\varepsilon^0$ and $u^*$ currently unknown. The elastic properties are the same for both states. Thus, one can write the following for stress and strain:

$$ \begin{aligned} \varepsilon_{ij} &= \frac{1}{2}(u_{i,j} + u_{j,i}) & \varepsilon_{ij}^* &= \frac{1}{2}(u_{i,j}^* + u_{j,i}^*) \\ \sigma_{ij} &= M_{ijkl}(\varepsilon_{kl} - \varepsilon_{kl}^0) & \sigma_{ij}^* &= M_{ijkl}\varepsilon_{kl}^* \end{aligned} \quad (10) $$

It is to be noticed, that $\varepsilon^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:

$$ \begin{aligned} -\int_{\Gamma} u_i^* \sigma_{ij} n_j d\Gamma + \int_{\Omega} f_i u_i^* d\Omega = & -\int_{\Gamma} u_i \sigma_{ij}^* n_j d\Gamma + \int_{\Omega} f_i^* u_i d\Omega \\ & -\int_{\Omega} \varepsilon_{ij}^0 \sigma_{ij}^* d\Omega \end{aligned} \quad (11) $$

**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 $\Omega$ whose boundary $\Gamma$ is loaded on a part $\Gamma_c$ and with a uniform coefficient of friction. Since the body forces are neglected ($f_i=0$) and $\sigma_{ij}n_j=-\mathbf{p}_i$ (values of $\mathbf{p}_i$ corresponding to either the pressure for $i=3$ or the shear for $i=1$ or 2), Eq. (11) becomes

$$ \int_{\Gamma_c} u_i^* \mathbf{p}_i d\Gamma = \int_{\Gamma} u_i \mathbf{p}_i^* d\Gamma + \int_{\Omega} f_i^* u_i d\Omega - \int_{\Omega_p} \varepsilon_{ij}^0 \sigma_{ij}^* d\Omega \quad (12) $$

where $\Omega_p$ is the initial volume occupied by the initial strains.

Let us consider the state $(u^*, \varepsilon^*, \sigma^*, 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 $\mathbf{p}^*=(A)=1/d\Gamma$ at point A within an elementary surface $d\Gamma=dx_1dx_2$, and $\mathbf{p}^*(M \neq 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)=\mu\mathbf{p}(A)$ and $t(M \neq 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,

$$ \begin{aligned} u_3(A) = & \int_{\Gamma} u_i(M) \mathbf{p}_i^*(M) d\Gamma = \int_{\Gamma_c} u_3^*[M, \mathbf{p}^*(A)] \mathbf{p}_3(M) d\Gamma \\ & + 2\mu_L \int_{\Omega_p} \varepsilon_{ij}^p(M) \varepsilon_{ij}^*[M, \mathbf{p}^*(A)] d\Omega \end{aligned} \quad (13) $$

where $M$ is point of the integration surface or volume and where the plastic strain $\varepsilon^p$ corresponds to the initial state $\varepsilon^0$. The notation $[M, \mathbf{p}^*(A)]$ means the value at point $M$ due to the application of the pressure $\mathbf{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,

$$ u_3(A) = u^e(A) + u^r(A) \quad (14) $$

**Application to Stress Calculation.** Consider the reciprocal theorem applied again to a half-space $\Omega$ whose boundary $\Gamma$ is loaded on a part $\Gamma_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^{**}, \varepsilon^{**}, \sigma^{**}, f_i^{**})$ corresponding to a body force applied at point B of the half-space in the direction $k$ and of magnitude 1,

$$ \begin{aligned} & \int_{\Omega} f_i^{**} u_i d\Omega = u_k(B) & \int_{\Gamma} u_{ni} \mathbf{p}_i^{**} d\Gamma = 0 \\ & \int_{\Gamma} u_{ii} t_i^{**} d\Gamma = 0 \quad \text{and } t = \mu \mathbf{p} \end{aligned} \quad (15) $$

Equation (11) then becomes

$$ \begin{aligned} u_k(B) = & \int_{\Gamma_c} u_{ni}^{**}(M, B) \mathbf{p}_i(M) d\Gamma + \int_{\Gamma_c} u_{ii}^{**}(M, B) t_i(M) d\Gamma \\ & + 2\mu_L \int_{\Omega_p} \varepsilon_{ij}^p(M) \varepsilon_{kij}^{**}(M, B) d\Omega. \end{aligned} \quad (16) $$

Hence,

$$ u_k(B) = u_k^n(B) + u_k^l(B) + u_k^r(B) \quad (17) $$

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,

$$ \sigma = \sigma^n + \sigma^l + \sigma^r \quad (18) $$

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.

Initial conditions:
$$W, h_i(x_1, x_2), p(x_1, x_2), \varepsilon^p, \text{hardening parameters} \quad (19)$$

Load balance:

$$W + \delta W = \int_{\Gamma_c} [p(x_1, x_2) + \delta p(x_1, x_2)] d\Gamma \quad (20)$$

Surface separation:

$$h(x_1, x_2) = h_i(x_1, x_2) + \delta + u_3^{e(B1+B2)}(x_1, x_2) + u_3^r(x_1, x_2) + \delta u_3^r(x_1, x_2) \quad (21)$$

with

$$u_3^e(A) = \int_{\Gamma_c} u_{3i}^*(M, A)(p_i + \delta p_i)(M) d\Gamma \quad (22)$$

the elastic displacement due to contact pressure and tangential loading,

$$u_3^r(A) = 2\mu_L \int_{\Omega_p} \varepsilon_{ij}^p(M) \varepsilon_{3ij}^*(M, A) d\Omega \quad (23)$$

the displacement due to plastic strain,

$$\delta u_3^r(A) = 2\mu_L \int_{\Omega_p} \delta \varepsilon_{ij}^p(M) \varepsilon_{3ij}^*(M, A) d\Omega \quad (24)$$

where the displacement due to the plastic strain increment is $\delta\varepsilon^p$. Plasticity model:

$$\delta\varepsilon^p = f(\sigma, \delta\sigma, \text{hardening parameters})$$

Stress calculation:

$$\sigma = \sigma^e(p) + \sigma^r(\varepsilon^p)$$

$$\delta\sigma = \delta\sigma^e(\delta p) + \delta\sigma^r(\delta\varepsilon^p) \quad (25)$$

Contact conditions:

$$h(x_1, x_2) \ge 0 \quad (26a)$$

$$\text{if } h(x_1, x_2) > 0 \text{ then } p(x_1, x_2) + \delta p(x_1, x_2) = 0 \quad (26b)$$

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 $\nu=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.

$$\sigma_{VM} = B(C + e^p)^n \quad (27)$$

Some results will be presented as normalized by the microyield stress $\sigma_c$, defined as the stress corresponding to a proof strain of $20 \times 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 \times 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=\pm 20a$, $y=\pm 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
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 1 Comparison of numerical results (SAM versus FEM, elastic solution, contact pressure distribution) for various friction coefficients: (a) $\mu=0$, (b) $\mu=0.2$, and (c) $\mu=0.4$</figcaption></figure>

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

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 2 Comparison of numerical results (SAM versus FEM, elastic solution, von Mises stress under load (profile at $x=y=0$)) for various friction coefficients: (a) $\mu=0$, (b) $\mu=0.2$, and (c) $\mu=0.4$</figcaption></figure>

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
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 3 Comparison of numerical results (SAM versus FEM, elastic-plastic solution, contact pressure distribution) for various friction coefficients: (a) $\mu=0$, (b) $\mu=0.2$, and (c) $\mu=0.4$</figcaption></figure>

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 these differences will be used in the near future to improve the SAM code. The advantage of the SAM

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 4 Comparison of numerical results (SAM versus FEM, elastic-plastic solution, von Mises stress under load (profile at $x=y=0$)) for various friction coefficients: (a) $\mu=0$, (b) $\mu=0.2$, and (c) $\mu=0.4$</figcaption></figure>

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 (Semi-analytical Method Results).** Numerical simulations are now performed to investigate the effect of normal and tangential loadings
<figure><img src="image_14.png" /><figcaption>Fig. 5 Comparison of numerical results (SAM versus FEM, elastic-plastic solution, equivalent plastic strain (profile at x=y=0)) for various friction coefficients: (a) $\mu=0$, (b) $\mu=0.2$, and (c) $\mu=0.4$</figcaption></figure>

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 $\omega/\omega_{c0}$ of 3.2. The critical interference $\omega_{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, $\omega_{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 $\omega_c/\omega_{c0}$ and $W_c/W_{c0}$, where $\omega_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.

<figure><img src="image_15.png" /><figcaption>Fig. 6 Influence of the friction coefficient on the critical load (circle symbols) and interference (square symbols) at the onset of yielding, normalized by the values for the frictionless contact</figcaption></figure>

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 $\mu$. The corresponding contact pressure distributions at the normal load of 5000 N are presented in Fig. 8. It should be noted that for $\mu=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 $\mu=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 $\mu=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

<figure><img src="image_16.png" /><figcaption>Fig. 7 Maximum of the equivalent plastic strain versus the corresponding Hertzian contact pressure normalized by the microyield stress for various friction coefficients</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 8 Pressure distribution for various friction coefficients. Normal load of 5000 N.</figcaption></figure>

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 $\mu$) 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

<figure><img src="image_18.png" /><figcaption>Fig. 10 Value and location of the maximum von Mises stress (under load) versus the friction coefficient</figcaption></figure>

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.

<figure><img src="image_19.png" /><figcaption>Fig. 9 von Mises stress profile at $x=y=0$. Normal load 5000 N. (a) Under load; (b) after unloading (residual).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 11 Residual von Mises stress profile in the plane $y=0$. Normal load of 5000 N; frictionless contact.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Fig. 12 Equivalent plastic strain profile versus depth at x=y=0. (a) Regular view. (b) Zoom of (a).</figcaption></figure>

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_{\text{hydr}}$, defined here as

$$P_{\text{hydr}} = - \frac{\sigma_{ii}}{3} \qquad (28)$$

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

<figure><img src="image_22.png" /><figcaption>Fig. 13 Hydrostatic pressure after unloading (residual), profile at x=y=0, frictionless contact, for various normal loads ranging from 1000 N to 5000 N</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 14 Hydrostatic pressure after unloading (residual), profile at x=y=0, normal load of 5000 N, for various friction coefficients ranging from 0 to 0.5</figcaption></figure>
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.

## Acknowledgment

The authors would like to acknowledge support for this research by the Agence Universitaire de la Francophonie (AUF) through the Ph.D. grant of E. Antaluca.

## Nomenclature

B, C, and n = Swift's law parameters

B1, B2 = body 1 and body 2

d = direction vector (CGM)

E = Young modulus, Pa

e$^p$ = equivalent plastic strain

f = body force, N m$^{-3}$

h = surface separation, m

h$_i$ = initial surface separation, m

M$_{ijkl}$ = elastic constant matrix

p = pressure, Pa

P$_0$ = initial pressure, Pa

P$_{hydr}$ = hydrostatic pressure, Pa

r = residue (CGM)

t = surface shear stress, Pa

U = influence coefficients for surface normal displacement

u = surface displacement, m

U$^*$ = internal complementary energy, J

U$_E$ = internal complementary energy, J

U$_E$ = elastic strain energy, J

u$^e$ = u$^n$+u$^t$=elastic normal displacement, m

u$^n$ = normal surface displacement due to normal loading, m

u$^t$ = residual displacement, m

u$^t$ = normal surface displacement due to tangential loading, m

V$^*$ = total complementary energy, J

W = normal load, N

W$_c$ = critical normal load at the onset of yielding, N

W$_{c0}$ = critical normal load at the onset of yielding for frictionless contact, N

Γ$_c$ = contact surface

δ = interference, m

ε$^0$ = initial inelastic strain tensor

ε$^p$ = plastic strain tensor

μ$_L$ = Lamé coefficient

μ = coefficient of friction

σ = total stress tensor, Pa

σ$^e$ = σ$^n$+σ$^t$=elastic stress tensor, Pa

σ$^n$ = elastic stress tensor due to normal loading, Pa

σ$^t$ = residual stress tensor, Pa

σ$^t$ = elastic stress tensor due to tangential loading, Pa

ω$_c$ = critical interference at the onset of yielding, m

ω$_{c0}$ = critical interference at the onset of yielding for frictionless contact, m

## Convention

$$a_{i,j} = \partial a_i / \partial 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 \times 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.

Normal load uniformly distributed over a rectangle. The stress components due to a uniform pressure $p$ over a rectangle $2a \times 2b$ are given by

$$\sigma_{ij} = \frac{p}{2\pi} [F_{ij}(x+a, y+b, z) - F_{ij}(x+a, y-b, z) + F_{ij}(x-a, y-b, z) - F_{ij}(x-a, y+b, z)]$$

with the following functions:

$$F_{xx}(x,y,z) = 2\nu \tan^{-1}\left(\frac{z^2 + y^2 - y\rho}{zx}\right) + 2(1 - 2\nu)\tan^{-1}\left(\frac{\rho - y + z}{x}\right) + \frac{xyz}{\rho(x^2 + z^2)}$$

$$F_{yy}(x,y,z) = 2\nu \tan^{-1}\left(\frac{z^2 + y^2 - y\rho}{xz}\right) + 2(1 - 2\nu)\tan^{-1}\left(\frac{\rho - x + z}{y}\right) + \frac{xyz}{\rho(y^2 + z^2)}$$

$$F_{zz}(x,y,z) = \tan^{-1}\left(\frac{y^2 + z^2 - y\rho}{xz}\right) - \frac{xyz}{\rho}\left(\frac{1}{x^2 + z^2} + \frac{1}{y^2 + z^2}\right),$$

$$F_{xy}(x,y,z) = \frac{-z}{\rho} - (1 - 2\nu)\ln(\rho + z)$$

$$F_{xz}(x,y,z) = \frac{z^2 y}{\rho(x^2 + z^2)}$$

$$F_{yz}(x,y,z) = \frac{z^2 x}{\rho(y^2 + z^2)}$$
Tangential load uniformly distributed along x over a rectangle. The stress components due to a uniform traction $t_x$ along x over a rectangle $2a \times 2b$ are given by

$$ \sigma_{ij} = \frac{t_x}{2\pi} [G_{ij}(x+a, y+b, z) - G_{ij}(x+a, y-b, z) + G_{ij}(x-a, y-b, z) - G_{ij}(x-a, y+b, z)] $$

with the following functions:

$$ G_{xx}(x,y,z) = \frac{-z}{\rho} \left\{ 1 + \frac{yz - x^2}{[(\rho + z)(\rho - y)]} \right\} + 2\nu \left[ \frac{y}{(\rho + z)} \right] - 2\ln(\rho - y) $$

$$ G_{yy}(x,y,z) = \frac{-yz}{[\rho(\rho + z)]} - 2\nu \left[ \frac{y}{\rho + z} + \ln(\rho - y) \right] $$

$$ G_{zz}(x,y,z) = \frac{yz^2}{\rho(x^2 + z^2)} $$

$$ G_{xy}(x,y,z) = \frac{-xz}{[\rho(\rho + z)]} - 2\nu \left[ \frac{x}{(\rho + z)} \right] - \ln(\rho - x) $$

$$ G_{xz}(x,y,z) = \frac{xyz}{\rho(x^2 + z^2)} + \tan^{-1} \left( \frac{z^2 + y^2 - y\rho}{xz} \right) $$

$$ G_{yz}(x,y,z) = \frac{-z}{\rho} $$

Tangential load uniformly distributed along y over a rectangle. The stress components due to a uniform traction $t_y$ along y applied over a rectangle $2a \times 2b$ are obviously related to the previous contribution of $t_x$ and corresponding functions G.

## References

[1] Lamagnère, P., Fougères, R., Lormand, G., Vincent, A., Girodin, D., Dudagne, G., and Vergne, F., 1998, "A Physically Based Model for Endurance Limit of Bearing Steels," ASME J. Tribol., 120, pp. 421-426.

[2] Nélias, D., Jacq, C., Lormand, G., Dudagne, G., and Vincent, A., 2005, "A New Methodology to Evaluate the Rolling Contact Fatigue Performance of Bearing Steels With Surface Dents—Application to 32CrMoV13 (Nitrided) and M50 Steels," ASME J. Tribol., 127, pp. 611-622.

[3] Vincent, A., Nélias, D., Jacq, C., Robin, Y., and Dudagne, G., 2006, "Comparison of Fatigue Performances of 32CrMoV13 and M50 Steels in Presence of Surface Indents," J. ASTM Int., 3(2), pp. 1-9.

[4] Jacq, C., Nélias, D., Lormand, G., and Girodin, D., 2002, "Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code," ASME J. Tribol., 124, pp. 653-667.

[5] Boucly, V., Nélias, D., Liu, S., Wang, Q. J., and Keer, L. M., 2005, "Contact Analyses for Bodies With Frictional Heating and Plastic Behavior," ASME J. Tribol., 127, pp. 355-364.

[6] Nélias, D., Boucly, V., and Brunet, M., 2006, "Elastic-Plastic Contact Between

Rough Surfaces: Proposal for a Wear or Running-in Model," ASME J. Tribol., 128, pp. 236-244.

[7] Wang, F., and Keer, L. M., 2005, "Numerical Simulation for Three Dimensional Elastic-Plastic Contact With Hardening Behavior," ASME J. Tribol., 127, pp. 494-502.

[8] Polonsky, I. A., and Keer, L. M., 1999, "A Numerical Method for Solving Rough Contact Problems Based on the Multi-Level Multi-Summation and Conjugate Gradient Techniques," Wear, 231, pp. 206-219.

[9] Allwood, J., 2005, "Survey and Performance Assessment of Solution Methods for Elastic Rough Contact Models," ASME J. Tribol., 127, pp. 10-23.

[10] Liu, S., Wang, Q., and Liu, G., 2000, "A Versatile Method of Discrete Convolution and FFT (DC-FFT) for Contact Analyses," Wear, 243, pp. 101-111.

[11] Johnson, K. L., 1985, Contact Mechanics, Cambridge University Press, London.

[12] Kalker, J. J., 1990, Three Dimensional Elastic Bodies in Rolling Contact, Kluwer Academic, Dordrecht.

[13] Chiu, Y. P., 1977, "On the Stress Field Due to Initial Strains in a Cuboid Surrounded by an Infinite Elastic Space," ASME J. Appl. Mech., 44, pp. 587-590.

[14] Chiu, Y. P., 1978, "On the Stress Field and Surface Deformation in a Half-space With a Cuboidal Zone in Which Initial Strains Are Uniform," ASME J. Appl. Mech., 45, pp. 302-306.

[15] Brebbia, C. A., Telles, J. C. F., and Wrobel, L. C., 1984, Boundary Element Techniques. Theory and Applications in Engineering, Springer-Verlag, Berlin.

[16] Telles, J. C. F., and Brebbia, C. A., 1979, "On the Application of the Boundary Element Method to Plasticity," Appl. Math. Model., 3, pp. 466-470.

[17] Telles, J. C. F., and Brebbia, C. A., 1979, "The Boundary Element Method in Plasticity," New Developments in Boundary Element Methods, C. A. Brebbia, ed., Computational Mechanics Centre, Southampton, pp. 295-317.

[18] Liu, S., and Wang, Q., 2003, "Transient Thermoelastic Stress Fields in a Half-Space," ASME J. Tribol., 125, pp. 33-43.

[19] Brebbia, C. A., 1980, The Boundary Element Method for Engineers, Pentech, London.

[20] Mayeur, C., Sainsot, P., and Flamand, L., 1995, "A Numerical Elastoplastic Model for Rough Contact," ASME J. Tribol., 117, pp. 422-429.

[21] Fotiu, P. A., and Nemat-Nasser, S., 1996, "A Universal Integration Algorithm for Rate-Dependant Elastoplasticity," Comput. Struct., 59, pp. 1173-1184.

[22] Stephens, L. S., Yan, L., and Meletis, E. I., 2000, "Finite Element Analysis of the Initial Yielding Behaviour of a Hard Coating/Substrate System With Functionally Graded Interface Under Indentation and Friction," ASME J. Tribol., 122, pp. 381-387.

[23] Chang, W. R., Etsion, I., and Bogy, D. B., 1987, "An Elastic-Plastic Model for the Contact of Rough Surfaces," ASME J. Tribol., 109, pp. 257-263.

[24] Jackson, R., Chupoisin, I., and Green, I., 2005, "A Finite Element Study of the Residual Stress and Deformation in Hemispherical Contacts," ASME J. Tribol., 127, pp. 484-493.

[25] Kadin, Y., Kligerman, Y., and Etsion, I., 2006, "Multiple Loading-Unloading of an Elastic-Plastic Spherical Contact," Int. J. Solids Struct., 43, pp. 7119-7127.

[26] Kogut, L., and Etsion, I., 2002, "Elastic-Plastic Contact Analysis of a Sphere and a Rigid Flat," ASME J. Appl. Mech., 69, pp. 657-662.

[27] Tallian, T. E., Chiu, Y. P., and Amerongen, E. V., 1978, "Prediction of Traction and Microgeometry Effects on Rolling Contact Fatigue Life," ASME J. Lubr. Technol., 100, pp. 156-166.

[28] Nélias, D., Dumont, M.-L., Champiot, F., Vincent, A., Girodin, D., Fougères, R., and Flamand, L., 1999, "Role of Inclusions, Surface Roughness and Operating Conditions on Rolling Contact Fatigue," ASME J. Tribol., 121, pp. 240-251.

[29] Howell, M. B., Rubin, C. A., and Hahn, G. T., 2004, "The Effect of Dent Size on the Pressure Distribution and Failure Location in Dry Point Frictionless Rolling Contacts," ASME J. Tribol., 126, pp. 413-421.