On the influence of an energy harvesting device on a dynamical system

Abstract

This article studies the motion of a three degrees-of-freedom (DOF) dynamical system, which consists of two parts, a dynamic part, and an electromagnetic harvesting device. The composition of the system depends on the magnet oscillating in the coil of this device. An electric potential is exerted because of the electromotive force to obtain energy harvesting (EH). Lagrange's equations (LE) are used to derive the system’s equations of motion (EOM), which are solved analytically up to a third-order of approximation using the multiple-scales method (MSM). These solutions are considered novel, in which the mentioned method is applied to a novel dynamical system. The gained are confirmed through comparison with the numerical solutions of the EOM to reveal a high degree of consistency. To establish the system’s solvable conditions and modulation equations (ME), three scenarios of primary external resonance are investigated simultaneously. Certain diagrams of time histories for the amplitudes and phases of the examined dynamical system are plotted to show their behaviors at any given instance of time. The Routh–Hurwitz criteria (RHC) are used to establish the stability and instability regions in light of the graphs of resonance response curves. The presence of an energy harvester in dynamic devices helps us to convert vibrational motion into electrical energy that can be exploited in several applications in a variety of daily tasks, including structural and environmental monitoring, military applications, space exploration, and remote medical sensing.

Introduction

It is known that non-renewable energy is running out rapidly and will be depleted over the coming decades. Therefore, it is very necessary to change the type of energy production and create renewable, clean energy sources. We can also use so-called EH, in which it is difficult to use large batteries to replace them regularly in portable and wireless systems, so we resort to compact and efficient power sources that enable us to increase the use of wireless devices. Therefore, the best possible idea for using energy in small systems as the main source is EH. Thus, using available energy resources, temporary power tanks can constantly be replaced with charged energy storage elements, including sunlight, heat, and vibration.^1–3 Vibration is particularly charming because of its many uses. Vibration sources are available in environments such as vehicles and human movements.

Energy is harvested from electromagnetic, piezoelectric, or electrostatic transmission. Also, these vibrations often occur at low frequencies, which means that the energy levels are low at the output. The amount of energy produced by such systems is in the micro to milliwatts order.^4–10 Electromagnetic energy originates from Faraday's law of induction ¹¹ since any change in the magnetic environment of a wire coil induces an electromotive force in the coil. This occurs according to the motion of the magnet relative to the coil. ¹² Capacitors are used in electrostatic transformation to convert mechanical energy into an electric field, which may be employed to provide current to an electric circuit. ¹³ When a piezoelectric material is deformed under stress, an electric field is generated, and a voltage is produced. ¹⁴ In Reference 15, the authors demonstrated the study of vibration reduction and EH in a dynamic system of a spring pendulum. In Reference 16, the spring pendulum system was connected to two different EH devices, each in a separate case, in which the controlling EOM were obtained. Moreover, the ME were determined in light of the classification of the examined resonance cases. It is generally acknowledged that mechanical systems, especially vibrational dynamics problems such as References 17–20 are regarded as the most significant problems in mechanics due to their broad spectrum of applications, particularly in engineering. In Reference 17, a general model composed of a nonlinear, damped, excited spring pendulum is studied. Its suspension point was restricted to moving in an elliptic trajectory with a stationary angular velocity, while its other end was connected to a rigid body. For both the numerical and approximate solutions, the time-varying histories of the resulting solutions are displayed. Whereas in Reference 18, a 3DOF nonlinear system was investigated, in which LE were used to derive the governing EOM and solved analytically by applying the MSM method. ¹⁹ Curves of frequency response were graphed to show the influence of the system’s parameters on the stability and instability areas. The chaotic motion of a spring pendulum under the action of a harmonically excited force was investigated in Reference 20. The authors considered that the planar motion of its pivot point was in a circular path. The positive effect of a viscous fluid on a pendulum spring was studied in Reference 21 when the suspension point was stationary, where buoyancy, drag, lift, and external forces were considered.

Systems for EH in dampening vibrations have been examined in References 22,23. It is an electromagnetic harvester attached to the automatic parametric pendulum that acts as the absorber. The overall performance of a combined damping vibration was studied, and it was shown that this combination did not reduce the vibration-damping efficiency. It was observed in Reference 22 that energy could be harvested for the rotation and chaotic motion of the pendulum. In Reference 23, a pendulum system with two independent electromagnetic harvesters is presented. A direct-current motor is added to the swing shaft to serve as a second harvester. The results showed that the proposed system significantly improves the capacity of EH, and the efficiency of vibration suppression for rotary harvesters is low. In Reference 24, vibration absorbers were developed as energy-harvesting devices and vibration-damping systems to suppress undesired vibrations and generate induced voltages that can be used in low-voltage applications such as wireless sensors. In Reference 25, a numerical analysis of a semi-active suspended vibration absorber, or harvesting system, was performed. The primary objective was used to dampen the vibrations of the primary structure, and the secondary objective was considered an absorber that simultaneously harvested energy from the vibrations. The pendulum-tuned mass absorber was regulated by a magnetic levitation collector. The influence of EH, and the control of vibrations by specially designed pendulum dampers, and the use of nonlinear dampers as energy harvesters were studied in Reference 26. The use of energy recovery devices has been shown to provide better vibration reduction in the low-frequency range. In Reference 27, the authors present a new pendulum absorber harvester structure that improves vibration suppression and EH characteristics. The damper configuration of the system refers to the use of magnets, whose performance was studied empirically and theoretically. The results showed that the new damper harvester has greater performance than the traditional device with both oscillation suppression and EH. The movement of a 3DOF dynamic system, which creates an auto-parametric model, is examined in Reference 28. In addition to evaluating stable and unstable zones, the authors used the nonlinear stability approach to investigate the system’s frequency response curves. In Reference 29 the authors analyzed a piezoelectric energy harvester under synchronous resonances. The obtained results showed that the generated voltage and power level are amplified under the simultaneous effects of two strong excitations such as super-harmonic and harmonic resonances. Energy harvesting and Wideband vibration isolation based on a coupled piezoelectric-electromagnetic structure were investigated in Reference 30. Using energy harvesting to power marine environment monitoring was studied in Reference 31. In Reference 32, the vibrational motion of a 3DOF dynamical system connected to an electromagnetic EH device is studied, which converts the vibrational motion to electric energy. The motion of a nonlinear damped dynamical system with a piezoelectric transducer is investigated in Reference 33. The bifurcation diagram, Poincaré diagram, and Lyapunov exponent spectrum are presented. Moreover, three different motion modes of the system, which are quasi-periodic motion, hyperchaotic motion, and chaotic motion, are examined. It is worth noting that high output voltage and power are obtained by controlling the values of natural frequency, excitation amplitude, capacitance, coupling coefficient, load resistance, and damping coefficient.

In this paper, we demonstrate a simple device based on spring vibrations to harvest energy from a vibrating environment, which converts the vibrating mechanical energy into an electrical one. The studied system is composed of two connected masses with two damped springs. One of these masses moves horizontally, and the second one has a planar motion in view of the existing EH device. This system is acted upon by two external forces and one moment. The MSM is utilized to obtain novel approximate solutions that fit the generalized coordinates of the regulating system of the EOM. The emerging resonance cases are grouped, in which three of them are examined at once, and then the ME are established based on the examined resonance situations. When the stability of the system is examined by using RHC, several areas close to resonance ranges are plotted to assess the dynamic frequency response curves for various values of the acting parameters. The realized solutions are plotted for the same values of these parameters to observe how the system’s motion behaves. The various values of the excitation amplitudes, damping coefficients, and load resistance are displayed as functions of the output current time series and power. The results of this work can be used in a variety of applications, including automotive, ³⁴ health monitoring devices, ³⁵ wireless sensor networks, ³⁶ and civil engineering. ³⁷

Description of the dynamic system

This section describes the investigated system, which consists of two parts: a dynamic part and a harvester device. The dynamical part consists of masses $m_k\hspace{0.17em}\left(k=\mathrm{1,2}\right)$ connected with two damped springs of normal lengths ${\ell }_1$ and $\ell$ linear stiffness coefficient $k_k$ , and damping coefficients $C_k$ besides $C_3$ , see Figure 1. One of these masses is directed horizontally, and the other has an angle of rotation $\theta$ . An EH device is located at the damped spring of mass $m_2$ . The construction of this system can be estimated according to the magnet's oscillation within the coil of this device. As a result, an electric potential is generated through electromotive force, which is harnessed to harvest energy. It is hypothesized that the dynamical motion of this system is impacted by the external harmonic forces $F_1=F_x\hspace{0.17em}\cos \left({\mathrm{\Omega }}_{1}t\right)$ (in the horizontal direction), $F_2=F_r\hspace{0.17em}\cos \left({\mathrm{\Omega }}_{1}t\right)$ , (along the pendulum’s arm) and the moment $M\left(t\right)=M_{\theta }\hspace{0.17em}\cos \left({\mathrm{\Omega }}_{3}t\right)$ (at the point $O$ ), in which $\left(F_x,F_r,F_{\theta }\right)$ and ${\mathrm{\Omega }}_j\hspace{0.17em}\left(j=\mathrm{1,2,3}\right)$ represent the amplitudes and frequencies of the acted forces and moment, respectively. The harvester is coupled with the spring of mass $m_2$ , which consists of two magnets that are fastened to the end of a cylindrical tube.

Figure 1.

The dynamical system.

Due to the fact that any dynamical system has a Lagrangian $L=T-V$ , where $T$ and $V$ represent the system's kinetic and potential energies, respectively. In accordance with the examined system, it is possible to write

$\begin{array}{l}L=\frac{1}{2}\{\left(m_1+m_2\right){x^{\prime }}^2-k_1{x}^2-k_2{r}^2+2\hspace{0.17em}{m}_{2}r\hspace{0.17em}{x}^{\prime }{\theta }^{\prime }\hspace{0.17em}\cos \hspace{0.17em}\theta +m_2{\left(\ell +r\right)}^2{{\theta }^{\prime }}^2+m_2[2r^{\prime }{x}^{\prime }\hspace{0.17em}\sin \hspace{0.17em}\theta \\ +{r^{\prime }}^2+2\ell \hspace{0.17em}\cos \hspace{0.17em}\theta \left(\mathrm{g}+x^{\prime }{\theta }^{\prime }\right)-2g\ell \left]\right\}.\end{array}$

(1)

Here, dash represents the derivatives regarding time $t$ while $x^{\prime },r^{\prime },$ and ${\theta }^{\prime }$ are the system’s generalized velocities.

Based on equation (1), the following Lagrange’s equations can be identified to determine the governing system EOM

$\frac{d}{dt}\left(\frac{\partial L}{\partial q_s^{\prime }}\right)-\frac{\partial L}{\partial q_s}=Q_{q_s};\hspace{0.17em}\left(q_s=x,r,\theta \right)$

(2)

where $Q_{q_s}$ are the forces that have been applied generally to the dynamical system under investigation, and they take the following forms

$\begin{array}{l}{Q}_x=F_x\hspace{0.17em}\cos \left({\mathrm{\Omega }}_{1}t\right)-C_1{x}^{\prime },\\ Q_r=F_r\hspace{0.17em}\cos \left({\mathrm{\Omega }}_{2}t\right)-C_2{r}^{\prime }-B{L}_{coil}I,\\ Q_{\theta }=M_{\theta }\hspace{0.17em}\cos \left({\mathrm{\Omega }}_{1}t\right)-C_3{\theta }^{\prime }.\end{array}$

(3)

Here, $I$ is the generating current, $L_{coil}$ refers the coil length and $B$ is the magnetic flux.

Substituting equations (1) and (3) into equation (2), we can write the EOM as follows

$\begin{array}{l}\left(m_1+m_2\right)x^{″}+k_{1}x+m_2\left\{\left[r^{″}-\left(\ell +r\right){{\theta }^{\prime }}^2\right]\sin \hspace{0.17em}\theta +\left[2r^{\prime }{\theta }^{\prime }+\left(\ell +r\right){\theta }^{″}\right]\cos \hspace{0.17em}\theta \right\}\\ +C_1{x}^{\prime }=F_x\hspace{0.17em}\cos \left({\mathrm{\Omega }}_{1}t\right),\\ m_2{r}^{″}+r\left(k_2-m_2{{\theta }^{\prime }}^2\right)+C_2{r}^{\prime }+m_2\left(x^{″}\hspace{0.17em}\sin \hspace{0.17em}\theta -\ell {{\theta }^{\prime }}^2\right)+B{L}_{coil}I=F_r\hspace{0.17em}\cos \left({\mathrm{\Omega }}_{2}t\right),\\ m_2\left(\ell +r\right){\theta }^{″}+m_2\left[g\ell \hspace{0.17em}\sin \hspace{0.17em}\theta +\left(\ell +r\right)\left(2r^{\prime }{\theta }^{\prime }+x^{″}\hspace{0.17em}\cos \hspace{0.17em}\theta \right)\right]+C_3{\theta }^{\prime }=M_{\theta }\hspace{0.17em}\cos \left({\mathrm{\Omega }}_{3}t\right).\end{array}$

(4)

For the electrical circuit, we can write the Kirchhoff Voltage law and Faraday's one of electromagnetic induction. ¹¹

$L_{ind}\hspace{0.17em}{I}^{\prime }+RI-B{L}_{coil}{r}^{\prime }=0,$

(5)

where $R$ denotes the total of the resistances of the coil $R_C$ and the load $R_L$ . Let us insert the below dimensionless parameters

$\begin{array}{l}\tau =\omega \hspace{0.17em}t,{\omega }^2=\frac{g}{\ell },W_1^2=\frac{k_1}{m_1+m_2},{\omega }_1=\frac{W_1}{\omega },W_2^2=\frac{k_2}{m_2},\\ {\omega }_2=\frac{W_2}{\omega },M=\frac{m_2}{m_1+m_2},c_1=\frac{C_1}{\left(m_1+m_2\right)\omega },\alpha =\frac{{\mathrm{{\rm B}}}^2{L}_{coil}^2}{m_2{\omega }^2{L}_{ind}},\\ c_2=\frac{C_2}{m_2\omega },c_3=\frac{C_3}{m_2{\ell }^2\omega },f_1=\frac{F_x}{\left(m_1+m_2\right)\ell {\omega }^2},f_2=\frac{F_r}{m_2\ell {\omega }^2},\\ f_3=\frac{M_{\theta }}{m_2{\ell }^2{\omega }^2},p_1=\frac{{\mathrm{\Omega }}_1}{\omega },p_2=\frac{{\mathrm{\Omega }}_2}{\omega },p_3=\frac{{\mathrm{\Omega }}_3}{\omega },\eta =\frac{R}{L_{ind}{w}_1},\end{array}$

(6)

besides the dimensionless coordinates $U=\frac{x}{\ell },Z=\frac{r}{\ell },V=\left(L_{ind}/B\hspace{0.17em}{L}_{coil}\ell \right)I,$ and $\theta =\theta$ .

Making use of the aforementioned parameters (6) and equation (5), one can rewrite the system (4) in their non-dimensional forms as follows

$\begin{array}{l}\ddot{U}+{\omega }_1^{2}U+c_1\dot{U}+M\left(\theta -\frac{{\theta }^3}{6}\right)\hspace{0.17em}\left[\hspace{0.17em}\ddot{Z}-\left(1+Z\right){\dot{\theta }}^2\right]+M\left(1-\frac{{\theta }^2}{2}\right)\left[2\dot{Z}\dot{\theta }+\left(1+Z\right)\ddot{\theta }\right]=f_1\hspace{0.17em}\cos \left(p_1\tau \right),\\ \ddot{Z}+{\omega }_2^{2}Z+c_2\dot{Z}-{\dot{\theta }}^2+\alpha \hspace{0.17em}V-Z{\dot{\theta }}^2+\left(\theta -\frac{{\theta }^3}{6}\right)\ddot{U}=f_2\hspace{0.17em}\cos \left(p_2\tau \right),\\ \ddot{\theta }+\theta +c_3\dot{\theta }-{\theta }^3/6+\ddot{U}+2\hspace{0.17em}\dot{Z}\dot{\theta }+2Z\dot{Z}\dot{\theta }+Z\left(1-\frac{{\theta }^2}{2}\right)\hspace{0.17em}\ddot{U}+2Z\hspace{0.17em}\ddot{\theta }+Z^2\ddot{\theta }=f_3\hspace{0.17em}\cos \left(p_3\tau \right),\\ \dot{V}+\mu V-\dot{Z}=0.\end{array}$

(7)

Here, dots are the derivatives regarding $\tau$ . The aforementioned system is categorized as a nonlinear second-order differential system.

The proposed perturbation method

The major purpose of this section is to use the MSM to achieve the essential approximate solutions for the EOM of the previous systems (7) up to higher-order approximations. To achieve this aim, the functions $\sin \hspace{0.17em}\theta$ and $\cos \hspace{0.17em}\theta$ in system (7) can be approximated using Taylor series until the third-order of their expansions, which are valid in an extremely small region surrounding the position of static equilibrium.^38,39 Therefore, system (7) can thus be rewritten as

$\begin{array}{l}\ddot{U}+{\omega }_1^2\hspace{0.17em}U+M\left(\theta -\frac{{\theta }^3}{6}\right)\left[\ddot{Z}-\left(1+Z\right){\dot{\theta }}^2\right]+M\left(1-\frac{{\theta }^2}{2}\right)\left[2\dot{Z}\dot{\theta }+\left(1+Z\right)\ddot{\theta }\right]=f_1\hspace{0.17em}\cos \left(p_1\tau \right)-c_1\dot{U},\\ Z+{\omega }_2^{2}Z+\alpha V-{\dot{\theta }}^2-Z{\dot{\theta }}^2+\left(\theta -\frac{{\theta }^3}{6}\right)\ddot{U}=f_2\hspace{0.17em}\cos \left(p_2\tau \right)-c_2\dot{Z},\\ \ddot{\theta }+\theta +2\dot{Z}\dot{\theta }-\frac{{\theta }^3}{6}+2Z\dot{Z}\dot{\theta }+\ddot{U}+Z\left(1-\frac{{\theta }^2}{2}\right)\ddot{U}+2Z\ddot{\theta }+Z^2\ddot{\theta }=f_3\hspace{0.17em}\cos \left(p_3\tau \right)-c_3\dot{\theta },\\ \dot{V}+\mu V-\dot{Z}=0.\end{array}$

(8)

We suppose that the generalized coordinates can be represented in terms of tiny parameter $0<\epsilon \ll 1$ . As a result, we look at three additional variables, $u,z,\phi ,$ and $v$ , as shown below

$U\left(\tau \right)=\epsilon u\left(\tau ;\epsilon \right),Z\left(\tau \right)=\epsilon z\left(\tau ;\epsilon \right)\hspace{0.17em}\theta \left(\tau \right)=\epsilon \hspace{0.17em}\phi \left(\tau ;\epsilon \right),V\left(\tau \right)=\epsilon v\left(\tau ;\epsilon \right).$

(9)

In accordance with the MSM, one may construct the functions $u,z,\phi ,$ and $v$ in powers of $\epsilon$ as ¹⁹

$\begin{array}{l}u={\displaystyle \sum }_{k=0}^2{\epsilon }^k{u}_{k+1}\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^3\right),\\ z={\displaystyle \sum }_{k=0}^2{\epsilon }^k{z}_{k+1}\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^3\right),\\ \phi ={\displaystyle \sum }_{k=0}^2{\epsilon }^k{\phi }_{k+1}\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^3\right),\\ v={\displaystyle \sum }_{k=0}^2{\epsilon }^k{v}_{k+1}\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^3\right).\end{array}$

(10)

Here ${\tau }_n={\epsilon }^n\tau \hspace{0.17em}\left(n=\mathrm{0,1,2}\right)$ symbolize various time scales, with ${\tau }_0$ and ${\tau }_1,{\tau }_2$ are, respectively, defined as fast and slow scales. Consequently, the chain rule shown below may be employed to transform the $\tau$ -varying derivatives of ${\tau }_n$ as

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

(11)

It is important to recognize that terms beginning with $O\left({\epsilon }^3\right)$ and higher are careless due to their smallness.

To achieve the required solutions, the amplitudes $f_j\hspace{0.17em}\left(j=\mathrm{1,2,3}\right)$ , damping coefficients $c_j$ , and the parameter $M$ might be specified as

$f_j={\epsilon }^3{\stackrel{\sim }{f}}_j,c_j={\epsilon }^2{\stackrel{\sim }{c}}_j,M=\epsilon m\hspace{0.17em}\left(j=\mathrm{1,2,3}\right)$

(12)

Substituting (9) - (12) into (8) and equaling coefficients of different powers of $\epsilon$ in both sides to get the below three groups of partial differential equations (PDE)

Order of $\epsilon$

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

(13)

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

(14)

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

(15)

$\frac{\partial v_1}{\partial {\tau }_0}+\mu \hspace{0.17em}{v}_1=\frac{\partial z_1}{\partial {\tau }_0}.$

(16)

Order of ${\epsilon }^2$

$\frac{{\partial }^2{u}_2}{\partial {\tau }_0^2}+{\omega }_1^2{u}_2=-2\frac{{\partial }^2{u}_1}{\partial {\tau }_0\partial {\tau }_1}-m\frac{{\partial }^2{\phi }_1}{\partial {\tau }_0^2},$

(17)

$\frac{{\partial }^2{z}_2}{\partial {\tau }_0^2}+{\omega }_2^2{z}_2={\left(\frac{\partial {\phi }_1}{\partial {\tau }_0}\right)}^2-2\frac{{\partial }^2{z}_1}{\partial {\tau }_0\partial {\tau }_1}-{\phi }_1\frac{{\partial }^2{u}_1}{\partial {\tau }_0^2},$

(18)

$\frac{{\partial }^2{\phi }_2}{\partial {\tau }_0^2}+{\phi }_2+\frac{{\partial }^2{u}_2}{\partial {\tau }_0^2}=-\frac{\partial z_1}{\partial {\tau }_0}\frac{\partial {\phi }_1}{\partial {\tau }_0}-2\left(\frac{{\partial }^2{u}_1}{\partial {\tau }_0\partial {\tau }_1}+\frac{{\partial }^2{\phi }_1}{\partial {\tau }_0\partial {\tau }_1}\right)-z_1\left(\frac{{\partial }^2{u}_1}{\partial {\tau }_0^2}+2\frac{{\partial }^2{\phi }_1}{\partial {\tau }_0^2}\right),$

(19)

$\frac{\partial v_2}{\partial {\tau }_0}+\mu v_2=\frac{\partial z_1}{\partial {\tau }_1}-\frac{\partial v_1}{\partial {\tau }_1}+\frac{\partial z_2}{\partial {\tau }_0}.$

(20)

Order of ${\epsilon }^3$

$\begin{array}{l}\frac{{\partial }^2{u}_3}{\partial {\tau }_0^2}+{\omega }_1^2{u}_3={\stackrel{\sim }{f}}_1\hspace{0.17em}\cos \left(p_1\tau \right)-\frac{{\partial }^2{u}_1}{\partial {\tau }_1^2}-{\stackrel{\sim }{c}}_1\frac{\partial u_1}{\partial {\tau }_0}-2m\left(\frac{\partial z_1}{\partial {\tau }_0}\frac{\partial {\phi }_1}{\partial {\tau }_0}+\frac{{\partial }^2{\phi }_1}{\partial {\tau }_0\partial {\tau }_1}\right)-2(\frac{{\partial }^2{u}_1}{\partial {\tau }_0\partial {\tau }_2}\\ +\frac{{\partial }^2{u}_2}{\partial {\tau }_0\partial {\tau }_1})-m\left({\phi }_1\frac{{\partial }^2{z}_1}{\partial {\tau }_0^2}+z_1\frac{{\partial }^2{\phi }_1}{\partial {\tau }_0^2}+\frac{{\partial }^2{\phi }_2}{\partial {\tau }_0^2}\right),\end{array}$

(21)

$\begin{array}{l}\frac{{\partial }^2{z}_3}{\partial {\tau }_0^2}+{\omega }_2^2{z}_3={\stackrel{\sim }{f}}_2\hspace{0.17em}\cos \left(p_2\tau \right)-\alpha v_1-\frac{{\partial }^2{z}_1}{\partial {\tau }_1^2}-{\stackrel{\sim }{c}}_2\frac{\partial z_1}{\partial {\tau }_0}+2\left(\frac{\partial {\phi }_1}{\partial {\tau }_1}\frac{\partial {\phi }_1}{\partial {\tau }_0}+\frac{\partial {\phi }_1}{\partial {\tau }_0}\frac{\partial {\phi }_2}{\partial {\tau }_0}\right)+z_1{\left(\frac{\partial {\phi }_1}{\partial {\tau }_0}\right)}^2\\ -2\left(\frac{{\partial }^2{z}_1}{\partial {\tau }_0\partial {\tau }_2}+\frac{{\partial }^2{z}_1}{\partial {\tau }_0\partial {\tau }_1}\right)-{\phi }_1\left(\frac{{\partial }^2{u}_2}{\partial {\tau }_0^2}+2\frac{{\partial }^2{u}_1}{\partial {\tau }_0\partial {\tau }_1}\right)-{\phi }_2\frac{{\partial }^2{u}_1}{\partial {\tau }_0^2},\end{array}$

(22)

$\begin{array}{l}\frac{{\partial }^2{\phi }_3}{\partial {\tau }_0^2}+{\phi }_3={\stackrel{\sim }{f}}_3\hspace{0.17em}\cos \hspace{0.17em}\left(p_3\hspace{0.17em}\tau \right)-{\stackrel{\sim }{c}}_3\frac{\partial {\phi }_1}{\partial {\tau }_0}+\frac{1}{6}{\phi }_1^3-\frac{{\partial }^2{u}_1}{\partial {\tau }_1^2}-\frac{{\partial }^2{\phi }_1}{\partial {\tau }_1^2}-2(\frac{\partial {\phi }_1}{\partial {\tau }_1}\frac{\partial z_1}{\partial {\tau }_0}+\frac{{\partial }_1}{\partial {\tau }_1}\frac{\partial {\phi }_1}{\partial {\tau }_0}\\ +\frac{\partial z_2}{\partial {\tau }_0}\frac{\partial {\phi }_1}{\partial {\tau }_0}+\frac{\partial z_1}{\partial {\tau }_0}\frac{\partial {\phi }_2}{\partial {\tau }_0}+\frac{{\partial }^2{u}_1}{\partial {\tau }_0\partial {\tau }_2}+\frac{{\partial }^2{\phi }_1}{\partial {\tau }_0\partial {\tau }_2}+\frac{{\partial }^2{u}_2}{\partial {\tau }_0\partial {\tau }_1}+\frac{{\partial }^2{\phi }_2}{\partial {\tau }_0\partial {\tau }_1})-z_2\frac{{\partial }^2{u}_1}{\partial {\tau }_0^2}\\ -2\hspace{0.17em}{z}_1\left(\frac{{\partial }^2{u}_1}{\partial {\tau }_0\partial {\tau }_1}+\frac{\partial z_1}{\partial {\tau }_0}\frac{\partial {\phi }_1}{\partial {\tau }_0}+\frac{{\partial }^2{\phi }_2}{\partial {\tau }_0^2}+\frac{{\partial }^2{\phi }_1}{\partial {\tau }_0\partial {\tau }_1}+\frac{1}{2}\frac{{\partial }^2{u}_2}{\partial {\tau }_0^2}\right),\end{array}$

(23)

$\frac{\partial v_3}{\partial {\tau }_0}+\mu v_3=\frac{\partial z_3}{\partial {\tau }_0}-\left(\frac{\partial v_2}{\partial {\tau }_1}+\frac{\partial v_1}{\partial {\tau }_2}\right)+\frac{\partial z_2}{\partial {\tau }_1}+\frac{\partial z_1}{\partial {\tau }_2}.$

(24)

It is interesting to note that the forgoing equations (13) to (24) are of orders $\epsilon$ comprise 12 PDEs that can be solved consecutively. To accomplish this goal, let us start with the following general solutions of the homogenous equations (13) to (16)

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

(25)

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

(26)

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

(27)

The substitution of (26) into equation (16), yields a PDE that can be solved to have the form

$v_1=B_1{e}^{-\mu {\tau }_0}+\frac{i\hspace{0.17em}{\omega }_2}{\mu +i\hspace{0.17em}{\omega }_2}{A}_2\hspace{0.17em}{e}^{i\hspace{0.17em}{\omega }_2{\tau }_0}-\frac{i{\omega }_2}{\mu -i\hspace{0.17em}{\omega }_2}{\overline{A}}_2\hspace{0.17em}{e}^{-i\hspace{0.17em}{\omega }_2{\tau }_0};\hspace{0.17em}i=\sqrt{-1}.$

(28)

Here $A_j\hspace{0.17em}\left(j=\mathrm{1,2,3}\right)$ and $B_1$ are undetermined complex functions, while ${\overline{A}}_j$ are the related complex conjugates.

Substituting the first-order solutions (25) to (28) into the PDE of the second group (17) and (20) and then cancelling the emerging secular terms to get the second-order solutions as follows

$u_2=-\frac{m\hspace{0.17em}{A}_3}{{\omega }_1^2-1}{e}^{{\tau }_0\dot{i}}+CC,$

(29)

$\begin{array}{l}{\mathrm{z}}_2=\frac{e^{i{\tau }_0\left(1+{\omega }_1\right)}}{{\omega }_2^2-{\left(1+{\omega }_1\right)}^2}{\omega }_1^2{A}_1{A}_3-\frac{2e^{i{\tau }_0\left(1+{\omega }_1\right)}{\omega }_1{A}_1{A}_3}{\left[{\omega }_2^2-{\left(1+{\omega }_1\right)}^2\right]\left(1-{\omega }_1^2\right)}-\frac{e^{2i{\tau }_0}}{{\omega }_2^2-4} A_3^2+\frac{2{\omega }_1^2{\overline{A}}_1{A}_1}{\left(1-{\omega }_1^2\right){\omega }_2^2}\\ +\frac{e^{i{\tau }_0\left(1-{\omega }_1\right)}}{{\omega }_2^2-{\left({\omega }_1-1\right)}^2}{\omega }_1^2{\overline{A}}_1{A}_3+\frac{2e^{i{\tau }_0\left({\omega }_1-1\right)}{\omega }_1{\overline{A}}_3{A}_1}{\left(1-{\omega }_1^2\right)\left[{\omega }_2^2-{\left({\omega }_1-1\right)}^2\right]}+\frac{2}{{\omega }_2^2}{\overline{A}}_3{A}_3+CC,\end{array}$

(30)

$\begin{array}{l}{\phi }_2=\frac{2{\omega }_1{\omega }_2{\overline{A}}_2{A}_1{e}^{i\left({\omega }_1-{\omega }_2\right){\tau }_0}}{\left[1-{\left({\omega }_1-{\omega }_2\right)}^2\right]\left(1-{\omega }_1^2\right)}-\frac{4{\omega }_2{A}_2{A}_3{e}^{i\left(1+{\omega }_2\right){\tau }_0}}{1-{\left(1+{\omega }_2\right)}^2}+\frac{\left(2-{\omega }_1^2\right)}{{\left(1-{\omega }_1^2\right)}^3}{\omega }_1^{2}m{A}_1{e}^{i{\tau }_0{\omega }_1}\\ -\frac{2{\omega }_1{\omega }_2{A}_1{A}_2{e}^{i\left({\omega }_1+{\omega }_2\right){\tau }_0}}{\left[1-{\left({\omega }_1+{\omega }_2\right)}^2\right]\left(1-{\omega }_1^2\right)}-\frac{{\omega }_1^2{\overline{A}}_1{A}_2{e}^{i\left({\omega }_2-{\omega }_1\right){\tau }_0}}{1-{\left({\omega }_1-{\omega }_2\right)}^2}-\frac{2{\omega }_1^2{\overline{A}}_2{A}_1{e}^{i\left({\omega }_1-{\omega }_2\right){\tau }_0}}{\left[1-{\left({\omega }_1-{\omega }_2\right)}^2\right]\left(1-{\omega }_1^2\right)}\\ -\frac{2\left({\omega }_2-1\right)}{\left({\omega }_2-2\right){\omega }_2}{\overline{A}}_3{A}_2{e}^{i\left({\omega }_2-1\right){{\tau }_0}_2}+\frac{{\omega }_1^2\left({\omega }_1^2-3\right)}{\left({\omega }_1^2-1\right)\left({\omega }_2^2+2{\omega }_1{\omega }_2+{\omega }_1^2-1\right)}{A}_1{A}_2{e}^{i\left({\omega }_1+{\omega }_2\right){\tau }_0}+c.c,\end{array}$

(31)

$\begin{array}{l}{v}_2=\frac{{\overline{A}}_1{\overline{A}}_{3}i\left(1+{\omega }_1\right)e^{-i{\tau }_0\left(1+{\omega }_1\right)}}{\left[\mu -i\left(1+{\omega }_1\right)\right]\left[{\omega }_2^2-{\left(1+{\omega }_1\right)}^2\right]}\left[{\omega }_1^2+\frac{2{\omega }_1}{\left({\omega }_1^2-1\right)}\right]+\frac{2i{e}^{2i{\tau }_0}{A_3}^2}{\left({\omega }_2^2-4\right)\left(\mu +2i\right)}\\ +\frac{2i{e}^{2i{\tau }_0{\omega }_1}{\omega }_1^5{A_1}^2}{{\left({\omega }_1^2-1\right)}^2\left({\omega }_2^2-4{\omega }_1^2\right)\left(\mu +2i{\omega }_1\right)}+\frac{i{e}^{-i{\tau }_0\left({\omega }_1-1\right)}\left({\omega }_1-1\right)}{\left[\mu -i\left({\omega }_1-1\right)\right]}\\ \times A_3\left\{\frac{{\omega }_1^2{\overline{A}}_1}{{\omega }_2^2-{\left({\omega }_1-1\right)}^2}+\frac{2{\omega }_1{A}_1}{\left({\omega }_1^2-1\right)\left[{\omega }_2^2-{\left({\omega }_1+1\right)}^2\right]}\right\}+CC,\end{array}$

(32)

where $CC$ means the complex conjugates of the preceding terms.

Using the four equations provided above, we can determine the prerequisites for deleting secular terms from the second approximation as follows

$\frac{m\hspace{0.17em}{\omega }_1{A}_1}{1-{\omega }_1^2}-2i\frac{\partial A_1}{\partial {\tau }_1}=0,$

(33)

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

(34)

$\frac{m\hspace{0.17em}{\omega }_1{A}_3}{1-{\omega }_1^2}-2i\frac{\partial A_3}{\partial {\tau }_1}=0,$

(35)

$\frac{\partial B_1}{\partial {\tau }_1}=0,$

(36)

where the functions $A_2$ and $B_1$ depend on ${\tau }_2$ only.

The third-order solutions can be obtained by substituting (25) and (28) and (29) and (36) into (21) and (24) besides the removal of secular terms. The following conditions must be met to cancel these terms

$2i{\omega }_1\frac{\partial A_1}{\partial {\tau }_2}+i{\stackrel{\sim }{c}}_1{\omega }_1{A}_1+\frac{m^2{\omega }_1^2\left(4{\omega }_1^4-5{\omega }_1^2-3\right)A_1}{4{\left({\omega }_1^2-1\right)}^3}=0,$

(37)

$2i{\omega }_2\frac{\partial A_2}{\partial {\tau }_2}+i{\stackrel{\sim }{c}}_2{\omega }_2{A}_2+\frac{i\alpha {\omega }_2}{\mu +i{\omega }_2}{A}_2-2A_3{\overline{A}}_3{A}_2+2D_1=0,$

(38)

$\begin{array}{l}2i\frac{\partial A_3}{\partial {\tau }_2}+i{\stackrel{\sim }{c}}_3{A}_3+\frac{m^2{A}_3}{{\left({\omega }_1^2-1\right)}^2}-\frac{\left({\omega }_2^4+8{\omega }_2^2-32\right)\overline{A_3}{A_3}^2}{2{\omega }_2^2\left({\omega }_2^2-4\right)}-D_2{A}_3{\overline{A}}_1{A}_1\\ +\frac{4\left({\omega }_2^2-2\right)\overline{A_2}{A}_2{A}_3}{{\omega }_2^2-4}=0,\end{array}$

(39)

$\frac{\partial B_1}{\partial {\tau }_2}-\frac{\stackrel{\sim }{\alpha }\mu B_1}{{\mu }^2+{\omega }_2^2}=0,$

(40)

where

$D_1=\frac{{\omega }_1^2\left[6{\omega }_1^4+{\omega }_1^8+{\omega }_1^6\left({\omega }_2^2-5\right)+{\omega }_1^2\left({\omega }_2^2-1\right)-{\left({\omega }_2^2-1\right)}^2\right]}{{\left({\omega }_1^2-1\right)}^2\left[{\omega }_1^4+{\left({\omega }_2^2-1\right)}^2-2{\omega }_1^2\left({\omega }_2^2+1\right)\right]},$

(41)

$D_2=\frac{\left[{\omega }_2^2{\left({\omega }_2^2-1\right)}^2-\left[2{\omega }_1^4\left(4{\omega }_2^4-17{\omega }_2^2-16\right)+2{\omega }_1^6\left({\omega }_2^2-4\right)\left({\omega }_2^2+1\right)-{\omega }_1^2\left(2{\omega }_2^4+61{\omega }_2^2+40\right)+2\left(11{\omega }_2^4-21{\omega }_2^2+8\right)\right]{\omega }_1^2-2{\omega }_1^{10}{\omega }_2^2\right]}{{\left({\omega }_1^2-1\right)}^2{\omega }_2^2\left[{\omega }_1^4+{\left({\omega }_2^2-1\right)}^2-2{\omega }_1^2\left({\omega }_2^2+1\right)\right]}$

(42)

Moreover, the third-order approximations have the forms

$\begin{array}{l}{u}_3=\frac{{\stackrel{\sim }{f}}_1}{2\left({\omega }_1^2-{p_1}^2\right)}{e}^{i\tau p_1}+\frac{e^{i{\tau }_0\left(1+{\omega }_2\right)}}{\left[{\omega }_1^2-{\left(1+{\omega }_2\right)}^2\right]}\left[\frac{m{\left(1+{\omega }_2\right)}^2\left({\omega }_2-2\right)A_2{A}_3}{{\omega }_2}\right]-\frac{e^{i{\tau }_0\left({\omega }_1+{\omega }_2\right)}}{\left[{\omega }_1^2-{\left({\omega }_1+{\omega }_2\right)}^2\right]}\\ \times \left[\frac{m{\left({\omega }_1+{\omega }_2\right)}^2\left({\omega }_1^4+{\omega }_2^2-2{\omega }_1^2-1\right)A_1{A}_2}{\left({\omega }_1^2-1\right)\left({\omega }_1+{\omega }_2-1\right)\left({\omega }_1+{\omega }_2+1\right)}\right]-\left[\frac{m{\left({\omega }_1-{\omega }_2\right)}^2\left({\omega }_1^4+{\omega }_2^2-2{\omega }_1^2-1\right)A_1{\overline{A}}_2}{\left({\omega }_1^2-1\right)\left[{\left({\omega }_1-{\omega }_2\right)}^2-1\right]}\right]\\ \times \frac{e^{i{\tau }_0\left({\omega }_1-{\omega }_2\right)}}{\left[{\omega }_1^2-{\left({\omega }_1-{\omega }_2\right)}^2\right]}+\frac{m{\overline{A}}_3{A}_2{\left({\omega }_2-1\right)}^2\left({\omega }_2^2-2\right)}{{\omega }_2\left[{\omega }_1^2-{\left({\omega }_2-1\right)}^2\right]}{e}^{i{\tau }_0\left({\omega }_2-1\right)}-\frac{m^2{A}_3{e}^{i{\tau }_0}}{{\left({\omega }_1^2-1\right)}^2}\left(1-\frac{1}{{\omega }_1^2-1}\right)+C.C,\end{array}$

(43)

$z_3=z_{31}+z_{32}+z_{33}+z_{34},$

(44)

${\phi }_3={\phi }_{31}+{\phi }_{32}+{\phi }_{33},$

(45)

$\begin{array}{l}{v}_3=\frac{e^{-i{p}_2{\tau }_0}i\stackrel{\sim }{f_2}{p}_2}{2\left({\omega }_2^2-p_2^2\right)\left(\mu -i{p}_2\right)}-\frac{e^{-\mu {\tau }_0}\alpha \mu B_1}{{\mu }^2+{\omega }_2^2}+\frac{i{e}^{-i{\tau }_0\left({\omega }_2-2\right)}\left(3{\omega }_2^2-6{\omega }_2+4\right){\overline{A}}_2{A_3}^2}{{\omega }_2{({\omega }_2^2-\left({\omega }_2-2\right)}^2)]\left(\mu -i\left({\omega }_2-2\right)\right)}\\ -\frac{i{e}^{-i{\tau }_0\left(2{\omega }_1+{\omega }_2\right)}{\omega }_1^2\left(2{\omega }_1+{\omega }_2\right){({\omega }_1^2\left({\omega }_1^2-1\right)}^2-6{\omega }_1{\omega }_2-3{\omega }_2^2-1){\overline{A}}_1^2{\overline{A}}_2}{{\left({\omega }_1^2-1\right)}^2\left({\omega }_1+{\omega }_2-1\right)\left({\omega }_1+{\omega }_2+1\right){({\omega }_2^2-\left(2{\omega }_1+{\omega }_2\right)}^2)\left(\mu -i\left(2{\omega }_1+{\omega }_2\right)\right)}\\ +\frac{i{e}^{-i{\tau }_0\left(2{\omega }_1-{\omega }_2\right)}{\omega }_1^2\hspace{0.17em}\left(2{\omega }_1-{\omega }_2\right)\left(2{\omega }_1^2+2{\omega }_1^4-{\omega }_1^6-6{\omega }_1{\omega }_2+3{\omega }_2^2+1\right){\overline{A}}_1^2{A}_2}{{\left({\omega }_1^2-1\right)}^2{(\left({\omega }_1-{\omega }_2\right)}^2-1){({\omega }_2^2-\left(2{\omega }_1-{\omega }_2\right)}^2)\left(\mu -i\left(2{\omega }_1-{\omega }_2\right)\right)}\\ +\frac{i{e}^{i{\tau }_0\left(2{\omega }_1+{\omega }_2\right)}{\omega }_1^2\left(2{\omega }_1+{\omega }_2\right)\left(2{\omega }_1^2+2{\omega }_1^4-{\omega }_1^6+6{\omega }_1{\omega }_2+3{\omega }_2^2+1\right){A_1}^2{A}_2}{{\left({\omega }_1^2-1\right)}^2\left({\omega }_1+{\omega }_2-1\right)\left({\omega }_1+{\omega }_2+1\right){({\omega }_2^2-\left(2{\omega }_1+{\omega }_2\right)}^2)\left(\mu +i\left(2{\omega }_1+{\omega }_2\right)\right)}\\ -\frac{i{e}^{-i{\tau }_0\left({\omega }_2+2\right)}\left(4+5{\omega }_2\left({\omega }_2+2\right)\right){\overline{A}}_2{\overline{A}}_3^2}{{\omega }_2{({\omega }_2^2-\left({\omega }_2+2\right)}^2)\left(\mu -i\left({\omega }_2+2\right)\right)},\end{array}$

(46)

where $z_{3n}\hspace{0.17em}\left(n=\mathrm{1,2,3,4}\right)$ and ${\phi }_{3j}\left(j=\mathrm{1,2,3}\right)$ are functions of $A_{j0}$ and ${\omega }_{j0}$ , see Appendix I.

Evaluation of resonance

The main topic of the current section is to classify the evolving resonance cases. These cases could arise when the denominators of higher approximations equal null. ³⁶ Therefore, according to second and third approximations (29) and (32) and (43) and (46), one may obtain the following resonance cases

(a). Internal resonance occurs when one of the next conditions is satisfied

${\omega }_1\approx \pm 1,{\omega }_2\approx \pm 2\hspace{0.17em}{\omega }_2\approx 0.$

• (b). When $p_1\approx {\omega }_1,p_2\approx {\omega }_2,p_3\approx 1,$ there is a chance for external (primary) resonance.

It should be noted that the behavior of the studied system will be quite challenging if any one of these conditions is met. Additionally, the approximate solutions achieved above will still hold if the oscillations have values that are not very similar to the resonance situations. To examine the stability of the system, the requirements of solvability and ME for this system must be obtained. Consequently, we will simultaneously examine the three cases of primary external resonance. Specifically, we take into account $p_1\approx {\omega }_1,p_2\approx {\omega }_2,$ and $p_3\approx 1,$ that expresses how close $p_1,p_2,$ and $p_3$ are to ${\omega }_1,{\omega }_2,$ and $1$ , respectively. Following that, we add the detuning parameters ${\sigma }_j=\epsilon \hspace{0.17em}{\stackrel{\sim }{\sigma }}_j\hspace{0.17em}\left(j=1,\mathrm{2,3}\right)$ as follows

$p_1={\omega }_1+{\sigma }_1,p_2={\omega }_2+{\sigma }_2,p_3=1+{\sigma }_3,{\sigma }_j=\epsilon \hspace{0.17em}{\stackrel{\sim }{\sigma }}_j\hspace{0.17em}\left(j=1,\mathrm{2,3}\right).$

(47)

We are now in the process of obtaining the solvability conditions. For this purpose, substitute (47) into the PDE equations (27) and (19) and (21) and (23), then eliminating the secular terms, the conditions of solvability can be gotten easily as follows

- For the PDE of the second-order of approximation

$\frac{m\hspace{0.17em}{\omega }_1{A}_1}{1-{\omega }_1^2}-2i\frac{\partial A_1}{\partial {\tau }_1}=0,$

(48)

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

(49)

$\frac{m\hspace{0.17em}{\omega }_1{A}_3}{1-{\omega }_1^2}-2i\frac{\partial A_3}{\partial {\tau }_1}=0,$

(50)

$\frac{\partial B_1}{\partial {\tau }_1}=0,$

(51)

• - For the PDE of the third-order of approximation

${\stackrel{\sim }{f}}_1\hspace{0.17em}{e}^{i{\tau }_1{\stackrel{\sim }{\sigma }}_1}-4\mathrm{i}{\omega }_1\frac{\partial A_1}{\partial {\tau }_2}-2i{\stackrel{\sim }{c}}_1{\omega }_1{A}_1-\frac{m^2{\omega }_1^2\left(4{\omega }_1^4-5{\omega }_1^2-3\right)A_1}{2{\left({\omega }_1^2-1\right)}^3}=0,$

(52)

${\stackrel{\sim }{f}}_2\hspace{0.17em}{e}^{i{\tau }_1{\stackrel{\sim }{\sigma }}_2}-4i{\omega }_2\frac{\partial A_2}{\partial {\tau }_2}-2i{\stackrel{\sim }{c}}_2{\omega }_2{A}_2-2\frac{i\alpha {\omega }_2}{\mu +i{\omega }_2}{A}_2+4A_3{\overline{A}}_3{A}_2-4D_1=0,$

(53)

$\begin{array}{l}{\stackrel{\sim }{f}}_3\hspace{0.17em}{e}^{i{\tau }_1{\stackrel{\sim }{\sigma }}_3}-4i\frac{\partial A_3}{\partial {\tau }_2}-2i{\stackrel{\sim }{c}}_3{A}_3-2\frac{m^2{A}_3}{{\left({\omega }_1^2-1\right)}^2}+\frac{\left({\omega }_2^4+8{\omega }_2^2-32\right){\overline{A}}_3{A_3}^2}{{\omega }_2^2\left({\omega }_2^2-4\right)}+2D_2{A}_3{\overline{A}}_1{A}_1\\ -\frac{8\left({\omega }_2^2-2\right){\overline{A}}_2{A}_2{A}_3}{{\omega }_2^2-4}=0,\end{array}$

(54)

$\frac{\partial B_1}{\partial {\tau }_2}-\frac{\stackrel{\sim }{\alpha }\mu B_1}{{\mu }^2+{\omega }_2^2}=0,$

(55)

From the equations (55) and (51), we can express $B_1$ in the form

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

(56)

where $B_2$ is an unknown complex function of the higher slow time scales than ${\tau }_2$ .

We may identify the functions $A_j$ , which depend exclusively on ${\tau }_2$ , by using the previously solvability requirements (52) to (54). As a result, we take into account the polar forms $A_j$ as follows

$A_j=\frac{1}{2}{\stackrel{\sim }{h}}_j\hspace{0.17em}{e}^{i\hspace{0.17em}{\stackrel{\sim }{\psi }}_j},h_j=\epsilon {\stackrel{\sim }{h}}_j\hspace{0.17em}\left(j=1,\mathrm{2,3}\right).$

(57)

Here $\stackrel{\sim }{h}$ and $\stackrel{\sim }{\psi }$ express the amplitude and phase of the solutions $u,z,$ and $\phi$ . Consider the below adjusted phases

${\theta }_j={\tau }_1\hspace{0.17em}{\stackrel{\sim }{\sigma }}_j-{\psi }_j;\hspace{0.17em}j=\left(\mathrm{1,2,3}\right).$

(58)

Making use of (57) and (58) into (52) to (54), and then separating the real portions and the imaginary ones to get the following ME

$\begin{array}{l}\frac{d{h}_1}{d\tau }=\frac{f_1}{2{\omega }_1}\sin \hspace{0.17em}{\theta }_1-\frac{c_1{h}_1}{2},\\ h_1\frac{d{\theta }_1}{d\tau }=\frac{f_1}{2{\omega }_1}\cos \hspace{0.17em}{\theta }_1+h_1{\sigma }_1+\frac{M{\omega }_1}{2{\left(1-{{\omega }_1}^2\right)}^2}{h}_1-\frac{M^2{\omega }_1\left(4{{\omega }_1}^4-5{{\omega }_1}^2-3\right)}{8{\left({{\omega }_1}^2-1\right)}^3}{h}_1,\\ \frac{d{h}_2}{d\tau }=\frac{f_2}{2{\omega }_2}\sin \hspace{0.17em}{\theta }_2-\frac{c_2{h}_2}{2}-\frac{\alpha \mu }{2\left({\mu }^2+{{\omega }_2}^2\right)}{h}_2,\\ h_2\frac{d{\theta }_2}{d\tau }=\frac{f_2}{2{\omega }_2}\cos \hspace{0.17em}{\theta }_2-\frac{\alpha {\omega }_2}{2\left({\mu }^2+{\omega }_2^2\right)}{h}_2+h_2{\sigma }_2+\frac{h_3^2{h}_2}{4{\omega }_2}-\frac{D_1{h}_1^2{h}_2}{4{\omega }_2},\\ \frac{d{h}_3}{d\tau }=\frac{f_3}{2}\sin \hspace{0.17em}{\theta }_3-\frac{c_3{h}_3}{2},\\ h_3\frac{d{\theta }_3}{d\tau }=\frac{f_3}{2}\cos \hspace{0.17em}{\theta }_3+\frac{M}{2{\left({{\omega }_1}^2-1\right)}^2}{h}_3-\frac{M^2}{2{\left({{\omega }_1}^2-1\right)}^2}{h}_3+h_3{\sigma }_3\\ +\frac{\left({{\omega }_2}^4+8{{\omega }_2}^2-32\right)h_3^3}{16{{\omega }_2}^2\left({{\omega }_2}^2-4\right)}+\frac{D_2{h}_1^2{h}_3}{8}-\frac{4\left({{\omega }_2}^2-2\right)h_2^2}{8\left({{\omega }_2}^2-4\right)}{h}_3.\end{array}$

(59)

The previous autonomous system ME (59) consists of six ordinary differential equations (ODE) of adjusted phases ${\theta }_j$ and amplitudes $h_j$ . These equations can be solved numerically through the planes $h_n\hspace{0.17em}{\theta }_n$ with the initial conditions $h_1\left(0\right)=0.001,$ ${\theta }_1\left(0\right)=0,$ $h_2\left(0\right)=0.001,$ ${\theta }_2\left(0\right)=0,$ $h_3\left(0\right)=0.001,$ and ${\theta }_3\left(0\right)=0$ . These solutions are graphed in Figures 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, and 13 for a few different values of the system's dimensional parameters as $R=100\hspace{0.17em}\mathrm{\Omega },$ $L_{ind}=0.1\mathrm{{\rm H}},$ $\ell =0.5m,$ $B=1\hspace{0.17em}\mathrm{{\rm T}},$ $L_{coil}=1m$ and dimensionless parameters, that in accordance with the following data

$\begin{array}{l}{f}_1=0.02\times {10}^{-3},f_2=0.034\times {10}^{-3},f_3=0.034\times {10}^{-3},c_1\left(=\mathrm{0.001,0.005,0.009}\right),\\ c_2\left(=\mathrm{0.04259,0.00212,0.02129}\right),c_3\left(=\mathrm{0.00358,0.0358,0.358}\right),p_1=12,\\ {\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right),{\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right),p_2=32,\\ p_3=1.1,{\sigma }_1=p_1-{\omega }_1,{\sigma }_2=p_2-{\omega }_2,{\sigma }_3=p_3-1,m_1=10,m_2=5.\end{array}$

Figure 2.

Presents the waves of $h_1\left(\tau \right)$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for various values of: (a) $c_1\left(=0.001,0.005,0.009\right)$ , (b) $c_2\left(=0.04259,0.0212,0.002129\right)$ , and (c) $c_3\left(=0.00358,0.0358,0.358\right)$ .

Figure 3.

Shows the waves of $h_1\left(\tau \right)$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for various values of: (a) ${\omega }_1\left(=5.5636,5.29112,4.73416\right)$ , and (b) ${\omega }_2\left(=13.924,12.0675,10.2109\right)$ .

Figure 4.

Presents the waves of $h_2\left(\tau \right)$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for various values of: (a) $c_1\left(=0.001,0.005,0.009\right)$ , (b) $c_2\left(=0.04259,0.0212,0.002129\right)$ , and (c) $c_3\left(=0.00358,0.0358,0.358\right)$ .

Figure 5.

Shows the waves of $h_3\left(\tau \right)$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for various values of: (a) ${\omega }_1\left(=5.5636,5.29112,4.73416\right)$ , and (b) ${\omega }_2\left(=13.924,12.0675,10.2109\right)$ .

Figure 6.

Illustrates the waves of $h_3\left(\tau \right)$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for various values of: (a) $c_1\left(=0.001,0.005,0.009\right)$ , (b) $c_2\left(=0.04259,0.0212,0.002129\right)$ , and (c) $c_3\left(=0.00358,0.0358,0.358\right)$ .

Figure 7.

Presents the waves of $h_3\left(\tau \right)$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for various values of: (a) ${\omega }_1\left(=5.5636,5.29112,4.73416\right)$ , and (b) ${\omega }_2\left(=13.924,12.0675,10.2109\right)$ .

Figure 8.

Describes the waves of ${\theta }_1\left(\tau \right)$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for various values of ${\theta }_1\left(\tau \right)$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for various values of: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 9.

Presents the waves of ${\theta }_1\left(\tau \right)$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for various values of: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

Figure 10.

Reveals the behavior of ${\theta }_2\left(\tau \right)$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for various values of: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 11.

Presents the behavior of ${\theta }_2\left(\tau \right)$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for various values of: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

Figure 12.

Presents the waves of ${\theta }_3\left(\tau \right)$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for various values of: (a) $c_1=\mathrm{0.001,0.005,0.009}$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 13.

Presents the waves of ${\theta }_3\left(\tau \right)$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for various values of: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

To examine the significance of acting parameters such as $c_j\left(j=\mathrm{1,2,3}\right)$ and ${\omega }_j\hspace{0.17em}\left(j=\mathrm{1,2}\right)$ on the behavior of the system, Figures 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, and 13 are plotted. These figures reveal the transient responses of the amplitudes $h_j\hspace{0.17em}\left(j=\mathrm{1,2,3}\right)$ and modified phases ${\theta }_j$ for system (59) of the above-mentioned specific values. Generally, one can note that the drawn curves in Figures 2, 3, 4, 5, 6, and 7 have decay manner over the examined time interval which reveals that the behaviors $h_j$ have stable procedures, while ${\theta }_j$ curves have various behaviors, as explored in Figures 8, 9, 10, 11, 12, and 13.

The curves in Figure 2 are drawn when the damping coefficients $c_j$ have different values and reveal the time variation of the amplitude $h_1$ over the time interval $\left[\mathrm{0,500}\right]$ . The decay behavior of these curves during this interval is expected from the structure of the dynamical system (59) of ME. It is observed that the plotted waves have low impulses with the variation of $c_1$ and $c_3$ , as seen in Figures 2(a) and (b)) in which the waves’ amplitudes decrease with the increase of $c_1$ and $c_3$ . On the other hand, there is no variation of the wave when $c_2$ has various values, as seen in Figure 2(b), which is due to the fact that the related equation does not depend on $c_2$ directly. A closer examination of the graphed curves in Figure 3 shows that they are calculated when of ${\omega }_1$ and ${\omega }_2$ take various values. It is noted that the values of ${\omega }_1$ have positive impact on the curves of modified amplitude $h_1$ , as seen in Figure 3(a), while the values of ${\omega }_2$ have no influence on the wave manner, as explored in Figure 3(b). The waves of the indicated figures have decay behavior during the whole-time interval; in addition to the amplitude of the waves decreases with the increase of ${\omega }_1$ values.

The curves of Figures 4 and 5 depict fascinating changes in the modified amplitude $h_2$ for variations of the damping coefficients and frequencies, respectively. In fact, as drawn in Figure 4, the time history of $h_2$ has a fast descent through a very small range of the amplitude during the investigated time interval, as illustrated in Figures 4(a) and (b). The amplitudes of the waves illustrating $h_2$ vary between increasing and decreasing, as indicated in Figure 4(a). The curves depicted in the portions of Figure 4 show that $c_2$ and $c_3$ have an excellent impact on the behavior of $h_2$ , while it has no influence on the variation of $c_1$ . The drawn waves in Figures 5(a) and (b) present the decay behavior of $h_2$ , in which it has no variation with ${\omega }_1$ , whereas a positive variation with ${\omega }_2$ is observed.

The main purpose of the graphed curves in Figures 6 and 7 is to examine the variation of the third amplitude $h_3\left(\tau \right)$ when $c_1,c_2,c_3,{\omega }_1,$ and ${\omega }_2$ have distinct values. These curves confirmed that the values of the parameters $c_1$ and $c_2$ have no influence on the wave of $h_3$ , as seen in Figures 6(a) and (b). On the contrary, the modified amplitude $h_3$ is impacted with the variation of the values of $c_3$ , as plotted in Figure 6(c). High oscillations are observed at the beginning of time interval when $\tau \in \left[\mathrm{0,25}\right]$ , and then steady behavior is noticed till the end of the time interval. The steady manner begins earlier when the value of $c_3$ increases. Figure 7 investigates the effect of various values of the frequencies ${\omega }_1$ and ${\omega }_2$ on the temporal course of $h_3$ . Decay waves have been observed with the change of ${\omega }_1$ , as indicated in Figure 7(a). In Figure 7(b), the waves rise to a small range and then oscillate periodically till the end of the interval $\left[\mathrm{0,500}\right]$ , in which the amplitude of these waves decreases with the increase of ${\omega }_2$ values.

Figures 8, 9, 10, 11, 12, and 13 illustrate the time histories of the modified phases ${\theta }_1,{\theta }_2,$ and ${\theta }_3$ when the influence of the aforementioned parameters is taken into account. It is seen that ${\theta }_1\left(\tau \right)$ is impacted by the variation of $c_1$ and ${\omega }_1$ , as seen in Figures 8(a) and 9(a), respectively. On the contrary, its behavior has no variation with the change of other parameters, as indicated in Figures 8(a) and (c) and 9(b). In general, the behavior of the explored waves increases with small oscillations during the whole time interval, which was expected from the equations of the system (59).

From the curves of Figures 10 and 11, one can notice that the process of changing the values of $c_2,c_3,{\omega }_1,$ and ${\omega }_2$ has a positive impact on the plotted waves of ${\theta }_2\left(\tau \right)$ . They increase gradually till a specified value and then have a stationary behavior over the entire examined interval, as seen in Figure 10(b), while the waves in Figure 10(c) and 11(a) and (b) increase to a certain value with the increase of time. Moreover, it is seen in Figure 10(a) that ${\theta }_2\left(\tau \right)$ has not been impacted by the change in $c_1$ value, which goes back to the structure of the system of equation (59).

The curves of Figures 12 and 13 clearly demonstrate the variation of the function ${\theta }_3\left(\tau \right)$ for the above mentioned chosen values of the used parameters $c_1,c_2,c_3,{\omega }_1,$ and ${\omega }_2$ . The oscillated waves in Figure 12(a) and (b) increase during the whole interval and have no variation with the various values of $c_1$ and $c_2$ . On the other hand, positive oscillations are noted, when $c_3$ and ${\omega }_1$ have various values, at the beginning of a time interval and then they have steady behavior till the end of this interval, as seen in Figures 12(c) and 13(a), respectively. The influence of ${\omega }_2$ values on the plotted waves in Figure 13(b) shows that they are raised until a certain range and then oscillate over the examined interval. According to the aforementioned figures, we can note that the behavior of the plotted waves that describe the modified amplitudes and phases has a stationary form and is free of chaos.

Referring to the plotted curves of the amplitudes $h_j\left(\tau \right)\hspace{0.17em}\left(j=\mathrm{1,2,3}\right)$ and the modified phases ${\theta }_j\left(\tau \right),$ one can present the projection of these curves in the planes ${\theta }_j{h}_j$ for the same considered values of the damping and frequency parameters, as seen in Figures 14, 15, 16, 17, 18, and 19. Decay waves are noted for the graphed curves in Figures 14, 15, 16, 17, periodic waves are observed in Figure 18(a) and (b), while spiral curves heading towards one point are drawn in Figures 18(a) and 19. We can say that the difference in the drawn types of curves is due to the complete consistency between the time trajectories of the amplitudes and their corresponding phases. Therefore, the dynamical motion of the system is stable over the studied time interval.

Figure 14.

Illustrates the phase plane ${\theta }_1{h}_1$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for various values of: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 15.

Shows the phase plane ${\theta }_1{h}_1$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for various values of: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

Figure 16.

Shows the phase plane ${\theta }_2{h}_2$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for various values of: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 17.

Presents the phase plane ${\theta }_2{h}_2$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for various values of: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

Figure 18.

Presents the phase plane ${\theta }_3{h}_3$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for various values of: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 19.

Illustrates the phase plane ${\theta }_3{h}_3$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for various values of: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

The main purpose of the explored curves in Figures 20, 21, 22, 23, 24, 25 is to investigate the time histories of the achieved asymptotic solutions $x,r,$ and $\theta$ , taking into account the above-selected values of the parameters $c_1,c_2,c_3,{\omega }_1,$ and ${\omega }_2$ . Decay waves are plotted in parts of Figure 20, in which the solution $x$ is impacted by the values of $c_1$ and the amplitude of the graphed waves decreases with the increase of $c_1$ values, as demonstrated in part (a). On the other side, this solution isn’t influenced by the selected values of $c_2$ and $c_3$ , as indicated in the other parts (b) and (c). Periodic waves describing the temporal behavior of this solution are graphed in Figures 21(a) and (b) when ${\omega }_1$ and ${\omega }_2$ have the same considered values above, in which their amplitudes increase with the decrease of ${\omega }_1$ and ${\omega }_2$ values. According to this simulation, one concludes that the plotted wave of the solution $x$ is stable during the whole-time interval.

Figure 20.

Reveals the solution $x\left(\tau \right)$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for various values of: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 21.

Shows the solution $x\left(\tau \right)$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for various values of: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

Figure 22.

The behavior of $r\left(\tau \right)$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for various values of: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 23.

The behavior of $r\left(\tau \right)$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for various values of: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

Figure 24.

The behavior of $\theta \left(\tau \right)$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for various values of: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ ,and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 25.

The behavior of $\theta \left(\tau \right)$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for various values of: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

The variation of the solution $r$ via the dimensionless time $\tau$ is graphed, respectively, in Figures 22 and 23 for various values of the damping and frequency parameters. The plotted waves behave in the form of decay wave packets, in which this solution is acted only by the values of $c_2$ . The decay reaches its maximum behavior when the damping value is high, as seen in parts of Figure 22. Periodic and quasi-periodic waves with a slight decay for the function $r\left(\tau \right)$ are graphed, respectively, in parts (a) and (b) of Figure 23. The wave’s amplitude decreases with the increase in ${\omega }_2$ values.

The time history of the asymptotic solution $\theta$ is drawn in Figures 24 and 25, when the damping parameters and frequency ones have various values. The drawn curves have stationary behavior over the examined interval of time, as seen in parts of these figures.

Figure 26 shows the time histories of both the achieved approximate and the numerical solutions over the interval $\tau \in \left[\mathrm{0,500}\right]$ . The represented waves of both solutions have periodic or quasi-periodic forms with little decay, which gives an induction about the stability of these waves. A closer look at this figure shows that there is a great deal of consistency between these solutions, which supports the importance of the used perturbation approach that produces accurate approximate solutions that are compatible with the numerical solutions of the original governing system.

Figure 26.

Describes the time history for both the numerical and approximate solutions.

The phase portraits curves for the obtained solutions $x,r,$ and $\theta$ are depicted in Figures 27, 28, 29, 30, 31 and 32. These curves are drawn for the same values of $c_j$ and ${\omega }_k$ to show closed curves whether symmetric or not that gives an induction about the stationary behavior of the obtained solutions. It is noticed that the phase portraits of the solution $\theta$ have deformation closed curves when $c_j$ and ${\omega }_k$ have different values. This deformation is due to that the equation of the first approximation of this solution is nonhomogeneous.

Figure 27.

Illustrates the phase portraits for the solution $x$ when ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 28.

Illustrates the phase portraits for the solution $x$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

Figure 29.

Describes the phase portraits for the solution $r$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 30.

Illustrates the phase portraits for the solution $r$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

Figure 31.

Presents the phase portraits for the solution $\theta$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 32.

Illustrates the phase portraits for the solution $\theta$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

Curves of the generated current $I$ by the EH device are represented graphically in Figures 33 and 34 for various values of the aforementioned parameters of damping and frequency. These curves have been graphed according to the current equation (5) to show their time behavior over the examined interval. Small amplitudes of the produced current are noted, resulting in low power levels (in milliwatts) at the output, which is considered suitable for sensor applications. The waveform of the current exhibits the existence of harmonic signals on the basic ones at the beginning of the wave's transient current caused by the RL circuit. The current then exhibits stationary behavior, which is necessary for usage. It is clear that the current has not been impacted by the change in $c_1$ , as seen in Figure 33(a), while it has been influenced by the values of $c_2$ and $c_3$ , as indicated in Figures 33(b) and (c). Based on the drawn curves, the current rises and oscillates at the beginning of the time interval until a certain time value, and then it has steady behavior. Increasing the value of the damping coefficient leads to faster wave stability than its lower value.

Figure 33.

Detects the time history of the current $I$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ for various values of: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 34.

Detects the time history of the current $I$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ for various values of: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

The good effect of changing the frequency values is clearly evident in Figure 34, in which periodic waves are graphed with a slight decay at the end of the time interval, and we find that this change seams slightly and clearly when ${\omega }_1$ and ${\omega }_2$ vary.

Steady-State solutions

The purpose of the present section is to investigate the model's vibrations in the steady-state scenario. In other words, when the transient process's vibrations disappear due to the presence of the system's damping. In these circumstances, the system of ME (59) plays a crucial role in defining this goal, in which the amplitudes $h_j\hspace{0.17em}\left(j=\mathrm{1,2,3}\right)$ and the adjusted phases ${\theta }_j$ have zero first derivatives. As a result, one can write

$\frac{d{h}_j}{dt}=0,\frac{d{\theta }_j}{dt}=0;\hspace{0.17em}\left(j=\mathrm{1,2,3}\right).$

(60)

According to the above conditions, equations of system (59) have the form

$\begin{array}{l}\frac{f_1}{2{\omega }_1}\sin \hspace{0.17em}{\theta }_1-\frac{c_1{h}_1}{2}=0,\\ \frac{f_1}{2{\omega }_1}\cos \hspace{0.17em}{\theta }_1+h_1{\sigma }_1+\frac{M{\omega }_1}{2{\left(1-{{\omega }_1}^2\right)}^2}{h}_1-\frac{M^2{\omega }_1\left(4{{\omega }_1}^4-5{{\omega }_1}^2-3\right)}{8{\left({{\omega }_1}^2-1\right)}^3}{h}_1=0,\\ \frac{f_2}{2{\omega }_2}\sin \hspace{0.17em}{\theta }_2-\frac{c_2{h}_2}{2}-\frac{\alpha \mu }{2\left({\mu }^2+{{\omega }_2}^2\right)}{h}_2=0,\\ \frac{f_2}{2{\omega }_2}\cos \hspace{0.17em}{\theta }_2-\frac{\alpha {\omega }_2}{2\left({\mu }^2+{\omega }_2^2\right)}{h}_2+h_2{\sigma }_2+\frac{h_3^2{h}_2}{4{\omega }_2}-\frac{D_1{h}_1^2{h}_2}{4{\omega }_2}=0,\\ \frac{f_3}{2}\sin \hspace{0.17em}{\theta }_3-\frac{c_3{h}_3}{2}=0,\\ \frac{f_3}{2}\cos \hspace{0.17em}{\theta }_3+\frac{M}{2{\left({{\omega }_1}^2-1\right)}^2}{h}_3-\frac{M^2}{2{\left({{\omega }_1}^2-1\right)}^2}{h}_3+h_3{\sigma }_3+\frac{D_2{h}_1^2{h}_3}{8}\\ +\frac{\left({{\omega }_2}^4+8{{\omega }_2}^2-32\right)h_3^3}{16{{\omega }_2}^2\left({{\omega }_2}^2-4\right)}-\frac{4\left({{\omega }_2}^2-2\right)h_2^2}{8\left({{\omega }_2}^2-4\right)}{h}_3=0.\end{array}$

(61)

After we get rid of the modified phases ${\theta }_j$ from equation (61), we can obtain the following nonlinear equations in terms of $h_j$ and ${\sigma }_j$

$\begin{array}{l}{f}_1^2={\omega }_1^2\left\{c_1^2{h}_1^2+{\left[\frac{M^2{\omega }_1\left(4{\omega }_1^4-5{\omega }_1^2-3\right)}{4{\left({{\omega }_1}^2-1\right)}^3}{h}_1-2h_1{\sigma }_1-\frac{M{\omega }_1}{{\left(1-{\omega }_1^2\right)}^2}{h}_1\right]}^2\right\}\\ f_2^2={\omega }_2^2\left\{{\left[c_2{h}_2+\frac{\alpha \mu }{\left({\mu }^2+{\omega }_2^2\right)}{h}_2\right]}^2+{\left[\frac{\alpha {\omega }_2}{\left({\mu }^2+{\omega }_2^2\right)}{h}_2-2h_2{\sigma }_2-\frac{h_3^2{h}_2}{2{\omega }_2}+\frac{D_1{h}_1^2{h}_2}{2{\omega }_2}\right]}^2\right\}\\ f_3^2=c_3^2{h}_3^2+{\left[\frac{4\left({{\omega }_2}^2-2\right)h_2^2}{4\left({\omega }_2^2-4\right)}{h}_3-\frac{M}{{\left({{\omega }_1}^2-1\right)}^2}{h}_3+\frac{M^2}{{\left({\omega }_1^2-1\right)}^2}{h}_3-2h_3{\sigma }_3-\frac{D_2{h}_1^2{h}_3}{4}-\frac{\left({\omega }_2^2+8{\omega }_2^2-32\right)h_3^3}{8{\omega }_2^2\left({\omega }_2^2-4\right)}\right]}^2\end{array}$

(62)

The oscillations in the steady-state scenario are regarded as crucial aspects for the stability analysis of the system of ME. We take into account the system's performance in an area that is relatively near the fixed locations to accomplish such a scenario. Consequently, the below expressions of $h_j$ and ${\theta }_j$ are taken into consideration

$h_j=h_{j0}+h_{j1},{\theta }_j={\theta }_{j0}+{\theta }_{j1};\hspace{0.17em}j=\left(\mathrm{1,2,3}\right)$

(63)

where $h_{j0}$ and ${\theta }_{j0}$ are the unperturbed steady-state solutions, while $h_{j1}$ and ${\theta }_{j1}$ are the suggested perturbed ones. The substitution of (63) into (59) yields

$\begin{array}{l}\frac{d{h}_{11}}{d\tau }=\frac{f_1}{2{\omega }_1}\cos \hspace{0.17em}{\theta }_{10}{\theta }_{11}-\frac{c_1{h}_{11}}{2},\\ h_{10}\frac{d{\theta }_{11}}{d\tau }=-\frac{f_1}{2{\omega }_1}\sin \hspace{0.17em}{\theta }_{10}{\theta }_{11}+h_{11}{\sigma }_1+\frac{M{\omega }_1}{2{\left(1-{{\omega }_1}^2\right)}^2}{h}_{11}-\frac{M^2{\omega }_1\left(4{{\omega }_1}^4-5{{\omega }_1}^2-3\right)}{8{\left({{\omega }_1}^2-1\right)}^3}{h}_{11},\\ \frac{d{h}_{21}}{d\tau }=\frac{f_2}{2{\omega }_2}\cos \hspace{0.17em}{\theta }_{20}{\theta }_{21}-\frac{c_2{h}_{21}}{2}-\frac{\alpha \mu }{2\left({\mu }^2+{{\omega }_2}^2\right)}{h}_{21},\\ h_{20}\frac{d{\theta }_{21}}{d\tau }=-\frac{f_2}{2{\omega }_2}\sin \hspace{0.17em}{\theta }_{20}{\theta }_{21}-\frac{\alpha {\omega }_2}{2\left({\mu }^2+{\omega }_2^2\right)}{h}_{21}+h_{21}{\sigma }_2+(h_{30}^2{h}_{21}\\ +2h_{20}{h}_{30}{h}_{31})\left(\frac{1}{4{\omega }_2}-\frac{D_1}{4{\omega }_2}\right),\\ \frac{d{h}_{31}}{d\tau }=\frac{f_3}{2}\cos \hspace{0.17em}{\theta }_{30}{\theta }_{31}-\frac{c_3{h}_{31}}{2},\\ h_{30}\frac{d{\theta }_{31}}{d\tau }=-\frac{f_3}{2}\sin \hspace{0.17em}{\theta }_{30}{\theta }_{31}+\frac{M}{2{\left({{\omega }_1}^2-1\right)}^2}{h}_{31}-\frac{M^2}{2{\left({{\omega }_1}^2-1\right)}^2}{h}_{31}+h_{31}{\sigma }_3\\ +\frac{3\left({{\omega }_2}^4+8{{\omega }_2}^2-32\right)h_{31}{h}_{30}^2}{16{{\omega }_2}^2\left({{\omega }_2}^2-4\right)}+\frac{D_2\left(2h_{10}{h}_{30}{h}_{11}+h_{10}^2{h}_{31}\right)}{8}\\ -\frac{4\left({{\omega }_2}^2-2\right)}{8\left({{\omega }_2}^2-4\right)}\left(2h_{20}{h}_{30}{h}_{21}+h_{20}^2{h}_{31}\right).\end{array}$

(64)

As mentioned above $h_{j1}$ and ${\theta }_{j1}$ refer to the un-identified functions. Then, we can express their solutions as a linear superposition of ${\mathrm{\hslash }}_m{e}^{\lambda \tau }\hspace{0.17em}\left(m=\mathrm{1,2},...,6\right)$ , where ${\mathrm{\hslash }}_m$ and $\lambda$ are constants and eigenvalues, respectively. If $h_{j0}$ and ${\theta }_{j0}$ of the steady-state solution are stable, then the real parts of the roots of the next characteristic equation of (64) will be negative

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

(65)

where ${\mathrm{\Gamma }}_m\hspace{0.17em}\left(m=\mathrm{1,2},...,6\right)$ are functions of $h_{j0}$ and ${\theta }_{j0}$ , see Appendix II

Based on the criteria of Routh–Hurwitz, ⁴⁰ the stability’s conditions of can be written in the following forms

$\begin{array}{l}{\mathrm{\Gamma }}_1>0,{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_2-{\mathrm{\Gamma }}_3>0,{\mathrm{\Gamma }}_3\left({\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_2-{\mathrm{\Gamma }}_3\right)+{\mathrm{\Gamma }}_1\left({\mathrm{\Gamma }}_5-{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_4\right)>0,\\ {\mathrm{\Gamma }}_1\left({\mathrm{\Gamma }}_2{\mathrm{\Gamma }}_3{\mathrm{\Gamma }}_4-{{\mathrm{\Gamma }}_2}^2{\mathrm{\Gamma }}_5+2{\mathrm{\Gamma }}_4{\mathrm{\Gamma }}_5-{\mathrm{\Gamma }}_3{\mathrm{\Gamma }}_6\right)+{\mathrm{\Gamma }}_3\left({\mathrm{\Gamma }}_2{\mathrm{\Gamma }}_5-{\mathrm{\Gamma }}_3{\mathrm{\Gamma }}_4\right)+{{\mathrm{\Gamma }}_1}^2\left({\mathrm{\Gamma }}_2{\mathrm{\Gamma }}_6-{{\mathrm{\Gamma }}_4}^2\right)-{{\mathrm{\Gamma }}_5}^2>0,\\ -{\mathrm{\Gamma }}_5\left({\mathrm{\Gamma }}_4\left(-{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_2{\mathrm{\Gamma }}_3+{{\mathrm{\Gamma }}_3}^2+{{\mathrm{\Gamma }}_1}^2{\mathrm{\Gamma }}_4\right)+{\mathrm{\Gamma }}_1\left({{\mathrm{\Gamma }}_2}^2-2{\mathrm{\Gamma }}_4\right){\mathrm{\Gamma }}_5-{\mathrm{\Gamma }}_2{\mathrm{\Gamma }}_3{\mathrm{\Gamma }}_5+{{\mathrm{\Gamma }}_5}^2\right)+[{{\mathrm{\Gamma }}_3}^3-{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_3\\ \times \left({\mathrm{\Gamma }}_2{\mathrm{\Gamma }}_3+3{\mathrm{\Gamma }}_5\right)+{{\mathrm{\Gamma }}_1}^2\left({\mathrm{\Gamma }}_3{\mathrm{\Gamma }}_4+2{\mathrm{\Gamma }}_2{\mathrm{\Gamma }}_5\right)]{\mathrm{\Gamma }}_6-{{\mathrm{\Gamma }}_1}^3{{\mathrm{\Gamma }}_6}^2>0,\\ {\mathrm{\Gamma }}_6(-{\mathrm{\Gamma }}_5\left({\mathrm{\Gamma }}_4\left(-{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_2{\mathrm{\Gamma }}_3+{{\mathrm{\Gamma }}_3}^2+{{\mathrm{\Gamma }}_1}^2{\mathrm{\Gamma }}_4\right)+{\mathrm{\Gamma }}_1\left({{\mathrm{\Gamma }}_2}^2-2{\mathrm{\Gamma }}_4\right){\mathrm{\Gamma }}_5-{\mathrm{\Gamma }}_2{\mathrm{\Gamma }}_3{\mathrm{\Gamma }}_5+{{\mathrm{\Gamma }}_5}^2\right)+[{{\mathrm{\Gamma }}_3}^3\\ -{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_3\left({\mathrm{\Gamma }}_2{\mathrm{\Gamma }}_3+3{\mathrm{\Gamma }}_5\right)+{{\mathrm{\Gamma }}_1}^2\left({\mathrm{\Gamma }}_3{\mathrm{\Gamma }}_4+2{\mathrm{\Gamma }}_2{\mathrm{\Gamma }}_5\right)]{\mathrm{\Gamma }}_6-{{\mathrm{\Gamma }}_1}^3{{\mathrm{\Gamma }}_6}^2)>0.\end{array}$

(66)

The analysis of stability

The topic of studying the stability of steady-state oscillations is especially important and unique. Therefore, the stability of the considered system is checked in the current section using the approach of nonlinear stability analysis, in which the stability simulation and the influence of the acted parameters are evaluated. It is noted that some important parameters like; the frequencies ${\omega }_j\hspace{0.17em}\left(j=\mathrm{1,2,3}\right)$ , detuning parameters ${\sigma }_j$ , and damping coefficients $c_j$ , play a significant role in the processes of stability. A specific action was used with different system parameters to produce the stability diagrams for the system of ME (59). For various parametrical zones, the altered amplitudes $h_j$ fluctuation with time are depicted, and their characteristics are shown throughout phase plane paths. The possible fixed points of the system are drawn. In Figures 35, 36, 37, 38, 39, and 40, the change of the potential fixed points with various detuning parameters is examined, in which its range is considered to be $-1\le {\sigma }_1\le 1$ . The dotted curve refers to the unstable range, on the contrast the continuous lines denote the stable one. Figures 35 and 36 show the variation of $h_1$ via ${\sigma }_1$ when the used different values of $c_j\hspace{0.17em}\left(j=\mathrm{1,2,3}\right)$ and ${\omega }_k\hspace{0.17em}\left(k=\mathrm{1,2}\right)$ are considered, respectively. One can see from the drawn frequency response curves (FRC) in Figures 35, 37, and 39 that the ranges of the unstable and stable fixed points are, respectively, $-1\le {\sigma }_1\le -0.08$ and $-0.09\le {\sigma }_1\le 1$ when $c_j$ have the considered various values. On the other hand, the fixed point reverses its stability at ${\sigma }_1=-0.09$ . Figures 36, 38, 39, and 40 show the influence of ${\omega }_k$ , on the plotted stability curves in the plane ${\sigma }_1{h}_j$ .It is found that the fixed point reverses its stability at ${\sigma }_1=-0.12$ and ${\sigma }_1=-0.09$ at the considered values of ${\omega }_1$ and ${\omega }_2$ respectively.

Figure 35.

Shows the FRC in the plane ${\sigma }_1{h}_1$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ when: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 36.

Explains the FRC in the plane ${\sigma }_1{h}_1$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ when: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

Figure 37.

Investigates the FRC in the plane ${\sigma }_1{h}_2$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ when: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 38.

Shows the FRC in the plane ${\sigma }_1{h}_2$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ when: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

Figure 39.

Presents stability and instability ranges in the plane ${\sigma }_1{h}_3$ at ${\omega }_1=5.5636$ and ${\omega }_2=13.924$ when: (a) $c_1\left(=\mathrm{0.001,0.005,0.009}\right)$ , (b) $c_2\left(=\mathrm{0.04259,0.0212,0.002129}\right)$ , and (c) $c_3\left(=\mathrm{0.00358,0.0358,0.358}\right)$ .

Figure 40.

Presents stability and instability ranges in the plane ${\sigma }_1{h}_3$ at $c_1=0.001,c_2=0.04259,$ and $c_3=0.00358$ when: (a) ${\omega }_1\left(=\mathrm{5.5636,5.29112,4.73416}\right)$ , and (b) ${\omega }_2\left(=\mathrm{13.924,12.0675,10.2109}\right)$ .

Based on the above simulation, the examined system has one critical fixed point for all values ${\sigma }_1$ , which means that it is under the influence of transcritical bifurcation. This demonstrates that there is no pattern of system behavior that is qualitative when the parameters are changed.

Conclusion

The motion of a 3DOF dynamical system composed of two coupled parts, which are a dynamic part and a HE device, has been studied. The magnet oscillating in the coil of this device determines how the system is constructed. The electromotive force produced an electric potential, which has been used to harvest the energy. The EOM of this system have been derived using LE and solved analytically up to a third-order using the MSM. The obtained novel solutions have been validated through comparison with the numerical solutions of the EOM, which demonstrate a high level of consistency. Three basic cases of external resonance have been simultaneously explored to determine the solvability conditions and the ME. To illustrate their behaviors at any given time, specific diagrams of time histories have been drawn and discussed for the obtained solutions and the current, taking into account certain values of the used parameters. The system of ME has been solved numerically to reveal the amplitudes and adjusted phases for the same data of the considered parameters. The stability criteria of Routh–Hurwitz have been used to discuss the stability and instability zones of this system and show the influence of the acting parameters on the FRC. Therefore, it is possible to convert vibrational motion into electrical energy by including an energy harvester in vibrating dynamical systems. This electrical energy can then be used for a variety of daily tasks, such as remote medical sensing, environmental and structural monitoring, concrete pouring, military applications, and space exploration.

Supplemental Material

Supplemental Material - On the influence of an energy harvesting device on a dynamical system

Supplemental Material for On the influence of an energy harvesting device on a dynamical system by TS Amer, A Arab and AA Galal in Journal of Low Frequency Noise, Vibration & Active Control