Stick-slip analysis of a circular point contact between a rigid sphere and a flat unidirectional composite with cylindrical fibers

International Journal of Solids and Structures 48 (2011) 3510–3520

Contents lists available at SciVerse ScienceDirect

International Journal of Solids and Structures journal homepage: www.elsevier.com/locate/ijsolstr

Stick-slip analysis of a circular point contact between a rigid sphere and a flat unidirectional composite with cylindrical fibers Julien Leroux, Daniel Nélias ⇑ Université de Lyon, CNRS INSA-Lyon, LaMCoS UMR5259, F-69621, France

a r t i c l e

i n f o

Article history: Received 2 June 2011 Received in revised form 16 August 2011 Available online 16 September 2011 Keywords: Contact mechanics Stick-slip Fretting Composite material Eshelby’s equivalent inclusion method

a b s t r a c t The stick-slip contact problem is investigated here when at least one of the contacting bodies behaves as an ideal composite material with long fibers perpendicular to the direction of movement. Cylindrical inhomogeneous inclusions within a homogeneous media and with axes parallel to the contact surface are considered. The Eshelby’s equivalent inclusion method is used to solve the problem numerically. Interactions between close inclusions are taken into account in the numerical procedure, as well as the coupling between the normal and tangential contact problems. It is found that the presence of heterogeneities in the vicinity of the surface contact affects significantly the contact pressure distribution and subsequently the distribution of shear and slip at the interface. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Damage, such as loss of material or/and cracking, induced by the small amplitude oscillatory relative movement between contacting surfaces is usually termed fretting. Three different fretting modes have been defined by Mohrbacher et al. (1995), depending on the forces and moments transmitted through the contact, as shown in Fig. 1. Slips are transversal in fretting mode I (tangential loading); slips are radial in fretting mode II (normal loading); and slips are circumferential in fretting mode III (torsional loading). For many industrial applications where two bodies are put into contact, fretting damage often occurs as the result of combined fretting modes (Ruiz et al., 1992). Based on experimental observations all along the life of a contact, Zhou and Vincent (1993) have established three fretting regimes: the partial slip regime, the gross slip regime and the mixed fretting regime which corresponds to sliding conditions evolving from gross slip to partial slip. The latter can be explained by a modification of the contact conditions: wear of the surfaces and superficial material properties change such as it leads to an increase of the friction coefficient. According to Vincent et al. (1992), the gross slip regime is confronted to wear while the partial slip regime to cracks. Cracks and wear compete on the mixed regime. Recently Gallego et al. (2010a) have presented a numerical technique based on semi-analytical methods in order to simulate fretting wear under gross, mixed and partial slip conditions. Each body

⇑ Corresponding author. Tel.: +33 4 72438490; fax: +33 4 72438913. E-mail address: [email protected] (D. Nélias). 0020-7683/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijsolstr.2011.09.007

was considered homogeneous and purely elastic. The efficiency of the method has been highlighted but the analysis limited to homogeneous materials so far. The intend is here to demonstrate that the presence of heterogeneous inclusions, also called inhomogeneities in the sense of Eshelby (1961), may significantly modify the displacement and stress fields within the contacting bodies and subsequently the slip and shear distributions at the interface. Many experimental investigations for the fretting behaviour of composite materials have been published in the recent years, see for example Xu et al. (2007), however very little is known about the physics involved in the contact problem. Some elements of response can be brought through numerical modelling. For instance Kuo (2008) and Leroux et al. (2010) have highlighted that the presence of inhomogeneities in the vicinity of the contact surface could significantly change the contact pressure distribution which subsequently affects the subsurface stress distribution. A first attempt is made here to describe how the presence of long fibers may affect the pressure, shear and slip distributions in a fretting contact.

2. Theoretical background The method proposed to solve the 3D contact problem in presence of inelastic strains due to the presence of heterogeneous inclusions is based on semi-analytical methods. Reusing the theoretical background and some of the numerical techniques initially proposed by Jacq et al. (2002) for solving the 3D contact problem for elastic–plastic materials, by Boucly et al. (2005) for thermalelastic–plastic contact and Boucly et al. (2007) for the collision between asperities, the method has been later extended to inelastic

J. Leroux, D. Nélias / International Journal of Solids and Structures 48 (2011) 3510–3520

3511

Fig. 1. Three different fretting modes: tangential fretting mode I, radial fretting mode II, and circumferential fretting mode III.

eigenstrains (Zhou et al., 2009; Leroux et al., 2010; Zhou et al., 2011) with, most of the time, little change in the formulation. Note that in Leroux et al. (2010) and Zhou et al. (2011) the presence of inhomogeneities – spherical for the former, cuboidal for the latter – modifies the contact pressure distribution. This has been already observed by Kuo (2007, 2008) but for the two-dimensional case (plane strain) in presence of cylindrical heterogeneous inclusions. One of the difficulties in presence of multiple inclusions (from now heterogeneous inclusions or inhomogeneities in the sense of Eshelby will be referred as inclusions for simplicity) is that interactions between them when they are getting closer should be carefully considered. Note that in Leroux et al. (2010) the size and location of each inclusion is independent of the mesh, which can be seen as an enrichment technique, whereas here and in the work of Liu and Wang (2005) and Zhou et al. (2009, 2011) each inclusion is discretized into unit cuboids following the mesh. The algorithm developed to solve the incremental contact problem is presented Fig. 2. The theory relies on the formulation of the Betti’s reciprocal theorem. The numerical scheme is based on the use of discrete convolution fast Fourier transform (DC-FFT) and conjugate gradient method (CGM), as proposed by Liu and Wang (2001). Step 1: The initial state may include residual strains (for example thermoelastic strains) and/or initial stresses (such as those due to the processes). The elastic–plastic algorithm requires first a solver for the elastic contact problem with any initial surface separation, in which the initial geometry may be modified to account for the permanent deformation of the surface due to residual strains. In a more conventional manner it requires also the knowledge of the elastic displacements of the surface, regularly updated. Steps 2–4: The Conjugate Gradient Method (CGM) is used to solve the normal and tangential contact problems, see Polonsky and Keer (1999) for the normal problem and Gallego et al. (2010b) for the normal and tangential problem, but both problems are solved separately. The discrete convolution and fast Fourier transforms technique (DC-FFT) are used to accelerate the numerical summations (Liu et al., 2000). Equilibrium equations are also useful and can be found in Gallego et al. (2006), and Gallego and Nélias (2007) when the center of pressure is not known in advance (or when both forces and moments are imposed to the contact). The process consists of first

solving the normal problem while the shear tractions are supposed nil and the normal residual displacements are enforced. Then the shear tractions can be solved while the tangential residual displacements are now enforced and pressure is known. The first step is repeated with the new shear tractions previously found, etc. until convergence is reached. This process is known as the Panagiotopoulos process (Panagiotopoulos, 1975) and is still valid for elastic–plastic contacts or the problem with inelastic eigenstrains. It is noted that the coupling between normal and tangential effects in stick/slip contact problems should be taken into account when materials are dissimilar (Gallego et al., 2010b). Once the contact problem has converged, the stress calculation is launched considering the contribution of residual (or inelastic) strains, contact pressure and shears distributions, producing new normal and tangential (residual) surface displacements that should be added to the initial surface geometry. Steps 12–14: The plasticity model is then used to calculate the plastic strain increment if any, based on the return mapping algorithm (Ortiz and Simo, 1986; Simo and Taylor, 1986; Fotiu and Nemat-Nasser, 1996) in the framework of semi-analytical methods (Nélias et al. (2006)). The increments of normal and tangential surface displacements are calculated analytically by the expressions given by Jacq et al. (2002) and Fulleringer and Nélias (2010) for normal and tangential displacements, respectively. This process is repeated until the residual displacement increment converges. Plastic strains, load, contact pressure, residual surface displacement and hardening parameters are then incrementally updated to define a new initial condition for the next loading step. Note that the effect of plasticity has not been investigated here. However if one considers that plastic strain is only one particular type of inelastic strain, it is interesting to see how both eigenstrains and plastic strains can be treated numerically in the same contact solver. Moreover this gives an overview on how the contact solver is organized when both inhomogeneous and inelastic materials have to be considered. Compared to the widely used finite element method (FEM) this numerical technique has two main advantages: first the robustness of

3512

J. Leroux, D. Nélias / International Journal of Solids and Structures 48 (2011) 3510–3520

Initial state 1 Normal Problem

Tangential Problem 3

2

Equivalent Inclusion model

Subsurface stresses Equivalent elastic problem

4 Determination of eigenstrainsε* ε*

For two mirror -image eigenstrains

5 Stresses on a extended surface

Disturbance stresses 7

6 Misfit displacements induced by inclusions 8

Misfit displacements induced by the pressure field 9 HalfHalf-space subject to pressure

Sum of stresses and displacements 10

Convergence on displacements

Interactions between inhomogeneities

Applied strain

11 16 Plastic strain and residual stresses 12

Interpolation 17

Residual displacements 13

For two mirror-image eigenstrains

Disturbance stresses

Stresses on a extended surface

18

Convergence on residual displacements

Relaxation

Eigenstrains

19

14

Convergence on eigenstrains

End of load

Load increment 15

20 HalfHalf-space subject to pressure

21 End Fig. 2. Flowchart of the contact algorithm.

the contact solver and second the low memory space required which, combined with an efficient numerical procedure, allows to solve transient and inelastic 3D problem in less time than for a 2D problem by FEM. The contact solver may be also coupled with FEM to account for the flexibility of the structure and used as a structural zoom when a finer mesh is needed for the contact (Gallego et al., 2010a). 3. Inhomogeneous aspects An inclusion module has been added to the inelastic contact solver in order to take into account the presence of inhomogeneous inclusions of arbitrary shape in an isotropic half space, hence modifying both the residual surface displacements and the subsurface stress and strain distributions. Step 5 – Determination of eigenstrainsThe equivalent inclusion method (EIM) proposed by Eshelby (1957, 1959, 1961) will be used here to model a matrix containing a single non-homogeneous inclusion. In this case the subdomain X that contains the inhomogeneity is occupied by a different elastic

material than the matrix. The stiffness tensors of the matrix and inhomogeneity materials are Cm and CI, respectively. The matrix is subjected to an applied strain e1, uniform far from the inhomogeneity (i.e. at infinity). The stress disturbance due to the presence of the inhomogeneity can be simulated by that of an appropriately chosen equivalent transformation strain eij (eigenstrain) distributed over an inclusion X which have the same elastic material properties than the half-space. It is straightforward to verify that the necessary and sufficient condition for the equivalence between the fields induced by the inhomogeneity and the inclusion is

Cm ijkl







I  1 e1 kl þ ekl  ekl ¼ C ijkl ekl þ ekl



in X

ð1Þ

The above equation is referred to as the ‘equivalency equation’ and can be rewritten as  1 DC ijkl ekl  C m ijkl ekl ¼ DC ijkl ekl

Cm ijkl

in X

ð2Þ

C Iijkl .

where DC ijkl ¼  The ‘equivalent inclusion method’ requires the determination of the eigenstrain eij which is related to compatibility strain eij by:

eij ¼ Sijkl ekl in X where Sijkl is the Eshelby tensor.

ð3Þ

3513

J. Leroux, D. Nélias / International Journal of Solids and Structures 48 (2011) 3510–3520

Eqs. 2 and 3 – which are valid only within the inclusion domain X – form a closed system of algebraic equations for the eigenstrain. The Eshelby’s tensor is only a function of the matrix Poisson’s ratio and the geometrical aspect ratios. The various components of the Eshelby’s tensor are well documented by Mura (1987, 1988) for the general ellipsoid, for various specific subsets of the ellipsoid (i.e., elliptic cylinder, flat ellipsoid, etc.), as well as other shapes such as the parallelepiped. In this paper, only formulations for cuboids will be used in order to break up an inhomogeneity of any shape into small elementary cuboids. When many inhomogeneities are considered, the eigenstrain eij will be estimated for each inhomogeneity whenever the applied strain is modified as followed:



eij ¼  DC ijkl  Sklmn  dkm  dln  C mijkl

1

 DC ijkl  e1 kl

ð4Þ

This approach is purely analytical and uses matrix computations, except for the inversion of the first term in Eq. (4) which is made numerically by the Gauss–Seidel method. However, the Eshelby’s theory considers that each inclusion is located within a uniform applied strain field, which prevents considering interactions among them. Therefore the method should be modified to account for high gradient strains due to contact loading or when two inhomogeneities are sufficiently close to interact each other. A polynomial distribution of eigenstrains was found to be useful in dealing with the elastic interactions between inhomogeneities (Sendeckyj, 1967; Moschovidis, 1975; Rodin and Hwang, 1991). In particular, Eshelby (1961) demonstrated that if eij is a polynomial of degree n in Cartesian coordinates, then the induced strain inside the inclusion is a sum of monomials of degree n, n  2, n  4, etc. Since there is no general analytical solution to this problem, Moschovidis and Mura (1975) proposed to make a quadratic expansion of the applied strain and the eigenstrain around the inclusion centers of coordinates xI:









1 1 1 I I I e1 ij ðxÞ ¼ Eij þ Eijk xk  xk þ Eijkl xk  xk xl  xl       eij ðxÞ ¼ Eij þ Eijk xk  xIk þ Eijkl xk  xIk xl  xIl





eij ðxÞ ¼ eij ðxI Þ þ eij;k ðxI Þ xk  xIk þ



   1 eij;kl ðxI Þ xk  xIk xl  xIl 2!

ð5Þ ð6Þ ð7Þ

The coefficients Eij..m and E1 ij::m are symmetrical constants with respect to the free indices i, j since the strain is symmetric in i and j, and having values independent of the order in which the summation indices appear, i.e., Eijkl = Eijlk, Eijklm = Eijkml, etc. After a strenuous development, the coefficients of the power series are equated. Since the E tensors are symmetric, the quadratic approximation of the solution for n inhomogeneities requires determining 60n unknowns per tensor. Although in this case the mathematical expressions are somewhat lengthy, the structure of the approach is essentially the same as for a uniform eigenstrain. This generalization for the inclusion problem immediately suggests that the inhomogeneity problem can also be extended to the general case of a polynomial applied strain. A very clear presentation of these topics is given elsewhere by Moschovidis (1975). Thus this method allows accounting for a non-uniform applied strain which generates a nonuniform eigenstrain. Nevertheless it is important to consider all mutual influences between multiple inhomogeneities (steps from 16 to 21). The implementation of the method which allows taking into account the mutual effects between inhomogeneities is as follows. The applied strain field e1 is determined from the stress calculation (step 4). When this one cannot be approximated by a uniform field (step 16) it can also be fitted by a second order expansion around the center of the inclusion (step 17). The individual eigenstrains and the induced strains are calculated from the system of Eqs. (4) and (6), respectively (steps 18–20). It is important to not

consider the influence of non-homogeneous inclusion on itself by cancelling the influence coefficients within each inclusion. These induced strains are added to the initial applied strain and the individual eigenstrains are then re-calculated. These eigenstrains are compared to the ones found during the previous step (step 21). This process is repeated until each individual eigenstrain converges. Once e⁄ is estimated (step 4), strains and stresses inside and outside the inhomogeneity can be found (steps from 6 to 10). Steps (6 and 7) – Determination of induced strainThen Eqs. (2) and (3) allow to determine the induced strains eij and/or stress rij inside the inclusion. Mura (1987) describes the mathematical operator based on Green’s functions that allow determining eij and rij inside and outside the inclusion as follows:

eij ðxÞ ¼ Dijkl ðxÞEkl þ Dijklq ðxÞEklq þ Dijklqr ðxÞEklqr

ð8Þ

where Dijkl..m is a complex influence coefficient derived from harmonic U and biharmonic w potentials, and constructed to be symmetric with respect to i, j and k, l.

  1  2mdkl /q::m;ij  ð1  mÞ /q::m;kj dil w 8pð1  mÞ q::m;klij  þ/q::m;ki djl þ /q::m;lj dik þ /q::m;li djk

Dijklq::m ¼

ð9Þ

where

/ij::k ðxÞ ¼

Z Z Z

0

x0i x0j ::x0k dx x0 j

X

wij::k ðxÞ ¼

jx  Z Z Z

and 0

x0i x0j ::x0k jx  x0jdx

ð10Þ

X

m is the Poisson’s ratio, dij is the Kronecker  delta  (dij = 1 for i = j = 1, 2, 3 and dij = 0 for i – j), and jx  x0 j2 ¼ xi  x0i xi  x0i . The potential functions U and w defined by Eq. (10) have been obtained by MacMillan (1958) when the inclusion domain X is cuboidal (see Appendix A). The functions Dijkl which have been calculated at the inclusion center xI constitute the Eshelby tensor Sijkl. After this large system has been solved, the induced strain and stress tensors can be computed at any point in the elastic medium. However, the formalism considered above is validated only for infinite spaces. The three-dimensional contact problem involves to use a halfspace which is bounded by the surface plane x3 = 0 in the Cartesian system (x1, x2, x3). Mura and Cheng (1977) and Mura and Seo (1979) worked on such application, trying to extend the equation given by Moschovidis – relating induced strains and eigenstrains – to infinite half-spaces. Despite of its complexity, those solutions consider only an inclusion with hydrostatic eigenstrains. It is obvious that hydrostatic eigenstrains are not necessarily encountered in inclusions under contact conditions. Even if those solutions are not fruitful here, Zhou et al. (2009) introduced a fast method for solving the problem of three-dimensional arbitrarily shaped inclusions in an isotropic half-space, as done by Jacq et al. (2002) for each plastic cuboid when numerically solving the elastic–plastic contact problem. The common method is to decompose an isotropic infinite half-space problem into three subproblems, as shown in Fig. 3. The first sub-problem considers  the same inclusion but in an infinite space e ¼ e11 ; e22 ; e33 ; e12 ;   e13 ; e23 Þ. The second one considers an image counterpart – an inclusion as seen in a mirror – in the same infinite space. Its eigenstrain is modified in order to respect the conditions of symmetry, i.e.   es ¼ e11 ; e22 ; e33 ; e12 ; e13 ; e23 . The solution for one inclusion in a uniform strain field and in an infinite space as discussed above is fully analytical, while when in a half-space the procedure is partially numerical. The solution for one inclusion in an infinite half space is the sum of the solutions for this inclusion and its mirror image in the infinite space. It is

3514

J. Leroux, D. Nélias / International Journal of Solids and Structures 48 (2011) 3510–3520

ε *s σn

x

x 1 ,x 1 ’

O

’ ,x 2

2

=

x 1 ,x 1 ’

O x

’ ,x 2

2

+

x

−

x 1 ,x 1 ’

O ’ ,x 2

2

x 1 ,x 1 ’ x x 2,

’

O

2

ε* x 3 ,x 3 ’

x 3 ,x 3 ’

x 3 ,x 3 ’

x 3 ,x 3 ’

(1)

(3)

(2)

Fig. 3. Decomposition of the problem for an inclusion containing eigenstrain e⁄ in an isotropic half-space into three subproblems.

now necessary to determine the induced stresses all over the infinite space for two inclusions of eigenstrains e⁄ and es . The summation of the solutions of problems (1) and (2) in Fig. 3 generates only normal traction rn with no tangential traction at the surface x3 = 0. Adding an opposite normal stress rn at the surface of the half-space leads to the creation of a free surface. This pressure field comes from the knowledge of r33 produced by both inclusions









rn ðx1 ; x2 ; 0Þ ¼ B33kl x1  xI1 ; x2  xI2 ; xI3 ekl xI1 ; xI2 ; xI3      B33kl x1  xI1 ; x2  xI2 ; xI3 eskl xi1 ; xi2 ; xi3

ð11Þ

sufficiently small. Note that, in the case of plastic strains (i.e. ekk = 0), the exact solution for the normal stress distribution is given by Jacq et al. (2002) in the form of Bessel’s functions. The same approach can be used when many inhomogeneities are present by simply superimposing the effects of each one. Note that a new class of inclusions – i.e. of the same size, shape, and elastic constants – will require the computation of a new set of influence coefficients for the calculation of eigenstrains. The induced stress field r becomes:

rij ðx1 ;x2 ;x3 Þ ¼

XXX x3

þ where Bijkl are the influence coefficients that relate the induced stress rij at point (x1, x2, x3) to the eigenstrain at the inclusion center   of coordinate point xI ¼ xI1 ; xI2 ; xI3 . Finally the stress components at each point of the infinite half-space is the sum of the stresses created by the first inclusion (1), the image counterpart (2) and the normal stress distribution at the surface x3 = 0 emulating a free surface (3). 

 

rij ðx1 ;x2 ;x3 Þ ¼ Bijkl x1  xI1 ;x2  xI2 ;x3  xI3 e xI1 ;xI2 ;xI3     þ Bijkl x1  xI1 ;x2  xI2 ;x3 þ xI3 es xI1 ;xI2 ;xI3 

Z

þ1 1

Z

þ1

1



 0 0    K ij x1  x01 ;x2  x02 ;x3 rn x01 ;x02 ;0 dx1 dx2 ;ði;j;k;l ¼ 1;2;3Þ:

ð12Þ

One of the problems encountered by Mura is the integration of this third term while rn is a complex function related to the influence of both inclusion on the stress r33 at the interface. For this reason, it was easier to solve this integral with hydrostatic eigenstrains. However the normal stress distribution rn can be approximated numerically using influence coefficients as for a contact problem and in a fast and efficient way using numerical methods such as a 2D-FFT. This constitutes a semi-analytical approach removing the restrictions on the eigenstrains mentioned above. The accuracy of the solution depends on the estimation of the normal stress field rn. This estimation is most likely accurate on the domain considered – i.e. the domain meshed explicitly – whereas outside of this domain the pressure is implicitly equal to zero. Therefore the surface of the half-space should be sufficiently enlarged to keep the numerical error due to the truncature



x2

x1

    Bijkl x1  x01 ;x2  x02 ;x3  x03 ekl x01 ;x02 ;x03

XXX x3

x2

x2

x1

XX

x1

    Bijkl x1  x01 ;x2  x02 ;x3 þ x03 eskl x01 ;x02 ;x03

    Nij x1  x01 ;x2  x02 ;x3 rn x01 ;x02 ;0

ð13Þ

where the triple summations are done on the domain considered. 2D-FFT is commonly used to transform a double summation into a convolution product for contact problems. Also, this method has been used extensively in previous plastic solvers in order to transform a triple summation into a simple summation of convolution products, used for residual stresses induced by plasticity (Wang and Keer, 2005). It is called the 3D-FFT, and transforms a triple summation into a convolution product. More details on how can be calculated Bijkl, Kij and Nij in Eqs. (11)–(13) can be found in Leroux et al. (2010). Steps (8, and 9) – Determination of misfit displacementsInserting the eigenstrain into the total strain, the displacement on the contact surface due to the eigenstrain (also called the eigendisplacement) can be calculated. The normal displacements are produced by the normal stress field rn, while the tangential displacements are due to the presence of both inclusions plus the effect of the normal stress field.

ui ðx1 ;x2 Þ ¼

XXX x3

þ þ

x

x

    Dsikl x1  x01 ;x2  x02 ;x3  x03 ekl x01 ;x02 ;x03

2 1 XX X

x

x

x2

x1

3 2 X X

x1

    Dsikl x1  x01 ;x2  x02 ;x3 þ x03 eskl x01 ;x02 ;x03

    K pi x1  x01 ;x2  x02 rn x01 ;x02 ;0 ; ði ¼ 1;2Þ

ð14Þ

3515

J. Leroux, D. Nélias / International Journal of Solids and Structures 48 (2011) 3510–3520

u3 ðx1 ; x2 Þ ¼

XX x1

 p

 0

K 3 x1  x01 ; x2  x2

r

 n

 x01 ; x02 ; 0

ð15Þ

x2

where

Dsikl ¼

   1 w  2mdkl /;i  2ð1  mÞ /;k dij þ /;j dik 8pð1  mÞ ;kli

ð16Þ

The effect of a uniform pressure acting on a rectangular area has been analysed in detail by Love (1927) and Johnson (1985). Here, Kp are the influence coefficients that relate the tangential and normal displacements at the surface point (x1, x2, 0) to the normal   stress (or pressure) rn at the surface point x01 ; x02 ; 0 , given in Appendix B. The surface contact geometry is modified by adding the normal displacement and the contact pressure distribution are updated (Steps 10 and 11). The elastic displacements are calculated from the updated contact pressure. This procedure is repeated until convergence is obtained. A close-loop links the variations of the surface geometry, contact pressure and total strain, and is repeated until the difference between the eigen-displacements of two following iterative steps is less than the prescribed error. A homogeneous equivalent problem is found. Once it is done, the plastic solver can be activated if yielding occurs. It should be noted that three loops are considered for each loading increment: (i) the normal and tangential contact problem, (ii) the homogeneous equivalent problem and (iii) the plastic solver. Those loops are time consuming but no other option has been considered yet.

homogeneous half-space, this load leads to a contact radius a = 3.01 lm and a maximum contact pressure P0 = 4213 MPa. The rigid body displacement is varying between ux = ±0.04 lm. The friction coefficient l is here taken equal to 0.3. The incremental load path to describe a full permanent fretting loop is given in Fig. 4. From now, cylindrical fibers with axes perpendicular to the direction of the tangential load as shown in Fig. 5a will be considered. These ones are discretized into many unit cuboidal elements with the same size as the numerical discretization mesh. Four elements are used along the fiber radius r (eight for the diameter), with r/a = 0.13, the fibers being equally spaced with a distance equals to r between them and with the top surface located at a distance r/2 from the free surface, see Fig. 5b. Each element is treated as a cuboidal inhomogeneity twice stiffer than the matrix but with the same Poisson’s ratio. The pressure and shear distributions are plotted in Fig. 6 for different points of the fretting loop, along the x-axis in the left column (a) and the y-axis in the right one (b). For comparison the distribution of contact pressure and shears for the homogeneous half-space are also plotted with dotted lines. The contact pressure and x-coordinate are normalized by the Hertz contact pressure P0 and the contact radius a, respectively. Shears are normalized by the maximum contact pressure P0 times the friction coefficient l. A fretting loop in mixed slip regime (i.e. partial then gross slip) will now be simulated and results analysed. Point A of the fretting

4. Application to fretting mode I This first fretting mode, which is also the most often studied starting from Cattaneo (1938) and Mindlin (1949), see also Hills et al. (1988), corresponds to an oscillatory tangential movement/ loading. Two elastic bodies are first brought into contact with a normal load W. A tangential displacement ux along the x direction is then applied. Accordingly Coulomb the tangential load is limited by the normal load times the friction coefficient, i.e. Fx 6 l  W. In the case of spheres, both authors demonstrated independently that the circular contact area is divided into an inner circular stick area and an outer annular slip area. The algorithm presented above is now used to get a numerical solution. A rigid sphere of radius 105 lm is brought into contact on a flat composite where in what follows the Young’s modulus and Poisson’s ratio of the matrix are chosen as Em = 210 GPa and mm = 0.3, respectively. The maximal normal load Wmax = 80 mN, is applied gradually every mN. For the

Fig. 4. Loading path to simulate an entire fretting loop. The normal load W is applied first and kept equal to the maximum load Wmax when reached. A tangential displacement along x is then applied and cycling between ux = +0.04 lm and 0.04 lm is imposed to reproduce an entire fretting loop.

x

R z

Ha lf-s pac e

y

0.5r

r Cylindrical fibers

(a)

r = 0.13a

r 3r

(b)

Fig. 5. Schematic view of the contact problem with a rigid sphere pressed against a half-space containing cylindrical fibers. (a) ux is the amplitude of the tangential displacement along the x axis. (b) Size and location of the cylindrical fibers within the homogeneous and isotropic matrix.

3516

J. Leroux, D. Nélias / International Journal of Solids and Structures 48 (2011) 3510–3520

(a)

(b)

Fig. 6. Dimensionless numerical contact pressure and shears. Results are plotted along the x-axis in (a) and the y-axis in (b) for the loading points A, B, C, and D of the fretting loop (represented by the vertical dash line on the top right figure incrusted in each figure, see also Fig. 4). The plots in dash line correspond to the homogeneous solution (i.e. without fibers). When Qy is not indicated, its value is nil.

loop corresponds to the increment when the maximal normal load Wmax is reached before the tangential displacement starts.

Considering coupling effects for dissimilar materials in contact all influence coefficients linking the local traction and the elastic

J. Leroux, D. Nélias / International Journal of Solids and Structures 48 (2011) 3510–3520

deflection in all directions are considered. It reveals that the Hertz solution is no more accurate when the elastic properties of the two contacting bodies differ. Friction acts to increase the maximum contact pressure while decreasing the normal rigid displacement. In presence of cylindrical fibers the contact pressure exhibits peaks of magnitude higher than the homogeneous solution and the contact area slightly decreases. An increase of more than 20% of the maximum contact pressure is found when the fibers are twice stiffer than the half-space. Note that conversely to sinusoidal roughness that adds a sinusoidal perturbation of the Hertz contact pressure, the presence of stiff fibers increases mostly the contact pressure while decreasing the contact area whereas that of soft fibers decreases mostly the pressure compared to the Hertz solution while slightly increasing the contact are to maintain the normal identical. In other words, the presence of a stiff inclusion close to the surface produces mostly a local rise of the contact pressure (see Fig. 6a), and those of a soft one mostly a local collapse of the contact pressure, whereas the presence of roughness implies at the same time peaks of pressure on the top of asperities and drop of pressure where valleys are located. The same trends have been discussed by Leroux et al. (2010) for spherical inclusions (see Fig. 6). Loading points B and D of the fretting loop are taken in partial slip, i.e. when the tangential load is lower than the friction coefficient times the normal load. Globally the contact is sticking, but a slipping annulus appears at the periphery of the contact area. In the sliding area the dimensionless shear matches exactly the dimensionless contact pressure, which means that there are no local sticking zones. In the central stick zone the dimensionless shear is still perturbed by the presence of fibers. The loading point C is taken during gross slip. The dimensionless contact pressure and shear along the x-axis are identical and highly perturbed by the presence of stiff fibers close to the contact surface. For fretting problems when wear is considered, it is clear that the changes in the contact pressure and shear distributions induced by the presence of heterogeneities close to the contact surface could significantly modify the wear progress. Fretting wear predictions for a cylinder on a flat surface and under gross slip regime have been highlighted by Gallego and Nélias (2007). The prediction of the local wear within a sliding contact requires the use (and before that the identification) of a wear coefficient and the knowledge of the energy dissipated. The energy dissipated by friction depends on the slip amplitude and the surface shear stress.

3517

Thus the disturbed elastic tangential behavior could also affect the wear progress and will now be investigated. More details could be seen in an animation of the dimensionless contact pressure and shears along the x-axis during the entire fretting loop (Leroux and Nélias (2011a) Movie 1). The slip distribution for the loading point B of the fretting loop is plotted along the x-axis in Fig. 7(a) and the y-axis in Fig. 7(b). For comparison the slip distribution for the homogeneous half-space is also plotted with dotted lines. The slip sx along the x-axis is slightly amplified due to the presence of stiff cylindrical fibers while its distribution along the y-axis is slightly reduced. Note that the

(a)

(b)

Fig. 7. Numerical slips sx and sy. Results are plotted along the x-axis in (a) and y-axis in (b). Results correspond to slips for the loading point B introduced in Fig. 6. The dash line corresponds to the solution without fibers. When sy is not indicated, its value is nil.

Fig. 8. Fretting loop. Results for the homogeneous case and the simulation with cylindrical fibers are superimposed.

3518

J. Leroux, D. Nélias / International Journal of Solids and Structures 48 (2011) 3510–3520

Fig. 9. Distribution of dimensionless shear qx on the contact surface (a) with cylindrical fibers and (b) for a homogeneous half-space. Results are given for different loading points of the fretting loop.

transverse slip, due to the fact that the materials in contact are elastically dissimilar, leads to a slip sy that could not be neglected anymore (see Fig. 7(b)).

The fretting loop can be obtained by plotting the tangential load Fx versus the rigid body displacement dx, cf. Fig. 8. First it can be observed that the starting point and the final point of the tangential

3519

J. Leroux, D. Nélias / International Journal of Solids and Structures 48 (2011) 3510–3520

where

x3

 1 1 n 2 R  c21 c2 c3 logðR þ c1 Þ HðcÞ ¼ c1 c2 c3 R þ 4 o 6   þ R2  c22 c3 c1 logðR þ c2 Þ þ R2  c23 c1 c2 logðR þ c3 Þ    

1 c2 c 3 c3 c 1 c1 c 2  c41 arctan þ c42 arctan þ c43 arctan 12 c1 R c2 R c3 R

2a3

 EðcÞ ¼ c1 c2 logðR þ c3 Þ þ c2 c3 logðR þ c1 Þ þ c3 c1 logðR þ c2 Þ    

1 2 c 2 c3 c3 c1 c1 c 2  c1 arctan þ c22 arctan þ c23 arctan ðA2Þ 2 c1 R c2 R c3 R  1 with c ¼ ðc1 ; c2 ; c3 Þ; R ¼ c21 þ c22 þ c23 2 and cn defined as:

x1 2a

2

O

c1 ¼ ðx1  a1 ;x2  a2 ;x3  a3 Þ; c2 ¼ ðx1 þ a1 ;x2  a2 ;x3  a3 Þ

2a1

x2

c3 ¼ ðx1 þ a1 ;x2 þ a2 ;x3  a3 Þ; c4 ¼ ðx1  a1 ;x2 þ a2 ;x3  a3 Þ c5 ¼ ðx1  a1 ;x2 þ a2 ;x3 þ a3 Þ; c6 ¼ ðx1  a1 ;x2  a2 ;x3 þ a3 Þ

Fig. 10. Right parallelepiped inclusion.

ðA3Þ

c7 ¼ ðx1 þ a1 ;x2  a2 ;x3 þ a3 Þ; c8 ¼ ðx1 þ a1 ;x2 þ a2 ;x3 þ a3 Þ loading cycle are not superimposed on the fretting loop, outlining a slip ratcheting behavior. More surprisingly it is found that the fretting loop is almost not disturbed by the presence of the cylindrical fibers. Lastly it is also interesting to represent the dimensionless shear qx on the contact surface compared to the results for a homogeneous half-space, cf. Fig. 9. In addition to these results, the animation of the dimensionless shear qx on the contact surface during the entire fretting loop (Leroux and Nélias (2011b) Movie 2) shows that the distribution of the dimensionless shear magnitude is clearly disturbed along the fiber axis. This should also play an important role in the prediction of wear.

Appendix B. Elastic displacement at the interface subject to normal pressure (Kp) Consider the contact between a sphere and an elastic half-space with elastic constants (E1, m1) and (E2, m2), respectively. The contact surface x3 = 0 is discretized into rectangular surface areas 2D1  2D2. The origin of the Cartesian coordinate system (x1, x2, x3) is set to be the initial contact point. The elastic displacement at an observation point P(n  1, n2, 0) is related to the pressure field at the center Q ¼ n01 ; n02 ; 0 by the function Kp below:

Kp1 ðc1 ; c2 Þ ¼



4 1  m 21 1  m 22 X þ Kx ðc ; c Þ pE1 pE2 i¼1 i 1 2

5. Conclusion

with

This paper presents a semi-analytical method based on the Equivalent Inclusion Method to model three-dimensional fretting contact problem with composite materials. A first attempt is made here with a rigid sphere pressed against an elastic half-space containing long cylindrical fibers. The efficiency of the method was investigated through an example corresponding to fretting mode I (tangential displacement). Partial and gross slip fretting regimes have been simulated in a single fretting loop (mixed slip regime). It is found that cylindrical fibers when located close to the surface – or by extension any type of heterogeneities present nearby the contact surface – affect very significantly the distribution of surface stresses. Interestingly the change in the fretting loop, which is at a larger scale of observation, appears almost negligible. Next step will be to investigate the effect of the presence of heterogeneous inclusions on wear.

B K x1 ðc1 ;c2 Þ ¼ 2ðc1 þ D1 Þarctan B @

Appendix A. Harmonic U and biharmonic w potentials for a cuboidal domain

wðxÞ ¼

8 P

ð1Þn Hðcn Þ

n¼1

/ðxÞ ¼

8 P

ð1Þn Eðcn Þ

n¼1

ðA1Þ

ðc2 þD2 Þþ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ðc2 þD2 Þ þðc 1 þD1 Þ ðc1 þD1 Þ

1

C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C 2 2 A þðc2 þ D2 Þlog ðc2 þ D2 Þ þ ðc1 þ D1 Þ 0 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2

B K x2 ðc1 ;c2 Þ ¼ 2ðc1  D1 Þarctan B @

ðc2 D2 Þþ

ðc1 D1 Þ þðc 2 D2 Þ ðc1 D1 Þ

C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C A þðc2  D2 Þlog ðc2  D2 Þ2 þ ðc1  D1 Þ2 0 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi B

ðc2 D2 Þþ

ðc2 D2 Þ2 þðc1 þD1 Þ2 ðc 1 þD1 Þ

ðc2 þD2 Þþ

ðc2 þD2 Þ2 þðc1 D1 Þ2 ðc 1 D1 Þ

C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C 2 2 A þðc2  D2 Þlog ðc2  D2 Þ þ ðc1 þ D1 Þ 0 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

K x3 ðc1 ;c2 Þ ¼ 2ðc1 þ D1 Þarctan B @

B

K x4 ðc1 ;c2 Þ ¼ 2ðc1  D1 Þarctan B @

C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C 2 2 A þðc2 þ D2 Þlog ðc2 þ D2 Þ þ ðc1  D1 Þ ðB2Þ

n01

n02 .

where, c1 ¼ n1  and c2 ¼ n2  The influence coefficient K p2 ðc1 ; c2 Þ is obtained from the influence coefficient K p1 ðc1 ; c2 Þ by permutation of the variables c1 and c2.

Kp1 ðc1 ; c2 Þ ¼ MacMillan (1958) derived the potential functions of a homogenous rectangular parallelepiped of size 2a1  2a2  2a3. Let the origin of a rectangular coordinate system be taken with the origin at the center of the parallelepiped and the axes parallel to its edges (Fig. 10). In this way the complete potential functions of the homogeneous right parallelopiped are found to be

0

ðB1Þ



4 1  m 21 1  m 22 X þ Kx ðc ; c Þ pE1 pE2 i¼1 i 1 2

with

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 2 2 þD2 Þ þðc 1 þD1 Þ pðc ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðc2 D2 Þþ ðc2 D2 Þ2 þðc1 þD1 Þ2  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðc þD Þþ ðc2 þD2 Þ2 þðc1 þD1 Þ2 ffi K n2 ðc1 ; c2 Þ ¼ ðc2 þ D2 Þ log 1 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðc1 D1 Þþ ðc2 þD2 Þ2 þðc1 D1 Þ2  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðc D Þþ ðc2 D2 Þ2 þðc1 D1 Þ2 ffi K n3 ðc1 ; c2 Þ ¼ ðc1  D1 Þ log 2 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðc2 þD2 Þþ ðc2 þD2 Þ2 þðc1 D1 Þ2  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðc D Þþ ðc2 D2 Þ2 þðc1 D1 Þ2 ffi K n4 ðc1 ; c2 Þ ¼ ðc2  D2 Þ log 1 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2

K n1 ðc1 ; c2 Þ ¼ ðc1 þ D1 Þ log



ðc2 þD2 Þþ

ðc1 þD1 Þþ

where, c1 ¼ n1 

n01

ðB1Þ

and c2 ¼ n2 

ðc2 D2 Þ þðc1 þD1 Þ

n02 .

ðB4Þ

3520

J. Leroux, D. Nélias / International Journal of Solids and Structures 48 (2011) 3510–3520

Appendix C. Supplementary material Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.ijsolstr.2011.09.007. References Boucly, V., Nélias, D., Liu, S., Wang, Q.J., Keer, L.M., 2005. Contact analyses for bodies with frictional heating and plastic behavior. ASME Journal of Applied Mechanics 127 (2), 355–364. Boucly, V., Nélias, D., Green, I., 2007. Modeling of the rolling and sliding contact between two asperities. ASME Journal of Applied Mechanics 129 (2), 235–245. Cattaneo, C., 1938. Sul contatto di due corpi elastici: distribuzione locale degli sforzi. Rendiconti dell Accademia Nazionale dei Lincei 27 (6), 342–348, 434– 436, 474–478. Eshelby, J.D., 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.D., 1959. The elastic field outside an elastic inclusion. Proceedings of the Royal Society of London A252, 561–569. Eshelby, J.D., 1961. Elastic inclusions and inhomogeneities. In: Sneddon, L.N., Hill, R. (Eds.), Progress in Solid Mechanics, vol. 2. North-Holland, Amsterdam, pp. 9– 140. Fotiu, P.A., Nemat-Nasser, S., 1996. A universal integration algorithm for ratedependant elastoplasticity. Computers and Structures 59 (6), 1173–1184. Fulleringer, B., Nélias, D., 2010. On the tangential displacement of a surface point due to a cuboid of uniform plastic strain in a half-space. ASME Journal of Applied Mechanics 77 (2), 21014-1_021014-7. Gallego, L., Nélias, D., Jacq, C., 2006. A comprehensive method to predict wear and to define the optimum geometry of fretting surfaces. ASME Journal of Tribology 128 (3), 476–485. Gallego, L., Nélias, D., 2007. Modeling of fretting wear under gross slip and partial slip conditions. ASME Journal of Tribology 129 (3), 528–535. Gallego, L., Fulleringer, B., Deyber, S., Nélias, D., 2010a. Multiscale computation of fretting wear at the blade/disk interface. Tribology International 43, 708–718. Gallego, L., Nélias, 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. Hills, D.A., Nowell, D., O’Connor, J.J., 1988. On the mechanics of fretting fatigue. Wear 125 (1–2), 129–146. Jacq, C., Nélias, D., Lormand, G., Girodin, D., 2002. Development of a threedimensional semi analytical elastic–plastic contact code. ASME Journal of Tribology 124 (4), 653–667. Johnson, K.L., 1985. Contact Mechanics. Cambridge University Press, London, UK. Kuo, C.H., 2007. Stress disturbances caused by the inhomogeneity in an elastic halfplane subjected to contact loading. International Journal of Solids and Structures 44, 860–873. Kuo, C.H., 2008. Contact stress analysis of an elastic half-plane containing multiple inclusions. International Journal of Solids and Structures 45, 4562–4573. Leroux, J., Fulleringer, B., Nélias, D., 2010. Contact analysis in presence of spherical inhomogeneities within a half-space. International Journal of Solids and Structures 47 (22–23), 3034–3049. Leroux, J., Nélias, D., 2011a. Movie 1: dimensionless contact pressure and shears along the x-axis during fretting loop. Leroux, J., Nélias, D., 2011b. Movie 2: dimensionless shear qx on the contact surface during fretting loop. Love, A.E.H., 1927. A Treatise on the Mathematical Theory of the Elasticity, fourth ed. Cambridge University Press, Cambridge, UK. Liu, S.B., Wang, Q., Liu, G., 2000. A versatile method of discrete convolution and FFT (DC-FFT) for contact analyses. Wear 243, 101–111.

Liu, S.B., Wang, Q.J., 2001. A three-dimensional thermomechanical model of contact between nonconforming rough surfaces. ASME Journal of Tribology 123, 17–26. Liu, S., Wang, Q., 2005. Elastic fields due to eigenstrains in a half-space. ASME Journal of Applied Mechanics 72, 871–878. MacMillan, W.D., 1958. The Theory of the Potential. Dover, New York. Mindlin, R.D., 1949. Compliance of elastic bodies in contact. ASME Journal of Applied Mechanics 16, 259–268. Mohrbacher, H., Celis, J.P., Roos, J.R., 1995. Laboratory testing of displacement and load induced fretting. Tribology International 28, 269–278. Moschovidis, Z.A., 1975. Two ellipsoidal inhomogeneities and related problems treated by the equivalent inclusion method. Ph.D. thesis, Northwestern University, Illinois. Moschovidis, Z.A., Mura, T., 1975. Two-ellipsoidal inhomogeneities by the equivalent inclusion method. ASME Journal of Applied Mechanics 42, 847–852. Mura, T., Cheng, D.C., 1977. The elastic field outside an ellipsoidal inclusion. ASME Journal of Applied Mechanics 44, 591–594. Mura, T., Seo, K., 1979. The elastic field in a half space due to ellipsoidal inclusions with uniform dilatational eigenstrains. ASME Journal of Applied Mechanics 46, 568–572. Mura, T., 1987. Micromechanics of Defects in Solids, second ed. Martinus Nijhoff, Dordrecht. Mura, T., 1988. Inclusion problem. ASME Journal of Applied Mechanics Review 41, 15–20. Nélias, D., Boucly, V., Brunet, M., 2006. Elastic–plastic contact between rough surfaces: proposal for a wear or running-in model. ASME Journal of Tribology 128 (2), 236–244. Ortiz, M., Simo, J.C., 1986. An analysis of a new class of integration algorithms for elasto-plastic constitutive relations. International Journal for Numerical Methods in Engineering 23, 353–366. Panagiotopoulos, P.D., 1975. A nonlinear programming approach to the unilateral contact and friction boundary value problem in the theory of elasticity. Inzenyr Architekt 44 (6), 421–432. Polonsky, I.A., Keer, L.M., 1999. A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient method. Wear 231, 206–219. Rodin, G.J., Hwang, Y.L., 1991. On the problem of linear elasticity for an infinite region containing a finite number of nonintersecting spherical inhomogeneities. International Journal of Solids and Structures 27 (2), 145–159. Ruiz, C., Wang, Z.P., Webb, P.H., 1992. Techniques for the characterization of fretting fatigue damage. ASTM STP 1159, 170–177. Sendeckyj, G.P., 1967. Ellipsoidal inhomogeneity problem. Ph.D. dissertation, Northwestern University, Evanston, IL. Simo, J.C., Taylor, R.L., 1986. A return mapping algorithm for plane stress elastoplasticity. International Journal for Numerical Methods in Engineering 22, 649–670. Vincent, L., Berthier, Y., Godet, M., 1992. Testing methods in fretting fatigue: a critical appraisal. ASTM STP 1159, 33–48. Wang, F., Keer, L.M., 2005. Numerical simulation for three dimensional elasticplastic contact with hardening behavior. ASME Journal of Tribology 127, 494–502. Xu, Y., Zhang, Y., Cheng, L., 2007. Preparation and friction behavior of carbon fiber reinforced silicon carbide matrix composites. Ceramics international 33, 439– 445. Zhou, K., Chen, W.W., Keer, L.M., Wang, Q.J., 2009. A fast method for solving threedimensional arbitrarily shaped inclusions in a half space. Computer Methods in Applied Mechanics and Engineering 198 (9–12), 885–892. Zhou, K., Chen, W.W., Keer, L.M., Ai, X., Sawamiphakdi, K., Glaws, P., Wang, Q.J., 2011. Multiple 3D inhomogeneous inclusions in a half space under contact loading. Mechanics of Materials 43 (8), 444–457. Zhou, Z.R., Vincent, M., 1993. Effect of external loading on wear maps of aluminium alloys. Wear, 619–623.