Dynamic analysis of a cantilever beam coupled with bistable nonlinear energy sink: application to vibration suppression

1 Introduction

The efficient suppression of unwanted detrimental vibrations is critically important for the safe and reliable operation in various engineering disciplines, from aerospace (Araujo et al., 2025) and automobile (Xue et al., 2022a; Xue et al., 2022b) to agricultural machinery (Ch et al., 2020; Xu et al., 2019; Qian et al., 2021; Faheem et al., 2021). Over the years, a variety of vibration mitigation strategies have been developed, including active, semi-active, and passive vibration control. Linear passive dynamic vibration absorbers have been widely adopted due to the structural simplicity and ease of implementation. However, the effective operation of a linear device relies on precise tuning of system parameters to achieve fundamental resonance between the primary system and the linear vibration absorber. This dependency adversely affects the vibration reduction performance when subjected to broadband or frequency-varying excitations, thereby restricting its practical applications. To address these limitations, the introduction of essential (nonlinearizable) nonlinearity can effectively dissipate the vibrational energy over a wider frequency range because of the unfixed natural frequencies of the nonlinear energy sink (NES), demonstrating significant potential for application in diverse vibration control scenarios (Xu et al., 2023; Haris et al., 2020). Additionally, nonlinear systems can exhibit fruitful dynamical behaviors, including bifurcation phenomena, sensitivity to system parameters, and multiscale wave dynamics (Nadeem and Iambor, 2025; Nadeem et al., 2026a; Nadeem et al., 2025a; Nadeem et al., 2026b; Nadeem et al., 2025b).

A typical nonlinear energy sink consists of a lightweight substructure coupled to the primary system via a nonlinear spring or damper. The novel phenomenon, targeted energy transfer (TET) mechanism is triggered by the transient resonance captures (Gourdon et al., 2007; Wei et al., 2019). Specifically, vibrational energy from the impulsively excited main structure (donor) is rapidly transferred to the nonlinear attachment (recipient), and the energy is then localized and ultimately dissipated by the damping element of the nonlinear substructure. When a primary system subjected to sustained harmonic excitation is coupled with a nonlinear energy sink, the strongly nonlinear system may exhibit a novel form of nonlinear dynamic response, namely, strongly modulated response (SMR) (Gendelman et al., 2008; Liu et al., 2019). The exact fundamental resonance leads to the relaxation oscillations generated by the rapid jumps between disjoint stable branches of the slow invariant manifold (SIM) (Liu and Wang, 2023). Therefore, the SMR regime is regarded as the steady-state analogue of targeted energy transfer between the linear oscillator and the nonlinear attachment. Accordingly, the ubiquitous SMR regimes yield efficient and broadband vibration suppression performance. Various types of NES configurations have been proposed to enhance vibration mitigation performance, for example, cubic (Li et al., 2022; Li et al., 2025), vibro-impact (Lizunov et al., 2025), non-smooth (Chen J. E. et al., 2020), multistable (Hu et al., 2025), multi-degree-of-freedom (Liu and Wang, 2023) and magnetic inerter (Sun et al., 2025) NESs. These diverse designs demonstrate effective and robust control of structural vibrations for various systems, and the performance of the NES is significantly improved by efficient and irreversible energy dissipation (Gendelman et al., 2000; Vakakis and Gendelman, 2000; Ding and Chen, 2020).

For the continuous systems, NES exhibits superior vibration mitigation performance in diverse practical engineering applications. Georgiades and Vakakis (2007) demonstrated that the NES can passively absorb and locally dissipate up to 87% of broadband shock-induced vibration energy through irreversible targeted energy transfer (TET). Philip et al. (2023) showed that the NES with nonlinear damping efficiently suppresses the vibration even at low excitation amplitudes through the SMR regimes for the harmonically excited fluid-conveying pipes. Ahmadabadi and Khadem (2012) studied nonlinear normal modes (NNMs) of a cantilever beam coupled with a grounded or ungrounded NES, and they further revealed that the topological structure and bifurcations of NNMs play a critical role in triggering one-way irreversible energy pumping in continuous systems. Zhang et al. (2024) proposed a one-third power nonlinear energy sink for micro-vibration suppression of satellite solar panel, demonstrating its high effectiveness when the system exhibits strongly modulated response. The bifurcation analysis showed that the presence of folding singularities yields the efficient targeted energy transfer for micro-vibration mitigation of flexible structures. Kumar and Kumar (2024) showed that the two-degree-of-freedom NES outperforms the traditional single-degree-of-freedom NES under high excitation levels, and the first resonant peak amplitude was reduced by up to 95%. Chen et al. (2018) demonstrated that the parallel NES presented superior vibration suppression performance in comparison to a single NES of equivalent mass, and the steady-state high amplitude isolated branches can be further eliminated for harmonic excitation. Liu et al. (2024) show that both centralized and distributed NES cells can effectively suppress the high-order modal responses of Timoshenko beams, and the attenuation efficiency enhances as the number of NES cells increases.

On the other hand, the external excitation should surpass a well-defined energy threshold (Dekemele et al., 2019; Romeo et al., 2015), otherwise effective energy transfer becomes unachievable for the NES with cubic stiffness. Accordingly, the energy threshold for the generation of SMR should be reduced to improve the vibration mitigation performance of NES over a wider range of external excitations. It is found that the bistable nonlinear energy sink (BNES) with a negative stiffness mechanism can effectively lower the SMR energy threshold for the lower excitation levels (Liao et al., 2023; Philip et al., 2022). Fang et al. (2017) proposed a continuous bistable nonlinear energy sink which was realized by a cantilever beam incorporating partial constrained layer damping and magnetostatic interaction, which significantly enhances the efficiency of energy transfer and dissipation under broadband shock excitations. Zeng et al. (2022) demonstrated that the energy threshold is effectively reduced and the robustness of the BNES is simultaneously enhanced by the buckling beam. Habib and Romeo (2017) investigated the effectiveness of BNES by the geometry and topological structure of the analytically obtained slow invariant manifold, and the tuned BNES demonstrates a robust capability to suppress oscillations of the primary system across a broadband impulsive excitations and multiple vibrational modes. Chen Y-Y. et al. (2020) developed a magnetic BNES that robustly protects the shear frame for seismic mitigation because of the immediate cascade of broadband transient internal resonance captures. Wang et al. (2024) proposed a novel track BNES that mitigates and presents robust seismic response mitigation despite the variation of impulsive energy levels.

More recently, Wang and Ding (2023) investigated TET in a nonlinear oscillator coupled with a BNES by nonlinear normal modes and frequency-energy plots and reveal that asymmetric backbone branches and multiple resonance captures contribute significantly to its superior energy dissipation performance. Li et al. (2023) investigated the Hamiltonian dynamics and TET of a grounded BNES and found that asymmetric potential well can enhance energy exchange and suppress reverse energy transfer from the NES to the primary oscillator. Nasar and AL-Shudeifat (2024) proposed a novel spring-loaded inverted pendulum-based BNES and demonstrated its robust vibration suppression and superior energy dissipation performance. Kumar and Kumar (2024) investigated the multi-mode vibration suppression of a cantilever beam carrying an unbalanced rotor by a NES and demonstrated substantial reductions in multiple resonant responses through optimized absorber parameters. Chen et al. (2025a) proposed a cellular bistable track nonlinear energy sink and demonstrated that the number of coupled cells and system parameters strongly influence the SMR activation threshold and vibration-suppression performance. Sun et al. (2026) investigated two-degree-of-freedom NESs and demonstrated their superior vibration-suppression performance over a relatively wide excitation range.

Adding to the body of the aforementioned study, in this work, we aim to investigate the complex multiscale nonlinear dynamic responses of a flexible Euler-Bernoulli beam coupled to a bistable NES and develop a unified analytical framework that combines the complexification-averaging method, the method of multiple scales, local bifurcation analysis, and slow invariant manifold analysis to reveal the formation of strongly modulated responses and the associated mechanism of steady-state targeted energy transfer. The numerical simulations and analytical asymptotic prediction were performed to investigate the vibration mitigation performance of the bistable NES when the slender beam is subjected to sustained harmonic excitations. The Galerkin method was applied to truncate the nonlinear governing equations of the beam by mode decomposition. The complexification-averaging method (CX-A) is subsequently implemented to analytically predict the global system behavior when the beam is subjected to harmonic excitation. Moreover, the method of multiple scales and bifurcation analysis were applied to show the complex dynamics of the nonlinear system at fast and slow time scales. The results indicate the efficient vibration mitigation performance by the proposed BNES. The outline of this paper is structured as follows, Section 2 presents the physical model of the continuous system and the corresponding reduced-order governing equations of a cantilever beam coupled with a bistable nonlinear energy sink (BNES). After that, in Section 3, the CX-A and multiple scales methods were applied to investigate the global system behavior at different time scales when the system is under sustained harmonic forcing, and the corresponding bifurcation analysis and mechanism of steady-state targeted energy transfer are revealed. Section 4 evaluates the vibration mitigation performance of the proposed BNES and the parametric studies are performed. Finally, the main conclusions of this work are summarized in the last section.

2 Model formulation and reduced-order modeling of the coupled system

2.1 Physical model of a cantilever beam with a bistable nonlinear energy sink

The schematic diagram of the continuous physical model is illustrated in Figure 1. The considered system consists of a slender linear cantilever beam and a coupled bistable nonlinear energy sink. The uniform and homogeneous beam is subjected to the external load F(t) at position a, and the BNES is coupled to the beam at position d. The length, density and cross-sectional area of the beam are denoted by L, ρ and A, respectively. With the assumption of distributed damping, the damping coefficient of the beam is μ. The Young’s modulus and moment of inertia of the beam are E and I respectively, and EI is the bending stiffness. The distributed mass ρA and bending stiffness EI determine the natural frequencies of the cantilever beam, while the damping coefficient μ characterizes the structural energy dissipation of the beam. The continuous system considered in this study is a slender linear cantilever beam. Since we focus on the transverse vibration of the beam around the first natural frequency, the Euler-Bernoulli beam theory is adopted by neglecting transverse shear deformation and rotary inertia. The governing equations of the cantilever beam coupled with BNES are derived by the Hamilton’s principle, which are expressed as.where the transverse displacement of the beam is represented by w (x,t), which is a function of time t and the spatial coordinate x along the axis in the initial configuration, α denotes the mass of the BNES, y the displacement of the NES, C the damping coefficient of BNES, K ₁ the negative stiffness of BNES, and K ₃ the cubic stiffness of BNES, the δ(x) Dirac delta function. The damping coefficient c characterizes the viscous energy dissipation in the BNES. It is noted that the bistability of the nonlinear attachment arises from the negative stiffness around the trivial equilibrium position. The positive nonlinear stiffness provides large-displacement stabilization, while the negative linear stiffness and positive nonlinear stiffness together determine the stable equilibrium positions and potential-barrier height of the BNES.

ρ A ∂ 2 w ∂ t 2 + μ ∂ w ∂ t + E I ∂ 4 w ∂ x 4 + K 3 w d , t − y 3 + K 1 w d , t − y + C w ˙ d , t − y ˙ δ x − d = F t δ x − a α y ¨ + K 3 y − w d , t 3 + K 1 y − w d , t + C y ˙ − w ˙ d , t = 0 (1)

2.2 Truncated model by Galerkin method

The transient dynamics of the continuous system are investigated by the projection of the complex responses in the infinite-dimensional orthonormal basis, i.e., the linear vibration modes of the linear cantilever beam. Accordingly, with the application of the method of separation of variables, the transverse displacement of the beam is approximately expressed as the superposition of intrinsic mode function, and therefore, can be written in the discretized formwhere ϕ _*n*(x) and q _*n*(t) denote the n^th vibration mode function and the corresponding generalized modal displacement of the cantilever beam. Additionally, the n^th vibration mode function and the natural frequency ω _*r* corresponding to it are given by where Equations 3, 4 define the vibration mode function and the corresponding natural frequency, respectively, and

w x , t = ∑ n = 1 ∞ ϕ n x q n t (2)

ϕ r x = cos λ r L x − cosh λ r L x + ς r sin λ r L x − sinh λ r L x (3)

ω r = λ r 2 E I ρ A L 4 (4)

ς r = sin ⁡ λ r − sinh ⁡ λ r cos ⁡ λ r + cosh ⁡ λ r (5)

where Equation 5 defines the coefficient λr used in the mode function. The eigenvalue λ _*r* of (Equation 3) can be obtained by solving the following algebraic equation

1 + cos ⁡ λ r ⁡ cosh ⁡ λ r = 0 (6)

The eigenvalue λr in Equation 3 can be obtained by solving Equation 6. Substituting Equation 2 into the governing equations of the continuous system Equation 1 and using orthonormality relations, the partial differential equations are transformed into an infinite set of ordinary differential equationswhere b _*fh* and b _*sh* denote the generalized modal mass and modal stiffness, respectively.

b f h q ¨ s t + μ ρ A b f h q ˙ s t + b s h q s t + { K 1 ∑ r = 1 ∞ q r t ϕ r d − y + K 3 ∑ r = 1 ∞ q r t ϕ r d − y 3 + C ∑ r = 1 ∞ q ˙ r t ϕ r d − y ˙ } ϕ s d = F t ϕ s a α y ¨ + K 1 y − ∑ r = 1 ∞ q r t ϕ r d + K 3 y − ∑ r = 1 ∞ q r t ϕ r d 3 + C y ˙ − ∑ r = 1 ∞ q ˙ r t ϕ r d = 0 (7)

b f h = ρ A ∫ 0 L ϕ s x 2 d x , b s h = E I ∫ 0 L ϕ s x ϕ s ⁗ x d x (8)

Equation 8 defines the generalized modal mass and modal stiffness of the cantilever beam. In this work, we aim to investigate the complex nonlinear dynamics of the continuous system under sustained harmonic excitation F(t) = A ₀cosωt. The fundamental resonance is of interest when the frequency of the harmonic excitation is in the neighborhood of the first natural frequency of the cantilever beam.

The governing Equation 7 are then numerically integrated by the fourth-order Runge-Kutta method with zero initial conditions and the system parameters are set as α = 0.1 kg, ρ = 7.85 × 10³ kg/m³, a = 0.33 m, A = 1.05 × 10⁻³ m², EI = 4750 Nm², L = 1.1m, d = 0.7m, A ₀ = 100 N, k ₁ = −100N/m, k ₃ = 1 × 10⁶ N/m³, c = 10 Ns/m, and μ = 0.5 Ns/m. The convergence of the Galerkin truncation is examined by retaining different approximate bending modes N = 1, 2, and 3, and Figure 2a illustrates the frequency-response curves in the vicinity of the first natural frequency of the cantilever beam. Figures 2b–d further compare the beam displacement at the BNES attachment point at three different excitation frequencies, corresponding to representative small-amplitude periodic, large-amplitude periodic, and strongly modulated responses, respectively. The nearly indistinguishable results obtained using the three modal truncations in both the time and frequency domains indicate that the contributions of higher-order modes are negligible within the considered frequency range around the fundamental resonance.

Therefore, the one-mode approximation of the continuous system can accurately predict the system steady-state dynamics. Accordingly, the one-mode approximation of the continuous system is adopted (Zhang et al., 2024; Mamaghani et al., 2016; Mamaghani et al., 2018; Chen et al., 2025b) and the one-mode approximate model of the governing equations of the nonlinear system (7) is expressed aswhere q ₁(t) and q ₁(t)ϕ ₁(d) denote the displacement of the first mode and the displacement of the cantilever beam at the location of the BNES.

b f h q ¨ 1 t + μ ρ A b f h q ˙ 1 t + b s h q 1 t + { K 1 q 1 t ϕ 1 d − y + K 3 q 1 t ϕ 1 d − y 3 + C q ˙ 1 t ϕ 1 d − y ˙ } ϕ 1 d = F t ϕ 1 a α y ¨ + K 1 y − q 1 t ϕ 1 d + K 3 y − q 1 t ϕ 1 d 3 + C y ˙ − q ˙ 1 t ϕ 1 d = 0 (9)

The first equation of Equation 9 is multiplied by ϕ ₁(d), additionally, by simple algebraic manipulations, the equations of motion of the approximate model are written in a simpler formwherewhere Equation 11 defines the transformed coordinate and the system parameters used in Equation 10, and z denotes the transverse displacement of the beam at the BNES attachment point.

z ¨ + ω 1 2 z + ε ζ z ˙ + ε k 3 z − y 3 + ε k 1 z − y + ε c z ˙ − y ˙ = ε f ⁡ cos ⁡ ω t ε y ¨ + ε k 3 y − z 3 + ε k 1 y − z + ε c y ˙ − z ˙ = 0 (10)

z = q 1 t ϕ 1 d , α k 1 = K 1 , α k 3 = K 3 , α c = C , ε = α ϕ 1 d 2 b f h , ζ = μ b f h ρ A α ϕ 1 d 2 , f = A 0 ϕ 1 a α ϕ 1 d (11)

For the sake of brevity, the following variables are introducedwhere v(t) of Equation 12 represents the displacement of the mass center of the nonlinear system and w(t) the relative displacement between the BNES and the beam at the location of the BNES. To analytically examine the steady-state responses of the nonlinear system, the complexification-averaging method (CX-A) developed by Manevitch (2001) is implemented to predict the steady-state varying amplitudes of the system responses. With the assumption that the continuous system is in the steady-state fundamental resonance regime with the sustained external harmonic forcing, then the dominant oscillating frequencies of the cantilever beam and BNES are identical to the frequency of harmonic excitation. Accordingly, the slow- and fast-varying parts of the steady-state forced response are separated and the global dynamics are then reduced to the important slow-flow manifold. Therefore, to apply the CX-A method, the following complex variables are introduced to perform the slow/fast partitionswhere the complex functions φ ₁(t) and φ ₂(t) represent the slowly varying amplitudes of the mass center and relative displacement of the cantilever beam coupled with BNES, respectively. In addition, the fast-oscillating terms e ^*iωt* are modulated by the corresponding slowly varying envelopes and the dominant fast oscillating frequency is identical to the frequency of applied harmonic excitation. On the other hand, in comparison to the cantilever beam, the mass of the lightweight BNES is sufficiently small (i.e., 0<ε ≤ 1), and the mass ratio ε can be considered as the perturbation factor for the following asymptotic analysis. Substituting Equation 13 into Equation 10 and keeping the terms containing e ^*iωt*, the resulting equations are subsequently averaged with respect to the dominant fast frequency ω. The governing equations given in Equation 10 for real variables can be equivalently expressed in the complex form in terms of the complex functions φ ₁(t) and φ ₂(t). Consequently, one obtains the complex differential equations governing the slow modulation of the nonlinear system

v t = z t + ε y t , w t = z t − y t (12)

ψ 1 t = v ˙ + i ω v = φ 1 t e i ω t , ψ 2 t = w ˙ + i ω w = φ 2 t e i ω t (13)

φ ˙ 1 + i ω 2 φ 1 + ζ ε ω − i ω 1 2 2 1 + ε ω φ 1 + ε φ 2 = ε f 2 φ ˙ 2 + ε 2 c + ζ + c 2 ε + 1 φ 2 + ε ζ φ 1 2 1 + ε − i ω 1 2 2 1 + ε ω φ 1 + ε φ 2 − i k 1 2 1 + ε ω φ 2 + 2 ε φ 2 + ε 2 φ 2 + i ω 2 φ 2 − 3 i k 3 1 + ε φ 2 2 φ 2 8 ω 3 = ε f 2 (14)

3 Nonlinear system behavior at fast and slow time scale

3.1 Stability and bifurcation analysis at fast time scale

Different types of steady-state responses of the original nonlinear continuous system, for example, periodic or quasi-periodic, can be analytically predicted by the slow flow (14). The asymptotic analysis of the forced responses of the nonlinear system Equation 14 starts by investigating the stationary points, which describes the steady-state periodic motion. Accordingly, the fixed points are obtained by setting the derivatives of the complex variables equal to zero (i.e., φ ˙ 1 = 0 and φ ˙ 2 = 0), then one has the following algebraic equations that describe the steady-state motion under the condition of fundamental resonancewhere φ ₁₀ and φ ₂₀ denote the fixed points of the complex functions φ ₁(t) and φ ₂(t) respectively, and provides predictions for the amplitudes of the mass center and relative motion between the cantilever beam and NES. The variation of the relative motion effectively affects the vibration absorption efficiency of BNES, and a larger amplitude of the relative motion synchronized with the small motion of mass center indicates a stronger capacity of the nonlinear vibration absorber. Accordingly, the relative motion between the continuous system and the BNES is of interest, and the bifurcation analysis of the steady-state amplitudes of the relative motion is performed. Referring to the first equation of Equation 15, expressing φ ₁₀ in terms of φ ₂₀ and substituting it into the second equation of Equation 15, one has

i ω 2 φ 10 + ζ ε ω − i ω 1 2 2 1 + ε ω φ 10 + ε φ 20 = ε f 2 ε 2 c + ζ + c 2 ε + 1 φ 20 + ε ζ φ 10 2 1 + ε − i ω 1 2 2 1 + ε ω φ 10 + ε φ 20 − i k 1 2 1 + ε ω φ 20 + 2 ε φ 20 + ε 2 φ 20 + i ω 2 φ 20 − 3 i k 3 1 + ε φ 20 2 φ 20 8 ω 3 = ε f 2 (15)

1 ω 3 ε ( i ω 4 − − c − ζ ε − c ω 3 − i c ζ + k 1 ε + ω 1 2 + k 1 ω 2 − c ω 1 2 + ζ ε k 1 ω + i k 1 ω 1 2 ) φ 20 + 3 k 3 4 ω 5 ε − i 1 + ε ω 2 − ζ ε ω + i ω 1 2 φ 20 2 φ 20 = f (16)

By performing simple algebraic manipulations, Equation 16 can be written in the following compact form,where

α 1 Z + α 2 Z 2 + α 3 Z 3 + α 4 = 0 (17)

Z = φ 20 2 , α 1 = a 3 2 + a 4 2 , α 2 = 2 a 1 a 3 + a 2 a 4 , α 3 = a 1 2 + a 2 2 , α 4 = − a 5 2 , a 1 = − 3 k 3 ζ 4 ω 4 , a 2 = − 3 k 3 4 ω 3 − 3 k 3 4 ω 3 ε + 3 k 3 ω 1 2 4 ω 5 ε , a 3 = c + ζ + c ε − c ω 1 2 ω 2 ε − k 1 ζ ω 2 , a 4 = ω ε − c ζ + k 1 ω − ω 1 2 + k 1 ω ε + k 1 ω 1 2 ω 3 ε , a 5 = f . (18)

The coefficients in Equation 17 are given by Equation 18. The steady-state amplitude of the relative motion between the cantilever beam and BNES can be obtained by solving the polynomial (17) in terms of the modulus of φ ₂₀. The solutions of algebraic cubic Equation 17 depends on the corresponding values of coefficients, and one or three real solutions may exist. Accordingly, at certain critical values of the parameters, two fixed points might collide and annihilate each other and the resulting saddle-node (SN) bifurcation is expected. To determine the loci of SN bifurcation points, the derivative of Equation 17 with respect to the square of modulus Z should be equal to zero

3 α 3 Z 2 + 2 α 2 Z + α 1 = 0 (19)

Accordingly, the SN bifurcation points should simultaneously satisfy Equations 17, 19 and can be computed in the parameter space by combining the two equations. Then the two-dimensional space plane of parameters (c,f) is expressed aswhere

β 1 a 5 4 + β 2 a 5 2 + β 3 = 0 (20)

β 1 = 3 9 a 1 2 + 9 a 2 2 2 a 1 2 + a 2 2 β 2 = 36 a 1 2 + a 2 2 a 1 a 3 + a 2 a 4 a 3 2 + 9 a 4 3 a 1 2 − 16 a 2 a 3 a 4 a 1 + a 2 2 9 a 3 2 + a 4 2 β 3 = 36 a 3 2 + a 4 2 2 a 1 2 + a 2 2 a 1 a 4 − a 2 a 3 2 (21)

Equations 20, 21 give the algebraic condition for the saddle-node bifurcation boundary in the parameter plane and define the corresponding coefficients, respectively. In this work, the fundamental resonance between the nonlinear system and external harmonic excitation is considered, therefore, the frequency of the harmonic excitation is in the neighborhood of the first natural frequency of the cantilever beam, i.e., ω = ω ₁ (1+εσ), where σ is the frequency detuning parameter. Figure 3 presents the representative SN bifurcation curve in the two-dimensional parameter plane (c,f), and the other system parameters are assigned as ε = 0.01kg, ζ = 5 Ns/m, k ₁ = -100N/m, k ₃ = 1 × 10⁶ N/m³, σ = 1. Clearly, the boundary of the SN bifurcation is characterized by a triangular shape, and the whole parameter plane is divided into two regions. There are three real solutions inside the region boundary, on the contrast, only one steady-state periodic solution exists outside the region. Figure 4 illustrates the SN bifurcation curves for varying positive and negative frequency detuning parameters. Obviously, the variation of frequency detuning parameter does not alter the triangular geometry of the SN bifurcation boundaries. Specifically, for the positive frequency detuning parameter, an increase in the frequency detuning parameter results in an expansion of the zone of existence for three real solutions, indicating that multiple steady-state periodic solutions can be achieved for a broader range of forcing amplitudes. On the other hand, the decrease of negative frequency detuning parameter leads to the upward shift of the bifurcation region.

While the SN bifurcation characterizes how external excitation or damping affect the number of real equilibrium points, it is the instability of an equilibrium point that triggers the Hopf bifurcation, subsequently leading to SMR, which is the favorable kind of steady-state responses for the efficient energy transfer and vibration absorption by the BNES. Therefore, a stability analysis of the equilibrium points is essential in the slow flow (14). To perform the analysis of Hopf bifurcation, the small complex perturbations are introduced near the fixed pointswhere δ ₁(t) and δ ₂(t) denote the distances around the equilibrium points φ ₁₀ and φ ₂₀ respectively. The linearized ordinary differential equations governing the slow modulation of the complex perturbations δ ₁(t) and δ ₂(t) around the fixed points are obtained by substituting (Equation 22) into (Equation 14) and neglecting high-order perturbation termswhere asterisk represents the conjugate. The characteristic polynomial of the Jacobian matrix of linearized system (Equation 23) is expressed aswhere µ is the eigenvalue and the coefficients γ i , i = 1 … 4 of Equation 24 are given in Supplementary Appendix A. At the Hopf bifurcation point, the characteristic polynomial should possess a pair of purely imaginary roots, which is expressed as

φ 1 t = φ 10 + δ 1 t , δ 1 t ≪ 1 φ 2 t = φ 20 + δ 2 t , δ 2 t ≪ 1 (22)

δ ˙ 1 = − i ω 2 + i ω 1 2 − ζ ε ω 2 1 + ε ω δ 1 + i ε ω 1 2 − ζ ω ε 2 2 1 + ε ω δ 2 δ ˙ 1 * = i ω 2 − i ω 1 2 + ζ ε ω 2 1 + ε ω δ 1 * − i ε ω 1 2 + ζ ω ε 2 2 1 + ε ω δ 2 * δ ˙ 2 = i ω 1 2 − ζ ω ε 2 1 + ε ω δ 1 + 3 i k 3 1 + ε ϕ 20 2 4 ω 3 δ 2 − i ω 2 δ 2 − c + ζ ε 2 + 2 ε c + c 2 1 + ε δ 2 + i k 1 ε 2 + ε ω 1 2 + 2 ε k 1 + k 1 2 1 + ε ω δ 2 + 3 i k 3 1 + ε ϕ 20 2 8 ω 3 δ 2 * δ ˙ 2 * = − i ω 1 2 + ζ ω ε 2 1 + ε ω δ 1 * − 3 i k 3 1 + ε ϕ 20 2 4 ω 3 δ 2 * + i ω 2 δ 2 * − c + ζ ε 2 + 2 ε c + c 2 1 + ε δ 2 * − i k 1 ε 2 + ε ω 1 2 + 2 ε k 1 + k 1 2 1 + ε ω δ 2 * − 3 i k 3 1 + ε ϕ 20 2 8 ω 3 δ 2 (23)

µ 4 + γ 1 µ 3 + γ 2 µ 2 + γ 3 µ + γ 4 = 0 (24)

μ = ± j Ω (25)

Subsequently, substituting Equation 25 into the characteristic Equation 24 and equating the real and imaginary parts to zero yields

Ω 4 − γ 2 Ω 2 + γ 4 = 0 Ω γ 1 Ω 2 + γ 3 = 0 (26)

By eliminating Ω in the Equation 26, one obtains the following algebraic equation

γ 3 2 − γ 2 γ 3 γ 1 + γ 4 γ 1 2 = 0 (27)

After simple algebraic manipulations, Equation 27 can be expressed in the polynomial form with respect to the amplitude Z, which is written aswhere the coefficients are given in Supplementary Appendix B.

v 5 Z 4 + v 4 Z 3 + v 3 Z 2 + v 2 Z + v 1 = 0 (28)

On the other hand, the Hopf bifurcation points should also satisfy Equation 17, for a specific damping coefficient c, the corresponding amplitude of the harmonic excitation f can be computed by combining Equation 17 and Equation 28. Accordingly, the Hopf bifurcation boundaries can also be computed in the two-dimensional parameter plane (c,f). It is noted that the occurrence of a Hopf bifurcation indicates the loss of stability of the steady-state responses, indicating a transition from periodic response to the quasi-periodic or strongly modulated response (SMR). Figure 5a illustrates the SN and Hopf bifurcation curves in the (c,f) plane for the system parameters ε = 0.01kg, ζ = 5 Ns/m, k ₁ = -100N/m, k ₃ = 1 × 10⁶ N/m³, σ = 1.2, and the solid and dashed curves denote the Hopf and SN bifurcation boundaries, respectively. The results show that the Hopf and saddle-node bifurcations coexist in certain forcing amplitude ranges. Moreover, Figure 5b presents the amplitude responses of the nonlinear system for c = 5 Ns/m, and the unstable branches are marked by red asterisks. Obviously, these two bifurcation types, both leading to a loss of stability, coincide over certain forcing amplitude ranges.

Figure 6 shows the frequency response curve and bifurcation behavior of relative motion, and the system parameters are assigned as ε = 0.01kg, ζ = 5 Ns/m, k ₁ = -100N/m, k ₃ = 1 × 10⁶ N/m³ and c = 5 Ns/m. The asterisk denotes the unstable steady-state quasi-periodic responses, which corresponds to the SMR. The SMR regime contributes effectively to vibration mitigation due to the irreversible energy transfer and efficient dissipation, and the underlying mechanism of energy transfer will be examined in details during the following analysis. The above local bifurcation analysis predicts the stability of equilibrium points and dynamic evolution of system responses with the variation of system parameters. In the following analysis, to fully understand the dynamic evolution and energy transfer of quasi-periodic responses, the global nonlinear analysis is further conducted to investigate the underlying mechanisms.

3.2 Strongly modulated responses at slow time scale

Strongly modulated responses (SMRs) have been identified as an effective mechanism for steady-state targeted energy transfer and also yields relaxation oscillations in the condition of 1:1 resonance between the system and external harmonic excitation. To investigate the energy transfer mechanism of the quasi-periodic response of the system under harmonic excitation, referring to the slow flow (14), φ ₁(t) can be expressed in terms of φ ₂(t). For the sake of brevity, the two first-order differential equations governing the slow modulations of envelopes φ ₁(t) and φ ₂(t) can be reduced into a single second-order complex modulation equation with respect to φ ₂(t), which is expressed aswhere

d 2 φ 2 d t 2 + d d t η 1 + η 3 φ 2 − 3 i k 3 1 + ε 8 ω 3 φ 2 2 φ 2 − 3 k 3 η 3 1 + ε j 8 ω 3 ϕ 2 2 ϕ 2 + η 2 2 ε − η 1 η 3 ϕ 2 = η 2 + η 3 ε f 2 (29)

η 1 = 1 2 1 + ε ω ( c + ζ ω − i k 1 ε 2 + i ω 2 − i ω 1 2 − 2 i k 1 + 2 c ω ε + i ω 2 − i k 1 + c ω ) , η 2 = i ω 1 2 − ε ω ζ 2 1 + ε ω , η 3 = i ω 2 1 + ε + ε ω ζ − i ω 1 2 2 1 + ε ω (30)

where Equation 30 defines the coefficients in the Equation 29. We recall the condition of the fundamental resonance between the nonlinear system and external harmonic excitation, and the frequency of the harmonic excitation is in the vicinity of the first natural frequency of the cantilever beam, i.e., ω = ω ₁ (1+εσ).

Proceeding to the next step, the method of multiple scales is implemented to study the dynamic evolution of the slow flow at different time scales, and a series of new fast and slow time scales are introduced as τ n = ε n t , n = 0 , 1 , where τ ₀ and τ ₁ represent the fast and slow time scales respectively. Therefore, φ ₂(t) can be expressed by the new independent time scales, which is written as

φ 2 = φ 2 τ 0 , τ 1 , … , d d t = ∂ ∂ τ 0 + ε ∂ ∂ τ 1 + ⋯ , d 2 d t 2 = ∂ 2 ∂ τ 0 2 + 2 ε ∂ 2 ∂ τ 0 ∂ τ 1 + ⋯ (31)

Substituting Equation 31 into Equation 29 and equating the coefficients of the terms at the same power ε, one obtains the hierarchy of subproblems at fast and slow time scales

O ε 0 : ∂ 2 φ 2 ∂ τ 0 2 + ∂ ∂ τ 0 c φ 2 2 + i φ 2 ω 1 2 − i k 1 φ 2 2 ω 1 − 3 i k 3 8 ω 1 3 φ 2 2 φ 2 = 0 O ε 1 : 2 ∂ 2 φ 2 ∂ τ 0 ∂ τ 1 + ∂ ∂ τ 1 c φ 2 2 + i ω 1 φ 2 2 − i k 1 φ 2 2 ω 1 − 3 i k 3 8 ω 1 3 φ 2 2 φ 2 + ∂ ∂ τ 0 c + ζ φ 2 2 + 3 i ω 1 σ φ 2 2 + i k 1 σ − 1 φ 2 2 ω 1 + 3 i k 3 3 σ − 1 8 ω 1 3 φ 2 2 φ 2 + 3 2 σ + 1 ω 1 − i ζ k 3 16 ω 1 3 φ 2 2 φ 2 + − ω 1 2 σ 2 + i 2 c σ + c + ζ ω 1 4 + c ζ + 2 k 1 σ + k 1 4 − i ζ k 1 4 ω 1 φ 2 = j 4 f ω 1 . . . (32)

The first equation of (32) describes the dynamic evolution of the nonlinear system at the fast time scale, and the integration of the governing equations is performed to obtain the leading-order approximation of the slow flowwhere C τ _ 1 ) , τ _ 2 , … denotes the integration constant with respect to the fast time scale. The equilibrium points of the first integral are obtained by setting the time derivative in Equation 33 equal to zero, which is expressed as

∂ ∂ τ 0 φ 2 + c φ 2 2 + i ω 1 φ 2 2 − i k 1 φ 2 2 ω 1 − 3 i k 3 8 ω 1 3 φ 2 2 φ 2 = C τ 1 , τ 2 , … (33)

c 2 φ 2 τ 1 + i ω 1 2 φ 2 τ 1 − i k 1 2 ω 1 φ 2 τ 1 − 3 i k 3 8 ω 1 3 φ 2 τ 1 2 φ 2 τ 1 = C τ 1 (34)

The complex algebraic equation provides the projection of slow invariant manifold (SIM), which captures the equilibrium points of the nonlinear system at the fast time scale. It is noted that the steady-state amplitude only depends on the slower time scales, i.e., φ _ ( 2 ) = φ _ ( 2 ) τ _ 1 , …. In addition, the polar coordinate transformation of the equilibrium points is written as φ 2 τ 1 = N τ 1 e i θ τ 1, and the substitution of the polar form of the fixed points into Equation 34 yields the real algebraic expression of SIMwhere

p 1 Z 3 + p 2 Z 2 + p 3 Z = 4 C 2 ω 1 2 (35)

Z τ 1 = N 2 τ 1 , p 1 = 9 k 3 2 16 ω 1 8 , p 2 = 3 k 1 k 3 − 3 k 3 ω 1 2 2 ω 1 6 , p 3 = ω 1 4 + ω 1 2 c 2 − 2 ω 1 2 k 1 2 + k 1 2 ω 1 4 (36)

The number of roots of Equation 35 depends on the coefficients of system parameters in Equation 36 and the energy-related integral constant |C (τ ₁)|. The number of real roots for the cubic Equation 35 is determined by the monotonicity of its homogeneous part. The polynomial is strictly monotonic and has only one root if the derivative of its homogeneous part has no real solutions. Conversely, if this derivative has real solutions, the system becomes non-monotonic and may exhibit either one or three real roots, depending on the constant of integral. Distinction between the two cases is investigated by examining whether the derivative of the homogeneous part of the polynomial with respect to Z has real roots, accordingly, the derivative of Equation 35 yields

3 p 1 Z 2 + 2 p 2 Z + p 3 = 0 (37)

The roots of Equation 37 are expressed as:

Z 1 , 2 = N 1 , 2 2 = − p 2 ∓ p 2 2 − 3 p 1 p 3 3 p 1 (38)

where Equation 38 provides the amplitudes associated with the two folding points of the SIM. If Equation 35 has a single real root, no folding singularities exist on the entirely stable SIM, thereby the steady-state response is finally attracted to the stable fixed points. When the polynomial function attains one local maximum and one local minimum, then the SIM is composed of two stable branches and one unstable branch with respect to the fast time scale. Then steady-state response can either be attracted to the stable orbit in phase space or exhibit persistent relaxation oscillations due to the jump phenomenon between disjoint stable branches on the SIM.

Figure 7a presents the typical topology of SIM, and the system parameters are assigned as ε = 0.01kg, ζ = 5 Ns/m, k ₁ = -100N/m, k ₃ = 1 × 10⁶ N/m³, c = 5 Ns/m. The geometry illustrates the energy transfer path of SMR at the slow time scale, and the red asterisk denotes the unstable branch. The SMR regime generates from a cyclic interplay between 1:1 transient resonance capture relying on the stable branch and rapid isoenergetic jumps from the folding point to the landing point. The energy transfer loop of the SMR regime starts from the rapid jump from the singular point N ₁ on the lower stable branch to the landing point N _*u* on the stable upper branch. Then the slow flow of the system moves along the SIM until reaching another folding point N ₂ due to the transient resonance capture, and then, another isoenergetic jump back to the lower branch is triggered. After that, the slow flow continues moving along the SIM to the initial departure point, thereby completing one full cycle of energy transfer path and initiating the next. Therefore, the persistent bifurcations of the slow flow between two disjoint stable branches facilitate relaxation oscillations and yields SMR in the sustained forced nonlinear system, which is favorable for the purpose of efficient vibration mitigation. As illustrated in Figure 7b, the topological geometry of the SIM varies with the damping coefficient of the BNES. A key observation is that the unstable branch vanishes as damping increases, which prevents the triggering of the SMR regime. Since the SMR is crucial for effective energy transfer and vibration mitigation, its suppression under high damping is undesirable.

To determine the landing point following a rapid jump phenomenon, the bifurcation is considered to be isoenergetic. For an exit singular point N ₁ on the lower stable branch, the corresponding landing point N _*u* on the upper stable branch is identified by enforcing the invariance of the energy-related parameter |C| in Equation 35, and the value of N _*u* is expressed in the closed form

Z u = N u 2 = − p 2 + 2 p 2 2 − 3 p 1 p 3 3 p 1 (39)

Combining with Equation 34, the phase of the landing point on the upper stable branch of the SIM can also be computed aswhere Equations 39, 40 give the amplitude and corresponding phase of the landing point on the upper stable branch after the rapid jump, respectively. P Z = 4 ω 1 4 − 4 k 1 ω 1 2 − 3 k 3 Z / 8 c ω 1 3 and θ ₁ denotes the phase of the folding point on the lower stable branch. Similarly, a jump originating from an exit folding point N ₂ on the upper stable branch leads to a landing point N _*d* on the lower stable branch, and the amplitude and its corresponding phase of the landing point are computed explicitly by the same conservation law, which are expressed as

θ u = θ 1 + tan − 1 P Z 1 − P Z u 1 + P Z 1 P Z u (40)

Z d = N d 2 = − p 2 − 2 p 2 2 − 3 p 1 p 3 3 p 1 θ d = θ 2 + tan − 1 P Z 2 − P Z d 1 + P Z 2 P Z d (41)

where Equation 41 gives the amplitude and phase of the landing point on the lower stable branch. To validate the accuracy of the analytical study, numerical simulations were conducted by directly integrating Equation 10 with the application of Runge–Kutta method. As shown in Figure 8a, the response of the cantilever beam eventually settles into a steady-state periodic oscillation with zero initial conditions, and the system parameters are set as ε = 0.01kg, ζ = 5 Ns/m, k ₁ = -100N/m, k ₃ = 1 × 10⁶ N/m³, c = 5 Ns/m, σ = 0, f = 100 N. As shown in Figure 8b, the numerically obtained responses of the nonlinear system are superimposed on the analytically obtained SIM, the trajectory of the numerically obtained responses closely follows the geometry of the SIM and finally settles to the fixed point on the stable branch, confirming that the global system behavior stabilizes at a single equilibrium.

Despite the aforementioned steady-state periodic oscillations, the steady-state SMR is rather ubiquitous in the strongly nonlinear system due to the presence of unstable branches of the SIM. As illustrated in Figure 9a, the numerical simulations of steady-state responses demonstrate that the system indeed exhibits the SMR regime by direct integration of the governing Equation 10. As presented in Figure 9b, the numerical simulations are in good accordance with the analytical predictions on the SIM, where the SMR regime is characterized by sustained jump bifurcations of the slow flow between the two disjoint stable branches, which also confirms the analytical model of the slow-flow dynamics.

Following the asymptotic analysis of the dynamic evolution at the fast time scale, to further study the jumping bifurcations on the SIM and the corresponding relaxation oscillations of SMR scenario at the slow time scale, the dynamic behavior of the super-slow flow on the analytically obtained SIM is investigated in the limit of τ ₀→+ ∞. Accordingly, the variable φ ₂ is independent of the fast time scale τ0 but depends on the slow time scale τ ₁, i.e., φ ₂ = φ ₂ (τ ₁). Considering the O(ε) order of the subproblem of Equation 32, then the derivatives of the terms in the second equation of Equation 32 with respect to τ ₀ are zero, thereby yielding complex differential equation governing the dynamic evolution of super-slow flow on the SIM

∂ ∂ τ 1 c φ 2 2 + i ω 1 φ 2 2 − i k 1 φ 2 2 ω 1 − 3 i k 3 8 ω 1 3 φ 2 2 φ 2 + 3 2 σ + 1 ω 1 − i ζ k 3 16 ω 1 3 φ 2 2 φ 2 + i 2 c σ + c + ζ ω 1 4 + c ζ + 2 k 1 σ + k 1 4 − 2 ω 1 3 σ + i ζ k 1 4 ω 1 φ 2 − i 4 f ω 1 = 0 (42)

By means of straightforward algebraic manipulations, Equation 42 can be written in a more compact form aswhere asterisk represents the complex conjugate, and

c 2 + i ω 1 2 − k 1 2 ω 1 − 3 k 3 φ 2 2 4 ω 1 3 ∂ φ 2 ∂ τ 1 − 3 i k 3 φ 2 2 8 ω 1 3 ∂ φ 2 * ∂ τ 1 = G (43)

G = − 3 2 σ + 1 ω 1 − i ζ k 3 16 ω 1 3 φ 2 2 φ 2 − i 2 c σ + c + ζ ω 1 4 + c ζ + 2 k 1 σ + k 1 4 − 2 ω 1 3 σ + i ζ k 1 4 ω 1 φ 2 + i 4 f ω 1 (44)

To study the dynamic behavior of super-slow flow at the slow time scale, a reduced complex governing equation is derived by combining the complex conjugate of Equations 43, 44 and performing the necessary algebraic operations.

3 p 1 φ 2 4 + 2 p 2 φ 2 2 + p 3 ∂ φ 2 ∂ τ 1 = c 2 − i ω 1 2 − k 1 2 ω 1 − 3 k 3 φ 2 2 4 ω 1 3 G + 3 i k 3 φ 2 2 8 ω 1 3 G * (45)

To derive a set of real ordinary differential equations, the complex variable φ ₂ of Equation 45 is transformed as polar form φ 2 τ 1 = N τ 1 e i θ τ 1. The separation of real and imaginary parts of the resulting equation yields where prime denotes the differentiation with respect to the slow time scale τ ₁, and g(N) = 3p ₁ N ⁴+2p ₂ N ²+p ₃.

∂ N 2 ∂ τ 1 = g − 1 N − 9 ζ k 3 2 128 ω 1 6 N 5 − 1 8 ω 1 2 ( c + ζ ω 1 4 + ζ c 2 − 2 k 1 ω 1 2 + ζ k 1 2 ) N + 3 ζ k 3 ω 1 2 − k 1 16 ω 1 4 N 3 − 3 f k 3 32 ω 1 2 cos ⁡ θ N 2 + f 8 ω 1 2 ⁡ cos ⁡ θ + c ω 1 ⁡ sin ⁡ θ − k 1 ⁡ cos ⁡ θ ≡ f 1 N , θ ⇒ N ′ = f 1 N , θ (46)

∂ θ ∂ τ 1 = g − 1 N − 27 k 3 2 2 σ + 1 128 ω 1 5 N 4 − 3 k 3 − 8 σ − 1 ω 1 2 + 2 c ζ + 8 k 1 σ + 4 k 1 32 ω 1 3 N 2 + 9 f k 3 32 ω 1 2 N ⁡ sin ⁡ θ + 1 8 ω 1 − 2 ω 1 4 σ + − 2 c 2 + 4 k 1 σ − c 2 + k 1 ω 1 2 − k 1 2 2 σ + 1 + f 8 − ω 1 2 ⁡ sin ⁡ θ + c ω 1 ⁡ cos ⁡ θ + k 1 ⁡ sin ⁡ θ ≡ f 2 N , θ ⇒ θ ′ = f 2 N , θ (47)

Both the ordinary stable equilibrium points and folding singularities are determined by Equation 46 and Equation 47. The regular fixed points relying on the SIM should satisfy the following conditions

f 1 N , θ = 0 , f 2 N , θ = 0 , g N ≠ 0 (48)

On the other hand, the conditions of folding singularities of the super-slow flow should meet

f 1 N , θ = 0 , f 2 N , θ = 0 , g N = 0 (49)

Equations 48, 49 specify the conditions for the ordinary equilibrium points on the SIM and the folding singularities of the super-slow flow, respectively. Accordingly, the denominator g(N) = 0 determines the folding singularities that trigger the SMR regime and relaxation oscillations. The folding points correspond to the folding lines in the phase portrait, marking the jump points between stable SIM branches. To compute both types of equilibrium points, i.e., the ordinary stable equilibrium points and folding singular points, the numerators of Equation 46 and Equation 47 should simultaneously be equal to zero

f 1 N , θ = 0 , f 2 N , θ = 0 (50)

The two algebraic equations of Equation 50 can be written in the matrixwhere

α 11 N α 12 N α 21 N α 22 N cos ⁡ θ sin ⁡ θ = β 1 N β 2 N (51)

α 11 N = − f − 4 ω 1 4 + 3 N 2 k 3 + 4 k 1 ω 1 2 32 ω 1 2 , α 12 N = c f ω 1 8 , α 21 N = c f ω 1 8 N , α 22 N = f − 4 ω 1 4 + 9 N 2 k 3 + 4 k 1 ω 1 2 32 ω 1 2 N β 1 N = − 9 ζ k 3 2 128 ω 1 6 N 5 − − 24 ζ k 3 ω 1 4 + 24 ζ k 3 k 1 ω 1 2 128 ω 1 6 N 3 − 16 c + ζ ω 1 8 + 16 ζ c 2 − 2 k 1 ω 1 6 + 16 ζ k 1 2 ω 1 4 128 ω 1 6 N β 2 N = − 27 k 3 2 2 σ + 1 128 ω 1 5 N 4 − 3 k 3 − 8 σ − 1 ω 1 2 + 2 c ζ + 8 σ k 1 + 4 k 1 32 ω 1 3 N 2 + − 2 σ ω 1 4 + − 2 c 2 + 4 k 1 σ − c 2 + k 1 ω 1 2 − k 1 2 2 σ + 1 8 ω 1 (52)

where Equation 52 defines the coefficients in the algebra equation given by Equation 51. The phase of the singular points on the upper and lower folding lines is determined from Equation 51 and is expressed aswhere Θ 1 , 2 f s and Θ 3 , 4 f s of Equation 53 represent the phase pairs of the folding singularities on the corresponding lower and higher folding lines, respectively. Obviously, the existence of folding singularities should satisfy the following conditionswhere Equation 54 provides the existence condition for the folding singularities and k = 1,2. Specifically, the existence of the folding singularities on the lower stable branch should satisfy following condition

Θ 1 , 2 = Θ 1 , 2 f s = sin − 1 α 12 N 1 α 11 N 1 2 + α 12 N 1 2 ± ⁡ cos − 1 β 1 N 1 α 11 N 1 2 + α 12 N 1 2 Θ 3 , 4 = Θ 3 , 4 f s = sin − 1 α 12 N 2 α 11 N 2 2 + α 12 N 2 2 ± ⁡ cos − 1 β 1 N 2 α 11 N 2 2 + α 12 N 2 2 (53)

β 1 N k α 11 N k 2 + α 12 N k 2 ≤ 1 (54)

f ≥ f 1 c r i t = 9 ζ k 3 2 N 1 5 + 24 ζ k 3 ω 1 2 k 1 − ω 1 2 N 1 3 + 16 ω 1 4 c + ζ ω 1 4 + ζ c 2 − 2 k 1 ω 1 2 + ζ k 1 2 N 1 128 ω 1 6 − 4 ω 1 4 + 3 N 1 2 k 3 + 4 k 1 ω 1 2 / 32 ω 1 2 2 + c ω 1 / 8 2 (55)

Similarly, if the solvability condition is satisfied, the second pair of folding singularities will exist,

f ≥ f 2 c r i t = 9 ζ k 3 2 N 2 5 + 24 ζ k 3 ω 1 2 k 1 − ω 1 2 N 2 3 + 16 ω 1 4 c + ζ ω 1 4 + ζ c 2 − 2 k 1 ω 1 2 + ζ k 1 2 N 2 128 ω 1 6 − 4 ω 1 4 + 3 N 2 2 k 3 + 4 k 1 ω 1 2 / 32 ω 1 2 2 + c ω 1 / 8 2 (56)

Equations 55, 56 define the critical forcing amplitudes for the existence of folding singularities on the lower and upper folding lines, respectively. The phase trajectories of the super-slow flow relying on the SIM are illustrated in Figure 10 for different external forcing amplitudes. Referring to the topological geometry of the SIM, each of the phase portrait is divided into three regions by the two folding lines N=N ₁ and N=N ₂, specifically, the lower, middle, and upper region correspond to lower stable branch, unstable branch and higher stable branch respectively. There is no trajectory in the middle blank region due to the jump bifurcation, and the arrows on the stable phase trajectory denote the directions of the corresponding super-slow flow at the slow time scale.

The system parameters are assigned as k ₁ = -100N/m, k ₃ = 1 × 10⁶ N/m³, c = 5 Ns/m. Accordingly, the critical forcing amplitudes required for the existence of folding singularities on the lower and higher folding lines, which trigger the SMR regime, are determined to be f _1*crit* = 34.7 N and f _2*crit* = 331.1N, respectively. As illustrated in Figure 10a, when the external forcing amplitude f = 30 N is lower than the first critical threshold, despite the possibility of a single jump from the upper to the lower branch for specific initial conditions, the super-slow flow is always attracted to a stable equilibrium point, which prevents the realization of the SMR regime. As the amplitude of the external harmonic excitation further increases to f = 200 N that is greater than the first critical amplitude, there exist two singular points on the lower folding line (cf. Figure 10b), and the rapid jump from the lower folding line to the upper folding line is possible. On the other hand, the forcing amplitude also exceeds the critical threshold of SN bifurcation, then the resulting cyclic rapid jump bifurcations between stable branches of the SIM become possible, resulting the relaxation oscillations.

Further increasing the forcing amplitude to f = 350N, as shown in Figure 10c, the amplitude of harmonic excitation exceeds both the second critical threshold and the Hopf bifurcation point. As a result, the sustained jump bifurcations of the slow flow between the two stable branches of the slow invariant manifold yields the instability of the steady-state strongly modulated response, which facilitates efficient and irreversible TET from the cantilever beam to the nonlinear energy sink. Specifically, following an isoenergetic jump bifurcation, the super-slow flow jumps to a landing point N _*u* on the stable upper branch. Then the super-slow flow evolves along this stable branch until reaching the upper folding line N=N ₂, at which point it jumps back to the lower branch and lands at the point N _*d* on the lower stable branch. A subsequent super-slow transition returns the trajectory to the original departure point N ₁, thereby completing one SMR cycle and initiating the next cyclic motion. For a further increase in the amplitude of the external harmonic excitation up to f = 500N, the system response transcends the Hopf and SN bifurcation regimes. In this scenario, as shown in Figure 10d, the phase trajectories that jump from the lower to the upper region are predominantly captured by a stable equilibrium. Consequently, the steady-state response converges to a periodic orbit as the forced dynamics are attracted to this equilibrium.

In summary, by the application of slow-fast decomposition and complexification-averaging method, the rich and complex dynamics of the nonlinear continuous system with BNES are justified. The SMR regime generates from the cyclic motion of the super-slow flow between the stable SIM branches, and the strongly modulated responses enable efficient and irreversible TET from the cantilever beam to the nonlinear energy sink. To achieve improved vibration reduction, it is necessary to broaden the frequency band of the SMR, which provides a guideline for the design of the BNES.

4 Vibration suppression performance of the BNES

This section investigates the vibration suppression efficacy of the BNES on the directly forced cantilever beam. In consistency to the preceding analysis, the fundamental resonance regimes are examined, accordingly, the oscillating frequencies of the cantilever beam and the BNES are identical to the frequency of external harmonic excitation. The objective of this section is to numerically verify the vibration mitigation performance of the BNES by parametric analysis, which plays a critical role in the vibration mitigation performance of the NES.

It is noted that multiple fundamental resonance regimes may coexist within the frequency range of interest for the essential nonlinear system, for example, periodic, quasi-periodic or chaotic responses. Accordingly, since the amplitude of non-periodic response is time dependent, the energy spectrum of the modal displacement is applied to evaluate the capacity of vibration attenuation by the proposed NES. The energy spectrum is computed by averaging the kinetic and potential energy of the beam across a specific time period when the harmonically forced nonlinear system reaches steady state, which is expressed as

E a = 1 2 x ˙ 1 2 + 1 2 ω 1 2 x 1 2 t (57)

Equation 57 defines the time-averaged energy used to evaluate the vibration suppression performance of the BNES. Figure 11a illustrates the effect of the cubic stiffness of the BNES on the energy spectrum of steady-state cantilever beam, and a smaller value of the energy spectrum corresponds to higher vibration mitigation efficiency. The other system parameters are fixed as ε = 0.01kg, ζ = 5 Ns/m, k ₁ = -100N/m, c = 8 Ns/m and f = 50 N. The total averaging energy of the cantilever beam exhibits a monotonic reduction with increasing cubic stiffness of the BNES. Figure 11b presents the effect of the damping coefficient of the BNES on the energy spectrum of steady-state cantilever beam, and the other system parameters are set as ε = 0.01kg, ζ = 5 Ns/m, k ₁ = -100N/m, k ₃ = 1 × 10⁶ N/m³, and f = 50 N. When the damping coefficient reaches 9 Ns/m, the system response becomes periodic, and the amplitude decreases with further increases in damping. In contrast, at low damping ratios, unstable steady-state strongly modulated responses emerge, indicating efficient energy transfer and dissipation within the BNES. This is reflected in a significantly lower energy spectrum compared to the high-damping case at the fundamental resonance. Thus, the occurrence of SMRs substantially enhances the vibration mitigation performance of the BNES.

Figure 11c illustrates the influences of the beam damping coefficient on the energy spectrum of steady-state cantilever beam, and the other system parameters are set as ε = 0.01kg, c = 8 Ns/m, k ₁ = −100N/m, k ₃ = 1 × 10⁶N/m³, and f = 50 N. The variations of damping coefficient of the beam shift and narrow the frequency interval of the strongly modulated response. Consequently, the higher damping ratio of the beam may attenuate nonlinear modulation and targeted energy transfer in the continuous system. Figure 11d presents the effects of the excitation amplitude on the energy spectrum of steady-state cantilever beam and the other system parameters are set as ε = 0.01kg, c = 8 Ns/m, ζ = 5 Ns/m, k ₁ = −100N/m, and k ₃ = 1 × 10⁶N/m³. As the excitation amplitude increases, the beam energy level also rises and the frequency range of SMR broadens around the fundamental resonance, which indicates that the stronger external harmonic excitation extends the range of effective targeted energy transfer.

5 Conclusion

This study systematically investigates the complex nonlinear behavior and vibration suppression performance of a flexible cantilever beam coupled with a bistable nonlinear energy sink under harmonic excitations through numerical simulations and theoretical analysis. The results demonstrate that the proposed BNES can achieve efficient and robust vibration suppression in typical excitation environments through the targeted energy transfer mechanism, revealing its significant potential for broadband vibration control in engineering structures.

Under sustained harmonic excitation, the continuous system presents fruitful nonlinear dynamical responses. When the excitation frequency is near the first natural frequency of the beam, the system exhibits various steady-state responses, including periodic and strongly modulated responses. The energy transfer mechanism of the SMR regime is studied analytically by the complexification-averaging method and multiple scales method. The analysis indicates that SMR regime generates from sustained jump bifurcations of the slow flow between two disjoint stable branches of the slow invariant manifold, resulting in relaxation oscillations. This cyclic capture-jump dynamic behavior provides an ideal path for steady-state, irreversible targeted energy transfer from the continuous structure to the BNES.

Furthermore, in-depth bifurcation analysis and parametric studies reveal that the performance of the BNES is highly sensitive to its parameter settings. Saddle-node and Hopf bifurcations collectively determine the existence and stability regions of the system’s steady-state responses. Parametric analysis further indicates that within a certain damping range, lower damping coefficients result in the SMR, thereby achieving optimal vibration suppression performance, while increasing the cubic stiffness generally helps reduce the vibration energy level of the primary structure. This provides clear theoretical guidance for the engineering design and parameter tuning of the BNES.

In summary, this work not only theoretically clarifies the complex dynamic mechanisms in the coupled system of a bistable nonlinear energy sink and a continuum structure, particularly the conditions for the generation of strongly modulated responses and the energy transfer path, but also verifies its excellent vibration suppression performance under harmonic excitations through detailed numerical simulations and theoretical analysis. This research lays a solid theoretical foundation for applying BNES to the broadband passive control of practical engineering structures.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmech.2026.1920583/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Roman Protasov, Slovenska technicka univerzita v Bratislave Strojnicka fakulta, Slovakia

Reviewed by: Muhammad Nadeem, Qujing Normal University, China

Omid Koochakianfard, Tarbiat Modares University, Iran