Koliann Mam

Michaël Peigney
email: michael.peigney@polytechnique.org

Dominique Siegert

Finite strain effects in piezoelectric energy harvesters under direct and parametric excitations

Keywords: piezoelectric materials, finite strains, energy harvesting, parametric resonance, nonlinearities

This paper addresses the dynamic behavior of piezoelectric cantilevers under base excitations. Such devices are frequently used for applications in energy harvesting. An Euler-Bernoulli model that accounts for large-deflection effects and piezoelectric nonlinearities is proposed. Closed-form expressions of the frequency response are derived, both for direct excitation (i.e. with a base acceleration transverse to the axis of the cantilever) and parametric excitation (i.e. with a base acceleration along the axis of the cantilever). Experimental results are reported and used for assessing the validity of the proposed model. Building on the model presented, some critical issues related to energy-harvesting are investigated, such as the influence of nonlinearities on the optimal load resistance, the limits of validity of linear models, and hysteresis effects in the electrical power. The efficiency of direct and parametric excitation is also compared in detail.

Introduction

Energy harvesting from ambient vibrations has become an increasingly active topic in recent years [START_REF] Priya | Energy harvesting technologies[END_REF][START_REF] Erturk | Piezoelectric energy harvesting[END_REF][START_REF] Daqaq | On the role of nonlinearities in vibratory energy harvesting: a critical review and discussion[END_REF][START_REF] Abdelkefi | Aeroelastic energy harvesting: A review[END_REF]. The overall idea is to use ambient vibrational energy as a source for operating low-power electronic devices. Energy harvesters lie on a transduction mechanism that converts mechanical energy into electric energy. Among the various options that are available, piezoelectric materials -such as PZT ceramics -are often used. Those dielectric materials have a non symmetric crystalline structure that results in a natural polarization. Through the deformation of the crystalline lattice, deforming a piezoelectric material entails a change of the polarization, which in turn produces a variation of the electric field and a current flow in a surrounding electrical circuit. Such a behavior makes for very simple mechanical designs of energy harvesters: The cantilever configuration, represented in Fig. 1, is the most frequently used design [START_REF] Roundy | A piezoelectric vibration based generator for wireless electronics[END_REF][START_REF] Du Toit | Design considerations for mems-scale piezoelectric mechanical vibration energy harvesters[END_REF][START_REF] Ajitsaria | Modeling and analysis of a bimorph piezoelectric cantilever beam for voltage generation[END_REF][START_REF] Peigney | Piezoelectric energy harvesting from traffic-induced bridge vibrations[END_REF]. It consists of a beam clamped on one side, with piezoelectric layers covering (all or part of) the top and bottom surfaces. Such a configuration is most commonly used in direct excitation, i.e. with a base acceleration that is transverse to the axis of the cantilever beam. The large majority of related studies uses a linear modeling approach, which rests on the underlying assumption that nonlinearities are sufficiently small to be neglected. Piezoelectric materials, however, have been reported to exhibit a nonlinear behavior, even at weak electric fields [START_REF] Von Wagner | Piezo-beam systems subjected to weak electric field: experiments and modelling of non-linearities[END_REF][START_REF] Stanton | Nonlinear nonconservative behavior and modeling of piezoelectric energy harvesters including proof mass effects[END_REF]. Assessing the validity of linear models thus needs to be clearly discussed, as pointed out notably in [START_REF] Erturk | An experimentally validated bimorph cantilever model for piezoelectric energy harvesting from base excitations[END_REF][START_REF] Abdelkefi | Modeling, validation, and performance of low-frequency piezoelectric energy harvesters[END_REF].

Besides direct excitation, parametric excitation is an other way of using the piezoelectric cantilever as an energy harvester [START_REF] Daqaq | Investigation of power harvesting via parametric excitations[END_REF][START_REF] Abdelkefi | Global nonlinear distributed-parameter model of parametrically excited piezoelectric energy harvesters[END_REF]. Parametric excitations is achieved by applying a base acceleration along the axis of the cantilever beam. In contrast with direct excitations, oscillations only occur if the excitation amplitude is sufficiently large. For excitations beyond that critical value, the amplitude of oscillations grows quadratically with the excitation amplitude (as far as nonlinear effects can be neglected) and can possibly reach large values [START_REF] Butikov | Parametric excitation of a linear oscillator[END_REF]. Such a behavior raised the interest of using parametric excitations for energy harvesting [START_REF] Jia | A parametrically excited vibration energy harvester[END_REF][START_REF] Jia | An auto-parametrically excited vibration energy harvester[END_REF]. Regarding piezoelectric cantilevers, parametric excitation has been far less studied than direct excitation. In particular, the comparison between direct and parametric excitations has not been fully addressed in the literature.

In order to address the issues mentioned above, we develop an Euler-Bernoulli model of piezoelectric cantilevers that accounts for large-deflection effects and piezoelectric elastic nonlinearities (Sect. 2). The predictions of that model are subsequently compared with experiments. A similar approach has been followed by Stanton et al. [START_REF] Stanton | Nonlinear piezoelectricity in electroelastic energy harvesters: Modeling and experimental identification[END_REF] for direct excitations, and by Daqaq and Stabler [START_REF] Daqaq | Investigation of power harvesting via parametric excitations[END_REF] for parametric excitations. Regarding direct excitations, Stanton et al. [START_REF] Stanton | Nonlinear piezoelectricity in electroelastic energy harvesters: Modeling and experimental identification[END_REF] took piezoelectric nonlinearities into account but ignored the geometric nonlinearities due to the large deflections. When comparing with experimental results, those authors needed to introduce a quadratic damping term to get a satisfactory match. But if the deflections are sufficiently large for quadratic damping to be significant, then there is no reason to neglect geometric nonlinearities in the first place. Geometric nonlinearities indeed give rise to a quadratic nonlinearity (through piezoelectric coupling). Regarding parametric excitations, Daqaq and Stabler [START_REF] Daqaq | Investigation of power harvesting via parametric excitations[END_REF] took geometric nonlinearities into account but ignored the piezoelectric nonlinearities. This might explain the noticeable differences obtained between the theoretical and experimental results in their paper.

The model developed in this paper accounts both for geometric and material nonlinearities. Abdelkefi et al. [START_REF] Abdelkefi | Global nonlinear distributed-parameter model of parametrically excited piezoelectric energy harvesters[END_REF][START_REF] Abdelkefi | Effects of nonlinear piezoelectric coupling on energy harvesters under direct excitation[END_REF] proposed a distributed-parameter model of piezoelectric harvesters that includes both piezoelectric and geometric nonlinearities. They performed a parametric study to investigate the effects of all the nonlinear parameters introduced. Abdelkefi et al. also argued that 3 modes (at least) are needed in the Galerkin projection in most cases. Comparison with experiments, however, is not discussed in the works by Abdelkefi et al. [START_REF] Abdelkefi | Global nonlinear distributed-parameter model of parametrically excited piezoelectric energy harvesters[END_REF][START_REF] Abdelkefi | Effects of nonlinear piezoelectric coupling on energy harvesters under direct excitation[END_REF]. Our modeling approach is similar to theirs. However, the comparison between the model and the experiments led us to make a different choice of material nonlinearities than those considered in [START_REF] Abdelkefi | Global nonlinear distributed-parameter model of parametrically excited piezoelectric energy harvesters[END_REF]: We have taken a cubic term in the expression of the stress into account, and neglected some nonlinear piezoelectric coupling terms. Our expansion of the piezoelectric enthalpy is based on the assumption that, in the range of excitations considered, mechanical nonlinearities (large strain effects) dominate electric nonlinearities (large electric field effects).

Our purpose is to derive a model that captures the main physics of the problem but still remains simple enough for analytical investigation to be tractable and parameter identification to remain simple. To that end, we use a Galerkin projection on the linear mode whose frequency is the closest to the excitation frequency. This procedure leads to a 2-degrees of freedom dynamical system. Analytical expressions of approximate steady-state solutions are derived in Sect. 3 by using a multiple time scale expansion [START_REF] Nayfeh | Nonlinear Oscillations, wiley classics library Edition[END_REF][START_REF] Nayfeh | Introduction to Perturbation Techniques[END_REF]. This enables us to obtain closed-form expressions of the frequency response, both for direct and parametric excitations. In Sect. 4 is reported an experimental study that allows the validity of the proposed model to be assessed. In Sect 5 we further explore some consequences of the model on some issues related to energy harvesting. The most common configuration used in energy harvesting is represented in Fig 1 . It consists of a cantilever beam (of length L) subjected to base excitations. Piezoelectric patches cover a portion [0, L b ] of the top and bottom surfaces of the beam. The goal of this Section is to derive a model of the piezoelectric-equipped beam that accounts for large deflections. In Sect. 2.1 we first study the purely mechanical problem of a heterogeneous beam. Specific features of the piezoelectric are introduced in Sect. 2.2. In particular, motivated by experimental observations, we consider nonlinear constitutive laws for the piezoelectric material.

Large-deflection model of a cantilever piezoelectric beam

Vibrations of an Euler-Bernoulli beam undergoing large displacements 2.1.1. Kinematics

Consider the in-plane motion of a beam under base excitations, first ignoring the piezoelectric patches. The axes (x, y) are attached to the beam and ξ(x, y, t) is the dis-placement (at time t) of a point located at (x, y) in the reference configuration. The axis y = 0 is chosen as the neutral axis of the beam. The displacement ξ is written as ξ(x, y, t) = ξ x (x, y, t)e x + ξ y (x, y, t)e y with ξ x (x, y, t) = u(x, t) -y sin Θ(x, t), ξ y (x, y, t) = v(x, t) + y cos Θ(x, t) -y.

(2.1)

In (2.1), u(x, t) is the axial displacement of the neutral axis, v(x, t) is the deflection of the neutral axis, and Θ(x, t) is the angle of rotation of the cross-sections. Eq. (2.1) corresponds to rigid body kinematics of the cross-sections. It follows from (2.1) that the Green-Lagrange strain tensor L = 1 2 (∇ξ+∇ T ξ+∇ T ξ•∇ξ) takes the form

L = L xx L xy L xy 0 (2.2) with L xx = 1 2 (e 2 -1) -yκ[(1 + u ) cos Θ + v sin Θ] + 1 2 y 2 κ 2 , L xy = 1 2 (-(1 + u ) sin Θ + v cos Θ).
(2.

3)

The superscript denotes differentiation with respect to the x direction. In (2.3), κ = Θ is the curvature of the beam and e is the axial dilatation of the beam, defined by

e = (1 + u ) 2 + v 2 .
Under the Euler-Bernoulli assumption that the cross-sections remain orthogonal to the neutral axis, the angle Θ satisfies the relations

cos Θ = 1 + u e , sin Θ = v e ;
(2.4) so that L xy = 0. We consider that the beam is inextensible, i.e. e = 1. Differentiating (2.4) with respect to x and using the inextensibility condition, we have

κ = v cos Θ = v √ 1 -v 2 .
Hence, up to the fourth-order in v 4 ,

κ = v (1 + 1 2 v 2 ). (2.5)
It follows that, up to the fourth-order in v 4 ,

L xx = -yκ = -yv (1 + 1 2 v 2 ), u = - v 2 2 (2.6)
Observe that, under the considered assumptions, the Green-Lagrange strain tensor L is equal to the infinitesimal strain tensor = 1 2 (∇ξ + ∇ T ξ), i.e. we have L = with = -yκ 0 0 0 .

(2.7)

Eq. (2.7) is formally identical to the expression used in the linear beam theory. The crucial difference is that the curvature κ now depends on the deflection v in a nonlinear fashion, as can be seen from Eq. (2.5).

Equation of motion

The equation of motion can be obtained using Hamilton's principle

δ T 0 L + W e = 0 (2.8)
where L is the lagrangian and W e is the external work. The relation (2.8) holds on any time interval [0, T ]. In the present case, L takes the form

L = U -W + 1 2 L 0 λ(x)(1 -(1 + u ) 2 + v 2 ) 2 dx
where U is the kinetic energy, W is the elastic energy, and λ(x) is a Lagrange multiplier associated with the inextensibility constraint. To proceed further, one has to specify the expressions of U , W and W e . The elastic energy W is taken in the form

W = 1 2 L 0 EI(x)κ 2 (x) dx (2.9)
where EI(x) is the flexural rigidity of the beam (not necessarily a constant). The kinetic energy U is written as

U = 1 2 L 0 µ(x)(( u + uB ) 2 + ( v + vB ) 2 ) dx (2.10)
where µ(x) is the linear mass density and (u B , v B ) is the base displacement. Substituting in (2.8) and omitting the detail of the calculation, the dynamic equation (in the case W e = 0) is found to be

µ(v + vB ) + EIv = -EIv v v -(EI) v v 2 - 1 2 v s L µ ∂ 2 ∂t 2 s 0 v 2 ds ds + v s L µds + µv üB (2.11)
which is complemented by the boundary conditions

v(0, t) = v (0, t) = 0, v (L, t) = (EIv (1 + v 2 )) (L, t) = 0.
(2.12)

Single mode approximation

In this paper, we are interested in deriving a simple model for the dynamics of the beam. This can be accomplished by using the single mode approximation, i.e. by writing v(x, t) as v(x, t) = φ(x)r(t) (2.13) where φ is the linear mode shape whose associated natural frequency is the closest to the main excitation frequency. Provided the mode shape φ is normalized in such fashion that φ(L) = 1, the variable r(t) can be interpreted as the tip deflection of the beam. Using (2.13), the kinetic energy and the elastic energy become

U = 1 2 L 0 µ ( u + uB ) 2 + (φ ṙ + vB ) 2 dx , W = 1 2 L 0 EIφ 2 r 2 1 + φ 2 r 2 dx.
To account for energy dissipation in the system (notably due to friction with the air), the external work δW e is taken as

δW e = -(c + d| ṙ|) ṙδr (2.14)
where c and d are positive parameters. The introduction of a quadratic damping coefficient d results from the fact that the deflections can be large. It is considered that the linear and quadratic damping parameters account for all the sources of energy dissipation, i.e. friction forces (with the air and the support) as well as material damping. Substituting in (2.8) and omitting the detail of the calculations, we obtain the dynamic equation (2.16)

In (2.15), the transverse excitation vB acts as a forcing term (direct excitation) whereas the axial excitation üB acts as a time-dependent variation of the stiffness (parametric excitation). The axial excitation üB appears in (2.15) as a consequence of the inextensibily condition, which relates a variation δv of the transverse displacement to a variation δu of the axial displacement. From (2.6) and (2.13) we have indeed δu = -v δv = -δr(t)r(t)φ 2 . Now when writing down Hamilon's principle, the term δW coming from the kinetic energy gives rise to the integral µ(x) uB δ udx, which after integration by part (both in space and time) leads to the term ξ üB r(t) that appears in (2.15).

Remark on the Galerkin approximation: One can use a more general approximation of v(x, t) in the form

v(x, t) = n i=1 φ i (x)r i (t)
where φ 1 ,• • • ,φ n are given mode shapes. In that case, Eq. (2.15) is replaced by a set of n coupled nonlinear equations as detailed in Appendix A.

Remark on the assumption of inextensibility: The equations (2.11) and (2.15) have been obtained using the condition of inextensibility of the beam. To justify that assumption, it is interesting to compare the first natural frequencies of transverse and axial vibrations. Considering for simplicity that the flexural rigidity EI, the mass density µ and the area A of the cross section are all independent on x, it is well known that the first natural frequency ω 0 of transverse vibrations is given by

ω 0 = k v L 2 EI µ (2.17)
and that the first natural frequency ω 0 of axial vibrations is given by

ω 0 = π 2L EA µ .
In (2.17), k v is the first positive solution to the equation cos

k v cosh k v = -1. Numerical calculations show that k v ≈ 1.875. It follows that ω 0 ω 0 2 ∼ L 2 A I ∼ L h 2 1
where h is the length of the beam in the y direction. Therefore, if the excitation frequency ω is of the order of ω 0 -which will be the case in this paper -then it makes sense to neglect axial vibrations and to consider that the beam is inextensible.

Piezoelectric behavior at large strain 2.2.1. Constitutive laws

We now detail the behavior of the piezoelectric. In accordance with common practice, the behavior of the piezoelectric is written in a local basis (e 1 , e 2 , e 3 ) where e 3 is the poling direction. For the patch located on the top surface y = h 2 , the poling direction is set equal to e y . For the patch located on the bottom surface y = -h 2 , the poling direction is set equal to -e y (so as to avoid cancellation effects when the patches are connected in series). Writing the strain in the local basis (e 1 , e 2 , e 3 ), the kinematics of the beam imply that 11 is the only non vanishing component. We furthermore assume that the electric field E and the electric displacement D are along the e 3 direction, i.e. E = E 3 e 3 and D = D 3 e 3 . In such conditions, the free energy density H of the piezoelectric can be regarded as a function of ( 11 , E 3 ). The constitutive laws of the piezoelectric are given by

σ 11 = ∂H ∂ 11 , D 3 = - ∂H ∂E 3 . (2.18)
The enthalpy density H is commonly taken as quadratic in ( 11 , E 3 ), in which case the constitutive relations (2.18) are linear. Such a choice is valid provided that 11 and E 3 remain small enough compared to some characteristic values of the material, henceforth denoted by 0 and E 0 . To account for large deflections, we need an expression of H holding for small -but not infinitesimally small -values of 11 / 0 . The dimensionless electric field E 3 /E 0 is assumed to remain small compared to 11 / 0 . Writing H as a function of the dimensionless variables ( 11 / 0 , E 3 /E 0 ), it is natural to consider a Taylor expansion of H, i.e. H(

11 0 , E 3 E 0 ) 1≤i+j≤n c ij ( 11
0 ) i ( E 3 E 0 ) j . (2.19)
where n is the order of the expansion. Assuming that E 3 /E 0 is of the order of ( 11 / 0 ) 2 and limiting the expansion to the order n = 4, we get an expression of the form

H( 11 , E 3 ) = 1 2 E p 2
11 -e 31 E 3 11 -

1 2 ε S 33 E 2 3 + 1 6 c 0 2 11 E 3 + 1 6 c 1 3 11 + 1 24 c 2 4
11 .

(2.20)

In (2.20), E p is the Young's modulus at infinitesimal strains, e 31 is the linear piezoelectric constant, and ε S 33 is the dielectric permittivity at constant strain. The additional parameters c 0 , c 1 and c 2 introduce nonlinearities in the constitutive laws. In particular, for short circuit conditions (E 3 = 0), the stress-strain relation obtained from (2.20) Thus the parameter c 1 (resp. c 2 ) models a linear (resp. quadratic) dependence of the secant modulus with respect to the applied strain. In particular, the parameter c 2 drives the stress-strain behavior at high strain: If c 2 < 0 then the stress as well as the secant modulus decrease for sufficiently high strain, i.e. the material exhibits stress softening. Conversely, if c 2 > 0 then the stress and the secant modulus increase for sufficiently high strain. In that case, the material exhibits stress hardening.

The parameter c 0 in (2.20) can be interpreted as a nonlinear piezoelectric coupling parameter. In the following, we make the simplifying assumption that c 0 = 0, i.e. we only take purely elastic nonlinearities of the piezoelectric into account. Some justification for dropping the parameter c 0 is provided by the results of Abdelkefi et al. [START_REF] Abdelkefi | Global nonlinear distributed-parameter model of parametrically excited piezoelectric energy harvesters[END_REF][START_REF] Abdelkefi | Effects of nonlinear piezoelectric coupling on energy harvesters under direct excitation[END_REF]. Those authors considered an expansion of the form (2.19) up to the order n = 3, retaining all the terms in the expansion and using numerical simulations to investigate the influence of the constitutive parameters on the frequency response of the piezoelectric. In particular, the primary effect of the parameter c 0 in (2.20) was found to reduce the amplitude of the response. In practice it turns out to be difficult to distinguish that effect from the damping contribution (2.14). In contrast, elastic nonlinearities will be shown in Sect. 3.2-4 to have a dramatic influence on the frequency response of the system (notably on the resonance frequency). Our objective being to keep the model as simple as possible while still capturing the main physics of the problem, we choose to set c 0 = 0. Ultimately the validity of that assumption will be validated in Sect.4 by comparison with the experiments.

Equation of motion for the piezoelectric-equipped beam

The equation of motion of the piezoelectric-equipped beam can be obtained using Hamilton's principle (2.8) where the lagrangian L is now given by

L = U -W -W p + 1 2 L 0 λ(x)(1 -(1 + u ) 2 + v 2 ) 2 dx. (2.23)
In (2.23), W p is the elastic energy of the piezoelectric. We have

δW p = Ωp ∂H ∂ 11 δ 11 dΩ
where Ω p is the domain occupied by the piezoelectric. Assuming the electric field to be constant in Ω p and provided that the patches are connected in series, we have

E 3 = -V (t) 2hp
where V (t) is the voltage across the piezoelectric. We note from (2.7) and (2.13) that

11 = -y r(t)φ + 1 2 r(t) 3 φ φ 2 (2.24)
Hence

δ 11 = -y φ + 3 2 r(t) 2 φ φ 2 δr.
Substituting in (2.23) and using the expression (2.20) for H, we obtain (omitting the details of the calculations)

δW p = Ωp E p 11 δ 11 dΩ -(θ + N r 2 )V δr + α p r 3 δr (2.25)
where

θ = χ Lp 0 φ (x)dx, N = 3 2 χ Lp 0 φ φ 2 dx, χ = 1 2 e 31 b(h + h p ), α p = 1 5 c 2 b([ h 2 + h p ] 5 -[ h 2 ] 5 ) Lp 0 φ 4 dx. (2.26)
In (2.26), b is the width of the patches.

Observe that the parameter c 1 does not appear in (2.25): This is due to the symmetric arrangement of the two patches. The term Ω E p 11 δ 11 dΩ is the linear contribution to the mechanical energy of the piezoelectric and can be put in the form (2.9). Using the expression (2.26), Hamilton's principle leads to the dynamic equation

M r + (K + ξ üB ) r + β rr + ṙ2 r + αr 3 + c ṙ + d ṙ| ṙ| -θV -N r 2 V = -mv B (2.27)
where it is assumed that the flexural rigidity EI(x) and the linear mass density µ(x) take the contribution of the piezoelectric into account. In (2.27), the parameter α is defined by

α = α b + α p .
(2.28)

The parameter α characterizes the cubic nonlinearity of the system. As can be seen on (2.28), α results from a mechanical contribution α b and a piezoelectric contribution α p . Whereas α b is always positive (see Eq. (2.16)), α p can either be positive or negative, depending on the sign of c 2 . Eq. (2.27 ) has been obtained under the assumption that the substrate material is linear. This is reasonable if the strain is not too large (below 0.1% for a steel substrate). In other words, piezoelectric nonlinearities are assumed to appear at the smaller strain level than the substrate nonlinearities.

Electrical circuit

v( h 2 + hp) piezoelectric v( h 2 ) v(-h 2 -hp) piezoelectric v(-h 2 ) R i Figure 3: Electrical circuit.
The equation (2.27) is complemented by the electrical equation of the circuit (Fig. 3). The piezoelectric patches are connected in series to a load resistance R. We denote by

v( h 2 + h p ) (resp. v( h 2 )) the voltage of the electrode at y = h 2 + h p (resp. y = h 2 ). Similarly we denote by v(-h 2 -h p ) (resp. v(-h 2 )) the voltage of the electrode at y = -h 2 -h p (resp. y = -h 2 ). We have V R = i where V = v( h 2 + h p ) -v(-h 2 -h p )
is the voltage across the resistance load. The current intensity i in the circuit is estimated as the surface integral of Ḋ3 on the mid-surface of a piezoelectric layer, i.e.

i = y= h+hp 2 Ḋ3 dxdz.
From the constitutive relation (2.18) and the expression (2.24) of the strain, we find

Ḋ3 = -ye 31 φ + 3 2 r(t) 2 φ φ 2 ṙ + ε S 33 Ė3 . (2.29) which after integration gives i = -θ ṙ -N r 2 ṙ -2C p v( h 2 + h p ) -v( h 2 ) with C p = ε S 33 bL p 2h p . ( 2

.30)

Conducting a similar reasoning for the piezoelectric on the bottom surface leads to

i = -θ ṙ -N r 2 ṙ -2C p v(- h 2 ) -v(- h 2 -h p ) . Since v( h 2 ) = v(-h
2 ) (see Fig. 3), we finally arrive at

C p V + V R + θ ṙ + N r 2 ṙ = 0. (2.

31)

Remark: Piezoelectric patches are often manufactured as bimorph, i.e. with two piezoelectric layers separated by an epoxy layer. The analysis detailed above can easily be extended to such configurations.

Coupled dynamical system

The two equations (2.27) and (2.31) govern the evolution of the variables (r, V ). In the following, we study the dynamics of the system under a harmonic base acceleration, in either direct or parametric excitation. For a direct excitation (u B = 0) with harmonic base acceleration (v B = B cos ωt), we obtain the dynamical system

M r + Kr + β rr + ṙ2 r + αr 3 + c ṙ + d ṙ| ṙ| -θV -N r 2 V = mBω 2 cos ωt, C p V + V R + θ ṙ + N r 2 ṙ = 0. (2.32)
For a parametric excitation (v B = 0) with harmonic base acceleration (u B = B cos ωt), we obtain the dynamical system

M r + K -ξBω 2 cos ωt r + β rr + ṙ2 r + αr 3 + c ṙ + d ṙ| ṙ| -θV -N r 2 V = 0, C p V + V R + θ ṙ + N r 2 ṙ = 0.
(2.33)

The linear version of (2.32) -obtained by setting α = d = N = 0 -is commonly used in studies related to piezoelectric energy harvesting. In such case, the steady state harmonic response of the system is of the form

r = Re(a e iωt ) , V = Re(b e iωt ) with b = -iRθωa/(1 + iRC p ω) and a = mω 2 B M (ω 2 0 -ω 2 ) + icω + iRθ 2 ω/(1 + iRC p ω)
.

Let P be the average value of the electrical power V (t) 2 /R delivered to the load resistance.

We have

P = 1 2 Rθ 2 ω 2 1 + R 2 C 2 p ω 2 a 2 .
(2.34)

At the first order in (c, θ 2 ), it can be verified that the resonance frequency ω is given by

ω = ω 0 + λ (2.35)
where

ω 0 = K/M and λ = R 2 C p θ 2 ω 0 2M (1 + R 2 C 2 p ω 2 0 )
.

Observe from (2.35) that the piezoelectric coupling results in a shift λ of the resonance frequency compared to the natural frequency ω 0 of the beam.

Crucial differences are found in the case of parametric excitation (2.33) compared to the case of direct excitation. In contrast with (2.32), r = V = 0 is always a solution of (2.33): The piezoelectric beam may remain at equilibrium even if excited. If, however, the amplitude B of the excitation is sufficiently large, then the equilibrium state r = V = 0 becomes unstable so that any small disturbance leads to oscillations with possibly large amplitudes. The limit of stability B s (ω) (i.e. the maximum amplitude B for which the equilibrium r = V = 0 is stable) can be estimated using Floquet theory, as detailed in Appendix B. In particular, for the undamped beam without piezoelectric coupling (θ = N = c = d = 0), B s (ω) vanishes for ω of the form

ω n = 2ω 0 n
where n ∈ N. The frequency ω n is called the n th parametric resonance. Those results are modified in presence of piezoelectric coupling: At the first order in (N, c, θ 2 , ω -2ω 0 ), the limit of stability of the piezoelectric beam is given by

4B 2 s (ω)ω 2 0 ξ 2 = c + Rθ 2 1 + R 2 C 2 p ω 2 0 2 + R 2 C p θ 2 ω 0 1 + R 2 C 2 p ω 2 0 -M (ω -2ω 0 ) 2 .
(2.36)

We refer to Appendix B for a derivation of (2.36) using Floquet theory. Observe that B s (ω) is minimum for

ω = 2(ω 0 + λ) (2.37)
Hence, in a similar way to direct excitation, piezoelectric coupling results in a shift of the parametric resonance. Observe that the relative shift (with respect to the relevant natural resonance frequency) is the same as for direct excitation and equal to λ/ω 0 . Also note that B s (ω) increases with the linear damping coefficient c: the more damping is present in the system, the higher the excitation needs to be for oscillations to occur.

First-order expansion of the nonlinear response

Steady-state solution

We are interested in studying the steady-state response of the system for an excitation frequency close to the main resonance ω R (i.e. ω R = ω 0 for direct excitation, and ω R = 2ω 0 for parametric excitation). The systems (2.32) and (2.33) are nonlinear and cannot be solved exactly. However, using perturbation methods, we can obtain a closed-form expansion of the solution in the case where the nonlinear terms as well as the excitation amplitude B are considered as small. More precisely, the parameters α, c, d, β, B, θ, N in (2.32)-(2.33) are assumed to vary as

B = O(ε 2 ), c = O(ε 2 ), d = O(ε 2 ), α = O(ε 2 ), β = O(ε 2 ), N = O(ε 2 ), θ = O(ε)
where ε is a small dimensionless parameter driving the expansion. Although they do not correspond to nonlinear terms, the linear damping coefficient c and the piezoelectric coupling parameter θ need to be considered as small quantities for the analysis to be tractable. Observe, however, that θ is allowed to remain one order of magnitude larger than the other parameters. The excitation frequency ω is written as

ω = ω R + σ
where the frequency shift σ is assumed to be of the order 2 in ε.

Following the method of multiple time scales [START_REF] Nayfeh | Nonlinear Oscillations, wiley classics library Edition[END_REF][START_REF] Nayfeh | Introduction to Perturbation Techniques[END_REF], the functions r(t) and V (t) are expanded as

r = r 0 (T 0 , T 1 , T 2 ) + εr 1 (T 0 , T 1 , T 2 ) + ε 2 r 2 (T 0 , T 1 , T 2 ) + O(ε 3 ), V = V 0 (T 0 , T 1 , T 2 ) + εV 1 (T 0 , T 1 , T 2 ) + ε 2 V 2 (T 0 , T 1 , T 2 ) + O(ε 3 ) (3.1)
with T i = ε i t. The time scale T 0 is the physical time t, and the time scale T i+1 is a slow time scale compared to T i . The expressions of r i and V i in (3.1) are obtained by expanding (2.32) and (2.33) in power of ε and collecting the terms with the same order. The analysis is similar to that performed by Abdelkefi et al. [START_REF] Abdelkefi | Global nonlinear distributed-parameter model of parametrically excited piezoelectric energy harvesters[END_REF][START_REF] Abdelkefi | Effects of nonlinear piezoelectric coupling on energy harvesters under direct excitation[END_REF], the main difference being in the definition of the various constants that appear in (2.32) and (2.33).

In the case of direct excitations, solutions of (2.33) are found to be of the form (in the first order in ε)

r = a cos (ωt -γ) , V = ω 0 θR 1 + R 2 C 2 p ω 2 0 a cos (ωt -γ + γ V ) (3.2)
where (a, γ) are functions of T 2 that satisfy the modulation equations

ω 0 M da dT 2 = -ω 0 c eff 2 a - 4ω 2 0 d 3π a 2 - mBω 2 0 2 sin γ, ω 0 M dγ dT 2 = ω 0 M (σ -λ) - 1 2 κ eff a 2 - mBω 2 0 2a cos γ. (3.3)
In (3.3), the parameters c eff and κ eff are defined by

c eff = c + Rθ 2 1 + R 2 C 2 p ω 2 0 , κ eff = 3α -2βω 2 0 4 . (3.4)
The parameter c eff can be interpreted as an effective linear damping parameter, that combines mechanical and electrical damping (represented by the terms c and Rθ 2 /1 + R 2 C 2 p ω 2 0 , respectively). The parameter κ eff can be interpreted as an effective nonlinear parameter.

Eq. (3.2) shows that r 0 is a slowly modulated harmonic function of pulsation ω 0 . In the expression of V in (3.7), there appears a additional phase shift γ V whose exact expression is not needed for our purpose. Steady-state solutions are obtained by enforcing that da/dT 2 = dγ/dT 2 = 0 in (3.3), i.e.

0 = -ω 0 c eff 2 a - 4ω 2 0 d 3π a 2 - mω 2 0 B 2 sin γ, 0 = ω 0 M (σ -λ) - 1 2 κ eff a 2 - mω 2 0 B 2a cos γ. (3.5)
The amplitude a (as well as the phase shift γ) depends nonlinearly on the excitation amplitude B and the frequency shift σ.

A similar analysis can be performed for parametric excitations. The modulation equations become

ω 0 M da dT 2 = -ω 0 c eff 2 a - 4ω 2 0 d 3π a 2 -ξω 2 0 Ba sin γ, ω 0 M dγ dT 2 = ω 0 M (σ -2λ) -κ eff a 2 -2ξω 2 0 B cos γ. (3.6)
In that case, steady-state solutions are found to be of the form

r = a cos ω 2 t - γ 2 , V = ω 0 θR 1 + R 2 C 2 p ω 2 0 a cos ω 2 t - γ 2 + γ V (3.7)
where (a, γ) satisfy

0 = -ω 0 c eff 2 - 4ω 2 0 d 3π a -ξω 2 0 B sin γ, 0 = ω 0 M (σ -2λ) -κ eff a 2 -2ξω 2 0 B cos γ. (3.8)
Observe from (3.7) that the system exhibits a harmonic response at half the frequency of the excitation.

Frequency-response

The expressions (3.5-3.8) can be used to study the influence of the constitutive parameters on the steady-state frequency response of the system. To that purpose, it is useful to express the various parameters in terms of dimensionless quantities (denoted by a tilde ˜) as follows:

c eff = ceff √ KM , d = dM/L , ξ = ξM/L , κ eff = κeff K/L 2 , a = Lã , B = L B.
It is also convenient to express the excitation frequency ω in terms of the relative frequency shift σ with respect to the linear resonance frequency, i.e.

σ = ω -ω R ω R - 1 2 R 2 C 2 p θ 2 M (1 + R 2 C 2 p ω 2 0 ) . (3.9)
For energy harvesting, the electric power delivered to the load resistance is the main quantity of interest. Since the instantaneous power is equal to V 2 (t)/R, it is clear from the expression (3.2-3.7) that the average power P delivered to the load resistance is

P = P 0 ã2 (3.10)
where ã corresponds to a steady-state solution and

P 0 = L 2 2 Rθ 2 ω 2 0 1 + R 2 C 2 p ω 2 0 .
Eliminating the phase γ from (3.5)-(3.8) leads to the equation

m M 2 B2 4ã 2 = ( 1 2 ceff + 4 d 3π ã) 2 + (σ - 1 2 κeff ã2 ) 2 for direct excitation (3.11)
and to the equation

ξ2 B2 = ( 1 2 ceff + 4 d 3π ã) 2 + (σ - 1 2 κeff ã2 ) 2
for parametric excitation.

(3.12)

In (3.11), the constitutive parameters only appear through the effective linear damping ceff , the quadratic damping parameter d and the effective nonlinear parameter κeff . In (3.12), the parametric coefficient ξ acts as an additional parameter.

In order to study the influence of the parameters (c eff , d, κeff ) on the power P , the dimensionless frequency responses obtained from (3.10-3.11) are plotted in Figs 4-6 both for direct (left) and parametric (right) excitations. In the parametric case, it can be observed in Figs 4-6 that P vanishes at two frequencies. It can easily be checked from (3.12) that those two frequencies are given by the relation

ξ2 B2 = 1 4 c2 eff + σ2 . (3.13)
Note that the relation (3.13) corresponds to the limit of stability of the equilibrium position, as obtained in Eq. (2.36) from Floquet theory.

Increasing the damping parameters ceff and d has the effect of decreasing the amplitude a and the electric power P (Figs. 45). In the case of parametric excitation, the linear damping ceff has the additional effect of reducing the frequency range over which oscillations occur. In contrast, the quadratic damping parameter d does not display such an influence. Whereas the damping parameters primarily act on the amplitude of the frequencyresponse, the nonlinear parameter κeff primarily affects the resonance frequency: Increasing κeff results in an increase of the resonance frequency, as shown in Fig. 6. Also observe the foldover effect (i.e. the bending of the resonance curve) that is typical of nonlinear oscillators. If κeff > 0, the foldover effect is of the hardening type, i.e. the resonance curve is bent towards high frequencies. Conversely, if κeff < 0 then the foldover effect is of the softening type, i.e. the resonance curve is bent towards low frequencies.

For later reference, we note that the expression of the peak power P M and the corresponding frequency shift σM can be obtained from (3.11-3.12). We have indeed

P M = P 0 ã2 M (3.14)
where ãM is the peak displacement amplitude. We obtain from (3.11)-(3.12) that σM = 1 2 κeff ã2 A piezoelectric cantilever has been designed in order to study the validity of the model presented in Sect. 2-3. The device is shown in Fig. 7: Two Midé QP20W bimorph piezoelectric patches were bonded to the upper and lower surfaces of a 222 x 0.8 x 40 mm 3 steel plate. The piezoelectric patches were bonded at the clamped end side of the plate. A 12 g tip mass was placed on the steel surface for tuning the resonance frequency and amplifying the amplitudes of vibrations. In Table 1 are reported the piezoelectric material parameters provided in the manufacturer's datasheet (E p and E e are respectively the Young's modulus of the piezoelectric material and of the epoxy, ρ p and ρ e are the mass densities, e 31 is the piezoelectric constant, ε S 33 is the relative permittivity at constant strain). Table 1: Material parameters for the piezoelectric patch.

M (3.15) where ãM =          3π 16 d   -c eff + c2 eff + 32 d 3π m M | B|   for direct excitation,
The experimental set-up is shown in Fig. 8: An electromagnetic shaker was used to apply prescribed harmonic excitations to the piezoelectric cantilever. Direct excitation was achieved by placing the cantilever horizontally (Fig. 9) while parametric excitation was achieved by placing the cantilever vertically (Fig. 10). The piezoelectric patches were connected in series, with a resistive load R = 28 kΩ. The voltage across the resistive load as well as the base acceleration were monitored. The particular value R = 28 kΩ was chosen because it was found experimentally to be reasonably close to the optimal resistance (i.e. the resistance that maximizes the electrical power) over all the excitations tested (both direct and parametric).

Resonance frequencies in the linear regime

In order to evaluate the resonance frequencies for direct excitations in the linear regime, a spectral analysis was performed on signals measured from free vibrations. Two resonance frequencies were found in the range 0 -150 Hz: the first resonance frequency ω 0 /2π is approximatively equal to 15.14 Hz, and the second resonance frequency is approximatively equal to 81.1 Hz (Fig. 11). As a comparison, the first two resonance frequencies obtained from a Finite Element Analysis (FEA) of the beam (carried out with the FEM software CAST3M [START_REF]CEA, Cast3m homepage[END_REF]) are approximatively equal to 15.17 Hz and 109.9 Hz. The corresponding mode shapes φ 1 and φ 2 are represented in Fig. 12. A finite-element beam model with piecewise-continuous properties was used so as to take the geometries of the piezoelectric patches and of the distributed proof mass into account. Regarding the first resonance ω 0 , the relative error between the experimental and numerical values is about 0.2%. For the second resonance, the relative error is about 34 %. Such a large value probably comes from the fact that the assumption of a perfect clamping -used in the FEA -is no longer valid for high vibration modes. For our purpose, the main observation is that the gap between the first two resonance frequencies is large, which justifies taking φ = φ 1 in the single mode calculated from the mode shape φ 0 in Fig. 12 and the material parameters in Table 1. The obtained values are reported in Table 2. At this point the parameters c, d and α p remain to be identified.

M K ξ β α b θ N C p m (g) (N.m -1 ) (kg.m -1 ) (kg.m -2 ) (N.m -3 ) (N.V -1 ) (N.m -2 .V -1 ) (nF) (g

Identification of the linear damping parameter c

The linear damping parameter c can be obtained from the small amplitude decay in free vibrations. For small amplitudes of the displacement, all the nonlinear terms in (2.32-2.33) indeed become negligible so that the voltage response V (t) in free vibrations takes the form

V (t) = ω 0 θR 1 + R 2 C 2 p ω 2 0 V 0 e c eff t/2M cos (ω 0 t -γ + γ V ) (4.1)
where V 0 is a constant. Eq. (4.1) can notably be obtained by using (3.2) and neglecting the quadratic term in (3.3). It follows from Eq. (4.1) that the ratio between the initial voltage V (0) and the n th peak of the voltage response is

V n V (0) = e -nπc eff /2M ω 0 . (4.2)
In Fig. 13(left) are shown some measurements of the voltage response V (t) in free vibrations at small amplitudes. Those measurements have been obtained with a resistive load R = 28 kΩ. The corresponding values of log V (0)/V n are shown in Fig. 13(right) and display a linear dependence with n. Using (4.2), the slope of that line gives the value of c eff , from which the linear damping c can be deduced by using Eq. (3.4). Carrying out such a procedure gives c = 0.034 N.s.m -1 .

Identification of the quadratic damping d and of the nonlinear piezoelectric parameter

α p Using the formula established in Sect. 3.2, the nonlinear parameters (d, α p ) can be obtained from the measurement of the peak power P M and the corresponding frequency shift σM at any given amplitude of excitation (either direct or parametric). In the case of a parametric excitation, we have indeed from (3.16)

d = 3π 4L |ξB| - c eff 2ω 0 P 0 P M , κ eff = 2K L 2 P 0 P M σM 2 . (4.3)
For a parametric excitation with acceleration amplitude U 0 = Bω 2 = 5.83 g, a peak power P M of 0.9605 mW was measured (see Fig. 14(a)). The corresponding frequency shift σM is equal to 6.6 10 -3 . Using those measurements, the relations (4.3) give

d = 0.018 kg.m -1 , κ eff = -2.29 10 4 N.m -3
In the identication procedure that has been just described, the parameters d and κ eff are found simultaneously. Since κ eff is directly related to the stress-strain nonlinearity of the piezoelectric, an other possible approach for identifying κ eff is to use (quasi-static) tests in short-circuit conditions. At this point all the model parameters have been found. Out of the 12 parameters of the model, 9 are directly available from a finite element analysis (Table 2). The remaining 3 are deduced from experimental measurements as detailed above (Table 3). Table 3: Model parameters (identified from experiments).

Frequency response in the nonlinear regime

The red dots in Fig. 14 show the average power P delivered to the load resistance for several amplitudes and frequencies of the excitation. Such measurements have been performed both for direct and parametric excitations. For parametric excitations, a foldover effect of the softening type can clearly be observed on Fig. 14(a)(c). A similar effect is also visible for direct resonance at sufficiently high excitation amplitude, see Fig. 14(b). Note that a softening behaviour under direct excitation has also been observed experimentally in [START_REF] Stanton | Nonlinear nonconservative behavior and modeling of piezoelectric energy harvesters including proof mass effects[END_REF][START_REF] Priya | High-power resonant measurements of piezoelectric materials: importance of elastic nonlinearities[END_REF]. From the definition (3.4) of κ eff , we can deduce that α p = -3.09 104 N.m -3 . Finding a negative value for κ eff is consistent with the parametric study detailed in Sect. 3.2: Softening resonance curves, as observed experimentally, correspond to negative values of the nonlinearity parameter κ eff . Recall from (2.28) and (3.4) that

κ eff = 1 4 (3α b -2βω 2 0 ) +
Using the values in Table 2, the purely mechanical contribution 1 4 (3α b -2βω 2 0 ) is found to positive. Hence the introduction of the nonlinear piezoelectric term α p is crucial to capture the softening behavior that is observed in the frequency response. Note from (2.26) that α p is directly related to the piezoelectric parameter c 2 which models the quadratic dependence of the Young modulus with respect to the applied strain. Incidentally, finding α p < 0 means that c 2 < 0, i.e. that the piezoelectric material exhibits stress softening at high strain. A similar result (i.e c 2 < 0) was obtained in [START_REF] Von Wagner | Piezo-beam systems subjected to weak electric field: experiments and modelling of non-linearities[END_REF]. Figure 14: Frequency responses for several acceleration amplitudes U 0 (g = 9.81 m.s -2 ): experimental results (red dots), analytical expressions (solid blue line), numerical simulations (pink solid line).

Model validation

The solid blue lines in Fig. 14 shows the frequency response predicted by the expressions (3.10-3.11-3.12). A very good agreement with the experimental results is observed over all the range of tested loadings, both in direct and parametric excitations. We emphasize that the same set of parameters (Tables 23) is kept for all the curves in Fig. 14. Recall that the resonance frequency ω M and peak power P M in Fig. 14(a) are the only experimental measurements used for identification purpose. The very good agreement between the model and the experiments thus indicate that the model gives a fairly accurate description of the physical system, at least close to the resonance and for moderately large amplitude of excitations.

The pink solid curves in Fig. 14 show the frequency response obtained by numerical integration of the dynamical systems (2.32) and (2.33) using a Runge-Kutta scheme. More precisely, starting from the initial condition (r(0), V (0)) = (0.0001, 0), the time integration is performed on a time sufficiently long for the response to stabilize on a steady-state cycle (60 s was found to be sufficient in all the considered cases). The results obtained in such fashion are very close to the results provided by the closed-form expressions (3.10-3.11-3.12), thus validating the analytical approach considered in Sec. 3. For parametric excitations, the sudden 'vertical' jump displayed by the numerical resonance curve is related to instability phenomena, as detailed later in Sect. 5.3.

The results in Fig. 14 rely on the single mode approximation ( 14) of the dynamics. As discussed in detail by Abdelkefi et al. [START_REF] Abdelkefi | Global nonlinear distributed-parameter model of parametrically excited piezoelectric energy harvesters[END_REF][START_REF] Abdelkefi | Effects of nonlinear piezoelectric coupling on energy harvesters under direct excitation[END_REF], the number of required modes in the Galerkin projection is strongly influenced by the value of the excitation acceleration. In order to assess the domain of validity of the single mode approximation for the problem at hand, numerical simulations have been performed using up to 3 modes in the Galerkin projection. The corresponding equations are reported in Appendix A, see Eqs. (A.2-A.3). The same linear and quadratic damping have been used for all modes. Some results are reported in Fig. 15 (for parametric excitations) and in Fig. 16 (for direct excitations). Both for direct and parametric excitations, it has been found that the 1-, 2-and 3-modes approximations all give similar results up to some value U max 0 of the amplitude of excitations. Above that value, significant differences appear between the different approximations considered so that the validity of the single mode approximation becomes questionable, see Figs 15(b)-16(b). The limit value U max 0 is about 0.22g for direct excitations and 5.9g for parametric excitations: Above those values, the relative difference between the peak power provided by the 1-, 2-and 3-modes approximations is significant and exceeds 10%.

Applications in energy harvesting

In this Section we further explore some consequences of the proposed model for applications in energy harvesting. Except stated otherwise, all the numerical curves presented in the following have been obtained using the model parameters in Tables 2-3.

Optimal resistance

For energy harvesting, a major concern is to estimate the maximum electrical power that can be extracted from a given excitation. This is achieved by assuming the perfect tuning condition ω = ω M and optimizing the load resistance R. To that purpose, it is convenient to introduce the dimensionless quantities defined as

R = RC p ω 0 , θ2 = θ 2 C p K , c = c √ KM .
From the relation (3.10), we obtain that the optimal resistance necessarily satisfies

0 = ãM ∂P 0 ∂R + 2P 0 ∂ã M ∂R i.e. 0 = (1 - R2 )ã M + 2 R(1 + R2 ) ∂ã M ∂R . (5.1)
At this point it is necessary to distinguish between direct and parametric excitations because the expression of ãM is not the same in both cases, as can be seen from Eq. (3.16). The case of parametric excitation is the easiest one to handle. In that case, (5.1) becomes

0 = (1 - R2 )(1 - 3 θ2 2| B ξ| - c R + R2 ). (5.2)
The value R = 1 (i.e. R = 1/C p ω 0 ) is always a solution of (5.2) and corresponds to a stationary point of P (R) (but not necessarily to a maximum). Depending on the excitation amplitude B, other solutions to (5.2) may exist. After some manipulations on (5.2), we obtain that:

• if | B ξ| -c/2 < 3 2
θ2 , then there are two values of the optimal resistance obtained as the solutions to the quadratic equation 0 = 1 -3 θ2 2| B ξ|-c R + R2 . The corresponding value of the electric power is

P = π 2 32 d2 L 2 Kω 0 (2| B ξ| -c) 3 3 . • if | B ξ| -c/2 > 3 2
θ2 , then the optimal resistance is given by R = 1. The corresponding electric power is

P = π 2 32 d2 L 2 Kω 0 9 θ2 (2| B ξ| -c -1 2 θ2 ) 2 8 .
For the device considered in Sec. 4, the threshold value B separating the two regimes is equal to 0.011. In Fig. 17 is plotted the optimal resistance as a function of the excitation amplitude.

Although similar in spirit, the calculations for direct excitations are more involved. Substituting the relevant expression of ãM in (5.1) leads to the equation

0 = (1 -R2 ) (1 + R2 ) c2 eff + 32 d 3π m M | B| -2 Rθ 2 .
Again R = 1 is always a solution. Depending on the excitation amplitude B, other solution may be found by solving the equation Recall that ãM is required to be positive, so (5.4) implies that

0 = (1 + R2 ) c2 eff + 32 d 3π m M | B| -2 Rθ 2 . ( 5
0 > R2 0 -R0 θ2 c + 1. ( 5.5)
Now the minimum value of the polynomial x → x 2 -x θ2 c + 1 is equal to 1 -θ4 /4c 2 . Hence (5.5) implies that θ2 ≥ 2c.

(5.6)

If the model parameters do not satisfy (5.6) -which happens to be the case for the device considered in Sect. 4 -then (5.3) cannot have any admissible solution. In such situation, the optimal resistance is given by R = 1 and does not depend on the excitation amplitude.

Limits of the linear model

As mentioned in Sect. 2.3, linear models of piezoelectric cantilever under direct excitation are frequently used in studies related to energy harvesting. In such models, the peak power and the resonance frequency are given by (2.34) and (2.35), respectively. The nonlinear model presented in this paper allows one to assess the limits of validity of the linear model: In Fig. 18(left) is plotted the relative error between the linearized expression (2.34) and the expression (3.15-3.16) of the peak power P M , as a function of the acceleration amplitude U 0 . The relative error remains small (<0.2%) for acceleration U 0 below 0.01 m.s -2 , but grows fast for acceleration amplitude above 0.01 m.s -2 . The relative error exceeds 10 % for acceleration amplitude above 0.16 m.s -2 . Note that acceleration amplitudes above 0.16 m.s -2 are common in low-level ambient vibration sources [START_REF] Du Toit | Design considerations for mems-scale piezoelectric mechanical vibration energy harvesters[END_REF], meaning that nonlinear effects are expected to be significant in a lot of applications related to energy harvesting. As illustrated in Fig. 18(right), neglecting nonlinearities leads to overestimate the harvested power. A similar conclusion has been obtained in [START_REF] Stanton | Nonlinear piezoelectricity in electroelastic energy harvesters: Modeling and experimental identification[END_REF] for negligible geometric nonlinearities. As can be observed in Fig. 18, the relative error on the prediction of the resonance frequency remains small.

Stability and hysteresis effects

The bending of the resonance curve -that is observed either for direct or parametric excitation -gives rise to some hysteretic behavior. To illustrate that effect, the frequency response is plotted in Fig. 19 for a parametric excitation with acceleration amplitude U 0 = 6.73 g. The equilibrium position is stable outside of the domain delimited by the points A and D in Fig. 19. The stability of oscillatory solutions can be estimated by studying the jacobian of the system (3.8): In the steady-state solution (3.7), the amplitude a and the phase shift γ are indeed obtained as equilibrium solutions to the dynamical system (3.8). Such equilibria are (asymptotically) stable if all the eigenvalues of the jacobian have a negative real part [START_REF] Viterbo | Equations différentielles et systèmes dynamiques[END_REF]. As detailed in Appendix C, there is a simple geometrical interpretation of the stability condition for the problem at hand: The limit of stability is reached when the tangent vector to the resonance curve is vertical (points B and E in Fig. 19). The portion of the resonance curve that is bounded by the points B and E (showed as a dashed line) corresponds to unstable solutions: Any small disturbance will drive the system either to a stable oscillatory solution or to the equilibrium position. Points on the resonance curve that are outside of the portion delimited by B and E correspond to stable oscillatory solutions. The frequency shifts corresponding to the points A, B, D, E are denoted by σ(A), σ(B), σ(D) and σ(E), respectively. For frequency shifts between σ(D) and σ(E), observe that there are two stable oscillatory solutions. Depending on the initial conditions, the steady-state response will be either on the upper branch BF or the lower branch DE of the resonance curve (see [START_REF] Daqaq | Investigation of power harvesting via parametric excitations[END_REF] for a more detailed discussion along those lines). Similarly, for frequency shifts between σ(B) and σ(D), there are two possible steady-state solutions, namely the equilibrium position and an oscillatory solution (branch BF ). In particular, it can be noted that the resonance frequency ω M falls in that range. As a consequence, the steady-state solution achieving peak power (point M in Fig. 19) is reached for some but not all initial conditions. The stable steady-state oscillatory solution is unique only for frequency shifts between σ(E) and σ(A). For frequencies in that range, the electrical power is maximum at the point labeled F in Fig. 19: That point correspond to the maximum power that can be extracted without setting provision on the initial conditions.

σ [s -1 ] -4 -3 -2 -1 0 1 2

Dependence of the harvested energy on the amplitude of excitation

The peak power P M as given in Eq. (3.14) is plotted in Fig. 20 (solid lines) as a function of the excitation amplitude B, both for direct and parametric excitations. The peak power P M is obtained by setting the excitation frequency equal to the relevant resonant frequency for the type of excitation considered (direct or parametric). The exact expression of the resonant frequency is given by Eqs. (3.15-3.16 )). Besides, the load resistance is set to its optimal value as discussed in Sect. 5.1. Hence, the curves in Fig. 20 (solid lines) shows the maximum power that can be harvested for a excitation of given amplitude. In the case of direct excitations, the B -P M curve (solid line) is almost linear. In contrast, in the case of parametric excitation, the B -P M curve has more of a parabolic shape: the harvested power starts slowly but grows more rapidly with the amplitude of excitation. Similar remarks apply to the the half-power bandwidth, as represented in Fig. 21. Those results suggest that the benefits of parametric excitation over direct excitation can only appear at high excitation amplitude.

B [m]

In that regard, it should be kept in mind that the distortion of the resonance curve gets more pronounced as the excitation amplitude B increases. As a result, for B large enough, there are two stable solutions at the resonance frequency ω M , as discussed previously in Sect. 5.3. In such case, the optimal oscillatory solution is reached only for certain initial conditions. The dotted lines in Fig. 20 shows the maximum power that can be obtained independently on the initial state. Both for direct and parametric excitation, those curves are almost linear with comparable slopes. A similar observation can be made for the half-power bandwidth (Fig. 21).

Conclusion

In this paper, a Galerkin projection has been used to study the behaviour of a distributedparameter model of piezoelectric harvesters that included both piezoelectric and geometric nonlinearities. That model holds for moderately large geometric nonlinearities. Regarding the piezoelectric, mechanical nonlinearities (large strain effects) are assumed to dominate electric nonlinearities (large electric field effects). A cubic term appears in the stress-strain relation as a result, but the piezoelectric coupling is assumed to remain linear. Material nonlinearities in the substrate are neglected. The beam is furthermore assumed to be inextensible, which is reasonable for slender beams and excitation frequencies of the order of the first natural frequency ω 0 . As demonstrated in Sect. 4, the proposed model gives a fairly accurate description of the dynamics near the main resonance frequency and for moderately large amplitude of excitation, both in direct and parametric excitations. Most of the model parameters are readily available from a linear finite element analysis of the cantilever beam. The remaining nonlinear parameters can be obtained by an identification procedure that requires measuring the peak power and corresponding resonance frequency at a single excitation amplitude.

Various issues related to energy harvesting have been discussed in detail. Under certain conditions given in Sect. 5.1, the optimal resistance can be made independent on the excitation amplitude, which is a welcome feature for effective predictions in applications. Concerning direct excitations, it appears that the limits of validity of a linear model are relatively low (0.16 m.s -2 for the tested device) -lower than what could intuitively be expected. As a consequence, nonlinearities effects can be significant even for low-level ambient vibrations. For the tested device, ignoring nonlinearities leads to overestimate the harvested energy, as detailed in Sect. 5.2.

Compared to direct excitation, the benefits of parametric excitation are likely to occur only at high amplitudes, far above the stability limit of the equilibrium. Moreover, at such levels of amplitude, several steady-state solutions are stable. Therefore, depending on the initial conditions, the response of the system may not converge towards the most favorable steady-state solution.

where

M ij = L 0 µφ i φ j dx K i,j = L 0 EIφ i φ j dx α ijkl = α b ijkl + α p ijkl α b ijkl = L 0 EI φ i φ j φ k φ l + φ i φ j φ k φ l dx α p ijkl = 1 5 c 2 b h 2 + h p 5 - h 2 5 Lp 0 φ i φ j φ k φ l dx β ijkl = L 0 µ x 0 φ i φ j ds x 0 φ k φ l ds dx ξ ij = L 0 φ i φ j x L µds dx θ i = χφ i (L p ) N ijk = χ 2 Lp 0 φ i φ j φ k + φ i φ j φ k + φ i φ j φ k dx m i = L 0 µφ i dx
and (c i , d i ) are respectively the linear and the quadratic damping for the mode i. In addition, the following electrical equation has to be satisfied:

C p V + V R + n i=1 θ i ṙi + n i,j,k=1
L ijk r i r j ṙk = 0 (A.3)

with

L ijk = χ Lp 0 1 2 φ i φ j φ k + φ i φ j φ k dx
The system (A.2-A.3) is a set of n+1 coupled nonlinear equations that governs the evolution of (r 1 , • • • , r n , V ).

Appendix B. Stability of equilibrium for the piezoelectric beam under parametric excitation

Let T = 2π/ω be the period of the excitation. Setting X = (r, ṙ, V ), the set of equations (2.33) can be put in the generic form Ẋ(t) = F(X(t), t)

(B.1)
where F is T -periodic. As noted in Sect.2.3, the equilibrium state X 0 (t) = 0 is a particular solution of the nonlinear differential equation (2.33). Using Routh criterion [START_REF] Viterbo | Equations différentielles et systèmes dynamiques[END_REF] In the present case, the linearized equation at X 0 (t) = 0 can be written as M r + (K + ξ üB ) r + c ṙ -θV = 0,

C p V + V R + θ ṙ = 0. (B.4)
We are interested in finding the maximum value of B for which r = V = 0 is an asymptotically stable solution to (B.4). From Floquet theory, the limit of stability is reached when one of the Floquet multipliers has a null real part, i.e. when (B.4) admits a non zero solution of the form r(t) = e iµt p(t) , V (t) = e iµt q(t). (B.5)

where µ ∈ R and (p, q) are T -periodic functions. Using Fourier decomposition, the functions p and q can be written as First consider the case θ = c = 0, in which case (2.33) reduces to the Mathieu equation [START_REF] Landau | Mechanics[END_REF]. Using (B.3), it can be proved that µ = k 0 ω 2 for some k 0 ∈ N. It follows from (B.6) that the limit of stability is equal to 0 if K -M (kω + k 0 ω 2 ) 2 = 0 for some k, i.e. if ω = 2ω 0 n for some n ∈ N. Now consider that θ and c are small parameters. More precisely, we assume that c = O(ε) and θ 2 = O(ε) where ε 1. We want to find the limit of stability for ω = 2ω 0 + σ with σ = O(ε). To that purpose, the coefficients a k and the exponent µ in (B.6) are expanded as

p(t) =
a k = a 0 k + a 1 k , µ = ω 0 2 + µ 1
where a 1 k = a 0 k O(ε) and µ 1 = O(ε). Substituting in (B.6) and collecting the terms of the main order in ε, we obtain that a 0 k = 0 for k ∈ {-1, 0}. Note that it is necessary that (a 0 -1 , a 0 0 ) = (0, 0) for the functions p and q in (B.5) to be non zero. At the next order in ε, we obtain the relation 0 = -2M (2k + 1)(kσ + µ 1 )ω 0 a 0 k + (K -M (2k + 1) 2 ω 2 0 )a 1 k + 2Bξω 2 0 (a 0 k-1 + a 0 k+1 ) +ic(2k + 1)ω 0 a 0 k + Rθ 2 a 0 k i(2kω 0 + µ 1 ) 1 + iRC p (2kω 0 + µ 1 ) .

(B.7) In particular, using (B.7) with k = -1 and k = 0 yields the two equations 0 = -icω 0 a 0 -1 -ω 0 R(i-RCpω 0 )θ 2 1+C 2 p ω 2 0 R 2 a 0 -1 + 2M ω 0 (µ 1 -σ)a 0 -1 + 2Bω 2 0 ξa 0 0 , 0 = 2Bω 2 0 ξa 0 -1 + icω 0 a 0 0 -2M ω 0 µ 1 a 0 0 +

ω 0 R(i+Cpω 0 R)θ 2 1+C 2 p ω 2 0 R 2 a 0 0 . (B.8)
For (a 0 -1 , a 0 0 ) to be different from (0, 0), the determinant of the linear system (B.8) must vanish. That requirement implies that µ 1 = ω 0 x and 4B 2 ω 2 0 ξ 2 = (c +

Rθ 2 1 + R 2 C 2 p ω 2 0 ) 2 + ( R 2 C p θ 2 ω 0 1 + R 2 C 2 p ω 2 0 -M (ω -2ω 0 )) 2 (B.9)
The value of B defined by the above relation is the limit of stability at the excitation frequency ω.

In order to obtain a geometrical interpretation of the condition (C.4), consider the derivative of the frequency-response P (σ). Recalling that P (σ) = In ). In the unshaded domain, dP /dσ has the same sign as det J. Therefore, unstable solutions correspond to point such that dP /dσ < 0, i.e. such that the frequency response is decreasing. Similarly, in the shaded domain, unstable solutions correspond to point such that dP /dσ > 0, i.e. such that the frequency response is increasing. As a conclusion, limits of stability correspond to points where a change of sign occur in dP /dσ (setting aside the resonance point M ). For the case depicted in

Figure 1 :

 1 Figure 1: Piezoelectric cantilever beam under base excitation.

M

  r + (K + ξ üB ) r + β rr + ṙ2 r + α b r 3 + c ṙ + d ṙ| ṙ| + mv B =

2 y = -h 2 Figure 2 :

 222 Figure 2: Poling direction in the piezolectric patches.

Figure 4 :

 4 Figure 4: Frequency response for several values of the linear effective damping parameter ceff : (left) direct excitation with B = 0.01, (right) parametric excitation with B = 0.02 ( d = 0.1, κeff = 0, ξ = -1).

Figure 5 :

 5 Figure 5: Frequency response for several values of the quadratic damping parameter d : (left) direct excitation with B = 0.01, (right) parametric excitation with B = 0.02 (c eff = 0.01, κeff = 0, ξ = -1).

Figure 6 :

 6 Figure 6: Frequency response for several values of the effective nonlinear parameter κeff : (left) direct excitation with B = 0.01, (right) parametric excitation with B = 0.02 (c eff = 0.01, d = 0.1, ξ = -1).

Figure 7 :

 7 Figure 7: Piezoelectric cantilever used for experiments.

E p E e ρ p ρ e e 31 ε S 33 (

 33 MPa) (MPa) (kg.m -3 ) (kg.m -3 ) (C.m -2 ) (F.m -1 )

Figure 8 :

 8 Figure 8: Overall view of the experimental set-up.

Figure 9 :

 9 Figure 9: Set-up for direct excitation.

Figure 10 :

 10 Figure 10: Set-up for parametric excitation.

FrequencyFigure 11 :Figure 12 :

 1112 Figure 11: Spectral analysis of free vibrations.

s a mp l e n u mb e rFigure 13 :

 13 Figure 13: (left) Measurements of the voltage response V (t) in free vibration at small displacement amplitude; (right) Logarithmic ratio V n /V (0) as a function of n.

  s.m -1 ) (kg.m -1 ) (N.m -3 )

  Direct excitation, U 0 = 0.079 g.

Figure 15 :

 15 Figure 15: Numerical simulations of the frequency response in parametric excitation.

Figure 16 :

 16 Figure 16: Numerical simulation of the frequency response in direct excitation.

. 3 )Figure 17 :

 317 Figure 17: Optimal resistance under parametric resonance.

Figure 18 :

 18 Figure 18: Comparison of the linear and nonlinear models: relative error on the peak power and resonance frequency (left), frequency response for a direct excitation with amplitude U 0 = 0.4 m.s -2 (right).

3 aFigure 19 :

 319 Figure 19: Frequency response for parametric excitation with acceleration amplitude U 0 = 6.73 g: Stability of oscillatory solutions and hysteretic behavior.

Figure 20 :

 20 Figure 20: Electrical power at the resonance frequency, as a function of the amplitude of excitation: (red) direct excitation, (blue) parametric excitation.

Figure 21 :

 21 Figure 21: Half-power bandwidth as a function of the amplitude of excitation: (red) direct excitation, (blue) parametric excitation.

k∈Z a k e(B. 6 )

 6 ikωt , q(t) = k∈Z b k e ikωt Substituting in (B.4) leads to the following linear relations between the coefficients (a k , b k ):-M (kω + µ) 2 a k + Ka k + 1 2 ξBω 2 (a k-1 + a k+1 ) + ic(kω + µ)a k -θb k = 0 iC p (kω + µ)b k + 1 R b k + iθ(kω + µ)a k = 0.Eliminating b k leads to the relation-M (kω +µ) 2 a k +Ka k + 1 2 ξBω 2 (a k-1 +a k+1 )+ic(kω +µ)a k +Rθ 2 a k i(kω + µ) 1 + iRC p (kω + µ)= 0 The relation (B.6) holds for any k ∈ Z. At the limit of stability, the set of equations (B.6) admits a non zero solution {a k } k∈Z .

  Fig. C.22, points in the unshaded (resp. shaded) domain satisfy f ,γ (a * , γ * ) < 0 (resp. f ,γ (a * , γ * ) > 0). Note that f ,γ (a * , γ * ) vanishes at the resonance point of the frequencyresponse (denoted by M in Fig. C.22

  Fig. C.22, such a change of sign occurs at two points labeled as B and E. The same analysis can be transposed to direct excitations.

Figure C. 22 :

 22 Figure C.22: Geometrical interpretation of the condition of stability for steady-state oscillatory solutions.

Table 2 :

 2 Model parameters (obtained from FEA).

	)

  , X 0 is (asymptotically) stable if Y 0 (t) = 0 is an (asymptotically) stable solution to the linearized Note that (B.2) is a linear differential equation with T -periodic coefficients. The stability of solutions to such equations can be studied using Floquet theory: Y 0 (t) is asymptotically stable if all the Floquet exponents have a non positive real part. We recall that the distinctive property of the Floquet exponents µ i is that (B.2) admits a non zero solutionY i of the form Y i (t) = e µ i t P i (t)where P i is T -periodic and takes values in C. For later reference, we also recall that the Floquet multipliers µ 1 ,• • • ,µ N satisfy the relationµ 1 . • • • .µ N =

	equation
	Ẏ (t) =	∂F ∂X	(X 0 (t), t).Y (t).	(B.2)
		0	T	tr	∂F ∂X	(X 0 (t), t) dt.	(B.3)

  1 2 P 0 L 2 a * ,2 where P 0 is a constant, we have dP dσ = P 0 L 2 a * da * dσ .

	Now differentiating (C.2) with respect to σ gives
		J •	da * /dσ dγ * /dσ	+	0 2ω 0	=	0 0	.
	It follows that
	da * /dσ dγ * /dσ	= -J -1 •	0 2ω 0	= -	1 det J	g ,γ -g ,a f ,a -f ,γ	•	0 2ω 0	.
	Hence		dP dσ	= -2ω 0	P 0 L 2 a * det J	f ,γ (a

* , γ * ).

α p .

On a final note, we stress that further work is needed to extend the limits of validity of the proposed model, both in terms of frequencies and amplitude of the excitation. In particular, for frequencies that are not close to the main resonance and large amplitude of excitations, additional modes need to be taken into account in the Galerkin projection. Moreover, higher order nonlinearities than those considered in this paper -such as the nonlinear piezoelectric coupling -may come at play.

Appendix A. Galerkin approximation with n modes

If the displacement v is approximated as

where φ 1 , • • • , φ n are given mode shapes, then the coefficients r i (t) satisfy the n equations

β ijkl (r j r k rl + r j ṙk ṙl )

Appendix C. Stability of steady-state oscillatory solutions

The equations of modulation (3.6) can be put in the form

where

In (C.1), σ is the dimensionless frequency shift, as defined in (3.9). Steady-state oscillatory solutions correspond to equilibria of the dynamical system (C.1), i.e. to values (a * , γ

To study the stability of such equilibria, consider the jacobian J of the system, i.e. the matrix

where the partial derivatives of f and g are evaluated at (a * , γ * , σ). Let λ 1 ∈ C and λ 2 ∈ C be the eigenvalues of J, i.e the solution of the quadratic equation

An equilibrium (a * , γ * ) is (asymptotically) stable if the real parts of λ 1 and λ 2 are negative [START_REF] Viterbo | Equations différentielles et systèmes dynamiques[END_REF]. In the present case, an important observation is that tr J < 0. Therefore, setting ∆ = (tr J) 2 -4 det J, three cases can occur:

-if ∆ < 0, then λ 2 = λ1 and consequently λ 1 + λ 2 = 2 Re(λ 1 ) = tr J < 0. Hence (a * , γ * ) is a stable equilibrium.

-if ∆ ≥ 0 and det J ≥ 0 then (λ 1 , λ 2 ) are real numbers and have the same sign. Since λ 1 + λ 2 = tr J < 0, λ 1 and λ 2 are both negative. It follows that (a * , γ * ) is a stable equilibrium.

-if ∆ ≥ 0 and det J < 0 then (λ 1 , λ 2 ) are real numbers and have opposite sign. Hence (a * , γ * ) is an unstable equilibrium.

Since ∆ = (tr J) 2 -4 det J, note that det J < 0 automatically implies that ∆ ≥ 0. For the problem at hand, the condition of instability thus reduces to det J < 0.

(C.4)