High accuracy comsol simulation method of bimorph cantilever for piezoelectric vibration energy harvesting Open Access

Piezoelectric vibration energy harvesting technology is promising in power self-supply for IoT (Internet of Things) devices. It uses piezoelectric inertia devices to vibrate resonantly with the ambient vibration frequency and convert the vibration energy into electrical energy based on the piezoelectric effect. Piezoelectric bimorph cantilever is a typical vibration energy harvesting structure. It consists of an upper and lower piezoelectric layer and an intermediate substrate layer. Under the excitation of ambient vibration, the transverse resonance of the composite cantilever produces large longitudinal stress in the upper and lower piezoelectric layers, and the piezoelectric material generates electricity in d₃₁ mode.

In order to design piezoelectric bimorph cantilever and analyze its electromechanical coupling characteristics, many modeling methods have been proposed, such as distributed parameter method (DPM),¹ finite element method,² and lumped parameter method.³ 

DPM is a more accurate modeling method, which usually uses Newton’s or Hamilton’s principles⁴ to obtain coupled vibration equation. The analytical solution⁵ can be obtained by combining the mode shapes of Euler-Bernoulli beams with piezoelectric constitutive equations, and the numerical solution⁶ can be obtained by Rayleigh Ritz method. Modeling process of DPM is complex for irregular beam and cantilever with additional mass.⁷ It is difficult to obtain analytical solutions of modal modes. Kim⁸ investigated the distributed parameter model of piezoelectric cantilever with additional mass. Shengxi Zhou⁹ investigated the distributed parameter model for the zigzag piezoelectric cantilever structure.

Finite element analysis (FEA) is a numerical simulation method, which is more suitable for solving complex beam model. Based on Kirchhoff plate theory, Carlos De Marqui Junior modeled and simulated the piezoelectric energy harvester plate by FEA.¹⁰ It is time-consuming and labor-consuming to compile the FEA program. At present, the main commercial FEA softwares, such as ANSYS, Abaqus and COMSOL, all support piezoelectric coupling analysis of structures.

FEA has the advantage in the electromechanical coupling modeling and simulation of complex topological structures. Vinayaga K.¹¹ discussed the ANSYS modeling and simulation method of smart piezoelectric beam. Yang¹² used FEA to optimize the structural parameters of a clamped piezoelectric circular diaphragm. Qian¹³ used ANSYS to carry out modal analysis of a torsional piezoelectric vibration energy harvesting system. Fey¹⁴ used ANSYS to model and simulate a piezoelectric auxetic PZT lattice structure. Yang¹⁵ used ANSYS to model and simulate a two-dimensional piezoelectric energy harvester. Zhou¹⁶ used ANSYS to simulate a flexible longitudinal zigzag structure. Gohari¹⁷ used quadratic piezoelectric multi-layer shell element for analysis of smart laminated composite plates in ABAQUS. Bath¹⁸ used COMSOL to model and simulate 3D folded zigzag piezoelectric energy harvester. Dauksevicius¹⁹ used COMSOL to model and simulate a nonlinear piezoelectric vibration energy harvester with frequency-tuned impacting resonators.

In addition to electromechanical coupling simulation, FEA can also simulate the coupling of magnetic field and piezoelectric structure vibration, which is conducive to the design of nonlinear vibration system. Upadrashta²⁰ used ANSYS to model and simulate nonlinear piezoelectric energy harvesters with magnetic interaction. Dauksevicius²¹ used COMSOL to analyse magnetic plucking dynamics in a frequency up-converting piezoelectric energy harvester.

Another advantage of FEA is the whole system modeling and simulation of electromechanical coupling structure and circuit. Wu²² used COMSOL to analyze electrically rectified piezoelectric energy harvesters. Elvin N. G⁶ established a coupled finite element-circuit simulation model for analyzing piezoelectric energy generators in ABAQUS.

Although FEA has been widely used in the modeling and simulation of piezoelectric vibration energy harvester (PVEH), and some studies claim that the COMSOL simulation results fit well with the experimental results, but the detailed parameter values and modeling methods are not given, and the simulation is difficult to reproduce.

It is still necessary to analyze some basic and detailed problems, such as getting the appropriate modeling parameters for FEA, improving the accuracy of calculation and reducing the calculation time. In addition, the results of FEA simulation should be compared with the theoretical model with distributed parameters and the experimental results to verify the reliability of FEA.

This study will discuss a FEA method for the piezoelectric bimorph cantilever by COMSOL Multiphysics, aiming at structural eigenfrequency analysis and electromechanical coupling analysis. Compared with the DPM simulation results of Erturk A,²³ the reliability of FEA in COMSOL is verified.

The typical piezoelectric bimorph cantilever for energy harvesting is shown in figure 1. The piezoelectric ceramic layers of the piezoelectric bimorph cantilever are tensioned and compressed axially (in the direction of piezoelectric material d₃₁) by inertia force under transverse vibration or rotational excitation of the fixed-end base. Mechanical energy is converted into electric energy based on the piezoelectric effect. Adding tip mass at the free end can effectively reduce the resonant frequency of the device, and increase inertia and piezoelectric output. The upper and lower piezoceramic layers of bimorph can be connected in series or in parallel. Series type can produce larger voltage output and parallel type can produce larger current output.

A distributed parameter electromechanical coupling model given by Erturk A²³ well described the vibration energy harvesting characteristics of piezoelectric bimorph cantilever. The electromechanical frequency response functions were given both in series type and in parallel type respectively. DPM model reflecting the internal relationship between the steady-state electrical and mechanical responses and the excitations of base translational and rotational accelerations. For the common first-order resonance excitation problem, the multi-modal solution is simplified to a single-modal form. When the base is only excited by translational acceleration, the single mode voltage frequency response function is:

where in this example the first modal mass term γ₁ is calculated as -0.0224. The first modal electromechanical coupling term χ₁ is calculated as 0.0842 in series type and 0.1684 in parallel type, respectively. The internal capacitance term C_p is calculated as 7.5661 nF in series type and 30.264 nF in parallel type, respectively. The calculation method can be seen in Ref. 23.

For the first-order eigenfrequency analysis of piezoelectric bimorph cantilever, the composite cantilever beam can be equivalent to a homogeneous cantilever beam. The theoretical derivation of the bending section coefficient YI of composite cantilever beam can be expressed as follows:

For a homogeneous cantilever beam, its inertial moment is:

Therefore, the equivalent elasticity modulus of homogeneous cantilever beam is:

The equivalent density of homogeneous cantilever beam is:

The first order eigenfrequency of the theoretical cantilever beam is:

Here, λ₁ = 1.875, m is the mass of the cantilever beam of unit length. The explanation of parameters in all above formulas can be seen in Tables I–III.

An example of piezoelectric bimorph cantilever for theoretical analysis and experimental verification has been given by Erturk A.²³ In this paper, the piezoelectric bimorph cantilever with the same geometry and material parameters as the example is used, then its 1D (dimension), 2D and 3D finite element models are built in COMSOL, respectively. The eigenfrequency and voltage frequency response functions of the model are simulated, and the model parameters are discussed.

PZT-5H is chosen as piezoelectric material, brass as substrate material, and nickel electrode as top electrode and bottom electrode of piezoceramics. The thickness of nickel electrode is negligible in the simulation. As shown in Table I, structural geometry modeling parameters can be parameterized in the global definition of COMSOL.

According to the above geometric design parameters, 1D, 2D and 3D finite element parametric modeling are carried out in COMSOL as shown in figure 2. A 1D homogeneous Euler-Bernoulli cantilever beam model is established using beam module. The equivalent parameters of composite materials are calculated by formula for material properties. The piezoelectric double crystal cantilever models of 2D and 3D are built by solid mechanics module. The piezoelectric material parameters of PZT-5H are from COMSOL and material parameters of brass are user-defined. For geometric model, fixed constraint boundary conditions are added to the fixed end of double crystal cantilever beam, and the rest are free. For mesh generation, line segment mesh is used in 1D model, triangle mesh is used in 2D model and tetrahedron mesh is used in 3D model. Routine mesh size is chosen for 1D and 2D models. Ultra-coarsening mesh size is chosen for 3D model to save the calculation time.

The material of piezoelectric bimorph cantilever is composed of the upper and lower piezoelectric layers of PZT-5H and the intermediate substrate layer of brass.

Brass substrate is a linear elastic material defined in 1D. The material parameters in Table II are suitable for distributed parameter modeling and 1D, 2D and 3D finite element modeling.

For the definition of PZT-5H piezoelectric material parameters, the parameters of distributed parameter model are simplified in the form of 1D by Euler-Bernoulli beam theory as shown in Table III.

However, in the piezoelectric material library of COMSOL, PZT-5H is defined with 3D material parameters. For the stress-charge constitutive relation, the elastic matrix of PZT-5H is: $127.2080.2184.6780.21127.2084.6784.6784.67117.4400000000000000000022.9900022.9900023.47$ GPa. The electromechanical coupling matrix of PZT-5H is: $000000−6.623−6.62323.240017.035017.03500000$ C/m². The relative permittivity matrix of PZT-5H is: $1704.40001704.40001433.6$⁠.

Definition of linear elastic materials in 1D is enough for eigenfrequency simulation in COMSOL. For the selection of physical field, 1D FEA model chooses beam module to build homogeneous Euler-Bernoulli cantilever beam. Its equivalent material properties are calculated by formulas (3-5). Solid mechanics modules are selected to build solid models for 2D and 3D FEA models.

When the piezoelectric parameters of PZT-5H material are defined in 3D as material library of COMSOL, the eigenfrequency simulation results of COMSOL, DPM and experiment are compared in Table IV.

The first order eigenfrequency simulation results in DPM and COMSOL 1D model have almost no error compared with the experimental results. The simulation results obtained by COMSOL 2D and COMSOL 3D models with PZT-5H parameters in 3D are higher than the experimental result.

In order to explore the reason why the simulation results in COMSOL 2D and COMSOL 3D models with PZT-5H parameters in 3D are too high, the piezoelectric material parameters of PZT-5H are modified and replaced with the PZT-5H parameters in 1D as shown in Table III, which are consistent in the DPM model. The exchange formula for the 1D piezoelectric material parameters from 3D piezoelectric material parameters are:²³ 

In 1D COMSOL piezoelectric material parameters setting, the elastic matrix and the relative permittivity matrix are all set to be isotropic. All the parameters in electromechanical coupling matrix is set to be 0 except for e₃₁: $000000−16.600000000000$ C/m².

The simulation results of the first order eigenfrequency are shown in Table V.

When the piezoelectric material parameters of PZT-5H is simplified in 1D, the simulation results in COMSOL 2D and 3D models are both consistent with those in the DPM model and experiment. So COMSOL models with PZT-5H defined in 1D are more credible. In addition, the error is also related to the meshing accuracy. The denser the meshing, the higher the simulation accuracy and the longer the computational time.

The electromechanical coupling simulation of piezoelectric bimorph cantilever using COMSOL requires three physical fields: solid mechanics, electrostatic and circuit fields. Therefore, 1D beam model cannot handle electromechanical coupling analysis, and only 2D and 3D FEA models are promised to do electromechanical coupling simulation.

In the physical field of solid mechanics, an inertia load is applied to the whole structure, with the acceleration of 1 g (9.8 m/s²). Isotropic loss factors are chosen for the damping of linear elastic materials on the intermediate substrate and mechanical damping of piezoelectric materials on the upper and lower piezoelectric layers. The polarization direction of piezoelectric materials in the upper and lower layers of bimorph in series type is opposite. The method is to rotate the relative coordinate system of piezoelectric materials in Z direction by 180⁰. In the electrostatic physical field, attention should be paid to define the top electrode of the upper piezoelectric layer as terminal 1, the connection between the bottom electrode of the upper piezoelectric layer and the top electrode of the lower piezoelectric layer as terminal 2. The type of the terminal should be chosen as circuit type, and the bottom electrode of the lower piezoelectric layer as grounding. In the circuit physical field, attention should be paid to select the potential of load resistance as the voltage of terminal 1. In the multi-physical field interface, a piezoelectric coupling interface is applied to the piezoelectric bimorph.

For bimorph in parallel type, the polarization directions of piezoelectric materials in the upper and lower layers are the same. In the electrostatic physical field, the top electrode of the upper piezoelectric layer and the bottom electrode of the lower piezoelectric layer are defined as terminal 1, and the bottom electrode of the upper piezoelectric layer and the top electrode of the lower piezoelectric layer are defined as grounding.

Voltage frequency response is simulated in frequency domain. Load resistance of 470 Ω is selected for short circuit condition and 1 MΩ is selected for open circuit condition. The excitation frequency is near the first order eigenfrequency and changed with the step length of 0.1 Hz. The voltage-frequency response curves of piezoelectric bimorph in series type built by 2D and 3D COMSOL models and DPM model are simulated and shown in figure 3.

The peak point of the piezoelectric voltage frequency response curve represents the peak voltage at resonant frequency. The simulation curves of 2D and 3D COMSOL models almost coincide. The short-circuit and open-circuit resonant frequencies and corresponding voltages are analyzed and compared with DPM and experimental values, as shown in Table VI.

The 2D &3D COMSOL simulation results of resonance frequency and peak voltage are little higher than those of the test, but the maximum relative error is less than 1.21%. It can be considered that COMSOL FEA electromechanical coupling simulation can well describe the voltage-frequency response characteristics of the system.

If the 2D &3D COMSOL models use the piezoelectric material parameters of PZT-5H in 3D from material library of COMSOL, the simulation results and relative error comparison with experimental results are shown in Table VII.

The largest relative error of peak voltage is calculated in 2D COMSOL model with the value of +28.33%, and the largest relative error of resonant frequency is obtained in 3D model with the value of +4.35% by using piezoelectric material parameters of PZT-5H in 3D. The eigenfrequency simulation and electromechanical coupling simulation both indicate that the 1D piezoelectric material parameters consistent in the DPM are suitable in COMOSL modelling to achieve the desired simulation accuracy.

1D piezoelectric material parameters of PZT-5H are applied to 2D and 3D COMSOL models of piezoelectric bimorph in parallel type. The voltage-frequency response curves of piezoelectric bimorph in parallel type built by 2D and 3D COMSOL models as well as DPM model are simulated and shown in figure 4.

The simulation curves of 2D and 3D COMSOL models in parallel type are almost coincident. The simulation time of 2D COMSOL model is shorter than that of 3D COMSOL model. Therefore, priority should be given to 2D modelling and simulation in COMSOL.

The short-circuit and open-circuit resonant frequencies and corresponding voltages are analyzed. Since Alper Erturk did not give the experimental values of bimorph in parallel type, the COMSOL simulation results and relative error are compared with the DPM values, as shown in Table VIII.

The maximum relative errors of resonant frequency and peak voltage based on COMSOL simulation results are both less than 1.24% compared with the DPM values. It can be concluded that COMSOL FEA electromechanical coupling simulation can well describe the voltage-frequency response characteristics of the piezoelectric bimorph both in series type and parallel type.

In addition, the resonant frequency of bimorph is the same whether in series type or parallel type. The open circuit peak voltage in series type is almost twice as much as that in parallel type, which is consistent with the theoretical results.

The output power and optimal load resistance of 2D COMSOL model in series type and parallel type are analysed separately. In frequency domain analysis, the excitation frequency range is varying from 480 Hz to 550 Hz with the step length of 0.5 Hz. The excitation acceleration is still 1 g. The load resistance parameters are scanned by exponential step. The curves of output electric power, excitation frequency and load resistance are observed, as shown in figure 5 and figure 6. The average output power is P = V²/2R_l, where R_l is the load resistance and V is the voltage of load resistance.

There are two peak power for bimorph in both series and parallel type from the simulation curves of 2D COMSOL model. The comparison results are shown in Table IX.

The peak power of bimorph in series type and parallel type are the same, and the corresponding resonant frequency are the same, but the optimal load resistance in series type is four times than that in parallel type, which is consistent with the theory. The above situations are also suitable for sub peak. Interestingly, the resonant frequencies of main peaks are close to the open-circuit resonant frequencies, and the resonant frequencies of sub peaks are close to the short-circuit resonant frequencies. Both can achieve peak power output due to impedance matching.

Although Erturk A and Kim have given the DPM model of piezoelectric bimorph cantilever with tip mass, it is very difficult for the theoretical model to consider the vibration mode and mode electromechanical coupling coefficient under the influence of mass location, shape and moment of inertia. However, the advantage of FEA lies in the modeling and simulation of complex topological structures.

As shown in figure 7, this study establishes 2D COMSOL models of three kinds of piezoelectric bimorph cantilever with tip mass. The structural damping is changed to 0.00845 and the loss factor is 0.0169. Structure A is identical to DPM model of Erturk A. The material of mass block is brass and is located at the tip of the cantilever. Structure B moves the mass block to the center of the cantilever. Structural C replaces the mass block material with silicon, and the bonding length of the mass block occupies half of the cantilever length. All the weights of tip mass in the three structures are 0.239 g. The first-order eigenfrequency, resonant frequency and peak voltage of the three structures are simulated in open-circuit state, as shown in Table X.

The simulation results of structure A in COSMOL are consistent with those of DPM model by Erturk A, and the relative error is small, which shows the reliability of FEA modeling and simulation of piezoelectric bimorph cantilever with tip mass. Comparing the simulation results of the first-order eigenfrequency, resonant frequency and peak voltage of the three structures in open-circuit state, the difference between the structure B and structure A shows that the position of mass block has the great influence on the vibration and power generation performance of the system. The difference between the structure C and structure A shows that the bonding length of the mass block has the great influence on the vibration and power generation performance of the system. In addition, the clamping style of cantilever will also affect the performance of PVEH, and the specific analysis can be seen in the literature.²⁴ The above analyses prove the advantages in modeling and simulation of complex topological structures by COMSOL FEA.

COMSOL modelling and simulation method of bimorph cantilever for piezoelectric vibration energy harvesting is put forward in this study. 1D, 2D and 3D geometry models of piezoelectric bimorph cantilever in series and parallel types are built in COMSOL Multiphysics software. The structural eigenfrequency is simulated by linear elastic model. The resonant frequency and corresponding voltage in open-circuit and short-circuit states are simulated by electromechanical coupling simulation both in 2D and 3D model. The COMSOL results are compared with the experiment data and DPM model simulation results. A conclusion can be drawn that the 1D piezoelectric material parameters of PZT-5H, as same in DPM model, instead of the 3D piezoelectric material parameters of COMSOL, is more suitable for 1D, 2D and 3D COMSOL modelling, due to more accuracy of simulation results.

The simulation results of 2D model and 3D model are almost the same, Therefore, priority should be given to 2D modelling and simulation in COMSOL for simulation time saving.

The comparative analyses of bimorph in series and parallel connection show that the resonant frequencies of the two types are the same, and the open-circuit voltage in series type is almost twice than that in parallel type.

The simulation results of electric power, excitation frequency and load resistance show that there are two peak powers of bimorph in both series and parallel type. The peak output powers in series and in parallel are the same, and the corresponding excitation frequency is the same. The optimal load resistance in series type is four times than that in parallel type.

Performance comparation of three kinds of piezoelectric bimorph cantilever with tip mass indicates the position of the mass block and the bonding length of the mass block have a great influence on the vibration and power generation performance of the system. COMSOL FEA proves the advantages in modeling and simulation of complex topological structures.

FEA modelling and simulation methods of piezoelectric bimorph cantilever in COMSOL proposed in this study are accurate and simple. It has practical reference value for electromechanical coupling analysis and piezoelectric vibration energy harvesting design.

This work was supported in part by the National Key Research & Development (R&D) Program of China (Grant No. 2016YFB1200100), the National Natural Science Foundation of China (Grant Nos. 51890884, 51875449, 51421004, 91748207), the Key Research and Development Project of Shaanxi Province (2018ZDXM-GM-103, 2018ZDCXL-GY-02-02), the Fundamental Research Funds for the Central Universities (Grant No. xjj2017165), and the 111 Program (Grant No. B12016).