Daniel Nélias

Efoe Rodrigue

Thibaut Chaise

Daniel Nelias
email: daniel.nelias@insa-lyon.fr

Rolling Contact on a Viscoelastic Multi-layered Half-Space

Keywords: Viscoelastic material, Multi-layered half-space, Semi-Analytical Method, Transient analysis

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Various motions can occur when two bodies get into contact. One of the most common among them is the rolling contact. When the bodies in contact are considered completely elastic, the pure rolling and the pure sliding contacts are the same if there is no Coulomb's friction at the contact interface. Otherwise, when friction is considered at the contact interface, under tangential forces, the rolling contact shows sliding zones beside adhesion zones. This kind of rolling, called tractive rolling, is the general case of rolling contact. It comes clear that when there is no tangential (or tractive) force, the tractive rolling case coincides with the pure rolling contact. And, when the tractive force is higher than the product of the normal force by the friction coefficient, the contact becomes a pure sliding contact (Coulomb's law).

The general tractive rolling contact has been tackled first by [START_REF] Carter | On the action of a locomotive driving wheel[END_REF], followed by a certain number of scientists (see for example [START_REF] Kalker | On the rolling contact of two elastic bodies in the presence of dry friction[END_REF], [START_REF] Bentall | Slip in the rolling contact of two dissimilar elastic rollers[END_REF], [START_REF] Johnson | Contact mechanics[END_REF], [START_REF] Nowell | Tractive rolling of dissimilar elastic cylinders[END_REF]) among whom [START_REF] Kalker | On the rolling contact of two elastic bodies in the presence of dry friction[END_REF][START_REF] Kalker | Three-Dimensional Elastic Bodies in Rolling Contact[END_REF][START_REF] Kalker | Viscoelastic multilayered cylinders rolling with dry friction[END_REF] has given the most complete theory. The Kalker's theory covers the rolling contact in its various aspects.

In more recent works, the tractive rolling is modelled taking advantage of improvements in numerical techniques. For instance, [START_REF] Guler | Tractive rolling contact mechanics of graded coatings[END_REF] developed a semi-analytical method for solving the tractive rolling contact but for graded coatings. [START_REF] Wang | A numerical approach for analyzing three-dimensional steady-state rolling contact including creep using a fast semi-analytical method[END_REF] developed a numerical approach for analysing three-dimensional steady-state rolling contact including creep using a semianalytical method. [START_REF] Manyo | Tire-pavement tractive rolling contact under turning conditions: towards pavement top-down cracking[END_REF] further designed a similar semi-analytical method based on the Kalker's theory coupled with Conjugate Gradient algorithms and Fast Fourier Transforms.

Besides, the modelling of multi-layered structures is certainly of great interest for those who work with coated bodies or composites for example. The latter structures are widely popular and used in many fields. For an elastic bilayer structure i.e. a single elastic layer on an elastic substrate, a first resolution has been proposed by O'Sullivan and [START_REF] O'sullivan | Sliding contact stress field due to a spherical indenter on a layered elastic half-space[END_REF] and followed by others [START_REF] Nogi | Influence of a hard surface layer on the limit of elastic contact-part I: Analysis using a real surface model[END_REF], [START_REF] Plumet | A 3-D model for a multilayered body loaded normally and tangentially against a rigid body: Application to specific coatings[END_REF], [START_REF] Polonsky | A fast and accurate method for numerical analysis of elastic layered contacts[END_REF], [START_REF] Liu | Studying contact stress fields caused by surface tractions with a discrete convolution and fast fourier transform algorithm[END_REF], [START_REF] Stepanov | Modeling of Sliding of a Smooth Indenter over a Viscoelastic Layer Coupled with a Rigid Base[END_REF], ...). For two elastic layers on an elastic substrate, [START_REF] Chen | Computation of stresses and displacements in a layered elastic medium[END_REF] performed an exact analysis for a two dimensional case. For a general case of L layers, L being any positive integer, most of the literature works are based on finite element methods [START_REF] Chen | Finite element analysis of contact induced adhesion failure in multilayer coatings with weak interfaces[END_REF], [START_REF] Djabella | Finite element analysis of elastic stresses in multilayered systems[END_REF], [START_REF] Gorishnyy | Optimization of wear-resistant coating architectures using finite element analysis[END_REF], [START_REF] Kot | Contact mechanics of coating-substrate systems: Monolayer and multilayer coatings[END_REF], [START_REF] Komvopoulos | Finite element analysis of a layered elastic solid in normal contact with a rigid surface[END_REF], [START_REF] Assogba | Numerical investigation of the mechanical response of semi-rigid base asphalt pavement under traffic load and nonlinear temperature gradient effect[END_REF]). However, [START_REF] Yu | Analytical frequency response functions for contact of multilayered materials[END_REF] solved the problem analytically in the Fourier frequency domain by constructing a matrix system and used numerical algorithms to find a semi analytical solution for the influence coefficients. In geology science, the study of the ground lead to the modelling of the multi-layered systems (see for instance [START_REF] Kuo | Static response of a multilayered medium under inclined surface loads[END_REF], [START_REF] Bufler | Theory of elasticity of a multilayered medium[END_REF], [START_REF] Farrell | Deformation of the Earth by surface loads[END_REF], [START_REF] Peltier | The impulse response of a Maxwell Earth[END_REF]) but in 2-D.

In addition to the multi-layered structure's influence, when some of the layers in a multi-layered structure have a non negligible viscoelastic behaviour, it changes the whole behaviour of the structure. This can easily be seen by observing roads in summer, when the top layer's viscoelastic properties are emphasized by higher temperature. Then, the whole road does not respond mechanically as in winter. The modelling of viscoelastic contacts has been treated by many authors from [START_REF] Lee | The contact problem for viscoelastic bodies[END_REF] to recently [START_REF] Xu | Influence of viscoelastic mechanical properties of rail pads on wheel and corrugated rail rolling contact at high speeds[END_REF], [START_REF] Zhang | Transient and steady-state viscoelastic contact responses of layer-substrate systems with interfacial imperfections[END_REF], [START_REF] Nguyen | Time-stepping method for frictional contact of anisotropic viscoelastic solids[END_REF], [START_REF] Wallace | Three-dimensional rolling/sliding contact on a viscoelastic layered halfspace[END_REF], [START_REF] Menga | Exploring the effect of geometric coupling on friction and energy dissipation in rough contacts of elastic and viscoelastic coatings[END_REF], through [START_REF] Kalker | Viscoelastic multilayered cylinders rolling with dry friction[END_REF], [START_REF] Goryacheva | Contact of elastic bodies with thin visco-elastic coatings under conditions of rolling or sliding friction[END_REF], Wayne [START_REF] Wayne Chen | Semi-analytical viscoelastic contact modeling of polymerbased materials[END_REF] [START_REF] Koumi | Modeling of the contact between a rigid indenter and a heterogeneous viscoelastic material[END_REF][START_REF] Koumi | Rolling contact of a rigid sphere/sliding of a spherical indenter upon a viscoelastic half-space containing an ellipsoidal inhomogeneity[END_REF], [START_REF] Putignano | A review of boundary elements methodologies for elastic and viscoelastic rough contact mechanics[END_REF]; [START_REF] Putignano | Mechanics of rough contacts in elastic and viscoelastic thin layers[END_REF][START_REF] Putignano | Theory of reciprocating contact for viscoelastic solids[END_REF], [START_REF] Menga | The sliding contact of a rigid wavy surface with a viscoelastic half-space[END_REF][START_REF] Menga | Effect of thickness and boundary conditions on the behavior of viscoelastic layers in sliding contact with wavy profiles[END_REF], Stepanov andTorskaya (2016, 2018), [START_REF] Kusche | Frictional force between a rotationally symmetric indenter and a viscoelastic half-space: Elastomer friction[END_REF], [START_REF] Scaraggi | Rough contact mechanics for viscoelastic graded materials: The role of small-scale wavelengths on rubber friction[END_REF], [START_REF] Putignano | Viscoelastic reciprocating contacts in presence of finite rough interfaces: A numerical investigation[END_REF], [START_REF] Goryacheva | Modelling of the viscoelastic layer effect in rolling contact[END_REF], [START_REF] Torskaya | Effect of surface layers in sliding contact of viscoelastic solids 3-D model of material[END_REF], and others. The modelling of viscoelastic contacts is now widely known and solved for various motions (rolling, sliding, reciprocating, fretting), various punches (in term of geometry and mechanical properties) and with different methods (Boundary Element Methods, Semi-Analytical Methods, Finite Element Methods).

For the first time a method is proposed to solve the 3D contact problem with a viscoelastic multi-layered half space. Some examples are also provided to illustrate the performances of the numerical techniques used to tackle the problem. The exposed method allows to solve the tractive rolling contact for any geometry of punch and in both transients and steady state regimes, however in this paper we only focus on the effect of viscosity for a multi-layered half-space and a frictionless contact with a rolling or sliding sphere. A further paper will complete this work with details on the tractive rolling modelling in viscoelasticity. In the following, firstly we present our problem, find the elastic influence coefficients for the multi-layered half-space in a second time, perform an Elastic/Viscoelastic correspondence to account for the viscoelastic behaviour of each layer, and then we make a parametric study in a pure rolling case in transient and steady-state regimes. The system being considered consists of two bodies (1 and 2) in contact over a region Γ c .

Formulation

Body 1 can be of any kind of geometry, rigid, elastic or even viscoelastic, and body 2 is the actual viscoelastic multi-layered half-space. In this framework, we consider that the layers are perfectly bonded to one another. It means that the displacements are continuous at the interfaces.

In addition, we add the traction continuity condition at the interfaces.

The classical definition of the contact for numerical applications begins with the definition of the gap between the bodies: the surface separation g(x, y, t) between the two bodies is defined at every point at the contact surface as the sum of the initial separation g i (x, y), the elastic and the viscoelastic normal displacements of both bodies under the load history u (B 1 +B 2 ) 3 (x, y, t) and, if there is some, the rigid body displacement δ(t) at every time t. Then, we have the following equation:

g(x, y, t) = g i (x, y) + δ(t) + u (B 1 +B 2 ) 3 (x, y, t) (1)
Once the gap is defined, let us call p(x, y, t) the local pressure at a point (x, y) on the surface at a time t. Then we can identify the contact conditions:

(x, y) ∈ Γ c ; g(x, y, t) = 0 and p(x, y, t) > 0 (2a) (x, y) Γ c ; g(x, y, t) > 0 and p(x, y, t) = 0 (2b)

The above contact conditions take into account the non-penetration and the pressure positivity conditions. We add the load balance: if a normal force F N (t) is applied on the system, it must be equal to the integral of the contact pressure distribution p(x, y, t) over the contact zone Γ c . This reads as follow:

F N (t) = Γ c (t) p(x, y, t)dS (3)
It is important to remark that the contact is solved at every time step. It means that whatever the load history is, it is taken into account at the present step. This allow us to solve transient contact problems, as shown later in the results section.

In the above contact definition, the elastic and viscoelastic displacements remain unknown. In the next part, we will find the elastic response of the multi-layered half-space under a unit pressure and a unit shear. This response is the influence coefficient that we will formulate in terms of Papkovich-Neuber potentials.

Influence coefficients for multi-layered elastic half-space

Consider a multilayered half-space with L linear elastic and homogeneous layers on a linear elastic and homogeneous substrate. From the top to the bottom, the layer j shear modulus is µ j , its Poisson's ratio is ν j and its thickness is h j ( j = 1, ..., L). For the substrate, the shear modulus is µ L+1 or µ sub and the Poisson's ration is ν L+1 or ν sub . We use the Papkovich-Neuber potentials φ j and ψ j i to express the elastic field. φ j and ψ j i are harmonic functions of (x, y), and when no body forces are present, [START_REF] Malvern | Introduction to the mechanics of a continuous medium[END_REF] showed that ψ j 2 can vanish. Then we can write displacements and stresses as functions of the potentials φ j , ψ j

x and ψ j z :

2µ j u j i = φ j ,i + xψ j x,i + z j ψ j z,i -(3 -4ν j )ψ j i (4a) σ j ik = φ j ,ik -2ν j (ψ j x,x + ψ j z,z )δ ik -(1 -2ν j )(ψ j i,k + ψ j k,i ) + xψ j x,ik + z j ψ j z,ik (4b)
where the Einstein's notation is used and j = 1, ...L + 1.

In the above equations, the coma stand for the partial derivative (eg: φ ,x = ∂φ ∂x ). To avoid performing those partial derivatives, a linearization of Eqs. (4a) and (4b) is done using Fourier Transform (FT). Recalling the properties of FT:

FT x ∂ ∂x φ(x) = im φ(m) (5a) FT x [-ixφ(x)] = ∂ ∂m φ(m). ( 5b
)
Furthermore, we define the transformed Papkovich-Neuber potentials using bi-harmonic func-

tions: φj = A j exp -αz j + A j exp αz j (6a) ψ j x = B j exp -αz j + B j exp αz j (6b) ψj z = C j exp -αz j + C j exp αz j . (6c)
Therefore, we have linear equations in the Fourier frequency domain to solve.

Let's introduce the following notations: w j + = e αz j and w j -= e -αz j , where α = √ m 2 + n 2 is the radius in the frequency domain; (m, n) represents the Fourier Transform of (x, y).

We can, then, infer the forms of displacements and stresses in the Fourier frequency domain (see Appendix A).

In those displacements and stresses, the following A j , A j ,B j , B j ,C j and C j terms are unknown.

Note that A L+1 = 0, B L+1 = 0 and C L+1 = 0 because the elastic field must vanish at infinity. Thus, we need a system of 6L + 3 equations for finding the unknowns.

To construct the system to solve, we have to set the boundary conditions first. They are:

• At the top surface, which is the contact surface of the whole multi-layered half-space i.e. at z 1 = 0, pressure and shear fields are imposed:

σ 1 zz (x, y, 0) = -δ p (x, y), σ 1 zx (x, y, 0) = -δ q (x, y), σ 1 zy (x, y, 0) = 0. ( 7
)
where δ p and δ q are Dirac distributions, with subscripts p and q for normal and tangential loading, respectively.

• At all the interfaces, tractions and displacements continuity is required: ( j, j+1) 4(1 -ν j+1 )µ ( j, j+1) 2α(1

σ j zx (x, y, h j ) = σ j+1 zx (x, y, 0) u j x (x, y, h j ) = u j+1 x (x, y, 0)      2α(1 -ν 1 ) -2α(1 -ν 1 ) . . . . . . 0 -4(1 -ν j )θ j - -4(1 -ν j )θ j + 4(1 -ν j+1 )µ
-ν j )θ j --2α(1 -ν j )θ j + -2α(1 -ν j+1 ) 2α(1 -ν j+1 ) . . . . . . . . . . . . 0 -4(1 -ν L )θ L - -4(1 -ν L )θ L + 4(1 -ν L+1 )µ (L,L+1) 2α(1 -ν L )θ L --2α(1 -ν L )θ L       α α 2(1 -ν 1 ) -2(1 -ν 1 ) α -α 1 -2ν 1 1 -2ν 1 0 . . . . . . X 1 j X 2 j . . . . . . α α(θ L + ) 2 αh L αh L (θ L + ) 2 -αµ (L,L+1) θ L + 0 α -α(θ L + ) 2 3 -4ν L + αh L (3 -4ν L + αh L )(θ L + ) 2 -αµ (L,L+1) θ L + -µ (L,L+1) (3 -4ν L+1 )θ L + 0 α -α(θ L + ) 2 1 -2ν L + αh L (1 -2ν L -αh L )(θ L + ) 2 -αθ L + -(1 -2ν L+1 )θ L + α α(θ L + ) 2 2(1 -ν L ) + αh L -[2(1 -ν L ) -αh L ](θ L + )      A 1 A 1 . . . C j C j A j+1 A j+1 . . . C L C L A L+1 C L+1 ( 10
)
where θ j + = e αh j , θ j -= e -αh j , µ ( j, j+1) = µ j µ j+1 and the sub-matrices X 1 j and X 2 j are written as follow:

X 1 j =                          α α(θ j + ) 2 αh j αh j (θ j + ) 2 α -α(θ j + ) 2 3 -4v j + αh j 3 -4v j -αh j (θ j + ) 2 α -α(θ j + ) 2 1 -2v j + αh j 1 -2v j -αh j (θ j + ) 2 α α(θ j + )       -αµ ( j, j+1) θ j + -αµ ( j, j+1) θ j + 0 0 -αµ ( j, j+1) θ j + αµ ( j, j+1) θ j + -µ ( j, j+1) 3 -4v j+1 θ j + -µ ( j, j+1) 3 -4v j+1 θ j + -αθ j + αθ j + -1 -2v j+1 θ j + -1 -2v j+1 θ j + -αθ j + -αθ j + -2 1 -v j+1 θ j + 2 1 -v j+1 θ j +                          (11)
In Eqs. ((10)), the right hand terms remain to be defined. Those S 0 1 , • • • , S L 4 are functions of the B j and B j found solving Eqs. ((9)) and their derivatives. All those functions are given in Appendix B.

The systems of equations ((9)) and ((10)) can be solved to find explicit expressions of the variables A j , A j , B j , B j , C j and C j . Those explicit solutions have been given by [START_REF] Yu | Analytical frequency response functions for contact of multilayered materials[END_REF]. In the present work, we have chosen to find the exact analytical solution only for the Eq. ( (9)). For the system of equations ((10)), a numerical inversion is performed. This will help us avoid errors that can be made since the solutions are long and complex. Thus, we use a Gauss-Jordan algorithm to make the inversion of the (4L + 2) × (4L + 2) matrix. The cost of the numerical inversion depends, of course, on the number of layers, and also on the size of the computational area.

At this point, the Green functions are found in the Fourier frequency domain. For finding them in the space domain, we have to perform an Inverse Fourier Transform (IFT) (see [START_REF] Liu | Studying contact stress fields caused by surface tractions with a discrete convolution and fast fourier transform algorithm[END_REF]). In the present work, the IFT is handled with Fast Fourier Transform (FFT) algorithms as proposed by [START_REF] Liu | A versatile method of discrete convolution and FFT (DC-FFT) for contact analyses[END_REF].

Once the solution is found, we have the solution of the elastic contact on a multi-layered half-space with any L layers.

In the next section, we find a way to turn this elastic field to a viscoelastic one, and, to transform the static solution in order to make a quasi-static rolling motion.

Viscoelastic rolling contact

𝝁 𝟎 𝝁 𝟏 𝝁 𝒏 𝜼 𝟏 𝜼 𝒏 Figure 2: Generalized Maxwell model
A viscoelastic material is a one which mechanical behaviour's changes in time. In this work we will consider a generalized Maxwell model (see Fig. 2) and its Prony series decomposition:

R(t) =        µ 0 + n i=0 µ i exp(-t/τ i )        H(t) (12)
where, µ i is the spring stiffness, η i the dashpot viscosity, The Relaxation function in Eq. (( 12)) is related to the creep function with the relation:

τ i = η i /µ i the
t 0 J(ξ)R(t -ξ)dξ = t. ( 13
)
For one relaxation time, the creep and relaxation functions can directly be written as: Further, we will use this model (see Fig. 3) to compute simple configurations in section 3.

J(t) = 1 µ 0 - 1 µ 1 exp - t τ = 1 µ ∞ + 1 µ 1 1 -exp - t τ (14) R(t) =        µ 0 + (µ ∞ -µ 0 ) exp        - t µ 0 µ ∞ τ               (

Elastic/Viscoelastic correspondence

The Elastic/Viscoelastic correspondence that we use is the same as explained in [START_REF] Wallace | Three-dimensional rolling/sliding contact on a viscoelastic layered halfspace[END_REF]. We first discretize the computation time interval into regular time steps. Then, at a time step t we (i) replace the elastic modulus 1 2µ j by the corresponding creep function J j (t) and (ii) integrate the history.

Let's recall first the correspondence in the case of homogeneous bodies. For a viscoelastic body, the displacement u i at a time t can be written as:

u i (x, y, t) = t 0 R 2 F i j (x -x ′ , y -y ′ , t -t ′ ) σ j (x ′ , y ′ , t ′ )dx ′ dy ′ dt ′ = t 0 R 2 F i j (x ′ , y ′ , t -t ′ ) σ j (x -x ′ , y -y ′ , t ′ )dx ′ dy ′ dt ′ (16)
where σ j are the components of the applied stress at the surface and F i j (x, y, t) the viscoelastic influence coefficients. The dot˙represents the partial derivative with respect to time.

Consider a homogeneous semi-infinite body, the viscoelastic influence coefficients F i j (x, y, t) can be written as:

F i j (x, y, t) = J(t)G i j (x, y), ( 17
)
where the term G i j is included in the elastic influence coefficient and does not depend on the elastic properties of the body. J(t) is the creep function defining the viscoelastic behaviour of the body.

The reader may find this correspondence method for a homogeneous half-space in Wayne Chen et al. ( 2011) for instance. Thus the above displacement can be written as:

u i (x, y, t) = t 0 J(t -t ′ ) R 2 G i j (x ′ , y ′ ) σ j (x -x ′ , y -y ′ , t ′ )dx ′ dy ′ dt ′ . ( 18
)
This means that a convolution is performed between the elastic influence coefficients G i j -which are constant -and the applied stress at every time step. And then the time integral, containing the creep function and the history of the loading, is calculated. Now, when moving to the case of a multi-layered half-space, the factorisation in Eq. ( 17) is not possible in that same way.

What we have done for this model is:

✧ Firstly, we have observed that the use of the potentials of Papkovich-Neuber implies the following forms of the displacements:

2µ r u r i = φ r ,i + xψ r x,i + zψ r z,i -(3 -4ν)ψ r i , ( 19
)
where µ is the shear modulus of the layer r and φ and ψ i are the potentials of Papkovich-Neuber.

In the left hand term, the factor 2µ is ideal for the correspondence since it has to be replaced by its viscoelastic equivalent J(t).

✧ Secondly, solving the boundary conditions, it has come that the right hand term depends on the ratios µ k /µ k+1 , k being the subscript corresponding to the layer k according to the convention in the paper. Since we have supposed that the Poisson's ratios do not change with time, the change of behaviour of the material is given only by those ratios µ k /µ k+1 .

✧ The last remark have lead us to propose the following equation for the viscoelastic displacements:

u i (x, y, t) = t 0 J(t -t ′ ) R 2 Ḡ1 i j (x ′ , y ′ , J k (t ′ )/J k+1 (t ′ )) σ j (x -x ′ , y -y ′ , t ′ )dx ′ dy ′ dt ′ . ( 20
)
This equation is an approximation of the exact solution. It meets the exact solution in steady state regime and in some other cases as discussed in Appendix C. In transient regime, we believe the error made is marginal with regard to the behaviour of the Fourier transform (see Appendix C for details).

Therefore, we can perform the correspondence but we need to recalculate the right hand term of Eq. ( 19) at every time step. This recalculation of the influence coefficients has an additional cost on the overall simulation time. However, the simulation remains straight and fast.

The resulting equation for the normal displacement at the contact surface that we are seeking for in Eq. ((1)) is:

u B 2 z (x, y, t) = t 0 J 1 (t -ξ) Ḡ1 (x, y, ξ) * ∂p(x, y, ξ) ∂ξ dξ (21)
where,

• * represents the convolution product,

• J 1 (t) is the creep function of the top layer,

• Ḡ1 is contained -as explained above -in the Green function (the influence coefficients) for the normal displacement at the surface z = 0 of the half-space, in which the ratios µ k /µ k+1 are replaced by their equivalent J k (t)/J k+1 (t).

This correspondence can be applied for the two bodies in contact and then a contact between two viscoelastic multi-layered bodies can be simulated. For simplicity, we will restrain the results in this paper to the case an elastic body 1 in contact with a viscoelastic multi-layered body 2.

Rolling/sliding approach

The rolling contact can be a pure rolling motion (i.e with no traction force, so there is no need to take into account the friction coefficient) or a tractive rolling contact (where a traction force is present). Note that a frictionless sliding motion is equivalent to pure rolling [START_REF] Koumi | Rolling contact of a rigid sphere/sliding of a spherical indenter upon a viscoelastic half-space containing an ellipsoidal inhomogeneity[END_REF]).

In this paper, we consider a free rolling motion for a body of revolution over the half-space. A further paper will present the tractive rolling contact.

Let's assume that the time range of simulation is divided into N t uniform time steps of size ∆t, the elementary displacement of the body 1 is then ∆ r = v × ∆t. From an initial framework where the fixed coordinate system is represented by the superscript 0, to the actual framework, the normal displacement is:

u B 2 z (x, y, t) = t 0 J 1 (t -ξ)       Ḡ1 (x, y, ξ) * ∂p(x 0 -α∆r 0 x , y 0 -α∆r 0 y , ξ) ∂ξ       dξ. ( 22
)
It shows that we make a translation of the actual pressure field and the relaxed pressure fields from the previous step in the fixed framework with a vector α∆ r 0 α∆r 0 x , α∆r 0 y , α being the actual time step index.

Numerical approach

The numerical implementation needs a discretization of both space and time. The contact surface is discretized into N 1 × N 2 elements, while the time domain is subdivided into N t time steps ∆t. Then, using the finite difference method for ∂p ∂ξ the discretized form of Eq. (( 22)) becomes:

u B 2 z (i, j, α) = α k=0 J 1 [(α -k)∆t] ×         N 1 i ′ =1 N 2 j ′ =1 Ḡ1 i -i ′ , j -j ′ , k∆t × p(i ′ , j ′ , k) -p(i ′ , j ′ , k -1)         , (23)
where one can infer that i, j and α are the indices for x, y and t respectively. For k = 0, the term with k -1 has to be zero.

Note that the translation process explained in the previous section does not clearly appear in the notation of the pressure here; so, the reader should make sure to perform the translation between the actual and the precedent time steps.

To find the solution of the contact problem, we apply the variational formulation and use the Conjugate Gradient Method (CGM) algorithm to solve it (see [START_REF] Gallego | A fast and efficient contact algorithm for fretting problems applied to fretting modes I, II and III[END_REF] for the com-plete method). The convolution product terms are evaluated using a Discrete Convolution-Fast Fourier Transform (DC-FFT) algorithm to reduce the number of operations from (N 1 × N 2 ) 2 to [START_REF] Liu | A versatile method of discrete convolution and FFT (DC-FFT) for contact analyses[END_REF]).

(N 1 × N 2 × log(N 1 × N 2 )) (see the algorithm of
At the time step α, the only unknown in the above equation is the actual pressure field. The implementation gives:

1. For α = 0, the elastic solution of the contact is found: p(i, j, 0).

2. For any α, the Green functions are recalculated Ḡ1 (ii ′ , jj ′ , (α)∆t). Then, the term

H(α) = α-1 k=0 J 1 [(α -k)∆t] ×         N 1 i ′ =1 N 2 j ′ =1 Ḡ1 i -i ′ , j -j ′ , k∆t × p(i ′ , j ′ , k) -p(i ′ , j ′ , k -1)         -J 1 [(α -α)∆t] ×         N 1 i ′ =1 N 2 j ′ =1 Ḡ1 i -i ′ , j -j ′ , (α)∆t × p(i ′ , j ′ , α -1)         (24)
is calculated and stored.

The remaining term

J 1 [(α -α)∆t] × N 1 i ′ =1 N 2 j ′ =1 Ḡ1 (i -i ′ , j -j ′ , (α)∆t) × (p(i ′ , j ′ , α)) con-
tains the unknown pressure field and so is solved with the CGM algorithm.

When the surface pressure is known at the time step α, the 3D stresses can be computed using the actual contact pressure and the elastic stresses influence coefficients.

For the calculation of the three-dimensional stresses, two possibilities can be considered:

➢ Applying the correspondence to the stresses.

In this case, the above correspondence for the displacements has to be performed also for the stresses. This will induce a big amount of time for the calculation of the stresses. We have not chosen this option because this process leads to a precision which we have judged not necessary compared to the demand of memory and time.

➢ Instead, we have chosen to use the contact field which is solved in viscoelasticity with no restriction, as input for the calculation of the stresses. The elastic influence coefficients are taken for the purpose. It is a modelling choice that has been made by evaluating the cost versus the accuracy brought.

Using the developed model, a parametric study is performed in order to highlight the effects of some of the viscoelastic parameters on the rolling contact. The results are presented in the next section.

Results

Framework

For the presentation of some illustrative cases, let's first describe the chosen framework.

We apply a normal force F n or impose a normal rigid displacement δ on a sphere of radius R = 10mm which goes in contact with a multi-layered viscoelastic half-space. The sphere is free-rolling at velocity v along the direction x. Since we want to focus on the effects related to the viscoelastic multi-layered half-space, we set for the sphere a shear modulus µ sphere = 10 20 MPa and its Poisson's ratio is ν sphere = 0.3; this very high shear modulus will free us of its elastic effect.

For the layered half-space, we use the notations in Fig. 1 and set L = 10 layers, each of thickness

h j = a 0 /10 ( j = 1, • • • , 10).
The viscoelastic creep and relaxation functions are the ones in Eqs.

(( 14)) and (( 15)). In those equations we will use the subscript j for the layer j (from top to bottom) and the subscript L + 1 = 11 or sub for the substrate. For all, layers and substrate, we keep µ ∞ j /µ 0 j = 10 ( j = 1, • • • , L + 1), and we do not change the substrate mechanical characteristics which are: µ ∞sub = µ ∞11 = 3.86MPa and τ sub = τ 11 = 10s. It comes obviously that the substrate which remains the same will serve us as the reference. Thus, the contact considering no layer gives the useful Hertz parameters p 0 (the maximum contact pressure) and a 0 (the radius of the contact zone) in the elastic case.

We will perform the studies by changing the instantaneous modulus of the layers µ ∞ j or/and its relaxation time τ j ( j = 1, • • • , 10). The instantaneous shear moduli will change following three main cases (see Fig. 4):

• The "Alternate Modulus": the even index layers take the shear modulus value 4 × µ ∞sub while the odd index layers take the shear modulus value µ ∞sub : • The "Decreasing Modulus": the modulus decreases from 6 × µ ∞sub with a constant step to reach the substrate modulus µ ∞sub :

µ ∞ j = µ ∞ j+1 4 = µ ∞sub ( j = 1, 3, 5, 7, 9). ( 25
µ ∞ j = 6 - j -1 2 × µ ∞sub ( j = 1, • • • , 10). (26)
• The "Increasing Modulus": the modulus increases from µ ∞sub to 5.5×µ ∞sub in the 10th layer:

µ ∞ j = 1 + j -1 2 × µ ∞sub ( j = 1, • • • , 10). ( 27
)
We can replicate the cases above to have similar changes of the relaxation times. Then we have,

τ j = τ j+1
2 = τ sub (with j = 1, 3, 5, 7, 9) for the "Alternate Relaxation" time;

τ j = 3 -j-1 5 × τ sub (with j = 1, • • • , 10
) for the "Decreasing Relaxation" time and τ j = 1 + j-1 5 × τ sub (with j = 1, • • • , 10) for the "Increasing Relaxation" time.

Applied Normal Force

Let's begin our analysis by applying a concentrated force F n = 1.48N on the sphere. The rolling velocity is given by vτ sub /a 0 = 0.8. Then, we compute in the time domain [0, 6τ sub ] subdivided into 600 time steps. The contact surface is discretized with elements of size 0.2a 0 × 0.2a 0 .

- Figure 5 shows the evolution of the contact pressure along the rolling direction axis x. Results are plotted at (a) t = 0, (b) t = τ sub /4, (c) t = τ sub , (d) t = 5 × τ sub for the homogeneous halfspace and the three non-homogeneous half-spaces where only the instantaneous shear moduli of the layers vary (see section 3.1).

The first remark about Fig. 5 is that at t = 0, the elastic response correlates well with the results of [START_REF] Yu | Analytical frequency response functions for contact of multilayered materials[END_REF] confirming that our method for solving the matrix system in Eq. ((10)) works well. For the viscoelasticity, the method described in section 2.2.1 has already been validated by [START_REF] Wallace | Three-dimensional rolling/sliding contact on a viscoelastic layered halfspace[END_REF], thus we will not show again a validation of it. In another hand, one can observe the significant decrease of the contact pressure peak (at the edge of the contact zone) at t = τ sub /4 . At the time t = τ sub , the pressure presents the shape of a steady-state regime but will continue decreasing to reach the steady-state regime (t = 5τ sub ). One may also notice that the more severe situation is found for a decreasing shear modulus, for which the maximum pressure is always the highest. Note also that, since both surfaces in contact experience the same contact pressure, it may impact the resistance of the counterface.

The calculations for an applied force is more time consuming as shown previously by [START_REF] Koumi | Rolling contact of a rigid sphere/sliding of a spherical indenter upon a viscoelastic half-space containing an ellipsoidal inhomogeneity[END_REF]. This, because we need more time to reach the steady-state regime and also because we have to perform the computation over a wider surface. Thus, we will perform most of the parametric study with a prescribed normal rigid body displacement. Keep in mind that it will activate a relaxation phenomenon.

Applied Normal Rigid Body Displacement

In this section, we prescribe a normal displacement δ = 0.1a 0 and compute in the time domain [0, 2τ sub ] with 80 time steps. We discretize the space with elements of size 0.2a 0 × 0.2a 0 × 0.05a 0 .

We still use the same rigid sphere, and the rolling velocity is still given by vτ sub /a 0 = 0.8. We, then, analyse the effects of the variation of the instantaneous shear moduli, the effects of the variation of the relaxation times, and the effects of both.

Variation of the instantaneous shear moduli

Despite the absence of Coulomb's friction coefficient, the viscoelastic dissipation introduces an assymetry of the contact pressure distribution along the rolling direction, which in fine produces a resisting force or a resisting torque. This phenomenon is also encoutered for frictionless elasticplastic rolling/sliding contact, as shown in [START_REF] Chaise | Contact Pressure and Residual Strain in 3D Elasto-Plastic Rolling Contact for a Circular or Elliptical Point Contact[END_REF]. For more details regarding frictionless rolling/sliding viscoelastic contact, the reader may refer to Fig. 3 What is called here the apparent friction coefficient µ app , sometimes also called the viscous friction elsewhere, is a representation of the hysteresis losses in the viscoelastic material. Again it is observed for a frictionless contact when inelastic materials are in contact. These losses induce a tangential resisting force given by:

F T (t) = Γ c (t) ∂u 3 (x 1 , x 2 , t) ∂x 1 p (x 1 , x 2 , t) dx 1 dx 2 (28)
where Γ c (t) is the real contact area at time t.

Then the apparent friction coefficient is µ app = F T /F N , where F N is the normal load.

Firstly, we focus on the effect of the change of the instantaneous shear moduli of the layers.

For this purpose, we keep the relaxation times the same for all the layers

τ j = τ sub ( j = 1 • • • , 10)
and modify the instantaneous shear moduli following the cases described in section 3.1. Doing so, we can plot various outputs such as the contact pressure at different times, the apparent friction coefficient which expresses the resistance to the rolling, and the three-dimensional stresses.

Figure 6 shows the contact pressure along the rolling axis x at times (a) t = 0, (b) t = τ sub /10, (c) t = τ sub /4 and (d) t = τ sub . We compare to a homogeneous viscoelastic half-space with the mechanical properties of the substrate, three non-homogeneous half-spaces with 10 layers having the same relaxation time but with either increasing modulus, decreasing modulus or alternate modulus (see Fig. 4).

In all the non-homogeneous cases, the overall stiffness of the layers is higher than the one of the substrate. This is the reason why the contact pressure of the homogeneous shows the lowest maximum with a wider profile. Moreover, among those cases, the top layer is stiffer for the decreasing modulus. As shown for a bilayer by O' Sullivan and King (1988), [START_REF] Wang | A numerical approach for analyzing three-dimensional steady-state rolling contact including creep using a fast semi-analytical method[END_REF],

and [START_REF] Wallace | Three-dimensional rolling/sliding contact on a viscoelastic layered halfspace[END_REF], the top surface has the greater direct influence on the contact. Under this remark, it comes comprehensive that the decreasing modulus shows a higher maximum pressure with the smallest contact zone. For the alternate modulus case and the increasing modulus case, the top layer is the same. Thus, the results allow us to say that the overall stiffness of the alternate modulus case is lower than the of the increasing modulus Further, when the time increases, the trend will stay the same since all the layers have the same relaxation time. In addition, we see the through a transient regime to achieve seems be a steady-state in sub-figure (c) and (d) of Fig. 6. This transition from the transient to the steady-state regime is confirmed by the apparent friction coefficient in Fig. 7 which shows that the steady-state regime is achieved quite quickly; far before t = 0.5τ sub . Moreover, one can see the small effect of the instantaneous shear moduli on the resistance to the rolling. This correlates with results found for only one layer lying on the top of the half-space ( Wallace et al.

(

)). 2020
Beside the contact results, the sub-surface stresses are of undeniable interest. Figure 8 shows the dimensionless second invariant of the stress tensor √ J 2 /p 0 , at t = 0 (left) which corresponds to the elastic case, and at t = 2τ sub which is in the steady-state regime.

Alternate Modulus

t = 0 t = 2τ
Decreasing Modulus

Increasing Modulus Again, the 3D-stresses distribution correlates well with the results of [START_REF] Yu | Analytical frequency response functions for contact of multilayered materials[END_REF] at t = 0.

One can clearly see that the effect of alternate modulus a jump of the stress from a layer to another. The higher stresses appear when we move (from top to bottom,) from a stiffer layer to a less stiffer one; the exception being near the surface where the modulus is low (the top layer has the modulus of the substrate). For the increasing and the decreasing modulus cases, there is no jump of stresses between the layers and in the decreasing modulus case, the stresses distribution is closer to a homogeneous modulus case distribution. This remark, can mean that for a decreasing modulus case, the overall half-space behaves like a homogeneous half-space with a greater instantaneous shear modulus. Otherwise, the increasing modulus case shows an increase of the stress from top to bottom in the layers, and reaches a maximum at the interface with the substrate. At that interface, there is a jump from a very high value (the maximum) to a very low value in the substrate.

At the steady-state regime, for the three plotted cases, the stresses relax (their value decrease)

and the maxima change location. For each, the maximum goes then under the leading edge of the contact and very close to the surface. This should be correlated with the contact pressure distribution (see Fig. 6).

Variation of relaxation times

Let's now focus on the effect of the change of the relaxation times of the layers. For this purpose, we keep the instantaneous shear moduli the same for all the layers µ ∞ j = µ ∞sub ( j = 1 • • • , 10) and modify the relaxation times following the cases describes in section 3.1.

Figure 9 shows the contact pressure distribution along the x axis at different time steps. The elastic response of the three cases is the same as expected since the instantaneous moduli are the same for all the layers and the substrate. For the alternate, increasing and homogeneous relaxation cases, the difference is not such significant. The reason lies in the fact that the top layer is the same for those three cases; and as discussed earlier, the top layer influences the most the contact behavior. Otherwise, for the decreasing relaxation case, the result is different. The top layer in that case has a greater time. Thus, the contact pressure relaxes less rapidly. And, even when the is achieved, the pressure profile is less wide and the maximum value with 10 layers with different relaxation times; normal rigid displacement δ = 0.1a 0 prescribed and the dimensionless rolling velocity is given by vτ/a 0 = 0.8 is greater.

When we cast a look at the apparent friction coefficient, see Fig. 10, we see the similarity of the behaviour for all the cases but the decreasing relaxation time one. For the latter, not only the steady-state regime is achieved later, but also the apparent friction coefficient is nearly 20% greater.

Figure 11 presents the 3D stress distribution. At time t = 0 there is no difference as expected.

In the steady-state regime, here for t = 2τ, the distribution of stresses is completely different due to viscoelastic effects, however the difference between the various configurations of layers is not very pronounced. We recall that for the calculation of the stresses, we perform a convolution of the actual contact pressure with the elastic influence coefficients. This can explain that we do not observe a huge change between the alternate and increasing relaxation cases since the contact pressures also are close (see Fig. 9). The decreasing relaxation case shows again a narrower distribution of stresses with a higher maximum stress. However, in all the cases, there is a change of the distribution in the steady-state regime; the stresses relax and the maximum moves under the leading edge of the surface and closer to the surface.

Variation of both instantaneous shear moduli and relaxation times

We now consider a mixed framework: the increasing modulus case coupled with the decreasing relaxation time case. The goal is to show a more complex case, under an applied force or an applied body displacement. For both cases, we discretize the surface with elements of size 0.2a 0 × 0.2a 0 .

And, we perform the computations in the time domain [0, 6τ sub ] for the applied force case and in [0, 2τ sub ] for the applied body displacement.

Figure 12 shows the contact pressure profile along the x axis at various time steps for an applied force (left) and for an applied body displacement (right). One can immediately see the huge difference between the profiles in the transient and in the steady-state regime when there is an applied force compared to an applied body displacement. Furthermore, we can see once again that the steady-state regime is achieved earlier for an applied body displacement. Confirmation of this is given with the apparent coefficient plotted in 13. latter figure shows also for an applied force, the to rolling that for an applied body displacement.

Conclusion

In summary, for the time, a viscoelastic multi-layered half-space is modelled for a rolling contact problem with a semi-analytical method. Drawn from the method using the Papkovich-Neuber potentials and Fourier transformation, we constructed matrix systems for the multi-layered problem and then we applied an elastic/viscoelastic correspondence. Further, with the numerical tools of semi-analytical methods (mainly CGM and FFT algorithms), we developed a code for solving the problem efficiently. For example, a three-dimensional computation for L = 10 layers, for 80 time steps as in section 3.3, for the calculation of both the contact pressure distribution and the stress field within the viscoelastic multi-layered half-space, takes 8017s CPU (≈ 2h15min) on a personnal computer.

Using the developed model, we highlighted some effects of the variation of the instantaneous shear moduli and of the variation of the relaxation time, in a pure rolling case. To do so, we took the Prony representation of a generalized Maxwell viscoelastic model, with one relaxation time. Then, setting a frictionless rolling motion, we have analyzed the evolution with time of the contact pressure distribution and the 3D stress field. Alongside the known results for the decreasing, increasing, and alternate elastic modulus within the layers, the main new observations adding viscoelasticity are:

✔ The variation of the shear moduli of the layers does not have a significant effect on the resistance to motion of the rolling contact; but its effect is undeniable on both the contact pressure profile and the 3D stress distribution.

✔ When the relaxation time is changing from one layer to the other one, the top layer which is interfering directly in the contact plays a key role on the contact pressure distribution and on the apparent friction coefficient.

✔ In the steady-state regime, the maxima of stress in the half-space move below the surface at the leading edge of the rolling sphere whatever are the properties of the viscoelastic layers.

Finally, this study has been voluntarily limited to the frictionless or pure rolling contact problem.

As perspective, a coupling of this model with the tractive rolling algorithm is to be done to cover the complete rolling contact in viscoelasticity for multi-layered bodies.

Appendix A.

With the notations of section 2.1, one can express the Fourier Transform form of the elastic displacements and stresses with the Papkovich-Neuber potentials as follow:

2µ j ũj x =im D j w j -+ D j w j + -4 1 -ν j B j w j -+ B j w j + + m 2 z j α -1 B j w j --B j w j + + imz j C j w j -+ C j w j + (A.1a) 2µ j ũ j y =in D j w j -+ D j w j + -mnz j α -1 B j w j --B j w j + + inz j C j w j -+ C j w j + (A.1b) 2µ j ũj z = -α D j w j --D j w j + -imα -1 B j w j --B j w j + -3 -4ν j C j w j -+ C j w j + + imz j B j w j -+ B j w j + ) + z j α -C j w j -+ C j w j + (A.1c) σ j xx = -m 2 D j w j -+ D j w j + -2im 2 -ν j B j w j -+ B j w j + ) -2αν j -C j w j -+ C j w j + + im 3 z j α -1 B j w j --B j w j + -z j m 2 C j w j -+ C j w j + (A.1d) σ j yy = -n 2 D j w j -+ D j w j + -2imν j B j w j -+ B j w j + -2αν j -C j w j -+ C j w j + + imn 2 z j α -1 B j w j --B j w j + -z j n 2 C j w j -+ C j w j + (A.1e) σ j zz =α 2 D j w j -+ D j w j + + 2im 1 -ν j B j w j -+ B j w j + -2α 1 -ν j -C j w j -+ C j w j + -imz j α B j w j --B j w j + + z j α 2 C j w j -+ C j w j + (A.1f) σ j xy = -mn D j w j -+ D j w j + -2in 1 -ν j B j w j -+ B j w j + + im 2 nz j α -1 B j w j --B j w j + -z j mn C j w j -+ C j w j + (A.1g) σ j xz = -imα D j w j --D j w j + + 2α 1 -ν j + m 2 α -1 × B j w j --B j w j + -im 1 -2ν j C j w j -+ C j w j + -m 2 z j B j w j -+ B j w j + + iz j αm -C j w j -+ C j w j + (A.1h) σ j yz = -inα D j w j --D j w j + + mnα -1 B j w j --B j w j + -in 1 -2ν j C j w j -+ C j w j + -mnz j B j w j -+ B j w j + + inz j α -C j w j -+ C j w j + (A.1i) + ( j = 2, • • • , L) (B.5) B L+1 = 1 -v L 1 -v L+1 B L θ L --B L θ L + (B.6)
To write the second members of Eq. ((10)), we need the derivatives of the above functions.

These derivatives about m are:

t (L) 0,m = 0 (B.7a) t ( j-1) 0,m = 4µ ( j-1, j) (θ j -) 2 t ( j) 0,m -t ( j) 0 + 1 t ( j) 0 -1 mh j α t ( j) 0 + 1 -t ( j) 0 -1 (θ j -) 2 2 ( j = L, • • • , 2) (B.7b) B 1 ,m = (θ 1 -) 2 t (1) 0 + 1 -t (1) 0 -1 (θ 1 -) 2 -1 2α (1 -v 1 )                  2t (1) 0,m -( 2mh 1 α ) t (1) 0 + 1 t (1) 0 -1 t (1) 0 + 1 -t (1) 0 -1 (θ 1 -) 2         - m α 2          (B.8) B 1 ,m = B 1 ,m + m 2α 3 (1 -ν 1 ) (B.9) B j ,m = 1 -v j-1 1 -v j t ( j) 0 -1 (θ j -) 2 t ( j) 0 + 1 -t ( j) 0 -1 (θ j -) 2              -2mh j α t ( j) 0 + 1 + 2t ( j) 0m t ( j) 0 -1 t ( j) 0 + 1 -t ( j) 0 -1 (θ j -) 2 B j-1 θ j-1 --B j-1 θ j-1 + +B j-1 ,m θ j-1 --B j-1 ,m θ j-1 + - mh j-1 α B j-1 θ j-1 -+ B j-1 θ j-1 + ( j = 2, • • • , L) (B.10) B j ,m = 1 -v j-1 1 -v j t ( j) 0 + 1 t ( j) 0 + 1 -t ( j) 0 -1 (θ j -) 2              2t ( j) 0,m t ( j) 0 +1 (θ j -) 2 + t ( j) 0 -1 -2mh j α (θ j -) 2 t ( j) 0 + 1 -t ( j) 0 -1 (θ j -) 2 B j-1 θ j-1 --B j-1 θ j-1 + +B ( j-1) ,m θ j-1 --B j-1 ,m θ j-1 + - mh j-1 α B j-1 θ j-1 -+ B j-1 θ j-1 + ( j = 2, • • • , L) (B.11) B L+1 ,m = 1 -v L 1 -v L+1 B L ,m θ L -- mh L α θ L -B L -B L ,m θ L + - mh L α θ L + B L (B.12)
Here, we can express the functions of the second member of Eq. ((10)):

S (0) 1 = -α -1 -2imα -1 (1 -v 1 ) B 1 + B 1 -iα B 1 ,m + B 1 ,m (B.13) S (0) 2 = -i mα -1 B 1 -B 1 + α B 1 ,m -B 1 ,m = 0 (B.14) and for j = 1, • • • , L: S ( j) 1 = αi µ ( j, j+1) B j+1 ,m θ j + + B j+1 ,m θ j + + 1 α mh j B j -B j (θ j + ) 2 -B j ,m -B j ,m (θ j + ) 2 (B.15) S ( j) 2 =i µ ( j, j+1) m α B j+1 θ j + -B j+1 θ j + + µ ( j, j+1) α B j+1 ,m θ j + -B j+1 ,m θ j + +mh j B j + B j (θ j + ) 2 - m α B j -B j (θ j + ) 2 -α B j ,m -B j ,m (θ j + ) 2 (B.16) S ( j) 3 = i α -1 m B j+1 θ j + -B j+1 θ j + + α B j+1 ,m θ j + -B j+1 ,m θ j + -α -1 -h j mB j + α -1 + h j m(θ j + ) 2 B j -α B j ,m -B j ,m (θ j + ) 2 (B.17) S ( j) 4 =i 2mα -1 1 -v j+1 B j+1 θ j + + B j+1 θ j + + α B j+1 ,m θ j + + B j+1 ,m θ j + -2mα -1 1 -v j B j + B j (θ j + ) 2 + mh j B j -B j (θ j + ) 2 -α B j ,m + B j ,m (θ j + ) 2 (B.18) Appendix C.
Recall that the displacement u z at the surface of a layered body is: Finding the exact solution by inverting the FT of the influence coefficients with respect to time would demand a significant amount of ressources. We propose an alternative way.

u z (x, y, t) =FT -1 {û z (x, y, ω)} = FT -1 { Ĵ1 (ω) Ĝ(x,
Basically, we need to find the function F(x, y, t) which can be written as:

F(x, y, t) = Ĝ(x, y, ω)iω p(ω) exp (iωt) dω (C.6)

In our work we have hypothesized that:

F(x, y, t) = FT -1 { Ĝ(x, y, Ĵk (ω) Ĵk+1 (ω)

)iω p(ω)} ≈ Ḡ x, y, J k (t) J k+1 (t)

∂ ∂t p(t). (C.7)
This assumption is not mathematically correct. However, it is exact under some conditions and we believe that it is a good approximation in the other cases of our studies. In particular, it is exact if the term Ĝ(x, y,

Ĵk (ω)
Ĵk+1 (ω) ) does not depend on ω, with regard to Eq. ((C.6)). In what follows next, we give details to justify this assumption.

We recall that we use the Maxwell generalised model for which the creep function of the layer k can be written as:

J k (t) = j k 0 - N k l=1 j k l exp -t/τ k l , t > 0. (C.8)
The Fourier transform of this function is given by: Ĵk (ω) = j k 0 δ(ω) - Ĵk+1 (ω) , it can be seen that for the times where we want to observe viscoelasticity, the ratio takes a constant value apart from the frequencies close to ω = 0. This phenomenon is explained by the behaviour of the Fourier transform: "Generally speaking, the more concentrated f (x) is, the more spread out its Fourier transform f (ξ) must be. In particular, the scaling property of the Fourier transform may be seen as saying: if we squeeze a function in x, its Fourier transform stretches out in ξ. It is not possible to arbitrarily concentrate both a function and its Fourier transform." (source: Wikipedia, "Uncertainty principle" in https://en.wikipedia.org/wiki/Fourier transform ) When we want to observe viscoelasticity, we need the regions where the creep functions spread. Therefore, their transforms are concentrated around the frequency 0. On the basis of this observation, we can write: This term corresponds to very low frequencies where we can physically assume that the viscoelastic effects are not important. Thus, those low frequencies are not expected to bring a significant effect on the viscoelastic response of the structure. In addition, the ignored points are in an area so small in the frequency domain that in a discretised form they might be missed without significant loss of informations. This explained the fact that our approximation allows to get good results when compared to FEM as shown in the previous paper (see [START_REF] Wallace | Three-dimensional rolling/sliding contact on a viscoelastic layered halfspace[END_REF]).

F(t) =

Conclusion:

. In summary, we give areas of validity of our approximation:

1. In steady state regime, the approximation meets the exact solution. Indeed, in the steady state regime, the behaviour in Ĝ(x, y, ω) is constant; and according to Eq. ((C.6)), it is the exact solution to write Eq. ((C.7)). Therefore, since all the final steps of our results are given for the steady state regime, they are exact.

2. When the ratio of creep functions is a constant, the ratio of their FT gives a constant. Thence, the term ε 0 Ĵ1 (ω) Ĝ(x, y, ω) exp (iωt) dω is equal to 0. Consequently, the applied solution is exact in the study of the effect of "Variation of the instantaneous shear moduli" in the paper (section 3.3.1).

3. For the other cases, the error made is contained in the term ε 0 Ĵ1 (ω) Ĝ(x, y, ω) exp (iωt) dω.

This term involves very low frequencies where the viscoleastic behaviour is not physically very significant. Thus, ignoring it implies an error which is not expected to by high.

Figure 1 :

 1 Figure 1: System being considered. A pressure and a shear distributions are applied at the surface of a viscoelastic multi-layered half-space, where the L layers and substrate have different viscoelastic behaviours

  relaxation time of one elementary model and H(t) the Heaviside function.

Figure 3 :

 3 Figure 3: Creep function (a) and relaxation function (b)

Figure 4 :

 4 Figure 4: Different cases of variation of the instantaneous shear moduli of the layers, with 10 layers.

Figure 5 :

 5 Figure 5: Normalized contact pressures along the rolling direction axis x, for a spherical contact on a multi-layered viscoelastic half-space with 10 layers with different instantaneous shear moduli; normal force F n = 1.48N prescribed and the dimensionless rolling velocity is vτ sub /a 0 = 0.8. The results are plotted at (a) t = 0, (b) t = τ sub /4, (c) t = τ sub , (d) t = 5 × τ sub .

  [START_REF] Koumi | Rolling contact of a rigid sphere/sliding of a spherical indenter upon a viscoelastic half-space containing an ellipsoidal inhomogeneity[END_REF] and Eqs.(22)(23)(24) in[START_REF] Wallace | Three-dimensional rolling/sliding contact on a viscoelastic layered halfspace[END_REF].

Figure 6 :Figure 7 :

 67 Figure 6: Normalized contact pressures along the rolling direction axis x, for a spherical contact on a multi-layered viscoelastic half-space with 10 layers with different instantaneous shear moduli; normal rigid displacement δ = 0.1a 0 prescribed and the dimensionless rolling velocity is vτ/a 0 = 0.8. The results are plotted at (a) t = 0, (b) t = τ sub /10, (c) t = τ sub /4, (d) t = τ sub .

Figure 8 :

 8 Figure 8: Normalized second invariant of the stress tensor √ J 2 /p 0 in the plane y = 0, for a spherical contact on a multilayered viscoelastic half-space with 10 layers with different instantaneous shear moduli; normal rigid displacement δ = 0.1a 0 prescribed and the dimensionless rolling velocity is vτ/a 0 = 0.8. Results are plotted at t = 0 (left) and t = 2τ (right)

Figure 9 :Figure 10 :

 910 Figure 9: Normalized contact pressures along the rolling direction, for a spherical contact on a multi-layered viscoelastic half-space with 10 layers with different relaxation times; normal rigid displacement δ = 0.1a 0 prescribed and the dimensionless rolling velocity is vτ/a 0 = 0.8. The results are plotted at (a) t = 0, (b) t = τ sub /10, (c) t = τ sub /4, (d) t = τ sub .

Figure 11 :Figure 12 :

 1112 Figure 11: Normalized second invariant of the stress tensor √ J 2 /p 0 in the plane y = 0, for a spherical contact on a multi-layered viscoelastic half-space with 10 layers with different relaxation times; normal rigid displacement δ = 0.1a 0 prescribed and the dimensionless rolling velocity is vτ/a 0 = 0.8. Results are plotted at t = 0 (left) and t = 2τ (right)

Figure 13 :

 13 Figure13: Apparent friction coefficient during time, for a spherical contact on a multi-layered viscoelastic half-space with 10 layers in the increasing modulus case coupled with the decreasing relaxation time case; the dimensionless rolling velocity is vτ/a 0 = 0.8. Results are plotted for an applied force (left) and for an applied body displacement (right)

  y, ω)iω p(ω)} (C.1) = G(x, y, t -ξ) ∂p(ξ) ∂ξ dξ. (C.2) = J 1 (t -ξ)F(x, y, ξ)dξ (C.3)whereF(x, y, t) = FT -1 { Ĝ(x, y, ω)iω p(ω)} (C.4) and G(x, y, t) = FT -1 { Ĵ1 (ω) Ĝ(x, y, ω)} = J 1 (t -η) Ḡ(x,y, η)dη. (C.5)

F

  y, ω)iω p(ω) exp (iωt) dω= ε 0 Ĝ(x, y, ω)iω p(ω) exp (iωt) dω + β ε Ĝ(x, y, ω)iω p(ω) exp (iωt) dω (C.10)where β is the equivalent of t in the frequency domain; and ε is a small frequency which indicates the frequencies over which the ratioĴk (ω) Ĵk+1 (ω) is not constant. Over [ε; β], the ratio Ĵk (ω)Ĵk+1 (ω) is constant. Thus, it makes sense to write: y, ω)iω p(ω) exp (iωt) dω + Ḡ x, ythe error in the calculation of F(x, y, t) lies in the term ε 0 Ĵ1 (ω) Ĝ(x, y, ω) exp (iωt) dω.

B j-1 θ j-1 --B j-1 θ j-1

Acknowledgements

The authors thank the Agence Nationale de la Recherche (ANR) for supporting this research work through the MACADAM project (ANR-17-CE22-0006-02). The authors would like to thank the reviewers for their insightful remarks.

Appendix B.

The analytical solution of Eq. ( ( 9)) is derived from the steps below. First let's use the following notations:

Then we can write the B functions as:

1v j t ( j) 0 + 1 t ( j) 0 + 1t ( j) 0 -1 (θ j -)