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 [110].

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. [1446] 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 [1427], extensions of the methodology to (i) structures composed of substructures [2831], (ii) structures with defects [3236], (iii) to address design issues [26,37], to include (iv) piezoelectric effects [12,38], (v) viscoelastic damping [39], and (vi) the effects of temperature [4046]. Applications/validations have been described in these papers and others, e.g., Refs. [4762], 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.

Fig. 1
Buckled beam piezoelectric energy harvester from Ref. [6] under (a) ground excitation y(t) and (b) pressure f
Fig. 1
Buckled beam piezoelectric energy harvester from Ref. [6] under (a) ground excitation y(t) and (b) pressure f
Close modal

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

The reduced order models considered here represent the vector of displacements of the degrees of freedom of a finite element model in a “modal expansion” form, i.e., as
w(t)=n=1Mqn(t)Ψ(n)
(1)

In this equation, the vectors Ψ(n) are the basis functions and qn(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.

The set of (nonlinear) ordinary differential equations governing the evolution of the generalized coordinates qn(t) is usually assumed to be nonlinear in stiffness only with that component being a cubic polynomial of these coordinates, i.e.,
Mijq¨j+Dijq˙j+Kij(1)qj+Kijl(2)qjql+Kijlp(3)qjqlqp=Fi,i=1,2,,M
(2)
where summation over repeated indices is assumed here and in the sequel so that 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 Mij, Dij, 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 Fi denotes the ith 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. [1462]. 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,3946,53,5559] 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)Ψ(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.

On that basis, the dual mode construction proceeds as follows. First, the loadings
FFE(m)=αj(m)KFE(1)Ψ(j)(nosummationonj)
(3)
and
FFE(m)=αj(m)KFE(1)Ψ(j)+αl(m)KFE(1)Ψ(l)(nosummationonj,l)
(4)
are applied statically to the structure for a series of values of αj(m) and αl(m) such that the corresponding deformations range from nearly linear to significantly nonlinear. Then, the displacement fields are projected along the linear modes kept in the basis and, finally, the residual of projections are processed through a proper orthogonal decomposition (POD) the eigenvectors of which form the dual modes, see Refs. [13,15,16,21,22,62] for details.

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.

While integral expressions for all coefficients of Eq. (2) are given in Ref. [64], it is more convenient to identify these coefficients directly from the finite element model (see Ref. [63] for an application of the formulas of Ref. [64]). In this regard, note that the terms Kijl(2)qjql and Kilj(2)qlqj for jl 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 ≤ jM, jlM, and lpM, 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
M=ΨTM(FE)ΨandK(1)=ΨTK(FE)Ψ
(5)
where Ψ denotes the matrix formed by the basis functions Ψ(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 qj. 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 Fi. Imposing that these displacements and modal forces satisfy the ROM governing equation, Eq. (2), leads to the conditions
Kij(1)qj+Kijl(2)qjql+Kijlp(3)qjqlqp=Fi,i=1,,M
(6)
to be satisfied by the ROM coefficients Kij(1), Kijl(2), and Kijlp(3).
Selecting appropriately the imposed displacements allows for a simple determination of each of these coefficients. Specifically, impose first displacements proportional to a single mode, say mode j, i.e.,
w1=qjΨ(j)(nosummationonj)
(7)

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 qj in Eq. (7) thus generates the 3M 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 qj 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.

In the second step of the approach, the imposed displacements are linear combinations of two modes, e.g.,
w2=qjΨ(j)+qlΨ(l)(nosummationonjorl)
(8)
Imposing three such combinations for every 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
w3=qjΨ(j)+qlΨ(l)+qpΨ(p)(nosummationonjorlorp)
(9)
for every 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 M3/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 M2/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.

Fig. 2
Rigid bar linkage model. NLk1, k3 denotes a nonlinear spring with linear (k1) and cubic (k3) stiffnesses.
Fig. 2
Rigid bar linkage model. NLk1, k3 denotes a nonlinear spring with linear (k1) and cubic (k3) stiffnesses.
Close modal
In deriving the equations of motion of the linkage model, let N1 and N2 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)
N1cosθ=P2k2x
(10)

Similarly, at point A one has

horizontally:
(N1N2)cosθ=k¯2x
(11)
vertically:
F=Fs(N1+N2)sinθ
(12)
where Fs = k1y + k3y3 is the force in the nonlinear spring acting vertically. Taking into account that x = L (1 − cos θ) and y = L sin θ, it is found from Eqs (10)(12) that
F=k1y+k3y3+(4k2+k¯2)(LL2y21)yP2yL2y2
(13)

Note that the linear buckling load of the linkage is Pbuck = k1L/2.

To connect the above analysis to a reduced order model of an actual structure (such as the beam of Fig. 1), assume that only two modes are selected to represent its motions as in Eqs. (1) and (2) with the first mode Ψ(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 q1 should be present in the transverse equation and even ones in the in-plane equation. Then, Eq. (2) becomes in static conditions
K11(1)q1+K112(2)q1q2+K1111(3)q13+K1122(3)q1q22=F
(14)
K22(1)q2+K211(2)q12+K222(2)q22+K2112(3)q12q2+K2222(3)q23=P
(15)
To reduce further Eqs. (14) and (15), note that the in-plane force P is involved linearly in Eq. (13) which can be matched by Eqs. (14) and (15) when only linear terms in q2 are present in these equations, i.e.,
K11(1)q1+K112(2)q1q2+K1111(3)q13=F
(16)
K22(1)q2+K211(2)q12+K2112(3)q12q2=P
(17)
Then, solving Eq. (17) for q2 yields
q2=PK211(2)q12K22(1)+K2112(3)q12
(18)
Inserting then this expression in Eq. (16) leads finally to
[K11(1)q1+K1111(3)q13K112(2)K211(2)q13K22(1)+K2112(3)q12]+[K112(2)q1K22(1)+K2112(3)q12]P=F
(19)
or equivalently
[K11(1)q1+K1111(3)q13+K112(2)K211(2)K2112(3)(K22(1)K22(1)+K2112(3)q121)q1]+[K112(2)q1K22(1)+K2112(3)q12]P=F
(20)

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 L2y2 for the linkage. However, if y is small enough that L2y2Ly2/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 L2y2 and Ly2/2L 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 q1 = y, both representing transverse displacements.

To confirm the above discussion, data from the linkage model, i.e., values of 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 q1 from Eq. (19) for the ensemble of values of P =Pi and F = Fj 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
Emod=i,j|y(Pi,Fj)q1(Pi,Fj)|
with y(Pi, Fj) and q1(Pi, Fj) obtained respectively from Eqs. (13) and (19). Once these coefficients were determined, the vertical displacement predicted from Eq. (19) as y = q1 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.
Fig. 3
Transverse force (F) versus transverse displacement (y) for different in-plane forces (P). Linkage data and its approximations by full and simplified ROMs: L = 1, k1 = 1, k3 = 3, and k2=k¯2=0.
Fig. 3
Transverse force (F) versus transverse displacement (y) for different in-plane forces (P). Linkage data and its approximations by full and simplified ROMs: L = 1, k1 = 1, k3 = 3, and k2=k¯2=0.
Close modal
Fig. 4
In-plane force (P) versus transverse displacement (y) for different transverse forces (F). Linkage data and its approximations by full and simplified ROMs: L = 1, k1 = 1, k3 = 3, and k2=k¯2=0.
Fig. 4
In-plane force (P) versus transverse displacement (y) for different transverse forces (F). Linkage data and its approximations by full and simplified ROMs: L = 1, k1 = 1, k3 = 3, and k2=k¯2=0.
Close modal
The model of Eq. (19) is often referred to as a “condensed” model as the in-plane displacements have been eliminated leaving only the transverse ones to describe the force-response behavior. Past applications of condensed ROMs (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
[K11(1)q1+K^1111(3)q13]+K112(2)K22(1)Pq1=F
(21)
where
K^1111(3)=K1111(3)K112(2)K211(2)K22(1)
(22)
Note that Eq. (21) correctly recovers the buckling load
Pbuck=K11(1)K22(1)(K112(2))
(23)
where K112(2) is negative.

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 q1. 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 q2 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.

Fig. 5
Finite element mesh of the beam harvester
Fig. 5
Finite element mesh of the beam harvester
Close modal
Table 1

Geometric and material properties of the beam

PropertySymbolValue
LengthL055 mm
Widthb11 mm
Thickness of steel shimhs0.1 mm
Thickness of piezo layerhp0.08 mm
Density of steel shimρs7850 kg/m3
Density of piezo layerρp4000 kg/m3
Young's modulus of steel shimEs2.03 × 1011 Pa
Young's modulus of piezo layerEp4.00 × 1010 Pa
Coupling coefficientD31−10 pC/N
Electrical permittivityε33T8.854 × 10−10 F/m
PropertySymbolValue
LengthL055 mm
Widthb11 mm
Thickness of steel shimhs0.1 mm
Thickness of piezo layerhp0.08 mm
Density of steel shimρs7850 kg/m3
Density of piezo layerρp4000 kg/m3
Young's modulus of steel shimEs2.03 × 1011 Pa
Young's modulus of piezo layerEp4.00 × 1010 Pa
Coupling coefficientD31−10 pC/N
Electrical permittivityε33T8.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).

Fig. 6
Static displacements induced by the combinations of transverse uniform pressures f = 2.8, 14, and 42 Pa and in-plane forces, from nastran and the ROM (see Sec. 4.3): (a) transverse displacement at the center point and (b) in-plane displacement of the right endpoint
Fig. 6
Static displacements induced by the combinations of transverse uniform pressures f = 2.8, 14, and 42 Pa and in-plane forces, from nastran and the ROM (see Sec. 4.3): (a) transverse displacement at the center point and (b) in-plane displacement of the right endpoint
Close modal
Fig. 7
Displacements along the beam for a uniform pressure f = 42 Pa and various in-plane forces predicted by nastran: (a) transverse and (b) in-plane displacements. Scaled by the peak value.
Fig. 7
Displacements along the beam for a uniform pressure f = 42 Pa and various in-plane forces predicted by nastran: (a) transverse and (b) in-plane displacements. Scaled by the peak value.
Close modal
Fig. 8
Buckling eigenvector, predicted by both nastran and the ROM
Fig. 8
Buckling eigenvector, predicted by both nastran and the ROM
Close modal

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 q1 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 F1. The resulting data are shown in Fig. 9.

Fig. 9
Modal force along the buckling mode versus in-plane force P for imposed transverse displacements along the buckling mode at various levels q1. The value q1 = 8.23 × 10−5 corresponds to a peak transverse displacement of 18 thicknesses.
Fig. 9
Modal force along the buckling mode versus in-plane force P for imposed transverse displacements along the buckling mode at various levels q1. The value q1 = 8.23 × 10−5 corresponds to a peak transverse displacement of 18 thicknesses.
Close modal

It should be noted, see Fig. 9(a), that the modal force F1 is negative for smaller values of q1 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.

Quite interestingly, it is seen from Fig. 9(b) that the necessary modal force F1 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 q1 (the value q1 = 8.23 × 10−5 corresponds to a peak transverse displacement of 18 thicknesses). From this observation, it is concluded that
F1(q1,P)=f0(q1)+PfP(q1)
(24)
The functions f0(q1) and fP(q1) are shown in Fig. 10 with the latter extracted as
fP(q1)=F1(q1,P)f0(q1)P
(25)
for each value of P considered generating the extremely thin band of values shown.
Fig. 10
Functions (a) f0(q1), (b) fP(q1), and (c) fP(q1)/q1 for the beam harvester and their approximations
Fig. 10
Functions (a) f0(q1), (b) fP(q1), and (c) fP(q1)/q1 for the beam harvester and their approximations
Close modal

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 f0(q1) 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 (q1 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.

Of particular interest here is the analysis of the force component fP(q1) associated with the in-plane force, see Fig. 10(b). At first, it appears essentially linear with respect to the generalized coordinate q1 as in Eq. (21) but when considering fP(q1)/q1, see Fig. 10(c), it is clearly seen that there is a definite deviation from it which increases with increasing value of q1. 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 q1 seems appropriate and would lead to
fP(q1)=K112(2)q1K22(1)+K2112(3)q12K112(2)K22(1)[1K2112(3)K22(1)q12]q1
(26)
so that
fP(q1)q1K112(2)K22(1)[1K2112(3)K22(1)q12]
(27)

As shown in Fig. 10(c), the quadratic fit of fP(q1)/q1 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 L2y2 of the linkage model for which the change of 2.5% in Fig. 10(c) would imply L/L2y2=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 2L 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.

The above findings provide the first supporting evidence for the need to include in the ROM cubic terms of the form of K2112(3) that couple transverse and in-plane coordinates and of the superiority of the ROM model of Eq. (19) over the one of Eq. (21). Having established the presence of these terms, it may be wondered whether their inclusion in the component f0(q1) associated with Eq. (21), i.e.,
f0(q1)=K11(1)q1+K1111(3)q13K112(2)K211(2)q13K22(1)+K2112(3)q12
(28)
would lead to a better match in Fig. 10(a) than obtained with the cubic fit. In this regard, note that K211(2)=K112(2)/2 given the symmetry of coefficients induced by the existence of a potential. Then, the introduction of the term K2112(3) in Eq. (28) would decrease the overall cubic term in this equation as q1 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.

This analogy is based on the similarity of the strains induced by piezoelectric and thermal effects. For the current beam model, the piezoelectric material has direction 1 defined along the in-plane direction, from one support to the other, and direction 3 defined as transverse to the beam. Then, the normal Green strain in the 1-direction induced by the piezoelectric effect, E11, due to an applied electrical field in the 3-direction, E3, is given as
E11=d31E3
(29)
where d31 is the piezoelectric coefficient. Assume next that the electrical field is uniform and given as E3 = V/hp where V is the voltage across the PZT element (between the electrodes on the top and bottom surfaces of the beam) and hp is the thickness of the PZT. Then, Eq. (29) can be rewritten as
E11=d31V/hp
(30)
On the other hand, the thermal strain induced by a temperature change is given as
E11=α(TTref)
(31)
where α is the coefficient of thermal expansion, T is the current temperature, and Tref is the reference temperature.
Comparing Eqs. (30) and (31) and setting Tref = 0 leads to the piezoelectric-thermal analogy
T=Vwhenα=d31/hp
(32)
that is, the piezoelectric effect can be equivalently modeled by a temperature change when the thermal expansion coefficient is set according to Eq. (32). This analogy effectively replaces the structure with PZT by the same structure without piezoelectric effect but with an applied temperature.
Thus, instead of developing a new ROM methodology to include piezoelectric effects (see Ref. [38] for such an approach), one can rely on the ROM work developed for heated structures [4046]. In these investigations, it has been assumed that the temperature can also be expressed in a modal form, i.e., as
T(t)=n=1μτn(t)ϕ(n)
(33)
where the vectors ϕ(n) and τ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
Mijq¨j+Dijq˙j+Kij(1)qjKijl(th)qjτl+Kijl(2)qjql+Kijlp(3)qjqlqp=Fi+Fil(th)τl
(34)
In the present context, there is a single thermal mode (µ = 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 ϕ(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
Mijq¨j+Dijq˙j+(Kij(1)Kij1(p)V)qj+Kijl(2)qjql+Kijlp(3)qjqlqp=Fi+Fi1(p)V
(35)
where the superscript (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ϕ(1) (uniform here) for which the linear stiffness terms are Kij(1)Kij1(p)T. Then, combining these values yields Kij1(p).

To identify the coefficients Fi1(p), the temperature distribution Tϕ(1) is applied again to the nastran mesh with all degrees of freedom blocked. The corresponding reaction forces of a linear static solution are then outputted and from them the corresponding modal forces Fi, i = 1, …, M, are computed. Finally,
Fi1(p)=Fi/T
(36)
With regard to ROM formulation, it remains only to discuss the handling of the ground excitation versus applied forces as shown in Eqs. (2) and (35). The standard discussion of this issue in the linear case applies here as well. Specifically, denote by u 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 w are the same as those for W without ground motion but with the applied pseudo force
Fgr(t)=M(FE)uy¨(t)
(37)

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 αj(m)=α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.

In anticipation of dynamic computations, it is necessary to also introduce damping in the harvester. Specifically, the Rayleigh damping model was adopted here because it can be introduced conveniently and consistently in the ROM and in the finite element model. The elements of the corresponding damping matrices are
Dij(FE)=αMij(FE)+βKij(FE)andDij=αMij+βKij(1)
(38)
in the finite element model and in the ROM, respectively. In these equations, the 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.
Given its intended application, the first validation of the ROM focused on the prediction of the buckling load and the associated eigenvector. The value obtained from nastran, using a SOL 105 command, was 10.75 N. The determination of its counterpart for the ROM is carried out from Eq. (6) in which nonlinear terms are neglected. For clarity of presentation, the vector of generalized coordinates q 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
Kij(tt)qj(t)+Kijl(ttd)qj(t)ql(d)=0,i=1,2,summationoverj=1,2andoverl=1,,4
(39)
and
Ksl(dd)ql(d)=Fs(d)=F^s(d)Ps=1,,4,summationoverl=1,,4
(40)
where Kij(tt) and Ksl(dd) are the elements of the transverse–transverse and in-plane–in-plane blocks of the linear ROM stiffness matrix. Moreover, Kijl(ttd) is an element of the part of the third-order tensor Kijl(2) that involves the transverse–transverse–in-plane components. Finally, in Eq. (40)Fs(d)=F^s(d)P are the modal forces associated with the in-plane force and thus are proportional to P.
Solving Eq. (40) for the vector q(d) yields
q(d)=[K(dd)]1F^(d)P
(41)
Introducing this expression in Eq. (39) finally yields the eigenvalue problem
Kij(tt)qj(t)={Kijl(ttd)[K(dd)]ls1F^s(d)}Pqj(t)i=1,2,sumoverj=1,2;l,s=1,,4
(42)

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×104 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.

Fig. 11
Components of static response along (a) buckling mode, (b) enrichment mode, and (c) first dual, predicted by nastran and the ROM for a uniform pressure f = 42 Pa
Fig. 11
Components of static response along (a) buckling mode, (b) enrichment mode, and (c) first dual, predicted by nastran and the ROM for a uniform pressure f = 42 Pa
Close modal

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×1016 and K4114(3)=6×1015 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×1017 and K4114(3)=2.3×1017 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¨(t) can be modeled as a bandlimited white noise in the range [0,1000] Hz and with a standard deviation of 3g. 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.

Fig. 12
Power spectral densities of the displacements of the beam center, tuned ROM, and nastran predictions. Standard deviation of transverse loading of 3g and various in-plane forces: (a) transverse and (b) in-plane.
Fig. 12
Power spectral densities of the displacements of the beam center, tuned ROM, and nastran predictions. Standard deviation of transverse loading of 3g and various in-plane forces: (a) transverse and (b) in-plane.
Close modal

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.

Fig. 13
Power spectral densities of the displacements of the beam center, cleaned ROM, and nastran predictions. Standard deviation of transverse loading of 3g and various in-plane forces: (a) transverse and (b) in-plane.
Fig. 13
Power spectral densities of the displacements of the beam center, cleaned ROM, and nastran predictions. Standard deviation of transverse loading of 3g and various in-plane forces: (a) transverse and (b) in-plane.
Close modal

4.4 Energy Harvesting.

The time histories of the structural generalized coordinates obtained from the constructed ROM are further used to predict the voltage output of the piezoelectric layers. This is done by assuming an open-loop circuit and ignoring the piezoelectric feedback effect on the structural vibration. Specifically, the formula derived in Ref. [12] is revised for the current uniform electric field, and the voltage output due to structural vibration is then
V(t)=[12Kij1(p)qi(t)qj(t)+Fi1(p)qi(t)]/Cp
(43)
where Cp 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.2g, 1g, and 3g 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 = 3g. 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).

Fig. 14
Time histories of (a) beam center structural displacement and (b) output voltage predicted by the ROM. Ground acceleration of standard deviation σa = 3g.
Fig. 14
Time histories of (a) beam center structural displacement and (b) output voltage predicted by the ROM. Ground acceleration of standard deviation σa = 3g.
Close modal

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.

Fig. 15
Standard deviation (STD) of voltage output induced by the ground motions of standard deviations σa as a function of the in-plane force, ROM computations
Fig. 15
Standard deviation (STD) of voltage output induced by the ground motions of standard deviations σa as a function of the in-plane force, ROM computations
Close modal

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

1.
Daqaq
,
M. F.
,
Masana
,
R.
,
Erturk
,
A.
, and
Quinn
,
D. D.
,
2014
, “
On the Role of Nonlinearities in Vibratory Energy Harvesting: A Critical Review and Discussion
,”
Appl. Mech. Rev.
,
66
(
4
), p.
040801
.
2.
Hu
,
N.
, and
Burgueno
,
R.
,
2015
, “
Buckling-Induced Smart Applications: Recent Advances and Trends
,”
Smart Mater. Struct.
,
24
(
6
), p.
063001
.
3.
Arrieta
,
A. F.
,
Hagedorn
,
P.
,
Erturk
,
A.
, and
Inman
,
D. J.
,
2010
, “
A Piezoelectric Bistable Plate for Nonlinear Broadband Energy Harvesting
,”
Appl. Phys. Lett.
,
97
(
10
), p.
104102
.
4.
Chen
,
L.
, and
Jiang
,
W.
,
2014
, “
Snap-Through Piezoelectric Energy Harvesting
,”
J. Sound Vib.
,
333
(
18
), pp.
4314
4325
.
5.
Wu
,
N.
,
He
,
Y.
, and
Fu
,
J.
,
2021
, “
Bistable Energy Harvester Using Easy Snap-Through Performance to Increase Output Power
,”
Energy
,
226
, p.
120414
.
6.
Cottone
,
F.
,
Gammaitoni
,
L.
,
Vocca
,
H.
,
Ferrari
,
M.
, and
Ferrari
,
V.
,
2012
, “
Piezoelectric Buckled Beams for Random Vibration Energy Harvesting
,”
Smart Mater. Struct.
,
21
(
3
), p.
035021
.
7.
Vocca
,
H.
,
Cottone
,
F.
,
Neri
,
I.
, and
Gammaitoni
,
L.
,
2013
, “
A Comparison Between Nonlinear Cantilever and Buckled Beam for Energy Harvesting
,”
Eur. Phys. J. Spec. Top.
,
222
(
7
), pp.
1699
1705
.
8.
Masana
,
R.
, and
Daqaq
,
M. F.
,
2011
, “
Electromechanical Modeling and Nonlinear Analysis of Axially Loaded Energy Harvesters
,”
ASME J. Vib. Acoust.
,
133
(
1
), p.
011007
.
9.
Pennisi
,
G.
,
Mann
,
B. P.
,
Naclerio
,
N.
,
Stephan
,
C.
, and
Michon
,
G.
,
2018
, “
Design and Experimental Study of a Nonlinear Energy Sink Coupled to an Electromagnetic Energy Harvester
,”
J. Sound Vib.
,
437
, pp.
340
357
.
10.
Fang
,
S.
,
Chen
,
K.
,
Xing
,
J.
,
Zhou
,
S.
, and
Liao
,
W.-H.
,
2021
, “
Tuned Bistable Nonlinear Energy Sink for Simultaneously Improved Vibration Suppression and Energy Harvesting
,”
Int. J. Mech. Sci.
,
212
, p.
106838
.
11.
Wang
,
X. Q.
,
Liao
,
Y.
, and
Mignolet
,
M. P.
,
2021
, “
Uncertainty Analysis of Piezoelectric Vibration Energy Harvesters Using a Finite Element Level-Based Maximum Entropy Approach
,”
ASCE-ASME J. Risk Uncertainty Eng. Syst. Part B: Mech. Eng.
,
7
(
1
), p.
010906
.
12.
Vyas
,
V.
,
Wang
,
X. Q.
,
Jain
,
A.
, and
Mignolet
,
M. P.
,
2015
, “
Nonlinear Geometric Reduced Order Model for the Response of a Beam With a Piezoelectric Actuator
,”
Proceedings of the 56th SDM Conference (AIAA SciTech)
,
Kissimmee, FL
,
Jan. 5–9
, AIAA Paper No. AIAA-2015-0692.
13.
Mignolet
,
M. P.
,
Przekop
,
A.
,
Rizzi
,
S. A.
, and
Spottswood
,
S. M.
,
2013
, “
A Review of Indirect/Non-intrusive Reduced Order Modeling of Nonlinear Geometric Structures
,”
J. Sound Vib.
,
332
(
10
), pp.
2437
2460
.
14.
Hollkamp
,
J. J.
,
Gordon
,
R. W.
, and
Spottswood
,
S. M.
,
2005
, “
Nonlinear Modal Models for Sonic Fatigue Response Prediction: A Comparison of Methods
,”
J. Sound Vib.
,
284
(
3–5
), pp.
1145
1163
.
15.
Kim
,
K.
,
Radu
,
A.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2013
, “
Nonlinear Reduced Order Modeling of Isotropic and Functionally Graded Plates
,”
Int. J. Non-Linear Mech.
,
49
, pp.
100
110
.
16.
Perez
,
R. A.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2014
, “
Non-Intrusive Structural Dynamic Reduced Order Modeling for Large Deformations: Enhancements for Complex Structures
,”
ASME J. Comput. Nonlinear Dyn.
,
9
(
3
), p.
031008
.
17.
Przekop
,
A.
,
Guo
,
X.
, and
Rizzi
,
S. A.
,
2012
, “
Alternative Modal Basis Selection Procedures for Nonlinear Random Response Simulation
,”
J. Sound Vib.
,
331
(
17
), pp.
4005
4024
.
18.
Givois
,
A.
,
Grolet
,
A.
,
Thomas
,
O.
, and
Deü
,
J.-F.
,
2019
, “
On the Frequency Response Computation of Geometrically Nonlinear Flat Structures Using Reduced-Order Finite Element Models
,”
Nonlinear Dyn.
,
97
(
2
), pp.
1747
1781
.
19.
Wang
,
X. Q.
,
Khanna
,
V.
,
Kim
,
K.
, and
Mignolet
,
M. P.
,
2021
, “
Nonlinear Reduced Order Modeling of Flat Cantilevered Structures: Challenges and Remedies
,”
ASCE J. Aerosp. Eng.
,
34
(
6
), p.
04021085
.
20.
Mahdiabadi
,
M. K.
,
Tiso
,
P.
,
Brandt
,
A.
, and
Rixen
,
D. J.
,
2021
, “
A Non-intrusive Model-Order Reduction of Geometrically Nonlinear Structural Dynamics Using Modal Derivatives
,”
Mech. Syst. Signal Process
,
147
, p.
107126
.
21.
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2021
, “
Discussion on “A Non-intrusive Model-Order Reduction of Geometrically Nonlinear Structural Dynamics Using Modal Derivatives
,”
Mech. Syst. Signal Process
,
159
, p.
107638
.
22.
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2020
, “
Toward a Systematic Construction of the Basis for Nonlinear Geometric Reduced Order Models
,”
Proceedings of the 11th International Conference on Structural Dynamics (EURODYN 2020)
,
Athens, Greece
, in virtual,
Nov. 23–25
.
23.
Muravyov
,
A. A.
, and
Rizzi
,
S. A.
,
2003
, “
Determination of Nonlinear Stiffness With Application to Random Vibration of Geometrically Nonlinear Structures
,”
Comput. Struct.
,
81
(
15
), pp.
1513
1523
.
24.
Vizzaccaro
,
A.
,
Givois
,
A.
,
Longobardi
,
P.
,
Shen
,
Y.
,
Deü
,
J.-F.
,
Salles
,
L.
,
Touzé
,
C.
, and
Thomas
,
O.
,
2020
, “
Non-intrusive Reduced Order Modelling for the Dynamics of Geometrically Nonlinear Flat Structures Using Three-Dimensional Finite Elements
,”
Comput. Mech.
,
66
(
6
), pp.
1293
1319
.
25.
Shen
,
Y.
,
Vizzaccaro
,
A.
,
Kesmia
,
N.
,
Yu
,
T.
,
Salles
,
L.
,
Thomas
,
O.
, and
Touzé
,
C.
,
2021
, “
Comparison of Reduction Methods for Finite Element Geometrically Nonlinear Beam Structures
,”
Vibration
,
4
(
1
), pp.
175
204
.
26.
Kim
,
E.
, and
Cho
,
M.
,
2017
, “
Equivalent Model Construction for a Non-linear Dynamic System Based on an Element-Wise Stiffness Evaluation Procedure and Reduced Analysis of the Equivalent System
,”
Comput. Mech.
,
60
(
5
), pp.
709
724
.
27.
Shen
,
Y.
,
Béreux
,
N.
,
Frangi
,
A.
, and
Touzé
,
C.
,
2021
, “
Reduced Order Models for Geometrically Nonlinear Structures: Assessment of Implicit Condensation in Comparison With Invariant Manifold Approach
,”
Eur. J. Mech. A/Solids
,
86
, p.
104165
.
28.
Wang
,
Y.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2018
, “
Component-Centric Reduced Order Modeling for the Prediction of the Nonlinear Geometric Response of a Part of a Stiffened Structure
,”
ASME J. Comput. Nonlinear Dyn.
,
13
(
13
), p.
121006
.
29.
Kuether
,
R. J.
,
Allen
,
M. S.
, and
Hollkamp
,
J. J.
,
2016
, “
Modal Substructuring of Geometrically Nonlinear Finite-Element Models
,”
AIAA J.
,
54
(
2
), pp.
691
702
.
30.
Kuether
,
R. J.
,
Allen
,
M. S.
, and
Hollkamp
,
J. J.
,
2017
, “
Modal Substructuring of Geometrically Nonlinear Finite Element Models With Interface Reduction
,”
AIAA J.
,
55
(
5
), pp.
1695
1706
.
31.
Mahdiabadi
,
M. K.
,
Bartl
,
A.
,
Xu
,
D.
,
Tiso
,
P.
, and
Rixen
,
D. J.
,
2019
, “
An Augmented Free-Interface-Based Modal Substructuring for Nonlinear Structural Dynamics Including Interface Reduction
,”
J. Sound Vib.
,
462
, p.
114915
.
32.
Perez
,
R. A.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2014
, “
Prediction of Displacement and Stress Fields of a Notched Panel With Geometric Nonlinearity by Reduced Order Modeling
,”
J. Sound Vib.
,
333
(
24
), pp.
6572
6589
.
33.
Wang
,
X. Q.
,
Phlipot
,
G. P.
,
Perez
,
R. A.
, and
Mignolet
,
M. P.
,
2018
, “
Locally Enhanced Reduced Order Modeling for the Nonlinear Geometric Response of Structures with Defects
,”
Int. J. Non-Linear Mech.
,
101
, pp.
1
7
.
34.
O’Hara
,
P. J.
, and
Hollkamp
,
J. J.
,
2014
, “
Modeling Vibratory Damage With Reduced-Order Models and the Generalized Finite Element Method
,”
J. Sound Vib.
,
333
(
24
), pp.
6637
6650
.
35.
Wang
,
X. Q.
,
O'Hara
,
P. J.
,
Mignolet
,
M. P.
, and
Hollkamp
,
J. J.
,
2019
, “
Reduced Order Modeling With Local Enrichment for the Nonlinear Geometric Response of a Cracked Panel
,”
AIAA J.
,
57
(
1
), pp.
421
436
.
36.
Marconi
,
J.
,
Tiso
,
P.
, and
Braghin
,
F.
,
2020
, “
A Nonlinear Reduced Order Model With Parametrized Shape Defects
,”
Comput. Methods Appl. Mech. Eng.
,
360
, p.
112785
.
37.
Hollkamp
,
J. J.
,
Perez
,
R. A.
, and
Spottswood
,
S. M.
,
2017
, “
Design Sensitivities of Components Using Nonlinear Reduced-Order Models and Complex Variables
,”
Proceedings of the International Modal Analysis Conference, IMAC XXXV
,
Garden Grove, CA
,
Jan. 30–Feb. 2
, also in: G. Kerschen, ed., Nonlinear Dynamics, Vol. 1,
Springer
.
38.
Givois
,
A. J.
,
Deü
,
J.-F.
, and
Thomas
,
O.
,
2021
, “
Dynamics of Piezoelectric Structures With Geometric Nonlinearities: A Non-intrusive Reduced Order Modelling Strategy
,”
Comput. Struct.
,
253
, p.
106575
.
39.
Wang
,
X. Q.
,
Song
,
P.
,
Mignolet
,
M. P.
, and
Chen
,
P. C.
,
2021
, “
Reduced Order Nonlinear Damping Model: Formulation and Application to Post-Flutter Aeroelastic Behavior
,”
AIAA J.
,
59
(
10
), pp.
4144
4154
.
40.
Perez
,
R.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2011
, “
Nonlinear Reduced Order Models for Thermoelastodynamic Response of Isotropic and FGM Panels
,”
AIAA J.
,
49
(
3
), pp.
630
641
.
41.
Matney
,
A.
,
Mignolet
,
M. P.
,
Culler
,
A. J.
,
McNamara
,
J. J.
, and
Spottswood
,
S. M.
,
2015
, “
Panel Response Prediction Through Reduced Order Models With Application to Hypersonic Aircraft
,”
Proceedings of the AIAA Science and Technology Forum and Exposition (SciTech2015)
,
Orlando, FL
,
Jan. 5–9
, AIAA Paper No. AIAA-2015-1630.
42.
Matney
,
A. K.
,
Mignolet
,
M. P.
,
Spottswood
,
S. M.
,
Culler
,
A. J.
, and
McNamara
,
J. J.
,
2014
, “
Thermal Reduced Order Model Adaptation to Aero-Thermo-Structural Interactions
,”
Proceedings of the AIAA Science and Technology Forum and Exposition (SciTech2014)
,
National Harbor, MD
,
Jan. 13–17
, Paper No. AIAA-2014-0493.
43.
Matney
,
A.
,
Perez
,
R. A.
,
Spottswood
,
S. M.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2012
, “
Nonlinear Structural Reduced Order Modeling Methods for Hypersonic Structures
,”
Proceedings of the 53rd Structures, Structural Dynamics and Materials Conference
,
Honolulu, HI
,
Apr. 23–26
, AIAA Paper No. AIAA-2012-1972.
44.
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2020
, “
“Enrichments of Structural Bases for the Reduced Order Modeling of Heated Structures Undergoing Nonlinear Geometric Response
,”
Proceedings of the International Modal Analysis Conference, IMAC XXXVIII
,
Houston, TX
,
Feb. 10–13
.
45.
Matney
,
A.
,
Perez
,
R.
, and
Mignolet
,
M. P.
, “
Nonlinear Unsteady Thermoelastodynamic Response of a Panel Subjected to an Oscillating Flux by Reduced Order Models
,”
Proceedings of the 52nd Structures, Structural Dynamics and Materials Conference
,
Denver, CO
,
Apr. 4–7 2011
, Paper No. AIAA-2011–2016.
46.
Perez
,
R.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2010
, “
Steady and Unsteady Nonlinear Thermoelastodynamic Response of Panels by Reduced Order Models
,”
Proceedings of the 51st Structures, Structural Dynamics, and Materials Conference
,
Orlando, FL
,
Apr. 12–15
, Paper No. AIAA-2010-2724.
47.
Kuether
,
R. J.
,
Deaner
,
B. J.
,
Hollkamp
,
J. J.
, and
Allen
,
M. S.
,
2015
, “
Evaluation of Geometrically Nonlinear Reduced-Order Models With Nonlinear Normal Modes
,”
AIAA J.
,
53
(
11
), pp.
3273
3285
.
48.
Gordon
,
R. W.
, and
Hollkamp
,
J. J.
,
2011
, “
Reduced-Order Models for Acoustic Response Prediction of a Curved Panel
,”
Proceedings of the 52nd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference
,
Denver, CO
,
Apr. 4–7
, AIAA Paper No. AIAA-2011-2081.
49.
Perez
,
R.
,
Bartram
,
G.
,
Beberniss
,
T.
,
Wiebe
,
R.
, and
Spottswood
,
S. M.
,
2017
, “
Calibration of Aero-Structural Reduced Order Models Using Full-Field Experimental Measurements
,”
Mech. Syst. Signal Process
,
86
(Part B), pp.
49
65
.
50.
Wiebe
,
R.
,
Perez
,
R. A.
, and
Spottswood
,
S. M.
,
2016
, “
Robust Simulation of Buckled Structures Using Reduced Order Modeling
,”
J. Phys.: Conf. Ser.
,
744
, p.
012118
.
51.
Van Damme
,
C. I.
,
Allen
,
M. S.
, and
Hollkamp
,
J. J.
,
2020
, “
Updating Geometrically Nonlinear Reduced-Order Models Using Nonlinear Modes and Harmonic Balance
,”
AIAA J.
,
58
(
8
), pp.
3553
3568
.
52.
Van Damme
,
C. I.
,
Allen
,
M. S.
, and
Hollkamp
,
J. J.
,
2020
, “
Evaluating Reduced Order Models of Curved Beams for Random Response Prediction Using Static Equilibrium Paths
,”
J. Sound Vib.
,
468
, p.
115018
.
53.
Phlipot
,
G.
,
Wang
,
X. Q.
,
Mignolet
,
M. P.
,
Demasi
,
L.
, and
Cavallaro
,
R.
,
2014
, “
Nonintrusive Reduced Order Modeling for the Nonlinear Geometric Response of Some Joined Wings
,”
Proceedings of the AIAA Science and Technology Forum and Exposition (SciTech2014)
,
National Harbor, MD
,
Jan. 13–17
, AIAA Paper No. AIAA-2014-0151.
54.
Lulf
,
F.
,
Tran
,
D. M.
,
Matthies
,
H. G.
, and
Ohayon
,
R.
,
2015
, “
An Integrated Method for the Transient Solution of Reduced Order Models of Geometrically Nonlinear Structures
,”
Comput. Mech.
,
55
(
2
), pp.
327
344
.
55.
Lin
,
J.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2018
, “
Nonlinear Reduced Order Modeling of Strongly Nonlinear Behavior: A Revisit of a Curved Beam Example
,”
Proceedings of the International Modal Analysis Conference, IMAC XXXVI
,
Orlando, FL
,
Feb. 12–15
.
56.
Lin
,
J.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2019
, “
Nonlinear Reduced Order Modeling of a Cylindrical Shell Exhibiting Mode Veering and Symmetry Breaking
,”
Proceedings of the International Modal Analysis Conference, IMAC XXXVII
,
Orlando, FL
,
Jan. 28–31
.
57.
Wainwright
,
B. A.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2020
, “
Investigation of Out-of-Band Response in Reduced Order Models of Nonlinear Geometric Response
,”
Proceedings of the International Modal Analysis Conference, IMAC XXXVIII
,
Houston, TX
,
Feb. 10–13
.
58.
Wainwright
,
B. A.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2019
, “
Nonlinear Reduced Order Modeling for the Dynamic Response of a Built-up Structure With Strong Asymmetry Through Thickness
,”
Proceedings of the International Modal Analysis Conference, IMAC XXXVII
,
Orlando, FL
,
Jan. 28–31
.
59.
Spottswood
,
S. M.
,
Eason
,
T. G.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2009
, “
Nonlinear Reduced Order Modeling of Curved Beams: A Comparison of Methods
,”
Proceedings of the 50th Structures, Structural Dynamics, and Materials Conference
,
Palm Springs, CA
,
May 4–7
, AIAA Paper No. AIAA-2009-2433.
60.
Spottswood
,
S. M.
,
Hollkamp
,
J. J.
, and
Eason
,
T. G.
,
2010
, “
Reduced-Order Models for a Shallow Curved Beam Under Combined Loading
,”
AIAA J.
,
48
(
1
), pp.
47
55
.
61.
Przekop
,
A.
, and
Rizzi
,
S. A.
,
2007
, “
Dynamic Snap-Through of Thin-Walled Structures by a Reduced-Order Method
,”
AIAA J.
,
45
(
10
), pp.
2510
2519
.
62.
Wang
,
X. Q.
,
Perez
,
R. A.
,
Wainwright
,
B.
,
Wang
,
Y.
, and
Mignolet
,
M. P.
,
2022
, “
Nonintrusive Nonlinear Reduced Order Models for Structures in Large Deformations: Validations to Atypical Structures and Basis Construction Aspects
,”
Vibration
,
5
(
1
), pp.
20
58
.
63.
Capiez-Lernout
,
E.
,
Soize
,
C.
, and
Mignolet
,
M. P.
,
2014
, “
Post-Buckling Nonlinear Static and Dynamical Analyses of Uncertain Cylindrical Shells and Experimental Validation
,”
Comput. Methods Appl. Mech. Eng.
,
271
, pp.
210
230
.
64.
Mignolet
,
M. P.
, and
Soize
,
C.
,
2008
, “
Stochastic Reduced Order Models for Uncertain Geometrically Nonlinear Dynamical Systems
,”
Comput. Methods Appl. Mech. Eng.
,
197
(
45–48
), pp.
3951
3963
.
65.
Wang
,
X. Q.
,
Perez
,
R.
,
Mignolet
,
M. P.
,
Capillon
,
R.
, and
Soize
,
C.
,
2013
, “
Nonlinear Reduced Order Modeling of Complex Wing Models
,”
Proceedings of the 54th Structures, Structural Dynamics and Materials Conference
,
Boston, MA
,
Apr. 8–11
, AIAA Paper No. AIAA-2013-1520.
66.
Lin
,
J.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2020
, “
Non-intrusive Identification of Nonlinear Reduced Order Models: Symmetry and Tuning
,”
Proceedings of the International Modal Analysis Conference, IMAC XXXVIII
,
Houston, TX
,
Feb. 10–13
.
67.
Lin
,
J.
,
2020
, “
Nonlinear Reduced Order Modeling of Structures Exhibiting a Strong Nonlinearity
,”
Ph.D. thesis
,
Arizona State University
, Tempe, AZ.