Contact analysis in the presence of an ellipsoidal inhomogeneity within a half space

Accepted Manuscript Contact analysis in the presence of an ellipsoidal inhomogeneity within a half space Koffi Espoir Koumi, Lv Zhao, Julien Leroux, Thibaut Chaise, Daniel Nelias PII: DOI: Reference:

S0020-7683(13)00516-7 http://dx.doi.org/10.1016/j.ijsolstr.2013.12.035 SAS 8238

To appear in:

International Journal of Solids and Structures

Received Date: Revised Date:

30 September 2013 25 November 2013

Please cite this article as: Koumi, K.E., Zhao, L., Leroux, J., Chaise, T., Nelias, D., Contact analysis in the presence of an ellipsoidal inhomogeneity within a half space, International Journal of Solids and Structures (2013), doi: http://dx.doi.org/10.1016/j.ijsolstr.2013.12.035

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

1

2

3

4 5

6

Contact analysis in the presence of an ellipsoidal inhomogeneity within a half space Koffi Espoir KOUMIa,b , Lv ZHAOa , Julien LEROUXb , Thibaut CHAISEa , Daniel NELIAS1a a Universit´ e

7

8

de Lyon, INSA-Lyon, LaMCoS UMR CNRS 5259, F69621 Villeurbanne, FRANCE b SNECMA, Centre de Villaroche, 77550 Moissy Cramayel, FRANCE

Abstract Many materials contain inhomogeneities or inclusions that may greatly affect their mechanical properties. Such inhomogeneities are for example encountered in the case of composite materials or materials containing precipitates. This paper presents an analysis of contact pressure and subsurface stress field for contact problems in the presence of anisotropic elastic inhomogeneities of ellipsoidal shape. Accounting for any orientation and material properties of the inhomogeneities are the major novelties of this work. The semi-analytical method proposed to solve the contact problem is based on Eshelby’s formalism and uses 2D and 3D Fast Fourier Transforms to speed up the computation. The time and memory necessary are greatly reduced in comparison with the classical finite element method. The model can be seen as an enrichment technique where the enrichment fields from the heterogeneous solution are superimposed to the homogeneous problem. The definition of complex geometries made by combination of inclusions can easily be achieved. A parametric analysis on the effect of elastic properties and geometrical features of the inhomogeneity (size, depth and orientation) 1 corresponding

author, [email protected]

1

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. It is shown that the presence of an inclusion below the contact surface affects significantly the contact pressure and subsurfaces stress distributions when located at a depth lower than 0.7 times the contact radius. The anisotropy directions and material data are also key elements that strongly affect the elastic contact solution. In the case of normal contact between a spherical indenter and an elastic half space containing a single inhomogeneity whose center is located straight below the contact center, the normal stress at the inhomogeneity/matrix interface is mostly compressive. Finally when the axes of the ellipsoidal inclusion do not coincide with the contact problem axes, the pressure distribution is not symmetrical. 9

Keywords: Inhomogeneity, Inclusion, Anisotropy, Eigenstrain, Eshelby’s Equivalent

10

Inclusion Method (EIM), Contact Mechanics

11

1. Introduction

12

From a fundamental point of view the presence of inhomogeneities in a matrix is

13

very interesting for contact problems, first because it modifies the subsurface stress

14

field, and also because it disturbes very significantly the contact pressure distribution.

15

On one hand, for metallic materials, aluminum alloys or ceramics, failure mechanisms

16

can be initiated by the presence of material impurities such as cavities, inclusions or

17

precipitates. Inclusions can raise internal stresses, at the origin of fatigue cracks which

18

in turn reduce the service life of mechanical components such as rolling elements ∗ [email protected]

Preprint submitted to IJSS

January 2, 2014

19

(Voskamp, 1985; Nelias et al., 1999; Nugent et al., 2000; Vignal et al., 2003; Chen

20

et al., 2008a). On the other hand, in heterogeneous materials such as composite mate-

21

rials, the presence of inhomogeneities that are added to improve structural properties

22

strongly affects the distribution of stresses including at the surface interface.

23

The disturbance field caused by the presence of an inhomogeneity embedded in

24

an infinite space has been investigated by many authors [Eshelby (1957, 1959, 1961),

25

Willis (1964), Walpole (1967), Asaro and Barnett (1975), Mura and Furuhashi (1984)].

26

Eshelby (1957, 1959) was the first to initiate a method dealing with an inhomogeneity

27

in an infinite space, which is called the ’Equivalent Inclusion Method’ (EIM). Moscho-

28

vidis and Mura (1975) extended this method to obtain the Eshelby’s tensor for the

29

exterior points of an ellipsoidal inclusion as a function of harmonic and bi-harmonic

30

potentials. All these authors limited their analysis to an infinite body.

31

The analytical solution for displacements, strains and stresses in a half-space are

32

difficult to obtain. Many authors investigated the effect of inhomogeneities in a half-

33

space subjected to a prescribed load. To solve the problem, some assumptions are

34

usually made. Mindlin and Cheng (1950) considered a half-space containing a hy-

35

drostatic eigenstrain; Aderogba (1976) assumed a spherical inclusion with an arbitrary

36

eigenstrain; Chiu (1977, 1978) limited his analysis to a cuboidal inclusion with an in-

37

compressible eigenstrain; Seo and Mura (1970) assumed an ellipsoidal inclusion with

38

a pure dilatational eigenstrain.

39

Meanwhile, for the exterior points of the ellipsoidal inclusion in an infinite space,

40

the Eshelby’s tensor was expressed in terms of the harmonic and bi-harmonic potentials

41

of the inclusion (Mura, 1987). The potentials can also be written in terms of elliptical

3

42

integrals.

43

Several authors already investigated stress concentration due to the presence of in-

44

clusions or inhomogeneities. However until very recently they did not solve the contact

45

problem, instead they assumed the hertzian contact pressure distribution as input, see

46

for example Kabo and Ekberg (2002, 2005) or Courbon et al. (2005). For the 2D con-

47

tact problem the effect of the presence of an elliptical heterogeneous inclusion on the

48

contact loading was first solved numerically by Kuo (2007, 2008) using the Boundary

49

Element Method (BEM). The effect of inhomogeneities on the three-dimensional con-

50

tact problem solution was first investigated by Nelias and co-workers (Leroux et al.,

51

2010; Leroux and Nelias, 2011) and then by Zhou et al. (2011a) based on the semi-

52

analytical method (SAM) initially proposed by Jacq et al. (2002) to numerically solve

53

3D elastic-plastic contact. For more details on the modeling of inclusions in contact

54

problems the reader may refer to the review paper by Zhou et al. (2013).

55

SAMs have been continuously developed since ten years, and then successfully ap-

56

plied to several leading edge problems such as thermo-elastic-plastic contact modeling

57

(Boucly et al., 2005), modeling of plasticity and accumulation of plastic strains (Boucly

58

et al., 2005; Wang and Keer, 2005), running-in (Nelias et al., 2007b) and wear modeling

59

(Gallego et al., 2006; Gallego and Nelias, 2007; Gallego et al., 2010b,a), simulation of

60

single impact (Chaise et al., 2011), shot peening (Chaise et al., 2012) and low plasticity

61

burnishing (Nelias et al., 2007a; Chen et al., 2008b; Chaise and Nelias, 2011), model-

62

ing of cuboidal inclusions (Liu and Wang, 2005; Zhou et al., 2009, 2011a,a,b, 2012), as

63

well as to account for material or coating anisotropy (Bagault et al., 2012, 2013). With

64

the same technique Wang et al. (2009) solved a contact between two joined quarter

4

65

spaces and a rigid sphere.

66

In this paper, the analysis will be limited to a single elastic ellipsoidal inhomogene-

67

ity which can be anisotropic and with any orientation. This is one of the key steps for

68

the modeling of contact for composite materials since fibers are mostly anisotropic.

69

It should be outlined that the intent here is not to simulate the macroscopic response

70

of the contact assuming homogenized and anisotropic material properties, as very re-

71

cently studied by Rodriguez-Tembleque et al. (2013). The effects of Young’s modulus,

72

bulk coefficient and shear modulus will be first investigated. It will be shown that

73

the presence of an heterogeneous inclusion will not only affect the local stress field,

74

but also very significantly the contact pressure distribution. Finally the presence of an

75

anisotropic and ellipsoidal inclusion located below the surface and with its axes not

76

coincident with the contact ones will be presented.

77

2. Contact problem formulation

78

79

Generally, the formulation of the normal contact between two finite bodies B1 and B2 (Fig. 1) consists in a set of equations and inequalities that are recalled below:

• The load balance. The applied load W and the integration of the contact pressure p(x1 , x2 ) in the contact region Γc must be strictly equal. W=

R Γc

p(x1 , x2 )dΓ

(1)

• The surface separation. The gap between the two contacting surfaces is: 1 +B2 ) h(x1 , x2 ) = hi (x1 , x2 ) + δ + u(B (x1 , x2 ) 3

5

(2)

80

1 +B2 ) where hi (x1 , x2 ) is the initial geometry, δ the rigid body displacement, and u(B (x1 , x2 ) 3

81

the sum of normal displacements of surfaces 1 and 2, that can be due to elastic

82

deflection (under loading only), plastic deformation or the presence of inhomo-

83

geneities. • The contact conditions. The distance h(x1 , x2 ) is always positive, because the

contacting bodies can not interpenetrate each other. The conditions are defined by the inequalities: h(x1 , x2 ) > 0

contact : hi (x1 , x2 ) = 0 and p(x1 , x2 ) > 0 separation : hi (x1 , x2 ) > 0 and p(x1 , x2 ) = 0

W M2

O

 x1

hi

 x2 M1

 x3 W

Figure 1: Contact problem description.

6

(3)

84

85

3. Theoretical Background - Eshelby’s equivalent inclusion method in Contact Mechanics

86

Given that the presence of inhomogeneities creates an incompatibility of deforma-

87

tion between the inhomogeneities and the matrix, the Eshelby’s equivalent inclusion

88

method is used.

89

3.1. Eshelby’s solution for an infinite space

90

An infinite matrix D with the elastic stiffness tensor CiMjkl containing an ellipsoidal

91

domain Ω with the elastic stiffness tensor CiIjkl is submitted at infinity to a uniform strain

92

ε0 . The strain field is disturbed by the presence of the inhomogeneity. The disturbed

93

contributions of the stresses and displacements are denoted σi j and ui respectively. The

94

total stress is σi j + σ0i j and the total displacement ui + u0i .

95

The stress components are in self-equilibrium if one neglects the body forces:

σi j, j = 0

96

(4)

and σi j = 0 at infinity. Under the hypothesis of linear isotropic elasticity, the stress components are calculated by the Hooke’s law:

σ0i j + σi j = CiIjkl (u0k,l + uk,l ) = CiIjkl (ε0kl + εkl ) σ0i j + σi j = CiMjkl (u0k,l + uk,l ) = CiMjkl (ε0kl + εkl )

in Ω in D − Ω

(5)

97

The Eshelby’s equivalent inclusion method (EIM) consists in representing the el-

98

lipsoidal inhomogeneity as an inclusion having the same elastic propreties CiMjkl as the 7

99

100

matrix but being subjected to an additional imaginary strain called eigenstrain ε∗ giving:

CiIjkl (ε0kl + εkl ) = CiMjkl (ε0kl + εkl − ε∗kl )

in Ω

(6)

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.(6). In particular, the eigenstrain ε∗i j is related to compatibility strain εi j by: εi j = S i jkl × ε∗kl ,

101

(7)

where S i jkl is the Eshelby’s tensor. Substitution of Eq.(7) into Eq.(6) leads to: ∆Ci jkl S klmn ε∗mn + CiMjkl ε∗kl = −∆Ci jkl ε0kl

(8)

where ∆Ci jkl = CiIjkl − CiMjkl

When the inhomogeneity contains an initial eigenstrain ε p (inhomogeneous inclusion), previous equations (Eq.6 - Eq.8) are modified as follows.

σ0i j + σi j = CiIjkl (ε0kl + εkl − εklp ) σ0i j + σi j = CiMjkl (ε0kl + εkl )

in Ω

in D − Ω

(9)

Using the EIM it yields, CiIjkl (ε0kl + εkl − ε p ) = CiMjkl (ε0kl + εkl − εklp − ε∗kl )

8

(10)

Hence, p M ∗∗ ∗∗ σi j = CiIjkl (S klmn ε∗∗ mn − ε ) = C i jkl (S klmn εmn − εkl )

(11)

where ε∗∗ = ε∗ + ε p

102

3.2. Determination of compatibility strain

103

The Eshelby’s solution is only valid for a uniformly applied strain. However, the

104

contact problem involves nonuniform strains. Moschovidis and Mura (1975) have ex-

105

tended the Eshelby’s EIM for nonuniformly applied strain. When the applied strain is a

106

n-order polynomial, the equivalent eigenstrain is also treated as a n-order polynomial. Considering an ellipsoidal inhomogeneity, if the applied strain is given by: ε0i j (x) = Ei0j + Ei0jk xk + Ei0jkl xk xl + · · · .

(12)

where the coefficients E 0 are constant, then the eigenstrain is assumed to be: ε∗i j (x) = Bi j + Bi jk xk + Bi jkl xk xl + · · · .

(13)

The strain associated to the eigenstrain is given as: εi j (x) = Di jkl (x)Bkl + Di jklq (x)Bklq + Di jklqr (x)Bklqr + · · · .

(14)

107

For the interior points of an ellipsoidal inhomogeneity, the coefficients Di jkl are constant

108

and Di jklq are linear in x. The expressions for Di jkl , Di jklq , Di jklqr are given in Mura (1987).

109

The results presented here are limited to a uniform eigenstrain, therefore, only the

110

calculation of the tensor Di jkl is performed.

9

Di jkl =

1 [Ψ,i jkl − 2νδkl φ,i j − (1 − ν)(δkl φil + δki φ, jl + δ jl φ,ik + δli φ, jk )] 8π(1 − ν) Z Ψ (x) = |x − x0 |dx0

(15)

Ω

Z

φ(x) =

Ω

1 dx0 |x − x0 |

The harmonic potential φ(x) and the biharmonic potential Ψ (x) can be expressed as a function of the elliptical integrals E(θ0 , k) and F(θ0 , k) (Gradshteyn and Ryzhik, 1965), where: E(θ , k) =

θ0

Z

0

(1 − k2 sinw)1/2 dw 0

F(θ0 , k) =

θ0

Z 0

1 (1 −

k2 sinw)1/2

θ0 = sin−1 1 − k=

a23

dw

(16)

!1/2

a21

3(a21 − a22 ) (a21 − a23 )

(17)

Assuming that a1 > a2 > a3 , with a1 ,a2 ,a3 the semi-axes of the ellipsoidal inclusion. The Eshelby’s tensor S i jkl is obtained from Eq. (15) as: S i jkl = Di jkl (xI )

111

where xI = (x1I , x2I , x3I ) represents the cartesian coordinates of the inclusion center.

112

3.3. Half-space solution

(18)

113

Three dimensional contact problems involve a half-space that is bounded by the

114

surface plane x3 = 0 in the cartesian coordinate system ( x1 , x2 , x3 ). Mura (1987) intro-

115

duced a method to determine the eigenstrain. The obtained equations, despite of their

116

complexity, are valid only in the case of hydrostatic eigenstrains. In order to get rid of

10

117

this restrictive hypothesis which is incompatible with the contact problem, Zhou et al.

118

(2009) proposed an ingenious method allowing to extend the previous solution, valid

119

only for infinite spaces, to half spaces. The solution for an isotropic half space con-

120

sists in decomposing the problem into three subproblems (Fig. 2), known as Chiu’s

121

decomposition (Chiu, 1978).

122

123

124

125

126

127

(1) An inclusion with the prescribed eigenstrain ε∗ = (ε∗11 ; ε∗22 ; ε∗33 ; ε∗12 ; ε∗13 ; ε∗23 ) in an infinite space.

(2) A symmetric inclusion with a mirror eigenstrain ε∗s = (ε∗11 ; ε∗22 ; ε∗33 ; ε∗12 ; −ε∗13 ; −ε∗23 ) in the same space.

(3) A normal traction distribution −σn at the surface of the half space ( x3 = 0) which is a function of the eigenstrains ε∗ and ε∗s .

𝜀∗ 𝜎

O

𝑥 ,𝑥

=

O

𝑥 ,𝑥

+

O

𝑥 ,𝑥

-

𝑥 ,𝑥

O

𝜀∗ 𝑥 ,𝑥

𝑥 ,𝑥

𝑥 ,𝑥

𝑥 ,𝑥

Figure 2: Decomposition of the half-space solution into three sub-problems.

128

The summation of the two solutions (1) and (2) leaves the plane of symmetry ( x3 =

129

0) free of shear tractions. By adding an opposite normal stress σn , the condition of

130

free surface traction is obtained. The stress at any point of the domain meshed with

11

131

n1 × n2 × n3 cuboids is given by:

σi j (x1 , x2 , x3 ) =

nX 3 −1 n 2 −1 n 1 −1 X X

Bi jkl (x1 − x1I , x2 − x2I , x3 − x3I )ε∗kl (x1I , x2I , x3I )

x3I =0 x2I =0 x1I =0

+

nX 3 −1 n 2 −1 n 1 −1 X X

Bi jkl (x1 − x1I , x2 − x2I , x3 + x3I )ε∗skl (x1I , x2I , −x3I )

(19)

x3I =0 x2I =0 x1I =0

−

nX 2 −1 n 1 −1 X

Mi j (x1 − x1I , x2 − x2I , x3 )σn (x1I , x2I , 0)

x2I =0 x1I =0 132

where Bi jkl are the influence coefficients that relate the constant eigenstrain at the

133

point (x1I , x2I , x3I ) which is the inclusion center in an infinite space to the stress σi j at the

134

point (x1 , x2 , x3 ). Mi j represent the influence coefficients relating the normal traction σn

135

within a discretized area centered at (x1I , x2I , 0) to the stress σi j at the point (x1 , x2 , x3 ).

Bi jkl (x) = CiMjmn Dmnkl (x)

for x in D − Ω

Bi jkl (x) = CiMjmn (Dmnkl (x) − Imnkl ) 136

137

138

for x in Ω

(20) (21)

where Ii jkl = 12 (δil δ jk + δik δ jl) is the fourth-order identity tensor. For a single inclusion centered at (x1I , x2I , x3I ) in the half-space, the normal traction σn at the surface point (x1, , x2, , 0) is obtained as:

σn (x10 , x20 , 0) = − B33kl (x10 − x1I , x20 − x2I , −x3I , )ε∗kl (x1I , x2I , x3I )

(22) −

B33kl (x10

−

x1I , x20

−

x2I , x3I , )ε∗skl (x1I , x2I , −x3I )

139

In Eq.(19), each component Mi j () is obtained by a double integration of the function

140

Fi j () over a discretized surface area 2∆x1 × 2∆x2 centered at (x1I , x2I , 0), see appendices A

141

and B.

12

Mi j (x1 − x1I , x2 − x2I , x3 ) =

Z

x1I +∆x1 x1I −∆x1

x2I +∆x2

Z

x2I −∆x2

Fi j (x1 − x10 , x2 − x20 , x3 )dx10 x20

(23)

142

The 3D-FFT is used to accelerate the calculation of the first (1) and second terms (2)

143

and the 2D-FFT for the third term (3). Wrap around order and zero-padding techniques

144

are used in order to remove the induced periodicity error (Liu et al., 2000).

145

3.4. Normal displacement of a surface point

146

The surface normal ’eigen-displacements’ can be obtained when inserting the eigen-

147

strain into the total strain. They are generated by the pressure field σn only. The normal

148

displacements are calculated as:

u3 (x1 , x2 ) =

nX 2 −1 n 1 −1 X

K n (x1 − x10 , x2 − x20 )σn (x10 , x20 )

(24)

x20 =0 x10 =0 149

To solve the equation above numerically, the surface in contact is discretized into

150

n1 × n2 rectangular elements of uniform size 2∆x1 × 2∆x2 . Then, pressure and displace-

151

ment within each discrete patch are treated as constant and their values located at the

152

center. The effect of an uniform pressure on a rectangular area has been given by Love

153

(1952) and Johnson (1985). K n denotes the influence coefficients that relate the normal

154

pressure at the surface point (x10 , x20 , 0) to the normal displacement at the surface point

155

(x1 , x2 , 0), recalled in Appendix C.

156

4. Integration of the inhomogeneity effects in the Contact Algorithm

157

In order to integrate the inhomogeneity effects in the contact algorithm, an equiva-

158

lent elastic algorithm is proposed, as shown in Fig. 3. The left frame in red presents the

13

Normal problem Pressure

Tangential problem Shear Equivalent inclusion model

Subsurface stresses (3D-FFT)

Equivalent elastic problem

For two mirror-image eigenstrains

Determination of eigenstrain ε∗ calculation Disturbance stresses

Stresses on extended surface

Misfit displacements induced by inclusions

Misfit displacements and stresses induced by the pressure field Ux,Uy,Uz Half-space subject to pressure

Sum of stresses and displacements Convergence on displacements

Figure 3: Flowchart of the algorithm for the equivalent elastic problem.

159

calculation of displacements due to the eigenstrain and the right frame in blue the ones

160

due to the contact pressure. The effect of an inclusion on the contact problem derives

161

from the fact that the surface contact geometry is modified by the eigen-displacement

162

produces by the eigenstrain. The contact pressure and shears are then updated, which

163

modifies the eigenstrain value. The elastic displacements are obtained from the up-

164

dated contact pressure via the resolution of the elastic contact problem. The algorithm

165

is repeated until convergence of the normal displacements is obtained.

166

It should be noted that in the frictional contact between dissimilar elastic materi-

167

als, the tangential displacements of the surface points, as analyzed by Fulleringer and

168

Nelias (2010) in analytical form for a cuboid of uniform eigenstrain, should be consid-

169

ered.

14

170

171

172

5. How to consider the orientation of an ellipsoidal inclusion In order to take into account the orientation of the inclusion (Fig. 4), it is necessary to respect the three rules described below.

173

• Chiu’s decomposition: Given that the orientations of the source inclusion and

174

the mirror inclusion are different, the influence coefficients of the two inclusions

175

are consequently different. However, the Chiu’s symmetry permits to offset the

176

effect of the shear stresses on the free surface.

177

• The Wrap Around: In the Wrap Around technique, the extension of the zone

178

should be considered. Based on the symmetries, this technique is only valid in

179

the inclusion’s coordinate system, which implies a development of new coeffi-

180

cients of symmetries. On each mesh point of the contact coordinate system, it is

181

necessary to associate a corresponding image point in the inclusion coordinate

182

system and the corresponding coefficient of symmetry.

183

• The third step consists in taking into account the correct change of coordinate

184

both for the source inclusion and the mirror inclusion. Hence, two transition

185

matrixes based on the ZXZ convention of Euler angles are introduced, one for

186

the source inclusion and the other one for the mirror inclusion.

187

188

It should be noted that this model is valid whatever the contacting surface geometry is, in other words it is not limited to the contact of a sphere and a flat.

15

A-A

P

P 𝑅

o

X2

A

X1

𝑎1 𝑎2

2𝑎 ∗

o

A

dx3 M

𝑎3

𝜉

Ɵ 2𝑎1

2𝑎3

X3

x1

x3

(a)

(b)

Figure 4: Contact of a sphere over a flat in the presence of an arbitrarily oriented inclusion. (a) 3D configuration and (b) 2D configuration.

189

190

191

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

192

The configuration is shown in Fig. 5. A contact between a sphere and a plane (solid)

193

containing an ellipsoidal inclusion at a certain depth (dx3 ) is described. In order to

194

correctly describe the interface, the inclusion is obtained by making a solid partition in

195

a local coordinate system related to the contact coordinate, via a rotation of 30 degrees

196

around the x2 axis (Fig. 6). The size of the sphere and the solid are recalled in Table 1.

197

The geometrical properties and position of the inclusion are normalized by the contact

198

half-width a∗ , as shown in Table 2.

199

Since the model is symmetric with respect to the plane x2 = 0, only half of the

200

bodies in contact are meshed. The solid is clamped on the lower surface, the sphere

201

can just move vertically. A point force is applied on the sphere generating a half-width 16

Sphere

Subsurface layer

Solid

Inclusion

Figure 5: Finite Element Model used for the validation.

17

30

𝑜

𝑥1

o

𝑥1

30o 𝑥2 = 𝑥2, 𝑥3,

𝑥3

Figure 6: Local coordinate system for the model partition.

Table 1: Geometry of the bodies in contact

Region

Geometry(mm)

Sphere

R=31

Solid

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∗

18

202

of contact of 1mm (from Hertzian theory). The latter is much smaller compared with

203

the radius of the sphere (1/31) such that the half-space assumption (i.e. the contact

204

size should be small compared to the dimensions of the bodies in contact) is satisfied.

205

A local area called ’subsurface layer’ is meshed by 30, 250 linear hexahedral elements

206

(C3D8) of size 0.02a∗ , while the rest of solid is meshed with 1, 088, 349 linear tetrahedral

207

(C3D4) elements. The (quarter) sphere contains 39, 671 quadratic tetrahedral (C3D10)

208

elements and the inclusion 11, 009 linear tetrahedral (C3D4) elements.

209

Three calculation are carried out. The first case corresponds to an homogeneous

210

body to be compared with the Hertz’s solution. The second case considers an isotropic

211

tilted inclusion and is used to validate the framework for the inclusion orientation. At

212

last, the third case considers an orthotropic elastic inclusion, for which the orthotropic

213

coordinate system coincides with the contact coordinate system.

214

215

The material properties for each region of the model and for each case are shown in Table 3.

216

Results from the Semi Analytical Method and the Final Element Model are pre-

217

sented in Fig. 7 for the three cases studied. The comparison focuses on pressure

218

distribution obtained by both methods. An excellent agreement is observed between

219

both methods thus validating the SAM framework.

220

7. Parametric study

221

In this section, a normal contact between a spherical indenter and an elastic half

222

space containing a single inhomogeneity is considered. The radius of the indenter is

223

R = 62mm and the normal force P = 10000N . The Young’s modulus and Poisson’s ratio

19

Table 3: Material properties for each region and for each case

Region

Material properties

Sphere

E S = 1010 GPa,

νS = 0.3

Solid

E M = 210 GPa,

ν M = 0.3

case 1 (Isotropic)

E I = 210 GPa,

νI = 0.3

case 2 (Isotropic)

E I = 840 GPa,

νI = 0.3

case 3 (Orthotropic)

E1I = 210 GPa,

E2I = 623 GPa,

E3I = 50 GPa

I µ12 = 83 GPa,

I µ13 = 400 GPa,

I µ23 = 20 GPa

I ν12 = 0.15,

I ν13 = 0.26,

I ν23 = 0.40

Inhomogeneity

20

1.25

1.25

1

1

0.75

0.75

P/P0

P/P0 0.5

0.5

0.25

0.25

SAM result FEM result 0 −1.5

−1

−0.5

0

0.5

SAM result FEM result 1

0 −1.5

1.5

−1

−0.5

x1 /a∗

0

0.5

1

1.5

x1 /a∗

(a)

(b) 1.25

1

0.75

P/P0 0.5

0.25

SAM result FEM result 0 −1.5

−1

−0.5

0

0.5

1

1.5

x1 /a∗

(c)

Figure 7: Comparaison between the Semi Analytical Method (SAM) and the FEM; (a) Hertzian contact, (b) isotropic tilted inhomogeneity, and (c) orthotropic tilted inhomogeneity.

21

224

for the half space are chosen as E M = 210GPa and ν M = 0.3 respectively (the equivalent

225

bulk modulus and shear modulus are k M = 175GPa and µ M = 80.77GPa). The resulting

226

contact width and maximum contact pressure for the homogeneous half-space, given

227

by Hertz theory, are: 2a∗ = 2mm and P0 = 4750MPa. An ellipsoidal inhomogeneity with

228

semi-axes a1 = 0.4a∗ , a2 = 0.1a∗ , a3 = 0.1a∗ and its center located at dx3 = 0.3a∗ will be

229

considered.

230

7.1. Isotropic inhomogeneity

231

The effect of an isotropic inhomogeneity on the contact pressure distribution and

232

the subsurface stress field is here investigated. The analysis will be subdivided into

233

three cases:

234

• Effect of the Young’s moduli ratio γ = E I /E M for νI = ν M = 0.3;

235

• Effect of the bulk moduli ratio δ = k I /k M for µI = µ M = 80.77GPa;

236

• Effect of the shear moduli ratio η = µI /µ M for k I = k M = 175GPa.

237

238

239

7.1.1. Effect of the Young’s modulus The Poisson’s ratio of the inhomogeneity is chosen as νI = 0.3. For this case, the elastic moduli of the inhomogeneity are given as follows:

CiIjkl =

E I νI EI δi j δkl + (δik δ jl + δil δ jk ) I (1 + − 2ν ) (1 + νI ) νI )(1

νI = ν M and E I = γE M

(25) (26)

CiIjkl = γCiMjkl ∆C = CiIjkl − CiMjkl = (γ − 1)CiMjkl

22

(27)

Equation (10) becomes: (γ − 1)CiMjkl S klmn ε∗mn + Ci jkl ε∗kl = (1 − γ)CiMjkl ε0kl

(28)

240

The dimensionless contact pressure distribution is presented for different values of

241

γ in Fig. 8. The case of an inclusion parallel to the surface is presented in Fig. 8a

242

(θ = 0). Figure 8b presents the case for θ = 45◦ .

243

It can be observed that the magnitude of the contact pressure increases when the

244

inclusion becomes stiffer, i.e. when γ > 1. Conversely when the inclusion is softer

245

than the matrix, i.e. γ < 1, the substrate material surrounding the inhomogeneity be-

246

comes more compliant and the contact pressure gets smaller than the Hertzian pressure.

247

Meanwhile the contact area increases. In Fig. 8a, γ → ∞ corresponds to a peak mag-

248

nitude augmentation of 24.3% while γ → 0 leads to a reduction of 32.12%, compared

249

with the homogeneous half space γ = 1. For an inclined rigid inclusion (γ → ∞ and

250

θ = 45◦ , see Fig. 8b) it should be noticed a sharp increase of the maximum contact

251

pressure which is here five times the Hertz’s solution. Note also that when the inclu-

252

sion becomes infinitely soft (γ → 0), as for a void or cavity, the pressure drops locally

253

to zero.

254

The influence of the orientation angle is observed for the case of a stiff inclusion

255

(γ = 4) in Fig. 9. When θ , 0 both the surface pressure and stress component lose their

256

symmetry in the plane x1 = 0. When θ approaches 90◦ , the inhomogeneity gets very

257

closed to the surface and the peak magnitude of the contact pressure increases by more

258

than 245.47%, compared with the case θ = 0.

259

In Figures 10a and 10b, the effect of the inhomogeneity’s center location on the

260

contact pressure distribution is shown. When its dimensionless depth (dx3 /a∗ ) is greater 23

1.5

γ γ γ γ γ

1.25

→0 = 0.5 =1 =2 →∞

1

P /P0

0.75

0.5

0.25

0 −1.5

−1

−0.5

0

0.5

1

1.5

x1 /a∗ (a) 5.5

γ γ γ γ γ

5 4.5

→0 = 0.5 =1 =2 →∞

4 3.5

P/P0

3 2.5 2 1.5 1 0.5 0 −1.5

−1

−0.5

0

0.5

1

1.5

x1 /a∗ (b)

Figure 8: Effect of the dimensionless Young’s modulus γ = E I /E M (with νI = ν M = 0.3) on the dimensionless contact pressure distribution for isotropic ellipsoidal inclusion (a1 = 0.4a∗ , a2 = a3 = 0.1a∗ , dx3 = 0.3a∗ ); (a) θ = 0 and (b) θ = 45◦ .

24

4

θ θ θ θ θ

3.5 3

= 00 = 300 = 600 = 900 = 1050

2.5

P /P0 2 1.5 1 0.5 0 −1.5

−1

−0.5

0

0.5

1

1.5

x1 /a∗

Figure 9: Dimensionless contact pressure distribution for different orientation angles θ for a stiff (γ = 4) isotropic ellipsoidal inclusion (a1 = 0.4a∗ , a2 = a3 = 0.1a∗ ,dx3 = 0.4a∗ )

261

than 0.7, the effect of the inhomogeneity on the contact pressure distribution becomes

262

negligible.

263

The effect of the inhomogeneity on the subsurface stresses is now investigated and

264

the results are presented in Figs. 11 and 12. As a general trend it is observed that

265

the stress components (σi j ) inside the inhomogeneity increase as the inhomogeneity

266

becomes stiffer and at the contrary decrease as the inhomogeneity becomes more com-

267

pliant.

268

The stress at the interface between the inclusion and the matrix is of great interest to

269

study crack initiation or debonding of composite’s fibers. When the point M describes

270

the inclusion matrix interface (see Fig. 4b), the angle ξ ranges from 0◦ to 360◦ . Figures

271

13 and 14 show the normal and shear stresses at the inhomogeneity / matrix interface as

272

a function of the angle ξ (the position of the point M at this interface) for three different

273

Young’s moduli ratios γ. Note that the normal stress is almost always negative which

25

1.25

dx3 /a∗ dx3 /a∗ dx3 /a∗ dx3 /a∗ dx3 /a∗ dx3 /a∗

1

P/P0

= 0.25 = 0.3 = 0.4 = 0.5 = 0.6 = 0.7

0.75

0.5

0.25

0 −2

−1

0

1

2

x1 /a∗ (a) 1.25

dx3 /a∗ dx3 /a∗ dx3 /a∗ dx3 /a∗ dx3 /a∗ dx3 /a∗

1

P/P0

= 0.25 = 0.3 = 0.4 = 0.5 = 0.6 = 0.7

0.75

0.5

0.25

0 −2

−1

0

1

2

x1 /a∗ (b)

Figure 10: Dimensionless contact pressure distribution for different dimensionless depths dx3 /a∗ for an ellipsoidal cavity (a1 = 0.3a∗ , a2 = a3 = 0.1a∗ , γ → 0); (a) θ = 0 and (b) θ = 45◦ .

26

(a)

(b)

(c)

Figure 11: σ33 /P0 in the plane x2 = 0 for various dimensionless Young’s moduli γ = E I /E M (with νI = ν M = 0.3) for isotropic ellipsoidal inclusion (a1 = 0.4a∗ , a2 = a3 = 0.1a∗ , dx3 /a∗ = 0.4, θ = 30◦ ); (a) γ = 0.25, (b) γ = 4, and (c) γ → ∞.

27

(a)

(b)

(c)

Figure 12: σ13 /P0 in the plane x2 = 0 for various dimensionless Young’s moduli γ = E I /E M (with νI = ν M = 0.3) for isotropic ellipsoidal inclusion (a1 = 0.4a∗ , a2 = a3 = 0.1a∗ , dx3 /a∗ = 0.4, θ = 30◦ ); (a) γ = 0.25, (b) γ = 4, and (c) γ → ∞.

28

274

means compressive. One can remark that the normal shear stress at the inhomogeneity /

275

matrix interface increases when the inhomogeneity becomes stiffer. For the asymptotic

276

case of a rigid inclusion (i.e. when γ → ∞) the maximum value of the interfacial shear

277

stress reached 0.69P0 in the case of the ellipsoidal inhomogeneity (Fig. 13) while it

278

becomes slightly higher (0.86P0 ) in the case of a spherical inhomogeneity (Fig. 14).

279

7.1.2. Effect of the bulk modulus

280

The effect of the bulk modulus on the contact problem is now investigated. The

281

bulk modulus characterizes the volume variation. The inhomogeneity and matrix shear

282

moduli are set equal (µI = µM = 80.77GPa). Note that in the case of spherical inhomo-

283

geneity, the eigenstrain is purely hydrostatic. δ = k I /k M is the ratio of the inhomogene-

284

ity’s bulk modulus to the matrix one. In particular, δ → ∞ and δ → 0 correspond to a

285

Poisson’s ratio of 0.5 (i.e. incompressible) and -1 for the inhomogeneity, respectively. Since Ci jkl is an isotropic tensor, C I can be decomposed in the base of isotropic tensors J and K . CiIjkl = 3k I Ji jkl + 2µI Ki jkl

286

(29)

where

287

Ji jkl = 13 δi j δkl and Ki jkl = Ii jkl − Ji jkl

288

Ii jkl =

289

Ji jkl εkl =

1 2



 δil δ jk + δik δ jl , is the fourth-order identity tensor εkk δ 3 ij

(spherical part), and Ki jkl εkl = ei j (deviatoric part)

µI = µ M , k I = δk M

29

(30)

0.5 0 −0.5 −1

σn /P0 −1.5 −2 −2.5 −3 0

γ = 0.25 γ=4 γ→∞ 60

120

180

240

300

360

300

360

ξ(degree) (a) 1

0.5

τ /P0

0

−0.5

−1 0

γ = 0.25 γ=4 γ→∞ 60

120

180

240

ξ(degree) (b)

Figure 13: Dimensionless normal and shear stresses (σn /P0 and τ/P0 ) in the plane x2 = 0 for various dimensionless Young’s moduli γ = E I /E M in presence of a single isotropic ellipsoidal inclusion (a1 = 0.4a∗ , a2 = a3 = 0.1a∗ , dx3 /a∗ = 0.4, θ = 30◦ ); (a) σn /P0 and (b) τ/P0 .

30

0.5 0 −0.5 −1

σn /P0 −1.5 −2 −2.5 −3 0

γ = 0.25 γ=4 γ→∞ 60

120

180

240

300

360

300

360

ξ(degree) (a) 1.5

1

0.5

τ /P0

0

−0.5

−1

−1.5 0

γ = 0.25 γ=4 γ→∞ 60

120

180

240

ξ(degree) (b)

Figure 14: Dimensionless normal and shear stresses (σn /P0 and τ/P0 ) in the plane x2 = 0 for various dimensionless Young’s moduli γ = E I /E M in presence of a single isotropic and spherical inclusion (a1 = a2 = a3 = 0.2a∗ , dx3 /a∗ = 0.3); (a) σn /P0 and (b) τ/P0 .

31

∆C = C I − C M = 3(k I − k M )Ji jkl = 3(δ − 1)k M Ji jkl

(31)

Then Eq. (10) becomes

290

3(δ − 1)k M Ji jkl S klmn ε∗mn + CiMjkl ε∗kl = −3(δ − 1)k M Ji jkl ε0kl

(32)

(δ − 1)k M δi j S kkmn ε∗mn + CiMjkl ε∗kl = −(δ − 1)k M δi j ε0kk

(33)

where ε0kk represents the spherical part of ε0kl .

291

Figures 15a and 15b present the dimensionless contact pressure distribution for

292

various bulk moduli ratios δ. It can be observed that the contact pressure locally in-

293

creases on the top of the inclusion when it tends to be incompressible (i.e. when k I or

294

δ → ∞). This peak of pressure is less marked when the ellipsoidal inclusion, which

295

center is located at the dimensionless depth dx3 /a∗ = 0.3, is oriented parallel to the sur-

296

face (θ = 0, Fig. 15a) whereas it becomes sharp for the inclined inclusion (Fig. 15b).

297

This can be explained by the fact that in the latter configuration (θ = 45◦ ) the surface

298

of the inclusion becomes closer to the contact interface. Note that there is no need to

299

have a stiff inclusion to increase the maximum contact pressure, a relatively soft but

300

incompressible one will have a very similar behavior.

301

7.1.3. Effect of the shear modulus

302

The bulk modulus of the inhomogeneity is chosen equal to the matrix one, k I = k M =

303

175GPa. When the inhomogeneity is spherical, the eigenstrain is purely deviatoric.

304

The ratio of the inhomogeneity’s shear modulus to the matrix one is defined by the

305

dimensionless parameter η = µI /µ M . 32

1.25

δ δ δ δ δ

1

→0 = 1/2 =1 =2 →∞

0.75

P /P0 0.5

0.25

0 −1.5

−1

−0.5

0

0.5

1

1.5

x1 /a∗ (a) 1.5

δ δ δ δ δ

1.25

→0 = 1/2 =1 =2 →∞

1

P /P0

0.75

0.5

0.25

0 −1.5

−1

−0.5

0

0.5

1

1.5

x1 /a∗ (b)

Figure 15: Effect of the dimensionless bulk modulus δ = k I /k M on the dimensionless contact pressure distribution for a single isotropic ellipsoidal inclusion (a1 = 0.4a∗ , a2 = a3 = 0.1a∗ , dx3 = 0.3a∗ ); (a) θ = 0 and (b) θ = 45◦ .

33

The elastic modulus of the inhomogeneity is given as follow: CiIjkl = 3k I Ji jkl + 2µI Ki jkl

(34)

k I = k M , µI = ηµ M

(35)

∆C = C I − C M = 2(η − 1)µ M Ki jkl

(36)

2(η − 1)µ M Ki jkl S klmn ε∗ mn + CiMjkl ε∗kl = −2(η − 1)µ M e0i j

(37)

Equation (10) becomes

306

where e0i j represents the deviatoric part of ε0i j .

307

The effect of the dimensionless shear modulus, η, on the contact pressure distri-

308

bution is shown in Figs. 16a and 16b. It can be observed that variations of contact

309

pressure for different values of γ (see Fig. 8) and η (see Fig. 16) are quite similar.

310

7.2. Anisotropic inhomogeneity

311

The material properties of the inhomogeneity are here taken anisotropic. An or-

312

thotropic material is considered. The elastic coefficients are given in the orthotropy

313

I axes Rorthotropic as follows: E1I = 210GPa, E2I = 623GPa, E3I = 50GPa, µ12 = 83GPa,

314

I I I I I µ13 = 400GPa µ23 = 20GPa, ν12 = 0.15, ν13 = 0.26, and ν23 = 0.4. The inclusion orienta-

315

tion θ is set equal to 45◦ .

316

Euler Angles of type ZXZ are introduced (Fig. 17) in order to define the orientation

317

of the orthotropy using a combination of three rotations (ϕ,θ,ψ) around the axes of the

318

contact reference frame ( x1 , x2 , x3 ).

319

Four material orientations are defined and studied this way:

34

1.5

η η η η η

1.25

→0 = 1/2 =1 =2 →∞

1

0.75

P /P0 0.5

0.25

0 −1.5

−1

−0.5

0

0.5

1

1.5

x1 /a∗ (a) 5

η η η η η

4.5 4

→0 = 1/2 =1 =2 →∞

3.5 3 2.5

P /P0

2 1.5 1 0.5 0 −1.5

−1

−0.5

0

0.5

1

1.5

x1 /a∗ (b)

Figure 16: Effect of the dimensionless shear modulus η = µI /µ M on the dimensionless contact pressure distribution for a single isotropic ellipsoidal inclusion (a1 = 0.4a∗ , a2 = a3 = 0.1a∗ , dx3 = 0.3a∗ ); (a) θ=0 and (b) θ=45◦ .

35

x3contact

x3contact x3’

x2’

𝜃

𝜑

x3orthotropic

x3contact orthotropic x2

𝜃

x2contact contact

x1contact x1’

x2’’

x2

𝜑

x2contact

x1contact x1’

𝜑

x1contact

𝜓

x1

orthotropic

x1c’

Figure 17: Euler Angles.

320

• Case 1: ϕ = 0◦ , θ = 0◦ , ψ = 0◦ , the orthotropy and contact axes coincide.

321

• Case 2: ϕ = 90◦ , θ = 45◦ , ψ = 90◦ , the orthotropy and ellipsoid axes coincide.

322

• Case 3: ϕ = 45◦ , θ = 60◦ , ψ = 30◦ .

323

• Case 4: ϕ = −30◦ , θ = 45◦ , ψ = 30◦ .

324

The contact pressure distributions corresponding to the four cases are shown in

325

Figs. 18a and 18b.

326

One can see that depending on the orthotropic orientation relative to the contact

327

orientation, the contact pressure value can either increase (Cases 1 and 4) or decrease

328

(Cases 2 and 3). This is due to the fact that the inhomogeneity stiffness tensor associ-

329

ated to the contact axes varies widely in the four cases.

330

8. Conclusion

331

A numerical method has been proposed to model the effect of an heterogeneous

332

ellipsoidal inclusion with arbitrary orientation on the solution of a three-dimensional

333

contact problem. The elastic properties of the inhomogeneity can be either isotropic,

36

1.5

Hertz Case 1 Case 2 Case 3 Case 4

1.25

1

P /P0 0.75

0.5

0.25

0 −1.5

−1

−0.5

0

0.5

1

1.5

x2 /a∗ (a) 1.5

Hertz Case 1 Case 2 Case 3 Case 4

1.25

1

P /P0 0.75

0.5

0.25

0 −1.5

−1

−0.5

0

0.5

1

1.5

x1 /a∗ (b)

Figure 18: Dimensionless contact pressure distribution for the heterogeneous solid in presence of an orthotropic tilted ellipsoidal inclusion (a1 = 0.4a∗ , a2 = a3 = 0.1a∗ , dx3 = 0.4a∗ , θ = 45◦ ); (a) in the plane x1 = 0 and (b) in the plane x2 = 0.

37

334

orthotropic or fully anisotropic. The proposed method has been validated by perform-

335

ing a comparison with the results of a finite element model.

336

It is found that the presence of such an inclusion located below the contact surface

337

modifies significantly to very significantly the pressure distribution when its center is

338

located at a depth up to 0.7 times the contact radius. An asymmetry of the pressure

339

distribution is also observed when the axes of the ellipsoidal inclusion do not coincide

340

with the contact problem axes. Still regarding the contact pressure distribution it was

341

shown that both the dimensionless Young’s modulus of the inclusion and the dimen-

342

sionless bulk modulus have a very strong effect on the local contact pressure. In other

343

terms the compressibility properties of the inclusion are as important than the dimen-

344

sionless Young’s modulus. A local peak of pressure up to 5 times the one in the absence

345

of inclusion can be observed in some extreme configurations.

346

The normal and shear stresses at the interface between the inclusion and the matrix

347

have also been investigated. In the configurations analyzed it was found that the normal

348

stress is mostly compressive, whatever the elastic properties of the inclusion including

349

when of ellipsoidal shape. Further work is now required to study more accurately what

350

happens at the inclusion/matrix interface, as a function of the elastic properties of both

351

materials, aiming at developing a decohesion model. Finally it should be pointed out

352

that i) the method can account for the interactions between close inclusions, and ii) any

353

geometry and material properties of the contacting bodies may be considered in the

354

modeling (in other words the indenter does not have to be rigid and spherical).

38

355

References

356

Aderogba, K., 1976. On eigenstresses in a semi-infinite solid. Mathematical Proceed-

357

ings of the Cambridge Philosophical Society 80(3), 555–562.

358

Asaro, R., Barnett, D., 1975. The nonuniform transformation strain problem for an

359

anisotropic ellipsoidal inclusion. Journal of the Mechanics of Physics and Solids 23,

360

77–83.

361

Bagault, C., Nelias, D., Baietto, M.C., 2012. Contact analyses for anisotropic half

362

space: effect of the anisotropy on the pressure distribution and contact area. ASME

363

Journal of Tribology 134, 031401 (8 pages).

364

Bagault, C., Nelias, D., Baietto, M.C., Ovaert, T.C., 2013. Contact analyses for

365

anisotropic half-space coated with an anisotropic layer: Effect of the anisotropy

366

on the pressure distribution and contact area. International Journal of Solids and

367

Structures 50, 743–754.

368

Boucly, V., Nelias, D., Liu, S., Wang, Q., Keer, L., 2005. Contact analyses for bodies

369

with frictional heating and plastic behaviour. ASME Journal of Tribology 127(2),

370

355–364.

371

Chaise, T., Li, J., Nelias, D., Kubler, R., Taheri, S., Douchet, G., Robin, V., Gilles, P.,

372

2012. Modelling of multiple impacts for the prediction of distortions and residual

373

stresses induced by ultrasonic shot peening (usp). Journal of Materials Processing

374

Technology 212, 2080–2090.

39

375

Chaise, T., Nelias, D., 2011. Contact pressure and residual strain in 3d elasto-plastic

376

rolling contact for a circular or elliptical point contact. Journal of Tribology 133,

377

041402–1.

378

Chaise, T., Nelias, D., Sadeghi, F., 2011. On the effect of isotropic hardening on

379

the coefficient of restitution for single or repeated impacts using a semi-analytical

380

method. Tribology Transactions 54, 714–722.

381

Chen, W., Liu, S., Wang, Q., 2008a. Fft-based numerical methods for elasto-plastic

382

contacts of nominally flat surfaces.

383

011022–1–11.

ASME Journal of Applied Mechanics 75,

384

Chen, W., Wang, Q., Wang, F., Keer, L., Cao, J., 2008b. Three-dimensional repeated

385

elasto-plastic point contact, rolling and sliding. ASME Journal of Applied Mechan-

386

ics 75, 021021–1–12.

387

388

Chiu, Y., 1977. On the stress field due to initial strains in a cuboid surrounded by an infinite elastic space. ASME Journal of Applied Mechanics 44, 587–590.

389

Chiu, Y., 1978. On the stress field and surface deformation in a half space with a

390

cuboidal zone in which initial strains are uniform. ASME Journal of Applied Me-

391

chanics 45, 302–306.

392

Courbon, J., Lormand, G., Dudragne, G., Daguier, P., Vincent, A., 2005. Influence of

393

inclusion pairs, clusters and stringers on the lower bound of the endurance limit of

394

bearing steels. Tribology International 36, 921–928.

40

395

396

397

398

399

400

Eshelby, J., 1957. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proceedings of the Royal Society of London A241, 376–396. Eshelby, J., 1959. The elastic field outside an elastic inclusion. Proceedings of the Royal Society of London A252, 561–569. Eshelby, J., 1961. Elastic inclusions and inhomogeneities. Progress in Solid Mechanics 2, 89–140.

401

Fulleringer, B., Nelias, D., 2010. On the tangential displacement of a surface point

402

due to a cuboid of uniform plastic strain in a half-space. ASME Journal of Applied

403

Mechanics 77(2), 021014–1/7.

404

405

406

407

408

409

Gallego, L., Fulleringer, B., Deyber, S., Nelias, D., 2010a. Multiscale computation of fretting wear at the blade/disk interface. Tribology International 43, 708–718. Gallego, L., Nelias, D., 2007. Modeling of fretting wear under gross slip and partial slip conditions. ASME Journal of Tribology 129(3), 528–535. Gallego, L., Nelias, D., Deyber, S., 2010b. A fast and efficient contact algorithm for fretting problems applied to fretting modes i, ii and iii. Wear 268, 208–222.

410

Gallego, L., Nelias, D., Jacq, C., 2006. A comprehensive method to predict wear and

411

to define the optimum geometry of fretting surfaces. ASME Journal of Tribology

412

128(3), 476–485.

413

Gradshteyn, Ryzhik, I., 1965. Table of Integrals, Series and Products. Academic Press.

41

414

Jacq, C., Nelias, D., Lormand, G., Girodin, D., 2002.

Development of a three-

415

dimensional semi-analytical elastic-plastic contact code. ASME Journal of Tribol-

416

ogy 124(4), 653–667.

417

Johnson, K.L., 1985. Contact Mechanics. Cambridge University Press.

418

Kabo, E., Ekberg, A., 2002. Fatigue initiation in railway wheels-a numerical study of

419

420

421

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.

422

Kuo, C., 2007. Stress disturbances caused by the inhomogeneity in an elastic half-

423

plane subjected to contact loading. International Journal of Solids and Structures 44,

424

860–873.

425

426

Kuo, C., 2008. Contact stress analysis of an elastic half-plane containing multiple inclusions. International Journal of Solids and Structures 45, 4562–4573.

427

Leroux, J., Fulleringer, B., Nelias, D., 2010. Contact analysis in presence of spherical

428

inhomogeneities within a half-space. International Journal of Solids and Structures

429

47, 3034–3049.

430

Leroux, J., Nelias, D., 2011. Stick-slip analysis of a circular point contact between a

431

rigid sphere and a flat unidirectional composite with cylindrical fibers. International

432

Journal of Solids and Structures 48, 3510–3520.

433

434

Liu, S., Wang, Q., 2005. Elastic fields due to eigenstrains in a half-space. ASME Journal of Tribology 72, 871–878. 42

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

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.E.H., 1952. A Treatise on the Mathematical Theory of Elasticity. 4th edition ed., Cambridge University Press. Mindlin, R., Cheng, D., 1950. Thermoelastic stress in the semi-infinite solid. Journal of Applied Physics 21, 931–933. Moschovidis, Z., Mura, T., 1975. Two-ellipsoidal inhomogeneities by the equivalent inclusion method. ASME Journal of Applied Mechanics 42, 847–852. Mura, T., 1987. Micromechanics of Defects in Solids. 2nd ed., Kluwer Academic Publishers. Mura, T., Furuhashi, R., 1984. The elastic inclusions with a sliding interface. ASME Journal of Applied Mechanics 51, 308–310. Nelias, D., Antaluca, E., Boucly, V., 2007a. Rolling of an elastic ellipsoid upon an elastic-plastic flat. Journal of Tribology 129, 791–800. Nelias, D., Antaluca, E., Boucly, V., Cretu, S., 2007b. A 3d semi-analytical model for elastic-plastic sliding contacts. ASME Journal of Tribology 129(4), 761–771.

451

Nelias, D., Dumont, M.L., Champiot, F., Vincent, A., Girodin, D., Fougres, R., Fla-

452

mand, L., 1999. Role of inclusions, surface roughness and operating conditions on

453

rolling contact fatigue. ASME Journal of Tribology 121(2), 240–251.

454

Nugent, E., Calhoun, R., Mortensen, A., 2000. Experimental investigation of stress

43

455

and strain fields in a ductile matrix surrounding an elastic inclusion. Acta Materiala

456

48, 1451–1467.

457

Rodriguez-Tembleque, L., Buroni, F., Abascal, R., Saez, A., 2013. Analysis of frp

458

composites under frictional contact conditions. International Journal of Solids and

459

Structures 50, 3947–3959.

460

Seo, T., Mura, T., 1970. The elastic field in half space due to ellipsoidal inclusions

461

with uniform dilatational eigenstrains. ASME Journal of Applied Mechanics 46,

462

568–572.

463

Vignal, V., Oltra, R., Josse, C., 2003. Local analysis of the mechanical behavior of

464

inclusions-containing stainless steels under strstrain conditions. Scripta Materiala

465

49, 779–784.

466

467

468

469

Voskamp, A., 1985. Material response to rolling contact loading. ASME Journal of Tribology 107, 359–366. Walpole, L., 1967. The elastic field of an inclusion in an anisotropic medium. Proceedings of the Royal Society of London A300, 270–288.

470

Wang, F., Block, J., Chen, W., Martini, A., Keer, L., Wang, Q., 2009. A multi-scale

471

model for the simulation and analysis of elasto-plastic contact of real machined sur-

472

faces. ASME Journal of Tribology 131, 021409–1–6.

473

474

Wang, F., Keer, L., 2005. Numerical simulation for three dimensional elastic-plastic contact with hardening behavior. ASME Journal of Tribology 127, 494–502.

44

475

476

Willis, J., 1964. Anisotropic elastic inclusion problems. The Quarterly Journal of Mechanics and Applied Mathematics 17(2), 157–174.

477

Zhou, K., Chen, W., Keer, L., Ai, X., Sawamiphakdi, K., Glaws, P., Wang, Q., 2011a.

478

Multiple 3d inhomogeneous inclusions in a half space under contact loading. Me-

479

chanics of Materials 43, 444–457.

480

Zhou, K., Chen, W., Keer, L., Wang, Q., 2009. A fast method for solving three-

481

dimensional arbitrarily shaped inclusions in a half space. Computer Methods in

482

Applied Mechanics and Engineering 198(9-12), 885–892.

483

484

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. Mechanics of Materials 60, 144–158.

485

Zhou, K., Keer, L., Wang, Q., 2011b. Semi-analytic solution for multiple interacting

486

three-dimensional inhomogeneous inclusions of arbitrary shape in an infinite space.

487

International Journal for Numerical Methods in Engineering 87, 617–638.

488

Zhou, K., Keer, L., Wang, Q., Ai, X., Sawamiphakdi, K., Glaws, P., Paire, M., Che, F.,

489

2012. Interaction of multiple inhomogeneous inclusions beneath a surface. Com-

490

puter Methods in Applied Mechanics and Engineering 217-220, 25–33.

491

Appendix A. Stress in a half-space due to a concentrated unit normal force at the

492

surface origin( Fi j )

F11 (x1 , x2 , x3 ) =

" # 1 1 − 2ν x3 x12 − x22 x3 x22 3x3 x12 (1 − ) + − , 2π r2 ρ r2 ρ3 ρ5

45

F22 (x1 , x2 , x3 ) = F11 (x2 , x1 , x3 ),

F33 (x1 , x2 , x3 ) = −

F12 (x1 , x2 , x3 ) =

3 x33 , 2π ρ5

" # 1 1 − 2ν x3 x1 x2 x3 x2 x1 3x3 x2 x1 (1 − ) + − , 2π r2 ρ r2 ρ3 ρ5

F13 (x1 , x2 , x3 ) = −

3 x1 x32 , 2π ρ5

F23 (x1 , x2 , x3 ) = F12 (x2 , x1 , x3 ),

where r2 = x12 + x22 , 493

494

ρ=

q

x12 + x32 + x32 ,

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

Appendix B. Stresses in a half-space subject to normal pressure (Mi j )

495

An isotropic half-space is submitted an uniform normal pressure σn in a discretized

496

surface area of 2∆x1 × 2∆x2 at the center point P(x10 , x20 , 0). The stress at an observation

497

point Q(x1 , x2 , x3 ) is given by Zhou et al. (2009) and Johnson (1985):

σi j (x1 , x2 , x3 ) = Mi j (x1 − x10 , x2 − x20 , x3 )σn (x1 , x2 ) σi j (x1 , x2 , x3 ) =

σn [hi j (ξ1 + ∆x1 , ξ2 + ∆x2 , ξ3 ) − hi j (ξ1 + ∆x1 , ξ2 − ∆x2 , ξ3 ) 2π + hi j (ξ1 − ∆x1 , ξ2 − ∆x2 , ξ3 ) − hi j (ξ1 − ∆x1 , ξ2 + ∆x2 , ξ3 )]

46

where ξi = xi − xi0 .

The functions hi j () in Eq.(B1) are defined by h11 (x1 , x2 , x3 ) = 2ν tan−1

x22 + x32 − ρx2 ρ − x2 + x3 x1 x2 x3 + 2(1 − ν) tan−1 + , x1 x 3 x1 ρ(x12 + x32 ) h22 (x1 , x2 , x3 ) = h11 (x2 , x1 , x3 ),

h33 (x1 , x2 , x3 ) = tan−1

! x22 + x32 − ρx2 x1 x2 x3 1 1 − + , x 1 x3 ρ x12 + x32 x22 + x32

h12 (x1 , x2 , x3 ) = −

x3 − (1 − 2ν) log(ρ + x3 ), ρ

h13 (x1 , x2 , x3 ) = −

x2 x32 ρ(x13 + x32 )

,

h23 (x1 , x2 , x3 ) = h13 (x2 , x1 , x3 ),

where ρ=

498

q

x12 + x22 + x32 .

Appendix C. Normal displacement at the surface subject to normal pressure (K n)

499

The contact between a sphere and an elastic half-space having respectively elastic

500

constants (E1 , ν1 ) and (E2 , ν2 ), where the surface x3 = 0 is discretized into rectangular

501

surface area of 2∆1 × 2∆2 , is now considered. The initial contact point coincides with

502

the origin of the Cartesian coordinate system (x1 , x2 , x3 ). The relationship between the

503

normal displacement at an observation point P(ξ1 , ξ2 , 0) and the pressure field at the

504

center Q(ξ10 , ξ20 , 0) is built using the function K n .

" K n (c1 , c2 ) =

# 4 1 − ν12 1 − ν22 X + K n (c1 , c2 ), πE1 πE2 p=1 p

47

  (c2 + ∆2 ) + K1n (c1 , c2 ) = (c1 + ∆1 ) log  (c2 − ∆2 ) +   (c1 + ∆1 ) + K2n (c1 , c2 ) = (c2 + ∆2 ) log  (c1 − ∆1 ) +   (c2 − ∆2 ) + K3n (c1 , c2 ) = (c1 − ∆1 ) log  (c2 + ∆2 ) +   (c1 − ∆1 ) + K4n (c1 , c2 ) = (c2 − ∆2 ) log  (c1 + ∆1 ) +

 p (c2 + ∆2 )2 + (c1 + ∆1 )2   , p (c2 − ∆2 )2 + (c1 + ∆1 )2  p (c2 + ∆2 )2 + (c1 + ∆1 )2   , p (c2 + ∆2 )2 + (c1 − ∆1 )2  p (c2 − ∆2 )2 + (c1 − ∆1 )2   , p (c2 + ∆2 )2 + (c1 − ∆1 )2  p (c2 − ∆2 )2 + (c1 − ∆1 )2   , p (c2 − ∆2 )2 + (c1 + ∆1 )2

where c1 = ξ1 − ξ10

505

and c2 = ξ2 − ξ20

Appendix D. Euler rotation matrix

   cos ϕ    =  sin ϕ RZ (ϕ)     0

− sin ϕ cos ϕ 0

   cos ψ    RZ (ψ) =  sin ψ     0    1    RX (α) =  0     0

− sin ψ cos ψ 0

0 cos α sin α

  0     0     1        0     1  0

   0    − sin α     cos α 

P(ϕ, α, ψ) = RZ (ϕ)RX (α)RZ (ψ)

48

Nomenclature

Letters a∗

contact radius for contact problem

a1 ,a2 ,a3

semi-axes of an ellipsoidal inhomogeneity

B∗i jkl

influence coefficients that relating the stress σi j at point (x13 , x2 , x3 ) to the constant eigenstrain at the point (x1k , x2k , x3k )

CiMjkl , CiIjkl

elastic constants of the matrix and the inhomogeneity

e0i j

deviatoric part of ε0i j

E M ,E I

Young’s modulus of the matrix and the inhomogeneity

E(θ0 , k), F(θ0 , k)

function of elliptic integral

h

distance between the two surfaces of the contacting bodies

Ii jkl

the fourth-order identity tensor

Kn

coefficients in the normal displacement at the contact surface due to the contact pressure

k M , kI

bulk modulus of the matrix and the inhomogeneity

L1 , L2 , L3

lengths of the three sides of the matrix in EF model

Mi j

influence coefficients relating the stress σi j at the point (x1 , x2 , x3 ) to the normal traction σn within a discretized area centered at (x1k , x2k , 0)

n1 , n2 , n3

grid number in the half-space along the Cartesian directions x1 , x2 , x3 , respectively

P

normal applied load

P0

maximum Hertzian pressure

p R

contact pressure distribution 49 indenter radius

S i jkl

components of the Eshelby’s tensor

u0i

displacements corresponding to the infinite applied strain ε0i j

ui

disturbed contribution of the displacements

Nomenclature

Letters W

applied exterior load

dx3

depth of the inclusion from the surface of the matrix in EF model

xI = (x1I , x2I , x3I )

Cartesian coordinates of the inclusion center

Greek letters ε0i j

infinite applied strain

εi j

strain due to eigenstrains

ε∗i j

eigenstrain due to the presence of inhomogeneities

ε0kk

spherical part of ε0i j

εp

initial eigenstrain of the inhomogeneity

σ0i j

stress corresponding to the infinite applied strain ε0i j

σi j

disturbed contribution of the stresses

φ, Ψ

harmonic and biharmonic potentials of mass density ε∗i j

φi j .., Ψi j ..

harmonic and biharmonic potentials of mass density xi x j . . .

δi j

Kronecker symbol

σn

normal pressure due to the summution of both symmetric inclusions

∆x1 , ∆x2

half size of the discretized surface area

ν M , νI

Poisson’s ratio of the matrix M and the inclusion I

µ M , µI

shear modulus of the matrix and the inclusion

γ

the ratio of the inhomogeneity’s Young’s modulus to the matrix’s

δ η

the ratio of the inhomogeneity’s bulk modulus to the matrix’s 50 the ratio of the inhomogeneity’s shear modulus to the matrix’s

θ

the tilted angle of the inhomogeneity in the x1 Ox3 plan

Nomenclature

Letters Acronyms and f ast Fourier trans f orms

2D-FFT

two-dimensional fast Fourier transform

3D-FFT

three-dimensional fast Fourier transform

FFT −1

inverse FFT operation

b Bi jkl

frequency response of coefficients Bi jkl in the frequency domain

bi j M

frequency response of coefficients Mi j in the frequency domain

51