Modeling, Simulation and Optimization of Piezoelectric Bimorph Transducer for Broadband Vibration Energy Harvesting

We demonstrate the detailed analysis for conversion of piezoelectric properties into compliance matrix and simulate a series bimorph configuration for vibration based energy generation. Commercially available software COMSOL Multiphysics was used to apply boundary conditions for optimization of geometric parameters such as length, width and thickness of piezoelectric layer to study voltage and power characteristics of the harvester. The resulting energy harvester was found to generate 1.73 mW at 53.4 Hz across a 3MΩ load with an energy density of 13.08mJ/cm3. We also investigated feasibility of this model by comparing it with existing experimental data of known piezoelectric ceramic compositions and found good correlation between the two.


Introduction
Recent advances in micro-sensors technology have increased the need for cheaper and more efficient energy harvesting devices to provide on-board power.Literature suggests that many piezoelectric energy harvesting power generators have emerged in the past decade: Fang et al. (2006) studied a six-layer composite cantilever beam design, which has a 2.16μw electric power output under a resonance frequency of 608Hz.The natural frequency can be tuned by adjusting the weight of the tip mass.To overcome the issue of low power output level, Liu et al. (2008) studied a micro array of a multi-layered composite cantilever beam generating 3.98μw of effective electric power, yet the study did not show results on the micro-array system's frequency band coverage (although natural frequency of three individual beams is presented in the study).
It is well known that the voltage and power output is maximum when the piezoelectric device operates at it fundamental frequency.In this paper, we use COMSOL Multiphysics to determine the fundamental mode of piezoelectric device by Eigen frequency analysis.Power density decreases significantly when vibration frequency deviates away from the resonant frequency in Figure 1.In this research, we present study of detailed conversion of any polycrystalline anisotropic piezoelectric material properties to anisotropic compliance matrix that were directly used in simulation to obtain mechanical and electrical responses.

Nomenclature
Figure 1.Frequency-voltage response of the beam (PZT5H) Figure 2. Shoe energy harvester (Mitcheson et al., 2008;Mitcheson, 2015) Figure2 shows an energy-harvesting shoe design created by embedding a piezoelectric bender inside the sole; two piezoelectric unimorphs are attached to a curved plate, which is mechanically deformed under the footpad when a human's foot strikes on the ground.The unimorph beam is defined as a single-layer structure as shown in Figure 3.
Research shows that "2-8mW power can be obtained from vibration of a walking human" (Mitcheson et al., 2008).
In that application, a piezoelectric bender is suitable in harvesting ambient walking vibration energy because the weight of the device is much lighter than other types of energy harvesting devices such as electromagnetic generator, so that the extra weight of the piezoelectric bender inside the shoes will not interfere with the way humans normally walk.Using a unimorph also helps to reduce the weight of the piezoelectric bender, because a unimorph consists of only one layer of piezoelectric material, therefore it takes half of the weight of a bimorph configuration.
A bimorph beam is defined as a two-layer structure as shown in Figure 4.It is well known that a bimorph configuration can be used in series as well as in parallel based on electrical connection between the piezoelectric layers and the direction of polarization.Priya and Inman (2009) suggested a formula to calculate the voltage of the bimorph beam in series.Fang et al. (2006) used peak voltage and resistance to calculate average power.Roundy, Wright, and Rabaey (2003) derived a formula to calculate the optimal resistance in an electrical circuit.
The piezoelectric cantilever beam design highly depends on its engineering application in terms of vibration mode and number of layers of piezoelectric material.In the energy-harvesting shoe application, the weight of the device is an important design factor; therefore, the unimorph design was adopted to the transducer as shown in Figure 2; whereas in this research, achieving broadband response and maximizing the electrical voltage and power performances are essential in designing an energy harvester, therefore a bimorph configuration is used.Choosing the appropriate vibration mode is crucial in designing a high-performance energy harvester.There are two modes of operation commonly used in energy harvesting applications, namely the 3-3 mode and the 3-1 mode.The naming convention of the "A"-"B" mode indicates that the force is applied in "B" direction and electric charge appears in "A" direction.The x-direction in the Cartesian coordination system is commonly known as the "1" direction.The y direction is commonly known as the "2" direction.The z direction is commonly known as the "3" direction.
The 3-3 mode in unimorphs is shown in Figure 3.An external force (green arrow) is applied on the upper surface of the beam in the -z direction, which compresses the material and deforms the upper surface of the beam.Because of the compressive top-down force, the distance from the negative side to the positive side of individual unit cell gets smaller around the deformed area.Therefore, according to Eqn. (5) and Eqn.(6) for calculating electric dipole moments, an additional positive electric charge q (+ in red) is accumulated around the deformed area, and an equal amount of negative electric charge (-in blue) is distributed on the bottom surface of the unimorph.

Theoretical Analysis
To The piezoelectric moduli d ijk is a multi-linear relation between stress and electric charge.The piezoelectric moduli d ijk is a 3 rd rank tensor, which is also known as the coupling matrix, has the form of: The piezoelectric charge constant d 33 is higher than d 31 as shown in the coupling matrix (2).The piezoelectric charge constant d 33 is three-fold of d 33 in term of the magnitude for the first sample and is shown as follows: Same holds true for the second sample as follows 0 0 0 0 0 0 0 0 0 0 0 0 6.02 6.02 19.8 0 0 0 Although the piezoelectric charge constant d in the 3-3 mode is higher than that of the 3-1 mode, the 3-1 mode is more suitable in this research because of the vibration motion of the beam: strain built inside piezoelectric material along the "1" direction and electric charge accumulated in the "3" direction.In this research, a high-performance energy harvester was designed to scavenge low frequency mechanical vibration from the surrounding environment and the 3-1 mode has "lower stiffness which is suitable to detect small mechanical vibration" (Blevins & Plunkett, 1980).The resonant frequency of the cantilever beam which operates at the 3-1 mode is much lower than that of the 3-3 mode.Therefore, the 3-1 mode is a more attractive choice to harvest low frequency mechanical vibrations.Typically, in the 3-1 mode, the fundamental frequency of the cantilever beam is under 200Hz.The exact value of first mode fundamental frequency depends on many factors such as material composition of the beam and the geometry of the beam (length, width and thickness).The fundamental frequency can be analytically determined by = ( . ) (Bedekar, Oliver, & Priya, 2010;Bedekar, Oliver, & Priya, 2010) In Eqn. 4, E p is Young's modulus, I p is the rotational inertia, L is the length of the beam, M t is the tip mass, M b is the mass of the beam, and f r is the fundamental frequency.
Besides choosing an appropriate vibration mode of piezoelectric cantilever beam, the polarization is also one of the most important manufacturing process steps for piezoelectric material.The process is called "polling", which requires applying a very large electrical field (>2000V/mm) to align all the electrical dipoles in the same direction in the piezoelectric material.In fact, the polarization step is mandatory before the material can display piezoelectric properties.After the polling process is completed, the piezoelectric material has remnant electrical polarization.In COMSOL simulations, setting up the direction of polarization is achieved by defining a base vector coordination system, which is indicative of the polling direction.The electric dipole moment can be described by

= | |
(5) is the electric charge and is the distance vector, the total electric n dipole moments in one system can be expressed by which is a tensor product of piezoelectric moduli and the stress tensor .Piezoelectric moduli is a tensor of third order, which has 27, 3 3 , elements.However, many finite element simulation software like COMSOL, ANSYS, etc., use a reduced matrix form , which has 3 rows and 6 columns (18 elements), because of index condensation (due to the symmetry property of piezoelectric crystals).
The Eqns. 8,9,10 can be expanded to describe the direct piezoelectric effect in the xyz (1-2-3) direction: The bimorph design was chosen over the unimorph configuration for the superior electrical voltage and power performance, because a bimorph is a structure of two layers of material laminated together with conductive epoxy material on the horizontal centerline, as shown in Figure4.Under vibration mode, one layer expands and the other layer contract, and the directions of polarities are reversed in series connection, therefore the voltages are added up instead of cancelling out each other, as shown in Figure 4.The electrical charge is marked in red.Each individual electrical dipole moment is represented in circle.The direction of the electric dipole moment is from negative to positive, which is indicated by two arrows in Figure 4.The polarization directions are opposite in a bimorph structures.The tensile, compressive stress directions applied on the upper/lower layer are also opposite to each other.Therefore, the generated electric charges are distributed in the same direction ("+" and"-" are marked in red in Figure4).Negative charges are accumulated at the bottom surface of the beam(s).A similar effect from the simulation is shown in Figure 5, in which, red arrows show the polarization vector pointing up from simulation results, which supports the concept from Figure 4.This effect was also confirmed by the cross section of potential from simulation results as shown in Figure 5.

Analytical Modeling for Conversion of Piezoelectric Properties
Mechanical properties and electrical properties are dependent on each other and cannot be separated therefore, coupled equations or constitutive equations in strain-charge form are shown in equation 11,12: (Moheimani & Fleming, 2006) where D is the electric charge displacement, d is the piezoelectric coefficient (moduli) matrix, T is stress, ε is the permittivity matrix, E is electric field; S is mechanical strain, s is the compliance.Specifically, compliance s ijkl is a 4th rank tensor, which would have 81 elements.Due to symmetry, the compliance matrix has 6 rows and 6 columns, which can be written in matrix form as: Where G is the shear modulus.Notice that the compliance matrix must be a positive definite matrix, and all elements on the diagonal must be non-zero.

G = ( )
All values to the compliance matrix were substituted for PZTPZN scheme 4 sample as follows: 1 Therefore, mechanical-electrical coupling coefficient k 31 can be found from Eqn. 20.Finding coupling coefficient k is the key to find more unknown elements in the compliance matrix.The piezoelectric charge constant d, dielectric constant and compliance s can be found in Table 2 and Table 3.Based on Hooke's law, the relation between strain and stress is established by: = s (14) is strain, s is compliance of the material which is inverse of stiffness, and is stress.The relation between the negative ratio transverse strain and axial strain is expressed by the definition of Poisson's ratio as follows: The material is electrically poled along the z axis ("3" direction).The piezoelectric ceramics material is transversely isotropic material (Blevins, & Plunkett, 1980), therefore the transversely isotropic property was used to find more elements in the compliance matrix using the following equations: ) There is a more general form of calculation of as follows: By Eqn.17 and the transversely isotropic property, more compliance matrix elements can be found as follows: To find another important mechanical-electrical coupling coefficient k 33 , mechanical work and electrical work of the cantilever beam are needed as follows: 2 is the conversion ratio of mechanical work to electrical work.It can be expressed as: Taking the square root on both sides and substituting and , and Q=C p V, the following equation is obtained: Δ is the defection of the cantilever beam, F is the exciting force, Q is the electric charge accumulated on the surface of the beam, C p is the capacitance of the cantilever beam between upper surface and lower surface, V is the voltage of the cantilever beam between the upper surface and lower surface of the beam.
In the 3-1 mode, the beam can be simplified by considering the compliance on lower right elements in the compliance matrix.
S 44 = S 55 = S 33 (22) So far, all piezoelectric elements of the compliance matrix are obtained for the anisotropic model simulations.A well-ordered form is tabulated as shown in Table 4 and Table 5.

Validation Hypothesis of the Theoretical Model of COMSOL Multiphysics Simulation
To build the model correctly, besides applying the material condition, two kinds of boundary conditions were applied carefully to our simulation model: electrical boundary condition and structural boundary condition.The electrical boundary condition includes: connecting ground reference voltage to the bottom surface of the bimorph cantilever beam shown in Figure7; coupling an electrical circuit together with an external load resistor and connecting it to the upper surface of the bimorph cantilever beam shown in Figure7.The structural boundary condition includes a gravitational force on the beam and an exciting force 0.001N to vibrate the beam structure in the z direction as a boundary load condition; the fixed constraint on the end of the beam which is attached to a fixed end as shown in Figure 6.The exciting force at the free-end needs to be calculated and validated to make sure the exciting force cannot exceed the buckling force.The beam would buckle if the applied force exceeds buckling force.The buckling force is defined as the maximum force that can apply on the beam before non-linear behavior of the beam appears.The calculation of the buckling force is expressed as follows: = (23) F b is the buckling force, is density, Y is Young's modulus, fr is fundamental frequency, is mechanical quality factor, V is the volume of the cantilever beam, which can be calculated as follows: V = l w t (24) Where, l is the length of the cantilever beam, w is the width of the cantilever beam, t is the thickness of the cantilever beam.To calculate voltage, a thin geometry parameter in length (40mm), width (2mm) and thickness (0.2mm) of the cantilever beam is considered, because as can be seen from Eqn.23, volume is positively proportional to the magnitude of the buckling force.Keeping all other parameters constant, the smaller the volume of the beam, the lower the buckling force of the cantilever beam.
Given the density, Young's modulus, fundamental frequency, mechanical quality factors of the first sample and the second sample, the buckling forces were obtained by Eqn.23 for both samples: 5.79N and 30.4N.The exciting force 0.001N is far less than the buckling force of the first sample and the second sample.Therefore, the magnitude of the exciting force is reasonable.In other words, the beam will not collapse under the load boundary condition.Alternatively, the bucking force can also be obtained by running "linear buckling" analysis in COMSOL Multiphysics.In summary, the appropriate vibration mode, the magnitude of loading force and numbers of layers of the cantilever beam, and the two new PZT-PZN material compositions are chosen to explore and design a low frequency piezoelectric energy harvester.The results reflect changes in the length of the beam from 40mm to 60mm with increments of 4mm, width of the beam from 2mm to 18mm with increments of 4mm, thickness from 0.2mm to 0.5mm with increments of 0.05mm.

Results and Discussion
The shift of open circuit voltage can be seen in Figure 8.The open circuit voltage is negatively proportional to the thickness of the bimorph on the upper right subplot; as the beam gets thicker the beam become stiffer.In other words, the maximum mechanical deformation of the beam decreases with increasing cross section area.
It can be seen in Figure 9 that the longer bimorph has a lower resonance frequency (left); the thicker bimorph has a higher resonance frequency (right) as the natural frequency of the beam is also inversely proportional to the thickness of the beam In the following section, it is shown how the analytical resonant frequency changed when length, width, and thickness were changed individually.
Analytical resonance frequency is determined by young modulus E p , rotational momentum I p , length L, the mass of beam M b and tip mass M t , (Bedekar, Oliver, & Priya, 2010;Kim, Tadesse, Priya, & Inman, 2009) It can be seen from the upper left subplot to the right subplot in Figure9 that the longer bimorph has lower analytical resonance frequency (left).As shown in Eqn.25, the resonance frequency is negatively correlated to the cube of the beam length-the denominator grows faster than the numerator as the beam gets longer.Therefore, the resonance frequency decreases as the length gets longer, and the thicker bimorph has higher resonance frequency (right subplot).A similar effect appears on the unimorph (Kim et al., 2008).The thickness of the beam has great linear impact on the frequency of the cantilever.It is concluded from Figure 9 that the resonance frequency is positively proportional to the thickness of the beam; the length is negatively proportional to the resonance frequency; the width of the beam has very little significance on the resonance frequency.It can be seen from the upper left subplot to the right subplot in Figure 11 that the longer bimorph had higher power output, the wider bimorph had lower power output, and the thicker bimorph had lower power output.When the length is 60mm, the width is 2mm, the thickness is 0.2 mm, the maximum power of the first sample obtained was 1.75× 10 watts; the maximum power of second sample obtained was 3.14× 10 watts.Figure 11 shows how power changed when length, width, and thickness changed individually.The maximum power density of the first PZT-PZN sample is 7.28mW/cm 3 and, that of second PZT-PZN sample is 13.08mW/cm 3 .The power density results agree with the published results of H. Kim et al.: power density of piezoelectric micro-energy-harvesting devices ranges between 10 to 400 μW/cm 3 in the acceleration range of 2.3 to 78.4 m/s 2 and the frequency range of 0.08 to 1 kHz (Kim & Priya, 2008).It is clear from the results that the power generated is directly proportional to the length of the beam whereas it's inversely proportional to the width and thickness of the beam.The ultimate purpose of designing a high-performance energy harvester is to provide power to electrical components, such as resistors.The optimum resistance was calculated based on maximum power generation.In practical scenario, the device will be operated at natural frequency of the beam at an optimized resistance for its geometry so that the device always generates "maximum" power.As the red curve shows (in Figure12), when the resistance reached 5× 10 ohms, the power reached maximum value (0.4 mW).As mentioned in the earlier section, a fixed length of the beam was set at 60mm, width changes from 2mm to 18mm, 4mm interval, thickness changes from 0.2mm to 0.5 with 0.05mm interval, thus making it a total of 35 permutations.These 35 cases were selected of distinct geometry of the cantilever beam and its corresponding vibration to 35 different geometry at resonant frequency.Therefore, 35 curves with 35 distinct colors were obtained and shown in Figure12 and Figure13.In Figure13, all thirty-five curves showed a consistent pattern: when the load resistance goes up, the voltage saturates at some value.In Figure12, power goes up at the beginning because the voltage goes up and when the voltage becomes saturated at the open circuit voltage, power goes down because the load resistance goes up continuously.The electric power generated by a vibrating beam can be approximated at natural frequency by (Kim & Priya, 2008): m is the mass of the beam, Y is the amplitude of the vibration, ω is the resonance frequency, is the damping constant (Kim & Priya, 2008), as the equation 26 shows that the power of the beam vibrating (P res ) at resonance frequency increases when the damping constant decreases.The damping constant is inversely proportional to mechanical quality factor: The damping constant and mechanical quality factor is tabulated in Table 5: The PZT-PZN samples' mechanical quality factors Q m is higher than that of PZT5A/H, which indicates that the PZT-PZN sample has lower damping, therefore the power of the beam vibrating (P res ) at resonance frequency is higher using PZT-PZN samples than that of traditional PZT5A/H piezoelectric material.
The efficiency η of mechanical energy to electrical energy conversion can be approximated by: Where, k is the electrical coupling factor that appears on the denominator and numerator of the equation 28: Q m is mechanical quality factor, when Q m increases, efficiency η increases.The product of piezoelectric voltage constant and charge constant g 31* d 31 is reported to be indicative of the piezoelectric energy harvesters' performance (Kim & Priya, 2008).However, the result in this research does not support the theory (g * d is indicative of power): the performance of PZT-PZT samples outperform PZT5A/H based on the result of this research is shown in Table 6 and Table 7.
Since to the method to find the "optimal resistance" is well-known by varying the external resistance, it is natural to consider the question: how to find the optimal resistance changes when length, width, and thickness are changed individually.The longer bimorph has higher optimal resistance; the wider bimorph has lower optimal resistance.The thicker bimorph has little effect on optimal resistance, as it is seen in Figure14, from the upper left subplot to right subplot.
Figure 14.Optimal resistance response in terms of variation of beam's dimension (l,w,t)

Experimental Validation
The experimental power and voltage data of PZT-PZN (scheme 4) piezoelectric 3-layer bimorph cantilever beam (25mm x 5.5mm x 0.4mm for single PZT-PZN layer, 0.05mm thickness Brass) is compared against the simulated average voltage and average power data, as shown in Figure 15 and Figure 16: As can be seen in Figure 15, the value of simulated voltage (using our model) and measured voltage are fitted well.It can also be seen in Figure16 that the maximum average power is at the same level between simulated average power and measured average power.Figure15 and Figure16 together show that our simulation model is validated by the experimental measured data.
To complement simulation work, the research group made PZT5A-brass (C22000)-PZT5A and PZT5H-brass(C22000)-PZT5H sandwiched beams.The width of the beam is 10mm and the length of the beam is 60mm, the thickness of the composite beam is shown in Figure 17:

Conclusions
Our main contributions in this research are to use newly developed PZT-PZN material and converting its piezoelectric material properties into the compliance matrix which is suitable for COMSOL Multiphysics simulations.A list of relevant piezoelectric material parameters of given samples are selected and converted to construct compliance matrix, and anisotropic simulation models are built.From Table 8, the results of simulations show that this research is accurate and promising in modeling piezoelectric energy harvesters by validating between the simulated and measured experimental average voltage and power performance of PZNPZN scheme 4 material.
It is interesting to note that when the length of the beam increases, the electric power output increases, whereas the power decreases with increasing width and thickness of the beam.In future work, the investigators will work on optimization of geometry and material properties to provide a broadband energy harvesting response which will enable the device to operate in off-resonance condition and still be able to generate sufficient power for charging small scale electronics such as cell-phone, smart watches.
mechanical work W E : electrical work G: Shear modulus Q: electric charge F: applied force

Figure 5 .
Figure 5. Polarization vector of bimorph cantilever bending downward simulation (3-1 mode) mechanical-electrical coupling coefficient k , piezoelectric charge constant d and dielectric constant and compliance s is established by = (Priya & Inman, 2009)

Figure 8 .
Figure 8. Sample1's open circuit voltage response in terms of variation of beam's dimension (l,w,t)

Figure 10 .
Figure 10.Sample1's analytical resonance frequency response in terms of variation of beam's dimension(l,w,t)

Figure 11 .
Figure 11.Sample1's power response in terms of variation of beam's dimension (l,w,t)

Figure 15 .
Figure 15.RMS voltage vs load resistance of PZT-PZN (scheme 4) cantilever beam in measured and simulated form at resonance frequency at 183Hz

Figure 17 .
Figure 17.Thickness view of the upper and lower PZT5 layer is 0.3mm, and middle layer (C22000 brass) is 0.025mm

Table 2 .
Piezoelectric Material properties values

Table 3 .
Elements of compliance matrix of sample one

Table 4 .
Elements of compliance matrix of sample two

Table 5 .
Comparison of four different piezoelectric materials between the quality factor, damping constant and product of voltage and charge constant

Table 7 .
Comparison of TWO different piezoelectric materials

Table 8 .
Maximum power Comparison of four different piezoelectric materials