**2. Formulation of continuous slip-stick elastic contact problem**

In contact problem formulation, it is convenient to describe the initial geometries ( ) 1 2 (,) *<sup>i</sup> hi x x* of the two contacting bodies *i* 1,2 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 1 2 **T**(,) *T T* , the bending (or flexing) moments 1 2 *M* , *M* and the torsional moment *M*<sup>3</sup> . 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 ( )*<sup>j</sup> <sup>i</sup> u* , and move as rigid-bodies with translations ( )*<sup>j</sup> i* and rotations ( )*<sup>j</sup> i* , with *i j* 1,2,3, 1,2 .

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 ( )*<sup>j</sup> p* in the normal direction and shear traction ( ) 1 2 (,) *<sup>j</sup>* **<sup>q</sup>** *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.

#### **2.1. The contact model in the normal direction**

Based on the works developed in (Johnson, 1985; Polonsky & Keer, 1999), the model for the contact in the normal direction consists in the following equations and inequalities:

1. The static force equilibrium:

130 Numerical Simulation – From Theory to Industry

required.

solutions.

load was applied in one step.

(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

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

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

**2. Formulation of continuous slip-stick elastic contact problem** 

In contact problem formulation, it is convenient to describe the initial geometries ( )

of the two contacting bodies *i* 1,2 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 1 2 **T**(,) *T T* , the bending

1 2 (,) *<sup>i</sup> hi x x*

$$\mathcal{W}(t) = \iint\limits\_{\Gamma\_{\mathbb{C}}(t)} p(\mathbf{x}\_1, \mathbf{x}\_2, t) d\mathbf{x}\_1 d\mathbf{x}\_2;\tag{1}$$

$$M\_1(t) = \iint\limits\_{\Gamma\_{\mathbb{C}}(t)} p(\mathbf{x}\_1, \mathbf{x}\_2, t) \mathbf{x}\_2 d\mathbf{x}\_1 d\mathbf{x}\_2; \ \ M\_2(t) = \iint\limits\_{\Gamma\_{\mathbb{C}}(t)} p(\mathbf{x}\_1, \mathbf{x}\_2, t) \mathbf{x}\_1 d\mathbf{x}\_1 d\mathbf{x}\_2. \tag{2}$$

2. The geometrical condition of deformation:

$$h(\mathbf{x}\_1, \mathbf{x}\_2, t) = h i(\mathbf{x}\_1, \mathbf{x}\_2) + u\_3(\mathbf{x}\_1, \mathbf{x}\_2, t) - a\_3(t) - \phi\_1(t)\mathbf{x}\_2 - \phi\_2(t)\mathbf{x}\_1, \quad (\mathbf{x}\_1, \mathbf{x}\_2) \in \Gamma\_C(t). \tag{3}$$

3. The contact complementarity conditions:

$$\begin{cases} p(\mathbf{x}\_1, \mathbf{x}\_2, t) > 0 \land h(\mathbf{x}\_1, \mathbf{x}\_2, t) = 0, & (\mathbf{x}\_1, \mathbf{x}\_2) \in \Gamma\_{\mathcal{C}} \{t\}; \\ p(\mathbf{x}\_1, \mathbf{x}\_2, t) = 0 \land h(\mathbf{x}\_1, \mathbf{x}\_2, t) > 0, & (\mathbf{x}\_1, \mathbf{x}\_2) \notin \Gamma\_{\mathcal{C}} \{t\}. \end{cases} \tag{4}$$

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, ( ) *<sup>C</sup> t* denotes the

contact area, and 1 2 *hx x t* ( , ,) the surface separation, both established at a specified point *t* on the loading curve. When 0 *t* , (1) (2) 12 12 12 12 *h x x hi x x hi x x hi x x* ( , ,0) ( , ) ( , ) ( , ) . The other relative (composite) quantities, lacking the superscript indicating the contacting body, are defined in a similar way: (1) (2) 33 3 () () () *ttt* , and (1) (2) ( ) ( ) ( ), 1,2 *ii i t t ti* . Computation of the relative normal displacement 312 *uxxt* ( , ,) will be discussed in Section 3.2.

Numerical Simulation of Slip-Stick Elastic Contact 133

(7)

(8)

, and points

 and (2) (1) 11 1 *uu u* ;

*S*

the frictional coefficient and 1 2 **s**(,) *s s*

. Rigid-body motion and elastic deformation due

<sup>1</sup> along the direction of 1 *x*

11 1 

  *C S*

**2.2. The contact model in the tangential direction** 

1. The static force equilibrium (at any time *t* ):

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,

12 12

(5)

(6)

32 31 12 2

*t t xx t*

( ) ( , , ) , 1,2;

3 212 1 112 2 1 2

2. The geometrical condition of deformation in the time frame 1 2 [,] *t t* , in which the

( , ,) ( , ,) ( , ,) ( , ,) ... ( , ,) ( , ,) ( , ,) ( , ,)

() () ( ) ( ) , ( , ) ( ). () () *<sup>C</sup>*

1 22 1 22 1 22 1 21 2 2

( , , ) ( , , ) ( , , ) ( , , ) 0, ( , ) ( ) ( ).

*x x t px x t x x t x x t ij t t*

the relative slip distances. The base for Eq. (7) is presented in Fig. 2, which depicts the slip-

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 *A*1 and *A*<sup>2</sup> , respectively, where the bodies deform due shear

Firstly, the above mentioned axis is assumed to pass through the initial point of contact, therefore 11 21 *A* () () *t At O* . In the considered time frame, all points in the contacting

1 and (2) 

 () () *t ut* , with (1) (2) 

*A*1 and *A*2 also undergo tangential displacements. If for these points the composite

therefore, 12 22 *A* () () *t At O* and consequently a stick regime is established in *O* .

 

1 22 1 22 1 22 1 21 2

( , , ) ( , , ) ( , , ) ( , , ) 0 , ( , ) ( );

*x x t px x t x x t x x t ij t*

12 11 2

*t t x*

22 21 1

*t t x*

1 1 22 1 1 21 1 1 22 1 1 21 2 1 22 2 1 21 2 1 22 2 1 21

*sxxt sxxt uxxt uxxt sxxt sxxt uxxt uxxt*

*M t q x x t x q x x t x dx dx*

() ( , ,) ( , ,) .

*T t q x x t dx dx i*

1985; Chen & Wang, 2008), consists in the following equations and inequalities:

( )

*C i i t*

( )

*<sup>C</sup> t*

tangential load is not allowed to change sign:

 

 

**q s s q s s**

3. The contact complementarity conditions:

tractions.

stick contact process in the direction of 1 *x*

bodies undergo rigid-body translations (1)

parameters cancel each other, i.e. 12 12

Here, *<sup>S</sup>* is the stick area, *C S* the slip region,

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) <sup>3</sup> *u t*( ) and (2) <sup>3</sup> *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.

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

#### **2.2. 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* ):

132 Numerical Simulation – From Theory to Industry

3.2.

3 

(1)

<sup>3</sup> *u t*( ) and (2)

body translations (1)

3 

( )*t* and (2)

defined in a similar way: (1) (2)

contact area, and 1 2 *hx x t* ( , ,) the surface separation, both established at a specified point *t* on

relative (composite) quantities, lacking the superscript indicating the contacting body, are

Computation of the relative normal displacement 312 *uxxt* ( , ,) will be discussed in Section

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.

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 rigid-

contacting bodies are compressed. In order to accommodate the interpenetration distance

 ( )*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

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.

<sup>3</sup> *u t*( ) pointing inward the corresponding body, as both normal contact tractions

 

33 3

12 12 12 12 *h x x hi x x hi x x hi x x* ( , ,0) ( , ) ( , ) ( , ) . The other

 () () () *ttt* , and (1) (2) ( ) ( ) ( ), 1,2 *ii i* 

( )*t* , which have opposite signs due to the fact that the

 *t t ti* .

the loading curve. When 0 *t* , (1) (2)

Therefore, study of adhesion effects is beyond the point of this study.

3 

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

$$T\_i(t) = \iint\_{\Gamma\_C(t)} q\_i(\mathbf{x}\_1, \mathbf{x}\_2, t) d\mathbf{x}\_1 d\mathbf{x}\_2, \quad i = 1, 2; \tag{5}$$

$$M\_3(t) = \iint\limits\_{\Gamma\_{\mathbb{C}}(t)} \left[ q\_2(\mathbf{x}\_1, \mathbf{x}\_2, t)\mathbf{x}\_1 - q\_1(\mathbf{x}\_1, \mathbf{x}\_2, t)\mathbf{x}\_2 \right] d\mathbf{x}\_1 d\mathbf{x}\_2. \tag{6}$$

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:

$$
\begin{bmatrix}
\mathbf{s}\_1(\mathbf{x}\_1, \mathbf{x}\_2, t\_2) - \mathbf{s}\_1(\mathbf{x}\_1, \mathbf{x}\_2, t\_1) \\
\mathbf{s}\_2(\mathbf{x}\_1, \mathbf{x}\_2, t\_2) - \mathbf{s}\_2(\mathbf{x}\_1, \mathbf{x}\_2, t\_1)
\end{bmatrix} = \begin{bmatrix}
u\_1(\mathbf{x}\_1, \mathbf{x}\_2, t\_2) - u\_1(\mathbf{x}\_1, \mathbf{x}\_2, t\_1) \\
u\_2(\mathbf{x}\_1, \mathbf{x}\_2, t\_2) - u\_2(\mathbf{x}\_1, \mathbf{x}\_2, t\_1)
\end{bmatrix} - \dots
$$

$$
\begin{bmatrix}
o\_1(t\_2) - o\_1(t\_1) \\o\_2(t\_2) - o\_2(t\_1)
\end{bmatrix} - \left(\phi\_3(t\_2) - \phi\_3(t\_1)\right) \begin{bmatrix}
\mathbf{x}\_2 \\
\mathbf{x}\_1
\end{bmatrix}, \quad \text{( $\mathbf{x}\_1, \mathbf{x}\_2$ ) \in \Gamma\_C(t\_2).} \tag{7}
$$

3. The contact complementarity conditions:

$$\begin{cases} \left\| \mathbf{q}(\mathbf{x}\_{1},\mathbf{x}\_{2},t\_{2}) \right\| \leq \mu p(\mathbf{x}\_{1},\mathbf{x}\_{2},t\_{2}) \land \left\| \mathbf{s}(\mathbf{x}\_{1},\mathbf{x}\_{2},t\_{2}) - \mathbf{s}(\mathbf{x}\_{1},\mathbf{x}\_{2},t\_{1}) \right\| = 0 \; \; \; \; \langle i,j \rangle \in \Gamma\_{\mathcal{S}}(t\_{2});\\ \left\| \mathbf{q}(\mathbf{x}\_{1},\mathbf{x}\_{2},t\_{2}) \right\| = \mu p(\mathbf{x}\_{1},\mathbf{x}\_{2},t\_{2}) \land \left\| \mathbf{s}(\mathbf{x}\_{1},\mathbf{x}\_{2},t\_{2}) - \mathbf{s}(\mathbf{x}\_{1},\mathbf{x}\_{2},t\_{1}) \right\| > 0, \; \; \langle i,j \rangle \in \Gamma\_{\mathcal{S}}(t\_{2}) - \Gamma\_{\mathcal{S}}(t\_{2}). \end{cases} \tag{8}$$

Here, *<sup>S</sup>* 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 *A*1 and *A*<sup>2</sup> , respectively, where the bodies deform due shear tractions.

Firstly, the above mentioned axis is assumed to pass through the initial point of contact, therefore 11 21 *A* () () *t At O* . In the considered time frame, all points in the contacting bodies undergo rigid-body translations (1) 1 and (2) <sup>1</sup> along the direction of 1 *x* , and points *A*1 and *A*2 also undergo tangential displacements. If for these points the composite parameters cancel each other, i.e. 12 12 () () *t ut* , with (1) (2) 11 1 and (2) (1) 11 1 *uu u* ; therefore, 12 22 *A* () () *t At O* 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 11 21 *A* () () *t At* , 1 2 *A* ( ) *t* diverge from 2 2 *A* ( ) *t* with a relative (composite) slip distance (1) (2) 11 1 *ss s* . 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.

Numerical Simulation of Slip-Stick Elastic Contact 135

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

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 *<sup>P</sup>* of the contacting bodies needs to be considered, which constitutes an important advantage of SAM over other numerical techniques. In *<sup>P</sup>* , 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 1 2 depends on the

elementary patches) are identified by a pair of indices (,) *i j* , with 1 1 *i N* , 2 1 *i N* . Any continuous distribution 1 2 *f*(,) *x 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 \* \*

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,

To simulate the loading history, additional discretization is performed in the temporal domain, meaning the load is imposed in small steps *k N* 1, . Consequently, *pi jk* (, , ) 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

Digitization of Eqs. (1) - (4) lead to the following numerical model, which will be referred

(,) ( ) ( ) ( , , ); *<sup>C</sup> ij A k W k pi jk* 

the slope of the initial separation should be small everywhere over *<sup>P</sup>* .

constant in the numerical program, e.g. the grid steps.

1 2 *x x*, are the coordinates of the control point of cell (,) *i j* .

, *<sup>i</sup>* 1,2 . The grid control points (centroids of rectangular

(9)

the theory of Digital Signal Processing (DSP) is also discussed.

**3.1. Principles of problem discretization** 

grid steps *<sup>i</sup>* in the direction of *<sup>i</sup> x*

1 2 *fij fx x* (,) ( , ) , where \* \*

from now on as the NC:

**Figure 2.** Geometrical condition of deformation in the tangential direction
