A Two-Step Hybrid Approach for Modeling the Nonlinear Dynamic Response of Piezoelectric Energy Harvesters

1Department of Innovation Engineering, University of Salento, Via Monteroni, 73100 Lecce, Italy 2Department of Structural and Geotechnical Engineering, Sapienza University of Rome, Via Eudossiana 18, 00184 Rome, Italy 3Department of Electrical and Information Engineering, Technical University of Bari, Via Orabona 4, 70125 Bari, Italy 4Department of Civil Engineering and Architecture, Technical University of Bari, Via Orabona 4, 70125 Bari, Italy


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][2][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 large-scale applications of this technology.
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 2 Shock and Vibration (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 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 FE-based 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.

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][42][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 B and S, respectively, where B, S ⊂ R 3 .When the electromechanical body deforms, the nonlinear mapping function  : B → S at time instant  maps the material point X ∈ B onto x ∈ S: x =  (X, ) . (1) The displacement vector u is obtained as the difference between the positions vectors of the current and initial configuration: whereas the deformation gradient F can be defined as a function of the displacement gradient H: where According to Faraday's law, where  ⇀ e is the electric field vector in the current configuration.Consequently, it is possible to define  ⇀ e 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 that where  represents the mechanical Cauchy stress tensor,  ⇀ d denotes the electric displacement,  ,V is the mechanical density,  ,V is free electric charge density, and  ,V b V indicates the volume force (in the current configuration).The balance of mass implies that  ,V =  , /  , where  , is the density in the initial configuration and   = det F. Local balance of angular momentum guarantees that  =   .Moreover, if  ⇀ p is the electric polarization vector in the current configuration and  0 is the vacuum permittivity, then The transformation between volume elements V,  and electric charges  ,V ,  , in the current and reference configurations, respectively, is based on the following equations: Furthermore, the condition   ̸ = 0 ensures that the tensor F 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: Shock and Vibration 3 where P and S are the total first and second Piola-Kirchoff stress tensors, respectively.Moreover,  , b  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, F, E, and C and the electric field vector  → E .All the quantities refer to the initial configuration.Here, E indicates the Green-Lagrange strain tensor that can be expressed as a function of the deformation gradient tensor F by where is the right Cauchy-Green tensor.Objectivity requires that where  and  are the Lame constants, C 1 and C 2 are further material constants to be calibrated, and I 1 , I 4 , and I 5 are computed for a transversely isotropic material according to [45,46] and they are equal to where I 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,  ⇀ D 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 to respectively.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 are where u and t 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, N is the outward unit normal to Γ.
The boundary conditions for the electric field are where  and 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 (F):

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 u = ∑   û 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 R(p) = 0, where the unknown variables p ∈ R  tp have to be determined ( tp is the total number of global unknowns of the problem).If p  ⊂ p is a subset of the global vector of unknowns p on which the th element depends explicitly, then the element R  and the Gauss point R  contributions to the global residuals R are explicit functions of p  , that is, R  (p  ) and R  (p  ).In particular, at the element level, the internal residuals R  are obtained using AceGen as follows [47]: where p  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 R are approximated as where ⋀ 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 R  = δ(p  )/ δp  .
The requirement of zero global residuals R(p) yields a set of nonlinear systems of equations, which is solved numerically by means of a Newton-Raphson algorithm.The system of equations R(p) = 0 is parameterized, introducing a load factor Λ in the form R(p, Λ) = 0. Thus, it is considered as [47] where R int denotes the contribution of internal forces to the nodal force vector and R ref 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 p () results in the expression where the upper index () denotes the quantities at the  iteration and Δp () = p () − p (−1) , whereas D indicates the directional derivative required for the linearization.In particular, the linearization of the vector R(p) yields the tangent matrix K = K () = DR(p () ).Within the FE procedure, the tangent operator K ≈ ⋀ =1   K  is formed from the Gauss point tangent operator K  : 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  em (C,  ⇀ E ) and refers to the electromechanical domain whereas the second term is   (F) 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 are where u  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 M  , C  , and K  are the structural mass, damping, and stiffness matrices, respectively.The matrices K  and K  are the stiffness matrices due to piezoelectric mechanical coupling.
The matrix K  is the stiffness matrix resulting from the electrical fields.The vectors f  and f  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 C  is considered: where  and  are constant multipliers that can be determined by assuming a constant damping ratio .

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,  eq  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 to where   () 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 reads where , , 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 T  is the modal mechanical forcing function coefficient that multiplies the base acceleration z  .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 is where , , and  PVDF are the width, length, and thickness of the piezoelectric layer, respectively, whereas  33 is the permittivity constant.The th modal electromechanical coupling term for unimorph configuration can be obtained as [48] where  31 is the piezoelectric constant and ℎ  is the position of the bottom PVDF layer from the neutral axis.It is understood that  =1 =  in (33).Following [48], ℎ  = ℎ pa −  PVDF , where ℎ pa is the distance from the top of the PVDF layer to the neutral axis, and it is equal to with ℵ =  mylar / PVDF , where  mylar and  PVDF 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,  mylar is the thickness of the substrate.

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 non-linear 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  tip ).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   =  2 , 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 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).

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  31 = 23 ⋅ 10 −12 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  with the other main material and geometrical data adopted in this numerical study.

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]  = where and   is the location of the axis of the th layer with respect to an arbitrary reference. layer 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  eq mylar 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  eq mylar 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.

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 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  2 /(), 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.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 ( tip ).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,   = 3/ 3 .Figure 7 provides the comparison between the FE analysis and the approximation obtained by means of the analytical Duffing model.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  PVDF whereas that of the substrate  mylar 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   = 9.8 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.

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 Voltage (V) Figure 10: FRFs of the maximum voltage for base acceleration amplitudes ranging from 1  to 10  , where (a) Here,  indicates the imaginary unit.Using the mode superposition technique, it is possible to express H  () in the form Frequency (rad/s) Tip displacement (mm)  where the SRth element of the matrix H  is obtained from the measured transfer function between points S and R. Υ  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

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.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.where  is a complex quantity expressed as follows:

Figure 1 :
Figure 1: Schematic representation of the unimorph electromechanical generator with a piezoelectric layer made of PVDF.

Figure 4 :Figure 5 :
Figure 4: Details of the FE mesh discretization.

4. 4 .Figure 6 :
Figure 6: Contour levels of vertical displacement components and Cauchy stresses (color bars refer to the final load step).

Figure 7 :
Figure 7: FE-based capacity curve and analytical approximation.

Frequency
ratios) are extracted based on the SDOF curve fitting of the FRF.Experimental modal shapes are derived considering the receptance matrix H  ():
) and 19(b) while the corresponding FFTs (see Figures19(c) and 19(d)) and time histories of the output voltage (see Figures 19(e) and 19(f)) are measured through a resistance   = 1 MΩ.

Figure 18 :Figure 19 :
Figure 18: Experimental characterization of device natural frequencies for chirp signals of 3 V and 10 V.

Table 1 :
Material and geometrical data used in the numerical study.

Table 3 :
Experimental versus numerical results.
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 Therefore, the steady-state amplitude  V of the output voltage can be computed as V (  ) =            =    tip (  ) tip  eq  [  2 + 1/ (