# Rolling contact on a viscoelastic multi-layered half-space

Efoé Rodrigue Wallace, Thibaut Chaise, Daniel Nélias

► To cite this version:

Efoé Rodrigue Wallace, Thibaut Chaise, Daniel Nélias. Rolling contact on a viscoelastic multi-layered half-space. International Journal of Solids and Structures, 2022, 239-240, 10.1016/j.ijsolstr.2021.111388. hal-03659660

HAL Id: hal-03659660

https://hal.science/hal-03659660

Submitted on 17 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Rolling Contact on a Viscoelastic Multi-layered Half-Space

Efoe Rodrigue WALLACEª, Thibaut CHAISEª, Daniel NELIAS<fnref n="1" />ª

ªUniv Lyon, INSA Lyon, CNRS, LaMCoS, UMR5259, 69621 Villeurbanne, France

## Abstract

This paper's purpose is the contact of a rolling body on a viscoelastic multi-layered half-space. Firstly, the contact is formulated with the classical definition of the gap and the pressure states in the contact zone and outside. L layers are considered on a substrate, L being any positive integer. Secondly, the influence coefficients related to an elastic multi-layered half-space are found, using the Papkovich-Neuber potentials. Traction and displacements continuity is assumed at the interfaces, and then a construction of matrix systems corresponding to a unit pressure and a unit shear imposed at the top surface, leads to the general problem to solve. Further, the solution for any kind of pressure and/or shear distributions at the contact surface is inferred by convolution with the influence coefficients found earlier, using the Fast Fourier Transform (FFT) algorithms. Throughout the steps, an Elastic/Viscoelastic correspondence is used in order to take into account not only the change of behaviour of the half-space during time, but also to superpose the load history. Conjugate Gradient Method (CGM) algorithms are used to solve the variational problem that yields from the contact problem definition. In the present paper, a parametric study is performed to highlight the effects of the elastic modulus and the relaxation time on some relatively simple cases of rolling contact.

**Key words:** Viscoelastic material, Multi-layered half-space, Semi-Analytical Method, Transient analysis.

<footnote n="1">¹corresponding author, daniel.nelias@insa-lyon.fr</footnote>
# 1. 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 Carter (1926), followed by a certain number of scientists (see for example Kalker (1967), Bentall and Johnson (1967), Johnson (1985), Nowell and Hills (1988)) among whom Kalker (1967, 1990, 1991) 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, Guler et al. (2012) developed a semi-analytical method for solving the tractive rolling contact but for graded coatings. Wang et al. (2012) developed a numerical approach for analysing three-dimensional steady-state rolling contact including creep using a semi-analytical method. Manyo et al. (2020) 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 King (1988) and followed by others (Nogi and Kato (1997), Plumet and Dubourg (1998), Polonsky and Keer (2000), Liu and Wang (2002), Stepanov and Torskaya (2018), ...). For two elastic layers on an elastic substrate, Chen (1971) 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 (Chen
and Bull (2009), Djabella and Arnell (1994), Gorishnyy et al. (2003), Kot (2012), Komvopoulos (1988), Assogba et al. (2020)). However, Yu et al. (2014) 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 Kuo (1969), Bufler (1971), Farrell (1972), Peltier (1974)) 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 Lee and Radok (1960) to recently Xu et al. (2020), Zhang et al. (2020), Nguyen and Hwu (2020), Wallace et al. (2020), Menga et al. (2021), through Kalker (1991), Goryacheva et al. (1995), Wayne Chen et al. (2011) Koumi et al. (2014, 2015), Putignano and Carbone (2014); Putignano et al. (2015, 2016), Menga et al. (2014, 2016), Stepanov and Torskaya (2016, 2018), Kusche (2017), Scaraggi and Comingio (2017), Putignano and Carbone (2018), Goryacheva and Miftakhova (2019), Torskaya and Stepanov (2019), and others. The modelling of viscoelastic contacts is now widely known and solved for various motions (rolling, sliding, reciprocating, fretting), various punches (in terms 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.

## 2. Formulation

<figure><img src="image_2.png" /><figcaption>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</figcaption></figure>

The system being considered consists of two bodies (1 and 2) in contact over a region $\Gamma_c$. 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_3^{(B_1+B_2)}(x, y, t)$ and, if there is some, the rigid body displacement $\delta(t)$ at every time $t$. Then, we have the following equation:

$$g(x, y, t) = g_i(x, y) + \delta(t) + u_3^{(B_1+B_2)}(x, y, t) \quad (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) \in \Gamma_c; \quad g(x, y, t) = 0 \quad \text{and} \quad p(x, y, t) > 0 \qquad (2a) $$

$$ (x, y) \notin \Gamma_c; \quad g(x, y, t) > 0 \quad \text{and} \quad p(x, y, t) = 0 \qquad (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 $\Gamma_c$. This reads as follow:

$$ F_N(t) = \int_{\Gamma_c(t)} p(x, y, t)dS \qquad (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.

## 2.1. 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 $\mu_j$, its Poisson's ratio is $\nu_j$ and its thickness is $h_j$ ($j = 1, ..., L$). For the substrate, the shear modulus is $\mu_{L+1}$ or $\mu_{sub}$ and the Poisson's ratio is $\nu_{L+1}$ or $\nu_{sub}$. We use the Papkovich-Neuber potentials $\phi^j$ and $\psi^j$ to express the elastic field. $\phi^j$ and $\psi^j$ are harmonic functions of $(x, y)$, and when no body forces are present, Malvern (1969) showed that $\psi^j$ can vanish. Then we can write displacements and stresses as functions of the potentials $\phi^j$, $\psi^j$ and $\psi^j$:

$$ 2\mu_j u_i^j = \phi_{,i}^j + x\psi_{,i}^j + z_j\psi_{,z,i}^j - (3 - 4\nu_j)\psi_i^j \qquad (4a) $$
$$ \sigma_{ik}^{j} = \phi_{,ik}^{j} - 2\nu_{j}(\psi_{x,x}^{j} + \psi_{z,z}^{j})\delta_{ik} - (1 - 2\nu_{j})(\psi_{i,k}^{j} + \psi_{k,i}^{j}) + x\psi_{x,ik}^{j} + z_{j}\psi_{z,ik}^{j} \quad (4b) $$

where the Einstein's notation is used and $j = 1, \dots, L + 1$.

In the above equations, the coma stand for the partial derivative (eg: $\phi_{,x} = \frac{\partial \phi}{\partial 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 \left[ \frac{\partial}{\partial x} \phi(x) \right] = im\tilde{\phi}(m) \quad (5a) $$

$$ FT_x[-ix\phi(x)] = \frac{\partial}{\partial m} \tilde{\phi}(m). \quad (5b) $$

Furthermore, we define the transformed Papkovich-Neuber potentials using bi-harmonic functions:

$$ \tilde{\phi}^j = A^j \exp(-\alpha z_j) + \bar{A}^j \exp(\alpha z_j) \quad (6a) $$

$$ \tilde{\psi}_x^j = B^j \exp(-\alpha z_j) + \bar{B}^j \exp(\alpha z_j) \quad (6b) $$

$$ \tilde{\psi}_z^j = C^j \exp(-\alpha z_j) + \bar{C}^j \exp(\alpha z_j). \quad (6c) $$

Therefore, we have linear equations in the Fourier frequency domain to solve.

Let's introduce the following notations: $w_+^j = e^{\alpha z_j}$ and $w_-^j = e^{-\alpha z_j}$, where $\alpha = \sqrt{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$, $\bar{A}^j$, $B^j$, $\bar{B}^j$, $C^j$ and $\bar{C}^j$ terms are unknown.
Note that $\bar{A}^{L+1} = 0$, $\bar{B}^{L+1} = 0$ and $\bar{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:

$$ \sigma_{zz}^{1}(x, y, 0) = -\delta_{p}(x, y), \quad \sigma_{zx}^{1}(x, y, 0) = -\delta_{q}(x, y), \quad \sigma_{zy}^{1}(x, y, 0) = 0. \quad (7) $$

where $\delta_p$ and $\delta_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:

$$ \begin{aligned} \sigma_{zx}^{j}(x, y, h_j) &= \sigma_{zx}^{j+1}(x, y, 0) & u_x^j(x, y, h_j) &= u_x^{j+1}(x, y, 0) \\ \sigma_{zy}^{j}(x, y, h_j) &= \sigma_{zy}^{j+1}(x, y, 0) & u_y^j(x, y, h_j) &= u_y^{j+1}(x, y, 0) \\ \sigma_{zz}^{j}(x, y, h_j) &= \sigma_{zz}^{j+1}(x, y, 0) & u_z^j(x, y, h_j) &= u_z^{j+1}(x, y, 0) \end{aligned} \quad (8) $$

Hence, applying double Fourier Transform (FT) to the above boundary conditions, a system of 6L + 3 equations with 6L + 3 unknowns comes. Following the steps of Yu et al. (2014), the system can be split into two matrix systems: one keeping the applied shear contribution and the other, the contribution of the applied pressure. Thus, the matrix systems to solve are:

$$ \left( \begin{array}{c} 2\alpha(1 - v_1) \\ \vdots \\ -4(1 - v_j)\theta_+^j \\ 2\alpha(1 - v_j)\theta_-^j \\ 0 \end{array} \right) \left( \begin{array}{c} -2\alpha(1 - v_1) \\ \vdots \\ -4(1 - v_j)\theta_+^j \\ -2\alpha(1 - v_j)\theta_-^j \\ 2\alpha(1 - v_j) \end{array} \right) = \left( \begin{array}{c} 0 \\ \vdots \\ 4(1 - v_{j+1})\mu^{(j,j+1)} \\ 4(1 - v_{j+1})\mu^{(j,j+1)} \\ 2\alpha(1 - v_{j+1}) \end{array} \right) \left( \begin{array}{c} B^1 \\ \vdots \\ \bar{B}^j \\ B^{j+1} \\ \vdots \\ \bar{B}^L \\ B^{L+1} \end{array} \right) = \left( \begin{array}{c} -\bar{q} \\ \vdots \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \end{array} \right) \quad (9) $$
and

$$
\begin{equation}
\begin{pmatrix}
\begin{array}{r@{\quad}l}
\begin{array}{cc} \alpha & \alpha \\ \alpha & -\alpha \end{array} &
\begin{array}{cc} 2(1-\nu_1) & -2(1-\nu_1) \\ 1-2\nu_1 & 1-2\nu_1 \end{array} \\
\multicolumn{2}{c}{\vdots} \\
X_{1j} & X_{2j} \\
\multicolumn{2}{c}{\vdots} \\
0 &
\end{array}
&
\begin{array}{@{}ccccccc@{}}
0 & & & & & & \\
\begin{array}{cccc} \alpha & \alpha(\theta_+^L)^2 & \alpha h_L & \alpha h_L(\theta_+^L)^2 & -\alpha \mu^{(L,L+1)} \theta_+^L & 0 & \\
\alpha & -\alpha(\theta_+^L)^2 & 3-4\nu_L + \alpha h_L & (3-4\nu_L + \alpha h_L)(\theta_+^L)^2 & -\alpha \mu^{(L,L+1)} \theta_+^L & -\mu^{(L,L+1)}(3-4\nu_{L+1})\theta_+^L & \\
\alpha & -\alpha(\theta_+^L)^2 & 1-2\nu_L + \alpha h_L & (1-2\nu_L - \alpha h_L)(\theta_+^L)^2 & -\alpha \theta_+^L & -(1-2\nu_{L+1})\theta_+^L & \\
\alpha & \alpha(\theta_+^L)^2 & 2(1-\nu_L) + \alpha h_L & -[2(1-\nu_L) - \alpha h_L](\theta_+^L)^2 & -\alpha \theta_+^L & -2(1-\nu_{L+1})\theta_+^L &
\end{array}
&
\begin{pmatrix} A^1 \\ \bar{A}^1 \\ \vdots \\ C^j \\ \bar{C}^j \\ A^{j+1} \\ \bar{A}^{j+1} \\ \vdots \\ C^L \\ \bar{C}^L \\ A^{L+1} \\ \bar{A}^{L+1} \\ \vdots \\ C^{L+1} \end{pmatrix}
=
\begin{pmatrix} S_1^0 \\ S_2^0 \\ \vdots \\ S_1^j \\ S_2^j \\ S_3^j \\ S_4^j \\ \vdots \\ S_1^L \\ S_2^L \\ S_3^L \\ S_4^L \end{pmatrix}
\end{pmatrix}
\end{equation}
$$

where $\theta_+^j = e^{\alpha h_j}$, $\theta_-^j = e^{-\alpha h_j}$, $\mu^{(j,j+1)} = \frac{\mu_j}{\mu_{j+1}}$ and the sub-matrices $X_{1j}$ and $X_{2j}$ are written as follow:

$$
X_{1j} = \begin{pmatrix} \alpha & \alpha(\theta_+^j)^2 & \alpha h_j & \alpha h_j(\theta_+^j)^2 \\ \alpha & -\alpha(\theta_+^j)^2 & 3-4\nu_j + \alpha h_j & (3-4\nu_j - \alpha h_j)(\theta_+^j)^2 \\ \alpha & -\alpha(\theta_+^j)^2 & 1-2\nu_j + \alpha h_j & (1-2\nu_j - \alpha h_j)(\theta_+^j)^2 \\ \alpha & \alpha(\theta_+^j)^2 & 2(1-\nu_j) + \alpha h_j & -[2(1-\nu_j) - \alpha h_j](\theta_+^j)^2 \end{pmatrix} \qquad X_{2j} = \begin{pmatrix} -\alpha\mu^{(j,j+1)}\theta_+^j & -\alpha\mu^{(j,j+1)}\theta_+^j & 0 & 0 \\ -\alpha\mu^{(j,j+1)}\theta_+^j & \alpha\mu^{(j,j+1)}\theta_+^j & -\mu^{(j,j+1)}(3-4\nu_{j+1})\theta_+^j & -\mu^{(j,j+1)}(3-4\nu_{j+1})\theta_+^j \\ -\alpha\theta_+^j & \alpha\theta_+^j & -(1-2\nu_{j+1})\theta_+^j & -(1-2\nu_{j+1})\theta_+^j \\ -\alpha\theta_+^j & -\alpha\theta_+^j & -2(1-\nu_{j+1})\theta_+^j & 2(1-\nu_{j+1})\theta_+^j \end{pmatrix} \tag{11}
$$

In Eqs. ((10)), the right hand terms remain to be defined. Those $S_1^0, \dots, S_4^L$ are functions of the $B^j$ and $\bar{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, \bar{A}^j, B^j, \bar{B}^j, C^j$ and $\bar{C}^j$. Those explicit solutions have been given by Yu et al. (2014). 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) \times (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 Liu and Wang (2002)). In the present work, the IFT is handled with Fast Fourier Transform (FFT) algorithms as proposed by Liu et al. (2000).

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.

## 2.2. Viscoelastic rolling contact

<figure><img src="image_3.png" /><figcaption>Figure 2: Generalized Maxwell model</figcaption></figure>

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) = \left[ \mu_0 + \sum_{i=0}^{n} \mu_i \exp(-t/\tau_i) \right] H(t) \quad (12)$$

where,

$\mu_i$ is the spring stiffness,

$\eta_i$ the dashpot viscosity,

$\tau_i = \eta_i/\mu_i$ the relaxation time of one elementary model and

$H(t)$ the Heaviside function.

The Relaxation function in Eq. (12) is related to the creep function with the relation:

$$\int_{0}^{t} J(\xi)R(t-\xi)d\xi = t. \quad (13)$$

For one relaxation time, the creep and relaxation functions can directly be written as:

$$J(t) = \left[ \frac{1}{\mu_0} - \frac{1}{\mu_1} \exp\left(-\frac{t}{\tau}\right) \right] = \left[ \frac{1}{\mu_\infty} + \frac{1}{\mu_1} \left( 1 - \exp\left(-\frac{t}{\tau}\right) \right) \right] \quad (14)$$
$$R(t) = \left[ \mu_0 + (\mu_{\infty} - \mu_0) \exp \left( - \frac{t}{\frac{\mu_0}{\mu_{\infty}} \tau} \right) \right] \quad (15)$$

<figure><img src="image_4.png" /><figcaption>Figure 3: Creep function (a) and relaxation function (b)</figcaption></figure>

Further, we will use this model (see Fig. 3) to compute simple configurations in section 3.

### 2.2.1. Elastic/Viscoelastic correspondence

The Elastic/Viscoelastic correspondence that we use is the same as explained in Wallace et al. (2020). We first discretize the computation time interval into regular time steps. Then, at a time step $t$ we (i) replace the elastic modulus $\frac{1}{2\mu_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:

$$ \begin{aligned} u_i(x, y, t) &= \int_0^t \int \int_{\mathbb{R}^2} \mathcal{F}_{ij}(x-x', y-y', t-t') \dot{\sigma}_j(x', y', t') dx'dy'dt' \\ &= \int_0^t \int \int_{\mathbb{R}^2} \mathcal{F}_{ij}(x', y', t-t') \dot{\sigma}_j(x-x', y-y', t') dx'dy'dt' \end{aligned} \quad (16) $$

where $\sigma_j$ are the components of the applied stress at the surface and $\mathcal{F}_{ij}(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 $\mathcal{F}_{ij}(x, y, t)$ can
be written as:

$$ \mathcal{F}_{ij}(x, y, t) = J(t)G_{ij}(x, y), \qquad (17) $$

where the term $G_{ij}$ 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) = \int_0^t J(t-t') \int_{\mathbb{R}^2} G_{ij}(x', y') \dot{\sigma}_j(x-x', y-y', t') dx'dy'dt'. \qquad (18) $$

This means that a convolution is performed between the elastic influence coefficients $G_{ij}$ - 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\mu^r u_i^r = \phi_{,i}^r + x\psi_{,x,i}^r + z\psi_{,z,i}^r - (3 - 4\nu)\psi_i^r, \qquad (19) $$

where $\mu$ is the shear modulus of the layer $r$ and $\phi$ and $\psi_i$ are the potentials of Papkovich-Neuber.

In the left hand term, the factor $2\mu$ 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 $\mu^k/\mu^{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 $\mu^k/\mu^{k+1}$.
◇ The last remark have lead us to propose the following equation for the viscoelastic displacements:

$$u_i(x, y, t) = \int_0^t J(t-t') \int_{\mathbb{R}^2} \bar{G}_{ij}^1(x', y', J^k(t')/J^{k+1}(t')) \dot{\sigma}_j(x-x', y-y', t') dx'dy'dt'. \quad (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. (11) is:

$$u_z^{B2}(x, y, t) = \int_0^t J^1(t-\xi) \left[ \bar{G}^1(x, y, \xi) * \frac{\partial p(x, y, \xi)}{\partial \xi} \right] d\xi \quad (21)$$

where,

* * represents the convolution product,

* $J^1(t)$ is the creep function of the top layer,

* $\bar{G}^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 $\mu^k/\mu^{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.
### 2.2.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 (Koumi et al. (2015)). 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 $\Delta t$, the elementary displacement of the body 1 is then $\Delta \vec{r} = \vec{v} \times \Delta 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_z^{B2}(x, y, t) = \int_0^t J^1(t - \xi) \left[ \bar{G}^1(x, y, \xi) * \frac{\partial p(x^0 - \alpha \Delta r_x^0, y^0 - \alpha \Delta r_y^0, \xi)}{\partial \xi} \right] d\xi. \quad (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 $\alpha \Delta \vec{r}^0(\alpha \Delta r_x^0, \alpha \Delta r_y^0)$, $\alpha$ being the actual time step index.

### 2.2.3. Numerical approach

The numerical implementation needs a discretization of both space and time. The contact surface is discretized into $N_1 \times N_2$ elements, while the time domain is subdivided into $N_t$ time steps $\Delta t$. Then, using the finite difference method for $\frac{\partial p}{\partial \xi}$ the discretized form of Eq. (22) becomes:

$$u_z^{B2}(i, j, \alpha) = \sum_{k=0}^{\alpha} J^1[(\alpha - k)\Delta t] \times \left[ \sum_{i'=1}^{N_1} \sum_{j'=1}^{N_2} \bar{G}^1(i - i', j - j', k\Delta t) \times (p(i', j', k) - p(i', j', k - 1)) \right], \quad (23)$$

where one can infer that $i, j$ and $\alpha$ 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 Gallego et al. (2010) 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 \times N_2)^2$ to $(N_1 \times N_2 \times \log(N_1 \times N_2))$ (see the algorithm of Liu et al. (2000)).

At the time step $\alpha$, the only unknown in the above equation is the actual pressure field. The implementation gives:

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

2. For any $\alpha$, the Green functions are recalculated $\tilde{G}^1(i - i', j - j', (\alpha)\Delta t)$. Then, the term

$$ \mathcal{H}(\alpha) = \sum_{k=0}^{\alpha-1} J^1[(\alpha - k)\Delta t] \times \left[ \sum_{i'=1}^{N_1} \sum_{j'=1}^{N_2} \tilde{G}^1(i - i', j - j', k\Delta t) \times (p(i', j', k) - p(i', j', k-1)) \right] \\ - J^1[(\alpha - \alpha)\Delta t] \times \left[ \sum_{i'=1}^{N_1} \sum_{j'=1}^{N_2} \tilde{G}^1(i - i', j - j', (\alpha)\Delta t) \times (p(i', j', \alpha-1)) \right] \quad (24) $$

is calculated and stored.

The remaining term $J^1[(\alpha - \alpha)\Delta t] \times [\sum_{i'=1}^{N_1} \sum_{j'=1}^{N_2} \tilde{G}^1(i - i', j - j', (\alpha)\Delta t) \times (p(i', j', \alpha))]$ contains the unknown pressure field and so is solved with the CGM algorithm.

When the surface pressure is known at the time step $\alpha$, 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.

# 3. Results

## 3.1. 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 $\delta$ on a sphere of radius $R = 10$mm 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 $\mu_{sphere} = 10^{20}$MPa and its Poisson's ratio is $\nu_{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, \dots, 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 $\mu_{\infty j}/\mu_{0j} = 10$ ($j = 1, \dots, L + 1$), and we do not change the substrate mechanical characteristics which are: $\mu_{\infty sub} = \mu_{\infty 11} = 3.86$MPa and $\tau_{sub} = \tau_{11} = 10$s. 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 $\mu_{\infty j}$ or/and its relaxation time $\tau_j$ ($j = 1, \dots, 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 \times \mu_{\infty sub}$ while the odd index layers take the shear modulus value $\mu_{\infty sub}$:

$$ \mu_{\infty j} = \frac{\mu_{\infty j+1}}{4} = \mu_{\infty sub} \quad (j = 1, 3, 5, 7, 9). \qquad (25) $$
<figure><img src="image_5.png" /><figcaption>Figure 4: Different cases of variation of the instantaneous shear moduli of the layers, with 10 layers.</figcaption></figure>

• The "Decreasing Modulus": the modulus decreases from $6 \times \mu_{\infty sub}$ with a constant step to reach the substrate modulus $\mu_{\infty sub}$:

$$ \mu_{\infty j} = \left[ 6 - \frac{j-1}{2} \right] \times \mu_{\infty sub} \quad (j = 1, \dots, 10). \qquad (26) $$

• The "Increasing Modulus": the modulus increases from $\mu_{\infty sub}$ to $5.5 \times \mu_{\infty sub}$ in the 10th layer:

$$ \mu_{\infty j} = \left[ 1 + \frac{j-1}{2} \right] \times \mu_{\infty sub} \quad (j = 1, \dots, 10). \qquad (27) $$

We can replicate the cases above to have similar changes of the relaxation times. Then we have, $\tau_j = \frac{\tau_{j+1}}{2} = \tau_{sub}$ (with $j = 1, 3, 5, 7, 9$) for the "Alternate Relaxation" time; $\tau_j = [3 - \frac{j-1}{5}] \times \tau_{sub}$ (with $j = 1, \dots, 10$) for the "Decreasing Relaxation" time and $\tau_j = [1 + \frac{j-1}{5}] \times \tau_{sub}$ (with $j = 1, \dots, 10$) for the "Increasing Relaxation" time.

## 3.2. 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\tau_{sub}/a_0 = 0.8$. Then, we compute in the time domain $[0, 6\tau_{sub}]$ subdivided into 600 time steps. The contact surface is discretized with elements of size $0.2a_0 \times 0.2a_0$.
<figure><img src="image_6.png" /><figcaption>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}$.</figcaption></figure>
Figure 5 shows the evolution of the contact pressure along the rolling direction axis x. Results are plotted at (a) $t = 0$, (b) $t = \tau_{sub}/4$, (c) $t = \tau_{sub}$, (d) $t = 5 \times \tau_{sub}$ for the homogeneous half-space 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 Yu et al. (2014) 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 Wallace et al. (2020), 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 = \tau_{sub}/4$. At the time $t = \tau_{sub}$, the pressure presents the shape of a steady-state regime but will continue decreasing to reach the steady-state regime ($t = 5\tau_{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 Koumi et al. (2015). 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.

### 3.3. Applied Normal Rigid Body Displacement

In this section, we prescribe a normal displacement $\delta = 0.1a_0$ and compute in the time domain $[0, 2\tau_{sub}]$ with 80 time steps. We discretize the space with elements of size $0.2a_0 \times 0.2a_0 \times 0.05a_0$. We still use the same rigid sphere, and the rolling velocity is still given by $v\tau_{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.

#### 3.3.1. Variation of the instantaneous shear moduli

Despite the absence of Coulomb's friction coefficient, the viscoelastic dissipation introduces an asymmetry of the contact pressure distribution along the rolling direction, which in fine produces
a resisting force or a resisting torque. This phenomenon is also encountered for frictionless elastic-plastic rolling/sliding contact, as shown in Chaise and Nelias (2011). For more details regarding frictionless rolling/sliding viscoelastic contact, the reader may refer to Fig. 3 in Koumi et al. (2015) and Eqs. (22-24) in Wallace et al. (2020).

What is called here the apparent friction coefficient $\mu_{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) = \int_{\Gamma_c(t)} \frac{\partial u_3(x_1, x_2, t)}{\partial x_1} p(x_1, x_2, t) dx_1 dx_2 \quad (28)$$

where $\Gamma_c(t)$ is the real contact area at time t.

Then the apparent friction coefficient is $\mu_{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 $\tau_j = \tau_{sub}$ ($j = 1, \dots, 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 = \tau_{sub}/10$, (c) $t = \tau_{sub}/4$ and (d) $t = \tau_{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), Wang et al. (2012), and Wallace et al. (2020), 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
<figure><img src="image_7.png" /><figcaption>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 $\delta = 0.1a_0$ prescribed and the dimensionless rolling velocity is $v\tau/a_0 = 0.8$. The results are plotted at (a) $t = 0$, (b) $t = \tau_{sub}/10$, (c) $t = \tau_{sub}/4$, (d) $t = \tau_{sub}$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Apparent friction coefficient, for a spherical contact on a multi-layered viscoelastic half-space with 10 layers with different instantaneous shear moduli; normal rigid displacement $\delta = 0.1a_0$ prescribed and the dimensionless rolling velocity is given by $v\tau/a_0 = 0.8$</figcaption></figure>

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 one of the increasing modulus case.

Further, when the time increases, the trend will stay the same since all the layers have the same relaxation time. In addition, we can see the evolution through a transient regime to achieve what seems to be a steady-state regime 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\tau_{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 $\sqrt{J_2}/p_0$, at $t = 0$ (left) which corresponds to the elastic case, and at $t = 2\tau_{sub}$ which is in the steady-state regime.
<figure><img src="image_9.png" /><figcaption>Figure 8: Normalized second invariant of the stress tensor $\sqrt{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 instantaneous shear moduli; normal rigid displacement $\delta = 0.1a_0$ prescribed and the dimensionless rolling velocity is $v\tau/a_0 = 0.8$. Results are plotted at $t = 0$ (left) and $t = 2\tau$ (right)</figcaption></figure>
Again, the 3D-stresses distribution correlates well with the results of Yu et al. (2014) at $t = 0$. One can clearly see that the effect of the alternate modulus induces 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).

### 3.3.2. 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 $\mu_{\infty j} = \mu_{\infty sub}$ ($j = 1, \dots, 10$) and modify the relaxation times following the cases described 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 relaxation time. Thus, the contact pressure relaxes less rapidly. And, even when the steady-state regime is achieved, the pressure profile is less wide and the maximum value
<figure><img src="image_10.png" /><figcaption>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 $\delta = 0.1a_0$ prescribed and the dimensionless rolling velocity is $v\tau/a_0 = 0.8$. The results are plotted at (a) $t = 0$, (b) $t = \tau_{sub}/10$, (c) $t = \tau_{sub}/4$, (d) $t = \tau_{sub}$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Apparent friction coefficient during time, for a spherical contact on a multi-layered viscoelastic half-space with 10 layers with different relaxation times; normal rigid displacement $\delta = 0.1a_0$ prescribed and the dimensionless rolling velocity is given by $v\tau/a_0 = 0.8$</figcaption></figure>

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\tau$, 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
<figure><img src="image_12.png" /><figcaption>Figure 11: Normalized second invariant of the stress tensor $\sqrt{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 $\delta = 0.1a_0$ prescribed and the dimensionless rolling velocity is $v\tau/a_0 = 0.8$. Results are plotted at $t = 0$ (left) and $t = 2\tau$ (right)</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: Normalized contact pressure along the rolling direction axis x, 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)</figcaption></figure>

leading edge of the surface and closer to the surface.

### 3.3.3. 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 \times 0.2a_0$. And, we perform the computations in the time domain $[0, 6\tau_{sub}]$ for the applied force case and in $[0, 2\tau_{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
<figure><img src="image_14.png" /><figcaption>Figure 13: 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)</figcaption></figure>

this is given with the apparent friction coefficient plotted in Fig. 13. The latter figure shows also that for an applied force, the resistance to rolling is higher than for an applied body displacement.

4. Conclusion

In summary, for the first 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 ($\approx 2h15min$) on a personal 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.

## 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 also like to thank the reviewers for their insightful remarks.
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\mu_j \tilde{u}_x^j = im(D^j w_-^j + \bar{D}^j w_+^j) - 4(1 - v_j)(B^j w_-^j + \bar{B}^j w_+^j) + m^2 z_j \alpha^{-1} (B^j w_-^j - \bar{B}^j w_+^j) + imz_j (C^j w_-^j + \bar{C}^j w_+^j) \quad (A.1a)$$

$$2\mu_j \tilde{u}_y^j = in(D^j w_-^j + \bar{D}^j w_+^j) - mnz_j \alpha^{-1} (B^j w_-^j - \bar{B}^j w_+^j) + inz_j (C^j w_-^j + \bar{C}^j w_+^j) \quad (A.1b)$$

$$2\mu_j \tilde{u}_z^j = -\alpha(D^j w_-^j - \bar{D}^j w_+^j) - im\alpha^{-1}(B^j w_-^j - \bar{B}^j w_+^j) - (3 - 4v_j)(C^j w_-^j + \bar{C}^j w_+^j) + imz_j(B^j w_-^j + \bar{B}^j w_+^j) + z_j\alpha(-C^j w_-^j + \bar{C}^j w_+^j) \quad (A.1c)$$

$$\tilde{\sigma}_{xx}^j = -m^2(D^j w_-^j + \bar{D}^j w_+^j) - 2im(2 - v_j)(B^j w_-^j + \bar{B}^j w_+^j) - 2\alpha v_j(-C^j w_-^j + \bar{C}^j w_+^j) + im^3 z_j \alpha^{-1}(B^j w_-^j - \bar{B}^j w_+^j) - z_j m^2(C^j w_-^j + \bar{C}^j w_+^j) \quad (A.1d)$$

$$\tilde{\sigma}_{yy}^j = -n^2(D^j w_-^j + \bar{D}^j w_+^j) - 2imv_j(B^j w_-^j + \bar{B}^j w_+^j) - 2\alpha v_j(-C^j w_-^j + \bar{C}^j w_+^j) + imn^2 z_j \alpha^{-1}(B^j w_-^j - \bar{B}^j w_+^j) - z_j m^2(C^j w_-^j + \bar{C}^j w_+^j) \quad (A.1e)$$

$$\tilde{\sigma}_{zz}^j = \alpha^2(D^j w_-^j + \bar{D}^j w_+^j) + 2im(1 - v_j)(B^j w_-^j + \bar{B}^j w_+^j) - 2\alpha(1 - v_j)(-C^j w_-^j + \bar{C}^j w_+^j) - imz_j \alpha(B^j w_-^j - \bar{B}^j w_+^j) + z_j \alpha^2(C^j w_-^j + \bar{C}^j w_+^j) \quad (A.1f)$$

$$\tilde{\sigma}_{xy}^j = -mn(D^j w_-^j + \bar{D}^j w_+^j) - 2in(1 - v_j)(B^j w_-^j + \bar{B}^j w_+^j) + im^2 nz_j \alpha^{-1}(B^j w_-^j - \bar{B}^j w_+^j) - z_j mn(C^j w_-^j + \bar{C}^j w_+^j) \quad (A.1g)$$

$$\tilde{\sigma}_{xz}^j = -im\alpha(D^j w_-^j - \bar{D}^j w_+^j) + (2\alpha(1 - v_j) + m^2\alpha^{-1}) \times (B^j w_-^j - \bar{B}^j w_+^j) - im(1 - 2v_j)(C^j w_-^j + \bar{C}^j w_+^j) - m^2 z_j(B^j w_-^j + \bar{B}^j w_+^j) + iz_j \alpha m(-C^j w_-^j + \bar{C}^j w_+^j) \quad (A.1h)$$

$$\tilde{\sigma}_{yz}^j = -in\alpha(D^j w_-^j - \bar{D}^j w_+^j) + mn\alpha^{-1}(B^j w_-^j - \bar{B}^j w_+^j) - in(1 - 2v_j)(C^j w_-^j + \bar{C}^j w_+^j) - mnz_j(B^j w_-^j + \bar{B}^j w_+^j) + inz_j\alpha(-C^j w_-^j + \bar{C}^j w_+^j) \quad (A.1i)$$
## Appendix B.

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

$$t_0^{(L+1)} = 1 \qquad (B.1a)$$

$$t_0^{(j)} = \mu^{(j,j+1)} \frac{t_0^{(j+1)} + 1 + [t_0^{(j+1)} - 1](\theta_-^{j+1})^2}{t_0^{(j+1)} + 1 - [t_0^{(j+1)} - 1](\theta_-^{j+1})^2} \quad (j = L, \dots, 1) \qquad (B.1b)$$

Then we can write the B functions as:

$$\bar{B}^1 = \frac{(\theta_-^1)^2}{t_0^{(1)} + 1 - [t_0^{(1)} - 1](\theta_-^1)^2} \frac{-1}{2\alpha(1 - \nu_1)} [t_0^{(1)} - 1] \qquad (B.2)$$

$$B^1 = \bar{B}^1 - \frac{1}{2\alpha(1 - \nu_1)} \qquad (B.3)$$

$$\bar{B}^j = \frac{1 - \nu_{j-1}}{1 - \nu_j} \frac{[t_0^{(j)} - 1](\theta_-^j)^2}{t_0^{(j)} + 1 - [t_0^{(j)} - 1](\theta_-^j)^2} [B^{j-1}\theta_-^{j-1} - \bar{B}^{j-1}\theta_+^{j-1}] \quad (j = 2, \dots, L) \qquad (B.4)$$

$$B^j = \frac{1 - \nu_{j-1}}{1 - \nu_j} \frac{t_0^{(j)} + 1}{t_0^{(j)} + 1 - [t_0^{(j)} - 1](\theta_-^j)^2} [B^{j-1}\theta_-^{j-1} - \bar{B}^{j-1}\theta_+^{j-1}] \quad (j = 2, \dots, L) \qquad (B.5)$$

$$B^{L+1} = \frac{1 - \nu_L}{1 - \nu_{L+1}} [B^L \theta_-^L - \bar{B}^L \theta_+^L] \qquad (B.6)$$

To write the second members of Eq. ((10)), we need the derivatives of the above functions. These derivatives about $m$ are:

$$t_{0,m}^{(L)} = 0 \qquad (B.7a)$$

$$t_{0,m}^{(j-1)} = 4\mu^{(j-1,j)} (\theta_-^j)^2 \frac{t_{0,m}^{(j)} - [t_0^{(j)} + 1][t_0^{(j)} - 1]^{\frac{mh_j}{\alpha}}}{\{t_0^{(j)} + 1 - [t_0^{(j)} - 1](\theta_-^j)^2\}^2} \quad (j = L, \dots, 2) \qquad (B.7b)$$

$$\bar{B}_{,m}^1 = \frac{(\theta_-^1)^2}{t_0^{(1)} + 1 - (t_0^{(1)} - 1)(\theta_-^1)^2} \frac{-1}{2\alpha(1 - \nu_1)} \left\{ \left[ \frac{2t_{0,m}^{(1)} - \left(\frac{2mh_1}{\alpha}\right)(t_0^{(1)} + 1)(t_0^{(1)} - 1)}{(t_0^{(1)} + 1) - (t_0^{(1)} - 1)(\theta_-^1)^2} \right] - \frac{m}{\alpha^2} \right\} \qquad (B.8)$$
$$B_{,m}^{1} = \bar{B}_{,m}^{1} + \frac{m}{2\alpha^{3}(1 - \nu_{1})} \qquad (B.9)$$

$$\bar{B}_{,m}^{j} = \frac{1 - \nu_{j-1}}{1 - \nu_{j}} \frac{[t_{0}^{(j)} - 1](\theta_{-}^{j})^{2}}{t_{0}^{(j)} + 1 - [t_{0}^{(j)} - 1](\theta_{-}^{j})^{2}} \left\{ \frac{-2mh_{j}}{\alpha} \frac{[t_{0}^{(j)} + 1] + \frac{2t_{0m}^{(j)}}{t_{0}^{(j)-1}}}{t_{0}^{(j)} + 1 - [t_{0}^{(j)} - 1](\theta_{-}^{j})^{2}} [\bar{B}^{j-1}\theta_{-}^{j-1} - \bar{B}^{j-1}\theta_{+}^{j-1}] \right. \\ \left. + B_{,m}^{j-1}\theta_{-}^{j-1} - \bar{B}_{,m}^{j-1}\theta_{+}^{j-1} - \frac{mh_{j-1}}{\alpha} [\bar{B}^{j-1}\theta_{-}^{j-1} + \bar{B}^{j-1}\theta_{+}^{j-1}] \right\} \quad (j = 2, \dots, L) \qquad (B.10)$$

$$B_{,m}^{j} = \frac{1 - \nu_{j-1}}{1 - \nu_{j}} \frac{t_{0}^{(j)} + 1}{t_{0}^{(j)} + 1 - [t_{0}^{(j)} - 1](\theta_{-}^{j})^{2}} \left\{ \frac{\frac{2t_{0m}^{(j)}}{t_{0}^{(j)+1}}(\theta_{-}^{j})^{2} + [t_{0}^{(j)} - 1]\frac{-2mh_{j}}{\alpha}(\theta_{-}^{j})^{2}}{t_{0}^{(j)} + 1 - [t_{0}^{(j)} - 1](\theta_{-}^{j})^{2}} [\bar{B}^{j-1}\theta_{-}^{j-1} - \bar{B}^{j-1}\theta_{+}^{j-1}] \right. \\ \left. + B_{,m}^{(j-1)}\theta_{-}^{j-1} - \bar{B}_{,m}^{(j-1)}\theta_{+}^{j-1} - \frac{mh_{j-1}}{\alpha} [\bar{B}^{j-1}\theta_{-}^{j-1} + \bar{B}^{j-1}\theta_{+}^{j-1}] \right\} \quad (j = 2, \dots, L) \qquad (B.11)$$

$$B_{,m}^{L+1} = \frac{1 - \nu_{L}}{1 - \nu_{L+1}} \left[ B_{,m}^{L}\theta_{-}^{L} - \frac{mh_{L}}{\alpha}\theta_{-}^{L}B^{L} - \bar{B}_{,m}^{L}\theta_{+}^{L} - \frac{mh_{L}}{\alpha}\theta_{+}^{L}\bar{B}^{L} \right] \qquad (B.12)$$

Here, we can express the functions of the second member of Eq. ((10)):

$$S_{1}^{(0)} = -\alpha^{-1} - 2im\alpha^{-1}(1 - \nu_{1})[B^{1} + \bar{B}^{1}] - i\alpha[B_{,m}^{1} + \bar{B}_{,m}^{1}] \qquad (B.13)$$

$$S_{2}^{(0)} = -i \left\{ m\alpha^{-1} [B^{1} - \bar{B}^{1}] + \alpha [B_{,m}^{1} - \bar{B}_{,m}^{1}] \right\} = 0 \qquad (B.14)$$

and for $j = 1, \dots, L$:

$$S_{1}^{(j)} = \alpha i \left\{ \mu^{(j,j+1)} [B_{,m}^{j+1}\theta_{+}^{j} + \bar{B}_{,m}^{j+1}\theta_{+}^{j}] + \frac{1}{\alpha} mh_{j} [B^{j} - \bar{B}^{j}(\theta_{+}^{j})^{2}] - B_{,m}^{j} - \bar{B}_{,m}^{j}(\theta_{+}^{j})^{2} \right\} \qquad (B.15)$$

$$S_{2}^{(j)} = i \left\{ \frac{\mu^{(j,j+1)}m}{\alpha} [B^{j+1}\theta_{+}^{j} - \bar{B}^{j+1}\theta_{+}^{j}] + \mu^{(j,j+1)}\alpha [B_{,m}^{j+1}\theta_{+}^{j} - \bar{B}_{,m}^{j+1}\theta_{+}^{j}] \right. \\ \left. + mh_{j} [B^{j} + \bar{B}^{j}(\theta_{+}^{j})^{2}] - \frac{m}{\alpha} [B^{j} - \bar{B}^{j}(\theta_{+}^{j})^{2}] - \alpha [B_{,m}^{j} - \bar{B}_{,m}^{j}(\theta_{+}^{j})^{2}] \right\} \qquad (B.16)$$

$$S_{3}^{(j)} = i \left\{ \alpha^{-1}m [B^{j+1}\theta_{+}^{j} - \bar{B}^{j+1}\theta_{+}^{j}] + \alpha [B_{,m}^{j+1}\theta_{+}^{j} - \bar{B}_{,m}^{j+1}\theta_{+}^{j}] \right. \\ \left. - (\alpha^{-1} - h_{j})mB^{j} + (\alpha^{-1} + h_{j})m(\theta_{+}^{j})^{2}\bar{B}^{j} - \alpha [B_{,m}^{j} - \bar{B}_{,m}^{j}(\theta_{+}^{j})^{2}] \right\} \qquad (B.17)$$
$$S_4^{(j)} = i \left\{ \begin{aligned} & 2m\alpha^{-1}(1 - v_{j+1}) \left[ B^{j+1}\theta_+^j + \bar{B}^{j+1}\theta_+^j \right] + \alpha \left[ B_{,m}^{j+1}\theta_+^j + \bar{B}_{,m}^{j+1}\theta_+^j \right] \\ & - 2m\alpha^{-1}(1 - v_j) \left[ B^j + \bar{B}^j(\theta_+^j)^2 \right] + mh_j \left[ B^j - \bar{B}^j(\theta_+^j)^2 \right] - \alpha \left[ B_{,m}^j + \bar{B}_{,m}^j(\theta_+^j)^2 \right] \end{aligned} \right\} \quad (B.18)$$

## Appendix C.

Recall that the displacement $u_z$ at the surface of a layered body is:

$$u_z(x, y, t) = FT^{-1}\{\hat{u}_z(x, y, \omega)\} = FT^{-1}\{\hat{J}^1(\omega)\hat{\bar{G}}(x, y, \omega)i\omega\hat{p}(\omega)\} \quad (C.1)$$

$$= \int G(x, y, t - \xi) \frac{\partial p(\xi)}{\partial \xi} d\xi. \quad (C.2)$$

$$= \int J^1(t - \xi) F(x, y, \xi) d\xi \quad (C.3)$$

where

$$F(x, y, t) = FT^{-1}\{\hat{\bar{G}}(x, y, \omega)i\omega\hat{p}(\omega)\} \quad (C.4)$$

and

$$G(x, y, t) = FT^{-1}\{\hat{J}^1(\omega)\hat{\bar{G}}(x, y, \omega)\} = \int J^1(t - \eta)\bar{G}(x, y, \eta)d\eta. \quad (C.5)$$

Finding the exact solution by inverting the FT of the influence coefficients with respect to time would demand a significant amount of resources. We propose an alternative way.

Basically, we need to find the function $F(x, y, t)$ which can be written as:

$$F(x, y, t) = \int \hat{\bar{G}}(x, y, \omega)i\omega\hat{p}(\omega) \exp(i\omega t) d\omega \quad (C.6)$$

In our work we have hypothesized that:

$$F(x, y, t) = FT^{-1}\{\hat{\bar{G}}(x, y, \frac{\hat{J}^k(\omega)}{\hat{J}^{k+1}(\omega)})i\omega\hat{p}(\omega)\} \approx \bar{G}\left(x, y, \frac{\hat{J}^k(t)}{\hat{J}^{k+1}(t)}\right) \frac{\partial}{\partial t}p(t). \quad (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 $\hat{\bar{G}}(x, y, \frac{\hat{J}^k(\omega)}{\hat{J}^{k+1}(\omega)})$ does not depend on $\omega$, 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_0^k - \sum_{l=1}^{N_k} j_l^k \exp(-t/\tau_l^k), \quad t > 0. \qquad (C.8)$$

The Fourier transform of this function is given by:

$$\hat{J}^k(\omega) = j_0^k \delta(\omega) - \sum_{l=1}^{N_k} \frac{j_l^k \tau_l^k}{1 + i\omega\tau_l^k}. \qquad (C.9)$$

If one plots the norm of the ratio $\frac{\hat{J}^k(\omega)}{\hat{J}^{k+1}(\omega)}$, 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 $\omega = 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 $\hat{f}(\xi)$ 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 $\xi$. It is not possible to arbitrarily concentrate both a function and its Fourier transform." (source: Wikipedia, "Uncertainty principle" in https://en.wikipedia.org/wiki/Fouriertransform)

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:

$$\begin{aligned} F(t) &= \int_{0}^{\beta} \hat{G}(x, y, \omega) i\omega \hat{p}(\omega) \exp(i\omega t) d\omega \\ &= \int_{0}^{\epsilon} \hat{G}(x, y, \omega) i\omega \hat{p}(\omega) \exp(i\omega t) d\omega + \int_{\epsilon}^{\beta} \hat{G}(x, y, \omega) i\omega \hat{p}(\omega) \exp(i\omega t) d\omega \end{aligned} \qquad (C.10)$$

where $\beta$ is the equivalent of $t$ in the frequency domain; and $\epsilon$ is a small frequency which indicates the frequencies over which the ratio $\frac{\hat{J}^k(\omega)}{\hat{J}^{k+1}(\omega)}$ is not constant.

Over $[\epsilon; \beta]$, the ratio $\frac{\hat{J}^k(\omega)}{\hat{J}^{k+1}(\omega)}$ is constant. Thus, it makes sense to write:

$$F(t) = \int_{0}^{\epsilon} \hat{G}(x, y, \omega) i\omega \hat{p}(\omega) \exp(i\omega t) d\omega + \bar{G}\left(x, y, \frac{J^k(t)}{j^{k+1}(t)}\right) \int_{\epsilon}^{\beta} i\omega \hat{p}(\omega) \exp(i\omega t) d\omega. \qquad (C.11)$$
It comes that the error in the calculation of $F(x, y, t)$ lies in the term $\int_0^\epsilon \hat{J}^1(\omega) \hat{G}(x, y, \omega) \exp(i\omega t) dw$.
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 Wallace et al. (2020)).

**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 $\hat{G}(x, y, \omega)$ 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 $\int_0^\epsilon \hat{J}^1(\omega) \hat{G}(x, y, \omega) \exp(i\omega t) dw$ 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 $\int_0^\epsilon \hat{J}^1(\omega) \hat{G}(x, y, \omega) \exp(i\omega t) dw$.
This term involves very low frequencies where the viscoelastic behaviour is not physically very significant. Thus, ignoring it implies an error which is not expected to be high.

## References

Assogba, O.C., Tan, Y., Zhou, X., Zhang, C., Anato, J.N., 2020. Numerical investigation of the mechanical response of semi-rigid base asphalt pavement under traffic load and nonlinear temperature gradient effect. Construction and Building Materials 235, 117406.

Bentall, R., Johnson, K., 1967. Slip in the rolling contact of two dissimilar elastic rollers. International Journal of Mechanical Sciences 9, 389–404.

Bufler, H., 1971. Theory of elasticity of a multilayered medium. Journal of Elasticity 1, 125–143.
Carter, F.W., 1926. On the action of a locomotive driving wheel. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 112, 151–157.

Chaise, T., Nelias, D., 2011. Contact Pressure and Residual Strain in 3D Elasto-Plastic Rolling Contact for a Circular or Elliptical Point Contact. Journal of Tribology 133, 041402.

Chen, J., Bull, S., 2009. Finite element analysis of contact induced adhesion failure in multilayer coatings with weak interfaces. Thin Solid Films 517, 3704–3711.

Chen, W., 1971. Computation of stresses and displacements in a layered elastic medium. International Journal of Engineering Science 9, 775–800.

Djabella, H., Arnell, R., 1994. Finite element analysis of elastic stresses in multilayered systems. Thin Solid Films 245, 27–33.

Farrell, W.E., 1972. Deformation of the Earth by surface loads. Reviews of Geophysics 10, 761.

Gallego, L., Nelias, D., Deyber, S., 2010. A fast and efficient contact algorithm for fretting problems applied to fretting modes I, II and III. Wear 268, 208–222.

Gorishnyy, T.Z., Olson, L.G., Oden, M., Aouadi, S.M., Rohde, S.L., 2003. Optimization of wear-resistant coating architectures using finite element analysis. Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films 21, 332–339.

Goryacheva, I., Goryachev, A., Sadegi, F., 1995. Contact of elastic bodies with thin visco-elastic coatings under conditions of rolling or sliding friction. Journal of Applied Mathematics and Mechanics 59, 607–614.

Goryacheva, I., Miftakhova, A., 2019. Modelling of the viscoelastic layer effect in rolling contact. Wear 430-431, 256–262.

Guler, M.A., Adibnazari, S., Alinia, Y., 2012. Tractive rolling contact mechanics of graded coatings. International Journal of Solids and Structures 49, 929–945.

Johnson, K.L., 1985. Contact mechanics. Cambridge University Press, Cambridge [Cambridgeshire] ; New York.

Kalker, J., 1967. On the rolling contact of two elastic bodies in the presence of dry friction. Ph.D. thesis. Delft University.

Kalker, J.J., 1990. Three-Dimensional Elastic Bodies in Rolling Contact. volume 2 of *Solid Mechanics and Its Applications*. Springer Netherlands, Dordrecht.

Kalker, J.J., 1991. Viscoelastic multilayered cylinders rolling with dry friction. Journal of Applied Mechanics 58, 666–679.

Komvopoulos, K., 1988. Finite element analysis of a layered elastic solid in normal contact with a rigid surface. Journal of Tribology 110, 477–485.

Kot, M., 2012. Contact mechanics of coating-substrate systems: Monolayer and multilayer coatings. Archives of Civil and Mechanical Engineering 12, 464–470.

Koumi, K.E., Chaise, T., Nelias, D., 2015. Rolling contact of a rigid sphere/sliding of a spherical indenter upon a
viscoelastic half-space containing an ellipsoidal inhomogeneity. Journal of the Mechanics and Physics of Solids 80, 1–25.

Koumi, K.E., Nelias, D., Chaise, T., Duval, A., 2014. Modeling of the contact between a rigid indenter and a heterogeneous viscoelastic material. Mechanics of Materials 77, 28–42.

Kuo, J.T., 1969. Static response of a multilayered medium under inclined surface loads. Journal of Geophysical Research 74, 3195–3207.

Kusche, S., 2017. Frictional force between a rotationally symmetric indenter and a viscoelastic half-space: Elastomer friction. ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik 97, 226–239.

Lee, E.H., Radok, J.R.M., 1960. The contact problem for viscoelastic bodies. Journal of Applied Mechanics 27, 438–444.

Liu, S., Wang, Q., 2002. Studying contact stress fields caused by surface tractions with a discrete convolution and fast fourier transform algorithm. Journal of Tribology 124, 36–45.

Liu, S., Wang, Q., Liu, G., 2000. A versatile method of discrete convolution and FFT (DC-FFT) for contact analyses. Wear 243, 101–111.

Malvern, L.E., 1969. Introduction to the mechanics of a continuous medium. Prentice-Hall series in engineering of the physical sciences, Prentice-Hall, Englewood Cliffs, N.J.

Manyo, E.Y., Reynaud, P., Picoux, B., Tautou, R., Allou, F., Petit, C., Nelias, D., 2020. Tire-pavement tractive rolling contact under turning conditions: towards pavement top-down cracking. International Journal of Pavement Engineering, 1–10.

Menga, N., Afferrante, L., Carbone, G., 2016. Effect of thickness and boundary conditions on the behavior of viscoelastic layers in sliding contact with wavy profiles. Journal of the Mechanics and Physics of Solids 95, 517–529.

Menga, N., Carbone, G., Dini, D., 2021. Exploring the effect of geometric coupling on friction and energy dissipation in rough contacts of elastic and viscoelastic coatings. Journal of the Mechanics and Physics of Solids 148, 104273.

Menga, N., Putignano, C., Carbone, G., Demelio, G.P., 2014. The sliding contact of a rigid wavy surface with a viscoelastic half-space. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 470, 20140392.

Nguyen, V.T., Hwu, C., 2020. Time-stepping method for frictional contact of anisotropic viscoelastic solids. International Journal of Mechanical Sciences 184, 105836.

Nogi, T., Kato, T., 1997. Influence of a hard surface layer on the limit of elastic contact-part I: Analysis using a real surface model. Journal of Tribology 119, 493–500.

Nowell, D., Hills, D., 1988. Tractive rolling of dissimilar elastic cylinders. International Journal of Mechanical Sciences 30, 427–439.

O'Sullivan, T.C., King, R.B., 1988. Sliding contact stress field due to a spherical indenter on a layered elastic half-
space. Journal of Tribology 110, 235-240.

Peltier, W.R., 1974. The impulse response of a Maxwell Earth. Reviews of Geophysics 12, 649.

Plumet, S., Dubourg, M.C., 1998. A 3-D model for a multilayered body loaded normally and tangentially against a rigid body: Application to specific coatings. Journal of Tribology 120, 668-676.

Polonsky, I.A., Keer, L.M., 2000. A fast and accurate method for numerical analysis of elastic layered contacts. Journal of Tribology 122, 30-35.

Putignano, C., Carbone, G., 2014. A review of boundary elements methodologies for elastic and viscoelastic rough contact mechanics. Physical Mesomechanics 17, 321-333.

Putignano, C., Carbone, G., 2018. Viscoelastic reciprocating contacts in presence of finite rough interfaces: A numerical investigation. Journal of the Mechanics and Physics of Solids 114, 185-193.

Putignano, C., Carbone, G., Dini, D., 2015. Mechanics of rough contacts in elastic and viscoelastic thin layers. International Journal of Solids and Structures 69-70, 507-517.

Putignano, C., Carbone, G., Dini, D., 2016. Theory of reciprocating contact for viscoelastic solids. Physical Review E 93, 043003.

Scaraggi, M., Comingio, D., 2017. Rough contact mechanics for viscoelastic graded materials: The role of small-scale wavelengths on rubber friction. International Journal of Solids and Structures 125, 276-296.

Stepanov, F.I., Torskaya, E.V., 2016. Study of stress state of viscoelastic half-space in sliding contact with smooth indenter. Journal of Friction and Wear 37, 101-106.

Stepanov, F.I., Torskaya, E.V., 2018. Modeling of Sliding of a Smooth Indenter over a Viscoelastic Layer Coupled with a Rigid Base. Mechanics of Solids 53, 60-67.

Torskaya, E.V., Stepanov, F.I., 2019. Effect of surface layers in sliding contact of viscoelastic solids 3-D model of material. Frontiers in Mechanical Engineering 5, 26.

Wallace, E.R., Chaise, T., Nelias, D., 2020. Three-dimensional rolling/sliding contact on a viscoelastic layered half-space. Journal of the Mechanics and Physics of Solids 143, 104067.

Wang, Z., Jin, X., Keer, L.M., Wang, Q., 2012. A numerical approach for analyzing three-dimensional steady-state rolling contact including creep using a fast semi-analytical method. Tribology Transactions 55, 446-457.

Wayne Chen, W., Jane Wang, Q., Huan, Z., Luo, X., 2011. Semi-analytical viscoelastic contact modeling of polymer-based materials. Journal of Tribology 133, 041404.

Xu, J., Wang, K., Liang, X., Gao, Y., Liu, Z., Chen, R., Wang, P., Xu, F., Wei, K., 2020. Influence of viscoelastic mechanical properties of rail pads on wheel and corrugated rail rolling contact at high speeds. Tribology International 151, 106523.

Yu, C., Wang, Z., Wang, Q.J., 2014. Analytical frequency response functions for contact of multilayered materials. Mechanics of Materials 76, 102-120.

Zhang, X., Wang, Q.J., He, T., 2020. Transient and steady-state viscoelastic contact responses of layer-substrate systems with interfacial imperfections. Journal of the Mechanics and Physics of Solids 145, 104170.