**3. Formulation of discrete slip-stick elastic contact problem**

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

### **3.1. Principles of problem discretization**

134 Numerical Simulation – From Theory to Industry

(1) (2)

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

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.

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

**3. Formulation of discrete slip-stick elastic contact problem** 

The main obstacles in solving the continuous contact problem consist in the fact that neither the contact and stick area, nor the tractions distributions are known in advance. A trial-anderror approach is therefore necessary, but this infer integration of arbitrary functions (contact tractions distributions) over irregular domains (contact or stick area) in displacements computation. As this cannot generally be performed analytically, numerical 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 grid steps *<sup>i</sup>* in the direction of *<sup>i</sup> x* , *<sup>i</sup>* 1,2 . The grid control points (centroids of rectangular 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 \* \* 1 2 *fij fx x* (,) ( , ) , where \* \* 1 2 *x x*, are the coordinates of the control point of cell (,) *i j* .

The limiting surfaces of the contacting bodies are sampled in two height arrays corresponding to grid control points. Such data can be obtained from an optical profilometer or can be generated numerically. The sum of the two heights at node (,) *i j* yields the composite initial surface height *hi i j* (,) . For the half-space approximation to remain valid, the slope of the initial separation should be small everywhere over *<sup>P</sup>* .

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 constant in the numerical program, e.g. the grid steps.

Digitization of Eqs. (1) - (4) lead to the following numerical model, which will be referred from now on as the NC:

$$\mathcal{W}(k) = \Delta \sum\_{(i,j) \in A\_{\mathbb{C}}(k)} p(i, j, k); \tag{9}$$

$$M\_1(k) = \Delta \sum\_{(i,j) \in A\_\mathbb{C}(k)} p(i, j, k) \mathbf{x}\_2(i, j); \quad M\_2(k) = \Delta \sum\_{(i,j) \in A\_\mathbb{C}(k)} p(i, j, k) \mathbf{x}\_1(i, j); \tag{10}$$

$$h(\mathbf{i}, \mathbf{j}, \mathbf{k}) = h\mathbf{i}(\mathbf{i}, \mathbf{j}) + \mu\_3(\mathbf{i}, \mathbf{j}, \mathbf{k}) - \alpha\_3(\mathbf{k}) - \phi\_1(\mathbf{k})\mathbf{x}\_2(\mathbf{i}, \mathbf{j}) - \phi\_2(\mathbf{k})\mathbf{x}\_1(\mathbf{i}, \mathbf{j}), \mathbf{(i}, \mathbf{j}) \in A\_\mathbf{p};$$

$$\begin{cases} p(i,j,k) > 0 \land h(i,j,k) = 0, \quad (i,j) \in A\_{\mathbb{C}}(k);\\ p(i,j,k) = 0 \land h(i,j,k) > 0, \quad (i,j) \in A\_{p} - A\_{\mathbb{C}}(k). \end{cases} \tag{12}$$

Numerical Simulation of Slip-Stick Elastic Contact 137

(18)

1 12 2 12

. (20)

. The influence

(1) 1 (1) 2

(21)

(22)

(23)

*q q*

(1) (1) (1) 32 33

*KK p*

.

cannot be performed analytically except for a few cases; however, its discrete counterpart

(,) ( , ) ( , ),

is the influence coefficient, expressing the displacement in the

2 21 1

 

(,) () , () .

, induced in the cell (,) *i j* of the body ( ) *<sup>n</sup>* by a unit contact traction,

(19)

can be computed numerically for any contact area and any distribution of tractions:

*C*

( ,)

uniformly distributed in the cell (,) *k* , acting along direction of *x*

*x xk n n*

*x xk*

( ) ( )

(2) (1) (2) (2) (2) (2) (1) <sup>3</sup> 3 3 31 32 33 <sup>31</sup>

influence coefficients ( ) *<sup>n</sup> <sup>K</sup>m* can be computed using the following relations:

*nn n nn n m m*

1

*nn n nn n m m*

*u uu KKK p K* 

2 21 1

() 2 () 2

() 2 () 2

equivalent manner using the symbol " " to denote discrete cyclic convolution:

*u uu KKK q K K K u uu KKK q K K K*

where ( )(,) *<sup>n</sup> K i kj* 

*nnn*

direction of *x*

*k A*

( ) ( ) ( )

*n nn*

*u ij K i kj p k*

coefficients yield from integration of Green functions over rectangular elementary patches:

*K i kj G x i x x j x dx dx*

Eq. (18) is in fact a two-dimensional convolution product, and can be written in an

() () () *nnn uK p*

 

The assembly of contributions of all contact tractions to displacements can be expressed as:

(2) (1) (2) (2) (2) (2) (1) (1) (1) 1 1 11 12 13 1 11 12 13 1 (2) (1) (2) (2) (2) (2) (1) (1) (1) 2 2 2 21 22 23 2 21 22 23

> (2) (1) (2) (1) (2) (1) 11 11 12 12 13 13 1 (2) (1) (2) (1) (2) (1) 21 21 22 22 23 23 2 (2) (1) (2) (1) (2) (1) 31 31 32 32 33 33

*KKKKKK q KKKKKK q KKKKKK p*

When the second subscript is omitted, the referred displacement is the sum of all contributions. The positive sign in computation of 3 *u* is related to the fact that normal displacements are computed in coordinate systems linked to each body, having the 3 *x* axes pointing inward. It should also be noted that from pressure definition (2) (1) *ppp* , but (2) (1) , 1,2 *i ii q q qi* , as the mutual actions of two bodies upon each other are equal, but directed to contrary parts. The

() () ()

*nn n m*

*K E ij*

( , ,,)

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

 

> 

, , , ln ln ;

 

( , , ( , ) , ( , ) ) ( , , ( , ) , ( , ) ); 2 2 2 2

() () () 2 2 2 2 33 1 2 ( ) 1 2 12 2 1 12

*<sup>n</sup> k E xx x x xx x x xx*

() () () 1 2 () () () 1 2 1 2 1 2 () () () 1 2 () () () 1 2 1 2 1 2

*k E x ij x ij k E x ij x ij*

<sup>2</sup> ( )

*n*

*E* 

*k E x ij x ij k E x ij x ij*

In a similar way, the continuous model consisting in Eqs. (5) - (8) has its discrete counterpart, referred to as the TC:

$$T\_n(k) = \Delta \sum\_{(i,j)\in A\_\mathbb{C}(k)} q\_n(i,j,k), n = 1,2;\tag{13}$$

$$M\_3(k) = \Delta \sum\_{(i,j) \in A\_{\mathbb{C}}(k)} \left[ q\_2(i,j,k) \mathbf{x}\_1(i,j) - q\_1(i,j,k) \mathbf{x}\_2(i,j) \right] \tag{14}$$

$$
\begin{aligned}
\begin{bmatrix}
s\_1(i,j,k) - s\_1(i,j,k-1) \\
s\_2(i,j,k) - s\_2(i,j,k-1)
\end{bmatrix} &= \begin{bmatrix}
u\_1(i,j,k) - u\_1(i,j,k-1) \\
u\_2(i,j,k) - u\_2(i,j,k-1)
\end{bmatrix} - \begin{bmatrix}
o\_1(k) - \alpha\_1(k-1) \\
o\_2(k) - \alpha\_2(k-1)
\end{bmatrix} - \dots \\
\begin{bmatrix}
\phi\_3(k) - \phi\_3(k-1)
\end{bmatrix} \begin{bmatrix}
x\_2(i,j) \\ x\_1(i,j)
\end{bmatrix}, \quad (i,j) \in A\_\mathbb{C}(k);
\end{aligned} \tag{15}
$$

$$\begin{cases} \left\| \mathbf{q}(i,j,k) \right\| \le \mu p(i,j,k) \land \left\| \mathbf{s}(i,j,k) - \mathbf{s}(i,j,k-1) \right\| = 0, \langle i,j \rangle \in A\_{\mathbb{S}}(k);\\ \left\| \mathbf{q}(i,j,k) \right\| = \mu p(i,j,k) \land \left\| \mathbf{s}(i,j,k) - \mathbf{s}(i,j,k-1) \right\| > 0, \langle i,j \rangle \in A\_{\mathbb{C}}(k) - A\_{\mathbb{S}}(k). \end{cases} \tag{16}$$

In these equations, , *AP C A* and *AS* are the discrete counterparts of *<sup>P</sup>* , *<sup>C</sup>* and *<sup>S</sup>* , respectively, consisting in sets of elementary patches. Consequently, in the numerical formulation, the contact or the stick area can only vary in equal increments . A fine mesh is thus required for accurate estimation of contact domains.

#### **3.2. Numerical computation of displacement fields**

In Contact Mechanics, it is common practice to assimilate the contacting bodies with elastic half-spaces. This assumption holds if the dimensions of contact area are small compared to significant dimensions of contacting bodies, and allows expressing displacements according to superposition principle applied to Green functions for the elastic half-space:

$$\mathbf{u}\_{\vec{\eta}}^{(k)}(\mathbf{x}\_{1},\mathbf{x}\_{2}) = \iint p\_{\vec{\eta}}^{(k)}(\mathbf{x}\_{1}^{\prime},\mathbf{x}\_{2}^{\prime}) \cdot \mathbf{G}\_{\vec{\eta}}^{(k)} \left(\mathbf{x}\_{1} - \mathbf{x}\_{1\prime}^{\prime} \cdot \mathbf{x}\_{2} - \mathbf{x}\_{2}^{\prime}\right) d\mathbf{x}\_{1}^{\prime} d\mathbf{x}\_{2}^{\prime}; \mathbf{i}, \mathbf{j} = 1,2,3; k = 1,2,\tag{17}$$

where ( ) 1 2 , *<sup>k</sup> G xx ij* denotes displacement of a point in body ( *<sup>k</sup>* ), in the direction of *<sup>i</sup> <sup>x</sup>* , induced by a unit point force applied in origin along direction of *<sup>j</sup> x* . The functions *G* for the elastic half-space, also referred to as fundamental solutions or Green functions, were derived by Boussinesq, (Boussinesq, 1969) and Cerruti (Cerruti, 1882). Integral in Eq. (17) cannot be performed analytically except for a few cases; however, its discrete counterpart can be computed numerically for any contact area and any distribution of tractions:

136 Numerical Simulation – From Theory to Industry

counterpart, referred to as the TC:

1 2 (,) ( ) ( ) ( , , ) ( , );

2 1

3 3 12 21 ( , , ) ( , ) ( , , ) ( ) ( ) ( , ) ( ) ( , ), ( , ) ; *<sup>P</sup> h i j k hi i j u i j k k k x i j k x i j i j A* 

In a similar way, the continuous model consisting in Eqs. (5) - (8) has its discrete

( ) ( , , ), 1,2;

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

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

2

1

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

*i jk pi jk i jk i jk i j A k A k*

In these equations, , *AP C A* and *AS* are the discrete counterparts of *<sup>P</sup>* , *<sup>C</sup>* and *<sup>S</sup>* , respectively, consisting in sets of elementary patches. Consequently, in the numerical formulation, the contact or the stick area can only vary in equal increments . A fine mesh

In Contact Mechanics, it is common practice to assimilate the contacting bodies with elastic half-spaces. This assumption holds if the dimensions of contact area are small compared to significant dimensions of contacting bodies, and allows expressing displacements according

1 2 12 1 1 2 2 1 2 (,) (,) , ; , 1,2,3; 1,2,

1 2 , *<sup>k</sup> G xx ij* denotes displacement of a point in body ( *<sup>k</sup>* ), in the direction of *<sup>i</sup> <sup>x</sup>*

the elastic half-space, also referred to as fundamental solutions or Green functions, were derived by Boussinesq, (Boussinesq, 1969) and Cerruti (Cerruti, 1882). Integral in Eq. (17)

(17)

*ij <sup>j</sup> ij u x x p x x G x x x x dx dx i j k*

*x ij*

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

*pi jk hi jk i j A k pi jk hi jk i j A A k*

(,) ( )

1 1 1 1 1 1 2 2 2 2 2 2

*s ijk s ijk u ijk u ijk k k s ijk s ijk u ijk u ijk k k*

*x ij k k*

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

*i jk pi jk i jk i jk i j A k*

(,) ( )

*<sup>C</sup> ij A k*

to superposition principle applied to Green functions for the elastic half-space:

( ) ( ) ( )

induced by a unit point force applied in origin along direction of *<sup>j</sup> x*

 

3 3

(,) ( ) ( 1) (,)

**q s s q s s**

is thus required for accurate estimation of contact domains.

**3.2. Numerical computation of displacement fields** 

*C k kk*

where ( )

*C n n ij A k T k q ijk n* 

3 2112

*M k q ijkx ij q ijkx ij*

 

(,) ( ) ( ) ( , , ) ( , ); *<sup>C</sup> ij A k M k pi jkx i j* 

> 

*C P C*

(13)

(14)

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

 

 

*S*

*C S*

(10)

(11)

(15)

(16)

,

. The functions *G* for

(12)

*<sup>C</sup> ij A k M k pi jkx i j* 

$$\ln\_{\zeta\bar{\xi}}^{(n)}(\mathbf{i},\mathbf{j}) = \sum\_{(k,\ell)\neq A\_{\mathbb{C}}} K^{(n)}\_{\zeta\bar{\xi}}(\mathbf{i}-k,\mathbf{j}-\ell) \cdot \mathbf{p}^{(n)}\_{\bar{\xi}}(\mathbf{k},\ell)\_{\mathbf{\prime}} \tag{18}$$

where ( )(,) *<sup>n</sup> K i kj* is the influence coefficient, expressing the displacement in the direction of *x* , induced in the cell (,) *i j* of the body ( ) *<sup>n</sup>* by a unit contact traction, uniformly distributed in the cell (,) *k* , acting along direction of *x* . The influence coefficients yield from integration of Green functions over rectangular elementary patches:

$$K^{(n)}\_{\\\zeta\bar{\xi}}(\mathbf{i}-\mathbf{k},j-\ell) = \int\_{\substack{\mathbf{x}\_{2}(\ell)-\Lambda\_{2}/2\\\mathbf{x}\_{1}(\ell)-\Lambda\_{2}/2}}^{\mathbf{x}\_{2}(\ell)+\Lambda\_{2}/2} \int\_{\substack{\mathbf{x}\_{1}(k)-\Lambda\_{1}/2\\\mathbf{x}\_{1}(k)-\Lambda\_{1}/2}}^{\mathbf{x}\_{1}(\ell)+\Lambda\_{1}/2} G^{(n)}\_{\zeta\bar{\xi}}\left(\mathbf{x}\_{1}(\ell)-\mathbf{x}\_{1}',\mathbf{x}\_{2}(j)-\mathbf{x}\_{2}'\right) d\mathbf{x}\_{1}'d\mathbf{x}\_{2}'.\tag{19}$$

Eq. (18) is in fact a two-dimensional convolution product, and can be written in an equivalent manner using the symbol " " to denote discrete cyclic convolution:

$$
\mu\_{\zeta\overline{\xi}}^{(n)} = K\_{\zeta\overline{\xi}}^{(n)} \otimes p\_{\zeta}^{(n)}.\tag{20}
$$

The assembly of contributions of all contact tractions to displacements can be expressed as:

$$
\begin{aligned}
\begin{bmatrix} u\_1\\ u\_2\\ u\_3 \end{bmatrix} = \begin{bmatrix} u\_1^{(2)} - u\_1^{(1)}\\ u\_2^{(2)} - u\_2^{(1)}\\ u\_3^{(2)} + u\_3^{(1)} \end{bmatrix} = \begin{bmatrix} K\_{11}^{(2)} & K\_{12}^{(2)} & K\_{13}^{(2)}\\ K\_{21}^{(2)} & K\_{22}^{(2)} & K\_{23}^{(2)}\\ K\_{31}^{(2)} & K\_{32}^{(2)} & K\_{33}^{(2)} \end{bmatrix} \otimes \begin{bmatrix} q\_1^{(2)}\\ q\_2^{(2)}\\ p^{(2)} \end{bmatrix} - \begin{bmatrix} K\_{11}^{(1)} & K\_{12}^{(1)} & K\_{13}^{(1)}\\ K\_{21}^{(1)} & K\_{22}^{(1)} & K\_{23}^{(1)}\\ -K\_{31}^{(1)} & -K\_{32}^{(1)} & -K\_{33}^{(1)} \end{bmatrix} \otimes \begin{bmatrix} q\_1^{(1)}\\ q\_2^{(1)}\\ p^{(1)} \end{bmatrix} = \\ = \begin{bmatrix} K\_{11}^{(2)} + K\_{11}^{(1)} & K\_{12}^{(2)} + K\_{12}^{(1)} & K\_{13}^{(2)} - K\_{13}^{(1)}\\ K\_{21}^{(2)} + K\_{21}^{(2)} & K\_{22}^{(2)} - K\_{23}^{(1)}\\ K\_{31}^{(2)} + K\_{31}^{(2)} + K\_{32}^{(2)} & K\_{33}^{(2)} - K\_{33}^{(2)} \end{bmatrix} \otimes \begin{bmatrix} q\_1\\ q\_1\\ p \end{bmatrix}.
\end{aligned} \tag{21}$$

When the second subscript is omitted, the referred displacement is the sum of all contributions. The positive sign in computation of 3 *u* is related to the fact that normal displacements are computed in coordinate systems linked to each body, having the 3 *x* axes pointing inward. It should also be noted that from pressure definition (2) (1) *ppp* , but (2) (1) , 1,2 *i ii q q qi* , as the mutual actions of two bodies upon each other are equal, but directed to contrary parts. The influence coefficients ( ) *<sup>n</sup> <sup>K</sup>m* can be computed using the following relations:

$$K\_{lm}^{(n)}(\nu^{(n)}, E^{(n)}, i, j) = $$

$$k\_{lm}^{(n)}(\nu^{(n)}, E^{(n)}, x\_1(i, j) + \frac{\Delta\_1}{2}, x\_2(i, j) + \frac{\Delta\_2}{2}) + k\_{lm}^{(n)}(\nu^{(n)}, E^{(n)}, x\_1(i, j) - \frac{\Delta\_1}{2}, x\_2(i, j) - \frac{\Delta\_2}{2}) - \dots \quad \text{(22)}$$

$$k\_{lm}^{(n)}(\nu^{(n)}, E^{(n)}, x\_1(i, j) + \frac{\Delta\_1}{2}, x\_2(i, j) - \frac{\Delta\_2}{2}) - k\_{lm}^{(n)}(\nu^{(n)}, E^{(n)}, x\_1(i, j) - \frac{\Delta\_1}{2}, x\_2(i, j) + \frac{\Delta\_2}{2});$$

$$k\_{33}^{(n)}(\nu^{(n)}, E^{(n)}, x\_1, x\_2) = \frac{1 - \left(\nu^{(n)}\right)^2}{\pi E^{(n)}} \left[x\_1 \ln\left(x\_2 + \sqrt{x\_1^2 + x\_2^2}\right) + x\_2 \ln\left(x\_1 + \sqrt{x\_1^2 + x\_2^2}\right)\right];\tag{23}$$

$$k\_{13}^{(n)}\left(\nu^{(n)}, E^{(n)}, \mathbf{x}\_1, \mathbf{x}\_2\right) = \frac{(1+\nu^{(n)})(1-2\nu^{(n)})}{2\pi E^{(n)}} \left[2\mathbf{x}\_1 \tan^{-1}\left(\frac{\sqrt{\mathbf{x}\_1^2 + \mathbf{x}\_2^2} - \mathbf{x}\_2}{\mathbf{x}\_1}\right) - \mathbf{x}\_2 \ln\left(\sqrt{\mathbf{x}\_1^2 + \mathbf{x}\_2^2}\right)\right];\tag{24}$$

$$k\_{11}^{(n)}\left(\nu^{(n)}, E^{(n)}, \mathbf{x}\_1, \mathbf{x}\_2\right) = \frac{1 + \nu^{(n)}}{\pi E^{(n)}} \left[ \mathbf{x}\_2 \ln\left(\mathbf{x}\_1 + \sqrt{\mathbf{x}\_1^2 + \mathbf{x}\_2^2}\right) - \mathbf{x}\_2 + (1 - \nu^{(n)}) \mathbf{x}\_1 \ln\left(\mathbf{x}\_2 + \sqrt{\mathbf{x}\_1^2 + \mathbf{x}\_2^2}\right) \right];\tag{25}$$

$$k\_{21}^{(n)}\left(\nu^{(n)}, E^{(n)}, \mathbf{x}\_1, \mathbf{x}\_2\right) = -\frac{\nu^{(n)}(1+\nu^{(n)})}{\pi E^{(n)}} \sqrt{\mathbf{x}\_1^2 + \mathbf{x}\_2^2};\tag{26}$$

Numerical Simulation of Slip-Stick Elastic Contact 139

from Eqs. (11) and (15) respectively, is to be solved, having as unknowns the normal or the tangential tractions. The main difficulties consist in the fact that the systems are essentially non-linear, due to presence of terms containing rigid-body motions, and their size, related to the number of cells in *AC* and *AS* respectively, is not known in advance. Moreover, resolution of a linear system with *N* unknowns has an order of computation of <sup>3</sup> *O N*( ) when direct methods like Gaussian elimination are used. In contact problems, grids with <sup>3</sup> *N* 10 nodes on the contact area are usually considered for nominal contact surfaces, in order to get an accurate estimation of contact area extents and to minimize the discretization error (i.e. the error introduced by digitization of problem parameters). When deterministic rough surfaces are studied, a minimum of <sup>6</sup> *N* 10 points is required to reproduce the detailed features of rough contact behaviour. Consequently, only fast numerical methods, based on iterative approaches, are suitable for this type of problem. According to a review paper by Allwood (Allwood, 2005), the Conjugate Gradient Method (CGM) has the fastest convergence. Implementation of CGM to this type of problem is authorised by the fact that the influence coefficients matrix is symmetrical and positive definite, therefore convergence

The algorithm used in this study is based on the modified CGM advanced by Polonsky and Keer (Polonsky & Keer, 1999) for the study of frictionless rough contact problems under normal loading. Existence and uniqueness of the solution is discussed by Kalker and van Randen (Kalker & van Randen, 1972). This algorithm was later enhanced by Spinu and Diaconescu (Spinu & Diaconescu, 2008) to allow for a bending moment on the contact area, and it was recently applied by Spinu and Frunza (Spinu & Frunza, 2011a) to solve

Eq. (21) suggests that, in case of dissimilarly elastic materials (i.e. (1) (2) *E E* and / or

It should be noted that in some cases, matching names will be used for variables having the same role in the algorithm structure, although their content might be different in the NC from that in the TC. Let us assume that a new loading increment *k* is applied. In both NC and TC, all parameters are specific to the achieved loading level; therefore, the third formal parameter, related to the loading history, will be omitted for brevity. An additional symbol

 ), computation of displacements in either direction require the knowledge of contact tractions on every direction. From this yields the coupling between the NC and the TC. However, in case of similarly elastic materials, Eq. (32) proves that normal displacements are independent of shear tractions, and tangential displacements are independent of pressure. In the latter case, solution of NC is decoupled from that of TC; however, pressure distribution is still needed in the TC for assessment of shear tractions in the slip zone, according to Coulomb's law. In the framework developed herein, we refer to the solution of uncoupled NC or TC as being the solution of NC when shear tractions are known, but otherwise arbitrarily distributed, as well as the solution of TC when pressure is

is assured.

(1) (2) 

numerically the Cattaneo-Mindlin problem.

known, but otherwise arbitrarily distributed.

**4.1. Numerical solution of uncoupled contact problems** 

$$k\_{23}^{(n)}\left(\nu^{(n)}, E^{(n)}, \mathbf{x}\_1, \mathbf{x}\_2\right) = k\_{13}^{(n)}\left(\nu^{(n)}, E^{(n)}, \mathbf{x}\_2, \mathbf{x}\_1\right);\tag{27}$$

$$\left(k\_{31}^{(n)}\left(\nu^{(n)}, E^{(n)}, \mathbf{x}\_1, \mathbf{x}\_2\right) = -k\_{13}^{(n)}\left(\nu^{(n)}, E^{(n)}, \mathbf{x}\_1, \mathbf{x}\_2\right);\tag{28}$$

$$k\_{12}^{(n)}\left(\nu^{(n)}, E^{(n)}, \mathbf{x}\_1, \mathbf{x}\_2\right) = k\_{21}^{(n)}\left(\nu^{(n)}, E^{(n)}, \mathbf{x}\_1, \mathbf{x}\_2\right);\tag{29}$$

$$k\_{22}^{(n)}\left(\nu^{(n)}, E^{(n)}, \mathbf{x}\_1, \mathbf{x}\_2\right) = k\_{11}^{(n)}\left(\nu^{(n)}, E^{(n)}, \mathbf{x}\_2, \mathbf{x}\_1\right);\tag{30}$$

$$k\_{32}^{(n)}\left(\nu^{(n)}, E^{(n)}, \mathbf{x}\_1, \mathbf{x}\_2\right) = -k\_{23}^{(n)}\left(\nu^{(n)}, E^{(n)}, \mathbf{x}\_{2'}, \mathbf{x}\_1\right),\tag{31}$$

where ( ) *<sup>n</sup> E* and ( ) *<sup>n</sup>* are the elastic constants (Young modulus and Poisson's ratio, respectively) of materials of the two contacting bodies. If (1) (2) *EEE* and (1) (2) , then (1) (2) , , 1,2,3, *K K K ij ij ij ij* and Eq. (21) takes a simplified form:

$$
\begin{bmatrix} u\_1 \\ u\_2 \\ u\_3 \end{bmatrix} = 2 \begin{bmatrix} \mathcal{K}\_{11} & \mathcal{K}\_{12} & \mathcal{0} \\ \mathcal{K}\_{21} & \mathcal{K}\_{22} & \mathcal{0} \\ \mathcal{0} & \mathcal{0} & \mathcal{K}\_{33} \end{bmatrix} \otimes \begin{bmatrix} q\_1 \\ q\_2 \\ p \end{bmatrix}. \tag{32}
$$

The most efficient way to compute the two-dimensional convolution products in Eqs. (21) or (32) is by using spectral methods. According to convolution theorem, the convolution of two signals, each having *N* samples, requires <sup>2</sup> *O N*( ) operations in time/space domain, but only *ON N* ( log ) operations in the frequency domain, where it reduces to element-wise product.

However, when transferring to frequency domain via fast Fourier transform (FFT), presumption of signal periodicity is automatically assumed, which leads to contamination of convolution result due to spurious neighbouring periods. Liu & al. (Liu & al., 2000) advanced a fast and robust algorithm to circumvent the periodicity error, which is also applied here. This algorithm is able to assess linear convolution by computing the discrete cyclic convolution of the two terms on an extended domain, with a special, "warp-around" arrangement of influence coefficients. The base for this approach is discussed in detail in (Press et al., 1992).

### **4. Algorithm description**

Considering the similarity between NC and TC, the same type of algorithm could be used to solve either model considered independently. Indeed, in both cases, a linear system arising from Eqs. (11) and (15) respectively, is to be solved, having as unknowns the normal or the tangential tractions. The main difficulties consist in the fact that the systems are essentially non-linear, due to presence of terms containing rigid-body motions, and their size, related to the number of cells in *AC* and *AS* respectively, is not known in advance. Moreover, resolution of a linear system with *N* unknowns has an order of computation of <sup>3</sup> *O N*( ) when direct methods like Gaussian elimination are used. In contact problems, grids with <sup>3</sup> *N* 10 nodes on the contact area are usually considered for nominal contact surfaces, in order to get an accurate estimation of contact area extents and to minimize the discretization error (i.e. the error introduced by digitization of problem parameters). When deterministic rough surfaces are studied, a minimum of <sup>6</sup> *N* 10 points is required to reproduce the detailed features of rough contact behaviour. Consequently, only fast numerical methods, based on iterative approaches, are suitable for this type of problem. According to a review paper by Allwood (Allwood, 2005), the Conjugate Gradient Method (CGM) has the fastest convergence. Implementation of CGM to this type of problem is authorised by the fact that the influence coefficients matrix is symmetrical and positive definite, therefore convergence is assured.

The algorithm used in this study is based on the modified CGM advanced by Polonsky and Keer (Polonsky & Keer, 1999) for the study of frictionless rough contact problems under normal loading. Existence and uniqueness of the solution is discussed by Kalker and van Randen (Kalker & van Randen, 1972). This algorithm was later enhanced by Spinu and Diaconescu (Spinu & Diaconescu, 2008) to allow for a bending moment on the contact area, and it was recently applied by Spinu and Frunza (Spinu & Frunza, 2011a) to solve numerically the Cattaneo-Mindlin problem.
