The stability analysis of a dynamical system equipped with a piezoelectric energy harvester device near resonance

Abstract

In this work, the influence of a piezoelectric device on the planar motion of a two-degrees-of-freedom (DOF) dynamical system is examined. This system comprises a nonlinear damping spring pendulum (SP), whose pivot point moves on a Lissajous curve at a constant angular velocity in an anticlockwise direction amidst the influence of external forces and moments. The second kind of Lagrange’s equations is used to derive the governing system of equations of motion (EOM), which is transformed into a dimensionless form. Up to the third approximation, the novel asymptotic solutions of this system are obtained using the multiple scales technique (MST). The numerical solutions of the EOM are obtained applying the fourth-order Runge–Kutta method. The comparison between both solutions, the asymptotic solutions and the numerical ones reveals great agreement between them, which in turn indicates the accuracy of the used perturbation technique. The solvability criteria are established once the generated secular terms are eliminated, and all resonance cases are subsequently categorized. In light of the system’s adjusted phases, the modulation equations for two of these cases are simultaneously investigated. Graphically, numerous plots depicting the time histories of the achieved solutions are examined, and the nonlinear stability analyses of the modulation equations are investigated in accordance with the steady-state solutions. The outputs of the energy harvesting (EH) device are scrutinized to investigate and analyze the influences of damping coefficients, excitation amplitudes, and different frequencies. Multiple regions of stability and instability are delineated within the frequency response curves (FRC), illustrating the impact of varying parameters on the dynamic behavior of the system. The obtained results are considered significant due to their practical applications in our lives. It can be utilized for medical applications, charging electrical gadgets, powering sensors, smartphones, and other gadgets. It also powers keyless entry systems, patient monitors, airbag sensors, fish and depth finders, and audible alarms like smoke alarms.

Introduction

Piezoelectric EH devices have become a popular way to capture energy from ambient vibrations in dynamic systems. The potential uses of this confluence of mechanical motion and electrical energy conversion in powering wireless networks, sensors, and small-scale electronic equipment in situations where conventional power sources are unavailable or unfeasible have generated a great deal of attention. A fascinating feature of this field is the complex interaction between the piezoelectric device’s EH properties and the motion of the underlying dynamic system. Research initiatives aiming at comprehending and enhancing the performance of such systems are built on this interconnectedness.

In this regard, a crucial field of research is the motion analysis of a coupled vibrating dynamical system using a piezoelectric EH device. This multidisciplinary endeavor explores the intricate dynamics that dictate the coupled behavior of piezoelectric transducers and mechanical structures. By means of theoretical examination, numerical models, and empirical confirmations, scientists strive to elucidate the complexities of interconnected oscillating dynamic systems enhanced by piezoelectric energy harvesters.

Over the past several years, there has been an increase in interest in transforming the vibrational motion of various dynamical systems into useful electrical energy using harvesting devices.^1–3 Applications of EH were studied in Refs 4,5 and include distributed wireless sensor nodes for monitoring the condition of structures, embedded and implanted sensor nodes for medical applications, load batteries for large systems, pressure monitoring of cars, power supply for unmanned vehicles, and home security systems. Now, wireless electronic devices or systems are growing in popularity ⁶ because they don’t require a power connection. Piezoelectric materials have been the subject of extensive research to create simple yet efficient EH devices. ⁷ The literature on piezoelectric EH contains a wide range of models ^8–14 to represent the behavior of harvesters connected electromechanically. These models explore the viability of using vibrations as a power source in situations where vibrations are prevalent, without arguing that it is the best or most flexible way to scavenge ambient power. ¹⁵ An overview of earlier work and design aspects for piezoelectric-based EH for MEMS-scale sensors are given in Ref 16. Harvested ambient vibration energy may be used to meet the power needs of sophisticated MEMS-scale autonomous sensors for a variety of applications, including structural health monitoring.

In Ref 17, the recommended energy harvesters under harmonic base excitations from the theory of flexoelectricity are developed using the Hamilton’s principle. The authors offered a simple device for collecting energy from oscillating surrounding environment that is based on spring oscillations. To obtain solutions that fit generalized coordinates, vibrating mathematical models are investigated in Refs 18–23. The stability of stable fixed points is verified, and the Routh-Hurwitz nonlinear stability analysis is examined in Ref 18. In Ref 19, the nonlinear dynamics of low-frequency stimulation in broadband piezoelectric energy harvesters is investigated. The fractional calculus is used and the impact of damping, excitation amplitude, and frequency on dynamic behaviors are examined. Chaos and periodic motion are all displayed by the system when the fractional order is changed. Whereas in Ref 20, two dynamical theories for transforming an SP vibrational motion into electrical energy were studied. Piezoelectric and electromagnetic EH devices were coupled independently. The primary controlling system was calculated using Lagrange’s equations, and its solutions were attained. The stability of fixed points was explored, as well as external resonance instances. The positive impact of excitation amplitudes and damping coefficients on the output voltage, current, and power is determined. The dynamical movements of damped SPs with 2DOF or 3DOF that are subjected to excitation forces and moments are examined in Refs 21–23. A general vibrating system of a three DOF double pendulum was investigated in Refs 21,22 wherein the pivot point follows an elliptic trajectory with a stationary angular velocity. Moreover, in Ref 23, the movement of this point was restricted to be on a Lissajous for the issue of a harmonically stimulated damped SP. To assess the accuracy of the analytical solutions, they were obtained up to a higher approximation and compared to the NS.

Additionally, a rigid rod is thought to relate with the pendulum’s motion, in which two separates massless upstretched rigid rods are used in Ref 24 to examine the dynamical behavior of a 2DOF double pendulum system. The motion of a 3DOF dynamical system consisting of three pendulums in a plane is examined in Refs 25,26 in which two of them are rigid bodies that are attached to another massless rigid arm, ²⁵ while a triple rigid body pendulum is investigated in. Ref 26. The movement of an unusual 3DOF auto-parametric dynamical system consisting of a forced nonlinear damped Duffing oscillator as the primary system and a damped SP as the secondary system is examined in. Ref 27,28 In Ref 29, the planar motion of an attached 2DOF damped system with an auto-parametric model is investigated. The extension of this research has been examined in Ref 30 for the case of an SP. The motion of a 3DOF dynamical system composed of a nonlinear Duffing oscillator linked to a nonlinear damped SP is investigated in Ref 31 It is taken into mind that this system is connected to an electromagnetic device of EH. It is found that at low frequencies, the EH devices produce large output power and current. In Ref 32, the authors focused their study on the motion of a novel 3DOF system comprised of two connected segments. A piezoelectric transducer is attached to this model to transform the resulting vibrational motion into electrical power. The bifurcation plots and the Lyapunov exponent diagrams are drawn to show the impact of various parameters on the system’s behavior.

In recent years, the study of nonlinear dynamical systems and their interactions with various physical phenomena has gained significant attention. The research has focused on understanding and harnessing these complex interactions for practical applications, such as energy harvesting and vibration control. In the following, some recent studies have been selected for these systems. In Ref. 33, the authors explored the innovative concept of utilizing vortex-induced vibrations to extract energy using piezoelectric materials. The internal resonance mechanisms within a cylinder subjected to these vibrations are investigated, aiming to enhance the efficiency of energy extraction. A detailed analysis is provided of the dynamic interactions and their potential for practical energy harvesting applications. In Ref. 34, the authors delved into the performance and dynamical behaviors of piezoelectric energy harvesters influenced by dissipative viscoelastic impacts. The nonlinear dynamics and the intricate interactions between the harvester’s components are emphasized. By providing insights into the system’s performance under various conditions, the optimization of energy harvesting devices has been improved. The complex dynamical interactions between a double pendulum system and nonlinear sloshing in cranes carrying a liquid tank are investigated in Ref. 35. It is aimed at understanding the coupled oscillations and their impact on the stability and control of the cranes. Such insights are crucial for designing safer and more efficient crane operations in industrial and maritime settings. The nonlinear static and dynamic behaviors of rod pendulums that are either partially or fully submerged in still water are examined in Ref. 36. This work provided valuable knowledge about the effects of submersion on the pendulum’s motion, contributing to a broader understanding of fluid-structure interactions and their implications for marine engineering applications. The analysis of a time-delayed excited damping Duffing oscillator is examined in Ref. 37. The complex dynamics introduced by time delays and their effects on the system’s stability and response are explored. Such investigations are pertinent for the design of advanced damping systems in mechanical and structural engineering. The use of active dampers to mitigate multistability vibrations in a 2DOF rotor model is examined in Ref. 38. The challenges posed by simultaneous multiparametric and external harmonic excitations are addressed, offering solutions to enhance the rotor’s stability and performance in various engineering applications. The vibrational dynamics of a system subjected to external torque and excitation forces are investigated in Ref. 39. This research aims to provide a comprehensive understanding of how these external influences affect the system’s behavior, with implications for the design and control of mechanical systems exposed to dynamic loading conditions. The stability of a 3DOF vibrating system in the vicinity of resonances is explored in Ref. 40. By examining the system’s behavior near resonance conditions, the authors provide insights that are critical for predicting and mitigating instability in engineering systems subjected to vibrational excitations.

In this study, a novel dynamical system of a SP is investigated for converting its vibrational movement into electricity by connecting a piezoelectric energy harvester device to this system. The main EOM are determined using Lagrange’s equations from second kind, and the piezoelectric circuit’s mechanism is employed to construct the equivalent equation. New analytical solutions of the system’s EOM up to the third approximation are achieved utilizing the MST, and their accuracies are verified by comparison with the numerical solutions using the aforementioned numerical method of the EOM. The conditions of solvability are established once the secular terms are removed, and all resonance cases are then categorized. The equations of modulation for two of these cases are simultaneously assessed in view of the system’s altered phases. Several plots showing the time histories of the obtained solutions are presented graphically using the Wolfram Mathematica 12.1. The stability analyses of the ME are examined in accordance with the solutions at the steady-state scenario. Moreover, the stability of the nonlinear amplitudes within the modulation equations system and the attributes of nonlinear analysis have been showcased through the utilization of complex transformations applied to the system’s amplitudes. The piezoelectric device is affected by the motion of the main dynamical system, so we must study the behavior of the main system. The impacts of damping coefficients, excitation amplitudes, and different frequencies are investigated and analyzed with regard to the outputs of EH devices, such as voltage and power. To illustrate how the acted-upon parameters affect the system’s behavior, a number of stable and unstable regions for the FRC are drawn. The acquired results are significant because they have real-world applicability in our daily lives. It can be used to power sensors, charge electrical devices, and perform medical procedures.

The dynamic approach modelling

This section introduces the required investigated system, in which the EOM can be obtained in main form and the corresponding dimensionless ones, besides the piezoelectric circuit equation. Therefore, let’s have a look at the motion of a 2DOF dynamical system that includes a nonlinear SP that has a linear stiffness coefficient $k_1$ and a nonlinear one with coefficient $k_2$ , respectively. The point of suspensions $O$ is limited to travel along a Lissajous curve with angular velocity $\Omega$ in the direction of anticlockwise direction with coordinates $X_O=R_X\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{X}t,$ and $Y_O=R_Y\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\Omega }_{Y}t$ , where $R_X,R_Y,{\Omega }_X,$ and ${\Omega }_Y$ are known parameters.

To follow the description of the resented system, we take into account the action of two harmonic external forces, $F_1\left(t\right)=F_1\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{1}t$ on the mass $m$ along the spring direction and $F_2\left(t\right)=F_2\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{2}t$ in the direction of increasing the rotation angle $\theta$ , in addition to the moment $M_{\theta }\left(t\right)=M_O\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{3}t$ acting at $O$ and $X\left(t\right)$ represents the spring elongation, as seen in Figure 1. Here, $\left({\Omega }_1,{\Omega }_2,{\Omega }_3\right)$ and $\left(F_1,F_2,M_O\right)$ are the frequencies and amplitudes of $F_1\left(t\right),F_2\left(t\right),$ and $M_{\theta }\left(t\right)$ , respectively. Let $C_1$ and $C_2$ denote the coefficients of damping of the damping force $C_1\dot{X}$ and damping moment $C_2\dot{\theta }$ , respectively, $l_0$ is the pendulum’s normal length and $g$ is the gravity’s acceleration. The resistive load of the piezoelectric circuit $R_p$ will now be taken into consideration in the connection of a piezoelectric device. Let $\gamma$ points to the linear coupling coefficient for the piezoelectric circuit, and $c_p$ signifies the piezoelectric capacitance.

Figure 1.

The dynamical system with piezoelectric EH devices.

We may express the mass’s Cartesian coordinates as follows

$X_A=R_X\hspace{0.17em}\mathit{cos}\left({\Omega }_{X}t\right)+\left(l_0+X\right)\mathit{cos}\hspace{0.17em}\theta ,Y_A=R_Y\hspace{0.17em}\mathit{sin}\left({\Omega }_{Y}t\right)+\left(l_0+X\right)\mathit{sin}\hspace{0.17em}\theta .$

(1)

The dynamical system’s potential and kinetic energy may now be expressed as

$V=-mg\left[R_X\hspace{0.17em}\mathit{cos}\left({\Omega }_{X}t\right)+\left(l_0+X\right)\mathit{cos}\hspace{0.17em}\theta \right]+\frac{1}{2}{k}_1{X}^2+\frac{1}{4}{k}_2{X}^4,T=\frac{1}{2}m\left\{{\left[R_Y{\Omega }_Y\hspace{0.17em}\mathit{cos}\left({\Omega }_{Y}t\right)+\dot{X}\mathit{sin}\hspace{0.17em}\theta +\left(l_0+X\right)\dot{\theta }\mathit{cos}\hspace{0.17em}\theta \right]}^2+{\left[\dot{X}\mathit{cos}\hspace{0.17em}\theta -R_X{\Omega }_X\hspace{0.17em}\mathit{sin}\left({\Omega }_{X}t\right)-\left(l_0+X\right)\dot{\theta }\mathit{sin}\hspace{0.17em}\theta \right]}^2\right\},$

(2)

where the differentiation over time $t$ is represented by the dot. The following second kind of Lagrange’s equations are applied to derive the main EOM for the examined system.

$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\eta }}\right)-\left(\frac{\partial L}{\partial \eta }\right)=Q_{\eta }; \eta =\left(X,\theta \right).$

(3)

Here, $L=T-V$ expresses the Lagrangian and $Q_{\eta }$ represents the related forces of the coordinates $\eta =\left(X,\theta \right)$ which have the forms

$Q_X=F_1\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{1}t-C_1\dot{X}-\gamma \upsilon ,Q_{\theta }=M_0\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{3}t+\left(l_0+X\right)F_2\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\Omega }_{2}t-C_2\dot{\theta },$

(4)

where $\upsilon$ is the voltage of the load’s resistance $R_p$ . Additionally, the mechanism’s equation for the circuit of piezoelectric may has the form ²⁰

$c_p\dot{\upsilon }+\frac{\upsilon }{R_p}-\gamma \dot{X}=0.$

(5)

Now, let’s have a look at the following dimensionless parameters

$\begin{array}{l}{\omega }_1^2=\frac{k_1}{m}, {\omega }_2^2=\frac{g}{l}, {\omega }^2=\frac{{\omega }_2^2}{{\omega }_1^2}, r_X=\frac{R_X}{l}, r_Y=\frac{R_Y}{l}, p_X=\frac{{\Omega }_X}{{\omega }_1}, p_Y=\frac{{\Omega }_Y}{{\omega }_1}, p_1=\frac{{\Omega }_1}{{\omega }_1}, p_2=\frac{{\Omega }_2}{{\omega }_1},\\ p_3=\frac{{\Omega }_3}{{\omega }_1}, \alpha =\frac{k_2{l}^2}{{\omega }_1^{2}m}, \mu =\frac{\gamma }{{\omega }_1^{2}ml}, c_1=\frac{C_1}{{\omega }_{1}m}, c_2=\frac{C_2}{l^2{\omega }_{1}m}, f_1=\frac{F_1}{l{\omega }_1^{2}m}, f_2=\frac{F_2}{l^2{\omega }_1^{2}m}, f_3=\frac{M_0}{{\omega }_1^{2}m{l}^2},\\ Z=\frac{X-Z_r}{l}, l=l_0+Z_r, {\zeta }_s=\frac{Z_r}{l}, k_1{Z}_r+k_2{Z}_r^3=mg, {\zeta }_s+\alpha {\zeta }_s^3={\omega }^2, \tau ={\omega }_{1}t,\end{array}$

(6)

where $Z_r$ is the spring’s static elongation. The substitution of (2), (4), and (6) into (3) produces the below governing EOM in dimensionless forms

$\ddot{Z}+c_1\dot{Z}+Z+\alpha Z^3+3\alpha {\zeta }_{s}Z\left({\zeta }_s+Z\right)+r_X{p}_X^2\left(\frac{1}{2}{\theta }^2-1\right)\mathit{cos}\left(p_X\tau \right)+r_Y{p}_Y^2\theta \times \left(\frac{1}{6}{\theta }^2-1\right)\mathit{sin}\left(p_Y\tau \right)+\frac{1}{2}{\omega }^2{\theta }^2+\mu \upsilon -\left(1+Z\right){\dot{\theta }}^2=f_1\hspace{0.17em}\mathit{cos}\left(p_1\tau \right),$

(7)

${\left(1+Z\right)}^2\ddot{\theta }+c_2\dot{\theta }+\left(1+Z\right){\omega }^2\left(\theta -\frac{{\theta }^3}{6}\right)+2\left(1+Z\right)\dot{Z}\dot{\theta }+\left(1+Z\right)\left(1-\frac{1}{6}{\theta }^2\right)\theta r_X{p}_X^2\hspace{0.17em}\mathit{cos}\left(p_X\tau \right)+\left(1+Z\right)\left(\frac{1}{2}{\theta }^2-1\right)r_Y{p}_Y^2\hspace{0.17em}\mathit{sin}\left(p_Y\tau \right)=f_3\hspace{0.17em}\mathit{cos}\left(p_3\tau \right)+l\left(1+Z\right)f_2\hspace{0.17em}\mathit{cos}\left(p_2\tau \right).$

(8)

The dimensionless mechanics equation of the piezoelectric circuit may be determined using the procedure above by substituting (6) into (5) as shown below.

$\dot{\upsilon }+\frac{\upsilon }{R_p{c}_p{\omega }_1}=\frac{l\gamma }{c_p}\dot{Z}.$

(9)

The suggested approach

In the current section, the MST is used to obtain the asymptotic analytic solutions of the aforementioned system of equations (7)–(9). Therefore, we focus our study on the dynamical behavior of this system in a tiny region that is closed to the point of static equilibrium point. ²³ Then, one can express its vibrational amplitudes as follows

$Z\left(\tau \right)=\epsilon \tilde{\chi }\left(\tau ;\epsilon \right), \theta \left(\tau \right)=\epsilon \tilde{\vartheta }\left(\tau ;\epsilon \right), \upsilon \left(\tau \right)=\epsilon \tilde{v}\left(\tau ;\epsilon \right),$

(10)

where $0<\epsilon \ll 1$ is a small parameter. The smallness of parameters and variables was taken into consideration in accordance with

$$\begin{gathered} r_X={\epsilon }^2{\tilde{r}}_X, r_Y={\epsilon }^2{\tilde{r}}_Y, \mu ={\epsilon }^2\tilde{\mu }, c_1={\epsilon }^2{\tilde{c}}_1, c_2={\epsilon }^2{\tilde{c}}_2, \alpha ={\epsilon }^2\tilde{\alpha }, \\ {f}_1={\epsilon }^3{\tilde{f}}_1, f_2={\epsilon }^3{\tilde{f}}_2, f_3={\epsilon }^3{\tilde{f}}_3, \gamma ={\epsilon }^2\tilde{\gamma }, R_p=\frac{{\tilde{R}}_p}{{\epsilon }^2}, c_p={\epsilon }^2{\tilde{c}}_p. \end{gathered}$$

(11)

According to the MST methodology, the solutions $\stackrel{\sim }{\chi },\stackrel{\sim }{\vartheta },$ and $\stackrel{\sim }{v}$ , can be expressed as as shown below ²⁶

$\begin{array}{ll}\stackrel{\sim }{\chi }={\displaystyle \sum }_{k=1}^3{\epsilon }^k{\chi }_k\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^4\right),& \\ \stackrel{\sim }{v}={\displaystyle \sum }_{k=1}^3{\epsilon }^k{v}_k\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^4\right).& \\ \stackrel{\sim }{\vartheta }={\displaystyle \sum }_{k=1}^3{\epsilon }^k{\vartheta }_k\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^4\right),& \end{array}$

(12)

Here ${\tau }_n={\epsilon }^n\tau \left(n=0,1,2\right)$ represent specific time scales, in which ${\tau }_1$ is a fast one while ${\tau }_2$ and ${\tau }_3$ are the slow ones. Now, the derivatives with regard to $\tau$ must be transformed to these scales. Accordingly, the use of the following operators fulfils this purpose

$\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).$

(13)

It is clear from these operators that terms $O\left({\epsilon }^3\right)$ and higher are not taken into account for their smallness. The substitution of (10)–(13) into equations (7)–(9) yields the below sets of partial differential equations that are linked to the various powers of $\epsilon$ :

(a) Order $\epsilon$

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

(14)

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

(15)

$\frac{\partial v_1}{\partial {\tau }_0}+\frac{v_1}{{\stackrel{\sim }{c}}_p{\stackrel{\sim }{R}}_p{\omega }_1}=\frac{l\stackrel{\sim }{\gamma }}{{\stackrel{\sim }{c}}_p}\frac{\partial {\chi }_1}{\partial {\tau }_0}.$

(16)

(b) Order of ${\epsilon }^2$

$\frac{{\partial }^2{\chi }_2}{\partial {\tau }_0^2}+{\chi }_2=\frac{1}{2}{p}_X^2{\stackrel{\sim }{r}}_X{e}^{i{p}_X{\tau }_0}-\frac{1}{2}{\omega }^2{\vartheta }_1^2+{\left(\frac{\partial {\vartheta }_1}{\partial {\tau }_0}\right)}^2-2\frac{{\partial }^2{\chi }_1}{\partial {\tau }_0\partial {\tau }_1},$

(17)

$\frac{{\partial }^2{\vartheta }_2}{\partial {\tau }_0^2}+{\omega }^2{\vartheta }_2=\frac{1}{2}i{\stackrel{\sim }{r}}_Y{p}_Y^2{e}^{i{p}_Y{\tau }_0}-2\left(\frac{{\partial }^2{\vartheta }_1}{\partial {\tau }_0{\tau }_1}+\frac{\partial {\chi }_1}{\partial {\tau }_0}\frac{\partial {\vartheta }_1}{\partial {\tau }_0}\right)-{\chi }_1\left({\omega }^2{\vartheta }_1+2\frac{{\partial }^2{\vartheta }_1}{\partial {\tau }_0^2}\right),$

(18)

$\frac{\partial v_2}{\partial {\tau }_0}+\frac{v_2}{{\stackrel{\sim }{c}}_p{\stackrel{\sim }{R}}_p{\omega }_1}=-\frac{\partial v_1}{\partial {\tau }_1}+\frac{l\stackrel{\sim }{\gamma }}{{\stackrel{\sim }{c}}_p}\left(\frac{\partial {\chi }_1}{\partial {\tau }_1}+\frac{\partial {\chi }_2}{\partial {\tau }_0}\right).$

(19)

(c) Order of ${\epsilon }^3$

$\frac{{\partial }^2{\chi }_3}{\partial {\tau }_0^2}+{\chi }_3=\frac{1}{2}{\stackrel{\sim }{f}}_1{e}^{i{p}_1{\tau }_0}-\frac{1}{2}i{\stackrel{\sim }{r}}_Y{p}_Y^2{\vartheta }_1{e}^{i{p}_Y{\tau }_0}-{\omega }^2{\vartheta }_1{\vartheta }_2-{\stackrel{\sim }{c}}_1\frac{\partial {\chi }_1}{\partial {\tau }_0}-\frac{{\partial }^2{\chi }_1}{\partial {\tau }_1^2}-3\stackrel{\sim }{\alpha }{{\zeta }_s}^2{\chi }_1+2\left(\frac{\partial {\vartheta }_1}{\partial {\tau }_0}\frac{\partial {\vartheta }_1}{\partial {\tau }_1}+\frac{\partial {\vartheta }_1}{\partial {\tau }_0}\frac{\partial {\vartheta }_2}{\partial {\tau }_0}\right)-2\left(\frac{{\partial }^2{\chi }_1}{\partial {\tau }_0\partial {\tau }_2}+\frac{{\partial }^2{\chi }_2}{\partial {\tau }_0\partial {\tau }_1}\right)+{\phi }_1{\left(\frac{\partial {\vartheta }_1}{\partial {\tau }_0}\right)}^2-\stackrel{\sim }{\mu }{v}_1,$

(20)

$\frac{{\partial }^2{\vartheta }_3}{\partial {\tau }_0^2}+{\omega }^2{\vartheta }_3=\frac{1}{2}\left(l{\stackrel{\sim }{f}}_2{e}^{i{p}_2{\tau }_0}+{\stackrel{\sim }{f}}_3{e}^{i{p}_3{\tau }_0}-i{\stackrel{\sim }{r}}_Y{p}_Y^2{\chi }_1{e}^{i{p}_Y{\tau }_0}-{\stackrel{\sim }{r}}_X{p}_X^2{\vartheta }_1{e}^{i{p}_X{\tau }_0}\right)-{\stackrel{\sim }{c}}_2\frac{\partial {\vartheta }_1}{\partial {\tau }_0}-{\omega }^2\left({\vartheta }_2{\chi }_1+{\vartheta }_1{\chi }_2\right)+\frac{1}{6}{\omega }^2{\vartheta }_1^3-2{\chi }_1\left(\frac{\partial {\chi }_1}{\partial {\tau }_0}\frac{\partial {\vartheta }_1}{\partial {\tau }_0}+2\frac{{\partial }^2{\vartheta }_1}{\partial {\tau }_0\partial {\tau }_1}+\frac{{\partial }^2{\vartheta }_2}{\partial {\tau }_0^2}\right)-2\frac{\partial {\chi }_1}{\partial {\tau }_0}\left(\frac{\partial {\vartheta }_1}{\partial {\tau }_1}+\frac{\partial {\vartheta }_2}{\partial {\tau }_0}\right)-2\frac{\partial {\vartheta }_1}{\partial {\tau }_0}\left(\frac{\partial {\chi }_1}{\partial {\tau }_1}+\frac{\partial {\chi }_2}{\partial {\tau }_0}\right)-\frac{{\partial }^2{\vartheta }_1}{\partial {\tau }_0^2}\left({\chi }_1^2+2{\chi }_2\right)-2\left(\frac{{\partial }^2{\vartheta }_1}{\partial {\tau }_0\partial {\tau }_2}+\frac{{\partial }^2{\vartheta }_2}{\partial {\tau }_0\partial {\tau }_1}\right)-\frac{{\partial }^2{\vartheta }_1}{\partial {\tau }_1^2},$

(21)

$\frac{\partial v_3}{\partial {\tau }_0}+\frac{v_3}{{\stackrel{\sim }{c}}_p{\stackrel{\sim }{R}}_p{\omega }_1}=-\frac{\partial v_2}{\partial {\tau }_1}-\frac{\partial v_1}{\partial {\tau }_2}+\frac{l\stackrel{\sim }{\gamma }}{{\stackrel{\sim }{c}}_p}\left(\frac{\partial {\chi }_3}{\partial {\tau }_0}+\frac{\partial {\chi }_2}{\partial {\tau }_1}+\frac{\partial {\chi }_1}{\partial {\tau }_2}\right).$

(22)

Within the framework of the previous equations (14)–(22), it can be stated that they can be solved sequentially which emphasizes the significance of the solutions that are provided by the first group of equations (14)–(16). Consequently, the general solutions of equations (13)–(15) have the forms

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

(23)

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

(24)

$v_1={\stackrel{\sim }{R}}_p\stackrel{\sim }{\gamma }l{\omega }_1\left(\frac{A_1{e}^{i{\tau }_0}}{{\stackrel{\sim }{R}}_p{\stackrel{\sim }{c}}_p{\omega }_1-i}+\frac{{\overline{A}}_1{e}^{-i{\tau }_0}}{{\stackrel{\sim }{R}}_p{\stackrel{\sim }{c}}_p{\omega }_1+i}\right),$

(25)

where $A_s\left(s=1,2\right)$ portray unclear, complex processes functioning at slow time scales ${\tau }_s$ , while ${\overline{A}}_s$ indicate their complex conjugate. It is known that the secular terms can be produced by substituting the previous solutions (23)–(25) into the higher order partial differential equations (17)–(19). The following conditions can be used to eliminate them

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

(26)

Consequently, the subsequent second-order solutions are as follows

${\chi }_2={\omega }^2{A}_2{\overline{A}}_2+\frac{{\stackrel{\sim }{r}}_X{p}_X^2}{2\left(1-p_X^2\right)}{e}^{i{p}_X{\tau }_0}+\frac{3{\omega }^2{A}_2^2}{2\left(4{\omega }^2-1\right)}{e}^{2i\omega {\tau }_0}+CC,$

(27)

${\vartheta }_2=\frac{\omega \left(\omega -2\right){\overline{A}}_1{A}_2}{\left(2\omega -1\right)}{e}^{i\left(\omega -1\right){\tau }_0}-\frac{\omega \left(\omega +2\right)A_1{A}_2}{\left(2\omega +1\right)}{e}^{i\left(\omega +1\right){\tau }_0}+\frac{i{\stackrel{\sim }{r}}_Y{p}_Y^2}{2\left(p_Y^2-{\omega }^2\right)}{e}^{i{p}_Y{\tau }_0}+CC,$

(28)

$v_2=\frac{i{\omega }_{1}l\stackrel{\sim }{\gamma }{\stackrel{\sim }{r}}_X{p}_X^3{\stackrel{\sim }{R}}_p}{2\left(p_X^2-1\right)\left({\omega }_1{\stackrel{\sim }{c}}_p{p}_X^2{\stackrel{\sim }{R}}_p-1\right)}{e}^{i{p}_X{\tau }_0}-\frac{3i{\omega }^3{A}_2^2{\stackrel{\sim }{R}}_p{\omega }_{1}l\stackrel{\sim }{\gamma }}{\left(4{\omega }^2-1\right)\left(4{\stackrel{\sim }{R}}_p{\omega }^2{\omega }_1{\stackrel{\sim }{c}}_p-1\right)}{e}^{2i\omega {\tau }_0}+CC,$

(29)

where $CC$ indicates the conjugations of the previous expressions. Similar to the earlier case, the canceling of the secular terms requested that $A_s$ depended on ${\tau }_2$ . Therefore, canceling these terms produces the following conditions

$i{\stackrel{\sim }{c}}_1{A}_1+2i\frac{\partial A_1}{\partial {\tau }_2}+3\stackrel{\sim }{\alpha }{\zeta }_s^2{A}_1+\frac{{\stackrel{\sim }{R}}_p{\omega }_{1}l\stackrel{\sim }{\gamma }\stackrel{\sim }{\mu }{A}_1}{\left({\omega }_1{\stackrel{\sim }{c}}_p{\stackrel{\sim }{R}}_p-i\right)}+\frac{6{\omega }^2\left({\omega }^2-1\right)A_1{A}_2{\overline{A}}_2}{\left(4{\omega }^2-1\right)}=0,$

(30)

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

(31)

Finally, we can get the below third-order solutions

$x_3=\frac{{\stackrel{\sim }{f}}_1{e}^{i{p}_1{\tau }_0}}{2\left(1-p_1^2\right)}+3A_1\left\{\frac{{\left(\omega +{\omega }^2\right)}^2{A}_2^2{e}^{i\left(1+2\omega \right){\tau }_0}}{\left(1+2\omega \right)\left[1-{\left(1+2\omega \right)}^2\right]}+\frac{{\left({\omega }^2-\omega \right)}^2{\overline{A}}_2^2{e}^{i\left(1-2\omega \right){\tau }_0}}{\left(1-2\omega \right)\left[1-{\left(1-2\omega \right)}^2\right]}\right\}-\frac{i{\stackrel{\sim }{r}}_Y{p}_Y^3}{2}\left\{\frac{\left(p_Y-2\omega \right){\overline{A}}_2{e}^{i\left(p_Y-\omega \right){\tau }_0}}{\left(p_Y^2-{\omega }^2\right)\left[1-{\left(p_Y-\omega \right)}^2\right]}+\frac{\left(p_Y+2\omega \right)A_2{e}^{i\left(p_Y+\omega \right){\tau }_0}}{\left(p_Y^2-{\omega }^2\right)\left[1-{\left(p_Y+\omega \right)}^2\right]}\right\}+CC,$

(32)

${\mathit{\phi }}_{\mathbf{3}}=\frac{{\stackrel{\sim }{\mathit{f}}}_{\mathbf{2}}{\mathit{e}}^{\mathit{i}{\mathit{p}}_{\mathbf{2}}{\mathit{\tau }}_{\mathbf{0}}}}{\mathbf{2}\left({\mathit{\omega }}^{\mathbf{2}}-{\mathit{p}}_{\mathbf{2}}^{\mathbf{2}}\right)}+\frac{{\stackrel{\sim }{\mathit{f}}}_{\mathbf{3}}{\mathit{e}}^{\mathit{i}{\mathit{p}}_{\mathbf{3}}{\mathit{\tau }}_{\mathbf{0}}}}{\mathbf{2}\left({\mathit{\omega }}^{\mathbf{2}}-{\mathit{p}}_{\mathbf{3}}^{\mathbf{2}}\right)}-\frac{\left(\mathbf{49}{\mathit{\omega }}^{\mathbf{2}}-\mathbf{1}\right){\mathit{A}}_{\mathbf{2}}^{\mathbf{3}}{\mathit{e}}^{\mathbf{3}\mathit{i}\mathit{\omega }{\mathit{\tau }}_{\mathbf{0}}}}{\mathbf{48}\left(\mathbf{4}{\mathit{\omega }}^{\mathbf{2}}-\mathbf{1}\right)}+\frac{{\mathit{p}}_{\mathit{X}}{\stackrel{\sim }{\mathit{r}}}_{\mathit{X}}}{\mathbf{2}\left({\mathit{p}}_{\mathit{X}}^{\mathbf{2}}-\mathbf{1}\right)}\left\{\frac{\left({\mathit{p}}_{\mathit{X}}-\mathbf{1}+\mathit{\omega }\right)}{\left({\mathit{p}}_{\mathit{X}}+\mathbf{2}\mathit{\omega }\right)}\times \left({\mathit{p}}_{\mathit{X}}+\mathbf{1}+\mathit{\omega }\right){\mathit{A}}_{\mathbf{2}}{\mathit{e}}^{\mathit{i}\left({\mathit{p}}_{\mathit{X}}+\mathit{\omega }\right){\mathit{\tau }}_{\mathbf{0}}}+\frac{\left({\mathit{p}}_{\mathit{X}}-\mathbf{1}-\mathit{\omega }\right)\left({\mathit{p}}_{\mathit{X}}+\mathbf{1}-\mathit{\omega }\right)}{\left({\mathit{p}}_{\mathit{X}}-\mathbf{2}\mathit{\omega }\right)}{\overline{\mathit{A}}}_{\mathbf{2}}{\mathit{e}}^{\mathit{i}\left({\mathit{p}}_{\mathit{X}}-\mathit{\omega }\right){\mathit{\tau }}_{\mathbf{0}}}\right\}-\frac{\mathit{i}{\mathit{p}}_{\mathit{Y}}^{\mathbf{3}}{\stackrel{\sim }{\mathit{r}}}_{\mathit{Y}}}{\mathbf{2}\left({\mathit{p}}_{\mathit{Y}}^{\mathbf{2}}-{\mathit{\omega }}^{\mathbf{2}}\right)}\times \left\{\frac{\left({\mathit{p}}_{\mathit{Y}}+\mathbf{2}\right){\mathit{A}}_{\mathbf{1}}{\mathit{e}}^{\mathit{i}\left({\mathit{p}}_{\mathit{Y}}+\mathbf{1}\right){\mathit{\tau }}_{\mathbf{0}}}}{\left[{\mathit{\omega }}^{\mathbf{2}}-{\left({\mathit{p}}_{\mathit{Y}}+\mathbf{1}\right)}^{\mathbf{2}}\right]}+\frac{\left({\mathit{p}}_{\mathit{Y}}-\mathbf{2}\right){\overline{\mathit{A}}}_{\mathbf{1}}{\mathit{e}}^{\mathit{i}\left({\mathit{p}}_{\mathit{Y}}-\mathbf{1}\right){\mathit{\tau }}_{\mathbf{0}}}}{\left[{\mathit{\omega }}^{\mathbf{2}}-{\left({\mathit{p}}_{\mathit{Y}}-\mathbf{1}\right)}^{\mathbf{2}}\right]}\right\}+\frac{\mathit{\omega }{\mathit{A}}_{\mathbf{2}}}{\mathbf{4}}\left\{\frac{\left({\mathit{\omega }}^{\mathbf{2}}+\mathbf{5}\mathit{\omega }+\mathbf{6}\right){\mathit{A}}_{\mathbf{1}}^{\mathbf{2}}{\mathit{e}}^{\mathit{i}\left(\mathit{\omega }+\mathbf{2}\right){\mathit{\tau }}_{\mathbf{0}}}}{\left(\mathbf{2}\mathit{\omega }+\mathbf{1}\right)}+\frac{\left({\mathit{\omega }}^{\mathbf{2}}-\mathbf{5}\mathit{\omega }+\mathbf{6}\right){\overline{\mathit{A}}}_{\mathbf{1}}^{\mathbf{2}}{\mathit{e}}^{\mathit{i}\left(\mathit{\omega }-\mathbf{2}\right){\mathit{\tau }}_{\mathbf{0}}}}{\left(\mathbf{2}\mathit{\omega }-\mathbf{1}\right)}\right\}+\mathit{C}\mathit{C},$

(33)

${\mathit{v}}_{\mathbf{3}}=-\mathit{i}{\mathit{\omega }}_{\mathbf{1}}{\stackrel{\sim }{\mathit{R}}}_{\mathit{p}}\mathit{l}\stackrel{\sim }{\mathit{\gamma }}{\stackrel{\sim }{\mathit{r}}}_{\mathit{Y}}{\mathit{p}}_{\mathit{Y}}^{\mathbf{3}}\left\{\frac{\left({\mathit{p}}_{\mathit{Y}}-\mathbf{2}\mathit{\omega }\right){\overline{\mathit{A}}}_{\mathbf{2}}{\mathit{e}}^{\mathit{i}{\mathit{\tau }}_{\mathbf{0}}\left({\mathit{p}}_{\mathit{Y}}-\mathit{\omega }\right)}}{\mathbf{2}\left[\mathbf{1}+\mathit{i}{\mathit{\omega }}_{\mathbf{1}}{\stackrel{\sim }{\mathit{c}}}_{\mathit{p}}{\stackrel{\sim }{\mathit{R}}}_{\mathit{p}}\left({\mathit{p}}_{\mathit{Y}}-\mathit{\omega }\right)\right]\left({\mathit{p}}_{\mathit{Y}}+\mathit{\omega }\right)\left({\mathit{p}}_{\mathit{Y}}-\mathit{\omega }-\mathbf{1}\right)\left({\mathit{p}}_{\mathit{Y}}-\mathit{\omega }+\mathbf{1}\right)}+\frac{\left({\mathit{p}}_{\mathit{Y}}+\mathbf{2}\mathit{\omega }\right){\mathit{p}}_{\mathit{Y}}{\mathit{A}}_{\mathbf{2}}{\mathit{e}}^{\mathit{i}{\mathit{\tau }}_{\mathbf{0}}\left({\mathit{p}}_{\mathit{Y}}+\mathit{\omega }\right)}}{\mathbf{2}\left[\mathbf{1}+\mathit{i}{\mathit{\omega }}_{\mathbf{1}}{\stackrel{\sim }{\mathit{c}}}_{\mathit{p}}{\stackrel{\sim }{\mathit{R}}}_{\mathit{p}}\left({\mathit{p}}_{\mathit{Y}}+\mathit{\omega }\right)\right]\left({\mathit{p}}_{\mathit{Y}}-\mathit{\omega }\right)\left({\mathit{p}}_{\mathit{Y}}+\mathit{\omega }-\mathbf{1}\right)\left({\mathit{p}}_{\mathit{Y}}+\mathit{\omega }+\mathbf{1}\right)}\right\}-\frac{\mathbf{3}}{\mathbf{4}}\mathit{i}{\mathit{\omega }}_{\mathbf{1}}{\stackrel{\sim }{\mathit{R}}}_{\mathit{p}}\mathit{l}\stackrel{\sim }{\mathit{\gamma }}\mathit{\omega }{\mathit{A}}_{\mathbf{1}}\left\{\frac{\left(\mathit{\omega }+\mathbf{1}\right){\mathit{A}}_{\mathbf{2}}^{\mathbf{2}}{\mathit{e}}^{\mathit{i}{\mathit{\tau }}_{\mathbf{0}}\left(\mathbf{2}\mathit{\omega }+\mathbf{1}\right)}}{\left[\mathbf{1}+\mathit{i}{\mathit{\omega }}_{\mathbf{1}}{\stackrel{\sim }{\mathit{c}}}_{\mathit{p}}{\stackrel{\sim }{\mathit{R}}}_{\mathit{p}}\left(\mathbf{2}\mathit{\omega }+\mathbf{1}\right)\right]}+\frac{\left(\mathit{\omega }-\mathbf{1}\right){\overline{\mathit{A}}}_{\mathbf{2}}^{\mathbf{2}}{\mathit{e}}^{\mathit{i}{\mathit{\tau }}_{\mathbf{0}}\left(-\mathbf{2}\mathit{\omega }+\mathbf{1}\right)}}{\left[\mathbf{1}+\mathit{i}{\mathit{\omega }}_{\mathbf{1}}{\stackrel{\sim }{\mathit{c}}}_{\mathit{p}}{\stackrel{\sim }{\mathit{R}}}_{\mathit{p}}\left(\mathbf{1}-\mathbf{2}\mathit{\omega }\right)\right]}\right\}+\frac{\mathit{i}{\mathit{\omega }}_{\mathbf{1}}{\stackrel{\sim }{\mathit{R}}}_{\mathit{p}}\mathit{l}\stackrel{\sim }{\mathit{\gamma }}{\mathit{p}}_{\mathbf{1}}{\stackrel{\sim }{\mathit{f}}}_{\mathbf{1}}{\mathit{e}}^{\mathit{i}{\mathit{p}}_{\mathbf{1}}{\mathit{\tau }}_{\mathbf{0}}}}{\mathbf{2}\left(\mathbf{1}+\mathit{i}{\mathit{\omega }}_{\mathbf{1}}{\mathit{p}}_{\mathbf{1}}{\stackrel{\sim }{\mathit{c}}}_{\mathit{p}}{\stackrel{\sim }{\mathit{R}}}_{\mathit{p}}\right)\left(\mathbf{1}-{\mathit{p}}_{\mathbf{1}}^{\mathbf{2}}\right)}+\mathit{C}\mathit{C}.$

(34)

Resonance categorization and modulation equations

This section’s goals are to categorize the resonance cases, address with one potential case, and obtain the modulation equations. When the denominators in the second and third order solutions disappear, these situations may appears. ⁴¹ Consequently, the following possibilities are obtained

1- The case of primary external resonance is obtained at $p_1=1,p_2=\omega ,p_3=\omega ,$

2- The case of internal resonance will be satisfied at $\omega =\pm 0.5,$ $p_Y=\omega ,$ $p_Y=\omega +1,$ $p_Y=1-\omega ,\omega =0,$ $p_X=\pm 1,p_X=\pm 2\omega .$

When one of the preceding prerequisites for resonance is satisfied, the examined system will behave in a very complicated manner. It’s crucial to note that even if the oscillations are far from resonance, the solutions produced are still appropriate. We’ll now look into the first case in greater depth, in which two of them occurred simultaneously, that is, $p_1\approx 1,p_2\approx \omega$ . Then, one can consider the detuning parameters ${\sigma }_s\left(s=1,2\right)$ that help achieve this aim, which express the nearness of $p_1$ and $p_2$ to $1$ and $\omega$ , respectively. Therefore, one writes

$p_1=1+{\sigma }_1,p_2=\omega +{\sigma }_2,$

(35)

where ${\sigma }_s={\epsilon }^2{\stackrel{\sim }{\sigma }}_s$ . Detuning parameters are seen as a distance from the vibrations and the strict resonance. Eliminating secular terms makes it possible to acquire the criteria of solvability. As a consequence, the next equations estimate that the necessary conditions are fulfilled

$$\begin{gathered} \frac{1}{2}{\tilde{f}}_1{e}^{i{\tau }_1{\tilde{\sigma }}_1}-i{\tilde{c}}_1{A}_1-2i\frac{\partial A_1}{\partial {\tau }_2}-\frac{{\tilde{R}}_p{\omega }_{1}l\tilde{\gamma }\tilde{\mu }{A}_1}{\left({\omega }_1{\tilde{c}}_p{\tilde{R}}_p-i\right)}-\frac{6{\omega }^2\left({\omega }^2-1\right)A_1{A}_2{\overline{A}}_2}{\left(4{\omega }^2-1\right)}-3\tilde{\alpha }{\zeta }_s^2{A}_1=0, \\ \frac{1}{2}{\tilde{f}}_2{e}^{i{\tau }_1{\tilde{\sigma }}_2}-i\omega {\tilde{c}}_2{A}_2-2i\omega \frac{\partial A_2}{\partial {\tau }_2}-\frac{{\omega }^2\left({\omega }^2-1\right)A_2}{\left(4{\omega }^2-1\right)}\left[6A_1{\overline{A}}_1-\frac{1}{2}\left(8{\omega }^2+1\right)A_2{\overline{A}}_2\right]=0. \end{gathered}$$

(36)

The established requirements of solvability (26) and (36) constitute four nonlinear partial differential equations regarding to $A_s$ that rely on the slow time scale ${\tau }_2$ . As a result, we may consider the polar form of these functions as:

${\mathit{A}}_{\mathit{s}}\left({\mathit{\tau }}_{\mathbf{2}}\right)=\frac{{\tilde{\mathit{h}}}_{\mathit{s}}\left({\mathit{\tau }}_{\mathbf{2}}\right)}{\mathbf{2}}{\mathit{e}}^{\mathit{i}{\mathit{\psi }}_{\mathit{s}}{\mathit{\tau }}_{\mathbf{2}}}; {\mathit{h}}_{\mathit{s}}=\mathit{\epsilon }{\tilde{\mathit{h}}}_{\mathit{s}}\left(\mathit{s}=\mathbf{1},\mathbf{2}\right).$

(37)

Since the functions $A_s$ don’t depend on the scales ${\tau }_0$ and ${\tau }_1$ , it is possible to simplify the first-order derivative operator as follows

$\frac{\partial A_s}{\partial {\tau }_0}=\frac{\partial A_s}{\partial {\tau }_1}=0,\frac{\partial A_s}{\partial \tau }={\epsilon }^2\frac{\partial A_s}{\partial {\tau }_2}.$

(38)

Referring to (38), the criteria of solvability (36) may be transformed into ordinary differential equations by using next modified phases

${\mathit{\theta }}_{\mathit{s}}\left({\mathit{\tau }}_{\mathbf{1}},{\mathit{\tau }}_{\mathbf{2}}\right)={\mathit{\tau }}_{\mathbf{1}}{\tilde{\mathit{\sigma }}}_{\mathit{s}}-{\mathit{\psi }}_{\mathit{s}}\left({\mathit{\tau }}_{\mathbf{2}}\right), \left(\mathit{s}=\mathbf{1},\mathbf{2}\right).$

(39)

Substituting (37)–(39) into (36) and separating the real and imaginary portions, one obtains the below required system

$\begin{array}{ll}{h}_1\frac{d{\theta }_1}{d\tau }=h_1\left({\sigma }_1-S_1\right)-\frac{3}{4}\frac{{\omega }^2\left({\omega }^2-1\right)}{\left(4{\omega }^2-1\right)}{h}_1{h}_2^2-\frac{3}{2}\alpha {\zeta }_s^2{h}_1+\frac{f_1}{2}\mathit{cos}\hspace{0.17em}{\theta }_1,& \\ \frac{d{h}_1}{d\tau }=-h_1\left(\frac{c_1}{2}+S_2\right)+\frac{f_1}{2}\mathit{sin}\hspace{0.17em}{\theta }_1,& \frac{d{h}_2}{d\tau }=-\frac{1}{2}{c}_2{h}_2+\frac{f_2}{2\omega }\mathit{sin}\hspace{0.17em}{\theta }_2.\\ h_2\frac{d{\theta }_2}{d\tau }=h_2{\sigma }_2+\frac{\omega \left({\omega }^2-1\right)}{4\left(4{\omega }^2-1\right)}\left(\frac{\left(8{\omega }^2+1\right)}{4}{h}_2^3-3h_2{h}_1^2\right)+\frac{f_2}{2\omega }\mathit{cos}\hspace{0.17em}{\theta }_2,& \end{array}$

(40)

Here, ${\mathit{S}}_{\mathbf{1}}=\frac{{\mathit{R}}_{\mathit{p}}^{\mathbf{2}}{\mathit{\omega }}_{\mathbf{1}}^{\mathbf{2}}\mathit{\mu }\mathit{l}\mathit{\gamma }{\mathit{c}}_{\mathit{p}}}{\mathbf{2}\left(\mathbf{1}+{\mathit{\omega }}_{\mathbf{1}}^{\mathbf{2}}{\mathit{c}}_{\mathit{p}}^{\mathbf{2}}{\mathit{R}}_{\mathit{p}}^{\mathbf{2}}\right)}$ and ${\mathit{S}}_{\mathbf{2}}=\frac{{\mathit{R}}_{\mathit{p}}{\mathit{\omega }}_{\mathbf{1}}\mathit{\mu }\mathit{l}\mathit{\gamma }}{\mathbf{2}\left(\mathbf{1}+{\mathit{\omega }}_{\mathbf{1}}^{\mathbf{2}}{\mathit{c}}_{\mathit{p}}^{\mathbf{2}}{\mathit{R}}_{\mathit{p}}^{\mathbf{2}}\right)}$ .

Numerical results

This section is devoted to examining the numerical solutions of the aforementioned system (40) using the Runge–Kutta method from fourth-order with the help of Wolfram Mathematica 13.2. It must be mentioned that these equations cannot be solved analytically to obtain the amplitudes and the modified phases $h_s\left(s=1,2\right)$ and ${\theta }_s$ . Therefore, the following listed initial circumstances are taken into account $h_1\left(0\right)=0.01,$ $h_2\left(0\right)=0.02,{\theta }_1\left(0\right)={\theta }_2\left(0\right)=0.$ In Figures 2 and 3, $h_s$ and ${\theta }_s$ are displayed to show the temporal histories of the amplitudes $h_s\left(s=1,2\right)$ and the adjusted phases ${\theta }_s$ . It has taken into account the below values of the used parameters for the sketched curves:

$\mathit{\alpha }=\mathbf{0.140625},\mathit{\mu }=\mathbf{0.06667},{\mathit{\zeta }}_{\mathit{s}}=\mathbf{0.032667},{\mathit{f}}_{\mathbf{1}}=\mathbf{0.0006},{\mathit{f}}_{\mathbf{2}}=\mathbf{0.0006222},{\mathit{f}}_{\mathbf{3}}=\mathbf{0.004444},{\mathit{c}}_{\mathbf{1}}=\mathbf{0.005},{\mathit{c}}_{\mathbf{2}}=\mathbf{0.012},{\mathit{p}}_{\mathit{X}}=\mathbf{0.01},{\mathit{p}}_{\mathit{Y}}=\mathbf{0.0125},{\mathit{r}}_{\mathit{X}}=\mathbf{0.05333},{\mathit{r}}_{\mathit{Y}}=\mathbf{0.08},{\mathit{p}}_{\mathbf{1}}={\mathit{p}}_{\mathbf{2}}={\mathit{p}}_{\mathbf{3}}=\mathbf{0.005},{\mathit{\sigma }}_{\mathbf{1}}=\mathbf{0.002},{\mathit{\sigma }}_{\mathbf{2}}=\mathbf{0.002}.$

Figure 2.

The amplitude $h_1$ and the modified phase ${\theta }_1$ as a function of the dimensionless time parameter $\tau$ at: (a), (b) ${\omega }_1\left(=15.5,17.9,20\right)$ , and (c), (d) ${\omega }_2\left(=2.9,3.2,3.6\right)$ .

Figure 3.

The amplitude $h_2$ and the modified phase ${\theta }_2$ as a function of the dimensionless time parameter $\tau$ at: (a), (b) ${\omega }_1\left(=15.5,17.9,20\right)$ and (c), (d) ${\omega }_2\left(=2.9,3.2,3.6\right)$ .

The calculations for these figures were performed for various values of ${\omega }_1\left(=15.5,17.9,20\right)$ and ${\omega }_2\left(=2.9,3.2,3.6\right)$ as indicated, respectively, in portions [(a), (b)] and [(c), (d)]. As time goes on, the oscillations of the waves described by the amplitude $h_1$ have decaying behaviors, as shown in Figure 2(a) and (c). On the other side, the behaviors of the $h_2$ waves increase to certain values, and then they behave in a stationary manner, as seen in Figure 3(a) and (c). Furthermore, the variance of the ${\theta }_s$ typically displays a stable behavior once it reaches a certain value over time, as seen in Figures 2(b) and (d) and 3 (b) and (d). Based on these simulations, one can conclude that $h_s$ and ${\theta }_s$ exhibit stable behaviors.

Figures 4 and 5 depict the time histories of the found approximative solutions $X$ and $\theta$ , up to the third degree of $\epsilon$ . These graphs are created using the previously utilized data at various settings ${\omega }_s$ . We may infer from the portrayed curves in these plots that the behavior of the solution $X$ as a function of the dimensionless time $\tau$ take the form of decaying waves, which are represented by wave packets with decreasing amplitudes until they reach a certain value, and then they have the form of stationary flow, see drawn in Figures 4(a) and 5(a). However, the amplitudes of waves illustrating the time behavior of the solution $\theta$ increase to specified values, then they have stationary till the end of the investigated time interval, as displayed in Figures 4(b) and 5(b).

Figure 4.

Time history of the asymptotic solutions $\mathit{X}\left(\mathit{\tau }\right)$ and $\mathit{\theta }\left(\mathit{\tau }\right)$ when the frequency ${\mathit{\omega }}_{\mathbf{1}}\left(=\mathbf{15.5},\mathbf{17.9},\mathbf{20}\right)$ has various values.

Figure 5.

Time history of asymptotic solutions $X\left(\tau \right)$ and $\theta \left(\tau \right)$ when the frequency ${\omega }_2\left(=2.9,3.2,3.6\right)$ has various values.

The analytic solutions and the numerical ones of the original EOM are plotted graphically, as seen in Figure 6. A deeper examination of this figure reveals excellent consistency between the two solutions, which is an indication of the approximations’ high accuracy.

Figure 6.

Comparison between the numerical solutions and the analytical ones at ${\mathit{\omega }}_{\mathbf{1}}=\mathbf{20}$ and ${\mathit{\omega }}_{\mathbf{2}}=\mathbf{3.6}$ when ${\mathit{c}}_{\mathbf{1}}=\mathbf{0.005}$ and ${\mathit{c}}_{\mathbf{2}}=\mathbf{0.012}$ .

The graphs of phase plane for the represented solutions in Figures 5 and 6 have been displayed in Figures 7 and 8. As may be seen from the drawn curves in the planes $\mathit{X}{\mathit{X}}^{\prime }$ and $\mathit{\theta }{\mathit{\theta }}^{\prime }$ that they have directed spiral paths into only a certain point, which indicates their stable behavior.

Figure 7.

Phase portrait in the plane $X{X}^{\prime }$ at different values of the frequency ${\omega }_1\left(=15.5,17.9,20\right)$ .

Figure 8.

Phase portrait in the plane $\theta {\theta }^{\prime }$ at different values of the frequency ${\omega }_1\left(=15.5,17.9,20\right)$ .

To continue the examination of the stability of the studied system, the relationships between the ${\mathit{h}}_{\mathit{s}}\left(\mathit{s}=\mathbf{1},\mathbf{2}\right)$ and the ${\mathit{\theta }}_{\mathit{s}}$ have been portrayed in Figures 9 and 10 according to the numerical outcomes of the system of equation (40), which is considered one of the best ways to clarify the behavior of the system. It is observed that all curves take the form of spiral trajectories towards one point, which indicates that this system has stable behavior during the examined time interval.

Figure 9.

The phase plane ${\theta }_1{h}_1$ at the mentioned values of the frequencies ${\omega }_1$ and ${\omega }_2$ .

Figure 10.

The phase plane ${\theta }_2{h}_2$ at the aforementioned values of the frequencies ${\omega }_1$ and ${\omega }_2$ .

Stable-state approaches

The examination of the solutions for the studied dynamical system’s steady-state situation is the main topic of this section. This scenario is where the solutions would appear when the transitory processes cease to exist. In the present scenario, we take into account the zero value of $\frac{\mathit{d}{\mathit{h}}_{\mathit{s}}}{\mathit{d}\mathit{\tau }}=\frac{\mathit{d}{\mathit{\theta }}_{\mathit{s}}}{\mathit{d}\mathit{\tau }}=\mathbf{0}\left(\mathit{s}=\mathbf{1},\mathbf{2}\right)$ . ²⁷ Therefore, by using the system of equation (40), we may generate the following form of algebraic equations:

$\begin{array}{ll}{f}_1\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_1+2h_1\left({\sigma }_1-S_1\right)-\frac{3}{2}\frac{{\omega }^2\left({\omega }^2-1\right)}{\left(4{\omega }^2-1\right)}{h}_1{h}_2^2-3\alpha {\zeta }_s^2{h}_1=0,& \\ f_1\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_1-h_1\left(c_1+2S_2\right)=0,& \frac{f_2}{\omega }\mathit{sin}\hspace{0.17em}{\theta }_2-c_2{h}_2=0.\\ \frac{f_2}{\omega }\mathit{cos}\hspace{0.17em}{\theta }_2+2h_2{\sigma }_2+\frac{\omega \left({\omega }^2-1\right)}{2\left(4{\omega }^2-1\right)}\left(\frac{\left(8{\omega }^2+1\right)}{4}{h}_2^3-3h_2{h}_1^2\right)=0,& \end{array}$

(41)

According to an examination of the previous system (41), the below frequency response functions are produced when the altered phases ${\theta }_s$ are removed from its equations.

$$\begin{gathered} f_1^2=h_1^2\left\{{\left[3\alpha {\zeta }_s^2+\frac{3}{2}\frac{{\omega }^2\left({\omega }^2-1\right)}{\left(4{\omega }^2-1\right)}{h}_2^2-2\left({\sigma }_1-S_1\right)\right]}^2+{\left(c_1+2S_2\right)}^2\right\}, \\ {f}_2^2={\omega }^2{h}_2^2\left\{{\left[2{\sigma }_2+\frac{\omega \left({\omega }^2-1\right)}{2\left(4{\omega }^2-1\right)}\left(\frac{\left(8{\omega }^2+1\right)}{4}{h}_2^2-3h_1^2\right)\right]}^2+c_2^2\right\}. \end{gathered}$$

(42)

It is crucial to remember that all potential steady-state solutions may be represented for the amplitudes $h_s$ at ${\sigma }_1=p_1-1$ and ${\sigma }_2=p_2-\omega$ , within the context of system (42), as depicted in Figure 11, in which it has been calculated for certain data of ${\sigma }_s\left(s=1,2\right)$ . The blue and dashed red curves show, respectively, the roots of the first and second equations of system (42). In other words, these curves show where the roots of these equations are located, and the spots where the two curves cross reflect the fixed points that are equivalent to the solutions at the steady-state scenario. These points provide an irrefutable formula for determining the axial amplitudes and regulating vibration in the steady-state. Furthermore, steady-state oscillations could be either stable or unstable. Small black circles are used to denote the stable fixed points, whereas small yellow ones are provided to denote the unstable case of fixed points.

We are going to examine the dynamical system’s behavior in a location close to the fixed points for further investigation of its stability. For this goal, we take into account very slight changes to the amplitudes and phases as follows.

$h_1=h_{10}+h_{11}, h_2=h_{20}+h_{21}, {\theta }_1={\theta }_{10}+{\theta }_{11}, {\theta }_2={\theta }_{20}+{\theta }_{21},$

(43)

where $h_{10},h_{20},{\theta }_{10},$ and ${\theta }_{20}$ indicate the steady-state solutions, while $h_{11},h_{21},{\theta }_{11},$ and ${\theta }_{21}$ signify the minor perturbations in comparison to $h_{10},h_{20},{\theta }_{10},$ and ${\theta }_{20}$ . After linearization, we may obtain the next system according to the substitution of (43) into (40).

$$\begin{gathered} {\mathit{h}}_{\mathbf{10}}\frac{\mathit{d}{\mathit{\theta }}_{\mathbf{11}}}{\mathit{d}\mathit{\tau }}={\mathit{h}}_{\mathbf{11}}\left({\mathit{\sigma }}_{\mathbf{1}}-{\mathit{S}}_{\mathbf{1}}\right)-\frac{\mathbf{3}}{\mathbf{4}}\frac{{\mathit{\omega }}^{\mathbf{2}}\left({\mathit{\omega }}^{\mathbf{2}}-\mathbf{1}\right)}{\left(\mathbf{4}{\mathit{\omega }}^{\mathbf{2}}-\mathbf{1}\right)}{\mathit{h}}_{\mathbf{20}}{\mathit{h}}_{\mathbf{11}}\left({\mathit{h}}_{\mathbf{20}}+\mathbf{2}{\mathit{h}}_{\mathbf{10}}\right)-\frac{\mathbf{3}}{\mathbf{2}}\mathit{\alpha }{\mathit{\zeta }}_{\mathit{s}}^{\mathbf{2}}{\mathit{h}}_{\mathbf{11}}-\frac{{\mathit{f}}_{\mathbf{1}}}{\mathbf{2}}{\mathit{\theta }}_{\mathbf{11}}\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\mathit{\theta }}_{\mathbf{10}}, \\ \frac{\mathit{d}{\mathit{h}}_{\mathbf{11}}}{\mathit{d}\mathit{\tau }}=-{\mathit{h}}_{\mathbf{11}}\left(\frac{{\mathit{c}}_{\mathbf{1}}}{\mathbf{2}}+{\mathit{S}}_{\mathbf{2}}\right)+\frac{{\mathit{f}}_{\mathbf{1}}}{\mathbf{2}}{\mathit{\theta }}_{\mathbf{11}}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\mathit{\theta }}_{\mathbf{10}}, \\ {\mathit{h}}_{\mathbf{20}}\frac{\mathit{d}{\mathit{\theta }}_{\mathbf{21}}}{\mathit{d}\mathit{\tau }}={\mathit{h}}_{\mathbf{21}}{\mathit{\sigma }}_{\mathbf{2}}+\frac{\mathit{\omega }\left({\mathit{\omega }}^{\mathbf{2}}-\mathbf{1}\right)}{\mathbf{4}\left(\mathbf{4}{\mathit{\omega }}^{\mathbf{2}}-\mathbf{1}\right)}{\mathit{h}}_{\mathbf{21}}\left\{\frac{\left(\mathbf{8}{\mathit{\omega }}^{\mathbf{2}}+\mathbf{1}\right)}{\mathbf{4}}{\mathit{h}}_{\mathbf{20}}^{\mathbf{2}}-\mathbf{3}{\mathit{h}}_{\mathbf{10}}^{\mathbf{2}}-\mathbf{6}{\mathit{h}}_{\mathbf{10}}{\mathit{h}}_{\mathbf{20}}\right\}-\frac{{\mathit{f}}_{\mathbf{2}}}{\mathbf{2}\mathit{\omega }}{\mathit{\theta }}_{\mathbf{21}}\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\mathit{\theta }}_{\mathbf{20}}, \\ \frac{\mathit{d}{\mathit{h}}_{\mathbf{21}}}{\mathit{d}\mathit{\tau }}=-\frac{\mathbf{1}}{\mathbf{2}}{\mathit{c}}_{\mathbf{2}}{\mathit{h}}_{\mathbf{21}}+\frac{{\mathit{f}}_{\mathbf{2}}}{\mathbf{2}\mathit{\omega }}{\mathit{\theta }}_{\mathbf{21}}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\mathit{\theta }}_{\mathbf{20}}. \end{gathered}$$

(44)

Figure 11.

Solution at the steady-state scenario: (a) three fixed points when ${\sigma }_1=0.01,{\sigma }_2=0.003$ , (b) five fixed points when ${\sigma }_1=0.02,{\sigma }_2=0.009$ and (c) one fixed point when ${\sigma }_1=0.001,{\sigma }_2=0.003$ .

We may express the solutions of this system as in a form of linear combination of $k_j{e}^{\lambda \tau }$ where $k_j\left(j=1,2,3,4\right)$ are constants and $\lambda$ is the unknown perturbation’s eigenvalue, respectively, that predicated on the tiny of the perturbations functions $h_{11},h_{21},{\theta }_{11},$ and ${\theta }_{21}$ . In light of the above, the fixed points in (44) are stable asymptotically when the real portions for the roots of the next characteristic equation are negative.

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

(45)

where ${\Gamma }_j$ are expressed in Appendix I.

According to Routh– Hurwitz criteria, ⁴² the following conditions must be met in order for fixed points to be stable:

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

(46)

Examination of stability

This section’s goal is to investigate the examined model’s nonlinear stability using the approach of nonlinear stability. ⁴³ On the basis of the foregoing, the dynamical motion is impacted by the external harmonic forces $F_s\left(t\right);\left(s=1,2\right)$ and moment $M\left(t\right)$ . Furthermore, to simulate the system of nonlinear equation (40), the stability conditions are used. It was discovered that a number of stability criteria factors, including frequency $\omega$ , damping coefficients $c_s$ , and detuning parameters ${\sigma }_s$ , significantly contribute to instability. To display the associated diagrams of the system’s stability, suitable regulations are used with a variety of system (40) parameters. The amplitudes $h_s$ for the waves describing the solutions time histories are shown in separate parametric regions, and their characteristics are graphed throughout the phase plane’s trajectories. Moreover, the potential fixed points have been plotted. Figures 12–20 examine how these points vary depending on the detuning parameters ${\sigma }_s$ .

Figure 12.

The FRC of the amplitudes $h_s\left({\sigma }_1\right);s=1,2$ at different values of $C_1$ at ${\sigma }_2=0$ and $C_2=0.08$ .

Figure 13.

The FRC of the amplitudes ${\mathit{h}}_{\mathit{s}}\left({\mathit{\sigma }}_{\mathbf{1}}\right);\mathit{s}=\mathbf{1},\mathbf{2}$ at different values of ${\mathit{C}}_{\mathbf{1}}$ when ${\mathit{\sigma }}_{\mathbf{2}}=\mathbf{0.001}$ and ${\mathit{C}}_{\mathbf{2}}=\mathbf{0.08}$ .

Figure 14.

The FRC of the amplitudes $h_s\left({\sigma }_1\right);s=1,2$ at different values of $C_1$ when ${\sigma }_2=-0.001$ and $C_2=0.08$ .

Figure 15.

The FRC of the amplitudes $h_s\left({\sigma }_2\right);s=1,2$ at different values of $C_1$ when ${\sigma }_1=0$ and $C_2=0.08$ .

Figure 16.

The FRC of the amplitudes $h_s\left({\sigma }_1\right);s=1,2$ at ${\sigma }_2=0$ and $C_1=0.01$ when $C_2$ has distinct values.

Figure 17.

The FRC of the amplitudes $h_s\left({\sigma }_1\right);s=1,2$ at ${\sigma }_2=0.001$ and $C_1=0.01$ when $C_2$ has distinct values.

Figure 18.

The FRC of the amplitudes $h_s\left({\sigma }_1\right);s=1,2$ at ${\sigma }_2=-0.001$ and $C_1=0.01$ when $C_2$ varies.

Figure 19.

The FRC of the adjusted amplitudes $h_s\left({\sigma }_2\right);s=1,2$ at ${\sigma }_1=0$ and $C_1=0.01$ when $C_2$ varies.

Figure 20.

The FRC of the modified amplitudes $h_s\left({\sigma }_1\right);s=1,2$ at ${\sigma }_2=0,$ $C_1=0.01$ and $C_2=0.08$ when ${\omega }_1$ varies.

The displayed curves in Figure 12 are calculated at the zero value of ${\sigma }_2$ and $C_2=0.08$ , in addition to the various values of $C_1$ . It must be noted that the dashed and solid curves are used to depict the unstable and stable fixed points. A more detailed examination for the parts of this figure shows that the amplitude $h_1$ vary with the change of the values of $C_1$ , as indicated in Figure 12(a), while the variation of $h_1$ according to the values of the same parameter seems to be slightly as observed in Figure 12(c)–(e). Moreover, the ranges containing the zones of stability and instability are ${\sigma }_1\le 0.001$ and $0.001<{\sigma }_1$ , respectively, when $C_1=0.01$ for the blue color. On the other hand, when $C_1=0.05$ and $C_1=0.09$ , the stability and instability areas are found at the ranges ${\sigma }_1\le 0$ , ${\sigma }_1\le -0.001$ and $0<{\sigma }_1$ , $-0.001<{\sigma }_1$ , respectively, for the red and black curves. This graph has eight peak fixed points (PFP) and six critical fixed points (CFP). The distribution of these points covered all portions of Figure 12 and their locations are stated in Table 1.

Table 1.

Displays the location of the PFP and the CFP for the FRC at ${\mathit{C}}_{\mathbf{2}}=\mathbf{0.08}$ for different values of ${\mathit{C}}_{\mathbf{1}}$ .

Figure	$C_1$	PFP	CFP	${\sigma }_s\left(s=1,2\right)$
Figure 12(a)	0.01	(0.001,0.0437)	(0.001,0.0437)	${\sigma }_2=0$
	0.05	(0,0.01003)	(0,0.01003)
	0.09	(‒0.001,0.0053)	(‒0.001,0.0053)
Figure 12(b)	0.01	(0.001,0.01294)	(0.001,0.01294)
Figure 12(c)	0.05	(‒0.0019,0.0142045)	(0,0.142032)
		(0,0.142032)
		(0.00383,0.014244)
Figure 12(d)	0.09	(‒0.001,0.0142045)	(‒0.001,0.0142045)
Figure 13(a)	0.01	(0.001,0.044)	(0.001,0.044)	${\sigma }_2=0.001$
	0.05	(0,0.01011)	(0,0.01011)
	0.09	(‒0.001,0.005708)	(‒0.001,0.005708)
Figure 13(b)	0.01	(0.0000685,0.01421)	(0.001,0.01367)
		(0.001,0.01367)
		(0.002177,0.01411)
	0.05	(0,0.01401)	(0,0.01401)
	0.09	(‒0.001,0.01397)	(‒0.001,0.01397)
Figure 14(a)	0.01	(0.001,0.04297)	(0.001,0.04297)	${\sigma }_2=-0.001$
	0.05	(0,0.0104)	(0,0.0104)
	0.09	(‒0.001,0.005561)	(‒0.001,0.005561)
Figure 14(b)	0.01	(0.001,0.01205)	(0.001,0.01205)
	0.05	(0,0.01394)	(0,0.01394)
	0.09	(‒0.001,0.01398)	(‒0.001,0.01398)
Figure 15(a)	0.01	(0.0012,0.02833)	(‒0.001,0.02873)	${\sigma }_1=0$
Figure 15(b)	0.05	(0,0.009602)	(‒0.002,0.00961)
Figure 15(c)	0.09	(‒0.00023,0.005484)	(‒0.003,0.00549)
Figure 15(d)	0.01	(0.0009941,0.01419)	(‒0.001,0.01327)
	0.05	(‒0.000249,0.01427)	(‒0.002,0.01368)
	0.09	(0.000376,0.01419)	(‒0.003,0.01226)

Moreover, Figures 13 and 14 show the FRC for both ${\mathit{h}}_{\mathbf{1}}\left({\mathit{\sigma }}_{\mathbf{1}}\right)$ and ${\mathit{h}}_{\mathbf{2}}\left({\mathit{\sigma }}_{\mathbf{1}}\right)$ at ${\mathit{\sigma }}_{\mathbf{2}}=\mathbf{0.001}$ and ${\mathit{\sigma }}_{\mathbf{2}}=-\mathbf{0.001}$ , respectively, with the same values of ${\mathit{C}}_{\mathbf{1}}$ and ${\mathit{C}}_{\mathbf{2}}=\mathbf{0.08}$ . These curves have the same areas of stability and instability of Figure 12. The good impact for various values of when ${\mathit{C}}_{\mathbf{1}}$ on the behavior of ${\mathit{h}}_{\mathit{s}}$ has been described by the plotted curves in Figure 15 when ${\mathit{\sigma }}_{\mathbf{1}}=\mathbf{0}$ and ${\mathit{C}}_{\mathbf{2}}=\mathbf{0.08}$ . It observed that at ${\mathit{C}}_{\mathbf{1}}=\mathbf{0.01},\mathbf{0.05},$ and $\mathbf{0.09}$ , the instability and stability regions are ${\mathit{\sigma }}_{\mathbf{2}}\le -\mathbf{0.001},$ ${\sigma }_2\le -0.002,{\mathit{\sigma }}_{\mathbf{2}}\le -\mathbf{0.003}$ and $-\mathbf{0.001}<{\mathit{\sigma }}_{\mathbf{2}},-\mathbf{0.002}<{\mathit{\sigma }}_{\mathbf{2}},-\mathbf{0.003}<{\mathit{\sigma }}_{\mathbf{2}}$ , respectively, with blue, red, and black colors. The observed PFP and CFP of Figures 13–15 are determined, as mentioned in Table 1. Altering the damping factor ${\mathit{C}}_{\mathbf{2}}$ produces the behavior of the FRC for the system (4). These variations are drawn in Figures 16–18 at ${\mathit{C}}_{\mathbf{1}}=\mathbf{0.01}$ and ${\mathit{\sigma }}_{\mathbf{2}}=\left(\mathbf{0},\mathbf{0.001},-\mathbf{0.001}\right)$ , respectively. The stability and instability areas are ${\mathit{\sigma }}_{\mathbf{1}}\le \mathbf{0.001},{\mathit{\sigma }}_{\mathbf{1}}\le \mathbf{0},{\mathit{\sigma }}_{\mathbf{1}}<-\mathbf{0.001}$ and $.\mathbf{001}<{\mathit{\sigma }}_{\mathbf{1}},\mathbf{0}<{\mathit{\sigma }}_{\mathbf{1}},$ $-\mathbf{0.001}<{\mathit{\sigma }}_{\mathbf{1}}$ , respectively, for the blue, red and black curves.

The graphed curves in Figure 19 shows $h_1\left({\sigma }_2\right)$ and $h_2\left({\sigma }_2\right)$ when ${\sigma }_1=0$ and $C_1=0.01$ , while the drawn ones in Figure 20 presents $h_1\left({\sigma }_1\right)$ and $h_2\left({\sigma }_1\right)$ at ${\sigma }_2=0$ besides the different values of the frequency ${\omega }_1\left(=15.5,17.9,20\right)$ . It is noted that the stable and unstable regions for the plotted curves in Figures 19 and 20 are $-0.001<{\sigma }_1,{\sigma }_1\le 0$ and ${\sigma }_1\le -0.001,0<{\sigma }_1$ , respectively. The position of gained PFP and CFP for Figures 16–20 are determined, in Tables 2 and 3, respectively.

Table 2.

Reveals the location of the PFP and the CFP for the FRC at ${\mathit{C}}_{\mathbf{1}}=\mathbf{0.01}$ when ${\mathit{C}}_{\mathbf{2}}$ has different values.

Figure	${\mathit{C}}_{\mathbf{2}}$	PFP	CFP	${\sigma }_s\left(s=1,2\right)$
Figure 16(a)	0.08	(0.001,0.04326)	(0.001,0.04326)	${\sigma }_2=0$
	0.25	(0.001073,0.03945)	(0,0.03226)
	0.41	(0.001476,0.03896)	(-0.001,0.01541)
Figure 16(b)	0.08	(0.001,0.01294)	(0.001,0.01294)
Figure 16(c)	0.25	(0.0006097,0.004696)	(0,0.004719)
Figure 16(d)	0.41	(0.001177,0.002832)	(-0.001,0.002841)
Figure 17(a)	0.08	(0.001,0.04385)	(0.001,0.04385)	${\mathit{\sigma }}_{\mathbf{2}}=\mathbf{0.001}$
	0.25	(0.001073,0.03901)	(0,0.03211)
	0.41	(0.0008714,0.03896)	(-0.001,0.01523)
Figure 17(b)	0.08	(-0.001,0.0142)	(0.001,0.01369)
		(0.001,0.01369)
		(0.002387,0.01412)
Figure 17(c)	0.25	(0,0.004734)	(0,0.004734)
		(0.001398,0.004724)
		(0.002283,0.004732)
Figure 17(d)	0.41	(0.0004135,0.002841)	(-0.001,0.00284)
		(0.0008667,0.002839)
		(0.002906,000284)
Figure 18(a)	0.08	(0.001,0.042971)	(0.001,0.042971)	${\mathit{\sigma }}_{\mathbf{2}}=-\mathbf{0.001}$
	0.25	(0.001274,0.003974)	(0,0.03241)
	0.41	(0.001073,0.03924)	(-0.001,0.01555)
Figure 18(b)	0.08	(0.001,0.01201)	(0.001,0.01201)
Figure 18(c)	0.25	(0.001041,0.004651)	(0,0.004686)
Figure 18(d)	0.41	(0.001177,0.002823)	(-0.001,0.002837)
Figure 19(a)	0.08	(0.0009941,0.0283)	(-0.001,0.02877)	${\sigma }_1=0$
	0.25	(0.0009941,0.03092)	(-0.001,0.031)
	0.41	(0.0009941,0.03143)	(-0.001,0.0315)
Figure 19(b)	0.08	(0.0009941,0.0142)	(-0.001,0.01324)
	0.25	(0.0009941,0.007212)	(-0.001,0.006839)
	0.41	(0.0009941,0.004865)	(-0.001,0.004545)

Table 3.

Outlines the PFP and CFP locations for the FRC at $C_1=0.01$ , $C_2=0.08$ , and when ${\omega }_1$ has various values.

Figure	${\omega }_1$	PFP	CFP	${\sigma }_s\left(s=1,2\right)$
Figure 20(a)	15.5	(0,0.005415)	(0,0.005415)	${\sigma }_2=0$
	17.9	(0,0.009815)	(0,0.009815)
	20.0	(0,0.03138)	(0,0.03138)
Figure 20(b)	15.5	---------	(0,0.04491)
	17.9	---------	(0,0.03154)
	20.0	---------	(0,0.02731)

Now, we will explore the stability of the nonlinear amplitudes of system (40) and demonstrate the characteristics of the nonlinear analysis. Consequently, let’s think about the following transformations ⁴⁴ .

$A_s=\left[P_s\left({\tau }_2\right)+i{Q}_s\left({\tau }_2\right)\right]e^{i{\stackrel{\sim }{\sigma }}_s{\tau }_1};\left(s=1,2\right),$

(47)

where

$P_s=\epsilon p_s,Q_s=\epsilon q_s.$

Making use of (40) and (47), then isolating the real imaginary terms to obtain

$\frac{d{p}_1}{d\tau }=q_1\left({\sigma }_1-\frac{3}{2}\alpha {\zeta }_s^2-S_1-S_2\right)-\frac{1}{2}{c}_1{p}_1-\frac{3{\omega }^2\left({\omega }^2-1\right)q_1}{4\left(4{\omega }^2-1\right)}\left(p_2^2+q_2^2\right),\frac{d{q}_1}{d\tau }=\left(S_1+S_2-{\sigma }_1\right)p_1-\frac{1}{2}{c}_1{q}_1+\frac{3{\omega }^2\left({\omega }^2-1\right)p_1}{4\left(4{\omega }^2-1\right)}\left(p_2^2+q_2^2\right)+\frac{3}{2}\alpha {\zeta }_s^2{p}_1-\frac{f_1}{4},\frac{d{p}_2}{d\tau }={\sigma }_2{q}_2-\frac{1}{2}{c}_2{p}_2+\frac{\omega \left({\omega }^2-1\right)q_2}{4\left(4{\omega }^2-1\right)}\left[\left(8{\omega }^2+1\right)\left(p_2^2+q_2^2\right)-4\left(p_1^2+q_1^2\right)\right],\frac{d{q}_2}{d\tau }=-{\sigma }_2{p}_2-\frac{1}{2}{c}_2{q}_2-\frac{\omega \left({\omega }^2-1\right)p_2}{4\left(4{\omega }^2-1\right)}\left[\left(8{\omega }^2+1\right)\left(p_2^2+q_2^2\right)-4\left(p_1^2+q_1^2\right)\right]-\frac{f_2}{4\omega }.$

(48)

Over the time span of the study, the adjusted amplitudes were confirmed in a variety of parametric areas, and the amplitudes’ characteristics have been plotted through the curves of the phase plane, as shown in Figures 21–24, in which the previous data have been used. The variation of the modified phases $p_s$ and $q_s$ via $\tau$ are presented in Figures 21 and 22 to show their fluctuations when ${\omega }_1$ and ${\omega }_2$ take various values. It is observed that the explored waves of $p_1\left(\tau \right)$ and $q_1\left(\tau \right)$ oscillates at the beginning of the examined time interval till to reach to steady behavior at the end of this interval. In other words, their behavior have the forms of decay procedure, as seen in Figures 21(a) and (b) and 22(a) and (b). On the other side, the plotted curves in Figures 21(c) and (d) and 22(c) and (d) for the functions $p_2\left(\tau \right)$ and $q_2\left(\tau \right)$ have a sharp decline through the first quartal of the examined interval, and then their behaviors have a steady manner till the end of this interval. On the other hand, the projections of the modulated trajectories on the planes $p_1{q}_1$ and $p_2{q}_2$ are shown in parts of Figures 23 and 24. These trajectories have directed spiral curves towards a fixed point, revealing that these amplitudes have stationary behavior.

Figure 21.

The behavior of ${\mathit{p}}_{\mathit{s}}\left(\mathit{\tau }\right)$ and ${\mathit{q}}_{\mathit{s}}\left(\mathit{\tau }\right);\left(\mathit{s}=\mathbf{1},\mathbf{2}\right)$ when ${\mathit{\omega }}_{\mathbf{1}}\left(=\mathbf{15.5},\mathbf{17.9},\mathbf{20}\right)$ has various values.

Figure 22.

The behavior of $p_s\left(\tau \right)$ and $q_s\left(\tau \right);\left(s=1,2\right)$ when ${\omega }_2\left(=2.9,3.2,3.6\right)$ has different values.

Figure 23.

The phase plane trajectories in the plans $p_1{q}_1$ and $p_2{q}_2$ at various values of ${\omega }_1\left(=15.5,17.9,20\right)$ .

Figure 24.

The phase plane trajectories in the plans $p_1{q}_1$ and $p_2{q}_2$ at various values of ${\omega }_2\left(=2.9,3.2,3.6\right)$ .

Performance of piezoelectric device

The main purpose of this part is to look into how the piezoelectric device affects the behavior of the considered dynamical system and to assess its importance in the generation of electrical energy. This device consists of dielectric materials, which is susceptible to polarization as a result of mechanical stress brought on by the dynamical model’s vibrations. In our scenario, polarization of these materials produces an electric field. As a result, the required electrical energy is created from mechanical energy utilizing the dynamical model and this device. The energy-collecting gadget produces electrical energy, which can be used for a variety of tasks, such as habitat monitoring (temperature, light, and humidity), medical remote sensing, structural monitoring, such as defibrillator monitoring and pacemaker, aerospace applications and military, and emergency medical response monitoring.

In order to get high performance, we’re going to look at how different values of the used parameters effect on the system. Figures 25–27 illustrate the time histories of the piezoelectric’s output voltages, powers, and energy, respectively. An inspection of the parts of Figure 25 reveals that we have periodic decay of standing waves throughout the whole-time interval. The amplitudes of these waves decrease with the increase of the damping parameter ${\mathit{C}}_{\mathbf{1}}$ and increase with the values of ${\mathit{F}}_{\mathbf{1}}$ , respectively, as explored in Figure 25(a) and (b). The impact of various values of ${\mathit{F}}_{\mathbf{2}},{\mathit{\omega }}_{\mathbf{1}},$ and ${\mathit{\omega }}_{\mathbf{2}}$ produces periodic decaying waves during the examined interval, as noted in Figure 25(c)–(e). The positive influence of the parameters ${\mathit{C}}_{\mathbf{1}},{\mathit{F}}_{\mathbf{1}},$ and ${\mathit{F}}_{\mathbf{2}}$ on the behavior of the produced power as a result of the piezoelectric device, is drawn in Figure 26. Analyzing the plotted curves in Figure 26 demonstrated that the output power increases and decreases with the increase of the values of ${\mathit{C}}_{\mathbf{1}}$ and $\left({\mathit{F}}_{\mathbf{1}},{\mathit{F}}_{\mathbf{2}}\right)$ , respectively. The same analyses can be applied on the produced curves in Figure 27 for the output energy $\mathit{E}$ , in which this energy decreases with the increase of ${\mathit{C}}_{\mathbf{1}}$ and increases with the raise values of ${\mathit{F}}_{\mathbf{1}}$ and ${\mathit{F}}_{\mathbf{2}}$ .

Figure 25.

The voltage’s behavior via $\tau$ at the different values of $C_1,F_1,F_2,{\omega }_1,$ and ${\omega }_2$ .

Figure 26.

The temporal behavior of the output power at different values of $C_1,F_1,$ and $F_2$ .

Figure 27.

The energy’s time histories at different values of $C_1,F_1,$ and $F_2$ .

Comparisons with benchmark work with other previous works, for example, Ref. 23 and Ref. 45 is provided in Table 4 to include the design, the state of the pivot point, the number of DOF, the obtained solutions, stability regions, fixed points, CFP, PFP, nonlinear analysis, as well as the EH device.

Table 4.

Comparison between the present work and previous ones in Refs. 23 and 45.

No.	Comparison items	The present study	The work in Ref. 23	The work in Ref. 45
1	Degrees-of-freedom	2DOF	3DOF	2DOF
2	Linearity of the spring	Nonlinear	Linear	Linear
3	Pivot point	In motion	In motion	Fixed
4	Analytic solution	The analytic solution has been obtained	The analytic solution has been obtained	The analytic solution has been obtained
5	Stability	The stability and instability areas have been investigated according to Routh-Hurwitz criteria	Not achieved	Not achieved
6	Fixed points	All fixed points have been examined whether stable or not	Not achieved	These pints have been checked
7	PFP and CFP	All arising PFP and CFP have been tabulated and classified	Not achieved	Not achieved
8	Nonlinear analysis	Nonlinear analysis has been examined	Not achieved	Not achieved
9	EH device	The influence of a piezoelectric energy harvesting device has been taken into account	Not achieved	Not achieved

Conclusions

A novel dynamical system with a 2DOF nonlinear damped oscillating SP system that is coupled to a piezoelectric EH device has been examined, in which its suspension point moves in a Lissajous route with a constant angular velocity. The following results have been listed:

1. The controlling motion’s system has been derived utilizing Lagrange’s equations and solved analytically applying using the MST up to higher approximation to obtain novel analytic solution.

2. The achieved analytic solution has been verified through comparison with the numerical one of the original system of motion, to reveal the great accuracy of the analytic solution and to demonstrate the validity of the MST.

3. The solvability conditions are established once the produced secular words have been removed. Consequently, every resonance case has been categorized. The modulation equations for two of these resonances have been explored concurrently in view of the system’s adjusted phases.

4. From a visual standpoint numerous diagrams showing the temporal histories of the achieved solutions have been analyzed. The gained solutions for the steady-state case have been used to guide the investigation of the nonlinear stability of the modulation equations.

5. The effects of multiple parameters, including damping coefficients, various frequencies, and excitation amplitudes, on the behavior of the system of modulation equations have been thoroughly investigated and analyzed on the outputs of the EH device.

6. To show how the acted-upon parameters influence the system’s behavior, many zones of stability and instability for the FRC have been depicted.

7. The stability and instability zones are ${\sigma }_1\le 0.001$ and $0.001<{\sigma }_1$ at ${\sigma }_2=0,C_2=0.08,$ and $C_1=0.01$ , as illustrated in Figure 12.

8. The effect of ${\mathit{\sigma }}_{\mathbf{2}}$ on ${\mathit{h}}_{\mathbf{1}}$ seems to be slightly due to the system of equation (42), as indicated in Figure 15(b) and (c).

9. The amplitudes are inversely correlated with the damping coefficients.

10. The piezoelectric device generates $0.02$ volt, $7*10^{-7}$ Watt, and $0.002$ Watt. hour at $F_1=0.09,$ $F_2=0.09,$ and $C_1=0.01$ , as graphed in Figures 25–27.

11. Furthermore, the stability of the nonlinear amplitudes of the modulation equations and the characteristics of the nonlinear analysis have been demonstrated in light of the use of complex transformations of the system’s amplitudes.

The piezoelectric device, which is connected to the vibrating system, produces electrical energy, as mentioned above, which is used in some practical applications, like as a source of electricity for sensors, for charging electronic devices, and for other medicinal uses.