Viscoelastic Damping in alternate reciprocating contacts

Reciprocating motion between viscoelastic solids occurs in a number of systems and, in particular, in all the dampers which exploits, as a physical principle, the viscoelastic dissipation. So far, any attempt to predict the behavour of this field of dampers relies on approximate methodologies and, often, on a steady-state approach, with a consequent poor understanding of the phenomenon. Here, we develop a methodology capable of simulating the actual mechanics of the problem and, in particular, we shed light on how the presence of not fully relaxed viscoelastic regions, during the punch motion, determine the viscoelastic dissipation. The latter is shown to be dependent ultimately on two dimensionless parameters, i.e. the maximum speed in the cycle and the frequency. Finally, the importance of considering a rough interface is enlightened.

in the Boundary Element method. Section III include the results in terms of deformation -the dissipated power is indeed proportional to deformed volume-and in terms of the tangential force in a range of different conditions of frequency and maximum speed. Furthermore, the role of the roughness is discussed. Final remarks close the paper and comment on the significance of the reciprocating regime to assess the dissipated power.

Formulation
By recalling the translational invariance and the elastic-viscoelastic correspondence principle 22 , the general linear-viscoelastic contact problem between a rigid indenter and a viscoelastic slab can be formulated by means of the following integral equations relating interfacial stresses and strains: where x is the in-plane position vector, t is the time, u(x, t) and σ(x, t) are respectively the normal surface displacement of the viscoelastic solid and the normal interfacial stress, x ( )  and t ( )  are the Green's function and the creep function, and the symbol '·' refers to the Lagrangian time derivative. The latter quantity is the creep function which satisfies causality, i.e.  < = t ( 0) 0, and, for a generic viscoelastic material 22 , can be written as: where t ( )  is the Heaviside step function, the real quantity E 0 is the low-frequency elastic modulus, τ ( )  is a strictly positive function usually defined as the creep (or retardation) spectrum 22,23 , and τ is the relaxation time, continuously distributed on the real axis. Now, solving directly Eq. (1) requires to discretize both the time and the space domains and this may be unfeasible in a number of cases, and, remarkably, when roughness has to be accounted for. Indeed, given a rough spectrum, covering typically several orders of magnitude, the number of boundary elements required would be too large for the computational resources currently available.
However, since we are going to study a reciprocating motion, given the periodicity of the system, we can try to simplify Eq. (1). To this end, let us start assuming that the interfacial normal stress distribution σ(x, t) is equal to : this means that the shape of interfacial distribution is time-indipendent and moves on the viscoelastic half-space with a sinusoidal law of amplitude ξ 0 and angular frequency ω. Now, because of the linearity and the translational invariance characterizing the problem, we can carry out the replacement , and therefore Eq. (1) can be rewritten as: 2 A significant advantage is embedded in Eq. (3): there exists only a parametric dependence on the time and, therefore, only the space domain has to be discretized to solve completely the problem. Now, G(x, t) has to be determined. To this end, as shown in refs 12 and 21, it is possible to recall the general relation between stress and displacement distributions: Furthermore, given the constant shape assumption previously stated, we can replace σ( in Eq. (3) and, then, obtain: can be re-written as: Replacing Eq. (10) into Eq. (7), we obtain i t k k 0  and, taking the inverse Fourier transform, we can rewrite Eq. (6) as: In order to efficiently estimate G(x, t), let us focus preliminarily on the term A k (x). Indeed, it is possible to apply the Chebishev-Gauss quadrature rule to the integral term ∫ at n nodes so that A k (x) can be numerically estimated as: . Now, we could directly calculate G(x, t) by means of a numerical estimation of the series in Eq. (15). However, in this paper, we propose a much more efficient approach. Indeed, what Eq. (15) says is that G(x, t) can be expressed as a Fourier series whose coefficient are U k . Consequently, if we generate an array = ... ...
, i.e. the Green's function for 2N + 1 values of the time t, simply by using the Inverse Fast Fourier Transform. In comparison with the estimation in the time domain of the series in Eq. (15), such an approach makes the search for the converged value of G(x, t) easier and speeds the entire calculation of orders of magnitude also when the frequency ω is small and computing in the time domain tends to be very slow. Incidentally, it should be observed that the final expression for G(x, t) relies on the factorization, performed in Eq. (1), of the integral equation kernel in two terms, that are  x ( ) and  t ( ). This is allowed due to the assumption on the rheology of the contacting solids, which are considered homogeneous linear viscoelastic. In the case of materials exhibiting a more general time-dependent behaviour, given the periodicity of the system, G(x, t) could be again expressed as a Fourier series, whose coefficients would depend on the particular conditions marking the bodies in contact. This is, however, out of the scope of the paper.
As suggested for the elastic case in refs 24 and 25. and for viscoelastic systems in ref. 12, once the Green's function G(x, t) is known, we can mesh the contact domain in M square cells and, then, discretize Eq. (3) by assuming that in each boundary element the normal stress σ is constant and equal to σ j . The displacement is, then, measured at the centre x i of the each i-th square and is equal to is reduced to the following linear system: . Indeed, the contact problem can be solved by using the iterative technique developed in ref. 24. for elastic contacts, thus providing contact areas, stresses and strains. Interestingly, the method does not require any discretization of the time domain as t is treated as a parameter. Furthermore, once the solution is known in terms of stresses and strains, following the approach commonly used in viscoelastic contact mechanics, stated in ref. 12, it is straightforward to calculate the viscoelastic interfacial tangential force as: Finally, we observe that G(x, t) has been calculated under the assumption that the shape of the stress field at the interface, which is generally equal to σ , where a 0 is the characteristic dimension of the contact region.

Results and Discussion
Let us start focusing briefly on the properties of the Green's function G(x, t), whose formulation has been shown to be parametrically dependant on the time. Indeed, given the scope of this paper, i.e. damping in viscoelastic reciprocating contacts, it is important to analyze the displacement field since the amount of dissipated energy is proportional to the region within the material, which undergoes time-dependent deformations 12,20 . The Green's function, being just the displacement produced by a concentrated force, represents, then, the simplest displacement distribution that can be obtained. Now, to carry out the analysis, we employ a simple one relaxation time material with the high frequency modulus E ∞ equal to  Figure 1 shows the evolution of the Green's function for different values of the time ω π π ∈ − t [ /2, /2]. An arrow refers, in each case, to the current position of the force. At ω π =− t /2 the force has just reached the left dead-point and starts moving from left to right: the speed grows and, consequently, the material stiffens and the displacement decreases. Later on, for time larger than ω = t 0, the speed starts to decrease until the force stops for ω π = t /2. We observe that the displacement distribution shows three peaks: one corresponds to the current position of the force; the other two are found at the dead-points and are due to a typically viscoelastic delay in the deformation recover. Ultimately, the material is not able to fully relax and, consequently, to recover the viscoelastic deformation during a time interval comparable to the oscillation period.
The occurrence of multi-peaked displacement distributions is intrinsically related to viscoelasticity and characterizes also more complicated loading conditions, as for example when a rigid punch moves over regions that are still relaxing. To show this, let us study the contact of a rigid sphere of radius R which is in reciprocating sliding over a viscoelastic material with the following motion law Coherently with what we observed before when looking at the Green's function, upon reversal of the sliding direction, the punch moves over a viscoelastic region that has not yet had the time to recover the deformation after the previous contact of the rigid body. Furthermore, as the speed increases, the material stiffening is not negligible anymore and a marked decrease of the penetration is noticed. Ultimately, three different deformation peaks are observed in the path covered between the two dead-points: one corresponds to the current position of the punch, and the other two, being close to the dead-points, occur because the material does not fully recover the deformation in the oscillation period. Interestingly,we notice that such an interplay between the displacement fields provoke stress distributions, which show, upon the motion inversion, the jagged trend observed in Fig. 3. Indeed, starting from ωt = −π/2, although the sliding speed goes to zero, the interfacial normal stress distribution has an asymmetric shape, marked by a peak on the left of the contact patch: this is due to a tipically viscoelastic time-delay which prevents the material to relax immediately when the sliding speed vanishes. Clearly, such an effect cannot be embedded with a steady-state approch, which, for a speed equal to zero, predicts the static Hertzian solution. Now, as the sphere starts moving to the right, the peak shown by the reciprocating model at ωt = −π/2 does not disappear immediately but has to gradually decrease. At the same time, the punch is moving towards the right and, as shown also in steady-state conditions, a peak at the leading edge has to appear in the pressure distribution. Furthermore, at the center of such distribution, in correspondence of the maximum of the displacement field in the contact area, the pressure must still preserve an Hertzian-like trend. As a result, we have the multi-peaked pressure distributions observed in Fig. 3: the deviation from the steady-state predictions is evident and, interestingly, is significant also at ωt = −0.28π, i.e. when the punch is far enough from the inversion point and the pressure has evolved towards a smoother distribution. Indeed, this a further evidence of the necessity of considering the viscoelastic reciprocating properties. Now, in order to relate such a mechanical behavior with viscoelastic tangential forces and with the energy dissipation, let us identify the parameters that govern the phenomenon. Naturally, the ratio ∞ E E / 0 plays an important role since it quantifies the degree of viscoelasticity: if = ∞ E E / 1 0 , it is straightforward to show that the material is perfectly elastic. Indeed, due to the sum rule 12 , we can write: However, more specifically related to the reciprocating motion and, therefore, to parameters like the stroke ξ 0 , the frequency ω and the contact area a 0 , we find two different dimensionless groups 21 , which govern the problem. The first one is Π = τ/t 0 , where ωξ = t a / 0 0 0 and a 0 is the Hertzian contact radius, and compares the relaxation time τ with the time t 0 needed by the sphere to cover a distance a 0 . Π can be seen, indeed, as the maximum dimensionless speed reached by the punch during its motion since it is equal to ωξ τ τ Π = = a v a / / 0 0 m ax 0 . When Π ≪ 1, we are then in soft elastic conditions; for very large values of Π, then, we are in the glassy elastic regime in the most part of the track, but close to the inversion points, where the speed vanishes, some viscoelastic effects may be present. The second dimensionless parameter, governing the phenomenon, is the dimensionless frequency ωτ πτ Ξ = = T 2 / . Given a regime where viscoelasticity plays an important role, i.e. Π ≈ 1, if Ξ < 1, the reciprocating motion occurs on time-scales longer than the material relaxation time τ and the system behaves like in the steady-state case 12 . If Ξ ≈ 1, the system will show the significant effects specific of the reciprocating motions, as previously enlightened in Fig. 2. Finally, for very large values of Ξ ≫ 1, the punch carries out very fast reciprocating oscillations with an almost negligible stroke and again the high frequency elastic behavior will be recovered. Incidentally, we observe that, because of the constant shape assumption, which requires ξ Ξ Π =  a / / 1 0 0 , our methodology should not be utilized to predict the behavior of the system when ξ when Ξ/Π ≫ 1 the behavior of the system is elastic, i.e. memory effects disappears, and the condition ξ Ξ Π =  a / / 1 0 0 becomes needless, as the purely elastic regime is intrinsically embedded in the proposed methodology. In fact, the Green's function G(x, t) in Eq. (3), for large Ξ/Π asymptotically converges to a time-indipendent purely elastic Green's function G(x) provided by the Boussinesq solution 26 .
Let us, now, study the influence of the parameters Ξ and Π on the tangential forces during the cycle. In Fig. 4, after fixing the maximum speed Π equal to Π = 23.5, the dimensionless tangential force F t /F n is plotted as a function of the dimensionless abscissa ξ x/ 0 , which identifies the position of the sphere center along the path, for different values of the frequency Ξ. In particular, blue lines refer to the simulations obtained by means of the reciprocating model, whereas values on the dotted red lines are computed by employing a steady-state model 12 given, in each instant t of the cycle, a rolling velocity . Note that the steady-state model neglects any reciprocating effect as the parameter Ξ does not play any role in this case.
As expected, for low values of Ξ, the reciprocating solution approaches the steady-state regime. A different behavior occurs, however, when Ξ is increased: upon the motion inversion, due to the interaction of the punch with a still deformed area, the actual reciprocating behavior strongly deviates from what is predicted by the steady-state model. In detail, because of the reciprocating interplay, when the motion is inverted, the tangential force F t /F n has the same direction as the sliding speed and the entire dissipation cycle appears twisted. Furthermore, as Ξ and Ξ/Π grow, the cycle area tends to reduce in agreement with what expected theoretically, i.e. that the system approaches an elastic behaviour. Now, after recalling that, in each point of the cycle, the punch has a speed v(t), let us introduce the dimensionless dissipated work per cycle δ: In Fig. 5, fixed Π, δ decreases with Ξ and goes to zero for Ξ ≫ 1. Therefore, employing an approach based on steady-state assumptions, as the one followed usually when designing dampers (see e.g. ref. 27.), may lead to a massive overestimation of the dissipated power and, then, to a poor design.
At this point, in order to study the influence of the other parameter, i.e. the maximum dimensionless speed in the cycle Π, it may be interesting to fix the frequency Ξ and sweep a range of values for Π. Let us preliminarly recall that, in steady-state conditions, where friction force and friction coefficient are functions only of the   . Now, in a reciprocating cycle, depending on the ratio ζ Π/ , we will find different shapes for the hysteretic curve. Indeed, as shown in Fig. 6, if ζ Π < , the viscoelastic tangential force F t /F n continuously increases when the speed grows in the cycle. The latter is, then, convex. In fact, since the instantaneously velocity v(t) is always smaller than the ζ τ a / 0 , the viscoelastic material is excited in a frequency range where the loss moduls of the material (i.e. ω Im E [ ( )]) always grows with the frequency 22,23 : this makes, indeed, friction an increasing function of the rolling velocity. On the contrary, when ζ Π > , F t /F n reaches again its maximum for ζ τ = v t a ( ) / 0 , but there will exist a time range where v(t) exceeds ζ τ a / 0 , thus leading the viscoelastic material to be exicited on the decreasing branch of the loss-modulus vs. frequency curve. As a matter of fact, the hysteretic curve tends to look spiked close to the dead-points. Such a trend is found both when employing the reciprocating method and in steady-state conditions, but significant quantitave differences occur between the two approaches. As previously discussed, this is due to the interplay between not fully relaxed regions of the viscoelastic half-space.
When looking at the dissipated work δ, from a quantitative point of view, we notice that δ decreases with Π and goes to zero once the system approaches the soft elastic regime. On the other side, for very large values of Π, the system is in the glassy elastic regime during the most part of the cycle, apart from very narrow regions close to each dead-point, where the punch arrests and inverts its motion. In Fig. 7, we can focus on δ and observe a bell-shaped curve: as expected, the real dissipation computed by the viscoelastic reciprocating model is smaller in comparison with the predictions carried out by means of the steady-state approach. Such a scenario confirms the necessity of the methodology introduced in this paper to account for the actual interactions in the reciprocating dynamics. This is even more true when the real rough spectra of the contacting surfaces is accounted for. The presence of a rough interface covering many orders of magnitude has, undoubtedly, a strong importance: as already shown in the steady-state case 17 , the viscoelastic dissipation grows by increasing the number of rough scales. In detail, given a self-affine fractal surface numerically generated, for example, by means of the spectral method described in ref. 25, viscoelastic friction grows significantly with the short-length cut-off wavevector q 1 = 2πN/L 0 , where L 0 is the roll-off Figure 6. The dimensionless tangential force F t /F n as a function of the dimensionless abscissa for several values of Π and Ξ = 0.5. Blue lines refer to results obtained by means of the viscoelastic reciprocating approach, whereas red dotted lines are obtained employing a classic steady-state model. length and N the number of scales. Such a trend is confirmed for the reciprocating conditions in Fig. 8, where the tangential force F t /F n is shown both for a smooth and for a rough sphere. This rough interface has been obtained spreading on the smooth sphere a fractal self-affine surface whose spectral components are in the range q 0 < q < q 1 , where q 0 = 2π/R and q 1 = Nq 0 . In particular, results in Fig. 8 are obtained with N = 32 and show, in addition to the increased dissipation, a different shape due to the presence of the roughness.

Conclusions
In this paper, we present a numerical methodology capable of providing accurate predictions in the case of reciprocating motion between viscoelastic solids. Implications for the damping of rubber-like components are enormous since reliable numerical predictions are finally possible. In particular, in the simple case of the reciprocating contact of a sphere over a viscoelastic half-space, the phenomenon is shown dependent just on two parameters, that are the maximum dimensionless speed Π reached during the motion and the dimensionless frequency Ξ. To this extent, it is crucial to notice that, if the analysis of the viscoelastic dissipation was conducted as usual done in literature with a steady-state approach (see e.g. ref. 27.), such a dependence on Ξ would be completely neglected. This issue is particularly critical for Ξ/Π ≫ 1, when the viscoelastic damping has to go to zero, as correctly predicted by the reciprocating model.
The damping overestimation, intrinsically related to the steady-state approach, is, indeed, related to the poor description that such stationary models produce on the reciprocating physics. Indeed, upon the motion inversion,  the punch gets into contact with a region of viscoelastic half-space, which is still deformed and has to recover the deformation yet. This produces the multi-peaked displacement distribution shown in the paper and, more interestingly, gives origin to the strong dependence on the frequency Ξ. All this is ignored from the steady-state point of view, which ultimately ignores all the reciprocating loading history of the material. This is even more significant when considering the surface roughness, which, as in any phenomenon related to viscoelastic contact mechanics, plays a crucial role. The viscoelastic material is stressed over a frequency range, whose width is strictly related to the rough spectra 17 , and, consequently, the interplay between regions not fully relaxed cannot be neglected.