Optimizing stability and characteristics of a vibrating rigid body pendulum with energy harvesting device

Abstract

The focus of this paper is on reducing vibrations and capturing energy from a three degrees-of-freedom (3DOF) dynamical system. The system consists of a damped spring-pendulum with an elliptical pivot path at a constant angular velocity and a connected rigid body experiencing external harmonic forces and moments. An independent electromagnetic energy harvesting system has been integrated into the structure, utilizing the motion of a magnet within a coil. The optimization efforts target both energy harvesting (EH) and vibration reduction capabilities. The second kind of Lagrange’s equations are used to formulate the governing kinematic equations, which are then solved asymptotically using the multiple-scales technique (MST), producing novel and precise results up to the third approximation. Various resonance cases are identified, with three examined concurrently. The time evolution of the solutions, as well as the modified phases and amplitudes, are analyzed through graphical representations, taking into account the influential parameters of the system. Graphical plots of phase portraits successfully depicted the system’s stable dynamics. The power output and current time profiles of the electromagnetic device are presented to demonstrate the impact of various parameters on system behavior. Analysis of stability and instability areas reveals that the system exhibits stable performance across a wide range of parameters. This model is gaining significance in contemporary applications, especially in control sensors used across various sectors such as industry, construction, automotive, transportation, and infrastructure.

Introduction

Over the past few decades, a great deal of study has been done on reducing damaging vibrations in a variety of architectural structures, including robotic limbs, high-rise buildings, bridges, spires, and turbines. ¹ Vibration energy is an intriguing new energy source that is both sustainable and renewable. The process involves converting vibrational kinetic energy into unstable electrical power by scavenging, or the use of large amplitude vibrations for energy harvesting. Vibrations, mechanical strain, and stress in mechanical systems can all be effectively captured by this process. Furthermore, EH has shown to be a perfect option for supplying power to tiny gadgets such as wireless sensors.^2,3

The limitation of reliable power sources has spurred the growing adoption of mobile and cordless technologies into daily routines. Depending on bulky batteries that require frequent replacement for certain devices poses significant challenges. Hence, integrating EH as the primary power source for small-scale systems emerges as a practical solution. Leveraging sources such as sunlight, vibrations, and heat enables the swapping of temporary energy reservoirs with consistently charged energy storage components for a continuous power supply.^4–6 Energy extraction from these origins employs electromagnetic conversion, electrostatic, or piezoelectric methods. However, the oscillations produced usually function within lower frequency ranges, resulting in comparatively modest power outputs. As a result, the EH from these mechanisms typically falls within the micro to milliwatt range.^7–13

Electromagnetic energy harvesters harness the principles of Faraday’s induction law, which dictates that any alteration in the magnetic field surrounding a wire coil will induce an electromotive force within the coil. This phenomenon occurs when a magnet moves relative to the coil. ¹⁴ Electrostatic transmutation, on the other hand, employs capacitors to convert mechanical energy into an electric field. By utilizing a variable capacitor, it functions as a current source capable of powering an electric circuit. ¹⁵ Meanwhile, piezoelectric materials generate an electric field when subjected to pressure, leading to voltage generation when deformed under stress. ¹⁶

In Ref. [17], the MST is employed to analyze an electromagnetic elastic-pendulum harvester model, focusing specifically on a two-to-one ratio of vibratory energy harvesters. In Ref. [18], the authors explore techniques for minimizing vibrations and capturing energy from a spring-pendulum system that moves along an elliptical trajectory. The pendulum’s structure is modified to integrate a separate electromagnetic harvesting system. In Ref. [19], the nonlinear movement of a connected damped pendulum to electromagnetic and piezoelectric devices is explored, employing approximate solutions and graphical representations to assess EH performance.

In Ref. [20], a novel strategy is introduced to enhance the efficiency of energy harvesting under low excitation conditions. This method involves incorporating a piezoelectric element attached to an inverted beam and pendulum system. Validation of this approach is achieved through an extensive finite element model, yielding valuable insights and simulations. Additionally, ²¹ investigates the application of a pendulum-based harvester to tap into the energy generated by vibrations from high-speed train motion. Expanding on this theme, Ref. [22] taps into the behavior of the base of a double pendulum under excitation to generate electric energy on a small level. Notably, there exists a qualitative alignment between numerical and experimental outcomes. Meanwhile, Ref. [23] introduces a magnetic rolling pendulum tailored for EH, showcasing remarkable energy production in various resonance modes. This system utilizes a magnetically levitated permanent magnet, effectively addressing issues related to mechanical damping.

To simulate EH within a particular excitation range, researchers create a nonlinear dielectric elastomer pendulum, ²⁴ and their results are rigorously analyzed by experimentation. Furthermore, Ref. [25] presents a simple double pendulum oscillator device that may be easily adjusted to change its inherent frequency and effectively transfer wave dynamics into electrical energy. Moreover, Ref. [26] performs a comprehensive analysis of experimental and theoretical aspects related to the dynamic motion of several pendulum systems in a magnetic field. Whereas Ref. [27] and Ref. [28] explore the stability of dynamical systems combined with EH mechanisms, demonstrating how vibrations can be transformed into useful electrical energy resources. Continuing with Ref. [29], the authors presented a meticulous numerical exploration of an auto-parametric pendulum-absorber paired with a composite EH arrangement. Their analysis extends to evaluating the electrical efficacy within the system and its adherence to various electrical standards. Meanwhile, in Ref. [30], an EH scheme centered on a dynamic pendulum-absorber coupled with an electromagnetic harvester is delved.

A novel method for simultaneous EH and vibration attenuation is presented in Ref. [31], where an electromagnetic harvester is linked to an auto-parametric pendulum that acts as the absorber. Examining this new hybrid harvester in detail demonstrates that the effectiveness of vibration mitigation is unaffected. Especially, the energy extracted from the chaotic and rotating motion of the pendulum turns out to be significant. In Ref. [32], authors unveil a pendulum configuration incorporating two separate electromagnetic harvesters, one of which integrates a direct current (DC) electric motor onto the pendulum’s axis. The obtained outcomes reveal a significant enhancement in EH capability, albeit with potential trade-offs in vibration suppression performance. Transitioning to Ref. [33], the vibration-damping device evolves into an EH apparatus and vibration-dampening system, adept at counteracting undesired vibrations and furnishing voltage for low-power applications like wireless sensors. In Ref. [34], an exhaustive numerical analysis is presented for a semi-active suspension system equipped with an absorber engineered for both vibration reduction in the primary structure and energy generation. This absorber, a pendulum-tuned mass absorber, is controlled using a magnetic levitation harvesting mechanism, illustrating a multifaceted solution to EH and vibration control.

In Ref. [35], the authors explored the synergy between EH and vibration control by employing a specially designed pendulum-absorber in conjunction with nonlinear vibration absorbers optimized for energy extraction. The obtained results highlighted on the effectiveness of integrating an energy recovery device in attenuating vibrations, especially at lower frequencies. Meanwhile, in Ref. [36], a novel configuration of a pendulum-absorber harvester is introduced, leveraging magnets to enhance both vibrations damping and EH capabilities. By combining theoretical insights with empirical investigations, this innovative setup outperforms conventional devices in terms of both vibration attenuation and energy extraction efficiency. Turning to Ref. [37], the harmonic balance method took center stage in exploring the response and stability of mechanical and electrical systems for EH. Here, a harvesting system comprising fixed magnets at one end of a tube, with a movable third magnet, undergoes thorough scrutiny regarding its parameters' impact on responsiveness and energy retrieval, particularly in proximity to the main resonance. The energy harvesting resources for railways, various harvesting methods with their pros and cons, and their applications in the industry were explored in Ref. [38]. The authors introduced and compared different models and prototypes of energy harvesters. The potential applications and future challenges are also discussed to highlight key principles for designing railway energy harvesters.

Investigating nonlinear vibrations in 3DOF systems enables engineers to better predict and control the behavior of structures and machinery under real-world conditions. In such a system, there are three independent generalized coordinates. They typically correspond to three different modes of vibration or motion. ³⁹

In linear systems, the relationship between forces and displacements follows Hooke’s Law and the equations of motion are linear differential equations. The system’s response can be superposed, meaning that if the system is subjected to multiple forces, the resulting motion is the sum of the responses to each force. Whereas in nonlinear systems the relationship between forces and displacements is not linear. This can result from nonlinearities in the springs, and dampers, or due to interactions between different parts of the system. Nonlinearities introduce complexities such as amplitude-dependent frequencies, bifurcations, and chaotic behavior. ⁴⁰ The equations governing this system are typically derived from Lagrangian mechanics or Newton’s laws. Nonlinear systems can have multiple equilibrium points. Stability analysis involves studying how small perturbations affect these equilibrium points. They can exhibit chaotic behavior, where small changes in initial conditions lead to drastically different outcomes. Tools such as Lyapunov exponents and Poincaré maps are used to analyze chaos. ⁴¹ Nonlinear vibrations are important in various engineering applications, including mechanical systems with large deflections, aerospace structures, and vehicle suspension systems. Understanding these vibrations helps in designing systems that can withstand or avoid undesirable resonances and ensure safety and performance.

Our work meticulously outlines the mathematical representation and operational principles, employing the MST to derive approximate solutions aligned with a novel 3DOF vibrating dynamical system according to its generalized coordinates. This system involves a spring-pendulum system where the pivot point travels an elliptical route at a constant angular speed, coupled with a rigid body connected to the spring. It is influenced by an external harmonic force and moments. Additionally, it features an independent electromagnetic EH device, which utilizes the movement of a magnet within a coil to enhance its functionality. Furthermore, equations of modulation are formulated to address resonance scenarios, while dynamic response curves near resonance regions inform assessments of system stability. Historical variations of solutions over time also undergo meticulous scrutiny, solidifying the device’s potential for transformative applications. Phase portraits serve as visual aids in depicting the system’s enduring behavior. Through plotted output power and current time histories of the electromagnetic device, the influence of various parameters on the system’s dynamics becomes evident. Detailed scrutiny of adjusted phases and amplitudes showcases the impact of the body’s physical parameters at each scenario. Moreover, phase planes are meticulously presented to validate the existence of stable fixed points (FPs) under steady-state conditions. In-depth nonlinear stability analysis employing the Routh-Hurwitz criteria (RHC) is undertaken and thoroughly deliberated upon. This rigorous examination sheds light on the system’s robustness and resilience in the face of dynamic perturbations, enriching our understanding of its intricate behavior. Through this analysis, regions of stability and instability are explored, indicating that the system maintains stability over a broad spectrum of parameters. This examined system holds growing significance in contemporary contexts, especially in applications concerning sensor control across diverse sectors such as industry, construction, automotive, transportation, and infrastructure.

Choosing an elliptical path for the pivot point of a pendulum in a vibrating system can offer diverse advantages and compelling dynamic behavior, such as:

This trajectory can more precisely mimic specific real-world movements than simple linear or circular paths, which is particularly beneficial in scenarios where the pivot point’s motion isn’t restricted to a single axis or a basic geometric pattern. This motion induces two-dimensional harmonic excitation, which can produce more complex resonance phenomena and mode interactions, useful for investigating the system’s response in more realistic conditions. ⁴²

Moreover, by adjusting the parameters of the elliptical path (such as the major and minor axes lengths and the frequency of motion), the energy input into the system can be precisely controlled. This allows for fine-tuning the excitation to achieve desired vibrational characteristics. The nonlinearities of the vibrating systems can be introduced through this motion, offering extensive opportunities to study nonlinear behavior, stability, bifurcations, and chaos in forced vibrating systems. In many practical engineering systems, the pivot point of a pendulum or similar mechanism may follow an elliptical trajectory due to design constraints or operational conditions. Studying such paths can provide insights that are directly applicable to real-world systems, such as robotics, aerospace, and mechanical engineering. In certain instances, an elliptical path can simulate the effects of gravitational forces in a varying field, with the changing orientation of the path representing the gravitational force components in different directions. ⁴³

Thus, employing an elliptical path can offer a more comprehensive understanding of the system’s dynamics and aid in designing more effective control strategies for practical applications.

Characterization of the problem

Let’s consider the nonlinear vibrations of a 3DOF dynamical system. Here, we examine the motion of a rigid body constrained to move within a plane, hanging from a point $\mathrm{A}$ to one end of a weightless spring possessing a specific stiffness value $k$ . The spring’s other end is fastened to a mobile location $N$ , tracing an elliptical path characterized by major and minor axes denoted by $2b$ and $2a$ , respectively. Considering a point $Q$ depicted in Figure (1), which serves as the projection of the moving point $N$ onto an auxiliary circle $b$ and maintaining a steady angular velocity $\Omega$ . Hence, the coordinates of $N$ can be formulated after a specific duration of time as follows:

${\mathit{X}}_{\mathit{N}}=\mathit{a}\hspace{0.17em}\mathit{cos}\hspace{0.17em}\mathbf{\Omega }\mathit{t},{\mathit{Y}}_{\mathit{N}}=\mathit{b}\hspace{0.17em}\mathit{sin}\hspace{0.17em}\mathbf{\Omega }\mathit{t}.$

(1)

Figure 1.

A schematic diagram of the examined dynamical system.

The rigid body’s characteristics are defined by its mass $m$ , its moment of inertia $U_A$ about the axis passing through the point $\mathrm{A}$ and perpendicular to the plane’s motion. Point $\mathrm{C}$ indicates the body’s center of mass, and the distance $p=\text{AC}$ representing its eccentricity. The spring’s damping coefficient is denoted by $C_1$ , while $C_2$ and $C_3$ are the rotatory damping coefficients at the points of rotation $N$ and $\mathrm{A}$ , respectively. Considering $\mathcal{l}={\mathcal{l}}_0+Z_{st}$ is the length of the spring at the static equilibrium position after time $t$ , where ${\mathcal{l}}_0$ is the spring’s natural length and $Z_{st}=\frac{mg}{k}$ represents the static elongation, the acceleration due to gravity is denoted as $g$ . All play important roles in the system.

Moreover, the system under study is connected to an electromagnetic device that enables the extraction of electrical power. This is achieved by harnessing the current $\mathrm{I}$ generated from the magnet’s oscillations within a coil. The electric component of the harvester is made up of various elements, including the load resistor $R_L$ , coil’s inductance $L_{ind}$ , the resistance of a coil $R_C$ , and coil length $L_{coil}$ , as shown in Figure (1).

The motion is considered under the influence of harmonic external force $F\left(t\right)=F_0\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{1}t$ , and moments $M_{\phi }\left(t\right)=M_1^{*}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{2}t$ and $M_{\psi }\left(t\right)=M_2^{*}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{3}t$ . Here, $F_0,M_s^{*}\left(s=1,2\right)$ and ${\Omega }_j\left(j=1,2,3\right)$ are the amplitudes and frequencies of the external force $F\left(t\right)$ and the moments $M_{\phi }\left(t\right),M_{\psi }\left(t\right)$ , respectively.

To derive the equations of motion (EOM), we scrutinize the variables representing the spring’s elongation $Z$ , as well as the angular displacements $\phi$ and $\psi$ at $N$ and $\mathrm{A}$ , respectively.

Thus, the system’s kinetic and potential energies are denoted as $T$ and $V$ , along with dissipation’s function $D$ may be expressed as follows

$\mathit{T}=\frac{\mathbf{1}}{\mathbf{2}}\mathit{m}{\mathit{\Omega }}^{\mathbf{2}}\left[{\mathit{a}}^{\mathbf{2}}\hspace{0.17em}{\mathit{sin}}^{\mathbf{2}}\hspace{0.17em}\mathit{\Omega }\mathit{t}+{\mathit{b}}^{\mathbf{2}}\hspace{0.17em}{\mathit{cos}}^{\mathbf{2}}\hspace{0.17em}\mathit{\Omega }\mathit{t}\right]+\mathit{m}\mathit{\Omega }(\mathit{b}\hspace{0.17em}\mathit{sin}\hspace{0.17em}\mathit{\phi }\hspace{0.17em}\mathit{cos}\hspace{0.17em}\mathit{\Omega }\mathit{t}-\mathit{a}\hspace{0.17em}\mathit{cos}\hspace{0.17em}\mathit{\phi }$

$\times \mathit{sin}\hspace{0.17em}\mathit{\Omega }\mathit{t})\dot{\mathit{Z}}+\mathit{m}\mathit{\Omega }\left({\ell }_{\mathbf{0}}+\mathit{Z}\right)\left(\mathit{b}\hspace{0.17em}\mathit{cos}\hspace{0.17em}\mathit{\phi }\hspace{0.17em}\mathit{cos}\hspace{0.17em}\mathit{\Omega }\mathit{t}+\mathit{a}\hspace{0.17em}\mathit{sin}\hspace{0.17em}\mathit{\phi }\hspace{0.17em}\mathit{sin}\hspace{0.17em}\mathit{\Omega }\mathit{t}\right)\dot{\mathit{\phi }}$

$+\mathit{m}\mathit{p}\mathit{\Omega }\left(\mathit{b}\hspace{0.17em}\mathit{cos}\hspace{0.17em}\mathit{\psi }\hspace{0.17em}\mathit{cos}\hspace{0.17em}\mathit{\Omega }\mathit{t}+\mathit{a}\hspace{0.17em}\mathit{sin}\hspace{0.17em}\mathit{\psi }\hspace{0.17em}\mathit{sin}\hspace{0.17em}\mathit{\Omega }\mathit{t}\right)\dot{\mathit{\psi }}+\mathit{m}\mathit{p}\dot{\mathit{Z}}\dot{\mathit{\psi }}\mathit{sin}\left(\mathit{\phi }-\mathit{\psi }\right)$

$+\mathit{m}\mathit{p}\left({\ell }_{\mathbf{0}}+\mathit{Z}\right)\dot{\mathit{\phi }}\dot{\mathit{\psi }}\mathit{cos}\left(\mathit{\phi }-\mathit{\psi }\right)+\frac{\mathbf{1}}{\mathbf{2}}\mathit{m}{\dot{\mathit{Z}}}^{\mathbf{2}}+\frac{\mathbf{1}}{\mathbf{2}}\mathit{m}{\left({\ell }_{\mathbf{0}}+\mathit{Z}\right)}^{\mathbf{2}}{\dot{\mathit{\phi }}}^{\mathbf{2}}$

$+\frac{\mathbf{1}}{\mathbf{2}}{\dot{\mathit{\psi }}}^{\mathbf{2}}\left(\mathit{m}{\mathit{p}}^{\mathbf{2}}+{\mathit{U}}_{\mathit{A}}\right),\mathit{D}=\frac{\mathbf{1}}{\mathbf{2}}\left({\mathit{C}}_{\mathbf{1}}{\dot{\mathit{Z}}}^{\mathbf{2}}+{\mathit{C}}_{\mathbf{2}}{\dot{\mathit{\phi }}}^{\mathbf{2}}+{\mathit{C}}_{\mathbf{3}}{\dot{\mathit{\psi }}}^{\mathbf{2}}\right),$

$\mathit{V}=-\mathit{m}\mathit{g}\left[\mathit{a}\hspace{0.17em}\mathit{cos}\hspace{0.17em}\mathit{\Omega }\mathit{t}+\left({\ell }_{\mathbf{0}}+\mathit{Z}\right)\mathit{cos}\hspace{0.17em}\mathit{\phi }+\mathit{p}\hspace{0.17em}\mathit{cos}\hspace{0.17em}\mathit{\psi }\right]+\frac{\mathbf{1}}{\mathbf{2}}\mathit{k}{\mathit{Z}}^{\mathbf{2}},$

(2)

where the dots represent the rates of change over time $t$ .

Referring to equation (2), obtaining the Lagrangian $S=T-V$ is straightforward. Subsequently, we can calculate the EOM according to the second kind of Lagrange’s equations ⁴⁴ as outlined below

$\frac{d}{dt}\left(\frac{\partial S}{\partial \dot{Z}}\right)-\left(\frac{\partial S}{\partial Z}\right)=Q_Z-\left(\frac{\partial D}{\partial \dot{Z}}\right)$

$\frac{d}{dt}\left(\frac{\partial S}{\partial \dot{\phi }}\right)-\left(\frac{\partial S}{\partial \phi }\right)=Q_{\phi }-\left(\frac{\partial D}{\partial \dot{\phi }}\right),$

$\frac{d}{dt}\left(\frac{\partial S}{\partial \dot{\psi }}\right)-\left(\frac{\partial S}{\partial \psi }\right)=Q_{\psi }-\left(\frac{\partial D}{\partial \dot{\psi }}\right),$

(3)

where $Q_Z,Q_{\phi },$ and $Q_{\psi }$ are the generalized forces, which take the forms

$Q_Z=F\left(t\right)=F_0\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{1}t$

$Q_{\phi }=M_{\phi }\left(t\right)=M_1^{*}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{2}t,$

$Q_{\psi }=M_{\psi }\left(t\right)=M_2^{*}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{3}t.$

(4)

The substitution of (2) and (4) into (3), produces EOM in the form

$m\ddot{Z}+C_1\dot{Z}+kZ-mg\hspace{0.17em}\mathit{cos}\hspace{0.17em}\phi -m\left({\ell }_0+Z\right){\dot{\phi }}^2-m{\Omega }^2(b\hspace{0.17em}\mathit{sin}\hspace{0.17em}\phi \hspace{0.17em}\mathit{sin}\hspace{0.17em}\Omega t$

$+a\hspace{0.17em}\mathit{cos}\hspace{0.17em}\phi \hspace{0.17em}\mathit{cos}\hspace{0.17em}\Omega t)+mp\ddot{\psi }\mathit{sin}\left(\phi -\psi \right)-mp{\dot{\psi }}^2\hspace{0.17em}\mathit{cos}\left(\phi -\psi \right)+B{L}_{coil}\mathrm{{\rm I}}=F_0\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{1}t,$

(5)

$m{\left({\ell }_0+Z\right)}^2\ddot{\phi }+C_2\dot{\phi }+mg\left({\ell }_0+Z\right)\mathit{sin}\hspace{0.17em}\phi +2m\left({\ell }_0+Z\right)\dot{Z}\dot{\phi }+mp\left({\mathcal{l}}_0+Z\right)\ddot{\psi }$

$\times \mathit{cos}\left(\phi -\psi \right)+mp\left({\ell }_0+Z\right){\dot{\psi }}^2\hspace{0.17em}\mathit{sin}\left(\phi -\psi \right)+m{\Omega }^2\left({\ell }_0+Z\right)(a\hspace{0.17em}\mathit{sin}\hspace{0.17em}\phi \hspace{0.17em}\mathit{cos}\hspace{0.17em}\Omega t$

$-b\hspace{0.17em}\mathit{cos}\hspace{0.17em}\phi \hspace{0.17em}\mathit{sin}\hspace{0.17em}\Omega t)=M_1^{*}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{2}t,$

(6)

$\ddot{\psi }\left(m{p}^2+U_A\right)+C_3\dot{\psi }+mgp\hspace{0.17em}\mathit{sin}\hspace{0.17em}\psi +\mathit{sin}\left(\phi -\psi \right)\left[mp\ddot{Z}-mp\left({\ell }_0+Z\right){\dot{\phi }}^2\right]$

$+\mathit{cos}\left(\phi -\psi \right)\left[mp\left({\ell }_0+Z\right)\ddot{\phi }+2mp\dot{\phi }\dot{Z}\right]+mp{\Omega }^2(a\hspace{0.17em}\mathit{sin}\hspace{0.17em}\psi \hspace{0.17em}\mathit{cos}\hspace{0.17em}\Omega t$

$-b\hspace{0.17em}\mathit{cos}\hspace{0.17em}\psi \hspace{0.17em}\mathit{sin}\hspace{0.17em}\Omega t)=M_2^{*}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{3}t,$

(7)

By utilizing Kirchhoff’s voltage law and Faraday’s law ¹⁸ of electromagnetic induction on the circuit shown in Figure (1), we are able to express the following equation

$L_{ind}\dot{\mathrm{{\rm I}}}+R\mathrm{{\rm I}}-B{L}_{coil}\dot{Z}=0,$

(8)

where the sum of the coil’s and load’s resistances $R_C$ and $R_L$ is represented by $R$ and $B$ is the magnetic flux density of the magnet.

The below dimensionless parameters may be used to convert the EOM into their dimensionless form

$w_1^2=\frac{k}{m},w_2^2=\frac{g}{\mathcal{l}},w_3^2=\frac{mpg}{U_A},{\omega }_1=\frac{w_1}{w_2},{\omega }_2=\frac{w_2}{w_3},\varpi =\frac{\Omega }{w_2},$

$s=\frac{p}{\mathcal{l}},p_1=\frac{{\Omega }_1}{w_2},p_2=\frac{{\Omega }_2}{w_2},p_3=\frac{{\Omega }_3}{w_2},\gamma =\frac{{{\rm B}}^2{L}_{coil}^2}{mg{L}_{ind}}\mathcal{l},\alpha =\frac{R}{L_{ind}{w}_2},$

$c_1=\frac{C_1}{m{w}_2},c_2=\frac{C_2}{m{w}_2{\mathcal{l}}^2},c_3=\frac{C_3}{m{w}_{2}p\mathcal{l}},f_1=\frac{F_0}{mg},f_2=\frac{M_1^{*}}{m\mathcal{l}g},$

$f_3=\frac{M_2^{*}}{mpg},z=\frac{Z-Z_{st}}{\mathcal{l}},x=\frac{L_{ind}}{B{L}_{coil}\mathcal{l}}\mathrm{{\rm I}},\tau =w_{2}t.$

(9)

Therefore, the substitution of (9) into (5)-(8) yields the following forms of the EOM

$z^{″}+c_1{z}^{\prime }+{\omega }_1^{2}z+\left(1-\mathit{cos}\hspace{0.17em}\phi \right)-\left(1+z\right){\phi }^{\prime 2}-s{\psi }^{\prime 2}\hspace{0.17em}\mathit{cos}\left(\phi -\psi \right)+s{\psi }^{″}\hspace{0.17em}\mathit{sin}\left(\phi -\psi \right)$

$-\frac{{\varpi }^2}{\mathcal{l}}\left(b\hspace{0.17em}\mathit{sin}\hspace{0.17em}\phi \hspace{0.17em}\mathit{sin}\hspace{0.17em}\varpi \tau +a\hspace{0.17em}\mathit{cos}\hspace{0.17em}\phi \hspace{0.17em}\mathit{cos}\hspace{0.17em}\varpi \tau \right)+\gamma x=f_1\hspace{0.17em}\mathit{cos}\hspace{0.17em}{p}_1\tau ,$

(10)

${\left(1+z\right)}^2{\phi }^{″}+c_2{\phi }^{\prime }+\left(1+z\right)\mathit{sin}\hspace{0.17em}\phi +2\left(1+z\right)z^{\prime }{\phi }^{\prime }+s\left(1+z\right){\psi }^{″}\hspace{0.17em}\mathit{cos}\left(\phi -\psi \right)$

$+s\left(1+z\right){\psi }^{\prime 2}\hspace{0.17em}\mathit{sin}\left(\phi -\psi \right)+\frac{{\varpi }^2}{\mathcal{l}}(1+z)(a\hspace{0.17em}\mathit{sin}\hspace{0.17em}\phi \hspace{0.17em}\mathit{cos}\hspace{0.17em}\varpi \tau -b\hspace{0.17em}\mathit{cos}\hspace{0.17em}\phi \hspace{0.17em}\mathit{sin}\hspace{0.17em}\varpi \tau )=f_2\hspace{0.17em}\mathit{cos}\hspace{0.17em}{p}_2\tau ,$

(11)

$\left(s+{\omega }_2^2\right){\psi }^{″}+c_3{\psi }^{\prime }+\mathit{sin}\hspace{0.17em}\psi +\mathit{sin}\left(\phi -\psi \right)\left[z^{″}-\left(1+z\right){\phi }^{\prime 2}\right]+\mathit{cos}(\phi -\psi )[\left(1+z\right){\phi }^{″}$

$+2z^{\prime }{\phi }^{\prime }]+\frac{{\varpi }^2}{\mathcal{l}}\left(a\hspace{0.17em}\mathit{sin}\hspace{0.17em}\psi \hspace{0.17em}\mathit{cos}\hspace{0.17em}\varpi \tau -b\hspace{0.17em}\mathit{cos}\hspace{0.17em}\psi \hspace{0.17em}\mathit{sin}\hspace{0.17em}\varpi \tau \right)=f_3\hspace{0.17em}\mathit{cos}\hspace{0.17em}{p}_3\tau ,$

(12)

$x^{\prime }+\alpha x-z^{\prime }=0.$

(13)

The prime represent differentiation concerning a variable $\tau$ . Equations (10)–(12) are considered nonlinear second-order ordinary differential equations (ODE).

It must be noted that in energy harvesting systems, continuous physical force is crucial for generating electricity efficiently. Here’s an explanation of how this force works and its role in energy harvesters:

Energy harvesters convert mechanical energy from physical forces into electrical energy. This process relies on a continuous external force to maintain the motion required for energy conversion. Many harvesters utilize vibrational forces, such as the motion of a vibrating rigid body. For instance, an electromagnetic device might be attached to a vibrating structure, as seen in Figure (1). To continuously generate electricity, the harvester must maintain consistent motion. For vibrational energy harvesters, this means that the vibrating source must produce ongoing oscillations according to the rotational harmonic moments at the points $O$ and $\mathrm{A}$ , and the harmonic physical force $F\left(t\right)$ in the spring’s direction. The applied force and moments cause the harvester (e.g., a mass-spring system in the electromagnetic harvester) to oscillate or vibrate.

Perturbation methodology

The primary goal of this section is to acquire the sought-after asymptotic solutions to the EOMs outlined in the provided system (10)-(13), specifically targeting up to the third-order of approximation. This will be achieved by leveraging the MST as a crucial perturbation technique. Furthermore, the systems’ resonance cases ⁴⁵ will be categorized. To achieve this, we will employ Taylor series expansion to approximate the trigonometric functions in $\phi$ and $\psi$ , as well as their derivatives up to the third-order, which is applicable in the vicinity of the static equilibrium position. Consequently, equations describing the system can be reformulated accordingly.

$z^{″}+c_1{z}^{\prime }+{\omega }_1^{2}z+\frac{{\phi }^2}{2}-\left(1+z\right){\phi }^{\prime 2}-\frac{s}{36}[\phi \psi \left({\phi }^2-6\right)\left({\psi }^2-6\right)+9\left({\phi }^2-2\right)({\psi }^2-2)]{\psi }^{\prime 2}$

$+\frac{s}{12}[\phi \left({\phi }^2-6\right)\left({\psi }^2-2\right)-\psi \left({\phi }^2-2\right)({\psi }^2-6)]{\psi }^{″}+\frac{{\varpi }^2}{\mathcal{l}}[b\left(\frac{{\phi }^3}{6}-\phi \right)\mathit{sin}\hspace{0.17em}\varpi \tau$

$+a\left(\frac{{\phi }^2}{2}-1\right)\mathit{cos}\hspace{0.17em}\varpi \tau ]+\gamma x=f_1\hspace{0.17em}\mathit{cos}\hspace{0.17em}{p}_1\tau ,$

(14)

${\left(z+1\right)}^2{\phi }^{″}+c_2{\phi }^{\prime }+\left(1+z\right)\left(\phi -\frac{{\phi }^3}{6}\right)+2\left(1+z\right)z^{\prime }{\phi }^{\prime }+\frac{s}{36}\left(1+z\right)[\phi \psi \left({\phi }^2-6\right)\left({\psi }^2-6\right)$

$+9\left({\phi }^2-2\right)({\psi }^2-2)]{\psi }^{″}+\frac{s}{12}\left(1+z\right)[\phi \left({\phi }^2-6\right)\left({\psi }^2-2\right)-\psi \left({\phi }^2-2\right)({\psi }^2-6)]{\psi }^{\prime 2}$

$+\frac{{\varpi }^2}{\mathcal{l}}\left(1+z\right)\left[b\left(\frac{{\phi }^2}{2}-1\right)\mathit{sin}\hspace{0.17em}\varpi \tau -a\left(\frac{{\phi }^3}{6}-\phi \right)\mathit{cos}\hspace{0.17em}\varpi \tau \right]=f_2\hspace{0.17em}\mathit{cos}\hspace{0.17em}{p}_2\tau ,$

(15)

$\left(s+{\omega }_2^2\right){\psi }^{″}+c_3{\psi }^{\prime }+\left(\psi -\frac{{\psi }^3}{6}\right)-\frac{1}{12}\left(1+z\right)[\phi \left({\phi }^2-6\right)\left({\psi }^2-2\right)-\psi \left({\phi }^2-2\right)({\psi }^2-6)]{\phi }^{\prime 2}$

$+\frac{1}{36}\left(1+z\right)[\phi \psi \left({\phi }^2-6\right)\left({\psi }^2-6\right)+9\left({\phi }^2-2\right)({\psi }^2-2)]{\phi }^{″}+\frac{1}{18}[\phi \psi \left({\phi }^2-6\right)\left({\psi }^2-6\right)$

$+9\left({\phi }^2-2\right)({\psi }^2-2)]z^{\prime }{\phi }^{\prime }+\frac{1}{12}[\phi \left({\phi }^2-6\right)\left({\psi }^2-2\right)-\psi \left({\phi }^2-2\right)({\psi }^2-6)]z^{″}+\frac{{\varpi }^2}{\mathcal{l}}\times$

$\left[b\left(\frac{{\psi }^2}{2}-1\right)\mathit{sin}\hspace{0.17em}\varpi \tau -a\left(\frac{{\psi }^3}{6}-\psi \right)\mathit{cos}\hspace{0.17em}\varpi \tau \right]=f_3\hspace{0.17em}\mathit{cos}\hspace{0.17em}{p}_3\tau ,$

(16)

$x^{\prime }+\alpha x-z^{\prime }=0.$

(17)

It is assumed that the generalized coordinates might be represented in relation to the order $0<\epsilon \ll 1$ . Hence, new variables $y,\varphi ,\zeta ,$ and $\chi$ are taken into consideration in the following manner:

$z\left(\tau \right)=\epsilon y\left(\tau ;\epsilon \right),\phi \left(\tau \right)=\epsilon \varphi \left(\tau ;\epsilon \right),\psi \left(\tau \right)=\epsilon \zeta \left(\tau ;\epsilon \right),x\left(\tau \right)=\epsilon \chi \left(\tau ;\epsilon \right).$

(18)

According to MST, we can present the function $y,\varphi ,\zeta ,$ and $\chi$ in terms of powers of $\epsilon$ in the following manner ⁴⁶ :

$y={\displaystyle \sum }_{k=0}^2{\epsilon }^k{y}_{k+1}\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^3\right),$

$\varphi ={\displaystyle \sum }_{k=0}^2{\epsilon }^k{\varphi }_{k+1}\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^3\right),$

$\zeta ={\displaystyle \sum }_{k=0}^2{\epsilon }^k{\zeta }_{k+1}\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^3\right)$

$\chi ={\displaystyle \sum }_{k=0}^2{\epsilon }^k{\chi }_{k+1}\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^3\right),$

(19)

where ${\tau }_n={\epsilon }^n\tau \left(n=0,1,2\right)$ perform different time scales, in which ${\tau }_0$ and ${\tau }_1,{\tau }_2$ act as the fast and slow time scales, respectively. In order to establish functions $y,\varphi ,\zeta ,$ and $\chi$ based on ${\tau }_n$ , we convert $\tau$ to ${\tau }_n$ using specific operators ⁴⁷

$\frac{d}{d\tau }=\frac{\partial }{\partial {\tau }_0}+\epsilon \frac{\partial }{\partial {\tau }_1}+{\epsilon }^2\frac{\partial }{\partial {\tau }_2},$

$\frac{d^2}{d{\tau }^2}=\frac{{\partial }^2}{\partial {\tau }_0^2}+2\epsilon \frac{{\partial }^2}{\partial {\tau }_0\partial {\tau }_1}+{\epsilon }^2\left(\frac{{\partial }^2}{\partial {\tau }_1^2}+2\frac{{\partial }^2}{\partial {\tau }_0\partial {\tau }_2}\right)+O\left({\epsilon }^3\right).$

(20)

It’s crucial to acknowledge that terms of $O\left({\epsilon }^3\right)$ and beyond are deemed insignificant due to their comparatively small magnitudes.

The eccentricity, the damping coefficients, and generalized forces can be represented in terms of $\epsilon$ as shown below:

$a={\epsilon }^2\stackrel{\sim }{a},b={\epsilon }^2\stackrel{\sim }{b},s={\epsilon }^2\stackrel{\sim }{s},\gamma ={\epsilon }^2\stackrel{\sim }{\gamma },$

$\alpha =\stackrel{\sim }{\alpha },c_j={\epsilon }^2{\stackrel{\sim }{c}}_j,f_j={\epsilon }^3{\stackrel{\sim }{f}}_j,\left(j=1,2,3\right)$

(21)

where the parameters $\stackrel{\sim }{a},\stackrel{\sim }{b},\stackrel{\sim }{s},\stackrel{\sim }{\gamma },\stackrel{\sim }{\alpha },{\stackrel{\sim }{c}}_j,$ and ${\stackrel{\sim }{f}}_j$ are of order 1.

Inserting equations (18)–(21) into equations (14)–(17) and subsequently comparing the coefficients of the same powers of $\epsilon ,{\epsilon }^2$ and ${\epsilon }^3$ , we can derive the below three distinct sets of partial differential equations (PDEs):Order of $\epsilon$

$\frac{{\partial }^2{y}_1}{\partial {\tau }_0^2}+{\omega }_1^2{y}_1=0,$

(22)

$\frac{{\partial }^2{\varphi }_1}{\partial {\tau }_0^2}+{\varphi }_1=0,$

(23)

$\frac{{\partial }^2{\zeta }_1}{\partial {\tau }_0^2}+\frac{1}{{\omega }_2^2}{\zeta }_1=-\frac{1}{{\omega }_2^2}\frac{{\partial }^2{\varphi }_1}{\partial {\tau }_0^2},$

(24)

$\frac{\partial {\chi }_1}{\partial {\tau }_0}+\stackrel{\sim }{\alpha }{\chi }_1=\frac{\partial y_1}{\partial {\tau }_0}.$

(25)

Order of ${\epsilon }^2$

$\frac{{\partial }^2{y}_2}{\partial {\tau }_0^2}+{\omega }_1^2{y}_2=\frac{{\varpi }^2}{\mathcal{l}}\stackrel{\sim }{a}\mathit{cos}\left(\varpi {\tau }_0\right)-2\frac{{\partial }^2{y}_1}{\partial {\tau }_0{\tau }_1}+{\left(\frac{\partial {\varphi }_1}{\partial {\tau }_0}\right)}^2-\frac{1}{2}{{\varphi }_1}^2,$

(26)

$\frac{{\partial }^2{\varphi }_2}{\partial {\tau }_0^2}+{\varphi }_2=\frac{{\varpi }^2}{\mathcal{l}}\stackrel{\sim }{b}\mathit{sin}\left(\varpi {\tau }_0\right)-2\left[\frac{{\partial }^2{\varphi }_1}{\partial {\tau }_0\partial {\tau }_1}+y_1\frac{{\partial }^2{\varphi }_1}{\partial {\tau }_0^2}+\frac{\partial y_1}{\partial {\tau }_0}\frac{\partial {\varphi }_1}{\partial {\tau }_0}\right]-y_1{\varphi }_1,$

(27)

$\frac{{\partial }^2{\zeta }_2}{\partial {\tau }_0^2}+\frac{1}{{\omega }_2^2}{\zeta }_2=\frac{{\varpi }^2}{\mathcal{l}{\omega }_2^2}\stackrel{\sim }{b}\mathit{sin}\left(\varpi {\tau }_0\right)-2\frac{{\partial }^2{\zeta }_1}{\partial {\tau }_0\partial {\tau }_1}-\frac{2}{{\omega }_2^2}\left(\frac{{\partial }^2{\varphi }_1}{\partial {\tau }_0\partial {\tau }_1}+\frac{\partial y_1}{\partial {\tau }_0}\frac{\partial {\varphi }_1}{\partial {\tau }_0}\right)$

$+\frac{1}{{\omega }_2^2}\left({\zeta }_1\frac{{\partial }^2{y}_1}{\partial {\tau }_0^2}-{\varphi }_1\frac{{\partial }^2{y}_1}{\partial {\tau }_0^2}-y_1\frac{{\partial }^2{\varphi }_1}{\partial {\tau }_0^2}\right)-\frac{1}{{\omega }_2^2}\frac{{\partial }^2{\varphi }_2}{\partial {\tau }_0^2},$

(28)

$\frac{\partial {\chi }_2}{\partial {\tau }_0}+\stackrel{\sim }{\alpha }{\chi }_2=\frac{\partial y_2}{\partial {\tau }_0}-\frac{\partial {\chi }_1}{\partial {\tau }_1}+\frac{\partial y_1}{\partial {\tau }_1}.$

(29)

Order of ${\epsilon }^3$

$\frac{{\partial }^2{y}_3}{\partial {\tau }_0^2}+{\omega }_1^2{y}_3={\stackrel{\sim }{f}}_1\hspace{0.17em}\mathit{cos}\left(p_1{\tau }_0\right)+\frac{{\varpi }^2}{\mathcal{l}}\stackrel{\sim }{b}{\varphi }_1\hspace{0.17em}\mathit{sin}\left(\varpi {\tau }_0\right)-\frac{{\partial }^2{y}_1}{\partial {\tau }_1^2}$

$-2\left(\frac{{\partial }^2{y}_1}{\partial {\tau }_0\partial {\tau }_2}+\frac{{\partial }^2{y}_2}{\partial {\tau }_0\partial {\tau }_1}\right)+2\left(\frac{\partial {\varphi }_1}{\partial {\tau }_0}\frac{\partial {\varphi }_2}{\partial {\tau }_0}+\frac{\partial {\varphi }_1}{\partial {\tau }_0}\frac{\partial {\varphi }_1}{\partial {\tau }_1}\right)$

$+y_1{\left(\frac{\partial {\varphi }_1}{\partial {\tau }_0}\right)}^2-{\stackrel{\sim }{c}}_1\frac{\partial y_1}{\partial {\tau }_0}-{\varphi }_1{\varphi }_2-\stackrel{\sim }{\gamma }{\chi }_1,$

(30)

$\frac{{\partial }^2{\varphi }_3}{\partial {\tau }_0^2}+{\varphi }_3={\stackrel{\sim }{f}}_2\hspace{0.17em}\mathit{cos}\left(p_2{\tau }_0\right)-\frac{{\varpi }^2}{\mathcal{l}}\stackrel{\sim }{a}{\varphi }_1\hspace{0.17em}\mathit{cos}\left(\varpi {\tau }_0\right)+\frac{{\varpi }^2}{\mathcal{l}}\stackrel{\sim }{b}{y}_1\hspace{0.17em}\mathit{sin}\left(\varpi {\tau }_0\right)$

$-2\left(\frac{{\partial }^2{\varphi }_1}{\partial {\tau }_0\partial {\tau }_2}+\frac{{\partial }^2{\varphi }_2}{\partial {\tau }_0\partial {\tau }_1}\right)-2\left(y_1\frac{{\partial }^2{\varphi }_2}{\partial {\tau }_0^2}+y_2\frac{{\partial }^2{\varphi }_1}{\partial {\tau }_0^2}\right)-4y_1\frac{{\partial }^2{\varphi }_1}{\partial {\tau }_0\partial {\tau }_1}$

$-2y_1\frac{\partial {\varphi }_1}{\partial {\tau }_0}\frac{\partial y_1}{\partial {\tau }_0}-2\left(\frac{\partial {\varphi }_2}{\partial {\tau }_0}\frac{\partial y_1}{\partial {\tau }_0}+\frac{\partial {\varphi }_1}{\partial {\tau }_0}\frac{\partial y_1}{\partial {\tau }_1}+\frac{\partial {\varphi }_1}{\partial {\tau }_0}\frac{\partial y_2}{\partial {\tau }_0}+\frac{\partial {\varphi }_1}{\partial {\tau }_1}\frac{\partial y_1}{\partial {\tau }_0}\right)$

$-y_1{\varphi }_2-y_2{\varphi }_1-y_1^2\frac{{\partial }^2{\varphi }_1}{\partial {{\tau }_0}^2}-{\stackrel{\sim }{c}}_2\frac{\partial {\varphi }_1}{\partial {\tau }_0}+\frac{1}{6}{\varphi }_1^3-\frac{{\partial }^2{\varphi }_1}{\partial {\tau }_1^2}-\stackrel{\sim }{s}\frac{{\partial }^2{\zeta }_1}{\partial {\tau }_0^2},$

(31)

$\frac{{\partial }^2{\zeta }_3}{\partial {\tau }_0^2}+\frac{1}{{\omega }_2^2}{\zeta }_3=\frac{{\stackrel{\sim }{f}}_3}{{\omega }_2^2}\mathit{cos}\left(p_3{\tau }_0\right)-\frac{{\varpi }^2}{\mathcal{l}{\omega }_2^2}\stackrel{\sim }{a}{\zeta }_1\hspace{0.17em}\mathit{cos}\left(\varpi {\tau }_0\right)-2\left(\frac{{\partial }^2{\zeta }_1}{\partial {\tau }_0\partial {\tau }_2}+\frac{{\partial }^2{\zeta }_2}{\partial {\tau }_0\partial {\tau }_1}\right)$

$-\frac{{\partial }^2{\zeta }_1}{\partial {\tau }_1^2}-\frac{1}{{\omega }_2^2}\frac{{\partial }^2{\varphi }_1}{\partial {\tau }_1^2}-\frac{1}{{\omega }_2^2}\frac{{\partial }^2{\varphi }_3}{\partial {\tau }_0^2}-\frac{{\stackrel{\sim }{c}}_3}{{\omega }_2^2}\frac{\partial {\zeta }_1}{\partial {\tau }_0}-\frac{\stackrel{\sim }{s}}{{\omega }_2^2}\frac{{\partial }^2{\zeta }_1}{\partial {\tau }_0^2}+\frac{1}{6{\omega }_2^2}{\zeta }_1^3$

$-\frac{1}{{\omega }_2^2}\left({\varphi }_1\frac{{\partial }^2{y}_2}{\partial {\tau }_0^2}+{\varphi }_2\frac{{\partial }^2{y}_1}{\partial {\tau }_0^2}-{\zeta }_1\frac{{\partial }^2{y}_2}{\partial {\tau }_0^2}-{\zeta }_2\frac{{\partial }^2{y}_1}{\partial {{\tau }_0}^2}+y_1\frac{{\partial }^2{\varphi }_2}{\partial {{\tau }_0}^2}+y_2\frac{{\partial }^2{\varphi }_1}{\partial {{\tau }_0}^2}\right)$

$-\frac{2}{{\omega }_2^2}\left(\frac{{\partial }^2{\varphi }_1}{\partial {\tau }_0\partial {\tau }_2}+\frac{{\partial }^2{\varphi }_2}{\partial {\tau }_0\partial {\tau }_1}\right)-\frac{2}{{\omega }_2^2}\left({\varphi }_1\frac{{\partial }^2{y}_1}{\partial {\tau }_0\partial {\tau }_1}+y_1\frac{{\partial }^2{\varphi }_1}{\partial {\tau }_0\partial {\tau }_1}-{\zeta }_1\frac{{\partial }^2{y}_1}{\partial {\tau }_0\partial {\tau }_1}\right)$

$-\frac{2}{{\omega }_2^2}\left(\frac{\partial {\varphi }_2}{\partial {\tau }_0}\frac{\partial y_1}{\partial {\tau }_0}+\frac{\partial {\varphi }_1}{\partial {\tau }_0}\frac{\partial y_2}{\partial {\tau }_0}+\frac{\partial {\varphi }_1}{\partial {\tau }_0}\frac{\partial y_1}{\partial {\tau }_1}+\frac{\partial {\varphi }_1}{\partial {\tau }_1}\frac{\partial y_1}{\partial {\tau }_0}\right)+\frac{1}{{\omega }_2^2}{\varphi }_1{\left(\frac{\partial {\varphi }_1}{\partial {\tau }_0}\right)}^2$

$-\frac{1}{{\omega }_2^2}{\zeta }_1{\left(\frac{\partial {\varphi }_1}{\partial {\tau }_0}\right)}^2-\frac{1}{{\omega }_2^2}{\zeta }_1{\varphi }_1\frac{{\partial }^2{\varphi }_1}{\partial {\tau }_0^2}+\frac{1}{2{\omega }_2^2}\left({\zeta }_1^2\frac{{\partial }^2{\varphi }_1}{\partial {{\tau }_0}^2}+{\varphi }_1^2\frac{{\partial }^2{\varphi }_1}{\partial {{\tau }_0}^2}\right),$

(32)

$\frac{\partial {\chi }_3}{\partial {\tau }_0}+\stackrel{\sim }{\alpha }{\chi }_3=\frac{\partial y_3}{\partial {\tau }_0}-\left(\frac{\partial {\chi }_2}{\partial {\tau }_1}+\frac{\partial {\chi }_1}{\partial {\tau }_2}\right)+\frac{\partial y_2}{\partial {\tau }_1}+\frac{\partial y_1}{\partial {\tau }_2}.$

(33)

It’s worth noting that the three systems of order of $\epsilon$ create a total of PDEs that can be solved one after the other. In order to accomplish this task, let’s begin with equations (22)–(25), which possess the following general solutions:

$y_1=A_1{e}^{i{\omega }_1{\tau }_0}+{\overline{A}}_1{e}^{-i{\omega }_1{\tau }_0},$

(34)

${\varphi }_1=A_2{e}^{i{\tau }_0}+{\overline{A}}_2{e}^{-i{\tau }_0},$

(35)

${\zeta }_1=A_3{e}^{\frac{i{\tau }_0}{{\omega }_2}}+{\overline{A}}_3{e}^{\frac{-i{\tau }_0}{{\omega }_2}}+\frac{1}{\left(1-{\omega }_2^2\right)}\left(A_2{e}^{i{\tau }_0}+{\overline{A}}_2{e}^{-i{\tau }_0}\right),$

(36)

where $A_j=A_j\left({\tau }_1,{\tau }_2\right)\left(j=1,2,3\right)$ are unidentified complex functions, while $\overline{A_j}$ indicate their complex conjugates.

In order to find the overall solution of ${\chi }_1$ , substituting equation (34) into equation (25) to derive a nonhomogeneous PDE that can be solved using a specific form

${\chi }_1=B_1{e}^{-\alpha {\tau }_0}+i{\omega }_1\left[\frac{A_1}{\left(\alpha +i{\omega }_1\right)}{e}^{i{\omega }_1{\tau }_0}-\frac{{\overline{A}}_1}{\left(\alpha -i{\omega }_1\right)}{e}^{-i{\omega }_1{\tau }_0}\right];i=\sqrt{-1},$

(37)

where, the definition of $B_1$ mirror the previous definition of $A_j$ .

Inserting equations (34)-(37) into (26)-(29) and eliminating terms resulting in secular behavior, we can derive a consistent asymptotic second-order of solutions as follows:

$y_2=-\frac{3A_2^2}{2\left({\omega }_1^2-4\right)}{e}^{2i{\tau }_0}+\frac{{\varpi }^2\stackrel{\sim }{a}}{2\mathcal{l}\left({\omega }_1^2-{\varpi }^2\right)}{e}^{i\varpi {\tau }_0}+\frac{1}{{\omega }_1^2}{A}_2{\overline{A}}_2+CC,$

(38)

${\varphi }_2=-\frac{\left(1+2{\omega }_1\right)}{{\omega }_1\left(2+{\omega }_1\right)}{A}_1{A}_2{e}^{i\left(1+{\omega }_1\right){\tau }_0}+\frac{\left(1-2{\omega }_1\right)}{{\omega }_1\left(2-{\omega }_1\right)}{\overline{A}}_1{A}_2{e}^{i\left(1-{\omega }_1\right){\tau }_0}-\frac{i{\varpi }^2\stackrel{\sim }{b}}{2\mathcal{l}\left(1-{\varpi }^2\right)}{e}^{i\varpi {\tau }_0}+CC,$

(39)

${\zeta }_2=\frac{1+2{\omega }_1+\left({\omega }_1-1\right){\left({\omega }_1+1\right)}^3{\omega }_2^2}{{\omega }_1\left({\omega }_1+2\right)\left({\omega }_2^2-1\right)\left[1-{\omega }_2^2{\left(1+{\omega }_1\right)}^2\right]}{A}_1{A}_2{e}^{i\left(1+{\omega }_1\right){\tau }_0}+{\overline{A}}_1{A}_2{e}^{i\left(1-{\omega }_1\right){\tau }_0}$

$\times \frac{1-2{\omega }_1+\left({\omega }_1+1\right){\left({\omega }_1-1\right)}^3{\omega }_2^2}{{\omega }_1\left({\omega }_1-2\right)\left({\omega }_2^2-1\right)\left[1-{\omega }_2^2{\left(1-{\omega }_1\right)}^2\right]}-\frac{{\omega }_1^2}{\left[1-{\left(1+{\omega }_1{\omega }_2\right)}^2\right]}{A}_1{A}_3{e}^{i\left(\frac{1+{\omega }_1{\omega }_2}{{\omega }_2}\right){\tau }_0}$

$-\frac{{\omega }_1^2}{\left[1-{\left(1-{\omega }_1{\omega }_2\right)}^2\right]}{\overline{A}}_1{A}_3{e}^{i\left(\frac{1-{\omega }_1{\omega }_2}{{\omega }_2}\right){\tau }_0}-\frac{i{\varpi }^2\stackrel{\sim }{b}}{2\mathcal{l}\left(1-{\varpi }^2\right)\left(1-{\varpi }^2{\omega }_2^2\right)}{e}^{i\varpi {\tau }_0}+CC,$

(40)

${\chi }_2=-\frac{3i}{\left(\alpha +2i\right)\left({\omega }_1^2-4\right)}{A}_2^2{e}^{2i{\tau }_0}+\frac{i{\varpi }^3\stackrel{\sim }{a}}{2\mathcal{l}\left({\omega }_1^2-{\varpi }^2\right)\left(\alpha +i\varpi \right)}{e}^{i\varpi {\tau }_0}+CC,$

(41)

where $CC$ represents the preceding terms’ conjugate.

Based on the provided information, it is possible to establish the criteria of second-order approximations for eliminating secular terms as in the following specified forms:

$\frac{\partial A_1}{\partial {\tau }_1}=0,\frac{\partial A_2}{\partial {\tau }_1}=0,\frac{\partial A_3}{\partial {\tau }_1}=0,\frac{\partial B_1}{\partial {\tau }_1}=0.$

(42)

Hence, $A_j=A_j\left({\tau }_2\right)$ and $B_1=B_1\left({\tau }_2\right)$ .

Making use of equations (34)–(37) and (38)-(42) into (30)-(33), the third-order approximate solutions can be expressed through the elimination of terms that lead to secular ones. To eliminate these terms, certain conditions need to be met

$2i{\omega }_1\frac{\partial A_1}{\partial {\tau }_2}+A_1\left[i{\stackrel{\sim }{c}}_1{\omega }_1+\frac{i\stackrel{\sim }{\gamma }{\omega }_1}{\left(\alpha +i{\omega }_1\right)}+\frac{6\left({\omega }_1^2-1\right)}{\left({\omega }_1^2-4\right)}{A}_2{\overline{A}}_2\right]=0,$

(43)

$2i\frac{\partial A_2}{\partial {\tau }_2}+A_2\left[i{\stackrel{\sim }{c}}_2-\frac{\stackrel{\sim }{s}}{\left(1-{\omega }_2^2\right)}+\frac{6\left({\omega }_1^2-1\right)}{\left({\omega }_1^2-4\right)}{A}_1{\overline{A}}_1-\frac{{\omega }_1^4+7{\omega }_1^2-8}{2{\omega }_1^2\left({\omega }_1^2-4\right)}{A}_2{\overline{A}}_2\right]=0,$

(44)

$\frac{2i}{{\omega }_2}\frac{\partial A_3}{\partial {\tau }_2}+A_3[\frac{i{\stackrel{\sim }{c}}_3}{{\omega }_2^3}-\frac{\stackrel{\sim }{s}}{{\omega }_2^2\left({\omega }_2^2-1\right)}-\frac{1}{2{\omega }_2^2}{A}_3{\overline{A}}_3+\frac{\left(1-2{\omega }_2^2\right)}{{\omega }_2^2{\left(1-{\omega }_2^2\right)}^2}{A}_2{\overline{A}}_2$

$-\frac{2{\omega }_1^4}{{\omega }_2^2\left(4-{\omega }_1^2{\omega }_2^2\right)}{A}_1{\overline{A}}_1]=0,$

(45)

$\frac{\partial B_1}{\partial {\tau }_2}-\frac{\alpha \stackrel{\sim }{\gamma }{B}_1}{\left({\alpha }^2+{\omega }_1^2\right)}=0.$

(46)

Subsequently, the solutions of the third approximations $y_3,{\varphi }_3,{\zeta }_3$ , and ${\chi }_3$ assume the following forms:

$y_3=\frac{{\stackrel{\sim }{f}}_1{e}^{i{p}_1{\tau }_0}}{2\left({\omega }_1^2-p_1^2\right)}+\frac{i{\varpi }^3\left(\varpi +2\right)\stackrel{\sim }{b}{A}_2{e}^{i\left(\varpi +1\right){\tau }_0}}{2\mathcal{l}\left[{\omega }_1^2-{\left(\varpi +1\right)}^2\right]\left(1-{\varpi }^2\right)}+\frac{i{\varpi }^3\left(\varpi -2\right)\stackrel{\sim }{b}{\overline{A}}_2{e}^{i\left(\varpi -1\right){\tau }_0}}{2\mathcal{l}\left[{\omega }_1^2-{\left(\varpi -1\right)}^2\right]\left(1-{\varpi }^2\right)}$

$-\frac{3\left({\omega }_1+1\right)A_1{A}_2^2{e}^{i\left(2+{\omega }_1\right){\tau }_0}}{4{\omega }_1\left({\omega }_1+2\right)}+\frac{3\left({\omega }_1-1\right){\overline{A}}_1{A}_2^2{e}^{i\left(2-{\omega }_1\right){\tau }_0}}{4{\omega }_1\left({\omega }_1-2\right)}-\frac{\stackrel{\sim }{\gamma }{B}_1{e}^{-\alpha {\tau }_0}}{\left({\alpha }^2+{\omega }_1^2\right)}+CC,$

(47)

${\varphi }_3=\frac{{\stackrel{\sim }{f}}_2{e}^{i{p}_2{\tau }_0}}{2\left(1-p_2^2\right)}-\frac{\left({\omega }_1^2-49\right)A_2^3{e}^{3i{\tau }_0}}{48\left({\omega }_1^2-4\right)}+\frac{\stackrel{\sim }{s}{A}_3{e}^{\left(\frac{i}{{\omega }_2}\right){\tau }_0}}{\left({\omega }_2^2-1\right)}-\frac{i{\varpi }^3\stackrel{\sim }{b}\left(2{\omega }_1+\varpi \right)A_1{e}^{i\left(\varpi +{\omega }_1\right){\tau }_0}}{2\mathcal{l}\left(1-{\varpi }^2\right)\left(1-{\left(\varpi +{\omega }_1\right)}^2\right)}$

$+\frac{i{\varpi }^3\stackrel{\sim }{b}\left(2{\omega }_1-\varpi \right){\overline{A}}_1{e}^{i\left(\varpi -{\omega }_1\right){\tau }_0}}{2\mathcal{l}\left(1-{\varpi }^2\right)\left(1-{\left(\varpi -{\omega }_1\right)}^2\right)}-\frac{\varpi \stackrel{\sim }{a}\left(1+2\varpi +{\varpi }^2-{\omega }_1^2\right)A_2{e}^{i\left(\varpi +1\right){\tau }_0}}{2\mathcal{l}\left({\omega }_1^2-{\varpi }^2\right)\left(2+\varpi \right)}$

$+\frac{\varpi \stackrel{\sim }{a}\left(1-2\varpi +{\varpi }^2-{\omega }_1^2\right){\overline{A}}_2{e}^{i\left(\varpi -1\right){\tau }_0}}{2\mathcal{l}\left({\omega }_1^2-{\varpi }^2\right)\left(2-\varpi \right)}+\frac{\left(6{\omega }_1^3+11{\omega }_1^2+6{\omega }_1+1\right)A_1^2{A}_2{e}^{i\left(1+2{\omega }_1\right){\tau }_0}}{4{\omega }_1^2\left({\omega }_1+2\right)\left(1+{\omega }_1\right)}$

$+\frac{\left(6{\omega }_1^3-11{\omega }_1^2+6{\omega }_1-1\right){\overline{A}}_1^2{A}_2{e}^{i\left(1-2{\omega }_1\right){\tau }_0}}{4{\omega }_1^2\left({\omega }_1-2\right)\left(1-{\omega }_1\right)}+CC,$

(48)

${\zeta }_3=\frac{{\stackrel{\sim }{f}}_3{e}^{i{p}_3{\tau }_0}}{2\left(1-{\omega }_2^2{p}_3^2\right)}+\frac{p_2^2{\stackrel{\sim }{f}}_2{e}^{i{p}_2{\tau }_0}}{2\left(1-p_2^2\right)\left(1-{\omega }_2^2{p}_2^2\right)}-\frac{A_3^3}{48}{e}^{\left(\frac{3i}{{\omega }_2}\right){\tau }_0}+E_1{A}_2^3{e}^{3i{\tau }_0}+E_2{A}_2{e}^{i{\tau }_0}$

$+E_3{A}_2^2{A}_3{e}^{i\left(\frac{2{\omega }_2+1}{{\omega }_2}\right){\tau }_0}+E_4{A}_2^2{\overline{A}}_3{e}^{i\left(\frac{2{\omega }_2-1}{{\omega }_2}\right){\tau }_0}+E_5{A}_2{A}_3^2{e}^{i\left(\frac{{\omega }_2+2}{{\omega }_2}\right){\tau }_0}+E_6{A}_2{\overline{A}}_3^2{e}^{i\left(\frac{{\omega }_2-2}{{\omega }_2}\right){\tau }_0}$

$+E_7{A}_2{e}^{i\left(\varpi +1\right){\tau }_0}+E_8{\overline{A}}_2{e}^{i\left(\varpi -1\right){\tau }_0}-E_9{A}_1^2{A}_2{e}^{i\left(1+2{\omega }_1\right){\tau }_0}-E_{10}{\overline{A}}_1^2{A}_2{e}^{i\left(1-2{\omega }_1\right){\tau }_0}$

$+E_{11}{A}_1{e}^{i\left(\varpi +{\omega }_1\right){\tau }_0}+E_{12}{\overline{A}}_1{e}^{i\left(\varpi -{\omega }_1\right){\tau }_0}+E_{13}{A}_3{e}^{i\left(\frac{\varpi {\omega }_2+1}{{\omega }_2}\right){\tau }_0}+E_{14}{\overline{A}}_3{e}^{i\left(\frac{\varpi {\omega }_2-1}{{\omega }_2}\right){\tau }_0}$

$+E_{15}{A}_1^2{A}_3{e}^{i\left(\frac{1+2{\omega }_1{\omega }_2}{{\omega }_2}\right){\tau }_0}+E_{16}{\overline{A}}_1^2{A}_3{e}^{i\left(\frac{1-2{\omega }_1{\omega }_2}{{\omega }_2}\right){\tau }_0}+CC,$

(49)

where $E_{\upsilon }\left(\upsilon =1,2,...,16\right)$ can be obtained easily (see Appendix I).

${\chi }_3=\frac{i{p}_1{\stackrel{\sim }{f}}_1{e}^{i{p}_1{\tau }_0}}{2\left({\omega }_1^2-p_1^2\right)\left(\alpha +i{p}_1\right)}+\frac{A_1{e}^{i{\omega }_1{\tau }_0}}{\left(\alpha +i{\omega }_1\right)}\left[1-\frac{i{\omega }_1}{\left(\alpha +i{\omega }_1\right)}\right][-\frac{{\stackrel{\sim }{c}}_1}{2}-\frac{\stackrel{\sim }{\gamma }}{2\left(\alpha +i{\omega }_1\right)}$

$+\frac{3i\left({\omega }_1^2-1\right)}{{\omega }_1\left({\omega }_1^2-4\right)}{A}_2{\overline{A}}_2]-\frac{{\varpi }^3\stackrel{\sim }{b}\left(\varpi +2\right)A_2{e}^{i\left(\varpi +1\right){\tau }_0}}{2\mathcal{l}\left(1-\varpi \right)\left[{\omega }_1^2-{\left(\varpi +1\right)}^2\right]\left[\alpha +i\left(\varpi +1\right)\right]}$

$+\frac{{\varpi }^3\stackrel{\sim }{b}\left(\varpi -2\right){\overline{A}}_2{e}^{i\left(\varpi -1\right){\tau }_0}}{2\mathcal{l}\left(\varpi +1\right)\left[{\omega }_1^2-{\left(\varpi -1\right)}^2\right]\left[\alpha +i\left(\varpi -1\right)\right]}+\frac{3i{\left({\omega }_1+1\right)}^2{A}_1{A}_2^2{e}^{i\left(2+{\omega }_1\right){\tau }_0}}{{\omega }_1\left[{\omega }_1^2-{\left(2+{\omega }_1\right)}^2\right]\left[\alpha +i\left(2+{\omega }_1\right)\right]}$

$-\frac{3i{\left({\omega }_1-1\right)}^2{\overline{A}}_1{A}_2^2{e}^{i\left(2-{\omega }_1\right){\tau }_0}}{{\omega }_1\left[{\omega }_1^2-{\left(2-{\omega }_1\right)}^2\right]\left[\alpha +i\left(2-{\omega }_1\right)\right]}+CC.$

(50)

Now, we can go back to determine the functions $A_j$ and $B_1$ by applying the elimination criteria for secular terms (42)-(46), in addition to the specified initial conditions.

$z\left(0\right)={\eta }_{01},\phi \left(0\right)={\eta }_{03},\psi \left(0\right)={\eta }_{05},x\left(0\right)={\eta }_{07},$

$z^{\prime }\left(0\right)={\eta }_{02},{\phi }^{\prime }\left(0\right)={\eta }_{04},{\psi }^{\prime }\left(0\right)={\eta }_{06},$

(51)

where ${\eta }_{01},{\eta }_{02},...,{\eta }_{07}$ are known quantities.

It must be noted that approximate solutions can provide insights into the behavior of complex systems without the need for overly complex mathematical analysis. Even if an exact solution is not feasible, approximate solutions can help us understand the system’s response to various inputs and initial conditions. These solutions can be used to validate analytical models developed for vibrating dynamical systems, predict the behavior of vibrating dynamical systems, and make timely adjustments to optimize performance or prevent failures. By providing insights into the behavior of vibrating dynamical systems under different conditions, it can help assess the risks associated with system failures or performance degradation. This is particularly important in safety-critical applications such as aerospace or automotive engineering.

Resonance scenarios and stability assessment

In this part, our attention is directed towards categorizing resonance scenarios. It is known that resonance can occur when the denominators of the obtained approximations, except the first one, reach zero. Hence, equations (38)–(41), and (47)-(50) unveil multiple resonance cases, neatly falling into three distinct partitions.

The first scenario is the primary external resonance that can be done at $p_1\approx {\omega }_1,p_2\approx 1,$ $p_3\approx \frac{1}{{\omega }_2}$ , the second case is known by internal resonance that has been found at ${\omega }_1\approx 0,{\omega }_1\approx 1,{\omega }_1\approx 2,$ ${\omega }_1\approx \varpi ,{\omega }_2\approx 0,$ ${\omega }_2\approx \pm 1,{\omega }_2\approx 3,$ ${\omega }_2\approx \frac{1}{3},\varpi =0,$ $\varpi =1,\varpi =2,\varpi =\frac{1}{{\omega }_2},\varpi =\frac{2}{{\omega }_2},$ whereas the third resonance scenario is a combination resonance which can arise at ${\omega }_1-\varpi \approx \pm 1,{\omega }_1+\varpi \approx 1,$

${\omega }_1\approx \pm \left[\left(\frac{1}{{\omega }_2}\right)-1\right],\varpi \approx \left(\frac{1}{{\omega }_2}\right)\pm 1,{\omega }_1\pm \varpi \approx \frac{1}{{\omega }_2}.$

It is worth noting that when any of these conditions are met, the system’s behavior becomes highly intricate. Additionally, if the vibrations' values deviate significantly from instances of resonance, the aforementioned estimated solutions will remain valid.^48,49

In order to assess the stability of the system, it is essential to determine the conditions under which solutions can be found and to derive the equations that account for any necessary modifications. Therefore, we will investigate three cases involving primary external resonance occurring concurrently, specifically when considering a combination of $p_1\approx {\omega }_1,$ $p_2\approx 1,$ and $p_3\approx \frac{1}{{\omega }_2}$ which describes the nearness of $p_1,p_2,$ and $p_3$ to ${\omega }_1,$ $1,$ and $\frac{1}{{\omega }_2}$ , respectively. Next, the parameters of detuning ${\sigma }_j\left(j=1,2,3\right)$ can be used as follows:

$p_1={\omega }_1+{\sigma }_1,p_2=1+{\sigma }_2,p_3=\left(\frac{1}{{\omega }_2}\right)+{\sigma }_3.$

(52)

These parameters serve as indicators of the degree to which the oscillations deviate from the tough resonance point. As a result, we represent them relation to ${\epsilon }^2$ according to the following notations

${\sigma }_j={\epsilon }^2{\stackrel{\sim }{\sigma }}_j,$

(53)

Inserting equations (52) and (53) into equations (26)–(29) and (30)-(33), and removing terms that yield secular ones, the solvability requirements can be derived as follows:

- For the approximation of second-order:

$\frac{\partial A_1}{\partial {\tau }_1}=0,\frac{\partial A_2}{\partial {\tau }_1}=0,\frac{\partial A_3}{\partial {\tau }_1}=0,\frac{\partial B_1}{\partial {\tau }_1}=0.$

(54)

- For the approximation of third-order:

${\stackrel{\sim }{f}}_1{e}^{i{\stackrel{\sim }{\sigma }}_1{\tau }_2}-4i{\omega }_1\frac{\partial A_1}{\partial {\tau }_2}-2A_1\left[i{\stackrel{\sim }{c}}_1{\omega }_1+\frac{i\stackrel{\sim }{\gamma }{\omega }_1}{\left(\alpha +i{\omega }_1\right)}+\frac{6\left({\omega }_1^2-1\right)}{\left({\omega }_1^2-4\right)}{A}_2{\overline{A}}_2\right]=0.$

(55)

${\stackrel{\sim }{f}}_2{e}^{i{\stackrel{\sim }{\sigma }}_2{\tau }_2}-4i\frac{\partial A_2}{\partial {\tau }_2}-2A_2[i{\stackrel{\sim }{c}}_2-\frac{\stackrel{\sim }{s}}{\left(1-{\omega }_2^2\right)}+\frac{6\left({\omega }_1^2-1\right)}{\left({\omega }_1^2-4\right)}{A}_1{\overline{A}}_1$

$-\frac{{\omega }_1^4+7{\omega }_1^2-8}{2{\omega }_1^2\left({\omega }_1^2-4\right)}{A}_2{\overline{A}}_2]=0.$

${\stackrel{\sim }{f}}_3{e}^{i{\stackrel{\sim }{\sigma }}_3{\tau }_2}-4i{\omega }_2\frac{\partial A_3}{\partial {\tau }_2}-A_3[\frac{2i{\stackrel{\sim }{c}}_3}{{\omega }_2}-\frac{2\stackrel{\sim }{s}}{\left({\omega }_2^2-1\right)}-A_3{\overline{A}}_3+\frac{2\left(1-2{\omega }_2^2\right)}{{\left(1-{\omega }_2^2\right)}^2}{A}_2{\overline{A}}_2$

$-\frac{4{\omega }_1^4}{\left(4-{\omega }_1^2{\omega }_2^2\right)}{A}_1{\overline{A}}_1]=0.$

$\frac{\partial B_1}{\partial {\tau }_2}-\frac{\alpha \stackrel{\sim }{\gamma }{B}_1}{\left({\alpha }^2+{\omega }_1^2\right)}=0.$

Using the last equation of (55) along with the conditions (54), the expression of $B_1$ can be obtained in the form

$B_1=e^{\left(\frac{\alpha \stackrel{\sim }{\gamma }}{{\alpha }^2+{\omega }_1^2}{\tau }_2+B_2\right)},$

(56)

where $B_2$ represents an undetermined complex function of the higher slow time scales compared to ${\tau }_2$ .

By taking into account the solvability conditions mentioned above, we are able to find the unspecified functions $A_j\left(j=1,2,3\right)$ where they are depending only on ${\tau }_2$ . Thus, we express the polar forms of $A_j$ as follows:

$A_j=\frac{1}{2}{\stackrel{\sim }{h}}_j\left({\tau }_2\right)e^{i{\mu }_j\left({\tau }_2\right)},h_j=\epsilon {\stackrel{\sim }{h}}_j\left(j=1,2,3\right),$

(57)

where ${\mu }_j$ and ${\stackrel{\sim }{h}}_j$ denote real functions of the phases and amplitudes for the solutions $y,\varphi ,$ and $\zeta$ , respectively.

Since $A_j=A_j\left({\tau }_2\right)$ , their first-order derivatives can be expressed as follows:

$\frac{d{A}_j}{d\tau }={\epsilon }^2\frac{\partial A_j}{\partial {\tau }_2}.$

(58)

After analyzing the formula presented in equation (58), the PDEs described in system (55) can be converted into a set of ordinary DEs according to the below substitution of the modified phases

${\theta }_j\left({\tau }_2\right)={\tau }_2{\stackrel{\sim }{\sigma }}_j-{\mu }_j\left({\tau }_2\right)j=\left(1,2,3\right).$

(59)

Substituting equations (57)–(59) into (55), and then splitting them into their real and imaginary components, we can get the required system of ordinary DEs in the form

$2{\omega }_1\frac{d{h}_1}{d\tau }=f_1\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_1-{\omega }_1{h}_1\left[c_1+\frac{\alpha \gamma }{\left({\alpha }^2+{\omega }_1^2\right)}\right],$

$2{\omega }_1{h}_1\frac{d{\theta }_1}{d\tau }=f_1\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_1+h_1\left[2{\omega }_1{\sigma }_1-\frac{\gamma {\omega }_1^2}{\left({\alpha }^2+{\omega }_1^2\right)}-\frac{3\left({\omega }_1^2-1\right)}{2\left({\omega }_1^2-4\right)}{h}_2^2\right],$

$2\frac{d{h}_2}{d\tau }=f_2\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_2-c_2{h}_2,$

$2h_2\frac{d{\theta }_2}{d\tau }=f_2\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_2+h_2\left[2{\sigma }_2+\frac{s}{\left(1-{\omega }_2^2\right)}-\frac{3\left({\omega }_1^2-1\right)}{2\left({\omega }_1^2-4\right)}{h}_1^2+\frac{\left({\omega }_1^4+7{\omega }_1^2-8\right)}{8{\omega }_1^2\left({\omega }_1^2-4\right)}{h}_2^2\right],$

$2{\omega }_2\frac{d{h}_3}{d\tau }=f_3\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_3-\frac{c_3{h}_3}{{\omega }_2},$

$2{\omega }_2{h}_3\frac{d{\theta }_3}{d\tau }=f_3\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_3+h_3[2{\omega }_2{\sigma }_3+\frac{s}{\left({\omega }_2^2-1\right)}+\frac{1}{8}{h}_3^2-\frac{\left(1-2{\omega }_2^2\right)}{4{\left(1-{\omega }_2^2\right)}^2}{h}_2^2$

$+\frac{{\omega }_1^4}{2\left(4-{\omega }_1^2{\omega }_2^2\right)}{h}_1^2].$

(60)

The equations of the aforementioned system depict the modulation equations for both modified phases ${\theta }_j$ and amplitudes $h_j$ . Specifically, they comprise six first-order ordinary DEs, representing the occurrence of three resonance cases simultaneously. The solutions to these equations are graphed for various selected values of the system’s parameters, as depicted in Figures (2)–(7), based on the following dataset

$F_0=0.07N,M_1=0.4\mathrm{N}.\mathrm{m},M_2=0.06\mathrm{N}.\mathrm{m},k=50\mathrm{N}/\mathrm{m},$

$c_1=0.05\mathrm{N}.\mathrm{s}/\mathrm{m},c_2=0.06\mathrm{N}.\mathrm{m}.\mathrm{s}/\text{rad},c_3=0.05\mathrm{N}.\mathrm{m}.\mathrm{s}/\text{rad},$

$a=0.4m,b=0.9m,p=0.6m,\Omega =2.5\text{rad}.{\mathrm{s}}^{‐1},{\Omega }_1=2.5\text{rad}.{\mathrm{s}}^{‐1},$

${\Omega }_2=2.5\text{rad}.{\mathrm{s}}^{‐1},{\Omega }_3=1.5\text{rad}.{\mathrm{s}}^{‐1},R=200\Omega ,B=2{\rm T},L_{ind}=0.4{\rm H},$

$L_{coil}=1m,p_1={\omega }_1+{\sigma }_1,p_2=1+{\sigma }_2,p_3=\frac{1}{{\omega }_2}+{\sigma }_3,\mathcal{l}=\left(0.45,0.50,0.55\right)m,$

$m=\left(2.5,2.8,3\right)\text{kg},U_A=\left(4,4.3,4.5\right)\text{kg}.{\mathrm{m}}^2,w_1=\left(4.47,4.23,4.08\right)s^{‐1},$

$w_2=\left(4.67,4.43,4.22\right)s^{‐1},w_3=\left(2.10,2.03,1.98\right)s^{‐1}.$

Figure 2.

Illustrates the variation in the adjusted amplitude $h_1$ versus $\tau$ at: (a) $w_1=\left(4.47,4.23,4.08\right)s^{‐1}$ , (b) $w_2=\left(4.67,4.43,4.22\right)s^{‐1}$ , (c) $w_3=\left(2.10,2.03,1.98\right)s^{‐1}$ .

Figure 3.

Reveals the time history of $h_2$ across different values of $w_j\left(j=1,2,3\right)$ .

Figure 4.

Illustrates the temporal history of $h_3$ at different values of $w_j\left(j=1,2,3\right)$ .

Figure 5.

Illustrates the behavior of ${\theta }_1\left(\tau \right)$ at: (a) $w_1=\left(4.47,4.23,4.08\right)s^{‐1}$ , (b) $w_2=\left(4.67,4.43,4.22\right)s^{‐1}$ , (c) $w_3=\left(2.10,2.03,1.98\right)s^{‐1}$ .

Figure 6.

Showcases the behavior of ${\theta }_2\left(\tau \right)$ at various values of $w_j\left(j=1,2,3\right)$ .

Figure 7.

Displays the behavior of ${\theta }_3\left(\tau \right)$ at various values of $w_j\left(j=1,2,3\right)$ .

These figures explore the influence of system’s parameters on the time behavior of the amplitudes $h_j\left(j=1,2,3\right)$ and modified phases ${\theta }_j$ for various values of $w_j$ . Generally, it must be mentioned that these values are associated with the pendulum’s standard length of the arm $\mathcal{l}$ , the rigid body mass $m$ , and the moment of inertia $U_A$ .

As observed from the depicted waveforms in these figures, they display diminishing patterns over time, suggesting that their behavior is stable. It is seen that $h_1$ has been influenced by the change of $w_1$ and $w_2$ values. The decaying manner of the waves decrease with the increase of $w_1$ and decrease with $w_2$ values, as evidenced in Figure 2(a) and (b), respectively. Conversely, $w_3$ values have no discernible effect on the waves’ behavior, as demonstrated in Figure 2(c). This observation aligns with the mathematical formulas governing the system of equation (60).

Parts of Figure (3) illustrate the changes in $h_2$ overtime as $w_1,w_2,$ and $w_3$ undergo changes. It’s clear that an elevation in $w_1$ values leads to a rise in the vertical extent of the wave represented by $h_2$ , whereas the decay of the waves follows the direction of the minimum value of $w_1$ , as demonstrated in Figure 3(a). Raising $w_2$ values, results in an upsurge of the wave initially drawn near and adjacent to the $h_2$ -axis, followed by a decay in the wave until it stabilizes after the first third of the examined time period, as drawn in Figure 3(b). On the contrary, as depicted in Figure 3(c), an augmentation in $w_3$ values initially causes a reduction in the adjusted oscillations towards the $h_2$ -axis, with the waves eventually attaining stationary behavior at a minimum value of $w_3$ . From this analysis, it can be concluded that the variation of $h_2$ with time exhibits oscillations during the initial third of the time interval being studied, but eventually stabilizes towards the conclusion of the time scale.

The temporal fluctuation of the amplitude $h_3$ is graphically depicted in Figure (4) when considering various values of $w_j$ . It’s observed that the oscillations’ number in the waves describing $h_3$ is fewer when compared to the fluctuations in $h_1$ and $h_2$ . The primary physical rationale for these variations is attributed to the systematic construction of the equations within the system (60). In general, the decaying behavior of the depicted waves in this figure is observed throughout the entire time interval.

The following Figures (5)-(7) illustrate how the modified phases ${\theta }_j$ vary over time $\tau$ as the values of $w_j$ differ. It was observed that altering the values of $w_j$ has a minimal impact on the waves representing the behavior of ${\theta }_1$ , as seen in Figure (5). Conversely, this effect becomes more pronounced as we move to the waves of ${\theta }_2$ and ${\theta }_3$ , as explored in Figures (6) and (7), respectively. Symmetric decaying curves are observed about a parallel horizontal axis to $\tau$ -axis. The adjusted oscillations to the ${\theta }_2$ -axis increase with the increase of $w_1$ and $w_2$ , as indicated in Figure 5(a) and (b), respectively. In a given context of the analysis of Figure 5(c) that is associated with changes in the variable $w_3$ , one can say that the plotted waves exhibit a decreasing trend as the values of $w_3$ increase. It is clear from the Figures (5) and (6) that the intensity of vibrations is greater than in Figure (7), where the reason is due to the mathematical formulas present in the system (60). The waves plotted in Figure (7) behave in a decaying manner over time. In part (a), the amplitudes of the plotted waves and their oscillations number decrease as $w_1$ increases, whereas in parts (b) and (c), one can observe an increasing trend in their amplitudes.

In summary, decaying waves play a crucial role in controlling vibrations, enhancing stability, improving performance, ensuring safety, and reducing noise in vibrating systems.

It’s notable that the patterns exhibited by the curves in Figures (2)-(7) can be interpreted graphically between adjusted amplitudes and phases. This can be done when the temporal parameter is excluded from the relationships that govern them, resulting in spiral curves, as seen in Figures (8)-(10). Each curve converges towards a single point in the planes $h_j{\theta }_j$ , when $w_j$ vary affirming the inherent stability of these waves. Overall, spiral plane curves provide valuable insights into the behavior of vibrating systems, aiding in stability analysis, mode shape representation, frequency response (FR) visualization, and damping analysis.

Figure 8.

Describes the curves in the ${\theta }_1{h}_1$ plane at different values of $w_j\left(j=1,2,3\right)$ .

Figure 9.

Illustrates the curves in the ${\theta }_2{h}_2$ plane at distinct values of $w_j\left(j=1,2,3\right)$ .

Figure 10.

Illustrates the phase curves in ${\theta }_3{h}_3$ plane for different values of $w_j\left(j=1,2,3\right)$ .

The curves depicted in Figures (11)-(13) illustrate the temporal evolution of solutions $z,\phi ,$ and $\psi$ across distinct time intervals, two different time ranges, highlighting the overarching shape of the outcomes. The plotted curves exhibit periodic wave patterns, signifying the stability of the obtained solutions. A detailed examination of Figure (11) reveals that the plotted waves in segments (a) and (d ) manifest as standing waves with discernible nodes. Notably, their amplitudes escalate alongside the increase in $w_1$ values, while maintaining a consistent number of oscillations and wavelengths. Figure 11(b) and (e) present progressive waves evolving over time in response to varying $w_2$ values, where their amplitudes diminish proportionally with the growth of these values. The periodic waves in segments (c) and ( f ) exhibit no alteration with the augmentation of $w_3$ .

Figure 11.

Showcases the efficiency of $w_j\left(j=1,2,3\right)$ values on the solution $z$ .

Figure 12.

Demonstrates the impact of $w_j\left(j=1,2,3\right)$ values on the solution $\phi$ .

Figure 13.

Shows the impact of changes in $w_j\left(j=1,2,3\right)$ on the solution $\psi$ .

A closer examination of the presented waves in Figure 12(a), (d) and (c), (f) shows that the solution $\phi$ does not vary with the change in $w_1$ and $w_3$ values, while it has been varied with the variation of $w_2$ values, as seen in Figure 12(b) and (e). The rationale behind this lies in the mathematical structure of this solution, where it explicitly depends on $w_2$ .

An examination of the graphed waves in Figure (13) shows that how the solution $\psi$ fluctuates in response to changes in $w_j$ values. As the mathematical structure of this solution explicitly relies on $w_j$ values, it’s reasonable to anticipate that it will be influenced by variations in $w_j$ values. Curves of this figure affirm this prediction, showcasing standing waves with discernible nodes are observed in Figure 13(a), (d)) and (c), (f), while progressive waves are noted in segments (b) and (e).

Generally, periodic waves serve as valuable indicators of stability in vibrating dynamical systems, reflecting the delicate balance between forces, resonance phenomena, damping effects, and control strategies. Analyzing the characteristics of these waves provides valuable insights into the behavior and performance of complex systems across various domains.

The phase plane curves corresponding to the achieved solutions are plotted in $z{z}^{\prime },\phi {\phi }^{\prime },$ and $\psi {\psi }^{\prime }$ planes as, depicted in Figures (14)-(16), with variations in $w_1,w_2,$ and $w_3$ . These curves take the form of closed loops, reinforcing the graphed curves of the aforementioned solutions and indicating their stability.

Figure 14.

The phase plane portrait of the solution $z$ at varies values of $w_j\left(j=1,2,3\right)$ .

Figure 15.

Explores the curves in the $\phi {\phi }^{\prime }$ plane at various values of $w_j\left(j=1,2,3\right)$ .

Figure 16.

Showcases the visual representations of the curves in $\psi {\psi }^{\prime }\to$ plane at the aforementioned values of $w_j\left(j=1,2,3\right)$ .

The periodic waves in Figures (17) and (18) illustrate the variation of the generated current $\mathrm{I}$ and power $P$ at various values of $w_j$ during short and large scales of time. Minimal amplitudes of the current generated by the harvesting circuit is demonstrated, leading to low-power levels at the output within the milliwatts range, which is ideal for certain applications. The current exhibits a positive correlation with the change in $w_1$ and $w_2$ values, as evidenced by the presence of standing and progressive waves observed in Figure 17(a), (d) and (b), (e). On the contrary, the current $\mathrm{I}$ shows no variation with the change in $w_3$ values, as noted in segments (c) and (e) of this figure.

Figure 17.

Identifies the chronological progression of the generated current $\mathrm{I}$ at various values of $w_j\left(j=1,2,3\right)$ .

Figure 18.

Illustrates the influences of $w_j\left(j=1,2,3\right)$ on the output power $P$ .

Observing the depicted waves in Figure (18), it’s apparent that the generated power $P$ exhibits periodic behavior, influenced by variations in both $w_1$ and $w_2$ values, with slight alterations noted in response to changes in $w_3$ values.

In the systems energy harvesting, the power output depends on both the current generated by the harvesting device and the voltage at which that current is produced. The relationship between current and power varies depending on the specific characteristics of the energy source, the efficiency of the harvesting technology, and the load connected to the system. Maximizing power output typically involves optimizing the harvesting device’s design and operating conditions to achieve the desired balance between current and voltage.

Solutions at the steady-state

This section addresses the oscillations of the mentioned model when it reaches a steady-state scenario. This means that the transient process of the oscillations disappears due to the system’s damping. In this scenario, equation (60) plays a crucial role in determining the steady-state conditions, where the amplitudes $h_j\left(j=1,2,3\right)$ and modified phases ${\theta }_j$ have zero derivatives.^50,51 Subsequently, we were able to write

$\frac{d{h}_j}{dt}=\frac{d{\theta }_j}{dt}=0.$

(61)

In accordance with these conditions, the below equations are produced

$f_1\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_1-{\omega }_1{h}_1\left[c_1+\frac{\alpha \gamma }{\left({\alpha }^2+{\omega }_1^2\right)}\right]=0,$

$f_1\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_1+h_1\left[2{\omega }_1{\sigma }_1-\frac{\gamma {\omega }_1^2}{\left({\alpha }^2+{\omega }_1^2\right)}-\frac{3\left({\omega }_1^2-1\right)}{2\left({\omega }_1^2-4\right)}{h}_2^2\right]=0,$

$f_2\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_2-c_2{h}_2=0,$

$f_2\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_2+h_2\left[2{\sigma }_2+\frac{s}{\left(1-{\omega }_2^2\right)}-\frac{3\left({\omega }_1^2-1\right)}{2\left({\omega }_1^2-4\right)}{h}_1^2+\frac{\left({\omega }_1^4+7{\omega }_1^2-8\right)}{8{\omega }_1^2\left({\omega }_1^2-4\right)}{h}_2^2\right]=0,$

$f_3\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_3-\frac{c_3{h}_3}{{\omega }_2}=0,$

$f_3\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_3+h_3\left[2{\omega }_2{\sigma }_3+\frac{s}{\left({\omega }_2^2-1\right)}+\frac{1}{8}{h}_3^2-\frac{\left(1-2{\omega }_2^2\right)}{4{\left(1-{\omega }_2^2\right)}^2}{h}_2^2+\frac{{\omega }_1^4}{2\left(4-{\omega }_1^2{\omega }_2^2\right)}{h}_1^2\right]=0.$

(62)

It is noteworthy that the earlier system consists of six algebraic equations involving the variables $h_j$ and ${\theta }_j$ . Upon eliminating ${\theta }_j$ , we directly obtain three equations relating the amplitudes $h_j$ and the detuning parameters ${\sigma }_j$

$f_1^2=4{\omega }_1^2{h}_1^2\{{\left[\frac{c_1}{2}+\frac{\alpha \gamma }{2\left({\alpha }^2+{\omega }_1^2\right)}\right]}^2+{\left[{\sigma }_1-\frac{\gamma {\omega }_1}{2\left({\alpha }^2+{\omega }_1^2\right)}-\frac{3\left({\omega }_1^2-1\right)}{4{\omega }_1\left({\omega }_1^2-4\right)}{h}_2^2\right]}^2\},$

$f_2^2=h_2^2\left\{c_2^2+4{\left[{\sigma }_2+\frac{s}{2\left(1-{\omega }_2^2\right)}-\frac{3\left({\omega }_1^2-1\right)}{4\left({\omega }_1^2-4\right)}{h}_1^2+\frac{\left({\omega }_1^4+7{\omega }_1^2-8\right)}{16{\omega }_1^2\left({\omega }_1^2-4\right)}{h}_2^2\right]}^2\right\},$

$f_3^2=h_3^2\{\frac{c_3^2}{{\omega }_2^2}+4{\omega }_2^2[{\sigma }_3+\frac{s}{2{\omega }_2\left({\omega }_2^2-1\right)}+\frac{1}{16{\omega }_2}{h}_3^2-\frac{\left(1-2{\omega }_2^2\right)}{8{\omega }_2{\left(1-{\omega }_2^2\right)}^2}{h}_2^2$

$+\frac{{\omega }_1^4}{4{\omega }_2(4-{\omega }_1^2{\omega }_2^2)}{h}_1^2{]}^2\}.$

(63)

For the system’s stability analysis, steady-state vibrations are considered essential. To assess this, we focus on the system’s behavior near its equilibrium points. In doing so, we examine the expressions for $h_j$ and ${\theta }_j$ .

$h_j=h_{j0}+h_{j1},$

${\theta }_j={\theta }_{j0}+{\theta }_{j1},j=\left(1,2,3\right).$

(64)

Here, we have the steady-state solutions represented by $h_{j0}$ and ${\theta }_{j0}$ , and the small perturbations are $h_{j1}$ and ${\theta }_{j1}$ which are minor in comparison with $h_{j0}$ and ${\theta }_{j0}$ . By substituting equation (64) into equation (60), we can derive the subsequent system

$2{\omega }_1\frac{d{h}_{11}}{d\tau }=f_1{\theta }_{11}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{10}-{\omega }_1{h}_{11}\left[c_1+\frac{\alpha \gamma }{\left({\alpha }^2+{\omega }_1^2\right)}\right],$

$2{\omega }_1{h}_{10}\frac{d{\theta }_{11}}{d\tau }=-f_1{\theta }_{11}\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_{10}+{\omega }_1{h}_{11}\left[2{\sigma }_1-\frac{\gamma {\omega }_1}{\left({\alpha }^2+{\omega }_1^2\right)}\right]-\frac{3\left({\omega }_1^2-1\right)}{2\left({{\omega }_1}^2-4\right)}\left(2h_{10}{h}_{20}{h}_{21}+h_{11}{h}_{20}^2\right),$

$2\frac{d{h}_{21}}{d\tau }=f_2{\theta }_{21}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{20}-c_2{h}_{21},$

$2h_{20}\frac{d{\theta }_{21}}{d\tau }=-f_2{\theta }_{21}\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_{20}+h_{21}\left[2{\sigma }_2+\frac{s}{\left(1-{\omega }_2^2\right)}+\frac{3\left({\omega }_1^4+7{\omega }_1^2-8\right)}{8{\omega }_1^2\left({\omega }_1^2-4\right)}{h}_{20}^2\right]$

$-\frac{3\left({\omega }_1^2-1\right)}{2\left({\omega }_1^2-4\right)}\left(2h_{10}{h}_{20}{h}_{11}+h_{10}^2{h}_{21}\right),$

$2{\omega }_2\frac{d{h}_{31}}{d\tau }=f_3{\theta }_{31}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{30}-\frac{c_3}{{\omega }_2}{h}_{31},$

$2{\omega }_2{h}_{30}\frac{d{\theta }_{31}}{d\tau }=-f_3{\theta }_{31}\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_{30}+h_{31}\left[2{\omega }_2{\sigma }_3+\frac{s}{\left({\omega }_2^2-1\right)}+\frac{3}{8}{h}_{30}^2\right]-\frac{\left(1-2{\omega }_2^2\right)}{4{\left(1-{\omega }_2^2\right)}^2}$

$\times \left(2h_{20}{h}_{21}{h}_{30}+h_{20}^2{h}_{31}\right)+\frac{{\omega }_1^4}{2\left(4-{\omega }_1^2{\omega }_2^2\right)}\left(2h_{10}{h}_{11}{h}_{30}+h_{10}^2{h}_{31}\right).$

(65)

The functions $h_{j1}$ and ${\theta }_{j1}$ are defined previously as perturbed functions for the amplitudes and modified phases. Thereafter, their solutions can be formulated as a linear combination of eigenvalues $\lambda$ and constants ${\hslash }_{\delta }$ in the form ${\hslash }_{\delta }{e}^{\lambda \tau }\left(\delta =1,2,...,6\right)$ . The real parts of the roots of the characteristic equation (66) should be negative for the steady-state solutions $h_{j1}$ and ${\theta }_{j1}$ to be asymptotically stable

${\lambda }^6+{\Gamma }_1{\lambda }^5+{\Gamma }_2{\lambda }^4+{\Gamma }_3{\lambda }^3+{\Gamma }_4{\lambda }^2+{\Gamma }_5\lambda +{\Gamma }_6=0,$

(66)

where ${\Gamma }_{\delta }\left(\delta =1,2,...,6\right)$ can be obtained easily (see Appendix II).

We will now outline the fundamental stability criteria that align with the RHC ⁵² for certain solutions in the scenario of steady-state as in the below forms

${\Gamma }_1>0,{\Gamma }_1{\Gamma }_2-{\Gamma }_3>0,$

${\Gamma }_3\left({\Gamma }_1{\Gamma }_2-{\Gamma }_3\right)+{\Gamma }_1\left({\Gamma }_5-{\Gamma }_1{\Gamma }_4\right)>0,$

${\Gamma }_1\left({\Gamma }_2{\Gamma }_3{\Gamma }_4-{{\Gamma }_2}^2{\Gamma }_5+2{\Gamma }_4{\Gamma }_5-{\Gamma }_3{\Gamma }_6\right)+{\Gamma }_3\left({\Gamma }_2{\Gamma }_5-{\Gamma }_3{\Gamma }_4\right)+{{\Gamma }_1}^2\left({\Gamma }_2{\Gamma }_6-{{\Gamma }_4}^2\right)-{{\Gamma }_5}^2>0,$

$-{\Gamma }_5\left[{\Gamma }_4\left(-{\Gamma }_1{\Gamma }_2{\Gamma }_3+{{\Gamma }_3}^2+{{\Gamma }_1}^2{\Gamma }_4\right)+{\Gamma }_1\left({{\Gamma }_2}^2-2{\Gamma }_4\right){\Gamma }_5-{\Gamma }_2{\Gamma }_3{\Gamma }_5+{{\Gamma }_5}^2\right]$

$+\left[{{\Gamma }_3}^3-{\Gamma }_1{\Gamma }_3\left({\Gamma }_2{\Gamma }_3+3{\Gamma }_5\right)+{{\Gamma }_1}^2\left({\Gamma }_3{\Gamma }_4+2{\Gamma }_2{\Gamma }_5\right)\right]{\Gamma }_6-{{\Gamma }_1}^3{{\Gamma }_6}^2>0,$

${\Gamma }_6\{-{\Gamma }_5\left[{\Gamma }_4\left(-{\Gamma }_1{\Gamma }_2{\Gamma }_3+{{\Gamma }_3}^2+{{\Gamma }_1}^2{\Gamma }_4\right)+{\Gamma }_1\left({{\Gamma }_2}^2-2{\Gamma }_4\right){\Gamma }_5-{\Gamma }_2{\Gamma }_3{\Gamma }_5+{{\Gamma }_5}^2\right]$

$+\left[{{\Gamma }_3}^3-{\Gamma }_1{\Gamma }_3\left({\Gamma }_2{\Gamma }_3+3{\Gamma }_5\right)+{{\Gamma }_1}^2\left({\Gamma }_3{\Gamma }_4+2{\Gamma }_2{\Gamma }_5\right)\right]{\Gamma }_6-{{\Gamma }_1}^3{{\Gamma }_6}^2\}>0.$

(67)

There are numerous important advantages of applying the RHC to a vibrating dynamic system, especially when it relates to stability analysis. Here are a few main benefits:

Stability determination

The main advantage is that the RHC offer a methodical way to assess the system’s stability without having to directly solve for the characteristic equation’s roots. This is especially helpful for higher-order systems, where locating the roots can be difficult and time-consuming.

Comprehensive analysis

It makes it possible to analyze every root of the characteristic equation simultaneously, providing a comprehensive view of the stability of the system. This is important because an unstable system might result from even one root having a positive real part.

Predictive power

By using the RHC, one can predict changes in stability due to parameter variations. This is useful when designing and managing vibrating systems whose characteristics could alter over time or as a result of outside factors.

Overall, the RHC are valuable tools in the analysis and design of vibrating dynamic systems, offering a robust, efficient, and insightful means to assess and ensure stability.

The analysis of stability

This section delves into the examination of the system’s stability under consideration in utilizing the RHC of nonlinear stability. The analysis involves determining the stability parameters and conducting simulations of the system. The parameters of detuning ${\sigma }_1,{\sigma }_2,$ and ${\sigma }_3$ along with the frequencies $w_1,w_2,$ and $w_3$ are crucial in the stabilization process. Based on chosen parameter values, stability plots for system (65) are produced. The adjusted amplitudes $h_1,h_2,$ and $h_3$ are drawn across distinct parametric regions, showcasing their characteristics through the trajectories of phase plane.

The potential stability and instability regions of FPs are illustrated in Figures (19)-(21) with varying values of ${\sigma }_1$ , showing that this range lies in the interval ${\sigma }_1\in \left[-0.2,0.2\right]$ . A closer look at Figure 19(a) reveals that at $w_1\left(=4.47,4.23\right)$ , the stable FPs fall within the area ${\sigma }_1\le 0.01,$ while the unstable zones of these points exist in the region $0.01<{\sigma }_1$ . Additionally, for $w_1=4.08$ , the stable areas of FPs span to be in the range ${\sigma }_1\le 0$ , while the unstable ones are found in the region of $0<{\sigma }_1$ . The drawn solid curve and dashed one depict the stable and unstable zones of FPs, respectively.

Figure 19.

Presents the FR curves in the plane ${\sigma }_1{h}_1$ at ${\sigma }_2=p_2-1,$ and ${\sigma }_3=p_3-\left(\frac{1}{{\omega }_2}\right)$ .

Figure 20.

Examines the fluctuation of $h_2$ via ${\sigma }_1$ at ${\sigma }_2=p_2-1,$ and ${\sigma }_3=p_3-\left(\frac{1}{{\omega }_2}\right)$ .

Figure 21.

Shows the FR curves in the plane ${\sigma }_1{h}_3$ at ${\sigma }_2=p_2-1,$ and ${\sigma }_3=p_3-\left(\frac{1}{{\omega }_2}\right)$ when: (a) $w_1=\left(4.47,4.23,4.08\right)s^{‐1}$ , (b) $w_2=\left(4.67,4.43,4.22\right)s^{‐1}$ , (c) $w_3=\left(2.10,2.03,1.98\right)s^{‐1}$ .

Figure 19(b) clearly demonstrates that the FPs maintain stability within the specified range of ${\sigma }_1\le 0$ , yet their stability loses once ${\sigma }_1$ surpasses $0$ . Figures 19(c) shows that at $w_3=2.10$ , the stable zones are found in the range ${\sigma }_1\le 0$ , while the unstable ones are existed in the range $0<{\sigma }_1$ , with a lying critical FP between the two regions. At $w_3\left(=2.03,1.98\right)$ , the stable FPs will be in the interval ${\sigma }_1\le 0.01$ , and the unstable points are within $0.01<{\sigma }_1$ . Based on the analysis of Figure (19), one can conclude that the amplitudes $h_2$ and $h_3$ via ${\sigma }_1$ demonstrate analogous patterns of stability and instability. These patterns are further illustrated in Figures (20) and (21). These figures are calculated at ${\sigma }_2=p_2-1,$ and ${\sigma }_3=p_3-\left(\frac{1}{{\omega }_2}\right)$ . The critical values of the presented wave, marking the occurrence of transcritical bifurcation due to a FP, can be determined by the intersections of these curves.

Figures (22)-(24) illustrate the changes in FR curves of $h_1,h_2,$ and $h_3$ through the detuning parameter ${\sigma }_2$ for various values of $w_j\left(j=1,2,3\right)$ ranging from $-0.2$ to $0.4$ . The purpose of these graphs is to showcase the impact of $w_j$ on the FR curves of $h_j$ in relation to ${\sigma }_2$ at ${\sigma }_1=p_1-{\omega }_1,$ and ${\sigma }_3=p_3-\left(\frac{1}{{\omega }_2}\right)$ . Each curve has a single critical FP that delineates the stable and unstable regions at $w_1\left(=4.47,4.23,4.08\right)$ , as depicted in Figure 22(a). The stable region ranges from $-0.2$ to $0.13$ for ${\sigma }_2$ , while the unstable region ranges from $0.13$ to $0.4$ . Figure 22(b) reveals a critical point and three peaks when $w_2=4.67$ , with stability and instability regions falling between $-0.2$ to $0.1$ and $0.1$ to $0.4$ for ${\sigma }_2$ respectively. Similarly, at $w_2=4.43$ , the stable region ranges from $-0.2$ to $0.13$ for ${\sigma }_2$ , while the unstable region ranges from $0.13$ to $0.4$ . At $w_2=4.22$ , the stable region is between $-0.2$ to $0.16$ for ${\sigma }_2$ , while the unstable region is between $0.16$ to $0.4$ .

Figure 22.

Represents the FR curves in the plane ${\sigma }_2{h}_1$ at ${\sigma }_1=p_1-{\omega }_1,$ and ${\sigma }_3=p_3-\left(\frac{1}{{\omega }_2}\right)$ when: (a) $w_1=\left(4.47,4.23,4.08\right)s^{‐1}$ , (b) $w_2=\left(4.67,4.43,4.22\right)s^{‐1}$ , (c) $w_3=\left(2.10,2.03,1.98\right)s^{‐1}$ .

Figure 23.

Shows the change in $h_2$ versus ${\sigma }_2$ at ${\sigma }_1=p_1-{\omega }_1,$ and ${\sigma }_3=p_3-\left(\frac{1}{{\omega }_2}\right)$ .

Figure 24.

Depicts how the modified amplitude $h_3$ varies as a function of the detuning parameter ${\sigma }_2$ at ${\sigma }_1=p_1-{\omega }_1,$ and ${\sigma }_3=p_3-\left(\frac{1}{{\omega }_2}\right)$ .

In Figure 22(c), a single critical FP is observed at $w_3=2.10$ , while the regions of stability and instability are found at ${\sigma }_2\le 0.13$ and ${\sigma }_2>0.13$ , respectively. On the other hand, at $w_3\left(=2.03,1.98\right)$ , the stable region zone are located at ${\sigma }_2\le 0.11$ , while the unstable one is found in the range ${\sigma }_2>0.11$ along with multiple peaks of FPs. Ultimately, both $h_2$ and $h_3$ display analogous stability and instability patterns concerning ${\sigma }_2$ , as depicted in Figures (23) and (24).

The adjusted magnitudes $h_j$ are plotted against ${\sigma }_3$ in Figures (25)-(27) to demonstrate the significant impact of $w_j$ on the dynamics of the system at ${\sigma }_1=p_1-{\omega }_1,$ and ${\sigma }_2=p_2-1$ . The analysis focuses on the stability and instability characteristics within the range $-0.2\le {\sigma }_3\le 0.2$ . The behavior of the FR curves, representing the variation of $h_1$ via ${\sigma }_3$ at different values of $w_j$ values, is depicted as straight lines in Figure (25). According to the calculations in part (a) of Figure (25), (a) single critical FP exists, with stable and unstable FP located in the intervals ${\sigma }_3\le -0.06,{\sigma }_3\le -0.07,$ ${\sigma }_3\le -0.08,$ and $-0.06<{\sigma }_3,$ $-0.07<{\sigma }_3,-0.08<{\sigma }_3$ , respectively. When $w_2\left(=4.67,4.43\right)$ , as in Figure 25(b), the FPs are stable within the range ${\sigma }_3\le -0.08,$ but lose stability in the region $-0.08<{\sigma }_3$ . Similarly, at $w_2=4.22$ , stable and unstable regions occur in ${\sigma }_3\le -0.09,$ and $-0.09<{\sigma }_3,$ respectively. At $w_3=2.10$ , as seen in Figure 25(c), the FPs are stable in the interval ${\sigma }_3\le -0.08,$ while these points lose their stabilities in the range $-0.08<{\sigma }_3$ . Furthermore, when $w_3\left(=2.03,1.98\right)$ , only one critical and one peak FPs are observed, with stable regions existing in ${\sigma }_3\le -0.07$ and unstable regions in $-0.07<{\sigma }_3$ , respectively. Figures (26) and (27) illustrate the FR curves of $h_2$ and $h_3$ for various values of $w_j$ . It is important to note that the stability and instability regions depicted in these figures are consistent with those in segments of Figure (25).

Figure 25.

Shows the FR curves in the plane ${\sigma }_3{h}_1$ at ${\sigma }_1=p_1-{\omega }_1,$ and ${\sigma }_2=p_2-1$ .

Figure 26.

Illustrates the FR curves for $h_2$ via ${\sigma }_3$ at ${\sigma }_1=p_1-{\omega }_1,$ and ${\sigma }_2=p_2-1$ .

Figure 27.

Displays the frequency characteristics of $h_3$ via ${\sigma }_3$ at: (a) $w_1=\left(4.47,4.23,4.08\right)s^{‐1}$ , (b) $w_2=\left(4.67,4.43,4.22\right)s^{‐1}$ , and (c) $w_3=\left(2.10,2.03,1.98\right)s^{‐1}$ .

Nonlinear analysis

The investigation in this section focuses on the stability of a dynamical model utilizing nonlinear amplitudes $A_j\left(j=1,2,3\right)$ . To examine the behavior of the nonlinear amplitude in system (55) and evaluate its stability, the following transformations will be introduced

$A_1=\frac{1}{2}\left({\stackrel{\sim }{U}}_1\left({\tau }_2\right)+i{\stackrel{\sim }{V}}_1\left({\tau }_2\right)\right)e^{i{\stackrel{\sim }{\sigma }}_1{\tau }_2},$

$A_2=\frac{1}{2}\left({\stackrel{\sim }{U}}_2\left({\tau }_2\right)+i{\stackrel{\sim }{V}}_2\left({\tau }_2\right)\right)e^{i{\stackrel{\sim }{\sigma }}_2{\tau }_2},$

$A_3=\frac{1}{2}\left({\stackrel{\sim }{U}}_3\left({\tau }_2\right)+i{\stackrel{\sim }{V}}_3\left({\tau }_2\right)\right)e^{i{\stackrel{\sim }{\sigma }}_3{\tau }_2},$

$U_j=\epsilon {\stackrel{\sim }{U}}_j,V_j=\epsilon {\stackrel{\sim }{V}}_j,{\sigma }_j={\epsilon }^2{\stackrel{\sim }{\sigma }}_j;\left(j=1,2,3\right).$

(68)

Through the application of (55) and (68) and the subsequent isolation of the real and imaginary terms, we get

$f_1+{\omega }_1\left(2\frac{d{V}_1}{d\tau }+2{\sigma }_1{U}_1+c_1{V}_1\right)-\frac{3\left({\omega }_1^2-1\right)}{2\left({\omega }_1^2-4\right)}{U}_1\left(U_2^2+V_2^2\right)$

$+\frac{{\omega }_1\gamma }{\left({\alpha }^2+{\omega }_1^2\right)}\left(\alpha V_1-{\omega }_1{U}_1\right)=0,$

${\omega }_1\left(2{\sigma }_1{V}_1-2\frac{d{U}_1}{d\tau }-c_1{U}_1\right)-\frac{3\left({\omega }_1^2-1\right)}{2\left({\omega }_1^2-4\right)}{V}_1\left(U_2^2+V_2^2\right)-\frac{{\omega }_1\gamma }{\left({\alpha }^2+{\omega }_1^2\right)}\left(\alpha U_1+{\omega }_1{V}_1\right)=0,$

$f_2+2\left(\frac{d{V}_2}{d\tau }+{\sigma }_2{U}_2\right)+c_2{V}_2+U_2[\frac{s}{\left(1-{\omega }_2^2\right)}-\frac{3\left({\omega }_1^2-1\right)}{2\left({\omega }_1^2-4\right)}\left(U_1^2+V_1^2\right)$

$+\frac{\left({\omega }_1^4+7{\omega }_1^2-8\right)}{8{\omega }_1^2\left({\omega }_1^2-4\right)}\left(U_2^2+V_2^2\right)]=0,$

$2\left({\sigma }_2{V}_2-\frac{d{U}_2}{d\tau }\right)-c_2{U}_2+V_2[\frac{s}{\left(1-{\omega }_2^2\right)}-\frac{3\left({\omega }_1^2-1\right)}{2\left({\omega }_1^2-4\right)}\left({U_1}^2+{V_1}^2\right)$

$+\frac{\left({\omega }_1^4+7{\omega }_1^2-8\right)}{8{\omega }_1^2\left({\omega }_1^2-4\right)}\left(U_2^2+V_2^2\right)]=0,$

$f_3+2{\omega }_2\left(\frac{d{V}_3}{d\tau }+{\sigma }_3{U}_3\right)+\frac{c_3{V}_3}{{\omega }_2}+U_3[\frac{s}{\left({\omega }_2^2-1\right)}+\frac{1}{8}\left(U_3^2+V_3^2\right)$

$-\frac{\left(1-2{\omega }_2^2\right)}{4{\left(1-{\omega }_2^2\right)}^2}\left(U_2^2+V_2^2\right)+\frac{{\omega }_1^4}{2\left(4-{\omega }_1^2{\omega }_2^2\right)}\left(U_1^2+V_1^2\right)]=0,$

${\sigma }_3{V}_3-\frac{d{U}_3}{d\tau }-\frac{1}{2{\omega }_2^2}{c}_3{U}_3+V_3[\frac{s}{2{\omega }_2\left({\omega }_2^2-1\right)}+\frac{1}{16{\omega }_2}\left(U_3^2+V_3^2\right)$

$-\frac{\left(1-2{\omega }_2^2\right)}{8{\omega }_2{\left(1-{\omega }_2^2\right)}^2}\left(U_2^2+V_2^2\right)+\frac{{\omega }_1^4}{4{\omega }_2\left(4-{\omega }_1^2{\omega }_2^2\right)}\left(U_1^2+V_1^2\right)]=0.$

(69)

Throughout the study, the adjusted amplitudes were validated in a range of parametric areas, and their characteristics have been plotted on the phase plane curves, as demonstrated in Figures (28)-(36), utilizing the earlier data.

Figure 28.

Shows the decaying behavior of the function $U_1\left(\tau \right)$ .

Figure 29.

Presents the decaying behavior of the function $U_2\left(\tau \right)$ .

Figure 30.

Describes the time history of the function $U_3$ .

Figure 31.

Describes the variation of the function $V_1\left(\tau \right)$ .

Figure 32.

Shows $V_2\left(\tau \right)$ variation at various values of $w_j$ .

Figure 33.

The time behavior of $V_3$ .

Figure 34.

Shows the trajectory in phase plane $U_1{V}_1$ when $w_j$ vary.

Figure 35.

Describes the paths in the plane $U_2{V}_2$ when $w_j$ vary.

Figure 36.

Illustrates the trajectories in the plane $U_3{V}_3$ when $w_j$ vary.

The presented waves in Figures (28)-(33) illustrate the temporal evolution of the phases $U_j$ and $V_j$ via $\tau$ . Thus, it can be observed that the time histories of such phases undergo fluctuations, first oscillating between growth and decline, and then stabilizing towards the end of the time period. In other words, these waves have decaying form over the examined time, where the amplitudes of these waves decrease and increase with the raise of $w_1$ and $w_2$ , respectively, as explored in Figure 28(a) and (b). Whereas there is no change in the waves’ amplitude when $w_3$ varies, as seen in Figure 28(c). This behavior signifies the stability of the modified amplitudes throughout the duration of the study.

One may observe an analogous behavior for the waves describing the function $U_2\left(\tau \right)$ at various values of $w_2$ , as in Figure (29), with the decay of these waves being greater than those depicted in Figure (28). This implies that these waves reach a stable stage more swiftly than before. Segments of Figure (30) show decaying curves for countable waves as the adjusted phase $U_j$ varies with different $w_j$ values. The waves illustrated in Figures (31)-(32) for the functions $V_j\left(\tau \right)$ demonstrate similar variations to those seen in Figures (28)-(30) when $w_1,w_2,$ and $w_3$ at their specified values.

Figures (34)-(36) depict the projections of the presented trajectories in Figures (28)-(33) onto phase planes $U_j{V}_j$ after removing the dimensionless time $\tau$ . The drawn paths in these planes have spiral shapes that converge towards a single point which showcases the stable behavior of the system under investigation.

Conclusion

This study has effectively addressed the dual objectives of reducing vibrations and harnessing energy in a 3DOF dynamical model. This model consists of a damped spring-pendulum with its pivot moving in an elliptical path at a constant angular velocity, coupled with a rigid body experiencing external harmonic forces and moments. The optimization of energy harvesting and vibration reduction has been achieved in the study by employing an electromagnetic system that utilizes the movement of a magnet within a coil. Novel, precise solutions up to the third approximation were derived in the analysis through the use of the second kind of Lagrange’s equations and the MST. Multiple resonance states were identified and scrutinized at the same time. The time evolution of solutions, along with modified phases and amplitudes, was detailed through comprehensive graphical representations, highlighting the major parameters that influence system behavior. The system’s consistent dynamics were depicted through phase portraits, and the power output and current time profiles of the electromagnetic device underscored the influence of various parameters. The stability analysis using RHC revealed a wide range of stable parameter configurations, highlighting the system’s robustness. The model’s relevance spans contemporary applications in control sensors, offering benefits across various sectors including industry, construction, automotive, transportation, and infrastructure. The contributions of this study provide important insights for enhancing the efficiency and stability of vibration reduction and energy harvesting systems in real-world scenarios.