An effective hybrid computational framework is described here in order to assess the nonlinear dynamic response of piezoelectric energy harvesting devices. The proposed strategy basically consists of two steps. First, fully coupled multiphysics finite element (FE) analyses are performed to evaluate the nonlinear static response of the device. An enhanced reduced-order model is then derived, where the global dynamic response is formulated in the state-space using lumped coefficients enriched with the information derived from the FE simulations. The electromechanical response of piezoelectric beams under forced vibrations is studied by means of the proposed approach, which is also validated by comparing numerical predictions with some experimental results. Such numerical and experimental investigations have been carried out with the main aim of studying the influence of material and geometrical parameters on the global nonlinear response. The advantage of the presented approach is that the overall computational and experimental efforts are significantly reduced while preserving a satisfactory accuracy in the assessment of the global behavior.

1. Introduction

The possibility of using electronic devices that do not require periodic replacement of batteries is nowadays the major challenge in several engineering fields. In particular, vibration-based piezoelectric energy harvesting devices are emerging as a valid technological option to power miniaturized electronic sensors in civil structural health monitoring applications [13]. It can be easily understood, therefore, that the implementation of an efficient computational framework for the analysis and design of these devices is of paramount importance in order to foster the future large-scale applications of this technology.

In this perspective, the electromechanical response of piezoelectric devices can be assessed through different numerical and analytical techniques, for example, reduced-order models, finite element (FE) models, and circuit analogy methods [2, 49]. Most of the efforts in this field are based on the hypothesis of linear behavior [1014], but nonlinear phenomena can greatly impact the final performances [1518]. The causes of the nonlinear response can be traced back to several mechanisms, such as instability phenomena [1922], nonlinear material constitutive law [2327], geometric effects [2830], impacts [3134], and damping [35]. Additionally, nonlinearities can arise from the electrical circuit (diodes, among others) [36].

As regards the applications in civil engineering, the dynamic response of most structures and infrastructures is characterized by low frequency content. In these cases, the frequency tuning of piezoelectric beams necessitates the use of flexible materials, low thickness-to-length ratios, and/or relatively heavy additional masses [37, 38]. Such expedients make piezoelectric energy harvesting devices more prone to exhibit nonlinear mechanical behaviors.

While FE simulations can be a viable computational strategy for nonlinear static analyses, reduced-order models (ROMs) are more attractive for nonlinear dynamic analyses because they allow saving elaboration time. Moreover, ROMs are useful in order to separate and identify the effects due to material and geometric nonlinearities [35]. With some exceptions [5], however, most reduced-order models assume a linear electromechanical response and the modal superposition principle is extensively adopted to derive the state-space representation of the system.

Therefore, in this paper, we propose an efficient multiscale hybrid approach to model accurately the nonlinear dynamic response of PVDF energy harvesters. The term multiscale refers specifically to the device and system scales. First, the FE method is employed to solve the equations governing the response under static loading (device scale). Hence, a global curve that provides the tip displacement evolution for increasing values of the external load is obtained (i.e., pushover curve). The global FE-based solutions are then used in place of experimental data to identify the values of linear and nonlinear lumped coefficients of the reduced-order model (system scale). Finally, the nonlinear differential equations governing the dynamic behavior are solved to estimate the frequency response functions of tip displacement and output voltage.

2. Numerical Modeling

2.1. A Short Review on Nonlinear Electroelasticity

Piezoelectric polymers like PVDF [39] represent a valid solution for the development of flexible energy harvesting devices [2, 40]. The main advantage of PVDF with respect to other piezoelectric materials (in particular piezoceramics) is the possibility of sustaining large displacement without failure or drastic reduction of the piezoelectric efficiency [4143]. Throughout this paper, therefore, it will be assumed that the piezoelectric layers are made of PVDF.

As soon as a piezoelectric solid undergoes large deformations and rotations, the classical small strain electromechanical constitutive equations lead to incorrect results. For the sake of completeness, we briefly review hereafter the equations for the continuum mechanical description of a piezoelectric solid under large strains. The interested reader can refer to [44] for a more complete discussion.

The reference and deformed configurations are denoted by and , respectively, where . When the electromechanical body deforms, the nonlinear mapping function at time instant maps the material point onto :The displacement vector is obtained as the difference between the positions vectors of the current and initial configuration:whereas the deformation gradient can be defined as a function of the displacement gradient :whereAccording to Faraday’s law,where is the electric field vector in the current configuration. Consequently, it is possible to define as the gradient of a scalar electric potential :Velocity and acceleration of a material point with respect to the reference configuration are defined, respectively, by the following material time derivatives:In the current configuration, the balance of momentum and Gauss’s law state thatwhere represents the mechanical Cauchy stress tensor, denotes the electric displacement, is the mechanical density, is free electric charge density, and indicates the volume force (in the current configuration). The balance of mass implies that , where is the density in the initial configuration and . Local balance of angular momentum guarantees that . Moreover, if   is the electric polarization vector in the current configuration and is the vacuum permittivity, thenThe transformation between volume elements , and electric charges , in the current and reference configurations, respectively, is based on the following equations:Furthermore, the condition ensures that the tensor is not singular and, as a consequence, the deformation process will be smooth.

In the initial configuration, the local balance of momentum given by (8) can be recast with respect to different stress and strain measures:where and are the total first and second Piola–Kirchoff stress tensors, respectively. Moreover, represents the body force in the initial configuration. Constitutive equations satisfying the material objectivity principle are also required. To this end, it is assumed that a strain energy density function exists for the electromechanical body that, in general, can be defined with respect to different kinematics tensors, namely, , , and and the electric field vector . All the quantities refer to the initial configuration. Here, indicates the Green–Lagrange strain tensor that can be expressed as a function of the deformation gradient tensor bywhereis the right Cauchy–Green tensor. Objectivity requires that . A compressible Neo–Hookean type material model with a total energy density is used in this work:where and are the Lame constants, and are further material constants to be calibrated, and , , and are computed for a transversely isotropic material according to [45, 46] and they are equal to where is the identity tensor. In this study, it is assumed that the interaction between electric fields and matter is mainly confined within the finite space occupied by the matter. Once is defined, it is possible to derive the following constitutive equations:Here, is the dielectric displacement vector computed in the initial configuration while the total first Piola–Kirchoff stress tensor and the total Cauchy stress tensor are equal torespectively. The transformation from the material to the current configuration is possible by means of the following relationships:The Dirichlet and Neumann boundary conditions for the mechanical field arewhere and are prescribed mechanical displacement and surface traction vectors in the reference configuration, respectively. The boundary of the domain is ( and are its Dirichlet and Neumann portions, resp.), with and . Moreover, is the outward unit normal to . The boundary conditions for the electric field arewhere and are prescribed values of electric potential on and electric charge flux on , respectively. Moreover, , .

A cantilever-type configuration is here assumed for the energy harvesting device (see Figure 1). It is made of a piezoelectric layer and a substrate layer used for the deposition of the polymeric mixture during the fabrication process before the polarization of dipoles. This second material layer is described here by a compressible Neo–Hookean type material with total energy densities :

2.2. Finite Element Discretization for Static Problems

The standard nodal FE discretization is employed following [46, 47]. In doing so, the advanced symbolic computational tools for FE analysis available in the AceGen/AceFEM are useful to facilitate the full automation of the linearization process. Therefore, let and be the discretized displacement and electric potential fields, respectively, where are the nodal displacements and are the nodal electric potentials ( are the shape functions). According to [47], since all quantities in (14) and (21) depend on the displacement field and the electric potential, the resulting system of nonlinear equations has the general form of , where the unknown variables have to be determined ( is the total number of global unknowns of the problem). If is a subset of the global vector of unknowns on which the th element depends explicitly, then the element and the Gauss point contributions to the global residuals are explicit functions of , that is, and . In particular, at the element level, the internal residuals are obtained using AceGen as follows [47]:where is the unknown vector related to the element that collects all nodal displacements and/or nodal electrical potentials . Within the FE procedure, the global residuals are approximated aswhere is the standard FE assembly operator, indicates the number of elements, is the number of Gauss points, and indicates a generic Gauss point ( and are the Gauss point weight and Jacobian determinant, resp.). The Gauss point contribution to the residuals is .

The requirement of zero global residuals yields a set of nonlinear systems of equations, which is solved numerically by means of a Newton–Raphson algorithm. The system of equations is parameterized, introducing a load factor in the form . Thus, it is considered as [47]where denotes the contribution of internal forces to the nodal force vector and is the reference load vector associated with the pattern of the applied nodal forces and/or electrical charges. Taylor series expansion of (24) at the (known) state results in the expressionwhere the upper index denotes the quantities at the iteration and , whereas indicates the directional derivative required for the linearization. In particular, the linearization of the vector yields the tangent matrix . Within the FE procedure, the tangent operator is formed from the Gauss point tangent operator :

2.3. Finite Element Discretization for Dynamic Problems

For our applications, it is useful to separate residual and stiffness matrix terms that refer to the mechanical and electrical unknowns and . To this end, the total energy of the discretized system is introduced as a function of two contributions: the first term is and refers to the electromechanical domain whereas the second term is and refers to the mechanical domain. Thus,Consequently, the following equations are derived: Using an implicit time integration scheme and Newton’s method iterations to determine the dynamic equilibrium at time , the coupled nonlinear FE equations arewhere and are the global time-dependent nodal displacement vector and electric potential vector, respectively, whereas indicates the time step increment. The upper dots indicate the time derivative. The matrices , , and are the structural mass, damping, and stiffness matrices, respectively. The matrices and are the stiffness matrices due to piezoelectric mechanical coupling. The matrix is the stiffness matrix resulting from the electrical fields. The vectors and are the unbalanced force vectors due to mechanical and electrical contributions, respectively. According to the classical Rayleigh damping approach, a mass- and stiffness-proportional damping matrix is considered:where and are constant multipliers that can be determined by assuming a constant damping ratio .

2.4. Reduced-Order Modeling

Reduced-order models allow mitigating the computational efforts required for FE-based dynamic analyses. According to [8, 48], the forced vibrations of linear elastic piezoelectric cantilever generators can be defined using the modal coordinates as follows:for each mode , where is the modal mechanical damping ratio, is the undamped natural frequency, and are modal electromechanical coupling terms, is the modal mechanical forcing function, is the capacitance, is the load resistance, and is the voltage response across the external resistive load. Having calculated the modal coordinates , the transverse displacement of the neutral axis relative to the moving base at position and time is equal towhere is a mass normalized eigenfunction (mode shape). Different analytical expressions can be derived for each lumped coefficient depending on the considered system (series or parallel connection of the piezoelectric layers, unimorph or bimorph layout, etc.). Using Hamilton’s principle and Galerkin’s method, Stanton et al. [49] extended the formulation [8] to the case of nonlinear piezoelectricity.

In this paper, we will focus on the numerical computations of the electromechanical response for piezoelectric unimorphs that exploit a nonlinear behavior to increase their bandwidth. In particular, following [49] and assuming a linear damping, the single-mode approximation for the piezoelectric cantilever in Figure 1 under base acceleration readswhere , , and are nonlinear coefficients that can be determined based on the orthogonal basis functions used to represent in the modal space the transverse deflection of the device, whereas is the modal mechanical forcing function coefficient that multiplies the base acceleration . Based on experimental evidence [17], it has been shown that the nonlinear response mainly depends on mechanical stiffening/softening effects; therefore, we will assume that both and are equal to zero in the following computations. Moreover, for a rectangular cross section, , whereas the equivalent capacitance iswhere , , and are the width, length, and thickness of the piezoelectric layer, respectively, whereas is the permittivity constant. The th modal electromechanical coupling term for unimorph configuration can be obtained as [48]where is the piezoelectric constant and is the position of the bottom PVDF layer from the neutral axis. It is understood that in (33). Following [48], , where is the distance from the top of the PVDF layer to the neutral axis, and it is equal towith , where and are the elastic modulus of the substrate (here assumed to be made of mylar, without loss of generality) and that of the PVDF layer, respectively. Moreover, is the thickness of the substrate.

3. Hybrid Computational Strategy

The solution of nonlinear dynamic equilibrium equations using mode superposition techniques was first studied in [50] and implemented successfully in [51] for mechanical problems. A review of reduced-order model techniques is reported in [52], including a description of the modal coordinate reduction. The use of modal derivatives for nonlinear model order reduction has been recently proposed in [53] in the context of isogeometric FE analysis. To the authors’ knowledge, few studies describe efficient numerical procedures for reduced-order model techniques and coupled domains (such as in electromechanical systems). The most important contributions are provided in [54, 55], whose strategy is implemented in [56].

In this framework, a hybrid multiscale computational approach is proposed here for modeling efficiently the nonlinear dynamic response of piezoelectric cantilever devices (see Figure 2). It is based on two steps. Specifically, we first solve (24) in the physical domain of the piezoelectric cantilever (device scale) using the FE method. In doing so, a pattern of forces () is statically applied to the structural model (including nonlinear effects) and the total reaction at the base is plotted against a reference displacement (i.e., the tip displacement ). This allows obtaining the capacity curve (also named pushover curve) of the device. This curve essentially allows reducing the nonlinear static problem at the device scale to an equivalent single degree of freedom (SDOF) system. The second step of the analysis consists in the derivation of a phenomenological law for the nonlinear stiffness from the estimated capacity curve. Hence, the goal is to compute the coefficients of a nonlinear spring suitable for reduced-order modeling. With reference to (33), these coefficients are and , where is the equivalent mass. Looking at (33), it is clear that the effects of material and geometrical nonlinearities are merged in a single coefficient . This assumption is meaningful for flexible PVDF EHs. In fact, the material nonlinear constitutive equations, derived from the potentials introduced in (14) and (21), allow for catching large deformations and, consequently, geometrical stiffening or softening effects. Therefore, it appears reasonable to consider the two sources of nonlinearity as coupled in the reduced-order model. For each case, the coefficient is evaluated by fitting the analytical approximation with the pushover curves. The obtained values are used to update the equations set (33)-(34), which is the form of the ordinary differential equations system that describes the dynamics of the piezoelectric energy harvester. Several patterns of forces have to be considered for a multimodal response analysis. Finally, the nonlinear state-space model equation is solved using a standard time discretization algorithm. This numerical strategy is implemented in Mathematica using the advanced symbolic computational tools available in the AceGen/AceFEM [47]. Plane-stress four-node large-displacement electromechanical elements with three degrees of freedom for node are implemented in AceGen/AceFEM. Analytical expressions can be also derived for harmonic base vibrations (see the Appendix).

4. Applications

4.1. Numerical Data

Energy harvesters similar to those tested by Elvin et al. [1, 17] are considered here. As shown in Figure 1(a), the harvester consists of a PVDF film with electrodes on both sides connected with wires to output the generated charge. Moreover, there are coatings to protect the device from damage. The thickness of the PVDF layer is 28 m, and the piezoelectric strain constant is  pC/N. According to Figure 1(a), five layers are considered in the real configuration of the device, namely, (i) a mylar substrate with a thickness equal to 6 m, (ii) a silver layer with a thickness equal to 8 m, (iii) a central layer of PVDF with a thickness equal to 28 m, (iv) another silver layer with a thickness equal to 8 m, and (v) a final layer of mylar with a thickness equal to 140 m. Density and Young’s modulus of PVDF and mylar are provided in Table 1, together with the other main material and geometrical data adopted in this numerical study.

4.2. Thickness and Stiffness Homogenization

According to [12], we use the multimorph structure concept to physically and mathematically describe the equivalent stiffness of the piezoelectric energy harvester. A fixed-free cantilever configuration is assumed. Moreover, each multiple thin material film layer has length , width , and thickness . It is assumed that the thin-film layer interfaces are smooth and continuous and do not slip with respect to each other. Each layer is considered uniform with Young’s modulus , rotational inertia , and cross-sectional areas = . The subscript denotes the th layer. The first four modal shapes and modal frequencies of this device have been calculated (see Figure 3). The homogenized flexural rigidity about the neutral axis located at is then given by [12]whereand is the location of the axis of the th layer with respect to an arbitrary reference. is the number of layers considered (5 for our device). The homogenized stiffness is given in Table 1. Since a lumped tip mass is assumed, it is observed that modal truncation to the first mode is enough to describe the overall device response. For broadband energy harvesting and/or for piezoelectric beams with different boundary constraints, however, several modes can be required. Moreover, it is not necessary in the present work to use FE discretization for each material layer. In fact, an equivalent mylar layer thickness is determined in such a way that the first resonance frequency predicted using the FE method and that obtained by means of homogenized stiffness given by (38) are close to each other. The equivalent mylar layer thickness is given in Table 1 (see Figure 1(b)). Although each layer can be modeled explicitly in the FE model, the use of an equivalent layer is useful to reduce the total elaboration time in nonlinear analyses because of the more refined mesh required by thin layers, as well as prevent distortion phenomena of the elements at the interface between layers having different thicknesses.

4.3. Capacity Curve

According to the proposed hybrid computational strategy, a FE analysis is first performed in order to derive the capacity curve of the device. The adopted FE mesh is illustrated in Figure 4. Open-circuit conditions are assumed. Figure 5 compares the results of the large-displacement analysis (performed by means of AceFEM) with the experimental evidence provided in [17]. Herein, it can be noted that the vertical tip displacements are normalized by the beam length while the load is normalized by , where is the homogenized stiffness coefficient. This comparison demonstrates a very good agreement between numerical predictions and experimental outcomes. The computer program Comsol Multiphysics is also used to further validate the results obtained with the developed FE codes. Finally, contour levels of the stress components and deformed shape are shown in Figure 6.

4.4. Calibration of the Analytical Model

The second step of the proposed hybrid computational procedure aims at deriving a phenomenological law for the nonlinear stiffness from the estimated capacity curve. FE analyses allow for evaluating the total reaction at the device clamped end as a function of the tip displacement . This facilitates the calibration of the nonlinear spring-type element of the reduced-order model. Hence, the capacity curve is now approximated using the analytical Duffing model by means of the relationship where and characterize the constitutive law of the nonlinear spring. In the static case, . Figure 7 provides the comparison between the FE analysis and the approximation obtained by means of the analytical Duffing model.

4.5. Nonlinear Frequency Response Functions

The influence of material and geometrical parameters on the nonlinear frequency response of the device is now investigated. To this end, four different values of , , , and are considered in the nonlinear static analysis. For each case, the capacity curves are obtained and these global FE solutions are used in place of experimental curves to identify the values of the linear and nonlinear lumped parameters. The resulting nonlinear dynamic equation system in (33) and (34) is solved in the frequency domain (see the Appendix). In particular, the frequency response functions (FRFs) of tip displacement and voltage through a resistance  MΩ are determined for four values of , , , and . The reference values for , , , and are those listed in Table 1 and are indicated as , , , and , respectively, throughout this parametric study. Figure 8 provides the capacity curves obtained through the nonlinear FE analyses while Table 2 reports the computed values of the nonlinear stiffness coefficients for each case.

Figures 9 and 10 highlight the role of the piezoelectric elastic modulus whereas that of the substrate can be inferred from Figures 11 and 12. The effects of the thickness layers are illustrated in Figures 13, 14, 15, and 16. Overall, it can be observed that the occurrence of geometric stiffening effects induces a significant variation of the frequency at which the device exhibits the maximum tip displacement and output voltage. For instance, as shown in Figure 9 for the original device configuration (case a), the peak displacement is achieved at 316 rad/s for small deformations (amplitude of the base excitation equal to or less than 1, with  m/s2). For very large deformations (intensity of the base excitation equal to 10), the peak displacement is achieved at 350 rad/s. In such a case, an increment of the elastic modulus up to 60% shifts this frequency to 420 rad/s.

5. Validation

5.1. Experimental Layout and Analysis

The proposed approach is further validated against new experimental data. The tested piezoelectric energy harvester is shown in Figure 1(a) (the corresponding data are given in Table 1). The adopted experimental testing equipment (see Figure 17) consists of a signal generator (exciter) that provides a controlled input voltage to the piezoelectric structure and an analysis system with tools for signal processing and modal characterization. Modal data (i.e., natural frequencies and damping ratios) are extracted based on the SDOF curve fitting of the FRF. Experimental modal shapes are derived considering the receptance matrix :Here, indicates the imaginary unit. Using the mode superposition technique, it is possible to express in the formwhere the th element of the matrix is obtained from the measured transfer function between points and . indicates experimental modal shapes.

The frequency response analysis is performed by evaluating the out-of-plane displacements of the piezoelectric cantilevers by a noncontact measurement system. An AC voltage signal generated by a signal generator embedded into the vibrometer system has been applied to the cantilever, thus resulting in the actuation and deflection of the tip. The cantilever is tested assuming a frequency sweep at a given voltage in order to measure the peak of the tip deflection. The measured first resonance frequency and damping ratio are 324 rad/sec and 4%, respectively (see Figure 18). Moreover, two base acceleration time histories are considered for the experimental characterization of the device response under large strains. In particular, the input acceleration values are measured directly by the LDV system, by focusing the laser beam on the fixed part of the cantilever. The harmonic dynamic inputs have amplitude equal to 3.7 and 10.5 while the frequency is 330 rad/sec. The FFTs of the input acceleration are provided in Figures 19(a) and 19(b) while the corresponding FFTs (see Figures 19(c) and 19(d)) and time histories of the output voltage (see Figures 19(e) and 19(f)) are measured through a resistance  M.

5.2. Comparison between Experimental Data and Numerical Predictions

The validation of the proposed hybrid computational strategy is now performed by means of a direct comparison between the predicted and the measured output voltage in the time domain. After estimating the value of the nonlinear term through the comparison between the reduced-order model and the static solution of (24), (33) and (34) are solved using a Runge–Kutta algorithm. Finally, the numerical time histories of tip displacement and/or velocity and voltage difference on the resistance are computed. Figures 19(e) and 19(f) provide the comparison between experimental data and predicted values of the voltage difference in the time domain for an amplitude of the base acceleration equal to 3.7 and 10.5. We can observe a good agreement between experimental evidence and numerical simulations, as confirmed by the results in Table 3. Consequently (as pointed out in the Appendix), based on experimental evidence, it is thus confirmed for PVDF EHs that a harmonic form can be assumed for tip displacement and output voltage in case of harmonic base acceleration, and the corresponding amplitudes have a nonlinear dependence with respect to the frequency of the input signal.

6. Conclusion

An innovative multiscale hybrid approach has been proposed to model accurately the nonlinear dynamic response of piezoelectric cantilever beams. Once the tip displacements are obtained as a function of the load increments by means of static nonlinear FE analyses, linear and nonlinear lumped coefficients of the spring element into the analytical Duffing model are calibrated in such a way so as to approximate as best as possible the reference nonlinear capacity curve. Finally, the nonlinear dynamic differential equations governing the response of the piezoelectric cantilever beam can be solved in order to estimate the frequency response functions of tip displacement and output voltage. The attractive features of this computational procedure are twofold. From a numerical standpoint, the FE analysis is performed to solve a static problem, which is less time-consuming than a dynamic problem. Since the nonlinear behavior is reflected into the reduced-order model adopted for the dynamic analyses, through the comparison with FE-based capacity curves, devices with a larger bandwidth and better performances in the frequency range of interest can be fabricated. The experimental work was limited to the estimation of damping and natural frequencies of the device without nonlinear effects in order to predict the electrical response in the time domain. A large parametric investigation has been also performed in order to assess the role of material and geometrical parameters in the nonlinear response of piezoelectric unimorphs for energy harvesting applications. Final results have shown that the proposed hybrid approach leads to satisfactory results while reducing the overall numerical and experimental efforts. Future work will concern the further validation of the proposed numerical procedure with the aim of identifying experimentally the nonlinear coefficients and the overall FRFs.


Let us consider (33) and (34) under the hypothesis of sinusoidal external acceleration acting at the base with amplitude , frequency , and phase . Moreover, based on experimental evidence (see Figure 18), let us assume that the electrical contribution to the mechanical equation is small enough to be neglected. For the single-mode approximation, (33) and (34) reduce toThe solution of the first equation is well established [57] and it allows correlating the amplitude of the tip oscillation with the amplitude of the external acceleration as a function of and :According to [57], provided that the following inequality is verified:the tip oscillation can be well approximated by a harmonic signal in the form . It is worth stressing that is a nonlinear function of the input excitation frequency . Therefore, the second equation can be recast in the formThe steady-state solution of (A.4) is given bywhere is a complex quantity expressed as follows:Therefore, the steady-state amplitude of the output voltage can be computed as

Conflicts of Interest

The authors declare that they have no conflicts of interest.


Claudio Maruccio acknowledges the support from the Italian MIUR through the project FIRB Futuro in Ricerca 2010 “Structural Mechanics Models for Renewable Energy Applications” (RBFR107AKG). Giuseppe Quaranta acknowledges the support from Sapienza University of Rome through the project “Smart Solutions for the Assessment of Structures in Seismic Areas.”