## Abstract

The focus of this investigation is first on assessing the validity to structures under in-plane forces, in particular near buckling, of the reduced order modeling approach for nonlinear geometric response that has been extensively developed in the last two decades. This focus is motivated by a class of piezoelectric energy harvesters that rely on strongly nonlinear behavior, such as large amplitude responses, to achieve broadband energy harvesting. A simple, two-rigid bars linkage that approximates a buckling beam is first considered to discover the features of the nonlinear force–displacement relationship induced by an in-plane loading. It is observed that the corresponding form of this relationship is not consistent with the one derived from a reduced order model (ROM) but can be closely approximated by it over a large displacement range. This analysis emphasizes in particular the role of a group of ROM coefficients that are usually considered unimportant. A similar study is performed next for the buckled harvester modeled within nastran and it is again found that a close match of the force–displacement relationship can be achieved. Based on that positive outlook, a six basis functions ROM of this beam harvester that includes piezoelectric effects is built and identified. It is found to provide a close match of nastran nonlinear predictions over a broad range of transverse and in-plane loadings in static and dynamic conditions. The ROM usefulness in predicting the open-circuit voltage is demonstrated.

## 1 Introduction

Piezoelectric energy harvesting from ambient vibrations has attracted extensive research interest as a potential solution to continuous and convenient energy supply for consumer electronics, wireless sensors, portable health monitors, etc. Earlier studies were focused on the harvesting of vibratory energy around a dominant natural frequency. This linear energy harvester performs well when the ambient vibration is narrow-banded close to the natural frequency of the harvester but is less efficient when the environmental vibration has a wide frequency spectrum. In recent years, a variety of harvester designs have been proposed that exhibit strongly nonlinear responses to broaden the effective bandwidth of energy harvesting [1–10].

Evaluating the efficiency of these nonlinear harvester designs requires the determination of their response under a variety of dynamic scenarios representative of the environment in which they operate. Such computations can certainly be carried out using finite element software with dynamic nonlinear solutions and piezoelectric capability, e.g., nastran (using the thermal analogy, e.g., see Refs. [11,12]) and abaqus, but such efforts can be quite time consuming. Accordingly, it would be desirable to have an alternative approach that is both accurate and computationally efficient, i.e., essentially the equivalent of modal analysis for linear systems in which the response involves fixed modes and time-varying modal coordinates.

Fortunately, such a direct extension already exists in the form of nonintrusive nonlinear reduced order models (ROMs) developed in the last couple of decades, see Ref. [13] for a somewhat dated review and Refs. [14–46] and references for other/newer investigations. These ROMs exhibit many features of modal models, e.g., they can be derived nonintrusively from finite element software, they express the response as a sum of components associated with different modes/basis functions, and the form of the governing equations is given with parameters depending on the structure and the modes selected. However, they involve more efforts in their construction as the basis involves more than the standard linear modes and the governing equations are nonlinear, but these properties are expected given the nonlinearity of the response.

The existing literature on these ROMs has addressed basic formulation issues of basis selection and coefficients identification [14–27], extensions of the methodology to (i) structures composed of substructures [28–31], (ii) structures with defects [32–36], (iii) to address design issues [26,37], to include (iv) piezoelectric effects [12,38], (v) viscoelastic damping [39], and (vi) the effects of temperature [40–46]. Applications/validations have been described in these papers and others, e.g., Refs. [47–62], ranging from simple beam/plate models to notably more complex structures, e.g., see Refs. [16,18,28,36,43,47,53,62]. Most of them focus on structures, often aircraft panels, that are subjected to transverse forces/pressures large enough to induce larger than linear responses. A few examples have been considered in the literature in which the response is strongly nonlinear, i.e., the snap-through of clamped–clamped curved beams [21,50,52,55,59,60], the thermal buckling and ensuing snap-through of straight beams [40,61], and a nonconventional buckling [53]. While using a somewhat different formulation, the shearing of a thin-walled cylinder has also been captured [63]. These ROMs have also demonstrated the capability to capture nonlinear features in the response such as symmetry breaking [56] and internal resonance [57]. Accordingly, they are potentially also applicable to harvester designs with strongly nonlinear features but this validity needs to be carefully confirmed.

The strongly nonlinear harvester design of interest here is the buckled beam introduced in Ref. [6] and shown in Fig. 1(a). The applied in-plane force *P* is large enough to induce the buckling of the beam and the ground excitation (representing the external vibration) leads to large displacements around the buckled position, potentially from one buckled position to the other. From the standpoint of energy harvesting, it would be very desirable to determine the response of the beam under dynamic loadings of different amplitudes, different frequency bands, and with different in-plane force *P* to assess its optimum value. These multiple parameters lead to a significant number of full nonlinear solutions that would have to be computed and thus such a study would strongly benefit from an accurate and expedient ROM.

However, it is unclear if a ROM would have these good properties as no prior application of this methodology has included in-plane forces, much less those large enough to induce buckling. In this regard, it should be recognized that thermal buckling and buckling through the application of a compressive force lead to mathematically different models. Temperature affects directly the linear stiffness of the ROM and buckling occurs when this matrix first achieves a zero eigenvalue. Mechanical buckling by an in-plane force is a much more complex process as this force affects first the in-plane displacement and is reflected in the transverse motions by appropriately capturing the nonlinear in-plane–transverse coupling.

In this light, the present investigation focuses first on an assessment of the potential for ROMs to capture the response of structures near buckling and then on an exemplification of their usefulness to study the energy harvested by the beam of Fig. 1 relying on a previously developed thermal analogy of piezoelectric effects. Prior to analyzing this structure, a simple illustrative example will be considered that provides a very useful perspective on the application range of the ROM. For completeness, a short review of the ROM construction is however first performed.

## 2 Reduced Order Model Formulation Review

In this equation, the vectors **Ψ**^{(n)} are the basis functions and *q*_{n}(*t*) are the corresponding generalized coordinates. The use of such a representation requires three key steps: (i) the derivation of the governing equations for the generalized coordinates, (ii) the selection of the basis, and (iii) the identification of the coefficients involved in the governing equations from an underlying finite element model. These three steps are reviewed in the ensuing sections.

### 2.1 Governing Equations.

*q*

_{n}(

*t*) is usually assumed to be nonlinear in stiffness only with that component being a cubic polynomial of these coordinates, i.e.,

*j*,

*l*, and

*p*in Eq. (2) are dummy indices ranging each from 1 to

*M*, see also discussion in Sec. 2.3. The terms

*M*

_{ij},

*D*

_{ij}, and $Kij(1)$ are recognized as the elements of the mass, damping (added to capture the dissipation effects), and linear stiffness matrices, while $Kijl(2)$ and $Kijlp(3)$ will be referred to as the quadratic and cubic stiffness coefficients. Note further that

*F*

_{i}denotes the

*i*th modal force originating from the external excitation pulled back to the undeformed configuration as necessary.

Equation (2) was derived in Ref. [64] from finite deformation elasticity under the assumption of a Saint-Venant/Kirchhoff constitutive behavior, i.e., with a linear relationship between the Green strain tensor and the second Piola–Kirchhoff stress tensor. It has also been derived [18] under the assumption of the von Karman strain definition.

The use of Eq. (2) has led to an excellent matching between ROM and full finite element predictions for a variety of structures, e.g., see Refs. [14–62]. However, the level of displacements and rotations in most of these cases has been modest, most often below four thicknesses peak transverse deflections. Notable exceptions to this situation have been the snap-through of curved clamped–clamped beam [21,50,52,55,59,60]—up to about ten thicknesses—and a study of cantilevered structures [19] in which the displacements were expectedly much larger and for which the ROMs have been successful for displacements up to 30% of span. For larger deflections of these same structures, the full finite element predictions (nastran) deviate from the ROM owing to differences in the nonlinear finite element formulation.

More importantly in the present context, none of the prior ROM investigations has considered the effects of an in-plane force, neither of small magnitude nor, as desired here, close to or exceeding the buckling limit.

### 2.2 Basis Functions.

The second key step in the formulation of the reduced order models is the selection of the basis functions **Ψ**^{(n)}. Assume first that only transverse loads act on the structure. Then, the corresponding dominant deformations are also transverse and have been found to be well represented by a set of the low-frequency linear modes. However, in-plane deformations will also arise due to the membrane stretching effect. These latter deformations are typically stiff and thus would be associated with high-frequency linear modes which could, with great care, be used to model them, see Ref. [13] and references therein. Other complements to the basis formed by the low frequency modes are (i) the dual modes of Refs. [13,15,16,21,22,62] for formulation details and Refs. [19,28,32,33,35,39–46,53,55–59] for additional validations and (ii) the modal derivatives of Refs. [20,21].

Based on the above discussion, the dual modes are constructed by exciting the structure along the linear modes retained in the basis and capturing the displacements that occur in addition to these modes when the structure is in the nonlinear regime. Given that the dual modes are expected to be associated with high frequencies, this analysis is done statically.

To induce a static response primarily along a linear mode **Ψ**^{(n)} retained in the basis, the structure is subjected to a force vector proportional to $KFE(1)\Psi (n)$ where $KFE(1)$ is the finite element linear stiffness matrix. If the structure responds linearly, then the displacements will be exactly proportional to **Ψ**^{(n)}. On the contrary, when nonlinearity is present, there will be an additional component to the response that includes both components along other modes in the basis resulting from nonlinear modal coupling and, more importantly, the membrane effects for that particular loading.

Note that the obtained dual modes are designed to capture the membrane effect, i.e., the in-plane displacements induced by the large transverse ones, not those resulting from in-plane loading. Accordingly, it may be necessary to also enrich the basis with functions that will properly capture the effects of the in-plane forces. If so, assuming again that the in-plane deformations are stiff, it is proposed here to select such enrichments as the static linear response to the in-plane loading.

### 2.3 Identification of Reduced Order Model Coefficients.

*j*≠

*l*cannot be distinguished from each other in Eq. (2) as they involve the same product of generalized coordinates. A similar situation occurs with cubic terms. Accordingly, when identifying the coefficients as below, it is sufficient to consider the range of the indices

*j*,

*l*, and

*p*in Eq. (2) to be 1 ≤

*j*≤

*M, j*≤

*l*≤

*M*, and

*l*≤

*p*≤

*M*, see Ref. [13]. Consider first the determination of the ROM mass and linear stiffness matrices. They can be obtained as in the linear case once their finite element counterparts

*M*

^{(FE)}and

*K*

^{(FE)}have been extracted from the finite element software. Specifically, it is found that

**Ψ**

^{(n)}stacked as columns. For finite element models developed in a commercial software, e.g., nastran, abaqus, it is necessary to proceed differently, nonintrusively, to identify the quadratic and cubic stiffness coefficients. The process initially devised in Ref. [23] is adopted here. It relies on imposing on the finite element model static displacement fields that are expressed exactly as linear combinations of the structural modes, i.e., with specified values

*q*

_{j}. The finite element software is then used to determine the force vectors needed on all degrees of freedom to generate these displacements. Projecting them on each of the

*M*basis functions (linear modes, duals, and enrichments) then yields the corresponding modal forces

*F*

_{i}. Imposing that these displacements and modal forces satisfy the ROM governing equation, Eq. (2), leads to the conditions

*j*, i.e.,

Then, Eq. (6) will only involve the coefficients $Kij(1)$, $Kijj(2)$, and $Kijjj(3)$ for each value of *i* independently. Proceeding with three different values of *q*_{j} in Eq. (7) thus generates the 3*M* equations necessary to evaluate all of the coefficients $Kij(1)$, $Kijj(2)$, and $Kijjj(3)$ for a given value of *j*. Repeating this process for every value of *j* = 1, …, *M* then yields all coefficients of these types. Note that the above computations could be reduced to using only two different values of *q*_{j} in Eq. (7) if one relies on the linear stiffness coefficients $Kij(1)$ determined from Eq. (5). However, it is usually recommended to not perform this simplification but rather to use the $Kij(1)$ values of Eq. (5) as a check that the process of Eq. (7) has been well carried out.

*j*and

*l*>

*j*leads to the determination of all ROM coefficients of the form $Kijl(2)$, $Kijll(3)$, and $Kijjl(3)$. Finally, imposing a displacement in the form of

*j*, and

*p*>

*l*>

*j*leads to the necessary equations to evaluate the remaining ROM coefficients $Kijlp(3)$.

The above process requires of the order of *M*^{3}/6 static solutions to complete the identification of the linear, quadratic, and cubic stiffness coefficients. However, an alternative approach has also been developed [16] in which the tangent stiffness matrix of the finite element model is outputted for each displacement field, then projected on the basis and finally equated to the ROM tangent stiffness matrix. In this strategy, only the displacements of the form of Eqs. (7) and (8) are needed and the number of static solutions needed is accordingly reduced to order of *M*^{2}/2.

## 3 A Simple Illustrative Example

Before focusing on the buckled harvester of Fig. 1, it would be desirable to develop a perspective on the effects of an in-plane force on the nonlinear geometric response of a prototypical structure. To this end, it is proposed here to consider the rigid bar linkage model of Fig. 2 which has clear similarities with the harvester but only involves two degrees of freedom, i.e., the horizontal and vertical positions *x* and *y* which are linked by the angle *θ*. The focus of its study here is to assess whether or over what range of displacements its response to the forces *P* and *F* is well captured by the ROM of Eq. (1) with two modes, one transverse and one in-plane.

*N*

_{1}and

*N*

_{2}denote the compressive forces in the bars

*AB*and

*AO*, respectively. The force balance at point

*B*then yields (see Fig. 2 for all variables)

Similarly, at point *A* one has

Note that the linear buckling load of the linkage is *P*_{buck} = *k*_{1}*L*/2.

**Ψ**

^{(1)}representing transverse motions and the second one

**Ψ**

^{(2)}modeling in-plane displacements. Assume next that the structure is symmetric with respect to the horizontal axis as the linkage model of Fig. 2. Then, the response in the transverse direction should change sign with

*F*changing sign but the in-plane displacements should stay the same. Accordingly, only odd powers of

*q*

_{1}should be present in the transverse equation and even ones in the in-plane equation. Then, Eq. (2) becomes in static conditions

Comparing Eqs. (13) and (20), it is clear that the ROM equation form does *not exactly* match the one for the rigid bar linkage. The difference lies in the denominator terms $K22(1)+K2112(3)q12$ for the ROM versus $L2\u2212y2$ for the linkage. However, if *y* is small enough that $L2\u2212y2\u2248L\u2212y2/2L$, its two-term Taylor expansion, then the forms of Eqs. (13) and (20) match. Note further that this approximation holds well for a quite large displacement, for *y*/*L* = 0.6—a very large displacement—the relative difference between $L2\u2212y2$ and *L* − *y*^{2}/2*L* is only 2.5%. On this basis, it is tentatively concluded that the ROM would provide a very close match of the exact linkage behavior for expected smaller displacements under the physical transformation *q*_{1} = *y*, both representing transverse displacements.

*P*given

*y*and

*F*as well as values of

*F*given

*y*and

*P*were generated and used to identify the ROM coefficients $K11(1)$ and $K1111(3)$ and the ratios $K112(2)K211(2)/K22(1)$ and $K2112(3)/K22(1)$ through a best fit minimization of the absolute value between the displacement

*y*of the linkage data and the one predicted as

*q*

_{1}from Eq. (19) for the ensemble of values of

*P =*

*P*and

_{i}*F = F*obtained above. That is, determine the ROM coefficients $K11(1)$, $K1111(3)$, $K112(2)K211(2)/K22(1)$, and $K2112(3)/K22(1)$ to minimize

_{j}*y*(

*P*

_{i},

*F*

_{j}) and

*q*

_{1}(

*P*

_{i},

*F*

_{j}) obtained respectively from Eqs. (13) and (19). Once these coefficients were determined, the vertical displacement predicted from Eq. (19) as

*y*=

*q*

_{1}could be compared to the one determined from Eq. (13). The results of this modeling, shown in Figs. 3 and 4, fully confirm the excellent match of the exact linkage data by the ROM for transverse displacements up to about 60% of the bar length.

*without*in-plane force) have typically assumed that the coefficient $K2112(3)$ can be neglected keeping only as in-plane–transverse coupling in the in-plane equation the term $K211(2)q12$. Performing this simplification in Eq. (19) then leads to the force–displacement relationship

The modeling of the linkage data of Figs. 3 and 4 was repeated but with the model of Eq. (21) in which the coefficients $K11(1)$ and $K^1111(3)$ and the ratio $K112(2)/K22(1)$ were identified by the same optimization as performed for Eq. (19). The results are also presented in Figs. 3 and 4 and they demonstrate that this approximation does not do as well as Eq. (19), it is only valid for *y*/*L* < 0.15 or so which may still be sufficient for practical, structural applications (see the buckled harvester discussion in Sec. 4). The weakness of Eq. (21) in comparison to Eq. (19) lies in the *P* related term. Based on the Taylor series arguments above, it could be expected that this term is correctly captured up to and including third order, (*y*/*L*)^{3}, in Eq. (19) but only up to and including first order, *y*/*L*, in Eq. (21). For the term independent of *P*, third order approximation is always achieved owing to the coefficients $K1111(3)$ and $K^1111(3)$.

It is thus concluded from this analysis that the term $K2112(3)$ does play a significant role in ROMs with nonzero in-plane loading. This result can be justified as well from order of magnitude arguments. Specifically, when *P* = 0, one concludes from Eq. (17) that $q2=O(q12)$ and the term $K2112(3)q12q2$ is of order 4 in *q*_{1}. However, when *P* is present and of order *O*(1) (i.e., of the same order of magnitude as the transverse loading), then its associated deflection *q*_{2} is also of order *O*(1) and $K2112(3)q12q2$ is of order $O(q12)$ as the term $K211(2)q12$. Thus, $K2112(3)$ ought to be retained in the modeling, even though that term is often difficult to identify accurately.

## 4 The Buckled Harvester

### 4.1 Geometry and Nonlinear Response Features.

The geometric and material properties of the beam are the same as those in Ref. [6], see Table 1, which is modeled by an array of 55 × 11 4-node shell elements in nastran (CQUAD4 elements) as shown in Fig. 5. A composite material card was employed to reflect that the beam is composed of a steel shim with two layers of PZT, one on top and one on bottom, that cover it completely. As shown in Fig. 1, the beam is clamped at both ends but its right end is free to move in the in-plane direction. A buckling analysis (SOL 105 in nastran) was first performed and yielded the buckling value of *P* = 10.75 N.

Property | Symbol | Value |
---|---|---|

Length | L_{0} | 55 mm |

Width | b | 11 mm |

Thickness of steel shim | h_{s} | 0.1 mm |

Thickness of piezo layer | h_{p} | 0.08 mm |

Density of steel shim | ρ_{s} | 7850 kg/m^{3} |

Density of piezo layer | ρ_{p} | 4000 kg/m^{3} |

Young's modulus of steel shim | E_{s} | 2.03 × 10^{11} Pa |

Young's modulus of piezo layer | E_{p} | 4.00 × 10^{10} Pa |

Coupling coefficient | D_{31} | −10 pC/N |

Electrical permittivity | $\epsilon 33T$ | 8.854 × 10^{−10} F/m |

Property | Symbol | Value |
---|---|---|

Length | L_{0} | 55 mm |

Width | b | 11 mm |

Thickness of steel shim | h_{s} | 0.1 mm |

Thickness of piezo layer | h_{p} | 0.08 mm |

Density of steel shim | ρ_{s} | 7850 kg/m^{3} |

Density of piezo layer | ρ_{p} | 4000 kg/m^{3} |

Young's modulus of steel shim | E_{s} | 2.03 × 10^{11} Pa |

Young's modulus of piezo layer | E_{p} | 4.00 × 10^{10} Pa |

Coupling coefficient | D_{31} | −10 pC/N |

Electrical permittivity | $\epsilon 33T$ | 8.854 × 10^{−10} F/m |

Given the uncertainty on whether the ROM approach of Sec. 2 would be applicable, a series of static responses without piezoelectric effects were first computed to understand the features of the response and assess the likelihood that they would be well captured by the ROM. Key in this study is the effect of the in-plane force and its interplay with the transverse loading. Accordingly, a first set of computations were performed with a varying in-plane force and a fixed, nonzero uniform pressure (*f* = 42 Pa) as in Fig. 1(b).

A summary of these computations is shown in Fig. 6 from which the very sharp increase in both transverse and in-plane displacements near the buckling force is observed. These figures are consistent with the linkage model results, see Fig. 4. Interestingly, as shown in Fig. 7(a), the shapes of the transverse displacement along the beam are very much independent of the loading conditions and are effectively nearly identical to the buckling mode shown in Fig. 8. The in-plane displacements show more variability but are dominated by the linear displacement field induced by the in-plane force. This situation reflects the small magnitude of the uniform pressure which, without axial force, leads to only very small in-plane displacements. Thus, at low values of the in-plane force, the displacement is nearly linear. For higher forces, the transverse displacements become larger and so do the nonlinearly induced membrane displacements which provide the modulation of the linear displacements shown in Fig. 7(b).

These last observations suggest that, at least phenomenologically, the beam could be analyzed in terms of just 1 mode in the transverse direction which would more specifically be the buckling mode of Fig. 8. This simplification would enable simple plots of displacements—transverse load—and in-plane force, as Figs. 3 and 4, to be made that would clarify the physical behavior of the beam. Figure 7(a) shows however that such a basis would not be sufficient to get *quantitatively* accurate predictions of the transverse response as it changes slightly with varying in-plane force.

To achieve this effort, the transverse displacements were imposed on the beam finite element model in the shape of the buckling mode (to avoid duplication of efforts, the shape of the buckling mode predicted by the ROM was used; it is nearly identical to the finite element one, see Fig. 8 and Sec. 4.3) with amplitudes (denoted here as *q*_{1} since they will be the first ROM generalized coordinates) ranging up to 18 thicknesses of the entire assembly, i.e., 0.26 mm, and with different in-plane forces spanning the domain [0,11.4] N. For each of these combinations, the in-plane displacements were left free and the necessary transverse forces and moments were outputted from nastran. They were then projected on the buckling mode to obtain the corresponding modal force *F*_{1}. The resulting data are shown in Fig. 9.

It should be noted, see Fig. 9(a), that the modal force *F*_{1} is negative for smaller values of *q*_{1} when the in-plane force *P* is larger than the buckling load, i.e., for *P* = 11.4 N. This is expected as the beam would on its own reach its corresponding buckled state which is larger than the imposed displacements. It is thus necessary to hold back the beam.

*F*

_{1}varies linearly with respect to the in-plane force

*P*consistently with the linkage model of Eq. (13) and the condensed ROM equation Eq. (19) even for large values of

*P*and

*q*

_{1}(the value

*q*

_{1}= 8.23 × 10

^{−5}corresponds to a peak transverse displacement of 18 thicknesses). From this observation, it is concluded that

*f*

_{0}(

*q*

_{1}) and

*f*

_{P}(

*q*

_{1}) are shown in Fig. 10 with the latter extracted as

*P*considered generating the extremely thin band of values shown.

The data of Fig. 10 provide an excellent framework to assess the adequacy of the ROM models of Eqs. (19) and (21). Consider first the function *f*_{0}(*q*_{1}) of Fig. 10(a) which represents the behavior of the beam without in-plane force. In Eq. (21), the force–displacement relationship should be a combination of linear and cubic terms only and it is seen in Fig. 10(a) that an approximation with these terms would indeed be very accurate till about 10–13 thicknesses (*q*_{1} values in the range of 4.5 × 10^{−5} to 6 × 10^{−5}). The deviation of the two curves for larger amplitudes could potentially reflect the known difference in the nonlinear modeling used by nastran and by the ROM.

*f*

_{P}(

*q*

_{1}) associated with the in-plane force, see Fig. 10(b). At first, it appears essentially linear with respect to the generalized coordinate

*q*

_{1}as in Eq. (21) but when considering

*f*

_{P}(

*q*

_{1})/

*q*

_{1}, see Fig. 10(c), it is clearly seen that there is a definite deviation from it which increases with increasing value of

*q*

_{1}. Since this deviation remains small through the current range of displacements, i.e., about 2.5% at 18 thicknesses, a Taylor series approximation of the term $1/(K22(1)+K2112(3)q12)$ with respect to

*q*

_{1}seems appropriate and would lead to

As shown in Fig. 10(c), the quadratic fit of *f*_{P}(*q*_{1})/*q*_{1} closely matches the modeling of the nastran data and is thus fully consistent with Eq. (27) with $K112(2)$ and $K2112(3)$ both negative. Interestingly, $K2112(3)$ negative with $K22(1)$ positive is similar to the *L*^{2} − *y*^{2} of the linkage model for which the change of 2.5% in Fig. 10(c) would imply $L/L2\u2212y2=1.025$ and thus *y*/*L* = 0.22. Since the 2.5% change occurs for a peak displacement of 18 thicknesses, one would estimate the length *L* for the harvester to be 82 thicknesses or 21.3 mm. Geometrically, the total length of 55 mm would equal to 2*L* or *L* = 27.5 mm which is surprisingly close to the above 21.3 mm and would seem to give credence to the linkage model as being simple but phenomenologically correct.

*f*

_{0}(

*q*

_{1}) associated with Eq. (21), i.e.,

*decrease*the overall cubic term in this equation as

*q*

_{1}increases which is the contrary of the trend necessary to obtain a better match in Fig. 10(a). Note that this discussion assumes that the coupling mechanism between in-plane and transverse motions is the same for the in-plane motions resulting from the in-plane force and for those induced by the membrane stretching which will be shown not to be the case for the harvester in Sec. 4.3.

### 4.2 Reduced Order Modeling With Piezoelectric Effects and Ground Motion.

Before proceeding with a ROM of the harvester of Fig. 1, it is necessary to extend the formulation of Sec. 2.1 to account for (i) the piezoelectric effects and (ii) the ground motion excitation versus mechanical forcing.

As in prior publications [11,12], the piezoelectric effects will be introduced in the ROM through the thermal analogy. In fact, this analogy must also be used with the nastran finite element model since piezoelectric modeling is not available in this structural software.

*E*

_{11}, due to an applied electrical field in the 3-direction, $E3$, is given as

*d*

_{31}is the piezoelectric coefficient. Assume next that the electrical field is uniform and given as $E3$ =

*V*/

*h*

_{p}where

*V*is the voltage across the PZT element (between the electrodes on the top and bottom surfaces of the beam) and

*h*is the thickness of the PZT. Then, Eq. (29) can be rewritten as

_{p}*T*is the current temperature, and

*T*

_{ref}is the reference temperature.

*T*= 0 leads to the piezoelectric-thermal analogy

_{ref}*τ*

_{n}(

*t*) are the temperature basis functions and corresponding generalized coordinates. Then, assuming that the material properties are independent of temperature (which is the case here since the piezoelectric properties are independent of the applied potential

*V*), it is found that Eq. (2) must be replaced by

*µ*= 1) to consider since the applied potential is the same everywhere on the PZT layers. Moreover, as the PZT covers the entire beam, it is concluded that $\varphi (1)=e$, i.e., the vector whose components are all equal to 1. Finally, from Eq. (32), it is concluded that

*τ*

_{1}=

*V*. Under these conditions, Eq. (34) simplifies as

*th*) of a thermal problem has been switched to (

*p*) to denote the piezoelectric effect. Note that a multithermal mode formulation would be necessary when there are different patches on the beam exhibiting different potentials.

In Eq. (35), $Kij1(p)$ and $Fi1(p)$ are the electromechanical coupling coefficients which need to be evaluated from the finite element model. This identification is again described here using the thermal analogy since nastran does not have a piezoelectric material card. For other software with piezoelectric capabilities, the process is identical but with the application of a voltage.

To obtain $Kij1(p)$, the imposed displacement approach of Eq. (7) is performed twice: once without a temperature to determine $Kij(1)$ and once with a temperature distribution $T\varphi (1)$ (uniform here) for which the linear stiffness terms are $Kij(1)\u2212Kij1(p)T$. Then, combining these values yields $Kij1(p)$.

*i*= 1, …,

*M*, are computed. Finally,

**the vector of components equal to 1 in the direction of the ground motion and 0 otherwise. The absolute displacements of the structure can then be expressed as $W(t)=y(t)u+w(t)$. Since**

*u***is here a rigid body mode for the structure, it induces no strain even in the nonlinear case and thus the equations of motion relating to the relative displacement**

*u***are the same as those for**

*w***without ground motion but with the applied pseudo force**

*W*With the above discussions, the ROM formulation is now complete.

### 4.3 Reduced Order Model Construction and Structural Only Validations.

The construction of the ROM of the buckled harvester of Fig. 1 with the properties of Table 1 started with the selection of the linear (transverse) modes to be included in the basis. Given the expected frequency range of the ground motions, up to 1000 Hz, and their uniformity along the beam, the first two linear transverse symmetric modes were selected with natural frequencies of 270.87 and 1469.51 Hz.

The selection of dual modes to capture the membrane stretching was carried out as in Sec. 2.2 by analyzing the nonlinear response to the loadings of Eq. (4) with the combinations of the linear modes (*j*, *l*) = (1, 1), (1, 2), and (2, 2). For each combination, ten different values *m* = 1, …, 10 of $\alpha j(m)=\alpha l(m)$ were selected so that the beam deformations span from nearly linear to a maximum transverse displacement of six thicknesses. Once these responses were obtained from nastran, their transverse components were first removed by making the displacements orthogonal to the two transverse basis functions. Next, a POD analysis was carried out on the residuals of projection and dominant POD eigenvectors were retained as dual modes. In this process, one dual mode was taken from each combination, and thus three dual modes are constructed.

The final aspect of the basis construction is capturing the displacements induced by the in-plane force. As suggested in Sec. 2.2, it was observed that the three dual modes did not lead to a good modeling of the in-plane linear displacements induced by the force *P*. Accordingly, the basis was further enriched by adding to it the static linear displacement field induced by this in-plane force. This step completed the basis selection which then contained six modes, two transverse and four in-plane.

Finally, all coefficients of the ROM, including the piezoelectric coupling terms, were identified using the tangent stiffness matrix version of the imposed displacement approach of Sec. 2.3 and the above procedure for the piezoelectric-related coefficients.

*α*and

*β*were selected as 84.14/s and 1.8876 10

^{−6}s which leads to damping ratios of 2.6% and 1.3% for the two transverse linear modes.

**is split here into the components associated with the two transverse modes, denoted by $q(t)$, and those corresponding to the in-plane motions (three dual modes and one enrichment) denoted by $q(d)$. Then, without transverse loading but with the in-plane force, Eq. (6) reduces in the linear static case to**

*q**P*.

The lowest eigenvalue was found to be *P* = 10.78 N which matches very well the nastran result of 10.75 N. The agreement of the predicted buckling eigenvector is also very good, see Fig. 8.

As discussed in Sec. 4.1, the transverse displacements of the harvester are most clearly analyzed in terms of the buckling mode which strongly dominates the response, although not entirely, see Fig. 7(a). Accordingly, the transverse basis, i.e., the two linear symmetric transverse modes, was replaced by the two eigenvectors of Eq. (42). In effect, this change is only a rotation of the transverse basis and thus should not change the ROM predictions. It does however bring to the forefront the most important nonlinear stiffness terms in the identification process and thus permits a more accurate estimation of them. The dual modes were not changed as they were based on all combinations of the two linear modes and thus remain relevant for the modeling of the membrane stretching effects.

The identification of the stiffness coefficients was performed again with the new buckling mode-related basis. In this process, it was observed that the condensed cubic stiffness coefficient associated with the buckling mode, $K^1111(3)$, see discussion of Eq. (22) and Ref. [19] for the generalization to multiple in-plane modes, was much less than its “regular” counterpart $K1111(3)$. Specifically, it was found that $K^1111(3)/K1111(3)=5.09\xd710\u22124$ which implies that the stiffening observed in the absence of in-plane motions $(K1111(3))$ is almost completely eliminated by the softening membrane stretching effects. This situation has been observed before in connection with structures that have some free boundary conditions, e.g., cantilevered structures [19] and the freely expanding plate of Ref. [62]. Clearly, the harvester does exhibit this feature as its right end is free to move in-plane. To accurately identify the net stiffening effect on the harvester, the decondensation approach of Ref. [19] was accordingly adopted.

As part of the decondensation, this ROM was “cleaned,” i.e., all coefficients that should be zero by symmetry, see discussion in Sec. 3, were set to that value as opposed to keeping the small but nonzero identified values. Moreover, all cubic stiffness terms coupling the in-plane (denoted by the index *i*) and transverse (denoted by the index *t*) modes and those relevant to the in-plane only modes were zeroed out leaving only the terms $Ktttt(3)$. This elimination includes in particular the terms $Kttii(3)$ where *t* = 1, 2 (mode 1 is the buckling mode from the ROM) and *i* = 3, 4, 5, 6 (mode 3 is the enrichment, mode 4 is the dominant (1,1) dual) some of which were suggested in Sec. 4.1 to play an important role when the deformations become large. This cleaning has been performed in particular for flat structures because the values identified by the algorithms of Refs. [16,23] do not appear consistent with nastran results, e.g., see discussion of Ref. [65].

The assessment of this cleaned ROM was carried out on the uniform static pressure cases of Fig. 6 and its corresponding predictions are also shown in this figure. To complement these results, shown in Fig. 11, are comparisons of the projections of the nastran responses on the buckling mode, enrichment mode, and first dual (modes 1, 3, and 4) with the corresponding generalized coordinates of the ROM. Note in these figures that the in-plane force is scaled by the relevant buckling force, i.e., 10.75 for nastran and 10.78 for the ROM, so that the buckling event occurs for both approaches at the same point. This slight difference and the scaling are the origin of the shift of the enrichment curves of Fig. 11(b) corresponding to the ROM and to nastran. In the absence of scaling, these curves are on top of each other through most of the in-plane force range.

Overall, the ROM performs quite well until the in-plane force becomes very close to the buckling point, i.e., till about 0.98 of that value, but appears to miss part of the physics above it. This conclusion is most easily drawn from the behavior of the enrichment mode, see Fig. 11(b), which the cleaned ROM predicts to continue to increase essentially linearly while the ROM clearly displays a rapid, nonlinear increasing trend. In Figs. 11(a) and 11(c), the change of behavior is reflected by a notable increase in the curvature very near the buckling point which is not present in the cleaned ROM results.

These results and the discussion of Fig. 10(c) suggest that some of the cubic stiffness coefficients coupling transverse and in-plane modes should be included. Modifying the ROM by reintroducing or changing the value of a few stiffnesses has been referred to as “tuning” in prior efforts [66,67] and is warranted when very small such changes have a large impact on the ROM predictions. Its use is intended to palliate the somewhat limited accuracy with which the nonlinear ROM coefficients can be estimated by the strategies of Refs. [16,23]. Such changes were found to be necessary when considering the snap-through behavior of a curved clamped–clamped beam because of the delicate balance of the terms occurring at that condition [66,67]. A similar situation occurs here and is most visible in the condensed ROM as in Eqs. (19) and (21): at buckling the overall linear terms of the buckling mode become zero but the net nonlinear stiffening effect of that mode, measured by $K^1111(3)$, is small. Thus, the presence of any extra term, even if small, is likely to affect the ROM predictions notably.

Following the discussion of Sec. 4.1, the tuning of the cleaned ROM focused on the coefficients $K3114(3)$ and $K4114(3)$ which appeared as the most important ones. The inclusion of the latter one was motivated by the dominance in the harvester response near buckling of mode 4 (first dual), its generalized coordinate being ten times as large as those of modes 5 and 6 and two orders of magnitude or so larger than those of the enrichment mode. The inclusion of $K3114(3)$ was motivated by the fact that only the generalized coordinates of the enrichment (mode 3) are visibly dependent on the in-plane force and by the clear lack of matching shown in Fig. 10(b).

As opposed to the formal, optimization-based tuning approach of Refs. [66,67], the selection of $K3114(3)$ and $K4114(3)$ and the coefficients they are directly related to by symmetry (i.e., $K4113(3)$, $K1134(3)$, and $K1133(3)$) was achieved by a manual trial and error approach which is appropriate here given the small number of such coefficients and the low cost of the ROM computations. Key in that effort was not simply getting a better prediction but actually matching the physical trends of the nastran response. This was achieved by selecting $K3114(3)=3\xd71016$ and $K4114(3)=\u22126\xd71015$ and the corresponding ROM predictions are shown as before in Figs. 6 and 11. Clearly, this tuned ROM performs extremely well, not only are the key physical displacements of Fig. 6 matching their nastran counterparts but also the trends of the generalized coordinates in Fig. 11 are closely matched demonstrating that these changes are genuinely needed to represent the physics of the harvester near buckling.

Note that the values of these coefficients identified by the approach of Ref. [16] were consistently $K3114(3)=1.38\xd71017$ and $K4114(3)=\u22122.3\xd71017$ which are very far from the above values. With these values, the ROM predictions deviate from the nastran ones quite early on and convergence problems are encountered when solving the nonlinear ROM equations near buckling. Similar observations made in the past with other structures, e.g., see Ref. [65], have led to the zeroing out of all $Kttii(3)$ coefficients as a simple, yet suboptimum strategy to improve ROM predictions.

On the basis of the excellent matching observed above, the validations of the tuned ROM focused finally on the dynamic response which is the source of energy harvesting. To this end, it was assumed that the ground acceleration $y\xa8(t)$ can be modeled as a bandlimited white noise in the range [0,1000] Hz and with a standard deviation of 3*g*. The time-step was selected as 5 × 10^{−5} s and the solution was determined for a total of 100,000 steps 20,000 of which were considered part of the transient phase and not used for spectra or standard deviation computations. Then, shown in Fig. 12 are the power spectral densities of the transverse and in-plane displacements at the beam center predicted by nastran and by the ROM for the three in-plane forces *P* = 0 N, 8.4 N, and 10.8 N. Note, as in Fig. 6, that the ROM predictions were obtained for in-plane forces scaled from their nastran counterparts by the factor of 10.78/10.75 to maintain the same ratio of in-plane force over respective buckling force. This scaling is performed in all ensuing plots as well.

When the in-plane force is zero, there is a clear frequency peak in the spectra which is associated with the dominant first transverse mode. When the force is increased to 8.4 N, the peak is shifted to a lower frequency due to the softening effect of the in-plane force. When the force is further increased to 10.8 N, i.e., just above the buckling threshold, the peak in the transverse spectrum almost disappears, suggesting that buckling is indeed occurring. For the first two in-plane force levels, the ROM predictions match the nastran results near perfectly. For the highest in-plane force, i.e., slightly post-buckled conditions, the matching of the transverse spectra is still very good both qualitatively and quantitatively but the in-plane ROM spectrum is slightly shifted from its nastran counterpart. The matching of the standard deviations of the transverse deflection at the center determined from the nastran and ROM time histories is consistent with the transverse spectra of Fig. 12(a), i.e., 0.15, 0.33, 2.18 for the ROM for the values *P* = 0 N, 8.4 N, and 10.8 N versus 0.15, 0.33, and 2.11 for nastran.

From the above results and the static ones of Figs. 6 and 11, it is concluded that the tuned ROM is very good to excellent up to and including the buckling threshold. For comparison, the results obtained with the cleaned ROM are shown in Fig. 13. As expected from the computations of Fig. 6, it is seen that the predictions from this ROM are significantly in error, most notably too small, with respect to those from the full finite element predictions when the in-plane force becomes close to or equal the buckling load.

### 4.4 Energy Harvesting.

*C*

_{p}is the total electric capacity of the piezoelectric patch. This effort is carried out for the ground acceleration of standard deviation

*σ*

_{a}equal to 0.2

*g*, 1

*g*, and 3

*g*and for the set of in-plane forces

*P*considered earlier, i.e.,

*P*= 0 N, 8.4 N, and 10.8 N.

Shown in Fig. 14 are the time histories of the beam center transverse displacement and the output voltage for the three values of *P* and *σ*_{a} = 3*g*. The time history plot of the transverse displacement clearly shows that buckling occurs for *P* = 10.78 N and that the vibration magnitude and the corresponding output voltage are significantly increased for that in-plane force versus the smaller ones owing to frequent excursions from one buckled position to its mirror image. Note that the time variations of the transverse displacement and voltage are very similar because the effect of the nonlinear coupling term $Kij1(p)qi(t)qj(t)$ is very small owing to a small value of $Kij1(p)$.

To clarify the dependence of the voltage output on the in-plane force, the above computations were repeated for a series of values of *P* in [0, 10.8] N and the standard deviation of the voltage was determined for each case. The results of these computations are shown in Fig. 15. It can be seen that when the in-plane force is small and well below the buckling value, the output voltage is small. However, when the in-plane force is increased, the output voltage becomes larger, increasing very rapidly when the in-plane force becomes close to the buckling value as previously observed [7]. For comparison, the corresponding results obtained with the cleaned ROM are also shown in Fig. 15. They parallel the discussion of Figs. 12 and 13, i.e., the power output is well computed until the in-plane force approaches the buckling load. In the neighborhood of that value, the power output is significantly underestimated by the cleaned ROM.

The strength of reduced order modeling is really highlighted by the computations of Fig. 15 as each data point takes 14,763 s for the full finite element computations but only 46 s for the ROM ones, i.e., 320 times faster, on a DELL Precision T7500 workstation with 48 GB of RAM. Moreover, the ROM provides a compact, easily codable structural model that can be used in the design of the electrical circuitry to optimize the power output, a task that could not be carried out with a finite element model in a commercial software such as nastran.

## 5 Summary

The focus of this investigation has been on assessing the applicability of nonintrusive ROMs for the nonlinear geometric response of structures subjected to in-plane forces near buckling such as an energy harvester recently proposed. The overarching conclusion is that these ROMs are indeed applicable and can provide excellent static and dynamic predictions of the response but the identification of their nonlinear stiffness coefficients must be done carefully to capture the delicate balance of the small linear terms, nonlinear terms, and transverse loading that occurs very near buckling. Specifically, it was proposed that the basis explicitly include the buckling mode and the linear response of the structure to the in-plane force (referred to as enrichment here) to reduce the number of dominant terms in the model. Next, as appropriate, a careful identification of the net stiffening effects of the structure, after the inclusion of the softening membrane effects, may need to be carried out following the approach referred to as decondensation. Finally, it was demonstrated that the cubic stiffness terms coupling the buckling mode, the enrichment, and the dominant dual mode play a significant role very near and at buckling. This finding was supported by the analysis of both the buckled harvester model but also a simple but phenomenologically appropriate linkage model. The positive importance of these cubic stiffness terms has not been demonstrated before and requires that they be appropriately tuned because their identification using existing strategies is not reliable.

The developed ROMs provide extremely efficient computational tools for the analysis of the piezoelectric energy harvesters as they dramatically reduce the computational effort necessary to produce the time histories of the dynamic response on which they are based. Moreover, the ROM is a compact, easily programmable structural model that can be used as a black box for fully coupled electromechanical analyses at the contrary of finite element models in commercial structural software.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## References

*Nonlinear Dynamics*, Vol. 1,