Critical thresholds for mode-coupling instability in viscoelastic sliding contacts

Mode-coupling instabilities are known to trigger self-excited vibrations in sliding contacts. Here, the conditions for mode-coupling (or “flutter”) instability in the contact between a spherical oscillator and a moving viscoelastic substrate are studied. The work extends the classical 2-Degrees-Of-Freedom conveyor belt model and accounts for viscoelastic dissipation in the substrate, adhesive friction at the interface and nonlinear normal contact stiffness as derived from numerical simulations based on a boundary element method capable of accounting for linear viscoelastic effects. The linear stability boundaries are analytically estimated in the limits of very low and very high substrate velocity, while in the intermediate range of velocity the eigenvalue problem is solved numerically. It is shown how the system stability depends on externally imposed parameters, such as the substrate velocity and the normal load applied, and on contact parameters such as the interfacial shear strength τ0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{0}$$\end{document} and the viscoelastic friction coefficient. In particular, for a given substrate velocity, there exist a critical shear strength τ0,crit\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{0,crit}$$\end{document} and normal load Fn,crit\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_{n,crit}$$\end{document}, which trigger mode-coupling instability: for shear stresses larger than τ0,crit\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{0,crit}$$\end{document} or normal load smaller than Fn,crit\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_{n,crit}$$\end{document}, self-excited vibrations have to be expected.


Introduction
Friction-induced vibrations (FIVs) are an incredibly widespread phenomenon, where friction plays an unexpected role: in fact, instead of acting as a damping source, attenuating the system dynamics, dissipative forces are the original cause that triggers this class of vibrations. The large interest in the analysis of FIVs is not only theoretical and aimed to assess the intricate correlation with friction [1][2][3], but it is also largely motivated by the variety and the importance of fields and applications, where FIVs occur: these range from automotive and railway systems [4][5][6][7][8][9][10][11] to the aerospace field [12][13][14] and include also areas far from industrial engineering, like bio-mechanics [15][16][17][18]. In all these different systems, FIVs may be related to unwanted acoustic emissions and to occurrence of additional stresses in the components affected by this phenomenon. Several mechanisms can originate FIVs. Among these: (i) the "sprag-slip effect", pioneeringly investigated by Spurr in Ref. [19], was related to "jamming" phenomena at the interface level, (ii) the "Stribeck's effect", which is due to a falling characteristic of the friction law [20][21][22][23][24][25] and (iii) the "mode-coupling" instability (also referred as flutter), which is the result of the coupling of two stable vibrational modes that originates one stable and one unstable mode [26]. The latter is due to friction-related non-symmetric entries in the stiffness matrix, which makes the system not conservative and allows energy to flow from the interface to the vibrational motion [4]. In the last two decades, a large research effort has been dedicated to characterize the origin of FIVs, developing a minimal, yet representative, model consisting in a 2-DOF (Degrees Of Freedom) oscillator, in contact with a moving frictional substrate [4] [11]. The basic model has been originally conceived by Hamabe et al. in 1999 [20] to investigate the insurgence of FIVs in drum brakes; later on, the scheme has been extensively employed for a number of different applications and, more generally, for improving our understanding of FIVs main features [22,23,[26][27][28][29][30][31][32][33]. It is important to notice that, although real structures may be affected by a much higher level of complexity, this 2-DOF dynamic system is a paradigmatic and illustrative scheme that allows us to explore the main features of FIVs at a fundamental level. By means of such a model, in Ref. [29] an intuitive explanation is provided for the insurgence of FIVs, by showing that out-of-phase variation of the frictional force and tangential displacement may lead to energy transfer from interface to the mechanical system. Furthermore, in Ref. [30], Hoffmann and Gaul investigated the effect of added linear damping to the insurgence of flutter instability, while in Ref. [31], an approximate estimation of the limit cycle oscillations has been proposed. These studies have provided fundamental contributions to the FIVs assessment, but, at the same time, most of them have been affected by a significant issue: the interface interactions between the oscillator and the moving substrate are treated by means of phenomenological laws, the simplest one being the classical Coulomb model, aimed at describing the contact interactions between the oscillator and the substrate [33]. However, this phenomenological approach leads to a strong simplification of the contact mechanics occurring at the interface: for example, contact stiffness is often assumed linear [34] and/or the frictional resistance is taken proportional to the contact normal force [22,[35][36][37][38].
Only very recently, some attempts have been made to account for a more realistic description of the interfacial phenomena, as in Ref. [11], where the problem of a spherical punch in contact with a viscoelastic sub-strate has been investigated. Nevertheless, Ref. [11] was specifically focused on enlightening the role that viscoelastic dissipation plays in triggering the instability, but the analysis did not consider any elastic crosscoupling in the mechanical system between the horizontal and vertical motion and neglected the Coulomb "adhesive" part of the friction coefficient. As a consequence, the possible mechanisms for instability were restricted just to the Stribeck's effect. However, as previously pointed out, and as schematically depicted in Fig. 1, in practical applications, cross-coupling may play a fundamental role. Among these applications, we recall the case of tip rubbing in turbomachinery, where self-excited vibrations are triggered by the contact between a rotating blade and the external (fixed) case, which is coated by a polymer viscoelastic layer. Self-sustained vibrations and internal resonances may originate unexpected large stresses, wear and damage [12][13][14]. Similar issues may arise when a rubbermade face seal is inserted between rotor and stator in a mechanical system [39]: because of the roughness at the microscales, there may occur friction-induced vibrations [40]. There exist similar conditions, again at a microscale level, in the windscreen wiper blades, where FIVs are correlated with tedious acoustic emissions [41], and, even in the bio-mechanical field, when fingers touch and slide over rough surfaces [15][16][17][18]. Finally, when dealing with dampers [42,43], where rigid punches are put into contact with viscoelastic layers, FIVs can occur and change the operational parameters according to which the system has been designed. All these systems are characterized by an elastic structure in contact with a moving soft substrate, marked by a time-dependent viscoelastic rheology. Consequently, we focus on the simple yet representative 2-DOF conveyor belt model (see Fig. 1). Such a scheme is a classical lumped model to study FIVs, but here we introduce two fundamentally new aspects, enabling us to obtain a better description of the interface interactions and, thus, of the conditions for instability. First of all, contact stiffness and viscoelastic dissipation are obtained from adhesiveless (steady-state) contact mechanics boundary element method numerical simulations [39,[44][45][46][47][48][49]. Secondly, in accordance with up-to-date contact mechanics experiments involving soft contacts (polymers and elastomers, [50][51][52][53]), the Coulomb-like fraction of the friction forces is assumed to be proportional, through a characteristic interfacial shear-strength τ 0 , to The paper is organized as follows. In Sect. 2, the problem is defined mathematically, by introducing the equations of motion and the form of the contact forces that will be used along the paper. In Sect. 3, a linear stability analysis is performed, and the stability boundaries are derived analytically in the limit of low and high sliding velocity. In Sect. 4, the results of the stability analysis are discussed and, finally, in the conclusions, final remarks are drawn.

Dynamic model
The mechanical system under investigation is constituted by a rigid spherical oscillator with radius R, mass m, horizontal and vertical stiffness k x and k y , offdiagonal spring k xy (inclined by an angle α, see Fig. 2), horizontal and vertical linear damping coefficients c x and c y . The oscillator is pressed by a normal force F n against a viscoelastic substrate moving at velocity v d in the horizontal direction; x is the parameter corresponding to the horizontal displacement of the center of mass of the spherical oscillator, whereas y is the vertical displacement. The latter is equal to 0 when the sphere first touch the substrate: thus, y corresponds to the sphere indentation into the moving viscoelastic substrate.
x + k x + k xy cos 2 α x + k xy sin α cos α y = F μ m ·· y+c y · y+ k y +k xy sin 2 α y+ k xy sin α cos α x=F n −F k (1) where a dot superposed means differentiation with respect to time. Regarding the external forces, F μ is the frictional force acting along the x-axis and is related to the hysteresis occurring due to deformation of the viscoelastic substrate (see Ref. [44] for more details); F k is the normal contact force in the vertical direction and is, thus, equal to the integral of the contact pressure over the contact area. F k is larger than 0 when the sphere indents the substrate, i.e., for y > 0. Consequently, both F μ and F k depend on the viscoelastic properties of the substrate and have been numerically obtained by means of boundary element method simulations.
The interaction between the spherical punch and the viscoelastic substrate determines the normal restoring force F k and the friction force F μ acting on the sphere. In general, as shown in detail in Ref. [11], these components depend on the kinematic parameters that determine the relative motion between the sphere and the substrate, hence for a viscoelastic material with a single relaxation time τ we have: Ultimately, this occurs as, in the mechanics between viscoelastic elements, the contact solution depends, at each time step, on the current deformation of the substrate and on the velocity at which the material is being deformed. Furthermore, in the contact mechanics analysis we have conducted via boundary element method, we have accounted for the material hysteretic dissipative behavior in the in-plane motion, i.e., along the x−direction, which plays a major role as determines the viscoelastic part of the friction coefficient, while the dissipation occurring in the vertical motion will be accounted through the linear damping coefficient c y . Now, when looking at Eq. (2), it is clear that viscoelastic materials respond differently when excited at different frequencies. Then, if the punch is oscillating in the horizontal direction with frequency ω and amplitude A x , then the dimensionless frequency of the punch oscillation, with respect to the moving substrate, is: If Ξ << 1, the rigid punch will get in contact with a material that had enough time to relax or that was never deformed; hence, we can drop the dependence on the punch relative position x−v d t R . Therefore, F μ and F k are written as a function of the following dimensionless groups: where we have introduced the groups v rel = v rel τ/R and y = y/R, being v rel = · x − v d the relative velocity between the punch and the substrate. When writing the friction force F μ , two contributions are considered: (i) the viscoelastic part of the friction coefficient μ visc ( v rel , y) originated by the internal dissipation in the viscoelastic material and (ii) the adhesive part of the frictional force that, as suggested by robust experimental outcomes [50][51][52][53], is assumed proportional to the apparent contact area A ( v rel , y) through a characteristic interfacial shear-strength τ 0 . Hence, the tangential force F μ ( v rel , y) is equal to: Notice that both μ visc ( v rel , y) and A ( v rel , y) are time varying quantities that depends at every time instant on the indentation of the spherical punch and on the relative horizontal velocity between the punch and the moving substrate [44,45,[47][48][49]. This is briefly recalled in Appendix A, where the main peculiarities of viscoelastic materials and the Boundary Element Method implemented to study the viscoelastic contact mechanics are summarized. Let us define, now, the dimensionless viscoelastic contact force F k = F k /E * 0 R 2 and the dimensionless contact area A = A/R 2 , where E * 0 = E 0 / 1 − ν 2 is the composite modulus, with E 0 the viscoelastic modulus of the substrate at zero frequency, i.e., E 0 = E (ω = 0), and ν the Poisson's ratio. As shown with more details in Appendix A, in fact, a generic linear viscoelastic material is characterized by the so-called rubbery modulus, that is E 0 , by the glassy modulus E ∞ , being the viscoelastic modulus for ω tending to infinite, and by a series of relaxation times τ . The simulations have been performed using a single relaxation time linear viscoelastic material with E 0 = 10 6 Pa , E ∞ = 10 7 Pa, τ = 0.01 s, R = 0.01 m and ν = 0.5. Such a choice is consistent with the approach adopted in the entire paper: the simple one-relaxation time material allows us to clearly see which are the possible mechanisms of instability in a viscoelastic material and how these depend on the relative velocity between the oscillator and the moving substrate, including the correct elastic limits at low and high relative speed. Although real materials have several relaxation times, we expect our results to be paradigmatic of the general behavior expected in these conditions. Furthermore, let us notice that, once the relaxation spectrum is known for a real viscoelastic material, the method presented here could be repeated to determine the stability boundaries for any given application.
Turning back to our simplified model, Fig. 3 shows: (a) the reaction force F k , (b) the contact area A and (c) the viscoelastic friction coefficient μ visc as a function of log 10 y and log 10 v rel . The black dots represent the numerical results obtained from the boundary element method analysis, while the shaded surfaces refer to the analytical functions used to fit the numerical results. Notice that the fitting laws for all the three quantities have been selected based on physical considerations and on the knowledge of the contact mechanics problem. Indeed, the advantage of using analytical fitting functions is twofold: (i) it guarantees a fast calling and evaluation of the numerical functions during the dynamical simulations and (ii) it allows any researcher to easily replicate or extend this study in the future.
Let us further comment on the physical ratio for the fitting equations employed for F k , A and μ visc (for more numerical detail on the fitting procedure the reader is referred to Appendix B). As for the viscoelastic friction coefficient, μ visc ( v rel , y) is proportional to the imaginary part of the viscoelastic modulus Im(E(ω)) and, thus, for a one relaxation time material, the dependence of μ visc on the speed v rel has a bell shape and is well described by a Gaussian function. The dependence on the indentation y can be easily recalled by using the Hertzian theory. The viscoelastic friction force F visc is proportional to the volume being deformed and, thus F visc ∝ a 2 y with a being the contact radius; in an Hertzian contact a = (Ry) 1/2 with R the radius of the spherical indenter, hence F visc ∝ y 2 . By recalling that the normal force is proportional to y 3/2 one obtains μ visc ∝ y 1/2 . Consequently, we have implemented the following fitting equation: which has provided a coefficient of determination r 2 = 0.9957 (see Fig. 3a). Notice that the viscoelastic friction monotonically increases with the indentation y/R as a larger penetration implies a larger volume that deforms and dissipates. Furthermore, with respect to the velocity dependence, the bell-shaped trend in the semi-logarithmic scale, where friction vanishes at very low and very high speed values, occurs as in the limit of very high and very low speeds, the material behaves elastically, whereas it reaches a maximum in the range of intermediate velocity.
With regard to the contact area A ( v rel , y), according to Hertz theory, the apparent contact area should depend linearly on the indentation y, whereas the dependence on velocity resemble a flipped Gaussian function. The fitting equation results: which gave a coefficient of determination r 2 = 0.9995 (see Fig. 3b). Indeed, the contact area A is constant at low and high velocity, where the material behaves elastically and the contact radius is only defined by the imposed indentation y as predicted by Hertzian theory, while, as shown in Ref. [44], due to viscoelasticity, it reaches a minimum in the range of intermediate velocities. Finally, the substrate reaction force should be proportional to ∝ y 3/2 , while the dependence on the velocity can be well captured by an erf (x) function. An excellent fit was obtained using the fit equation which provided a coefficient of determination r 2 = 0.9999 (see Fig. 3c). For a given penetration, the elastic reaction force increases moving from the rubbery to the glassy region, i.e., from low to high velocity.

Dimensionless formulation
After introducing the equilibrium equations in the horizontal and vertical directions and the equations that determine the contact forces, we develop here a dimensionless formulation, that will be used along the rest of the paper. Let us introduce the damping ratio ξ x and ξ y , referred to the motion components, respectively, along the −x axis and the −y axis, and, similarly, the natural frequencies ω nx , ω ny and ω nx y : By defining the dimensionless time θ and a reference displacement x 0 as we introduce the following dimensionless parameters: Using (9,11) and d dt = ω nx d dθ , the equilibrium equations can be written in dimensionless form as Notice that, in case of lift-off of the spherical punch, no contact between the punch and the substrate occurs, hence, in such conditions, F k ( v, y) =0 and μ visc ( v, y) = 0; Eqs. (12) will be modified accordingly: If not differently stated, in the next sections, we will assume the following model parameters: where the parameters in the second row may be useful if one wants to go back from dimensionless to dimensional quantities. Here, we are focusing our attention on the particular case Ω = η = 1: such a choice is carried out as we aim at investigating how Stribeck and mode-coupling instabilities co-exist and interact. Clearly, setting a vanishing (or infinite) value for Ω or η would give the system two far apart modes, thus preventing the mode-coupling instability to take place.

Linear stability analysis
Let us assume that the substrate is sliding at a given velocity v d and a steady state solution exists so that (12), the static equilibrium position ( x e , y e ) can be obtained by solving the following algebraic system of equations: (15) where y e = y e / R. By defining u (t) = x (t) − x e and q (t) = y (t) − y e , the system of equations (12) is linearized and written in matrix form as: Equation (16) can be rewritten in the classical form: where M, C, K are, respectively, the mass, damping and stiffness matrix, U T = { u, q} and c 22 = 2ξ y ; c 12 = 0; k 11 = 1 + η 2 cos 2 α; k 12 = k 21 = η 2 sin α cos α; The stability of the equilibrium state ( x e , y e ) can be determined by solving the eigenvalue problem det(A − λI) = 0 (21) being A the state matrix If all eigenvalues λ i have negative real parts (Re (λ i ) < 0, ∀ λ i ) the equilibrium solution ( x e , y e ) is linearly stable against small perturbations.

Equivalent damping in the horizontal direction
Let us preliminary focus on the damping matrix C and, in particular, on the search for possible sources of negative damping. As well known in literature, in fact, negative damping is one of the elements that trigger FIVs. When looking at matrix C, we realize that the component c 12 is null, while the term c 22 is positive and depends on the external structural source of damping. As for c 21 , this correlates the damping in the horizontal direction to the oscillations in the vertical motion: c 21 is always positive as F k increases monotonically with the sliding velocity (see Fig. 3c). Finally, the remaining term c 11 contains two contributions, one originated by the external damping ξ x , the other related to the dissipative processes occurring in the viscoelastic substrate. Specifically, we can define the system equivalent damping ratio ξ eq along the x-direction as ξ eq = c 11 /2, and, then, write ξ eq as: Now, in Fig. 4, for the model parameters previously fixed, the equivalent damping ratio ξ eq is plotted as a function of the dimensionless substrate velocity v d , for different values of the normal force F n = 7.5 * 10 −2 , 10 −3 , 10 −4 , 10 −5 and the tangential stress τ 0 = [0, 0.01, 0.05, 0.1].
Let us recall that the friction μ visc vanishes at low and high values of the sliding substrate velocity, i.e., when the substrate acts as elastic, respectively, in the rubbery and the glassy regimes; in these limits, the contact area A assumes the constant values given by the elastic Hertzian solution. Thus, in the limits of vanishing or infinite sliding velocity, the derivatives ∂μ/∂ v rel and ∂ A/∂ v rel go to zero and, thus, the equivalent damping coefficient is very small and equal to the linear damping coefficient ξ eq = ξ x , which is pretty small when compared to the interfacial damping. On the other hand, as soon as the speed is slightly increased, the viscoelastic friction markedly increases. This is the region where the transition from "stick to slip" takes place, which implies a large value for the derivative term ∂μ/∂ v rel and ultimately a large equivalent damping. The latter behavior is well-known in tribology, particularly for the rate and state friction models developed mostly by the geophysics community [54][55][56].
For very small values of τ 0 , the equivalent damping is uniquely defined by the viscoelastic friction (see, in all the panels in Fig. 4, the pale blue curves): starting from a value asymptotically tending to ξ x , ξ eq increases with the speed, reaches a maximum and then decreases, eventually becoming negative for v d being approximately equal to approximately v d ≈ 10 −2 as visible in Fig. 4 insets. Finally, for very large values of v d , ξ eq tends to ξ x .
Regarding the dependence on the interfacial shear stresses, by increasing τ 0 , the effect of ∂ A/∂ v rel becomes dominant. We employ τ 0 ∈ [0.01, 0.1] as this is a representative range for soft materials. Figure 4 shows that increasing τ 0 reduces ξ eq in the range of intermediate velocity. For low normal loads, this reduction may be strong enough to bring ξ eq in the negative half-plane. Although this may resemble counterintuitive, in Fig. 3b, it can be seen that the contact area A decreases with velocity, reaches a minimum and, only after that, increases again. This means that there is a speed range, where ∂ A/∂ v rel is negative and there occurs a negative damping effect. Ultimately, the faster the sphere goes, the smaller the contact area gets, the smaller the "adhesive" shear resistance is. Conversely, in the range of high velocity (see Fig. 4), where the contact area provides a positive damping (∂ A/∂ v rel > 0, see Fig. 3b), the higher τ 0 the higher ξ eq .

Linear stability for high/low velocity
Before proceeding with the stability map, it is interesting to consider the system under investigation in the limit of very low or very high driving velocity. In this limit, all the dissipative interfacial contributions vanishes as the material behaves elastically, and ∂μ/∂ v rel and ∂ A/∂ v rel goes to zero; there exist the externally imposed damping ratios ξ x and ξ y , but, as typical in most of the mechanical structures, these can be assumed very small (in the order of 0.01). It is hence licit to neglect the contribution of the damping matrix in the stability analysis (C = 0) to determine the system eigenvalues. With this approximation, it is possible to obtain, from the state matrix A, two pairs of complexconjugate eigenvalues: where λ 1,2 are purely imaginary, while λ 3,4 can be real or complex depending on the terms in the stiffness matrix. Incidentally, notice that we have neglected the contribution of the damping matrix, but the effect of the viscoelastic substrate still remains through the stiffness matrix K. While the terms k 11 and k 12 depend only on the structural stiffness, the terms k 12c and k 22 depend on the substrate properties and its deformation: although the layer, in this conservative limit of small and large speeds, behaves elastically, contact area and substrate restoring force still depend on the indentation. Figure 5 shows the system eigenvalues as a function of the interfacial shear strength τ 0 for v d = 10 −4 and F n = 7.5 * 10 −6 . At low interfacial shear strength, two purely imaginary eigenvalues exist, and then, at about τ 0crit ≈ 3.5 * 10 −2 the two modes merge and a pair of eigenvalues develop: one of these has a positive real part, showing the classical picture of a flutter type instability. In Fig. 6, the eigenvalues λ i are plotted as a function of the normal force F n for a given τ 0 = 0.01 and v d = 1.58 * 10 −4 . One recognizes that for F n 7.5 * 10 −7 there exist two purely imaginary eigenvalues (recall we set C = 0), while for a smaller normal force the two modes merge and a pair of eigenvalues, one with positive real part, appear. We conclude that both low normal forces and high interfacial shear stresses will promote flutter instability.
Looking at Fig. 5 and 6, it could be interesting to determine the critical interfacial shear strength τ 0,crit and the critical normal force F n,crit at which the flutter instability originates: we will, be able to define the "critical point", where the two modes merges. From Eq. (25), one obtains which can be cast in terms of critical shear strength τ 0,crit and critical normal load F n,crit : for shear stresses larger than τ 0,crit and normal load smaller than F n , self-excited vibrations will be triggered. These critical values are, respectively: where we introduce the dimensionless tangential stiffness k x = k x / E * 0 R . Figure 7 shows the critical shear strength τ 0,crit , above which the system is linearly unstable, as a function of the speed v d for several values of the applied normal force F n,crit . This critical shear stress τ 0,crit increases both with the substrate velocity and with the applied normal force. Notice that for a soft material a realistic value of τ 0,crit is in the range [0.01, 0.1]: hence, low normal forces can easily lead to flutter. Figure 8 shows the critical normal force F n,crit below which there exists flutter instability as a function of the substrate velocity v d and for different values of τ 0 = [0.01, 0.05, 0.1]. Ultimately, it is shown that F n,crit decreases with v d and increases with τ 0 : thus, for a given τ 0 a smaller velocity promotes instability.

Stability map for the full system
The dynamical system considered is linearly stable to small perturbations if max [Re (λ i )] < 0, being λ i the  27)), while the red solid line is the stability boundary obtained with the full system but setting η = 0. Panel (a) for F n = 10 −6 while panel (b) for F n = 10 −3 i−th eigenvalue obtained from Eq. (21). Here, the stability maps of the full system are shown for a given set of governing parameters: this has been carried out by accounting for the damping matrix C. Indeed, Fig. 9 shows the stability maps as a function of the driving velocity v d and of the interfacial shear stress τ 0 , for different values of the normal force F n . Now, for small values of the normal force and, specifically, for F n = 10 −6 , panel (a) shows that the stability boundary of the full system (obtained from Eq. (21), black dashed line in 9) is in good agreement with the estimate obtained neglecting the damping matrix C over all the velocity range, except in the range between v d 10 −6 ∼ 10 −4 . The agreement seems to worsen for larger normal forces, as clear in panel (b), where the stability analysis is carried out for F n = 10 −3 . We easily understand this observing that dissipation in the viscoelastic substrate is proportional to the volume of material deformed during the normal sliding, thus, at low normal forces, the only damping remaining in the system is related to ξ x and ξ y whose values are indeed quite small. For a higher normal force, (panel (b)) Eq. (27) gives an accurate result both at high and low driving velocity while it overestimates the stability range at intermediate velocity. This is immediately explained when one remembers the rheology of viscoelastic materials: they behave elastically (no dissipation occurs) at low and high frequency of excitation (glassy and rubbery region) while they dissipate energy for intermediate frequency (transition regime). It may result unexpected that the full system shows an instability region that is larger than the conservative counter-part. Nevertheless, Fig. 4 shows that the system equivalent damping (x-direction) is strongly affected by the contact area evolution with the driving velocity. In particular, negative damping is achieved for relatively small normal forces (see Fig. 4 bottom panels), which clearly foster instability. Furthermore, at very low interfacial shear strength ( τ 0 10 −3 ∼ 10 −4 ), the terms related to the contact area variation in the damping matrix C vanish and the only remaining contribution is related to the variation in the viscoelastic friction coefficient μ visc . Clearly, two different mechanisms of instability are supervening here: one related to the mode coupling instability, the other to the so-called Stribeck effect. To isolate the effect of the negative damping induced by the variation in the viscoelastic friction coefficient and of the contact area with the relative velocity, we show in Fig. 9 the stability boundary (red solid line) that one would obtain by setting η = 0. Indeed, for η = 0 the two modes are decoupled, which prevents the development of flutter instability. In the region of intermediate velocity, the prevailing mechanism for instability is due to a Stribeck effect. In this respect, the accuracy of Eq. (27) can be reconsidered: the analysis in the limit of low or high velocities, where viscoelastic dissipation vanishes, correctly predicts the stability boundary due to flutter in the region where the latter is the dominant mechanism for instability. Figure 10 shows the stability maps as a function of the driving velocity v d and of the applied normal force F n for two interfacial shear stresses τ 0 = [0.01, 0.1], respectively, in panel (a) and panel (b). As shown in Fig. 9, the stability boundary obtained by employing Eq. (28) for the conservative system (dashed black line) matches the full system stability boundary only at low and high sliding velocity. On the other end in the transition region the viscoelastic properties of the substrate strongly affect the system stability introducing the positive/negative damping contributions originated by the contact area variation and the friction coefficient evolution with the relative velocity. By comparing panel (a) and panel (b), one notices that the unstable region enlarges for materials with a higher interfacial shear strength: this is clearly determined by the terms related to the contact area variation in the damping matrix. This is demonstrate by the stability boundary we have obtained by setting η = 0: this again shows how, in the region of intermediate velocities, the prevailing mechanism for instability is due to a Stribeck's effect.

Post-instability behavior
To better understand the importance of accounting for the effect of a time varying contact area in determining the dynamical response of a mechanical system in relative motion with a viscoelastic substrate, here we focus on the post-instability response of the dynamical system by means of time integration results. An illustrative case is considered for F n = 7.5 * 10 −4 , v d = 5 * 10 −3 and a varying characteristic interfacial shear strength 10 −4 < τ 0 < 5 * 10 −2 . Since the Coulomb part of the frictional resistance is proportional to the interfacial shear strength, a vanishing τ 0 "suppresses" the effect of the varying contact area, while this gets stronger for high values of τ 0 . Figure 11 shows the root-mean- Increasing τ 0 completely changes the stability of the system. Indeed, for τ 0 2 * 10 −3 , the system is stable against small perturbations and the vibration amplitude is about 0. For larger τ 0 the system becomes unstable and the vibration amplitude X rms strongly increases with τ 0 . Clearly, accounting for the time varying contact area strongly influences the system dynamics. Figure 12 shows the dimensionless (a) horizontal displacement, (b) relative velocity, and (c) contact area as a function of the dimensionless time θ for the simulations in Fig. 11 with τ 0 = 10 −4 and τ 0 = 0.05 (respectively, dashed blue and solid black lines). For τ 0 = 10 −4 the system is linearly stable and { x, A, τ 0 } are constant. For τ 0 = 5 * 10 −2 a high amplitude limit cycle develops, characterized by "stick-slip" like motion. Notice that, during the "slip" phase ( v rel < 0), a rapid reduction in the contact area takes place, and this comes along with a sudden reduction in the Coulomb part of the frictional resistance, further strengthening the instability effect.

Conclusions
In this work, the linear stability analysis of a spherical oscillator moving on a viscoelastic substrate has been studied, by determining the contact interactions between the substrate and the punch through simulations based on a boundary element formulation. The frictional dissipation has been modeled adding up the viscoelastic contribution due to bulk hysteresis in the viscoelastic material to a Coulomb "adhesive" contribution that, consistently with recent measurements, has been assumed to be proportional to the apparent contact area via a characteristic interfacial shear strength.
Such a simple yet representative model has allowed us to investigate how two different instability mechanisms, that is, the Stribeck's one related to the falling characteristic of the friction law with the sliding velocity and the mode-coupling one, co-exist and compete in the context of viscoelastic contact. We have performed a linear stability analysis, which in the limit of low and high velocity, has provided analytical results in closed form, thus enabling us to estimate the critical point for flutter instability in terms of critical shear strength (normal force) above (below) which instability occurs. Comparisons with the full stability maps have shown that the critical points determined analytically are accurate in the limit of high/low velocity or when the applied normal force is small. Conversely, in the range of intermediate velocity, the analytical estimates are not accurate as we have shown that the dominant mechanism for instability is due to a Stribeck effect, i.e., it is dominated by the viscoelastic damping (positive or negative) introduced by the variation in friction coefficient and contact area with the relative velocity is.
Our results show how, the evolution of the apparent contact area at the interface may play a fundamental role in the developing of unstable vibrations and show once more how common simplifications of adopting a constant friction coefficient or a linear contact stiffness, can largely oversimplify the problem of detecting the parametric regions where mechanical systems experience self-sustained vibrations. This calls for further development of physics-based theoretical/numerical approaches capable of correctly accounting for the interactions that take place at the interface between contacting bodies. Availability of data and materials All data available upon request at the authors' email address.

Author contributions
Code availability Custom code available upon request at AP email address.

Declarations
Conflict of interest All the authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

Viscoelastic contact mechanics
A linear viscoelastic material responds to the following integral equation [57]: which relates the time-dependent strain ε (t) to the time derivative of the stress σ (t) by means of the creep function J (t). The latter satisfies causality, i.e., J (t < 0) = 0 and can be written as with H (t) being the Heaviside step function, while E 0 is the rubber elastic modulus of the material at zero-frequency, C (τ ) is a strictly positive function defined as the creep spectrum [57,58] and τ is the relaxation time, continuously distributed on the real axis. Obviously, in order to employ practically the creep function in a numerical formulation, we have to discretize the aforementioned relation and, so, rewrite J (t) as J (t) = H (t) 1/E 0 − n k=1 C k exp (−t/τ k ) . Now, it is possible to Fourier-transform the equations presented above by writing ε (ω) = σ (ω) /E (ω) with E (ω) being the viscoelastic modulus. The latter, defined as E (ω) = [iω J (ω)] −1 , can be written as: which, in the discretized form, becomes 1/E (ω) = 1/E ∞ + n k=1 C k /(1 + iωτ k ) . As shown in Fig. 13, at low frequencies, the viscoelastic material is in the elastic 'rubbery' region, where E(ω) becomes real and tends to E 0 , whilst for very large frequencies, E (ω) is again real, but this time tends to E ∞ and the material is in the so-called glassy region. Finally, in the intermediate frequency range, we have the proper viscoelastic region: the loss tangent ImE (ω) /ReE (ω) is, then, very large [see Fig. 13b]. Energy dissipation in the sliding contact will be, indeed, originated in this region of the viscoelastic spectrum . Now, let us briefly describe the Boundary Element formulation necessary to tackle the viscoelastic contact mechanics problem and determine the solution in terms of stresses, strains and, ultimately, friction. Basically, by recalling the translational invariance and the elasticviscoelastic correspondence principle [57], we can formulate the general linear-viscoelastic contact problem between a three-dimensional rigid indenter and a viscoelastic slab in terms of an integral Equation relating the normal surface displacement u of the viscoelastic solid and the normal interfacial stress σ : where η is the in-plane position vector referred to a frame moving with the rigid indenter and t is the time, while the integral kernel is given by J (t) and by the elastic Green's function G (η). Now, if we assume that the sliding occurs at constant velocity, thus implying that the rigid punch always slides over completely relaxed material regions, Eq. (31) can be re-written in the following form: where v is the speed by which the viscoelastic solid is deformed. In the contact configuration adopted in this paper (see Fig. 2), , the velocity v is then equal to v =v rel i with i being the unit vector along the x− direction. By employing the expression for the viscoelastic Green's function G (η, v) explicitly given in Ref. [44], the viscoelastic problem requires to solve Eq. [44] in the contact area. The latter can be be found by following the iterative numerical scheme aimed at determining the contact solution [49]: we discretize the contact domain in N square cells and assume that, in each square element, the normal stress σ is constantly equal to σ k = σ (η k ) where η k is the position vector of the center of the square cell D k . Thus, the normal displacement u i = u (η i ) at the center of the i-th square cell can be related to σ k with the following linear system: where each element of the response matrix L ik (v) is equal to L ik (v) = G η i −η k , v and parametrically depends on the velocity v. Eq. (33) can be easily solved and, then, to determine the contact area, the iterative procedure described in details in Ref. [49] can be adopted. Once the contact stresses and displacements are known, by recalling that v = vi , it is possible to compute the friction force as F f ric = Ω d 2 ησ (η)∂u/∂η.

Appendix -B
In this paper, as detailed in Sect. 2, we have fitted the contact solution in terms of friction force, normal force and contact area. Here, we report the fitting coefficients obtained from the fitting procedure (we used MAT-LAB). For the viscoelastic reaction force, the fitting function F k = a 1 y 3/2 1 2 a 2 + erf a 3 log 10 v + a 4 log 10 y + a 5 (34) with coefficients a 1 = 11.890; a 2 = 1.199; a 3 = 0.873; a 4 = −0.449; a 5 = 0.412; gave a coefficient of determination r 2 = 0.9999 (see Fig. 3a). For the contact area, the fitting function gave a coefficient of determination r 2 = 0.9957 (see Fig. 3c). Figure 14 compares the numerical results with the predictions obtained using the fitting equations proposed above. In each panel, the coefficient of determination r 2 is reported.