Research Article  Open Access
A TwoStep Hybrid Approach for Modeling the Nonlinear Dynamic Response of Piezoelectric Energy Harvesters
Abstract
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 reducedorder model is then derived, where the global dynamic response is formulated in the statespace 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, vibrationbased piezoelectric energy harvesting devices are emerging as a valid technological option to power miniaturized electronic sensors in civil structural health monitoring applications [1–3]. 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 largescale applications of this technology.
In this perspective, the electromechanical response of piezoelectric devices can be assessed through different numerical and analytical techniques, for example, reducedorder models, finite element (FE) models, and circuit analogy methods [2, 4–9]. Most of the efforts in this field are based on the hypothesis of linear behavior [10–14], but nonlinear phenomena can greatly impact the final performances [15–18]. The causes of the nonlinear response can be traced back to several mechanisms, such as instability phenomena [19–22], nonlinear material constitutive law [23–27], geometric effects [28–30], impacts [31–34], 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 thicknesstolength 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, reducedorder 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 reducedorder models assume a linear electromechanical response and the modal superposition principle is extensively adopted to derive the statespace 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 FEbased solutions are then used in place of experimental data to identify the values of linear and nonlinear lumped coefficients of the reducedorder 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 [41–43]. 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 cantilevertype 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 :
(a)
(b)
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 timedependent 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 stiffnessproportional damping matrix is considered:where and are constant multipliers that can be determined by assuming a constant damping ratio .
2.4. ReducedOrder Modeling
Reducedorder models allow mitigating the computational efforts required for FEbased 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 singlemode 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 reducedorder 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 reducedorder 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 reducedorder 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 reducedorder 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 statespace 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]. Planestress fournode largedisplacement 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 fixedfree cantilever configuration is assumed. Moreover, each multiple thin material film layer has length , width , and thickness . It is assumed that the thinfilm 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 crosssectional 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.
(a)
(b)
(c)
(d)
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. Opencircuit conditions are assumed. Figure 5 compares the results of the largedisplacement 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 springtype element of the reducedorder 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.
(a)
(b)
(c)
(d)
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/s^{2}). 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
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 outofplane 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.
(a)
(b)
(c)
(d)
(e)
(f)
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 reducedorder 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 timeconsuming than a dynamic problem. Since the nonlinear behavior is reflected into the reducedorder model adopted for the dynamic analyses, through the comparison with FEbased 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.
Appendix
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 singlemode 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 steadystate solution of (A.4) is given bywhere is a complex quantity expressed as follows:Therefore, the steadystate amplitude of the output voltage can be computed as
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
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.”
References
 N. G. Elvin, N. Lajnef, and A. A. Elvin, “Feasibility of structural monitoring with vibration powered sensors,” Smart Materials and Structures, vol. 15, no. 4, article no. 011, pp. 977–986, 2006. View at: Publisher Site  Google Scholar
 C. Maruccio, G. Quaranta, L. D. Lorenzis, and G. Monti, “Energy harvesting from electrospun piezoelectric nanofibers for structural health monitoring of a cablestayed bridge,” Smart Materials and Structures, vol. 25, no. 8, Article ID 085040, 2016. View at: Publisher Site  Google Scholar
 G. Quaranta, F. Trentadue, C. Maruccio, and G. C. Marano, “Analysis of piezoelectric energy harvester under modulated and filtered white Gaussian noise,” Mechanical Systems and Signal Processing, vol. 104, pp. 134–144, 2018. View at: Publisher Site  Google Scholar
 J. G. Smits, S. I. Dalke, and T. K. Cooney, “The constituent equations of piezoelectric bimorphs,” Sensors and Actuators A: Physical, vol. 28, no. 1, pp. 41–61, 1991. View at: Publisher Site  Google Scholar
 A. H. Nayfeh, M. I. Younis, and E. M. AbdelRahman, “Reducedorder models for MEMS applications,” Nonlinear Dynamics, vol. 41, no. 13, pp. 211–236, 2005. View at: Publisher Site  Google Scholar
 M. Staworko and T. Uhl, “Modeling and simulation of piezoelectric elements  comparison of available methods and tools,” Mechanics, vol. 27, no. 4, pp. 161–171, 2008. View at: Google Scholar
 Y. Yang and L. Tang, “Equivalent circuit modeling of piezoelectric energy harvesters,” Journal of Intelligent Material Systems and Structures, vol. 20, no. 18, pp. 2223–2235, 2009. View at: Publisher Site  Google Scholar
 A. Erturk and D. J. Inman, “An experimentally validated bimorph cantilever model for piezoelectric energy harvesting from base excitations,” Smart Materials and Structures, vol. 18, no. 2, Article ID 025009, 2009. View at: Publisher Site  Google Scholar
 M. Pozzi, “Impulse excitation of piezoelectric bimorphs for energy harvesting: A dimensionless model,” Smart Materials and Structures, vol. 23, no. 4, Article ID 045044, 2014. View at: Publisher Site  Google Scholar
 S. Adhikari, M. I. Friswell, and D. J. Inman, “Piezoelectric energy harvesting from broadband random vibrations,” Smart Materials and Structures, vol. 18, no. 11, Article ID 115005, 2009. View at: Publisher Site  Google Scholar
 M. I. Younis, Elements of LumpedParameter Modeling in MEMS, Springer Science and Business Media, Berlin, Germany, 2011.
 R. Andosca, T. G. McDonald, V. Genova et al., “Experimental and theoretical studies on MEMS piezoelectric vibrational energy harvesters with mass loading,” Sensors and Actuators A: Physical, vol. 178, pp. 76–87, 2012. View at: Publisher Site  Google Scholar
 S. Zhao and A. Erturk, “Electroelastic modeling and experimental validations of piezoelectric energy harvesting from broadband random vibrations of cantilevered bimorphs,” Smart Materials and Structures, vol. 22, no. 1, Article ID 015002, 2013. View at: Publisher Site  Google Scholar
 M. R. A. Nabawy and W. J. Crowther, “Dynamic electromechanical coupling of piezoelectric bending actuators,” Micromachines, vol. 7, no. 1, pp. 1–10, 2016. View at: Publisher Site  Google Scholar
 M. Lallart and D. Guyomar, “An optimized selfpowered switching circuit for nonlinear energy harvesting with low voltage output,” Smart Materials and Structures, vol. 17, no. 3, Article ID 035030, 2008. View at: Publisher Site  Google Scholar
 J. M. Renno, M. F. Daqaq, and D. J. Inman, “On the optimal energy harvesting from a vibration source,” Journal of Sound and Vibration, vol. 320, no. 12, pp. 386–405, 2009. View at: Publisher Site  Google Scholar
 N. G. Elvin and A. A. Elvin, “Large deflection effects in flexible energy harvesters,” Journal of Intelligent Material Systems and Structures, vol. 23, no. 13, pp. 1475–1484, 2012. View at: Publisher Site  Google Scholar
 R. Patel, S. McWilliam, and A. A. Popov, “Optimization of piezoelectric cantilever energy harvesters including nonlinear effects,” Smart Materials and Structures, vol. 23, no. 8, p. 085002, 2014. View at: Publisher Site  Google Scholar
 K. Das and R. C. Batra, “Pullin and snapthrough instabilities in transient deformations of microelectromechanical systems,” Journal of Micromechanics and Microengineering, vol. 19, no. 3, p. 035008, 2009. View at: Publisher Site  Google Scholar
 A. Erturk and D. J. Inman, “Broadband piezoelectric power generation on highenergy orbits of the bistable Duffing oscillator with electromechanical coupling,” Journal of Sound and Vibration, vol. 330, no. 10, pp. 2339–2353, 2011. View at: Publisher Site  Google Scholar
 H. Vocca, I. Neri, F. Travasso, and L. Gammaitoni, “Kinetic energy harvesting with bistable oscillators,” Applied Energy, vol. 97, no. 9, pp. 771–776, 2012. View at: Publisher Site  Google Scholar
 B. Ando, S. Baglio, G. L'Episcopo, and C. Trigona, “Investigation on mechanically bistable MEMS devices for energy harvesting from vibrations,” Journal of Microelectromechanical Systems, vol. 21, no. 4, Article ID 6189362, pp. 779–790, 2012. View at: Publisher Site  Google Scholar
 S. C. Stanton, C. C. McGehee, and B. P. Mann, “Reversible hysteresis for broadband magnetopiezoelastic energy harvesting,” Applied Physics Letters, vol. 95, no. 17, Article ID 174103, 2009. View at: Publisher Site  Google Scholar
 R. Ramlan, M. J. Brennan, B. R. Mace, and I. Kovacic, “Potential benefits of a nonlinear stiffness in an energy harvesting device,” Nonlinear Dynamics, vol. 59, no. 4, pp. 545–558, 2010. View at: Publisher Site  Google Scholar
 A. M. Elshurafa, K. Khirallah, H. H. Tawfik, A. Emira, A. K. S. Abdel Aziz, and S. M. Sedky, “Nonlinear dynamics of spring softening and hardening in foldedmems comb drive resonators,” Journal of Microelectromechanical Systems, vol. 20, no. 4, Article ID 5772897, pp. 943–958, 2011. View at: Publisher Site  Google Scholar
 I. Neri, F. Travasso, H. Vocca, and L. Gammaitoni, “Nonlinear noise harvesters for nanosensors,” Nano Communication Networks, vol. 2, no. 4, pp. 230–234, 2011. View at: Publisher Site  Google Scholar
 P. Harte, E. Blokhina, O. Feely, D. FournierPrunaret, and D. Galayko, “Electrostatic vibration energy harvesters with linear and nonlinear resonators,” International Journal of Bifurcation and Chaos, vol. 24, no. 11, Article ID 1430030, 2014. View at: Publisher Site  Google Scholar
 H. Li, S. Preidikman, B. Balachandran, and C. D. Mote Jr., “Nonlinear free and forced oscillations of piezoelectric microresonators,” Journal of Micromechanics and Microengineering, vol. 16, no. 2, pp. 356–367, 2006. View at: Publisher Site  Google Scholar
 L. Gammaitoni, I. Neri, and H. Vocca, “Nonlinear oscillators for vibration energy harvesting,” Applied Physics Letters, vol. 94, no. 16, Article ID 164102, 2009. View at: Publisher Site  Google Scholar
 S. D. Nguyen and E. Halvorsen, “Nonlinear springs for bandwidthtolerant vibration energy harvesting,” Journal of Microelectromechanical Systems, vol. 20, no. 6, Article ID 6062373, pp. 1225–1227, 2011. View at: Publisher Site  Google Scholar
 L.C. Julin Blystad and E. Halvorsen, “A piezoelectric energy harvester with a mechanical end stop on one side,” Microsystem Technologies, vol. 17, no. 4, pp. 505–511, 2011. View at: Publisher Site  Google Scholar
 H. Liu, C. Lee, T. Kobayashi, C. J. Tay, and C. Quan, “Investigation of a MEMS piezoelectric energy harvester system with a frequencywidenedbandwidth mechanism introduced by mechanical stoppers,” Smart Materials and Structures, vol. 21, no. 3, Article ID 035005, 2012. View at: Publisher Site  Google Scholar
 C. P. Le, E. Halvorsen, O. Søråsen, and E. M. Yeatman, “Microscale electrostatic energy harvester using internal impacts,” Journal of Intelligent Material Systems and Structures, vol. 23, no. 13, pp. 1409–1421, 2012. View at: Publisher Site  Google Scholar
 P. Basset, D. Galayko, F. Cottone et al., “Electrostatic vibration energy harvester with combined effect of electrical nonlinearities and mechanical impact,” Journal of Micromechanics and Microengineering, vol. 24, no. 3, Article ID 035001, 2014. View at: Publisher Site  Google Scholar
 Y. Yang and D. Upadrashta, “Modeling of geometric, material and damping nonlinearities in piezoelectric energy harvesters,” Nonlinear Dynamics, vol. 84, no. 4, pp. 2487–2504, 2016. View at: Publisher Site  Google Scholar
 K. Kanda, T. Saito, Y. Iga, K. Higuchi, and K. Maenaka, “Influence of parasitic capacitance on output voltage for seriesconnected thinfilm piezoelectric devices,” Sensors, vol. 12, no. 12, pp. 16673–16684, 2012. View at: Publisher Site  Google Scholar
 L. G. W. Tvedt, D. S. Nguyen, and E. Halvorsen, “Nonlinear behavior of an electrostatic energy harvester under wide and narrowband excitation,” Journal of Microelectromechanical Systems, vol. 19, no. 2, Article ID 5404427, pp. 305–316, 2010. View at: Publisher Site  Google Scholar
 R. Perez, M. Kral, and H. Bleuler, “Study of polyvinylidene fluoride (PVDF) based bimorph actuators for laser scanning actuation at kHz frequency range,” Sensors and Actuators A: Physical, vol. 183, pp. 84–94, 2012. View at: Publisher Site  Google Scholar
 P. Ueberschlag, “PVDF piezoelectric polymer,” Sensor Review, vol. 21, no. 2, pp. 118–125, 2001. View at: Publisher Site  Google Scholar
 M. Dietze and M. EsSouni, “Structural and functional properties of screenprinted PZTPVDFTrFE composites,” Sensors and Actuators A: Physical, vol. 143, no. 2, pp. 329–334, 2008. View at: Publisher Site  Google Scholar
 C. Maruccio and L. De Lorenzis, “Numerical homogenization of piezoelectric textiles for energy harvesting,” Frattura ed Integritá Strutturale, vol. 29, pp. 49–60, 2014. View at: Google Scholar
 L. Persano, C. Dagdeviren, C. Maruccio, L. De Lorenzis, and D. Pisignano, “Cooperativity in the enhanced piezoelectric response of polymer nanowires,” Advanced Materials, vol. 26, no. 45, pp. 7574–7580, 2014. View at: Publisher Site  Google Scholar
 C. Maruccio, L. De Lorenzis, L. Persano, and D. Pisignano, “Computational homogenization of fibrous piezoelectric materials,” Computational Mechanics, vol. 55, no. 5, pp. 983–998, 2015. View at: Publisher Site  Google Scholar  MathSciNet
 J. Yang, An Introduction to The Theory of Piezoelectricity, Springer Science and Business Media, Berlin, Germany, 2005, Anintroductiontothetheoryofpiezoelectricity. Springer Science Business Media.
 G. Brandstetter and S. Govindjee, “Cyclic steady states of nonlinear electromechanical devices excited at resonance,” International Journal for Numerical Methods in Engineering, vol. 110, no. 13, pp. 1227–1246, 2016. View at: Publisher Site  Google Scholar  MathSciNet
 P. Steinmann, “Computational nonlinear electroelasticity getting Started,” CISM International Centre for Mechanical Sciences, vol. 527, pp. 181–230, 2011. View at: Google Scholar
 J. Korelc and P. Wriggers, Automation of Finite Element Methods: Basic Equations of Continuum Mechanics, Springer International Publishing, 1 edition, 2016. View at: Publisher Site
 A. Erturk and D. J. Inman, “A distributed parameter electromechanical model for cantilevered piezoelectric energy harvesters,” Journal of Vibration and Acoustics, vol. 130, Article ID 041002, pp. 1–15, 2008. View at: Publisher Site  Google Scholar
 S. C. Stanton, A. Erturk, B. P. Mann, and D. J. Inman, “Nonlinear piezoelectricity in electroelastic energy harvesters: modeling and experimental identification,” Journal of Applied Physics, vol. 108, no. 7, Article ID 074903, 2010. View at: Publisher Site  Google Scholar
 K.J. Bathe and S. Gracewski, “On nonlinear dynamic analysis using substructuring and mode superposition,” Computers & Structures, vol. 13, no. 56, pp. 699–707, 1981. View at: Publisher Site  Google Scholar
 ADINA Theory and Modeling Guide. ADINA Inc.71, Elton Avenue, Watertown, MA 02472 USA.
 Z.Q. Qu, Model Order Reduction Techniques with Applications in Finite Element Analysis, SpringerVerlag, London, UK, 2004.
 O. Weeger, U. Wever, and B. Simeon, “On the use of modal derivatives for nonlinear model order reduction,” International Journal for Numerical Methods in Engineering, vol. 108, no. 13, pp. 1579–1602, 2015. View at: Publisher Site  Google Scholar  MathSciNet
 L. D. Gabbay, J. E. Mehner, and S. D. Senturia, “Computeraided generation of nonlinear reducedorder dynamic macromodelsI: nonstressstiffened case,” Journal of Microelectromechanical Systems, vol. 9, no. 2, pp. 262–269, 2000. View at: Publisher Site  Google Scholar
 J. E. Mehner, L. D. Gabbay, and S. D. Senturia, “Computeraided generation of nonlinear reducedorder dynamic macromodelsII: stressstiffened case,” Journal of Microelectromechanical Systems, vol. 9, no. 2, pp. 270–278, 2000. View at: Publisher Site  Google Scholar
 AnsysCoupled field analysis guide. ANSYS INC., Canonsburg (PA).
 R. H. Enns and G. C. McGuire, Nonlinear Physics with Mathematica for Scientists and Engineers, Springer International Publishing, 2001. View at: Publisher Site
Copyright
Copyright © 2018 Claudio Maruccio et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.