Perspective Chapter: Maximizing Energy Collection from Nonlinear Harvesting System through Optimization and Control Techniques with Induced Time Delays

Abstract

This investigation explores the utilization of periodic and Quasi-Periodic (QP) vibrations for Energy Harvesting (EH) in a delayed nonlinear oscillator system. The system consists of a delayed Duffing-van der Pol oscillator and a delayed piezoelectric coupling medium, with a focus on the occurrence of delay parametric resonance. This occurs when the frequency of the delay width modulation in the mechanical component is close to twice the oscillator’s natural frequency. The double-step stress system is used to approximate the QP delay width, which is then harnessed for power generation. The results suggest that only periodic vibration can be used for energy harvesting with small, unmodulated delay width in the mechanical component. However, larger values result in instability of periodic vibration, and only QP vibration can be used for better performance. The impact of delay in the electrical circuit on the performance of periodic and quasi-periodic vibration-based EH is also analyzed, and it is found that suitable values of delay and delay frequency in the electrical component can improve energy harvesting performance within a specific range of parameters. These findings are supported by numerical simulations.

1. Introduction

The limitations of linear attachments in Energy Harvesting (EH) systems can be overcome by using nonlinear stiffness in the mechanical part. This leads to improved EH performance, either in the case of monostable harvesters with hardening characteristics [1, 2] or bistable ones [3]. However, the utilization of nonlinear attachments may result in instability and unexpected oscillations near the limits of the stable frequency response range, as stated in Ref. [4]. Under self-excitation and linear stiffness, energy harvesting in the range of large-amplitude oscillations can be achieved through Limit-Cycle (LC) oscillations in the harvester. However, these LC oscillations can become unstable through a secondary Hopf bifurcation, leading to Quasi-Periodic (QP) vibrations over a wide frequency range, which are not near resonance [5, 6].

However, in certain Energy Harvesting systems subjected to aerodynamic and base excitations, it has been observed that Quasi-Periodic (QP) vibrations lead to a significant decrease in harvested power beyond the flutter speed. To optimize energy extraction, it is necessary to avoid these QP vibrations [7, 8].

Despite this, research has shown that when there is a time delay, the vibrations from QP result in large, broad-spectrum amplitude swings [9]. Inspired by these results, in [10] a van der Pol-type energy harvester system was created that uses delay amplitude modulation to harness the large-amplitude, broad-spectrum vibrations from QP. This resulted in efficient energy extraction. The study in [10] introduced time delay in the mechanical component, while [11] explored the case where delay was present in both the mechanical and electrical subsystems. The findings in [11] revealed that the highest output power of the harvester does not always coincide with the maximum system response amplitude, depending on the values of time delay frequency and amplitude.

It was emphasized that in a delayed Duffing-type oscillator system, which is subjected to harmonic excitation and connected to a piezoelectric circuit harvester, the large-amplitude QP vibrations generated can be utilized for energy extraction across a broad range of excitation frequencies that are away from the resonance. This approach allows for the avoidance of hysteresis and instability that may occur near resonance [12].

Recently, in [13] it has been shown that a delayed Duffing harvester device, which is connected to an electrical circuit through a piezoelectric mechanism and has a modulated delay amplitude, is capable of producing large-amplitude Quasi-Periodic (QP) vibrations due to the modulation of the delay amplitude. These vibrations can be effectively utilized for energy extraction, away from resonance, with desirable performance.

Building upon previous studies that have explored the use of QP vibrations for energy harvesting in systems with time delay in the mechanical component [10, 12, 13], and those that have included time delay in both the mechanical and electrical components as a control mechanism [11], this study aims to broaden the scope of the investigation by examining the energy harvesting performance of a delayed Duffing-van der Pol oscillator coupled with a delayed piezoelectric harvester circuit. The focus is on the impact of time delay in both the mechanical and electrical components on the system’s energy harvesting performance.

Additionally, the time delay in the piezoelectric coupling can be utilized as a means to control and optimize the harvester’s output power. On the other hand, the time delay in the mechanical subsystem of the harvester is a common occurrence in milling and turning operations, resulting in inherent delay in the position [14, 15]. As a result, it is not accounted for as an extra source of energy for the harvester.

In the following section, the harvester system will be introduced. The multiple scales method will be used to determine the periodic response and harvested power near a delay parametric resonance. The effect of the delay parameters on the periodic response and resultant power output will also be explored. The second-step multiple scales method will be utilized in Section 3 to determine the QP response and evaluate the harvested power. The impact of different harvester system parameters on the energy harvesting performance will be examined, with a summary of the results presented in the conclusion section.

2. Explanation of the model and energy harvesting utilizing periodic vibrations

The energy harvesting device under consideration in this study consists of a piezoelectric device connecting a delayed Duffing-van der Pol oscillator to an electrical circuit, as shown in the schematic illustration in Figure 1. The mechanics and electronics of the harvester are impacted by a lag in feedback, leading to a system equation that can be expressed dimensionlessly

Figure 1.

Schematic description of the EH system.

$$\ddot{\mathrm{x}}\left(\mathrm{t}\right)+\mathrm{\delta }\dot{\mathrm{x}}\left(\mathrm{t}\right)+\mathrm{\lambda }\dot{\mathrm{x}}\left(\mathrm{t}\right)\mathrm{x}{\left(\mathrm{t}\right)}^2+\mathrm{x}\left(\mathrm{t}\right)+\mathrm{\gamma x}{\left(\mathrm{t}\right)}^3-\mathrm{\chi v}\left(\mathrm{t}\right)=\mathrm{\alpha }\left(\mathrm{t}\right)\mathrm{x}\left(\mathrm{t}-{\mathrm{\tau }}_1\right)$$E1

$$\dot{\mathrm{v}}\left(\mathrm{t}\right)+\mathrm{\beta v}\left(\mathrm{t}\right)+\mathrm{\kappa }\dot{\mathrm{x}}\left(\mathrm{t}\right)={\mathrm{\alpha }}_3\mathrm{v}\left(\mathrm{t}-{\mathrm{\tau }}_2\right)$$E2

The relative displacement of the rigid mass m is represented by$\mathrm{x}\left(\mathrm{t}\right)$, the voltage across the load resistance is $\mathrm{v}\left(\mathrm{t}\right)$, and various system parameters are denoted by $\mathrm{\delta }$, $\mathrm{\lambda }$, $\mathrm{\gamma }$,$\mathrm{\chi }$, $\mathrm{\kappa }$, $\mathrm{\beta }$, ${\mathrm{\alpha }}_3$, ${\mathrm{\tau }}_2$, $\mathrm{\alpha }\left(\mathrm{t}\right)$, and ${\mathrm{\tau }}_1$. These parameters correspond to the mechanical damping ratio, stiffness parameter, piezoelectric coupling in the mechanical attachment, piezoelectric coupling in the electrical circuit, reciprocal of the time constant of the electrical circuit, feedback gain and time delay in the electric circuit, and feedback gain and time delay in the mechanical component. It is worth mentioning that the delay in the mechanical component originates from milling and turning operations, as described in Refs. [14, 15]. In contrast, the delay feedback in the electrical circuit is incorporated to enhance the harvester’s output power, as indicated in Ref. [11]. The delay control with modulated delay amplitude in the position is represented by $\mathrm{\alpha }\left(\mathrm{t}\right)$

$$\mathrm{\alpha }\left(\mathrm{t}\right)={\mathrm{\alpha }}_1+{\mathrm{\alpha }}_{2}cos\left(\mathrm{\omega t}\right)$$E3

The unmodulated delay amplitude is denoted by ${\mathrm{\alpha }}_1$, while the amplitude and frequency of modulation are represented by ${\mathrm{\alpha }}_2$ and $\mathrm{\omega }$ respectively. The use of modulated delay amplitude has been extensively utilized to enhance energy harvesting performance [10, 11, 13].

This study expands upon previous investigations, where the effect of time delay on the energy harvesting performance was examined. For instance, the impact of a uniform delay in both the mechanical and electrical components without a nonlinear stiffness ($\mathrm{\gamma }$= 0) was analyzed in [11], and the case of a linear damper and unmodulated delay was studied in [12]. Here, the focus is on examining the optimization of the EH performance through the introduction of different time delays and amplitudes in the mechanical and electrical components, as described by Eqs. (1) and (2). In this study, we examine the system’s behavior at the delay parametric resonance by applying the resonance condition $1=\frac{{\mathrm{\omega }}^2}{4}+\mathrm{\sigma }$ where $\mathrm{\sigma }$ is a detuning parameter. The multiple scales method [16] is applied by introducing a bookkeeping parameter, $\mathrm{ϵ}$, to allow for scaling of parameters as $\mathrm{\delta }=\mathrm{ϵ}\stackrel{\sim }{\mathrm{\delta }},\mathrm{\lambda }=\mathrm{ϵ}\stackrel{\sim }{\mathrm{\lambda }},\mathrm{\gamma }=\mathrm{ϵ}\stackrel{\sim }{\mathrm{\gamma }},\mathrm{\chi }=\mathrm{ϵ}\stackrel{\sim }{\mathrm{\chi }},{\mathrm{\alpha }}_1=\mathrm{ϵ}{\stackrel{\sim }{\mathrm{\alpha }}}_1,{\mathrm{\alpha }}_2={\stackrel{\sim }{\mathrm{\alpha }}}_2$, $\mathrm{\sigma }=\mathrm{ϵ}\stackrel{\sim }{\mathrm{\sigma }}$. Consequently, Eqs. (1) and (2) are expressed as

$$\ddot{\mathrm{x}}\left(\mathrm{t}\right)+\frac{{\mathrm{\omega }}^2}{4}\mathrm{x}=\mathrm{ϵ}\left[-\stackrel{\sim }{\mathrm{\delta }}\dot{\mathrm{x}}\left(\mathrm{t}\right)-\stackrel{\sim }{\mathrm{\lambda }}\dot{\mathrm{x}}\left(\mathrm{t}\right)\mathrm{x}{\left(\mathrm{t}\right)}^2-\stackrel{\sim }{\mathrm{\gamma }}\mathrm{x}{\left(\mathrm{t}\right)}^3+\stackrel{\sim }{\mathrm{\chi }}\mathrm{v}\left(\mathrm{t}\right)-\stackrel{\sim }{\mathrm{\sigma }}\mathrm{x}\left(\mathrm{t}\right)+\left({\stackrel{\sim }{\mathrm{\alpha }}}_1+{\stackrel{\sim }{\mathrm{\alpha }}}_{2}cos\left(\mathrm{\omega t}\right)\right)\mathrm{x}\left(\mathrm{t}-{\mathrm{\tau }}_1\right)\right]$$E4

$$\dot{\mathrm{v}}\left(\mathrm{t}\right)+\mathrm{\beta v}\left(\mathrm{t}\right)+\mathrm{\kappa }\dot{\mathrm{x}}\left(\mathrm{t}\right)={\mathrm{\alpha }}_3\mathrm{v}\left(\mathrm{t}-{\mathrm{\tau }}_2\right)$$E5

A solution to Eqs. (4) and (5) can be sought in the form

$$\mathrm{x}\left(\mathrm{t}\right)={\mathrm{x}}_0\left({\mathrm{T}}_0,{\mathrm{T}}_1\right)+\mathrm{ϵ}{\mathrm{x}}_1\left({\mathrm{T}}_0,{\mathrm{T}}_1\right)+\mathrm{O}\left({\mathrm{ϵ}}^2\right)$$E6

$$\mathrm{v}\left(\mathrm{t}\right)={\mathrm{v}}_0\left({\mathrm{T}}_0,{\mathrm{T}}_1\right)+\mathrm{ϵ}{\mathrm{v}}_1\left({\mathrm{T}}_0,{\mathrm{T}}_1\right)+\mathrm{O}\left({\mathrm{ϵ}}^2\right)$$E7

where ${\mathrm{T}}_0=\mathrm{t}$, and ${\mathrm{T}}_1=\mathrm{ϵt}$. In terms of the variables ${\mathrm{T}}_{\mathrm{i}}$, the time derivatives become $\frac{\mathrm{d}}{\mathrm{dt}}={\mathrm{D}}_0+\mathrm{ϵ}{\mathrm{D}}_1+\mathrm{O}\left({\mathrm{ϵ}}^2\right)$ and $\frac{{\mathrm{d}}^2}{\mathrm{d}{\mathrm{t}}^2}={\mathrm{D}}_0^2+{\mathrm{ϵ}}^2{\mathrm{D}}_1^2+2\mathrm{ϵ}{\mathrm{D}}_0{\mathrm{D}}_1+\mathrm{O}\left({\mathrm{ϵ}}^2\right)$ where ${\mathrm{D}}_{\mathrm{i}}^{\mathrm{j}}=\frac{{\partial }^{\mathrm{j}}}{{\partial }^{\mathrm{j}}{\mathrm{T}}_{\mathrm{i}}}$. By substituting (6) and (7) into (4) and (5) and equating the coefficients of like powers of $\mathrm{ϵ}$, the following hierarchy of equations to second order is obtained:

$${\mathrm{D}}_0^2{\mathrm{x}}_0+\frac{{\mathrm{\omega }}^2}{4}{\mathrm{x}}_0=0$$E8

$${\mathrm{D}}_0{\mathrm{v}}_0+\mathrm{\beta }{\mathrm{v}}_0+\mathrm{\kappa }{\mathrm{D}}_0{\mathrm{x}}_0={\mathrm{\alpha }}_3{\mathrm{v}}_{0{\mathrm{\tau }}_2}$$E9

$${\mathrm{D}}_0^2{\mathrm{x}}_1+\frac{{\mathrm{\omega }}^2}{4}{\mathrm{x}}_1=-2{\mathrm{D}}_0{\mathrm{D}}_1{\mathrm{x}}_0-\stackrel{\sim }{\mathrm{\delta }}{\mathrm{D}}_0{\mathrm{x}}_0-\stackrel{\sim }{\mathrm{\lambda }}{\mathrm{x}}_0^2{\mathrm{D}}_0{\mathrm{x}}_0-\stackrel{\sim }{\mathrm{\sigma }}{\mathrm{x}}_0-\stackrel{\sim }{\mathrm{\gamma }}{\mathrm{x}}_0^3+\stackrel{\sim }{\mathrm{\chi }}{\mathrm{v}}_0+\left({\stackrel{\sim }{\mathrm{\alpha }}}_1+{\stackrel{\sim }{\mathrm{\alpha }}}_{2}cos\left(\mathrm{\omega t}\right)\right){\mathrm{x}}_{0{\mathrm{\tau }}_1}$$E10

$${\mathrm{D}}_0{\mathrm{v}}_1+\mathrm{\beta }{\mathrm{v}}_1=-{\mathrm{D}}_1{\mathrm{v}}_0-\mathrm{\kappa }{\mathrm{D}}_0{\mathrm{x}}_1-\mathrm{\kappa }{\mathrm{D}}_1{\mathrm{x}}_0+{\mathrm{\alpha }}_3{\mathrm{v}}_{1{\mathrm{\tau }}_2}$$E11

Initially, the solution can be characterized by

$${\mathrm{x}}_0\left({\mathrm{T}}_0,{\mathrm{T}}_1\right)=\mathrm{A}\left({\mathrm{T}}_1\right){\mathrm{e}}^{\mathrm{i}\frac{\mathrm{\omega }{\mathrm{T}}_0}{2}}+\overline{\mathrm{A}}\left({\mathrm{T}}_1\right){\mathrm{e}}^{-\mathrm{i}\frac{\mathrm{\omega }{\mathrm{T}}_0}{2}}$$E12

$${\mathrm{v}}_0\left({\mathrm{T}}_0,{\mathrm{T}}_1\right)=\frac{-\text{κiωA}\left({\mathrm{T}}_1\right)}{2\mathrm{\beta }+\mathrm{i\omega }-2{\mathrm{\alpha }}_3{\mathrm{e}}^{-\mathrm{i}\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}}}{\mathrm{e}}^{\mathrm{i}\frac{\mathrm{\omega }{\mathrm{T}}_0}{2}}+\frac{\mathrm{\kappa i\omega }\overline{\mathrm{A}}\left({\mathrm{T}}_1\right)}{2\mathrm{\beta }-\mathrm{i\omega }-2{\mathrm{\alpha }}_3{\mathrm{e}}^{\mathrm{i}\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}}}{\mathrm{e}}^{-\mathrm{i}\frac{\mathrm{\omega }{\mathrm{T}}_0}{2}}$$E13

where $\mathrm{A}\left({\mathrm{T}}_1\right)$ and $\overline{\mathrm{A}}\left({\mathrm{T}}_1\right)$ are unknown complex conjugates. Substituting eqs. (12) and (13) into (10) and (11) and deleting the secular terms results in:

$$\begin{array}{l}-\mathrm{i\omega }\left({\mathrm{D}}_1\mathrm{A}\right)-\frac{\mathrm{i}\stackrel{\sim }{\mathrm{\delta }}\mathrm{\omega A}}{2}-\frac{\mathrm{i}\stackrel{\sim }{\mathrm{\lambda }}\mathrm{\omega }}{2}{\mathrm{A}}^2\overline{\mathrm{A}}-\stackrel{\sim }{\mathrm{\sigma }}\mathrm{A}-3\stackrel{\sim }{\mathrm{\gamma }}{\mathrm{A}}^2\overline{\mathrm{A}}-\frac{\mathrm{\kappa i\omega }\stackrel{\sim }{\mathrm{\chi }}\mathrm{A}}{2\mathrm{\beta }-2{\mathrm{\alpha }}_{3}cos\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)+\mathrm{i}\left(\mathrm{\omega }+2{\mathrm{\alpha }}_{3}sin\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}\\ +{\stackrel{\sim }{\mathrm{\alpha }}}_1\mathrm{A}{\mathrm{e}}^{-\frac{\mathrm{i\omega }{\mathrm{\tau }}_1}{2}}+\frac{{\stackrel{\sim }{\mathrm{\alpha }}}_2}{2}\overline{\mathrm{A}}{\mathrm{e}}^{\frac{\mathrm{i\omega }{\mathrm{\tau }}_1}{2}}\\ =0\end{array}$$E14

We can obtain modulation equations up to first order by expressing $\mathrm{A}=\frac{1}{2}\mathrm{a}{\mathrm{e}}^{\mathrm{i\theta }}$, where $\mathrm{a}$ and $\mathrm{\theta }$ are the amplitude and phase, respectively.

$$\left\{\begin{array}{ll}\frac{\mathrm{da}}{\mathrm{dt}}={\mathrm{S}}_1\mathrm{a}+{\mathrm{S}}_2{\mathrm{a}}^3+{\mathrm{S}}_3\text{acos}\left(2\mathrm{\theta }\right)+{\mathrm{S}}_4\text{asin}\left(2\mathrm{\theta }\right)& \\ & \\ \mathrm{a}\frac{\mathrm{d\theta }}{\mathrm{dt}}={\mathrm{S}}_5\mathrm{a}+{\mathrm{S}}_6{\mathrm{a}}^3+{\mathrm{S}}_4\text{acos}\left(2\mathrm{\theta }\right)-{\mathrm{S}}_3\text{asin}\left(2\mathrm{\theta }\right)& \end{array}\right.$$E15

where ${\mathrm{S}}_{\mathrm{i}}\left(\mathrm{i}=1,\dots ,6\right)$ are given in Appendix. The solution to the first order given by (12) and (13) is as follows ${\mathrm{x}}_0\left({\mathrm{T}}_0,{\mathrm{T}}_1\right)=\text{acos}\left(\frac{\mathrm{\omega t}}{2}+\mathrm{\theta }\right)$ and ${\mathrm{v}}_0\left({\mathrm{T}}_0,{\mathrm{T}}_1\right)=\text{Vcos}\left(\frac{\mathrm{\omega t}}{2}+\mathrm{\theta }+arctan\frac{2\mathrm{\beta }-2{\mathrm{\alpha }}_{3}cos\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)}{\mathrm{\omega }+2{\mathrm{\alpha }}_{3}sin\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)}\right)$ such that the condition $\mathrm{\omega }+2{\mathrm{\alpha }}_{3}sin\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\ne 0$ is satisfied. Moreover the voltage amplitude $\mathrm{V}$ is given by

$$\mathrm{V}=\frac{\mathrm{\kappa \omega }}{\sqrt{{\left(2\mathrm{\beta }-2{\mathrm{\alpha }}_{3}cos\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}^2+{\left(\mathrm{\omega }+2{\mathrm{\alpha }}_{3}sin\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}^2}}\mathrm{a}$$E16

To find the steady-state response of system (15), which represents the periodic oscillations of Eqs. (4) and (5), we set $\frac{\mathrm{da}}{\mathrm{dt}}$ and $\frac{\mathrm{d\theta }}{\mathrm{dt}}$ equal to zero. After eliminating the phase, we are left with an algebraic equation in amplitude $\mathrm{a}$

$${\left({\mathrm{S}}_1\mathrm{a}+{\mathrm{S}}_2{\mathrm{a}}^3\right)}^2+{\left({\mathrm{S}}_5\mathrm{a}+{\mathrm{S}}_6{\mathrm{a}}^3\right)}^2=\left({\mathrm{S}}_3^2+{\mathrm{S}}_4^2\right){\mathrm{a}}^2$$E17

To find the average power, we integrate the dimensionless form of the instantaneous power ($\mathrm{P}\left(\mathrm{t}\right)=\mathrm{\beta v}{\left(\mathrm{t}\right)}^2$) over one period of the delay modulation $\mathrm{T}$. This results in the following expression:

$${\mathrm{P}}_{\mathrm{av}}=\frac{1}{\mathrm{T}}{\int }_0^{\mathrm{T}}‍\mathrm{\beta }{\mathrm{v}}^2\mathrm{dt}$$E18

where $\mathrm{T}=\frac{4\mathrm{\pi }}{\mathrm{\omega }}$. Then, the average power expressed by ${\mathrm{P}}_{\mathrm{av}}=\frac{\mathrm{\beta }{\mathrm{V}}^2}{2}$ reads

$${\mathrm{P}}_{\mathrm{av}}=\frac{1}{2}\left[\frac{\mathrm{\beta }{\mathrm{\kappa }}^2{\mathrm{\omega }}^2}{{\left(2\mathrm{\beta }-2{\mathrm{\alpha }}_{3}cos\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}^2+{\left(\mathrm{\omega }+2{\mathrm{\alpha }}_{3}sin\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}^2}\right]{\mathrm{a}}^2$$E19

The value of amplitude “$\mathrm{a}$” is derived from Eq. (17). We determine the highest power response by employing a maximizing process, which yields

$${\mathrm{P}}_{\max }=\frac{\mathrm{\beta }{\mathrm{\kappa }}^2{\mathrm{\omega }}^2{\mathrm{a}}^2}{{\left(2\mathrm{\beta }-2{\mathrm{\alpha }}_{3}cos\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}^2+{\left(\mathrm{\omega }+2{\mathrm{\alpha }}_{3}sin\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}^2}$$E20

The impact of various system parameters on both the steady-state response and maximum power output of the harvester is evaluated using Eqs. (17) and (20).

3. Vibration-based electromagnetic harvesting with quasi-periodic behavior

We investigate the response of the harvester system under quick perturbations by using the second-step perturbation method [17]. To make the calculation easier, we transform the polar form (15) into a Cartesian system using a variable change $\mathrm{u}=\text{acosθ}$ and $\mathrm{w}=-\text{asinθ}$

$$\left\{\begin{array}{ll}\frac{\mathrm{du}}{\mathrm{dt}}=\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)\mathrm{w}+\mathrm{\eta }\left\{\left({\mathrm{S}}_1+{\mathrm{S}}_3\right)\mathrm{u}+\left({\mathrm{S}}_2\mathrm{u}+{\mathrm{S}}_6\mathrm{w}\right)\left({\mathrm{u}}^2+{\mathrm{w}}^2\right)\right\}& \\ & \\ \frac{\mathrm{dw}}{\mathrm{dt}}=-\left({\mathrm{S}}_5+{\mathrm{S}}_4\right)\mathrm{u}+\mathrm{\eta }\left\{\left({\mathrm{S}}_1-{\mathrm{S}}_3\right)\mathrm{w}+\left({\mathrm{S}}_2\mathrm{w}-{\mathrm{S}}_6\mathrm{u}\right)\left({\mathrm{u}}^2+{\mathrm{w}}^2\right)\right\}& \end{array}\right.$$E21

where $\mathrm{\eta }$ is a record-keeping factor included to account for damping and nonlinearity. The periodic solution of the slow flow (21) can be represented in the form of

$$\mathrm{u}\left(\mathrm{t}\right)={\mathrm{u}}_0\left({\mathrm{T}}_1,{\mathrm{T}}_2\right)+\mathrm{\eta }{\mathrm{u}}_1\left({\mathrm{T}}_1,{\mathrm{T}}_2\right)+\mathrm{O}\left({\mathrm{\eta }}^2\right)$$E22

$$\mathrm{w}\left(\mathrm{t}\right)={\mathrm{w}}_0\left({\mathrm{T}}_1,{\mathrm{T}}_2\right)+\mathrm{\eta }{\mathrm{w}}_1\left({\mathrm{T}}_1,{\mathrm{T}}_2\right)+\mathrm{O}\left({\mathrm{\eta }}^2\right)$$E23

where ${\mathrm{T}}_1=\mathrm{t}$ and ${\mathrm{T}}_2=\mathrm{\eta t}$. In terms of the variables ${\mathrm{T}}_{\mathrm{i}}$, the time derivatives become $\frac{\mathrm{d}}{\mathrm{dt}}={\mathrm{D}}_1+\mathrm{\eta }{\mathrm{D}}_2+\mathrm{O}\left({\mathrm{\eta }}^2\right)$ where ${\mathrm{D}}_{\mathrm{i}}^{\mathrm{j}}=\frac{{\partial }^{\mathrm{j}}}{{\partial }^{\mathrm{j}}{\mathrm{T}}_{\mathrm{i}}}$. By substituting expressions (22) and (23) into (21), and comparing the coefficients of equivalent powers of $\mathrm{\eta }$, we arrive at the following set of equations

$${\mathrm{D}}_1^2{\mathrm{u}}_0+{\mathrm{\nu }}^2{\mathrm{u}}_0=0$$E24

$$\left({\mathrm{S}}_5-{\mathrm{S}}_4\right){\mathrm{w}}_0={\mathrm{D}}_1{\mathrm{u}}_0$$E25

$$\begin{array}{ll}{\mathrm{D}}_1^2{\mathrm{u}}_1+{\mathrm{\nu }}^2{\mathrm{u}}_1=& -{\mathrm{D}}_1{\mathrm{D}}_2{\mathrm{u}}_0+\left({\mathrm{S}}_1+{\mathrm{S}}_3\right){\mathrm{D}}_1{\mathrm{u}}_0+\left({\mathrm{S}}_2{\mathrm{D}}_1{\mathrm{u}}_0+{\mathrm{S}}_6{\mathrm{D}}_1{\mathrm{w}}_0\right)\left({\mathrm{u}}_0^2+{\mathrm{w}}_0^2\right)\\ & -\left({\mathrm{S}}_5-{\mathrm{S}}_4\right){\mathrm{D}}_2{\mathrm{w}}_0+\left({\mathrm{S}}_2{\mathrm{u}}_0+{\mathrm{S}}_6{\mathrm{w}}_0\right){\mathrm{D}}_1\left({\mathrm{u}}_0^2+{\mathrm{w}}_0^2\right)+\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)\left({\mathrm{S}}_1-{\mathrm{S}}_3\right){\mathrm{w}}_0\\ & +\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)\left({\mathrm{S}}_2{\mathrm{w}}_0-{\mathrm{S}}_6{\mathrm{u}}_0\right)\left({\mathrm{u}}_0^2+{\mathrm{w}}_0^2\right)\end{array}$$E26

$$\left({\mathrm{S}}_5-{\mathrm{S}}_4\right){\mathrm{w}}_1={\mathrm{D}}_1{\mathrm{u}}_1+{\mathrm{D}}_2{\mathrm{u}}_0-\left({\mathrm{S}}_1+{\mathrm{S}}_3\right){\mathrm{u}}_0-\left({\mathrm{S}}_2{\mathrm{u}}_0+{\mathrm{S}}_6{\mathrm{w}}_0\right)\left({\mathrm{u}}_0^2+{\mathrm{w}}_0^2\right)$$E27

where $\mathrm{\nu }=\sqrt{{\mathrm{S}}_5^2-{\mathrm{S}}_4^2}$ represents the frequency of the limit cycle of the slow flow, which is equal to the frequency of the quasi-periodic modulation. The initial solution can be represented as

$${\mathrm{u}}_0\left({\mathrm{T}}_1,{\mathrm{T}}_2\right)=\mathrm{R}\left({\mathrm{T}}_2\right)cos\left(\mathrm{\nu }{\mathrm{T}}_1+\mathrm{\psi }\left({\mathrm{T}}_2\right)\right)$$E28

$${\mathrm{w}}_0\left({\mathrm{T}}_1,{\mathrm{T}}_2\right)=\frac{-\mathrm{R}\left({\mathrm{T}}_2\right)\mathrm{\nu }}{\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)}sin\left(\mathrm{\nu }{\mathrm{T}}_1+\mathrm{\psi }\left({\mathrm{T}}_2\right)\right)$$E29

$\mathrm{R}$ and $\mathrm{\psi }$ represent the amplitude and phase, respectively, of the slow flow limit cycle. By substituting (28) and (29) into (26) and eliminating non-oscillatory terms, we arrive at the following autonomous slow-slow flow system for $\mathrm{R}$ and $\mathrm{\psi }$

$$\left\{\begin{array}{ll}\frac{\mathrm{dR}}{\mathrm{dt}}={\mathrm{S}}_1\mathrm{R}+\left[\frac{{\mathrm{S}}_2}{2}+\frac{{\mathrm{S}}_2{\mathrm{\nu }}^2}{2{\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)}^2}\right]{\mathrm{R}}^3& \\ & \\ \mathrm{R}\frac{\mathrm{d\psi }}{\mathrm{dt}}=\left[\frac{\mathrm{\nu }{\mathrm{S}}_6}{4\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)}+\frac{3{\mathrm{\nu }}^3{\mathrm{S}}_6}{8{\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)}^3}+\frac{3}{8\mathrm{\nu }}{\mathrm{S}}_6\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)\right]{\mathrm{R}}^3& \end{array}\right.$$E30

The equilibria of this slow-slow flow system determine the QP solutions of the original system (1), (2). To find the non-trivial equilibrium, we set $\frac{\mathrm{dR}}{\mathrm{dt}}$ equal to zero and obtain the following expression:

$$\mathrm{R}=\sqrt{-\frac{2{\mathrm{S}}_1{\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)}^2}{{\mathrm{S}}_2{\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)}^2+{\mathrm{S}}_2{\mathrm{\nu }}^2}}$$E31

Therefore, the approximate periodic solution for the slow flow (21) is given by

$$\mathrm{u}\left(\mathrm{t}\right)=\text{Rcos}\left(\mathrm{\theta t}\right)$$E32

$$\mathrm{w}\left(\mathrm{t}\right)=\frac{-\mathrm{R\nu }}{\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)}sin\left(\mathrm{\theta t}\right)$$E33

The approximate amplitude $\mathrm{a}\left(\mathrm{t}\right)$ of the QP response can be expressed as

$$\mathrm{a}\left(\mathrm{t}\right)=\sqrt{\frac{{\mathrm{R}}^2}{2}+\frac{{\mathrm{R}}^2{\mathrm{\nu }}^2}{2{\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)}^2}+\left[\frac{{\mathrm{R}}^2}{2}-\frac{{\mathrm{R}}^2{\mathrm{\nu }}^2}{2{\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)}^2}\right]cos\left(2\mathrm{\theta t}\right)}$$E34

and the envelope of the QP modulation is defined by ${\mathrm{a}}_{\min }$ and ${\mathrm{a}}_{\max }$ given by

$${\mathrm{a}}_{\min }=min\left\{\sqrt{\frac{{\mathrm{R}}^2}{2}+\frac{{\mathrm{R}}^2{\mathrm{\nu }}^2}{2{\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)}^2}\pm \left[\frac{{\mathrm{R}}^2}{2}-\frac{{\mathrm{R}}^2{\mathrm{\nu }}^2}{2{\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)}^2}\right]}\right\}$$E35

$${\mathrm{a}}_{\max }=max\left\{\sqrt{\frac{{\mathrm{R}}^2}{2}+\frac{{\mathrm{R}}^2{\mathrm{\nu }}^2}{2{\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)}^2}\pm \left[\frac{{\mathrm{R}}^2}{2}-\frac{{\mathrm{R}}^2{\mathrm{\nu }}^2}{2{\left({\mathrm{S}}_5-{\mathrm{S}}_4\right)}^2}\right]}\right\}$$E36

Therefore, the QP response of the original Eq. (2) can be explicitly expressed as

$$\mathrm{x}\left(\mathrm{t}\right)=\mathrm{u}\left(\mathrm{t}\right)cos\left(\frac{\mathrm{\omega t}}{2}\right)+\mathrm{w}\left(\mathrm{t}\right)sin\left(\frac{\mathrm{\omega t}}{2}\right)$$E37

The QP solution for the voltage v(t) can be found by inserting the expression for $\mathrm{x}\left(\mathrm{t}\right)$ from Eq. (37) into the second equation of system (2), and then evaluating the convolution integral with the boundary condition $\mathrm{v}\left(0\right)=\mathrm{v}\left(\mathrm{T}\right)$, where $\mathrm{T}=\frac{2\mathrm{\pi }}{\mathrm{\nu }}$. This results in

$$\mathrm{v}\left(\mathrm{t}\right)=-\mathrm{\kappa }{\mathrm{e}}^{\left({\mathrm{\alpha }}_3{\mathrm{e}}^{\mathrm{\beta }{\mathrm{\tau }}_2}-\mathrm{\beta }\right)\mathrm{t}}{\int }_0^{\mathrm{t}}‍\dot{\mathrm{x}}\left({\mathrm{t}}^{\prime }\right){\mathrm{e}}^{\left(\mathrm{\beta }-{\mathrm{\alpha }}_3{\mathrm{e}}^{\mathrm{\beta }{\mathrm{\tau }}_2}\right){\mathrm{t}}^{\prime }}{\mathrm{dt}}^{\prime }$$E38

As a result, the QP modulation region yields the following values for power: average, maximum, and output power

$${\mathrm{P}}_{\mathrm{QP}}\left(\mathrm{t}\right)=\mathrm{\beta }{\left(\mathrm{\kappa }{\mathrm{e}}^{\left({\mathrm{\alpha }}_3{\mathrm{e}}^{\mathrm{\beta }{\mathrm{\tau }}_2}-\mathrm{\beta }\right)\mathrm{t}}{\int }_0^{\mathrm{t}}‍\dot{\mathrm{x}}\left({\mathrm{t}}^{\prime }\right){\mathrm{e}}^{\left(\mathrm{\beta }-{\mathrm{\alpha }}_3{\mathrm{e}}^{\mathrm{\beta }{\mathrm{\tau }}_2}\right){\mathrm{t}}^{\prime }}{\mathrm{dt}}^{\prime }\right)}^2$$E39

$${\mathrm{P}}_{\text{avQP}}=\frac{\mathrm{\beta }{\mathrm{\kappa }}^2{\mathrm{\nu }}^2}{2\left[{\left(\mathrm{\beta }-2{\mathrm{\alpha }}_{3}cos\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}^2+{\left(\mathrm{\nu }+2{\mathrm{\alpha }}_{3}sin\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}^2\right]}{\mathrm{a}}^2$$E40

$${\mathrm{P}}_{\text{maxQP}}=\frac{\mathrm{\beta }{\mathrm{\kappa }}^2{\mathrm{\nu }}^2}{\left[{\left(\mathrm{\beta }-2{\mathrm{\alpha }}_{3}cos\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}^2+{\left(\mathrm{\nu }+2{\mathrm{\alpha }}_{3}sin\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}^2\right]}{\mathrm{a}}^2$$E41

The amplitude $\mathrm{a}$ in Eqs. (40) and (41) can be calculated using Eqs. (35) and (36).

In this work, the response of the energy harvester system near the delay parametric resonance is studied. The effect of different system parameters, including the unmodulated delay amplitude, on the steady-state response and the maximum output power of the harvester is examined. To obtain the quasi-periodic response and corresponding power, a second-step perturbation method is used. The variation of the amplitude of the periodic and quasi-periodic responses and the maximum output power is shown in Figure 2. The results of the analytical prediction are compared to numerical simulations obtained by using dde23 algorithm. The box inset in the figures presents time histories of the amplitude and power responses.

Figure 2.

Vibration and power amplitudes as a function of ${\mathrm{\alpha }}_1$ for ${\mathrm{\alpha }}_2=0.25$, ${\mathrm{\tau }}_1=5.2$, ${\mathrm{\tau }}_2=4.2$, $\mathrm{\chi }=0.05$, $\mathrm{\beta }=0.05$, $\mathrm{\lambda }=0.2$, $\mathrm{\delta }=-0.1$, $\mathrm{\gamma }=0.05$, $\mathrm{\omega }=2,$ and $\mathrm{\kappa }=0.5$. The black (gray) line represents the delayed (undelayed) electric circuit with ${\mathrm{\alpha }}_3$ equal to β (${\mathrm{\alpha }}_3$ equal to 0). The graph shows both the analytical predictions (solid lines for stable and dashed line for unstable) and numerical simulation results (represented by circles).

Figure 2 displays the change in amplitude of both the periodic and QP responses and the maximum output power (${\mathrm{P}}_{\max },{\mathrm{P}}_{\text{maxQP}}$) as the unmodulated delay amplitude ${\mathrm{\alpha }}_1$ increases. It compares the analytical predictions (solid lines for stable and dashed lines for unstable) with numerical simulations (represented by circles) from the dde23 algorithm [18]. The plots show that only periodic vibration-based energy can be extracted at a small ${\mathrm{\alpha }}_1$ value, but as it increases, the periodic solution becomes unstable and energy can only be extracted from QP vibration with improved performance compared to the periodic output power. The figure also demonstrates that the introduction of delay in the electrical circuit (${\mathrm{\alpha }}_3$≠ 0) decreases the amplitude of both the periodic and QP modulation (Figure 2a, black line) while increasing the harvested power (Figure 2b, black line), indicating that the maximum power output does not always correspond with the maximum oscillation amplitude.

4. Conclusions

This research evaluated the effectiveness of energy harvesting (EH) in a system comprising a delayed Duffing-van der Pol oscillator and a piezoelectric harvester with a delayed response. The delay in the mechanical and piezoelectric elements was assumed to have differences in both timing and magnitude. The examination centered around the vicinity of delay parametric resonance where the frequency of modulation was close to twice the natural frequency of the oscillator. Approximation methods were used to calculate the periodic and quasi-periodic vibrations for EH purposes. The influence of the delay in the piezoelectric subsystem on the EH performance of the delayed Duffing-van der Pol harvester was analyzed. The findings showed that the presence of a modulated delay in the mechanical subsystem leads to an ideal set of system parameters that maximizes the quasi-periodic vibration amplitude and corresponding output power. Moreover, the results indicated that the presence of a delay in the electrical circuit improves the output power.

Appendices

${\mathrm{S}}_1=\frac{-\mathrm{\delta }}{2}-\frac{\mathrm{\chi \kappa }\left(2\mathrm{\beta }-2{\mathrm{\alpha }}_{3}cos\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}{{\left(2\mathrm{\beta }-2{\mathrm{\alpha }}_{3}cos\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}^2+{\left(\mathrm{\omega }+2{\mathrm{\alpha }}_{3}sin\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}^2}-\frac{{\mathrm{\alpha }}_1}{\mathrm{\omega }}sin\frac{\mathrm{\omega }{\mathrm{\tau }}_1}{2}$, ${\mathrm{S}}_2=\frac{-\mathrm{\lambda }}{8}$, ${\mathrm{S}}_3=\frac{{\mathrm{\alpha }}_2}{2\mathrm{\omega }}sin\frac{\mathrm{\omega }{\mathrm{\tau }}_1}{2.}$

${\mathrm{S}}_4=-\frac{{\mathrm{\alpha }}_2}{2\mathrm{\omega }}cos\frac{\mathrm{\omega }{\mathrm{\tau }}_1}{2}$, ${\mathrm{S}}_5=\frac{\mathrm{\sigma }}{\mathrm{\omega }}+\frac{\mathrm{\chi \kappa }\left(\mathrm{\omega }+2{\mathrm{\alpha }}_{3}sin\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}{{\left(2\mathrm{\beta }-2{\mathrm{\alpha }}_{3}cos\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}^2+{\left(\mathrm{\omega }+2{\mathrm{\alpha }}_{3}sin\left(\frac{\mathrm{\omega }{\mathrm{\tau }}_2}{2}\right)\right)}^2}-\frac{{\mathrm{\alpha }}_1}{\mathrm{\omega }}cos\frac{\mathrm{\omega }{\mathrm{\tau }}_1}{2}$, ${\mathrm{S}}_6=\frac{3\mathrm{\gamma }}{4\mathrm{\omega }}$