Sage Journals: Discover world-class research

Abstract

This article describes an analytical model of the cantilever-type energy harvester based on Euler–Bernoulli’s beam theory. Starting from the Hamiltonian form of total energy equation, the bending mode shapes and electromechanical dynamic equations are derived. By solving the constitutive electromechanical dynamic equation, the frequency transfer function of output voltage and power can be obtained. Through a case study of a unimorph piezoelectric energy harvester, this analytical modeling method has been validated by the finite element method.

Keywords

Energy harvesting piezoelectric Euler–Bernoulli’s beam finite element method

Introduction

With the rapid advancement of information technology in the past 20 years, we are facing a new era of Internet of things (IOTs). However, this new technology is greatly dependent on various kinds of wireless sensor nodes (WSNs),¹ portable electronics,² and micro-electromechanical system (MEMS) devices.³ Currently, most of them are powered only by batteries, which may limit the whole system’s operational time and performance by periodically replacing or recharging these batteries. Since MEMS are scaled down in the range of 100 nm–1 mm and WSNs operate cyclically, the required power for one separate device is limited into micro-Watt level, which makes it possible to develop regenerative micro power supply such as solar cells, thermoelectric energy harvester, and mechanical vibration energy harvester⁴ to power these devices. However, compared with conventional battery technology, all of these complementary schemes have their corresponding advantages and disadvantages. Cook-Chennault et al.⁵ had ever reviewed these power solutions for portable MEMS devices, where the author made an emphasis on piezoelectric vibration generator for its structure simplicity, easy fabrication, high output voltage, and long lifespan duration. In this article, we mainly focus on the analytical modeling of piezoelectric vibration energy harvester.

The typical structure of piezoelectric scavenger is the rectangular cantilever beam,⁶ which mainly consists of an elastic substrate and a piezoelectric layer as shown in Figure 1, where one end of the beam is fixed with host structure and the other end is kept free. When the host structure is vertically vibrating, the fixed end will rigidly translate along with the fixture, and the free end will move with a lagged phase, by which mechanical strains are induced in the dynamic vibrating process. If the piezoelectric patch is coated with electrodes on top and bottom surfaces, electric charges across the electrode surfaces on the basis of direct piezoelectric effect can be induced from the mechanical strain.

Figure 1.

Schematic figure showing the dynamic motion of the unimorph piezoelectric energy harvester.

In this article, the structure is outlined as follows. In section “Analytical model of piezoelectric unimorph cantilever energy harvester,” the analytical model of piezoelectric unimorph cantilever energy harvester is described. Section “Multimode frequency response of the unimorph cantilever beam” provides multimode frequency response of the unimorph cantilever beam. In section “Case study of piezoelectric unimorph cantilever scavenger in bending mode,” both analytical and numerical results are compared to prove the validity of the analytical method. Beside this, the transient dynamic response of the device is also characterized. The conclusion remarks are finally made in section “Conclusion.”

Analytical model of piezoelectric unimorph cantilever energy harvester

When the composite cantilever beam vibrates with host structure vertically, arbitrary point P will translate from the original position to terminal position $P ″$ as denoted in Figure 1, where the vector $\vec{PP ″}$ is resolved into $\vec{PP'}$ and $\vec{P' P ″}$ . For the vector $\vec{PP'}$ , we assume that the cantilever beam moves rigidly with a displacement of $W_{base}$ . In terms of vector $\vec{P' P ″}$ , the cantilever along the longitudinal direction will deflect with a relative transverse displacement of $w_{rel}$ . Based on Euler–Bernoulli’s beam theory, the plane normal to the neutral plane before deformation is still perpendicular to it after a rotation. Therefore, the rotation angle $θ$ of point $P ″$ is defined as $\partial w_{rel} / \partial x$ . The position vectors of point P and $P ″$ in the OXZ coordinate system can be expressed as equations (1) and (2), respectively

\vec{OP} = x {\vec{e}}_{1} + z {\vec{e}}_{3}

(1)

\vec{OP ″} = (x - z \frac{\partial w_{rel}}{\partial x}) {\vec{e}}_{1} + (z + w_{base} + w_{rel}) {\vec{e}}_{3}

(2)

In order to characterize the dynamic motion with respect to reference coordinate, the absolute velocity of point $P ″$ can be derived based on equations (1) and (2)

\begin{matrix} \frac{d \vec{PP ″}}{dt} = & \frac{d \vec{OP ″}}{dt} - \frac{d \vec{OP}}{dt} = - z \frac{\partial {\overset{\cdot}{w}}_{rel} (x, t)}{\partial x} {\vec{e}}_{1} \\ + ({\overset{\cdot}{w}}_{base} (t) + {\overset{\cdot}{w}}_{rel} (x, t)) {\vec{e}}_{3} \end{matrix}

(3)

And the relative displacement of $P ″$ with respect to $P'$ is derived as

\vec{P' P ″} = \vec{PP ″} - \vec{PP'} = - z \frac{\partial w_{rel} (x, t)}{\partial x} {\vec{e}}_{1} + w_{rel} (x, t) {\vec{e}}_{3}

(4)

Based on the continuum thermodynamics, the strain tensor can be denoted as

ε_{ij} = \frac{1}{2} (u_{i, j} + u_{j, i})

(5)

where $ε_{ij}$ is the strain tensor, and $u_{i, j}$ denotes the displacement derivative with respect to coordinate j. By differentiating equation (4) with respect to x, the stretch strain field across the whole composite beam structure can be expressed as

ε_{xx} = ε_{1} = \frac{\partial (\vec{P' P ″} \cdot {\vec{e}}_{1})}{\partial x} = - z \frac{\partial^{2} w_{rel} (x, t)}{\partial x^{2}} = - z ε_{1}^{(1)}

(6)

As the linear electrical enthalpy of piezoelectric material can be approximated by^7,8

H (ε_{kl}, E_{k}) = \frac{1}{2} c_{ijkl}^{E} ε_{ij} ε_{kl} - e_{ijk} E_{i} ε_{jk} - \frac{1}{2} ζ_{ij}^{ε} E_{i} E_{j}

(7)

Then, the constitutive relation derived from equation (7) can be obtained as

σ_{ij} = \frac{\partial H}{\partial ε_{ij}} = c_{ijkl}^{E} ε_{kl} - e_{kij} E_{k}

(8)

D_{i} = - \frac{\partial H}{\partial E_{i}} = e_{ikl} ε_{kl} + ζ_{ik}^{ε} E_{k}

(9)

In equations (7)–(9), the terms $c_{ijkl}^{E}$ , $e_{ijk}$ , and $ζ_{ij}^{E}$ are the elastic stiffness at constant electric field, piezoelectric coupling coefficient, and permittivity at constant strain, respectively. $σ_{ij}$ , $D_{i}$ , and $E_{k}$ represent the stress tensor, electric displacement vector, and electric field vector that are expressed in Voigt’s notation. The total energy equation of the piezoelectric composite beam is expressed in the Hamiltonian form shown as^7,9

\int_{t_{1}}^{t_{2}} (δ KE - δ PE + δ WE + δ W_{f}) dt = 0

(10)

where KE, PE, WE, and $W_{f}$ are the kinetic energy, potential energy, electric energy derived by the direct piezoelectric effect from piezoelectric patches, and mechanical energy from the external host structure vibration, respectively. If the piezoelectric layer is polarized upward in thickness direction and electrodes are set at major surfaces of the beam, then electric field will be induced in the vertical direction. Considering the assumption of Euler–Bernoulli’s beam theory that only the strain component $ε_{xx}$ is nontrivial while others are all zero, equations (8) and (9) will be reduced as

σ_{xx} = c_{11}^{E} ε_{1} - e_{31} E_{3}

(11)

D_{3} = e_{31} ε_{1} + ζ_{33}^{ε} E_{3}

(12)

where $E_{3}$ is the electric field in vertical direction. Combining equations (3)–(12), the Lagrangian $L_{a}$ of the total energy that takes both external mechanical and electrical works into consideration can be deduced as

\begin{matrix} \int_{t_{1}}^{t_{2}} δ (L_{a} + W_{f}) dt = \int_{t_{1}}^{t_{2}} (δ KE - δ PE + δ WE + δ W_{f}) dt \\ = \int_{t_{1}}^{t_{2}} {\int_{S} [I^{(C, k)} \frac{\partial^{2} {\overset{\cdot\cdot}{w}}_{rel}}{\partial x^{2}} δ w_{rel}] ds} dt + \int_{t_{1}}^{t_{2}} {\int_{S} [- I^{(A, k)} {\overset{\cdot\cdot}{w}}_{rel} δ w_{rel} - I^{(A, k)} {\overset{\cdot\cdot}{w}}_{base} δ w_{rel}] dS} dt \\ - \int_{t_{1}}^{t_{2}} {\int_{S} [C_{11}^{(F, k)} \frac{\partial^{4} w_{rel}}{\partial x^{4}} δ w_{rel} - R_{31}^{(H, k)} \frac{\partial^{2} v (t)}{\partial x^{2}} δ w_{rel}] dS} dt + \int_{t_{1}}^{t_{2}} {\int_{S} [R_{31}^{(H, k)} \frac{\partial^{2} w_{rel}}{\partial x^{2}} + S_{33}^{(k)} v (t)] δ v (t) dS} dt + \int_{t_{1}}^{t_{2}} q δ v (t) dt \\ - \int_{t_{1}}^{t_{2}} {\oint_{S} n_{x} (C_{11}^{(F, k)} \frac{\partial^{2} w_{rel}}{\partial x^{2}}) δ (\frac{\partial w_{rel}}{\partial x}) dS} dt + \int_{t_{1}}^{t_{2}} {\oint_{S} n_{x} C_{11}^{(F, k)} \frac{\partial^{3} w_{rel}}{\partial x^{3}} δ w_{rel} dS} dt + \int_{t_{1}}^{t_{2}} {\oint_{S} n_{x} R_{31}^{(H, k)} v (t) δ (\frac{\partial w_{rel}}{\partial x}) dS} dt \\ - \int_{t_{1}}^{t_{2}} {\oint_{S} n_{x} R_{31}^{(H, k)} \frac{\partial v (t)}{\partial x} δ w_{rel} dS} dt - \int_{t_{1}}^{t_{2}} {\oint_{S} I^{(C, k)} n_{x} \frac{\partial {\overset{\cdot\cdot}{w}}_{rel}}{\partial x} δ w_{rel} ds} dt \end{matrix}

(13)

where $I^{(A, k)}$ and $I^{(C, k)}$ are the zeroth and second mass moments of inertia per unit area with respect to the characteristic material property and cross-sectional area of piezoelectric layer and substrate layer, respectively. q and $v (t)$ are the electric charges and electric potential, respectively. The coefficients $C_{11}^{(F, k)}$ , $R_{31}^{(H, k)}$ , and $S_{33}^{(k)}$ are defined as follows

\begin{matrix} C_{11}^{(F, k)} & = \int_{- (t_{s} + t_{p} - z_{s})}^{z_{s} - t_{p}} z^{2} E_{s} dz + \int_{z_{s} - t_{p}}^{z_{s}} z^{2} c_{11}^{E} dz \\ = \frac{1}{3} E_{s} [t_{s}^{3} - 3 t_{s} (z_{s} - t_{p}) (t_{s} - z_{s} + t_{p})] \\ + \frac{1}{3} c_{11}^{E} (t_{p}^{3} + 3 z_{s}^{2} t_{p} - 3 z_{s} t_{p}^{2}) \end{matrix}

(14)

R_{31}^{(H, k)} = \int_{z_{s} - t_{p}}^{z_{s}} z e_{31} (\frac{1}{t_{p}}) dz = e_{31} (z_{s} - \frac{t_{p}}{2})

(15)

S_{33}^{(k)} = \int_{z_{s} - t_{p}}^{z_{s}} \frac{ζ_{33}^{ε}}{t_{p}^{2}} dz = \frac{ζ_{33}^{ε}}{t_{p}}

(16)

In the equations, $t_{p}$ , $t_{s}$ , $z_{s}$ , and b represent the piezoelectric layer thickness, substrate layer thickness, distance between neutral plane and top surface, and composite beam width, respectively. $c_{11}^{E}$ and $E_{s}$ denote the stiffness coefficient for piezoelectric layer and Young’s modulus for substrate, respectively. $ζ_{33}$ is the relative permittivity. From equation (13), it can be seen that the Lagrangian theorem can be reformulated in terms of virtual relative displacement and electrical potential for the basic constitutive electromechanical dynamic equations of piezoelectric unimorph cantilever beam. With the parameter of virtual relative transverse displacement and electrical potential being set as zero, the corresponding electromechanical dynamic equations can be deduced as follows

\begin{matrix} δ w_{rel} : {\hat{I}}^{(C, k)} \frac{\partial^{2} {\overset{\cdot\cdot}{w}}_{rel}}{\partial x^{2}} - {\hat{I}}^{(A, k)} {\overset{\cdot\cdot}{w}}_{rel} - {\hat{I}}^{(A, k)} {\overset{\cdot\cdot}{w}}_{base} \\ - {\hat{C}}_{11}^{(F, k)} \frac{\partial^{4} w_{rel}}{\partial x^{4}} + {\hat{R}}_{31}^{(H, k)} \frac{\partial^{2} v (t)}{\partial x^{2}} = 0 \end{matrix}

(17)

δ v : \int_{0}^{L} [{\hat{R}}_{31}^{(H, k)} \frac{\partial^{2} w_{rel}}{\partial x^{2}} + {\hat{S}}_{33}^{(k)} v (t)] dx + q = 0

(18)

where ${\hat{I}}^{(A, k)}$ , ${\hat{I}}^{(C, k)}$ , ${\hat{C}}_{11}^{(F, k)}$ , ${\hat{R}}_{31}^{(H, k)}$ , and ${\hat{S}}_{33}^{(k)}$ denote the integration of the corresponding terms with respect to y. In order to characterize the dynamic equation, it should first derive the natural frequency and its corresponding mode shapes, where external excitations including both terms ${\overset{\cdot\cdot}{w}}_{base}$ and $v (t)$ will not be considered. Therefore, it is assumed that $w_{rel} = Φ (x) e^{i ω t}$ . And the corresponding eigenfunction can be derived as

\frac{\partial^{4} Φ (x)}{\partial x^{4}} + \frac{{\hat{I}}^{(C, k)} ω^{2}}{{\hat{C}}_{11}^{(F, k)}} \frac{\partial^{2} Φ (x)}{\partial x^{2}} - \frac{{\hat{I}}^{(A, k)} ω^{2}}{{\hat{C}}_{11}^{(F, k)}} Φ (x) = 0

(19)

The solution to the ordinary differential equation (19) is obtained as

Φ (x) = C_{1} \cosh μ x + C_{2} \sinh μ x + C_{3} \cos ν x + C_{4} \sin ν x

(20)

where

μ = \sqrt{- \frac{{\hat{I}}^{(C, k)} ω^{2}}{2 {\hat{C}}_{11}^{(F, k)}} + \frac{1}{2} \sqrt{{(\frac{{\hat{I}}^{(C, k)} ω^{2}}{{\hat{C}}_{11}^{(F, k)}})}^{2} + 4 \frac{{\hat{I}}^{(A, k)} ω^{2}}{{\hat{C}}_{11}^{(F, k)}}}}

ν = \sqrt{\frac{{\hat{I}}^{(C, k)} ω^{2}}{2 {\hat{C}}_{11}^{(F, k)}} + \frac{1}{2} \sqrt{{(\frac{{\hat{I}}^{(C, k)} ω^{2}}{{\hat{C}}_{11}^{(F, k)}})}^{2} + 4 \frac{{\hat{I}}^{(A, k)} ω^{2}}{{\hat{C}}_{11}^{(F, k)}}}}

The associated boundary conditions are given as

\begin{matrix} Φ (0) = & 0, \frac{\partial Φ (0)}{\partial x} = 0, \frac{\partial^{2} Φ (L)}{\partial x^{2}} = 0 and \frac{\partial^{3} Φ (L)}{\partial x^{3}} \\ + \frac{{\hat{I}}^{(C, k)} ω^{2}}{{\hat{C}}_{11}^{(F, k)}} \frac{\partial Φ (L)}{\partial x} = 0 \end{matrix}

(2)

Combining equations (20) and (21), the coefficients $C_{1}$ to $C_{4}$ can be determined according to the equation below

{\begin{matrix} Φ (0) = C_{1} + C_{3} = 0 \\ \frac{\partial Φ (0)}{\partial x} = μ C_{2} + ν C_{4} = 0 \\ \frac{\partial^{2} Φ (L)}{\partial x^{2}} = μ^{2} \cosh μ L C_{1} + μ^{2} \sinh μ L C_{2} - ν^{2} \cos ν L C_{3} - ν^{2} \sin ν L C_{4} = 0 \\ \frac{\partial^{3} Φ (L)}{\partial x^{3}} + \frac{{\hat{I}}^{(C, k)} ω^{2}}{{\hat{C}}_{11}^{(F, k)}} \frac{\partial Φ (L)}{\partial x} = (μ^{3} \sinh μ L + \frac{{\hat{I}}^{(C, k)} ω^{2}}{{\hat{C}}_{11}^{(F, k)}} μ \sinh μ L) C_{1} + (μ^{3} \cosh μ L + \frac{{\hat{I}}^{(C, k)} ω^{2}}{{\hat{C}}_{11}^{(F, k)}} μ \cosh μ L) C_{2} \\ + (ν^{3} \sin ν L - \frac{{\hat{I}}^{(C, k)} ω^{2}}{{\hat{C}}_{11}^{(F, k)}} ν \sin ν L) C_{3} - (ν^{3} \cos ν L - \frac{{\hat{I}}^{(C, k)} ω^{2}}{{\hat{C}}_{11}^{(F, k)}} ν \cos ν L) C_{4} = 0 \end{matrix}

(22)

In order to obtain a nontrivial value of $C_{1}$ to $C_{4}$ , the determinant of coefficient matrix should be equal to zero, and thus, the values for natural frequencies $ω_{i}$ can be determined in numerical methods. Based on the obtained result, it is easy to derive the mode shape function as equation (23)

Φ (x) = C_{1} (\cosh μ x - a \frac{ν}{μ} \sinh μ x - \cos ν x + a \sin ν x)

(23)

where

a = \frac{μ^{2} \cosh μ L + ν^{2} \cos ν L}{μ ν \sinh μ L + ν^{2} \sin ν L}

Multimode frequency response of the unimorph cantilever beam

If Ritz method is utilized to represent the transverse displacement, the solution form can be expressed as

w_{rel} (x, t) = \sum_{r = 1}^{m} c_{r} Φ_{r} (x) e^{i ω t}

(24)

where $c_{r}$ is the unknown Ritz coefficients that can be determined by eigenfunction

\sum_{r = 1}^{m} [K_{qr} - ω^{2} M_{qr}] c_{r} = 0, q = 1, 2, \dots, m

(25)

As the values of $ω$ can be determined by equation (22), each value of $ω$ can determine a corresponding eigenvector of ${c_{r}}$ based on equation (25). Therefore, the arbitrary Ritz mode shapes can conversely be expressed in terms of both eigenvector $c_{kr}$ and truncated mode shapes $Φ_{k} (x)$

Φ_{r} (x) = \sum_{k = 1}^{m} c_{kr} Φ_{k} (x) r = 1, 2, \dots, m

(26)

If $Φ_{r} (x)$ is normalized with respect to the mass

\begin{matrix} {\hat{Φ}}_{r} (x) = \frac{Φ_{r} (x)}{{(\int_{0}^{L} {\hat{I}}^{(A, k)} Φ_{r}^{2} (x) dx + \int_{0}^{L} {\hat{I}}^{(C, k)} {(\frac{\partial Φ_{r} (x)}{\partial x})}^{2} dx)}^{1 / 2}} \\ r = 1, 2, \dots, m \end{matrix}

(27)

Then, equation (24) can be rearranged as

w_{rel} (x, t) = \sum_{r = 1}^{m} {\hat{Φ}}_{r} (x) w_{r} (t)

(28)

By substituting equation (28) into equations (17) and (18), the dynamic equations can be rewritten as

{\overset{\cdot\cdot}{w}}_{r} + 2 ς_{r} ω_{r} {\overset{\cdot}{w}}_{r} + ω_{r}^{2} w_{r} - {\hat{P}}_{r} v = - {\hat{Q}}_{r} {\overset{\cdot\cdot}{w}}_{base}

(29)

\sum_{r = 1}^{m} {\hat{P}}_{r} {\overset{\cdot}{w}}_{r} + P_{D} \overset{\cdot}{ν} + R_{L} ν = 0

(30)

where

\begin{matrix} 2 ς_{r} ω_{r} = α + β ω_{r}^{2} \\ {\hat{P}}_{r} = \int_{0}^{L} {\hat{R}}_{31}^{(H, k)} \frac{\partial^{2} {\hat{Φ}}_{r} (x)}{\partial x^{2}} dx \\ {\hat{Q}}_{r} = \int_{0}^{L} {\hat{I}}^{(A, k)} {\hat{Φ}}_{r} (x) dx \end{matrix}

If Laplace transforms are applied into equations (29) and (30), the transfer functions of $w_{abs} (x, t)$ and ${\overset{\cdot}{w}}_{abs} (x, t)$ related to input acceleration of ${\overset{\cdot\cdot}{w}}_{base} (t)$ will be deduced as equations (31) and (32)

\begin{matrix} H_{1}^{abs} (x, j ω) = \sum_{r = 1}^{m} \\ [\frac{{\hat{Φ}}_{r} (x)}{ω_{r}^{2} - ω^{2} + j 2 ζ_{r} ω_{r} ω} (\frac{{\hat{P}}_{r} \sum_{r = 1}^{m} \frac{j {\hat{P}}_{r} {\hat{Q}}_{r} ω}{ω_{r}^{2} - ω^{2} + j 2 ζ_{r} ω_{r} ω}}{j ω P_{D} + R_{L} + \sum_{r = 1}^{m} \frac{j ω P_{r}^{2}}{ω_{r}^{2} - ω^{2} + j 2 ζ_{r} ω_{r} ω}} - {\hat{Q}}_{r})] - \frac{1}{ω^{2}} \end{matrix}

(31)

\begin{matrix} H_{4}^{abs} (x, j ω) = j ω \sum_{r = 1}^{m} \\ [\frac{{\hat{Φ}}_{r} (x)}{ω_{r}^{2} - ω^{2} + j 2 ζ_{r} ω_{r} ω} (\frac{{\hat{P}}_{r} \sum_{r = 1}^{m} \frac{j {\hat{P}}_{r} {\hat{Q}}_{r} ω}{ω_{r}^{2} - ω^{2} + j 2 ζ_{r} ω_{r} ω}}{j ω P_{D} + R_{L} + \sum_{r = 1}^{m} \frac{j ω P_{r}^{2}}{ω_{r}^{2} - ω^{2} + j 2 ζ_{r} ω_{r} ω}} - {\hat{Q}}_{r})] + \frac{1}{j ω} \end{matrix}

(32)

The time-dependent dynamical responses of voltage and power outputs are deduced as equations (33) and (34), respectively

v (t) = - ω^{2} w_{base} e^{j ω t} \frac{\sum_{r = 1}^{m} \frac{j ω {\hat{P}}_{r} {\hat{Q}}_{r}}{ω_{r}^{2} - ω^{2} + j 2 ζ_{r} ω_{r} ω}}{j ω P_{D} + R_{L} + \sum_{r = 1}^{m} \frac{j ω {\hat{P}}_{r}^{2}}{ω_{r}^{2} - ω^{2} + j 2 ζ_{r} ω_{r} ω}}

(33)

P (t) = \frac{ω^{4} w_{base}^{2} e^{j 2 ω t}}{R_{load}} {[\frac{\sum_{r = 1}^{m} \frac{j ω {\hat{P}}_{r} {\hat{Q}}_{r}}{ω_{r}^{2} - ω^{2} + j 2 ζ_{r} ω_{r} ω}}{j ω P_{D} + R_{L} + \sum_{r = 1}^{m} \frac{j ω {\hat{P}}_{r}^{2}}{ω_{r}^{2} - ω^{2} + j 2 ζ_{r} ω_{r} ω}}]}^{2}

(34)

Correspondingly, the frequency response functions (FRFs) of $v (t)$ and $P (t)$ related to input acceleration of ${\overset{\cdot\cdot}{w}}_{base} (t)$ will be deduced as

H_{2} (j ω) = \frac{\sum_{r = 1}^{m} \frac{j ω {\hat{P}}_{r} {\hat{Q}}_{r}}{ω_{r}^{2} - ω^{2} + j 2 ζ_{r} ω_{r} ω}}{j ω P_{D} + R_{L} + \sum_{r = 1}^{m} \frac{j ω {\hat{P}}_{r}^{2}}{ω_{r}^{2} - ω^{2} + j 2 ζ_{r} ω_{r} ω}}

(35)

H_{3} (ω) = \frac{1}{R_{load}} {[\frac{\sum_{r = 1}^{m} \frac{j ω {\hat{P}}_{r} {\hat{Q}}_{r}}{ω_{r}^{2} - ω^{2} + j 2 ζ_{r} ω_{r} ω}}{j ω P_{D} + R_{L} + \sum_{r = 1}^{m} \frac{j ω {\hat{P}}_{r}^{2}}{ω_{r}^{2} - ω^{2} + j 2 ζ_{r} ω_{r} ω}}]}^{2}

(36)

From equation (34), the optimal harvested power will be determined by setting the external load resistance $R_{load}$ as

\begin{matrix} R_{load} = \\ \frac{| \sum_{r = 1}^{m} \frac{ω {\hat{P}}_{r}^{2} (2 ζ_{r} ω_{r} ω)}{{(ω_{r}^{2} - ω^{2})}^{2} + {(2 ζ_{r} ω_{r} ω)}^{2}} - j [ω P_{D} + \sum_{r = 1}^{m} \frac{(ω_{r}^{2} - ω^{2}) ω {\hat{P}}_{r}^{2}}{{(ω_{r}^{2} - ω^{2})}^{2} + {(2 ζ_{r} ω_{r} ω)}^{2}}] |}{{[\sum_{r = 1}^{m} \frac{ω {\hat{P}}_{r}^{2} (2 ζ_{r} ω_{r} ω)}{{(ω_{r}^{2} - ω^{2})}^{2} + {(2 ζ_{r} ω_{r} ω)}^{2}}]}^{2} + {[ω P_{D} + \sum_{r = 1}^{m} \frac{(ω_{r}^{2} - ω^{2}) ω {\hat{P}}_{r}^{2}}{{(ω_{r}^{2} - ω^{2})}^{2} + {(2 ζ_{r} ω_{r} ω)}^{2}}]}^{2}} \end{matrix}

(37)

The optimal resistance is obtained by differentiating equation (36) with respect to $R_{load}$ and setting the differentiation function to be zero.

Case study of piezoelectric unimorph cantilever scavenger in bending mode

Analytical simulation

In order to prove the feasibility of the model, the basic unimorph cantilever structure for the vibration energy harvesting device is proposed. Piezoelectric material of PZT-5H^10,11 is used in the top layer where the polarized direction is set upward in the vertical direction. Parameters for PZT-5H including stiffness constant $c_{11}$ , mass density $ρ_{p}$ , piezoelectric coupling coefficient $d_{31}$ , and relative permittivity $ζ_{33}$ are $126 GPa$ , $7500 kg / m^{3}$ , $- 274.8 pm / V$ , $3200 ζ_{0}$ , and $8.854 \times 10^{- 12} F / m$ , respectively. The material for elastic substrate is selected to be 45 steel whose Young’s modulus, Poisson’s ratio, and mass density are $205 GPa$ , 0.28, and $7850 kg / m^{3}$ , respectively. The dimension of the geometric structure and its corresponding inertial moment of mass are listed in Table 1.

Table 1.

Dimensions of geometric structure and corresponding coefficients.

Geometric properties	Piezoelectric layer	Elastic substrate
Length, L (mm)	40	40
Thickness, t (mm)	0.11	0.24
Width, b (mm)	8	8
$z_{s} (mm)$	0.1915
${\hat{I}}^{(A, k)}$	$2.1672 \times 10^{- 2}$
${\hat{I}}^{(C, k)}$	$4.0769 \times 10^{- 10}$
${\hat{C}}_{11}^{(F, k)}$	$4.6504 \times 10^{- 3}$
${\hat{R}}_{31}^{(H, k)}$	$- 7.098 \times 10^{- 6}$
${\hat{S}}_{33}^{(k)}$	$2.06057 \times 10^{- 6}$

Since the order of ${\hat{I}}^{(C, k)}$ is small compared with other parameters, the value of ${\hat{I}}^{(C, k)}$ can be treated as 0. Therefore, the determinant of coefficient matrix of equation can be reduced as

1 + \cosh μ L \cos μ L = 0

(38)

Through the numerical calculations, the first three resonant frequencies are determined to be 162.0125, 1015.3149, and 2842.9110 Hz, respectively. The first three bending mode shapes can be obtained and they are plotted as shown in Figure 2.

Figure 2.

Bending mode shapes of the cantilever beam at 162.0125, 1015.3149, and 2842.9110 Hz.

By further substitution of these parameters and normalized mode shapes into equation (35), the FRF of output voltage is determined. From the results shown in Figure 3(a), the highest output voltages around the first three modes are 0.526 V/g, 0.00437 V/g, and 0.1944 mV/g in open circuit. With the load resistance being decreased, the output voltage monotonically diminishes and finally reaches a minimum value of 0.00471 V/g at 100 Ω. However, when it comes to power output response, the maximum power output at 162 Hz is 13.6 µW/g², where the impedance of the piezoelectric energy harvester is matched with a resistance of 11.043 kΩ. Similarly, the output power generated across the load resistances at 100 Ω, 679 Ω, 1.9 kΩ, and 30 kΩ all present peak values at 160 Hz, while for peak value of output power at 10 MΩ, the corresponding frequency is shifted to 162.5 Hz. This frequency shift is due to the fact that the response functions of voltage and power are related to load resistance and piezoelectric power; the frequency of the peak value of output power is shifted according to formulas (35) and (36). For second mode at 1015 Hz, the maximum power output is 5.13 nW/g², where the optimum load resistance of 1.9 kΩ is proved. The resonant frequencies for devices with load resistances of 100 Ω, 679 Ω, and 1.9 kΩ are 1015 Hz, while for that of 11.043 kΩ, 30 kΩ, and 10 MΩ, the frequencies are shifted to be 1014.5 Hz. In the third mode, the maximum power output is 27.614 pW/g² at the load resistance of 679 Ω. The resonant frequencies for the load resistance of 100 Ω, 679 Ω, and 1.9 kΩ are 2824, 2817.5, and 2811.5 Hz, respectively. The third resonant frequencies in other cases are all 2810 Hz. The effective factors to the resonant frequency include the geometric dimensions, Young’s modulus, stiffness, and mass densities, which can be found from the combination of formula (38), expression for µ, and $I^{(A, k)}$ and $C_{11}^{(F, k)}$ . As the resonant frequency is determined, the effective factors to the response function of voltage and power are $P_{D}$ , $R_{load}$ , $ζ_{r}$ ${\hat{P}}_{r}$ , and ${\hat{Q}}_{r}$ , which are related with mass densities of substrate layer and piezoelectric layer, piezoelectric coupling coefficient, geometric dimensions, and shape functions. As equation (20) expresses, the coefficients of the shape function depends on the boundary conditions.

Figure 3.

Frequency response functions (FRFs) of (a) output voltage and (b) power.

Numerical analysis

In this section, we use the finite element analysis package ANSYS 12.0 to validate the proposed analytical method. Although the geometric structure dimensions are still the same as described in section “Multimode frequency response of the unimorph cantilever beam,” there are some differences on material constants that need to be clarified. Finite element method (FEM) has been widely used in modeling mechanical energy harvesting devices.¹² In the analytical method, as the beam model is simplified to Euler–Bernoulli type and only transverse displacement is considered, beam width plays no effect on the whole device output performance. Therefore, the geometric model can be built in a two-dimensional domain in ANSYS software. Here, plane elements plane223 and plane82 are used for modeling the piezoelectric and substrate layers. Since plane223 is used in a piezoelectric-coupled field analysis, it works either in a plane stress or plane strain condition, where more material constants such as $c_{22}$ , $c_{33}$ , and $c_{44}$ should be considered for the positive definite matrix of $[c^{E}]$ . Piezoelectric material constants of $ζ_{11}$ and $e_{33}$ should also be considered in the FEM for proper calculations. It is noted that since the coupling effect among the strain components is not considered in the analytical assumption, the terms $c_{12}$ and $c_{13}$ are treated to be 0 in $[c^{E}]$ . The geometric dimensions are exactly the same as shown in Table 1. Additional material parameters that are defined in the FEM are shown in Table 2.

Table 2.

Material parameters for PZT-5H.

Elastic stiffness (GPa)	Coupling coefficients (c/m²)	Dielectric constants (c/m²)
$c_{22} = 126$	$e_{31} = - 6.5$	$ζ_{31} = 3200$
$c_{33} = 117$	$e_{32} = - 6.5$	$ζ_{32} = 3200$
$c_{44} = 23.3$	$e_{33} = 23.3$

In the modal analysis, the composite beams are modeled by two glued rectangles where one end is fixed and the other end is kept free. On the top and bottom surfaces of piezoelectric layer, the electrodes are modeled using the coupled commands. The mode shapes of the geometric structure in displacement contour are shown in Figure 4. The mode shape of the rectangular cantilever beam in FEM is the same as that in analytical simulation, which can be proved by the calculated resonant frequencies and the shapes. Although the shapes at mode 2 between analytical and FEM methods are in opposite directions, they are attributed to the negative value of coefficient $c_{kr}$ been chosen in equation (23). The calculated resonant frequencies of 162.26, 1015.53, and 2841.02 Hz coincide well with analytical results of 162.0125, 1015.3149, and 2842.9110 Hz, respectively.

Figure 4.

Mode shapes of the rectangular cantilever beam in FEM.

As the first three mode shapes and resonant frequencies are determined, the dynamic response of output voltage and power should be characterized as well. In analytical method, the optimum resistance around the first three resonant modes is determined to be 11.043 kΩ, 1.9 kΩ, and 679 Ω according to equation (37). In FEM, we calculate the optimum external load resistance according to

R_{load} = \frac{1}{2 π f_{r} C}

(39)

where $f_{r}$ and C denote the resonant frequency and piezoelectric capacitance, respectively. The capacitance of piezoelectric layer is determined by the extracted charges in the static analysis by setting a voltage of 1 V across the electrodes. Here, as the geometric structure is in two-dimensional domain, the obtained value of capacitance is $1.218791812 \times 10^{- 5} F / m$ . Correspondingly, the optimum load resistances at resonant frequencies of 162.26, 1015.53, and 2841.02 Hz are solved to be 80.481, 12.859, and 4.597 Ωm, respectively. In actual three-dimensional models, the beam width should be considered. The optimum resistances are calculated to be 10.06 kΩ, 1.607 kΩ, and 574.565 Ω by dividing 80.481, 12.859, and 4.597 Ωm with the beam width of 8 mm. In addition, the resistances of 100 Ω, 30 kΩ, and 10 MΩ are also included for a direct comparison with analytical output results. In Figure 5(a), the output voltage is characterized by sweeping the driven frequency from 0 to 3000 Hz with a step of 1 Hz. From the results, it is found that the output voltage is decreased from mode 1 to mode 3, but increased when the external load resistance varies from 100 Ω to 10 MΩ. The peak output voltages across 10 MΩ resistor around the three resonant modes are 0.93739, 0.04958, and 0.00606 V/g at 162, 1015, and 2829 Hz, respectively, which are all a bit overestimated compared with analytical results. The reason to this phenomenon is that the coupling effect of $e_{32}$ and $e_{33}$ on the output voltage is considered. However, the predicted varying trends of the output voltage responses are still similar with each other between analytical and finite element models.

Figure 5.

Frequency response functions (FRFs) of (a) output voltage and (b) power in FEM.

The maximum power output at 162 Hz is 35.4697 µW/g², where the optimum resistance of 10.06 kΩ is matched to the piezoelectric energy harvester. The generated peak output power at 162 Hz across 100 Ω, 574.565 Ω, 1.607 kΩ, 30 kΩ, and 10 MΩ are 0.98979, 5.49292, 14.0783, 22.6027, and 0.008787 µW/g², respectively. For the second mode at 1015 Hz, the maximum power output is 743.688 nW/g², where the optimum resistance is changed to be 1.607 kΩ. Correspondingly, the induced power across 100 Ω, 574.565 Ω, 10.06 kΩ, 30 kΩ, and 10 MΩ are 93.7168, 474.044, 236.42, 81.4825, and 0.24578 nW/g², respectively. In the third mode, the peak output power is 31.0336 nW/g² at 2836 Hz, where the matched external load resistance is 574.565 Ω. Similar to previous analytical results, the resonant frequencies across 100 Ω, 1.607 kΩ, 10.06 kΩ, 30 kΩ, and 10 MΩ are shifted to be 2843, 2831, 2830, 2829, and 2829 Hz, respectively. Through the comparison of the results between analytical method and FEM, the power output around the first resonant mode is comparable with each other, while the difference in the outputs at second and third modes is quite large. The reason to this result is attributed to the assumption of plane strain in FEM. In FEM, not only $c_{11}$ , $ζ_{11}$ , and $e_{33}$ are considered but also the coefficients $c_{22}$ , $c_{33}$ , $c_{44}$ , $ζ_{11}$ , and $e_{33}$ are taken into account in calculating the values of $K_{qr}$ , $M_{qr}$ , $Q_{q}$ , $P_{q}$ , $P_{r}$ , and $P_{D}$ . Regarding the simulation in ANSYS software, the two-dimensional simulation of the bending mode cantilever energy harvester is based on the plane strain assumption, where the constitutive equation is expressed as the following formula

\begin{matrix} σ_{1} = & (c_{11} - \frac{c_{13}^{2}}{c_{33}}) ε_{1} + (c_{12} - \frac{c_{13} c_{32}}{c_{33}}) ε_{2} \\ + (e_{13} - \frac{c_{13}}{c_{33}} e_{33}) E_{3} \end{matrix}

(40)

\begin{matrix} σ_{2} = & (c_{21} - \frac{c_{23} c_{31}}{c_{33}}) ε_{1} + (c_{22} - \frac{c_{23} c_{32}}{c_{33}}) ε_{2} \\ + (e_{23} - \frac{c_{23}}{c_{33}} e_{33}) E_{3} \end{matrix}

(41)

σ_{44} = c_{44} ε_{4}

(42)

Comparing the modified constitutive equations with formulas (11) and (12), the terms $c_{12}$ , $c_{13}$ , $c_{23}$ , and $c_{32}$ and $e_{13}$ , $e_{23}$ , and $e_{33}$ play important roles in the dynamic responses. Although the terms $c_{12}$ , $c_{13}$ , $c_{23}$ , and $c_{32}$ are treated to be 0, the coefficients $c_{22}$ , $c_{33}$ , $c_{44}$ , $ζ_{11}$ , and $e_{33}$ still make extra impact on the output voltage and power responses at corresponding resonant modes compared with the analytical method. However, $c_{11}$ , $ζ_{11}$ , and $e_{33}$ are dominant coefficients in this calculation in the first resonant mode. While for higher modes, the effects of $c_{22}$ , $c_{33}$ , and $c_{44}$ should not be neglected. When the driven frequency is increased up to higher resonant frequencies, the deviation of the output voltage and power with respect to the analytical results can be augmented according to equations (35) and (36).

Transient analysis of the dynamic response

In the transient analysis using FEM, we focus on the dynamic response at the first mode, where the external driving frequency is set constant at 162.26 Hz. The vibration displacement amplitude at the fixed end is set with a value of 3.96 µm, which means that the acceleration amplitude is defined to be 0.104 m/s² (approximately 0.01 g). External load resistor of 10.06 kΩ is connected to the piezoelectric layer. When the device is working, the fixed end of the cantilever vibrates rigidly with the external excitation. But for the free end, it will induce relative transverse displacement with respect to the rigid moving base because of the transverse bending motion. Through the finite element analysis, the displacement response of the tip is shown in Figure 6(a) and (b). From the sinusoidal wave, it is found that the displacement oscillates cyclically around the equilibrium point with a time period of 6.16 ms. Initially, the tip motion is in a transient status in 0.72152 s. After it transcends this time period, the oscillation is kept stable at amplitude of 0.543 mm, which can be observed in the zoomed out Figure 6(b).

Figure 6.

(a) Transient dynamic response of tip displacement in 2 s and (b) zoomed out figure of transient dynamic response of the tip displacement at stable state.

The transient dynamic response of the output voltage across the load resistor is shown in Figure 7(a) and (b), where the output voltage response enters a stable state from 0.72383 s. The amplitude of the output voltage is 0.9669 V with a frequency of 162.26 Hz. However, when it comes to the generated instantaneous power across 10.06 kΩ, the frequency is double of 162.26 Hz, as the instantaneous power is calculated by $P = (1 / 2) (V^{2} / R)$ , where both the positive and negative voltages could be utilized to power the external electric load. This can be observed in Figure 8(b). The dynamic response of the generated power is shown in Figure 8(a), where the maximum instantaneous output power of 92.828 µW is first reached at 0.68069 s. In the time period of 2 s, the induced energy from the bending mode motion can be calculated. Figure 9 shows the relationship between the induced electric energy and the external excitation. When the base is excited with an acceleration of 0.001 g, the dissipated energy across the external load resistance is 0.817 µJ. With the excitation increased from 0.001 to 0.27 g, the maximum converted energy of 59.6 mJ can be obtained. As for the amplitude of output voltage response, it is almost in a linear relationship with the output excitation.

Figure 7.

(a) Transient dynamic response of output voltage in 2 s and (b) zoomed out figure of transient dynamic response of the output voltage at stable state.

Figure 8.

(a) Transient dynamic response of output power in 2 s and (b) zoomed out figure of transient dynamic response of the output power at stable state.

Figure 9.

Open voltage and dissipated electrical energy versus base excitation.

Conclusion

In this article, the analytical model of piezoelectric unimorph cantilever energy harvester has been provided based on Euler–Bernoulli’s beam theory. By expressing the Hamiltonian form of the total energy equation into strong form and weak form, the mode shape functions and electromechanical dynamic equations have been obtained. Utilizing the Ritz method and Laplace transforms, the FRFs of output voltage and generated power across load resistor are provided. To prove the validity of analytical model, a specific unimorph piezoelectric energy harvester has been studied. By substituting the related parameters in the derived formulations, bending mode shapes as well as first three resonant frequencies of 162, 1015, and 2842 Hz are obtained. In order to characterize the dynamic behavior of the device, the varying trends of output voltage and generated power versus driven frequencies are predicted. In this way, the optimum load resistances at corresponding resonant modes are determined. For confirming the effectiveness of analytical model, the proposed piezoelectric energy harvester is also simulated in commercial finite element software ANSYS. Through the specified definition of the material property and modeling elements, the calculated first three bending mode resonant frequencies—162, 1015, and 2841 Hz—coincide well with analytical results. Through the static analysis of piezoelectric energy harvester, the resistance of the optimum power output corresponding to each mode is found to be comparable with analytical results as well. By conducting the harmonic analysis, the varying trends of FRFs of output voltage and generated power also agree well with each other. Only one issue that should be noted is that the output values at second and third modes are overestimated compared with the analytical model. This phenomenon is attributed to the fact that the additional coefficients $c_{22}$ , $c_{33}$ , $c_{44}$ , $ζ_{11}$ , and $e_{33}$ have been taken into account. When the swept frequency is increased from the first mode to higher modes, the minor difference between the two methods around the first mode is amplified by the increased driving frequency. However, the results at first mode matched well with each other. Finally, the transient dynamic behavior of the cantilever energy harvester has been examined in FEM. As the external excitation is increased from 0.001 to 0.27 g, the output voltage across the external load is linearly increased. The maximum converted energy of 59.6 mJ can be derived in 2 s under the vibration amplitude of 0.27 g.

Footnotes

Academic Editor: Xiaotun Qiu

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Natural Science Foundation of Hubei Province (2014CFB835).

References

Jin

Zhang

One-to-many chaotic synchronization with application in wireless sensor network. IEEE Commun Lett 2013; 17: 1782–1785.

Chew

Printed circuit board based memristor in adaptive lowpass filter. Electron Lett 2012; 48: 1610–1611.

Vilela

Forgie

. Miniature humidity micro-sensor based on organic conductive polymer—poly(3,4-ethylenedioxythiophene). Micro Nano Lett 2009; 4: 84–87.

Chew

Design and characterisation of a piezoelectric scavenging device with multiple resonant frequencies. Sensor Actuat A: Phys 2010; 162: 82–92.

Cook-Chennault

Thambi

Sastry

AM.

Powering MEMS portable devices—a review of non-regenerative and regenerative power supply systems with special emphasis on piezoelectric energy harvesting systems. Smart Mater Struct 2008; 17: 043001.

Roundy

Leland

Baker

. Improving power output for vibration-based energy scavengers. IEEE Pervas Comput 2005; 4: 28–36.

Lumentut

MF.

Mathematical dynamics of electromechanical piezoelectric energy harvesters. PhD Thesis, Curtin University, Bentley, WA, Australia, 2011.

Yang

Mechanics of piezoelectric structures. Singapore: World Scientific Publishing Company, 2006.

Yang

JS.

Variational formulation of the equations for small fields superposed on finite biasing fields in an electroelastic body. IEEE T Ultrason Ferr 2004; 51: 1030–1034.

10.

Pozzi

Zhu

ML.

Plucked piezoelectric bimorphs for knee-joint energy harvesting: modelling and experimental validation. Smart Mater Struct 2011; 20: 055007.

11.

Chidambaram

Mazzalai

Muralt

Measurement of effective piezoelectric coefficients of PZT thin films for energy harvesting application with interdigitated electrodes. IEEE T Ultrason Ferr 2012; 59: 1624–1631.

12.

Abas

Kim

Zhai

. Finite element analysis of vibration-driven electro-active paper energy harvester with experimental verification. Adv Mech Eng 2015; 7. DOI: 10.1177/1687814015571231.

Analytical simulation of the cantilever-type energy harvester

Abstract

Keywords

Introduction

Analytical model of piezoelectric unimorph cantilever energy harvester

Multimode frequency response of the unimorph cantilever beam

Case study of piezoelectric unimorph cantilever scavenger in bending mode

Analytical simulation

Numerical analysis

Transient analysis of the dynamic response

Conclusion

Footnotes

Declaration of conflicting interests

Funding

References