# Finite Element Modeling of a Piezoelectric Composite Beam and Comparative Performance Study of Piezoelectric Materials for Voltage Generation

Authors: Action Nechibvute; Albert Chawanda; Pearson Luhanga
Journal: ISRN Materials Science
 (2012)

---

## Abstract

A comparative study of the traditional PZT ceramics and new single crystals is critical in selecting the best material and optimization of transducer design for applications such as conversion of ambient vibrations into useful electrical energy. However, due to material and fabrication costs and the need for rapid prototyping while optimizing transducer design, primary comparisons can be based on simulation. In this paper, the COMSOL Multiphysics finite element package was used to study the direct piezoelectric effect when an external load is applied at the free end of a piezoelectric composite beam. The primary output parameters such as electric potential and electric field were studied as a function of the input strain and stress. The modeling is presented for the relatively new single crystal lead magnesium niobate-lead titanate (PMN32) and three different lead zirconate titanate ceramics (PZT-5A, PZT-5H, and PZT-4). Material performance was assessed by using a common geometry and identical excitation conditions for the different piezoelectric materials. For each material, there are three analyses performed, namely, static, eigenfrequency, and transient/time-dependent analysis. Comparative results clearly suggest that the new crystal material PMN32 is capable of outperforming presently useing piezoelectric ceramics for voltage generation.


---

## Body

## 1. Introduction

Piezoelectric materials have the novel ability of transferring from electrical to mechanical energy and vice-versa. This property is observable in many crystalline materials such as lead zirconate (PZT) ceramics where the phenomenon has found practical use in sensors and actuators [1–3]. Recently, the direct piezoelectric effect has been applied in energy harvesting where mechanical deformation on the piezoelectric material caused by ambient vibrations is converted to useful electrical energy [4–6]. The electrical energy is used to power ultralow power electronics such as wireless sensor nodes and implantable biomedical devices [5–7]. The challenge facing piezoelectric energy harvesting is the low power output of the energy generators. One way of improving the direct piezoelectric effect and subsequently the energy harvesting capabilities of piezoelectric generators is the development of single crystal materials with high voltage generation abilities under low mechanical excitation compared to traditional PZT ceramics [8]. Understanding the direct piezoelectric performance of the PZT ceramic compared to new single crystal material such as PMN32 is very critical in the selection of the best material for a particular application. A detailed comparative study of different material performances is primarily inhibited by high material and fabrication costs and hence simulations in virtual design environments are required [9–11]. Finite element modeling (FEM) is an enabling tool that can allow detailed analysis of models and can predict the behaviour of electromechanical structures under real world conditions while enabling faster and cheaper prototype development [12, 13]. In this study, we model the direct piezoelectric effect on a composite beam consisting of a stainless steel substrate sandwiched between two piezoelectric layers using the commercial FEM software package COMSOL Multiphysics version 3.5, MEMS modules. With the geometry of the composite beam fixed, four different piezoelectric materials are used, namely, PZT-5A, PZT-5H, PZT-4 and PMN32. In each case, three analyses were used: static, eigenfrequency, and transient/time-dependent. The static analysis was used to find the magnitudes and locations of maximum stress, strain, and electrical potential of the composite cantilever beam when a static load was applied to the beam’s free end. The eigenfrequency analysis was then performed to find the first six modes of vibration and their associated mode shapes. Finally, time-dependent analysis was carried out to solve for the transient solution when the applied load was time-dependent with a frequency matching and off-the-beam’s first resonance frequency.The paper is organized as follows. Section2 presents the modeling of the composite beam, highlighting the equations governing the operation of the piezoelectric materials, the geometry, and procedures used in the COMSOL Multiphysics modeling environment. Section 3 gives a summary of the simulation results. Section 4 presents a brief discussion of the results and comparison of performance of PZT and PMN32 materials. Section 5 concludes the paper.


## 2. Modeling

### 2.1. Electromechanical Model of a Linear Piezoelectric Material and FEM

The stress charge form of the electromechanical constitutive equations for linear piezoelectricity are given by [14] as follows:
(1)T=cES-eE,D=eTS+εSE,
where T is the stress vector; D is the electric flux density vector; S is the strain vector; E is the electric field vector, cE is the elasticity matrix (evaluated at constant electric field); e is the piezoelectric stress matrix; εS is the dielectric matrix (evaluated at constant mechanical strain).The equations in (1) represent the material behaviour for which the FEM software solves. The finite element discretization is performed by establishing nodal solution variables and the element shape functions over an element domain which approximates the following solution [15]:
(2)uc=NuT·u,Vc=NVT·V,
where uc is the displacement within element domain in the x, y, z directions, Vc is the electrical potential within element domain, Nu is the matrix of displacement shape functions, NV is the vector of the electrical potential shape function, u is the vector of nodal displacements, and V is the vector of nodal electrical potential. Using (2), the strain S and electric field E are thus related to the displacements and potentials by (3) and (4), respectively. Consider
(3)S=Bu·u,(4)E=-BV·V,
where
(5)Bu=[∂∂x00∂∂y0∂∂z0∂∂y0∂∂x∂∂z000∂∂z0∂∂y∂∂x]T,BV=[∂∂x∂∂y∂∂z]T.
Upon the application of the variational principle and the finite element discretization, the coupled finite element matrix equation is
(6)[M000][u¨V¨]+[C000][u˙V˙]+[KKzKzTKd][uV]=[FL],
where M is the structural mass given by M=∫ρNuNuTdv; K is the structural stiffness given by K=∫BuTcBudv; Kz is the piezoelectric coupling matrix given by Kz=-∫BuTeBVdv; Kd is the dielectric conductivity given by K=∫BVTεBVdv; C is the structural damping matrix; ρ is the mass density; F is the structural load vector (a vector of nodal forces, surface forces, and body forces); L is the electrical load vector (a vector of nodal, surface, and body charges). (The integration is over the whole element.)


### 2.2. Geometry and Modeling Definition in COMSOL Multiphysics

The composite beam and its associated dimensions are shown in Figure1 and Table 1, respectively. Throughout the modeling experiments in COMSOL Multiphysics-MEMS modules, the geometry and the dimensions shown in Figure 1 were not changed; only different piezoelectric materials were changed one type at a time. The material properties of the different piezoelectric materials are shown in Appendix A.Table 1
Dimensions of the composite beam.
Layer
Length (mm)
Width (mm)
Thickness (mm)
Top piezoelectric
25
7.5
0.1
Middle stainless steel
25
7.5
0.05
Bottom piezoelectric
15
7.5
0.1Figure 1
Cantilever composite beam.


### 2.3. Model Procedures

To simulate the geometry in Figure1, the piezo plane stress and plane stress were chosen as the Multiphysics problem in 2D. Appendix A shows the properties of different materials retrieved from the COMSOL library [16] and used in the models. The geometric model is implemented for a case where there is perfect bonding existing at the interfaces, and hence no slip occurs between the layers that make up the composite beam. The material properties from the COMSOL library were considered linear for the range of voltages used. However, in actual materials there are nonlinear and hysteresis effects and this may accept the results in this study. Nevertheless, the nonlinear effects are not prominent in the voltage range investigated in this study and hence use of the linear properties was considered very valid. It is apparent that three-dimensional (3D) modeling might be used; however, in this study 2D modeling was used. The small thickness of the piezoelectric layers and the substrate layer in the composite beam imply that in 3D, the unavoidable number of mesh elements caused the solver to consume a longer time to obtain a solution and in some cases failed to converge in the transient analysis [17, 18]. In fact, when we tried to run the 3D models for the time-dependent analysis, the solver took too long to converge to a solution and in some cases failed to converge. In addition, in problems where there is symmetry about one axis, 2D analysis is more preferable compared to 3D. In principle, the use of the FEM analysis is to find an approximate solution to the problem at hand which is close to the analytical analysis (if available). In other words, while the 2D solution is not exactly the same as the 3D solution, it reasonably approximates to the 3D solution and in turn a very close approximate to the analytical solution. Critical to finding a solution that is very close to the analytical solution is the elements meshing size. It is for this reason that the results of 2D analysis would be more accurate than the results of 3D analysis, in which one can use a very fine meshing size, whereas one can only use a coarse mesh in the 3D analysis to avoid having a tremendous number of elements that the FEM software could not solve. In the static analysis, the beam’s free end is subjected to an arbitrarily chosen external load of load FX=0 and FY=-105 N/m2. Inertial and damping effects are ignored. In the eigenfrequency analysis, the external load was removed and the model was submitted for analysis enabling the first six modes of vibration and their associated frequencies to be realized. Time-dependent analysis was performed to find the transient response to a harmonic load with the same amplitude as the static load (i.e., a harmonic load FX=0 and FY=-105sin(2πft)). The excitation frequency (f) was set to the first natural frequency found in the eigenfrequency analysis. In order to find the off-resonance performance of the composite beam, the analysis was also performed with the excitation frequencies below and above the first resonance values. Damping is very important in the analysis of time-dependent in order to get results that are as close to the real world as possible. The transient analysis model in COMSOL Multiphysics uses the Rayleigh damping given by C=αdMM+βdKK, where C is the damping matrix, M is the mass matrix, and K is the stiffness matrix [16]. Throughout the transient analysis modeling, the damping coefficients αdM and βdK were set at 94.25 and 0.0001, respectively.The piezoelectric beam was meshed using the COMSOL Multiphysics standard meshing tool at 11160 triangular elements and 65118 degrees of freedom. The optimal mesh density was determined by gradually increasing the mesh density starting with a coarse mesh and finding the corresponding final results. By following this method, it was found that as the mesh density increases the final results change as well by a considerable value. At the mesh density of 11160 elements, it was found that any increasing in the mesh density had very small impact on the final results (a mesh density of about 65000 elements changes the final results by ~0.01%).


## 2.1. Electromechanical Model of a Linear Piezoelectric Material and FEM

The stress charge form of the electromechanical constitutive equations for linear piezoelectricity are given by [14] as follows:
(1)T=cES-eE,D=eTS+εSE,
where T is the stress vector; D is the electric flux density vector; S is the strain vector; E is the electric field vector, cE is the elasticity matrix (evaluated at constant electric field); e is the piezoelectric stress matrix; εS is the dielectric matrix (evaluated at constant mechanical strain).The equations in (1) represent the material behaviour for which the FEM software solves. The finite element discretization is performed by establishing nodal solution variables and the element shape functions over an element domain which approximates the following solution [15]:
(2)uc=NuT·u,Vc=NVT·V,
where uc is the displacement within element domain in the x, y, z directions, Vc is the electrical potential within element domain, Nu is the matrix of displacement shape functions, NV is the vector of the electrical potential shape function, u is the vector of nodal displacements, and V is the vector of nodal electrical potential. Using (2), the strain S and electric field E are thus related to the displacements and potentials by (3) and (4), respectively. Consider
(3)S=Bu·u,(4)E=-BV·V,
where
(5)Bu=[∂∂x00∂∂y0∂∂z0∂∂y0∂∂x∂∂z000∂∂z0∂∂y∂∂x]T,BV=[∂∂x∂∂y∂∂z]T.
Upon the application of the variational principle and the finite element discretization, the coupled finite element matrix equation is
(6)[M000][u¨V¨]+[C000][u˙V˙]+[KKzKzTKd][uV]=[FL],
where M is the structural mass given by M=∫ρNuNuTdv; K is the structural stiffness given by K=∫BuTcBudv; Kz is the piezoelectric coupling matrix given by Kz=-∫BuTeBVdv; Kd is the dielectric conductivity given by K=∫BVTεBVdv; C is the structural damping matrix; ρ is the mass density; F is the structural load vector (a vector of nodal forces, surface forces, and body forces); L is the electrical load vector (a vector of nodal, surface, and body charges). (The integration is over the whole element.)


## 2.2. Geometry and Modeling Definition in COMSOL Multiphysics

The composite beam and its associated dimensions are shown in Figure1 and Table 1, respectively. Throughout the modeling experiments in COMSOL Multiphysics-MEMS modules, the geometry and the dimensions shown in Figure 1 were not changed; only different piezoelectric materials were changed one type at a time. The material properties of the different piezoelectric materials are shown in Appendix A.Table 1
Dimensions of the composite beam.
Layer
Length (mm)
Width (mm)
Thickness (mm)
Top piezoelectric
25
7.5
0.1
Middle stainless steel
25
7.5
0.05
Bottom piezoelectric
15
7.5
0.1Figure 1
Cantilever composite beam.


## 2.3. Model Procedures

To simulate the geometry in Figure1, the piezo plane stress and plane stress were chosen as the Multiphysics problem in 2D. Appendix A shows the properties of different materials retrieved from the COMSOL library [16] and used in the models. The geometric model is implemented for a case where there is perfect bonding existing at the interfaces, and hence no slip occurs between the layers that make up the composite beam. The material properties from the COMSOL library were considered linear for the range of voltages used. However, in actual materials there are nonlinear and hysteresis effects and this may accept the results in this study. Nevertheless, the nonlinear effects are not prominent in the voltage range investigated in this study and hence use of the linear properties was considered very valid. It is apparent that three-dimensional (3D) modeling might be used; however, in this study 2D modeling was used. The small thickness of the piezoelectric layers and the substrate layer in the composite beam imply that in 3D, the unavoidable number of mesh elements caused the solver to consume a longer time to obtain a solution and in some cases failed to converge in the transient analysis [17, 18]. In fact, when we tried to run the 3D models for the time-dependent analysis, the solver took too long to converge to a solution and in some cases failed to converge. In addition, in problems where there is symmetry about one axis, 2D analysis is more preferable compared to 3D. In principle, the use of the FEM analysis is to find an approximate solution to the problem at hand which is close to the analytical analysis (if available). In other words, while the 2D solution is not exactly the same as the 3D solution, it reasonably approximates to the 3D solution and in turn a very close approximate to the analytical solution. Critical to finding a solution that is very close to the analytical solution is the elements meshing size. It is for this reason that the results of 2D analysis would be more accurate than the results of 3D analysis, in which one can use a very fine meshing size, whereas one can only use a coarse mesh in the 3D analysis to avoid having a tremendous number of elements that the FEM software could not solve. In the static analysis, the beam’s free end is subjected to an arbitrarily chosen external load of load FX=0 and FY=-105 N/m2. Inertial and damping effects are ignored. In the eigenfrequency analysis, the external load was removed and the model was submitted for analysis enabling the first six modes of vibration and their associated frequencies to be realized. Time-dependent analysis was performed to find the transient response to a harmonic load with the same amplitude as the static load (i.e., a harmonic load FX=0 and FY=-105sin(2πft)). The excitation frequency (f) was set to the first natural frequency found in the eigenfrequency analysis. In order to find the off-resonance performance of the composite beam, the analysis was also performed with the excitation frequencies below and above the first resonance values. Damping is very important in the analysis of time-dependent in order to get results that are as close to the real world as possible. The transient analysis model in COMSOL Multiphysics uses the Rayleigh damping given by C=αdMM+βdKK, where C is the damping matrix, M is the mass matrix, and K is the stiffness matrix [16]. Throughout the transient analysis modeling, the damping coefficients αdM and βdK were set at 94.25 and 0.0001, respectively.The piezoelectric beam was meshed using the COMSOL Multiphysics standard meshing tool at 11160 triangular elements and 65118 degrees of freedom. The optimal mesh density was determined by gradually increasing the mesh density starting with a coarse mesh and finding the corresponding final results. By following this method, it was found that as the mesh density increases the final results change as well by a considerable value. At the mesh density of 11160 elements, it was found that any increasing in the mesh density had very small impact on the final results (a mesh density of about 65000 elements changes the final results by ~0.01%).


## 3. Results Summary

See Figures2, 3, 4, 5, 6, and 7 and Tables 2 and 3Table 2
Summary of static analysis results.
Displacement (mm)
Von Mises stress (MPa)
Strain ×103
Electric potential (V)
Electric field norm (MV/m)
Electric energy density (J/m3)
PZT-5A
0.6434
38.33
0.5455
15.8335
0.6627
1508.2
PZT-5H
0.625
39.96
0.5455
12.702
0.5225
1748.524
PZT-4
0.5
37.85
0.4129
14.9615
0.6146
1098.351
PMN32
1.59
38.26
1.382
18.0223
0.5399
3785.961Table 3
Summary of eigenfrequency analysis results.
f1/Hz
f2/Hz
f3/Hz
f4/Hz
f5/Hz
f6/Hz
PZT-5A
237.503744
1117.15114
2991.352491
5945.001334
9627.396607
14705.579508
PZT-5H
243.85559
1145.492784
3067.62181
6096.261632
9873.353056
15082.978781
PZT-4
274.079525
1263.173713
3391.705788
6723.378674
10903.812945
16655.978892
PMN32
146.877937
728.170014
1952.253729
3878.945263
6304.482396
9575.683795The first six deformation modes of the composite piezoelectric, using PZT-5A.
(a)
(b)
(c)
(d)
(e)
(f)Plots of von Mises stress and displacement as a function of excitation frequency.
(a)
(b)
(c)
(d)Plots of von Mises stress and electric field norm as a function of excitation frequency.
(a)
(b)
(c)
(d)Plots of electric potential and strain as function of excitation frequency.
(a)
(b)
(c)
(d)Figure 6
Displacement, strain and electric potential profiles for the different material models (at first resonance frequency).Time-dependent electric potential profiles for the different material models (at first resonance).
(a)
(b)
(c)
(d)


## 4. Results Analysis, Discussion, and Comparison of Performance

The consistency in the values of the applied load, composite beam geometry, and boundary conditions used in the FEM study laid a strong basis to allow fair comparison of results from the static and time-dependent analyses for the four types of piezoelectric material models. The results of static analysis in Table1 show that with a static load of 100 kN/m2, the model which consists of PMN32 material gave the highest electric potential of 18.0223 V. The model with PZT-5A gave the second highest electric potential of 15.8335 V, and that consisting of PZT-4 was third with 14.9615 V whereas PZT-5H model ranked fourth with a potential of 12.702 V. As may be expected from the direct piezoelectric effect, the material model with the highest electric potential is also the one with the highest electric energy density, von Mises stress, first principal strain, and beam displacement. In addition, the static analysis confirmed that strain and electric potential are maximum at the beam’s fixed end and are reduced in value as one traverses towards the free end of the beam (Appendix B). Generally, the static analysis suggests that for the studied geometry, PMN32 exhibits strongest direct piezoelectric performance relative to PZT ceramics in terms of voltage generation.The results of eigenfrequency analysis presented in Table2 were used to set up the excitation frequency of the time-dependent load in the transient analysis. The eigenfrequency analysis reveals that the PZT-4 material model has the highest first resonance frequency of 274.08 Hz, followed by PZT-5H with 243.86 Hz and PZT-5A with 237.50 Hz, and PMN32 material model has the lowest first resonance frequency of 146.88 Hz. Guided by the values of the first resonance frequencies, time-dependent analysis using the same magnitude for the dynamic load was performed for the different material models and the plots in Figures 3–5 show the strong frequency dependence of the key electrical and mechanical parameters in the direct piezoelectric property of the materials. Time-dependent results show that stress, strain, electric field, and electrical potential are at their maximum values at the first resonance; the values of these quantities decrease as the excitation frequency deviates from the resonance frequency. The model which consists of PM32 material gave the highest electric potential of 111.7425 V, followed by PZT-5A with 92.982 V, and that consisting of PZT-4 was third highest with 82.51, and PZT-5H model gave the lowest potential of 74.2325 V. Figure 7 shows the electric potential distribution along the composite beam’s top surface. From the figure, it is observed that the electric potential is the highest at the beam’s fixed end and changes from maximum to minimum with respect to time. This is in excellent agreement with expectation from the direct piezoelectricity relationship: the maximum electric potential exists at the area that has maximum strain. The relative performance comparison is presented in Figure 8. These results confirm the superior performance of PMN32 over PZT ceramics as has been reported by the work of other researchers [10, 19, 20].Relative performance comparison of the four piezoelectric material models.
(a)
(b)


## 5. Conclusion

In this paper, FEM using the commercial COMSOL Multiphysics package was used to study performance of a composite piezoelectric beam. For both static and time-dependent mechanical excitations, simulation results have shown that PM32 has the highest greatest voltage generation ability. PZT-5A was a close second, followed by PZT-4 in the third rank while PZT-5H ranked forth with the lowest electric potential. In addition to the largest values of the electrical potential, the PMN32 material model also gave the highest average values of electric field, electric energy density, stress, strain, and displacement. On the other hand, PMN32 material model has the lowest first resonance frequency of compared to the PZT models. For future research, the strong piezoelectric effect demonstrated by PMN32 needs to be further studied, particularly in the context of vibration to electric power generation to enable optimization of the material properties and design geometry of the cantilever device.



---

*Source: 921361-2012-09-18.xml*