Mechanics of Materials 40 (2008) 936–948

Contents lists available at ScienceDirect

Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat

A numerical model for the point contact of dissimilar materials considering tangential tractions W. Wayne Chen *, Q. Jane Wang Department of Mechanical Engineering, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA

a r t i c l e

i n f o

Article history: Received 31 May 2007 Received in revised form 14 June 2008

Keywords: Contact mechanics Stress analysis Dissimilar-material contact

a b s t r a c t This paper presents a three-dimensional numerical model for the simulation of the contacts of elastically dissimilar materials. This model is based on the Boussinesq–Cerruti integral equations relating normal and tangential surface tractions to surface displacements and employs the static friction law to identify the onset of local micro-slip. An iterative conjugate gradient method (CGM) is utilized to determine the unknown contact and stick area efficiently and accurately with the assistance of the discrete convolution and fast Fourier transform (DC-FFT) algorithm. The model leads to the solutions of surface real contact and stick areas, pressure, tangential tractions, and the subsurface stress field. It is verified through the comparison of the numerical results with analytical solutions. Simulations are performed for a ball and half-space contact with the normal loading alone first, and then the results for the contact with increasing tangential loading are also presented. The effects of shear tractions on the contact area, the stick zone, pressure and the intensity of vonMises stress are discussed. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction The classical analysis of frictionless contacts between two counterformal elastic bodies with quadratic surface geometry presented by Hertz facilitates investigations of the interaction between mechanical components (Johnson, 1985). However, the frictionless assumption is only valid when contact bodies have the same material properties or can not sustain shear tractions. Shear tractions will take place at the surface when bodies with dissimilar materials are brought into contact. If the shear traction magnitude in a contact area exceeds the static friction limit, the contact area may experience micro-slip opposite to the friction direction; otherwise it remains stick. When an increasing external tangential loading is applied, more and more contact areas are involved in the relative slip until a gross sliding begins. It is important to study the stick–slip phenomenon because it is seen in many mechanical systems. * Corresponding author. Tel.: +1 847 467 6961; fax: +1 847 491 3915. E-mail address: [email protected] (W.W. Chen). 0167-6636/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechmat.2008.06.002

The exact partial slip solution of cylindrical contacts of similar materials subjected to a global tangential loading was derived by Cattaneo (1938) and Mindlin (1949) independently. The Cattaneo–Mindlin solution has been proven to be valid for a general elliptical or circular contact (Johnson, 1985). Nowell et al. (1988) derived closedform solutions for interfacial tractions induced by the normal contact between dissimilar elastic cylinders when a tangential force, less than that causing a gross sliding, was applied. The coupling effect of shear tractions with normal pressure was included in the formulation. Spence (1973, 1975) gave analytical solutions of interfacial tractions due to the indentation of an elastic half-space by an axisymmetric punch. Spence first solved the problem with a flat-ended punch, and then used a transformation technique to obtain the solution for a spherical indenter. Numerical simulation is an effective approach to revealing the contact situation in the normal indentation with friction. Kosior and his co-workers (1999, 2000) modeled the normal elastic contact with Coulomb friction by means of the finite element method (FEM) coupled with a domain

W.W. Chen, Q.J. Wang / Mechanics of Materials 40 (2008) 936–948

937

Nomenclature Letters a0 Hertz contact radius, mm Ac, Ast,Asl element sets of contact zone, stick zone, and slip zone c stick zone radius D, B the influence coefficient E1, E2 Young’s moduli of two contact bodies, GPa 0 E , E* equivalent Young’s moduli, GPa g square summation of slip distance in the stick zone G, T the Green’s function h, hi surface gap, initial surface clearance, mm p, qx, qy pressure, shear tractions, MPa p0 Hertz contact peak pressure, MPa sx, sy the slip distances along x and y axes, mm tx, ty conjugate direction ux, uy, uz surface displacements, mm W, Fx, Fy normal load, tangential load components, N x, y, z space coordinates, mm

decomposition method and the boundary element method (BEM). Brizmer et al. (2006, 2007) used FEM to investigate the normally loaded contact of an elastic–plastic sphere and a rigid flat with perfect slip and full stick conditions under combined normal and tangential loading. Kogut and Etsion (2003) investigated the real contact between a spherical and a half-space with FEM and treated the sliding inception as the moment when the entire contact zone yields. FEM is applicable for simulations with complicated material behaviors and contact shapes. However, the fixed boundary condition used in FEM needs a large simulation domain, which may increase the computational burden; and the surface elements in the potential contact region have to be treated in a special way. Semi-analytical methods (SAM) (Bjorklund and Andersson, 1994; Polonsky and Keer, 1999; Jacq et al., 2002; Li and Berger, 2003; Liu et al., 2007; Nelias et al., 2007) were developed to investigate the contact problems. The present model also utilizes SAM to solve the stick–slip contact problems under combined normal and tangential loadings. SAM is more efficient than FEM because only the contact area of interest needs to be meshed and closed-form displacement solutions of surface tractions are used. Due to the coupling effect between normal pressure and shear tractions, an iterative approach is used to determine the stick and slip status of the contact zone. It is noteworthy that contact pressure is assumed to be parallel everywhere, thus the present model may have difficulty simulating the contact of bodies with large geometry slope, such as a sharp conic tip. The contact problems studied in this paper (a sphere on a half-space) comply the above assumption and can be simulated by the present model without any difficulty. The investigation of the subsurface stress field is a necessary step leading to understanding of plasticity and failure. Hamilton (1983) developed explicit formulae

Greek letters b Dunders constant dx, dy, dz the rigid body deformation along three axes, mm D1, D2, D3 element sizes along three axes, mm lf0 friction coefficient l , l* equivalent shear moduli, GPa m1, m2 poisson ratio of two contact bodies rmn Cauchy stress components (m, n = x, y, z), MPa s updating length along conjugate direction Special marks * continuous convolution () single underline, the vector containing variable values at all elements () double underlines, the influence coefficient matrix

for stresses beneath a sliding Hertz contact, where shear traction was treated as a product of pressure and a specified friction coefficient. Using the solutions of Hertz pressure and Spence’s shear tractions as boundary conditions, Hills and Sackfield (1987) derived the equations for the stress field induced by the normal contact between two spheres having different elastic constants. Numerical models (Munisamy and Hills, 1992; Liu and Wang, 2002) have been developed to analyze the subsurface stress field caused by arbitrarily distributed pressure based on the generic stress solutions of uniform tractions in a surface rectangle (Hills et al., 1993). These models provide powerful approaches to evaluating the subsurface stress field in a half-space if the surface contact pressure or shear traction distributions are known. However, the surface tractions caused by the spherical contact subjected to both normal and tangential loads under the stick–slip condition are not well studied, and no analytical solution has been derived so far. This study uses SAM to investigate the stress field in a half-space under a stick–slip contact boundary. The purpose of this paper is to develop a three-dimensional numerical model for contacts of dissimilar materials. The model is applied to analyze a point stick–slip contact subjected to both normal and tangential loads. The potential contact area is discretized into small rectangular elements. General solutions of displacements due to uniform tractions over the element at the origin are derived. The interaction between pressure and shear tractions due to dissimilar elastic properties is also included. The CGM algorithm applied by Polonsky and Keer (1999) in normal contact simulations is extended here to determine the real stick zone. The discrete convolution and fast Fourier transform algorithm (DC-FFT) (Liu et al., 2000) is also used to improve the computation efficiency. The resulting surface tractions and subsurface stress field are presented.

938

W.W. Chen, Q.J. Wang / Mechanics of Materials 40 (2008) 936–948

where

2. Theoretical background 2.1. Problem formulations The geometry of the contact problem between a ball (body 2, E2 and m2) and a half-space (body 1, E1 and m1) is shown in Fig. 1. The x and y axes are on the surface, while the z axis points inwards the half-space. The ball indenter is compressed onto the half-space by a normal load, W, along the z axis. Tangential loads, Fx and Fy, are applied on the ball parallel to the x–y plane. The contact interaction results in normal pressure p and shear tractions qx and qy in the interface. Note that the tractions acting at a point on one surface should act on the corresponding point on the other surface with the same size but at the opposite direction. Based on the Boussinesq–Cerruti integral equations (Johnson, 1985), the relative surface displacements at any point (x, y) can be expressed as,

ux ðx; yÞ ¼

Z

þ1

1

Z

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 þ y 2 ;

r¼

þ1

fGxx ðx  x0 ; y  y0 Þ  qx ðx0 ; y0 Þ

1

þ Gxy ðx  x0 ; y  y0 Þ  qy ðx0 ; y0 Þ þ Gxz ðx  x0 ; y  y0 Þ  pðx0 ; y0 Þg dx0 dy0

ð1Þ

1=E ¼ ð1  m21 Þ=E1 þ ð1  m22 Þ=E2 ; 1=l0 ¼ ½ð1 þ m1 Þð1  2m1 Þ=2E1  ½ð1 þ m2 Þð1  2m2 Þ=2E2 ;

1=l ¼ ½m1 ð1 þ m1 Þ=E1 þ ½m2 ð1 þ m2 Þ=E2

  1 x2 y 2 xy ; Gxy ¼ þ ; pr3 E0 E pl r3 x xy Gxz ¼  0 2 ; Gyx ¼ ; pl r pl r3  2  1 x y2 y Gyy ¼ 3  þ 0 ; Gyz ¼  0 2 ; pr E E pl r x y 1 ; Gzy ¼ ; Gzz ¼  Gzx ¼ pl0 r2 pl0 r2 pE r Gxx ¼

W

NX 1 1 N 2 1  X n¼0

Here, symbol ‘*’ means the continuous convolution, Gmn (m, n = x, y, and z) are the Green’s functions, which are listed as follows,

my mz qx½n;w Dmx ½an;bw þ qy½n;w D½an;bw þ p½n;w D½an;bw

ð4Þ

where Dmn ½a;b are the discrete influence coefficients (ICs), which is the displacement um at (2aD1, 2bD2) due to a unit surface traction along the ‘n’ direction at the element centered at the origin.

Z

D1

D1

Z

D2

Gmn ð2aD1  x0 ; 2bD2  y0 Þ dx0 dy0

ð5Þ

D2

The closed-form expressions of indefinite integrals of these Green’s functions have been derived by James and Busby (1995). One has,

E 2 , v2 δx

Body 2

Fx

δz

x

E1 , v1



w¼0

0 6 a; n 6 N 1  1; 0 6 b; w 6 N 2  1; m; n ¼ x; y; z

Dmn ½a;b ¼

ð2Þ

ð3Þ

In order to perform numerical calculation, the contact surface of interest is meshed into small rectangular elements, which have the uniform size of 2D1  2D2. N1, N2 are the numbers of elements along the x and y directions. Tractions and displacements in each discrete patch are treated as constant, where the values at the center of a patch are used. qn[a,b] (n = x, y, z) and um[a,b] (m = x, y, z) are surface traction qn and elastic displacement um of the element centered at (2aD1, 2bD2), respectively. Here, qz = p. The displacement can then be written in the form of discrete convolution,

um½a;b ¼

uz ðx; yÞ ¼ Gzx  qx þ Gzy  qy þ Gzz  p

and



¼ Gxx  qx þ Gxy  qy þ Gxz  p uy ðx; yÞ ¼ Gyx  qx þ Gyy  qy þ Gyz  p

1=E0 ¼ ð1 þ m1 Þ=E1 þ ð1 þ m2 Þ=E2 ;

hi

Fy

uz

δx

δz

x

h

Body 1

y

z z

Fig. 1. Geometry of the contact problem under investigation, W is the normal load, and Fx and Fy are the tangential loads.

ux

sx

Fig. 2. Contact variables in the x–z plane, dx and dz are the rigid body approaches, ux and uz the normal elastic displacement, and h is the surface gap.

939

W.W. Chen, Q.J. Wang / Mechanics of Materials 40 (2008) 936–948

Dmn ½a;b ¼ fmn ð2aD1 þ D1 ; 2bD2 þ D2 Þ þ fmn ð2aD1 a

ux  dx ¼ sx ;

 fmn ð2aD1  D1 ; 2bD2 þ D2 Þ 0 6 a 6 N1  1; 0 6 b 6 N2  1; m; n ¼ x; y; z

ð6Þ

Gmn.

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 þ y2 Þ  x y lnðx þ x2 þ y2 Þ þ ; pE pE0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 þ y 2 ; f xy ðx; yÞ ¼  

fxx ðx; yÞ ¼

x lnðy þ

pl

f yx ðx; yÞ ¼ fxy ðy; xÞ;

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi y x2 þ y2  y þ x tan1 ; x f yy ðx; yÞ ¼ fxx ðy; xÞ;

f yz ðx; yÞ ¼ fxz ðy; xÞ;

f zx ðx; yÞ ¼ fxz ðx; yÞ;

f xz ðx; yÞ ¼ 

1 

pl0

y ln

Ac ¼ f½a; bjh½a;b ¼ 0; p½a;b > 0g Ac ¼ f½a; bjh½a;b > 0; p½a;b ¼ 0g

f zy ðx; yÞ ¼ fxz ðy; xÞ; and pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x lnðy þ x2 þ y2 Þ þ y lnðx þ x2 þ y2 Þ f zz ðx; yÞ ¼ pE

ð7Þ

P 1 1 PN2 1 mn We use um ¼ Dmn  qn ¼ Nn¼0 w¼0 D½an;bw qn½n;w to represent the basic discrete convolutions. Symbols with single underlines are vectors containing N1  N2 values of corresponding variables at all elements, and those with double underlines are the matrices of influence coefficients. Thus, Eq. (4) can be rewritten as,

2

3

2

32

3

Dxx

Dxy

Dxz

6 7 6 Dyx 4 uy 5 ¼ 4 Dzx uz

Dyy

6 7 Dyz 7 54 qy 5 zz p D

un

D

zy

qx

Here, dx, dy, and dz are the rigid body translations of N1  N2 elements due to external forces along three axes, respectively, sx and sy the relative slip distance parallel to the x and y axes, hi is the initial surface clearance, and h the surface gap after loading. The meanings of variables are illustrated in Fig. 2 graphically. The element centered at (2aD1, 2bD2) is represented by the notation [a, b]. The element in contact has a positive pressure and zero surface gap, while the non-contact element has zero pressure and an open gap. Thus contact element set Ac and non-contact element set Ac can be defined as,

ð8Þ

The general contact model (Johnson, 1985) is summarized as follows,

Ast ¼ f½a; bjjqj½a;b 6 lf p½a;b ; jsj½a;b ¼ 0g Asl ¼ f½a; bjjqj½a;b ¼ lf p½a;b ; jsj½a;b 6¼ 0g where jqj½a;b ¼

lf is the static friction coefficient. In the slip region, the

Normal Contact analysis p(x,y)

The stick-slip loop (CGM) qx(x,y), qy(x,y) ∇ ∇

Fx‘=

1

2

∑qx(x,y); Fy‘=

1

2

∑ qy(x,y)

Rigid body displacement adjustments δ x = δ x Fx /Fx‘; δ y = δ y Fy /Fy‘

No

ð11Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q2x½a;b þ q2y½a;b ; jsj½a;b ¼ s2x½a;b þ s2y½a;b , and

Normal and tangential loads W, Fx, and Fy

Surface geometry update by adding the normal displacements due to qx(x,y) and qy(x,y)

ð10Þ

In this study, a quasi-static process is assumed considering a slow-loading process. Motion initiation is simulated, where the micro-slips at the surface are relatively small as compared to contact dimensions. The static friction is used as the limit of the shear traction. In the stick region, the shear traction magnitude is less than the static friction, and the slip distance should vanish. We have the stick element set, Ast, and the slip element set, Asl, as,

∇ ∇

Here, fmn =

ð9Þ

uz þ hi  dz ¼ h

 D1 ; 2bD2  D2 Þ  fmn ð2aD1 þ D1 ; 2bD2  D2 Þ

RR

uy  dy ¼ sy

Yes Does pressure p(x,y) converge?

Fig. 3. Flowchart for the numerical contact simulation.

End

940

W.W. Chen, Q.J. Wang / Mechanics of Materials 40 (2008) 936–948

orientation of the micro-slip vector should be opposite to that of the shear traction vector. That is,

ðqx½a;b ; qy½a;b Þ  ðsx½a;b ; sy½a;b Þ ¼ jqj½a;b  jsj½a;b

rmn ðx; y; zÞ ¼

Z

1

1

½T xmn ðx  x0 ; y  y0 ; zÞqx ðx0 ; y0 Þ

þ T ymn ðx  x0 ; y  y0 ; zÞqy ðx0 ; y0 Þ

T

0

0

þ T zmn ðx  x0 ; y  y0 ; zÞpðx0 ; y0 Þdx dy ;

if ½a; b 2 Asl

ð12Þ

Asl \ Ast ¼ 0

m; n ¼ x; y; z

ð14Þ

Here, T xmn ; T ymn , and T zmn are stress solutions of the surface unit concentrated point force along x, y, and z directions, respectively. The exact formulas of Tmn can be found in Chapter 3 in Johnson (1985). The simulation domain is discretized into N3 layers along the z direction, and the space between adjacent layers is D3. All subsurface layers are meshed using the same rectangle element system as the

ð13Þ

The subsurface elastic stresses, rmn, at (x, y, z) can be associated with the surface tractions as,

a

1

1

In addition, the contact area is compose of the stick and slip regions,

Asl [ Ast ¼ Ac ;

Z

0.18 Fx=0.9μ fW

0.15

qx/p0

0.12 Fx=0.5μfW

0.09 0.06 0.03

Numerical solution Analytical solution

0.00 -1.0

-0.5

0.00

.5

1.0

x/a0

b 1.0

0.8

1.0

0.8

δx/δ0

δ x /δ 0

0.6

0.6 c/a0 Numerical solution Analytical solution

0.4

0.4

0.2

0.2

0.0

0.0 0.0

c/a0

Numerical solution Analytical solution

0.5

1.0

Fx /μ fW Fig. 4. Comparisons of the numerical results with analytical solutions, (a) shear traction qx along the x axis, and (b) variations of the stick zone radius and tangential rigid body deformation as a function of the tangential force.

941

W.W. Chen, Q.J. Wang / Mechanics of Materials 40 (2008) 936–948

b

a

0.4 1.0 μf /β =1.8

|q/p0

p/p0

Hertz solution μf /β =0.5

0.3

0.8 0.6 0.4

0.2

Hertz solution μ /β =0.5 f

0.2

0.1

μf /β =1.8

0.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0

0.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0

x/a0

x/a0

Fig. 5. Dissimilar-material contact under a normal load alone, (a) surface pressure along the x axis, and (b) shear traction magnitude jqj along the x axis.

where

Table 1 Comparisons of contact variables under different values of lf/b

lf/b

0.0 (Hz) 1.0000 1.0000 0.0000

pmax/p0 a/a0 c/a0

0.5 1.0811 0.9915 0.1853

1.0 1.0789 0.9839 0.4887

1.4 1.0770 0.9792 0.6248

Bkmn½a;b;c ¼

1.8 1.0761 0.9773 0.7060

Z

D1

Z

D1

k; m; n ¼ x; y; z; 0 6 b;

D2

D2

T kmn ð2aD1  x0 ; 2bD2  y0 ; cD3 Þ dx0 dy0 ;

0 6 a;

w 6 N2  1;

n 6 N1  1;

0 6 c 6 N3  1

ð15Þ

The closed-form expressions of indefinite double integrals of Tmn are listed in Appendix. topmost surface layer. rmn[a,b,c] are the elastic stresses of the element centered at the point of (2aD1, 2bD2, cD3), which can be evaluated as,

rmn½a;b;c ¼

2.2. Numerical procedure The inputs of the contact analysis include the initial surface geometry, hi, material properties, normal and tangential loads, and the static friction coefficient. First of all, the normal frictionless contact without tangential loads can be solved by the algorithm proposed in Polonsky and Keer

NX 1 1 N 2 1  X

qx½n;w Bxmn½an;bw;c þ qy½n;w Bymn½an;bw;c

n¼0

w¼0

þp½n;w Bzmn½an;bw;c

a



b 1

0.5

0.5

0.15

1

-0

0.0

-0.3

y a0

0.12

0.14

.02

0.3

-0.06 4 -0.0

0.05

-0.5

0.2

-0.1 -0.15 -0.2

0.1 0.08 0.06 0.04

- 0. 08

0

-0.05

-0.14

-0.1

-0.25

0

2 -0.1

y a0

0.1

-0.5

2

-1

-1 -1

-0.5

0

x a0

0.5

1

-1

-0.5

0

0.5

1

x a0

Fig. 6. Contours of the dimensionless shear traction qx on the surface for the contacts under a normal load alone, (a) lf/b = 0.5, and (b) lf/b = 1.8. (The regions inside the bold dashed lines are the stick zones, and those between the bold solid lines and the bold dashed lines are the slip zones.)

W.W. Chen, Q.J. Wang / Mechanics of Materials 40 (2008) 936–948

m ¼ x; y

ð16Þ

ð17Þ

Eq. (17) is actually a linear equation system with respect to unknowns qx and qy, which can be solved by the CG algorithm iteratively (Press et al., 1992). At the nth iterative ðnÞ ðnÞ step, the residual values of the equations, sx and sy , (i.e. the relative slip distances) are calculated from Eq. (17) ðnÞ ðnÞ using tractions in the current step, qx and qy . The DCFFT algorithm (Liu et al., 2000) is utilized to evaluate the linear convolutions existing in Eq. (17) efficiently. Then the square sum of the residual vectors is calculated, ðnÞ

sx½a;b

 ðnÞ þ sy½a;b

for ½a; b 2 Ast

ð18Þ

5 0.5

1

0.45

-1

b

0.64 p0

0.42a0

0.6 0. 55

0.5

0.5 0.3

-1

c

ðnþ1Þ

qy

"

# ðnÞ qx ðnÞ

qy

"  sðnÞ 

ðnÞ

tx

0

0.72a0

1

0.68 p0

0.74 p0

0 0.23a0

0. 12

ð20Þ

5

x a0

0. 62

0.5

0.57

z a0

i h iT ðnÞ ðnÞ ðnÞ ðnÞ sx sy  tx ty " # i Dxx Dxy h iT ðnÞ ðnÞ ðnÞ ty tx ty yx yy D D

0.4

0.2

0.45

0.52 0.47 0.42 7 0.3 2 0.3

0.

1

¼

0.63 p 0

0.6

1

Next, the shear tractions are updated by using a vector whose orientation is along the conjugate directions, as indicated in Eq. (21).

qx

0.18a0

67 0.

h

# ðnþ1Þ

1

1 0.

The updating length, s(n), along the conjugate directions is then computed as,

"

0

0

ðnÞ

ð19Þ

ðnÞ

5

b¼0 ðnÞ

tx

0.3

x a0

We can calculate the conjugate directions, t x and ty , for solution searching as, ( ðnÞ ðnÞ ðn1Þ sk½a;b þ ggðn1Þ t k½a;b for ½a; b 2 Ast ðnÞ t k½a;b ¼ ; where k ¼ x; y 0 for ½a; b 62 Ast

sðnÞ ¼ h

0.4

5

a¼0

2

2 

0.5

2

g ðnÞ ¼

 NX 1 1 N 2 1  X

0.5

0.2

Dyy

2 3 # qx       Dxz 6 7 dx sx 0 q  ¼ ¼ 4 5 y Dyz dy sy 0 p

0. 3

Dyx

0.48a0

0 .2

Dxy

1 0.

27

Dxx

0.62 p0

0

z a0

"

a

6 0.

where d0 = 3lfW  [(2  m1)(1 + m1)/E1 + (2  m2)(1 + m2)/E2]/ 8a0. Because the part of the contact area in stick is unknown in advance, an iterative method based on the conjugate gradient (CG) algorithm, which was used to solve the normal contact by Polonsky and Keer (1999), is employed here to determine the real stick zone. We initially assume that the entire contact area is in stick, thus the relative slip distances of all elements are zero. Plugging Eq. (8) into the first two equations in Eq. (9), we have,

0.2

dm ¼ d0 ½1  ð1  F m =lf WÞ2=3 ;

The boundary conditions in Eqs. (11) and (12) are checked to determine the stick and slip status for each element. If jqj[a,b] > lfp[a,b], the element [a, b] is removed from the stick element set Ast, and the shear traction magnitude over this element is set equal to the static friction, lfp[a,b]. Another examination is to check the orientations of the shear traction and micro-slip for the element in the slip set Asl. The element is drawn back to the stick element set Ast if the micro-slip is not opposite to the shear traction on this element. The (n + 1)th iterative step is performed based on the new stick element set and the new element shear tractions,

0. 3 0.2 5

(1999). The contact area and contact pressure, p, are then obtained. The initial rigid body horizontal translations, dx and dy, are estimated from specified tangential loads using the exact formula in Johnson (1985), which is derived under several restricted assumptions. However, it is a fairly good estimation to start the numerical iteration.

z a0

942

-1

0

1

x a0

# ð21Þ

ðnÞ

ty

ð0Þ

ð0Þ

Here, the values are initialized by setting qx ¼ qy ¼ 0; ð0Þ ð0Þ tx ¼ t y ¼ 0, and g(0) = 1.

Fig. 7. Contours of the dimensionless von-Mises stress in the y = 0 plane for the contacts under a normal load alone, (a) frictionless Hertz solution, (b) lf/b = 0.5, and (c) lf/b = 1.8. (Triangles mark the maximum subsurface stresses, and circles the maximum surface stresses. The maximum values are shown in the rectangular boxes.)

943

W.W. Chen, Q.J. Wang / Mechanics of Materials 40 (2008) 936–948 ðnþ1Þ

ðnþ1Þ

A new stick–slip loop needs to be performed using the updated rigid body translations. The computation ends until both normal pressure and shear tractions converge. Fig. 3 shows the flowchart of this numerical procedure.

qx and qy . The entire procedure from Eqs. (18)–(21) is repeated until the square sum of residual values g(n) is less than a specified tolerance error. Because of the dissimilarity of material properties, the tangential tractions will induce relative normal displacements and change the normal contact geometry and pressure distribution. Therefore, the normal contact needs to be re-simulated with the updated surface geometry. The total tangential forces in this step are calculated as,

F 0x

¼ D1 D2

NX 1 1 N 2 1 X

a¼0

qx½a;b ;

F 0y

¼ D1 D2

NX 1 1 N 2 1 X

a¼0

b¼0

qy½a;b

3. Verification In order to verify the current model, the numerical results of the point contact of two elastic materials are compared with the analytical solution given in Johnson (1985). An increasing tangential force, Fx, less than lfW that leads to a gross sliding, is applied along the x axis. The analytical solution is derived based on the assumptions that shear tractions have no mutual interaction with pressure and are parallel to the tangential force everywhere. The

ð22Þ

b¼0

The rigid body translations are adjusted based on the differences between the computed tangential forces ðF 0x ; F 0y Þ and the specified forces (Fx, Fy).

dx ¼ dx F x =F 0x ;

dy ¼ dy F y =F 0y

ð23Þ

Fx

a 1

0.7

0.6

0.2 Fx μ f W = 0.0

0.5

b

0.25

Ast /Ac

y a0

0.20 0

0.15 0.10

-0.5 0.05

0.9

0.8

0.00

-1 -1

-0.5

0

0.5

0.0

1

0.2

0.4

0.6

0.8

1.0

0.945

x a0

Fx /μfW

Fig. 8. (a) Evolution of the stick zone with the increasing shear force (The dashed lines are the boundaries of the stick zones, and the regions between the solid line and the dashed lines are the slip zones.), and (b) variation of the stick area as a function of the shear force. (The dimensionless shear force leading to a zero stick area is given in the rectangular box.)

b

a

0.3

1.0 0.2 0.1

0.6 0.4

Hertz solution Fx /μfW=0.2 Fx /μ fW=0.6

0.2

Fx /μ fW=0.9 Fx

0.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0

x/a0

q x /p0

p/p0

0.8

0.0 Fx /μ fW=0.0

-0.1

Fx /μ fW=0.2 Fx /μ fW=0.6

-0.2 -0.3

Fx /μfW=0.9 Fx

-1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0

x/a0

Fig. 9. Contact results affected by the increasing shear force, (a) surface pressure distribution along the x axis, and (b) shear traction qx distribution along the x axis.

944

W.W. Chen, Q.J. Wang / Mechanics of Materials 40 (2008) 936–948

solution of rigid body displacement is given in Eq. (16), and that of shear traction along the x axis is shown in Eq. (24).

4. Results and discussion

2sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi3  2 r 2 r c 5  Hðc  rÞ  qx ðx; yÞ ¼ lf p0 4 1  1 a0 a0 c

The current model is used to analyze the point contact of dissimilarly elastic materials with shear tractions on the interface. Suppose that a rigid ball indenter with a radius of R = 18 mm is compressed on an elastic half-space by a normal load P = 20 N. The half-space is assumed to be steel with elastic properties of E1 = 210 GPa, and m1 = 0.3. The resulting Hertz solutions are a0 = 0.105 mm and p0 = 860 MPa. A 256  256  32 mesh system is used here to discretize the physical contact domain with dimensions of 3a0  3a0  1.5a0. The contact initially takes place at the origin point. The directions of surface tractions correspond to those of tractions applied on the halfspace. In the following analyses, the stresses are normalized by p0, and the coordinates are normalized by a0.

ð24Þ Here, H(r) is Heaviside’s step function, c/a0 = (1  Fx/lfW)1/3 the stick zone radius, and a0 and p0 are the Hertz contact radius and peak pressure, respectively. The numerical verification is conducted using the same assumptions as those used in the analytical derivation. Note that these assumptions are released in Section 4. Fig. 4a presents the longitudinal traction, qx, along the x axis, and Fig. 4b shows variations of the stick radius, c, and the rigid body displacement, dx, with the increasing tangential load. Excellent agreements can be found between the current numerical results and the analytical solutions.

b

a

1 8

0.13

0.5

0.5 8

0.3

0

0

0.2

-0.12

-0.5

-0.05

0.23

0

0.28

0.1

2 -0.2 -0.17

y a0

0. 05

0.0

y a0

0.0

3

1

0.2 5 0.1

-0.5

-0.07

.0 -0

5

2

0.1

-1

-1 -1

-0.5

0

0.5

1

-1

-0.5

0

x a0

c

0.5

1

x a0 1

0. 01 05 0.

0.5

0.13

0.17 0.13

0

0.21

y a0

0.09

0.29 0.25

-0.5

0.2

1 0.17

-1 -1

-0.5

0

0.5

1

x a0 Fig. 10. Contours of the dimensionless shear traction along the x direction, qx, on the surface, (a) Fx/lfW = 0.2, (b) Fx/lfW = 0.6, and (c) Fx/lfW = 0.9. (The bold dashed lines enclose the stick zones, and the regions between the bold solid lines and the bold dashed lines are the slip zones.)

945

W.W. Chen, Q.J. Wang / Mechanics of Materials 40 (2008) 936–948

4.1. Contact under a normal load alone Because the contact bodies have dissimilar elastic constants, the normal contact pressure can cause the relative tangential deformation at the interface. Shear tractions are generated to oppose the tendency of this relative tangential deformation. This kind of interaction between the normal and shear tractions relies on the degree of dissimilarity of material properties, which can be quantified by the Dundurs constant, b.

b¼

ð1 þ m1 Þð1  2m1 Þ=2E1  ð1 þ m2 Þð1  2m2 Þ=2E2 ð1  m21 Þ=E1 þ ð1  m22 Þ=E2

ð25Þ

Here, the properties of contact bodies result in b = 0.2857. Note that for the contact of the same materials, b = 0. The coefficient of friction, lf, has a significant influence on the stick–slip phenomenon; however, this effect virtually depends on the discrepancy of material elastic properties. Therefore, the ratio of friction coefficient over the Dunders

b

1

0.16

0.1

0.08

0.2 0.2 0.24

0.5

y a

0

y a

0. 12

0.04

0

0.1

5

0.25

0.5

0.2

1

0.05

a

constant, lf/b, is used to characterize contact performances. Fig. 5 presents the profiles of pressure and shear traction magnitude along the x axis for lf/b = 0 (Hertz contact), 0.5, 1.8. The presence of shear tractions enhances the peak pressure, but the increase in lf/b reduces the peak pressure. This phenomenon can be explained by the distribution of the shear traction magnitude shown in Fig. 5b. A higher value of lf/b leads to a smaller shear traction near the center of the contact area, which has a smaller effect on the peak pressure. However, an interface with larger lf/b induces a larger shear traction in the peripheral zone of the contact area. The characteristic values of contact, including the peak pressure, the contact radius, and the stick zone radius, are listed in Table 1 with the increasing lf/b. As the value of lf/b increases, the contact radius decreases slightly, while the stick zone expands significantly. The contours of dimensionless shear traction along the x axis, qx, in the interface are plotted in Fig. 6; they are

0

0

-0.05 04 -0 .

-0.24 -0.2

-0.5

2

-0.15

-0.1

.1

-0.25 -0.2

-0 .1 6

-0

-0.5

-0.08

-1

-1 -1

-0.5

0

0.5

1

-1

-0.5

0

x a0

0.5

1

x a0

c

1

0.06

0.09

12 0.

0.5

0.03

y a

0

0.15

0 -0

-0.12

6 .0 -0

-0.5

09 -0 .

-0.15

.0 3

-1 -1

-0.5

0

0.5

1

x a0 Fig. 11. Contours of the dimensionless shear traction along the y direction, qy, on the surface, (a) Fx/lfW = 0.2, (b) Fx/lfW = 0.6, and (c) Fx/lfW = 0.9. (The bold dashed lines enclose the stick zones, and the regions between the bold solid lines and the bold dashed lines are the slip zones.)

W.W. Chen, Q.J. Wang / Mechanics of Materials 40 (2008) 936–948

0.72 p0

Stick zone

0

0.6

0.57

7

0. 12

0.0

7

0.19a0

0.62 0.65 p0

0.5

0.5

7 0.52

0.470.42 0.37

0. 3 0.2 2 7

1 0.17

-1

0

0.22

1

x a0

b

Stick zone

0.79 p0

0

0.68 0.65 p0

3

0.5

0.18

0.5

0.1a0

0.6

0.08

1

8 0.53 0.48

28 0.

0.23

0.43

0.3 0.38 3

-1

0

1

x a0

c

Stick zone 0.76 p0

0

0.64 0.38a0

1 29 0. 4

0.19

-1

0.64 0.59 0.54 0.49 4 0.4 9 0.3

0. 34

0.5

0.14 0.09

0.65 p0

0.2

Once the understanding of normal contact is established, the contact subjected to tangential loading can be readily investigated. This section assumes that an increasing tangential load, Fx, pushes the spherical indenter along the x axis, but the magnitude of Fx is not large enough to induce a gross sliding. In this part of study, the ratio of friction coefficient over the Dundurs constant, lf/b, is selected to be 1.0. Fig. 8a presents the evolution of the stick part of the contact area as the tangential load, Fx, increases gradually, and Fig. 8b shows the variation of the stick to contact area ratio as a function of the increasing shear force. When Fx/ lfW is less than 0.6, the stick zone is dragged by the tangential force towards the trailing edge of friction (it touches the contact border as Fx/lfW is about 0.6), and the stick area reduces steadily. Further increase in the tangential force makes the stick zone shrinks dramatically at the rear part of contact area until the stick zone vanishes completely. For the contact of the same materials, the gross sliding takes place at Fx = lfW because only the x-direction shear traction, qx, exists. Fig. 8b indicates that the stick zone reduces to zero at Fx/lfW = 0.945, thus a tangential force less than lfW will lead to a gross sliding due to the occurrence of the y-direction shear traction, qy. The profiles of pressure p and the x-direction shear traction, qx, along the x axis are plotted in Fig. 9 for different magnitudes of the tangential force. As shown in Fig. 9a, the increase in the tangential force reduces pressure at the contact edge in the friction leading direction and enhances that at the edge in the friction trailing direction. Variation of the peak pressure is trivial when Fx/ lfW < 0.6. After Fx increases beyond 0.6lfW, the peak pressure decreases, and its position is shifted opposite to the tangential force. The x-direction shear traction, qx, has both positive and negative values when the tangential load is less than 0.6lfW. The value of the negative part of qx

a

z a0

4.2. Contact with increasing tangential load

decreases with the increasing tangential force, and qx becomes entirely positive when the tangential force is larger than the transitional value of 0.6lfW. Figs. 10 and 11 give contours of the shear traction components, qx and qy, respectively. As indicated in these figures, qx contours are symmetric and while qy contours anti-symmetric with respect to the x axis. As the tangential force increases, the zone containing positive qx expands and that containing negative qx shrinks until qx in the entire contact area become positive (having the same direction of friction). The maximum positive and negative values of qx locate at the boundary of the stick zone, and their positions move along with the stick zone towards the trailing edge of friction. Fig. 11 validates the existence

z a0

anti-symmetric with respect to the y axis. qx tends to resist the x-direction tangential deformations induced by the normal pressure; these tangential deformations have opposite directions on left and right half planes. The maximum absolute values of the shear traction locate on the boundary of the stick zone. Considering the symmetry of point contact problems, the contour of qy can be obtained through folding that of qx with respect to the line of y = x. Fig. 7 shows the dimensionless von-Mises stress contours in the vertical cross section (y = 0) for different values of lf/b. The maximum equivalent stress intensity of 0.62p0 locates at the z = 0.48a0 below the origin for the Hertz pressure distribution. After taking into account the shear tractions induced by normal load W, the global maximums of the von-Mises stress occur on the surface instead of under the surface, which are marked by circles in Fig. 7. The increase in lf/b enhances the maximum von-Mises stresses and moves their positions away further from the center of the contact area. In addition, the subsurface maximum von-Mises stress (indicated by triangles) increases with lf/b, and its location is also uplifted towards the surface.

z a0

946

0

1

x a0 Fig. 12. Contours of the dimensionless von-Mises stress in the y = 0 plane, (a) Fx/lfW = 0.2, (b) Fx/lfW = 0.6, and (c) Fx/lfW = 0.9.

W.W. Chen, Q.J. Wang / Mechanics of Materials 40 (2008) 936–948

947

of the y-direction shear traction, qy, in dissimilar-material contacts subjected to the x-direction tangential load, Fx. The maximum qy decreases with the increasing tangential load. The positions of maximum qy are also at the interface of stick and slip zones, and they are shifted opposite to friction and close to the x axis. The dimensionless von-Mises stress contours in the vertical cross section (y = 0) are presented in Fig. 12 for the increasing tangential force, Fx. The local maximum stresses under the surface are marked by triangles, and those on the surface by circles. The corresponding stick zones are labeled on the figures. Fig. 12 reveals that the global maximum stress is at the interface of stick and slip zones and locates at the leading edge of stick zone. When Fx is below 0.6lfW, the increase in the tangential force enhances the surface maximum von-Mises stress and lifts the position of subsurface maximum stress. However, the increase in Fx beyond 0.6lfW actually decreases the intensity of surface maximum stress and lowers the position of subsurface maximum stress. In addition, the increase in the tangential force has a negligible effect on the intensity of the maximum subsurface von-Mises stress in the plane of symmetry.

concentrated point forces at the surface of half-space and is capable to identify the stick and slip status of the contact area. A good agreement with the analytical solution of a frictional point contact case well verifies the numerical accuracy of the present model. The analyses for the contacts between dissimilar materials under a pure normal load show the existence of shear tractions at the interface. The peak pressure is larger than that from the Hertz frictionless contact solution. Increase in the static friction coefficient leads to a higher peak shear traction, a larger stick region, and a higher maximum vonMises stress. An increasing tangential force, Fx, is applied to investigate the effect of friction on the contacts of dissimilar materials. The stick zone is dragged opposite to friction as Fx increases and eventually vanishes at the rear part of the contact area. The global maximum von-Mises stress is on the surface and locates at the leading interface of stick and slip zones. Because the y-direction shear traction, qy, is also induced at the surface, the gross sliding in the x direction can be reached at a lower tangential force than that predicted by the classic Coulomb friction law.

5. Conclusions

Acknowledgements

A three-dimensional numerical model for the contact of elastically dissimilar materials has been developed, where the coupling effect of normal pressure and shear tractions is considered. The model is based on the exact solutions of

The authors acknowledge research supports from US National Science Foundation, Office of Naval Research, and Department of Energy. The authors also thank Prof. Leon M. Keer for helpful discussions.

Appendix. Indefinite double integrals of the Green’s functions, Tmn

Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z Z

    1 y z 2mx2  2 lnðq þ yÞ þ zð1  2mÞ þ 2p qðq þ zÞ qðq þ yÞ qðq þ yÞ   Z Z 1 y 2 m y T xyy ðx; y; zÞ dx dy ¼ T yxx ðy; x; zÞ dx dy ¼ 2m lnðq þ yÞ  zð1  2mÞ  2p qðq þ zÞ q Z Z 1 T yzz ðy; x; zÞ dx dy ¼ ½z2 =qðq þ yÞ T xzz ðx; y; zÞ dx dy ¼ 2p   Z Z 1 x 2m x T yxy ðy; x; zÞ dx dy ¼ lnðq þ xÞ  zð1  2mÞ T xxy ðx; y; zÞ dx dy ¼  2p qðq þ zÞ q    Z Z 1 xz xy x y 1 T yz ðy; x; zÞ dx dy ¼  T xz ðx; y; zÞ dx dy ¼  tan 2p zq qðq þ yÞ Z Z 1 T yxz ðy; x; zÞ dx dy ¼ ðz=qÞ T xyz ðx; y; zÞ dx dy ¼ 2p       1 xy x xz  þ 2ð1  2mÞ tan1 T zxx ðx; y; zÞ dx dy ¼ 2m tan1 2p zq qþyþz qðq þ yÞ       1 xy y yz z 1 1 T yy ðx; y; zÞ dx dy ¼ þ 2ð1  2mÞ tan 2m tan  2p zq qþxþz qðq þ xÞ     1 xy xz yz z 1 þ  tan þ T zz ðx; y; zÞ dx dy ¼ 2p zq qðq þ yÞ qðq þ xÞ 1 ½ð2m  1Þ lnðq þ zÞ  z=q T zxy ðx; y; zÞ dx dy ¼ 2p 1 T zxz ðx; y; zÞ dx dy ¼ ½z2 =qðq þ yÞ 2p 1 ½z2 =qðq þ xÞ T zyz ðx; y; zÞ dx dy ¼ 2p

T xxx ðx; y; zÞ dx dy ¼

where q ¼

Z Z

T yyy ðy; x; zÞ dx dy ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 þ y2 þ z2 ; tan1 ðxy=zqÞjz¼0 ¼ sgnðxyÞp=2.

948

W.W. Chen, Q.J. Wang / Mechanics of Materials 40 (2008) 936–948

References Bjorklund, S., Andersson, S., 1994. A numerical method for real elastic contacts subjected to normal and tangential loading. Wear 179, 117– 122. Brizmer, V., Kligerman, Y., Etsion, I., 2006. The effect of contact conditions and material properties on the elasticity terminus of a spherical contact. Int. J. Solids Struct. 43, 5736–5749. Brizmer, V., Kligerman, Y., Etsion, I., 2007.Elastic–plastic spherical contact under combined normal and tangential loading in full stick. Tribol. Lett. 25 (1), 61–70. Cattaneo, C., 1938. Sul contatto di due corpi elstici: Distribuzione locale degli sforzi. Reconditi dell Accademia natzionale dei Lincei 27, 342– 348. 474–478, 434–436. Guyot, N., Kosior, F., Maurice, G., 2000. Coupling of finite elements and boundary elements methods for study of the frictional contact problem. Comput. Methods Appl. Mech. Eng. 181, 147–159. Hamilton, G.M., 1983. Explicit equations for the stresses beneath a sliding spherical contact. Proc. Inst. Mech. Eng. 197C, 53–59. Hills, D.A., Sackfield, A., 1987. The stress field induced by normal contact between dissimilar spheres. ASME J. Appl. Mech. 54, 8–14. Hills, D.A., Nowell, D., Sackfield, A., 1993. Mechanics of Elastic Contacts. Butterworth Heinemann Ltd., Oxford. Jacq, C., Nelias, D., Lormand, G., Girodin, D., 2002. Development of athreedimensional semi-analytical elastic–plastic contact code. ASME J. Tribol. 124, 653–667. James, R.D., Busby, H.R., 1995. Elasticity solutions for constant and linearly varying applied to a rectangular surface patch on the elastic half-space. J. Elasticity 38, 153–163. Johnson, K.L., 1985. Contact Mechanics. Cambridge University Press, London. Kogut, L., Etsion, I., 2003. A semi-analytical solution for the sliding inception of a spherical contact. ASME J. Tribol. 125, 499–506.

Kosior, F., Guyot, N., Maurice, G., 1999. Analysis of frictional contact problem using boundary element method and domain decomposition method. Int. J. Numer. Meth. Eng. 46, 65–82. Li, J., Berger, E.J., 2003. A semi-analytical approach to three-dimensional normal contact problems with friction. Comput. Mech. 30, 310–322. Liu, C.H., Lin, Y.-H., Lin, P.-H., 2007. A numerical analysis of partial slip problems under Hertzian contacts. Meccanica 42, 197–206. Liu, S.B., Wang, Q., Liu, G., 2000. A versatile method of discrete convolution and FFT (DC-FFT) for contact analyses. Wear 243, 101–111. Liu, S.B., Wang, Q., 2002. Study contact stress fields caused by surface tractions with a discrete convolution and fast Fourier transform algorithm. J. Tribol. 124, 36–45. Mindlin, R.D., 1949. Compliance of elastic bodies in contact. ASME J. Appl. Mech. 16, 259–268. Munisamy, R.L., Hills, D.A., 1992. A numerical analysis of an elastically dissimilar three-dimensional sliding contact. Proc. Inst. Mech. Eng. 206, 203–211. Nelias, D., Antaluca, E., Boucly, V., Cretu, S., 2007. A three-dimensional semianalytical model for elastic–plastic sliding contacts. ASME J. Tribol. 129, 761–771. Nowell, D., Hills, D.A., Sackfield, A., 1988. Contact of dissimilar elastic cylinders under normal and tangential loading. J. Mech. Phys. Solids 36 (1), 59–75. Polonsky, I.A., Keer, L.M., 1999. A numerical method for solving rough contact problems based on multi-level multi-summation and conjugate gradient techniques. Wear 231, 206–219. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes in FORTRAN: The Art of Scientific Computing, second ed. Cambridge University Press, London. Spence, D.A., 1973. Eigenvalue problem for elastic contact with finite friction. Proc. Camb. Philol. Soc. 73, 249–268. Spence, D.A., 1975. The Hertz contact problem with finite friction. J. Elasticity 5, 297–319.