Open access peer-reviewed chapter

# Numerical Simulation of Slip-Stick Elastic Contact

By Sergiu Spinu and Dumitru Amarandei

Submitted: November 29th 2011Reviewed: March 29th 2012Published: September 19th 2012

DOI: 10.5772/48451

## 1. Introduction

Fretting defines a condition in which mechanical contacts are subjected to alternating tangential displacements, small compared to dimensions of contact area, due to oscillating loading conditions. Fretting wear and fretting fatigue are between the most important factors responsible for contact failure, especially when high loads are transmitted through non-conforming contacts, leading to highly localized stress concentrators in the vicinity of the contact region. Prediction of life span of machine elements working in such conditions requires assessment of stress and strain in the contacting bodies, which is the main subject of Contact Mechanics. Although fretting is intrinsically a multidisciplinary process, involving adhesion, oxidation, abrasion and pitting, modern approach suggests that contact stresses play a chief role.

While analytic solutions in this research fieldlead to complex mathematical models, many without closed-form solution, numerical approach reveals itself as a useful engineering tool, capable of extending the few existing analytical results to technologically important contact scenarios. A numerical study may advance the understanding of fretting contact and provide assistance to the design of contacts with improved load-carrying capacity.

Elastic contact analysis considering interfacial friction and slip-stick behaviour originated in the works of Cattaneo (Cattaneo, 1938) and Mindlin(Mindlin, 1949). They proved independently that, even when the contacting bodies are globally sticking, a peripheral region of slip is to be assumed in order to remain in the frame of Linear Theory of Elasticity and to obey the Coulomb’s law of friction. Based on these results, Johnson (Johnson, 1985) advanced the closed-form solution for the contact between similarly elastic materials undergoing a fretting loop.

In case of dissimilarly elastic materials, when the effects of normal and tangential tractions are coupled, an iterative solution has been achieved (Hills et al., 1993) only for the plane (i.e. cylindrical) contact. Many authors employ the so-called Goodman approximation (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 threedimensional 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 semi-analytical (SAM) in a review paper by Renauf et al. (Renauf et al., 2011). The strongpoint 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 required.

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 elasticlayered half-spaces. However, the loading history was not accounted forin these studies, i.e. the full load was applied in one step.

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

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

In contact problem formulation, it is convenient to describe the initial geometries hi(i)(x1,x2)of the two contacting bodies i=1,2in a Cartezian coordinate systems having its x1and x2axes 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 x3-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 T(T1,T2), the bending (or flexing) moments M1,M2and the torsional moment M3. 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 ui(j), and move as rigid-bodies with translations ωi(j)and rotations ϕi(j), with i=1,2,3,j=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 p(j)in the normal direction and shear traction q(j)(q1,q2)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:

W(t)=ΓC(t)p(x1,x2,t)dx1dx2;E1
M1(t)=ΓC(t)p(x1,x2,t)x2dx1dx2;M2(t)=ΓC(t)p(x1,x2,t)x1dx1dx2.E2
1. The geometrical condition of deformation:

h(x1,x2,t)=hi(x1,x2)+u3(x1,x2,t)ω3(t)ϕ1(t)x2ϕ2(t)x1,(x1,x2)ΓC(t).E3
1. The contact complementarity conditions:

{p(x1,x2,t)>0h(x1,x2,t)=0,(x1,x2)ΓC(t);p(x1,x2,t)=0h(x1,x2,t)>0, (x1,x2)ΓC(t).E4

The temporal dimension tis included in this model along the spatial dimensions x1,x2to provide basis for reproduction of the loading history. Consequently, ΓC(t)denotes the contact area, and h(x1,x2,t)the surface separation, both established at a specified point ton the loading curve. When t=0, h(x1,x2,0)=hi(x1,x2)=hi(1)(x1,x2)+hi(2)(x1,x2). The other relative (composite) quantities, lacking the superscript indicating the contacting body, are defined in a similar way: ω3(t)=ω3(1)(t)+ω3(2)(t), and ϕi(t)=ϕi(1)(t)+ϕi(2)(t),i=1,2. Computation of the relative normal displacement u3(x1,x2,t)will be discussed in Section 3.2.

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 t=0, in undeformed state) profile of the contacting bodies, but in positions (relative to the initial point of contact O) corresponding to rigid-body translations ω3(1)(t)and ω3(2)(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 u3(1)(t)and u3(2)(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):

Ti(t)=ΓC(t)qi(x1,x2,t)dx1dx2,i=1,2;E5
M3(t)=ΓC(t)[q2(x1,x2,t)x1q1(x1,x2,t)x2]dx1dx2.E6
1. The geometrical condition of deformation in the timeframe [t1,t2], in which the tangential load is not allowed to change sign:

[s1(x1,x2,t2)s1(x1,x2,t1)s2(x1,x2,t2)s2(x1,x2,t1)]=[u1(x1,x2,t2)u1(x1,x2,t1)u2(x1,x2,t2)u2(x1,x2,t1)]...[ω1(t2)ω1(t1)ω2(t2)ω2(t1)](ϕ3(t2)ϕ3(t1))[x2x1],(x1,x2)ΓC(t2).E7
1. The contact complementarity conditions:

{q(x1,x2,t2)μp(x1,x2,t2)s(x1,x2,t2)s(x1,x2,t1)=0,(i,j)ΓS(t2);q(x1,x2,t2)=μp(x1,x2,t2)s(x1,x2,t2)s(x1,x2,t1)>0,(i,j)ΓC(t2)ΓS(t2).E8

Here, ΓSis the stick area, ΓCΓSthe slip region, μthe frictional coefficient and s(s1,s2)the relative slip distances. The base for Eq. (7) is presented in Fig. 2, which depicts the slip-stick contact process in the direction of x1. 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 [t1,t2]. Let us consider two points P1and P2on the contacting bodies (1) and (2), respectively, located at t1on 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 A1and A2, respectively, where the bodies deform due shear tractions.

Firstly, the above mentioned axis is assumed to pass through the initial point of contact, thereforeA1(t1)A2(t1)O. In the considered timeframe, all points in the contacting bodies undergo rigid-body translations ω1(1)and ω1(2)along the direction of x1, and points A1and A2also undergo tangential displacements. If for these points the composite parameters cancel each other, i.e. ω1(t2)=u1(t2), with ω1=ω1(1)+ω1(2)and u1=u1(2)u1(1); therefore,A1(t2)A2(t2)Oand consequently a stick regime is established in O.

Secondly, the P1P2axis intersects the contact area in the peripheral region, where, althoughA1(t1)A2(t1), A1(t2)diverge from A2(t2)with a relative (composite) slip distance s1=s1(1)+s1(2). This position corresponds to a region in relative slip (also referred to as micro-slip). 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-and- error 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

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 ΓPof the contacting bodies needs to be considered, which constitutes an important advantage of SAM over other numerical techniques. In ΓP, 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Δ2depends on the grid steps Δiin the direction of xi, i=1,2. The grid control points (centroids of rectangular elementary patches) are identified by a pair of indices (i,j), with 1i<N1, 1i<N2. Any continuous distributionf(x1,x2)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 f(i,j)=f(x1*,x2*), where x1*,x2*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 ΓP.

To simulate the loading history, additional discretization is performed in the temporal domain, meaning the load is imposed in small steps k=1,N¯. Consequently, p(i,j,k)denotes the nodal (elementary) pressure at the intersection of the line iwith the column jof the rectangular grid, achieved after application of kloading 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:

W(k)=Δ(i,j)AC(k)p(i,j,k);E9
M1(k)=Δ(i,j)AC(k)p(i,j,k)x2(i,j);M2(k)=Δ(i,j)AC(k)p(i,j,k)x1(i,j);E10
h(i,j,k)=hi(i,j)+u3(i,j,k)ω3(k)ϕ1(k)x2(i,j)ϕ2(k)x1(i,j),(i,j)AP;E11
{p(i,j,k)>0h(i,j,k)=0,(i,j)AC(k);p(i,j,k)=0h(i,j,k)>0,(i,j)APAC(k).E12

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

Tn(k)=Δ(i,j)AC(k)qn(i,j,k),n=1,2;E13
M3(k)=Δ(i,j)AC(k)[q2(i,j,k)x1(i,j)q1(i,j,k)x2(i,j)];E14
[s1(i,j,k)s1(i,j,k1)s2(i,j,k)s2(i,j,k1)]=[u1(i,j,k)u1(i,j,k1)u2(i,j,k)u2(i,j,k1)][ω1(k)ω1(k1)ω2(k)ω2(k1)]...                                  (ϕ3(k)ϕ3(k1))[x2(i,j)x1(i,j)],(i,j)AC(k);E15
{q(i,j,k)μp(i,j,k)s(i,j,k)s(i,j,k1)=0,(i,j)AS(k);q(i,j,k)=μp(i,j,k)s(i,j,k)s(i,j,k1)>0,(i,j)AC(k)AS(k).E16

In these equations, AP,ACand ASare the discrete counterparts of ΓP, ΓCand ΓS, 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 isthus 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:

uij(k)(x1,x2)=ΓCpj(k)(x1,x2)Gij(k)(x1x1,x2x2)dx1dx2;i,j=1,2,3;k=1,2,E17

where Gij(k)(x1,x2)denotes displacement of a point in body (k), in the direction of xi, induced by a unit point force applied in origin along direction of xj. The functions Gfor 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:

uζξ(n)(i,j)=(k,)ACKζξ(n)(ik,j)pξ(n)(k,),E18

where Kζξ(n)(ik,j)is the influence coefficient, expressing the displacement in the direction of xζ, induced in the cell (i,j)of the body (n)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)(ik,j)=x2()Δ2/2x2()+Δ2/2x1(k)Δ1/2x1(k)+Δ1/2Gζξ(n)(x1(i)x1,x2(j)x2)dx1dx2.E19

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

uζξ(n)=Kζξ(n)pξ(n)E20

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

[u1u2u3]=[u1(2)u1(1)u2(2)u2(1)u3(2)+u3(1)]=[K11(2)K12(2)K13(2)K21(2)K22(2)K23(2)K31(2)K32(2)K33(2)][q1(2)q2(2)p(2)][K11(1)K12(1)K13(1)K21(1)K22(1)K23(1)K31(1)K32(1)K33(1)][q1(1)q2(1)p(1)]==[K11(2)+K11(1)K12(2)+K12(1)K13(2)K13(1)K21(2)+K21(1)K22(2)+K22(1)K23(2)K23(1)K31(2)K31(1)K32(2)K32(1)K33(2)+K33(1)][q1q2p].E21

When the second subscript is omitted, the referred displacement is the sum of all contributions. The positive sign in computation of u3is related to the fact that normal displacements are computed in coordinate systems linked to each body, having the x3axes pointing inward. It should also be noted that from pressure definition p(2)=p(1)=p, but qi(2)=qi(1)=qi,i=1,2, as the mutual actions of two bodies upon each other are equal, but directed to contrary parts. The influence coefficients Km(n)can be computed using the following relations:

Km(n)(ν(n),E(n),i,j)=km(n)(ν(n),E(n),x1(i,j)+Δ12,x2(i,j)+Δ22)+km(n)(ν(n),E(n),x1(i,j)Δ12,x2(i,j)Δ22)...km(n)(ν(n),E(n),x1(i,j)+Δ12,x2(i,j)Δ22)km(n)(ν(n),E(n),x1(i,j)Δ12,x2(i,j)+Δ22);E22
k33(n)(ν(n),E(n),x1,x2)=1(ν(n))2πE(n)[x1ln(x2+x12+x22)+x2ln(x1+x12+x22)];E23
)
k13(n)(ν(n),E(n),x1,x2)=(1+ν(n))(12ν(n))2πE(n)[2x1tan1(x12+x22x2x1)x2ln(x12+x22)];E24
k11(n)(ν(n),E(n),x1,x2)=1+ν(n)πE(n)[x2ln(x1+x12+x22)x2+(1ν(n))x1ln(x2+x12+x22)];E25
k21(n)(ν(n),E(n),x1,x2)=ν(n)(1+ν(n))πE(n)x12+x22;E26
k23(n)(ν(n),E(n),x1,x2)=k13(n)(ν(n),E(n),x2,x1);E27
k31(n)(ν(n),E(n),x1,x2)=k13(n)(ν(n),E(n),x1,x2);E28
k12(n)(ν(n),E(n),x1,x2)=k21(n)(ν(n),E(n),x1,x2);E29
k22(n)(ν(n),E(n),x1,x2)=k11(n)(ν(n),E(n),x2,x1);E30
k32(n)(ν(n),E(n),x1,x2)=k23(n)(ν(n),E(n),x2,x1),E31

where E(n)and ν(n)are the elastic constants (Young modulus and Poisson’s ratio, respectively) of materials of the two contacting bodies. If E(1)=E(2)=Eand ν(1)=ν(2)=ν, then Kij(1)=Kij(2)=Kij,i,j=1,2,3,and Eq. (21) takes a simplified form:

[u1u2u3]=2[K11K120K21K22000K33][q1q2p].E32

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 Nsamples, requires O(N2)operations in time/space domain, but only O(NlogN)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 ACand ASrespectively, is not known in advance. Moreover, resolution of a linear system with Nunknowns has an order of computation of O(N3)when direct methods like Gaussian elimination are used. In contact problems, grids with N=103nodes 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 N=106points 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.

### 4.1. Numerical solution of uncoupled contact problems

Eq. (21) suggests that, in case of dissimilarly elastic materials (i.e. E(1)E(2)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 kis 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. δu1(i,j)=u1(i,j,k)u1(i,j,k1).

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:

NC:θ=0;Rold=1;d3(i,j)=0,pold(i,j)=0,(i,j)AP;TC:θ=0;Rold=1;d1(i,j)=0,d2(i,j)=0,qold(i,j)=0,(i,j)AP.E33

Here, di,i=1,2,3is the initial descend direction in the CGM. The role of the other variables will be discussed later.

1. 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=AP, and in the TC that all cells in the contact area are in stick, i.e. AS=AC. The contact tractions are sought in the following form:

NC:p(i,j)=a3+a1x2(i,j)+a2x1(i,j),(i,j)AP;TC:{q1(i,j)=a3x2(i,j)+a1;q2(i,j)=a3x1(i,j)+a2;(i,j)AC,E34

where the unknown coefficients ai(k),i=1,2,3are computed by plugging Eq. (34) into Eqs. (9), (10) and (13), (14), respectively, yielding:

NC:[a3a1a2]T=Μ1[WM1M2]T/Δ;TC:[a1a2a3]T=Μ1[T1T2M3]T/Δ;E35
NC:Μ=[(i,j)AC1(i,j)ACx2(i,j)(i,j)ACx1(i,j)(i,j)ACx2(i,j)(i,j)ACx22(i,j)(i,j)ACx2(i,j)x1(i,j)(i,j)ACx1(i,j)(i,j)ACx1(i,j)x2(i,j)(i,j)ACx12(i,j)];TC:M=[(i,j)AS1(i,j)AS0(i,j)ASx2(i,j)(i,j)AS0(i,j)AS1(i,j)ASx1(i,j)(i,j)ASx2(i,j)(i,j)ASx1(i,j)(i,j)AS[x12(i,j)+x22(i,j)]].E36
1. 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, δu1(i,j),δu2(i,j), is needed.

2. 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:

NC:hi(i,j)+u3(i,j)ω3ϕ1x2(i,j)ϕ2x1(i,j)=0,(i,j)AC;TC:[δu1(i,j)δu2(i,j)][δω1δω2]+δϕ3[x2(i,j)x1(i,j)]=,(i,j)AS,E37

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

NC:(ω3,ϕ1,ϕ2)=(i,j)AC(hi(i,j)+u3(i,j)ω3ϕ1x2(i,j)ϕ2x1(i,j))2;TC:(δω1,δω2,δϕ3)=(i,j)AS((δu1(i,j)δω1δϕ3x2(i,j))2+(δu2(i,j)δω2δϕ3x1(i,j))2).E38

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:

NC:[ω3ϕ1ϕ2]=Μ1[(i,j)AC[u3(i,j)+hi(i,j)](i,j)AC[(u3(i,j)+hi(i,j))x2(i,j)](i,j)AC[(u3(i,j)+hi(i,j))x1(i,j)]];TC:[δω1δω2δϕ3]=Μ1[(i,j)ASδu1(i,j)(i,j)ASδu2(i,j)(i,j)AS[δu1(i,j)x2(i,j)+δu2(i,j)x1(i,j)]].E39

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 systemsare solved on successive contact and stick areas, until static equilibrium equations and complementarity conditions are fully satisfied for every cell in the grid.

1. Start the conjugate gradient loop by computing the residuals r(i,j),=1,2,3:

NC:r3(i,j)=hi(i,j)+u3(i,j)ω3ϕ1x2(i,j)ϕ2x1(i,j),(i,j)AP;TC:[r1(i,j)r2(i,j)]=[δu1(i,j)δu2(i,j)][δω1δω2]+δϕ3[x2(i,j)x1(i,j)],(i,j)AS,E40

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

NC:R=(i,j)AC(r3(i,j))2;TC:R=(i,j)AS[(r1(i,j))2+(r2(i,j))2].E41
1. Compute the descend direction di,i=1,2,3in 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:

NC:d3(i,j){r3(i,j)(Rθ/Rold)d3(i,j),(i,j)AC;0,(i,j)APAC;TC:[d1(i,j)d2(i,j)]{[r1(i,j)r2(i,j)]+RθRold[d1(i,j)d2(i,j)],(i,j)AS;T,(i,j)ACAS,E42

and overwrite the content of Rold: RoldR.

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

NC:α=(i,j)ACr3(i,j)d3(i,j)/(i,j)ACc3(i,j)d3(i,j);TC:α=(i,j)ASr1(i,j)d1(i,j)+(i,j)ASr2(i,j)d2(i,j)(i,j)ASc1(i,j)d1(i,j)+(i,j)ASc2(i,j)d2(i,j),E43

whereci,i=1,2,3denotes the convolution of the influence coefficients matrix with the descend direction:

NC:c3=[K33(2)+K33(1)]d3;TC:[c1c2]=[K11(2)+K11(1)K12(2)+K12(1)K21(2)+K21(1)K22(2)+K22(1)][d1d2].E44
1. Memorize the current contact tractions in a new variable for subsequent error computation: poldp,qoldq.

2. Update the system solution by making a step of length α, computed using Eq.(43), along direction di,i=1,2,3, derived in Eq. (42):

NC:p(i,j)p(i,j)αd3(i,j),(i,j)AC;TC:[q1(i,j)q2(i,j)][q1(i,j)q2(i,j)]α[d1(i,j)d2(i,j)],(i,j)AS.E45
1. 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 pressuresare 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 h(i,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. αr(i,j).

{ACAC{(i,j)},p(i,j)=0,(i,j){(i,j):p(i,j)<0};ACAC{(i,j)},p(i,j)p(i,j)αh(i,j),(i,j){(i,j):p(i,j)=0h(i,j)<0}.E46

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:

{ASAS{(i,j)},q(i,j)μp(i,j)q(i,j)q(i,j),(i,j){(i,j):q(i,j)>μp(i,j)};ASAS{(i,j):q(i,j)s(i,j)>0}.E47

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

1. 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:

NC:p(i,j)p(i,j)+a3+a1x2(i,j)+a2x1(i,j),(i,j)AP;TC:[q1(i,j)q2(i,j)][q1(i,j)q2(i,j)]+[a3x2(i,j)+a1a3x1(i,j)+a2];(i,j)AS,E48

yielding the following correcting parameters:

NC:[a3a1a2]=Μ1[W/Δ(i,j)ACp(i,j)M1/Δ(i,j)ACp(i,j)x2(i,j)M2/Δ(i,j)ACp(i,j)x1(i,j)];TC:[a1a2a3]=Μ1[T1/Δ(i,j)ACq1(i,j)T2/Δ(i,j)ACq2(i,j)M3/Δ(i,j)AC[q2(i,j)x1(i,j)+q1(i,j)x2(i,j)]].E49
1. Verify convergence to the imposed precision ε:

NC:(i,j)AC|p(i,j)pold(i,j)|/(i,j)ACpold(i,j)ε;TC:(i,j)ACq(i,j)qold(i,j)/(i,j)ACp(i,j)ε,E50

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
1. Acquire the state after k1loading increments and establish the new increment.

2. Adopt vanishing guess value for the incrementof shear tractions δq(i,j), i.e. q(i,j,k)=q(i,j,k1),(i,j)AC(k1).

3. Solve the uncoupled NC using the algorithm described in the previous section, and obtain pressure p(i,j,k),(i,j)AC(k).

4. Memorize the obtained pressure for subsequent error computation: pold(i,j,k)p(i,j,k).

5. Solve the uncoupled TC having as input the previously computed pressure, and obtain the shear tractions q(i,j,k),(i,j)AC(k). The latter constitutes a better approximation to coupled problem solutioncompared to shear tractions adopted in step 2.

6. Solve the uncoupled NC again, using as input the shear tractions computed in step 5. The resulting pressure distribution p(i,j,k)is in its turn a better approximation to coupled problem solution than the one obtained in step 3.

7. Verify if pressure distributions resulted in two subsequent computations, p(i,j,k)and pold(i,j,k), vary within an imposed precision goal. If convergence is not met, return to step 4.

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.

## 5. Results and discussions

The newly advanced algorithmis validated in this section by comparison with existing closed-form solutions for the contact of similarly elastic materials. The results are then extended by simulating the three-dimensional contact between dissimilarly elastic materials, proving the capacity of the program to solve a large variety of problems incorporating the slip-stick processes in elastic mechanical contacts.

### 5.1. Simulation of a fretting loop

The newly proposed algorithm is first validated by comparison with the closed-form solution advanced in (Johnson, 1985) for the spherical contact undergoing a fretting loop. A steel ball of radius R=18mmis pressed against an elastic half-space having the same elastic properties, with a normal force W=1kN. In this case, the NC is uncoupled from the TC, as shown by Eq. (32). A tangential force T1oscillating between two limiting values Tlimand Tlim, where Tlim=0.9μW, is subsequently applied. During a fretting contact process, it is expected that friction vary on the contact area, as well as with accumulation of debris particles resulted from additional wear. However, for validation purposes, a frictional coefficient uniform over all contact area and constant during load application is assumed in this study. The numerical approach can equally handle mapped distributions of μ, if such information is available. The loading history for fretting processes is depicted in Fig. 4, in which the load is applied incrementally in N=700equal steps. Figure 4.The loading history in simulation of a fretting loop; all moments are assumed to vanish

Nis chosen differently in the following simulations. When studying the contact between similarly elastic materials, accurate (i.e. concurring with the existing closed-form solution) results can be obtained with even a small number of increments (N=42increments in simulations depicted in Figs. 5 and 6). Moreover, it is found that every different state on the loading path can be obtained by simulating only the states when the tangential load changes its sign of, e.g. every state M on the trajectory DF require computation of three states: the ones corresponding to points B and D, as well as the current one (corresponding to point M). This is not the case when simulating the contact between dissimilarly elastic materials. It is found that a large number of increments is required to obtain the detailed contact behaviour, due to non-overlapping stick regions from one loading step to the next. This feature was also observed by Gallego et al. (Gallego et al., 2010), who pointed out that a waved shear tractions profile is predicted numerically when the number of increments is small. Although some perturbations still appear, the large number of increments used in the current simulations N=700led to well converged numerical solutions.

The Hertz frictionless theory for the similarly elastic contact scenario predicts a central pressure pH=1.996GPaand a contact radius aH=0.489mm, which are used as normalizers. Shear tractions profiles corresponding to equal loading levels but laying on different trajectories are depicted in Fig. 5, proving that the current state depends not only on the present load level, but also on its past evolution, thus showing a hysteretic behaviour or memory effect. Only profiles corresponding to states on AD are shown in Fig. 5. States on the DF trajectory can be obtained from those on AD by symmetry with respect to the x1-axis, and states on FH overlap those on BD corresponding to the same loading level. Figure 5.Profiles of shear tractions in the plane x2=0 corresponding to different points on the loading curve, μ=0.1

As depicted in Fig. 6, the force-displacement curve forms a hysteretic loop, also referred to as a fretting loop. The rigid-body tangential displacement is normalized by its value ωlimcorresponding to Tlim. The analytical model predicts that the states corresponding to points B and F on the loading curve should overlap, which is also obtained through numerical simulation, within an imposed precision.

From a mechanical point of view, regions undergoing the most severe stresses are of interest, leading to prediction of yield inception or crack nucleation in the contacting bodies. An algorithm to assess stresses due to known, but otherwise arbitrarily distributed contact tractions is readily available (Liu & Wang, 2002). A typical distribution of von Mises equivalent stress σVM/pHin the plane x2=0, corresponding to point B on the loading path, is depicted in Fig. 7. Evolution of magnitude and of depth of von Mises maximum stress in a fretting loop is depicted in Figs. 8-9, for various frictional coefficients.

It is found through numerical simulation that shear tractions have a weak impact on the maximum von Mises stress when μ<0.2. However, larger friction processes lead to a competition between the in-depth maximum, which fluctuates around its position in the frictionless contact (corresponding to point A in Fig. 9), and a second extremum developing on the surface, at the trailing edge of the contact. The latter can seriously diminish the load-carrying capacity of the contact, especially in case of surfaces with poorly controlled surface quality, due to superimposition of microtopography-induced stress perturbations.

### 5.2. Simulation of torsional contact

Program validation is subsequently performed in case of a spherical contact undergoing torsion applied simultaneously to a normal constant force. Based on the results of this author, (Mindlin, 1949), Johnson (Johnson, 1985) presents the closed-form solution for this contact scenario when a partial slip regime is established on the contact area (i.e. when the torsional moment M3is smaller than a limiting value M3lim=3πμWaH/16). The solution is later reviewed and enhanced for the case of viscoelastic materials by Dintwa et al. (Dintwa et al., 2005).

A spherical indenter of radius R=18mmis pressed with a normal force W=1kNagainst an elastic half-space, having the same elastic parameters, and subsequently an increasing torsional moment M3M3limis applied, leading to shear tractions q1and q2depicted in Figs. 10a and b, respectively. The norm of these tractionsq(i,j)=q1(i,j)+q1(i,j)is compared in Fig. 11 with analytical results,using cylindrical coordinates.

Figure 12 depicts the dimensionless stick radius when the torsional moment M3is varied in the domain corresponding to partial slip. In all investigated cases, a good agreement between analytical results and numerical predictions is found, giving further confidence in the newly advanced computer program. Figure 10.Distribution of shear tractions in torsional contact, μ=0.1, M3=0.9M3 lim Figure 11.Profiles of radial shear tractions, μ=0.1

### 5.3. Extension of results

The numerical program advanced in this study is subsequently used to simulate the fretting contact between dissimilarly elastic materials. To our best knowledge, an analytical solution to this contact process has not been advanced yet. The ball in the previous example is considered rigid, and the loading history is simulated using N=700equal increments. Figure 13.Shear tractions profiles in the fretting contact between dissimilarly elastic materials, corresponding to different points on the loading path, μ=0.3 Figure 14.Shear tractions profiles in the fretting contact between dissimilarly elastic materials, corresponding to different points on the loading path, μ=0.6

Numerical simulations predict that the investigated fretting contact exhibits a unique behaviour in the first two loading cycles, after which it stabilizes to a fixed path. Figures 13 and 14 suggest that the shear tractions profiles corresponding to the BD and FH trajectories no longer match, as in case of similarly elastic materials. However, states D and H overlap, beside a few perturbations induced by the former boundaries of the stick area. Presumably, these perturbations are related to the discretization error, and the D and H states can be considered as concurring. Subsequent oscillating loading is expected to lead to states following the same fixed path, as in case of similarly elastic materials. Ananalogous behaviour was observed by Wang et al. (Wang et al., 2010) when studying numerically the partial slip contact of elastic layered half-spaces.

## 6. Conclusions

The work reported herein advances a numerical model for the fretting contact between dissimilarly elastic materials. A numerical approach is required to simulate this type of contact process, as analytical models can incorporate neither the loading history, which must be reproduced when friction is accounted for, nor the coupling of normal and tangential effects.

The implemented algorithm is based on three levels of iterations, fully incorporating the interconnectivity between normal and tangential tractions. The innermost level solves, in a conjugate gradient type loop, either the normal, or the tangential uncoupled contact problem, while the intermediate loop iterates between these solutions, until pressure convergence is reached. The outer level reproduces the loading history, and is based on the assumption that irreversibility of friction requiressimulation of all previous states in an incremental load process.

The strong points of this algorithm consist in reduced computational complexity, compared to finite element analysis, as well as in its ability to handle arbitrarily shaped contact geometries or imposed frictional coefficient maps. Comparison with existing closed-form solutions for the spherical contact undergoing a fretting loop or torsion, when a uniform frictional coefficient is assumed, gives confidence in the newly advanced model.

Evolution of stress state during a fretting loop in the spherical contact between similarly elastic materials is assessed. It is found that for specific combinations of loading levels and frictional coefficients, the most severely stressed region can be found on the bounding surfaces, at the trailing edge of the contact.

Numerical simulations suggest that the fretting contact between dissimilarly elastic materials exhibits a unique path in the first two loading cycles, which stabilizes with subsequent oscillating loading to a fixed trajectory, as in case of similarly elastic materials.

The newly advanced algorithm is expected to solve a large variety of contact problems involving interfacial friction, leading to a better understanding of complex multidisciplinary phenomena like fretting wear and fretting fatigue, in which transient contact tractions and induced subsurface stresses play an important role, as well as to a more accurate prediction of contact failure through yield inception or crack nucleation. Study of partial slip elastic-plastic contact is anticipated for future contributions, by addition of the residual term, related to development of plastic strains.

## Acknowledgement

This paper was supported by the project “Progress and development through post-doctoral research and innovation in engineering and applied sciences – PRiDE - Contract no. POSDRU/89/1.5/S/57083”, project co-funded from European Social Fund through Sectorial Operational Program Human Resources 2007-2013.

## How to cite and reference

### Cite this chapter Copy to clipboard

Sergiu Spinu and Dumitru Amarandei (September 19th 2012). Numerical Simulation of Slip-Stick Elastic Contact, Numerical Simulation - From Theory to Industry, Mykhaylo Andriychuk, IntechOpen, DOI: 10.5772/48451. Available from:

### Embed this chapter on your site Copy to clipboard

<iframe src="http://www.intechopen.com/embed/numerical-simulation-from-theory-to-industry/numerical-simulation-of-slip-stick-elastic-contact" />

Embed this code snippet in the HTML of your website to show this chapter

### More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

### Related Content

Next chapter

#### Numerical Simulations of Radiation and Scattering Characteristics of Dipole and LOOP Antennas

By Oleg A. Yurtsev and Grigory V. Ptashinsky

#### Two Phase Flow, Phase Change and Numerical Modeling

Edited by Amimul Ahsan

First chapter

#### Modeling the Physical Phenomena Involved by Laser Beam – Substance Interaction

By Marian Pearsica, Stefan Nedelcu, Cristian-George Constantinescu, Constantin Strimbu, Marius Benta and Catalin Mihai

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

View all books