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

138 Numerical Simulation – From Theory to Industry

where ( ) *<sup>n</sup> E* and ( ) *<sup>n</sup>*

**4. Algorithm description** 

*n n nnn*

 () () 2 2 () () () 1 122 2 2 13 1 2 ( ) 1 2 12

 ( ) () () () 2 2 ( ) 2 2 11 1 2 ( ) 2 1 12 2 1 2 12

> () () () () () 2 2 21 1 2 ( ) 1 2 (1 ) , , , ;

 () () () () () () <sup>23</sup> 1 2 13 2 1 , , , , ,, ; *nnn nnn k E xx k E xx*

 () () () () () () <sup>31</sup> 1 2 13 1 2 , , , , ,, ; *nnn nnn k E xx k E xx*

 () () () () () () <sup>12</sup> 1 2 21 1 2 , , , , ,, ; *nnn nnn k E xx k E xx*

 () () () () () () <sup>22</sup> 1 2 11 2 1 , , , , ,, ; *nnn nnn k E xx k E xx*

 () () () () () () <sup>32</sup> 1 2 23 2 1 , , , , ,, , *nnn nnn k E xx k E xx*

> 1 11 12 1 2 21 22 2

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).

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

*u KK q u KK q u Kp*

respectively) of materials of the two contacting bodies. If (1) (2) *EEE* and (1) (2)

3 33

0 0

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

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

*k E xx x x xx E x*

2

*E* 

then (1) (2) , , 1,2,3, *K K K ij ij ij ij* and Eq. (21) takes a simplified form:

*n*

*<sup>n</sup> nnn <sup>n</sup>*

 

(1 )(1 2 ) , , , 2 tan ln ;

<sup>1</sup> , , , ln (1 ) ln ;

*E*

 

> 

 

 

> 

0 2 0.

are the elastic constants (Young modulus and Poisson's ratio,

 

1

 

(26)

(27)

(28)

(29)

(30)

(31)

(32)

 

 ,

(25)

*xxx*

(24)

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

It should be noted that in some cases, matching names will be used for variables having the same role in the algorithm structure, although their content might be different in the NC from that in the TC. Let us assume that a new loading increment *k* is applied. In both NC and TC, all parameters are specific to the achieved loading level; therefore, the third formal parameter, related to the loading history, will be omitted for brevity. An additional symbol " " is used to denote the increment of the referred quantity after application of the *k*-th loading increment, i.e. 11 1 *u ij u ijk u ijk* ( , ) ( , , ) ( , , 1) .

The contact tractions corresponding to the new loading level *p k*( ) or **q**( ) *k* can be assessed using the same type of algorithm, consisting in the sequences described in the following paragraphs.

1. Initialize auxiliary variables:

$$\begin{aligned} \text{NC}: \quad \theta = 0; \quad \mathcal{R}\_{\text{old}} = 1; d\_3(i, j) = 0, p\_{\text{old}}(i, j) = 0, (i, j) \in A\_p; \\ \text{TC}: \quad \theta = 0; \mathcal{R}\_{\text{old}} = 1; d\_1(i, j) = 0, d\_2(i, j) = 0, \mathbf{q}\_{\text{old}}(i, j) = \mathbf{0}, (i, j) \in A\_p. \end{aligned} \tag{33}$$

Numerical Simulation of Slip-Stick Elastic Contact 141

*C*

*S*

2

2 2

.

(39)

 

 *u ij u ij* ( , ), ( , ) , is

(37)

3. Compute the relative displacement field over the grid domain using Eq. (21). It should be remembered that in the NC, pressure was adopted in step 2 and shear tractions are assumed as known, but arbitrarily distributed, while in the TC, shear tractions were adopted in step 2, and pressure is assumed known but arbitrarily distributed. It should

also be noted that in case of TC, the increment of displacement field, 1 2

3 3 12 21

*NC hi i j u i j x i j x i j i j A*

3

(,) (,) 0 : ,(,) , (,) (,) <sup>0</sup>

*u ij x ij TC ij A u ij x ij*

1 12

 

 

(,)

*ij A*

(,)

2

3

*ij A*

*S*

*C*

2 21

312 3 3 12 21

*NC hi i j u i j x i j x i j*

(,) <sup>3</sup> 1

**Μ**

1

**Μ**

(,)

*ij A*

*ij A*

(,)

(,) 1

*TC u ij*

(,)

*ij A*

2 2

*S*

: (,)

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

4. Estimate the rigid-body motions in order to linearize Eqs. (11) and (15). According to complementarity conditions in Eqs. (12) and (16), the following relations hold:

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

resulting in linear systems having a number of equations equal the number of cells in the contact and in the stick area, respectively. When convergence is reached, the equations in every of these two systems are not all independent; however, during iterations, they appear to be independent. To get the best possible estimates of rigidbody motions, the systems are treated as over-determined and the best fit is sought using least square minimization. The functional is defined as the square sum of the

1 23 1 1 32 2 2 31

*TC u ij x ij u ij x ij*

1 32

*ij A*

*C*

*NC u i j hi i j x i j*

*C*

: (,) (,) (,) ;

*C*

 

the partial derivatives of to zero, leading to the following solution:

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

The rigid-body motions that minimize the goal function in Eq. (38) are found by setting

3

*ij A*

(,)

*i j*

3 1

*u i j hi i j x i j*

*S*

*S*

*A*

 

(,) (,) (,)

(,) (,)

*u i j hi i j*

1

12 21

*u ij*

(,)

(,) (,) (,) (,)

 *u ijx ij u ijx ij*

(38)

 

 

needed.

residuals:

 

Here, , 1,2,3 *<sup>i</sup> d i* is the initial descend direction in the CGM. The role of the other variables will be discussed later.

2. Adopt the initial guess for contact tractions, using the static force equilibrium. In a first approximation, it will be considered in the NC that all cells are in contact, i.e. *AC P A* , and in the TC that all cells in the contact area are in stick, i.e. *AS C A* . The contact tractions are sought in the following form:

$$\begin{aligned} \text{NC}: \quad &p(i,j) = a\_3 + a\_1 \mathbf{x}\_2(i,j) + a\_2 \mathbf{x}\_1(i,j), \quad (i,j) \in A\_P; \\ \text{TC}: \quad &\begin{cases} q\_1(i,j) = a\_3 \mathbf{x}\_2(i,j) + a\_1; \\ q\_2(i,j) = a\_3 \mathbf{x}\_1(i,j) + a\_2; \end{cases} \quad (i,j) \in A\_C. \end{aligned} \tag{34}$$

where the unknown coefficients ( ), 1,2,3 *<sup>i</sup> aki* are computed by plugging Eq. (34) into Eqs. (9), (10) and (13) , (14), respectively, yielding:

$$\begin{aligned} \text{NC}: \quad \begin{bmatrix} a\_3 & a\_1 & a\_2 \end{bmatrix}^T &= \mathbf{M}^{-1} \cdot \begin{bmatrix} \mathbf{W} & \mathbf{M}\_1 & \mathbf{M}\_2 \end{bmatrix}^T \Big/ \Delta; \\ \text{TC}: \quad \begin{bmatrix} a\_1 & a\_2 & a\_3 \end{bmatrix}^T &= \mathbf{M}^{-1} \cdot \begin{bmatrix} T\_1 & T\_2 & \mathbf{M}\_3 \end{bmatrix}^T \Big/ \Delta; \end{aligned} \tag{35}$$

2 1 (,) (,) (,) 2 2 2 2 1 (,) (,) (,) 2 1 1 2 1 (,) (,) (,) (,) ( , 1 (,) (,) : (,) (,) (,) (,) ; (,) (,) (,) (,) 1 0 : *C C C CCC C C C S ij A ij A ij A ij A ij A ij A ij A ij A ij A ij A i j x ij x ij NC x ij x ij x ijx ij x ij x ijx ij x ij TC* **Μ M** 2 ) (,) 1 (,) (,) (,) 2 2 2 1 12 (,) (,) (,) (,) 0 1 (,) . (,) (,) (,) (,) *S S SS S S SS A ij A ij A ij A ij A ij A ij A ij A x ij x ij x ij x ij x ij x ij* (36)

3. Compute the relative displacement field over the grid domain using Eq. (21). It should be remembered that in the NC, pressure was adopted in step 2 and shear tractions are assumed as known, but arbitrarily distributed, while in the TC, shear tractions were adopted in step 2, and pressure is assumed known but arbitrarily distributed. It should also be noted that in case of TC, the increment of displacement field, 1 2 *u ij u ij* ( , ), ( , ) , is needed.

140 Numerical Simulation – From Theory to Industry

1. Initialize auxiliary variables:

will be discussed later.

loading increment, i.e. 11 1 

tractions are sought in the following form:

Eqs. (9), (10) and (13) , (14), respectively, yielding:

" is used to denote the increment of the referred quantity after application of the *k*-th

The contact tractions corresponding to the new loading level *p k*( ) or **q**( ) *k* can be assessed using the same type of algorithm, consisting in the sequences described in the following

*u ij u ijk u ijk* ( , ) ( , , ) ( , , 1) .

3

1 2

*NC R d i j p i j i j A*

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

: 0; 1; ( , ) 0, ( , ) 0, ( , ) ,( , ) . *old old P*

*TC R d i j d i j i j i j A*

Here, , 1,2,3 *<sup>i</sup> d i* is the initial descend direction in the CGM. The role of the other variables

2. Adopt the initial guess for contact tractions, using the static force equilibrium. In a first approximation, it will be considered in the NC that all cells are in contact, i.e. *AC P A* , and in the TC that all cells in the contact area are in stick, i.e. *AS C A* . The contact

3 12 21

*NC p i j a a x i j a x i j i j A*

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

where the unknown coefficients ( ), 1,2,3 *<sup>i</sup> aki* are computed by plugging Eq. (34) into

1 312 1 2 1 123 12 3

*T T*

1 (,) (,)

(,) (,) (,) (,)

*x ij x ijx ij x ij*

*T T*

1 1 2 1

) (,)

0 1 (,) .

*S S*

*A ij A*

2 1 12

*x ij x ij x ij x ij*

(,) (,) (,) (,)

: ;

**Μ**

**Μ**

: ;

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

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

*ij A ij A ij A*

*C C C*

 

*CCC*

*C C C*

(,) (,) (,)

(,) (,) (,)

*ij A ij A ij A*

*S SS*

1 0

(,) (,) (,)

*ij A ij A ij A*

(,) ( ,

*ij A i j*

*S*

:

**M**

**Μ**

*TC*

(,) (,) (,)

*ij A ij A ij A*

*SS S*

*ij A ij A ij A*

*NC x ij x ij x ijx ij*

 

1 32 1 2 31 2

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

*q ij ax ij a TC ij A q ij ax ij a*

*NC a a a WM M*

*TC a a a TT M*

*old old P*

**q 0** (33)

*C*

2 1

*x ij x ij*

*P*

2

2

*x ij*

*x ij*

(,)

1

2 2

(34)

(35)

(36)

" 

paragraphs.

4. Estimate the rigid-body motions in order to linearize Eqs. (11) and (15). According to complementarity conditions in Eqs. (12) and (16), the following relations hold:

$$\begin{aligned} \text{NC}: \quad &hi(i,j) + u\_3(i,j) - \alpha\_3 - \phi\_1 \mathbf{x}\_2(i,j) - \phi\_2 \mathbf{x}\_1(i,j) = 0, (i,j) \in A\_\mathbb{C};\\ \text{TC}: \quad &\begin{bmatrix} \delta u\_1(i,j) \\ \delta u\_2(i,j) \end{bmatrix} - \begin{bmatrix} \delta \phi\_1 \\ \delta \phi\_2 \end{bmatrix} + \delta \phi\_3 \begin{bmatrix} \mathbf{x}\_2(i,j) \\ \mathbf{x}\_1(i,j) \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}, (i,j) \in A\_{S'} \end{aligned} \tag{37}$$

resulting in linear systems having a number of equations equal the number of cells in the contact and in the stick area, respectively. When convergence is reached, the equations in every of these two systems are not all independent; however, during iterations, they appear to be independent. To get the best possible estimates of rigidbody motions, the systems are treated as over-determined and the best fit is sought using least square minimization. The functional is defined as the square sum of the residuals:

$$\begin{split} \text{NC}: & \mathcal{N}(a\_{3},\phi\_{1},\phi\_{2}) = \sum\_{(i,j) \in A\_{C}} \left( h\mathbf{i}(\mathbf{i},j) + u\_{3}(\mathbf{i},j) - \alpha\_{3} - \phi\_{1}\mathbf{x}\_{2}(\mathbf{i},j) - \phi\_{2}\mathbf{x}\_{1}(\mathbf{i},j) \right)^{2}; \\ & \text{TC}: & \mathcal{N}(\delta a\_{1},\delta a\_{2},\delta \phi\_{3}) = \sum\_{(i,j) \in A\_{S}} \left( \left( \delta u\_{1}(\mathbf{i},j) - \delta a\_{1} - \delta \phi\_{3}\mathbf{x}\_{2}(\mathbf{i},j) \right)^{2} + \left( \delta u\_{2}(\mathbf{i},j) - \delta a\_{2} - \delta \phi\_{3}\mathbf{x}\_{1}(\mathbf{i},j) \right)^{2} \right). \end{split} \tag{38}$$

The rigid-body motions that minimize the goal function in Eq. (38) are found by setting the partial derivatives of to zero, leading to the following solution:

$$\begin{aligned} \text{NC}: \quad \begin{bmatrix} \alpha\_{3} \\ \phi\_{1} \\ \phi\_{2} \end{bmatrix} = \mathbf{M}^{-1} \cdot \begin{bmatrix} \sum\_{(i,j) \in A\_{C}} \left[ u\_{3}(i,j) + h\dot{u}(i,j) \right] \\ \sum\_{(i,j) \in A\_{C}} \left[ \left( u\_{4}(i,j) + h\dot{u}(i,j) \right) x\_{2}(i,j) \right] \\ \sum\_{(i,j) \in A\_{C}} \left[ \left( u\_{3}(i,j) + h\dot{u}(i,j) \right) x\_{1}(i,j) \right] \end{bmatrix}; \\\ \text{TC}: \begin{bmatrix} \delta\alpha\_{1} \\ \delta\alpha\_{2} \\ \delta\phi\_{3} \end{bmatrix} = \mathbf{M}^{-1} \cdot \begin{bmatrix} \sum\_{(i,j) \in A\_{S}} \delta u\_{1}(i,j) \\ \sum\_{(i,j) \in A\_{S}} \delta u\_{2}(i,j) \\ \sum\_{(i,j) \in A\_{S}} \delta u\_{2}(i,j) \\ \sum\_{(i,j) \in A\_{S}} \left[ \delta u\_{1}(i,j) x\_{2}(i,j) + \delta u\_{2}(i,j) x\_{1}(i,j) \right] \end{bmatrix}. \end{aligned} \tag{39}$$

With this approximation, Eqs. (11) and (15) finally reduce to linear systems having the nodal contact tractions as unknowns. The systems sizes are given by the extents of the contact and of the stick area, respectively, both a priori unknown. A trial-and-error approach is used, in which the systems are solved on successive contact and stick areas, until static equilibrium equations and complementarity conditions are fully satisfied for every cell in the grid.

5. Start the conjugate gradient loop by computing the residuals *r ij* ( , ), 1,2,3 :

$$\begin{aligned} \text{NC}: \quad r\_3(i, j) &= \text{hi}(i, j) + u\_3(i, j) - \alpha\_3 - \phi\_1 \mathbf{x}\_2(i, j) - \phi\_2 \mathbf{x}\_1(i, j), (i, j) \in A\_P; \\ \text{TC}: \quad \begin{bmatrix} r\_1(i, j) \\ r\_2(i, j) \end{bmatrix} &= \begin{bmatrix} \delta u\_1(i, j) \\ \delta u\_2(i, j) \end{bmatrix} - \begin{bmatrix} \delta \alpha\_1 \\ \delta \alpha\_2 \end{bmatrix} + \delta \phi\_3 \begin{bmatrix} \mathbf{x}\_2(i, j) \\ \mathbf{x}\_1(i, j) \end{bmatrix}, \text{ (i, j)} \in A\_{S'} \end{aligned} \tag{40}$$

and the square sum of the residuals on the indicated domains:

$$\begin{aligned} \text{NC}: \quad & R = \sum\_{(i,j) \in A\_{\mathbb{C}}} \left( r\_3(i,j) \right)^2; \\ \text{TC}: \quad & R = \sum\_{(i,j) \in A\_{\mathbb{S}}} \left[ \left( r\_1(i,j) \right)^2 + \left( r\_2(i,j) \right)^2 \right]. \end{aligned} \tag{41}$$

Numerical Simulation of Slip-Stick Elastic Contact 143

*S*

*C*

(44)

, computed using Eq. (43),

(45)

(46)

(47)

 

(2) (1) 3 33 33 3

: ;

*NC c K K d*

*TC*

9. Update the system solution by making a step of length

along direction , 1,2,3 *<sup>i</sup> d i* , derived in Eq. (42):

computation: , *old p p* **q q** *old* .

times, i.e.

area, auxiliary variable

*rij* (,) .

*C C C C*

*S S*

*S S*

stick region, yielding an adjusted stick area:

*A A ij ij ij*

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

: .

8. Memorize the current contact tractions in a new variable for subsequent error

3

 

11 1 22 2

zone to the non-contact area, in order to obey the non-adhesion hypothesis.

*A A i j pi j i j i j pi j*

**q s**

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

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

*A A i j pi j pi j hi j i j i j pi j hi j* 

In the TC, cells for which Coulomb's friction law is not verified are removed from the stick region, and the corresponding shear tractions are set to the value of limiting friction. Also, cells having micro-slip **s**(,) *i j* not opposite to shear tractions **q**(,) *i j* pass from the slip to the

*i j A A i j i j pi j i j i j i j pi j i j*

In either the NC or the TC, if any cell is removed from or re-enters the contact or the stick

resetting the conjugate directions *<sup>i</sup> d* once new cells enter or leave the system. Indeed, these

**<sup>q</sup> q q q**

*NC p i j p i j d i j i j A*

: ( , ) ( , ) ( , ),( , ) ; (,) (,) (,) : ,( , ) . (,) (,) (,)

*q ij q ij d ij TC ij A q ij q ij d ij*

10. Impose complementarity conditions to adjust the size of the system. In the NC, the contact or non-contact status of every cell in the contact domain is to be determined, while in the TC, the stick or slip regime of every cell in the contact area must be assessed.

In the NC, cells having negative pressure are removed from the current contact area, and the corresponding nodal pressures are set to zero. It should be remembered that the current model cannot incorporate adhesion; therefore, normal contact tractions can only be compressive, not tensile. This assumption leads to specific cells passing from the contact

On the other hand, as bodies are considered impenetrable in the frame of Linear Theory of Elasticity, specific cells having vanishing pressure but negative gap *hi j* (,) (which equals the residual), are reinstated in the current contact area. To this end, the corresponding nodal pressures are adjusted with a quantity proven (Polonsky & Keer, 1999) to be positive at all

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

is set to zero, otherwise it is set to unity. This variable allows

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

*c d KKKK*

*c KKKK d*

6. Compute the descend direction , 1,2,3 *<sup>i</sup> d i* in the CGM algorithm. In the multidimensional space of contact tractions, every new descend direction is constructed from the residual to be **K** - orthogonal (Shewchuk, 1994) to all previous residuals and search directions:

$$\begin{aligned} N\mathbf{C}: \quad d\_3(i,j) &\leftarrow \begin{cases} r\_3(i,j) - \left(R\theta \Big/ R\_{old}\right) d\_3(i,j), & (i,j) \in A\_{\curlyvee}; \\ 0, & (i,j) \in A\_{\lozenge} - A\_{\curlyvee}; \end{cases} \\ TC: \quad \begin{bmatrix} d\_1(i,j) \\ d\_2(i,j) \end{bmatrix} &\leftarrow \begin{cases} \left[r\_1(i,j)\right] \\ r\_2(i,j) \end{cases} + \frac{R\theta}{R\_{old}} \begin{bmatrix} d\_1(i,j) \\ d\_2(i,j) \end{bmatrix}, & (i,j) \in A\_{\lozenge}; \\ \left[\left[0\quad 0\right]\right]^T, (i,j) \in A\_{\lozenge} - A\_{\lozenge}. \end{aligned} \tag{42}$$

and overwrite the content of *Rold* : *R R old* .

7. Assess the length of the step to be made in the computed descend direction:

$$\begin{aligned} \text{NC}: \quad \alpha &= \sum\_{\langle i,j\rangle \in A\_{\mathbb{C}}} r\_3(i,j) d\_3(i,j) \Big/ \sum\_{\langle i,j\rangle \in A\_{\mathbb{C}}} c\_3(i,j) d\_3(i,j); \\ \sum\_{\mathbf{C}} r\_1(i,j) d\_1(i,j) &+ \sum\_{\langle i,j\rangle \in A\_{\mathbb{S}}} r\_2(i,j) d\_2(i,j) \\ \text{TC}: \quad \alpha &= \frac{\langle i,j\rangle \in A\_{\mathbb{S}}}{\sum\_{\langle i,j\rangle \in A\_{\mathbb{S}}} c\_1(i,j) d\_1(i,j) + \sum\_{\langle i,j\rangle \in A\_{\mathbb{S}}} c\_2(i,j) d\_2(i,j)}{\sum\_{\langle i,j\rangle \in A\_{\mathbb{S}}} c\_2(i,j) d\_2(i,j)}. \end{aligned} \tag{43}$$

where , 1,2,3 *<sup>i</sup> c i* denotes the convolution of the influence coefficients matrix with the descend direction:

$$\begin{aligned} \text{NC}: \quad c\_3 &= \left| K\_{33}^{(2)} + K\_{33}^{(1)} \right| \otimes d\_3; \\ \text{TC}: \quad \begin{bmatrix} c\_1 \\ c\_2 \end{bmatrix} &= \begin{bmatrix} K\_{11}^{(2)} + K\_{11}^{(1)} & K\_{12}^{(2)} + K\_{12}^{(1)} \\ K\_{21}^{(2)} + K\_{21}^{(1)} & K\_{22}^{(2)} + K\_{22}^{(1)} \end{bmatrix} \otimes \begin{bmatrix} d\_1 \\ d\_2 \end{bmatrix}. \end{aligned} \tag{44}$$

8. Memorize the current contact tractions in a new variable for subsequent error computation: , *old p p* **q q** *old* .

142 Numerical Simulation – From Theory to Industry

search directions:

With this approximation, Eqs. (11) and (15) finally reduce to linear systems having the nodal contact tractions as unknowns. The systems sizes are given by the extents of the contact and of the stick area, respectively, both a priori unknown. A trial-and-error approach is used, in which the systems are solved on successive contact and stick areas, until static equilibrium equations and complementarity conditions are fully satisfied for every cell in the grid.

5. Start the conjugate gradient loop by computing the residuals *r ij* ( , ), 1,2,3 :

3 3 3 12 21

*r ij u ij x ij TC ij A r ij u ij x ij*

 

: (,) ;

*C*

*TC R r i j r i j*

*S*

 

(,) (,) (,) : ,(,) , (,) (,) (,)

 

3

2

1 2

2 2

*old C*

*C S*

*S*

(42)

(43)

(41)

 

*P*

(40)

*S*

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

: (,) (,) .

6. Compute the descend direction , 1,2,3 *<sup>i</sup> d i* in the CGM algorithm. In the multidimensional space of contact tractions, every new descend direction is constructed from the residual to be **K** - orthogonal (Shewchuk, 1994) to all previous residuals and

3 3

*r ij R R d ij ij A NC d i j ij A A*

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

1 1

*r ij d ij R d ij ij A*

 

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

2 2

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

*C C*

: , (,) (,) (,) (,)

where , 1,2,3 *<sup>i</sup> c i* denotes the convolution of the influence coefficients matrix with the

*S S*

*S S*

*T*

*TC r ij d ij R*

7. Assess the length of the step to be made in the computed descend direction:

(,) (,)

*ij A ij A*

*NC r ijd ij c ijd ij*

*ij A ij A*

*ij A ij A*

(,) (,)

(,) (,)

0 0 ,( , ) ,

3 3 3 3

11 22

(,) (,) (,) (,)

*r ijd ij r ijd ij*

*c ijd ij c i jd ij*

1 1 2 2

*ij A A*

*old*

*P C*

3

*NC r i j hi i j u i j x i j x i j i j A*

1 1 12

3

1

2

and overwrite the content of *Rold* : *R R old* .

*TC*

descend direction:

*d ij*

(,)

2 2 21

(,)

*NC R r i j*

*ij A*

(,)

*ij A*

and the square sum of the residuals on the indicated domains:

9. Update the system solution by making a step of length , computed using Eq. (43), along direction , 1,2,3 *<sup>i</sup> d i* , derived in Eq. (42):

$$\begin{aligned} \text{NC}: \quad &p(i,j) \leftarrow p(i,j) - \alpha d\_3(i,j), (i,j) \in A\_{\text{C}};\\ \text{TC}: \quad &\begin{bmatrix} q\_1(i,j) \\ q\_2(i,j) \end{bmatrix} \leftarrow \begin{bmatrix} q\_1(i,j) \\ q\_2(i,j) \end{bmatrix} - \alpha \begin{bmatrix} d\_1(i,j) \\ d\_2(i,j) \end{bmatrix} \mathbf{(}i,j) \in A\_{\text{S}}.\end{aligned} \tag{45}$$

10. Impose complementarity conditions to adjust the size of the system. In the NC, the contact or non-contact status of every cell in the contact domain is to be determined, while in the TC, the stick or slip regime of every cell in the contact area must be assessed.

In the NC, cells having negative pressure are removed from the current contact area, and the corresponding nodal pressures are set to zero. It should be remembered that the current model cannot incorporate adhesion; therefore, normal contact tractions can only be compressive, not tensile. This assumption leads to specific cells passing from the contact zone to the non-contact area, in order to obey the non-adhesion hypothesis.

On the other hand, as bodies are considered impenetrable in the frame of Linear Theory of Elasticity, specific cells having vanishing pressure but negative gap *hi j* (,) (which equals the residual), are reinstated in the current contact area. To this end, the corresponding nodal pressures are adjusted with a quantity proven (Polonsky & Keer, 1999) to be positive at all times, i.e. *rij* (,) .

$$\begin{cases} A\_{\mathbb{C}} \leftarrow A\_{\mathbb{C}} - \left\{ (\mathbf{i}, \mathbf{j}) \right\}, p(\mathbf{i}, \mathbf{j}) = 0, (\mathbf{i}, \mathbf{j}) \in \left\{ (\mathbf{i}, \mathbf{j}) : p(\mathbf{i}, \mathbf{j}) < 0 \right\};\\ A\_{\mathbb{C}} \leftarrow A\_{\mathbb{C}} \cup \left\{ (\mathbf{i}, \mathbf{j}) \right\}, p(\mathbf{i}, \mathbf{j}) \leftarrow p(\mathbf{i}, \mathbf{j}) - ah(\mathbf{i}, \mathbf{j}), (\mathbf{i}, \mathbf{j}) \in \left\{ (\mathbf{i}, \mathbf{j}) : p(\mathbf{i}, \mathbf{j}) = 0 \land h(\mathbf{i}, \mathbf{j}) < 0 \right\}. \end{cases} \tag{46}$$

In the TC, cells for which Coulomb's friction law is not verified are removed from the stick region, and the corresponding shear tractions are set to the value of limiting friction. Also, cells having micro-slip **s**(,) *i j* not opposite to shear tractions **q**(,) *i j* pass from the slip to the stick region, yielding an adjusted stick area:

$$\begin{cases} A\_S \leftarrow A\_S - \{ (i, j) \}, \mathbf{q}(i, j) \leftarrow \mu p(i, j) \frac{\mathbf{q}(i, j)}{\| \mathbf{q}(i, j) \|}, (i, j) \in \left\{ (i, j) : \left\| \mathbf{q}(i, j) \right\| > \mu p(i, j) \right\}; \\ A\_S \leftarrow A\_S \cup \left\{ (i, j) : \mathbf{q}(i, j) \cdot \mathbf{s}(i, j) > 0 \right\}. \end{cases} \tag{47}$$

In either the NC or the TC, if any cell is removed from or re-enters the contact or the stick area, auxiliary variable is set to zero, otherwise it is set to unity. This variable allows resetting the conjugate directions *<sup>i</sup> d* once new cells enter or leave the system. Indeed, these

new cells have no precedent in the minimization process and therefore a new search for the optimal path in the CGM algorithm must be conducted.

11. Adjust contact tractions (i.e. system solution) according to the static equilibrium equations. The correction is sought in the same form as in step 1 of the algorithm:

$$\begin{aligned} \text{NC}: \quad &p(i,j) \leftarrow p(i,j) + a\_3 + a\_1 \mathbf{x}\_2(i,j) + a\_2 \mathbf{x}\_1(i,j), \quad (i,j) \in A\_P; \\ \text{TC}: \quad &\begin{bmatrix} q\_1(i,j) \\ q\_2(i,j) \end{bmatrix} \leftarrow \begin{bmatrix} q\_1(i,j) \\ q\_2(i,j) \end{bmatrix} + \begin{bmatrix} a\_3 \mathbf{x}\_2(i,j) + a\_1 \\ a\_3 \mathbf{x}\_1(i,j) + a\_2 \end{bmatrix}; \quad (i,j) \in A\_{S'} \end{aligned} \tag{48}$$

Numerical Simulation of Slip-Stick Elastic Contact 145

**q**(,) *i j* , i.e.

**Figure 3.** Flow-chart for the uncoupled NC or TC

( , , ) ( , , 1), ( , ) ( 1) *<sup>C</sup>* **q q** *ijk ijk ij A k* .

obtain pressure ( , , ),( , ) ( ) *<sup>C</sup> pi jk i j A k* .

than when the full load is applied in one step.

**5. Results and discussions** 

(, , ) (, , ) *old p i jk pi jk* .

step 4.

1. Acquire the state after 1 *k* loading increments and establish the new increment.

3. Solve the uncoupled NC using the algorithm described in the previous section, and

4. Memorize the obtained pressure for subsequent error computation:

5. Solve the uncoupled TC having as input the previously computed pressure, and obtain the shear tractions ( , , ),( , ) ( ) *<sup>C</sup>* **q** *ijk ij A k* . The latter constitutes a better approximation to

6. Solve the uncoupled NC again, using as input the shear tractions computed in step 5. The resulting pressure distribution *pi jk* (, , ) is in its turn a better approximation to

7. Verify if pressure distributions resulted in two subsequent computations, *pi jk* (, , ) and (, , ) *old p ijk* , vary within an imposed precision goal. If convergence is not met, return to

This algorithm can be used to simulate any loading history in the fully coupled three dimensional elastic contact with slip and stick. Essentially, three levels of iterations are employed. The inner level is responsible for solving either the normal (NC), or the tangential (TC) unconnected contact problem. The intermediate level mutually adjusts these solutions, based on the contribution of tractions in each direction to displacement fields. The outer level is employed to reproduce the loading history, and is related to friction being a dissipative process. Spinu and Frunza (Spinu & Frunza, 2011b) recently proved that results obtained by simulating the loading history are a closer match to existing analytical solutions

The newly advanced algorithm is validated in this section by comparison with existing closed-form solutions for the contact of similarly elastic materials. The results are then

2. Adopt vanishing guess value for the increment of shear tractions

coupled problem solution compared to shear tractions adopted in step 2.

coupled problem solution than the one obtained in step 3.

yielding the following correcting parameters:

$$\begin{aligned} \text{NC}: \begin{bmatrix} a\_3\\ a\_1\\ a\_2 \end{bmatrix} = \mathbf{M}^{-1} \cdot \begin{bmatrix} \text{W}/\text{\AA} - \sum\_{(i,j) \in A\_c} p(i,j) \\\ M\_1/\text{\AA} - \sum\_{(i,j) \in A\_c} p(i,j)x\_2(i,j) \\\ M\_2/\text{\AA} - \sum\_{(i,j) \in A\_c} p(i,j)x\_1(i,j) \end{bmatrix};\\ \text{TC}: \begin{bmatrix} a\_1\\ a\_2\\ a\_3 \end{bmatrix} = \mathbf{M}^{-1} \cdot \begin{bmatrix} \text{T}\_1/\text{\AA} - \sum\_{(i,j) \in A\_c} q\_1(i,j) \\\ \text{T}\_2/\text{\AA} - \sum\_{(i,j) \in A\_c} q\_2(i,j) \\\ \text{M}\_3/\text{\AA} - \sum\_{(i,j) \in A\_c} \left[ q\_2(i,j)x\_1(i,j) + q\_1(i,j)x\_2(i,j) \right] \end{bmatrix}. \end{aligned} \tag{49}$$

#### 12. Verify convergence to the imposed precision :

$$\begin{aligned} \text{NC}: \sum\_{\{i,j\} \neq A\_{\mathbb{C}}} \left| p(i,j) - p\_{\text{old}}(i,j) \right| \Big/ \sum\_{\{i,j\} \neq A\_{\mathbb{C}}} p\_{\text{old}}(i,j) \leq \varepsilon; \\ \text{TC}: \sum\_{\{i,j\} \neq A\_{\mathbb{C}}} \left\| \mathbf{q}(i,j) - \mathbf{q}\_{\text{old}}(i,j) \right\| \Big/ \sum\_{\{i,j\} \neq A\_{\mathbb{C}}} p(i,j) \leq \varepsilon. \end{aligned} \tag{50}$$

and return to step 2 if the precision goals are not met. The algorithm can be summarized in the flow-chart depicted in Fig. 3.

#### **4.2. Numerical solution of coupled contact problems**

In the general case, the NC and the TC cannot be solved independently, as the two problems are coupled, i.e. solution of each one requires resolution of the other. While the NC can be decoupled from the TC by neglecting the normal displacements induced by tangential tractions, as in case of Goodman approximation, solution of NC is always required to solve the TC, as shear tractions are linked to pressure in the slip area through Coulomb's law. The solution of the partial slip-stick elastic contact problem, involving the coupled NC and TC, can be obtained in an iterative manner, using the algorithm consisting in the following steps.

**Figure 3.** Flow-chart for the uncoupled NC or TC

optimal path in the CGM algorithm must be conducted.

yielding the following correcting parameters:

2

*a*

*a*

3

12. Verify convergence to the imposed precision

*a*

*a*

new cells have no precedent in the minimization process and therefore a new search for the

11. Adjust contact tractions (i.e. system solution) according to the static equilibrium equations. The correction is sought in the same form as in step 1 of the algorithm:

> 1 1 32 1 2 2 31 2

*q ij q ij ax ij a TC ij A q ij q ij ax ij a*

(,) (,) (,)

(,) 3 1

*NC a M p i j x i j*

11 2 (,)

: (,) (,) ;

(,) 1

*TC a T q ij*

2 2 2

(,)

(,) (,)

*ij A ij A old ij A ij A*

*TC i j i j p i j*

*NC p i j p i j p i j*

(,) (,)

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

: (,)

1

**Μ**

summarized in the flow-chart depicted in Fig. 3.

**4.2. Numerical solution of coupled contact problems** 

in the following steps.

**Μ**

3 12 21

*P*

(48)

(49)

*S*

2

**q q** (50)

, ) (,)

*jx ij*

.

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

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

2 1 (,)

*M pi jx i j*

*ij A*

*ij A*

*C*

*ij A*

*W pi j*

 

*C*

1 1

*ij A*

*T q ij*

*ij A*

(,)

(,)

*C*

(,) (,)

(,)

(,) (,) (

*C*

*C*

3 21 1

:

*old old*

: (,) (,) (,) ;

: (,) (,) (,) , *C C*

and return to step 2 if the precision goals are not met. The algorithm can be

*C C*

In the general case, the NC and the TC cannot be solved independently, as the two problems are coupled, i.e. solution of each one requires resolution of the other. While the NC can be decoupled from the TC by neglecting the normal displacements induced by tangential tractions, as in case of Goodman approximation, solution of NC is always required to solve the TC, as shear tractions are linked to pressure in the slip area through Coulomb's law. The solution of the partial slip-stick elastic contact problem, involving the coupled NC and TC, can be obtained in an iterative manner, using the algorithm consisting

*M q ijx ij q i*

*NC p i j p i j a a x i j a x i j i j A*


This algorithm can be used to simulate any loading history in the fully coupled three dimensional elastic contact with slip and stick. Essentially, three levels of iterations are employed. The inner level is responsible for solving either the normal (NC), or the tangential (TC) unconnected contact problem. The intermediate level mutually adjusts these solutions, based on the contribution of tractions in each direction to displacement fields. The outer level is employed to reproduce the loading history, and is related to friction being a dissipative process. Spinu and Frunza (Spinu & Frunza, 2011b) recently proved that results obtained by simulating the loading history are a closer match to existing analytical solutions than when the full load is applied in one step.
