# Optimisation and performance comparison of direct and parametric piezoelectric energy harvester with geometrical nonlinearities

Hugo Fayolle, Christophe Giraud-Audine, Olivier Thomas

► To cite this version:

Hugo Fayolle, Christophe Giraud-Audine, Olivier Thomas. Optimisation and performance comparison of direct and parametric piezoelectric energy harvester with geometrical nonlinearities. Journal of Sound and Vibration, 2025, 619, pp.119343. 10.1016/j.jsv.2025.119343. hal-05394631

HAL Id: hal-05394631

https://hal.science/hal-05394631v1

Submitted on 3 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ELSEVIER

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Optimisation and performance comparison of direct and parametric piezoelectric energy harvester with geometrical nonlinearities

Hugo Fayolle<sup>a ID</sup>, Christophe Giraud-Audine<sup>b ID</sup>, Olivier Thomas<sup>a ID,*</sup>

## ARTICLE INFO

### Keywords:

- Piezoelectric energy harvester
- Geometrical nonlinearities
- Arbitrary laminated model
- Large rotation
- Parametric excitation
- Perturbation methods

## ABSTRACT

The behaviour and performance of a flexible beam piezoelectric energy harvester (PEH) in direct and parametric excitation, including the effect of geometrical nonlinearities, are addressed in this paper. First, the electromechanical modelling of the harvester is addressed. A Timoshenko geometrically exact model of a laminated piezoelectric beam in large rotation, including parametric excitation, is proposed, extending previous results of the literature. Then, considering cantilever boundary conditions, it is simplified under Euler-Bernoulli, inextensible assumptions and third order Taylor expansion, suitable for modal expansion. Two perturbation methods are tested to compute vibratory response under direct and parametric excitations. They are both found inaccurate for large amplitude oscillations, leading to a preference for numerical solving by continuation of periodic solutions. PEH behaviours and performances are finally carefully estimated regarding the optimal harvested power in a shunted resistor at resonance. Comparisons between direct and parametric forcing are proposed and the effect of geometrical nonlinearities is estimated. An interesting result is that the harvested power under parametric excitation cannot compete with the one under direct excitation.

## 1. Introduction

Ambient energy harvesting is a large studied field in order to replace or assist batteries for small electronic devices supply and data sharing for Internet of Things [1], in vivo applications [2] or wireless sensor networks [3] applications. The purpose is to avoid battery change operation for the system to be autonomous all along its lifetime. An important research area concerns vibration energy harvesters and more precisely piezoelectric energy harvesters (PEH), which have been the subject of numerous literature during the past decades (see e.g. the recent reviews on the subject [4–6]). The basic architecture is composed of three components: a mechanical resonator, a piezoelectric transducer and an electronic circuit. A popular design of the resonator is a thin flexible cantilever beam onto which one or several piezoelectric patches are glued, to take advantage of bending deformations [4]. The behaviour and the performance of these systems are well known in the case of linear resonators [7], that are known to be limited in terms of frequency bandwidth, such that taking advantage of nonlinearities is an active subject of research to enhance the performances [8,9]. Current strategies, among others, include the use of resonators with high piezoelectric coupling and quality factor with an adapted electronics [10], Duffing-like response [11,12], bistable resonators [13,14], resonators with internal resonances [15] and chaos [12] and harvesters for rotating machinery, taking advantage of an excitation due to gravity [16].

The present study focuses on the investigation of the energy harvesting performance and behaviour of a resonator under parametric excitation. In this situation, the cantilevered resonator is orientated such that the main input vibration direction is parallel

<footnote><sup>a</sup>Arts et Métiers Institute of Technology, LISPEN, F-59000 Lille, France</footnote>

<footnote><sup>b</sup>Arts et Métiers Institute of Technology, L2EP, Univ. Lille, Centrale Lille, JUNIA, F-59000 Lille, France</footnote>

<footnote>* Corresponding author.
E-mail addresses: hugo.fayolle_1@ensam.eu (H. Fayolle), christophe.giraud-audine@ensam.eu (C. Giraud-Audine), olivier.thomas@ensam.eu (O. Thomas).</footnote>
Nomenclature

<table><tr><td>ō, ō</td><td>Dimensionless version of ◦ variable</td></tr><tr><td>e&#773;<sup>k</sup><sub>31</sub></td><td>Piezoelectric coefficient of the kth layer</td></tr><tr><td>ε, S</td><td>Strain and stress tensor fields</td></tr><tr><td>E, D</td><td>Electric field and displacement</td></tr><tr><td>e<sub>x</sub>, e<sub>y</sub>, e<sub>z</sub></td><td>Basis vectors</td></tr><tr><td>χ, χ&#773;</td><td>Coefficients of the linear and nonlinear piezoelectric coupling terms</td></tr><tr><td>ū<sub>0</sub>, īū<sub>0</sub></td><td>Axial and transverse prescribed accelerations</td></tr><tr><td>Δ</td><td>Dirac delta functions difference</td></tr><tr><td>δ</td><td>Parametric modal forcing</td></tr><tr><td>δ<sub>cr</sub></td><td>Critical parametric forcing</td></tr><tr><td>ε<sup>k</sup><sub>33</sub></td><td>Dielectric constant of the kth layer</td></tr><tr><td>γ</td><td>Amplitude of direct/parametric external acceleration</td></tr><tr><td>Γ, Π</td><td>Coefficients of the cubic nonlinear terms in short circuit</td></tr><tr><td>Î</td><td>Coefficient of the cubic nonlinear term in open circuit</td></tr><tr><td>ω&#775;<sub>l</sub></td><td>Natural frequency of the ith mode in open circuit</td></tr><tr><td>Θ&#775;<sub>p</sub></td><td>Effective bending piezoelectric constant for a non-symmetric laminae</td></tr><tr><td>Θ&#775;<sub>p</sub></td><td>Effective capacitance for a non-symmetric laminae</td></tr><tr><td>D&#775;</td><td>Effective bending stiffness for a non-symmetric laminae</td></tr><tr><td>κ</td><td>Electromechanical coupling factor</td></tr><tr><td>P<sub>η</sub></td><td>Power</td></tr><tr><td>T<sub>η</sub></td><td>Kinetic energy</td></tr><tr><td>γ<sub>η</sub></td><td>Potential energy</td></tr><tr><td>ν<sub>k</sub></td><td>Poisson's ratio of the kth layer</td></tr><tr><td>Ω</td><td>Driving frequency</td></tr><tr><td>ω<sub>l</sub></td><td>Natural frequency of the ith mode in short circuit</td></tr><tr><td>Φ<sub>l</sub></td><td>Deformed shape of the ith mode in short circuit</td></tr><tr><td>ρ<sub>k</sub></td><td>Mass density of the kth layer</td></tr><tr><td>τ<sub>e</sub></td><td>Electrical time constant</td></tr><tr><td>θ</td><td>Rotation of cross section</td></tr><tr><td>ξ</td><td>Modal damping ratio</td></tr><tr><td>Ξ<sub>p</sub>, Θ<sub>p</sub></td><td>Piezoelectric constants of the laminated beam and the pth patch</td></tr><tr><td>A, B, D</td><td>Stiffness constant of the laminated beam</td></tr><tr><td>a, α</td><td>amplitude/phase pair of q(t)</td></tr><tr><td>a<sub>h</sub>, ā<sub>h</sub></td><td>Cosine and sine Fourier coefficients of q(t)</td></tr><tr><td>b, β</td><td>amplitude/phase pair of Q(t)</td></tr><tr><td>C</td><td>Equivalent capacitance of the piezoelectric patches</td></tr><tr><td>e, Ę, Ē</td><td>Axial, curvature and shear strains</td></tr><tr><td>E<sub>z</sub>, D<sub>z</sub></td><td>Transverse components of E and D</td></tr><tr><td>F</td><td>Direct modal forcing</td></tr><tr><td>H, Ę</td><td>Shear stiffness and shera correction coefficient of the laminated beam</td></tr><tr><td>K<sub>l</sub></td><td>Modal stiffness of the ith mode</td></tr><tr><td>L<sub>p</sub>, h<sub>p</sub>, b<sub>p</sub></td><td>Length, thickness and width of the pth piezoelectric layer</td></tr><tr><td>m, I</td><td>Inertia constants of the laminated beam</td></tr><tr><td>M<sub>l</sub></td><td>Modal mass of the ith mode</td></tr><tr><td>N, M, T</td><td>Beam generalised stresses (axial force, bending moment and shear force)</td></tr><tr><td>n, p, q</td><td>External loading: axial force, transverse force and bending moment per unit length</td></tr><tr><td>P<sub>η</sub></td><td>Averaged power over one period</td></tr><tr><td>q</td><td>Modal coordinate</td></tr><tr><td>R</td><td>Resistance of the electrical charge</td></tr><tr><td>S, Λ</td><td>Correction coefficients for the direct and parametric drive</td></tr><tr><td>S<sub>xx</sub>, S<sub>xz</sub></td><td>Axial and shear components of S</td></tr><tr><td>t</td><td>Time</td></tr><tr><td>u, w</td><td>Axial and transverse displacement fields</td></tr></table>
<table><tr><td>V, Q</td><td>Voltage and charge at the terminal of the electric charge</td></tr><tr><td>V<sub>p</sub>, Q<sub>p</sub></td><td>Voltage and charge pair of the pth piezoelectric layer</td></tr><tr><td>x, y, z</td><td>Coordinates of any point M in the beam continuum</td></tr><tr><td>Y<sub>k</sub>, G<sub>k</sub></td><td>Young's and shear moduli of the kth layer</td></tr><tr><td>z<sub>n</sub></td><td>Transverse location of the neutral line of the laminated beam</td></tr><tr><td>SC, OC</td><td>Short and open circuit</td></tr></table>

<figure><img src="image_5.png" /><figcaption>Fig. 1. Parametric energy harvester with large rotation oscillations under (a) direct or (b) parametric excitation.</figcaption></figure>

to the beam axis, such that the bending response appears thanks to a dynamic instability, above a threshold. In this case, a linear model for the resonator predicts exponentially growing vibrations, which are bounded in practice by geometrical nonlinearities, which have to be mandatory included in the model to predict the harvested power. This article thus addresses the impact of the geometrical nonlinearities on the harvested power, for a PEH under parametric excitation, but also in the case of direct excitation, in order to systematically compare the two forcing situations, direct and parametric, sketched in Fig. 1.

Whereas direct excitation of the resonator is addressed in the majority of contributions about PEH, parametric energy harvesting was addressed only scarcely. The pioneering work was proposed in [17], where the authors study a cantilever beam axially moved with a one-mode model with a quadratic damping. In [18] the same problem is studied with a more refined model that includes geometrical nonlinearities and the piezoelectric coupling. A particular practical configuration is proposed in [19], for which the cantilever beam is included in a rotating wind turbine and parametrically excited thanks to magnets. Experimental proofs of concept are proposed in [20,21]. A theoretical short study of the harvested power of a Mathieu-Duffing type oscillator is proposed in [22]. Then, in a series of recent articles [23–26], the common effect of simultaneous direct and parametric excitation is investigated. Finally, the combination of a parametric PEH to other sources of nonlinearities is explored in [27] with motion limiters, in [28] with time-varying snap-through effects and in [29] with a 1:3 internal resonance. It is worth citing [30] that explores the same concept of parametric energy harvesting with an electromagnetic transduction, as well as [31,32] in which the parametric amplification concept is explored, to enhance the performances.

Appropriate PEH modelling requires a well suited model of laminated cantilever with piezoelectric patches. Cantilever beams are very flexible structures in bending, such that large amplitude vibrations are observed for a low excitation level, leading to geometrical nonlinearities. These nonlinearities are due to large rotation and are different from those due to axial/bending coupling, that is observed when the axial motion is restrained by the boundary conditions (see e.g. [33]). For slender resonators, beam models are well adapted and so-called geometrically exact models naturally include the nonlinear effect of the large rotations of the cross sections, without restriction (see the pioneering works [34] for a beam kinematics without shear deformation and [35] that includes it, [36,37] for overviews on this model and [38–41] for very large deformation simulations). An interesting simplification of this model, thanks to an inextensible assumption and a Taylor expansion, was proposed in [42] and also widely used ever since. However, most of those studies consider a homogeneous material, in contradiction with piezoelectric resonators which are often laminated, including several layers of different materials. For macro structures, the most common architecture is an elastic thin beam/plate/shell onto which one or several piezoelectric thin patches are glued [43,44]. For micro/nano electromechanical systems (M/NEMS), more layers of different materials are often encountered [45,46]. The most simple theories replace the laminated cross-section of the slender structure by an equivalent single homogeneous layer [47]. Whereas most studies do not include geometrical nonlinearities [48], some of them take into account membrane/bending nonlinear coupling using von Kármán assumptions [47,49,50], but large rotation laminated beam theories including piezoelectric layers are scarce. In [51], the authors propose an extension of the classical inextensible model of [42] to the laminated piezoelectric case, with an attempt to extend it to a geometrically exact theory in [52]. In [53], a geometrically exact piezoelectric beam model is solved by numerical integration to estimate the harvested power under strong nonlinear dynamics and a similar model is proposed in [54] for energy harvesting under fluid structure interactions.

In this context, the main goal of this article is to optimise the energy harvesting of a parametrically excited cantilever beam with piezoelectric patches and to compare its performance to the one of the same resonator under direct excitation, considered as a reference. A geometrically exact model of a laminated beam with piezoelectric layers, under both direct and parametric excitation, coupled to an electric circuit is proposed. This model is then rigorously simplified with an inextensible assumption, a third order
<figure><img src="image_6.png" /><figcaption>Fig. 2. Piezoelectric laminated beam under study: (a) a typical configuration with detail of (b) the Timoshenko kinematics and (c) of a piezoelectric layer.</figcaption></figure>

Taylor expansion and a one eigenmode projection, to obtain a reduced order analytical model of the system. The derivation, reported in Section 2, includes all the refinements, some of them being overlooked in previous works of the literature with similar models [17,20,52,54,55]. Then, in the attempt of obtaining closed-form expressions for the vibration amplitude and energy, two first order classical perturbative solutions, stemming from a one-harmonic balance method and a multiple scale method, are compared in Section 3. Because no quadratic drag-like damping is included in the model (on the contrary to previous works [17,20,55]), it is found that those approximated solutions are far from reference numerical solutions (obtained with numerical continuation), which are thus preferred in the remaining of this work. The optimisation of the parametric PEH is addressed in Section 4 and compared to direct PEH, considered as a reference. This optimisation is done by (1) carefully selecting a reduced set of independent parameters using proper scaling of the model equations and (2) performing an exhaustive search on this reduced set of parameters to understand their optimal values that maximise the harvested power. The effect of geometrical nonlinearities on a cantilever resonator under direct excitation is also addressed. Finally, an analogy with a pendulum is done in Section 5, to give more insights in the parametric PEH performances.

## 2. Governing equations

The target of this section is to obtain a model of a laminated piezoelectric beam allowing large rotation vibrations, including both direct or parametric forcing. The proposed model is based on the geometrically exact beam model of [35], as written in [37], which is extended to include a laminated cross section with piezoelectric layers, based on an equivalent single layer theory, as proposed in [47,48]. A Timoshenko kinematics, allowing for shear deformations, is selected for the sake of generality. Then, this general model is simplified to allow an analytical treatment, namely by (i) neglecting shear with an Euler-Bernoulli kinematics, (ii) considering cantilever boundary conditions and an inextensible middle line and (iii) performing a third order Taylor expansion in the transverse displacement field. This allows using modal expansions to obtain a discretised reduced order model in order to study the nonlinear dynamics and performance of the resonator for PEH, in the next sections.

### 2.1. A piezoelectric laminated geometrically exact beam model

#### 2.1.1. Beam deformations and strains

Following the derivations in [37], we consider a straight cantilever beam of length L composed of several layers of different materials (Fig. 2(a)). Following a Timoshenko kinematics restricted to the plane ($e_x, e_z$) (with $e_x$ the unit vector in the axial direction and ($e_x, e_y, e_z$) an orthonormal basis), we assume that any cross section initially plane undergoes a rigid body motion (Fig. 2(b)). If x is the axial location of the centre G of the cross section (due to the laminated nature of the cross section, this point may not be the geometrical centre, as discussed in Section 2.2) and if $u(x,t)e_x + w(x,t)e_z$ denotes its displacement field, the one of any point M in the cross section such that $GM = ye_y + ze_z$ is:

$$u = [u - z \sin \theta]e_x + [w + z(\cos \theta - 1)]e_z, \quad (1)$$

where $\theta(x,t)$ is the rotation of the cross section around $e_y$.

The strains are obtained using the Green-Lagrange tensor:

$$\epsilon = \frac{1}{2}(\nabla u + \nabla^T u + \nabla u^T \nabla u), \quad (2)$$
where $\nabla \circ$ is the tensor gradient of vector $\circ$. Since the beam is assumed to allow large displacements $w(x,t)$ and rotations $\theta(x,t)$ but small local strains, it is possible to simplify the Green–Lagrange strains $\epsilon$ by linearising the local strains while keeping the rigid body rotations of the cross-sections exact, using a so-called consistent linearisation procedure, described in [37]. This simplified strains, equivalent to those of the standard Cosserat rod theory [35,36,56], reads:

$$ \bar{\epsilon} = \begin{pmatrix} e - z\bar{\kappa} & \bar{\gamma}/2 & 0 \\ \bar{\gamma}/2 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \qquad (3) $$

in basis $(e_x, e_y, e_z)$, with

$$ \left\{ \begin{aligned} e &= (1 + u') \cos \theta + w' \sin \theta - 1 \\ \bar{\kappa} &= \theta' \\ \bar{\gamma} &= -(1 + u') \sin \theta + w' \cos \theta, \end{aligned} \right. \qquad (4) $$

being respectively the axial strain, curvature and shear strain and $\circ' = \partial\circ/\partial x$. One must remark that the shear strain $\bar{\gamma}$ is close, but not exactly equal, to the shear angle $\bar{\gamma} = \theta - \varphi$ defined in Fig. 2(b), with $\varphi$ the rotation angle of the normal to the middle line of the beam (refer to [37] for details).

### 2.1.2. Laminated structure and constitutive law

We consider that the beam is composed of several parts in the $x$ (axial) direction, as sketched in Fig. 2(a). In each part, the material and geometrical properties of the cross sections are uniform in the $x$-direction. In the $z$ (transverse) direction, which is the direction of bending, each cross section is constituted of $K$ plane layers ($K \in \mathbb{N}^*$), each one being made of a transverse isotropic material, eventually piezoelectric. The $k$th layer ($k = 1, \dots, K$) is situated between heights $z_{k-1}$ and $z_k$, as shown in Fig. 2(b); its thickness (along $e_z$) is $h_k = z_k - z_{k-1}$ and width (along $e_y$) is $b_k$. The coordinate system origin $O$ is situated at one end of the beam in the $x$-direction and the structure of the beam is symmetrical with respect to the plane $(O, e_x, e_z)$. In the $z$-direction, the origin location ($z=0$) is arbitrary and will be discussed in Section 2.2.

The strains $\bar{\epsilon}$ are work conjugated<fnref n="1" /> to a stress tensor $S$ (which is, because of the consistent linearisation, a Biot stress tensor instead of a second Piola–Kirchhoff one [58]). In addition, the electric state of the piezoelectric layers is defined by the electric $E$ and electric displacement $D$ fields, also work conjugated [60,61]. According to classical beam theories, the axial stress is dominant with respect to the two transverse stresses $S_{xx} \gg S_{yy}, S_{zz} \Rightarrow S_{yy} \approx S_{zz} \approx 0$. Moreover, since the piezoelectric patches are polarised along the $z$-direction, the $z$ component of the electric field is dominant and $E_z \gg E_x, E_y \Rightarrow E_x \approx E_y \approx 0$.

We focus in this article on a linear piezoelectric constitutive law, neglecting any nonlinear piezoelectric effects that may appear for large levels of electric field [18,62,63]. Including those effects in our model is left for upcoming studies, even if they might have a first order effect by changing the hardening/softening behaviour of some resonances [64–66]. Then, since the local strains are small, a linear relationship between stresses and strains is sufficient (a so-called Saint-Venant-Kirchhoff law). With those assumptions, the three-dimensional piezoelectric constitutive law in the $k$th layer reduces to, with Voigt notations [50,60,61]:

$$ \left\{ \begin{aligned} S_{xx} &= Y_k(e - z\bar{\kappa}) - \bar{\epsilon}_{31}^k E_z \\ S_{xz} &= G_k \bar{\gamma} \\ D_z &= \bar{\epsilon}_{31}^k (e - z\bar{\kappa}) + \epsilon_{33}^k E_z \end{aligned} \right. \qquad (5) $$

with $D_z$ the electric displacement along the $(z,3)$ direction, work conjugated to $E_z$. The material properties of the $k$th layer are its Young’s modulus $Y_k$ in the $e_x$ direction and its shear modulus $G_k$ in the $(x, z)$ plane (for an isotropic layer, $G_k = Y_k/[2(1 + v_k)]$) with $v_k$ the Poisson’s ratio; however, standard piezoelectric ceramics are transverse isotropic, such that $G_k$ is not directly related to $Y_k$ in this case), its piezoelectric constant $\bar{\epsilon}_{31}^k$ and its dielectric constant $\epsilon_{33}^k$. For any purely elastic layer, the piezoelectric constant vanishes, $\bar{\epsilon}_{31}^k = 0$.

### 2.1.3. Generalised forces

Among all the elastic layers in a given part of the beam, some of them are piezoelectric, their indices being gathered in a subset $P$ of $\{1, \dots, K\}$. For any of them, as sketched in Fig. 2(c), the upper and lower surfaces are considered electrically conducting. This is realised in practice by the upper and lower layers, called electrodes. For macro applications, these electrodes are often very thin and made of a conducting material, with their mechanical stiffness and inertia often neglected [65,67,68]. For M/NEMS, dedicated upper and lower conducting layers are used (see examples in [45,46]). These conducting layers are naturally non-piezoelectric, since piezoelectric materials are by essence dielectric [69]. Since the piezoelectric layers are considered thin (their thickness is small with respect to their length, $h_k \ll L_p$), we assume that the electric field is uniform in the domain it occupies [43], such that the electric potential in the $p$th piezoelectric layer is:

$$ V_p(t) = -E_p(t)h_p, \qquad (6) $$

with $V_p(t)$ the potential difference between its upper and lower electrodes.

<footnote n="1">¹ A pair of variables that are work conjugated means that their product gives a thermodynamic potential energy, see e.g. [57, §3.5.2] or [58,59]</footnote>
The generalised forces are defined by integrating the stresses of Eqs. (5) across the surface S of the cross section, such that:

$$N = \int_S S_{xx} d\sigma = Ae - B\bar{\kappa} + \sum_{p \in P} \Xi_p V_p, \quad (7a)$$

$$M = \int_S zS_{xx} d\sigma = Be - D\bar{\kappa} + \sum_{p \in P} \Theta_p V_p, \quad (7b)$$

$$T = \int_S S_{xz} d\sigma = \bar{k}H\bar{\gamma}, \quad (7c)$$

with N the axial force, M the bending moment and T the shearing force, $\bar{k}$ the shear correction factor [47] and $d\sigma$ is an infinitesimal surface element. The lamination dependent factors are [48]:

$$A = \sum_{k=1}^{K} b_k h_k Y_k, \quad B = \sum_{k=1}^{K} b_k h_k \bar{\bar{\kappa}}_k Y_k, \quad D = \sum_{k=1}^{K} b_k \left( \frac{h_k^3}{12} + h_k \bar{\bar{\kappa}}_k^2 \right) Y_k, \quad (8a)$$

$$\Xi_p = b_p \bar{\bar{e}}_{31}^p, \quad \Theta_p = b_p \bar{\bar{\zeta}}_p \bar{\bar{e}}_{31}^p, \quad H = \sum_{k=1}^{K} b_k h_k G_k, \quad (8b)$$

with $\bar{\bar{\zeta}}_k = (z_k + z_{k+1})/2$ the mean height of the kth layer. (A, D, B, H) are respectively the axial, bending, axial-bending coupling and shear stiffnesses of the laminated beam and $(\Xi_p, \Theta_p)$ are the axial and bending piezoelectric coupling coefficients. All those coefficients are piecewise functions of x because of the different parts of the beam in the axial direction.

In addition, the potential $V_p(t)$ of the pth piezoelectric layer is work conjugated to the electric charge $Q_p(t)$, that is the sum of the free electric charge contained in the upper electrode (by electrical influence, the charge in the lower electrode is $-Q_p(t)$). These free electric charges are contained in the surface $S_{el}$ separating the pth piezoelectric layer and its upper electrode, and are linked to the integral of the electric displacement jump across the surface:

$$Q_p = \int_{S_{el}} [D_z \frac{z_k^+}{\bar{\bar{\kappa}}} d\sigma = b_p \int_{x_-}^{x_+} -D_z(x, z_k) dx = -\Xi_p \int_{x_-}^{x_+} e(x) dx + \Theta_p \int_{x_-}^{x_+} \bar{\bar{\kappa}}(x) dx + C_p V_p, \quad (9)$$

since $D_z = 0$ inside the (conducting) electrode. One thus obtains:

$$Q_p = -\Xi_p L_p \bar{\bar{e}}_{x_-}^{x_+} + \Theta_p [\theta_{x_-}^{x_+} + C_p V_p], \quad (10)$$

where $C_p$ is the electric capacitance of the undeformed (blocked) piezoelectric patch, defined by:

$$C_p = c_{33}^p \frac{b_p L_p}{h_p}, \quad (11)$$

with $x_-$, $x_+$ the piezoelectric patch ends axial positions and $L_p = x_+ - x_-$ its length. One observes that the electric charge created in the electrodes of the pth piezoelectric patch by the mechanical deformation is composed of two terms, proportional to: (i) the averaged axial strain $e_{x_-}^{x_+} = 1/L_p \int_{x_-}^{x_+} e(x) dx$ and (ii) the difference of rotation $[\theta]_{x_-}^{x_+} = \theta(x_+, t) - \theta(x_-, t)$ of its two end cross sections.

### 2.1.4. Equations of motion

The equations of motion are derived from the virtual work principle and are identical to those of a purely elastic and homogeneous beam. They read [37,40]:

$$\left\{ \begin{array}{l} m\ddot{u} = (N \cos \theta - T \sin \theta)' + n \\ m\ddot{w} = (N \sin \theta + T \cos \theta)' + p \\ I\ddot{\theta} = T(1+e) - N\bar{\gamma} - M' + q, \end{array} \right. \quad (12a)$$

$$m\ddot{w} = (N \sin \theta + T \cos \theta)' + p \quad (12b)$$

$$I\ddot{\theta} = T(1+e) - N\bar{\gamma} - M' + q, \quad (12c)$$

with

$$m = \sum_{k=1}^{K} b_k h_k \rho_k, \quad I = \sum_{k=1}^{K} b_k \left( \frac{h_k^3}{12} + h_k \bar{\bar{z}}_k^2 \right) \rho_k, \quad (13)$$

the mass per unit length and the moment of inertia of the cross section, with $\rho_k$ the density of the kth layer. One must notice that the inertial part of Eq. (12c) is under this form only if the origin for the lamination heights $z_k$ (the location where $z = 0$) is chosen such that $\int_S z\rho d\sigma = 0$. Since the rotatory inertia will be neglected in Section 2.2 and after, this point will not be an issue. In addition, $n(x,t)$, $p(x,t)$ and $q(x,t)$ are the external axial force, transverse force and moment per unit length, respectively, and $\dot{o} = \partial\circ/\partial t$.

Eqs. (12) along with the electrical equation (10) fully governs the electromechanical behaviour of the piezoelectric laminated beam, in deformations with arbitrary cross section rotation $\theta(x, t)$. In the general case, finding solutions for the three kinematics fields $(w, u, \theta)$ and the electrical quantities $(V_p, Q_p)$ relies on a numerical solving, for instance with the finite-element discretisation method (see [40,70] for the purely elastic case). In the sequel, those equations are simplified to allow a semi-analytical treatment.
<figure><img src="image_7.png" /><figcaption>Fig. 3. Several particular cases: (a,c): symmetric lamination; (b) asymmetric lamination with one part; (d) asymmetric lamination with several parts. Piezoelectric layers are shown in red and the middle lines for each configurations are shown with a dashed-dotted line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 2.2. Third order inextensible model

This section aims at extending the inextensible model of [37,41,42], dedicated to beams made of a homogeneous material, to the present case of a laminated structure including piezoelectric layers. The assumptions are (i) an Euler–Bernoulli kinematics (shear strain neglected, $\bar{\gamma} = 0$) in addition to a negligible rotatory inertia ($I = 0$), (ii) and inextensible assumption ($e = 0$), valid with free boundary conditions at $x = L$, that allows eliminating the axial motion $u$ in favour of the transverse motion $w$, and (iii) a third order truncation in the transverse displacement field $w$.

Since energy harvesting is under concern, the considered external excitation is an acceleration field $\Gamma_0(t) = \ddot{u}_0(t)e_x + \ddot{w}_0(t)e_z$. To include this effect in the equations of motion, one has to replace in Eqs. (12) ($u, w$) by $(\ddot{u} + u_0, \ddot{w} + w_0)$, with $(\ddot{u}, \ddot{w})$ being the displacement fields with respect to the moving accelerated frame [40]. To simplify the derivations, $(\ddot{u}, \ddot{w})$ are simply denoted by $(u, w)$ in the following, showing that the external axial and transverse forces per unit length are uniform and reads:

$$n(x, t) = -m\ddot{u}_0(t), \quad p(x, t) = -m\ddot{w}_0(t). \qquad (14)$$

Since no external moment per unit length is applied, $q = 0$.

The core assumption of the present inextensible model is assumption (ii), that enforces the curvilinear length of the beam to remain constant during deformation. This may be valid only for the neutral line, i.e. the line onto which the axial strain $\bar{\epsilon}_{xx}$ is zero during bending deformation. In the present case of a laminated structure, considering first no piezoelectric loading ($V_p = 0 \forall p$), the location $z_n$ of the neutral line verifies zero axial strain ($\epsilon_{xx} = 0$) under bending without tension ($N = 0$):

$$\left\{ \begin{array}{l} \bar{\epsilon}_{xx} = e - z_n \bar{\kappa} = 0 \\ N = Ae - B \bar{\kappa} = 0 \end{array} \right. \Rightarrow z_n = \frac{B}{A}. \qquad (15)$$

As a consequence, one may verify that by choosing $z_n$ as the origin of the z-coordinates, i.e. the location of G the cross-section centre, one obtains $B = 0$ with Eq. (8b): the bending and axial motion are uncoupled in Eqs. (7) (i.e. ($N, M$) depends only on ($e, \bar{\kappa}$), respectively).

Then, several situations may be considered, as sketched in Fig. 3. In the case of a symmetric lamination (two examples are shown in Figs. 3(a,c)), it is possible to find a single neutral line, common to all parts of the beam, in the middle of the cross sections in the z-direction. It is also the case for an asymmetric lamination uniform in the x direction (Fig. 3(b)), with the z location of the neutral line being defined by Eq. (15). On the contrary, in the case of an asymmetric lamination with several different parts in the x-direction (Fig. 3(d)), the neutral line would be discontinuous, such that there is no continuous line for which $e = 0$ in bending for all $x \in [0, L]$. A way of managing these particularities in the model is first to eliminate $e$ between Eqs. (7a,b) and (9), as done in the linear case in [48]:

$$M = \frac{B}{A}N - \hat{D}\bar{\kappa} + \sum_{p \in P} \hat{\Theta}_p V_p, \qquad (16)$$

$$Q_p = -\frac{\bar{\epsilon}_p L_p}{A} \bar{N}_{x-}^{x+} + \hat{\Theta}_p [\theta]_{x-}^{x+} + \hat{C}_p V_p + \sum_{j \neq p} C_{pj} V_j, \qquad (17)$$

with the effective bending stiffness, piezoelectric constants and capacitances:

$$\hat{D} = D - \frac{B^2}{A}, \quad \hat{\Theta}_p = \Theta_p - \frac{B \bar{\epsilon}_p}{A}, \quad \hat{C}_p = C_p + \frac{L_p \bar{\epsilon}_p^2}{A}, \quad C_{pj} = \frac{L_p \bar{\epsilon}_p \bar{\epsilon}_j}{A}, \qquad (18)$$
and $\bar{N}_{-}^{x+} = 1/L_p \int_{-L_p}^{L_p} e(x) dx$ the averaged axial force. In the two above Eqs. (16), (17), instead of writing the variables $(N, M, Q_p)$ as a function of $(e, \bar{\kappa}, V_p)$, the role of $e$ and $N$ is exchanged such that $(e, M, Q_p)$ are now written as a function of $(N, \bar{\kappa}, V_p)$.

Then, the third order inextensible model is written following the procedure of [37,41,42]. Assumption (i) leads to neglect the shear strain ($\bar{\gamma} = 0$) as well as the rotatory inertia ($I = 0$), such that Eq. (12c) reduces to $T = M'(1+e)$. Then, $N$ is obtained with (12a) and assumption (iv) $(N(L,t) = T(L,t) = 0, \forall t \ge 0)$ as:

$$N = \frac{1}{\cos \theta} \int_{L}^{x} m(\xi)(\ddot{u} + \ddot{u}_{0}) d\xi + \frac{M' \tan \theta}{1+e}, \qquad (19)$$

such that it can be eliminated in Eqs. (12b), yielding:

$$m(x)\ddot{w} - \left(\frac{M'}{(1+e)\cos\theta}\right)' - \left[\tan\theta \int_{L}^{x} m(\xi)(\ddot{u} + \ddot{u}_{0}) d\xi\right]' = -m(x)\ddot{w}_{0}. \qquad (20)$$

Finally, assumption (i) ($\theta = \varphi$) leads to the following kinematical relations, obtained considering the displacement of a beam element of length $dx$:

$$\sin \theta = \frac{w'}{1+e}, \quad \cos \theta = \frac{1+w'}{1+e}. \qquad (21)$$

The end of the derivation needs the inextensibility assumption (i), that must be considered differently depending on the symmetry of the lamination.

### 2.2.1. The case of a symmetric lamination

If the lamination is symmetric (Figs. 3(a,c)), as discussed earlier, choosing the origin of the z-coordinate in the middle of the cross section leads to $B = 0$, canceling the material axial/bending coupling in Eqs. (7a,b) and enforcing a neutral line at $z = 0$ onto which $e = 0$ in linear bending. In the nonlinear case, $N(t)$ is non-zero because of the geometrical nonlinearities due to Eq. (19) and $e$ is consequently not strictly zero (see Eq. (7a)). However, it is reasonable to consider $e$ small such that $1+e \approx 1$ and, equivalently, that the beam is inextensible, yielding with Eq. (21) to:

$$1 + u' = \sqrt{1 - w'^2}. \qquad (22)$$

This last relationship allows $u$ to be eliminated in favour of $w$ in Eq. (20). Finally, considering Eq. (16) with $B = 0$ to eliminate $M$ in favour of $w$ in Eqs. (20) and subsequently performing the third order Taylor expansions:

$$\theta = \arcsin w' \approx w' + \frac{w'^3}{6}, \quad \tan \theta \approx w' + \frac{w'^3}{2}, \quad \frac{1}{\cos \theta} = \frac{1}{\sqrt{1-w'^2}} \approx 1 + \frac{w'^2}{2}, \quad u' \approx -\frac{w'^2}{2}, \qquad (23)$$

one obtains the following Partial Differential Equations (PDE) of motion:

$$\left\{ \begin{array}{l} m(x)\ddot{w} + [\hat{D}(x)w'']'' + (w'[w'w''\hat{D}(x)]')' + \left[\frac{w'}{2}\int_{x}^{L} m(\xi)\left(\int_{\xi}^{0} w'^2 d\eta\right) d\xi\right]' \\ \qquad + \sum_{p \in P} \hat{\Theta}_p \left[\Delta_p(x)\left(1+\frac{w'^2}{2}\right)\right]' V_p + \left[(w' + \frac{w'^3}{2}) \int_{x}^{L} m(\xi) d\xi\right]' \ddot{u}_0 = -m(x)\ddot{w}_0, \end{array} \right. \qquad (24a)$$

$$Q_p = \hat{\Theta}_p \left[ w' + \frac{w'^3}{6} \right]_{-}^{x+} + \hat{C}_p V_p - \frac{\bar{\Xi}_p L_p}{A} \bar{N}_{-}^{x+} + \sum_{j \neq p} C_{pj} V_j, \qquad (24b)$$

with $\Delta_p(x) = \delta(x - x_+) - \delta(x - x_-)$ with $\delta$ the Dirac delta function.

The above equations are the direct extension of the model of [37,42] to the laminated and piezoelectric case with parametric forcing. The first line of Eq. (24a) takes into account the laminated structure of the beam, through coefficients $\hat{D}(x)$ and $m(x)$ which are piecewise constant. As in the case of a homogeneous beam [37,42], the geometrical nonlinearities classically appear as a nonlinear curvature term (the second term) and a nonlinear inertial term (the last term), the latter being the consequence of the condensation of the axial motion $u(x,t)$ into the transverse motion. Then, the second line shows the axial acceleration $\ddot{u}_0(t)$ multiplied by $w'$, characteristic of a parametric forcing term. It also includes the direct excitation in its second member (proportional to $\ddot{w}_0(t)$) as well as the piezoelectric coupling proportional to voltage $V_p(t)$ in the $p$th piezoelectric patch. The associated electric charge $Q_p(t)$ (Eq. (24b)) classically includes an electromechanical coupling term proportional to $\hat{\Theta}_p$ and a purely capacitive one. Some additional and less standard terms, shown with an *underline*, are a direct consequence of the non-zero value of $e$ and $N$.

The electrical part of Eqs. (24a,b) may be simplified by considering colocalised piezoelectric patches as in Fig. 4. In the series connection, the charge and voltage viewed from the electric circuit are $Q = Q_1 = Q_2$ and $V = V_1 + V_2$, and the polarisation of the two patches must be opposed to obtain a coupling with the bending motion, such that $\bar{e}_{31}^1 = -\bar{e}_{31}^2$. In the parallel connection,
<figure><img src="image_8.png" /><figcaption>Fig. 4. Possible connections between two colocalised piezoelectric patches and the shunt resistor. (a) Series connection, (b) Parallel connection.</figcaption></figure>

$Q = Q_1 - Q_2$, $V = V_1 = -V_2$ and $\bar{e}_{31}^1 = \bar{e}_{31}^2$. In both cases, the additional term in Eq. (24b) vanishes, such that one obtains:

$$ \left\{ \begin{aligned} & m(x)\ddot{w} + [\hat{D}(x)w'']'' + (w'[\hat{D}(x)w'w'')]' + \left[ \frac{w'}{2} \int_x^L m(\xi) \left( \int_\xi^0 w'^2 d\eta \right) d\xi \right]' \\ & \qquad + \hat{\Theta} \left[ \Delta(x) \left(1 + \frac{w'^2}{2}\right) \right]' V + \left[ \left(w' + \frac{w'^3}{2}\right) \int_x^L m(\xi) d\xi \right]' \ddot{u}_0 = -m(x)\ddot{w}_0, \end{aligned} \right. \quad (25a) $$

$$ Q = \hat{\Theta} \left[ w' + \frac{w'^3}{6} \right]_{x_-}^{x_+} + CV, \quad (25b) $$

with, for the series connection $\hat{\Theta} = \Theta_0$ and $C = C_0/2$, for the parallel connection $\hat{\Theta} = 2\Theta_0$ and $C = 2C_0$, where

$$ \Theta_0 = \frac{b_p(h_b + h_p)}{2} \bar{e}_{31}, \quad C_0 = \epsilon_{33} \frac{b_p L_p}{h_p}, \quad (26) $$

the constants of one of the layers, with $(h_b, h_p)$ the thicknesses of the elastic and piezoelectric layers, respectively (see Fig. 3(c)) and $b_p$ the width of the piezoelectric layers. One may remark that Eq. (25b) may also be derived directly from Eq. (10).

It is interesting to notice that the geometrical nonlinearities add cubic nonlinear piezoelectric terms of the form "$w^2V/2$" in Eq. (25a) and "$w^3/6$" in Eq. (25b), of coefficients different in both equations. They look similar to terms obtained in the case of axial/bending coupling nonlinearities, of the form "$wV$" and "$w^2/2$" (see [50]), but are different because of the order of the nonlinearity and the coefficients. The particular value of the coefficients is explained by the existence of an electromechanical potential energy $P_c$ (see Eq. (51)).

### 2.2.2. The case of an asymmetric lamination

We consider here the case of a non-symmetric lamination, for which the term $BN/A$ may not vanish in Eq. (16). Indeed, if the beam has a uniform cross section along the x-direction (Fig. 3(b)), choosing $z_n$ (Eq. (15)) as the origin of the z-coordinates leads to $B = 0$ for all $x \in [0, L]$ such that Eq. (25a) is still exact in this case.

On the contrary, if the beam is composed of several parts of different lamination (Fig. 3(d)), one can choose the origin of the z-coordinate such that $B = 0$ in parts I and III, but one will have $B \neq 0$ in part II. If Eq. (19) is considered to eliminate $N$ in (16), one would obtain a differential equation in $M$ such that its direct substitution in Eq. (20) to obtain a single equation in $w(x, t)$ would not be possible. This break of symmetry would also add some quadratic nonlinearities in the problem since $N$ is quadratic in $w$ (see Eq. (19)), whereas Eqs. (25) in the symmetric case are purely cubic (see the Appendix of [52] for further details).

For the electric part, the same effect arises in the form of the term proportional to $\tilde{N}_x^{\lambda+}$ in Eq. (25b), leading to a slight quadratic nonlinear effect in the piezoelectric charge, due to the geometrical nonlinearities. This effect is always present in the case of an asymmetric lamination, such that a part of the electric charge $Q_p(t)$ comes from $N(t)$.

One may however think of neglecting these effect of non-zero $N(t)$, to arrive at the same governing equations (25a,b) with $\hat{D}$, $\hat{\Theta}$ and $\hat{C}$ obtained with Eqs. (18), whatever be the lamination.

## 2.3. Reduced order model

### 2.3.1. Modal expansion

The electromechanical problem (25a,b) is discretised with the help of the following single mode modal expansion:

$$ w(x, t) \approx \Phi_l(x)q(t), \quad (27) $$

with $q(t)$ the modal coordinate and $\Phi_l(x)$ the $i$th mode shape. For the sake of generality, a more complete modal expansion, not needed in the remaining of this article, is given in Appendix B. Indeed, only resonant excitations in the vicinity of the $i$th resonance will be considered, neglecting any nonlinear mode coupling. This approximation will be validated in Section 3. The eigenmodes are the free vibration eigen-solutions $(\omega_l, \Phi_l(x))$ of Eq. (25a) linearised and considered in short-circuit ($V = 0$), that reads:

$$ [\hat{D}(x)\Phi_l(x)'' - \omega_l^2 m(x)\Phi_l(x) = 0, \quad (28) $$
where $(\omega_i, \Phi_i(x))$ are the ith natural frequency (in [rad/s]) and mode shape in short circuit. The modes are normalised as follows:

$$ \int_{0}^{L} m(x) \Phi_{i}(x)^{2} dx = M_{i}, \quad \forall i = 1, \dots, N, \qquad (29) $$

with $M_i$ the ith modal mass, expressed in [kg], such that $\Phi_i(x)$ has no dimension and $q_i(t)$ is expressed in [m]. Notice that the value of $M_i$ depends on the scaling of the mode shape $\Phi_i$.

Then, Eq. (27) is introduced into Eq. (25a). Multiplying by $\Phi_i(x)$ both members of the equation and integrating the result over the length of the beam, one obtains:

$$ \left\{ \begin{aligned} & \ddot{q} + 2\xi\omega_i\dot{q} + \omega_i^2 q + \Gamma q^3 + \Pi(q\ddot{q} + \dot{q}^2)q - \frac{\chi}{M_i}V - \frac{\bar{\chi}}{2M_i}q^2V + (\Lambda + \bar{\Lambda}q^2)q\dot{u}_0 = -S\dot{w}_0, \\ & Q = \chi q + \frac{\bar{\chi}}{6}q^3 + CV. \end{aligned} \right. \qquad (30a) \qquad (30b) $$

In the above equations, $(\Gamma, \Pi)$ are the coefficients of the cubic geometrically nonlinear curvature and inertia terms, $(\chi, \bar{\chi})$ are the coefficients of the linear and nonlinear piezoelectric coupling terms and $(S, \Lambda)$ are correction factors of the direct and parametric forcing that take into account the considered mode shape $\Phi_i$. All those coefficients are defined in Appendix A (Eqs. (A.1)) and their value depend on the scaling of the mode shape $\Phi_i$. $S$ has no units and $\Lambda$ has unit [m⁻¹] (see Eqs. (A.1) and (B.9)). In addition, a linear viscous damping term has been added, with $\xi$ the modal damping factor.

Considering an electronic circuit reduced to a simple resistor $R$, Ohm's law leads to $V = -R\dot{Q}$. In addition, it may be convenient to rewrite Eq. (30a) by eliminating $V$ with the help of Eq. (30b), thus selecting the electric charge $Q(t)$ as main electrical unknown. Then, not considering the nonlinear parametric forcing term $\bar{\Lambda}q^3\dot{u}_0$ (see Appendix C) leads to:

$$ \left\{ \begin{aligned} & \ddot{q} + 2\xi\omega_i\dot{q} + \omega_i^2 q + \Gamma q^3 + \Pi(q\ddot{q} + \dot{q}^2)q + \frac{\chi R}{M_i C} \dot{Q} - \frac{\bar{\chi} R}{2M_i C} q^2 Q + \Lambda q\dot{u}_0 = -S\dot{w}_0, \\ & \tau_e \dot{Q} + Q - \chi q - \frac{\bar{\chi}}{6} q^3 = 0, \end{aligned} \right. \qquad (31a) \qquad (31b) $$

where $\hat{\omega}_i$ is the natural frequency in open circuit (when $Q=0$), $\hat{\Gamma}$ the associated cubic coefficient and $\tau_e$ the electrical circuit time constant, defined as:

$$ \hat{\omega}_i = \omega_i + \frac{\chi^2}{M_i C}, \quad \hat{\Gamma} = \Gamma + \frac{2\chi\bar{\chi}}{3M_i C}, \quad \tau_e = RC. \qquad (32) $$

Since in our case the electronic circuit is a simple resistor, another way of writing Eq. (31a) is simply to replace $V$ in Eq. (30a) by $-R\dot{Q}$, yielding to the following equation, equivalent to Eq. (31a):

$$ \ddot{q} + 2\xi\omega_i\dot{q} + \omega_i^2 q + \Gamma q^3 + \Pi(q\ddot{q} + \dot{q}^2)q + \frac{\chi R}{M_i} \dot{Q} + \frac{\bar{\chi} R}{2M_i} q^2 \dot{Q} + \Lambda q\dot{u}_0 = -S\dot{w}_0, \qquad (33) $$

The above electro-mechanical model of the system (Eqs. (31a,b) or Eqs. (33), (31b)) is equivalent to a classical lumped parameter model for which the cantilever beam is replaced by an equivalent mass/spring system of mass $M_i$ and stiffness $K_i = \omega_i^2 M_i$ in short-circuit. In this case, $q(t)$ is the displacement of the mass, that depends on the scaling of the mode shape $\Phi_i$. Choosing $\Phi_1(L) = 1$ would give $q(t) = w(L, t)$ (Eq. (27)), such that the displacement of the equivalent mass $M_i$ would be the one of the tip of the cantilever.

### 2.3.2. Dimensionless problem

It is convenient to render Eqs. (31a,b) and (33) dimensionless to identify the independent parameters of the problem. Therefore, the following dimensionless variables, denoted with an overbar, are introduced:

$$ \bar{\tau} = \omega_i t, \quad \hat{\omega}_i = \frac{\hat{\omega}_i}{\omega_i}, \quad \bar{q} = \frac{q}{L}, \quad \bar{Q} = \frac{Q}{L\omega_i\sqrt{M_i C}}, \quad \kappa = \frac{\chi}{\omega_i\sqrt{M_i C}}, \quad \bar{\kappa} = \frac{\bar{\chi} L^2}{\omega_i\sqrt{M_i C}}, \qquad (34a) $$

$$ \{\bar{\Gamma}, \hat{\Gamma}\} = \frac{\{\Gamma, \hat{\Gamma}\} L^2}{\omega_i^2}, \quad \bar{\Pi} = \Pi L^2, \quad \bar{\Lambda} = L\Lambda, \quad \bar{\delta}(t) = \bar{\Lambda} \frac{\dot{\bar{u}}_0(t)}{L\omega_i^2}, \quad \bar{F}(t) = S \frac{\dot{\bar{w}}_0(t)}{L\omega_i^2}, \quad \bar{\tau}_e = \omega_i \tau_e = \omega_i RC, \qquad (34b) $$

where the dimensionless coefficients $(\bar{\Gamma}, \bar{\Pi}, \bar{\Lambda}, \kappa, \bar{\kappa}, S)$ are written in Appendix B (Eqs. (B.9)). The dimensionless equations are then:

$$ \left\{ \begin{aligned} & \ddot{q} + 2\xi\dot{q} + \hat{\omega}_i^2 \bar{q} + \hat{\Gamma} \bar{q}^3 + \bar{\Pi}(\bar{q}\ddot{q} + \dot{q}^2)\bar{q} - \kappa \bar{Q} - \frac{\bar{\kappa}}{2} \bar{q}^2 \bar{Q} + \bar{\delta} \bar{q} = - \bar{F}, \\ & \ddot{\bar{q}} + 2\xi \dot{\bar{q}} + \bar{q} + \bar{\Gamma} \bar{q}^3 + \bar{\Pi}(\bar{q}\ddot{q} + \dot{q}^2)\bar{q} + \kappa \bar{\tau}_e \bar{Q} + \frac{\bar{\kappa} \bar{\tau}_e}{2} \bar{q}^2 \bar{Q} + \bar{\delta} \bar{q} = - \bar{F}, \end{aligned} \right. \qquad (35a) \qquad (35b) $$

$$ \bar{\tau}_e \dot{Q} + \bar{Q} - \kappa \bar{q} - \frac{\bar{\kappa}}{6} \bar{q}^3 = 0. \qquad (35c) $$

$$ \bar{\tau}_e \dot{Q} + \bar{Q} - \kappa \bar{q} - \frac{\bar{\kappa}}{6} \bar{q}^3 = 0. \qquad (35c) $$

These equations were also obtained in [20]. They bring out parameter $\kappa$, which is the classical electromechanical modal coupling factor often used as a performance indicator [67,71,72]. It relates the natural frequencies in short and open circuits (see Eq. (32)).
<table><caption>**Table 1**
Numerical values of the material parameters, with $\epsilon_0 = 8.854 \cdot 10^{-12}$ F/m the free space permittivity [43,73].</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Steel (elastic parts)</th><th colspan="2">PIC151 (piezoelectric parts)</th></tr><tr><th>ρ<sub>b</sub></th><th>6850</th><th>ρ<sub>p</sub></th><th>8500</th></tr></thead><tbody><tr><td>Density</td><td>kg/m<sup>3</sup></td><td>ρ<sub>b</sub></td><td>6850</td><td>ρ<sub>p</sub></td><td>8500</td></tr><tr><td>Young modulus</td><td>GPa</td><td>Y<sub>b</sub></td><td>210</td><td>Y<sub>p</sub></td><td>66.7</td></tr><tr><td>Piezoelectric constant</td><td>C/m<sup>2</sup></td><td>–</td><td>–</td><td>ε&#773;<sub>31</sub></td><td>-14</td></tr><tr><td>Dielectric constant</td><td>F/m</td><td>–</td><td>–</td><td>ε<sub>33</sub></td><td>2068ε<sub>0</sub></td></tr></tbody></table>

<table><caption>**Table 2**
Values of the reduced order model dimensionless parameters for the three configurations (a,b,c) under study. In all cases, $h_p = 1.01h_b$ for the steel/PZT couple of materials (see [48]).</caption><thead><tr><th></th><th>L<sub>p</sub></th><th>κ</th><th>κ&#773;</th><th>Γ&#773;</th><th>Γ&#773;&#773;</th><th>Λ&#773;</th><th>Ξ&#773;</th><th>S</th><th>Φ<sub>1</sub>(1)</th><th>Φ'<sub>1</sub>(1)</th></tr></thead><tbody><tr><td>Single layer beam</td><td>–</td><td>–</td><td>–</td><td>3.271</td><td>4.597</td><td>-1.571</td><td>-4.043</td><td>0.783</td><td>2</td><td>2.753</td></tr><tr><td>Config. (a)</td><td>0.635L</td><td>0.316</td><td>4.054κ</td><td>2.423</td><td>3.313</td><td>-1.811</td><td>-3.357</td><td>1.159</td><td>1.833</td><td>3.034</td></tr><tr><td>Config. (b)</td><td>L</td><td>0.222</td><td>3.352κ</td><td>1.447</td><td>2.033</td><td>-1.571</td><td>-1.788</td><td>1.778</td><td>1.330</td><td>1.831</td></tr><tr><td>Config. (c)</td><td>0.6L</td><td>0.254</td><td>4.941κ</td><td>2.938</td><td>4.060</td><td>-1.703</td><td>-3.864</td><td>0.958</td><td>1.949</td><td>3.004</td></tr></tbody></table>

<figure><img src="image_9.png" /></figure>

and [60]), as well as the cubic coefficient, as follows:

$$ \kappa^2 = \frac{\omega_i^2 - \omega_t^2}{\omega_i^2} \Rightarrow \hat{\omega}_i^2 = 1 + \kappa^2, \quad \hat{\Gamma} = \bar{\Gamma} + \frac{2\kappa\bar{\kappa}}{3}. \qquad (36) $$

Because of the geometrical nonlinearities, a nonlinear electromechanical coupling factor $\bar{\kappa}$ also appears in the equations. A complete dimensionless treatment from the initial PDE (25) is possible and proposed in Appendix B.

### 2.3.3. Practical configurations under study

To decrease the number of free parameters for optimisation purpose, some choices are made, without loss of generality. We focus on the piezoelectric beam configurations sketched in Fig. 3. The elastic and piezoelectric parts are made in standard steel and PIC151 from PI Ceramic GmbH, respectively (Table 1). Then, the placements and dimensions of the piezoelectric patches are chosen to maximise the coupling factor $\kappa$, as proposed in [48]. Qualitatively, as it can be guessed from Eq. (A.1g), one has to maximise the mode shape slope difference at the ends of the piezoelectric patch: $\Phi'_i(x_+) - \Phi'_i(x_-)$. The focus is also put on the first mode of vibration (i = 1) for the remainder of the paper. Therefore, it leads to $x_- = 0$, meaning aligning the patch to the clamped end of the beam. Then, since the added inertia/stiffness of the piezoelectric patches are taken into account, the thickness $h_p/h_b$ and length $L_p/L_b$ are optimised following the results of [48]. As a consequence, the three configurations of Table 2 are considered, with the numerical values of the coefficients obtained with Eqs. (B.9). For this, $\Phi_i(x)$ is written in terms of sin, cos, sinh and cosh functions, by assembling the exact solutions on the different parts of the beam, fully taking into account the piecewise nature of $m(x)$ and $\hat{D}(x)$, as described in the Appendix of [48]. It must also be noted that the values of the coefficients in Table 2 depend on the scaling of the modes, indicated in the last column with the value of $\Phi_1(L)$ and that has been chosen such that the dimensionless modal mass $\bar{m}_1 = 1$ (see Eq. (B.6)).

Considering the values of the coefficients of the model in Table 2, one can observe that they share the same order of magnitude, whatever be the considered configuration, which further justifies the dimensionless writing of the equations.

## 3. Single mode dynamics and solution

In this section, three different methods of solving Eqs. (35) are studied and compared. Two are analytical: the Method of Multiple Scales (MMS) and the Harmonic Balance Method (HBM) with $H = 1$ harmonic (HBM-H1). These methods of solution are fully described in [74,75] for several direct and parametric nonlinear oscillators, and HBM-H1 and MMS are used in [23–25] for similar PEH. Finally, numerical continuation of periodic solutions is used as a third method for reference. In practice, the software MANLAB is used, based on a HBM with an arbitrary number of harmonics and an asymptotic numerical method for efficient continuation [76,77]. It allows also the computation of bifurcations and the stability of the periodic solutions [78]. This reference numerical method will be denoted as HBM-H30 in the following, as a number $H = 30$ of harmonics was found largely sufficient to reach convergence.

The comparison between the results of the different methods is carried out on a simple cantilever beam model, without piezoelectric patches, adequate to highlight the limit of the analytical solutions. Moreover, it also enables to validate the result of the single-mode truncation, by comparison to a finite-element solving of the PDE found in [40,41], called FEHBM in the following.
### 3.1. Direct excitation

The oscillations of the problem (35a,b,c) under direct only harmonic excitation, i.e. without parametric drive ($\delta = 0$), and with $\bar{F}(t) = F \cos \bar{\Omega}\bar{t}$ with $F$ the amplitude of the forcing and $\bar{\Omega} = \Omega/\omega_i$ the excitation frequency ($\bar{\Omega}$ is dimensionless whereas $\Omega$ is expressed in [rad/s]), are addressed here. The analytical expression obtained with both the HBM-H1 and the MMS (restricted to the first order) are both:

$$ \bar{q}(\bar{t}) = a \cos(\bar{\Omega}\bar{t} + \alpha), \quad \bar{Q}(t) = b \cos(\bar{\Omega}\bar{t} + \beta), \qquad (37) $$

where $(a, \alpha)$ and $(b, \beta)$ are the amplitude/phase pairs of $q(t)$ and $Q(t)$.

Introducing the above ansatz in Eqs. (35b,c) and applying a standard HBM procedure restricted to one harmonic leads to the following equations:

$$ -(T_1 + T_3 a^2)a = F \cos \alpha, \quad (T_2 + T_4 a^2)a = F \sin \alpha, \qquad (38) $$

that can be combined to give:

$$ (T_1 a + T_3 a^3)^2 + (T_2 a + T_4 a^3)^2 = F^2, \quad \tan(\alpha) = -\frac{T_2 + T_4 a^2}{T_1 + T_3 a^2}, \qquad (39) $$

where

$$ T_1 = 1 - \bar{\Omega}^2 + \kappa^2 \zeta_1, \quad T_2 = 2\xi \bar{\Omega} + \kappa^2 \zeta_2, \quad \zeta_1 = \frac{(\bar{\Omega} \bar{\tau}_e)^2}{1 + (\bar{\Omega} \bar{\tau}_e)^2}, \qquad (40a) $$

$$ T_3 = \frac{3\bar{\Gamma} - 2\bar{\Pi}\bar{\Omega}^2}{4} + \frac{\kappa\bar{\kappa}}{2}\zeta_1, \quad T_4 = \frac{\kappa\bar{\kappa}}{4}\zeta_2, \quad \zeta_2 = \frac{\bar{\Omega}\bar{\tau}_e}{1 + (\bar{\Omega}\bar{\tau}_e)^2}. \qquad (40b) $$

It has been verified that starting from the alternate formulation Eqs. (35a,c) instead of (35b,c), one obtains the same results. The first equation above is a third order polynomial in $a^2$ which has no analytical solution in the general case.

Similarly, starting from Eq. (35a,b,c) and applying the Method of Multiple Scales, similar relationships than Eqs. (39) are obtained for $(a^2, \alpha)$, with different values for the coefficients $(T_i, \zeta_j)$, $i = 1, ..., 4$, $j = 1, 2$, which are also slightly different if Eqs. (35a,c) or Eqs. (35b,c) are considered. As it will be seen in Section 3.3, the former is better suited for the PEH study around the short-circuit conditions and the latter around the open-circuit conditions, especially if $\kappa$ is large. With Eq. (35b), the coefficients are denoted by $T_i = \bar{T}_i^b$ and $\zeta_j = \bar{\zeta}_j^a$ and take the form:

$$ \bar{T}_1^b = 2(1 - \bar{\Omega}) + \kappa^2 \bar{\zeta}_1^b, \quad \bar{T}_2^b = 2\xi + \kappa^2 \bar{\zeta}_2^b, \quad \bar{\zeta}_1^b = \frac{\bar{\tau}_e^2}{1 + \bar{\tau}_e^2}. \qquad (41a) $$

$$ \bar{T}_3^b = \frac{3\bar{\Gamma} - 2\bar{\Pi}}{4} + \frac{\kappa\bar{\kappa}}{2}\bar{\zeta}_1^b, \quad \bar{T}_4^b = \frac{\kappa\bar{\kappa}}{4}\bar{\zeta}_2^b, \quad \bar{\zeta}_2^b = \frac{\bar{\tau}_e}{1 + \bar{\tau}_e^2}. \qquad (41b) $$

Besides, Eq. (35a) leads to write $T_i = \bar{T}_i^a$ and $\zeta_j = \bar{\zeta}_j^a$, with

$$ \bar{T}_1^a = 2(\hat{\omega}_i - \bar{\Omega})\hat{\omega}_i - \kappa^2 \bar{\zeta}_1^a, \quad \bar{T}_2^a = 2\xi\hat{\omega}_i + \kappa^2 \bar{\zeta}_2^a, \quad \bar{\zeta}_1^a = \frac{1}{1 + \bar{\tau}_e^2 \hat{\omega}_i^2}. \qquad (42a) $$

$$ \bar{T}_3^a = \frac{3\bar{\Gamma} - 2\bar{\Pi}\hat{\omega}_i^2}{4} + \frac{\kappa\bar{\kappa}\bar{\tau}_e\hat{\omega}_i}{2}\bar{\zeta}_2^a, \quad \bar{T}_4^a = \frac{\kappa\bar{\kappa}}{4}\bar{\zeta}_2^a, \quad \bar{\zeta}_2^a = \frac{\bar{\tau}_e\hat{\omega}_i}{1 + \bar{\tau}_e^2\hat{\omega}_i^2}. \qquad (42b) $$

### 3.2. Parametric excitation

Under parametric excitation only ($F(t) = 0$), the main resonance is obtained by forcing at twice the frequency $\Omega$ of the response, such that $\delta(\bar{t}) = \delta \cos(2\bar{\Omega}\bar{t})$ with $\delta$ the amplitude of the parametric forcing. The same solution (37) is obtained for both the HBM-H1 and the MMS for the PEH dynamics (Eqs. (35)). Both methods lead to:

$$ -2(T_1 + T_3 a^2) = \delta \cos 2\alpha, \quad 2(T_2 + T_4 a^2) = \delta \sin 2\alpha, \qquad (43) $$

which gives:

$$ (T_3^2 + T_4^2)a^4 + 2(T_1 T_3 + T_2 T_4)a^2 + (T_1^2 + T_2^2) - \delta^2/4 = 0, \quad \tan(2\alpha) = -\frac{T_2 + T_4 a^2}{T_1 + T_3 a^2}. \qquad (44) $$

Solving the first equation for $a$ leads to:

$$ a^2 = \frac{-2(T_1 T_3 + T_2 T_4) \pm \sqrt{\delta^2(T_3^2 + T_4^2) - (2T_1 T_4 + 2T_2 T_3)^2}}{2(T_3^2 + T_4^2)}. \qquad (45) $$

In the above equations, coefficients $T_i$ are defined by Eqs. (40a,b) for the HBM-H1 and by Eqs. (41a,b) or (42a,b) for the MMS.
As a well known result (see [75,79,80]), a non-zero response under parametric forcing is obtained only above a given threshold $\delta_{cr}$, i.e for $\delta > \delta_{cr}$. This threshold is obtained by enforcing the term under the square root of Eq. (45) to be positive, such that:

$$ \delta_{cr} = 2 \frac{T_1 T_4 + T_2 T_3}{\sqrt{T_3^2 + T_4^2}} \qquad (46) $$

As this expression depends on $\bar{\Omega}$ (see Eqs. (40a,b), (41a,b) and (42a,b)), an additional constraint is needed. Numerically, it is shown that the birth of the parametric response, for $\delta = \delta_{cr}$ and $a = 0$, appears on the frequency-axis at a frequency $\bar{\Omega}_{cr}$ very close to the resonance frequency of the linear and direct response of the system, which is the recommended choice in practice.<fnref n="2" /> The analyse of $\delta_{cr}$ will be considered in Section 4.4.

The parametric threshold takes particular values in the two extreme cases of short circuit (SC) and open-circuit (OC). The SC condition is obtained with a zero resistance: $R = 0 \Rightarrow \bar{\tau}_e = 0$, such that $\zeta_1 = \zeta_2 = 0 \Rightarrow T_4 = 0 \Rightarrow \delta_{cr} = 2T_2 = 4\xi\bar{\Omega}$ and $\bar{\Omega}_{cr} = 1$, the natural frequency in SC. The OC condition is obtained with an infinite resistance: $R \to +\infty \Rightarrow \bar{\tau}_e \to +\infty$, such that $\zeta_1 = 1$, $\zeta_2 = 0 \Rightarrow T_4 = 0$, leading to the same relationship $\delta_{cr} = 2T_2 = 4\xi\bar{\Omega}$. The difference with the SC condition lies in the value of $\bar{\Omega}_{cr} = \sqrt{1 + \kappa^2} = \hat{\omega}_i$. In summary:

$$ \delta_{cr}^{\text{sc}} = 4\xi, \quad \delta_{cr}^{\text{oc}} = 4\xi\sqrt{1 + \kappa^2}. \qquad (47a,b) $$

Those values where validated by comparison to the reference numerical HBM-H30 method.

With the MMS, the SC conditions $\bar{\tau}_e = 0$ in Eqs. (41) and the OC condition $\bar{\tau}_e \to +\infty$ in Eqs. (42) leads to the same values of $\delta_{cr}^{\text{sc}}$ and $\delta_{cr}^{\text{oc}}$ than with the HBM-H1 (Eqs. (47a,b)). However, the SC conditions in Eqs. (42) and the OC condition in Eqs. (41) lead to slightly different results, showing that the validity domain of the MMS method is limited to oscillations in the vicinity of the linear natural frequency of the nonlinear oscillator under concern, which are 1 for Eq. (35a) and $\hat{\omega}_i$ for Eq. (35b). This also confirms that the validity domain of the HBM-H1 is often larger than the one of the MMS, a result already pointed out in the literature, in [80] for instance.

### 3.3. Backbone curves in short- and open-circuit

The free periodic oscillations of Eqs. (35a,b,c), obtained if the system is conservative, are addressed in this section. Both structural and electrical damping must be canceled, leading to $\xi = 0$ and either $\bar{\tau}_e = 0$ (SC) or $\bar{\tau}_e \to +\infty$ (OC), as well as the forcing ($F = \delta = 0$). Considering the HBM-H1 solution (Eqs. (40a,b)), in both SC and OC cases, $T_2 = T_4 = 0$ such that Eqs (39) or (44) degenerate to $a = 0$ or $(T_3a^2 + T_1)^2 = 0$. Two backbone curves $a^2 = -T_1/T_3 = f(\bar{\Omega})$ are obtained in the two extreme SC and OC cases, that read:

$$ \text{SC: } a^2 = \frac{4(\bar{\Omega}^2 - 1)}{3\bar{\Gamma} - 2\bar{\Pi}\bar{\Omega}^2}, \quad \text{OC: } a^2 = \frac{4(\bar{\Omega}^2 - \hat{\omega}_i^2)}{3\hat{\Gamma} - 2\hat{\Pi}\bar{\Omega}^2} = \frac{4(\bar{\Omega}^2 - 1 - \kappa^2)}{3\bar{\Gamma} + 2\kappa\bar{\kappa} - 2\hat{\Pi}\bar{\Omega}^2}. \qquad (48) $$

For the MMS solution, following the same procedure in SC (starting from Eq. (35b) and taking $\bar{\tau}_e = 0$) and in OC (starting from Eq. (35a) and taking $\bar{\tau}_e \to +\infty$), leads to:

$$ \text{SC: } a^2 = \frac{8(\bar{\Omega} - 1)}{3\bar{\Gamma} - 2\bar{\Pi}}, \quad \text{OC: } a^2 = \frac{8(\bar{\Omega} - \hat{\omega}_i)\hat{\omega}_i}{3\hat{\Gamma} - 2\bar{\Pi}\hat{\omega}_i^2} = \frac{8(\bar{\Omega}\sqrt{1 + \kappa^2} - 1 - \kappa^2)}{3\bar{\Gamma} + 2\kappa\bar{\kappa} - 2\hat{\Pi}(1 + \kappa^2)}. \qquad (49) $$

The same situation than observed for the computation of $\delta_{cr}$ in the previous Section 3.2 is obtained here: considering Eq. (35a) in SC and Eq. (35b) in OC leads to erroneous results. In particular, the foot of the backbone curves, for $a = 0$, do not coincide with the respective natural frequencies in SC ($\bar{\Omega} = 1$) and in OC ($\bar{\Omega} = \hat{\omega}_i = \sqrt{1 + \kappa^2}$).

It is interesting to notice that those backbone curves can also be obtained under a phase resonance condition, as proved in [81]. In the case of a direct excitation, if the response is in phase quadrature with the excitation, the forcing term exactly cancels the damping term at any time $t$, such that the system behaves as if it was in free conservative oscillations, i.e. on the backbone curve. In the present case, having $a^2 = -T_1/T_3$ on the backbone curve (Eqs. (48) or (49)) is equivalent to having a phase lag of $\alpha = +\pi/2$, thanks to Eq. (39).<fnref n="3" /> Analogously, the same situation occurs in parametric excitation for $2\alpha = \pi/2 \mod 2\pi$, obtained for $\alpha = \pi/4$ or $\alpha = -3\pi/4$. Both values are possible and lead to two distinct solutions with the same amplitude (see [82]).

Eqs. (48) and (49) show that the hardening/softening behaviour of the backbone curves is the result of the combined values of $\bar{\Gamma}$ and $\bar{\Pi}$, but may also be changed by the piezoelectric coupling, through $(\kappa, \bar{\kappa})$ in the OC case. This effect is very slight in practice. Moreover, in the case of the HBM-H1 method, there is a singularity for $\Omega = \sqrt{3\bar{\Gamma}/2\bar{\Pi}}$ (in SC) or for $\Omega = \sqrt{3\hat{\Gamma}/2\hat{\Pi}}$ (in OC), already

<footnote n="2">² We also tested to enforce $a = 0$ in Eq. (45), giving $T_1T_3 + T_2T_4 = 0 \Rightarrow \bar{\Omega} = \bar{\Omega}_{cr}$, but it was found that the resonance frequency choice gives results closer to the reference numerical solution for large values of $\kappa$ and for $\bar{\tau}_e \approx 1$.</footnote>

<footnote n="3">³ Notice that the phase lag between the displacement and the direct forcing is often $-\pi/2$ [81]. The + sign obtained here comes from the - sign of the forcing in the second member of Eqs. (35a,b).</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 5.** Backbone curves of Eqs. (35b) for a single layer cantilever beam ($\Gamma = 3.271$, no piezoelectric coupling $\kappa = \bar{\kappa} = 0$). (a) Backbone curves for different values of $\bar{\Pi}$ (starting from a simple Duffing oscillator with $\bar{\Pi} = 0$ to a single layer cantilever beam with $\bar{\Pi} = 4.597$) and obtained with different methods: first order MMS (dashed line), HBM-H1 (dash-dotted line), HBM-H30 (solid line) and the geometric exact beam model FEHBM of [41] (thick solid line). (b) Deformed shape of the cantilever beam at the validity limit ($w_{H1} = 0.7L$) of the HBM-H3 model.</figcaption></figure>

observed in [83]. For the homogeneous single layer cantilever beam considered in Fig. 5, this singularity is observed for $\bar{\Omega} = 1.0331$, far beyond the validity limit of the analytical solutions, as addressed in the next section.

### 3.4. Validity limit of the models

To choose the most suitable solving method, the case of a homogeneous single layer cantilever beam is considered (the values of the parameters are gathered in the first row of Table 2). First, a backbone curve comparison between the analytical first order MMS and HBM-H1 solutions and the numerical HBM-H30 method is proposed in Fig. 5 for different values of the inertial nonlinearities (the term $\bar{\Pi}(\bar{q}\ddot{q} + \bar{q}^2)\bar{q}$ in Eq. (35a,b), by artificially changing the value of $\bar{\Pi}$, from the Duffing case ($\bar{\Pi} = 0$) to the single layer cantilever case ($\bar{\Pi} = 4.597$). For the latter, a backbone curve obtained with a finite element discretisation of the geometrically exact beam model (12) (FEHBM, i.e. without the 3rd order expansion of Eq. (25) and the single mode expansion of Eq. (27)), solved by numerical continuation and taken from [40,41], is added as a reference.

To obtain a clear physical insight of the results, the y-axis corresponds to the vibration amplitude of the first harmonic of the tip of the beam scaled by its length, $w_{H1}/L$. For the two analytical solutions, Eqs. (37), (34a), (27) and the value of $\Phi(L)$ in Table 2 leads to $w(L,t)/L = \Phi_1(L)a \cos \bar{\Omega}\bar{t} \Rightarrow w_{H1}/L = 2a$. For the HBM-H30 solution and the FEHBM, only the first harmonic is considered in the plot.

An analysis of Fig. 5 shows that the inertial nonlinearity term (proportional to $\bar{\Pi}$) has a dramatic effect on the quality of the two analytical solutions MMS and HBM-H1: for $\bar{\Pi} = 0$ (Duffing oscillator), these two solutions match the reference HBM-H30 solution on a wide range of $\bar{\Omega} \in [0, 1.05]$ and $w_{H1}/L \in [0, 0.4] \Leftrightarrow a \in [0, 0.2]$. On the contrary, for the case of the first mode of a cantilever beam ($\bar{\Pi} = 4.597$), the match occurs on a smaller amplitude range ($w_{H1}/L \in [0, 0.3]$), which corresponds to a very small frequency range $\bar{\Omega} \in [1, 1.0017]$, since the mode is less hardening.

On the contrary, the reduced order model (35a,b) (based on a 3rd. order Taylor expansion and a one mode expansion) has a wider validity range, up to an amplitude $w_{H1} \approx 0.7L$, assessed by comparing the HBM-H30 solution with the FEHBM result. This limit corresponds to $\bar{\Omega} = 1.011$ and to a rotation of the cantilever tip of about $0.3\pi$, as illustrated in Fig. 5(b) (see also Fog. 17(b)). This is mainly a result of the third-order Taylor expansions (23).

The discrepancy between the results of HBM-H1, MMS and HBM-H30 is also observed on the frequency response curves, shown for direct and parametric excitation in Fig. 6, for a response amplitude close to the HMB-H30 validity limit $w_{H1} \approx 0.7L$. The estimations of the resonance amplitude under direct excitation are similar for the three methods, with a slight difference in the resonance frequencies explained by the discrepancy of the backbone curves. On the contrary, a huge difference is observed in the case of parametric forcing, since the first order analytical MMS and HBM-H1 are unable to predict the amplitude at resonance. In particular, those analytical solutions are unbounded in amplitude and follow the backbone curve to infinity. Similar results are obtained in the case of parametrically driven Duffing oscillator in [79]. On the contrary, the threshold of parametric response is well predicted by the analytical expressions (47), which is logical since it is associated to a near zero motion amplitude $a$.

It is interesting to mention that the reduced validity range of the analytical methods (HBM-H1 and MMS) in the present case of the first mode of a cantilever beam was not noticed in previous publications on the subject since most of them consider a quadratic damping term (of the form $|\dot{q}|\dot{q}$) to take into account the drag effect of the surrounding air (see [17,18,20,84]). This kind of damping model, though very realistic [81], is not considered in the present study.

To validate the single mode expansion used in the HBM-H30 model, Fig. 7 shows a comparison with the geometrically exact beam FEHBM model [40,41], in which all modes are kept. A correct agreement is obtained, with a slight overestimation of the
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 6. Single layer cantilever beam frequency response curves obtained with three solving methods: HBM-H1 (dash-dotted line), first order MMS (dashed line) and numerical HBM-H30 (solid line) ($\bar{F} = 3.271, \bar{N} = 4.597, \kappa = \bar{\kappa} = 0$), for a damping ratio $\xi = 0.01$, for $F = 0.007$ and for $\delta = 0.042$. (a) direct excitation; (b) parametric excitation.</figcaption></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 7. Single layer cantilever beam frequency response curves comparison obtained with 2 models: the geometric exact beam model FEHBM of [41] (orange) and the truncated model (Eq. (35b)) solved with HBM-H30 (blue), for a damping ratio $\xi = 0.01$, for $F = [0.0021; 0.0057; 0.0089; 0.0158]$ and for $\delta = [0.0402; 0.0413; 0.0425; 0.0444]$. The backbone curves are indicated in dashed lines. (a) direct excitation; (b) parametric excitation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

amplitude at resonance by the HBM-H30 results in the validity range of $w_{H1} = 0.7L$ defined above. This may be explained by the one mode truncation of the HBM-H30 model (Eq. (35a,b)), that neglects the geometrically nonlinear effect of the other modes. In particular, it was shown in [81] that the linear mass proportional damping model used in the FEHBM computation can lead to a nonlinear damping effect at large amplitude, that is not captured by the one-mode model HBM-H30.

As a conclusion, it appears that the validity range of the analytical solutions HBM-H1 and MMS, of the order of $w_{H1} \approx 0.3L$ for the backbone curves, is too small for the present application of PEH since it will be shown in Section 4.4 that large amplitude of the cantilever is required to obtain a competitive energy harvesting. As a consequence, the numerical HBM-H30 numerical solving will be preferred in the optimisations of the next section.

## 4. Optimisation for energy harvesting

The purpose of this section is to analyse the performance of the PEH under direct excitation (called a DPEH) and parametric excitation (called PPEH) and optimise it, i.e. finding the values of the independent parameters of the system that guarantee a maximum value of the power dissipated in the load resistance $R$. The goals are (i) to analyse the effect of the geometrical nonlinearities on a DPEH, (ii) to investigate the performance of a PPEH and (iii) to systematically compare their performance to the reference case of a linear DPEH. Without loss of generality, all the numerical results of this section are obtained with a configuration (a) PEH (see Table 2).

### 4.1. Maximum power analysis

To highlight the harvested power, we conduct an energy assessment of the one-mode model (31a,b). By multiplying Eqs. (31a,b) respectively with $M_i\dot{q}$ and $\dot{Q}/C$ and adding the two obtained equations, one arrives to:

$$ \frac{d}{dt} [\tau_m + \nu_m + \nu_e + \nu_c] + P_d + P_h = P_e + P_p, \qquad (50) $$
with

$$ \mathcal{T}_m = \frac{1}{2} M_i \dot{q}^2, \quad \mathcal{V}_m = \frac{1}{2} M_i \dot{\omega}_i^2 q^2 + \frac{1}{4} M_i \dot{\Gamma} q^4 + \frac{1}{2} M_i \Pi \dot{q}^2 q^2, \quad \mathcal{V}_e = \frac{Q^2}{2C}, \quad \mathcal{V}_c = -\frac{\chi}{C} Q q - \frac{\tilde{\chi}}{6C} Q q^3, \quad (51a) $$

$$ P_d = 2\xi M_i \omega_i \dot{q}^2, \quad P_h = R\dot{Q}^2, \quad P_e = -M_i S \ddot{w}_0 \dot{q}, \quad P_p = -M_i \Lambda \ddot{u}_0 \dot{q}, \quad (51b) $$

which are the instantaneous mechanical kinetic and potential energies $\mathcal{T}_m$ and $\mathcal{V}_m$, the electrical potential energy $\mathcal{V}_e$, the piezoelectric coupling energy $\mathcal{V}_c$, the power dissipated by mechanical damping $P_d$, the harvested power (dissipated in the resistor) $P_h$ and the external power, direct $P_e$ and parametric $P_p$, all functions of time $t$.

We are interested in the solution of Eqs. (31a,b) under a harmonic forcing, either in direct or parametric excitations:

$$ \ddot{w}_0 = \gamma \cos \Omega t, \quad \text{or} \quad \ddot{u}_0 = \gamma \cos 2\Omega t, \quad (52) $$

where $\gamma$ is the amplitude of the prescribed base acceleration. Consequently, $q(t)$ would also be harmonic if no nonlinearities were present. In our case, the nonlinearities are responsible for the creation of harmonics such that $(q(t), Q(t))$ are T-periodic in the general case, with $T = 2\pi/\Omega$. They could be more complex (quasi-periodic or chaotic, see [33]), but this case is not considered in the present study. Then, integrating Eq. (50) over one period $T$ of the motion leads to define the average powers as:

$$ P_{\eta} = \frac{1}{T} \int_{0}^{T} P_{\eta}(t) dt, \quad \eta = h, d, e, p. \quad (53) $$

The average harvested power is then:

$$ P_h = P_e + P_p - P_d = -\frac{M_i}{T} \int_0^T (S\ddot{w}_0\dot{q} + \Lambda\ddot{u}_0q\dot{q} + 2\xi\omega_i\dot{q}^2) dt, \quad (54) $$

since the integral of the term $d/dt[-]$ in Eq. (50), linked to potential and kinetic energies, is zero over one period of the periodic motion.

### 4.1.1. Maximum power analysis in direct forcing

From Eq. (54), it is possible to estimate the maximum harvested power in direct excitation only ($\ddot{u}_0 = 0$), as proposed in [85]. Since the motion is periodic, one expands $q(t)$ in Fourier series:

$$ q(t) = a_0 + \sum_{h=1}^{H} (a_h \cos h\Omega t + \bar{a}_h \sin h\Omega t), \quad (55) $$

where $(a_0, a_h, \bar{a}_h)$, $h = 1, \dots, H$, are real numbers and $H \in \mathbb{N}^*$ is the number of harmonics retained. Then, Eqs. (54) along with the Parseval's identity leads to:

$$ P_h = \frac{M_i \Omega}{2} \underbrace{\left[ -S\gamma\bar{a}_1 - 2\xi\omega_i\Omega \sum_{h=1}^{H} h^2 (a_h^2 + \bar{a}_h^2) \right]}_{P_e} \underbrace{\left. \right\} \quad (56) $$

The maximum extracted power $P_{\max}$ is obtained by rendering $P_h$ extremal, such that:

$$ \frac{\partial P_h}{\partial a_h} = \frac{\partial P_h}{\partial \bar{a}_h} = 0, \quad \forall h = 1, \dots, H \quad \Rightarrow \quad \bar{a}_1 = -\frac{S\gamma}{4\xi\omega_i\Omega}, \quad a_1 = a_h = \bar{a}_h = 0, \quad \forall h \ge 2. \quad (57) $$

This leads to:

$$ P_{\max} = \frac{M_i(S\gamma)^2}{16\xi\omega_i\Omega}, \quad q(t) = \bar{a}_1 \sin \Omega t, \quad P_{\max} = P_d = P_e/2. \quad (58) $$

Notice that the term $M_i S^2$ is independent of the mode shape scaling (see Eqs. (29) and (A.1e)) and that it represents the equivalent mass of the energy harvester.

The above equation, already obtained in [85], is interesting since it proves several general results. It first gives an upper bound of the harvested power under direct mono-harmonic forcing, this upper bound being obtained when the external power $P_e$ is separated in two equal parts, half of it being dissipated in the mechanical damper. This shows that in the best condition, only half of the external power is harvested, the other half being lost in mechanical damping. Then, the maximum power is harvested when the displacement $q(t)$ of the resonator is in phase quadrature with the input acceleration $\ddot{w}_0(t)$. Since $\bar{a}_1 < 0$, this is equivalent to having a phase shift $\alpha = \pi/2$, as if the system were at its phase resonance in open or short circuit (see Section 3.3). Moreover, higher harmonics created by the nonlinearities must be zero to harvest the maximum power. In other words, even if the nonlinearities are here conservative (associated to the potential energy $\mathcal{V}_m$) and thus play no role in the power balance of Eq. (54), they may create harmonics of order $h \ge 2$ in the response that would increase the dissipation in the mechanical damper and reduce the harvested power (see the term $a_h^2 + \bar{a}_h^2$ in Eq. (56)). Finally, as a well known result, Eq. (58) shows that the harvested power may be increased by (i) increasing its mass or (ii) decreasing its natural frequency or (iii) decreasing its damping factor.
### 4.1.2. Maximum power analysis in parametric forcing

We now consider the case where $\dot{w} = 0$ and a parametric excitation $\ddot{u}_0$ is applied. Considering Eqs. (54) and (55) in the same way as in the previous section, and using trigonometric identities, the harvested power under parametric excitation reads:

$$P_h = \frac{M_h \Omega}{2} \left[ \underbrace{-\Lambda\gamma \left( a_1 \ddot{a}_1 + 2a_0 \ddot{a}_2 + \sum_{h=1}^{H} (a_h \ddot{a}_{h+2} - \ddot{a}_h a_{h+2}) \right)}_{P_p} - \underbrace{2\xi \omega_l \Omega \sum_{h=1}^{H} h^2 (a_h^2 + \ddot{a}_h^2)}_{P_d} \right]. \quad (59)$$

The condition for the power to be a local extremum is now written as:

$$\frac{\partial P_h}{\partial a_0} = \frac{\partial P_h}{\partial a_h} = \frac{\partial P_h}{\partial \ddot{a}_h} = 0, \quad \forall h = 1, \dots, H, \quad (60)$$

which leads to $(2H + 1)$ equations in the $(2H + 1)$ unknowns $(a_0, a_h, \ddot{a}_h)$, that can be put under the following linear system:

$$M_H a = 0 \quad (61)$$

with $a = [a_1, \ddot{a}_1, a_2, \ddot{a}_2, \dots, a_H, \ddot{a}_H]^T$ and $M_H$ a matrix that is a sole functions of the parameter $Y = \Lambda\gamma/4\xi\omega_l\Omega$, as shown in Appendix E. It is a system of linear homogeneous equations, such that if $\det(M_H) \neq 0$, the unique solution is the trivial one, $a = 0 \Rightarrow q(t) = 0 \forall t$. It corresponds to the lowest possible bound, without any motion and harvested power, giving no interesting information.

It is thus worse exploring the case for which the system is singular, verifying $\det(M_H) = 0$. If $H = 1$ (it is the single harmonic solution for $q(t)$), Eq. (61) simply reads:

$$\begin{bmatrix} 1 & Y \\ Y & 1 \end{bmatrix} \begin{bmatrix} a_1 \\ \ddot{a}_1 \end{bmatrix} = 0 \quad (62)$$

such that $\det(M_1) = 1 - Y^2 = 0$ leads to $Y = 1$ and:

$$\Lambda\gamma = 4\xi\omega_l\Omega \Leftrightarrow \delta = \delta_{cr}^{sc} = 4\xi, \quad \ddot{a}_1 = -a_1 \Rightarrow q(t) = a_1(\cos\Omega t - \sin\Omega t), \quad P_h = 0. \quad (63)$$

As a consequence, on the contrary to the DPEH case of Section 4.1.1, considering Eq. (61) with $H = 1$ harmonic only gives, in the present PPEH case, a lower bound in power: the acceleration input is equal to the parametric excitation threshold, for which naturally no power is harvested ($P_h = 0$). It is interesting to note that in this case, the phase shift between the displacement and the input acceleration is obtained as $\alpha \in \{\pi/4, -3\pi/4\}$, depending on the sign of $a_1$, as if the system were at its phase resonance in open or short circuit (see Section 3.3).

With $H \ge 2$, as explained in Appendix E, solving $\det(M_H) = 0$ with symbolic computation gives polynomials in $Y^2$ with $H - 1$ roots, equivalent to values of $Y$ that render Eq. (61) singular. The smaller value is around $Y \approx 0.91$ (instead of $Y = 1 \Leftrightarrow \delta = \delta_{cr}^{sc}$ for $H = 1$) and the other values are larger (2.9, 7.5 etc., see Fig. E.1), with no clear interpretation.

As a conclusion, in the present case of a parametric excitation, the extremal harvested power analysis, using Eq. (54), fails to give an upper bound of the harvested power, as it was the case in Section 4.1.1 for direct excitation. This may be a consequence of the fact that the interesting parametric response appears after a bifurcation, and is unbounded if no nonlinearities ($H = 1$) are considered.

### 4.2. Reference solution: a linear DPEH

Even if a linear and directly excited PEH is well known in literature [7,86,87], some of its useful characteristics concerning its performance in term of harvested power and its optimisation are reminded in this part, as a reference for the two upcoming parts about the effect of the nonlinearities and of a parametric excitation.

Considering Eqs. (35a,c) without nonlinearities and parametric forcing ($\bar{\Gamma} = \bar{I} = \bar{k} = \bar{\delta} = 0$), the exact solution is the mono-harmonic solutions $\bar{q}(t) = \Re(\hat{a}e^{j\bar{\Omega}t})$, $\bar{Q}(t) = \Re(\hat{b}e^{j\bar{\Omega}t})$ for the dimensionless displacement and electric charge, with $\hat{a}, \hat{b} \in \mathbb{C}$ their complex amplitudes, such that:

$$\hat{b} = -\frac{\kappa}{\bar{D}(\bar{\Omega})} F, \quad \hat{a} = -\frac{1 + j\tau_e \bar{\Omega}}{\bar{D}(\bar{\Omega})} F, \quad (64)$$

with $\bar{D}(\bar{\Omega}) = 1 - (1 + 2\bar{\tau}_e\xi)\bar{\Omega}^2 + j\bar{\Omega}[\bar{\tau}_e(1 + \kappa^2) + 2\xi - \bar{\tau}_e\bar{\Omega}^2]$. Then, considering the definition (51), (53) of the average harvested power $P_h$, Parseval's identity and the definition of the dimensionless variables (34), one can show that:

$$P_h = P_{\max} \cdot \bar{P}(\xi, \kappa, \bar{\tau}_e, \bar{\Omega}), \quad \bar{P} = \frac{8\xi\kappa^2\bar{\tau}_e\bar{\Omega}^2}{|\bar{D}(\bar{\Omega})|^2}, \quad (65)$$

where $\bar{P} \le 1$ is a dimensionless power, which equals 1 for $P = P_{\max}$, and $P_{\max}$ is defined in Eq. (58).

Eq. (65) shows that, apart from the maximum value $P_{\max}$ introduced in Section 4.1, the harvested power depends on only four dimensionless independent parameters: $(\xi, \kappa, \bar{\tau}_e, \bar{\Omega})$, through $\bar{P}$, which is especially considered in the following. Plotting $\bar{P}$ as a function of $\bar{\Omega}$ for several values of $\bar{\tau}_e$ from SC to OC conditions in Fig. 8(e,f), one observes that the power response of the system is resonant and that it exists optimal values of $\bar{\tau}_e$ that render $\bar{P}$ maximum. Because of this resonant behaviour, it is useful, from the power frequency response, to extract the resonant value of the power $\bar{P}_{\text{res}}(\xi, \kappa, \bar{\tau}_e) = \bar{P}(\xi, \kappa, \bar{\tau}_e, \bar{\Omega}_{\text{res}})$, where $\bar{\Omega}_{\text{res}}$ is the power
<figure><img src="image_15.png" /><figcaption>**Fig. 8.** Frequency response of a linear DPEH for several resistive loads $R$ ($\bar{\tau}_e$ takes 30 values from 0.001 (≈SC, blue curve) to 1000 (≈OC, yellow curve)), constant damping $\xi = 0.01$ and for two coupling factors ($\kappa = 0.16 < \kappa_{lim}$ for the left column, plots (a,c,e) and $\kappa = 0.316 > \kappa_{lim}$ for the right column, plots (b,d,f)). (a,b) tip amplitude of the PEH normalised by its resonant SC value $F\Phi_1(L)/(2\xi)$; (c,d) phase lag $\alpha$ between the displacement and the external acceleration; (e,f) dimensionless power $P$. The × markers show the power resonance and the red/orange curves are the optimal ones (with maximal power at resonance). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

resonance frequency. The points ($\bar{\Omega}_{res}$, $\bar{P}_{res}$) are shown in Fig. 8(e,f) by × markers. This resonant power $\bar{P}_{res}$ is thus a function of ($\xi$, $\kappa$, $\bar{\tau}_e$) only. Its evolution as a function of ($\bar{\tau}_e$, $\kappa$) is illustrated in Fig. 9(e) for $\xi = 0.01$. The obtained surface is very symmetrical with respect to $\bar{\tau}_e$ in logarithmic scale, with a mean value of $\bar{\tau}_e = 1$. It is characterised by a plateau at $\bar{P}_{res} = 1$, obtained for values of ($\bar{\tau}_e$, $\kappa$) in the form of a horseshoe, shown in blue in Fig. 9(e). Interestingly, plotting this surface for other values of the damping ratio $\xi$ leads to a similar shape, with the plateau being translated toward the positive $\kappa$ when $\xi$ increases. It can then be concluded that there exist a value $\kappa_{lim}(\xi)$, shown by a blue dot in Fig. 9(e), that conducts to two situations. If the coupling factor is $\kappa < \kappa_{lim}$, then the optimal value of the electrical time constant is $\bar{\tau}_e = RC\omega_l = 1$ and the maximal power $P_{max}$ cannot be harvested: $P_h < P_{max}$ ($\bar{P}_{res} = \bar{P} < 1$). This situation is the one of Figs. 8(a,c,e), 9(a,b) and the red dot in Fig. 9(e). Otherwise, for any value $\kappa > \kappa_{lim}$, two optimal values of $\bar{\tau}_e$ exist and are associated to the maximal power that can be harvested, $P_h = P_{max}$ ($\bar{P} = 1$). This situation is the one of Figs. 8(b,d,f), 9(c,d) and the orange dots in Fig. 9(e). It means that if $\kappa > \kappa_{lim}(\xi)$, provided a suitable choice of the resistive load $R$ ($\bar{\tau} = RC\omega_l$), the maximal power that can be harvested. The evolution of $\kappa_{lim}(\xi)$ as a function of $\xi$ is shown in Fig. 12(a), showing that the larger the damping ratio is, the larger is the $\kappa$ threshold to harvest the maximum power. A simple analytical expression, $\kappa_{lim} = 2\sqrt{\xi + \xi^2}$, is provided in [86,87]. It exactly matches our numerical results of Fig. 12(a).

For physical insights, Fig. 8 also shows the vibration response of the PEH. Its amplitude, in Fig. 8(a,b), is plotted scaled by $F\Phi(L)/(2\xi)$, for reference to the amplitude at resonance in the SC case. Indeed, setting $\bar{\tau}_e = 0$ and $\bar{\Omega} = 1$ in Eq. (64) gives $|\hat{a}| = F/(2\xi)$. With $\Phi_1(L) = 1.833$ (config. (a) of Table 2), the tip displacement of the PEH is $w_{H1}/L = \Phi_1(L)|\hat{a}|$, such that $w_{H1}/L \times 2\xi/(\Phi_1(L)F) = 1$ at the SC resonance, whatever be the forcing $F$ and the damping $\xi$. Fig. 8(b) shows that the maximal power is harvested with a PEH displacement close to half its value at resonance if $\kappa > \kappa_{lim}$ whereas it is larger in the case $\kappa < \kappa_{lim}$, as shown in Fig. 8(a). The surfaces of Fig. 9(a–d) show another view of those results.
<figure><img src="image_16.png" /><figcaption>**Fig. 9.** 3D plots of the response and performance of a linear DPEH. (a,c) tip amplitude of the PEH normalised by its resonant SC value $F\Phi_1(L)/(2\xi)$, as a function of electrical time constant $\bar{\tau}_e$ and forcing frequency $\bar{\Omega}$, for $\xi = 0.01$ and for $\kappa = 0.16 < \kappa_{\text{lim}}$ (a) and $\kappa = 0.316 > \kappa_{\text{lim}}$ (c); (b,d) same plot with the dimensionless power $\bar{P}$; (e) dimensionless power $\bar{P}_{\text{res}}$ at power resonance as a function of $\bar{\tau}_e$ and $\kappa$ for $\xi = 0.01$. The red/orange curves in (a,b/c,d) are the optimal ones in term of harvested power, shown as red/orange bullet in (e). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Fig. 8(c,d) show that the power resonance (for $\bar{\Omega} = \bar{\Omega}_{\text{res}}$) is obtained at a phase lag of $+\pi/2$, whatever be the value of $\bar{\tau}_e$. This is an interesting result, extending the concept of phase resonance to a damped system, since a finite $\bar{\tau}_e$ means a finite resistance $R$. This is an opening to the use of a phase locked loop technique (see e.g. [64,81]) to directly measure the behaviour of the system at its power resonance.

As a conclusion, in addition to the increase of $P_{\text{max}}$ already discussed in Section 4.1, the above discussion about $\bar{P}$ gives insights on the effect of the three parameters ($\bar{\tau}_e, \kappa, \xi$). For any value of ($\kappa, \xi$), there are optimal values of $\bar{\tau}_e$ that maximise the harvested power. They are (i) $\bar{\tau}_e = 1$ if $\kappa < \kappa_{\text{lim}}$ or (b) two possible values linked to the horseshoe curve of Fig. 9(e) if $\kappa > \kappa_{\text{lim}}$. In practice, one must adjust the equivalent impedance $R = \bar{\tau}_e/(\omega_1 C)$ of the electronic circuit accordingly. Then, decreasing the damping factor $\xi$, in addition to increasing the maximum harvested power $P_{\text{max}}$, also decreases the limit value $\kappa_{\text{lim}}$ above which the maximum power $P_{\text{max}}$ may be harvested. This threshold behaviour also shows that in term of harvested power, it is not necessary to increase as much as possible the piezoelectric coupling factor: it is only necessary to have it larger than $\kappa_{\text{lim}}(\xi)$. This behaviour is different than what is observed for piezoelectric shunt damping applications, whose performances are increasing functions of $\kappa$ [67]. However, increasing $\kappa$ above $\kappa_{\text{lim}}$ can be interesting to increase the frequency bandwidth of the energy harvesting, linked to the width of the horseshoe, if a special synchronised switched electronics is used [88].

## 4.3. A nonlinear DPEH

In this section, the effect of the geometrical nonlinearities on the energy harvesting performance and optimisation of a PEH under direct forcing (DPEH) are studied. The main difference with respect to a linear PEH addressed in Section 4.2 is that all results are here dependent on an additional parameter: the amplitude of the forcing, because of the nonlinearities. The considered model is Eqs. (35a,c) without parametric forcing ($\bar{\delta} = 0$). Because of the reduced validity limit of the analytical models (see Section 3.4 and Fig. 5), all results of this section are obtained thanks to the numerical continuation of periodic solutions, with the software Manlab [77,78], based on a HBM method with $H = 30$ harmonics.

To analyse the effect of nonlinearities, Figs. 10 and 11 are the analogues of Figs. 8 and 9, but now with the geometrical nonlinearities included. Because of the amplitude dependence due to the nonlinearities, the PEH response is shown for the largest amplitude obtained within the validity range of the model, to comfortably observe the effect of those nonlinearities. The validity limit is a vibration amplitude of the tip of the PEH of $w_{H1} = 0.7L$ at the SC resonance (see Fig. 5), which corresponds to a forcing amplitude $F = 0.0076$ with a damping ratio $\xi = 0.01$.
<figure><img src="image_17.png" /><figcaption>**Fig. 10.** Frequency response of a nonlinear DPEH several resistive loads $R$ ($\bar{\tau}_e$ takes 20 values from 0.001 (≈ SC, blue curve) to 1000 (≈ OC, yellow curve)) and constant damping $\xi = 0.01$ and coupling factor $\kappa = 0.316$. The forcing $F = 0.0076$ has been chosen to obtain a maximal tip amplitude of the harvester of $w_{H1} = 0.7L$, close to the validity limit of the model. (a) tip amplitude (first harmonic scaled by the PEH length $w_{H1}/L$); (b) phase lag $\alpha$ between the displacement and the external acceleration; (c) dimensionless power $\bar{P}$. The × markers show the power resonance and the red curves are the optimal ones (with maximal power at resonance). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 11.** 3D plots of the response and performance of a nonlinear DPEH. (a) tip amplitude of the PEH (first harmonic scaled by the PEH length $w_{H1}/L$) as a function of electrical time constant $\bar{\tau}_e$ and forcing frequency $\bar{\Omega}$, for $\xi = 0.01$, $\kappa = 0.316$ and $F = 0.0076$; (b) same plot with the dimensionless power $\bar{P}$. The red curves are the optimal ones in term of harvested power.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 12. For linear and nonlinear DPEHs: (a) effect of the normalised direct forcing $F/\xi$ on the law $\kappa_{lim} = f(\xi)$; (b) tip amplitude of vibration of the PEH (first harmonic scaled by the PEH length $w_{H1}/L$) at the short circuit (SC) resonance and at the optimal power resonance (any location on the power plateau of Fig. 9(e)), as a function of normalised direct forcing $F/\xi$ and for several damping ratios $\xi$; (c) normalised power $\bar{P}_{pla}$ on the “horseshoe” plateau of Fig. 9(c) as a function of normalised direct forcing $F/\xi$ and for several damping ratios $\xi$.</figcaption></figure>

The main effect of nonlinearities, visible in Figs. 10(a-c) and 11(a), is the hardening behaviour of the vibration mode, which bends the resonance curves to the high frequencies. Because the first mode of the cantilever beam is considered, this effect is very slight (the resonance frequency change is 1% for $w_{H1} = 0.7L$).

Apart from this hardening effect, the main characteristics of a linear PEH, explored in Section 4.2, are still valid in the present nonlinear case. In particular, the harvested power is written in the same way as in the linear case (Eq. (65)) by:

$$P_h = P_{max} \cdot \bar{P}(\xi, \kappa, \bar{\tau}_e, \bar{\Omega}, F), \quad (66)$$

in which the dimensionless power $\bar{P}$ is now also a function of the forcing amplitude $F$. There is no closed form expression for $\bar{P}$, but it is computed numerically (directly as an auxiliary variable in Manlab, see Appendix D). Figs. 10(c) and 11(b) show $\bar{P}$ as a function of the parameters, and its maximum value is very close to 1, leading to think that the dependence of $\bar{P}$ (and thus of $P_h$) on $F$ is slight. Moreover, plotting the power at resonance $\bar{P}_{res}$ as a function of $(\kappa, \bar{\tau}_e)$ in the present nonlinear case gives the same horseshoe shaped surface plot than the one in Fig. 9(e) (this plot is not shown here since it is visually identical to Fig. 9(e)), leading to the same behaviour: a plateau of harvested power, provided the piezoelectric coupling factor $\kappa$ is larger than a limit value $\kappa_{lim}$ and that $\bar{\tau}_e$ is adjusted at one of its two optimal values.

Fig. 12 quantitatively illustrates the effect of $F$ on the system's behaviour. First, Fig. 12(a) shows the limit value $\kappa_{lim}(\xi, F)$ (it is now a function of $F$ in addition to $\xi$) as a function of the damping ratio for several values of $F/\xi$. Then, Fig. 12(b) shows the cantilever tip vibration amplitude as a function of $F/\xi$ for several values of $\xi$. The two cases of the SC resonance (the resonance on the left of Fig. 10(a)) and the optimal power resonance (the left red 'o' marker of Fig. 10(a)) are considered, the latter being half smaller than the former. Finally, Fig. 12(c) shows the dimensionless harvested power $\bar{P}$ as a function of $F/\xi$ for several values of $\xi$. For all those cases, it is observed that the nonlinear curves are almost merged with the curves from the linear PEH, under a large vibration amplitude range of $w_{H1} \le 0.7L$, such that the dependence of the system's behaviour on the forcing amplitude $F$, in term of response amplitude, harvested power and $\kappa_{lim}(\xi)$, can be neglected in practice: the geometrical nonlinearities do not affect significantly the performance of the PEH. Only the resonance frequency is slightly changed.

Interestingly, Fig. 10(b) shows that the power resonance still appears at a phase lag of $+\pi/2$, despite the hardening behaviour, a result in accordance with [81].

### 4.4. A (nonlinear) PPEH

In this section, the energy harvesting performance and the optimisation of a PEH under parametric forcing (PPEH) are studied. The main difference with respect to a DPEH (addressed in the previous sections Sections 4.2 and 4.3) is linked to the parametric response of the resonator. In our case, the considered model is Eqs. (35a,c) without direct forcing ($F = 0$). As already addressed in Section 3.2, under harmonic parametric forcing $\delta(t) = \delta \cos(2\Omega\bar{t})$ at a frequency close to twice the natural frequency of the PEH ($\Omega \approx \omega_l \Leftrightarrow \bar{\Omega} \approx 1$), a resonant response at frequency $\Omega$ is observed above a forcing threshold $\delta_{cr}$ (i.e. for $\delta > \delta_{cr}$). Unlike in the direct cases, the vibration amplitude at resonance in the present parametric case is not a linear function of the forcing amplitude $\delta$ as it grows rapidly after the threshold $\delta_{cr}$, an a priori interesting property to harvest energy (an illustration is given in Fig. 16(b), that will be analysed hereafter).

Figs. 13(a) and 14(a) illustrate the frequency response of the PPEH for various resistive loads $R$ and for a constant forcing amplitude $\delta = 1.048 \delta_{cr}^{sc} (\delta_{cr}^{sc} = 4\xi$ is the threshold in SC), chosen to obtain, at the SC resonance, the largest amplitude within the validity range of the model (a vibration amplitude of the tip of the PEH of $w_{H1} = 0.7L$, see Fig. 5). The behaviour resembles that of
<figure><img src="image_20.png" /><figcaption>Fig. 13. Frequency response of a PPEH for several resistive loads R ($\bar{\tau}_e$ takes 20 values from 0.01 (≈ SC, blue curve) to 100 (≈ OC, yellow curve)) and constant damping $\xi = 0.01$ and coupling factor $\kappa = 0.158$. The forcing $\delta/(4\xi) = 1.048$ has been chosen to obtain a maximal tip amplitude of the harvester of $w_{H1} = 0.7L$, close to the validity limit of the model. (a) tip amplitude (first harmonic scaled by the PEH length $w_{H1}/L$) of the PEH; (b) phase lag $\alpha$ between the displacement and the external acceleration; (c) dimensionless power $\bar{P}$. The × markers show the power resonance and the red curves are the optimal ones (with maximal power at resonance). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

a DPEH as two large resonances are obtained in the two extreme cases of SC ($R = 0$) and OC ($R \to +\infty$). In between, the dissipative effect of the resistive load decreases the resonant response, but, unlike in the DPEH case, it also increases the threshold $\delta_{cr}$, which makes the parametric response disappearing in a given resistance range, for which $\delta < \delta_{cr}$. This question, illustrated in Fig. 15, will be addressed hereafter. That is why the coupling factor was chosen $\kappa = 0.158$ in this section (half the value of config(a), Table 2), to obtain a noticeable OC resonance. Indeed, $\kappa = 0.316$ leads to a higher threshold value such that OC resonance disappears. One may notice that the OC parametric resonance is smaller than the SC one. This is explained because its natural frequency is slightly larger, thus increasing the threshold ($\delta_{cr}^{oc} = 4\xi\sqrt{1+\kappa^2} > \delta_{cr}^{sc} = 4\xi$, see Eq. (47)).

To consider quantitatively the harvested power $P_h$, we normalise it in the same way as for the nonlinear DPEH (Eq. (66)), by the maximal value $P_{max}$ (Eq. (58)) with $\gamma$ now the parametric input acceleration, as:

$$P_h = P_{max} \cdot \bar{P}(\xi, \kappa, \bar{\tau}_e, \bar{\Omega}, \delta), \qquad (67)$$

where $\bar{P}$ is now a function of $\delta$ (in place of $F$). Figs. 13(b) and 14(b) show that, analogously to what is observed for a DPEH, there exist optimal values of $R \leftrightarrow \bar{\tau}_e$ that maximise $\bar{P}$. However, in the present parametric case, since the OC parametric resonance is lower than the SC one (see above), the harvested power is smaller for the OC resonance than for the SC resonance. Consequently, in the following, only the optimal harvested power close to the SC resonance will be considered.

Following the same process than for a DPEH, Figs. 14(c,d) (analogous to Fig. 11(e)) shows the dimensionless power $\bar{P}_{res}(\xi, \kappa, \bar{\tau}_e, \delta)$ at resonance as a function of ($\bar{\tau}_e, \kappa$), for $\xi = 0.01$ and $\delta = 1.048 \delta_{cr}^{sc}$. Interestingly, the obtained surface has still a “horseshoe” shape, with its right part truncated (for large values of $\bar{\tau}_e$, on the side of the OC resonance), because, as already observed, $\delta_{cr}^{oc} > \delta_{cr}^{sc}$. There is
<figure><img src="image_21.png" /><figcaption>**Fig. 14.** Various plots of the response and performance of a PPEH. (a) tip amplitude (first harmonic scaled by the PEH length $w_{H1}/L$) of the PEH as a function of electrical time constant $\bar{\tau}_e$ and forcing frequency $\bar{\Omega}$, for $\xi = 0.01$, $\kappa = 0.158$ and $\delta/(4\xi) = 1.048$; (b) same plot with the dimensionless power $\bar{P}$; (c,d) dimensionless power $\bar{P}_{res}$ at resonance as a function of $\bar{\tau}_e$ and $\kappa$ for $\xi = 0.01$. The red curve in (a,b) is the optimal one in term of harvested power, shown by a red bullet in (c,d). The blue curve in (c,d) corresponds to the power plateau ($\bar{P}_{res} = \bar{P}_{pla}$). The red curve in (d) corresponds to the parametric threshold above which a non-zero response is obtained. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 15.** Excitation threshold $\delta_{cr}$ normalised by $\delta_{cr}^{sc} = 4\xi$, as a function of $\bar{\tau}_e$ and $\kappa$, for $\xi = 0.01$. The blue curve corresponds to the threshold $\delta_{cr}^{pla}$ for the half “horseshoe” curve on the power plateau of Fig. 14(c). The red curve is the level set $\delta/(4\xi) = 1.048$, corresponding to a maximal tip amplitude of $w_{H1} = 0.7L$. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

also a plateau, only on the left part of the horseshoe (on the side of the SC resonance), which is related to a constant¹ (and optimal) value $\bar{P}_{pla}$ of $\bar{P}$. The corresponding “half horseshoe” curve of optimal values of $(\kappa, \bar{\tau}_e)$ is shown in blue in Figs. 14(c,d).

<footnote n="4">⁴ A slight decrease of around 3% is observed between $\bar{P}_{pla}(\kappa = 0.5)$ and $\bar{P}_{pla}(\kappa = \kappa_{lim})$ for $\xi = 0.01$. Increasing $\xi$ tends to increase the gap between $\bar{P}_{pla}(\kappa = 0.5)$ and $\bar{P}(\kappa = \kappa_{lim})$ ($\approx 11\%$ for $\xi = 0.1$).</footnote>
<figure><img src="image_23.png" /><figcaption>**Fig. 16.** For a PPEH: (a) effect of the normalised parametric forcing $\delta/(4\xi)$ on the law $\kappa_{lim} = f(\xi)$; (b) tip amplitude of vibration of the PEH (first harmonic scaled by the PEH length $w_{H1}/L$) at the short circuit (SC) resonance and at the optimal power resonance (any location on the power plateau of Fig. 14(c)), as a function of normalised parametric forcing $\delta/(4\xi)$ and for several damping ratios $\xi$; (c) normalised power $\bar{P}_{pla}$ on the “horseshoe” plateau of Fig. 9(e) as a function of normalised parametric forcing $\delta/(4\xi)$ and for several damping ratios $\xi$.</figcaption></figure>

Unlike for a DPEH, the parametric response disappears when the forcing is beyond the threshold, $\delta < \delta_{cr}$. As shown by Eq. (46), $\delta_{cr}$ is a function of $(\xi, \kappa, \bar{\tau}_e)$. Plotting it as a function of $(\kappa, \bar{\tau}_e)$ in Fig. 15 shows that it is an increasing function of $\kappa$ (which is logical since the damping effect of the electrical resistance is applied to the mechanical system through the piezoelectric coupling) and that it is the larger for a matched electrical impedance, around $\bar{\tau}_e = 1$. Fig. 15 also shows the level set $\delta = 1.048 \times 4\xi$ (in red): for all values of $(\kappa, \bar{\tau}_e)$ for which the surface is above this level, no parametric response is obtained for this value of $\delta$. This red curve is also superimposed in Fig. 14(d), explaining quantitatively why the harvested power is zero in the middle area of Figs. 14(c,d).

Fig. 15 also shows the $(\bar{\tau}_e, \kappa)$ line (in blue) that corresponds to the “half horseshoe” power plateau curve of Figs. 14(c,d). It was verified numerically that this blue curve is indistinguishable from one of the level set of the surface, such that $\delta_{cr}$ on this line is a constant, denoted as $\delta_{cr}^{pla}$. Since the location of the power plateau is a function of $(\xi, \delta)$ (refer to the upcoming discussion on $\kappa_{lim} = f(\xi, \delta)$), $\delta_{cr}^{pla}$ is also a function of $(\xi, \delta)$. However, even if the threshold surface of Fig. 15 is a function of $\xi$ (see Eq. (46)), since the “half horseshoe” curve always lies in the lower part of the threshold surface, it was numerically observed that $\delta_{cr}^{pla}/4\xi$ apparently does not depend on $\xi$ (it thus depends only on $\delta$). We obtained $\delta_{cr}^{pla}(1.048\delta_{cr}^{sc}) = 1.021\delta_{cr}^{sc}$, whereas $\delta_{cr}^{pla}(1.093\delta_{cr}^{sc}) = 1.038\delta_{cr}^{sc}$.

Interestingly, Fig. 13(b) shows that the power resonance (when $\bar{P} = \bar{P}_{res}$) is obtained at a given phase lag of $\alpha = \pi/4$ or $\alpha = -3\pi/4$ (as for the conservative backbone curves, see Section 3.2), whatever be the values of the parameters. As for the DPEH (for which $\alpha = \pi/2$), this result extends the concept of phase resonance to a damped parametric system and paves the way to using a phase locked loop controller to directly measure the behaviour of the system at its power resonance.

The “half horseshoe” surface of Figs. 14(c,d) depends on $\xi$ and $\delta$ and shows, as for a DPEH, a limit value $\kappa_{lim}$ of $\kappa$, above which the maximal power at resonance (denoted $\bar{P}_{pla}$) can be harvested. Obviously, $\kappa_{lim}$ is a function of $\xi$ and $\delta$ and this dependence was numerically evaluated. The result is shown in Fig. 16(a). As in the DPEH case, $\kappa_{lim}$ is an increasing function of $\xi$, such that the lower $\xi$ is, the better it is for this $\kappa$ limit. However, unlike the nonlinear DPEH case, $\kappa_{lim}$ depends on the forcing ($\delta$ in our case, whereas it was $F$ in the DPEH case, see Fig. 12(a)): $\kappa_{lim}$ is an increasing function of $\delta$, such that increasing the parametric forcing also increases the $\kappa$ limit, thus going in the wrong way. On the contrary, Fig. 16(a) shows that $\kappa_{lim}$ for the PPEH is almost one order of magnitude smaller than $\kappa_{lim}$ for the DPEH, leading to think that poorly coupled piezoelectric elements ($\kappa < 0.04$ for $\xi < 0.01$) may be sufficient to obtain the best power harvesting under parametric forcing.

Finally, Figs. 16(b,c) show the effect of the parametric forcing amplitude $\delta$ and the damping ratio $\xi$ on (i) the PEH vibration amplitude and (ii) the maximum harvested power at the plateau of Fig. 14(c). The first result is that, due to the intrinsic nonlinear response of a parametric resonator, both of those quantities nonlinearly depend on $\delta$. As the vibration amplitude of the PEH grows rapidly above the threshold $\delta_{cr}$, it is worth noticing that the interval of excitation amplitude in which the vibration amplitude remains below the large amplitude limit $w_{H1} = 0.7L$ is relatively tight ($1 < \delta/(4\xi) < 1.048$ for the SC resonance). As for the maximal power, unlike for a DPEH, $\bar{P}_{pla}$ is a growing function of $\delta$, close to a parabolic shape. In addition, it is observed that the effect of $\xi$ on those two quantities can be neglected (as for the nonlinear DPEH).

However, the main result, visible in Fig. 14(c), is that the value of the normalised power $\bar{P}_{pla}$ is much lower than its value for a DPEH, almost two orders of magnitude lower (of the order of 0.013 instead of 1 for large amplitude of vibration of the PEH, of the order of $w_{H1} = 0.7L$). This result will be analysed and explained in the next section.

# 5. Direct and parametric energy harvesting mismatch: analogy with a pendulum system and power transfer

## 5.1. A pendulum

In order to offer more general results illustrating the performance differences between the direct and parametric responses of the cantilever PEH of the previous sections, we consider the multidirectional excited pendulum shown in Fig. 17(a). Pendulum systems are considered elsewhere for energy harvesting in both parametric and direct driving (see e.g [89,90]), but it is here selected as the
simplest rigid body – and qualitative – model of a cantilever beam. The pendulum is of length l, with all of its inertia concentrated at its tip. To model the geometrical nonlinearities of a cantilever beam, the pendulum is connected to the ground with a nonlinear torsional spring, that creates a restoring moment with linear and cubic parts, with no gravity considered.

For this system, the equation of motion governing the position angle $\theta(t)$, with additional linear damping, reads:

$$ \ddot{\theta} + 2\xi\omega_0\dot{\theta} + \omega_0^2\dot{\theta} + \Gamma\theta^3 - \sin\theta\frac{\gamma_x}{l} = -\cos\theta\frac{\gamma_y}{l}, \quad (68) $$

where $(\gamma_x(t), \gamma_y(t))$ are the two components of the applied ground acceleration, $\omega_0$ is the natural frequency and $\Gamma$ the cubic nonlinearity coefficients related to the torsional spring.

We render this equation dimensionless with the following variables and harmonic forcing is considered:

$$ \bar{t} = \omega_0 t, \quad \bar{\Gamma} = \frac{\Gamma}{\omega_0^2}, \quad \gamma_x = l\omega_0^2\bar{\gamma}_x \cos 2\bar{\Omega}\bar{t}, \quad \gamma_y = l\omega_0^2\bar{\gamma}_y \cos \bar{\Omega}\bar{t}, \quad (69) $$

where $\bar{\Omega}$ is the dimensionless excitation frequency. As a simplification, and also to use the analytical results of Section 3, we linearise the sine and cosine terms, such that Eq. (68) becomes:

$$ \ddot{\theta} + 2\xi\dot{\theta} + (1 - \bar{\gamma}_x \cos 2\bar{\Omega}\bar{t})\theta + \bar{\Gamma}\theta^3 = -\bar{\gamma}_y \cos \bar{\Omega}\bar{t}. \quad (70) $$

The oscillations of the pendulum, neglecting higher harmonics, are $\theta(t) \approx \theta_m \cos(\bar{\Omega}\bar{t} + \alpha)$. The parametric threshold is introduced as $\bar{\gamma}_{cr} = 4\xi$ (see Eq. (47), such that parametric oscillations are observed only if $\bar{\gamma}_y > \bar{\gamma}_{cr}$). We are interested in the amplitude at resonance of the pendulum.

Firstly, the pendulum behaviour is linearised and considered under direct excitation ($\bar{\gamma}_x = \bar{\Gamma} = 0$). In this situation, the amplitude of oscillations at resonance is:

$$ \theta_{\text{lin}} \approx \frac{\bar{\gamma}_y}{2\xi} = 2\frac{\bar{\gamma}_y}{\bar{\gamma}_{cr}}. \quad (71) $$

Secondly, nonlinear response is considered. As proved in Section 3 and in [79], since here the inertial nonlinearities are zero ($\bar{\Gamma} = 0$), the HBM-H1 solutions of Eqs. (39) and (44) are acceptable. For direct excitation ($\bar{\gamma}_x = 0$), the resonance is obtained for $\alpha = \pi/2$, such that Eqs. (38) and (48) with $\kappa = \bar{\kappa} = 0$ give:

$$ \bar{\gamma}_y^2 = T_2^2 \theta_m^2 = 4\xi^2 \bar{\Omega}^2 \theta_m^2, \quad \bar{\Omega}^2 = 1 + \frac{3\bar{\Gamma}\theta_m^2}{4}, \quad (72) $$

$$ \Rightarrow \quad \theta_D = \sqrt{\frac{2}{3\Gamma}\left(\sqrt{1 + \frac{3\Gamma\bar{\gamma}_y^2}{4\xi^2} - 1}\right)} = \sqrt{\frac{2}{3\Gamma}\left(\sqrt{1 + 12\Gamma\left(\frac{\gamma_y}{\gamma_{cr}}\right)^2} - 1\right)} \approx \theta_{\text{lin}}. \quad (73) $$

Finally, under parametric excitation ($\bar{\gamma}_y = 0$), the phase lag at resonance is $\alpha \in \{\pi/4, -3\pi/4\}$ and one obtains with Eqs. (43) and (48):

$$ \bar{\gamma}_x^2 = 4T_2^2 = 16\xi^2 \bar{\Omega}^2 \Rightarrow \theta_p = \frac{2}{\sqrt{3\Gamma}} \sqrt{\left(\frac{\bar{\gamma}_x}{\bar{\gamma}_{cr}}\right)^2 - 1}. \quad (74) $$

One may observe that the amplitude at resonance, in all the three above cases, are sole function of $\bar{\gamma}/\bar{\gamma}_{cr}$ and $\bar{\Gamma}$, where $\bar{\gamma}$ is either $\bar{\gamma}_x$ (for parametric excitation) or $\bar{\gamma}_y$ (for direct excitation), producing a convenient scaling. Fig. 17(b) shows $\theta_{\text{lin}}, \theta_D, \theta_p = f(\bar{\gamma}/\bar{\gamma}_{cr})$ and will be analysed hereafter.

## 5.2. Comparison to a PEH

The pendulum angle $\theta(t)$ is compared to the rotation angle of the tip cross-section of the PEH. A single layer elastic beam (without piezoelectric patches) as well as a configuration (a) PEH are considered (see Table 2). For this, a HBM-H30 Manlab solving of Eq. (35a) for several amplitudes of $F$ in direct resonance (with $\bar{\delta} = 0$) and $\bar{\delta}$ in parametric resonance (with $\bar{\Gamma} = 0$). The first harmonic amplitude $a_{H1}$ of $\bar{q}(t)$ is identified, at resonance for the elastic beam, and at the optimal power resonance (on the plateau of Figs. 9(e) and 14(c)) for the PEH. In each case, the first harmonics of the beam's tip rotation angle is estimated as:

$$ \theta(1, \bar{t}) \approx \bar{w}'(1, \bar{t}) = \Phi'_1(1)\bar{q}(\bar{t}) \Rightarrow \theta_{H1} \approx \Phi'_1(1)a_{H1}. \quad (75) $$

Then, $\theta_{H1}$ is plotted in Fig. 17(b) (and denoted as the truncated model) as a function of:

$$ \frac{\gamma_y}{\gamma_{cr}} = \frac{\bar{\Lambda}F}{4\xi^2S}, \quad \text{or} \quad \frac{\gamma_x}{\gamma_{cr}} = \frac{\bar{\delta}}{4\xi^2}. \quad (76) $$

respectively for direct and parametric excitations, where $S, \bar{\Lambda}$ and $\Phi'_1(1)$ are read in Table 2. The direct and parametric response of a homogeneous cantilever beam, computed with the geometrically exact beam model of [40,41] (FEHBM)), are also shown in Fig. 17(b).

First of all, considering Figs. 12(b) and 16(b), the vibration amplitude of the PEH at the optimal power resonance, as a function of $\gamma/\gamma_{cr}$, is almost independent of the damping ratio $\xi$. As a consequence, the curves of Fig. 17(b) for the elastic cantilever and the PEH are master curves, that depend only on the mechanical structure of the PEH (linear and nonlinear stiffness).
<figure><img src="image_24.png" /><figcaption>**Fig. 17.** Pendulum under direct and parametric excitation with a torsional spring stiffness. (a) Scheme of the system. (b) Amplitude of the first harmonic of the oscillations $\theta(t)$ of the pendulum as a function of the scaled input acceleration $\bar{\gamma}/\bar{\gamma}_{cr}$, for several values of $\Gamma$ in the direct (solid curves) and parametric (dashed curves) cases (from blue to yellow). The additional lines with markers are the tip rotation at resonance of a cantilever PEH in direct (triangle) and parametric (circle) excitation: (purple) config (a) of Table 2 PEH, 3rd. order truncated model HBM-H30 for ($\kappa, \tau$) on the plateau; (red) elastic single layer cantilever, 3rd. order truncated model HBM-H30; (black) elastic single layer cantilever, geometrically exact model FEHBM. The maximum amplitude ($\approx 0.3\pi$ for the single layer beam and $\approx 0.5\pi$ for the PEH, different because of the mode shapes) of the truncated HBM-H30 model (purple, red) correspond to its validity range ($w_{H1} = 0.7L$), see Fig. 5. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. 18.** Energy harvesting response of the config (a) PEH (see Table 2) in direct (yellow) and parametric (blue) excitation at optimal power harvesting for $\bar{\tau}_e = 1$; (a) amplitude response; (b,c) dimensionless powers (direct and parametric input powers $\bar{P}_e, \bar{P}_p$, dissipated power $\bar{P}_d$ and harvested power $\bar{P}_h$) as a function of the input acceleration $\gamma/\gamma_{cr}$ and the tip vibration amplitude of the PEH $w_{H1}/L$. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Interestingly, Fig. 17(b) shows that the direct and parametric responses of the pendulum are, in this $(\theta, \bar{\gamma}/\bar{\gamma}_{cr})$ plot, qualitatively very similar. The main result is that the parametric pendulum can only overcome the amplitude of the direct pendulum, for the same input acceleration, if the amplitude of vibration is very large, of the order of $5\pi/8$ for the chosen cubic coefficients. An analogue result is observed for a cantilever beam, for which the FEHBM model predicts a match between direct and parametric response for a tip angle of $\approx 3\pi/4$.

Consequently, Fig. 17(b) gives a first insight into the discrepancy in harvested power between a DPEH and a PPEH observed in Section 4.4: to obtain a relatively large amplitude of vibration of the PEH ($w_{H1} = 0.7L \Leftrightarrow \theta_{H1} \approx 0.5\pi$, illustrated in Fig. 5(b)), the necessary input acceleration is 1.79 times higher in parametric excitation than in direct one. Since the harvested power is proportional to the square of the input acceleration, this 1.79 ratio in acceleration would give a 3.2 ratio in power, in favour of the DPEH with respect to the DPEH. However, this does not explain the large discrepancy in power observed by comparing the horseshoe amplitudes of Figs. 12(c) and 16(c), for which the DPEH is superior in power by a factor $1/0.013 \approx 77$.

To understand this discrepancy, Figs. 18(b,c) show the evolution of the several powers (direct $P_e$ and parametric $P_p$ input powers, dissipated $P_d$ and harvested $P_h$) as a function of the input acceleration $\gamma/\gamma_{cr}$ and the tip vibration amplitude of the PEH $w_{H1}/L$. Instead of scaling them by $P_{max}$, which includes $\gamma^2$, the considered powers are here made dimensionless by Eqs. (34) ($\bar{P}_\eta = P_\eta/(M_i L^2 \omega_i^3)$, $\eta = h, d, e, p$), to observe there increase with the input acceleration. One observes in Fig. 18(b,c) that in direct excitation, the harvested and dissipated powers are exactly half of the input power: $P_h = P_d = P_e/2$, as predicted in Section 4.1.1. Notice that the nonlinearities do not change significantly the result (see Fig. 12(c)). On the contrary, for parametric excitation, one first observes that the input power $P_p$ for the PPEH is almost half the one ($P_e$) for the DPEH, at equivalent vibration amplitude. In addition, the dissipated power is almost the same for DPEH and PPEH, at equivalent vibration amplitude, since it is proportional
to the squared vibration amplitude ($\tilde{P}_d = 2\xi\dot{q}^2$). A slight difference is noticed, due to the power resonance frequency, which is slightly lower for the PPEH than for the DPEH (see Fig. 18(a)), such that $\tilde{P}_d$ for the PPEH is slightly smaller. The power resonance frequencies difference is explained by the fact the resonance depends on the values of $\kappa$ and $\tau_e$, which are different for the DPEH and PPEH plateaus. At the end, since the harvested power is the difference of input and dissipated powers ($P_h = P_p - P_d$), it is very small in parametric excitation. This explains the large ratio (77) in terms of harvested power between the PPEH and the DPEH.

## 6. Conclusions

In this work, a complete study on the modelling, behaviour and performance of a direct (DPEH) and a parametric (PPEH) cantilever PEH in large rotation is proposed, with three groups of results.

Firstly, the inextensible model of [42], for a cantilever beam in large rotation, is extended to the case of (i) a laminated beam including (ii) piezoelectric layers and (iii) their coupling with an electric circuit. It is discretised with the help of an eigenmode expansion, giving rise to a single degree of freedom reduced order model (ROM). Its validity range is investigated with comparison to a full finite-element numerical method, leading to a maximal tip amplitude displacement of ≈ 0.7 times the length of the beam.

Secondly, the ROM is solved by classical perturbation methods (the method of multiple scale (MMS) and the harmonic balance method, reduced to a single harmonic (HBM-H1)). The validity range of those two methods is investigated and found to be too small for the large-amplitude application required for energy harvesting. This is due to the presence of the inertial nonlinear terms as well as the parametric excitation. A numerical method (a continuation method based on a HBM with $H = 30$ harmonics) is then proposed and recommended.

Thirdly, the optimisation of a linear DPEH (as a reference), compared to the one of a nonlinear DPEH and a PPEH, is investigated. The free parameters are the piezoelectric coupling factor $\kappa$, the impedance of the electrical circuit load, the mechanical damping ratio and, in addition, the input acceleration in case of nonlinearities. It is found that the effect of the geometrical nonlinearities on the DPEH can be neglected. On the contrary, the PPEH is intrinsically nonlinear, with a threshold effect and a nonlinear evolution of the vibration amplitude as a function of the excitation. However, a behaviour qualitatively similar to that of the DPEH is found in terms of maximal power with a horseshoe shape. In particular, a limit value of $\kappa$ is obtained, above which an optimum of harvested power is achievable, provided the electric circuit impedance is adjusted.

Finally, the main (and very disappointing) result is that the harvested power under parametric excitation is found almost 1.9 order of magnitude lower than under direct excitation, for equivalent input accelerations. Some physical explanations are raised, by comparison to a multibody pendulum equivalent system. The main reason is that the input power under parametric excitation is almost half the one in direct excitation, at equivalent vibration amplitude of the resonator, whereas the dissipated power in the mechanical damper is the same. An opening, left for future work, would be to try to maximise the ratio of the harvested power and the dissipated power $P_h/P_d$ for the PPEH, that may lead to exhibit better performances in other areas of the parameter space.

In contrast, it turns out that the limit value of $\kappa$ is almost one order of magnitude smaller in the parametric case in comparison to the direct case, leading to think that poorly coupled piezoelectric elements may be sufficient to obtain the best power harvesting under parametric forcing. In addition, testing a multibody harvester in the shape of a pendulum, with electromagnetic transduction, may be a promising way for obtaining competitive parametric harvesters.

### CRediT authorship contribution statement

**Hugo Fayolle:** Writing – review & editing, Writing – original draft, Validation, Investigation, Conceptualization.
**Christophe Giraud-Audine:** Writing – review & editing, Validation, Supervision, Methodology, Investigation.
**Olivier Thomas:** Writing – review & editing, Writing – original draft, Validation, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

The French région Hauts-de-France is thanked for participating to the PhD grant of the first author. This article is dedicated to the memory of Simon Benacchio who left us far too early. So long, Simon !
## Appendix A. Coefficients of the modal expansion

The coefficients with dimension of the single modal expansion (30) are given by:

$$ \Gamma = \frac{2}{M_i} \int_{0}^{L} \hat{D}(\Phi_i' \Phi_i'')^2 dx, \qquad (A.1a) $$

$$ \Pi = \frac{1}{M_i} \int_{0}^{L} \Phi_{i}^{\prime 2} \int_{L}^{x} m \int_{x}^{0} \Phi_{i}^{\prime 2} dx dx dx, \qquad (A.1b) $$

$$ \Lambda = \frac{1}{M_i} \int_{0}^{L} \Phi_{i}^{\prime 2} \int_{L}^{x} m dx dx, \qquad (A.1c) $$

$$ \tilde{\Lambda} = \frac{1}{M_i} \int_{0}^{L} \Phi_{i}^{\prime 4} \int_{L}^{x} m dx dx, \qquad (A.1d) $$

$$ S = \frac{1}{M_i} \int_{0}^{L} m \Phi_{i} dx, \qquad (A.1e) $$

$$ \chi = \hat{\Theta}[\Phi_{i}^{\prime 1}]_{x_{-}}^{x_{+}} \qquad (A.1f) $$

$$ \tilde{\chi} = \hat{\Theta}[\Phi_{i}^{\prime 3}]_{x_{-}}^{x_{+}}. \qquad (A.1g) $$

Notice that some integration by parts have been performed, taking into account the clamped-free boundary conditions ($\Phi_i(0) = \Phi_i'(0) = \Phi_i''(L) = \Phi_i'''(L) = 0$).

## Appendix B. Multi-mode modal expansion

In Section 2.3.1, the PDEs (24) are reduced to a single mode model, in order to write a reduced order model adapted to the energy harvesting optimisations of Section 4. In the process, they are made dimensionless after the modal expansion. Here, we present a multi-mode expansion of the PDEs, for which it is convenient to make them dimensionless before the modal expansion.

To this end, we define the following dimensionless variables:

$$ \bar{i} = \omega_0 t, \quad \bar{x} = \frac{x}{L}, \quad \{\bar{w}, \bar{u}\} = \frac{\{w, u\}}{L}, \quad \bar{m} = \frac{m}{m_b}, \quad \bar{D} = \frac{\hat{D}}{D_b}, \quad \bar{\Delta} = L\Delta, \qquad (B.1a) $$

$$ \bar{Q} = \frac{Q}{Q_0} = \frac{Q}{\omega_0 \sqrt{L^3 m_b C}}, \quad \bar{V} = \frac{C}{Q_0} V, \qquad (B.1b) $$

such that the EDPs (25) become:

$$ \left\{ \begin{array}{l} \bar{m}\ddot{w} + [\bar{D}\bar{w}^{\prime\prime}]'' + (\bar{w}^{\prime} [\bar{D}\bar{w}^{\prime}\bar{w}^{\prime\prime}]')' + \left[ \frac{\bar{w}'}{2} \int_x^1 \bar{m}(\xi) \left( \int_\xi^0 \bar{w}^{\prime 2} d\eta \right)^{\cdot\cdot} d\xi \right]' \\ \qquad + \alpha \left[ \bar{\Delta} \left( 1 + \frac{\bar{w}^{\prime 2}}{2} \right) \right]' \bar{V} + \left[ \left( \bar{w}^{\prime} + \frac{\bar{w}^{\prime 3}}{2} \right) \int_x^1 \bar{m}(\xi) d\xi \right]' \ddot{u}_0 = -\bar{m}\ddot{w}_0, \end{array} \right. \qquad (B.2a) $$

$$ \bar{Q} - \bar{V} - \alpha \left[ \bar{w}^{\prime} + \frac{\bar{w}^{\prime 3}}{6} \right]_{x_{-}}^{x_{+}} = 0, \qquad (B.2b) $$

where $\circ'$ and $\odot$ denote now derivatives with respect to $\bar{x}$ and $\bar{i}$ and with

$$ \omega_0 = \frac{1}{L^2} \sqrt{\frac{D_b}{m_b}}, \quad \alpha = \frac{\hat{\Theta} \sqrt{L}}{\sqrt{C D_b}}. \qquad (B.3) $$

In the above equations, the structural constants ($m, \hat{D}$) of the beam have been scaled by nominal values ($m_b, D_b$) that are arbitrary and that have to be chosen to simplify the problem. In beam configurations with a uniform lamination, such as those of Figs. 3(a,b), one may choose ($m_b, D_b$) as the mass per unit length and stiffness of any cross section, such that $\bar{m} = \bar{D} = 1, \forall \bar{x} \in [0, 1]$. In cases with a non-uniform lamination, such as those of Figs. 3(c,d), one may choose the constants of parts I and III for ($m_b, D_b$).

The dynamic electromechanical problem (B.2a,b) is discretised with the help of the following modal expansion on K short-circuit eigenmodes:

$$ \bar{w}(\bar{x}, t) = \sum_{i=1}^{K} \Phi_{i}(\bar{x}) q_{i}(t), \qquad (B.4) $$

with {$q_i$}$_{i \in N^+}$ the modal coordinates. The modes are solutions of the dimensionless version of Eq. (28):

$$ [\bar{D}(\bar{x})\Phi_{i}(\bar{x})]'' - \bar{\omega}_{i}^{2}\bar{m}(\bar{x})\Phi_{i}(\bar{x}) = 0, \qquad (B.5) $$

with $\bar{\omega}_i = \omega_i / \omega_0$ and $\Phi_i(\bar{x})$ the same mode shapes than in (27), with $x$ simply replaced by $\bar{x}$. With the following dimensionless modal mass:

$$ \int_{0}^{1} \bar{m}(\bar{x}) \Phi_{i}(\bar{x})^{2} d\bar{x} = \bar{m}_{i}, \forall i = 1, \dots, K, \qquad (B.6) $$
introducing Eq. (B.4) into Eq. (B.2a), multiplying the result by $\Phi_l(x)$ and subsequently using the orthogonality properties of the modes, one obtains:

$$ \ddot{q}_k + 2\xi_k \bar{\omega}_k \dot{q}_k + \bar{\omega}_k^2 q_k + \sum_{i,j,l=1}^{K} \left[ \bar{\Gamma}_{ijl}^k q_l q_j q_l + \bar{\Pi}_{ijl}^k q_l (q_j \ddot{q}_l + \dot{q}_j \dot{q}_l) \right] - \frac{\bar{\omega}_k}{\sqrt{\bar{m}_k}} \left( \kappa_k + \sum_{i,j=1}^{K} \frac{\bar{\kappa}_{ij}^k}{2} q_l q_j \right) \bar{V} \quad (B.7a) $$

$$ + u_0 \left( \sum_{l=1}^{K} \bar{\Lambda}_l^k q_l + \sum_{i,j,l=1}^{K} \bar{\Lambda}_{ijl}^k q_l q_j q_k \right) = -F_k $$

$$ \bar{Q} - \bar{V} - \sum_{i=1}^{K} \sqrt{\bar{m}_i} \bar{\omega}_i \left( \kappa_i q_i - \frac{1}{6} \sum_{j,l=1}^{K} \bar{\kappa}_{ij}^k q_l q_j q_k \right) = 0, \quad (B.7b) $$

with

$$ \bar{\Gamma}_{ijl}^k = \frac{1}{\bar{m}_k} \int_0^1 \bar{D} \Phi_l' \Phi_l'' (\Phi_l' \Phi_k')' d\bar{x}, \quad (B.8a) $$

$$ \bar{\Pi}_{ijl}^k = \frac{1}{\bar{m}_k} \int_0^1 \Phi_k' \Phi_l' \int_1^{\bar{x}} \bar{m} \int_{\bar{x}}^0 \Phi_j' \Phi_k' d\bar{x} d\bar{x}, \quad (B.8b) $$

$$ \bar{\Lambda}_l^k = \frac{1}{\bar{m}_k} \int_0^1 \Phi_k' \Phi_l' \int_1^{\bar{x}} \bar{m} d\bar{x} d\bar{x}, \quad (B.8c) $$

$$ \bar{\Lambda}_{ijl}^k = \frac{1}{\bar{m}_n} \int_0^1 \Phi_k' \Phi_l' \Phi_j' \Phi_l' \int_1^{\bar{x}} \bar{m} d\bar{x}, \quad (B.8d) $$

$$ S_k = \frac{1}{\bar{m}_k} \int_0^1 \bar{m} \Phi_k d\bar{x}, \quad (B.8e) $$

$$ F_k = S_k \ddot{w}_0, \quad (B.8f) $$

$$ \kappa_k = \frac{\alpha}{\bar{\omega}_k \sqrt{\bar{m}_k}} [\Phi_k' \Phi_k^{\prime\prime}]_{\bar{x}_{-}}^{\bar{x}_{+}}, \quad (B.8g) $$

$$ \bar{\kappa}_{ij}^k = \frac{\alpha}{\bar{\omega}_k \sqrt{\bar{m}_k}} [\Phi_k' \Phi_l' \Phi_j' \Phi_l']_{\bar{x}_{-}}^{\bar{x}_{+}}. \quad (B.8h) $$

In the above equations, the $\sqrt{\bar{m}_k}$ factor in $\kappa_k$ was chosen to obtain a coefficient independent of the modal scaling. In particular, if the deformed shape is multiplied by a factor $a$, the modal mass $\bar{m}_k$ is multiplied by $a^2$ according to Eq. (B.6). All the other coefficients are dependent on the modal scaling of the mode shapes.

Since each mode has its own natural frequency, the above scaling of the time was chosen according to the time constant $1/\omega_0$ defined by Eq. (B.3). As a consequence, the dimensionless frequency $\omega_l$ appears in Eqs. (B.7a). In the case of the one mode expansion (on the ith mode) of Section 2.3.1, the time scaling is chosen as $1/\omega_l$ to further simplify the equations. The corresponding coefficients are then:

$$ \bar{\Gamma} = \frac{\bar{\Gamma}_{iii}^i}{\bar{\omega}_i^2} = \frac{2}{\bar{m}_i \bar{\omega}_i^2} \int_0^1 \bar{D} (\Phi_i' \Phi_i'')^2 d\bar{x} \quad (B.9a) $$

$$ \bar{\Pi} = \bar{\Pi}_{iii}^i = \frac{1}{\bar{m}_i} \int_0^1 \Phi_i^{\prime 2} \int_1^x \bar{m} \int_x^0 \Phi_i^{\prime 2} d\bar{x} d\bar{x}, \quad (B.9b) $$

$$ \bar{\Lambda} = \bar{\Lambda}_i^l = \int_0^1 \Phi_i^{\prime 2} \int_1^x \bar{m} d\bar{x} d\bar{x}, \quad (B.9c) $$

$$ S = \bar{S}_i = \frac{1}{\bar{m}_i} \int_0^1 \bar{m} \Phi_i d\bar{x}, \quad (B.9d) $$

$$ \kappa = \kappa_i = \frac{\alpha}{\bar{\omega}_i \sqrt{\bar{m}_i}} [\Phi_i' \Phi_i^{\prime\prime}]_{\bar{x}_{-}}^{\bar{x}_{+}}, \quad (B.9e) $$

$$ \bar{\kappa} = \bar{\kappa}_{ii}^l = \kappa_i \frac{[\Phi_i^{\prime 3}]_{\bar{x}_{-}}^{\bar{x}_{+}}}{[\Phi_i']_{\bar{x}_{-}}^{\bar{x}_{+}}}. \quad (B.9f) $$

Notice that changing the time scaling changes also the scaling of $V$ and $Q$.

## Appendix C. Effect of parametric cubic excitation term

A cubic nonlinear parametric excitation term $\bar{\Lambda}q^3\dot{u}_0$ naturally appears in Eq. (30), after the third order truncation. Available references including this term are very scarce. To the knowledge of the authors, it can only be found in [91], which deals with parametric amplification and attenuation and claims that it is two orders of magnitude smaller than comparable terms, which is contrary to our results (see Table 2).
<figure><img src="image_26.png" /><figcaption>Fig. C.1. Frequency response curves of a single layer cantilever beam obtained with Eq. (35b) for 3 different values of $\tilde{\Lambda}$: no cubic parametric term ($\tilde{\Lambda} = 0$), 10% of its calculated value ($\tilde{\Lambda} = 0.4043$) and its calculated value ($\tilde{\Lambda} = 4.043$). The HBM-H200 results (solid lines) are compared with numerical integration results (black circles) and the reference solution obtained with the exact geometrically beam model (red squares).</figcaption></figure>

The effects of this term on the PEH response are shown in Fig. C.1 for a single layer cantilever beam (First row of Table 2). Eq. (35a) with $\kappa = \tilde{\kappa} = 0$ is solved with numerical continuation (with $H = 200$ harmonics, HBM-H200) and the results are compared with numerical integration and the reference solution obtained with the exact geometrically beam model (FEHBM). The HBM-H200 and numerical integration are very similar, with a slight difference at high amplitude related to the fact that only the first harmonic amplitude is plotted for the HBM-H200 results, whereas the maximum amplitude over one period is considered for time integration. However, those two solutions can be considered equivalent.

First, considering $\tilde{\Lambda} = 0$, it appears that the truncated and the geometrically exact model strongly coincide. On the contrary, increasing the value of $\tilde{\Lambda}$ also increases the value of amplitude, leading to dramatic discrepancies. This unusual phenomenon may be explained by the one mode truncation, the effect of the other mode having potentially a compensating effect. Investigations on this question are left for future work. In the present work, this cubic term is neglected, the corresponding results being convincing enough within the validity range of the model.

## Appendix D. Manlab coding of the solved system

The HBM-H30 numerical computations are realised with the software Manlab, which requires a first order dynamical system with quadratic only nonlinearities [76–78]. This is obtained by adding auxiliary variables and additional algebraic equations to the dimensionless Eqs. (35a,c):

$$ \dot{q} = v \qquad (D.1a) $$

$$ \dot{v} = a \qquad (D.1b) $$

$$ \dot{Q} = I \qquad (D.1c) $$

$$ 0 = a - 2\xi v - (1 + \kappa_2)v - (\bar{\Gamma} + 2c/3\kappa_2)g - \bar{\Pi}(wa + rq) + \kappa Q + c\kappa/2h + f - p \qquad (D.1d) $$

$$ 0 = \bar{\tau}_e I + Q - \kappa q - c\kappa/6g \qquad (D.1e) $$

$$ 0 = w - q^2 \qquad (D.1f) $$

$$ 0 = g - wq \qquad (D.1g) $$

$$ 0 = r - v^2 \qquad (D.1h) $$

$$ 0 = h - wQ \qquad (D.1i) $$

$$ 0 = I_2 - I^2 \qquad (D.1j) $$

$$ 0 = \kappa_2 - \kappa^2 \qquad (D.1k) $$

$$ 0 = x - \delta \sin 2\bar{\Omega}\bar{t} \qquad (D.1l) $$

$$ 0 = p - qx \qquad (D.1m) $$

$$ 0 = f - F \sin \bar{\Omega}\bar{t} \qquad (D.1n) $$
$$0 = \tilde{P}_h - \bar{\tau}_e I_2 \qquad (D.1o)$$

$$0 = \tilde{P}_d - 2\xi r \qquad (D.1p)$$

$$0 = \tilde{P}_e - f\nu \qquad (D.1q)$$

$$0 = \tilde{P}_p - p\nu \qquad (D.1r)$$

where $[q, v, Q]$ are the main variables, $\kappa$ may be a parameter or a variable (see below), $[a, I, w, g, r, h, I_2, \kappa_2, x, p, f]$, $[\tilde{P}_h, \tilde{P}_d, \tilde{P}_e, \tilde{P}_p]$ are auxiliary variables and $(\xi, \tilde{\Gamma}, \tilde{\Pi}, \bar{\tau}_e, c, F, \delta)$ are parameters, with $c = \bar{\kappa}/\kappa$, defined in Table 2.

For the direct and parametric excitation cases, $\delta = 0$ and $F = 0$ are set, respectively. Then, a power correcting factor $P_{corr}$ is adjusted to compute $\tilde{P} = P_h/P_{max} = \tilde{P}_h/P_{corr}$, as defined by Eqs. (66) and (67). Considering that $\tilde{P}_h = \bar{\tau}_e \tilde{Q}^2$ in Eq. (D.1m) and the definition of the dimensionless variables (34), $P_{corr} = F^2/(16\xi)$ in direct excitation and $P_{corr} = S^2\delta^2/\bar{\Lambda}^2/(16\xi)$ in parametric excitation, such that $\tilde{F}(t) \equiv S\bar{\delta}(t)/\bar{\Lambda}$ leads to the same acceleration amplitude $\gamma$ (see Eq. (52)).

Several kinds of simulations are performed:

* Firstly, standard frequency responses can be obtained by leaving free $\bar{\Omega}$ and keeping $\kappa$ constant, leading to Figs. 6, 7 (with $\kappa = 0$) and Figs. 10, 11, 13, 14(a,b);

* Secondly, backbone curves of Fig. 5 are obtained by considering Hamiltonian free oscillator periodic solutions, i.e. only Eqs. (D.1a,b,d,f-h) without forcing ($F = \delta = 0$) and with the damping term $2\xi\dot{q}$ replaced by $\lambda\dot{q}$, with $\lambda$ a free bifurcation parameter, this latter term being used for numerical stability, to equilibrate an additional phase condition (see [78]). For short circuit, $\kappa = 0$, and for open circuit, the two terms $\kappa Q + c\kappa h/2$ are cancelled and $\kappa$ is a fixed parameter;

* Thirdly, less standard simulations are also performed following the power phase resonance conditions observed in Section 4.2, Sections 4.3 and 4.4. For the plots of Figs. 9(c), 12(b,c), 14(c,d), 16(b,c), F (for DPEH) or $\delta$ (for PPEH) are set to a constant values but $\bar{\Omega}$ and $\kappa$ are left free, the system being closed by a phase condition, such that the first harmonic of $q(t)$ has a prescribed phase lag $\alpha$ with $\tilde{F}(t)$ and $\bar{\delta}(t)$. For a DPEH, $\alpha = \pi/2$ and for a PPEH, $\alpha = -\pi/4$ or $3\pi/4$. In this case, $\kappa$ is considered as a variable and the system (D.1) must also be quadratic in $\kappa$, leading to the definition of the auxiliary variable $\kappa_2 = \kappa^2$. In this context, simulations with constant $\bar{\tau}_e$ help obtain sections $P_h = f(\kappa)$ of the horseshoe surfaces efficiently.

For the surface plots of Figs. 11(c) and 14(c,d), since the step-length of the arclength continuation is automatically adapted, some post-processing interpolations are performed with scattered data interpolation numerical algorithms (like griddata in Matlab or Python).

## Appendix E. Power extremum in parametric excitation

Differentiating Eq. (59) with respect to $(a_0, a_h, \tilde{a}_h)$ and dividing the results by $4\xi\omega_1\Omega$, $M_H$ of Eq. (61) can be split into two submatrices:

$$\begin{bmatrix} M_H^o & 0 \\ 0 & M_H^e \end{bmatrix} \begin{bmatrix} a^o \\ a^e \end{bmatrix} = 0, \quad \Rightarrow \quad \begin{cases} M_H^o a^o = 0 \\ M_H^e a^e = 0 \end{cases} \qquad (E.1)$$

with $a^o = [a_1, a_3, \tilde{a}_3, \dots, a_H, \tilde{a}_H]^T$ and $a^e = [a_0, a_2, \tilde{a}_2, a_4, \tilde{a}_4, \dots, a_H, \tilde{a}_H]^T$, that collects respectively the odd and even harmonics coefficients, and $M_H^o, M_H^e$ defined as:

$$M_H^o = \begin{bmatrix} 1 & Y & 0 & Y & 0 & 0 & 0 & 0 & \dots \\ Y & 1 & -Y & 0 & 0 & 0 & 0 & 0 & \dots \\ 0 & -Y & 3^2 & 0 & 0 & Y & 0 & 0 & \dots \\ Y & 0 & 0 & 3^2 & -Y & 0 & 0 & 0 & \dots \\ 0 & 0 & 0 & -Y & 5^2 & 0 & 0 & Y & \dots \\ 0 & 0 & Y & 0 & 0 & 5^2 & -Y & 0 & \dots \\ 0 & 0 & 0 & 0 & 0 & -Y & 7^2 & 0 & \dots \\ 0 & 0 & 0 & 0 & Y & 0 & 0 & 7^2 & \dots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix} \qquad (E.2)$$

and

$$M_H^e = \begin{bmatrix} 0 & 0 & 2Y & 0 & 0 & 0 & 0 & \dots \\ 0 & 2^2 & 0 & 0 & Y & 0 & 0 & \dots \\ 2Y & 0 & 2^2 & -Y & 0 & 0 & 0 & \dots \\ 0 & 0 & -Y & 4^2 & 0 & 0 & Y & \dots \\ 0 & Y & 0 & 0 & 4^2 & -Y & 0 & \dots \\ 0 & 0 & 0 & 0 & -Y & 6^2 & 0 & \dots \\ 0 & 0 & 0 & Y & 0 & 0 & 6^2 & \dots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix} \qquad (E.3)$$

There diagonals comes from the $P_d$ terms and are filled with $h^2$, $h$ being the number of the corresponding harmonic. Then, $P_p$ is responsible of the off diagonal terms, adding $4 \times 4$ antidiagonals $[Y, -Y, -Y, Y]$. As a consequence, two independent problems are obtained, one for the odd harmonics and one for the even harmonics.
<figure><img src="image_27.png" /><figcaption>Fig. E.1. Roots Y of det($M_H^o$) = det($M_H^e$) = 0 as a function of H.</figcaption></figure>

With $H \ge 3$ odd, solving det($M_H^o$) = 0 with symbolic computation gives a polynomial in $Y^2$ of order $(H+1)/2$, whereas with $H \ge 2$ even, det($M_H^o$) = 0 gives $Y=0$ and a polynomial in $Y^2$ of order $H/2-1$. Each polynomial has roots giving particular values of $Y$. Thus, gathering the contributions of both odd and even harmonics gives $H-1$ values of $Y$ that render Eq. (61) singular. The polynomial are obtained by symbolic computations and the numerical values of their roots, and thus of $Y$, are shown in Fig. E.1.

## Data availability

Data will be made available on request.

## References

[1] M. Shirvanimoghaddam, K. Shirvanimoghaddam, M.M. Abolhasan, M. Farhangi, V.Z. Barsari, H. Liu, M. Dohler, M. Naebe, Towards a green and self-powered internet of things using piezoelectric energy harvesting, IEEE Access 7 (2019) 94533-94556, http://dx.doi.org/10.1109/ACCESS.2019.2928523.

[2] N. Sezer, M. Koç, A comprehensive review on the state-of-the-art of piezoelectric energy harvesting, Nano Energy (80) (2021) 105567, http://dx.doi.org/10.1016/j.nanoen.2020.105567.

[3] A.J. Williams, M.F. Torquato, I.M. Cameron, A.A. Fahmy, J. Sienz, Survey of energy harvesting technologies for wireless sensor networks, IEEE Access 9 (2021) 77493-77510, http://dx.doi.org/10.1109/ACCESS.2021.3083697.

[4] A. Sadaf, R. Ahmed, H. Ahmed, Cantilever configurations in vibration-based piezoelectric energy harvesting: a comprehensive review on beam shapes and multi-beam formations, Smart Mater. Struct. 33 (2024) 123001, http://dx.doi.org/10.1088/1361-665X/ad8b87.

[5] M. Safaei, H.A. Sodano, S.R. Anton, A review of energy harvesting using piezoelectric materials: state-of-the-art a decade later (2008-2018), Smart Mater. Struct. 28 (2019) 113001, http://dx.doi.org/10.1088/1361-665X/ab36e4.

[6] H. Liang, G. Haoa, O.Z. Olszewskia, A review on vibration-based piezoelectric energy harvesting from the aspect of compliant mechanisms, Sensors Actuators A: Phys. 331 (2021) 112743, http://dx.doi.org/10.1016/j.sna.2021.112743.

[7] A. Erturk, D. Inman, Piezoelectric Energy Harvesting, John Wiley & sons, inc., USA, 2011, http://dx.doi.org/10.1002/9781119991151.

[8] M. F. Daqaq, R. Masana, A. Erturk, D.D. Quinn, On the role of nonlinearities in vibratory energy harvesting: A critical review and discussion, Appl. Mech. Rev. 66 (4) (2014) 040801, http://dx.doi.org/10.1115/1.4026278.

[9] Y. Jia, Review of nonlinear vibration energy harvesting: Duffing, bistability, parametric, stochastic and others, J. Intell. Mater. Syst. Struct. 31 (7) (2020) 921-944, http://dx.doi.org/10.1177/1045389X20905989.

[10] D. Gibus, A. Morel, P. Gasnier, A. Ameye, A. Badel, High performance piezoelectric vibration energy harvesting by electrical resonant frequency tuning, Smart Mater. Struct. 31 (12) (2022) 125012, http://dx.doi.org/10.1088/1361-665X/ac9d74.

[11] S. Leadenham, A. Erturk, Nonlinear M-shaped broadband piezoelectric energy harvester for very low base accelerations: primary and secondary resonances, Smart Mater. Struct. 24 (2015) 055021, http://dx.doi.org/10.1088/0964-1726/24/5/055021.

[12] J. Chen, B. Bao, J. Liu, Y. Wu, Q. Wang, Piezoelectric energy harvester featuring a magnetic chaotic pendulum, Energy Convers. Manage. 269 (2022) 116155, http://dx.doi.org/10.1016/j.enconman.2022.116155.

[13] C.R. McInnes, D.G. Gorman, M.P. Cartmell, Enhanced vibrational energy harvesting using nonlinear stochastic resonance, J. Sound Vib. 318 (2008) 655-662, http://dx.doi.org/10.1016/j.jsv.2008.07.017.

[14] A. Morel, L. Charleux, Q. Demouron, A. Benhemou, D. Gibus, C. Saint-Martin, A. Carré, E. Roux, T. Huguet, A. Badel, Simple analytical models and analysis of bistable vibration energy harvesters, Smart Mater. Struct. 31 (10) (2022) 105016, http://dx.doi.org/10.1016/j.ymssp.2024.112013.

[15] Y. Fan, M.H. Ghayesh, T.-F. Lu, A proof-of-concept study on an internal-resonance-based piezoelectric energy harvester with coupled three-dimensional bending-torsion motions, J. Vib. Acoust. 144 (2022) 061004, http://dx.doi.org/10.1115/1.4055720.

[16] H. Fu, X. Mei, D. Yurchenko, S. Zhou, S. Theodossiades, K. Nakano, E.M. Yeatman, Rotational energy harvesting for self-powered sensing, Joule 5 (2021) 1074-1118, http://dx.doi.org/10.1016/j.joule.2021.03.006.

[17] M.F. Daqaq, C. Stabler, Y.Q.T. Seuaciuc-Osório, Investigation of power harvesting via parametric excitations, J. Intell. Mater. Syst. Struct. 20 (5) (2009) 545-557, http://dx.doi.org/10.1177/1045389X08100978.

[18] A. Abdelkefi, A.H. Nayfeh, M.R. Hajj, Effects of nonlinear piezoelectric coupling on energy harvesters under direct excitation, Nonlinear Dynam. 67 (2) (2012) 1221-1232, http://dx.doi.org/10.1007/s11071-011-0064-9.

[19] M.A. Karami, J.R. Farmer, D.J. Inman, Parametrically excited nonlinear piezoelectric compact wind turbine, Renew. Energy 50 (2013) 977-987, http://dx.doi.org/10.1016/j.renene.2012.07.037.

[20] K. Mam, M. Peignay, D. Siegert, Finite strain effects in piezoelectric energy harvesters under direct and parametric excitations, J. Sound Vib. 389 (2017) 411-437, http://dx.doi.org/10.1016/j.jsv.2016.11.022.
[21] A. Garg, S.K. Dwivedy, Piezoelectric energy harvester under parametric excitation: A theoretical and experimental investigation, J. Intell. Mater. Syst. Struct. 31 (4) (2020) 612-631, http://dx.doi.org/10.1177/1045389X19891523.

[22] P. Alevras, S. Theodossiades, H. Rahnejat, Broadband energy harvesting from parametric vibrations of a class of nonlinear mathieu systems, Appl. Phys. Lett. 110 (23) (2017) 233901, http://dx.doi.org/10.1063/1.4984059.

[23] F. Fang, G. Xia, J. Wang, Nonlinear dynamic analysis of cantilevered piezoelectric energy harvesters under simultaneous parametric and external excitations, Acta Mech. Sin. 34 (3) (2018) 561-577, http://dx.doi.org/10.1007/s10409-017-0743-y.

[24] G. Xia, F. Fang, M. Zhang, Q. Wang, J. Wang, Performance analysis of parametrically and directly excited nonlinear piezoelectric energy harvester, Arch. Appl. Mech. 89 (10) (2019) 2147-2166, http://dx.doi.org/10.1007/s00419-019-01568-3.

[25] G. Xia, F. Fang, Q. Wang, M. Zhang, J. Wang, Performance analysis of piezoelectric energy harvesters with a tip mass and nonlinearities of geometry and damping under parametric and external excitations, Arch. Appl. Mech. 90 (10) (2020) 2297-2318, http://dx.doi.org/10.1007/s00419-020-01721-3.

[26] G. Xia, X. Kang, C.W. Lim, Y. Liu, D. Chen, L. Tang, T. Han, Parametric excitation analysis for system performance of piezoelectric energy harvesters, Appl. Math. Model. 21 (2023) 321-338, http://dx.doi.org/10.1016/j.apm.2023.05.003.

[27] Y. Fan, M.H. Ghayesh, T.-F. Lu, M. Amabili, Design, development, and theoretical and experimental tests of a nonlinear energy harvester via piezoelectric arrays and motion limiters, Int. J. Non-Linear Mech. 142 (2022) 103974, http://dx.doi.org/10.1016/j.ijnonlinmec.2022.103974.

[28] Y. Fan, M.-Q. Niu, M.H. Ghayesh, M. Amabili, L.-Q. Chen, Nonlinear vibration energy harvesting via parametric excitation: Snap-through with time-varying potential wells, Mech. Syst. Signal Process. 220 (2024) 111625, http://dx.doi.org/10.1016/j.ymssp.2024.111625.

[29] A. Garg, S.K. Dwivedy, Dynamic analysis of piezoelectric energy harvester under combination parametric and internal resonance: a theoretical and experimental study, Nonlinear Dynam. 101 (4) (2020) 2107-2129, http://dx.doi.org/10.1007/s11071-020-05931-w.

[30] Y. Jia, J. Yan, K. Soga, A.A. Seshia, A parametrically excited vibration energy harvester, J. Intell. Mater. Syst. Struct. 25 (3) (2014) 278-289, http://dx.doi.org/10.1177/1045389X13491637.

[31] M.F. Daqaq, D. Bode, Exploring the parametric amplification phenomenon for energy harvesting, Proc. Inst. Mech. Eng. Part I: J. Syst. Control. Eng. 225 (4) (2011) 456-466, http://dx.doi.org/10.1177/2041304110401145.

[32] U. Nabholz, L. Lamprecht, J.E. Mehner, A. Zimmermann, P. Degenfeld-Schonburg, Parametric amplification of broadband vibrational energy harvesters for energy-autonomous sensors enabled by fieldinduced striction, Mech. Syst. Signal Process. 139 (2020) 106642, http://dx.doi.org/10.1016/j.ymssp.2020.106642.

[33] O. Thomas, Understanding, computing and identifying the nonlinear dynamics of elastic and piezoelectric structures thanks to nonlinear modes, in: C. Touzé, A. Frangi (Eds.), Model Order Reduction for Design, Analysis and Control of Nonlinear Vibratory Systems, Springer Series CISM courses and lectures, New York, NY, 2024.

[34] G. Kirchhoff, Ueber das gleichgewicht und die bewegung eines unendlich dünnen elastischen stabes, J. Für Die Reine Und Angew. Math. 56 (1859) 285-313.

[35] E. Reissner, On one-dimensional finite-strain beam theory: The plane problem, J. Appl. Math. Phys. 23 (1972) 795-804, http://dx.doi.org/10.1007/BF01602645.

[36] C. Meier, A. Popp, W.A. Wall, Geometrically exact finite element formulations for slender beams: Kirchhoff-Love theory versus simo-Reissner theory, Arch. Comput. Methods Eng. 26 (2019) 163-243, http://dx.doi.org/10.1007/s11831-017-9232-5.

[37] O. Thomas, A. Sénéchal, J.-F. Deü, Hardening/softening behavior and reduced order modeling of nonlinear vibrations of rotating cantilever beams, Nonlinear Dynam. 86 (2) (2016) 1293-1318, http://dx.doi.org/10.1007/s11071-016-2965-0.

[38] H. Farokhi, Y. Xia, A. Erturk, Experimentally validated geometrically exact model for extreme nonlinear motions of cantilevers, Nonlinear Dynam. 107 (2022) 457-475, http://dx.doi.org/10.1007/s11071-021-07023-9.

[39] H. Farokhi, E. Kohtanen, A. Erturk, Extreme parametric resonance oscillations of a cantilever: An exact theory and experimental validation, Mech. Syst. Signal Process. 196 (2023) 110342, http://dx.doi.org/10.1016/j.ymssp.2023.110342.

[40] M. Debeurre, A. Grolet, B. Cochelin, O. Thomas, Finite element computation of nonlinear modes and frequency response of geometrically exact beam structures, J. Sound Vib. 548 (2023) 117534, http://dx.doi.org/10.1016/j.jsv.2022.117534.

[41] M. Debeurre, A. Grolet, B. Cochelin, O. Thomas, Extreme nonlinear dynamics of cantilever beams: effect of gravity and slenderness on the nonlinear modes, Nonlinear Dynam. 11 (14) (2023) 12787-12815, http://dx.doi.org/10.1007/s11071-023-08637-x.

[42] M.R.M.C. da Silva, C.C. Glynn, Nonlinear flexural-flexural-torsional dynamics of extensional beams. II. Forced motions, J. Struct. Mech. 6 (4) (1978) 449-461, http://dx.doi.org/10.1080/03601217808907349.

[43] O. Thomas, J.-F. Deü, J. Ducarne, Vibrations of an elastic structure with shunted piezoelectric patches: efficient finite element formulation and electromechanical coupling coefficients, Internat. J. Numer. Methods Engrg. 80 (2) (2009) 235-268, http://dx.doi.org/10.1002/nme.2632.

[44] A. Sénéchal, O. Thomas, J.-F. Deü, Optimization of shunted piezoelectric patches for vibration reduction of complex structures - application to a turbojet fan blade, in: Proceedings of the ASME 2010 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, IDETC/CIE 2010, Montreal, Canada, 2010.

[45] D. Dezest, O. Thomas, F. Mathieu, L. Mazenq, C. Soyer, J. Costecalde, D. Remiens, J.-F. Deü, L. Nicu, Wafer-scale fabrication of self-actuated piezoelectric nanoelectromechanical resonators based on lead zirconate titanate (PZT), J. Micromech. Microeng. 25 (3) (2015) 035002, http://dx.doi.org/10.1088/0960-1317/25/3/035002.

[46] D. Saya, D. Dezest, A.J. Welsh, F. Mathieu, O. Thomas, T. Leïchlé, S. Trolier-McKinstry, L. Nicu, Piezoelectric nanoelectromechanical systems integrating microcontact printed lead zirconate titanate films, J. Micromech. Microeng. 30 (2020) 035004, http://dx.doi.org/10.1088/1361-6439/ab60bf.

[47] J.N. Reddy, Mechanics of Laminated Composite Plates, CRC Press, 1997, http://dx.doi.org/10.1201/b12409.

[48] J. Ducarne, O. Thomas, J.-F. Deü, Placement and dimension optimization of shunted piezoelectric patches for vibration reductions, J. Vib. Control 331 (14) (2012) 3286-3303, http://dx.doi.org/10.1016/j.jsv.2012.03.002.

[49] A. Lazarus, O. Thomas, J.-F. Deü, Finite elements reduced order models for nonlinear vibrations of piezoelectric layered beams with applications to NEMS, Finite Elem. Anal. Des. 49 (1) (2012) 35-51, http://dx.doi.org/10.1016/j.finel.2011.08.019.

[50] A. Givois, J.-F. Deü, O. Thomas, Dynamics of piezoelectric structures with geometric nonlinearities: A non-intrusive reduced order modelling strategy, Comput. Struct. 253 (2021) 106575, doi:10.1016/j.compstruc.2021.106575.

[51] V. Guillot, A.T. Savadkoohi, C.-H. Lamarque, Analysis of a reduced-order nonlinear model of a multi-physics beam, Nonlinear Dynam. 97 (2019) 1311-1401, http://dx.doi.org/10.1007/s11071-019-05054-x.

[52] V. Guillot, A. Givois, M. Colin, O. Thomas, Theoretical and experimental investigation of a 1:3 internal resonance in a beam with piezoelectric patches, J. Vib. Control 26 (13-14) (2020) 1119-1132, http://dx.doi.org/10.1177/1077546320910536.

[53] S. Fallahpasand, M. Dardel, Piezoelectric energy harvesting from highly flexible cantilever beam, Proc IMechE Part K: J Multi-Body Dyn. 233 (1) (2019) 71-92, http://dx.doi.org/10.1177/1464419318791104.

[54] L. Shang, C. Hoareau, A. Zilian, A geometrically nonlinear shear deformable beam model for piezoelectric energy harvesters, Acta Mech. 232 (12) (2021) 4847-4866, http://dx.doi.org/10.1007/s00707-021-03083-5.

[55] A. Abdelkefi, A.H. Nayfeh, M.R. Hajj, Global nonlinear distributed-parameter model of parametrically excited piezoelectric energy harvesters, Nonlinear Dynam. 67 (2) (2012) 1147-1160, http://dx.doi.org/10.1007/s11071-011-0059-6.
[56] M. Géradin, A. Cardona, Flexible Multibody Dynamics: A Finite Element Approach, John Wiley & Sons, Ltd, Chichester, 2001, http://dx.doi.org/10.1016/0045-7949(94)90300-X.

[57] R.W. Ogden, Non-Linear Elastic Deformations, Dover, 1984, http://dx.doi.org/10.1007/BF00047287.

[58] P. Wriggers, Nonlinear Finite Element Methods, Springer, 2008, http://dx.doi.org/10.1007/978-3-540-71001-1.

[59] G.A. Holzapfel, Nonlinear Solid Mechanics: A Continuum Approach for Engineering, Wiley, 2001, http://dx.doi.org/10.1023/A:1020843529530.

[60] ANSI/IEEE Std 176-1987, IEEE Standard on Piezoelectricity, The Institute of Electrical and Electronics Engineers, Inc., 1988.

[61] H.F. Tiersten, Linear Piezoelectric Plate Vibrations, Springer US, Boston MA, 1969.

[62] U.V. Wagner, P. Hagedorn, Piezo-beam systems subjected to weak electric field: experiments and modelling of nonlinearities, J. Sound Vib. 256 (5) (2002) 861-872, http://dx.doi.org/10.1006/jsvi.2002.5024.

[63] S. Leadenham, A. Erturk, Unified nonlinear electroelastic dynamics of a bimorph piezoelectric cantilever for energy harvesting, sensing, and actuation, Nonlinear Dynam. 79 (3) (2015) 1727-1743, http://dx.doi.org/10.1007/s11071-014-1770-x.

[64] V. Denis, M. Jossic, C. Giraud-Audine, B. Chomette, A. Renault, O. Thomas, Identification of nonlinear modes using phase-locked-loop experimental continuation and normal form, Mech. Syst. Signal Process. 106 (2018) 430-452, http://dx.doi.org/10.1016/j.ymssp.2018.01.014.

[65] A. Givois, C. Giraud-Audine, J.-F. Deü, O. Thomas, Experimental analysis of nonlinear resonances in piezoelectric plates with geometric nonlinearities, Nonlinear Dynam. 102 (2020) 1451-1462, http://dx.doi.org/10.1007/s11071-020-05997-6.

[66] D. Gibus, P. Gasnier, A. Morel, N. Garraud, A. Badel, Non-linear losses study in strongly coupled piezoelectric device for broadband energy harvesting, Mech. Syst. Signal Process. 165 (2022) 108370, http://dx.doi.org/10.1016/j.ymssp.2021.108370.

[67] O. Thomas, J. Ducarne, J.-F. Deü, Performance of piezoelectric shunts for vibration reduction, Smart Mater. Struct. 21 (1) (2012) 015008, http://dx.doi.org/10.1088/0964-1726/21/1/015008.

[68] B. Lossouarn, L. Rouleau, R. Darleux, J.-F. Deü, Comparison of passive damping treatments based on constrained viscoelastic layers and multi-resonant piezoelectric networks, J. Struct. Dyn. 1 (2021) 30-48, http://dx.doi.org/10.25518/2684-6500.63.

[69] C.T. Johnk, Engineering Electromagnetic Fields and Waves, Wiley, 1988, http://dx.doi.org/10.1007/978-3-319-90847-2.

[70] E. Cottanceau, O. Thomas, P. Véron, M. Alochet, R. Deligny, A finite element/quaternion/asymptotic numerical method for the 3D simulation of flexible cables, Finite Elem. Anal. Des. 139 (2018) 14-34.

[71] O. Thomas, B. Legrand, C. Fuinel, Optimization of length and thickness of smart transduction layers on beam structures for control and m/NEMS applications, in: Proceedings of SMASIS 2015 (ASME 2015 Conference on Smart Materials Adaptive Structures and Intelligent Systems, Colorado Springs, USA, 2015, Paper 8857.

[72] B. Lossouarn, Optimal placement of thin piezoelectric patches for maximum modal coupling, J. Intell. Mater. Syst. Struct. 35 (20) (2024) 1591-1605, http://dx.doi.org/10.1177/1045389X241282054.

[73] P. Polytec, Piezoelectric ceramic products, 2008, URL https://www.piceramic.com/en/knowledge-center/downloads.

[74] A.H. Nayfeh, D.T. Mook, Nonlinear Oscillations, John Wiley & sons, inc., New-York, 1979, http://dx.doi.org/10.1002/9783527617586.

[75] J.J. Thomsen, Vibrations and Stability, McGraw-Hill, Denmark, 1997, http://dx.doi.org/10.1007/978-3-030-68045-9.

[76] L. Guillot, B. Cochelin, C. Vergez, A taylor series-based continuation method for solutions of dynamical systems, Nonlinear Dynam. 98 (2019) 2827-2845, http://dx.doi.org/10.1007/s11071-019-04989-5.

[77] L. Guillot, A. Lazarus, O. Thomas, C. Vergez, B. Cochelin, Manlab 4.0: An Interactive Path-Following and Bifurcation Analysis Software, Tech. rep., CNRS, Laboratoire de Mécanique et d'Acoustique, 2018, URL http://manlab.lma.cnsrs-mrs.fr.

[78] L. Guillot, A. Lazarus, O. Thomas, C. Vergez, B. Cochelin, A purely frequency based floquet-hill formulation for the efficient stability computation of periodic solutions of ordinary differential systems, J. Comput. Phys. 416 (2020) 109477, http://dx.doi.org/10.1016/j.jcp.2020.109477.

[79] S. Benacchio, C. Giraud-Audine, O. Thomas, Effect of dry friction on a parametric nonlinear oscillator, Nonlinear Dynam. 108 (2) (2022) 1005-1026, http://dx.doi.org/10.1007/s11071-022-07233-9.

[80] E. Abboud, A. Grolet, H. Mahé, O. Thomas, On the combination of multiple harmonics in direct and parametric forcing of an oscillator: new results on parametric amplification, Meccanica 57 (2022) 2447-2471, http://dx.doi.org/10.1007/s11012-022-01583-9.

[81] M. Debeurre, S. Benacchio, A. Grolet, C. Grenat, C. Giraud-Audine, O. Thomas, Phase resonance testing of highly flexible structures: measurement of conservative nonlinear modes and nonlinear damping identification, Mech. Syst. Signal Process. 215 (2024) 111423, http://dx.doi.org/10.1016/j.ymssp.2024.111423.

[82] I. Mahboob, H. Yamaguchi, Bit storage and bit flip operations in an electromechanical oscillator, Nature Nanotechnology 3 (2008) 275-279, http://dx.doi.org/10.1038/nano.2008.84.

[83] M. Aghamohammadi, V. Sorokin, B. Mace, Nonlinear dynamics of parametrically excited cantilever beams with a tip mass considering nonlinear inertia and duffing-type nonlinearity, Nonlinear Dynam. 111 (8) (2023) 7251-7269, http://dx.doi.org/10.1007/s11071-023-08236-w.

[84] H. Yabuno, Y. Ide, N. Aoshima, Nonlinear analysis of a parametrically excited cantilever beam: Effect of the tip mass on stationary response, JSME Int. J. Ser. C Mech. Syst. Mach. Elements Manuf. 41 (1998) 555-562, http://dx.doi.org/10.1299/jsmec.41.555.

[85] A. Brenes, A. Morel, J. Juillard, E. Lefeuvre, A. Badel, Maximum power point of piezoelectric energy harvesters: a review of optimality condition for electrical tuning, Smart Mater. Struct. 29 (3) (2020) 033001, http://dx.doi.org/10.1088/1361-665X/ab6484.

[86] J.M. Renno, M.F. Daqaq, D.J. Inman, On the optimal energy harvesting from a vibration source, J. Sound Vib. 320 (1-2) (2009) 386-405, http://dx.doi.org/10.1016/j.jsv.2008.07.029.

[87] Y. Liao, H. Sodano, Optimal power, power limit and damping of vibration based piezoelectric power harvesters, Smart Mater. Struct. 27 (7) (2017) 075057, http://dx.doi.org/10.1088/1361-665X/aabf4a.

[88] A. Brenes, A. Morel, D. Gibus, C.-S. Yood, P. Gasnier, E. Lefeuvre, A. Badel, Large-bandwidth piezoelectric energy harvesting with frequency-tuning synchronized electric charge extraction, Sensors Actuators A: Phys. 302 (2020) 111759, http://dx.doi.org/10.1016/j.sna.2019.111759.

[89] M. Marszal, B. Witkowski, K. Jankowski, P. Perlikowski, T. Kapitaniak, Energy harvesting from pendulum oscillations, Int. J. Non-Linear Mech. 94 (2017) 251-256, http://dx.doi.org/10.1016/j.ijnonlinmec.2017.03.022.

[90] P. Sommermann, M.P. Cartmell, Dynamics of an omnidirectional pendulum energy harvester: A comparative analysis between numerical and experimental results, Int. J. Non-Linear Mech. 159 (2024) 104588, http://dx.doi.org/10.1016/j.ijnonlinmec.2023.104588.

[91] V. Kumar, J.K. Miller, J.F. Rhoads, Nonlinear parametric amplification and attenuation in a base-excited cantilever beam, J. Sound Vib. 330 (22) (2011) 5401-5409, http://dx.doi.org/10.1016/j.jsv.2011.06.006.