Friction-Induced Vibration Suppression via the Tuned Mass Damper: Optimal Tuning Strategy

Abstract

Friction-induced vibrations are a significant problem in various engineering applications, while dynamic vibration absorbers are an economical and effective tool for suppressing various kinds of vibrations. In this study, the archetypal mass-on-moving-belt model with an attached dynamic vibration absorber was considered. By adopting an analytical procedure, the optimal tuning of the absorber’s parameters was defined. Furthermore, the bifurcations occurring at the loss of stability were analytically investigated; this analysis illustrated that a properly chosen nonlinearity in the absorber’s stiffness permits controlling the supercritical or subcritical character of the bifurcation. However, a numerical analysis of the system’s dynamics, despite confirming the analytical results, also illustrated that the system’s global behavior is only slightly affected by the bifurcation character. Indeed, a dynamic vibration absorber possessing a perfectly linear restoring force function seems to provide the optimal performance; namely, it minimizes the velocity range for which stick–slip oscillations exists.

1. Introduction

Friction-induced vibrations (FIVs) are a peculiar type of oscillations generated by the friction acting between two bodies in relative motion. They consist of either the successions of stick and slip phases between the two bodies [1], or of quasi-harmonic oscillations having an approximately sinusoidal displacement–time relation [2,3]. Although for some specific applications these kinds of vibrations are intentionally generated, such in the case of violin strings [4] and singing wine glasses [5], typically they are seen as a detrimental phenomenon, as in the cases of brake squeal [6] or earthquakes [7,8].

Several methods exist for suppressing FIVs. One possibility is to reduce the friction force at the interface utilizing a lubricant. This method is efficient if friction is not required for the device to operate, such as in the case of hinge squeaking; however, it cannot be adopted for brake squeal mitigation, where high friction is strictly required. Most brake squeal suppression methods consist of increasing the system damping, which is obtained with various techniques [6,9]. Experimental observations also illustrated that isolating the natural frequencies of the brake system’s components at low frequencies tends to reduce the occurrence of audible brake squeal [10]; however, in many cases, this strategy is not effective [11]. For active methods for suppressing FIV, Cunefare and Graf [12] proposed adopting a dither exciting the system at non-audible frequencies, which can suppress brake squeal. Papangelo and Ciavarella [13] proposed to mitigate FIVs by normal load variation, for which they provided a closed-form solution.

The dynamic vibration absorber (DVA) is a practical tool for suppressing undesired vibrations in several engineering applications. Its classical design [14] consists of a mass attached to the host structure through a spring and a damper. By tuning its natural frequency in correspondence of the frequency to be damped, it is able to dynamically interact with the host structure dissipating vibration energy. It is successfully employed in several engineering fields for the suppression of various kinds of vibrations, such as flutter instabilities [15,16], rolling motions in ships [17], helicopter rotor oscillations [18] and machine tool vibrations [19,20]. Although DVAs are a mature technology, which was first proposed more than one hundred years ago [17], to the authors’ knowledge, there are only a few and relatively recent studies addressing its implementation to suppress FIVs. Popp and Rudolph [21] numerically and experimentally analyzed the performance of a DVA for FIV suppression; by utilizing a single-degree-of-freedom (DoF) primary system, they illustrated its beneficial effect. Chatterjee [22] studied the stability properties of an undamped DVA attached to a two-DoF primary system. Very recently, Niknam and Farhang [23] proposed a study similar to that of Chatterjee [22], where they also provided some numerical simulations of the full system, missing in [22]. Despite the promising results obtained in [21,22,23], a clear tuning strategy of the DVA’s parameters for maximizing its performance is still missing. This paper aims to fill this gap by providing a precise tuning of the absorber parameters for optimizing stability properties and studying the behavior of the host system with the attached DVA while stability is lost.

The rest of the paper is organized as follows. In Section 2, the mechanical model, consisting of the host mass-on-moving-belt system and the attached DVA, is introduced. In Section 3, the stability analysis of the host system, without and with the DVA, is performed, providing explicit equations for the optimal tuning of the absorber parameters. In Section 4 and Section 5, the bifurcations occurring at the loss of stability of the host system, without and with absorber, are analytically studied. Furthermore, the effect of the addition of a cubic term in the absorber’s restoring force is analytically investigated; analytical results are integrated by numerical simulation, illustrating the system’s behavior at high amplitudes. In Section 6, conclusions about the benefits and limitations of the DVA are presented.

2. Mechanical Model

2.1. Primary System

The primary system considered in this study is the classical mass-on-moving-belt model, which is an archetypal system for studying FIVs [1,4,21]. As shown in Figure 1, this single-DoF system consists of a mass $m_1$, a linear spring $k_1$ and a linear damper $c_1$. The mass of the system is in contact with the belt, which moves at a constant driving speed v, while the friction coefficient $\mu \left(v_{\mathrm{rel}}\right)$ of the contact is a function of the relative velocity $v_{\mathrm{rel}}=v-{\dot{x}}_1$.

Figure 1b illustrates the free body diagram (FBD) of the host system. The forces acting upon the lumped mass $m_1$ are the normal forces $F_{\mathrm{N}}$ that cancel each other; the damping and spring forces $F_{c1}$ and $F_{k1}$, respectively; and the friction force $F_{\mathrm{R}}$. (For the represented FBD, it is assumed that $v_{\mathrm{rel}}\ge 0$.) Based on Newton’s second law of motion and considering that $F_{\mathrm{R}}=\mu \left(v\mathrm{rel}\right)F_{\mathrm{N}}$, $F_{c1}=c_1{\dot{x}}_1$ and $F_{k1}=k_1{x}_1$, we obtain that the second order differential equation describing the motion of the one DoF system is

$$m_1{\ddot{x}}_1+c_1{\dot{x}}_1+k_1{x}_1=F_{\mathrm{R}},$$

(1)

with

$$\left\{\begin{array}{cc}{F}_{\mathrm{R}}=\mu \left(v_{\mathrm{rel}}\right)F_{\mathrm{N}}\hfill & v_{\mathrm{rel}}\ne 0\hfill \\ |F_{\mathrm{R}}|\le {\mu }_{\mathrm{s}}{F}_{\mathrm{N}}\hfill & v_{\mathrm{rel}}=0,\hfill \end{array}\right.$$

(2)

where the overdot indicates derivation with respect to the time t. The system is in stick condition when the relative velocity is zero ($v_{\mathrm{rel}}=0$). In this case, the friction force is smaller or equal to ${\mu }_{\mathrm{s}}{F}_{\mathrm{N}}$, where ${\mu }_{\mathrm{s}}$ is the static friction. In sliding condition ($v_{\mathrm{rel}}\ne 0$), the direction of the friction force $F_{\mathrm{R}}$ depends on the sign of the relative velocity $v_{\mathrm{rel}}$, which is included in the mathematical formulation of $\mu \left(v_{\mathrm{rel}}\right)$. Additional details about the friction coefficient utilized are provided in Section 2.3. Let us introduce the following expressions

$${\zeta }_1=\frac{c_1}{2\sqrt{m_1{k}_1}}; {\omega }_{\mathrm{n}1}=\sqrt{\frac{k_1}{m_1}} ; x_0=\frac{F_{\mathrm{N}}}{k_1} ; \tau ={\omega }_{\mathrm{n}1}t .$$

(3)

By dividing Equation (1) by the mass $m_1$, applying the expressions from Equation (3) and dividing it by ${\omega }_{\mathrm{n}1}^2$, we obtain

$$x_1^{″}+2{\zeta }_1{x}_1^{\prime }+x_1=\frac{F_{\mathrm{R}}}{k_1} ,$$

(4)

where prime ${}^{\prime }$ indicates derivation with respect to the dimensionless time $\tau$. Then, introducing the dimensionless displacement ${\tilde{x}}_1=x_1/x_0$, the system is eventually reduced to

$${\tilde{x}}_1^{″}+2{\zeta }_1{\tilde{x}}_1^{\prime }+{\tilde{x}}_1={\tilde{F}}_{\mathrm{R}} ,$$

(5)

where ${\tilde{F}}_{\mathrm{R}}=F_{\mathrm{R}}/F_{\mathrm{N}}$.

2.2. Mechanical Model of the Host System with the Absorber

We now attach a DVA to the host system. The basis of the model is the same as the one mentioned in the previous subsection. The only additional element is the absorber mass $m_2$, which is attached to the primary system through a spring and a linear damper. The schematic depiction of this two-DoF system is provided in Figure 2.

The absorber’s spring encompasses a linear term and a cubic term. The differential equations describing the dynamics of the system are

$$\begin{array}{cc}& m_1{\ddot{x}}_1+c_1{\dot{x}}_1+k_1{x}_1+c_2({\dot{x}}_1-{\dot{x}}_2)+k_2(x_1-x_2)+k_{\mathrm{nl}2}{(x_1-x_2)}^3=F_{\mathrm{R}}\hfill \\ & m_2{\ddot{x}}_2+c_2({\dot{x}}_2-{\dot{x}}_1)+k_2(x_2-x_1)+k_{\mathrm{nl}2}{(x_2-x_1)}^3=0 ,\hfill \end{array}$$

(6)

where $k_2$ and $k_{\mathrm{nl}2}$ are the linear and cubic coefficients of the absorber stiffness, respectively, while $c_2$ is the linear coefficient of the absorber damping. Let us introduce the following expressions

$${\zeta }_2=\frac{c_2}{2\sqrt{m_2{k}_2}} ; {\omega }_{\mathrm{n}2}=\sqrt{\frac{k_2}{m_2}} ; \epsilon =\frac{m_2}{m_1} ; \gamma =\frac{{\omega }_{\mathrm{n}2}}{{\omega }_{\mathrm{n}1}} .$$

(7)

By diving Equation (6) by $m_1$, applying the expressions from (7) and (3), dividing by ${\omega }_{\mathrm{n}1}^2$ and utilizing the dimensionless time $\tau$ and dimensionless displacements ${\tilde{x}}_1$ and ${\tilde{x}}_2=x_2/x_0$, we obtain

$$\begin{array}{cc}& {\tilde{x}}_1^{″}+2{\zeta }_1{\tilde{x}}_1^{\prime }+{\tilde{x}}_1+2\epsilon {\zeta }_2\gamma ({\tilde{x}}_1^{\prime }-{\tilde{x}}_2^{\prime })+\epsilon {\gamma }^2({\tilde{x}}_1-{\tilde{x}}_2)+\epsilon {\kappa }_{\mathrm{nl}2}{({\tilde{x}}_1-{\tilde{x}}_2)}^3={\tilde{F}}_{\mathrm{R}}\hfill \\ & \epsilon {\tilde{x}}_2^{″}+2\epsilon {\zeta }_2\gamma ({\tilde{x}}_2^{\prime }-{\tilde{x}}_1^{\prime })+\epsilon {\gamma }^2({\tilde{x}}_2-{\tilde{x}}_1)+\epsilon {\kappa }_{\mathrm{nl}2}{({\tilde{x}}_2-{\tilde{x}}_1)}^3=0 ,\hfill \end{array}$$

(8)

where ${\kappa }_{\mathrm{nl}2}=k_{\mathrm{nl}2}{x}_0^2/\left(k_1\epsilon \right)$. A variable change is performed, where ${\tilde{x}}_3={\tilde{x}}_1-{\tilde{x}}_2$ (relative displacement of $m_2$), hence Equation (8) transforms into

$$\begin{array}{cc}& {\tilde{x}}_1^{″}+2{\zeta }_1{\tilde{x}}_1^{\prime }+{\tilde{x}}_1+2\epsilon {\zeta }_2\gamma {\tilde{x}}_3^{\prime }+\epsilon {\gamma }^2{\tilde{x}}_3+\epsilon {\kappa }_{\mathrm{nl}2}{\tilde{x}}_3^3={\tilde{F}}_{\mathrm{R}}\hfill \\ & \epsilon \left(\right({\tilde{x}}_1^{″}-{\tilde{x}}_3^{″})-2{\zeta }_2\gamma {\tilde{x}}_3^{\prime }-{\gamma }^2{\tilde{x}}_3-{\kappa }_{\mathrm{nl}2}{\tilde{x}}_3^3)=0 .\hfill \end{array}$$

(9)

2.3. Friction Force

The applied friction law is described by an exponential decaying function, as done for instance in [1], i.e.,

$$\mu \left(v_{\mathrm{rel}}\right)=({\mu }_{\mathrm{d}}+({\mu }_{\mathrm{s}}-{\mu }_{\mathrm{d}}){\mathrm{e}}^{-{\displaystyle \frac{|v_{\mathrm{rel}}|}{{\tilde{v}}_0}}})\mathrm{sgn}(v_{\mathrm{rel}}) ,$$

(10)

where the relative velocity is $v_{\mathrm{rel}}=v-{\tilde{x}}_1^{\prime }$. The values assumed by $\mu$ for a range of relative velocities for ${\mu }_{\mathrm{s}}=1$, ${\mu }_{\mathrm{d}}=0.5$ and $v_0=0.5$ are represented in Figure 3. The values for the friction law adopted in the present study are the same utilized in [1]. Nevertheless, as illustrated below, the optimization of the absorber parameters does not strictly depend on the considered friction law.

3. Linear Stability Analysis

We now aim at analyzing the system’s behavior for small perturbations around the equilibrium. Thus, we linearize it around its trivial solution and we study the stability of this trivial solution. The analysis is performed for both the systems, with and without DVA, in order to assess the beneficial effects of the DVA.

3.1. Linear Stability of the Host System without the DVA

We first linearize Equation (5) around its equilibrium point $x_{1\mathrm{e}}=\mu (v_{\mathrm{rel}}=v)$, obtaining

$$z_1^{″}-2\psi z_1^{\prime }+z_1=0 ,$$

(11)

where $z_1=x_{1\mathrm{e}}+{\tilde{x}}_1$ and $2\psi =\partial \mu /\partial z_1^{\prime }{|}_{z_1^{\prime }=0}-2{\zeta }_1$. Considering the friction law adopted, we have that

$$\frac{\partial \mu }{\partial z_1^{\prime }}{|}_{z_1^{\prime }=0}=\frac{{\mu }_s-{\mu }_d}{v_0}{\mathrm{e}}^{-\frac{v}{v_0}} ,$$

(12)

(valid for $v>0$), which is the slope of the friction force coefficient at the belt velocity v.

Equation (11) corresponds to a linear oscillator, whose trivial solution is asymptotically stable if and only if $\psi <0$. According to the friction law utilized, and considering that $\psi$ is a monotonically decreasing function of v, the trivial solution is stable for

$$v>v_{\mathrm{h},\mathrm{cr}}=v_{0}ln\left({\displaystyle \frac{{\mu }_{\mathrm{s}}-{\mu }_{\mathrm{d}}}{2v_0{\zeta }_1}}\right) .$$

(13)

This result is well-known and better discussed, for instance, in [1]. The practical consequences are that, if the belt moves at a speed lower than $v_{\mathrm{h},\mathrm{cr}}$, the equilibrium of the system is unstable and stick–slip oscillations occur. More details about these stick–slip oscillations are provided below.

3.2. Linear Stability of the Host System with DVA

To study the stability of the system with the DVA, we linearize Equation (9) around the equilibrium ${\tilde{x}}_1=x_{1\mathrm{e}}$ and ${\tilde{x}}_3=0$. By reformulating Equation (9) in explicit form with respect to ${\tilde{x}}_1^{″}$ and ${\tilde{x}}_3^{″}$ and by utilizing the variables and parameters introduced in the previous subsection, we obtain

$$\left[\begin{array}{c}{z}_1^{\prime }\\ z_2^{\prime }\\ z_3^{\prime }\\ z_4^{\prime }\end{array}\right]=\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ -1& -{\gamma }^2\epsilon & 2\psi & -2\epsilon {\zeta }_2\gamma \\ -1& -{\gamma }^2(\epsilon +1)& 2\psi & -2(\epsilon +1)\gamma {\zeta }_2\end{array}\right]\left[\begin{array}{c}{z}_1\\ z_2\\ z_3\\ z_4\end{array}\right]=\mathit{Az} ,$$

(14)

where $z_1=x_{1\mathrm{e}}+{\tilde{x}}_1$, $z_2=x_3$, $z_3=z_1^{\prime }$ and $z_4=z_2^{\prime }$.

We analyze the characteristic exponents of the system to determine the stability. The characteristic polynomial is

$$p\left(\lambda \right)=det(\mathit{A}-\lambda \mathit{I})=0 .$$

(15)

The characteristic polynomial is in the form of

$$p\left(\lambda \right)=a_0{\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4=0 (a_0>0) .$$

(16)

The stability is determined based on the Liénard–Chipart conditions (LCC) [24]. For the polynomial $p\left(\lambda \right)$ to have all roots with negative real parts, it is necessary and sufficient that the coefficients of the polynomial are all positive ($a_i>0,i=1,\cdots ,4$) and that the determinant inequalities ${\Delta }_1>0$ and ${\Delta }_3>0$ (Hurwitz determinantal inequalities) are verified, which in our case means

$$\begin{array}{cc}\hfill a_0& =1>0 ,\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill a_1& =2\left(\gamma \right(1+\epsilon ){\zeta }_2-\psi )>0⟺\psi <(1+\epsilon )\gamma {\zeta }_2 ,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill a_2& =1+\gamma (\gamma (1+\epsilon )-4{\zeta }_2\psi )>0⟺\psi <{\displaystyle \frac{1+{\gamma }^2(1+\epsilon )}{4\gamma {\zeta }_2}} ,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill a_3& =2\gamma ({\zeta }_2-\gamma \psi )>0⟺\psi <{\displaystyle \frac{{\zeta }_2}{\gamma }} ,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill a_4& ={\gamma }^2>0 ,\hfill \end{array}$$

(21)

$${\Delta }_1=a_1=2\gamma (1+\epsilon ){\zeta }_2-2\psi >0 (\mathrm{Already} \mathrm{present} \mathrm{in} (18\left)\right) ,$$

(22)

$$\begin{array}{cc}\hfill {\Delta }_3& =a_1{a}_2{a}_3-a_0{a}_3^2-a_1^2{a}_4=-4\gamma \left({\gamma }^4{(\epsilon +1)}^2{\zeta }_2\psi \right.-{\gamma }^3{\psi }^2(4(\epsilon +1){\zeta }_2^2+\epsilon )\hfill \\ & +2{\gamma }^2{\zeta }_2\psi (2(\epsilon +1){\zeta }_2^2+2{\psi }^2-1)\left.-\gamma {\zeta }_2^2(\epsilon +4{\psi }^2)+{\zeta }_2\psi \right) >0 .\hfill \end{array}$$

(23)

3.2.1. Analytical Optimal Solution

Considering the stability analysis performed in the previous section, we aim at finding the parameter values of the absorber which maximize the stable region. The linear system in Equation (14) has the same mathematical form as the one studied in [25]. Thus, we can follow the same steps in order to optimize the absorber.

First, we look at the curves where the coefficients and the Hurwitz determinants are zeros; these are the boundaries where certain roots change. At certain boundaries, the stable/unstable transition takes place. Following the steps discussed in [25], we can define specific points on these curves, which helps us find the optimal parameters. Figure 4 shows the stability regions for different values of ${\zeta }_2$; the boundary curves of the LCC are also depicted. The gray shaded region is the stable region; we can observe that the stability boundary is not at a constant value of $\psi$ as it is the case for the host system without the DVA; instead, it is a function of $\gamma$, with a pronounced peak for $\gamma \approx 1$. This is an expected feature, considering that the DVA usually needs to be tuned at a frequency close to the natural frequency of the primary system [14]. For different values of ${\zeta }_2$, the maximum value of $\psi$ also changes, thus we need to find the optimal combination of $({\zeta }_2,\gamma )$ such that the value of $\psi$ generating instabilities is maximized. For the optimization, the mass ratio $\epsilon$ is assumed constant; however, the results of the analysis show that larger values of $\epsilon$ increase the stable region.

As we can see in Figure 4a, the intersection of the curves $a_1={\Delta }_1=0$ and $a_3=0$ defines a point that we denote with C. Its coordinates in the $(\psi ,\gamma )$ plane are C $=({\zeta }_2\sqrt{1+\epsilon },1/\sqrt{1+\epsilon })$. For low values of ${\zeta }_2$, Point C marks the maximal value of $\psi$ providing stability, which we call ${\psi }^{*}$. Increasing the value of ${\zeta }_2$, Figure 4c illustrates that Point C does not identify any more ${\psi }^{*}$ and that the stability boundary is defined by Equation ${\Delta }_3=0$. Nevertheless, also for large damping ${\zeta }_2$, the value of $\gamma$ which maximizes ${\psi }^{*}$ is very close to the $\gamma$ coordinate of Point C. Therefore, we define Point P as the point of the curve ${\Delta }_3=0$ having the same $\gamma$ coordinate as C. Substituting the $\gamma$ coordinate of C, $\gamma =1/\sqrt{1+\epsilon }$, into ${\Delta }_3=0$, we obtain Point P $=(\epsilon /(4{\zeta }_2\sqrt{1+\epsilon }),1/\sqrt{1+\epsilon })$.

Since Points C and (approximately) P alternatively mark ${\psi }^{*}$, depending on the value of the absorber damping ${\zeta }_2$, by choosing ${\zeta }_2$ such that P and C coincide, we can maximize ${\psi }^{*}$. Imposing equality between the $\psi$ coordinates of C and P, we attain

$${\zeta }_2\sqrt{1+\epsilon }=\frac{\epsilon }{\left(4\sqrt{1+\epsilon }{\zeta }_2\right)} ,$$

(24)

therefore, ${\psi }^{*}$ is maximized for

$$\gamma ={\gamma }_{\mathrm{opt}}={\displaystyle \frac{1}{\sqrt{1+\epsilon }}} \mathrm{and} {\zeta }_2={\zeta }_{2\mathrm{opt}}={\displaystyle \frac{1}{2}}\sqrt{{\displaystyle \frac{\epsilon }{1+\epsilon }}} ,$$

(25)

and the corresponding maximal value of ${\psi }^{*}$ is

$${\psi }_{\max }={\displaystyle \frac{\sqrt{\epsilon }}{2}} .$$

(26)

Accordingly, if $\gamma ={\gamma }_{\mathrm{opt}}$ and ${\zeta }_2={\zeta }_{2\mathrm{opt}}$, the system is stable if $2\psi <\sqrt{\epsilon }$, or, in other words, if

$${\displaystyle \frac{{\mu }_{\mathrm{s}}-{\mu }_{\mathrm{d}}}{v_0}} {\mathrm{e}}^{-\frac{v}{v_0}}-2{\zeta }_1<\sqrt{\epsilon } .$$

(27)

Considering that $\psi$ is a monotonically decreasing function of v, the equilibrium of the system is stable if

$$v>v_{\mathrm{cr}}=v_{0}log\left({\displaystyle \frac{{\mu }_{\mathrm{s}}-{\mu }_{\mathrm{d}}}{v_0(\sqrt{\epsilon }+2{\zeta }_1)}}\right) .$$

(28)

The stability chart corresponding to this case is illustrated in Figure 4b. We remark that the optimal tuning proposed here does not strictly depend on the friction law considered, which could be modeled with alternative functions without varying ${\gamma }_{\mathrm{opt}}$ and ${\zeta }_{2\mathrm{opt}}$. This represent a clear advantage in the case of real engineering applications.

3.2.2. Numerical Validation

The optimization procedure utilized in the previous section is based on a heuristic approach, which does not prove that ${\gamma }_{\mathrm{opt}}$ and ${\zeta }_{2\mathrm{opt}}$ provide the maximal possible value of ${\psi }^{*}$. Therefore, its validity should be verified numerically. The numerical analysis is performed by directly computing eigenvalues of matrix $\mathit{A}$ on a grid of the $({\zeta }_2,\gamma ,\psi )$ space for a fixed $\epsilon$ value ($\epsilon =0.05$) and identifying the couple of values ${\zeta }_2,\gamma$ which provides the maximal ${\psi }^{*}$. After several trials, the analysis is finally performed on the grid described in

Table 1. Parameter grid for the optimum search.

Parameter	min	max	Step
${\zeta }_2$	0.1	0.11	${10}^{-5}$
$\gamma$	0.97589	0.97591	${10}^{-6}$
$\psi$	0.1124	0.1126	${10}^{-5}$

.

This analysis provided the optimal values for ${\zeta }_2$ and $\gamma$, which are indicated in

Table 2. Comparison of numerical and analytical results.

Parameter	Numerical	Analytical	Relative Error [%]
${\zeta }_2$	0.10977	0.109109	0.605867 [%]
$\gamma$	0.975899	0.975900	0.000109945 [%]
$\psi$	0.11247	0.111803	0.596226 [%]
$v_{\mathrm{cr}}$	1.36678	1.36883	0.364409 [%]

and directly compared with the optimal values obtained analytically.

Although numerical and analytical optimal parameters do not coincide, their difference is minimal and negligible for most engineering applications. In particular, the optimal $\gamma$ value is practically the same in both cases. The critical velocity, computed utilizing the parameter values indicated in

Table 3. Numerical input data.

${\mathit{\mu }}_{\mathbf{s}}$	${\mathit{\mu }}_{\mathbf{d}}$	${\mathit{v}}_0$	$\mathit{\epsilon }$	${\mathit{\zeta }}_1$
1	0.5	0.5	0.05	0.05

, has a difference of less than 0.4% in the two cases. We remark that the difference between numerical and analytical computation is not related to the inaccuracy of stability estimation through the LCC, which exactly predicts an equilibrium’s stability, but to the heuristic approach utilized for the optimization.

Figure 5, illustrating the curve ${\Delta }_3=0$ for various values of ${\zeta }_2$, enables us to understand the reason for the difference between the results obtained with the numerical and analytical approach. The blue line in the figure corresponds to the analytical optimization, while the green line to the numerical one. The yellow and red curves refer to values of ${\zeta }_2$ slightly higher and lower than the optimal ones, respectively. The inaccuracy of the analytical procedure is due to the fact that the peak of the ${\Delta }_3=0$ curve does not exactly lie on Point P (which has a fixed $\gamma$ value and it is not represented in the figure). However, considering the minimal difference found and the practical compactness of Equations (24) and (25), these will be utilized in the continuation of the paper.

3.3. Evaluation of the Absorber’s Performance

We now aim at evaluating the performance of the DVA by performing various comparison between the system with and without DVA. As mentioned above, the critical velocities for both cases are

$$\begin{array}{cc}\hfill v_{\mathrm{h},\mathrm{cr}}& =v_{0}ln\left({\displaystyle \frac{{\mu }_{\mathrm{s}}-{\mu }_{\mathrm{d}}}{2v_0{\zeta }_1}}\right) ,\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill v_{\mathrm{cr}}& =v_{0}ln\left({\displaystyle \frac{{\mu }_{\mathrm{s}}-{\mu }_{\mathrm{d}}}{v_0(\sqrt{\epsilon }+2{\zeta }_1)}}\right) .\hfill \end{array}$$

(30)

Considering the characteristics of a logarithm function, we can state that, if the argument becomes 1, the logarithm function yields 0. Thus, there is a certain parameter set for which the critical velocity becomes 0 (inherent stability). We define ${\zeta }_{1,\mathrm{cr}}$ as the critical primary damping parameter, for which the critical velocity becomes 0. To obtain ${\zeta }_{1,\mathrm{cr}}$, we solve the arguments of the logarithm functions for 1; these yield

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\mu }_{\mathrm{s}}-{\mu }_{\mathrm{d}}}{2v_0{\zeta }_{1\mathrm{h},\mathrm{cr}}}}& =1⟶{\zeta }_{1\mathrm{h},\mathrm{cr}}={\displaystyle \frac{{\mu }_{\mathrm{s}}-{\mu }_{\mathrm{d}}}{2v_0}} ,\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\mu }_{\mathrm{s}}-{\mu }_{\mathrm{d}}}{v_0(\sqrt{\epsilon }+2{\zeta }_{1,\mathrm{cr}})}}& =1⟶{\zeta }_{1,\mathrm{cr}}={\displaystyle \frac{{\mu }_{\mathrm{s}}-{\mu }_{\mathrm{d}}-v_0\sqrt{\epsilon }}{2v_0}} .\hfill \end{array}$$

(32)

Utilizing values in Table 3, the numerical values for the critical primary damping are

$$\begin{array}{cc}\hfill {\zeta }_{1\mathrm{h},\mathrm{cr}}& =0.5 ,\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill {\zeta }_{1,\mathrm{cr}}& =0.388197 .\hfill \end{array}$$

(34)

This shows that the application of an optimally tuned DVA with a mass of only 5% of the host system mass enables to reduce the critical primary damping of 22%.

Considering, instead, the critical velocity as a base of comparison, we define the improvement factor $\phi$ such that

$$\phi :={\displaystyle \frac{{\tilde{v}}_{\mathrm{h},\mathrm{cr}}-{\tilde{v}}_{\mathrm{cr}}}{{\tilde{v}}_{\mathrm{h},\mathrm{cr}}}}\times 100\% .$$

(35)

The critical velocities for both cases and the improvement factor are plotted against the varying ${\zeta }_1$ in Figure 6, utilizing the parameter values in Table 3. We can observe that the difference in critical velocity is more significant for smaller values of ${\zeta }_1$, i.e., for a slightly damped host system. We also notice that, if the host system is completely undamped, then the critical velocity is undefined, meaning that the equilibrium is always unstable. For any value of the host system damping ${\zeta }_1$, the improvement factor is almost always above 50%.

Let us also observe what happens if we vary the value of the mass ratio $\epsilon$. Similar to before, the critical velocities for both systems and the improvement factor curve are illustrated in Figure 7. The parameter values are those indicated in Table 3. The critical velocity of the host system is constant because it does not depend on $\epsilon$; however, for the system with the DVA, the critical velocity monotonously decreases with $\epsilon$. Utilizing the parameter values in Table 3, the critical velocities for both the host system and the system with DVA are $v_{\mathrm{h},\mathrm{cr}}=1.151$ and $v_{\mathrm{cr}}=0.5641$, hence the improvement provided by the DVA is of $\phi =51\%$.

4. Bifurcation Analysis of the Host System without the DVA

The analysis performed in Section 3 refers to the system linearized around its equilibrium. Therefore, it is able to describe its dynamics only in the vicinity of the equilibrium, while phenomena occurring when the stability is lost are overlooked. Additionally, it provides no information about the stable equilibrium’s robustness if the system is subject to non-small perturbations. To investigate the system behavior at the loss of stability and correctly evaluate the DVA performance, we reintroduce the nonlinear terms and analytically perform a bifurcation analysis of the system without and with DVA.

Considering the system in Equation (5), we first center the system around its equilibrium point $x_{1\mathrm{e}}$ by introducing the variable $z_1=x_{1\mathrm{e}}+{\tilde{x}}_1$, and then we expand it in Taylor series up to the third order, obtaining

$$\left[\begin{array}{c}{z}_1^{\prime }\\ z_3^{\prime }\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -1& \frac{{\mu }_d-{\mu }_s}{v_0}{\mathrm{e}}^{-\frac{v}{v_0}}-2{\zeta }_1\end{array}\right]\left[\begin{array}{c}{z}_1\\ z_3\end{array}\right]+\left[\begin{array}{c}0\\ -\left(\frac{{\mu }_d-{\mu }_s}{2v_0^2}{\mathrm{e}}^{-\frac{v}{v_0}}\right)z_3^2-\left(\frac{{\mu }_d-{\mu }_s}{6v_0^3}{\mathrm{e}}^{-\frac{v}{v_0}}\right)z_3^3\end{array}\right]={\mathit{A}}_{\mathrm{h}}{\mathit{z}}_{\mathrm{h}}+{\mathit{b}}_{\mathrm{h}} .$$

(36)

For $v\approx v_{\mathrm{h},\mathrm{cr}}$, matrix ${\mathit{A}}_{\mathrm{h}}$ has complex conjugate eigenvalues ${\lambda }_{1,2\mathrm{h}}={\alpha }_{1\mathrm{h}}\pm \mathrm{i}{\omega }_{1\mathrm{h}}$ and eigenvectors ${\mathit{s}}_1={\overline{\mathit{s}}}_2$, which are reduced to

$${\lambda }_{1\mathrm{h}}=\mathrm{i}, {\lambda }_{2\mathrm{h}}=-\mathrm{i}, {\mathit{s}}_1=\left[\begin{array}{c}-\mathrm{i}\\ 1\end{array}\right], {\mathit{s}}_2=\left[\begin{array}{c}\mathrm{i}\\ 1\end{array}\right]$$

(37)

for $v=v_{\mathrm{h},\mathrm{cr}}$.

We then define the transformation matrix

$${\mathit{T}}_{\mathrm{h}}=\left[\begin{array}{cc}\mathrm{Re}\left({\mathit{s}}_2\right)& \mathrm{Im}\left({\mathit{s}}_2\right)\end{array}\right]{|}_{v=v_{\mathrm{h},\mathrm{cr}}} ,$$

(38)

we apply the transformation $\mathit{z}={\mathit{T}}_{\mathrm{h}}{\mathit{y}}_{\mathrm{h}}$ and we pre-multiply Equation (36) by ${\mathit{T}}_{\mathrm{h}}^{-1}$, leading to

$${\mathit{y}}_{\mathrm{h}}^{\prime }={\mathit{T}}_{\mathrm{h}}^{-1}{\mathit{A}}_{\mathrm{h}}{\mathit{T}}_{\mathrm{h}}{\mathit{y}}_{\mathrm{h}}+{\mathit{T}}_{\mathrm{h}}^{-1}{\mathit{b}}_{\mathrm{h}}={\mathit{W}}_{\mathrm{h}}{\mathit{y}}_{\mathrm{h}}+{\tilde{\mathit{b}}}_{\mathrm{h}} ,$$

(39)

where

$${\mathit{W}}_{\mathrm{h}}=\left[\begin{array}{cc}{\alpha }_{1\mathrm{h}}\left(v\right)& -1\\ 1& {\alpha }_{1\mathrm{h}}\left(v\right)\end{array}\right] \mathrm{and} {\tilde{\mathit{b}}}_{\mathrm{h}}=\left[\begin{array}{c}-\left(\frac{{\mu }_d-{\mu }_s}{2v_0^2}{\mathrm{e}}^{-\frac{v}{v_0}}\right)y_2^2-\left(\frac{{\mu }_d-{\mu }_s}{6v_0^3}{\mathrm{e}}^{-\frac{v}{v_0}}\right)y_2^3\\ 0\end{array}\right].$$

(40)

For $v=v_{\mathrm{h},\mathrm{cr}}$, ${\alpha }_{1\mathrm{h}}=0$, ${\alpha }_{1\mathrm{h}}$ is kept as a generic function of v, since ${\alpha }_{1\mathrm{h}}$ is the critical term causing the instability. The system in Equation (39) is in the so-called Jordan normal form.

By performing several transformations, namely transformation in complex form, near-identity transformation and transformation in polar coordinates, the bifurcation can be characterized through its normal form

$$r^{\prime }={\alpha }_{1\mathrm{h}}\left(v\right)r+{\delta }_{\mathrm{h}}{r}^3 ,$$

(41)

where

$${\delta }_{\mathrm{h}}=\frac{{\mathrm{e}}^{-\frac{v}{v_0}}({\mu }_s-{\mu }_d)}{16v_0^3} .$$

(42)

Details of this standard procedure can be found in [26]. Non-zero real equilibrium solutions of Equation (41) correspond to periodic motion of the system in Equation (5). Linearizing ${\alpha }_{1\mathrm{h}}\left(v\right)$ around $v=v_{\mathrm{h},\mathrm{cr}}$, we obtain

$$r^{\prime }=r({\alpha }_{1\mathrm{h}}^{*}(v-v_{\mathrm{h},\mathrm{cr}})+{\delta }_{\mathrm{h}}{r}^2), \mathrm{where} {\alpha }_{1\mathrm{h}}^{*}=\frac{\mathrm{d}{\alpha }_{1\mathrm{h}}}{\mathrm{d}v}{|}_{v=v_{\mathrm{h},\mathrm{cr}}}=-{\mathrm{e}}^{-\frac{v}{v_0}}\frac{{\mu }_s-{\mu }_d}{2v_0^2} ,$$

(43)

which has solutions

$$r=r_{\mathrm{h}0}=0 \mathrm{and} r=r_{\mathrm{h}}^{*}=\sqrt{-\frac{{\alpha }_{1\mathrm{h}}^{*}(v-v_{\mathrm{h},\mathrm{cr}})}{{\delta }_{\mathrm{h}}}}=2\sqrt{2v_0(v-v_{\mathrm{h},\mathrm{cr}})} .$$

(44)

The trivial solution $r_{\mathrm{h}0}$ exists for any value v and it is stable for $v>v_{\mathrm{h},\mathrm{cr}}$. Differently, $r_{\mathrm{h}}^{*}$ is real only if the argument of the square root in Equation (44) is non-negative, which occurs for $v>v_{\mathrm{h},\mathrm{cr}}$. Since ${\mu }_{\mathrm{s}}>{\mu }_{\mathrm{d}}$, in all relevant cases $\delta$ is positive (see Equation (42)), which, as clearly explained in [26], means that the bifurcation is subcritical. This implies that $r_{\mathrm{h}}^{*}$ corresponds to unstable solutions of Equation (41). This result is in accordance with [1].

A practical consequence of the subcritical character of the bifurcation is that the system, even within the stable region of the equilibrium ($v>v_{\mathrm{h},\mathrm{cr}}$), can experience large oscillations. If the system, while in equilibrium, is subject to a sufficiently large perturbation, which makes it cross the unstable periodic solution in the phase space, it will leave its basin of attraction and it will reach another attractor, which in this case consists of stick–slip oscillations.

The bifurcation diagram in Figure 8a clearly illustrates this scenario. The dashed line indicates a branch of unstable periodic solutions generated at the bifurcation (this branch was obtained through time reverse numerical simulations). The solid line, instead, marks the branch of stick–slip oscillation. The thin solid red line represents the branch of unstable periodic solutions obtained from the analytical computation. We remark on the excellent agreement of the analytically computed solution with the numerical one at low amplitudes. For $v\in (1.15, 1.83)$, the system presents two stable solutions, the trivial one and a stick–slip periodic solution, and an unstable periodic solution, as illustrated in Figure 8b for $v=1.3$. Depending on the initial conditions, the system will either converge towards the trivial solution (red curve in Figure 8c) or will undergo stick–slip oscillations (blue curve in Figure 8c). Numerical solutions were computed utilizing the switch model proposed in [4].

5. Bifurcation Analysis of the Host System with the DVA

To evaluate the DVA performance when stability is lost, we analyze the bifurcation behavior of the system with an attached DVA. The analysis is performed assuming that $\gamma$ and ${\zeta }_2$ are tuned approximately according to Equation (25). An analysis of the eigenvalues of matrix $\mathit{A}$ illustrates that at the loss of stability, if $\gamma$ and ${\zeta }_2$ are tuned approximately according to Equation (25), a couple of complex conjugate eigenvalues leaves the left-hand side of the complex plane, meaning that their real parts become positive. This scenario corresponds to the occurrence of a Hopf bifurcation. We also notice that, if ${\zeta }_2\le {\zeta }_{2\mathrm{opt}}$ and $\gamma ={\gamma }_{\mathrm{opt}}$, not one, but two couples of complex conjugate eigenvalues leave the left-hand side of the complex plane. Referring to the stability chart in Figure 4a, the entire unstable region matrix $\mathit{A}$ has only one couple of eigenvalues with positive real part, except in the loop delimited by Points C and P, where all four eigenvalues have positive real part. This scenario corresponds to a Hopf–Hopf (or double Hopf) bifurcation. In the following, the case of a single Hopf bifurcation is analyzed.

The first step of the analysis consists of transforming the system in Equation (9) into first-order form, similar to Equation (14), but including nonlinear terms up to the third order, which leads to

$$\begin{array}{cc}\hfill & \left[\begin{array}{c}{z}_1^{\prime }\\ z_2^{\prime }\\ z_3^{\prime }\\ z_4^{\prime }\end{array}\right]=\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ -1& -{\gamma }^2\epsilon & 2\psi & -2\epsilon {\zeta }_2\gamma \\ -1& -{\gamma }^2(\epsilon +1)& 2\psi & -2(\epsilon +1)\gamma {\zeta }_2\end{array}\right]\left[\begin{array}{c}{z}_1\\ z_2\\ z_3\\ z_4\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ -\frac{z_3^2({\mu }_{\mathrm{d}}-{\mu }_{\mathrm{s}})e^{-\frac{v}{v_0}}(3v_0+z_3)}{6v_0^3}-\epsilon k_{\mathrm{nl}2}{z}_2^3\\ -\frac{z_3^2({\mu }_{\mathrm{d}}-{\mu }_{\mathrm{s}})e^{-\frac{v}{v_0}}(3v_0+z_3)}{6v_0^3}-(\epsilon +1)\epsilon k_{\mathrm{nl}2}{z}_2^3\end{array}\right]\hfill \\ \hfill & =\mathit{A}\mathit{z}+\mathit{b} .\hfill \end{array}$$

(45)

In the vicinity of the loss of stability, $\mathit{A}$ has two couples of complex conjugate eigenvalues ${\lambda }_{1,2}={\alpha }_1\pm \mathrm{i}{\omega }_1$ and ${\lambda }_{3,4}={\alpha }_2\pm \mathrm{i}{\omega }_2$. To decouple the linear part of the system, we define the transformation matrix

$$\mathit{T}=\left[\begin{array}{cccc}\mathrm{Re}\left({\mathit{s}}_1\right)& \mathrm{Im}\left({\mathit{s}}_1\right)& \mathrm{Re}\left({\mathit{s}}_3\right)& \mathrm{Im}\left({\mathit{s}}_3\right)\end{array}\right] ,$$

(46)

where ${\mathit{s}}_1$ and ${\mathit{s}}_3$ are two of the eigenvectors of $\mathit{A}$, and we apply the transformation $\mathit{z}=\mathit{Ty}$, obtaining

$$\dot{\mathit{y}}={\mathit{T}}^{-1}\mathit{Ay}+{\mathit{T}}^{-1}\mathit{b}=\mathit{Wy}+\tilde{\mathit{b}} ,$$

(47)

where

$$\mathit{W}=\left[\begin{array}{cccc}{\alpha }_1& {\omega }_1& 0& 0\\ -{\omega }_1& {\alpha }_1& 0& 0\\ 0& 0& {\alpha }_2& {\omega }_2\\ 0& 0& -{\omega }_2& {\alpha }_2\end{array}\right] .$$

(48)

For the sake of brevity, the explicit formulation of $\tilde{\mathit{b}}$ is omitted here.

In the case of a single Hopf bifurcation, only ${\alpha }_1$ becomes positive at the loss of stability, while ${\alpha }_2$ remains negative. Therefore, only the first two equations of Equation (47) are linearly related to the bifurcation, while $y_3$ and $y_4$ have minor local effect at the bifurcation. Next, we aim at reducing the dynamics of the system to the so-called center manifold, which is a two-dimensional surface tangent at the bifurcation point to the subspace spanned by the two eigenvectors ${\mathit{s}}_1$ and ${\mathit{s}}_2$ related to the bifurcation. To do so, we approximate $y_3$ and $y_4$ by $y_3={\eta }_{320}{y}_1^2+{\eta }_{311}{y}_1{y}_2+{\eta }_{302}{y}_2^2$ and $y_4={\eta }_{420}{y}_1^2+{\eta }_{411}{y}_1{y}_2+{\eta }_{402}{y}_2^2$, reducing the system to

$$\begin{array}{cc}\hfill & y_1={\alpha }_1{y}_1+{\omega }_1{y}_2+\sum _{j+k=2,3}{a}_{jk}{y}_1^j{y}_2^k+\mathrm{h}.\mathrm{o}.\mathrm{t}.\hfill \\ \hfill & y_2=-{\omega }_1{y}_1+{\alpha }_1{y}_2+\sum _{j+k=2,3}{b}_{jk}{y}_1^j{y}_2^k+\mathrm{h}.\mathrm{o}.\mathrm{t}. ,\hfill \end{array}$$

(49)

where j and k are non-negative integers (more details on this procedure can be found in [26]) and h.o.t. stands for higher order terms.

The system in Equation (49) has the same form as Equation (39); therefore, exactly the same steps can be performed to reduce the system to its normal form, that is

$$r^{\prime }={\alpha }_1\left(v\right)r+\delta r^3 ,$$

(50)

where [26]

$$\delta =\frac{1}{8}(\frac{1}{{\omega }_1}((a_{20}+a_{02})(-a_{11}+b_{20}-b_{02})+(b_{20}+b_{02})(a_{20}-a_{02}+b_{11}))+(3a_{30}+a_{12}+b_{21}+3b_{03})).$$

(51)

Imposing $\epsilon =0.05$, $\gamma ={\gamma }_{\mathrm{opt}}$ and ${\zeta }_2=1.05 {\zeta }_{2\mathrm{opt}}$, we obtain

$$\delta =0.00474+0.68 k_{\mathrm{nl}2} .$$

(52)

Proceeding as done for the host system without DVA, we have that the non-trivial solutions of Equation (50) is given by

$$r=r^{*}=\sqrt{-\frac{{\alpha }_1^{*}(v-v_{\mathrm{cr}})}{\delta }} ,$$

(53)

where ${\alpha }_1^{*}=\mathrm{d}{\alpha }_1{/\mathrm{d}v|}_v=v_{\mathrm{cr}}$. We notice that $\delta$ is positive if the DVA is linear ($k_{\mathrm{nl}2}=0$), which means that also in this case the bifurcation is subcritical, and it generates unstable periodic solutions. Analyzing other values of ${\zeta }_2$ and $\epsilon$, we verified that the subcritical characteristic persists for a relatively large parameter value range. The corresponding bifurcation diagram is illustrated in Figure 9a. Comparing Figure 9a and Figure 8a, we notice that, although the linear DVA does not change the characteristic of the bifurcation, the advantages in terms of vibrations suppression persist also in the nonlinear range. In fact, for the considered parameter values, in the host system without DVA, stick–slip oscillations exist for $v\in \left(0, 1.83\right]$, while, with the addition of the absorber, they are limited to the range $v\in \left(0, 0.768\right]$.

Equation (52) suggests that, if $k_{\mathrm{nl}2}<-0.00697$, $\delta$ becomes negative, making therefore the bifurcation supercritical. This scenario is confirmed by the bifurcation diagram depicted in Figure 9b, for $k_{\mathrm{nl}2}=-0.01$. Although at first sight it seems that the bifurcation is subcritical, the inset illustrates that the bifurcation is indeed supercritical; however, the branch of periodic solutions bends rapidly to the right in correspondence of a fold bifurcation, making the overall scenario similar to the case of $k_{\mathrm{nl}2}=0$. The figure confirms the correctness of the analytical computation; nevertheless, it also points out that the performed local analysis is unable to capture the global behavior of the system, which is not qualitatively affected by the variation of the nonlinear characteristic of the DVA’s spring. Furthermore, we notice that the addition of the softening nonlinear spring enlarges the bistable range, making stick–slip oscillations exist up to $v<0.86$, instead of 0.768 as in the case of $k_{\mathrm{nl}2}=0$.

Figure 9c illustrates the bifurcation diagram obtained for a hardening absorber’s spring ($k_{\mathrm{nl}2}=0.01$). In this case, the range of existence of stick–slip oscillations is further enlarged, persisting up to $v<1.035$. We also remark that increasing the value of $k_{\mathrm{nl}2}$ above 0.01 or decreasing it below −0.01 provided only worse performance than those illustrated in Figure 9. This result suggests that any low order nonlinearity of the absorber’s stiffness is detrimental concerning the DVA effectiveness. This finding is somehow surprising, considering that in similar applications the addition of a properly tuned nonlinear term in the DVA’s stiffness provided some advantages [20,25,27].

Regarding Figure 9b, we notice that the branch of stick–slip oscillations presents two folds for $v\approx 0.48$. However, an analysis of the system’s steady state solutions before and after the folds did not reveal any particular detail relevant from an engineering point of view; therefore, the phenomenon was not analyzed in further detail. We also remark that, in Figure 9b,c, the branches of stable and unstable solutions do not encounter each other at a well defined point, as happens in Figure 8a, for instance. This is probably related to the fact that the branches of unstable solutions in Figure 9 were obtained adopting the shooting method (employing MatCont [28], a MATLAB-based toolbox for numerical continuation) of the system smoothed assuming that $v_{\mathrm{rel}}$ is always positive. This assumption makes the considered system unable to exhibit stick–slip oscillations, but keeps it equivalent to the original system for $v>z_3$. In contrast, the stable branches were obtained from direct numerical simulations of the full system. Therefore, inaccuracies of the smoothed system in the proximity of the onset of stick–slip motions are possible.

As mentioned at the beginning of this section, for ${\zeta }_2<{\zeta }_{2\mathrm{opt}}$ and $\gamma ={\gamma }_{\mathrm{opt}}$, the system undergoes a Hopf–Hopf bifurcation. However, acknowledging the fact that the bifurcation analysis seems to be an inefficient tool for investigating the post-bifurcation behavior of the system, which is dominated by large amplitude oscillations, and considering that the analysis of such a bifurcation requires a significant analytical effort, the detailed investigation of this case is omitted in this study.

6. Conclusions

In this study, the problem of suppressing FIVs through a DVA was addressed. Possibly the simplest system exhibiting FIVs was considered, i.e., the mass-on-moving-belt system, to which a classical DVA was attached. The optimal tuning of the absorber parameters was defined through an analytical procedure, which enabled us to reduce the critical velocity by approximately 50%, with an additional mass of only 5% of the primary system’s mass.

The post-bifurcation behavior analysis illustrated that, although a linear DVA is unable to change the bifurcation character at the loss of stability, it can still significantly reduce the extent of the bistable region. Globally, the area of existence of stick–slip oscillations is reduced by 58%, with a DVA mass of only 5%. The bifurcation analysis proved that it is possible to change the bifurcation character if a small softening term is included in the absorber. However, this has only a local beneficial effect, while, globally, it enlarges the region of existence of stick–slip motions. The performance also worsens if an additional hardening term is introduced, suggesting that the spring characteristic should be maintained as linear as possible. Large order nonlinearities, such as non-smoothness, might have beneficial effects; nevertheless, their analysis was not addressed in this study, and it is left for future developments. Other possible future developments of the present study include the analysis of the Hopf–Hopf bifurcation occurring at the loss of stability for ${\zeta }_2<{\zeta }_{2\mathrm{opt}}$ and the analysis of the performance of the DVA if the primary system has two DoF, encompassing, therefore, coupling instabilities as well [22].