Geometrical nonlinearities and shape effects in electromechanical models of piezoelectric bridge structures

We consider nonlinear shape effects appearing in the lumped electromechanical model of a bimorph piezoelectric bridge structure due to the interaction between the electromechanical constitutive model and the geometry of the structure. At finite proof-mass displacement and electrode voltage, the shape of the beams is no longer given by Euler-Bernoulli theory which implies that shape effects enter in both the electrical and mechanical domains and in the coupling between them. Accounting for such effects is important for the accurate modelling of, e.g., piezoelectrical energy harvesters and actuators in the regime of large deflections and voltages. We present a general method, based on a variational approach minimizing the Gibbs enthalpy of the system, for computing corrections to the nominal shape function and the associated corrections to the lumped model. The lowest order correction is derived explicitly and is shown to produce significant improvements in model accuracy, both in terms of the Gibbs enthalpy and the shape function itself, over a large range of displacements and voltages. Furthermore, we validate the theoretical model using large deflection finite element simulations of the bridge structure and conclude that the lowest order correction substantially improve the model, obtaining a level of accuracy expected to be sufficient for most applications. Finally, we derive the equations of motion for the lowest order corrected model and show how the coupling between the electromechanical properties and the geometry of the bridge structure introduces nonlinear interaction terms.


Introduction
The application of nonlinear mechanical structures in inertial energy harvesting using piezoelectric transduction is a growing field of research motivated by the potential advantages over linear resonator systems in terms of bandwidth, durability and amplitude [1][2][3][4][5]. In addition, manifestly nonlinear modes of operation have been suggested to target excitation regimes where linear systems typically produce insufficient power [1,[6][7][8]. Conversely, microelectromechanical (MEMS) piezoelectric actuators and switches are characterized by fast response and high bandwidth [9] and able to produce large amplitude actuation. Accurate operation in the large deflection regime requires accounting for nonlinear effects both in terms of the displacement and the voltage applied to produce it [10,11].
A common structure in applications, which features mechanical nonlinearities, is the doubly clamped bridge with a central proof-mass, as illustrated in Fig. 1. At infinitesimal deflections d of the proof-mass, the deformation of the beams is described by the Euler-Bernoulli theory of bending strain and the restoring force is linear in d. However, at finite deflections d, the boundary conditions and symmetry of the problem induce stretching strain in the beams, corresponding to nonlinear contributions to the restoring force. At deflections d which are large compared to the thickness of the beams, nonlinear effects must be accounted for to accurately describe the system [1,5]. In addition, the shape of the beams, i.e., the functional coordinate dependence of the displacement field, depends on both the proof-mass displacement d and the voltage v which must be accounted for.
The purpose of the present paper is the derivation and validation of a theoretical model describing the shape of the beams in the piezoelectric bridge structure as a function of the displacement d and voltage v, as well as a lumped model for the piezoelectric bridge structure, that is a system of differential equations governing the displacement d and voltage v. Such a theoretical model provides a way to understand how shape effects enter into both the static and dynamical properties of the system and is essential for accurately designing structures for various applications in the domain of large deflections and voltages.
A variational approach was used in [12] to compute large d corrections to the nominal shape function (valid at infinitesimal d) in a purely mechanical bridge structure. In the present paper, we extend the analysis to a piezoelectric bridge structure by including electromechanical effects in the description of the system. Consequently, the shape function acquires a dependence on the electric voltage v across the piezoelectric layer in addition to the deflection d. We construct a model for the electromechanical bridge based on a power series Ansatz for the shape function and perform a variational analysis to derive conditions corresponding to the minimization of the Gibbs enthalpy. We also compute explicitly the lowest order corrections to the nominal shape function and validate the resulting model using finite element model (FEM) simulations. Finally, we derive the equations of motion for the model including the lowest order correction by incorporating the electromechanical shape effects, and discuss their impact on the dynamics of the system. The variational approach results in an improved electromechanical model of piezoelectric bridge structures accounting for electromechanical shape effects, extending the work by Gafforelli et al. [6]. The importance of including the shape effects in order to accurately model the operation of bridge type harvesters, sensors and actuators in realistic applications where nonlinear effects are significant is demonstrated.
In addition to clearly displaying the the functional relationships between design parameters and nonlinear effects in the system, the lumped model offers significant advantages over FEM simulations in terms of computational cost and possibilities to explore parameter space. In particular, this is true for the simulation of system dynamics where the nonlinearities require computationally intensive timedependent FEM simulations where the time step required for accuracy is typically several orders of magnitude smaller than the time scale of the dynamics of the system and transient phenomena.
Lumped models of the bridge structure (as well as other common structures appearing in MEMS applications) have been derived before and have proven extremely useful [5,6,13]. In [6], a Rayleigh-Ritz method was employed to derive a nonlinear model of the piezoelectric bridge based on a time-dependent rescaling of the Euler-Bernoulli displacement profile valid a infinitesimal amplitudes. A similar approach was previously used to derive a nonlinear model of microcantilever actuators [14] and its extension to a nonlinear modal analysis was considered in [15,16]. Furthermore, the approach utilized in the derivation of our model has been successfully applied on numerous occasions in the form of a Galerkin method used to model geometrical nonlinearities in the context of efficient full spatio-temporal numerical simulations of MEMS structures [17,18]. However, to the best of our knowledge, the model presented in the present paper is the first lumped model of the piezoelectric bridge system accounting for geometrical nonlinearities and incorporating their effect on the shape of the bridge.

Electromechanical model
In this section, we consider a nonlinear electromechanical model of a piezoelectric bimorph operating in the e 31 mode, and the boundary conditions associated with a doubly clamped bridge structure. At large deflections d, the mechanical strain according to Green-Lagrange is given as the sum of a bending strain and a stretching strain where w 1 and w 3 are the displacement field in the x 1 -and x 3 -directions, respectively. The stretching strain, corresponding to axial forces in the harvester, is due to the moderate rotation effects and can be neglected for infinitesimal d where the bending strain dominates. For a device operating in the e 31 -mode, the electric potential is taken to depend only on the x 1 coordinate perpendicular to the piezoelectric layer and the bimorph electrodes, and the corresponding electric field is given by the gradient The constitutive piezoelectric relations in the e 31 mode are given by (3) Fig. 1 Illustration of the doubly clamped bridge structure with definitions of proof mass displacement d, beam length L and coordinates x 1 and x 3 where T 33 is the stress, D 1 is the displacement field, c 33 is the elastic modulus at constant electric field strength, 11 is the permittivity at constant mechanical strain and e 31 is the piezoelectric coupling coefficient 1 . The elastic modulus for the substrate and piezoelectric materials of the bimorph are denoted C and c, respectively.

Boundary conditions
The mechanical boundary conditions for a single beam are In addition, we assume that the w 1 displacement field is independent of both the x 1 and the x 2 coordinates. In this situation, the symmetry of the problem implies that the function w 1 is antisymmetric with respect to x 3 = L∕2 , which imposes additional constraints on the derivatives of the displacement w 1 according to where ℕ 1 is the set of all positive integers 2 . We will discuss the properties of the structure required for the assumption w 1 ∕x 2 = 0 to be valid in more detail below. Furthermore, symmetry of the bridge configuration preventing displacement of the proof mass in the x 3 -direction imposes zero x 3 -displacement conditions at x 3 = 0 and x 3 = L implying that the deformation w 1 is associated with a change in length of the harvester and that the inplane displacement field w 3 and its derivatives can be safely neglected.
The electrical boundary conditions depend on the details of the piezoelectric layer configuration. We will consider a bimorph configuration illustrated in Fig. 2. In order to utilize the stretching strain contribution dominating at large deflection, the polarized piezoelectric region and the corresponding electrode extend the full length of the beam. We note that using this electrode configuration, the bending strain contributions to the electromechanical energy will cancel identically.

Gibbs enthalpy
Under the symmetry conditions discussed above implying that w 3 and its derivatives can be neglected, the total Green-Lagrange strain reduces to We note that the stretching strain contribution is manifestly positive, meaning that at large deflection operation the harvester voltage will not alternate with the deflection d as in the case of bending strain harvesting at small deflections.
The electric potential is assumed to be linear in the x 1 -coordinate across the piezoelectric layer, giving a electric field where v is the electrode voltage and h is the thickness of the piezoelectric layer.
The internal potential energy of the system is the energy stored in the electroelastic field, which in the situation of constant E and S is given by the Gibbs enthalpy [19], where V refers to the combined volume of both beams and A to the cross-sectional area in the directions transverse to the The shape function To simplify the notation, we now drop the subscript on w 1 and introduce the dimensionless coordinate The functional form of the deflection profile of the beam is then described by the shape function w( |d, v) where we have made the dependence of w on the proof-mass deflection d and the electrode voltage v explicit. The constraints on w( |d, v) , corresponding to (5), (6) and (7) are and where w ′ and w (n) denote derivatives with respect to .
We will now proceed to compute expressions for the various parts of the Gibbs enthalpy in terms of dimensionless integrals of the shape function, which will subsequently be (11) used in the variational analysis in the new section. For brevity, the dependence of the shape function on and v will usually be left implicit below.

Shape integral formulation
Considering first the mechanical part of the Gibbs enthalpy, we can integrate out the dependence on the transverse directions and enforce the boundary condition (16) yielding where V b and V s are the shape integrals corresponding to bending and stretching We have also introduced the constants and containing all dependence on geometry and mechanical properties for the bending and stretching strain energies, respectively, in terms of the moment of area integrals.
The electromechanical contribution to the Gibbs enthalpy is similarly given by where the shape integral is and we have introduced the constant accounting for the geometry and electromechanical properties of the structure. As mentioned above, the symmetry of the bimorph configuration implies the identical cancellation of the contribution from the coupling of the bending strain to the electromechanical enthalpy. (20) Finally, integrating over all coordinate directions in the electrical contribution gives where we have introduced the constant for consistency in the treatment of the different parts of V . We recognize V E as (minus) the ordinary capacitive energy stored in the piezoelectric layers and note that it is independent of the shape function due to our assumption (9) that the electric potential is linear 3 .
Collecting the parts of the Gibbs enthalpy in terms of the dimensionless integrals containing its dependence on the shape function w, we have This formulation, where the dependence on the potential v is explicit and all the dependence on the proof-mass displacement d is contained in the dimensionless integrals will be convenient in the next section when we consider the variational analysis of the system which is the main focus of this paper.

Variational analysis
The defining property of the shape function w( ) is that it minimizes the Gibbs enthalpy V subject to the boundary conditions. The shape function that extremizes the functional V[w] can be found as the solution to the variational equation 4 .
The fact that V e is independent of the shape function (under our assumptions on the electrical field) implies that it will not contribute to the functional derivative.
We will use a power series Ansatz for the shape function in the dimensionless length parameter where the coefficients a n are determined by the constraints (15)-(17) on w and the minimization condition (28). In practice, we will truncate the power series at some power N, in order to simplify the problem and determine a suitable approximation for the true shape function. It should be noted that this is different from an ordinary perturbation theory expansion of a scalar quantity in terms of a small parameter, since is actually a coordinate and the minimization is performed over the entire domain. However, due to symmetry, the shape function is uniquely determined by its values over the interval ∈ [0, 1∕2] where higher powers of are suppressed. Furthermore, the approximate strain energy V is a decreasing function of N since as we add higher powers to w( ) a solution with unchanged energy is obtained by simply setting all added coefficients a n = 0 . Therefore, in the limit N → ∞ , we expect to approach the exact solution.
In the first part of this section, we will provide some general results and describe the method used to find the approximate shape function w for arbitrary N. Subsequently, we will exemplify the procedure for N = 5 which represents the first correction to the trivial (in a sense to be made precise below) shape function.

General method and results
Imposing boundary and symmetry constraints on the shape function power series results in dependencies among the coefficients a n . We first enforce the boundary conditions (15) and (16) which further reduces the number of independent coefficients a n . The requirement (28) that the strain energy functional be minimized by w, allows the determination of the remaining coefficients a * n through the relations We note the importance of enforcing the boundary conditions and symmetry constrains before applying the variational principle, since this guarantees that the procedure extremizes the strain energy in the space of shape functions compatible with the corresponding geometrical constraints.
Truncating the power series expression at some order N simply amounts to replacing the infinite sums with finite ones by setting a n = 0 for n > N . Since the shape function w is antisymmetric with respect to = 1∕2 , it is natural to restrict truncations to N odd. In particular, the boundary and symmetry conditions then reduce the number of independent coefficients to be determined by the variational Eq. (35) to (N + 1 − 4)∕2.

Lowest order shape function
From the above expressions, it is clear that the lowest order of truncation allowed in order to satisfy the boundary conditions (30) and (31) is N = 3 . Is this case (34) and (35) impose no additional constraints and the solution is simply given by In should be noted that (36) is the exact solution for a doubly clamped beam in the Euler-Bernoulli theory 5 of bending, applicable for infinitesimal deformations where bending provides the dominant contribution to the strain energy. The shape function (36), therefore, constitutes the lowest order shape function in the sense of accounting for large deformations of the doubly clamped beam.

N = 5 corrections
Truncating the general expressions for boundary and symmetry conditions at N = 5 , we can now compute the first corrections to (36). For N = 5 , there is only one independent coefficient, which we take to be a 5 = a . The boundary conditions and the n = 2 symmetry constraint allow us to express the dependent coefficients as  The purpose of computing corrections to the shape function is to obtain an improved estimate compared to the lowest order N = 3 function, which we have seen describes the pure bending situation at infinitesimal displacements and applied voltages and is therefore not expected to properly account for large deformations of the doubly clamped structure or large electrical fields. In the form derived above, the variational problem depends on the electromechanical properties and the geometry of the bridge structure.
In order to obtain results which are independent of the precise properties and geometry of the structure, we introduce dimensionless quantities where the characteristic length and voltage scales of the problem are given by We also introduce the dimensionless shape enthalpy which contains only shape-dependent terms which will contribute to the variation with respect to w. The variational condition determining the bridge shape is then finally given as This is a cubic equation which we can solve numerically for given the dimensionless displacement and voltage , yielding the results in Fig. 3 for different ranges of the dimensionless displacement and voltage. We note that since the enthalpy is not manifestly positive, the relative reduction of W in the N = 5 model compared to the N = 3 model is not a meaningful (or even regular) measure. Therefore, the dimensionless shape enthalpy reduction is used to quantify the improvement in model accuracy in Fig. 4. We note that in all cases the reduction is non-negative implying that the shape enthalpy W 3 in the N = 3 model is larger than the corresponding quantity W 5 in the N = 5 model.
As a special case, setting = 0 reproduces the results in the purely mechanical case [12]. In this case, the enthalpy reduces to potential strain energy and the relative reduction becomes well-defined. The parameter and the relative reduction in strain energy ΔW are show in Fig. 5.
It is straightforward to verify numerically that for the range of and considered here, the discriminant of the cubic Eq. (44) is negative, implying that the existence of a single real root which guarantees the uniqueness of the numerical solution.

Model validation
While it is clear from the results presented in the previous section that the N = 5 model provides a (substantial) improvement compared to the N = 3 model, it remains to validate the variational approach in comparison with the physical systems it is meant to describe, and to quantify the improvement obtained by including the N = 5 model corrections in this context. To this end, the theoretical models are evaluated for a validation model and compared to the

Simulation models
We consider a bridge structure with geometry defined by

Validation results
The primary quantification of the improvement obtained by accounting for the N = 5 corrections is the Gibbs enthalpy. As above, we are primarily concerned with the part that depends on the shape function, that is the shape enthalpy W . For the validation ranges specified above, the shape enthalpy W is shown as a function of both d and v in Fig. 6, and as a function of d for constant voltage v in Fig. 7. It is clear from the results that the N = 5 corrections compensate for essentially the entire enthalpy error caused by the incorrect assumption of the form of the shape function in the N = 3 model. Furthermore, the improvement is significant in the range considered in the validation model. For example, at d = 30 m and v = 200 mV, the error relative to the simulation is reduced from 30% in the N = 3 model to 2% in the N = 5 model.
The shape function w( ) can also be validated using the simulation model. Figure 8 shows results obtained from the theoretical models and the simulation model for proof-mass displacement d = 30 m. While it is not obvious how to quantify the improvement in terms of the shape function in a physically relevant way, the N = 5 model appears to produce a better qualitative agreement with the simulated results than does the N = 3 model. We also note that the FEM simulation validates the boundary conditions (5)- (7).

System dynamics
In this section, we incorporate the shape effects considered above into the equations of motion for the electromechanical bridge system, in order to account for their influence on the dynamics of the system. To accomplish this, we use the Euler-Lagrange equations which are conveniently formulated in suitable generalized coordinates of the system. The Euler-Lagrange approach yielding equations of motion for the lumped variables are applicable since the time scale of spatial strain and voltage equilibration in the structure is typically orders of magnitude smaller than the time scale of the change of global variables in applications.
For the bridge structure, the fundamental quantity in the electrical domain is the charge q separated across the piezoelectric layer. The rate of change in q is related to the voltage v across the layer by the continuity equation in the electrical domain, the precise form of which depends on the effective external impedance between the top and bottom electrode. We will consider the case where the circuit connecting the electrodes is a simple resistive load R, which is a common choice in, e.g., the energy harvesting literature and makes it easy to consider the open and closed circuit dynamics as the limiting cases R → ∞ and R → 0 , but a more general terminating circuit can be considered in an analogous way. The equation of continuity is then where the voltage difference between the top and bottom electrodes is 2v. Since it is not possible to absorb all Here, the Lagrangian is L = K − V , where V is the Gibbs enthalpy and K is the kinetic energy in the system, and P is the power dissipating from the system. For a purely resistive load, the only contribution to the kinetic energy is mechanical where M is the proof-mass and the kinetic energy in the beams is neglected. Further neglecting mechanical damping effects, the dissipated power P is simply

Equations of motion for N = 5
We have seen that the Gibbs enthalpy for the N = 5 model of the electromechanical bridge system takes the general form for the coupling between the domains. As for the static properties and the shape function, the dynamics of the N = 3 model is obtained by setting the model coefficient a = 0 in the N = 5 model. We note that this implies that while the N = 3 shape function is that of the Euler-Bernoulli theory, the dynamics of the N = 3 model includes the cubic term accounting for the inclusion of stretching in the model [1,5]. Furthermore, the N = 3 model reproduces the structure of the equations of motion derived in [6,8], numerical factors differing due to the configuration of the piezoelectric layers.

Dynamical shape effects
The inclusion of shape effects in the analysis of the bridge system causes two different kinds of modifications to the equations of motion in the N = 5 model compared to the N = 3 model. The first is a straightforward addition of terms accounting for mechanical and electromechanical couplings to the shape of the beams. We also note that the linear coefficient k 1 in the mechanical domain receives a correction from the stretching strain when a ≠ 0 . The second is the appearance in (56) of terms proportional to the derivative ȧ . In other words, the voltage v couples to the shape not only through the dependence of the electromechanical interactions on a, but directly to the rate of change a with respect to time.

Conclusion
In this paper, we have shown how nonlinear shape effects appear in the lumped electromechanical model of a bimorph piezoelectric bridge structure. Under broadly applicable assumptions regarding the form of the Green-Lagrange strain, these effects arise from the interaction between the electromechanical constitutive model and the geometry of the structure. At finite proof-mass displacement and electrode voltage, the shape of the beams minimizing the Gibbs enthalpy is no longer the solution of Euler-Bernoulli theory, which influences the coupling between the electrical and mechanical domains.
Adding higher order terms in to the Euler-Bernoulli solution, we use a variational approach to derive corrections to the shape function of any order. The lowest order N = 5 corrections are parametrized by a single dimensionless model coefficient and we derive explicitly the equilibrium equation determining as a function of the dimensionless displacement and voltage . Numerical solution of the governing Eq. (44) reveals an intricate nonlinear dependence of both the model coefficient and the (shape-dependent part of) the enthalpy on and . For a large range of displacements and voltages, we show that the corrections to the shape function produce a reduction in the shape enthalpy corresponding to an improved model accuracy.
Subsequently, we validate the theoretical model using large deflection FEM simulations of a bridge structure in the MEMS regime. The validation model and the simulated ranges of displacements and voltages are selected to be representative of both the energy harvesting (or sensing) and the actuation modes of operating the structure. The validation results indicate that the N = 5 model represents a substantial improvement of the lowest order N = 3 model. In fact, the accuracy achieved by the N = 5 model in terms of the Gibbs enthalpy and shape function is expected to be sufficient for most applications.
Finally, we derive explicitly the equations of motion for the N = 5 model and show how the coupling between the electromechanical properties and the geometry of the bridge structure introduces nonlinear interaction terms. The equations of motion follow directly from the Euler-Lagrange equations by inserting the expression for the N = 5 Gibbs enthalpy W , which has been independently validated.
As always, there are limitations associated with the assumptions made in the derivation of the model presented here. Such limitations include the restriction to a linear constitutive relation (3) and a linear relation (9) between the voltage and electric field as well as ignoring any spatial electromechanical fringe effects. In addition, no estimates of the convergence of the power series approximation as N increases are provided. While the numerical validation suggests that already the N = 5 model is sufficient for most applications, it is important to keep such theoretical limitations in mind, in particular when considering extending the scope of the model.
Care must also be taken extending the ranges of (dimensionless) displacement and voltage in the model. The discriminant of the third degree polynomial Eq. (44) can become non-negative for some values of and , which corresponds to multiple real roots, and multiple local extrema of the enthalpy, which must be investigated to determine the global minimum and the corresponding shape function.
The theoretical modelling of nonlinear piezoelectric structures can be motivated as a tool to make accurate predictions of the behaviour of structures appearing in applications to sensing, energy harvesting or actuation at a fraction of the computational cost of FEM models. However, perhaps the most important aspect of the theoretical lumped model is its potential to produce a deeper understanding of the system under investigation by identifying the underlying mathematical structures. More specifically, this includes the functional dependence of the coupling constants appearing in both the static and dynamical model on the geometrical and electromechanical parameters of the system. Such an understanding is essential for a successful and systematic design process in practical applications. In conclusion, the present paper elucidates the interaction between the geometry and the electromechanical properties of a nonlinear bridge structure and shows that this interaction must be accounted for in order to obtain accurate models in several common applications of piezoelectric devices.
Funding Open access funding provided by RISE Research Institutes of Sweden.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.