A geometrically nonlinear shear deformable beam model for piezoelectric energy harvesters

An electromechanical model for beam-like piezoelectric energy harvesters based on Reissner’s beam theory is developed in this paper. The proposed model captures first-order shear deformation and large displacement/rotation, which distinguishes this model from other models reported in the literature. All governing equations are presented in detail, making the associated framework extensible to investigate various piezoelectric energy harvesters. The weak formulation is then derived to obtain the approximate solution to the governing equations by the finite element method. This solution scheme is completely coupled, and thus allows for two-way interaction between mechanical and electrical fields. To validate this model, extensive numerical examples are implemented in the linear and nonlinear regime. In the linear limit, this model produces results in excellent agreement with reference data. In the nonlinear regime, the large amplitude response of the piezoelectric beam induced by strong base excitation or fluid flow is considered, and the comparison of results with literature data is encouraging. The ability of this nonlinear model to predict limit cycle oscillations in axial flow is demonstrated.


Introduction
The growing demand for small-sized and low-power electronic devices has led to a focused research effort on the technology of energy harvesting, by which a permanent and autonomous power generator is possible due to the extraction of a usable form of energy from ambient energy sources. The locally available energy sources include mechanical, thermal, solar energy, and so forth. They can be converted to electrical energy by a particular transduction mechanism of either fundamental physical interactions (e.g., Faraday's law, from magnetic to electric) or material properties (e.g., piezoelectric media, from mechanical to electrical).
The interest of this study is piezoelectric energy harvesters (PEHs). One of the commonly investigated configurations of PEHs is a cantilevered beam made of a flexible and conducting material as a substrate layer with one (unimorph) or two (bimorph) piezoelectric layer(s) and a harvesting circuit attached. Similar laminated piezoelectric structures may also be used as sensors/actuators [1,2], if the main focus is on control applications. For energy harvesting applications, the external excitation source can be a vibrating host structure [3][4][5] or fluid flow [4,[6][7][8][9]. In the latter case, significant attention has been paid to PEHs in axial flow [8,10], where the effective and sustained extraction of flow energy is expected from limit cycle oscillations (LCOs). LCOs are large-amplitude self-excited and self-limiting vibrations resulting from, e.g., aeroelastic nonlinearity [10].
Accurate mathematical modeling is crucial to the design and optimization of PEHs. To date, lumped parameter [11] and distributed parameter [12] models have been developed in literature, and detailed discussion of issues and corrections of these models can be found in [13]. Among distributed parameter beam models, Euler-Bernoulli beam assumptions together with small displacement/rotation are most commonly used. It is well known that Euler-Bernoulli beam theory works perfectly for thin beams. However, for thick beams with low length-to-thickness aspect ratio, which leads to non-negligible rotary inertia effects and shear deformation [14], the Timoshenko beam theory would be a better choice. In 2010, Dietl et al. [15] proposed a Timoshenko beam model for a cantilevered bimorph PEH using the methodology of force and moment balance. Based on this paper, Zhu et al. [16] performed a more concrete analysis of the difference between the two PEH models. It is reported that at low length-to-thickness aspect ratio, the Euler-Bernoulli model tends to over predict the resonance frequency as well as the electrical output of the energy harvester. Energy methods are also adopted to build a Timoshenko beam model for PEHs, such as the work by Erturk [14]. More recently, Zhao et al. [17] also presented a Timoshenko beam model for a unimorph PEH. The unique contribution of this investigation consists in the solution: The steady-state Green's function method and Laplace transform method are firstly used to obtain the closed-form solution to the electromechanical PEH model. All these models are geometrically linear and thus only feasible in small deformation regimes.
When it comes to large deformation regimes, for instance, harvesting energy from LCOs, as mentioned above, geometrical nonlinearity of the structure must be taken into account. The most popular nonlinear structural model in this field is the inextensible beam model [8,10,18,19]. This model is derived from Euler-Bernoulli beam assumptions [20]. Therefore, it also suffers from the limitations indicated above in [16]. Further, since the above investigations [14][15][16][17] were only performed in the linear regime, more efforts are needed to explore the difference between PEH models based on Euler-Bernoulli assumptions and Timoshenko assumptions in the nonlinear regime.
Apart from analytic approaches with closed-form solutions generally used in the aforementioned research, the approximate solutions using the finite element (FE) method are also proposed, e.g., models based on Euler-Bernoulli assumptions by [6,21,22], a linear three-dimensional (3D) model by [23]. Finite element modeling is particularly valuable when studying the electromechanical behavior of complex PEHs because closed-form solutions are only available to PEHs with simple configurations [23].
This paper presents a nonlinear model of PEHs based on the geometrically exact beam theory [24]. With exact kinematic relations, this model incorporates first-order shear deformation and allows for large displacement/rotation. These two benefits are validated in the linear regime and the nonlinear regime, respectively.

General governing equations
In this Section, the equations governing the electromechanical behavior of a beam-like PEH are presented. For the sake of brevity, only a symmetric bimorph with a parallel-connection circuit is considered, but the basic idea of this method can be extended to consider PEHs in other configurations, as seen in the various numerical examples in Sect. 4 and the discussion on partial coverage of electrodes in Sect. 2.3.

Strong form equations
The geometrically exact beam theory dealing with the plane deformation of an originally straight beam [24] is used. The plane deformation may be stretching, bending, and shearing, which is linked to the axial strain ε, curvature κ, and shear strain γ , respectively. The cross section of the composite beam is rectangular and assumed to perform only rigid body motions during the deformation.
The electrodes are assumed to continuously cover the entire top and bottom surfaces of the upper and lower piezoelectric layers. In this case, the electric potential on each surface does not change in the axial direction of the beam, which is also called the equipotential condition [23]. The mechanical effects and the resistance of the electrodes are negligible.
A global Cartesian coordinate system O − XY Z is used to describe beam deformation, and the reference configuration A is for the undeformed beam while the current configuration A * is for the deformed beam. For each layer, the axis through the geometric center of the cross section is selected as the x-axis. The position of an arbitrary point P on each layer after deformation is determined by the axial displacement u, the transverse displacement w, and the rotation angle of the beam cross section ψ, as seen in Fig. 1a. The voltage output of the PEH, denoted by electric potential φ, is the voltage across the resistive load R in the circuit, as seen in Fig. 1b.
The mechanical variables u, w, ψ are spatially and temporally dependent while the electric variable φ is only temporally dependent due to the equipotential condition. From the perspective of mechanics, on account of the feature that the bimorph has only one dominant dimension, it is advantageous with little loss of precision to use the fiber aligned with the x-axis to represent the corresponding layer. Accordingly, the variation of mechanical quantities along the beam thickness is not considered unless otherwise specified. As a result, a set of 10 variables is under consideration: nine mechanical unknowns U m = {u i (x, t), w i (x, t), ψ i (x, t)} with ∀i ∈ {u, s, l} and one electric unknown U e = {φ(t)}, where and hereafter the subscripts (·) u , (·) s , (·) l indicate quantities related to the upper piezoelectric layer, the substrate structure, and the lower piezoelectric layer, respectively. When there is no need to distinguish the upper and lower piezoelectric layers, the subscripts (·) u and (·) l can be reduced to (·) p .
The three layers of the bimorph have the same length L and width b, and the thickness is denoted by h u , h s , h l (h u = h l = h p due to the symmetry). Both the piezoelectric and the substrate materials are assumed to have homogeneous mechanical properties with density ρ u , ρ s , ρ l (ρ u = ρ l = ρ p due to the symmetry). For this symmetric bimorph, the x-axis of the substrate structure coincides with the neutral axis of the composite beam.

Kinematic relations
According to the geometrically exact beam theory, the apparent deformation-induced strains (also called beam strain measures in classical beam theories [25], or more specifically, Reissner's generalized strains [26]) in the reference configuration are where (·) refers to the differential with respect to (w.r.t.) the local coordinate x. Because the fiber aligned with the x-axis of each layer is employed as the representative of the corresponding layer and thus z = 0, the beam strain measures of interest can be written as for ∀i ∈ {u, s, l}. In the above expressions, no assumptions are made for the rotation angle ψ i in addition to that the beam cross section itself is rigid (equivalent to a first-order shear deformation theory). For the piezoelectric layers, the non-negligible electric field E 3 can be seen as the electric counterpart to the mechanical strain. E 3 is assumed to be uniform across the thickness of the piezoelectric layers and is given in terms of φ [15] as It is pertinent to mention here that the relations between the electric field and the output voltage are dependent on the specific configurations of PEHs, see, e.g., for more details on the expression of E 3 in [3] for the series connection case. Another point worth attention is that, as indicated in [27], Maxwell's equations governing the electrodynamic behavior of a continuum usually refer to its current configuration, in which the physical quantities are measured. It means, Eq. (7) is stated w.r.t. the current configuration.

Constitutive relations
Linear material laws are used in this part to build the relations between the strain measures and the resultant loads for the PEH. For the substrate layer, the following linear relations in the reference configuration are commonly used [28]: where the stress resultants N s , V s , and M s are the normal force, shear force, and bending moment at the beam cross section, respectively; c Y,s and c G,s are Young's modulus or shear modulus; ν is the shear correction factor; A s is the area of the beam cross section; I s is the second area moment. For the piezoelectric layer, using the stress-electric displacement form of the linear constitutive equations for an X Z-plane Timoshenko beam model with the initial poling axis in 3-direction [29], the stress resultants in the reference configuration can be expressed as and additionally, the non-negligible electric displacement, D 3 , can be regarded as the electric counterpart to the mechanical force, given by for ∀i ∈ {u, l}. In Eqs. (9), (10), c E Y,i , c E G,i , S 33,i , e 31,i are material constants, and the superscripts (·) E , (·) S denote that the constants are evaluated at constant electric field and constant strain, respectively. Using Eq. (7) to replace E 3,i in Eqs. (9), (10) by φ, the stress resultants N i , V i , M i , D 3,i in the reference configuration w.r.t. variables {ε i , γ i , κ i , φ} with ∀i ∈ {u, l} can be formulated. It is noteworthy that Eq. (7) is defined on the current configuration while Eqs. (9), (10) on the reference configuration, so the above naive replacement leads to inconsistency in the configurations. This inconsistency is inconsequential in small deformation settings, and may not matter even in some large deformation cases, which can be inferred from [8,22], where this inconsistency is ignored in the theoretical model but good agreement with the experimental results is still achieved by the theoretical prediction results. To eliminate this inconsistency, one can refer to [27] to transform the electric quantities from the current configuration to the reference configuration.
Following the principle of virtual work as in [26], for ∀i ∈ {u, s, l}, the mechanical stress resultants in the current configuration are If the linear constitutive models are not appropriate (e.g., when dealing with large strains), one can select other constitutive models such as examples given in [25,26,30].

Equilibrium relations
The equilibrium relations are established at the deformed segment (the current configuration) by the force balance in 1-and 3-directions as well as the moment balance about the segment center point, which are (see [24]) with ∀i ∈ {u, s, l} in the dynamic context. In Eqs. (12), (14), f 1,i , f 3,i , and m 5,i are external distributed loads w.r.t. 1-, 3-, and 5-direction, respectively; the superscript( ·) means the second-order derivative w.r.t. time t, and (·) is still the derivative w.r.t. x. For the piezoelectric layer, according to the differential form of Gauss' law, the following equation: holds for ∀i ∈ {u, l} also in the current configuration. The area of the upper and lower surfaces of the piezoelectric layers after deformation is assumed to be identical to the initial area when integration of D 3,i is involved.

Boundary conditions
In line with the physical significance, the Dirichlet and the Neumann boundary conditions are called the kinematic and dynamic boundary conditions hereinafter, respectively. The mechanical kinematic boundary conditions are given by whereŪ m is a prescribed displacement or rotation angle in different directions as shown in Fig. 2a, and m is the corresponding kinematic boundary. The electric counterpart to the mechanical displacement, namely φ used in this work, is actually the difference between the electric potential (but not the electric potential itself) of the upper and lower electrodes attached on each piezoelectric layer, so no boundary values can be prescribed for φ, meaning there is no electric kinematic boundary condition. The mechanical dynamic equilibrium boundary conditions are for ∀i ∈ {u, s, l}, whereF 1 ,F 3 ,M 5 are prescribed external force or bending moment in different directions, [m] and[ J ] are the mass and the mass moment of inertia on the corresponding boundaries 1 i , 3 i , and 5 i , respectively, as depicted in Fig. 2b.
The electric equilibrium boundary conditions are for ∀i ∈ {u, l}, whereq is the prescribed free charge per unit electrode surface area.

Initial conditions
The mechanical and electric initial conditions are given as the zeroth-, first-, and second-order time derivatives of U m and U e : where α(x), β with the superscripts (·) O , (·) I , (·) II are a set of known time-independent functions. Only two time derivatives of U m and U e are required as the initial conditions for each unknown variable.

Coupling conditions
The three layers of the PEH are independent hitherto, and ten unknowns are concerned. By enforcing proper coupling conditions between them, the ten unknowns can be reduced to four: For the beam model under consideration, only the axial strain varies in the beam thickness direction as indicated by Eq. (1). Assuming that all strains at interfaces are continuous and the cross section of the whole laminate beam remains rigid during deformation, the following relations between mechanical unknowns of each layer can be obtained using Eqs. (1)- (3): where C k with ∀k ∈ {1, 2, 3, 4} are integration constants. Not knowing the values of C k is not problematic because only derivatives will be used in the final weak formulation. It should be noted that by this coupling the continuity of transverse shear stresses at interfaces computed from the constitutive laws is not satisfied. This issue may be alleviated by layerwise theories or zig-zag theories [31]. However, for thin and moderately thick laminated structures, the model resulting from the above coupling conditions is still able to predict adequately accurate global results, which is discussed further in Appendix A. From Gauss' law, the net charge in a parallel-connection circuit [15] is According to the definition of electric current and Ohm's law, the current in the PEH circuit, I C , has two different expressions as below, Combining them, the coupling condition of the circuit is 2.3 Discussion on partial coverage of electrodes For simplicity, in above derivation, full coverage of the electrodes is assumed. On account that the power output of PEHs is highly relevant to the electrode configuration due to the effect of strain nodes [32], it is beneficial here to make the above governing equations applicable to partially covered PEHs to track the state-of-the-art investigations [5,8].
With the help of Eqs. (1), (10), and (15), the following exact expression of E 3 can be formulated: where C 5 is the integration constant. Assuming that the coordinate of the interface between the covered and uncovered part is x = x p and the value of E 3 at the covered part is E 3 = E p , then the electric field E 3 of the If we stick to use the fiber aligned with the x-axis as the representative of the corresponding piezoelectric layer, which means z = 0 in Eq. (32), E 3 is still spatially independent, and the value of E 3 is identical to the value in the full coverage case. Therefore, all governing relations proposed above are still valid and only the integration interval has to be adapted to match the new length of the piezoelectric layer and the electrodes when deriving the weak formulation, which is presented in Sect. 3. This conclusion is consistent with the method used in other literature such as [5], where only the integration interval in the electrical equation is adapted when considering a bimorph with partial electrode coverage. Another interesting point of the exact expression E 3 (x, z) is that the distribution of E 3 is not uniform throughout the thickness of the piezoelectric layer, which is in contrast to the assumption about the electric field in the previous text (this assumption is commonly used in literature, such as [3,5,15]). If the variation of E 3 in the thickness direction is taken into account, the bending moment of the piezoelectric layer is S 33 I p κ p , so clearly the bending moment has one more term associated with the electric property of the material when compared with the simplified bending moment expression Eq. (9). Actually, if more complete constitutive laws for the piezoelectric materials are employed, as seen in paper [33], where the electric displacement in 1-direction is also involved, another term related to e 15 will also be present. The two terms including e 31 and e 15 together are called the induced electric bending moment, and may not be neglected in some special cases [33].

Summary
In this Section, a unified framework including kinematic, constitutive, equilibrium relations as well as boundary conditions and initial conditions is applied separately to the different layers of a beam-like PEH; in particular, for the piezoelectric layers, the electric quantities are treated as electric counterparts to mechanical quantities. The different layers of a PEH and the attached circuit finally form one interdependent whole by the coupling conditions. For convenience, the governing equations presented in this Section are simplified to a certain extent, which, in addition to making the electromechanical model simpler, also provides the space to improve the model precision in the future. As indicated above, the potential directions include: to select appropriate constitutive laws when large strains exist; to make the reference/current configuration of electric quantities consistent; to fulfill the continuity condition of transverse shear stresses at interfaces; to incorporate the bending moment contributed by the electric property. Nonlinearities activated by physical excitation phenomena [34,35] (the structural models of both are based on the Euler-Bernoulli assumptions) can be accounted for in the presented framework in order to study the overall behavior and performance of associated beam-like PEH systems.

Finite element model
Applying the method of weighted residuals to the equilibrium equations (12)- (15) and the boundary conditions (17)- (20), also introducing the electric coupling equation (31), with test functions matched with virtual work as {δε i , δγ i , δκ i , δφ}, for ∀i ∈ {u, s, l}, yields the weak formulation of the electromechanical system as follows: Equation (33) is consistent with the principle of virtual work (in the form of the principle of virtual displacement) with lines (33a-33e) devoted to the mechanical virtual work and lines (33f-33g) devoted to the electric virtual work. Line (33h) can be regarded as the method of Lagrange multipliers applied to Eq. (31). The last step to formulate the final weak formulation exclusively expressed with unknowns U = {u s (x, t), w s (x, t), ψ s (x, t), φ(t)} is to impose the mechanical coupling conditions Eqs. (23)-(28) on Eq. (33). The final weak formulation is not given in this paper, considering it is only concerned with simple mathematical substitution. Hereafter, the subscript (·) s linked to the displacement variables is omitted for convenience.
To obtain the approximate solution, the trial function for any unknown variable a(x, t) is introduced: In Eq. (34), n N is the number of ansatz functions, S a j (x) are ansatz functions of local or global support and fulfill the Dirichlet boundary conditions, q a j (t) are the associated displacement/voltage degrees of freedom, which are functions of time and are to be determined. In this work, the finite element (FE) method is employed, i.e., using ansatz functions of local support to discretize the electromechanical system in space. The semi-discrete matrix form can be written as: where q is the vector of displacement/voltage unknowns, q = (q u , q w , q ψ , q φ ) T ; M lin , C lin , and K lin are the mass matrix, damping matrix, and stiffness matrix associated with linear terms in Eq. (33), respectively; F nl (q,q, q) is associated with nonlinear terms, andq,q in F nl are induced by the coupling equations (23)- (28); f is the vector of externally applied load. Nonlinearity associated withq,q can normally be neglected, because these terms include a factor h s +h p 2 , which, for typical beam structures, is much smaller than other length measurements.

Model validation
This section is devoted to validate the model proposed in Sect. 3. All numerical simulations are implemented using the FEniCS framework [36,37]. The choice of function space is: continuous Galerkin of degree 2 for S u j (x) and S w j (x), continuous Galerkin of degree 1 for S ψ j (x), and S φ j (x) = 1. By this choice, shear locking effects are avoided [38]. The mass-proportional damping [23] is used in simulations when needed. The convergence of simulation results w.r.t. the number of finite elements and time step (in nonlinear problems) is examined for every example. All geometric and material parameters are given in Appendix B.

Linear regime
In this part, the linearized weak formulation of the model proposed in Sect. 3 is derived to facilitate the frequency domain analysis. Considering the main objective is to capture shear deformation effects, a numerical example aiming at thick beams (non-symmetric unimorph) given by [17] is implemented in addition to an experimentally validated example of thin beams (symmetric bimorph) given by [3].

Linearized weak formulation
If small deformation is concerned, the linear weak formulation of a symmetric bimorph in a parallel configuration can be achieved by linearizing the nonlinear kinematic relations Eqs.
with equivalent parameters being In Eq. (36), the terms of boundary conditions and external loads are omitted for convenience. It can be found that the equivalent parameters Eqs. (38)-(42) are the same as those presented in [15], while Eq. (37) is special to the present formulation. The expanded matrix form with sub-matrices resulting from FE discretization is ⎡ Equation (43) is a linear second-order ordinary differential equation (ODE) system that can be solved in frequency or time domain.

Numerical examples
Thin bimorph This example is from [3], where the single-mode electromechanical frequency response functions (FRFs) that relate the voltage output to translational base excitation are presented for a bimorph with a tip mass both in series and parallel configurations. All physical parameters in the numerical simulations are kept the same as the setup given in [3]. The thickness-to-length aspect ratio of the bimorph is (2h p +h s )/L = 1.3%. The response of voltage output φ in the frequency domain for 3 different resistive load values (R = 1 k , R = 33 k , R = 470 k ) is computed and then compared with the analytic results, as shown in Fig. 3, where the results of both parallel and series (in this case, ϑ = − 1 2 e 31 b(h s + h p ), C 0 = 1 2 S 33 bL/ h p ) configurations are given.
From Fig. 3, it can be seen that the numerical FRFs agree well with the analytic FRFs. Furthermore, the precise values of voltage output for different resistive loads at their respective resonant excitation frequencies are also compared, as shown in Table 1. Although discrepancies are observed between the voltage output from the present work and the experimental results from [3], the present work obtains almost the same results as that from the single-mode FRF given by [3] and the fully three-dimensional finite element model proposed by [23]. The fully three-dimensional model is apparently powerful when dealing with complicated PEHs, but the current beam model is particularly suitable for thin-walled PEHs because it saves much more computational resource compared to the three-dimensional model.

Thick unimorph
This example is from [17], where the first natural frequencies of a Timoshenko beam model-based unimorph are computed using the Green's function method and Laplace transform technique. The equivalent parameters of a unimorph can be obtained either directly from [12] or derived following the procedure presented in Sect. 2 with the neutral axis of the composite beam as the x-axis of the substrate structure. Table 2 shows the comparison of the first natural frequencies when the unimorph is moderately thick with different thickness-to-length aspect ratio under short-circuit (R = 1×10 2 ) and open-circuit (R = 1×10 6 ) conditions. The analytic solutions are the closed-form solutions based on the steady-state Green's function method and Laplace transform technology presented in [17], where data are available only under the shortcircuit condition. The 1D FE data are obtained from the linear beam model Eq. (36), while the 3D FE data are obtained from a continuum-based three dimensional model. It can be seen that good agreement is achieved in all comparative cases. Therefore, at this point, it is persuasive to say that the current electromechanical beam model is capable of capturing shear deformation effects.   [3] are calculated from the single-mode FRF, while all other reference voltage data are from [3,23]. Unavailable reference data are indicated by -

Nonlinear regime
In this part, two large deformation examples are implemented based on the nonlinear model proposed in Sect. 3. The first example is a highly flexible unimorph with a tip mass under strong base excitation that leads to an extreme deflection of the beam. The second example is a preliminary examination of [8], where encouraging prediction of the voltage output is achieved from a flow-driven bimorph although the fluid model is highly simplified. The generalized-α scheme with parameters ensuring unconditional stability and secondorder accuracy [39] is used to update fields in the dynamic system over time, and the nonlinear equations are solved by the Newton-Raphson method.

PEH under strong base excitation
This example is from [22], where the computational model is based on the Euler-Bernoulli assumptions under small strains and large deflections, and a decoupled solution scheme is employed, namely, solving the structural model firstly in a commercial FE program (ABAQUS) and then substituting the structural results into the circuit equation to obtain electric results. According to the specific circuit equation used in [22] and the circuit model of a unimorph with a resistive load R in [40], the circuit equation Eq. (31) in simulations is adapted to be where C p is the piezoelectric capacitance, φ oc is the open-circuit voltage across the piezoelectric layer, given by [22] where L p is the length of the piezoelectric layer (L p = 31 mm), h c is the distance of the top of the piezoelectric layer from the neutral axis of the composite beam, and ψ(L p ) is the rotation angle of the beam at L p . In the numerical simulations, all parameters of the setup are from Table 1 in Ref. [22]. Figure 4 shows the transverse displacement and the voltage at different excitation frequencies in the range 4.0-5.4 Hz when the base acceleration is 4 ms −2 . The reference data are taken from Fig. 9 in Ref. [22] (only data in the range 4.0-4.9 Hz are given), and the predicted data are the maximum steady responses in the simulations. Although there is a small phase shift between the reference and computed resonant frequency, i.e., from 4.64 Hz to 5.0 Hz, the values of both simulated tip displacement and voltage output coincide well with the reference data. The maximum tip displacement is about 25 mm, nearly 80% of the beam length (L = 32 mm). As for the Fig. 4 Comparison of reference and simulated displacement and voltage at different excitation frequencies voltage, the maximum value in simulations is 62 V, 7.5% lower than the reference voltage (67 V). This voltage discrepancy could be attributed to the different solution schemes used in our simulations and the reference paper [22]. When large deformation is concerned, the electromechanical models are normally highly nonlinear and thus difficult to obtain solutions, so decoupled solution schemes are commonly used, such as [8,22,41] (all of them use Euler-Bernoulli assumptions). In all numerical simulations of the present work, a completely coupled solution scheme is employed, i.e., solving the mechanical and electric unknowns simultaneously. By a numerical comparison, paper [22] indicates that the coupled scheme tends to predict lower voltage output than the decoupled scheme, especially when the electromechanical coupling is strong.
The dynamic response of this PEH at the resonant excitation frequency is shown in Fig. 5. The small saddles in Fig. 5c are caused by bending backward, which can be clearly seen from the deformed shapes of the unimorph (Fig. 5d). The unique bending backward phenomenon captured by nonlinear models in extremely large deflection cases is also reported recently in [41,42].
Additionally, a numerical test is made to show the difference between the electromechanical results predicted by the linear model Eq. (36) and the nonlinear model Eq. (33) under the base excitation condition. The thin bimorph in a parallel configuration presented in Sect. 4.1.2 is used. To ensure this PEH appropriately flexible to undergo large displacement/rotation, the thickness of the substrate structure and the piezoelectric layers is reduced by half. All other parameters, including the damping ratio, remain the same. In this test, the resistive load is 33 k , and the frequency of the base excitation is the linear resonant frequency when R = 33 k (17.2 Hz). The comparison is given in Fig. 6, where the normalized tip displacement (NTD) is the percentage of the maximum steady tip displacement over the beam length. It can be seen that when the tip displacement is about 20% of the beam length, the linear displacement prediction starts to deviate from the nonlinear prediction; as to the power output, the deviation point presented by NTD is about 30%. Fig. 6 illustrates that for a PEH operating in a large deformation regime, e.g., in this particular example, when the NTD is larger than 30%, it is necessary to employ a geometrically nonlinear model to obtain accurate electromechanical response.

PEH harvesting energy from axial air flow
In principle, to validate the proposed nonlinear model for flow-driven PEHs, e.g., harvesting energy from LCOs, an accurate fluid model for the axial flow is essential. Although considerable work has been accomplished in this area, as evident from the introduction sections of [8,43], it remains difficult to mathematically express and to numerically compute the fluid force acting on the deforming structure. Hence, detailed fluid modeling is excluded from the scope of the current paper.
Considering that the electric output of PEHs is directly related to the structural deformation, the electric output can be determined for a given deformation pattern of the piezoelectric beam, no matter what (external force) induces this deformation pattern. Accordingly, the basic idea of the model validation in this case is to manufacture the deformation pattern present in [8] and then compare the voltage output. Reproducing the deformation pattern to a highly precise extent is only possible in a limited way due to the absence of a proper fluid model as well as other physical parameters. However, effective comparison can still be expected if the main characteristics of the deformation pattern are preserved. It is reported in [8] that the modal content of the LCOs of the specimen is predominantly comprised of the first and second vibration modes; and the transverse displacement at 80% of the length from the clamped end is approximately 0.012 m (namely, w| x=0.8L = 0.012 m) when the airflow speed is 34 ms −1 (data estimated from Fig. 7c in Ref. [8]). These two points are taken as the main characteristics of the deformation pattern to be reproduced. A simplified linear fluid model [43] is employed to manufacture the above deformation pattern, which reads where p is the fluid pressure, ρ air is the airflow density, and U ∞ is the free-stream velocity in x-direction. A bimorph of continuous piezoelectric layer coverage in a series configuration with the same parameters as in [8] is used in the numerical simulations. Epoxy layers between the substrate and each piezoelectric layer are neglected. The external force model Eq. (46) is integrated in the nonlinear model Eq. (33), and an artificial initial velocity in the transverse directionẇ| t=0 = 0.1x is applied to start the simulations. Figure 7 shows the electromechanical response of this bimorph when U ∞ = 40 ms −1 and R = 10 M . At this flow velocity, Fig. 6 Comparison of linear and nonlinear PEH models under base excitation the flutter amplitude at x = 0.8L is 11.9 mm (Fig. 7a), the modal form involves mainly the second-order vibration modes (Fig. 7c), and the maximum voltage output in LCOs is 47.0 V (Fig. 7b). When R = 100 M , the corresponding voltage is 50.9 V. In [8], the voltage output of the specimen is between 30 V and 40 V (data estimated from Fig. 8c of Ref. [8]) when w| x=0.8L ≈ 0.012 m, R ∈ [10 M , 100 M ]. Therefore, the identical magnitude of voltage to the reference data is predicted by the present model. The difference between the exact values can be explained by the roughness of the fluid model used in the current simulations. With such a simplified fluid model, the manufactured deformation pattern is not the same as that in [8].
In this flutter case, the solution convergence of Eq. (33) is observed to deteriorate severely with the longitudinal inertia term ρ i A iüi δu i . Since the motion in this direction is not of interest, this term is omitted in the above numerical simulations. It is also noteworthy that no physical damping (e.g., material viscous damping) or numerical damping from the generalized-α scheme [39] is introduced into the above simulations. Therefore, the nonlinear restoring force, which is necessary to keep solutions bounded in time in the postcritical regime [44], as seen in Fig. 7, is due to the geometrical nonlinearity of the structural model. Contrasting this with a linear structural model such as Eq. (36), if the same linear fluid model Eq. (46) is used, the dynamic response of the system will be unbounded in time, and thus leading to non-physical results (Fig. 8).
Remark In the strict sense, the above validation is rather rough since the fluid model is highly simplified, but it is still clear that the proposed model has the potential to predict reasonable electromechanical response for a fluid-driven PEH. What's more, the significance of introducing nonlinearity in the structural model for the fluttering case is elucidated by the comparison with a linear model for the same setup.

Conclusions
The beam-like configuration is commonly used in piezoelectric energy harvesting devices. This article proposes a model for PEHs based on the geometrically exact beam theory and appropriate solutions to the governing equations using the finite element method. With the help of the exact kinematics, this model is able to capture both shear deformation and large displacement/rotation. Various numerical examples covering the symmetric bimorph/non-symmetric unimorph, series/parallel configuration, thin/thick beams, and (strongly) base-excited/flow-driven PEHs are implemented, and good results are obtained, clearly exhibiting the broad applicability of the framework presented in the current work. In the nonlinear regime, the comparative investigation of a PEH under base excitation shows that geometrical nonlinearity should be considered when the tip displacement is over 30% of the beam length; and for the flow-driven PEHs, geometrical nonlinearity alone of the structural model is sufficient to provide the nonlinear restoring force necessary for the occurrence of limit cycle oscillations in axial flow.

Appendix A
The aim of this part is to discuss the accuracy of the proposed beam model, which does not fulfill the stress compatibility between the layers due to the lack of warping displacements through the beam thickness. A static problem of a cantilevered laminate comprising one core layer and two face layers is investigated using the above beam model and a 2D continuum model. For simplicity, we choose only two different materials without piezoelectric effects, one for the core layer and the other one for the two face layers. The geometric and mechanical parameters are given in Table 3. The governing equations for a 3D continuum are where F = I +∇d is the deformation gradient tensor, S is the 2nd Piola-Kirchhoff stress tensor, E is the Green-Lagrange strain tensor, g is the body force vector, and λ, μ are Lamé constants. In the 2D implementation, plane strain assumptions are used; in the 1D (beam) implementation, the equivalent parameter (E I ) e in Eq. (38) is multiplied with the coefficient 1 1−η 2 to obtain the appropriate bending stiffness. We consider two laminates A and B, where each layer has a thickness of 0.01 m and 0.03 m, respectively, and thus the thickness-to-length ratio α is 3% (case A) and 9% (case B). For both A and B, the core layer may be the strong material (c Y = 100 GPa) or the weak material (c Y = 20 GPa). For laminate A, the static body force is applied in the transverse direction (3-direction) of the structure as g = (0; 3 × 10 7 ) Nm −3 , while g = (0; 5 × 10 8 ) Nm −3 for B.
The deflection of the centerline of laminate A and B is obtained from an FE solution of the proposed beam model and the above 2D model. The results are shown in Fig. 9.
In Fig. 9, the solid lines with circle markers are the 1D results when the shear correction factor (SCF) 5 6 is used. In Fig. 9b, the triangle and square markers indicate the SCF is 1 2 and 1, respectively. For laminate A, the difference between the 1D and 2D results cannot be distinguished, so only one SCF is employed. For laminate B, the SCF does influence the 1D results, although rather slightly in this case. This may be explained by the role of the SCF in the first-order shear deformation theory (FSDT). In general, FSDT fails to represent the high-order variation of transverse shear stresses through the thickness. To amend this, the SCF is introduced to match the global response obtained from FSDT with the elasticity solutions. For laminated structures, an ideal SCF could be computed using some sophisticated methods [45], where the configuration, geometry, and   mechanical parameters are all of interest. Therefore, to further improve the accuracy of the beam model, an optimal SCF should be carefully determined case by case.
We also compare the strain and stress distribution over the middle cross section of laminate B when a weak core is used, as shown in Fig. 10. The axial strain for the 2D model is the axial component of E, i.e., E 11 ; for the 1D model, the axial strain is ε(x, z). Similarly, the shear strain is E 13 in 2D and γ (x, z) in 1D; the axial and shear stresses in 2D are S 11 and S 13 , respectively; the axial stress in 1D is c Y ε(x, z), and the shear stress is c Y γ (x,z) 2(1+η) . It can be seen that the proposed geometrically nonlinear 1D model predicts perfect axial strain/stress for the setup under consideration. This provides a justification for the observed adequately accurate global response of the 1D model, as demonstrated by Fig. 9, although the correct shear strain/stress is missing. This point may be illustrated from the perspective of strain energy. Defining the component of strain energy as below, the percentage of axial strain energy W 11 and shear strain energy W 13 in the total energy W 2D = W 11 + W 13 + W 33 can be computed from the 2D model for a laminate with a weak core and various thicknesses, as shown in Table 4. For the 1D model, W 11 and W 13 can also be defined in a similar manner, e.g., replacing E 11 by ε(x, z), and the total energy is W 1D = W 11 + W 13 . The proportion of shear strain energy increases monotonically as the laminate becomes thicker, but for all the given parameters, i.e., α 12%, the axial strain energy accounts for more than 95% share (referring to the 2D results). Hence, the accuracy of the model for global response is dominated by the prediction of axial strain/stress when the laminate is thin or moderately thick. In summary, it is clear from the above discussion that the 1D model in this work can be used to obtain effective global response for common beam-like piezoelectric energy harvesters, where the structures are usually not very thick and the material parameters of the substrate and the piezoelectric patches do not deviate from each other too significantly.

Appendix B
Geometric and material parameters used in the numerical examples of this work are given in Tables 5, 6, 7, and 8, where the length is denoted by L, width by b, thickness by h, density by ρ, Young's modulus by c Y , Poisson ratio by η, piezoelectric constants by d 31 and e 31 , and permittivity by S 33 . Other parameters: Tip mass M t = 0.012 kg, 0 = 8.854 pF m −1 , mechanical damping ratio of the first vibration mode ζ 1 = 0.027, shear correction factor is 0.83 Other parameters: 0 = 8.854 pF m −1 , shear correction factor is 0.87. For the 3D FE implementation: The elasticity matrix, permittivity matrix, and piezoelectric moduli matrix are taken from Appendix B of [33], i.e., C 11    Other parameters: Shear correction factor is 0.87, air density ρ air = 1.225 kg m −3