## Abstract

The inerter has been integrated into various vibration mitigation devices, whose mass amplification effect could enhance the suppression capabilities of these devices. In the current study, the inerter is integrated with a pendulum vibration absorber, referred to as inerter pendulum vibration absorber (IPVA). To demonstrate its efficacy, the IPVA is integrated with a linear, harmonically forced oscillator seeking vibration mitigation. A theoretical investigation is conducted to understand the nonlinear response of the IPVA. It is shown that the IPVA operates based on a nonlinear energy transfer phenomenon wherein the energy of the linear oscillator transfers to the pendulum vibration absorber as a result of parametric resonance of the pendulum. The parametric instability is predicted by the harmonic balance method along with the Floquet theory. A perturbation analysis shows that a pitchfork bifurcation and period doubling bifurcation are necessary and sufficient conditions for the parametric resonance to occur. An arc-length continuation scheme is used to predict the boundary of parametric instability in the parameter space and verify the perturbation analysis. The effects of various system parameters on the parametric instability are examined. Finally, the IPVA is compared with a linear benchmark and an autoparametric vibration absorber and shows more efficacious vibration suppression.

## 1 Introduction

The inerter is a mechanical device with two terminals, each of which exerts an equal and opposite inertial force proportional to the relative acceleration between the terminals [1]. The inerter amplifies the inertial effects of a small mass by using motion transmission mechanisms, fluids, and levers [1]. By virtue of its mass amplification effect, the inerter has been studied to improve the performance of various passive vibration mitigation techniques in the last decade. Ikago et al. [2] developed the tuned viscous mass damper (TVMD), which consisted of a tuning spring in series with the inerter and a viscous damper in parallel. It was shown that the TVMD outperformed the viscous damper alone when applied to a seismically excited single degree-of-freedom (SDOF) structure. Furthermore, Lazar et al. [3] proposed the tuned inerter damper (TID), wherein the inerter was substituted for the oscillating mass of a tuned mass damper (TMD). The TID and TMD were compared in seismically excited multiple-degree-of-freedom structures and demonstrated similar effectiveness. Later, Lazar et al. [4] considered the TID in suppressing the midspan vibration of cables and showed that the TID outperformed the optimal viscous damper. Moreover, Qian et al. [5] studied serial and parallel connections between the TID and a base-isolation system and concluded that the serial TID outperformed the parallel TID for practical structures.

The inerter has also been applied to enhance the inertial effects of dynamic vibration absorbers (DVAs). Marian and Giaralis [6] proposed the tuned mass damper inerter (TMDI), which consisted of a TMD and the inerter in series. In a 3DOF structure simulation, they showed that for achieving similar vibration control performance, the weight of the TMDI was four times lighter than the TMD. Furthermore, De Domenico and Ricciardi [7] incorporated the TMDI in a base-isolation system and demonstrated that the displacement demand of the base-isolated structure could be significantly reduced. Moreover, Joubaneh and Barry [8] studied the performance of four models of electromagnetic resonant shunt TMDI (ERS-TMDI) on both vibration suppression and energy harvesting and identified the best model. Their parametric studies showed that increasing the inertance enhances the performance of the best model in terms of both vibration mitigation and energy harvesting. On the other hand, Tai [9] proposed the tuned inerter-torsional-mass damper (TITMD), which integrated the inerter and a torsional mass damper. In comparison with the TMDI, the TITMD achieved 20–70% improvement when having identical weights.

The rest of this article is organized as follows. In Sec. 2, we will present the design of the IPVA and derive the equations of motion. In Sec. 3, we study the stability of periodic solutions of the system and conduct a bifurcation analysis to determine the boundary of parametric instability. In Sec. 4, direct numerical integration is used to verify the stability boundary obtained from Sec. 3. Section 5 shows the effects of various parameters on the stability boundary. Section 6 compares the proposed system with a linear benchmark and an autoparametric vibration absorber in Ref. [17]. This study is concluded in Sec. 7.

## 2 Inerter Pendulum Vibration Absorber

In this section, we present the design of the IPVA and derive the corresponding equations of motion. We present two designs: rack-pinion based IPVA and ball screw based IPVA, as shown in Fig. 1. Specifically, the rack-pinion and ball screw IPVA are used to suppress vibrations in the horizontal and vertical directions, respectively. Their working principles are briefly described as follows. Figure 1(a) shows a SDOF linear oscillator (primary structure) that moves in the horizontal direction and is excited by an external force F, which is modeled by a linear spring of stiffness k, a mass M, and a viscous damper of damping coefficient c. The external force is assumed to be harmonic, i.e., F = F0 sin Ωt, where F0 and Ω are the force amplitude and excitation frequency, respectively. Denote by x the displacement of primary structure. To suppress the vibration of primary structure, a rack-pinion based pendulum vibration absorber is considered, which consists of a rack-pinion of radius R, a carrier of radius Rp, and a pendulum vibration absorber of mass m and length r. The rack-pinion is installed between the primary mass and the ground to convert the linear oscillation of the primary mass into rotation of the pinion. As a result, the linear displacement of primary structure and the pinion’s angular displacement are related through x = . Note that because the pendulum oscillates in the horizontal plane, the gravity is neglected. Moreover, the carrier is fixed to the pinion, and the pendulum is pivoted on a point on the circumference of the carrier. Figure 1(c) shows a primary structure that moves in the vertical direction and incorporates the ball screw IPVA. The ball screw is attached between the ground and the primary structure such that linear motion x of the mass is converted to angular motion θ of the ball screw with effective radius R = L/2π. L refers to the lead value associated with the ball screw as a measure of the ratio of linear displacement of the nut to full rotation of the screw. More specifically, x = , which is the same as the rack-pinion design. A pendulum, of length r and mass m, is also attached perpendicular to the screw, at a radius of Rp from the center of the screw and with an angular displacement of ϕ with respect to the attachment point. Note that two pendulums are shown, but one is to be fixed in place to avoid rotating unbalance. Also, the pendulum moves in the horizontal plane; thus, the gravity is also neglected in the ball screw IPVA.

Fig. 1
Fig. 1
Close modal

Both systems consist of two degrees-of-freedom, one associated with the pinion’s angular displacement (θ) and the other with the pendulum’s angular displacement relative to the pinion (ϕ).

### 2.1 Equations of Motion.

Although different mechanisms are used, the working principle of both systems are identical. Therefore, their equations of motion are identical. Lagrange’s equations are used to derive the equations of motion. First, the total kinetic energy of the system is derived as follows:
$T=TM+Tc+Tp$
(1)
where
$TM=12M(Rθ˙)2,Tc=12Jθ˙2,Tp=12Jp(θ˙+ϕ˙)2+12m(Rp2θ˙2+r2(θ˙+ϕ˙)2+2Rprcos(ϕ)θ˙(θ˙+ϕ˙))$
(2)
are the kinetic energy of the structure, carrier, and pendulum, respectively. Here, J is the moment of inertia of the carrier-pinion composite and Jp is the moment of inertia of pendulum with respect to its center of mass
$V=12kx2=12kR2θ2$
(3)
To account for energy loss at the pivot point of the pendulum, a torsional viscous damping coefficient cp is introduced. The virtual work done by the force F = F0 sin(Ωt), the damping torque in the pendulum, and the damping force in the primary mass can be derived as Fδx, $−cpϕ˙δϕ˙,$ and $−cx˙δx˙$, respectively, where cp and c are the torque damping coefficient in the pendulum and damping coefficient of the viscous damper between primary mass and ground, respectively. Then, the virtual work done by the force F, the damping torque (due to cp and viscous damping c) are derived as follows:
$δW=FRδθ−cpϕ˙δϕ−cR2θ˙δθ$
(4)
Therefore, the equations of motion of the system obtained using the Lagrange’s equations are written as follows:
$(MR2+J+mRp2+mr2+2mRprcosϕ)θ¨+(mr2+mRprcosϕ+Jp)ϕ¨+cR2θ˙+kR2θ−2mRprϕ˙θ˙sinϕ−mRprϕ˙2sinϕ=F0sin(Ωt)R,(mr2+Jp)ϕ¨+m(r2+Rprcosϕ)θ¨+cpϕ˙+mRprθ˙2sinϕ=0$
(5)
It is assumed that the pendulum is made of a point mass such that its moment of inertia with respect to the pivot point is much larger than the pendulum’s moment of inertia with respect to its center of mass, i.e., mr2 > > Jp. Furthermore, $J=(mp+mc)Rg2$, where mp and mc are the pinion mass and carrier mass, respectively, and Rg is the radius of gyration. As the primary mass M is much larger than the sum of mp and mc, and as the pinion radius R and the radius of gyration Rg have the same order of magnitude, it is assumed that MR2 > >J. Without loss of generality, J and Jp are neglected. We rescale the time and convert Eq. (5) into a dimensionless form for further analysis using the following parameters:
$μr=mRp2MR2,ω0=kM,ω=Ωω0,τ=ω0t,η=rRpξ=c2ω0M,ξp=cp2ω0MR2,f=F0MRω02,()′=d()dτ$
(6)
Denote x = [θ, ϕ]T and f = [f sin ωτ, 0]T. The dimensionless equations of motion are obtained as follows:
$Mx″+Cx′+Kx+g(x,x′,x″)=f$
(7)
where
$M=[1+μr(1+η2)μrη2μrη2μrη2],C=[2ξ002ξp],K=[1000],g(x,x′,x″)=μrη[(2θ″+ϕ″)cosϕ−ϕ′(2θ′+ϕ′)sinϕθ″cosϕ+θ′2sinϕ]$
(8)

It is worth noting that the strength of the nonlinear inertial terms g(x, x′, x″) is proportional to μr and η. The moment of inertia ratio μr can be readily magnified by adjusting the ratio Rp/R, thereby creating strong nonlinear inertial effects with a small pendulum mass. For example, for a mass ratio $mM=3%$, a ratio $Rp/R=10$ leads to μr = 0.3, indicating that the inertia effect is magnified by a factor of ten. Furthermore, the pendulum length ratio η is proportional to the length of the pendulum. Therefore, a long pendulum leads to strong nonlinear inertial effects.

## 3 Parametric Resonance of IPVA

According to the studies on autoparametric resonance, parametric resonance plays an essential role in transferring the kinetic energy of a primary structure to the pendulum vibration absorber [1618]. As will be demonstrated in Sec. 6, when parametric resonance occurs to the IPVA, a similar energy transfer phenomenon is observed, resulting in vibration mitigation of the primary structure. In this section, we will determine the conditions for which parametric resonance will occur to the IPVA. To this end, we use the harmonic balance method to determine the parametric instability of the system.

### 3.1 Harmonic Balance Method.

By virtue of the harmonic balance method, periodic solutions of the system are assumed to take the following form:
$θp(τ)=∑p=1P(Θpccos(pωτν)+Θpssin(pωτν))ϕp(τ)=Φ0+∑p=1P(Φpccos(pωτν)+Φpssin(pωτν))$
(9)
where Θp, Φp, and Φ0 are unknown Fourier coefficients to be determined. Note that $ν∈N$ accounts for subharmonics. Furthermore, a constant Φ0 is included to consider asymmetric oscillation of the centrifugal pendulum [24]. Denote by xp = [θp, ϕp]T the vector of assumed periodic solutions. After substituting Eq. (9) into the equations of motion (7), we obtain the following residue term:
$R(τ)=Mxp″+Cxp′+Kxp−g(xp,xp′,xp″)−f$
(10)
To obtain an expression relating the Fourier coefficients, a Galerkin procedure [25] is used to project (10) on the orthogonal trigonometric basis, yielding 2P + 1 nonlinear algebraic equations
$h0(x^)=∫02πν/ωR(τ)dτ=0,hps(x^)=∫02πν/ωR(τ)sin(pωτν)dτ=0hpc(x^)=∫02πν/ωR(τ)cos(pωτν)dτ=0$
(11)
where , $Θ=[Θ1c,…,ΘPc,Θ1s,…,ΘPs]T$, and $Φ=[Φ1c,…,ΦPc,Φ1s,…,ΦPs]T$. Note that $g(xp,xp′,xp″)$ will result in composite trigonometric terms such as $cos(Φpssin(pωτ/ν))$. These terms can be expanded using the Jacobi–Anger expansion, namely, an infinite series of products of Bessel functions and trigonometric functions [26] (see also Appendix  B for the expansion formulas). For the current study, the Jacobi–Anger expansion is truncated at Bessel functions of order up to third to capture the necessary nonlinear effects.

We solve Eq. (11) for the Fourier coefficients using the Newton–Raphson method. Substitution of the Fourier coefficients into Eq. (9) will lead to the periodic solutions. The stability of the periodic solutions will be determined in the next section.

### 3.2 Stability.

To determine the stability of the periodic solutions, small perturbations are introduced into Eq. (9) as follows:
$θ(τ)=θp(τ)+δθ(τ)andϕ(τ)=ϕp(τ)+δϕ(τ)$
(12)
where $|δθ(τ)|<<1$ and $|δϕ(τ)|<<1$. Denote by $δ=[δθ,δϕ]T$ the vector of small perturbations. Substitution of Eq. (12) into Eq. (7) and linearization with respect to θp(τ) and ϕp(τ) yield
$(M+∂g∂x″)δ″+(C+∂g∂x′)δ′+(K+∂g∂x)δ=0$
(13)
where the Jacobian matrices ∂g/∂x″, ∂g/∂x′, and ∂g/∂x are evaluated at x = xp, $x′=xp′$, and $x″=xp″$, respectively, wherever appropriate. Note that the Jacobian matrices are periodic functions of period T = 2πν/ω, e.g., ∂g/∂x(τ) = ∂g/∂x(τ + T). Because Eq. (13) have periodic coefficients, one can use Floquet theory to determine the stability [27]. To this end, Eq. (13) is transformed into the state-space form and numerically integrated using matlab’s ode45 (based on an explicit Runge–Kutta integration method) over one period T to obtain the fundamental matrix. The absolute and relative tolerance were taken to be 10−9 to ensure accuracy. If any eigenvalue of the fundamental matrix has a magnitude greater than unity, the periodic solutions are unstable.
Although Eq. (13) can be used to determine the stability of arbitrary periodic solutions, it is a numerical approach; hence, it is hard to understand how parametric resonance occurs to the IPVA. To gain physical insights, we also use a semi-analytical approach to determine the stability. To that end, we apply a multiple-scale approach to Eq. (13) as follows. Because a compact and lightweight design of the IPVA system is preferred in practical applications, we consider μr < < 1. Assuming that the parameters μr, ξ, ξp are small quantities, we set $μr=ϵμ^r,ξ=ϵξ^,ξp=ϵξ^p$ and introduce the following asymptotic expansions:
$δθ,ϕ(τ)=δθ,ϕ(0)(τ0,τ1,…)+ϵδθ,ϕ(1)(τ0,τ1,…)+⋯τk=ϵkτ,k=0,1,…ddτ=∂∂τ0+ϵ∂∂τ1+⋯$
(14)
where |ε| < <1 is a small bookkeeping parameter.
After substituting Eq. (14) into Eq. (13) and collecting terms that will lead to parametric instabilities, the equation obtained in order O(ε0) is expressed as follows:
$∂2δθ(0)∂τ02+δθ(0)=0,∂2δϕ(0)∂τ02+2ξ^pμ^rη2∂δϕ(0)∂τ0+A(xp)ηδϕ(0)=0$
(15)
where
$A(xp)=cos(ϕp)(θp′)2−sin(ϕp)θp″$
(16)
is a periodic coefficient of period T. It is worth noting that the first equation in Eq. (15) shows that $δθ(0)$ are stable harmonic functions. Therefore, the stability of the periodic solutions are determined by the second equation in Eq. (15).
When the nonlinearity is weak, periodic solutions (9) are dominated by primary harmonics, i.e., P = 1 and ν = 1. As Eq. (15) is derived by the multiple-scale approach, it is accurate when the nonlinearity is weak. Therefore, in addition to Eq. (13) (Floquet theory), Eq. (15) is used to determine the boundary of parametric instability for periodic solutions of primary harmonics, which will explain how parametric resonance occurs to the IPVA. Thus, P = 1 and ν = 1 are substituted in Eq. (9) to obtain $θp=Θ1ccos(ωτ)+Θ1ssin(ωτ)$ and $ϕp=Φ0+Φ1ccos(ωτ)+Φ1ssin(ωτ)$. After substituting these expressions of Θ and Φ into Eq. (15) and expanding in terms of Bessel functions up to third order, we arrive at a damped Mathieu equation as follows:
$∂2δϕ(0)∂τ02+2ξpμrη2∂δϕ(0)∂τ0+u(x^)ηδϕ(0)+v(x^)ηcos(ωτ−γ)δϕ(0)=0$
(17)
where
$u(x^)=ω2Θ1c(Φ0)2{Θ1[J0(Φ1)+J2(Φ1)s(2α)]+2J1(Φ1)s(α)},v(x^)2=[b0b1+2b2(b1+b3)]c(α)+(2b0b2+2b1b3−b12)c(2α)+(2b0b3−b1b2)c(3α)−b1b3c(4α)+b02+5b124+b22+b32,b0=s(Φ0)Θ1ω2J0(Φ1),b1=−s(Φ0)Θ12ω2J1(Φ1),b2=−s(Φ0)Θ1ω2J2(Φ1),b3=−s(Φ0)Θ12ω2J3(Φ1)Θ1=(Θ1c)2+(Θ1s)2,Φ1=(Φ1c)2+(Φ1s)2,$
(18)
$α=tan−1(Φ1sΘ1c−Φ1cΘ1sΦ1cΘ1c+Φ1sΘ1s)$
(19)
where c(·) = cos(·), s(·) = sin(·), and Jn(·) denotes Bessel functions of the first kind of order n. Note that the detail of phase angle γ is not provided because γ is irrelevant to stability. There are two things worth noting in Eqs. (17) and (19). First, b0 = b1 = b2 = b3 = 0 or $v(x^)=0$ when Φ0 = 0. Because $v(x^)$ is the magnitude of parametric excitation, no parametric instabilities can occur when $v(x^)=0$. In other words, nonzero asymmetric oscillation, i.e., Φ0 ≠ 0, is a necessary condition for parametric instabilities. Second, the linear stiffness term $u(x^)$ is composed of nonlinear inertial coupling induced by the carrier motion Θ1. As the nonlinear inertial coupling results in linear stiffness per se, the pendulum can have parametric resonance without having any linear stiffness. Compared to the autoparametric vibration absorbers, which would need to have low linear stiffness to tune their natural frequency around half the natural frequency of the primary structure [1618], the nonlinear inertial coupling of the IPVA enables compact designs.
To determine the boundary of parametric instability, Eq. (17) is transformed to the standard form of Mathieu equation [13]:
$∂2Ψ∂w2+p(x^)Ψ−2q(x^)cos(2w)Ψ=0$
(20)
where
$Ψ=δϕ(0)exp(2ξ^pwμ^rη2ω),2w=ωτ−γp(x^)=4u(x^)ηω2−4ξ^p2μ^r2η4ω2,q(x^)=−2v(x^)ηω2$
(21)
The boundary of parametric instability for Eq. (20) corresponds to the transition curves in the pq plane [28]. Because we seek the boundary that occurs with low force magnitudes, we compute the transition curve that occurs with the lowest p and q values. Mathematically, this transition curve is expressed as $p=A1(q)$, where $A1(q)$ are the characteristic values for even Mathieu functions with characteristic exponent 1 and parameter q [28]. In this article, $A1(q)$ is computed by the “MathieuCharacteristicA” function of wolfram mathematica 11.3. Note that p and q are functions of f and ω. Therefore, $p=A1(q)$ is solved with Eq. (11) simultaneously to yield the transition curves in the fω plane.

### 3.3 Pitchfork Bifurcation.

As mentioned in Sec. 3.2, nonzero asymmetric oscillation, i.e., Φ0 ≠ 0, is necessary to induce parametric instabilities. To determine when it occurs, $θp=Θ1ccos(ωτ)+Θ1ssin(ωτ)$ and $ϕp=Φ0+Φ1ccos(ωτ)+Φ1ssin(ωτ)$ are substituted into Eq. (13) to solve for stable periodic solutions with Φ0 ≠ 0. We use the pendulum length ratio η as the bifurcation parameter to obtain a bifurcation diagram that shows the parameter space, wherein Φ0 ≠ 0 will occur. To track the bifurcation points with varying η, a bifurcation tracking algorithm that is based on arclength continuation is used with Eq. (13); see Appendix  A for the detail. Figure 2 shows a bifurcation diagram of Φ0 with varying η. Three branches of the bifurcation were obtained using three different sets of initial conditions (one corresponding to each branch, namely, the lower, middle, and upper). It can be observed that Φ0 undergoes a supercritical pitchfork bifurcation at a critical value of η. After this critical value of η, Φ0 ≠ 0 and parametric instabilities become possible. For the rest of this article, we will only explore parametric instabilities with the parameters that lead to Φ0 ≠ 0.

Fig. 2
Fig. 2
Close modal

### 3.4 Period Doubling Bifurcation.

Within the parameter space wherein Φ0 ≠ 0 exist, the boundary of parametric instability is computed in the fω plane. To find an initial bifurcation point for the bifurcation tracking algorithm described in Appendix  A, ω = 0.8 is set and Eq. (13) is repeatedly used to compute the Floquet multipliers as f decreases until the maximum magnitude of the Floquet multipliers becomes unity. Afterward, the bifurcation tracking algorithm will generate the boundary as described in Appendix  A. To verify whether the boundary is indeed of parametric instability, the Mathieu equation (20) is used to generate the transition curve as described in Sec. 3.2. The boundary and transition curve for a set of parameters are shown in Fig. 3. As shown, although the transition curve underestimates the boundary, they are in qualitative agreement. Specifically, the discrepancy between the two curves increases as the force magnitude f increases. Because the transition curve is predicted by the perturbation method, it is reasonable that it is more accurate for small force magnitudes. Thus, the comparison verifies the claim that the boundary indicates parametric instability. To gain more insight, the Floquet exponents corresponding to a few points on the boundary are computed and found equal to ±/T, where $i=−1$. According to Ref. [29], this indicates period doubling bifurcation. Since periodic doubling bifurcation is a co-dimension one bifurcation, it is a curve in a parameter plane [30]. Therefore, the parametric instability boundary is in fact a boundary of period doubling bifurcation. When this bifurcation occurs, the pendulum oscillation will have subharmonics of ω/2, i.e., ν = 2 in Eq. (9). It is worth noting that the autoparametric vibration absorbers also have a similar bifurcation behavior, that is, subharmonics of half excitation frequency induced by parametric instabilities [1618].

Fig. 3
Fig. 3
Close modal

Within the parameter space wherein subharmonics of ω/2 exist, the stability of the subharmonics can be further investigated. Preliminary investigations indicates the presence of another period doubling bifurcation, implying that subharmonics of quarter frequency will appear. Therefore, it is hypothesized that there exists a cascade of period doubling bifurcations in the fω space, which eventually leads to chaotic motions of the system. Determination of the boundary of this additional period doubling bifurcation, however, is out of the scope of this paper.

Figure 3 is a bifurcation diagram that shows the parameter space for qualitatively different solutions, defined by the instability boundary. By locating the parameters in Fig. 3, the qualitative behavior of the corresponding solutions can be predicted. For example, $×2$ resides in the parameter space just above the boundary. Accordingly, periodic solutions of primary harmonics along with subharmonics of excitation frequency are predicted at $×2$. Next, we verify the predictions by Fig. 3 by direct numerical integration.

## 4 Numerical Demonstration

To verify the bifurcation analysis in Sec. 3.4, numerical integration (matlab’s ode45) is used to obtain the solutions of Eq. (7) at three representative points in Fig. 3 (denoted by markers “$×$” followed by numbers, e.g., $×2$). Among these three points, points $×1$ and $×2$ lead to periodic solutions, whereas point $×3$ leads to nonperiodic solutions. The fast Fourier transform (FFT) of the periodic solutions are computed to reveal the frequency components, which are shown in Figs. 4(a), 5(a), and 6(a). On the other hand, a time series of the solutions are presented to show the dynamical behaviors, shown in Figs. 4(a), 5(b), 6(b), 4(c), 5(c), and 6(c). Note that the frequencies $ω^$ of the FFT are normalized with respect to the excitation frequency. It follows that primary harmonics correspond to components at $ω^=1$, subharmonics of half excitation frequency correspond to components at $ω^=0.5$, etc.

Fig. 4
Fig. 4
Close modal
Fig. 5
Fig. 5
Close modal
Fig. 6
Fig. 6
Close modal

There are several things worth noting in Figs. 36. First, the prediction at point $×1$ is in good agreement with the numerical solutions. As shown in Fig. 3, point $×1$ is below the instability boundary. It is expected that the periodic solutions are dominated by primary harmonics. This prediction is verified by Fig. 4, which shows that the periodic solutions at point $×1$ have the largest components at $ω^=1$, corresponding to primary harmonics. Furthermore, in Fig. 3, as we increase the value of f and reach point $×2$, the primary harmonics undergo a period doubling bifurcation. As a result, subharmonics of half excitation frequency should arise. As shown in Fig. 5, subharmonics of half excitation frequency indeed exist, which verifies the prediction in Fig. 3. Second, the parameters at $×3$ lead to strong nonperiodic solutions composed of both oscillation and intermittent rotations of the pendulum, as shown in Figs. 6(c) and 7(c). Similar nonperiodic solutions are also observed in autoparametric resonance vibration absorbers [18].

Fig. 7
Fig. 7
Close modal
In addition to FFT, the Poincaré sections are used to demonstrate the period doubling bifurcations predicted by Fig. 3. The Poincaré sections are computed by the Hénon trick [31], which are defined as follows:
$Pn(x0)=xn(τ0+2nπ/ω;xn−1,τ0),n=1,2,…$
(22)
where xn−1 and xn are the solutions of the system (7), which pass through the Poincaré section at time τ = τ0 + 2(n − 1)π/ω and τ = τ0 + 2/ω, respectively. Successively, the points x0, x1 = P(x0), x2 = P2(x0), … correspond to the intersection of the trajectory x(τ; x0, τ0) with the sections at τ = τ0, τ0 + 2π/ω, τ0 + 4π/ω, …, respectively.

To demonstrate the period doubling bifurcation, Poincaré sections corresponding to point $×1$ and $×2$ are plotted in Fig. 7. As shown in Fig. 7(a), $×1$ leads to a fixed point on the Poincaré section, corresponding to period-1 solutions. Figure 7(b), on the other hand, shows two fixed points, corresponding to period-2 solutions. Therefore, it is clear that the system has undergone a period doubling bifurcation when moving from $×1$ to $×2$.

## 5 Parametric Studies

In this section, we analyze the effect of the parameters on the instability boundary. We consider four parameters in Eq. (6), namely, μr, η, ξ, and ξp. It can be seen that these parameters can be varied independently of each other. Therefore, we will observe the effect by varying one parameter while keeping the other parameters constant. We start by increasing the value of η while keeping the others constant. From Fig. 8, it can be observed that increasing η does not make any significant change in the lowest f value for parametric instability to occur, which corresponds to the vertex of the boundaries. That means that value of η should not influence the energy transfer capabilities of the system by a lot. However, a minimum threshold value of η is required for period-doubling bifurcation to occur as discussed in Sec. 3.3. We next vary μr. From Fig. 9, it can be observed that the required values of f to attain parametric instability decrease as μr increases. This can be attributed to the fact that the inertia supplied by the pendulum vibration absorber increases as the mass amplification factor μr is increased. The value of μr can be controlled by changing the ratio Rp/R, which can be adjusted by changing the carrier radius (Rp).

Fig. 8
Fig. 8
Close modal
Fig. 9
Fig. 9
Close modal

While keeping the other parameters constant, we vary the viscous damping ratio ξ and observe its effects. From Fig. 10, it can be seen that the requirement of f to achieve nonlinear energy transfer increases with the increase in the viscous damping. In a similar fashion, we vary the ξp while keeping the other parameters constant. We see that the values of f required to achieve parametric instability increase as the viscous damping increases. The observations on the effects of both viscous damping match well with the effect of viscous damping on parametric instability—the larger the viscous damping, the larger the force it takes to cause parametric instability [28].

Fig. 10
Fig. 10
Close modal

Last but not least, the parameter μr has a significant influence on the instability boundary, as demonstrated in Sec. 5. As shown in Fig. 9, a larger μr not only leads to lower force magnitudes required to cause parametric resonance but a wider frequency bandwidth of parametric resonance, which is beneficial in terms of vibration mitigation. The parameter μr can be readily increased by changing the ratio of Rp/R without incurring the large weight to the system, which is attributed to the mass amplification effect of the inerter.

## 6 Discussion

In the beginning of this study, it was proposed that a nonlinear energy transfer phenomenon similar to autoparametric resonance takes place when the parametric instability occurs. To demonstrate this, we compare the proposed system with two systems, a linear benchmark and an autoparametric vibration absorber with parametrically excited pendulum [17]. The linear system here is characterized by locking the pendulum at its initial position (ϕ = 0), effectively removing all the nonlinearities in the system. By setting ϕ = ϕ′ = ϕ″ = 0 in Eq. (7), the equation of motion of the linear system is written as follows:

$[1+μr(1+η2+2η)]θ″+2ξθ′+2ξpθ′+θ=fsinωτ$
(23)
By using $θ=12Θeiωτ+c.c.$, the equation can be solved to obtain
$Θ=|f2iω(ξ+ξp)+1−ω2[1+μr(1+η2+2η)]|$
(24)
We computed the root-mean-square (rms) of the IPVA system and compared it with the linear system. The comparison is shown in Fig. 12. The response from the 2401th to 3000th cycle was used to compute the rms to eliminate transient effects. The IPVA parameters used in Figs. 12(a) and 12(b) correspond to Figs. 3 and 11, respectively.
Fig. 11
Fig. 11
Close modal
Fig. 12
Fig. 12
Close modal

Several things are worth noting in Fig. 12. First, it is shown that the response of the primary structure flattens for a range of excitation frequencies. In comparison with the response of the linear system, the IPVA shows significant vibration suppression with the flattening region. For example, as shown in Fig. 12(b), θ for f0 = 0.025 flattens for ω ∈ [0.81, 0.87]. In comparison with Figs. 3 and 11, it is clear that the flattening occurs when the system is within the parametric instability boundary. For example, as shown in Fig. 11, when ξp = 0.015 and f0 = 0.025, the system is within the parametric instability boundary for ω ∈ [0.81, 0.87]. This observation agrees with Fig. 12(b). Second, within the flattening region, the response of the primary structure barely increases despite an increase in the force magnitude, suggesting a saturation phenomenon similar to autoparametric vibration absorbers [23] and nonlinear vibration absorbers with quadratic nonlinearities [32]. Note that the response of the IPVA system in Fig. 12(b) is nonperiodic for f = 0.035 within a range of ω; thus, different initial conditions may lead to different rms responses. To examine the effect of initial conditions, ten rms responses were computed and plotted at six discrete ω values, each corresponding to a different initial condition vector $[θ,ϕ,θ˙,ϕ˙]T$ that was randomly chosen from a standard normal distribution with zero mean and a unit standard deviation; see the inset in Fig. 12(b).

Next, we compare IPVA with the autoparametric vibration absorber shown in Fig. 13. Because the autoparametric system oscillates in the vertical direction, the ball screw IPVA is considered hereinafter. The equation of motion of the autoparametric system is written as follows [17]:
$(M+m)x¨+cx˙+kx+ml(ϕ¨sinϕ+ϕ˙2cosϕ)=F0sin(Ωt)ml2ϕ¨+caϕ˙+(mgl+mlx¨)sinϕ=0$
(25)
where x and ϕ represent the primary structure displacement and pendulum angular displacement, respectively, M, k, and c are the mass, stiffness, and viscous damping coefficient of the primary structure, respectively, and m and l are the pendulum mass and length, respectively. Furthermore, a viscous damping coefficient ca is introduced to account for energy loss at the pivot point of the pendulum.
Fig. 13
Fig. 13
Close modal

According to Ref. [17], when the natural frequency of the pendulum is tuned around half of the natural frequency of the primary structure, the system shows autoparametric resonance for a certain set of parameters when excited harmonically. This autoparametric resonance results in energy transfer from the primary structure to the pendulum, thereby achieving vibration suppression of the primary structure. Because the IPVA system and the autoparametric system achieve vibration suppression in a similar way, the latter is an ideal benchmark system for comparison. For a fair comparison, the primary structure parameters, excitation force magnitude, and pendulum mass are kept identical in both systems. Specifically, primary mass M = 5 kg, natural frequency of the primary structure $ω0=k/M=4πrad/s$, pendulum mass m = 0.5 kg, force magnitude $F0=0.491N$, and ξ = 0.005 or ξ = 0.01. Note that two values of ξ are considered to examine the performance of the IPVA when the damping ratio of the primary structure changes. The remaining parameters pertaining to the autoparametric system are ca = 8.68 × 10−5N · m · s and l = 12.42 cm, which were taken from Ref. [17]. Specifically, the pendulum length was chosen to achieve autoparametric resonance, and the pendulum damping coefficient was determined from ξa = ca/(2ml2ωp) = 0.05, where $ωp=g/l$ is the natural frequency of the pendulum. The remaining parameters pertaining to the IPVA system are R, Rp, and r. Three sets of R, Rp, and r were chosen as follows: (a) R = 2.49 cm, Rp = 4.97 cm, and r = 1.99 cm; (b) R = 1.78 cm, Rp = 3.55 cm, and r = 1.42 cm; and (c) R = 2.07 cm, Rp = 4.14 cm, and r = 1.66 cm. These three sets are labeled as IPVA (a), IPVA (b), and IPVA (c), respectively, in Figs. 14 and 15. These three sets all lead to μr = 0.4 and η = 0.4 and lead to f = 0.025, f = 0.035, and f = 0.030, respectively. In this way, the dependence of the IPVA on different values of f will be examined. The rms response of the IPVA system and autoparametric system were computed using the same direct numerical integration scheme with the same settings that were used to generate Fig. 12. The effects of initial conditions were examined for both systems when their responses were observed to be nonperiodic using the same method used to obtain Fig. 12(b).

Fig. 14
Fig. 14
Close modal
Fig. 15
Fig. 15
Close modal

There are several things worth noting in Fig. 14. First, it can be clearly seen that the response curve of the primary structure displacement x for both IPVA (a) and IPVA (b) flattens for a range of ω, which demonstrates the energy transfer phenomenon for two different sets of parameters. Specifically, IPVA (b) has a more compact design (smaller R, Rp, and r) and shows better performance. Although the autoparametric system shows similar vibration suppression, both IPVA (a) and IPVA (b) outperform it. Second, let us examine the pendulum response in Fig. 14(b). As seen, for both IPVA (a) and IPVA (b), the pendulum response significantly increases within the ω range of parametric instability (ω ∈ [0.81, 0.87] for IPVA (a) and ω ∈ [0.80, 0.89] for IPVA (b)), indicating that the kinetic energy of the primary structure transfers to the pendulum, resulting in the response flattening observed in Fig. 14(a). It is noteworthy that both IPVA (a) and IPVA (b) have a larger pendulum angular velocity than the autoparametric system. Similarly, Figs. 15(a) and 15(b) show the comparison of IPVA (c) with the autoparametric system for ξ = 0.01. As can be observed, IPVA (c) outperforms the autoparametric system in terms of vibration suppression. Furthermore, it also leads to a larger pendulum angular velocity.

In addition to better vibration suppression, the IPVA system has two other advantages in comparison with the autoparametric system. First, it generates higher pendulum angular velocities, as shown in Figs. 14(b) and 15(b). Kecik and Boroweic [33] proposed an energy harvesting autoparametric system where they installed an electromagnetic generator at the pendulum pivot point to convert the pendulum angular motion into electricity. As the larger angular velocity, the larger electricity can be generated, and the larger angular velocity in the IPVA system may lead to better performance in terms of energy harvesting, which remains to be explored in the future. Second, the IPVA system leads to a more compact design. The largest length in the IPVA system is the sum Rp + r of the carrier radius and pendulum length. IPVA (a) has Rp + r = 6.69 cm, which is the maximum among the three. On the contrary, the autoparametric system requires a long pendulum (l = 12.42 cm) as it needs this length to tune the natural frequency.

## 7 Conclusion

This study analyzes the IPVA system proposed in Ref. [34] with a focus on vibration suppression of a linear oscillator subject to single harmonic excitation. It is shown that for a given excitation force magnitude, the pendulum parameters can be chosen such that parametric resonance occurs to the pendulum vibration absorber for a specific range of excitation frequencies. It is also shown that a pitchfork bifurcation and period doubling bifurcation of the pendulum response are necessary and sufficient conditions for the parametric resonance. Furthermore, when parametric resonance occurs, the kinetic energy of the linear oscillator transfers to the pendulum, resulting in vibration suppression of the linear oscillator. A saturation phenomenon similar to autoparametric vibration absorbers and nonlinear vibration absorbers is observed in the IPVA system; that is, the response of the linear oscillator saturates despite the increase in the force magnitude. Meanwhile, the increased energy due to the increase in the force magnitude seems to transfer to the pendulum, resulting in increased pendulum response. Furthermore, the system is compared to the autoparametric vibration absorber. It is shown to outperform the autoparametric vibration absorber in terms of vibration absorption and energy harvesting capabilities.

## Acknowledgment

This material is based upon work supported by the start-up funding in Michigan State University and by the National Science Foundation under Grant No. 2127495. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of Michigan State University and the National Science Foundation.

## Conflict of Interest

There are no conflicts of interest.

### Appendix A: Bifurcation Tracking Algorithm

#### Algorithm Formulation

A bifurcation tracking algorithm based on Ref. [35] is formulated as follows. First, define a vector function consisting of all the algebraic equations in Eq. (11) as follows:
$f(x^;λ)=(h0(x^;λ)h1c(x^;λ)⋮hPc(x^;λ)h1s(x^;λ)⋮hPs(x^;λ))$
(A1)
where $x^$ is the vector consisting of all the Fourier coefficients in Eq. (11) and $λ$ are bifurcation parameters of interest, e.g., $λ=η$ in Sec. 3.3 and $λ=[f,ω]T$ in Sec. 3.4. Second, define a scalar function $g(x^;λ)$ that outputs the maximum magnitude of eigenvalues of the fundamental matrix of Eq. (13). A bifurcation point will satisfy the following equations:
$h(x^;λ)=(f(x^;λ)g(x^;λ)−1)=0$
(A2)
Now suppose a solution $x^0$ to Eq. (A2) is found for a bifurcation point $λ=λ0$, e.g., by numerical integration. To obtain a neighboring bifurcation point, we use an arc length continuation method. Assume that there exists a neighboring solution $x^1=x^0+δx^0$ for a neighboring bifurcation point $λ1=λ0+δλ0$ satisfying Eq. (A2) with $‖δx^0‖2+‖δλ0‖2<<1$, where ‖.‖ depicting the ℓ2 norm.
Imposing the constancy of arc length constraint, we obtain
$‖δx^0‖2+‖δλ0‖2=‖x^1−x^0‖2+‖λ1−λ0‖2=s2$
(A3)
where s is sufficiently the small arc length. Therefore, solving Eqs. (A2) and (A3) together will give a neighboring bifurcation point.

#### Numerical Implementation

The bifurcation tracking algorithm is implemented numerically using Newton–Raphson iterations as follows. Suppose a solution $x^0$ for a bifurcation point $λ0$ is known. Define a new vector function:
$p(x^,λ;x^0,λ0)={h(x^;λ)‖x^−x^0‖2+‖λ−λ0‖2−s2}$
(A4)
which needs to equal 0 to obtain a neighboring bifurcation point. Suppose two neighboring solutions $x^k−1$ and $x^k$ are found for two bifurcation points $λk−1$ and $λk$, respectively. Denote by $y=[x^T,λT]T$ the solution to Eq. (A4). A neighboring solution is initially guessed to be
$yk+1g=yk+(yk−yk−1)=2yk−yk−1$
(A5)
Next, Eq. (A4) is solved using Newton–Raphson iterations with the initial guess. The correction to the initial guess is given by
$yk+1(n+1)=yk+1(n)−J−1p(yk+1(n);x^k,λk),n=0,1,…$
with $yk+1(0)=yk+1g$ being the initial guess and J is the Jacobian matrix. To calculate J, we use the forward difference method, i.e., the vth column of J is expressed as follows:
$J(:,v)=p(yk+1(n)+ϵev;x^k,λk)−p(yk+1(n);x^k,λk)ϵ$
(A6)
where ev is the vth column in the N × N identity matrix and N being the dimension of y and 0 < ε < <1. For the calculations in this study, ε = 10−5 was taken and the convergence criterion $‖yk+1(n+1)−yk+1(n)‖<10−10$ was used. Newton–Raphson iterations will give the value of yk+1, this (along with yk) can be further used to calculate yk+2, thereby tracking the bifurcation points.

### Appendix B: Jacobi–Anger Expansion

The Jacobi-Anger formulas used for expansion are as follows [26]:
$cos(ϕ0sinψ)=J0(ϕ0)+2J2(ϕ0)cos(2ψ)+2J4(ϕ0)cos(4ψ)+…,$
(B1)
$sin(ϕ0sinψ)=2J1(ϕ0)sin(ψ)+2J3(ϕ0)sin(3ψ)+2J5(ϕ0)sin(5ψ)+…$
(B2)
where Jm(ϕ0) is a Bessel function of the first kind of order m. For ϕ0 = 2.0 radians (115 deg), J4(ϕ0) = 0.034 and J5(ϕ0) = 0.007. Therefore, only the first two terms are required in the expansion for a good accuracy over a wide range of pendulum oscillations.

## References

1.
Smith
,
M. C.
,
2020
, “
The Inerter: A Retrospective
,”
Annu. Rev. Control, Rob., Auton. Syst.
,
3
, pp.
361
391
.
2.
Ikago
,
K.
,
Saito
,
K.
, and
Inoue
,
N.
,
2012
, “
Seismic Control of Single-Degree-of-Freedom Structure Using Tuned Viscous Mass Damper
,”
Earthquake Eng. Struct. Dyn.
,
41
(
3
), pp.
453
474
.
3.
Lazar
,
I.
,
Neild
,
S.
, and
Wagg
,
D.
,
2014
, “
Using an Inerter-Based Device for Structural Vibration Suppression
,”
Earthquake Eng. Struct. Dyn.
,
43
(
8
), pp.
1129
1147
.
4.
Lazar
,
I.
,
Neild
,
S.
, and
Wagg
,
D.
,
2016
, “
Vibration Suppression of Cables Using Tuned Inerter Dampers
,”
Eng. Struct.
,
122
, pp.
62
71
.
5.
Qian
,
F.
,
Luo
,
Y.
,
Sun
,
H.
,
Tai
,
W. C.
, and
Zuo
,
L.
,
2019
, “
Optimal Tuned Inerter Dampers for Performance Enhancement of Vibration Isolation
,”
Eng. Struct.
,
198
, p.
109464
.
6.
Marian
,
L.
, and
Giaralis
,
A.
,
2014
, “
Optimal Design of a Novel Tuned Mass-Damper–Inerter (TMDI) Passive Vibration Control Configuration for Stochastically Support-Excited Structural Systems
,”
Probabilistic. Eng. Mech.
,
38
, pp.
156
164
.
7.
De Domenico
,
D.
, and
Ricciardi
,
G.
,
2018
, “
An Enhanced Base Isolation System Equipped With Optimal Tuned Mass Damper Inerter (TMDI)
,”
Earthquake Eng. Struct. Dyn.
,
47
(
5
), pp.
1169
1192
.
8.
Joubaneh
,
E. F.
, and
Barry
,
O. R.
,
2019
, “
On the Improvement of Vibration Mitigation and Energy Harvesting Using Electromagnetic Vibration Absorber-Inerter: Exact H2 Optimization
,”
ASME J. Vib. Acoust.
,
141
(
6
), p.
061007
.
9.
Tai
,
W.-C.
,
2020
, “
Optimum Design of a New Tuned Inerter-Torsional-Mass-Damper Passive Vibration Control for Stochastically Motion-Excited Structures
,”
ASME J. Vib. Acoust.
,
142
(
1
), p.
011015
.
10.
Qian
,
F.
, and
Zuo
,
L.
,
2021
, “
Tuned Nonlinear Spring-Inerter-Damper Vibration Absorber for Beam Vibration Reduction Based on the Exact Nonlinear Dynamics Model
,”
J. Sound. Vib.
,
509
, p.
116246
.
11.
Kakou
,
P.
, and
Barry
,
O.
,
2021
, “
Simultaneous Vibration Reduction and Energy Harvesting of a Nonlinear Oscillator Using a Nonlinear Electromagnetic Vibration Absorber-inerter
,”
Mech. Syst. Signal. Process.
,
156
, p.
107607
.
12.
Yang
,
J.
,
Jiang
,
J. Z.
, and
Neild
,
S. A.
,
2020
, “
Dynamic Analysis and Performance Evaluation of Nonlinear Inerter-Based Vibration Isolators
,”
Nonlinear Dyn.
,
99
, pp.
1823
1839
.
13.
Hatwal
,
H.
,
Mallik
,
A.
, and
Ghosh
,
A.
,
1983
, “
Forced Nonlinear Oscillations of An Autoparametric System–Part 1: Periodic Responses
,”
ASME J. Appl. Mech.
,
50
(
3
), pp.
657
662
.
14.
Vyas
,
A.
, and
Bajaj
,
A.
,
2001
, “
Dynamics of Autoparametric Vibration Absorbers Using Multiple Pendulums
,”
J. Sound. Vib.
,
246
(
1
), pp.
115
135
.
15.
Bajaj
,
A.
,
Chang
,
S.
, and
Johnson
,
J.
,
1994
, “
Amplitude Modulated Dynamics of a Resonantly Excited Autoparametric Two Degree-of-Freedom System
,”
Nonlinear Dyn.
,
5
(
4
), pp.
433
457
.
16.
Hatwal
,
H.
,
Mallik
,
A.
, and
Ghosh
,
A.
,
1982
, “
Non-Linear Vibrations of a Harmonically Excited Autoparametric System
,”
J. Sound. Vib.
,
81
(
2
), pp.
153
164
.
17.
Song
,
Y.
,
Sato
,
H.
,
Iwata
,
Y.
, and
Komatsuzaki
,
T.
,
2003
, “
The Response of a Dynamic Vibration Absorber System With a Parametrically Excited Pendulum
,”
J. Sound. Vib.
,
259
(
4
), pp.
747
759
.
18.
Warminski
,
J.
, and
Kecik
,
K.
,
2009
, “
Instabilities in the Main Parametric Resonance Area of a Mechanical System With a Pendulum
,”
J. Sound. Vib.
,
322
(
3
), pp.
612
628
.
19.
Yan
,
Z.
, and
Hajj
,
M. R.
,
2015
, “
Energy Harvesting From an Autoparametric Vibration Absorber
,”
Smart Mater. Struct.
,
24
(
11
), p.
115012
.
20.
Yan
,
Z.
, and
Hajj
,
M. R.
,
2017
, “
Nonlinear Performances of an Autoparametric Vibration-Based Piezoelastic Energy Harvester
,”
J. Intell. Mater. Syst. Struct.
,
28
(
2
), pp.
254
271
.
21.
Kecik
,
K.
,
2018
, “
Assessment of Energy Harvesting and Vibration Mitigation of a Pendulum Dynamic Absorber
,”
Mech. Syst. Signal. Process.
,
106
, pp.
198
209
.
22.
Felix
,
J. L. P.
,
Balthazar
,
J. M.
,
Rocha
,
R. T.
,
Tusset
,
A. M.
, and
Janzen
,
F. C.
,
2018
, “
On Vibration Mitigation and Energy Harvesting of a Non-Ideal System With Autoparametric Vibration Absorber System
,”
Meccanica
,
53
(
13
), pp.
3177
3188
.
23.
Tan
,
T.
,
Yan
,
Z.
,
Zou
,
Y.
, and
Zhang
,
W.
,
2019
, “
Optimal Dual-Functional Design for a Piezoelectric Autoparametric Vibration Absorber
,”
Mech. Syst. Signal. Process.
,
123
, pp.
513
532
.
24.
Sharif-Bakhtiar
,
M.
, and
Shaw
,
S.
,
1992
, “
Effects of Nonlinearities and Damping on the Dynamic Response of a Centrifugal Pendulum Vibration Absorber
,”
ASME J. Vib. Acoust.
,
114
(
3
), pp.
305
311
.
25.
Detroux
,
T.
,
Renson
,
L.
,
Masset
,
L.
, and
Kerschen
,
G.
,
2015
, “
The Harmonic Balance Method for Bifurcation Analysis of Large-Scale Nonlinear Mechanical Systems
,”
Comput. Methods. Appl. Mech. Eng.
,
296
, pp.
18
38
.
26.
Newland
,
D. E.
,
1964
, “
Nonlinear Aspects of the Performance of Centrifugal Pendulum Vibration Absorbers
,”
ASME J. Manuf. Sci. Eng.
,
86
(
3
), pp.
257
263
.
27.
Hamdan
,
M.
, and
Burton
,
T.
,
1993
, “
On the Steady State Response and Stability of Non-Linear Oscillators Using Harmonic Balance
,”
J. Sound. Vib.
,
166
(
2
), pp.
255
266
.
28.
Kovacic
,
I.
,
Rand
,
R.
, and
Mohamed Sah
,
S.
,
2018
, “
Mathieu’s Equation and Its Generalizations: Overview of Stability Charts and Their Features
,”
Appl. Mech. Rev.
,
70
(
2
), p.
020802
.
29.
Xie
,
L.
,
Baguet
,
S.
,
Prabel
,
B.
, and
Dufour
,
R.
,
2017
, “
Bifurcation Tracking by Harmonic Balance Method for Performance Tuning of Nonlinear Dynamical Systems
,”
Mech. Syst. Signal. Process.
,
88
, pp.
445
461
.
30.
Carroll
,
T. L.
, and
Pecora
,
L. M.
,
1995
,
Nonlinear Dynamics in Circuits
,
World Scientific
,
Singapore
.
31.
Gourdon
,
E.
,
Alexander
,
N. A.
,
Taylor
,
C. A.
,
Lamarque
,
C.-H.
, and
Pernot
,
S.
,
2007
, “
Nonlinear Energy Pumping Under Transient Forcing With Strongly Nonlinear Coupling: Theoretical and Experimental Results
,”
J. Sound. Vib.
,
300
(
3–5
), pp.
522
551
.
32.
Oueini
,
S. S.
,
Nayfeh
,
A. H.
, and
Pratt
,
J. R.
,
1998
, “
A Nonlinear Vibration Absorber for Flexible Structures
,”
Nonlinear Dyn.
,
15
(
3
), pp.
259
282
.
33.
Kecik
,
K.
, and
Borowiec
,
M.
,
2013
, “
An Autoparametric Energy Harvester
,”
Eur. Phys. J. Spec. Top.
,
222
(
7
), pp.
1597
1605
.
34.
Gupta
,
A.
, and
Tai
,
W.-C.
,
2020
, “
Broadband and Enhanced Energy Harvesting Using Inerter Pendulum Vibration Absorber
,”
International Design Engineering Technical Conferences and Computers and Information in Engineering Conference
,
Virtual Conference
,
Aug. 17–19
, p. V007T07A007.
35.
Marathe Amol
,
C. A.
,
2006
, “
Asymmetric Mathieu Equations
,”
Proc. R. Soc. A
,
462
, pp.
1643
1659
.