A piezoelectric energy harvester for human body motion subjected to two different transversal reciprocating excitations

Harvesting energy from human body motion to supply electricity for wearable devices is focused on in this paper. Based on the fact that the frequency of human body motion is lower and the motions of different human body parts are variable, a piezoelectric energy harvester subjected to two different transversal reciprocating excitations is studied in this paper. Each excitation is treated as a transverse rheonomic constraint. The dynamics equation of the beam is established using the Hamiltonian principle. Expressing the transverse rheonomic constraint as a periodic function, closed-form solutions of the dynamics equation are obtained. And the characteristics of energy harvesters are investigated based on the closed-form solutions. The results show that the difference between the two excitations will certainly cause the energy harvester to generate more output power at lower frequencies of excitations, and the larger the difference, the more the output power will be generated. This unusual characteristic at the lower frequency enables the proposed harvester to be quite suitable to harvest energy from the motions of the human body.

1 Introduction

The market for wearable devices, such as pacemakers, hearing aids, and so on, has been rapidly growing for 3 decades. However, the energy density of batteries for these devices limits their usability, which forces the user to carry more battery packs or recharge their devices on a power grid somewhere. Among the various methods, harvesting energy from human movements may be a convenient solution to the problem.

Compared with other ambient mechanical vibrations, human movements have two features. One is that the frequency of human movement is much lower. It has been acknowledged that the frequency of human motion is less than 10 Hz (Green et al., 2013) and that of human upper limbs is 0.5–3 Hz (George et al., 2020). The other is that at different points in the human body, the motions are variable in amplitude, frequency, or even in regularity (Hallal and Pratt, 2020; Ren et al., 2020).

One efficient way to convert human movements into electronic energy is by using piezoelectric materials for the reason that piezoelectric ceramics have large power densities, ease of application, and feasibility of fabrication (Roundy and Wright, 2004; Sodano et al., 2004; Maamer et al., 2019; Blad and Tolou, 2019).

A normal piezoelectric energy harvester is a cantilever metallic beam attached to piezoelectric ceramics layers. The vibrating host structure that the metallic beam is fixed on is named an excitation or a vibrator (Sodano et al., 2020; Maamer et al., 2019). The operating principle of the normal piezoelectric energy harvester is based on the resonant interaction between the excitation and the cantilever beam. Therefore, only when the model frequency of the beam is equal, or very close, to the frequency of the excitation can the normal piezoelectric energy harvester generate maximum energy transduction.

Thus, for the vibration of a frequency lower than 30 Hz, the performance of the normal piezoelectric energy harvester is not as good as expected. Therefore, a large variety of configurations to modify the normal piezoelectric energy harvester have been proposed (Blad and Tolou, 2019). Liu et al. (2012) designed an S-shaped piezoelectric PZT (PbZr_xTi_(1−x)O₃) cantilever for harvesting vibration energy at lower than 30 Hz. Dechant et al. (2017) utilized mechanical impacts on two mechanical stoppers to convert low-frequency mechanical vibrations, the best at 7–25 Hz, into high-frequency resonance oscillations of the transducer. Leadenham and Erturk (2015) developed an M-shaped piezoelectric energy harvester configuration that exhibits a nonlinear frequency response under very low vibration levels. The secondary resonance of their harvester is as low as 4.5 Hz. Halim and Park (2014) used a low-frequency flexible beam with an extended tip mass driven by the vibrator to hit two high-frequency rigid piezoelectric generating beams. The best-performing bandwidth of this harvester is 7–14.5 Hz.

In this paper, we propose a piezoelectric energy harvester for human body motion. In contrast with the harvester whose configuration is changed from a normal piezoelectric energy harvester with one end of a metallic beam fixed on the excitation, the metallic beam of the harvester proposed is connected to the human body at two points, acting as two vibrators (see Fig. 1). Its proposition utilizes the feature that the motions of the two human body points are distinct. Its configuration is simpler than existing configurations of piezoelectric energy harvesters for lower-frequency vibration; consequently, the cost of fabrication may be reduced greatly but at the expense of losing a little comfort. The excitations are treated as transverse rheonomic constraints acting on the beam for the derivation of the dynamics equation of the energy harvester. Based on the analysis of dynamics equations, we find such a piezoelectric energy harvester performs very well, especially at frequencies lower than 30 Hz. We are inspired by the fact that a metallic beam of an energy harvester subjected to two different transversals reciprocating is quite suitable for harvesting energy from human body motions.

Figure 1 Schematic diagram of a piezoelectric energy harvester for human body motion subjected to two excitations.

The rest of the paper is organized as follows. In Sect. 2, the dynamics equation of the harvester with two excitations is derived, and its closed-form solution is obtained in the case when the excitations can be expressed as periodic functions. In Sect. 3, the performance of the energy harvester is analyzed when the two vibrators are with the same frequency but different amplitudes and/or phases, to illustrate the abilities of the harvester, especially at low frequencies. And in Sect. 4, the conclusions are presented.

2 Modeling of the harvester for human body motion subjected to two excitations

Figure 1 is a schematic diagram of the piezoelectric energy harvester for human body motion. The beam is connected at two points of the human body, which are referred as vibrator 1 and vibrator 2, respectively. As shown in Fig. 1, one end of the beam is clamped on vibrator 1, i.e., the back of the human body, and the other end of the beam is free, while vibrator 2, i.e., the arm, excites the beam at the middle part of the beam. It can be observed that the motion of vibrator 1 will be dissimilar than the motion of vibrator 2 when a human being is walking, jogging, or conducting other movements. Under this circumstance, we develop the piezoelectric energy harvester. Of course, the two vibrators may be on other parts of the human body, provided that the motions of the two vibrators are different. Otherwise, the beam will vibrate under excitation of its inertia force, similar to the normal piezoelectric energy harvester, which will also be discussed for comparison.

2.1 Derivation of the electromechanical model

The physical model and the simplified model for analysis of the harvester are presented in Figs. 2 and 3, respectively.

Figure 2 Physical model of the piezoelectric energy harvester for human body motion subjected to two excitations.

Figure 3 Simplified model for analysis.

In Fig. 2, the beam is clamped onto vibrator 1. For convenience, the effects of the clamping on the beam in Fig. 2 are represented by two elements in Fig. 3; one is a sliding constraint and the other one is a transverse rheonomic constraint, P₁(t) at a₁l. The former implies that the bending deformation of the beam at the clamping end is neglected, the end is tolerated to move in the direction of w, and the latter allows us to consider a more general case to establish the dynamics equation. In this paper, a₁ will be set to zero. Vibrator 2 is also a transverse rheonomic constraint, P₂(t) acting at a₂l ($\mathrm{0}\le a_{\mathrm{2}}\le \mathrm{1}$).

We assume the beam is an Euler–Bernoulli beam; both transverse displacements, P₁(t) and P₂(t), can be represented by modal coordinates. We ignore the inherent piezoelectric nonlinearities of the harvester for the reason that the acceleration is very small (Erturk and Inman, 2011; Wagner and Hagedorn, 2002).

Let the transverse displacement of the beam be w(x,t). Then, vibrators 1 and 2 can be expressed as a constraint equation as follows:

$$\begin{array}{}\text{(1)}& {\stackrel{\mathrm{̃}}{p}}_h\left(t\right)=w\left(a_{h}l,t\right)-p_h\left(t\right)=\mathrm{0}(h=\mathrm{1},\mathrm{2}).\end{array}$$

Let the rth order modal function be ϕ_r(x) and the modal coordinate be q_r(t); the transverse displacement of the beam is

$$\begin{array}{}\text{(2)}& w\left(x,t\right)=\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\varphi }}_r\left(x\right)q_r\left(t\right).\end{array}$$

With the constraints depicted in Fig. 3, the boundary conditions of the beam is ${\mathit{\varphi }}^{\prime }\left(\mathrm{0}\right)=\mathrm{0},{\mathit{\varphi }}^{\prime \prime \prime }\left(\mathrm{0}\right)=\mathrm{0},{\mathit{\varphi }}^{\prime \prime }\left(l\right)=\mathrm{0},and{\mathit{\varphi }}^{\prime \prime \prime }\left(l\right)=\mathrm{0}$. Thus, the rth normalized model function is

$$\begin{array}{}\text{(3)}& {\mathit{\varphi }}_r\left(x\right)=\left\{\begin{array}{l}\sqrt{\frac{\mathrm{2}}{ml}}\frac{\cos \left(\frac{{\mathit{\lambda }}_r}{l}x\right)+\frac{\cos \left({\mathit{\lambda }}_r\right)}{\cosh \left({\mathit{\lambda }}_r\right)}\cosh \left(\frac{{\mathit{\lambda }}_r}{l}x\right)}{\sqrt{\mathrm{1}+\text{csch}{\left({\mathit{\lambda }}_r\right)}^{\mathrm{2}}\sin {\left({\mathit{\lambda }}_r\right)}^{\mathrm{2}}}}(r\ge \mathrm{2})\\ \sqrt{{\displaystyle \frac{\mathrm{1}}{ml}}}(r=\mathrm{1})\end{array}\right.,\end{array}$$

where λ_r can be obtained by solving the following equation:

$$\begin{array}{}\text{(4)}& \cosh \left(\mathit{\lambda }\right)\sin \left(\mathit{\lambda }\right)+\cos \left(\mathit{\lambda }\right)\sinh \left(\mathit{\lambda }\right)=\mathrm{0}.\end{array}$$

The kinetic energy of the beam can be written as

$$\begin{array}{}\text{(5)}& T={\displaystyle \frac{m}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{l}{\int }}{\dot{w}}^{\mathrm{2}}\mathrm{d}x,\end{array}$$

where $m=b({\mathit{\rho }}_s{h}_s+\mathrm{2}{\mathit{\rho }}_{\mathrm{p}}{h}_{\mathrm{p}})$.

The potential energy of the beam can be written as

$$\begin{array}{}\text{(6)}& U={\displaystyle \frac{\text{YI}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{l}{\int }}{w^{\prime \prime }}^{\mathrm{2}}\mathrm{d}x,\end{array}$$

where YI $=c_{\mathrm{11}}{I}_p+E{I}_s$.

The complementary energy of the beam can be written as

$$\begin{array}{}\text{(7)}& W_e=-U+\underset{\mathrm{0}}{\overset{l}{\int }}{\displaystyle \frac{{\mathit{\alpha }}^{\mathrm{2}}{\mathit{\epsilon }}_{\mathrm{33}}b{v}^{\mathrm{2}}}{\mathrm{4}{h}_{\mathrm{p}}}}-\mathit{\alpha }{e}_{\mathrm{31}}{h}_{\mathrm{pc}}bv{w}^{\prime \prime }\mathrm{d}x,\end{array}$$

where $h_{\mathrm{pc}}=(h_s+h_{\mathrm{p}})/\mathrm{2}$ and α=2 when the piezoelectric layers are connected in parallel, and α=1 when the piezoelectric layers are connected in series.

The virtual work of air damping and material damping can be written as

$$\begin{array}{}\text{(8)}& Q_s=-\underset{\mathrm{0}}{\overset{l}{\int }}{c}_s{I}_{\mathrm{b}}{\dot{w}}^{\prime \prime }\mathit{\delta }{w}^{\prime \prime }-c_{\mathrm{a}}\dot{w}\mathit{\delta }w\mathrm{d}x.\end{array}$$

As a conservative element, the virtual work of the load resistance is

$$\begin{array}{}\text{(9)}& Q_R=-{\displaystyle \frac{v}{R}}\mathit{\delta }v.\end{array}$$

Using κ₁(t) and κ₂(t) to represent Lagrange multipliers, by the Hamiltonian principle, the dynamic model can be transformed into a functional extremum problem as follows:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& \underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\delta }\left[T+W_e+\sum _{i=\mathrm{1}}^{\mathrm{2}}{\mathit{\kappa }}_i\left(t\right){\stackrel{\mathrm{̃}}{p}}_i\left(t\right)\right]\mathrm{d}t+\underset{\mathrm{0}}{\overset{t}{\int }}{Q}_s\mathrm{d}t\\ & +\underset{\mathrm{0}}{\overset{t}{\int }}{Q}_R\mathrm{d}t=\mathrm{0}.\end{array}\end{array}$$

With the procedure of derivation presented in Appendix B, the governing equations of the energy harvester can be written as

$$\begin{array}{}\text{(11)}& \left\{\begin{array}{l}{\ddot{q}}_r\left(t\right)+\mathrm{2}{\mathit{\zeta }}_r{\mathit{\omega }}_r{\dot{q}}_r\left(t\right)+{\mathit{\omega }}_r^{\mathrm{2}}{q}_r\left(t\right)-{\mathit{\chi }}_{r}v\left(t\right)\\ =\sum _{i=\mathrm{1}}^{\mathrm{2}}{\mathit{\kappa }}_i\left(t\right)\frac{{\mathit{\varphi }}_r\left(a_{i}l\right)}{m}\\ \sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right)q_r\left(t\right)-p_{\mathrm{1}}\left(t\right)=\mathrm{0}\\ \sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right)q_r\left(t\right)-p_{\mathrm{2}}\left(t\right)=\mathrm{0}\\ C_{p}v\left(t\right)+\frac{v\left(t\right)}{R}+\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\chi }}_r{\dot{q}}_r\left(t\right)=\mathrm{0}\left\{r=\mathrm{1},\mathrm{\dots },\mathrm{\infty }\right\}\end{array}\right..\end{array}$$

It can be seen that Eq. (11) is an infinite-order differential algebraic equation (DAE). From Eq. (11), we obtain the undamped natural frequency of the beam ω_r as follows:

$$\begin{array}{}\text{(12)}& {\mathit{\omega }}_r={\mathit{\lambda }}_r^{\mathrm{2}}\sqrt{{\displaystyle \frac{\text{YI}}{m{l}^{\mathrm{4}}}}},\end{array}$$

and the expressions of χ_r and C_p as follows:

$$\begin{array}{}\text{(13)}& {\mathit{\chi }}_r=\mathit{\alpha }{e}_{\mathrm{31}}{h}_{\mathrm{p}}b{\mathit{\varphi }}_r^{\prime }(x{)}_{x=l},\end{array}$$

$$\begin{array}{}\text{(14)}& C_p={\displaystyle \frac{\mathit{\alpha }{\mathit{\epsilon }}_{\mathrm{33}}bl}{\mathrm{2}{h}_{\mathrm{p}}}}.\end{array}$$

The damping ratio ζ_r for λ_r can be determined, according to the formulas of the damping ratio in Erturk and Inman (2011), by

$$\begin{array}{}\text{(15)}& {\mathit{\zeta }}_r={\displaystyle \frac{{\mathit{\zeta }}_j{\mathit{\lambda }}_j^{\mathrm{2}}\left({\mathit{\lambda }}_i^{\mathrm{4}}-{\mathit{\lambda }}_r^{\mathrm{4}}\right)+{\mathit{\zeta }}_i{\mathit{\lambda }}_i^{\mathrm{2}}\left({\mathit{\lambda }}_r^{\mathrm{4}}-{\mathit{\lambda }}_j^{\mathrm{4}}\right)}{\left({\mathit{\lambda }}_i^{\mathrm{4}}-{\mathit{\lambda }}_j^{\mathrm{4}}\right){\mathit{\lambda }}_r^{\mathrm{2}}}}.\end{array}$$

It should be noted that the subscripts i, j, and r can be chosen independently in Eq. (15), provided that neither i nor j is equal to zero.

2.2 Closed-form solution of modeling under a periodical displacement excitation

When a human being is walking or jogging, the motions of their back and upper arm are periodic, and the experiment results of Stansfield et al. (2001) indicate that the frequencies of both parts of the human body are the same. In this paper, for simplification, we suppose that both transverse rheonomic constraints, i.e., both displacement functions of excitation P_h(t) and $h=\mathrm{1},\mathrm{2}$ are periodic functions with the same frequency but different amplitude, so they can be expressed as

$$\begin{array}{}\text{(16a)}& {\displaystyle }{P}_{\mathrm{1}}\left(t\right)=f_{\mathrm{1}}{e}^{I\mathit{\varpi }t}\text{(16b)}& {\displaystyle }{P}_{\mathrm{2}}\left(t\right)=f_{\mathrm{2}}{e}^{I\mathit{\varpi }t}.\end{array}$$

Then the form of the rth modal coordinates can be written as

$$\begin{array}{}\text{(17)}& q_r\left(t\right)=H_r{e}^{I\mathit{\varpi }t}.\end{array}$$

Similarly, the ith Lagrange multiplier can be written as

$$\begin{array}{}\text{(18)}& {\mathit{\kappa }}_i\left(t\right)={\mathit{\kappa }}_i{e}^{I\mathit{\varpi }t},\end{array}$$

and the output voltage of the energy harvester is in the form of

$$\begin{array}{}\text{(19)}& v\left(t\right)=v{e}^{I\mathit{\varpi }t}.\end{array}$$

Substituting Eqs. (16) to (19) into Eq. (11), we get

$$\begin{array}{}\text{(20)}& \left\{\begin{array}{l}-{\mathit{\varpi }}^{\mathrm{2}}{H}_r{e}^{I\mathit{\varpi }t}+\mathrm{2}I\mathit{\varpi }{\mathit{\zeta }}_r{\mathit{\omega }}_r{H}_r{e}^{I\mathit{\varpi }t}\\ +{\mathit{\omega }}_r^{\mathrm{2}}{H}_r{e}^{I\mathit{\varpi }t}-{\mathit{\chi }}_{r}v{e}^{I\mathit{\varpi }t}\\ ={\mathit{\kappa }}_{\mathrm{1}}\frac{{\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right)}{m}{e}^{I\mathit{\varpi }t}+{\mathit{\kappa }}_{\mathrm{2}}\frac{{\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right)}{m}{e}^{I\mathit{\varpi }t}\\ \sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right)H_r{e}^{I\mathit{\varpi }t}-f_{\mathrm{1}}{e}^{I\mathit{\varpi }t}=\mathrm{0}\\ \sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right)H_r{e}^{I\mathit{\varpi }t}-f_{\mathrm{2}}{e}^{I\mathit{\varpi }t}=\mathrm{0}\\ \left(I\mathit{\varpi }{C}_p+\frac{\mathrm{1}}{R}\right)v{e}^{I\mathit{\varpi }t}+I\mathit{\varpi }\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\chi }}_r{H}_r{e}^{I\mathit{\varpi }t}=\mathrm{0}\\ \left\{r=\mathrm{1},\mathrm{\dots },\mathrm{\infty }\right\}\end{array}\right..\end{array}$$

The rth order modal equation in Eq. (20) yields

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& \left({\mathit{\omega }}_r^{\mathrm{2}}-{\mathit{\varpi }}^{\mathrm{2}}+I\mathit{\varpi }\mathrm{2}{\mathit{\zeta }}_r{\mathit{\omega }}_r\right)H_r\\ & =\left({\mathit{\kappa }}_{\mathrm{1}}{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right)}{m}}+{\mathit{\kappa }}_{\mathrm{2}}{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right)}{m}}+{\mathit{\chi }}_{r}v\right).\end{array}\end{array}$$

Let ${\mathrm{\Gamma }}_r=({\mathit{\omega }}_r^{\mathrm{2}}-{\mathit{\varpi }}^{\mathrm{2}}+I\mathit{\varpi }\mathrm{2}{\mathit{\zeta }}_r{\mathit{\omega }}_r{)}^{-\mathrm{1}}$; from Eq. (21), the modal coordinate can be written as

$$\begin{array}{}\text{(22)}& H_r={\mathrm{\Gamma }}_r({\mathit{\kappa }}_{\mathrm{1}}{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right)}{m}}+{\mathit{\kappa }}_{\mathrm{2}}{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right)}{m}}+{\mathit{\chi }}_{r}v).\end{array}$$

And from Eq. (17), we have

$$\begin{array}{}\text{(23)}& q_r\left(t\right)={\mathrm{\Gamma }}_r({\mathit{\kappa }}_{\mathrm{1}}{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right)}{m}}+{\mathit{\kappa }}_{\mathrm{2}}{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right)}{m}}+{\mathit{\chi }}_{r}v)e^{I\mathit{\varpi }t}.\end{array}$$

Substituting Eq. (23) into the algebraic equations in Eq. (20), we get

$$\begin{array}{}\text{(24)}& \begin{array}{rl}\sum _{r=\mathrm{1}}^{\mathrm{\infty }}& {\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right){\mathrm{\Gamma }}_r({\mathit{\kappa }}_{\mathrm{1}}{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right)}{m}}+{\mathit{\kappa }}_{\mathrm{2}}{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right)}{m}}+{\mathit{\chi }}_{r}v)e^{I\mathit{\varpi }t}\\ & -f_{\mathrm{1}}{e}^{I\mathit{\varpi }t}=\mathrm{0}\\ \sum _{r=\mathrm{1}}^{\mathrm{\infty }}& {\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right){\mathrm{\Gamma }}_r({\mathit{\kappa }}_{\mathrm{1}}{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right)}{m}}+{\mathit{\kappa }}_{\mathrm{2}}{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right)}{m}}+{\mathit{\chi }}_{r}v)e^{I\mathit{\varpi }t}\\ & -f_{\mathrm{2}}{e}^{I\mathit{\varpi }t}=\mathrm{0}.\end{array}\end{array}$$

From Eq. (24), it can be seen that the coefficients of e^Iϖt in Eq. (24) should satisfy the following equations:

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{\mathit{\kappa }}_{\mathrm{1}}\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathrm{\Gamma }}_r{\displaystyle \frac{{\mathit{\varphi }}_r(a_{\mathrm{1}}l{)}^{\mathrm{2}}}{m}}& +{\mathit{\kappa }}_{\mathrm{2}}\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathrm{\Gamma }}_r{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right){\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right)}{m}}\\ & +\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right){\mathrm{\Gamma }}_r{\mathit{\chi }}_{r}v=f_{\mathrm{1}}\\ {\mathit{\kappa }}_{\mathrm{1}}\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathrm{\Gamma }}_r{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right){\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right)}{m}}& +{\mathit{\kappa }}_{\mathrm{2}}\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathrm{\Gamma }}_r{\displaystyle \frac{{\mathit{\varphi }}_r(a_{\mathrm{2}}l{)}^{\mathrm{2}}}{m}}\\ & +\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right){\mathrm{\Gamma }}_r{\mathit{\chi }}_{r}v=f_{\mathrm{2}}.\end{array}\end{array}$$

Substituting Eq. (23) into the electronic circuit equation in Eq. (20), we have

$$\begin{array}{}\text{(26)}& \begin{array}{rl}& \left(I\mathit{\varpi }{C}_p+{\displaystyle \frac{\mathrm{1}}{R}}+I\mathit{\varpi }\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\chi }}_r^{\mathrm{2}}{\mathrm{\Gamma }}_r\right)v\\ & +I\mathit{\varpi }{\mathit{\kappa }}_{\mathrm{1}}\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\chi }}_r{\mathrm{\Gamma }}_r{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right)}{m}}\\ & +I\mathit{\varpi }{\mathit{\kappa }}_{\mathrm{2}}\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\chi }}_r{\mathrm{\Gamma }}_r{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right)}{m}}=\mathrm{0}.\end{array}\end{array}$$

Let $d=(I\mathit{\varpi }{C}_p+\frac{\mathrm{1}}{R}+I\mathit{\varpi }\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\chi }}_r^{\mathrm{2}}{\mathrm{\Gamma }}_r{)}^{-\mathrm{1}}$, then the formula to calculate the output voltage is

$$\begin{array}{}\text{(27)}& v=-dI\mathit{\varpi }\left({\mathit{\kappa }}_{\mathrm{1}}\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\chi }}_r{\mathrm{\Gamma }}_r{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right)}{m}}+{\mathit{\kappa }}_{\mathrm{2}}\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\chi }}_r{\mathrm{\Gamma }}_r{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right)}{m}}\right).\end{array}$$

Substituting Eq. (27) into Eq. (24), we obtain

$$\begin{array}{}\text{(28)}& \begin{array}{l}{\mathit{\kappa }}_{\mathrm{1}}{\mathit{\xi }}_{\mathrm{11}}+{\mathit{\kappa }}_{\mathrm{2}}{\mathit{\xi }}_{\mathrm{12}}=f_{\mathrm{1}}\\ {\mathit{\kappa }}_{\mathrm{1}}{\mathit{\xi }}_{\mathrm{21}}+{\mathit{\kappa }}_{\mathrm{2}}{\mathit{\xi }}_{\mathrm{22}}=f_{\mathrm{2}},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(29)}& {\mathit{\xi }}_{ij}=\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathrm{\Gamma }}_r{\mathit{\varphi }}_r\left(a_{i}l\right)\left[{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{j}l\right)}{m}}-{\mathit{\chi }}_{r}dI\mathit{\varpi }\sum _{i=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\chi }}_i{\mathrm{\Gamma }}_i{\displaystyle \frac{{\mathit{\varphi }}_i\left(a_{j}l\right)}{m}}\right],\end{array}$$

we have

$$\begin{array}{}\text{(30a)}& {\displaystyle }{\mathit{\kappa }}_{\mathrm{1}}={\displaystyle \frac{f_{\mathrm{2}}{\mathit{\xi }}_{\mathrm{12}}-f_{\mathrm{1}}{\mathit{\xi }}_{\mathrm{22}}}{{\mathit{\xi }}_{\mathrm{12}}{\mathit{\xi }}_{\mathrm{21}}-{\mathit{\xi }}_{\mathrm{11}}{\mathit{\xi }}_{\mathrm{22}}}}\text{(30b)}& {\displaystyle }{\mathit{\kappa }}_{\mathrm{2}}={\displaystyle \frac{f_{\mathrm{1}}{\mathit{\xi }}_{\mathrm{21}}-f_{\mathrm{2}}{\mathit{\xi }}_{\mathrm{11}}}{{\mathit{\xi }}_{\mathrm{12}}{\mathit{\xi }}_{\mathrm{21}}-{\mathit{\xi }}_{\mathrm{11}}{\mathit{\xi }}_{\mathrm{22}}}}.\end{array}$$

Substituting Eq. (30) into Eq. (27), we get the formula to calculate the amplitude of the output voltage as follows:

$$\begin{array}{}\text{(31)}& \begin{array}{rl}v& =-dI\mathit{\varpi }{\displaystyle \frac{f_{\mathrm{2}}{\mathit{\xi }}_{\mathrm{12}}-f_{\mathrm{1}}{\mathit{\xi }}_{\mathrm{22}}}{{\mathit{\xi }}_{\mathrm{12}}{\mathit{\xi }}_{\mathrm{21}}-{\mathit{\xi }}_{\mathrm{12}}{\mathit{\xi }}_{\mathrm{22}}}}\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\chi }}_r{\mathrm{\Gamma }}_r{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{1}}l\right)}{m}}\\ & -dI\mathit{\varpi }{\displaystyle \frac{f_{\mathrm{1}}{\mathit{\xi }}_{\mathrm{21}}-f_{\mathrm{2}}{\mathit{\xi }}_{\mathrm{11}}}{{\mathit{\xi }}_{\mathrm{12}}{\mathit{\xi }}_{\mathrm{21}}-{\mathit{\xi }}_{\mathrm{11}}{\mathit{\xi }}_{\mathrm{22}}}}\sum _{r=\mathrm{1}}^{\mathrm{\infty }}{\mathit{\chi }}_r{\mathrm{\Gamma }}_r{\displaystyle \frac{{\mathit{\varphi }}_r\left(a_{\mathrm{2}}l\right)}{m}}.\end{array}\end{array}$$

Furthermore, the output power of the harvester can be expressed as

$$\begin{array}{}\text{(32)}& \text{Po}\left(t\right)={\displaystyle \frac{(v{e}^{I\mathit{\varpi }t}{)}^{\mathrm{2}}}{R}}.\end{array}$$

Next, we suppose that both transverse rheonomic constraints are quite distinct in amplitude and phase. Let $\mathit{\epsilon }=f_{\mathrm{2}}/f_{\mathrm{1}}$ be the ratio of the amplitudes and $\mathit{\beta }\in [\mathrm{0},\mathrm{2}\mathit{\pi })$ be the phase difference, respectively. Then the two transverse rheonomic constraints can be written as

$$\begin{array}{}\text{(33)}& P_{\mathrm{2}}\left(t\right)=\mathit{\epsilon }{e}^{\mathit{\beta }I}{P}_{\mathrm{1}}\left(t\right).\end{array}$$

The power frequency response function (FRF, Erturk and Inman, 2009) of the harvester is

$$\begin{array}{}\text{(34)}& \begin{array}{ll}{P}^{\prime }o=\frac{\text{Po}\left(t\right)}{{\left(-f_{\mathrm{1}}{\mathit{\varpi }}^{\mathrm{2}}{e}^{I\mathit{\varpi }t}\right)}^{\mathrm{2}}}& (f_{\mathrm{1}}\ne \mathrm{0})\\ P^{\prime }o={\displaystyle \frac{\text{Po}\left(t\right)}{{\left(-f_{\mathrm{2}}{\mathit{\varpi }}^{\mathrm{2}}{e}^{I\mathit{\varpi }t}\right)}^{\mathrm{2}}}}& (f_{\mathrm{1}}=\mathrm{0}).\end{array}\end{array}$$

3 The performance of the harvester and discussions

In this section, we use the formulas derived in the last section to analyze the influence of the main parameters of the harvester, ε and β, as expressed in Eq. (33), on the harvesting characteristics. In the numerical calculation, we select the first eight modal functions.

To analyze the performance of the harvester at low frequencies, we expand Eq. (31) into a power series at ϖ=0; then, we have

$$\begin{array}{}\text{(35)}& \begin{array}{rl}v& =\mathit{\varpi }I{f}_{\mathrm{1}}R\left(-\mathrm{1}+e^{I\mathit{\beta }}\mathit{\epsilon }\right)\\ & \cdot {\displaystyle \frac{{\mathit{\chi }}_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{4}}{\mathit{\lambda }}_{\mathrm{4}}^{\mathrm{4}}\left[{\mathit{\varphi }}_{\mathrm{2}}\left(\mathrm{0}\right)-{\mathit{\varphi }}_{\mathrm{2}}\left(a_{\mathrm{2}}l\right)\right]+{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{4}}{\mathit{\lambda }}_{\mathrm{4}}^{\mathrm{4}}{\mathit{\chi }}_{\mathrm{3}}\left[{\mathit{\varphi }}_{\mathrm{3}}\left(\mathrm{0}\right)-{\mathit{\varphi }}_{\mathrm{3}}\left(a_{\mathrm{2}}l\right)\right]+\mathrm{\dots }}{{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{4}}{\mathit{\lambda }}_{\mathrm{4}}^{\mathrm{4}}{\left[{\mathit{\varphi }}_{\mathrm{2}}\left(\mathrm{0}\right)-{\mathit{\varphi }}_{\mathrm{2}}\left(a_{\mathrm{2}}l\right)\right]}^{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{4}}{\mathit{\lambda }}_{\mathrm{4}}^{\mathrm{4}}{\left[{\mathit{\varphi }}_{\mathrm{3}}\left(\mathrm{0}\right)-{\mathit{\varphi }}_{\mathrm{3}}\left(a_{\mathrm{2}}l\right)\right]}^{\mathrm{2}}+\mathrm{\dots }}}\\ & +{\mathit{\varpi }}^{\mathrm{2}}{f}_{\mathrm{1}}R\left(-\mathrm{1}+e^{I\mathit{\beta }}\mathit{\epsilon }\right)\\ & \cdot {\displaystyle \frac{l^{\mathrm{4}}mR{\mathit{\chi }}_{\mathrm{3}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{4}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{4}}{\left[{\mathit{\varphi }}_{\mathrm{2}}\left(\mathrm{0}\right)-{\mathit{\varphi }}_{\mathrm{2}}\left(a_{\mathrm{2}}l\right)\right]}^{\mathrm{2}}+\mathrm{\dots }}{\text{YI}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{4}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{4}}{\left({\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{4}}{\left[{\mathit{\varphi }}_{\mathrm{2}}\left(\mathrm{0}\right)-{\mathit{\varphi }}_{\mathrm{2}}\left(a_{\mathrm{2}}l\right)\right]}^{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{4}}{\left[{\mathit{\varphi }}_{\mathrm{3}}\left(\mathrm{0}\right)-{\mathit{\varphi }}_{\mathrm{3}}\left(a_{\mathrm{2}}l\right)\right]}^{\mathrm{2}}\right)}^{\mathrm{2}}+\mathrm{\dots }}}\\ & +{\mathit{\varpi }}^{\mathrm{3}}I{f}_{\mathrm{1}}R\\ & \cdot {\displaystyle \frac{\begin{array}{c}<l^{\mathrm{4}}mYI{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{6}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}{\left[{\mathit{\varphi }}_{\mathrm{3}}\left(\mathrm{0}\right)-{\mathit{\varphi }}_{\mathrm{3}}\left(a_{\mathrm{2}}l\right)\right]}^{\mathrm{2}}\cdot \mathit{\{}lm{\mathit{\chi }}_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{4}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{4}}\\ \left[{\mathrm{e}}^{i\mathit{\beta }}\mathit{ϵ}{\mathit{\varphi }}_{\mathrm{3}}\left(\mathrm{0}\right)-{\mathit{\varphi }}_{\mathrm{3}}\left(a_{\mathrm{2}}l\right)\right]\left({\mathit{\varphi }}_{\mathrm{2}}\left(a_{\mathrm{2}}l\right){\mathit{\varphi }}_{\mathrm{3}}\left(\mathrm{0}\right)-{\mathit{\varphi }}_{\mathrm{2}}\left(\mathrm{0}\right){\mathit{\varphi }}_{\mathrm{3}}\left(a_{\mathrm{2}}l\right)\right)\mathit{\}}>\end{array}}{{\text{YI}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{6}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{6}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{8}}{\left[{\mathit{\varphi }}_{\mathrm{2}}\left(\mathrm{0}\right)-{\mathit{\varphi }}_{\mathrm{2}}\left(a_{\mathrm{2}}l\right)\right]}^{\mathrm{4}}+\mathrm{\dots }}}\\ & +\mathrm{\dots }+{\mathit{\varpi }}^{\mathrm{3}}\left(-\mathrm{1}+e^{I\mathit{\beta }}\mathit{\epsilon }\right)\left(\mathrm{\dots }\right)+O\left({\mathit{\varpi }}^{\mathrm{4}}\right).\end{array}\end{array}$$

From Eq. (35), it can be found that the output voltage v approaches zero as the frequency ϖ approaches zero. But the decreasing rate is changed with the variation of the parameters of ε and β. If the two excitations are different in the phase β≠0 or β=0 but ε≠1, then $\left(-\mathrm{1}+e^{I\mathit{\beta }}\mathit{\epsilon }\right)\ne \mathrm{0}$. The term O(ϖ) dominates the changing rate of the output voltage, v, of the energy harvester; the term O(ϖ²) dominates the changing rate of the output power, Po(t), in Eq. (32); and further, the term O(ϖ⁻²) dominates the changing rate of power FRF , P^′o, in Eq. (34), where “−” means that the value of power FRF will increase as the frequency ϖ approaches zero.

If the two excitations are the same, i.e., both β=0 and ε=1, the coefficients of ϖ and ϖ² in Eq. (35) are equal to zero, i.e., $\left(-\mathrm{1}+{\mathrm{e}}^{I\mathit{\beta }}\mathit{\epsilon }\right)=\mathrm{0}$, while the coefficient of ϖ³ is not equal to zero. Then the term O(ϖ³) dominates the changing rate of the output voltage v of the energy harvester, the term O(ϖ⁶) dominates the changing rate of output power Po(t), and the term O(ϖ⁺²) dominates the changing rate of power FRF P^′o, where “+” means that the value of power FRF will decrease as frequency ϖ approaches zero. Therefore, the decreasing rate of the power FRF of the harvester subjected to two different excitations is ϖ⁴ times faster than that of the harvester subjected to the two same excitations. In the latter case, the beam vibrates under the excitation of its inertial force, similar to the case of the normal piezoelectric energy harvester (Erturk and Inman, 2011), in which the metallic beam is fixed to the vibrator on one point, and its output voltage can be expanded into power series as

$$\begin{array}{}\text{(36)}& \begin{array}{rl}{v}_{\mathrm{normal}}& =-IR{f}_{\mathrm{1}}\left[\sum _{i=\mathrm{1}}^{\mathrm{\infty }}{\displaystyle \frac{l^{\mathrm{4}}{m}^{\mathrm{2}}{\mathit{\chi }}_i}{\text{YI}{\mathit{\lambda }}_i^{\mathrm{4}}}}\underset{\mathrm{0}}{\overset{l}{\int }}{\mathit{\varphi }}_i\left(x\right)\mathrm{d}x\right]{\mathit{\varpi }}^{\mathrm{3}}\\ & +O\left({\mathit{\varpi }}^{\mathrm{4}}\right).\end{array}\end{array}$$

In Eqs. (35) and (36), the lowest power of ϖ is 1 and 3, respectively. As a whole, the characteristic of the harvester subjected to two different excitations is quite different from that of the normal piezoelectric energy harvester or that of the harvester subjected to the two same excitations. To present intuitive comparisons in the characteristic of the harvesters, we provide some results of a numerical analysis in which the parameters of the beam, piezoelectric material, and electronic circuit are the same as that in the classical experiment of an energy harvester with PZT-5A (Erturk and Inman, 2011).

Figure 4a presents the output power $\left|\text{Po}\left(t\right)/e^{\mathrm{2}I\mathit{\varpi }t}\right|\times {\mathit{\varpi }}^{-n}$ $n=\mathrm{1},\mathrm{2},\mathrm{3}$ generated by the harvester subjected to two different excitations. In this figure, as $n=\mathrm{1},\mathrm{3}$, the value of $\left|\text{Po}\left(t\right)/e^{\mathrm{2}I\mathit{\varpi }t}\right|\times {\mathit{\varpi }}^{-n}$ decreases or increases, respectively, as the frequency ϖ approaches zero. However, as n=2, the value of $\left|\text{Po}\left(t\right)/e^{\mathrm{2}I\mathit{\varpi }t}\right|\times {\mathit{\varpi }}^{-\mathrm{2}}$ is a constant concerning the variations of frequency ϖ. It can be inferred that the value of $\left|\text{Po}\left(t\right)/e^{\mathrm{2}I\mathit{\varpi }t}\right|$ is equivalent infinitesimal to ϖ². Similarly, Fig. 4b and c show that $\left|\text{Po}\left(t\right)/e^{\mathrm{2}I\mathit{\varpi }t}\right|$ is equivalent infinitesimal to ϖ⁶ for the harvester subjected to the two same excitations and the normal harvester, respectively. It is confirmed that for the three kinds of harvesters, the output powers approach zero as the frequency ϖ approaches zero, but the decreasing rate of the harvester subjected to the two different excitations is slower ϖ⁴ times than that of the normal harvester and the harvester subjected to the two same excitations. It means that at the low frequency, the output power of the harvester subjected to the two different excitations is certainly higher than that of the two other harvesters. In other words, as the frequency ϖ approaches zero, the harvester subjected to the two different excitations performs better than the two other harvesters.

Figure 4 $\left|\text{Po}\left(t\right)/e^{\mathrm{2}I\mathit{\varpi }t}\right|\times {\mathit{\varpi }}^{-n}(n\in \mathbb{N})$ of the three kinds of piezoelectric energy harvesters.

Undoubtedly, such good performance is influenced by the parameters of the system.

Figure 5a shows the variations of power FRF P^′o concerning the phase difference β and the frequency ϖ of the excitations as parameters R=10⁵Ω, a₂=0.5, and ε=1. From this figure, it can be seen that the variation trends of the power FRF for the increasing frequency are variant when β is changed. Figure 5b shows the curves of the power FRF as β=0, $\mathit{\beta }=\mathit{\pi }/\mathrm{2}$, and β=π, along with the curve of the normal piezoelectric energy harvester reproduced from research work by Erturk and Inman (2011).

Figure 5 Variation of power FRF for frequency and phase difference (R=10⁵Ω, a₂=0.5, and ε=1).

In Fig. 5b, the shapes of the curves for the harvesters subjected to the two different excitations are not the same as that for the normal piezoelectric energy harvesters and those with the two same excitations. For harvesters subjected to two different excitations such as $\mathit{\beta }=\mathit{\pi }/\mathrm{2}$ or β=π, with the frequency increasing from 0 to 3000 Hz, the curve firstly tends downward, then upward to a local maximum around 1400 Hz, and after that downward continually. However, the curve for the normal piezoelectric energy harvester tends firstly upward, reaches the global maximum around 500Hz, then downward continually. The curve for the harvester subjected to the two same excitations tends firstly upward to one of its local maximums, then downward followed by upward to reach another local maximum, and afterwards it goes downward continually.

We are interested in the behaviors of the harvesters at low frequency. As shown in Fig. 5b, in the low-frequency interval (0–30 Hz), for the harvester subjected to two different excitations, i.e., $\mathit{\beta }=\mathit{\pi }/\mathrm{2}$ or β=π, the values of the power FRF are monotonically decreasing with the frequency increase of excitation. In other words, the lower the frequency of excitation, the better the performance of the harvester. This indeed is an unusual characteristic of the piezoelectric energy harvester. And the minimum value is $\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{6}}$ (W g⁻²) at 30 Hz. While in the same interval, the values of the power FRF for both the normal piezoelectric energy harvester and the harvester subjected to the two same excitations are monotonically increasing. And the former is larger than the latter; at a frequency of 30 Hz, the former is $\mathrm{5.8}\times {\mathrm{10}}^{-\mathrm{11}}$ (W g⁻²), much lower than that of the harvester subjected to the two different excitations.

Now we turn to discuss the effects of the parameter ε on power FRF at a specific low frequency. We set the frequency to 2 Hz. Figure 6 shows the variation of the power FRF for a variation of ε for some phase difference β as R=10⁵Ω, a₂=0.5, and ϖ=2 Hz. It can be seen that for the phase differences β>0.4π, the values of the power FRF will monotonically increase with the increase of ε, while for β<0.4π, with the increase of ε, the values of the power FRF firstly decrease and then increase. The curve for β=0 is the lowest one, and the minimum value, $\mathrm{1.12}\times {\mathrm{10}}^{-\mathrm{14}}$ (W g⁻²), occurs when ε=1, i.e., the harvester subjected to the two same excitations.

Figure 6 Power FRF with respect to the variation of ε (R=10⁵Ω, a₂=0.5, ϖ=2 Hz).

It must be noted that “power FRF” is not the actual power. The analysis of “power FRF” shows the transformation ability of the harvester and the classic harvester to the excitation motion. We will use an example to show how this difference in transformation ability is reflected in the actual power.

To compare the performance of the harvester subjected to the two different excitations with that of a normal piezoelectric energy harvester, we set f₁=0.01l, ε=1, β=0.3π, and a₂=0.5. The output power Po(t) of the harvester together with that of the normal piezoelectric energy harvester are depicted in Fig. 7, and four frequencies at which the comparisons are made are listed in Table 1. It should be pointed out that resistances corresponding to the maximum output power are not the same according to Eq. (32).

Figure 7 The comparison under low-frequency excitation (a₂=0.5, ε=1, and β=0.3π).

Table 1 Comparisons of maximum output power at four frequencies.

As shown in Table 1, the output power of the harvester subjected to the two excitations as ε=1 and β=0.3π is much higher than that of the normal piezoelectric energy harvester: at 30 Hz, the former is ×10⁶ times higher than the latter; at 10 Hz, it is 8×10⁶ times; and at 1 Hz, it is 5×10¹² times. In short, at lower frequencies, the energy harvester subjected to two different excitations will generate more times power than a normal harvester.

Figure 8 Power FRF (R=10⁵Ω and ϖ=2 Hz).

Finally, we discuss the influence of the parameter a₂ on the performance of the harvester subjected to the two different excitations. We set $\mathit{\beta }=\left\{\mathrm{0},\mathrm{0.2}\mathit{\pi },\mathrm{0.4}\mathit{\pi },\mathrm{0.6}\mathit{\pi }\right\}$, $\mathit{\epsilon }=\left\{\mathrm{0.5},\mathrm{1}\right\}$, R=10⁵Ω, and ϖ=2 Hz. With Eq. (34), the power FRF P^′o for the variation of parameters a₂ are shown in Fig. 8. It can be seen from Fig. 8a that the trends of all the curves for β≠0 are almost the same: a smaller a₂, in general, results in greater output power. But there is an exception: the curve for the case when β=0 and ε=1, i.e., the harvester subjected to the two same excitations, not only for its trend but for the value of the output power which is lower. However, in Fig. 8b, for β=0, these distinctions do not exist. It also indicates that differences between the two excitations are crucial for obtaining a better performance of the energy harvester with two excitations.

4 Conclusion

For harvesting energy from human body motion, it is proposed in this paper to connect two parts of the human body with different transversal reciprocating excitations using a metallic beam attached with piezoelectric ceramics layers to form the energy harvester. In this paper, we emphasize the analysis of the characteristics of the energy harvester subjected to two different excitations.

Employing the Hamiltonian principle, the governing equation of the piezoelectric energy harvester is established. Assuming that both the two excitations are periodic functions, the closed-form solutions of the dynamics equation are obtained. Utilizing the closed-form solutions, the influences of the differences in the amplitude and phase difference of the excitations on the performance of the harvester are analyzed together with numerical examples.

The results show that as long as there exists some differences between the two excitations, either in amplitude and/or in phase difference, the energy harvester can generate more output power at low frequencies, and the larger the differences are, the more the output power will be generated. Contrasting the normal piezoelectric energy harvester and the harvester subjected to two of the same excitations, the harvester subjected to two different excitations exhibits a characteristic: the output power increases with the decrease of frequency. This unusual feature enables the proposed harvester to be quite suitable to harvest energy from human body motions.

In this paper, the analysis is based on a mathematical model. In future research, a prototype device will be built to validate the results experimentally. Additionally, the same kind of harvester subjected to more than two excitations with different excitations will be studied.

Appendix A: Nomenclature

l	length of beam
a1, a2	position parameters of vibrator 1 and
	vibrator 2
t	time
Ph(t)	the hth displacement functions
${\stackrel{\mathrm{̃}}{P}}_h\left(t\right)$	the hth constraint equation
w(x,t)	transverse displacement of the beam
ϕr(x)	the rth order modal function of the beam
qr(t)	the rth modal coordinate
I	$\sqrt{-\mathrm{1}}$
ρs	density of the metal layer
ρp	density of the piezoelectric layer
hs	thickness of the metal layer
hp	thickness of the piezoelectric layers
b	width of the beam
c11	elastic constant of the piezoelectric layers
E	elastic modulus of the metal layer
m	the mass per unit length of the beam
YI	bending stiffness
Ip, Is, Ib	cross-sectional area moment of inertia of
	the metal layer, inertia of piezoelectric
	layers, the whole beam
ε33	dielectric constant
e31	piezoelectric constant
Qs	virtual work of air damping and
	material damping
cs	average material damping coefficient
ca	air damping coefficient
QR	virtual work of the load resistance
v(t)	voltage
κ(t)	Lagrange multiplier
ζr	rth order damping ratio
ωr	natural frequency of the beam
Po(t)	output power
P′o	power FRF

fr	amplitude of the rth vibrators
ϖ	frequency of the vibrators
χr	electromechanical coupling coefficient
Cp	capacitance
Hr	amplitude of rth modal coordinates
κi	amplitude of ith Lagrange multipliers
Γr	grouped term 1
ξij	grouped term 2
β	phase difference

Appendix B: Derivation of Eq. (11)

Let

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}& L\left(w^{\prime \prime },\dot{w},v,{\mathit{\kappa }}_{\mathrm{1}},{\mathit{\kappa }}_{\mathrm{2}}\right)={\displaystyle \frac{m}{\mathrm{2}}}{\dot{w}}^{\mathrm{2}}-{\displaystyle \frac{\text{YI}}{\mathrm{2}}}{w^{\prime \prime }}^{\mathrm{2}}+{\displaystyle \frac{{\mathit{\alpha }}^{\mathrm{2}}{\mathit{\epsilon }}_{\mathrm{33}}b{v}^{\mathrm{2}}}{\mathrm{4}{h}_{\mathrm{p}}}}\\ & -\mathit{\alpha }{e}_{\mathrm{31}}{h}_{\mathrm{pc}}bv{w}^{\prime \prime }+{\mathit{\kappa }}_{\mathrm{1}}{\stackrel{\mathrm{̃}}{p}}_{\mathrm{1}}\left(t\right)w\mathit{\delta }i\left(x-a_{\mathrm{1}}l\right)\\ & +{\mathit{\kappa }}_{\mathrm{2}}{\stackrel{\mathrm{̃}}{p}}_{\mathrm{2}}\left(t\right)w\mathit{\delta }i\left(x-a_{\mathrm{2}}l\right),\end{array}\end{array}$$

where δi is the Dirac Dirichlet function. Using Eq. (B1) in Eq. (10), it becomes

$$\begin{array}{}\text{(B2)}& \underset{\mathrm{0}}{\overset{t}{\int }}\underset{\mathrm{0}}{\overset{l}{\int }}\mathit{\delta }L\mathrm{d}x\mathrm{d}t+\underset{\mathrm{0}}{\overset{t}{\int }}{Q}_s\mathrm{d}t+\underset{\mathrm{0}}{\overset{t}{\int }}{Q}_R\mathrm{d}t=\mathrm{0}.\end{array}$$

Integrate Eq. (B2) by parts, and it becomes

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}& \underset{\mathrm{0}}{\overset{t}{\int }}\underset{\mathrm{0}}{\overset{l}{\int }}\mathit{\delta }L\mathrm{d}x\mathrm{d}t+\underset{\mathrm{0}}{\overset{t}{\int }}{Q}_s\mathrm{d}t+\underset{\mathrm{0}}{\overset{t}{\int }}{Q}_R\mathrm{d}t\\ & =\underset{\mathrm{0}}{\overset{t}{\int }}\underset{\mathrm{0}}{\overset{l}{\int }}[\left({\displaystyle \frac{d^{\mathrm{2}}}{\mathrm{d}{x}^{\mathrm{2}}}}{\displaystyle \frac{\partial L}{\partial w^{\prime \prime }}}-{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\displaystyle \frac{\partial L}{\dot{\partial w}}}-c_s{I}_{\mathrm{b}}{\dot{w}}^{\prime \prime \prime \prime }-c_{\mathrm{a}}\dot{w}\right)\mathit{\delta }w\\ & +\left({\displaystyle \frac{\partial L}{\partial v}}\right)\mathit{\delta }v+\left({\displaystyle \frac{\partial L}{\partial {\mathit{\kappa }}_{\mathrm{1}}}}\right)\mathit{\delta }{\mathit{\kappa }}_{\mathrm{1}}+\left({\displaystyle \frac{\partial L}{\partial {\mathit{\kappa }}_{\mathrm{2}}}}\right)\mathit{\delta }{\mathit{\kappa }}_{\mathrm{2}}]\mathrm{d}x\mathrm{d}t.\end{array}\end{array}$$

In Eq. (B3), the coefficients of δw, δv, δκ₁, and δκ₂ equal 0. Notice that $v=v\left[H\left(x\right)-H(x-l)\right]$, where H is the Heaviside function. Using Eq. (B1) in Eq. (B3), we get

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}m{\dot{w}}^{\prime \prime }& +\text{YI}{w}^{\prime \prime \prime \prime }+c_s{I}_{\mathrm{b}}{\dot{w}}^{\prime \prime \prime \prime }+c_{\mathrm{a}}\dot{w}\\ & -{\mathit{\chi }}_{r}v\left[\mathit{\delta }{i}^{\prime }\left(x\right)-\mathit{\delta }{i}^{\prime }(x-l)\right]\\ & =\sum _{i=\mathrm{1}}^{\mathrm{2}}{\mathit{\kappa }}_i\left(t\right)\left(w(x,t)\mathit{\delta }i\left(x-a_{i}l\right)\right),\end{array}\end{array}$$

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}& C_{p}v\left(t\right)+{\displaystyle \frac{v\left(t\right)}{R}}+e_{\mathrm{31}}{\displaystyle \frac{h_{\mathrm{p}}+h_s}{\mathrm{2}}}b\underset{\mathrm{0}}{\overset{l}{\int }}{\dot{w}}^{\prime \prime }\mathrm{d}x=\mathrm{0}\\ & \left\{r=\mathrm{1},\mathrm{\dots },\mathrm{\infty }\right\},\end{array}\end{array}$$

$$\begin{array}{}\text{(B6)}& {\stackrel{\mathrm{̃}}{p}}_{\mathrm{1}}\left(t\right)={\stackrel{\mathrm{̃}}{p}}_{\mathrm{2}}\left(t\right)=\mathrm{0}.\end{array}$$

Use Eq. (2) in Eqs. (B4)–(B6) and multiply ϕ_r(x) at both sides of the equations. Integrate equations from 0 to l, and using ${\int }_{\mathrm{0}}^{l}m{\mathit{\varphi }}_s\left(x\right){\mathit{\varphi }}_r\left(x\right)={\mathit{\delta }}_{rs}$, one can obtain Eq. (11).