1 Introduction

The energy harvesting from ambient environmental vibration sources has been widely studied over the past decades for its promising energy density and ease of acquisition (Beeby et al., 2006; MacCurdy et al., 2008; Yang et al., 2010; Wang et al., 2012). Piezoelectric and electromagnetic are the most effective way of turning mechanical vibrations into electrical energy (Roundy et al., 2004; Wacharasindhu and Kwon, 2008; Tadesse et al., 2009; Yuan et al., 2009; Wang et al., 2013). To enhance the energy conversion efficiency, researchers proposed a kind of hybrid energy harvesting technology by coupling the piezoelectric and the electromagnetic mechanism together.

The hybrid energy harvester has received great attention for its important advantages, such as being clean, stable, having high energy density, and a wide frequency band. For example, a general spring-mass-damper model was presented by Williams and Yates (1996), where the electrical energy harvested was equivalent to the energy dissipated in the electrical damping. The spring-mass-damper model was widely used by simplifying the effect of two complex mechanisms into viscous damping. Challa et al. (2009) generated 332 μW at 21.6 Hz from a cantilever beam hybrid energy harvester. This harvester produced a 30% increase in power output compared to 257 μW and 244 μW generated from the stand-alone piezoelectric and electromagnetic energy harvesting devices, respectively. Wischke et al. (2009) presented a model of a hybrid generator which was based on a beam suspended on both sides. Gonsalez et al. (2010) conducted a model considering the coupling influences between the piezoelectric and electromagnetic systems. They found that the mechanical system was strongly influenced when an electrical circuit was coupled to the transducer.

From the above descriptions, it can be found that the previous devices can have effective outputs only when they work near their optimal frequencies. However, it is difficult to achieve a high power output from a low frequency vibration in the environment. Therefore, a new hybrid energy harvester is presented to broaden the frequency bandwidth and coordinate ambient low-frequency vibration. Firstly, the mathematical model of the hybrid energy harvester will be established. Secondly, the corresponding numerical simulation will be performed. Finally, the validity of the mathematical model will be verified through a finite element simulation and supporting experiments.

2 Mathematical modeling

Fig. 1 illustrates a schematic of the newly presented hybrid energy harvester. This harvester consists of a Pb(Zr,Ti)O3 (PZT) bimorph cantilever beam mounted on the vibrating shaker and a magnetic levitation system fixed on the PZT beam. The magnetic levitation system used two outer magnets that were mechanically attached to a support. A center magnet was placed between the two outer magnets. The magnetic poles were oriented to repel the center magnet. Thus, the center magnet was suspended by a nonlinear restoring force (Mann and Owens, 2010). Spiral coils were wound around the outer casing to transform the mechanical energy of the center magnet to electrical energy.

Fig. 1
figure 1

Schematic of a hybrid energy harvester

Full size image

2.1 Spring-mass-damper model

In this section, the hybrid energy harvester device can be simplified as a spring-mass-damper system (Fig. 2), consisting of a mechanical damping C m (which accounts for the energy losses due to structural and viscous damping), a piezoelectric damping C p, and an electrical damping C e.

Fig. 2
figure 2

Schematic of the spring-mass-damper system

Full size image

Here the vibrating structure is represented as a spring with a stiffness K and a mass m. A source vibration of w(t)=Acos(ωt) (where A and ω are source vibration amplitude and frequency, respectively) generates a relative deflection z(t) in the vibrating structure. A differential equation using the above-described model is given as

$$\ddot z\left( t \right) + 2\zeta {\omega _n}\dot z\left( t \right) + \omega _n^2z\left( t \right) = A{\omega ^2}{\text{cos(}}\omega t{\text{),}}$$
((1))

where ω n is the undamped natural frequency of the vibrating structure, \({\omega _n} = \sqrt {K/m}\), and ζ is the total damping ratio, ζ = C / (2mω n ).

The solution of Eq. (1) can be divided into two parts, i.e., the transient solution and the steady-state solution. The transient solution is caused by the beginning unstable vibration or disturbance of the system itself. So the transient solution is ignored, by retaining the steady-state solution.

$$z(t) = \frac{{A{{(\omega /{\omega _n})}^2}}}{{\sqrt {{{\left( {1 - {{\left( {\omega /{\omega _n}} \right)}^2}} \right)}^2} + {{\left( {2\zeta \left( {\omega /{\omega _n}} \right)} \right)}^2}} }}{\text{cos}}\left( {\omega t - \theta } \right),$$
((2))

where θ is the phase difference and it can be ignored for having almost no effect on the power output.

In the hybrid power generation device, z(t) represents the vibration of the piezoelectric beam, that is, the electromagnetic energy harvester device input vibration.

2.2 Power output from a stand-alone electromagnetic energy harvester

Before starting to build the mathematical model of the magnetic suspension energy harvester, one should know the expression of the acting force between a fixed permanent magnet and a vibrating magnet. There are some empirical formulas for solving the above acting force. However, these empirical formulas can only provide an approximate result. To establish a more accurate model for a magnetic suspension system, the software MAXWELL was used to obtain the acting force, and to efficiently and accurately solve the electromagnetic field simulation directly from the actual model.

Fig. 3 illustrates the relationship between the restoring force and the separation distance of two magnets simulated by MAXWELL. A polynomial curve-fitting algorithm in MATLAB was employed based on the above relationship. It can be found that the cubic polynomial fitting can effectively meet the accuracy requirements (Fig. 3).

Fig. 3
figure 3

Relationship between the restoring force and the separation distance of two magnets and the cubic polynomial fitting

Full size image

Therefore, according to the fitting results, the force between the two magnets can be written as

$${F_{\text{e}}} = {b_3}{s^3} + {b_2}{s^2} + {b_1}s + {b_0},$$
((3))

where s is the separation distance between two magnets, and the coefficients b 0, b 1, b 2, and b 3 are obtained using MATLAB

To determine the restoring forces of the suspension electromagnetic system, Fig. 4 illustrates a schematic diagram of the magnetic levitation system.

Fig. 4
figure 6

Schematic of the magnetic levitation system

Full size image

Hence, when the suspended magnet deviates as displacement y from the equilibrium position, the total magnetic restoring force can be calculated as

$$\begin{array}{*{20}{c}} {F\left( y \right) = {b_3}\left[ {{{\left( {h - y} \right)}^3} - {{\left( {h + y} \right)}^3}} \right]} \\ { + {b_2}\left[ {{{\left( {h - y} \right)}^2} - {{\left( {h + y} \right)}^2}} \right] + 2{b_1}y} \\ { = 2{b_3}{y^3} - \left( {6{h^2}{b_3} + 4h{b_2} + 2{b_1}} \right)y} \\ { = {k_3}{y^3} + {k_1}y,} \end{array}$$
((4))

where h stands for the separation distance between the fixed magnets and the suspension magnet, and the coefficients are k 1=39.5 N/m, k 3=97320 N/m3.

According to (Mann and Sims, 2009), a governing equation is given as

$${m_0}\ddot x + {C_{{\text{em}}}}\left( {\dot x - \dot z} \right) + {k_1}\left( {x - z} \right) + {k_3}{\left( {x - z} \right)^3} + \alpha i = - mg,$$
((5))

where m 0 is the magnet mass; g is the gravitational constant; C em is the damping coefficient used to describing the mechanical damping; x is the motion of the center magnet; i is the electrical current; and α is the electromechanical coupling coefficient. This term is commonly expressed as α=NBL, where N is the number of coil turns, B is the average magnetic field strength, and L is the coil length.

According to (Stephen, 2006), the coil magnetic force is assumed to be proportional to the relative velocity between the center magnet and outer housing, \(i = \frac{{\alpha \left( {\dot x - \dot z} \right)}}{{{R_{\text{e}}} + {R_1}}}\), where R e is the internal resistance of the coil, and R 1 is the resistance of the external load. To obtain the power output from a stand-alone electromagnetic energy harvester, the input acceleration A 0 remains unchanged, so \(z\left( t \right) = \frac{{{A_0}}}{{{\omega ^2}}}{\text{cos}}\left( {\omega t} \right)\). C e describes the combined mechanical and electrical damping of the magnetic levitation system, C e=C et+C em, where C et is the expression for the electrical damping coefficient, \({C_{{\text{et}}}} = \frac{{{\alpha^2}}}{{{R_{\text{e}}} + {R_1}}}\). Substituting y(t)=x(t)−z(t) into Eq. (5), we can obtain a governing equation written in terms of the relative displacement between the center magnet and outer housing.

$$\ddot y + 2\mu \varepsilon \dot y + \omega _0^2y + \varepsilon \beta {y^3} = {F_0} + 2\varepsilon {F_1}{\text{ cos}}\left( {\omega t} \right),$$
((6))

where the coefficients of Eq. (5) have been rewritten as A 0=2εF 1, \({A_0} = 2\varepsilon {F_1}\), \({\beta _0} = \varepsilon \beta = \tfrac{{{k_3}}}{m}\), \(\omega _0^2 = \tfrac{{{k_1}}}{m}\), \( - g = {F_0}\), \(\frac{{{C_{\text{e}}}}}{m} = 2\mu \varepsilon\) and ε has been introduced as a parameter (Mann and Sims, 2009).

Eq. (6) is a second-order nonlinear differential equation. The following is the process of solving this equation using a multiple time scale. The assumed solution to Eq. (6) is written as a first-order expansion.

$$y\left( t \right) = {y_0}\left( t \right) + \varepsilon {y_1}\left( t \right).$$
((7))

The independent time scales are defined as τ0= t, τ1=εt. It follows that the derivatives with respect to time become partial derivatives with respect to the corresponding time scale after setting \({D_0}y = \frac{{\partial y}}{{\partial {\tau _0}}}\) and \({D_1}y = \frac{{\partial y}}{{\partial {\tau _1}}}\).

$$\frac{{{\text{d}}y}}{{{\text{d}}t}} = {D_0}{y_0} + \varepsilon \left( {{D_0}{y_0} + {D_1}{y_0}} \right) + {\varepsilon ^2}{D_1}{y_1},$$
((8))
$$\begin{array}{*{20}{l}} {\frac{{{{\text{d}}^2}y}}{{{\text{d}}{t^2}}} = D_0^2{y_0} + \varepsilon \left( {D_0^2{y_1} + {D_0}{D_1}{y_0}} \right)} \\ {{\text{ }} + {\varepsilon ^2}\left( {2{D_0}{D_1}{y_1} + {D_1}{y_0}} \right) + {\varepsilon ^3}{D_1}{y_1}} \end{array}$$
((9))

Substituting Eqs. (8) and (9) into Eq. (6) and separating the result into the orders of ε, the following two linear expressions can be obtained:

$${\varepsilon ^0}:D_0^2{\ddot y_0} + \omega _0^2{y_0} = {F_0},$$
((10))
$$\begin{array}{*{20}{l}} {{\varepsilon ^1}:D_0^2{{\ddot y}_0} + \omega _0^2{y_1} = 2{F_1}{\text{ cos}}\left( {\omega t} \right) - \beta y_0^3} \\ {{\text{ }} - 2\mu {D_0}{y_0} - 2{D_0}{D_1}{y_0}.} \end{array}$$
((11))

The solution of Eq. (10) can be written as

$${y_0} = A\left( {{\tau _1}} \right){{\text{e}}^{{\text{i}}{\omega _{\text{0}}}{\tau _{\text{0}}}}} + \bar A\left( {{\tau _1}} \right){{\text{e}}^{ - {\text{i}}{\omega _{\text{0}}}{\tau _{\text{0}}}}} + \frac{{{F_0}}}{{\omega _0^2}},$$
((12))

where A(τ1 and Ā(τ1) are complex conjugates that are functions of τ1.

After setting ω=ω0+εσ, and substituting Eq. (12) into Eq. (11), the following expression can be obtained:

((13))

where O(einω 0 τ 0) is the polynomial that n ε1, σ is called a detuning parameter for describing the nearness of ω and ω0. In the method of multiple time scale, the secular terms equal 0.

((14))

According to (Nayfeh and Balachandran, 2008), the polar forms, \(A\left( {{\tau _1}} \right) = \frac{1}{2}\alpha \left( {{\tau _1}} \right){{\text{e}}^{{\text{i}}\Phi \left( {{\tau _1}} \right)}}\) and \(\overline A \left( {{\tau _1}} \right) = \frac{1}{2}\alpha \left( {{\tau _1}} \right){{\text{e}}^{{\text{ - i}}\Phi \left( {{\tau _1}} \right)}}\), are now introduced into the above equation. After substituting γ=στ1−Φ=(ω−ω0)τ0−Φ and separating into the real and imaginary components, the resulting equations are

$${F_1}{\text{cos(}}\gamma {\text{)}} = - {\omega _0}{\alpha _0}\sigma + \frac{{3\beta }}{8}\left( {a_0^3 + 4{a_0}\frac{{F_0^2}}{{\omega _0^4}}} \right),$$
((15))
$${F_1}{\text{sin}}\gamma {\text{ = }}{\omega _{\text{0}}}\mu {\alpha _0}.$$
((16))

The equilibrium solutions of Eqs. (15) and (16) represent the steady-state periodic solutions of the system. They can be found by setting \(\frac{{\partial \alpha }}{{\partial {\tau _1}}} = 0\) and \(\frac{{\partial \gamma }}{{\partial {\tau _1}}} = 0\). Squaring and adding these relationships results into the expression of the frequency response for the nonlinear system, the initial parameters are introduced into the expression. Then, the modified frequency response equation is as follows:

$$\begin{array}{*{20}{l}} {\frac{{9k_3^2}}{{64{m_0}{k_1}}}\alpha _0^6 + \left[ {\frac{{9{m_0}{g^2}k_3^2}}{{8k_1^3}} + \frac{{3{k_3}}}{{4{m_0}}}\left( {1 - \sqrt {\frac{{{m_0}}}{{{k_1}}}} \omega } \right)} \right]\alpha _0^4} \\ { + \left[ {\frac{{9{m_0}^3{g^4}k_3^2}}{{4k_1^5}} + \frac{{3{m_0}{g^2}{k_3}}}{{k_1^2}}\left( {1 - \sqrt {\frac{{{m_0}}}{{{k_1}}}} \omega } \right)} \right.} \\ { + \left. {{{\left( {\sqrt {\frac{{{k_1}}}{{{m_0}}}} - \omega } \right)}^2} + \frac{{{c^2}}}{{4{m_0}^2}}} \right]\alpha _0^2 - \frac{{{m_0}A_0^2}}{{4{k_1}}} = 0.} \end{array}$$
((17))

The approximate solution of Eq. (6) is

$${y_0} = {a_0}({\text{cos(}}\omega t) - \gamma ) - \frac{g}{{\omega _0^2}}.$$
((18))

Thus, the output power of the electromagnetic energy harvester can be expressed as

$${P_{\text{e}}} = {\left[ {\frac{{\alpha \omega {\alpha _0}}}{{\sqrt 2 \left( {{R_{\text{e}}} + {R_1}} \right)}}} \right]^2}{R_1}.$$
((19))

2.3 Power output from a stand-alone piezoelectricity energy harvester

In this section, a double-chip piezoelectric cantilever is used. The open voltage U p generated by the piezoelectric effect depends on the stress σstruc and the material properties of PZT, which was given by Challa et al. (2009)

$${U_{\text{p}}} = \frac{{2{t_{\text{p}}}{d_{31}}{\sigma _{{\text{struc}}}}}}{\varepsilon },$$
((20))

where T p is the thickness of the piezoelectric ceramic, d 31 is the piezoelectric strain constant, and ε is the dielectric constant of the piezoelectric material.

As mentioned in the context, z(t) represents the vibration of the piezoelectric beam. After setting A1 and the acceleration of z(t), the equivalent force exerting on the piezoelectric beam can be expressed as

$$\begin{array}{*{20}{l}} {F = K \cdot z(t) = m\omega _n^2z(t)} \\ {{\text{ }} = \frac{{m{A_1}}}{{\sqrt {{{\left[ {1 - {{\left( {\omega /{\omega _n}} \right)}^2}} \right]}^2} + {{\left( {2\zeta \omega /{\omega _n}} \right)}^2}} }}{\text{cos(}}\omega t + \theta {\text{)}}{\text{.}}} \end{array}$$
((21))

The average bending stress of the piezoelectric film can be written as

$$\begin{array}{*{20}{l}} {{\sigma _{{\text{struc}}}} = \frac{1}{{{L_{\text{p}}}}}\int_{L - {L_{\text{p}}}}^L {\left( {\frac{1}{{{t_{\text{p}}}}}\int_{0.5{h_{\text{p}}}}^{0.5{h_{\text{p}}} + {t_{\text{p}}}} {\frac{{Fxy}}{I}{\text{d}}y} } \right){\text{d}}x} } \\ {{\text{ }} = \frac{{3m(2L - {L_{\text{p}}})({h_{\text{p}}} + {t_{\text{p}}}){A_1}{\text{ cos(}}\omega t + \theta {\text{)}}}}{{b{{({h_{\text{p}}} + 2{t_{\text{p}}})}^3}\sqrt {{{\left[ {1 - {{\left( {\omega /{\omega _n}} \right)}^2}} \right]}^2} + {{\left( {2\zeta \omega /{\omega _n}} \right)}^2}} }},} \end{array}$$
((22))

where L is the thickness of the piezoelectric ceramic, L p and h p are the length and thicknesses of the substrate, respectively, and I is the equivalent area moment, I = b(h p + 2T p)3 /12.

By inserting Eq. (22) into Eq. (20), the piezoelectric power output can be written as

$$\begin{array}{*{20}{l}} {{P_{\text{P}}}{\text{ = }}\frac{1}{T}\int_0^T {\frac{{U_{\text{p}}^2{R_2}}}{{{{({R_{\text{p}}} + {R_2})}^2}}}} {\text{d}}t} \\ { = \frac{{18{R_2}}}{{{{({R_{\text{p}}}{\text{ + }}{R_2})}^2}}}\frac{{t_{\text{p}}^2d_{31}^2{m^2}{{(2L - {L_{\text{p}}})}^2}{{({h_{\text{p}}} + {t_{\text{p}}})}^2}A_1^2}}{{{\varepsilon ^2}{b^2}{{({h_{\text{p}}} + 2{t_{\text{p}}})}^6}\left\{ {{{\left[ {1 - {{\left( {\frac{\omega }{{{\omega _n}}}} \right)}^2}} \right]}^2} + {{\left( {2\zeta \frac{\omega }{{{\omega _n}}}} \right)}^2}} \right\}}},} \end{array}$$
((23))

where R p is the internal resistance of the piezoelectric ceramic, and R 2 is the resistance of the external load.

2.4 Power output from a hybrid energy harvester

It is critical to note that this total power output is not equal to the sum of the individual stand-alone power outputs generated from the piezoelectric and electromagnetic techniques. The electromagnetic damping C e should be loaded to the piezoelectric system, and the vibration of the piezoelectric beam should be used as the input vibration of the electromagnetic system. After adding the coupled outputs of the piezoelectric system and the electromagnetic system, the total power output can be expressed as

$$P = {P_{\text{p}}} + {P_{\text{e}}} = {P_{\text{p}}} + \frac{{{\alpha ^2}{\omega ^2}\alpha _0^2{R_1}}}{{2{{({R_{\text{e}}} + {R_1})}^2}}}.$$
((24))

3 Numerical simulation

The materials of the magnet and the piezoelectric ceramic are NdFe35 and PZT-5H, respectively. Table 1 lists the dimension and material property parameters of the system.

Table 1 Dimension and material property parameters of the system
Full size table

The optimal external resistances of the electromagnetic and piezoelectric systems are 1 kΩ and 178 kΩ, respectively. In the case of the input acceleration of 5 m/s2, the power output of the hybrid energy harvester and individual energy harvester are shown in Fig. 5.

Fig. 5
figure 4

Comparison of the single device and the coupled power output (the input acceleration equals 5 m/s 2)

Full size image

The piezoelectric-electromagnetic vibration energy harvester has two output peak frequencies. The peak frequency of the electromagnetic occurs at 8.373 Hz, while the peak frequency of the piezoelectric occurs at 14.83 Hz. So it has a considerable power output in the range of 6 to 20 Hz, which is three times the bandwidth of the stand-alone, single harvesting mode devices. Furthermore, it increases the efficiency of energy conversion and the power output. The two peak power outputs are 11.4 mW and 21.6 mW, respectively.

4 Finite element simulation and experiments

The model of the electromagnetic system was established using MAXWELL (Fig. 6a) and the model of the piezoelectric system was established in ANSYS (Fig. 6b). Fig. 7 shows the experimental system. Fig. 8 shows the comparison of the hybrid output powers obtained from the numerical simulation, the finite element simulation, and the experiments. The experimental results show that the peak power outputs of the hybrid harvester are 7.2 mW and 16.4 mW at the resonance frequencies of 8.5 Hz and 16 Hz, respectively. Compared to the previous studies (Shen et al., 2008; Torah et al., 2008; Yang and Tang, 2009), this hybrid harvester effectively broadens the frequency bandwidth. The average error rates of the peak frequency of the numerical analysis and finite element simulation to the experimental results are very small, but that of the power output are slightly larger. The reasons for these slightly high error rates are related to the following aspects:

  1. 1.

    Simplified models are used in the mathematical modeling.

  2. 2.

    Even though the multiple time scale method can solve the governing equation, finally, the solution is an approximate one.

  3. 3.

    Lateral force is ignored in the mathematical modeling. Lateral force increases the mechanical damping, so it can affect the finite element result and numerical result.

Fig. 6
figure 5

Finite element models of the electromagnetic system (a) and piezoelectricity system (b)

Full size image
Fig. 7
figure 7

Experimental system

Full size image
Fig. 8
figure 8

Comparison of the numerical, the finite element simulation, and the experimental results

Full size image

5 Conclusions

In this paper, a new hybrid energy harvester using two power generation mechanisms (piezoelectric and electromagnetic) is presented. The theoretical power output models were developed based on the piezoelectric equation and the nonlinear theory. The experimental and finite element simulation results of this coupled device show a very good agreement with the numerical results of the developed mathematical model. It demonstrates that the theoretical model of the hybrid harvester is valid. The present work shows that the coupling of energy harvesting techniques is a simple and effective way to broaden the frequency bandwidth and increase the energy harvesting power output. Further, developing multiple power output peaks from different power generation mechanisms is a promising method to obtain an even better output effect.