Modeling of the contact between a rigid indenter and a heterogeneous viscoelastic material

Mechanics of Materials 77 (2014) 28–42

Contents lists available at ScienceDirect

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

Modeling of the contact between a rigid indenter and a heterogeneous viscoelastic material Koffi Espoir Koumi a,b, Daniel Nelias a,⇑, Thibaut Chaise a, Arnaud Duval a a b

Université de Lyon, INSA-Lyon, LaMCoS UMR CNRS 5259, F69621 Villeurbanne, France SNECMA, Centre de Villaroche, 77550 Moissy Cramayel, France

a r t i c l e

i n f o

Article history: Received 10 April 2014 Received in revised form 30 June 2014 Available online 10 July 2014 Keywords: Viscoelasticity Inhomogeneity Anisotropy Eigenstrain Eshelby’s equivalent inclusion method (EIM) Contact Mechanics

a b s t r a c t In this paper, the contact problem between a rigid indenter and a viscoelastic half space containing either isotropic or anisotropic elastic inhomogeneities is solved. The model presented here is 3D and based on semi-analytical methods. To take into account the viscoelastic properties of the matrix, contact and subsurface problem equations are discretized in the spatial and temporal dimensions. A conjugate gradient method and the fast Fourier transform are used to solve the normal problem, contact pressure, subsurface problem and real contact area simultaneously. The Eshelby’s formalism is applied at each step of the temporal discretization to account for the effect of the inhomogeneity on pressure distribution and subsurface stresses. This method can be seen as an enrichment technique where the enrichment fields from heterogeneous solutions are superimposed to the homogeneous viscoelastic problem solution. Note that both problems are fully coupled. The model is validated by comparison with a Finite Element Model. A parametric analysis of the effect of elastic properties and geometrical features of the inhomogeneity is proposed. The model allows to obtain the contact pressure distribution disturbed by the presence of inhomogeneities as well as subsurface and matrix/inhomogeneity interface stresses at every step of the temporal discretization. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Composite materials are more and more used in several domains of the engineering: aeronautics, automobiles . . .Some of these composite materials are made of by viscoelastic matrix and generally reinforced with anisotropic elastic fibers. Considering the complexity of these materials, it is difficult to develop a mathematical model allowing to analyze or predict their behavior within the framework of contact mechanics. Several authors have shown interest to the problem of contact between homogeneous viscoelastic materials. Lee and Radok (1960) obtained the solution of contact pressure and area for the spherical ⇑ Corresponding author. E-mail address: [email protected] (D. Nelias). http://dx.doi.org/10.1016/j.mechmat.2014.07.001 0167-6636/Ó 2014 Elsevier Ltd. All rights reserved.

indentation of a linear viscoelastic material. This solution is valid only for monotonical increase of the contact area and leads to a negative pressure when the contact area decreases. Since then, several approaches were proposed in the literature to overcome this issue. Hunter (1960) and Graham (1967) introduced methods for viscoelastic indentation test which are able to handle the case when the contact radius possesses a single maximum. Ting (1966, 1968) dealt with the indentation problem in which the time dependent contact area is an arbitrary function of time, and the indenter has an axisymmetric profile. Most of the solutions presented above are based on complex analytical formalisms, and limited to particular geometries (cone, sphere) and to the ideal viscoelastic material with only one relaxation time. Several other authors were interested in the indentation problem, rolling, sliding and

K.E. Koumi et al. / Mechanics of Materials 77 (2014) 28–42

29

Nomenclature

Letters a contact radius a1 ; a2 ; a3 semi-axes of an ellipsoidal inhomogeneity Bijkl influence coefficients that relating the stress rij at point ðx31 ; x2 ; x3 Þ to the constant eigenstrain at the point ðxk1 ; xk2 ; xk3 Þ I CM ; C ijkl ijkl elastic constants of the matrix and the inhomogeneity EI Young’s modulus of the inhomogeneity h distance between the two surfaces of the contacting bodies HðtÞ the Heaviside step function Iijkl the fourth-order identity tensor JðtÞ viscoelastic creep function Kn coefficients in the normal displacement at the contact surface due to the contact pressure L1 ; L2 ; L3 lengths of the three sides of the matrix in Finite Element Model M ij influence coefficients relating the stress rij at the point ðx1 ; x2 ; x3 Þ to the normal traction rn within a discretized area centered at ðxk1 ; xk2 ; 0Þ n1 ; n2 ; n3 grid sizes in the half-space along the Cartesian directions x1 ; x2 ; x3 , respectively P normal applied load P0 maximum Hertzian pressure p contact pressure distribution D indenter diameter RðtÞ viscoelastic relaxation function Sijkl components of the Eshelby’s tensor u0i ui W dx3

displacements corresponding to the infinite applied strain e0ij disturbed contribution of the displacements applied exterior load depth of the inclusion from the surface of the matrix in EF model

rolling friction for viscoelastic materials. This is the case of Argatov (2012) who gave the analytical solution of the rebound indentation problem for an isotropic linear viscoelastic layer loaded with a spherical punch. His model is only valid for the indentation phase with a monotonic loading. A contact between an axisymmetric indenter and a viscoelastic half-space is presented by Greenwood (2010). Considering the sliding contact or rolling friction in viscoelastic contact problem; Hunter (1961) and later, Goriacheva (1973) proposed a solution for the rolling contact between a rigid indenter and a viscoelastic half-space. Both approaches have in common to limit themselves to particular shape (cylindrical) of the indenter and to an ideal viscoelastic material. To get rid of the indenter shape limitation, Vollebregt (2011) proposed a model based on boundary element method. His model can only account for viscoelastic material with one single relaxation time. Recently, Carbone and Putignano (2013) introduced a novel methodology to investigate steady-state viscoelastic sliding or rolling contact based on boundary element

xI ¼ ðxI1 ; xI2 ; xI3 Þ Cartesian coordinates of the inclusion center Greek letters e0ij infinite applied strain eij strain due to eigenstrains eij eigenstrain due to the presence of inhomogeneities r0ij stress corresponding to the infinite applied strain e0ij rij disturbed contribution of the stresses /; W harmonic and biharmonic potentials of mass density eij dij Kronecker symbol rn normal pressure due to the summation of both symmetric inclusions Dx1 ; Dx2 half size of the discretized surface area mM ; mI Poisson’s ratio of the matrix M and the inclusion I c the ratio of the inhomogeneity Young’s modulus to the matrix g the dashpot viscosity s the relaxation time h the tilted angle of the inhomogeneity in the x1 Ox3 plan Acronyms and fast Fourier transforms 2D-FFT two-dimensional fast Fourier transform 3D-FFT three-dimensional fast Fourier transform FFT 1 inverse FFT operation b ijkl frequency response of coefficients Bijkl in the B frequency domain b ij M frequency response of coefficients M ij in the frequency domain

method. Their approach is able to deal with a large spectra of relaxation times for the viscoelastic materials. Other authors (Nasdala et al., 1998; Le and Rahler, 1994; Nackenhorst, 2004; Padovan et al., 1992, 1984) solved the sliding/rolling friction problem of viscoelastic contact using the Finite Element Method. This method is able to account for the real viscoelastic materials and any geometry of contacting bodies. However using the Finite Element Method, the accuracy of the contact solution in terms of pressure and subsurface stresses is often insufficient. A robust semi-analytic method for contact between a rigid indenter and a viscoelastic half-space has first been introduced by Chen et al. (2011). This semi-analytical approach can account: for a wide spectra of relaxation times for linearly viscoelastic materials, an arbitrary loading profile (with possibility of decreasing the contact area); it can also simulate contact with a rough surface by incorporating the asperity heights into initial surface gap. The model presented in this paper is based on the semi-analytical method introduced by Chen et al. (2011) for viscoelastic

30

K.E. Koumi et al. / Mechanics of Materials 77 (2014) 28–42

materials, also accounting for the presence of isotropic or anisotropic ellipsoidal inhomogeneities of any orientation within the matrix (Koumi et al., 2014). Many authors already investigated stress concentration due to the presence of inclusions or inhomogeneities embedded in purely elastic domain. However until very recently they did not solve the contact problem, instead they assumed the hertzian contact pressure distribution as input, see for example Kabo and Ekberg (2002, 2005) or Courbon et al. (2005). For a 2D contact problem the effect of the presence of an elliptical heterogeneous inclusion on the contact loading was first solved numerically by Kuo (2007, 2008) using the boundary element method (BEM). The effect of inhomogeneities on the three-dimensional contact problem solution was first investigated by Nelias and co-workers (Leroux et al., 2010; Leroux and Nelias, 2011) and then by Zhou et al. (2009, 2011a,b, 2012) based on the semi-analytical method (SAM) initially proposed by Jacq et al. (2002) to numerically solve 3D elastic–plastic contact. For more details on the modeling of inclusions in contact problems the reader may refer to the review paper by Zhou et al. (2013). Recently, Koumi et al. (2014) presented a model to account for the effect of ellipsoidal anisotropic inhomogeneities on contact problem. Besides field pressure and the subsurface stresses, the normal and shear stresses at the interface between the inclusion and the matrix have also been investigated. As sus-mentioned the approaches introduced by these authors deal with a purely elastic half or infinite space. To the best of the authors knowledge, no theoretical studies have been published on the contact between a rigid indenter and an heterogeneous viscoelastic material containing one or two inhomogeneities. The model presented here, will propose a new approach to solve the contact problem in the presence of an heterogeneous viscoelastic half-space. The approach introduced by Chen et al. (2011) for the contact between a rigid sphere and an homogeneous viscoelastic solid will be used. A temporal and a spatial discretization of equations governing the contact is performed. The Eshelby’s formalism is applied at every step of the temporal discretization in order to take into account the effect of inhomogeneities. The analysis will be limited to the contact between a spherical indenter and an heterogeneous viscoelastic half space containing one or two interacting inhomogeneities for comprehension. This model is able to take into account any number of relaxations times and to manage a wide range of semi-infinite geometries, potentially including rough surfaces, isotropic and anisotropic inhomogeneities.

in time of these types of materials, the relaxation function defined by Eq. (1) is introduced.

rðtÞ e0

RðtÞ ¼

ð1Þ

Likewise, the creep compliance function JðtÞ is defined to be the ratio of the strain history over the constant stress.

JðtÞ ¼

eðtÞ r0

ð2Þ

The creep and relaxation function are linked by the following expression:

Z

t

JðnÞRðt  nÞdn ¼ t

ð3Þ

0

Based on the generalized Maxwell model (see Fig. 1), the relaxation function can be written as follows:

"

l0 þ

RðtÞ ¼

#

n X

li expðt=si Þ HðtÞ

Eq. (4) is also known as the Prony series decomposition. Under the assumption of linear viscoelasticity, the stress and deformation fields can be written as follows:

rðtÞ ¼

Z

t

eðtÞ ¼

deðnÞ dn dn

ð5Þ

drðnÞ dn dn

ð6Þ

Rðt  nÞ 0

Z

t

Jðt  nÞ

0

2.2. Contact between a rigid indenter and an homogeneous viscoelastic material The theory presented in this part was developed by Chen et al. (2011). It is resumed briefly here in order to facilitate the first understanding of the rest of the paper. A normal contact problem between two contacting bodies (see Fig. 2) can be described by a set of equations that must be solved simultaneously. These equations are:  The load balance. The applied load WðtÞ and the integration of the contact pressure pðx1 ; x2 ; tÞ in the real contact area at time t Cc ðtÞ must be strictly equal.

WðtÞ ¼

Z

pðx1 ; x2 ; tÞdx1 dx2

Cc ðtÞ

2. Theoretical background and model description 2.1. Some basic principles on viscoelastic materials Viscoelasticity is the property of certain materials that exhibit both viscous and elastic behaviors under mechanical solicitation. Whether under creep or relaxation solicitations, the mechanical response of a viscoelastic material is time-dependent. In order to characterize the dependence

ð4Þ

i¼0

Fig. 1. Generalized Maxwell model.

ð7Þ

K.E. Koumi et al. / Mechanics of Materials 77 (2014) 28–42

P

31

3. Theoretical Background – Eshelby’s equivalent inclusion method in contact mechanics Accounting for the presence of inhomogeneities in the contacting bodies (see Fig. 3) consists in adding in Eq. (8) the eigendisplacement u3 induced by inhomogeneities. Details on the calculation procedure will be given later in Section 3.3. Eq. (8) is then modified as follows:

o

X2

X1

hðx1 ; x2 ; tÞ ¼ hi ðx1 ; x2 Þ þ d þ u3 ðx1 ; x2 ; tÞ þ u3 ðx1 ; x2 Þ

( ),

X3

Fig. 2. Contact of rigid indenter over a viscoelastic half-space.

 The surface separation. The gap between the two contacting surfaces hðx1 ; x2 ; tÞ is defined by the initial geometry hi ðx1 ; x2 Þ, the rigid body displacement d and by the normal surface displacement of the viscoelastic half space u3 ðx1 ; x2 ; tÞ

hðx1 ; x2 ; tÞ ¼ hi ðx1 ; x2 Þ þ d þ u3 ðx1 ; x2 ; tÞ

Z

Z

1

1

1

1

t  nÞ

Z

t 0

Gðx1  x01 ; x2  x02 ; x3 ;

@pðx01 ; x02 ; tÞ 0 0 dx1 x2 dn @n

ð9Þ

where, Gðx1 ; x2 ; tÞ is the viscoelastic Green’s function. Gðx1 ; x2 ; tÞ can be decomposed into a spatial part and a temporal part.

ð1  mÞ Gðx1 ; x2 ; tÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi JðtÞ ¼ Gðx1 ; x2 ÞJðtÞ p x21 þ x22 u3 ðx1 ; x2 ; aDtÞ ¼

Z

1

Z

1



1 1

ð10Þ

0

a X J½ða  kÞDt  ½pðx01 ; x02 ; kÞ

An infinite matrix D with the elastic stiffness tensor CM ijkl ðaDtÞ containing an ellipsoidal domain X with the elastic stiffness tensor C Iijkl is submitted at infinity to a uniform strain e0 . The strain field is disturbed by the presence of the inhomogeneity. The Eshelby’s equivalent inclusion method (EIM) consists in representing the ellipsoidal inhomogeneity as an inclusion having the same elastic properties C M ijkl ðaDtÞ as the matrix but being subjected to an additional imaginary strain called eigenstrain e giving: 0  C Iijkl ðe0kl þ ekl Þ ¼ C M in X ijkl ðaDtÞ  ðekl þ ekl  ekl Þ

eij ¼ Sijkl ðaDtÞ  ekl ;

k¼0

ð11Þ

 The contact conditions. The distance hðx1 ; x2 ; tÞ is always positive, because the contacting bodies can not interpenetrate each other. The conditions are defined by the inequalities:

ð14Þ

ð15Þ

where Sijkl ðaDtÞ is the Eshelby’s tensor at t ¼ aDt. The Eshelby’s tensor depends only on the Poisson’s ratio of the matrix. In the model, this Poisson’s ratio is assumed to be constant. It is shown that

Sijkl ðaDtÞ ¼ Sijkl ðt ¼ 0Þ ¼ Sijkl

ð16Þ

Substitution of Eq. (15) into Eq. (14) leads to:

contact : hðx1 ; x2 ; tÞ ¼ 0 and pðx1 ; x2 ; tÞ > 0 in Cc ðtÞ

 0 DC ijkl ðaDtÞSklmn emn þ C M ijkl ðaDtÞekl ¼ DC ijkl ðaDtÞekl

where

separation : hðx1 ; x2 ; tÞ > 0 and ðx1 ; x2 ; tÞ ¼ 0 out Cc ðtÞ

3.1. Eshelby’s solution for an infinite space

The necessary and sufficient condition for the equivalence of the stresses and strains in the two above problems of inhomogeneity and inclusion is provided by Eq. (14). In particular, the eigenstrain eij is related to compatibility strain eij by:

Gðx1  x01 ; x2  x02 Þdx1 x02

 pðx01 ; x02 ; k  1Þ

The algorithm presented in Fig. 4 is used to couple the contact problem with the subsurface inhomogeneous problem on the field u3 . The basics of the method in the case of a semi-infinite elastic half-space are detailed extensively in Leroux et al. (2010), Leroux and Nelias (2011) and Koumi et al. (2014). Sections 3.1–3.3 present the various steps of calculation of u3 in the case of a viscoelastic body containing either isotropic or anisotropic elastic inhomogeneities. The calculation of the eigenstrain e in an infinite space is presented first. The calculation of the eigenstrain e in a semi-infinite space is then detailed. And finally the calculation of the eigendisplacement u3 is presented. Eshelby’s formalism is applied at every step of the temporal discretization in order to take into account the effect of inhomogeneities. C M ijkl ðaDtÞ represents the elastic stiffness tensor of the viscoelastic matrix at time t ¼ aDt.

ð8Þ

where hi ðx1 ; x2 Þ is the initial geometry, d the rigid body displacement. The surface normal displacement, u3 ðx1 ; x2 ; tÞ, at time t of a viscoelastic surface can be expressed as:

u3 ðx1 ; x2 ; tÞ ¼

ð13Þ

ð12Þ

DC ijkl ðaDtÞ ¼ C Iijkl  C M ijkl ðaDtÞ

ð17Þ

32

K.E. Koumi et al. / Mechanics of Materials 77 (2014) 28–42

A-A

(b)

(a)

P

P

o

X2

A

A 1

( ), 2

∗

2

X1

dx3

o

M

x1 θ

3

2

X3

2

3

1

( ),

x3

Fig. 3. Contact of a rigid spherical indenter over a viscoelastic half-space in the presence of an arbitrarily oriented inclusion. (a) 3D configuration and (b) 2D configuration.

3.2. Half-space solution Three dimensional contact problems involve a halfspace that is bounded by the surface plane x3 ¼ 0 in the cartesian coordinate system (x1 ; x2 ; x3 ). Jacq et al. (2002) and later Zhou et al. (2009) proposed a method allowing to extend the previous solution, valid only for infinite spaces, to half spaces. The solution consists in decomposing the problem into three subproblems (see Fig. 5), known as Chiu’s decomposition (Chiu, 1978). (1) An inclusion with the prescribed eigenstrain e ¼ ðe11 ; e22 ; e33 ; e12 ; e13 ; e23 Þ in an infinite space. (2) A symmetric inclusion with a mirror eigenstrain es ¼ ðe11 ; e22 ; e33 ; e12 ; e13 ; e23 Þ in the same space. (3) A normal traction distribution rn at the surface of the half space (x3 ¼ 0) which is a function of the eigenstrains e and es . The summation of the two solutions (1) and (2) leaves the plane of symmetry (x3 ¼ 0) free of shear tractions. By adding an opposite normal stress rn , the condition of free surface traction is obtained. The stress at any point of the domain meshed with n1  n2  n3 cuboids is given by:

rij ðx1 ; x2 ; x3 ; aDtÞ ¼

nX 3 1n 2 1n 1 1 X X

Bijkl ðx1  xI1 ; x2  xI2 ; x3  xI3 ; aDtÞ

where Bijkl are the influence coefficients that relate the constant eigenstrain at the point ðxI1 ; xI2 ; xI3 Þ which is the inclusion center in an infinite space to the stress rij at the point ðx1 ; x2 ; x3 Þ. M ij represents the influence coefficients relating the normal traction rn within a discretized area centered at ðxI1 ; xI2 ; 0Þ to the stress rij at the point ðx1 ; x2 ; x3 Þ.

Bijkl ðx; aDtÞ ¼ C M for x in  X ijmn ðaDtÞDmnkl ðx; aDtÞ

Bijkl ðx; aDtÞ ¼ C M for x in X ijmn ðaDtÞðDmnkl ðx; aDtÞ  Imnkl Þ ð20Þ where Iijkl ¼ 12 ðdil djk þ dik djlÞ is the fourth-order identity tensor. The second-order tensor M ij (see Appendix B) and the fourth-order tensor Dijkl depends only of the Poisson’s ratio which is assumed to be constant. Thus

Dijkl ðx; aDtÞ ¼ Dijkl ðxÞ

xI3 ¼0 xI2 ¼0 xI1 ¼0

Mij ðx; aDtÞ ¼ M ij ðxÞ

Dijkl ¼

1 ½W;ijkl  2mdkl /;ij  ð1  mÞðdkl /il 8pð1  mÞ þ dki /;jl þ djl /;ik þ dli /;jk Þ

WðxÞ ¼

Z

ð18Þ

ð23Þ

0

jx  x0 jdx

X

Z X

xI2 ¼0 xI1 ¼0

 rn ðxI1 ; xI2 ; 0; aDtÞ

ð22Þ

The expression for Dijkl is given in Mura (1987).

/ðxÞ ¼

 eskl ðxI1 ; xI2 ; xI3 ; aDtÞ nX 2 1n 1 1 X Mij ðx1  xI1 ; x2  xI2 ; x3 ; aDtÞ 

ð21Þ

and

xI3 ¼0 xI2 ¼0 xI1 ¼0

 ekl ðxI1 ; xI2 ; xI3 ; aDtÞ nX 3 1n 2 1n 1 1 X X Bijkl ðx1  xI1 ; x2  xI2 ; x3 þ xI3 ; aDtÞ þ

ð19Þ

1 0 dx jx  x0 j

For a single inclusion centered at ðxI1 ; xI2 ; xI3 Þ in the halfspace, the normal traction rn at the surface point ðx01 ; x02 ; 0Þ is obtained as:

K.E. Koumi et al. / Mechanics of Materials 77 (2014) 28–42

33

Fig. 4. Flowchart of the algorithm for the coupling between viscoelastic contact problem and inhomogeneity problem.

rn ðx01 ; x02 ; 0; aDtÞ ¼ B33kl ðx01  xI1 ; x02  xI2 ; xI3 ; aDtÞekl ðxI1 ; xI2 ; xI3 ; aDtÞ

In Eq. (18), each component Mij ðÞ is obtained by a double integration of the function F ij ðÞ over a discretized surface area 2Dx1  2Dx2 centered at ðxI1 ; xI2 ; 0Þ, see Appendices A and B.

 B33kl ðx01  xI1 ; x02  xI2 ; xI3 ; aDtÞeskl ðxI1 ; xI2 ; xI3 ; aDtÞ

ð24Þ

34

K.E. Koumi et al. / Mechanics of Materials 77 (2014) 28–42

∗

,

O

,

O

=

O

+

,

-

O

,

∗

,

,

,

,

Fig. 5. Decomposition of the half-space solution into three sub-problems.

Mij ðx1  xI1 ; x2  xI2 ; x3 Þ ¼

Z

xI1 þDx1

xI1 Dx1

Z

xI2 þDx2

xI2 Dx2

4.1. Material properties of the viscoelastic half-space

F ij ðx1 0

 x01 ; x2  x02 ; x3 Þdx1 x02

ð25Þ

The 3D-FFT is used to accelerate the calculation of the first (1) and second terms (2) and the 2D-FFT for the third term (3). Wrap around order and zero-padding techniques are used in order to remove the induced periodicity error (Liu et al., 2000).

The surface normal eigendisplacement u3 can be obtained when inserting the eigenstrain into the total strain. They are generated by the pressure field rn only. The normal displacements are calculated as:

u3 ðx1 ; x2 ; aDtÞ ¼

JðtÞ ¼

1

þ

t

ð27Þ

l g

and:

3.3. Normal displacement of a surface point

nX 2 1n 1 1 X

The same material properties are considered for the viscoelastic half-space in Sections 4.2 and 4.3. The basic Maxwell viscoelastic model with one relaxation time is considered (see Figs. 6 and 7). The creep and relaxation function are respectively:

K n ðx1  x01 ; x2  x02 Þrn ðx01 ; x02 ; aDtÞ

x02 ¼0 x01 ¼0

RðtÞ ¼ l expðt=sÞ

ð28Þ

where l is the spring stiffness, g the dashpot viscosity, and s ¼ g=l the relaxation time. The viscoelastic half space is assumed to have a timeindependent Poisson’s ratio mM =0.3. The instantaneous shear modulus is Rð0Þ ¼ l ¼ 80:77 Gpa and the relaxation time s ¼ 25 s.

ð26Þ To solve the equation above numerically, the surface in contact is discretized into n1  n2 rectangular elements of uniform size 2Dx1  2Dx2 . Then, pressure and displacement within each discrete patch are treated as constant and their values located at the center. The effect of an uniform pressure on a rectangular area has been given by Love (1952) and Johnson (1985). K n ðx1  x01 ; x2  x02 Þ denotes the influence coefficients that relate the normal pressure at the surface point ðx01 ; x02 ; 0Þ to the normal displacement at the surface point ðx1 ; x2 ; 0Þ, recalled in Appendix C. 4. Validation and results A normal contact between a spherical rigid indenter and a viscoelastic half space containing one or two inhomogeneities is considered. The diameter of the rigid indenter is D ¼ 62 mm and a normal force P ¼ 10000 N is applied. a and P0 are the elastic solutions of contact radius and peak pressure for the equivalent homogeneous elastic problem when the instantaneous shear modulus Rð0Þ ¼ l is used.

4.2. Validation by Finite Element Analysis In order to validate the semi-analytical method, a comparison with a finite element (FE) model was performed using the commercial FE package Abaqus v6.11–2. The configuration is presented in Fig. 8. A contact between a rigid sphere and a plane (solid) containing an ellipsoidal inclusion at a certain depth (dx3 ) is described. In order to correctly describe the interface, the inclusion is obtained by making a solid partition in a local coordinate system related to the contact coordinate, via a rotation of 30° around the x2 axis (see Fig. 9). The size of the sphere and the solid are recalled in Table 1. The geometrical properties and position of the inclusion are normalized by the contact half-width a , as shown in Table 2.

Fig. 6. Maxwell viscoelastic model.

35

K.E. Koumi et al. / Mechanics of Materials 77 (2014) 28–42

(a)

(b)

2

1

1.75

0.9 0.8

1.5

0.7 0.6

1.25

0.5 0.4

1

0.3 0.2

0.75

0.1 0 0

0.25

0.5

0.75

1

0.5 0

0.25

0.5

0.75

1

Fig. 7. Relaxation function (a) and Creep function (b).

 Mesh of the Finite Element Model. Since the model is symmetric with respect to the plane x2 ¼ 0, only half of the bodies in contact are meshed. A local area called ’subsurface layer’ is meshed by 30; 250 linear hexahedral elements (C3D8) of size 0:02a to correctly describe the contact. The rest of the solid is meshed with 1; 088; 349 linear tetrahedral (C3D4) elements. The (quarter) sphere contains 39; 671 quadratic tetrahedral (C3D10) elements and the inclusion 11; 009 linear tetrahedral (C3D4) elements.  Mesh of the Semi-Analytical Model. The mesh is limited here to the area of interest ½2a; 2a  ½2a; 2a  ½0; 2a. The meshing of the boundary conditions is not necessary. Note that it is one of the main advantages of the semi analytical methods. The area of interest (ie ½2a; 2a ½2a; 2a  ½0; 2a) is meshed with 281,775 cuboidal elements of same dimension. Two calculations are carried out. The first case corresponds to a contact between a rigid indenter and an homogeneous viscoelastic half-space. The second case considers a contact between a rigid indenter and an heterogeneous viscoelastic half space containing an elastic isotropic tilted inclusion. The inhomogeneity material properties are shown in Table 3. Results from the Semi Analytical Method (SAM) and the Finite Element Model (FEM) are presented in Fig. 10 for the two cases studied. The comparison focuses on the pressure distribution obtained by both methods for t ¼ 0; t ¼ s=2; t ¼ s. An excellent agreement is observed between both methods thus validating the SAM framework.

4.3. Results The simulation in the time domain varies in ½0; s and is divided into 40 time steps. a ðtÞ represents the

instantaneous contact radius corresponding to the homogeneous viscoelastic case. a ðtÞ increases with t. The instantaneous Young’s modulus of the viscoelastic half space at the initial time is given by Eq. (29).

EM ð0Þ ¼ 2  ð1 þ mM Þ  Rð0Þ

ð29Þ

The ratio c can be expressed as:

c ¼ EI =EM ð0Þ where EI is the inhomogeneity.

ð30Þ Young’s

modulus

of

the

elastic

4.3.1. Isotropic inhomogeneity Figs. 11–14 present the normalized contact pressure distribution respectively for a ratio c equal to 1; 2; 0:5 and 0:2. An ellipsoidal tilted inhomogeneity with semiaxes a1 ¼ 0:4a ; a2 ¼ a3 ¼ 0:1a located at dx3 ¼ 0:4a with an angle of h ¼ 30 is considered. Results are plotted in the plane x2 ¼ 0. In the first case (see Fig. 11) when c ¼ 1 it can be observed that at t ¼ 0 the inhomogeneity has no effect on the contact pressure distribution. However at t ¼ s=2 and at t ¼ s the peak magnitude pressure increases respectively by 6% and 18% compared with the homogeneous viscoelastic case. This result is consistent due to the fact that the substrate material surrounding the inhomogeneity becomes softer when the time increases while the inhomogeneity remains elastic. The case of a stiff inhomogeneity (c ¼ 2) is presented in Fig. 12. It leads to a peak magnitude augmentation respectively of 6%; 14% and 24% at t ¼ 0; t ¼ s=2; t ¼ s. When c ¼ 0:5 (Fig. 13), the contact pressure fields obtained in the homogeneous viscoelastic half-space and in the heterogeneous viscoelastic half-space are the same for t ¼ s=2. At t = 0 a peak pressure magnitude reduction of 9:5% is obtained while t ¼ s corresponds to a pressure augmentation of 10:5%.

36

K.E. Koumi et al. / Mechanics of Materials 77 (2014) 28–42

Sphere

Subsurface layer

Solid

Inclusion

Fig. 8. Finite Element Model used for the validation.

30o

Table 3 Inhomogeneity material properties.

1

1

Region

Material properties

Inhomogeneity

EI ¼ 840 GPa

mI ¼ 0:3

30o 2

=

, 2 3

, 3

Fig. 9. Local coordinate system for the model partition; Reference coordinate system ðx1 ; x2 ; x3 Þ and inhomogeneity coordinate system ðx01 ; x02 ; x03 Þ.

Table 1 Geometry of the bodies in contact. Region

Geometry (mm)

Sphere Solid

R = 31 L1 ¼ L2 ¼ L3 ¼ 60

Table 2 Size and location of the ellipsoidal inclusion (reference case). Geometry

Position

a1 ¼ 0:4a ; a2 ¼ a3 ¼ 0:1a

dx3 ¼ 0:4a

For c ¼ 0:2 (Fig. 14) the peak pressure decreases by 17:15% when t ¼ 0 and by 12% At t ¼ s the contact pressure distribution is same as in the case of the homogeneous half-space and in the case of heterogeneous half-space.

magnitude for t ¼ s=2. almost the viscoelastic viscoelastic

4.3.2. Two spherical inhomogeneities In this section the case of two spherical inhomogeneities S1 and S2 centered respectively at (dx1 ¼ 0:3a ; dx2 ¼ 0; dx3 ¼ 0:4a ) and (dx1 ¼ 0:3a ; dx2 ¼ 0; dx3 ¼ 0:4a ) is presented. The two inhomogeneities have the same radii r ¼ 0:14a . The elastic properties are such that c ¼ 0:6 for S1 and 0:8 for S2 . The conjugate gradient method (CGM) algorithm introduced by Zhou et al. (2011) is used to take into account the mutual influence between the two inhomogeneities. One can remark that when t ¼ 0 both S1 and S2 behave as soft inhomogeneities. At t ¼ s=2; S1 behaves as a soft inhomogeneity and S2 as a hard inhomogeneity. When

37

K.E. Koumi et al. / Mechanics of Materials 77 (2014) 28–42

(a)

(b)

1.2

1.2

t=0

t=0

1

1

t=τ /2

0.8 0.6

t=τ /2

0.8

t=τ

0.6

0.4

0.4

0.2

0.2

0 −2

−1.5

−1

−0.5

0

0.5

1

1.5

0 −2

2

t=τ

−1.5

−1

−0.5

0

0.5

1

1.5

2

Fig. 10. Comparison of the normalized pressure distribution between the Semi Analytical Method (SAM) and the Finite Element Model (FEM); (a) Homogeneous viscoelastic half-space, (b) Heterogeneous viscoelastic half-space containing isotropic tilted inhomogeneity (dx3 =a ¼ 0:4; a1 ¼ 0:4a ; a2 ¼ a3 ¼ 0:1a ; h=30°).

(a)

(b)

1.2

1.2

t=0

t=0

1

1 t=τ /2

0.8 0.6

0.6 t=τ

0.4

t=τ

0.4

0.2 0 −2

t=τ /2

0.8

0.2

−1.5

−1

−0.5

0

0.5

1

1.5

2

0 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Fig. 11. Distribution of normalized contact pressure for c ¼ EI =EM ð0Þ ¼ 1; mI ¼ 0:3 (a1 ¼ 0:4a ; a2 ¼ a3 ¼ 0:1a ; dx3 ¼ 0:4a ; h ¼ 30 ); (a) x1 ¼ 0 and (b) x2 ¼ 0.

t ¼ s both inhomogeneities have the same effect as two hard inhomogeneities. Fig. 15.

In the plane x2 ¼ 0 the peak magnitude increases from 2:1% (at t ¼ 0) to 49% (at t ¼ s).

4.3.3. Anisotropic inhomogeneity The case of an anisotropic inhomogeneity is also presented (see Fig. 16). The material properties are given as follows: EI1 ¼ 210 Gpa, EI2 ¼ 623 Gpa, EI3 ¼ 50 Gpa, lI12 ¼ 83 Gpa, lI13 ¼ 400 Gpa, lI23 ¼ 20 Gpa, mI12 ¼ 0:15; mI13 ¼ 0:26; mI23 ¼ 0:40. The orthotropic orientation relative to the contact orientation chosen using Euler Angles ZXZ conventions is: u ¼ 45 ; h ¼ 60 ; w ¼ 30 . One can note that, at the initial time,the inhomogeneity behaves as a hard inhomogeneity in the plane x2 ¼ 0 and as a soft inhomogeneity in the plane x1 ¼ 0. This is mainly due to the anisotropic material property and the orientation of the inclusion.

4.3.4. Effect of the inhomogeneity location Now, a rigid inhomogeneity located at dx3 ¼ 0:7a is studied. It was shown in Koumi et al. (2014) that the effect of the inhomogeneity on the contact pressure distribution can be neglected if the inclusion center is located at a depth below 0.7 times the contact radius. It can be observed in Fig. 17 that this conclusion is no longer valid in the case of a viscoelastic heterogeneous half-space. The peak magnitude increases from 2:8% (when t ¼ 0) to 6:4% (when t ¼ s). The effect of the inhomogeneity on contact pressure distribution increases with the simulation time t.

38

K.E. Koumi et al. / Mechanics of Materials 77 (2014) 28–42

(a)

(b)

1.2 1

1

t=0

t=0 t=τ /2

0.8

t=τ /2

0.8

0.6

0.6 t=τ

0.4

t=τ

0.4 0.2

0.2 0 −2

1.2

−1.5

−1

−0.5

0

0.5

1

1.5

0 −2

2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Fig. 12. Distribution of normalized contact pressure for c ¼ EI =EM ð0Þ ¼ 2; mI ¼ 0:3 (a1 ¼ 0:4a ; a2 ¼ a3 ¼ 0:1a ; dx3 ¼ 0:4a ; h ¼ 30 ); (a) x1 ¼ 0 and (b) x2 ¼ 0.

(b)

(a) 1.2

1.2 1

1

t=0

t=0 t=τ /2

0.8

t=τ /2

0.8

0.6 0.6

0.2

0.2 0 −2

t=τ

0.4

t=τ

0.4

−1.5

−1

−0.5

0

0.5

1

1.5

2

0 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Fig. 13. Distribution of contact pressure for c ¼ EI =EM ð0Þ ¼ 0:5; mI ¼ 0:3 (a1 ¼ 0:4a ; a2 ¼ a3 ¼ 0:1a ; dx3 ¼ 0:4a ; h ¼ 30 ); (a) x1 ¼ 0 and (b) x2 ¼ 0.

4.3.5. Subsurface stresses Whereas most models found in the literature, as those described in the introduction, are limited to the calculation of the pressure distribution and contact area (for homogeneous bodies), the model presented in this paper allows to obtain the subsurface stresses. In the context of this paper, we will limit ourselves to the normal/shear stress at the inclusion matrix interface and to the von Mises stress. These scalar variables are quite important when developing a decohesion or damage model. In this section, a spherical inhomogeneity with radius r ¼ 0:15a and c ¼ 0:8 centered at (dx1 ¼ 0; dx2 ¼ 0; dx3 =a ¼ 0:4) is considered.  Normal and shear stress at inhomogeneity/matrix interface. The normal and shear stresses at the inhomogeneity/ matrix interface are now investigated (see Fig. 18).

When the point M describes the inclusion matrix interface (see Fig. 3(b)), the angle n ranges from 0 to 360 . Fig. 18(a) and (b) show the normal and shear stresses at the inhomogeneity/matrix interface as a function of the angle n (the position of the point M at this interface). In our case the normal stress is mostly compressive. The peak values normal stress corresponds respectively to 0:83P 0 ; 0:74P 0 ; 0:66P0 when t ¼ 0; t ¼ s=2; t ¼ s. The absolute value of the peak values shear stress at the inclusion matrix interface are equal to 0:28P 0 ; 0:237P0 ; 0:20P0 when t ¼ 0; t ¼ s=2; t ¼ s, respectively.  Distribution of the von Mises stress. The von Mises stress is presented in Fig. 19. One can remark that the position and the value of the peak magnitude of the von Mises stress vary widely with the simulation time. At the t ¼ 0 the maximum value of the von Mises stress corresponds to 0:687P 0 and it located at the

39

K.E. Koumi et al. / Mechanics of Materials 77 (2014) 28–42

(a)

(b)

1.2 1

1

t=0

0.8

t=0

0.8

t=τ /2

t=τ /2

0.6

0.6 t=τ

0.4

t=τ

0.4 0.2

0.2 0 −2

1.2

−1.5

−1

−0.5

0

0.5

1

1.5

0 −2

2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Fig. 14. Distribution of contact pressure for c ¼ EI =EM ð0Þ ¼ 0:2; mI ¼ 0:3 (a1 ¼ 0:4a ; a2 ¼ a3 ¼ 0:1a ; dx3 ¼ 0:4a ; h ¼ 30 ); (a) x1 ¼ 0 and (b) x2 ¼ 0.

(a)

(b)

1.2 1

1

t=0

t=τ /2

0.6

0.6 t=τ

0.4

t=τ

0.4 0.2

0.2 0 −2

t=0

0.8

t=τ /2

0.8

1.2

−1.5

−1

−0.5

0

0.5

1

1.5

2

0 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Fig. 15. Distribution of normalized contact pressure for viscoelastic half-space containing two spherical inhomogeneities with c1 ¼ EI =EM ð0Þ ¼ 0:6 and c2 ¼ EI =EM ð0Þ ¼ 0:8 (with mI ¼ mM ¼ 0:3) on the dimensionless contact pressure distribution for isotropic ellipsoidal inclusion (a1 ¼ 0:4a ; a2 ¼ a3 ¼ 0:1a ; dx3 ¼ 0:3a ); (a) x1 ¼ 0 and (b) x2 ¼ 0.

equator of the inclusion (n ¼ 0 and n ¼ 180 , see Fig. 3(b)) of the inhomogeneity/matrix interface. While when t ¼ s=2; t ¼ s the maximum of von Mises stress corresponds respectively to 0:502P 0 and 0:45P0 ; located at the south pole (n ¼ 270 , see Fig. 3(b)) of the inclusion matrix interface. 5. Conclusion In this paper a numerical method has been proposed to model the effect of an ellipsoidal inhomogeneity embedded into a viscoelastic half- space. The viscoelastic contact problem is fully coupled with the inhomogeneity problem that account for mutual influence between inhomogeneities. The proposed method has been validated by performing a comparison with the results of a finite element model. Contact pressure distribution, subsurface stresses,

normal and shear stresses at the interface between the inclusion and the matrix have been investigated. It can been observed that the effect of the inhomogeneity becomes more pronounced when the simulation time t increases. Two explanations can be given: (i) the substrate material surrounding the inhomogeneity becomes softer when time increases while the inhomogeneity remains elastic (ii) dx3 (the position of the center of the inhomogeneity) is fixed but a ðtÞ (the instantaneous contact radius corresponding to the homogeneous viscoelastic case) increases with time, therefore the relative depth of the inclusion compared to the contact radius dx3 =a ðtÞ decreases. It is shown that the position and the value of the peak magnitude of the von Mises stress vary widely with the simulation time. It can also been observed that the normal and shear stresses at the inhomogeneity/matrix interface decrease when the simulation time t increases.

40

K.E. Koumi et al. / Mechanics of Materials 77 (2014) 28–42

(a)

(b)

1.2 1

1.2 1

t=0 t=τ /2

0.8

0.6 t=τ

0.4

0.2

−1.5

−1

−0.5

0

0.5

1

1.5

0 −2

2

Fig. 16. Dimensionless contact pressure distribution in the EI2 ¼ 623 Gpa, (a1 ¼ 0:4a ; a2 ¼ a3 ¼ 0:1a ; dx3 ¼ 0:4a ; h ¼ 30 ; EI1 ¼ 210 Gpa, I I I m12 ¼ 0:15; m13 ¼ 0:26; m23 ¼ 0:40); (a) x1 ¼ 0 and (b) x2 ¼ 0.

1

−0.5

an

0

ellipsoidal

lI12 ¼ 83 Gpa,

0.5

1

1.5

anisotropic

lI13 ¼ 400 Gpa,

1

t=τ /2

2

inhomogeneity lI23 ¼ 20 Gpa,

t=0 t=τ /2

0.8

0.6

0.6 t=τ

0.4

t=τ

0.4

0.2

0.2

−1.5

−1

−0.5

0

0.5

1

1.5

0 −2

2

−1.5

Fig. 17. Dimensionless contact pressure distribution in the case of an ellipsoidal (a1 ¼ 0:4a ; a2 ¼ a3 ¼ 0:1a ; h ¼ 30 ; c ¼ EI =EM ð0Þ ! 1; mI ¼ 0:3); (a) x1 ¼ 0 and (b) x2 ¼ 0.

(a)

−1

1.2

t=0

0.8

−1.5

case of EI3 ¼ 50 Gpa,

(b)

1.2

0 −2

t=τ

0.4

0.2

(a)

t=τ /2

0.8

0.6

0 −2

t=0

0.25

(b)

0

−1

−0.5

hard

0

0.5

inhomogeneity

1

1.5

located

2

at

dx3 =a ¼ 0:7

0.5

0.25

−0.25 0 −0.5 −0.25

−0.75

−1 0

60

120

180

240

300

Fig. 18. Dimensionless normal and shear stresses (rn =P 0 and ða1 ¼ a2 ¼ a3 ¼ 0:15a ; dx3 =a ¼ 0:4Þ; (a) rn =P 0 and (b) s=P 0 .

360

−0.5 0

60

120

180

240

300

360

s=P0 ) in the plane x2 ¼ 0 in the presence of a single isotropic and spherical inclusion

41

K.E. Koumi et al. / Mechanics of Materials 77 (2014) 28–42

Appendix B. Stresses in a half-space subject to normal pressure ðM ij Þ An isotropic half-space is submitted an uniform normal pressure rn in a discretized surface area of 2Dx1  2Dx2 at the center point Pðx01 ; x02 ; 0Þ. The stress at an observation point Q ðx1 ; x2 ; x3 Þ is given by Zhou et al. (2009) and Johnson (1985):

rij ðx1 ; x2 ; x3 Þ ¼ Mij ðx1  x01 ; x2  x02 ; x3 Þrn ðx1 ; x2 Þ rij ðx1 ; x2 ; x3 Þ ¼

rn ½hij ðn1 þ Dx1 ; n2 þ Dx2 ; n3 Þ 2p  hij ðn1 þ Dx1 ; n2  Dx2 ; n3 Þ þ hij ðn1  Dx1 ; n2  Dx2 ; n3 Þ  hij ðn1  Dx1 ; n2 þ Dx2 ; n3 Þ

where

ni ¼ xi  x0i The functions hij ðÞ above are defined by

x22 þ x23  qx2 þ 2ð1  mÞ tan1 x1 x3 q  x2 þ x3 x1 x2 x3  þ x1 qðx21 þ x23 Þ

h11 ðx1 ; x2 ; x3 Þ ¼ 2m tan1

h22 ðx1 ; x2 ; x3 Þ ¼ h11 ðx2 ; x1 ; x3 Þ x22 þ x23  qx2 x1 x3   x1 x2 x3 1 1  þ 2 2 2 2 q x1 þ x3 x2 þ x3

h33 ðx1 ; x2 ; x3 Þ ¼ tan1 Fig. 19. Dimensionless von Mises stress (rVM =P 0 ) in the plane x2 ¼ 0 in the presence of a single isotropic and spherical inclusion ða1 ¼ a2 ¼ a3 ¼ 0:15a ; dx3 =a ¼ 0:4Þ; (a) t ¼ 0, (b) t ¼ s=2 and (c) t ¼ s.

Appendix A. Stress in a half-space due to a concentrated unit normal force at the surface origin (F ij )

h12 ðx1 ; x2 ; x3 Þ ¼ 

h13 ðx1 ; x2 ; x3 Þ ¼  F 11 ðx1 ; x2 ; x3 Þ ¼

    1 1  2m x3 x21  x22 x3 x22 3x3 x21 1 þ 3  2 2 5 2p r q r q q

x3

q

 ð1  2mÞ logðq þ x3 Þ

x2 x23 qðx31 þ x23 Þ

h23 ðx1 ; x2 ; x3 Þ ¼ h13 ðx2 ; x1 ; x3 Þ where

F 22 ðx1 ; x2 ; x3 Þ ¼ F 11 ðx2 ; x1 ; x3 Þ F 33 ðx1 ; x2 ; x3 Þ ¼ 

F 12 ðx1 ; x2 ; x3 Þ ¼

3 x33 2p q5

    1 1  2m x3 x1 x2 x3 x2 x1 3x3 x2 x1 1  þ  2p r2 q r2 q3 q5

F 13 ðx1 ; x2 ; x3 Þ ¼ 

3 x1 x23 2p q5

F 23 ðx1 ; x2 ; x3 Þ ¼ F 12 ðx2 ; x1 ; x3 Þ where

r 2 ¼ x21 þ x22 ; q ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x21 þ x23 þ x23

with m, the Poisson’s ratio of the isotropic half-space.

q¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x21 þ x22 þ x23

Appendix C. Normal displacement at the surface subject to normal pressure ðK n Þ The contact between a sphere and an elastic half-space having respectively elastic constants ðE1 ; m1 Þ and ðE2 ; m2 Þ, where the surface x3 ¼ 0 is discretized into rectangular surface area of 2D1  2D2 , is now considered. The initial contact point coincides with the origin of the Cartesian coordinate system ðx1 ; x2 ; x3 Þ. The relationship between the normal displacement at an observation point Pðn1 ; n2 ; 0Þ and the pressure field at the center Q ðn01 ; n02 ; 0Þ is built using the function K n .

K n ðc1 ; c2 Þ ¼

  4 1  m21 1  m22 X þ K n ðc ; c Þ pE1 pE2 p¼1 p 1 2

42

K.E. Koumi et al. / Mechanics of Materials 77 (2014) 28–42

K n1 ðc1 ; c2 Þ ¼ ðc1 þ D1 Þ 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 ðc þ D Þ þ ðc2 þ D2 Þ2 þ ðc1 þ D1 Þ2 C 2 B 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA  log @ ðc2  D2 Þ þ ðc2  D2 Þ2 þ ðc1 þ D1 Þ2 K n2 ðc1 ; c2 Þ ¼ ðc2 þ D2 Þ 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 2 2 Bðc1 þ D1 Þ þ ðc2 þ D2 Þ þ ðc1 þ D1 Þ C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA  log @ ðc1  D1 Þ þ ðc2 þ D2 Þ2 þ ðc1  D1 Þ2 K n3 ðc1 ; c2 Þ ¼ ðc1  D1 Þ 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 ðc  D Þ þ ðc2  D2 Þ2 þ ðc1  D1 Þ2 C 2 2 B qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA  log @ ðc2 þ D2 Þ þ ðc2 þ D2 Þ2 þ ðc1  D1 Þ2 K n4 ðc1 ; c2 Þ ¼ ðc2  D2 Þ 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 ðc  D Þ þ ðc2  D2 Þ2 þ ðc1  D1 Þ2 C 1 B 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiA  log @ ðc1 þ D1 Þ þ ðc2  D2 Þ2 þ ðc1 þ D1 Þ2 where

c1 ¼ n1  n01

and c2 ¼ n2  n02

References Argatov, I., 2012. An analytical solution of the rebound indentation problem for an isotropic linear viscoelastic layer loaded with a spherical punch. Acta Mech. 223, 1441–1453. Carbone, G., Putignano, C., 2013. A novel methodology to predict sliding/ rolling friction in viscoelastic materials: theory and experiments. J. Mech. Phys. Solid. 61, 1822–1834. Chen, W., Wang, Q., Huan, Z., Luo, X., 2011. Semi-analytical viscoelastic contact modeling of polymer-based materials. J. Tribol. 133, 041404. Chiu, Y., 1978. On the stress field and surface deformation in a half space with a cuboidal zone in which initial strains are uniform. J. Appl. Mech. 45, 302–306. Courbon, J., Lormand, G., Dudragne, G., Daguier, P., Vincent, A., 2005. Influence of inclusion pairs, clusters and stringers on the lower bound of the endurance limit of bearing steels. Tribol. Int. 36, 921–928. Goriacheva, G., 1973. Contact problem of rolling of a viscoelastic cylinder on a base of the same material. J. Appl. Math. Mech. 37 (5), 925–933. Graham, G., 1967. The contact problem in the linear theory of viscoelasticity when the time dependent contact area has any number of maxima and minima. Int. J. Eng. Sci. 5, 495–514. Greenwood, J., 2010. Contact between an axisymmetric indenter and a viscoelastic half-space. Int. J. Mech. Sci. 52, 829–835. Hunter, S., 1960. The hertz problem for a rigid spherical indenter and viscoelastic half-space. J. Mech. Phys. Solids 8, 219–234. Hunter, S., 1961. The rolling contact of a rigid cylinder with a viscoelastic half space. J. Appl. Mech. 28, 611–617. Jacq, C., Nelias, D., Lormand, G., Girodin, D., 2002. Development of a threedimensional semi-analytical elastic-plastic contact code. J. Tribol. 124 (4), 653–667.

Johnson, K., 1985. Contact Mechanics. Cambridge University Press. Kabo, E., Ekberg, A., 2002. Fatigue initiation in railway wheels – a numerical study of the influence of defects. Wear 253, 26–34. Kabo, E., Ekberg, A., 2005. Material defects in rolling contact fatigue of railway wheels-the influence of defect size. Wear 258, 1194–1200. Koumi, K., Zhao, L., Leroux, J., Chaise, T., Nelias, D., 2014. Contact analysis in the presence of an ellipsoidal inhomogeneity within a half space. Int. J. Solids Struct. 51, 1390–1402. Kuo, C., 2007. Stress disturbances caused by the inhomogeneity in an elastic half-plane subjected to contact loading. Int. J. Solids Struct. 44, 860–873. Kuo, C., 2008. Contact stress analysis of an elastic half-plane containing multiple inclusions. Int. J. Solids Struct. 45, 4562–4573. Lee, E., Radok, J., 1960. The contact problem for viscoelastic bodies. J. Appl. Mech. 27, 438–444. Leroux, J., Nelias, D., 2011. Stick-slip analysis of a circular point contact between a rigid sphere and a flat unidirectional composite with cylindrical fibers. Int. J. Solids Struct. 48, 3510–3520. Leroux, J., Fulleringer, B., Nelias, D., 2010. Contact analysis in presence of spherical inhomogeneities within a half-space. Int. J. Solids Struct. 47, 3034–3049. Le Tallec, P., Rahler, C., 1994. Numerical models of steady rolling for nonlinear viscoelastic structures in finite deformations. Int. J. Numer. Methods Eng. 37, 1159–1186. Liu, S., Wang, Q., Liu, G., 2000. A versatile method of discrete convolution and fft (dc-fft) for contact analyses. Wear 243 (1–2), 101–111. Love, A., 1952. A Treatise on the Mathematical Theory of Elasticity, fourth ed. Cambridge University Press. Mura, T., 1987. Micromechanics of Defects in Solids, second ed. Kluwer Academic Publishers. Nackenhorst, U., 2004. The ale-formulation of bodies in rolling contact: theoretical foundations and finite element approach. Comput. Methods Appl. M 193, 4299–4322. Nasdala, L., Kaliske, M., Becker, A., Rothert, H., 1998. An efficient viscoelastic formulation for steady-state rolling structures. Comput. Mech. 22, 395–403. Padovan, J., Paramadilok, O., 1984. Transient and steady state viscoelastic rolling contact. Comput. Struct. 20, 545–553. Padovan, J., Kazempour, A., Tabaddor, F., Brockman, B., 1992. Alternative formulations of rolling contact problems. Finite Elem. Anal. Design 11, 275–284. Ting, T., 1966. The contact stresses between a rigid indenter and a viscoelastic half-space. J. Appl. Mech. 33, 845–854. Ting, T., 1968. Contact problems in the linear theory of viscoelasticity. J. Appl. Mech. 35, 248–254. Vollebregt, E., 2011. User guide for contact. J.J. Kalker’s variational contact model. Tech. rep., Technical Report TR09-03. Zhou, K., Chen, W., Keer, L., Wang, Q., 2009. A fast method for solving three-dimensional arbitrarily shaped inclusions in a half space. Comput. Methods Appl. Mech. Eng. 198 (9–12), 885–892. Zhou, K., Chen, W., Keer, L., Ai, X., Sawamiphakdi, K., Glaws, P., Wang, Q., 2011a. Multiple 3d inhomogeneous inclusions in a half space under contact loading. Mech. Mater. 43, 444–457. Zhou, K., Keer, L., Wang, Q., 2011b. Semi-analytic solution for multiple interacting three-dimensional inhomogeneous inclusions of arbitrary shape in an infinite space. Int. J. Numer. Methods Eng. 87, 617–638. Zhou, K., Keer, L., Wang, Q.J., 2011c. Semi-analytic solution for multiple interacting three-dimensional inhomogeneous inclusions of arbitrary shape in an infinite space. Int. J. Numer. Methods Eng. 87, 617–638. Zhou, K., Keer, L., Wang, Q., Ai, X., Sawamiphakdi, K., Glaws, P., Paire, M., Che, F., 2012. Interaction of multiple inhomogeneous inclusions beneath a surface. Comput. Methods Appl. Mech. Eng. 217–220, 25– 33. Zhou, K., Hoh, H.J., Wang, X., Keer, L.M., Pang, J.H.L., Song, B., Wang, Q.J., 2013. A review of recent works on inclusions. Mech. Mater. 60, 144– 158.