Design optimization of vibration energy harvesters fabricated by lamination of thinned bulk-PZT on polymeric substrates

Andrés Vásquez Quintero; Nadine Besse; Pattanaphong Janphuang; Robert Lockhart; Danick Briand; Nico F de Rooij

doi:10.1088/0964-1726/23/4/045041

1. Introduction

The increasing use of portable and ultra-low power smart wireless devices consuming less than 100 μW continues to boost interest in energy harvesters that could provide autonomy for extended periods of time. They could either recharge or completely replace electro-chemical energy sources such as batteries, thus reducing maintenance/replacement costs, and enabling the realization of smart systems placed at remote and/or isolated locations. The piezoelectric principle for energy conversion is favored for micro-sized devices, due to its achievable small dimensions and its high power density [1–5]. Among different transduction mechanisms, resonant vibration harvesters are widely investigated in the literature due to a greater power generation from ambient sources [5]. Compact and functional MEMS piezoelectric energy harvesters based on thin and thick piezoelectric materials have been demonstrated.

Thanks to MEMS fabrication processes and the relatively high electromechanical coupling, different active piezoelectric materials have been integrated onto rigid silicon cantilevers; of those, one can cite sputtered aluminum nitride (AlN) [6–8], deposited thin film lead zirconate titanate (PZT) [9–12], screen-printed thick PZT [13, 14], and thinned bulk-PZT [15, 16]. Besides the silicon industry, the plastic and printed electronics industry is starting to establish itself. A new class of smart sensing objects or large area surfaces is emerging thanks to the use of polymeric foil with unique characteristics. Energy components such as flexible solar cells, batteries, and supercapacitors are currently under development [17–19], with some commercial examples already available. The realization of mechanical energy harvesters on flexible or elastomeric substrates has gained a lot of attention recently. The fabrication processes of vibration harvesters on flexible substrates using screen-printed PZT on polyimide [20], coated polyvinylidene difluoride (PVDF) on PET substrates [21, 22], and the transfer of bulk-PZT onto PET [23] have been reported. The latter takes advantage of the excellent piezoelectric properties of commercially available bulk piezoelectric ceramics. More recently, interest has increased in fully flexible energy harvesters employing bending modes instead of resonant modes, to be used in low-frequency and on-body applications. Such configurations typically use zinc oxide (ZnO) and PZT as vertically aligned nanowires (NWs) [24–26], transfer-printed nano-ribbons [27], and transfer-printed PZT onto polydimethylsiloxane (PDMS) substrates [28].

In this work, we combine the superior piezoelectric properties of bulk-PZT [16, 29, 30] with PET polymeric substrates towards the realization of vibration energy harvesters to achieve high power densities at low frequencies (<50 Hz). A cost-effective, foil-based, and low-temperature fabrication process was developed which is based on the lamination of thinned bulk-PZT (50 μm) onto a flexible PET substrate at low temperatures. Models are required to geometrically optimize the harvester structure in order to maximize its output power and further reduce its resonance frequency. Some groups have already presented analytical models describing the mechanical and electrical behavior of vibration cantilever harvesters [31–37], accounting for point proof masses and the back-coupling effect. However, the combination of the back-coupling effect with a distributed model of the proof mass to determine the mechanical and the electrical harvester performance has not been addressed yet. The two models presented here account for a distributed mass at the free end, the back-coupling effect, and the additional stress-charges generated in the anchor region. The first model, referred to as the hybrid model, is based on the Euler–Bernoulli beam theory used by Ertuk [31–33]. Additionally, it is further extended by considering the distributed mass used in the multi-span theory [38], and combining it with the back-coupling effect and FEM simulations proposed individually by other authors [34, 37]. The presented model allows for optimization of the mass/beam ratio, which was shown to have an important influence in the output power and resonance frequency of the harvester.

The second, referred to as the numerical model, was developed entirely using FEM simulations in COMSOL, integrating the piezoelectric and the electric circuit modules. Both models (hybrid and numerical) were validated with experimental data. Finally, in addition to the geometrical optimization of the structure for maximum output power generation at a minimized frequency, the electrical and mechanical characterization of fabricated optimal devices are presented and compared to the developed models.

2. Design and methodology

The vibration energy harvester presented is a unimorph cantilever beam with a rectangular cross-sectional area and a proof mass located at the free end. It consists of a piece of thinned bulk-PZT (d₃₁ mode) laminated onto a flexible polymeric substrate (polyethylene terephthalate—PET) at low temperature (<85 °C). Figure 1 provides the design and the dimension nomenclature which will be used throughout this document. Two modeling approaches were followed to design and optimize the harvesters in terms of power at both low frequency and low acceleration.

**Figure 1.** Schematic representation of the vibration energy harvester, illustrating the multi-layer structure and the dimension nomenclature.
Download figure:
Standard image High-resolution image

The hybrid model is a combination of analytical equations calculated using MATLAB and a simulated correction factor (α) from COMSOL simulations. The correction factor is included to account for the extra stress-charges generated in the anchor electrode area (A) that are excluded by the hybrid model, which assumes perfect clamping at one end of the beam (figure 2(a)).

**Figure 2.** Schematic representation of (a) the *hybrid model* and (b) the *numerical model* using analytical and FEM modeling, respectively.
Download figure:
Standard image High-resolution image

The numerical model calculates all of the output variables using only COMSOL simulations, which allows evaluation of more complex geometries and clamping conditions. This model implements clamping boundaries that are more realistic by including an anchor region connected to a rigid substrate (figure 2(b)). The anchor electrode area (A) directly includes both the extra capacitance and the additional stress-charges generated.

Following the fabrication of the first generation of harvesters used to validate the models and to quantify the damping coefficients, a geometrical optimization was performed that considered the length, width, and thicknesses of each of the layers to achieve high power at low frequencies. Based on the optimal dimensions, a second generation of devices was fabricated and compared once more to the developed models to corroborate their validity and the design approach.

3. Materials

The vibration energy harvester is composed of an active layer of bulk lead zirconate titanate (PZT) laminated onto a flexible polyethylene terephthalate (PET) substrate (Melinex^® ST506, 125 μm thick from DuPont^©). Both layers were assembled together using the dry film photoresist (DP) PerMX3014 (14 μm thick) from DuPont. The PZT layer is PSI-5H4E^® (127 μm thick), from Piezo Systems^©, Inc. (Massachusetts, USA), which comes in the form of a sheet (7 × 7 cm²). Vacuum sputtered nickel (Ni) electrodes were pre-deposited on both sides (100 nm thick). The piezoelectric ceramic presents extremely low current leakage and low magnetic permeability, with a relative dielectric constant of 3800 (when mechanically free). The PZT thickness was thinned down from 127 μm (initial value) to 75 μm and finally 50 μm in order to enhance the output power and to reduce the resonance frequency for the same lateral device dimensions. The dry adhesive ARClear 8932^®, from Adhesives Research Inc., was used both to attach the fabricated harvesters onto rigid copper-clad glass-reinforced epoxy (FR-4) substrates for characterization purposes, and to fix the proof mass (a tungsten (W) sheet 500 μm thick) to the beam. The electrical contacts to the Ni electrodes were made using silver-filled epoxy glue (EPO-TEK H70E^®), from Epoxy Technology Inc., and wire soldering.

4. Modeling and experimental methods

4.1. Hybrid model

The hybrid model combines analytical equations and a correction factor from COMSOL simulations. The analytical side associates already existing theories into a more complete vibration energy harvester model. The equations describe the mechanical and electrical behavior of the harvester, combining the back-coupling effect [31] and the distributed mass model at the free end [38]. This implementation is not trivial, since the mass locally modifies the neutral plane, the moment of inertia, and thus the deflection profile of the beam [38].

The analytical development is based on Euler–Bernoulli beam theory [31, 32], which describes the deflection of the beam as a function of position and time in response to an external force [33, 39–41]. The cantilever is excited by both the movement of its base and the inertia of the beam and mass combined [31, 42]. Besides enlarging the beam deflection, the mass creates a high stress region along its inner edge, improving the charge generation. All the equations were implemented and numerically solved using MATLAB scripts for static and harmonic acceleration regimes in the vertical direction (z-axis).

To begin with, the multi-layer cantilever was transformed into an equivalent two-span beam system which accounts for the distributed mass at the free end [38], as shown in figure 3. This transformation is the core of the analytical approach dividing the beam into two blocks (b1, b2) in the vertical direction at the point where the mass starts. Each segment is represented by its equivalent bending stiffness (EI_1,2), calculated from the material properties and their respective dimensions, as follows [43]:

$\begin{eqnarray}\label {eq1} h_{1,2}&=&\frac {\sum {E_{{\rm MAT}} z_{{\rm MAT}} A_{{\rm MAT}}} }{\sum {E_{{\rm MAT}} A_{{\rm MAT}}} }\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray} \label {eq2} I_{{\rm MAT}1,2}&=&\frac {w_{{\rm MAT}} t_{\rm MAT^{3}}}{12} + (h_{1,2}-z_{{\rm MAT}})^{2} A_{{\rm MAT}}\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray} \label {eq3} EI_{1,2}&=&\sum {E_{{\rm MAT}1,2} I_{{\rm MAT}1,2}}, \end{eqnarray} \tag{ 3 }$

where h_1,2 is the neutral plane with respect to the z-axis for each block, E_MAT, A_MAT, w_MAT, and t_MAT are the Young's modulus, the cross-sectional area, the width, and the thickness of the respective materials, z_MAT is the distance from the bottom plane to the center of each material, and I_MAT1,2 is the inertia for each material in each of the blocks.

**Figure 3.** Schematic of the equivalent Young's modulus and the moment of inertia between the harvester layers and the modeled multi-span beam.
Download figure:
Standard image High-resolution image

Using the multi-span beam system, the general Euler–Bernoulli motion equation set is presented as follows [31]:

$\begin{eqnarray}&&\label {eq4} EI_{1,2}\frac {\partial ^{4}z_{1,2}(x,t)}{\partial x^{4}} + m_{1,2}\frac {\partial ^{2}z_{1,2}(x,t)}{\partial t^{2}} \nonumber \\ &&+\ c_{1,2}\frac {\partial z_{1,2}(x,t)}{\partial t} + K_{1,2} U(t)=f_{1,2}(x,t), \end{eqnarray} \tag{ 4 }$

where z_1,2(x, t), c_1,2, m_1,2, f_1,2(x, t), K_1,2U(t) are the vertical displacement as a function of horizontal position and time, the damping, the equivalent mass, the external forces, and the back-coupling effect, for b1 and b2, respectively. Although other models have already been proposed in the literature[41, 44], in the present work we have developed a model which combines a distributed mass approach with the back-coupling effect. The back-coupling effect describes the extra stiffness of the PZT beam during operation, due to an opposing bending force caused by the stress-charges generated in the active layer [45].

The method of separation of variables was used to solve equation (4) for the vertical displacement (Z_1,2(x, t)), as shown in equation (5).

$\begin{eqnarray}&&\label {eq5} {Z}_{1,2}(x,t)=X_{1,2}(x ) T(t), \end{eqnarray} \tag{ 5 }$

where X_1,2(x) and T(t) are the vertical displacement functions with respect to space and time, respectively.

The homogeneous Euler–Bernoulli equation (neglecting the damping, the back-coupling, and the external forces) was solved in order to find the natural frequency of the system and later to calculate the general displacement equation of the beam. The spatial coordinate and the characteristic equations are then defined as follows:

$\begin{eqnarray}&&\label {eq6} \frac {\partial ^{4}X_{1,2}(x)}{\partial x^{4}}-\beta _{1,2}^{4}{ X}_{1,2}(x)=0\end{eqnarray} \tag{ 6 }$

$\begin{eqnarray} &&\label {eq7} \beta _{1,2}^{4}=\left (\frac {\lambda _{1,2}}{l_{B1,2}} \right )^{4}=\frac {m_{1,2}}{EI_{1,2}} \omega _{{\rm n}}^{2}\end{eqnarray} \tag{ 7 }$

$\begin{eqnarray} &&\label {eq8} \omega _{{\rm n}}^{2}=\left (\frac {\lambda _{1,2}}{l_{B1,2}} \right )^{4}\frac {EI_{1,2}}{m_{1,2}}, \end{eqnarray} \tag{ 8 }$

where β_1,2 is the wavenumber, ω_n is the circular frequency, l_B1,2 is the length of the beam, and λ_1,2 is the eigenvalue, for b1 and b2, respectively.

The characteristic equation (6) is solved for the spatial composition (X_1,2(x)), which gives the generic equation (9), including eight coefficients: A_1,2, B_1,2, C_1,2, and D_1,2. The coefficients were found using the boundary conditions for the system [46], shown in the set of equation (10).

$\begin{eqnarray}&&\label {eq9} X_{1,2}(x)=A_{1,2}(\cos \langle \beta _{1,2}x \rangle +\cosh \langle \beta _{1,2}x \rangle ) \nonumber \\ &&+\ B_{1,2}(\cos \langle \beta _{1,2}x \rangle -\cosh \langle \beta _{1,2}x \rangle )\nonumber \\ &&+\ C_{1,2}(\sin \langle \beta _{1,2}x \rangle +\sinh \langle \beta _{1,2}x \rangle )\nonumber \\ &&+\ D_{1,2}(\sin \langle \beta _{1,2}x \rangle -\sinh \langle \beta _{1,2}x \rangle ) \end{eqnarray} \tag{ 9 }$

$\begin{eqnarray} &&\begin{array}{l}\label {eq10} X_{1}(0,t)=0 \quad\qquad X_{1}^{\prime }(0,t)=0 \\ X_{2}^{\prime \prime }(l_{B2},t )=0 \quad\qquad X_{2}^{\prime \prime \prime } (l_{B2},t )=0 \\ X_{1}(l_{B1},t )=X_{2}(0,t) \quad\qquad X_{1}^{\prime }(l_{B1},t )=X_{2}^{\prime }(0,t) \\ EI_{1}X_{1}^{\prime \prime }(l_{B1},t )=EI_{2}X_{2}^{\prime \prime }(0,t)\\ EI_{1}X_{1}^{\prime \prime \prime }(l_{B1},t )= EI_{2}X_{2}^{\prime \prime \prime }(0,t). \end{array} \end{eqnarray} \tag{ 10 }$

Substituting the boundary conditions in the spatial equation (9), the set of eight coefficients can be found for each block: b1 equations (11) and b2 equations (12).

$\begin{eqnarray} &&\begin{array}{l}\label {eq11} A_{1}=0; \quad\qquad B_{1}=1; \quad\qquad C_{1}=0 \\ D_{1}={-B}_{1}\left [(\cos \langle \beta _{1}l_{B1} \rangle {\vphantom {\frac {\beta _{1}^{2}EI_{1}}{\beta _{2}^{2}EI_{2}}}}- \cosh \langle \beta _{1}l_{B1}\rangle )\right .\\ \quad \times \ (\sin \langle \beta _{2}l_{B2}\rangle - \sinh \langle \beta _{2}l_{B2}\rangle )\\ \quad +\ \displaystyle \ \frac {\beta _{1}^{2}EI_{1}}{\beta _{2}^{2}EI_{2}} (\cos \langle \beta _{1}l_{B1}\rangle +\ \cosh \langle \beta _{1}l_{B1}\rangle )\\ \quad \times \ (\sin \langle \beta _{2}l_{B2}\rangle -\sinh \langle \beta _{2}l_{B2}\rangle )\\ \quad +\ \displaystyle \frac {\beta _{1}}{\beta _{2}}(\sin \langle \beta _{1}l_{B1}\rangle +\sinh \langle \beta _{1}l_{B1}\rangle )\\ \quad \times \ (\cos \langle \beta _{2}l_{B2}\rangle -\cosh \langle \beta _{2}l_{B2}\rangle )\\ \quad \displaystyle +\ \frac {\beta _{1}^{3}EI_{1}}{\beta _{2}^{3}EI_{2}} (\sin \langle \beta _{1}l_{B1}\rangle -\sinh \langle \beta _{1}l_{B1}\rangle )\\ \quad \left .\times \ (\cos \langle \beta _{2}l_{B2}\rangle + \cosh \langle \beta _{2}l_{B2}\rangle ){\vphantom {\frac {\beta _{1}^{2}EI_{1}}{\beta _{2}^{2}EI_{2}}}} \right ]\\ \quad \times \ \left [(\sin \langle \beta _{1}l_{B1}\rangle {\vphantom {\frac {\beta _{1}^{2}EI_{1}}{\beta _{2}^{2}EI_{2}}}} -\sinh \langle \beta _{1}l_{B1}\rangle )\right .\\ \quad \times \ (\sin \langle \beta _{2}l_{B2}\rangle +\sinh \langle \beta _{2}l_{B2}\rangle )\\ \quad \displaystyle +\ \frac {\beta _{1}^{2}EI_{1}}{\beta _{2}^{2}EI_{2}} (\sin \langle \beta _{1}l_{B1}\rangle + \sinh \langle \beta _{1}l_{B1}\rangle )\\ \quad \times \ (\sin \langle \beta _{2}l_{B2}\rangle -\sinh \langle \beta _{2}l_{B2}\rangle )\\ \quad \displaystyle -\ \frac {\beta _{1}}{\beta _{2}}(\cos \langle \beta _{1}l_{B1}\rangle -\cosh \langle \beta _{1}l_{B1}\rangle )\\ \quad \times \ (\cos \langle \beta _{2}l_{B2}\rangle -\cosh \langle \beta _{2}l_{B2}\rangle )\\ \quad \displaystyle -\ \frac {\beta _{1}^{3}EI_{1}}{\beta _{2}^{3}EI_{2}} (\cos \langle \beta _{1}l_{B1}\rangle +\cosh \langle \beta _{1}l_{B1}\rangle )\\ \quad \left .\times \ (\cos \langle \beta _{2}l_{B2}\rangle +\cosh \langle \beta _{2}l_{B2}\rangle ){\vphantom {\frac {\beta _{1}^{2}EI_{1}}{\beta _{2}^{2}EI_{2}}}} \right ]^{-1} \end{array}\end{eqnarray} \tag{ 11 }$

$\begin{eqnarray} &&\begin{array}{l}\label {eq12} \displaystyle A_{2}=\left (\frac {1}{2}\right )[B_{1}(\cos \langle \beta _{1}l_{B1}\rangle -\cosh \langle \beta _{1}l_{B1}\rangle )\\ \quad +\ D_{1}(\sinh \langle \beta _{1}l_{B1}\rangle -\sinh \langle \beta _{1}l_{B1}\rangle ) ] \\ \displaystyle B_{2}=\left (\frac {1}{2}\right )\frac {\beta _{1}^{2}EI_{1}}{\beta _{2}^{2}EI_{2}}[B_{1}(\cos \langle \beta _{1}l_{B1}\rangle +\cosh \langle \beta _{1}l_{B1}\rangle )\\ \quad +\ D_{1}(\sinh \langle \beta _{1}l_{B1}\rangle +\sinh \langle \beta _{1}l_{B1}\rangle ) ] \\ \displaystyle C_{2}=\left (\frac {1}{2}\right )\frac {\beta _{1}}{\beta _{2}}[{-B}_{1}(\sin \langle \beta _{1}l_{B1}\rangle +\sinh \langle \beta _{1}l_{B1}\rangle )\\ \quad +\ D_{1}(\cos \langle \beta _{1}l_{B1}\rangle -\cosh \langle \beta _{1}l_{B1}\rangle ) ]\\ \displaystyle D_{2}=\left (\frac {1}{2}\right )\frac {\beta _{1}^{3}EI_{1}}{\beta _{2}^{3}EI_{2}}[{-B}_{1}(\sin \langle \beta _{1}l_{B1}\rangle -\sinh \langle \beta _{1}l_{B1}\rangle )\\ \quad +\ D_{1}(\cos \langle \beta _{1}l_{B1}\rangle +\cosh \langle \beta _{1}l_{B1}\rangle )]. \end{array} \end{eqnarray} \tag{ 12 }$

Using the Euler–Bernoulli equation along with the orthogonality condition for natural modes [38], and assuming a sinusoidal acceleration a applied on the piezoelectric energy harvester, the amplitude T of the vertical displacement with respect to time T(t) is given by equation (13) [31].

$\begin{eqnarray}&&\label {eq13} T= \frac {\gamma -\varphi U}{\omega _{{\rm n}}^{2}-\omega _{{\rm e}}^{2} +j 2 \zeta \omega _{{\rm n}} \omega _{{\rm e}}}, \end{eqnarray} \tag{ 13 }$

where

$\begin{equation}\label {eq14} \gamma =\frac {\sum {m_{1,2} a} \int _0^{l_{B1,2}} {X_{1,2}(x) \, {\rm d} x} }{\sum {m_{1,2} } \int _0^{l_{B1,2}} {X_{1,2}^{2}(x ) \,{\rm d} x} } \end{equation} \tag{ 14 }$

$\begin{equation} \label {eq15} \varphi =\frac {\sum {K_{1,2} \int _0^{l_{B1,2}} {\frac {{\rm d} ^{2}X_{1,2}(x)}{{{\rm d} x}^{2}}\, {\rm d} x} } }{\sum {m_{1,2} \int _0^{l_{B1,2}} {X_{1,2}^{2}(x) \,{\rm d} x} } } \end{equation} \tag{ 15 }$

$\begin{equation} \label {eq16} K_{1,2}=-E_{{\rm PZT}} d_{31} \left (\frac {t_{{\rm PZT}}}{2}-h_{1,2}\right ). \end{equation} \tag{ 16 }$

Considering the electrical circuit as a capacitor (representing a piezoelectric energy harvester having no internal resistance) connected to an external resistor R, and assuming that all the current I_c generated by the PZT flows into the resistor R, the voltage U across the resistive load can be calculated by equation (17).

$\begin{eqnarray}\label {eq17} U&=&R I_{{\rm c}}=j{ \omega }_{{\rm e}} R{ Q}_{{\rm c}}\nonumber \\ &=&j{ \omega }_{{\rm e}} R [(Q_{d1}+Q_{d2} )-C U], \end{eqnarray} \tag{ 17 }$

where Q_d1,2 are stress-charges generated along the surface of the PZT and C is the internal capacitance of the PZT layer, given by equations (18) and (19).

$\begin{eqnarray}&&\label {eq18} Q_{d_{1.2}}=-\frac {\alpha }{t_{{\rm PZT}}} \int _0^{l_{B1,2}} \int _0^{w_{B}} \int _{-h_{1,2}}^{t_{{\rm PZT}}-h_{1,2}} d_{31}\nonumber \\ &&\times \ E_{{\rm PZT}} z \frac {{\rm d} ^{2}Z_{1,2}(x,t)}{{\rm d} x^{2}} \,{\rm d} x \,{\rm d} y \,{\rm d} z\end{eqnarray} \tag{ 18 }$

$\begin{eqnarray} &&\label {eq19} C=\varepsilon _{0} \varepsilon _{33} \frac {w_{B} (l_{A}+l_{B})}{t_{{\rm PZT}}}, \end{eqnarray} \tag{ 19 }$

where α is the correction factor extracted from COMSOL simulations. It is defined as the sum of all the stress-charges generated within the PZT (including the extra charges) divided by the sum of the stress-charges created in the cantilever (without the anchor) (refer to section 6.1.2).

Finally, inserting equation (13) into equation (17) and rearranging the terms, the average generated power was obtained by solving equation (20).

$\begin{eqnarray}&&\label {eq20} P=\frac {\vert U \vert ^{2}}{2 R}\nonumber \\ &&=(R{ \omega }_{{\rm e}}^{2} (q_{d1}+q_{d2} )^{2} \gamma ^{2}/2)\nonumber \\ &&\quad \times ([\omega _{{\rm n}}^{2}-\omega _{{\rm e}}^{2}-2 \zeta \omega _{{\rm n}} \omega _{{\rm e}}^{2} R C ]^{2}\nonumber \\ &&\quad +\ [\omega _{b} R (q_{d1}+q_{d2} ) \varphi \nonumber \\ &&\quad +\ (\omega _{{\rm n}}^{2}-\omega _{{\rm e}}^{2} ) \omega _{{\rm e}} R C+2 \zeta \omega _{{\rm n}} \omega _{{\rm e}} ]^{2})^{-1}, \end{eqnarray} \tag{ 20 }$

where

$\begin{eqnarray}\label {eq21} q_{d_{1.2}}&=&\frac {Q_{d_{1.2}}}{T} =-\frac {\alpha }{t_{{\rm PZT}}} \int _0^{l_{B1,2}} \int _0^{w_{B}} \int _{-h_{1,2}}^{t_{{\rm PZT}}-h_{1,2}} d_{31}\nonumber \\ &&\times \ E_{{\rm PZT}} z \frac {{\rm d} ^{2}X_{1,2}(x)}{{\rm d} x^{2}} \, {\rm d} x\, {\rm d} y \, {\rm d} z. \end{eqnarray} \tag{ 21 }$

Equation (20) uses the calculated stress-charges for each block q_d1 and q_d2 with respect to space, the estimated optimal load R (found at the maximum power by sweeping a resistive load and the resonance frequency), the measured damping coefficient ζ, the theoretical natural frequency ω_n, and the excitation frequency of the base ω_e for a given acceleration a.

4.2. FEM simulations

Two finite element simulations were implemented in COMSOL Multiphysics v4.2 using the clamping conditions shown in figure 2(b). The first uses the piezoelectric module of COMSOL to calculate the correction factor α emulating the total generated stress-charges in the cantilever beam. As discussed previously, this factor is then implemented in the hybrid model. The second simulation defines the numerical model, combining the piezoelectric and electric circuit modules to introduce a resistive load connected across the electrodes of the PZT material. This combination enables the simulation of the output power generated by the system as a function of load, acceleration, damping, and excitation frequency, as well as geometrical parameters.

In both cases, the material models used were as follows: piezoelectric material for the PZT layer and linear elastic material for the others (PET substrate, dry film photoresist, tungsten mass, and rigid fixation substrate). In order to get accurate results within a reasonable calculation time, the mesh was optimized depending on the different materials and dimensions. The mesh was defined with the mapped and distributed functions of COMSOL using square-shaped elements. The PZT layer has the smallest mesh size (0.1 mm × 0.2 mm × 1 μm for the length, width, and thickness) due to its electromechanical behavior, while the proof mass and the rest of the structure have a mesh size of 0.1 mm × 0.2 mm × 20 μm. The designs were simulated using the 3D geometry representation, and the mechanical boundary conditions were set such that the base was rigidly anchored while the rest of the structure remained free. In the second simulation, only the PZT bottom electrode was electrically set as ground while the top electrode was connected to a variable resistor through the terminal boundary condition. The metal conductors and the adhesion between layers were assumed to be ideal. The type of materials and their corresponding mechanical and electrical properties used for the analytical and FEM approaches are shown in tables 1 and 2, respectively.

Table 1. Mechanical properties of the materials used for the analytical calculations and the FEM simulations.

Mechanical properties

Material	E (GPa)	Density (kg m⁻³)	Poisson's ratio (—)
PZT-5H	62	7 800	0.27
PET-Melinex	4.6	1 430	0.33
DP-PerMX	2	1 430	0.33
Tungsten	411	19 350	0.28
PCB	22	1 900	0.28

Table 2. Electrical properties of the PZT-5H used for the analytical calculations and the FEM simulations.

Electrical properties (PZT)

Piezo. coefficients (pC N⁻¹)		Rel. perm.
d₃₁	−320	ε₁₁	3130
d₃₃	650	ε₂₂	3130
d₁₅	741	ε₃₃	3800

4.3. Validation procedures

The same validation procedure was followed for the hybrid model and the numerical model. After defining the geometry, the mechanical/electrical properties, and the boundary conditions, the first mode eigen-frequency for the respective design was calculated/simulated. Next, a frequency sweep was performed at short- and open-circuit conditions within a range that includes the eigen-frequency previously found. In [47], it was shown that the losses due to damping are mainly caused by internal structural friction with low influence of the surrounding air (no significant differences were observed between measurements under vacuum and air environments). The effect of air damping was tested by measuring the damping coefficient of ten devices with different dimensions (areas ranging from 30 to 65 mm²). The damping coefficient was shown to be 0.012 ± 0.002, which represents a relatively small variation, and therefore it was assumed to be constant over this range of dimensions. In addition, a linear relationship between the damping coefficient and the input acceleration (range from 0.3 g to 1 g) was observed with an incremental ratio of around 0.0037 g⁻¹. Finally, a parametric sweep study was implemented to find the optimal resistive load value that gives the maximum output power. Using equations (1)–(21), and in parallel the frequency domain sweep simulation in COMSOL, the stress-charges at the surface, the voltage, and the power at the optimal load were calculated or simulated, respectively. It should be noted that this procedure can be performed including or excluding the proof mass, a process which allows one to compare its influence on the mechanical and electrical dynamic behavior of the cantilever stack.

The characterization procedure followed for the fabricated devices is based on [48], where an automatic test bench is used to measure values such as resonance and anti-resonance frequency, coupling coefficient, power as a function of frequency (short and open circuit), power as a function of resistive load (at fixed frequency and acceleration), optimal load, power, voltage, bandwidth, and tip displacement.

After validation, the analytical model was used to optimize the geometry of the harvester in order to obtain the maximum output power at low frequencies. Due to the relatively fast feedback given by the MATLAB scripts, the beam length (0–30 mm), the beam width (2–5 mm), the PZT/PET length ratio (0–1), the PZT thickness (10−150 μm), and the mass/beam length ratio (0–1) were swept, and optimized values were determined.

5. Fabrication

An optimized vibration energy harvester fabrication is presented. The process is based on lamination steps in order to bond bulk-PZT sheets to PET substrates and printing techniques to make the electrical contacts [23]. Besides targeting low complexity, the process is developed at low temperature due to the relatively low glass transition temperature of PET, ≈80 °C. Figure 4 shows the main fabrication steps, indicating some typically used parameters for the geometry and the type of material.

To enhance the output power and to achieve low resonance frequencies, the PZT layer was thinned down at wafer level using a grinder DAG810^® from Disco^© [15]. The latter process uses a rotating wheel with diamond sweeps that removes material from the PZT sheet at very precise vertical displacements (figure 4(a)). The grinding speed was 30 μm min⁻¹, with a total thickness variation of around 1 μm and a final roughness at the PZT surface of less than 100 nm. Next, the PZT was diced with a diamond saw according to the required dimensions. In parallel, the PET substrates were laser ablated to the respective dimensions with a Speedy 300 fiber^TM CO₂ laser from Trotec^® (figure 4(b)).

The transferring process starts with the lamination of the dry film photoresist (DP) PerMX (14 μm thick) onto the diced PZT beams (figure 4(c)). The lamination was performed using a semi-automatic double-roller laminator at 85 °C, 2 bar, and 2 m min⁻¹. After a post-lamination bake (5 min @ 85 °C), the stack was exposed to UV light (250 mJ cm⁻²) from the PET side in order to polymerize the DP [49]. Finally, a post-exposure bake (10 min @ 60 °C) was performed in order to finalize the cross-linking process. Next, the PZT/DP stack was aligned and bonded to the PET cantilever substrate using 8 kPa of pressure at 85 °C for 10 min (figure 4(d)). In order to obtain a good adhesion between the layers and to achieve the flexible mechanical properties of the DP, the stack was hard cured at 85 °C for 1 h in a convection oven [49]. To facilitate the characterization experiments, the PZT/DP/PET stack was fixed to a FR-4 rigid substrate using the dry adhesive AR8932. The top and bottom contacts were connected using stencil-printed silver-filled epoxy adhesive (EPO-TEK H70E), from the top to the copper-clad PCB and through the dry adhesive vias, respectively (figures 4(e) and (f)). The conductive adhesive was cured in a convection oven for 2 h at 85 °C. Finally, a laser-ablated tungsten mass was fixed at the free end of the beam using the dry adhesive.

6. Results and discussion

6.1. Hybrid model (analytical and FEM based)

6.1.1. Two-span beam displacement continuity

The displacement of each block (b1, b2) was calculated from equation (9), solving for X₁(x) and X₂(x), the deflections of the two blocks along the x-axis. The solution of these two spatial displacements was obtained using the boundary conditions presented in equation (10), which describes the interaction of the two blocks with respect to each other and the clamping. Using the mentioned displacement profiles, X₁(x) and X₂(x), the resonance frequency, the strain, the stress-charges, the voltage, and the power were obtained according to equations (13)–(21).

Figure 5 displays the calculated displacement values along the x-axis for each block, in order to show the continuity of the equations at their interface (x = 9.5 mm). The insets of figure 5 show a schematic and the tip displacement as a function of frequency of a device with a distributed mass at the free end (distributed on 50% of the total beam length) and at the optimal load. The continuity corroborates the validity of the boundary conditions for the equivalent blocks.

6.1.2. Additional charge simulation

The hybrid model uses a mechanically fixed boundary condition along the beam thickness (at the PCB frame) and not at the base of the anchor point (as in reality) due to the relatively lower complexity. As shown in equation (21), the additional stress-charges are corrected with the factor α from COMSOL simulations. This factor depends on the mechanical properties (i.e. flexural rigidity) of the beam. Since the same materials, cross-section, and fixation technique were used in all of our devices, the correction factor is almost constant. Figure 6 shows the first resonant mode (eigen-frequency simulation at optimal load) solved for the von Mises stress of the structure and for the current density of the PZT at the anchor point and at the interface with the mass (close-up view). It should be noted that the maximum stress region (at the anchor point) corresponds to the maximum current density. The latter is an indication of the stress-induced charges in the PZT. These simulations show that the stress-charges generated increase by 10–15% when the anchor region is included and a flexible PET substrate is used.

**Figure 6.** First resonance mode of the multi-stack harvester. The color legend of the whole structure represents the simulated von Mises stress (red to blue range), while the color legend of the maximum stress regions represents the simulated current density (black to white range).
Download figure:
Standard image High-resolution image

6.1.3. Model validation

The validity of the hybrid model for a range of accelerations (0.3–1 g) and cantilever areas (30−70 mm², with a linear increase of the damping coefficient with respect to the acceleration) was proven by comparing calculated and experimental data of the power generated at different vibrating frequencies. The optimal load used for the power calculation was found by sweeping it from 100 Ω (short circuit) to 1 MΩ (open circuit). Figures 7(a)–(d) compare the same device: with (A) and without (B) mass at 0.3 g, 0.5 g, 0.7 g, and 1 g, respectively. The solid lines represent the values obtained from the hybrid model, while the dotted lines denote the experimental data. The measured dimensions of devices A and B are given in table 3. Figure 7 shows a good correlation between the calculated and measured curves, validating the model at least for the mentioned ranges of acceleration and dimensions. The average differences for all accelerations (0.3 g, 0.5 g, 0.7 g, and 1 g) between the calculated and measured values for output power and resonance frequency were 13 ± 5%, 1.5 ± 0.5% with mass and 15 ± 6%, 4.5 ± 0.4% without mass, respectively. By simulating and comparing the output power versus frequency curves with and without the back-coupling effect (neglecting the coupling factor K_1,2 in equations (13)–(16)), it was found that the power and resonance frequency were just slightly decreased and up-shifted, respectively. In [47], it was shown by our group that the back-coupling effect in Si-based devices was relatively higher than in polymer-based ones, probably due to the higher electromechanical coupling coefficient.

**Figure 7.** Comparison of calculated (solid lines) and measured (dotted lines) power versus frequency for a device with (A) and without (B) mass, at (a) 0.3 g, (b) 0.5 g, (c) 0.7 g, and (d) 1 g.
Download figure:
Standard image High-resolution image

Table 3. Geometrical parameters used for the simulated and fabricated devices with the hybrid model and the numerical model.

Material	Width (mm)	Length (mm)	Thickness (μm)
PZT-5H	2.0	19.5	70
PET-Melinex	2.0	19.5	125
DP-PerMX	2.0	19.5	14
Tungsten	2.2	10.0	500
Extra PZT (l_a)	2.0	4.0	70

6.1.4. Geometrical optimization

The geometrical optimization was performed using the relatively fast feedback from MATLAB scripts of the hybrid model. Figure 8 shows the curves of power (left axis) and resonance frequency (right axis) versus beam length, beam width, PZT thickness, the PZT/PET length ratio, and the mass/beam length ratio.

A longer beam increases the power and reduces the frequency (equation (8)). However, it could affect the mechanical integrity of the beam by causing fractures due to both an increase in mass and larger deflections. Furthermore, a wider beam augments the output power due to an increase in area and mass of the device, whereas the resonance frequency remains constant provided that geometry of the structure is a cantilever and not a plate (i.e. length/width ratio ≫ 1). Here, a trade-off is identified between the size of the device and the output power, and therefore the optimum value depends on the application of interest (figures 8(a) and (b)).

Full coverage of PZT over the PET cantilever (i.e. PZT/PET ratio = 1) was found to give the maximum power for constant thicknesses, as shown in figure 8(c). The frequency presents a maximum value at about 60% coverage due to the higher relative stiffness of the stack (i.e. the Young's modulus of PZT is higher than that of PET). For larger coverage ratios, the frequency decreases, since the effect of the additional mass is more significant than the effect of its rigidity.

In addition to the PZT/PET length ratio, their respective thicknesses play an important role in the harvester optimization due to the location of the neutral plane. In [23], it was shown that the power can be maximized by placing the neutral plane as close as possible to the interface between the PZT and the PET layers. Appropriate manipulation of the neutral plane can be achieved by either increasing the thickness of the PET or by decreasing the thickness of the PZT layer, due to their respective Young's moduli (E_PZT > E_PET). In this work, the PZT thickness was decreased, following the grinding procedure described previously, in order to keep the overall dimensions of the design as small as possible. Additionally, the resonance frequency is directly proportional to the thickness of the PZT layer; therefore, reducing the thickness reduces the resonance frequency of the cantilever.

Using the hybrid model, the neutral plane, the power, and the resonance frequency were calculated while sweeping the PZT thickness from 10 to 150 μm and keeping the PET thickness constant at 125 μm. The maximum output power was found for a PZT thickness of 50 μm (figures 8(d) and (e)). Here, also, the resonance frequency tends to increase for thicker layers, due to the increased stiffness of the beam.

Finally, the mass/beam length ratio was increased from the free end of the beam to the anchor point fully covering the beam (i.e. mass/beam ratio = 1). The power and the resonance frequency of the harvester were calculated over this range, as shown in figure 8(f). The output power reaches a maximum value and the resonance frequency falls to a minimum value for a ratio of about 0.5, due to the trade-off between the stiffness and the mass of the cantilever.

6.2. Numerical model (entirely FEM based)

6.2.1. Model validation

The validity of the numerical model implemented in COMSOL for the range of accelerations (0.3–1 g) and cantilever areas (30−70 mm² with a linear increase of the damping coefficient with respect to the acceleration) was proven by comparing the simulations to the experimental data representing the output power at different frequencies. The optimal load used for the power simulations was found after sweeping the connected resistor from 100 Ω (short circuit) to 1 MΩ (open circuit). Figures 9(a)–(d) compare the same device with (A) and without (B) mass at 0.3 g, 0.5 g, 0.7 g, and 1 g, respectively. The solid lines represent the values simulated with the numerical model, while the dotted lines denote the experimental data. The dimensions and measured characteristics of the devices are found in table 3. Figure 9 shows a good correlation between the simulated and measured curves, validating the model at least for the mentioned ranges of accelerations and dimensions. The average differences for all accelerations (0.3 g, 0.5 g, 0.7 g, and 1 g) between the calculated and measured values for power and resonance frequency were 16 ± 8%, 1.3 ± 0.6% with mass and 7 ± 4%, 1.4 ± 0.3% without mass, respectively. It should be noted that for a device with a mass, excited at relatively high accelerations (e.g. 1 g), the measured power curve (figure 9(d)) is not symmetric, as a result of the nonlinear response of the PZT under large stress [50]. The use of nonlinear elastic mechanical properties and damping coefficients has been proposed to overcome these limitations [51]. In its present state, these effects are not taken into account in the hybrid model or in the numerical model, which limits their respective accuracy.

**Figure 9.** Comparison of FEM-simulated (solid lines) and measured (dotted lines) power versus frequency for a device with (A) and without (B) mass, at (a) 0.3 g, (b) 0.5 g, (c) 0.7 g, and (d) 1 g.
Download figure:
Standard image High-resolution image

6.2.2. Optimized harvester

Two optimal devices, known as harvester I (long) and harvester II (short), were fabricated following the results and design recommendations from the harvester geometrical optimization (figure 8). The dimensions were as follows: length 19.3 mm (device I) and 15 mm (device II), width 4 mm, PZT/PET length ratio 1, PZT thickness 50 μm, and mass/beam length ratio 0.5. Harvester I without mass is depicted in figure 10(a), while the optimized one (including the tungsten proof mass) is shown in figure 10(b). A cross-sectional view of the multi-layer stack is presented in figure 10(c), indicating the respective materials and thicknesses. The final dimensions of the devices are reported in table 4.

**Figure 10.** Pictures of fabricated optimal harvester I: (a) side view of the harvester mounted on a PCB without the mass, (b) side view of the harvester with the mass compared to a button battery, and (c) cross-sectional view of the harvester showing the PET, DP, PZT, and W layers.
Download figure:
Standard image High-resolution image

Table 4. Geometrical parameters of the optimized and fabricated harvester.

Material	Width (mm)	Length I/II (mm)	Thickness (μm)
PZT-5H	4.3	19.3/15	52
PET-Melinex	4.3	19.3/15	125
DP-PerMX	4.3	19.3/15	14
Tungsten	4.3	9.7	505
Extra elec.	4.0	4.0	70

The numerical model simulations (solid lines) and measured values (dotted lines) for output power were compared and are found to be in good agreement at the optimal load for optimized harvesters I and II: with (figures 11(a), 12(a)) and without (figures 11(b), 12(b)) the proof mass. Figure 11(a) shows the output power curves for accelerations from 0.1 to 0.6 g. The input acceleration was limited to 0.6 g, since above this value the beam deflection was greater than the allocated space (≈2.5 mm) and could induce damage of the device. The damping value at 0.1 g was calculated from the linear regression of damping versus acceleration (0.0037 g⁻¹) mentioned previously. Figure 11(a) shows output powers of 88.7 μW at 48.6 Hz and 6.7 μW at 49.8 Hz for optimized harvester I with an optimal load of 75 kΩ at 0.5 g and 0.1 g, respectively. As pointed out previously, the model excludes the nonlinear effects which appear for large beam deflections. The calculated dynamic electromechanical coupling coefficient $k_{2}=(f_{{\rm oc}}^{2}-f_{{\rm sc}}^{2})/f_{{\rm sc}}^{2}$ (where f_oc and f_sc are the open-circuit and short-circuit resonance frequencies) was calculated for the optimized harvesters I and II, resulting in values of 2.03% and 1.8%, respectively.

**Figure 11.** Comparison of numerical FEM-simulated (solid lines) and measured (dotted lines) power versus frequency for optimized harvester I with (a) and without (b) mass at different accelerations.
Download figure:
Standard image High-resolution image

Relatively high values for normalized power density (NPD) [52] of 11 683 μW g⁻² cm⁻³ and volume displacement figure of merit (FOM) [48] of 227 μW g⁻² cm⁻³ were found for optimized harvester I. The volume displacement FOM was calculated using FOM = P × g⁻² × (V_Σ)⁻¹ × Q⁻¹, where P is the output power, g is the input acceleration, Q is the quality factor, and V_Σ = w_B × l_B × (t_PET + t_DR + t_PZT + t_W + d_max), with d_max as the maximum peak-to-peak displacement of the cantilever beam. Optimized harvester I outperformed the previous version [47] in terms of output power (at 0.1 g), NPD (including the empty volume between the proof mass and the anchor region of the beam), and FOM (including the beam displacement) by a factor of 5.1, 2.46, and 2.3, respectively. Figure 11(b) shows the maximum generated power of 12 μW, a resonance frequency of 185 Hz, and an optimal load of 15 kΩ at 1 g without proof mass.

Figure 12 shows a good agreement between the measured harvester (harvester II) and its numerical simulations, resulting in a maximum power of 48 μW at 110.5 Hz and 32 kΩ of optimal load.

Additionally, figure 13 shows the average current and voltage curves for optimized harvester I (with and without mass). It is pointed out that the average current and voltage values at the resonance frequency, optimal load, and including the proof mass, are 40 μA and 3.1 V, respectively.

A comparison of the power as a function of load curves between the numerical model simulations (solid lines) and measured data (dotted lines) was performed to completely corroborate the validity of the model. Figure 14(a) shows the output power as a function of load at 0.1 g and a fixed frequency of 49.8 Hz (the natural frequency of the harvester), while figure 14(b) shows the power curves and frequency shifts (additional stiffness due to the back-coupling effect) for several resistive loads (from 1 to 500 kΩ). For the measured data as well as for the FEM simulations, the optimal load was found to be between 80 and 100 kΩ, providing a maximum power of 6.7 μW at 0.1 g. It should be noted that the system generates considerable power levels at relatively low frequencies (≈50 Hz) and at low input accelerations, which was one of the main reasons to use the polymeric PET substrate.

**Figure 14.** Comparison of FEM-simulated (solid lines) and measured (dotted lines) power versus load for an optimized device with mass at 0.1 g of acceleration for optimized harvester I; (a) power versus load and (b) power versus frequency at different loads (from 1 to 500 kΩ).
Download figure:
Standard image High-resolution image

The overview and comparison of optimized harvester I with other published harvesters is summarized in table 5. As we can see, the NPD of harvester I outperformed the others due to the relative higher displacements thus higher power per acceleration. However, higher flexibility causes reduction of the FOM, which takes into account tip displacements.

Table 5. Performance comparison between the optimized harvester presented in this work (harvester I) and two other relevant devices from literature.

Device	PZT type	Freq. (Hz)	Power (μW g⁻²)	Current (μA g⁻¹)	Voltage (V g⁻¹)	R_l(kΩ)	Disp. (mm)	Total volume (mm³)	k²(%)	NPD (μW g⁻² cm⁻³)	FOM (μW g⁻² cm⁻³)
This work (I)	Thinned bulk PZT-5H	49.8	670	100	5	75	2	57.3	2.03	11 683	227.4
[16] (Michigan)	Thinned bulk PZT-5A	167	274	—	11.7	500	—	27	—	10 150	—
[53] (DTU)	Screen-printed PZT	235	14	—	1.5	100	—	12.9	—	1 087	—

7. Conclusion

The modeling and the fabrication of an optimized thinned bulk-PZT-based vibration energy harvester on a plastic polymeric substrate (PET) is presented in this paper. A hybrid model (analytical + FEM) and a numerical model (only FEM) were validated and shown to be in good agreement with the experimental data collected from fabricated devices. The hybrid model, based on the Euler–Bernoulli beam theory, the multi-span theory and enhanced with COMSOL simulations, was used to geometrically optimize the harvester. The power was maximized for a low-frequency operation (below 100 Hz) through a parametric study. The optimized device was in good agreement with the developed models, generating a power of 6.7 μW at 49.8 Hz, 0.1 g across a load of 75 kΩ. It exhibited a normalized power density NPD of 11 683 μW g⁻² cm⁻³ and a figure of merit (FOM) of 227 μW g⁻² cm⁻³ (including displacement). In addition, an output power higher than 100 μW was obtained for accelerations above 0.5 g.

Acknowledgments

This work was partly funded by the EU FP7 project FlexSmell, a Marie Curie Initial Training Network (ITN), under grant no. 238454. The authors thank Dr Rolanas Dauksevicius at EPFL-IMT SAMLAB for fruitful discussions, and Dr Bill A McDonald from DuPont^® Teijin Films and Dr Nikki Birkbeck from Adhesives Research, for the supply of Melinex^® ST506 and ARClear 8923 samples, respectively.

Design optimization of vibration energy harvesters fabricated by lamination of thinned bulk-PZT on polymeric substrates

Article metrics

Submit

Permissions

Author e-mails

Author affiliations

Dates

Peer review information

Abstract

1. Introduction

2. Design and methodology

3. Materials