Modeling and analysis of a piezoelectric transducer embedded in a nonlinear damped dynamical system

Abstract

This paper focuses on the dynamical analysis of the motion of a new three-degree-of-freedom (DOF) system consisting of two segments that are attached together. External harmonic forces energize this system. The equations of motion (EOM) are derived utilizing Lagrangian equations, and the approximate solutions up to the third order are investigated using the methodology of multiple scales. A comparison between these solutions and numerical ones is constructed to confirm the validity of the analytic solutions. The modulation equations (ME) are acquired from the investigation of the resonance cases and the solvability conditions. The bifurcation diagrams and spectrums of Lyapunov exponent are presented to reveal the different types of the system’s motion and to represent Poincaré maps. The piezoelectric transducer is connected to the dynamical system to convert the vibrational motion into electricity; it is one of the energy harvesting devices which have various applications in our practical life like environmental and structural monitoring, medical remote sensing, military applications, and aerospace. The influences of excitation amplitude, natural frequency, coupling coefficient, damping coefficient, capacitance, and load resistance on the output voltage and power are performed graphically. The steady-state solutions and stability analysis are discussed through the resonance curves.

1 Introduction

The primary energy-producing resources are nonrenewable fossil fuels, but they are swiftly running out and will be exhausted within the next several decades. It is essential to alter the way energy is produced and create clean, renewable energy sources. The most promising form of renewable energy is energy harvesting, which involves gathering wasted ambient energy and transforming it into a more usable form of energy.

One of the used energy harvesting tools, to convert mechanical stress and vibrations into electrical power, is a piezoelectric device [1,2,3]. An innovative piezoelectric energy (PE) harvesting model is examined theoretically and experimentally in [4]. It is based on coupled transverse and interference galloping to enhance the energy's efficiency of the traditional one. In [5], the authors analyzed theoretically the nonlinear dynamic behavior of the energy harvester that is consisting of buckled piezoelectric ribbon layer. The EOM are derived using Lagrange’s equations and the analytical solutions are obtained utilizing the harmonic balance method.

Hamlehdar et al. [6] offered a comprehensive analysis for the methods of capturing energy, from fluid flow, utilizing piezoelectric materials. In [7], the authors suggested a thin-walled PE harvester actuated by fluid pressure to get homogenous solutions, in which the weighted residual approach was used. Two numerical examples are discussed to forecast how the variations and output power would be responded. In [8], a linked approach for modeling numerous PE harvesters depending on vortex-induced vibrations phenomena at random positions was presented and experimentally validated. The Navier–Stokes equations [9] for incompressible fluid flow, piezoelectric equations of Gauss law, and a damper system of mass-spring were coupled to achieve the mathematical formulation.

Wu et al. [10] examined several unique concepts for PE harvesting from natural resources and environmental vibration. The authors provided a detailed summary and a discussion of the approaches and the employed mechanisms to improve the efficiency of piezoelectric power production and energy harvesting. The outcomes of modeling a novel nonlinear multi-stable system for obtaining energy from vibrating mechanical devices were introduced in [11]. In addition, the model was composed of a beam and a system of springs to calculate the potential of a quasi-flat well. They used the pattern of Lyapunov distribution to demonstrate that energy harvesting was diminished in the chaotic motion zones.

A wind energy harvester based on flutter that uses a hybrid of piezoelectric and electromagnetic devices is investigated in [12]. A scientific foundation for the creation of remarkably effective wind energy harvesters is provided in this study. A spring pendulum system was linked to two different energy harvesting devices in [13] each in a separate case, the authors obtained the EOM, and the ME, and they classified all various cases of resonance. The temporal histories of motion and the resonance curves were graphed in their study. The authors in [14], a tri-stable magneto-piezoelastic absorber was used to examine concurrently energy harvesting and vibration isolation. Concepts of Hamilton and Euler–Bernoulli beam, in addition to a magnetic force are used to derive the nonlinear EOM. A one DOF dynamical system of a spring coupled with an energy harvester was investigated in [15]. The outcomes of the optimization revealed that the new framework can output several times more power than the traditional single DOF system. A reduction of vibration and an energy harvesting of a spring pendulum dynamical system is examined in [16]. The pendulum's structure was updated using a separate electromagnetic harvesting device. All zones of stability and instability are examined and discussed.

It is well recognized that mechanical systems, in particular vibrational dynamical problems, e.g. [17,18,19], are among the most significant issues in mechanics due to their wide range of applications, particularly in engineering. In [20], the authors performed a dynamic simulation and bifurcation evaluation for a blade-disk rotor system backed by a rolling bearing. It is found that the slender blades will cause a chaotic behavior at high speeds, and blade stiffness and mass will definitely expand the rotational speed range of quasi-periodic motion as the system's motion sequence gets complicated. In [21], the authors provided a method for identifying Hopf bifurcations [22] in 2D dynamical systems that are characterized by ordinary differential equations (ODE) using an analogue active network. Additionally, they discovered that the experimental technique is capable of detecting global impacts that cannot be expected from a local study, whereas in [23], the authors examined the motion of a damped system attached to an automatic parametric pendulum to constitute a 3DOF dynamical model. They acquired the governing equations and solved them analytically, in which the stability and instability regions are analyzed. The topic of variable-length pendulums was covered in [24]. Using mathematical modeling, dynamical analysis, and innovative computer simulations, an attempt is made to evaluate current developments in this subject in a novel way.

Our purpose of this work is to investigate the vibration analysis of a 3DOF new dynamical system. This system consists of two connected parts; the first part is linked to a piezoelectric transducer to convert the vibrations and stresses into electrical power, and the second part is a nonlinear damped pendulum. The regulating EOM are solved analytically using the methodology of multiple scales till the third order of approximations. Various cases of resonance are assorted, in which three cases of external resonances are investigated at the same time. Therefore, the ME are acquired by utilizing the solvability conditions. To confirm the reliability of the analytic solutions, the numerical ones are compared with them. The bifurcation diagrams are depicted. The spectrum of Lyapunov exponents and Poincaré maps are graphed to illustrate several types of system movements. In the electrical part concerning the PE harvesting, various parameters studied the influence of the output voltage and power graphically. The stability of the model is explored with the help of resonance figures. The value of this work is due to the investigation of the vibration analysis of a 3DOF dynamical system, to acquire beneficial electrical energy for several life applications in the presence of a piezoelectric device.

2 Characterization of the dynamical system

The dynamical system under consideration is composed of two segments linked together. The first portion consists of a mass $m_{1}$, which is coupled to a piezoelectric transducer laden with the load resistance $R_{L}$, and furthermore, it is linked to a linear spring with length $l_{1}$ and stiffness coefficient $k_{1}$. The second portion is a nonlinear spring pendulum with mass $m_{2}$, length $l_{2}$, and the stiffness coefficients $k_{2}$ and $k_{3}$. The pendulum’s system is hanged on the mass $m_{1}$ from point $O$ as seen in Fig. 1. The mentioned dynamical system is energized by the external harmonic forces $F_{1} (t) = F_{1} \,\cos \Omega_{1} t,\,\,F_{2} (t) = F_{2} \,\cos \Omega_{2} t$ and the moment $M(t) = F_{3} \,\cos \Omega_{3} t.$ The frequencies and amplitudes of the applied moment and external forces, respectively, are $\Omega_{n} \,\,(n = 1,2,3)$ and $F_{n}$.

Fig. 1

The dynamical system's depiction

The described dynamical system has three degrees of freedom in addition to three generalized coordinates $X(t),\,\,Z(t)$ (representing the springs’ elongations) and $\theta (t)$ (indicating the rotational angle). It is considered that $C_{n} \,\,(n = 1,2,3)$ are the viscous damping coefficients.

Since Lagrange’s function of any dynamical system has the form $L = T - V$, where $T$ and $V$ are the kinetic and potential energies of the system. According to the investigated system, one can write

$$ \begin{aligned} L & = \frac{1}{2}\{ (m_{1} + m_{2} )X^{{\prime}{2}} + 2m_{2} \,X^{\prime}[Z^{\prime}\sin \theta + (\ell_{2} + Z)\theta ^{\prime}\cos \theta ] + m_{2} [Z^{{\prime}{2}} + (\ell_{2} + Z)^{2} \theta ^{{\prime}{2}} ]\} \\ \,\,\,\,\, - \frac{1}{2}(k_{1} X^{2} + k_{2} Z^{2} + \frac{1}{2}k_{3} Z^{4} ) + m_{2} g(\ell_{2} + Z)\cos \theta . \\ \end{aligned} $$

(1)

Here, $(X^{\prime},Z^{\prime},\theta^{\prime})$ are the system’s generalized velocities that correspond to their generalized coordinates.

The governing system of motion can be obtained according to the following Lagrange’s equations

$$ \frac{d}{dt}\left( {\frac{\partial L}{{\partial q^{\prime}_{j} }}} \right) - \left( {\frac{\partial L}{{\partial q_{j} }}} \right) = Q_{{q_{j} }} ;\,\,\,\,\,\,q_{j} ( = X,Z,\theta ) $$

(2)

where $Q_{{q_{j} }}$ are the applied generalized forces on the examined dynamical system and they have the forms

$$ \begin{aligned} &Q_{X} = F_{1} \,\cos \Omega_{1} t+ F_{2} \,\sin\theta\cos \Omega_{2} t\, - C_{1} \dot{X} - \gamma v, \hfill \\ &Q_{Z} = F_{2} \,\cos \Omega_{2} t\, - C_{2} \dot{Z}, \hfill \\ &Q_{\theta } = F_{3} \,\cos \Omega_{3} t\, - C_{3} \dot{\theta }\,. \hfill \\ \end{aligned} $$

(3)

Here $v$ is the piezoelectric transducer's output voltage and $\gamma$ refers the coupling coefficient.

Now, we are going to propose the following dimensionless parameters to follow the system's inquiry approach

$$ \begin{aligned} X = x\,\ell + \frac{{(m_{1} + m_{2} )g}}{{k_{1} }},\quad \,Z = z\,\ell + \frac{{m_{2} \,g}}{{k_{2} }},\quad \ell_{1} = \ell - \frac{{(m_{1} + m_{2} )g}}{{k_{1} }},\,\,\,\,\,\ell_{2} = \ell - \frac{{m_{2} \,g}}{{k_{2} }},\,\,\,\,\,\, \hfill \\ \omega_{1}^{2} = \frac{{k_{1} \,}}{{(m_{1} + m_{2} )}},\quad \omega_{2}^{2} = \frac{g\,}{\ell }\,,\,\,\,\,\,\,\omega_{3}^{2} = \frac{{k_{2} \,}}{{m_{2} }},\,\,\,\,\,\,\beta = \,\,\frac{{m_{2} \,}}{{(m_{1} + m_{2} )}},\,\,\,\,\alpha = \frac{{k_{3} \,\ell^{2} }}{{\omega_{1}^{2} m_{2} }},\, \hfill \\ \eta = \frac{{m_{2} \,g}}{{k_{2} \,\ell }},\,\,\,\,\,\varpi^{2} = \frac{{\omega_{3}^{2} }}{{\omega_{1}^{2} }},\,\,\,\,\,\,\omega^{2} = \frac{{\omega_{2}^{2} }}{{\omega_{1}^{2} }},\,\,\,\,p_{1} = \frac{{\Omega_{1} }}{{\omega_{1} }},\,\,\,\,\,\,p_{2} = \frac{{\Omega_{2} }}{{\omega_{1} }},\,\,\,\,\,\,p_{2} = \frac{{\Omega_{3} }}{{\omega_{1} }},\,\, \hfill \\ \mu = \frac{\gamma }{{\ell \,(m_{1} + m_{2} )\,\omega_{1}^{2} }},\,\,\,\,\,c_{1} = \frac{{C_{1} }}{{\omega_{1} (m_{1} + m_{2} )}},\,\,\,\,\,\,c_{2} = \frac{{C_{2} }}{{\omega_{1} m_{2} }},\,\,\,\,\,\,c_{3} = \frac{{C_{3} }}{{\ell^{2} \omega_{1} m_{2} }}, \hfill \\ f_{1} = \frac{{F_{1} }}{{\ell \,\omega_{1}^{2} (m_{1} + m_{2} )}},\,\,\,\,f_{2} = \frac{{F_{2} }}{{\ell \,\omega_{1}^{2} m_{2} }},\,\,\,\,\,\,f_{3} = \frac{{F_{3} }}{{\ell^{2} \omega_{1}^{2} m_{2} }},\,\,\,\,\,\,\tau = \omega_{1} t, \hfill \\ \end{aligned} $$

(4)

where the gravity’s acceleration is indicated by $g$, the natural frequencies are symbolized by $\omega_{n} \,\,(n = 1,2,3)$, and $\tau$ is the dimensionless time parameter.

Substituting (1), (3), and (4) into (2) to obtain the following dimensionless system of the EOM

$$ \begin{aligned} \ddot{x} + \beta \,\ddot{z}\,\sin \theta \, + \beta \,(1 + z)\,\ddot{\theta }\cos \theta + 2\beta \,\dot{z}\,\dot{\theta }\cos \theta - \beta \,(1 + z)\,\dot{\theta }^{2} \sin \theta \, \hfill \\ \,\,\,\, + c_{1} \dot{x} + x + \mu \,v + \omega^{2} = f_{1} \cos p_{1} \tau + f_{2}\beta \sin\theta\cos p_{2} \tau\,,\hfill \\ \end{aligned} $$

(5)

$$ \begin{aligned} \ddot{z} + \ddot{x}\sin \theta - (1 + z)\dot{\theta }^{2} + c_{2} \,\dot{z} + \varpi^{2} z + 3\,\eta^{2} \alpha \,z + 3\,\eta \,\alpha \,z^{2} + \alpha \,z^{3} \hfill \\ \,\,\,\, + \omega^{2} (1 - \cos \theta ) = f_{2} \cos p_{2} \tau \,, \hfill \\ \end{aligned} $$

(6)

$$ (1 + z)^{2} \ddot{\theta } + (1 + z)\ddot{x}\cos \theta + 2(1 + z)\,\dot{z}\,\dot{\theta } + c_{3} \,\dot{\theta } + \omega^{2} (1 + z)\sin \theta = f_{3} \cos p_{3} \tau \,, $$

(7)

Moreover, the formula of the piezoelectric transducer's equation can be according to its electrical circuit mechanism, as follows

$$ \dot{v} + \,\frac{v\,}{{c_{p} \,R_{L} \,\omega_{1} }} = \frac{\ell \,\gamma }{{c_{p} \,}}\,\dot{x}, $$

(8)

where $c_{p}$ is the capacitor's equivalent capacity.

3 The approximate methodology

In this part, we are going to solve the abovementioned system of Eqs. (5)–(8) analytically using one of the perturbation techniques. The multiple scale methodology [25] is adopted to obtain the solutions up to the third approximation. Accordingly, we assumed the solutions in terms of the minuscule parameter $\varepsilon$, as follows

$$ \begin{aligned} x(\tau ) = \varepsilon \,\chi (\tau ;\varepsilon ),\quad z(\tau ) = \varepsilon \,\zeta (\tau ;\varepsilon ), \hfill \\ \theta (\tau ) = \varepsilon \,\varphi (\tau ;\varepsilon ),\quad v(\tau ) = \varepsilon \,\upsilon (\tau ;\varepsilon ). \hfill \\ \end{aligned} $$

(9)

The multiple time scales are represented by $\tau_{n} = \varepsilon^{n} \tau \,\,(n = 0,\,1\,,\,2)$, the solutions $\chi ,\,\zeta ,\,\varphi$, and $\upsilon$ can be written as follows

$$ \begin{aligned} \chi = \sum\limits_{k = 1}^{3} {\varepsilon^{k} \chi_{k} (\tau_{0} ,\tau_{1} ,\tau_{2} )} + O(\varepsilon^{4} ), \hfill \\ \zeta = \sum\limits_{k = 1}^{3} {\varepsilon^{k} \zeta_{k} (\tau_{0} ,\tau_{1} ,\tau_{2} )} + O(\varepsilon^{4} ), \hfill \\ \varphi = \sum\limits_{k = 1}^{3} {\varepsilon^{k} \varphi_{k} (\tau_{0} ,\tau_{1} ,\tau_{2} )} + O(\varepsilon^{4} ), \hfill \\ \upsilon = \sum\limits_{k = 1}^{3} {\varepsilon^{k} \upsilon_{k} (\tau_{0} ,\tau_{1} ,\tau_{2} )} + O(\varepsilon^{4} ). \hfill \\ \end{aligned} $$

(10)

Considering the narrowness of the following variables and parameters

$$ \begin{aligned} &\mu = \varepsilon^{2} \,\tilde{\mu },\,\,\,\,\,\,\,\gamma = \varepsilon^{2} \,\tilde{\gamma },\,\,\,\,\,\,\,\alpha = \varepsilon^{2} \,\tilde{\alpha },\,\,\,\,\,\,\,\beta = \varepsilon \,\tilde{\beta }, \hfill \\ c_{n} = \,\varepsilon^{2} \,\tilde{c}_{n} ,\,\,\,\,\,\,\,f_{n} = \,\varepsilon^{3} \,\tilde{f}_{n} ;\,\,\,\,(n = 1,2,3),\,\,\,\,\,\, \hfill \\ &R_{L} = \frac{{\tilde{R}_{L} }}{{\varepsilon^{2} }},\,\,\,\,\,\,\,c_{p} \, = \varepsilon^{2} \tilde{c}_{p} . \hfill \\ \end{aligned} $$

(11)

The next operators are employed to transform the derivatives regarding $\tau$, to the corresponding ones regarding the time scales $\tau_{n}$

$$ \begin{aligned} &\frac{d}{d\tau } = \frac{\partial }{{\partial \tau_{0} }} + \varepsilon \frac{\partial }{{\partial \tau_{1} }} + \varepsilon^{2} \frac{\partial }{{\partial \tau_{2} }}, \hfill \\ &\frac{{d^{2} }}{{d\tau^{2} }} = \frac{{\partial^{2} }}{{\partial \tau_{0}^{2} }} + 2\varepsilon \frac{{\partial^{2} }}{{\partial \tau_{0} \partial \tau_{1} }} + \varepsilon^{2} \left( {\frac{{\partial^{2} }}{{\partial \tau_{1}^{2} }} + 2\frac{{\partial^{2} }}{{\partial \tau_{0} \partial \tau_{2} }}} \right) + O(\varepsilon^{3} ). \hfill \\ \end{aligned} $$

(12)

To acquire the following the groups of partial differential equations (PDE) in powers of $\varepsilon$, substituting (9)–(12) into the EOM (5)–(8), and neglecting the terms of $O(\varepsilon^{3} )$ due to their smallness.

• (i) Order of $(\varepsilon )$

  $$ \frac{{\partial^{2} \chi_{1} }}{{\partial \tau_{0}^{2} }} + \chi_{1} = 0, $$

  (13)
  $$ \frac{{\partial^{2} \zeta_{1} }}{{\partial \tau_{0}^{2} }} + \varpi^{2} \zeta_{1} = 0, $$

  (14)
  $$ \frac{{\partial^{2} \varphi_{1} }}{{\partial \tau_{0}^{2} }} + \omega^{2} \varphi_{1} + \frac{{\partial^{2} \chi_{1} }}{{\partial \tau_{0}^{2} }} = 0, $$

  (15)
  $$ \frac{{\partial \upsilon_{1} }}{{\partial \tau_{0} }} + \frac{{\upsilon_{1} }}{{\omega_{1} \,\tilde{c}_{p} \,\tilde{R}_{L} }} = \frac{{\ell \,\tilde{\gamma }}}{{\tilde{c}_{p} }}\frac{{\partial \chi_{1} }}{{\partial \tau_{0} }}. $$

  (16)

• (ii) Order of $(\varepsilon^{2} )$

  $$ \frac{{\partial^{2} \chi_{2} }}{{\partial \tau_{0}^{2} }} + \chi_{2} = - \tilde{\beta }\frac{{\partial^{2} \varphi_{1} }}{{\partial \tau_{0}^{2} }} - 2\frac{{\partial^{2} \chi_{1} }}{{\partial \tau_{0} \partial \tau_{1} }}, $$

  (17)
  $$ \frac{{\partial^{2} \zeta_{2} }}{{\partial \tau_{0}^{2} }} + \varpi^{2} \zeta_{2} = - \frac{1}{2}\omega^{2} \varphi_{1}^{2} + \left( {\frac{{\partial \varphi_{1} }}{{\partial \tau_{0} }}} \right)^{2} - 2\frac{{\partial^{2} \zeta_{1} }}{{\partial \tau_{0} \partial \tau_{1} }} - \varphi_{1} \frac{{\partial^{2} \chi_{1} }}{{\partial \tau_{0}^{2} }}, $$

  (18)
  $$ \begin{aligned} \frac{{\partial^{2} \varphi_{2} }}{{\partial \tau_{0}^{2} }} + \omega^{2} \varphi_{2} + \frac{{\partial^{2} \chi_{2} }}{{\partial \tau_{0}^{2} }} &= - \zeta_{1} \left( {\omega^{2} \varphi_{1} + 2\frac{{\partial^{2} \varphi_{1} }}{{\partial \tau_{0}^{2} }} + \frac{{\partial^{2} \chi_{1} }}{{\partial \tau_{0}^{2} }}} \right) \hfill \\ &\quad - 2\left( {\frac{{\partial \zeta_{1} }}{{\partial \tau_{0} }}\frac{{\partial \varphi_{1} }}{{\partial \tau_{0} }} + \frac{{\partial^{2} \varphi_{1} }}{{\partial \tau_{0} \partial \tau_{1} }} + \frac{{\partial^{2} \chi_{1} }}{{\partial \tau_{0} \partial \tau_{1} }}} \right), \hfill \\ \end{aligned} $$

  (19)
  $$ \frac{{\partial \upsilon_{2} }}{{\partial \tau_{0} }} + \frac{{\upsilon_{2} }}{{\omega_{1} \,\tilde{c}_{p} \,\tilde{R}_{L} }} = \frac{{\ell \,\tilde{\gamma }}}{{\tilde{c}_{p} }}(\frac{{\partial \chi_{1} }}{{\partial \tau_{1} }} + \frac{{\partial \chi_{2} }}{{\partial \tau_{0} }}) - \frac{{\partial \upsilon_{1} }}{{\partial \tau_{1} }}. $$

  (20)

• (iii) Order of $(\varepsilon^{3} )$

  $$ \begin{aligned} \frac{{\partial^{2} \chi_{3} }}{{\partial \tau_{0}^{2} }} + \chi_{3} &= \frac{1}{2}\tilde{f}_{1} e^{{i\,p_{1} \,\tau_{0} }} - \tilde{\mu }\,\upsilon_{1} \\ &\quad -\, \frac{{\partial^{2} \chi_{1} }}{{\partial \tau_{1}^{2} }} - \tilde{c}_{1} \frac{{\partial \chi_{1} }}{{\partial \tau_{0} }} - 2\left( {\frac{{\partial^{2} \chi_{1} }}{{\partial \tau_{0} \partial \tau_{2} }} + \frac{{\partial^{2} \chi_{2} }}{{\partial \tau_{0} \partial \tau_{1} }}} \right) \hfill \\ &\quad -\, \tilde{\beta }\left( {2\frac{{\partial \zeta_{1} }}{{\partial \tau_{0} }}\frac{{\partial \varphi_{1} }}{{\partial \tau_{0} }} + \varphi_{1} \frac{{\partial^{2} \zeta_{1} }}{{\partial \tau_{0}^{2} }} + \zeta_{1} \frac{{\partial^{2} \varphi_{1} }}{{\partial \tau_{0}^{2} }} + 2\frac{{\partial^{2} \varphi_{1} }}{{\partial \tau_{0} \partial \tau_{1} }} + \frac{{\partial^{2} \varphi_{2} }}{{\partial \tau_{0}^{2} }}} \right), \hfill \\ \end{aligned} $$

  (21)
  $$ \begin{aligned} \frac{{\partial^{2} \zeta_{3} }}{{\partial \tau_{0}^{2} }} + \varpi^{2} \zeta_{3} &= \frac{1}{2}\tilde{f}_{2} e^{{i\,p_{2} \,\tau_{0} }} - 3\eta^{2} \tilde{\alpha }\,\zeta_{1} - \omega^{2} \varphi_{1} \varphi_{2} - \frac{{\partial^{2} \zeta_{1} }}{{\partial \tau_{1}^{2} }} - \tilde{c}_{2} \frac{{\partial \zeta_{1} }}{{\partial \tau_{0} }} \hfill \\ &\quad +\, 2\frac{{\partial \varphi_{1} }}{{\partial \tau_{0} }}\left( {\frac{1}{2}\zeta_{1} \frac{{\partial \varphi_{1} }}{{\partial \tau_{0} }} + \frac{{\partial \varphi_{2} }}{{\partial \tau_{0} }} + \frac{{\partial \varphi_{1} }}{{\partial \tau_{1} }}} \right) + 2\left( {\frac{{\partial^{2} \zeta_{1} }}{{\partial \tau_{0} \partial \tau_{2} }} + \frac{{\partial^{2} \zeta_{2} }}{{\partial \tau_{0} \partial \tau_{1} }}} \right) \hfill \\ &\quad -\, \varphi_{1} \left( {2\frac{{\partial^{2} \chi_{1} }}{{\partial \tau_{0} \partial \tau_{1} }} + \frac{{\partial^{2} \chi_{2} }}{{\partial \tau_{0}^{2} }}} \right) - \varphi_{2} \frac{{\partial^{2} \chi_{1} }}{{\partial \tau_{0}^{2} }}\,, \hfill \\ \end{aligned} $$

  (22)
  $$ \begin{aligned} \frac{{\partial^{2} \varphi_{3} }}{{\partial \tau_{0}^{2} }} + \omega^{2} \varphi_{3} &= \frac{1}{2}\tilde{f}_{3} e^{{i\,p_{3} \,\tau_{0} }} - \omega^{2} \left( {\zeta_{2} \varphi_{1} + \zeta_{1} \varphi_{2} - \frac{{\varphi_{1}^{3} }}{6}} \right) - \frac{{\partial^{2} \varphi_{1} }}{{\partial \tau_{1}^{2} }} - \frac{{\partial^{2} \chi_{1} }}{{\partial \tau_{1}^{2} }} - \tilde{c}_{3} \frac{{\partial \varphi_{1} }}{{\partial \tau_{0} }} \hfill \\ &\quad - 2\frac{{\partial \zeta_{1} }}{{\partial \tau_{0} }}\left( {\frac{{\partial \varphi_{1} }}{{\partial \tau_{1} }} + \frac{{\partial \varphi_{2} }}{{\partial \tau_{0} }}} \right) - 2\frac{{\partial \varphi_{1} }}{{\partial \tau_{0} }}\left( {\frac{{\partial \zeta_{1} }}{{\partial \tau_{1} }} + \zeta_{1} \frac{{\partial \zeta_{1} }}{{\partial \tau_{0} }} + \frac{{\partial \zeta_{2} }}{{\partial \tau_{0} }}} \right) - 2\left( {\frac{{\partial^{2} \varphi_{1} }}{{\partial \tau_{0} \partial \tau_{2} }} + \frac{{\partial^{2} \chi_{1} }}{{\partial \tau_{0} \partial \tau_{2} }}} \right) \hfill \\&\quad - 2\left\{ {\zeta_{1} \left[ {2\frac{{\partial^{2} \varphi_{1} }}{{\partial \tau_{0} \partial \tau_{1} }} + \frac{{\partial^{2} \chi_{1} }}{{\partial \tau_{0} \partial \tau_{1} }} + \frac{{\partial^{2} \varphi_{2} }}{{\partial \tau_{0}^{2} }} + \frac{1}{2}\frac{{\partial^{2} \chi_{2} }}{{\partial \tau_{0}^{2} }}} \right] + \frac{{\partial^{2} \varphi_{2} }}{{\partial \tau_{0} \partial \tau_{1} }} + \frac{{\partial^{2} \chi_{2} }}{{\partial \tau_{0} \partial \tau_{1} }}} \right\} \hfill \\&\quad - \frac{{\partial^{2} \varphi_{1} }}{{\partial \tau_{0}^{2} }}(\zeta_{1}^{2} + 2\zeta_{2} ) - \frac{{\partial^{2} \chi_{1} }}{{\partial \tau_{0}^{2} }}\left( {\zeta_{2} - \frac{{\varphi_{1}^{2} }}{2}} \right), \hfill \\ \end{aligned} $$

  (23)
  $$ \frac{{\partial \upsilon_{3} }}{{\partial \tau_{0} }} + \frac{{\upsilon_{3} }}{{\omega_{1} \,\tilde{c}_{p} \,\tilde{R}_{L} }} = \frac{{\ell \,\tilde{\gamma }}}{{\tilde{c}_{p} }}\left( {\frac{{\partial \chi_{1} }}{{\partial \tau_{2} }} + \frac{{\partial \chi_{2} }}{{\partial \tau_{1} }} + \frac{{\partial \chi_{3} }}{{\partial \tau_{0} }}} \right) - \frac{{\partial \upsilon_{1} }}{{\partial \tau_{2} }} - \frac{{\partial \upsilon_{2} }}{{\partial \tau_{1} }}. $$

  (24)

The PDE of these groups can be solved sequentially. Therefore, we begin with the general solutions of the first group of Eqs. (13)–(16) as follows

$$ \chi_{1} \, = A_{1} \,e^{{i\tau_{0} }} + CC, $$

(25)

$$ \zeta_{1} \, = A_{2} \,e^{{i\,\varpi \,\tau_{0} }} + CC, $$

(26)

$$ \varphi_{1} \, = A_{3} \,e^{{i\,\omega \,\tau_{0} }} + \frac{{A_{1} \,e^{{i\tau_{0} }} }}{{1 - \omega^{2} }} + CC, $$

(27)

$$ \upsilon_{1} \, = \frac{{\ell \,\omega_{1} \tilde{\gamma }\,\tilde{R}_{L} A_{1} \,e^{{i\,\,\tau_{0} }} }}{{(\omega_{1} \,\tilde{c}_{p} \,\tilde{R}_{L} - i\,)}} + CC, $$

(28)

Here, $CC$ denotes the complex conjugates of the foregoing terms, and $A_{n} \,(n = 1,\,2,\,3)$ and $\overline{A}_{n}$ are unknown complex functions of the scales $(\tau_{1} ,\tau_{2} )$ and their conjugates, respectively.

Substituting the solutions (25)–(28) into the second group of PDE (17)–(20), and then eliminating the produced secular terms, to get the following conditions

$$ \frac{{\partial A_{n} }}{{\partial \tau_{1} }} = 0,\,\quad \,n = 1,\,2,\,3. $$

(29)

As a corollary of these conditions, the solutions of Eqs. (17)–(20) are

$$ \chi_{2} \, = \frac{{\omega^{2} \tilde{\beta }A_{3} \,e^{{i\,\omega \,\tau_{0} }} }}{{(1 - \omega^{2} )}} + CC, $$

(30)

$$ \begin{aligned} \zeta_{2} \, & = \frac{{(4 - 3\omega^{2} )A_{1} \overline{A}_{1} }}{{\varpi^{2} (\omega^{2} - 1)^{2} }} + \frac{{\omega^{2} A_{3} \overline{A}_{3} }}{{\varpi^{2} }}\, - \frac{{3\omega^{2} A_{1}^{2} \,e^{{2i\,\tau_{0} }} }}{{2(\varpi^{2} - 4)(\omega^{2} - 1)^{2} }} - \frac{{3\omega^{2} A_{3}^{2} \,e^{{2i\,\omega \,\tau_{0} }} }}{{2(\varpi^{2} - 4\omega^{2} )}} \\ &\quad + \frac{{[2\omega (\omega + 1) - 1]\omega^{2} A_{1} A_{3} \,e^{{i(\,\omega + 1)\,\tau_{0} }} }}{{(\omega^{2} - 1)[\varpi^{2} - (\omega + 1)^{2} ]}} + \frac{{[2\omega (\omega - 1) - 1]\omega^{2} A_{3} \overline{A}_{1} \,e^{{i(\,\omega - 1)\,\tau_{0} }} }}{{(\omega^{2} - 1)[\varpi^{2} - (\omega - 1)^{2} ]}} + CC, \\ \end{aligned} $$

(31)

$$ \begin{aligned} \varphi_{2} & = \frac{{(2\omega^{2} - 2\varpi - 3)A_{1} A_{2} e^{{i(\varpi + 1)\tau_{0} }} }}{{(\omega^{2} - 1)[\omega^{2} - (\varpi + 1)^{2} ]}} + \frac{{(2\omega^{2} + 2\varpi - 3)A_{2} \overline{A}_{1} e^{{i(\varpi - 1)\tau_{0} }} }}{{(\omega^{2} - 1)[\omega^{2} - (\varpi - 1)^{2} ]}} \\ &\quad - \frac{{\omega (2\varpi + \omega )A_{2} A_{3} e^{{i(\varpi + \omega )\tau_{0} }} }}{\varpi (\varpi + 2\omega )} + \frac{{\omega (2\varpi - \omega )A_{2} \overline{A}_{3} e^{{i(\varpi - \omega )\tau_{0} }} }}{\varpi (\varpi - 2\omega )} + CC, \\ \end{aligned} $$

(32)

$$ \upsilon_{2} = - \frac{{\ell \,\omega^{3} \omega_{1} \tilde{R}_{L} \tilde{\gamma }A_{3} \,e^{{i\,\omega \tau_{0} }} }}{{(\omega \,\omega_{1} \tilde{c}_{p} \tilde{R}_{L} - i\,)(\omega^{2} - 1)}} + CC, $$

(33)

The requirements of removing the secular terms from the third group of PDE (21)–(24) can be obtained after inserting the solutions (25)–(28) and (30)–(33) into these equations. Therefore, one obtains

$$ i\,\tilde{c}_{1} A_{1} + 2i\frac{{\partial A_{1} }}{{\partial \tau_{2} }} + D_{1} A_{1} = 0, $$

(34)

$$ i\,\tilde{c}_{2} \varpi A_{2} \,\, + 3\,\eta^{2} \tilde{\alpha }A_{2} + 2\,i\,\varpi \frac{{\partial A_{2} }}{{\partial \tau_{2} }} + D_{2} A_{2} A_{3} \overline{A}_{3} + D_{3} A_{1} A_{2} \overline{A}_{1} = 0, $$

(35)

$$ i\,\tilde{c}_{3} \omega A_{3} + 2\,i\,\omega \frac{{\partial A_{3} }}{{\partial \tau_{2} }} - D_{2} A_{3} (D_{4} A_{3} \overline{A}_{3} - A_{2} \overline{A}_{2} ) + D_{5} A_{1} A_{3} \overline{A}_{1} = 0. $$

(36)

Ultimately, the solutions of Eqs. (21)–(24) can be found in the forms

$$ \begin{aligned} \chi_{3} &= \frac{{\tilde{f}_{1} e^{{i\,p_{1} \,\tau_{0} }} }}{{2(1 - p_{1}^{2} )}} + \frac{{(\varpi + 1)^{2} (\omega^{2} + \varpi^{2} - 2)\tilde{\beta }A_{1} A_{2} e^{{i\,(\varpi + 1)\tau_{0} }} }}{{\varpi (\varpi + 2)(\varpi - \omega + 1)(\varpi + \omega + 1)(\omega^{2} - 1)}} \hfill \\ &\quad + \frac{{(\varpi - \omega )(\varpi + \omega )^{3} \tilde{\beta }A_{2} A_{3} e^{{i\,(\varpi + \omega )\tau_{0} }} }}{{\varpi (\varpi + 2\omega )[1 - (\varpi + \omega )^{2} ]}} + \frac{{(\varpi - 1)^{2} (\omega^{2} + \varpi^{2} - 2)\tilde{\beta }A_{2} \overline{A}_{1} e^{{i\,(\varpi - 1)\tau_{0} }} }}{{\varpi (\varpi - 2)(\omega^{2} - 1)[(\varpi - 1)^{2} - \omega^{2} ]}} \hfill \\ &\quad + \frac{{(\varpi + \omega )(\varpi - \omega )^{3} \tilde{\beta }A_{2} \overline{A}_{3} e^{{i\,(\varpi - \omega )\tau_{0} }} }}{{\varpi (\varpi - 2\omega )[1 - (\varpi - \omega )^{2} ]}} + CC, \hfill \\ \end{aligned} $$

(37)

$$ \zeta_{3} = \zeta_{31} + \zeta_{32} + \zeta_{33} , $$

(38)

$$ \varphi_{3} = \varphi_{31} + \varphi_{32} + \varphi_{33} , $$

(39)

where $\zeta_{3n}$ and $\varphi_{3n} ,\,(n = 1,2,3)$ are defined as in the appendix (I).

$$ \begin{aligned} \upsilon_{3} & = \frac{{\ell \,p_{1} \omega_{1} \tilde{R}_{L} \,\tilde{\gamma }\,\tilde{f}_{1} e^{{i\,p_{1} \tau_{0} }} }}{{2(1 - p_{1}^{2} )(p_{1} \omega_{1} \tilde{c}_{p} \tilde{R}_{L} - i\,)}} + \frac{{\ell \,\omega_{1} \tilde{R}_{L} \,\tilde{\gamma }\,\tilde{\beta }(\omega^{2} + \varpi^{2} - 2)}}{{\varpi (\omega^{2} - 1)}} \\ &\quad \times \left\{ {\frac{{(\varpi + 1)^{3} A_{1} A_{2} e^{{i(\varpi + 1)\tau_{0} }} }}{{(\varpi + 2)[(\varpi + 1)^{2} - \omega^{2} ][\omega_{1} \tilde{c}_{p} \tilde{R}_{L} (\varpi + 1) - i\,]}}} \right. \\ &\quad \left. { + \frac{{(\varpi - 1)^{3} \overline{A}_{1} A_{2} e^{{i(\varpi - 1)\tau_{0} }} }}{{(\varpi - 2)[(\varpi - 1)^{2} - \omega^{2} ][\omega_{1} \tilde{c}_{p} \tilde{R}_{L} (\varpi - 1) - i\,]}}} \right\} \\ &\quad - \frac{{\ell \,\omega_{1} \tilde{R}_{L} \,\tilde{\gamma }\,\tilde{\beta }}}{\varpi }\left\{ {\frac{{(\varpi - \omega )(\varpi + \omega )^{4} A_{2} A_{3} e^{{i(\varpi + \omega )\tau_{0} }} }}{{(\varpi + 2\omega )[(\varpi + \omega )^{2} - 1][\omega_{1} \tilde{c}_{p} \tilde{R}_{L} (\varpi + \omega ) - i\,]}}} \right. \\ &\quad \left. { + \frac{{(\varpi + \omega )(\varpi - \omega )^{4} A_{2} \overline{A}_{3} e^{{i(\varpi - \omega )\tau_{0} }} }}{{(\varpi - 2\omega )[(\varpi - \omega )^{2} - 1][\omega_{1} \tilde{c}_{p} \tilde{R}_{L} (\varpi - \omega ) - i\,]}}} \right\} + CC. \\ \end{aligned} $$

(40)

4 Resonance investigation

In this section, different aspects of resonance [26] will be identified which affect the behavior of the system. Moreover, the solvability conditions that generate the modulation equations will be determined. Meanwhile, the frequency response curves will be used to highlight the zones of stability and instability according to the system’s physical parameters.

The resonance conditions are fulfilled based on the denominators of the second and third orders of solutions (30)–(33) and (37)–(40) when they have zero values. Here, we have two divisions of resonance.

• a. Internal resonance takes place, once the following terms are met

  $$ \begin{aligned} \varpi \approx 0,\quad \,\varpi \approx - 1,\quad \varpi \approx \pm 2,\quad \varpi \approx \pm \omega ,\quad \varpi \approx \pm 2\omega ,\quad \varpi \approx \pm (\omega + 1), \hfill \\ \varpi \approx \pm (\omega - 1),\quad \omega \approx \pm 1,\quad \omega \approx \pm (2\varpi + 1),\quad \omega \approx \pm (2\varpi - 1),\quad \omega \approx \pm 3,\quad \omega \approx \pm {1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-0pt} 3}. \hfill \\ \end{aligned} $$

• b. External primary resonance can be achieved when $p_{1} \approx 1,$ $p_{2} \approx \varpi ,$ and $p_{3} \approx \omega$.

Now, we will explore the cases of external primary resonance at the same time. Then, we introduce the detuning parameter $\sigma_{n} \,\,(n = 1,2,3)$ according to

$$ p_{1} = 1 + \sigma_{1} ,\quad p_{2} = \varpi + \sigma_{2} ,\quad p_{3} = \omega + \sigma_{3} ;\,\quad \sigma_{n} = \varepsilon^{2} \tilde{\sigma }_{n} . $$

(41)

Substituting (41) into the third-order solutions (37)–(40), and making use of the conditions (34)–(36), the solvability conditions become

$$ \begin{aligned} & \frac{1}{2}\tilde{f}_{1} \,e^{{i\,\tau_{2} \tilde{\sigma }_{1} }} - i\,\tilde{c}_{1} A_{1} - 2i\frac{{\partial A_{1} }}{{\partial \tau_{2} }} - D_{1} A_{1} = 0, \\ & \frac{1}{2}\tilde{f}_{2} \,e^{{i\,\tau_{2} \tilde{\sigma }_{2} }} - i\,\tilde{c}_{2} \varpi A_{2} \, - 3\,\eta^{2} \tilde{\alpha }A_{2} \\ &\quad - 2\,i\,\varpi \frac{{\partial A_{2} }}{{\partial \tau_{2} }} - D_{2} A_{2} A_{3} \overline{A}_{3} - D_{3} A_{1} A_{2} \overline{A}_{1} = 0, \\ &\frac{1}{2}\tilde{f}_{3} \,e^{{i\,\tau_{2} \tilde{\sigma }_{3} }} - i\,\tilde{c}_{3} \omega A_{3} - 2\,i\,\omega \frac{{\partial A_{3} }}{{\partial \tau_{2} }} \\ &\quad + D_{2} A_{3} [D_{4} A_{3} \overline{A}_{3} - A_{2} \overline{A}_{2} ] - D_{5} A_{1} A_{3} \overline{A}_{1} = 0, \\ \end{aligned} $$

(42)

Let us consider $\tilde{a}_{n}$ and $\tilde{\psi }_{n}$ to express the amplitudes and phases of the solutions $\chi ,\zeta ,$ and $\varphi$. Therefore, we can specify the functions $A_{n}$ in their polar forms, as follows

$$ A_{n} = \frac{{\tilde{a}_{n} (\tau_{2} )}}{2}e^{{i\tilde{\psi }_{n} \tau_{2} }} ;\quad \,a_{n} = \varepsilon \tilde{a}_{n} \quad (n = 1,2,3). $$

(43)

Henceforth, their first derivatives can be stated as follows

$$ \frac{{\partial A_{j} }}{\partial \tau } = \varepsilon^{2} \frac{{\partial A_{j} }}{{\partial \tau_{2} }}.\,\,\,\,\, $$

(44)

For simplicity, the above solvability conditions of PDE (42) can be transformed into the ordinary form. To gain this objective, the following modified phases can be utilized

$$ \begin{aligned} \theta_{1} (\tau_{1} ,\tau_{2} ) = \tau_{2} \,\tilde{\sigma }_{1} - \psi_{1} (\tau_{2} ), \hfill \\ \theta_{2} (\tau_{1} ,\tau_{2} ) = \tau_{2} \,\tilde{\sigma }_{2} - \psi_{2} (\tau_{2} ), \hfill \\ \theta_{3} (\tau_{1} ,\tau_{2} ) = \tau_{2} \,\tilde{\sigma }_{3} - \psi_{3} (\tau_{2} ). \hfill \\ \end{aligned} $$

(45)

Substituting (43)–(45) into (42) and isolating the real and imaginary portions to obtain the following ME

$$ \begin{aligned} a_{1} \frac{{d\theta_{1} }}{d\tau } = \frac{1}{2}[f_{1} \cos \theta_{1} - a_{1} (D_{6} - 2\sigma_{1} )], \hfill \\ \frac{{da_{1} }}{d\tau } = \frac{1}{2}[f_{1} \sin \theta_{1} - a_{1} (D_{7} + c_{1} )], \hfill \\ a_{2} \frac{{d\theta_{2} }}{d\tau } = \frac{1}{2\varpi }\left\{ {f_{2} \cos \theta_{2} - a_{2} \left[ {3\alpha \,\eta^{2} - 2\varpi \sigma_{2} + \frac{1}{4}(D_{2} a_{3}^{2} + D_{3} a_{1}^{2} )} \right]} \right\}, \hfill \\ \frac{{da_{2} }}{d\tau } = \frac{1}{2\varpi }(f_{2} \sin \theta_{2} - \varpi \,a_{2} c_{2} ), \hfill \\ a_{3} \frac{{d\theta_{3} }}{d\tau } = \frac{1}{2\omega }\left\{ {f_{3} \cos \theta_{3} + \frac{{\omega \,a_{3} }}{2}\left[ {4\sigma_{3} + \frac{1}{2}D_{2} a_{3} (D_{4} a_{3}^{2} - a_{2}^{2} )\,\, - D_{5} a_{1}^{2} } \right]} \right\}, \hfill \\ \frac{{da_{3} }}{d\tau } = \frac{1}{2\omega }(f_{3} \sin \theta_{3} - \omega \,a_{3} c_{3} ), \hfill \\ \end{aligned} $$

(46)

where the terms $D_{n} \,\,(n = 1,2,...,7)$ are given in Appendix (I).

The modulation equations system (46) may be numerically solved through the planes $a_{n} \theta_{n}$ and according to the arbitrary initial conditions $a_{1} (0) = 0.01,$ $a_{2} (0) = 0.1,$ $a_{3} (0) = 0.08,$ and $\theta_{1} (0) = \theta_{2} (0) = \theta_{3} (0) = 0$.

A comparison between the numerical solutions and the approximate ones is graphed in Fig. 2. The portions of this figure are drawn according to the following data of the system’s parameter [27]

$$ \begin{aligned} f_{1} = f_{2} = f_{3} = 0.01 \times 10^{ - 6} ,\quad c_{1} = 0.002,\quad \,c_{2} = 0.003,\quad c_{3} = 0.0012, \hfill \\ \varpi = 0.01,\quad \,\omega = 0.15,\quad \mu = 0.000003,\quad \,\eta = 0.0001,\,\quad \beta = 0.01 \hfill \\ \alpha = 0.01,\,\quad \gamma = 0.00003,\quad \,\sigma_{1} = \sigma_{2} = - 0.003,\,\quad \sigma_{3} = 0.003. \hfill \\ \end{aligned} $$

Fig. 2

The numerical and approximate solutions comparison

The inference from the figure’s portions reveals high consistency among these solutions; meanwhile, it confirms an excellent accuracy of the approximate solutions. The subsequent parameters were considered for the comparison as shown in Fig. 2.

5 Chaotic motion

In the mentioned dynamical system, tiny adjustments to the system's parameters can have an impact on the system's behavior. Hence, the system offered three different types of motion as we will describe below. The bifurcation diagrams of the variables $X$ and $Z$, and the spectrum of Lyapunov exponents (LEs of $\lambda_{n} ,\,n = 1,2,...,6$) versus the excitation amplitude $f_{1}$ were simulated to display the kinds of these motions [28]. Figure 3a, b shows the bifurcation diagrams of $X$ and $Z$, and Fig. 3c displays the spectrum LEs that describe different types of motion. According to the below three figures, we have three ranges of the excitation amplitude $f_{1}$, every range refers to different kind of system’s motion. The first range at $f_{1} \in [0,0.04065]$, we noticed in this range of values, the spectrum of LEs lines are approximately under the axis $\lambda = 0$ or equal to zero, which indicated that the motion is in a quasi-periodic sequence and it is seen also in the bifurcation diagram of $X,\,\,Z$ for the same range of $f_{1}$ especially in Fig. 3a, the system exhibits nearly as a line. The second range is founded at $f_{1} \in (0.04065,0.23]$. It is noticed that, there is more than one LEs spectrum that has a positive value like the blue and red lines. Hence, the system behavior is a hyper-chaotic motion as it is obviously observed in the bifurcation diagrams for the aforementioned second range. Finally, we can observe from the third range of $f_{1} \in (0.23,0.5]$ that only one spectrum has a positive value. Therefore, the system is in the chaotic phase as seen at the same time in the bifurcation diagrams of $X$ and $Z$.

Fig. 3

The bifurcation diagrams of $X,\,\,Z$, and their spectrum of LEs

The phase portraits and Poincaré maps are plotted in Fig. 4, where the blue curves denote the phase portraits and the red dots symbolize the Poincaré maps. The figures were depicted for different values of the excitation amplitude $f_{1}$ to show the various three motions of the system.

Fig. 4

The phase portraits and Poincare maps of $X,\,\,Z$

In Fig. 4a, we consider a value $f_{1} = 0.0005$ that lies in the first range of the excitation amplitude. The red dots pattern of the Poincaré map shows approximately a closed curve for both variables $X$ and $Z$, which leads to the previous conclusion of the quasi-periodic movement of the system.

The red dots at $f_{1} = 0.24$ as in Fig. 4b start to diverge randomly because of the chaotic state which appeared in the third range of values of $f_{1}$. In the same way, in Fig. 4c, the system at $f_{1} = 0.11$ behaves as the hyper-chaotic motion as in the second range of $f_{1}$. Therefore, the red dots are in a very messy distribution unlike in the previous two figures.

6 The piezoelectric transducer

In this portion of the study, we pursued to take the advantage of producing electrical energy by integrating a piezoelectric device with our mentioned dynamical model, which is one of the energy harvesting devices. The piezoelectric device consists of dielectric materials, which can be polarized due to the mechanical stress that comes from the vibrations of the dynamical model in our case when these materials are polarized; an electric field is generated. This process is called the piezoelectric effect. So, electrical energy from mechanical energy using the merge between the piezoelectric device and the dynamical model is produced. The electrical energy is generated by the energy harvesting device which may be used for many different purposes, including environmental monitoring like habitat monitoring (light, temperature, and humidity), structural monitoring, and medical remote sensing such as emergency medical response monitoring, pacemaker, defibrillator, military applications, and aerospace.

The piezoelectric transducer is merged with our dynamical model specifically with the mass $m_{1}$ and loaded with the resistance $R_{L}$, in which the influence of various system parameters on the energy production is checked. In Fig. 5, the effect of different values of the damping coefficient $c_{1} ( = 0.002,\,0.005,\,0.008)$, capacitance $c_{p} ( = 0.001,\,0.004,\,0.007)$, and the load resistance $R_{L} ( = 500,\,700,\,900)$ on the time histories of output voltage and power is drawn. We logically concluded that the output voltage and the power decrease with the increasing of $c_{1}$ and $c_{p}$. It is observed that the output voltage rises as the load resistance increases in accordance with Ohm's Law [29], and in contrast, as $R_{L}$ grows, the output power falls.

Fig. 5

The time histories of the output voltage, and power at different values of damping coefficient $c_{1}$, capacitance $c_{p}$, and load resistance $R_{L}$

Simultaneously, the effect of other parameters which, increase the productivity of electrical energy is studied, as shown in Fig. 6. Moreover, the impact of various values of the excitation amplitude $f_{1} ( = 0.0005,\,0.0007,\,0.0009)$, natural frequency $\omega_{1} ( = 2.5,\,3.5,\,4.5)$, and the coupling coefficient $\gamma ( = 0.0003,\,0.0009,\,0.0015)$ is examined. The rising of the mentioned parameters yields, the output voltage and power growth as time passed.

Fig. 6

The behavior of $v$ and $P$ versus $\tau$, when varying of excitation amplitude $f_{1}$, natural frequency $\omega_{1}$, and coupling coefficient $\gamma$

7 Steady-state solutions and stability analysis

This section of the research covers the steady-state solutions and the stability analysis of the dynamical model according to the resonance response curves. When transient processes are eliminated as a result of the system's dampening, steady-state vibrations are apparent [30]. We presumed that the amplitudes’ and phases’ time derivatives are equal to zero in modulation Eq. (46). Therefore, these equations are transformed from a system of ODE to a system of algebraic equations as follows

$$ \begin{aligned} f_{1} \cos \theta_{1} - a_{1} (D_{6} - 2\sigma_{1} ) = 0, \hfill \\ f_{1} \sin \theta_{1} - a_{1} (D_{7} + c_{1} ) = 0, \hfill \\ f_{2} \cos \theta_{2} - a_{2} \left[ {3\alpha \,\eta^{2} - 2\varpi \sigma_{2} + \frac{1}{4}(D_{2} a_{3}^{2} + D_{3} a_{1}^{2} )} \right] = 0, \hfill \\ f_{2} \sin \theta_{2} - \varpi \,a_{2} c_{2} = 0, \hfill \\ f_{3} \cos \theta_{3} + \frac{{\omega a_{3} }}{2}\left[ {4\sigma_{3} + \frac{1}{2}D_{2} a_{3} (D_{4} a_{3}^{2} - a_{2}^{2} )\, - D_{5} a_{1}^{2} } \right] = 0, \hfill \\ f_{3} \sin \theta_{3} - \omega \,a_{3} c_{3} = 0. \hfill \\ \end{aligned} $$

(47)

The following set of nonlinear algebraic equations is a relation between the amplitudes $a_{n}$ and the detuning parameters $\sigma_{n}$ after getting rid of the modified phases $\theta_{n}$ from Eq. (47)

$$ \begin{aligned} f_{1}^{2} = a_{1}^{2} \{ [D_{6} - 2\sigma_{1} ]^{2} + [D_{7} + c_{1} ]^{2} \} , \hfill \\ f_{2}^{2} = a_{2}^{2} \left\{ {\varpi^{2} c_{2}^{2} + \left[ {3\alpha \,\eta^{2} - 2\varpi \sigma_{2} + \frac{1}{4}(D_{2} a_{3}^{2} + D_{3} a_{1}^{2} )} \right]^{2} } \right\}, \hfill \\ f_{3}^{2} = \omega^{2} a_{3}^{2} \left\{ {\left[ {2\sigma_{3} + \frac{1}{4}D_{2} a_{3} (D_{4} a_{3}^{2} - a_{2}^{2} ) - D_{5} a_{1}^{2} } \right]^{2} + c_{3}^{2} } \right\}. \hfill \\ \end{aligned} $$

(48)

The stability of steady-state vibrations is a unique and significant feature of the topic of this study. Consequently, we will examine the system's behavior in the immediate vicinity of the fixed locations to understand this issue. The fixed point analogous to the steady-state solution is asymptotically stable in the perspective of Lyapunov if and only if the real components of all eigenvalues are negative [31].

We will illustrate the stability of the fixed points through the resonance response curves which can be plotted using the numerical solutions of Eq. (48). In Figs. 7 and 8, we noticed that all fixed points for the ranges of $\sigma_{1}$ and $\sigma_{2}$ are stable. In contrast and as shown in Fig. 9, the influence of $\sigma_{3}$ on the stability of the fixed points is obvious, in which the fixed points of some areas lost their stability and turn into unstable fixed points. It must be noted that the stable fixed points are represented with blue points while the red ones are used to refer to the unstable ones.

Fig. 7

Resonance curves of the amplitudes $a_{n} ,\,n = 1,2,3$ versus $\sigma_{1} ,$ and $\sigma_{2} = - 0.002,\,\sigma_{3} = - 0.001$

Fig. 8

The amplitudes $a_{n}$ vs $\sigma_{2} ,$ at $\sigma_{1} = 0.001,\,\sigma_{3} = - 0.001$

Fig. 9

Representation of the resonance curves of the amplitudes $a_{n}$ vs $\sigma_{3} ,$ at $\sigma_{1} = 0.001,\,\sigma_{2} = - 0.002$

8 Conclusion

The nonlinear motion of a 3DOf dynamical system consisting of two linked components: The first one is connected to a PE device, while the other is a nonlinear damped pendulum, has been investigated. The governing equations are obtained using Lagrange’s equations and are solved analytically using the multiple scale methodology till the third-order approximation. These solutions are confirmed numerically and supported with graphs. All the obtained external resonance cases are examined and the modulation equations are achieved. The bifurcation diagrams, Poincaré maps, and the spectrum of Lyapunov exponents are depicted and explored that the system acted in three different movements: a quasi-periodic, hyper-chaotic, and chaotic motion. It is noted that we can get high output voltage and power by controlling the values of excitation amplitude, natural frequency, coupling coefficient, capacitance, damping coefficient, and load resistance. The influences of these parameters are represented graphically to be more comprehensible. Moreover, the stability of the system using the resonance response curves is drawn and analyzed to clarify the stable and unstable fixed points. The significance of the investigated dynamical system is due to its uses in diverse applications including medical remote sensing like emergency medical response monitoring, and environmental applications; for example, temperature and humidity monitoring.

Appendix I

$$ \begin{aligned} \zeta_{31} & = \frac{{\tilde{f}_{2} e^{{i\,p_{2} \tau_{0} }} }}{{2(\varpi^{2} - p_{2}^{2} )}} - \frac{{\tilde{\beta }\,\omega^{4} A_{3}^{2} e^{{2i\omega \tau_{0} }} }}{{(\varpi^{2} - 4\omega^{2} )^{2} (\omega^{2} - 1)}} - \frac{{2\tilde{\beta }\,\omega^{4} A_{3} \overline{A}_{3} }}{{\varpi^{2} (\omega^{2} - 1)}} + \frac{{\tilde{\beta }\,\omega^{4} A_{3} \overline{A}_{1} e^{{i(\omega - 1)\tau_{0} }} }}{{(\omega^{2} - 1)^{2} [\varpi^{2} - (\omega^{2} - 1)^{2} ]}} \\ \,\,\, + \frac{{\tilde{\beta }\,\omega^{4} A_{3} \overline{A}_{1} e^{{i(\omega + 1)\tau_{0} }} }}{{(\omega^{2} - 1)^{2} [\varpi^{2} - (\omega^{2} + 1)^{2} ]}} + \frac{{[4\omega^{4} - 5\omega^{2} - 3\varpi (\varpi + 2) - 2]A_{2} A_{1}^{2} e^{{i(\varpi + 2)\tau_{0} }} }}{{4(\varpi + 1)(\varpi - \omega + 1)(\varpi + \omega + 1)(\omega^{2} - 1)^{2} }} \\ \,\,\,\, + \frac{{[4\omega^{4} - 5\omega^{2} - 3\varpi (\varpi - 2) - 2]A_{2} \overline{A}_{1}^{2} e^{{i(\varpi - 2)\tau_{0} }} }}{{4[(\varpi - 1)^{3} - \omega^{2} (\varpi - 1)](\omega^{2} - 1)^{2} }} - \frac{{3\omega (\varpi + \omega )A_{2} A_{3}^{2} e^{{i(\varpi + 2\omega )\tau_{0} }} }}{4\varpi (\varpi + 2\omega )} \\ &\quad + \frac{{3\omega (\varpi - \omega )A_{2} \overline{A}_{3}^{2} e^{{i(\varpi - 2\omega )\tau_{0} }} }}{4\varpi (\varpi - 2\omega )}, \\ \end{aligned} $$

$$ \begin{aligned} \zeta_{32} & = \omega \,A_{1} A_{2} A_{3} e^{{i(\varpi + \omega + 1)\tau_{0} }} [6\varpi^{4} + 12\varpi^{3} (\omega + 1) - \omega (2\omega (\omega + 1) - 1)(\omega^{2} - 1) \\ &\quad + \varpi^{2} (2\omega (13 - 4\omega (\omega - 1)) + 4) - 2\varpi (\omega + 1)(\omega \,(\omega (4\omega - 1) - 7) \\ &\quad + 1)]/[\varpi (\varpi - \omega + 1)(\omega - 1)(\omega + 1)^{2} (\varpi + \omega + 1)(2\varpi + \omega + 1)(\varpi + 2\omega )] \\&\quad + \omega \,A_{2} A_{3} \overline{A}_{1} e^{{i(\varpi + \omega - 1)\tau_{0} }} [6\varpi^{4} + 12\varpi^{3} (\omega - 1)\, + \omega (2\omega (\omega - 1) - 1)(\omega^{2} - 1) \\ &\quad + \varpi^{2} (2\omega (4\omega (\omega + 1) + 13) + 4) + 2\varpi (\omega - 1)(\omega (4\omega^{2} + \omega - 7) - 1)]/[\varpi (\varpi \\ &\quad + 2\omega )(\omega^{2} - 1)(\varpi^{2} - (\varpi + \omega - 1)^{2} )(\omega^{2} - (\varpi - 1)^{2} )], \\ \end{aligned} $$

$$ \begin{aligned} \zeta_{33} & = \,\omega \,A_{2} \overline{A}_{1} \overline{A}_{3} e^{{i(\varpi - \omega - 1)\tau_{0} }} [ - 6\varpi^{4} + 12\varpi^{3} (\omega + 1)\, + \omega (2\omega (\omega + 1) - 1)(\omega^{2} - 1)\, \\ &\quad + \varpi^{2} (2\omega (4\omega (\omega - 1) - 13) - 4) - 2\varpi (\omega + 1)(\omega (\omega (4\omega - 1) - 7) + 1)]\,/[\varpi (\varpi \\ &\quad - 2\omega )(\varpi - \omega - 1)(2\varpi - \omega - 1)(\omega - 1)(\omega + 1)^{2} (\varpi + \omega - 1)] \\&\quad + \omega \,A_{1} A_{2} \overline{A}_{3} e^{{i(\varpi - \omega + 1)\tau_{0} }} [2\varpi (\varpi + 1)(3\varpi (\varpi + 1) - 1) + \omega - 2\varpi \,\omega (2\varpi + 3)(3\varpi \\ &\quad , + 2) + 2\omega^{2} (4\varpi (\varpi + 2) + 1) + \omega^{3} (8\varpi^{2} + 6\varpi - 3) - 2\omega^{4} (4\varpi + 1) + 2\omega^{5} ] \\ &a\quad /[\varpi (\varpi - 2\omega )(\varpi - \omega + 1)(\varpi + \omega + 1)(\omega^{2} - 1)(\varpi^{2} - (\varpi - \omega + 1)^{2} )] + CC, \\ \end{aligned} $$

$$ \begin{aligned} \varphi_{31} &= \frac{{\tilde{f}_{3} e^{{i\,p_{3} \,\tau_{0} }} }}{{2(\omega^{2} - p_{3}^{2} )}} + \frac{{p_{1}^{2} \tilde{f}_{1} e^{{i\,p_{1} \,\tau_{0} }} }}{{2(1 - p_{1}^{2} )(\omega^{2} - p_{1}^{2} )}} - \frac{{[\varpi^{2} (4\omega^{2} - 3) + 18\omega^{4} - 79\omega^{2} + 12]A_{1}^{3} e^{{3\,i\,\tau_{0} }} }}{{6(\omega^{2} - 9)(\varpi^{2} - 4)(\omega^{2} - 1)^{3} }} \hfill \\ &\quad - \frac{{(\varpi^{2} + \omega^{2} - 2)\tilde{\beta }A_{2} }}{{\varpi (\omega^{2} - 1)}}\{ \frac{{(\varpi + 1)^{4} A_{1} e^{{i(\varpi + 1)\tau_{0} }} }}{{(\varpi + 2)[(\varpi + 1)^{2} - \omega^{2} ]^{2} }} + \frac{{(\varpi - 1)^{4} \overline{A}_{1} e^{{i(\varpi - 1)\tau_{0} }} }}{{(\varpi - 2)[(\varpi - 1)^{2} - \omega^{2} ]^{2} }}\} \hfill \\ &\quad - \tilde{\beta }A_{2} A_{3} e^{{i(\varpi + \omega )\tau_{0} }} [4\varpi^{2} \omega^{4} + 4\varpi^{3} \omega^{5} + \omega^{6} - \omega^{8} + \varpi (\omega^{2} - 1)(\varpi^{5} + 4\varpi^{4} \omega - 2\omega^{5} ) \hfill \\ &\quad + \varpi^{4} \omega^{2} (6\omega^{2} - 5)]/\{ \varpi (\omega^{2} - 1)(\varpi + 2\omega )[\omega^{2} - (\varpi + \omega )^{2} ][(\varpi + \omega )^{2} - 1]\} \hfill \\ &\quad - \tilde{\beta }A_{2} \overline{A}_{3} e^{{i(\varpi - \omega )\tau_{0} }} \,[4\varpi^{2} \omega^{4} - 4\varpi^{3} \omega^{5} + \omega^{6} - \omega^{8} + \varpi (\omega^{2} - 1)(\varpi^{5} - 4\varpi^{4} \omega + 2\omega^{5} ) \hfill \\ &\quad + \varpi^{4} \omega^{2} (6\omega^{2} - 5)]/\{ \varpi (\omega^{2} - 1)(\varpi - 2\omega )[\omega^{2} - (\varpi - \omega )^{2} ][(\varpi - \omega )^{2} - 1]\} \hfill \\ &\quad - \frac{{(\varpi^{2} - 49\omega^{2} )A_{3}^{3} e^{{3\,i\,\,\omega \tau_{0} }} }}{{48(\varpi^{2} - 4\omega^{2} )}}, \hfill \\ \end{aligned} $$

$$ \begin{aligned} \varphi_{32} & = - A_{1}^{2} A_{3} e^{{i(\omega + 2)\tau_{0} }} \{ \varpi^{4} (3\omega^{2} - 2) + \varpi^{2} [\omega (\omega (2\omega (\omega - 9) - 45) - 12) + 20] + \omega [\omega (\omega (\omega \\ &\quad \times (3\omega (\omega + 6) + 7) + 36) + 132) + 48] - 48\} /\{ 2(\varpi^{2} - 4)(\omega^{2} - 1)^{2} [(\omega + 1)^{2} - \varpi^{2} ][\omega^{2} \\ &\quad - (\omega + 2)^{2} ]\} + \overline{A}_{1}^{2} A_{3} e^{{i(\omega - 2)\tau_{0} }} \{ \varpi^{4} (3\omega^{2} - 2) + \varpi^{2} [\omega (\omega (2\omega (\omega + 9) - 45) + 12) + 20] \\ &\quad + \omega [\omega (\omega (\omega (3\omega (\omega - 6) + 7) - 36) + 132) - 48] - 48\} /\{ 2(\varpi^{2} - 4)(\omega^{2} - 1)^{2} [\omega^{2} - (\omega \\ &\quad - 2)^{2} (\varpi - \omega + 1)(\varpi + \omega - 1)]\} - A_{1} A_{3}^{2} e^{{i(2\omega + 1)\tau_{0} }} \{ \varpi^{4} (2\omega^{2} - 1) + \varpi^{2} (1 + 6\omega - 4\omega^{2} (4\omega \\ &\quad + 1)(\omega + 2) + \omega^{2} [\omega (\omega (50\omega^{2} + 96\omega + 39) + 6) + 5]\} /\{ 2(\omega^{2} - 1)(4\omega^{2} - \varpi^{2} )[(\omega + 1)^{2} \\ &\quad - \varpi^{2} ][\omega^{2} - (2\omega + 1)^{2} ]\} - \overline{A}_{1} A_{3}^{2} e^{{i(2\omega - 1)\tau_{0} }} \{ \varpi^{4} (2\omega^{2} - 1) + \varpi^{2} (1 - 6\omega - 4\omega^{2} (4\omega \, - 1)(\omega \\ &\quad - 2) + \omega^{2} [\omega (\omega (50\omega^{2} - 96\omega + 39) - 6) + 5]\} /\{ 2(\omega^{2} - 1)(4\omega^{2} - \varpi^{2} )[(\omega - 1)^{2} - \varpi^{2} ] \\ &\quad \times [\omega^{2} - (2\omega - 1)^{2} ]\} , \\ \end{aligned} $$

$$ \begin{aligned} \varphi_{33} = - \frac{{[2\omega^{4} + 18\varpi + 19\varpi^{2} + 6\varpi^{3} - 2\omega^{2} (4\varpi^{2} + 6\varpi + 3) + 5]A_{1} A_{2}^{2} e^{{i(2\varpi + 1)\tau_{0} }} }}{{(\omega^{2} - 1)[\omega^{2} - (\varpi + 1)^{2} ][\omega^{2} - (2\varpi + 1)^{2} ]}} \hfill \\ &\quad - \frac{{[2\omega^{4} - 18\varpi + 19\varpi^{2} - 6\varpi^{3} - 2\omega^{2} (4\varpi^{2} - 6\varpi + 3) + 5]\overline{A}_{1} A_{2}^{2} e^{{i(2\varpi - 1)\tau_{0} }} }}{{(\omega^{2} - 1)[\omega^{2} - (\varpi - 1)^{2} ][\omega^{2} - (2\varpi - 1)^{2} ]}} \hfill \\ &\quad - \frac{{\omega (\varpi + \omega )(2\varpi + \omega )(3\varpi + \omega )A_{2}^{2} A_{3} e^{{i(2\varpi + \omega )\tau_{0} }} }}{{\varpi (\varpi + 2\omega )[\omega^{2} - (2\varpi + \omega )^{2} ]}} + \frac{{\omega (2\varpi - \omega )(3\varpi^{2} - 4\varpi \omega + \omega^{2} )A_{2}^{2} \overline{A}_{3} e^{{i(2\varpi - \omega )\tau_{0} }} }}{{\varpi (\varpi - 2\omega )[\omega^{2} - (2\varpi - \omega )^{2} ]}} \hfill \\ &\quad + \frac{{A_{1} e^{{i\,\tau_{0} }} }}{{2\varpi^{2} (\omega^{2} - 1)^{4} }}\{ \frac{{[3(\varpi^{4} - 12\varpi^{2} + 32) - \omega^{2} (4\varpi^{4} - 47\varpi^{2} + 136) - 6\omega^{4} (3\varpi^{2} - 8)]A_{1} \overline{A}_{1} }}{{(\varpi^{2} - 4)}} \hfill \\ &\quad + ([2(\omega^{2} - 1)^{2} (2\varpi^{2} (5 - 8\varpi^{2} + 3\varpi^{4} + \omega^{2} (3\varpi^{2} - 11) + 2\omega^{4} (\varpi^{2} + 4) - 2\omega^{6} ]A_{2} \overline{A}_{2} \hfill \\ &\quad /[\omega^{4} + (\varpi^{2} - 1)^{2} - 2\omega^{2} (\varpi^{2} + 1)]) + ([\varpi^{2} (\varpi^{2} - 1)^{2} - \omega^{2} (2\varpi^{6} + 7\varpi^{4} - 16\varpi^{2} - 3) \hfill \\ &\quad + \omega^{4} (2\varpi^{4} - 7\varpi^{2} + 8) - \omega^{6} (2\varpi^{2} + 7) + 2\omega^{8} ]A_{3} \overline{A}_{3} )\} + CC, \hfill \\ \end{aligned} $$

$$ \begin{aligned} D_{1} = \frac{{\ell \,\omega_{1} \tilde{R}_{L} \,\tilde{\gamma }\,\tilde{\mu }}}{{\omega_{1} \,\tilde{c}_{p} \tilde{R}_{L} - i}},\quad D_{2} = \frac{{6\omega^{2} (\varpi^{2} - \omega^{2} )}}{{(\varpi^{2} - 4\omega^{2} )}}, \hfill \\ D_{3} = \frac{{2[8 + 3\varpi^{4} - 19\omega^{2} + 15\omega^{4} - 4\omega^{6} + \varpi^{2} (4\omega^{4} + 2\omega^{2} - 9)]}}{{(\omega^{2} - 1)^{2} [\varpi^{4} + (\omega^{2} - 1)^{2} - 2\varpi^{2} (\omega^{2} + 1)]}}, \hfill \\ D_{4} = \frac{{(\varpi^{2} + 8\omega^{2} )}}{{12\varpi^{2} }}, \hfill \\ D_{5} = \omega^{2} (\omega^{2} - 1)^{2} (3\omega^{2} - 4) - \varpi^{6} (3\omega^{2} - 2) + \varpi^{4} (\omega^{4} - 2\omega^{2} - 6) \hfill \\ &\quad - \varpi^{2} (\omega^{6} + \,14\omega^{4} - 19\omega^{2} - 4)]/\varpi^{2} (\omega^{2} - 1)^{2} [(\varpi^{2} - 1)^{2} - 2\omega^{2} (\varpi^{2} - 1) + \omega^{4} , \hfill \\ D_{6} = \frac{{\ell \,\omega_{1}^{2} R_{L}^{2} \,\gamma \,\mu \,}}{{(\omega_{1}^{2} c_{p}^{2} R_{L}^{2} + 1)}},\,\quad D_{7} = \frac{{\ell \,\omega_{1} R_{L} \,\gamma \,\mu \,}}{{(\omega_{1}^{2} c_{p}^{2} R_{L}^{2} + 1)}}. \hfill \\ \end{aligned} $$