Influence of the Motion of a Spring Pendulum on Energy-Harvesting Devices

Abstract

Energy harvesting is becoming more and more essential in the mechanical vibration application of many devices. Appropriate devices can convert the vibrations into electrical energy, which can be used as a power supply instead of ordinary ones. This study investigated a dynamical system that correlates with two devices, namely a piezoelectric device and an electromagnetic one, to produce two novel models. These devices are connected to a nonlinear damping spring pendulum with two degrees of freedom. The damping spring pendulum is supported by a point moving in a circular orbit. Lagrange’s equations of the second kind were utilized to obtain the equations of motion. The asymptotic solutions of these equations were acquired up to the third approximation using the approach of multiple scales. The comparison between the approximate and the numerical solutions reveals high consistency between them. The steady-state solutions were investigated, and their stabilities were checked. The influences of excitation amplitudes, damping coefficients, and the different frequencies on energy-harvesting device outputs are examined and discussed. Finally, the nonlinear stability analysis of the modulation equations is discussed through the stability and instability ranges of the frequency response curves. The work is significant due to its real-life applications, such as a power supply of sensors, charging electronic devices, and medical applications.

1. Introduction

Energy harvesting (EH) has been a substantial aspect of research work in the previous few years. It transforms the surrounding energy presented on Earth into electrical power to drive autonomous electronic devices or circuits [1,2]. This energy can be harnessed from solar energy, thermal energy, and the most vital source for harnessing, kinetic energy, especially from vibrational motion [3]. EH depending on vibration can drive devices like sensors and wireless electronics as well as monitoring applications [4]. Its mechanism’s models depend on the uses of electromagnetic, piezoelectric, and hybrid electromagnetic devices.

In [5], Glynne-Jones et al. designed an electromagnetic generator to convert the environment vibrational energy into an electrical one to power sensor devices. The response of a harvester of energy is investigated in [6]. In [7], the authors computed the coupling due to the electromagnetic energy harvesters and deduced the resonant harvesters’ output voltage. A prototype was designed in [8] to harvest energy from the speed of passing vehicles and control their speed simultaneously.

In [9], the authors investigated a self-powered and self-operating tracking system based on an electromagnetic hybrid friction blue energy-harvesting module, theoretically and experimentally. The unique spin of the designed gyro enables the system to search for the energy of low-frequency wave and non-uniform vibration with high efficiency. A hybrid designed propeller wheel was reported in [10] with a nano-coupled triboelectric generator and an electromagnetic generator that can harvest the wind’s energy even in a weak airflow. At the same time, Fan et al. [11] harvested energy from low frequency. They used a piezoelectric energy harvester (PEH) with stoppers and tuned the device at monostable mode to avoid problems in a bistable or tristable system. A multi-mode PEH with analytical modeling is investigated in [12] to produce different close peaks of the voltage output from the vibrational sources. In [13], theoretical and experimental broadband PEH studies depend on kinetic energy with a triple well potential generated by the magnetic field. The mathematical model is examined to show the response of tristable energy generators.

A comparison between three nonlinear planar springs, namely bi-leg, quad-leg, and pent-leg designs, is presented in [14] to examine the influence of the interaction of a vibration mode to extend the operating frequency range of the PEH. In [15], the authors answered the question: How can the electricity be generated from a mechanical structure driven by vapor momentum? They investigated a design of an energy harvester that experimentally depends on vapor-induced vibration for the heat pipe application. The velocity of the vapor impacted the induced power of the device.

A report on a cost-effective wearable PEH using electrophoresis, conductive fabrics, and coated with nickel-copper-polyester is presented in [16]. The finite element method (FEM) was used in [17] to obtain the numerical results of the effectiveness of the piezoelectric layer’s length and thickness on the von Mises stress. The influence of increasing the seismic mass’s thickness on the cantilever structure’s resonance frequency was studied in [18], in which the FEM was applied to achieve numerical results. It is observed that the resonance frequency of the cantilever decreased with increasing size of the seismic mass. The authors [19] presented the design, manufacture, and testing of an underground double-axis PEH based on the piezoelectric stack.

A modified numerical method was utilized in [20] to obtain the solutions of second-order boundary value problems of ordinary differential equations (ODEs) directly for Dirichlet and mixed cases. In addition, a hybrid single-step implicit block method was used to solve the nonlinear differential equation of a circular sector oscillator in [21]. Significant success was achieved in solving the stiff ODEs using the hybrid single-step implicit block method.

A new flexoelectric energy harvester and broadband piezoelectric in the presence of an applied axial compressive load with an excitation base is investigated in [22], where the energy harvester operates in either a pre-torsion configuration or a post-configuration, depending on the axial load’s amplitude. Moreover, a converter is used in [23] to transform the piezoelectric vibration into electricity as a power source for wireless electronics. A simple piezoelectric spring configuration depending on a standard binder clip structure is studied in [24]. The harvester consists of a spring’s pendulum that converts the mass’s energy into an electrical one in the piezoelectric transducer. Two degrees of freedom (DOF) mechanical model with developed glands was examined in [25] by simulating the impact behavior as linear stiffness piecewise. A new horizontal asymmetric U-shaped vibration-based PEH is presented in [26], in which the nonlinearity is combined with multimodality.

Erturk et al. [27] examined the high-energy orbits in the piezo-magneto-elastic energy harvester theoretically and experimentally over the excitation frequencies area. A comparison between an electromagnetic PEH and a PEH with distinct interface circuits is found in [28]. Furthermore, it is shown in [29] that the effect of the backward electromechanical coupling on the mechanical vibrations of the energy harvesters can be approximated using the different values of the excitation frequency and damping coefficient.

The periodic and quasi-periodic vibration of a dynamical model connected to a piezoelectric circuit was studied in [30]. The approximate solutions were obtained using the complexification-averaging method. A special EH designed system was examined asymptotically using the harmonic balance method in [31], in which the attained results are verified numerically and experimentally. In contrast, the author in [32] designed another nonlinear system that suppressed the vibration and EH at the same time. It is concluded that the starting area of resonance is the best area for suppression.

The dynamical motion of a damped spring pendulum with two or three DOF in the presence of excitation forces and moments is investigated in several works, e.g., [33,34,35,36,37,38,39,40]. The damped motion of a spring pendulum, in which its suspended point moves in a circular trajectory, is investigated in [33]. The authors provided a numerical example for the model’s physical parameters to study the examined model’s chaotic behavior. The planar motion of a linear spring-rigid-body pendulum with fixed pivot is studied in [34], and its generalization is found in [35] and [36] when the fixed point is constrained to move on an ellipse with a constant angular velocity for a linear and nonlinear spring stiffness, respectively. The authors provided other motions as additional examples in [36] for the motion of the fixed point in horizontal and vertical paths. The chaotic motion of an elastic pendulum is explored in [37], the movement of its suspension point in a circular route. It is observed that the approximated system has a bifurcation that leads to chaotic motion via a series of period-doubling bifurcations. The nonlinear motion of a double pendulum system with 3DOF in the presence of one external force and two moments is examined in [38]. The effects of absorber motion on the 3DOF dynamical system is studied in [39]. On the other hand, nonlinear movements of an elastic pendulum under the influence of two perpendicular excitation forces are investigated in [40]. The pivot points of the pendulums considered in [38,39,40] are considered to move in elliptic trajectories. The asymptotic solutions are obtained up to higher order of approximation using the approach of multiple scales (AMS). The solvability conditions and the modulation equations (ME) are achieved in line with the examined resonance cases. The stability and instability areas of the investigated models in [39,40] are explored and represented graphically for suitably selected values of the physical parameters of these models.

The present study investigated two novel dynamical models to convert the vibrational motion of a spring pendulum into electrical energy. Two EH devices, namely piezoelectric and electromagnetic, were connected with the spring separately as two models. The main governing system of motion was derived using Lagrange’s equations, by which the mechanism of the piezoelectric and electromagnetic circuits were used to obtain their corresponding equations [29]. The asymptotic solutions up to the third approximation were obtained utilizing the AMS. External resonance cases between the classified resonance cases were examined. Therefore, ME were obtained to explore the solutions at the steady state and examine the fixed point stability. The numerical results were compared with the approximate ones to reveal the accuracy between them. The influence of damping coefficients and excitation amplitudes on the output voltage, current, and power are represented graphically. Moreover, the response of the excitation frequency on the output power of the considered two dynamical systems was obtained from the energy-harvesting devices. Moreover, the nonlinear stability approach was used to examine the ME through some plots of stability and instability areas for the frequency response of these equations. The importance of this work is due to its significant applications in practical life, such as in a power supply of sensors, charging electronic devices, and various medical applications.

2. Dynamical Modeling

This section is devoted to introducing a complete description of the considered models. Therefore, let us consider the motion of 2DOF dynamical models consisting of nonlinear spring pendulums with the same linear and nonlinear stiffness coefficients $k_1$ and $k_2$, respectively. It is considered that the point of suspensions $O$ moves in a circular path of radius $R$, with an angular velocity $\mathsf{\Omega }$, in which it is subjected to a harmonic torque $M(t)=M_0\hspace{0.17em}\cos {\mathsf{\Omega }}_{2}t$. The pendulum’s other end $Q$ is excited by an external harmonic force $F(t)=F_0\hspace{0.17em}\cos {\mathsf{\Omega }}_{1}t$ that drives the mass $m$ along the pendulum length; see Figure 1. Here, $F_0,M_0$, and ${\mathsf{\Omega }}_1,{\mathsf{\Omega }}_2$ are the amplitudes and frequencies of $F(t)$ and $M(t)$. Let $C_1$ and $C_2$ represent the damping coefficients of the viscous force $C_1\hspace{0.17em}\dot{x}$ and viscous moment $C_2\hspace{0.17em}\dot{\theta }\hspace{0.17em}$, respectively; $l_0$ is the spring’s normal length and $g$ is the gravitational acceleration.

Now, we consider a connection of two devices (piezoelectric and electromagnetic) separately with the considered model, in which the resistive loads of the piezoelectric and electromagnetic circuits are $R_p$ and $R_m$, respectively. Let $c_p$ denote the capacitance of the piezoelectric, $l_m$ the inductance of the coil, and ${\gamma }_j\hspace{0.17em}(j=1,2)$ the linear coupling coefficients for the piezoelectric circuit and electromagnetic one.

Based on the portrayed Figure 1, we can write the cartesian coordinates of the point $Q$ as follows

$$\begin{array}{l}X=R\hspace{0.17em}\cos \hspace{0.17em}\mathsf{\Omega }t\hspace{0.17em}+(l_0+x)\cos \hspace{0.17em}\theta ,\\ Y=R\hspace{0.17em}\sin \hspace{0.17em}\mathsf{\Omega }t\hspace{0.17em}+(l_0+x)\sin \hspace{0.17em}\theta .\end{array}$$

(1)

Henceforth, the potential and kinetic energies of the dynamical system can be written as

$$\begin{array}{l}V=\frac{1}{2}{k}_1\hspace{0.17em}{x}^2+\frac{1}{4}{k}_2\hspace{0.17em}{x}^4-m\hspace{0.17em}g\hspace{0.17em}[R\hspace{0.17em}\cos \hspace{0.17em}\mathsf{\Omega }t\hspace{0.17em}+(l_0+x)\cos \hspace{0.17em}\theta ],\\ T=m\{\frac{1}{2}\hspace{0.17em}[R^2\hspace{0.17em}{\mathsf{\Omega }}^2+{\dot{x}}^2+{(l_0+x)}^2{\dot{\theta }}^2]-\hspace{0.17em}R\hspace{0.17em}\mathsf{\Omega }\hspace{0.17em}\dot{x}\hspace{0.17em}\sin (\mathsf{\Omega }t-\theta )\\ \hspace{1em}+\hspace{0.17em}R\hspace{0.17em}\mathsf{\Omega }\hspace{0.17em}(l_0+x)\dot{\theta }\cos (\mathsf{\Omega }t-\theta )\},\end{array}$$

(2)

where the dot represents the differentiation with respect to time $t$.

In order to obtain the governing equations of motion (EOM) for the mentioned systems, the following Lagrange’s equations are used

$$\begin{array}{l}\frac{d}{dt}(\frac{\partial L}{\partial \dot{x}})-(\frac{\partial L}{\partial x})=Q_x,\\ \frac{d}{dt}(\frac{\partial L}{\partial \dot{\theta }})-(\frac{\partial L}{\partial \theta })=Q_{\theta }.\end{array}$$

(3)

Here $L=T-V$ is the Lagrangian function, and $Q_{\hspace{0.17em}x}\hspace{0.17em}$ and $Q_{\hspace{0.17em}\theta }\hspace{0.17em}$ are the generalized forces corresponding to the generalized standing coordinates $x$ and $\theta$, respectively. Moreover, the mechanism’s equations for the piezoelectric circuit and electromagnetic one can be expressed as follows

$$\begin{array}{l}{c}_p\hspace{0.17em}\dot{v}+\frac{v}{R_p}-{\gamma }_1\hspace{0.17em}\dot{x}=0\\ l_m\hspace{0.17em}\dot{q}+R_m\hspace{0.17em}{l}_m\hspace{0.17em}q-{\gamma }_2\hspace{0.17em}\dot{x}=0\end{array}$$

(4)

where $v$ and $q$ are the voltage and the current of the load resistances $R_p$ and $R_m$, respectively. The generalized forces $Q_{\hspace{0.17em}x}\hspace{0.17em}$ and $Q_{\hspace{0.17em}\theta }\hspace{0.17em}$ have the forms

$$\begin{array}{l}{Q}_x=F_0\hspace{0.17em}\cos {\mathsf{\Omega }}_{1}t\hspace{0.17em}-C_1\dot{x}\hspace{0.17em}-PM,\\ Q_{\theta }=M_0\hspace{0.17em}\cos {\mathsf{\Omega }}_{2}t\hspace{0.17em}-C_2\dot{\theta }\hspace{0.17em},\end{array}$$

where $PM$ represents a coupling term of the two models equals ${\gamma }_{1}v$ for the piezoelectric circuit (model a) or ${\gamma }_{2}q$ for the electromagnetic circuit (model b).

Let us consider the following dimensionless parameters

$$\begin{array}{ll}{\omega }_1^2& =\frac{k_1\hspace{0.17em}}{m},\hspace{0.17em}{\omega }_2^2=\frac{g\hspace{0.17em}}{l}\hspace{0.17em},\hspace{0.17em}{W}^2=\frac{{\omega }_2^2}{{\omega }_1^2}\hspace{0.17em},\hspace{0.17em}r=\hspace{0.17em}\frac{R\hspace{0.17em}}{l},\hspace{0.17em}p=\frac{\mathsf{\Omega }}{{\omega }_1},\\ p_1& =\frac{{\mathsf{\Omega }}_1}{{\omega }_1},\hspace{0.17em}{p}_2=\frac{{\mathsf{\Omega }}_2}{{\omega }_1},\hspace{0.17em}\alpha =\frac{k_2\hspace{0.17em}{l}^2}{{\omega }_1^{2}m},\hspace{0.17em}{\mu }_1=\frac{{\gamma }_1}{{\omega }_1^{2}m\hspace{0.17em}l},\\ {\mu }_2& =\frac{{\gamma }_2}{{\omega }_1^{2}m\hspace{0.17em}l},\hspace{0.17em}PM(={\mu }_1\hspace{0.17em}v\hspace{0.17em},{\mu }_2\hspace{0.17em}q\hspace{0.17em}),\hspace{0.17em}{c}_1=\frac{C_1}{{\omega }_{1}m},\\ c_2& =\frac{C_2}{l^2{\omega }_{1}m},\hspace{0.17em}{f}_1=\frac{F_0}{l{\omega }_1^{2}m},\hspace{0.17em}{f}_2=\frac{M_0}{{\omega }_1^{2}m\hspace{0.17em}{l}^2},\hspace{0.17em}\tau ={\omega }_{1}t,\end{array}$$

(5)

where

$$z_r\hspace{0.17em}+\alpha \hspace{0.17em}{z}_r^3=W^2,\hspace{0.17em}{z}_r={\zeta }_s\hspace{0.17em}l,\hspace{0.17em}x=z\hspace{0.17em}l+z_r,\hspace{0.17em}{l}_0=l-z_r.$$

Making use of (2), (5), and (6) into (3) yields the following two dimensionless EOM

$$\begin{array}{l}\ddot{z}+c_1\hspace{0.17em}\dot{z}+z+3\hspace{0.17em}\alpha \hspace{0.17em}{\zeta }_{s}z(\hspace{0.17em}{\zeta }_s+z)+\alpha \hspace{0.17em}{z}^3-r\hspace{0.17em}{p}^2(\cos p\tau +\theta \hspace{0.17em}\sin \hspace{0.17em}p\tau )\\ \hspace{0.17em}+\frac{1}{2}{W}^2{\theta }^2+PM-(1+z){\dot{\theta }}^2=f_1\cos p_1\tau \hspace{0.17em},\end{array}$$

(6)

$$\begin{array}{l}{(1+z)}^2\ddot{\theta }+c_2\hspace{0.17em}\dot{\theta }+(1+z)[W^2(\theta -\frac{{\theta }^3}{6})-r\hspace{0.17em}{p}^2(\sin p\tau -\theta \hspace{0.17em}\cos \hspace{0.17em}p\tau )]\\ \hspace{0.17em}+2(1+z)\dot{z}\dot{\theta }=f_2\cos p_2\tau .\end{array}$$

(7)

According to the above procedure, the dimensionless mechanism equations of the piezoelectric circuit and electromagnetic one can be obtained from the substitution of (6) into (4) as follows

$$\begin{array}{l}\hspace{0.17em}\dot{v}+\frac{v}{R_p\hspace{0.17em}{c}_p\hspace{0.17em}{\omega }_1}=\frac{l\hspace{0.17em}{\gamma }_1}{c_p}\dot{z}\\ \dot{q}+\frac{q\hspace{0.17em}{R}_m}{\hspace{0.17em}{l}_m\hspace{0.17em}{\omega }_1}=\frac{l\hspace{0.17em}{\gamma }_2}{l_m}\dot{z}\end{array}$$

3. The Proposed Method

In the present section, the AMS is applied to achieve the asymptotic solutions of the above Equations (7)–(9). Therefore, we examine the dynamics of considered systems in a small neighborhood of their static equilibrium position [41]. Consequently, we express the amplitudes of all oscillations in terms of a small parameter $0<\epsilon <<1$ as follows

$$\begin{array}{l}z(\tau )=\epsilon \hspace{0.17em}\tilde{x}(\tau ;\epsilon ),\hspace{0.17em}\theta (\tau )=\epsilon \hspace{0.17em}\tilde{\phi }(\tau ;\epsilon ),\hspace{0.17em}\\ v(\tau )=\epsilon \hspace{0.17em}\tilde{v}(\tau ;\epsilon )\hspace{0.17em},\hspace{0.17em}q(\tau )=\epsilon \hspace{0.17em}\tilde{q}(\tau ;\epsilon ).\hspace{0.17em}\end{array}$$

(8)

Consider employment of the following variables and parameters

$$\begin{array}{l}r={\epsilon }^2\hspace{0.17em}\tilde{r},\hspace{0.17em}{\mu }_1={\epsilon }^2\hspace{0.17em}{\tilde{\mu }}_1,\hspace{0.17em}{\mu }_2={\epsilon }^2\hspace{0.17em}{\tilde{\mu }}_2,\hspace{0.17em}{l}_m\hspace{0.17em}={\epsilon }^2{\tilde{l}}_m\\ {\gamma }_1={\epsilon }^2\hspace{0.17em}{\tilde{\gamma }}_1,\hspace{0.17em}{\gamma }_2={\epsilon }^2\hspace{0.17em}{\tilde{\gamma }}_2,\hspace{0.17em}\alpha ={\epsilon }^2\hspace{0.17em}\tilde{\alpha },\\ c_1=\hspace{0.17em}{\epsilon }^2\hspace{0.17em}{\tilde{c}}_1,\hspace{0.17em}{f}_1=\hspace{0.17em}{\epsilon }^3\hspace{0.17em}{\tilde{f}}_1\hspace{0.17em},\hspace{0.17em}{f}_2=\hspace{0.17em}{\epsilon }^3\hspace{0.17em}{\tilde{f}}_2,\\ R_p={\tilde{R}}_p/{\epsilon }^2\hspace{0.17em},\hspace{0.17em}{c}_p\hspace{0.17em}={\epsilon }^2{\tilde{c}}_p,\hspace{0.17em}{R}_m={\epsilon }^2{\tilde{R}}_m.\end{array}$$

We looked for the desired solutions $\tilde{x},\hspace{0.17em}\tilde{\phi },\tilde{v},$ and $\tilde{q}$ as power series of $\epsilon$ according to

$$\begin{array}{l}\tilde{x}={\displaystyle {\sum }_{k=1}^3{\epsilon }^k{x}_k({\tau }_0,\hspace{0.17em}{\tau }_1,\hspace{0.17em}{\tau }_2)}+O({\epsilon }^4),\\ \tilde{\phi }={\displaystyle {\sum }_{k=1}^3{\epsilon }^k{\phi }_k({\tau }_0,\hspace{0.17em}{\tau }_1,\hspace{0.17em}{\tau }_2)}+O({\epsilon }^4),\\ \tilde{v}={\displaystyle {\sum }_{k=1}^3{\epsilon }^k{v}_k({\tau }_0,\hspace{0.17em}{\tau }_1,\hspace{0.17em}{\tau }_2)}+O({\epsilon }^4),\\ \tilde{q}={\displaystyle {\sum }_{k=1}^3{\epsilon }^k{q}_k({\tau }_0,\hspace{0.17em}{\tau }_1,\hspace{0.17em}{\tau }_2)}+O({\epsilon }^4),\end{array}$$

(9)

where ${\tau }_n={\epsilon }^n\tau \hspace{0.17em}(n=0,\hspace{0.17em}1\hspace{0.17em},\hspace{0.17em}2)$ are different time scales.

The system of Equations (7)–(9) was transformed to a set of partial differential equations (PDE), i.e., we used the following operators of time derivatives

$$\begin{array}{l}\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(\frac{{\partial }^2}{\partial {\tau }_1^2}+2\frac{{\partial }^2}{\partial {\tau }_0\partial {\tau }_2})+O({\epsilon }^3).\end{array}$$

(10)

Terms of $O({\epsilon }^3)$ and higher orders in the operators (12) are neglected. Employing a procedure of splitting with respect to perturbation parameter $\epsilon$, the following nine linear PDE are obtained

(i) Order of $(\epsilon )$

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

(11)

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

(12)

$$\begin{array}{l}\frac{\partial v_1}{\partial {\tau }_0}+\frac{v_1}{{\omega }_1{\tilde{c}}_p{\tilde{R}}_p}=\frac{l\hspace{0.17em}{\tilde{\gamma }}_1}{{\tilde{c}}_p}\frac{\partial x_1}{\partial {\tau }_0},\\ \frac{\partial q_1}{\partial {\tau }_0}+\frac{q_1\hspace{0.17em}{\tilde{R}}_m}{{\omega }_1{\tilde{l}}_m}=\frac{l\hspace{0.17em}{\tilde{\gamma }}_2}{{\tilde{l}}_m}\frac{\partial x_1}{\partial {\tau }_0}.\end{array}$$

(13)

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

$$\frac{{\partial }^2{x}_2}{\partial {\tau }_0^2}+x_2=\frac{1}{2}{p}^2\hspace{0.17em}\tilde{r}{e}^{i\hspace{0.17em}p\hspace{0.17em}{\tau }_0}-\frac{1}{2}{W}^2{\phi }_1^2+{(\frac{\partial {\phi }_1}{\partial {\tau }_0})}^2-2\frac{{\partial }^2{x}_1}{\partial {\tau }_0\partial {\tau }_1},\hspace{0.17em}$$

(14)

$$\begin{array}{l}\frac{{\partial }^2{\phi }_2}{\partial {\tau }_0^2}+W^2{\phi }_2=\frac{1}{2}i\hspace{0.17em}{p}^2\hspace{0.17em}\tilde{r}{e}^{i\hspace{0.17em}p\hspace{0.17em}{\tau }_0}-W^2{x}_1{\phi }_1-2\frac{\partial x_1}{\partial {\tau }_0}\frac{\partial {\phi }_1}{\partial {\tau }_0}\\ \hspace{0.17em}-2\frac{{\partial }^2{\phi }_1}{\partial {\tau }_0{\tau }_1}\hspace{0.17em}-2\hspace{0.17em}{x}_1\hspace{0.17em}\frac{{\partial }^2{\phi }_1}{\partial {\tau }_0^2},\end{array}$$

(15)

$$\begin{array}{l}\frac{\partial v_2}{\partial {\tau }_0}+\frac{v_2}{{\omega }_1{\tilde{c}}_p{\tilde{R}}_p}=\frac{l\hspace{0.17em}{\tilde{\gamma }}_1}{{\tilde{c}}_p}(\frac{\partial x_1}{\partial {\tau }_1}+\frac{\partial x_2}{\partial {\tau }_0})-\frac{\partial v_1}{\partial {\tau }_1},\\ \frac{\partial q_2}{\partial {\tau }_0}+\frac{q_2\hspace{0.17em}{\tilde{R}}_m}{{\omega }_1{\tilde{l}}_m}=\frac{l\hspace{0.17em}{\tilde{\gamma }}_2}{{\tilde{l}}_m}(\frac{\partial x_1}{\partial {\tau }_1}+\frac{\partial x_2}{\partial {\tau }_0})-\frac{\partial q_1}{\partial {\tau }_1}.\end{array}$$

(16)

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

$$\begin{array}{rr}\frac{{\partial }^2{x}_3}{\partial {\tau }_0^2}+x_3& =\frac{1}{2}{\tilde{f}}_1{e}^{i\hspace{0.17em}{p}_1\hspace{0.17em}{\tau }_0}-\frac{1}{2}i\hspace{0.17em}{p}^2\hspace{0.17em}\tilde{r}{\phi }_1{e}^{i\hspace{0.17em}p\hspace{0.17em}{\tau }_0}-W^2{\phi }_1{\phi }_2-\frac{{\partial }^2{x}_1}{\partial {\tau }_1^2}-{\tilde{c}}_1\frac{\partial x_1}{\partial {\tau }_0}\\ & -2(\frac{{\partial }^2{x}_1}{\partial {\tau }_0\partial {\tau }_2}+\frac{{\partial }^2{x}_2}{\partial {\tau }_0\partial {\tau }_1})+2(\frac{\partial {\phi }_1}{\partial {\tau }_0}\frac{\partial {\phi }_1}{\partial {\tau }_1}+\frac{\partial {\phi }_1}{\partial {\tau }_0}\frac{\partial {\phi }_2}{\partial {\tau }_0})\hfill \\ & +\hspace{0.17em}{x}_1{(\frac{\partial {\phi }_1}{\partial {\tau }_0})}^2-{\tilde{\mu }}_1\hspace{0.17em}{v}_1(or\hspace{0.17em}{\tilde{\mu }}_2\hspace{0.17em}{q}_1\hspace{0.17em})-3\hspace{0.17em}\tilde{\alpha }\hspace{0.17em}{\zeta }_s{}^2{x}_1,\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill \frac{{\partial }^2{\phi }_3}{\partial {\tau }_0^2}+W^2{\phi }_3& =\frac{1}{2}({\tilde{f}}_2{e}^{i\hspace{0.17em}{p}_2\hspace{0.17em}{\tau }_0}-i\hspace{0.17em}{p}^2\hspace{0.17em}\tilde{r}{x}_1{e}^{i\hspace{0.17em}p\hspace{0.17em}{\tau }_0}-p^2\hspace{0.17em}\tilde{r}{\phi }_1{e}^{i\hspace{0.17em}p\hspace{0.17em}{\tau }_0})\hfill \\ & -W^2({\phi }_2{x}_1+{\phi }_1{x}_2)\hspace{0.17em}+\frac{1}{6}{W}^2{\phi }_1^3-\frac{{\partial }^2{\phi }_1}{\partial {\tau }_1^2}-{\tilde{c}}_2\frac{\partial {\phi }_1}{\partial {\tau }_0}\hfill \\ & -2\frac{\partial x_1}{\partial {\tau }_0}(\frac{\partial {\phi }_1}{\partial {\tau }_1}+\frac{\partial {\phi }_2}{\partial {\tau }_0})-2\frac{\partial {\phi }_1}{\partial {\tau }_0}(\frac{\partial x_1}{\partial {\tau }_1}+\frac{\partial x_2}{\partial {\tau }_0})\hfill \\ & -2x_1(\frac{\partial x_1}{\partial {\tau }_0}\frac{\partial {\phi }_1}{\partial {\tau }_0}+2\frac{{\partial }^2{\phi }_1}{\partial {\tau }_0\partial {\tau }_1}+\frac{{\partial }^2{\phi }_2}{\partial {\tau }_0^2})-\frac{{\partial }^2{\phi }_1}{\partial {\tau }_0^2}\hfill \\ & \times (x_1^2+2x_2)-2(\frac{{\partial }^2{\phi }_1}{\partial {\tau }_0\partial {\tau }_2}+\frac{{\partial }^2{\phi }_2}{\partial {\tau }_0\partial {\tau }_1}),\hfill \end{array}$$

(18)

$$\begin{array}{l}\frac{\partial v_3}{\partial {\tau }_0}+\frac{v_3}{{\omega }_1{\tilde{c}}_p{\tilde{R}}_p}=\frac{l\hspace{0.17em}{\tilde{\gamma }}_1}{{\tilde{c}}_p}(\frac{\partial x_1}{\partial {\tau }_2}+\frac{\partial x_2}{\partial {\tau }_1}+\frac{\partial x_3}{\partial {\tau }_0})-\frac{\partial v_1}{\partial {\tau }_2}-\frac{\partial v_2}{\partial {\tau }_1},\\ \frac{\partial q_3}{\partial {\tau }_0}+\frac{q_3\hspace{0.17em}{\tilde{R}}_m}{{\omega }_1{\tilde{l}}_m}=\frac{l\hspace{0.17em}{\tilde{\gamma }}_2}{{\tilde{l}}_m}(\frac{\partial x_1}{\partial {\tau }_2}+\frac{\partial x_2}{\partial {\tau }_1}+\frac{\partial x_3}{\partial {\tau }_0})-\frac{\partial q_1}{\partial {\tau }_2}-\frac{\partial q_2}{\partial {\tau }_1}.\end{array}$$

(19)

These equations can be solved sequentially. Therefore, the solutions of Equations (13)–(15) are as follows

$$x_1\hspace{0.17em}=A\hspace{0.17em}{e}^{i{\tau }_0}+\overline{A}\hspace{0.17em}{e}^{-i{\tau }_0},$$

(20)

$${\varphi }_1=B\hspace{0.17em}{e}^{iW{\tau }_0}+\overline{B}{e}^{-iW{\tau }_0},$$

(21)

$$\begin{array}{l}{v}_1\hspace{0.17em}=\frac{{\tilde{R}}_p{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_1}{{\tilde{R}}_p{\omega }_1{\tilde{c}}_p-i}A\hspace{0.17em}{e}^{i{\tau }_0}+\frac{{\tilde{R}}_p{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_1}{{\tilde{R}}_p{\omega }_1{\tilde{c}}_p+i}\overline{A}\hspace{0.17em}{e}^{-i{\tau }_0},\\ q_1\hspace{0.17em}=\frac{{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_2}{{\omega }_1{\tilde{l}}_m-i\hspace{0.17em}{\tilde{R}}_m}A\hspace{0.17em}{e}^{i{\tau }_0}+\frac{{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_2}{{\omega }_1{\tilde{l}}_m+i\hspace{0.17em}{\tilde{R}}_m}\overline{A}\hspace{0.17em}{e}^{-i{\tau }_0},\end{array}$$

(22)

where $A$ and $B$ represent unknown complex functions of the slow time scales ${\tau }_1$ and ${\tau }_2$, while $\overline{A}$ and $\overline{B}$ denote their counterpoint complex conjugates.

Substitution of the solutions (22)–(24) into the higher-order Equations (16)–(18) yields secular terms. To eliminate these terms, the following conditions are employed

$$\frac{\partial A}{\partial {\tau }_1}=0,\hspace{0.17em}\frac{\partial B}{\partial {\tau }_1}=0.\hspace{0.17em}$$

(23)

Therefore, the following second-order solutions are yielded

$$x_2=\hspace{0.17em}{W}^{2}B\hspace{0.17em}\overline{B}+\frac{3\hspace{0.17em}{W}^2{B}^2{e}^{2iW{\tau }_0}\hspace{0.17em}}{2(4\hspace{0.17em}{W}^2-1)}+\frac{\tilde{r}\hspace{0.17em}{p}^2{e}^{i\hspace{0.17em}p{\tau }_0}\hspace{0.17em}}{2(1-\hspace{0.17em}{p}^2)}+CC,\hspace{0.17em}$$

(24)

$$\begin{array}{l}{\phi }_2=-\frac{W(W+2)A\hspace{0.17em}B\hspace{0.17em}{e}^{i(W+1){\tau }_0}}{(2\hspace{0.17em}W+1)}+\frac{W(W-2)\overline{A}\hspace{0.17em}B\hspace{0.17em}{e}^{i(W-1){\tau }_0}}{(2\hspace{0.17em}W-1)}\\ \hspace{0.17em}+\frac{i\hspace{0.17em}\tilde{r}\hspace{0.17em}{p}^2{e}^{i\hspace{0.17em}p{\tau }_0}\hspace{0.17em}}{2(\hspace{0.17em}{p}^2-W^2)}+CC,\end{array}$$

(25)

$$\begin{array}{l}{v}_2\hspace{0.17em}=-\frac{\hspace{0.17em}{p}^3\tilde{r}{\tilde{R}}_p{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_1\hspace{0.17em}{e}^{i\hspace{0.17em}p{\tau }_0}}{2(p^2-1)({\tilde{R}}_{p}p{\omega }_1{\tilde{c}}_p-i)}+\frac{3\hspace{0.17em}{W}^3{B}^2{\tilde{R}}_p\hspace{0.17em}{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_1{e}^{2iW{\tau }_0}}{(4\hspace{0.17em}{W}^2-1)(2{\tilde{R}}_{p}W{\omega }_1{\tilde{c}}_p-i)}+CC\hspace{0.17em},\\ q_2\hspace{0.17em}=-\frac{\hspace{0.17em}{p}^3\tilde{r}{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_2\hspace{0.17em}{e}^{i\hspace{0.17em}p{\tau }_0}}{2(p^2-1)(p{\omega }_1{\tilde{l}}_m-i\hspace{0.17em}{\tilde{R}}_m)}+\frac{3\hspace{0.17em}{W}^3{B}^2{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_2{e}^{2iW{\tau }_0}}{(4\hspace{0.17em}{W}^2-1)(2W{\omega }_1{\tilde{l}}_m-i\hspace{0.17em}{\tilde{R}}_m)}+CC\hspace{0.17em},\end{array}$$

(26)

where $CC$ denotes the conjugates of the previous terms.

As per the previous method, the elimination of the secular terms required that $A$ and $B$ depended just upon ${\tau }_2$. The solutions for Equations (19) and (20) can be obtained by elimination of the secular terms, which yields the following equations

$$\begin{array}{l}i\hspace{0.17em}{\tilde{c}}_{1}A+2i\frac{\partial A}{\partial {\tau }_2}+3\tilde{\alpha }{\zeta }_s^{2}A+\frac{{\tilde{R}}_p{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_1\hspace{0.17em}{\tilde{\mu }}_{1}A}{({\tilde{R}}_p{\omega }_1{\tilde{c}}_p-i)}+\frac{6\hspace{0.17em}{W}^2(W^2-1)A\hspace{0.17em}B\hspace{0.17em}\overline{B}}{(4\hspace{0.17em}{W}^2-1)}=0.\hspace{0.17em}\\ i\hspace{0.17em}{\tilde{c}}_{1}A+2i\frac{\partial A}{\partial {\tau }_2}+3\tilde{\alpha }{\zeta }_s^{2}A+\frac{{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_2\hspace{0.17em}{\tilde{\mu }}_{2}A}{({\omega }_1{\tilde{l}}_m-i\hspace{0.17em}{\tilde{R}}_m)}+\frac{6\hspace{0.17em}{W}^2(W^2-1)AB\hspace{0.17em}\overline{B}}{(4\hspace{0.17em}{W}^2-1)}=0.\hspace{0.17em}\end{array}$$

(27)

$$iW\hspace{0.17em}{\tilde{c}}_{2}B+2iW\frac{\partial B}{\partial {\tau }_2}+\frac{W^2(W^2-1)}{2(4\hspace{0.17em}{W}^2-1)}B[12\hspace{0.17em}A\overline{A}-(8W^2+1)B\overline{B}]=0.\hspace{0.17em}$$

(28)

Finally, we obtain the third-order solutions as follows

$$\begin{array}{l}{x}_3=\frac{{\tilde{f}}_1{e}^{i\hspace{0.17em}{p}_1{\tau }_0}}{2(1-p_1^2)}-\frac{i\hspace{0.17em}{p}^3\tilde{r}}{2(p^2-W^2)}\{\frac{\hspace{0.17em}(p+2W)B\hspace{0.17em}{e}^{i{\tau }_0(p+W)}}{[1-{(p+W)}^2]}+\frac{(p-2W)\overline{B}\hspace{0.17em}{e}^{i{\tau }_0(p-W)}}{[1-{(p-W)}^2]}\}\\ \hspace{0.17em}-\frac{3\hspace{0.17em}WA}{4}\{\frac{(1+W)B^2{e}^{i{\tau }_0(1+2\hspace{0.17em}W)}}{(2W+1)}-\frac{(W-1)\hspace{0.17em}{\overline{B}}^2{e}^{i{\tau }_0(1-2\hspace{0.17em}W)}}{(2W-1)}\}+CC,\end{array}$$

(29)

$$\begin{array}{cc}\hfill {\phi }_3& =\frac{{\tilde{f}}_2{e}^{i\hspace{0.17em}{p}_2{\tau }_0}}{2(W^2-p_2^2)}-\frac{i\hspace{0.17em}{p}^3\tilde{r}}{2(p^2-W^2)}\{\frac{\hspace{0.17em}(p+2)A\hspace{0.17em}{e}^{i{\tau }_0(p+1)}}{[W^2-{(p+1)}^2]}+\frac{\hspace{0.17em}(p-2)\overline{A}{e}^{i{\tau }_0(p-1)}}{[W^2-{(p-1)}^2]}\}\hfill \\ & \hspace{0.17em}+\frac{pA\hspace{0.17em}\tilde{r}\hspace{0.17em}}{2(p^2-1)}\{\frac{\hspace{0.17em}{e}^{i{\tau }_0(p+W)}}{(p+2W)}(p-1+W)(p+1+W)+\frac{\hspace{0.17em}{e}^{i{\tau }_0(p-W)}}{(p-2W)}\hfill \\ & \times (p-1-W)(p+1-W)\}+\frac{WB\hspace{0.17em}}{4}\{\frac{\hspace{0.17em}(W^2+5W+6)\hspace{0.17em}{A}^2{e}^{i{\tau }_0(W+2)}}{(2W+1)}\hfill \\ & +\frac{(W^2-5W+6)\hspace{0.17em}{\overline{A}}^2{e}^{i{\tau }_0(W-2)}}{(2W-1)}\}-\frac{\hspace{0.17em}(49\hspace{0.17em}{W}^2-1)B^3{e}^{3\hspace{0.17em}i\hspace{0.17em}W{\tau }_0}}{48\hspace{0.17em}(4\hspace{0.17em}{W}^2-1)}\hspace{0.17em}+CC,\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill v_3& =\frac{i\hspace{0.17em}{\omega }_1\hspace{0.17em}{\tilde{R}}_p\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_1\hspace{0.17em}{p}_1{\tilde{f}}_1{e}^{i\hspace{0.17em}{p}_1{\tau }_0}}{2(1+i\hspace{0.17em}{\omega }_1\hspace{0.17em}{\tilde{c}}_p\hspace{0.17em}{\tilde{R}}_p\hspace{0.17em}{p}_1)(1-p_1^2)}-\frac{3\hspace{0.17em}i{\omega }_1\hspace{0.17em}{\tilde{R}}_p\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_1\hspace{0.17em}W\hspace{0.17em}(W+1)A{B}^2\hspace{0.17em}{e}^{i{\tau }_0(2W+1)}}{4[1+i\hspace{0.17em}{\omega }_1\hspace{0.17em}{\tilde{c}}_p\hspace{0.17em}{\tilde{R}}_p\hspace{0.17em}(2W+1)]}\\ & -\frac{i\hspace{0.17em}{\omega }_1\hspace{0.17em}{\tilde{R}}_p\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_1\hspace{0.17em}{p}^3\hspace{0.17em}(p-2W)\hspace{0.17em}\tilde{r}\hspace{0.17em}\overline{B}{e}^{i{\tau }_0(p-W)}}{2[1+i\hspace{0.17em}{\omega }_1\hspace{0.17em}{\tilde{c}}_p\hspace{0.17em}{\tilde{R}}_p\hspace{0.17em}(p-W)](p-W-1)(p-W+1)(p+W)}\hfill \\ & -\frac{i\hspace{0.17em}{\omega }_1\hspace{0.17em}{\tilde{R}}_p\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_1\hspace{0.17em}{p}^3\hspace{0.17em}(p+2W)\hspace{0.17em}\tilde{r}\hspace{0.17em}B\hspace{0.17em}{e}^{i{\tau }_0(p+W)}}{2[1+i\hspace{0.17em}{\omega }_1\hspace{0.17em}{\tilde{c}}_p\hspace{0.17em}{\tilde{R}}_p\hspace{0.17em}(p+W)](p+W-1)(p+W+1)(p-W)}\hfill \\ & -\frac{3\hspace{0.17em}i{\omega }_1\hspace{0.17em}{\tilde{R}}_p\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_1\hspace{0.17em}W\hspace{0.17em}(W-1)A{\overline{B}}^2\hspace{0.17em}{e}^{i{\tau }_0(-2W+1)}}{4[1+i\hspace{0.17em}{\omega }_1\hspace{0.17em}{\tilde{c}}_p\hspace{0.17em}{\tilde{R}}_p\hspace{0.17em}(-2W+1)]}+CC,\hfill \\ \hfill q_3& =\frac{i\hspace{0.17em}{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_2\hspace{0.17em}{p}_1{\tilde{f}}_1{e}^{i\hspace{0.17em}{p}_1{\tau }_0}}{2({\tilde{R}}_m+i\hspace{0.17em}{\omega }_1\hspace{0.17em}{\tilde{l}}_m\hspace{0.17em}{p}_1)(1-p_1^2)}-\frac{3\hspace{0.17em}i{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_2\hspace{0.17em}W\hspace{0.17em}(W+1)A{B}^2\hspace{0.17em}{e}^{i{\tau }_0(2W+1)}}{4[{\tilde{R}}_m+i\hspace{0.17em}{\omega }_1\hspace{0.17em}{\tilde{l}}_m\hspace{0.17em}(2W+1)]}\hfill \\ & -\frac{i\hspace{0.17em}{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_2\hspace{0.17em}{p}^3\hspace{0.17em}(p-2W)\hspace{0.17em}\tilde{r}\hspace{0.17em}\overline{B}{e}^{i{\tau }_0(p-W)}}{2[{\tilde{R}}_m+i\hspace{0.17em}{\omega }_1\hspace{0.17em}{\tilde{l}}_m\hspace{0.17em}(p-W)](p-W-1)(p-W+1)(p+W)}\hfill \\ & -\frac{i\hspace{0.17em}{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_2\hspace{0.17em}{p}^3\hspace{0.17em}(p+2W)\hspace{0.17em}\tilde{r}\hspace{0.17em}B\hspace{0.17em}{e}^{i{\tau }_0(p+W)}}{2[{\tilde{R}}_m+i\hspace{0.17em}{\omega }_1\hspace{0.17em}{\tilde{l}}_m\hspace{0.17em}(p+W)](p+W-1)(p+W+1)(p-W)}\hfill \\ & -\frac{3\hspace{0.17em}i{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_2\hspace{0.17em}W\hspace{0.17em}(W-1)A{\overline{B}}^2\hspace{0.17em}{e}^{i{\tau }_0(-2W+1)}}{4[{\tilde{R}}_m+i\hspace{0.17em}{\omega }_1\hspace{0.17em}{\tilde{l}}_m\hspace{0.17em}(1-2W)]}+CC.\hfill \end{array}$$

(31)

4. Vibrations and Resonance Conditions

This section is introduced to review and classify the resonance cases and to obtain the ME. It is known that these cases arise when the denominator approaches zero [42]. One may consider the following cases:

Primary external resonance, at $p_1\approx 1,\hspace{0.17em}{p}_2\approx W$,

Internal resonance occurs if $W\approx \pm 0.5$.

If one of the above resonance conditions is satisfied, the behavior of the system is extremely complex. It is notable that the achieved solutions are still valid when the oscillations are clear of resonance. Now, we highlight the two cases of primary external resonance if they occur together, i.e., $p_1\approx 1,\hspace{0.17em}{p}_2\approx W.$ For this purpose, we introduce the following detuning parameters ${\sigma }_j\hspace{0.17em}(j=1,\hspace{0.17em}2)$, which express the nearness $p_1$ and $p_2$ to $1$ and $W$

$$\begin{array}{l}{p}_1=1+{\sigma }_1,\\ p_2=W+{\sigma }_2,\end{array}$$

(32)

where ${\sigma }_j={\epsilon }^2{\tilde{\sigma }}_j$.

At that point, the effectiveness of resonance reverberates in secular terms. The parameters of detuning are viewed as a proportion of separation of the vibrations from resonance. The conditions of solvability can be acquired as a result of eliminating secular terms; therefore, the following equations, which estimate the required conditions, are satisfied

$$\begin{array}{l}\frac{1}{2}{\tilde{f}}_1\hspace{0.17em}{e}^{i{\tau }_2{\tilde{\sigma }}_1}-\hspace{0.17em}i\hspace{0.17em}{\tilde{c}}_1\hspace{0.17em}A-2i\frac{\partial A}{\partial {\tau }_2}-3\tilde{\alpha }{\zeta }_s^{2}A-\tilde{H}-\frac{6\hspace{0.17em}{W}^2(W^2-1)A\hspace{0.17em}B\hspace{0.17em}\overline{B}}{(4\hspace{0.17em}{W}^2-1)}=0,\\ \frac{1}{2}{\tilde{f}}_2\hspace{0.17em}{e}^{i{\tau }_2{\tilde{\sigma }}_2}-iW\hspace{0.17em}{\tilde{c}}_2\hspace{0.17em}B-\frac{\hspace{0.17em}{W}^2(W^2-1)\hspace{0.17em}}{2(4\hspace{0.17em}{W}^2-1)}B[12A\overline{A}\hspace{0.17em}-(8W^2+1)B\overline{B}]\\ -2iW\frac{\partial \hspace{0.17em}B}{\partial {\tau }_2}=0.\hspace{0.17em}\end{array}$$

(33)

where

$\tilde{H}=\frac{{\tilde{R}}_p{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_1\hspace{0.17em}{\tilde{\mu }}_1\hspace{0.17em}A}{({\tilde{R}}_p{\omega }_1{\tilde{c}}_p-i)}$ for model (a) or $\tilde{H}=\frac{{\omega }_1\hspace{0.17em}l\hspace{0.17em}{\tilde{\gamma }}_2\hspace{0.17em}{\tilde{\mu }}_2\hspace{0.17em}A}{({\omega }_1{\tilde{l}}_m-i\hspace{0.17em}{\tilde{R}}_m)}$ for model (b).

The conditions of solvability (25) and (35) establish a system consisting of four nonlinear PDE in terms of the unknown functions $A$ and $B$ that depend on the slow time scale ${\tau }_2$. Therefore, we can assume these functions in the polar form as

$$\begin{array}{l}A=\frac{\tilde{a}({\tau }_2)}{2}{e}^{i{\tilde{\psi }}_1{\tau }_2};\hspace{0.17em}a=\epsilon \tilde{a},\hspace{0.17em}\\ B=\frac{\tilde{b}({\tau }_2)}{2}{e}^{i{\tilde{\psi }}_2{\tau }_2};\hspace{0.17em}b=\epsilon \tilde{b}.\end{array}$$

(34)

Since $A$ and $B$ are independent functions of variables ${\tau }_0$ and ${\tau }_1$, then the first-order derivative operator can be simplified as the form

$$\frac{\partial A}{\partial \tau }={\epsilon }^2\frac{\partial A}{\partial {\tau }_2};\hspace{0.17em}\frac{\partial B}{\partial \tau }={\epsilon }^2\frac{\partial B}{\partial {\tau }_2}\hspace{0.17em}.\hspace{0.17em}$$

(35)

Based on these conditions, the conditions of solvability (35) can be transformed into PDE according to the following counterpoint modified phases

$$\begin{array}{l}{\theta }_1({\tau }_1,{\tau }_2)={\tau }_1\hspace{0.17em}{\tilde{\sigma }}_1-{\psi }_1({\tau }_2);\hspace{0.17em}\\ {\theta }_2({\tau }_1,{\tau }_2)={\tau }_1\hspace{0.17em}{\tilde{\sigma }}_2-{\psi }_2({\tau }_2).\end{array}$$

(36)

Substituting (36)–(38) into (35) and separating the real and imaginary portions, the following system is obtained

$$\begin{array}{l}a\frac{d{\theta }_1}{d\tau }=\frac{f_1}{2}\cos {\theta }_1-\frac{a}{2}[\frac{3W^2(W^2-1)\hspace{0.17em}{b}^2\hspace{0.17em}}{2(4W^2-1)\hspace{0.17em}}+3\alpha {\zeta }_s^2+\hspace{0.17em}2(H_1-{\sigma }_1)],\\ \frac{da}{d\tau }=\frac{f_1}{2}\sin {\theta }_1-\frac{a}{2}\hspace{0.17em}[c_1+\hspace{0.17em}2H_2],\\ b\frac{d{\theta }_2}{d\tau }=\frac{f_2}{2W}\cos {\theta }_2+b[{\sigma }_2+\frac{W(W^2-1)}{16(4W^2-1)}((8W^2+1)b^2-12a^2)\hspace{0.17em}],\\ \frac{db}{d\tau }=\frac{f_2}{2W}\sin {\theta }_2-\frac{1}{2}{c}_{2}b.\\ \end{array}$$

(37)

Here, $H_1=\frac{R_p^2{\omega }_1^2\hspace{0.17em}{\mu }_1\hspace{0.17em}l\hspace{0.17em}{\gamma }_1\hspace{0.17em}{c}_p\hspace{0.17em}}{2(R_p^2{\omega }_1^2{c}_p^2+1)\hspace{0.17em}}$ and $H_2=\frac{R_p\hspace{0.17em}{\omega }_1\hspace{0.17em}{\mu }_1\hspace{0.17em}l\hspace{0.17em}{\gamma }_1\hspace{0.17em}}{2(R_p^2{\omega }_1^2{c}_p^2+1)\hspace{0.17em}}$ for the model (a) while $H_1=\frac{{\omega }_1^2\hspace{0.17em}{\mu }_2\hspace{0.17em}l\hspace{0.17em}{\gamma }_2\hspace{0.17em}{l}_m\hspace{0.17em}}{2({\omega }_1^2{l}_m^2+R_m^2)\hspace{0.17em}}$ and $H_2=\frac{R_m\hspace{0.17em}{\omega }_1\hspace{0.17em}{\mu }_2\hspace{0.17em}l\hspace{0.17em}{\gamma }_2\hspace{0.17em}{l}_m\hspace{0.17em}}{2({\omega }_1^2{l}_m^2+R_m^2)\hspace{0.17em}}$ for the model (b).

Equation (39) can be solved numerically to obtain the amplitudes $a,\hspace{0.17em}b$ and the modified phases ${\theta }_j\hspace{0.17em}(j=1,2)$, through which it is impossible to acquire their solutions analytically. Therefore, the following initial conditions are considered

$$a(0)=0.01,\hspace{0.17em}b(0)=0.02\hspace{0.17em},\hspace{0.17em}{\theta }_j(0)=0.$$

The solutions of these equations are presented in Figure 2 and Figure 3 for the two energy-harvesting models. Curves of these figures explore the time histories of the amplitudes and the modified phases, in which its oscillations gradually decrease over time, as shown in Figure 2I, while the oscillations of Figure 3II behave in a stationary manner after a period of time. Model (a) was elaborated for the piezoelectric device, while model (b) was elaborated for the electromagnetic device. These curves take into account the following fixed parameters

$$\begin{array}{l}{f}_1=1.67*{10}^{-6},\hspace{0.17em}{f}_2=2.22*{10}^{-6},\hspace{0.17em}{c}_1=0.0002,\hspace{0.17em}{c}_2=0.0004,\hspace{0.17em}\\ W=0.057,\hspace{0.17em}p=0.003,p_1=0.0016,\hspace{0.17em}{p}_2=0.0016,\hspace{0.17em}r=0.05,\hspace{0.17em}\\ \alpha =0.00007,{\mu }_1={\mu }_2=0.0125,\hspace{0.17em}{\sigma }_1=0.002,\hspace{0.17em}{\sigma }_2=0.002.\end{array}$$

(38)

It is notable that the difference between the two models becomes clear for the amplitude $a$, the modified phase ${\theta }_1$, and the elongation, $x$, as shown in parts (I) of Figure 2, Figure 3 and Figure 4. The reason is due to the presence of energy-harvesting devices, which is predicted from the first two equations of (39) and the approximate solution of $x$. On the other hand, there is no variation regarding the behaviors of the amplitude $b$, the modified phase ${\theta }_2$, or the angle $\theta$, as observed from parts (II) of the same figures, which is expected from the last two equations of (39) besides the asymptotic solution of the angle $\theta$.

The phase portraits of the obtained solutions $x$ and $\theta$ are displayed in parts of Figure 5. It is noted that the included curves of Figure 5 have closed trajectories, which expresses the steady behavior of these solutions.

A comparison between the numerical results of the original governing equations and the asymptotic solutions is represented graphically in Figure 6I for the elongation $x$ and Figure 6II for the angle $\theta$. A closer look at this figure shows high consistency between both solutions, which reflects the good accuracy of the obtained approximate solutions.

5. Steady-State Solutions

In this section, we focus our attention on investigating the solutions at the steady-state case of the considered dynamical models. This case corresponds to the zero values of $\frac{d{\theta }_1}{d\tau },\frac{d{\theta }_2}{d\tau },\hspace{0.17em}\frac{da}{d\tau }$, and $\frac{db}{d\tau }$ [43]. Therefore, we can use Equation (39) to obtain the following form of algebraic equations

$$\begin{array}{l}{f}_1\cos {\theta }_1-\frac{a}{2}\{\frac{3W^2(W^2-1)\hspace{0.17em}{b}^2\hspace{0.17em}}{(4W^2-1)\hspace{0.17em}}+2[3\alpha {\zeta }_s^2+2(\hspace{0.17em}{H}_1-{\sigma }_1)]\}=0,\\ f_1\sin {\theta }_1-a\hspace{0.17em}(c_1+\hspace{0.17em}2H_2)=0,\\ \frac{f_2}{W}\cos {\theta }_2+b\{2{\sigma }_2+\frac{W(W^2-1)}{8(4W^2-1)}[(8W^2+1)b^2-12a^2]\hspace{0.17em}\}=0,\\ \frac{f_2}{W}\sin {\theta }_2-c_{2}b=0.\end{array}$$

(39)

Removing the modified phases ${\theta }_1$ and ${\theta }_2$ from Equation (40) produces the relationships between both amplitudes and the frequency through the introduced detuning parameters

$$\begin{array}{l}{f}_1^2=a^2{\{\frac{3W^2(W^2-1)\hspace{0.17em}{b}^2\hspace{0.17em}}{2(4W^2-1)\hspace{0.17em}}+3\alpha {\zeta }_s^2+2(H_1-{\sigma }_1)]}^2+{[c_1+2\hspace{0.17em}{H}_2]}^2\},\\ f_2^2=W^2{b}^2\{{[2{\sigma }_2+\frac{W(W^2-1)}{8(4W^2-1)}((8W^2+1)b^2-12a^2)\hspace{0.17em}]}^2+c_2^2\}.\end{array}$$

(40)

Investigation of stability considers one of the fatal parts of the vibrations at a steady state. The behavior of the system was explored in the neighborhood region to the fixed points. To fulfil this target, let us assume

$$\begin{array}{l}a=a_{10}+a_{11},\hspace{0.17em}b=b_{10}+b_{11},\\ {\theta }_1={\theta }_{10}+{\theta }_{11},\hspace{0.17em}{\theta }_2={\theta }_{20}+{\theta }_{21},\end{array}$$

(41)

where $a_{10},{\theta }_{10},b_{10},$ and ${\theta }_{20}$ represent the steady-state solutions, while $a_{11},{\theta }_{11},b_{11},$ and ${\theta }_{21}$ denote the tiny perturbations comparing with $a_{10},{\theta }_{10},b_{10},$ and ${\theta }_{20}$.

Substituting (42) into (39), we can obtain the following system of ordinary differential equations after linearization

$$\begin{gathered} \begin{array}{l}{a}_{10}\frac{d{\theta }_{11}}{d\tau }=-\frac{f_1}{2}\sin {\theta }_{10}{\theta }_{11}-\frac{a_{11}}{2}[\frac{3W^2(W^2-1)\hspace{0.17em}{b}_{10}^2\hspace{0.17em}}{2(4W^2-1)\hspace{0.17em}}+3\alpha {\zeta }_s^2+2\hspace{0.17em}(H_1-{\sigma }_1)]\\ \\ \end{array} \\ \begin{array}{l}\hspace{0.17em}-\frac{3W^2(W^2-1)\hspace{0.17em}{a}_{10}{b}_{10}\hspace{0.17em}{b}_{11}\hspace{0.17em}}{2(4W^2-1)\hspace{0.17em}},\\ \end{array} \\ \begin{array}{l}\frac{d{a}_{11}}{d\tau }=\frac{f_1}{2}\cos {\theta }_{10}{\theta }_{11}-\frac{a_{11}}{2}\hspace{0.17em}[c_1+\hspace{0.17em}2H_2],\\ \end{array} \\ \begin{array}{l}{b}_{10}\frac{d{\theta }_{21}}{d\tau }=-\frac{f_2}{2W}\sin {\theta }_{20}{\theta }_{21}+b_{11}[{\sigma }_2+\frac{3W(W^2-1)}{16(4W^2-1)}((8W^2+1)b_{10}^2-4a_{10}^2)\hspace{0.17em}]\\ \end{array} \\ \begin{array}{l}\hspace{0.17em}-\frac{3W(W^2-1)\hspace{0.17em}{a}_{10}{b}_{10}\hspace{0.17em}{a}_{11}\hspace{0.17em}}{2(4W^2-1)\hspace{0.17em}},\\ \end{array} \\ \frac{d{b}_{11}}{d\tau }=\frac{f_2}{2W}\cos {\theta }_{20}{\theta }_{21}-\frac{1}{2}{c}_2{b}_{11}. \end{gathered}$$

(42)

Considering that $a_{11},{\theta }_{11},b_{11},$ and ${\theta }_{21}$ are the unknown perturbed functions of the linear system (39), we can express their solutions in the form $\hspace{0.17em}{k}_s\hspace{0.17em}{e}^{\lambda \tau }$ in which $k_s\hspace{0.17em}(s=1,2,3,4)$ are constants, and $\lambda$ is the eigenvalue of the unknown perturbation. Consequently, the fixed points in (43) are asymptotically stable when the real parts of the roots of the following characteristic equation of (43),

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

(43)

are negative, Here ${\mathsf{\Gamma }}_1,{\mathsf{\Gamma }}_2,{\mathsf{\Gamma }}_3,$ and ${\mathsf{\Gamma }}_4$ take the forms

$$\begin{array}{l}{\mathsf{\Gamma }}_1=\frac{1}{2}(c_1+c_2+H_2+\frac{f_1}{a_{10}}\sin {\theta }_{10}+\frac{f_2}{b_{10}\hspace{0.17em}W}\sin {\theta }_{20}),\\ {\mathsf{\Gamma }}_2=\frac{c_2}{4}(c_1+H_2)+\frac{f_1\hspace{0.17em}{f}_2\hspace{0.17em}\sin {\theta }_{10}\sin {\theta }_{20}}{4a_{10}{b}_{10}W}+\frac{f_1\hspace{0.17em}\cos {\theta }_{10}}{4a_{10}}[\frac{3W^2(W^2-1)\hspace{0.17em}\hspace{0.17em}{b}_{10}^2\hspace{0.17em}}{4(4W^2-1)\hspace{0.17em}}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+6\alpha {\zeta }_s^2-4({\sigma }_1-H_1)]+\frac{f_2\hspace{0.17em}\cos {\theta }_{20}}{64b_{10}W}[\frac{3W^2(W^2-1)(-8a_{10}^2+b_{10}^2(8W^2+1))}{(4W^2-1)}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+32{\sigma }_2]+\frac{f_2\hspace{0.17em}\sin {\theta }_{10}}{4a_{10}}[c_1+c_2+H_2]+\frac{f_2\hspace{0.17em}{H}_2\sin {\theta }_{20}}{4b_{10}W}+\frac{f_2\hspace{0.17em}(c_1+c_2)\sin {\theta }_{20}}{4b_{10}W},\end{array}$$

$$\begin{array}{l}{\mathsf{\Gamma }}_3=\frac{c_2{f}_1\hspace{0.17em}\cos {\theta }_{10}}{16a_{10}}[\frac{3W^2(W^2-1)\hspace{0.17em}\hspace{0.17em}{b}_{10}^2\hspace{0.17em}}{(4W^2-1)\hspace{0.17em}}+6\alpha {\zeta }_s^2-4({\sigma }_1-H_1)]+\frac{c_1\hspace{0.17em}{f}_2\hspace{0.17em}\cos {\theta }_{20}}{4b_{10}W}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+(c_1+H_2)[\frac{c_1{f}_2\hspace{0.17em}\cos {\theta }_{20}}{128b_{10}}[\frac{3(W^2-1)(-8a_{10}^2+b_{10}^2(8W^2+1))}{(4W^2-1)}]+\frac{c_2{f}_1\hspace{0.17em}\sin {\theta }_{10}}{8a_{10}}]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{f_1\hspace{0.17em}{f}_2\hspace{0.17em}\cos {\theta }_{20}\sin {\theta }_{10}}{128a_{10}{b}_{10}W}[\frac{3W(W^2-1)(-8a_{10}^2+b_{10}^2(8W^2+1))}{(4W^2-1)}+32{\sigma }_2]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{f_1{f}_2\hspace{0.17em}\cos {\theta }_{10}\sin {\theta }_{20}}{16a_{10}{b}_{10}W}[\frac{3W^2(W^2-1)\hspace{0.17em}\hspace{0.17em}{b}_{10}^2\hspace{0.17em}}{(4W^2-1)\hspace{0.17em}}+6\alpha {\zeta }_s^2-4({\sigma }_1-H_1)]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{c_2{f}_2\hspace{0.17em}\sin {\theta }_{20}}{8b_{10}W}(c_1+H_2)+\frac{f_1{f}_2\hspace{0.17em}{H}_2\sin {\theta }_{10}\sin {\theta }_{20}}{8a_{10}{b}_{10}W}(c_1+c_2+H_2)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{f_2\hspace{0.17em}{H}_2\cos {\theta }_{20}}{4b_{10}W},\end{array}$$

$$\begin{array}{l}{\mathsf{\Gamma }}_4=\frac{3f_1\hspace{0.17em}{f}_2\hspace{0.17em}(W^2-1)\cos {\theta }_{10}\cos {\theta }_{20}}{16a_{10}{b}_{10}(4W^2-1)}[\frac{3}{4}(W^2-1)[3W^2{a}_{10}^2{b}_{10}^2+\frac{1}{4}(-8a_{10}^2+b_{10}^2(8W^2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+1))\hspace{0.17em}\alpha {\zeta }_s^2]+{\sigma }_1[\frac{1}{8}{b}_{10}^2(8W^2+1)+a_{10}^2]]+\frac{f_1\hspace{0.17em}{f}_2\hspace{0.17em}{\sigma }_2\cos {\theta }_{10}\cos {\theta }_{20}}{8a_{10}{b}_{10}W}(3\alpha {\zeta }_s^2-2{\sigma }_1)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{3f_1\hspace{0.17em}{f}_2\hspace{0.17em}{b}_{10}W(W^2-1)\cos {\theta }_{10}\cos {\theta }_{20}}{16a_{10}(4W^2-1)}[{\sigma }_2+\frac{3W(W^2-1)(8W^2+1)\hspace{0.17em}\hspace{0.17em}{b}_{10}^2\hspace{0.17em}}{32(4W^2-1)\hspace{0.17em}}]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{f_1\hspace{0.17em}{f}_2\hspace{0.17em}{H}_1\cos {\theta }_{10}\cos {\theta }_{20}}{128a_{10}{b}_{10}W}[\frac{3W(W^2-1)(-8a_{10}^2+b_{10}^2(8W^2+1))}{(4W^2-1)}+32{\sigma }_2]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{f_1\hspace{0.17em}{f}_2\hspace{0.17em}\cos {\theta }_{20}\sin {\theta }_{10}}{256a_{10}{b}_{10}W}[32{\sigma }_2(c_1+H_2)+\frac{3H_{2}W(W^2-1)(-8a_{10}^2+b_{10}^2(8W^2+1))}{(4W^2-1)}]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{c_2{f}_1\hspace{0.17em}{f}_2\hspace{0.17em}\cos {\theta }_{20}\sin {\theta }_{10}}{32a_{10}{b}_{10}W}[\frac{3b_{10}^{2}W(W^2-1)}{(4W^2-1)}+2(3\alpha {\zeta }_s^2-2({\sigma }_1-H_1))]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{c_2{f}_1\hspace{0.17em}{f}_2\hspace{0.17em}\sin {\theta }_{10}\sin {\theta }_{20}}{16a_{10}{b}_{10}W}(c_1+H_2).\end{array}$$

The appropriate and necessary conditions of the stability for the fixed points, as indicated by Routh–Hurwitz criteria [44], are

$$\begin{array}{l}{\mathsf{\Gamma }}_1>0,\hspace{0.17em}{\mathsf{\Gamma }}_3({\mathsf{\Gamma }}_1\hspace{0.17em}{\mathsf{\Gamma }}_2-{\mathsf{\Gamma }}_3)-{\mathsf{\Gamma }}_4\hspace{0.17em}{\mathsf{\Gamma }}_1^2>0,\\ {\mathsf{\Gamma }}_1\hspace{0.17em}{\mathsf{\Gamma }}_2-{\mathsf{\Gamma }}_3>0,\hspace{0.17em}{\mathsf{\Gamma }}_4>0.\end{array}$$

(44)

It is worth mentioning that all possible steady-state solutions can be described on the plane of coordinates $a,b$ as shown in Figure 7. The solutions of the system of Equation (41) are represented graphically to have the curves of these figures at selected data of ${\sigma }_j\hspace{0.17em}(j=1,2)$. These curves indicate the locus of the roots of these equations, where the intersection points of both curves represent the fixed points that correspond to the steady-state solutions. If these points are marked by black dots, then it expresses a stable steady-state case; otherwise, it is an unstable steady state. The number of possible intersection points depends on the parameters of the dynamical models. The minimum number of possible amplitudes is one, as in Figure 7I, while the obtained maximum number equals five according to the used data of ${\sigma }_j\hspace{0.17em}(j=1,2)$ as explored in Figure 7V.

6. Stability Analysis

The investigated models were analyzed using the nonlinear stability approach [45]. Based on the above, the dynamical motion was considered under the influence of external harmonic force $F(t)$ and moment $M(t)$. The criteria of stability were implemented besides the simulations of the nonlinear system of Equation (39). It is observed that various parameters in the stability criterion, such as the frequency $W$, coefficients of damping $c_j\hspace{0.17em}(j=1,2)$, and the detuning parameters ${\sigma }_j$, play a crucial destabilizing impact. Appropriate policies with varied system (39) parameters were considered to plot the corresponding diagrams of the system’s stability.

For distinct parametric regions, time histories of the amplitudes $a$ and $b$ are displayed, and their attributes are presented through the trajectories of the phase plane. The probably fixed points were drawn. Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18, Figure 19, Figure 20, Figure 21, Figure 22 and Figure 23 explore the variation of these points with the detuning parameters ${\sigma }_j$. In other words, Figure 8 is calculated when ${\sigma }_2=0,W=0.448,$ and $c_2=0.011$ at various values of $c_1$, in which the stable fixed points are represented by solid curves while the unstable ones are drawn by the dashed curves. The stability and instability areas lie in the ranges $-0.05\le {\sigma }_1\hspace{0.17em}<0.001$ and $0.001\le {\sigma }_1\hspace{0.17em}\le 0.05$, respectively. Here, one critical fixed point is observed at the points $(0.0014,0.0157),$$(0.0014,0.0236),$$(0.0014,0.0414)$ and $(0.00149,0.0142),$$(0.00149,0.0141),$ $(0.00149,0.013)$ for Figure 8I,II, respectively. Moreover, Figure 9 and Figure 10 illustrate the frequency response of both $a({\sigma }_1)$ and $b({\sigma }_1)$ at ${\sigma }_2=0.001$ and ${\sigma }_2=-0.001$, respectively. The comparison between Figure 8, Figure 9 and Figure 10 shows the slight variation of the critical fixed points as seen in parts (I) of these figures, while the change of the locus of the corresponding fixed points becomes very clear as observed in parts (II) of the same figures.

Parts of Figure 11 illustrate the variation of $c_1$ values on the behavior of the amplitudes $a$ and $b$ at ${\sigma }_1=0,W=0.448,$ and $c_2=0.011$. It is shown that the stability and the instability areas are $-0.002<{\sigma }_2\hspace{0.17em}\le 0.05$ and $-0.05\le {\sigma }_2\hspace{0.17em}\le -0.002$, respectively, with one critical stable fixed point. Figure 11 is replotted when ${\sigma }_1=-0.001$ to reveal the variation of ${\sigma }_1$ values on the frequency response of each $a({\sigma }_2)$ and $b({\sigma }_2)$ with the same other data as shown in Figure 12. On the other hand, according to the used data, Routh–Hurwitz conditions (45) are not satisfied when ${\sigma }_1=0.001$. When we compared Figure 11 with Figure 12, we concluded that the range of $a$ and $b$ axes varied toward decreasing while the number of critical fixed points remained unchanged.

By changing the damping coefficient, Figure 13 and Figure 14 show at and, respectively. The influence on the behavior and of the frequency response curves becomes clearer without sharp peaks than on the behavior of and, which included sharp peaks. The reason goes back to the system of Equation (39). Only one critical fixed point is observed regardless of the values of and on the whole used domain. The stability and instability $c_2$ areas in this domain of the curves of Figure 13 are $-0.05\le {\sigma }_1\hspace{0.17em}<0.001$ and $0.001\le {\sigma }_1\hspace{0.17em}\le 0.05$, respectively. The corresponding regions of Figure 14 are $-0.001<{\sigma }_2\hspace{0.17em}\le 0.05$ and $-0.05\le {\sigma }_2\hspace{0.17em}<-0.001$.

A closer look at Figure 15 and Figure 16 show that they are calculated at $W=0.448$, $c_1=0.001$ when ${\sigma }_2=0.001$ and ${\sigma }_2=-0.001$, respectively. The comparison between these figures with Figure 13 shows that the corresponding parts of Figure 15 and Figure 16 have critical fixed points more than parts of Figure 13. Therefore, the change of the values ${\sigma }_2$ has a significant impact on the stability areas and on the number of critical points. On the other hand, Figure 14 is regraphed when ${\sigma }_1=-0.001$ with the same other values of Figure 14 to produce Figure 17. It clearly shows that the change ${\sigma }_1$ has a good impact on the behavior of $a$ compared to $b$ as shown in parts (I) and (II) of Figure 17, in which the locus of critical points is changed as in part (I). The positive value ${\sigma }_1=0.001$ produces unfulfilled conditions of Routh–Hurwitz (45).

The influence of different values of frequency $W$ on the frequency response curves is shown in Figure 18 and Figure 19 when $W=0.1$, $c_1=0.001$, $c_2=0.011$ at ${\sigma }_2=0.001$ at and ${\sigma }_1=-0.001$, respectively. The stable fixed points exist in the range $-0.05\le {\sigma }_1\hspace{0.17em}\le 0.001$, while the unstable ones occur in the range $0.001<{\sigma }_1\hspace{0.17em}\le 0.05$, as indicated in parts of Figure 18. It is noted that this figure contains only one sharp peak for each curve directed upward. Parts of Figure 19 explore the stability and instability regions of $a({\sigma }_2)$ and $b({\sigma }_2)$, which lie in the ranges $0.001<{\sigma }_2\hspace{0.17em}\le 0.05$ and $-0.05\le {\sigma }_2\hspace{0.17em}\le 0.001$, respectively. The entire domain contains one critical fixed point for each curve that is directed downward for $a({\sigma }_2)$ as seen in Figure 19I, while it takes an upward direction for $b({\sigma }_2)$ as shown in Figure 19II.

Figure 18 was redrawn at ${\sigma }_2=0$ and ${\sigma }_2=-0.001$ when $W=0.1$, $c_1=0.001$,$c_2=0.011$ to obtain Figure 20 and Figure 21, respectively. In contrast, Figure 19 was replotted again at ${\sigma }_1=0$ to yield Figure 22 when $W=0.1$, $c_1=0.001$ and $c_2=0.011$. The impact of ${\sigma }_2$ values becomes clear on the modified amplitude $b$ than the modified amplitude $a$ as seen in parts of Figure 20 and Figure 21, while the variation with respect to the amplitude $a$ becomes more evident than the change of $b$ as graphed in parts of Figure 22.

Figure 23 shows the frequency response curves of $a({\sigma }_1)$ and $b({\sigma }_1)$ at ${\sigma }_2=0.001,W=0.448,$ $c_1=0.001,$ and ${\sigma }_2(=0.001,\hspace{0.17em}0,\hspace{0.17em}-0.001)$ when. Solid and dashed curves describe the ranges of stable and unstable fixed points, respectively.

We demonstrated the properties of the nonlinear analysis for the nonlinear amplitudes of the equations of the system (39) and explored their stabilities. Therefore, we considered the following transformation

$$\begin{array}{l}A=({\tilde{u}}_1+i\hspace{0.17em}{\tilde{v}}_1)e^{i\hspace{0.17em}{\tilde{\sigma }}_1\hspace{0.17em}{\tau }_2},\\ B=({\tilde{u}}_2+i\hspace{0.17em}{\tilde{v}}_2)e^{i\hspace{0.17em}{\tilde{\sigma }}_2\hspace{0.17em}{\tau }_2},\hspace{0.17em}\end{array}$$

(45)

where $u_j=\epsilon \hspace{0.17em}{\tilde{u}}_j\hspace{0.17em},\hspace{0.17em}{v}_j=\epsilon \hspace{0.17em}{\tilde{v}}_j\hspace{0.17em};\hspace{0.17em}(j=1,2).$ Making use of (39) and (46), the real and the imaginary parts were separated to obtain

$$\begin{gathered} \begin{array}{l}\frac{d{u}_1}{d\tau }=-\frac{3W^2(W^2-1)\hspace{0.17em}{v}_1\hspace{0.17em}}{4(4W^2-1)\hspace{0.17em}}(u_2^2+v_2^2)+v_1({\sigma }_1-\frac{3}{2}\alpha {\zeta }_s^2-\hspace{0.17em}{H}_1-H_2)\\ \end{array} \\ \begin{array}{l}\hspace{0.17em}-\frac{1}{2}{c}_1\hspace{0.17em}{u}_1,\\ \end{array} \\ \begin{array}{l}\frac{d{v}_1}{d\tau }=\frac{3W^2(W^2-1)\hspace{0.17em}{u}_1\hspace{0.17em}}{4(4W^2-1)\hspace{0.17em}}(u_2^2+v_2^2)+\frac{u_1}{2}\hspace{0.17em}[3\alpha {\zeta }_s^2+2\hspace{0.17em}(H_1+H_2-{\sigma }_1)]\\ \end{array} \\ \begin{array}{l}\hspace{0.17em}-\frac{1}{4}(2c_1\hspace{0.17em}{v}_1+f_1),\\ \end{array} \\ \begin{array}{l}\frac{d{u}_2}{d\tau }=\frac{W(W^2-1)\hspace{0.17em}{v}_2\hspace{0.17em}}{4(4W^2-1)\hspace{0.17em}}[(8W^2+1)(u_2^2+v_2^2)-4(u_1^2+v_1^2)]\\ \end{array} \\ \begin{array}{l}\hspace{0.17em}+\frac{1}{2}(2{\sigma }_2{v}_2-c_2\hspace{0.17em}{u}_2),\\ \end{array} \\ \begin{array}{l}\frac{d{v}_1}{d\tau }=-\frac{W(W^2-1)\hspace{0.17em}{u}_2\hspace{0.17em}}{4(4W^2-1)\hspace{0.17em}}[(8W^2+1)(u_2^2+v_2^2)-4(u_1^2+v_1^2)]\\ \end{array} \\ \hspace{0.17em}-\frac{1}{4W}[2W(2{\sigma }_2{u}_2+c_2{v}_2)+f_2], \end{gathered}$$

(46)

The amplitudes are confirmed through the examined interval of time in distinct parametric areas, and the characteristics of the amplitudes are presented in the phase plane curves as seen in Figure 24, Figure 25 and Figure 26 for the dynamical system, which is connected to the piezoelectric device in the model (a) and Figure 27, Figure 28 and Figure 29 for model (b) when the system is attached to the electromagnetic device. The following data are used

$$\begin{array}{l}{f}_1=1.67\times {10}^{-6},\hspace{0.17em}{f}_2=2.22\times {10}^{-6},\hspace{0.17em}{c}_1=0.0002,\hspace{0.17em}{c}_2=0.0004,\hspace{0.17em}\\ W=0.057,\hspace{0.17em}{\mu }_1={\mu }_2==0.0125,\hspace{0.17em}{\zeta }_s=0.003,\hspace{0.17em}\alpha =7\times {10}^{-4},\hspace{0.17em}\\ {\gamma }_1={\gamma }_2=10,\hspace{0.17em}{\sigma }_1=0.002,\hspace{0.17em}{\sigma }_2=0.002.\end{array}$$

Figure 24 and Figure 25 express the variation of new modified phases $u_1,v_1$ and $u_2,v_2$ that is observed in the system of Equation (45) via time $\tau$. Decay curves are sketched over time until the end of the investigated time interval, as seen in parts of Figure 24, while the plotted curves in Figure 25 behave in a steady manner as of the second half of the interval of time. On the other hand, parts of Figure 26 represent the projections of the equations of modulation trajectories on the planes $u_1{v}_1$ and $u_2{v}_2$ in which circular spiral curves are drawn. These curves approach one point, which indicates the stationary behavior of the investigated models. These simulations can be applied to the plotted curves in Figure 27, Figure 28 and Figure 29. The difference between the curves included in these figures and the corresponding ones in Figure 24, Figure 25 and Figure 26 is very slight because it lacks the considerable influence of energy-harvesting devices connected with the dynamical systems on its stability.

The trajectory behavior is an intrinsic indicator of a spiral through the other paths. If these paths appear to be circular, the spiral is stiff. Otherwise, spiral spindles produce a Hopf bifurcation in dramatic systems. To modify the spiral tip, one can use any external constraint or other separated approaches.

7. Energy-Harvesting Device Performance

Based on the mentioned dynamical models, when piezoelectric and electromagnetic devices are connected with the 2DOF spring pendulum systems, one can expect to obtain an energy harvester with enhanced power output.

EH devices can employ vibration energy and produce electricity, as we mentioned before. In our work, we used the considered dynamical system to obtain the necessary mechanical vibrations to power the two different models (a) and (b). As shown in model (a), the piezoelectric device is a transducer, which can convert mechanical stress into electricity using the advantages of piezoelectric materials. The electromagnetic device mechanism in the model (b) is made up of a magnet that is wrapped in a coil to create a magnetic flux. When a dynamical system vibrates, it cuts the magnetic flow, resulting in the generation of electricity.

Moreover, we examined the impact of different selected values of the physical parameters on the system to obtain high performance. The time histories of the output voltages and output powers of the two models of piezoelectric and electromagnetic devices are plotted in Figure 30 and Figure 31. It is noticed that the maximum value of voltage, current, and even the output power decreases with the increase of the damping coefficients $(c=c_1=c_2)$. The increasing excitation amplitude increases the output voltage, current, and power, as indicated in Figure 32 and Figure 33. The piezoelectric device’s output power is much greater than the electromagnetic one of the model (b), as shown in Figure 31 and Figure 33.

The response of the excitation frequencies is represented in the curves of Figure 34 and Figure 35 for different selected values of the damping coefficients. High voltage, current, and power are generated at a frequency of approximately 10 Hz. It is observed that when the value of the damping coefficient increases, the generated voltage and power of the piezoelectric device decrease, as seen in parts (a) of Figure 34 and Figure 35. The current and power of the electromagnetic device remain unchanged due to the current’s approximate equations, as seen in parts (b) of Figure 34 and Figure 35.

8. Conclusions

1.

A 2DOF nonlinear damped vibrating spring pendulum system moving in a circular path with a stationary angular velocity, in which it is connected with energy-harvesting devices, was investigated as a novel model.

2.

The governing equations were derived using Lagrange’s equation and solved asymptotically using the AMS to obtain results with high accuracy.

3.

The solvability conditions were acquired in light of the resonance cases.

4.

The achieved ME was solved to obtain two nonlinear algebraic equations in terms of the amplitudes and the detuning parameters.

5.

The time histories of the dynamical motion, the resonance curves responses, and the solutions at the steady-state cases were graphed to reveal the excellent impact of the selected values of the model parameters on the motion.

6.

The achieved asymptotic solutions were verified compared with the numerical results that reveal the high solidity between them.

7.

The damping coefficients and excitation amplitudes that influence the output voltage, current, and power were checked. Moreover, the output power of the systems was reviewed according to the response of excitation frequency.

8.

The nonlinear stability of the ME was examined through the stability and instability areas for the frequency response curves. In addition, the nonlinear analysis for the nonlinear amplitudes of these equations is presented.

9.

Electrical energy was generated from piezoelectric and electromagnetic devices.