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

Authors: Claudio Maruccio; Giuseppe Quaranta; Pasquale Montegiglio; Francesco Trentadue; Giuseppe Acciani
Journal: Shock and Vibration
 (2018)

---

## 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 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.


---

## Body

## 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 [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 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, 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 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 [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 byB and S, respectively, where B,S⊂R3. When the electromechanical body deforms, the nonlinear mapping function φ:B→S at time instant t maps the material point X∈B onto x∈S:(1)x=φX,t.The displacement vector u is obtained as the difference between the positions vectors of the current and initial configuration:(2)uX,t=φX,t-X,whereas the deformation gradient F can be defined as a function of the displacement gradient H:(3)F=gradφX,t=∂x∂X=1+H,where(4)H=gradu.According to Faraday’s law,(5)curle⃑=0,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 ϕ:(6)e⃑=-gradϕ=-∂ϕ∂x.Velocity and acceleration of a material point with respect to the reference configuration are defined, respectively, by the following material time derivatives:(7)VX,t=φ˙X,t,AX,t=φ¨X,t=V˙X,t.In the current configuration, the balance of momentum and Gauss’s law state that(8)div⁡σ+ρm,vb¯v=ρm,vA,div⁡d⃑=ϱe,v,where σ represents the mechanical Cauchy stress tensor, d⃑ denotes the electric displacement, ρm,v is the mechanical density, ϱe,v is free electric charge density, and ρm,vb¯v indicates the volume force (in the current configuration). The balance of mass implies that ρm,v=ρm,V/JF, where ρm,V is the density in the initial configuration and JF=det⁡F. Local balance of angular momentum guarantees that σ=σT. Moreover, if  p⃑ is the electric polarization vector in the current configuration and ϵ0 is the vacuum permittivity, then(9)d⃑=ϵ0e⃑+p⃑.The transformation between volume elements dv, dV and electric charges ϱe,v, ϱe,V in the current and reference configurations, respectively, is based on the following equations:(10)ϱe,v=ϱe,VJF,dv=JFdV.Furthermore, the condition JF≠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:(11)div⁡P+ρm,Vb¯V=ρm,VA,div⁡FS+ρm,Vb¯V=ρm,VA,where P and S are the total first and second Piola–Kirchoff stress tensors, respectively. Moreover, ρm,Vb¯V 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(12)E=12FFT-I=C-I2,where(13)C=FTFis the right Cauchy–Green tensor. Objectivity requires that ψ~F,E→=objψC,E→=ψ-E,E→. A compressible Neo–Hookean type material model with a total energy density ψC,E→ is used in this work:(14)ψC,E→=μ12I1-3-ln⁡JF+λ2ln⁡JF2+C1I4+C2I5,where λ and μ are the Lame constants, C1 and C2 are further material constants to be calibrated, and I1, I4, and I5 are computed for a transversely isotropic material according to [45, 46] and they are equal to (15)I1=C:I=TrC,I4=E→⊗E→:I=TrE→⊗E→,I5=E→⊗E→:C=TrCE→⊗E→,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:(16)S=∂ψ∂E=2∂ψ∂C,D⃑=-∂ψ∂E⃑.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(17)P=FS,σ=FSFTJF,respectively. The transformation from the material to the current configuration is possible by means of the following relationships:(18)e⃑=F-TE⃑,d⃑=FD⃑JF.The Dirichlet and Neumann boundary conditions for the mechanical field are(19)u=u¯onΓut=P·N=t¯onΓt,where u¯ and t¯ are prescribed mechanical displacement and surface traction vectors in the reference configuration, respectively. The boundary of the domain is Γ (Γu and Γt are its Dirichlet and Neumann portions, resp.), with Γ=Γu∪Γt and Γu∩Γt=⊘. Moreover, N is the outward unit normal to Γ. The boundary conditions for the electric field are(20)ϕ=ϕ¯onΓϕ,d→=D→·N=d¯→onΓd→,where ϕ¯ and d¯→ are prescribed values of electric potential on Γϕ and electric charge flux on Γd→, respectively. Moreover, Γ=Γϕ∪Γd→, Γϕ∩Γd→=⊘.A cantilever-type configuration is here assumed for the energy harvesting device (see Figure1). 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:(21)ψF=12λJF-12+μ12TrC-3-ln⁡JF.Figure 1
Schematic representation of the unimorph electromechanical generator with a piezoelectric layer made of PVDF.
(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 u=∑Niu^i and ϕ=∑Niϕ^i be the discretized displacement and electric potential fields, respectively, where u^i are the nodal displacements and ϕ^i are the nodal electric potentials (Ni 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∈Rntp have to be determined (ntp is the total number of global unknowns of the problem). If pe⊂p is a subset of the global vector of unknowns p on which the eth element depends explicitly, then the element Re and the Gauss point Rg contributions to the global residuals R are explicit functions of pe, that is, Re(pe) and Rg(pe). In particular, at the element level, the internal residuals Re are obtained using AceGen as follows [47]:(22)Re=∂∂pe∫VψpedV=0,where pe is the unknown vector related to the element that collects all nodal displacements u^i and/or nodal electrical potentials ϕ^i. Within the FE procedure, the global residuals R are approximated as(23)R≈⋀e=1neRe=⋀e=1ne∑g=1ngwgJgRg=0,where ⋀ is the standard FE assembly operator, ne indicates the number of elements, ng is the number of Gauss points, and g indicates a generic Gauss point (wg and Jg are the Gauss point weight and Jacobian determinant, resp.). The Gauss point contribution to the residuals is Rg=δ^ψ(pe)/δ^pe.The requirement of zero global residualsR(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](24)Rp,Λ=Rintp-ΛRref=0,where Rint denotes the contribution of internal forces to the nodal force vector and Rref 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(25)Rpι+Δpι=Rpι+DRpι·Δpι+rpι,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(ι)=DRpι. Within the FE procedure, the tangent operator K≈⋀e=1neKe=⋀e=1ne∑g=1ngwgKg is formed from the Gauss point tangent operator Kg:(26)Kg=δ^Rgδ^pe.


### 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 unknownsu^i and ϕ^i. 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 ψm(F) and refers to the mechanical domain. Thus,(27)Π=⋃ψemC,E⃑+ψmF.Consequently, the following equations are derived: (28)Rui=δΠδu^i,Rϕi=δΠδϕ^i,Kuuij=δRuiδu^j,Kϕϕij=δRϕiδϕ^j,Kuϕij=δRuiδϕ^j,Kϕuij=δRϕiδu^j.Using an implicit time integration scheme and Newton’s method iterations to determine the dynamic equilibrium at time t+Δt, the coupled nonlinear FE equations are(29)Muut+Δt000u¨dιϕ¨ι+Cuut+Δt000u˙dιϕ˙ι+KuutKuϕtKϕut-KϕϕtΔudιΔϕι=fut+Δtι-1fϕt+Δtι-1,where ud 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 Muu, Cuu, and Kuu are the structural mass, damping, and stiffness matrices, respectively. The matrices Kuϕ and Kϕu are the stiffness matrices due to piezoelectric mechanical coupling. The matrix Kϕϕ is the stiffness matrix resulting from the electrical fields. The vectors fu 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 Cuu is considered:(30)Cuu=χMuu+βKuu,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 ηr as follows:(31)η¨r+2ζrωrη˙r+ωr2ηr-θrυ=fr,Cpeqυ˙+υRL+∑r=1∞τrη˙r=0,for each mode r, where ζr is the modal mechanical damping ratio, ωr is the undamped natural frequency, θr and τr are modal electromechanical coupling terms, fr is the modal mechanical forcing function, Cpeq is the capacitance, RL is the load resistance, and υ is the voltage response across the external resistive load. Having calculated the modal coordinates ηr(t), the transverse displacement of the neutral axis u(ξ,t) relative to the moving base at position ξ=x/L and time t is equal to(32)uξ,t=∑r=1∞ϕrξηrt,where ϕr(ξ) 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(33)η¨+2ζωη˙+ω2η+αη3-θ+φη2υ=Tnz¨g,(34)Cpeqυ˙+υRL+τ+κη2η˙=0,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 Tn is the modal mechanical forcing function coefficient that multiplies the base acceleration z¨g. 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(35)Cpeq=ϵ33BLTPVDF,where B, L, and TPVDF are the width, length, and thickness of the piezoelectric layer, respectively, whereas ϵ33 is the permittivity constant. The rth modal electromechanical coupling term for unimorph configuration can be obtained as [48](36)θr=EPVDFd31B2hb+TPVDFϕr′L2,where d31 is the piezoelectric constant and hb is the position of the bottom PVDF layer from the neutral axis. It is understood that θr=1=θ in (33). Following [48], hb=hpa-TPVDF, where hpa is the distance from the top of the PVDF layer to the neutral axis, and it is equal to(37)hpa=TPVDF2+2ℵTmylarTPVDF+ℵTPVDF22Tmylar+ℵTPVDF,with ℵ=Emylar/EPVDF, where Emylar and EPVDF 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, Tmylar is the thickness of the substrate.


## 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 byB and S, respectively, where B,S⊂R3. When the electromechanical body deforms, the nonlinear mapping function φ:B→S at time instant t maps the material point X∈B onto x∈S:(1)x=φX,t.The displacement vector u is obtained as the difference between the positions vectors of the current and initial configuration:(2)uX,t=φX,t-X,whereas the deformation gradient F can be defined as a function of the displacement gradient H:(3)F=gradφX,t=∂x∂X=1+H,where(4)H=gradu.According to Faraday’s law,(5)curle⃑=0,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 ϕ:(6)e⃑=-gradϕ=-∂ϕ∂x.Velocity and acceleration of a material point with respect to the reference configuration are defined, respectively, by the following material time derivatives:(7)VX,t=φ˙X,t,AX,t=φ¨X,t=V˙X,t.In the current configuration, the balance of momentum and Gauss’s law state that(8)div⁡σ+ρm,vb¯v=ρm,vA,div⁡d⃑=ϱe,v,where σ represents the mechanical Cauchy stress tensor, d⃑ denotes the electric displacement, ρm,v is the mechanical density, ϱe,v is free electric charge density, and ρm,vb¯v indicates the volume force (in the current configuration). The balance of mass implies that ρm,v=ρm,V/JF, where ρm,V is the density in the initial configuration and JF=det⁡F. Local balance of angular momentum guarantees that σ=σT. Moreover, if  p⃑ is the electric polarization vector in the current configuration and ϵ0 is the vacuum permittivity, then(9)d⃑=ϵ0e⃑+p⃑.The transformation between volume elements dv, dV and electric charges ϱe,v, ϱe,V in the current and reference configurations, respectively, is based on the following equations:(10)ϱe,v=ϱe,VJF,dv=JFdV.Furthermore, the condition JF≠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:(11)div⁡P+ρm,Vb¯V=ρm,VA,div⁡FS+ρm,Vb¯V=ρm,VA,where P and S are the total first and second Piola–Kirchoff stress tensors, respectively. Moreover, ρm,Vb¯V 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(12)E=12FFT-I=C-I2,where(13)C=FTFis the right Cauchy–Green tensor. Objectivity requires that ψ~F,E→=objψC,E→=ψ-E,E→. A compressible Neo–Hookean type material model with a total energy density ψC,E→ is used in this work:(14)ψC,E→=μ12I1-3-ln⁡JF+λ2ln⁡JF2+C1I4+C2I5,where λ and μ are the Lame constants, C1 and C2 are further material constants to be calibrated, and I1, I4, and I5 are computed for a transversely isotropic material according to [45, 46] and they are equal to (15)I1=C:I=TrC,I4=E→⊗E→:I=TrE→⊗E→,I5=E→⊗E→:C=TrCE→⊗E→,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:(16)S=∂ψ∂E=2∂ψ∂C,D⃑=-∂ψ∂E⃑.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(17)P=FS,σ=FSFTJF,respectively. The transformation from the material to the current configuration is possible by means of the following relationships:(18)e⃑=F-TE⃑,d⃑=FD⃑JF.The Dirichlet and Neumann boundary conditions for the mechanical field are(19)u=u¯onΓut=P·N=t¯onΓt,where u¯ and t¯ are prescribed mechanical displacement and surface traction vectors in the reference configuration, respectively. The boundary of the domain is Γ (Γu and Γt are its Dirichlet and Neumann portions, resp.), with Γ=Γu∪Γt and Γu∩Γt=⊘. Moreover, N is the outward unit normal to Γ. The boundary conditions for the electric field are(20)ϕ=ϕ¯onΓϕ,d→=D→·N=d¯→onΓd→,where ϕ¯ and d¯→ are prescribed values of electric potential on Γϕ and electric charge flux on Γd→, respectively. Moreover, Γ=Γϕ∪Γd→, Γϕ∩Γd→=⊘.A cantilever-type configuration is here assumed for the energy harvesting device (see Figure1). 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:(21)ψF=12λJF-12+μ12TrC-3-ln⁡JF.Figure 1
Schematic representation of the unimorph electromechanical generator with a piezoelectric layer made of PVDF.
(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 u=∑Niu^i and ϕ=∑Niϕ^i be the discretized displacement and electric potential fields, respectively, where u^i are the nodal displacements and ϕ^i are the nodal electric potentials (Ni 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∈Rntp have to be determined (ntp is the total number of global unknowns of the problem). If pe⊂p is a subset of the global vector of unknowns p on which the eth element depends explicitly, then the element Re and the Gauss point Rg contributions to the global residuals R are explicit functions of pe, that is, Re(pe) and Rg(pe). In particular, at the element level, the internal residuals Re are obtained using AceGen as follows [47]:(22)Re=∂∂pe∫VψpedV=0,where pe is the unknown vector related to the element that collects all nodal displacements u^i and/or nodal electrical potentials ϕ^i. Within the FE procedure, the global residuals R are approximated as(23)R≈⋀e=1neRe=⋀e=1ne∑g=1ngwgJgRg=0,where ⋀ is the standard FE assembly operator, ne indicates the number of elements, ng is the number of Gauss points, and g indicates a generic Gauss point (wg and Jg are the Gauss point weight and Jacobian determinant, resp.). The Gauss point contribution to the residuals is Rg=δ^ψ(pe)/δ^pe.The requirement of zero global residualsR(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](24)Rp,Λ=Rintp-ΛRref=0,where Rint denotes the contribution of internal forces to the nodal force vector and Rref 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(25)Rpι+Δpι=Rpι+DRpι·Δpι+rpι,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(ι)=DRpι. Within the FE procedure, the tangent operator K≈⋀e=1neKe=⋀e=1ne∑g=1ngwgKg is formed from the Gauss point tangent operator Kg:(26)Kg=δ^Rgδ^pe.


## 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 unknownsu^i and ϕ^i. 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 ψm(F) and refers to the mechanical domain. Thus,(27)Π=⋃ψemC,E⃑+ψmF.Consequently, the following equations are derived: (28)Rui=δΠδu^i,Rϕi=δΠδϕ^i,Kuuij=δRuiδu^j,Kϕϕij=δRϕiδϕ^j,Kuϕij=δRuiδϕ^j,Kϕuij=δRϕiδu^j.Using an implicit time integration scheme and Newton’s method iterations to determine the dynamic equilibrium at time t+Δt, the coupled nonlinear FE equations are(29)Muut+Δt000u¨dιϕ¨ι+Cuut+Δt000u˙dιϕ˙ι+KuutKuϕtKϕut-KϕϕtΔudιΔϕι=fut+Δtι-1fϕt+Δtι-1,where ud 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 Muu, Cuu, and Kuu are the structural mass, damping, and stiffness matrices, respectively. The matrices Kuϕ and Kϕu are the stiffness matrices due to piezoelectric mechanical coupling. The matrix Kϕϕ is the stiffness matrix resulting from the electrical fields. The vectors fu 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 Cuu is considered:(30)Cuu=χMuu+βKuu,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 ηr as follows:(31)η¨r+2ζrωrη˙r+ωr2ηr-θrυ=fr,Cpeqυ˙+υRL+∑r=1∞τrη˙r=0,for each mode r, where ζr is the modal mechanical damping ratio, ωr is the undamped natural frequency, θr and τr are modal electromechanical coupling terms, fr is the modal mechanical forcing function, Cpeq is the capacitance, RL is the load resistance, and υ is the voltage response across the external resistive load. Having calculated the modal coordinates ηr(t), the transverse displacement of the neutral axis u(ξ,t) relative to the moving base at position ξ=x/L and time t is equal to(32)uξ,t=∑r=1∞ϕrξηrt,where ϕr(ξ) 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(33)η¨+2ζωη˙+ω2η+αη3-θ+φη2υ=Tnz¨g,(34)Cpeqυ˙+υRL+τ+κη2η˙=0,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 Tn is the modal mechanical forcing function coefficient that multiplies the base acceleration z¨g. 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(35)Cpeq=ϵ33BLTPVDF,where B, L, and TPVDF are the width, length, and thickness of the piezoelectric layer, respectively, whereas ϵ33 is the permittivity constant. The rth modal electromechanical coupling term for unimorph configuration can be obtained as [48](36)θr=EPVDFd31B2hb+TPVDFϕr′L2,where d31 is the piezoelectric constant and hb is the position of the bottom PVDF layer from the neutral axis. It is understood that θr=1=θ in (33). Following [48], hb=hpa-TPVDF, where hpa is the distance from the top of the PVDF layer to the neutral axis, and it is equal to(37)hpa=TPVDF2+2ℵTmylarTPVDF+ℵTPVDF22Tmylar+ℵTPVDF,with ℵ=Emylar/EPVDF, where Emylar and EPVDF 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, Tmylar 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 Figure2). 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 (Fz) is statically applied to the structural model (including nonlinear effects) and the total reaction Rz 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 αD=αm and κD=ω2m, where m 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 αD 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).Figure 2
Implemented hybrid computational strategy.


## 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 d31=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 Table 1, together with the other main material and geometrical data adopted in this numerical study.Table 1
Material and geometrical data used in the numerical study.
Property Symbol Value Elastic modulus of the piezoelectric film E P V D F 3.0 GPa Elastic modulus of the substrate E m y l a r 3.79 GPa Poisson ratio of the piezoelectric film ν P V D F 0.3 Poisson ratio of the substrate ν m y l a r 0.35 Damping ratio ζ 4.0% Density of the piezoelectric film ρ P V D F 1780 kg/m3 Density of the substrate ρ m y l a r 1390 kg/m3 Lumped tip mass M l u m p 0.032 g Piezoelectric coupling coefficient d 31 23 ⋅ 1 0 - 12 pC/N Electrical permittivity ϵ 33 106e-12 F/m Unimorph length L 31.7 mm Unimorph width B 16.0 mm Circuit resistance R L 10 MΩ Thickness of the piezoelectric film T P V D F 28μm Homogenized stiffness E I 37.7 Nmm2 Equivalent mylar layer thickness T m y l a r e q 172μm


### 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 L, width B, and thickness Ts. 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 Es, rotational inertia Is, and cross-sectional areas As = B·Ts. The subscript s denotes the sth 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 zN is then given by [12](38)EI=∑s=1NlayerAsEszs-zN2+Ts212,where(39)zN=∑s=1NlayerEsTszs∑s=1NlayerEsTsand zs is the location of the axis of the sth layer with respect to an arbitrary reference. Nlayer is the number of layers considered (5 for our device). The homogenized stiffness EI 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 Tmylareq 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 Tmylareq 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.Figure 3
First four modal shapes and modal frequencies of the considered piezoelectric device: (a)ω1 = 316 rad/sec, (b) ω2 = 1980 rad/sec, (c) ω3 = 5545 rad/sec, and (d) ω4 = 10866 rad/sec.
(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 Figure4. 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 L while the load Fz is normalized by L2/(EI), where EI 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 σyy and deformed shape are shown in Figure 6.Figure 4
Details of the FE mesh discretization.Figure 5
Vertical tip deflection of the cantilever unimorph: comparison between FE solutions and experimental data.Figure 6
Contour levels of vertical displacement components and Cauchy stresses (color bars refer to the final load step).


### 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 endRz 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 (40)Rz=κDδtip+αDδtip3,where κD and αD characterize the constitutive law of the nonlinear spring. In the static case, κD=3EI/L3. Figure 7 provides the comparison between the FE analysis and the approximation obtained by means of the analytical Duffing model.Figure 7
FE-based capacity curve and analytical approximation.


### 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 ofEPVDF, Emylar, TPVDF, and Tmylar 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 RL=10 MΩ are determined for four values of EPVDF, Emylar, TPVDF, and Tmylar. The reference values for EPVDF, Emylar, TPVDF, and Tmylar are those listed in Table 1 and are indicated as EPVDF0, Emylar0, TPVDF0, and Tmylar0, 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 αD for each case.Table 2
Computation ofαD coefficients [10-6 N/mm3].
E m y l a r 0 % → α D E P V D F 0 % → α D T m y l a r 0 % → α D T P V D F 0 % → α D 50 → 6.0 120 → 7.5 80 → 4.0 50 → 5.5 100 → 8.5 100 → 8.5 100 → 8.5 100 → 8.5 150 → 8.5 140 → 8.0 90 → 5.5 120 → 6.5 200 → 8.5 160 → 8.5 120 → 11.0 140  → 8.0Figure 8
FE-based capacity curves: parametric analysis.
(a) (b) (c) (d)Figures9 and 10 highlight the role of the piezoelectric elastic modulus EPVDF whereas that of the substrate Emylar 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 1ag, with ag=9.8 m/s2). For very large deformations (intensity of the base excitation equal to 10ag), 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.Figure 9
FRFs of the maximum tip displacement for base acceleration amplitudes ranging from 1ag to 10ag, where (a) EPVDF=EPVDF0, (b) EPVDF = 1.2EPVDF0, (c) EPVDF = 1.4EPVDF0, and (d) EPVDF = 1.6EPVDF0.
(a) (b) (c) (d)Figure 10
FRFs of the maximum voltage for base acceleration amplitudes ranging from 1ag to 10ag, where (a) EPVDF=EPVDF0, (b) EPVDF = 1.2EPVDF0, (c) EPVDF = 1.4EPVDF0, and (d) EPVDF = 1.6EPVDF0.
(a) (b) (c) (d)Figure 11
FRFs of the maximum tip displacement for base acceleration amplitudes ranging from 1ag to 10ag, where (a) Emylar=Emylar0, (b) Emylar = 0.5Emylar0, (c) Emylar = 1.5Emylar0, and (d) Emylar = 2.0Emylar0.
(a) (b) (c) (d)Figure 12
FRFs of the maximum voltage for base acceleration amplitudes ranging from 1ag to 10ag, where (a) Emylar=Emylar0, (b) Emylar = 0.5Emylar0, (c) Emylar = 1.5Emylar0, and (d) Emylar = 2.0Emylar0.
(a) (b) (c) (d)Figure 13
FRFs of the maximum tip displacement for base acceleration amplitudes ranging from 1ag to 10ag, where (a) TPVDF=TPVDF0, (b) TPVDF=0.5TPVDF0, (c) TPVDF=1.2TPVDF0, and (d) TPVDF=1.4TPVDF0.
(a) (b) (c) (d)Figure 14
FRFs of the maximum voltage for base acceleration amplitudes ranging from 1ag to 10ag, where (a) TPVDF=TPVDF0, (b) TPVDF=0.5TPVDF0, (c) TPVDF=1.2TPVDF0, and (d) TPVDF=1.4TPVDF0.
(a) (b) (c) (d)Figure 15
FRFs of the maximum tip displacement for base acceleration amplitudes ranging from 1ag to 10ag, where (a) Tmylar=Tmylar0, (b) Tmylar=0.8Tmylar0, (c) Tmylar=0.9Tmylar0, and (d) Tmylar=1.2Tmylar0.
(a) (b) (c) (d)Figure 16
FRFs of the maximum voltage for base acceleration amplitudes ranging from 1ag to 10ag, where (a) Tmylar=Tmylar0, (b) Tmylar=0.8Tmylar0, (c) Tmylar=0.9Tmylar0, and (d) Tmylar=1.2Tmylar0.
(a) (b) (c) (d)


## 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 d31=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 Table 1, together with the other main material and geometrical data adopted in this numerical study.Table 1
Material and geometrical data used in the numerical study.
Property Symbol Value Elastic modulus of the piezoelectric film E P V D F 3.0 GPa Elastic modulus of the substrate E m y l a r 3.79 GPa Poisson ratio of the piezoelectric film ν P V D F 0.3 Poisson ratio of the substrate ν m y l a r 0.35 Damping ratio ζ 4.0% Density of the piezoelectric film ρ P V D F 1780 kg/m3 Density of the substrate ρ m y l a r 1390 kg/m3 Lumped tip mass M l u m p 0.032 g Piezoelectric coupling coefficient d 31 23 ⋅ 1 0 - 12 pC/N Electrical permittivity ϵ 33 106e-12 F/m Unimorph length L 31.7 mm Unimorph width B 16.0 mm Circuit resistance R L 10 MΩ Thickness of the piezoelectric film T P V D F 28μm Homogenized stiffness E I 37.7 Nmm2 Equivalent mylar layer thickness T m y l a r e q 172μm


## 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 L, width B, and thickness Ts. 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 Es, rotational inertia Is, and cross-sectional areas As = B·Ts. The subscript s denotes the sth 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 zN is then given by [12](38)EI=∑s=1NlayerAsEszs-zN2+Ts212,where(39)zN=∑s=1NlayerEsTszs∑s=1NlayerEsTsand zs is the location of the axis of the sth layer with respect to an arbitrary reference. Nlayer is the number of layers considered (5 for our device). The homogenized stiffness EI 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 Tmylareq 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 Tmylareq 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.Figure 3
First four modal shapes and modal frequencies of the considered piezoelectric device: (a)ω1 = 316 rad/sec, (b) ω2 = 1980 rad/sec, (c) ω3 = 5545 rad/sec, and (d) ω4 = 10866 rad/sec.
(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 Figure4. 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 L while the load Fz is normalized by L2/(EI), where EI 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 σyy and deformed shape are shown in Figure 6.Figure 4
Details of the FE mesh discretization.Figure 5
Vertical tip deflection of the cantilever unimorph: comparison between FE solutions and experimental data.Figure 6
Contour levels of vertical displacement components and Cauchy stresses (color bars refer to the final load step).


## 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 endRz 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 (40)Rz=κDδtip+αDδtip3,where κD and αD characterize the constitutive law of the nonlinear spring. In the static case, κD=3EI/L3. Figure 7 provides the comparison between the FE analysis and the approximation obtained by means of the analytical Duffing model.Figure 7
FE-based capacity curve and analytical approximation.


## 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 ofEPVDF, Emylar, TPVDF, and Tmylar 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 RL=10 MΩ are determined for four values of EPVDF, Emylar, TPVDF, and Tmylar. The reference values for EPVDF, Emylar, TPVDF, and Tmylar are those listed in Table 1 and are indicated as EPVDF0, Emylar0, TPVDF0, and Tmylar0, 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 αD for each case.Table 2
Computation ofαD coefficients [10-6 N/mm3].
E m y l a r 0 % → α D E P V D F 0 % → α D T m y l a r 0 % → α D T P V D F 0 % → α D 50 → 6.0 120 → 7.5 80 → 4.0 50 → 5.5 100 → 8.5 100 → 8.5 100 → 8.5 100 → 8.5 150 → 8.5 140 → 8.0 90 → 5.5 120 → 6.5 200 → 8.5 160 → 8.5 120 → 11.0 140  → 8.0Figure 8
FE-based capacity curves: parametric analysis.
(a) (b) (c) (d)Figures9 and 10 highlight the role of the piezoelectric elastic modulus EPVDF whereas that of the substrate Emylar 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 1ag, with ag=9.8 m/s2). For very large deformations (intensity of the base excitation equal to 10ag), 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.Figure 9
FRFs of the maximum tip displacement for base acceleration amplitudes ranging from 1ag to 10ag, where (a) EPVDF=EPVDF0, (b) EPVDF = 1.2EPVDF0, (c) EPVDF = 1.4EPVDF0, and (d) EPVDF = 1.6EPVDF0.
(a) (b) (c) (d)Figure 10
FRFs of the maximum voltage for base acceleration amplitudes ranging from 1ag to 10ag, where (a) EPVDF=EPVDF0, (b) EPVDF = 1.2EPVDF0, (c) EPVDF = 1.4EPVDF0, and (d) EPVDF = 1.6EPVDF0.
(a) (b) (c) (d)Figure 11
FRFs of the maximum tip displacement for base acceleration amplitudes ranging from 1ag to 10ag, where (a) Emylar=Emylar0, (b) Emylar = 0.5Emylar0, (c) Emylar = 1.5Emylar0, and (d) Emylar = 2.0Emylar0.
(a) (b) (c) (d)Figure 12
FRFs of the maximum voltage for base acceleration amplitudes ranging from 1ag to 10ag, where (a) Emylar=Emylar0, (b) Emylar = 0.5Emylar0, (c) Emylar = 1.5Emylar0, and (d) Emylar = 2.0Emylar0.
(a) (b) (c) (d)Figure 13
FRFs of the maximum tip displacement for base acceleration amplitudes ranging from 1ag to 10ag, where (a) TPVDF=TPVDF0, (b) TPVDF=0.5TPVDF0, (c) TPVDF=1.2TPVDF0, and (d) TPVDF=1.4TPVDF0.
(a) (b) (c) (d)Figure 14
FRFs of the maximum voltage for base acceleration amplitudes ranging from 1ag to 10ag, where (a) TPVDF=TPVDF0, (b) TPVDF=0.5TPVDF0, (c) TPVDF=1.2TPVDF0, and (d) TPVDF=1.4TPVDF0.
(a) (b) (c) (d)Figure 15
FRFs of the maximum tip displacement for base acceleration amplitudes ranging from 1ag to 10ag, where (a) Tmylar=Tmylar0, (b) Tmylar=0.8Tmylar0, (c) Tmylar=0.9Tmylar0, and (d) Tmylar=1.2Tmylar0.
(a) (b) (c) (d)Figure 16
FRFs of the maximum voltage for base acceleration amplitudes ranging from 1ag to 10ag, where (a) Tmylar=Tmylar0, (b) Tmylar=0.8Tmylar0, (c) Tmylar=0.9Tmylar0, and (d) Tmylar=1.2Tmylar0.
(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 Figure1(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 Huu(ω):(41)Huuω=Kuu-ω2Muu+ιCuuω-1.Here, ι indicates the imaginary unit. Using the mode superposition technique, it is possible to express Huu(ω) in the form(42)Huuω=∑r=1nrmΥrΥrTωr2-ω2+ι2ζrωr,where the SRth element of the matrix Huu is obtained from the measured transfer function between points S and R. Υr indicates experimental modal shapes.Figure 17
Block diagram of the experimental setup.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 Figure18). 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.7ag and 10.5ag 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 RL=1 MΩ.Figure 18
Experimental characterization of device natural frequencies for chirp signals of 3 V and 10 V.Figure 19
FFT of 3.7ag and 10.5ag base acceleration (a-b); FFT of the measured output voltages (c-d); comparison among experimental (solid blue lines) and numerical (dashed red lines) output voltages for 3.7ag and 10.5ag base acceleration (e-f).
(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 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 RL 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.7ag and 10.5ag. 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.Table 3
Experimental versus numerical results.
Experimental voltage peak [V] Numerical voltage peak [V] 3.7ag 5.2 5.2 10.5ag 10.3 10.1


## 5.1. Experimental Layout and Analysis

The proposed approach is further validated against new experimental data. The tested piezoelectric energy harvester is shown in Figure1(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 Huu(ω):(41)Huuω=Kuu-ω2Muu+ιCuuω-1.Here, ι indicates the imaginary unit. Using the mode superposition technique, it is possible to express Huu(ω) in the form(42)Huuω=∑r=1nrmΥrΥrTωr2-ω2+ι2ζrωr,where the SRth element of the matrix Huu is obtained from the measured transfer function between points S and R. Υr indicates experimental modal shapes.Figure 17
Block diagram of the experimental setup.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 Figure18). 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.7ag and 10.5ag 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 RL=1 MΩ.Figure 18
Experimental characterization of device natural frequencies for chirp signals of 3 V and 10 V.Figure 19
FFT of 3.7ag and 10.5ag base acceleration (a-b); FFT of the measured output voltages (c-d); comparison among experimental (solid blue lines) and numerical (dashed red lines) output voltages for 3.7ag and 10.5ag base acceleration (e-f).
(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 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 RL 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.7ag and 10.5ag. 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.Table 3
Experimental versus numerical results.
Experimental voltage peak [V] Numerical voltage peak [V] 3.7ag 5.2 5.2 10.5ag 10.3 10.1


## 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.



---

*Source: 2054873-2018-02-26.xml*