sliding of a spherical indenter upon a viscoelastic half-space containing an ellipsoidal inhomogeneity

Journal of the Mechanics and Physics of Solids 80 (2015) 1–25

Contents lists available at ScienceDirect

Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps

Rolling contact of a rigid sphere/sliding of a spherical indenter upon a viscoelastic half-space containing an ellipsoidal inhomogeneity Koffi Espoir Koumi a,b, Thibaut Chaise a, Daniel Nelias a,n a b

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

a r t i c l e i n f o

abstract

Article history: Received 13 August 2014 Received in revised form 12 February 2015 Accepted 2 April 2015 Available online 10 April 2015

In this paper, the frictionless rolling contact problem between a rigid sphere and a viscoelastic half-space containing one elastic inhomogeneity is solved. The problem is equivalent to the frictionless sliding of a spherical tip over a viscoelastic body. The inhomogeneity may be of spherical or ellipsoidal shape, the later being of any orientation relatively to the contact surface. The model presented here is three dimensional and based on semi-analytical methods. In order to take into account the viscoelastic aspect of the problem, contact equations are discretized in the spatial and temporal dimensions. The frictionless rolling of the sphere, assumed rigid here for the sake of simplicity, is taken into account by translating the subsurface viscoelastic fields related to the contact problem. Eshelby's formalism is applied at each step of the temporal discretization to account for the effect of the inhomogeneity on the contact pressure distribution, subsurface stresses, rolling friction and the resulting torque. A Conjugate Gradient Method and the Fast Fourier Transforms are used to reduce the computation cost. The model is validated by a finite element model of a rigid sphere rolling upon a homogeneous vciscoelastic half-space, as well as through comparison with reference solutions from the literature. A parametric analysis of the effect of elastic properties and geometrical features of the inhomogeneity is performed. Transient and steadystate solutions are obtained. Numerical results about the contact pressure distribution, the deformed surface geometry, the apparent friction coefficient as well as subsurface stresses are presented, with or without heterogeneous inclusion. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Contact Mechanics Semi-analytical methods (SAM) Viscoelasticity Rolling torque Apparent friction coefficient Eshelby's equivalent inclusion method (EIM) Inhomogeneity Eigenstrain

1. Introduction Polymer based materials are extensively used in several engineering domains due to their numerous advantages such as low cost of raw materials and less complex manufacturing, low weight, and compatibility with most liquids and lubricants. They are also sometimes biocompatible or biodegradable. In order to increase some of their mechanical properties these materials can be reinforced by adding small particles or fibers. Since Hunter (1961), it is well known that (frictionless) rolling over a viscoelastic material induces an apparent friction coefficient. This apparent friction has been studied by several authors in the case of homogeneous viscoelastic solid, mostly in the steady-state regime.

n

Corresponding author. Fax: þ 33 4 72 43 89 13. E-mail address: [email protected] (D. Nelias).

http://dx.doi.org/10.1016/j.jmps.2015.04.001 0022-5096/& 2015 Elsevier Ltd. All rights reserved.

2

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

center

Nomenclature Letters contact radius an a1, a2, a3 semi-axes of an ellipsoidal inhomogeneity influence coefficients that relating the stress Bijkl sij at point (x1, x2, x 3 ) to the constant eigenstrain at the point (x1k , x2k , x 3k ) M I Cijkl, Cijkl elastic constants of the matrix and the inhomogeneity EI Young's modulus of the inhomogeneity h distance between the two surfaces of the contacting bodies H(t) the Heaviside step function Iijkl the fourth-order identity tensor J(t) viscoelastic creep function Kn coefficients in the normal displacement at the contact surface due to the contact pressure Mij influence coefficients relating the stress sij at the point (x1, x2, x 3 ) to the normal traction sn within a discretized area centered at (x1k , x2k , 0) n1, n2, n3 grid sizes in the half-space along the Cartesian directions x1, x2, x3, respectively P0 maximum Hertzian pressure p contact pressure distribution D indenter diameter R(t) viscoelastic relaxation function Sijkl components of Eshelby's tensor 0 ui displacements corresponding to the infinite 0 applied strain εij ui disturbed contribution of the displacements W applied exterior load dx3 depth of the inclusion from the surface of the matrix in EF model x I = (x1I , x2I , x 3I ) Cartesian coordinates of the inclusion

Greek letters

εij εij

0

εij⁎ 0

sij

sij ϕ, Ψ

δij sn

Δx1, Δx2 νM, νI

γ η τ θ

infinite applied strain strain due to eigenstrains eigenstrain due to the presence of inhomogeneities stress corresponding to the infinite applied 0 strain εij disturbed contribution of the stresses harmonic and biharmonic potentials of mass density εij⁎ Kronecker symbol normal pressure due to the summation of both symmetric inclusions half-size of the discretized surface area Poisson's ratio of the matrix M and the inclusion I the ratio of the inhomogeneity Young's modulus to the matrix the dashpot viscosity the relaxation time the tilted angle of the inhomogeneity in the x1Ox 3 plane

Acronyms and fast Fourier transforms 2D-FFT 3D-FFT FFT  1 Blijkl

mij M

two-dimensional fast Fourier transform three-dimensional fast Fourier transform inverse FFT operation frequency response of coefficients Bijkl in the frequency domain frequency response of coefficients Mij in the frequency domain

The first works on viscoelastic contact focused on the indentation problem between a rigid indenter and a viscoelastic solid. Lee and Radok (1960) obtained the contact pressure distribution for the spherical indentation of a linear viscoelastic material for a monotonic increase of the contact area. Their model has been extended by Hunter (1960) and Graham (1967) to the indentation of viscoelastic materials when the contact radius possesses a single maximum. More recently Greenwood (2010) introduced a model to solve the contact problem between an axisymmetric indenter and a viscoelastic half-space. The analytical solution of the rebound indentation problem for a linear viscoelastic layer has been given by Argatov (2012). One of the main limitations of all the approaches mentioned above is their restriction to ideal linear viscoelastic materials with one relaxation time. Chen et al. (2011) recently introduced a robust semi-analytical approach to solve indentation problems between a rigid indenter and a homogeneous viscoelastic half-space. The semi-analytical approach allows to account for a wide spectra of relaxation times for linear viscoelastic materials, an arbitrary loading profile and can also be used to simulate the contact between rough surfaces. The semi-analytical approach was recently extended to solve the indentation problem between a rigid indenter and a heterogeneous viscoelastic material (Koumi et al., 2014a). The model can account for the presence of isotropic or anisotropic elastic ellipsoidal inhomogeneities of any orientation within the viscoelastic matrix. Other authors investigated the steady-state response of the frictionless rolling contact problem when one of the bodies in contact is homogeneous and viscoelastic. The two-dimensional contact problem of a rigid cylinder rolling over a viscoelastic halfspace was first solved by Hunter (1961). His plane strain model was limited to an ideal viscoelastic material with one relaxation time. Later, Panek and Kalker (1980) extended Hunter's approach to three-dimensional problems by using an approximation based on the elastic line contact theory (Kalker, 1972, 1977; Panek and Kalker, 1977). This assumption represents a very strong approximation. Goriacheva (1973) used another approach to solve the rolling contact problem between a cylindrical rigid indenter and viscoelastic halfspace. The three-dimensional sliding contact problem of a smooth indenter and a viscoelastic halfspace has been solved later by Aleksandrov et al. (2010). The contact pressure distribution and the resulting torque are presented in the case of an ideal viscoelastic material with one relaxation time. Persson (2010) presented a new analytical theory for the rolling contact

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

3

of a hard cylinder or a sphere over a homogeneous viscoelastic solid. This approach assumes that, when the rolling speed increases, the viscoelastic material stiffens, but the pressure distribution keeps on being of a Hertzian type. Although Persson's theory provides a good estimate of the apparent friction coefficient, his model is not able to predict the contact pressure distribution for a force driven problem. Nowadays many researchers are employing the finite element (FE) method (Padovan and Paramadilok, 1984; Padovan et al., 1992; Le Tallec and Rahler, 1994; Nasdala et al., 1998; Nackenhorst, 2004). The FE method can account for any viscoelastic material and any geometry of the contacting bodies. For three-dimensional problems and a moving load the computation cost is very high, almost unaffordable when a fine mesh is required. Besides FE model users should be very careful with the way the contact is solved since a correct contact pressure distribution is one of the most difficult outputs of the FE method. Recently, Carbone and Putignano (2013) introduced a novel methodology to investigate steady-state viscoelastic sliding or rolling contacts based on the boundary element method. This approach is able to deal with large spectra of relaxation times for the viscoelastic materials. Steady-state solutions in terms of apparent friction coefficient and contact pressure distribution have been reached for homogeneous viscoelastic bodies. This paper proposes a new approach, based on semi-analytical methods, to solve the three-dimensional rolling/sliding contact problem with heterogeneous viscoelastic bodies. The model permits to get the solution in terms of contact pressure distribution, subsurface stresses, apparent friction coefficient for frictionless or frictional rolling/sliding contact problems, whatever the geometry is, both in the transient and then steady-state regimes. The analysis will be limited here for the sake of simplicity to the frictionless contact between a spherical and rigid indenter rolling and/or sliding (since frictionless) over a heterogeneous viscoelastic solid containing mostly one elastic isotropic inhomogeneity. It should be noted that the model may account for any number of relaxation times for the viscoelastic behavior, surface roughness if needed, and may also consider many either isotropic or anisotropic inhomogeneities of any shape including when lying at the surface. The model developed in this paper is originally intended to composite materials. Some of these composite materials are made of a viscoelastic matrix and often reinforced with elastic fibers. Besides the apparent friction coefficient, the knowledge of the contact pressure distribution and subsurface stresses is also of main importance. The contact pressure distribution will be seen by the counterface, sometimes metallic such as in some aeronautical applications (composite blade against steel disk for example). Subsurface stresses are quite important when developing a decohesion or damage model, to better understand the damage mechanisms in composite materials when two surfaces are rubbing. The multiple phenomena involved being quite complex, the first step is to look at the effect of one single elastic inhomogeneity embedded into a viscoelastic matrix. The effect of a few interacting heterogeneities will be also investigated.

2. Theoretical background and model description 2.1. Contact between a rigid indenter and a homogeneous viscoelastic material The theory presented below is limited to the framework of linear viscoelasticity. Based on the generalized Maxwell model (see Fig. 1), the relaxation function (for t ≥ 0) can be written as follows:

R (t) =

⎡ σ (t) = ⎢μ 0 + ε 0 H (t) ⎣⎢

n

⎤

i=0

⎥⎦

∑ μi exp( − t /τi ) ⎥ H (t)

(1)

where ε0 H (t) is the product of the constant applied strain and the unit step function, μi the spring stiffness, ηi the dashpot viscosity, and τi = ηi /μi the relaxation time of one elementary model. Eq. (1) is also known as the Prony series decomposition. Likewise, the creep compliance function J(t) (for t ≥ 0) is defined as the ratio of the strain history over the product of the constant applied stress and the unit step function:

Fig. 1. Generalized Maxwell model.

4

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

Fig. 2. Rolling of a sphere or sliding of a spherical indenter over a viscoelastic half-space.

Fig. 3. Schematic view of the deformed geometry (the vertical deformations are amplified) and contact pressure distribution (dash lines) for asymptotic solutions of the frictionless viscoelastic (VE) contact problem between a sphere and a flat with the same imposed normal displacement: rigid sphere against viscoelastic flat (left column), viscoelastic sphere against rigid flat (right column). Comparison between indentation (top line), pure sliding (middle line) and pure rolling (bottom line). F is the normal load, T the tangential (resisting) force, and M the rolling (resisting) moment. V is the velocity of the sphere center and ω the rotational speed of the sphere when rolling. δ is the shift of the center of pressure from the sphere center, and R the sphere radius. Note that the pressure is always normal to the surface, i.e. towards the sphere center when rigid (left column) or vertical when the flat is rigid (right column).

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

J (t) =

ε (t) σ0 H (t)

5

(2)

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

∫0

t

J (ξ ) R (t − ξ ) dξ = t

(3)

A geometrical description of a rolling/sliding sphere or a sliding (spherical) indenter over a viscoelastic body is presented in Fig. 2. Six different configurations can be found when one of the bodies in contact behaves viscoelastically whereas the other one is rigid, depending on the kinematics, see Fig. 3. In pure rolling the angular velocity of the rolling sphere is defined as ω. A frictionless contact is considered. The rolling problem of a rigid spherical indenter over a viscoelastic half-space is equivalent to simple sliding of spherical indenter at velocity:

(4)

v=R×ω

where v is the sliding velocity and R the sphere or indenter tip radius. For the sake of brevity, this paper focuses on the rolling of a sphere or sliding of an indenter moving at a constant velocity. It is however important to note that the model may account for any kind of indenter's geometries moving at non-constant velocity. The reference frames used in the modeling are defined as follows: R0 (O0, X10 , X20 , X30 ) corresponds to the initial contact set of axes at the initial time t¼ 0. For t ≠ 0 the current contact axes are identified by R (O, X1, X2 , X3 ), O being the projection of the indenter's center on the matted surface. R (O, X1, X2 , X3 ) follows the indenter movement. The contact equations are given and solved in this reference frame R (O, X1, X2 , X3 ). The total time of simulation is divided into Nt uniform time steps Δt . The sphere or indenter movement at each time step is equal to

→ → Δ r = Δt . v

(5)

⎯⎯⎯→ ⎯⎯⎯→ ⎯⎯⎯→ X10 , X20 , X30 are the unit vectors along the (OX10 ) , (OX20 ) , (OX30 ) axes, respectively. A rolling contact problem between two contacting bodies can be described by a set of equations that must be solved simultaneously. These equations are:



Load balance: The applied load W(t) and the integration of the contact pressure p (x1, x2, t) in the real contact area Γc (t) at time t must be strictly equal:

W (t) =

∫Γ (t) p (x1, x2, t) dx1 dx2

(6)

c

 Surface separation: The gap between the two contacting surfaces in the moving reference frame linked to the contact h (x1, x2, t) is defined by the initial geometry hi (x1, x2 ), the rigid body displacement δ and by the normal surface displacement of the viscoelastic half-space u3 (x1, x2, t)

h (x1, x2 , t) = hi (x1, x2 ) + δ (t) + u3 (x1, x2 , t)

(7)

where hi (x1, x2 ) is the initial geometry, and δ the rigid body displacement. The model is able to account for any kind of sliding indenter geometry. The information about the indenter geometry is contained in hi (x1, x2 ) which is the gap between the contacting surfaces. The surface normal displacement, u3 (x1, x2, t), at time t of a surface point of the viscoelastic half-space can be expressed as

u3 (x1, x2 , t) =

∞

∞

∫−∞ ∫−∞ ∫0

t

G (x1 − x1′ , x2 − x2′ , t − ξ)

∂p (x1′ , x2′ , t) dx1′ x2′ dξ ∂ξ

(8)

where G (x1, x2, t) are the viscoelastic Green's functions. G (x1, x2, t) can be decomposed into space and time functions:

G (x1, x2 , t) =

Therefore it yields

(1 − ν) π x12 + x22

J (t) = G (x1, x2 ) J (t) (9)

6

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

u3 (x1, x2 , αΔt) =

∞

∞

∫−∞ ∫−∞ G (x1 − x1′, x2 − x2′ ) dx1′ x2′

α

∑ J [(α − k) Δt] k=0

→ ⎯⎯⎯→ → ⎯⎯⎯→ × [p (x1′0 − α·Δ r · X10 , x2′0 − α·Δ r · X20 , k) → ⎯⎯⎯→ → ⎯⎯⎯→ − p (x1′0 − α·Δ r · X10 , x2′0 − α·Δ r · X20 , k − 1)]

(10)

 The contact conditions: The distance h (x1, x2, t) should not be negative, because the contacting bodies cannot interpenetrate each other. These conditions are defined by the following inequalities:

contact:

h (x1, x2 , t) = 0

separation:

and

h (x1, x2 , t) > 0

p (x1, x2 , t) > 0 in Γc (t)

and

p (x1, x2 , t) = 0 out Γc (t)

(11)

These equations are solved simultaneously using the Conjugate Gradient Method (CGM). 2.2. Theoretical background – Eshelby's equivalent inclusion method in Contact Mechanics The presence of a single or multiple inhomogeneities within one of the bodies in contact (see Figs. 4 and 5) is taken into account by adding in Eq. (7) the eigendisplacement u3⁎ induced by these inhomogeneities. Eq. (7) is then modified as follows:

h (x1, x2 , t) = hi (x1, x2 ) + δ (t) + u3 (x1, x2 , t) + u3* (x1, x2 , t)

(12)

The algorithm presented in Fig. 6 is used to couple the contact problem with the subsurface inhomogeneous problem through the additional displacement field u3⁎. The bases of the method in the case of an elastic or elastic–plastic half-space are detailed extensively in Jacq et al. (2002), Nelias et al. (2006), Chen et al. (2008a), Leroux et al. (2010), Leroux and Nelias

Fig. 4. 3D view of a moving sphere over a viscoelastic half-space containing one ellipsoidal inclusion, the sphere is moving along the X1-axis.

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

7

Fig. 5. 2D view of a moving sphere over a viscoelastic half-space containing one ellipsoidal inclusion, the sphere is moving along the X1-axis.

(2011), and Koumi et al. (2014b) and for a viscoelastic half-space (Chen et al., 2011; Koumi et al., 2014a). Note that the problem with a moving load is quite similar to the elastic–plastic one, except that the permanent print on the elastic–plastic surface (Nelias et al., 2007; Chen et al., 2008b; Chaise and Nelias, 2011) is recovered for a viscoelastic material. Sections 2.2.1–2.2.3 present the various steps for the calculation of u3⁎ in the case of a viscoelastic body containing elastic inhomogeneities. The calculation of the eigenstrain ε⁎ is presented first. The calculation of the subsurface stresses is then detailed. Finally the calculation of the eigendisplacement u3⁎ is reported. Eshelby's formalism, also known as the Equivalent Inclusion Method (EIM), is applied at every step of the temporal M discretization in order to take into account the presence of inhomogeneities. Cijkl (αΔt) represents the elastic stiffness tensor of the viscoelastic matrix at time t = αΔt .

2.2.1. Determination of the eigenstrain ε⁎ M Considering an infinite matrix D with the elastic stiffness tensor Cijkl (αΔt) containing an ellipsoidal domain Ω with the I elastic stiffness tensor Cijkl submitted at infinity to a uniform strain ε0. Using the EIM method, the eigenstrain induced by the presence of the inhomogeneity can be calculated by the following equation: M ⁎ ΔCijkl (αΔt) Sklmn (αΔt) εmn + Cijkl (αΔt) ε kl⁎ = − ΔCijkl (αΔt) ε kl0

(13)

8

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

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

where I M ΔCijkl (αΔt) = Cijkl − Cijkl (αΔt)

Sijkl (αΔt) is Eshelby's tensor at time t = αΔt . Eshelby's tensor depends solely on Poisson's ratio of the matrix. In the model, this Poisson's ratio is assumed to be constant. It can be shown that

Sijkl (αΔt) = Sijkl (t = 0) = Sijkl

(14)

2.2.2. Calculation of the subsurface stresses Applying Eshelby's EIM and Chiu's decomposition (Chiu, 1978) to the contact problem, the stress components at any point of the domain meshed with n1 × n2 × n3 cuboids are given by

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

9

n 3 − 1 n 2 − 1 n1 − 1

σ ij (x1, x 2, x 3, αΔt) =

∑ ∑ ∑ x 3I = 0 x 2I = 0 x1I = 0

Bijkl (x1 − x1I , x 2 − x 2I , x 3 − x 3I , αΔt) ε kl⁎ (x1I , x 2I , x 3I , αΔt)

n 3 − 1 n 2 − 1 n1 − 1

+

∑ ∑ ∑ x 3I = 0 x 2I = 0 x1I = 0

⁎ (x I , x I , − x I , αΔt) Bijkl (x1 − x1I , x 2 − x 2I , x 3 + x 3I , αΔt) ε skl 1 2 3

n 2 − 1 n1 − 1

−

Mij (x1 − x1I , x 2 − x 2I , x 3, αΔt) σ n (x1I , x 2I , 0, αΔt)

∑ ∑ x 2I = 0 x1I = 0

(15)

where Bijkl are the influence coefficients that relate the constant eigenstrain at the point (x1I , x2I , x 3I ) which is the inclusion center in an infinite space to the stress sij at the point (x1, x2, x 3 ). Mij represents the influence coefficients relating the normal traction sn within a discretized area centered at (x1I , x2I , 0) to the stress sij at the point (x1, x2, x 3 ): M Bijkl (x, αΔt) = Cijmn (αΔt) Dmnkl (x, αΔt)

for x in D − Ω

M Bijkl (x, αΔt) = Cijmn (αΔt)(Dmnkl (x, αΔt) − Imnkl ) 1 2

(

for x in Ω

(16) (17)

)

where Iijkl = δ il δ jk + δ ik δ jl is the fourth-order identity tensor. The second-order tensor Mij (see Appendix B) and the fourth-order tensor Dijkl depends only on Poisson's ratio which is assumed to be constant here for simplicity. It yields

Dijkl (x, αΔt) = Dijkl (x)

(18)

Mij (x, αΔt) = Mij (x)

(19)

and

The expression of Dijkl is given by (Mura, 1987)

Dijkl =

1 [Ψ, ijkl − 2νδ kl ϕ, ij − (1 − ν)(δ kl ϕil + δ ki ϕ, jl + δ jl ϕ, ik + δli ϕ, jk )] 8π (1 − ν)

Ψ (x) =

∫Ω |x − x′| dx′

ϕ (x) =

∫Ω

1 |x − x′|

(20)

dx′

For a single inclusion centered at (x1I , x2I , x 3I ) within the half-space, the normal traction sn at the surface point (x1′ , x2′ , 0) is derived:

σ n (x1′ , x2′ , 0, αΔt) = − B33kl (x1′ − x1I , x2′ − x2I , − x 3I , αΔt) ε kl⁎ (x1I , x2I , x 3I , αΔt) ⁎ (x I , x I , − x I , α Δt) − B33kl (x1′ − x1I , x2′ − x2I , x 3I , αΔt) ε skl 1 2 3

(21)

In Eq. (15), each component Mij () is obtained by a double integration of the functions Fij () over a discretized surface area 2Δx1 × 2Δx2 centered at (x1I , x2I , 0), see Appendices A and B:

Mij (x1 − x1I , x2 − x2I , x 3 ) =

x1I +Δx1

x 2I +Δx 2

I 1

I 2

∫x −Δx ∫x −Δx 1

2

Fij (x1 − x1′ , x2 − x2′ , x 3 ) dx1′ x2′

(22)

3D-FFT are employed to accelerate the calculation of the first and second terms in Eq. (15), and 2D-FFT for the third term. Wrap around order and zero-padding techniques are used in order to remove the induced periodicity error (Liu et al., 2000). 2.2.3. Calculation of the normal displacement u3⁎ The surface normal eigendisplacement u3⁎ can be obtained when inserting the eigenstrain into the total strain. They are generated by the pressure field sn only. The normal displacements are calculated as n2 − 1 n1 − 1

u3⁎ (x1, x2 , αΔt) =

∑ ∑ x 2′ = 0 x1′= 0

K n (x1 − x1′ , x2 − x2′ ) σ n (x1′ , x2′ , αΔt) (23)

To solve the equation above numerically, the surface in contact is discretized into n1 × n2 rectangular elements of uniform size 2Δx1 × 2Δx2. Then, pressure and displacement within each surface element are treated as constant and their values located at the center. The effect of a uniform pressure on a rectangular area has been given by Love (1952) and Johnson (1985). K n (x1 − x1′ , x2 − x2′ ) denotes the influence coefficients that relate the normal pressure at the surface point (x1′ , x2′ , 0) to the normal displacement at the surface point (x1, x2, 0) , recalled in Appendix C.

10

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

2.3. Determination of the apparent friction coefficient and the resulting torque The apparent friction coefficient and the resulting torque are the two aspects of the same phenomenon: the hysteresis losses in the viscoelastic material. Considering a rolling contact along the X1-axis, the apparent friction coefficient and the resulting torque can be obtained using Eqs. (24) and (25).

 Apparent friction coefficient: The tangential force induced by the hysteresis losses within a viscoelastic material is given by

FT (t) =



∫Γ (t) c

∂u3 (x1, x2 , t) p (x1, x2 , t) dx1 dx2 ∂x1

(24)

where Γc (t) is the real contact area at time t. The friction coefficient is then calculated as μapp = FT /FN , where FN is the normal load. The finite difference method is used to calculate the expression ∂u3 (x1, x2, t)/∂x1, and the rectangle rule is used to evaluate the value of the integral FT(t). Resulting torque: For a rolling sphere that creates a resulting torque, acting opposite to the rotational movement, given by

M (t) =

∫Γ (t) x1p (x1, x2, t) dx1x2 c

(25)

The rectangle rule is used to evaluate the value of the integral M(t). As mentioned above the resulting torque and the apparent friction coefficient are induced by the hysteresis losses of viscoelastic materials. At the macroscopic scale these two scalar variables are related by the equation below:

M (t) = μ app × F N × R

(26)

3. Validation The semi-analytical model is validated by two different approaches: firstly by comparison with a finite element model using the commercial FE package Abaqus v6.9 and secondly by comparison with the model developed by Carbone and Putignano (2013) based on the boundary element method.

 Validation by finite element analysis: The configuration is shown in Fig. 7. A rigid sphere and a viscoelastic solid are meshed. The sphere may translate and the contact is assumed frictionless. Only half of the bodies in contact are meshed since the problem is symmetric with respect to the plane Y = 0. A local area called ‘subsurface layer’ is meshed by 30,250 linear hexahedral elements of size 0.02a⁎ to describe the contact problem as accurately as possible, an being the contact radius of the equivalent Hertzian solution. The rest of the solid is meshed with 1,088,349 linear tetrahedral (C3D4) elements. The (quarter) sphere contains 39,671 quadratic tetrahedral (C3D10) elements and the inclusion 11,009 linear tetrahedral elements. The total simulation time is 2τ , divided into 100 fixed time steps:

ΔtFEM = τ /50



(27)

The calculation is performed in three steps. Step 1 consists in translating the sphere of a distance corresponding to −v × τ relative to the contact center. The second step consists in applying the maximum load assuming that the material behaves elastically only. The last step consists in combining the translation of the rigid spherical indenter of ΔxFEM = vΔtFEM and the resolution of the viscoelastic contact. Validation using Carbone and Putignano (2013) approach: Carbone and Putignano (2013) introduced a novel methodology to predict sliding and rolling friction of viscoelastic materials. Their model is based on the boundary element method and is limited to homogeneous viscoelastic materials. They obtained the steady-state solution, e.g. the transient solution observed at the beginning of the movement has not been modeled.

In order to facilitate the comparison of our model with the numerical simulations presented in Carbone and Putignano (2013), the same material properties and geometry are used. A rigid spherical indenter of radius R ¼10 mm moving at constant velocity on a viscoelastic half-space is considered. The viscoelastic properties of the semi-infinite body are given by Eqs. (28) and (29) and plotted in Fig. 8. Creep function:

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

11

Fig. 7. Finite element model used for the validation.

Fig. 8. Relaxation function (a) and creep function (b).

⎡ ⎛ t ⎞⎤ ⎡ 1 ⎛ t ⎞⎞⎤ 1 1 1⎛ ⎜⎜1 − exp ⎜− ⎟ ⎟⎟ ⎥ J (t) = ⎢ exp ⎜⎜− ⎟⎟ ⎥ = ⎢ − + ⎝ τ ⎠ ⎠ ⎥⎦ μ1 μ1 ⎝ ⎢⎣ μ 0 ⎝ τ ⎠ ⎥⎦ ⎣⎢ μ∞

(28)

Relaxation function:

⎡ ⎛ ⎞⎤ ⎢ ⎜ ⎟⎥ t R (t) = ⎢μ 0 + (μ∞ − μ 0 ) exp ⎜− μ ⎟ ⎥ 0 ⎟⎥ ⎢ ⎜ τ⎟ ⎜ ⎢⎣ ⎝ μ∞ ⎠ ⎥⎦

(29)

where μ is the spring stiffness with μ∞ = 3.86 MPa , μ∞ /μ0 = 10, and τ the relaxation time. The viscoelastic half-space is assumed to have a time-independent Poisson's ratio νM ¼0.3. The initial shear modulus is R (0) = μ∞ and the relaxation time τ = 0.01 s.

12

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

Fig. 9. Comparison of the normalized pressure distribution between the semi-analytical method (SAM), the finite element model (FEM), and the boundary element method (BEM) (see Carbone and Putignano, 2013 for BEM) in steady-state regime and for a fixed value of penetration δ = 0.1 mm .

3.1. Contact pressure distribution In this subsection comparison focuses on the pressure distribution obtained by both methods. Comparison is presented first (Fig. 9) for a fixed value of the prescribed normal displacement δ ¼0.1 mm, and second (Fig. 10) in the case of a fixed value of normal force FN = 0.15 N. Results are presented for the dimensionless moving velocity vτ /a⁎ = 0.8. an and P0 are the contact radius and peak pressure of the equivalent homogeneous elastic problem when the instantaneous shear modulus R (0) = μ∞ is used. An excellent agreement is observed between the three methods.

3.2. Apparent friction coefficient Carbone and Putignano (2013) presented the apparent friction coefficient as an output of their model. In this section the apparent coefficient of friction obtained with SAM and BEM are compared. It is however necessary to remind that the model presented in Carbone and Putignano (2013) gives the value of the apparent friction coefficient in the steady-state regime only, which corresponds to the asymptotic solution. Meanwhile the semi-analytical method allows to obtain this value in the transient regime as well, before to reach the steady-state solution. Fig. 11 presents the apparent friction coefficient for a constant dimensionless velocity vτ /a⁎ = 0.6. A good agreement is observed between the two methods where the stationary regime is reached thus validating the SAM framework. It is interesting to note that the steady-state solution is reached when t > τ /2, approximately.

Fig. 10. Comparison of the normalized pressure distribution between the semi-analytical method (SAM), the finite element model (FEM), and the boundary element method (BEM) (see Carbone and Putignano, 2013 for BEM) in steady-state regime and for a constant applied load FN = 0.15 N .

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

13

Fig. 11. Comparison of the apparent friction coefficient between the semi-analytical method (SAM) and the boundary element method (BEM) (see Carbone and Putignano, 2013) for a constant applied load FN = 0.15 N and dimensionless speed vτ /a⁎ = 0.6 .

4. Results 4.1. Description of the problem In this section the rolling contact between a rigid sphere of radius R¼10 mm moving at constant velocity on a heterogeneous viscoelastic half-space is considered. The creep function of the viscoelastic half-space is given by Eq. (28), with μ∞ = 3.86 MPa , μ∞ /μ0 = 10, and τ = 0.01 s. Results are presented first in the case of a viscoelastic homogeneous half-space. Afterwards a heterogeneous viscoelastic half-space containing one ellipsoidal inhomogeneity with different orientations and elastic properties is considered. The inhomogeneity is assumed here to have an isotropic elastic behavior. The simulation in the time domain varies between [0, 2τ], divided into 60 time steps. a⁎ (t) represents the instantaneous contact radius corresponding to the homogeneous viscoelastic case. The instantaneous Young's modulus of the viscoelastic half-space at the initial time is given by

E M (0) = 2 × (1 + ν M) × R (0) The ratio

γ=

(30)

γ can be expressed as

E I /E M (0)

(31)

I

where E is Young's modulus of the elastic inhomogeneity.

ν I = ν M = 0.3

(32)

ν I, νM are respectively Poisson's ratio of the inhomogeneity and the viscoelastic half-space. 4.2. Case of a homogeneous viscoelastic half-space (transient regime response) In this section two different cases are presented. The first case corresponds to a fixed value of the normal displacement

δ ¼0.1 mm, and the second one to a fixed value of the external normal force FN = 1.48 N. For the elastic case (t¼0) the first and second cases are equivalent in terms of contact pressure distribution and subsurface stresses. The main purpose of this section is to study the effect of the external load on the response in the transient regime.

 Prescribed normal displacement: The external load corresponds to a fixed value of the penetration δ = 0.1 mm and the simulation is done for different values of the dimensionless speed vτ /a0 = 0.4, 0.8, and 1.2. The initial solution at t ¼0 is not plotted since it corresponds to the classical Hertz solution. Figs. 12, 13a and b present the contact pressure distribution, apparent friction coefficient and resulting torque at different times, respectively. Again it is observed that the steady-state regime is reached at time t ≅ τ /2. It is also observed that the friction force acting against the movement of the sliding indenter (Fig. 13a), or the rolling torque in case of a rolling sphere (Fig. 13b), is strongly dependent on the moving velocity. The apparent friction coefficient reaches 0.045 for the dimensionless speed vτ /a0 = 1.2, the normal displacement imposed and the viscoelastic properties chosen. Fig. 14 presents the deformation of the viscoelastic halfspace induced by the rolling of the rigid sphere. The viscoelastic effect at the trailing edge of the contact is clearly visible as soon the movement starts, producing an asymmetric point which reaches an asymptotic shape in the steady-state

14

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

Fig. 12. Distribution of normalized contact pressure for homogeneous viscoelastic half-space at different time steps for different values of the dimensionless speed vτ /a⁎ for prescribed normal displacement δ = 0.1 mm . (a) vτ /a⁎ = 0.4 , (b) vτ /a⁎ = 0.8, and (c) vτ /a⁎ = 1.2.

regime.

 Prescribed normal force: The external load is a fixed normal force FN = 1.48 N. Again the Hertz solution which corresponds to the initial time t ¼0 is not plotted. Figs. 15, 16a and b present the contact pressure distribution, resulting torque and apparent friction coefficient at different times, respectively.

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

15

Fig. 13. Apparent friction coefficient (a) and resulting torque (b) versus dimensionless time for prescribed normal displacement δ = 0.1 mm .

The contact pressure profile is very different compared to the situation where the normal displacement is prescribed, for the same values of dimensionless speed. It is also interesting to observe that, this time, it takes more time to reach the steady-state regime, obtained here when t ≅ 4τ . This observation can be explained by going back to the material properties plotted in Fig. 8. One can remark that the relaxation function stabilizes at τ /2 while the creep function stabilizes near 4τ . For a prescribed normal displacement, the material response corresponds to a relaxation phenomenon while for a prescribed normal force this response corresponds to a creep phenomenon. It can be concluded that the time required to reach the steady-state regime for a rolling or sliding contact depends on the viscoelastic material properties, mostly on the relaxation function for prescribed normal displacement, and on the creep function for prescribed normal load. 4.3. Heterogeneous viscoelastic half-space The rigid sphere is now rolling on a viscoelastic solid containing one elastic inhomogeneity of ellipsoidal shape, with a translation velocity vτ /a⁎ = 0.6. Considering that dX I = (dx1I , dx2I , dx 3I ) is the center of the inhomogeneity in the R (O, X1, X2 , X3 ) coordinate system, it yields dX I = (0.6a⁎, 0.0, dx 3I ) when t¼0; dX I = (0.0, 0.0, dx 3I ) at t = τ ; and dX I = ( − 0.6a⁎, 0.0, dx 3I ) for t = 2τ . In this section the effects of the elastic properties of a spherical inhomogeneity (inclusion) and of the orientation of one ellipsoidal inhomogeneity are investigated. In what follows Young's modulus of the elastic inhomogeneity over the one of the viscoelastic matrix is defined by the dimensionless parameter γ = E I /E M (0), whereas θ gives the orientation of the ellipsoid. For all simulations Poisson's ratio of both inhomogeneity and matrix will be identical: ν I = ν M = 0.3. The contact pressure distribution will be plotted in the both coordinate systems R0 (O0, X10 , X20 , X30 ) and R (O, X1, X2 , X3 ), the first one being the fixed reference frame and the second one the reference frame linked to the contact (moving with). The configuration of the problem is recalled in Figs. 4 and 5. 4.3.1. Spherical inhomogeneity The case of a spherical inhomogeneity located at dX I = ( − 0.6a⁎, 0.0, 0.3a⁎) of radius r = 0.15a⁎ is considered. Young's modulus of the inhomogeneity is chosen in such a manner that it is identical to the viscoelastic matrix one at the beginning of the loading, i.e. γ = E I /E M (0) = 1. The results in terms of contact pressure distribution and apparent friction coefficient are presented in Figs. 17 and 18, respectively, for a prescribed normal penetration δ = 0.1 mm and dimensionless speed vτ /a0 = 0.6. In Fig. 17 for t ¼0 the presence of the inhomogeneity has no effect on the contact pressure distribution, and the result

Fig. 14. Deformation of the viscoelastic half-space at different time steps with a prescribed normal displacement δ = 0.1 mm , for the dimensionless speed vτ /a⁎ = 1.2.

16

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

Fig. 15. Distribution of the normalized contact pressure for a homogeneous viscoelastic half-space at different time steps and for different values of the dimensionless speed vτ /a⁎ in the case of a prescribed normal force FN = 1.48 N . (a) vτ /a⁎ = 0.4 , (b) vτ /a⁎ = 0.8, and (c) vτ /a⁎ = 1.2.

corresponds to the Hertzian solution, which is coherent with the choice of the elastic properties γ = E I /E M (0) = 1). Then, when t increases until one time the relaxation time τ approximately, one may observe a peak of pressure above the inhomogeneity. Finally the effect of the inhomogeneity on the contact pressure decreases with time and becomes negligible for t = 2τ . Two explications can be given: (i) the hysteresis losses induced by the rolling contact are stronger than the elastic

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

17

Fig. 16. Apparent friction coefficient (a) and resulting torque (c) versus dimensionless time in the case of a prescribed normal force FN = 1.48 N .

Fig. 17. Distribution of normalized contact pressure for a viscoelastic half-space containing a spherical inhomogeneity γ = E I /E M (0) = 1 (a1 = a2 = a3 = r = 0.25a⁎, dx 3 = 0.3a⁎); the value of the penetration is δ = 0.1 mm and the dimensionless speed is vτ /a0 = 0.6 ; (a) in the (moving) contact reference frame R (O, X1, X2, X3 ) , (b) in the fixed reference frame R0 (O0, X10 , X20 , X30 ) .

Fig. 18. Apparent friction coefficient versus dimensionless time for a viscoelastic half-space containing a spherical inhomogeneity γ = E I /E M (0) = 1 (a1 = a2 = a3 = r = 0.25a⁎, dx 3 = 0.3a⁎); with prescribed normal penetration δ = 0.1 mm and dimensionless speed vτ /a0 = 0.6 .

18

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

Fig. 19. Distribution of the dimensionless contact pressure for different orientations θ when the viscoelastic half-space contains an ellipsoidal inhomogeneity γ = E I /E M (0) = 1 (a1 = 0.4a⁎ , a2 = 0.5a⁎, a3 = 0.2a⁎, dx 3 = 0.4a⁎); with prescribed normal penetration δ = 0.1 mm and dimensionless speed vτ /a0 = 0.6 ; θ = 30° ((a) in R (O, X1, X2, X3 ) , (b) in R0 (O0, X10 , X20 , X30 ) ); θ = 45° ((c) in R (O, X1, X2, X3 ) , (d) in R0 (O0, X10 , X20 , X30 ) ); θ = 60° ((e) in R (O, X1, X2, X3 ) , (f) in R0 (O0, X10 , X20 , X30 ) ); θ = 80° ((g) in R (O, X1, X2, X3 ) , (h) in R0 (O0, X10 , X20 , X30 ) ).

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

19

Fig. 20. Apparent friction coefficient for different orientations θ when the viscoelastic half-space contains a tilted ellipsoidal inhomogeneity γ = E I /E M (0) = 1 (a1 = 0.4a⁎, a2 = 0.5a⁎, a3 = 0.2a⁎ , dx 3 = 0.4a⁎); with prescribed normal penetration δ = 0.1 mm and dimensionless speed vτ /a0 = 0.6 .

contribution of the inhomogeneity and (ii) for t = 2τ the center of the inhomogeneity is located at dX I = ( − 0.6a⁎, 0.0, 0.3a⁎) which is out of the actual contact area [ − 0.5a⁎, a⁎] × [ − a⁎, a⁎]. It is obvious that a spherical inhomogeneity has an effect on the contact pressure distribution (see Leroux et al., 2010; Leroux and Nelias, 2011). However from Fig. 18 a spherical inhomogeneity has no effect on the apparent friction coefficient (or the resisting torque). This is due to the symmetry of the inhomogeneity regarding the normal to the contact surface. According to Eshelby's formalism the (additional) contribution of a single spherical inhomogeneity on eigenstrains, then on the surface displacement and the contact pressure distribution, is symmetric. We will show now that only titled ellipsoidal inhomogeneity or non-symmetric inhomogeneities may modify the apparent friction coefficient, which is a quite interesting result for practical applications. 4.3.2. Ellipsoidal inhomogeneity In this section the effect of a single tilted ellipsoidal inhomogeneity, of orientation angle θ, is considered (see Figs. 4 and 5). The effects of its orientation and Young's modulus are studied. Results are plotted in the plane x2 = 0.

 Effect of the orientation: Results are presented for different values of θ. The influence of the orientation angle is observed



for the case γ = 1. Results in terms of contact pressure distribution and apparent friction coefficient are presented in Figs. 19 and 20 for different values of the orientation angle θ. When θ tends to 90°, the inhomogeneity gets very close to the surface and the peak of pressure due to the presence of the inhomogeneity increases. The apparent friction coefficient presents an increase of up to 7% compared to the homogeneous viscoelastic solution in the transient regime, which appears when the observation time corresponds to t = τ /2 i.e. when the location of the heterogeneity coincides with the sphere center. This increase vanishes later when the inhomogeneity moves away from the contact. Effect of Young's modulus: Figs. 21 and 22 present the results for different elastic properties of the inhomogeneity. The relatively soft inhomogeneity with γ = 0.1 corresponds to E I = 2 × (1 + ν) × μ0 . The magnitude of the contact pressure increases when the inhomogeneity becomes stiffer (γ = 1, γ = 4 ). Conversely when the inhomogeneity is softer than the matrix, the contact pressure drops locally compared to the solution with a homogeneous viscoelastic half-space. The apparent friction coefficient increases respectively up to 4% and 6% for γ = 1 and γ = 4 in the transient regime. It should also be noted that in the steady-state regime (t > τ /2), a reduction of the apparent friction coefficient is observed for relatively soft tilted inhomogeneities (γ = 0.1) whereas stiffer inhomogeneities (γ = 1 or γ = 4 ) tend to increase it.

4.3.3. Multiple interacting inhomogeneities It has been shown that the presence of a tilted ellipsoidal inhomogeneity has an effect on the apparent friction coefficient when it goes through the contact, however limited in time. The idea here is to investigate what happens when the viscoelastic body contains several tilted ellipsoidal inhomogeneities, that are continuously moving through the contact. It should be mentioned that the interaction between close inhomogeneities is considered, as already described by Leroux et al. (2010). As in Zhou et al. (2011b), a modified conjugate gradient algorithm is used to account for the interaction between inhomogeneities. The same approach was also employed by Zhou (2012) to predict the effective moduli of periodic composites with arbitrary inhomogeneity distribution. The bulk and shear moduli obtained are well located within the Hashin– Shtrikman bounds which demonstrate the accuracy of the method. The case of three ellipsoidal inhomogeneities of semiaxes (a1 = 0.4a⁎, a2 = 0.5a⁎, a3 = 0.2a⁎) with γ = 0.5 and θ = 60° is considered. The center of these inhomogeneities is located

20

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

Fig. 21. Distribution of the dimensionless contact pressure for different values of γ when the viscoelastic half-space contains a single tilted ellipsoidal inhomogeneity (a1 = 0.4a⁎, a2 = 0.5a⁎, a3 = 0.2a⁎, dx 3 = 0.4a⁎, θ = 30°); with prescribed normal penetration δ = 0.1 mm and dimensionless speed vτ /a0 = 0.6 ; γ = 0.1 ((a) in R (O, X1, X2, X3 ) ), (b) in R0 (O0, X10 , X20 , X30 ) ); γ = 1 ((c) in R (O, X1, X2, X3 ) , (d) in R0 (O0, X10 , X20 , X30 ) ); γ = 4 ((e) in R (O, X1, X2, X3 ) , (f) in R0 (O0, X10 , X20 , X30 ) ).

at depth dx 3 = 0.4a⁎, and at abscissa dx1= −0.1a⁎, 0.4a⁎, and 0.9a⁎, respectively. Results in terms of contact pressure distribution are presented in Fig. 23a for the homogeneous case and in Fig. 23b for the heterogeneous one. One can remark that the presence of the three inhomogeneities modifies very significantly the contact pressure distribution. The effect of inhomogeneities on the apparent friction coefficient, see Fig. 24, is this time more pronounced than when just one is present (see Figs. 20 and 22 for comparison). It can be explained by the fact that the volume fraction occupied by the elastic

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

21

Fig. 22. Apparent friction coefficient for different values of γ when the viscoelastic half-space contains an ellipsoidal inhomogeneity (a1 = 0.4a⁎, a2 = 0.5a⁎, a3 = 0.2a⁎, dx 3 = 0.4a⁎, θ = 30°); with prescribed normal penetration δ = 0.1 mm and dimensionless speed vτ /a0 = 0.6 .

Fig. 23. Distribution of normalized contact pressure for a viscoelastic half-space containing three interacting tilted ellipsoidal inhomogeneities γ = E I /E M (0) = 0.5 (a1 = 0.4a⁎, a2 = 0.5a⁎, a3 = 0.2a⁎ , dx 3 = 0.4a⁎, θ = 60° ); the value of the penetration is δ = 0.1 mm and the dimensionless speed is vτ /a0 = 0.6 : (a) homogeneous, (b) heterogeneous.

Fig. 24. Apparent friction coefficient when the viscoelastic half-space contains three interacting tilted ellipsoidal inhomogeneities γ = E I /E M (0) = 0.5 (a1 = 0.4a⁎, a2 = 0.5a⁎, a3 = 0.2a⁎, dx 3 = 0.4a⁎, θ = 60°); with prescribed normal penetration δ = 0.1 mm and dimensionless speed vτ /a0 = 0.6 .

22

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

inhomogeneities – in the vicinity of the contact – is this time higher. 4.3.4. Subsurface stresses In addition, the model proposed here allows to obtain the subsurface stresses. As an illustration the von Mises stress is plotted for a spherical inhomogeneity of radius r = 0.15a⁎ located at depth dx 3I = 0.3a⁎ and with γ = 0.4 . Note that the inhomogeneity is initially softer (relatively) than the viscoelastic matrix (E I = 0.4E ∞) and then becomes progressively harder than the matrix at very low strain rate (E I = 4E 0 ). Results of the numerical simulation are presented in Fig. 25 at different time steps when the inhomogeneity is moving through the contact. The von Mises stress is plotted in the fixed coordinate system R0 (O0, X10 , X20 , X30 ). At the initial time t¼0 the maximum value of the von Mises stress is maximum at the equator of the spherical inhomogeneity/matrix interface, reaching 0.8P0 , significantly higher than the homogeneous solution (0.6P0 ). It remains maximum in the vicinity of the inhomogeneity as far as this heterogeneity remains within the contact, however the absolute value drops quickly to 0.14P0 for t = τ /2 while moving to the bottom side of the inhomogeneity/matrix interface. When the inhomogeneity moves away from the contact (t = 2τ ) the maximum of the von Mises stress corresponds to 0.1a⁎ located at x1/a⁎ = 1 in the viscoelastic matrix. This observation is consistent since the peak of pressure for a homogeneous viscoelastic body is located at the same position. Interacting inhomogeneities may also be considered. Fig. 26 presents the dimensionless subsurface stress σ33/P0 at t ¼0 and t = τ /3, with three tilted ellipsoidal inhomogeneities as for Figs. 23b and 24. The subsurface stresses are disturbed by the presence of the inhomogeneities. It is also possible to account for a cluster of inhomogeneities (Leroux et al., 2010) or to describe a real composite material (Leroux et al., 2013).

Fig. 25. Von Mises stress field at different time steps when the viscoelastic half-space contains a spherical inhomogeneity of radius r = 0.15a⁎ located at depth dx 3I = 0.3a⁎ (γ = 0.4 ); prescribed normal penetration δ = 0.1 mm and dimensionless speed vτ /a0 = 0.6 ; (a) t = 0 , (b) t = τ /2 , (c) t = τ , (d) t = 3⁎τ /2, (e) t = 2⁎τ .

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

23

Fig. 26. σ33/P0 in the plane x2 = 0 at different time steps when the viscoelastic half-space contains three interacting tilted ellipsoidal inhomogeneities γ = E I /E M (0) = 0.5 (a1 = 0.4a⁎, a2 = 0.5a⁎, a3 = 0.2a⁎ , dx 3 = 0.4a⁎); prescribed normal penetration δ = 0.1 mm and dimensionless speed vτ /a0 = 0.6 ; (a) t = 0 , (b) t = τ /3.

5. Conclusion A numerical method has been proposed to model the rolling contact problem between a rigid sphere – or a sliding indenter – and a heterogeneous viscoelastic half-space. The model has been first validated by performing a comparison with the results of a finite element model and with the model recently presented by Carbone and Putignano (2013) based on the boundary element method. Contact pressure distribution, subsurface stresses, apparent friction coefficient and the resulting torque induced by the viscoelastic hysteresis losses have been analyzed in the presence of a single spherical or ellipsoidal heterogeneity. The main results are summarized below: 1. The contact pressure distribution in the steady-state regime depends on the way the load is imposed, either by a prescribed normal penetration or a fixed normal load. 2. The time laps required to reach the steady-state regime differs also if the load is imposed through a prescribed normal penetration or via a constant normal load. It depends mostly on the relaxation function for prescribed normal displacement, and on the creep function for prescribed normal load. 3. Whereas the presence of heterogeneous inhomogeneities nearby the contact surface modifies the contact pressure distribution, it has been observed in transient regime that spherical inhomogeneities have no effect on the apparent friction coefficient and subsequently on the rolling torque. Only titled ellipsoidal or asymmetric inhomogeneities modify the apparent friction coefficient in the transient regime. 4. As a consequence in the steady-state regime a reduction of the apparent friction coefficient is observed for relatively soft inhomogeneities whereas stiffer inhomogeneities tend to increase it, which is a very interesting outcome for practical applications.

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

F11 (x1, x2 , x 3 ) =

⎡ x ⎞ x 2 − x22 x 3 x22 3x 3 x12 ⎤ 1 ⎢ 1 − 2ν ⎛ ⎥, ⎜⎜1 − 3 ⎟⎟ 1 + − 2 2 3 ρ⎠ 2π ⎢⎣ r r ρ ρ5 ⎥⎦ ⎝

F22 (x1, x2 , x 3 ) = F11 (x2 , x1, x 3 ),

F33 (x1, x2 , x 3 ) = −

F12 (x1, x2 , x 3 ) =

3 x 33 , 2π ρ5

x ⎞x x x x x 3x 3 x2 x1 ⎤ 1 ⎡⎢ 1 − 2ν ⎛ ⎥, ⎜1 − 3 ⎟ 1 2 + 3 2 1 − 2π ⎢⎣ r 2 ⎝ ρ ⎠ r2 ρ3 ρ5 ⎥⎦

F13 (x1, x2 , x 3 ) = −

3 x1x 32 , 2π ρ5

F23 (x1, x2 , x 3 ) = F12 (x2 , x1, x 3 ),

24

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

where

r 2 = x12 + x22 , with

ρ=

x12 + x 32 + x 32 ,

ν being Poisson's ratio of the isotropic half-space.

Appendix B. Stresses in a half-space subject to normal pressure (Mij ) An isotropic half-space is submitted a uniform normal pressure sn in a discretized surface area of 2Δx1 × 2Δx2 at the center point P (x1′ , x2′ , 0) . The stress at an observation point Q (x1, x2, x 3 ) is given by (Zhou et al., 2009a; Johnson, 1985)

σij (x1, x2 , x 3 ) = Mij (x1 − x1′ , x2 − x2′ , x 3 ) σ n (x1, x2 ) σij (x1, x2 , x 3 ) =

σn [hij (ξ1 + Δx1, ξ2 + Δx2 , ξ3 ) − hij (ξ1 + Δx1, ξ2 − Δx2 , ξ3 ) 2π

+ hij (ξ1 − Δx1, ξ2 − Δx2 , ξ3 ) − hij (ξ1 − Δx1, ξ2 + Δx2 , ξ3 )]

(B.1)

where

ξ i = xi − xi′. The functions hij () in Eq. (B.1) are defined by

h11 (x1, x2 , x 3 ) = 2ν tan−1

x22 + x 32 − ρx2 ρ − x2 + x 3 x1x2 x 3 + 2(1 − ν) tan−1 + , x1x 3 x1 ρ (x12 + x 32 )

h22 (x1, x2 , x 3 ) = h11 (x2 , x1, x 3 ),

h33 (x1, x2 , x 3 ) = tan−1

h12 (x1, x2 , x 3 ) = −

h13 (x1, x2 , x 3 ) = −

⎞ x22 + x 32 − ρx2 x1x2 x 3 ⎛ 1 1 ⎜⎜ ⎟⎟, − + 2 2 2 2 ρ ⎝ x1 + x 3 x1x 3 x2 + x 3 ⎠

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

x2 x 32 ρ (x13 + x 32 )

,

h23 (x1, x2 , x 3 ) = h13 (x2 , x1, x 3 ), where

ρ=

x12 + x22 + x 32 .

Appendix C. Normal displacement at the surface subjected to normal pressure (K n) The contact between a sphere and an elastic half-space having respectively elastic constants (E1, ν1) and (E2, ν2 ) , where the surface x 3 = 0 is discretized into rectangular surface area of 2Δ1 × 2Δ2, is now considered. The initial contact point coincides with the origin of the Cartesian coordinate system (x1, x2, x 3 ). The relationship between the normal displacement at an observation point P (ξ1, ξ2, 0) and the pressure field at the center Q (ξ1′, ξ2′ , 0) is built using the function Kn:

⎤ ⎡ 1 − ν12 1 − ν22 ⎥ K n (c1, c2 ) = ⎢ + ⎢⎣ πE1 πE2 ⎥⎦

4

∑ K pn (c1, c2 ), p=1

⎛ (c2 + Δ2 ) + K1n (c1, c2 ) = (c1 + Δ1 ) log ⎜ ⎜ ⎝ (c2 − Δ2 ) +

⎞ (c2 + Δ2 )2 + (c1 + Δ1 )2 ⎟ , ⎟ (c2 − Δ2 )2 + (c1 + Δ1 )2 ⎠

K.E. Koumi et al. / J. Mech. Phys. Solids 80 (2015) 1–25

⎛ (c1 + Δ1 ) + K2n (c1, c2 ) = (c2 + Δ2 ) log ⎜ ⎜ ⎝ (c1 − Δ1 ) +

⎞ (c2 + Δ2 )2 + (c1 + Δ1 )2 ⎟ , ⎟ (c2 + Δ2 )2 + (c1 − Δ1 )2 ⎠

⎛ (c2 − Δ2 ) + K3n (c1, c2 ) = (c1 − Δ1 ) log ⎜ ⎜ ⎝ (c2 + Δ2 ) +

⎞ (c2 − Δ2 )2 + (c1 − Δ1 )2 ⎟ , ⎟ (c2 + Δ2 )2 + (c1 − Δ1 )2 ⎠

⎛ (c1 − Δ1 ) + K4n (c1, c2 ) = (c2 − Δ2 ) log ⎜ ⎜ ⎝ (c1 + Δ1 ) +

⎞ (c2 − Δ2 )2 + (c1 − Δ1 )2 ⎟ , ⎟ (c2 − Δ2 )2 + (c1 + Δ1 )2 ⎠

25

where

c1 = ξ1 − ξ1′

and

c2 = ξ2 − ξ2′.

References Aleksandrov, V., Goryacheva, I., Torskaya, E., 2010. Sliding contact of a smooth indenter and a viscoelastic half-space (3d problem). Dokl. Phys. 55, 77–80. Argatov, I., 2012. An analytical solution of the rebound indentation problem for an isotropic linear viscoelastic layer loaded with a spherical punch. Acta Mech. 223, 1441–1453. Carbone, G., Putignano, C., 2013. A novel methodology to predict sliding/rolling friction in viscoelastic materials: theory and experiments. J. Mech. Phys. Solids 61, 1822–1834. Chaise, T., Nelias, D., 2011. Contact pressure and residual strain in 3d elasto-plastic rolling contact for a circular or elliptical point contact. J. Tribol. 133 (4), 041402. Chen, W., Liu, S., Wang, Q., 2008a. Fft-based numerical methods for elasto-plastic contacts of nominally flat surfaces. J. Appl. Mech. 75. 011022-1–11. Chen, W., Wang, Q., Huan, Z., Luo, X., 2011. Semi-analytical viscoelastic contact modeling of polymer-based materials. J. Tribol. 133, 041404. Chen, W., Wang, Q., Wang, F., Keer, L., Cao, J., 2008b. Three-dimensional repeated elasto-plastic point contact, rolling and sliding. J. Appl. Mech. 75. 021021-1–12. Chiu, Y., 1978. On the stress field and surface deformation in a half space with a cuboidal zone in which initial strains are uniform. J. Appl. Mech. 45, 302–306. Goriacheva, G., 1973. Contact problem of rolling of a viscoelastic cylinder on a base of the same material. J. Appl. Math. Mech. 37 (5), 925–933. Graham, G., 1967. The contact problem in the linear theory of viscoelasticity when the time dependent contact area has any number of maxima and minima. Int. J. Eng. Sci. 5, 495–514. Greenwood, J., 2010. Contact between an axisymmetric indenter and a viscoelastic half-space. Int. J. Mech. Sci. 52, 829–835. Hunter, S., 1960. The Hertz problem for a rigid spherical indenter and viscoelastic half-space. J. Mech. Phys. Solids 8, 219–234. Hunter, S., 1961. The rolling contact of a rigid cylinder with a viscoelastic half space. J. Appl. Mech. 28, 611–617. Jacq, C., Nelias, D., Lormand, G., Girodin, D., 2002. Development of a three-dimensional semi-analytical elastic–plastic contact code. J. Tribol. 124 (4), 653–667. Johnson, K., 1985. Contact Mechanics. Cambridge University Press, Cambridge. Kalker, J., 1972. On elastic line contact. ASME J. Appl. Mech. 39, 1125–1132. Kalker, J., 1977. The surface displacement of an elastic half space loaded in a slender, bounded, curved surface region with application to the calculation of the contact pressure under a roller. J. Inst. Math. Appl. 19, 127–144. Koumi, K., Nelias, D., Chaise, T., Duval, A., 2014a. Modeling of the contact between a rigid indenter and a heterogeneous viscoelastic material. Mech. Mater. 77, 28–42. Koumi, K., Zhao, L., Leroux, J., Chaise, T., Nelias, D., 2014b. Contact analysis in the presence of an ellipsoidal inhomogeneity within a half space. Int. J. Solids Struct. 51, 1390–1402. Le Tallec, P., Rahler, C., 1994. Numerical models of steady rolling for non-linear viscoelastic structures in finite deformations. Int. J. Numer. Methods Eng. 37, 1159–1186. Lee, E., Radok, J., 1960. The contact problem for viscoelastic bodies. J. Appl. Mech. 27, 438–444. Leroux, J., Fulleringer, B., Nelias, D., 2010. Contact analysis in presence of spherical inhomogeneities within a half-space. Int. J. Solids Struct. 47, 3034–3049. Leroux, J., Nelias, D., 2011. Stick-slip analysis of a circular point contact between a rigid sphere and a flat unidirectional composite with cylindrical fibers. Int. J. Solids Struct. 48, 3510–3520. Leroux, J., Nelias, D., Ruiz-Sabariego, J.-A., 2013. Modeling a frictional contact for composite materials [modelisation d’un contact frottant pour materiaux composites]. Mater. Techn. 101 (2). 205-1/7 (in French). Liu, S., Wang, Q., Liu, G., 2000. A versatile method of discrete convolution and fft (dc-fft) for contact analyses. Wear 243 (1–2), 101–111. Love, A., 1952. A Treatise on the Mathematical Theory of Elasticity, 4th edition. Cambridge University Press, Cambridge. Mura, T., 1987. Micromechanics of Defects in Solids, 2nd edition. Kluwer Academic Publishers, Dordrecht. Nackenhorst, U., 2004. The ale-formulation of bodies in rolling contact: theoretical foundations and finite element approach. Comput. Method Appl. Mech. Eng. 193, 4299–4322. Nasdala, L., Kaliske, M., Becker, A., Rothert, H., 1998. An efficient viscoelastic formulation for steady-state rolling structures. Comput. Mech. 22, 395–403. Nelias, D., Antaluca, E., Boucly, V., 2007. Rolling of an elastic ellipsoid upon an elastic–plastic flat. J. Tribol. 129 (4), 791–800. Nelias, D., Boucly, V., Brunet, M., 2006. Elastic–plastic contact between rough surfaces: proposal for a wear or running-in model. J. Tribol. 128 (2), 236–244. Padovan, J., Kazempour, A., Tabaddor, F., Brockman, B., 1992. Alternative formulations of rolling contact problems. Finite Elem. Anal. Des. 11, 275–284. Padovan, J., Paramadilok, O., 1984. Transient and steady state viscoelastic rolling contact. Comput. Struct. 20, 545–553. Panek, C., Kalker, J., 1977. A solution for the narrow rectangular punch. J. Elast. 2, 213–215. Panek, C., Kalker, J., 1980. Three-dimensional contact of a rigid roller traversing a viscoelastic half space. J. Inst. Math. Appl. 26, 299–313. Persson, B., 2010. Rolling friction for hard cylinder and sphere on viscoelastic solid. Eur. Phys. J. E 33, 327–333. Zhou, K., 2012. Elastic field and effective moduli of periodic composites with arbitrary inhomogeneity distribution. Acta Mech. 223, 293–308. Zhou, K., Chen, W., Keer, L., Wang, Q., 2009a. A fast method for solving three-dimensional arbitrarily shaped inclusions in a half space. Comput. Methods Appl. Mech. Eng. 198 (9–12), 885–892. Zhou, K., Keer, L., Wang, Q., 2011b. Semi-analytic solution for multiple interacting three-dimensional inhomogeneous inclusions of arbitrary shape in an infinite space. Int. J. Numer. Methods Eng. 87, 617–638.