Analysis and Optimization of a Piezoelectric Energy Harvester

. The paper aims to assess and improve the performances of a multilayer piezoelectric MEMS device for vibrations harnessing. Two operating modes are possible: at resonance and outside resonance. In some applications it is not possible to operate at resonance, functioning being mostly at low frequencies in a quasi-static regime. An Euler-Bernoulli classic beam theory mathematical model was studied for estimating the behaviour of multilayer piezoelectric generators, in terms of deflection and voltage, at functioning under resonance frequency. The analytical results were compared with the finite element method simulation in COMSOL Multiphysics. The main goal of this study is to obtain an accurate model for engineering design purposes, with simple analytical equations and ease of use, but with predictable errors. The study proved the usefulness of the derived model but also its limitations. It also proves the need to improve the model using plate theory, for sensors with high width/height ratio.


Introduction
Different transduction mechanisms (piezoelectric [1][2][3][4], electromagnetic [4][5][6], electrostatic [7,8]) have been studied for converting vibrations into electricity, but piezoelectric energy harvesting has received the most attention due to the high power density and ease of application [1]. Hybrid harvesters using combinations of piezoelectricity with either of the other two methods have also been studied [9,10]. Piezoelectric materials have the property that, when subjected to mechanical strain, they produce an electrical charge proportional to the stress applied, unlike electromagnetic devices, which exhibit hysteretic responses.
Piezoelectric harvesters use the direct piezoelectric effect for converting mechanical vibrations into electrical energy, and tap this energy by connecting the harvester electrodes with an electric circuit [3]. For obtaining maximum electrical response, it is preferable to excite a given harvester at its fundamental frequency (first bending mode). At resonance, a small excitation can produce large vibration response, a sharp voltage peak occurring at that frequency. Operating at resonance enables energy harvesters to be placed on equipment producing high vibrations (engines, compressors, turbines, etc.).
Usually, a cantilever structure is preferred due to high electrical energy to transducer mass and volume ratio. To improve the output, an inertial mass is placed on the free end. Designing such a transducer involves deriving a lumped parameters model.
In other applications, the mechanical energy results from exerting a force on a point on the transducer. Such systems are currently integrated in suspensions and vibration isolation elements, under stairs, in highways, etc. The transducer does not support the external mass and force, this being the role of a flexible mechanical suspension, but is actuated during suspension movement. In this case the regime is quasi-static, the operating frequency being very low comparing with the resonance frequency.

Mathematical model
Simple, accurate and easily developable mathematical models are necessary for analysing piezoelectric structures. The best choice is performing a finite element method (FEM) simulation of a 3D model. This option gives very accurate results, however the main drawback is that electronic circuitry integration is not simple. Often it is more convenient to elaborate a concentrated parameters model of the assembly integrating the piezoelectric devices [2]. Previous studies have helped to derive simple and precise analytical models for a piezoelectric actuator and have also shown the modelling errors [11].
From engineering design point of view, approximate formulas are necessary, for ease of use during design and optimisation stages. These models can be validated by FEM simulations and experimental studies, as we did in the present paper.
The computations did not consider the materials yield strengths when applying the variable force. The purpose of the work herein was only of assessing the general behaviour of the harvester. The maximum voltages and deflections that can be obtained are limited by material yield strengths.
The simplest models are based on theory of elasticity hypotheses, which are applied to homogenous multilayered cantilever structures (Fig. 1 a)), with stress and strain exerted only on beam direction [12]:  Bending takes place on a normal plane to the maximum dimension of the beam;  The points on median line displace only on normal direction to this line. Euler-Bernoulli Beam Theory assumptions are also considered [12]:  A straight line or a cross-section, which is initially perpendicular on the median line of a layer, also remains perpendicular after deformation;  The perpendiculars on the median lines have the same rotation angles for all layers, resulting in a complete transmission of the specific deformations to the separation surfaces between layers. The usual equations are presented hereinafter.

A) Layers' constitutive piezoelectric equations
The system of equations for mechanical stresses and electrical displacements is: where: H (k-1) ≤ z ≤ H k ; T -Mechanical stress; Y -Young's modulus; S -Strain; s -Elastic constant; E -Electric field intensity; D -Electrical displacement; d -Piezoelectric charge/strain constant; ϵ -Permittivity; k -Layer number; x, z -Coordinates on x, z axes. Fig. 1. a) Multilayer element [12]; b) Deformation of infinitesimal element [12] Since only x direction components were taken into consideration for vectors and tensors, it is assumed that = .

B) Specific deformations and section displacements
For a beam with length dx, its cross-surface deformation occurs as in the schematic representation in Fig. 1 b). At coordinate, there is a line whose length dx stays constant after deformation (neutral fibre), the sectioning surfaces having the relative rotations .
We note = ∑ ℎ . The strain and stress in the fibre at z coordinate are:

C) Cross-surface equilibrium condition
The forces equilibrium condition on axial direction and the momentums equilibrium condition in report to the neutral axis, within the section at x coordinate, are: where: N -External reduced forces in cross-section; M -External reduced bending momentum in cross-section; ∆ ( ) -Electric potential of layer ( ); ( ) -Electrical potential on layers separation surface ( ).
Knowing the external axial force and bending momentum, the expressions for and can be inferred from equations (3) and (4).
In this case, equation (4) becomes: The momentum-(inverse of bending radius) constant K M and the piezoelectric momentum M piezo , calculated in relation to the neutral axis, are given by: Using the expressions for and , the electric displacement D and the electric charge Q accumulated on the electrodes can be calculated.
The charge accumulated on an electrode, for constant beam width b, and the electric displacement are: where: -electrical charge; -electrical displacement; -beam width; -charge/strain constant; -elastic constant; -Electrical field intensity; -Strain; -Permittivity. Superscript * has the significance of the value measured on the electrode.  [12] The stiffness constants are:

Trimorph using direct piezoelectric effect
where: ℎ = ℎ = ℎ , ℎ = ℎ and ( ) = ( ) = , ( ) = ; = 0, = ℎ , = ℎ + ℎ , = 2ℎ + ℎ . For trimorph, the electric induction in the vicinity of the superior electrode of layer 3 is: where: ℎ -height of piezoelectric layer; ℎ -height of central support layer; -rotation angle; -piezoelectric permittivity. The electric induction in the vicinity of the inferior electrode of layer 3 is: The overall electrical charge generated in superior layer 3 is: For trimorph with fixed constraint at one end and free at the other end, we obtain: where: l -beam length; b -beam width; -charge constant; Y -Young's Modulus; p -uniformly distributed pressure; F -Concentrated force at free end; M -concentrated bending momentum at free end; K -equivalent beam stiffness, given by equation (6). The total charge measured on the terminals depends on how the layers are connected, being = 2 for parallel connection (Fig. 2a) and = for series connection (Fig. 2b). For usual case = 0, = 0.

Finite element method simulation
In order to evaluate the accuracy of the models and to estimate the modelling errors, FEM simulations were performed in COMSOL Multiphysics, for a three-layered cantilever structure, in the configuration from Fig. 2 b). The FEM 3D model designed consists of a beam with two PZT-5H (Lead Zirconate Titanate) piezoelectric layers of 0.5 mm thickness, and a 1 mm thick central metallic layer, of a high-strength alloy steel.
Two sets of simulations were performed, static and dynamic. In both cases a cantilever structure was considered, with fixed constraint at one end and free in rest. In the case of static simulations, a variable force was applied, between 0 ÷ 30 N, uniformly distributed on the surface of the free end (values of 0 ÷ 200 N/m), perpendicular on the surface of the composite structure. We evaluated the displacement, the total voltage in the two piezoelectric layers and the maximum mechanical stress (von Mises stress and the component of the stress tensor on x direction), calculated as functions of force.
Numerical simulations using finite element method were conducted, assuming: -complete anisotropic 3D model; -plane stress: thin plate (less likely, assumes width much smaller than height); -plane strain: thick plate (more realistic, assuming width much larger than height).
Several options were also considered for declaring how the bar was clamped, both at the fixed and free end. The possible influence of motion and deformation blocking are relevant both for a correct design of structure coupling with neighbouring elements in the vibrating equipment, and also for simplifying the simulation model. The complete declaration of all components leads to a very high number of finite elements, and therefore to higher computational resources (in terms of processor, drive speed and very high memory requirement), leading to a high computing time.
The simulations were performed considering two intertwined physics (Fig. 3), namely:  Solid Mechanics, where layers materials are assigned, the material of layers 1 and 3 being declared piezoelectric. As well, fixed constraints and loads (force) are declared.  Electrostatics, where boundary conditions and initial conditions are declared and a conservation condition of the electric charge is introduced for layers 1 and 3. The two physics are coupled by activating the Multiphysics -Piezoelectric effect option, which allows simultaneous numerical solving of equations from the theory of elasticity, electrostatics and constitutive equations of the piezoelectric material.
In the case of 2D simulation in plane stress and plane strain hypotheses, similar models were used but the mathematical model of Solid Mechanics physics is simplified assuming a plane state of mechanical stresses (specific to a very thin structure) or a plane state of specific deformations (specific to a structure of infinite width). Hereinafter, considering a trimorph with ratio b/(2hp+hs)=10, the influence of clamping declaration was analysed for the fixed constrained end. Fastening between clamping plates (with condition of null displacements on exterior superior and inferior faces) was compared with the case when declaring null displacements for the cross-sectional surface of the three layered beam. The two possibilities were analysed, combined with the cases of cantilever with free or blocked tip.
We evaluated, for all the considered values of the force applied at the free end on a direction perpendicular to the surface (z-direction): the electric field in the piezoelectric layers, the voltage differences between central layer declared as ground and the upper and lower layer declared as terminals (based on which the total output voltage was computed), the deformations (based on which we computed the average deflection of the upper free end line on the superior edge), the distribution of stress tensor component on x direction evaluated in the piezoceramic layers, and the maximum overall stress in the piezoceramic layers, computed with von Mises relations (in the following figures, a force of 100 N/m).    between the analytic model and the 3D models (under 3%), higher differences appear when using 2D models (the plane strain model overrates by about 30% and the plane stress model underrates by about 10%). In the case of von Mises stresses (figure 8), the plane stress model overrates negligibly (below 2%) and the plane strain underrates by about 15%. For x-direction component of the mechanical stress tensor (figure 9), the 2D models (plane stress and plane strain) give similar results and very close to the theoretical model, but they are all significantly underrated comparing with the 3D models. These are grouped into models with a mobile tip subjected to deformations, with blocked curvature and respectively free. The last category leads to mechanical stress having very close value to the mechanical anisotropic model, about 20% smaller than the first group (which coincides with the mechanical isotropic model), and followed by about 20% smaller with 2D models.  These results lead to the idea of an underrating of the mechanical stresses in isotropic or theoretical models and in the 2D models, which are consistent with previous studies on modelling of piezoelectric actuators [11].
Differences can occur due to blocking free end deformation, but also due to not taking into account the plate assumptions. The current use of increasingly large piezoelectric structures, possible to manufacture as a result of technological advances, may amplify the errors due to neglecting the plate type models. Consequently, we analysed the influence of harvester width on the dimensions evaluated above, presented in the subsequent section.

Validity of beam models for large width structures
The first question is the independency of displacement in relation to beam width, assumed in the theoretical model. In the first case, the cantilever beam was considered fastened between two rigid parts in the clamping section, with free bending of the rest of the structure, on the in-plane directions y and x. The second case supposes constraining the bending of the free end, by means of other two bars. Mechanical load is defined as the constant per width (total force will increase with increasing width, but it should lead, according to the theoretical model, to constant values of tip deflection, voltage and mechanical stress).
The ratio between width b and the overall thickness (2hp+hs) was varied. A variation is obtained for the output voltage ( figure 10 a), the deflection at the free end ( figure 10 b), the maximum von Mises stress in the piezoelectric layers (figure 10 c) and the maximum x-axis component of the stress tensor in the piezoelectric layers ( figure 10 d). Theoretically, according to the linear model, all these parameters should be independent of beam width b.
Analysing the output voltage variation, we observed that the theoretical model is a good approximation for smaller width beams. For high width beams the voltages increase, with a higher value for blocked curvature of free end. The values of output voltage tend to stabilize as the width grows. The maximum voltages and deflections are limited by the lowest material yield strength.
The deflections of the free end on perpendicular direction on the trimorph were evaluated as average deflection of the superior edge. The deflections for constricted end were evaluated on the superior edge of the mobile end. We notice that for high b/(2hp+hs) ratio, in case of free curvature tip, the deflection value for the superior edge is higher than for blocked curvature tip. For small b/(2hp+hs), both are higher than the theoretical values, but tend towards these for a very wide beam.  High von Mises stresses occur mainly not due to the piezoelectric bending in xOz plane, but due to the clamping that constrains the bending in yOz plane [11].
There are significant differences between FEM simulations for small b/(2hp+hs) ratios. At high ratios, differences do exist, but very small, negligible, for von Mises stress, and smallest but not negligible for x-axis stress component. Nevertheless, both FEM simulations show differences more than 50-60% higher than the theoretical values calculated analytically with Euler-Bernoulli beam theory model.

Conclusions
The voltage values for 3D modelling are not dependent on the mechanical clamping mode and are comparable to the Euler-Bernoulli values for small width/height ratios, but however increase for large ratios, tending to stabilize as the ratio increases. Displacements for 3D modelling are higher than those for the Euler-Bernoulli model.
The mechanical von Mises maximum stresses in the piezoelectric layers for 3D simulations are comparable to the theoretical ones, but the component in the x-direction of the mechanical stress tensor depends on the free end clamping mode and is higher for 3D modelling than for the analytical model. Simplified 2D models lead to overrated stress (for plane strain) or significant stress (for plane stress). 2D models can be used for mechanical computations (displacement and von Mises stress), mainly for plane stress model.
At large beam widths, the behaviour of the harvester changes significantly: voltage increases and tends to stabilize, displacement decreases and tends to stabilize and the mechanical stress decreases with width. The theoretical model offers a good approximation of voltages at small widths but not for displacements and mechanical stresses at large widths. When not properly considering the sizes, the theoretical model underrates these parameters. An improved model, considering the plate theory, should be analysed.
The computations did not consider the yield strength limitations when applying the variable force, the purpose only being of assessing the general behaviour of the trimorph.