*3.1. Bulk Stresses*

Figure 5 shows the results of FEM simulations for the boundary value problem shown in Figure 4. They refer to the single arm model, but, from a qualitative point of view, the considerations that are going to be drawn below for the bulk stresses also apply to the other two models herein considered.

Figure 5a,b display the stresses developing in the bulk at the end of the normal loading stage, and they display three distinct areas with high stresses where the harmonic profile comes into contact. Because of the presence of friction and, since coupling effect are fully included, an anti-symmetric distribution of *τxz* arises, even in the pure normal loading stage, see Figure 5b. The following two figures represent the same quantities at a subsequent load stage, where the normal imposed displacement has reached its maximum, and the indenter slides at constant velocity. Figure 5c,d show the stresses during the next stage of sliding, corresponding to a lateral shift of the harmonic profile of about half of its wavelength. The advantage of the finite element method is evident from the possibility to consider any finite-size problem geometry and boundary conditions, like, in this case, the output automatically including not only contact tractions, but also bulk stresses.

**Figure 5.** Model predictions: bulk stresses during the normal approach, (**a**,**b**), and during full sliding, (**c**,**d**), all scaled by a reference elastic modulus *E<sup>f</sup>* ,0 = 8.147 × 10 <sup>2</sup> Pa.

## *3.2. Interface Tractions*

Figure 6 highlights the evolution of contact tractions in time for the single arm model and selected stages of the contact simulation. The curves in Figure 6a correspond to the purely normal loading sequence, where normal contact tractions progressively increase along with the value of the applied normal displacement, which linearly rises from zero up to the final value of 2*g*0. Black curves denote the symmetric distribution of normal contact tractions *pz*(*x*) divided by *E<sup>f</sup>* ,0 , while red curves represent the anti-symmetric distribution of tangential contact tractions *qx*(*x*), scaled by *f E<sup>f</sup>* ,0 . Points along the interface, where k*qx*(*x*)k/(*f E<sup>f</sup>* ,0 ) equals *pz*(*x*)/*E<sup>f</sup>* ,0 , are in a state of slip, while, when the inequality k*qx*(*x*)k/(*f E<sup>f</sup>* ,0 ) < *pz*(*x*)/*E<sup>f</sup>* ,0 holds, then there is a state of stick.

Figure 6b refers to the next stage of the contact problem when, keeping the normal displacement constant, a far-field displacement linearly increasing with time is applied in the tangential direction. While, for the given rheological model, the results that are shown in Figure 6a are evaluated in a condition of zero tangential velocity, Figure 6b–d are referred to *<sup>v</sup>* <sup>=</sup> 2.154 <sup>×</sup> <sup>10</sup>−<sup>2</sup> m/s. This specific value has been chosen amid the other entries of Table 2, because it is in the middle of the range, determining the highest viscoelastic effects, and it is also low enough for analysing the transition from *stick/slip* to *full sliding*, Figure 6b. Here, tangential traction distributions change their shape from the classical anti-symmetric form towards a state of increasing slip, which terminates in the full slip condition. The transition from *stick-slip* to *full slip* is strongly affected by the velocity of the horizontal displacement: the faster the slip, the more abrupt such a transition.

Figure 6c refers to the situation of sliding after full slip (*gross sliding*) and, in particular, it shows the evolution over time of the normal contact tractions. We see a transition from the symmetric contact traction distribution along the whole interface at the onset of full slip, as shown in black, towards other distributions in different scales of grey shifted along the interface to the right, as long as the tangential displacement increases. A certain degree of relaxation is observed after the onset of full slip. As the sliding proceeds in time, virgin material is perturbed, and a recovery in stiffness takes place.

(**a**) *p<sup>z</sup>* and *q<sup>x</sup>* during the normal loading stage. (**b**) Tractions in the partial-slip regime.

an already stressed portion of the interface.

Finally, Figure 6d captures the first overlapping of a new contact zone with a previously loaded portion of the interface. Here, the role of the relaxation time is important, since viscoelastic effects do alter the solution that corresponds to a linear elastic material that has no memory effects.

The resultant tangential force *Qx*, integral of tangential contact tractions along the interface, is plotted vs. time in Figure 7a–c for the three viscoelastic models investigated herein. In each subfigure, different curves correspond to different far-field horizontal displacement velocities. Darker curves correspond to slower velocities.

**Figure 7.** Time evolution of the resultant tangential force *Q<sup>x</sup>* for different rheological models.

In all of the cases, for *t*/*t*<sup>0</sup> ≤ 1, tangential tractions are vanishing, since, in that stage, the imposed displacement is only acting in the normal direction. Therefore, tangential contact tractions are due to frictional coupling effects and their sum over the whole contact zones is vanishing by definition, since they correspond to self-equilibrated distributions. For *t*/*t*<sup>0</sup> > 1, the indenter starts sliding and we assist to a transition from *stick-slip* to *full slip* with an oscillatory behaviour when the contact profile enters in contact with unrelaxed material portions. When the velocity is low, no rate effects are evident, and the mechanical response is smooth. On the other hand, by increasing the applied velocity, the importance of viscoelasticity increases and oscillating responses do appear.

*Lubricants* **2020**, *8*, 107

The integral of tangential tractions related to two linear elastic models that are characterised by short and long term modulus are also plotted in Figure 7; for comparison, see black dash-dotted lines. The elastic moduli are evaluated as:

$$E^{\text{el},\infty} = \lim\_{t \to 0} E(t) = E^{\infty} \qquad \qquad \qquad E^{\text{el},0} = \lim\_{t \to \infty} E(t) = E^{\infty} (1 - \sum\_{i=1}^{n} \mu\_i) \tag{13}$$

The curves *Q* el,∞ *<sup>x</sup>* and *Q* el,0 *<sup>x</sup>* are evaluated under the assumption of linear elasticity, neglecting the dynamic effects. For this reason, they lead to constant values as soon as the horizontal far-field displacement is applied, without any oscillation. The only factor that plays a role is the velocity, which governs the transition from *stick/slip* to full sliding. In the figures, only the curves that correspond to the highest value of *v* are plotted. In all three models, the instantaneous (higher) and long term (lower) curves are extreme bounds to the values that are related to viscoelastic simulations, with a gap increasing from the single arm to the three arms model, consistent with their respective stiffness.

The steady-state solution strongly depends on the rheological properties of the material, as shown in Figure 7d. In general, for the present case study, the higher the number of arms, the higher the total tangential force. In all cases, the highest velocity determines the highest value of the steady state *Q*<sup>0</sup> *x* . This is in accordance with the fact that, in a condition of *gross slip*, *Q<sup>x</sup>* = *f Nz*, and for high velocities, the material is excited in its high frequency region, thus resulting in a vertical response that is governed by the higher *glassy* Young's modulus. The increased stiffness leads to higher *N*<sup>0</sup> *<sup>z</sup>* and, in turn, higher *Q*0 *<sup>x</sup>* values.

## **4. Conclusions**

In this study, a novel finite element procedure has been proposed, which allows for investigating transient and steady state sliding of a rigid indenter over a viscoelastic continuum. In particular, the representative problem of an indenter with harmonic profile sliding over a viscoelastic layer of finite depth has been analysed, employing different sliding velocities together with three different rheological models, which are characterised by Prony series with one, two, and three arms, respectively. A regularised version of the classic Coulomb friction law has been employed for the evaluation of the interface tangential tractions.

Numerical results pinpoint a strong dependence of the mechanical response in terms of steady-state forces *N<sup>z</sup>* and *Q<sup>x</sup>* on both the velocity and rheological model employed, obtaining increasing forces for higher velocities and more relaxation terms that are involved in the rheological approximation.

It is worth mentioning that the proposed methodology appears to be suitable for the investigation of a class of problems for which a solution could be difficult to be found while using other techniques. The proposed approach is capable of overcoming the limitations of other solution schemes thanks to the capability of FEM of solving linear and nonlinear boundary value problems with arbitrary material constitutive laws and geometries. Moreover, the use of the recently developed interface finite element [6–8] has further advantages. First of all, the possibility of taking into account arbitrary shapes for the indenting profile as analytical functions that are embedded into the interface element. The ability of simulating partial slip scenarios involving finite sliding of the indenter should also be mentioned.

Moreover, as a key advantage when compared to other models that are available in the literature that neglect the effect of Coulomb friction, focusing on viscoelastic dissipation only, here viscoelastic effects and frictional effects can be simultaneously investigated, since they are inherently coupled in the formulation. Neglecting interface tangential tractions, together with their related coupling affecting the distribution of normal tractions, could be reasonable when incompressibility conditions are approached. On the other hand, several evidences can be found that, as the Young's modulus of a viscoelastic material changes with time, so does the Poisson's ratio. Because the latter quantity governs the coupling between normal and tangential tractions, a fully coupled model is worth study for fine precision engineering applications. As a final remark, the proposed interface finite element has the further advantage of being easily extended for taking thermal effects into account. These could be relevant not only for the analysis of temperature transfer across the interface, but also to simulate frictional heat generation, thus leading to a thermodynamically accurate model that is capable of investigating a wide class of realistic viscoelastic dissipative phenomena.

**Author Contributions:** Conceptualization, J.B. and M.P.; methodology, J.B. and M.P.; software, J.B.; validation, J.B.; formal analysis, J.B. and M.P.; investigation, J.B. and M.P.; resources, M.P.; data curation, J.B.; writing—original draft preparation, J.B.; writing—review and editing, M.P.; visualization, J.B.; supervision, M.P.; project administration, M.P.; funding acquisition, M.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received funding from the Italian Ministry of University and Research through the Research Project of National Interest (PRIN 2017) "XFAST-SIMS: Extra fast and accurate simulation of complex structural systems" (grant no. D68D19001260001).

**Conflicts of Interest:** The authors declare no conflict of interest.
